参考的例子是实现mvtnorm::dmvnorm()

代码如下：

<br />
#include < RcppArmadillo.h ><br />
// [[Rcpp::depends(RcppArmadillo)]]<br />
using namespace Rcpp ;<br />
const double log2pi = std::log(2.0 * M_PI);</p>
<p>// [[Rcpp::export]]<br />
arma::vec dmvnrm_arma1(arma::mat x,<br />
                      arma::rowvec mean,<br />
                      arma::mat sigma) {<br />
    int n = x.n_rows;<br />
    int xdim = x.n_cols;<br />
    arma::vec out(n);<br />
    arma::mat rooti = arma::trans(arma::inv(trimatu(arma::chol(sigma))));<br />
    double rootisum = arma::sum(log(rooti.diag()));<br />
    double constants = -(static_cast<double>(xdim)/2.0) * log2pi;</p>
<p>    for (int i=0; i < n; i++) {<br />
        arma::vec z = rooti * arma::trans( x.row(i) - mean) ;<br />
        out(i)      = constants - 0.5 * arma::sum(z%z) + rootisum;<br />
    }<br />
    out = exp(out);<br />
    return(out);<br />
}<br />
//[[Rcpp::export]]<br />
arma::rowvec dmvnrm_arma1(arma::mat X,arma::vec mu,arma::mat Sigma){<br />
   int k = X.n_cols;<br />
   int n = X.n_rows;<br />
   arma::mat rooti = arma::inv(arma::chol(Sigma));<br />
   arma::mat MU = repmat(mu,1,n);<br />
   arma::mat QUAD = (rooti.t()*(X.t()-MU)) % (rooti.t()*(X.t()-MU));<br />
   arma::rowvec quads = sum(QUAD);<br />
   arma::rowvec FD(n);<br />
   FD.fill(-(k/2)*std::log(2*M_PI)+arma::accu(arma::log(rooti.diag())));<br />
   return arma::exp(FD- 0.5*quads);<br />
}<br />
/*** R<br />
dMvn <- function(X,mu,Sigma) {<br />
  k <- ncol(X)<br />
  rooti <- backsolve(chol(Sigma),diag(k))<br />
  quads <- colSums((crossprod(rooti,(t(X)-mu)))^2)<br />
  return(exp(-(k/2)*log(2*pi) + sum(log(diag(rooti))) - .5*quads))<br />
}<br />
*/<br />

这个

dmvnrm_arma1

改写了

dmvnrm_arma

，有点向量化的味道，可是实验结果：

<br />
sigma <- bayesm::rwishart(10,diag(8))$IW<br />
means <- rnorm(8)<br />
X     <- mvtnorm::rmvnorm(900000, means, sigma)</p>
<p>benchmark(mvtnorm::dmvnorm(X,means,sigma,log=F),<br />
          dmvnrm_arma(X,means,sigma),<br />
          dmvnrm_arma1(X,means,sigma),<br />
          dMvn(X,means,sigma),<br />
          order="relative", replications=100)[,1:4]<br />
                                        test replications elapsed relative<br />
2               dmvnrm_arma(X, means, sigma)          100   40.51    1.000<br />
4                      dMvn(X, means, sigma)          100   53.96    1.332<br />
1 mvtnorm::dmvnorm(X, means, sigma, log = F)          100   68.99    1.703<br />
3              dmvnrm_arma1(X, means, sigma)          100  104.76    2.586</p>
<p>

<br />
 dmvnrm_arma1(X, means, sigma)#有点类似sugar语法，怎么变得这么慢。<br />

同样的使用

RcppEigen

<br />
#include < RcppEigen.h ><br />
using namespace Rcpp;<br />
// [[Rcpp::depends(RcppEigen)]]<br />
//[[Rcpp::export]]<br />
using Eigen::Map;<br />
using Eigen::MatrixXd;<br />
using Eigen::VectorXd;<br />
using Eigen::LLT;<br />
VectorXd dmvnrm_Eigen(Map<MatrixXd> X,Map<VectorXd> mu,Map<MatrixXd> Sigma){<br />
  const int k = X.cols();<br />
   int n = X.rows();<br />
   MatrixXd D(k,k);<br />
   LLT<MatrixXd>lltOfsigma(Sigma);<br />
   MatrixXd L = lltOfsigma.matrixL();<br />
   MatrixXd rooti = L.inverse();<br />
   MatrixXd MU = mu.replicate(1,n);<br />
   MatrixXd Q = rooti.transpose()*(X.transpose()-MU);<br />
   MatrixXd QUAD = Q.cwiseProduct(Q);<br />
   VectorXd quads = QUAD.colwise().sum();<br />
   VectorXd FD;<br />
   double Con = -(k/2)*std::log(2*M_PI)+(rooti.diagonal().log().sum());<br />
   FD.setConstant(1, n, Con);<br />
   return (FD - 0.5*quads).exp();<br />
}<br />

运算时间基本上呵呵了。

看来是代码水平太着急啦。
</p>