# test runmin, runmax and runmed
k=15; n=200;
x = rnorm(n,sd=30) + abs(seq(n)-n/4)
col = c("black", "red", "green", "blue", "magenta", "cyan")
plot(x, col=col[1], main = "Moving Window Analysis Functions")
lines(runmin(x,k), col=col[2])
lines(runmed(x,k), col=col[3])
lines(runmax(x,k), col=col[4])
legend(0,.9*n, c("data", "runmin", "runmed", "runmax"), col=col, lty=1 )
#test runmean and runquantile
y=runquantile(x, k, probs=c(0, 0.5, 1, 0.25, 0.75), endrule="constant")
plot(x, col=col[1], main = "Moving Window Quantile")
lines(runmean(y[,1],k), col=col[2])
lines(y[,2], col=col[3])
lines(runmean(y[,3],k), col=col[4])
lines(y[,4], col=col[5])
lines(y[,5], col=col[6])
lab = c("data", "runmean(runquantile(0))", "runquantile(0.5)",
"runmean(runquantile(1))", "runquantile(.25)", "runquantile(.75)")
legend(0,0.9*n, lab, col=col, lty=1 )
#test runmean and runquantile
k =25
m=runmed(x, k)
y=runmad(x, k, center=m)
plot(x, col=col[1], main = "Moving Window Analysis Functions")
lines(m , col=col[2])
lines(m-y/2, col=col[3])
lines(m+y/2, col=col[4])
lab = c("data", "runmed", "runmed-runmad/2", "runmed+runmad/2")
legend(0,1.8*n, lab, col=col, lty=1 )
# numeric comparison between different algorithms
numeric.test = function (n, k) {
eps = .Machine$double.eps ^ 0.5
x = rnorm(n,sd=30) + abs(seq(n)-n/4)
# numeric comparison : runmean
a = runmean(x,k)
b = runmean(x,k, alg="R")
d = runmean(x,k, alg="exact")
e = filter(x, rep(1/k,k))
stopifnot(all(abs(a-b)<eps, na.rm=TRUE));
stopifnot(all(abs(a-d)<eps, na.rm=TRUE));
stopifnot(all(abs(a-e)<eps, na.rm=TRUE));
# numeric comparison : runmin
a = runmin(x,k, endrule="trim")
b = runmin(x,k, endrule="trim", alg="R")
c = apply(embed(x,k), 1, min)
stopifnot(all(a==b, na.rm=TRUE));
stopifnot(all(a==c, na.rm=TRUE));
# numeric comparison : runmax
a = runmax(x,k, endrule="trim")
b = runmax(x,k, endrule="trim", alg="R")
c = apply(embed(x,k), 1, max)
stopifnot(all(a==b, na.rm=TRUE));
stopifnot(all(a==c, na.rm=TRUE));
# numeric comparison : runmad
a = runmad(x,k, endrule="trim")
b = apply(embed(x,k), 1, mad)
stopifnot(all(a==b, na.rm=TRUE));
# numeric comparison : runquantile
a = runquantile(x,k, c(0.3, 0.7), endrule="trim")
b = t(apply(embed(x,k), 1, quantile, probs = c(0.3, 0.7)))
stopifnot(all(abs(a-b)<eps));
}
numeric.test(50, 3) # test different window size vs. vector ...
numeric.test(50,15) # ... length combinations
numeric.test(50,49)
numeric.test(49,49)
# speed comparison
x=runif(100000); k=991;
system.time(runmean(x,k))
system.time(runmean(x,k, alg="R"))
system.time(runmean(x,k, alg="exact"))
system.time(filter(x, rep(1/k,k), sides=2)) #the fastest alternative I know
k=91;
system.time(runmad(x,k))
system.time(apply(embed(x,k), 1, mad)) #the fastest alternative I know
# numerical comparison of round-off error handling
test.runmean = function (x, k) {
a = k*runmean(x,k, alg="exact")
b = k*runmean(x,k, alg="C")
d = k*runmean(x,k, alg="R")
e = k*filter(x, rep(1/k,k))
f = k* c(NA, NA, apply(embed(x,k), 1, mean), NA, NA)
x = cbind(x, a, b, d, e, f)
colnames(x) = c("x","runmean(alg=exact)","runmean(alg=C)",
"runmean(alg=R)","filter","apply")
return(x)
}
a = rep( c(1, 10, -10, -1, 0, 0, 0), 3) # nice-behaving array
b = rep( c(1, 10^20, -10^20, -1, 0, 0, 0), 3) # round-off error prone array
d = rep( c(1, 10^20, 10^40, -10^40, -10^20, -1, 0), 3)
test.runmean(a, 5) #all runmean algorithms give the same result
test.runmean(b, 5) #runmean(alg=R) gives wrong result
test.runmean(d, 5) #only runmean(alg=exact) gives correct result
Run the code above in your browser using DataLab