引言

貝葉斯方法是天生用來做推斷的方法钓账，然而它常隱藏在課本的數(shù)學(xué)分析的背后。
隨著近年來貝葉斯方法在機(jī)器學(xué)習(xí)競(jìng)賽中成功應(yīng)用，其重要性又引起了學(xué)習(xí)者的興趣疏魏，但是其難點(diǎn)在于貝葉斯數(shù)學(xué)和概率編程之間的銜接。《Probabilistic Programming and Bayesian Methods for Hackers》一書試圖彌補(bǔ)以上的遺憾晤愧。
有關(guān)概率編程和貝葉斯方法的實(shí)驗(yàn)大莫，我將以該在線書籍作為學(xué)習(xí)資料，希望從中能將理論知識(shí)和編程實(shí)踐想結(jié)合官份，以加深對(duì)理論的理解只厘，以提高解決具體問題和數(shù)學(xué)建模的能力。

我們期望使用貝葉斯推斷來進(jìn)行統(tǒng)計(jì)分析舅巷，利用計(jì)算機(jī)的計(jì)算能力通過概率編程來解決問題羔味。該方法去除了原本數(shù)學(xué)分析的干預(yù)。PyMC是一個(gè)馬爾可夫鏈蒙特卡羅采樣(Markov chain Monte Carlo Sampling)Python工具包钠右，它實(shí)現(xiàn)了Metropolis-Hastings算法赋元，包含了畫圖、擬合優(yōu)度(goodness-of-fit)和聚合診斷(convergence diagnostics)的方法。

貝葉斯思想

比如測(cè)試一個(gè)算法的正確性搁凸，我們通過一個(gè)個(gè)測(cè)試用例來驗(yàn)證算法的正確性媚值，每次測(cè)試都讓我們對(duì)算法的準(zhǔn)確性更加有信心，但是我們無法保證算法是沒有bug的护糖，除非測(cè)試了所有可能的問題杂腰，但這常常是不實(shí)際的。
貝葉斯推斷是同樣的一個(gè)道理椅文，通過成功的結(jié)果來更新信任度喂很，在沒有驗(yàn)證所有可能的結(jié)果之前，我們無法絕對(duì)肯定皆刺。

現(xiàn)在少辣，我們來對(duì)比一下貝葉斯和頻率派對(duì)于概率的解釋。
從貝葉斯的世界觀中羡蛾，概率是一個(gè)事件發(fā)生的可信度(believability)的度量漓帅，也就是說，我們對(duì)一件事情發(fā)生是多么有信心痴怨。
而頻率派的觀點(diǎn)忙干，他們?cè)O(shè)定概率是事件的長(zhǎng)期頻率。比如飛機(jī)發(fā)生事故的概率浪藻，頻率派解釋為飛機(jī)事故的長(zhǎng)期頻率捐迫。這在很多事件的概率問題上看上去確實(shí)有邏輯合理性，但當(dāng)一個(gè)事件不會(huì)經(jīng)常出現(xiàn)爱葵，無法統(tǒng)計(jì)長(zhǎng)期的發(fā)生頻率的話施戴，就很難解釋通了。比如總統(tǒng)選舉的概率萌丈，因?yàn)樵撨x舉只發(fā)生一次赞哗。
貝葉斯觀點(diǎn)認(rèn)為，概率是事件發(fā)生的可信度或者信心的測(cè)度辆雾。對(duì)于一個(gè)事件的發(fā)生肪笋，如果將這種可信度設(shè)為0，那么說明我們對(duì)該事件發(fā)生沒有信心度迂；相反藤乙，如果將這種可信度設(shè)為1，那么說明我們十分肯定該事件一定會(huì)發(fā)生英岭。0到1之間的可信度是對(duì)結(jié)果發(fā)生的權(quán)重的考慮湾盒。這樣的定義與飛機(jī)發(fā)生事故的概率的解釋是一致的：通過觀測(cè)了飛機(jī)事故的發(fā)生頻率，在沒有附加信息的情況下诅妹，飛機(jī)事故的可信度應(yīng)該等于該頻率。

