当前位置：首页 > news >正文

机器学习 [白板推导]（十）[马尔可夫链蒙特卡洛法]

news 2025/8/13 8:29:27

13. 马尔可夫链蒙特卡洛法（Markov Chain & Monte Carlo，MCMC）

13.1. 蒙特卡洛方法简介

在推断中的近似推断问题中，常常会使用随机近似推断方法，例如最经典的MCMC，核心思想是蒙特卡洛方法，即基于采样的随机近似方法。

例如，要计算推断中的期望 $Ez[f(z)]E_{z}\left[f(z) \right ]$ ，假设随机变量 $z$ 是连续/无限取值的离散，则对其求期望常常是困难的（连续-积分/无限离散-无穷级数），此时可以随机采样 $N$ 个样本 $,zN}\left \{ z_1,z_2,\cdots,z_N \right \}$ ，则期望的估计值为 $,f(zN)}E_{z}\left[f(z) \right ]=\text{mean}\left \{ f(z_1),f(z_2),\cdots,f(z_N) \right \}$ 。

蒙特卡洛方法中常常使用以下采样方法：

概率分布采样：利用变量的概率分布 $p (z)$ 进行随机数生成，但这个方法仅限概率分布简单的情况；
Adjection Sampling （接收-拒绝采样法）：
- 当随机变量的概率分布 $p (z)$ 较为复杂时，可以设一个简单分布 $q (z)$ ，以及寻找一个常数 $M$ ，使得定义域上的任意 $M⋅q(z)⩾p(z)\forall \ z_i,\ \ M\cdot q(z)\geqslant p(z)$ 。
- 引入一个新的参数 $α(z)=p(z)M⋅q(z)\alpha(z)=\frac{p(z)}{M\cdot q(z)}$ ，称为接收率，其值在0-1之间，且越接近1表示 $p (z)$ 和 $M⋅q(z)M\cdot q(z)$ 越接近；
- 此时可以从简单分布 $q (z)$ 中进行采样，每采样一个 $z_i$ 的同时也从均匀分布 $U (0, 1)$ 中采样一个 $u_i$ ，当 $ui⩽α(zi)u_i\leqslant \alpha(z_i)$ 时，接受采样结果，否则拒绝，需要重新采样；
- 这种采样方法又叫接受-拒绝法，其优点是容易实现，缺点是效率低，如果 $p (z)$ 所含面积占 $M⋅q(z)M\cdot q(z)$ 所含面积比例过小，则采样很容易被拒绝，需要重新采样，类似的，如果待采样变量维度过高，也很容易采样被拒绝；
重要度采样：同样设一个简单分布 $q (z)$ ，则原本难算的期望 $Ez[f(x)]=∫zp(z)f(z)dz=∫zp(z)q(z)f(z)⋅q(z)dz=Eq(z)[p(z)q(z)f(z)]E_z\left[f(x) \right ]=\int_zp(z)f(z)dz=\int_z\frac{p(z)}{q(z)}f(z)\cdot q(z)dz=E_{q(z)}\left[\frac{p(z)}{q(z)}f(z) \right ]$ 可以化为另外一个期望，使得原本从 $p (z)$ 中采样的困难转为从简单分布 $q (z)$ 中采样，而 $p(z)q(z)\frac{p(z)}{q(z)}$ 称为重要度；

13.2. 马尔科夫链简介

马尔可夫链： $,xt−1)p(x_t)=p(x|x_1, x_2, \cdots, x_{t-1})$ ，其中 $x_t$ 表示序列上某个步/时间的样本，若一阶马尔可夫，则 $p(x_t)=p(x|x_{t-1})$ ，序列上所有随机变量的取值集合相同，被称为状态空间 $S$ ， $p(x_t)=p(x|x_{t-1})$ 被称为 t 时刻的状态转移函数。

