贝叶斯方法和采样技术结合-MCMC-后验分布

贝叶斯方法与采样技术的结合

贝叶斯方法和核心就是计算某个未知参数的后验分布,当参数的后验分布知道了，就可以利用分布获得你想要知道的信息。在这里假设观测数据为$Y$,未知参数为$theta$,似然函数$f(y|theta)$,参数的先验分布为$p(theta)$

目标：估计$theta$
主要想法：求$p(theta|y)$，由于$$p(theta|y)=\frac{p(theta)p(y|theta)}{p(y)}$$
分母是观测数据的边际概率，设计到了求和或者求积分，而且分母部分是一个常数，所以我们可以认为$p(theta|y)$正比于$p(theta)p(y|theta)$,我们把$p(theta)p(y|theta)$看作未归一化的未知参数的后验分布，记为$\breve{p(theta|y)}$。
我们一般就基于这个为归一化的后验分母进行一些推断，把它也叫做目标分布，那我当初是非常不理解的，我 以前认为那这样的话，分布不就变了吗，但是你要记住它是一个常数，一个确定的值，所以 不管参数是什么，$p(y)$都是一样的，所以你会发现$$\frac{p(thet{a}^{\ast }|y)}{p(theta|y)}=\frac{\breve{p(thet{a}^{\ast }|y)}}{\breve{p(theta|y)}}$$
这个比值与$p(y)$没有任何的关系，所以我们在进行采样技术时，获得不了归一化的目标分布也是没有关系的，我们可以先采样获得来自未归一化的目标分布的样本，然后再求分布，因为分布就是求积分或者求和，求积分或者求和可以转化成求期望，求期望我们就可以利用大数定律的思想。

我之前的博客当中有非常详细的介绍过近似推断方法，近似推断法包括直接采样法和间接采样法，其中直接采样法我主要介绍的是利用均匀分布获得来自目标分布的样本，间接采样法里面我讨论了拒绝采样，重要性采样，MH算法和吉布斯采样算法。

拒绝采样法–核心思想

我们想要获得来自目标分布$p(x)$的样本，在这里我们就是想要获得来自目标分布$\breve{p(theta|y)}$的样本，如果这个分布比较复杂没办法直接抽样，那么怎么办呢，我们引入一个好抽样的且接近$p(x)$的分布$q(x)$，把这个容易抽样的分布叫做建议分布，使得存在一个常数$k$,对于任意从建议分布中抽取的样本$x^{\ast }$,都有$kq(x^{\ast })>p(x^{\ast })$，那么令$\alpha =\frac{p(x^{\ast })}{kq(x^{\ast })}$，$\alpha$为抽取的样本$x^x$被接受的比例，我们从均匀分布$u(0,1)$中任意取一个随机数$u$,当$u<\alpha$则表明$x^{\ast }$被接受了，被接受是来自目标分布$p(x)$的样本了，再对$p(x)$作一些推理时，$x^{\ast }$就被考虑在内了，否则不能把$x^{\ast }$看做来自目标分布的样本

注意：这种方法需要确定k,k的选取很重要，如果k太大，那么必然拒绝的比例会增大，那么如果想要获得足够的来自目标分布的样本就要增加采样次数，而且在高维空间中，值接近，分布可能差很多

重要性采样–核心思想
重要性采样算法的思想是根据拒绝采样来的，拒绝采样那里我们的建议分布要是一个接近目标分布的分布，但是当我们的最终目的是估计一个函数$f(x)$在分布$p(x)$下的期望，我们其实不必要求最终获得的样本一定是来自目标分布$p(x)$，主要是因为$f(x)p(x)=q(x)\frac{p(x)}{q(x)}f(x)$,我们可以令$\frac{p(x)}{q(x)}=w(x)$,$w(x)$为重要性权重，所以$E_{p(x)}f(x)=E_{q(x)}w(x)f(x)$,所以同样大数定律，从$q(x)$中抽取样本点，估计$wf(x)$的期望即可

最开始我说了当我们只能知道未归一化的目标分布$\breve{p(x)}$时,然后你可能会好奇我们采用拒绝采样法获得的接受的样本$x^{\ast }$它是服从未归一化的目标的分布的把，如果我们只是目标是求期望的话，我们自然利用重要性采样的思想求
$$E_{p(x)}f(x)={\int }_{x}f(x)\frac{\breve{p(x)}}{Z}dx=\frac{{\int }_{x}f(x)\breve{p(x)}dx}{{\int }_x\breve{p(x)}dx}$$
而对于分子和分母就可以利用建议分布$q(x)$抽取样本，利用重要性权重求期望。

