http://hyry.dip.jp/tech/slice/slice.html/42

本文通過實例介紹二項分布和Beta分布的含義，并使用pymc對拋硬幣進行模擬實驗，從而獲得Beta分布。

二項分布和Beta分布

In?[15]: %pylab inline import pylab as pl import numpy as np from scipy import stats Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.

二項分布

在概率論和統計學中，二項分布是n個獨立的[是/非]試驗中成功的次數的離散概率分布，其中每次試驗的成功概率為p。舉兩個例子就很容易理解二項分布的含義了：

• 拋一次硬幣出現正面的概率是0.5(p)，拋10(n)次硬幣，出現k次正面的概率。

• 擲一次骰子出現六點的概率是1/6，投擲6次骰子出現k次六點的概率。

在上面的兩個例子中，每次拋硬幣或者擲骰子都和上次的結果無關，所以每次實驗都是獨立的。二項分布是一個離散分布，k的取值范圍為從0到n，只有n+1種可能的結果。

scipy.stats.binom為二項分布，下面用它計算拋十次硬幣，出現k次正面的概率分布。

In?[16]: n = 10 k = np.arange(n+1) pcoin = stats.binom.pmf(k, n, 0.5) pcoin Out[16]: array([ 0.00097656, 0.00976563, 0.04394531, 0.1171875 , 0.20507813,0.24609375, 0.20507813, 0.1171875 , 0.04394531, 0.00976563,0.00097656]) In?[17]: pl.stem(k, pcoin, basefmt="k-") pl.margins(0.1)

下面是投擲6次骰子，出現6點的概率分布。

In?[18]: n = 6 k = np.arange(n+1) pdice = stats.binom.pmf(k, n, 1.0/6) pdice Out[18]: array([ 3.34897977e-01, 4.01877572e-01, 2.00938786e-01,5.35836763e-02, 8.03755144e-03, 6.43004115e-04,2.14334705e-05]) In?[19]: pl.stem(k, pdice, basefmt="k-") pl.margins(0.1)

Beta分布

對于硬幣或者骰子這樣的簡單實驗，我們事先能很準確地掌握系統成功的概率。然而通常情況下，系統成功的概率是未知的。為了測試系統的成功概率p，我們做n次試驗，統計成功的次數k，于是很直觀地就可以計算出p=k/n。然而由于系統成功的概率是未知的，這個公式計算出的p只是系統成功概率的最佳估計。也就是說實際上p也可能為其它的值，只是為其它的值的概率較小。

例如有某種特殊的硬幣，我們事先完全無法確定它出現正面的概率。然后拋10次硬幣，出現5次正面，于是我們認為硬幣出現正面的概率最可能是0.5。但是即使硬幣出現正面的概率為0.4，也會出現拋10次出現5次正面的情況。因此我們并不能完全確定硬幣出現正面的概率就是0.5，所以p也是一個隨機變量，它符合Beta分布。

Beta分布是一個連續分布，由于它描述概率p的分布，因此其取值范圍為0到1。 Beta分布有α和β兩個參數，其中α為成功次數加1，β為失敗次數加1。

連續分布用概率密度函數描述，下面繪制實驗10次，成功4次和5次時，系統成功概率p的分布情況。可以看出k=5時，曲線的峰值在p=0.5處，而k=4時，曲線的峰值在p=0.4處。

In?[20]: n = 10 k = 5 p = np.linspace(0, 1, 100) pbeta = stats.beta.pdf(p, k+1, n-k+1) plot(p, pbeta, label="k=5", lw=2)k = 4 pbeta = stats.beta.pdf(p, k+1, n-k+1) plot(p, pbeta, label="k=4", lw=2) xlabel("$p$") legend(loc="best");

下面繪制n=10,k=4和n=20,k=8的概率分布。可以看出峰值都在p=0.4處，但是n=20的山峰更陡峭。也就是說隨著實驗次數的增加，p取其它值的可能就越小，對p的估計就更有信心，因此山峰也就更陡峭了。

In?[30]: n = 10 k = 4 p = np.linspace(0, 1, 100) pbeta = stats.beta.pdf(p, k+1, n-k+1) plot(p, pbeta, label="n=10", lw=2)n = 20 k = 8 pbeta = stats.beta.pdf(p, k+1, n-k+1) plot(p, pbeta, label="n=20", lw=2) xlabel("$p$") legend(loc="best");

用pymc模擬

假設我們的知識庫中沒有Beta分布，如何通過模擬實驗找出p的概率分布呢？pymc是一個用于統計估計的庫，它可以通過?先驗概率和?觀測值?模擬出?后驗概率?的分布。下面先解釋一下這兩個概率：

• 先驗概率：在貝葉斯統計中，某一不確定量p的先驗概率分布是在考慮"觀測數據"前，能表達p不確定性的概率分布。

• 后驗概率：在考慮相關證據或數據后所得到的不確定量p的概率分布。

拿前面拋硬幣的實驗來說，如果在做實驗之前能確信硬幣出現正面的概率大概在0.5附近的話，那么它的先驗概率就是一個以0.5為中心的山峰波形。而如果是某種特殊的硬幣，我們對其出現正面的概率完全不了解，那么它的先驗概率就是一個從0到1的平均分布。為了估計這個特殊硬幣出現正面的概率，我們做了20次實驗，其中出現了8次正面。通過這個實驗，硬幣出現正面的可能性的后驗概率就如上圖中的綠色曲線所示。

pymc庫可以通過先驗概率和觀測值模擬出后驗概率的分布，這對于一些復雜的系統的估計是很有用的。下面我們看看如何用pymc來對這個特殊硬幣出現正面的可能性進行估計：

• 首先pcoin是這個特殊硬幣出現正面的概率，由于我們沒有任何先驗知識，因此它的先驗概率是一個從0到1的平均分布(Uniform)。

• 假設我們做了20次實驗，其中8次為正面。根據前面的介紹可知，出現正面的次數符合二項分布(Binomial)，并且這個二項分布的概率p為pcoin。這個通過value參數指定了實驗的結果。因此experiment雖然是一個二項分布，但是它已經不能取其它值了。

In?[32]: import pymc pcoin = pymc.Uniform("pcoin", 0, 1) experiment = pymc.Binomial("experiment", 20, pcoin, value=8)

接下來通過MCMC對象模擬pcoin的后驗概率。MCMC是Markov chain Monte Carlo(馬爾科夫蒙特卡洛)的縮寫，它是一種用馬爾可夫鏈從隨機分布取樣的算法。通過調用MCMC對象的sample()，可以對pcoin的后驗概率分布進行取樣。這里30000為取樣次數，5000表示不保存頭5000次取樣值。這時因為MCMC算法通常有一個收斂過程，我們希望只保留收斂之后的取樣值。

In?[33]: mc = pymc.MCMC([pcoin]) mc.sample(30000, 5000) [****************100%******************] 30000 of 30000 complete

通過MCMC對象trace()可以獲得某個不確定量的取樣值。下面的程序獲得pcoin的25000次取樣值，并用hist()顯示其分布情況。由結果可知pcoin的分布與前面介紹的Beta分布一致。

In?[31]: pcoin_trace = mc.trace("pcoin")[:] hist(pcoin_trace, normed=True, bins=30); plot(p, pbeta, "r", label="n=20", lw=2) Out[31]: [<matplotlib.lines.Line2D at 0x5182190>] In?[34]: pcoin_trace.shape Out[34]: (25000,)

總結

以上是生活随笔為你收集整理的二项分布和Beta分布的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。