如何产生正态分布的随机数？

Milo Yip ，游戏程序员、《游戏引擎架构》译者

我为了这个问题做了个开源项目 miloyip/normaldist-benchmark · GitHub。

这个项目比较不同正态随机数生成算法的性能，它测试的函数原型是：

void normaldistf(float* data, size_t n);
void normaldist(double* data, size_t n);

它测试了以下实现的性能，并验证其正确性（mean, SD, Skewness, Kurtosis）。

boxmuller: Box-Muller Transform。
cpp11normal: C++11 的 std::normal_distribution。（在 clang/LLVM 中的实现是 Box-Muller transform。）
clt：利用中心极限定理，把m个均匀分布随机数相加。
inverse: 利用正态分布 CDF 的反函数做 Inverse transform sampling。
marsagliapolar: Marsaglia polar method。
ziggurat: Ziggurat algorithm，n=128。
null: 作为比较基准，生成均匀分布的随机数。

一些算法包含 SSE2 和 AVX 指令集的版本。

注意 clt 算法的 abs(kurtosis) >= 0.01，未能通过正确性验证。其他除了 null 外，都通过正确性验证。

先看看 float 版本的结果（iMac Corei5-3330S@2.70GHz clang 6.1 64-bit）：

http://rawgit.com/miloyip/normaldist-benchmark/master/result/Corei5-3330S@2.70GHz_mac64_clang6.1.html

可能很多人以为 Zigguart 算法是最快的，但查表操作、不定数量的分支都不利于 SIMD 并行化。而 Marsaglia polar method 也要做 rejection sampling，所以也会造成分支。Marsaglia polar method 原意是为了改善 Box-Muller 的性能的，但在此测试中反而较慢。Inverse 虽然简单，但该反函数也需要分支处理中央和尾的部分，形成分支。

在加入 SSE2 和 AVX 加速后（其 PRNG 为简单 SIMD 并行的 LCG），Box-Muller 有很明显优势，SSE2 已是 ziggurat 的 1.79x，AVX 是 2.99x。AVX 实现中使用了SSE2/SSE4.1的指令分拆来做一些操作，要在 AVX2 中才能完全实现 8-way 计算。

相对于维基上 Box-Muller 的例子代码：

double generateGaussianNoise(double mu, double sigma)
{
	const double epsilon = std::numeric_limits<double>::min();
	const double two_pi = 2.0*3.14159265358979323846;

	static double z0, z1;
	static bool generate;
	generate = !generate;

	if (!generate)
	   return z1 * sigma + mu;

	double u1, u2;
	do
	{
	   u1 = rand() * (1.0 / RAND_MAX);
	   u2 = rand() * (1.0 / RAND_MAX);
	}
	while ( u1 <= epsilon );

	z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
	z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
	return z0 * sigma + mu;
}

我实现的Box-Muller 做了一些优化：

template <class T>
void boxmuller(T* data, size_t count) {
    assert(count % 2 == 0);
    static const T twopi = T(2.0 * 3.14159265358979323846);

    LCG<T> r;
    for (size_t i = 0; i < count; i += 2) {
        T u1 = 1.0f - r(); // [0, 1) -> (0, 1]
        T u2 = r();
        T radius = std::sqrt(-2 * std::log(u1));
        T theta = twopi * u2;
        data[i    ] = radius * std::cos(theta);
        data[i + 1] = radius * std::sin(theta);
    }
}

直接写入两个结果（这个与 API 设计相关）
把 [0, 1) 的随机数转换成 (0, 1]，避免了 log(0) 的出现，不需要做 rejection。
在 SSE2/AVX 版本中，利用 sincos() 函数同时计算 sin 和 cos。

如果目标平台有 SIMD 指令集，建议使用 Box-Muller，否则可考虑 Zigguart。

如果有其他改进或算法，欢迎 pull request。

编辑于 2015-07-09 8 条评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

何史提，从一到无穷大

小小、知乎用户、知乎用户等人赞同

谢邀。

先假设 $Z \sim \mathcal{N}(0,1)$ 。

用Box-Muller方法，随机抽出两个从

均匀分布的数字

和

。然后
$z_1 = \sqrt{-2 \log u} \cos 2\pi v$
$z_2 = \sqrt{-2 \log u} \sin 2\pi v$
那

和

都是正态分布的。
证明可用极坐标，请参考教科书中的Box-Muller方法。

另外使用反函数，先随机抽出一个从

均匀分布的数字

，然后
$z=\sqrt{2} \text{erf}^{-1}(2u-1)$
那

