在嘈杂的时间序列中检测周期最大值（峰值）（使用R语言？）

Question

这个问题是关于一种算法来确定数字序列中最大值的数量和位置。因此，这个问题有统计学味道，但更偏向编程，因为我不关心具体的统计特性，解决方案需要用R语言。使用统计学来回答这个问题是可以的，但不是必须的。 我想要提取时间序列数据中周期极值（即有序数列）。这样的数据的一个例子是太阳耀斑时间序列（约11年周期，在9到14年之间变化）。这些周期不是完美间隔重复的，峰值也不总是相同高度的。 我发现一篇最近的论文描述了一种算法，并且这篇论文实际上以太阳耀斑作为一个例子（图5，Scholkmann等人2012，Algorithms）。我希望这个算法或同样有效的算法可以用作R包。 Scholkmann的“自动多尺度峰值检测”的链接： http://www.mdpi.com/1999-4893/5/4/588 我已经尝试了“pastecs”包中的“turningpoints”函数，但它似乎太敏感了（即检测到太多的峰值）。我想尝试先平滑时间序列，但我不确定这是否是最好的方法（我不是专家）。 感谢您的任何指导。

- rbatt

---

由于这个问题与 R 语言不太相关，更多地涉及统计方法，我认为将其提问在 CrossValidated 而不是 Stack Overflow 更为合适。 - undefined

统计数据可能对回答这个问题有用，但并不一定是必需的。我同意这个问题处于界面上。然而，对这个问题的回答必须包括编程，并且可以（或者不可以）包括统计学。因此，我认为在SO上更合适。 - undefined

Answer 1

这里提供一种使用R语言中的wmtsa包解决问题的方法。我添加了一个自己编写的小函数，以便在wmtsa::wavCWTPeaks将其接近最大值后更容易搜索。

PeakCycle <- function(Data=as.vector(sunspots), SearchFrac=0.02){
    # using package "wmtsa"
    #the SearchFrac parameter just controls how much to look to either side 
    #of wavCWTPeaks()'s estimated maxima for a bigger value
    #see dRange
    Wave <- wavCWT(Data)
    WaveTree <- wavCWTTree(Wave)
    WavePeaks <- wavCWTPeaks(WaveTree, snr.min=5)
    WavePeaks_Times <- attr(WavePeaks, which="peaks")[,"iendtime"]

    NewPeakTimes <- c()
    dRange <- round(SearchFrac*length(Data))
    for(i in 1:length(WavePeaks_Times)){
        NewRange <- max(c(WavePeaks_Times[i]-dRange, 1)):min(c(WavePeaks_Times[i]+dRange, length(Data)))
        NewPeakTimes[i] <- which.max(Data[NewRange])+NewRange[1]-1
    }

    return(matrix(c(NewPeakTimes, Data[NewPeakTimes]), ncol=2, dimnames=list(NULL, c("PeakIndices", "Peaks"))))
}

dev.new(width=6, height=4)
par(mar=c(4,4,0.5,0.5))
plot(seq_along(as.vector(sunspots)), as.vector(sunspots), type="l")
Sunspot_Ext <- PeakCycle()
points(Sunspot_Ext, col="blue", pch=20)

- rbatt

Answer 2

如果峰值几乎是周期性的（具有缓慢波动的周期），例如太阳黑子的例子， 您可以使用 Hilbert变换或经验模态分解来平滑时间序列。

library(EMD)
x <- as.vector(sunspots)
r <- emd(x)
# Keep 5 components -- you may need more, or less.
y <- apply( r$imf[,5:10], 1, sum ) + mean(r$residue)
plot(x, type="l", col="grey")
lines( y, type="l", lwd=2)
n <- length(y)
i <- y[2:(n-1)] > y[1:(n-2)] & y[2:(n-1)] > y[3:n]
points( which(i), y[i], pch=15 )

- Vincent Zoonekynd

---

谢谢你的建议。最终我找到了“wmsta”包，非常有帮助。 - undefined