12 马尔可夫链蒙特卡洛

MCMC 是一种随机的近似推断，其核心就是基于采样的随机近似方法蒙特卡洛方法。对于采样任务来说，有下面一些常用的场景：

采样作为任务，用于生成新的样本
求和/求积分

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

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

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

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

因此需要借助其他手段，如蒙特卡洛方法中的拒绝采样，重要性采样和 MCMC。

12.1 蒙特卡洛方法

蒙特卡洛方法旨在求得复杂概率分布下的期望值：$\mathbb{E}_{z|x}[f(z)]=\int p(z|x)f(z)dz\simeq\frac{1}{N}\sum\limits_{i=1}^Nf(z_i)$，也就是说，从概率分布中取 $N$ 个点，从而近似计算这个积分。采样方法有：

概率分布采样，首先求得概率密度的累积密度函数 CDF，然后求得 CDF 的反函数，在0到1之间均匀采样，代入反函数，就得到了采样点。但是实际大部分概率分布不能得到 CDF。
Rejection Sampling 拒绝采样：对于概率分布 $p(z)$，引入简单的提议分布 $q(z)$，使得 $\forall z_i,Mq(z_i)\ge p(z_i)$。我们先在 $ q(z)$ 中采样，定义接受率：$\alpha=\frac{p(z^i)}{Mq(z^i)}\le1$。算法描述为：
1. 取 $z^i\sim q(z)$。
2. 在均匀分布中选取 $u$。
3. 如果 $u\le\alpha$，则接受 $z^i$，否则，拒绝这个值。
Importance Sampling：直接对期望：$\mathbb{E}_{p(z)}[f(z)]$ 进行采样。 \[ \mathbb{E}_{p(z)}[f(z)]=\int p(z)f(z)dz=\int \frac{p(z)}{q(z)}f(z)q(z)dz\simeq\frac{1}{N}\sum\limits_{i=1}^Nf(z_i)\frac{p(z_i)}{q(z_i)} \] 于是采样在 $ q(z)$ 中采样，并通过权重计算和。重要值采样对于权重非常小的时候，效率非常低。重要性采样有一个变种 Sampling-Importance-Resampling，这种方法，首先和上面一样进行采样，然后在采样出来的 $N$ 个样本中，重新采样，这个重新采样，使用每个样本点的权重作为概率分布进行采样。

12.2 MCMC

马尔可夫链式一种时间状态都是离散的随机变量序列。我们关注的主要是齐次的一阶马尔可夫链。马尔可夫链满足：$p(X_{t+1}|X_1,X_2,\cdots,X_t)=p(X_{t+1}|X_t)$。这个式子可以写成转移矩阵的形式 $p_{ij}=p(X_{t+1}=j|X_t=i)$。我们有： \[ \pi_{t+1}(x^*)=\int\pi_i(x)p_{x\to x^*}dx \] 如果存在 $\pi=(\pi(1),\pi(2),\cdots),\sum\limits_{i=1}^{+\infin}\pi(i)=1$，有上式成立，这个序列就叫马尔可夫链 $X_t$ 的平稳分布，平稳分布就是表示在某一个时刻后，分布不再改变。MCMC 就是通过构建马尔可夫链概率序列，使其收敛到平稳分布 $p(z)$。引入细致平衡：$\pi(x)p_{x\to x^*}=\pi(x^*)p_{x^*\to x}$。如果一个分布满足细致平衡，那么一定满足平稳分布（反之不成立）： \[ \int\pi(x)p_{x\to x^*}dx=\int\pi(x^*)p_{x^*\to x}dx=\pi(x^*) \] 细致平衡条件将平稳分布的序列和马尔可夫链的转移矩阵联系在一起了，通过转移矩阵可以不断生成样本点。假定随机取一个转移矩阵 $(Q=Q_{ij})$，作为一个提议矩阵。我们有： \[ p(z)\cdot Q_{z\to z^*}\alpha(z,z^*)=p(z^*)\cdot Q_{z^*\to z}\alpha(z^*,z) \] 取： \[ \alpha(z,z^*)=\min\{1,\frac{p(z^*)Q_{z^*\to z}}{p(z)Q_{z\to z^*}}\} \] 则 \[ p(z)\cdot Q_{z\to z^*}\alpha(z,z^*)=\min\{p(z)Q_{z\to z^*},p(z^*)Q_{z^*\to z}\}=p(z^*)\cdot Q_{z^*\to z}\alpha(z^*,z) \] 于是，迭代就得到了序列，这个算法叫做 Metropolis-Hastings 算法：

