在這篇文章中，我從非技術層面來介紹一下Markov chain Monte Carlo，通常簡稱為MCMC。MCMC被廣泛應用于貝葉斯統計模型擬合。MCMC有很多不同的方法，這里我將主要介紹Metropolis–Hastings 算法。為了簡便起見,我將 忽略一些細節，強烈推薦大家在練習使用MCMC之前閱讀[BAYES] 手冊。
我們繼續前一篇文章貝葉斯統計介紹Part 1：基本概念中提到的硬幣投擲例子。我們對參數θ的后驗分布感興趣，也就是投擲硬幣時，“頭”在上面的概率。我們的先驗分布是扁平的,無信息Beta分布參數1和1。然后我們將使用二項式似然函數對我們的實驗數據進行量化，10次投擲硬幣4次“頭”朝上的結果。我們可以使用MCMC的M-H算法來生成后驗分布θ的樣本。然后可以使用這個樣本來估算如后驗分布的均值等問題。
這種技術有三個基本部分：
1. Monte Carlo
2. Markov chains
3. M–H algorithm
Monte Carlo methods
蒙特卡羅是指依靠偽隨機數的生成方法(簡稱偽隨機數)。圖1說明了蒙特卡羅實驗的一些特點。
Figure 1: Proposal distributions, trace plots, and density plots
在這個例子中，我將從正態分布的均值0.5和任意方差σ中生成一系列隨機數。正態分布稱為建議分布。
圖中右上角部分稱為軌跡圖，它會根據繪圖的順序來顯示θ值。它也顯示建議分布順時針旋轉90度的情況，我每畫一個θ值就會將它向右平移。綠點軌跡圖顯示了θ的當前值。
左上角的密度圖是逆時針旋轉90度的直方圖樣本。我旋轉了坐標軸，因此θ軸線是平行的。同樣的，綠色的點密度圖顯示了θ的當前值。
讓我們繪制10,000個θ隨機值來看看是如何工作的。
Animation 1: A Monte Carlo experiment
動畫圖1說明了幾個重要的特點。首先建議分布不會從一個迭代中更改到另一個中。其次，隨著樣本容量的增加，密度圖看起來越來越像建議分布。第三，軌跡圖是平穩的--所有迭代的變量，變化形式看起來都是一樣的。
蒙特卡羅仿真生成的樣本看起來像建議分布，很多時候也正是我們所需的。但是我們需要一個額外的工具從后驗分布中生成一個樣本。
Markov chain Monte Carlo methods
Markov chain是一個數字序列,每個數字都依賴于序列中的前一個數。例如，我們從一個正常的建議分布與平均等于θ的前一個值來繪制θ值。圖2顯示了 軌跡圖和密度圖，當前的θ值(θt=0.530)是從一個平均等于θ前值的建議分布(θt−1=0.712)來繪制的。
Figure 2: A draw from a MCMC experiment
下面圖3顯示序列中的下一次迭代的軌跡圖和密度圖。建議密度的均值現在是θt−1=0.530，并且我們從新的建議分布來繪制出了隨機數值θt=0.411。
Figure 3: The next iteration of the MCMC experiment
讓我們使用MCMC繪制10,000個θ隨機值來看看是如何工作的。
Animation 2: A MCMC experiment
動畫圖2顯示了Monte Carlo實驗和MCMC實驗直接的差異。首先，建議分布是隨著每個迭代變化而變化的。這就生成一個隨機游走模式的軌跡圖：所有迭代的變化是不一樣的。其次，產生的密度圖看上去不像建議分布或其他任何有用的分布，肯定也不是一個后驗分布。
MCMC從后驗分布獲取樣本失敗的原因是我們還沒在生成樣本的過程中輸入過任何后驗分布的信息。我們可以通過保持θ值來提高樣本，θ值更可能是在后驗分布下的值而不太可能是丟棄值。
但是基于后驗分布，θ值接受和拒絕建議值最明顯的困難是我們通常不知道后驗分布的函數形式。如果我們知道它的函數形式，就可以直接計算概率，而不用生成一個隨機樣本。那么我們如何不用知道函數形式就可以接受或拒絕建議θ值呢? 答案就是M-H算法!
MCMC和M–H 算法
M-H算法可以用在我們不知道后驗分布函數形式的情況下來判斷是否接受θ建議值。讓我們將M-H算法分成幾步，然后經過幾次迭代來看看它是怎么實現的。
Figure 4: MCMC with the M–H algorithm: Iteration 1
圖4顯示了等于(θt−1=0.0.517)的建議分布軌跡圖和密度圖。我們已經繪制了一個預估值θ(θnew=0.380)，我將這個值定義為θnew是因為還沒有確定是否要用這個值。
我們從步驟1計算比率開始使用M-H算法
我們不知道后驗分布的函數形式，但是在這個例子中，我們可以替代先驗分布和似然函數的乘積。比率的計算并不是都那么簡單，但是我們可以把事情變得簡單些。
圖4步驟1顯示了(θnew=0.380) 和(θt−1=0.0.517) 的后驗概率比等于1.307。步驟2我們計算接受概率α(θnew,θt−1)，也就是最低比值r(θnew, θt -1)和1。步驟2是必須的，因為比率肯定在0和1之間。
接受概率等于1，所以我們接受(θnew=0.380),建議分布的平均值在下個迭代中就變成了0.380.
Figure 5: MCMC with the M–H algorithm: Iteration 2
圖5顯示了下一個迭代(θnew=0.286)是根據平均值為(θt−1=0.380)后驗分布繪制的。步驟1
計算的比率r(θnew,θt−1)等于0.747，小于1。步驟2中計算的接受概率α(θnew,θt−1)等于0.747。
我們不會自動拒絕θnew，因為接受概率小于1。相反，我們可以繪制來自步驟3中的均勻分布(0，1)的隨機數U。如果U小于接受概率，我們就接受建議值θnew，反之，我們拒絕θnew并保留θt−1。
在這里u=0.094小于了接受概率(0.747)，所以我們就接受(θnew=0.286)。我在圖5中將θnew和0.286以綠色來表示它們已經被接受。
Figure 6: MCMC with the M–H algorithm: Iteration 3
圖6顯示了下一個迭代(θnew=0.088)是根據平均值為(θt−1=0.286)后驗分布繪制的。步驟1計算的比率r(θnew,θt−1)等于0.039，小于1. 步驟2中計算的接受概率α(θnew,θt−1)等于0.039.步驟3中計算的U 的值是0.247，大于接受概率。因此在步驟4中，我們拒絕θnew=0.088并保留θt−1=0.286。
讓我們使用M-H算法的MCMC繪制10,000個θ隨機值來看看是如何工作的。
Animation 3: MCMC with the M–H algorithm
動畫圖3表明了幾個問題。首先，建議分布會隨著大部分迭代變化而變化。需要注意的是, 如果θnew值被接受，它們會以綠色表示，如果被拒絕會以紅色表示。其次，我們只使用MCMC的時候，軌跡圖不顯示隨機游走模式，所有迭代中的變化是相似的。最后，密度圖看起來像一個有用的分布。
旋轉一下密度圖結果，更仔細的看看。
Figure 7: A sample from the posterior distribution generated with MCMC with the M–H algorithm
圖7顯示了我們使用M-H算法MCMC生成的樣本直方圖。在這個例子中，我們知道后驗分布的參數是Beta 5和7. 紅色線條顯示后驗分布分析，藍色線條是樣本的核密度圖。核密度圖跟Beta(5,7)分布相似，這說明我們的樣本是一個很好的近似的理論后驗分布。
我們可以使用樣本來計算后驗分布均值或中位數，95%的置信區間的概率，或θ落在任意區間的概率。
Stata中的MCMC和M-H算法
讓我們回到使用bayesmh的投擲硬幣實例中。隨著二項式的可能性，我們指定一個beta分布參數1和1。
Example 1: Using bayesmh with a Beta(1,1) prior
輸出結果告訴我們bayesmh跑完了12,500個MCMC迭代。“Burn-in”告訴我們為了減少鏈中隨機起始值的影響，迭代中的2500個被丟棄。下面那行告訴我們，最終的MCMC樣本大小為10,000, 再下面一行表明我們的數據集有10個觀測值。接受概率θ建議值比例包含在最終的MCMC樣本中。我建議您閱讀Stata Bayesian Analysis Reference Manual 并了解效率的定義，但是高效率表明低自相關，而低效率表明高自相關。Monte Carlostandard error (MCSE)顯示在系數表中，估算后驗均值的近似誤差。
檢查chain的收斂性
“convergence”這個詞在MCMC和極大似然中的表述有不同的含義。用于最大似然估計算法迭代直至收斂到最大。MCMC鏈不迭代直到確定最佳值。外鏈簡單迭代直到達到要求的樣本量大小，任何算法停止。事實上，外鏈停止運行并不表明后驗分布中最佳樣本的生成。我們必須檢查樣本以制止問題的出現。我們可以用bayesgraph diagnostics圖形形式來檢查樣本。
圖8顯示包含MCMC樣本的軌跡圖、直方圖和密度圖的診斷圖和correlegram。軌跡圖有一個平穩模式，這也正是我們希望看到的。動畫圖2中的隨機游走模式顯示了chain存在的問題。直方圖沒有任何特殊特征如多模式。完整樣本的核密度圖，chain的前半部分，chain的后半部分都看起來很相似并沒有任何特殊的特征如chain的前半部分和后半部分密度不同。使用Markov chain來生成樣本并產生一個自相關，但是自相關會隨著較大的滯后值迅速減小。這些圖形并不能說明我們的樣本有任何問題。
總結
這篇文章介紹了使用M-H算法MCMC背后的想法。請注意我省略了一些細節，忽略了一些假定條件以便讓事情變得簡單，順著我們的感覺往下發展。Stata的bayesmh命令實際上執行了一個更難的算法，我們稱之為帶M-H算法的自適應MCMC。但是基本概念是一樣的，希望能給大家一些啟發。
