正態(tài)分布隨機數(shù)生成算法

1、概率論與數(shù)理統(tǒng)計課程設(shè)計題目：正態(tài)分布隨機數(shù)生成算法要編程得到服從均勻分布的偽隨機數(shù)是容易的。C語言、Java語言等都提供了相應(yīng)的函數(shù)。但是要想生成服從正態(tài)分布的隨機數(shù)就沒那么容易了。得到服從正態(tài)分布的隨機數(shù)的基本思想是先得到服從均勻分布的隨機數(shù)，再將服從均勻分布的隨機數(shù)轉(zhuǎn)變?yōu)榉恼龖B(tài)分布。接下來就先分析三個從均勻分布到正態(tài)分布轉(zhuǎn)變的方 法。然后編程實現(xiàn)其中的兩個方法并對程序?qū)崿F(xiàn)運作的效果進行統(tǒng)計分析。1、方法分析(1)利用分布函數(shù)的反函數(shù)若要得到分布函數(shù)為F(x)的隨機變量Y?？闪钛?=f -(u)，其中u是服從均勻分布的隨機變量，有p (丫 乞 y)1二 p(u < F -(y)二

2、F (y)因而，對于任意的分布函數(shù)，只要求出它的反函數(shù)，就可以由服從均勻分布的隨機變量實例來生成服從該分布函數(shù)的隨機變量實例。現(xiàn)在來看正態(tài)分布的分布函數(shù)，對于XN(,；2)，其分布函數(shù)為：1x顯然，要想求其反函數(shù)是相當困難的，同時要想編程實現(xiàn)也很復(fù)雜。可見，用此種方法來生成服從正態(tài)分布的隨機變量實例并不可取。(2)利用中心極限定理第二種方法利用林德伯格一萊維(Lindeberg Levi)中心極限定理：如果隨機變量序列X- X2,，獨立同分布，并且具有有限的數(shù)學(xué)期望和方差2E Xi =，D Xi =：.:0(i=1,2,)，則對一切R 有因此，對于服從均勻分布的隨機變量 Xi，只要n充分大，隨

3、機變量- Xi - v J 就服從N 0,1。我們將實現(xiàn)這一方法。(3)使用Box Muller方法 2x先證明e 2 dx 二 2 二:x2 oO令 I2 dx，則2 _xI2oO e 2QO dx e2y2dy 二e2 dxdy令 x = r cos y 二 r sin v，則有-石=2 二 o e rdr = 2 二。接下來再來得出 Box Muller方法:設(shè)X ,Y為一對相互獨立的服從正態(tài)分布的隨機變量。則有概率密度函數(shù)x2：；y2令x = R cos h y二R sin二，其中v 0, 2，貝U R有分布函數(shù):P(R _r)=r o3122_ue 2 udud v2 _r2 -：e

4、 2 udu0丄2令Fr r =1 e 2 RZpln 1 Z如果X服從均勻分布，則R的分布函數(shù)即為Fr r。最后，可以用代替1 -Z，令二為2二U 2，其中6 U 0,1 , U2 0,1，得:X = R cos v -、. -2 ln U t cos 2c. U 2 ，丫 = R sin v -2 ln U t sin 2c. U 2從而，只需要有兩個服從均勻分布的隨機變量U-U2，就能通過公式X = R cos =2 l n U 勺 cos 2 二U 2來得到一個服從正態(tài)分布的隨機變量X。用Box Muller方法來生成服從正態(tài)分布的隨機數(shù)是十分快捷方便的。我們也將實現(xiàn)這一方法。2、實現(xiàn)

5、與分析（1）利用中心極限定理方法的實現(xiàn)與分析利用中心極限定理來生成隨機數(shù)的函數(shù)（C+語言）編寫如下：const int N = 200;double getRand（）double s = 0;for ( int i = 0; i != N; +i)s += double (rand() % 1000) / 1000;return s;函數(shù)生成的隨機數(shù)是N個0,1間服從均勻分布的隨機數(shù)的和。這里N為200。從而理論上產(chǎn)生的隨機數(shù)應(yīng)近似服從N(n、n；2),其中n為N,即200,為0.5,；二2為1/12。程序生成了 200個隨機數(shù)，并求出樣本均值與樣本方差，也即與二2的最大似然估計：/生成隨機

6、數(shù)并存儲double sum, store200, xi, su = 0, sb = 0, ssb = 0;int cnt = 0;sum = 0;for ( int i = 0; i != 200; +i)xi = getRand();sum += xi;storei = xi;/得到樣本均勻與樣本方差su = sum / 200;for ( int i = 0; i != 200; +i)sb += (storei - su) * (storei - su);sb /= 200;ssb = sqrt(sb);此次選取yt=90, y2 =92,y94,y4=96,y5 =98, y6=10

7、0,y7 =102,y8=104,y9 =106, yg =108，它們將實軸分成11個互不相交的區(qū)間，從而將樣本值分成11組。程序統(tǒng)計了每組中的樣本數(shù)量。為方便計算，程序還計算出了山- ?)int segments12, m = 2;double x1 = 90, x10 = 108;memset(segments, 0, sizeof (segments);for ( int i = 0; i != 200; +i)if (storei <= x1)+segments0;else if (storei > x10)+segments10;else+segments int (s

8、torei - x1) / m + 1);cout << 'i'<< 't'<< "ni" << endl;for ( int i = 0; i != 11; +i)cout << i + 1 <<'t'<< segmentsi;if (i < 10)cout << 't'<< fixed << setprecision(2) << (90 + i * 2 - su) / s

