6MATLAB2014a数理统计

MATLAB数理统计引言•数理统计研究的对象是受随机因素影响的数据•数据样本少则几个，多则成千上万，人们希望能用少数几个包含其最多相关信息的数值来体现数据样本总体的规律。•面对一批数据如何进行描述与分析，需要掌握参数估计和假设检验这两个数理统计的最基本方法。•我们将用MATLAB的统计工具箱(StatisticsToolbox)来实现数据的统计描述和分析。频数表和直方图•一组数据（样本）往往是杂乱无章的，作出它的频数表和直方图，可以看作是对这组数据的一个初步整理和直观描述。•将数据的取值范围划分为若干个区间，然后统计这组数据在每个区间中出现的次数，称为频数，由此得到一个频数表。以数据的取值为横坐标，频数为纵坐标，画出一个阶梯形的图，称为直方图，或频数分布图。•若样本容量不大，能够手工作出频数表和直方图，当样本容量较大时则可以借助MATLAB这样的软件了。例1学生的身高和体重学校随机抽取100名学生，测量他们的身高和体重，所得数据如表身高体重17275171621666216055155571735816655170631675317360身高体重16955168671686517567176641685016149169631716117864身高体重16964165521645917374172691695217357173611667016357身高体重17165169621705817264169581677217576164591666316954身高体重16747168651656416857176571705715851165621725316966数据输入•两种方法•一种是在交互环境中直接输入•如果在统计中数据量比较大，这样作不太方便；•另一种办法是先把数据写入一个纯文本数据文件data.txt中，数据列之间用空格键或Tab键分割，该数据文件data.txt存放在matlab\work子目录下，在MATLAB中用load命令读入数据•具体作法是：loaddata.txt•这样在内存中建立了一个变量data，它是一个数据矩阵。作频数表及直方图•用hist命令实现，其用法是：[N,X]=hist(Y,M)•数组（行、列均可）Y的频数表•它将区间[min(Y),max(Y)]等分为M份（缺省时M设定为10），N返回M个小区间的频数，X返回M个小区间的中点。•hist(Y,M)•数组Y的直方图。例1的程序loaddata.txt;high=data(:,1:2:9);high=high(:);weight=data(:,2:2:10);weight=weight(:);[n1,x1]=hist(high)[n2,x2]=hist(weight)subplot(1,2,1)hist(high)subplot(1,2,2)hist(weight)例2统计字符出现的频数统计下列五行字符串中字符a、g、c、t出现的频数1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggagg2.cggaggacaaacgggatggcggtattggaggtggcggactgttcgggga3.gggacggatacggattctggccacggacggaaaggaggacacggcggacataca4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagctta5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggc例2的程序clcfid1=fopen('shuju.txt','r');i=1;while(~feof(fid1))data=fgetl(fid1);a=length(find(data==97));b=length(find(data==99));c=length(find(data==103));d=length(find(data==116));e=length(find(data=97&data=122));f(i,:)=[abcdea+b+c+d];i=i+1;endfhe=[sum(f(:,1))sum(f(:,2))sum(f(:,3))sum(f(:,4))...sum(f(:,5))sum(f(:,6))]fid2=fopen('pinshu.txt','w');fprintf(fid2,'%8d%8d%8d%8d%8d%8d\n',f');fclose(fid1);fclose(fid2);把五行字符复制到一个纯文本数据文件shuju.txt中，放在matlab\work子目录下统计结果最后写到一个纯文本文件pinshu.txt中。在程序中多引进了几个变量，是为了检验字符串是否只包含a、g、c、t四个字符概率分布离散型随机变量:离散均匀分布二项分布泊松分布几何分布超几何分布负二项分布连续型随机变量:连续均匀分布指数分布正态分布对数正态分布χ2分布非中心χ2分布t分布非中心t分布F分布非中心F分布β分布γ分布Rayleigh分布Weibull分布统计中常用的几个概率分布),(12N正态分布））分布（）Chisquare223）t分布（student）4）F分布MATLAB统计工具箱(toolbox\stats)中的概率分布normchi2tf正态分布χ2分布t分布F分布pdfcdfinvstatrnd概率密度概率分布逆概率分布均值与方差随机数生成使用方法：将分布命令字符与函数命令字符接起来，并输入自变量（可以是标量、数组或矩阵）和参数p=normpdf(x,mu,sigma);均值mu，标准差sigma的正态分布在x的密度函数p=p(x)（mu=0，sigma=1时可缺省）p=tcdf(x,n);t分布(自由度n)在x的分布函数p=F(x)x=chi2inv(p,n);χ2分布(自由度n)使分布函数F(x)=p的x（即p分位数）x=norminv(0.9,0,2)x=2.5631p=normcdf(2.5631,0,2)p=0.9000[m,v]=fstat(n1,n2);F分布(自由度n1,n2)的均值m和方差v概率分布实例计算统计量的函数样本均值:mean调用格式:(1)Y=mean(X)(2)Y=mean(X,dim)指定维例:X=[012;345];Y=mean(X)Y=1.50002.50003.5000Y=mean(X,2)Y=14计算统计量的函数样本方差:var调用格式:(1)Y=var(X)(2)Y=var(X,1)用n-1标准化用n标准化X=[4-21;957];Y=var(X)Y=12.500024.500018.0000Y=var(X,1)Y=6.250012.25009.0000计算统计量的函数样本标准差:std调用格式:(1)Y=std(X)(2)Y=std(X,1)X=[4-21;957];Y=std(X)Y=3.53554.94974.2426Y=std(X,1)Y=2.50003.50003.0000计算统计量的函数协方差矩阵:cov调用格式:(1)C=cov(X)(2)C=cov(X,Y)(3)C=cov(X,1)(4)C=cov(X,Y,1)X=[-112;-231;403];Y=cov(X)Y=10.3333-4.16673.0000-4.16672.3333-1.50003.0000-1.50001.0000计算统计量的函数相关系数:corrcoef调用格式:R=corrcoef(X)中心矩:moment调用格式:M=moment(X,order)X=[-112;-231;403];Y=corrcoef(X)Y=1.0000-0.84860.9333-0.84861.0000-0.98200.9333-0.98201.0000Y=moment(X,3)Y=11.40740.74070参数估计与假设检验最大似然估计和区间估计:mle调用格式:(1)phat=mle(dist,data)(2)[phat,pci]=mle(dist,data,alpha)data=[0.9501,0.2311,0.6068,0.4860,...0.8913,0.7621,0.4565,0.0185,0.8214,0.4447];phat=mle('normal',data)phat=0.56690.2835正态分布的均值和均方差的最大似然估计参数估计与假设检验单总体的U检验:ztest调用格式:(1)H=ztest(X,m,sigma)(2)H=ztest(X,m,sigma,alpha)(3)H=ztest(X,m,sigma,alpha,tail)(4)[H,p]=ztest(...)(5)[H,p,ci]=ztest(...)例3某车间用一台包装机包装糖果。包得的袋装糖重是一个随机变量，它服从正态分布。当机器正常时，其均值为0.5公斤，标准差为0.015公斤。某日开工后为检验包装机是否正常，随机地抽取它所包装的糖9袋，称得净重为（公斤）：0.4970.5060.5180.5240.4980.5110.5200.5150.512,问机器是否正常？解：总体σ已知，X~N(μ,.0152)，μ未知。于是提出假设H0:μ=0.500和H1:μ≠0.5。Matlab实现如下：x=[0.4970.5060.5180.5240.498...0.5110.5200.5150.512];[h,p,ci]=ztest(x,0.5,0.015)求得h=1，p=0.0248，说明在0.05的水平下，可拒绝原假设，即认为这天包装机工作不正常。参数估计与假设检验单总体的t检验:ttest调用格式:(1)H=ttest(X)(2)H=ttest(X,m)(3)H=ttest(...,alpha)(4)H=ttest(...,tail)(5)[H,p]=ttest(...)(6)[H,p,ci]=ttest(...)例4某种电子元件的寿命x(以小时计)服从正态分布,μ,σ2均未知.现得16只元件的寿命如下:159280101212224379179264222362168250149260485170问是否有理由认为元件的平均寿命大于225(小时)?解按题意需检验H0:μμ0=225,H1:μ225，取α=0.05。Matlab实现如下：x=[159280101212224379179264...222362168250149260485170];[h,p,ci]=ttest(x,225,0.05,1)求得h=0，p=0.2570，说明在显著水平为0.05的情况下，不能拒绝原假设，认为元件的平均寿命不大于225小时。参数估计与假设检验双总体的t检验:ttest2调用格式:(1)H=ttest2(X,Y)(2)H=ttest2(X,Y,alpha)(3)H=ttest2(X,Y,alpha,tail)(4)[H,p]=ttest2(...)(5)[H,p,ci]=ttest2(...)例5在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的得率,试验是在同一平炉上进行的。每炼一炉钢时除操作方法外,其它条件都可能做到相同。先用标准方法炼一炉,然后用建议的新方法炼一炉,以后交换进行,各炼了10炉,其得率分别为1°标准方法78.172.476.274.377.478.476.075.676.777.32°新方法79.181.077.379.180.079.179.177.380.282.1设这两个样本相互独立且分别来自正态总体N(μ1,σ2)和N(μ2,σ2)，μ1,μ2,σ2均未知，问建议的新方法能否提高得率?(取α=0.05)解（i）需要检验假设H0:μ1=μ2，H1:μ1μ2(ii)Matlab实现x=[78.172.476.274.377.478.476.075.6