(Monte_Carlo)蒙特卡罗模拟

Monte Carlo名字的由来：
• Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名
Nicholas Metropolis (1915-1999)
Monte-Carlo, Monaco
Monte Carlo方法的基本思想
蒙特卡罗方法，或称计算机随机模拟方法，是一种基 于‚随机数‛的计算方法。源于美国在第二次世界大战研 制原子弹的‚曼哈顿计划‛，该计划的主持人之一数学家 冯· 诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来 命名这种方法，为它蒙上了一层神秘色彩。 蒙特卡罗方法的基本思想很早以前就被人们所发现和 利用。早在17世纪，人们就知道用事件发生的‚频率‛来 决定事件的‚概率‛。19世纪人们用蒲丰投针的方法来计 算圆周率π，上世纪40年代电子计算机的出现，特别是近 年来高速电子计算机的出现，使得用数学方法在计算机上 大量、快速地模拟这样的试验成为可能。
问题分析
需要模拟出以下两件事：
[1] 观察所对目标的指示正确与否 模拟试验有两种结果，每一种结果出现的概率都是1/2． 因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为 指示正确，反之为不正确． [2] 当指示正确时，我方火力单位的射击结果情况 模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)， 毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)． 这时可用投掷骰子的方法来确定： 如果出现的是１、２、３三个点：则认为没能击中敌人； 如果出现的是４、５点：则认为毁伤敌人一门火炮； 若出现的是６点：则认为毁伤敌人两门火炮．
非常见分布的随机数的产生
• 逆变换方法
由 定 理 1 ， 要 产 生 来 自 F (x) 的 随 机 数 ， 只 要 先 产 生 来 自 U ( 0 ,1) 随 机 数 u ， 然 后 计 算 F 可。具体步骤如下：
1
(u ) 即
(1 ) 生 成 (0 ,1) 上 均 匀 分 布 的 随 机 数 U 。
蒲丰投针实验： 法国科学家蒲丰(Buffon)在1777年提出的蒲丰 投针实验是早期几何概率一个非常著名的例子。蒲丰 投针实验的重要性并非是为了求得比其它方法更精确 的π值，而是它开创了使用随机数处理确定性数学问 题的先河，是用偶然性方法去解决确定性计算的前导。 由此可以领略到从‚概率土壤‛上开出的一朵瑰丽的鲜 花-蒙特卡罗方法(MC) 蒲丰投针实验可归结为下面的数学问题：平面上画 有距离为a的一些平行线，向平面上任意投一根长为l (l<a)的针，假设针落在任意位置的可能性相同，试求针 与平行线相交的概率P(从而求π)
蒲丰投针实验：
如右图所示，以M表示针落下 后的中点，以x表示M到最近一条平行
线的距离，以φ表示针与此线的交角：
针落地的所有可能结果满足： 0 x a / 2 , 0 其样本空间视作矩形区域Ω, 面积是: S a / 2 针与平行线相交的条件： x l sin / 2 , 0 0 它是样本空间Ω子集A，面积是：S A l sin / 2 d l 0 syms l phi; int('l/2*sin(phi)',phi,0,pi) %ans=l 因此，针与平行线相交的概率为： p S ( A ) / S ( ) 2 l / a
举例
例1 在我方某前沿防守地域，敌人以一个炮排（含两门火炮） 为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地 进行了伪装并经常变换射击地点． 经过长期观察发现，我方指挥所对敌方目标的指示有50％是 准确的，而我方火力单位，在指示正确时，有1/3的射击效果能 毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮． 现在希望能用某种方式把我方将要对敌人实施的20次打击结 果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。 分析：这是一个概率问题，可以通过理论计算得到相应的概 率和期望值.但这样只能给出作战行动的最终静态结果，而显 示不出作战行动的动态过程. 为了能显示我方20次射击的过程，现采用模拟的方式。
•马氏链蒙特卡罗方法 为要模拟服从给定分布的随机变量,用生成 一个易于实现的不可约遍历链 X { X , n 0 } 作为随机样本，使其平稳分布为给定分布 的方法，称为马氏链蒙特卡罗方法.
n
[1] 连续型随机变量(以指数分布为例)：
e f (t ) 0
t
1 -e f ( t ) d t F ( x ) 0 t 0 t 0
4，5
k1=k1+1
k2=k2+1
k3=k3+1 Y
k1=k1+1
i＜20? N
E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20 停止
模拟结果
试验 序号 １ ２ ３ ４ ５ ６ ７ ８ ９ １０ 投硬币 结 果 正 正 反 正 正 反 正 正 反 反 ∨ ∨ ∨ ∨ ∨ ∨ ∨ ３ ６ ∨ ∨ 指示 正确 ∨ ∨ ∨ １ ２ 指 示 不正确 掷骰子 结 果 ４ ４ ∨ ∨ ∨ ∨ ∨ ∨ ０ 消灭敌人火炮数 １ ∨ ∨ ２
一些人进行了实验，其结果列于下表 ：
实验者 沃尔弗(Wolf) 年份 1850 投计次数 5000 π的实验值 3.1596
斯密思(Smith)
福克斯(Fox) 拉查里尼 (Lazzarini)
1855
1894 1901
3204
1120 3408
3.1553
3.1419 3.1415929
蒙特卡罗投点法是蒲丰投针实验的推广： 在一个边长为a的正方形内随机投点， 该点落在此正方形的内切圆中的概率 y (a/2,a/2) 应为该内切圆与正方形的面积比值， 即 π a/2 2 : a 2 π/4 n=10000; a=2; m=0; o x for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)^2+(y-a/2)^2 <= (a/2)^2 ) m=m+1; end end disp(['投点法近似计算的π为: ',num2str(4*m/n)]);
rand('seed',0.1); rand(1) ％每次运行程序产生的值是相同的 rand('state',sum(100*clock)*rand);
rand(1) %每次重新启动matlab时，输出的随机数不一样
注意: 产生一个参数为λ的指数分布的随机数应输入 exprnd(1/λ)
产生m×n阶参数为A1,A2,A3的指定分布'name'的随机数矩阵 random('name',A1,A2,A3,m,n) 举例: 产生2×4阶的均值为0方差为1的正态分布的随机数矩阵 random('Normal',0,1,2,4) 'name'的取值可以是(详情参见help random)： 'norm' or 'Normal' / 'unif' or 'Uniform' 'poiss' or 'Poisson' / 'beta' or 'Beta' 'exp' or 'Exponential' / 'gam' or 'Gamma' 'geo' or 'Geometric' / 'unid' or 'Discrete Uniform' ……
基本思想
由上面的例子可以看出，当所求问题的解 是某个事件的概率，或者是某个随机变量的数 学期望，或者是与之有关的量时，通过某种试 验的方法，得出该事件发生的频率，再通过它 得到问题的解。这就是蒙特卡罗方法的基本思 想。 蒙特卡罗方法的关键步骤在于随机数的产 生，计算机产生的随机数都不是真正的随机数 (由算法确定的缘故)，如果伪随机数能够通过 一系列统计检验，我们也可以将其当作真正的 随机数使用。
符号说明
i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数； k3：击中敌人两门火炮的射击总数； E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数．
模拟框图
初始化:i=0,k1=0,k2=0,k3=0 i=i+1 Y
1，2，3
硬币正面?
6
N
骰子点数?
(2 ) 计 算 X F
-1
(U ) , 则 X 为 来 自 F ( x )
分布的随机数.
Acceptance-Rejection 方法
最早由 Von Neumann提出，现在已经广 泛应用于各种随机数的生成。 基本思路： 通过一个容易生成的概率分布 g 和一个取舍 准则生成另一个与 g 相近的概率分布 f 。
E=7/20=0.35
E1 0
13 20
1
4 20
2
3 20
=0.5
示不正 观察所对目标指示正确
确
A 0： 射 中 敌 方 火 炮 的 事 件 ； A 1 ： 射 中 敌 方 一 门 火 炮 的 事 件 ； A 2： 射 中 敌 方 两 门 火 炮 的 事 件 ． 则由全概率公式： E = P ( A 0) = P ( j = 0 ) P ( A 0∣ j = 0 ) + P ( j = 1 ) P ( A 0∣ j = 1 )
试验 序号 １１ １２ １３ １４ １５ １６ １７ １８ １９ ２０
投硬币 结 果 正 反 正 反 正 正 正 正 反 正
指示 正确 ∨
指 示 不正确
掷骰子 结 果 ２ ０ ∨ ∨ ３ ∨ ∨ ６ ４ ２ ４ ∨
消灭敌人火炮数 １ ２
∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ６
∨ ∨
∨ ∨ ∨
从以上模拟结果可计算出：
[2] 离散型随机变量(以离散分布为例)：
x=[2,4,6,8]; px=[0.1,0.4,0.3,0.2]; %以下为程序片段 Fx=0; for n=1:length(px),Fx=[Fx,sum(px(1:n))];end r=rand; index=find(r<Fx); x(index(1)-1) %通用离散分布随机数产生程序
x
x
x 0 x 0
syms t x lambda; Fx=int('lambda*exp(-lambda*t)',t,0,x) %分布函数 syms r; Fxinv=finverse(Fx,x); %求反函数 Fxinv=subs(Fxinv,x,r) %替换反函数变量x为r Fxinv=inline(Fxinv) x=Fxinv(3,rand) %产生参数 lambda=3 指数分布的随机数 %指数分布随机数产生函数已经提供 exprnd(1/3,1,1)
Monte Carlo 模拟
内容提纲
1.引言 2.Monte Carlo模拟基本思想 3.随机数生成函数 4.应用实例举例 5.排队论模拟 6.Monte Carlo模拟求解规划问题
引言(Introduction)
Monte Carlo方法：
蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四 十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法， statistical simulation method 利用随机数进行数值模拟的方法
从而有： 2 l / ap
蒲丰投针实验的计算机模拟： format long; %设置15位显示精度 a=1; l=0.6; %两平行线间的宽度和针长 figure; axis([0,pi,0,a/2]); %初始化绘图板 set(gca,'nextplot','add'); %初始化绘图方式为叠加 counter=0; n=2010; %初始化计数器和设定投针次数 x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n); %样本空间Ω for i=1:n if x(i)<l*sin(phi(i))/2 %满足此条件表示针与线的相交 plot(phi(i),x(i),‘r.’); counter=counter+1; %统计针与线相交的次数 frame(counter)=getframe; %描点并取帧 end end fren=counter/n; pihat=2*l/(a*fren) %用频率近似计算π figure(2) movie(frame,1) %播放帧动画1次