我正在编写 rcpp 代码，我想在包“invgamma”中使用函数 dinvgamma(rinvgamma)。以下是我的全部代码。我尝试将包“invgamma”放入环境中，然后将其中的函数作为 Rcpp::Function 调用。
#include <Rcpp.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <R_ext/Utils.h>
using namespace Rcpp;
// [[Rcpp::export]]
RcppExport SEXP updatesigama2_mu(SEXP sigma2_mu,
SEXP mu,
SEXP u0,
SEXP v0,
SEXP K,
SEXP SS,
SEXP acc,
SEXP sigma2_mu_list)
{
BEGIN_RCPP
Rcpp::Environment invgamma("package:invgamma");
Rcpp::Function dinvgamma = invgamma["dinvgamma"];
Rcpp::Function rinvgamma = invgamma["rinvgamma"];
double xacc = Rcpp::as<double>(acc);
Rcpp::NumericVector xsigma2_mu_list(sigma2_mu_list);
Rcpp::NumericVector xmu(mu);//vector mu
double xsigma2_mu = Rcpp::as<double>(sigma2_mu);
int xK = Rcpp::as<int>(K);
int xSS = Rcpp::as<int>(SS);// time for irrecation
double xu0 = Rcpp::as<double>(u0);
double xv0 = Rcpp::as<double>(v0);
Rcpp::RNGScope scope;
int c = 0; int d = 0;
c = xu0 + 0.5*xK + 1;
d = xv0 + 0.5*sum(xmu);
for (int ss = 0; ss<xSS; ss++){//iteration
Rcpp::NumericVector tmp = rinvgamma(1,0,1);//proposal distribution Normal(0,10)
Rcpp::NumericVector u = Rcpp::runif(1);
Rcpp::NumericVector a = dinvgamma(tmp[0], c, pow(d,-1),d, false ) * dinvgamma(xsigma2_mu,1,0,1,false)/
(dinvgamma(xsigma2_mu,c,pow(d,-1),d,false)*dinvgamma(tmp[0],1,0,1,false))
xsigma2_mu_list[1] = tmp[0];
xsigma2_mu_list[2] = a[0];
if ( u[0] <= a[0] ){
xsigma2_mu = tmp[0];
xacc += 1;
}
}
return Rcpp::List::create(Rcpp::Named("sigma2_mu") = xsigma2_mu,
Rcpp::Named("acc") = xacc,
Rcpp::Named("sigma2_mu_list") = xsigma2_mu_list);
END_RCPP
}
我将它用作以下表格，但它不起作用。它是否错过了一些东西？
Rcpp::NumericVector a = dinvgamma(tmp[0], c, pow(d,-1),d, false ) * dinvgamma(xsigma2_mu,1,0,1,false)/
(dinvgamma(xsigma2_mu,c,pow(d,-1),d,false)*dinvgamma(tmp[0],1,0,1,false))
加载某个包中定义的 R 函数没有原则问题。但是，必须附加该软件包才能使其环境可用。请参阅rfunc()示例中的函数。在逆 Gamma 的情况下，根据 Gamma 函数定义您自己的函数会更容易。请参阅sugar()示例中的函数。
例子：
#include <Rcpp.h>
// [[Rcpp::export]]
Rcpp::List rfunc() {
Rcpp::Environment invgamma("package:invgamma");
Rcpp::Function dinvgamma = invgamma["dinvgamma"];
Rcpp::Function rinvgamma = invgamma["rinvgamma"];
Rcpp::NumericVector tmp = rinvgamma(5, 1);
Rcpp::NumericVector a = dinvgamma(tmp, 1);
return Rcpp::List::create(Rcpp::Named("tmp") = tmp,
Rcpp::Named("a") = a);
}
Rcpp::NumericVector rinvgamma(R_xlen_t n,
double shape,
double rate = 1.0) {
return 1.0/Rcpp::rgamma(n, shape, rate);
}
Rcpp::NumericVector dinvgamma(Rcpp::NumericVector x,
double shape,
double rate = 1.0,
bool log = false) {
Rcpp::NumericVector log_f = Rcpp::dgamma(1.0/x, shape, rate, true) - 2 * Rcpp::log(x);
if (log)
return log_f;
return Rcpp::exp(log_f);
}
// [[Rcpp::export]]
Rcpp::List sugar() {
Rcpp::NumericVector tmp = rinvgamma(5, 1);
Rcpp::NumericVector a = dinvgamma(tmp, 1);
return Rcpp::List::create(Rcpp::Named("tmp") = tmp,
Rcpp::Named("a") = a);
}
/*** R
library(invgamma)
set.seed(42)
rfunc()
set.seed(42)
sugar()
microbenchmark::microbenchmark(rfunc(), sugar())
*/
输出：
> library(invgamma)
> set.seed(42)
> rfunc()
$tmp
[1] 0.5156511 5.5426504 1.8711424 41.7271256 2.3376817
$a
[1] 0.5408323347 0.0271775313 0.1673728698 0.0005607317 0.1193024224
> set.seed(42)
> sugar()
$tmp
[1] 0.5156511 5.5426504 1.8711424 41.7271256 2.3376817
$a
[1] 0.5408323347 0.0271775313 0.1673728698 0.0005607317 0.1193024224
> microbenchmark::microbenchmark(rfunc(), sugar())
Unit: microseconds
expr min lq mean median uq max neval
rfunc() 115.098 117.1595 130.80325 117.9270 119.429 1342.420 100
sugar() 7.333 8.3810 26.03649 9.2195 10.023 1657.404 100
使用 Rcpp 糖的性能提升令人印象深刻！
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句