MCMC采样和M-H采样

本篇主要解决一个问题，对于给定的概率平稳分布 $\pi$, 如何找到与之对应的概率转移矩阵 $P$。这个问题的解决主要是两个方法。MCMC 采样 和 M-H 采样。

1: 马尔可夫链的细致平稳条件
目标是找到一个矩阵 $P$, 使得 $\pi P=\pi$，记此式为式子1。
上式的左端的第 $j$ 元素是，${\sum }_{i=0}^{\infty }{\pi }_i{P}_{i,j}$； 右端第 $j$ 元素是 ${\pi }_j$。
如果说，可以找到一个矩阵 $P$, 使得 ${\pi }_i{P}_{i,j}={\pi }_j{P}_{j,i}$，记此式为式子2。
则式子1成立。所以我们只需要找到一个满足式子2的矩阵P即可。
式子2被称为马尔可夫链的细致平衡条件。

2: MCMC 采样
由于一般情况下，目标平稳分布𝜋(𝑥)和某一个马尔科夫链状态转移矩阵𝑄不满足细致平稳条件，即 $\pi (i)Q_{i,j}\ne \pi (j)Q_{j,i}$。
那么怎么使得该式子成立？可以引进一个参数 $\alpha (i,j)$, 使得
$\pi (i)Q_{i,j}\alpha (i,j)=\pi (j)Q_{j,i}\alpha (j,i)$。记该式子为式子3。

只需要使得 $\alpha (i,j)=\pi (j)Q_{j,i}$, 那么，$\alpha (j,i)=\pi (i)Q_{i,j}$，从而式子3成立。

所以分布 $\pi (x)$ 对应的马尔可夫链转移矩阵 $P$, 满足
$P(i,j)=Q_{i,j}\alpha (i,j)$。

“”“也就是说，我们的目标矩阵𝑃可以通过任意一个马尔科夫链状态转移矩阵𝑄乘以𝛼(𝑖,𝑗)得到。𝛼(𝑖,𝑗)我们有一般称之为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵𝑃可以通过任意一个马尔科夫链状态转移矩阵𝑄以一定的接受率获得。这个很像我们在MCMC(一)蒙特卡罗方法第4节讲到的接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵𝑄通过一定的接受-拒绝概率得到目标转移矩阵𝑃,两者的解决问题思路是类似的。“”“

“”“好了，现在我们来总结下MCMC的采样过程。

1）输入我们任意选定的马尔科夫链状态转移矩阵𝑄，平稳分布𝜋(𝑥)，设定状态转移次数阈值𝑛1，需要的样本个数𝑛2
　　　　2）从任意简单概率分布采样得到初始状态值𝑥0
　　　　3）for 𝑡=0 to 𝑛1+𝑛2−1:

a) 从条件概率分布𝑄(𝑥|𝑥𝑡)中采样得到样本𝑥∗
　　　　　　b) 从均匀分布采样𝑢∼𝑢𝑛𝑖𝑓𝑜𝑟𝑚[0,1]
　　　　　　c) 如果𝑢<𝛼(𝑥𝑡,𝑥∗)=𝜋(𝑥∗)𝑄(𝑥∗,𝑥𝑡), 则接受转移𝑥𝑡→𝑥∗，即𝑥𝑡+1=𝑥∗
　　　　　　d) 否则不接受转移，即𝑥𝑡+1=𝑥𝑡
　　　　样本集(𝑥𝑛1,𝑥𝑛1+1,…,𝑥𝑛1+𝑛2−1)即为我们需要的平稳分布对应的样本集。

上面这个过程基本上就是MCMC采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于𝛼(𝑥𝑡,𝑥∗)可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个𝑛1要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的M-H采样出场了。”“”
　　　　
3: M-H采样
Metropolis-Hastings采样或M-H采样
M-H采样解决了我们上一节MCMC采样接受率过低的问题。
“”“
我们回到MCMC采样的细致平稳条件：
𝜋(𝑖)𝑄(𝑖,𝑗)𝛼(𝑖,𝑗)=𝜋(𝑗)𝑄(𝑗,𝑖)𝛼(𝑗,𝑖)
　　　　我们采样效率低的原因是𝛼(𝑖,𝑗)太小了，比如为0.1，而𝛼(𝑗,𝑖)为0.2。即：
𝜋(𝑖)𝑄(𝑖,𝑗)×0.1=𝜋(𝑗)𝑄(𝑗,𝑖)×0.2
　　　　这时我们可以看到，如果两边同时扩大五倍，接受率提高到了0.5，但是细致平稳条件却仍然是满足的，即：
𝜋(𝑖)𝑄(𝑖,𝑗)×0.5=𝜋(𝑗)𝑄(𝑗,𝑖)×1
　　　　这样我们的接受率可以做如下改进，即：
𝛼(𝑖,𝑗)=𝑚𝑖𝑛{𝜋(𝑗)𝑄(𝑗,𝑖)𝜋(𝑖)𝑄(𝑖,𝑗),1}
　　　　通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

1）输入我们任意选定的马尔科夫链状态转移矩阵𝑄，平稳分布𝜋(𝑥)，设定状态转移次数阈值𝑛1，需要的样本个数𝑛2
　　　　2）从任意简单概率分布采样得到初始状态值𝑥0
　　　　3）for 𝑡=0 to 𝑛1+𝑛2−1:

a) 从条件概率分布𝑄(𝑥|𝑥𝑡)中采样得到样本𝑥∗
　　　　　　b) 从均匀分布采样𝑢∼𝑢𝑛𝑖𝑓𝑜𝑟𝑚[0,1]
　　　　　　c) 如果𝑢<𝛼(𝑥𝑡,𝑥∗)=𝑚𝑖𝑛{𝜋(𝑗)𝑄(𝑗,𝑖)𝜋(𝑖)𝑄(𝑖,𝑗),1}, 则接受转移𝑥𝑡→𝑥∗，即𝑥𝑡+1=𝑥∗
　　　　　　d) 否则不接受转移，即𝑥𝑡+1=𝑥𝑡
　　　　样本集(𝑥𝑛1,𝑥𝑛1+1,…,𝑥𝑛1+𝑛2−1)即为我们需要的平稳分布对应的样本集。

很多时候，我们选择的马尔科夫链状态转移矩阵𝑄如果是对称的，即满足𝑄(𝑖,𝑗)=𝑄(𝑗,𝑖),这时我们的接受率可以进一步简化为：
𝛼(𝑖,𝑗)=𝑚𝑖𝑛{𝜋(𝑗)𝜋(𝑖),1}
“”“

4: M-H采样实例

略

5: M-H采样总结

(1): 效率低
(2): 如果只知道各个特征之间的条件概率分布，怎么采样？

Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下，下一篇我们就来讨论Gibbs采样。