原文連結：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 被歸一化——這意味著如果我們將 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，建立漂亮的圖。

概率

2. def _int(sef, daa):
3. ""二進位制資料的概率。""
4. elfus = {s:0 for s in ['0', '1']}
6. def proes_data(sef, ata):
7. ""處理資料。""
8. for s in ['0', '1']:
9. unts[s] =cont(s)
12. ""處理概率。""
13. n0 = couts['0']
14. n1 = couts['1']
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)
38. def prob(sef, p0):
39. ""獲取資料的概率。"""
40. prdat,_ = ._procs_proabiti(p0)
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. ""二進位制資料的β先驗。""
4. elfa0 = alha0
5. sl.1 = alha1
6. slf.0r = ba(self.0, sef.a1)
8. def intrval(slf, pb):
9. ""包含`prob'的pdf區域的端點。
11. 例如：interval(0.95)
12. """
14. return sef.pvitervl(rob)
16. def mean(self):
17. """返回先驗平均數。"""
20. def pdf(self, p0):
21. """p0處的概率密度。"""
23. return self.p.pdf(p0)
25. def plot(self):
26. ""顯示平均值和95%置信區間的圖。""
28. fig, ax = tsuplots(1, 1)
29. x = np.rae(0., 1., 0.01)
31. # 獲得先驗平均數 p0
32. mean = mean()
34. # 獲得包含95%概率的低/高點
35. low_p0, high_p0 = l.nterval(.95)
36. xo =nang(_0, hih_p0, 0.01)
38. # 繪製pdf
39. ax.plot(x, self.pdf(x), 'r-')
41. # 填充95%的區域
42. between(x_prob 0, sel.pdf(x_pob)
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. ""後驗。
4. data: 一個數據樣本的列表
5. 先驗：β先驗類的一個例項
6. """
7. elflilihod = liliod(dta)
8. sefprir = rio
10. self._css_steror()
11. def _pces_posteror(slf):
12. ""使用傳遞的資料和先驗來處理後驗。""
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
20. el0rv= beta(selfa0 + sfn0,
21. sef.a1 + slf.1)
22. def interalself, prob):
23. ""含`prob`的pdf區域的端點。
26. 例如：interval(0.95)
27. """
30. def mean(sef):
31. """返回後驗平均數。"""
34. def pdf(sef, p0):
35. """p0處的概率密度。"""
38. def plot(slf):
39. ""顯示先驗、似然和後驗的圖。""
41. ## 先驗
42. # 得到先驗平均數 p0
43. pri_mean =eorman()
45. #得到包含95%概率的低/高分值
46. plo_p0,= interval(0.95)
47. prob = arange(low_p0, i_hgh_0, 0.01)
49. # 繪製pdf
50. plot(prior.pdf(x)
52. # 填充95%的區域
53. x.ll_between(pri_p )
55. # 平均值
56. astm([pri_mean])
59. ## 似然
60. # 繪製似然圖
61. li = [sel.likliood.pro]
64. # ##後驗
66. #獲得後驗平均數p0
67. ostmen = mean()
69. #得到包含95%概率的低/高點
70. ow_p0, _high_p0 = interval(0.95)
71. prob = np.rngest_low_p0po_highp0 0.01)
73. # 繪製pdf
74. plot(x, slf.pd(x)
76. # 填充95%的區域
77. fil_etween(pos_xob, 0,self.df(pt_pr )
78. # 平均值
79. ax2].t_mean])

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

例子

讓我們從資料和統一先驗的示例開始。

1. # 資料
2. data1 = [0,0,0,0,1,1,0,0,0,1]
4. # 先驗
5. pro1 = prir(1, 1)
7. # 後驗
8. ot1 = psteior(daa1, ior1)
9. plot()

這裡需要注意的事項：

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

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

1. # 先驗
2. prir(4, 6)
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])
8. # 先驗
9. prior(1,1)
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模型實現