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

**R on Ralf Stubner**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 <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.

**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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.