机器学习中的概率模型和概率密度估计方法及VAE生成式模型详解之四（第2章）

Posted 2020-11-19

第2章  分布密度估算的采样法

   

2.1 如何采样

   

有了概率分布P(X)，我们要生成符合这个分布的样本。这个过程也叫采样(sampling)。

   

这其实并不是一个简单的任务。我们知道，计算机可以生成符合某些指定分布的伪随机数（伪随机数即样本），但并不支持任意的分布。对于所支持的非均匀分布的伪随机数x的产生，我们需要先产生均匀分布的伪随机数y，然后通过某种变换，把y转换成指定的非均匀分布的x。例如下图，y轴上是一些均匀分布的点，通过某个函数Fx映射到x轴，就得到非均匀分布的点了。如果函数Fx足够复杂，我们可以把任意一个分布转换成另一个分布。如何在两个分布之间找到这样的函数，这是VAE的关键点之一。

   

均匀分布的伪随机数的生成常用两个算法，它们都是用数列迭代的方法产生伪随机数。一个叫线性同余序列(Linear Congruential Generator, LCG)，公式是

   

   

来自https://en.wikipedia.org/wiki/Inverse_transform_sampling

   

产生[0, M]之间的一个随机数。a或b，以及M要足够大才能有好的结果。早期的Matlab使用的是这个算法：a=7⁵，b=0， M=2³¹ – 1。

   

n₀是初始数，称为"种子"(seed)。若要重现实验结果，在产生随机数时要提供使用的seed。

   

另一个叫Lagged Fibonacci Generator，公式是

   

   

新版的Matlab (>R5.0) 使用的是这个算法 ：a=1，b=1，? =17，?=5。

   

均匀分布到非均匀分布的转换方法也有两种，一种可以称为CDF Transformation，另一种可以称为PDF Transformation。

   

前者需要知道随机变量的CDF(Cumulative Distribution Function)，即离散变量

   

   

连续变量

   

先在均匀分布上采集一个随机变量u，则

   

   

就是符合Fx分布的随机变量，即Z具有与X相同的分布。F_x^-1是F_x的逆函数。

   

证明过程很简单。Z的CDF是F_Z(z)，

   

   

大括号内的不等式两边加上F_X，则左边变成u，右边变成F_X(z)。CDF具有单调性，所以上面的不等式加上F_X 也不改变不等号的方向。因此，

   

证毕。

   

显然，只有当F_x^-1有解析解时，CDF Transformation才是可行的。

   

比如，

则

   

如果我们知道的是x的概率密度函数f(x)，为了使用上述方法，我们要先计算它的CDF：

   

这一步也无法保证有解析解。下面是一个有解析解的例子：

   

则

   

第二种方法是PDF Transformation，直接利用概率密度函数（PDF）推导出函数h，使得X = h(U)把均匀分布的样本U转换成需要的X的分布。这个方法避免求积分，但需要求h^-1的雅可比矩阵，看上去也不是一个容易使用的方法。我们常用的标准高斯分布?(x) ~ N(0, 1)的伪随机数的产生方法是Box-Muller公式，其推导就使用了PDF Transformation，但用了一个巧妙的方法，即极坐标方法，避开了实际求雅可比矩阵，这里就不展开了。下图直观地解释了这种方法的原理：每个圈（从浅绿色到深绿色）里的数据点的数量基本一样（在极坐标里是均匀分布的），但在x-y坐标系里就是高斯分布的了（点的密集度从中央到边缘越来越低，其实在x-轴和y-轴上的分布都是高斯分布）。

   

Polar method is sometimes more convenient to think of generating a pair of random variates (x, y) from the target distribution p(z). We can map this pair to a representation in polar form (rcosθ, rsinθ), which exposes other mechanisms for sampling,   e.g., the famous Box-Muller transform is  derived in this way.

   

来源：http://blog.shakirm.com/2015/10/machine-learning-trick-of-the-day-4-reparameterisation-tricks/

   

对于标准正态（高斯）分布N(0, 1)，另一种采样方法是利用中心极限定理，从二值分布（Bernoulli分布）大量地采集样本x₁，x₂，…，x_n，然后计算出一个高斯随机数：

   

   

问题是，这个方法需要很长的计算时间，因为每个高斯随机数的产生都要先产生n个Bernoulli分布的随机数。因此一般我们还是用Box-Muller公式，直接把均匀分布转换成标准高斯分布。

   

