常见采样方法

sampling学习小结

本文是在阅读prml的第11章后，写的一个总结。

---

  首先要说说为什么需要采样。通常没我们会遇到很多问题无法用分析的方法求得精确解，例如待求问题很复杂。遇到这种情况的时候，人们通常会采用一些方法得到近似解，其中随机模拟就是一类近似求解的方法，采样就是随机模拟的核心。

-基本方法：
从基本概率分布产生新变量分布,假设有某分布的分布函数是$F(x)$，若$y$是$(0,1)$上的均匀分布随机变量，那么随机变量$F^{-1}(y)$的分布函数是$F(x)$。
问题：$F^{-1}(y)$难解

-基于proposal distribution方法：
从分布$p(z)$采样很困难，但是可以找到一个相对更容易的分布$q(z)$，这个$q(z)$称之为proposal distribution。

• Rejection sampling
  

具体操作如下，设定一个方便抽样的函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。（参考上图）
- x 轴方向：从 q(x) 分布抽样得到 a。
- y 轴方向：从均匀分布（0, kq(a)) 中抽样得到 u。
- 如果刚好落到灰色区域： u > p(a), 拒绝， 否则接受这次抽样
- 重复以上过程

为什么需要乘上k呢？
对于两个分布，其积分面积是一致的，因为都等于概率1， 不能时刻保证任意一个分布大于另一个分布，所以乘上k后能够保证这一点。

在选择k的过程中，需要选择尽量小，这样能保证拒绝率小，从而提高采样的的效率。

• Adaptive Rejective Sampling

  拒绝采样的弱点在于当被拒绝的点很多时，采样的效率会非常不理想。同时我如果能够找到一个跟目标分布函数非常接近的参考函数，那么就可以保证被接受的点占大多数（被拒绝的点很少）。这样一来便克服了拒绝采样效率不高的弱点。如果函数是 log-concave 的话，那么我们就可以采样自适应的拒绝采样方法。

什么是log-concave呢？
就是对任意的p(x) 如果他的log p(x)是凹函数，那么p(x) 就是log-concave.

在对数图像上找一些点做图像的切线，如下图所示。因为对数图像是凹函数，所以每个切线都相当于一个超平面，而且对数图像只会位于超平面的一侧。

因为我们之前对p(x) 取了对数，还原回去就是做一个指数运算。如下图中的一些曲线，这些曲线会被原函数的图形紧紧包裹住。特别是当这些的指数函数变得很多很稠密时，以彼此的交点作为分界线，我们其实相当于得到了一个分段函数。这个分段函数是原函数的一个逼近。用这个分段函数来作为参考函数再执行Reject Sampling，自然就完美的解决了我们之前的问题。

- Importance sampling
假设对$p(z)$的采样困难，但是对一个给定的$z$，却可以容易计算其值。若现在不是要计算分布的采样，而是计算一个函数$f(z)$在这个分布下的期望$E(f)=\int p(z)f(z)dz$。

$E(f) =\int p(z)f(z)dz = \frac{p(z)f(z)}{q(z)}q(z)dz \approx \frac{1}{L}\sum_{l=1}^L \frac{p(z^l)}{q(z^l)}f(z^l)$
其中，$\frac{p(z^l)}{q(z^l)}$看作是important weight。
问题：同样在于在高维空间中寻找到合适的$q$很难。

• Metropolis-Hastings

在说到Metropolis-Hastings算法之前，需要提一个很重要的平稳分布：$p^*(z)=\sum_{z'}T(z,z')p^*(z^{'})$，表示转以后分布不再发生变化，这样一来，我们可以控制马尔科夫收敛到我们指定的分布$p^*(z)$。这也就是说，通过稳态的Markov Chain进行转移计算，等效于从$p^*(z)$中采样。再者，马氏链的收敛性质是由转移矩阵决定的，所以基于马氏链的采样问题就转化为如何构造矩阵，使得平稳分布恰好是我们要的分布$p^*(z)$。

为了得到平稳分布，需要满足细致平稳条件，细致平稳充分不必要条件公式是：$p^*(z)T(z,z') =p^*(z')T(z',z)$，由此，我们可以通过设置接受率的方式来满足条件，接受率$A(z^*,z^{(\tau)}) = min(1, \frac{\widetilde{p}(z^*)q(z^{\tau}|z^*)}{\widetilde{p}(z^{\tau})q(z^{*}|z^{\tau})}$)。

最后，给出M-H算法的描述：
有一个易采样的分布$q(z)$，假设当前已经采样出的样本是$z^{(\tau)}$，那么下一个样本从分布$q(z|z^{(\tau)})$中取得：记新采样的样本为$z^*$，以概率$A(z^*,z^{\tau})$接受该样本。如下图所示：

综上，Metropolis方法需要先选一个比较容易取样的proposal distribution，从这个分布里取样，然后通过接受率决定是否采用这个样本。一个简单的例子就是对于proposal distribution我们可以采用Gaussian centred on the current state，其实很好理解，就是上一步节点的值可以做下一步节点需要采样的proposal distribution即高斯分布的均值，这样下一步节点的状态由上一步完全决定，这就是一个马尔科夫链。但是，高斯分布的方差是固定的，其实方差就是步长，如何选择步长这是一个state of the art问题，步子太小扩散太慢，步子太大，拒绝率会很高，原地踏步。

• gibbs sampling
  gibbs采样其实是每次只对一个维度的变量进行采样，固定住其他维度的变量，然后迭代，可以看作是接受率为1的Metropolis-Hastings的特例。

  

考察$x$坐标相同的两个点$A(x_1,y_1),B(x_1,y_2)$，我们可以的得到：
$p(x_1,y_1)p(y_2|x_1) = p(x_1)p(y_1|x_1)p(y_2|x_1)$ $p(x_1,y_2)p(y_1|x_1) = p(x_1)p(y_2|x_1)p(y_1|x_1)$
因此，$p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1)$

由此等式，我们发现，在$x=x_1$这条平行于$y$轴的直线上，如果使用条件分布$p(y|x_1)$作为任意两点的转移概率，那么任意两点之间的转移满足细致平稳条件。于是可以构造如下的概率转移矩阵Q:

为了更好说明gibbs采样时M-H的特例，我们从M-H的角度重新认识下Gibbs采样。考虑对$z_k$的更新采样，在得到$z_k^{\tau +1}$后的$z$记为$z^*$，那么有$q(z^*|z) = p(z_k|z{''})$,其中$z''$表示集合中去掉$z_k$。在$z$和$z^*$之间，只有$z_k$是不同的，因此$z''=z^{*''}$

在M-H算法中的接受率始终为1：

将该算法推广到高维，就可以得到高维的Gibbs采样：

需要指出的是，算法在收敛后，得到的分布$p(x_1, x_2,...,x_n)$的样本，这些样本不独立，但是在所要求的采样目标就是采样得到符合给定概率分布的样本，并不要求各样本独立，因此，仍符合要求。

---

参考
【1】https://wenku.baidu.com/view/cc4c8eaca58da0116c1749dc.html
【2】http://www.cnblogs.com/daniel-D/p/3388724.html
【3】http://blog.youkuaiyun.com/losteng/article/details/51098647
【4】http://www.52nlp.cn/