我們?cè)O(shè)定可信度的測(cè)量是針對(duì)一個(gè)個(gè)體的，而不是對(duì)于自然而言的吭狡。這是因?yàn)榧庋辏煌膫€(gè)體具有不同的信息，所以其對(duì)于事件發(fā)生的信任程度也不同划煮。
將可信度認(rèn)為是一種概率的哲學(xué)觀點(diǎn)對(duì)人類是很自然的送丰，我們持續(xù)不斷地與世界萬物互動(dòng)接觸，總是看到一部分真相弛秋，卻收集了一些證據(jù)來構(gòu)成這種可信度器躏。
貝葉斯派通過觀察到更多的證據(jù)和現(xiàn)象來更正自己的信任度，尤其是當(dāng)證據(jù)與初始認(rèn)為的一致時(shí)蟹略，該證據(jù)就更不能被忽視了登失。我們將更新的可信度設(shè)為P(A|X)，解釋為在給定證據(jù)X時(shí)挖炬，事件A發(fā)生的概率揽浙。這被稱作后驗(yàn)概率。

1. 先驗(yàn)概率P(A)：硬幣有50%的機(jī)會(huì)向上意敛；后驗(yàn)概率P(A|X)：你看到了硬幣落到地上時(shí)正面朝上馅巷，這就是信息X，于是我們就將硬幣朝上的概率設(shè)為1.0草姻，硬幣朝下的概率設(shè)為0.0钓猬。
2. 先驗(yàn)概率P(A)：大型復(fù)雜的程序代碼可能bug；后驗(yàn)概率P(A|X)：如果代碼通過了所有的測(cè)試用例X撩独，盡管仍然有bug存在的可能性逗噩，但是該可能性已經(jīng)減小了。

貝葉斯推斷實(shí)戰(zhàn)

如果是頻率派和貝葉斯派編寫一個(gè)函數(shù)跌榔，輸入一個(gè)統(tǒng)計(jì)問題异雁，而返回的結(jié)果是不同的。頻率派的推斷方程會(huì)返回一個(gè)數(shù)來表示一個(gè)估計(jì)僧须，而貝葉斯派的方程將返回一個(gè)概率纲刀。
以之前代碼調(diào)試問題為例，當(dāng)詢問“代碼通過了所有測(cè)試担平，那么該代碼是不是沒有bug了”示绊，頻率派的方程會(huì)返回“是”；而貝葉斯派的方程將返回“是”和“否”的概率暂论。

令N為我們手上數(shù)據(jù)證據(jù)的數(shù)量面褐，當(dāng)我們收集了無限多個(gè)證據(jù)，即N→∞取胎，貝葉斯的結(jié)果和頻率派的結(jié)果是一致的展哭。所以湃窍，對(duì)于很大的N，統(tǒng)計(jì)推斷或多或少有一定客觀性匪傍；而對(duì)于較小的N您市，該推斷就不太穩(wěn)定了：頻率派估計(jì)就有較大的方差和較大的置信區(qū)間(confidence intervals)。這卻是貝葉斯分析擅長(zhǎng)的役衡，貝葉斯派引入了先驗(yàn)概率茵休，返回一個(gè)概率，這樣就保有了對(duì)于小數(shù)據(jù)集情況下統(tǒng)計(jì)推斷所反映出來的不確定性(uncertainty)手蝎。

有人會(huì)認(rèn)為榕莺，當(dāng)N很大的時(shí)候，人們會(huì)不太在意這兩種技術(shù)的區(qū)別棵介，因?yàn)榇藭r(shí)它們提供了類似的推斷钉鸯，甚至更傾向于計(jì)算簡(jiǎn)單的頻率派方法。
而有專家[Andrew Gelman (2005)]認(rèn)為鞍时，樣本數(shù)量永遠(yuǎn)不會(huì)很大亏拉。如果N很小，為了得到一個(gè)足夠精確的推論逆巍，你需要獲得更多的數(shù)據(jù)及塘。但是一旦N足夠大，你便開始細(xì)分?jǐn)?shù)據(jù)以學(xué)習(xí)更多（比如锐极，在一個(gè)公開的民意投票中笙僚，一旦你對(duì)于全國(guó)的情況有一個(gè)很好的估計(jì)，你便會(huì)劃分不同的人群來進(jìn)行更多的估計(jì)灵再，比如按照性別肋层，不同的年齡，不同的地域等）翎迁。N不會(huì)很充足栋猖，因?yàn)橐坏┠阏J(rèn)為有充足的數(shù)據(jù)，你便轉(zhuǎn)向更加細(xì)致的問題汪榔，這樣你便需要更多的數(shù)據(jù)蒲拉。

