正态分布随机数生成算法

概率论与数理统计课程设计

题目：正态分布随机数生成算法

要编程得到服从均匀分布的伪随机数是容易的。C语言、Java语言等都提供了相应的函

数。但是要想生成服从正态分布的随机数就没那么容易了。

得到服从正态分布的随机数的基本思想是先得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。接下来就先分析三个从均匀分布到正态分布转变的方法。然后编程实现其中的两个方法并对程序实现运作的效果进行统计分析。

1、 方法分析

（1） 利用分布函数的反函数

若要得到分布函数为F(x)的随机变量Y 。

可令1()Y F u -=， 其中u 是服从均匀分布的随机变量，有

1

()(())()

P Y y P U F y F y -≤=≤=

因而，对于任意的分布函数，只要求出它的反函数，就可以由服从均匀分布的随机变量实例来生成服从该分布函数的随机变量实例。

现在来看正态分布的分布函数，对于2

~(,)X N μσ，其分布函数为：

2

2()21

()t x F x e μσ

---∞

=

⎰

显然，要想求其反函数是相当困难的，同时要想编程实现也很复杂。可见，用此种方法来生成服从正态分布的随机变量实例并不可取。

（2） 利用中心极限定理

第二种方法利用林德伯格—莱维(Lindeberg —Levi)中心极限定理：如果随机变量序列

12,,,,n X X X 独立同分布，并且具有有限的数学期望和方差

()()2

,0(1,2,),i i E X D X i μσ

==>= 则对一切x R ∈有

2

2

1lim t

n

x i n i P X n x dt μ-

-∞

→∞

=⎛

⎫⎫

-≤=

⎪⎪⎭⎝⎭

∑⎰

因此，对于服从均匀分布的随机变量i X ，只要n 充分大，

11

n i i X n μ=⎫

-⎪

⎭∑就服从()0,1N 。我们将实现这一方法。

（3） 使用Box Muller 方法 先证明2

2

2x

e

dx π-∞-∞

=⎰

：

令2

2

x

I e

dx -∞-∞

=

⎰

，则

2

2

22

2

2

2

2

x

y

x y I

e

dx e

dy e

dxdy --+∞∞∞-

-∞

-∞

-∞

=

=

⎰

⎰

⎰

令cos ,sin x r y r θθ==，则有

2

2

22

2

2

22r

r

I

e

rdrd e

rdr πθπ

π∞∞-

-

=

==⎰⎰

⎰

。

接下来再来得出Box Muller 方法：

设(),X Y 为一对相互独立的服从正态分布的随机变量。则有概率密度函数

()()2

2

2

,1,2x y X Y f x y e π

+-=

令cos ,sin x R y R θθ==,其中[]~0,2θπ，则R 有分布函数：

2

2

2

22

2

2

1()12u

u

r

r r P R r e

udud e

udu e

π

θπ

---≤=

=

=-⎰⎰

⎰

令()2

2

1r

R F r e

-=-

()1

R

R F Z -==

如果X 服从均匀分布，则R 的分布函数即为()R F r 。

最后，可以1U 用代替()1Z -，令θ为22U π，其中()1~0,1U U ，()2~0,1U ，得:

()()22cos 2,sin 2X R U Y R U θπθπ==

==

从而，只需要有两个服从均匀分布的随机变量12,U U ，就能通过公式

()2cos 2X R U θπ==

来得到一个服从正态分布的随机变量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。从而理论上产生的随机数应近似服从2(,)N n n μσ，其中n 为N ，即200，μ为0.5，2

σ为1/12。程序生成了200个随机数，并求出样本均值与样本方差，也即μ与2

σ的最大似然估计：

//生成随机数并存储 double sum, store[200], xi, su = 0, sb = 0, ssb = 0; int cnt = 0; sum = 0;

for (int i = 0; i != 200; ++i) { xi = getRand(); sum += xi; store[i] = xi;

}

//得到样本均匀与样本方差 su = sum / 200;

for (int i = 0; i != 200; ++i)

sb += (store[i] - su) * (store[i] - su); sb /= 200;

ssb = sqrt(sb);

此次选取1234567890,92,94,96,98,100,102,104,y y y y y y y y ========

910106,108y y ==，它们将实轴分成11个互不相交的区间，从而将样本值分成11组。

程序统计了每组中的样本数量。为方便计算，程序还计算出了

ˆ()ˆi y μσ

-：

int segments[12], m = 2; double x1 = 90, x10 = 108;

memset(segments, 0, sizeof (segments)); for (int i = 0; i != 200; ++i) { if (store[i] <= x1)

++segments[0];

else if (store[i] > x10)