One-dimensional numerical integration in R with values > 2^1024?

user_8926425 0
Viewed 144

I am trying to use R to evaluate a definite integral of a univariate function that contains an expression of the form exp(exp(x)), with the upper bound on x exceeding 100. When using the base integrate() function, I get the error "non-finite function value" in this case, as the values involved exceed the largest number that can be represented using R's double precision arithmetic (2^1024 or ~10^300).

The Brobdingnag package is very helpful when it comes to dealing with very large numbers, but integrate() coerces all values to doubles internally, so if I try to define the integrand as, say, exp(as.brob(exp(x))) (an expression that can be evaluated), I just get a different error ("evaluation of function gave a result of wrong length"). I also tried using the integrateR() function from the Rmpfr package, but on my set-up (package version 0.8-40, R version 4.1.0), even just attempting to run the example code given in the documentation (integrateR(exp, mpfr(0,200), 1, rel.tol=1e-25, verbose=TRUE)) would abort my R session.

Is there an alternative to R's integrate() that can deal with the very large numbers arising in the integrand?

Edit from 2021-09-15: Minimal reproducible example of an attempted solution added below:

if(!require("Brobdingnag") {install.packages("Brobdingnag")}
library(Brobdingnag)

f <- function(x) {
  term1 <- 0.8361913 * exp(0.1516063*x) * exp(0.1788979 * exp(8.577809*x))
  term2 <- 0.9512496 * exp(8.577809*x) + 0.04875068
  return(term1 * term2)
}

log_f <- function(x) {
  x <- as.brob(x)
  term1 <- 0.1788979 * (exp(8.577809*x) - 1) + 0.1516063*x 
  term2 <- log(0.9512496 * exp(8.577809*x) + 0.04875068)
  return(term1 + term2)
}

u <- 131.6
log_S <- log(u) + log_f(u)    # +exp(1127.1)
integrate(function(x) {as.double(exp(log_f(x) - log_S))}, lower = 0, upper = u)
# 0 with absolute error < 0
numerical-integration precision r
1 Answers

A strategy that will probably work:

Pick a scaling value S that will prevent overflow of the full integral. For example, if the function is increasing in x (and positive) and we integrate from 0 to u, the integral is bounded by u*max(f(x)) = u*f(u). (If f is non-monotonic, you just have to have some reasonable way of guessing the order of magnitude of its maximum.) Then scale your integrand by this value:

log_S <- log(u) + log_f(u)
i0 <- integrate(function(x) { exp(log_f(x) - log_S) }, from = 0, to u)
log_result <- log(i0) + u*log_S

Where log_f() is the log of the original function (e.g. exp(x) if f is exp(exp(x)).

(If this log_S is too big, you can shrink it; the idea is to pick a scaling function that prevents the full integral from overflowing but doesn't allow the largest values of the integrand to overflow.)

Either (1) all the values are of large magnitude, in which case they won't be lost by underflow, or (2) the values that are not of large magnitude (and are thus lost by underflow) do not make up an important component of the integral.

A clever analyst could certainly come up with a function that will break this scheme (e.g. a function with a Dirac delta function or other narrow, tall spike in it), but this should be OK for a large class of 'not too nasty' functions.

You could also do a brute-force Riemann sum with Brobdingnag or Rmpfr or matrixStats::logSumExp (i.e., pick a sensible dx and compute sum(f(x_i)*dx)), but if you wanted fancy stuff like adaptive choice of dx you'd have to re-implement it yourself ...

user_1902770
Related