AI算法 - 抽样方法

1 采样介绍

   在粒子滤波等算法中都提到采样，即从后验概率分布中采样，也有介绍抽样分布。从概率分布p(x) 中获取样本的过程叫采样，

大多数时候遇到的问题模型是采样的逆过程: 给一定数量的数据或样本，计算他们服从的概率分布。采样主要研究的是给定一
个概率分布p(x) ，如何生成满足条件的样本？
更专业一点就是随机抽样，随机现象包括随机变量和随机过程，随机抽样定义：从随机变量或随机过程服从的概率分布中获取其样本值的数学方法。因此需要思考3个问题：
问题1: 概率分布有不同的形式，有些随机变量和过程的概率分布是显示的，有些是未知的(不知道具体的函数)，该如何获取样本。
问题2: 生成的样本需要满足什么样的条件，比如既然生成的样本来自P(x)分布，这些样本是否可以重新计算出新的P’(x)，两者逼近。
问题3: 生成样本的量化指标，比如样本数量，计算时间，抽样的精度。
本篇针对项目中问题，主要讨论概率分布已知的形式(如均匀分布、泊松分布等)的采样，概率分布的未知情况下的采样，这种采样的分布是近似的采样方法，得到的样本也是近似的，因此需要讨论抽样的精度和效率。

2 基本采样思想

本节先从随机变量的抽样引出基本的抽样方法，然后将问题扩展到随机向量和随机过程的采样，抽样方法的技巧是不同的。
2.1 直接抽样
原理就是直接从完全己知概率分布出发, 利用均匀分布U(0 1) 的随机数,使用严格精确的数学方法构造抽样算法，产生随机变量和随机
过程的样本值 x, 使得样本具有独立同分布。主要包括3中方法：
方法1：逆变方法(CDF)，f(x)概率分布已知，可以直接经过初等变换得到累计概率分布的反函数。常见的如下：

   方法2：取舍方法(Rejecting Sampling)，从建议概率分布中直接抽样，建议概率分布解析式可求。

随机变量X的概率分布为f(x),建议概率分布为g(x), 对样本X，f(x) < cg(x) . 取舍方法是从g(x)中抽样产生样本X, 由上图可知部分样本落在法f(x)分布中，以f(x)/cg(x) 接受样本X，以1- f(x)/c*g(x)拒绝样本X。
方法3：复合方法，随机变量X服从的条件概率分布g(x|y), 而随机变量Y的概率分布服从h(y), 复合分布的形式如下：

复合分布的算法流程主要是：
(1) 从概率分布h(y)抽样产生Y
(2) 从条件概率分布g(x|Y)抽样产生X

直接抽样思想中逆变方法没有判别和取舍问题，一次计算就可以得到样本，对于取舍和复合方法需要进行多次搜索比较，抽样过程也使用了随机数，无法确定需要搜索多少次，无规律可循，因此每次计算不要不同的运算量，抽样的效率很低。是否可以找到更高效的抽样方法了？

2.2 典型采样技巧

2.2.1 变换算法

将复杂的概率分布抽样变换为简单的概率分布抽样，同时给出两者的变换关系。即知道随机变量概率分布，对应的随机变量函数也是随机变量。将需要从复杂的概率分布f(x)抽样产生X，变成先从简单的概率分布法f(y),抽样产生样本Y，由变换关系x = g(y), 将样本值Y带入得到最终样本值X。 抽样算法如下：
 (1) f(y)抽样产生Y 
 (2) 样本值X= g(Y)
 将一维变换方法推广到二维，二维联合概率分布f(x,y)抽样比较复杂，可以先从f(v,w)抽样产生V，W, 然后由变换关系x = g1(V,W) 和y = g2(V,W)变换得到最终的样本值X,Y。
 标准正态分布采样使用了变换算法的思想技巧，也就是著名的Box-Muller变换。

2.2.2 随机向量抽样

以上主要是二维随机变量，对于多维随机变量构成的随机向量联合概率密度函数f(x1, x2,…,xn)，分为两种情况：
（1）随机向量相互独立：联合f(x1, x2,…,xn) = f1(x1)*f2(x2)*f3(x3)…*fn(xn), 可独立的对概率分布抽样。
（2）随机向量相关：联合概率分布可以分解为单个随机变量的边缘分布与多个随机变量的条件概率分布乘积。

首先从边缘概率分布 f1(x1)抽样，得到随机变量X1的样本值，然后依次地从各个条件概率分布f(x2|x1,x2…xn)抽样，得到随机变量X2,X3,…Xn的样本值。
随机向量抽样主要问题：是条件概率密度函数是多维的积分计算，无解析形式，理论上可行但是实际很难做，导致抽样失效。

2.2.3 随机过程抽样

	随机过程简记为X(t)，随机过程又称为随机函数，是定义于基本概率空间上的一族随机变量, 这些随机变量是相关的, 不是独立的。对于完全已知的概率分布随机过程抽样包括需要掌握以下两类：

