机器学习中的采样方法
首先，我们需要了解为什么需要采样，也就是说采样的目的是什么。一般来说，对于模型f(z)，z服从概率分布p(z)，我们需要求该模型的期望： \(\mathbb{E}[f]=\int f(z)p(z)dz\)
然而，$\mathbb{E}[f]$的积分往往难以求得，但如果我们可以从p(z)中独立采样出$L$个样本$z^{(l)}, l=1,2,…,L$， 那么由大数定律我们可以近似的认为： \(\hat{f}=\frac{1}{L}\sum_{l=1}^Lf(z^{(l)})\)
这就是我们常说的蒙特卡洛积分，基于该原理（用多个样本的均值来近似积分）的采样都称为蒙特卡洛采样。 那么问题就来了，我们如何获取从p(z)中独立采样出的$L$个样本$z^{(l)}$呢？这就是sample method要解决的问题了。
how can we obtain independent samples from a distribution p(z) we do not know how to sample from?
目录： （一）基本的概率分布采样方法 （二）马尔可夫链蒙特卡洛 （三）吉布斯采样
（一）基本的概率分布采样方法
1. 转换采样 transformation method
- 假设我们可以轻松的获得服从均匀分布（0，1）的样本x；
-
我们可以用概率密度函数（pdf）p(z)求得z的累计分布函数（cdf）h(z)，h(z)的取值范围必定为[0,1]的均匀分布，如下图（把y换成z = =）：
- 那么这时我们就可以认为x是h(z)的独立采样样本，即x=h(z), 则可以求得h(z)的反函数$z=h^{-1}(x)$
- 即 从均匀分布的[0,1]中采样x，通过反函数计算得出z，则z是 从分布p(z)中独立采样出的样本。
2. 向前采样 （ancestral/forward sampling）
对于一个取值为k的离散变量v，其取值概率分别为$\theta_1,…,\theta_k$。我们要如何利用均匀分布来对k进行采样呢？
很简单，我们可以将均匀分布按照k的概率分布进行如下切片：
这样我们就可以从[0,1]均匀分布中随机采出一个值x来，判断x的值属于哪个区间，进而将x转换为变量v的采样样本。
由此，对于包含n个变量的样本，我们很容易的想到使用贝叶斯网络进行采样，这种方法称为ancestral sampling (或forward sampling)。对于一个已经构建好的贝叶斯网络（如下图），我们可以从根节点变量开始采样，之后采用各个节点下的条件概率分布（CPD）表进行采样，重复这个过程直到n个变量都被采到。
在我们的图例中，若我们要对样本$(\text{student’s grade, exam difficulty, intelligence level}) \rightarrow (g, d, i)$进行采样。首先，我们先分别采样出$d$和$i$，然后根据CPD$p(g|d,i)$ 采样出$g$得到我们的样本$(g, d, i)$。
3. 拒绝采样 （Rejection sampling）
当我们难以直接从$p(z)$ 进行采样，我们可以假设$p(z)$ 可以通过任意个样本$z$来估计，即 \(p(z)=\frac{1}{Z_p}\tilde{p}(z)\)
并且我们假设已知一个分布$q(z)$ 存在一个常量k满足 \(kq(z)\geq\tilde{p}(z)\)
则我们可以使用拒绝采样来获取服从$p(z)$分布的独立样本。
- 从$q(z)$分布中采样$z_0$点
- 从均匀分布$[0, kq(z_0)]$中随机采样$u_0$点
- 若$u_0>\tilde{p}(z_0)$, 则拒绝样本$z_0$，否则接受样本$z_0$
- 重复1～3步骤，所有接受样本$z$的集合可以认为是从分布$p(z)$独立采样出的样本集合。
Note: 为了提高采样效率（降低拒绝样本的出现概率），我们应该尽可能地选取靠近$p(z)$的分布$q(z)$
然而，一个样本被接受的概率是多少呢？ \(p(accept) = \int \frac{\tilde{p}(z)}{kq(z)}q(z)dz=\frac{1}{k}\int\tilde{p}(z)dz\) 即在符合条件$kq(z)\geq \tilde{p}(z)$的条件下， k越小，样本被接受的概率越大。
然而对于高维空间中的样本，其接受率随着空间维度指数下降（维度诅咒），参数k的确定也变成了一种挑战。
4. 重要性采样 (Importance sampling)
拒绝采样的缺点来源于拒绝样本这个动作会带来采样过程中额外的开销。于是重要性采样就弥补了这一缺点：对于通过$q(z)$得到的样本，全部接受。然而，通过$q(z)$采出的样本显然不符合原本的$p(z)$分布，为了矫正这个偏差，我们给每个样本附一个重要性权重，比如使$\frac{p(z_0)}{q(z_0)}=1$的样本权重较大，使$\frac{p(z_0)}{q(z_0)}=0.1$的样本权重较小。很容易的，我们就能想到，其实$\frac{p(z_0)}{q(z_0)}$本身就是一个很好的权重标识：
\begin{align}
\mathbb{E}[f]=& \int f(z)p(z)dz
=& \int f(z)\frac{p(z)}{q(z)}q(z)dz
\approx& \frac{1}{L}\sum_{l=1}^{L}f(z^{(l)})\frac{p(z^{(l)})}{q(z^{(l)})}
=& \sum_{l=1}^L w_lf(z^{(l)})
\end{align}
5. 重要性重采样 (Sampling-importance-resampling)
我们可以看到，重要性采样是用来近似函数f(z)的期望，而非样本分布本身。我们可以通过重要性的重采样来获得服从p(z)分布的独立样本。
- 从分布q(z)中采样L个样本$z^{(1)}, z^{(2)},…,z^{(L)}$
- 计算每个样本点$z^{(l)}$的权重$\frac{p(z^{(l)})}{q(z^{(l)})}$
- 将计算出的权重归一化（normalize）为$w_1, w_2,…,w_L$, 其中$w_l=\frac{p(z^{(l)})}{q(z^{(l)})}\frac{1}{\sum_m\frac{p(z^{(m)})}{q(z^{(m)})}}$
- 根据归一化后的权重$w_1, w_2,…,w_L$从样本集$z^{(1)}, z^{(2)},…,z^{(L)}$中随机采样L个样本$\hat{z}^{(1)}, \hat{z}^{(2)},…,\hat{z}^{(L)}$
- 当$L\rightarrow \infty$时，样本$\hat{z}^{(1)}, \hat{z}^{(2)},…,\hat{z}^{(L)}$为服从p(z)分布的样本集。
（二）马尔可夫链蒙特卡洛
马尔科夫链蒙特卡洛 （Marko chain Monte Carlo）同样也是采用建议分布（proposal distribution），即q(z)，来进行采样。与前面介绍的方法不同的是，MCMC的采样过程与采样序列相关，即其proposal distribution 与当前采样样本相关：$q(z|z^{T})$。
1. 马尔可夫链
很容易的我们想到了马尔可夫链：
\(P(X_{t+1}=x | X_t, X_{t-1}, ... X_0)= P(X_{t+1}=x|X_t)\)
马尔可夫链有很好的收敛性质：
马氏链的收敛定理非常重要，MCMC都是以这个定理作为基础理论的。
NOTE:
- 该定理中马氏链的状态可以是无穷多个的
- 定理中的“非周期”可以不过多理解（实际上是我没找到更详细资料），因为我们遇到的问题绝大多数是非周期的马尔可夫链 连通的含义是值i可以通过有限n步转移到j；马氏链的任何两个状态是连通的含义是指存在一个n使得矩阵$P^n$中的任何一个元素值都大于0。
2. Metropolis-Hastings
如果我们能够构造一个转移矩阵为$P$的马氏链，使得马氏链的平稳分布恰好是$p(z)$，那么我们从任何一个初始状态$x_0$出发经过N次P转移，得到一个转移序列$z_0, z_1, z_2,…,z_n,z_{n+1}$，假设马氏链在第n步已经收敛，则$z_n,z_{n+1},…$为$\pi(x)$的样本。
这个idea在1953年被Metropolis首次提出，因此也称为Metropolis算法。虽然Metropolis算法是为了考虑物理中的玻尔兹曼分布的采样问题而提出，但Metripolis算法是首个普适但采样方法，并随之引发了一系列MCMC变种。Metropolis-Hastings就是最常见的一种MCMC改进算法。
转移矩阵P决定了马氏链的收敛性质，因此我们的问题核心就是如何构造P使得平稳分布恰好是p(z)分布，这就要使用到“细致平稳条件”。
细致平稳条件：如果非周期马氏链的转移矩阵P和分布$\pi(z)$满足\(\pi(i)P_{ij}=\pi(j)P_{ji} \hspace{3ex}\text{ for all i, j}\)则$\pi(z)$是马氏链的平稳分布，上式被称为细致平稳条件（detailed balance condition）
通常我们的目标$p(z)$不满足细致平稳条件， 也就是说对于转移矩阵P：\(p(z_i)P_{z_i,z_j}\neq p(z_j)P_{z_j,z_i}\)，因此$p(z)$不是马氏链的平稳分布。但我们可以引入一个因子$\alpha(z_i,z_j)$，使其满足：
\(p(z_i)P_{z_i,z_j}\alpha(z_i,z_j) = p(z_j)P_{z_j,z_i}\alpha(z_j,z_i)\) 例如，我们可以让$\alpha(z_i,z_j) = p(z_j)P_{z_j,z_i}$, $\alpha(z_j,z_i)=p(z_i)P_{z_i,z_j}$，这样就使得上式成立。
于是，我们有\(p(z_i)\underbrace{P_{z_i,z_j}\alpha(z_i,z_j)}_{Q(z_i,z_j)}=p(z_j)\underbrace{P_{z_j,z_i}\alpha(z_j,z_i)}_{Q(z_j,z_i)}\) 即，马氏链Q的平稳分布为p(z)。
实际上我们将引入的因子$\alpha{z_i,z_j}$称为接受率，物理意义上为：对于马氏链P，从状态$z_i$以$P(z_i,z_j)$的概率跳转到状态$z_j$时，以$\alpha{z_i,z_j}$的概率接受这个转移。由此，MCMC采样算法如下：
- 初始化马氏链初始状态$Z_0$
- 对$t=0,1,2,…,$循环以下采样过程
-
第t个时刻马氏链状态为$Z_t=z_t$，采样 $\hat{z} \sim P(z z_t)$ - 从均匀分布[0,1]采样u
-
若$u<\alpha(z_t,\hat{z})=p(\hat{z})P(z_t \hat{z})$ 则接受转移 $z_t\rightarrow\hat{z}$， 即$Z_{t+1}=\hat{z}$ - 否则不接受转移，即$Z_{t+1}=z_t$ 在上述过程中，接受率可能存在偏小的情况，造成马氏链状态拒绝跳转，停留在原地，收敛到平稳分布$p(z)$的速度太慢。我们可以在平稳细致条件的等式两边乘上一个扩大系数M，即： \(p(z_i)P_{z_i,z_j}\alpha(z_i,z_j)\times M= p(z_j)P_{z_j,z_i}\alpha(z_j,z_i)\times M\)
-
通常我们取\(\alpha(z_i,z_j)=\min\{\frac{p(z_j)P(z_j,z_i)}{p(i)P(z_i,z_j)}, 1\}\) 使得两边的接受率同比例放大到其中一边为1。这就是Metropolis-Hastings算法。
（三）吉布斯采样
Gibbs采样是基于Metropolis-Hastings的一种采样方法，适用于高维空间中的样本采样。先上Gibbs采样过程（懒得翻译了，直接截图了）：
证明过程在文献【3】中的“3.4 Gibbs sampling”中写的灰常清楚了，我就不摘抄了。
参考文献
【1】Pattern Recognition and Machine Learning Chapter 11: Sampling Methods, Elise Arnaud and Jakob Verbeek 【2】Sampling methods 【3】LDA-math-MCMC 和 Gibbs Sampling