頻率派方法是錯(cuò)誤的嗎？

頻率派方法依然很有用痴腌，并且在很多領(lǐng)域都是先進(jìn)水平(state-of-the-art)雌团。像最小二乘線性回歸、LASSO回歸士聪、期望最大算法都很有用锦援、很快速。貝葉斯方法補(bǔ)充了這些技術(shù)的不足剥悟，或者用彈性建模(flexible modeling)來闡述潛在的系統(tǒng)灵寺。

關(guān)于大數(shù)據(jù)

反常的是曼库，大數(shù)據(jù)的預(yù)測(cè)分析問題常常使用相對(duì)簡(jiǎn)單的算法，因此替久，我們辯解大數(shù)據(jù)預(yù)測(cè)的困難點(diǎn)不在于算法的使用凉泄，而是存儲(chǔ)和使用大數(shù)據(jù)的難度躏尉。
更多的難以分析的問題包含中數(shù)據(jù)(medium data)以至于小數(shù)據(jù)(small data)蚯根。

貝葉斯框架

入門實(shí)驗(yàn)

拋硬幣實(shí)驗(yàn)

假設(shè)你不知道拋硬幣正面朝上的的概率，你定義該真實(shí)比率為p胀糜，但是并沒有關(guān)于p的任何先驗(yàn)知識(shí)颅拦。 于是，我們嘗試拋硬幣并記錄其觀測(cè)結(jié)果教藻。有趣的問題來了距帅，隨著我們觀察到越來越多的數(shù)據(jù)，我們的推斷是如何發(fā)生變化的括堤？

%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(11, 9)

import scipy.stats as stats


dist = stats.beta
#多次試驗(yàn)
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
#構(gòu)造500個(gè)符合伯努利分布的隨機(jī)抽樣
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)


#使用伯努利分布的共軛先驗(yàn)——Beta分布
for k, N in enumerate(n_trials):
    sx = plt.subplot(len(n_trials) / 2, 2, k + 1)
    plt.xlabel("$p$, probability of heads") \
        if k in [0, len(n_trials) - 1] else None
    plt.setp(sx.get_yticklabels(), visible=False)#不顯示y軸刻度
    heads = data[:N].sum()#統(tǒng)計(jì)1的個(gè)數(shù)
    y = dist.pdf(x, 1 + heads, 1 + N - heads)#構(gòu)造Beta分布
    plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads))
    plt.fill_between(x, 0, y, color="#FF99CC", alpha=0.4)#填充色彩到曲線下區(qū)域
    plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)#在0.5處畫虛線

    leg = plt.legend()
    leg.get_frame().set_alpha(0.4)
    plt.autoscale(tight=True)
    
plt.suptitle("Bayesian updating of posterior probabilities",
             y=1.02,
             fontsize=14)

plt.tight_layout()

二項(xiàng)分布與Beta分布

二項(xiàng)分布
在概率論和統(tǒng)計(jì)學(xué)中，二項(xiàng)分布是n個(gè)獨(dú)立的[是/非]試驗(yàn)中成功的次數(shù)的離散概率分布去枷，其中每次試驗(yàn)的成功概率為p怖辆。
比如：

• 拋一次硬幣出現(xiàn)正面的概率是0.5是复，拋10次硬幣，出現(xiàn)k次正面的概率竖螃。
• 擲一次骰子出現(xiàn)六點(diǎn)的概率是1/6淑廊，擲 6次骰子，出現(xiàn)k次六點(diǎn)的概率特咆。
  每次拋硬幣或者擲骰子都和上次的結(jié)果無關(guān)季惩，所以每次實(shí)驗(yàn)都是獨(dú)立的。二項(xiàng)分布是一個(gè)離散分布腻格，k的取值范圍為從0到n画拾，只有n+1種可能的結(jié)果。

scipy.stats.binom為二項(xiàng)分布菜职，下面用它計(jì)算拋十次硬幣青抛，出現(xiàn)k次正面的概率分布。

n = 10
k = np.arange(n+1)
pcoin = stats.binom.pmf(k, n, 0.5)#k -> quantiles, n and p -> shape parameters
plt.stem(k, pcoin)#繪制火柴棍圖
plt.margins(0.1,0.05)#設(shè)置自動(dòng)縮放的邊緣

