Book Notes on Pattern Recognition and Machine Learning (PRML)
我们首先来考虑一个最基本的问题,在给定某范围内均匀分布的随机变量后,我们如何用其得到任意分布的随机变量。
假设我们已知均匀分布在 (0,1) 上的随机变量 z,我们对 z 进行变换得到 y,即 \(y=f(z)\) ,那么 y 的分布可以写成
\[p(y) = p(z)|\frac{dz}{dy}|\]在这个例子中,z 是均匀分布的,所以 \(p(z) = 1\) 。我们的目标是找到这样一个变换函数 \(f(\cdot)\) 使得 y 的分布为 \(p(y)\) 。通过积分可以将这两者联系起来
\[z=h(y) \equiv \int_{-\infty}^y p(\hat{y})d\hat{y}\]所以对应的变换函数就是上式的反函数,即 \(y = h^{-1}(z)\)
这里通过一个例子来说明,假设我们期望的 y 服从指数分布,即 \(p(y) = \lambda exp(-\lambda y)\) ,这里 y 为非负数,那么此时我们可以得到 z 和 y 的关系
\[z=h(y) \equiv \int_{-\infty}^y \lambda exp(-\lambda \hat{y})d\hat{y} = 1-exp(-\lambda y)\]我们对其求反函数,即可得到变换函数
\[y = h^{-1}(z) = -\frac{\ln (1-z)}{\lambda}\]所以,通过上式我们实现了从均匀分布中得到指数分布
拒绝采样思想最典型的应用是求圆周率,我们通过随机生成点计算到原点的距离,接受落在圆内的点,以此来进行估计。
假设我们要从一个和普通标准分布不同的不常见分布 \(p(z)\) 中采样,直接从该分布中取样不太现实。但我们通常可以对任意给定的 z 值计算 \(p(z)\)
给定一个概率分布 \(p(z) = \tilde{p}(z)/Z_p\) ,其中 \(\tilde{p}(z)\) 已知,\(Z_p\) 为未知的归一化常数。我们想对 \(p(z)\) 进行拒绝采样,需要引入一个参考分布 \(q(z)\) ,这是一个人为选择的容易采样的分布。接下来我们找一个常数值 k 对任意值都满足 \(kq(z)\geq \tilde{p}(z)\)
现在我们就可以借助 \(q(z)\) 来进行拒绝采样。我们首先从 \(q(z)\) 中随机生成一个数 \(z_0\) ,接着从区间 \([0,k(z_0)]\) 中均匀随机生成另一个数 \(\mu_0\) 。如果 \(\mu_0 > \tilde{p}(z_0)\) 则拒绝该采样;否则接受。最后剩下的接受样本就是我们的采样结果。
为了保证采样被接受的概率,我们需要选择满足条件的尽可能小的 k 值,否则将会有非常多的样本被拒绝浪费。但现在仍然存在一个问题,那就是参考分布 \(q(z)\) 的选取。我们希望参考分布和目标分布之间尽可能的相似,但通常我们很难根据先验知识直接构建一个合适的 \(q(z)\) ,所以这里考虑利用已知的 \(p(z)\) 测定值构建包络函数
当 \(p(z)\) 是对数凹分布时,我们可以选择 \(\ln p(z)\) 的若干处的切线及交点组成包络函数。其本身由一系列分段指数分布组成,即
\[q(z) = k_i\lambda_iexp\{-\lambda_i(z-z_i)\}, \quad z_{i-1} < z \leq z_i\]此时我们在拒绝采样的过程中,也可以根据被拒绝的样本逐步更新优化我们的包络函数 \(q(z)\) ,这种方法也被称为 Adaptive rejection sampling
但拒绝采样只适用于一维二维的低维数据。当维度变高时,k 的值会逐渐增大,这意味着接受率会随着维度增加而指数级减小。同时参考分布也会随着维度增加变得更难确定。
我们采样的根本目的还是为了计算期望,即
\[E[f] = \int f(z)p(z)dz\]而重要采样可以直接用来近似期望,跳过了从分布 \(p(z)\) 中采样的环节。
当我们可以直接从 \(p(z)\) 中采样时,我们可以用有限项求和代替积分来近似期望,即
\[E[f] \approx \frac{1}{L}\sum_{i=1}^L f(z^{(i)})\]但通常从 \(p(z)\) 中直接采样时不现实的,但给定 \(z_0\) 的值则可以容易地计算 \(p(z_0)\),一种替代方法是
\[E[f] \approx \sum_{i=1}^L f(z^{(i)})p(z^{(i)})\]我们可以发现,在这种替代方法中,求和项随着 z 维度的增加而指数级增加。而且从 z 空间中均匀采样的效率很低,我们希望从 \(p(z)\) 较大的区域或者说乘积 \(f(z)p(z)\) 较大的区域采样,这就引出了重要采样。
在重要采样中,我们还是引入一个参考分布 \(q(z)\),不过和拒绝采样不同,在这里我们假设 \(q(z)\) 和 \(p(z)\) 具有相同的性质,即 \(q(z) = \tilde{q}(z)/Z_q\) ,其中 \(\tilde{q}(z)\) 已知,\(Z_q\) 为未知的归一化常数,此时我们再写出期望的计算式
\[E[f] = \int f(z)p(z)dz = \frac{Z_q}{Z_p} \int f(z)\frac{\tilde{p}(z)}{\tilde{q}(z)}q(z)dz\]我们用 \(\tilde{r}_i =\tilde{p}(z^{(i)})/\tilde{q}(z^{(i)})\) 代替,用求和代替积分近似。可以从上式变换得到
\[E[f] \approx \frac{Z_q}{Z_p}\frac{1}{L} \sum_{i=1}^L \tilde{r}_if(z^{(i)})\]这个 \(\tilde{r}_i\) 也被称为重要权重 (importance weights),用于纠正从错误分布里采样引入的 bias。不同于拒绝采样,重要采样中所有的样本都会被保留。我们用保留的样本集计算比值 \(Z_p/Z_q\),有
\[\frac{Z_p}{Z_q} = \frac{1}{Z_q}\int \tilde{p}(z)dz = \int \frac{\tilde{p}(z)}{\tilde{q}(z)}q(z)dz \approx \frac{1}{L}\sum_{i=1}^L\tilde{r}_i\]将这个比值代入上面期望的计算式中,我们可以得到重要采样最终的期望计算式
\[E[f] \approx \sum_{i=1}^L w_if(z^{(i)}),\quad where \ w_i = \frac{\tilde{r}_i}{\sum_m \tilde{r}_m} = \frac{\tilde{p}(z^{(i)})/q(z^{(i)})}{\sum_m \tilde{p}(z^{(m)})/q(z^{(m)})}\]和拒绝采样类似,重要采样的近似程度取决于 \(q(z)\) 和 \(p(z)\) 之间的相似程度。通常 \(p(z)f(z)\) 集中分布在 z 空间下的很小一部分区域,对应我们的重要权重 \(\{r_i\}\) 中部分元素值很大,剩下的权重则没那么重要。因此有效的采样数实际上是远小于总采样数 L 的。重要采样还有一个缺点,即使期望估计效果非常差,也没有指标能看出来。所以这就要求参考分布 \(q(z)\) 不应该在 \(p(z)\) 较大的地方取值较小或为零
在拒绝采样中,我们需要找一个 k 值来约束 \(kq(z)\geq \tilde{p}(z)\) 。而在很多情况下,合适的 k 值也很难找。因为一旦 k 值过大就会导致低接受率。所以提出了 sampling-importance-resampling (SIR) 的方法,结合拒绝采样和重要采样来避免寻找合适 k 值
SIR 方法分为三个步骤
此时我们重采样得到的样本就可以近似认为是对 \(p(z)\) 的采样。可以证明,当 \(L \rightarrow \infty\) 时,重采样等同于对 \(p(z)\) 的采样。
我们首先写出重采样对应的分布
\[p(z\leq a) = \sum_{i:z_i\leq a} w_i = \frac{\sum_i I(z^{(i)}\leq a) \tilde{p}(z^{(i)})/q(z^{(i)})}{\sum_m \tilde{p}(z^{(m)})/q(z^{(m)})}\]当 \(L \rightarrow \infty\) 时,我们就可以用积分来代替求和,此时再写出分布
\[p(z\leq a) = \frac{\int I(z\leq a) \{\tilde{p}(z)/q(z)\} q(z)dz}{\int \{\tilde{p}(z)/q(z)\} q(z)dz} = \frac{\int I(z\leq a)\tilde{p}(z)dz}{\int \tilde{p}(z)dz} = \int I(z\leq a)p(z)dz\]可以发现这就是 \(p(z)\) 的累积分布函数
在上面的基础采样方法中,我们介绍了用于估计函数期望的两种采样方法。但无论是拒绝采样还是重要采样,都很难应用在高维度的数据上。所以我们在这里引入 Markov chain Monte Carlo (MCMC) 框架。
在之前的采样方法中,我们都是不断地从参考分布中独立采样。在 MCMC 中,我们维护对当前状态 \(z^{(\tau)}\) 的记录,参考分布 \(q(z\vert z^{(\tau)})\) 取决于当前的状态。所以我们的一系列采样 \(z^{(1)},z^{(2)},...,\) 构成了马尔科夫链。这里的采样分布和之前一致,\(p(z) = \tilde{p}(z)/Z_p\) ,其中 \(\tilde{p}(z)\) 已知,\(Z_p\) 为未知的归一化常数,我们首先来看最基础的 Metropolis 算法
使用该算法时,我们假设参考分布具有对称性,即 \(q(z_a\vert z_b) = q(z_b\vert z_a)\) 。在每轮算法迭代中,我们首先从参考分布 \(q(z\vert z^{(\tau)})\) 中生成采样 \(z^\star\) ,然后计算接受概率
\[A(z^\star,z^{(\tau)}) = min(1,\frac{\tilde{p}(z^\star)}{\tilde{p}(z^{(\tau)})})\]我们在 (0,1) 区间上取随机数 \(\mu\) ,如果 \(A(z^\star,z^{(\tau)}) > \mu\) ,那么我们就接受这个样本 \(z^\star\) 。
我们可以发现,如果从状态 \(z^{(\tau)}\) 转移到状态 \(z^\star\) 可以提升 \(p(z)\),那么这个样本点 \(z^\star\) 一定会被接受。根据马尔科夫的规则,如果我们接受了样本 \(z^\star\) ,则下一个状态会被设置为该样本,即 \(z^{(\tau+1)} = z^\star\) ;否则被设置为和当前状态一致,即 \(z^{(\tau+1)} = z^{(\tau)}\) 。接着迭代到了下一轮,从参考分布 \(q(z\vert z^{(\tau+1)})\) 中采样并以此类推
当我们对问题进行泛化,不再假设参考分布具有对称性时,此时的算法称为 Metropolis-Hastings 算法,A 的计算式发生了改变
\[A(z^\star,z^{(\tau)}) = min(1,\frac{\tilde{p}(z^\star)q(z^{(\tau)}\vert z^\star)}{\tilde{p}(z^{(\tau)})q(z^\star \vert z^{(\tau)})})\]Gibbs 采样可以看成是 Metropolis-Hastings 算法的一个特例。我们假设要采样的分布为 \(p(\mathbf{z}) = p(z_1,..,z_M)\) 。
我们首先选取初始状态。之后在算法的迭代中,每次将一个变量的值变成基于剩下变量的采样值。这里直接给出采样流程
记当前状态为 \(\tau\),在算法的一次更新中
接下来我们来看一下 Gibbs 采样与 Metropolis-Hastings 算法之间的联系。
在 Gibbs 采样中,我们每次根据剩下的状态元素来更新某个状态元素。对应的参考分布可以写为 \(q_k(z^\star \vert z) = p(z_k^\star \vert z_{\backslash k})\) ,此时我们写出对应的 A
\[A(z^\star,z) = \frac{p(z^\star)q_k(z\vert z^\star)}{p(z)q_k(z^\star \vert z)} = \frac{p(z_k^\star \vert z^\star_{\backslash k}) p(z^\star_{\backslash k}) q_k(z\vert z^\star) }{p(z_k \vert z_{\backslash k}) p(z_{\backslash k}) q_k(z^\star \vert z)} = \frac{p(z_k^\star \vert z^\star_{\backslash k}) p(z^\star_{\backslash k}) p(z_k\vert z^\star_{\backslash k}) }{p(z_k \vert z_{\backslash k}) p(z_{\backslash k}) p(z^\star_k \vert z_{\backslash k})}\]在 Gibbs 采样中,每次更新状态元素时,剩下元素没有发生改变。也就是说 \(z^\star_{\backslash k} = z_{\backslash k}\) ,代入上式可以发现
\[A(z^\star,z) = 1\]说明 Metropolis-Hastings 每一步迭代得到的样本都会被接受,对应 Gibbs 采样中的每一次更新
这里待采样分布和之前一致,\(p(z) = \tilde{p}(z)/Z_p\) ,其中 \(\tilde{p}(z)\) 已知,\(Z_p\) 为未知的归一化常数。
Slice sampling 的思想是引入另一个变量 u,抽取联合分布样本,从概率密度函数 p(x) 所围成的面积里均匀采样,即
\[\hat{p}(z,u) = 1/Z_p \ \ if \ 0\leq u \leq \tilde{p}(z); \quad else\ \ 0\]此时上式的边缘分布就是待采样分布,所以我们可以通过对 \(\hat{p}(z,u)\) 采样并忽略 u 来实现对 \(p(z)\) 的采样,证明如下
\[\int \hat{p}(z,u)du = \int_0^{\tilde{p}(z)} \frac{1}{Z_p}du = \frac{\tilde{p}(z)}{Z_p} = p(z)\]在进行采样时,可以采用 Gibbs 采样方法,即
从直观上来描述,我们首先在上一步选取了样本 \(z^{(t-1)}\) ,然后由均匀分布得到 \(u^{(t)} \sim U(0,\tilde{p}(z^{(t−1)})\) 。此后,再从 \(z^{(t)}\sim U(\tilde{p}(z)>u^{(t)})\) 采样。
我们可以发现,\(u\) 是很容易得到的,因为它的分布是确定的。而 \(z\) 的采样则不那么容易,它的分布甚至是不连续的。原论文中给出了两种启发式的解决方案来探寻不连续的 \(z\) 采样(这也就是叫做slice sampling的原因)。第一种方式叫 stepping out and shrinkage procedure,第二种方式叫 doubling procedure