随机数产生原理

bi (mod M )( i 1,2,)
则有：
a1 b1 a2 b2 (mod M ) a1a2 b1b2 (mod M )
（2）如果
ca cb(mod M ),
M a b mod (c , M )
则
其中（c，M）是c，M的最大公约数。
证：因为分布函数F（x）是在（0，1）上取值单调递增的连续 函数，所以当X在（－∞，x）内取值时，随机变量R则在（0， F(x)）上取值，且对应于（0，1）上的一个R值，至少有一个 x满足，见图4.2.1
r F ( x ) P x
(4.2.1)
以
F1 (r ) 表示随机变量 R 的分布函数，则有
当 M=8 时，取初值得“菲波那西”数列。
的特例。
(4.3.6)
0，1，1，2，3，5，8，13，21，34，55，89，144， 253 ……
对上述数列取模得： 0，1，1，2，3，5，0，5，5，7，1，1…… 再除以模 M 我们可得到 (0,1) 之间的序列 。 (4.3.7)
我们知道对于一个来自均匀分布的随机序列，应该满足独立性、 均匀性等统计特性，但伪随机数往往有一些缺陷，例如 (4.3.7) 序列到一定长度后，又开始重复下面的序列，这称为
计算结果为： Q = 0.5613 JF_Sin_xy = 0.5612 当抽样数很大时结果很接近。我们可以从Monte-Carlo方法中看出， 如果维数增加实际计算难度并不增加，因此是处理高维问题的有效 方法。
这里 x 是积分定义域中的均匀分布随机数，这是革命性的一个视角。 从这个视角，我们把繁难的积分计算变成了简单的算术平均，而大 量的抽样计算又是计算机的拿手好戏，更重要的是当维数增加时并
例4.2.1 产生1000个均匀分布随机数，利用变换产生λ=6的指数 分布并进行拟合优度检验。 clc,clear x=linspace(0,20,100); R=rand(1,1000); % 产生1000个（0，1）均匀随机数 ex=-6*log(1-R); % 变换为指数分布随机数 ex=ex'; m=mean(ex) v=var(ex) subplot(1,2,1),cdfplot(ex) subplot(1,2,2),hist(ex) kstest(ex, [ex expcdf(ex, 6)]) % 拟合优度检验 结果为：
是10的好几次方。而链的构象又是千差万别，而且是随机游 动的。如何在中微观上几乎是无规律的现象中去判断其宏观 的性质？用数学的解析式来说明这样的现象是苍白无力的。 Monte-Carlo 方 法 是 一 个 很 好 的 工 具 ， 它 使 得 科 学 家 用
Monte-Carlo方法去探索高分子运动的规律。一个典型的例
H=0,
接受原假设，变换后的确为λ=6的指数分布
三、 （0，1）均匀分布随机数的产生


1、算法要求 （ 1） 产生的数值序列要具有均匀总体简 单子样的一些概率统计特性，通常包括分布 的均匀性，抽样的随机性，试验的独立性和 前后的一致性。 （2） 产生的随机数要有足够长的周期。 （3） 产生随机数速度快，占用内存小。
-0.5
0.5 三根链的起点 (0,0,0)
0

末端距=（X2+Y2+Z2)1/2 这里X，Y，Z为链的末端点 的坐标。 显然末端距随链的不同而不
红链的末端距
1
0

绿链的末端距
同，即为随机变量。
二、伪随机数产生原理

前面 Monte-Carlo 方法的例子是以高质量的随机数为 基础的。通过完全的随机抽样或调查可以产生随机序 列。
周期性是一种明显的规律，与随机性矛盾。通常我们只能选用
一个周期内的序列作为我们的伪随机数。因此研究一种算法， 使得其产生的随机序列的周期尽可能长，我们可以通过调节
（4.3.1）的参数来实现。
因此如何来获得一个周期比较长的序列，就成了我们研究的一 个内容：有关伪随机数序列的周期有如下的一些结论：
定理 4.3.1 混合同余法产生序列达最大周期 M 的充要条件：
E[ f ( x )] f ( x )dx
即是f（x）的均值，对于均值我们有一个很好的估计，即
N
ˆ [ f ( x )] E
f (x )
i 1 i
N
【例4.1.1】 用Monte-Carlo 对

1
0
x dx
2
积分
解：将积分区域和值域看成是一个边长为一的正方形。利用均匀分 布随机数将点撒在正方形中，计算小于函数的个数并除全部点数。 这就是积分的近似值。
其中a，b，M 以及初值 y 都是正整数，容易看出 x 满足 0≤x≤1。其中 mod M 运算定义为：任一整数 y 可唯一表示为 公式
y nM z
则
y(mod M ) z
乘同余法 当 b = 0 时，有
yn ayn1 xn yn (mod M ) / M
加同余法 以下形式的同余法称为加同余法
面积计算结果为： s = 0.3482
【例4.1.2】 利用Monte-Carlo方法计算定积分。
0 x 1, 0 y 1

