拓端tecdat|PYTHON貝葉斯推斷計算：用BETA先驗分布推斷概率和可視化案例

原文鏈接：http://tecdat.cn/?p=24084 

原文出處：拓端數據部落公眾號

在這篇文章中，我將擴展從數據推斷概率的示例，考慮 0 和 1之間的所有（連續）值，而不是考慮一組離散的候選概率。這意味着我們的先驗（和后驗）現在是一個 probability density function (pdf) 而不是 probability mass function (pmf)。

我考慮了從數據序列推斷 p0，即零的概率：

我使用 p0的不同先驗來解決相同的問題，該先驗允許 0 和 1 之間的連續值而不是一組離散的候選值。

概率

我們推理問題的起點是 似然 ——觀察到的數據序列的概率，寫成就像我知道 p0 的值一樣：

 

為了清楚起見，我們可以插入 p0=0.6，並找到給定未知概率值的指定數據序列的概率：

概率的更一般形式，不是特定於所考慮的數據序列，是

其中 n0 是零的數量，n1 是考慮的任何數據序列 D 中的 1 的數量。

先驗 - Beta 分布

我們使用貝塔分布來表示我們的先驗假設/信息。其數學形式是

其中 α0 和 α1 是我們必須設置的超參數，反映我們關於 p0 值的假設/信息。但是，只需將p0 視為我們想要推斷的參數——忽略參數是概率。

請注意 后驗 pdf 也將是 Beta Distribution，因此值得努力適應 pdf。

先驗均值—— 大多數人想要一個數字或點估計來表示推理結果或先驗中包含的信息。然而，在貝葉斯推理方法中，先驗和后驗都是 pdf 或 pmf。獲得點估計的一種方法是取 相關參數相對於先驗或后驗的 平均值。例如，對於 Beta 先驗我們得到：

pdf 被歸一化—— 這意味着如果我們將 p0 從 0 積分到 1，我們會得到一個：

因為以下關系：

就我們而言，最重要的信息是 b 在 0 到 1 區間上進行歸一化，這對於像 p0 這樣的概率是必要的。

先驗假設和信息可以通過設置超參數來反映--超參數α0α0和α1α1影響pdf的形狀，使先驗信息的編碼更加靈活。

例如，使用α0=1, α1=1 不能反映p0 的優選值。這個pdf看起來像

另一個先驗可以指定 α0=5, α1=5，它更接近 p0=1/2 附近的值

最后，我們可以用α0≠α1得到非對稱的先驗，可以看到α0=2，α1=8。

 

關於設置超參數需要記住的一些事情：

• 如果 α0=α1，先驗將是對稱的，先驗平均值等於 Eprior[p0]=1/2。
• 如果 α0≠α1，則先驗將是不對稱的，先驗平均值不同於 1/21/2。
• 先驗的強度與α0+α1的總和有關。將α總和與數據中的n0+n1進行比較，將α的視為假數。這些和的相對大小控制着先驗和似然對后驗形狀的影響。這在下面的Python例子中會變得很清楚。

累積分布函數 (cdf)  beta累積分布函數 (cdf)讓我們計算 p0 小於或等於值 x 的概率。具體來說，cdf定義為：

該積分也被稱為不完全Beta ing積分，並表示為Ix(α0,α1)。

如果我們想知道p0在數值xl和xh之間的概率，我們可以用cdf來計算。
 

不完全 Beta 積分或 cdf 及其逆積分允許從先驗或后驗計算置信區間。使用這些工具，可以說p0 的值有 95% 的概率在某個范圍內——同樣，我們將使用 Python 代碼在下面繪制它。

貝塔分布是這個問題的共軛先驗--這意味着后驗將具有與先驗相同的數學形式（它也是一個貝塔分布），並更新了超參數。這種數學上的 "共鳴 "真的很好，讓我們不用MCMC就能做完整的貝葉斯推斷。

現在我們談談貝葉斯定理和這個問題的后驗pdf。

貝葉斯定理和后驗

我們的最終目標是后驗概率密度函數，結合似然和先驗，在考慮數據后對我們對p0的知識做一個更新的反映。后驗pdf的形式是（在這種情況下）。