那如果我们是想要得到来自归一化的目标分布的样本时，在拒绝采样当中，我们首先获得的是来自未归一化的目标分布的样本，然后你要记得我们为什么不好获得归一化的目标分布，首先是因为分母不好求，分母涉及到了边缘化的思想（会进行多重积分），那么如果我们获得了未归一化的样本，那么分母是未归一化函数的定积分，既然是定积分那么就可以利用蒙塔卡洛的思想进行抽样近似计算，样本也知道了，就相当于求$f(x)=1$这个函数的期望。所以利用来自未归一化分布的样本可以计算得到归一化的目标分布。

然后你要知道在拒绝抽样当中$q(x)$是接近$\breve{p(x)}$的，那么现在$\frac{q(x)}{z}$是接近$\frac{\breve{p(x)}}{z}$的，你很容易知道新的建议分布中抽取 的样本与原来的建议分布的样本相比，只是一个倍数的问题，也就是归一化目标分布的抽取的样本值是未归一化抽取的样本值的$\frac{1}{z}$倍。

所以进行总结：

• 当我们已知归一化的目标分布时，我们利用拒绝采样接受的样本经过部分处理之后应该是服从未归一化的目标分布的
• 当我们已知的是未归一化的目标分布时，同样我们利用拒绝采样获得的样本再进行一些处理，最后得到的样本是服从未归一化的目标分布的，那么如果我们确实想要获得归一化的目标分布的样本，那就利用取得的样本计算归一化常数$z$，然后再拿之前的样本除以$z$即可。

上面讲解的都是简单的近似采样的方法，虽然简单但是其思想是非常重要的，近似的思想，求期望求定积分的思想都是需要掌握的

下面介绍更加常用和更加符合实际情况的采样方法

MCMC算法

MCMC算法我在上次博客介绍的很多了，这里再来做个概括，回想一下我们在拒绝采样当中是不是有一个建议分布，从建议分布中抽取样本，但是这里的建议分布我们是说接近目标分布且容易抽样的分布，t+1时刻抽取的样本与t时刻的是没有关系的，都是从给定的这个建议分布中抽样，那么MCMC（蒙塔卡洛马尔科夫链）的建议分布就做出了改进，我们假设建议分布是服从马尔科夫性质，也就是从建议分布中抽出来的样本是一个马尔科夫链

马尔科夫链

马尔科夫链也就是当状态转移分布确定时，t+1时刻抽取出来的样本与t时刻抽取的样本有关，而与其他时刻抽取出来的样本都没有关系，也就是 说假设我们观测的一个人的心情，我今天的心情只与昨天的心情有关，和从昨天的心情转移到今天的心情的概率。关于马尔科夫链的详细讲解大家可以看李航的统计学习方法那本书。

平稳分布

我们的目标是想要获得来自目标分布的样本，那么根据马尔科夫链的性质，对于不可约的非周期且正常返的马尔科夫链有唯一的平稳分布。平稳分布是假设初始的状态分布为${\pi }_0$,拿我的心情来举例子，给定状态转移矩阵$p$，其中$p_{ij}$表示$p(i_{t+1}=q_j|i_t=q_i)$。今天我开心，那么明天我的心情服从的分布为$p{\pi }_0$，注意我现在已经知道了今天的心情，那么我今天的心情服从的分布就是只有开心对应的元素为1,别的心情对应的元素都为0，我明天的心情服从的分布是$p{\pi }_0$，这是一个列向量，也就是我明天的心情服从这样的一个离散的状态分布，那么我明天的心情就是从这个离散的状态分布进行随机抽样。那么后天我的心情服从的分布就为$p^2{\pi }_0$,一直是这样的过程，直到$p^n{\pi }_0$不发生变化，即是直到有一天当期服从的分布为$\pi$，下一期我的心情服从的分布按理说应该为$p\pi$，但是此时$p\pi =\pi$,那么就说明了此时我的心情服从的分布是一个平稳的分布，就把满足$p\pi =\pi$的$\pi$叫做我的心情的平稳分布

所以在这里我们的想法就是构造一个马尔科夫链，使得这个马尔科夫链的平稳分布就是我们的目标分布，那么关键也就是找到状态转移分布$q(x)$,当变量是离散型随机变量时，我们想找到一个状态转移概率矩阵，当变量是连续型随机变量时，我们想要找到状态转移核函数。

细致平稳方程

如果状态分布为$\pi$,其中状态$i$的概率为${\pi }_i$,状态$j$的概率为${\pi }_j$,状态转移分布为$p$,如果有$p_{ij}{\pi }_i=p_{ji}{\pi }_j$成立，那么我们称$\pi$是这个马尔科夫链的平稳分布，当得到平稳分布时，继续按照这种模式今年采样，那么可以认为状态平稳时以后获得的样本就是来自平稳分布的。

MH算法

