采样方法【1】

对于大多数的概率模型，直接推导求取它的参数，如均值，积分等，通常是很棘手的，比如我们要求函数$f(z)$在概率分布$p(z)$下的期望：

$$E(f) = \int f(z)p(z)dz$$
但是如果我们可以从模型的概率分布采样到足够多的数据 $z^{(l)},l=1,...,L$，根据大数定理，期望可以用样本的均值来逼近
$$E(\hat f)=1/L\sum_{l=1}^Lf(z^(l))$$
而且
$$var(\hat f) = \frac{1}{L}E[(f-E(f)]^2$$

1. 原始采样法(Ancestral Sampling)

有向图模型的概率分布可以表示成条件概率的积：

$$p(z) = \prod_{i=1}^Mp(z_i|pa_i)$$
$z_i$是与节点i相关联的变量， $pa_i$指节点i的父节点相关联的变量。
根据这些变量在有向图中的拓扑顺序依次采样；如果有些变量已经有观察到的值，那么在采样过程中需要把采样的值与已知值做比较，相等则保留采样的值，不相等则这一轮采样的值全部丢弃，从头再来。这种方法的性能随着已观测变量的数目增加急剧降低，所以这种方法在实际中很少会采用。

2. 基本采样法（Basic Sampling）

假设$z$服从某种简单的标准分布，如均匀分布，且$z=f(y)$，那么

$$p(y) = p(z)\left|\frac{dz}{dy}\right|$$
对上式求积分，即
$$z = h(y) = \int_{-\inf}^{y}p(\hat y)d\hat y$$

直接举个简单的例子吧，比如指数函数

$$p(y) =\lambda exp(-\lambda y)$$
对上面的概率密度函数积分，得到
$$z = h(y) = 1 - exp(-\lambda y)$$
$z$的取值范围是0到1，这样子我们可以从0到1的均匀分布采样得到z的值【很多语言都能产生从0到1的伪随机数】，并通过$y=h^{-1}(z)$求出y

Box-Muller method

这个方法可以从高斯函数中采样，见PRML p527

3. 拒绝采样（Rejection Sampling）

在实际应用中，$p(z)$并不服从一些简单分布，无法直接从$p(z)$采样的时候。通常，

$$p(z)=\frac{1}{Z_p}\tilde p(z)$$
$\tilde p(z)$可以很容易求得，但是归一化项 $Z_p$的值是不知道的
这时候我们需要一个相对简单的，能直接采样的分布 $q(z)$，叫做建议分布（proposal distribution），使得存在 $k$，对于所有的$z$右边的不等式成立： $kq(z) \ge \tilde p(z)$
对z进行采样的方法如下：
1) 从 $q(z)$采样得到 $z_0$
2) 从 $[0, kq(z)]$的均匀分布中采样得到 $u_0$
3) 如果 $u_0>\tilde p(z_0)$，拒绝采样，否则接受 $z_0$作为采样值
所以一个采样值会被接受的概率为
$$\begin{align} p(accept) & = \int{\tilde p(z)/kq(z)}q(z)dz \\ & = \frac{1}{k}\int \tilde p(z)dz \end{align}$$
也即下图中非阴影部分表示的区域

所以为了使接受的概率比较大，k应该尽可能小

4. 自适应拒绝采样（Adaptive Rejection Sampling

在实际中，合适的$q(z)$也是很难确定的
当$p(z)$为log凹函数的时候，可以采用ARS：在$lnp(z)$上的某些均匀间隔的点$z_i$计算函数值和梯度可以得到相应的切线方程，切线方程的交可以形成一个包络函数(envelope function)，如下图

所以 $q(z)$的方程由下式给出：
$$q(z) = k_i\lambda_iexp(-\lambda_i(z-z_i)) \qquad \hat z_{i-1,i} \le z \le \hat z_{i,i+1}$$
$\hat z_{i-1,i}$是点 $z_{i-1}$和 $z_i$的切线方程的交点， $\lambda_i$是斜率， $k_i$是截距
可以应用和拒绝采样一样的算法，不同的是当一个采样被拒绝时，我们可以再细分化网格点，用更多的切线方程来逼近原函数分布。
拒绝采样在低维【1，2维】的时候是一个有用的方法，但是并不适用于高维的情况。

