随机数产生原理

显然乘同余法产生的序列达不到周期 M（不满足（1））。当
取 m2k （k为任意整数）时，因为 M 只有一个素因子2，
且4是 M 的因子，则由条件（2）、（3）有
，
从而a混合1(同m 余发o4生d)器达到最大周期的算法为：
x yii (y4ic/ 2k1)yi1(2d1)(m2ok)d （3.4.8)
这里 x 是积分定义域中的均匀分布随机数，这是革命性的一个视角。 从这个视角，我们把繁难的积分计算变成了简单的算术平均，而大 量的抽样计算又是计算机的拿手好戏，更重要的是当维数增加时并 不增加计算难度，从而用 Monte-Carlo 方法研究高维积分问题已 是当今计算数学界的热门课题。
2、管理科学的系统仿真研究
我们知道对于一个来自均匀分布的随机序列，应该满足独立性、 均匀性等统计特性，但伪随机数往往有一些缺陷，例如 (4.3.7) 序列到一定长度后，又开始重复下面的序列，这称为 周期性是一种明显的规律，与随机性矛盾。通常我们只能选用 一个周期内的序列作为我们的伪随机数。因此研究一种算法， 使得其产生的随机序列的周期尽可能长，我们可以通过调节 （4.3.1）的参数来实现。
F1(r)
(4.2.3)
例4.2.1 求指数分布的随机数。令
Rxetd t(ex1) 0
1Rexl1 n R ) ( x x l1 n R ) (
从而我们用服从[0，1]上的随机数R，通过上面的公式就可以得到 指数分布的随机数了。
例4.2.1 产生1000个均匀分布随机数，利用变换产生λ=6的指数 分布并进行拟合优度检验。
的特例。
xynn
yn1 yn2 yn(moM d)/
M
(4.3.6)
当 M=8 时，取初值得“菲波那西”数列。
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…… (4.3.7)
再除以模 M 我们可得到 (0,1) 之间的序列 。
的坐标。
显然末端距随链的不同而不 同，即为随机变量。
2
红链的末端距
1 0
绿链的末端距
二、伪随机数产生原理
前面Monte-Carlo方法的例子是以高质量的随机数为 基础的。通过完全的随机抽样或调查可以产生随机序 列。
x1,x2,,xn
当我们用Monte-Carlo方法研究一个实际问题时，我 们需要快速地获得大量的随机数。用计算机产生这样 的随机数是非常方便的，用数学方法在计算机上产生 的随机数称为伪随机数。
四、其他分布随机数的产生
1、 直接抽样法
由基本定理我们知道，对于有些随机变量可以建立 与R的函数关系，因此只需对R进行抽样，利用函数的 映射关系我们就可以方便地得到该随机变量的抽样了。 如前面的指数分布随机数。
例：用计算机模拟高分子链
链的末端距
末端距：空间一链的末端与 始端的距离为末端距，由于 我们将始端放在坐标原点， 所以末端距的 计算公式为：
蓝链的末端距
-1.5
0 -0.5 -1
-0.5
0 三根0.链5 的起点 1 (0,0,0)
末端距=（X2+Y2+Z2)1/2 这里X，Y，Z为链的末端点
其中c，d为任意整数。混合同余发生器是否达到最大周期M 与初始值无关。
对于乘同余发生器，由同余运算的定义，知其由如下性质
(1) 如果 a ib i(m M )o i( d 1 ,2 , ) 则有： a1b1a2b2(moMd) a1a2b1b2(moMd)
（2）如果 cac(bmoMd), 则
abmod(c,MM)
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)
即是f（x）的均值，对于均值我们有一个很好的估计，即
N
f (xi )
Eˆ[ f (x)] i1 N
【例4.1.1】 用Monte-Carlo 对 1 x 2dx 积分 0
解：将积分区域和值域看成是一个边长为一的正方形。利用均匀分 布随机数将点撒在正方形中，计算小于函数的个数并除全部点数。 这就是积分的近似值。
其中（c，M）是c，M的最大公约数。
利用这些性质可得到以下定理。
定理4.3.2 对乘同余发生器，若 (x0,M)1 ，则使
aV 1(moM d)
成立得最小正整数 V 就是此发生器得周期。
在数论中称 V 为 a 关于 M 的阶数，对于乘同余发生器，若种
子与 M 互素，则其周期就是关于 M 的阶数。这样一来，寻找 达到最大周期的同余发生器的问题就转化为数论方面寻求 M达 到最大阶数 a 的问题了。Knuth 对这一问题的研究作了总结。
% 利用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
为了达到快速的要求，一般采用递推公式
x n f ( a 0 ,a 1 , ,a k |x n 1 ,x n 2 , x n k 1 )（4.3.1)
目前最常用的方法是ห้องสมุดไป่ตู้述方法的一个特例：
混合同余法
xynn
ayn1 b yn(moM d)/
M
(4.3.2)
其中a，b，M 以及初值 y 都是正整数，容易看出 x 满足 0≤x≤1。其中 mod M 运算定义为：任一整数 y 可唯一表示为
r F (x ) P x (4.2.1)
以 F1(r) 表示随机变量 R 的分布函数，则有
F 1 ( x ) P R r P f ( ) r
0
= PXF1(r) r
1
证毕
当r0 当0r1 当r1
(4.2.2)
图4.2.1
基本定理给出了任一随机变量和均匀分布R之间的关系。而有 些随机变量可以通过分布函数的逆变换来获得，因此如果我们 可以产生高质量的均匀分布，我们就可以通过变换获得高质量 的其他分布。见公式 (4.2.3)
管理科学中的系统仿真研究，如服务系统、库存系统等。其共性就 是研究的对象是随机数，如顾客到达时间一般是一个服从指数分布 的随机数，而服务时间也可以看成是服从某种分布的随机数，当一 个系统是多队多服务体系时，问题就变的相当复杂了。我们很难用 数学的解析式来表达。这时Monte-Carlo方法也是有利的武器。对 于 这 个 领 域 的 已 有 各 种 比 较 成 熟 的 专 用 软 件 如 GPSS 、 SIMULATION等可以使用。
因此如何来获得一个周期比较长的序列，就成了我们研究的一 个内容：有关伪随机数序列的周期有如下的一些结论：
定理 4.3.1 混合同余法产生序列达最大周期 M 的充要条件：
（1） b 与 M 互素
（2） 对于 M 的任意素因子 p，有 a1(mopd)
（3） 如果 4 是 M 的因子，则 a1(mo4d)
3. 物理化学中的分子领域
50年代科学家已经在高分子领域使用Monte-Carlo方法了。 这一领域所研究的问题更加复杂，计算量非常之大。高分子 材料是由几乎是无穷的高分子链组成，而每一个链的长度又 是10的好几次方。而链的构象又是千差万别，而且是随机游 动的。如何在中微观上几乎是无规律的现象中去判断其宏观 的性质？用数学的解析式来说明这样的现象是苍白无力的。 Monte-Carlo 方 法 是 一 个 很 好 的 工 具 ， 它 使 得 科 学 家 用 Monte-Carlo方法去探索高分子运动的规律。一个典型的例 子是：对于高分子多链体的研究这是一个很复杂的问题，直 到标度理论和重正化群理论方法的引入，才使得单链构象统 计问题得到了较好的解决。
f(x)dx这里 x{x1,x2, ,xn}
这是一个众所周知的积分公式，我们当然也可以把它一般的 看为是一个高维积分，如果从传统的数值计算方法来看待， 则高维问题是随着维数的增加计算量成指数增加的，计算很 快就失去控制。但是如果我们换一个角度来看待这个问题， 从概率论的角度，它实际是：
E[f(x)] f(x)dx
基本定理：
如果随机变量X的分布函数F（x）连续，则R = F（x）是
[0，1]上的均匀分布的随机变量。
证：因为分布函数F（x）是在（0，1）上取值单调递增的连续
函数，所以当X在（－∞，x）内取值时，随机变量R则在（0， F(x)）上取值，且对应于（0，1）上的一个R值，至少有一个 x满足，见图4.2.1
1、数值计算的研究
数值计算的研究可以说是一切Monte-Carlo应用的基础，在 计算数学领域我们遇到了很多的复杂计算，一个典型的例子是 计算积分。对于一维、二维的问题早已获得解决。但是当遇到 高维积分问题时，我们传统的数值方法都由于计算量太大而陷 于了困境。但是高维积分问题又偏偏是物理、高分子化学、运 筹学和最优化问题迫切而必须解决的问题。我们来看一个例子。
公式
ynM z 则 y(mM o)d z
乘同余法
当 b = 0 时，有
xynn
ayn1 yn(moM d)/
M
加同余法 以下形式的同余法称为加同余法
xyn n yyn n(1m ynM o2 )d/ M ynk1
(4.3.4) (3.4.5)
例4.3.1 历史上比较有名的称为“菲波那西”数列为加同余 法
第四章 随机数产生原理
一、引言 二、伪随机数产生原理 三、[0，1]均匀分布随机数的算法 四、其他分布随机数的产生 五、正态分布随机数的产生 六、MATLAB统计库中的随机数发生器 七、随机数的检验 八、案例3 九、习题
一、引言
以随机数产生为基础的Monte-Carlo方法已成 为现代科研的重要手段之一。其意义早以超出 了概率论与数理统计的范畴。广泛应用于计算 方法、随机数规划、管理科学、物理化学中的 高分子结构的研究等领域。我们来看一些例子。
x=rand(2,100000);
% 产生两组随机数
Sin_xy= sin(x(1,:).^2+x(2,:).^2); % 代入函数
JF_Sin_xy=mean(Sin_xy) % 用Monte-Carlo方法求积分值
计算结果为：
Q = 0.5613 JF_Sin_xy = 0.5612 当抽样数很大时结果很接近。我们可以从Monte-Carlo方法中看出， 如果维数增加实际计算难度并不增加，因此是处理高维问题的有效 方法。
从算法上我们知道公式是递推的，因此一般的随机发生器程 序都要预先赋初值，这种初值为种（Seed），有些统计软件 如SPSS要求用户给出Seed.
以均匀分布（0，1）随机变量R变换成的随机变量。以r,ε,u, 分别表示 (0,1) 均匀分布,指数分布，N(0,1) 标准正态分布。 其他常见的分布如卡方分布、F分布等的抽样方法见表4.3.1。
面积计算结果为： s = 0.3482
【例4.1.2】 利用Monte-Carlo方法计算定积分。
sinx(2 y2)dxdy
0x1,0y1
解：抽两组随机数，求每组元素的平方代入给定的函数，然后求平
均值即得积分的近似值。
% Monte-Carlo方法积分二重积分并与数值方法的结果进行比较
Q = dblquad('sin(x.^2+y.^2)',0,1,0,1) % 数值求积分命令
kstest(ex, [ex expcdf(ex, 6)]) % 拟合优度检验
结果为：
H=0, 接受原假设，变换后的确为λ=6的指数分布
三、 （0，1）均匀分布随机数的产生
1、算法要求 （1） 产生的数值序列要具有均匀总体简
单子样的一些概率统计特性，通常包括分布 的均匀性，抽样的随机性，试验的独立性和 前后的一致性。 （2） 产生的随机数要有足够长的周期。 （3） 产生随机数速度快，占用内存小。