// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
namespace rstub {
// [...]
}
class exp4: public Numer::Func {
private:
double mean;
public:
(double mean_) : mean(mean_) {}
exp4
double operator()(const double& x) const {
return exp(-pow(x - mean, 4) / 2);
}
};
// [[Rcpp::export]]
::NumericVector integrate_exp4(const double &mean, double lower, double upper) {
Rcpp(mean);
exp4 functiondouble err_est;
int err_code;
double result = rstub::integrate(function, lower, upper, err_est, err_code);
return Rcpp::NumericVector::create(Rcpp::Named("result") = result,
::Named("error") = err_est);
Rcpp}
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
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 funcdouble lower;
double upper;
public:
(T _func, double _lower, double _upper) :
transform_infinite(_func), lower(_lower), upper(_upper) {}
func
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) {
::stop("At least on limit must be infinite.");
Rcpp} 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) {
= 0.0;
err_est = 0;
err_code 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);
= -1.0;
sign }
<T> g(f, lower, upper);
transform_infinitereturn 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.