蒙特卡洛-马尔科夫链(MCMC)初步

。

解决什么问题？
我们常常遇到这样的问题：模型构建好之后，有一个概率p(x)p(x)（称为目标分布），不能显式的给出其表达，只能生成一系列符合这个分布的xx。这种问题称为“采样”。

特别地，在贝叶斯方法中，关注的是后验概率p(x|D)p(x|D)。在给定观测DD的情况下，需要估计系统参数xx。如果后验概率没有明确表达，或者由于多重积分难以计算，则无法直接求解xx，但是可以生成一系列符合此概率的分布。
常见的情况是，后验概率难以明确表达。在已知观测数据x的情况下，利用p(x|D)=p(x)⋅p(D|x)/p(D)∝p(D|x)p(x|D)=p(x)⋅p(D|x)/p(D)∝p(D|x)。对似然函数进行采样。

名字解析
Monte Carlo方法：用多次随机求平均的方法来逼近一个值。实际是采样方法的核心

设计MCMC方法的一个难处，在于如何设计合理的转移概率函数，使得MC链的稳态分布等于要求的概率分布。
MCMC的另一个问题是混合速度，在这篇博客中有所讨论。

MCMC是一个大类，有许多种具体算法，以下举例几种最为著名的。

Gibbs Sampling
Gibbs Sampling处理这样的问题：对于一个高维随机变量x=[x1,x2,x3]x=[x1,x2,x3]，不能写出其各分量的联合概率p(x)p(x)，但是可以写出各个分量之间的条件概率。

首先任取初始值x0x0。在已知当前采样xtxt时，按照如下方法生成t+1时刻采样。

根据x2x2,x3x3采样x1x1：x1t+1∝p(x1|x2t,x3t)xt+11∝p(x1|xt2,xt3)
根据x1x1,x3x3采样x2x2：x2t+1∝p(x2|x1t+1,x3t)xt+12∝p(x2|xt+11,xt3)
根据x1x1,x2x2采样x3x3：x3t+1∝p(x3|x1t+1,x2t+1)xt+13∝p(x3|xt+11,xt+12)
每一次采样，都尽可能利用其他分量的最新结果。Gibbs Sampling特别适用于贝叶斯网络的采样，因为该网络本身就包含一系列条件概率。

Metropolis-Hastings(MH)
MH方法的一个优点是，即使不能写出概率密度函数p(x)p(x)，可以用一个和其成正比的函数f(x)f(x)来采样。如前所述，这在贝叶斯方法中非常方便。
首先任取初始值x0x0。在已知当前采样xtxt时，按照如下方法生成t+1时刻采样。

根据一个概率分布Q(xt+1|xt)Q(xt+1|xt)生成一个候选采样xt+1xt+1。
其中QQ称为proposal density或者jumping distribution。可以取以xtxt为中心的高斯分布。
比较新旧采样的概率密度函数
如果f(xt+1)>f(xt)f(xt+1)>f(xt)，接受xt+1xt+1为新采样；
否则，按照p=f(xt+1)/f(xt)p=f(xt+1)/f(xt)选择xt+1xt+1为新采样，或者维持xtxt为新采样。
一种常见的情况，如果待求概率是一个后验p(x|D)p(x|D)，且其先验是高斯，则可以做如下变换：
p(x|D)∝p(x)⋅p(D|x)=N(x;0,Σ)⋅p(D|x)
p(x|D)∝p(x)⋅p(D|x)=N(x;0,Σ)⋅p(D|x)

直接使用以下方法生成候选样本
xt+1=1−ϵ2−−−−−√xt+ϵν
xt+1=1−ϵ2xt+ϵν
其中ν∼N(0,Σ)ν∼N(0,Σ)，ϵ∈[−1,+1]ϵ∈[−1,+1]是扰动步长。新样本等于当前样本和先验的加权和。另一种表达方法是
xt+1=cosθ⋅xt+sinθ⋅ν
xt+1=cos⁡θ⋅xt+sin⁡θ⋅ν
不同的ϵϵ对应的新采样的轨迹是半个椭圆圆周。椭圆的两轴分别为xtxt和νν。如果ϵ=0ϵ=0，则新采样和旧采样相同（红色）。ϵ或者θϵ或者θ控制扰动幅度。