二項分布和Beta分布

原文為： 二項分布和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分布有$\alpha$和$\beta$兩個參數，其中$\alpha$為成功次數加1，$\beta$為失敗次數加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,)