hulu机器学习问题与解答系列|十四：如何对高斯分布进行采样

欢迎回到“采样”系列~

今天的内容是

【如何对高斯分布进行采样】

场景描述

高斯分布，又称正态分布，是一个在数学、物理及工程领域都非常重要的概率分布。在实际应用中，我们经常需要对高斯分布进行采样。虽然在很多编程语言中，直接调用一个函数就可以生成高斯分布随机数，但了解其中的具体算法能够加深我们对相关概率统计知识的理解；此外，高斯分布的采样方法有多种，通过展示不同的采样方法在高斯分布上的具体操作以及性能对比，我们会对这些采样方法有更直观的印象。

问题描述

如果让你来实现一个高斯分布随机数生成器，你会怎么做？

背景知识：概率统计

解答与分析

首先，假设随机变量z服从标准正态分布N(0,1)，令

x = σ·z + μ

则x服从均值为μ、方差为σ2的高斯分布N(μ, σ2)。因此，任意高斯分布都可以由标准正态分布通过拉伸和平移得到，所以这里我们只考虑标准正态分布的采样。另外，几乎所有的采样方法都是以均匀分布随机数作为基本操作，因此这里假设我们已经有均匀分布随机数生成器了。均匀分布随机数一般用线性同余法来生成（伪随机数），具体参见文献[1]。

常见的采样方法有逆变换法 (Inverse Transform Method)、拒绝采样法 (Rejection Sampling)、重要性采样及其重采样 (Importance Sampling, Sampling-Importance-Resampling)、马尔科夫蒙特卡洛采样法 (Markov Chain Monte Carlo) 等。具体到高斯分布，我们需要如何采样呢？ 

如果直接用逆变换法，基本操作如下：

上述逆变换法需要求解erf(x)的逆函数，这并不是一个初等函数，没有显式解，计算起来比较麻烦。为了避免这种非初等函数的求逆操作，Box-Muller算法采用如下解决方案：既然单个高斯分布的累计分布函数不好求逆，那么两个独立的高斯分布的联合分布呢？假设x, y是两个服从标准正态分布的独立随机变量，它们的联合概率密度为：

考虑(x, y)在圆盘(x, y) | x2 + y2 ≤ R2上的概率：

通过极坐标变换将(x, y)转化为(r, θ)，可以很容易求得上述二重积分：

这里F(R)可以看成是极坐标中r的累积分布函数。由于F(R)的计算公式比较简单，逆函数也很容易求得，所以可以利用逆变换法来对r进行采样；对于θ，在[0, 2π]上进行均匀采样即可。这样就得到了(r,θ)，经过坐标变化即可得到符合标准正态分布的(x, y)。具体采样过程如下：

Box–Muller算法由于需要计算三角函数，相对来说还是比较耗时，而Marsaglia polar method则避开了三角函数的计算，因而更快，其具体采样操作如下：

除了逆变化法，我们还可以利用拒绝采样法，选择一个比较好计算累积分布逆函数的参考分布来覆盖当前正态分布（可以乘以一个常数倍），进而转化为对参考分布的采样以及对样本点的拒绝/接收操作。考虑到高斯分布的特性，这里可以用指数分布来作为参考分布。指数分布的累积分布及其逆函数都比较容易求解。由于指数分布的样本空间为x≥0，而标准正态分布的样本空间为(-∞, +∞)，因此还需要利用正态分布的对称性来在半坐标轴和全坐标轴之间转化。具体来说，取λ=1的指数分布作为参考分布，

实际应用时，M需要尽可能小，这样每次的接受概率大，采样效率更高。因此，可以取

因此，具体的采样过程如下：

拒绝采样法的效率取决于接受概率的大小：参考分布与目标分布越接近，则采样效率越高。有没有更高效的拒绝采样算法呢？这就是Ziggurat算法，该算法本质也是拒绝采样，但采用多个阶梯矩形来逼近目标分布（如图1所示）。Ziggurat算法虽然看起来稍微繁琐，但实现起来并不复杂，操作也非常高效，其具体采样步骤可以参考文献[2]。

图1  Ziggurat算法

（图片来源于Numerical Computing with MATLAB第九章Random Numbers）

总结与扩展

高斯分布的采样方法还有很多，我们只列举了几种最常见的。具体面试时，候选人不需要回答所有的方法，知道其中一两种即可，面试官可以针对这一两种方法深入提问，如理论证明、优缺点、性能等。如果候选人没有思路，面试官可以引导其回忆那些通用的采样方法，如何将那些策略用到高斯分布这个具体案例上。另外，本题还可以适当扩展，把一般的高斯分布换成截尾高斯分布 (Truncated Gaussian Distribution) ，如何采样？如果是高维随机变量，拒绝采样法会存在什么问题？怎么解决呢？

参考文献：

[1] Linear congruential generator,

https://en.wikipedia.org/wiki/Linear_congruential_generator

[2] Ziggurat algorithm,

https://en.wikipedia.org/wiki/Ziggurat_algorithm

下一题预告

【多层感知机与布尔函数】

场景描述

神经网络概念的诞生很大程度上受到了神经科学的启发。生物学研究表明，大脑皮层的感知与计算功能是通过分多层实现的，例如视觉图像，首先光信号进入大脑皮层的V1区，即初级视皮层，之后依次通过V2层，V4层，即纹外皮层，进入下颞叶参与物体识别。深度神经网络，除了其模拟人脑功能的多层结构，最大的优势在于能够以紧凑简洁的方式来表达比浅层网络更复杂的函数集合（这里的“简洁”可定义为隐层单元的数目与输入单元的数目呈多项式关系）我们的问题将从一个简单的例子引出，已知神经网络中每个节点都可以进行“逻辑与／或／非”的运算，如何构造一个多层感知机 (Multi-Layer Perceptron, MLP) 网络实现n个输入比特的奇偶校验码（任意布尔函数）？

问题描述

1. 如何用多层感知机实现一个异或逻辑（仅考虑二元输入）？

2. 如果只使用一个隐层，需要多少隐节点能够实现包含n元输入的任意布尔函数？

3. 上面的问题中，由单隐层变为多隐层，需要多少节点？

4. 合理配置后所需的最少网络层数是多少？