文章目录

前言The Box–Muller transformThe Ziggurat algorithm(金字形神塔)附录：Inverse transform sampling直观解释

前言

在Monte Carlo模拟技术中，许多地方都需要用到符合标准正态分布(高斯)的随机数来设计采样方案，因此了解如何用均匀分布随机数(实际上是均匀分布的伪随机数)来生成标准正态分布的随机数十分重要。本文将对这个最基本的问题做讨论，并提供c++11代码。

在讨论更高效的算法之前，我们先来看看能不能基于中心极限定理来设计算法。中心极限定理告诉我们，对于一组i.i.di.i.di.i.d的随机数{xk}∼U(μ,σ2)\{x_k\}\sim U(\mu,\sigma^2){xk​}∼U(μ,σ2)，有1n∑i=1nxi→N(μ,σ2/n)\frac{1}{n}\sum_{i=1}^n x_i \rightarrow N(\mu,\sigma^2/n)n1​∑i=1n​xi​→N(μ,σ2/n)。这个算法有两个问题：

计算量大。生成一个数需要用n个均匀分布随机数。无法准确刻画正态分布的末端效应。生成的数均有界，不会是很大的数。

此外，如果用“拒绝采样(rejection sampling)“的思路，在覆盖f(x)=ce−x2/2f(x)=ce^{-x^2/2}f(x)=ce−x2/2的矩形内均匀投点，保留曲线下的点，则计算较复杂(exp函数)，而且舍弃点多代价大。

我们下面介绍两种更高效的算法：The Box–Muller transform 和 The Ziggurat algorithm。它们一个是对分布函数做了变换，另一个还是使用了“拒绝采样“的思路，但并不是简单的仅用一个大矩形覆盖。

The Box–Muller transform

The Box–Muller transform 把一对均匀分布随机数映射到一对标准正态分布随机数。它有两种形式：

基本形式：用(0,1)(0,1)(0,1)均匀分布随机数，需要计算三角函数sin⁡\sinsin和cos⁡\coscos；极坐标形式：用(−1,1)(-1,1)(−1,1)均匀分布随机数，且不需要计算三角函数。

我们希望计算积分I=∫−∞∞e−x2/2dxI=\int_{-\infty}^{\infty}e^{-x^2/2}dxI=∫−∞∞​e−x2/2dx，可以先取它的平方并用极坐标表示，I2=∫−∞∞∫−∞∞e−(x2+y2)/2dxdy=∫02π∫0∞rer2/2drdθ.I^2=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-(x^2+y^2)/2}dxdy=\int_{0}^{2\pi}\int_{0}^{\infty}re^{r^2/2}drd\theta.I2=∫−∞∞​∫−∞∞​e−(x2+y2)/2dxdy=∫02π​∫0∞​rer2/2drdθ.可以看到极角θ\thetaθ服从(0,2π)(0,2\pi)(0,2π)均匀分布，径向距离rrr服从rer2/2re^{r^2/2}rer2/2分布函数(即r2r^2r2服从χ2分布\chi^2分布χ2分布)。如果把rrr的累积分布函数(cumulative distribution function)写出来F(x)=∫0xrer2/2dr=1−e−x2/2,F(x)=\int_0^xre^{r^2/2}dr=1-e^{-x^2/2},F(x)=∫0x​rer2/2dr=1−e−x2/2,我们就可以通过FFF的逆函数F−1(u)F^{-1}(u)F−1(u)，将均匀分布u∼U(0,1)u\sim U(0,1)u∼U(0,1)映射到目标分布，即x=−2ln⁡(1−u)x=\sqrt{-2\ln (1-u)}x=−2ln(1−u)​，也等价于x=−2ln⁡ux=\sqrt{-2\ln u}x=−2lnu​ (注意到如果u∼U(0,1)u\sim U(0,1)u∼U(0,1)，那么1−u1-u1−u也是)。

基本形式：

假设U1U_1U1​和U2U_2U2​是区间(0,1)(0,1)(0,1)上的均匀分布随机数，令

Z1=Rcos⁡(θ)=−2ln⁡U1cos⁡(2πU2)Z_1=R\cos(\theta)=\sqrt{-2\ln U_1}\cos(2\pi U_2)Z1​=Rcos(θ)=−2lnU1​​cos(2πU2​)，

Z2=Rsin⁡(θ)=−2ln⁡U1sin⁡(2πU2)Z_2=R\sin(\theta)=\sqrt{-2\ln U_1}\sin(2\pi U_2)Z2​=Rsin(θ)=−2lnU1​​sin(2πU2​)，

则Z1Z_1Z1​和Z2Z_2Z2​是两个独立的标准正态分布随机变量。

极坐标形式：

转换到极坐标的方法叫做Marsaglia polar method。该方法对水平和竖直方向(−1,1)(-1,1)(−1,1)正方形区域内均匀投点，把落在单位圆内的点(x,y)(x,y)(x,y)通过s=x2+y2s=x^2+y^2s=x2+y2映射到单位圆周上，即(xs,ys)(\frac{x}{\sqrt{s}},\frac{y}{\sqrt{s}})(s​x​,s​y​)，这样就可以不用计算三角函数。由于s≡r2s\equiv r^2s≡r2独立于极角，且在区间(0,1)(0,1)(0,1)上均匀分布（因为∫∫dxdy=∫∫rdrdθ=∫∫12dr2dθ\int\int dxdy=\int\int rdrd\theta=\int\int \frac{1}{2}dr^2d\theta∫∫dxdy=∫∫rdrdθ=∫∫21​dr2dθ），因此可以用这个sss继续得到径向距离−2ln⁡(s)\sqrt{-2\ln(s)}−2ln(s)​。于是我们有：

