第12章 其他建模方法
12.1 计算机模拟
蒙特卡罗（Monte Carlo）方法，也称为计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹的“曼哈顿计划”。该计划的主持人之一，数学家冯·诺依曼用驰名世界的赌城—摩纳哥的蒙特卡罗来命名这种方法，使它蒙上了一层神秘的色彩。
早在17世纪，人们就已经知道用事件发生的“频率”作为事件的“概率”的近似值。只要设计一个随机试验，使一个事件的概率与某未知数有关，然后通过重复试验，以频率近似表示概率，即可求得该未知数的近似解。显然，利用随机试验求近似解，试验次数要相当多才行。随着20世纪40年代电子计算机的出现，人们便开始利用计算机来模拟所设计的随机试验，使得这种方法得到迅速的发展和广泛的应用。
12.1.1 随机变量的模拟
利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟。
1.离散随机变量的模拟
设随机变量X的分布律为P{X?xi}?pi,i?1,2,L，令
p(0)?0,p(i)??pj,i?1,2,L，
j?1i将p(i)作为分点，将区间(0,1)分为一系列小区间(p(i?1),p(i))(i?1,2,L)。对于均匀分布的随机变量R~U(0,1)，则有
P{p(i?1)?R?p(i)}?p(i)?p(i?1)?pi,i?1,2,L，
由此可知，事件{p(i?1)?R?p(i)}和事件{X?xi}有相同的发生概率。因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况。
例12.1 已知在一次随机试验中，事件A,B,C发生的概率分别为0.4，0.5，0.1。模拟1000次随机试验，计算事件A,B,C发生的频率。
在一次随机试验中，事件发生的概率分布见表12.1。
表12.1 事件发生的概率分布
事件 概率 累积概率 A 0.4 0.4 B 0.5 0.9 C 0.1 1 我们用产生[0,1]区间上均匀分布的随机数，来模拟事件A,B,C的发生。由表12.1的数据和几何概率的知识，可以认为如果产生的随机数在区间[0,0.4]上，事件A发生了；产生的随机数在区间(0.4,0.9]上，事件B发生了；产生的随机数在区间(0.9,1]上，事件C发生了。产生1000个[0,1]区间上均匀分布的随机数，统计随机数落在相应区间上的次数，就是在这1000次模拟中事件A,B,C发生的次数，再除以总的试验次数1000，即得到事件A,B,C发生的频率。
293
模拟的MATLAB程序如下：
clc, clear, n=1000;
a=rand(1,n); %产生n个[0,1]区间上的随机数
n1=sum(a<=0.4), n2=sum(a>0.4 & a<=0.9), n3=sum(a>=0.9) f=[n1,n2,n3]/n %计算各事件发生的频率
例12.2 事件Ai(i?1,2,L,10)发生的概率如表12.2所示，求在10000次模拟中，事件Ai(i?1,2,L,10)发生的频数。
表12.2 事件Ai发生的概率分布
事件 概率 A1 0.2 A2 0.05 A3 0.01 A4 0.06 A5 0.08 A6 0.1 A7 0.3 A8 0.05 A9 0.03 A10 0.12 模拟的MATLAB程序如下：
clc, clear, n=10000;
a=[0.2,0.05,0.01,0.06,0.08,0.1,0.3,0.05,0.03,0.12]; b=cumsum(a) %求累加和
c=[0,b(1:end-1)]; %所有小区间的左端点 m=zeros(1,10); %计数器初始化 d=rand(1,n); %产生n个随机数
for i=1:n
ind=find(d(i)>c,1,'last'); %从后往前找满足d(i)大于c中第一个数的地址 m(ind)=m(ind)+1; end
m %显示频数数据 2.连续型随机变量的模拟
利用[0,1]区间上的均匀分布随机数可以产生具有给定分布的随机变量数列。
我们知道，若随机变量?的概率密度函数和分布函数分布为f(x),F(x)，则随机变量
??F(?)的分布就是区间[0,1]上的均匀分布。因此，若Ri是[0,1]中均匀分布的随机数，那么方程
?xi??f(x)dx?Ri
(12.1)
的解xi就是所求的具有概率密度函数为f(x)的随机抽样。这可简单解释如下。
若某个连续型随机变量?的分布函数为
F(x)??x??f(x)dx，
不失一般性，设F(x)是严格单调增函数，存在反函数x?F?1(y)，下面我们证明随机变量
??F(?)服从[0,1]上的均匀分布，记?的分布函数为G(y)，由于F(x)是分布函数，它的取值在[0,1]上，从而当0?y?1时
G(y)?P{??y}?P{F(?)?y}?P{??F?1(y)}?F(F?1(y))?y，
因而?的分布函数为
294
?0,y?0,?G(y)??y,0?y?1,
?1,y?1.??服从[0,1]上的均匀分布。
R为[0,1]区间上均匀分布的随机变量，则随机变量??F?1(R)的分布函数为F(x)，概率
密度函数为f(x)，这里F?1(x)是F(x)的反函数。
所以，只要分布函数F(x)的反函数F?1(x)存在，由[0,1]区间上均匀分布的随机数Rt，求xt?F?1(Rt)，即解方程
F(xt)?Rt，
就可得到分布函数为F(x)的随机抽样xt。
例12.3 求具有指数分布
??e??x,x?0,f(x)??
0,x?0.?的随机抽样。
设Ri是[0,1]区间上均匀分布的随机数，利用(12.1)式，得
Ri??xi??xif(x)dx???e??xdx?1?e??xi.
0所以
1xi??ln(1?Ri)
?就是所求的随机抽样。
由于1?Ri也服从均匀分布，所以上式又可简化为
1xi??lnRi.
? MATLAB工具箱中常用分布的随机数都可以产生，实际上不需要我们去生成各种分布
的随机数。
12.1.2 随机模拟的例子
例12.4 设计随机试验求?的近似值。
在单位正方形中取1000000个随机点(xi,yi)，i?1,2,L,1000000，统计点落在x2?y2?1内的频数n。则由几何概率知，任取单位正方形内一点，落在单位圆内部（第一象限部分）的概率为p??n?4n，由于试验次数充分多，频率近似于概率，有所以??。 ?，4100000041000000y1O295
1x
数学建模简明教程课件-第12章-其他建模方法
