4随机模拟

延安大学数学与计算机科学学院
零件参数
[例3] 粒子分离器的某关键参数（记作y）由7个零件的 参数（记作x1，x2，…，x7）决定，
y=174.42( x3 x1 )( )0.85 x 5 x2 x1 x4 0.56 x4 1.16 1 2.62 1 0.36( ) ( ) x2 x2 x6 x7
a=D1 b=D1.*D2./3
x1=normrnd(a(1),b(1),1,1); …… x7=normrnd(a(7),b(7),1,1)
零件期望a 零件标准差b 按指定分布产生随机 数，作为零件参数
y=funY(x1,x2,x3,x4,x5,x6,x7) 计算该产品的关键参 Q=funQ(y） 数y值和损失的大小 产生一组零件参数，相当制作一个产品。重复N＝1000次 （即生产N个产品），求出损失费用的平均.(代码略)
生成N=20天的需求量的matlab代码可以为
延安大学数学与计算机科学学院
报童诀窍的简化版
生成N=20天的需求量的matlab代码（定义函数） function y=randfun1(N) 文件名为randfun1.m y=zeros(1,N); for i=1:N t=rand(1,1); 产生1行N列的全0矩阵，目的 if t<=0.34 分配好矩阵大小，可以省略 X=500; elseif t<=0.70 根据随机数t的范围，确定 X=510; 需求量 X 值，并保存到数 else 组的相应位置中（关键部 X=520; 分） end y(i)=X; end
延安大学数学与计算机科学学院
报童诀窍的简化版
问题2：如何从需求量计算利润？售出价a=1， 购进价b=0.75，退回价c=0.6，购进数n=510 函数文件名fun2.m n(a b) xn 模拟程序代码 function y=fun2(x) y N=1000; (a b) (n x)(b c) x n x a=1; %售出价 x=randfun1(N); b=0.75; %购进价 y=0; c=0.6; %退回价 for i=1:N n=510; %购进数量 y=y+fun2(x(i)); if(x>n) end y=n*(a-b); y/N else y=x*(a-b)-(n-x)*(b-c) end
由于y的表达式很复杂，要想计算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] D2=[5,5,5,15,15,5,5,]./100 %标定值 %容差等级
2
xi～N (ai , bi )
零件损失Q是一个随机变量，平均损失就是期望E(Q) 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)
延安大学数学与计算机科学学院
报童诀窍的简化版
问题1：如何产生以下分布规律的随机数据？ 需求量X 500 510 520 注:rand(m,n)可以生成 概率P 0.34 0.36 0.30 [0,1]上均匀分布随机数
• 把 [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 • 重复多次就可以若干天的需求量了
延安大学数学与计算机科学学院
两船到达码头时刻服从[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(随机数) if((x<=y & y<=x+2) |(y<=x & x<=y+1)) c=c+1; %如果能相遇，则计数器加1 end end P=c/N %显示相遇的概率近似值
概率P
0.34 0.36 0.30
按照需求量的分布规律，随机生成N=20个数据: 510 520 500 510 520 500 500 500 500 510 500 500 500 520 510 520 510 510 代表20天的需求量，计算出报童在这20天的总利 润和平均利润，用平均利润来近似报童的平均收 入。 这也是Monte Carlo方法. 问 •如何按分布规律产生随机数据？ 题 •随机数据很多时，如何编程？
延安大学数学与计算机科学学院
报童诀窍的简化版
也可完整模拟程序：
售出价a 购进价b 退回价c a=1;b=0.75;c=0.6;n=510; 购进数n 总利润y n=510; N=1000; y=0; 模拟天数N for i=1:N t=rand(1,1); if (t<=0.34) x=500; 根据随机数t， elseif (t<=0.70) x=510; 计算需求量x值 else x=520; end 根据需求量x， if(x>n) y=y+n*(a-b); 计算利润并累 else y=y+x*(a-b)-(n-x)*(b-c); 加到y中。 end End fprintf(1,'平均利润=%.3f',y/N); 显示平均利润
延安大学数学与计算机科学学院
二、Matlab中的部分随机数产生命令
注：以下都是产生不同分布m×n 阶随机矩阵
命令 rand(m,n) unifrnd(a,b,m,n) unifrnd(N,m,n) randn(m,n) exprnd(λ,m,n) poissrnd(λ,m,n) normrnd(μ,σ,m,n)
P{ X z } 1 2

z

e xp[ ( t )2 / 2 2 ]dt p
计算函数:z = norminv(p，mu，sigma)
延安大学数学与计算机科学学院
统计直方图matlab函数hist
2000
• 直方图绘图函数: hist(data，n) 1000 0 其中,data是需要处理的数据块, 1 2 3 4 5 绘图原理:利用data中最小数和最大数构成一区间, 将区间等分为n个小区间，统计落入每个小区间的 数据量,以数据量为高度绘小矩形，形成直方图。 如果省略参数n，MATLAB将n的默认值取为10。 • 直方图统计计算函数：N=hist(data，n) 计算结果N是n个数的一维数组，分别表示data中各 个小区间的数据量。这种方式只计算而不绘图。
延安大学数学与计算机科学学院
四、练习
一列火车在下午1点后离开A站前往B站，火车 离开A站的规律如下：
离站时间 13:00 13 : 05 13 :10 概率 0.7 0.2 0.1
火车从A到B所需时间平均为30分钟，标准差 为2分钟。如张三到达B站的时间分布为: 到达时间 13： 28 13： 30 13： 32 13： 34
延安大学数学与计算机科学学院
三、一些随机模拟的例子
[例1](相遇问题) 甲、乙两船在24小时内独立 地随机到达码头. 设两船到达码头时刻都服从 [0,24]上的均匀分布，如果甲船到达码头后停 留2小时,乙船到达码头后停留1小时.问两船相 遇的概率有多大？（可用几何概率，此处略） 分析：如果知道甲、乙两船到达的时刻x和y，两 船能相遇的条件就是：x y x 2或y x y 1 两船到达的时刻x和y可以随机生成，生成一组到 达时刻，可以确定是否能相遇。 如果重复很多次，统计相遇的比例就可近似为相 遇的概率。 这也是Monte Carlo思想.
均值=标定值 标准差=容差等级*标定值÷3
延安大学数学与计算机科学学院
零件参数
x7 0.75 B
标 定 值 容差等级
x1 0. 2 B
x2 0.3 B
x3 x4 x5 x6 0.1 0.1 1.5 Байду номын сангаас6 B C C B
容差等级 A=1% B=5% C=15%
均值=标定值
标准差=容差等级*标定值÷3
3 2
关键参数y的目标值是1.50，当偏离为±0.1 时，产品 为次品，损失为1000元；当偏离为±0.3时，产品为废品， 损失为 9000 元。由于工艺原因， 7 个零件参数可以看着 是正态随机变量，在后面的标定值及容差等级情况下， 求产品的平均损失？ 容差等级 x1 x2 x3 x4 x5 x6 x7 A=1% B=5% 标 定 值 0. 2 0.3 0.1 0.1 1.5 16 0.75 容差等级 B B B C C B B C=15%
延安大学数学与计算机科学学院
计算机随机模拟(Monte Carlo)
（建模培训）
吕佳
延安大学数学与计算机科学学院
一、简介
有些问题，由于随机因素很多，用概率论 的方法进行求解可能很难很复杂，这时就需要 借助随机模拟方法得到近似解答。 随机模拟法也叫蒙特卡罗(Monte Carlo) 方法，也称为计算机随机模拟方法。 由于Monte Carlo法计算量大，精度不高， 因而需要借助计算机，并仅适合一些用解析方 法或常规数值方法难以解决问题的低精度求解 和验证