这里只讨论一阶马尔可夫链，因为更高阶的马尔可夫链可以推导转换为一阶马尔科夫链，同时若任意的 $s$ ，有 $P(x_{t+s}|x_{t+s-1})=P(x_t|x_{t-1})$ ，则称为时间齐次的马尔可夫链，这里也只讨论时间齐次的马尔科夫链。

对于状态空间 $S$ 全部是离散值的马尔可夫链，其所有的状态转移函数可以写作一个矩阵，即 $}\left \{ \ p_{ij}=P(X_t=i|X_{t-1}=j) \ \right \}$ ，称为状态转移矩阵，满足 $∑ipij=1\sum_{i}p_{ij}=1$ 。

基于状态转移函数，马尔可夫链可以表示为，对 $t$ 时刻，取值为 $i$ 的概率为 $P(xt=i)=∑jpij⋅P(xt−1=j)=∑jpij⋅∑kpjkP(xt−2=k)=⋯P(x_t=i)=\sum_{j}p_{ij}\cdot P(x_{t-1}=j)=\sum_{j}p_{ij}\cdot \sum_kp_{jk} P(x_{t-2}=k)=\cdots$ ，为了方便表示，设矩阵 $ΠS×T={πi(t)}=P(xt=i)\Pi_{S\times T}=\{ \pi_i(t) \}=P(x_t=i)$ ，即为将每一步的计算值存储下来便于后面的计算，其中 $]T\vec{\pi}(0) =\left [ \pi_{S_0}(0), \pi_{S_1}(0), \cdots \right ]^T$ 被称为初始状态分布，常常是个独热向量，表示马尔可夫链从某个固定状态开始。

如果马尔可夫链的状态空间上存在某个分布 $]T\vec{\pi}=\left [ \pi_{S_0}, \pi_{S_1}, \cdots \right ]^T$ ，使得 $π=Pπ\pi=P\pi$ ，则称该分布为平稳分布，如果马尔可夫链进入平稳分布，则其接下来的状态分布较为平稳，变化很小，同时一个马尔可夫链可能存在唯一的平稳分布，无穷多个平稳分步，或不存在平稳分布。

马尔可夫链的性质：

可约性：若对于任意初始状态 $x_0=i$ 和任意状态 $j$ ，都至少存在一个时刻 $t > 0$ ，使得 $P(x_t=j|x_0=i)>0$ ，则说明任意初始状态经过一定时间可以到达任意状态，这种情况称马尔可夫链式不可约的，否则，说明从某个状态出发会永远到不了其他某些状态，这个马尔可夫链称为可约的。
周期性：对任意初始状态 $x_0=i$ ，若其经过一定时间后可以回到状态 $x_t=i$ ，且所有最短途径时长 ${ti}\left \{ t_i \right \}$ 的最大公约数不为1，则称该马尔科夫链具有周期性，否则是非周期的。
正常返：定义 $p_{ij}^t$ 表示 $P(xt=i,xt−m≠i∣x0=j)P(x_t=i,x_{t-m}\neq i|x_0=j)$ ，即从状态 $j$ 出发，经过 $t$ 时间首次来到状态 $i$ 的概率，若对任意 $j$ 和 $i$ ，都有 $lim⁡t→∞pijt>0\lim_{t \to \infty }p_{ij}^t>0$ ，则称为正常返的马尔科夫链，否则称为不是正常返。
可逆性：对于一个马尔可夫链，若存在一种状态分布 $]\vec{\pi}=[\pi_0, \pi_1, \cdots]$ ，使得任意两个状态 $i, j$ 都有 $P(xt=i∣xt−1=j)⋅P(xt−1=j)=P(xt=j∣xt−1=i)⋅P(xt−1=i)P(x_t=i|x_{t-1}=j)\cdot P(x_{t-1}=j)=P(x_t=j|x_{t-1}=i)\cdot P(x_{t-1}=i)$ ，即 $pijπj=pjiπip_{ij}\pi_j=p_{ji}\pi_i$ （该式被称作细致平衡方程，detailed balance），则称该马尔可夫链是可逆的。

