马尔科夫链蒙特卡罗法

作者: 单调不减 | 来源:发表于2019-06-22 21:50 被阅读0次

蒙特卡罗法（Monte Carlo Method）也称为统计模拟方法，是通过概率模型的随机抽样进行近似数值计算的方法。

马尔科夫链蒙特卡罗法（Markov Chain Monte Carlo，MCMC）则是以马尔科夫链为概率模型的蒙特卡罗法。MCMC构建一个马尔科夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔科夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算。

1、蒙特卡罗法

在我浅显的认知里，蒙特卡罗法就是在正方形里随机取点用概率近似求内接圆面积的近似方法。看过《统计学习方法》后的简单介绍后发现自己还是naive，蒙特卡罗法还是很强大的。

蒙特卡罗法要解决的问题是：假设概率分布定义已知，通过抽样获得概率分布的随机样本，并通过得到的随机样本对概率分布的特征进行分析。蒙特卡罗法的核心是随机抽样。

一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。后两者适用于概率密度复杂、不能直接抽样的情况。

接受-拒绝法的基本思想如下：假设 $p(x)$ 不可直接抽样，找到一个可以直接抽样的分布，称为建议分布。假设 $q(x)$ 为建议分布的概率密度函数，并且有 $cq(x)\geq p(x)$ ，其中 $c>0$ 。按照 $q(x)$ 进行抽样，假设得到结果是 $x^*$ ，再按照 $\frac{p(x^*)}{cq(x^*)}$ 的比例决定是否接受 $x^*$ 。如下图所示，落到 $p(x^*)$ 范围内就接受（绿色），否则拒绝（红色）。

具体算法流程如下：

（1）选择概率密度函数为 $q(x)$ 的概率函数作为建议分布，使其对任一 $x$ 满足 $cq(x)\geq p(x)$ ，其中 $c>0$ 。

（2）按照建议分布 $q(x)$ 随机抽样得到样本 $x^*$ ，再按照均匀分布在 $(0,1)$ 范围内抽样得到 $u$ 。

（3）若 $u\leq\frac{p(x^*)}{cq(x^*)}$ ，则将 $x*$ 作为抽样结果，否则返回（2）。

（4）直至得到 $n$ 个随机样本，结束。

接受-拒绝法实际是按照 $p(x)$ 的涵盖面积占 $cq(x)$ 的涵盖面积的比例进行抽样。其优点是容易实现，缺点是效率可能不高。若 $p(x)$ 涵盖面积占 $cq(x)$ 的涵盖面积的比例很低，就会导致拒绝比例很高，抽样效率很低。

一般蒙特卡罗法可以用于数学期望估计。按概率分布 $p(x)$ 独立地抽取 $n$ 个样本 $x_1,x_2,\dots,x_n$ ，之后计算函数 $f(x)$ 的样本均值 $\hat{f}_n$ ：

$\hat{f}_n=\frac{1}{n}\sum_{i=1}^n f(x_i)$

作为数学期望 $E_{p(x)}[f(x)]$ 的近似值。这里可以通过离散的情况简单的理解，比如 $x=x_i$ 被抽到 $m$ 次，把 $m$ 和 $\frac{1}{n}$ 相乘就得到 $x=x_i$ 的频率，可用于估计概率 $p(x_i)$ ，从而上式可视为 $\sum_x f(x)p(x)$ 即期望 $E_{p(x)}[f(x)]$ 的近似值。

此外，一般蒙特卡罗法还可以用于积分计算。假设有一个函数 $h(x)$ ，目标是计算该函数积分：

$\int_{x}h(x)dx$

将 $h(x)$ 分解成一个函数 $f(x)$ 和一个概率密度函数 $p(x)$ 的乘积，则有：

$\int_x h(x)dx=\int_x f(x)p(x)dx=E_{p(x)}[f(x)]$

