# Numerical integration over an infinite interval in Rcpp

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

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
```
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:

`integrate_exp4(4, 0, 4)`

```
## result error
## 1.077900e+00 9.252237e-08
```

However, it produces `NA`

for infinite ones:

`integrate_exp4(4, -Inf, Inf)`

```
## result error
## NaN NaN
```

This is disappointing, since base R’s `integrate()`

handles this without problems:

```
exp4 <- function(x, mean) exp(-(x - mean)^4 / 2)
integrate(exp4, 0, 4, mean = 4)
```

`## 1.0779 with absolute error < 1.3e-07`

`integrate(exp4, -Inf, Inf, mean = 4)`

`## 2.155801 with absolute error < 7.9e-06`

In this particular case the problem can be easily solved in two different ways. First, the integral can be expressed in terms of the Gamma function:

\[

\int_{-\infty}^{\infty} \mathrm{d}x \exp\left(-\frac{(x-\mu)^4}{2}\right) =

2^{-\frac{3}{4}} \Gamma\left(\frac{1}{4}\right) \approx 2.155801

\]

Second, the integrand is almost zero almost everywhere:

It is therefore sufficient to integrate over a small region around `mean`

to get a reasonable approximation for the integral over the infinite range:

`integrate_exp4(4, 1, 7)`

```
## result error
## 2.155801e+00 9.926448e-13
```

However, the trick to approximate the integral over an infinite range with an integral over a (possibly large) finite range does not work for functions that approach zero more slowly. The help page for `integrate()`

has a nice example for this effect:

```
## a slowly-convergent integral
integrand <- function(x) {1/((x+1)*sqrt(x))}
integrate(integrand, lower = 0, upper = Inf)
```

`## 3.141593 with absolute error < 2.7e-05`

```
## don't do this if you really want the integral from 0 to Inf
integrate(integrand, lower = 0, upper = 10)
```

`## 2.529038 with absolute error < 3e-04`

`integrate(integrand, lower = 0, upper = 100000)`

`## 3.135268 with absolute error < 4.2e-07`

`integrate(integrand, lower = 0, upper = 1000000, stop.on.error = FALSE)`

`## failed with message 'the integral is probably divergent'`

How does `integrate()`

handle the infinite range and can we replicate this in Rcpp? The help page states:

If one or both limits are infinite, the infinite range is mapped onto a finite interval.

This is in fact done by a different function from R’s C-API: `Rdqagi()`

instead of `Rdqags()`

. In principle one could call `Rdqagi()`

via Rcpp, but this is not straightforward. Fortunately, there are at least two other solutions.

The GNU Scientific Library provides a function to integrate over the infinte interval \((-\infty, \infty)\), which can be used via the RcppGSL package:

```
// [[Rcpp::depends(RcppGSL)]]
#include
```
#include
double exp4 (double x, void * params) {
double mean = *(double *) params;
return exp(-pow(x-mean, 4) / 2);
}
// [[Rcpp::export]]
Rcpp::NumericVector normalize_exp4_gsl(double &mean) {
gsl_integration_workspace *w = gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
F.function = &exp4;
F.params = &mean;
gsl_integration_qagi(&F, 0, 1e-7, 1000, w, &result, &error);
gsl_integration_workspace_free (w);
return Rcpp::NumericVector::create(Rcpp::Named("result") = result,
Rcpp::Named("error") = error);
}

`normalize_exp4_gsl(4)`

```
## result error
## 2.155801e+00 3.718126e-08
```

Alternatively, one can apply the transformation used by GSL (and probably R) also in conjunction with RcppNumerical. To do so, one has to substitute \(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}

\]

Now one could write the code for the transformed function directly, but it is of course nicer to have a general solution, i.e. use a class template that can transform *any function* in the desired fashion

```
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include
```
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::plugins(cpp11)]]
template class trans_func: public T {
public:
using T::T;
double operator()(const double& t) const {
double x = (1-t)/t;
return (T::operator()(x) + T::operator()(-x))/pow(t, 2);
}
};
// [[Rcpp::export]]
Rcpp::NumericVector normalize_exp4(const double &mean) {
trans_func f(mean);
double err_est;
int err_code;
const double result = Numer::integrate(f, 0, 1, err_est, err_code);
return Rcpp::NumericVector::create(Rcpp::Named("result") = result,
Rcpp::Named("error") = err_est);
}

`normalize_exp4(4)`

```
## result error
## 2.155801e+00 1.439771e-06
```

Note that the `exp4`

class is identical to the one from the initial example. This means one can use the same class to calculate integrals over a finite range and after transformation over an infinite range.

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