在Python中正常化後驗分佈時遇到問題

Question

我正在編寫一篇文章中Dirichlet-Multinomial posterior的推導，我已閱讀過該文章，並且遇到了將分佈總結爲1的問題。下面是代碼的非簡化形式：

def pcn(X, n, N, c, alpha): 
    pnc = np.math.factorial(np.sum([n[i] for i in range(len(n))]))/ \ 
      np.product([np.math.factorial(n[i]) for i in range(len(n))])* \ 
      np.product([c[i]**n[i] for i in range(len(n))]) 
    pc = G(len(X)*alpha)/ \ 
     np.product([G(alpha) for i in range(len(n)) if i in X])* \ 
     np.product([(c[i])**(alpha - 1) for i in range(len(n)) if i in X]) 
    pn = np.math.factorial(N)/ \ 
     np.product([np.math.factorial(n[i]) for i in range(len(n)) if i in X])* \ 
     G(len(X)*alpha)/ \ 
     G(len(X)*alpha + N)* \ 
     np.product([G(alpha + n[i])/G(alpha) for i in range(len(n)) if i in X]) 
    return pnc 

，我在這裏簡單去除得到劃分出來的零件：

def pcns(X, n, N, c, alpha): 
    pnc = np.product([c[i]**n[i] for i in range(len(n))]) 

    pc = np.product([(c[i])**(alpha - 1) for i in range(len(n))])/ \ 
     np.product([G(alpha) for i in range(len(n))]) 

    pn = np.product([G(alpha + n[i])/G(alpha) for i in range(len(n))])/ \ 
     G(len(X)*alpha + N) 

    return pnc * pc/pn 

我設置的輸入變量和C的數組初始化輸入：

X = [0,1] 
n = [6.0, 3.0] 
N = sum(n) 
alpha = 20 
c = np.linspace(0, 1, 1000) 

，然後我遍歷C和在每個C和情節評估我的後路功能：

dist = [] 
for i in c: 
    dist.append(pcns(X, n, N, [i, 1-i], alpha)) 

plt.plot(c,dist) 

當我總結了dist是999或len(c) - 1我得到的值。會有人碰巧知道爲什麼它不是1？

Answer 1

tl; dr：你正在計算一個離散的近似到一個確定的積分，並且忘記把這些片段乘以dx〜= delta_x。

---

即使這是一個正確的歸一化概率分佈，爲什麼要dist總和爲1？由於我們沒有G，我無法準確再現您的結果，所以我們來考慮一下簡單的高斯。我們可以設置平均0.5，將STDDEV的東西少，所以從0到1的範圍應包含大部分的概率：

>>> import scipy.stats 
>>> c = np.linspace(0, 1, 1000) 
>>> p = scipy.stats.norm(0.5,0.02).pdf(c) 
>>> sum(p) 
999.00000000000045 

此後又有這麼999。但是，爲什麼應該這是1？應該是一個的數量是總的綜合概率。這裏我們只是將某些點的概率分佈函數的值和它們加在一起。

更易例如：我們從0到1知道的x^2的定積分爲1/3，但

>>> x = np.linspace(0, 1, 1000) 
>>> sum(x**2) 
333.50016683350009 

當真正我們想要的東西更象（粗矩形近似，高於所需一點答案，因爲我們包括從x = 1點貢獻，盒子實際上是出了積分區）：

>>> sum(x**2) * (x[1]-x[0]) 
0.3338340008343344 

或者在你以前的情況

>>> sum(p) * (c[1]-c[0]) 
1.0000000000000004 

請注意，在這裏我們可以寫* (c[1]-c[0])，因爲c是等間隔的，所以Int(p(x) dx) ~= sum(p[x] * delta_x)中的每個「dx」都是相同的。一般來說，我們希望更多的東西像sum(p[:-1] * np.diff(c))。