Boost 1.66 has tanh-sinh quadrature
https://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/double_exponential.html
@randommm @Bob_Carpenter Is this better than what got halfway merged in?
It works well. You can save this as a .cpp file and click on the Source icon in the top right of the RStudio editor pane to play around with it.
#include <Rcpp.h>
#include <boost/math/quadrature/tanh_sinh.hpp>
using namespace Rcpp;
// [[Rcpp::plugins(cpp14)]]
// [[Rcpp::depends(BH)]]
// [[Rcpp::export]]
double test_beta_integral(double alpha, double beta) {
boost::math::quadrature::tanh_sinh<double> integrator;
auto f = [&](double x, double xc) {
return x <= 0.5 ? std::pow(x, alpha - 1) * std::pow(1 - x, beta - 1) :
std::pow(x, alpha - 1) * std::pow(xc, beta - 1);
};
double termination = std::sqrt(std::numeric_limits<double>::epsilon());
double error;
double L1;
size_t levels;
double Q = integrator.integrate(f, 0.0, 1.0, termination, &error, &L1, &levels);
double condition_number = L1/std::abs(Q);
std::cout << "error = " << error << std::endl
<< "L1 = " << L1 << std::endl
<< "levels = " << levels << std::endl <<
"condition_number = " << condition_number << std::endl;
return Q;
}
It is accurate down to about alpha = 0.025
and beta = 0.025
. It also comes with sinh-sinh quadrature for integration over the whole real line and exp-sinh quadrature for integration to infinity on one side.
All else being equal, I’d much rather go with Boost-maintained numerics code.