貝葉斯推斷簡介

你是一名程序員榕栏，今天你實(shí)現(xiàn)了一個(gè)非常復(fù)雜的算法扒磁。你不知道這段代碼有沒有 bug。根據(jù)你以往的經(jīng)驗(yàn)缸榛，寫出 bug 是在所難免的内颗，畢竟你不是 Jeff Dean。你先拿一些簡單的測(cè)試用例 run 了一下找前，通過了纸厉。你又試了一些復(fù)雜的測(cè)試用例，也通過了沃缘。接下來，你連續(xù)試了一大堆用例水慨，都通過了。于是你松了一口氣朝抖，開始覺得這段程序可能沒有 bug 了。

如果你這樣思考問題砌滞，那么恭喜你绊茧，你是一個(gè)貝葉斯主義者华畏。所謂貝葉斯推斷，說白了胧卤，就是在觀察到新的證據(jù)之后唯绍，更新你的信念。貝葉斯主義者往往不會(huì)對(duì)某件事確信無疑枝誊，但他可以對(duì)某件事很有信心况芒。拿上面的例子來說，我們無法保證代碼 100% 沒有 bug（除非窮盡所有可能的情況）叶撒。但當(dāng)代碼通過了大量的測(cè)試用例之后绝骚，我們對(duì)它更有信心了。

我們經(jīng)常用到“概率”這個(gè)詞，但對(duì)什么是“概率”，其實(shí)是有不同看法的。

在貝葉斯學(xué)派看來，一個(gè)事件的概率是人們根據(jù)已有的知識(shí)和經(jīng)驗(yàn)對(duì)該事件發(fā)生的可能性所給出的個(gè)人信念^[1]洒宝，這種信念用區(qū)間 $[0,1]$ 中的一個(gè)數(shù)來表示将宪，可能性大的對(duì)應(yīng)較大的數(shù)。
與之相對(duì)的，在頻率學(xué)派看來，一個(gè)事件在獨(dú)立重復(fù)實(shí)驗(yàn)中發(fā)生的頻率會(huì)趨于一個(gè)極限，這個(gè)極限就是該事件的概率。

乍一看，頻率學(xué)派的觀點(diǎn)好像更嚴(yán)謹(jǐn)。舉例來說仗谆，你想知道拋一枚硬幣國徽朝上的概率阔涉，你就不停地拋這枚硬幣暖侨，統(tǒng)計(jì)國徽朝上的次數(shù)的占比葫掉，隨著你拋的次數(shù)的增多，你會(huì)發(fā)現(xiàn)這個(gè)占比漸漸趨于穩(wěn)定，這個(gè)極限值就是國徽朝上的概率。問題是，有些情形下，你是無法做獨(dú)立重復(fù)實(shí)驗(yàn)的窍奋。比如老板問你目前負(fù)責(zé)的某個(gè)項(xiàng)目有多大的概率能拿到業(yè)務(wù)成果窖逗。這個(gè)項(xiàng)目只能做一次，要么成功，要么失敗，怎么計(jì)算成功的頻率呢？對(duì)此，頻率學(xué)派可能會(huì)告訴你，假設(shè)有很多個(gè)平行世界，這個(gè)項(xiàng)目在一些平行世界里成功了扯俱，而在另一些平行世界里失敗了为流。這個(gè)項(xiàng)目成功的概率就是在平行世界里成功的頻率的極限莲祸。但作為貝葉斯主義者，我們不需要訴諸虛無縹緲的平行世界∈勇基于以往的工作經(jīng)驗(yàn)和對(duì)這個(gè)項(xiàng)目現(xiàn)狀的分析，我們依然能夠給出我們的信念，“老板鞠评，我有 7 成把握，這個(gè)項(xiàng)目能成÷匦裕”

我們把對(duì)事件 $A$ 發(fā)生的可能性的信念記為 $P(A)$ ，稱作事件 $A$ 的先驗(yàn)概率（prior probability）巩剖。當(dāng)觀察到證據(jù) $E$ 之后，我們更新了對(duì)事件 $A$ 的信念语稠。更新后的信念記為 $P(A|E)$ 衣式，稱作事件 $A$ 的后驗(yàn)概率（posterior probability）^[2]。那么，在觀察到證據(jù) $E$ 之后，我們?cè)趺床拍芨孪闰?yàn)概率來得到后驗(yàn)概率呢碧绞？很簡單照激，我們使用貝葉斯定理^[3]：
$\begin{aligned} P(A|E) &= \frac{P(E|A)P(A)}{P(E)} \\ &\propto P(E|A)P(A) \end{aligned}$
$P(E|A)$ 是當(dāng)事件 $A$ 發(fā)生時(shí)六水，我們觀察到證據(jù) $E$ 的概率港准。它是關(guān)于證據(jù) $E$ 的函數(shù)毛萌，稱作證據(jù) $E$ 的似然函數(shù)（likelihood）。

