我一直在編寫一個R函數來計算關於某些分佈的積分,參見下面的代碼。Rcpp功能比相同R功能
EVofPsi = function(psi, probabilityMeasure, eps=0.01, ...){
distFun = function(u){
probabilityMeasure(u, ...)
}
xx = yy = seq(0,1,length=1/eps+1)
summand=0
for(i in 1:(length(xx)-1)){
for(j in 1:(length(yy)-1)){
signPlus = distFun(c(xx[i+1],yy[j+1]))+distFun(c(xx[i],yy[j]))
signMinus = distFun(c(xx[i+1],yy[j]))+distFun(c(xx[i],yy[j+1]))
summand = c(summand, psi(c(xx[i],yy[j]))*(signPlus-signMinus))
}
}
sum(summand)
}
它工作正常,但它是相當緩慢的。通常聽說用C++等編譯語言對函數進行重新編程會加快速度,特別是因爲上面的R代碼涉及雙循環。我也是,使用Rcpp:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double EVofPsiCPP(Function distFun, Function psi, int n, double eps) {
NumericVector xx(n+1);
NumericVector yy(n+1);
xx[0] = 0;
yy[0] = 0;
// discretize [0,1]^2
for(int i = 1; i < n+1; i++) {
xx[i] = xx[i-1] + eps;
yy[i] = yy[i-1] + eps;
}
Function psiCPP(psi);
Function distFunCPP(distFun);
double signPlus;
double signMinus;
double summand = 0;
NumericVector topRight(2);
NumericVector bottomLeft(2);
NumericVector bottomRight(2);
NumericVector topLeft(2);
// compute the integral
for(int i=0; i<n; i++){
//printf("i:%d \n",i);
for(int j=0; j<n; j++){
//printf("j:%d \n",j);
topRight[0] = xx[i+1];
topRight[1] = yy[j+1];
bottomLeft[0] = xx[i];
bottomLeft[1] = yy[j];
bottomRight[0] = xx[i+1];
bottomRight[1] = yy[j];
topLeft[0] = xx[i];
topLeft[1] = yy[j+1];
signPlus = NumericVector(distFunCPP(topRight))[0] + NumericVector(distFunCPP(bottomLeft))[0];
signMinus = NumericVector(distFunCPP(bottomRight))[0] + NumericVector(distFunCPP(topLeft))[0];
summand = summand + NumericVector(psiCPP(bottomLeft))[0]*(signPlus-signMinus);
//printf("summand:%f \n",summand);
}
}
return summand;
}
我很開心,因爲這個C++函數可以正常工作。然而,當我測試了這兩種功能,C++的一個跑慢:
sourceCpp("EVofPsiCPP.cpp")
pFGM = function(u,theta){
u[1]*u[2] + theta*u[1]*u[2]*(1-u[1])*(1-u[2])
}
psi = function(u){
u[1]*u[2]
}
print(system.time(
for(i in 1:10){
test = EVofPsi(psi, pFGM, 1/100, 0.2)
}
))
test
print(system.time(
for(i in 1:10){
test = EVofPsiCPP(psi, function(u){pFGM(u,0.2)}, 100, 1/100)
}
))
那麼,有沒有某種專家圍繞願意解釋我這個?我是否像猴子一樣編碼,並且有辦法加快這一功能?此外,我還有第二個問題。事實上,我可以用SEXP替換輸出類型double,並且參數類型Function by SEXP也是如此,它似乎沒有改變任何東西。那麼區別是什麼呢?
非常感謝你提前, 吉爾達斯
雖然沒有'NumericVector'的代碼,我猜想創建這麼多臨時對象可能是瓶頸。 – lapk
將R函數傳遞給C++函數很慢。如果你希望事情變得快速,你需要在C++/Rcpp代碼中重新定義這些R函數。 –