二項分布和Beta分布
二項分布
在概率論和統計學中，二項分布是n個獨立的[是/非]試驗中成功的次數的離散概率分布，其中每次試驗的成功概率為p。舉兩個例子就很容易理解二項分布的含義了：
- 拋一次硬幣出現正面的概率是0.5(p)，拋10(n)次硬幣，出現k次正面的概率。
- 擲一次骰子出現六點的概率是1/6，投擲6次骰子出現k次六點的概率。
在上面的兩個例子中，每次拋硬幣或者擲骰子都和上次的結果無關，所以每次實驗都是獨立的。二項分布是一個離散分布，k的取值范圍為從0到n，只有n+1種可能的結果。
n = 10
k = np.arange(n+1)
pcoin = stats.binom.pmf(k, n, 0.5)
[ 0.00097656, 0.00976563, 0.04394531, 0.1171875 , 0.20507813, 0.24609375, 0.20507813, 0.1171875 , 0.04394531, 0.00976563, 0.00097656 ]
下面是投擲6次骰子，出現6點的概率分布。
n = 6
k = np.arange(n+1)
pdice = stats.binom.pmf(k, n, 1.0/6)
[ 3.34897977e-01, 4.01877572e-01, 2.00938786e-01, 5.35836763e-02, 8.03755144e-03, 6.43004115e-04, 2.14334705e-05]
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處。
下面繪制n=10,k=4和n=20,k=8的概率分布。可以看出峰值都在p=0.4處，但是n=20的山峰更陡峭。也就是說隨着實驗次數的增加，p取其它值的可能就越小，對p的估計就更有信心，因此山峰也就更陡峭了。
MCMC
假設我們的知識庫中沒有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雖然是一個二項分布，但是它已經不能取其它值了。
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算法通常有一個收斂過程，我們希望只保留收斂之后的取樣值。
mc = pymc.MCMC([pcoin])
mc.sample(30000, 5000)
[****************100%******************] 30000 of 30000 complete
通過MCMC對象trace()可以獲得某個不確定量的取樣值。下面的程序獲得pcoin的25000次取樣值，並用hist()顯示其分布情況。由結果可知pcoin的分布與前面介紹的Beta分布一致。
pcoin_trace = mc.trace("pcoin")[:]
hist(pcoin_trace, normed=True, bins=30);
plot(p, pbeta, "r", label="n=20", lw=2)
pcoin_trace.shape