1D numerical integration (again)


#1

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?


#2

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.


#3

All else being equal, I’d much rather go with Boost-maintained numerics code.