# 1D numerical integration (again)

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>
using namespace Rcpp;

// [[Rcpp::plugins(cpp14)]]
// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double test_beta_integral(double alpha, double beta) {
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.