分开写或者不分开写我的windows下有溢出(Revolution6),我是这么改的:
<br />
/*File: C_dmvnrm_Eigen.cpp*/<br />
#include R.h<br />
#include Rcpp.h<br />
#include RcppEigen.h</p>
<p> RcppExport SEXP C_dmvnrm_Eigen(SEXP X_,SEXP mu_,SEXP Sigma_){<br />
using Eigen::Map;<br />
using Eigen::MatrixXd;<br />
using Eigen::VectorXd;<br />
using Eigen::LLT;<br />
typedef Eigen::Map<Eigen::MatrixXd> MapMatd;<br />
typedef Eigen::Map<Eigen::VectorXd> MapVecd;<br />
const MapMatd X(Rcpp::as<MapMatd>(X_));<br />
const MapVecd mu(Rcpp::as<MapVecd>(mu_));<br />
const MapMatd Sigma(Rcpp::as<MapMatd>(Sigma_));<br />
int k = X.cols();<br />
int n = X.rows();<br />
MatrixXd D(k,k);<br />
LLT<MatrixXd>lltOfsigma(Sigma);<br />
MatrixXd L = lltOfsigma.matrixL();<br />
MatrixXd rooti = L.inverse()*Eigen::MatrixXd::Identity(k,k);<br />
double Con=0.;<br />
for (int i=0;i<k;i++){<br />
Con+=log(rooti(i,i));<br />
Rprintf("[%d,%d]:=%f\t",i,i,rooti(i,i));<br />
}<br />
Rprintf("\nRe:\t%f\n",Con);</p>
<p> return(Rcpp::wrap(Con));<br />
}</p>
<p>/*** R<br />
sigma <- bayesm::rwishart(10,diag(8))$IW<br />
means <- rnorm(8)<br />
X <- mvtnorm::rmvnorm(90, means, sigma)<br />
dmvnrm_Eigen(X,means,sigma)<br />
*/<br />
</p>
Eigen 下
rooti.diagonal() = array
和
array=rooti.diagonal()
的语法都是允许的,可能这就是问题吧。
<br />
dyn.load("C_dmvnrm_Eigen.dll")<br />
library(rbenchmark)</p>
<p>sigma <- as.matrix(<br />
bayesm::rwishart(10,diag(8))$IW<br />
)<br />
means <- rnorm(8)<br />
library(RcppEigen)</p>
<p>X <- as.matrix(mvtnorm::rmvnorm(90, means, sigma))</p>
<p>.Call("C_dmvnrm_Eigen",X,means,sigma)<br />
dyn.unload("C_dmvnrm_Eigen.dll")</p>
<p>> .Call("C_dmvnrm_Eigen",X,means,sigma)<br />
[0,0]:=1.887568 [1,1]:=1.976924 [2,2]:=2.268791 [3,3]:=2.510085 [4,4]:=2.274426 [5,5]:=3.133264 [6,6]:=3.331032 [7,7]:=3.274877<br />
Re: 7.409761<br />
[1] 7.409761<br />
> dyn.unload("C_dmvnrm_Eigen.dll")<br />
</p>