（《机器学习》完整版系列）第14章 概率图模型——14.7 MH采样算法的详细推导（提议分布，“拒绝采样”）

基于“拒绝采样”的 MH 算法是 MCMC 的重要“艺术”品。
常常有人误解算法中的“拒绝”，这里有必要澄清。 “拒绝”是丢弃了样本，但不是说本次没有产生样本从而浪费了计算，而应理解为： 对本次采样进行了修正，即丢弃新采到的样本，而复制一份先前一步的样本作为本次的样本。

MH算法

上篇的 MCMC 算法并没有给出构造状态转移概率的方案，因此，状态转移概率的设计就成了门“艺术”，其中，基于“拒绝采样”的 MH 算法就是 MCMC 的重要“艺术”品。

设有一提议分布（先验的）$Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$，对其采样得${\mathrm{x}}^{\ast }$，然后再对样本进行修正，而修正策略是“简单粗暴”：

（1）按一定的概率接受${\mathrm{x}}^{\ast }$为${\mathrm{x}}^t$，设接受的概率服从分布$A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$；

（2）当拒绝${\mathrm{x}}^{\ast }$为${\mathrm{x}}^t$时，${\mathrm{x}}^t:={\mathrm{x}}^{t-1}$，即重复一次值。

即转移概率可定义为
$$\begin{array}{rlr}& T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})& \\ & =\{\begin{array}{ll}Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})& (\text{若}{\mathrm{x}}^{\ast }\ne {\mathrm{x}}^{t-1},{\mathrm{x}}^t={\mathrm{x}}^{\ast })\\ Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})& (\text{若}{\mathrm{x}}^{\ast }={\mathrm{x}}^{t-1},{\mathrm{x}}^t={\mathrm{x}}^{t-1})\\ Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})(1-A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1}))& (\text{若}{\mathrm{x}}^{\ast }\ne {\mathrm{x}}^{t-1},{\mathrm{x}}^t={\mathrm{x}}^{t-1})\end{array}& \\ & =\{\begin{array}{ll}Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1}{)}_{{\mathrm{x}}^{\ast }\ne {\mathrm{x}}^{t-1}}A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1}{)}_{{\mathrm{x}}^{\ast }\ne {\mathrm{x}}^{t-1}}& (\text{若}{\mathrm{x}}^t\ne {\mathrm{x}}^{t-1})\\ Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1}{)}_{{\mathrm{x}}^{\ast }={\mathrm{x}}^{t-1}}& \\ +{\left[Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})(1-A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1}))\right]}_{{\mathrm{x}}^{\ast }\ne {\mathrm{x}}^{t-1}}& (\text{若}{\mathrm{x}}^t={\mathrm{x}}^{t-1})\end{array}& \end{array}\text{\hspace{1em}(14.52)}$$

我们需要转移概率满足平稳性（细致平稳）条件式(14.49)，依式(14.52)分两种情况讨论。

I. 当${\mathrm{x}}^t\ne {\mathrm{x}}^{t-1}$时

