如何用均匀分布随机数生成正态分布随机数

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

在讨论更高效的算法之前，我们先来看看能不能基于中心极限定理来设计算法。中心极限定理告诉我们，对于一组 i . i . d i.i.d i.i.d的随机数 { x k } ∼ U ( μ , σ 2 ) \{x_k\}\sim U(\mu,\sigma^2) {xk​}∼U(μ,σ2)，有 1 n ∑ i = 1 n x i → 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)。这个算法有两个问题：

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

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

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

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

基本形式： 假设 U 1 U_1 U1​和 U 2 U_2 U2​是区间 ( 0 , 1 ) (0,1) (0,1)上的均匀分布随机数，令 Z 1 = R cos ⁡ ( θ ) = − 2 ln ⁡ U 1 cos ⁡ ( 2 π U 2 ) Z_1=R\cos(\theta)=\sqrt{-2\ln U_1}\cos(2\pi U_2) Z1​=Rcos(θ)=−2lnU1​ ​cos(2πU2​)， Z 2 = R sin ⁡ ( θ ) = − 2 ln ⁡ U 1 sin ⁡ ( 2 π U 2 ) Z_2=R\sin(\theta)=\sqrt{-2\ln U_1}\sin(2\pi U_2) Z2​=Rsin(θ)=−2lnU1​ ​sin(2πU2​)， 则 Z 1 Z_1 Z1​和 Z 2 Z_2 Z2​是两个独立的标准正态分布随机变量。

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

算法局限性： 尾部截断(Tails truncation) 由于计算机有数值精度限制，其表示的数字不可能无线接近0。比如，使用32bit精度，那么最接近0的数字是 2 − 32 2^{-32} 2−32。因此，当 u 1 u_1 u1​和 u 2 u_2 u2​都取 2 − 32 2^{-32} 2−32时， x x x的最大值是 − 2 ln ⁡ ( 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.6 6.6 6.6个标准差，这意味着会有 2 × 1 0 − 11 2×10^{-11} 2×10−11的比例(尾部)会因为计算机精度问题无法采样到，这在大规模生成高斯随机数和研究rare events的时候要特别注意。

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

数学解释： 假设 F ( x ) = P r [ X ; x ] F(x)=Pr[X;x] F(x)=Pr[X