集成（）给出完全错误的数字答案

【问题标题】：integrate() gives totally wrong number集成（）给出完全错误的数字
【发布时间】：2020-06-01 12:05:38
【问题描述】：

integrate() 给出了可怕的错误答案：

integrate(function (x) dnorm(x, -5, 0.07), -Inf, Inf, subdivisions = 10000L)
# 2.127372e-23 with absolute error < 3.8e-23

返回值显然应该是1（正常分布积分为1），但是integrate()返回的数字小得离谱，报错错误，没有警告......

有什么想法吗？

这似乎默认的integrate() 有可怕的错误......我只是偶然发现了这个！是否有任何可靠的 R 包来计算数值积分？

编辑：我尝试了包pracma，我发现同样的问题！：

require(pracma)
integral(function (x) dnorm(x, -5, 0.07), -Inf, Inf)
# For infinite domains Gauss integration is applied!
# [1] 0

编辑：嗯...深入挖掘，似乎他很难找到数值大于 0 的函数的非常狭窄的域。当我将限制设置为某个（非常接近到 0, 1) 分位数，它开始工作：

integral(function (x) dnorm(x, -5, 0.07), qnorm(1e-10, -5, 0.07), qnorm(1 - 1e-10, -5, 0.07))

但无论如何，这是一个非常可怕的问题......想知道是否有任何补救措施。

【问题讨论】：

这在某种意义上是混乱的，因为dnorm(x, -5.5, 0.07) 是“1”，但-6 又很小。如果您查看测试的x 值（可能在dnorm 之前添加message(paste(x, collapse="\n"))），它可能会提供一些关于它的估计的见解。

标签： r

【解决方案1】：

来自在线文档：“与所有数值积分例程一样，这些例程在有限的点集上评估函数。如果函数在几乎所有范围内近似恒定（特别是零），则结果和错误估计可能严重错误。”

我认为这是“买者自负”的意思。我注意到在您的示例中，绝对误差大于积分值。鉴于您知道所有 x 的 f(x) > 0，至少它让您有机会发现出了问题。抓住这个机会就看你自己了。

integrate( function(x) dnorm(x, -5, 0.07), -20, 10, subdivisions=1000L)

给予

1 with absolute error < 9.8e-07

在线文档中的警告告诉我，鉴于您对 buggy 的明显定义，您的问题的答案是“不，没有可靠的数值积分方法。不是 R 或任何其他语言”。不应盲目使用数值积分技术。用户需要检查他们的输入是否合理，输出是否合理。仅仅因为计算机给了你答案，你就相信答案是不好的。

另见this post。

【讨论】：

【解决方案2】：

进一步扩展 @r2evan 和 @Limey 的 cmets：

@Limey：对于这样的非常普遍的问题，根本没有办法保证通用的解决方案。

解决此类问题的一种方法是使用更多关于被积函数属性的知识（@r2evans 的回答）； answer referenced by @Limey 详细介绍了不同的问题。

您可能没有想到的一个“陷阱”是尝试一堆通用方法、调整设置等可能会误导您得出结论认为某些设置/方法一般优于您尝试的第一个未能得到正确答案。（有效的方法可能效果更好，因为它们通常更好，但在一个例子上尝试并不能证明这一点！）

例如，pcubature() 的描述（?cubature::pcubature 中说

