Rcpp: Error with numerical integration. No matching constructor for initialization

Question

I have this function in R, that I am trying to replicate into Rcpp but it gives me some problems. I get an error when trying to use the RcppNumerical integrate function. It gives me the following error code, when I try to source the function into R: No matching constructor for initialization of 'BetaFunc'.

I am a bit unsure how to use the class public Func in Rcpp that is needed in order to use the RcppNumerical integrate function.

Here is the Rcpp function that I am having trouble with:

// [[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;

}

And here is the R function that is working as it should:

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)
}

This is just a minimal example that I am using to test the functions:

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))

Help is greatly appreciated.

Answer 1

As @coatless mentioned in the comments, the function signature is off. To be precise, the last argument of the constructor for BetaFunc is defined as a NumericalVector but is called with an IntegerVector. From my point of view the latter is more appropriate. Changing the last argument of BetaFunc() to IntegerVector makes the example to compile. However, it is actually slower than your R code. And the results are not the same, either.