如果我正在使用的值是一个异常值,那么我需要计算一个ts,以前的30个值.
我正在使用的数据的维度为600列到200000行.所以我想利用数据表速度的好处.
我的功能是:
es_outlier<-function(vect){
qq =quantile(vect, prob=c(0.25,0.75), na.rm=T)
q3=qq[2]
IC=q3-qq[1]
limSup=q3+IC*1.5
vector_final=abs(vect)>limSup
return(vector_final[length(vect)] )
}
示例表将是:
library(data.table)
dt<-data.table(x1=runif(50000), x2=runif(50000))
dt$x1[555]<-2000
dt$x2[556]<-2000
我可以用动物园包解决这个问题:
zoo::rollapply(dt,30,es_outlier, fill=NA,align='right')
但它需要花费很多时间,而且比实际数据要少.
我想要像:
dt[, (nom):=lapply(.SD,function, n=30)]
我尝试使用Rcpp,但它没有分位数功能.
是否有更快的方式来应用我的功能?
PS:对于一个小表,函数返回:
x<-data.frame(x1=1:8, x2=c(1:7,2000))
x_dt<-data.table(x)
zoo::rollapply(x_dt,5,es_outlier, fill=NA,align='right')
x1 x2
NA NA
NA NA
NA NA
NA NA
FALSE FALSE
FALSE FALSE
FALSE FALSE
FALSE TRUE
最佳答案 建议存储排序后的矢量,这样当从一个窗口移动到另一个窗口时,只需要添加1个新元素.虽然不是很好的加速…
set.seed(25L)
N <- 50000
dt <- data.frame(x1=runif(N), x2=runif(N))
dt$x1[555] <- 2000
dt$x2[556] <- 2000
wl <- 30
####################################################################################################
#' Calculate IQR for a sorted vector with 30 observations
#'
#' @details assume that sorted is sorted. using type 7 in ?quantile.
#'
#' @param sorted sorted numeric vector
#'
#' @return the interquartile range
#'
iqr30obs <- function(sorted) {
c(sorted[8] + 0.25 * (sorted[9] - sorted[8]), sorted[22] + 0.75 * (sorted[23] - sorted[22]))
} #iqr30obs
es_outlier2 <- function(vect) {
start <- 1
end <- start + wl - 1
sorted <- sort(structure(vect[start:end], names=start:end))
i <- 0
res <- rep(NA, nrow(dt))
while (end < nrow(dt)) {
locFirstObs <- which(names(sorted)==start)
if (!(i > 9 && i < 22 && locFirstObs > 9 && locFirstObs < 22)) {
#changes in the 8th. 9th, 22th and 23th positions after removing first obs
#and adding new observation
qt <- iqr30obs(sorted)
iqr1.5 <- 1.5 * (qt[2] - qt[1])
}
res[end] <- sorted[as.character(end)] < qt[1] - iqr1.5 |
sorted[as.character(end)] > qt[2] + iqr1.5
#moving to next window ----
#remove the first observation in the window
sorted <- sorted[-locFirstObs]
#create the new observation to add to window
toAdd <- structure(vect[end+1], names=end+1)
#insert this new observation into the sorted vector while maintaining order
for (i in seq_along(sorted)) {
if (toAdd < sorted[i]) {
sorted <- c(sorted[seq_len(i-1)], toAdd, sorted[i:(wl-1)])
break
}
}
if (i == length(sorted)) {
sorted <- c(sorted, toAdd)
}
#increment indices
start <- start + 1
end <- end + 1
} #while
res
} #es_outlier2
es_outlier<-function(vect){
qq =quantile(vect, prob=c(0.25,0.75), na.rm=T)
q3=qq[2]
IC=q3-qq[1]
limSup=q3+IC*1.5
vector_final=abs(vect)>limSup
return(vector_final[length(vect)] )
}
结果:
system.time(es_outlier2(dt$x1))
# user system elapsed
# 4.62 0.00 4.67
system.time(es_outlier2(dt$x2))
# user system elapsed
# 4.56 0.00 4.83
system.time(zoo::rollapply(dt, 30, es_outlier, fill=NA, align='right'))
# user system elapsed
# 17.59 0.01 17.69