例 1：2020 年 12 月 23 日吨拗，你看完勇士 VS 籃網(wǎng)的比賽哈恰，庫里三分球 10 投 2 中。你想預(yù)測(cè)一下庫里在本賽季的三分球命中率遍希。因?yàn)檫@是庫里本賽季的第一場(chǎng)比賽拥诡，直接用 $2/10=20\%$ 顯然是不合理的没讲。那么該怎么做呢潘靖？

查閱資料誓军，我發(fā)現(xiàn)庫里在最近的兩個(gè)賽季三分球一共是 859 投 366 中。用 $\Theta$ 表示庫里本賽季三分球的命中率。根據(jù)歷史數(shù)據(jù)鹰贵，我認(rèn)為 $\Theta$ 服從參數(shù)為 $\alpha=366$ 和 $\beta=859-366=493$ 的 Beta 分布。即我對(duì)庫里本賽季三分球命中率 $\Theta=\theta$ 的信念為
$P(\Theta=\theta)=\mathrm{Beta}(\theta;\alpha=366, \beta=493)$
Beta 分布的期望為
$\mathbb E(\Theta)=\frac{\alpha}{\alpha+\beta}$
所以，在開始的時(shí)候电爹，我覺得庫里本賽季的三分球命中率會(huì)在 $366/(366+493)\approx42.61\%$ 附近波動(dòng)立哑。

已知庫里的三分命中率為 $\theta$ 塞俱，那么他投 $n$ 次三分球命中的次數(shù) $X$ 應(yīng)該服從參數(shù)為 $n$ 和 $\theta$ 的二項(xiàng)分布，即似然為
$P(X=x|\Theta=\theta)=\mathrm{Binomial}(x;n,\theta)$
在看到本場(chǎng)比賽庫里 10 投 2 中的結(jié)果后，我就可以更新我的信念了柔纵。我們有
$\mathrm{Beta}(\theta;\alpha,\beta) = \frac{1}{\mathrm{B}(\alpha, \beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}$
其中
$\mathrm B(\alpha, \beta) = \int_0^1\theta^{\alpha-1}(1-\theta)^{\beta-1}\mathrm d\theta$
又有
$\mathrm{Binomial}(x;n,\theta) = \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x}$
因此
$\begin{aligned} P(\Theta=\theta|X=x) &= \frac{P(X=x|\Theta=\theta)P(\Theta=\theta)}{P(X=x)} \\ &= \frac{P(X=x|\Theta=\theta)P(\Theta=\theta)}{\int_0^1 P(X=x|\Theta=\theta)P(\Theta=\theta)\mathrm d\theta} \\ &= \frac{\mathrm{Binomial}(x;n,\theta)\mathrm{Beta}(\theta;\alpha,\beta)}{\int_0^1 \mathrm{Binomial}(x;n,\theta)\mathrm{Beta}(\theta;\alpha,\beta)\mathrm d\theta} \\ &=\frac{\frac{1}{\mathrm B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\cdot\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x}}{\int_0^1 \frac{1}{\mathrm B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\cdot\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x} \mathrm d\theta }\\ &=\frac{\theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}}{\int_0^1 \theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}\mathrm d\theta} \\ &=\frac{1}{\mathrm B(\alpha+x, \beta+n-x)}\theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}\\ &\equiv \mathrm{Beta}(\theta;\alpha+x,\beta+n -x) \end{aligned}$
于是揭蜒，我的信念更新為
$\begin{aligned} P(\Theta=\theta|X=2) &= \mathrm{Beta}(\theta;\alpha=366+2,\beta=493+10-2)\\ &= \mathrm{Beta}(\theta;\alpha=368,\beta=501) \end{aligned}$
也就是說隐轩，在觀察到庫里首場(chǎng)比賽三分球 10 投 2 中之后，我更新了我的信念，認(rèn)為他本賽季的三分球命中率將會(huì)在 $368/(368+501)\approx42.35\%$ 附近波動(dòng)。

