马尔科夫链蒙特卡罗⽅法（MCMC）
⼀.蒙特卡罗法的缺陷
通常的蒙特卡罗⽅法可以模拟⽣成满⾜某个分布的随机向量，但是蒙特卡罗⽅法的缺陷就是难以对⾼维分布进⾏模拟。对于⾼维分布的模拟，最受欢迎的算法当属马尔科夫链蒙特卡罗算法(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,
opcns
现在我们来看⼀下，这样模拟后得到的转移分布函数：
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
Loading [MathJax]/jax/output/HTML-CSS/fonts/TeX/fontdata.js
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议