（1） 正态过程抽样
（2） 马尔可夫过程抽样
给出马尔可夫过程的X(t)的联合概率密度函数：
f(x(t0),x(t1),x(t2)…x(tn)) = f(x(t0)) * f(x(t1)|*(t0)） …f(x(t m))*f(x(tm-1))
当t =0，从初始概率分布f(x(t0))产生样本X0； 从条件概率分布f(x(t+1)|x(t))抽样产生样本x(t+1) ,循环抽样。

2.3 马尔科夫链蒙特卡洛技巧

	直接抽样方法在下面两种情况遇到了困难:一是由于抽样效率太低，直接抽样算法失效;二是概率分布是不完全己知概率分布，归一化常数是不知道的，直接抽样方法不适用。本节主要讨论MCMC方法的基本理论。首先介绍一下马尔科夫链的重要性质，MC(马尔科夫链)与MC(蒙特卡罗)有啥关系？

2.3. 1马尔可夫链性质

     马尔可夫链的平稳分布，来源于向量转移概率定义：对于非零向量X ={X1，X2,....,Xn},使得：XP = X，其中P为一概率矩阵，则称X为P的固定概率向量。
	 设X ={X1，X2,....,Xn}为一状态概率向量，P为状态转移概率矩阵，若XP=X，则称X为该马尔科夫链的一个平稳分布。系统的平稳状态指的是系统再经过一步状态转移其状态(平稳状态概率向量)概率保持不变。马尔可夫链的细致平衡条件如下：状态转移矩阵P 和概率分布 π ( x ) 对于所有的 i , j 满足：
	 																		π(i)P(i,j)=π(j)P(j,i) 
	 则称概率分布 π ( x )是状态转移矩阵P的平稳分布。当给定P(习题集会给定P计算平稳分布)，为已知转移矩阵，通过联合方程或者特征值分解，得到平稳分布 π ( x )。实际问题是状态转移矩阵P大多数是不太好找的，即很难随机找到一个状态转移矩阵Q满足细致平衡条件：
	 																       π(i)Q(i,j) ≠ (j)Q(j,i)

这个时候要靠MCMC采样来解决这个问题，引入一个参数α(i,j)，使等式成立：
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
那么具体的α(i,j)需要满足：
α(i,j)=π(j)Q(j,i)
α(j,i)=π(i)Q(i,j)
因此通过构造一个随机转移矩阵Q，和一个参数α(i,j)，得到了π ( x )对应的状态转移矩阵P。这里给Q和α(i,j)进行了命名：
(1) 转移矩阵Q：建议转移分布，具体是条件概率分布Q(x=j|x=i)，在状态x=i的条件下状态x=j的条件概率
（2）参数α(i,j) ： 接受概率， 0<α(i,j)<1
根据前面的拒绝-接受思想，在原来的马氏链上，从状态i以建议转移矩阵Q(i,j)的概率转移到状态j的时候，以α(i,j)的概率接受。 即目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q 一定的接受率获得。

因此，可以总结马尔可夫链采样方式:
(1) 首先选择一个建议概率分布Q，从建议概率分布抽样产生一个候选样本值
(2) 然后建立可操作的状态转移规则，己知概率分布在状态转移规则下是不变的
(3) 根据状态转移规则，用接受概率α(i,j)判断状态是否转移，重复这个过程，将产生一条马尔可夫链
(4) 如果所生成的马尔可夫链是不可约的、非周期的和遍历的，经过足够多步的状态转移后，马尔可夫链的平稳分布将渐近已知概率分布
(5) 马尔可夫链的各个状态就是随机变量的样本值，有了样本值可以采用蒙特卡罗模拟计算复杂积分。

将以上5条转化为数学符号描述，得到MCMC过程如下：

样本集(xn1,xn1+1,…,xn1+n2−1)即为我们需要的平稳分布对应的样本集。由于α(xt,x∗)可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。因此引入M-H算法。

2.3.2 M-H采样

	把建议概率分布从对称形式推广到非对称形式，从而把梅特罗波利斯算法(M)推广到一般情况，这种算法称为黑斯廷斯算法（H），也就是通用梅特罗波利斯算法(MH)。

黑斯廷斯算法只是把梅特罗波利斯算法的建议概率分布从对称形式推广到非对称形式，并没有改进梅特罗波利斯算法。

2.3.3 改进方法

马尔可夫链蒙特卡罗抽样方法的其他算法都是通用梅特罗波利斯算法的改进和发展，针对上述问题，提出如下两类方法做为改进方向：
（1）改进建议概率分布Q：无论是梅特罗波利斯算法，还是黑斯廷斯算法，对于如何选择建议概率分布都很笼统，梅特罗波利斯算法的建议
概率分布是对称形式，黑斯廷斯算法把建议概率分布推广到非对称形式。建议概率分布可任意选取，任意性很大，建议概率分布的选择并不是
唯一，不同的建议概率分布，收敛速度差别很大。因此寻找一个理想的建议概率分布犹如一门艺术，要找到一个好的建议概率分布相当困难。典
型的早期方法如建议概率分布采用具体的独立游走分布。

目前比较重要的方法是建议概率分布是独立建议概率分布：
1. 单个独立抽样算法
2. 多个独立抽样算法
3. SIS独立抽样算法 ：粒子滤波应用

    （2）全条件概率分布抽样方法：全条件概率分布抽样方法使得抽样遵循己知概率分布的动态性，其显著特点是每步法代都用全条件概率分布来
    构建马尔可夫链，这些全条件概率分布是通过将己知概率分布限制在一定子空间产生的。共同点是马尔可夫链转移规则建立在全条件概率分
    布抽样的基础上，目的是为了提高马尔可夫链蒙特卡罗方法的收敛速度，条件概率分布抽样方法最大特点是接受概率为1，不用进行状态转移
    判断，抽样效率高。

最典型的抽样是Gibbs采样(吉布斯抽样),以下重点介绍上面改进的方法。

2.3.3.1 独立抽样

独立抽样算法的特点是建议概率分布选择为独立建议概率分布，独立建议概率分布不是条件概率分布，只与状态j有关，与状态i 无关，独立建议概率分布Q(x=j,x =i) = Q(x=j), 独立建议概率分布只有一个转移状态x =j, 对应的算法称为单个独立抽样算法。在当前状态x=i,从Q(x=j)抽样得到样本值 j ，接受概率变为：

其中，w称为权重，w(j) = 概率分布 π / 建议概率分布Q , 独立抽样算法的状态转移规则是：从独立建议概率分布Q(j)抽样产生候选样本值X*，若u<α(i,j)，则选取候选样本值X*，否则选取当前样本值X。
独立抽样算法的抽样效率依赖于独立建议概率分布Q(j)与己知概率分布π(j)的接近程度。

多个独立抽样马尔可夫链步长越小，马尔可夫链移动速度就越慢，但是如果马尔可夫链步长变大，则将导致接受概率很低，拒绝率很高。这两种情况都使得梅特罗波利斯算法收敛速度慢，算法效率低。为了加速收敛，提高抽样效率，Frenkel和Smit(1996)提出多个独立抽样算法有多个转移状态使得抽样在不降低接受概率的情况下，进行大步长跳跃转移，加快收敛速度。这里不做深入总结。

SIS独立抽样序贯抽样是一个迭代过程， 在序贯抽样基础上加上重要抽样， 构成序贯重要抽样。
我们知道上面讨论的目标平稳分布为π(X), X={x1,x2,x3,…,xm },在抽样的过程中目标概率分布的样本容量m是固定不变的,叫做固定样本抽样。而序贯抽样，不用规定样本容量，先抽取少量样本，根据抽样结果，持续的抽样，这种样本数量变化的叫做序贯样本抽样。序贯抽样的目标平稳分布可以表示为：

根据独立抽样算法：独立建议概率分布Q(x=j,x =i) = Q(x=j)，即Q(y,x) =Q(y)或Q(y|x) =Q(y),更换变量名。则序贯抽样的独立建议概率分布表示为：

序贯重要抽样权重的递推公式为：

SIS独立抽样算法的接受概率为：

SIS 独立抽样算法的状态转移判断准则为:从Q(Y) 抽样产生候选样本值Y ， 若u< α(X,Y)， 则选取候选样本值Y ， 否则选取当前样本值X 。随机变量的概率分布是变化的， 每次法代有不同概率分布每次迭代都要进行多次蒙特卡罗模拟， 所以很花费时间， 计算量很大， 但是随着迭代次数增加，有可能得到的估计值方差较小， 从而减小误差提高估计精度和计算效率。此处可以思考两个问题：
(1) 递推权重是否会带来样本退化，怎么解决？
(2) 蒙特卡罗模拟的估计值方差如何减小？

2.3.3.2 Gibbs抽样

Gibbs抽样属于全条件概率分布抽样方法，与以上抽样方法有哪些改进了，先看目标平稳概率分布：

单随机变量x的全条件概率分布：

吉布斯算法首先将己知概率分布分解为多个单变量全条件概率分布然后通过对多个单变量全条件概率分布抽样来实现对己知概率分布抽样， 称为多级吉布斯抽样， 若只有两个单变量，则称为两级吉布斯抽样。和随机向量抽样有点相似，有啥区别了：
(1) 随机向量直接抽样算法需要计算边缘概率分布和条件概率分布的高维积分， 其积分计算和抽样方法都是很困难的。
(2) 吉布斯算法只需要计算全条件概率分布的一维积分， 其计算和抽样则比较容易。

再看建议概率分布 从状态X=i到状态X=j建议概率分布Q(i|j), 与X=j到状态X=i的建议概率分布Q(j|i) 满足：

代入接受概率为：

Gibbs的接受概率为1，马尔可夫链转移规则不需要进行状态接受判断，没有状态拒绝这一步，所以收敛速度比较快，抽样效率较高。
具体采用过程可参照相关博客。

3 参考文献

[1]: https://www.cnblogs.com/pinard/p/6638955.html
[2]: https://blog.csdn.net/anshuai_aw1/article/details/84875521
[3]: 《蒙特卡罗方法理论和应用》康崇禄