Numerics

Numerical integration over an infinite interval in Rcpp (part 2)

In a previous post I have shown that without intervention RcppNumerical does not handle integration over infinite ranges. In this post I want to generalize the method to integrals where only one of the limits is infinite. In addition, I want to make …

Numerical integration over an infinite interval in Rcpp

On Stack Overflow the question was asked how to numerically integrate a function over a infinite range in Rcpp, e.g. by using RcppNumerical. As an example, the integral \[ \int_{-\infty}^{\infty} \mathrm{d}x \exp\left(-\frac{(x-\mu)^4}{2}\right) \] was given. Using RcppNumerical is straight forward. One defines a class that extends Numer::Func for the function and an interface function that calls Numer::integrate on it: // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::depends(RcppNumerical)]] #include <RcppNumerical.h> class exp4: public Numer::Func { private: double mean; public: exp4(double mean_) : mean(mean_) {} double operator()(const double& x) const { return exp(-pow(x-mean, 4) / 2); } }; // [[Rcpp::export]] Rcpp::NumericVector integrate_exp4(const double &mean, const double &lower, const double &upper) { exp4 function(mean); double err_est; int err_code; const double result = Numer::integrate(function, lower, upper, err_est, err_code); return Rcpp::NumericVector::create(Rcpp::Named("result") = result, Rcpp::Named("error") = err_est); } This works fine for finite ranges: