——皮埃爾—西蒙·拉普拉斯
本文我們將學(xué)習(xí)貝葉斯統(tǒng)計中的核心概念以及一些用于貝葉斯分析的基本工具。大部分內(nèi)容都是一些理論介紹，其中會涉及一些Python代碼，盡管本文內(nèi)容有點偏理論，可能會讓習(xí)慣代碼的你感到有點不安，不過這會讓你在后面應(yīng)用貝葉斯統(tǒng)計方法解決問題時容易一些。
本文包含以下主題：
統(tǒng)計模型；
概率及不確定性；
貝葉斯理論及統(tǒng)計推斷；
單參數(shù)推斷以及經(jīng)典的拋硬幣問題；
如何選擇先驗；
如何報告貝葉斯分析結(jié)果；
安裝所有相關(guān)的Python庫。
統(tǒng)計學(xué)主要是收集、組織、分析并解釋數(shù)據(jù)，因此，統(tǒng)計學(xué)的基礎(chǔ)知識對數(shù)據(jù)分析來說至關(guān)重要。分析數(shù)據(jù)時一個非常有用的技巧是知道如何運用某種編程語言（如Python）編寫代碼。真實世界里充斥著復(fù)雜而雜亂的數(shù)據(jù)，因此對數(shù)據(jù)做一些預(yù)處理操作必不可少。即便你的數(shù)據(jù)已經(jīng)是整理好的了，掌握一定的編程技巧仍然會給你帶來很大幫助，因為如今的貝葉斯統(tǒng)計絕大多數(shù)都是計算統(tǒng)計學(xué)。
數(shù)據(jù)是統(tǒng)計學(xué)最基本的組成部分。數(shù)據(jù)的來源多，比如實驗、計算機模擬、調(diào)查以及觀測等。假如我們是數(shù)據(jù)生成或收集人員，首先要考慮的是要解決什么樣的問題以及打算采用什么方法，然后再去著手準(zhǔn)備數(shù)據(jù)。事實上，統(tǒng)計學(xué)有一個叫做實驗設(shè)計的分支專門研究如何獲取數(shù)據(jù)。在這個數(shù)據(jù)泛濫的年代，我們有時候會忘了獲取數(shù)據(jù)并非總是很便宜。例如，盡管大型強子碰撞加速裝置一天能產(chǎn)生上百TB的數(shù)據(jù)，但其建造卻要花費數(shù)年的人力和智力。本文已經(jīng)獲取了數(shù)據(jù)并且數(shù)據(jù)是整理好的（這在現(xiàn)實中通常很少見），以便關(guān)注到本文的主題上來。
假設(shè)我們已經(jīng)有了數(shù)據(jù)集，通常的做法是先對其探索并可視化，這樣我們就能對手頭的數(shù)據(jù)有個直觀的認(rèn)識?？梢酝ㄟ^如下兩步完成所謂的探索式數(shù)據(jù)分析過程：
描述性統(tǒng)計；
數(shù)據(jù)可視化。
其中，描述性統(tǒng)計是指如何用一些指標(biāo)或統(tǒng)計值來定量地總結(jié)或刻畫數(shù)據(jù)，例如你已經(jīng)知道了如何用均值、眾數(shù)、標(biāo)準(zhǔn)差、四分位差等指標(biāo)來描述數(shù)據(jù)。數(shù) 據(jù)可視化是指用生動形象的方式表述數(shù)據(jù)，你大概對直方圖、散點圖等表現(xiàn)形式比較熟悉。乍看起來，探索式數(shù)據(jù)分析似乎是在復(fù)雜分析之前的一些準(zhǔn)備工作，或者是作為一些復(fù)雜分析方法的替代品，不過探索式數(shù)據(jù)分析在理解、解釋、檢查、總結(jié)及交流貝葉斯分析結(jié)果等過程中依然有用。
有時候，畫畫圖、對數(shù)據(jù)做些簡單的計算（比如求均值）就夠了。另外一些時候，我們希望從數(shù)據(jù)中挖掘出一些更一般性的結(jié)論。我們可能希望了解數(shù)據(jù)是怎么生成的，也可能是想對未來還未觀測到的數(shù)據(jù)做出預(yù)測，又或者是希望從多個對觀測值的解釋中找出最合理的一個，這些正是統(tǒng)計推斷所做的事情。模型分為許多種，統(tǒng)計推斷依賴的是概率模型，許多科學(xué)研究（以及我們對真實世界的認(rèn)識）也都是基于模型的， 大腦不過是對現(xiàn)實進(jìn)行建模的一臺機器，可以觀看相關(guān)的TED演講了解大腦是如何對現(xiàn)實進(jìn)行建模的，網(wǎng)址為http://www.tedxriodelaplata.org/videos/m%C3%A1quina-construye-realidad。
什么是模型？模型是對給定系統(tǒng)或過程的一種簡化描述。這些描述只關(guān)注系統(tǒng)中某些重要的部分，因此，大多數(shù)模型的目的并不是解釋整個系統(tǒng)。此外，假如我們有兩個模型能用來解釋同一份數(shù)據(jù)并且效果差不多，其中一個簡單點，另一個復(fù)雜一些，通常我們傾向于更簡單的模型，這稱作奧卡姆剃刀，我們會在第6章模型比較部分討論貝葉斯分析與其之間的聯(lián)系。
不管你打算構(gòu)建哪種模型，模型構(gòu)建都遵循一些相似的基本準(zhǔn)則，我們把貝葉斯模型的構(gòu)建過程總結(jié)為如下3步。
（1）給定一些數(shù)據(jù)以及這些數(shù)據(jù)是如何生成的假設(shè)，然后構(gòu)建模型。通常，這里的模型都是一些很粗略的近似，不過大多時候也夠用了。
（2）利用貝葉斯理論將數(shù)據(jù)和模型結(jié)合起來，根據(jù)數(shù)據(jù)和假設(shè)推導(dǎo)出邏輯結(jié)論，我們稱之為經(jīng)數(shù)據(jù)擬合后的模型。
（3）根據(jù)多種標(biāo)準(zhǔn)，包括真實數(shù)據(jù)和對研究問題的先驗知識，判斷模型擬合得是否合理。
通常，我們會發(fā)現(xiàn)實際的建模過程并非嚴(yán)格按照該順序進(jìn)行的，有時候我們有可能跳到其中任何一步，原因可能是編寫的程序出錯了，也可能是找到了某種改進(jìn)模型的方式，又或者是我們需要增加更多的數(shù)據(jù)。
貝葉斯模型是基于概率構(gòu)建的，因此也稱作概率模型。為什么基于概率呢？因為概率這個數(shù)學(xué)工具能夠很好地描述數(shù)據(jù)中的不確定性，接下來我們將對其進(jìn)行深入了解。
盡管概率論是數(shù)學(xué)中一個相當(dāng)成熟和完善的分支，但關(guān)于概率的詮釋仍然有不止一種。對于貝葉斯派而言，概率是對某一命題不確定性的衡量。假設(shè)我們對硬幣一無所知，同時沒有與拋硬幣相關(guān)的任何數(shù)據(jù)，那么可以認(rèn)為正面朝上的概率介于0到1之間，也就是說，在缺少信息的情況下，所有情況都是有可能發(fā)生的，此時不確定性也最大。假設(shè)現(xiàn)在我們知道硬幣是公平的，那么我們可以認(rèn)為正面朝上的概率是0.5或者是0.5附近的某個值（假如硬幣不是絕對公平的話），如果此時收集數(shù)據(jù)，我們可以根據(jù)觀測值進(jìn)一步更新前面的先驗假設(shè)，從而降低對該硬幣偏差的不確定性。按照這種定義，提出以下問題都是自然而且合理的：火星上有多大可能存在生命？電子的質(zhì)量是9.1×10?31kg的概率是多大？1816年7月9號是晴天的概率是多少？值得注意的是，類似火星上是否有生命這種問題的答案是二值化的，但我們關(guān)心的是，基于現(xiàn)有數(shù)據(jù)以及我們對火星物理條件和生物條件的了解，火星上存在生命的概率有多大。該命題取決于我們當(dāng)前所掌握的信息而非客觀的自然屬性。我們使用概率是因為我們對事件不確定，而不代表事件本身是不確定的。由于這種概率的定義取決于我們的認(rèn)知水平，有時也被稱為概率的主觀定義，這也就解釋了為什么貝葉斯派總被稱作主觀統(tǒng)計。然而，這種定義并不是說所有命題都是同等有意義的，僅是承認(rèn)我們對世界的理解是基于現(xiàn)有的數(shù)據(jù)和模型，是不完備的。不基于模型或理論去理解世界是不可能的。即使我們能脫離社會預(yù)設(shè)，最終也會受到生物上的限制：受進(jìn)化過程影響，我們的大腦已經(jīng)與世界上的模型建立了聯(lián)系。我們注定要以人類的方式去思考，永遠(yuǎn)不可能像蝙蝠或其他動物那樣思考。而且，宇宙是不確定的，總的來說我們能做的就是對其進(jìn)行概率性描述。需要注意的是，不管世界的本原是確定的還是隨機的，我們都將概率作為衡量不確定性的工具。
邏輯是指如何有效地推論，在亞里士多德學(xué)派或者經(jīng)典的邏輯學(xué)中，一個命題只能是對或者錯，而在貝葉斯學(xué)派定義的概率中，確定性只是概率上的一種特殊情況：一個正確命題的概率值為1，而一個錯誤命題的概率值為0。只有在我們擁有充分的數(shù)據(jù)表明火星上存在生物生長和繁殖時，我們才認(rèn)為“火星上存在生命”這一命題為真的概率值為1。不過，需要注意的是，認(rèn)定一個命題的概率值為0通常是比較困難的，這是因為火星上可能存在某些地方還沒被探測到，或者是我們的實驗有問題，又或者是某些其他原因?qū)е挛覀冨e誤地認(rèn)為火星上沒有生命而實際上有。與此相關(guān)的是克倫威爾準(zhǔn)則（Cromwell’s Rule），其含義是指在對邏輯上正確或錯誤的命題賦予概率值時，應(yīng)當(dāng)避免使用0或者1。有意思的是，考克斯（Cox）在數(shù)學(xué)上證明了如果想在邏輯推理中加入不確定性，我們就必須使用概率論的知識，由此來說，貝葉斯定理可以視為概率論邏輯上的結(jié)果。因此，從另一個角度來看，貝葉斯統(tǒng)計是對邏輯學(xué)處理不確定性問題的一種擴充，當(dāng)然這里毫無對主觀推理的輕蔑?，F(xiàn)在我們已經(jīng)熟悉了貝葉斯學(xué)派對概率的理解，接下來就了解下概率相關(guān)的數(shù)學(xué)特性吧。如果你想深入學(xué)習(xí)概率論，可以參考閱讀Joseph K Blitzstein和Jessica Hwang寫的《概率導(dǎo)論》（Introduction to Probability）。
概率值介于[0,1]之間（包括0和1），其計算遵循一些法則，其中之一是乘法法則：
上式中，A和B同時發(fā)生的概率值等于B發(fā)生的概率值乘以在B發(fā)生的條件下A也發(fā)生的概率值，其中，表示A和B的聯(lián)合概率， 表示條件概率，二者的現(xiàn)實意義是不同的，例如，路面是濕的概率跟下雨時候路面是濕的概率是不同的。條件概率可能比原來的概率高，也可能低。如果B并不能提供任何關(guān)于A的信息，那么，
條件概率是統(tǒng)計學(xué)中的一個核心概念，接下來我們將看到，理解條件概率對于理解貝葉斯理論至關(guān)重要。這里我們換個角度來看條件概率，假如我們把前面的公式調(diào)整下順序，就可以得到下面的公式：
需要注意的是，我們不對0概率事件計算條件概率，因此，分母的取值范圍是(0,1)，從而可以看出條件概率大于或等于聯(lián)合概率。為什么要除以p(B)呢？因為在已經(jīng)知道事件B發(fā)生的條件下，我們考慮的可能性空間就縮小到了事件B發(fā)生的范圍內(nèi)，然后將該范圍內(nèi)A發(fā)生的可能性除以B發(fā)生的可能性便得到了條件概率 。需要強調(diào)的是，所有的概率本質(zhì)上都是條件概率，世間并沒有絕對的概率。不管我們是否留意，在談到概率時總是存在這樣或那樣的模型、假設(shè)或條件。比如，當(dāng)我們討論明天下雨的概率時，在地球上、火星上甚至宇宙中其他地方明天下雨的概率是不同的。同樣，硬幣正面朝上的概率取決于我們對硬幣有偏性的假設(shè)。在理解概率的含義之后，接下來我們討論另一個話題：概率分布。
概率分布是數(shù)學(xué)中的一個概念，用來描述不同事件發(fā)生的可能性，通常這些事件限定在一個集合內(nèi)，代表了所有可能發(fā)生的事件。在統(tǒng)計學(xué)里可以這么理解：數(shù)據(jù)是從某種參數(shù)未知的概率分布中生成的。由于并不知道具體的參數(shù)，我們只能借用貝葉斯定理從僅有的數(shù)據(jù)中反推參數(shù)。概率分布是構(gòu)建貝葉斯模型的基礎(chǔ)，不同分布組合在一起之后可以得到一些很有用的復(fù)雜模型。
本文會介紹一些概率分布，在第一次介紹某個概率分布時，我們會先花點時間理解它。最常見的一種概率分布是高斯分布，又叫正態(tài)分布，其數(shù)學(xué)公式描述如下：
上式中，μ和σ是高斯分布的兩個參數(shù)。第1個參數(shù)μ是該分布的均值（同時也是中位數(shù)和眾數(shù)），其取值范圍是任意實數(shù)，即
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import seaborn as sns
mu_params = [-1, 0, 1]
sd_params = [0.5, 1, 1.5]
x = np.linspace(-7, 7, 100)
f, ax = plt.subplots(len(mu_params), len(sd_params), sharex=True, sharey=True)
for i in range(3):
for j in range(3):
mu = mu_params[i]
sd = sd_params[j]
y = stats.norm(mu, sd).pdf(x)
ax[i,j].plot(x, y)
ax[i,j].plot(0, 0,
label='$\\mu$ = {:3.2f}\n$\\sigma$ = {:3.2f}'.format (mu, sd), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$x$', fontsize=16)
ax[1,0].set_ylabel('$pdf(x)$', fontsize=16)
plt.tight_layout()
上面代碼的輸出結(jié)果如下：
由概率分布生成的變量（例如x）稱作隨機變量，當(dāng)然這并不是說該變量可以取任意值，相反，我們觀測到該變量的數(shù)值受到概率分布的約束，而其隨機性源于我們只知道變量的分布卻無法準(zhǔn)確預(yù)測該變量的值。通常，如果一個隨機變量服從在參數(shù)μ和σ下的高斯分布，我們可以這樣表示該變量：
其中，符號～讀作服從于某種分布。
隨機變量分為兩種：連續(xù)變量和離散變量。連續(xù)隨機變量可以從某個區(qū)間內(nèi)取任意值（我們可以用Python中的浮點型數(shù)據(jù)來表示），而離散隨機變量只能取某些特定的值（我們可以用Python中的整型數(shù)據(jù)來表示）。
許多模型都假設(shè)，如果對服從于同一個分布的多個隨機變量進(jìn)行連續(xù)采樣，那么各個變量的采樣值之間相互獨立，我們稱這些隨機變量是獨立同分布的。用數(shù)學(xué)語言描述就是，如果兩個隨機變量x和y對于所有可能的取值都滿足
時間序列是不滿足獨立同分布的一個典型例子。在時間序列中，需要對時間維度的變量多加留心。下面的例子是從http://cdiac.esd.ornl.gov中獲取的數(shù)據(jù)。這份數(shù)據(jù)記錄了從1959年到1997年大氣中二氧化碳的含量。我們用以下代碼將它寫出來：
data = np.genfromtxt('mauna_loa_CO2.csv', delimiter=',')
plt.plot(data[:,0], data[:,1])
plt.xlabel('$year$', fontsize=16)
plt.ylabel('$CO_2 (ppmv)$', fontsize=16)
圖中每個點表示每個月空氣中二氧化碳含量的測量值，可以看到測量值是與時間相關(guān)的。我們可以觀察到兩個趨勢：一個是季節(jié)性的波動趨勢（這與植物周期性生長和衰敗有關(guān)）；另一個是二氧化碳含量整體性的上升趨勢。
到目前為止，我們已經(jīng)學(xué)習(xí)了一些統(tǒng)計學(xué)中的基本概念和詞匯，接下來讓我們首先看看神奇的貝葉斯定理：
看起來稀松平常，似乎跟小學(xué)課本里的公式差不多，不過這就是關(guān)于貝葉斯統(tǒng)計你所需要掌握的全部。首先看看貝葉斯定理是怎么來的，這對我們理解它會很有幫助。事實上，我們已經(jīng)掌握了如何推導(dǎo)它所需要的全部概率論知識。
根據(jù)前面提到的概率論中的乘法準(zhǔn)則，我們有以下式子：
上式還可以寫成如下形式：
由于以上式子的左邊相等，于是可以得到：
對上式調(diào)整下順序，便得到了貝葉斯定理：
現(xiàn)在，讓我們看看這個式子的含義及其重要性。首先，上式表明
在前面的式子中，如果我們將H理解為假設(shè)，D理解為數(shù)據(jù)，那么貝葉斯定理告訴我們的就是，在給定數(shù)據(jù)的條件下如何計算假設(shè)成立的概率。不過，如何把假設(shè)融入貝葉斯定理中去呢？答案是概率分布。換句話說，H是一種狹義上的假設(shè)，我們所做的實際上是尋找模型的參數(shù)（更準(zhǔn)確地說是參數(shù)的分布）。因此，與其稱H為假設(shè)，不如稱之為模型，這樣能避免歧義。
貝葉斯定理非常重要，后面會反復(fù)用到，這里我們先熟悉下其各個部分的名稱：
p(H )：先驗；
p(D | H )：似然；
p(H | D)：后驗；
p(D)：證據(jù)。
先驗分布反映的是在觀測到數(shù)據(jù)之前我們對參數(shù)的了解，如果我們對參數(shù)一無所知（就跟《權(quán)力的游戲》中的雪諾一樣），那么可以用一個不包含太多信息的均勻分布來表示。由于引入了先驗，有些人會認(rèn)為貝葉斯統(tǒng)計是偏主觀的，然而，這些先驗不過是構(gòu)建模型時的一些假設(shè)罷了，其主觀性跟似然差不多。
似然是指如何在實驗分析中引入觀測數(shù)據(jù)，反映的是在給定參數(shù)下得到某組觀測數(shù)據(jù)的可信度。
后驗分布是貝葉斯分析的結(jié)果，反映的是在給定數(shù)據(jù)和模型的條件下我們對問題的全部認(rèn)知。需要注意，后驗指的是我們模型中參數(shù)的概率分布而不是某個值，該分布正比于先驗乘以似然。有這么個笑話：貝葉斯學(xué)派就像是這樣一類人，心里隱約期待著一匹馬，偶然間瞥見了一頭驢，結(jié)果堅信他看到的是一頭騾子。當(dāng)然了，如果要刻意糾正這個笑話的話，在先驗和似然都比較含糊的情況下，我們會得到一個（模糊的）“騾子”后驗。不過，這個笑話也講出了這樣一個道理，后驗其實是對先驗和似然的某種折中。從概念上講，后驗可以看做是在觀測到數(shù)據(jù)之后對先驗的更新。事實上，一次分析中的后驗，在收集到新的數(shù)據(jù)之后，也可以看做是下一次分析中的先驗。這使得貝葉斯分析特別適合于序列化的數(shù)據(jù)分析，比如通過實時處理來自氣象站和衛(wèi)星的數(shù)據(jù)從而提前預(yù)警災(zāi)害，更詳細(xì)的內(nèi)容可以閱讀在線機器學(xué)習(xí)方面的算法。
最后一個概念是證據(jù)，也稱作邊緣似然。正式地講，證據(jù)是在模型的參數(shù)取遍所有可能值的條件下得到指定觀測值的概率的平均。不過，本文的大部分內(nèi)容并不關(guān)心這個概念，我們可以簡單地把它當(dāng)作歸一化系數(shù)。我們只關(guān)心參數(shù)的相對值而非絕對值。把證據(jù)這一項忽略掉之后，貝葉斯定理可以表示成如下正比例形式：
前面，我們學(xué)習(xí)了幾個重要概念，其中有兩個是貝葉斯統(tǒng)計的核心概念，這里我們用一句話再重新強調(diào)下：概率是用來衡量參數(shù)不確定性的，貝葉斯定理就是用來在觀測到新的數(shù)據(jù)時正確更新這些概率以期降低我們的不確定性。
現(xiàn)在我們已經(jīng)知道什么是貝葉斯統(tǒng)計了，接下來就從一個簡單的例子入手，通過推斷單個未知參數(shù)來學(xué)習(xí)如何進(jìn)行貝葉斯統(tǒng)計。
拋硬幣是統(tǒng)計學(xué)中的一個經(jīng)典問題，其描述如下：我們隨機拋一枚硬幣，重復(fù)一定次數(shù)，記錄正面朝上和反面朝上的次數(shù)，根據(jù)這些數(shù)據(jù)，我們需要回答諸如這枚硬幣是否公平，以及更進(jìn)一步這枚硬幣有多不公平等問題。拋硬幣是一個學(xué)習(xí)貝葉斯統(tǒng)計非常好的例子，一方面是因為幾乎人人都熟悉拋硬幣這一過程，另一方面是因為這個模型很簡單，我們可以很容易計算并解決這個問題。此外，許多真實問題都包含兩個互斥的結(jié)果，例如0或者1、正或者負(fù)、奇數(shù)或者偶數(shù)、垃圾郵件或者正常郵件、安全或者不安全、健康或者不健康等。因此，即便我們討論的是硬幣，該模型也同樣適用于前面這些問題。
為了估計硬幣的偏差，或者更廣泛地說，想要用貝葉斯學(xué)派理論解決問題，我們需要數(shù)據(jù)和一個概率模型。對于拋硬幣這個問題，假設(shè)我們已試驗了一定次數(shù)并且記錄了正面朝上的次數(shù)，也就是說數(shù)據(jù)部分已經(jīng)準(zhǔn)備好了，剩下的就是模型部分了。考慮到這是第一個模型，我們會列出所有必要的數(shù)學(xué)公式，并且一步一步推導(dǎo)。
首先，我們要抽象出偏差的概念。我們稱，如果一枚硬幣總是正面朝上，那么它的偏差就是1，反之，如果總是反面朝上，那么它的偏差就是0，如果正面朝上和反面朝上的次數(shù)各占一半，那么它的偏差就是0.5。這里用參數(shù)θ來表示偏差，用y表示N次拋硬幣實驗中正面朝上的次數(shù)。根據(jù)貝葉斯定理，我們有如下公式：
這里需要指定我們將要使用的先驗
假設(shè)多次拋硬幣的結(jié)果相互之間沒有影響，也就是說每次拋硬幣都是相互獨立的，同時還假設(shè)結(jié)果只有兩種可能：正面朝上或者反面朝上。但愿你能認(rèn)同我們對這個問題做出的合理假設(shè)?；谶@些假設(shè)，一個不錯的似然候選是二項分布：
這是一個離散分布，表示N次拋硬幣實驗中y次正面朝上的概率（或者更通俗地描述是，N次實驗中，y次成功的概率）。下面的代碼生成了9個二項分布，每個子圖中的標(biāo)簽顯示了對應(yīng)的參數(shù)：
n_params = [1, 2, 4]
p_params = [0.25, 0.5, 0.75]
x = np.arange(0, max(n_params)+1)
f, ax = plt.subplots(len(n_params), len(p_params), sharex=True, sharey=True)
for i in range(3):
for j in range(3):
n = n_params[i]
p = p_params[j]
y = stats.binom(n=n, p=p).pmf(x)
ax[i,j].vlines(x, 0, y, colors='b', lw=5)
ax[i,j].set_ylim(0, 1)
ax[i,j].plot(0, 0, label='n = {:3.2f}\np =
{:3.2f}'.format(n, p), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$\\theta$', fontsize=14)
ax[1,0].set_ylabel('$p(y|\\theta)$', fontsize=14)
ax[0,0].set_xticks(x)
二項分布是似然的一個合理選擇，直觀上講，θ可以看作拋一次硬幣時正面朝上的可能性，并且該過程發(fā)生了y次。類似地，我們可以把“1?θ”看作拋一次硬幣時反面朝上的概率，并且該過程發(fā)生了“N?y”次。
假如我們知道了θ，那么就可以從二項分布得出硬幣正面朝上的分布。如果我們不知道θ，也別灰心，在貝葉斯統(tǒng)計中，當(dāng)我們不知道某個參數(shù)的時候，就對其賦予一個先驗。接下來繼續(xù)選擇先驗。
這里我們選用貝葉斯統(tǒng)計中最常見的beta分布，作為先驗，其數(shù)學(xué)形式如下：
仔細(xì)觀察上面的式子可以看出，除了Γ部分之外，beta分布和二項分布看起來很像。Γ是希臘字母中大寫的伽馬，用來表示伽馬函數(shù)?，F(xiàn)在我們只需要知道，用分?jǐn)?shù)表示的第一項是一個正則化常量，用來保證該分布的積分為1，此外，
params = [0.5, 1, 2, 3]
x = np.linspace(0, 1, 100)
f, ax = plt.subplots(len(params), len(params), sharex=True, sharey=True)
for i in range(4):
for j in range(4):
a = params[i]
b = params[j]
y = stats.beta(a, b).pdf(x)
ax[i,j].plot(x, y)
ax[i,j].plot(0, 0, label='$\\alpha$ = {:3.2f}\n$\\beta$ = {:3.2f}'.format(a, b), alpha=0)
ax[i,j].legend(fontsize=12)
ax[3,0].set_xlabel('$\\theta$', fontsize=14)
ax[0,0].set_ylabel('$p(\\theta)$', fontsize=14)
plt.savefig('B04958_01_04.png', dpi=300, figsize=(5.5, 5.5))
為什么要在模型中使用beta分布呢？在拋硬幣以及一些其他問題中使用beta分布的原因之一是：beta分布的范圍限制在0到1之間，這跟我們的參數(shù)
首先回憶一下貝葉斯定理：后驗正比于似然乘以先驗。
對于我們的問題而言，需要將二項分布乘以beta分布：
現(xiàn)在，對上式進(jìn)行簡化。針對我們的實際問題，可以先把與θ不相關(guān)的項去掉而不影響結(jié)果，于是得到下式：
重新整理之后得到：
可以看出，上式和beta分布的形式很像（除了歸一化部分），其對應(yīng)的參數(shù)分別為
計算后驗并畫圖
現(xiàn)在已經(jīng)有了后驗的表達(dá)式，我們可以用Python對其計算并畫出結(jié)果。下面的代碼中，其實只有一行是用來計算后驗結(jié)果的，其余的代碼都是用來畫圖的：
theta_real = 0.35
trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
beta_params = [(1, 1), (0.5, 0.5), (20, 20)]
dist = stats.beta
x = np.linspace(0, 1, 100)
for idx, N in enumerate(trials):
if idx == 0:
plt.subplot(4,3, 2)
else:
plt.subplot(4,3, idx+3)
y = data[idx]
for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')):
p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)
plt.plot(x, p_theta_given_y, c)
plt.fill_between(x, 0, p_theta_given_y, color=c, alpha=0.6)
plt.axvline(theta_real, ymax=0.3, color='k')
plt.plot(0, 0, label='{:d} experiments\n{:d} heads'.format(N, y), alpha=0)
plt.xlim(0,1)
plt.ylim(0,12)
plt.xlabel(r'$\theta$')
plt.legend()
plt.gca().axes.get_yaxis().set_visible(False)
plt.tight_layout()
plt.savefig('B04958_01_05.png', dpi=300, figsize=(5.5, 5.5))
在上圖的第一行中，實驗的次數(shù)為0，因此第一個圖中的曲線描繪的是先驗分布，其中有3條曲線，每條曲線分別表示一種先驗。
藍(lán)色的線是一個均勻分布先驗，其含義是：偏差的所有可能取值都是等概率的。
紅色的線與均勻分布有點類似，對拋硬幣這個例子而言可以理解為：偏差等于0或者1的概率要比其他值更大一些。
最后一條綠色的線集中在中間值0.5附近，該分布反映了通常硬幣正面朝上和反面朝上的概率大致是差不多的。我們也可以稱，該先驗與大多數(shù)硬幣都是公平的這一信念是兼容的。“兼容”這個詞在貝葉斯相關(guān)的討論中會經(jīng)常用到，特別是在提及受到數(shù)據(jù)啟發(fā)的模型時。
剩余的子圖描繪了后續(xù)實驗的后驗分布
貝葉斯分析的結(jié)果是后驗分布而不是某個值，該分布描述了根據(jù)給定數(shù)據(jù)和模型得到的不同數(shù)值的可能性。
后驗最可能的值是根據(jù)后驗分布的形態(tài)決定的（也就是后驗分布的峰值）。
后驗分布的離散程度與我們對參數(shù)的不確定性相關(guān)；分布越離散，不確定性越大。
盡管1/2=4/8=0.5，但從圖中可以看出，前者的不確定性要比后者大。這是因為我們有了更多的數(shù)據(jù)來支撐我們的推斷，該直覺也同時反映在了后驗分布上。
在給定足夠多的數(shù)據(jù)時，兩個或多個不同先驗的貝葉斯模型會趨近于收斂到相同的結(jié)果。在極限情況下，如果有無限多的數(shù)據(jù)，不論我們使用的是怎樣的先驗，最終都會得到相同的后驗。注意這里說的無限多是指某種程度而非某個具體的數(shù)量，也就是說，從實際的角度來講，某些情況下無限多的數(shù)據(jù)可以通過比較少量的數(shù)據(jù)近似。
不同后驗收斂到相同分布的速度取決于數(shù)據(jù)和模型。從前面的圖中可以看出，藍(lán)色和紅色的后驗在經(jīng)過8次實驗之后就很難看出區(qū)別了，而紅色的曲線則一直到150次實驗之后才與另外兩個后驗看起來比較接近。
有一點從前面的圖中不太容易看出來，如果我們一步一步地更新后驗，最后得到的結(jié)果跟一次性計算得到的結(jié)果是一樣的。換句話說，我們可以對后驗分150次計算，每次增加一個新的觀測數(shù)據(jù)并將得到的后驗作為下一次計算的先驗，也可以在得到150次拋硬幣的結(jié)果之后一次性計算出后驗，而這兩種計算方式得到的結(jié)果是完全一樣的。這個特點非常有意義，在許多數(shù)據(jù)分析的問題中，每當(dāng)我們得到新的數(shù)據(jù)時可以更新估計值。
從前面的例子可以看出，先驗對分析的結(jié)果會有影響。一些貝葉斯分析的新手（以及一些詆毀該方法的人）會對如何選擇先驗感到茫然，因為他們不希望先驗起到?jīng)Q定性作用，而是更希望數(shù)據(jù)本身替自己說話！有這樣的想法很正常，不過我們得牢記，數(shù)據(jù)并不會真的“說話”，只有在模型中才會有意義，包括數(shù)學(xué)上的和腦海中的模型。面對同一主題下的同一份數(shù)據(jù)，不同人會有不同的看法，這類例子在科學(xué)史上有許多，查看以下鏈接可以了解《紐約時報》最近一次實驗的例子：http://www.nytimes.com/interactive/2016/09/20/upshot/the-error-the-polling-world-rarely-talks-about.html?_r=0。
有些人青睞于使用沒有信息量的先驗（也稱作均勻的、含糊的或者發(fā)散的先驗），這類先驗對分析過程的影響最小。本文將遵循Gelman、McElreath和Kruschke 3人的建議[1]，更傾向于使用帶有較弱信息量的先驗。在許多問題中，我們對參數(shù)可以取的值一般都會有些了解，比如，參數(shù)只能是正數(shù)，或者知道參數(shù)近似的取值范圍，又或者是希望該值接近0或大于/小于某值。這種情況下，我們可以給模型加入一些微弱的先驗信息而不必?fù)?dān)心該先驗會掩蓋數(shù)據(jù)本身的信息。由于這類先驗會讓后驗近似位于某一合理的邊界內(nèi)，因此也被稱作正則化先驗。
當(dāng)然，使用帶有較多信息量的強先驗也是可行的。視具體的問題不同，有可能很容易或者很難找到這類先驗，例如在我工作的領(lǐng)域（結(jié)構(gòu)生物信息學(xué)），人們會盡可能地利用先驗信息，通過貝葉斯或者非貝葉斯的方式來了解和預(yù)測蛋白質(zhì)的結(jié)構(gòu)。這樣做是合理的，原因是我們在數(shù)十年間已經(jīng)從上千次精心設(shè)計的實驗中收集了數(shù)據(jù)，因而有大量可信的先驗信息可供使用。如果你有可信的先驗信息，完全沒有理由不去使用。試想一下，如果一個汽車工程師每次設(shè)計新車的時候，他都要重新發(fā)明內(nèi)燃機、輪子乃至整個汽車，這顯然不是正確的方式。
現(xiàn)在我們知道了先驗有許多種，不過這并不能緩解我們選擇先驗時的焦慮?；蛟S，最好是沒有先驗，這樣事情就簡單了。不過，不論是否基于貝葉斯，模型都在某種程度上擁有先驗，即使這里的先驗并沒有明確表示出來。事實上，許多頻率統(tǒng)計學(xué)方面的結(jié)果可以看做是貝葉斯模型在一定條件下的特例，比如均勻先驗。讓我們再仔細(xì)看看前面那幅圖，可以看到藍(lán)色后驗分布的峰值與頻率學(xué)分析中θ的期望值是一致的：
注意，這里
在特定分析任務(wù)中，如果我們對某個先驗或者似然不確定，可以自由使用多個先驗或者似然進(jìn)行嘗試。模型構(gòu)建過程中的一個環(huán)節(jié)就是質(zhì)疑假設(shè)，而先驗就是質(zhì)疑的對象之一。不同的假設(shè)會得到不同的模型，根據(jù)數(shù)據(jù)和與問題相關(guān)的領(lǐng)域知識，我們可以對這些模型進(jìn)行比較。
由于先驗是貝葉斯統(tǒng)計中的一個核心內(nèi)容，在接下來遇到新的問題時我們還會反復(fù)討論它，因此，如果你對前面的討論內(nèi)容感到有些疑惑，別太擔(dān)心，要知道人們在這個問題上已經(jīng)困惑了數(shù)十年并且相關(guān)的討論一直在繼續(xù)。
現(xiàn)在我們已經(jīng)有了后驗，相關(guān)的分析也就結(jié)束了。下面我們可能還需要對分析結(jié)果進(jìn)行總結(jié)，將分析結(jié)果與別人分享，或者記錄下來以備日后使用。
根據(jù)受眾不同，你可能在交流分析結(jié)果的同時還需要交流模型。以下是一種簡單表示概率模型的常見方式：
這是我們拋硬幣例子中用到的模型。符號～表示左邊隨機變量的分布服從右邊的分布形式，也就是說，這里θ服從于參數(shù)為α和β的Beta分布，而y服從于參數(shù)為n = 1和p = θ的二項分布。該模型還可以用類似Kruschke文中的圖表示成如下形式：
在第一層，根據(jù)先驗生成了θ，然后通過似然生成最下面的數(shù)據(jù)。圖中的箭頭表示變量之間的依賴關(guān)系，符號～表示變量的隨機性。
本書中用到的類似Kruschke中的圖都是由Rasmus B??th（http://www.sumsar.net/blog/2013/10/diy-kruschke-style-diagrams/）提供的模板生成的。
貝葉斯分析的結(jié)果是后驗分布，該分布包含了有關(guān)參數(shù)在給定數(shù)據(jù)和模型下的所有信息。如果可能的話，我們只需要將后驗分布展示給觀眾即可。通常，一個不錯的做法是：同時給出后驗分布的均值（或者眾數(shù)、中位數(shù)），這樣能讓我們了解該分布的中心，此外還可以給出一些描述該分布的衡量指標(biāo)，如標(biāo)準(zhǔn)差，這樣人們能對我們估計的離散程度和不確定性有一個大致的了解。拿標(biāo)準(zhǔn)差衡量類似正態(tài)分布的后驗分布很合適，不過對于一些其他類型的分布（如偏態(tài)分布）卻可能得出誤導(dǎo)性結(jié)論，因此，我們還可以采用下面的方式衡量。
一個經(jīng)常用來描述后驗分布分散程度的概念是最大后驗密度（ Highest Posterior Density，HPD）區(qū)間。一個HPD區(qū)間是指包含一定比例概率密度的最小區(qū)間，最常見的比例是95%HPD或98%HPD，通常還伴隨著一個50%HPD。如果我們說某個分析的HPD區(qū)間是[2, 5]，其含義是指：根據(jù)我們的模型和數(shù)據(jù)，參數(shù)位于2～5的概率是0.95。這是一個非常直觀的解釋，以至于人們經(jīng)常會將頻率學(xué)中的置信區(qū)間與貝葉斯方法中的可信區(qū)間弄混淆。如果你對頻率學(xué)的范式比較熟悉，請注意這兩種區(qū)間的區(qū)別。貝葉斯學(xué)派的分析告訴我們的是參數(shù)取值的概率，這在頻率學(xué)的框架中是不可能的，因為頻率學(xué)中的參數(shù)是固定值，頻率學(xué)中的置信區(qū)間只能包含或不包含參數(shù)的真實值。在繼續(xù)深入之前，有一點需要注意：選擇95%還是50%或者其他什么值作為HPD區(qū)間的概率密度比例并沒有什么特殊的地方，這些不過是經(jīng)常使用的值罷了。比如，我們完全可以選用比例為91.37%的HPD區(qū)間。如果你選的是95%，這完全沒問題，只是要記住這只是個默認(rèn)值，究竟選擇多大比例仍然需要具體問題具體分析。
對單峰分布計算95%HPD很簡單，只需要計算出2.5%和97.5%處的值即可：
def naive_hpd(post):
sns.kdeplot(post)
HPD = np.percentile(post, [2.5, 97.5])
plt.plot(HPD, [0, 0], label='HPD {:.2f} {:.2f}'.format(*HPD),
linewidth=8, color='k')
plt.legend(fontsize=16);
plt.xlabel(r'$\theta$', fontsize=14)
plt.gca().axes.get_yaxis().set_ticks([])
np.random.seed(1)
post = stats.beta.rvs(5, 11, size=1000)
naive_hpd(post)
plt.xlim(0, 1)
對于多峰分布而言，計算HPD要稍微復(fù)雜些。如果把對HPD的原始定義應(yīng)用到混合高斯分布上，我們可以得到：
np.random.seed(1)
gauss_a = stats.norm.rvs(loc=4, scale=0.9, size=3000)
gauss_b = stats.norm.rvs(loc=-2, scale=1, size=2000)
mix_norm = np.concatenate((gauss_a, gauss_b))
naive_hpd(mix_norm)
plt.savefig('B04958_01_08.png', dpi=300, figsize=(5.5, 5.5))
從上圖可以看出，通過原始HPD定義計算出的可信區(qū)間包含了一部分概率較低的區(qū)間，位于[0,2]。為了正確計算出HPD，這里我們使用了plot_post函數(shù)，你可以從本書附帶的代碼中下載對應(yīng)的源碼：
from plot_post import plot_post
plot_post(mix_norm, roundto=2, alpha=0.05)
plt.legend(loc=0, fontsize=16)
plt.xlabel(r'$\theta$', fontsize=14)
從上圖可以看出，95%HPD包含兩個區(qū)間，同時plot_post函數(shù)也返回了兩個眾數(shù)。
貝葉斯方法的一個優(yōu)勢是：一旦得到了后驗分布，我們可以根據(jù)該后驗生成未來的數(shù)據(jù)y，即用來做預(yù)測。后驗預(yù)測檢查主要是對觀測數(shù)據(jù)和預(yù)測數(shù)據(jù)進(jìn)行比較從而發(fā)現(xiàn)這兩個集合的不同之處，其目的是進(jìn)行一致性檢查。生成的數(shù)據(jù)和觀測的數(shù)據(jù)應(yīng)該看起來差不多，否則有可能是建模出現(xiàn)了問題或者輸入數(shù)據(jù)到模型時出了問題，不過就算我們沒有出錯，兩個集合仍然有可能出現(xiàn)不同。嘗試去理解其中的偏差有助于我們改進(jìn)模型，或者至少能知道我們模型的極限。即使我們并不知道如何去改進(jìn)模型，但是理解模型捕捉到了問題或數(shù)據(jù)的哪些方面以及沒能捕捉到哪些方面也是非常有用的信息。也許模型能夠很好地捕捉到數(shù)據(jù)中的均值但卻沒法預(yù)測出罕見值，這可能是個問題，不過如果我們只關(guān)心均值，這個模型對我們而言也還是可用的。通常我們的目的不是去聲稱一個模型是錯誤的，我們更愿意遵循George Box的建議，即所有模型都是錯的，但某些是有用的。我們只想知道模型的哪個部分是值得信任的，并測試該模型是否在特定方面符合我們的預(yù)期。不同學(xué)科對模型的信任程度顯然是不同的，物理學(xué)中研究的系統(tǒng)是在高可控條件下依據(jù)高深理論運行的，因而模型可以看做是對現(xiàn)實的不錯描述，而在一些其他學(xué)科如社會學(xué)和生物學(xué)中，研究的是錯綜復(fù)雜的孤立系統(tǒng)，因而模型對系統(tǒng)的認(rèn)知較弱。盡管如此，不論你研究的是哪一門學(xué)科，都需要對模型進(jìn)行檢查，利用后驗預(yù)測和本文學(xué)到的探索式數(shù)據(jù)分析中的方法去檢查模型。
本書用到的代碼是用Python 3.5寫的，建議使用Python 3的最新版本運行，盡管大多數(shù)代碼也能在更老的Python版本（包括Python 2.7）上運行，不過可能會需要做些微調(diào)。
本書建議使用Anaconda安裝Python及相關(guān)庫，Anaconda是一個用于科學(xué)計算的軟件分發(fā)，你可以從以下鏈接下載并了解更多：https://www.continuum.io/downloads。在系統(tǒng)上裝好Anaconda之后，就可以通過以下方式安裝Python庫了：
conda install NamePackage
我們會用到以下Python庫：
Ipython 5.0；
NumPy 1.11.1；
SciPy 0.18.1；
Pandas 0.18.1；
Matplotlib 1.5.3；
Seaborn 0.7.1；
PyMC3 3.0。
在命令行中執(zhí)行以下命令即可安裝最新穩(wěn)定版的PyMC3：
pip install pymc3
我們的貝葉斯之旅首先圍繞統(tǒng)計建模、概率論和貝葉斯定理做了一些簡短討論，然后用拋硬幣的例子介紹了貝葉斯建模和數(shù)據(jù)分析，借用這個經(jīng)典例子傳達(dá)了貝葉斯統(tǒng)計中的一些最重要的思想，比如用概率分布構(gòu)建模型并用概率分布來表示不確定性。此外我們嘗試揭示了如何選擇先驗，并將其與數(shù)據(jù)分析中的一些其他問題置于同等地位（怎么選擇似然，為什么要解決該問題等）。本文的最后討論了如何解釋和報告貝葉斯分析的結(jié)果。本文我們對貝葉斯分析的一些主要方面做了簡要總結(jié)，后面還會重新回顧這些內(nèi)容，從而充分理解和吸收，并為后面理解更高級的內(nèi)容打下基礎(chǔ)。下面我們會重點關(guān)注一些構(gòu)建和分析更復(fù)雜模型的技巧，此外，還會介紹PyMC3并將其用于實現(xiàn)和分析貝葉斯模型。
我們尚不清楚大腦是如何運作的，是按照貝葉斯方式？還是類似貝葉斯的某種方式？又或許是進(jìn)化過程中形成的某種啟發(fā)式的方式？不管如何，我們至少知道自己是通過數(shù)據(jù)、例子和練習(xí)來學(xué)習(xí)的，盡管你可能對此有不同的意見，不過我仍然強烈建議你完成以下練習(xí)。
（1）修改生成本文的第3個圖的代碼，在圖中增加一條豎線用來表示觀測到的正面朝上的比例（正面朝上的次數(shù)/拋硬幣的次數(shù)），將該豎線的位置與每個子圖中后驗的眾數(shù)進(jìn)行比較。
（2）嘗試用不同的先驗參數(shù)（beta值）和不同的實驗數(shù)據(jù)（正面朝上的次數(shù)和實驗次數(shù)）重新繪制本文的第3個圖。
（3）閱讀維基百科上有關(guān)Cromwell準(zhǔn)則的內(nèi)容：https://en.wikipedia.org/wiki/Cromwell%27s_rule。
（4）探索不同參數(shù)下高斯分布、二項分布和beta分布的圖像，你可以為每個分布單獨畫一個圖而不是全都畫在一個網(wǎng)格中。
PyMOL社區(qū)活躍者傾情奉獻(xiàn)！
發(fā)現(xiàn)Python貝葉斯分析的力量！
本書介紹了貝葉斯統(tǒng)計中的主要概念，以及將其應(yīng)用于數(shù)據(jù)分析的方法。本書采用編程計算的實用方法介紹了貝葉斯建模的基礎(chǔ)，使用一些手工構(gòu)造的數(shù)據(jù)和一部分簡單的真實數(shù)據(jù)來解釋和探索貝葉斯框架中的核心概念，然后在本書涉及的模型中，抽象出了線性模型用于解決回歸和分類問題，此外還詳細(xì)解釋了混合模型和分層模型，并單獨用一章討論了如何做模型選擇，最后還簡單介紹了非參模型和高斯過程。
本書所有的貝葉斯模型都用PyMC3實現(xiàn)。PyMC3是一個用于概率編程的Python庫，其許多特性都在書中有介紹。在本書和PyMC3的幫助下，讀者將學(xué)會實現(xiàn)、檢查和擴展貝葉斯統(tǒng)計模型，從而解決一系列數(shù)據(jù)分析的問題。
聯(lián)系客服
