正確評估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中的雙積分

回答

相關問題