高斯随机数生成算法对比研究

张志军,刘行兵,段新涛

（河南师范大学 计算机与信息工程学院,河南 新乡 543007）

高斯随机数生成算法对比研究

张志军,刘行兵,段新涛

（河南师范大学 计算机与信息工程学院,河南 新乡 543007）

针对应用学科仿真实验中所需高斯随机数质量要求日趋严格的问题,对比研究了高斯随机数常用的中心极限生成算法、Box-Muller算法和极化判决算法,研究了各算法在尾部区域内生成高斯随机数的质量.仿真结果表明：极化判决算法为最佳选择,在增加一定量运算代价下,换取高斯随机数的高尾部精度.需注意极化判决算法中由于用到判决语句,导致高斯随机数的产生速度不恒定,硬件实现时需用先入先出缓冲器解决该问题.

高斯随机数；模拟

需要高斯分布随机数模拟实验的领域包括通信系统、金融建模等.在高斯随机数生成算法中,通常将服从均匀分布的随机数经算术运算转换得到,在可接受精度要求下,产生出满足特定应用环境下质量要求的高斯随机数[1-6].

高斯随机数生成算法研究的目标是为特定的应用提供足够精确的、统计独立的、取自服从理想高斯分布的随机变量的样本点.根据构造方法可将高斯随机数产生算法划分为“精确”法和“逼近”法.在“理想”环境下,精确法将会产生精准的高斯随机数,例如Box-Muller法[2-3]和极化判决算法[2,4].与精确法相比,不管计算精度多么高,逼近法生成的随机数是近似地服从高斯分布.逼近法的一个例子就是中心极限生成算法[1-2].使用该方法产生高斯随机数时,只有当无穷多个均匀随机数经线性组合后生成的高斯随机数才是精确的,当然也是不现实的.通常用该方法生成的高斯随机数都是近似的.

本文研究上述3种高斯随机数生成算法,对比这3种算法所需计算类型与计算次数,重点研究各算法在尾部区域内生成的高斯随机数质量,为硬件实现高斯随机数生成器提供参考.

1 算法描述

1.1 中心极限生成算法

多个相互独立的均匀随机数之和生成的随机数,其概率密度函数可通过每个均匀随机数密度函数逐个卷积得到.然而,中心极限理论提供了上述问题的另一种求法.由中心极限定理（Central Limit Theorem,CLT）可知[1],n个相互独立的、在（-0.5,0.5）区间上服从均匀分布的随机数之和的密度函数逼近于均值为零、方差为n/12的高斯分布概率密度函数,当n趋于无穷多个时,逼近的误差趋近于零.该方法原理比较简单,主要缺点是收敛速度过慢.

1.2 Box-Muller算法

Box-Muller算法是最早的精确转换方法之一,该方法将一对均匀随机数转换生成一对高斯随机数. Box-Muller算法基于下述事实：两个相互独立的、均值均为零、方差相同的高斯随机数联合二维分布是呈放射状对称的.

Box-Muller算法输出的高斯随机数可以认为是二维平面上某随机点坐标,该随机点的振幅由（0,1）区间上服从均匀分布的随机数转换而来,其相位为（0,1）区间上一均匀随机数乘2π得到.将该随机点映射到笛卡尔坐标轴上,其对应的坐标点就是服从高斯分布的随机数.文献[6]基于优化的函数硬件实现了固定精度的Box-Muller高斯随机数生成器,该硬件中正弦、余弦的计算可以在一步操作中同时实现.

Box-Muller算法描述：

（1）a2←-2ln（U1）,b←2πU2；

（2）返回一对相互独立的高斯随机数（a sinb,a cosb）.

1.3 极化判决算法

极化判决算法[2,4]作为Box-Muller方法一个扩充,以不同的方式获得二维高斯随机变量.在极化方法中,需要两个相互独立的（-1,1）区间内均匀随机数产生一高斯随机向量,但该向量的坐标需要判断：当该向量的模大于1时,该均匀随机数对将被舍弃；如果其模小于1（发生的概率为π/4）,将该随机向量点以一定的方式转换成两个高斯随机数输出.极化判决算法与Box-Muller算法不同之处在于其不需正/余弦三角函数计算,但需增加一次除法运算及两次乘法运算.

极化判决算法描述：

（1）x←V1,y←V2；

（2）d←x2+y2；

（3）如果0＜d＜1,转到（4）,否则,转到（1）；

（4）f2←-2（lnd）/d；

（5）返回（f×x,f×y）.

2 随机数样本检验

用概率密度函数图和拟合优度χ2检验对比上述3种高斯随机数生成算法的精准度.

图1示出由中心极限生成算法生成的随机数所对应的归一化概率密度函数曲线随变量个数n变化情况.实验中所需的n个相互独立的随机变量均为取值于（-0.5,0.5）区间服从均匀分布的随机数,每个随机数取109个样本.图1上半部分示出该随机数线性尺度下的密度函数曲线,当变量个数n大于2时,生成的概率密度函数曲线看似“非常”地接近于“理想”高斯分布密度函数,这是由于随机变量大于2σ时,随机事件发生的概率非常低以至于在线性尺度下区分不开.图1下半部分示出了该密度函数对数尺度下的曲线.很显然,当n等于12时所对应概率密度函数,在变量值大于3σ时密度函数曲线已“脱离”于理想的高斯分布密度函数曲线.

图1 中心极限生成算法产生的随机数对应密度函数对比Fig.1 Probability density function with different n using the central limit algorithm

