2013-07-30 202 views
3

我一直在編寫一個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也是如此,它似乎沒有改變任何東西。那麼區別是什麼呢?

非常感謝你提前, 吉爾達斯

+1

雖然沒有'NumericVector'的代碼,我猜想創建這麼多臨時對象可能是瓶頸。 – lapk

+2

將R函數傳遞給C++函數很慢。如果你希望事情變得快速,你需要在C++/Rcpp代碼中重新定義這些R函數。 –

回答

8

其他人已經回答了意見。所以我只強調一點:回調R函數的代價很高,因爲我們需要對錯誤處理格外小心。只用C++循環並調用R函數不會在C++中重寫您的代碼。嘗試將psipFGM改寫爲C++函數,並在此處向我們報告發生的情況。

你可能會認爲你失去了一些靈活性,你不能再使用任何R函數。對於這樣的情況,我建議使用某種混合解決方案,其中您已經在C++中實現了最常見的情況,否則將回退到R解決方案。

至於另一個問題,SEXP是一個R對象。這是R API的一部分。它可以是任何東西。當你從它創建一個Function(當創建一個函數需要一個參數Function時隱含地爲你完成),你可以保證這確實是一個R函數。開銷非常小,但代碼的表現力方面的收益是巨大的。

+7

+1 - 我們之前已經回答了類似的虛榮希望。如果僅僅在C++代碼片段中加入一個'R'函數可以加快速度,那麼我們都會一直這樣做。但是不可用午餐定理仍然成立...... –