Time-Series

如何突出顯示時間序列中的嘈雜補丁?

  • February 10, 2017

我有很多時間序列數據——水位和速度與時間的關係。它是水力模型模擬的輸出。作為確認模型按預期執行的審查過程的一部分,我必須繪製每個時間序列以確保數據中沒有“擺動”(參見下面的小擺動示例)。使用建模軟件的 UI 是檢查這些數據的一種非常緩慢且費力的方法。因此,我編寫了一個簡短的 VBA 宏來將模型中的各種數據(包括結果)導入 Excel,並一次將它們全部繪製出來。我希望編寫另一個簡短的 VBA 宏來分析時間序列數據並突出顯示任何可疑的部分。

到目前為止,我唯一的想法是我可以對數據的斜率進行一些分析。在給定的搜索窗口內,斜率多次從正數快速變為負數的任何地方都可以歸類為不穩定的。我錯過了任何更簡單的技巧嗎?本質上,“穩定”模擬應該提供非常平滑的曲線。任何突然的變化都可能是計算不穩定的結果。

輕微不穩定性示例

為簡單起見,我建議分析殘差相對於數據的穩健平滑的大小(絕對值)。對於自動檢測,請考慮用指標替換這些大小: 1 當它們超過某個高分位數時,例如在級別, 否則為 0。平滑此指標並突出顯示任何超過的平滑值.

數字

左邊的圖形藍色的數據點以及黑色的穩健的局部平滑。右圖顯示了該平滑殘差的大小。黑色虛線是它們的第 80 個百分位數(對應於)。紅色曲線的構造如上所述,但已按比例縮放(從和) 到用於繪圖的絕對殘差的中間範圍。

變化允許控制精度。在這種情況下,設置少於識別 22 小時左右的噪音中的短暫間隙,同時設置比…更棒0小時附近也出現了快速變化。

光滑的細節並不重要。在這個例子中,使用了一個 loess smooth (在R本地化loessspan=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))

引用自:https://stats.stackexchange.com/questions/261115

comments powered by Disqus