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

July 11, 2019
By

(This article was first published on R on Ralf Stubner, and kindly contributed to R-bloggers)

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, since I would like to write something like

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include 

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 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::infinity());
        bool lower_finite = (lower > -std::numeric_limits::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
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::QuadratureRule rule = Integrator::GaussKronrod41) {
    
    if (upper == lower) {
        err_est = 0.0;
        err_code = 0;
        return 0.0;
    }
    if (std::abs(upper) < std::numeric_limits::infinity() && 
        std::abs(lower) < std::numeric_limits::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 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.

To leave a comment for the author, please follow the link and comment on their blog: R on Ralf Stubner.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)