概率统计的数值实验-MATLAB在概率统计教学中的应用

引言

而MATLAB软件具有简单易学、易操作和绘图功能强等特点，利用MATLAB软件的图形可视功能将概率统计的内容用图形表示出来，通过图形让学生加深理解，以达到事半功倍的效果。

概率论与数理统计知识比较抽象，逻辑性较强。因此，建议让学生结合理论和公式推导，进行数值试验和相关调查，直观地感受数学概念和理论，从而提高学生解决实际问题的信心和能力。概率论

1.rand(m,n)：生成m×n的随机矩阵，每个元素都在(0,1)

间，生成方式为均匀分布。

2.randn(m,n)：生成m×n的随机矩阵，每个元素都在(0,1)

间，生成方式为正态分布。

3.randperm(m)：生成一个1~m的随机整数排列。

4.perms(1:n)：生成一个1~n的全排列，共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.91340.0975

0.63240.2785

>>randperm(5)%生成一个1~5的随机整数排列

ans=41523

>>a=[1242332];

unique(a)

ans=1234

例1

随机投掷均匀硬币，观察国徽朝上与国徽

朝下的频率。解>>n=3000~100000000;m=0;fori=1:nt=randperm(2);%生成一个1~2的随机整数排列x=t-1;%生成一个0~1的随机整数排列y=x(1);ify==0;m=m+1;endendp1=m/np2=1-p1

试验次数n300050001万2万3万国徽朝上频率0.50400.50060.48790.49990.5046国徽朝下频率0.49600.49940.51210.50010.4954试验次数n5万10万100万100万1亿国徽朝上频率0.50210.49990.49990.50010.5000国徽朝下频率0.49790.50010.50010.49990.5000可见当时，解记事件为第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;forj=1:mk=0;sui=randperm(n);fori=1:nifsui(i)==ik=k+1;elsek=k;endendifk==0p=p+1;elsep=p;endende=m/pe=2.7313模拟次数m400004000040000人数n100020005000e2.71552.70822.7202模拟次数m400040000400000人数n505050e2.73792.73132.7194

设针与平行线的夹角为，针的中心与最近直线的距离为。针与平行线相交的充要条件是，则所求概率为

故可得的近似计算公式，其中n为随机试验次数，m为针与平行线相交的次数。例3Buffon投针实验

在画有许多间距为的等距平行线的白纸上，随机投掷一根长为的均匀直针，求针与平行线相交的概率，并计算的近似值。解

>>clear,clfn=10000000;l=0.5;m=0;d=1;fori=1:nx=l/2*sin(rand(1)*pi);y=rand(1)*d/2;ifx>=ym=m+1;endendp1=m/npai=2*n*l/(m*d)试验次数n5千1万10万100万1000万针长l/平行间距d3/103/103/103/103/10相交频率0.18360.19710.18870.19050.1912π的近似值3.26803.04413.17983.14983.1387试验次数n5千1万10万100万1000万针长l/平行间距d2/52/52/52/52/5相交频率0.24960.25620.25490.25440.2543π的近似值3.20513.12263.13863.14513.1433试验次数n5千1万10万100万1000万针长l/平行间距d1/21/21/21/21/2相交频率0.32540.31480.31580.31780.3183π的近似值3.07313.17663.16673.14703.1417试验次数n5千1万10万100万1000万针长l/平行间距d4/54/54/54/54/5相交频率0.51420.51340.50860.50930.5093π的近似值3.11163.11653.14603.14183.1418试验次数n5千1万10万100万1000万针长l/平行间距d17/2017/2017/2017/2017/20相交频率0.54320.54520.54200.54120.5410π的近似值3.12963.11813.13663.14133.1426试验次数n5千1万10万100万1000万针长l/平行间距d9/109/109/109/109/10相交频率0.58600.57000.57560.57330.5731π的近似值3.07173.15793.12723.13953.1410例4在100个人的团体中，不考虑年龄差异，研究是否有两个以上的人生日相同。假设每人的生日在一年365天中的任意一天是等可能的，那么随机找n个人(不超过365人)。

(1)求这n个人生日各不相同的概率是多少？从而求这n个人中至少有两个人生日相同这一随机事件发生的概率是多少？

(2)近似计算在30名学生的一个班中至少有两个人生日相同的概率是多少？解:(1)>>clear,clfforn=1:100p0(n)=prod(365:-1:365-n+1)/365^n;p1(n)=1-p0(n);endp1=ones(1,100)-p0;n=1:100;plot(n,p0,n,p1,'--')xlabel('人数'),ylabel('概率')legend('生日各不相同的概率','至少两人生日相同的概率')axis([0100-0.11.199]),gridonp1(30)=0.7063,p1(60)=0.9941

