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