R中值> 2 ^ 1024的一维数值积分？

Question

我正在尝试使用 R 来评估包含 exp(exp(x)) 形式的表达式的单变量函数的定积分，x 的上限超过 100。使用基本 integrate() 函数时，我在这种情况下得到错误“非有限函数值”，因为所涉及的值超过了使用 R 的双精度算术可以表示的最大数字（2^1024 或 ~10^300）。

Brobdingnag 包在处理非常大的数字时非常有用，但 integrate() 在内部强制所有值加倍，所以如果我尝试将被积函数定义为，例如，exp(as.brob(exp(x)))（一个表达式可以评估），我只是得到一个不同的错误（“函数评估给出了错误长度的结果”）。我还尝试使用 Rmpfr 包中的 integrateR() 函数，但在我的设置（包版本 0.8-40，R 版本 4.1.0）中，甚至只是尝试运行文档中给出的示例代码（@ 987654329@) 将中止我的 R 会话。

是否有替代 R 的 integrate() 可以处理被积函数中出现的非常大的数字？

从 2021 年 9 月 15 日开始编辑：下面添加的尝试解决方案的最小可重复示例：

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

Answer 1

一个可能会奏效的策略：

选择一个缩放值S，以防止完整积分溢出。例如，如果函数在x（和正数）中递增，并且我们从 0 积分到 u，则积分以 u*max(f(x)) = u*f(u) 为界。 （如果f 是非单调的，你只需要有一些合理的方法来猜测它的最大值的数量级。）然后用这个值缩放你的被积函数：

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

其中log_f() 是原始函数的日志（例如exp(x) 如果f 是exp(exp(x))。

（如果这个log_S太大，你可以缩小它；想法是选择一个缩放函数，防止全积分溢出，但不允许被积函数的最大值溢出。）

要么 (1) 所有值都很大，在这种情况下，它们不会因下溢而丢失，或者 (2) 不是 大的值（因此下溢丢失）不构成积分的重要组成部分。

一个聪明的分析师当然可以想出一个可以打破这个方案的函数（例如，一个带有狄拉克三角函数或其他狭窄、高尖峰的函数），但这对于一大类“不太讨厌的功能。

您也可以使用 Brobdingnag 或 Rmpfr 或 matrixStats::logSumExp 进行蛮力黎曼求和（即，选择合理的 dx 并计算 sum(f(x_i)*dx)），但如果您想要诸如自适应选择之类的花哨的东西dx 你必须自己重新实现它......