于是函数 $h(x)$ 的积分就可以表示为函数 $f(x)$ 关于 $p(x)$ 的期望。实际上给定一个概率密度函数 $p(x)$ ，只要取 $f(x)=\frac{h(x)}{p(x)}$ ，即可得到上述形式。接下来，我们用期望的估计即可计算积分：

$\int_x h(x)dx=E_{p(x)}[f(x)]\approx \frac{1}{n}\sum_{i=1}^n f(x_i)$

2、马尔科夫链

马尔科夫链定义：

考虑一个随机变量的序列 $X=\{X_0,X_1,\dots,X_t,\dots\}$ ，这里 $X_t$ 表示时刻 $t$ 的随机变量， $t=0,1,2,\dots$ 。每个随机变量 $X_t(t=0,1,2,\dots)$ 的取值集合相同，称为状态空间，表示为 $S$ 。随机变量可以是离散的，也可以是连续的。以上随机变量的序列构成随机过程。

假设在时刻 $0$ 的随机变量 $X_0$ 遵循概率分布 $P(X_0)=\pi_0$ ，称为初始状态分布。在某时刻 $t\geq 1$ 的随机变量 $X_t$ 与前一时刻的随机变量 $X_{t-1}$ 有条件分布 $P(X_t|X_{t-1})$ 。若 $X_t$ 只依赖于 $X_{t-1}$ 而不依赖于过去的随机变量 $\{X_0,X_1,\dots,X_{t-2} \}$ ，这一性质称为马尔可夫性，即：

$P(X_t|X_0,X_1,\dots,X_{t-1})=P(X_t|X_{t-1}),\quad t=1,2,\dots$

具有马尔科夫性的随机序列 $X=\{X_0,X_1,\dots,X_t,\dots\}$ 称为马尔科夫链或马尔科夫过程。条件概率分布 $P(X_t|X_{t-1})$ 称为马尔科夫链的转移概率分布。

若马尔科夫链在时刻 $(t-1)$ 处于状态 $j$ ，在时刻 $t$ 移动到状态 $i$ ，将转移概率记作：

$p_{ij}=(X_t=i|X_{t-1}=j),\quad i=1,2,\dots,\quad j=1,2,\dots$

满足：

$p_{ij}\geq 0,\quad \sum_i p_{ij}=1$

马尔科夫链的转移概率 $p_{ij}$ 可用矩阵表示：

$P= \begin{bmatrix} p_{11} & p_{12} & p_{13} & \dots\\ p_{21} & p_{22} & p_{23} & \dots\\ p_{31} & p_{32} & p_{33} & \dots\\ \dots & \dots & \dots & \dots\\ \end{bmatrix} \qquad$

考虑马尔科夫链 $X=\{X_0,X_1,\dots,X_t,\dots\}$ 在时刻 $t$ 的概率分布，称为时刻 $t$ 的状态分布：

$\pi(t)=\begin{bmatrix} \pi_1(t)\\ \pi_2(t)\\ \dots\\ \dots\\ \end{bmatrix}$

其中 $\pi_i(t)$ 表示时刻 $t$ 状态为 $i$ 的概率 $P(X_t=i)$ 。

马尔科夫链 $X$ 在时刻 $t$ 的状态分布可以由在时刻 $(t-1)$ 的状态分布以及转移概率分布决定：

$\pi(t)=P\pi(t-1)$

递推得到：

$\pi(t)=P^t \pi(0)$

若状态空间 $S$ 上存在一个分布 $\pi$ 使得：

$\pi=P\pi$

则称 $\pi$ 为马尔科夫链 $X=\{X_0,X_1,\dots,X_t,\dots\}$ 的平稳分布。

下面给出几个定义及定理：

（1）不可约：对于任意状态 $i,j\in S$ ，存在时刻 $t$ ，使得时刻 $0$ 从状态 $j$ 出发，时刻 $t$ 到达状态 $i$ 的概率大于 $0$ ，则称此马尔科夫链 $X$ 是不可约的。否则称马尔科夫链是可约的。

