随机模拟

计算机随机模拟(Monte Carlo)
（建模培训）
吕佳
一、简介
有些问题，由于随机因素很多，用概率论的方法进行求解 可能很难和复杂，这时就需要借助随机模拟方法得到近似解答 。
随机模拟法也叫蒙特卡罗(Monte Carlo)方法，也称为计 算机随机模拟方法。
由于Monte Carlo法计算量大，精度不高，因而需要借助 计算机，并仅适合一些用解析方法或常规数值方法难以解决问 题的低精度求解和验证
if(x>n) y=y+n*(a-b); else y=y+x*(a-b)-(n-x)*(b-c); end
根据随机数t， 计算需求量x值
根据需求量x， 计算利润并累 加到y中。
End
显示平均利润（记作y）由7个零件的 参数（记作x1，x2，…，x7）决定，
均值为λ的指数分布
poissrnd(λ,m,n)
均值为λ的泊松分布
normrnd(μ,σ,m,n) 正态分布N(μ,σ2)
有关正态分布matlab计算函数
• 正态分布变量X的数学期望，方差 2 ，密度函数
f ( x) 1 exp[ ( x )2 / 2 2 ) 2
x (, )
三、一些随机模拟的例子
[例1](相遇问题) 甲、乙两船在24小时内独立地随机到达码头. 设两船到达码头时刻都服从[0,24]上的均匀分布，如果甲船到 达码头后停留2小时,乙船到达码头后停留1小时.问两船相遇的 概率有多大？（可用几何概率，此处略）
分析：如果知道甲、乙两船到达的时刻x和y，两船能相遇的条件就
模拟：通过产生指定分布的随机数，来代表7个零件的 参数值，计算y值，确定损失大小。多做几次可得均值
零件参数
%编写计算y值的函数：y=funY(x1,x2,x3,x4,x5,x6,x7) %编写计算损失费函数:Q=funQ(y)
D1=[0.1,0.3,0.1,0.1,1.5,16,0.75] %标定值
计算该产品的关键参 数y值和损失的大小
产生一组零件参数，相当制作一个产品。重复N＝1000次
（即生产N个产品），求出损失费用的平均.(代码略)
四、练习
一列火车在下午1点后离开A站前往B站，火车离开A站的规律 如下：
离站时间 13:00 13 : 05 13 :10 火车从概A到率B所需时间平0.均7为30分钟0.，2标准差为02.1分钟。如张
E(Q) 0 P(Q 0) 1000 P(Q 1000) 9000 P(Q 9000)
1000 P(0.1 | y 1.5 | 0.3) 9000 P(| y 1.5 | 0.3)
由于y的表达式很复杂，要想计算y的分布和上述概 率很困难，我们必须寻找较为有效的近似方法。
看着是正态随机变量，在后面的标定值及容差等级情况
下，求产品的平均损失？
x1 x2 x3 x4 x5 x6 x7 标 定 值 0.2 0.3 0.1 0.1 1.5 16 0.75 容差等级 B B B C C B B
容差等级 A=1% B=5% C=15%
均值=标定值 标准差=容差等级*标定值÷3
零件参数
• 把 [0,1]分成长度为0.34、0.36、0.30的三个区间 [0,0.34]、 (0.34,0.70]、(0.70,1]
• 用rand(1,1)产生1个[0,1]上均匀分布随机数， 如该数在[0,0.34]、 (0.34,0.70]或(0.70,1]内，相当于该天的需 求量相应为500、510和520
二、Matlab中的部分随机数产生命令
注：以下都是产生不同分布m×n 阶随机矩阵
命令
说明
rand(m,n) unifrnd(a,b,m,n)
[0,1]上均匀分布 [a,b]上均匀分布
unifrnd(N,m,n)
1，…,N的等概率分布
randn(m,n)
标准正态分布N(0,1)
exprnd(λ,m,n)
if((x<=y & y<=x+2) |(y<=x & x<=y+1))
c=c+1;
%如果能相遇，则计数器加1
end
end P=c/N %显示相遇的概率近似值
注意：每次运行的结果一般都不一样
报童诀窍的简化版
[例2]（报童诀窍的简化版）报童每天清晨从报社购进报纸零
售，晚上退回没有卖掉的报纸.若每份报纸的购进价为b=0.75元 ，售出价为a=1元，退回价为c=0.6元.每天需求量X是离散型随 机变量，其分布为
X 500 510 520 P 0.34 0.36 0.30
问：如果报童每天购进报纸为n=510份，每天的平均利润是多少 ？
方法一：概率方法（略） 方法二：如果我们知道每天的需求量，可直接计算利润。 而每天需求量可以按分布生成（随机模拟思想）
售出价a=1，购进价b=0.75，退回价c=0.6，
报童诀窍的简化版
计算密度值函数:normpdf(x，mu，sigma)
• 累积分布函数，即积分上限函数
P{ X x} 1 x exp[ (t )2 / 2 2 ]dt
2
计算概率值函数:normcdf(x，mu，sigma)
• 逆累积分布函数值，即已知概率值p，求z使得
P{ X z} 1 z exp[ (t )2 / 2 2 ]dt p
，用平均利润来近似报童的平均收入。
这也是Monte Carlo方法.
问 •如何按分布规律产生随机数据？ 题 •随机数据很多时，如何编程？
报童诀窍的简化版 问题1：如何产生以下分布规律的随机数据？
需求量X 500 510 520 注:rand(m,n)可以生成 概率P 0.34 0.36 0.30 [0,1]上均匀分布随机数
三到达B站的时间分布为:
问：到张达三时能赶间上该1趟3：火28车的1概3：率30是多1少3：32 13：34 概率 0.3 0.4 0.2 0.1
函数文件名fun2.m function y=fun2(x)
a=1; %售出价 b=0.75; %购进价 c=0.6; %退回价 n=510; %购进数量
if(x>n) y=n*(a-b);
else y=x*(a-b)-(n-x)*(b-c)
end
模拟程n序(a代 b码)
xn
y xNx(==a1r0a0nb0d) ;fu(nn1(xN))(;b c) x n
两船到达码头时刻服从[0,24]上的均匀分布，甲船停留2小时,乙 船停留1小时,相遇概率？
方法二 %计算机模拟程序
clc; N=1000; c=0; %模拟次数N，相遇次数c清零
for i=1:N
%重复N次到达时间
x=unifrnd(0,24,1,1); %甲船到达时间x(随机数)
y=unifrnd(0,24,1,1); %乙船到达时间y(随机数)
• 重复多次就可以若干天的需求量了
生成N=20天的需求量的matlab代码可以为
报童诀窍的简化版
生成N=20天的需求量的matlab代码（定义函数）
function y=randfun1(N)
y=zeros(1,N);
for i=1:N
t=rand(1,1);
if t<=0.34 X=500;
elseif t<=0.70 X=510;
是两：船到达的时刻x和y可以随机生x 成 ，y生成x一组2且到达y时刻x，可y以确1定
是否能相遇。 如果重复很多次，统计相遇的比例就可近似为相遇的概率。
这也是Monte Carlo思想.
两船到达码头时刻服从[0,24]上的均匀分布，甲船停留2小时,乙 船停留1小时,相遇概率？
方法一 可以使用软件（excel、C等实现）
y=174.42( x1 )( x3 )0.85 x5 x2 x1
3
1
2.62
1
0.36(
x4 x2
)0.56
2
(
x4 x2
)1.16
x6 x7
关键参数y的目标值是1.50，当偏离为±0.1时，产品
为次品，损失为1000元；当偏离为±0.3时，产品为废
品，损失为9000元。由于工艺原因，7个零件参数可以
else X=520;
end
y(i)=X;
end
文件名为randfun1.m
产生1行N列的全0矩阵，目的 分配好矩阵大小，可以省略
根据随机数t的范围，确定需求量X 值，并保存到数组的相应位置中（ 关键部分）
报童诀窍的简化版 问题2：如何从需求量计算利润？售出价a=1， 购进价b=0.75，退回价c=0.6，购进数n=510
y=0;
for i=1:N
y=y+fun2(x(i));
end
y/N
报童诀窍的简化版
也可完整模拟程序：
售出价a 购进价b 退回价c
a=1;b=0.75;c=0.6;n=510; n=510; N=1000; y=0;
购进数n
总利润y
for i=1:N
模拟天数N
t=rand(1,1); if (t<=0.34) x=500; elseif (t<=0.70) x=510; else x=520; end
D2=[5,5,5,15,15,5,5,]./100
%容差等级
a=D1 b=D1.*D2./3
零件期望a 零件标准差b
x1=normrnd(a(1),b(1),1,1); 按指定分布产生随机
……
数，作为零件参数
x7=normrnd(a(7),b(7),1,1)
y=funY(x1,x2,x3,x4,x5,x6,x7) Q=funQ(y）