MATLAB数据输入与分析.ppt

2019年12月5日,MATLAB和R软件,1,MATLAB数据输入与分析,2019年12月5日,MATLAB和R软件,2,引言,统计研究的对象是受随机因素影响的数据数据样本少则几个，多则成千上万，人们希望能用少数几个包含其最多相关信息的数值来体现数据样本总体的规律。面对一批数据如何进行描述与分析，需要掌握一些数理统计的最基本方法。我们将用MATLAB及其统计工具箱(StatisticsToolbox)来实现数据的统计描述和分析。,2019年12月5日,MATLAB和R软件,3,数据输入,两种方法一种是在交互环境中直接输入如果在统计中数据量比较大，这样作不太方便；另一种办法是先把数据写入一个数据文件中，在MATLAB中用相关命令读入数据具体作法：.mat文件.dat或.txt文件.xls文件,2019年12月5日,MATLAB和R软件,4,数据预处理,查找错误查找异常值,数据分析,描述性统计参数估计假设检验方差分析回归分析.,2019年12月5日,MATLAB和R软件,5,数据分析操作实例,数据集2*2因子行为实验数据在一给定试验中按键的反应时间和特性输入数据查找错误与异常值计算描述性统计量进一步数据分析,2019年12月5日,MATLAB和R软件,6,2019年12月5日,MATLAB和R软件,7,计算统计量的函数,样本均值: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,2019年12月5日,MATLAB和R软件,8,计算统计量的函数,样本方差: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,2019年12月5日,MATLAB和R软件,9,计算统计量的函数,样本标准差: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,2019年12月5日,MATLAB和R软件,10,计算统计量的函数,协方差矩阵: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,2019年12月5日,MATLAB和R软件,11,计算统计量的函数,相关系数: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,2019年12月5日,MATLAB和R软件,12,参数估计与假设检验,最大似然估计和区间估计: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,正态分布的均值和均方差的最大似然估计,2019年12月5日,MATLAB和R软件,13,参数估计与假设检验,单总体的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(.),例某车间用一台包装机包装糖果。包得的袋装糖重是一个随机变量，它服从正态分布。当机器正常时，其均值为0.5公斤，标准差为0.015公斤。某日开工后为检验包装机是否正常，随机地抽取它所包装的糖9袋，称得净重为（公斤）：0.4970.5060.5180.5240.4980.5110.5200.5150.512,问机器是否正常？,解：总体已知，XN(,.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的水平下，可拒绝原假设，即认为这天包装机工作不正常。,2019年12月5日,MATLAB和R软件,14,参数估计与假设检验,单总体的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(.),例某种电子元件的寿命x(以小时计)服从正态分布,2均未知.现得16只元件的寿命如下:159280101212224379179264222362168250149260485170问是否有理由认为元件的平均寿命大于225(小时)?,解按题意需检验H0:225，取=0.05。Matlab实现如下：x=159280101212224379179264.222362168250149260485170;h,p,ci=ttest(x,225,0.05,1)求得h=0，p=0.2570，说明在显著水平为0.05的情况下，不能拒绝原假设，认为元件的平均寿命不大于225小时。,2019年12月5日,MATLAB和R软件,15,参数估计与假设检验,双总体的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(.),例在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的得率,试验是在同一平炉上进行的。每炼一炉钢时除操作方法外,其它条件都可能做到相同。先用标准方法炼一炉,然后用建议的新方法炼一炉,以后交换进行,各炼了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:1X=binocdf(0:20,20,0.5);H,p,k,c=kstest(X,0.05,0)H=1p=2.8173e-005k=0.5000c=0.2873,X服从二项分布,不服从标准正态分布,拒绝原假设,2019年12月5日,MATLAB和R软件,17,非参数假设检验,例下面列出了84个伊特拉斯坎（Etruscan）人男子的头颅的最大宽度（mm），试检验这些数据是否来自正态总体（取=0.1)。141148132138154142150146155158150140147148144150149145149158143141144144126140144142141140145135147146141136140146142137148154137139143140131143141149148135148152143144141143147146150132142142143153149146149138142149142137134144146147140142140137152145,min(x),max(x)%求数据中的最小数和最大数hist(x,8)%画直方图fi=length(find(x=135,求得皮尔逊统计量chisum=1.9723，0.12(7-2-1)=0.12(4)=7.7794，故在水平0.1下接受H0，即认为数据来自正态分布总体。,2019年12月5日,MATLAB和R软件,18,非参数假设检验,双样本K-S检验:kstest2调用格式:(1)H=kstest2(X1,X2)(2)H=kstest2(X1,X2,alpha)(3)H=kstest2(X1,X2,alpha,tail)(4)H,p,ksstat,cv=kstest2(.),X=-1:1:5;Y=randn(20,1);H,p,k=kstest2(X,Y)H=1p=0.0219k=0.6143,Y是由randn生成的正态分布随机数,拒绝X和Y具有相同的分布的假设,2019年12月5日,MATLAB和R软件,19,非参数假设检验,符号检验:signtest调用格式:(1)p=signtest(X)(2)p=signtest(X,m)(3)p=signtest(X,Y)(4)p=signtest(.,alpha)(5)p,H=signtest(.)(6)p,H,stats=signtest(.),X=normrnd(0,1,20,1);Y=normrnd(0,2,20,1);p,H=signtest(X,Y,0.05)p=0.8238H=0,不能拒绝原假设,2019年12月5日,MATLAB和R软件,20,非参数假设检验,秩和检验:ranksum调用格式:(1)p=ranksum(X,Y)(2)p=ranksum(X,Y,alpha)(3)p,H=ranksum(.)(4)p,H,stats=ranksum(.),例7某商店为了确定向公司A或公司B购买某种产品，将A,B公司以往各次进货的次品率进行比较，数据如下所示，设两样本独立。问两公司的商品的质量有无显著差异。设两公司的商品的次品的密度最多只差一个平移，取=0.05。A：7.03.59.68.16.25.110.44.02.010.5B：5.73.24.211.09.76.93.64.85.68.410.15.512.3,解分别以A、B记公司A、B的商品次品率总体的均值。所需检验的假设是H0:A=B，H1:AB.Matlab实现如下：a=7.03.59.68.16.25.110.44.02.010.5;b=5.73.24.211.09.76.93.64.85.68.410.15.512.3;p,h=ranksum(a,b)求得p=0.8041，h=0，表明两样本总体均值相等的概率为0.8041，并不很接近于零，且h=0说明可以接受原假设，即认为两个公司的商品的质量无明显差异。,2019年12月5日,MATLAB和R软件,21,方差分析,单因素方差分析:anova1调用格式:(1)p=anova1(X)(2)p=anova1(X,group)(3)p=anova1(X,group,displayopt)(4)p,table=anova1(.)(5)p,table,stats=anova1(.),X=2.16503.69611.55383.64004.95511.62682.05912.29883.86444.20111.07513.79714.24602.65074.23481.35162.26412.36102.72965.86170.30352.87173.57744.98464.9438;p=anova1(X)p=5.9952e-005,2019年12月5日,MATLAB和R软件,22,方差分析,方差分析,2019年12月5日,MATLAB和R软件,23,例用4种工艺生产灯泡，从各种工艺制成的灯泡中各抽出了若干个测量其寿命，结果如下表，试推断这几种工艺制成的灯泡寿命是否有显著差异。,解编写程序如下：x=1620158014601500167016001540155017001640162016101750172016801800;x=x(1:4),x(16),x(5:8),x(9:11),x(12:15);g=ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4);p=anova1(x,g),求得0.01X=5.54.53.55.54.54.06.04.03.06.55.04.07.05.55.07.05.04.5;p=anova2(X,3)p=0.00000.00010.7462,2019年12月5日,MATLAB和R软件,25,方差分析,方差分析,2019年12月5日,MATLAB和R软件,26,例一火箭使用了4种燃料，3种推进器作射程试验，每种燃料与每种推进器的组合各发射火箭2次，得到结果如下：,试在水平0.05下，检验不同燃料（因素A）、不同推进器（因素B）下的射程是否有显著差异？交互作用是否显著？,方差分析,2019年12月5日,MATLAB和R软件,27,解编写程序如下：clc,clearx0=58.2,52.656.2,41.265.3,60.849.1,42.854.1,50.551.6,48.460.1,58.370.9,73.239.2,40.775.8,71.558.2,51.048.7,41.4;x1=x0(:,1:2:5);x2=x0(:,2:2:6);fori=1:4x(2*i-1,:)=x1(i,:);x(2*i,:)=x2(i,:);endp=anova2(x,2)求得p=0.00350.02600.001，表明各试验均值相等的概率都为小概率，故可拒绝均值相等假设。即认为不同燃料（因素A）、不同推进器（因素B）下的射程有显著差异，交互作用也是显著的。,2019年12月5日,MATLAB和R软件,28,回归分析,线性回归:regress调用格式:(1)B=regress(Y,X)(2)B,BINT=regress(Y,X)(3)B,BINT,R=regress(Y,X)(4)B,BINT,R,RINT=regress(Y,X)(5)B,BINT,R,RINT,stats=regress(Y,X)(6).=regress(Y,X,alpha),2019年12月5日,MATLAB和R软件,29,回归分析,X=111213141516171819110;Y=11.116512.062713.007514.0352,14.930316.169617.005918.179719.026420.0872;b,bint=regress(Y,X,0.05)b=10.04561.0030bint=9.916510.17470.98221.0238,回归分析,2019年12月5日,MATLAB和R软件,30,例合金的强度y与其中的碳含量x有比较密切的关系，今从生产中收集了一批数据如下表：x0.100.110.120.130.140.150.160.170.18y42.041.545.045.545.047.549.055.050.0试先拟合一个函数y(x)，再用回归分析对它进行检验。,解先画出散点图：x=0.1:0.01:0.18;y=42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0;plot(x,y,+)可知y与x大致上为线性关系。,用regress和rcoplot编程如下：clc,clearx1=0.1:0.01:0.18;y=42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0;x=ones(9,1),x1;b,bint,r,rint,stats=regress(y,x);b,bint,stats,rcoplot(r,rint),设回归模型为y=0+1x,得到b=27.4722137.5000bint=18.685136.259475.7755199.2245stats=0.798527.74690.0012即0=27.4722，1=137.5000，0的置信区间是18.6851,36.2594，1的置信区间是75.7755,199.2245；R2=0.7985，F=27.7469，p=0.0012。,可知模型成立,观察命令rcoplot(r,rint)所画的残差分布，除第8个数据外其余残差的置信区间均包含零点，第8个点应视为异常点，将其剔除后重新计算，可得b=30.7820109.3985bint=26.280535.283476.9014141.8955stats=0.918867.85340.0002应该用修改后的这个结果。,2019年12月5日,MATLAB和R软件,31,例:若x=2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,y=6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76;用进行拟合.,回归分析,非线性回归:nlinfit调用格式:(1)beta=nlinfit(X,Y,fun,beta0)(2)beta,R,J=nlinfit(X,Y,fun,beta0),先创建函数:functiony=f(beta,x)a=beta(1);b=beta(2);y=x./(a.*x+b);,2019年12月5日,MATLAB和R软件,32,回归分析,在命令窗口输入:x=2:16;y=6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,.10.59,10.6,10.8,10.6,10.9,10.76;beta0=0.1,0.1;beta=nlinfit(x,y,f,beta0)beta=0.08450.1152,拟合函数是,回归分析,2019年12月5日,MATLAB和R软件,33,例在研究化学动力学反应过程中，建立了一个反应速度和反应物含量的数学模型，形式为,其中1,5是未知的参数，x1,x2,x3是三种反应物（氢，n戊烷，异构戊烷）的含量，y是反应速度。今测得一组数据如下表，试由此确定参数1,5，并给出其置信区间。1,5的参考值为（0.1，0.05，0.02，1，2）。,序号反应速度y氢x1n戊烷x2异构戊烷x318.554703001023.79285801034.8247030012040.024708012052.754708010614.391001901072.541008065,序号反应速度y氢x1n戊烷x2异构戊烷x384.3547019065913.0010030054108.50100300120110.05100801201211.3228530010133回归分析,2019年12月5日,MATLAB和R软件,34,解首先，以回归系数和自变量为输入变量，将要拟合的模型写成函数文件huaxue.m：functionyhat=huaxue(beta,x);yhat=(beta(4)*x(2)-x(3)/beta(5)./(1+beta(1)*x(1)+.beta(2)*x(2)+beta(3)*x(3);然后，用nlinfit计算回归系数，用nlparci计算回归系数的置信区间，用nlpredci计算预测值及其置信区间,clc,clearx0=18.554703001023.79285801034.8247030012040.024708012052.754708010614.391001901072.54100806584.3547019065913.0010030054108.50100300120110.05100801201211.3228530010133x=x0(:,3:5);y=x0(:,2);beta=0.1,0.05,0.02,1,2;%回归系数的初值betahat,f,j=nlinfit(x,y,huaxue,beta);%f,j是下面命令用的信息betaci=nlparci(betahat,f,j);betaa=betahat,betaci%回归系数及其置信区间yhat,delta=nlpredci(huaxue,x,betahat,f,j),用nlintool得到一个交互式画面，左下方的Export可向工作区传送数据，如剩余标准差等。使用命令nlintool(x,y,huaxue,beta)可看到画面，并传出剩余标准差rmse=0.1933。,回归分析,2019年12月5日,MATLAB和R软件,35,例某厂生产的一种电器的销售量y与竞争对手的价格x1和本厂的价格x2有关。下表是该商品在10个城市的销售记录。,x1(元)120140190130155175125145180150 x2(元)10011090150210150250270300250y(个)10210012077469326696585,试根据这些数据建立y与x1和x2的关系式，对得到的模型和系数进行检验。若某市本厂产品售价160（元），竞争对手售价170（元），预测商品在该市的销售量。,解分别画出y关于x1和y关于x2的散点图，可以看出y与x2有较明显的线性关系，而y与x1之间的关系则难以确定，我们将作几种尝试，用统计分析决定优劣。设回归模型为y=0+1x1+2x2,回归分析,2019年12月5日,MATLAB和R软件,36,编写如下程序：x1=120140190130155175125145180150;x2=10011090150210150250270300250;y=10210012077469326696585;x=ones(10,1),x1,x2;b,bint,r,rint,stats=regress(y,x);b,bint,stats得到b=66.51760.4139-0.2698bint=-32.5060165.5411-0.20181.0296-0.461