来自 R 包 cubature 的 adaptIntegrate 得到一个巨大的积分错误

Posted 2023-03-24

【中文标题】来自 R 包 cubature 的 adaptIntegrate 得到一个巨大的积分错误

【英文标题】：adaptIntegrate from R package cubature gets an integral wrong by a huge margin

【发布时间】：2021-09-21 12:13:32

【问题描述】：

我对各种 2D 求积工具进行了一些测试，cubature 中的 adaptIntegrate 一直是最精确的工具之一，当它工作时。

麻烦的是，在某些情况下它完全出错了，但不是小数点，真的超出了规模，我不明白为什么。

我试图在矩形（实际上甚至是正方形）二维域上集成的功能是合乎逻辑的，形式如下：

条件_1 & (条件_2 | 条件_3 | ...)

这对我尝试过的任何工具都没有问题（pracma::integral2、pracma::quad2d、pracma::simpson2d...）。adaptIntegrate 然而，在大多数情况下提供最佳结果的同时，偶尔会完全失败。

例子：

require(cubature)

# Intersection of a circle of radius 4 centred in (0,0) with a circle of radius 1 centred in (0,0)

adaptIntegrate(function(x) ( (x[1]^2+x[2]^2 <= 16) & ((x[1]-0)^2+(x[2]-0)^2 <= 1) ), c(-4, -4), c(4, 4), absError = 1.e-2 )

#$integral
#[1] 3.141522
#
#$error
#[1] 0.009982141
#
#$functionEvaluations
#[1] 24089
#
#$returnCode
#[1] 0

正确：积分应该是 pi*1^2，大约是 3.142。 现在将小圆圈的中心移动到 (1,0)。它仍然完全包含在较大的圆圈中，因此交叉点仍然具有相同的区域。

# Intersection of a circle of radius 4 centred in (0,0) with a circle of radius 1 centred in (1,0)

adaptIntegrate(function(x) ( (x[1]^2+x[2]^2 <= 16) & ((x[1]-1)^2+(x[2]-0)^2 <= 1) ), c(-4, -4), c(4, 4), absError = 1.e-2 )

#$integral
#[1] 0
#
#$error
#[1] 0
#
#$functionEvaluations
#[1] 51
#
#$returnCode
#[1] 0

我无法弄清楚为什么会发生这种情况。如果我对小圆圈的 x 使用 1.1 而不是 1，它会回到一个非常令人满意的估计值。

如果我做错了什么，或者这实际上是一个错误，有没有人有任何想法？

谢谢！

PS 2 个注释：

是的，我知道这些例子很简单，我可以简单地使用小圆圈的面积。我需要解决的真正问题不那么微不足道（例如，小圆圈可以部分重叠大圆圈的周边，并且会有很多小圆圈，不仅是一个，它们之间也会重叠）。如果我不能依靠它来处理一个简单的案例，我就不能依靠它来处理更复杂的案例，可以吗？ 是的，我看到了一些与adaptIntegrate 不一致的帖子，例如this one。 我在那里看到了这个答案：

Cucture C 库给出的结果与我在上述问题中报告的结果相同，而且这不太可能是错误；相反，在某些情况下，h 自适应 cubature 例程（R 包是其接口）不如 Cubature 的 p 自适应例程准确，后者将适当区域的采样点数加倍。

这有什么帮助？当工具实际上给出错误答案时，我说 “这不是错误” 对我来说似乎是不正确的。此外，在我的情况下，pcubature 一直不太好，并且给出了相同类型的错误，所以......

【参考方案1】：

我无法解释为什么adaptIntegrate无法返回合理的结果。你可以写信给cubature的维护者或者这个软件的原作者。

但是，如果您将积分域限制为围绕较小/最小圆的矩形，您可以获得合理的结果！

最好将函数定义为返回数字 价值观。这也更容易写成几个圈子。

