MCMC-1|机器学习推导系列（十五）

一、蒙特卡洛方法

Monte Carlo Method也就是基于采样的随机近似方法。该方法旨在求得复杂概率分布下的期望值：$E_{z|x}[f(z)]=\int p(z|x)f(z)\mathrm{d}z\approx \frac{1}{N}{\sum }_{i=1}^{N}f(z_i)$，其中$z_i$是从概率分布$p(z|x)$中取的样本，也就是说从概率分布中取$N$个点，从而近似计算这个积分。

这里介绍三种采样方法：

1. 概率分布采样

首先要求得概率密度函数PDF的累计密度函数CDF，然后求CDF得反函数，在0-1之间均匀取样，代入反函数，就得到了取样点。这个方法的缺点就是大部分PDF很难求得CDF：

2. 拒绝采样（Rejection Sampling）

对于较复杂的概率分布$p(z)$，引入简单的提议分布（proposal distribution）$q(z)$，使得任意的$Mq(z_i)\ge p(z_i)$，然后对$q(z)$进行采样获得样本。具体的采样方法步骤为：
①选择概率密度函数为$q(z)$，作为提议分布，使其对任一$z$满足$Mq(z_i)\ge p(z_i)$，其中$M>0$；
②按照提议分布$q(z)$随机抽样得到样本$z_i$，再按照均匀分布在$(0,1)$范围内抽样得到$u_i$；
③如果$u_i\le \frac{p(z_i)}{Mq(z_i)}$，则将$z_i$作为抽样结果；否则，返回步骤②；
④直到获得$N$个样本，结束。

拒绝采样的优点是容易实现，缺点是采样效率可能不高。如果$p(z)$的涵盖体积占$Mq(z)$的涵盖体积的比例很低，就会导致拒绝的比例很高，抽样效率很低。注意，一般是在高维空间抽样，会遇到维度灾难的问题，即使$p(z)$与$Mq(z)$很接近，两者涵盖体积的差异也可能很大。

3. 重要性采样（Importance Sampling）

直接对期望$E_{p(z)}[f(z)]$进行采样。这里引入另一个分布$q(z)$：

$$\begin{gathered} E_{p(z)}[f(z)]=\int p(z)q(z)\mathrm{d}z \\ =\int \frac{p(z)}{q(z)}\cdot q(z)\cdot f(z)\mathrm{d}z \\ =\int f(z)\cdot \frac{p(z)}{q(z)}\cdot q(z)\mathrm{d}z \\ \approx \frac{1}{N}\sum _{i=1}^{N}f(z_i)\cdot \underset{weight}{\underbrace{\frac{p(z_i)}{q(z_i)}}} \\ {z}_i\sim q(z),i=1,2,\cdots ,N \end{gathered}$$

于是采样在$q(z)$中采样，并通过权重计算和。重要值采样对于权重⾮常⼩的时候，效率非常低。

重要性采样有⼀个变种 Sampling-Importance-Resampling，这种方法，首先和上面⼀样进行采样，然后在采样出来的N个样本中，重新采样，这个重新采样，使⽤每个样本点的权重作为概率分布进行采样。

二、马尔可夫链

1. 齐次马尔科夫链

考虑一个随机变量的序列 X = { X 0 , X 1 , ⋯   , X t , ⋯   } X=\left \{X_{0},X_{1},\cdots ,X_{t},\cdots \right \} X={X0​,X1​,⋯,Xt​,⋯}，这里的$X_t$表示$t$时刻的随机变量，每个随机变量的取值空间相同。

如果$X_t$只依赖于$X_{t-1}$，而不依赖于过去的随机变量 { X 0 , X 1 , ⋯   , X t − 2 } \left \{X_{0},X_{1},\cdots ,X_{t-2}\right \} {X0​,X1​,⋯,Xt−2​}，这一性质称为马尔可夫性，即

$$P(X_t|X_1,X_2,\cdots X_{t-1})=P(X_t|X_{t-1}),t=1,2,\cdots$$

