Wie kann ich verrauschte Patches in einer Zeitreihe hervorheben?

9

Ich habe viele Zeitreihendaten - Wasserstände und Geschwindigkeiten gegen die Zeit. Es ist die Ausgabe einer hydraulischen Modellsimulation. Als Teil des Überprüfungsprozesses, um zu bestätigen, dass das Modell die erwartete Leistung erbringt, muss ich jede Zeitreihe zeichnen, um sicherzustellen, dass die Daten keine "Wackelbewegungen" enthalten (siehe Beispiel für geringfügiges Wackeln unten). Die Verwendung der Benutzeroberfläche der Modellierungssoftware ist eine ziemlich langsame und mühsame Methode, um diese Daten zu überprüfen. Ich habe daher ein kurzes VBA-Makro geschrieben, um verschiedene Datenbits aus dem Modell einschließlich der Ergebnisse in Excel zu importieren und alle gleichzeitig zu zeichnen. Ich hoffe, ein weiteres kurzes VBA-Makro schreiben zu können, um die Zeitreihendaten zu analysieren und verdächtige Abschnitte hervorzuheben.

Mein einziger Gedanke ist bisher, dass ich eine Analyse der Steigung der Daten durchführen könnte. Überall dort, wo sich die Steigung innerhalb eines bestimmten Suchfensters schnell schnell von positiv zu negativ ändert, kann sie als instabil eingestuft werden. Vermisse ich einfachere Tricks? Im Wesentlichen sollte eine "stabile" Simulation eine sehr glatte Kurve liefern. Jegliche plötzlichen Änderungen sind wahrscheinlich auf eine Instabilität der Berechnungen zurückzuführen.

Beispiel kleine Instabilität

davehughes87
quelle
1
Lesen Sie Tukeys Buch EDA für eine Reihe einfacher Methoden. Zu Beginn des Buches beschreibt er beispielsweise einfache Glätter und ihre Verwendung, um Residuen zu erhalten. Eine nachfolgende Glättung der absoluten Residuen würde die lokale Variabilität Ihrer Kurven darstellen, hoch gehen, wenn Sie schnelle, plötzliche oder äußere Änderungen haben, und ansonsten niedrig bleiben. Viele viel ausgefeiltere Methoden sind möglich, aber vielleicht würde dies ausreichen. Tukeys Smoothers sind in VBA relativ einfach zu codieren: Ich habe es geschafft .
whuber
@whuber Dies ist im Wesentlichen die Leistung des gleitenden Hochpassfilters?
Amöbe
@ Amöbe Vielleicht. Mein Verständnis solcher Filter ist, dass sie nicht vollständig lokal und definitiv nicht robust sind, während Tukeys Glätter diese beiden wichtigen Eigenschaften haben. (Heutzutage verwenden die Leute Löss oder GAMs zum Glätten, was in Ordnung ist, aber diese sind viel weniger einfach zu implementieren.)
whuber

Antworten:

10

1- -αα

Zahl

1201α=0,201

αα0,20α0,20

Die Details des Glatten spielen keine Rolle. In diesem Beispiel wurde eine Lössglättung (implementiert Rwie loessbei span=0.05der Lokalisierung) verwendet, aber selbst ein Mittelwert mit Fenster hätte gut funktioniert. Um die absoluten Residuen zu glätten, habe ich einen Fenstermittelwert der Breite 17 (ca. 24 Minuten) gefolgt von einem Fenstermedian ausgeführt. Diese Fensterglättungen sind relativ einfach in Excel zu implementieren. Eine effiziente VBA-Implementierung (für ältere Versionen von Excel, aber der Quellcode sollte auch in neuen Versionen funktionieren) ist unter http://www.quantdec.com/Excel/smoothing.htm verfügbar .


R Code

#
# Emulate the data in the plot.
#
xy <- matrix(c(0, 96.35,  0.3, 96.6, 0.7, 96.7, 1, 96.73, 1.5, 96.74, 2.5, 96.75, 
               4, 96.9, 5, 97.05, 7, 97.5, 10, 98.5, 12, 99.3, 12.5, 99.35, 
               13, 99.355, 13.5, 99.36, 14.5, 99.365, 15, 99.37, 15.5, 99.375, 
               15.6, 99.4, 15.7, 99.41, 20, 99.5, 25, 99.4, 27, 99.37),
             ncol=2, byrow=TRUE)
n <- 401
set.seed(17)
noise.x <- cumsum(rexp(n, n/max(xy[,1])))
noise.y <- rep(c(-1,1), ceiling(n/2))[1:n]
noise.amp <- runif(n, 0.8, 1.2) * 0.04
noise.amp <- noise.amp * ifelse(noise.x < 16 | noise.x > 24.5, 0.05, 1)
noise.y <- noise.y * noise.amp

g <- approxfun(noise.x, noise.y)
f <- splinefun(xy[,1], xy[,2])
x <- seq(0, max(xy[,1]), length.out=1201)
y <- f(x) + g(x)
#
# Plot the data and a smooth.
#
par(mfrow=c(1,2))
plot(range(xy[,1]), range(xy[,2]), type="n", main="Data", sub="With Smooth",
     xlab="Time (hours)", ylab="Water Level")
abline(h=seq(96, 100, by=0.5), col="#e0e0e0")
abline(v=seq(0, 30, by=5), col="#e0e0e0")
#curve(f(x) + g(x), xlim=range(xy[,1]), col="#2070c0", lwd=2, add=TRUE, n=1201)
lines(x,y, type="l", col="#2070c0", lwd=2)

span <- 0.05
fit <- loess(y ~ x, span=span)
y.hat <- predict(fit)
lines(fit$x, y.hat)
#
# Plot the absolute residuals to the smooth.
#
r <-  abs(resid(fit))
plot(fit$x, r, type="l", col="#808080",
     main="Absolute Residuals", sub="With Smooth and a Threshold",
     xlab="Time hours", ylab="Residual Water Level")
#
# Smooth plot an indicator of the smoothed residuals.
#
library(zoo)
smooth <- function(x, window=17) {
  x.1 <- rollapply(ts(x), window, mean)
  x.2 <- rollapply(x.1, window, median)
  return(as.vector(x.2))
}
alpha <- 0.2
threshold <- quantile(r, 1-alpha)
abline(h=threshold, lwd=2, lty=3)
r.hat <- smooth(r >threshold)
x.hat <- smooth(fit$x)
z <- max(r)/2 * (r.hat > alpha)
lines(x.hat, z, lwd=2, col="#c02020")
par(mfrow=c(1,1))
whuber
quelle
1
+1. Haben Sie die Daten irgendwie aus dem Plot des OP herausgekratzt?
Amöbe
2
@Amoeba Das wäre zu viel Mühe, besonders für die wackeligen Stellen nach 15 Stunden. Ich musterte ein Dutzend Punkte auf der Kurve, zeichnete einen Spline, fügte einige Zwischenpunkte ein, um die seltsamen Spitzen zu beseitigen, die ein Spline erzeugen kann, und fügte einen stark negativ heteroskedastischen korrelierten Fehler hinzu. Der gesamte Vorgang dauerte nur wenige Minuten und führte zu einem Datensatz, der qualitativ dem in der Frage gezeigten entspricht.
whuber
Ich habe mich gefragt, wie Sie die Daten aus meiner Handlung erhalten haben! Prost! Ich werde es versuchen.
Davehughes87
FWIW, ich habe den Code gepostet, mit dem ich die Illustration gemacht habe. Auch wenn es nicht VBA ist, wird es vielleicht die Details klären. (cc @amoeba)
whuber