马尔可夫链性质的定理：

不可约且非周期的有限状态马尔可夫链一定有唯一平稳分布存在；
不可约且非周期，正常返的马尔可夫链一定有唯一平稳分布存在；
遍历定理：对于一个不可约，非周期，正常返的马尔可夫链，设其唯一的平稳分布为 $π⃗\vec{\pi}$ ，则对于任意初始状态 $x_0=j$ ，有 $lim⁡t→∞P(xt=i∣x0=j)=πi\lim_{t \to \infty } P(x_t=i|x_0=j)=\pi_i$ ，也就是无论初始状态，未来经过无限时刻，状态分布都会趋近于平稳分布，同时对于一个定义在状态空间上的函数 $f(x),x∈Sf(x),x\in S$ ，则 $[f(x)]=∑if(Si)πilim_{t \to \infty}\frac{1}{t}\sum_{k=1}^tf(x_k) \to E_{\vec{\pi}}\ [f(x)]=\sum _if(S_i)\pi_i$ 。
满足细致平稳方程的分布就是该马尔可夫链的平稳分布。

13.3. MCMC算法

MCMC适用于随机变量高维且相互不独立，密度函数复杂的情况，在这样的分布中求某个函数的期望 $Ep(x)[f(x)]E_{p(x)}\left [f(x) \right ]$ ，主要思想是在随机变量的状态空间（所有取值）上定义一个满足便利定理的马尔可夫链，则其平稳分布下函数值的加权平均即为函数的数学期望。

因此MCMC计算期望的方法为 $E^[f]=1n−m∑i=m+1nf(xi)\hat{E}[f]=\frac{1}{n-m}\sum_{i=m+1}^{n}f(x_i)$ ，其中 $m$ 是一个较大的数，使得 $t > m$ 时的 $x_t$ 分布趋于平稳分布，在 $m$ 之前的遍历称为燃烧期。

因为马尔科夫链满足遍历定理，所以随机游走的起点不影响最终结果，因为都会收敛到平稳分布，这个方法比接收-拒绝采样法效率高；

13.3.1. Metropolis-Hastings算法

在MCMC算法中，第一个难点是定义一个马尔可夫链用于随机游走，使其满足遍历定理，无论从任何起点出发都将收敛到平稳分布，对于离散分布，要定义状态转移矩阵 ${pij=P(xt=i∣xt−1=j)}\left \{ p_{ij}=P(x_t=i|x_{t-1}=j) \right \}$ ，对于连续分布，要定义状态转移核函数 $p(i,j)=P(x_t=i|x_{t-1}=j)$ 。