分析：在30名学生中至少两人生日相同的概率为70.63％。下面进行计算机仿真。

随机产生30个正整数，代表一个班30名学生的生日，然后观察是否有两人以上生日相同。当30个人中有两人生日相同时，输出“1”，否则输出“0”。如此重复观察100次，计算出这一事件发生的频率。

(2)>>clear,clfn=0;form=1:100%做100次随机试验

y=0;x=1+fix(365*rand(1,30));%产生30个随机数

fori=1:29%用二重循环寻找30个随机数中是否有相同数

forj=i+1:30ifx(i)==x(j)y=1;break;endendendn=n+y;%累计有两人生日相同的试验次数endf=n/m%计算频率f=0.6900f=0.7900f=0.6700f=0.7300f=0.7500f=0.6900f=0.7200f=0.6700f=0.6800……重复观察，数据如下：例5Galton钉板模型和二项分布

Galton钉板试验是由英国生物统计学家和人类学家Galton设计的。故而得名。

通过模拟Calton钉板试验，观察和体会二项分布概率分布列的意义、形象地理解DeMoivre-Laplace中心极限定理。共15层小钉Ox-8-7-6-5-4-3-2-1

12345678高尔顿钉板试验小球最后落入的格数?记小球向右落下的次数为

则记小球向左落下的次数为

则符号函数,大于0返回1,小于0返回-1,等于0返回0

高尔顿(FrancisGalton,1822-1911)英国人类学家和气象学家Ox-8-7

-6

-5

-4

-3

-2

-1

12345678记则近似高尔顿钉板试验？什么曲线共15层小钉小球碰第层钉后向右落下小球碰第层钉后向左落下

模拟Galton钉板试验的步骤：

(1)确定钉子的位置：将钉子的横、纵坐标存储在两个矩阵X和Y中。

(2)在Galton钉板试验中，小球每碰到钉子下落时都具有两种可能性，设向右的概率为p，向左的概率为q＝1-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.51290.46050.35040.09500.43370.7092而且再次运行该指令时结果保持不变。除非重设种子seed的值，如>>rand('seed',2),u=rand(1,6)u=0.02580.92100.70080.19010.86730.4185这样结果才会产生变化。

将[0,1]区间分成两段，区间[0,p)和[p,1]。如果随机数u属于[0,p)，让小球向右落下；若u属于[p,1]，让小球向左落下。将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落人某一格子的过程。

