1. 從“猜”到“算”為什么我們需要貝葉斯推斷如果你做過數(shù)據(jù)分析或者模型預(yù)測一定遇到過這樣的場景你有一個模型它基于一些初始假設(shè)比如“用戶點(diǎn)擊率是5%”來預(yù)測未來。然后你拿到了一批新的觀測數(shù)據(jù)比如實(shí)際投放廣告點(diǎn)擊率是7%。這時候你該怎么辦是固執(zhí)地堅(jiān)持最初的5%假設(shè)還是完全相信這7%的觀測結(jié)果傳統(tǒng)的頻率學(xué)派統(tǒng)計(jì)可能會告訴你用這7%的數(shù)據(jù)重新估計(jì)一個參數(shù)。但這似乎有點(diǎn)“喜新厭舊”完全拋棄了我們在拿到新數(shù)據(jù)之前的所有認(rèn)知那5%的假設(shè)并非憑空而來可能是基于歷史經(jīng)驗(yàn)或行業(yè)基準(zhǔn)。而貝葉斯推斷的核心思想恰恰是優(yōu)雅地解決了這個問題它讓我們能夠用新的觀測數(shù)據(jù)來“更新”我們舊有的信念從而得到一個融合了先驗(yàn)知識與新證據(jù)的、更靠譜的新信念。這個“先驗(yàn) - 觀測 - 后驗(yàn)”的過程就是貝葉斯推斷的骨架。公式P(θ|Data) ∝ P(θ) * P(Data|θ)看似簡單卻威力無窮。P(θ)是你的先驗(yàn)信念比如你認(rèn)為點(diǎn)擊率θ最可能在4%-6%之間P(Data|θ)是似然函數(shù)在給定某個點(diǎn)擊率θ的情況下觀測到當(dāng)前這批數(shù)據(jù)的可能性而P(θ|Data)就是后驗(yàn)分布——它綜合了你的先驗(yàn)和當(dāng)前數(shù)據(jù)告訴你現(xiàn)在應(yīng)該相信什么。那么問題來了。這個后驗(yàn)分布P(θ|Data)對于稍微復(fù)雜一點(diǎn)的模型其解析解往往異常復(fù)雜甚至不存在。比如你的參數(shù)θ不止一個或者你的模型結(jié)構(gòu)是非線性的。這時候我們就需要一種強(qiáng)大的計(jì)算方法來繞過數(shù)學(xué)上的高山直接窺探后驗(yàn)分布的模樣。這就是蒙特卡洛模擬登場的時候。它不跟你玩復(fù)雜的積分它用“暴力”但有效的方式——通過大量隨機(jī)抽樣來近似這個我們無法直接計(jì)算的后驗(yàn)分布。你可以把它想象成與其苦苦推導(dǎo)一個復(fù)雜函數(shù)的具體表達(dá)式不如讓計(jì)算機(jī)隨機(jī)生成成千上萬個可能的參數(shù)值然后看看這些值在考慮了數(shù)據(jù)和先驗(yàn)之后出現(xiàn)的頻率分布是怎樣的。這個頻率分布就是我們對后驗(yàn)分布的最佳近似。所以“蒙特卡洛模擬的貝葉斯推斷模型”這個標(biāo)題指向的正是現(xiàn)代數(shù)據(jù)分析中一個非常核心且實(shí)用的技術(shù)棧利用隨機(jī)抽樣的數(shù)值方法來實(shí)現(xiàn)對貝葉斯后驗(yàn)分布的推斷與探索。它特別適合那些模型假設(shè)相對靈活、參數(shù)關(guān)系復(fù)雜、或者你對不確定性量化有極高要求的場景比如金融風(fēng)險評估、流行病學(xué)模型預(yù)測、A/B測試的深入分析以及任何你覺得“傳統(tǒng)方法說不清楚”的復(fù)雜系統(tǒng)建模。2. 核心組件拆解先驗(yàn)、似然與蒙特卡洛采樣器要搭建這個模型我們必須先理解它的三個核心部件先驗(yàn)分布、似然函數(shù)和蒙特卡洛采樣算法。這三者環(huán)環(huán)相扣缺一不可。2.1 先驗(yàn)分布如何科學(xué)地“拍腦袋”先驗(yàn)分布P(θ)代表了我們在看到數(shù)據(jù)之前對模型參數(shù)θ的認(rèn)知。選擇一個合適的先驗(yàn)是貝葉斯分析的藝術(shù)也是容易引起爭議的地方。但它的選擇并非玄學(xué)而是有章可循的。無信息先驗(yàn)當(dāng)你對參數(shù)一無所知或者希望數(shù)據(jù)完全主導(dǎo)推斷時使用。例如對于一個介于0和1之間的概率參數(shù)如點(diǎn)擊率一個常見的無信息先驗(yàn)是Beta(1, 1)分布它等價于在[0,1]區(qū)間上的均勻分布。對于均值參數(shù)使用方差極大的正態(tài)分布如Normal(0, 1000)也可以近似表達(dá)“我什么都不知道”的態(tài)度。弱信息先驗(yàn)這是更推薦的做法。你基于領(lǐng)域知識或常識給參數(shù)一個合理的、但范圍較寬的限制。比如對于網(wǎng)頁點(diǎn)擊率你根據(jù)經(jīng)驗(yàn)知道它不太可能超過20%也不太可能低于0.1%。那么你可以設(shè)定一個Beta(2, 20)這樣的先驗(yàn)它的概率質(zhì)量大部分集中在0到0.2之間但又留有其他可能性的余地。這能防止模型在數(shù)據(jù)量極少時得出荒謬的結(jié)論。共軛先驗(yàn)這是一類數(shù)學(xué)上非常友好的先驗(yàn)。當(dāng)先驗(yàn)分布與似然函數(shù)屬于同一個分布族時后驗(yàn)分布也會屬于該族并且有解析解。例如二項(xiàng)分布的似然配合Beta先驗(yàn)后驗(yàn)依然是Beta分布。在蒙特卡洛方法普及前共軛先驗(yàn)是貝葉斯推斷的主要工具?，F(xiàn)在雖然我們不再依賴它來求解析解但選擇共軛先驗(yàn)作為起點(diǎn)依然能讓采樣更高效、更容易理解。注意先驗(yàn)的選擇會影響后驗(yàn)尤其是在數(shù)據(jù)量小的時候。一個基本原則是先驗(yàn)的信息強(qiáng)度應(yīng)該與你實(shí)際擁有的先驗(yàn)知識相匹配。不要用非常強(qiáng)的先驗(yàn)如方差極小的正態(tài)分布去扭曲數(shù)據(jù)本身傳達(dá)的信號除非你有極其充分的理由。2.2 似然函數(shù)連接參數(shù)與數(shù)據(jù)的橋梁似然函數(shù)P(Data|θ)衡量的是在給定參數(shù)值θ的條件下我們觀測到當(dāng)前這批數(shù)據(jù)的概率或概率密度。它反映了模型對數(shù)據(jù)的擬合程度。對于不同的數(shù)據(jù)類型和模型我們需要選擇不同的似然函數(shù)對于計(jì)數(shù)數(shù)據(jù)如點(diǎn)擊次數(shù)、失敗次數(shù)常用泊松分布或二項(xiàng)分布。對于連續(xù)測量數(shù)據(jù)如身高、溫度常用正態(tài)分布高斯分布。對于生存時間或等待時間數(shù)據(jù)常用指數(shù)分布或威布爾分布。構(gòu)建似然函數(shù)時一個關(guān)鍵假設(shè)是條件獨(dú)立性即各數(shù)據(jù)點(diǎn)在給定參數(shù)θ的條件下是相互獨(dú)立的。這使得聯(lián)合似然可以寫成單個數(shù)據(jù)點(diǎn)似然的乘積P(Data|θ) Π P(Data_i|θ)。這個假設(shè)簡化了計(jì)算也是很多模型的基礎(chǔ)。2.3 蒙特卡洛采樣器從后驗(yàn)分布中“撈”樣本這是將理論變?yōu)閷?shí)踐的關(guān)鍵一步。我們無法直接寫出后驗(yàn)分布P(θ|Data)的公式但我們可以設(shè)計(jì)一種算法讓它生成一系列隨機(jī)樣本{θ^(1), θ^(2), ..., θ^(N)}使得這些樣本的分布近似于后驗(yàn)分布。當(dāng)N足夠大時我們就可以用這些樣本的統(tǒng)計(jì)特性如均值、中位數(shù)、分位數(shù)來估計(jì)后驗(yàn)分布的特性。最經(jīng)典、應(yīng)用最廣的采樣器是馬爾可夫鏈蒙特卡洛。MCMC的核心思想是構(gòu)造一條馬爾可夫鏈?zhǔn)蛊淦椒€(wěn)分布恰好就是我們想要的后驗(yàn)分布。然后讓這條鏈運(yùn)行足夠長的時間“老化”階段之后產(chǎn)生的樣本就近似來自后驗(yàn)分布。Metropolis-Hastings算法是MCMC的基石。它的步驟非常直觀從一個初始參數(shù)值θ^(0)開始。對于每一次迭代t1, 2, ... a.提議根據(jù)一個提議分布q(θ* | θ^(t-1))比如一個以當(dāng)前值為中心的正態(tài)分布生成一個候選新值θ*。 b.計(jì)算接受率α min(1, [P(θ*)P(Data|θ*)] / [P(θ^(t-1))P(Data|θ^(t-1))] * [q(θ^(t-1)|θ*) / q(θ*|θ^(t-1))])。后一項(xiàng)是提議分布的修正項(xiàng)如果提議分布對稱如正態(tài)分布則此項(xiàng)為1。 c.決定以概率α接受θ*令θ^(t) θ*否則拒絕令θ^(t) θ^(t-1)。Gibbs抽樣是MH算法的一個特例適用于參數(shù)可以分成多個塊且每個塊在給定其他塊的條件后驗(yàn)分布易于直接抽樣的情況。它輪流對每個參數(shù)塊進(jìn)行抽樣接受率恒為1因此效率通常更高。如今像Stan、PyMC、JAGS這樣的概率編程語言已經(jīng)將復(fù)雜的MCMC采樣過程封裝起來。用戶只需要用類似數(shù)學(xué)公式的語言定義先驗(yàn)和似然軟件會自動選擇高效的采樣算法如NUTS No-U-Turn Sampler并完成抽樣。這極大地降低了貝葉斯建模的門檻。3. 一個完整的實(shí)戰(zhàn)案例估計(jì)廣告點(diǎn)擊率讓我們通過一個具體的例子將上述所有概念串聯(lián)起來。假設(shè)我們運(yùn)營一個網(wǎng)站新上線了一個廣告位。我們根據(jù)行業(yè)經(jīng)驗(yàn)認(rèn)為這個位置的點(diǎn)擊率CTR大概在2%左右但不確定。我們設(shè)定一個先驗(yàn)CTR ~ Beta(α2, β100)。這個Beta分布的均值是α/(αβ)2/102≈1.96%與我們2%的認(rèn)知接近同時它有較寬的分布方差較大表達(dá)了我們的不確定性。然后我們進(jìn)行了一次小規(guī)模測試展示了廣告1000次獲得了15次點(diǎn)擊。我們的數(shù)據(jù)是展示次數(shù)N1000點(diǎn)擊次數(shù)k15。似然函數(shù)很自然地我們假設(shè)每次展示是否點(diǎn)擊是一個伯努利試驗(yàn)?zāi)敲纯偟狞c(diǎn)擊次數(shù)k服從二項(xiàng)分布k ~ Binomial(nN, pCTR)?，F(xiàn)在我們的后驗(yàn)分布是P(CTR | k, N) ∝ P(CTR) * P(k | CTR, N)即Beta(α, β)乘以Binomial(N, CTR)。由于Beta分布是二項(xiàng)分布的共軛先驗(yàn)我們可以直接得到后驗(yàn)分布的解析解CTR | data ~ Beta(α_post α k, β_post β N - k) Beta(215, 1001000-15) Beta(17, 1085)。這個后驗(yàn)分布的均值是17/(171085) ≈ 1.57%?？梢钥吹皆谙闰?yàn)~1.96%和樣本均值15/10001.5%之間后驗(yàn)均值1.57%做了一個加權(quán)平均更靠近數(shù)據(jù)提供的證據(jù)因?yàn)閿?shù)據(jù)量N1000比先驗(yàn)的“等效樣本量”αβ102要大。但是讓我們假裝不知道共軛這個“捷徑”用蒙特卡洛方法來模擬一下。我們使用PyMC庫來實(shí)現(xiàn)import pymc as pm import arviz as az import numpy as np import matplotlib.pyplot as plt # 定義觀測數(shù)據(jù) impressions 1000 clicks 15 # 使用PyMC構(gòu)建模型 with pm.Model() as ctr_model: # 先驗(yàn)分布CTR ~ Beta(2, 100) ctr pm.Beta(ctr, alpha2, beta100) # 似然函數(shù)觀測到的點(diǎn)擊次數(shù) ~ Binomial(impressions, ctr) obs pm.Binomial(obs, nimpressions, pctr, observedclicks) # 使用MCMC采樣默認(rèn)使用NUTS算法 trace pm.sample(draws5000, tune1000, chains4, return_inferencedataTrue) # 后驗(yàn)分析 az.summary(trace) # 查看后驗(yàn)統(tǒng)計(jì)摘要均值、標(biāo)準(zhǔn)差、分位數(shù)等 az.plot_trace(trace) # 繪制軌跡圖檢查收斂性 az.plot_posterior(trace[ctr], ref_val0.015) # 繪制后驗(yàn)分布并與樣本均值(1.5%)比較 plt.show()運(yùn)行這段代碼MCMC采樣器會為我們從后驗(yàn)分布P(CTR | data)中抽取大量樣本。az.summary給出的后驗(yàn)均值、94%最高密度區(qū)間HDI等會與解析解Beta(17, 1085)的結(jié)果非常接近。通過plot_trace我們可以檢查馬爾可夫鏈?zhǔn)欠褚呀?jīng)收斂不同鏈混合良好軌跡像“毛毛蟲”。這個簡單的例子展示了完整的工作流定義先驗(yàn)和似然 - MCMC采樣 - 后驗(yàn)診斷與分析。即使模型復(fù)雜到?jīng)]有共軛解這個流程也完全適用。4. 模型診斷與收斂性判斷你的采樣結(jié)果可信嗎MCMC采樣不是魔法它可能失敗。如果采樣沒有收斂我們得到的樣本就不能代表真正的后驗(yàn)分布基于此做出的任何推斷都是危險的。因此模型診斷是貝葉斯工作流中至關(guān)重要、不可省略的一環(huán)。軌跡圖這是最直觀的診斷工具。將每條馬爾可夫鏈的采樣值按迭代次數(shù)畫出來。一個健康的軌跡圖應(yīng)該看起來像“平穩(wěn)的、毛茸茸的毛毛蟲”——沒有明顯的趨勢在均值附近隨機(jī)波動并且多條鏈的軌跡高度混合、重疊在一起。如果看到明顯的趨勢、周期性或幾條鏈分離很遠(yuǎn)說明采樣沒有收斂或混合性差。自相關(guān)圖它檢查樣本之間的自相關(guān)性。理想情況下隨著滯后階數(shù)的增加自相關(guān)系數(shù)應(yīng)迅速下降到0附近。如果自相關(guān)性很高且衰減很慢意味著采樣效率低下樣本中包含的信息量少你可能需要增加采樣次數(shù)或調(diào)整采樣算法如增加NUTS算法的目標(biāo)接受率。Gelman-Rubin診斷統(tǒng)計(jì)量這是一個量化指標(biāo)。它比較鏈間方差和鏈內(nèi)方差。R?R-hat統(tǒng)計(jì)量越接近1越好。通常R? 1.01被認(rèn)為是收斂的良好標(biāo)志。現(xiàn)代貝葉斯計(jì)算庫如ArviZ會為你計(jì)算這個值。有效樣本量由于MCMC樣本存在自相關(guān)N個樣本并不等同于N個獨(dú)立樣本。ESS衡量的是這些相關(guān)樣本相當(dāng)于多少個獨(dú)立樣本。你至少需要幾百個有效樣本才能對后驗(yàn)做出可靠的估計(jì)。對于尾部概率的估計(jì)則需要更多。實(shí)操心得不要只看默認(rèn)的總結(jié)報告。務(wù)必繪制并仔細(xì)查看軌跡圖和自相關(guān)圖。我曾在一個層次模型上忽略了診斷結(jié)果后驗(yàn)區(qū)間異常地窄導(dǎo)致過于自信的錯誤結(jié)論。后來檢查軌跡圖發(fā)現(xiàn)有一條鏈卡在了低概率區(qū)域。增加老化迭代次數(shù)并調(diào)整參數(shù)化方式后問題才得以解決。一個常見的技巧是如果自相關(guān)性太高可以嘗試對采樣結(jié)果進(jìn)行“稀釋”比如每隔10個樣本取一個但這會浪費(fèi)計(jì)算資源更好的辦法是優(yōu)化模型或使用更高效的采樣器。5. 超越點(diǎn)估計(jì)貝葉斯推斷的威力在于量化不確定性傳統(tǒng)頻率統(tǒng)計(jì)通常給出一個點(diǎn)估計(jì)如均值和一個置信區(qū)間。貝葉斯推斷則直接給出了參數(shù)的整個后驗(yàn)分布。這帶來了幾個無可比擬的優(yōu)勢直接的概率陳述我們可以直接從后驗(yàn)分布中計(jì)算任何感興趣的概率。例如在點(diǎn)擊率案例中我們可以輕松回答“點(diǎn)擊率超過2%的概率是多少” 只需計(jì)算后驗(yàn)樣本中大于0.02的比例。np.mean(trace[ctr] 0.02)。這是頻率學(xué)派的置信區(qū)間無法直接提供的。最高密度區(qū)間HDI是后驗(yàn)分布中一個區(qū)間它包含了指定概率質(zhì)量如94%的參數(shù)值并且區(qū)間內(nèi)的任一點(diǎn)的后驗(yàn)密度都高于區(qū)間外的點(diǎn)。它比基于分位數(shù)的等尾區(qū)間更能反映后驗(yàn)分布的“最可能”范圍尤其當(dāng)后驗(yàn)分布不對稱時。預(yù)測分布貝葉斯的真正目標(biāo)是預(yù)測新數(shù)據(jù)而不是僅僅估計(jì)參數(shù)。我們可以通過后驗(yàn)預(yù)測檢查來模擬新數(shù)據(jù)對于后驗(yàn)分布中的每一個參數(shù)樣本根據(jù)模型生成一組假想的新數(shù)據(jù)。所有這些生成的數(shù)據(jù)的分布就是預(yù)測分布。通過比較預(yù)測分布與實(shí)際觀測數(shù)據(jù)的分布我們可以檢驗(yàn)?zāi)Ｐ偷臄M合優(yōu)度。# 在PyMC中進(jìn)行后驗(yàn)預(yù)測檢查 with ctr_model: # 從后驗(yàn)中抽取預(yù)測樣本 ppc pm.sample_posterior_predictive(trace, extend_inferencedataTrue) # 檢查預(yù)測的點(diǎn)擊次數(shù)分布 ppc_obs ppc.posterior_predictive[obs].values.flatten() plt.hist(ppc_obs, bins30, alpha0.5, densityTrue, labelPredicted) plt.axvline(clicks, colorred, linestyle--, labelObserved (15)) plt.xlabel(Number of Clicks) plt.ylabel(Density) plt.legend() plt.show()如果觀測值紅色虛線落在預(yù)測分布的高概率區(qū)域說明模型能很好地解釋現(xiàn)有數(shù)據(jù)。反之則可能意味著模型設(shè)定有誤。決策分析在商業(yè)環(huán)境中我們最終要基于模型做決策。貝葉斯框架天然地將參數(shù)不確定性傳導(dǎo)至決策結(jié)果。例如我們可以計(jì)算在不同廣告點(diǎn)擊率假設(shè)下預(yù)期收益的后驗(yàn)分布從而選擇期望收益最高的方案并同時了解這個決策所伴隨的風(fēng)險收益的方差或尾部風(fēng)險。6. 進(jìn)階話題與常見陷阱當(dāng)你掌握了基礎(chǔ)后可能會遇到更復(fù)雜的場景和挑戰(zhàn)。層次模型當(dāng)數(shù)據(jù)存在分組結(jié)構(gòu)時如不同用戶、不同城市、不同時間點(diǎn)層次模型或多水平模型非常強(qiáng)大。它假設(shè)每個組有自己的參數(shù)但這些參數(shù)又來自一個共同的群體分布。這能在組間進(jìn)行“部分池化”讓數(shù)據(jù)量少的組從數(shù)據(jù)量大的組“借用”信息從而得到更穩(wěn)健的估計(jì)。在PyMC中這通常通過定義超先驗(yàn)來實(shí)現(xiàn)。先驗(yàn)的敏感性分析你的結(jié)論在多大程度上依賴于先驗(yàn)的選擇一個好的實(shí)踐是進(jìn)行敏感性分析嘗試幾種不同的、合理的先驗(yàn)如無信息先驗(yàn)、弱信息先驗(yàn)、信息性更強(qiáng)的先驗(yàn)觀察后驗(yàn)推斷如均值、HDI是否發(fā)生本質(zhì)變化。如果結(jié)論穩(wěn)定那么你的推斷對先驗(yàn)選擇是穩(wěn)健的如果變化劇烈則需要更謹(jǐn)慎地論證先驗(yàn)的合理性或者承認(rèn)數(shù)據(jù)本身的信息不足以得出強(qiáng)結(jié)論。計(jì)算挑戰(zhàn)與優(yōu)化維度災(zāi)難隨著參數(shù)數(shù)量增加后驗(yàn)分布的空間呈指數(shù)級增長采樣會變得極其困難。對策包括重新參數(shù)化模型如使用非中心化參數(shù)化、使用變分推斷等近似方法作為快速原型或者尋求更專業(yè)的采樣算法。模型錯誤指定如果似然函數(shù)嚴(yán)重偏離數(shù)據(jù)的真實(shí)生成過程那么再精巧的推斷也是徒勞。后驗(yàn)預(yù)測檢查是診斷模型誤設(shè)的重要工具。此外可以嘗試不同的分布族如用學(xué)生t分布替代正態(tài)分布來處理厚尾數(shù)據(jù)。初始化問題MCMC鏈的初始值如果選在低概率區(qū)域可能導(dǎo)致收斂緩慢甚至失敗。一個好的策略是使用從先驗(yàn)分布中抽取的多個隨機(jī)點(diǎn)進(jìn)行多次初始化或者使用變分推斷的結(jié)果作為MCMC的初始值。最后我想分享一個深刻的體會貝葉斯蒙特卡洛建模與其說是一門精確的科學(xué)不如說是一種迭代探索的藝術(shù)。你很少能第一次就設(shè)定出完美的模型。通常的流程是構(gòu)建一個簡單模型 - 運(yùn)行推斷 - 診斷收斂性、預(yù)測檢查- 發(fā)現(xiàn)模型缺陷 - 改進(jìn)模型增加層次結(jié)構(gòu)、改變似然、調(diào)整先驗(yàn)- 再次運(yùn)行推斷。這個循環(huán)可能會重復(fù)多次。擁抱這種不確定性利用后驗(yàn)分布提供的豐富信息進(jìn)行思考和決策才是貝葉斯思維帶給我們的最大價值。它迫使你明確地陳述你的假設(shè)先驗(yàn)誠實(shí)地面對數(shù)據(jù)帶來的更新似然并最終量化你所有認(rèn)知中的不確定性后驗(yàn)。