正确评估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，所以它不应该归因于积分方法。

我一直在与此战斗一段时间，并感谢任何帮助。谢谢。

编辑：纠正错字

来源

2010-10-27 Jason

你看过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

通过积分法给出的错误仅仅是一个数字，告诉你收敛行为是多么好。你有没有试图计算被积函数的显式值？

顺便说一句：你整合PDF的？如果是：你确定你的整合限制？

来源

2010-10-27 17:10:28 Philip

@ 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的例子 - 这将使我们能够真正尝试你的代码。

来源

2010-11-07 12:34:33

正确评估Python中的双积分

回答

相关问题