換句話說，這是 給定數據序列 D 和先驗假設的 p0 的概率密度，由具有超參數 (α0,α1)的 Beta pdf 反映。

在這種情況下，貝葉斯定理采用以下形式：

其中后驗 P(p0|D,α0,α1) 為藍色，似然 P(D|p0)為黑色，先驗 P(p0| α0,α1)是紅色的。請注意，歸一化 邊際似然 （上述等式中的分母）現在是一個積分。

嘗試將貝葉斯定理視為關於 p0 從 假設 (α0,α1) 更新到 假設 + 數據 (D,α0,α1) 的信息：

試着把貝葉斯定理看作是關於p0的信息被從假設（α0,α1）更新為假設+數據（D,α0,α1）。

為了得到后驗pdf，我們必須在貝葉斯定理的分母上做積分。在這種情況下，利用貝塔分布的特性，就可以進行計算。該積分如下。

最后一行的積分定義了一個貝塔函數，在關於先驗的一節中討論過，並且有一個已知的結果。
 

這意味着分母，也叫邊際似然，等於。 

同樣，我們得到這個結果是因為Beta分布是我們所考慮的伯努利過程概率的共軛先驗。請注意，來自先驗的超參數已經被計數數據所更新。

這與人們預期的完全一樣，不需要做所有的數學計算。在任何情況下，在用Python實現這一點之前，有幾個注意事項。
 

• 后驗pdf在0到1的區間內被歸一化，就像我們推斷p0這樣的概率時需要的那樣。
• 后驗平均數，是對我們的推斷給出一個點估計的方法是

• 后驗的cdf和先驗的一樣，因為我們仍然有一個Beta分布--只不過，現在的參數是用數據更新的。在任何情況下，我們都可以用不完整的Beta積分和它的逆向找到置信區間，如上所述。

Python 中的推理代碼

首先，我們導入一些包，使用這些包來計算和繪制先驗、似然和后驗。此外，使用 matplotlib，在本例中為 ggplot，創建漂亮的圖。

概率

1.  
    
2.  
   def _int(sef, daa):
3.  
   ""二進制數據的概率。""
4.  
   elfus = {s:0 for s in ['0', '1']}
5.  
    
6.  
   def proes_data(sef, ata):
7.  
   ""處理數據。""
8.  
   for s in ['0', '1']:
9.  
   unts[s] =cont(s)
10.  
     
11.  
     
12.  
    ""處理概率。""
13.  
    n0 = couts['0']
14.  
    n1 = couts['1']
15.  
     
16.  
    if p0 != 0 and p0 != 1:
17.  
    # 典型情況
18.  
    lgpdaa = n0*log(p0) + \
19.  
    n1*log(1.-p0)
20.  
    rdta = exp(lgpata)
21.  
    elif p0 == 0 and n0 != 0:
22.  
    # 如果n0不是0，p0就不可能是0
23.  
    lpr_dta = -np.inf
24.  
    p_dta = exp(lgrata)
25.  
    elif p0 == 0 and n0 == 0:
26.  
    # 數據與p0=0一致
27.  
    lgpr_at = n1*np.log(1.-p0)
28.  
    p_ata =.exp(lo_dta)
29.  
    elif p0 == 1 and n1 != 0:
30.  
    # 如果n1不是0，p0就不可能是1
31.  
    lpdta = -inf
32.  
    praa = exp(lor_ta)
33.  
    elif p0 == 1 and n1 == 0:
34.  
    # 數據與p0=1一致
35.  
    log_dta = n0*log(p0)
36.  
     
37.  
     
38.  
    def prob(sef, p0):
39.  
    ""獲取數據的概率。"""
40.  
    prdat,_ = ._procs_proabiti(p0)
41.  
     
42.  
     
43.  
    def logprob(sef, p0):
44.  
    ""獲取數據的概率對數。"""
45.  
    _, lgpr_a = sel_presplies(p0)

先驗分布

我們的先驗類基本上是一個圍繞 scipy的包，有一個繪圖方法。注意 plot() 方法得到了 Beta 分布的平均值，並使用  interval() 方法得到了一個概率為 95% 的區域--這是使用不完整的 Beta 積分和上面討論的它的逆值完成的。
 

