Metropolis-Hastings 算法

19.4 Metropolis-Hastings 算法

本节叙述 Metropolis-Hastings 算法。该算法是马尔可夫链蒙特卡罗法的代表算法。

---

19.4.1 基本原理

1. 马尔可夫链

假设要抽样的概率分布为 $p(x)$。Metropolis-Hastings 算法采用转移核为 $p(x,x^{\prime })$ 的马尔可夫链：
$$p(x,x^{\prime })=q(x,x^{\prime })\alpha (x,x^{\prime })\text{\hspace{1em}(19.38)}$$

其中，$q(x,x^{\prime })$ 和 $\alpha (x,x^{\prime })$ 分别称为建议分布（proposal distribution）和接受分布（acceptance distribution）。

建议分布 $q(x,x^{\prime })$ 是另一个马尔可夫链的转移核，并且 $q(x,x^{\prime })$ 是不可约的，即其概率值恒不为 0，同时是一个容易抽样的分布。接受分布 $\alpha (x,x^{\prime })$ 是：
$$\alpha (x,x^{\prime })=\min \left\{1,\frac{p(x^{\prime })q(x^{\prime },x)}{p(x)q(x,x^{\prime })}\right\}\text{\hspace{1em}(19.39)}$$

这样，转移核 $p(x,x^{\prime })$ 可以写成：
$$p(x,x^{\prime })=\{\begin{array}{ll}q(x,x^{\prime }),& p(x^{\prime })q(x^{\prime },x)\ge p(x)q(x,x^{\prime })\\ q(x^{\prime },x)\frac{p(x^{\prime })}{p(x)},& p(x^{\prime })q(x^{\prime },x)<p(x)q(x,x^{\prime })\end{array}\text{\hspace{1em}(19.40)}$$

转移核为 $p(x,x^{\prime })$ 的马尔可夫链上的随机游走以以下方式进行：如果在时刻 $t-1$ 处于状态 $x$, 即 $x_{t-1}=x$，则按照建议分布 $q(x,x^{\prime })$ 抽样产生一个候选状态 $x^{\prime }$，然后按照接受分布 $\alpha (x,x^{\prime })$ 抽样决定是否接受状态 $x^{\prime }$。以概率 $\alpha (x,x^{\prime })$ 接受 $x^{\prime }$，决定时刻 $t$ 转移到状态 $x^{\prime }$，而以概率 $1-\alpha (x,x^{\prime })$ 拒绝 $x^{\prime }$，决定时刻 $t$ 仍停留在状态 $x$。具体地，从区间 $[0,1]$ 的均匀分布中抽取一个随机数 $u$，决定时刻 $t$ 的状态为：
$$x_t=\{\begin{array}{ll}{x}^{\prime },& u\le \alpha (x,x^{\prime })\\ x,& u>\alpha (x,x^{\prime })\end{array}$$

可以证明，转移核为 $p(x,x^{\prime })$ 的马尔可夫链是可逆马尔可夫链（满足遍历定理），其平稳分布就是 $p(x)$，即要抽样的目标分布。也就是说这是马尔可夫链蒙特卡罗法的一个具体实现。

定理 19.6 由转移核 (19.38) ～ (19.40) 构成的马尔可夫链是可逆的，即：
$$p(x)p(x,x^{\prime })=p(x^{\prime })p(x^{\prime },x)\text{\hspace{1em}(19.41)}$$

并且 $p(x)$ 是该马尔可夫链的平稳分布。

证明 若 $x=x^{\prime }$，则式 (19.41) 显然成立。

设 $x\ne x^{\prime }$，则：
p ( x ) p ( x , x ′ ) = p ( x ) q ( x , x ′ ) min ⁡ { 1 , p ( x ′ ) q ( x ′ , x ) p ( x ) q ( x , x ′ ) } = min ⁡ { p ( x ) q ( x , x ′ ) , p ( x ′ ) q ( x ′ , x ) } = p ( x ′ ) q ( x ′ , x ) min ⁡ { p ( x ) q ( x , x ′ ) p ( x ′ ) q ( x ′ , x ) , 1 } = p ( x ′ ) q ( x ′ , x ) \begin{aligned} p(x)p(x, x') &= p(x)q(x, x') \min \left\{ 1, \frac{p(x')q(x', x)}{p(x)q(x, x')} \right\} \notag \\ &= \min \{p(x)q(x, x'), p(x')q(x', x)\} \notag \\ &= p(x')q(x', x) \min \left\{ \frac{p(x)q(x, x')}{p(x')q(x', x)}, 1 \right\} \notag \\ &= p(x')q(x', x) \end{aligned} p(x)p(x,x′)​=p(x)q(x,x′)min{1,p(x)q(x,x′)p(x′)q(x′,x)​}=min{p(x)q(x,x′),p(x′)q(x′,x)}=p(x′)q(x′,x)min{p(x′)q(x′,x)p(x)q(x,x′)​,1}=p(x′)q(x′,x)​
式 (19.41) 成立。

由式 (19.41) 知：
$$\begin{array}{rl}\int p(x)p(x,x^{\prime })dx& =\int p(x^{\prime })p(x^{\prime },x)dx\\ & =p(x^{\prime })\int p(x^{\prime },x)dx\\ & =p(x^{\prime })\end{array}$$

根据平稳分布的定义（式 (19.21)），$p(x)$ 是马尔可夫链的平稳分布。

2. 建议分布

建议分布 $q(x,x^{\prime })$ 有多种可能的形式，这里介绍两种常用形式。

第一种形式：假设建议分布是对称的，即对任意的 $x$ 和 $x^{\prime }$ 有：
$$q(x,x^{\prime })=q(x^{\prime },x)\text{\hspace{1em}(19.42)}$$

这样的建议分布称为 Metropolis 选择，也是 Metropolis-Hastings 算法最初采用的建议分布。
此时，接受分布 $\alpha (x,x^{\prime })$ 简化为：
$$\alpha (x,x^{\prime })=\min \left\{1,\frac{p(x^{\prime })}{p(x)}\right\}\text{\hspace{1em}(19.43)}$$

Metropolis 选择的一个特例是 $q(x,x^{\prime })$ 取条件概率分布 $p(x^{\prime }|x)$，定义为多元正态分布，其均值是 $x$，协方差矩阵是常数矩阵。

Metropolis 选择的另一个特例是令 $q(x,x^{\prime })=q(|x-x^{\prime }|)$，这时算法称为随机游走 Metropolis 算法。例如：
$$q(x,x^{\prime })\propto \exp \left[-\frac{(x^{\prime }-x{)}^2}{2}\right]$$

Metropolis 选择的特点是当 $x^{\prime }$ 与 $x$ 接近时，$q(x,x^{\prime })$ 的概率值高；否则，$q(x,x^{\prime })$ 的概率值低。状态转移在附近点的可能性更大。

第二种形式称为独立抽样。假设 $q(x,x^{\prime })=q(x^{\prime })$ 与当前状态 $x$ 无关，即 $q(x,x^{\prime })=q(x^{\prime })$。建议分布的计算按照 $q(x^{\prime })$ 独立抽样进行。此时，接受分布 $\alpha (x,x^{\prime })$ 可以写成：
$$\alpha (x,x^{\prime })=\min \left\{1,\frac{w(x^{\prime })}{w(x)}\right\},\text{\hspace{1em}(19.44)}$$

其中：$w(x^{\prime })=p(x^{\prime })/q(x^{\prime }),w(x)=p(x)/q(x).$

独立抽样实现简单，但可能收敛速度慢，通常选择接近目标分布 $p(x)$ 的分布作为建议分布 $q(x)$。

3. 满条件分布

马尔可夫链蒙特卡罗法的目标分布通常是多元联合概率分布 $p(x)=p(x_1,x_2,\cdots ,x_k)$，其中 $x=(x_1,x_2,\cdots ,x_k{)}^T$ 为 $k$ 维随机变量。如果条件概率分布 $p(x_I|x_{-I})$ 中所有 $k$ 个变量全部出现，其中 x I = { x i , i ∈ I } , x − I = { x i , i ∉ I } , I ⊆ K = { 1 , 2 , ⋯   , k } x_{I} = \{x_i, i \in I\}, x_{-I} = \{x_i, i \notin I\}, I \subseteq K = \{1, 2, \cdots, k\} xI​={xi​,i∈I},x−I​={xi​,i∈/I},I⊆K={1,2,⋯,k}，那么称这种条件概率分布为满条件分布（full conditional distribution）。

满条件分布有以下性质：对任意的 $x\in X$ 和任意的 $I\subseteq K$，有：
$$p(x_I|x_{-I})=\frac{p(x)}{\int p(x)d{x}_I}\propto p(x)\text{\hspace{1em}(19.45)}$$

而且，对任意的 $x,x^{\prime }\in X$ 和任意的 $I\subseteq K$，有：

$$\frac{p(x_I^{\prime }|x_{-I}^{\prime })}{p(x_I|x_{-I})}=\frac{p(x^{\prime })}{p(x)}\text{\hspace{1em}(19.46)}$$

Metropolis-Hastings 算法中，可以利用性质 (19.46) 简化计算，提高计算效率。具体地，通过满条件分布概率的比 $\frac{p(x_I^{\prime }|x_{-I}^{\prime })}{p(x_I|x_{-I})}$ 计算联合概率的比 $\frac{p(x^{\prime })}{p(x)}$ 而前者更容易计算。

例 19.9 设 $x_1$ 和 $x_2$ 的联合概率分布的密度函数为：

$$p(x_1,x_2)\propto \exp \left[-\frac{1}{2}(x_1-1{)}^2(x_2-1{)}^2\right]$$

求其满条件分布。

解 由满条件分布的定义有：
$$\begin{array}{rl}p(x_1\mid x_2)& \propto p(x_1,x_2)\\ & \propto \exp \left[-\frac{1}{2}(x_1-1{)}^2(x_2-1{)}^2\right]\\ & \propto N\left(1,(x_2-1{)}^{-2}\right)\end{array}$$

这里 $N(1,(x_2-1{)}^{-2})$ 是均值为 1，方差为 $(x_2-1{)}^{-2}$ 的正态分布。这时 $x_1$ 是变量，$x_2$ 是参数。同样可得：
$$\begin{array}{rl}p(x_2\mid x_1)& \propto p(x_1,x_2)\\ & \propto \exp \left[-\frac{1}{2}(x_2-1{)}^2(x_1-1{)}^2\right]\\ & \propto N\left(1,(x_1-1{)}^{-2}\right)\end{array}$$

19.4.2 Metropolis-Hastings 算法

算法 19.2 (Metropolis-Hastings 算法)
输入：抽样的目标分布的密度函数 $p(x)$，函数 $f(x)$。
输出：$p(x)$ 的随机样本 $x_{m+1},x_{m+2},\cdots ,x_n$，函数样本均值 $f_{mn}$。
参数：收敛步数 $m$，迭代步数 $n$。

1. 任意选择一个初始值 $x_0$。

2. 对 $i=1,2,\cdots ,n$ 循环执行：
   (a) 设状态 $x_{i-1}=x$，按照建议分布 $q(x,x^{\prime })$ 随机抽取一个候选状态 $x^{\prime }$。
   (b) 计算接受概率：
   $$\alpha (x,x^{\prime })=\min \left\{1,\frac{p(x^{\prime })q(x^{\prime },x)}{p(x)q(x,x^{\prime })}\right\}$$

   (c) 从区间 $(0,1)$ 中按均匀分布随机抽取一个数 $u$。若 $u\le \alpha (x,x^{\prime })$，则状态 $x_i=x^{\prime }$；否则，状态 $x_i=x$。

3. 得到样本集合 { x m + 1 , x m + 2 , ⋯   , x n } \{x_{m+1}, x_{m+2}, \cdots, x_n\} {xm+1​,xm+2​,⋯,xn​}，计算
   $$f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

19.4.3 单分量 Metropolis-Hastings 算法

在 Metropolis-Hastings 算法中，通常需要对多元变量分布进行抽样。有时对多元变量分布的抽样是困难的，可以对多元变量的每一变量的条件分布依次分别进行抽样，从而实现对整个多元变量的一次抽样，这就是单分量 Metropolis-Hastings（single-component Metropolis-Hastings）算法。

假设马尔可夫链的状态由 $k$ 维随机变量表示：
$$x=(x_1,x_2,\cdots ,x_k{)}^T$$

其中，$x_j$ 表示随机变量 $x$ 的第 $j$ 个分量，$j=1,2,\cdots ,k$，而 $x^{(i)}$ 表示马尔可夫链在时刻 $i$ 的状态：
$$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T,i=1,2,\cdots ,n$$

其中，$x_j^{(i)}$ 是随机变量 $x^{(i)}$ 的第 $j$ 个分量，$j=1,2,\cdots ,k$。

为了生成容量为 $n$ 的样本集合 { x ( 1 ) , x ( 2 ) , ⋯   , x ( n ) } \{x^{(1)}, x^{(2)}, \cdots, x^{(n)}\} {x(1),x(2),⋯,x(n)}，单分量 Metropolis-Hastings 算法由下面的 $k$ 步迭代来实现 Metropolis-Hastings 算法的一次迭代。

设在第 $i-1$ 次迭代结束时分量 $x_j$ 的取值为 $x_j^{(i-1)}$，在第 $i$ 次迭代的第 $j$ 步，对分量 $x_j$ 根据 Metropolis-Hastings 算法更新，得到其新的取值 $x_j^{(i)}$。首先，由建议分布 $q(x_j^{(i-1)},x_j|x_{-j}^{(i)})$ 抽样产生 $x_j$ 的候选值 $x_j^{\prime (i)}$，这里 $x_{-j}^{(i)}$ 表示在第 $i$ 次迭代的第 $j-1$ 步后的 $x^{(i)}$ 除去 $x_j^{(i-1)}$ 的所有值，即：
$$x_j^{(i-1)}=(x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\cdots ,x_k^{(i-1)}{)}^T$$

其中分量 $1,2,\cdots ,j-1$ 已更新。然后，按照接受概率：
$$\alpha (x_j^{(i-1)},x_j^{\prime (i)}|x_{-j}^{(i)})=\min \left\{1,\frac{p(x_j^{\prime (i)}|x_{-j}^{(i)})q(x_j^{\prime (i)},x_j^{(i-1)}|x_{-j}^{(i)})}{p(x_j^{(i-1)}|x_{-j}^{(i)})q(x_j^{(i-1)},x_j^{\prime (i)}|x_{-j}^{(i)})}\right\}\text{\hspace{1em}(19.47)}$$

抽样决定是否接受候选值 $x_j^{\prime (i)}$。如果 $x_j^{\prime (i)}$ 被接受，则令 $x_j^{(i)}=x_j^{\prime (i)}$；否则，令 $x_j^{(i)}=x_j^{(i-1)}$。其余分量在第 $j$ 步不改变。马尔可夫链的转移核概率为：
$$p\left(x_j^{(i-1)},x_j^{\prime (i)}|x_{-j}^{(i)}\right)=\alpha \left(x_j^{(i-1)},x_j^{\prime (i)}|x_{-j}^{(i)}\right)q\left(x_j^{(i-1)},x_j^{\prime (i)}|x_{-j}^{(i)}\right)\text{\hspace{1em}(19.48)}$$

图 19.10 示意了单分量 Metropolis-Hastings 算法的迭代过程。目标是对含有两个变量的随机变量 $x$ 进行抽样。如果变量 $x_1$ 或 $x_2$ 更新，那么在水平方向或垂直方向产生一个移动，连续水平移动和垂直移动产生一个新的样本点。注意由于建议分布可能不被接受，Metropolis-Hastings 算法可能在一些相邻的时刻不产生移动。

---

基本原理 1.马尔可夫链 理解
19.39当前状态x转移到候选状态x′的接受概率 公式解析
定理19.6 证明马尔可夫链可逆 公式解析
定理19.6 证明p(x)是马尔可夫链可逆的平稳分布 公式解析
基本原理 2.建议分布 理解
基本原理 3.满条件分布 理解
单分量Metropolis-Hastings算法 理解