看到這些計(jì)算過程竣况，有些讀者可能已經(jīng)崩潰了。然而這已經(jīng)算是最簡單的情形之一了悴灵。但是別著急，我們并不一定要把后驗(yàn)分布計(jì)算出來骂蓖！我們所關(guān)心的無非是后驗(yàn)分布的期望积瞒、方差、偏度登下、峰度之類的統(tǒng)計(jì)量茫孔，如果我們能夠設(shè)法從后驗(yàn)分布中采樣，就可以利用這些樣本來估計(jì)這些統(tǒng)計(jì)量被芳。那么缰贝，有什么方法可以從一個(gè)復(fù)雜的概率分布中采樣呢？非常幸運(yùn)的是畔濒，馬爾可夫鏈蒙特卡洛（Markov chain Monte Carlo, MCMC）方法可以做到這一點(diǎn)揩瞪。關(guān)于這個(gè)方法，后續(xù)我會(huì)專門介紹篓冲。這里我們只講一下實(shí)際操作。

在 Python 下宠哄，我們可以使用 PyMC3 輕松實(shí)現(xiàn) MCMC 采樣壹将。拿上面的例子來說，

import pymc3 as pm
import seaborn as sns

with pm.Model() as model:
    # theta 的先驗(yàn)分布為 Beta 分布
    theta = pm.Beta(
        name='theta', 
        alpha=366, 
        beta=493
    )
    # 似然函數(shù)為二項(xiàng)分布
    likelihood = pm.Binomial(
        name='likelihood', 
        p=theta,
        n=10,
        observed=[2]
    )
    # 利用 MCMC 方法從后驗(yàn)分布中采樣
    posterior = pm.sample()

# 利用后驗(yàn)分布的樣本估計(jì)后驗(yàn)分布的期望
print(posterior.theta.mean())
# 0.4234174176263356

# 樣本的直方圖
sns.distplot(posterior.theta)

采樣結(jié)果

例 2： 下圖展示了某商品在一段時(shí)間內(nèi)每天的銷量毛嫉。請(qǐng)問在這段時(shí)間內(nèi)銷量的平均水平是否發(fā)生過變化诽俯？

商品的銷量

用 $S_t$ 表示第 $t$ 天的銷量，這種計(jì)數(shù)類型的隨機(jī)變量適合用 Poisson 分布來刻畫，即似然為
$S_t\sim\mathrm{Poisson}(\lambda)$
假設(shè)在第 $\tau$ 天暴区，銷量的水平發(fā)生了變化闯团，則
$\lambda=\begin{cases} \lambda_1,&t<\tau\\ \lambda_2,&t\geq\tau \end{cases}$
我們只需要給出 $\lambda_1$ 、 $\lambda_2$ 和 $\tau$ 的先驗(yàn)分布仙粱，然后利用 MCMC 方法從他們的后驗(yàn)分布中采樣即可房交。
關(guān)于 $\lambda_1$ 和 $\lambda_2$ ，我們只知道他們一定是正數(shù)伐割，且和這段時(shí)間的平均銷量應(yīng)該相差不多候味。因此，我們假設(shè)他們的先驗(yàn)分布為指數(shù)分布
$\begin{aligned} \lambda_1&\sim\mathrm{Exp}(\alpha)\\ \lambda_2&\sim\mathrm{Exp}(\alpha) \end{aligned}$
其中 $\frac{1}{\alpha}$ 為這段時(shí)間的平均銷量隔心。
關(guān)于 $\tau$ 白群，我們所知甚少，只能假設(shè)它的先驗(yàn)分布是離散均勻分布
$\tau\sim\mathrm{DiscreteUniform}(0, 89)$

同樣地硬霍，我們使用 PyMC3 從后驗(yàn)分布中采樣

with pm.Model() as model:
    alpha = 1.0 / sales.mean()
    lambda_1 = pm.Exponential('lambda_1', alpha)
    lambda_2 = pm.Exponential('lambda_2', alpha)
    
    tau = pm.DiscreteUniform('tau', lower=0, upper=len(sales)-1)
    
    t = np.arange(len(sales))
    lambda_ = pm.math.switch(t < tau, lambda_1, lambda_2)
    
    likelihood = pm.Poisson('likelihood', lambda_, observed=sales)
    
    posterior = pm.sample(20000, tune=8000, step=pm.Metropolis())

結(jié)果如下：

采樣結(jié)果

根據(jù)采樣的結(jié)果帜慢，

\lambda_1

、

\lambda_2

和

\tau

的均值分別是

\begin{aligned} \lambda_1 &\approx 17.96\\ \lambda_2 &\approx 24.00\\ \tau &\approx 50.40 \end{aligned}

銷量的變化

這幾乎就是正確答案了唯卖，因?yàn)樗^的“銷量”實(shí)際上是我用下面的代碼隨機(jī)生成的

from scipy.stats import poisson
import numpy as np

sales = []

for _ in range(52):
    sales.append(poisson.rvs(mu=18))
