采样算法:马尔科夫链中的各类蒙特卡洛采样算法「2.1」

最近学习了机器学习中的马尔科夫链蒙特卡洛(Markov Chain Monte Carlo, 简称MCMC) 相关的知识。

主要内容包括：

【1】蒙特卡洛原则，及其应用于采样的必要性

【2】用于求解最大似然、近似推断、期望问题的经典采样算法：Metropolis-Hastings，Rejection，Importan，Metropolis和Gibbs算法。

【3】马尔可夫链各个性质在蒙特卡洛采样问题中的应用，包括同质性，平移不变性

—————【2】—————

采样过程的本质，可以理解为实现目的【已知目标函数的形状，不能求得函数的真实形式，但又要获取服从此目标函数分布的样本集合】而采用的获取样本的方法。我们的目标，是一个集合{x1,...,xN}。

【拒绝采样算法 Rejection】

假设要从分布p(z)中采样，p(z)不能直接求，但根据贝叶斯公式分解为 P_(z)*(1/Z), 这个假设即意为：分子项P_(z)类似 P(y|x)*P(x)，为已知；分母项为 P(y)=积分x(P(y|x)*P(x) dx), 分母项为不可直接求的未知正态化常量。

拒绝采样算法的思路是这样的：

下图中，红色曲线及其覆盖的白色区域，是我们要求的目标概率密度函数，即目标分布，同时也是上面提到的P_(z)，蓝色曲线及其与红色曲线之间的灰色区域，是我们假设出来的一个函数q(z)，它的性质，我们要求它乘以常数k（某个一定存在的常数）后，在z的定义域上每一点，都使得 k*q(z) >= P_(z) 成立。

由于这里的q(z)即途中蓝色曲线是我们所假定的某个函数，我们要求它具有易于完成采样操作的性质，例如，这里我们假设它是高斯分布。而目标函数红色曲线P_(z),我们知道它在整个定义域上对每个z的函数值，但不知道它的具体函数式。

因此接下来重复N次，获取N个样本：

1、在【0，1】的区间内均匀随机取一个小数u

2、从蓝色曲线的假定分布q(z)随机采样zi和q(zi)

3、检验 k*q(zi)*u > p_(z)，则拒绝；否则保留此样本作为最终采样结果集合的一个元素

这样就得到了满足真实分布p(z)的计量分布p_(z)的样本集合，即使未求出此分布的具体表达式。

拒绝采样算法使用的假设分布这样取是为了优化速度，如果使用更大范围的q(z)比如矩形，最终结果也不会改变，只是运行更慢。

它的问题在于：1、对整个定义域z，这个假设分布q(z)，不一定存在，因为它需要满足kq(z) > p_(z) 的性质，若k非常大才能满足之，则每个随机采样的样本被接受的概率太小，在高维情况下此方法会难以运算。

【重要性采样算法 Importance】

重要性采样的思路，同样假设一个易于采样的分布q(z)，但对其样本不是保留一部分，而是全部保留，但使用重要性参数调整其权重，来求目标概率密度函数的期望，此期望由蒙特卡洛原则证明趋近于真实目标函数的概率密度（见【1】）。

目标概率密度函数f(z)的期望，通过引入假设概率q(z)，转化为对 （f(z)*p/q）的期望，根据蒙特卡洛原则转化为L个样本的重要性加权函数值的和。重要性系数为p(zi)/q(zi)，记为ri。这些重要性系数修正了从假设分布q(z)中采样产生的误差，保留了所有采样样本。

通常，目标函数的后验估计P_(z)容易估计，但正态化系数Zp未知。（【1】中已经讨论），同样的方式定义假设分布 q(z)=q_(z)/Zq ,带入上面的推导，则 E[f]= Zq/Zp *(1/L)* sum(r_*f(z)), 这里 r_= p_(z)/q_(z)。

因此，求 Zp/Zq= P_(z)对所有z求积分 /Zq= 1/L * sum(r_),带入求得 E[f]=sum(w*f(z)),其中w为重要性权重w= r_i / sum(r_L)，r_=p_(z)/q_(z)可以轻松求得（p_,q_可观测）。

重要性采样算法的不足：

1、算法成功依赖于假设分布q(z)和p(z)的相似性

2、快速波动的p(z)*f(z)会导致大量样在本在z上集中在小区域

3、重要性权重r可能被一小部分值很大的权重主宰

这些因素的存在会导致采样结果失真于真实分布。

【在贝叶斯网络中应用重要性采样】

目标是求P(xf | xe)分布的样本,即为图中概率。

假设不能直接求此条件概率，使用重要性采样获取满足此条件概率的样本。首先将条件概率写成 p(x1,x4,x5,x2=1,x3=1)/p(x2=1,x3=1)=p(xf,xe)/p(xe)=1/Zp * p_(x)，其中p_(x)为未正态化的目标分布。

假设分布，可以设为q(x1,x4,x5)=p(x1)p(x4|x3=1)p(x5|x2=1),对此假设分布采样：

从p(x1)中采样x1，从p(x4|x3=1)中采样x4，从假设分布的边缘分布采样x1，x4，x5。边缘分布概率图中已经给出。

接下来求每个样本的重要性权重，wi= r_i / sum(r_L)。例如样本{x1=0,x4=1,x5=1},求出p_(xi)=0.006912,直接根据贝叶斯图给出的关系和概率值计算。q(xi)=0.288，根据我们之前假设的假设分布的概率计算关系q(x1,x4,x5)=p(x1)p(x4|x3=1)p(x5|x2=1。这样就求出了每个样本的权重。

最后，我们真正要求的条件概率 p(x1,x4,x5|x2=1,x3=1) 就是 【满足x1，x4，x5某个条件状态的所有样本权重之和】/【所有状态样本权重之和】。这样就证明了我们采集的样本服从于我们的目标分布。

回顾整个采样过程，避开了直接求条件概率，做了以下几件事：

1、将目标条件概率转换为 联合概率（易求）/ 全概率（回避），联合概率已知

2、应用蒙特卡洛原则，将期望公式变形为假设分布、联合概率相关的权重值公式

3、假设分布已知，联合概率已知，则可求出每个样本的权重值。

4、执行采样，对假设分布均匀采样，经过上述处理求出目标条件概率的数值。

且蒙特卡洛原则证明了随着样本量的扩大，此数值解逼近于真实积分值。