至此，我们知道如何对一些常见的有公式表达的概率分布进行采样了。我们记住这些方法的思路，就是先得到均匀分布的样本，再通过一个函数把这些样本转换成所需的分布的样本。进一步延伸，我们应该可以把从已知分布（不一定是均匀分布）采样的样本，通过一个函数，转换成我们需要的任意分布。这里，我们的任务就是找到这个函数。

   

2.2    蒙特卡罗方法

   

那么，如果概率分布f(x)是任意的，同时CDF Transformation和PDF Transformation都很难或无法找到解析解呢？这时我们一般使用一个直截了当的采样方法，名为MCMC（Markov Chain Monte Carlo）的采样方法，中文翻译成马尔科夫链-蒙特卡罗方法，就是在马尔科夫链上跑蒙特卡罗方法，所以它有两个部分，即蒙特卡罗方法和马尔科夫链。

   

先简要介绍一下蒙特卡罗方法（这时你的头脑中是否闪现过一个念头：也许还有拉斯×××方法？）。这个方法的最初目的是求定积分：

   

   

其中，?是定积分的x取值域。这个积分问题的难处在于函数f(x)很复杂，或者x的维度很高，使得这个定积分没有或很难找到解析解，常规的数值方法也不适用。这时我们可以用采样的方法来求I的近似值。先看个x为一维的简单例子（下图来自Wikipedia） ：

我们来求?的近似值。方法是在上图的正方形里撒均匀分布的N个点（红+蓝），计算出落入1/4圆里面有C个点（红）。正方形的面积是1，则圆的面积是

   

4C/N =? r² = ?

   

f(x)

对单一的一维变量x的函数f(x)做定积分，就是求f(x)的曲线下的面积。无论f(x)的曲线多么复杂，我们都可以把需要积分的这段放在一个正方形里，再用上述方法计算曲线下的面积。假如正方形的面积是S，N个均匀采样点{(x_i, y_i), i=1..N}中有C个点满足

   

   

则f(x)下的面积约是

   

   

当N足够大时，这个面积就近似于f(x)的定积分。这就是基本的蒙特卡罗方法。

   

这个方法需要在坐标系的横轴和纵轴方向以正方形为界，采集均匀分布的多对点{(x_i, y_i), i=1..N}。如果x是多维变量，那就需同时在各个维上采集均匀样本点。均匀分布的采样我们在上面已经介绍了。

   

这种计算定积分的近似值的方法称为Rejection Sampling。意思是，N个采样点中，只有落在f(x)曲线下的区域里的C个样本点才被接受（上图中红色的点），在这个区域之外的点（蓝色的点）被拒绝。

   

那么，定积分的这种近似解法和任意分布f(x)的采样有什么关系 ？因为CDF的计算公式是

   

我们可以用"拒绝"的方法把

内的均匀采样得到的样本点{(x_i, y_i), i=1..N}检查一遍，只留下满足

的那些x_i点作为符合f(x)分布的样本点，其它的点都被拒绝。这样，在f(x)越大的地方，留下的样本点越多，因为这些地方的f(x)下的面积也大。最后留下的{x_i}就是合法的样本点。这也算是一种从均匀分布到任意分布的转换方法。

   

这个方法的缺点是，如果那个正方形与f(x)之间的间隙比较大，比如当f(x)的曲线有 "陡峭的山峰"形状，则被拒绝的采样点也很多，造成较大的浪费（获得这些采样点所花的计算资源和时间）。

   

提高效率的一个方法是Importance Sampling，意思是，f(x)取值大的地方（即"尖锐的山峰"处）所取样本的点数应该多于f(x)取值小的地方，即越重要的地方，样本点应该越多，越密集。这也是概率密度的本意。下面是一个例子，来自http://bjlkeng.github.io/posts/markov-chain-monte-carlo-mcmc-and-the-metropolis-hastings-algorithm/。

   

图中蓝线的Target（即f(x)）是双伽马分布(Double Gamma Distribution)。可以看到，它有尖锐的"山峰"。如果我们围绕它画一个矩形（假如我们知道f(x)的峰值），并在这个矩形里进行均匀采样，那么被拒绝的点会非常多（落入蓝线条以上的区域的点）。但如果我们能用一个Envelope函数q(x)覆盖f(x)，然后在横坐标方向以q(x)为分布采样x_i，在纵坐标方向继续以均匀分布采样y_i，那么被拒绝的点就会少很多，因为大多数x_i被q(x)集中在中间"山峰"的位置。当然，q(x)的选择是个关键，它必须满足两个特点，一是它足够简单，容易采样，二是能大致符合f(x)的形状，使q(x)和 f(x)之间的空隙尽可能地小。图中，黑色虚线q(x)是高斯分布N (0, 2)。在f(x)的取值大的区域，从q(x)生成的样本x_i也聚集得多。前面的内容已经讲过如何从高斯分布采样（Box-Muller公式）。这样，我们可以用一个简短的程序就能实现双伽马分布的Importance Sampling。

   