当高斯分布密度函数尾部精度在6σ以上时,其随机事件发生的概率小于10-9,为此需产生随机样本数应至少为1011.由于线性尺度下该事件概率的差别几乎不可观测到,图2仅示出3种算法所生成的随机数对数尺度下的密度函数曲线.从图2可以看出,中心极限算法生成随机数质量最差,当大于3σ时,生成的密度函数曲线已“脱离”于理想的高斯分布密度函数曲线；Box-Muller算法“脱离”于理想的高斯分布密度函数曲线发生在约4.6σ左右；在这3种算法中,极化判决算法生成高斯随机数的质量最好,该算法生成高斯随机数的密度函数与理想的高斯分布密度函数开始“分离”发生在约5.4σ左右.

图2 3种算法生成的密度函数曲线对比Fig.2 Probability density function with different algorithms

皮尔逊拟合优度χ2检验中[7-8],将样本的观测频率与理论频率进行比较检验.其中样本观测频率是通过在γ个分割子区域内对样本进行统计.检验时,对不同的样本区间采用不同的分割方法：0≤X≤3.0σ区间分割100个子区间,3.0σ＜X≤4.5σ区间分割50个子区间,4.5σ＜X≤6.0σ区间分割成30个子区间.检验的χ2统计量用下式计算

式（1）中O（·）为子区间内样本观测到的实际频数,E（·）为子区间内理论频数,α为显著性水平,γ为分割的子区间数目.表1示出了极化判决算法生成随机数的χ2检验结果.

表1 极化判决算法的χ2检验结果Tab.1 Chi-Square test results for Polar-Rejection algorithm

表1中χ2obs表示对应区间内拟合优度χ2检验拒绝限,在显著水平为α下,计算式（1）得统计量若大于该值时,通过拟合检验,否则,未通过.从表1可得,极化判决算法在区间0≤X≤3.0σ和3.0σ＜X≤4.5σ均能通过检验；但在4.5σ＜X≤6.0σ区间内只有在显著水平α=0.1下通过,即极化判决算法生成的高斯随机数尾部拟合精度在[4.5σ,6.0σ]之间,和前面的分析是一致的.

表2示出高斯随机数3种生成算法所需运算对比情况,相对Box-Muller算法,性能最好的极化判决算法增加了加法、除法和比较运算,这是高尾部精度的代价,硬件实现时增加的运算量在可承受范围之内.

表2 3种算法所需运算与资源对比Tab.2 Comparison of computing resources with different algorithms

3 小结

本文研究了高斯随机数3种常用生成算法,包括中心极限生成算法、Box-Muller算法和极化判决算法,对比了各算法生成的高斯随机数在尾部区域内的质量以及这3种算法所需的运算类型与运算次数.在硬件实现时,极化判决算法不失为最佳的选择,即以增加一定量运算代价下,换取了高斯随机数的高尾部精度.同时需注意的是,极化判决算法由于需要用判决语句,导致产生高斯随机数的速度不恒定,在硬件实现时需用先入先出缓冲器解决该问题.

[1] Fischer H.A history of the central limit theorem：Fromclassical to modern probability theory[M].NewYork：Springer press,2010.

[2] Thomas D B,Luk W,Leong P H W,et al.Gaussian randomnumber generators[J].ACMComputing Surveys（CSUR）,2007,39（4）：11.

[3] Box G,Muller ME.A note on the Generation of randomnormal deviates[J].Annals of Math.Stat.,1958,29（2）：610-611.

[4] Muller ME.A comparison of methods for generating normal deviates on digital computers[J].Journal of the ACM（JACM）,1959,6（3）：376-383.

[5] Lee D U,Villasenor J D,Luk W,et al.A hardware Gaussian noise generator using the Box-Muller method and its error analysis[J]. IEEE Transactions on Computers,2006,55（6）：659-671.

[6] Alimohammad A,Fard S F,Cockburn B F,et al.A compact and accurate Gaussian variate generator[J].IEEE Transactions on VLSI Systems,2008,16（5）：517-527.

[7] D'Agostino R B.Goodness-of-fit-techniques[M].NewYork：CRC press,1986.

[8] 王松桂,张忠占,程维虎,等.概率论与数理统计[M].北京：科学出版社,2011：165-170.

（责任编辑：卢奇）

Algorithmcomparison on Gaussian randomnumber generators

Zhang Zhijun,Liu Xingbing,Duan Xintao
（College of Computer and Information Engineering,Henan Normal University,Xinxiang 453007,China）

Aiming at the quality of Gaussian randomnumbers in simulation experiments for various applied sciences,comparative study was made on the generate algorithms used to generating Gaussian randomnumber, focusing on the quality of each algorithmto generate Gaussian randomnumbers in the tail region.Simulation results showthat,namely to increase the operational costs,Polar-Rejection algorithmmay well be the best option with hightail accuracy.It should be noted that,due to using if-else statements in Polar-Rejection algorithm,resulting the rate of Gaussian randomnumbers is not constant.It should be using first-in first-out buffer to solve the problemin hardware implementation.

Gaussian randomnumber；simulation

TP301.6

A

1008-7516（2014）03-0056-04

10.3969/j.issn.1008-7516.2014.03.013

2014-03-24

国家自然科学基金资助项目（U1204606）；河南省基础与前沿技术研究计划项目（142300410004）；河南省教育厅科学技术研究重点研究项目（14B510020）

张志军（1979-）,男,河南辉县人,硕士,讲师.主要从事无线通信调制、信道编码研究.