c++产生均匀分布随机数赋值_不随机的随机数：高斯随机数生成器综述

随机数的使用非常广泛，例如在从统计总体中抽取有代表性的样本时，或者在将实验动物分配到不同的试验组的过程中，或者在进行蒙特卡洛模拟法计算的时候等等。事实上，这些统计技术中使用的随机数均为“伪随机数”，由各种各样的随机数生成器生成。

本文将高斯随机数生成器分为四类：一、直接利用累积分布函数的反函数生成随机数；二、转换变形方法，将均匀分布直接转换为高斯分布；三、拒绝采样，在转换方法中增加了按条件拒绝一些转换值；四、递归迭代方法，对已产生的高斯随机数进行线性变换，从而产生新的随机数。

01 反函数算法

基本原理：假设P(X)是一个累积分布函数，P^(-1)是它的反函数，若U是一个服从(0,1)均匀分布的随机变量，则P^(-1)(U)服从函数P给出的分布。

例如要生成一个服从指数分布的随机变量，我们知道指数分布的概率密度函数为

则累积分布函数为

算法的核心是要找到指数随机变量的值，使得Prob(V≤x)是[0,1]区间上的均匀分布。即给定一个满足均匀分布的变量u之后，要使得

通过求解 x，得到(1)式：

上式的作用就是把输入的均匀随机数 u 转换为指数随机数。

类似地，对正态分布来说，累积分布函数如下图所示，

图1：正态分布累积分布图

在y轴上产生服从(0,1)均匀分布的随机数，水平向右投影到曲线上，再垂直向下投影到x轴，这样在x轴上就得到了正态分布。

然而，正态分布的累积分布函数不能直接用数学形式写出，求解需要采用数值方法或者一些近似估计，例如多项式近似，这使得CDF反函数算法产生随机数的效率降低，关于CDF反函数的近似方法可参考Muller(1958), Gebhardt(1964),Wichura(1988)等。

02 转换变形算法

Box-Muller 方法

Box-Muller 方法(Box & Muller 1958b; Pike 1965)由一对均匀分布随机数得到一对高斯随机数，是一种最早的精确转换方法，并且在很长时间内该方法都是生成正态分布随机数的“标准”算法。算法特点是效率高，计算过程比较简单。

二维标准正态分布的概率密度函数为

在极坐标系下又可写成

由于 p(r) 是关于极坐标轴对称的，我们只需要考虑半径 r , 而角度 θ 可以用 [0,2π] 区间上的均匀分布，由于 p(r) 是对于变量 r^2 的指数分布，其中参数 a=1/2, r 可以通过公式 (1) 生成

其中 u1 是均匀分布随机数。进而生成两个满足标准正态分布且完全独立的随机数

其中 u2 是另一个满足均匀分布的随机数。

Box-Muller 转换算法还有一种极坐标形式，使得运算过程不包含三角函数运算(三角函数运算一般效率较低)，可以参考(Marsaglia and Ray, 1964)。

中心极限定理

中心极限定理：设 X1，X2,...,Xn 为独立同分布的随机变量序列，均值为 μ ，方差为 σ^2，则

具有渐近分布 N(0,1) , 也就是说，当 n 趋近于无穷时，

即n个独立同分布的随机变量之和的分布近似于正态分布，n 越大，近似程度越好。

根据中心极限定理，只要生成 n 个独立的均匀分布，加和就可以近似得到正态分布。这个算法的缺点是随着 n 的增大，收敛到正态分布的累积分布函数的速度会越来越慢。

除了常见的 Box-Muller 方法、中心极限定理方法，转换变形算法中还有用三角形分布分段线性逼近(Kabal，2000)，Monty Python方法(Marsaglia and Tsang，1998)等。另外，(Thomas & Luk， 2006)将大量三角形分布和中心极限定理相结合，提供了一种有效的生成高斯随机数的方法。

03 拒绝采样算法

最简单的拒绝采样(Rejection Sampling)如下图所示，首先在蓝色矩形框中，生成均匀分布的随机点(x, y)，然后拒绝所有未落在正态分布钟形曲线以下的点，剩下点的横坐标 x 就构成一个服从正态分布的随机数序列。

图2：正态分布曲线图

但是显然，直接运用均匀分布的拒绝率很高，会导致算法效率低下，若能有一种方法直接产生在钟形曲线下方面积内均匀分布的随机点，那么就会大大降低拒绝率，提升算法效率，这正是Ziggurat算法(Marsaglia and Tsang 1984a, 2000)的基本思路。

Ziggurat 算法

图3：Ziggurat算法示意图

由于曲线左右对称，故只考虑右半边图形，具体算法如下：

① 我们用 Ri (i=1,...,n) 来表示 n 个矩形，Ri 的右边界为 xi，上边界为 yi. Si 表示一个分割，当 i=n 时，Sn(最下方矩形)=Rn+tail，其他情况 Si=Ri.

