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

Question

集成（）给出了非常错误的答案：

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

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

有任何想法吗？

这似乎默认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))

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

Answer 1

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

我认为这意味着“警告清空者”。我注意到在您的示例中，绝对误差大于积分值。鉴于您知道所有 x 的 f(x) > 0，至少它让您有机会发现某些地方出了问题。抓住机会取决于你。

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

给

1 with absolute error < 9.8e-07

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

另见这篇文章。

Answer 2

试试包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")，部分介绍。我的重点。

这个 Rcubature包公开了底层 C cubature 库的 hcubature 和 pcubature 例程，包括矢量化接口。

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

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

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

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

Answer 3

进一步扩展@r2evan 和@Limey 的评论：

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

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

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

例如，pcubature()(in?cubature::pcubature说

对于少数 (<=3) 维的平滑被积函数，该算法通常优于 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")

白色方块是成功的（积分=1），黑色方块是坏的（积分=0）。

Answer 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

我认识到这是启发式的缓解措施，假设在集成之前了解有关您的发行版的一些信息，并且不是完美的“通用”解决方案。只是提供另一个视角。