检测嘈杂时间序列中的周期最大值(峰值)(在 R 中?) [英] Detecting cycle maxima (peaks) in noisy time series (In R?)

问题描述

这个问题是关于一种算法，用于确定数字序列中最大值的数量和位置.因此，这个问题有统计的味道，但更倾向于编程，因为我对具体的统计属性不感兴趣，解决方案需要在R中.使用统计来回答这个问题是可以的，但不是必需的.

我想提取时间序列数据(即有序的数字序列)中的最大周期数.此类数据的一个示例是太阳耀斑时间序列(~11 年周期，在 9 和 14 年之间).循环不会以完美的间隔重复，峰的高度也不总是相同.

我发现一篇最近的论文描述了一个算法，该论文实际上以太阳耀斑为例(图 5，Scholkmann 等人，2012，算法).我希望这个算法，或者同样有效的算法，可以作为 R 包提供.

链接到 Scholkmann 有关基于多尺度的自动峰值检测"的论文

This question is about an algorithm for determining the number and location of maxima in a sequence of numbers. Thus, there is a statistical flavor to the question, but it is more leaning towards programming, because I am not interested in the specific statistical properties, and the solution needs to be in R. The use of statistics to answer this question is OK, but not a requirement.

I want to extract maxima of cycles in time series data (i.e., an ordered sequence of numbers). An example of such data is the solar flare time series (~11 year cycle, between 9 & 14 years). The cycles don't repeat at a perfect interval, and the peaks aren't always the same height.

I found a recent paper describing an algorithm for this, and the paper actually uses solar flares as an example (Figure 5, Scholkmann et al. 2012, Algorithms). I was hoping that this algorithm, or an equally effective algorithm, was available as an R package.

Link to Scholkmann paper on "automatic multiscale-based peak detection" http://www.mdpi.com/1999-4893/5/4/588

I've tried the "turningpoints" function in the "pastecs" package but it seemed to be too sensitive (i.e., detected too many peaks). I thought of trying to smooth the time series first, but I'm not sure if this is the best approach (I'm no expert).

Thanks for any pointers.

解决方案

Here is a solution involving the wmtsa package in R. I added my own little function to facilitate the searching of maxima once the wmtsa::wavCWTPeaks got it close.

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)

这篇关于检测嘈杂时间序列中的周期最大值(峰值)(在 R 中?)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！