第3章
CHAPTER3
蒙特卡罗法
蒙特卡罗法正式出现在20世纪40年代,但实际上其起源可上溯到17世纪。蒙特卡罗法
的发展大体上经历了三个阶段。
蒙特卡罗法的第一阶段主要发展直接仿真的方法,也就是对具有随机性质的问题进行直
接仿真,没有形成一种普遍的和有效的解决实际问题的方法。这段时间主要的成果是发现和
论证随机变量的统计性质和对随机变量进行定量分析的理论与方法。
蒙特卡罗法的第二阶段主要发展了数学仿真方法,是蒙特卡罗法真正确立的阶段。数学
仿真方法就是针对具体的随机性质问题,把对随机现象进行直接仿真改为用数学方法进行仿
真求得解答。例如,18世纪用随机投针法计算π值,1889年研究了随机游动问题,1908年著
名统计学家Gossett推导了t 分布和一种χ2 分布,1928年Courant、Friedrichs和Lewy讨论
了二维和三维布朗运动。这种方法显然比直接仿真法有很多优越性。
蒙特卡罗法的第三阶段主要发展了拟仿真方法,该方法使蒙特卡罗法得到迅速发展和广
泛应用。为了更加有效地解决实际问题,可以人为地将原来的随机现象变为另一种与原问题
答案有一定关系的人为的“随机现象”,然后利用人为随机现象与原问题答案间的关系,再利用
计算机仿真方法得到原问题的解答,拟仿真方法的关键就是建立替代的人为随机现象。拟仿
真方法使蒙特卡罗法进入了迅速发展的阶段,可以说,没有拟仿真方法就没有蒙特卡罗法的应
用和发展。
现在,蒙特卡罗法仍然在进一步的发展中,特别是随着非线性物理学的发展、混沌现象、生
命科学和经济学的深入研究,蒙特卡罗法已经成为必备的定量分析研究的工具和仿真模拟的
理论基础。蒙特卡罗法学习的意义已经大大超出计算电磁学的领域。
3.1 蒙特卡罗法的基本原理
在自然界有相当多的物理量是具有一定程度不确定性的量。实际上,任何物理现象都是
在不断变动中的,环境因素也在不断地变化,就必然具有随机性。随机的物理量可以用概率论
和随机过程的理论研究和表达。但是,有些物理量,例如分形、混沌等非线性过程的参量,根本
无法用解析式表达而必须用计算机经不断的迭代和产生随机参量才能获得解答。目前,处理
这类问题最有效的方法就是计算机模拟,采用的工具就是蒙特卡罗法。从工程设计的角度,如
果不考虑物理量的不确定性就不会得到成功的、可靠的工程解决方案。
从科学的角度看,随机现象总是伴随着确定的物理规律,当改换角度看同一物理现象时,
65
确定的物理现象的背景很可能是随机现象。例如,加热封闭空间中的气体,温度和气压满足热
力学定律是确定不变的,但是从分子运动的角度理解加热封闭空间中的气体现象时,我们看到
的是大量的随机运动的分子,热力学规律不过是大量分子随机运动的统计结果。从自然哲学
的角度看,确定性与随机性是互为依存的同一运动的两方面,只看运动的确定性或者说只看运
动的随机性都是不完全的,都曲解了运动的本来面目。因此,本章的内容除了要介绍蒙特卡罗
法的基本原理和基本处理方法,还希望读者能够牢固地建立起基于随机性的理解和分析各种
现象的能力。如果读者能够主动地站在随机性的角度观察和理解世界,将能够看到同以往完
全不同的非常奇妙的图景,将能够发掘出随机现象中的确定性和规律性,将能够解释以前认为
是不可理解的、匪夷所思的现象。
在具体介绍蒙特卡罗法之前,笔者首先要强调,进行计算机模拟并不是最终的目的,计算
机模拟只是工具,一定要解决实际提出的工程和理论问题。因此,计算机模拟完成之后,还必
须回到原来提出的问题去,要对计算机模拟的结果进行分析研究并回到现实世界去,给出对现
实问题和现象的物理解释,否则其他人无法使用所得出的数值分析结果。这个过程的难度和
重要意义不亚于计算机模拟的全部工作,中国有句古语:行百里者半九十,说的就是这种情
形,最后的十里没有完成或没有走对,前面的九十里就毫无意义。在学习蒙特卡罗法时,这一
点尤其重要,拟仿真方法实际上是采用已有的随机模型的样板解决实际的问题,得到的解答就
必须回到原来的问题对得到的结果作出合理的物理解释。
3.1.1 蒙特卡罗法的基本过程
随机变量的全体称为总体,随机变量总数称为容量,容量可以是无限也可以是有限的。当
某个随机变量与存在的另一个随机变量之间有一定的函数关系时,可以定义后一个随机变量
为统计量。但随机变量本身并无确定的规律可循,又不可能每次都对所有的随机变量进行研
究。因此,为了研究随机变量之间的函数关系,为了在有限范围内研究在无限统计分析时才具
有的规律性,就必须首先构成样本空间。
在随机变量总体中抽出一定数量的个体来进行观测的过程称为抽样。抽样要采用科学的
方法使抽出来的个体能反映总体的统计或随机性质,反映以总体为单位时参与物理变化的情
况。例如,一棵树可以是一片森林的抽样,但是树叶就不是森林的科学抽样,一片森林对调节
地球或区域的生态有作用但是却无法讨论一棵树对生态的作用。
设随机变量总体为Ω,在保证条件无任何变化的情况下对随机变量Ω 进行n 次重复独立
观测,得到的结果称为n 次简单随机抽样,n 次抽样所得的结果记为ζ=(ξ1,ξ2,…,ξn ),称为
来自总体Ω 的随机样本,抽样次数n 称为样本容量。
样本应该尽可能多地反映出随机变量总体的特性,即每个ζi 都具有随机变量总体的特
征,也就是ζi 与Ω 具有相同的随机分布。此外,每次抽样都要保证是独立进行的,互不影响,
即(ξ1,ξ2,…,ξn)=ζi 是相互独立的随机变量。
定义3.1.1 如果随机变量(ξ1,ξ2,…,ξn)=ζi 相互独立,并且每一个ζ 与总体Ω 有相同
的概率分布,则随机变量(ξ1,ξ2,…,ξn)=ζi 称为总体Ω 的容量为n 的样本。
抽样之后的量ζi,仍是独立的随机变量,但ξ1,ξ2,…,ξn 则具有确切的数值。此时,ζi 可
以代表Ω 参加函数运算,并且由此定义出一类由随机变量的样本产生的随机变量,即统计量。
定义3.1.2 设ζi =(ξ1,ξ2,…,ξn )是总体Ω 的一个样本,G (ξ1,ξ2,…,ξn )是样本ζi =
(ξ1,ξ2,…,ξn )的一个函数,则称G (ξ1,ξ2,…,ξn )为一个统计量。因为G 是随机变量ξ1,
66
ξ2,…,ξn 的函数,所以统计量也是一个随机变量,也有自己的概率分布,称为抽样分布。求抽
样分布是数理统计的基本问题之一。
设X =(x1,x2,…,xn)是来自该统计结果的一个样本,样本的均值及方差与随机变量的
均值及方差有如下的关系。
设(x1,x2,…,xn)是来自正态分布N (μ,σ2)的一个样本,μ 和σ2 是未知的,现在要用样
本来求取μ 和σ2 的近似值。
定义样本均值为
..x =1 nΣn
i=1
xi
样本方差为
S2n
=1 nΣn
i=1 (xi -..x)2
设E(..x)=μ,即样本均值..x 是原随机变量μ 的无偏估值。由于样本方差的均值为
E(S2n)=σ2
nEnS2n
σ2
.
è .
.
. ÷
=n -1 n σ2 (3.1.1)
所以,S2n
并不是原随机变量方差σ2 的无偏估值。但是对S2n
略加修正即可写为
S2 = 1 n -1Σn
i=1 (xi -..x)2 (3.1.2)
S2 就是原随机变量方差σ2 的无偏估值。当n 趋于无穷大时S2n
趋于S2,因此称S2n
称为σ2 的
渐进无偏估计。通常使用时,由于样本容量不够大,使用S2n
会使得过小估计σ2,因此实际中
采用S2 的比较多。
蒙特卡罗法应用时首先要构成一个概率空间,然后在该概率空间中确定一个依赖随机变
量X ,其分布函数为F(x)的统计量g(x),使得g(x)的数学期望
E(g)=∫g(x)dF(x) (3.1.3)
正好等于所要求取的值G,式(3.1.3)中F(x)是x 的分布函数。最后,产生随机变量x 的简单
子样x1,x2,…,xN ,取其相应的统计量g(x1),g(x2),…,g(xN )的算术平均值N 为统计量
g(x)的无偏估值,作为G(ξi)的近似估计。令
N =1 NΣN
n=1
g(xn) (3.1.4)
若用N (X )表示对G(ξi)的估计量,则Eξ( N (X ))=G(X )称为G(ξi)的无偏估计。实
际上不可能进行无穷次实验获得渐近无偏估计G (X ),因此一般使用渐进无偏估计作为
G(X )的近似估计。
显然,蒙特卡罗法应用的关键是要人为地确定一个统计量,使其数学期望正好等于所要求
的值,称这样的统计量为无偏估计量。
蒙特卡罗法的最低要求:能确定一个与计算步数N 有关的无偏估计量N 使得当N →∞
时, N (X )依概率收敛于所要求的值G(X ),即对任意小的ε>0,有
lim N →∞
P(| N -G|<ε)=1 (3.1.5)
式中,P(*)表示事件*的概率。
蒙特卡罗法可以用来解决两类问题。
67
第一类是确定性问题,例如,计算重积分、求逆矩阵、解线性代数方程组、解积分方程、解偏
微分方程的边值问题、计算微分算子特征值等。
第二类是随机性问题,例如,中子扩散问题、运筹学的库存问题、随机服务系统中的排队问
题、信息和传输的算法、随机游动问题、动物生态竞争问题、传染病蔓延问题等。
解第一类问题时要假设一个概率事件作为拟仿真建模,使概率事件的统计结果就是要求
的值。首先,构成一个概率空间,也就是建立一个相应于所要求解问题的概率模型,使所要求
的解就是建立概率模型的某一随机变量的数学期望值。然后,在该概率空间中确定一个随机
变量x 和这个随机变量的统计量g(x),并对g(x)进行抽样和试验。通常在计算机上试验
时,要按x 的分布函数F(x)抽取随机变量x 的简单子样x1,x2,…,xN 。最后,这个统计量
的数学期望值E(g)的近似值或估计量g(x1),g(x2),…,g(xN )的算术平均值
N =1 nΣN
n=1
g(xn) (3.1.6)
就是要求的解G。
求解第二类问题时,通常没有或很难求得确定的表达式,一般都采用直接模拟的方法,即
根据实际问题的随机性用计算机进行抽样,再对试验得到的随机变量值进行统计平均,求得所
要的结果。
从以上两种问题可看出,蒙特卡罗法最关键的一步是确定一个概率空间中的统计量,即使
相应的无偏估计量的数学期望恰好等于所求的值。
3.1.2 蒙特卡罗法的基本问题
使用蒙特卡罗法一般要讨论如下一些基本问题,这些是不同于其他计算方法之处。
1.蒙特卡罗法的收敛性
蒙特卡罗法尽管可以进行相当多次的抽样,但是总不能达到无穷多次,因此通常用G 的
近似估计值N 代替数学期望值G。这就有个收敛的问题,因此蒙特卡罗法的最基本要求是
能确定这样一个与计算步数N 有关的统计估计值N ,使当N →∞时, N →期望值G。
蒙特卡罗法的近似估值N 依概率1收敛于G,即
P lim n→∞
( N =G ) =1 (3.1.7)
要满足这一点,无偏估计量g(x)的均值应满足
E(|g|)=∫|g(x)|dF(x)<+ ∞ (3.1.8)
的条件,也就是要求这个无偏估计量的绝对数学期望值存在,式中F(x)为x 的分布函数。
收敛速度:如果无偏估计量满足条件E(|g|r)=∫|g(x)|rdF(x)<+ ∞,其中1≤
r ≤2,则有
P lim N →∞
Nr-1 ( r ( N -G)=0) =1 (3.1.9)
N 依概率1收敛于G 的速度为N1-r
r ,因为使分式(1-r)/r取得最大值的r 为2,所以收敛速
度不会超过N -0.5。
2.蒙特卡罗法的误差
根据中心极限定理,只要所确定的无偏估计量具有有限的不等于零的方差σ2,则对于任
意非负的X 值均有
68
lim N →∞
P N
σ | N -G|<X
.
è .
.
. ÷
= 1
2π∫X
-Xe-t2
2dt
当N 足够大时,有
P | N -G|< Xσ
N
.
è .
.
. ÷
≈ 1
2π∫X
-Xe-t2
2dt=1-α (3.1.10)
1-α 称为置信水平,α 称为置信度,式(3.1.10)说明发生事件| n -G|<(Xσ/ N )的概
率由问题要求确定的置信水平决定。实际上可以由问题的要求得到1-α,再按正态分布确定
X ,则可以由X 求出真实值与估值的偏差为| n -G|<(Xσ/ N )。为计算方便,通常取置
信水平为0.5,0.95,0.997对应的X 为0.6745,1.96,3。置信度与正态分布取值区间关系如表
3.1所示。
表3.1 置信度与正态分布取值区间关系
α 0.0063 0.0027 0.0455 0.134 0.3173 0.617 0.50 0.05 0.02 0.01
X 4 3 2 1.5 1 0.5 0.6745 1.9600 2.3263 2.5758
3.蒙特卡罗法的费用
根据上述讨论,设问题要求误差为ε,置信水平为1-α,则可以利用正态积分表确定X ,再
由置信水平和(Xσ/ N )≤ε 确定N 。也就是要求误差小于ε 时,就必须要求样本容量N 为
N ≥ Xσ
ε
.
è .
.
. ÷
2 (3.1.11)
设每计算一次无偏估计量所需的费用为C,则误差为ε 的蒙特卡罗法总花费为
NC ≈ Xσ
ε
.
è .
.
. ÷
2
C = X
ε
.
è .
.
. ÷
2
σ2C ∝σ2C (3.1.12)
即蒙特卡罗法的费用与无偏估计量方差σ2、费用C 的乘积σ2C 成正比。在同一置信水平下,
σ2C 值越小就表示该方法越好。
3.1.3 蒙特卡罗法的特点
蒙特卡罗法有如下几个特点。
(1)收敛速度与问题维数无关。由式(3.1.12)可以看出,当置信水平确定之后,蒙特卡罗
法的误差由无偏估计量的方差σ2 和样本容量N 唯一确定。因此,蒙特卡罗法的收敛速度与
问题的维数无关,维数变化只是引起一次观察所需要的花费C 的大小变化。
(2)受问题条件限制的影响不大。采用蒙特卡罗法求积分问题或边值问题时,一旦确定
了拟仿真的概率模型边界条件和区域形状只影响概率事件的判决,问题的区域形状的限制对
蒙特卡罗法影响不大,不会明显增加问题的难度和计算花费。
(3)不必进行离散化处理。一般计算方法都首先要把问题进行离散化,进行区域切分要
确定网格,网格一旦确定后,问题只能针对网格结点上的值展开,网格以外点的值只能采用插
值方法近似求得。蒙特卡罗法完全不同,蒙特卡罗法不需要进行任何离散化处理,其计算针对
连续变化的量。这一优点可使蒙特卡罗法更加适用于高维问题,使精度进一步增加,可以节省
大量存储单元,对建立蒙特卡罗通用程序很有利。
(4)蒙特卡罗法是一种直接解决问题的方法。特别是那些具有随机性质的问题,蒙特卡
69
罗法可以直接建立随机模型并在计算机上实现针对求解问题的随机过程的实验,因此是一种
试验方法,是一种可靠的客观的方法,其结果常作为标准去衡量其他方法的正确性。
(5)误差容易确定。一般计算方法要估计计算结果与真值的误差是十分困难的事,计算
结果的误差与离散方法、计算机的有效位等因素都有关。而蒙特卡罗法的误差则只与X 、N 、
σ 有关,X 由置信度确定,N 为试验次数,σ 可以通过N 来给出:
σ^N ≈ 1 NΣN
n=1
g2(xn)- 2N
é
. êê
ù
. úú
12
(3.1.13)
由式(3.1.13)可知,蒙特卡罗法的误差很容易确定。
(6)蒙特卡罗法的缺点是对小维数的问题不如其他方法好,求得的结果一般都比真实
值低。
3.1.4 蒙特卡罗法待研究的若干问题
1.随机数
具有均匀分布的总体中所产生的简单子样称为随机数序列,每一个体称为随机数。伪随
机数则是由数学递推公式所产生的随机数。
产生伪随机数的方法有平方取中法、乘同余方法和乘加同余方法等。后两者很好地解决
了出现周期性的问题,但是仍然有许多需要研究的问题。例如,20世纪60年代初人们提出了
伪随机数序列的独立性问题。
2.已知分布的随机抽样
冯·诺依曼早在1951年就完成了随机抽样的基本理论,1954年Kahn对冯·诺依曼的理
论系统进行了补充研究并扩充了部分方法,但是此后在随机抽样的基本原理方面几乎没有什
么进展,特别是如何将原理用于若干重要分布的随机抽样得到理想的随机抽样方法的研究。
3.非归一问题的随机抽样
有些问题分布是已知的但不是归一的,其抽样方法是1953年由Metropolis确定的,但这
种方法有很大的缺点,特别对连续型分布的情况,需要提出更加理想的抽样方法。
4.蒙特卡罗法的基本技巧
蒙特卡罗法的基本技巧主要是讨论如何确定合适的无偏估计量使它相应的σ2 尽量小的
技巧。蒙特卡罗法的基本技巧随问题不同而变化,随着具体问题的应用的迅速增长,需要应用
者进一步发展各种不同类型的应用方法。
5.蒙特卡罗法的并行化计算方法
略。
3.1.5 随机变量的基本规律
1.随机变量
如果每次试验的结果都可以用ξ 来表示,ξ 为实数而且对任意实数X ,“ξ<X ”事件的发
生有着确定的概率,则称ξ 为随机变量。
事件ξ<X 对应的概率分布为P {ξ<X }=F(x),其中的随机变量可以是离散随机变量,
也可以是连续随机变量。连续随机变量ξ 可以用概率分布函数或概率分布密度p(z)来描写,
概率分布函数F(x)可以表示为
70
F(x)=∫x
-∞
p(z)dz (3.1.14)
表示随机变量ξ 出现在z 到z+dz 范围内的概率为p(z)dz,其中p(z)称为随机变量ξ
的分布密度函数,有如下性质:
p(z)≥0
∫∞
-∞
p(z)dz =1
P{a ≤ξ ≤b}=F(b)-F(a)=∫b
ap(z)dz
ì
.
í
...
...
(3.1.15)
2.数学期望值
数学期望值定义为随机变量出现概率最大的值,用E表示,即
E(ξ)=∫∞
-∞
XdF(x)=∫∞
-∞
Xp(x)dx (3.1.16)
若Fξ(x)是ξ 的分布函数,f(x)是连续函数,则有
E(f(ξ))=∫∞
-∞
f(x)dFξ(x)=∫∞
-∞
f(x)p(x)dx
3.方差
方差表示随机变量值偏离数学期望值的特性,定义为
D(ξ)=E(ξ-E(ξ))2 =∫∞
-∞ (ξ-E(ξ))2dFξ(x) (3.1.17)
D(ξ)=E(ξ2)- (E(ξ))2 (3.1.18)
4.特征函数
设随机变量ξ 的分布函数为F(x),则称eitξ 的数学期望f(t)=E(eitξ)=∫eitxdF(x)为ξ
的特征函数。
唯一性定理:分布函数由特征函数唯一决定,两者有一一对应关系。
5.中心极限定理
无论单个随机变量的分布如何,大量的独立随机变量之和总是满足正态分布的,称为高斯
分布,该分布可以由给定的数学期望值和方差完全确定下来,其分布密度函数为
f(x)= 1
σ 2πexp - (x -μ)2
2σ2
é
. êê
ù
. úú
(3.1.19)
6.分布函数的基本性质
(1)单调不减:若x1<x2 则F(x1)≤F(x2)。
(2)定义F(-∞)=lim x→-∞
F(x),F(+∞)=lim x→+∞
F(x),则有F(-∞)=0,F(+∞)=1。
(3)左连续:对于-∞<x0<+∞,有lim x→x0-0
F(x)=F(x0)。
这三条基本性质不但是分布函数要满足的必要条件,也构成随机变量分布函数的充要条
件。下面用符号N (μ,σ2)表示正态分布的概率函数,其概率密度函数为
p(x)= 1
2πσexp - (x -μ)2
2σ2
式中,σ>0。
71
7.随机变量序列的收敛性
1)依分布收敛
设随机变量ξn (n=1,2,…)和随机变量ξ 的分布函数分别为Fn (x)(n=1,2,…)和
F(x),若F(x)的所有连续点x 上都有lim n→∞
Fn(x)=F(x),则称随机变量序列{ξn}依分布收敛
于随机变量ξ,简记为ξn
L→ξ。
2)依概率收敛
若对任意给定的ε>0有lim n→∞
P{|ξn -ξ|<ε}=1或lim n→∞
P {|ξn -ξ|≥ε}=0,则称随机变量
序列{ξn}依概率收敛于随机变量ξ,简记为ξn
P→ξ。
3)r 阶收敛
若有E|ξ|r<∞和E|ξn|r<∞,且lim n→∞E|ξn -ξ|r =0(r>0,n=1,2,…),则称随机变量序
列{ξn}为依r 阶收敛于随机变量ξ,简记为ξn
r→ξ,r=1时称{ξn }为平均收敛,r=2时称
{ξn}为均方收敛。
4)依概率1收敛
若有P{ω:lim n→∞ξn(ω)=ξ(ω)}=1,简记为P {lim n→∞ξn =ξ}=1,则称随机变量序列{ξn }以概
率1(或几乎处处)收敛于随机变量ξ,记为ξn
a.e→ξ。
上述的几种收敛之间的关系可以用图3.1形象化表示,其中由左向右依次满足包含关系。
图3.1 几种收敛的关系
5)切比雪夫(Чебыщев)不等式
设随机变量ξ 的方差Dξ 为有限值,则对任意的ε>0,有
P{|ξ-Eξ|≥ε}≤ Dξ
ε2 或 P{|ξ-Eξ|<ε}≥1-Dξ
ε2
3.1.6 大数定律及中心极限定理的一般形式
1.大数定律
设{ξn}是随机变量序列,若存在常数序列{an}对任给的ε>0,有
lim n→∞
P 1 nΣn
k=1
ξk -an <ε =1
则称随机变量序列{ξn}服从大数定律。
3个常用的大数定律形式如下所述。
1)伯努利(Bernoulli)大数定律
n 重独立试验中事件A 出现的频率νnn
依概率收敛于事件A 在每次试验中出现的概率p
(0<p<1),即对任给的ε>0,有lim n→∞
P νnn
-p <ε =1。
72
伯努利大数定律以数学形式表达了随机事件出现的频率的稳定性,即当n 很大时,事件
A 发生的频率以很大的概率趋向于每次试验中出现的概率,即νnn
→p(A)=p。
2)切比雪夫大数定律
设{ξn}是相互独立的随机变量序列,又有Dξn ≤C,C 为常数,n =1,2,…,则对任意的
ε >0,有lim n→∞
P 1 nΣn
k=1
ξk -1 nΣn
k=1Eξk <ε =1。
3)辛钦(Χинчин)大数定律
设{ξn}是相互独立同分布的随机变量序列,且有有限的数学期望Eξn =μ,n=1,2,…,则
对任意的ε>0,有
lim n→∞
P 1 nΣn
k=1
ξk -μ <ε =1 即 1 nΣn
k=1
ξk
p→μ
辛钦大数定律可理解为:对一个随机变量ξ在进行n 次重复独立观测时,得到的算术平均
1 nΣn
k=1
ξk 在n 很大时会以很大的概率接近一个常数,即接近数学期望值。
2.中心极限定理
设相互独立的随机变量序列{ξn}有有限的数学期望Eξn 与方差Dξn >0(n=1,2,…),并
设ηn =Σn
k=1
ξk - Σn
k=1Eξk
Σn
k=1Dξk
,则ηn 称为{ξn}前n 项和的标准化(规范化),中心极限定理就是要寻
找{ξn}满足什么条件时,ηn 的渐近分布是标准正态分布,即
lim n→∞
P{ηn <x}=∫x
-∞
1
2πe-t2
2dt (3.1.20)
3.1.7 4 个常见的中心极限定理
1.勒维-林德伯格(Lévy-Lindeberg)中心极限定理
设{ξn }是相互独立、同分布的随机变量序列,且有有限的数学期望与方差Eξn =μ,
Dξn =σ2>0(n=1,2,…),则对任意的实数x,有
lim n→∞
P
ì
.
í .. ..
Σn
k=1
ξk -nμ
σ/n <x
ü
t
y ..
..
=∫x
-∞
1
2πe-t2
2dt (3.1.21)
或
Σn
k=1
ξk -E.
è .
Σn
k=1
ξk
.
. ÷
D.
è .
Σn
k=1
ξk
.
. ÷
L→N (0,1), n → ∞ (3.1.22)
或
1 nΣn
k=1
ξk -μ
σ/ n
L→N (0,1), n → ∞ (3.1.23)
73
这个定理表明:当n 充分大时,1 nΣn
k=1
ξk 近似服从正态分布N μ,σ2
n
.
è .
.
. ÷
。
2.棣莫弗-拉普拉斯(DeMoivre-Laplace)中心极限定理
设{ξn}是相互独立同分布的随机变量序列,且ξn ~B (1,p)(0<p<1,n=1,2,…),则对
任意实数x,有
lim n→∞
P
ì
.
í ..
..
Σn
k=1
ξk -np
np(1-p)<x
ü
t
y ..
..
=∫x
-∞
1
2πe-t2
2dt (3.1.24)
式(3.1.24)常常等价地表达如下:设νn 表示n 重独立试验中事件A 出现的次数,
p(0<p<1)是事件A 在每次试验中出现的概率,则对任意区间[a,b),有
lim n→∞
P a ≤ Vn -np
np(1-p)<b =∫b
a
1
2πe-t2
2dt (3.1.25)
此处,B(n,p)为二项分布,Eξ=np,Dξ=np(1-p)。
3.李雅普诺夫中心极限定理
设{ξn}是相互独立的随机变量序列,如果对某个δ >0有0<E|ξk -Eξk|2+δ < ∞,k=
1,2,…,且满足条件lim n→∞
1 B2+δ
n Σn
k=1E|ξk -Eξk|2+δ =0,其中Bn = Σn
k=1Dξk ,则对任意实数x,有
lim n→∞
P
ì
.
í ..
..
Σn
k=1
ξk - Σn
k=1Eξk
Bn <x
ü
t
y ..
..
=∫x
-∞
1
2πe-t2
2dt (3.1.26)
4.林德伯格中心极限定理
设{ξn}是独立随机变量序列,且有有限的方差Dξn >0,n=1,2,…,如果对任意实数τ>0
满足条件
lim n→∞
1 B2n
Σn
k=∫1 |x-Eξk|>τBn
|x -Eξk|2dFk(x)=0
式中,Bn = Σn
k=1Dξn ,Fk(x)是ξk 的分布函数,则对任意函数x ,有
lim n→∞
P
ì
.
í ..
..
Σn
k=1
ξk - Σn
k=1Eξk
Bn <x
ü
t
y ..
..
=∫x
-∞
1
2πe-t2
2dt (3.1.27)
此定理可解释为:若一个随机变量是由大量相互独立的随机因素的影响所造成的,而每一个
别因素在总影响中所起的作用都不是很大,则这种随机变量通常都服从或近似服从正态分布。
大数法则和中心极限定理是蒙特卡罗法的基础,应该通过各种形式的表述和应用来加深
理解。
3.1.8 几种常见的概率模型和分布
常见的概率模型(以下简称为概型)和分布有以下15种。
1.伯努利概型———二项分布
F(x)=P(μ <x)=ΣK <x
CK
npKqn-K , q =1-p (3.1.28)
74
2.泊松(Poisson)分布
F(x)=ΣK <x
(λt)K
K ! e-λt (3.1.29)
式中,λ 为常数,在电信科学中泊松分布可以用来表示单位时间内K 次呼叫的占线率。
3.均匀分布
令Beta分布中a=b=1,即得均匀分布
F(x)=
0, x <a
x -a
b-a , a ≤x ≤b
1, x >b
ì
.
í
..
.
..
(3.1.30)
4.正态分布
F(x)= 1
2πσ∫x
-∞e-(t-μ)2
2σ2 dt (3.1.31)
5.指数分布
F(x)= 0, x <0
1-e-λx , x ≤0 (3.1.32)
6.Gamma分布
Gamma分布是用如下概率密度表示的概率分布,记为Ga(α,λ)。
p(x;α,λ)= λα
Γ(α)xα-1e-λx (3.1.33)
7.Beta分布
Beta分布是用如下概率密度表示的概率分布,记为Be(a,b),其中a、b 是两个正值参数。
p(a;a,b)= Γ(a +b)
Γ(a)Γ(b)xa-1(1-x)b-1 (3.1.34)
8.t 分布
t 分布是用如下概率密度表示的概率分布,记为t(α),其中α 是正实数。
p(x,α)= Γα +1
2
.
è .
.
. ÷
απΓ α2
.
è .
.
. ÷
1+x2
α
.
è .
.
. ÷
-(α+1)
2 (3.1.35)
9.z 分布
z 分布是用如下概率密度表示的概率分布,记为z(a,b),其中a、b 是两个正值参数。
p(x;a,b)= Γ(a +b)
Γ(a)Γ(b)
xa-1
(1+x)a+b (3.1.36)
10.χ 2分布
在Gamma分布中,令α=n/2,λ=1/2,即得自由度为n 的χ2 分布,记为χ2(n)。
p(x;α,λ)= 1
2n/2Γ(n/2)xn2
-1e-x2
, x >0 (3.1.37)
11.指数分布
在Gamma分布中,令α=1,即得指数分布,记为exp(λ)。
p(x;λ)=λe-λx , x >0 (3.1.38)
75
12.反余弦分布
在Beta分布中,令a=b=1/2,即得反余弦分布。
13.多项分布
P(X1 =n1,…,Xr =nr)= n!
n1!…nr !pn11 …prnr (3.1.39)
14.非中心Gamma分布
非中心Gamma分布是如下概率密度表示的概率分布,记为Ga(α,λ,γ)。
p(x;α,λ,γ)=e-γγm
m ! PG (x;α +m ,λ) (3.1.40)
15.非中心t 分布
非中心t 分布是用如下概率密度表示的概率分布,记为t(n,γ),其中γ 是非负实数。
p(x,n,γ)= nn 2e-r2
2
πΓ n2
.
è .
.
. ÷
(n +x2)n+1
2 Σ∞
m =0Γm +n +1
2
.
è .
.
. ÷
λm
m !!
2x
n +x2
.
è ..
.
. ÷÷
m
(3.1.41)
3.1.9 蒙特卡罗法简单应用举例
为了有一个直观的基本概念,举例如下。
例3.1 Buffon投针试验求π的近似值(1777年)。
平面上有相距2a 的两条平行线,向这一平面随机投放一个长为2l 的针,设a>l>0,针
的中点用M 表示,y 是中点M 与最近一条平行线的距离,. 是针与平行线的交角。显然有如
下的条件(参见图3.2)。
(1)0≤y≤a,0≤.≤π;
(2)针与平行线相交的充要条件为y≤lsin.,于是,针与平行线相交概率为
p = 1
πa∫π
∫0 lsin.
0 dyd. = 1
πa∫π
0
lsin.d. =2l
πa (3.1.42)
此处,y 的变化范围是由0到a ,. 则只能由0变到π,. 和y 变化范围扫过的总面积为πa,总的
概率为1,参见图3.3,则所有可能发生的事件的总值为πa,而投针与平行线相交的概率面积为
∫π
∫0 lsin.
0 dyd.。用N 表示投针的次数,ν 是针与平行线相交的次数。当试验次数足够多时,相交
概率为p =ν/N ,由式(3.1.42)有ν
N ≈ 2l
πa,于是有π≈2lN
aν 的关系。
图3.2 Buffon投针试验示意图图3.3 投针试验中针与线相交概率
Buffon投针试验可以编计算机程序进行计算。Wolf在1853年取a=45,l=36,N =
76
5000进行试验时,得到v=2532,计算出的π值为3.1596。表3.2列出一些试验结果。
表3.2 历史上投针试验获得的π的近似值
试验者时 间针长(l/a) 投针次数相交次数π的估计值
Wolf 1850年0.80 5000 2532 3.15956
Smith 1855年0.60 3204 1218 3.15665
Fox 1884年0.75 1030 489 3.15951
Lazzarini 1925年0.83 3408 1808 3.14159292
实际上,该试验的精度并不高。设l=a,则p=2/π,由中心极限定理可知:当试验次数为
图3.4 随机投点求积分值
n 时,如果要求置信度为0.95,同时要求π的估计值达到小
数点后三位有效数字,则试验次数n 应该达到8.87×105,显
然表3.2中Lazzarini的结果纯属巧合(该结果接近祖冲之的
密率355/113)。
例3.2 见图3.4,求定积分值I =∫1
0
f(x)dx。
设f(x)是[0,1]上的连续函数且0≤f(x)≤1,那么当
向正方形[1,1]内均匀投点(ξ,η)时,随机点落入曲线
y=f(x)下面的概率就是要求的积分值。
Pr {y ≤f(x)} =∫1
∫0 f(x)
0 dydx =∫1
0
f(x)dx =I
例3.3 应用泊松过程估算地面移动设备雷击概率。
(1)泊松分布与二项分布。
泊松分布源于二项分布,因为泊松分布是某段时间间隔或某个给定区域内发生的结果是
否达到了某个数量,本质上还是达到与没达到的选择。二项分布为
P {y =k} =
n
k
.
è .
.
.÷pkqn-k p +q =1 k =0,1,2,…,n
实际情况中,当p 或q 很小,而试验次数n 很大时,n→∞,p→0,两者是同级的无限大/小,
所以np→λ,一定趋于一个常数λ。设试验时间无限长,每次试验的间隔时间为D,在区间D中
发生k 次呼叫的概率满足如下的二项分布:
pn(k)= n!
(n -k)!k! pk(1-p)n-k k =0,1,2,…,n
泊松定理可推导出二项分布的极限分布。
泊松定理如下:如果n→∞p→0,np→λ,则有
n!
(n -k)!k! pk(1-p)n-k → n→∞e-λ λk
k! k =0,1,2,…
(2)首次发生的小概率事件的泊松分布。
把泊松定理应用于p 很小,而试验次数n 很大时的二项分布问题中,也就是用于随机泊
松点的问题:随机放置n 个点到(-T/2,T/2)区间并且用P {k 在ta 内}表示这些点中的k
个将落在(t1,t2)内的概率。设问题要考虑的区间长度为ta =t2-t1,p=ta/T ,则在n.1并
且ta =t2-t1.T 时,有
77
P {k 在ta 内} ≈e-
nta
T
nta
T
.
è .
.
. ÷
k
k! =e-λta (λta ) k
k!
用X 表示在一个泊松试验中得到的某个随机变量结果的数量,称为泊松随机变量,在公
式中表示为X >x,其概率分布为泊松分布,结果的均值为μ=λt,t 是具体的时间、距离、面积
或体积;λ 是结果的发生率,表示在单位时间、长度、面积或体积内得到结果的平均数量。某段
时间间隔或某个给定区域内某事件发生的次数满足泊松分布:
p (x;λt) =e-λt (λt)x
x! x =0,1,2,…
其中,X (x)表示在给定的时间间隔或指定的区域t 内结果的发生数量;X 为泊松随机变量;λ
是在单位时间、长度、面积或体积内得到结果的平均数量。
例如,一个试验中,1μs内通过计数器的平均辐射粒子数为4个。在某个1μs中有6个辐
射粒子通过计数器的概率是多少?
按所给条件知道,x=6,μ=λt=4,有
p(6;4)=e-446
6! =Σ6
x=0
p(x;4)- Σ5
x=0
p(x;4)=0.8893-0.7851=0.1042
又例如,在12小时内,有180个电话随机打来,问4小时中打入电话数介于50~70个的
概率有多大? 此时p=4/12,x=180,k=50~70。
泊松分布为离散分布,用来计算在某时间间隔或空间范围内一定数目的泊松事件发生的
概率,就变成了连续分布的泊松问题。在很多应用中,时间或空间可以是随机变量,而泊松分
布是一个单参数的分布,参数λ 可以理解为单位时间内事件发生的平均数目。
在许多情况下,更关心的是问题中的事件发生没发生,可以等价地考虑首次发生该事件的
时间是否落在我们关心的时间或区间段内。当我们考虑描述事件首次发生所需要的等待时间
的随机变量时,可以研究在时间t 时,事件还没有发生时的概率,于是互补的概率问题是只要
发生事件的概率。于是,有
(首次发生事件所需要的等待时间的概率)=1- (只要发生事件的概率)
泊松事件首达时刻,即在某段时间间隔或某个给定区域内还没有发生概率事件的概率等
价于发生所有事件的概率,即
P (X ≥1) ==1-p (0;λt) =1-e-λt (λt) 0
0! =1-e-λt =e-λt
(3)根据问题条件构建概率密度函数。
为了构造移动设备雷击概率函数,首先要回顾一下构建概率函数的理论。针对随便选取
或构造的一些函数,怎么判断它是否为概率函数? 从数学上讲,应该满足随机函数的存在性
定理。
存在性定理:给定一个函数f x( ) 或它的积分F x( ) =∫x
-∞
f u( )du,如果要求能够构造一
个随机变量x 和一个随机试验,具有分布密度函数f x( ) 或分布函数F x( ) ,则要求该函数必
须具有如下性质:函数f (x ) 或分布函数F (x ) 必须是非负的并且在(- ∞,+ ∞) 上的积分
F x( ) =1;F x( ) 必须是右连续的,并且当x 从- ∞ 增加到∞ 时,函数f x( ) 必须单调地从0
增加到1。
在构建关于随机变量的概率分布函数时,其均值是最重要的信息,因为这是无限次概率试
78
验得到的必然性的结果。对于泊松概率分布来讲,其均值为μ=λt,
没有事件发生时x=0,泊松分布变成指数分布,则x 的累积分布为
P =1-e-λx =e-λx
指数分布函数:如果随机变量x 的概率密度函数为fx x( ) = λe-λx , x≥0
0, 其他{ ,则称随机变
量x 是服从具有参数λ 的指数分布。如果在互不相交的区间上事件发生是相互独立的,例如
电话呼叫的到达时间、公共汽车到达一个车站的时间或类似的事件的等待时间,都可以用指数
分布来描述。
因此,在构建移动设备雷击概率函数时,首先要合理地假设其数学期望值或均值。在我们
的问题中,讨论的是近地面发生的雷击事件,尽管描述车辆、公路、雷云的参数复杂各异,但是
有意义的最主要的是地面物体的引雷面积。因此,取地面移动目标的等效静止区域面积与雷
云区行车路线决定的公路的面积之比作为发生雷击概率的均值μ=λt。定义概率pem 主要由
处于雷云区的公路的总的引雷面积和移动目标的等效静止引雷面积决定。
设移动目标的等效静止区域面积为σem;雷云区的行车路线的公路总面积的下确界,也
就是公路的总建筑面积,仍然用σcm 表示;用σ'cm表示由移动出发点到发生雷击点的公路的
总建筑面积;公路的总长度用Lcm表示;由移动出发点到发生雷击点的公路的总长度用L'cm
表示;移动目标的等效静止引雷面积用σem表示,而移动目标的等效静止引雷面的长度用则
Lem表示。
从实际情况考虑,公路的宽度与移动车辆的宽度有一定的比例关系,例如两车道、四车道、
六车道等。假设最简单公路情况,也就是两车道的情况,就可以把公路宽度D 设为移动车辆
的2倍。根据地面物体引雷面与高度的关系,考虑车辆比地面有一定的高度,所以可以忽略这
个差异,则有
σem =Lem ×d车宽.Lem ×D (3.1.43)
σcm =Lcm ×D (3.1.44)
在我们的问题中,λ 是事件的发生率或均值,可以认为λ 等于移动目标的等效静止区域面积与
雷云区的行车路线的公路总面积之比,也就是公路的总长度与移动目标的等效静止引雷长度
之比,即λ=Lem
Lcm 。
问题中,首次发生雷击事件的位置是不确定的随机数,就是近地面发生的雷击地面移动设
备事件中的随机变量x。不管雷击发生在何处,都属于地面移动目标的“首次发生雷击”这个
事件,设L'cm是“首次发生雷击”的随机位置,就可以把随机变量x 用x=Lem
L'cm
表示。
经过推导得到
pem = Lem
Lcm
.
è .
.
. ÷
e-
Lem Lcmx =λe
-
Lem
Lcm 1-
Lem L'cm
( ) →λe-λ 1
1-x =λe-λ (x+x2+x3+… ) ≈λe-λx
在这个概率表达式中,选取e-
Lem
Lcm 1-Lem ( L'cm ) 是为了使概率函数满足存在性定理的要求:F(x)
必须是右连续的,并且当x 从-∞增加到∞时,函数f x( ) 必须单调地从0增加到1。
当L'cm=Lem时,由于公路非常短,不可能发生雷击,分母为0,则e-∞ =0;当Lcm为无穷长
时,表示公路非常长,其上移动的车辆一定会遭受雷击,也就是e-
Lem Lcm →e-Lem
∞ =e0=1。综上所
79
述,我们设计的首次发生雷击事件概率函数符合存在性定理的要求条件。
设λ=Lem
Lcm ,而首次发生雷击事件的位置L'cm是不确定的随机数,但是不管发生在何处,都
属于地面移动目标的“首次发生雷击”这个事件,设x=Lem
L'cm ,有
e
-
Lem
Lcm 1-
Lem L'cm
( ) →e-
Lem Lcm
( ) 1
1-x =e-
Lem Lcm
( ) (x+x2+x3+… ) →e-
Lem Lcm
( ) x =e-λx
以典型数值为例,设地面移动目标长度的线度L移动目标为10~20m,再加上地面移动目标
运动产生的等效长度40~60m,设移动目标等效静止引雷长度Lem 为50~80m。选取行驶的
公路长度为1km,按此数据计算可得到pem为0.0475~0.0784,pem随着等效静止引雷长度和
公路长度取值的不同而变化。
根据上述演算,大型运输车运输中被雷击是个小概率事件。但是,这个举例只是配合教学
的简单情况,实际情况中,雷云高度、厚度,云-地温度差决定的雷云带电量,以及运输货物的
易燃易爆性质等因素都对结果产生很大影响,并不能得出大型运输车运输中被雷击是小概率
事件的结论。
3.2 伪随机数
3.2.1 简单子样
产生随机数是执行蒙特卡罗法的基本条件,通常用随机数序列的方法实现。具有单位均
匀分布的总体中所产生的简单子样称为随机数序列,其中每一个体称为随机数。
简单子样指的是这样的随机子样{xn}N
n=1={x1,…,xN },它们之间相互独立,对于任意的
x 和所有的n=1,2,…,N ,满足
P(x <xn <x +dx)=f(x)dx (3.2.1)
简单子样的性质如下。
(1)对于任意的统计量g(x),{g(xn)}N
n=1具有同一分布且相互独立,数学期望均等于
G =E(g)=∫g(x)dF(x)=∫g(x)f(x)dx (3.2.2)
F(x)与f(x)分别为随机变量ξ 的分布函数和分布密度函数。
(2)只要E(|g|)=∫|g(x)|f(x)dx 存在,则蒙特卡罗近似估计值S
N =1 NΣN
n=1
g(xn)
便依概率1收敛于G。
(3)只要统计量g(x)的方差σ2 存在,且异于0,则N| S
N -G|/σ 就渐进于正态分布。
3.2.2 随机数与伪随机数
前面提到,由具有已知分布的总体中产生的简单子样在蒙特卡罗法的问题中占有非常重
要的地位。形式多样的简单子样可以按照一定的规则利用基础的随机变量产生,显然最合适
的基础随机变量是单位均匀分布的随机变量,就如同可以用黑点在白纸上构造出任意形状的
图案的情形一样。在蒙特卡罗法中,把这种分布的简单子样作为基础量来看待,而针对其他分
布产生简单子样时,就把这种具有单位均匀分布的随机变量看成已知条件。通常人们采用随
80
机数作为计算机上的单位均匀分布的随机变量。因此,生产随机数成为应用蒙特卡罗法的大
前提。
1.随机数
均匀分布也称为矩形分布,其中最基本的是单位均匀分布,其密度函数为
f(x)=1(0≤x≤1),而从具有单位均匀分布的总体中产生的简单子样x1,x2,…,xN 称为随
机数序列,其中的每一个体称为随机数,并且用ξ1,ξ2,…,ξN 等符号表示。
2.随机数的性质
设随机数ξ1,ξ2,…是相互独立的具有相同单位均匀分布的随机变量序列,则对任意自然
数S,由S 个随机数所组成的S 维空间上的点(ξn+1,ξn+2,…,ξn+s)也在S 维空间的单位方体
GS 上均匀分布,也就是对于任意小的ai,0≤ai≤1,i=1,2,…,s,有下式成立:
P(ξn+i ≤ai,i=1,2,…,s)=Πs
i=1
ai (3.2.3)
客观上存在许多真正的随机变量,因而用物理的方法产生随机数应当不成问题。例如,可
以用放射性同位素的衰变来产生随机数。1942年,RAND 公司就是以随机脉冲源,用电子旋
转轮产生随机数表的,但要产生又快又准的随机数是很困难的。例如,1975—1978年,Frigerio
做了一个高分率的计数器,每小时能产生6000个31位的真随机数。但是,随机过程一去不复
返,而且设备的维修费用昂贵,显然产生的随机数的数量少而且慢。
3.伪随机数
由于采用真随机数遇到困难,人们普遍用数学方法产生随机数,但这种随机数是由递推公
式和初值推出的,也就不满足随机数相互独立的要求。由于计算机所能表示的[0,1]上的数量
是有限多个,而随机数又由递推公式产生,所以就会出现周期现象。因此,一般称用数学方法
产生的随机数为伪随机数。
尽管存在这些问题,但在一定范围的应用中,只要这些伪随机数序列能通过一系列的统计
检验,如均匀性和独立性的检验,就可以把它当作随机数来使用。伪随机数的应用范围与伪随
机数的总体大小有直接关系,在应用中伪随机数的数量应能满足所设概率模型和问题精度的
要求。实际上,随机数无边无际变化的性质总是体现在有限逼近的过程中。
4.伪随机数的最大容量
由伪随机数组成的子样有一个最大限度,超出这个限度将出现周期循环现象,则从伪随机
数的初始值开始到出现循环现象为止所产生的伪随机数的个数称为伪随机数的最大容量。
3.2.3 产生伪随机数的几种方法
1.平方取中法
这是最早的一种方法,由冯·诺依曼提出,是历史上第一个产生伪随机数的方法。
以十进制为例,平方取中法是把一个2s 位的十进制数自乘,去头截尾保留中间的2s 位数
字,再用102s来除使结果落在[0,1]区间,得到的结果就是[0,1]上的伪随机数,平方取中法可
用公式表示为
十进制:
Xn+1≡[10-sX2n](mod102s)
ξn+1=Xn+1
102s (3.2.4)
81
二进制:
Xn+1≡[2-sX2n](mod22s)
ξn+1=Xn+1
22s (3.2.5)
式中,用10-s乘X2n
表示“截尾”,(mod102s)运算表示“去头”,A mod(N )表示取余运算,即A
除以N 所剩的余数,[*]符号(方括号)表示取整运算,取不超过“*”的最大整数。
与此类似,Forsythe提出类似的方法如下:
Xn+2≡[2-sXnXn+1](mod22s)
ξn+2=Xn+2
22s (3.2.6)
平方取中法的最大容量随初值的不同而变化,情形非常复杂。
练习题:取十进制数,用平方取中法产生伪随机数序列并比较序列长度,找到周期发生部
位(建议数字:2s=4,X0=1111和X0=6400)。
当用数学方法产生伪随机数时,如果计算的伪随机数产生了零值或产生了重复的元素,则
伪随机数序列中断,并决定了此随机数最大序列长度。平方取中法缺点很多,现已很少采用,
目前使用比较多的是乘同余法和乘加同余法。
2.加同余法
任取两个数作为初值,用下面的公式产生后续随机数:
Xn+2≡Xn +Xn+1(modM )
ξn+2=Xn+2
M
(3.2.7)
当采用X1=X2=1时,所产生的序列就是著名的斐波那契(Fibonacci)数列,该数列在生
物等多种学科中有广泛的应用。一般性地讨论加同余法的最大容量很困难,但是当取
M =2s 时,该随机序列的最大容量为3×2s-1。
3.乘同余法
乘同余法是Lehmer首先提出的。其方法为:对任意取定的初值再用式(3.2.8)决定后续
随机数。
Xn+1≡aXn(modM )
ξn+1=Xn+1
M
ì
.
í ..
..
(3.2.8)
式中,a 为选定的常数,在式中最好取X1 为4q+1型的数,q 为任意整数,a 取52K +1型的数,
K 是使52K +1达到计算机允许的最大奇数时所确定的数。
其最大容量主要取决于产生周期的时刻,与a 的取值有关,通常与初值X1 的选取无关,
与M 的具体值关系不大,但一般要取X1 与M 为互素数。
4.乘加同余法
乘加同余法由Rotenbery提出,其方法如下:对任意初值X1,由下式确定伪随机数。
Xn+1≡aXn +c (modM )
ξn+1=Xn+1
M , a,c 为常数(3.2.9)
乘子a 为小于模M 的正整数,c 为非负整数,常选为0。当c=0、M 取素数、a 是M 的一
个原根时用乘同余法产生的随机数的周期为M -1,而M 通常取为M =231-1,a=16087或
82
630360016。当取M =2K 、K 为整数、a 为模8余3或5的整数时所产生的随机数周期
为2K -2。
1968年,Marsaglia从理论上指出乘同余法的缺点是:给定长度K ,设(ξi1,ξi2,…,ξik )为
乘同余法所产生的随机数空间上的一点,这些点的分布具有栅格结构,也就是,这些点落在具
有确定数量的等距离的平行超平面上,也就是说是有规律的,破坏了应有的随机性。
5.移位寄存器法———Tausworthe法
一般形式为
βj =β0xTj
ξj =Σn
k=1
b(j-1)×n+kx2k-1, j=1,2,… { (3.2.10)
式中,bj∈Z(Z 表示只含数0和1的二元有限域),βj=(b1,j,b2,j,…,bn,j)是1×n 的二进制
矩阵。
T =(I+Rs)(I+Lt)
I=
1 0 0 … 0
0 1 0 … 0
0 0 1 … 0
. . . . .
0 0 0 … 1
é
.
êêêêêêê
ù
.
úúúúúúú
n×n
R =
0 1 0 … 0
0 0 1 … 0
. . . . .
0 0 0 … 1
0 0 0 … 0
é
.
êêêêêêê
ù
.
úúúúúúú
n×n
L =
0 0 … 0 0
1 0 … 0 0
0 1 … 0 0
. . . . .
0 0 … 1 0
é
.
êêêêêêê
ù
.
úúúúúúú
n×n
t 为整数,I、R、L 为n×n 矩阵。
移位寄存器法产生的伪随机数周期为2n-1,但这种伪随机数发生器可能具有严重的高位
序列相关性。
6.斐波那契法
该法的一般形式为Xi=(Xi-p ..Xi-q)mod2m ,其中..是一种运算操作,可以定义为加、
减、乘、逻辑异或(p、q 为延迟步,且p>q)。
斐波那契法最大的优点是可以直接进行浮点运算,而不像其他算法需要经整数转换。
7.混合法
混合法于1965年由Donald和Marsaglia提出,该方法将两个或三个乘同余方法简单组
合。步骤如下。
(1)给出两个单独的乘同余发生器:
{Xi =(αXi-1)modMx
Yi =(βYi-1)modMy
(2)组合规则:
Zi =Xi ..Yi ={Xi -Yi +Mx -1, Xi -Yi ≤1
Xi -Yi, 其他
ξi =Zi/Mx , i=1,2,3,…
ì
.
í
..
..
83
如果取α=40014,β=40692,Mx =2147483563,My =2147483399,产生的随机数周
期为(Mx -1)(My -1)/2。
8.复杂组合法
复杂组合法于1990年由Marsaglia、Zaman和Tsang设计。
基本思想:采用斐波那契序列和乘同余方法设置初值,再利用一个简单的算术序列和斐
波那契序列给出随机数。步骤如下。
(1)产生初始种子序列{Yn}{Zn}{bn}:
Yn =Yn-3*Yn-2*Yn-1 mod179
Zn =53*In-1 mod169
bn ={0, Yn *Z0 mod64<32
1, 其他
初值Y0,Y1,Y2 为[1,178]区间内的任意整数,不能同时全为1;初值Z0为[0,168]区间的
任意整数。由此产生M ×N 个bi,即i =1,2,…,M ×N 。M 表示所需的种子数,N 为二进
制小数的位数。把种子转换为实型数:
Xi =ΣN
j=1
b(i-1)*N+j*2-j, i=1,2,…,M
(2)产生简单算术序列{Cn}:
Cn =Cn-1 ..d ={Cn-1 -d, Cn-1 ≥d
Cn-1 -d +e, 其他
式中,d 等于任意整数除以2N 所得的数,e 等于对小于2N 的最大素数除以2N 所得的数,C0等
于任意整数除以2N 所得的数,序列{Cn}的周期为小于2N 的最大素数。
(3)斐波那契序列及组合规则。
斐波那契序列为
Xn =Xn-ip ..Xn-jp
定义运算..为
X ..Y ={x -y, x ≥y
x -y +1, 其他
(4)组合产生随机数序列{Un}:
Un =Xn ..Cn
这种方法产生的伪随机数序列有很长的周期且能通过各种统计检验,如果对(3)的斐波那
契序列稍做改进即可得到周期更长的伪随机数序列。
3.2.4 伪随机数的检验
只有经过一系列检验,所产生的伪随机数序列才能作为随机数使用。重要的检验内容简
述如下。
1.均匀性检验
一般伪随机数序列可表示为{ξn}N
n=1,均匀偏度定义为
Δ(N )=0s≤uxp≤1
1 NΣN
n=1
η(ξn <X )-X (3.2.11)