是正态分布的。

Python实现：

import numpy as np
from scipy.special import erfinv

def boxmullersampling(mu=0, sigma=1, size=1):
    u = np.random.uniform(size=size)
    v = np.random.uniform(size=size)
    z = np.sqrt(-2*np.log(u))*np.cos(2*np.pi*v)
    return mu+z*sigma
    
def inverfsampling(mu=0, sigma=1, size=1):
    z = np.sqrt(2)*erfinv(2*np.random.uniform(size=size)-1)
    return mu+z*sigma

发布于 2015-05-05 3 条评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

孙晓博，重度网瘾需电击，晚期懒癌无药医

知乎用户、Porphyah、Xenophen Li 等人赞同

最简单的：rejection sampling，思路很简单，也很容易实现，但效率较差
较复杂的：inverse CDF，直接利用累积分布函数（CDF）的反函数生成随机数，但计算中牵扯到比较复杂的误差函数erf（非初等函数）
更好的：Box-Muller算法，在很长时间内都是生成正态分布随机数的"标准"算法。Box-Muller算法的特点是效率高，并且计算过程比较简单（只用到了初等函数）。参见：Box-Muller transform
目前最好的（相较于其它实用算法）：ziggurat算法，效率很高，很多现代的编程语言都使用了这一算法。ziggurat并不是人名，其含义是“金字形神塔”，不是埃及那个金字塔，而是古代苏美尔人建造的类金字塔结构的神坛：神坛由多层平台构成，每层平台都呈矩形、卵形或正方形，且自下而上面积逐渐减小。ziggurat算法实际上是一种改进的、包含查表操作的rejection sampling。参见：Ziggurat algorithm

--------------------------------------------------------------------------------------------------------------------------------------

楼上知友 @何史提的答案中已经涵盖了 inverse CDF和 Box-Muller算法，这里只简要介绍一下 ziggurat算法：

先考虑最简单的 rejection sampling，其方法如下图所示，首先生成在蓝色矩形框中均匀分布的随机点

，然后舍弃（即拒绝，reject）所有未落在正态分布钟形曲线以下的点，剩余点的横坐标

就构成一个服从相应正态分布的随机数序列。rejection sampling的效率之所以低下，就是由于很多生成的随机点因为落在钟形曲线以上而被白白舍弃。

如果有一种方法能够直接产生在钟形曲线下方面积内均匀分布的随机点，那么就不存在舍弃随机点的情况，也就能够避免上面所说的“浪费”问题，这正是ziggurat算法的基本思路。具体来说，ziggurat算法的内容如下：

% 以下内容均遵从MATLAB语言的规定，数组及向量的下标索引从1开始编号

正态分布的钟形曲线左右对称，这里只取右半边考虑，如下图所示。首先将钟形曲线以下的面积等分为

份，包括

个水平矩形和最下方一直延伸至无限大的尾巴，即 个矩形及最下方尾巴的面积都相等。这些层叠的矩形看起来像古代苏美尔人建造的金字形神塔，这就是ziggurat这个名称的由来。图中取

，实践中

可能达到

以上。设各水平矩形右边界的横坐标（即各圆圈标记的横坐标）为 $z_k\ (k=1,2,\dots,n)$ ，其中

。当

确定后，

也就确定了，可以通过数值方法求解超越方程得到

的值并制表记录。这里称长度方向上与上一层矩形重叠的部分为当前矩形的 核心区，下图中各核心区的右边界就是虚线边，显见最上面的矩形没有核心区。
图片出自 Numerical Computing with MATLAB一书的第九章 Random Numbers。顺便提一下，此书作者就是MATLAB的创始者、MathWorks公司的主席和首席数学家 Cleve B. Moler。

定义 $\sigma_k\equiv z_{k-1}/z_k\ (k=2,3,\dots,n)$ ，即相邻两层矩形水平方向长度的比值，规定 $\sigma_1\equiv0$ 。显然，依照所要求的参数（ $\mu,\,\sigma$ ）生成正态分布时，所有的

、 $\sigma_k$ 都只需要在初始化时计算一次并制表记录即可，之后不需要再重复计算，除非所要求生成的正态分布的参数（ $\mu,\,\sigma$ ）发生改变。

初始化结束后就可以开始生成正态分布随机数了：首先生成一个

区间内均匀分布的随机整数和一个

区间内均匀分布的随机 浮点数。然后检查 $|u|<\sigma_j$ 是否成立，若成立则

就在第

段核心区内，显然也就在钟形曲线下方，故可将

作为正态分布的一个采样点返回。相应的伪代码为

