R - 가변 간격으로 롤링 통계를 계산하는 더 빠른 방법

Question

변수를 통해 롤링 통계 (롤링 평균, 중간 값, 백분위 수 등)를 계산하는 방법이 생길지 궁금합니다. 시간 간격 (윈도우 잉).

즉, 무작위 적으로 시간이 지정된 관찰 (즉, 일별 데이터 또는 주간 데이터가 아닌 관측치에는 틱 데이터와 같이 시간 스탬프가있는 것)이 주어 졌다고 가정하고 그 중앙 및 분산 통계를보고 싶다고 가정 해 봅시다. 이 통계가 계산되는 시간 간격을 넓히고 조이면됩니다.

이 작업을 수행하는 간단한 for 루프를 만들었습니다. 하지만 분명히 매우 느리게 실행됩니다 (실제로 루프가 속도를 테스트하기 위해 설정 한 작은 데이터 샘플을 통해 계속 실행되고 있다고 생각합니다). 나는 이것을하기 위해 ddply와 같은 것을 얻으려고 노력해왔다. - 그것은 일일 통계를 위해 달리는 것이 곤란한 것처럼 보이지만, 나는 그것을 벗어날 수 없다.

예 :

샘플 세트 업!

df <- data.frame(Date = runif(1000,0,30)) 
df$Price <- I((df$Date)^0.5 * (rnorm(1000,30,4))) 
df$Date <- as.Date(df$Date, origin = "1970-01-01") 