具有马尔可夫性的随机序列 X = { X 0 , X 1 , ⋯   , X t , ⋯   } X=\left \{X_{0},X_{1},\cdots ,X_{t},\cdots \right \} X={X0​,X1​,⋯,Xt​,⋯}称为马尔可夫链（Markov Chain），或马尔可夫过程（Markov Process）。条件概率分布$P(X_t|X_{t-1})$称为马尔可夫链的转移概率分布。

当转移概率分布$P(X_t|X_{t-1})$与$t$无关，也就是说不同时刻的转移概率是相同的，则称该马尔可夫链为时间齐次的马尔可夫链（Time Homogenous Markov Chain），形式化的表达是：

$$P(X_{t+s}|X_{t-1+s})=P(X_t|X_{t-1}),t=1,2,\cdots ;s=1,2,\cdots$$

2. 转移概率矩阵和状态分布

• 转移概率矩阵

如果马尔可夫链的随机变量$X_t(t=0,1,2,\cdots )$定义在离散空间，则转移概率分布可以由矩阵表示。若马尔可夫链在时刻$t-1$处于状态$j$，在时刻$t$移动到状态$i$，将转移概率记作：

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

满足：

$$p_{ij}\ge 0,\sum _i{p}_{ij}=1$$

马尔可夫链的转移概率可以由矩阵表示：

$$\begin{gathered} P=\left[\begin{array}{cccc}{p}_{11}& p_{12}& p_{13}& \cdots \\ p_{21}& p_{22}& p_{23}& \cdots \\ p_{31}& p_{32}& p_{33}& \cdots \\ \cdots & \cdots & \cdots & \cdots \end{array}\right] \\ {p}_{ij}\ge 0,\sum _i{p}_{ij}=1 \end{gathered}$$

• 状态分布

考虑马尔可夫链 X = { X 0 , X 1 , ⋯   , X t , ⋯   } X=\left \{X_{0},X_{1},\cdots ,X_{t},\cdots \right \} X={X0​,X1​,⋯,Xt​,⋯}在时刻$X_t(t=0,1,2,\cdots )$的概率分布，称为时刻$t$的状态分布，记作：

$$\pi (t)=\left[\begin{array}{c}{\pi }_1(t)\\ {\pi }_2(t)\\ ⋮\end{array}\right]$$

其中${\pi }_i(t)$表示时刻$t$状态为$i$的概率$P(X_i=i)$，即：

$${\pi }_i(t)=P(X_i=i),i=1,2,\cdots$$

对于马尔可夫链的初始状态分布：

$$\pi (0)=\left[\begin{array}{c}{\pi }_1(0)\\ {\pi }_2(0)\\ ⋮\end{array}\right]$$

其中${\pi }_i(0)$表示时刻$0$状态为$i$的概率$P(X_0=i)$。通常初始分布${\pi }_i(0)$的向量只有一个分量是$1$，其余分量为$0$，表示马尔可夫链是从一个具体状态开始的。

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

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

其中：

$$\begin{gathered} {\pi }_i(t)=P(X_t=i) \\ =\sum _{m}P(X_t=i|X_{t-1}=m)P(X_{t-1}=m) \\ =\sum _m{p}_{im}{\pi }_m(t-1) \end{gathered}$$

马尔可夫链在时刻$t$的状态分布，可以通过递推得到：

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

递推得到：

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

这里的递推式表明马尔可夫链的状态分布由初始分布和转移概率分布决定。

这里的$P^t$称为$t$步转移概率矩阵：

$$P_{ij}^t=P(X_t=i|X_0=j)$$

3. 平稳分布

• 定义

对于一个马尔可夫链$X$，如果其状态空间上存在一个分布：

$$\pi =\left[\begin{array}{c}{\pi }_1\\ {\pi }_2\\ ⋮\end{array}\right]$$

使得：

$$\pi =P\pi$$

则称$\pi$为马尔可夫链$X$的平稳分布。

