0

Rcpp を使用して一部の R コードの速度を改善したいと考えています。ただし、C++ に関する私の知識はほとんどありません。そこで、Rcpp で提供されているドキュメントと、Dirk Eddelbuttel のサイトで提供されているその他のドキュメントを確認しました。すべてを読んだ後、R で書いた単純なループを実行しようとしましたが、残念ながら実行できませんでした。R関数は次のとおりです。

Inverse Wishart

beta = matrix(rnorm(15),ncol=3)

a = rnorm(3) 

InW = function(beta,a) {

    n = nrow(beta)
    p = ncol(beta)
    I = diag(rep(1,times = p))
    H = matrix(0,nrow=p,ncol=p)
    for(i in 1:n){
    subBi = beta[i,]
          H = H + tcrossprod(a - subBi)
        }
    H = H + p * I

    T = t(chol(chol2inv(chol(H))))
    S = 0
    for(i in 1:(n+p)){
        u <- rnorm(p)
        S = S + tcrossprod(T %*% u)
        }
    D = chol2inv(chol((S)))
    ans = list(Dinv = S,D=D)
}

Rcppの学習の出発点として役立つので、誰かが私を助けてくれれば本当に感謝しています.

4

2 に答える 2

5

RcppArmadillo の基本的な例は次のようになります。

require(RcppArmadillo)
require(inline)

code <- '
  arma::mat beta = Rcpp::as<arma::mat>(beta_);
  int n = beta.n_rows; int p = beta.n_cols;
  arma::mat Ip = arma::eye<arma::mat>( p, p );
  int ii;
  double S=0;
  for (ii=0; ii<(n+p); ii++) {
    S += ii; // dummy calculation
  }
  return Rcpp::wrap(S);
 '

fun <- cxxfunction(signature(beta_ ="matrix"),
                       code, plugin="RcppArmadillo")

m <- matrix(1:9,3)
fun(m)

また、アルマジロのドキュメントを参照して、より高度な部分を見つけることができます。

于 2012-02-06T06:06:10.067 に答える
1

私の最初の質問に対する答えを以下に示します。効率的な方法ではないかもしれませんが、Rcpp コードは R コードと同じ結果を返します。バティストの助けに感謝します。

code <- '<br/>
arma::mat beta = Rcpp::as<arma::mat>(beta_);
arma::rowvec y = Rcpp::as<arma::rowvec>(y_);
int n = beta.n_rows; int p = beta.n_cols;
arma::mat Ip = arma::eye<arma::mat>( p, p );
int ii;
arma::mat H1 = beta,  d;
arma::mat H2=H1.zeros(p,p);
arma::rowvec S;
for (ii=0;ii<n;ii++){
S= beta.row(ii);
d = trans(y - S)*(y-S);
H2 = H2 + d ;
}
arma::mat H = chol(H2+p*Ip);
arma::mat Q , R;
 qr(Q,R,H);
arma::mat RR = R;
arma::mat TT = trans(chol(solve(trans(RR)*RR,Ip)));
int jj;
arma::mat SS = H1.zeros(p,p);
arma::colvec u;
arma::colvec V;
for(jj=0;jj<(n+p);jj++) {
       u = rnorm(p);
       V = TT*u;
      SS = SS + V * trans(V);
      }
arma::mat SS1 = chol(SS);
arma::mat Q1 , R1;
qr(Q1,R1,SS1);
arma::mat SS2 = R1;
arma::mat D = solve(trans(SS2)*SS2,Ip);
return Rcpp::List::create(Rcpp::Named("Dinv")=SS,Rcpp::Named("D")=D);
'
fun = cxxfunction(signature(beta_ ="matrix",y_="numeric"),code, plugin="RcppArmadillo")
m = matrix(rnorm(100),ncol=5)
vec = rnorm(5)
fun(m,vec)
于 2012-02-23T06:10:51.860 に答える