MetropolisHastings 算法的稳态证明:深入浅出
MetropolisHastings (MH) 算法是蒙特卡洛马尔可夫链 (MCMC) 方法中的一颗璀璨明珠,它能够从一个复杂的目标概率分布中抽取样本,而无需直接计算该分布的归一化常数。其核心思想是构建一个马尔可夫链,使得链的稳态分布恰好是我们的目标分布。那么,我们如何确信MH算法能够达到我们期望的稳态呢?这需要我们深入探究其背后的数学原理。
要证明MH算法能达到稳态,我们需要证明两点:
1. 不可约性 (Irreducibility): 从一个状态出发,能够以非零概率到达所有其他状态。
2. 非周期性 (Aperiodicity): 马尔可夫链不会陷入某种周期性的状态转移模式。
一旦这两个条件满足,理论上我们就能保证一个马尔可夫链存在唯一的稳态分布,并且链会收敛到它。而MH算法的设计正是为了满足这些条件,并确保其稳态分布就是我们的目标分布。
MH算法的核心机制与稳态的联系
在深入证明之前,我们先回顾一下MH算法的工作流程:
1. 初始化: 从一个任意的初始状态 $x^{(0)}$ 开始。
2. 提议 (Proposal): 基于当前状态 $x^{(t)}$,从一个提议分布 $q(x' | x^{(t)})$ 中抽取一个候选状态 $x'$。这个提议分布可以是我们自己选择的。
3. 接受率计算: 计算接受候选状态 $x'$ 的概率,即接受率 $alpha(x' | x^{(t)})$。其公式为:
$$ alpha(x' | x^{(t)}) = minleft(1, frac{pi(x') q(x^{(t)} | x')}{pi(x^{(t)}) q(x' | x^{(t)})}
ight) $$
其中 $pi(cdot)$ 是我们的目标分布。
4. 接受或拒绝: 以概率 $alpha(x' | x^{(t)})$ 接受候选状态 $x'$,即 $x^{(t+1)} = x'$;否则,拒绝候选状态,并保持当前状态,即 $x^{(t+1)} = x^{(t)}$。
正是这个精心设计的接受率,保证了MH算法的马尔可夫链能够收敛到目标分布 $pi(cdot)$。
证明的关键:细致平衡条件 (Detailed Balance Condition)
证明MH算法稳态分布是 $pi(cdot)$ 的核心在于证明其满足细致平衡条件。细致平衡条件是稳态分布的一个充分但不一定必要的条件。它指的是对于马尔可夫链中的任意两个状态 $i$ 和 $j$,在达到稳态后,从状态 $i$ 转移到状态 $j$ 的概率乘以从状态 $j$ 转移到状态 $i$ 的概率,等于从状态 $j$ 转移到状态 $i$ 的概率乘以从状态 $i$ 转移到状态 $j$ 的概率。更具体地说,如果我们用 $P_{ij}$ 表示从状态 $i$ 转移到状态 $j$ 的一步转移概率,并且 $pi_i$ 是状态 $i$ 的稳态概率,那么细致平衡条件可以表示为:
$$ pi_i P_{ij} = pi_j P_{ji} quad ext{for all } i, j $$
如果一个马尔可夫链满足细致平衡条件,并且这个链是不可约且非周期的,那么它一定收敛到以 $pi_i$ 为概率的稳态分布。
验证MH算法的细致平衡条件
现在,让我们来验证MH算法的转移概率 $P_{ij}$ 是否满足细致平衡条件,其中 $i$ 代表状态 $x^{(t)}$,$j$ 代表状态 $x'$。
MH算法的转移概率 $P_{x, x'}$ 可以分解为两个部分:
1. 提议并接受: 从状态 $x$ 提议到状态 $x'$,然后被接受。其概率为 $q(x'|x) alpha(x'|x)$。
2. 提议到其他地方并回到x (或者提议到x自己): 从状态 $x$ 提议到某个状态 $y
eq x'$,然后被拒绝,或者提议到 $x$ 自己,并被接受。由于我们关心的是从 $x$ 到 $x'$ 的转移,而回到 $x$ 的情况属于 $P_{xx}$ 的计算,所以我们主要关注提议到 $x'$ 并接受的情况。
更严谨地说,MH算法的转移概率 $P_{x, x'}$ 可以表示为:
$$ P_{x, x'} = q(x'|x) alpha(x'|x) quad ext{if } x'
eq x $$
$$ P_{x, x} = 1 sum_{x'
eq x} q(x'|x) alpha(x'|x) quad ext{if } x' = x $$
有了这些定义,我们开始验证细致平衡条件:
$$ pi(x) P_{x, x'} = pi(x') P_{x', x} $$
代入MH算法的转移概率公式:
$$ pi(x) left[ q(x'|x) alpha(x'|x)
ight] = pi(x') left[ q(x|x') alpha(x|x')
ight] $$
现在,我们代入接受率 $alpha(x'|x)$ 的定义:
$$ alpha(x'|x) = minleft(1, frac{pi(x') q(x|x')}{pi(x) q(x'|x)}
ight) $$
以及 $alpha(x|x')$ 的定义:
$$ alpha(x|x') = minleft(1, frac{pi(x) q(x'|x)}{pi(x') q(x|x')}
ight) $$
我们发现一个有趣的对称性。让 $A = frac{pi(x') q(x|x')}{pi(x) q(x'|x)}$。那么,
$$ alpha(x'|x) = min(1, A) $$
$$ alpha(x|x') = min(1, 1/A) $$
我们分两种情况讨论:
情况 1:$A ge 1$
在这种情况下,$A ge 1 implies frac{pi(x') q(x|x')}{pi(x) q(x'|x)} ge 1$。
那么:
$alpha(x'|x) = min(1, A) = 1$
$alpha(x|x') = min(1, 1/A) = 1/A$
代入细致平衡条件左侧:
$$ pi(x) left[ q(x'|x) alpha(x'|x)
ight] = pi(x) left[ q(x'|x) cdot 1
ight] = pi(x) q(x'|x) $$
代入细致平衡条件右侧:
$$ pi(x') left[ q(x|x') alpha(x|x')
ight] = pi(x') left[ q(x|x') cdot frac{1}{A}
ight] = pi(x') left[ q(x|x') cdot frac{pi(x) q(x'|x)}{pi(x') q(x|x')}
ight] $$
$$ = pi(x') cdot frac{pi(x) q(x'|x)}{pi(x')} = pi(x) q(x'|x) $$
可以看到,在这种情况下,细致平衡条件 $pi(x) P_{x, x'} = pi(x') P_{x', x}$ 成立。
情况 2:$A < 1$
在这种情况下,$A < 1 implies frac{pi(x') q(x|x')}{pi(x) q(x'|x)} < 1$。
那么:
$alpha(x'|x) = min(1, A) = A$
$alpha(x|x') = min(1, 1/A) = 1$
代入细致平衡条件左侧:
$$ pi(x) left[ q(x'|x) alpha(x'|x)
ight] = pi(x) left[ q(x'|x) cdot A
ight] = pi(x) left[ q(x'|x) cdot frac{pi(x') q(x|x')}{pi(x) q(x'|x)}
ight] $$
$$ = pi(x) cdot frac{pi(x') q(x|x')}{pi(x)} = pi(x') q(x|x') $$
代入细致平衡条件右侧:
$$ pi(x') left[ q(x|x') alpha(x|x')
ight] = pi(x') left[ q(x|x') cdot 1
ight] = pi(x') q(x|x') $$
同样可以看到,在这种情况下,细致平衡条件 $pi(x) P_{x, x'} = pi(x') P_{x', x}$ 也成立。
综合以上两种情况,我们证明了MH算法的转移概率满足细致平衡条件:
$$ pi(x) q(x'|x) alpha(x'|x) = pi(x') q(x|x') alpha(x|x') $$
这仅仅是证明了 MH 算法从状态 $x$ 到 $x'$ 的转移与从状态 $x'$ 到 $x$ 的转移在稳态下是平衡的。我们还需要确保整个链是可逆的,并且其稳态分布是目标分布 $pi(cdot)$。
可约性和非周期性:保证收敛
尽管细致平衡条件保证了我们选择的目标分布 $pi$ 是一个稳态分布,但要证明它是唯一的稳态分布,并且链会收敛到它,还需要满足链的不可约性和非周期性。
不可约性 (Irreducibility):
要保证MH算法的链是不可约的,我们需要确保对于任意两个状态 $x$ 和 $x'$,从 $x$ 可以通过有限步转移到达 $x'$,且概率大于零。
这通常依赖于提议分布 $q(cdot|cdot)$ 的选择。如果提议分布是全空间支持的,也就是说,对于任意的 $x$ 和 $x'$,都有 $q(x'|x) > 0$,那么MH算法生成的链就是不可约的。
例如,高斯分布提议机制(步长可调)在连续状态空间中通常是全空间支持的。在离散状态空间中,如果提议分布能够覆盖所有可能的下一个状态,也能实现不可约性。
为什么不可约性很重要?如果链不是不可约的,那么它可能会被分割成几个独立的子集,从一个子集中的状态可能永远无法到达另一个子集中的状态。这样一来,我们从一个子集开始的链将永远无法探索到全局的概率分布。
非周期性 (Aperiodicity):
周期性是指马尔可夫链在转移过程中可能陷入一个长度为 $d > 1$ 的循环。例如,一个状态只能在偶数步后回到自身。
要保证MH算法的链是非周期的,通常只要链不是完全确定性的(即,存在非零概率的自循环或转移到其他状态),并且提议分布允许有一定的“随机性”,就可以实现。
例如,即使提议分布 $q(x'|x)$ 允许从 $x$ 提议到 $x$,并接受(转移到 $x'$),如果提议分布 $q(x|x')$ 也允许从 $x'$ 提议到 $x$,那么这个转移就不一定是确定的。
具体来说,如果存在一个状态 $x$ 使得 $q(x|x) = 1$ 并且接受率 $alpha(x|x) = 1$,那么这个状态就可能成为一个周期性循环的一部分。但是,MH算法的接受率允许我们保持当前状态(即 $x^{(t+1)} = x^{(t)}$),这种“自循环”的发生概率由 $1 sum_{x'
eq x} q(x'|x) alpha(x'|x)$ 决定。只要这个概率不为零,并且在链的某个点上我们有机会从一个状态 $x$ 转移到另一个状态 $x'$,我们就可以打破周期性。
许多常见的提议分布,如高斯随机游走,即使在离散化或某些特殊情况下可能表现出周期性,但通过允许采样到当前状态本身(即提议 $x' = x$),通常可以确保非周期性。或者,更简单地说,只要存在一条从状态 $x$ 到达状态 $x$ 的路径,其长度不一定是固定的,就可能打破周期性。
结论:迈向稳态
一旦MH算法的马尔可夫链被证明是不可约且非周期性的,并且我们已经证明了它满足细致平衡条件,那么根据马尔可夫链的遍历性理论,我们可以得出以下结论:
1. 唯一稳态分布: 存在一个唯一的稳态分布 $pi_{steady}$,它满足 $pi_{steady} P = pi_{steady}$,其中 $P$ 是转移矩阵。
2. 收敛性: 对于任何初始分布,该马尔可夫链都将收敛到这个唯一的稳态分布 $pi_{steady}$。
3. 目标分布: 由于细致平衡条件 $pi(x) P_{x, x'} = pi(x') P_{x', x}$ 是在假设 $pi$ 是稳态分布的情况下推导出来的,并且MH算法正是按照这个规则构建转移概率的,这就意味着我们构造的这个唯一的稳态分布 $pi_{steady}$ 就是我们最初的目标分布 $pi$。
换句话说,MH算法通过精心设计的提议和接受机制,构建了一个满足细致平衡条件的马尔可夫链。只要提议分布保证了链的不可约性和非周期性,那么这个链的稳态分布就必然是我们想要抽样的目标分布 $pi$。从这个稳态分布中抽取的样本,就代表了目标分布的特征。
因此,MH算法之所以能够达到马尔可夫稳态,并且该稳态是我们的目标分布,其根本在于 细致平衡条件的满足,以及提议分布所赋予的 不可约性 和 非周期性。这是理解MH算法强大之处的关键所在。