j = randint(1, n);  % [1,n]区间内均匀分布的随机整数
u = 2 * rand() - 1;  % （-1,1）区间内均匀分布的随机浮点数
if abs(u) < sigma[j]
    return u * z[j];  % u * z[j]在第j段核心区内
end

由此可见，使用ziggurat算法时只需要生成一个随机整数和一个随机浮点数，做一次条件判断和一次乘法运算，外加两次查表即可得到一个服从正态分布的随机数，整个过程中没有开根号、求对数、算三角函数等复杂运算。

在上面伪代码中，那个if条件判断在绝大多数情况下都为true，剩下为false的就是被reject的情况，具体又包含三种可能：(1).

，最上面的矩形没有核心区；(2). $j=2\sim(n-1)$ ，

恰好落在核心区之外（即右侧包含曲线的小矩形中）；(3).

，

落在底部核心区之外的尾巴中。在这三种情况下就需要依据Box-Muller算法利用已经生成的均匀分布随机数进行额外运算以生成正态分布随机数。容易看出，

越大，出现这三种需要额外运算的情况的可能性就越低。 Numerical Computing with MATLAB一书中给出的数据是当

时，需要额外运算的可能性就小于3%，所以这部分额外运算对ziggurat算法的总体效率影响很小。

ziggurat算法适用于连续生成大量服从同一正态分布的随机数的情形，在这种情况下ziggurat算法的效率远高于Box-Muller算法。反之，如果只是需要少量正态分布随机数，那么考虑到计算

、 $\sigma_k$ 的开销，ziggurat算法的总体表现未必优于Box-Muller算法。

ziggurat算法的基本思想最早由美国数学家兼计算机科学家George Marsaglia于20世纪60年代提出，在之后的岁月中该算法不断得到改进和发展，目前常用的版本大多源自于George Marsaglia和香港数学家曾衛寰（Wai Wan Tsang）在2000年发表的论文 The Ziggurat Method for Generating Random Variables。

最后需要说明两点：

实际使用的ziggurat算法会比上面介绍的版本更加复杂。如果要在实际项目（特别是与加密有关的项目）中使用正态分布随机数，强烈建议直接使用标准库或是经广泛检验的第三方库所提供的函数，而不是照着上面的介绍自己随便写一个凑合着用。
ziggurat算法实际上并不局限于生成正态分布随机数，只要一个连续分布的概率密度函数（PDF）图象是一条下降的曲线或是像正态分布这样的钟形曲线，理论上都可以使用ziggurat算法生成。

编辑于 2015-07-04 10 条评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

知乎用户，摸鱼失败身败名裂

落菌、知乎用户赞同

如何生成随机数（上） « Free Mind

发布于 2015-07-02 添加评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

慧航，公众号:ecopaper，Keep calm and DO your…

Light Spell、毕成、知乎用户赞同

刚好我的专栏里面有一个现成的code： Julia call C - 码农经济学 - 知乎专栏

编辑于 2015-05-11 添加评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

九铭书记，知道一点微末的数学知识，一直想还原有趣…

李修赞同

最简单粗暴的方法：下载r，写一行r代码：y<--dnorm（x，0，1）。这是标准正态分布。把0和1可以替换成想要的mean和standard deviation

发布于 2015-07-04 2 条评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

匿名用户

知乎用户赞同

在获得均匀分布随机数的基础上, 进行 Box-Muller 变换, 一次计算可获得两个正态分布的随机数.
Box-Muller transform: https://en.wikipedia.org/wiki/Box-Muller_transform

下面是 C# 代码的实现:

/// <summary>
/// 产生正态分布的随机数
/// </summary>
/// <param name="averageValue"> 正态分布的平均值, 即 N(μ, σ^2) 中的 μ </param>
/// <param name="standardDeviation"> 正态分布的标准差, 即 N(μ, σ^2) 中的 σ </param>
/// <returns> 返回正态分布的随机数. 理论值域是 μ±∞; 落在 μ±σ, μ±2σ, μ±3σ 的概率依次为 68.26%, 95.44%, 99.74% </returns>
public float Normal(float averageValue, float standardDeviation)
{
	return averageValue + standardDeviation * (float)
		(
			Math.Sqrt(-2.0 * Math.Log(1.0 - NextDouble()))
			* Math.Sin((Math.PI + Math.PI) * NextDouble())
		);
}

其中 NextDouble 是产生一个 [0, 1) 区间的均匀分布随机数. 这个方法丢弃了另一个结果, 并且一次调用导致随机数种子变化两次.