MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第10章-Matlab在数理统计应用

【例10-1】设总体密度函数为，试从该总体中抽取容量为1000的简单随机样本。利用MATLAB编辑窗口保存以下程序，保存为ex11.m。n=1000;x=zeros(1,n);k=0;whilek<na=rand*pi-pi/2;b=rand/2;ifb<(cos(a)/2)k=k+1;x(k)=a;endend注意理解其原理。保存完成之后，在命令窗口执行ex11，则x被赋值。再执行下列命令，就可以得到这个容量为1000的样本的直方图，如图10-1所示。hist(x,-pi/2:0.2:pi/2)这里我们又用到了“hist”函数，其作用是在直角坐标系下绘制统计直方图，如图10-2所示，其格式为：hist(x,n)说明：x为统计数据，n表示直方图的区间数，缺省值n=10。同时，我们再介绍下“rose”函数，其作用是在极坐标系下绘制角度直方图，如图10-3所示，其用法为：rose(x,n)说明：n是在[0,2]范围内所分区域数，默认值n=20；x为指定的弧度数据。图10-1容量为1000的样本直方图图10-2直角坐标系下绘制的统计直方图图10-3极坐标系下绘制角度直方图【例10-2】某食品厂为加强质量管理，对生产的罐头重量X进行测试，在某天生产的罐头中抽取了100个，其重量测试数据记录如下：342340348346343342346341344348346346340344342344345340344344343344342343345339350337345349336348344345332342342340350343347340344353340340356346345346340339342352342350348344350335340338345345349336342338343343341347341347344339347348343347346344345350341338343339343346342339343350341346341345344342试根据以上数据作出X的频率直方图。解：在MATLAB编辑器中建立M文件LX0832.m：>>X=[342340348346343342346341344348...346346340344342344345340344344...343344342343345339350337345349...336348344345332342342340350343...347340344353340340356346345346...340339342352342350348344350335...340338345345349336342338343343...341347341347344339347348343347...346344345350341338343339343346...342339343350341346341345344342];>>hist(X,13)结果为：图10-4X的频率直方图【例10-3】某机床出次品的概率为0.01，求生产100件产品中：（1）恰有1件次品的概率；（2）至少有1件次品的概率。解：此问题可看作是100次独立重复试验，每次试验出次品的概率为0.01。（1）恰有1件次品的概率在MATLAB命令窗口键入：>>p=pdf('bino',1,100,0.01)p=0.3697或在MATLAB命令窗口键入：>>p=binopdf(1,100,0.01)p=0.3697（2）至少有1件次品的概率在MATLAB命令窗口键入：>>p=1-cdf('bino',0,100,0.01)p=0.6340或在MATLAB命令窗口键入：>>p=1-binocdf(0,100,0.01)p=0.6340【例10-4】自1875年到1955年中的某63年间，某城市夏季（5—9月间）共发生暴雨180次，试求在一个夏季中发生k次（k=0,1,2,…,8）暴雨的概率Pk（设每次暴雨以1天计算）。解：一年夏天共有天数为n=31+30+31+31+30=153故可知夏天每天发生暴雨的概率约为，很小，n=153较大，可用Poisson分布近似=np=。在MATLAB编辑器中编写M文件：LX0802.mp=input('inputp=')n=input('inputn=')lambda=n*pfork=1:9%循环变量的最小取值是从k=1开始。p_k(k)=poisspdf(k-1,lambda);endp_k在MATLAB的命令窗口键入LX0802，回车后按提示输入p和n的值，显示如下：inputp=180/(63*153)p=0.0187inputn=153n=153lamda=2.8571p_k=0.05740.16410.23440.22330.15950.09110.04340.01770.0063注意：在MATLAB中，p_k(0)被认为非法，因此应避免。【例10-5】某市公安局在长度为t的时间间隔内收到的呼叫次数服从参数为t/2的Poisson分布，且与时间间隔的起点无关（时间以小时计）。求：（1）在某一天中午12时至下午3时没有收到呼叫的概率；（2）某一天中午12时至下午5时至少收到1次呼叫的概率。解：在此题中，Lamda=t/2设呼叫次数X为随机变量，则该问题转化为：（1）求P{X=0}；（2）求1-P{X≤0}。解法一：在MATLAB命令窗口键入：>>poisscdf(0,1.5)%X=0表示0次呼叫，Lambda=t/2=1.5ans=0.2231即问题（1）中所求的概率为0.2231。>>1-poisscdf(0,2.5)ans=0.9179即问题（2）中所求的概率为0.9179。解法二：由于呼叫次数X≤0就是呼叫0次，即X=0。因此，此题也可用poisspdf求解。即：poisspdf(0,1.5)和1-poisspdf(0,2.5)。【例10-6】设盒中有5件同样的产品，其中3件正品，2件次品，从中任取3件，求不能取得次品的概率。在MATLAB编辑器中编辑M文件：LX0802.mN=input('inputN=')M=input('inputM=')n=input('inputn=')fork=1:M+1p_k=hygepdf(k-1,N,M,n)end在MATLAB的命令窗口键入LX0804，回车后按提示输入N、M、n的值，显示如下：inputN=5N=5inputM=2M=2inputn=3n=3p_k=0.1000p_k=0.6000p_k=0.3000这里，p_k=(0.10000.60000.3000)表示取到次品数分别为X=0,1,2的概率。【例10-7】计算正态分布N(0,1)在点0.7733的概率密度值。解：在MATLAB命令窗口键入：>>pdf('norm',0.7733,0,1)%利用通用函数ans=0.2958>>normpdf(0.7733,0,1)%利用专用函数ans=0.2958两者计算结果完全相同。【例10-8】绘制卡方分布密度函数在n分别等于1，5，15时的图形。解：在MATLAB编辑器中编辑M文件：LX0806.mx=0:0.1:30;y1=chi2pdf(x,1);plot(x,y1,':')holdony2=chi2pdf(x,5);plot(x,y2,'+')y3=chi2pdf(x,15);plot(x,y3,'o')axis([0,30,0,0.2])在命令窗口键入LX0806，回车后得结果为图10-8。图10-8卡方分布密度函数图【例10-9】某公共汽车站从上午7:00起每15分钟来一班车。若某乘客在7:00到7:30间的任何时刻到达此站是等可能的，试求他侯车的时间不到5分钟的概率。解：设乘客7点过X分钟到达此站，则X在[0，30]内服从均匀分布，当且仅当他在时间间隔（7:10，7:15）或（7:25，7:30）内到达车站时，侯车时间不到5分钟。故其概率为：p=P{10<X<15}+P{25<X<30}在MATLAB编辑器中建立M文件LX0807.m如下：formatratp1=unifcdf(15,0,30)-unifcdf(10,0,30);p2=unifcdf(30,0,30)-unifcdf(25,0,30);p=p1+p2运行结果为：p=1/3【例10-10】设X~N(3,22)，求P{2<X<5}，P{-4<X<10}，P{|X|>2}，P{X>3}。解：在MATLAB编辑器中编辑M文件LX0808.m如下：p1=normcdf(5,3,2)-normcdf(2,3,2)p2=normcdf(10,3,2)-normcdf(-4,3,2)p3=1-normcdf(2,3,2)+normcdf(-2,3,2)p4=1-normcdf(3,3,2)运行结果为：p1=1080/2027p2=2148/2149p3=7073/10138p4=1/2【例10-11】设随机变量X的概率密度为（1）确定常数c；（2）求X落在区间内的概率；（3）求X的分布函数F(x)。解：（1）在MATLAB编辑器中建立M文件LX08091.m如下：symscxp_x=c/sqrt(1-x^2);F_x=int(p_x,x,-1,1)运行结果为：F_x=pi*c由pi*c=1得c=1/pi（2）在MATLAB编辑器中建立M文件LX08092.m如下：symsxc=str2sym('1/pi');p_x=c/sqrt(1-x^2);formatratp1=int(p_x,x,-1/2,1/2)运行结果为：p1=1/3（3）在MATLAB编辑器中建立M文件LX08093.m如下：symsxtc=str2sym('1/pi');p_t=c/sqrt(1-t^2);F_x=int(p_t,t,-1,x)运行结果:F_x=asin(x)/pi+1/2>>simplify(F_x)ans=asin(x)/pi+1/2所以X的分布函数为：【例10-12】设lnX~N(1,22)，求P{}。解：利用对数分布累积专用函数，在MATLAB命令窗口键入：>>p=logncdf(2,1,2)-logncdf(1/2,1,2)p=232/965【例10-13】设X~N(3,22)，确定c，使得P{X>c}=P{X<c}。解：若要P{X>c}=P{X<c}，只需P{X>c}=P{X<c}=0.5在MATLAB命令窗口键入：>>c=icdf('norm',0.5,3,2)运行结果为：c=3【例10-14】在假设检验中，求临界值问题：已知：，查自由度为10的双边界检验t分布临界值。解：在MATLAB命令窗口键入：>>t0=icdf('t',0.025,10)运行结果为：t0=-586/263【例10-15】公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X（单位：cm）~N(175,36)，求车门的最低高度。解：设h为车门高度，X为男子身高，求满足条件P{X>h}≤0.01的h，即P{X<h}≥0.99。在MATLAB命令窗口键入：>>h=norminv(0.99,175,6)h=31556/167【例10-16】设二维随机向量(X,Y)的联合密度为求：（1）P{0<X<1,0<Y<1}；（2）（X，Y）落在x+y=1,x=0,y=0所围成区域G内的概率。解：在MATLAB编辑器中编辑M文件LX0814.m如下：symsxyf=exp(-x-y);p_XY=int(int(f,y,0,1),x,0,1)p_G=int(int(f,y,0,1-x),x,0,1)运行结果为：p_XY=exp(-2)*(exp(1)-1)^2p_G=1-2*exp(-1)【例10-17】设随机变量X的分布律为：X-2-1012p0.30.10.20.10.3求EX，E(X2-1)。解：在MATLAB编辑器中建立M文件LX0815.m：X=[-2-1012];p=[0.30.10.20.10.3];EX=sum(X*p')Y=X.^2-1;EY=sum(Y*p')运行结果：EX=0EY=8/5LX0815.m也可以写成如下形式：X=[-2-1012];p=[0.30.10.20.10.3];EX=sum(X.*p)Y=X.^2-1;EY=sum(Y.*p)运行结果：EX=0EY=8/5【例10-18】随机抽取6个滚珠测得直径如下：（单位：mm）14.7015.2114.9014.9115.3215.32试求样本均值。解：在MATLAB命令窗口键入：>>X=[14.7015.2114.9014.9115.3215.32];>>mean(X)ans=15.0600【例10-19】已知随机变量X的概率密度求EX和E(4X-1)。解：在MATLAB编辑器中建立M文件LX0817.m：symsxp_x=3*x^2;EX=int(x*p_x,0,1)EY=int((4*x-1)*p_x,0,1)运行结果为：EX=3/4EY=2【例10-20】设随机变量X的概率密度为求EX。解：在MATLAB编辑器中建立M文件LX0818.m：symsxp_x=1/2*exp(-abs(x));EX=int(x*p_x,-inf,inf)运行结果为：EX=0【例10-21】设随机变量X的分布律为X-2-1012p0.30.10.20.10.3求DX，D(X2-1)。解：在MATLAB编辑器中建立M文件LX0819.m：X=[-2-1012];p=[0.30.10.20.10.3];EX=sum(X.*p)Y=X.^2-1EY=sum(Y.*p)DX=sum(X.^2.*p)-EX.^2DY=sum(Y.^2.*p)-EY.^2运行结果：EX=0Y=30-103EY=8/5DX=13/5DY=76/25【例10-22】求下列样本的样本方差和样本标准差，方差和标准差14.7015.2114.9014.9115.3215.32解：在MATLAB编辑器中建立M文件LX0820.m：X=[14.7015.2114.9014.9015.3215.32];DX=var(X,1)%方差sigma=std(X,1)%标准差DX1=var(X)%样本方差sigma1=std(X)%样本标准差运行结果为：DX=241/4272sigma=371/1562DX1=371/1562sigma1=364/1399【例10-23】设X的密度函数为求DX，D(2X+1)。解：在MATLAB编辑器中建立M文件LX0821.m：symsxpx=1/(pi*sqrt(1-x^2));EX=int(x*px,-1,1)DX=int(x^2*px,-1,1)-EX^2y=2*x+1;EY=int(y*px,-1,1)DY=int(y^2*px,-1,1)-EY^2运行结果为：EX=0DX=1/2EY=1DY=2 【例10-24】求参数为0.12和0.34的分布的期望和方差。解：在MATLAB命令窗口键入：>>[m,v]=betastat(0.12,0.34)结果为：m=6/23v=271/2052【例10-25】按规定，某型号的电子元件的使用寿命超过1500小时为一级品，已知一样品20只，一级品率为0.2。问这样品中一级品元件的期望和方差为多少？解：分析可知此电子元件中一级品元件分布为二项分布，可使用binostat函数求解。在MATLAB命令窗口键入：>>[m,v]=binostat(20,.2)结果为：m=4v=16/5结果说明一级品元件的期望为4，方差为3.2。【例10-26】求参数为8的Poisson分布的期望和方差。解：>>[m,v]=poisstat(8)m=8v=8由此可见Poisson分布参数的值与它的期望和方差是相同的。【例10-27】设(X,Y)的联合分布为YX-112-12Z=X-Y，求EZ。解：在MATLAB编辑器中建立M文件LX0825.m：X=[-12];Y=[-112];fori=1:2forj=1:3Z(i,j)=X(i)-Y(j);endend%该循环计算X-Y的值Zp=[5/202/206/20;3/203/201/20];EZ=sum(sum(Z.*p))%将Z与p对应相乘相加运行结果为：EZ=-1/2【例10-28】射击试验中，在靶平面建立以靶心为原点的直角坐标系，设X、Y分别为弹着点的横坐标和纵坐标，它们相互独立且均服从N(0,1)，求弹着点到靶心距离的均值。解：弹着点到靶心的距离为，求EZ。其联合分布密度为在MATLAB编辑器中建立M文件LX0826.m：symsxyrtpxy=1/(2*pi)*exp(-1/2*(x.^2+y.^2));EZ=int(int(r*1/(2*pi)*exp(-1/2*r^2)*r,r,0,inf),t,0,2*pi)%利用极坐标计算较简单运行结果为：EZ=(2^(1/2)*pi^(1/2))/2即EZ=【例10-29】设(X,Y)的联合密度为求DX，DY和cov(X,Y)。解：在MATLAB编辑器中建立M文件LX0827.m：symsxypxy=1/8*(x+y);EX=int(int(x*pxy,y,0,2),0,2)EY=int(int(y*pxy,x,0,2),0,2)EXX=int(int(x^2*pxy,y,0,2),0,2)EYY=int(int(y^2*pxy,x,0,2),0,2)EXY=int(int(x*y*pxy,x,0,2),0,2)DX=EXX-EX^2DY=EYY-EY^2DXY=EXY-EX*EY运行结果为：EX=7/6EY=7/6EXX=5/3EYY=5/3EXY=4/3DX=11/36DY=11/36DXY=-1/36【例10-30】求向量a=[121221]的协方差。解：在MATLAB命令窗口键入：>>a=[121221];>>cov(a)ans=3/10【例10-31】求矩阵的协方差。解：在MATLAB命令窗口键入：>>d=rand(2,6)d=527/3095837/1562974/3643636/14751137/1448411/8003585/8811145/1381362/2055127/267393/30082660/4239>>cov1=cov(d)cov1=348/2855299/4130-97/4308617/56203-685/4239817/5744299/4130735/17093-191/1428263/9662-253/2636178/2107-97/4308-191/1428247/11300-110/542431375/46063-50/1903617/5620363/9662-110/5424368/68775-100/6871129/10070-685/4239-253/26361375/46063-100/6871277/1293-531/2816817/5744178/2107-50/1903129/10070-531/28162599/15659【例10-32】设(X,Y)的联合分布律为YX-112-12求X与Y的协方差及相关系数。解：在MATLAB编辑器中建立M文件LX0830.m：formatrat%有理格式输出X=[-12];Y=[-112];PXY=[5/202/206/20;3/203/201/20];%X、Y的联合分布PX=sum(PXY')%X的边缘分布PY=sum(PXY)%Y的边缘分布EX=sum(X.*PX)%X的期望EY=sum(Y.*PY)EXX=sum(X.^2.*PX)%计算EX2EYY=sum(Y.^2.*PY)DX=EXX-EX^2%X的方差DY=EYY-EY^2XY=[1-1-2;-224];%XY的取值EXY=sum(sum(XY.*PXY))DXY=EXY-EX*EY%X与Y的协方差ro_XY=DXY/sqrt(DX*DY)%X与Y的相关系数运行结果为：PX=13/207/20PY=2/51/47/20EX=1/20EY=11/20EXX=41/20EYY=41/20DX=819/400DY=699/400EXY=-1/4DXY=-111/400ro_XY=-365/2488即：=-111/400，=-365/2488【例10-33】设(X,Y)在单位圆G={(x,y)|x2+y2≤1}上服从均匀分布，即有联合密度求，，，。解：在MATLAB编辑器中建立M文件LX0831.m：symsxyrtpxy=str2sym('1/pi');

EX=int(int(r^2*cos(t)*pxy,r,0,1),0,2*pi)EY=int(int(r^2*sin(t)*pxy,r,0,1),0,2*pi)EXX=int(int(r^3*cos(t)^2*pxy,r,0,1),0,2*pi)EYY=int(int(r^3*sin(t)^2*pxy,r,0,1),0,2*pi)EXY=int(int(r^3*cos(t)*sin(t)*pxy,r,0,1),0,2*pi)DX=EXX-EX^2DY=EYY-EY^2DXY=EXY-EX*EYro_XY=DXY/sqrt(DX*DY)运行结果为：EX=0EY=0EXX=1/4EYY=1/4EXY=0DX=1/4DY=1/4DXY=0ro_XY=0【例10-34】设总体，为来自总体的简单随机样本，欲估计总体均值（注意未知），比较以下三个点估计量的好坏：，，解：本例题给出了利用MSE评价点估计量的随机模拟方法。由于的总体均值为，因此我们可以先取定一个固定值，例如，然后在这个参数已知且固定的总体中抽取容量为20的样本，分别用样本值依照三种方法分别计算估计值（注意谁也别偷看底牌），看看哪种方法误差大，哪种方法误差小。一次估计的比较一般不能说明问题，正如低手射击也可能命中10环，高手射击也可能命中9环。如果连续射击1万次，比较总环数（或平均环数），多者一定是高手。同理，如果抽取容量为20的样本次，分别计算小者为好。N=10000;m=5;n=20;mse1=0;mse2=0;mse3=0;fork=1:Nx=chi2rnd(m,1,n);m1=101*x(1)-100*x(2);m2=median(x);m3=mean(x);mes1=mse1+(m1-m)^2;mes2=mse2+(m2-m)^2;mes3=mse3+(m3-m)^2;endmse1=mes1/Nmse2=mes2/Nmse3=mes3/N以上程序保存为ex21.m，命令窗口中键入ex21，运算结果为mse1=2197/355mse2=1/2242736mse3=5/166436可见第一个虽为无偏估计量，但MSE极大，表现很差。第二个虽为有偏估计，但表现与第三个相差不多，也是较好的估计量。另外，重复运行ex21，每次的结果是不同的，但优劣表现几乎是一致的。【例10-35】设为来自上服从均匀分布的总体的简单随机样本，容易得到未知参数的矩估计量，最大似然估计量，试用随机模拟的方法比较两者的优劣。解：不妨设，以下程序给出了两者的评价。s=5;N=10000;mse1=0;mse2=0;fork=1:Nx=5.*rand(1,50);s1=2*mean(x);s2=max(x);mse1=mse1+(s1-s)^2;mse2=mse2+(s2-s)^2;endmse1=mse1/N;mse2=mse2/N;[mse1,mse2]参考运行结果：176/1063157/8123【例10-36】设1.1,2.2,3,3,4.4,5.5为来自正太总体的简单随机样本，求的置信水平为95%的置信区间。解：以下用MATLAB命令计算：>>x=[1.1,2.2,3.3,4.4,5.5];>>m=mean(x);>>c=2.3;>>d=c*norminv(0.975);>>a=m-d;>>b=m+d;>>[a,b]ans=-3265/27031382/177【例10-37】设两台车床加工同一零件，各加工8件，长度的误差为：A：-0.12-0.80-0.05-0.04-0.010.050.070.21B：-1.50-0.80-0.40-0.100.200.610.821.24求方差比的置信区间。解：用MATLAB计算如下：>>x=[-0.12,-0.80,-0.05,-0.04,-0.01,0.05,0.07,0.21];>>y=[-1.50,-0.80,-0.40,-0.10,0.20,0.61,0.82,1.24];>>v1=var(x);>>v2=var(y);>>c1=finv(0.025,7,7);>>c2=finv(0.975,7,7);>>a=(v1/v2)/c2;>>b=(v1/v2)/c1;>>[a,b]ans=121/5278910/1591【例10-38】随机产生100个分布数据，相应的分布参数真值为4和3。求4和3的最大似然估计值和置信度为99%的置信区间。解：在MATLAB编辑器中建立M文件LX0833.m：>>X=betarnd(4,3,100,1);%随机产生100个分布数据，参数为4和3。>>reshape(X,25,4)%这些数据只有一列，这里为了节约版面而改为4列数据[phat,pci]=betafit(X,0.01)运行结果：ans=909/1640319/41687/113153/269703/1298309/349869/3999647/13681173/20271133/2418512/707585/1457728/9211041/15611451/2319534/713449/1252559/8721031/1746661/2291117/5231030/1987674/863360/1067841/1330308/1129374/15173163/3982271/560925/1228626/1513211/287783/13871369/2065483/1256641/8101219/1837397/8471174/21833103/55451117/1953551/990563/7442353/40391079/1906897/994916/1173665/858768/8992097/2662863/2154569/1268541/896283/551769/18431120/4177220/8411953/23921135/2264488/541299/842755/1264496/7072119/2651292/1165337/643544/887511/773954/1199345/956498/9292843/6016281/3952549/3040353/1128595/6891670/2239682/1319474/899310/1111424/6237534/10389519/6321422/20952593/3141551/9071023/252449/115652/12651517/1920499/9571285/1712444/931393/11111247/1561478/1243797/110688/2631917/1996279/485>>[phat,pci]=betafit(X,0.01)phat=3634/9871378/531pci=1902/7852411/13748605/1538971/253说明：数据4.6613和3.5719为参数4和3的估计值；pci的第1列为参数4的置信区间，第2列为参数3的置信区间。随机产生的数据不同，其估计值和置信区间就不一样。【例10-39】设某种油漆的9个样品，其干燥时间（以小时计）分别为6.05.75.86.57.06.35.66.15.0设干燥时间总体服从正态分布N()，求和的置信度为0.95的置信区间（未知）。解：在MATLAB命令窗口键入：>>X=[6.05.75.86.57.06.35.66.15.0];>>[muhat,sigmahat,muci,sigmaci]=normfit(X,0.05)muhat=6%的最大似然估计值sigmahat=1030/1793%的最大似然估计值muci=%的置信区间9227/166010693/1660sigmaci=%的置信区间149/3842299/2089此解说明的最大似然估计值为6，置信区间为[9227/1660，10693/1660]；的最大似然估计值为1030/1793，置信区间为[149/384，2299/2089]。【例10-40】分别使用金球和铂球测定引力常数（1）用金球测定观察值为：6.6836.6816.6766.6786.6796.672（2）用铂球测定观察值为：6.6616.6616.6676.6676.664设测定值总体服从N()，和为未知。对（1）、（2）两种情况分别求和为的置信度为0.9的置信区间。解：在MATLAB命令窗口键入：>>j=[6.6836.6816.6766.6786.6796.672];>>b=[6.6616.6616.6676.6676.664];>>[muhat,sigmahat,muci,sigmaci]=normfit(j,0.1)%金球测定的估计muhat=581/87sigmahat=37/9564muci=10433/15632579/386sigmaci=93/3577079/9774说明金球测定数据的置信度为0.9的和置信区间为：[10433/1563，2579/386]：[93/35770，79/9774]>>[muhat,sigmahat,muci,sigmaci]=normfit(b,0.1)%铂球测定的估计muhat=833/125sigmahat=3/1000muci=2811/42211487/1723sigmaci=65/3336969/9695说明铂球测定数据的置信度为0.9的和置信区间为：[2811/422，11487/1723]：[65/33369，69/9695]【例10-41】已知小麦亩产服从正态分布，传统小麦品种平均亩产800斤，现有新品种产量未知，试种10块，每块一亩，产量为：775,816,834,836,858,863,873,877,885,901新产品亩产是否超过了800斤？假设检验就是概率意义上的反证法。要证明命题H1：，可以首先假设H0：。本体中容易计算样本均值超过800了，有没有可能超过800的原因是由于抽样的随机性引起的？是否总体均值根本没有变化？我们看如下的统计量：容易看出，如果新品种确有增产效应，应偏大，不利于H0，取，查表求临界值，使得，即构造不利于H0，有利于H1的小概率事件，如果在一次试验中该小概率事件发生了，就有理由拒绝H0，认为H1成立。严格逻辑意义上的反证法思路如下：欲证H1成立，先假设其否命题H0成立，然后找出逻辑意义上的矛盾，从而推翻H0成立，严格证明H1成立。假设检验的思路类似，只不过引出的不是矛盾，而是小概率事件在一次实验中发生。我们称想要证明的命题H1为备择假设，对立的命题H0称为原假设，面对样本，我们必须表态是接受原假设还是拒绝原假设，这有可能出现两类错误。如果客观上原假设的确成立，面对样本的异常我们拒绝了原假设，这种“以真为假”的错误我们称为第一类错误，发生的概率用表示；如果客观上备择假设成立，我们却接受了原假设，这种“以假为真”的错误我们称为第二类错误，用发生的概率用表示。假设假设检验一般首先控制第一类错误，即：当我们拒绝原假设时有比较充足的理由，犯错误的概率不超过预设的，称为显著性水平。常用的显著性水平有这种预设显著性水平的假设检验也称为显著性检验，以后我们提到的假设检验都是显著性检验。对于显著性检验，当接受原假设时，可以认为是拒绝的证据不足。对于例10-41的问题，取，当时拒绝原假设。这里称为检验统计量，所确定的的取值范围称为拒绝域。>>x=[775,816,834,836,858,863,873,877,885,901];>>T=(mean(x)-800)/(std(x)/sqrt(9))T=3146/755>>ta=tinv(0.95,9)ta=1384/755计算结果T=3146/755>ta=1384/755，故拒绝原假设，认为确有增产。【例10-42】某车间用一台包装机包装葡萄糖，包得的袋装糖重是一个随机变量，它服从正态分布。当机器正常时，其均值为0.5kg，标准差为0.015。某日开工后检验包装机是否正常，随机地抽取所包装的糖9袋，称得净重为：（kg）0.4970.5060.5180.5240.4980.5110.520.5150.512问机器工作是否正常？解：总体和已知，则可设样本的=0.015，于是X~N(,0.0152)，问题就化为根据样本值来判断=0.5还是≠0.5。为此，提出假设原假设：H0:0.5备择假设：H1:MATLAB实现>>X=[0.4970.5060.5180.5240.4980.5110.520.5150.512];>>[H,sig]=ztest(X,0.5,0.015,0.05,0)%注意：此处数据X只能为向量而非矩阵H=1sig=275/11087结果H=1，说明在0.05的水平下，可拒绝原假设，即认为这天包装机工作不正常。方差未知情形原假设H0：此时，检验统计量为，H0成立时，依据备择假设的不同提法，分三种情况分别给出拒绝域。双侧检验。备择假设H1：拒绝域：单侧检验（右侧）。备择假设H1：拒绝域：单侧检验（左侧）。备择假设H1：拒绝域：下面我们介绍一下在MATLAB中进行检验的函数。函数：ttest格式：H=ttest(X,m,alpha)[H,sig,ci]=ttest(X,m,alpha,tail)说明：X是样本m期望值alpha是检验水平（默认为0.05）tail——备选假设的选项，有三种情况：tail=0（默认）：tail=1：tail=-1：即：tail=0为双边检验，其余为单边检验问题。H是检验结果，两种情况：H=0是在水平下，接受原假设，或假设相容；H=1是在水平下，拒绝原假设，或假设不相容。sig是当原假设为真时（即成立），得到观察值的概率，当sig为小概率时，则对原假设提出质疑。ci是均值的置信度为的置信区间。【例10-43】某种电子元件的寿命X（以小时计）服从正态分布，和均未知。现测得16只元件的寿命如下：159280101212224379179264222362168250149260485170问是否有理由认为元件的平均寿命大于225小时？解：未知，按题意作如下假设H0：，H1：取。在MATLAB实现>>X=[159280101212224379179264222362168250149260485170];注意：此处数据X只能为向量而非矩阵>>[h,sig]=ttest(X,225,0.05,1)h=0sig=451/1755结果表明，h=0，即在显著水平为0.05的情况下，不能拒绝原假设，认为元件的平均寿命不大于225小时。【例10-44】在平炉上进行一项试验以确定改变操作的建议是否会增加钢的得率，试验是在同一只平炉上进行的。每炼一炉钢时除操作方法外，其它条件都尽可能做到相同。先用标准方法炼一炉，然后用建议的新方法炼一炉，以后交替进行，各炼10炉，其得率分别为：（1）标准方法：78.172.476.274.377.478.476.075.576.777.3（2）新方法：79.181.077.379.180.079.179.177.380.282.1设这两种方法相互独立，且分别来自正态总体N()和N()，、、均未知。问建议的新操作方法能否提高得率？（取）解：建立假设H0：=H1：<MATLAB实现如下：>>X=[78.172.476.274.377.478.476.075.576.777.3];>>Y=[79.181.077.379.180.079.179.177.380.282.1];>>[h,sig,ci]=ttest2(X,Y,0.05,-1)h=1sig=31/142468ci=-1/0-208/109结果h=1，表明在的显著水平下，可以拒绝原假设，即认为建议的新操作方法较原方法优。【例10-45】已知原来工艺下生产的某种灯泡的中位数为800小时，现改进生产工艺，试产10只灯泡，实验得到每只寿命为：775,816,834,836,858,863,873,877,885,901新工艺生产的灯泡寿命中位数是否超过了800小时？H0：一般情况下，灯泡寿命不是正态分布的。符号检验使用的是计数统计量，先设则有即记录样本点中大于800的个数。若H0成立，应该大约占样本容量的一半左右，若异常的大，说明备择假设H1：成立。H0成立时，，可以利用二项分布构造拒绝域：使得若H0成立时，，利用二项分布的分布律可以计算出临界值，用如下MATLAB函数文件计算：functiont=bt(n,a)SS=2^n*a;S=0;c=1;k=n+1;whileS<=SSk=k-1;S=S+c;c=c*k/(n-k+1);endt=k+1;以上自定义函数扩展了MATLAB的功能，可以替代教科书上的“符号检验临界值表”，并且可以使用任意的n及。注意：以上代码为自定义函数，不同于本书中之前新建文件操作(如例10-38的LX0833.m)，不能直接运行，需在命令窗口中调用。在例10-45中，，对于，使用命令t=bt(10,0.05)可以得到临界值9，临界值9，落在拒绝域内，故拒绝原假设，认为新工艺生产的灯泡寿命中位数超过了800小时。只要去代替，也可以进行双侧符号检验。【例10-46】20个品酒师对A、B两种白酒进行品尝，有17个品酒师认为A品质好，3个品酒师认为B品质好，在的显著性水平下，检验两种白酒品质是否存在差异？解：，设原假设为：H0：两种白酒品质无差异令表示认为A品质好的品酒师的人数，则H0成立时应该在10左右取值，如果值异常大，或者异常小，都说明两种白酒品质有差异。取临界值与，使得，，由于关于对称，故有，因此可用水平为的单侧检验求出临界值。命令>>t2=bt(20,0.05/2)t2=15得到，因此，此例中拒绝域为：，或者落在拒绝域内，可以认为两种白酒品质有显著差异。【例10-47】某班级共15名同学，某次英语水平考试，分数如下：男：53，55，59，65，71，77，81女：56，62，68，76，84，86，90，96在显著性水平下，能否认为女生英语水平高于男生？要求采用Wilcoxon秩和检验。解：注意这是一个单侧检验问题，使用MATLAB命令：>>x=[53,55,59,65,71,77,81];>>y=[56,62,68,76,84,86,90,96];>>rsum(x,y)ans=78>>c=wr(7,8,0.05)c=78上述计算中，注意到rsum(x,y)=78，而临界值为c=78,的值落在拒绝域内，故可拒绝原假设，认为女生成绩显著高于男生。【例10-48】某班级共15名同学，某次英语水平考试，分数如下：53，55，59，65，71，77，81，56，62，68，76，84，86，90，96在显著性水平下，能否认为平均成绩高于60分？要求分别用：(1)符号检验；(2)Wilcoxon符号秩检验。解：注意这是一个单侧检验问题：H0：H1：使用MATLAB命令：>>x=[53,55,59,65,71,77,81,56,62,68,76,84,86,90,96];(1)符号检验注意这里n=15,B=11,利用前面自定义的bt.m函数计算：>>t=bt(15,0.05)t=12得到临界值，B=11<，没有落入拒绝域，故接受H0，认为平均成绩没有明显高于60分。(2)Wilcoxon符号秩检验>>n=15;>>E=n*(n+1)/4E=60>>wp=rpsum(x,60)wp=106计算结果发现wp=106>E=60，满足单侧检验条件一，再计算>>p=signrank(x,60)p=0.0071结果得p=0.0071<2，故拒绝原假设，认为平均成绩明显高于60分。【例10-49】旱灾土地总面积问题。下表10-6是从1990年至2010年全国因干旱而受灾的土地总面积（单位：千公顷）数。（数据来源于全国统计局官网）表10-6各个年份全国受灾土地总面积年份受灾面积年份受灾面积年份受灾面积199018175199733516200417253199124914199814236200516028199232981199930156200620738199321097200040541200729386199430423200138472200812137199523455200222124200929259199620152200324852201013259解决以下问题：（1）计算所给样本的均值与标准差；（2）检验在显著水平为0.05的情况下，全国每年因干旱而受灾的土地总面积是否服从正态分布；（3）如果服从正态分布，用极大似然估计法对未知参数μ和σ作出估计；（4）若年受旱灾总面积大于35000千公顷即为重灾年，根据估计出的μ值和σ值，计算当年为重灾年的概率。分析问题：这是一个样本均值和标准差的计算以及正态性检验和计算的一系列问题。对于此类问题可以应用数学软件MATLAB进行处理，应用MATLAB可以很容易的计算出均值及标准差，此外，采用Jarque-Beran检验即可知道其是否服从正态分布，并估计出总体的均值μ和标准差σ。解决问题：下面计算样本的均值和标准差>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>[h,stats]=cdfplot(X)h=Line-属性:Color:[00.44700.7410]LineStyle:'-'LineWidth:0.5000Marker:'none'MarkerSize:6MarkerFaceColor:'none'XData:[1×44double]YData:[1×44double]ZData:[1×0double]显示所有属性stats=包含以下字段的struct:min:12137max:40541mean:2.4436e+04median:23455std:8.1234e+03图10-17样本经验分布函数图从输出结果可看出，样本的最小值为12137，最大值为40541，中值为23455，均值为24436，标准差为8123.4。验证下面检验其是否服从正态分布。MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>normplot(X);>>[h,P,Jbstat,CV]=jbtest(X,0.05)>>title('正态概率图')>>xlabel('数据');>>ylable('概率')运行程序后，输出如下：h=0P=0.4574Jbstat=0.9230CV=3.8801图10-18正态概率图由输出结果h=0且Jbstat<CV可得出结论，在置信度α=0.05下，受灾面积（原始数据）服从正态分布，且在正态概率图中，各点均落在直线两侧，也可说明这一结论是成立的。再用极大似然估计法对未知参数μ和σ作出估计：MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>phat=mle(X,'distribution','norm','alpha',0.05)运行程序后，输出如下phat=1.0e+04*2.44360.7928即受灾面积的μ估计值为24436，σ估计值为7928。最后根据估计出的μ值和σ值，计算出每年的受灾面积大于35000千公顷的概率：MATLAB程序代码如下>>clearall;p=normspec([35000inf],2443