0

最初に私の問題を説明したいと思います: 私がやりたいのは、30 分ごとのデータを持っている間に、24 時間枠での価格の急上昇の数を計算することです。

たとえば、次のような Stackoverflow の投稿をすべて見ました: Rollapply for time series

（もっと関連性のあるものがある場合は、私に知らせてください;））

私は自分のデータをアップロードできないし、おそらくアップロードすべきでもないので、最小限の例を次に示します。ランダム変数をシミュレートし、それを xts オブジェクトに変換し、ユーザー定義関数を使用して「スパイク」を検出します (もちろん、この場合はかなりばかげています。しかし、エラーを示しています)。

library(xts)
##########Simulate y as a random variable
y <- rnorm(n=100)
##########Add a date variable so i can convert it to a xts object later on
yDate <- as.Date(1:100)
##########bind both variables together and convert to a xts object
z <- cbind(yDate,y)
z <- xts(x=z, order.by=yDate)
##########use the rollapply function on the xts object:
x <- rollapply(z, width=10, FUN=mean)

関数は想定どおりに機能します。前の 10 個の値を取り、平均を計算します。

次に、ピークを見つけるための独自の関数を定義しました。ピークは極大値 (周囲の m ポイントよりも高い) であり、少なくとも時系列の平均 + h と同じ大きさです。これはにつながります：

find_peaks <- function (x, m,h){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])&x[i+1]>mean(x)+h) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}

そして正常に動作します: 例に戻ります:

plot(yDate,y)
#Is supposed to find the points which are higher than 3 points around them
#and higher than the average:
#Does so, so works.
points(yDate[find_peaks(y,3,0)],y[find_peaks(y,3,0)],col="red")

ただし、rollapply()関数を使用すると、次のことが起こります。

x <- rollapply(z,width = 10,FUN=function(x) find_peaks(x,3,0))
#Error in `[.xts`(x, c(z:i, (i + 2):w)) : subscript out of bounds

mパラメータのために、最初のポイントの負のインデックスが int で実行される可能性があるため、エラーが発生する可能性があると最初に考えました。残念ながら、mゼロに設定してもエラーは変わりません。

私もこのエラーを追跡しようとしましたが、ソースが見つかりません。誰か助けてくれませんか？

編集: スパイクの写真:オーストラリアの電力市場のスパイク。find_peaks(20,50) は赤い点がスパイクであると判断し、find_peaks(0,50) はさらに青い点がスパイクであると判断します (したがって、2 番目のパラメーター h が重要です。青い点は明らかに分析したいものではないためです)スパイクについて話すとき）。

4

1 に答える 1

0

あなたが何を求めているのか、私にはまだ完全にはわかりません。データのウィンドウが与えられ、その中心がウィンドウの残りの部分よりも大きいかどうか、同時にウィンドウの平均よりも大きいかどうかを識別したい場合+ h、次のように実行できます。

peakfinder = function(x,h = 0){
  xdat = as.numeric(x)
  meandat = mean(xdat)
  center = xdat[ceiling(length(xdat)/2)]
  ifelse(all(center >= xdat) & center >= (meandat + h),center,NA)
}

y <- rnorm(n=100)
z = xts(y, order.by = as.Date(1:100))
plot(z)
points(rollapply(z,width = 7, FUN = peakfinder, align = "center"), col = "red", pch = 19)

中心点が隣接点よりも大きい場合、必然的に局所平均よりも大きいように見えるかもしれませんが、関数のこの部分はh >= 0. 時系列のグローバル平均を使用する場合は、 の計算をmeandat、 に引数として渡された事前に計算されたグローバル平均に置き換えるだけpeakfinderです。

于 2016-02-04T13:53:56.670 に答える