有效计算双积分
I使用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积分或类似的东西？
我知道python中有用于monte积分的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
... )正如你所看到的，内部积分的极限不依赖于外部积分。有人能推荐库或方法来有效地计算这个积分吗？
回答 2
Stack Overflow用户
发布于 2015-05-24 13:55:07
在scikit-摩纳哥中处理非立方集成卷的最简单方法是重新定义您的集成函数，以便在集成区域之外返回0(请参阅文档的这部分)：
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所说，如果一体化区域的大部分区域为零，这将是相当低效的。为了解决这个问题，您可以使用守财奴算法或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的集成区域，直到积分收敛为止。我肯定会推荐使用守财奴算法(mcmiser在科学，摩纳哥)，而不是均匀抽样。
最后，从您使用过的函数的名称来判断，看起来您正在进行某种形式的贝叶斯更新。如果是这样的话，您可以考虑使用马尔可夫链蒙特卡罗的PyMC库，这可能比通用MC集成库更合适。
Stack Overflow用户
发布于 2015-05-24 13:45:30
蒙特卡罗集成是一个好主意，也不是很难实现。然而，相对于其他方法，均匀采样点的收敛速度较慢，这些方法迭代采样，在每一步，根据到该点的体积结果进行采样。
请检查递归分层抽样是否有此问题。
很可能很难适应您的自定义区域(这显然不是一个超立方体)。
不幸的是，我并不认为有Python库可以开箱即用。
https://stackoverflow.com/questions/30424074
复制相似问题
腾讯云开发者
