二、随机变量及其分布 主讲教师:董庆宽 副教授 研究方向:密码学与信息安全 电子邮件:qkdong@xidian.edu.cn 基于MATLAB的概率统计数值实验 二、随机变量及其分布 主讲教师:董庆宽 副教授 研究方向:密码学与信息安全 电子邮件:qkdong@xidian.edu.cn 个人主页:http://web.xidian.edu.cn/qkdong/
内容介绍 二、随机变量及其分布 1. MATLAB中概率分布函数 2. 二项分布实验 3. 泊松分布实验 4. 二项分布与泊松分布关系实验 5. 连续型随机变量分布实验 6. 随机变量的均值与方差 7. 逆累积分布函数实验 8. 中心极限定理实验
1. MATLAB中概率分布函数 MATLAB为常见自然概率分布提供了下列5类函数 ①概率密度函数(pdf),求随机变量X在x点处的概率密度值 ②累积分布函数(cdf),求随机变量X在x点处的分布函数值 ③逆累积分布函数(inv),求随机变量X在概率点处的分布函数反函数值 ④均值与方差计算函数(stat),求给定分布的随机变量X的数学期望E(X)和方差var(X) ⑤随机数生成函数(rnd),模拟生成指定分布的样本数据(调用格式:x=分布rnd(分布参数),如x=normrnd(0,1))
1. MATLAB中概率分布函数 常见的分布类型名如下 分布类型 MATLAB名称 正态分布 norm 二项分布 bino 指数分布 exp Poisson分布 poiss 均匀分布 unif 几何分布 geo β分布 beta 超几何分布 hyge Γ分布 gam 离散均匀分布 unid 对数正态分布 logn 连续均匀分布 rayleigh分布 rayl 负二项分布 nbin weibull 分布 weib 2分布 chi2 F分布 f 学生氏t分布 t
1. MATLAB中概率分布函数 具体函数的命名规则是: 函数名=分布类型名称+函数类型名称(pdf、cdf、inv、stat、rnd) 例如,normpdf、normcdf、norminv、normstat和normrnd分别是正态分布的概率密度、累积分布、逆累积分布、数字特征和随机数生成函数。 关于这5类函数的语法,请详见有关书籍 快捷的学习可借助MATLAB的系统帮助,通过指令doc获得具体函数的详细信息,语法是 doc <函数名>
2. 二项分布实验 已知Y~b(20, 0.2)求Y分布率的值,并划出图形 在Matlab中输入以下命令: binopdf(10,20,0.2) x=0:1:20; y=binopdf(x,20,0.2) plot(x,y, ‘r.’) 结果: ans = 0.0020 y =0.0115 0.0576 0.1369 0.2054 0.2182 0.1746 0.1091 0.0545 0.0222 0.0074 0.0020 0.0005 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2. 二项分布实验 已知Y~b(20, 0.3)求Y分布函数的值,画出函数图像 在Matlab中输入以下命令: binocdf(10,20,0.3) x=0:1:20; y=binocdf(x,20,0.3) ezplot('binocdf(t,20,0.3)',[0,20]) 结果: ans = 0.9994 y = 0.0115 0.0692 0.2061 0.4114 0.6296 0.8042 0.9133 0.9679 0.9900 0.9974 0.9994 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
2. 二项分布实验
2. 二项分布实验 到某服务机构办事总是要排队等待的。设等待时间T是服从指数分布的随机变量(单位:分钟),概率密度为 设某人一个月内要到此办事10次,若等待时间超过15分钟,他就离去。求: (1)恰好有两次离去的概率; (2)最多有两次离去的概率; (3)至少有两次离去的概率; (4)离去的次数占多数的概率。
2. 二项分布实验 解 首先求任一次离去的概率,依题意 设10次中离去的次数为X,则X~b(10, p) 解 首先求任一次离去的概率,依题意 设10次中离去的次数为X,则X~b(10, p) >> p=1-expcdf(15,10) %任一次离去的概率 p1=binopdf(2,10,p) %恰有两次离去的概率 q=binopdf([0:2],10,p);p2=sum(q) %最多有两次离去的概率 q=binopdf([0:1],10,p);p3=1-sum(q) %最少有两次离去的概率 q=binopdf([0:5],10,p);p4=1-sum(q) %离去的次数占多数的概率 p = 0.2231 p1 = 0.2972 p2 = 0.6073 p3 = 0.6899 p4 = 0.0112
3. 泊松分布实验 假设电话交换台每小时接到的呼叫次数X服从参数=3的泊松分布,求 在Matlab中输入以下命令: (1) 每小时恰有4次呼叫的概率 (2) 一小时内呼叫不超过5次的概率 (3) 画出分布律图像 在Matlab中输入以下命令: (1)p1= poisspdf(4,3) (2)p2= poisscdf(5,3) (3)x=0:1:20;y=poisspdf(x,3);plot(x,y)
3. 泊松分布实验
4. 二项分布与泊松分布关系实验 二项分布与泊松分布的关系 例7:X~b(200,0.02),Y 服从参数为4的泊松分布,划出分布率图像 y1=binopdf(x,200,0.02); y2=poisspdf(x,4); plot(x,y1,’r.’,x,y2,’b.’)
4. 二项分布与泊松分布关系实验 泊松定理 设λ>0是一个常数,n是任意正整数,设npn=λ,则对于任意固定的非负整数k,有 (用泊松分布来逼近二项分布的定理) 设λ>0是一个常数,n是任意正整数,设npn=λ,则对于任意固定的非负整数k,有 例9 某种重大疾病的医疗险种,每份每年需交保险费100元,若在这一年中,投保人得了这种疾病,则每份可以得到索赔额10000元,假设该地区这种疾病的患病率为0.0002,现该险种共有10000份保单,问: (1)保险公司亏本的概率是多少? (2)保险公司获利不少于80万元的概率是多少?
解 设 表示这一年中发生索赔的份数,依题意, 的统计规律可用二项分布 来描述。由二项分布与泊松分布的近似计算关系有 近似服从参数为2的泊松分布。 当索赔份数超过100份时,则保险公司发生亏本,亏本的概率为 当索赔份数不超过20份时,则保险公司获利就不少于80万元,其概率为
>> [p]=poisspdf([0:100],2);%计算101个泊松分布概率值 或 [p]=binopdf([0:100],10000,0.0002); %按二项分布计算 p1=1-sum(p) %求出保险公司亏本的概率 p1 = 0.0000 >> [p]=poisspdf([0:19],2);%计算出20个泊松分布概率值 或 [p]=binopdf([0:19],10000,0.0002); %按二项分布计算 p2=sum(p) %求出保险公司获利不少于80万元的概率 p2 = 1.0000
5. 连续型随机变量分布实验 离散均匀分布的概率密度函数和累积分布函数 unidpdf(X,N) unidcdf(X,N) 随机变量X在1到N上的N各自然数之间等可能取值 在Matlab中输入以下命令: x=1:1:10; y=unidpdf(x,10) 结果:y = 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 x=0:1:10; y=unidcdf(x,10) 结果:y = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
5. 连续型随机变量分布实验 连续均匀分布 例: 画出均匀分布U(2,5)的概率密度函数和分布函数的图形. 在Matlab中输入以下命令: 密度函数:f=unifpdf(x,a,b) 分布函数:f=unifcdf(x,a,b) 例: 画出均匀分布U(2,5)的概率密度函数和分布函数的图形. 在Matlab中输入以下命令: x=0:0.01:7; y=unifpdf(x, 2, 5); z=unifcdf(x, 2, 5); plot(x,y,x,z)
5. 连续型随机变量分布实验 (2) 指数分布 密度函数:f=exppdf(x,) 分布函数:F=expcdf(x,) 例: 画出指数分布E(1)的概率密度函数和分布函数的图形. 求P(0<X<5) P(0<X<20). 在Matlab中输入以下命令: x=0:0.1:5; y=exppdf(x,2); z=expcdf(x,2); plot(x,y,x,z) result1=expcdf(5,2)-expcdf(0,2) result2=expcdf(20,2)-expcdf(0,2)
结果:result1 = 0.91791500137610 result2 = 0.99995460007024
5. 连续型随机变量分布实验 (3) 正态分布 密度函数:f=normpdf(x,,) 分布函数:F=normcdf(x,,) 例: 画出正态分布N(1,4)的概率密度函数和分布函数的图形. 求P(1<X<6). 在Matlab中输入以下命令: x=-5:0.1:6; y=normpdf(x,1,2); z=normcdf(x,1,2); plot(x,y,x,z) result=normcdf(6,1,2)-normcdf(1,1,2)
5. 连续型随机变量分布实验 结果:Result =0.4938
5. 连续型随机变量分布实验 例11:在同一坐标下,画下列正态分布的密度函数图像 (1) μ=3, σ=0.5, 0.7, 1, 1.5, 2 (2) σ=0.5, μ=1,2,3,4 (1)命令: x=-6:0.1:6; y1=normpdf(x,3,0.5); y2=normpdf(x,3,0.7); y3=normpdf(x,3,1); y4=normpdf(x,3,1.5); y5=normpdf(x,3,2); plot(x,y1,'.',x,y2,'+',x,y3,'*',x,y4,'d',x,y5)
例 观察正态分布参数对密度曲线的影响。 解:>> clear mu1=2.5;mu2=3;sigma1=0.5;sigma2=0.6; x=(mu2-4*sigma2):0.01:(mu2+4*sigma2); y1=normpdf(x,mu1,sigma1); %考察均值的影响 y2=normpdf(x,mu2,sigma1); y3=normpdf(x,mu1,sigma1); %考察方差的影响 y4=normpdf(x,mu1,sigma2); subplot(1,2,1) %考察结果的可视化 plot(x,y1,'-g',x,y2,'-b') xlabel('\fontsize{12}μ1<μ2,σ1=σ2' ) legend('μ1','μ2') subplot(1,2,2) plot(x,y3,'-g',x,y4,'-b') xlabel('\fontsize{12}μ1=μ2,σ1<σ2' ) legend('σ1','σ2')
5. 连续型随机变量分布实验 计算正态分布的累积概率值 例,设X~N(4,32), P{3<X<6}, P{X>3} 解: 调用函数normcdf(x,μ,σ) 返回函数值 解: >> p1=normcdf(6,4,3)-normcdf(3,4,3) p1 = 0.3781 >> p2=1-normcdf(3,4,3) p2 = 0.6306
例 正态分布参数μ和σ对变量x取值规律的约束——3σ准则。 解:>> clear,clf %(标准)正态分布密度曲线下的面积 X=linspace(-5,5,100); Y=normpdf(X,0,1); yy=normpdf([-3,-2,-1,0,1,2,3],0,1); plot(X,Y,'k-',[0,0],[0,yy(4)],'c-.') hold on plot([-2,-2],[0,yy(2)],'m:',[2,2],[0,yy(6)],'m:',[-2,-0.5],[yy(6),yy(6)],'m:',[0.5,2],[yy(6),yy(6)],'m:') plot([-1,-1],[0,yy(3)],'g:',[1,1],[0,yy(5)],'g:',[-1,-0.5],[yy(5), yy(5)],'g:',[0.5,1],[yy(5),yy(5)],'g:') plot([-3,-3],[0,yy(1)],'b:',[3,3],[0,yy(7)],'b:',[-3,-0.5],[yy(7), yy(7)],'b:',[0.5,3],[yy(7),yy(7)],'b:')
5. 连续型随机变量分布实验 hold off text(-0.5,yy(6)+0.005,'\fontsize{14}95.44%') text(-3.2,-0.03,'\fontsize{10}μ-3σ') text(-2.2,-0.03,'\fontsize{10}μ-2σ') text(-1.2,-0.03,'\fontsize{10}μ-σ') text(-0.05,-0.03,'\fontsize{10}μ') text(0.8,-0.03,'\fontsize{10}μ+σ') text(1.8,-0.03,'\fontsize{10}μ+2σ') text(2.8,-0.03,'\fontsize{10}μ+3σ')
6. 随机变量的数学期望和方差 对于任意的分布,可用Matlab中的函数和运算编程实现 对于给定的分布,只需给出分布的参数,即可调用stat族函数,得出数学期望和方差,调用格式 [E,D]=分布+stat(参数) 例:求二项分布参数n=100,p=0.2的数学期望和方差: 解:>>n=100; >>p=0.2; >>[E,D]=binostat(n,p); 结果显示:E= 20 D= 16
例 绘制正态分布的密度函数、分布函数曲线,并求均值与方差 6. 随机变量的数学期望和方差 例 绘制正态分布的密度函数、分布函数曲线,并求均值与方差 解: >> clear mu=2.5;sigma=0.6; x=(mu-4*sigma):0.005:(mu+4*sigma); y=normpdf(x,mu,sigma); f=normcdf(x,mu,sigma); plot(x,y,'-g',x,f,':b') [M,V]=normstat(mu,sigma) legend('pdf','cdf',-1)
M=2.5000 V=0.3600 从图中可以看出,正态密度曲线是关于x=μ对称的钟形曲线(两侧在μ±σ处各有一个拐点),正态累积分布曲线当x=μ时F(x)=0.5。
7. 逆累积分布函数 逆累积分布函数就是返回给定概率条件下的自变量的临界值,实际上是分布函数的逆函数。 icdf(Inverse Cumulative Distribution Function) 即:在分布函数F(x)=p中已知p求其相对应的x的值 调用:在分布函数名后加inv 如:X=norminv(p,mu,sgm) 也有2)X=icdf(‘name’,p,A1,A2,A3),其中name为相应的函数名,如‘normal’;p为给定的概率值; A1,A2,A3为相应的参数
7. 逆累积分布函数 例、计算标准正态分布N(0,1)概率值0.1,0.3, 0.5,0.7,0.9,所对应的x的值 命令: 结果: y=0.1:0.2:0.9; x=norminv(y,0,1) 结果: x=-1.2816 -0.5244 0 0.5244 1.2816 检验:y1=normcdf(x,0,1); y1=0.1000 0.3000 0.5000 0.7000 0.9000
7. 逆累积分布函数 例、计算二项分布b(10,0.5)概率值0.1,0.3, 0.5,0.7,0.9,所对应的x的值 命令: p=0.1:0.2:0.9; x=binoinv(p,10,0.5) 结果: x=3 4 5 6 7 检验:y1=binocdf(x,10,0.5); 结果: y1=0.1719 0.3770 0.6230 0.8281 0.9453
7. 逆累积分布函数 在离散分布情形下,icdf 返回使cdf(x)p的第一个值x 上例中,对p=0.1,对应cdf(x)0.1的第一个值为3,故返回值为3 B(10,0.5)的分布函数图像
7.逆累积分布函数-上分位点 定义:上 分位点:设随机变量X的分布函数为: F(x),如果实数 满足P(X> )= ,则称 为上 例14、计算自由度为8的卡方分布的上 分位点, 其中α=0.1,0.05,0.025 命令: x=[0.1,0.05,0.025]; y=chi2inv(1-x,8) 结果: y=13.3616 15.5073 17.5345
例 标准正态分布α分位数的概念图示。 解 >> %α分位数示意图(标准正态分布,α=0.05) clear,clf 例 标准正态分布α分位数的概念图示。 解 >> %α分位数示意图(标准正态分布,α=0.05) clear,clf data=normrnd(0,1,300,1); xalpha1=norminv(0.05,0,1); xalpha2=norminv(0.95,0,1); xalpha3=norminv(0.025,0,1); xalpha4=norminv(0.975,0,1); subplot(3,1,1) capaplot(data,[-inf,xalpha1]);axis([-3,3,0,0.45]) subplot(3,1,2) capaplot(data,[xalpha2,inf]);axis([-3,3,0,0.45]) subplot(3,1,3) capaplot(data,[-inf,xalpha3]);axis([-3,3,0,0.45]) hold on capaplot(data,[xalpha4,inf]);axis([-3,3,0,0.45]) hold off xalpha1 xalpha2 xalpha3 xalpha4
xalpha1 = -1.6449 xalpha2 = 1.6449 xalpha3 = -1.9600 xalpha4 = 1.9600
8. 中心极限定理 例1利用随机数样本验证中心极限定理 独立同分布的随机变量的和的极限分布服从正态分布,通过产生容量为n的poiss分布和exp分布的样本,研究其和的渐近分布。 算法如下: ① 产生容量为n的独立同分布的随机数样本,得其均值和标准差; ② 将随机数样本和标准化; ③ 重复①、②; ④ 验证所得标准化的随机数样本和是否服从标准正态分布
8. 中心极限定理 >> clear n=2000; means=0; s=0; y=[]; lamda=4; a=lamda; for i=1:n r=poissrnd(a,n,1);%可换成r=exprnd(a,n,1); means=mean(r);%计算样本均值 s=std(r);%计算样本标准差 y(i)=sqrt(n).*(means-a)./sqrt(s); end normplot(y);%分布的正态性检验 title('poiss分布,中心极限定理')
8. 中心极限定理 棣莫弗-拉普拉斯定理的应用 Galton钉板模型和二项分布 Galton钉板试验是由英国生物统计学家和人类学家Galton设计的。故而得名。 通过模拟Calton钉板试验,观察和体会二项分布概率分布列的意义、形象地理解De Moivre -Laplace中心极限定理。
8. 中心极限定理 高尔顿钉板试验 高尔顿( Francis Galton,1822-1911) 英国人类学家和气象学家 落下的位置为 15层中向右的次数减向左的次数 O x -8 -7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 小球最后落入的格数 ? 共15层小钉 符号函数,大于0返回1,小于0返回-1,等于0返回0 记小球向右落下的次数为 则 记小球向左落下的次数为 则 W取值从-8到8
8. 中心极限定理 高尔顿钉板试验 什么曲线 ? 共15层小钉 记 则 小球碰第 层钉后向右落下 小球碰第 层钉后向左落下 近似 x O -8 -7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 共15层小钉 ? 什么曲线 O 小球碰第 层钉后向右落下 小球碰第 层钉后向左落下 记 近似 则
8. 中心极限定理 模拟Galton钉板试验的步骤: (1) 确定钉子的位置:将钉子的横、纵坐标存储在两个矩阵X和Y中。 (2) 在Galton钉板试验中,小球每碰到钉子下落时都具有两种可能性,设向右的概率为p,向左的概率为q=1-p,这里p=0.5,表示向左向右的机会是相同的。
8. 中心极限定理 模拟过程如下:首先产生一均匀随机数u,这只需调用随机数发生器指令rand(m,n)。 rand(m,n)指令:用来产生m×n个(0,1)区间中的随机数,并将这些随机数存于一个m×n矩阵中,每次调用rand(m,n)的结果都会不同。如果想保持结果一致,可与rand(‘seed’,s)配合使用,这里s是一个正整数,例如 >> rand('seed',1),u=rand(1,6) u = 0.5129 0.4605 0.3504 0.0950 0.4337 0.7092 而且再次运行该指令时结果保持不变。除非重设种子seed的值,如 >> rand('seed',2),u=rand(1,6) u = 0.0258 0.9210 0.7008 0.1901 0.8673 0.4185 这样结果才会产生变化。
8. 中心极限定理 将[0,1]区间分成两段,区间[0,p)和[p,1]。如果随机数u属于[0,p),让小球向右落下;若u属于[p,1] ,让小球向左落下。将这一过程重复n次,并用直线连接小球落下时所经过的点,这样就模拟了小球从顶端随机地落人某一格子的过程。 (3) 模拟小球堆积的形状。输入扔球次数m(例如m=50、100、500等等),计算落在第i个格子的小球数在总球数m中所占的比例,这样当模拟结束时,就得到了频率 用频率反映小球的堆积形状
8. 中心极限定理 (4)用如下动画指令制作动画: movien(n):创建动画矩阵;制作动画矩阵数据; Getframe:拷贝动画矩阵; movie(Mat, m):播放动画矩阵m次。 M文件如下:
8. 中心极限定理 解:>> clear,clf,m=100;n=5;y0=2;%设置参数 ballnum=zeros(1,n+1); p=0.5;q=1-p; for i=n+1:-1:1 %创建钉子的坐标x,y x(i,1)=0.5*(n-i+1); y(i,1)=(n-i+1)+y0; for j=2:i x(i,j)=x(i,1)+(j-1)*1; y(i,j)=y(i,1); end mm=moviein(m); %动画开始,模拟小球下落路径 for i=1:m s=rand(1,n); %产生n个随机数 xi=x(1,1);yi=y(1,1);k=1;l=1; %小球遇到第一个钉子 for j=1:n plot(x(1:n,:),y(1:n,:),‘o’,x(n+1,:),y(n+1,:),‘.-’),%画钉子的位置 axis([-2 n+2 0 y0+n+1]),hold on
k=k+1; %小球下落一格 if s(j)>p l=l+0;%小球左移 else l=l+1;%小球右移 end xt=x(k,l);yt=y(k,l);%小球下落点的坐标 h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1]) %画小球运动轨迹 xi=xt;yi=yt; ballnum(l)=ballnum(l)+1; %计数 ballnum1=3*ballnum./m; bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1]) %画各格子的频率 mm(i)=getframe; %存储动画数据 hold off movie(mm,1) %播放动画一次
……
作业 1. 已知二项分布Xb(15,0.2) 求: (1) 分布率和分布函数值,并画出曲线 (2) 求该分布的数学期望和方差 (4) 验证50个该分布的随机样本的和的标准化变量服从标准正态分布 2. 对于正态分布: (1)给出N(2,9)的密度曲线和分布函数曲线 (2)在同一图上划出均值为2,标准差为0.5,0.7,1,2的密度曲线 (3)已知=0.05,求N(0,1)的上分位点.
谢谢!