У меня есть эта функция в R, которую я пытаюсь скопировать в Rcpp, но она вызывает у меня некоторые проблемы. Я получаю ошибку при попытке использовать функцию интегрирования RcppNumeric. Это дает мне следующий код ошибки, когда я пытаюсь добавить функцию в R: No matching constructor for initialization of 'BetaFunc'.
Я немного не уверен, как использовать класс public Func в Rcpp, который необходим для того, чтобы использоватьRcppNumeric интегрировать функцию.
Вот функция Rcpp, с которой у меня возникли проблемы:
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppArmadillo.h>
#include <RcppNumerical.h>
#include <RcppEigen.h>
#include <RcppArmadilloExtensions/sample.h>
using namespace Numer;
using namespace arma;
using namespace Rcpp;
typedef Eigen::Map<Eigen::MatrixXd> MapMat;
typedef Eigen::Map<Eigen::VectorXd> MapVec;
// Define the function to integrate over
class BetaFunc: public Func {
private:
const double alpha;
const double beta;
const NumericVector vX;
const int i;
const NumericVector n;
const NumericVector indx;
public:
BetaFunc(const double alpha_, const double beta_, int i_,
const NumericVector vX_,
const NumericVector n_,
const NumericVector indx_) : alpha(alpha_), beta(beta_), i(i_), vX(vX_), n(n_), indx(indx_) {}
double operator()(const double& dZ) const {
double Vr = Rf_dbeta(dZ, vX[i] + alpha, n[i] - vX[i] + beta, 0);
int int_size = indx.length();
for (int j = 0; j < int_size; j++) {
int di = indx[j];
Vr *= Rf_pbeta(dZ, vX[di] + alpha, n[di] - vX[di] + beta, 1, 0);
}
return Vr;
}
};
// [[Rcpp::export]]
// Estimate the Bayesian posterior probability of each alternative
// being the best binomial bandit
NumericVector best_binomial_bandit(NumericVector vX, NumericVector n, double alpha = 1,
double beta = 1) {
int iK = vX.size();
NumericVector vAns(iK);
IntegerVector vSeq = seq(0, iK-1);
for (int i = 0; i < iK; i++) {
// Vector of zeros
LogicalVector vLogi(iK);
vLogi[i] = 1;
// Create integervector with all values except value i
IntegerVector indx = vSeq[vLogi < 1];
const double lower = 0, upper = 1;
BetaFunc f(alpha, beta, i, vX, n, indx);
double err_est;
int err_code;
vAns[i] = integrate(f, lower, upper, err_est, err_code);
}
return vAns;
}
А вот функция R, которая работает как следует:
bbb <-
function(x, n, alpha=1, beta=1) {
k <- length(x)
ans <- numeric(k)
for (i in (1:k)) {
indx <- (1:k)[-i]
f <- function(z) {
r <- dbeta(z, x[i] + alpha, n[i] - x[i] + beta)
for (j in indx) {
r <- r * pbeta(z, x[j] + alpha, n[j] - x[j] + beta)
}
return(r)
}
ans[i] = integrate(f, 0, 1)$value
}
return(ans)
}
ЭтоЭто лишь минимальный пример, который я использую для проверки функций:
library(Rcpp)
sourceCpp("./best_binomial_bandit.cpp")
set.seed(1)
x <- c(0.1, 0.2, 0.3, 0.4, 0.5)
n <- c(1, 2, 3, 4, 5)
bb_res <- bbb(x, n, 1, 1)
bb_cpp_res <- best_binomial_bandit(x, n, 0.2, 0.4)
library(microbenchmark)
microbenchmark(bbb(x, n, 1, 1),
best_binomial_bandit(x, n, 1, 1))
Помощь очень важна.