scipy.integrate 和 mpmath.quad 的积分结果不同

Question

我正在集成一个函数来实现本文中的功能，并且我观察到不同的scipy.integrate.quad 和 mpmath.quad 的输出如下：

import mpmath as mp
from scipy import integrate

print(mp.quad(lambda x: x**58 * (1-x)**167 * (x+0.07053941908706418)**86, [0,1]))
print(integrate.quad(lambda x: x**58 * (1-x)**167 * (x+0.07053941908706418)**86, 0, 1))

# Results
# 3.99873917020827e-90
# (1.0583370200366968e-89, 1.5341117982755304e-89)

基于外部考虑到这一点，我有理由相信 SciPy 结果是正确的（我对返回的 PDF 进行了检查，SciPy 结果返回的值是一个有效的 PDF，在支持上积分为 1）。

我是否错误地调用了任一函数，或者具体是 mp.quad ？

Answer 1

您正确地调用了这些函数，但是 mpmath.quad 是为任意精度编写的，因此当您将精度设置得更高时，它是一个更可靠的参考。另外，您可以尝试它的两种方法；我通常会相信结果彼此一致。

import numpy as np
from scipy import integrate
from scipy.integrate._tanhsinh import _tanhsinh  # another reference from SciPy
from mpmath import mp  # note: NOT `import mpmath as mp`
mp.dps = 150  # set to 150 decimal digits of precision

def f(x):
    return x**58 * (1 - x)**167 * (x + 0.07053941908706418)**86

F1 = float(mp.quad(f, (0, 1), method='tanh-sinh'))
F2 = float(mp.quad(f, (0, 1), method='gauss-legendre'))
F3 = _tanhsinh(f, 0, 1).integral
F4 = integrate.quad(f, 0, 1)[0]

1.0583376645182504e-89
1.0583376645182504e-89
1.0583376645182522e-89
1.0583370200366968e-89

所以 scipy.integrate.quad 与默认设置接近，但仅与其他前七位一致。请注意，它对其结果不是很有信心：

integrate.quad(f, 0, 1)
# (1.0583370200366968e-89, 1.5341117982755304e-89)

第二个数字是绝对误差估计，因此即使它相当准确，但它并没有声称如此。

您可以尝试使用 epsabs 和 epsrel 收紧其容差，但在这种情况下它们没有任何效果。