② 除了 Rn 以外，其他每个矩形面积相等，设为定值 A. Rn 的面积 = A-tail 的面积。这样保证从任何一个分割中抽样 (x,y) 的概率相等.

③ 当任意选定一个 Ri 在其中抽样 (x,y)，若 x

这里为了方便解释，只用了几个矩形，在程序实现的时候，一般使用128或256个矩形。

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

另外，ziggurat 算法实际上并不局限于生成正态分布随机数，只要一个连续分布的概率密度函数(PDF)图象是一条下降的曲线或是像正态分布这样的钟形曲线，理论上都可以使用 ziggurat 算法生成。

在拒绝采样算法中，除了流行的ziggurat 算法，还有 Polar 方法(Bell 1968; Knop 1969)、Marsaglia-Bray Rejection 方法(Marsaglia & Bray 1964) 、Ratio of Uniforms (Kinderman and Monahan ，1977)、Ahrens-Dieter Table-Free方法(Ahrens and Dieter，1988)、GRAND(Brent，1974)等。

04 递归迭代算法

Wallace随机数生成方法充分利用了高斯分布随机数的线性组合仍然是高斯分布这一性质(例如我们最常用的镜像法，如果u服从N(0, 1)那么-u也服从N(0,1)，就是这一方法的简单特例)，从而完全避免了对基本函数的求值，因此该方法有着最快的速率，但是也存在着相关性问题，具体算法可参考(Wallace，1996)，简要步骤如下：

① 初始化，建立一个有 N = KL 个独立高斯随机数的初始池，并将其归一化(也称之为二阶动量匹配)；

② 线性变换，把每 K 个值看作一个向量 X，令X' = AX，这里 A 是一个K维空间里的旋转矩阵。由于原始的 K 个值是高斯分布的，故通过旋转变换，新得到的 K个值仍为高斯分布的随机数。

经过这一步骤后，我们得到了一个新的高斯分布随机数池，由于 A 是正交矩阵，因此新的随机数池的平方和均值仍为1，这就需要进一步处理。如果x1, x2, ... , xN是高斯分布的独立样本，那么我们期望 x1^2+x2^2+...+xN^2 是一个卡方分布。

③ 校正，从得到的随机数池中取出一个值逐一除剩余值，这些剩余值被除后的 N-1 个数形成最终的高斯随机数池。

Wallace 用原始高斯随机数得到新的高斯随机数的方法，意味着输出样本间不可避免地将会存在一些相关性，但是我们可以通过仔细选择系统参数，将相关性效应降低，使得得到的高斯随机数满足应用中对随机数质量的要求。

05 总结

本文共介绍了4类生成高斯随机数的算法，其中反函数算法由于正态分布累积分布函数求解困难导致算法效率低下；用中心极限定理生成高斯随机数原理简洁，但是收敛速度慢；Box-Muller 方法和 ZIggurat 是目前较流行的两种算法，计算效率高。如果只需生成少量高斯随机数，建议使用 Box-Muller 方法；若要生成大量高斯随机数，则应使用 Ziggurat 算法，这种情形下，该算法效率远高于 Box-Muller 方法。此外，Wallace 方法有着最快的计算速率，但是同时存在相关性问题。

参考文献

[1] B OX , G. E. P. AND M ULLER , M. E. 1958b. A note on the generation of random normal deviates. Annals Math.Stat. 29, 2, 610–611.

[2] MARSAGLIA , G. AND T SANG , W. W. 1984a. A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions. SIAM J. Sci. Statis. Comput. 5, 349–359.

[3] MARSAGLIA , G. AND B RAY , T. A. 1964. A convenient method for generating normal variables. SIAM Rev. 6, 3,260–264.

[4] W ALLACE , C. S. 1996. Fast pseudorandom generators for normal and exponential variates. ACM Trans.Math. Softw. 22, 1, 119–127.

[5] K ABAL , P. 2000. Generating Gaussian pseudo-random deviates. Tech. Rep., Department of Electrical and Computer Engineering, McGill University.

[6] T HOMAS , D. B. AND L UK , W. 2006. Non-uniform random number generation through piecewise linear approximations. In International Conference on Field Programmable Logic and Applications.

[7] A HRENS , J. H. AND D IETER , U. 1988. Efficient table-free sampling methods for the exponential, Cauchy, andnormal distributions. Comm. ACM 31, 11, 1330–1337.

[8] K INDERMAN , A. J. AND M ONAHAN , J. F. 1977. Computer generation of random variables using the ratio of uniform deviates. ACM Trans. Math. Softw. 3, 3, 257–260.

关于镒链

上海镒链金融科技公司专注于风险管理、衍生品估值、金融数据、交易等前台核心技术，是全球FinTech的建设者，中国金融技术的领先者之一。