代码资源精通matlab科学计算13-3rd

13章概率统计计算统计工具统计工具箱介中的统计工具箱是一套建立在数值计算环境下的统计分析工具，该函数可以通过命令行或用户的应用程序调用，多数为中的M文件表13-1给出了统计工具箱中函数的主要分类。13-1文件输入/随量的数随量的数字特量，但能描述随量在某些方面的重要特征。本节主要讲述用求解随量属于统计工具的描述统计学（DescriptiveStatistics）部分。设离散性随量的分布律为P{Xxk}pk,k1,XE(X)xkk而对于来自总体X的一个样本设其样本值为x(x1,x2 ,xn)则定义样本均值为n1nx 可以证明x依概率收敛于X的均值因此通常用x来近似代替X的均值中提供了函数mean()来计算这种较为简单的均值，其调用格式如下。1．mXmean(X)2．mXdim认值为1。】解：通过在中编写程序ex1301.m来实现%ex1301.m用mean（）计算样本均值>>m=15.6】X1X2X3解：用矩阵的均值形式计算，在中编写程序ex1302.m来实现%ex1302.m1.8];%3.0];%2.8];% %将三个样本作为X %计算X %将三个样本作为X >m1=m2 X12.98X23.38X3的均值为3.38。【例13-3】 样本均值计算实例3。设随量X、Y的联合分布律如表13-2所示，分别计算E(X)，E(Y)，E(X+3)以及E(X+Y)。13-2X、YYX123400001020001 %ex1303.m %E(Xsum >>mX mY mX3 mXY 中还给出了求其他形式均值的一些函数，其命令形式与mean相同。各命令具体含义详见表13-3。13-3 如需了解表13-3中函数的详细用法，可到 的Help中的相关部分检索查询方差是用来刻画随量X取值分散程度的一个量。其一般用下式表达D(X)E{[x在应用上还引入与随量X具有相同量纲的量Xk

)E{[XE(X)]k},k1n1

S (xixnn1n1 nS (xixn1n1nn(xix2SSS或k

S Sn1Sn12n(xix1knBk (xix),knvar()std()来计算样本的标准差，用矩函数moment()计算样本的各阶中心矩，具体用法如下。var()Vvar(X) XX为矩阵则返回矩阵各列向Vvar(X,1) X1/n。var(X,0)等同于var(X)，其采用无偏式计算方差，前置因子为1/(n-1)。V XwX，wX，wXVvar(X,flag,dim)X在特定维上的方差，dim指定维数，flag指定选择的计算式。即flag0，选择无偏式计算；flag1，选择有偏式计算。的方差计算函数，调用格式与s 函数返回向量（矩阵）X的标准差（1/n-sstd(X flag01/n-1；flag11/nsstd(X,flag,dim)X在特定维上的标准差，dim指定维数，flag指定选函数用法同13-4】x＝[-1,-1,1,2]的方差和标准差以及后的方差和标准差，其中向量为w[1,2,3,4]（方差和标准差值指的解：通过 中编写程序ex1304.m来实现，具体代码如下所示%ex1304.m计算xx=[-111 %输入向量w=[123 %权向量>>v1 v2 v3 v4 v5 s1 s2 s3 moment() mmoment(Xorder dimX13-5】3解：通过在中编写程序ex1305.m来实现%ex1305.m3X=randn([5 %5 %计算矩阵X3 %计算矩阵X32>>==随量X、Y的协方差和相关系数的定义式如下cov(x,y)E{[xE(x)][yE(cof(x,y)

cov(x,D(x)D(对于n维随量，通常用协方差矩阵描述它的2阶中心矩。如对于二维随

c12c22

c11E{[xc12E{[xE(x)][yE(c21E{[yE(y)][xE(x)]}c22E{[yE(y)]2}中用cov()和corrcoef()计算协方差和相关系数。具体用法如下。Ccov(X)X为向量时，函数返回此向量的方差。X为矩阵时，矩阵的每一行表示一元素是X矩阵的列向量的方差值。Ccov(X,Y)X、YX、Y 计算协方差矩阵时前置系数取1/n。cov(X,0)与cov(X)相同，都是取前置系数为1/n-1，此用法可参考var函数。corrcoef()Rcorrcoef(X)X的相关系数矩阵，其各点值对应于相关矩阵的各点值除以Rcorrcoef(x,y) x、yx、y分别为列向量，则该命令等同于Rcorrcoef([xy])。[R,P Pp[R,P,RLO,RUP]corrcoef(...) 对于每一个R值返回的95%置信区间为[RLO,RUP]。[...]corrcoef(...,'param1',val1,'param2',val2,...) corrcoef13-6xy的协方差和相关系数。向量x与向量y分别如下所示：xy解：通过 中编写程序ex1306.m来实现，具体代码如下所示%ex1306.m %计算向量x的协方差，可知其等于x >>ex1306 cx=vx=cxy cor 【例13-7】随量数字特征综合计算实例。根据表13-4中所示的数据计算各变量13-4解：通过 中编写程序ex1307.m来实现，具体代码如下所示%ex1307.mVOL=[899292928992948989HP=[495570537055807366MPG=[65.455.949.046.546.245.443.439.339.6SP=[96971059610597107103100 %计算X %计算X >>===cv=corX随量x的偏斜度随量x的峰度

E(xy[D(x)]3/E(xk在中分别用函数skewness()、kurtosis()来实现上述计算skewnessyskewness(X,flag,dim)，在矩阵X的dimkurtosiskkurtosis(X)Xkkurtosis(Xflag)，flagskewnesskkurtosis(Xflagdim)skewness通过上面各式可以计算出向量（矩阵）Xflag1表示不进行偏移修正；flag0，表示进行偏移修正。dim指定在矩阵X的特定维上求相应统计量。下面是一段 命令窗口中输入的计算偏斜度和峰度的例子 %输入向量y=[10.2,9.4,8.5;12.1,10.3,11.0;12.1,13.4,11.4]输入矩阵>>ans%求x>>ans>>ans>>%求xans=>>%求yans >>skewness(y,1,2) ans=>>kurtosis(y) %求y峰度，flag1，dim1ans= 除了上面介绍的几类重要的数字特征外，的统计工具箱中还有很多对数据13-613-5 求随量的中 75%25%13-6 计算jackknife通过对数据进行bootstrap计算bootstrap13.3特殊分布的概率计正态分布和泊松分布等。中也有专门的对这些特殊分布进行操作处理的函数。20多种特殊概率分布。对于每种分布，该工具箱都提供表13-7给出了支持的20多种分布以及它们名称的字母缩写。可表13-7支持的分布类BetaF2t2Fftt23Ypdf('name',X,A1,A2,其中，name13-7中分布的字母缩写，X为样本矩阵，A1、A2A3是分布参数矩阵，Y为概率密度矩阵。对于某些分布，有些参数矩阵可以不入。X、A1、A2和A3必须是具有同样的大，另外对于每一种分布还有一个的求其概率密度的函数它们的基本用法和含义见表13-8。，13-8Beta2FYmvnpdf(X);Ymvnpdf(X,mu)tFt2tX13-8】概率密度计算实例。计算坐标[0,2]0.4，0.6，1.4，1.6处解：可以分别用通用函数pdf和函数unifpdf求解。在中编写程ex1308.mclearall 用用>>f1 f2 解：通过 中编写程序ex1309.m来实现，具体代码如下所示%ex1309.mclear 13-1

0

13-1若随量X的概率密度函数为f(x)，则累积分布函数的定义如下F(x)P(X

x)

f与求pdf类似， 含义可参见pdf。表13-9给出了各分布的的累积分布函数。其具体参数含义和对应的pdf函数13-9Beta2经验分布（kan-Meier布[f,FFt2ticdf，具体调用格式与cdf相同。表13-10给出了中的的逆累积分布函数。13-10Beta2FFt2t【例13-10】 P(X

，P(1X

，P(X6)，P(|X|P(Xc0.5c解：通过在中编写程序ex1310.m来实现，具体代码如下所示%ex1310.mcleary1=normcdf(4,1,3%P(X<4)y2=normcdf(5,1,3)-normcdf(1,1,3)%P(1<X<5)y3=1-normcdf(6,1,3)%P(X>6)y4=normcdf(4,1,3)-normcdf(-4,1,3)%P(|X|<4)%%%%%xc=norminv(0.5,1,3)>y1=y2=y3=y4=c1由计算结果可知，PX

P(|X|40.7936。当c1P(Xc0.5【例13-11】概率计算应用实例。设在一次实验中某发生的概率为0.4，重复该实验100次求出现该次数不小于50次的概率若该刚好发生30次求其概率ex1311.m%ex1311.mclear >>p1=p2=由以上计算结果可知，出现该次数不小于50的概率为0.0271，该刚好发300.01除了上节介绍的简单的均值函数mean以及方差函数var外，中还为每一种特殊分布提供了相应的计算其均值和方差的函数，详细函数说明见表13-11。13-11Beta2FFt2t【例13-12】 特殊分布的期望和方差计算实例1。求参数为0.4和0.8的分布的期解：在命令窗口中输入>>mv=【例13-13】 特殊分布的期望和方差计算实例2。求参数为5的泊松分布的期望和解：在命令窗口中输入>>m5v5由此例可知泊松分布的均值和方差相同，且都等于参数【例13-14】 特殊分布的期望和方差计算实例3。分别计算参数为1～6的离散均匀解：在命令窗口中输入>>mv013-11中的函数还可以同时计算具有不同参数的相同类型分布的均值random和特定分布的随机数产生函数（以rnd结尾。可以直接调用这些函数来获得所需的随机数。13-12给出了几13-12Beta2F13-15】0434解：在命令窗口中输入 %34列的N（0，4）n2=normrnd(0,2,[3 %34列的N（0，4） ====13.4参数估的统计工具箱中对一些常用的特殊分布也提供了相应的极大似然估计函数（fit结尾13-1313-13phat为Betaa、b为alphaphat为二项分布中发生概率的极大似然估计为alphaphatevfit(x)phat为极值分布的均值和标准差极大似然估计；pcialphaphat的极大似然估计；pcialphaphat分布参数的极大似然估计；pcialpha[phat,pci]phat计；pci为alpha[phat,pci]phat为负二项分布参数的极大似然估计；pcimuhat、sigmahat分别为正态分布的均值和标准差的极大似然估计；muci，sigmaci为相应的alphalambdshatlambdaci为alphaphat为瑞利分布参数的极大似然估计；pcialphaahat,bhata、b估计；aci、bcialphaphat为分布参数的极大似然估计；pci为alpha水信区间，其中α默认为0.05。例如：[pha,pci]pci95%mle（1）phat（2）[phat,pci]95%（3）[phat,pci]mle(x,'distribution',dist95%置信区间。13-16】1X在[a,b]上服从均匀分布，a，ba，b95%解：下面给出分别用函数unifit以及通用函数mle求解的过程，在 >>[ahat,bhat,aci,bci]=unifit(x)%unifitahatbhat5acibci= 命令窗口中继续输入>> %mlephat pci 1.2000由此解可得a的估计值为1.2，置信区间为[-71,1.2]；b的估计值为5[5,77.2]13-17】2。一地质学家为研究密执根湖湖滩地区的岩石成分，随机地从该地区取100个样品，每个样品有10块石子，记录了每个样品中属石灰分布。pp的极大似然估计值。该地质学家所得数据如表13-14所列。表 90310解：在命令窗口中输入 %100 %每一个样本的p %总体的pp1=13-18】3。已知下列数据为指数分布，求它们的点估计值和97%的置信区间。数据为：1，6，8，25，28，15，1，3，8。解：在命令窗口输入>> >>parmhat=parmci13.5假设检H0H1H0，进而得到H0的域；根据统计量的值是否落入域，做出或保留H0的判断。13-15所示的几种假设检验函数。13-152HPcdfLilliefors0ptxmt2FXXzsigma表中X为输入样本；H为检验结果，其值为0支持原假设，其值为1原假设；显著水平α默认值为0.05，也可以在参数中进行设定。p值是原假设的最小的α值。N(,22的检验（Z检验Zztest来实现，其调用格式有如下几种。hztest(x,m,sigma)在α0.05sigmam。xxhztest(xmsigmaalpha)alpha[hsig,ci]ztest(xmsigmaalpha,tail)tail（13-16）对三种备择假设形式进行选择，sig为p值，ci为1-alpha置信区间。13-16单样本均值检验：tailtailm(xmm(xmm(xm13-19】1。某糖厂有一台自动打包机打包，额定标准每包质量为100kg。设每包质量服从正态分布，且根据以往经验，其方差2(0.4)2。9998.5102.51019899102102.10.05解：H0100；H1100z统计量检验假设。具体通过下面在命令窗口中输入>>x=[9998.5102.51019899102102.1>>>>h1在命令窗口中继续输入>>h1sig=ci 由h1知原假设，即这天打包机打包的质量与规定质量有显著差异。另外可得0.030395%置信区间为[100.276,100.5502]2的检验（t检验在中用ttest来实现。ttest具体调用格式如下（1）h（2）hm（3）[h,sig]ttest(x,m,alpha,alphamztest13-20】26块进行设砖的抗断强度服从正态分布，问这批砖的平均抗断强度是否不大于32.50(kg/cm2)解：H0：„H1：t统计量进行检验。具体在命令窗口中输入>>x=[32.5629.6632.6430.0031.87>>>>h0断强度不大于32.5kg/cm2。两个正态总体均值差的检验（t在中用ttest2来检验两个正态样本的均值是否相同。这两个正态样本方httest2(x,[h,sig,ci]ttest2(x,y,alpha,x、y为两个输入样本，其长度可以不同。alpha为显著性水平，tail13-1713-17两样本均值检验：tailtailxxx13-21】313-18分别给出两个文学家马克·8篇小品文以及特·格拉斯的10篇小品文中由3个字母组成的单词的比例。问两个作家所写的小品文中包含由3个字母组成的单词的比例是否有显著差异。13-183解：取α0.05

H0：1H1：1ttest2在命令窗口中输入>>h1基于成对数据的检验（tttest的另一用法便是基于成对数据的检验，其具体调用格式如下：httest(x,y),x、yx-y构成的正态分布的均13-22】4Ix，Iy，用来测量材料中某种金属9（13-1913-19问：能否认为这两台仪器的测量结果有显著的差异（取解：DH0：DH1：Dtttest计算如下：在命令窗口输入：>>x=[0.20>>y=[0.10>>hHvartest(X,V,Hvartest(X,V,X为输入样本，V为待检验方差值。alpha为显著性水平，tail13-2013-20单样本方差检验：tailtail【例13-23】 螺钉的直径服从方差为0.0002(cm2)的正态分布。今从产品中随机抽取10只进量，得1.180.0002(cm2取α解：H020.0002vartest计算如下。在命令窗口中输入：>>x=[1.191.211.211.181.171.201.201.171.19>>>>h通常采用F统计量检验两独立样本的方差是否相同。中采用vartest2实现，Hvartest2(X,Hvartest2(X,Y,Hvartest2(X,Y,alpha,tail13-2113-21两样本方差检验：tailtailVxVxVx13-24】6。有两台机床加工同一种零件，这两台机床生产的零件尺寸服从正态分布今从两台机床生产的零件中分别抽取11个和9个零件进量，问甲机床的加工精度是否比乙机床的加工精度差？（取解：H0：12H0：12用vartest2计算，在命令窗口中输入>>x=[6.25.76.56.06.35.85.76.06.05.8>>y=[5.65.95.65.75.86.05.55.7>>h0.0513.6方差分单因素试验的方差分析是指试验中只有一个因素发生改变。在中，单因anova1anova1来自相同的总体或来自具有相同的均值的不同总体的p值，具体调用格式如下。回p值，ANOVA表格和矩阵X各列的箱线图。groupX样本的箱线图。groupXgroupXpanova1(X,group,'dis用dis yopt确定是否要作出ANOVA的表格和样本的箱线图。当dis yopt'on'(默 yept'off'时不给出表格和箱线图。sourcedfFMS2）X【例13-25】方差分析综合实例1。将抗生素注入会产生抗生素与血浆蛋白结合13-225种常用的抗生素注入到牛的体内时，抗生素与13-22解： 命令窗口中输入group=['青霉素';'四环素';'链霉素';'红霉素';'氯霉素>>pp0.0020.01，故可以判断这些抗生素的均值具有显著差异。ANOVA表格和X箱线图如图13-2和图13-3所示。13-2ANOVA5青霉 四环 链霉 红霉 氯霉13-3Xanova2Panova2(X,一个的观察点，则变量reps每一单元观察点的数目，每一单元包含reps行。repsreps默认值为1，此时anova2返回的p值向量包含两个值，分别为原ApB作用的行向量来自一个总体的p值。reps1p33ABp13-262（因解：用anova2函数求解， 命令窗口中输入>>>>p 13-4ANOVA13-4ANOVAppB因素（温度因素）的影响。由于p(1)很小，故可得B因素对产品强度影响显著。p(2)>0.05A（p(3)<0.0113.7回归分 ysis）是研究一个变量Y与其他若干变量X之间相关关系的yx呈现较明显的二次（或高次）函数关系，则可以选polyfit是多项式系数估计函数，其中用参数n设定多项式的最高次幂。polyvalpolyconfpolyfit p返回多项式系数的估计值，n设定多项式的最高次数。x，y Sx的范德蒙矩阵的QR分解的R用来作为polyval函数的输入，从而估计误差。polyval YXp[Y,DELTA]polyval(p,X,S) p，S为polyfit的输出，DELTA为误差估计。性回归模型中，Y±DELTA以50%的概率包含函数在X处的真值。polyconf[YDELTApolyconf(PXS)95%Y[YDELTApolyconf(PXSalphavalue1)value1设定显著性水平αpolytool(x,y,n)polytool(x,y,n,alpha)默认值为0.05。13-27】17297两人测量其旋转定向能力以对这种运动能力的影响先得到的一组数据如表24表13-24与旋转定向能力数24 中编写程序ex1327.m来实现，具体代码如下所示。%ex1327.m%%%%%y的拟合值及预测值y的置信半径deltaypx0>>p 8.9782-yColumns1through22.524326.058227.989628.318627.045024.168922.524326.0582Columns11through28.318627.045024.1689polytool13-5

13-5 统计工具箱中使用函数regress实现多元线性回归。具体调用格式为bregress(Y,[b,bint,r,rint,stats]regress(Y,X,YXXXx2kx1ky1Yy2 1xnk nyk=1alpha为显著性水平(默认时设定为0.05)；输出向量为b，bint为回归系数估计R2RF统计量值，第三个是与统计量F对应的概率P，当P<α时H0，回归模型成立。rcoplot（r，rint）13-28COD（y）（x1（x2（x3y13-25y解：通过在中编写程序ex1328.m来实现，具体代码如下所示%ex1328.mx1=[1.376,x2=[0.450,x3=[2.170x4=[0.8922,1.1610,0.5346,0.9589,1.0239,1.0499,1.1065,y=[5.19,5.30,5.60,5.82,6.00,6.06,6.45,saveex1328_datax1x2x3x4load >>bbintstats0.9846 statsR20.9846，F47.9654，P0.0123中采用nlinfit实现，其调用格式如下。（1）betanlinfit(X,Y,fun,返回非线性回归方程系数的最小二乘估计值。非线性方程由fun给定，fun为用户提供形如Yf(beta,X)的函数。beta0为系数初值。（2）[beta,r,J]nlinfit(X,Y,fun,betaRJacobiJnlparcinlinfit函数的输出beta，r，J。函数具体调用格式为：cinlparci(beta,r,beta95%beta100(1-alpha)%nlpredci计算，具体调用格式如（1）[ypred,delta]nlpredci(FUN,inputs,beta,r,inputsypred95%的置信区间为[ypred-delta,（2）[ypred,delta]nlpredci(FUN,inputs,beta,r,J,返回模型在对应inputs处输出预测值ypred以及其100(1alpha)%的置信区间[ypreddelta,ypreddelta]。nlintoolnlintool(X,Y,fun,X、Y95%nlintool(X,Y,fun,beta0,100(1alpha)%nlintool(X,Y,fun,beta0,alpha,'xname',x、y13-29】yx…8时，yxya(0.4944a、b95%置信区间。并显示出拟合曲线。此处设a、b的初值分别为0.30，0.02。解：funex1329.mfunctionyhat=funex1329(beta0,x)M文件ex1329.m%ex1329.mclear8.0010.0010.0010.0010.0012.0012.0012.0014.0020.0020.0020.0022.0028.0030.0030.0030.000.490.480.470.480.470.460.450.430.450.430.440.430.430.460.420.420.430.410.410.400.420.400.400.410.410.410.400.400.400.380.410.400.400.410.380.400.40beta0=[0.30 %95%生产时间氯气级分>>beta ci a、b0.3896，0.101195%置信区间分别为[0.3805,0.3988]和[0.0771,0.1251]。所以用最小二乘法对该模型系数进行估计，最终的模型表达式y0.38960.1004e013-6 氯