注意，来自q(x)的样本x_i被接受的概率等于f(x_i)/q(x_i)，即f(x_i)和q(x_i)的高度之比。这是因为在纵向的采样是在0 ~ q(x_i)之间的均匀采样，样本点落入f(x_i)以下部分的概率是f(x_i)/q(x_i)。请仔细观察图中的绿虚线和红虚线之间的关系。Importance Sampling也是一种Rejection Sampling，只不过在横轴方向的采样不再是均匀分布的，而是按照q(x)的分布采样。

   

让我们再用数学语言来描述一下Importance Sampling。

   

我们要求的积分

   

最后这个等号来自期望值E的定义。我们设

   

这样上式简记成

   

   

类似地，如果x是离散随机变量，可以求得

   

   

这样，我们把求积分（或求和）转化为求期望值。注意E的下标，它说明x的分布必须符合q(x)，或者说，x来自q(x)。那么如何求期望值呢？

   

根据大数定理，如果我们通过q(x)采集大量的i.i.d.样本x_i，i =1..N，那么我们可以用g(x_i) 的均值近似g(x) 的期望值：

   

其中，

   

   

而f(x)和q(x)都是已知的，计算所有样本点的g(x_i)的均值就可以得到I的近似解。这就是Importance Sampling的数学表达。

   

上面的公式中，g(x) = f(x)/q(x)。前面说过，g(x)的意义是，它是来自q(x)的样本x被接受的概率。

   

我们也可以解释

   

的意义：x_i 被选中的概率是q(x_i)，它被接受的概率是g(x_i)，那么，x_i 最终成为合法的样本点的概率f(x)是它被选中并被接受的概率，即这两个概率的乘积q(x_i)g(x_i) 。

   

同时，大数定理还告诉我们，在N是有限值时，均值到期望值之间的误差是：

其中，

此项即方差（variance）；而N(0,1)是标准高斯分布。

   

注：大数定理（Law of Large Numbers）的另一个名字是我们前面提到过的中心极限定理（Central Limit Theorem），前者可能在信息论里较常使用，后者在概率论里较常使用。

这个定理是概率论里最本质、最核心的定理，但它的证明比较复杂，思路是先把多个随机变量之和用其特征函数变换一下（关于特征函数，可参阅https://en.wikipedia.org/wiki/Normal_distribution里面的相关内容），再用泰勒级数展开特征函数，最后利用高斯函数的傅立叶变换就是其本身这个性质完成证明。这里就不再进一步展开了。

   

显然，当g(x) = I 时，误差趋于0。此时，

   

此q(x)即全部符合条件的q(x)中的最佳q(x)函数，记为q*(x)我们在前面说过，q(x)最好能大致符合f(x)的形状，现在这句话得到了数学解释 --- 这样的q(x)会使方差尽量小，亦即使误差尽量小。

   

如果能找到q*(x)，则我们只需一个样本（即N=1）就能得到误差为0的估计值。之所以如此，是因为q*(x)已经解决了最初的求积分的问题，而这是无法做到的（所以才需要估值的方法）。既然q*(x)是无法得到的，我们就要专注找其它能使方差减小的q(x)。在后面讲MCMC算法和Variational Inference算法时，我们还要回到找最佳q(x)这个话题。

   

q(x)对提高采样效率具有重大影响，选的不好会使采样效率变低。

   

从g(x)= f(x)/q(x)可以看出，如果q(x)产生的样本使g(x)很大，则方差也很大，而这种情况会发生在当q(x)很小时。所以我们希望q(x)大一些。但另一方面，如果

   

   

则样本的拒绝率会很高。所以q(x)也不能太大。总的来说，x的维数越高，好的q(x)越难找，近似算法的误差也越大，为减小误差所需的样本也越多。这里从另一个角度揭示了x的维数与模型学习所需样本数之间的关系。

   

即便找q(x)存在风险，但Importance Sampling让我们以较高的效率采集符合已知概率分布的样本。而且，它利用大数定理，把求积分（或求和）转换成求期望值，再转换成求大量样本的均值这个思想，在机器学习的很多算法，包括深度学习的很多算法中被广泛使用，而且上面对q(x)的讨论可以帮助读者深刻理解VAE和其它生成式模型。

   

2.3   马尔科夫链

   

在随机变量的维数比较高的时候，好的q(x)是很难找到的，这时马尔科夫链就出场了。

