Ошибка в rollapply: индекс за пределами

Question

Прежде всего я хотел бы описать свою проблему: Что я хочу сделать, так это рассчитать количество всплесков цен в 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 

я первый подумал, ну, возможно, ошибка возникает из-за того, что она может запускать int отрицательный индекс для первых точек из-за параметра m. К сожалению, установка m на ноль не изменяет ошибку.

Я тоже попытался проследить эту ошибку, но не нашел источник. Может ли кто-нибудь помочь мне здесь?

Edit: Изображение шипов: Spikes on the australian Electricity Market. find_peaks(20,50) determines the red points to be spikes, find_peaks(0,50) additionally finds the blue ones to be spikes (therefore, the second parameter h is important, because the blue points are clearly not what we want to analyse when we talk about spikes).

Answer 1

Я до сих пор не совсем уверен, что это то, что вы после этого. Исходя из предположения, что данное окно данных, которые вы хотите, чтобы определить, является ли больше, чем остальная частью окна, в то же время его центр как больше, чем среднее значение окна + 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.