高斯隨機數(shù)生成算法對比研究

張志軍,劉行兵,段新濤

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

高斯隨機數(shù)生成算法對比研究

張志軍,劉行兵,段新濤

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

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

高斯隨機數(shù)；模擬

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

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

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

1 算法描述

1.1 中心極限生成算法

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

1.2 Box-Muller算法

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

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

Box-Muller算法描述：

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

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

1.3 極化判決算法

極化判決算法[2,4]作為Box-Muller方法一個擴充,以不同的方式獲得二維高斯隨機變量.在極化方法中,需要兩個相互獨立的（-1,1）區(qū)間內(nèi)均勻隨機數(shù)產(chǎn)生一高斯隨機向量,但該向量的坐標(biāo)需要判斷：當(dāng)該向量的模大于1時,該均勻隨機數(shù)對將被舍棄；如果其模小于1（發(fā)生的概率為π/4）,將該隨機向量點以一定的方式轉(zhuǎn)換成兩個高斯隨機數(shù)輸出.極化判決算法與Box-Muller算法不同之處在于其不需正/余弦三角函數(shù)計算,但需增加一次除法運算及兩次乘法運算.

極化判決算法描述：

（1）x←V1,y←V2；

（2）d←x2+y2；

（3）如果0＜d＜1,轉(zhuǎn)到（4）,否則,轉(zhuǎn)到（1）；

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

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

2 隨機數(shù)樣本檢驗

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

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

圖1 中心極限生成算法產(chǎn)生的隨機數(shù)對應(yīng)密度函數(shù)對比Fig.1 Probability density function with different n using the central limit algorithm

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

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

皮爾遜擬合優(yōu)度χ2檢驗中[7-8],將樣本的觀測頻率與理論頻率進行比較檢驗.其中樣本觀測頻率是通過在γ個分割子區(qū)域內(nèi)對樣本進行統(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ū)間內(nèi)樣本觀測到的實際頻數(shù),E（·）為子區(qū)間內(nèi)理論頻數(shù),α為顯著性水平,γ為分割的子區(qū)間數(shù)目.表1示出了極化判決算法生成隨機數(shù)的χ2檢驗結(jié)果.

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

表1中χ2obs表示對應(yīng)區(qū)間內(nèi)擬合優(yōu)度χ2檢驗拒絕限,在顯著水平為α下,計算式（1）得統(tǒng)計量若大于該值時,通過擬合檢驗,否則,未通過.從表1可得,極化判決算法在區(qū)間0≤X≤3.0σ和3.0σ＜X≤4.5σ均能通過檢驗；但在4.5σ＜X≤6.0σ區(qū)間內(nèi)只有在顯著水平α=0.1下通過,即極化判決算法生成的高斯隨機數(shù)尾部擬合精度在[4.5σ,6.0σ]之間,和前面的分析是一致的.

表2示出高斯隨機數(shù)3種生成算法所需運算對比情況,相對Box-Muller算法,性能最好的極化判決算法增加了加法、除法和比較運算,這是高尾部精度的代價,硬件實現(xiàn)時增加的運算量在可承受范圍之內(nèi).

表2 3種算法所需運算與資源對比Tab.2 Comparison of computing resources with different algorithms

3 小結(jié)

本文研究了高斯隨機數(shù)3種常用生成算法,包括中心極限生成算法、Box-Muller算法和極化判決算法,對比了各算法生成的高斯隨機數(shù)在尾部區(qū)域內(nèi)的質(zhì)量以及這3種算法所需的運算類型與運算次數(shù).在硬件實現(xiàn)時,極化判決算法不失為最佳的選擇,即以增加一定量運算代價下,換取了高斯隨機數(shù)的高尾部精度.同時需注意的是,極化判決算法由于需要用判決語句,導(dǎo)致產(chǎn)生高斯隨機數(shù)的速度不恒定,在硬件實現(xiàn)時需用先入先出緩沖器解決該問題.

[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] 王松桂,張忠占,程維虎,等.概率論與數(shù)理統(tǒng)計[M].北京：科學(xué)出版社,2011：165-170.

（責(zé)任編輯：盧奇）

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

國家自然科學(xué)基金資助項目（U1204606）；河南省基礎(chǔ)與前沿技術(shù)研究計劃項目（142300410004）；河南省教育廳科學(xué)技術(shù)研究重點研究項目（14B510020）

張志軍（1979-）,男,河南輝縣人,碩士,講師.主要從事無線通信調(diào)制、信道編碼研究.