まえから気になっていた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倍ぐらいの時間がかかった....
何がボトルネックなんだろう...
状況はかなり悪く、山域ほぼ食い荒らされている状況。
どこもかしこも....