(3)模拟小球堆积的形状。输入扔球次数m(例如m＝50、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;fori=n+1:-1:1%创建钉子的坐标x，yx(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;forj=2:ix(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);endendmm=moviein(m);%动画开始，模拟小球下落路径fori=1:ms=rand(1,n);%产生n个随机数

xi=x(1,1);yi=y(1,1);k=1;l=1;%小球遇到第一个钉子

forj=1:nplot(x(1:n,:),y(1:n,:),‘o’,x(n+1,:),y(n+1,:),‘.-’),%画钉子的位置axis([-2n+20y0+n+1]),holdon

k=k+1;%小球下落一格

ifs(j)>pl=l+0;%小球左移

elsel=l+1;%小球右移

endxt=x(k,l);yt=y(k,l);%小球下落点的坐标

h=plot([xi,xt],[yi,yt]);axis([-2n+20y0+n+1])%画小球运动轨迹

xi=xt;yi=yt;endballnum(l)=ballnum(l)+1;%计数

ballnum1=3*ballnum./m;bar([0:n],ballnum1),axis([-2n+20y0+n+1])%画各格子的频率

mm(i)=getframe;%存储动画数据

holdoffendmovie(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的时间间隔内收到的紧急呼救次数服从参数为t／2的泊松分布，而与时间间隔的起点无关(时间以小时计)，求：

(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;whilep<=0.9q=poisspdf([0:m],3);p=sum(q);m=m+1;endmm=6

计算结果表明当m=6时，p>0.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('\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')例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-.')holdonplot([-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:')holdofftext(-0.5,yy(6)+0.005,'\fontsize{14}95.44%')text(-0.5,yy(5)+0.005,'\fontsize{14}68.26%')text(-0.5,yy(7)+0.005,'\fontsize{14}99.74%')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σ')例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])holdoncapaplot(data,[xalpha4,inf]);axis([-3,3,0,0.45])holdoffxalpha1

xalpha2

xalpha3

xalpha4xalpha1=-1.6449xalpha2=1.6449xalpha3=-1.9600xalpha4=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次故障纪录，故障出现时该刀具完成的零件数如下：459362624542509584433748815505

612452434982640742565706593680926653164487734608428115359384452755251378147438882453886265977585975549697515628954771609402960885610292837473677358638699634555570844166061062484120447654564339280246687539790581621724531512577496468499544645764558378765666763217715310851

试观察该刀具出现故障时完成的零件数属于哪种分布？>>%数据输入x1=[459362624542509584433748815505];x2=[612452434982640742565706593680];x3=[9266531644877346084281153593844];x4=[527552513781474388824538862659];x5=[77585975549697515628954771609];x6=[402960885610292837473677358638];x7=[699634555570844166061062484120];x8=[447654564339280246687539790581];x9=[621724531512577496468499544645];x10=[764558378765666763217715310851];x=[x1x2x3x4x5x6x7x8x9x10];%作频数直方图hist(x,10)[N,X]=hist(x,10)%分布的正态性检验normplot(x)N=33714242214832X=1.0e+003*0.10420.21460.32500.43540.54580.65620.76660.87700.98741.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、方差等于25／16的正态分布，即解>>%1、用随机数生成的方法模拟简单随机抽样x=[];%生成一个存放样本数据的空表（维数可变的动态矩阵）forbyk=1:500%循环控制，循环执行下面的指令500次，本例中相当于500次抽样

xx=normrnd(60,5,16,1);%生成一个来自N(60,25)的容量为16的样本（列向量）

x=[x,xx];%将样本数据逐列存入数表x，可从matlab的变量浏览器（workspace）中观察这个数表end

%2、计算每个样本的样本均值（1~500）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)%求（1~500）样本的样本均值的均值V=var(xmean)%求（1~500）样本的样本均值的方差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:n]rand('seed',3)R=binornd(n,p,1,m);%模拟服从二项分布的随机数，相当于模拟投球m次forI=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([-1601])%画频率图title('\fontsize{18}\fontname{华文新魏}5000次投球小球堆积的频率图')>>f=binopdf(x,n,p),bar(x,f),axis([-1601])title('\fontsize{18}\fontname{华文新魏}B(5,0.5)理论分布图')例18

利用随机数样本验证中心极限定理。

独立同分布的随机变量的和的极限分布服从正态分布，通过产生容量为n的poiss分布和exp分布的样本，研究其和的渐近分布。

算法如下：①产生容量为n的独立同分布的随机数样本，得其均值和标准差；②将随机数样本和标准化；③重复①、②；④验让所得标准化的随机数样本和是否服从标准正态分布。

具体程序如下：>>clearn=2000;means=0;s=0;y=[];lamda=4;a=lamda;fori=1:nr=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分布，中心极限定理')例19在同一坐标轴上画box图，并对两个班的成绩进行初步的分析比较。

两个教学班各30名同学，在数学课程上，A班用新教学方法组织教学，B班用传统方法组织教学，现得期末考试成绩如下。

A：82，92，77，62，70，36，80，100，74，64，63，56，72，78，68，65，72．70，58，92，79,92，65，56，85，73，61，71，42，89

B：57，67，64，54，77，65，71，58，59，69，67，84，63，95，81，46，49,60,64，66，74，55，58，63，65，68，76，72，48，72解>>clear

x=[82,92,77,62,70,36,80,100,74,64,63,56,72,78,68,65,72,70,58,92,79,92,65,56,85,73,61,71,42,89;57,67,64,54,77,65,71,58,59,69,67,84,63,95,81,46,49,60,64,66,74,55,58,63,65,68,76,72,48,72];

boxplot(x')

从图中直观地看出，两个班成绩的分布是正态（对称）的，A班成绩较为分散(方差大)，B班成绩则较集中(方差小)。A班成绩明显高于B班(均值比较．并且A班25％低分段上限接近B班中值线，A班中值线接近B班25％高分段下限)。A班的平均成绩约为70分(中值)，B班约为65分(中值)。A班有一名同学的成绩过低(离群)，而B班成绩优秀的只有一人(离群)。

需要注意的是，从图中我们不能得出新教学方法一定优于传统教学方法的结论，因为我们并不知道两个班级原有的数学基础是怎样的。三、MATLAB也为常用的三大统计分布提供了相应的pdf、cdf、inv、stat、rnd类函数，具体分布类型函数名称如下：

分布类型MATLAB名称

分布chi2t分布tF分布f非中心分布ncx2非中心t分布

nct非中心F分布ncf例20分布的密度函数曲线。解：

>>%绘制不同自由度的卡方分布概率密度曲线clear,clfX=linspace(0,20,100);Y1=chi2pdf(X,1);%自由度等于1Y2=chi2pdf(X,3);%自由度等于3Y3=chi2pdf(X,6);%自由度等于6plot(X,Y1,'-g',X,Y2,'-b',X,Y3,'-k')title('\fontsize{18}\fontname{华文新魏}不同自由度的{\chi}^2分布概率密度曲线的比较')text(0.6,0.6,'\fontsize{12}df:n=1')text(2.6