貝葉斯統計的概念很簡單，有一些固定的數據（固定的意思是指我們無法改變觀測值），和一些感興趣的參數，剩下要做的就是探索這些參數可能的取值，其中所有的不確定性都通過概率進行建模。
說句白話，能夠用已有的資料做基礎，嘗試預測未來。
在別的統計學範式中，未知量有多種不同的表示方式，而在貝葉斯的框架中，所有未知量都是同等對待的，如果不知道某個變量，我們就給其賦予一個概率分布。因此，貝葉斯定理就是將先驗概率分布
（在觀測到數據之前我們對問題的理解）轉化成後驗分布
（觀測到數據之後所得到的信息），換句話說，貝葉斯統計就是一種機器學習的過程。
以上公式看不懂沒關係，其實就是說的由已經知道的，可以推導出大概率的結論。比如，小王以前偷盜成性，以後在監獄裡吃牢飯的機率就很高，俗話說的「江山易改，本性難移」就說的這個道理。
儘管概念上很簡單，純粹的概率模型得到的表達式通常分析起來很棘手。許多年來，這都是一個問題，大概也是影響貝葉斯方法廣泛應用的最大原因之一。隨著計算時代的到來，數值化方法的發展使得計算幾乎任意模型的後驗成為了可能，這極大地促進了貝葉斯方法的應用。我們可以把這些數值化方法當作通用引擎，按照PyMC3的核心開發者之一Thomas Wiecki的說法，只需要按下按鈕，推斷部分就可以自動完成了。
自動化推斷促進了概率程式語言的發展，從而使得模型構建和推斷相分離。在概率程式語言的框架中，用戶只需要寥寥數行代碼描述概率模型，後面的推斷過程就能自動完成了。概率編程使得人們能夠更快速地構建複雜的概率模型並減少出錯的可能，可以預見，這將給數據科學和其他學科帶來極大的影響。
我認為，程式語言對科學計算的影響可以與60多年前Fortran語言的問世相對比。儘管如今Fortran語言風光不再，不過在當時Fortran語言被認為是相當革命性的。科學家們第一次從計算的細節中解放出來，開始用一種更自然的方式關注構建數值化的方法、模型和仿真系統。類似地，概率編程將處理概率和推斷的過程對用戶隱藏起來，從而使得用戶更多的去關注模型構建和結果分析。
推斷引擎
即便某些情況下不太可能從分析的角度得到後驗，我們也有辦法將後驗計算出來，其中一些方法列舉如下。
（1）非馬爾科夫方法：
網格計算；二次近似；變分方法。（2）馬爾科夫方法：
Metropolis-Hastings算法；漢密爾頓蒙特卡洛方法（Hamiltonian Monte Carlo）/不掉向採樣（No U- Turn Sampler，NUTS）。如今，貝葉斯分析主要是通過馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法進行，同時變分方法也越來越流行，特別是在一些較大的數據集上。學習貝葉斯分析並不需要完全掌握這些方法，不過從概念層面上對它們的工作原理有一定了解會很有幫助，比如調試模型。
非馬爾科夫方法
我們首先討論基於非馬爾科夫方法的推斷引擎。非馬爾科夫方法在低維的問題上要比馬爾科夫方法更快，對於某些問題，這類方法非常有用，而對另外一些問題，這類方法只能提供真實後驗的粗略近似，不過沒關係，它們可以為馬爾科夫方法提供一個不錯的初始點，後面我們會介紹馬爾科夫的含義。
1．網格計算
網格計算是一種暴力窮舉的方法。即便你無法計算出整個後驗，你也可以根據一些點計算出先驗和似然。假設我們要計算某個單參數模型的後驗，網格近似可以按照如下方式進行：
確定參數的一個合理區間（先驗會給你點提示）；在以上區間確定一些網格點（通常是等距離的）；對於網格中的每個點計算先驗和似然。視情況，我們可能會對計算結果進行歸一化（把每個點的計算結果除以所有點的計算結果之和）。
很容易看出，選的點越多（網格越密）近似的結果就越好。事實上，如果使用無限多的點，我們可以得到準確的後驗。網格計算的方法不能很好地適用於多參數（又或者稱多維度）的場景，隨著參數的增加，採樣空間相比後驗空間會急劇增加，換言之，我們花費了大量時間計算後驗值，但對於估計後驗卻幾乎沒有幫助，因而使得該方法對於大多數統計學和數據科學的問題都不太實用
。
下面的代碼用網格計算的方法解決第一章中的拋硬幣問題：
def posterior_grid(grid_points=100, heads=6, tosses=9): """ A grid implementation for the coin-flip problem """ grid = np.linspace(0, 1, grid_points) prior = np.repeat(5, grid_points) likelihood = stats.binom.pmf(heads, tosses, grid) unstd_posterior = likelihood * prior posterior = unstd_posterior / unstd_posterior.sum() return grid, posterior假設4次拋硬幣實驗中只有一次觀測到正面朝上，那麼就有：
points = 15h, n = 1, 4grid, posterior = posterior_grid(points, h, n)plt.plot(grid, posterior, 'o-', label='heads = {}\ntosses = {}'.format(h, n))plt.xlabel(r'$\theta$')plt.legend(loc=0)
2．二次近似
二次近似，也稱 拉普拉斯方法 或者正態近似，利用的是高斯分布來近似後驗。該方法通常很有效，原因是後驗分布眾數附近的區域通常是服從正態分布的，事實上很多情況下就是高斯分布。二次近似的計算過程分為兩步。首先，找出後驗分布的峰值，該值可以通過一些最優化算法得到（有許多開箱即用的算法用來求解函數的最大值或最小值），這樣我們得到了用來近似的高斯分布的均值；然後估計函數在眾數附近的曲率，根據該曲率可以得出近似高斯分布的標準差。在介紹完PyMC3之後，我們會講解如何使用該方法。
3．變分方法
現代的貝葉斯統計學大多都採用馬爾科夫方法（下一節會介紹），不過該方法對某些問題求解很慢而且不能很好地並行計算。一種簡單的做法是同時運行多個馬爾科夫鏈，然後將結果合併，不過對大多數問題而言這並不是一個合適的解決方案，如何找到有效的並行計算方式是當前的一個研究熱點。
對於較大的數據集或是某些計算量很大的似然而言，變分方法是一個更好的選擇。此外，這類方法能快速得到後驗的近似，為MCMC方法提供初始點。
變分方法的基本思想是用一個更簡單的分布去近似後驗，這聽起來有點像拉普拉斯近似，不過深入細節你會發現二者有很大不同。變分方法的最大缺點是我們必須對每個模型設計一個特定的算法，因而變分方法並不是一個通用的推斷引擎，而是與模型相關的。
當然，許多人都在嘗試將變分方法自動化。最近提出的一個方法是 自動差分變分推斷 （Automatic Differentiation Variational Inference，ADVI）。從概念層面上講，ADVI是按下面步驟工作的。
（1）對參數進行變換從而使得它們位於整個實軸上。例如，假設一個參數限定為正數，對其求log後得到的值位於無界區間
內。
（2）用高斯分布對無界參數近似。需要注意，轉換後參數空間裡的高斯分布在原來參數空間內是非高斯的，這也是與拉普拉斯近似方法的不同點。
（3）採用某種優化算法使得高斯近似的結果儘可能接近後驗，該過程通過最大化 證據下界 （Evidence Lower Bound，ELBO）實現。如何衡量兩個分布的相似性以及證據下界的具體含義都是一些數學細節。
ADVI算法已經在PyMC3中實現了，本書後面會使用到。
馬爾科夫方法
有許多相關的方法都稱MCMC方法。對於網格計算近似方法而言，我們需要根據給定的點計算出似然和先驗，並且近似出整個後驗分布。MCMC方法要比網格近似更好，這是因為其設計思想是將更多的時間用於高概率區域的計算，而在較低概率區間花費較少時間。事實上，MCMC方法會根據相對概率訪問參數空間中的不同區間。如果區間 A 的概率是區間 B 的兩倍，那麼我們從區間 A 採樣的次數可能是從區間 B 採樣次數的兩倍。因此，即使我們無法從分析的角度得到整個後驗，我們也可以採用MCMC的方法從中採樣，並且採樣數越多，效果越好。
要理解MCMC方法，我們先將其拆分成兩個MC部分，即 蒙特卡洛 部分和 馬爾科夫鏈 部分。
1．蒙特卡洛
蒙特卡洛這部分可以用隨機數的應用來解釋。蒙特卡洛方法是一系列應用非常廣泛的算法，其思想是通過隨機採樣來計算或模擬給定過程。
蒙特卡洛是位於摩納哥公國的一個非常有名的城市，蒙特卡洛方法的開發者之一是Stanislaw Ulam。Stan正是基於這一核心思想——儘管很多問題都難以求解甚至無法準確用公式表達，但我們可以通過採樣或者模擬來有效地研究。
據傳，起因是為了計算單輪遊戲中的概率，解決該問題的方法之一是組合分析法；另一種是Stan聲稱的，進行多次單輪遊戲，最後計算其中有多少次是我們感興趣的。這聽起來似乎是顯而易見的，或者至少是相當合理的，比如，你可能已經用重採樣的方法來解決統計問題。不過這個實驗是早在70年前進行的，當時，第一臺計算機才剛開始研發。該方法首次應用於一個核物理問題。如今，個人計算機都已經足夠強大，用蒙特卡洛方法可以解決許多有趣的問題，因此，該方法也被廣泛應用到科學、工程、工業以及藝術等各個方面。
在使用蒙特卡洛方法計算數值的例子中，一個教科書上非常經典的是估計π。實際使用中有更好的方法來計算π，不過這個例子仍然具有教學意義。我們可以通過以下過程估計π：
（1）在邊長為2 R 的正方形內隨機撒 N 個點。
（2）在正方形內畫一個半徑為 R 的圓，計算在圓內的點的個數。
（3）得出
的估計值
。
需要注意：如果一個點滿足以下條件，我們認為該點位於圓內：
正方形的面積是(2 R ) 2 ，圓的面積是π R 2 ，因此二者的面積之比是4/π，而圓和正方形的面積分別正比於圓內的點的個數和總的點數 N 。我們可以通過幾行簡單的代碼來模擬該蒙特卡洛過程計算π值，同時計算出估計值與實際值之間的相對誤差：
N = 10000x, y = np.random.uniform(-1, 1, size=(2, N))inside = (x**2 + y**2) <= 1pi = inside.sum()*4/Nerror = abs((pi - np.pi)/pi)* 100outside = np.invert(inside)plt.plot(x[inside], y[inside], 'b.')plt.plot(x[outside], y[outside], 'r.')plt.plot(0, 0, label='$\hat \pi$ = {:4.3f}\nerror = {:4.3f}%'.format(pi, error), alpha=0)plt.axis('square')plt.legend(frameon=True, framealpha=0.9, fontsize=16);
上面的代碼中，outside是用來繪圖的，在計算的過程中沒有用到。另外一點需要澄清的是，由於這裡用的是單位圓，因此在判斷一個點是否在圓內時沒有計算平方根。
2．馬爾科夫鏈
馬爾科夫鏈是一個數學對象，包含一系列狀態以及狀態之間的轉移概率，如果每個狀態轉移到其他狀態的概率只與當前狀態相關，那麼這個狀態鏈就稱為馬爾科夫鏈。有這樣一個馬爾科夫鏈之後，我們可以任取一個初始點，然後根據狀態轉移概率進行隨機遊走。假設能夠找到這樣一個馬爾科夫鏈，其狀態轉移概率正比於我們想要採樣的分布（如貝葉斯分析中的後驗分布），採樣過程就變成了簡單地在該狀態鏈上移動的過程。那麼，如何在不知道後驗分布的情況下找到這樣的狀態鏈呢？有一個概念叫做 細節平衡條件 （Detailed Balance Condition），直觀上講，這個條件是說，我們需要採用一種可逆的方式移動（可逆過程是物理學中的一個常見的近似）。也就是說，從狀態 i 轉移到狀態 j 的概率必須和狀態 j 轉移到狀態 i 的概率相等。
總的來說就是，如果我們能夠找到滿足細節平衡條件的馬爾科夫鏈，就可以保證從中採樣得到的樣本來自正確的分布。保證細節平衡的最流行的算法是 Metropolis-Hasting算法 。
3．Metropolis-Hasting算法
為了更形象地理解這個算法，我們用下面這個例子來類比。假設我們想知道某個湖的水容量以及這個湖中最深的點，湖水很渾濁以至於沒法通過肉眼來估計深度，而且這個湖相當大，網格近似法顯然不是個好辦法。為了找到一個採樣策略，我們請來了兩個好朋友小馬和小萌。經過討論之後想出了如下辦法，我們需要一個船（當然，也可以是竹筏）和一個很長的棍子，這比聲吶可便宜多了，而且我們已經把有限的錢都花在了船上。
（1）隨機選一個點，然後將船開過去。
（2）用棍子測量湖的深度。
（3）將船移到另一個地點並重新測量。
（4）按如下方式比較兩點的測量結果。
如果新的地點比舊的地點水位深，那麼在筆記本上記錄下新的測量值並重複過程（2）。如果新的地點比舊的地點水位淺，那麼我們有兩個選擇：接受或者拒絕。接受意味著記錄下新的測量值並重複過程（2）；拒絕意味著重新回到上一個點，再次記錄下上一個點的測量值。如何決定接受還是拒絕新的測量值呢？這裡的一個技巧便是使用Metropolis-Hastings準則，即接受新的測量值的概率正比於新舊兩點的測量值之比。
按照以上過程迭代下去，我們不僅可以得到整個湖的水容量和最深的點，而且可以得到整個湖底的近似曲率。你也許已經猜到了，在這個類比中，湖底的曲率其實就是後驗的分布，而最深的點就是後驗的眾數。根據小馬的說法，迭代的次數越多，近似的效果越好。
事實上，理論保證了在這種情形下，如果我們能採樣無數次，最終能得到完整的後驗。幸運地是，實際上對於很多問題而言，我們只需要相對較少地採樣就可以得到一個相當準確的近似。
現在讓我們從更正式的角度來看看該算法。對於很多分布而言（如高斯分布），我們有相當高效的算法對其採樣，但對於一些其他分布，情況就變了。Metropolis-Hastings算法使得我們能夠從任意分布中以概率 p ( x )得到採樣值，只要我們能算出某個與 p ( x )成比例的值。這一點很有用，因為在類似貝葉斯統計的許多問題中，最難的部分是計算歸一化因子，也就是貝葉斯定理中的分母。Metropolis-Hastings算法的步驟如下。
（1）給參數 x i 賦一個初始值，通常是隨機初始化或者使用某些經驗值。
（2）從某個簡單的分布
中選一個新的值
，如高斯分布或者均勻分布。這一步可以看做是對狀態
的擾動。
（3）根據Metropolis-Hastings準則計算接受一個新的參數值的概率：
（4）從位於區間[0,1]內的均勻分布中隨機選一個值，如果第（3）步中得到的概率值比該值大，那麼就接受新的值，否則仍保持原來的值。
（5）回到第（2）步重新迭代，直到我們有足夠多的樣本，稍後會解釋什麼叫足夠多。
有幾點需要注意。
如果選取的分布 是對稱的，那麼可以得到 ，這通常稱為 Metropolis準則 。步驟（3）和步驟（4）表明：我們總是會轉移到一個比當前狀態（或參數）概率更大的狀態（或參數），對於概率更小的，則會以 x i +1 與 x i 之比的概率接受。該準則中的接受步驟使得採樣過程相比網格近似方法更高效，同時保證了採樣的準確性。目標分布（貝葉斯統計中的後驗分布）是通過記錄下來的採樣值來近似的。如果我們接受轉移到新的狀態 x i +1 ，那麼我們就記錄該採樣值 x i +1 。如果我們拒絕轉移到 x i +1 ，那麼我們就記錄 x i 。最後，我們會得到一連串記錄值，有時候也稱 採樣鏈 或者 跡 。如果一切都正常進行，那麼這些採樣值就是後驗的近似。在採樣鏈中出現次數最多的值就是對應後驗中最可能的值。該過程的一個優點是：後驗分析很簡單，我們把對後驗求積分的過程轉化成了對採樣鏈所構成的向量求和的過程。
下面的代碼展示了Metropolis算法的一個基本實現。這段代碼並不是為了解決什麼實際問題，只是在這裡用來演示，如果我們知道怎麼計算給定點的函數值，我們就可能得到該函數的採樣。需要注意下面的代碼中不包含貝葉斯相關的部分，既沒有先驗也沒有數據。要知道，MCMC是一類能夠用於解決很多問題的通用方法。例如，在一個（非貝葉斯的）分子模型中，我們可能需要一個函數來計算在某個狀態 x 下系統的能量而不是簡單地調用func.pdf(x)函數。
metropolis函數的第一個參數是一個SciPy的分布，假設我們不知道如何從中直接採樣。
def metropolis(func, steps=10000): """A very simple Metropolis implementation""" samples = np.zeros(steps) old_x = func.mean() old_prob = func.pdf(old_x) for i in range(steps): new_x = old_x + np.random.normal(0, 0.5) new_prob = func.pdf(new_x) acceptance = new_prob/old_prob if acceptance >= np.random.random(): samples[i] = new_x old_x = new_x old_prob = new_prob else: samples[i] = old_x return samples
下面的例子中，我們把func定義成一個beta函數，因為beta函數可以通過改變參數調整分布的形狀。我們把metropolis函數的結果用直方圖畫出來，同時用紅色的線表示真實的分布。
func = stats.beta(0.4, 2)samples = metropolis(func=func)x = np.linspace(0.01, .99, 100)y = func.pdf(x)plt.xlim(0, 1)plt.plot(x, y, 'r-', lw=3, label='True distribution')plt.hist(samples, bins=30, normed=True, label='Estimated distribution')plt.xlabel('$x$', fontsize=14)plt.ylabel('$pdf(x)$', fontsize=14)plt.legend(fontsize=14)
現在你應該從概念上掌握了Metropolis-Hastings算法。也許你需要回過頭去重新閱讀前面幾頁才能完全消化。此外，我還強烈建議閱讀PyMC3核心作者之一寫的博文
。他用一個簡單的例子實現了metropolis方法，並將其用於求解後驗分布，文中用非常好看的圖展示了採樣的過程，同時簡單討論了最初選取的步長是如何影響結果的。
4．漢密爾頓蒙特卡洛方法/不掉向採樣
MCMC方法，包括Metropolis-Hastings，都在理論上保證如果採樣次數足夠多，最終會得到後驗分布的準確近似。不過，實際中想要採樣足夠多次可能需要相當長的時間，因此，人們提出了一些Metropolis-Hastings算法的替代方案。這些替代方案，包括Metropolis-Hastings算法本身，最初都是用來解決統計力學中的問題。
統計力學是物理學的一個分支，主要研究原子和分子系統的特性。漢密爾頓蒙特卡洛方法，又稱混合蒙特卡洛（Hybrid Monte Carlo，HMC），是這類改進方案之一。簡單來說，漢密爾頓這個詞描述的是物理系統的總能量，而另外一個名稱中的「混合」是指將Metropolis-Hastings算法與分子力學（分子系統中廣泛應用的一種仿真技巧）相結合。
HMC方法本質上和Metropolis-Hastings是一樣的，改進的地方在於：原來是隨機放置小船，現在有了一個更聰明的辦法，將小船沿著湖底方向放置。為什麼這個做法更聰明？因為這樣做避免了Metropolis-Hastings算法的主要問題之一：探索得太慢而且採樣結果自相關（因為大多數採樣結果都被拒絕了）。
那麼，如何才能不必深入其數學細節而理解漢密爾頓蒙特卡洛方法呢？假設我們還是在湖面上坐著船，為了決定下一步將要去哪，我們從當前位置往湖底扔了一個球，受「球狀奶牛」的啟發，我們假設球面是理想的，沒有摩擦，因而不會被泥巴和水減速。扔下球之後，讓它滾一小會兒，然後把船劃到球所在的位置。現在利用Metropolis-Hastings算法中提到的Metropolis準則來選擇接受或者拒絕，重複整個過程一定次數。改進後的過程有更高的概率接受新的位置，即使它們的位置相比前一位置距離較遠。
現在跳出我們的思維實驗，回到現實中來。基於漢密爾頓的方法需要計算函數的梯度。梯度是在多個維度上導數的推廣。我們可以用梯度信息來模擬球在曲面上移動的過程。因此，我們面臨一個權衡；HMC計算過程要比Metropolis-Hastings更複雜，但是被接受概率更高。對於一些複雜的問題，HMC方法更合適一些。HMC方法的另一個缺點是：想要得到很好的採樣需要指定一些參數。如果手動指定，需要反覆嘗試，這要求使用者有一定的經驗。幸運地是，PyMC3中有一個相對較新的不掉向採樣算法，該方法被證實可以有效提升HMC方法的採樣效率，同時不必手動調整參數。
5．其他MCMC方法
MCMC的方法很多，而且人們還在不斷提出新的方法。如果你認為你能提升採樣效率，會有很多人對你的想法感興趣。討論所有這類方法及其優缺點顯然超出了本書的範圍。不過，有些方法仍值得一提，因為你可能經常會聽到人們討論它們，所以最好是了解下他們都在討論些什麼。
在分子系統模擬中廣泛應用的另外一種採樣方法是 副本交換 （Replica Exchange），也稱 並行退火 （Parallel Tempering）或者 Metropolis Coupled MCMC （MC3；好多MC……）。
該方法的基本思想是並行模擬多個不同的副本，每個副本都按照Metropolis-Hastings算法執行。多個副本之間的唯一不同是一個叫做溫度的參數（又受到物理學的啟發！），該參數用來控制接受低概率採樣點的可能性。
每隔一段時間，該方法都嘗試在多個副本之間進行切換，切換過程同時遵循Metropolis-Hastings 準則來接受或拒絕，只不過現在考慮的是不同副本之間的溫度。切換過程可以在狀態鏈上進行隨機切換，不過，通常更傾向於在相鄰副本之間切換，也就是說，具有相似溫度的副本會有更高的接受概率。
該方法的直觀理解是：如果我們提高溫度，接受新位置的概率也會提升，反之則會降低。溫度更高的副本可以更自由地探索系統，因為這些副本的表面會變得相當平坦從而更容易探索。對於溫度無限高的副本，所有狀態都是等概率的。副本之間的切換避免了較低溫度的副本陷在局部最低點，因而該方法很適合探索有多個最低點的系統。
看了這麼多對Python和數學感興趣嗎？可以點擊下面購買對應的書籍。