張志軍,劉行兵,段新濤

（河南師范大學 計算機與信息工程學院,河南 新鄉(xiāng) 543007）

高斯隨機數生成算法對比研究

張志軍,劉行兵,段新濤

（河南師范大學 計算機與信息工程學院,河南 新鄉(xiāng) 543007）

針對應用學科仿真實驗中所需高斯隨機數質量要求日趨嚴格的問題,對比研究了高斯隨機數常用的中心極限生成算法、Box-Muller算法和極化判決算法,研究了各算法在尾部區(qū)域內生成高斯隨機數的質量.仿真結果表明：極化判決算法為最佳選擇,在增加一定量運算代價下,換取高斯隨機數的高尾部精度.需注意極化判決算法中由于用到判決語句,導致高斯隨機數的產生速度不恒定,硬件實現時需用先入先出緩沖器解決該問題.

高斯隨機數；模擬

需要高斯分布隨機數模擬實驗的領域包括通信系統(tǒng)、金融建模等.在高斯隨機數生成算法中,通常將服從均勻分布的隨機數經算術運算轉換得到,在可接受精度要求下,產生出滿足特定應用環(huán)境下質量要求的高斯隨機數[1-6].

高斯隨機數生成算法研究的目標是為特定的應用提供足夠精確的、統(tǒng)計獨立的、取自服從理想高斯分布的隨機變量的樣本點.根據構造方法可將高斯隨機數產生算法劃分為“精確”法和“逼近”法.在“理想”環(huán)境下,精確法將會產生精準的高斯隨機數,例如Box-Muller法[2-3]和極化判決算法[2,4].與精確法相比,不管計算精度多么高,逼近法生成的隨機數是近似地服從高斯分布.逼近法的一個例子就是中心極限生成算法[1-2].使用該方法產生高斯隨機數時,只有當無窮多個均勻隨機數經線性組合后生成的高斯隨機數才是精確的,當然也是不現實的.通常用該方法生成的高斯隨機數都是近似的.

本文研究上述3種高斯隨機數生成算法,對比這3種算法所需計算類型與計算次數,重點研究各算法在尾部區(qū)域內生成的高斯隨機數質量,為硬件實現高斯隨機數生成器提供參考.

1 算法描述

1.1 中心極限生成算法

多個相互獨立的均勻隨機數之和生成的隨機數,其概率密度函數可通過每個均勻隨機數密度函數逐個卷積得到.然而,中心極限理論提供了上述問題的另一種求法.由中心極限定理（Central Limit Theorem,CLT）可知[1],n個相互獨立的、在（-0.5,0.5）區(qū)間上服從均勻分布的隨機數之和的密度函數逼近于均值為零、方差為n/12的高斯分布概率密度函數,當n趨于無窮多個時,逼近的誤差趨近于零.該方法原理比較簡單,主要缺點是收斂速度過慢.

1.2 Box-Muller算法

Box-Muller算法是最早的精確轉換方法之一,該方法將一對均勻隨機數轉換生成一對高斯隨機數. Box-Muller算法基于下述事實：兩個相互獨立的、均值均為零、方差相同的高斯隨機數聯合二維分布是呈放射狀對稱的.

Box-Muller算法輸出的高斯隨機數可以認為是二維平面上某隨機點坐標,該隨機點的振幅由（0,1）區(qū)間上服從均勻分布的隨機數轉換而來,其相位為（0,1）區(qū)間上一均勻隨機數乘2π得到.將該隨機點映射到笛卡爾坐標軸上,其對應的坐標點就是服從高斯分布的隨機數.文獻[6]基于優(yōu)化的函數硬件實現了固定精度的Box-Muller高斯隨機數生成器,該硬件中正弦、余弦的計算可以在一步操作中同時實現.

Box-Muller算法描述：

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

（2）返回一對相互獨立的高斯隨機數（a sinb,a cosb）.

1.3 極化判決算法

