On 29 October 2014 at 16:08, 李晓辉 wrote: | 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 can't replicate that. | I do not know what is wrong with | double Con = -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); For me it also segfaults when I comment out the vector operation and use the loop, suggesting your error is somewhere else. The Rcpp Gallery has a working dmvnorm using RcppArmadillo at http://gallery.rcpp.org/articles/dmvnorm_arma/ If you want to use RcppEigen, I suggest you study the Eigen documentation carefully and examine your program, possibly with the help of a debugger. Dirk -- http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org _______________________________________________ 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