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。请注意，滑动窗口末端的索引存储在ledge 和redge 中。这些本质上是指向x 和y 中各自元素的指针。这些索引对于调用其他将向量以及开始和结束索引作为输入的 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