蒙特卡洛方法初探
实验目 的
实验内容
学习计算机模拟的基本过程与方法。
1 、模拟的概念。
4 、实验作业。
3 、计算机模拟实例。
2 、产生随机数的计算机命令
。
模拟的概念
模拟就是利用物理的、数学的模型来类比、模
仿现实系统及其演变过程，以寻求过程规律的一 种方法。
模拟的基本思想是建立一个试验模型，这个模 型包含所研究系统的主要特点．通过对这个实验 模型的运行，获得所要研究系统的必要信息
模拟的方法
1 、物理模拟：
对实际系统及其过程用功能相似的实物
系统去模仿。例如，军事演习、船艇实验、
沙盘作业等。
物理模拟通常花费较大、周期较长，且
在物理模型上改变系统结构和系数都较困难。
而且，许多系统无法进行物理模拟，如社会
经济系统、生态系统等。
在实际问题中，面对一些带随机因素的复杂系 统，用分析方法建模常常需要作许多简化假设， 与面临的实际问题可能相差甚远，以致解答根本 无法应用。这时，计算机模拟几乎成为唯一的选 择。
在一定的假设条件下，运用数学运算模拟系
统的运行，称为数学模拟。现代的数学模拟都是 在计算机上进行的，称为计算机模拟。
2 、数学模拟
计算机模拟可以反复进行，改变系统的结构
和系数都比较容易。
蒙特卡洛（ Monte Carlo ）方法是一种 应用随机数来进行计算机模拟的方法．此方 法对研究的系统进行随机观察抽样，通过对 样本值的统计分析，求得所研究系统的某些 参数．
用蒙特卡洛方法进行计算机模拟的步骤： [1] 设计一个逻辑框图，即模拟模型．这个框 图要正确反映系统各部分运行时的逻辑关系。 [2] 模拟随机现象．可通过具有各种概率分布 的模拟随机数来模拟随机现象．
产生模拟随机数的计算机命令
在 Matlab 软件中，可以直接产生满足各种分 布的随机数，命令如下：
1 ．产生 m*n 阶 (a ， b) 均匀分布 U （ a ， b ）的随机数矩阵： unifrnd (a,b,m, n); 产 生一个 [a ， b] 均匀分布的随机数： unifrnd (a,b)
当只知道一个随机变量取值在（ a ，
b ）内，但不知道（也没理由假设）它在何处 取值的概率大，在何处取值的概率小，就只 好用 U （ a ， b ）来模拟它。
2 ．产生 mm*nn 阶离散均匀分布的随 机数矩阵：
R = unidrnd(N)
R = unidrnd(N,mm,nn)
•当研究对象视为大量相互独立的随机变量之和， 且其中每一种变量对总和的影响都很小时，可以 认为该对象服从正态分布。
•若连续型随机变量 X 的概率密度函数为
其中 >0 为常数，则称 X 服从参数为 的指数 分布。
0 0
1 0 )
( /
x x x e
f x
•指数分布的期望值为
•排队服务系统中顾客到达间隔、质量与可靠性 中电子元件的寿命通常服从指数分布。
例 顾客到达某商店的间隔时间服从参数为 1 0( 分钟 ) 的指数分布 ( 指数分布的均值为 1
0) --- 指两个顾客到达商店的平均间隔时间
是 10 分钟 . 即平均 10 分钟到达 1 个顾客 . 顾客到达的间隔时间可用 exprnd(10) 模拟
。
•设离散型随机变量 X 的所有可能取值为 0,1,2,…, 且 取各个值的概率为
其中 >0 为常数，则称 X 服从参数为 的帕松分布。 ,
, 2 , 1 , 0
! , )
(
k k k e
X P
k
•帕松分布在排队系统、产品检验、天文、物理等 领域有广泛应用。
•帕松分布的期望值为
6 产生 1 个参数为 n,p 的二项分布的随机数 binorn d(n,p) ，产生 mn 个参数为 n,p 的二项分布的随机 数 binornd(n,p,m,n) 。
掷一枚均匀硬币，正面朝上的次数
X 服从参数为１， p 的二项分布 ,X~B
(1,p)
频率的稳定性模拟
1. 事件的频率
在一组不变的条件下，重复作 n 次试
验，记 m 是 n 次试验中事件 A 发生的次数
。 频率 f=m/n
2. 频率的稳定性
掷一枚均匀硬币，记录掷硬币试验中
频率 P* 的波动情况。
function liti1(p,mm)
pro=zeros(1,mm);
randnum = binornd(1,p,1,mm)
a=0;
for i=1:mm
a=a+randnum(1,i);
pro(i)=a/i;
end
pro=pro
num=1:mm;
plot(num,pro)
在 Matlab 中编辑 .m 文件输入以下命令：
在 Matlab 命令行中输入以下命令：
liti1(0.5,1000)
在 Matlab 命令行中输入以下命令：
liti1(0.5,10000)
练习 掷一枚不均匀硬币，正面出现概率
为 0.3 ，记录前 1000 次掷硬币试验中
正面频率的波动情况，并画图。
在 Matlab 命令行中输入以下命令：
liti1(0.3,1000)
例 2 掷两枚不均匀硬币，每枚正面出现概
率为 0.4 ，记录前 1000 次掷硬币试验中
两枚都为正面频率的波动情况，并画图。
在 Matlab 中编辑 .m 文件输入以下命令：
function liti2(p,mm)
pro=zeros(1,mm);
randnum = binornd(1,p,2,mm) ； a=0;
for i=1:mm
a=a+randnum(1,i)*randnum(2,i); pro(i)=a/i
;
end
pro=pro ， num=1:mm;plot(num,pro)
liti2(0.4,100)
liti2(0.4,10000)
例 3: 在一袋中有 10 个相同的球，分别标有 号码 1,2,…,10 。每次任取一个球，记录其号 码后放回袋中，再任取下一个。这种取法叫 做“有放回抽取”。今有放回抽取 3 个球，求这 3 个球的号码均为偶数的概率。（用频率估计 概率） 解：有放回取 3 个球 , 所有取法有 103 种 ; 有
放回取 3 个偶数号码的球 , 所有取法有 53 种 . 所以
8 1 10
) 5
( 3
3
A
P 古典概率模拟
function proguji=liti3(n,mm)
frq=0;
randnum=unidrnd(n,mm,3);proguji=0;
for i=1:mm a=(randnum(i,1)+1)*(randnum(i,2)+
1)*(randnum(i,3)+1) ；
if mod(a,2)==1
frq=frq+1
end
end;
proguji=frq/mm
例 4 两盒火柴，每盒 20 根。每次随机
在任一盒中取出一根火柴。问其中一盒
中火柴被取完而另一盒中至少还有 5 根
火 柴 的 概 率 有 多 大 ？ （ 用 频 率 估 计 概
率）
function proguji=liti40(mm)
%mm 是随机实验次数
frq=0; randnum=binornd(1,0.5,mm,2*20);proguji=0; for i=1:mm
a1=0;a2=0;j=1;
while (a1<20)&(a2<20)
if randnum(i,j)==1 a1=a1+1; else a2=a2+1; end
j=j+1; end
if abs(a1-a2)>=5 frq=frq+1; end end; proguji=frq/mm
>> liti40(100)
proguji = 0.4800
>> liti40(1000)
proguji = 0.4970
>> liti40(10000)
proguji = 0.4910
>> liti40(100000)
proguji = 0.4984
例 4 两盒火柴，每盒 n 根。每次随机在
任一盒中取出一根火柴。问其中一盒中
火柴被取完而另一盒中至少还有 k 根火
柴的概率有多大？（用频率估计概率）
function proguji=liti4(n,k,mm)
%n 是每盒中的火柴数 %k 是剩余的火柴数
%mm 是随机实验次数
frq=0; randnum=binornd(1,0.5,mm,2*n);proguji=0; for i=1:mm
a1=0;a2=0;j=1;
while (a1<n)&(a2<n)
if randnum(i,j)==1 a1=a1+1; else a2=a2+1; end
j=j+1; end
if abs(a1-a2)>=k , frq=frq+1; end
% a1=a1,a2=a2,frq % pause end; proguji=frq/mm
>> liti4(20,5,100)
proguji = 0.4800
>> liti4(20,5,1000)
proguji = 0.4970
>> liti4(20,5,10000)
proguji = 0.4910
>> liti4(20,5,100000)
proguji = 0.4984
几 何 概 率 模 拟1. 定义
向任一可度量区域 G 内投一点，如果所投的 点落在 G 中任意可度量区域 g 内的可能性与 g 的度量成正比，而与 g 的位置和形状无关
，则称这个随机试验为几何型随机试验。或 简称为几何概型。
2. 概率计 算
P(A)=[A 的度量 ]/[S 的度量 ] 例 5 两人约定于 12 点到 1 点到某地会面，先 到者等 20 分钟后离去，试求两人能会面的概 率？ 解：设 x, y 分别为甲、乙到达时刻 ( 分
钟 )令 A={ 两人能会面 }={(x,y)||x-y|≤20 ， x≤60,y≤6
0}
P(A)=A 的面积 /S 的面积
= （ 602-402 ） /602=5/9=0.5556
function proguji=liti5(mm)
%mm 是随机实验次数
frq=0;
randnum1=unifrnd(0,60,mm,1); randnum2=unifrnd(0,60,mm,1); randnum=randnum1-randnum2; proguji=0;
for ii=1:mm
if abs(randnum(ii,1))<=20 frq=frq+1;
end end
proguji=frq/mm
liti5(10000)
proguji = 0.5557
例 6 在我方某前沿防守地域，敌人以一个炮排
（含两门火炮）为单位对我方进行干扰和破坏． 为躲避我方打击，敌方对其阵地进行了伪装并经 常变换射击地点．
经过长期观察发现，我方指挥所对敌方目标的 指示有 50 ％是准确的，而我方火力单位，在指 示正确时，有 1/3 的概率能毁伤敌人一门火炮， 有 1/6 的概率能全部消灭敌人．
现在希望能用某种方式把我方将要对敌人实施
的 1 次打击结果显现出来，利用频率稳定性，确 定有效射击 ( 毁伤一门炮或全部消灭 ) 的概率 . 复杂概率模拟
分析： 这是一个复杂概率问题，可以通过理论 计算得到相应的概率 .
为了直观地显示我方射击的过程，现采用模拟 的方式。
需要模拟出以下两件事： 1. 问题分析
[1] 观察所对目标的指示正确与否
模拟试验有两种结果，每一种结果出现的概 率都是 1/2 ．
因此，可用投掷一枚硬币的方式予以确定 ,
当硬币出现正面时为指示正确，反之为不正确
．
[2] 当指示正确时，我方火力单位的射击结果 情况
模拟试验有三种结果：毁伤一门火炮的可能性为 1/3( 即 2/6) ，毁伤两门的可能性为 1/6 ，没能毁 伤敌火炮的可能性为 1/2( 即 3/6) ．
这时可用投掷骰子的方法来确定：
如果出现的是１、２、３三个点：则认为没能击 中敌人；
如果出现的是４、５点：则认为毁伤敌人一门火 炮；
若出现的是６点：则认为毁伤敌人两门火炮．
2. 符号假设
i ：要模拟的打击次数；
k1 ：没击中敌人火炮的射击总数； k2 ：击中敌人一门火炮的射击总数； k3 ：击中敌人两门火炮的射击总数．
E ：有效射击 ( 毁伤一门炮或两门炮 ) 的概率
3. 模拟框图
初始化 :i=0,k1=0,k2=0,k3=0 i=i+1
骰子点数 ?
k1=k1+1 k2=k2+1 k3=k3+1 k1=k1+1
i ＜ mm?
E=(k2+k3)/mm 停止
硬币正面 ?
Y N
N
Y
1， 2 ， 3
4， 5 6
function liti6(p,mm) efreq=zeros(1,mm);
randnum1 = binornd(1,p,1,mm);
randnum2 = unidrnd(6,1,mm);k1=0;k2=0;k3=0; for i=1:mm
if randnum1(i)==0 k1=k1+1; else
if randnum2(i)<=3 k1=k1+1;
elseif randnum2(i)==6 k3=k3+1; else k2=k2+1;
end end
efreq(i)=(k2+k3)/i; end
num=1:mm;plot(num,efreq)
在 Matlab 中编辑 .m 文件输入以下命令：
在 Matlab 命令行中输入以下命令：
liti6(0.5,2000)
在 Matlab 命令行中输入以下命令：
liti6(0.5,20000)
5. 理论计算
6. 结果比较
模拟结果与理论计算近似一致，能更加真 实地表达实际战斗动态过程．
三 . 蒙特卡洛模拟的理论基础与模拟结果的误 差
大数定律
中心极限定理
大数定律
贝努里（ Bernoulli ） 大数定律
设 nA 是 n 次独立重复试验中事件 A 发生的 次数 , p 是每次试验中 A 发生的概率，则
0
有
0 lim
n p
P nA
n
或 lim 1
n p
P nA
n
在概率的统计定义中，事件 A 发生的频率
“ 稳定于”事件 A 在一次试验中发生的概率是 指：
n nA
n nA
频率 与 p 有较大偏差
p n
nA
是
小概率事件 , 因而在 n 足够大时 , 可以用频率 近似代替 p . 这种稳定称为依概率稳定 .
贝努里（ Bernoulli ） 大数定律的意义：
设 X1 ， X2 ，…， XN 是来自总体 X (EX<
)
的简单随机样本，即 X1 ， X2 ，…， XN 独立同 分
布，则
EX N X
X N P
i i
N
1
1
1 )
lim
(
X E X
P N
N
辛钦 大数定律
即
中心极限定理
设随机变量序列 X1, X 2,, X n, 相互
独立，服从同一分布，且有期望和方差： , 2 , 1 ,
0 )
( ,
)
(X D X 2 k
E k k
则对于任意实数 x ,
x t
n
k k
n x e dt
n
n X
P 1 2
2
2 lim 1
若令
) 1 , 0 (
1 ~ N
n
n
n X
k k
n
k k
n X
X n
1
1 则
等价于
) 1 , 0 ( / n ~ N Xn
于是
(u ) 1 n
X u P N
这表明，不等式
近似地以概率 1 成立。上式也表明，
n Xn u
Xn 收敛到
的阶为 O(n -1/2) 。通常，蒙特卡罗方法的误差 ε 定义为
n u
关于蒙特卡罗方法的误差需说明两点：
（ 1 ），蒙特卡罗方法的误差为概率误差， 也即蒙特卡罗方法的收敛是概率意义下的 收敛，虽然不能断言其误差不超过某个
值，但能指出其误差以接近 1 的概率不超过 某个界限。
如 =0.5 ，误差 0.6745 / n. 此时，误差 超过 ε 的概率与小于 ε 的概率 1- 相等，都 等于 0.5 。
来代替，在计算所求量的同时，可计算出 。
2 1
1
2 (1 )
ˆ 1
n
i i
n
i i
n X n X
ˆ
关于蒙特卡罗方法的误差需说明两点：
（ 2 ）误差中的均方差 是未知的，必须 使用其估计值
显然，当给定置信度后，误差由 和 n 决定。要减小 ，或者是增大 n ，或者是减 小方差 2 。在 固定的情况下，要把精度提 高一个数量级，试验次数 n 需增加两个数量 级。因此，单纯增大 n 不是一个有效的办法。
四、减小方差的各种技巧
n u
另一方面，如能减小估计的均方差 ， 比如降低一半，那误差就减小一半，这相
当于 n 增大四倍的效果 (n=(u /)2) 。因此 降低方差的各种技巧，引起了人们的普遍 注意。
n u
一般来说，降低方差的技巧，往往会使观 察一个子样的时间增加。在固定时间内，使 观察的样本数减少。所以，一种方法的优 劣，需要由方差和观察一个子样的费用（使 用计算机的时间）两者来衡量。这就是蒙特
卡罗方法中效率的概念。它定义为 nc ，其中 c 是观察一个子样的平均费用。显然
nc= (u /)2 2c 它与 2c 成正比。
四、效率
总而言之，作为提高蒙特卡洛方法效率的 重要方向，是在减小标准差的同时兼顾考虑费 用大小，使 2c 尽可能地小。
nc= (u /)2 2c)
五、蒙特卡洛方法的特点
1) Monte Carlo 方法及其程序结构简单
—— 产生随机数，计算相应的值。即通过大 量的简单重复抽样和简单计算实现该方法。
2 ）收敛速度与问题维数无
—— Monte Carlo关 方法的收敛速度为 O(n -1/2) ， 与一般数值方法相比很慢。因此，用 Monte
Carlo 方法不能解决精确度要求很高的问题
—— 从
n u
可见， Monte Carlo 方法的误差只与标准差 和样本容量 n 有关，而与样本所在空间无关， 即 Monte Carlo 方法的收敛速度与问题维数无 关。而其他数值方法则不然。因此，这就决定 了 Monte Carlo 方法对多维问题的适用性。
—— 在解题时受问题条线限制的影响较小，例 如要计算 s 维空间中的任一区域 Ds 上的积分
s
D g x x xs dx dx dx
g
s
( 1, 2, , ) 1 2
3 ） Monte Carlo 方法的适用性强
时，无论区域 Ds 的形状多么特殊，只要能给出 描述 Ds的几何特征的条件，就可以从 Ds中均匀 产生 n 个点 , 得到积分的近似值(x1(i), x2(i),, xs(i) )
N i
i s i
s i
N g x x x
N g D
1
) ( )
( 2 ) (
1 , , , )
(
其中 Ds 为区域 Ds 的体积。这是数值方法难以作 到的。
3. 某厂生产的灯泡能用 1000 小时的概率为 0.8, 能用 1500 小时的概率为 0.4 , 求已用 1000 小时 的灯泡能用到 1500 小时的概率（频率估计概率）。
2. 在一袋中有 10 个相同的球，分别标有号码 1,2,…, 10 。今有放回任取两个球，求取得的第一个球号码 为奇数，第二个球的号码为偶数的概率（频率估计概率） 1.掷三枚不均匀硬币，每枚正面出现概率为 0.3 ，记 录前 1000 次掷硬币试验中至少两枚都为正面频率的 波动情况，并画图。
作业 :
4 : 两船欲停靠同一个码头 , 设两船到达码头的 时间各不相干，而且到达码头的时间在一昼夜内 是等可能的 . 如果两船到达码头后需在码头停 留的时间分别是 1 小时与 2 小 时，试求在一昼 夜内，任一船到达时，需要等待空出码头的概率 . （频率估计概率）
5 : 在 0,1,2,3,…..,9 中不重复地任取 4 个数，求 它们能排成首位非零的四位偶数的概率 . （频 率估计概率）
6: 从 1,2,….. ,10 十个数字中有放回地任取 5 个 数字 , 求取出的 5 个数字中按由小到大排列 , 中间的那个数等于 4 的概率 . （频率估计概 率）
