Wie kann eine robuste Schrittfunktion an eine Zeitreihe angepasst werden?

7

Ich habe eine etwas laute Zeitreihe, die auf verschiedenen Ebenen schwebt.

Zum Beispiel die folgenden Daten:

Geben Sie hier die Bildbeschreibung ein

Ich habe die durchgezogenen Liniendaten zur Verfügung und möchte eine Schätzung für die gestrichelte Linie erhalten. Es sollte stückweise konstant sein.

Welche Algorithmen sollten Sie hier ausprobieren?

Meine bisherigen Ideen drehen sich um 0-Grad-P-Splines (aber wie finde ich heraus, wo die Knoten platziert werden sollen?) Oder Strukturbruchmodelle. Ein Regressionsbaum ist die beste Idee, die ich derzeit habe, aber im Idealfall würde ich nach einer Methode suchen, die die Tatsache berücksichtigt, dass die beiden Ebenen bei y = 250 gleiche y-Werte haben. Wenn ich das richtig verstehe, würde ein Regressionsbaum diese beiden Intervalle in zwei verschiedene Gruppen mit jeweils unterschiedlichen Mittelwerten aufteilen.

Der R-Code, der es generiert hat, ist folgender:

set.seed(20181118)
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
x = 1:400
y = true_fct(x) + rt(length(x), df=1)
plot(x, y, type="l")
lines(x, true_fct(x), lty=2, lwd=3)
Alexander Engelhardt
quelle
2
Wenn Ihre Daten wirklich wie die simulierten aussehen, können Sie kaum etwas Besseres tun, als einen Fenstermedian mit einem sehr kleinen Fenster zu berechnen: Das würde alle Sprünge zuverlässig erkennen. Schätzen Sie die Werte anhand der Mediane der Antworten in jedem dieser erkannten Intervalle. Können Sie daher angeben, ob die impliziten Annahmen der Simulation - große Sprünge, stückweise konstante Mediane und Student t-Fehler - genau die Annahmen sind, die wir treffen sollten?
whuber
1
Vielen Dank für Ihren Kommentar! Ich habe zwei Bemerkungen: (1) Wie würde ich die Intervalle aus dem Fenstermedian erhalten? (2) Die Annahmen sind stückweise konstante Mediane und merkliche Sprünge, aber ich weiß nichts über die Fehlerverteilung, außer der Tatsache, dass große Ausreißer auftreten können.
Alexander Engelhardt
Manchmal funktionieren einfache nichtparametrische Methoden, wenn das Problem einfach ist. Ich möchte, dass Sie einen anspruchsvolleren / realistischeren Datensatz simulieren, in dem eine eingebettete Arima-Struktur und möglicherweise ein oder zwei saisonale Impulse vorhanden sind. Umfassende Ansätze für solche Probleme müssen autoregressive Strukturen und Anomalien während der Verarbeitung berücksichtigen und isolieren. Sie können eine andere Frage stellen und den etwas realistischeren Datensatz hinzufügen.
IrishStat
Ich sollte auch hinzufügen, wenn die Pegel- / Schrittverschiebungen gegenüber dem Fehlerprozess so groß sind, dass nichtparametrische Methoden eine nützliche Rolle spielen können und weniger, wenn dieses Verhältnis kleiner wird
IrishStat

Antworten:

7

Eine einfache, robuste Methode, um mit solchen Störungen umzugehen, ist die Berechnung von Medianen.

Ein rollierender Median über ein kurzes Fenster erkennt alle bis auf die kleinsten Sprünge, während Mediane der Antwort innerhalb von Intervallen zwischen erkannten Sprüngen ihre Pegel zuverlässig schätzen. (Sie können diese letztere Schätzung durch eine robuste Schätzung ersetzen, die von den Ausreißern nicht beeinflusst wird.)

Sie sollten diesen Ansatz mit realen oder simulierten Daten abstimmen, um akzeptable Fehlerraten zu erzielen. Zum Beispiel fand ich es für die Simulation in der Frage gut, das zweite und das 98. Perzentil zu verwenden, um Schwellenwerte für die Erkennung der Sprünge festzulegen. Unter anderen Umständen - beispielsweise wenn viele Sprünge auftreten könnten - würden zentralere Perzentile besser funktionieren.