所以我们想要找到$q(x,x^{\prime })$让它满足$q(x,x^{\prime })p(x)=q(x^{\prime },x)p(x^{\prime })$,但是实际直接找出一个这样的$q(x,x^{\prime })$是非常难的，想想也知道挺难的，所以我们给出一个试错的过程，我们不必要一定要找一个$q(x,x^{\prime })$，使得根据这个状态转移分布得到的马尔科夫链的平稳分布就是我们的目标分布$p(x)$,我们利用拒绝采样的思想，我们利用一个好抽样的状态转移分布$q(x,x^{\prime })$，抽取样本，但是这个样本需要一个接受率进而和来自均匀分布的随机数$u$进行比较决定接不接受这个样本，具体来说$$\alpha (x,x^{\prime })=min(1,\frac{p(x^{\prime })q(x^{\prime },x)}{p(x)q(x,x^{\prime })})$$
那么真正的状态转移分布为$q(x,x^{\prime }{)}^{{}^{\prime }}=\alpha (x,x^{\prime })q(x,x^{\prime })$，实际上是根据$q(x,x^{\prime }{)}^{{}^{\prime }}$进行抽取样本的

注意：这里的不接受不代表舍弃这个样本，不代表这个样本不是来自目标分布，在MH算法当中，不接受只是表明不更新，t+1时刻的状态仍然等于t时刻的状态。不要与单纯的拒绝采样混淆了。

上面说的最主要的就是状态转移分布$q(x,x^{\prime })$,它不同，MH算法又可以分为不同的形式。推断算法博客里面写过了,最常用的是设置$q(x,x^{\prime })=q(x‘,x)$，也就是设置对称的建议分布

• 最简答的一种是$x^{\prime }=x+u(-a,a)$,下一时刻的状态取值是在上一时刻的状态取值加上一个来自对称的均匀分布样本值。
• 其次就是$q(x,x^{\prime })$取条件概率分布$p(x^{\prime }|x)$，但是注意这个$p$当然不是我们的目标分布，只是表示的是条件概率分布，最常用的是当$x$已知时，$x^{\prime }$服从均值为$x$,协方差矩阵是常数矩阵的多元正态分布，即$x^{\prime }服从N(x,\sum )$
• 还有一种随机游走Metropolis算法，这个跟第一种的原理一样，只是取的游走的形式不一样，mcmc程序包里有metrop()函数使用的是标准正态变量的随机游走，即$x^{\prime }=x+scale\ast z$,z是来自标准正态分布的随机变量，sclae是一个尺度，是常数。

你可能同样会思考这个问题，当我们不知道归一化的目标分布的概率密度函数形式，怎么在MH算法中得到来自目标分布的样本呢？这里就跟之前说的方法不一样了，因为我们这里虽然最初不知道归一化的目标分布$p(x)$,但是我们的$$p(x^{\prime })q(x|x^{\prime })\alpha (x^{\prime },x)=p(x)q(x^{\prime }|x)\alpha (x,x^{\prime })$$
$\breve{p(x^{\prime })}q(x|x^{\prime })\alpha (x^{\prime },x)=\breve{p(x)}q(x^{\prime }|x)\alpha (x,x^{\prime })$,
所以我们抽取的样本就是来自归一化的目标分布的

其实归一化与未归一化相比也就是只相差一个常数而已。（这里我不知道理解的对不对：所以我们之前在拒绝采样里面讨论的还要除以这个常数我现在感觉没有必要了，大家都没有除以这个比例，因为对应的建议分布也是一个常数比例关系，所以我把建议分布直接看成接近归一化目标分布的容易抽要的分布也是可以的，我们就可以把抽取的样本经过处理之后看成来自归一化的目标分布的。）

小提醒：
在MH算法，我们经常使用对数化的目标分布，为了额编程时把乘除号变成加减号，我们一般都会进行取对数处理

现在来讨论最后一种采样方法

Gibbs采样

这里简单说一句就是gibbs采样，也就是吉布斯采样算法，它也是MCMC采样的一种，最初思想一样。都是为了找到$q(x,x^{\prime })p(x)=q(x^{\prime },x)p(x^{\prime })$，MH算法和Gibbs算法都是为了取得这个结果，采取的不同形式。

在吉布斯采样把满条件分布看作是建议分布，**注意这个满条件分布是目标分布的满条件分布，**不要与单分量的metroplolis算法搞混了，后者按照建议分布的满条件分布进行逐个维度上的采样，而前者建议分布直接就是目标分布的满条件分布，满条件分布就是在其他维度上的取值不变的条件下，每次对一个维度上进行采样。

当对下一个维度进行按照满条件分布进行采样时，之前采样过的维度值变成最新的。
Gibbs采样没有拒绝更新的这个过程，也就是接受概率为1

以上就是贝叶斯方法与采样算法的结合