一.蒙特卡羅法的缺陷
通常的蒙特卡羅方法可以模擬生成滿足某個分布的隨機向量，但是蒙特卡羅方法的缺陷就是難以對高維分布進行模擬。對於高維分布的模擬，最受歡迎的算法當屬馬爾科夫鏈蒙特卡羅算法(MCMC)，他通過構造一條馬爾科夫鏈來分步生成隨機向量來逼近制定的分布，以達到減小運算量的目的。
二.馬爾科夫鏈方法概要
馬爾科夫鏈蒙特卡羅方法的基本思路就是想辦法構造一個馬爾科夫鏈，使得其平穩分布是給定的某分布，再逐步生模擬該馬爾科夫鏈產生隨機向量序列。其基本思路如下。就像是普通的蒙特卡羅方法本質上依賴於概率論中的大數定理，蒙特卡羅方法的理論支撐是具有遍歷性的馬爾科夫鏈的大數定理。馬爾科夫鏈蒙特卡羅方法的大體思路如下：
(1)給定某個分布$p(x)$, 構造某個馬爾科夫鏈$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$使得$p$是其平穩分布，且滿足一定的特殊條件；
(2)從一點$x_{0}$出發，依照馬爾科夫鏈$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$隨機生成向量序列$x_{0},x_{1},...$；
(3)蒙特卡羅積分估計：計算$E_{p}(f)\approx\sum_{t=1}^{N}f(x_{t})$
三.MCMC的數學基礎——馬爾科夫鏈的遍歷性，大數定理
MCMC為什么可以近似計算積分? 其實在數學上這是不太平凡的，下面簡要介紹一下其數學理論依據。
3.1 馬爾科夫鏈與其遍歷性, 馬爾科夫鏈的大數定理：
所謂馬爾科夫鏈通俗的說就是一個隨機過程，其滿足，t時刻的狀態和$t-1$之前的狀態無關。我們用嚴格的測度論語言說就是：
定義3.1：定義於概率空間$(\Omega,\mathcal{G},P)$, 取值於$\mathcal{Y}\in\mathbb{R}^{K}$的隨機向量序列$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$稱為離散時間馬爾科夫鏈(Markov Chain of discrete time)如果其滿足：
對於任意$\mathcal{Y}$的Borel集$B\in \mathcal{B}_{\mathcal{Y}}$
$P(X_{t+1}^{-1}(B)\mid X_{t},...,X_{1})=P(X_{t+1}^{-1}(B)\mid X_{t})$
進一步的，如果$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$還滿足：
\begin{equation}P(X_{t+1}^{-1}(B)\mid X_{t})=P(X_{1}^{-1}(B)\mid X_{0})\end{equation}
我們稱馬爾科夫鏈$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$為時間齊次（time homogeneous)的，這時我們定義該馬爾科夫鏈的轉移核(transition kernel)$$P_{t}: \mathbb{N}\times\mathcal{B}_{\mathcal{Y}}\longrightarrow [0,1]:$$
$$P_{t}(y,A)\triangleq P(X_{t}\in A\mid X_{0}=y),$$
對任意$t\in\mathbb{N}$, 並且我們直接簡記$P(y,A)=P_{1}(y,A)$, 對$y\in\mathcal{Y}$, $A\in\mathcal{B}_{\mathcal{Y}}$。
為了方便起見，我們在這里的所有“馬爾科夫鏈”均為離散時間，時間齊次的馬爾科夫鏈。
如果這時候$\mathcal{Y}$上已經有一個定義於$\mathcal{Y}$上所有Borel集構成的$\sigma$代數$\mathcal{B}_{\mathcal{Y}}$上的$\sigma$有限測度$\nu$以及某個概率分布函數$p$, $\int_{\mathcal{Y}}pd\nu=1$, 則馬爾科夫鏈$\lbrace X_{t}\rbrace_{t=1}^{\infty}$稱為具有以$p$為平穩分布的遍歷性（ergodicity），如果其滿足：
1)非周期性(Aperiodicity)：不存整數$T>1$以及互不相交的Borel集$B_{i}\in\mathcal{B}_{\mathcal{Y}}$, $i=1,...,T$使得:
對任意$y_{i}\in B_{i}$, $j\equiv i+1(\text{mod } T)$, 有:
\begin{equation}P(y_{i},B_{j})=1,\end{equation}
2)$p-$平穩性($p$-Invariance): 對任意的$A\in\mathcal{B}_{\mathcal{Y}}$
\begin{equation}\int_{\mathcal{Y}}P(y,A)p(y)d\nu(y)=\int_{\mathcal{Y}}p(y)d\nu(y)\end{equation}
3)$p-$不可約性($p$-Irreducibility):
對任意的$y\in\mathcal{Y}$, $A\in\mathcal{B}_{\mathcal{Y}}$, $\int_{A}p d\nu>0$, 存在$t\in\mathbb{N}$使得
\begin{equation}P_{t}(y,A)>0\end{equation}
4)Harris 回歸性(Harris Recurrence): 對任意$A\in\mathcal{B}_{\mathcal{Y}}$, $\int_{A}p(y)d\nu(y)>0$, 馬爾科夫鏈$X$以概率為1，無數次經過$A$:
\begin{equation}P(\sum_{t=1}^{\infty}\text{I}_{\lbrace X_{t}\in A\rbrace}=+\infty\mid X_{0}=y)=1\end{equation}
3.2 馬爾科夫鏈的大數定理
我們知道，經典的大數定理要求是隨機序列$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$滿足獨立同分布假設。但是問題來了，對於馬爾科夫鏈，自然獨立同分布假設一般是不滿足的，那么憑什么大數定理成立，也就是我們可以生成馬爾科夫鏈的樣本來近似計算積分? 下面的馬爾科夫鏈的大數定理(large number theorem of markov chains)是MCMC算法的數學基礎，回答了以上問題：
定理3.1（馬爾科夫鏈的大數定理）：如果馬爾科夫鏈$\lbrace X_{t}\rbrace$具有以概率分布函數$p$為平穩分布的遍歷性，$g\in L^{1}_{pd\nu}(\mathcal{Y})$，那么我們有幾乎處處收斂：
\begin{equation}\frac{1}{m}\sum_{t=1}^{m}g(X_{t})\longrightarrow\int_{\mathcal{Y}}g(y)p(y)d\nu,\end{equation}
並且：
\begin{equation}P_{t}(y, A)\longrightarrow \int_{A}pd\nu\end{equation}
對於p-幾乎處處的$y\in\mathcal{Y}$.
馬爾科夫鏈的大數定理是概率論中極其深刻而優雅的定理，其證明相當復雜，完整的證明也許得寫本小書專門講，參見[3]的定理4.3。
馬爾科夫鏈的大數定理告訴我們以下的事實：
如果馬爾科夫鏈有以$p$為平穩分布的遍歷性，則對幾乎處處（幾乎所有）的$y$, 從$y$狀態出發，$X_{t}$漸進趨於均衡分布$p$，這是非常有趣的性質，那就是漸進狀態和初始值無關，這也自然的提供我們模擬生成服從$p$分布隨機數的思路。另外，（6）保證了估算積分值的合法性。
四. Metropolis Hastings 算法：
現在我們考慮如何實現一個MCMC算法。現在我們給定了某個概率分布函數$p$，其實整個算法實施的關鍵僅僅在於如何構造一個具有以$p$為平穩分布遍歷性的馬爾科夫鏈，下面的Metropolis-Hastings算法給出一種方法。
Metropolis-Hatings算法的歷史很耐人尋味，最初由物理學家在第二次世界大戰期間於洛斯阿拉莫斯實驗室研制原子彈時發現，於1953年首次提出於由Nicholas Metropolis, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, Edward Teller五人署名，短短四頁的一篇發表於某化學雜志的文章里[6]，他們在這篇文章里模擬了Boltzmann分布的采樣，一個推廣工作，也是現今更加廣泛采用的版本是由W.Hatings於1970年給出[7]，一個特殊的版本於1984年在S.Geman, D.Geman於伊辛模型的研究中獨立提出,現在被稱為Gibbs采樣算法[8]。但是這類算法並沒有得到統計學界的足夠關注，一直到1990年Gelfand, Smith的工作[9]。該算法名聲大噪之后，最初的五位提出者陷入曠日持久的名譽爭奪戰，因為其他人對該算法僅僅命名Metropolis-Hastings算法感到不滿。
我們在這里沿用之前的記號和約定。這時如果存在可測函數$p(\cdot,\cdot):\mathcal{Y}\times\mathcal{Y}\longrightarrow \mathbb{R}^{+}$, 使得馬爾科夫鏈$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$的轉移核滿足：
\begin{equation}P(y, A)=\int_{A}p(y,z)d\nu(z),\end{equation}
對任意$(y, A)\in\mathcal{Y}\times\mathcal{B}_{\mathcal{Y}}$，則稱該函數為馬爾科夫鏈$\lbrace X_{t}\rbrace$的轉移分布函數(triansition distribution function)。
Metropolis Hastings算法的具體思路是模擬一個轉移分布函數$p:\mathcal{Y}\times\mathcal{Y}\longrightarrow \mathbb{R}$滿足：
\begin{equation}p(y)p(y,z)=p(z)p(z,y)\end{equation}
的馬爾科夫鏈，以上條件稱作細致平衡(detailed balance)。我們知道，如果細致平衡條件滿足，則對任意的$A\in\mathcal{B}_{\mathcal{Y}}$兩邊取積分$\int_{\lbrace (y,z)\in\mathcal{Y}\times A\rbrace} d\nu(y)\times d\nu(z)$我們得到：
\begin{equation}\int_{\lbrace (y,z)\in\mathcal{Y}\times A\rbrace} p(y)p(y,z)d\nu(y)\times d\nu(z)=\int_{\lbrace (y,z)\in \mathcal{Y}\times A\rbrace}p(y)p(y,z)d\nu(y)\times d\nu(z),\end{equation}
由Fubini定理容易計算上式左邊為$\int_{A}p(y)P(y,A)d\nu(y)$, 右邊為$\int_{A}p d\nu$, 所以：
\begin{equation}\int_{A}p d\nu=\int_{A}p(y)P(y,A)d\nu(y),\end{equation}
所以這時相應的以$p(\cdot,\cdot)$為轉移分布函數的馬爾科夫鏈必然是$p-$不變的，也就是說細致平衡條件是比$p$-不變性更加強的限制條件。（物理意義是什么？以后慢慢了解。）
現在的問題是我們如何構造這樣的一個$p(\cdot,\cdot)$?下面我們探討一下。我們先做一個如下的模擬實驗：
和常見的蒙特卡羅方法一樣，假如我們已經有一個提議分布(proposal distribution) $q(y, z)$, 是一個常見的，容易進行模擬的分布，那么這時候由狀態$y$出發下一步模擬以分布$q(y, \cdot)$生成一個向量$z$, 一般來說$z\neq y$, 這時如果我們繼續做模擬，也就是模擬$z$ 以某概率(待定)$\alpha(y, z)$ 生成 $z$，而以概率$1-\alpha(y, z)$生成$y$,
現在我們來看一下，這樣模擬后得到的轉移分布函數：
\[p(y,z)=\begin{cases}q(y,z)\alpha(y,z)&\text{if } z\neq y,\\g(y)&\text{if }z=y,\end{cases}\]
其中$q(y)$是某個關於$y$的函數使得$\int_{\mathcal{Y}}p(y,z)d\nu(z)=1$, 例如我們可以取$q(y)=(1-\int_{\mathcal{\lbrace z\mid z\in\mathcal{Y},z\neq y\rbrace}}q(y,z)\alpha(y,z)d\nu(z) )/\nu(\lbrace y\rbrace)$, 當$\nu(\lbrace y\rbrace)>0$, 而直接取$0$當$\nu(\lbrace y\rbrace)=0.$
這時候我們要求$p(\cdot,\cdot)$滿足細致平衡條件，也就是要求：
\begin{equation}p(y)q(y,z)\alpha(y,z)=p(z)q(z,y)\alpha(z,y)\end{equation}
由上面的公式，我們知道，如果能夠找到某個對稱函數$\lambda(y,z)$使得：
$$\alpha(y,z)=\lambda(y,z)p(z)q(z,y),$$
則（12）自動滿足。而想其成為概率值，必然有$\lambda(y,z)\leq \frac{1}{p(z)q(z,y)}$, 再由對稱性$\lambda(y,z)=\lambda(z,y)\leq \frac{1}{p(y)p(y,z)}$, 於是乎我們只需要取$\lambda(y,z)=\min\lbrace \frac{1}{p(z)q(z,y)}, \frac{1}{p(y)p(y,z)}\rbrace,$
這時候自然有：
\begin{equation}\alpha(y,z)=\min\lbrace \frac{p(z)q(z,y)}{p(y)q(y,z)},1\rbrace\end{equation}
所以依照以上方法，我們確實可以得到某個$p(\cdot,\cdot)$滿足細致平衡條件，然后我們就可以模擬產生以此為轉移分布函數的馬爾科夫鏈的樣本。
我們總結Metropolis Hasting算法如下：
算法4.1(Metropolis Hasting)：
任意選擇某個$x_{0}$, 我們假設現在已經生成了$x_{0},...,x_{t}$, $t\in\mathbb{N}$, 現在我們迭代生成$x_{t+1}$, 以如下方式：
(1)按照提議分布$q(x_{t},\cdot)$生成某數$y_{t}$
(2)計算出：
\begin{equation}\alpha(x_{t},y_{t})=\min\lbrace \frac{p(y_{t})q(y_{t},x_{t})}{p(x_{t})q(x_{t},y_{t})},1\rbrace\end{equation}
(3)隨機生成某數$u\sim U([0,1])$，如果$u\leq \alpha(x_{t},y_{t})$，令$x_{t+1}=y_{t}$，否則令$x_{t+1}=x_{t}$.
以上是Metropolis Hastings算法的直觀推導，但是問題來了，這時候構造的馬爾科夫鏈是不是真的是遍歷的? 如果是的話我們才可以用該算法估算積分值。幸運的是，可以證明該種算法能確保相應馬爾科夫鏈的遍歷性：
定理4.1：如果馬爾科夫鏈$\lbrace X_{t}\rbrace_{t\in\mathbb{N}}$的轉移分布函數$p(\cdot,\cdot)$為上述Metropolis-Hastings算法中的方法構造的函數, 則該馬爾科夫具有遍歷性。
定理的證明在這里不贅述，具體來說在[5]中2.4節可以找到非周期性的證明，而在[4]中可以找到Harris回歸性的證明。
五.參考文獻
[1]Jun Shao, Mathematical Statistics, second edition, Springer Texts in Statistics, 2003 Springer Science+Business Media, LLC.
[2]Larry Wasserman, All of Statistics: A Concise Course in Statistical Inference,
[3]D. Revuz, Markov Chains, North-Holland Mathematical Library. Vol 11, 1984, North-Holland, Amsterdam, New York, Oxford.
[4]Luke Tierney, Markov Chains for Exploring Posterior Distributions, The Annals of Statistics, Vol. 22, No. 4. (Dec., 1994), pp. 1701-1728.
[5]E.Nummelin, General Irreducible Markov Chains and non-negative operators, Cambridge Tracks in Mathematics, Cambridge University Press,1984.
[6]N.Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller (1953). Equation of state calculations by fast computing machines. J. of Chemical Physics 21, 1087–1092.
[7]W.Hastings,(1970). Monte carlo sampling methods using markov chains and their applications. Biometrika 57(1), 97–109.
[8]S.Geman. and D. Geman (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans. on Pattern Analysis and Machine Intelligence 6(6).
[9]A.Gelfand, and A. Smith (1990). Sampling-based approaches to calculating marginal densities. J. of the Am. Stat. Assoc. 85, 385–409.
[10]Kevin P. Murphy, Machine Learning A Probabilistic Perspective,2012 Massachusetts Institute of Technology