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

Posted 2023-03-24

【中文标题】R中值> 2 ^ 1024的一维数值积分？

【英文标题】：One-dimensional numerical integration in R with values > 2^1024?

【发布时间】：2021-11-09 04:34:29

【问题描述】：

我正在尝试使用 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

【参考方案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 你必须自己重新实现它......

【讨论】：

感谢您的巧妙解决方法；看起来它可能会有所帮助。你介意澄清一下你的符号吗？ log_f(x) 是指原被积函数的对数吗？ （幸运的是，我感兴趣的函数可以表示为几个更简单函数的乘积，因此它的对数很容易通过分析获得 - 一般情况可能更困难。）另外，我是否正确理解这只会起作用如果被积函数在 [0, u] 上单调递增？

我现在已经实现了您的解决方案，虽然它在大多数情况下可能有效，但它并不能解决我的问题。部分原因是我的错：我的f(x) 包括exp(exp(A*x)) 而不仅仅是exp(exp(x))，对于某些A，甚至log_f(x) 超过2^1024。我通过使用Brobdingnag 并将结果强制为双精度来解决此问题（按照您的逻辑，下溢丢失的值对积分贡献不大），但遇到了一个更基本的问题：f(x) 增长如此之快这使得任何x != u 的log_f(x) - log_S 等于0，所以积分的对数就是u*log_S。