matlab

概率统计应用实验随机数与统计直方图 蒙特卡罗方法

matlab

均匀分布随机数

O

1

MATLAB产生均匀随机数方法 rand(m,n) 产生均匀随机数方法: 产生均匀随机数方法 ， 产生m× 个 之间均匀随机数 均匀随机数。 产生 ×n个 0,1 之间均匀随机数。随机数等可能落入 区间[0, 内长度相等子区间中 内长度相等子区间中。 区间 ，1]内长度相等子区间中。 观察1000个随机点分布情况 例1. 观察 个随机点分布情况 P=rand(2,1000); x=P(1,:);y=P(2,:); plot(x,y,'b.')

2/18

matlab

统计直方图2000

直方图绘图命令: 直方图绘图命令 hist(data,n) ,

1000 0

1 2 3 4 5 其中,data是需要处理的数据块 是需要处理的数据块, 其中 是需要处理的数据块 绘图原理:利用 利用data中最小数和最大数构成一区间 将 中最小数和最大数构成一区间,将 绘图原理 利用 中最小数和最大数构成一区间 区间等分为n个小区间 个小区间， 区间等分为 个小区间，统计落入每个小区间的数据 以数据量为高度绘小矩形，形成直方图。 量。以数据量为高度绘小矩形，形成直方图。如果省 略参数n, 的默认值取为10。 略参数 ，MATLAB将n的默认值取为 。 将 的默认值取为 直方图也可以用于统计计算

N=hist(data,n) , 计算结果N是 个数的一维数组 分别表示data中各个 个数的一维数组， 计算结果 是n个数的一维数组，分别表示 中各个 小区间的数据量。这种方式只计算而不绘图。 小区间的数据量。这种方式只计算而不绘图。3/18

matlab

统计10000个均匀随机数在五个 例2 统计 个均匀随机数在五个 小区间的分布 。 data=rand(10000,1); hist(data,5) N5=hist(data,5) N5 = 1969 2010 2018 1999 1 条形图是根据数据绘小矩形或 小柱体。使用格式: 小柱体。使用格式 bar(data) 0.5 或bar3(data) 0 x=linspace(0,pi,10); y=sin(x); bar(y,'r') bar3(y,'r')

2500 2000 1500 1000 500 0 0 0.2 0.4 0.6 0.8 1

2004

1 2 3 4 5 6 7 8 9 10

4/18

matlab

3000

产生正态分布随机数的函数为 randn(),使用格式为 ， R=randn(m,n) ,

2000 1000 0 -4

-2

0

2

4

产生m× 阶矩阵 矩阵中元素都是区间( , ) 阶矩阵R, 产生 ×n阶矩阵 ,矩阵中元素都是区间(– 3,3)内的 正态随机数。 正态随机数。 创建10000个正态随机数，将区间 ，3]分为十三 个正态随机数， 例3 创建 个正态随机数 将区间[–3, 分为十三 个小区间，分别绘频数和频率直方图。 个小区间，分别绘频数和频率直方图。 data=randn(10000,1); 0.3 N=hist(data,13); 0.2 figure(1),bar([-3:0.5:3],N,'r') 0.1 figure(2),M=N/10000; 0 -3 -2 -1 0 1 2 3 bar([-3:0.5:3],M,'r')5/18

matlab

均值

mean( )

方差 var( )

[a,b,c,d]=normfit(X,Pi) 计算随机变量X的均值和方差， 计算随机变量X的均值和方差，统计均值和方差 置信度为P 置信度为Pi的置信区间 例4 data=randn(

1,100); [a,b,c,d]=normfit(data,0.95)

6/18

matlab

蒙特卡罗方法，或称计算机随机模拟方法， 蒙特卡罗方法，或称计算机随机模拟方法，是一种基 随机统计”的计算方法。 于“随机统计”的计算方法。方法源于美国在第二次 世界大战中研制原子弹的“曼哈顿计划” 世界大战中研制原子弹的“曼哈顿计划”。 所围图形的面积. 例5 计算两条抛物线 y =x2 ,x = y 2 所围图形的面积 在正方形区域D内投入 个点 在正方形区域 内投入N个点，统计坐标满足 内投入 个点，x2 ≤ y ≤ x

的点P(x, 的数目 的数目M。 的点 ，y)的数目 。面积近似 计算公式为： 计算公式为：S=M/N 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.32767/18

matlab

相遇问题: 乙两船在24小时内独立地随机到 例6 相遇问题 甲、乙两船在 小时内独立地随机到 达码头. 达码头 设两船到达码头时刻分别为 X 和 Y 均匀分布随机变量 X ~ U(0 , 24), Y ~ U(0 , 24) 如果甲船到达码头后停留2小时 小时,乙船到达码头后停留 如果甲船到达码头后停留 小时 乙船到达码头后停留 1小时.问两船相遇的概率有多大？ Y 小时. 小时 问两船相遇的概率有多大？

24 S1 S 2 P {( X , Y ) ∈ D} = 24 22

24

S1 S2O 24

S1 = 0.5 × 22 2 S 2 = 0.5 × 23 2

X

D = {( x , y ) | x 1 ≤ y ≤ x + 2,0 < x , y < 24}

8/18

matlab

例7 相遇问题的统计试验 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 F= 0.1185

D1 = {( x , y ) | x ≤ y & y ≤ x + 2} D2 = {( x , y ) | y ≤ x & x ≤ y + 1}

S1 = 0.5 × 22 2

S 2 = 0.5 × 23 2

24 2 S1 S 2 P{( X , Y ) ∈ D} = = 0.1207 2 249/18

matlab

Y 服从二项分布

Y ~ B( n, p ) k =0,1,2,…,n

k P{Y = k} = Cn p k (1 p ) n k

10/18

matlab

计算二项分布随机变量X=k的命令使用格式为 的命令使用格式为 计算二项分布随机变量 Pk=binopdf(k,n,p) , , 其中， 是随机变量取值 是随机变量取值， 是贝努里试验的重数 是贝努里试验的重数， 其中，k是随机变量取值，n是贝努里试验的重数，p 重贝努里试验中事件A发生的概率 为n重贝努里试验中事件 发生的概率。 重贝努里试验中事件 发生的概率。 对于二项分布随机变量X,计算累加概率P{X ≤ k}的 对于二项分布随机变量 ，计算累加概率 的 MATLAB命令使用格式为 MATLAB命令使用格式为 P = binocdf(k,n,p) , , MATLAB的二项分布随机数发生器使用格式为 的二项分布随机数发生器使用格式为 R= binornd(n,p,L,M) , , , 二项分布随机数。 产生 L×M 个二项分布随机数。 ×

11/18

matlab

例8 有一千名以上的小学生参加保险公司的平安保险, 有一千名以上的小学生参加保险

公司的平安保险, 参加保险的小学生每人一年交保险费50元，若一年内 参加保险的小学生每人一年交保险费 元 出现意外事故，保险公司赔付一万元。统计表明， 出现意外事故，保险公司赔付一万元。统计表明，每 年一千名小学生中平均有两名学生出事故。 年一千名小学生中平均有两名学生出事故。保险公司 赔本的概率有多大？利用二项分布随机数进行模拟， 赔本的概率有多大？利用二项分布随机数进行模拟， 分析：小学生出意外事故的概率为p=0.002,设随机变 分析：小学生出意外事故的概率为 ， 为一年内出事故的小学生人数。 服从二项分布 量X为一年内出事故的小学生人数。X服从二项分布 为一年内出事故的小学生人数 B(n,p),其中 为投保人数。由于对出事故的小学生， 为投保人数。 , ,其中n为投保人数 由于对出事故的小学生， 保险公司一次性赔付一万元， 保险公司一次性赔付一万元，所以每年保险公司赔付 费为： (万元)。 )。一年中保险公司赔付费不超过总 费为：X(万元)。一年中保险公司赔付费不超过总 的保险收费则会获利， 的保险收费则会获利，如果赔付费超过总的保险收费 将会赔本。 将会赔本。每年保险公司所获利润为总保险收费减去 总的赔付费。 总的赔付费。12/18

