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

| 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

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 ??? <2009lxh at 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 at lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

Yixuan Qiu <yixuan.qiu at cos.name>
Department of Statistics,
Purdue University
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141029/6a3189d3/attachment-0001.html>