fn10 <- function(xy) 
    x <- xy[1]; y <- xy[2]
    (x^2+y^2 <= 16) * ((x-1)^2 + y^2 <= 1)


fn40 <- function(xy) 
    x <- xy[1]; y <- xy[2]
    (x^2+y^2 <= 16) * ((x-4)^2 + y^2 <= 1)


library(cubature)

adaptIntegrate(fn10, lowerLimit = c(0, -1), upperLimit = c(2, 1))
## $integral
## [1] 3.141828

adaptIntegrate(fn40, lowerLimit = c(3, -1), upperLimit = c(5, 1))
## $integral
## [1] 1.490316

我产生了兴趣并尝试通过不同的集成例程为f40 区域找到更好的价值。

library(pracma)

fun40 <- function(x, y)
    (x^2+y^2 <= 16) * ((x-4)^2 + y^2 <= 1)

quad2d(fun40, 3, 5, -1, 1, n = 512)                ## 200 ms
## [1] 1.489193

integral2(fun40,  3, 5, -1, 1)                     ##   3 ms
## [1] 1.494751

simpson2d(fun40, 3, 5, -1, 1, nx = 512, ny = 512)  ##  16 ms
## [1] 1.487213

那么这个积分的真实值是多少？ 两个圆相交于x0 = 31/8，因此一维 在区间[3, 4] 上积分的（矢量化）函数是：

x0 <- 31/8.0
fun <- function(x) 
    ifelse(x <= x0, sqrt(1 - (x-4)^2), sqrt(16 - x^2))


2 * quadgk(fun, 3, 4)
## [1] 1.487332

我们将积分加倍，因为我们只在 x 轴上方积分。 这应该精确到大约 0.00025 的误差。我们看到 简单的 Simpson 2D 求积非常接近。

注意：我们也可以利用三角公式来计算面积，即两个圆的面积，因此

# circle at (0.0)               circle at (4,0)
th1 = 2*acos(x0/4);             th2 = 2*acos((4-x0)/1)
A1 = 0.5*(th1 - sin(th1))*16;   A2 = 0.5*(th2 - sin(th2))*1

A1 + A2
## [1] 1.487332

这个值与之前在一维积分的帮助下找到的值一致。

【讨论】：

谢谢 Hans W。我有不止一两个圆圈，所以我不能限制矩形。要整合的领域是一个相当复杂的领域，由许多小圆圈组成，部分重叠，我想要其中的并集（你如何将其表示为数字？）然后这个并集必须与大圆圈相交。 This 是原始问题。是的，纯几何是一种选择，但在这种情况下实现起来并不容易。

至于给作者写信，这不是我经常做的事情，在这种情况下，我看不到任何明显的联系方式，至少在 CRAN 的包裹官方网页上：cran.r-project.org/web/packages/cubature/index.html。鉴于其他帖子给出的答案指出了 adaptIntegrate 的问题，不确定它是否值得。如果我的代码没有做错任何事情，并且问题确实是 cubature 函数的错误，那么我会接受它并使用其他工具。

这一切都可以处理。你提供了一个简单的例子，你得到了一个简单的答案。这就是生活。如果您想要更复杂的答案，请打开一个新线程，其中包含您应用程序中的真实示例。顺便说一句：我认为如果您使用的是开源软件，您会以某种方式提交错误报告。这将有助于我们所有人的生活变得更好。

好的，谢谢 Hans W.，那么我想您的回答应该停在第一句话。我确实做了一个注释（n. 1），实际问题更复杂，为了简洁起见，我保持简单。事实上，我的主要问题是我的代码是否有任何问题，或者更确切地说是我使用的工具中的错误。现在我们确定是后者。承诺提交错误报告，好的，我接受了，并且我确实这样做了（例如对于包 proxyC）。如果没有我可以联系任何人的联系方式，那么......（？）

维护者：B. Narasimhan 。这就是它在描述文件中所说的，很明显。它还链接到 Github 页面，在那里您可以找到 Cubature C 库作者的电子邮件地址。