一个马尔科夫链定义了一个状态转换模型

   

其中x和x' 表示状态，x→x'  表示状态x转化到另一个状态x' ，状态转换概率是T。所以，对所有的x，必须有

   

   

下面这个例子和图来自Wikipedia，是一个股市预测模型。

   

   

这相当于有一个随机变量market，它有三个取值，即用圆圈表示的三个状态：x⁰='Bull market' ，x¹='Bear market'和x²='Stagnant market'。而圆圈之间的有向连线表示状态之间的转换关系，旁边标注的数字是该转换的概率T。一个状态可能保持不变，所以，一个圆圈到自己也有"状态转换"概率。

   

特别提醒一下，圆圈代表的不是随机变量，而是随机变量的状态（即可能的取值）。如果一个模型里包含多个随机变量，或随机变量是多维向量，那么一个状态就是这些变量或向量各维的某个联合取值。

   

用X表示状态，x表示状态的取值，P表示某个状态出现的概率，T表示状态转换概率，t表示时间步数（time step）。每个time step状态转换一次。则新状态的概率更新为：

   

   

比如上图，

   

   

这就是说，当前时刻(t+1时)的状态只由上一时刻(t时)的状态所决定，与更早或更晚的状态无关。这是马尔科夫链的关键思想。

   

注：虽然这里的定义和例子是针对离散随机变量的，但也同样适用于连续随机变量，只不过求和变成求积分。以后不再特别说明。

   

在马尔科夫链上，如果在有限的步数内，一个任意状态到另一个任意状态的转换概率大于0，则称这个马尔科夫链是Regular Markov Chain，也叫各态历经的马尔科夫链。这种马尔科夫链有个性质，就是状态不断地转换，最终会进入一个唯一的稳态（stationary state），此时P(X)不再变化。我们用f(X)取代P(X)来表示稳态时的各个状态出现的概率。如果T选择得合适，那么可以证明，f(X)正是我们要从中采样的目标概率。此时继续运行马尔科夫链，即继续进行状态转换，那么，每个状态出现的次数与它的概率f(X)成正比。也就是说，每次的状态转换相当于一次从f(X)的采样。这就是MCMC的核心。

   

模型处于稳态时，满足条件：

   

   

不过，我们会利用一个更强的称为Time Reversibility的条件：

   

   

即状态从x' 转换到x与从x 转换到x' 具有相同的可能。这种状态也称为Detailed Balance。要想使稳态时的各状态的概率（或分布）是我们希望的f(x)，关键是找到合适的状态转换概率T。

   

2.4    MCMC的Metropolis-Hastings算法

   

如果我们任凭马尔科夫链完全跟随概率T更新状态，那么，在X的整个状态空间中，从一个状态转换到相距较远的另一个状态可能是非常困难的，因为状态转换都有本地化的倾向，即相邻的状态之间较容易转换，结果是，出现相距当前状态较远的某些状态很难出现，或者说，要经过大量的采样才会出现，无论它们的概率是什么。因此，采样的效率不高。解决办法是，让Importance Sampling运行在马尔科夫链上。这就是MCMC的Metropolis-Hastings算法。我们将看到，这个算法一举两得，即解决了状态转换概率T的选择问题，又解决了下一个状态的选择问题（即找好的q(x)的问题）。

   

Metropolis-Hastings算法的关键之处在于，它把状态转换概率T分解成两部分：

   

   

其中，q是proposal概率（当x是连续的，q是分布函数，下同），比如q(x→x' )表示以当前状态x为起点，建议跳到下一状态x' 的概率，或者说x' 被选中成为新状态的概率。q的作用与Importance Sampling中的q(x)是一样的；g是接受q的建议的概率，也与Importance Sampling中的q(x)有同样的作用。状态转换概率T(x→x')等于x' 被选中的概率乘上x' 被接受的概率，这和Importance Sampling同样是一致的。因此，找合适的T变成找合适的q，而g可由q和希望的概率分布f决定。

   

把上面两个等式代入Detailed Balance的条件

   

并整理，得到

   

   

结果是，T消失了，而我们可以自由地选择q(x)或g(x)，只要保证上面的式子成立就行。

   

上式如果分子大于或等于分母，则等号右边≥1，此时

   