5. 重要性采样（Importance Sampling）

重要性采样其实是一种计算期望的方法，算法本身并不能从$p(z)$产生样本。本文一开头提到的用大数定理逼近求期望的方法在高维的情况下并不适用【网格化，维数灾难】；而且概率密度函数的“质量”经常只局限在相对小的某个区域内，在高维情况下，很小部分在这个区域内采样的样本都将会对结果产生巨大的影响。重要性采样即是要在$p(z)$，或者更准确得说是在$p(z)f(z)$比较大的地方采集更多的样本，而在比较小的地方采集相对少的样本，以此减少工作量【所谓“重要性”】

【p(z)存在归一化项不可求的情况也是类似的，这里将不再赘述，可以参考PRML第11章的内容】
同拒绝采样，我们需要一个容易采样的分布$q(z)$：

$$\begin{align} E(f) & = \int f(z)p(z)dz \\ & = \int f(z)\frac{p(z)}{q(z)}q(z)dz \\ & \simeq \frac{1}{L} \sum_{l=1}^L\frac{p(z^{(l)}}{q(z^{(l)}}f(z^(l)) \end{align}$$
似然比 $w(z^(l)) = p(z^{(l)}/q(z^{(l)}$就是重要性权重。
和拒绝采样一样，这种方法依赖于 $q(z)$和 $p(z)$的相似程度。
如果有多个 $q(z)$可供选择，那么要选那个最小化估计期望方差的那个：
$$min\quad var_{q(z)}(f(z)w(z)) = E_q(z)(f^2(z)w^2(z)) - E^2(f)$$
第二项并不依赖于q(x)所以只需要最小化第一项，根据琴生不等式可以得到它的下限：
$$\begin{align} E_{q(z)}(f^2(x)w^2(x)) & \ge (E_{q(z)}(|f(z)|pz)))^2 \\ & = \left(\int |f(z)|w(z)q(z)dz\right)^2 \\ & = \left(\int |f(z)|p(z)dz\right)^2 \\ \end{align}$$
所以最优重要性分布（optimal importance distribution）为：
$$q^*(x) = \frac{|f(z)|w(z)}{\int |f(z)|w(z)dz}$$
虽然我们无法从 $|f(z)|w(z)$中直接采样，但是上面的结果表明，在 $|f(z)|w(z)$比较大的区域中采样比较高效；重要性采样也可以当成是蒙特卡罗算法中降低方差的一种方法【蒙特卡罗算法中用的是 $q(z)=p(z)$】，这个性质使得可以它可以在实际中用来估计尾概率【tail probability】，即用 $q(z) = I(z)p(z)【I是指示性函数，如果z属于重要区域则为1，否则为0】$，如下图：

重要性重采样

如果想要从重要性采样估计的分布中采样，可以这样子做：首先从$q(z)$中采样$z^{(l)},l=1,...,L$，然后计算$w^{(l)}$，最后从$z^{(l)}$中根据概率$w^{(l)}$采样产生新的$L$个样本。当$L$趋向于无穷大时，$L$个样本趋向于遵从p(z)的分布。
证明：考虑重采样的累积分布函数【$I$是指示性函数】

$$\begin{align} p(z \le a) & = \sum_{l:z^{(l)} \le a} w^{(l)} \\ & = \frac{\sum_lI(z^{(l)} \le a)p(z^{(l)})/q(z^{(l)})}{\sum_lp(z^{(l)})/q(z^{(l)})} \end{align}$$
当 $L \to \infty$时，求和可以转化为相对于 $q(z)$的积分：
$$\begin{align} p(z \le a) & = \frac{\int I(z \le a)\{p(z)/q(z)\}q(z)dz}{\int \{p(z)/q(z)\}q(z)dz} \\ & = \frac{\int I(z \le a)p(z)dz}{\int p(z)dz} \\ & = \int I(z \le a)p(z)dz \end{align}$$
证明结束。

6. 采样与EM算法

【哎，找个时间把EM算法也写一遍再来填坑吧】

参考资料：
[1] Pattern Recognition and Machine Learning
[2] An Introduction to MCMC for Machine Learning
[3] https://en.wikipedia.org/wiki/Importance_sampling