Rcpp+inlineを使ってみた

まえから気になっていたRcpp+inlineを使ってみた

単純な繰り返しループの集計計算がかなり時間がかかるので、Rcppでどれくらい改善されるか挑戦

まずは通常のコード

非常に単純で以下のxとBeta,Alphaの線形結合のquantileを求めるというもの。

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
res1 <- sapply(1:length(x), function(i) 
               quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025))

次にRcpp+inline

library(Rcpp)
library(inline)
src <- "NumericVector xx(x);
  NumericVector beta(b);
  NumericVector alpha(a);
  Function quantile(\"quantile\");
  int xsize=xx.size(), bsize=beta.size();
  NumericVector mu(bsize);
  NumericVector res(xsize);
  for(unsigned long int i=0; i < xsize; i++){
   for(unsigned long int j=0; j < bsize; j++){
    double dummy=exp(beta[j]*xx[i]+alpha[j]);
    mu[j]=0.0;
    mu[j]=dummy/(1.0+dummy);
   }
   res[i]=0.0;
   res[i]=as(quantile(mu, _[\"probs\"]=0.025));
  }
  return res;"

fun.R <- cxxfunction(signature(x="numeric",a="numeric",b="numeric"),
       src,plugin="Rcpp")

結果

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
system.time(res1 <- sapply(1:length(x), function(i) 
    quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025)))
   ユーザ   システム       経過  
    19.162      1.008     20.169 
system.time(res2 <- fun.R(x,Beta,Alpha))
   ユーザ   システム       経過  
     8.936      0.932      9.870 

計算結果も同じだった

summary(res1-res2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

RcppでRの関数を使うのではなく、GNU Scientific Library (GSL)のquantile(gsl_stats_quantile_from_sorted_data)も試してみたけどこちらは通常のRの1.5倍ぐらいの時間がかかった....
何がボトルネックなんだろう...