第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|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 (00,有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 0(n=1,2,…),则对任意的实数x,有 lim n→∞ P ì . í .. .. Σn k=1 ξk -nμ σ/n 0有00,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 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 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