在噪声时间序列中检测循环最大值（峰值）（In R？）

问题描述：

这个问题是关于确定数字序列中最大值数量和位置的算法。因此，这个问题有一个统计口味，但它更倾向于编程，因为我对特定的统计特性不感兴趣，并且解决方案需要在R.使用统计来回答这个问题是可以的，但不是要求。在噪声时间序列中检测循环最大值（峰值）（In R？）

我想提取时间序列数据中循环的最大值（即一个有序的数字序列）。这种数据的一个例子是太阳耀斑时间序列（〜11年周期，14年间的）。这些周期不会以完美的间隔重复，并且峰值并不总是相同的高度。

我发现了一篇最近的论文，描述了一个算法，并且该论文实际上使用了太阳耀斑作为例子（图5，Scholkmann et al。2012，算法）。我希望这个算法或者一个同样有效的算法可以作为一个R包提供。

链接到Scholkmann一篇关于“自动的基于多尺度峰值检测” http://www.mdpi.com/1999-4893/5/4/588

我试过在“pastecs”包“turningpoints”的功能，但它似乎过于敏感（即检测到过许多高峰）。我想先尝试平滑时间系列，但我不确定这是否是最好的方法（我不是专家）。

感谢您的指点。

由于这个问题是没有这么多关于R而是关于统计方法，我认为这是比SO更适合交叉验证。 – RoyalTS 2013-05-02 16:36:34

统计数据可能对回答问题很有用，但它们不一定是必需的。不过，我同意这个问题存在于界面上。但是，这个问题的答案*必须包括编程，*可能*（或不可以）包含统计。因此，我认为它更适合于SO。 – rbatt 2015-08-11 18:27:05

答

这里是涉及wmtsa包R.我说我自己的小功能，方便最大值的搜索，一旦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)

答

如果峰几乎周期（以缓慢波动时期），如太阳黑子例如，可以使用Hilbert transform或empirical mode decomposition平稳时间序列。

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)

Sunspots

感谢您的建议。我最终找到了“wmsta”包，这非常有帮助。 – rbatt 2013-05-02 20:35:59

在噪声时间序列中检测循环最大值（峰值）（In R？）

相关推荐