直观地来看，如果以该平稳分布作为初始分布，面向未来进行随机状态转移，之后任何一个时刻的状态分布都是该平稳分布。

• 定理

给定一个马尔可夫链$X$，状态分布为$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$为$X$的充要条件
为$\pi$是下列方程组的解：

$$\{\begin{array}{c}①x_i={\sum }_j{p}_{ij}{x}_j,i=1,2,\cdots \\ ②x_i\ge 0,i=1,2,\cdots \\ ③{\sum }_i{x}_i=1\end{array}$$

证明如下：

必要性证明
假设$\pi$是平稳分布，显然满足①和②，又因为${\pi }_i={\sum }_j{p}_{ij}{\pi }_j$，所以满足③。
充分性证明
由②和③知$\pi$是一概率分布。假设$\pi$是$X_t$的分布，则有：
$$P(X_t=i)={\pi }_i=\sum _j{p}_{ij}P(X_{t-1}=j)$$
又因为$\pi$满足①，所以有：
$$P(X_t=i)={\pi }_i=\sum _j{p}_{ij}{\pi }_j$$
综合两式则可得出$P(X_{t-1}=j)={\pi }_j,j=1,2,\cdots$，即$\pi$也是$X_{t-1}$的概率分布。事实上这对任意$t$都成立，所以$\pi$是马尔可夫链的平稳分布。

这个定理给出了一个求马尔可夫链平稳分布的方法。

4. 连续状态马尔可夫链

• 概率转移核

连续状态马尔可夫链的转移概率分布由概率转移核或转移核（transition kernel）表示。

在连续状态空间$S$上，对任意的$x\in S,A\subset S$（$A$可以理解为一个区间），转移核$P(x,A)$定义为：

$$P(x,A)={\int }_{A}p(x,y)\mathrm{d}y$$

其中$p(x,\cdot )$为概率密度函数，满足$p(x,\cdot )\ge 0,P(x,S)={\int }_{S}p(x,y)\mathrm{d}y=1$。转移核$P(x,A)$表示从$x\sim A$的转移概率：

$$P(X_t=A|X_{t-1}=x)=P(x,A)$$

有时也将概率密度函数$p(x,\cdot )$称为转移核。

• 平稳分布

若马尔可夫链的状态空间$S$上的概率分布$\pi (x)$满足条件：

$$\pi (y)=\int p(x,y)\pi (x)\mathrm{d}x,\forall y\in S$$

则$\pi (x)$为该马尔可夫链的平稳分布。等价地：

$$\pi (A)=\int p(x,A)\pi (x)\mathrm{d}x,\forall A\in S$$

或简写为：

$$\pi =P\pi$$

三、马尔可夫链的性质

以下通过离散状态马尔可夫链介绍马尔可夫链的性质，可以推广到连续状态马尔可夫链。

1. 不可约

在状态空间$S$中对于任意状态$i,j\in S$，如果存在一个时刻$t(t>0)$满足：

$$P(X_t=i|X_0=j)>0$$

也就是说，时刻$0$从状态$j$出发，时刻$t$到达状态$i$的概率大于$0$，则称此马尔可夫链是不可约的（irreducible），否则称马尔可夫链是可约的（reducible）。

直观上，一个不可约的马尔可夫链，从任意状态出发，当经过充分长时间后，可以到达任意状态。

举例：

不可约：

可约：

2. 非周期

在状态空间$S$中对于任意状态$i\in S$，如果时刻$0$从状态$i$出发，$t$时刻返回状态的所有时间长 { t : P ( X t = i ∣ X 0 = i ) > 0 } \left \{t:P(X_{t}=i|X_{0}=i)> 0\right \} {t:P(Xt​=i∣X0​=i)>0}的最大公约数是$1$，则称此马尔可夫链是非周期的（aperiodic），否则称马尔可夫链是周期的（periodic）。

直观上，一个非周期性的马尔可夫链，不存在一个状态，从这一个状态出发，再返回到这个状态时所经历的时间长呈一定的周期性，也就是说非周期性的马尔可夫链的任何状态都不具有周期性。

举例：

非周期：

周期：

3. 正常返

对于任意状态$i,j\in S$，定义概率$p_{ij}^t$为时刻$0$从状态$j$出发，时刻$t$首次转移到状态$i$的概率，即$p_{ij}^t=P(X_t=i,X_s\ne i,s=1,2,\cdots ,t-1|X_0=j),t=1,2,\cdots$。若对所有状态$i,j$都满足${\lim }_{t\to \infty }{p}_{ij}^t>0$，则称马尔可夫链是正常返的（positive recurrent）。

直观上，一个正常返的马尔可夫链，其中任意一个状态，从其他任意一个状态出发，当时间趋于无穷时，首次转移到这个状态的概率不为$0$。

定理：

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

4. 遍历定理

若马尔可夫链是不可约、非周期且正常返的，则该马尔可夫链有唯一平稳分布$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$，并且转移概率的极限分布是马尔可夫链的平稳分布：

$$\underset{t\to \infty }{\lim }P(X_t=i|X_0=j)={\pi }_i,i=1,2,\cdots ,j=1,2,\cdots$$

也就是：

$$\underset{t\to \infty }{\lim }P(X_t=i)={\pi }_i,i=1,2,\cdots$$

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

$$P\left\{{\hat{f}}_t\to E_{\pi }[f(X)]\right\}=1$$

这里：

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

$E_{\pi }[f(X)]={\sum }_{i}f(i){\pi }_i$是$f(X)$关于平稳分布$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$的数学期望，$P\left\{{\hat{f}}_t\to E_{\pi }[f(X)]\right\}=1$表示${\hat{f}}_t\to E_{\pi }[f(X)],t\to \infty$几乎处处成立或以概率$1$成立。

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

样本均值可以认为是时间均值，数学期望是空间均值。遍历定理表述了遍历性的含义：当时间趋于无穷时，时间均值等于空间均值。

遍历定理的三个条件：不可约、非周期、正常返，保证了当时间趋于无穷时达到任意一个状态的概率不为$0$。

理论上并不知道经过多少次迭代，马尔可夫链的状态分布才能接近平稳分布，在实际应用遍历定理时，取一个足够大的整数$m$，经过$m$次迭代之后认为状态分布就是平稳分布，这时计算从第$m+1$次迭代到第$n$次迭代的均值，即：

$$E_{\pi }[f(X)]=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

称为遍历均值。

5. 可逆马尔可夫链

• 定义

对于任意状态$i,j\in S$，对任意一个时刻$t$满足：

$$P(X_t=i|X_{t-1}=j){\pi }_j=P(X_{t-1}=j|X_t=i){\pi }_i,i,j=1,2,\cdots$$

或简写为：

$$p_{ij}{\pi }_j=p_{ji}{\pi }_i,i,j=1,2,\cdots$$

则称此马尔可夫链为可逆马尔可夫链（reversible Markov chain），上式又被称作细致平衡方程（detailed balance equation）。

直观上，如果有可逆马尔可夫链，那么以该马尔可夫链的平稳分布作为初始分布，进行随机状态转移，无论是面向过去还是面向未来，任何一个时刻的状态分布都是该平稳分布。

• 定理

满足细致平衡方程的状态分布$\pi$就是该马尔可夫链的平稳分布，即满足$P\pi =\pi 。$

证明：

$$(P\pi {)}_i=\sum _j{p}_{ij}{\pi }_j=\sum _j{p}_{ji}{\pi }_i={\pi }_i\underset{=1}{\underbrace{\sum _j{p}_{ji}}}={\pi }_i,i=1,2,\cdots$$

该定理说明，可逆马尔可夫链一定有唯一平稳分布，给出了一个马尔可夫链有平稳分布的充分条件（不是必要条件）。也就是说，可逆马尔可夫链满足遍历定理的条件。

参考资料

ref:李航《统计学习方法》