for _ in range(38):
    sales.append(poisson.rvs(mu=23))

sales = np.array(sales)

參考資料

Probabilistic Programming & Bayesian Methods for Hackers
Beta distribution - Wikipedia
Binomial distribution - Wikipedia
斯蒂芬-庫里NBA球員虎撲籃球

個(gè)人信念粱玲？那豈不是每個(gè)人都不同？的確如此耐床，因?yàn)槊總€(gè)人擁有的知識(shí)和經(jīng)驗(yàn)都不同密幔，對(duì)同一件事的判斷肯定不一樣。舉例來說撩轰，股(jiu)民(cai)和知道內(nèi)部消息的大佬對(duì)某只股票漲跌的可能性就有不同的信念胯甩，正所謂“時(shí)間的朋友不如領(lǐng)導(dǎo)的朋友”。 ?
注意到貝葉斯推斷可以迭代使用：后驗(yàn)概率可以當(dāng)作新的先驗(yàn)概率堪嫂，從而在觀察到更新的證據(jù)之后偎箫，再次更新我們的信念。 ?
需要說明的是皆串，這個(gè)公式并不是貝葉斯推斷所獨(dú)有的淹办。學(xué)習(xí)過概率論的朋友一眼就能看出來，這個(gè)就是課本上的“逆概公式”恶复。 ?

最后編輯于：2021.09.14 18:45:03

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末怜森，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子谤牡，更是在濱河造成了極大的恐慌副硅，老刑警劉巖，帶你破解...
沈念sama閱讀 218,755評(píng)論 6贊 507
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件翅萤，死亡現(xiàn)場(chǎng)離奇詭異恐疲，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,305評(píng)論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門培己，熙熙樓的掌柜王于貴愁眉苦臉地迎上來碳蛋，“玉大人，你說我怎么就攤上這事省咨∷嗟埽” “怎么了？”我有些...
開封第一講書人閱讀 165,138評(píng)論 0贊 355
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵茸炒，是天一觀的道長愕乎。經(jīng)常有香客問我，道長壁公，這世上最難降的妖魔是什么感论？我笑而不...
開封第一講書人閱讀 58,791評(píng)論 1贊 295
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮紊册，結(jié)果婚禮上比肄，老公的妹妹穿的比我還像新娘。我一直安慰自己囊陡，他們只是感情好芳绩，可當(dāng)我...
茶點(diǎn)故事閱讀 67,794評(píng)論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著撞反，像睡著了一般妥色。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上遏片，一...
開封第一講書人閱讀 51,631評(píng)論 1贊 305
城市分裂傳說
那天嘹害，我揣著相機(jī)與錄音，去河邊找鬼吮便。笑死笔呀，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的髓需。我是一名探鬼主播许师，決...
沈念sama閱讀 40,362評(píng)論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼僚匆！你這毒婦竟也來了微渠？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 39,264評(píng)論 0贊 276
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤咧擂，失蹤者是張志新（化名）和其女友劉穎逞盆，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體屋确，經(jīng)...
沈念sama閱讀 45,724評(píng)論 1贊 315
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 37,900評(píng)論 3贊 336
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了攻臀。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片焕数。...
茶點(diǎn)故事閱讀 40,040評(píng)論 1贊 350
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖刨啸，靈堂內(nèi)的尸體忽然破棺而出堡赔，到底是詐尸還是另有隱情，我是刑警寧澤设联，帶...
沈念sama閱讀 35,742評(píng)論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布善已，位于F島的核電站，受9級(jí)特大地震影響离例，放射性物質(zhì)發(fā)生泄漏换团。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,364評(píng)論 3贊 330
男人毒藥：我在死后第九天來索命
文/蒙蒙一宫蛆、第九天我趴在偏房一處隱蔽的房頂上張望艘包。院中可真熱鬧，春花似錦耀盗、人聲如沸想虎。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,944評(píng)論 0贊 22
一樁弒父案叛拷，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽舌厨。三九已至，卻和暖如春忿薇，著一層夾襖步出監(jiān)牢的瞬間裙椭，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 33,060評(píng)論 1贊 270
情欲美人皮
我被黑心中介騙來泰國打工煌恢，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留骇陈，地道東北人。一個(gè)月前我還...
沈念sama閱讀 48,247評(píng)論 3贊 371
代替公主和親
正文我出身青樓瑰抵，卻偏偏與公主長得像你雌，于是被迫代替她去往敵國和親。傳聞我的和親對(duì)象是個(gè)殘疾皇子二汛，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 44,979評(píng)論 2贊 355

貝葉斯推斷簡介

推薦閱讀更多精彩內(nèi)容