高效计算二重积分

Question

我使用 scipy 计算以下积分：

from scipy.stats import norm
def integrand(y, x):
    # print "y: %s  x: %s" % (y,x)
    return (du(y)*measurment_outcome_belief(x, 3)(y))*fv_belief(item.mean, item.var)(x)
    return dblquad(
        integrand, norm.ppf(0.001, item.mean, item.var), 
        norm.ppf(0.999, item.mean, item.var),
        lambda x: norm.ppf(0.001, x, 3), 
        lambda x: norm.ppf(0.999, x, 3))[0]

我有项目的信念状态（正态分布）和以项目的实际价值为条件的测量。（也是正态分布）使用这个积分我计算信息的价值（衡量这个项目有多有用）。

问题是，计算这个积分需要很多时间。是否有更有效的计算方法（我不需要 100% 的精度），比如 monte carlo 积分或类似的方法？

我知道 python 中有 skmonaco 库用于 monte carlo 积分，但积分的限制必须是数字，不像 scipy，内部积分限制，取决于外部（例如从上面

lambda x: norm.ppf(0.001, x, 3)

) 这里如何使用 skmonaco

计算二重积分

>>> from skmonaco import mcquad
>>> mcquad(lambda x_y: x_y[0]*x_y[1], # integrand
...     xl=[0.,0.],xu=[1.,1.], # lower and upper limits of integration
...     npoints=100000 # number of points
...     )

如您所见，此处的内部积分限制不依赖于外部积分。有人可以推荐图书馆或有效计算这个积分的方法吗？

Answer 1

蒙特卡洛积分是个好主意，实施起来也不难。但是，相对于其他方法，以统一方式对点进行采样收敛速度较慢，其他方法会迭代采样，并且在每一步都根据到该点的体积结果进行采样。

为此检查 recursive stratified sampling。

适应您的自定义区域（显然不是超立方体）可能不会很困难。

不幸的是，我不认为有 Python 库可以开箱即用。

Answer 2

在 scikit-monaco 中处理非三次积分体积的最简单方法是重新定义积分函数，使其 returns 0 在积分区域之外（参见 this文档部分):

def modified_integrand(xs):
    x, y = xs
    if norm.ppf(0.001, x, 3) < y < norm.ppf(0.999, x, 3):
        return integrand(xs) # the actual integrand
    else:
        return 0.0

正如 Ami Tavori 所说，如果积分区域的大部分为零，这将非常低效。为了解决这个问题，您可以使用 MISER 或 VEGAS 算法：这两种算法 'learn' 积分的形状，因为它们运行在感兴趣的区域更有效地分布点。

也就是说，您的积分区域基本上是一个旋转的矩形：

import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np

xs = np.linspace(-10, 10)
ys = np.linspace(-10, 10)

# Plot integration domain
# Red regions are in the domain of integration, blue
# regions are outside
plt.contourf(xs, ys, 
     [ [ 1 if norm.ppf(0.001, x, 3) < y < norm.ppf(0.999, x, 3) 
         else 0 for x in xs ] for y in ys ])

在这种情况下，您最好旋转积分坐标。例如，定义

r = x - y
R = (x + y)/2.0

你的被积函数是：

def rotated_integrand(rs):
    R, r = rs
    x = R + r/2.0
    y = R - r/2.0
    return integrand(np.array([x,y]))

沿 r 的积分限制现在是一个常数（-9.27..9.27，在您给出的示例中）。沿 R 的积分限制仍然是 (-inf, inf)，因此您需要逐渐增加沿 R 的积分区域，直到积分收敛。我肯定会建议为此使用 MISER 算法（mcmiser in scikit-monaco），而不是均匀采样。

最后，从您使用的函数名称来看，您似乎在进行某种形式的贝叶斯更新。如果是这种情况，您可以考虑使用 PyMC 马尔可夫链 Monte Carlo 库，这可能比通用 MC 集成库更合适。

高效计算二重积分

efficiently calculating double integral

python

scipy

montecarlo