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

前言

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

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

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

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

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

The Box–Muller transform

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

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

我们希望计算积分$I={\int }_{-\infty }^{\infty }{e}^{-x^2/2}dx$，可以先取它的平方并用极坐标表示，$$I^2={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-(x^2+y^2)/2}dxdy={\int }_0^{2\pi }{\int }_0^{\infty }r{e}^{r^2/2}drd\theta .$$可以看到极角 θ \theta