9、sb;cout << endl;程序的最終運行輸出如圖2- 1所示。樣本均值:99.2773 樣本方差匕16?8302kill1ni1215-2.26 -1.7?315-1.29424-0.80529-0316360.187370.668291.1591?1.641032.13114圖2 1最終輸出結(jié)果對結(jié)果的統(tǒng)計如表 2 1所示。由表2- 1中可見E2 =1 3.3 8 0,今k=11,m = 2,并令2 2a =0.05 ,則3a(k -m T尸0.05 (8 )=15.507.由于<，故可認為產(chǎn)生的隨機數(shù)服從正態(tài)分布。表2 1i厲亠yJni?i25 -200 ?)200

10、 ?1(_oo,9020.0 1190.0612(90, 9220.0 2652.0553(92, 941 00.0 6010.0344(94, 961 20.1 1 345.0295(96, 98370.1 6640.4166(98,100460.2 0690.5167(100,102420.1 7401.4898(102,104230.1 2950.3259(104,1061 50.0 7460.00010(106,1081 00.0 3391.52911(108，七左10.0 1 661.62113.380利用中心極限定理的方法雖然可以得到服從正態(tài)分布的隨機數(shù)樣本，其思想也較為簡 單，容

11、易想到。但是這種方法每次都要先產(chǎn)生若干個服從均勻分布的隨機數(shù)樣本并求它們的 和，因而算法的時間復(fù)雜度高。(2) Box Muller方法的實現(xiàn)與分析 使用Box Muller方法得到隨機數(shù)的函數(shù)如下： double getRand()double u1 = double (rand() % 10000) / 10000, u2 =double (rand() % 10000) / 10000, r;r = 20 + 5 * sqrt(-2.0 * (log(u1) / log(e) * cos(2 * pi * u2);return r;用此函數(shù)得到的隨機數(shù)樣本理論上服從N(20, 25)。所

12、實現(xiàn)的程序產(chǎn)生了 500個隨機變量的樣本其他與利用中心極限定理的實現(xiàn)基本相同。程序得到的最終結(jié)果如圖2-2所示。樣本均值：20.0566 樣本方差：2石3即k：ll m：2Bi1ni117211-1.57335-1,1044?558-0.40664-0-917920,388610-7?9551-1610301-S51130請按任意鍵繼續(xù)-圖2-2對結(jié)果的統(tǒng)計如表 2-2所示。表2 2i宀一i, yJ?i2(rii -500 ?i )500 ?i1(q,101 70.0 2501.622(10,121 10.0 3321.893(12,14350.0 6080.694(14,16470.0 9580.015(16,18580.1 2980.736(18,20640.1 5141.817(20,22920.1 5203.368(22,24610.1 3140.339(24,26550.0 9860.6510(26,28300.0 7241.0611(28,-He300.0 6060.0012.19z由表 2 2 可見 7-2 =1 2.1 9 今 k =11, m = 2,并令 a =0 . 0,則2 2