（2）非周期：对任意状态 $i\in S$ ，若时刻 $0$ 从状态 $i$ 出发， $t$ 时刻返回状态的所有时间长 $\{ t:P(X_t=i|X_0=i)>0 \}$ 的最大公约数为1，则称此马尔科夫链 $X$ 是非周期的，否则称马尔科夫链是周期的。

（3）正常返：对于任意状态 $i,j\in S$ ，定义概率 $p_{ij}^t$ 为时刻 $0$ 从状态 $j$ 出发，时刻 $t$ 首次转移到状态 $i$ 的概率，若对任意状态 $i,j\in S$ 都满足： $\lim_{t\rightarrow \infty}p_{ij}^t>0$ ，则称马尔科夫链 $X$ 是正常返的。

（4）可逆：若有一个状态分布 $\pi=(\pi_1,\pi_2,\dots)^T$ ，对任意状态 $i,j\in S$ ，对任意时刻 $t$ ，满足：

$p_{ji}\pi_i=p_{ij}\pi_j$

则称此马尔科夫链 $X$ 为可逆马尔科夫链。上式称为细致平衡方程。

定理1：不可约且非周期的有限状态马尔科夫链有唯一平稳分布存在。

定理2：不可约、非周期且正常返的马尔科夫链有唯一平稳分布存在。

定理3（遍历定理）：不可约、非周期且正常返的马尔科夫链有唯一平稳分布 $\pi=(\pi_1,\pi_2,\dots)^T$ ，且转移概率的极限分布式马尔科夫链的平稳分布：

$\lim_{t\rightarrow \infty}P(X_t=i|X_0=j)=\pi_i,\quad i=1,2,\dots,\quad j=1,2,\dots$

若 $f(X)$ 是定义在状态空间上的函数， $E_{\pi}[|f(X)|]<\infty$ ，则：

$P \{\hat{f}_t\rightarrow E_{\pi}[|f(X)|] \}=1$

这里：

$\hat{f}_t=\frac{1}{t}\sum_{s=1}^t f(x_s)$

遍历定理是说，满足相应条件的马尔科夫链，当时间趋于无穷时，马尔科夫链的状态分布趋近于平稳分布，随机变量的函数的样本均值以概率1收敛于该函数的数学期望。

定理4：满足细致平衡方程的状态分布 $\pi$ 就是该马尔科夫链的平稳分布，即：

$P\pi=\pi$

证明：

$(P\pi)_i=\sum_j p_{ij}\pi_j=\sum_j p_{ji}\pi_i=\pi_i\sum_j p_{ji}=\pi_i$

3、马尔科夫链蒙特卡罗法

假设多元随机变量 $x$ 概率密度函数为 $p(x)$ ， $f(x)$ 为 $x$ 的函数，目标是获得概率分布 $p(x)$ 的样本集合，以及函数 $f(x)$ 的数学期望 $E_{p(x)}[f(x)]$ 。

应用马尔科夫链蒙特卡罗法解决这个问题的基本想法是：在随机变量 $x$ 的状态空间 $S$ 上定义一个满足遍历定理的马尔科夫链 $X=\{X_0,X_1,\dots,X_t,\dots \}$ ，使其平稳分布就是抽样的目标分布 $p(x)$ 。然后在这个马尔科夫链上进行随机游走，每个时刻得到一个样本。根据遍历定理，当时间趋于无穷时，样本的分布趋于平稳分布，样本的函数均值趋近于函数的数学期望。

所以，当时间足够长时（时刻大于某个正整数 $m$ ， $0\sim m$ 的时间段称为燃烧期），在之后的时间（时刻小于某个给定正整数 $n$ ， $n>m$ ）里随机游走得到的样本集合 $\{x_{m+1},x_{m+2},\dots,x_n\}$ 就是目标概率分布的抽样结果，得到的函数均值就是要计算的期望值：