極化判決算法[2,4]作為Box-Muller方法一個擴充,以不同的方式獲得二維高斯隨機變量.在極化方法中,需要兩個相互獨立的（-1,1）區(qū)間內均勻隨機數產生一高斯隨機向量,但該向量的坐標需要判斷：當該向量的模大于1時,該均勻隨機數對將被舍棄；如果其模小于1（發(fā)生的概率為π/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 隨機數樣本檢驗

用概率密度函數圖和擬合優(yōu)度χ2檢驗對比上述3種高斯隨機數生成算法的精準度.

圖1示出由中心極限生成算法生成的隨機數所對應的歸一化概率密度函數曲線隨變量個數n變化情況.實驗中所需的n個相互獨立的隨機變量均為取值于（-0.5,0.5）區(qū)間服從均勻分布的隨機數,每個隨機數取109個樣本.圖1上半部分示出該隨機數線性尺度下的密度函數曲線,當變量個數n大于2時,生成的概率密度函數曲線看似“非?！钡亟咏凇袄硐搿备咚狗植济芏群瘮?這是由于隨機變量大于2σ時,隨機事件發(fā)生的概率非常低以至于在線性尺度下區(qū)分不開.圖1下半部分示出了該密度函數對數尺度下的曲線.很顯然,當n等于12時所對應概率密度函數,在變量值大于3σ時密度函數曲線已“脫離”于理想的高斯分布密度函數曲線.

圖1 中心極限生成算法產生的隨機數對應密度函數對比Fig.1 Probability density function with different n using the central limit algorithm

當高斯分布密度函數尾部精度在6σ以上時,其隨機事件發(fā)生的概率小于10-9,為此需產生隨機樣本數應至少為1011.由于線性尺度下該事件概率的差別幾乎不可觀測到,圖2僅示出3種算法所生成的隨機數對數尺度下的密度函數曲線.從圖2可以看出,中心極限算法生成隨機數質量最差,當大于3σ時,生成的密度函數曲線已“脫離”于理想的高斯分布密度函數曲線；Box-Muller算法“脫離”于理想的高斯分布密度函數曲線發(fā)生在約4.6σ左右；在這3種算法中,極化判決算法生成高斯隨機數的質量最好,該算法生成高斯隨機數的密度函數與理想的高斯分布密度函數開始“分離”發(fā)生在約5.4σ左右.

圖2 3種算法生成的密度函數曲線對比Fig.2 Probability density function with different algorithms

皮爾遜擬合優(yōu)度χ2檢驗中[7-8],將樣本的觀測頻率與理論頻率進行比較檢驗.其中樣本觀測頻率是通過在γ個分割子區(qū)域內對樣本進行統(tǒng)計.檢驗時,對不同的樣本區(qū)間采用不同的分割方法：0≤X≤3.0σ區(qū)間分割100個子區(qū)間,3.0σ＜X≤4.5σ區(qū)間分割50個子區(qū)間,4.5σ＜X≤6.0σ區(qū)間分割成30個子區(qū)間.檢驗的χ2統(tǒng)計量用下式計算

式（1）中O（·）為子區(qū)間內樣本觀測到的實際頻數,E（·）為子區(qū)間內理論頻數,α為顯著性水平,γ為分割的子區(qū)間數目.表1示出了極化判決算法生成隨機數的χ2檢驗結果.

表1 極化判決算法的χ2檢驗結果Tab.1 Chi-Square test results for Polar-Rejection algorithm

表1中χ2obs表示對應區(qū)間內擬合優(yōu)度χ2檢驗拒絕限,在顯著水平為α下,計算式（1）得統(tǒng)計量若大于該值時,通過擬合檢驗,否則,未通過.從表1可得,極化判決算法在區(qū)間0≤X≤3.0σ和3.0σ＜X≤4.5σ均能通過檢驗；但在4.5σ＜X≤6.0σ區(qū)間內只有在顯著水平α=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算法和極化判決算法,對比了各算法生成的高斯隨機數在尾部區(qū)域內的質量以及這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] 王松桂,張忠占,程維虎,等.概率論與數理統(tǒng)計[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-）,男,河南輝縣人,碩士,講師.主要從事無線通信調制、信道編碼研究.