MCMC 方法是从特定分布中生成随机数的方法.  不同于之前说过的逆变换采样,  MCMC 方法可以适用于多维度分布,  并且不存在过多额外的数学过程.

Markov Chain Monte Carlo 方法,  顾名思义是一种依赖 Markov Chain (马尔科夫链) 的随机采样方法,  所以首先有必要来认识一下马尔科夫链.

马尔科夫链是指从一个状态转移到另一个状态的随机过程,  并且转移概率只与当前状态有关,  不与之前状态有关.  对于这篇专栏来说了解到这样子就可以了,  详细可以参考wiki.

至于 Monte Cario (蒙特卡洛) 涉及到积分的计算,  然而这里不会计算积分所以就忽略掉了.

关于符号,  这篇专栏里统一使用偏向数学风格的符合记述,  所以向量默认列向量;  P(x=x₁) 表示随即状态 x 为 x₁ 的概率,  简记为 P(x₁);  P(x=x₁|y=y₁) 表示 x 为 y₁ 时 x 为 x₁的概率,  简记为 P(x₁|y₁);  P(x=x₁,y=y₁) 表示 x 为 x₁ 并且 y 为 y₁ 的概率,  简记为 P(x₁,y₁); etc.

随机状态的转移和分布

因为马尔科夫链里状态转移的概率只与当前状态有关,  那么状态 i 转移到状态 j 的概率可以表示为 ,  可以使用 转移矩阵(transition matrix) 表示 .

使用 Dₜ 表示在随机状态在时刻 t 的分布(distibution),  Dₜ(i) 表示随机状态在时刻 t 取得状态 i 的概率,  那么下一时刻的分布 Dₜ₊₁ 可以由  得出,  或者写为向量形式: .

当有一个特定的分布 π,  经过状态转移后仍然符合这个分布,  即 ,  那么则这个分布称为平稳分布(stationary distribution).  以向量形式来说为 ,  即 π 为转移矩阵 T 的特征向量,  并且特征值为1.

如果转移矩阵 T 不可约并且是非周期的,  那么以任意分布 D₀ 开始,  执行无限次转移后,  都会得到平稳分布,  即 ,  这个结论由 Perron–Frobenius 理论给出.  下面写一个小程序对这个结论进行粗略的验证:

可以看到只要 T 不改变,  无论是什么起始分布,  最终都会收敛到 [0.625  0.3125 0.0625]ᵀ.

既然无论任何分布,  施加多次转移矩阵 T 都会使得分布趋向 T 对应的平稳分布 π,  那么可以使用任意起始状态的马尔科夫链作为对平稳分布的采样.  但需要注意的是,  马尔科夫链起始几项不太符合平稳分布,  所以采样需要在分布趋向平稳后开始.  下面的以上面的转移矩阵作的例子

可以看到采样分布比较符合平稳分布.

MH 采样

在一般情况下,  我们是需要对某特定分布进行采样,  而不是从转移里得到分布,  所以有必要从特定分布里求出转移.

当分布 π 为转移 T 的平稳分布时,  某状态 i 的概率密度流出量等于流入量,  即满足等式 ,  满足这个等式的一个充分条件为 ,  后一条件被称为细致平衡(detailed balance)条件.

对于特定分布 π,  选取任意转移 Q,  则存在 α 使得  成立.  不难知道,   是一种方便的取法.  于是转移  的平稳分布即为 π.

使用接受-拒绝(acceptance-rejection)采样的思路,  生成相应的马尔科夫链步骤如下:  1) 选取随机起始状态 x₀;  2) 进入循环 t = 1, 2, ... :  2.a) 从条件分布 Q(x'|xₜ₋₁) 中得到候选状态 x';  2.b) 计算转移接受概率 α = Q(xₜ₋₁|x') π(x');  2.c) 选取分布在区间 [0,1) 的均匀随机数 s;  2.d) 如果 s < α, 那么接受转移 xₜ = x',  否则不接受转移 xₜ = xₜ₋₁.  3) 最后截断起始还没符合分布 π 的链,  得到可以代表分布 π 的采样集.

[虽然这里可以插入一个例子,  但是我真的很懒]

上述算法有两个很严重的问题:  看到接受概率 α,  在离散情况时这是两个小于 1 的数字相乘,  这意味着接受率很小,  从而导致需要很长的链才能得到足够符合分布的采样.  另外,  当状态不是离散而是连续时,  描述状态分布由向量变为函数,  即概率密度函数(probability density function),  根据 pdf 的定义,  pdf 的值允许大于1,  所以接受概率 α 可能大于1,  大于 1 的概率是没有意义的.

为此,  引入变量 b 使得 ,  可以看到无论 b 取何值,  细致平衡条件都是成立的.  为了解决在离散时 α 过小以及在连续时 α 可能无意义两个问题,  选取 b 使得 ,  不难得出 ,  又记 ,  那么 β 即是修改后的转移接受概率.  这种算法被称为 Metropolis-Hastings 算法.

其他 MCMC 方法都可以看作是 MH 算法的特例.  如果转移 Q 的对称的,  即 ,  MH 算法退化为 Metropolis 算法.

