正态分布随机数生成算法

1、概率论与数理统计课程设计题目：正态分布随机数生成算法要编程得到服从均匀分布的伪随机数是容易的。c语言、java语言等都提供了相应的函数。但是要想生成服从正态分布的随机数就没那么容易了。得到服从正态分布的随机数的基本思想是先得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。接下来就先分析三个从均匀分布到正态分布转变的方 法。然后编程实现其中的两个方法并对程序实现运作的效果进行统计分析。1、方法分析(1) 利用分布函数的反函数若要得到分布函数为 f(x)的随机变量y。可令y = f(u),其中u是服从均匀分布的随机变量，有p(y 0(i =1,2|)，则对一切 xwr有21 1

2、 e)广 1 二lim p l x xi -n x =1 -=e 2dtvvncr)2 -gcv21/.因此，对于服从均匀分布的随机变量xi,只要n充分大，随机变量xi-nn i、n 0.i m就服从n(0,1 )。我们将实现这一方法。(3) 使用box muller方法- -x2先证明广ehdx=2n：2二令 i = j e 2 dx ,则i2200e 2 1.joo2:-=y-dx e 200dy =匚edxdy令 x = r cos9, y = r sin 8,则有r2i22二二干00e rdrdur2=2冗 krdr =2n。接下来再来得出 box muller方法:设(x,y )为一

3、对相互独立的服从正态分布的随机变量。则有概率密度函数1fxy x，y =ex2 y2,则r有分布函数:r 10冗2 )令 x = rcos 工 y = rsin ，其中 i e 2 udud【-e 2 udu = 1 - e令 fr r =1 -e 2r = frz =2ln 1-z如果x服从均匀分布，则 r的分布函数即为 fr(r)。最后，可以ui用代替(1-z),令日为2叫2,其中uiu(0,1), u2(0,1),得:x =rcosu -、，-2ln u1 cos 2二u2 ,丫 =rsim -2-3 sin 2二u2从而，只需要有两个服从均匀分布的随机变量u1,u2,就能通过公式x =

4、rcos? -2lnu1 cos 2 ：u 2来得到一个服从正态分布的随机变量x。用box muller方法来生成服从正态分布的随机数是十分快捷方便的。我们也将实现这一方法。2、实现与分析(1)利用中心极限定理方法的实现与分析利用中心极限定理来生成随机数的函数(c+语言)编写如下:const int n = 200;double getrand() double s = 0;for ( int i = 0; i != n; +i) s += double (rand() % 1000) / 1000;return s; 函数生成的随机数是n个0,1间服从均匀分布的随机数的和。这里 n为200。

5、从而理论上产生的随机数应近似服从n(nn,n2),其中n为n,即200, n为0.5,仃2为1/12。程序生成了 200个随机数，并求出样本均值与样本方差，也即n与o 2的最大似然估计：/生成随机数并存储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)

6、sb += (storei - su) * (storei - su);sb /= 200;ssb = sqrt(sb);此次选取 y =90, y2 =92, =94, y4=96, y =98, y6 =100, y? =102, y8 =104,y9 =106,y10 =108,它们将实轴分成11个互不相交的区间，从而将样本值分成11组。程序统计了每组中的样本数量。为方便计算，程序还计算出了int segments12, m = 2;double x1 = 90, x10 = 108;memset(segments, 0, sizeof (segments);for ( int i =

7、0; i != 200; +i)if (storei x10)+segments10;else+segments int (storei - 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) / ssb;cout endl;程序的最终运彳r输出如图2 1所示。样本均值:993样本方差二1678302kill1n i11-2.2625-1.7?315-1.29424-0.8

8、0529-0.316360.187378291.1591?1.641032.13114图2 1最终输出结果对结果的统计如表 2 1所示。由表21中可见 * =13.380 ,今k =11,m= 2,并令豆=0.05,则 %(k-m-1 尸405(8 )=15.507.由于 13.380 15.507,故可认为产生的 随机数服从正态分布。表21 i(yi,yni?i2(n-200 ?)200 ?i1s9020.01190.0612(90,9220.02652.0553(92,94100.06010.0344(94,96120.11345.0295(96,98370.16640.4166(98,1

9、00460.20690.5167(100,102420.17401.4898(102,104230.12950.3259(104,106150.07460.00010(106,108100.03391.52911(108, f10.01661.62113.380利用中心极限定理的方法虽然可以得到服从正态分布的随机数样本，其思想也较为简单，容易想到。但是这种方法每次都要先产生若干个服从均匀分布的随机数样本并求它们的 和，因而算法的时间复杂度高。(2)box muller方法的实现与分析使用box muller方法得到随机数的函数如下：double getrand()double u1 = dou

10、ble (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;用此函数得到的随机数样本理论上服从n (20,25)。所实现的程序产生了500个随机变量的样本其他与利用中心极限定理的实现基本相同。程序得到的最终结果如图2 2所示。样本均值：200566 样本方差：26.387milm 二 2i，1ni117-1.96211-1-57335-1,1944?-0.7?5580 .40

11、664-0.01?920.388610.7?9551.16103g1.551130图2-2对结果的统计如表 2-2所示。表22 i(ynyjni?i- 2(n 500?)500?1s10170.02501.622(10,12110.03321.893(12,14350.06080.694(14,16470.09580.015(16,18580.12980.736(18,20640.15141.817(20, 22920.15203.368(22, 24610.13140.339(24, 26550.09860.6510(26, 28300.07241.0611(28,收300.06060.0012.19由表 2-2 可见 * =12.1 9 今 k=11,m = 2,并令 口=0.0,5 则 /(k-m-1 )=22.05 (8 ) =