4

特定の分布に関する積分を計算する 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)
}

正常に動作しますが、かなり遅いです。特に上記の R コードには二重ループが含まれているため、C++ などのコンパイル済み言語で関数を再プログラミングすると速度が向上するとよく耳にします。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)  
  }
))

それで、私にこれを説明してくれる親切な専門家はいますか？私は猿のようにコーディングしましたか?その機能を高速化する方法はありますか? さらに、2 つ目の質問があります。実際、出力型 double を SEXP に置き換え、引数型 Function を SEXP に置き換えることもできましたが、何も変わらないようです。違いは何ですか？

よろしくお願いします、ギルダス

4

1 に答える 1

13

他の人はすでにコメントで答えています。したがって、ここで強調しておきたいのは、R 関数のコールバックは、エラー処理に特に注意する必要があるため、コストがかかるということです。C++ でループして R 関数を呼び出すだけでは、C++ でコードを書き直すことにはなりません。psiandpFGMを C++ 関数として書き直してみて、何が起こるかをここに報告してください。

ある程度の柔軟性が失われ、R 関数を使用できなくなったと主張するかもしれません。このような状況では、最も一般的なケースを C++ で実装し、それ以外の場合は R ソリューションにフォールバックする、ある種のハイブリッド ソリューションを使用することをお勧めします。

他の質問については、 aSEXPは R オブジェクトです。これは R API の一部です。それは何でもかまいません。それからa を作成するFunctionと (引数をとる関数を作成するときに暗黙的に行われるようにFunction)、これが実際に R 関数であることが保証されます。オーバーヘッドは非常に小さいですが、コードの表現力の点で大きなメリットがあります。

于 2013-07-31T10:16:02.743 に答える