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和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)
}

现在解释一下rollmean_cpp.x并且y是数据.xout是请求滚动统计的点向量.width是滚动窗口的宽度*2.请注意,滑动窗口末端的indeces存储在ledge和中redge.这些基本上是指针到相应的元件x和y.这些indeces可能非常有利于调用其他C++函数(例如,中位数等),这些函数将向量以及开始和结束indeces作为输入.

对于那些想要"冗长"版本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