在Python中正确计算双重积分

Question

我正在尝试使用scipy计算一个确定的双重积分。这个积分的被积函数有点复杂，因为它包含了一些概率分布，用来给每个x和y的值加权（就像混合模型一样）。我写的代码计算出来的结果是一个负数，但它应该在[0,1]之间。此外，这个计算大约花了半个小时。

我有两个问题。

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的数量级，所以这不应该是积分方法的问题。

我已经为这个问题挣扎了一段时间，非常感谢任何帮助。谢谢。

编辑：修正了一个拼写错误

Answer 1

这个积分方法给出的错误信息其实就是一个数字，它告诉你收敛的情况怎么样。你有没有试着计算一下被积函数的具体值呢？

顺便问一下：你是在做概率密度函数的积分吗？如果是的话，你确定你的积分范围是正确的吗？

Answer 2

有几种方法可以提高计算的速度。

1. 你可以使用 epsabs 和 epsrel 这两个参数来调整 dblquad 的容忍度，这样可以加快积分的速度。当然，这样做的结果可能会不那么准确，但在调试的时候这样做是可以的。

2. 你可以通过重新排列代码来大幅减少 integrand 中函数的调用次数（注意，这段代码未经测试）。

   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
   

   不过要注意，Python 在调用函数时会有一些额外的开销，所以如果你只是用纯 Python 来写，这样的优化效果不会太明显（使用像 Cython 这样的工具可能会稍微有帮助）。

我不太明白为什么积分会变成负数。如果你能给我一个 p 的例子，我或许能帮你找出原因，这样我们就可以实际运行你的代码了。