一、概述

1、动机

MCMC是一种随机的近似推断，其核心就是基于采样的随机近似方法蒙特卡洛方法。对于采样的动机而言，有这样的使用需求：

1. 采样作为任务，用于生成新的样本本身就是机器学习的常见任务
2. 求和/求积分

采样结束后，我们需要评价采样出来的样本点是不是好的样本集：

1. 样本趋向于高概率的区域
2. 样本之间必须独立

具体采样中，采样是一个困难的过程：

1. 无法采样得到归一化因子，即无法直接对概率 $p(x)=\frac{1}{Z}\hat{p}(x)$ 采样，常常需要对 CDF 采样，但复杂的情况不行
2. 如果归一化因子可以求得，但是对高维数据依然不能均匀采样（维度灾难），这是由于对 $p$ 维空间，总的状态空间是 $K^p$ 这么大，于是在这种情况下，直接采样也不行

二、平稳分布

定义状态转移矩阵（随机矩阵）

$$Q={\left(\begin{array}{cccc}{Q}_{11}& Q_{12}& \cdots & Q_{1K}\\ ⋮& ⋮& ⋮& ⋮\\ Q_{k1}& Q_{k2}& \cdots & Q_{KK}\end{array}\right)}_{k\times k}$$

这个矩阵每一行代表这一个时刻$q^{(i)}x$的K(状态空间）种取值可能概率，一行之和为1,不过到底是行等于1还是列等于1，得看矩阵中每个元素即转移概率的定义。是从行上的状态转移到列上的状态，还是从列上的状态转移到行上的状态。此外，随机矩阵的特征值 ≤ 1。假设只有一个特征值为 ${\lambda }_i=1$。于是在马尔可夫过程中：

$$\begin{gathered} q^{t+1}(x=j)=\sum _{i=1}^K{q}^t(x=i)Q_{ij} \\ 这种转换的目的是将状态用矩阵表示 \end{gathered}$$

具体过程如下：

$$\begin{array}{rl}令q^{(t+1)}& ={\left(q^{(t+1}(x=1)q^{(t+1)}(x=2)\cdots q^{(t+1)}(x=k)\right)}_{1\times k}\\ \text{ 而 }& q^{(t+1)}(x=j)=\sum _{i=1}^k{q}^{(t)}(x=i)\cdot Q_{ij}\\ \therefore q^{(t+1)}& ={\left(\sum _{i=1}^k{q}^{(t)}(x=i)Q_{i1}\sum _{i=1}^k{q}^{(t)}(x=i)\cdot Q_{i2}\cdots \sum _{i=1}^k{q}^{(t)}(x=i)\cdot Q_{ik}\right)}_{i\times k}\\ & =q^{(t)}\cdot Q\\ \iff & q^{(t)}=\left(q^{(t)}(x=1)q^{(t)}(x=2)\cdots q^{(t)}(x=k)\right)\end{array}$$

$$\begin{gathered} 该过程可持续迭代 \\ \Rightarrow q^{(t+1)}=q^{(t)}\cdot Q=q^{(1)}{Q}^t \end{gathered}$$

对于Q而言：

∴Q=A⋅Λ⋅A−1Λ=(λ1…λk)⋅∣λi∣⩽1, for i=1,2,,…,K则:q(t+1)=q(1)⋅(A⋅∧⋅A−1)t=q(1)⋅A∧tA−1存在足够大的m,s.t.Λm=(01…0)由于随机矩阵特征值的绝对值均≤1，因此可以设λi为1(简化运算)q(m+1)=q(1)⋅A⋅ΛmA−1q(m+2)=q(m+1)⋅A⋅Λ⋅A−1=q(1)⋅A∧mA−1A∧A−1=q(1)∧m+1A−1=q(m+1)∴q(m+2)=q(m+1)当t>m时,q(m+1)=q(m+2)=⋯=q(∞)=⋯begin{aligned}
& therefore Q=A cdot Lambda cdot A^{-1}
&
Lambda=left(begin{array}{lll}
lambda_{1} &
& … &
& & lambda_k
end{array}right) cdotleft|lambda_iright| leqslant 1, text { for } i=1,2,,…,K
& 则: q^{(t+1)}=q^{(1)} cdotleft(A cdot wedge cdot A^{-1}right)^t=q^{(1)} cdot A wedge^t A^{-1}
& 存在足够大的m, s.t. Lambda^m = left(begin{array}{lll}
0 & &&
& 1 & &
& & …&
& & &0

end{array}right)

& 由于随机矩阵特征值的绝对值均≤1，因此可以设lambda_i为1(简化运算)
& q^{(m+1)}=q^{(1)} cdot A cdot Lambda^m A^{-1}
& q^{(m+2)}=q^{(m+1)} cdot A cdot Lambda cdot A^{-1}
& =q^{(1)} cdot A wedge^{m} A^{-1} A wedge A^{-1}
& =q^{(1)} wedge^{m+1} A^{-1}=q^{(m+1)}
& therefore q^{(m+2)}=q^{(m+1)}
& begin{array}{l}
当t > m时,
q^{(m+1)}=q^{(m+2)}=cdots=q^{(∞)}=cdots
end{array}
&
end{aligned}​∴Q=A⋅Λ⋅A−1Λ=⎝⎛​λ1​​…​λk​​⎠⎞​⋅∣λi​∣⩽1, for i=1,2,,…,K则:q(t+1)=q(1)⋅(A⋅∧⋅A−1)t=q(1)⋅A∧tA−1存在足够大的m,s.t.Λm=⎝⎛​0​1​…​0​⎠⎞​由于随机矩阵特征值的绝对值均≤1，因此可以设λi​为1(简化运算)q(m+1)=q(1)⋅A⋅ΛmA−1q(m+2)=q(m+1)⋅A⋅Λ⋅A−1=q(1)⋅A∧mA−1A∧A−1=q(1)∧m+1A−1=q(m+1)∴q(m+2)=q(m+1)当t>m时,q(m+1)=q(m+2)=⋯=q(∞)=⋯​​

正因为马尔科夫链有这样一个能够收敛的性质，因此我们可以联想到构造一条马尔科夫链，让最后收敛的平稳分布等于或者逼近我们想要的那种分布，后可从已经收敛的状态中，选择一个状态取出需要数量的样本，达到合理采样的目的

三、蒙特卡洛方法

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\left(z_i\right)\ge p\left(z_i\right)$ ，然后对 $q(z)$ 进行采样获得样本。具体的采样方法步骤为:

①选择概率密度函数为 $q(z)$ ，作为提议分布，使其对任一 $z_i$ 满足 $Mq\left(z_i\right)\ge p\left(z_i\right)$ ，其中 $M>0$ ；

②按照提议分布 $q(z)$ 随机抽样得到样本 $z_i$ ，再按照均匀分布在 $(0,1)$ 范围内抽样得到 $u_i$ ；

③如果 $u_i\le \frac{p\left(z_i\right)}{Mq\left(z_i\right)}$ ，则将 $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 f(z)\cdot p(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}{\frac{p(z_i)}{q(z_i)}} \\ {z}_i\sim q(z),i=1,2,\cdots ,N \end{gathered}$$

重要性采样并不像接受-拒绝采样一样是一种采样方法，只是一种让估计更加精确的积分策略。接受拒绝采样是一种生成我们想要的分布的方法，而重要性采样是利用某种分布（这种分布能区分出被积函数不同点的重要性）更好地积分。本质上二者做的是两件不同的事情。

例如，如果分布相当不同，采用这种方法采样是非常不合理的，越差的拟合效果反而有更大的权重了。

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

4. Sampling-Importance-Resampling

重要性采样有⼀个变种

这种方法，首先和上面⼀样进行采样，然后在采样出来的N个样本中，重新采样，这个重新采样，使⽤每个样本点的权重作为概率分布进行采样。在后续章节粒子滤波会详细介绍。

四、马尔可夫链的类型

1. 齐次(一阶)马尔科夫链

考虑一个随机变量的序列$X=\left\{X_0,X_1,\cdots ,X_t,\cdots \right\}$，这里的$X_t$表示$t$时刻的随机变量，每个随机变量的取值空间相同。

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

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

具有马尔可夫性的随机序列$X=\left\{X_0,X_1,\cdots ,X_t,\cdots \right\}$称为马尔可夫链（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. 转移概率矩阵和状态分布

PS: 以下内容结合这个图会容易理解一点，注意状态分布的定义和大佬的不一样

• 转移概率矩阵

如果马尔可夫链的随机变量$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$$

满足每个行的和为$1$（具体是行还是列根据自己的定义）：

$$p_{ij}\ge 0,\sum _i{p}_{ij}=1,p_{ij}=p(X_{t+1}=j|X_t=i)$$

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

$$\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=\left\{X_0,X_1,\cdots ,X_t,\cdots \right\}$在时刻$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_t=i)$，即：

$${\pi }_i(t)=P(X_t=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.Detail Balance

那么我们如何构造一个马氏链，使得它能收敛到平稳分布呢？Detail Balance，他是平稳分布的充分非必要条件（Detail Balance能推导到平稳分布，平稳分布不能推导到它）

Detail Balance：

$$p\left(x\mapsto x^{\ast }\right)\cdot \pi (x)=P\left(x^{\ast }\mapsto x\right)\cdot \pi \left(x^{\ast }\right)$$

其中：$x\mapsto x^{\ast }为x^{\ast }|x$

$$\begin{array}{rl}& \int p\left(x\mapsto x^{\ast }\right)\cdot \pi (x)dx\\ & =\int P\left(x^{\ast }\mapsto x\right)\cdot \pi \left(x^{\ast }\right)dx\\ & =\pi \left(x^{\ast }\right)\cdot \underset{1}{\int P\left(x^{\ast }\mapsto x\right)dx}\\ & =\pi \left(x^{\ast }\right)\end{array}$$

五、马尔可夫链的性质

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

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$时刻返回状态的所有时间长$\left\{t:P(X_t=i|X_0=i)>0\right\}$的最大公约数是$1$，则称此马尔可夫链是非周期的（aperiodic），否则称马尔可夫链是周期的（periodic)。

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

举例：

非周期：

周期：

3. 正常返

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

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

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

4.遍历定理
设有马尔可夫链 $X=\left\{X_0,X_1,\cdots ,X_t,\cdots \right\}$ ，若这个马尔可夫链是不可约、非周期且正常返的，则该马尔可夫链有唯一平稳分布 $\pi ={\left({\pi }_1,{\pi }_2,\cdots \right)}^T$ ，并且转移概率的极限分布是马尔可夫链的平稳分布:

$$\underset{t\to \infty }{\lim }P\left(X_t=i\mid X_0=j\right)={\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\left(x_s\right)$$

$E_{\pi }[f(X)]={\sum }_{i}f(i){\pi }_i$ 是 $f(X)$ 关于平稳分布 $\pi ={\left({\pi }_1,{\pi }_2,\cdots \right)}^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 成立。

遍历定理的直观解释: 满足相应条件的马尔可夫链，当时间趋于无穷时，马尔可夫链的状态分布趋近于平稳分布，随机变量的函数的样本均值 ${\hat{f}}_t$ 以概率 1 收敛于该函数的数学期望 $E_{\pi }[f(X)]$ 。

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

遍历定理的三个条件：不可约、非周期、正常返，保证了当时间趋于无穷时达到任意一个状态的概率不为 0 。
理论上并不知道经过多少次迭代，马尔可夫链的状态分布才能接近平稳分布，在实际应用遍历定理时，取一个足够大的整数 $m$ ，经过 $m$ 次迭代之后认为状态分布就是平稳分布，这时计算从第 $m+1$ 次迭代到第 $n$ 次迭代的均值，即:

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

称为遍历均值。

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}{\sum _j{p}_{ji}}={\pi }_i,i=1,2,\cdots$$

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

“开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 11 天，点击查看活动详情”