通过在0，1之间均匀分布取点 $u$
生成 $z^*\sim Q(z^*|z^{i-1})$
计算 $\alpha$ 值
如果 $\alpha\ge u$，则 $z^i=z^*$，否则 $z^{i}=z^{i-1}$

这样取的样本就服从 $p(z)=\frac{\hat{p}(z)}{z_p}\sim \hat{p}(z)$。

下面介绍另一种采样方式 Gibbs 采样，如果 $z$ 的维度非常高，那么通过固定被采样的维度其余的维度来简化采样过程：$z_i\sim p(z_i|z_{-i})$：

给定初始值 $z_1^0,z_2^0,\cdots$
在 $t+1$ 时刻，采样 $z_i^{t+1}\sim p(z_i|z_{-i})$，从第一个维度一个个采样。

Gibbs 采样方法是一种特殊的 MH 采样，可以计算 Gibbs 采样的接受率： \[ \frac{p(z^*)Q_{z^*\to z}}{p(z)Q_{z\to z^*}}=\frac{p(z_i^*|z^*_{-i})p(z^*_{-i})p(z_i|z_{-i}^*)}{p(z_i|z_{-i})p(z_{-i})p(z_i^*|z_{-i})} \] 对于每个 Gibbs 采样步骤，$z_{-i}=z_{-i}^*$，这是由于每个维度 $i$ 采样的时候，其余的参量保持不变。所以上式为1。于是 Gibbs 采样过程中，相当于找到了一个步骤，使得所有的接受率为 1。

12.3 平稳分布

定义随机矩阵： \[ Q=\begin{pmatrix}Q_{11}&Q_{12}&\cdots&Q_{1K}\\\vdots&\vdots&\vdots&\vdots\\Q_{k1}&Q_{k2}&\cdots&Q_{KK}\end{pmatrix} \] 这个矩阵每一行或者每一列的和都是1。随机矩阵的特征值都小于等于1。假设只有一个特征值为 $\lambda_i=1$。于是在马尔可夫过程中： \[ q^{t+1}(x=j)=\sum\limits_{i=1}^Kq^t(x=i)Q_{ij}\\ \Rightarrow q^{t+1}=q^t\cdot Q=q^1Q^t \] 于是有： \[ q^{t+1}=q^1A\Lambda^t A^{-1} \] 如果 $m$ 足够大，那么，$\Lambda^m=diag(0,0,\cdots,1,\cdots,0)$，则：$q^{m+1}=q^{m}$ ，则趋于平稳分布了。马尔可夫链可能具有平稳分布的性质，所以我们可以构建马尔可夫链使其平稳分布收敛于需要的概率分布（设计转移矩阵）。

在采样过程中，需要经历一定的时间（燃烧期/混合时间）才能达到平稳分布。但是 MCMC 方法有一些问题：

无法判断是否已经收敛
燃烧期过长（维度太高，并且维度之间有关，可能无法采样到某些维度），例如在 GMM 中，可能无法采样到某些峰。于是在一些模型中，需要对隐变量之间的关系作出约束，如 RBM 假设隐变量之间无关。
样本之间一定是有相关性的，如果每个时刻都取一个点，那么每个样本一定和前一个相关，这可以通过间隔一段时间采样。

12.4 隐马尔可夫模型

隐马尔可夫模型是一种概率图模型。我们知道，机器学习模型可以从频率派和贝叶斯派两个方向考虑，在频率派的方法中的核心是优化问题，而在贝叶斯派的方法中，核心是积分问题，也发展出来了一系列的积分方法如变分推断，MCMC 等。概率图模型最基本的模型可以分为有向图（贝叶斯网络）和无向图（马尔可夫随机场）两个方面，例如 GMM，在这些基本的模型上，如果样本之间存在关联，可以认为样本中附带了时序信息，从而样本之间不独立全同分布的，这种模型就叫做动态模型，隐变量随着时间发生变化，于是观测变量也发生变化：