예 기능 (즉

SummaryStats <- function(dataframe, interval){ 
    # Returns daily simple summary stats, 
    # at varying intervals 
    # dataframe is the data frame in question, with Date and Price obs 
    # interval is the width of time to be treated as a day 

    firstDay <- min(dataframe$Date) 
    lastDay <- max(dataframe$Date) 
    result <- data.frame(Date = NULL, 
         Average = NULL, Median = NULL, 
         Count = NULL, 
         Percentile25 = NULL, Percentile75 = NULL) 

    for (Day in firstDay:lastDay){ 

    dataframe.sub = subset(dataframe, 
       Date > (Day - (interval/2)) 
       & Date < (Day + (interval/2))) 

    nu = data.frame(Date = Day, 
        Average = mean(dataframe.sub$Price), 
        Median = median(dataframe.sub$Price), 
        Count = length(dataframe.sub$Price), 
        P25 = quantile(dataframe.sub$Price, 0.25), 
        P75 = quantile(dataframe.sub$Price, 0.75)) 

    result = rbind(result,nu) 

    } 

    return(result) 

} 

귀하의 조언은 환영받을 많은 관찰 정말 느린 실행

Answer 1

Rcpp 속도가 가장 중요한 문제인 경우 좋은 방법입니다. 롤링 평균 통계를 사용하여 예제로 설명하겠습니다.

벤치 마크 : Rcpp R

x = sort(runif(25000,0,4*pi)) 
y = sin(x) + rnorm(length(x),0.5,0.5) 
system.time(rollmean_r(x,y,xout=x,width=1.1)) # ~60 seconds 
system.time(rollmean_cpp(x,y,xout=x,width=1.1)) # ~0.0007 seconds 

코드 Rcpp 및 rollmean_cpp의 explantion 지금 R 기능

cppFunction(' 
    NumericVector rollmean_cpp(NumericVector x, NumericVector y, 
           NumericVector xout, double width) { 
    double total=0; 
    unsigned int n=x.size(), nout=xout.size(), i, ledge=0, redge=0; 
    NumericVector out(nout); 

    for(i=0; i<nout; i++) { 
     while(x[ redge ] - xout[i] <= width && redge<n) 
     total += y[redge++]; 
     while(xout[i] - x[ ledge ] > width && ledge<n) 
     total -= y[ledge++]; 
     if(ledge==redge) { out[i]=NAN; total=0; continue; } 
     out[i] = total/(redge-ledge); 
    } 
    return out; 
    }') 

rollmean_r = function(x,y,xout,width) { 
    out = numeric(length(xout)) 
    for(i in seq_along(xout)) { 
    window = x >= (xout[i]-width) & x <= (xout[i]+width) 
    out[i] = .Internal(mean(y[window])) 
    } 
    return(out) 
} 

대. x 및 y이 데이터입니다. xout은 롤링 통계가 요청되는 지점의 벡터입니다. width은 롤링 창의 너비 * 2입니다. 슬라이딩 윈도우의 끝 부분에 대한 indeces는 ledge 및 redge에 저장됩니다. 이것들은 기본적으로 x과 y의 각 요소에 대한 포인터입니다. 이 indeces는 벡터를 취하고 indents를 입력으로 시작하고 끝내는 다른 C++ 함수 (예 : 중앙값 등)를 호출하는 데 매우 유용 할 수 있습니다. (긴) 디버깅을 위해 rollmean_cpp의 "자세한"버전을 원하는 사람들을 위해

: 내 질문에 "케빈"에 응답에서

cppFunction(' 
    NumericVector rollmean_cpp(NumericVector x, NumericVector y, 
           NumericVector xout, double width) { 

    double total=0, oldtotal=0; 
    unsigned int n=x.size(), nout=xout.size(), i, ledge=0, redge=0; 
    NumericVector out(nout); 


    for(i=0; i<nout; i++) { 
     Rcout << "Finding window "<< i << " for x=" << xout[i] << "..." << std::endl; 
     total = 0; 

     // numbers to push into window 
     while(x[ redge ] - xout[i] <= width && redge<n) { 
     Rcout << "Adding (x,y) = (" << x[redge] << "," << y[redge] << ")" ; 
     Rcout << "; edges=[" << ledge << "," << redge << "]" << std::endl; 
     total += y[redge++]; 
     } 

     // numbers to pop off window 
     while(xout[i] - x[ ledge ] > width && ledge<n) { 
     Rcout << "Removing (x,y) = (" << x[ledge] << "," << y[ledge] << ")"; 
     Rcout << "; edges=[" << ledge+1 << "," << redge-1 << "]" << std::endl; 
     total -= y[ledge++]; 
     } 
     if(ledge==n) Rcout << " OVER "; 
     if(ledge==redge) { 
     Rcout<<" NO DATA IN INTERVAL " << std::endl << std::endl; 
     oldtotal=total=0; out[i]=NAN; continue;} 

     Rcout << "For interval [" << xout[i]-width << "," << 
       xout[i]+width << "], all points in interval [" << x[ledge] << 
       ", " << x[redge-1] << "]" << std::endl ; 
     Rcout << std::endl; 

     out[i] = (oldtotal + total)/(redge-ledge); 
     oldtotal=total+oldtotal; 
    } 
    return out; 
    }') 

x = c(1,2,3,6,90,91) 
y = c(9,8,7,5.2,2,1) 
xout = c(1,2,2,3,6,6.1,13,90,100) 
a = rollmean_cpp(x,y,xout=xout,2) 
# Finding window 0 for x=1... 
# Adding (x,y) = (1,9); edges=[0,0] 
# Adding (x,y) = (2,8); edges=[0,1] 
# Adding (x,y) = (3,7); edges=[0,2] 
# For interval [-1,3], all points in interval [1, 3] 
# 
# Finding window 1 for x=2... 
# For interval [0,4], all points in interval [1, 3] 
# 
# Finding window 2 for x=2... 
# For interval [0,4], all points in interval [1, 3] 
# 
# Finding window 3 for x=3... 
# For interval [1,5], all points in interval [1, 3] 
# 
# Finding window 4 for x=6... 
# Adding (x,y) = (6,5.2); edges=[0,3] 
# Removing (x,y) = (1,9); edges=[1,3] 
# Removing (x,y) = (2,8); edges=[2,3] 
# Removing (x,y) = (3,7); edges=[3,3] 
# For interval [4,8], all points in interval [6, 6] 
# 
# Finding window 5 for x=6.1... 
# For interval [4.1,8.1], all points in interval [6, 6] 
# 
# Finding window 6 for x=13... 
# Removing (x,y) = (6,5.2); edges=[4,3] 
# NO DATA IN INTERVAL 
# 
# Finding window 7 for x=90... 
# Adding (x,y) = (90,2); edges=[4,4] 
# Adding (x,y) = (91,1); edges=[4,5] 
# For interval [88,92], all points in interval [90, 91] 
# 
# Finding window 8 for x=100... 
# Removing (x,y) = (90,2); edges=[5,5] 
# Removing (x,y) = (91,1); edges=[6,5] 
# OVER NO DATA IN INTERVAL 

print(a) 
# [1] 8.0 8.0 8.0 8.0 5.2 5.2 NaN 1.5 NaN 

Answer 2

자 ... 루프 (매우 느리게 R)를하고 있습니다. 서브 세트 생성시 불필요한 데이터 사본을 만들고 rbind을 사용하여 데이터 세트를 축적합니다. 이를 피하면 일이 상당히 빨라질 것입니다. 이것을 시도하십시오 ...

Summary_Stats <- function(Day, dataframe, interval){ 
    c1 <- dataframe$Date > Day - interval/2 & 
     dataframe$Date < Day + interval/2 
    c(
     as.numeric(Day), 
     mean(dataframe$Price[c1]), 
     median(dataframe$Price[c1]), 
     sum(c1), 
     quantile(dataframe$Price[c1], 0.25), 
     quantile(dataframe$Price[c1], 0.75) 
    ) 
} 
Summary_Stats(df$Date[2],dataframe=df, interval=20) 
firstDay <- min(df$Date) 
lastDay <- max(df$Date) 
system.time({ 
    x <- sapply(firstDay:lastDay, Summary_Stats, dataframe=df, interval=20) 
    x <- as.data.frame(t(x)) 
    names(x) <- c("Date","Average","Median","Count","P25","P75") 
    x$Date <- as.Date(x$Date) 
}) 
dim(x) 
head(x) 

Answer 3

위, 나는 밖으로 아래 뭔가를 생각 생각합니다.

이 함수는 틱 데이터 (시간 관측치가 무작위 간격으로 나타나고 타임 스탬프로 표시됨)를 취하여 일정 간격 동안 평균을 계산합니다. 사슬로 연결된 모든 포인트의

library(Rcpp) 

cppFunction(' 
    NumericVector rollmean_c2(NumericVector x, NumericVector y, double width, 
           double Min, double Max) { 

double total = 0, redge,center; 
unsigned int n = (Max - Min) + 1, 
        i, j=0, k, ledge=0, redgeIndex; 
NumericVector out(n); 


for (i = 0; i < n; i++){ 
    center = Min + i + 0.5; 
    redge = center - width/2; 
    redgeIndex = 0; 
    total = 0; 

    while (x[redgeIndex] < redge){ 
    redgeIndex++; 
    } 
    j = redgeIndex; 

    while (x[j] < redge + width){ 
    total += y[j++]; 

    } 

    out[i] = total/(j - redgeIndex); 
} 
return out; 

    }') 

# Set up example data 
x = seq(0,4*pi,length.out=2500) 
y = sin(x) + rnorm(length(x),0.5,0.5) 
plot(x,y,pch=20,col="black", 
    main="Sliding window mean; width=1", 
    sub="rollmean_c in red  rollmean_r overlaid in white.") 


c.out = rollmean_c2(x,y,width=1,Min = min(x), Max = max(x)) 
lines(0.5:12.5,c.out,col="red",lwd=3) 

Answer 4

생각한다. 이 체인을 각 데이터 포인트가 노드 인 그래프로 생각하십시오. 그런 다음 각 노드에 대해 거리가 w 또는 그 이하인 다른 모든 노드를 찾고자합니다. 이렇게하기 위해서, 먼저 pairwise 거리를주는 행렬을 생성합니다. n 번째 행은 노드 n 노드의 거리를 구분합니다.

# First, some data 
x = sort(runif(25000,0,4*pi)) 
y = sin(x) + rnorm(length(x),0,0.5) 

# calculate the rows of the matrix one by one 
# until the distance between the two closest nodes is greater than w 
# This algorithm is actually faster than `dist` because it usually stops 
# much sooner 
dl = list() 
dl[[1]] = diff(x) 
i = 1 
while(min(dl[[i]]) <= w) { 
    pdl = dl[[i]] 
    dl[[i+1]] = pdl[-length(pdl)] + dl[[1]][-(1:i)] 
    i = i+1 
} 

# turn the list of the rows into matrices 
rarray = do.call(rbind, lapply(dl,inf.pad,length(x))) 
larray = do.call(rbind, lapply(dl,inf.pad,length(x),"right")) 

# extra function 
inf.pad = function(x,size,side="left") { 
    if(side=="left") { 
    x = c(x, rep(Inf, size-length(x))) 
    } else { 
    x = c(rep(Inf, size-length(x)), x) 
    } 
    x 
} 

그런 다음 매트릭스를 사용하여 각 창의 가장자리를 결정하십시오. 이 예에서는 w=2을 설정합니다. 정의 된 창으로

# How many data points to look left or right at each data point 
lookr = colSums(rarray <= w) 
lookl = colSums(larray <= w) 

# convert these "look" variables to indeces of the input vector 
ri = 1:length(x) + lookr 
li = 1:length(x) - lookl 

, 그것은 최종 답변을 얻을 수있는 *apply 기능을 사용하기 매우 간단합니다.

rolling.mean = vapply(mapply(':',li,ri), function(i) .Internal(mean(y[i])), 1) 

위의 코드는 모두 내 컴퓨터에서 약 50 초가 걸렸습니다. 이것은 내 대답에 rollmean_r 함수보다 약간 빠릅니다. 그러나 여기서 특히 좋은 점은 indeces가 제공된다는 것입니다. 그런 다음 *apply 함수를 사용하여 원하는 R 함수를 사용할 수 있습니다. 예를 들어,

rolling.mean = vapply(mapply(':',li,ri), 
             function(i) .Internal(mean(y[i])), 1) 

은 약 5 초가 걸립니다. 그리고,

rolling.median = vapply(mapply(':',li,ri), 
             function(i) median(y[i]), 1) 

이 약 14 초 걸린다. 원한다면 다른 답변에서 Rcpp 함수를 사용하여 indeces를 얻을 수 있습니다.