随机模拟实验

计算结果N是n个数的一维数组，分别表示data中各个小区间的数据量。 这种方式只计算而不绘图。
数学实验
例： 统计10000个均匀随机数在五个小区间的分布。
即观察10000 个随机数在[0,0.2],[0.2,0.4], [0.4,0.6],[0.6,0.8],[0.8,1]的分布情况。 data=rand(10000,1); figure(1),hist(data,5) N5=hist(data,5) figure(2),bar(N5,'r')
(1)甲比乙先到码头：x
数学实验
function F=shipmeet(N) if nargin==0, N=2000; end P=24*rand(2,N); X=P(1,:);Y=P(2,:); I=find(X<=Y&Y<=X+2); J=find(Y<=X&X<=Y+1); F=(length(I)+length(J))/N plot(X,Y,’b.’) ,hold on axis([0 24 0 24]) h1=plot([0 24],[0,24],’k’) h2=plot([1 24],[0,23],’r’) h3=plot([0 22],[2,24],’b’) set([h1,h2,h3],’linewidth’,2)
0, 0 r 0.5 f(r ) 0.5 r 1 1,
如果 0.5<r<1，认为出现数字面X=1，否则出现图案面
数学实验
算法设计： 输入：抛硬币次数n
输出：抛硬币出现数字面的概率估计值 p
1. 2. 3. 4. 5. 6.
初始化count=0 对于i=1,2,…,n，执行第3~4步 产生[0,1]上的随机数r 如果0.5<r<1,则count=count+1 计算p=count/n 输出p
人之一数学家冯· 诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种 方法，为它蒙上了一层神秘色彩。 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人 们就知道用事件发生的“频率”来决定事件的“概率”。 1777年，法国数学家
布丰（ Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法 求圆周率π，这被认为是蒙特卡罗方法的起源。上世纪 40年代电子计算机的出现， 特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速 地模拟这样的试验成为可能。
数学实验
统计直方图的绘制：hist(data，n)
其中data是需要处理的数据块。 绘图原理：利用data中最小数和最大数构成一区间 ,将区间等分为n个 小区间，统计落入每个小区间的数据量。以数据量为高度绘小矩形， 形成直方图。 如果省略参数n，MATLAB将n的默认值取为10。
直方图也可以用于统计计 N=hist(data，n)
>> A=[0 0 2 3 4 5 6 7 8];
>> [a,b]=find(A,1,'last')
a=
1
b=
9
数学实验
条形图的绘制：bar(data)或 bar3(data)
x=linspace(0,pi,10); y=sin(x); bar(y,'r') bar3(y,'r')
1
0.5
0
1 2 3 4 5 6 7 8 9 10
随机模拟实验
电子科技大学 数学科学学院
数学实验
随机模拟实验
背景 引例：模拟抛硬币实验
Matlab随机函数
蒙特卡罗方法
数学实验
（一）背景
随机模拟实验----应用广泛
解定积分、二重积分、三重积分 解优化模型
物流系统、雷达系统 模拟工厂的运行、股票交易、城市交通、军事演习
O
1
问题： [a，b]上的均匀随机数如何产生？ 一般区间 [a，b]上的均匀随机数 产生方法
R=a+(b-a)*rand
观察均匀性：
例 . 观察 1000 个随机数在[0,0.5],[0.5,1]分布情况
function F=myrand(n) if nargin==0,n=1000;end
X=rand(1,n); Index=find(X<0.5);
（1）返回向量中非零元素的位置 find(A)
>> A = [1 0 4 -3 0 0 0 8 6]; >> X = find(A)
X=
13489
数学实验
（2）返回矩阵中非零元素的位置 find(A) >> A= [3 2 0; -5 0 7; 0 0 1]; >> find(A) ans = 1 2 4 8 9
数学实验
版本2 R=rand(1, 1000); P=sum(R>0.5)/1000
数学实验
（三）Matlab随机函数
数学实验
均匀分布随机数
[0,1]内均匀随机数产生方法 : rand( ) rand(m) rand(m,n) rand([m,n]) 产生 m× m的均匀随机数 产生 m× n的均匀随机数 等效于rand(m,n)
目的：验证、发现、求解、探究
数学实验
（二）引例：模拟抛硬币实验
实验目的：重复多次抛一枚硬币，观察出现数字面的次数，如何模拟？
X 表示每次出现的数字面的次数，则
X=0 代表图案面； X=1 代表数字面；
理论上： P{X=0}=P{X=1}=1/2
分析：定义函数f(r), 这里 r为[0,1]内的随机数
x2 y x
模拟过程 : (1)在矩形区域内投 N个点 (2)统计落入所围区域内点的个数M (3)所围区域的面积= M/N
程序： N=1000; data=rand(N,2); x=data(:,1); y=data(:,2); II=find(y<=sqrt(x)&y>=x.ˆ2); M=length(II); S=M/N 运行结果 S= 0.3276
Monte Carlo名字的由来：
• Monte Carlo是摩纳哥（monaco)的首都， 该城以赌博闻名
Monte-Carlo, Monaco
Monte Carlo方法的提出：
蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方 法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持
运行结果：
F= 0.1175
20
15
10
5
0
0
5
10
15
20
概率值：P = 0.1207
平面多边形填充图方法 fill( )
x1=0:.01:1;y1=sqrt(x1); x2=1:-.01:0;y2=x2.^2; fill([x1,x2],[y1,y2],'r')
y =x2 , x = y 2 所围区域
数学实验
例3：相遇问题
甲 、乙两船在24小时内独立地随机到达码头 。设两船到达码头时刻分 别为 X 和 Y，均匀分布随机变量 X ~ U(0 , 24), Y ~ U(0 , 24)。如果甲船到
达码头后停留 2小时，乙船到达码头后停留
1小时 。问两船相遇的概率有多大？
相遇条件：
y且 y x 2 (2)乙比甲先到码头：y x且 x y 1
要求：模拟抛1000次硬币的过程，统计出现两面的频数。
问题1 ：写程序前需要设计算法吗？你如何做的？
问题2 ：这个算法有没有发挥Matlab语言特色？如果没有，如何修改？
数学实验
编写程序：
版本1 c=0; n=1000; for i=1:n r=rand; % 模拟一次投币 if r>0.5 %{0.5<r<1}={X=1}={出现数字面} c=c+1 end end p=c/n % 出现数字面的频率 注：版本1并未能发挥matlab的特点
为 16。积分值是一个三维图形所围体积 (V1)，该三维图形位于立方体区域
{(x，y，z) |0≤ x ≤4，–1≤ y ≤2，0 ≤ z ≤16 }
该立方体区域的体积(V2)为 192. 2.蒙特卡罗方法求解： 向V2内投N个随机点，统计落在V1内的个数M
程序:
function testmain N=100000; for k=1:7 % 多次模拟 V1(k)=mysim(N); end V1 function v =mysim(N) V2=192; %V2=4*3*16 d =rand(N,3); x =4*d(:,1); y =-1+3*d(:,2); z =16*d(:,3); % 下面表达式可结合find和 l.^2) & (x<=y+2) & (z<=x.*y.^2) ); v =V2*M/N;
V1 =
7.5859 7.5898 7.6262
7.5629 7.7894 7.6781
7.4304
符号结果 ：7.5857
例2：曲线所围区域面积
用随机模拟方法估算两条抛物线 y =x2 , x = y 2所围图形的面积 分析 : (1)两条曲线交点， [0,0],[1,1] (2)用一个矩形区域包含所围区域 0 y 1 0 x 1 (3)所围区域内的点满足
y1=-1:.1:2;y2=2:-.1:-1; x11=y1.*y1;x22=y2+2; fill([x11,x22],[y1,y2],'r')
y= x – 2 与 y2 = x 所围区域
x1=-1:0.1:1; y1=x1.^2.^(1/3); x2=1:-0.1:-1; y2=2-x2.^2; fill([x1,x2],[y1,y2],'c')