那么从当前状态x转换至q(x→x' )建议的状态x' 应该被接受，因为f(x' )也在上式的分子里，它的值很可能处于f(x)之上的高处，而我们以前说过，在概率密度高的地方最好也更密集地采集样本。由于接受概率不能大于1，我们取g(x→x' )=1。

   

如果分子小于分母，则f(x' )很可能但不肯定处于f(x )之下的低处。我们希望低处的样本少于高处的样本，但我们不能完全拒绝x' ，而是以一个概率接受x' 。我们取分母部分g(x' →x)=1，则接受x' 的概率

   

这时Detailed Balance仍然成立，同时状态转换得最快（请你想想为什么）。

   

归纳一下上面的两种情况，Metropolis-Hastings算法使用如下的x→x' 的接受概率：

   

   

这里需要问一个问题，在马尔科夫链上跑Importance Sampling算法与常规的Importance Sampling有什么不同或者说改进吗？毕竟，MCMC也没有告诉我们如何选择最关键的q(x)。

   

它们之间重要的差别在于它们的g(x)，当然q(x)也有差别，但这个差别比较隐秘。

   

如果我们为MCMC选q(x)为均匀分布，就如同在一般的Rejection Sampling里那样，那么，

   

   

因为

   

如果q(x)选为高斯分布，注意，它与在一般的Importance Sampling里使用的均值为0的高斯分布有所不同，它是以当前状态为均值（即中心值），以方差作为超参数的非标准高斯分布。此时，由于高斯分布也具有对称性，

   

   

所以，

总的来说，当q(x)具有对称性时，

   

   

这是Metropolis-Hastings算法的核心结果。

   

可见，新的样本被接受的概率取决于它与当前样本点的相对高度，即f(x' )与f(x)之比，比较的基准是f(x)。

   

而在一般的Importance Sampling中，g(x' )= f(x' )/q(x' ) ，比较的基准是q(x' )，且与当前样本点没有关系，对每个新样本点都是独立地进行接受与否的判断。

   

显然，从比较的基准角度看，如果q(x)选为高斯分布，则多数情况下x→x' 跳的不太远，MCMC的效率高于一般的Importance Sampling。同时，x→x' 也有一定的可能跳得较远，即仍然有可能遍历整个状态空间，不至于总是在某些相邻状态之间徘徊，有助于提高采样效率。而且，MCMC的最大的好处是，我们不必知道f(x)具体是什么，而只要知道两个样本点之间的相对高度就可以决定是否接受新的样本点。

   

MCMC的问题是，它必须从初始状态运行一段时间后才能进入稳态，进入稳态后才能开始采样。这段时间叫mixing time。我们无法知道MCMC系统什么时候进入了稳态。虽然有一些迹象可作为提示，但获得这些迹象也并非简单明了，而且也不能完全确定。

   

MCMC的另一个问题是，即使进入稳态，相邻的采样样本之间很可能是有关联的，而我们希望采集的样本都是i.i.d.的。解决的方法是，每间隔多次（比如1000次）采样再取一个样本，或者同时跑多个（比如100个）独立的马尔科夫链，在它们上面轮流采样。

   

无论用什么解决办法，这两个问题都会影响MCMC的采样效率和准确性。

   

总结一下采样，我们走了"均匀分布的采样→高斯分布的采样→任意分布的采样"这样一条路。

   

2.5    Bayesian Inference

   

现在，我们把MCMC用于Bayesian Inference，以获得分布参数的后验概率。

   

根据MCMC方法，我们先用先验概率P(?)以?_current为中心值在其附近产生一个新的参数?_new。P(?)相当于MCMC里的proposal分布。正像以前那样，我们一般选P(?)为对称的分布，而且一般用高斯分布。它的均值是当前参数?_current，方差则是一个超参数。方差的选择对最后结果的误差和MCMC算法的收敛速度有很大的影响，不过不是我们在这里要讨论的话题。

   

根据Metropolis-Hastings算法，这个新的参数是否被接受取决于目标分布（此处即后验概率）在?_new和?_current时的比例：

   

结果就是在?_new和?_current下的似然值之比（对数似然值之差）。如果这个比值大于等于1（对数似然值之差大于等于0），那么?_new就被接受；如果小于1，则以这个比值为概率接受?_new，未被接受时?_new = ?_current。接着使?_new成为新的?_current。

   

以?_new作为新的?_current，不断重复上述过程足够多次，最后当MCMC进入稳态时，我们可以开始为?采样。因为我们把节点状态概率设定为后验概率P(?|D)，所以在链上采集的样本符合这个概率。统计这些样本的empirical estimate就可以算出近似的后验概率。如果运行的时间足够长，MCMC也可以获得任意精度的后验概率。

   

但是，正如前面说过，采样法的计算时间比较长，尤其是当数据集很大或者要近似的模型比较复杂时（模型复杂时后验概率的状态空间会很大，遍历需更长时间）。下一章介绍的优化法，可以求得概率分布（或其参数）的近似值，以精度换取时间。