关于 Q 的选取一般有两种:  1) 直接选取一个与目标分布接近且容易生成的分布;  2) 以当前状态为中心向附近跳跃一段距离 (正态分布是好文明).

下面以第二种选取 Q 的方法作为例子:

可以看到生成的马尔科夫链分布是非常接近目标分布的,  并且转移接受率有 0.64.

不难知道,  MH 采样极度依赖于转移 Q,  如果 Q 选取不当,  接受转移率降低会导致:  1) 马尔科夫链需要抛弃起始极长一段;  2) 生成的随机采样质量差;  3) 计算成本增大.

选取合适的转移最简单就是观察目标分布的特征,  然后选取一个接近且容易计算的分布.  以上面展示的分布为例,  转移取 "均值=0, 标准差=7.5"的正态分布的右半部分,  转移接受率可以达到 0.47.  对于 MH 算法来说转移接受率大于 0.2 ~ 0.45 就已经可以使用了,  具体还是要看应用场景.

另外还有一个问题,  跟离散情况一样,  马尔科夫链开头几项认为是不符合目标分布,  所以需要抛弃的,  但是究竟需要抛弃多长的链?  可以同时运行许多许多不同起始状态的马尔科夫链,  抽取同一时刻的项跟目标分布进行对比,  这样可以大概地确定需要抛弃的数目.

高维采样

需要注意的是,  "状态" 并不是一个限制在一维的数字,  实际上就像数学上定义 "集合" 为一堆东西组成的对象,  并没有对 "东西" 作定义,  而 "状态" 可以理解为集合里的元素,  所以状态实际上也是任意的.  以下篇幅主要讨论一下二维连续分布,  更高维度的也可以作类比.

下面用 MH 采样对分布  进行采样,  其中 N 是归一化因子,  众所周知归一化因子不重要好吧.

可以看到采样的分布比较符合目标分布.  但是有一个问题,  转移接受率只有 0.27,  实际上这个数字对于高维的 MH 采样来说已经算不错的了,  普遍来说转移接受率都在 0.1 左右.  对此,  这里提出 Gibbs 采样.

在二维空间里进行采样,  采样步骤可以分解为:  1) 先对 x 分量进行采样,  2) 在取得 x 后再对 y 进行采样,  即是 .  但 x, y 的采样顺序是可以对换的,  即有 .  实际上联合上面两式可以得到贝叶斯定理:  .

考虑只有一个坐标不相同的两点,  比如说 x 不一样:  (x₁,y) 和 (x₂,y),  根据上面提出的采样分解步骤,  不难知道等式  成立,  对比细致平衡条件,  可以知道在 x 方向上转移矩阵可以为 ,  类似地 y 方向也有相同的结论:  .

于是 Gibbs 采样步骤为:  (起始状态和循环什么的就不写了)  1) 已知状态 (xₜ,yₜ);  2) 在 x 方向上根据 Pr(x|yₜ) 采样得到 xₜ₊₁;  3) 在 y 方向上根据 Pr(y|xₜ₊₁) 采样得到 yₜ₊₁;  4) 得到下一个状态 (xₜ₊₁,yₜ₊₁).

在有些地方见到的 Gibbs 采样在细节上会有不一样的地方,  比如说有些地方会把中间产生的状态 (xₜ₊₁,yₜ) 也当作一个合理的采样;  还有些地方不是坐标轴轮换采样,  而是随机选取一个坐标轴进行采样.  这些变化都是没问题的,  完全ok.

以上面的分布  为例,  可以求得条件分布: ,  其中  分别是与 x 无关的归一化因子,  可以看到 x 对于 y 的条件分布是 "均值为 , 标准差为 " 的正态分布,  y 对于 x 的分布也有相同的结论.  于是得到 Gibbs 采样的程序:

可以看到在 Gibbs 采样里,  因为接受-拒绝环节被取消了,  生成的样本比 MH 采样更加均匀.  实际上如果把中途生成的状态 (xₜ₊₁,yₜ) 也当作采样,  就可以看到状态是如何转移的:

Gbbis 采样相比与 MH 采样,  优势就是以 100% 概率接受状态转移,  并且不需要知道多个变量之间的联合分布,  只需要知道单个变量对于其他变量的条件分布,  在一些情况下联合分布比条件分布难求得多.

不过 Gibbs 采样也不是万能的,  如果条件分布也是复杂分布的话,  虽然可以使用简单的接受-拒绝采样,  但是从成本上来说说不定比 MH 采样更差.  对于条件分布极难甚至不能求出的情况,  Gbbis 采样就完全不适用,  这时候可以采用更高级的采样方法,  这里举几个例子:  自适应 MH 采样,  Hamiltonian 方法,  Langevin 规则,  切片采样,  可逆跳跃采样.

采样方法多种多样,  不过还是那句话

所以更高级的采样方法我也懒得去翻了,  这篇专栏只是学量子力学途中遇上一点难题,  才顺便拉出来的.

日常推销瑟图群:  274767696

封面pid: 66081482