在MH算法中，可以定义一个建议分布 $q (i, j)$ 较为简单易于操作，以及一个接收分布 $α(i,j)=min⁡{1,p(j)q(i,j)p(i)q(j,i)}\alpha (i,j)=\min \{ 1,\frac{p(j)q(i,j)}{p(i)q(j,i)} \}$ ，则转移核 $p(i,j)=q(i,j)α(i,j)={q(i,j)p(j)q(i,j)⩾p(i)q(j,i)q(j,i)⋅p(i)p(j)p(j)q(i,j)<p(i)q(j,i)p(i,j)=q(i,j)\alpha (i,j)=\left\{\begin{matrix} q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i)\cdot \frac{p(i)}{p(j)} & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.$ 。

通俗理解，MH算法的流程为，定义一个简易分布 $q (i, j)$ ，作为转移函数进行随机游走，每次从状态 $x_{t-1}=j$ 转移到状态 $i$ 时，计算接受分布 $α(i,j)=min⁡{1,p(j)q(i,j)p(i)q(j,i)}\alpha (i,j)=\min \{ 1,\frac{p(j)q(i,j)}{p(i)q(j,i)} \}$ ，并在0-1之间采样随机数 $u$ ，若 $u<α(i,j)u<\alpha (i,j)$ 则状态转移成功， $x_{t}=i$ ，否则停留在原状态 $x_{t}=j$ 。

由：
$\begin{aligned} p(j)p(i,j)&=p(j)\left\{\begin{matrix} q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i)\cdot \frac{p(i)}{p(j)} & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=\left\{\begin{matrix} p(j)q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ p(i) q(j,i) & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=p(i)\left\{\begin{matrix} q(i,j)\cdot \frac{p(j)}{p(i)} & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i) & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=p(i)p(j,i) \end{aligned}$
满足细致平衡方程，因此上面所定义的转移核函数是可逆的，即存在某个分布 $p (x)$ 使得 $p (i) p (j, i) = p (j) p (i, j)$ ，同时该分布为平稳分布。

Metropolis选择：将建议分布设为可逆分布，即 $q (i, j) = q (j, i)$ ，此时接收分布简化为 $α(i,j)=min⁡{1,p(j)p(i)}\alpha(i,j)=\min\{ 1,\frac{p(j)}{p(i)} \}$ ，这是MH算法最初采用的方法，将 $q (i, j)$ 定义为一个以 $j$ 为均值，方差为常数矩阵的正态分布，或是其他在 $j$ 附近概率较高的简单分布。

单分量Metropolis-Hastings算法：当多元变量分布抽样困难时，可以逐步对变量的每一个维度的分量进行抽样，此时建议分布变为 $q(x_t^{(m)}|x_t^{(-m)},x_{t-1}^{(m)})$ ，其中 $x_t^{(m)}$ 为t时刻状态的第 $m$ 个分量， $,xt−1(K)}x_t^{(-m)}=\{ x_t^{(1)}, x_t^{(2)}, \cdots, x_t^{(m-1)}, x_{t-1}^{(m+1)}, \cdots, x_{t-1}^{(K)} \}$ ，类似于每次固定 $K - 1$ 个维度，更新一个维度的过程，同时接受分布 $α(xt(m)∣xt−1(m),xt(−m))=min⁡{1,p(xt(m)∣xt(−m))q(xt−1(m)∣xt(m),xt(−m))p(xt−1(m)∣xt(−m))q(xt(m)∣xt−1(m),xt(−m))}\alpha (x_t^{(m)}|x_{t-1}^{(m)}, x_t^{(-m)})=\min\{ 1,\frac{p(x_{t}^{(m)}|x_{t}^{(-m)})q(x_{t-1}^{(m)}|x_{t}^{(m)}, x_{t}^{(-m)})}{p(x_{t-1}^{(m)}|x_{t}^{(-m)})q(x_t^{(m)}|x_{t-1}^{(m)}, x_{t}^{(-m)})} \}$ 。

13.3.2. 吉布斯抽样

满条件分布：对于一个多元概率联合分布 $,x(K)]Tp(x),x=[x^{(1)}, x^{(2)}, \cdots, x^{(K)}]^T$ ，定义一个条件概率分布 $p(x^I|x^{-I})$ ，其中 $I$ 是变量维度的一个子集，即 $xI={x(m),m∈I},x−I={x(m),m∉I}x^I=\{ x^{(m)},m\in I \}, x^{-I}=\{ x^{(m)},m\notin I \}$ ，也就是所有维度都在括号里出现，无论在 | 之前还是之后，这个条件概率分布叫做满条件分布，满条件分布满足性质： $p(xt−1I∣(xt−1−I)p(xtI∣(xt−I)=p(xt−1)p(xt)\frac{p(x^I_{t-1}|(x^{-I}_{t-1})}{p(x^I_{t}|(x^{-I}_{t})}=\frac{p(x_{t-1})}{p(x_{t})}$ ，左边计算大大方便，因此可以用满条件分布简化一些计算。