Z1=−2ln⁡U1cos⁡(2πU2)=−2ln⁡s(us)=u⋅−2ln⁡ss,Z_1=\sqrt{-2\ln U_1}\cos(2\pi U_2)= \sqrt{-2 \ln s} \left(\frac{u}{\sqrt{s}}\right) = u \cdot \sqrt{\frac{-2 \ln s}{s}},Z1​=−2lnU1​​cos(2πU2​)=−2lns​(s​u​)=u⋅s−2lns​​,

Z2=−2ln⁡U1sin⁡(2πU2)=−2ln⁡s(vs)=v⋅−2ln⁡ss.Z_2=\sqrt{-2\ln U_1}\sin(2\pi U_2)= \sqrt{-2 \ln s}\left( \frac{v}{\sqrt{s}}\right) = v \cdot \sqrt{\frac{-2 \ln s}{s}}.Z2​=−2lnU1​​sin(2πU2​)=−2lns​(s​v​)=v⋅s−2lns​​.

极坐标形式对应的c++代码如下（顺带说一下，c++标准函数库libstdc++里std::normal_distribution用的就是Marsaglia polar method，见这里）：

//// Marsaglia polar method. C++11 代码// 由于生成的是一对i.i.d.高斯变量，每次生成的两个样本都可以利用起来！//bool _M_saved_available=false;if (_M_saved_available) //输出下面暂时储存的x*mult{_M_saved_available = false;return = _M_saved;}else{double x, y, r2;do{x = 2.0 * rand() - 1.0;y = 2.0 * rand() - 1.0;r2 = x * x + y * y;}while (r2 > 1.0 || r2 == 0.0);const double mult = std::sqrt(-2 * std::log(r2) / r2);_M_saved = x * mult;_M_saved_available = true; //暂时储存x*mult，输出y*multreturn = y * mult;}

算法局限性：尾部截断(Tails truncation)

由于计算机有数值精度限制，其表示的数字不可能无线接近0。比如，使用32bit精度，那么最接近0的数字是2−322^{-32}2−32。因此，当u1u_1u1​和u2u_2u2​都取2−322^{-32}2−32时，xxx的最大值是−2ln⁡(2−32)cos⁡(2−32)≈6.6\sqrt{-2\ln(2^{-32})}\cos(2^{-32})\approx6.6−2ln(2−32)​cos(2−32)≈6.6，即产生的随机数最大不会超过6.66.66.6个标准差，这意味着会有2×10−112×10^{-11}2×10−11的比例(尾部)会因为计算机精度问题无法采样到，这在大规模生成高斯随机数和研究rare events的时候要特别注意。

The Ziggurat algorithm(金字形神塔)

The Ziggurat algorithm 属于rejection sampling的一种。这个方法适用于分布函数是单调递减的情况，而正态分布的正半轴就属于这种情况。

附录：Inverse transform sampling直观解释

数学解释：假设F(x)=Pr[X&lt;x]F(x)=Pr[X&lt;x]F(x)=Pr[X<x]是CDF，并假设u∼U(0,1)u\sim U(0,1)u∼U(0,1)，定义映射X=f(x)X=f(x)X=f(x)，这个映射是[0,1]→R[0,1]\rightarrow R[0,1]→R。我们希望找到一个映射使得XXX恰好服从F(x)F(x)F(x)的CDF分布，即Pr[f(U)&lt;x]=F(x)Pr[f(U)&lt;x]=F(x)Pr[f(U)<x]=F(x)。显然，f=F−1f=F^{-1}f=F−1，因为Pr[X&lt;x]=Pr[f(u)&lt;x]=Pr[u&lt;F(x)]=F(x)Pr[X&lt;x]=Pr[f(u)&lt;x]=Pr[u&lt;F(x)]=F(x)Pr[X<x]=Pr[f(u)<x]=Pr[u<F(x)]=F(x)。

图像解释：我们假设随机变量取值是离散的，比如一共有5种取值。右边是CDF，相当于把左边的概率“盒子“逐渐堆积起来，直到最后一列对应的FFF正好等于1。这时候对F(x)F(x)F(x)均匀采样的话，相当于在最后一列按照各个盒子面积的相对大小（各部分的概率大小）来采样，F(x)F(x)F(x)的逆函数只是用来回溯到对应的xxx而已。随机变量取连续情况同样适用上述论证过程。

[1]:

[2]: /jmcmanus/pagedown-extra “Pagedown Extra”

[3]: http://meta./questions/5020/mathjax-basic-tutorial-and-quick-reference

[4]: http://bramp.github.io/js-sequence-diagrams/

[5]: http://adrai.github.io/flowchart.js/

[6]: /benweet/stackedit

如果觉得《如何用均匀分布随机数生成正态分布随机数》对你有帮助，请点赞、收藏，并留下你的观点哦！