随机模拟实验

要求：模拟抛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的特点
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')
（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
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);
运行结果：
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 所围区域
数学实验
统计直方图的绘制：hist(data，n)
其中data是需要处理的数据块。 绘图原理：利用data中最小数和最大数构成一区间 ,将区间等分为n个 小区间，统计落入每个小区间的数据量。以数据量为高度绘小矩形， 形成直方图。 如果省略参数n，MATLAB将n的默认值取为10。
直方图也可以用于统计计 N=hist(data，n)
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
数学实验
例3：相遇问题
甲 、乙两船在24小时内独立地随机到达码头 。设两船到达码头时刻分 别为 X 和 Y，均匀分布随机变量 X ~ U(0 , 24), Y ~ U(0 , 24)。如果甲船到
达码头后停留 2小时，乙船到达码头后停留
1小时 。问两船相遇的概率有多大？
相遇条件：
y且 y x 2 (2)乙比甲先到码头：y x且 x y 1
6
7
8
>> find(A,3) ans = 3 4 5
数学实验
（5）返回最后一个非零值的位置 find(A,1,,'last') >> A=[0 0 2 3 4 5 6 7 8] A=
0 0 2 3 4 5 6 7 8
>> find(A,1,'last')
ans = 9
数学实验
（6）返回最后一个非零值的行列位置
N5 = 1969 2010 2018 1999 2004
2500 2000 1500 1000 500 0 0 0.2 0.4 0.6 0.8 1
2500 2000 1500
1000 500 0 1 2 3 4 5
（四） 蒙特卡罗方法
蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪 四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方 法，statistical simulation method 利用随机数进行数值模拟的方法
为 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和 length完成 M =sum((x>=y.^2) & (x<=y+2) & (z<=x.*y.^2) ); v =V2*M/N;
数学实验
A=
3 -5 0
2 0 0
0 7 1
（3）返回满足某条件的值find(A**) >> A= [3 2 0; -5 0 7; 0 0 1];
>> find(A>2)
ans =
1 8
数学实验
（4）返回前N个非零元素的位置find(A,X)
X必须为正整数，我们以A=[0 0 2 3 4 5 6 7 8]为例，在MTALAB主窗口中输 入如上向量，然后输入find(A,1)，回车，返回第一个非零值的位置，然后继 续输入 find(A,3)，返回前三个非零值的位置。 >> A=[0 0 2 3 4 5 6 7 8] A= 0 0 2 3 4 5 >> find(A,1) ans = 3
Monte Carlo名字的由来：
• Monte Carlo是摩纳哥（monaco)的首都， 该城以赌博闻名
Monte-Carlo, Monaco
Monte Carlo方法的提出：
蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方 法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持
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
数学实验
版本2 R=rand(1, 1000); P=sum(R>0.5)/1000
数学实验
（三）Matlab随机函数
数学Hale Waihona Puke Baidu验
均匀分布随机数
[0,1]内均匀随机数产生方法 : rand( ) rand(m) rand(m,n) rand([m,n]) 产生 m× m的均匀随机数 产生 m× n的均匀随机数 等效于rand(m,n)
运计算、量子热力学计算、空气动力学计算、核工程)等领域应用广泛。
蒙特卡罗方法在数学中最常见的应用就是蒙特卡罗积分。
数学实验
例1：估算二重积分
计算二重积分
D
xy2 dxdy
其中 D为 y= x – 2 与 y2 = x 所围区域 。
分析 ： 1.二重积分的几何意义是计算体积。
由于 D的边界曲线交点为 ：(–1,1),(4,2) 被积函数在求积区域内的最大值
f1=length(Index);F=[f1,n-f1]; 第一次实验: 490 510 第三次实验: 508 492
观察结果，你能得到什么结论？
数学实验
第二次实验: 497 503 第四次实验: 511 489
回顾：find函数的用法
find函数用于返回所需要元素的所在位置。
位置的判定：在矩阵中，第一列开始，自上而下，依次为1，2 ，3...,然后再从第二列，第三列依次往后数
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)所围区域内的点满足
>> 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
人之一数学家冯· 诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种 方法，为它蒙上了一层神秘色彩。 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人 们就知道用事件发生的“频率”来决定事件的“概率”。 1777年，法国数学家
布丰（ Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法 求圆周率π，这被认为是蒙特卡罗方法的起源。上世纪 40年代电子计算机的出现， 特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速 地模拟这样的试验成为可能。
目的：验证、发现、求解、探究
数学实验
（二）引例：模拟抛硬币实验
实验目的：重复多次抛一枚硬币，观察出现数字面的次数，如何模拟？
X 表示每次出现的数字面的次数，则
X=0 代表图案面； X=1 代表数字面；
理论上： P{X=0}=P{X=1}=1/2
分析：定义函数f(r), 这里 r为[0,1]内的随机数
Monte Carlo方法的基本思想：
当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，
通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率， 或者得到这个随机变量的某些数字特征，并将其作为问题的解。
Monte Carlo方法的应用：
蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学(如粒子输
随机模拟实验
电子科技大学 数学科学学院
数学实验
随机模拟实验
背景 引例：模拟抛硬币实验
Matlab随机函数
蒙特卡罗方法
数学实验
（一）背景
随机模拟实验----应用广泛
解定积分、二重积分、三重积分 解优化模型
物流系统、雷达系统 模拟工厂的运行、股票交易、城市交通、军事演习
(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)
计算结果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')