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
