Jak mogę wyróżnić hałaśliwe łatki w szeregu czasowym?

9

Mam wiele danych szeregów czasowych - poziomy wody i prędkości w funkcji czasu. Jest to wynik symulacji modelu hydraulicznego. W ramach procesu przeglądu, aby potwierdzić, że model działa zgodnie z oczekiwaniami, muszę wykreślić każdy szereg czasowy, aby upewnić się, że w danych nie ma żadnych „wahnięć” (patrz przykład mniejszego wahania poniżej). Korzystanie z interfejsu użytkownika oprogramowania do modelowania jest dość powolnym i pracochłonnym sposobem sprawdzania tych danych. Dlatego napisałem krótkie makro VBA, aby zaimportować różne bity danych z modelu, w tym wyniki do Excela i wykreślić je wszystkie naraz. Mam nadzieję napisać kolejne krótkie makro VBA, aby przeanalizować dane szeregów czasowych i wyróżnić wszelkie podejrzane sekcje.

Do tej pory myślałem tylko, że mógłbym przeprowadzić analizę nachylenia danych. Wszędzie tam, gdzie nachylenie gwałtownie zmienia się od wielokrotnego dodatniego do ujemnego w danym oknie wyszukiwania, można je zaklasyfikować jako niestabilne. Czy brakuje mi prostszych sztuczek? Zasadniczo „stabilna” symulacja powinna zapewniać bardzo płynną krzywą. Wszelkie nagłe zmiany mogą wynikać z niestabilności obliczeń.

Przykład niewielkiej niestabilności

davehughes87
źródło
1
Przeczytaj książkę EDA Tukeya, aby zapoznać się z zestawem prostych metod. Na przykład we wczesnej części książki opisuje proste wygładzacze i ich zastosowanie do uzyskiwania pozostałości. Kolejne wygładzenie wartości bezwzględnych wykreśli lokalną zmienność twoich krzywych, osiągając wysoki poziom tam, gdzie masz szybkie, nagłe lub odległe zmiany, i w przeciwnym razie pozostawałby niski. Możliwych jest wiele bardziej wyrafinowanych metod, ale być może to wystarczy. Wygładzacze Tukeya są stosunkowo łatwe do kodowania w VBA: Zrobiłem to .
whuber
@whuber Jest to zasadniczo moc przesuwanego filtra górnoprzepustowego?
ameba
@amoeba Może. Rozumiem, że takie filtry są takie, że nie są one całkowicie lokalne i zdecydowanie nie są niezawodne, podczas gdy wygładzacze Tukeya mają obie te ważne właściwości. (W dzisiejszych czasach ludzie używają Loess lub GAM do wygładzania, co jest w porządku, ale są one znacznie mniej łatwe do wdrożenia.)
whuber

Odpowiedzi:

11

Dla uproszczenia sugerowałbym przeanalizowanie wielkości (wartości bezwzględnych) reszt w stosunku do solidnego wygładzenia danych. W przypadku automatycznego wykrywania rozważ zastąpienie tych rozmiarów wskaźnikiem: 1, gdy przekroczą jakiś wysoki kwantyl, powiedzmy na poziomie , a 0 w przeciwnym razie. Wygładź ten wskaźnik i wyróżnij wygładzone wartości przekraczające .1αα

Postać

Grafika po lewej stronie przedstawia punktów danych w kolorze niebieskim wraz z solidnym, lokalnym wygładzeniem w kolorze czarnym. Grafika po prawej stronie pokazuje rozmiary resztek tej gładkości. Czarna linia przerywana to ich 80 percentyl (odpowiadający ). Czerwona krzywa jest skonstruowana jak opisano powyżej, ale została przeskalowana (od wartości i ) do średnicy bezwzględnych reszt do wykreślenia.1201α=0.201

Zmienianie pozwala kontrolować precyzję. W tym przypadku ustawienie poniżej identyfikuje krótką przerwę w hałasie około 22 godzin, podczas gdy ustawienie większe niż powoduje także szybką zmianę w pobliżu 0 godzin.αα0.20α0.20

Szczegóły gładkości nie mają większego znaczenia. W tym przykładzie less gładka (realizowane w Rjak loessze span=0.05aby go zlokalizować) użyto, ale nawet z okienkiem średnią zrobiłby grzywny. Aby wygładzić resztki bezwzględne, przeprowadziłem okienkową średnią szerokość 17 (około 24 minut), a następnie środkową okienkową. Te wygładzone okna są stosunkowo łatwe do wdrożenia w programie Excel. Wydajna implementacja VBA (dla starszych wersji programu Excel, ale kod źródłowy powinien działać nawet w nowych wersjach) jest dostępna pod adresem http://www.quantdec.com/Excel/smoothing.htm .


R Kod

#
# 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
źródło
1
+1. Czy w jakiś sposób zeskrobałeś dane z fabuły OP?
ameba
2
@Amoeba Byłoby to zbyt wielkim problemem, szczególnie dla krętych kawałków po 15 godzinach. Spojrzałem kilkanaście punktów na krzywą, narysowałem splajn, wstawiłem kilka punktów pośrednich, aby pozbyć się dziwnych kolców, które może wytworzyć splajn, i dodałem silnie negatywnie heteroscedastyczny błąd skorelowany. Cały proces zajął zaledwie kilka minut i zaowocował zestawem danych jakościowo takim jak ten pokazany w pytaniu.
whuber
Zastanawiałem się, skąd masz dane z mojej fabuły! Twoje zdrowie! Spróbuję.
davehughes87
FWIW, opublikowałem kod, którego użyłem do zrobienia ilustracji. Chociaż nie jest to VBA, może wyjaśni szczegóły. (cc @amoeba)
whuber