преобразование цикла из R в C ++ с использованием Rcpp

Question

Я хочу повысить скорость некоторых моих R-кодов с помощью Rcpp.Тем не менее, мои знания C ++ очень мало.Итак, я проверил документацию, предоставленную Rcpp, и другие документы, предоставленные на сайте Дирка Эддельбюттеля.Прочитав все это, я попытался выполнить простой цикл, который написал в 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.

Answer 1

Базовый пример 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)

, и вы можете просмотреть документ броненосца , чтобы найти более продвинутые биты и кусочки.

Answer 2

ответ на мой первый вопрос показан ниже. Это может быть не эффективным способом, но код 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)