该算法通常优于 h 自适应集成在几个 (

但是，请回想一下，pcubature() 恰好在您的示例中失败，这是一个平滑的低维情况 - 正是 pcubature() 应该执行更好的地方 - 这表明它可能是在这种情况下，hcubature() 有效，pcubature() 无效。

说明结果对参数的敏感程度（在本例中为下限/上限）：

library(emdbook)
cc <- curve3d(integrate( dnorm, mean=-5, sd=0.07,
        lower=x, upper=y, subdivisions=1000L)$value,
       xlim=c(-30,-10), ylim=c(0,30), n = c(61, 61),
       sys3d="image", col=c("black", "white"),
    xlab="lower", ylab="upper")

白色方块表示成功（integral=1），黑色方块表示错误（integral=0）。

【讨论】：

【解决方案3】：

试用包cubature。

library(cubature)

hcubature(function (x) dnorm(x, -5, 0.07), -Inf, Inf)
#$integral
#[1] 1
#
#$error
#[1] 9.963875e-06
#
#$functionEvaluations
#[1] 405
#
#$returnCode
#[1] 0

请注意，同一包中的函数 pcubature 也返回 0。

来自vignette("cubature")，简介部分。我的重点。

这个 R cubature 包暴露了 hcubature 和 pcubature 底层 C 立方库的例程，包括矢量化接口。

根据文档，建议使用 pcubature 仅用于平滑最多三个维度的被积函数。事实上，pcubature 例程的性能明显低于矢量化hcubature 在不适当的情况下。 因此，当您有疑问时，最好使用 hcubature.

由于在这种情况下，被积函数是正常密度、平滑的一维函数，因此有理由更喜欢pcubature。但它没有给出正确的结果。小插曲总结如下。

矢量化hcubature 似乎是一个很好的起点。

对于低维 (≤3) 的平滑被积函数，pcubature 可能值得一试。在生产包中使用之前进行试验。

【讨论】：

哇！！所以有可以处理这个的包，太棒了！但是你知道那么 hcubature 和 pcubature 有什么区别，为什么另一个返回 0 吗？只是想知道我可以依赖这个...
@TMS 我从文档中添加了一个引用。最后一点。

【解决方案4】：

有趣的解决方法：毫不奇怪，integrate 在采样值（在(-Inf,Inf) 上，不少于）更接近数据的“中心”时表现良好。您可以通过使用您的功能但暗示中心来减少这种情况：

无需调整：

t(sapply(-10:10, function(i) integrate(function (x) dnorm(x, i, 0.07), -Inf, Inf, subdivisions = 10000L)))
#       value        abs.error    subdivisions message call      
#  [1,] 0            0            1            "OK"    Expression
#  [2,] 1            4.611403e-05 10           "OK"    Expression
#  [3,] 6.619713e-19 1.212066e-18 2            "OK"    Expression
#  [4,] 7.344551e-71 0            2            "OK"    Expression
#  [5,] 3.389557e-06 6.086176e-06 3            "OK"    Expression
#  [6,] 2.127372e-23 3.849798e-23 2            "OK"    Expression
#  [7,] 1            3.483439e-05 8            "OK"    Expression
#  [8,] 1            6.338078e-07 11           "OK"    Expression
#  [9,] 1            3.408389e-06 7            "OK"    Expression
# [10,] 1            6.414833e-07 8            "OK"    Expression
# [11,] 1            7.578907e-06 3            "OK"    Expression
# [12,] 1            6.414833e-07 8            "OK"    Expression
# [13,] 1            3.408389e-06 7            "OK"    Expression
# [14,] 1            6.338078e-07 11           "OK"    Expression
# [15,] 1            3.483439e-05 8            "OK"    Expression
# [16,] 2.127372e-23 3.849798e-23 2            "OK"    Expression
# [17,] 3.389557e-06 6.086176e-06 3            "OK"    Expression
# [18,] 7.344551e-71 0            2            "OK"    Expression
# [19,] 6.619713e-19 1.212066e-18 2            "OK"    Expression
# [20,] 1            4.611403e-05 10           "OK"    Expression
# [21,] 0            0            1            "OK"    Expression

如果我们添加一个“居中”提示，我们会得到更一致的结果：

t(sapply(-10:10, function(i) integrate(function (x, offset) dnorm(x + offset, i, 0.07), -Inf, Inf, subdivisions = 10000L, offset = i)))
#       value abs.error    subdivisions message call      
#  [1,] 1     7.578907e-06 3            "OK"    Expression
#  [2,] 1     7.578907e-06 3            "OK"    Expression
#  [3,] 1     7.578907e-06 3            "OK"    Expression
#  [4,] 1     7.578907e-06 3            "OK"    Expression
#  [5,] 1     7.578907e-06 3            "OK"    Expression
#  [6,] 1     7.578907e-06 3            "OK"    Expression
#  [7,] 1     7.578907e-06 3            "OK"    Expression
#  [8,] 1     7.578907e-06 3            "OK"    Expression
#  [9,] 1     7.578907e-06 3            "OK"    Expression
# [10,] 1     7.578907e-06 3            "OK"    Expression
# [11,] 1     7.578907e-06 3            "OK"    Expression
# [12,] 1     7.578907e-06 3            "OK"    Expression
# [13,] 1     7.578907e-06 3            "OK"    Expression
# [14,] 1     7.578907e-06 3            "OK"    Expression
# [15,] 1     7.578907e-06 3            "OK"    Expression
# [16,] 1     7.578907e-06 3            "OK"    Expression
# [17,] 1     7.578907e-06 3            "OK"    Expression
# [18,] 1     7.578907e-06 3            "OK"    Expression
# [19,] 1     7.578907e-06 3            "OK"    Expression
# [20,] 1     7.578907e-06 3            "OK"    Expression
# [21,] 1     7.578907e-06 3            "OK"    Expression

我认识到这是对启发式的缓解，假设在集成之前了解您的分布，并且不是一个完美的“通用”解决方案。只是提供另一个视角。

【讨论】：

哇，这个抵消的想法不错！谢谢！
谢谢。我认为这是一种明智的黑客攻击，可以弥补不完美的启发式方法。从一开始就拥有更具弹性的启发式方法会更好。