我使用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％的精度），例如蒙特卡洛积分或类似的东西？
我知道python中有
skmonaco库用于蒙特卡洛积分，但是积分的极限必须是数字，与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
... )
如您所见，内部积分的极限不取决于外部积分。
有人可以推荐库或有效地计算此积分的方法吗？
最佳答案
在scikit-monaco中处理非三次积分量的最简单方法是重新定义积分函数，以使其在积分区域之外返回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算法：这两种算法都可以“学习”积分的形状，因为它们可以在感兴趣的区域中更有效地分布点。
也就是说，您的集成区域基本上是一个旋转的矩形：
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算法（在scikit-monaco中为mcmiser），而不是统一采样。最后，根据您使用过的函数的名称来判断，似乎您正在进行某种形式的贝叶斯更新。如果是这样，您可以考虑将PyMC库用于Markov Chain Monte Carlo，它可能比通用的MC集成库更合适。
关于python - 有效地计算双积分，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/30424074/