参考的例子是实现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>