Мы делаем пиковое обнаружение (и обнаружение долины), используя текущие медианы и медианное абсолютное отклонение Вы можете указать, какое отклонение от бегущей медианы означает пик.
На следующем шаге мы используем биномиальную модель, чтобы проверить, какие регионы содержат больше «экстремальных» значений, чем можно ожидать. Эта модель (в основном тест оценки) приводит к «пиковым областям» вместо отдельных пиков. Поворот вокруг, чтобы получить "области долины", тривиален.
Текущая медиана рассчитывается с использованием функции weightedMedian из пакета aroma.light. Мы используем функцию embed (), чтобы создать список «окон» и применить к нему функцию ядра.
Применение взвешенной медианы:
center <- apply(embed(tmp,wdw),1,weightedMedian,w=weights,na.rm=T)
Здесь tmp - это временный вектор данных и wdw размер окна (должен быть неравномерным). tmp создается путем сложения (wdw-1) / 2 значений NA на каждой стороне вектора данных. веса построены с использованием настраиваемой функции. Для сумасшедших мы используем ту же процедуру, но затем на diff (data) вместо самих данных.
Выполнение Пример кода:
require(aroma.light)
# make.weights : function to make weights on basis of a normal distribution
# n is window size !!!!!!
make.weights <- function(n,
type=c("gaussian","epanechnikov","biweight","triweight","cosinus")){
type <- match.arg(type)
x <- seq(-1,1,length.out=n)
out <-switch(type,
gaussian=(1/sqrt(2*pi)*exp(-0.5*(3*x)^2)),
epanechnikov=0.75*(1-x^2),
biweight=15/16*(1-x^2)^2,
triweight=35/32*(1-x^2)^3,
cosinus=pi/4*cos(x*pi/2),
)
out <- out/sum(out)*n
return(out)
}
# score.test : function to become a p-value based on the score test
# uses normal approximation, but is still quite correct when p0 is
# pretty small.
# This test is one-sided, and tests whether the observed proportion
# is bigger than the hypothesized proportion
score.test <- function(x,p0,w){
n <- length(x)
if(missing(w)) w<-rep(1,n)
w <- w[!is.na(x)]
x <- x[!is.na(x)]
if(sum(w)!=n) w <- w/sum(w)*n
phat <- sum(x*w)/n
z <- (phat-p0)/sqrt(p0*(1-p0)/n)
p <- 1-pnorm(z)
return(p)
}
# embed.na is a modification of embed, adding NA strings
# to the beginning and end of x. window size= 2n+1
embed.na <- function(x,n){
extra <- rep(NA,n)
x <- c(extra,x,extra)
out <- embed(x,2*n+1)
return(out)
}
# running.score : function to calculate the weighted p-value for the chance of being in
# a run of peaks. This chance is based on the weighted proportion of the neighbourhood
# the null hypothesis is calculated by taking the weighted proportion
# of detected peaks in the whole dataset.
# This lessens the need for adjusting parameters and makes the
# method more automatic.
# for a correct calculation, the weights have to sum up to n
running.score <- function(sel,n=20,w,p0){
if(missing(w)) w<- rep(1,2*n+1)
if(missing(p0))p0 <- sum(sel,na.rm=T)/length(sel[!is.na(sel)]) # null hypothesis
out <- apply(embed.na(sel,n),1,score.test,p0=p0,w=w)
return(out)
}
# running.med : function to calculate the running median and mad
# for a dataset. Window size = 2n+1
running.med <- function(x,w,n,cte=1.4826){
wdw <- 2*n+1
if(missing(w)) w <- rep(1,wdw)
center <- apply(embed.na(x,n),1,weightedMedian,w=w,na.rm=T)
mad <- median(abs(x-center))*cte
return(list(med=center,mad=mad))
}
##############################################
#
# Create series
set.seed(100)
n = 1000
series <- diffinv(rnorm(20000),lag=1)
peaks <- apply(embed.na(series,n),1,function(x) x[n+1] < quantile(x,probs=0.05,na.rm=T))
pweight <- make.weights(0.2*n+1)
p.val <- running.score(peaks,n=n/10,w=pweight)
plot(series,type="l")
points((1:length(series))[p.val<0.05],series[p.val<0.05],col="red")
points((1:length(series))[peaks],series[peaks],col="blue")
Пример кода выше разработан для поиска областей с большими колебаниями, а не долинами. Я немного адаптировал, но это не оптимально. Вдобавок ко всему, для серий более 20000 значений вам требуется много памяти, я больше не могу запускать ее на своем компьютере.
В качестве альтернативы, вы можете работать с приближением числовой производной и второй производной для определения долин. В вашем случае это может даже работать лучше. Прагматичный способ вычисления производных и минимумов / максимумов первой производной:
#first derivative
f.deriv <- diff(lowess(series,f=n/length(series),delta=1)$y)
#second derivative
f.sec.deriv <- diff(f.deriv)
#minima and maxima defined by where f.sec.deriv changes sign :
minmax <- cumsum(rle(sign(f.sec.deriv))$lengths)
op <- par(mfrow=c(2,1))
plot(series,type="l")
plot(f.deriv,type="l")
points((1:length(f.deriv))[minmax],f.deriv[minmax],col="red")
par(op)