如何将R包中的函数添加到rcpp代码中

Question

我正在编写 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))

Answer 1

加载某个包中定义的 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 糖带来的性能提升非常可观！

如何将R包中的函数添加到rcpp代码中

how to add function in R package into rcpp code

r

rcpp