$T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$为式(14.52)第一式，代入式(14.49)，有
$$\begin{array}{rl}p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)A({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})A({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\end{array}\text{\hspace{1em}(14.53)}$$
即
$$\begin{array}{rl}\frac{A({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}{A({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}& =\frac{p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}\end{array}\text{\hspace{1em}(14.54)}$$
我们找满足式(14.54)的充分条件即可：

i.当$\frac{p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}⩽1$时，取
$$\begin{array}{r}\{\begin{array}{ll}A({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =1\\ A({\mathrm{x}}^t|{\mathrm{x}}^{t-1})& =\frac{p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}\end{array}\end{array}\text{\hspace{1em}(14.55)}$$

ii.当$\frac{p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}>1$时，取
$$\begin{array}{r}\{\begin{array}{ll}A({\mathrm{x}}^t|{\mathrm{x}}^{t-1})& =1\\ A({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =\frac{p({\mathrm{x}}^t)Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^t|{\mathrm{x}}^{t-1})}\end{array}\end{array}\text{\hspace{1em}(14.56)}$$
现在，我们只关心$A({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$，故由式(14.55)、式(14.56)可定义
$$\begin{array}{rl}A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})& =\min \left(1,\frac{p({\mathrm{x}}^{\ast })Q({\mathrm{x}}^{t-1}|{\mathrm{x}}^{\ast })}{p({\mathrm{x}}^{t-1})Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})}\right)\end{array}\text{\hspace{1em}(14.57)}$$
其中，将${\mathrm{x}}^t$改为${\mathrm{x}}^{\ast }$是为了强调此时还未去确定${\mathrm{x}}^t$，此即【西瓜书式(14.28)】。

II. 当${\mathrm{x}}^t={\mathrm{x}}^{t-1}$时

$T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$为式(14.52)的第二式，也可以按上述方法处理，但用不着这么麻烦，因为“常数”序列肯定是平稳的：

设${\mathrm{x}}^t={\mathrm{x}}^{t-1}=\mathrm{x}$，有
$$\begin{array}{r}\{\begin{array}{l}T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})=T(\mathrm{x}|\mathrm{x})=T({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)\\ p({\mathrm{x}}^t)=p(\mathrm{x})=p({\mathrm{x}}^{t-1})\end{array}\end{array}\text{\hspace{1em}(14.58)}$$
式(14.58)代入式(14.49)验证知平稳性（细致平稳）成立，故，此时$T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$无论如何（当然包括式(14.52)第二式）都不影响平稳性（细致平稳）的成立。

我们再证明上述两个平稳马尔可夫链（I.与II.）的“混合体”也平稳：

记式(14.52)的第一式和第二式分别为$T_1({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$和$T_2({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$，则由指示函数及应用（将分段函数表达成一个式子的技术）中的式(B4)
、式(14.52)有
$$\begin{array}{r}T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})=T_1({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\mathbb{I}({\mathrm{x}}^t\ne {\mathrm{x}}^{t-1})+T_2({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\mathbb{I}({\mathrm{x}}^t={\mathrm{x}}^{t-1})\end{array}\text{\hspace{1em}(14.59)}$$
式(14.59)中交换${\mathrm{x}}^t$与${\mathrm{x}}^{t-1}$次序有
$$\begin{array}{r}T({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)=T_1({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)\mathbb{I}({\mathrm{x}}^t\ne {\mathrm{x}}^{t-1})+T_2({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)\mathbb{I}({\mathrm{x}}^t={\mathrm{x}}^{t-1})\end{array}\text{\hspace{1em}(14.60)}$$
由I.与II.的平稳性有
$$\begin{array}{r}\{\begin{array}{ll}p({\mathrm{x}}^t)T_1({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =p({\mathrm{x}}^{t-1})T_1({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\\ p({\mathrm{x}}^t)T_2({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =p({\mathrm{x}}^{t-1})T_2({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\end{array}\end{array}\text{\hspace{1em}(14.61)}$$
故
$$\begin{array}{r}\{\begin{array}{ll}p({\mathrm{x}}^t)T_1({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)\mathbb{I}({\mathrm{x}}^t\ne {\mathrm{x}}^{t-1})& =p({\mathrm{x}}^{t-1})T_1({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\mathbb{I}({\mathrm{x}}^t\ne {\mathrm{x}}^{t-1})\\ p({\mathrm{x}}^t)T_2({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)\mathbb{I}({\mathrm{x}}^t={\mathrm{x}}^{t-1})& =p({\mathrm{x}}^{t-1})T_2({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\mathbb{I}({\mathrm{x}}^t={\mathrm{x}}^{t-1})\end{array}\end{array}\text{\hspace{1em}(14.62)}$$
式(14.62)中两式相加，并应用式(14.60)，有
$$\begin{array}{rl}p({\mathrm{x}}^t)T({\mathrm{x}}^{t-1}|{\mathrm{x}}^t)& =p({\mathrm{x}}^{t-1})T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})\end{array}\text{\hspace{1em}(14.63)}$$
式(14.63)说明：“混合”转移概率$T$满足平稳性（细致平稳）条件式(14.49)。

综上，取式(14.57)的$A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$为接受采样函数（对应有“拒绝采样”），将状态转移概率$T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$定义为式(14.52)，则“混合体”马尔可夫链具有平稳性（细致平稳），从而当足够大$t$时，${\mathrm{x}}^t$近似服从分布$p(\mathrm{x})$，因此，可以在上面实现分布$p(\mathrm{x})$的采样。

将上述状态转移（概率）$T({\mathrm{x}}^t|{\mathrm{x}}^{t-1})$以算法形式实现即为MH算法【西瓜书图14.9】：

$Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$是拍拍脑袋给出的先验分布，作为已知，要采样的分布$p(\mathrm{x})$也是已知。

从先验分布$Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$中采样出${\mathrm{x}}^{\ast }$（第3句）；

利用已知的分布（$Q({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$和$p(\mathrm{x})$）由式(14.57)【西瓜书式(14.28)】计算接收概率$A({\mathrm{x}}^{\ast }|{\mathrm{x}}^{t-1})$，依此值将区间[0,1]分为接受区间$[0,A]$和拒绝区间$(A,1]$，掷骰子$u$碰运气，确定本次是接受还是拒绝（第4-9句）。

需要说明的是：常常有人误解算法中的“拒绝”，这里有必要澄清。 “拒绝”是丢弃了样本，但不是说本次没有产生样本从而浪费了计算，而应理解为： 对本次采样进行了修正，即丢弃新采到的样本，而复制一份先前一步的样本作为本次的样本（算法第8句）。 为什么会出现这种情况呢？这样理解：假定先前一步采样${\mathrm{x}}^{t-1}$是在$p(\mathrm{x})$的高密度区（如，正态分布的中心点附近），即应有较大概率再采到它，而当前步采样${\mathrm{x}}^t$实际并没有从$p(\mathrm{x})$中采，而是由${\mathrm{x}}^{t-1}$转移出${\mathrm{x}}^t$，这之间会有差距，因此需要修正，即调整使之适时地重复出现${\mathrm{x}}^{t-1}$，以保证采样集合中${\mathrm{x}}^{t-1}$的重复度符合分布的要求。

让满足细致平稳条件马尔可夫链充分运行后，得到的值系列为（式(14.43)截去头部）
$$\begin{array}{r}{\mathrm{x}}^T,{\mathrm{x}}^{T+1},{\mathrm{x}}^{T+2},\cdots \end{array}\text{\hspace{1em}(14.64)}$$
则序列(14.64)中的值都可以视为从分布$p(\mathrm{x})$中的采样，从中任取$N$个，则得到采样集，将采样集用于其它的计算，如，【西瓜书式(14.25)】计算函数期望的无偏估计。

取采样集需要注意的是不要在序列中连续取样，因为从马尔可夫链的生成过程知，后项是根据前项生成的，即并不是独立。 为保证采样的相对独立性，实践中通常的做法是：给定一个$n$值（$n$越大，采样样本的越相对独立性），在序列中每隔$n$个取一个，直至取足需要的$N$个样本。

一般的马尔可夫链的转移概率为$T(j|i)$并不满足细致平稳条件，即
$$\begin{array}{r}P(i)T(j|i)\ne P(j)T(i|j)\end{array}\text{\hspace{1em}(14.65)}$$
为表述方便，这里$i={\mathrm{x}}^{t-1},j={\mathrm{x}}^t$。

如何进行改造？

将不满足细致平稳条件的转移概率记为$Q(j|i)$，即
$$\begin{array}{r}P(i)Q(j|i)\ne P(j)Q(i|j)\end{array}\text{\hspace{1em}(14.66)}$$
式(14.66)中的四项相乘，并引入$\alpha (j|i)=P(j)Q(i|j)$，则
$$\begin{array}{r}P(i)P(j)Q(j|i)Q(i|j)=\text{式(14.65)左}\times \alpha (j|i)=\text{式(14.65)右}\times \alpha (i|j)\end{array}$$
即
$$\begin{array}{r}P(i)Q(j|i)\alpha (j|i)=P(j)Q(i|j)\alpha (i|j)\end{array}\text{\hspace{1em}(14.67)}$$
记$T(i|j)=Q(i|j)\alpha (i|j)$，则有
$$\begin{array}{r}P(i)T(j|i)=P(j)T(i|j)\end{array}\text{\hspace{1em}(14.68)}$$

由式(14.68)知，改造后的马尔可夫链满足细致平稳条件，即满足【西瓜书式(14.26)】，而MH算法的【西瓜书式(14.27)】变成了式(14.67)，$\alpha (j|i)$充当了接受概率。

本文为原创，您可以：

• 点赞（支持博主）
• 收藏（待以后看）
• 转发（他考研或学习，正需要）
• 评论（或讨论）
• 引用（支持原创）
• 不侵权

上一篇：14.6 采样：马尔可夫链蒙特卡罗MCMC方法（为何要采样？如何变通采样？）
下一篇：14.8 吉布斯采样算法的详细推导（将“多变量”联合采样变为交替地“单变量”采样）