Hi Xiaohui, There are several problems in the code I have found: 1. 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
MatrixXd rooti = L.inverse(); 2. k is an integer, so you need to use k / 2.0 in calculating Con. 3. Most importantly, the reason to cause segfault is the misuse of .log() and .exp() in the "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 double Con = -(k/2.0)*std::log(2*M_PI)+rooti.diagonal().array().log().sum(); and return (FD - 0.5*quads).array().exp(); 4. To correctly computing the quadratic terms, you should write MatrixXd Q = rooti*(X.transpose()-MU); instead of MatrixXd Q = rooti.transpose()*(X.transpose()-MU); Best, Yixuan 2014-10-29 4:08 GMT-04:00 李晓辉 <2009...@gmail.com>: > These days i rewrote the dmvnorm() function from the mvtnorm package > with RcppEigen which is publised in > http://gallery.rcpp.org/articles/dmvnorm_arma/. > The code as follows: > >> #include <RcppEigen.h> >> using namespace Rcpp; >> // [[Rcpp::depends(RcppEigen)]] >> //[[Rcpp::export]] >> using Eigen::Map; >> using Eigen::MatrixXd; >> using Eigen::VectorXd; >> using Eigen::LLT; >> VectorXd dmvnrm_Eigen(Map<MatrixXd> X,Map<VectorXd> mu,Map<MatrixXd> >> Sigma){ >> const int k = X.cols(); >> int n = X.rows(); >> MatrixXd D(k,k); >> LLT<MatrixXd>lltOfsigma(Sigma); >> MatrixXd L = lltOfsigma.matrixL(); >> MatrixXd rooti = L.inverse()*Eigen::MatrixXd::Identity(k,k); >> MatrixXd MU = mu.replicate(1,n); >> MatrixXd Q = rooti.transpose()*(X.transpose()-MU); >> MatrixXd QUAD = Q.cwiseProduct(Q); >> VectorXd quads = QUAD.colwise().sum(); >> VectorXd FD(n); >> double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); >> FD.setConstant(Con); >> return (FD - 0.5*quads).exp(); >> } > > /*** R >> sigma <- bayesm::rwishart(10,diag(8))$IW >> means <- rnorm(8) >> X <- mvtnorm::rmvnorm(90, means, sigma) >> dmvnrm_Eigen(X,means,sigma) >> */ > > > When i run above code with RStudio, it often got an error message *:R > Session Aborted* > it seems that* double Con = > -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); *got a > problem.Actually,if i rewrite these code as follow: > >> double Con=0.; >> for (int i=0;i<k;i++){ >> Con+=log(rooti(i,i)); >> } >> Con+=-(k/2)*std::log(2*M_PI) > > > it work well. > I do not know what is wrong with > *double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum();* > > > _______________________________________________ > Rcpp-devel mailing list > Rcpp-devel@lists.r-forge.r-project.org > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > -- Yixuan Qiu <yixuan....@cos.name> Department of Statistics, Purdue University
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
