什麼是貝葉斯理論？
貝葉斯理論就是允許我們根據抽樣分佈（或稱似然度）和先驗分佈推匯出後驗分佈。
什麼是抽樣分佈？
抽樣分佈指的是在給定引數θ的條件下，X發生的概率。記做p(X|θ).
舉例說明，我們記錄了1000次拋硬幣的資料。如果1表示硬幣正面。可以用python這麼表示：
輸入 [2]:
|
1 2 3 |
import numpy as np data_coin_flips = np.random.randint(2, size=1000) np.mean(data_coin_flips) |
輸出[2]:
0.46800000000000003
抽樣分佈允許我們指定資料的生成方式。對拋硬幣來說，我們可以假設總體資料服從以 p 為引數的伯努利分佈。這個引數 p 是丟擲1的概率（拋硬幣得到正面）。返回結果為1的概率p以及結果為0的概率（1-p）。
可以看出這種分佈完美的切合拋硬幣事件。如果硬幣是均勻的情況下，我們知道p值為0.5，因為拋得1（正面）或0（反面）的機會是平均的。在python中我們可以通過下面兩條語句建立這樣的分佈：
輸入[3]:
|
1 2 |
bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000) np.mean(bernoulli_flips) |
輸出[3]:
0.46500000000000002
既然已經定義了我們所認為的資料生成方式，就可以根據給定引數來計算資料的似然概率P(X|θ)。 因為我們已經選定了伯努利分佈，只需要一個引數p就可以了。
我們可以使用伯努利分佈的概率質量函式（PMF）來計算單次拋硬幣的期望概率。PMF根據一個資料值和給定的引數（本例中為p），返回相應的概率。對伯努利分佈來說，這個pmf函式非常簡單： 如果資料值為1， 那麼它的概率返回p， 如果資料值為0， 那麼返回概率為(1-p)。 我們可以定義如下函式：
In [4]:
|
1 2 3 4 5 6 7 |
def bern_pmf(x, p): if (x == 1): return p elif (x == 0): return 1 - p else: return "Value Not in Support of Distribution" |
現在我們使用這個函式，根據輸入輸引數來獲取每個資料點的概率。如果引數p的值為0.5，那麼返回值也總是0.5.
In [5]:
|
1 2 |
print(bern_pmf(1, .5)) print(bern_pmf(0, .5)) |
0.5
0.5
這是一種非常簡單的 PMF，不過其他分佈的 PMF 會複雜得多。所以呢，我們應該瞭解這個模組 Scipy，它將這些都作為內建函式，直接呼叫就可以了。像下面這樣描述PMF：
In [6]:
|
1 2 3 |
import scipy.stats as st print(st.bernoulli.pmf(1, .5)) print(st.bernoulli.pmf(0, .5)) |
0.5
0.5
這個模組很給力，但是我們實際上想知道的是這1000個樣本資料的總體概率。怎麼才能知道呢？這裡有個小技巧，假設我們的樣本資料是獨立同分布的。那麼我們就可以簡單的認為所有資料的概率就是每個個體的獨立概率的乘積：
p(x1,…,xn|β)=p(x1|β)∗…∗p(xn|β)。Python中可以這樣實現：
In [7]:
|
1 |
np.product(st.bernoulli.pmf(data_coin_flips, .5)) |
Out[7]:
9.3326361850321888e-302
得出這樣的資料能做什麼用呢？ 說實話，單個資料本身其實並沒有什麼卵用。我們現在需要做的是給我們的抽樣模型生成更多的分佈。現在，我們僅僅使用概率0.5來測試我們的模型，如果概率為0.8呢？ 或者0.2呢？我們資料的概率又會怎樣？想要驗證這些情況， 可以定義一個概率p 的網格。下面我將構造100個0到1之間的資料組成的網格（因為概率值都在0和1之間），之後呢，我會在這些值的基礎上來計算我們抽樣資料的概率。
In [8]:
|
1 2 3 4 5 6 7 8 |
import matplotlib.pyplot as plt import seaborn as sns %matplotlib inlines ns.set(style='ticks', palette='Set2')params = np.linspace(0, 1, 100) p_x = [np.product(st.bernoulli.pmf(data_coin_flips, p)) for p in params] plt.plot(params, p_x) sns.despine() |
現在可以看出一點端倪了。 從圖中很明顯看出，當概率p=0.5的時候，資料分佈的概率達到峰值，概率0.4到0.6之間也還算靠譜。Nice。假設我們的資料服從伯努利分佈的話，我們可以很清楚的看出資料服從的引數 p 的值了。
先驗分佈
貝葉斯理論認為我們不僅要考慮抽樣分佈還得考慮先驗分佈。什麼是先驗分佈呢？就是p(θ)，或者說為每個引數θ指定一個明確的概率。在我們的抽樣分佈中，我們給引數p定義了100個0到1之間的資料。現在我們需要對每個值都設定一個先驗概率。這個概率是得到新資料前我們假設的概率。大多數情況下， 我們假設硬幣是均勻的。看起來就像上圖那樣的分佈。我們可以這樣做：
In [9]:
|
1 2 3 4 5 |
fair_flips = bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000) p_fair = np.array([np.product(st.bernoulli.pmf(fair_flips, p)) for p in params]) p_fair = p_fair / np.sum(p_fair) plt.plot(params, p_fair) sns.despine() |
基本上，我們就像之前一樣建立1000次拋硬幣資料並生成抽樣分佈（不過我們將數值除以抽樣分佈的總和使得所有值加起來等於1）.這樣我們就得到了一個“均勻分佈的硬幣”的先驗概率。這也就意味著在得到任何實際資料之前，我們主觀上認定硬幣是均勻的。我們也可以從上圖中直接看出先驗概率，0.5時達到頂峰，並且絕大多數處於0.4和0.6區間之內。
我知道你在想什麼——你覺得這相當無聊對吧？抽樣分佈和先驗分佈看起來也沒差啊。所以呢， 我們將它們組合在一起，我們仍然假設先驗概率是均勻的，而將抽樣資料換成非均勻的硬幣：
In [27]:
|
1 2 3 4 5 6 7 8 |
unfair_flips = bernoulli_flips = np.random.binomial(n=1, p=.8, size=1000) p_unfair = np.array([np.product(st.bernoulli.pmf(unfair_flips, p)) for p in params]) fig, axes = plt.subplots(2, 1, sharex=True) axes[0].plot(params, p_unfair) axes[0].set_title("Sampling Distribution") axes[1].plot(params, p_fair) axes[1].set_title("Prior Distribution") sns.despine()plt.tight_layout() |
哎呀，事情變得有趣了呢。我們的資料證明硬幣不是均勻分佈的（因為我們根據非均勻的概率0.8來生成的資料）， 但是我們先驗概率由堅稱硬幣是均勻的。現在我們該怎麼辦呢？
Bayes Theorem (Posterior Distribution)
貝葉斯理論允許我們從抽樣和先驗概率來推匯出後驗概率。後驗概率 P(θ|x)，通俗一點地說，就是在給定資料情況下，我們要計算的該假設的概率。仔細想一想就會知道，我們一直想要的就是這個後驗概率。這個給定的資料可能來自於問卷調查，也可能來自於網際網路，在這些給定資料下，我們需要找出什麼樣的假設引數最有可能生成我們的樣本資料。怎樣得到後驗分佈呢？首先我們需要一些數學理論的說明（別擔心，很簡單）。
根據定義，我們可知（如果你信不過我的數學水平， 這個頁面可以幫你複習）：
- P(A|B)=P(A,B)/P(B) 。通俗一點說，給定資料B的前提下A發生的概率等於他們的聯合概率除以B的概率。
- P(B|A)=P(A,B)/P(A) 。與之類似的， 給定資料A的前提下B發生的概率等於他們的聯合概率除以A的概率。
可能你也注意到了， 這倆公式有著相同的因子，所以呢：
- P(A,B)=P(A|B)∗P(B)
- P(A,B)=P(B|A)∗P(A)
因此：
P(A|B)∗P(B)=P(B|A)∗P(A)
即：
P(A|B)=P(B|A)∗P(A)/P(B)
這裡就很明顯啦，將θ和X帶入公式，替換A、B，可知：
P(θ|X)=P(X|θ)∗P(θ)/P(X)
完美！現在們就可以帶入我們知道的術語啦：
後驗概率 = 似然度 * 先驗概率P(X)
不過問題是：P(X)是個啥？ 換句話說，我們資料的概率是多少？有點詭異……
好吧，我們回到數學，繼續使用B和A的公式：
我們知道P(B)=∑AP(A,B) (如果忘了，這個頁面可以幫你複習 )
從剛才的討論，我們知道：
P(A,B)=P(A|B)∗P(A)
因此：
P(B)=∑AP(A|B)∗P(A)
帶入 θ 和 X ：
P(X)=∑θP(θ|X)∗P(θ)
帶入術語：
P(X)=∑θ*似然度 * 先驗概率
簡直不能更贊！但是 ∑θ 是啥意思呢？ 這是指我們所有引數值的總和。在我們拋硬幣的例子中， 我們給引數 p 定義了一百個數值，所以我們需要對所有值計算它的似然度*先驗概率並求和。這個值也就是貝葉斯理論的分母。因此，我們得出最後的公式：
後驗概率 = 似然度 * 先驗概率 / ∑θ * 似然度 * 先驗概率
囉嗦了一大堆。上程式碼：
In [29]:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
def bern_post(n_params=100, n_sample=100, true_p=.8, prior_p=.5, n_prior=100): params = np.linspace(0, 1, n_params) sample = np.random.binomial(n=1, p=true_p, size=n_sample) likelihood = np.array([np.product(st.bernoulli.pmf(sample, p)) for p in params]) #likelihood = likelihood / np.sum(likelihood) prior_sample = np.random.binomial(n=1, p=prior_p, size=n_prior) prior = np.array([np.product(st.bernoulli.pmf(prior_sample, p)) for p in params]) prior = prior / np.sum(prior) posterior = [prior[i] * likelihood[i] for i in range(prior.shape[0])] posterior = posterior / np.sum(posterior) fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8)) axes[0].plot(params, likelihood) axes[0].set_title("Sampling Distribution") axes[1].plot(params, prior) axes[1].set_title("Prior Distribution") axes[2].plot(params, posterior) axes[2].set_title("Posterior Distribution") sns.despine() plt.tight_layout() return posterior |
In [25]:
example_post = bern_post()
你可能注意到了，我為先驗概率和似然度設定了100個觀察資料。增加了分佈的方差。資料越多，分佈資料表現的越集中。同樣，當你有更多的資料來估計似然度時，先驗分佈的作用就會變小了。
In [30]:
moredata_post = bern_post(n_sample=1000)
從上圖也可以看出資料量大小對結果的影響。因為我們使用更多的資料來估計似然度，從而我們的後驗分佈概率看起來與似然度十分接近。非常酷。
Conclusion
到這裡，你對貝葉斯理論就有了一個大概的瞭解了。如果你曾經懷疑一枚硬幣是否均勻，現在你就知道怎麼來分析問題了。又或者你想估計投票結果的概率？或者其他是或否的問題。我知道你在想什麼——你覺得這些都不夠酷。你想要做出一些預言或者跑一個非常牛叉的演算法…… 未來有可能。因為我確實想要寫一個系列來介紹貝葉斯線性迴歸。希望這篇文章是這個系列的第一篇 :)
一些附註
1. 你可能注意到貝葉斯理論的分母是個常數。所以如果你想要獲得後驗概率最大化，不需要計算出這個常數的實際值。因為這個緣故，你常常可以看到後驗概率正比例於（數學符號 ∝ ）似然度 * 先驗概率。
2. 頻率論統計關心的是似然度。或者可以說頻率論就是無資訊先驗的貝葉斯（類似於均勻分佈）。但也不要太討厭頻率論； 因為在實際應用過程中許多貝葉斯理論依賴於頻率論統計。
3. 現在你知道頻率論統計關注點在似然度，那麼也就明白了為什麼人們經常會誤解頻率論的置信區間。似然度是 P(X|θ) ——或者說在假設前提下得到我們資料的概率。有那麼一點迷惑人，因為我們給定的是資料，而不是我們的假設。大多數頻率論模型都是取似然函式分佈的最大值（也稱為最大似然估計）。主要是找出什麼樣的引數使得頻率達到極值。這裡有個關鍵點是：他們將資料視為隨機的，所以頻率論的置信區間意味著：如果你不斷地有新的資料（可能是一些問卷調查）加入樣本中，並對每個新樣本來計算置信區間，那麼其中 95% 的置信區間都會包含你正在估計的總體引數。