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

[Yixuan Qiu]

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>