[Rcpp-devel] The problem i encountered about RcppEigen double -(k/2)*std::log(2*M_PI)+rooti.diagonal().log().sum(); can not run .

[Dirk Eddelbuettel]

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 | edd at debian.org