<div dir="ltr"><div><div><div><div>Hi Xiaohui,<br></div>There are several problems in the code I have found:<br>1. <font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">MatrixXd rooti = L.inverse()*Eigen::MatrixXd::Identity(k,k) is not necessary. You are multiplying a matrix with an identity matrix, so it can be simplified as<br></font><br><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><span style="font-family:courier new,monospace">MatrixXd rooti = L.inverse();</span><br><br></font></div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">2. k is an integer, so you need to use k / 2.0 in calculating Con.<br><br></font></div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">3. Most importantly, the reason to cause segfault is the misuse of .log() and .exp() in the </font><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">"double Con = " and "return" statements. When you want to do coefficient-wise operations on a vector or a matrix, you should first convert it to Array class. So these two lines should be<br><br><span style="font-family:courier new,monospace">double Con = -(k/2.0)*std::log(2*M_PI)+rooti.diagonal().array().log().sum();<br></span><br></font></div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">and<br><br><span style="font-family:courier new,monospace">return (FD - 0.5*quads).array().exp();</span><br></font><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><br></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">4. To correctly computing the quadratic terms, you should write<br><br><span style="font-family:courier new,monospace">MatrixXd Q = rooti*(X.transpose()-MU);</span><br><br></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">instead of<br></font><br><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><span style="font-family:courier new,monospace">MatrixXd Q = rooti.transpose()*(X.</span><span style="font-family:courier new,monospace">transpose()-MU);</span><br><br><br><br></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">Best,<br></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif">Yixuan<br></font></div><div><font color="#333333" face="Helvetica Neue, Helvetica, Arial, sans-serif"><br></font></div></div><div class="gmail_extra"><br><div class="gmail_quote">2014-10-29 4:08 GMT-04:00 李晓辉 <span dir="ltr"><<a href="mailto:2009lxh@gmail.com" target="_blank">2009lxh@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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/." target="_blank">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>
<br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br><br clear="all"><br>-- <br>Yixuan Qiu <<a href="mailto:yixuan.qiu@cos.name" target="_blank">yixuan.qiu@cos.name</a>><br>Department of Statistics,<br>Purdue University<br>
</div>