：從貝葉斯推斷到不確定性量化)
不確定性這個東西做數(shù)據(jù)建模的人天天掛在嘴邊但真正把它當回事、拿進模型里去算的其實少得可憐。大多數(shù)業(yè)務分析還是停留在“點估計”的層面——均值多少、預測值多少、系數(shù)多大頂多再給個置信區(qū)間應付一下。問題是現(xiàn)實世界從來不給你一個確定答案。渠道投放效果會有波動用戶對價格的反應有隨機性供應鏈缺貨風險更是層層疊加。你要是只用單一數(shù)值去描述這些現(xiàn)象等于把所有風險都壓在一個毫無彈性的點上一旦真實值偏離整個決策鏈跟著遭殃。概率編程解決的就是這件事。它讓你把不確定性當成模型的一等公民不是事后貼上去的誤差條而是從數(shù)據(jù)生成機制里就帶出來的東西。用概率編程做不確定性推理本質上是在回答一個更誠實的問題基于當前掌握的數(shù)據(jù)和先驗知識某個未知量到底可能落在哪些值上各自的可能性有多大。這個思路在營銷歸因、推薦系統(tǒng)、風控建模、實驗分析這些場景里都適用而且一旦跑通你會發(fā)現(xiàn)它能帶來的信息量遠不止一個預測數(shù)字那么簡單。這篇文章寫給正在用傳統(tǒng)統(tǒng)計或機器學習方法做建模、想嘗試概率編程但不知道從哪下手的人。我會從概念拆解講到工具選型再用一個完整的實操案例帶你走一遍從數(shù)據(jù)到后驗分布的建模流程最后把我踩過的坑和排查經驗一并整理出來。適合作為你入坑概率編程的第一篇參考。1. 為什么常規(guī)建模方法在不確定性面前不夠用1.1 你平時用的模型其實是在“假裝確定”先做一個思想實驗。你在做某產品的付費轉化率預估用歷史數(shù)據(jù)訓練出一個模型算出轉化率是 2.3%。這個數(shù)字看起來精確但它背后省略了太多信息。當你說轉化率是 2.3% 的時候你其實是在說“我手頭所有數(shù)據(jù)指向的最佳猜測是 2.3%?！笨蓡栴}是這個猜測的置信度有多高如果樣本量只有 500 個用戶這個 2.3% 的可信區(qū)間可能橫跨 1.5% 到 3.5%如果樣本量是 50 萬區(qū)間可能收窄到 2.25% 到 2.35%。同樣的點估計背后代表的不確定性完全不同。傳統(tǒng)頻率派方法不是不能描述這種差異它用置信區(qū)間、標準誤來刻畫。但這里有幾個脫節(jié)的地方。第一置信區(qū)間的解釋方式極其繞它的嚴格含義是“重復抽樣下區(qū)間包含真實值的頻率”而不是“真實值落在這個區(qū)間的概率是 95%”——后者才是業(yè)務決策者真正想問的問題。第二當你引入先驗信息的時候頻率派沒有一套自然的框架來做這件事。第三復雜的層級結構、缺失數(shù)據(jù)、測量誤差這些實際建模中常見的問題用傳統(tǒng)方法處理起來要么得做大量近似假設要么計算復雜到難以下手。概率編程換了一個問法。它直接問你你對某個未知量先驗上有什么判斷觀察到的數(shù)據(jù)長什么樣然后通過貝葉斯公式把這兩部分融合起來輸出一個完整的后驗分布。這個分布就是你要的答案全貌——不光告訴你最大概率的值在哪還告訴你其他值有多大可能。1.2 偶然不確定性和認知不確定性兩個必須分清的概念進入概率編程之前先要把兩種性質完全不同的不確定性區(qū)分清楚因為它們對模型設計的影響是決定性的。偶然不確定性Aleatoric Uncertainty這類不確定性來自系統(tǒng)本身的隨機性是數(shù)據(jù)生成機制自帶的。比如拋硬幣的結果、用戶訪問時長的波動、傳感器讀數(shù)的噪聲。就算你收集一億個樣本這類不確定性不會消失它只是被刻畫得更精確而已。在概率編程里它通常體現(xiàn)在似然函數(shù)的噪聲項上。比如你假設觀測值服從正態(tài)分布那個方差參數(shù)就是在描述偶然不確定性的大小。認知不確定性Epistemic Uncertainty這類不確定性來自你認知的局限是“不知道”造成的。比如你只有 100 個樣本對模型參數(shù)的真實值心里沒底再比如你完全不確定某個渠道的投放效果究竟是正還是負。認知不確定性的特點是可以隨著數(shù)據(jù)積累而減小——樣本越多你對參數(shù)的認識越清晰后驗分布越集中。這兩種不確定性在業(yè)務上有著完全不同的含義。拿庫存預測來說偶然不確定性對應的是需求本身的隨機波動你再怎么優(yōu)化模型也改變不了消費者行為的隨機性能做的是在安全庫存里給它留出緩沖認知不確定性對應的是你對“需求均值到底是多少”這個參數(shù)的認識不清這個可以通過增加數(shù)據(jù)量、做更細致的分層來縮小。如果你在建模時把這兩種不確定性混為一談模型給出的風險評估就是失真的。概率編程天然適合區(qū)分這兩者。通過為模型參數(shù)設置先驗分布并做后驗推斷你可以直觀地看到哪些不確定性在數(shù)據(jù)進來后大幅縮小了——那是認知不確定性在消退哪些沒有縮小——那些就是系統(tǒng)本身的偶然不確定性需要你接受并適應它。這個視角上的轉變對實際決策的影響非常大。2. 概率編程里的核心概念從先驗到后驗的完整鏈路2.1 貝葉斯公式概率編程的地基概率編程的本質是貝葉斯推斷在計算層面的工程化實現(xiàn)。貝葉斯公式表達起來非常簡單$$P(\theta | D) \frac{P(D | \theta) P(\theta)}{P(D)}$$翻譯成大白話$P(\theta)$ 是先驗分布在看到數(shù)據(jù)之前你對參數(shù) $\theta$ 的認知。它可以是信息性的比如“根據(jù)過往經驗某渠道的轉化率大概率在 1%-3% 之間”也可以是弱信息性的比如“轉化率不會小于 0 也不會大于 100%”就行。$P(D | \theta)$ 是似然函數(shù)在給定參數(shù) $\theta$ 的情況下觀測到當前數(shù)據(jù)集 $D$ 的概率。它是模型與數(shù)據(jù)的連接點描述了數(shù)據(jù)生成機制。$P(\theta | D)$ 是后驗分布在看到數(shù)據(jù)之后你更新過的、對 $\theta$ 的最終認知。這是你做推斷、做預測、做決策的依據(jù)。理解貝葉斯公式的關鍵在于“更新”這個詞。先驗是你最初的判斷數(shù)據(jù)像秤砣一樣改變著判斷的天平后驗是了這個過程之后的結果。數(shù)據(jù)量越大、數(shù)據(jù)越“硬”后驗就越被數(shù)據(jù)主導先驗的影響越弱數(shù)據(jù)量小的時候后驗則主要是先驗在支撐——這也是貝葉斯方法在小樣本場景下依然能給出合理推斷的原因。概率編程做的事情用一句話概括就是你定義先驗和似然它幫你算后驗。這里的核心難點是$P(D)$ 這個歸一化常數(shù)在絕大多數(shù)現(xiàn)實模型里是沒有解析解的你需要用數(shù)值方法來逼近后驗分布。這就是 MCMC 和變分推斷登場的時刻。2.2 后驗推斷的兩條路MCMC 與變分推斷**MCMC馬爾可夫鏈蒙特卡洛**是目前概率編程中最主流的后驗推斷方法。它的思路很巧妙既然無法直接算出后驗的解析形式那就構造一種采樣策略讓采樣出來的樣本服從的分布逐漸逼近后驗分布。具體實現(xiàn)的算法很多比較常用的是 No-U-Turn SamplerNUTS它是 Hamilton Monte CarloHMC的一種自適應變體。NUTS 能在參數(shù)空間里利用梯度信息高效地探索后驗分布在參數(shù)維度不高、后驗形態(tài)不算太扭曲的情況下收斂速度相當理想。使用 MCMC 得到的是一堆從后驗分布中抽取的樣本。通過對這些樣本做統(tǒng)計你可以得到后驗均值、分位數(shù)、可信區(qū)間等一切你關心的量。樣本量越大對后驗的刻畫就越精細。它的優(yōu)點是準確——沒有系統(tǒng)性偏差只要鏈長足夠、診斷通過樣本分布就趨近于真實后驗。缺點是需要大量計算尤其在高維參數(shù)空間或大數(shù)據(jù)集場景下采樣成本可能高到不可接受。**變分推斷Variational Inference**是另一個方向。它不是采樣而是把后驗推斷轉化成一個最優(yōu)化問題設定一個參數(shù)化的分布族 $q(\theta; \phi)$然后調整變分參數(shù) $\phi$讓 $q$ 盡可能接近真實后驗 $P(\theta | D)$通常用 KL 散度來衡量兩個分布的距離。優(yōu)化完成后$q$ 就作為后驗的近似輸出。變分推斷的優(yōu)勢是速度快能處理大規(guī)模數(shù)據(jù)和復雜模型。缺點是有偏——$q$ 的分布族選擇如果表達能力不夠近似結果就和真實后驗有偏差?，F(xiàn)實使用中一般是這樣取舍的先期探索模型結構、需要快速迭代的時候用變分推斷最終關鍵結論的確定、當模型維度不算太高的時候用 MCMC。我在實際項目里經常兩個一起跑變分推斷給出的結果當“預報員”MCMC 的結果當“裁判員”。2.3 概率編程語言到底“編程”了什么很多人第一次看到 PyMC 或 Stan 的代碼時會有個困惑這不就是寫個模型公式然后調 sample() 嗎這算什么編程這個問題的答案在于概率編程語言提供的是一種生成式模型的表達框架。你寫的模型代碼實際上是在描述“數(shù)據(jù)是如何被生成出來的”——先定義參數(shù)的先驗分布再定義中間變量之間的確定性關系最后定義觀測變量的似然分布?？蚣茇撠熢诤笈_把這個描述轉換成對應的計算圖然后自動推導采樣所需的梯度信息、自動完成推斷過程。這種“聲明式建模”的方式帶來兩個巨大的好處。第一模型和推斷解耦你可以把精力全部放在模型結構的合理性上而不需要自己實現(xiàn)采樣算法、推導梯度公式。第二模型靈活度極高你可以輕松地構建多層級的回歸模型、帶有測量誤差的模型、混合分布模型等復雜結構這些用傳統(tǒng)統(tǒng)計軟件做會非常痛苦但在概率編程框架里只是多寫幾行代碼的事。當然這并不意味著你可以完全不懂推斷原理。實際上對 MCMC 的收斂性、采樣器行為、先驗敏感性有基本了解是排查模型問題的前提——這個問題我在后面的章節(jié)會詳細展開。3. 工具選型主流概率編程框架怎么選3.1 PyMC、Stan、NumPyro 三足鼎立概率編程的生態(tài)經過這么多年的發(fā)展基本形成了三足鼎立的局面老牌的 PyMC前身是 PyMC2/PyMC3、工業(yè)級口碑的 Stan尤其它的 R 接口、以及 JAX 生態(tài)下的新銳 NumPyro。選型沒有絕對的正確答案關鍵看你的技術棧和需求。PyMCPython 生態(tài)中最成熟、社區(qū)最活躍的選擇。API 設計對新手非常友好內置的pm.Model()上下文管理器讓模型定義像寫公式一樣自然。它封裝了多種采樣器和變分推斷方法pm.sample()一把梭就能完成推斷。配合arviz庫做后驗診斷和可視化體驗非常順滑。如果你的主力語言是 Python想快速上手、快速驗證想法PyMC 基本是首選。Stan以 C 為基礎實現(xiàn)通過cmdstanpy或 R 的rstan接口調用。Stan 的建模語言是獨立的 Stan Language需要在.stan文件中寫模型然后編譯運行。它的 HMC/NUTS 實現(xiàn)非常扎實在某些復雜模型和后驗幾何相對崎嶇的情況下采樣效率往往優(yōu)于其他框架。代價是建模和調試的迭代周期比 Python 原生環(huán)境要長一些。適合需要極致采樣性能、模型復雜度高的場景。NumPyro基于 JAX 構建用函數(shù)式編程風格來描述模型。它的最大優(yōu)勢是支持 JIT 編譯和 GPU 加速在對性能極度敏感的大規(guī)模任務上有明顯優(yōu)勢。如果你是 JAX 生態(tài)的玩家或者模型大到采樣不得不靠 GPU 加速NumPyro 非常值得嘗試。學習曲線相對陡——函數(shù)式風格、numpyro.sample的調用方式和plate與numpyro.handlers的交互邏輯都需要一些適應時間。3.2 我的選型建議按場景對號入座這里給你一個比較實際的參考場景推薦工具理由Python 新手 / 快速原型PyMC上手快、社區(qū)資料多、可視化配套完善復雜模型 / 需要最高采樣質量StanNUTS 實現(xiàn)成熟對崎嶇后驗的適應性強大數(shù)據(jù)量 / GPU 加速 / JAX 生態(tài)NumPyroJIT 向量化性能天花板最高R 用戶做統(tǒng)計建模rstan與 R 的統(tǒng)計生態(tài)銜接自然生產環(huán)境部署NumPyro / PyMC 序列化Python 棧部署路徑最順還有兩個輕量級選手值得一提。sampyl和emcee都是輕量級 MCMC 庫體量小、依賴少適合教學和小型項目。但它們功能覆蓋面有限生產級建模我還是建議回到三大件上來。我的個人經驗是不要貪多主攻一個熟悉另一個。比如你日常用 PyMC 建模遇到采樣性能瓶頸時能切到 Stan 去對比驗證或者你在 PyMC 里原型驗證沒問題之后用 NumPyro 在生產環(huán)境里重寫一遍上 GPU?？蚣苤g模型表達的思路是相通的重點是把核心概念吃透。4. 完整實操用概率編程做一個轉化率不確定性分析4.1 業(yè)務場景與數(shù)據(jù)構造我們來走一個完整的建模流程。場景設定某產品線做了一個新版本的功能產品經理想知道新版上線后購買轉化率到底比舊版高多少。收集到的數(shù)據(jù)形式很簡單舊版觀察了 2000 個用戶其中 96 人完成購買新版觀察了 2140 個用戶其中 132 人完成購買。傳統(tǒng)做法是算轉換率、做個卡方檢驗或者兩比例 z 檢驗給個 p 值。但概率編程能給出的信息豐富得多你能直接得到轉化率差異的完整后驗分布并基于概率做出決策比如“新版轉化率比舊版高的概率是 93.7%”。為了演示我們用 Python 的 PyMC 框架。核心代碼如下import pymc as pm import arviz as az import numpy as np # 觀測數(shù)據(jù) n_old, x_old 2000, 96 n_new, x_new 2140, 132 with pm.Model() as conversion_model: # 先驗舊版本轉化率 ~ Beta(1, 1) 等價于 Uniform(0, 1) p_old pm.Beta(p_old, alpha1, beta1) p_new pm.Beta(p_new, alpha1, beta1) # 似然觀測到的購買人數(shù)服從二項分布 obs_old pm.Binomial(obs_old, nn_old, pp_old, observedx_old) obs_new pm.Binomial(obs_new, nn_new, pp_new, observedx_new) # 派生量轉化率差異 diff pm.Deterministic(diff, p_new - p_old) # 執(zhí)行 MCMC 采樣 trace pm.sample(draws2000, tune2000, chains4, random_seed42)這段代碼一共做了四件事定義轉化率的先驗Beta(1,1) 表示完全沒有任何先驗傾向、定義觀測數(shù)據(jù)的似然二項分布、定義關心的派生量轉化率差異、執(zhí)行 NUTS 采樣。跑完之后trace里存的就是所有參數(shù)的 MCMC 樣本。4.2 后驗分析從樣本到業(yè)務決策采樣完成后我們把后驗樣本展示出來az.summary(trace, var_names[p_old, p_new, diff])你會得到一張表格每一行一個變量包含均值、標準差、HDI最高密度區(qū)間、有效樣本量等指標。假設輸出結果是mean sd hdi_3% hdi_97% ess_bulk r_hat p_old 0.048 0.003 0.042 0.054 6234 1.00 p_new 0.062 0.004 0.055 0.069 6157 1.00 diff 0.013 0.005 0.004 0.022 6128 1.00這份結果信息的密度遠超一個 p 值。首先p_old和p_new的 HDI 區(qū)間都是兩側遠離 0 的表明數(shù)據(jù)充分支持這兩個轉化率確實在 4%-7% 這個量級。其次diff的 94% HDI 區(qū)間是 [0.004, 0.022]不包含 0——這說明新版轉化率顯著高于舊版等等先別急著說“顯著”在貝葉斯框架里我們換一種更直白的表述。更直觀的做法是直接計算概率prob_improvement (trace.posterior[diff] 0).mean().item() print(f新版轉化率高于舊版的概率: {prob_improvement:.1%})假設輸出為 98.3%那這個信息可以直接扔給業(yè)務方我們不需要糾結“顯著性水平 0.05”或“p0.03”這種解釋易產生歧義的表述而是可以說“根據(jù)數(shù)據(jù)表現(xiàn)新版轉化率高于舊版的概率為 98.3%”。這個概率表述業(yè)務方一聽就懂決策者可以直接基于它來評估是否全量上線。4.3 先驗敏感性分析給結論上雙保險有經驗的建模者不會只看一組先驗跑出來的結果就下結論。貝葉斯方法的先驗選擇在一定程度上是主觀判斷所以要做敏感性分析換不同的先驗設置看結論會不會被顯著改變。我們把代碼改成用信息性強的先驗再跑一遍with pm.Model() as conversion_model_prior: # 更強的先驗根據(jù)歷史經驗轉化率集中在 3%-7% 之間 p_old pm.Beta(p_old, alpha10, beta170) p_new pm.Beta(p_new, alpha10, beta170) obs_old pm.Binomial(obs_old, nn_old, pp_old, observedx_old) obs_new pm.Binomial(obs_new, nn_new, pp_new, observedx_new) diff pm.Deterministic(diff, p_new - p_old) trace_prior pm.sample(draws2000, tune2000, chains4, random_seed42)對比兩次運行的diff后驗均值如果結論基本一致說明你的結論受先驗影響不大、是數(shù)據(jù)本身在驅動如果差異明顯說明當前數(shù)據(jù)量還不足以克服先驗的影響。這時候就要慎重對結論的表述要加上“基于當前先驗判斷”的限定條件。這是貝葉斯建模的好習慣也是規(guī)避“先驗挑數(shù)據(jù)”這個爭議的有效方式。4.4 層級模型的擴展不滿足于單一體轉化率上述模型是極其簡化的——它假設所有用戶同質忽略用戶來源渠道、訪問設備、時間區(qū)間等因素。一個更接近真實業(yè)務的建模方式是把這些變量吸納進模型。設想數(shù)據(jù)升級為用戶來源渠道有 A、B、C 三個每個渠道既有舊版也有新版訪問。你可以在模型里加入渠道作為分組變量構建一個層級模型# 假設 channels 是 0,1,2 編碼的渠道標簽 import pandas as pd # 這里簡化為構造一個示例 DataFrame np.random.seed(42) df pd.DataFrame({ channel: np.repeat([0, 1, 2], [800, 700, 640]), version: 0 # 占位 }) # 實際使用時數(shù)據(jù)是每個 [channel, version] 組合的訪問量和購買量 with pm.Model() as hierarchical_model: # 渠道間共享的超先驗 mu pm.Normal(mu, mu0, sigma1) sigma pm.HalfCauchy(sigma, beta0.5) # 渠道效應收縮到共同均值 z_channel pm.Normal(z_channel, mu0, sigma1, shape3) beta_channel pm.Deterministic(beta_channel, mu sigma * z_channel) # 版本效應按渠道建模時也能分層處理 p_logit mu beta_channel * ... # 根據(jù)具體編碼設計 # 似然觀測銷售數(shù)據(jù)服從伯努利/二項分布 ...層級模型的好處在于部分池化樣本量少的渠道會參考整體分布的信息避免因樣本稀疏產生極端的估計樣本量大的渠道則保持自己的特性。這在轉化率、客單價等業(yè)務指標的多維度分析中非常實用也是傳統(tǒng)頻率派方法處理起來最棘手、而概率編程最舒展的場景之一。5. 建模實操中的常見問題與排查經驗5.1 采樣不收斂最讓人抓狂的 MCMC 問題跑pm.sample()之后第一件事永遠是看r_hat和有效樣本量ess_bulk。r_hat大于 1.01 基本就意味著采樣鏈沒有收斂。常見原因和排查手段我整理成了一張速查表癥狀可能原因排查方向r_hat 恒大于 1.01鏈初始值差異太大burn-in 不夠增加 tune 長度提供更合理的初始值鏈呈“毛蟲”狀、游走緩慢后驗分布相關性高 / 參數(shù)化不佳嘗試標準化參數(shù)、非中心化參數(shù)化高維模型頻繁發(fā)散后驗幾何太崎嶇調整步長、換用pm.Metropolis對比、簡化模型結構eff 樣本量極低500自相關性過高增加chains或draws檢查模型是否有冗余參數(shù)警告出現(xiàn)很多divergencesNUTS 步長自適應遇到問題通常因先驗與數(shù)據(jù)沖突檢查先驗是否有不合理取值并調整或降低target_accept到 0.9有參數(shù)后驗分布“掛在邊界”先驗約束不夠或者數(shù)據(jù)對參數(shù)的識別力弱考慮更換先驗分布族例如在比例參數(shù)上使用 Beta 而不是 Uniform我遇到最多的還是參數(shù)化不當導致的無效采樣。一個典型的例子是方差參數(shù)用了HalfNormal(sigma10)這樣過寬的弱先驗。表面上看起來“不預設立場很客觀”但實際上過寬的先驗會在高參數(shù)空間分配大量概率質量導致采樣器在極端值之間掙扎表現(xiàn)為大量發(fā)散。改用HalfCauchy或適度窄的先驗分布之后采樣立刻恢復順暢。這里想多說一句弱先驗不等于打死不設先驗。在你對參數(shù)的合理范圍多少有數(shù)的情況下給一個符合業(yè)務常識的約束或弱信息先驗往往既不會扭曲后驗結論又能顯著改善采樣性能。這也是貝葉斯工程師的經驗之談。5.2 模型結果怪得離譜先檢查似然和數(shù)據(jù)的對應關系有次我做一個客單價預估模型后驗均值明顯偏離業(yè)務直覺。排查了半天最后發(fā)現(xiàn)是“數(shù)據(jù)對齊”的問題訓練數(shù)據(jù)里的一列應是當日訂單量結果因為某次處理腳本的 bug混入了前一天的訂單數(shù)據(jù)。這類“模型以外”的錯誤在統(tǒng)計推斷里比模型本身的錯誤更隱蔽。排查順序建議為檢查數(shù)據(jù)描述性統(tǒng)計——均值、分位數(shù)、缺失率是否合理。畫出關鍵變量在模型下的生成模擬數(shù)據(jù)posterior_predictive對比真實數(shù)據(jù)的分布。如果模擬數(shù)據(jù)和觀測數(shù)據(jù)差別巨大說明模型的一個或多個核心假設有問題。逐一簡化模型先跑一個只含主要效應、不含交互項和層級結構的模型看看能否得到大概合理的后驗再逐步增加復雜度。對連續(xù)型變量檢查數(shù)據(jù)縮放是否均勻。數(shù)值范圍跨度過大的特征會讓 NUTS 的步長自適應失靈建議對輸入數(shù)據(jù)做標準化同時可以在模型的確定性關系中使用相應的偏差。曾有很長一段時間我在做這類分析時從不做第 2 步——posterior_predictive檢查——后來一個實際案例教育了我模型對某個分組變量的層級方差估得很小把分組間差異全都“解釋”成了偶然噪聲完全不符合業(yè)務預期。當時的模型層設計沒有包含組內樣本量對先前方差收縮作用的權衡簡單來說就是層級結構設置不當。如果沒有事后預測檢查這個問題會一直潛伏到部署階段才爆發(fā)。5.3 變分推斷的結果不靠譜別急著怪算法使用變分推斷時最容易掉進的坑是q的分布族表達能力不足導致求出的近似后驗和真實后驗差異巨大但收斂過程看起來完全正常。一種驗證方法是把變分推斷結果和 MCMC 結果做對比。在數(shù)據(jù)量不大、模型不太復雜的時候這個對比成本有限卻能一針見血地發(fā)現(xiàn)問題。如果兩者差異在可接受范圍內就用變分推斷如果差異巨大檢查是不是原始模型的非光滑環(huán)節(jié)比如離散變量、硬截斷導致變分分布族嚴重失配。另外變分推斷對初始化也比較敏感多試幾組隨機種子如果結果的不確定性變動異常說明優(yōu)化過程本身可能有問題。從工程實踐角度看我的策略是變分推斷用在探索階段MCMC 用在發(fā)布階段。每次項目保質期內的最終結論都以 MCMC 驗證過的版本為準。5.4 業(yè)務落地時的不確定性傳播最后分享一個經常被忽略但實際影響很大的問題——概率編程的輸出如何進入下游系統(tǒng)。一種常見做法是只把后驗均值傳到下游這等于又把不確定性丟掉了。較好的做法是把后驗分布的完整樣本傳遞過去。具體而言如果你在做庫存優(yōu)化不是傳一個“預測需求 1300 件”而是傳一組后驗樣本“需求可能分布在 [1180, 1420]”下游的庫存優(yōu)化代碼可以對每一個樣本跑一遍最優(yōu)庫存決策得到一整套決策的分布最終選擇覆蓋某個分位的庫存量。如果你在做 A/B 實驗的后續(xù)部署判斷不是簡單用“贏了的版本全量上”而是用后驗差異的分布來確定上線節(jié)奏和回滾閾值。在做預算分配時可以把每個渠道的收益參數(shù)后驗樣本輸入優(yōu)化器這樣得到的分配方案本身就是帶不確定性的決策者可以看到最壞情況和最好情況之間的差異。很多業(yè)務場景本身并不缺少優(yōu)化方法缺少的是給優(yōu)化方法提供不確定性輸入的意識。概率編程給出了這個輸入它是從純粹統(tǒng)計推斷到實際業(yè)務決策的橋梁。6. 實戰(zhàn)進階從簡單模型走向生產級應用6.1 從單模型到貝葉斯工作流把概率編程真正用起來之后你會發(fā)現(xiàn)自己需要的不只是一個采樣函數(shù)而是一整套工作流?；谖业膶崙?zhàn)經驗一個可復用的貝葉斯建模流程應該是這樣的問題形式化把業(yè)務問題翻譯成一個可以計算概率的可量化問題。明確哪些是觀測變量、哪些是待推斷參數(shù)、哪些是派生業(yè)務指標。數(shù)據(jù)審查做描述性統(tǒng)計和可視化識別數(shù)據(jù)缺失、異常值和量綱問題。先驗設立基于業(yè)務常識設定合理的弱信息先驗。向業(yè)務方了解“過去的歷史水平”讓它成為先驗均值的重要參考。模型搭建與迭代從簡單模型開始驗證核心結構正確后再加復雜度。每次增加結構后都用后驗預測檢查來檢驗改進效果。推斷與診斷執(zhí)行 MCMC或變分推斷檢查r_hat、ess_bulk、發(fā)散情況。必要時做參數(shù)化調整。敏感性分析改變先驗設定、修改數(shù)據(jù)子集做重跑確認結論的穩(wěn)定性。產出與溝通把后驗分布翻譯成業(yè)務語言輸出“差異為 X 的概率是多少”、“落在哪個區(qū)間的把握有多大”這類可操作的結論。部署與監(jiān)控把模型固化為可重復執(zhí)行的推斷管道對新數(shù)據(jù)周期性地更新后驗并監(jiān)控分布漂移。這八步聽起來多但每一步的時間消耗并不大。熟練之后一個中等復雜度的模型從開始到出結論基本在一兩天內就能完成。6.2 讓模型可復現(xiàn)的幾個工程習慣概率編程的模型天然具有隨機性因此可復現(xiàn)性是工程化繞不開的問題。沒有固定隨機種子的模型運行每次結果都有細微差異這在調試時會讓人分不清改動模型前的效果差異是真實的還是隨機波動造成的。我現(xiàn)在在項目里嚴格遵守以下三條規(guī)則固定隨機種子所有采樣、變分推斷調用都顯式傳入random_seed并使它在每次運行中保持一致。記錄模型版本和依賴環(huán)境把模型代碼、依賴庫版本號、數(shù)據(jù)文件哈希寫進自動生成的模型元數(shù)據(jù)文件。后期排查問題這一條能救命。用緩存機制保存推斷結果同一個模型、同一份數(shù)據(jù)跑過一次后把 trace 序列化為文件。后續(xù)所有人的分析都基于這份緩存結果保證溝通討論的前后一致。這些習慣沒有任何高深之處但它們直接決定了概率編程能否從小項目走向正式工作流。7. 一些關于概率編程的實話做概率編程這些年最大的感受是這項技術真正改變人的地方不在采樣算法的效率而在思維方式的轉換。當你開始用分布去思考所有問題你會意識到以前那些點估計的結論有多少是被“平均”掩埋掉的。但也要說句公道話概率編程不是萬能的。它的計算成本遠高于傳統(tǒng)統(tǒng)計方法在超大樣本場景下需要投入相當?shù)臅r間去做工程優(yōu)化。如果業(yè)務場景的數(shù)據(jù)量足夠大、模型相對簡單、只需要一個穩(wěn)定輸出傳統(tǒng)方法可能更實惠。概率編程適合的是對不確定性本身有決策價值、數(shù)據(jù)規(guī)?？煽?、模型復雜度中上的場景——用一句直白的話說當你發(fā)現(xiàn)“上線后預測總偏差大、又不知道偏差從哪來”的時候概率編程大概能給你答案。另外概率編程也不該被你當成一個“萬能答案機”。它給出的后驗分布依然建立在模型假設之上。假設錯了后驗再漂亮也不可靠。這也是為什么做概率編程的人一半時間在建模另一半時間在驗證假設、檢查模型與數(shù)據(jù)的匹配度。這些約束和代價不是說概率編程不好而是說你有必要想清楚我到底需要的是不確定性信息還是僅僅需要結論的確定性去完成眼前這一步。想明白了再上手這個工具給你的回報會遠高于你的預期。 概率編程實戰(zhàn)：從貝葉斯推斷到不確定性量化
不確定性這個東西做數(shù)據(jù)建模的人天天掛在嘴邊但真正把它當回事、拿進模型里去算的其實少得可憐。大多數(shù)業(yè)務分析還是停留在“點估計”的層面——均值多少、預測值多少、系數(shù)多大頂多再給個置信區(qū)間應付一下。問題是現(xiàn)實世界從來不給你一個確定答案。渠道投放效果會有波動用戶對價格的反應有隨機性供應鏈缺貨風險更是層層疊加。你要是只用單一數(shù)值去描述這些現(xiàn)象等于把所有風險都壓在一個毫無彈性的點上一旦真實值偏離整個決策鏈跟著遭殃。概率編程解決的就是這件事。它讓你把不確定性當成模型的一等公民不是事后貼上去的誤差條而是從數(shù)據(jù)生成機制里就帶出來的東西。用概率編程做不確定性推理本質上是在回答一個更誠實的問題基于當前掌握的數(shù)據(jù)和先驗知識某個未知量到底可能落在哪些值上各自的可能性有多大。這個思路在營銷歸因、推薦系統(tǒng)、風控建模、實驗分析這些場景里都適用而且一旦跑通你會發(fā)現(xiàn)它能帶來的信息量遠不止一個預測數(shù)字那么簡單。這篇文章寫給正在用傳統(tǒng)統(tǒng)計或機器學習方法做建模、想嘗試概率編程但不知道從哪下手的人。我會從概念拆解講到工具選型再用一個完整的實操案例帶你走一遍從數(shù)據(jù)到后驗分布的建模流程最后把我踩過的坑和排查經驗一并整理出來。適合作為你入坑概率編程的第一篇參考。1. 為什么常規(guī)建模方法在不確定性面前不夠用1.1 你平時用的模型其實是在“假裝確定”先做一個思想實驗。你在做某產品的付費轉化率預估用歷史數(shù)據(jù)訓練出一個模型算出轉化率是 2.3%。這個數(shù)字看起來精確但它背后省略了太多信息。當你說轉化率是 2.3% 的時候你其實是在說“我手頭所有數(shù)據(jù)指向的最佳猜測是 2.3%?！笨蓡栴}是這個猜測的置信度有多高如果樣本量只有 500 個用戶這個 2.3% 的可信區(qū)間可能橫跨 1.5% 到 3.5%如果樣本量是 50 萬區(qū)間可能收窄到 2.25% 到 2.35%。同樣的點估計背后代表的不確定性完全不同。傳統(tǒng)頻率派方法不是不能描述這種差異它用置信區(qū)間、標準誤來刻畫。但這里有幾個脫節(jié)的地方。第一置信區(qū)間的解釋方式極其繞它的嚴格含義是“重復抽樣下區(qū)間包含真實值的頻率”而不是“真實值落在這個區(qū)間的概率是 95%”——后者才是業(yè)務決策者真正想問的問題。第二當你引入先驗信息的時候頻率派沒有一套自然的框架來做這件事。第三復雜的層級結構、缺失數(shù)據(jù)、測量誤差這些實際建模中常見的問題用傳統(tǒng)方法處理起來要么得做大量近似假設要么計算復雜到難以下手。概率編程換了一個問法。它直接問你你對某個未知量先驗上有什么判斷觀察到的數(shù)據(jù)長什么樣然后通過貝葉斯公式把這兩部分融合起來輸出一個完整的后驗分布。這個分布就是你要的答案全貌——不光告訴你最大概率的值在哪還告訴你其他值有多大可能。1.2 偶然不確定性和認知不確定性兩個必須分清的概念進入概率編程之前先要把兩種性質完全不同的不確定性區(qū)分清楚因為它們對模型設計的影響是決定性的。偶然不確定性Aleatoric Uncertainty這類不確定性來自系統(tǒng)本身的隨機性是數(shù)據(jù)生成機制自帶的。比如拋硬幣的結果、用戶訪問時長的波動、傳感器讀數(shù)的噪聲。就算你收集一億個樣本這類不確定性不會消失它只是被刻畫得更精確而已。在概率編程里它通常體現(xiàn)在似然函數(shù)的噪聲項上。比如你假設觀測值服從正態(tài)分布那個方差參數(shù)就是在描述偶然不確定性的大小。認知不確定性Epistemic Uncertainty這類不確定性來自你認知的局限是“不知道”造成的。比如你只有 100 個樣本對模型參數(shù)的真實值心里沒底再比如你完全不確定某個渠道的投放效果究竟是正還是負。認知不確定性的特點是可以隨著數(shù)據(jù)積累而減小——樣本越多你對參數(shù)的認識越清晰后驗分布越集中。這兩種不確定性在業(yè)務上有著完全不同的含義。拿庫存預測來說偶然不確定性對應的是需求本身的隨機波動你再怎么優(yōu)化模型也改變不了消費者行為的隨機性能做的是在安全庫存里給它留出緩沖認知不確定性對應的是你對“需求均值到底是多少”這個參數(shù)的認識不清這個可以通過增加數(shù)據(jù)量、做更細致的分層來縮小。如果你在建模時把這兩種不確定性混為一談模型給出的風險評估就是失真的。概率編程天然適合區(qū)分這兩者。通過為模型參數(shù)設置先驗分布并做后驗推斷你可以直觀地看到哪些不確定性在數(shù)據(jù)進來后大幅縮小了——那是認知不確定性在消退哪些沒有縮小——那些就是系統(tǒng)本身的偶然不確定性需要你接受并適應它。這個視角上的轉變對實際決策的影響非常大。2. 概率編程里的核心概念從先驗到后驗的完整鏈路2.1 貝葉斯公式概率編程的地基概率編程的本質是貝葉斯推斷在計算層面的工程化實現(xiàn)。貝葉斯公式表達起來非常簡單$$P(\theta | D) \frac{P(D | \theta) P(\theta)}{P(D)}$$翻譯成大白話$P(\theta)$ 是先驗分布在看到數(shù)據(jù)之前你對參數(shù) $\theta$ 的認知。它可以是信息性的比如“根據(jù)過往經驗某渠道的轉化率大概率在 1%-3% 之間”也可以是弱信息性的比如“轉化率不會小于 0 也不會大于 100%”就行。$P(D | \theta)$ 是似然函數(shù)在給定參數(shù) $\theta$ 的情況下觀測到當前數(shù)據(jù)集 $D$ 的概率。它是模型與數(shù)據(jù)的連接點描述了數(shù)據(jù)生成機制。$P(\theta | D)$ 是后驗分布在看到數(shù)據(jù)之后你更新過的、對 $\theta$ 的最終認知。這是你做推斷、做預測、做決策的依據(jù)。理解貝葉斯公式的關鍵在于“更新”這個詞。先驗是你最初的判斷數(shù)據(jù)像秤砣一樣改變著判斷的天平后驗是了這個過程之后的結果。數(shù)據(jù)量越大、數(shù)據(jù)越“硬”后驗就越被數(shù)據(jù)主導先驗的影響越弱數(shù)據(jù)量小的時候后驗則主要是先驗在支撐——這也是貝葉斯方法在小樣本場景下依然能給出合理推斷的原因。概率編程做的事情用一句話概括就是你定義先驗和似然它幫你算后驗。這里的核心難點是$P(D)$ 這個歸一化常數(shù)在絕大多數(shù)現(xiàn)實模型里是沒有解析解的你需要用數(shù)值方法來逼近后驗分布。這就是 MCMC 和變分推斷登場的時刻。2.2 后驗推斷的兩條路MCMC 與變分推斷**MCMC馬爾可夫鏈蒙特卡洛**是目前概率編程中最主流的后驗推斷方法。它的思路很巧妙既然無法直接算出后驗的解析形式那就構造一種采樣策略讓采樣出來的樣本服從的分布逐漸逼近后驗分布。具體實現(xiàn)的算法很多比較常用的是 No-U-Turn SamplerNUTS它是 Hamilton Monte CarloHMC的一種自適應變體。NUTS 能在參數(shù)空間里利用梯度信息高效地探索后驗分布在參數(shù)維度不高、后驗形態(tài)不算太扭曲的情況下收斂速度相當理想。使用 MCMC 得到的是一堆從后驗分布中抽取的樣本。通過對這些樣本做統(tǒng)計你可以得到后驗均值、分位數(shù)、可信區(qū)間等一切你關心的量。樣本量越大對后驗的刻畫就越精細。它的優(yōu)點是準確——沒有系統(tǒng)性偏差只要鏈長足夠、診斷通過樣本分布就趨近于真實后驗。缺點是需要大量計算尤其在高維參數(shù)空間或大數(shù)據(jù)集場景下采樣成本可能高到不可接受。**變分推斷Variational Inference**是另一個方向。它不是采樣而是把后驗推斷轉化成一個最優(yōu)化問題設定一個參數(shù)化的分布族 $q(\theta; \phi)$然后調整變分參數(shù) $\phi$讓 $q$ 盡可能接近真實后驗 $P(\theta | D)$通常用 KL 散度來衡量兩個分布的距離。優(yōu)化完成后$q$ 就作為后驗的近似輸出。變分推斷的優(yōu)勢是速度快能處理大規(guī)模數(shù)據(jù)和復雜模型。缺點是有偏——$q$ 的分布族選擇如果表達能力不夠近似結果就和真實后驗有偏差?，F(xiàn)實使用中一般是這樣取舍的先期探索模型結構、需要快速迭代的時候用變分推斷最終關鍵結論的確定、當模型維度不算太高的時候用 MCMC。我在實際項目里經常兩個一起跑變分推斷給出的結果當“預報員”MCMC 的結果當“裁判員”。2.3 概率編程語言到底“編程”了什么很多人第一次看到 PyMC 或 Stan 的代碼時會有個困惑這不就是寫個模型公式然后調 sample() 嗎這算什么編程這個問題的答案在于概率編程語言提供的是一種生成式模型的表達框架。你寫的模型代碼實際上是在描述“數(shù)據(jù)是如何被生成出來的”——先定義參數(shù)的先驗分布再定義中間變量之間的確定性關系最后定義觀測變量的似然分布?？蚣茇撠熢诤笈_把這個描述轉換成對應的計算圖然后自動推導采樣所需的梯度信息、自動完成推斷過程。這種“聲明式建模”的方式帶來兩個巨大的好處。第一模型和推斷解耦你可以把精力全部放在模型結構的合理性上而不需要自己實現(xiàn)采樣算法、推導梯度公式。第二模型靈活度極高你可以輕松地構建多層級的回歸模型、帶有測量誤差的模型、混合分布模型等復雜結構這些用傳統(tǒng)統(tǒng)計軟件做會非常痛苦但在概率編程框架里只是多寫幾行代碼的事。當然這并不意味著你可以完全不懂推斷原理。實際上對 MCMC 的收斂性、采樣器行為、先驗敏感性有基本了解是排查模型問題的前提——這個問題我在后面的章節(jié)會詳細展開。3. 工具選型主流概率編程框架怎么選3.1 PyMC、Stan、NumPyro 三足鼎立概率編程的生態(tài)經過這么多年的發(fā)展基本形成了三足鼎立的局面老牌的 PyMC前身是 PyMC2/PyMC3、工業(yè)級口碑的 Stan尤其它的 R 接口、以及 JAX 生態(tài)下的新銳 NumPyro。選型沒有絕對的正確答案關鍵看你的技術棧和需求。PyMCPython 生態(tài)中最成熟、社區(qū)最活躍的選擇。API 設計對新手非常友好內置的pm.Model()上下文管理器讓模型定義像寫公式一樣自然。它封裝了多種采樣器和變分推斷方法pm.sample()一把梭就能完成推斷。配合arviz庫做后驗診斷和可視化體驗非常順滑。如果你的主力語言是 Python想快速上手、快速驗證想法PyMC 基本是首選。Stan以 C 為基礎實現(xiàn)通過cmdstanpy或 R 的rstan接口調用。Stan 的建模語言是獨立的 Stan Language需要在.stan文件中寫模型然后編譯運行。它的 HMC/NUTS 實現(xiàn)非常扎實在某些復雜模型和后驗幾何相對崎嶇的情況下采樣效率往往優(yōu)于其他框架。代價是建模和調試的迭代周期比 Python 原生環(huán)境要長一些。適合需要極致采樣性能、模型復雜度高的場景。NumPyro基于 JAX 構建用函數(shù)式編程風格來描述模型。它的最大優(yōu)勢是支持 JIT 編譯和 GPU 加速在對性能極度敏感的大規(guī)模任務上有明顯優(yōu)勢。如果你是 JAX 生態(tài)的玩家或者模型大到采樣不得不靠 GPU 加速NumPyro 非常值得嘗試。學習曲線相對陡——函數(shù)式風格、numpyro.sample的調用方式和plate與numpyro.handlers的交互邏輯都需要一些適應時間。3.2 我的選型建議按場景對號入座這里給你一個比較實際的參考場景推薦工具理由Python 新手 / 快速原型PyMC上手快、社區(qū)資料多、可視化配套完善復雜模型 / 需要最高采樣質量StanNUTS 實現(xiàn)成熟對崎嶇后驗的適應性強大數(shù)據(jù)量 / GPU 加速 / JAX 生態(tài)NumPyroJIT 向量化性能天花板最高R 用戶做統(tǒng)計建模rstan與 R 的統(tǒng)計生態(tài)銜接自然生產環(huán)境部署NumPyro / PyMC 序列化Python 棧部署路徑最順還有兩個輕量級選手值得一提。sampyl和emcee都是輕量級 MCMC 庫體量小、依賴少適合教學和小型項目。但它們功能覆蓋面有限生產級建模我還是建議回到三大件上來。我的個人經驗是不要貪多主攻一個熟悉另一個。比如你日常用 PyMC 建模遇到采樣性能瓶頸時能切到 Stan 去對比驗證或者你在 PyMC 里原型驗證沒問題之后用 NumPyro 在生產環(huán)境里重寫一遍上 GPU?？蚣苤g模型表達的思路是相通的重點是把核心概念吃透。4. 完整實操用概率編程做一個轉化率不確定性分析4.1 業(yè)務場景與數(shù)據(jù)構造我們來走一個完整的建模流程。場景設定某產品線做了一個新版本的功能產品經理想知道新版上線后購買轉化率到底比舊版高多少。收集到的數(shù)據(jù)形式很簡單舊版觀察了 2000 個用戶其中 96 人完成購買新版觀察了 2140 個用戶其中 132 人完成購買。傳統(tǒng)做法是算轉換率、做個卡方檢驗或者兩比例 z 檢驗給個 p 值。但概率編程能給出的信息豐富得多你能直接得到轉化率差異的完整后驗分布并基于概率做出決策比如“新版轉化率比舊版高的概率是 93.7%”。為了演示我們用 Python 的 PyMC 框架。核心代碼如下import pymc as pm import arviz as az import numpy as np # 觀測數(shù)據(jù) n_old, x_old 2000, 96 n_new, x_new 2140, 132 with pm.Model() as conversion_model: # 先驗舊版本轉化率 ~ Beta(1, 1) 等價于 Uniform(0, 1) p_old pm.Beta(p_old, alpha1, beta1) p_new pm.Beta(p_new, alpha1, beta1) # 似然觀測到的購買人數(shù)服從二項分布 obs_old pm.Binomial(obs_old, nn_old, pp_old, observedx_old) obs_new pm.Binomial(obs_new, nn_new, pp_new, observedx_new) # 派生量轉化率差異 diff pm.Deterministic(diff, p_new - p_old) # 執(zhí)行 MCMC 采樣 trace pm.sample(draws2000, tune2000, chains4, random_seed42)這段代碼一共做了四件事定義轉化率的先驗Beta(1,1) 表示完全沒有任何先驗傾向、定義觀測數(shù)據(jù)的似然二項分布、定義關心的派生量轉化率差異、執(zhí)行 NUTS 采樣。跑完之后trace里存的就是所有參數(shù)的 MCMC 樣本。4.2 后驗分析從樣本到業(yè)務決策采樣完成后我們把后驗樣本展示出來az.summary(trace, var_names[p_old, p_new, diff])你會得到一張表格每一行一個變量包含均值、標準差、HDI最高密度區(qū)間、有效樣本量等指標。假設輸出結果是mean sd hdi_3% hdi_97% ess_bulk r_hat p_old 0.048 0.003 0.042 0.054 6234 1.00 p_new 0.062 0.004 0.055 0.069 6157 1.00 diff 0.013 0.005 0.004 0.022 6128 1.00這份結果信息的密度遠超一個 p 值。首先p_old和p_new的 HDI 區(qū)間都是兩側遠離 0 的表明數(shù)據(jù)充分支持這兩個轉化率確實在 4%-7% 這個量級。其次diff的 94% HDI 區(qū)間是 [0.004, 0.022]不包含 0——這說明新版轉化率顯著高于舊版等等先別急著說“顯著”在貝葉斯框架里我們換一種更直白的表述。更直觀的做法是直接計算概率prob_improvement (trace.posterior[diff] 0).mean().item() print(f新版轉化率高于舊版的概率: {prob_improvement:.1%})假設輸出為 98.3%那這個信息可以直接扔給業(yè)務方我們不需要糾結“顯著性水平 0.05”或“p0.03”這種解釋易產生歧義的表述而是可以說“根據(jù)數(shù)據(jù)表現(xiàn)新版轉化率高于舊版的概率為 98.3%”。這個概率表述業(yè)務方一聽就懂決策者可以直接基于它來評估是否全量上線。4.3 先驗敏感性分析給結論上雙保險有經驗的建模者不會只看一組先驗跑出來的結果就下結論。貝葉斯方法的先驗選擇在一定程度上是主觀判斷所以要做敏感性分析換不同的先驗設置看結論會不會被顯著改變。我們把代碼改成用信息性強的先驗再跑一遍with pm.Model() as conversion_model_prior: # 更強的先驗根據(jù)歷史經驗轉化率集中在 3%-7% 之間 p_old pm.Beta(p_old, alpha10, beta170) p_new pm.Beta(p_new, alpha10, beta170) obs_old pm.Binomial(obs_old, nn_old, pp_old, observedx_old) obs_new pm.Binomial(obs_new, nn_new, pp_new, observedx_new) diff pm.Deterministic(diff, p_new - p_old) trace_prior pm.sample(draws2000, tune2000, chains4, random_seed42)對比兩次運行的diff后驗均值如果結論基本一致說明你的結論受先驗影響不大、是數(shù)據(jù)本身在驅動如果差異明顯說明當前數(shù)據(jù)量還不足以克服先驗的影響。這時候就要慎重對結論的表述要加上“基于當前先驗判斷”的限定條件。這是貝葉斯建模的好習慣也是規(guī)避“先驗挑數(shù)據(jù)”這個爭議的有效方式。4.4 層級模型的擴展不滿足于單一體轉化率上述模型是極其簡化的——它假設所有用戶同質忽略用戶來源渠道、訪問設備、時間區(qū)間等因素。一個更接近真實業(yè)務的建模方式是把這些變量吸納進模型。設想數(shù)據(jù)升級為用戶來源渠道有 A、B、C 三個每個渠道既有舊版也有新版訪問。你可以在模型里加入渠道作為分組變量構建一個層級模型# 假設 channels 是 0,1,2 編碼的渠道標簽 import pandas as pd # 這里簡化為構造一個示例 DataFrame np.random.seed(42) df pd.DataFrame({ channel: np.repeat([0, 1, 2], [800, 700, 640]), version: 0 # 占位 }) # 實際使用時數(shù)據(jù)是每個 [channel, version] 組合的訪問量和購買量 with pm.Model() as hierarchical_model: # 渠道間共享的超先驗 mu pm.Normal(mu, mu0, sigma1) sigma pm.HalfCauchy(sigma, beta0.5) # 渠道效應收縮到共同均值 z_channel pm.Normal(z_channel, mu0, sigma1, shape3) beta_channel pm.Deterministic(beta_channel, mu sigma * z_channel) # 版本效應按渠道建模時也能分層處理 p_logit mu beta_channel * ... # 根據(jù)具體編碼設計 # 似然觀測銷售數(shù)據(jù)服從伯努利/二項分布 ...層級模型的好處在于部分池化樣本量少的渠道會參考整體分布的信息避免因樣本稀疏產生極端的估計樣本量大的渠道則保持自己的特性。這在轉化率、客單價等業(yè)務指標的多維度分析中非常實用也是傳統(tǒng)頻率派方法處理起來最棘手、而概率編程最舒展的場景之一。5. 建模實操中的常見問題與排查經驗5.1 采樣不收斂最讓人抓狂的 MCMC 問題跑pm.sample()之后第一件事永遠是看r_hat和有效樣本量ess_bulk。r_hat大于 1.01 基本就意味著采樣鏈沒有收斂。常見原因和排查手段我整理成了一張速查表癥狀可能原因排查方向r_hat 恒大于 1.01鏈初始值差異太大burn-in 不夠增加 tune 長度提供更合理的初始值鏈呈“毛蟲”狀、游走緩慢后驗分布相關性高 / 參數(shù)化不佳嘗試標準化參數(shù)、非中心化參數(shù)化高維模型頻繁發(fā)散后驗幾何太崎嶇調整步長、換用pm.Metropolis對比、簡化模型結構eff 樣本量極低500自相關性過高增加chains或draws檢查模型是否有冗余參數(shù)警告出現(xiàn)很多divergencesNUTS 步長自適應遇到問題通常因先驗與數(shù)據(jù)沖突檢查先驗是否有不合理取值并調整或降低target_accept到 0.9有參數(shù)后驗分布“掛在邊界”先驗約束不夠或者數(shù)據(jù)對參數(shù)的識別力弱考慮更換先驗分布族例如在比例參數(shù)上使用 Beta 而不是 Uniform我遇到最多的還是參數(shù)化不當導致的無效采樣。一個典型的例子是方差參數(shù)用了HalfNormal(sigma10)這樣過寬的弱先驗。表面上看起來“不預設立場很客觀”但實際上過寬的先驗會在高參數(shù)空間分配大量概率質量導致采樣器在極端值之間掙扎表現(xiàn)為大量發(fā)散。改用HalfCauchy或適度窄的先驗分布之后采樣立刻恢復順暢。這里想多說一句弱先驗不等于打死不設先驗。在你對參數(shù)的合理范圍多少有數(shù)的情況下給一個符合業(yè)務常識的約束或弱信息先驗往往既不會扭曲后驗結論又能顯著改善采樣性能。這也是貝葉斯工程師的經驗之談。5.2 模型結果怪得離譜先檢查似然和數(shù)據(jù)的對應關系有次我做一個客單價預估模型后驗均值明顯偏離業(yè)務直覺。排查了半天最后發(fā)現(xiàn)是“數(shù)據(jù)對齊”的問題訓練數(shù)據(jù)里的一列應是當日訂單量結果因為某次處理腳本的 bug混入了前一天的訂單數(shù)據(jù)。這類“模型以外”的錯誤在統(tǒng)計推斷里比模型本身的錯誤更隱蔽。排查順序建議為檢查數(shù)據(jù)描述性統(tǒng)計——均值、分位數(shù)、缺失率是否合理。畫出關鍵變量在模型下的生成模擬數(shù)據(jù)posterior_predictive對比真實數(shù)據(jù)的分布。如果模擬數(shù)據(jù)和觀測數(shù)據(jù)差別巨大說明模型的一個或多個核心假設有問題。逐一簡化模型先跑一個只含主要效應、不含交互項和層級結構的模型看看能否得到大概合理的后驗再逐步增加復雜度。對連續(xù)型變量檢查數(shù)據(jù)縮放是否均勻。數(shù)值范圍跨度過大的特征會讓 NUTS 的步長自適應失靈建議對輸入數(shù)據(jù)做標準化同時可以在模型的確定性關系中使用相應的偏差。曾有很長一段時間我在做這類分析時從不做第 2 步——posterior_predictive檢查——后來一個實際案例教育了我模型對某個分組變量的層級方差估得很小把分組間差異全都“解釋”成了偶然噪聲完全不符合業(yè)務預期。當時的模型層設計沒有包含組內樣本量對先前方差收縮作用的權衡簡單來說就是層級結構設置不當。如果沒有事后預測檢查這個問題會一直潛伏到部署階段才爆發(fā)。5.3 變分推斷的結果不靠譜別急著怪算法使用變分推斷時最容易掉進的坑是q的分布族表達能力不足導致求出的近似后驗和真實后驗差異巨大但收斂過程看起來完全正常。一種驗證方法是把變分推斷結果和 MCMC 結果做對比。在數(shù)據(jù)量不大、模型不太復雜的時候這個對比成本有限卻能一針見血地發(fā)現(xiàn)問題。如果兩者差異在可接受范圍內就用變分推斷如果差異巨大檢查是不是原始模型的非光滑環(huán)節(jié)比如離散變量、硬截斷導致變分分布族嚴重失配。另外變分推斷對初始化也比較敏感多試幾組隨機種子如果結果的不確定性變動異常說明優(yōu)化過程本身可能有問題。從工程實踐角度看我的策略是變分推斷用在探索階段MCMC 用在發(fā)布階段。每次項目保質期內的最終結論都以 MCMC 驗證過的版本為準。5.4 業(yè)務落地時的不確定性傳播最后分享一個經常被忽略但實際影響很大的問題——概率編程的輸出如何進入下游系統(tǒng)。一種常見做法是只把后驗均值傳到下游這等于又把不確定性丟掉了。較好的做法是把后驗分布的完整樣本傳遞過去。具體而言如果你在做庫存優(yōu)化不是傳一個“預測需求 1300 件”而是傳一組后驗樣本“需求可能分布在 [1180, 1420]”下游的庫存優(yōu)化代碼可以對每一個樣本跑一遍最優(yōu)庫存決策得到一整套決策的分布最終選擇覆蓋某個分位的庫存量。如果你在做 A/B 實驗的后續(xù)部署判斷不是簡單用“贏了的版本全量上”而是用后驗差異的分布來確定上線節(jié)奏和回滾閾值。在做預算分配時可以把每個渠道的收益參數(shù)后驗樣本輸入優(yōu)化器這樣得到的分配方案本身就是帶不確定性的決策者可以看到最壞情況和最好情況之間的差異。很多業(yè)務場景本身并不缺少優(yōu)化方法缺少的是給優(yōu)化方法提供不確定性輸入的意識。概率編程給出了這個輸入它是從純粹統(tǒng)計推斷到實際業(yè)務決策的橋梁。6. 實戰(zhàn)進階從簡單模型走向生產級應用6.1 從單模型到貝葉斯工作流把概率編程真正用起來之后你會發(fā)現(xiàn)自己需要的不只是一個采樣函數(shù)而是一整套工作流?；谖业膶崙?zhàn)經驗一個可復用的貝葉斯建模流程應該是這樣的問題形式化把業(yè)務問題翻譯成一個可以計算概率的可量化問題。明確哪些是觀測變量、哪些是待推斷參數(shù)、哪些是派生業(yè)務指標。數(shù)據(jù)審查做描述性統(tǒng)計和可視化識別數(shù)據(jù)缺失、異常值和量綱問題。先驗設立基于業(yè)務常識設定合理的弱信息先驗。向業(yè)務方了解“過去的歷史水平”讓它成為先驗均值的重要參考。模型搭建與迭代從簡單模型開始驗證核心結構正確后再加復雜度。每次增加結構后都用后驗預測檢查來檢驗改進效果。推斷與診斷執(zhí)行 MCMC或變分推斷檢查r_hat、ess_bulk、發(fā)散情況。必要時做參數(shù)化調整。敏感性分析改變先驗設定、修改數(shù)據(jù)子集做重跑確認結論的穩(wěn)定性。產出與溝通把后驗分布翻譯成業(yè)務語言輸出“差異為 X 的概率是多少”、“落在哪個區(qū)間的把握有多大”這類可操作的結論。部署與監(jiān)控把模型固化為可重復執(zhí)行的推斷管道對新數(shù)據(jù)周期性地更新后驗并監(jiān)控分布漂移。這八步聽起來多但每一步的時間消耗并不大。熟練之后一個中等復雜度的模型從開始到出結論基本在一兩天內就能完成。6.2 讓模型可復現(xiàn)的幾個工程習慣概率編程的模型天然具有隨機性因此可復現(xiàn)性是工程化繞不開的問題。沒有固定隨機種子的模型運行每次結果都有細微差異這在調試時會讓人分不清改動模型前的效果差異是真實的還是隨機波動造成的。我現(xiàn)在在項目里嚴格遵守以下三條規(guī)則固定隨機種子所有采樣、變分推斷調用都顯式傳入random_seed并使它在每次運行中保持一致。記錄模型版本和依賴環(huán)境把模型代碼、依賴庫版本號、數(shù)據(jù)文件哈希寫進自動生成的模型元數(shù)據(jù)文件。后期排查問題這一條能救命。用緩存機制保存推斷結果同一個模型、同一份數(shù)據(jù)跑過一次后把 trace 序列化為文件。后續(xù)所有人的分析都基于這份緩存結果保證溝通討論的前后一致。這些習慣沒有任何高深之處但它們直接決定了概率編程能否從小項目走向正式工作流。7. 一些關于概率編程的實話做概率編程這些年最大的感受是這項技術真正改變人的地方不在采樣算法的效率而在思維方式的轉換。當你開始用分布去思考所有問題你會意識到以前那些點估計的結論有多少是被“平均”掩埋掉的。但也要說句公道話概率編程不是萬能的。它的計算成本遠高于傳統(tǒng)統(tǒng)計方法在超大樣本場景下需要投入相當?shù)臅r間去做工程優(yōu)化。如果業(yè)務場景的數(shù)據(jù)量足夠大、模型相對簡單、只需要一個穩(wěn)定輸出傳統(tǒng)方法可能更實惠。概率編程適合的是對不確定性本身有決策價值、數(shù)據(jù)規(guī)?？煽?、模型復雜度中上的場景——用一句直白的話說當你發(fā)現(xiàn)“上線后預測總偏差大、又不知道偏差從哪來”的時候概率編程大概能給你答案。另外概率編程也不該被你當成一個“萬能答案機”。它給出的后驗分布依然建立在模型假設之上。假設錯了后驗再漂亮也不可靠。這也是為什么做概率編程的人一半時間在建模另一半時間在驗證假設、檢查模型與數(shù)據(jù)的匹配度。這些約束和代價不是說概率編程不好而是說你有必要想清楚我到底需要的是不確定性信息還是僅僅需要結論的確定性去完成眼前這一步。想明白了再上手這個工具給你的回報會遠高于你的預期。