概率统计的数值实验 ——Matlab.ppt

概率统计的数值实验MATLAB在概率统计教学中的应用,崔明涛2012年10月11日,引言,而MATLAB 软件具有简单易学、易操作和绘图功能强等特点， 利用MATLAB 软件的图形可视功能将概率统计的内容用图形表示出来，通过图形让学生加深理解，以达到事半功倍的效果。,概率论与数理统计知识比较抽象，逻辑性较强。因此，建议让学生结合理论和公式推导，进行数值试验和相关调查，直观地感受数学概念和理论，从而提高学生解决实际问题的信心和能力。,概率论,1.rand(m,n)：生成m×n的随机矩阵，每个元素都在(0,1) 间，生成方式为均匀分布。2.randn(m,n)：生成m×n的随机矩阵，每个元素都在(0,1) 间，生成方式为正态分布。3.randperm(m)：生成一个1m的随机整数排列。4.perms(1:n)：生成一个1n的全排列，共n!个。5.取整函数系列：（1）fix(x)：截尾法取整；（2）floor(x)：退一法取整（不超过x的最大整数）；（3）ceil(x)：进一法取整（= floor(x)+1）；（4）round(x)：四舍五入法取整。6.unique(a)：合并a中相同的项。7.prod(x)：向量x的所有分量元素的积。,一、MATLAB常用的与随机数产生相关的函数：,示例：>> rand(1) %生成一个(0,1)间的随机数 ans = 0.8147>> rand(2,2) %生成一个2×2阶(0,1)间的随机数矩阵ans = 0.9134 0.0975 0.6324 0.2785>> randperm(5) %生成一个15的随机整数排列ans = 4 1 5 2 3>> a=1 2 4 2 3 3 2;unique(a)ans = 1 2 3 4,例1 随机投掷均匀硬币，观察国徽朝上与国徽 朝下的频率。,解>> n= 3000100000000;m=0; for i=1:n t=randperm(2);%生成一个12的随机整数排列 x=t-1;%生成一个01的随机整数排列 y=x(1); if y=0; m=m+1; endendp1=m/np2=1-p1,可见当 时，,解 记事件 为第i个人拿到自已枪，事件 为第i个人没拿到自己枪，易知：又记 为没有一个人拿到自己枪的概率。有乘法公式可知：,例2 某班有n个人，每人各有一支枪，这些枪外形一样。某次夜间紧急集合，若每人随机地取走一支枪，问没有一个人拿到自己枪的概率是多少？,于是 所以特别地，当n较大时， 。 因此，可随机模拟出没有人拿到自己枪的频率，根据频率的稳定性，近似当做概率，然后去估计自然常数e。算法如下：,1、产生n个随机数的随机序列；2、检验随机列与自然列是否至少有一个配对；3、对没有一个配对的序列进行累积 p；4、重复1、2、3步 m 次；5、估计 。,具体程序及相关结果为（注：自然常数 e 2.7183）：,>> m=40000;n=50;p=0;for j=1:m k=0; sui=randperm(n); for i=1:n if sui(i)=i k=k+1; else k=k; end end if k=0 p=p+1; else p=p; endende=m/pe = 2.7313,设针与平行线的夹角为 ，针的中心与最近直线的距离为 。针与平行线相交的充要条件是 ，则所求概率为故可得 的近似计算公式 ，其中n为随机试验次数，m为针与平行线相交的次数。,例3 Buffon投针实验 在画有许多间距为 的等距平行线的白纸上，随机投掷一根长为 的均匀直针，求针与平行线相交的概率，并计算 的近似值。,解 >> clear,clfn=10000000;l=0.5;m=0;d=1;for i=1:n x=l/2*sin(rand(1)*pi);y=rand(1)*d/2; if x>=y m=m+1; endendp1=m/npai=2*n*l/(m*d),例4 在100个人的团体中，不考虑年龄差异，研究是否有两个以上的人生日相同。假设每人的生日在一年365天中的任意一天是等可能的，那么随机找n个人(不超过365人)。 (1)求这n个人生日各不相同的概率是多少？从而求这n个人中至少有两个人生日相同这一随机事件发生的概率是多少？ (2)近似计算在30名学生的一个班中至少有两个人生日相同的概率是多少？,解: (1)>> clear,clffor n=1:100 p0(n)=prod(365:-1:365-n+1)/365n; p1(n)=1-p0(n);endp1=ones(1,100)-p0;n=1:100;plot(n,p0,n,p1,'-')xlabel('人数'),ylabel('概率')legend('生日各不相同的概率','至少两人生日相同的概率')axis(0 100 -0.1 1.199),grid on,p1(30)=0.7063, p1(60)= 0.9941,分析：在30名学生中至少两人生日相同的概率为70.63。下面进行计算机仿真。 随机产生30个正整数，代表一个班30名学生的生日，然后观察是否有两人以上生日相同。当30个人中有两人生日相同时，输出“1”，否则输出“0”。如此重复观察100次，计算出这一事件发生的频率 。,(2)>> clear,clfn=0;for m=1:100 %做100次随机试验 y=0; x=1+fix(365*rand(1,30); %产生30个随机数 for i=1:29 %用二重循环寻找30个随机数 中是否有相同数 for j=i+1:30 if x(i)=x(j) y=1;break; end end end n=n+y; %累计有两人生日相同的试验次数endf=n/m %计算频率,f =0.6900f =0.7900f =0.6700f =0.7300f =0.7500f =0.6900f =0.7200f =0.6700f =0.6800,重复观察，数据如下：,例5 Galton钉板模型和二项分布 Galton钉板试验是由英国生物统计学家和人类学家Galton设计的。故而得名。 通过模拟Calton钉板试验，观察和体会二项分布概率分布列的意义、形象地理解De Moivre -Laplace中心极限定理。,共15层小钉,高尔顿钉板试验,小球最后落入的格数 ?,记小球向右落下的次数为 则,记小球向左落下的次数为 则,符号函数,大于0返回1,小于0返回-1,等于0返回0,高尔顿( Francis Galton,1822-1911) 英国人类学家和气象学家,O,记,则,近似,高尔顿钉板试验,共15层小钉,模拟Galton钉板试验的步骤： (1) 确定钉子的位置：将钉子的横、纵坐标存储在两个矩阵X和Y中。 (2) 在Galton钉板试验中，小球每碰到钉子下落时都具有两种可能性，设向右的概率为p，向左的概率为q1-p，这里p=0.5，表示向左向右的机会是相同的。 模拟过程如下：首先产生一均匀随机数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这样结果才会产生变化。,将0,1区间分成两段，区间0,p)和p,1。如果随机数u属于0,p)，让小球向右落下；若u属于p,1 ，让小球向左落下。将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落人某一格子的过程。 (3) 模拟小球堆积的形状。输入扔球次数m(例如m50、 100、500等等)，计算落在第i个格子的小球数在总球数m中所占的比例，这样当模拟结束时，就得到了频率 用频率反映小球的堆积形状。 (4)用如下动画指令制作动画： movien(n)：创建动画矩阵；制作动画矩阵数据； Getframe：拷贝动画矩阵； movie(Mat, m)：播放动画矩阵m次。 M文件如下：,解：>> 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:ix(i,j)=x(i,1)+(j-1)*1; y(i,j)=y(i,1); endendmm=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:nplot(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; end 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 offendmovie(mm,1) %播放动画一次,概率密度函数（pdf），求随机变量X在x点处的概率密度值累积分布函数（cdf），求随机变量X在x点处的分布函数值逆累积分布函数（inv），求随机变量X在概率点 处的分 布函数反函数值均值与方差计算函数（stat），求给定分布的随机变量X的 数学期望E（X）和方差var（X）。随机数生成函数（rnd），模拟生成指定分布的样本数据。,二、MATLAB为常见自然概率分布提供了下列5类函数：,具体函数的命名规则是：函数名分布类型名称+函数类型名称(pdf、cdf、inv、stat、rnd) 其中，分布类型名称如下： 分布类型 MATLAB名称,正态分布 norm指数分布 exp均匀分布 unif 分布 beta 分布 gam对数正态分布 lognrayleigh分布 raylweibull 分布 weib二项分布 binoPoisson分布 poiss几何分布 geo超几何分布 hyge离散均匀分布 unid负二项分布 nbin,例如，normpdf、normcdf、norminv、normstat和normrnd分别是正态分布的概率密度、累积分布、逆累积分布、数字特征和随机数生成函数。,关于这5类函数的语法，请详见有关书籍。 快捷的学习可借助MATLAB的系统帮助，通过指令doc获得具体函数的详细信息，语法是 doc ,例6 到某服务机构办事总是要排队等待的。设等待时间T是服从指数分布的随机变量(单位：分钟)，概率密度为设某人一个月内要到此办事10次，若等待时间超过15分钟，他就离去。求： (1)恰好有两次离去的概率； (2)最多有两次离去的概率； (3)至少有两次离去的概率； (4)离去的次数占多数的概率。,解 首先求任一次离去的概率，依题意 设10次中离去的次数为X，则 。,>> 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.2231p1 = 0.2972p2 = 0.6073p3 = 0.6899p4 = 0.0112,例7 某一急救中心在长度为t的时间间隔内收到的紧急呼救次数服从参数为t2的泊松分布，而与时间间隔的起点无关(时间以小时计)，求： (1) 在某一天中午12时至下午3时没有收到紧急呼救的概率； (2) 某一天中午12时至下午5时至少收到1次紧急呼救的概率。,（1）>> P1=poisscdf(0,3/2)P1 =0.2231或者>> P1=poisspdf(0,3/2)P1 =0.2231中午12时到下午3时没有收到紧急呼救的概率为0.2231。 （2）>> P2=1-poisscdf(0,5/2)P2 = 0.9179中午12时至下午5时至少收到1次紧急呼救的概率为0.9179。,解 本题计算需调用函数poisscdf，其格式为poisscdf(x,)， 返回 的值。,例8 某厂研发了一种新产品，现要设计它的包装箱，要求每箱至少装100件产品，且开箱验货时，每箱至少装有100件合格产品的概率不应小于0.9，假设随机装箱时每箱中的不合格产品数服从参数为3的泊松分布。 问：要设计的这种包装箱，每箱至少应装多少件产品才能满足要求？,解 设每箱至少装100+m件产品，X表示每箱中的不合格品数，则X 服从参数为3的泊松分布，即 ，依题意，即要求按下面的不等式确定m。,>> clear;clf,m=0;p=0;while p0.9。即设计的包装箱每箱至少应装106件产品。,例9 某种重大疾病的医疗险种，每份每年需交保险费100元，若在这一年中，投保人得了这种疾病，则每份可以得到索赔额10000元，假设该地区这种疾病的患病率为0.0002，现该险种共有10000份保单，问：(1)保险公司亏本的概率是多少?(2)保险公司获利不少于80万元的概率是多少?,解 设 表示这一年中发生索赔的份数，依题意， 的统计规律可用二项分布 来描述。由二项分布与泊松分布的近似计算关系有 近似服从参数为2的泊松分布。 当索赔份数超过100份时，则保险公司发生亏本，亏本的概率为 当索赔份数不超过20份时，则保险公司获利就不少于80万元，其概率为,>> p=poisspdf(0:19,2);%计算出20个泊松分布概率值 或 p=binopdf(0:19,10000,0.0002); %按二项分布计算 p2=sum(p) %求出保险公司获利不少于80万元的概率 p2 = 1.0000,>> p=poisspdf(0:100,2);%计算101个泊松分布概率值或 p=binopdf(0:100,10000,0.0002); %按二项分布计算 p1=1-sum(p) %求出保险公司亏本的概率 p1 = 0.0000,例10 设 ，求 ， 。,本题计算正态分布的累积概率值，调用函数normcdf，其格式为normcdf(x,)，返回 的值。解：>> p1=normcdf(6,4,3)-normcdf(3,4,3)p1 = 0.3781>> p2=1-normcdf(3,4,3)p2 = 0.6306,例11 绘制正态分布的密度函数、分布函数曲线，并求均值与方差。,解：>> clearmu=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.5000V=0.3600,从图中可以看出，正态密度曲线是关于x对称的钟形曲线(两侧在±处各有一个拐点)，正态累积分布曲线当x时F(x)0.5。,例12 观察正态分布参数对密度曲线的影响。,解：>> clearmu1=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('fontsize121<2,1=2' )legend('1','2')subplot(1,2,2)plot(x,y3,'-g',x,y4,'-b')xlabel('fontsize121=2,1<2' )legend('1','2'),例13 正态分布参数和对变量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 onplot(-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:'),hold offtext(-0.5,yy(6)+0.005,'fontsize1495.44%')text(-0.5,yy(5)+0.005,'fontsize1468.26%')text(-0.5,yy(7)+0.005,'fontsize1499.74%')text(-3.2,-0.03,'fontsize10-3')text(-2.2,-0.03,'fontsize10-2')text(-1.2,-0.03,'fontsize10-')text(-0.05,-0.03,'fontsize10')text(0.8,-0.03,'fontsize10+')text(1.8,-0.03,'fontsize10+2')text(2.8,-0.03,'fontsize10+3'),例14 标准正态分布分位数的概念图示。,解 >> %分位数示意图（标准正态分布，=0.05）clear,clfdata=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 oncapaplot(data,xalpha4,inf);axis(-3,3,0,0.45)hold offxalpha1 xalpha2 xalpha3 xalpha4,xalpha1 = -1.6449xalpha2 = 1.6449xalpha3 = -1.9600 xalpha4 = 1.9600,数理统计基础,Matlab统计工具箱中常见的统计命令,1、基本统计量对于随机变量x，计算其基本统计量的命令如下：均值：mean(x) 标准差：std(x)中位数：median(x) 方差：var(x)偏度：skewness(x) 峰度：kurtosis(x)2、频数直方图的描绘A、给出数组data的频数表的命令为： N,X=hist(data,k) 此命令将区间min(data),max(data)分为k个小区间(缺省为10)，返回数组data落在每一个小区间的频数N和每一个小区间的中点X。B、描绘数组data的频数直方图的命令为： hist(data,k),3、参数估计A、对于正态总体，点估计和区间估计可同时由以下命令获得：muhat,sigmahat,muci,sigmaci=normfit(x,alpha)此命令在显著性水平alpha下估计x的参数（alpha缺省值为5%），返回值muhat是均值的点估计值，sigmahat是标准差的点估计值，muci是均值的区间估计，sigmaci是标准差的区间估计。B、对其他分布总体，两种处理办法：一是取容量充分大的样本，按中心极限定理，它近似服从正态分布，仍可用上面估计公式计算；二是使用特定分布总体的估计命令，常用的命令如：muhat,muci=expfit(x,alpha)lambdahat, lambdaci=poissfit(x,alpha)phat, pci=weibfit(x,alpha),4、正态总体假设检验A、单总体均值的z检验： h,sig,ci=ztest(x,m,sigma,alpha,tail)检验数据x关于总体均值的某一假设是否成立，其中sigma为已知方差，alpha为显著性水平，究竟检验什么假设取决于tail的取值：tail=0，检验假设“x的均值等于m”tail=1，检验假设“x的均值大于m”tail=-1，检验假设“x的均值小于m”tail的缺省值为0， alpha的缺省值为5%。返回值h为一个布尔值，h=1表示可拒绝原假设， h=0表示不可拒绝原假设，sig为假设成立的概率，ci为均值的1- alpha置信区间。B、单总体均值的t检验： h,sig,ci=ttest(x,m,alpha,tail)C、双总体均值的t检验： h,sig,ci=ttest2(x,y,alpha,tail),5、非参数检验：总体分布的检验Matlab统计工具箱提供了两个对总体分布进行检验的命令：A、 h=normplot(x)此命令显示数据矩阵x的正态概率图，如果数据来自于正态分布，则图形显示出直线形态，而其他概率分布函数显示出曲线形态。B、h=weibplot(x)此命令显示数据矩阵x的Weibull概率图，如果数据来自于Weibull分布，则图形显示出直线形态，而其他概率分布函数显示出曲线形态。,例15 一道工序用自动化车床连续加工某种零件，由于刀具损坏等会出现故障。故障是完全随机的，并假定生产任一零件时出现故障机会均相同，工作人员是通过检查零件来确定工序是否出现故障的。现积累有100次故障纪录，故障出现时该刀具完成的零件数如下：459 362 624 542 509 584 433 748 815 505 612 452 434 982 640 742 565 706 593 680 926 653 164 487 734 608 428 1153 593 844 527 552 513 781 474 388 824 538 862 659 775 859 755 49 697 515 628 954 771 609 402 960 885 610 292 837 473 677 358 638 699 634 555 570 84 416 606 1062 484 120 447 654 564 339 280 246 687 539 790 581 621 724 531 512 577 496 468 499 544 645 764 558 378 765 666 763 217 715 310 851试观察该刀具出现故障时完成的零件数属于哪种分布？,>> %数据输入x1=459 362 624 542 509 584 433 748 815 505;x2=612 452 434 982 640 742 565 706 593 680;x3=926 653 164 487 734 608 428 1153 593 844;x4=527 552 513 781 474 388 824 538 862 659;x5=775 859 755 49 697 515 628 954 771 609;x6=402 960 885 610 292 837 473 677 358 638;x7=699 634 555 570 84 416 606 1062 484 120;x8=447 654 564 339 280 246 687 539 790 581;x9=621 724 531 512 577 496 468 499 544 645;x10=764 558 378 765 666 763 217 715 310 851;x=x1 x2 x3 x4 x5 x6 x7 x8 x9 x10;%作频数直方图hist(x,10)N,X=hist(x,10)%分布的正态性检验normplot(x),N = 3 3 7 14 24 22 14 8 3 2X = 1.0e+003 * 0.1042 0.2146 0.3250 0.4354 0.5458 0.6562 0.7666 0.8770 0.9874 1.0978,>> %参数估计muhat,sigmahat,muci,sigmaci=normfit(x)muhat = 594sigmahat = 204.1301muci = 553.4962 634.5038sigmaci = 179.2276 237.1329刀具寿命服从正态分布，均值估计值为594，方差估计值为204.1301，均值的95%置信区间为553.4962，634.5038，方差的95%置信区间为179.2276，237.1329,>> %假设检验h,sig,ci=ttest(x,594)%已知刀具寿命服从正态分布，方差未知的情况下，检验寿命均值是否等于594。h = 0sig = 1ci = 553.4962 634.5038检验结果：布尔变量h=0，表示不可拒绝原假设，说明假设寿命均值等于594是合理的。 95%置信区间为553.4962，634.5038完全包括594，估计精度较高。sig = 1远超过0.05，不可拒绝原假设所以可以认为刀具平均寿命为594（件）,例16 用模拟试验的方法直观地验证教材§6.2抽样分布定理一的结论。 假定变量 ，用随机数生成的方法模拟对 的500次简单随机抽样，每个样本的容量为16。利用这500×16个样本数据直观地验证样本均值 的抽样分布为均值等于60、方差等于2516的正态分布，即,解 >> %1、用随机数生成的方法模拟简单随机抽样x=;%生成一个存放样本数据的空表（维数可变的动态矩阵）for byk=1:500 %循环控制，循环执行下面的指令500次，本例中相当于500次抽样 xx=normrnd(60,5,16,1);%生成一个来自N(60,25)的容量为16 的样本（列向量） x=x,xx; %将样本数据逐列存入数表x，可从matlab的变量浏览器（workspace）中观察这个数表end %2、计算每个样本的样本均值（1500）xmean=mean(x);%可从变量浏览器中观察这500个数据 %3、绘制500个样本均值数据的直方图k=ceil(1.87*(length(x)-1)(2/5);%确定分组数h=histfit(xmean,k);%绘制附正态参考曲线的数据直方图set(h(1),'FaceColor','c','EdgeColor','w')%修饰，设置直方图线条颜色与填充色 %4、用这500个样本均值数据验证样本均值的均值和方差M=mean(xmean) %求（1500）样本的样本均值的均值V=var(xmean)%求（1500）样本的样本均值的方差,M = 59.9879V = 1.4129M = 60.0117V = 1.3900M = 59.9749V = 1.5158M = 59.9929V = 1.5757M = 59.8809V = 1.6855,例17 观察：用binornd模拟5000次投球过程，观察小球堆积的情况。,>> clear;clf, n=5;p=0.5;m=5000;x=0:1:nrand('seed',3)R=binornd(n,p,1,m);%模拟服从二项分布的随机数，相当于模拟 投球m次for I=1:n+1 %开始计数 k=; k=find(R=(I-1);%find是一个有用的指令，本语句的作用是找出R中等于(I-1)元素下标，并赋予向量k中 h(I)=length(k)/m;%计算落于编号(I-1)的格子中的小球频率endbar(x,h),axis(-1 6 0 1)%画频率图title('fontsize18fontname华文新魏5000次投球小球堆积的频率图'),>> f=binopdf(x,n,p), bar(x,f), axis(-1 6 0 1) title('fontsize18fontname华文新魏B(5,0.5)理论分布图'),例18 利用随机数样本验证中心极限定理。,独立同分布的随机变量的和的极限分布服从正态分布，通过产生容量为n的poiss分布和exp分布的样本，研究其和的渐近分布。 算法如下： 产生容量为n的独立同分布的随机数样本，得其均值和标准差； 将随机数样本和标准化； 重复、； 验让所得标准化的随机数样本和是否服从标准正态分布。 具体程序如下：,>> clearn=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);endnormplot(y);%分布的正态性检验title('poiss分布，中心极限定理'),