Hier ist das Ergebnis, das (a) die drei Sprünge als rote Punkte und (b) die vier geschätzten Ebenen als hellblaue Linien zeigt.

Zahl

Es wird geschätzt, dass die Sprünge bei den Indizes 100, 200, 250 auftreten (genau dort, wo sie durch die Simulation auftreten), und die resultierenden Werte werden auf 199,6, 249,8, 300,0 und 250,2 geschätzt: alle innerhalb von 0,4 der tatsächlichen zugrunde liegenden Werte.

Dieses hervorragende Verhalten bleibt bei wiederholten Simulationen bestehen (Entfernen des set.seedBefehls am Anfang).

Hier ist der RCode.

#
# Rolling medians.
#
rollmed <- function(x, k=3) {
  n <- length(x)
  x.med <- sapply(1:(n-k+10), function(i) median(x[i + 0:(k-1)]))
  l <- floor(k/2)
  c(rep(NA, l), x.med, rep(NA, k-l))
}
y.med <- rollmed(y, k=5)
#
# Changepoint analysis.
#
dy <- diff(y.med)
fourths <- quantile(dy, c(1,49)/50, na.rm=TRUE)
thresholds <- fourths + diff(fourths)*2.5*c(-1,1)
jumps <- which(dy < thresholds[1] | dy > thresholds[2]) + 1

points(jumps, y.med[jumps], pch=21, bg="Red")
#
# Plotting.
#
limits <- c(1, jumps, length(y)+1)
y.hat <- rep(NA, length(jumps)+1)
for (i in 1:(length(jumps)+1)) {
  j0 <- limits[i]
  j1 <- limits[i+1]-1
  y.hat[i] <- median(y[j0:j1])
  lines(x[j0:j1], rep(y.hat[i], j1-j0+1), col="skyblue", lwd=2)
}
whuber
quelle
+1, aber der Teil des Codes "Änderungspunktanalyse" ist für einige Benutzer möglicherweise nicht ganz klar. Vielleicht können Sie also kommentieren, was dort passiert?
Tim
@ Tim Danke für den Vorschlag. Der Zweck des ersten Absatzes besteht darin, diesen Algorithmus zu erläutern. Ich möchte die Details der Implementierung herunterspielen, da sie unwichtig sind: Es sollte ausreichen, eine robuste Ausreißer-Screening-Methode auf die Residuen anzuwenden.
whuber
Möglicherweise möchten Sie zoo::rollmedianeine ähnliche Funktion in Betracht ziehen , um Ihren Code zu vereinfachen.
usεr11852
@ usεr11852 Vielen Dank. Ich bin mir dessen bewusst, habe michzoo aber entschieden, es nicht zu benutzen, weil ich faul bin! Es war schneller und einfacher zu schreiben, rollmedals die Argumentaufrufe für jede Funktion zu überprüfen, die möglicherweise bereits verfügbar ist. Außerdem gefällt mir, wie rollmeddeutlich zeigt, was ich tue, anstatt die Details hinter einer Black Box zu verstecken.
whuber
Kein Problem. :) (Ich war mir sicher, dass Sie davon wussten zoo, ich war mir nicht sicher, ob Sie es nicht
freiwillig
3

Wenn Sie immer noch an einer Glättung mit L0-Strafen interessiert sind, würde ich einen Blick auf die folgende Referenz werfen: "Visualisierung genomischer Veränderungen durch segmentierte Glättung mit einer L0-Strafe" - DOI: 10.1371 / journal.pone.0038230 (eine nette Einführung in die Whittaker Smoother finden Sie in P. Eilers Papier "A Perfect Smoother" (DOI: 10.1021 / ac034173t). Um Ihr Ziel zu erreichen, müssen Sie natürlich ein wenig an der Methode arbeiten.

Grundsätzlich benötigen Sie 3 Zutaten:

  1. Je glatter - ich würde den Whittaker-Glätter verwenden. Außerdem werde ich die Matrixvergrößerung verwenden (siehe Eilers und Marx, 1996 - "Flexibles Glätten mit B-Splines und Strafen", S.101).
  2. Quantile Regression - Ich werde das R-Paket Quantreg (Rho = 0,5) für Faulheit verwenden :-)
  3. L0-Strafe - Ich werde der erwähnten "Visualisierung genomischer Veränderungen durch segmentierte Glättung unter Verwendung einer L0-Strafe" folgen - DOI: 10.1371 / journal.pone.0038230

Natürlich müssten Sie auch die optimale Glättungsmenge auswählen. Dies wird von meinen Zimmermannsaugen für dieses Beispiel gemacht. Sie können die Kriterien in DOI verwenden: 10.1371 / journal.pone.0038230 (S. 5, aber ich habe es in Ihrem Beispiel nicht versucht).

Unten finden Sie einen kleinen Code. Ich habe einige Kommentare als Leitfaden hinterlassen.

# Cross Validated example
rm(list = ls()); graphics.off(); cat("\014")

library(splines)
library(Matrix)
library(quantreg)

# The data
set.seed(20181118)
n = 400
x = 1:n
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
y = true_fct(x) + rt(length(x), df = 1)

# Prepare bases - Identity matrix (Whittaker)
# Can be changed for B-splines
B = diag(1, n, n)

# Prepare penalty - lambda parameter fix
nb = ncol(B)
D = diff(diag(1, nb, nb), diff = 1)
lambda = 1e2

# Solve standard Whittaker - for initial values
a = solve(t(B) %*% B + crossprod(D), t(B) %*% y, tol = 1e-50)    

# est. loop with L0-Diff penalty as in DOI: 10.1371/journal.pone.0038230
p = 1e-6
nit = 100
beta = 1e-5

for (it in 1:nit) {
  ao = a

  # Penalty weights
  w = (c(D %*% a) ^ 2  + beta ^ 2) ^ ((p - 2)/2)
  W = diag(c(w))

  # Matrix augmentation
  cD = lambda * sqrt(W) %*% D
  Bp = rbind(B, cD)
  yp =  c(y, 1:nrow(cD)*0)

  # Update coefficients - rq.fit from quantreg
  a = rq.fit(Bp, yp, tau = 0.5)$coef

  # Check convergence and update
  da = max(abs((a - ao)/ao))
  cat(it, da, '\n')
  if (da < 1e-6) break
}

# Fit 
v = B %*% a

# Show results
plot(x, y, pch = 16, cex = 0.5)
lines(x, y, col = 8, lwd = 0.5)
lines(x, v, col = 'blue', lwd = 2)
lines(x, true_fct(x), col = 'red', lty = 2, lwd = 2)
legend("topright", legend = c("True Signal", "Smoothed signal"), 
       col = c("red", "blue"), lty = c(2, 1))

Geben Sie hier die Bildbeschreibung ein PS. Dies ist meine erste Antwort auf Cross Validated. Ich hoffe es ist nützlich und klar genug :-)