下面是投擲6次骰子酬核，出現(xiàn)6點(diǎn)的概率分布蜜另。

n = 6
k = np.arange(n+1)
pdice = stats.binom.pmf(k, n, 1.0/6)
plt.stem(k, pdice)
plt.margins(0.1)

Beta分布
貝塔分布是一個(gè)作為伯努利分布和二項(xiàng)式分布的共軛先驗(yàn)分布的密度函數(shù)，在機(jī)器學(xué)習(xí)和數(shù)理統(tǒng)計(jì)學(xué)中有重要應(yīng)用嫡意。貝塔分布中的參數(shù)可以理解為偽計(jì)數(shù)举瑰，伯努利分布的似然函數(shù)可以表示為，表示一次事件發(fā)生的概率鹅很，它為貝塔有相同的形式嘶居，因此可以用貝塔分布作為其先驗(yàn)分布。

對(duì)于硬幣或者骰子這樣的簡(jiǎn)單實(shí)驗(yàn)促煮，我們事先能很準(zhǔn)確地掌握系統(tǒng)成功的概率邮屁。然而通常情況下，系統(tǒng)成功的概率是未知的菠齿。為了測(cè)試系統(tǒng)的成功概率p佑吝，我們做n次試驗(yàn)，統(tǒng)計(jì)成功的次數(shù)k绳匀，于是很直觀地就可以計(jì)算出p = k/n芋忿。然而由于系統(tǒng)成功的概率是未知的，這個(gè)公式計(jì)算出的p只是系統(tǒng)成功概率的最佳估計(jì)疾棵。也就是說實(shí)際上p也可能為其它的值戈钢，只是為其它的值的概率較小。
例如有某種特殊的硬幣是尔，我們事先完全無法確定它出現(xiàn)正面的概率殉了。然后拋10次硬幣，出現(xiàn)5次正面拟枚，于是我們認(rèn)為硬幣出現(xiàn)正面的概率最可能是0.5薪铜。但是即使硬幣出現(xiàn)正面的概率為0.4众弓，也會(huì)出現(xiàn)拋10次出現(xiàn)5次正面的情況。因此我們并不能完全確定硬幣出現(xiàn)正面的概率就是0.5隔箍，所以p也是一個(gè)隨機(jī)變量谓娃，它符合Beta分布。

Beta分布是一個(gè)連續(xù)分布蜒滩，由于它描述概率$p$的分布滨达，因此其取值范圍為0到1。 Beta分布有alpha和beta兩個(gè)參數(shù)帮掉，其中alpha為成功次數(shù)加1弦悉，beta為失敗次數(shù)加1窒典。

連續(xù)分布用概率密度函數(shù)描述蟆炊，下面繪制實(shí)驗(yàn)10次，成功4次和5次時(shí)瀑志，系統(tǒng)成功概率p的分布情況涩搓。可以看出k=5時(shí)劈猪，曲線的峰值在p=0.5處昧甘，而k=4時(shí)，曲線的峰值在p=0.4處战得。
當(dāng)n=30充边，k=12時(shí)，曲線的峰值還在p=0.4處常侦，但是因?yàn)殡S著實(shí)驗(yàn)次數(shù)的增多浇冰，p取其他值的可能就會(huì)變小，對(duì)p的估計(jì)就更有信心聋亡，因此山峰就更陡峭了肘习。

n = 10
k = 5
p = np.linspace(0, 1, 100)
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plt.plot(p, pbeta, label="k=5", lw=2)
k = 4
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plt.plot(p, pbeta, label="k=4", lw=2)
n = 30
k = 12
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plt.plot(p, pbeta, label="k=3", lw=2)

plt.xlabel("$p$")
plt.legend(loc="best")

小結(jié)

這一小節(jié)，我們先對(duì)貝葉斯概念有個(gè)初步的了解坡倔，在下一小節(jié)中漂佩，我將進(jìn)一步地介紹不同的數(shù)學(xué)分布，并引入pymc庫的使用罪塔。

轉(zhuǎn)載請(qǐng)注明作者Jason Ding及其出處
Github博客主頁(http://jasonding1354.github.io/)
CSDN博客(http://blog.csdn.net/jasonding1354)
簡(jiǎn)書主頁(http://www.reibang.com/users/2bd9b48f6ea8/latest_articles)
百度搜索jasonding1354進(jìn)入我的博客主頁