<div dir="ltr">These days i rewrote the <span style="color:rgb(51,51,51);font-family:'Helvetica Neue',Helvetica,Arial,sans-serif;font-size:13px;line-height:18px"> </span><code style="padding:2px 4px;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12px;color:rgb(221,17,68);border-top-left-radius:3px;border-top-right-radius:3px;border-bottom-right-radius:3px;border-bottom-left-radius:3px;border:1px solid rgb(225,225,232);line-height:18px;background-color:rgb(247,247,249)">dmvnorm()</code><span style="color:rgb(51,51,51);font-family:'Helvetica Neue',Helvetica,Arial,sans-serif;font-size:13px;line-height:18px"> function from the </span><code style="padding:2px 4px;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12px;color:rgb(221,17,68);border-top-left-radius:3px;border-top-right-radius:3px;border-bottom-right-radius:3px;border-bottom-left-radius:3px;border:1px solid rgb(225,225,232);line-height:18px;background-color:rgb(247,247,249)">mvtnorm</code><span style="color:rgb(51,51,51);font-family:'Helvetica Neue',Helvetica,Arial,sans-serif;font-size:13px;line-height:18px"> package with RcppEigen which is publised in <a href="http://gallery.rcpp.org/articles/dmvnorm_arma/.">http://gallery.rcpp.org/articles/dmvnorm_arma/.</a></span><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><span style="line-height:18px">The code as follows:</span></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><blockquote style="line-height:18px;margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex" class="gmail_quote">#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()*Eigen::MatrixXd::Identity(k,k);<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(n);<br>   double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum();<br>   FD.setConstant(Con);<br>   return (FD - 0.5*quads).exp();<br>}</blockquote><blockquote style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex" class="gmail_quote"><span style="line-height:18px">/*** R<br></span><span style="line-height:18px">sigma <- bayesm::rwishart(10,diag(8))$IW<br></span><span style="line-height:18px">means <- rnorm(8)<br></span><span style="line-height:18px">X     <- mvtnorm::rmvnorm(90, means, sigma)<br></span><span style="line-height:18px">dmvnrm_Eigen(X,means,sigma)<br></span><span style="line-height:18px">*/ </span></blockquote><div style="line-height:18px"><br></div><div style="line-height:18px">When i run above code with RStudio, it often got an error message <b>:R Session Aborted</b></div><div style="line-height:18px"> it seems that<b><u> double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); </u></b>got a problem.Actually,if i rewrite these code as follow:</div><div style="line-height:18px"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">double Con=0.;<br>   for (int i=0;i<k;i++){<br>   Con+=log(rooti(i,i));<br>   }<br>Con+=-(k/2)*std::log(2*M_PI)</blockquote></div><div style="line-height:18px"><br></div><div style="line-height:18px">it work well.</div><div style="line-height:18px">I do not know what is wrong with </div><div style="line-height:18px"><b><u>double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum();</u></b><br></div><div style="line-height:18px"><b><u><br></u></b></div></font></div></div>