graph LR;
    z1-->z2-->z3;

根据状态变量的特点，可以分为：

HMM，状态变量（隐变量）是离散的
Kalman 滤波，状态变量是连续的，线性的
粒子滤波，状态变量是连续，非线性的

12.5 HMM

HMM 用概率图表示为：

graph TD;
t1-->t2;
subgraph four
    t4-->x4((x4))
end
subgraph three
    t3-->x3((x3))
end
subgraph two
    t2-->x2((x2))
end
subgraph one
    t1-->x1((x1))
end

t2-->t3;
t3-->t4;

上图表示了四个时刻的隐变量变化。用参数 $\lambda=(\pi,A,B)$ 来表示，其中 $\pi$ 是开始的概率分布，$A$ 为状态转移矩阵，$B$ 为发射矩阵。

下面使用 $ o_t$ 来表示观测变量，$O$ 为观测序列，$V=\{v_1,v_2,\cdots,v_M\}$ 表示观测的值域，$i_t$ 表示状态变量，$I$ 为状态序列，$Q=\{q_1,q_2,\cdots,q_N\}$ 表示状态变量的值域。定义 $A=(a_{ij}=p(i_{t+1}=q_j|i_t=q_i))$ 表示状态转移矩阵，$B=(b_j(k)=p(o_t=v_k|i_t=q_j))$ 表示发射矩阵。

在 HMM 中，有两个基本假设：

齐次 Markov 假设（未来只依赖于当前）： \[ p(i_{t+1}|i_t,i_{t-1},\cdots,i_1,o_t,o_{t-1},\cdots,o_1)=p(i_{t+1}|i_t) \]
观测独立假设： \[ p(o_t|i_t,i_{t-1},\cdots,i_1,o_{t-1},\cdots,o_1)=p(o_t|i_t) \]

HMM 要解决三个问题：

Evaluation：$p(O|\lambda)$，Forward-Backward 算法
Learning：$\lambda=\mathop{argmax}\limits_{\lambda}p(O|\lambda)$，EM 算法（Baum-Welch）
Decoding：$I=\mathop{argmax}\limits_{I}p(I|O,\lambda)$，Vierbi 算法
1. 预测问题：$p(i_{t+1}|o_1,o_2,\cdots,o_t)$
2. 滤波问题：$p(i_t|o_1,o_2,\cdots,o_t)$

12.5.1 Evaluation

\[ p(O|\lambda)=\sum\limits_{I}p(I,O|\lambda)=\sum\limits_{I}p(O|I,\lambda)p(I|\lambda) \]

\[ p(I|\lambda)=p(i_1,i_2,\cdots,i_t|\lambda)=p(i_t|i_1,i_2,\cdots,i_{t-1},\lambda)p(i_1,i_2,\cdots,i_{t-1}|\lambda) \]

根据齐次 Markov 假设： \[ p(i_t|i_1,i_2,\cdots,i_{t-1},\lambda)=p(i_t|i_{t-1})=a_{i_{t-1}i_t} \] 所以： \[ p(I|\lambda)=\pi_1\prod\limits_{t=2}^Ta_{i_{t-1},i_t} \] 又由于： \[ p(O|I,\lambda)=\prod\limits_{t=1}^Tb_{i_t}(o_t) \] 于是： \[ p(O|\lambda)=\sum\limits_{I}\pi_{i_1}\prod\limits_{t=2}^Ta_{i_{t-1},i_t}\prod\limits_{t=1}^Tb_{i_t}(o_t) \] 我们看到，上面的式子中的求和符号是对所有的观测变量求和，于是复杂度为 $O(N^T)$。

