3 MCMC 采样与 M-H 采样

MCMC(三)MCMC采样和M-H采样 - 刘建平Pinard - 博客园

3.1 细致平稳条件：

$$\pi(i)P(i,j) = \pi(j)P(j,i).$$

其中 $\pi\in\mathbb{R}^{1\times n}$ 为平稳分布，$P(i,j)$ 表示从状态 i 到 j 的概率。

细致平衡条件表明，在平稳分布下，任意两个状态之间交换的概率"流量"是相等的。这意味着在稳态下，从状态 i 到 j 的概率流恰好等于从 j 到 i 的概率流，形成微观层面的平衡。

证明：左右两端对 i 求和，有

$$\sum_{i=1}^{n}\pi(i)P(i,j) = \sum_{i=1}^{n}\pi(j)P(j,i)$$

假设当前时刻为 t，则左端等于 $\pi_{t+1}(j)$，右端等于 $\pi(j)\sum_i P(j,i) = \pi(j) = \pi_t(j)$。

由于 j 是任意值，则我们有 $\pi_{t+1} = \pi_t$。

3.2 MCMC 采样

一般情况下，我们并不能具有 3.1 所述的细致平稳条件，即假设我们有一个目标平稳分布 $\pi$ 和一个马尔科夫链状态转移矩阵 $Q$：

$$\pi(i)Q(i,j)\neq \pi(j)Q(j,i),$$

这个时候我们就需要对上述条件进行改造，使其成立：

$$\begin{aligned} &\pi(i)\alpha(i,j)Q(i,j)= \pi(j)\alpha(j,i)Q(j,i),\\ \text{where }& \alpha(i,j)=\pi(j)Q(j,i),\\ &\alpha(j,i)=\pi(i)Q(i,j) \end{aligned}$$

MCMC 采样过程

1. 输入我们任意选定的马尔科夫链状态转移矩阵$Q$，平稳分布$\pi(x)$，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$

2. 从任意简单概率分布采样得到初始状态值$x_0$

3. for $t = 0$ to $n_1 + n_2 - 1$:

   a) 从$Q(x_*|x_t)$中采样得到样本$x_*$

   b) 从均匀分布采样$u \sim U[0,1]$

   c) 如果$u < \alpha(x_t,x_*) = \pi(x_*)Q(x_*,x_t)$，则接受转移$x_t \to x_*$，即$x_{t+1} = x_*$

   d) 否则不接受转移，即$x_{t+1} = x_t$

样本集$(x_{n_1},x_{n_1+1},\ldots,x_{n_1+n_2-1})$即为我们需要的平稳分布对应的样本集。

---

对于转移核$P(i,j)=Q(i,j)\cdot\alpha(i,j)$ 这个导出来的表达式，可以用下面这个想法去理解：

1. 目标分布 $\pi$ 是马尔科夫链 $A$ 的稳定分布，$P(i,j)$是 $A$ 的转移核；
2. $Q(i,j)$ 是另一个马尔科夫链 $B$ 的转移核，$A$ 和 $B$ 有相同的状态空间；
3. 对于 $P(i,j)$，它等于事件“i到j能够在B转移”和事件“A接受这种转移”的联合概率，
   • 第1个事件的概率就是等于 $Q(i,j)$，第2个事件的概率等价于均匀分布随机变量 $u$ 的 CDF，即 $P(u\leq \alpha(i,j)) = \alpha(i,j)$；
4. 由于这两个事件是相互独立的，所以联合事件的概率等于这两个概率的乘积，即 $A$ 的转移核 $P(i,j)=Q(i,j)\cdot\alpha(i,j)$.
5. 对应 $\alpha(i,j)$ 经过精心设计，转移核 $P(i,j)$ 也满足在马尔科夫链 $A$ 下的 detailed balance 条件，也就可以满足 $A$ 最后的平稳分布可以是 $\pi$。

3.3 Metropolis-Hastings 采样

上述分析过程并未考虑 $\alpha(i,j)$ 的实际大小，事实上，这个值一般会比较小，导致 MCMC 采样效率比较低，因此 Metropolis-Hastings 对其作了如下修正，使采样效率提高。

$$\alpha(i,j) = \min\{\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}, 1\}$$

注：选择 min 是因为我们判断过程中的 $u \sim U[0,1]$.

在进行 MH 采样时，有以下几点需要注意：

• 初始值 $x_0$ 可以随意给定；
• 基于前一个状态的值 $x_{\text{last}}$ 和提议分布（即条件分布） $Q(x_{\text{next}} | x_{\text{last}})$，生成候选参数；
• 基于前一个状态的值、当前候选参数的值，以及目标平稳分布 $\pi$，我们可以计算接受概率 $\alpha$；
• 从均匀分布采样一个随机数 $u\sim U(0,1)$：
  • 如果 $u<\alpha$，设定当前状态为候选参数；
  • 否则，拒绝候选，并将当前状态设为前一个状态。

实例（参考./MCMC.ipynb）：

• 目标平稳分布 $\pi$ 是一个均值3，标准差2的正态分布——$\pi(x)$即为 $\pi$ 在位置 $x$ 处的概率密度值；
• $Q(i,j)$ 是以 $i$ 为均值，1为方差的正态分布在位置 $j$ 的概率密度值
  • $Q$：提议分布 / 条件分布
  • 我们将会用到 $Q(x_t, x^*)$，即从状态 $x_t$ 转移到 $x^*$ 的概率。
• $x_0 = 0$
  • $x^*$ 是一个服从 $N(x_t,1)$ 的随机数
  • 由于 $Q(i,j) = Q(j,i)$，则接受概率的计算简化为 $\alpha(x_t,x^*) = \min\{1, \frac{\pi(x^*)Q(x^*, x_t)}{\pi(x_t)Q(x_t,x^*)} \} = \min\{1, \frac{\pi(x^*)}{\pi(x_t)} \}$