Gi_F.
quelle
1

Ich würde in Betracht ziehen, Ruey Tsays Artikel Ausreißer, Pegelverschiebungen und Varianzänderungen in Zeitreihen-Differenzierungsmodellen mit AR1- und 21-Ausreißern zu verwenden.

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Wir haben die Differenzierung deaktiviert und die Pegelverschiebungen werden speziell angezeigt.

Geben Sie hier die Bildbeschreibung ein

Tom Reilly
quelle
1
Ich frage mich, ob Sie die Betonung auf "robust" im Titel der Frage übersehen haben, da jede Methode, die zusätzlich zu den 3 tatsächlichen Sprüngen 18 Störparameter (entsprechend den in der Simulation eingeführten Ausreißern) identifiziert, kaum als robust (oder sparsam) angesehen werden kann diese Angelegenheit).
whuber
Das ist eine robuste Lösung. Ich bin mir nicht sicher, warum Sie gegen die Identifizierung und Anpassung von Ausreißern sind, aber es gibt eine Welt der Forschung, die dies unterstützt, und natürlich unsere Erfahrungen. Diese anderen Variablen sind Ausreißer. Ich habe ein Diagramm hinzugefügt, das die historischen Daten und eine bereinigte Version zeigt, um den Unterschied gegenüberzustellen.
Tom Reilly
1
Könnten Sie explizit angeben, wie hoch Ihre Schätzung der Schrittfunktion ist?
whuber
1
In der Periode 100 (x3), 200 (x2), 250 (x4) gibt es ein Flag, das den Schritt anzeigt. Der Differenzierungsoperator macht es etwas schwieriger zu sehen, aber der Effekt ist der gleiche. Ich habe ein Modell ohne Differenzierung hinzugefügt.
Tom Reilly