如何突出顯示時間序列中的嘈雜補丁?
我有很多時間序列數據——水位和速度與時間的關係。它是水力模型模擬的輸出。作為確認模型按預期執行的審查過程的一部分,我必須繪製每個時間序列以確保數據中沒有“擺動”(參見下面的小擺動示例)。使用建模軟件的 UI 是檢查這些數據的一種非常緩慢且費力的方法。因此,我編寫了一個簡短的 VBA 宏來將模型中的各種數據(包括結果)導入 Excel,並一次將它們全部繪製出來。我希望編寫另一個簡短的 VBA 宏來分析時間序列數據並突出顯示任何可疑的部分。
到目前為止,我唯一的想法是我可以對數據的斜率進行一些分析。在給定的搜索窗口內,斜率多次從正數快速變為負數的任何地方都可以歸類為不穩定的。我錯過了任何更簡單的技巧嗎?本質上,“穩定”模擬應該提供非常平滑的曲線。任何突然的變化都可能是計算不穩定的結果。
為簡單起見,我建議分析殘差相對於數據的穩健平滑的大小(絕對值)。對於自動檢測,請考慮用指標替換這些大小: 1 當它們超過某個高分位數時,例如在級別, 否則為 0。平滑此指標並突出顯示任何超過的平滑值.
左邊的圖形藍色的數據點以及黑色的穩健的局部平滑。右圖顯示了該平滑殘差的大小。黑色虛線是它們的第 80 個百分位數(對應於)。紅色曲線的構造如上所述,但已按比例縮放(從和) 到用於繪圖的絕對殘差的中間範圍。
變化允許控制精度。在這種情況下,設置少於識別 22 小時左右的噪音中的短暫間隙,同時設置比…更棒0小時附近也出現了快速變化。
光滑的細節並不重要。在這個例子中,使用了一個 loess smooth (在
R
本地化loess
中span=0.05
實現),但即使是窗口平均值也可以做得很好。為了平滑絕對殘差,我運行了寬度為 17(約 24 分鐘)的窗口平均值,然後是窗口中位數。這些窗口平滑在 Excel 中相對容易實現。http://www.quantdec.com/Excel/smoothing.htm提供了一個高效的 VBA 實現(適用於舊版本的 Excel,但源代碼即使在新版本中也應該可以工作)。
R
代碼# # 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))