相关疑难解决方法(0)

我需要评估一个需要长循环的函数(后验分布).显然我不想在R本身内这样做,因此我使用"inline"和"Rcpp"来实现C++.但是,我发现在每个循环使用R函数的情况下,cxx函数的运行速度与运行R代码一样慢(参见下面的代码和输出).特别是,我需要在每个循环中使用多元正态累积分布函数,因此我使用mvtnorm包中的pmvnorm().

如何在cxx函数中使用此R函数并加快速度？我想了解为什么会这样,所以我将来可以在cxxfunction中使用其他R函数.

谢谢.

test <- cxxfunction(
  signature(Num="integer",MU="numeric",Sigma="numeric"),
  body='
  RNGScope scope;

  Environment stats("package:mvtnorm");
  Function pmvnorm = stats["pmvnorm"];

  int num = Rcpp::as<int>(Num);
  NumericVector Ret(1);
  NumericMatrix sigma(Sigma);
  NumericVector mu(MU);
  NumericVector zeros(2);

for(int i = 0; i < num; i++)
{
  Ret = pmvnorm(Named("upper",zeros),Named("mean",MU),Named("sigma",sigma));
}
return Ret;
',plugin="Rcpp"
)

system.time(
test(10000,c(1,2),diag(2))
)
    user  system elapsed 
    5.64    0.00    5.75 

system.time(
for(i in 1:10000){
pmvnorm(upper=c(0,0),mean=c(1,2),sigma=diag(2))
}
)
   user  system elapsed 
   5.46    0.00    5.57 

我正在尝试在R建筑中rcpp实现多元正态分布的一阶导数,以实现此处和此处发布的多元正态分布.

这是一个快速R实现

mvnormDeriv = function(..., mu=rep(0,length(list(...))), sigma=diag(length(list(...)))) {
    if(sd(laply(list(...),length))!=0)
        stop("The vectors not same length.")
    fn = function(x) -1 * c((1/sqrt(det(2*pi*sigma))) * exp(-0.5*t(x-mu)%*%solve(sigma)%*%(x-mu))) * solve(sigma,(x-mu))
    out = t(apply(cbind(...),1,fn))
    colnames(out) = c('x', 'y')
    return(out[,1])
}

以及一些基准测试数据:

set.seed(123456789)
sigma = rWishart(1, 2, diag(2))
means = rnorm(2)
X     = rmvnorm(10000, means, sigma[,,1])
x1    = X[,1]
x2    = X[,2]
benchmark(mvnormDeriv(x1,x2,mu=means,sigma=sigma),
    order="relative", replications=5)[,1:4]

该公式可以在矩阵食谱(2012),公式346中找到.

我没rcpp能从这里修改多变量法线的实现.这是一些代码,我曾经尝试过

// [[Rcpp::export]]
arma::vec dmvnormDeriv_arma(arma::mat x,  SEXP mu_sexp, arma::mat sigma, bool log …