下面，记 $\alpha_t(i)=p(o_1,o_2,\cdots,o_t,i_t=q_i|\lambda)$，所以，$\alpha_T(i)=p(O,i_T=q_i|\lambda)$。我们看到： \[ p(O|\lambda)=\sum\limits_{i=1}^Np(O,i_T=q_i|\lambda)=\sum\limits_{i=1}^N\alpha_T(i) \] 对 $\alpha_{t+1}(j)$： \[ \begin{align}\alpha_{t+1}(j)&=p(o_1,o_2,\cdots,o_{t+1},i_{t+1}=q_j|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_{t+1},i_{t+1}=q_j,i_t=q_i|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Np(o_{t+1}|o_1,o_2,\cdots,i_{t+1}=q_j,i_t=q_i|\lambda)p(o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j|\lambda) \end{align} \] 利用观测独立假设： \[ \begin{align}\alpha_{t+1}(j)&=\sum\limits_{i=1}^Np(o_{t+1}|i_{t+1}=q_j)p(o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Np(o_{t+1}|i_{t+1}=q_j)p(i_{t+1}=q_j|o_1,\cdots,o_t,i_t=q_i,\lambda)p(o_1,\cdots,o_t,i_t=q_i|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Nb_{j}(o_t)a_{ij}\alpha_t(i) \end{align} \] 上面利用了齐次 Markov 假设得到了一个递推公式，这个算法叫做前向算法。

还有一种算法叫做后向算法，定义 $\beta_t(i)=p(o_{t+1},o_{t+1},\cdots，o_T|i_t=i,\lambda)$： \[ \begin{align}p(O|\lambda)&=p(o_1,\cdots,o_T|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_T,i_1=q_i|\lambda)\nonumber\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_T|i_1=q_i,\lambda)\pi_i\nonumber\\ &=\sum\limits_{i=1}^Np(o_1|o_2,\cdots,o_T,i_1=q_i,\lambda)p(o_2,\cdots,o_T|i_1=q_i,\lambda)\pi_i\nonumber\\ &=\sum\limits_{i=1}^Nb_i(o_1)\pi_i\beta_1(i) \end{align} \] 对于这个 $\beta_1(i)$： \[ \begin{align}\beta_t(i)&=p(o_{t+1},\cdots,o_T|i_t=q_i)\nonumber\\ &=\sum\limits_{j=1}^Np(o_{t+1},o_{t+2},\cdots,o_T,i_{t+1}=q_j|i_t=q_i)\nonumber\\ &=\sum\limits_{j=1}^Np(o_{t+1},\cdots,o_T|i_{t+1}=q_j,i_t=q_i)p(i_{t+1}=q_j|i_t=q_i)\nonumber\\ &=\sum\limits_{j=1}^Np(o_{t+1},\cdots,o_T|i_{t+1}=q_j)a_{ij}\nonumber\\ &=\sum\limits_{j=1}^Np(o_{t+1}|o_{t+2},\cdots,o_T,i_{t+1}=q_j)p(o_{t+2},\cdots,o_T|i_{t+1}=q_j)a_{ij}\nonumber\\ &=\sum\limits_{j=1}^Nb_j(o_{t+1})a_{ij}\beta_{t+1}(j) \end{align} \] 于是后向地得到了第一项。

12.5.2 Learning

为了学习得到参数的最优值，在 MLE 中： \[ \lambda_{MLE}=\mathop{argmax}_\lambda p(O|\lambda) \] 我们采用 EM 算法（在这里也叫 Baum Welch 算法），用上标表示迭代： \[ \theta^{t+1}=\mathop{argmax}_{\theta}\int_z\log p(X,Z|\theta)p(Z|X,\theta^t)dz \] 其中，$X$ 是观测变量，$Z$ 是隐变量序列。于是： \[ \lambda^{t+1}=\mathop{argmax}_\lambda\sum\limits_I\log p(O,I|\lambda)p(I|O,\lambda^t)\\ =\mathop{argmax}_\lambda\sum\limits_I\log p(O,I|\lambda)p(O,I|\lambda^t) \] 这里利用了 $p(O|\lambda^t)$ 和$\lambda$ 无关。将 Evaluation 中的式子代入： \[ \sum\limits_I\log p(O,I|\lambda)p(O,I|\lambda^t)=\sum\limits_I[\log \pi_{i_1}+\sum\limits_{t=2}^T\log a_{i_{t-1},i_t}+\sum\limits_{t=1}^T\log b_{i_t}(o_t)]p(O,I|\lambda^t) \] 对 $\pi^{t+1}$： \[ \begin{align}\pi^{t+1}&=\mathop{argmax}_\pi\sum\limits_I[\log \pi_{i_1}p(O,I|\lambda^t)]\nonumber\\ &=\mathop{argmax}_\pi\sum\limits_I[\log \pi_{i_1}\cdot p(O,i_1,i_2,\cdots,i_T|\lambda^t)] \end{align} \] 上面的式子中，对 $i_2,i_2,\cdots,i_T$ 求和可以将这些参数消掉： \[ \pi^{t+1}=\mathop{argmax}_\pi\sum\limits_{i_1}[\log \pi_{i_1}\cdot p(O,i_1|\lambda^t)] \] 上面的式子还有对 $\pi$ 的约束 $\sum\limits_i\pi_i=1$。定义 Lagrange 函数： \[ L(\pi,\eta)=\sum\limits_{i=1}^N\log \pi_i\cdot p(O,i_1=q_i|\lambda^t)+\eta(\sum\limits_{i=1}^N\pi_i-1) \] 于是： \[ \frac{\partial L}{\partial\pi_i}=\frac{1}{\pi_i}p(O,i_1=q_i|\lambda^t)+\eta=0 \] 对上式求和： \[ \sum\limits_{i=1}^Np(O,i_1=q_i|\lambda^t)+\pi_i\eta=0\Rightarrow\eta=-p(O|\lambda^t) \] 所以： \[ \pi_i^{t+1}=\frac{p(O,i_1=q_i|\lambda^t)}{p(O|\lambda^t)} \]

12.5.3 Decoding

Decoding 问题表述为： \[ I=\mathop{argmax}\limits_{I}p(I|O,\lambda) \] 我们需要找到一个序列，其概率最大，这个序列就是在参数空间中的一个路径，可以采用动态规划的思想。

定义： \[ \delta_{t}(j)=\max\limits_{i_1,\cdots,i_{t-1}}p(o_1,\cdots,o_t,i_1,\cdots,i_{t-1},i_t=q_i) \] 于是： \[ \delta_{t+1}(j)=\max\limits_{1\le i\le N}\delta_t(i)a_{ij}b_j(o_{t+1}) \] 这个式子就是从上一步到下一步的概率再求最大值。记这个路径为： \[ \psi_{t+1}(j)=\mathop{argmax}\limits_{1\le i\le N}\delta_t(i)a_{ij} \]

12.6 小结

HMM 是一种动态模型，是由混合树形模型和时序结合起来的一种模型（类似 GMM + Time）。对于类似 HMM 的这种状态空间模型，普遍的除了学习任务（采用 EM ）外，还有推断任务，推断任务包括：

译码 Decoding：$p(z_1,z_2,\cdots,z_t|x_1,x_2,\cdots,x_t)$
似然概率：$p(X|\theta)$
滤波：$p(z_t|x_1,\cdots,x_t)$，Online \[ p(z_t|x_{1:t})=\frac{p(x_{1:t},z_t)}{p(x_{1:t})}=C\alpha_t(z_t) \]
平滑：$p(z_t|x_1,\cdots,x_T)$，Offline \[ p(z_t|x_{1:T})=\frac{p(x_{1:T},z_t)}{p(x_{1:T})}=\frac{\alpha_t(z_t)p(x_{t+1:T}|x_{1:t},z_t)}{p(x_{1:T})} \] 根据概率图的条件独立性，有： \[ p(z_t|x_{1:T})=\frac{\alpha_t(z_t)p(x_{t+1:T}|z_t)}{p(x_{1:T})}=C\alpha_t(z_t)\beta_t(z_t) \] 这个算法叫做前向后向算法。
预测：$p(z_{t+1},z_{t+2}|x_1,\cdots,x_t),p(x_{t+1},x_{t+2}|x_1,\cdots,x_t)$ \[ p(z_{t+1}|x_{1:t})=\sum_{z_t}p(z_{t+1},z_t|x_{1:t})=\sum\limits_{z_t}p(z_{t+1}|z_t)p(z_t|x_{1:t}) \]

\[ p(x_{t+1}|x_{1:t})=\sum\limits_{z_{t+1}}p(x_{t+1},z_{t+1}|x_{1:t})=\sum\limits_{z_{t+1}}p(x_{t+1}|z_{t+1})p(z_{t+1}|x_{1:t}) \]