matlab

N=input(‘input N:=’); p=0.002; join=50;pay=10000; all=join*N X1=fix(all/pay) P1=1-binocdf(X1,N,p) puples=binornd(N,p,1,8) Pays=pay*puples profits=all-Pays

%赔付最大承受人数 %赔偿概率 %八年出事故人数模拟 %八年赔付金模拟 %八年利润模拟

13/18

matlab

正态分布变量X的数学期望 正态分布变量 的数学期望µ,方差σ 2 ,密度函数 1 x ∈ ( ∞ , + ∞ ) f ( x) = exp[ ( x µ ) 2 / 2σ 2 ) 2π σ 计算命令： 计算命令：y = normpdf(x,mu,sigma) , , 累积分布函数， 累积分布函数，即积分上限函数P{ X < x } = 1 2π σ

∫

x

∞

exp[ ( t µ ) 2 / 2σ 2 ]dt

计算命令 ：p = normcdf(x,mu,sigma) , , 逆累积分布函数值，即已知概率值p, 逆累积分布函数值，即已知概率值 ，求z 使得P{ X < z } = 1 2π σ

∫

z

∞

exp[ ( t µ ) 2 / 2σ 2 ]dt = p

计算命令 ：z = norminv(p,mu,sigma) , ,14/18

matlab

某城市中99%男子身高介于 男子身高介于1.52米到 米到1.88米，如果 例9 某城市中 男子身高介于 米到 米 男子上公交车时头与车门相碰的概率小于5%,公交车 男子上公交车时头与车门相碰的概率小于 ， 门的高度应该是多少？ 门的高度应该是多少？ 分析：设身高为正态分布随机变量X, 分析：设身高为正态分布随机变量 ，1.70为X的数学 为 的数学 期望，方差取为36。计算z,使得P{ 期望，方差取为 。计算 ，使得 X > z} = 0.05。

即 。 求逆累积函数在 x=0.95 处的值0.1 mu=170;sigam=6; z=norminv(0.95,mu,sigam) 0.05 data=mu+sigam*randn(1000,1); 0 II=find(data>=z); 150 F=length(II)/1000

160

170

180

190

z = 179.8691 F= 0.004715/18

matlab

练习与思考题1.美国总统选举前民意测验的抽样调查与计算面积的 美国总统选举前民意测验的抽样调查与计算面积的 美国总统选举前民意测验的抽样调查 蒙特卡罗方法有何相同之处? 蒙特卡罗方法有何相同之处 2.甲、乙两人在下午 点到 点之间独立地随机到达 乙两人在下午1点到 点到2点之间独立地随机到达 . 汽车站,这段时间内有四趟班车 开车时间分别为1:15, 这段时间内有四趟班车,开车时间分别为 汽车站 这段时间内有四趟班车 开车时间分别为 1:30, 1:45, 2:00;问在 (1)见车就乘 (2)最多等一趟 见车就乘, 最多等一趟 ；问在: 见车就乘 两种情况下, 车;两种情况下 两人同乘一辆车的概率多大 两种情况下 两人同乘一辆车的概率多大?

3. 贝努里概型取 p = 0.6 如何设计随机实验

X P

0 1–p

1 p16/18