1.  
   def _int_(sef, alpa0=1, alph1=1):
2.  
   ""二進制數據的β先驗。""
3.  
    
4.  
   elfa0 = alha0
5.  
   sl.1 = alha1
6.  
   slf.0r = ba(self.0, sef.a1)
7.  
    
8.  
   def intrval(slf, pb):
9.  
   ""包含`prob'的pdf區域的端點。
10.  
     
11.  
    例如：interval(0.95)
12.  
    """
13.  
     
14.  
    return sef.pvitervl(rob)
15.  
     
16.  
    def mean(self):
17.  
    """返回先驗平均數。"""
18.  
     
19.  
     
20.  
    def pdf(self, p0):
21.  
    """p0處的概率密度。"""
22.  
     
23.  
    return self.p.pdf(p0)
24.  
     
25.  
    def plot(self):
26.  
    ""顯示平均值和95%置信區間的圖。""
27.  
     
28.  
    fig, ax = tsuplots(1, 1)
29.  
    x = np.rae(0., 1., 0.01)
30.  
     
31.  
    # 獲得先驗平均數 p0
32.  
    mean = mean()
33.  
     
34.  
    # 獲得包含95%概率的低/高點
35.  
    low_p0, high_p0 = l.nterval(.95)
36.  
    xo =nang(_0, hih_p0, 0.01)
37.  
     
38.  
    # 繪制pdf
39.  
    ax.plot(x, self.pdf(x), 'r-')
40.  
     
41.  
    # 填充95%的區域
42.  
    between(x_prob 0, sel.pdf(x_pob)
43.  
     
44.  
    # 平均值
45.  
    ax.stem([mean], [sf.df(mea)

讓我們使用新代碼繪制一些具有一序列參數的 Beta pdf。

統一先驗

帶點的垂直線顯示 pdf 均值的位置。陰影區域表示對於給定的 α0 和 α1 值，概率為 95% 的（對稱）區域。如果您想要平均值和置信區間的實際值，也可以獲取這些值：

print("先驗均值: {}".format(pri.mean()))

上面的其他先前示例也有效：

1.  
   prior(5, 5)
2.  
   plot()

和

1.  
   por(2, 8)
2.  
   iplot()

 

了解超參數所反映的先前假設的均值和不確定性很有用。

后驗

最后，我們為后驗構建類。正如您所料，我將數據和先驗作為參數，並從這些元素中提取后驗所需的參數。

1.  
   def __int_(sf datprior):
2.  
   ""后驗。
3.  
    
4.  
   data: 一個數據樣本的列表
5.  
   先驗：β先驗類的一個實例
6.  
   """
7.  
   elflilihod = liliod(dta)
8.  
   sefprir = rio
9.  
    
10.  
    self._css_steror()
11.  
    def _pces_posteror(slf):
12.  
    ""使用傳遞的數據和先驗來處理后驗。""
13.  
     
14.  
    # 從似然和先驗中提取n0, n1, a0, a1
15.  
    seln = slfliklihod.counts['0']
16.  
    sel.n1 = elf.lkelihodcnts['1']
17.  
    lfa0 = sf.prir.0
18.  
    self.a = sef.priora
19.  
     
20.  
    el0rv= beta(selfa0 + sfn0,
21.  
    sef.a1 + slf.1)
22.  
    def interalself, prob):
23.  
    ""含`prob`的pdf區域的端點。
24.  
     
25.  
     
26.  
    例如：interval(0.95)
27.  
    """
28.  
     
29.  
     
30.  
    def mean(sef):
31.  
    """返回后驗平均數。"""
32.  
     
33.  
     
34.  
    def pdf(sef, p0):
35.  
    """p0處的概率密度。"""
36.  
     
37.  
     
38.  
    def plot(slf):
39.  
    ""顯示先驗、似然和后驗的圖。""
40.  
     
41.  
    ## 先驗
42.  
    # 得到先驗平均數 p0
43.  
    pri_mean =eorman()
44.  
     
45.  
    #得到包含95%概率的低/高分值
46.  
    plo_p0,= interval(0.95)
47.  
    prob = arange(low_p0, i_hgh_0, 0.01)
48.  
     
49.  
    # 繪制pdf
50.  
    plot(prior.pdf(x)
51.  
     
52.  
    # 填充95%的區域
53.  
    x.ll_between(pri_p )
54.  
     
55.  
    # 平均值
56.  
    astm([pri_mean])
57.  
     
58.  
     
59.  
    ## 似然
60.  
    # 繪制似然圖
61.  
    li = [sel.likliood.pro]
62.  
     
63.  
     
64.  
    # ##后驗
65.  
     
66.  
    #獲得后驗平均數p0
67.  
    ostmen = mean()
68.  
     
69.  
    #得到包含95%概率的低/高點
70.  
    ow_p0, _high_p0 = interval(0.95)
71.  
    prob = np.rngest_low_p0po_highp0 0.01)
72.  
     
73.  
    # 繪制pdf
74.  
    plot(x, slf.pd(x)
75.  
     
76.  
    # 填充95%的區域
77.  
    fil_etween(pos_xob, 0,self.df(pt_pr )
78.  
    # 平均值
79.  
    ax2].t_mean])
80.  
     

基本代碼就是這樣，讓我們​​做一些例子。

例子

讓我們從數據和統一先驗的示例開始。

1.  
   # 數據
2.  
   data1 = [0,0,0,0,1,1,0,0,0,1]
3.  
    
4.  
   # 先驗
5.  
   pro1 = prir(1, 1)
6.  
    
7.  
   # 后驗
8.  
   ot1 = psteior(daa1, ior1)
9.  
   plot()

這里需要注意的事項： 

• 先驗是統一的。這意味着概率和后驗具有相同的形狀。
• 95%的置信區間同時顯示在先驗和后驗中。

接下來，讓我們考慮具有不統一先驗的相同數據。數據集長度為 10，因此 n0+n1=10。讓我們用 α0+α1=10 ，設置先驗，但先驗在與似然不同的位置達到峰值（也許有專家說這應該是先驗設置）：

1.  
   # 先驗
2.  
   prir(4, 6)
3.  
    
4.  
   # 后驗
5.  
   ps2 = postior(da1, pio2)
6.  
   plot()

顯然數據和專家在這一點上存在分歧。然而，因為先驗的權重設置為 10 並且數據序列的長度為 10，所以后驗峰值位於先驗峰值和似然峰值的中間。嘗試使用這種效果來更好地理解先驗超參數、數據集長度和結果后驗之間的相互作用。

作為最后一個例子，我們考慮最后一個例子的兩個變體， 首先我們使用統一先驗：

1.  
   # 設定概率為0
2.  
   p0 = 0.2
3.  
   #設置rng種子為2
4.  
   n.rand.sed(2)
5.  
   # 產生數據
6.  
   dta2 =rando.i([0,1], 500, p=[p0, 1.-p0])
7.  
    
8.  
   # 先驗
9.  
   prior(1,1)
10.  
     
11.  
    # 后驗
12.  
    poteir(daa2, pio3)

請注意，概率和后驗峰值在同一個地方，正如我們所期望的那樣。但是，由於數據集較長（500 個值），峰值要強得多。

最后，我們在同一數據集上使用“錯誤先驗”。在這種情況下，我們將保持先驗強度為 10，即 α0+α1=10：

1.  
   # 先驗
2.  
   prior(6,4)
3.  
   # 后驗
4.  
   poseor(data, pior4)

請注意，盡管先驗在錯誤的位置達到峰值，但概率和后驗非常相似。這個例子表明，如果先驗沒有設置得太強，合理數量的數據應該產生不錯的推理。一般來說，最好讓 n0+n1>α0+α1 並考慮先驗和后驗的形狀。

---

最受歡迎的見解

1.matlab使用貝葉斯優化的深度學習

2.matlab貝葉斯隱馬爾可夫hmm模型實現

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實現貝葉斯線性回歸模型

7.R語言使用貝葉斯 層次模型進行空間數據分析

8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實現