正态分布随机数生成算法

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

2、F (y)因而，对于任意的分布函数，只要求出它的反函数，就可以由服从均匀分布的随机变量实例来生成服从该分布函数的随机变量实例。现在来看正态分布的分布函数，对于XN(,；2)，其分布函数为：1x显然，要想求其反函数是相当困难的，同时要想编程实现也很复杂。可见，用此种方法来生成服从正态分布的随机变量实例并不可取。(2)利用中心极限定理第二种方法利用林德伯格一莱维(Lindeberg Levi)中心极限定理：如果随机变量序列X- X2,，独立同分布，并且具有有限的数学期望和方差2E Xi =，D Xi =：.:0(i=1,2,)，则对一切R 有因此，对于服从均匀分布的随机变量 Xi，只要n充分大，随

3、机变量- Xi - v J 就服从N 0,1。我们将实现这一方法。(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方法:设X ,Y为一对相互独立的服从正态分布的随机变量。则有概率密度函数x2：；y2令x = R cos h y二R sin二，其中v 0, 2，贝U R有分布函数: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的分布函数即为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来得到一个服从正态分布的随机变量X。用Box Muller方法来生成服从正态分布的随机数是十分快捷方便的。我们也将实现这一方法。2、实现

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

6、数并存储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个互不相交的区间，从而将样本值分成11组。程序统计了每组中的样本数量。为方便计算，程序还计算出了山- ?)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最终输出结果对结果的统计如表 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.由于<，故可认为产生的随机数服从正态分布。表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利用中心极限定理的方法虽然可以得到服从正态分布的随机数样本，其思想也较为简 单，容

11、易想到。但是这种方法每次都要先产生若干个服从均匀分布的随机数样本并求它们的 和，因而算法的时间复杂度高。(2) Box Muller方法的实现与分析 使用Box Muller方法得到随机数的函数如下： 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;用此函数得到的随机数样本理论上服从N(20, 25)。所

12、实现的程序产生了 500个随机变量的样本其他与利用中心极限定理的实现基本相同。程序得到的最终结果如图2-2所示。样本均值：20.0566 样本方差：2石3即k：ll m：2Bi1ni117211-1.57335-1,1044?558-0.40664-0-917920,388610-7?9551-1610301-S51130请按任意键继续-图2-2对结果的统计如表 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