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

Abstract

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 it more user friendly.

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>

namespace rstub {
// [...]
} 

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, double lower, double upper) {
    exp4 function(mean);
    double err_est;
    int err_code;
    double result = rstub::integrate(function, lower, upper, err_est, err_code);
    return Rcpp::NumericVector::create(Rcpp::Named("result") = result,
                                       Rcpp::Named("error") = err_est);
}

and have it correctly handle different input:

rbind(
  integrate_exp4(4,    0,   4),
  integrate_exp4(4, -Inf, Inf),
  integrate_exp4(4,    3, Inf),
  integrate_exp4(4, -Inf,   3)
)

        result        error
[1,] 1.0779003 9.252237e-08
[2,] 2.1558005 1.439771e-06
[3,] 1.9903282 4.250105e-11
[4,] 0.1654723 6.251315e-14

The only differences in the above code to the sample code from the previous post is the usage or rstub::integrate instead of Numer:integrate and the as yet unspecified rstub namespace. What is needed in that namespace? First, we will need a template class that does the necessary variable substitutions. In the case where both limits are infinite, we use as before \(x = (1-t)/t\) resulting in

\[ \int_{-\infty}^{\infty} \mathrm{d}x f(x) = \int_0^1 \mathrm{d}t \frac{f((1-t)/t) + f(-(1-t)/t)}{t^2} \]

If only one of the limits is infinite, we use the substitutions \(x = a + (1-t)/t\) and \(x = b - (1-t)/t\) resulting in

\[ \int_{a}^{\infty} \mathrm{d}x f(x) = \int_0^1 \mathrm{d}t \frac{f(a+(1-t)/t)}{t^2} \]

and

\[ \int_{-\infty}^{b} \mathrm{d}x f(x) = \int_0^1 \mathrm{d}t \frac{f(b-(1-t)/t)}{t^2} \]

For the C++ template class aggregation is used instead of inheritance, allowing to easily specify the limits:

template<class T> 
class transform_infinite: public Numer::Func {
private:
    T func;
    double lower;
    double upper;
public:
    transform_infinite(T _func, double _lower, double _upper) : 
                      func(_func), lower(_lower), upper(_upper) {}

    double operator() (const double& t) const {
        double x = (1 - t) / t;
        bool upper_finite = (upper <  std::numeric_limits<double>::infinity());
        bool lower_finite = (lower > -std::numeric_limits<double>::infinity());
        if (upper_finite && lower_finite) {
            Rcpp::stop("At least on limit must be infinite.");
        } else if (lower_finite) {
            return func(lower + x) / pow(t, 2);
        } else if (upper_finite) {
            return func(upper - x) / pow(t, 2);
        } else {
            return (func(x) + func(-x)) / pow(t, 2);
        }
    }
};

Finally we need a wrapper function for Numer::integrate which checks if both limits are finite or not:

using Numer::Integrator;
template<class T>
double integrate(const T& f, double lower, double upper,
                 double& err_est, int& err_code,
                 const int subdiv = 100, const double& eps_abs = 1e-8, const double& eps_rel = 1e-6,
                 const Integrator<double>::QuadratureRule rule = Integrator<double>::GaussKronrod41) {
    
    if (upper == lower) {
        err_est = 0.0;
        err_code = 0;
        return 0.0;
    }
    if (std::abs(upper) < std::numeric_limits<double>::infinity() && 
        std::abs(lower) < std::numeric_limits<double>::infinity()) {
        return Numer::integrate(f, lower, upper, err_est, err_code, subdiv, eps_abs, eps_rel, rule);
    } else {
        double sign = 1.0;
        if (upper < lower) {
            std::swap(upper, lower);
            sign = -1.0;
        }
        transform_infinite<T> g(f, lower, upper);
        return sign * Numer::integrate(g, 0.0, 1.0, err_est, err_code, subdiv, eps_abs, eps_rel, rule);
    }
}

If both limits are finite, Numer::integrate is used directly. Otherwise the function is transformed and Numer::integrate is used with adjusted range. In addition, it is first checked that the upper limit is actually larger than the lower limit. If this is not the case, one of the properties of integration is used to swap the limits and change the sign:

\[ \int_{a}^{b} \mathrm{d}x f(x) = -\int_{b}^{a} \mathrm{d}x f(x) \]

Thereby we get the correct result even when the limits have been exchanged:

rbind(
  integrate_exp4(4,    3, Inf),
  integrate_exp4(4,  Inf,   3)
)

        result        error
[1,]  1.990328 4.250105e-11
[2,] -1.990328 4.250105e-11

In the end we needed on template class and one template function, which could be put into a separate header file, to generalize Numer::integrate for integration over an infinite interval.