正确评估Python中的双积分
我想用scipy来计算一个确定的双积分。被积函数有点复杂,因为它包含一些概率分布来给出x和y的每个值(像混合模型)的可能性有多大。下面的代码评估为负数,但它应该被[0,1]绑定。此外,计算需要大约半小时。正确评估Python中的双积分
我有两个问题。
1)有没有更好的方法来计算这个积分?
2)这个负值来自哪里?对我来说,最大的问题是如何加快计算速度,因为我可以在我的代码中发现后来自行导致负面的错误。
from scipy import stats
from scipy.integrate import dblquad
import itertools
p= [list whose entries are each different stats.beta(a,b) distributions]
def integrand(x,y):
delta=x-y
marg=0
for distA,distB in itertools.permutations(p,2):
first=distA.pdf(x)
second=distB.pdf(y)
weight1=0
weight2=0
for distC in p:
if distC == distA:
continue
w1=distC.cdf(x)-distC.cdf(y)
if weight1 == 0:
weight1=w1
else:
weight1=weight1*w1
marg+=(first*weight1*second)
I=delta*marg
return I
expect=dblquad(integrand,0,1,lambda x: 0, lambda x: x)
这实质上是要求两点之间的最大距离的期望值在分布向量中。积分的极限是yε[0,x]和xε[0,1]。这给了我大约-49,估计的积分误差为10e-10,所以它不应该归因于积分方法。
我一直在与此战斗一段时间,并感谢任何帮助。谢谢。
编辑:纠正错字
通过积分法给出的错误仅仅是一个数字,告诉你收敛行为是多么好。你有没有试图计算被积函数的显式值?
顺便说一句:你整合PDF的?如果是:你确定你的整合限制?
@ user485185:是的,被积函数包含pdf。换言之,它是(X-Y)* P(X-Y)。 P(X-Y)是计算x-y的概率,如下所示:对于给定的一对分布,加权该分布的概率给出值x或y(取决于你正在查看的变量时刻;评估为P_i(x)),其概率为分布集合的最小值或最大值(否则,您将不计算使用该特定分布的最大距离;评估为CDF_i(x)-CDF_i Y))。这我整合了x = [0,1]和y = [0,x]。 – Jason 2010-10-27 17:37:39
有几种方法可以提高计算速度。
您可以使用
epsabs
和epsrel
参数dblquad
增加你的集成tolreance。当然,你的结果不太准确,但是对于调试来说这很好。-
可以大大地重新排序的代码一样(警告,未经测试的代码)
def integrand(x, y): marg = 0.0 cdf = dict((id(distC), distC.cdf(x) - distC.cdf(y)) for distC in p) for distA in p: weight = numpy.prod(cdf[id(distC)] for distC in p if distC is not distA) marg += weight * distA.pdf(x) * sum( distB.pdf(y) for distB in p if distB is not distA) return (x-y) * marg
减少功能评估的数量
integrand
但是请注意,Python有函数调用相当的开销,在这么写这个纯Python不会让你太过分(使用类似Cython这个问题可能会有所帮助)。
我不知道为什么积分变为负值。也许我可以告诉你,如果你会举一个p
的例子 - 这将使我们能够真正尝试你的代码。
你看过http://code.google.com/p/mpmath/和http://code.google.com/p/sympy/ – pyfunc 2010-10-27 16:56:11
@pyfunc:我之前看过他们。 Sympy似乎不喜欢我的双重积分。 MPMath我认为使用一种类似的方法来评估积分,因为它是scipy所做的,所以它目前需要相当长的一段时间,上面的p矢量只包含三个分布。 – Jason 2010-10-27 19:34:29
我在任何地方都看不到psi1和psi2的定义,除非psi2总是小于psi1,否则不保证重量distC.cdf(psi1)-distC.cdf(psi2)不是负值。我不明白算法,不应该有像随机变量向量的维数(大于2)那么多的积分。如果太乱了,我会转向蒙特卡罗整合。 – user333700 2010-10-29 03:14:01