sin( x 2 y 2 )dxdy
解：抽两组随机数，求每组元素的平方代入给定的函数，然后求平 均值即得积分的近似值。 % Monte-Carlo方法积分二重积分并与数值方法的结果进行比较 Q = dblquad('sin(x.^2+y.^2)',0,1,0,1) % 数值求积分命令 x=rand(2,100000); % 产生两组随机数 Sin_xy= sin(x(1,:).^2+x(2,:).^2); % 代入函数 JF_Sin_xy=mean(Sin_xy) % 用Monte-Carlo方法求积分值
f ( x)dx
这里
x { x1 , x2 ,, xn }
这是一个众所周知的积分公式，我们当然也可以把它一般的
看为是一个高维积分，如果从传统的数值计算方法来看待， 则高维问题是随着维数的增加计算量成指数增加的，计算很
快就失去控制。但是如果我们换一个角度来看待这个问题，
从概率论的角度，它实际是：
于 这 个 领 域 的 已 有 各 种 比 较 成 熟 的 专 用 软 件 如 GPSS 、
SIMULATION等可以使用。
3. 物理化学中的分子领域
50年代科学家已经在高分子领域使用 Monte-Carlo方法了。 这一领域所研究的问题更加复杂，计算量非常之大。高分子
材料是由几乎是无穷的高分子链组成，而每一个链的长度又
a 1(mod 4)
yi (4c 1) yi 1 ( 2d 1)(mod 2 k ) k x y / 2 i i
（3.4.8)
其中c，d为任意整数。混合同余发生器是否达到最大周期M
与初始值无关。 对于乘同余发生器，由同余运算的定义，知其由如下性质 (1) 如果 ai
第四章 随机数产生原理
一、引言 二、伪随机数产生原理 三、[0，1]均匀分布随机数的算法 四、其他分布随机数的产生 五、正态分布随机数的产生 六、MATLAB统计库中的随机数发生器 七、随机数的检验 八、案例3 九、习题
一、引言

以随机数产生为基础的 Monte-Carlo 方法已成 为现代科研的重要手段之一。其意义早以超出
% 利用Monte-Carlo方法计算定积分 x=rand(1,1000); x_2=x.^2; JF=mean(x_2) % 作Monte-Carlo积分示意图 for i=1:1000 xx=rand(1,100); yy=rand(1,100); end x1=linspace(0,1,50); y1=x1.^2; plot(xx,yy,'.',x1,y1,'linewidth',2) axis equal h = legend('x-y','x^2'); JF = 0.3346
为了达到快速的要求，一般采用递推公式
xn f (a0 , a1 ,, ak | xn1 , xn 2 , xn k 1 )
目前最常用的方法是上述方法的一个特例： 混合同余法
（4.3.1)
yn ayn1 b xn yn (mod M ) / M
(4.3.2)
的其他分布。见公式 (4.2.3)
F (r )
1
(4.2.3)
例4.2.1 求指数分布的随机数。令
R e t dt (e x 1)
0
x
1 R e
x
ln( 1 R) x x ln( 1 R)
从而我们用服从[0，1]上的随机数R，通过上面的公式就可以得到 指数分布的随机数了。
x1 , x2 ,, xn

当我们用 Monte-Carlo方法研究一个实际问题时，我 们需要快速地获得大量的随机数。用计算机产生这样 的随机数是非常方便的，用数学方法在计算机上产生 的随机数称为伪随机数。
基本定理：
如果随机变量X的分布函数F（x）连续，则R = F（x）是 [0，1]上的均匀分布的随机变量。
（1） b 与 M 互素
（2） 对于 M 的任意素因子 p，有 （3） 如果 4 是 M 的因子，则
a 1(mod p) a 1(mod 4)
显然乘同余法产生的序列达不到周期 M（不满足（1））。当 取
m 2k
（k为任意整数）时，因为 M 只有一个素因子2， ，
且4是 M 的因子，则由条件（2）、（3）有 从而混合同余发生器达到最大周期的算法为：
不增加计算难度，从而用 Monte-Carlo 方法研究高维积分问题已
是当今计算数学界的热门课题。
2、管理科学的系统仿真研究
管理科学中的系统仿真研究，如服务系统、库存系统等。其共性就 是研究的对象是随机数，如顾客到达时间一般是一个服从指数分布 的随机数，而服务时间也可以看成是服从某种分布的随机数，当一 个系统是多队多服务体系时，问题就变的相当复杂了。我们很难用 数学的解析式来表达。这时Monte-Carlo方法也是有利的武器。对
(4.3.4)
y n y n 1 y n 2 y n k 1 xn yn (mod M ) / M
(3.4.5)
例4.3.1 历史上比较有名的称为“菲波那西”数列为加同余 法
y n y n 1 y n 2 xn yn (mod M ) / M
了概率论与数理统计的范畴。广泛应用于计算
方法、随机数规划、管理科学、物理化学中的
wk.baidu.com
高分子结构的研究等领域。我们来看一些例子。
1、数值计算的研究
数值计算的研究可以说是一切 Monte-Carlo 应用的基础，在
计算数学领域我们遇到了很多的复杂计算，一个典型的例子是
计算积分。对于一维、二维的问题早已获得解决。但是当遇到 高维积分问题时，我们传统的数值方法都由于计算量太大而陷 于了困境。但是高维积分问题又偏偏是物理、高分子化学、运 筹学和最优化问题迫切而必须解决的问题。我们来看一个例子。
子是：对于高分子多链体的研究这是一个很复杂的问题，直 到标度理论和重正化群理论方法的引入，才使得单链构象统
计问题得到了较好的解决。
例：用计算机模拟高分子链
链的末端距
蓝链的末端距

末端距：空间一链的末端与
始端的距离为末端距，由于 我们将始端放在坐标原点， 所以末端距的 计算公式为：
2
0 -0.5 -1 -1.5
F1 ( x ) PR r P f ( ) r
0 1 = P X F (r ) r 1


当r 0 当0 r 1 当r 1
(4.2.2)
证毕
图4.2.1
基本定理给出了任一随机变量和均匀分布R之间的关系。而有 些随机变量可以通过分布函数的逆变换来获得，因此如果我们 可以产生高质量的均匀分布，我们就可以通过变换获得高质量
利用这些性质可得到以下定理。 定理4.3.2 对乘同余发生器，若 ( x0 , M ) 1 ，则使
aV 1(mod M )
成立得最小正整数 V 就是此发生器得周期。
在数论中称 V 为 a 关于 M 的阶数，对于乘同余发生器，若种