matlab与统计回归分析

1、一Matlab作方差分析方差分析是分析试验（或观测）数据的一种统计方法。在工农业生产和科学研究中，经常要分析各种因素及因素之间的交互作用对研究对象某些指标值的影响。在方差分析中，把试验数据的总波动（总变差或总方差）分解为由所考虑因素引起的波动（各因素的变差）和随机因素引起的波动（误差的变差），然后通过分析比较这些变差来推断哪些因素对所考察指标的影响是显著的，哪些是不显著的。【例1】（单因素方差分析）一位教师想要检查3种不同的教学方法的效果，为此随机地选取水平相当的15位学生。把他们分为3组，每组5人，每一组用一种方法教学，一段时间以后，这位教师给15位学生进行统考，成绩见下表1。问这3种教学方

2、法的效果有没有显著差异。表1学生统考成绩表方法成绩甲7562715873乙7185689290丙7379607581Matlab中可用函数anova1（.）函数进行单因子方差分析。调用格式：p=anova1（X）含义：比较样本mxn的矩阵X中两列或多列数据的均值。其中，每一列表示一个具有m个相互独立测量的独立样本。返回：它返回X中所有样本取自同一总体（或者取自均值相等的不同总体）的零假设成立的概率p。解释：若p值接近0(接近程度有解释这自己设定)，则认为零假设可疑并认为至少有一个样本均值与其它样本均值存在显著差异。Matlab程序：Score=7562715873;8185689290;737

3、9607581;P=anova1(Score)输出结果：方差分析表和箱形图ANOVATableSourceSSdfMSFProbFColumns604.93332302.46674.25610.040088Error852.81271.0667Total1457.733314ColumnNumber由于p值小于0.05，拒绝零假设，认为3种教学方法存在显著差异。例2（双因素方差分析）为了考察4种不同燃料与3种不同型号的推进器对火箭射程（单位：海里）的影响，做了12次试验，得数据如表2所示。表2燃料-推进器-射程数据表推进器1推进器2推进器3燃料158.256.265.3燃料249.154.15

4、1.6燃料360.170.939.2燃料475.858.248.7在Matlab中利用函数anova2函数进行双因素方差分析。调用格式：p=anova2（X,reps）含义：比较样本X中两列或两列以上和两行或两行以上数据的均值。不同列的数据代表因素A的变化，不同行的数据代表因素B的变化。若在每个行-列匹配点上有一个以上的观测量，则参数reps指示每个单元中观测量的个数。返回：当reps=1（默认值）时，anova2将两个p值返回到向量p中。HOA:因素A的所有样本（X中的所有列样本）取自相同的总体；HOB:因素B的所有样本（X中的所有行样本）取自相同的总体。当reps1时，anova2还返回第

5、三个p值：H0AB:因素A与因素B没有交互效应。解释：如果任意一个p值接近于0,则认为相关的零假设不成立。Matlab程序：disp1=58.256.265.3;;60.170.939.2;75.858.248.7;p=anova2（disp1,1）输出结果：方差分析表ANOVATableSourceSSdfMSFProbFColumns157.59352.530.430590.73875Rows223.84672111.92330.917430.44912Error731.986121.9967Total1113.416711由于燃料和推进器对应的p值均大于0.05，

6、所以可以接受零假设H0A和H0B,认为燃料和推进器对火箭的射程没有显著影响。例3（双因素方差分析）设火箭的射程在其它条件基本相同时与燃料种类及推进器型号有关。现在考虑4种不同的燃料及3种不同型号的推进器，对于每种搭配个发射了火箭两次，得数据见表3。问各自变量和自变量的交互效应是否对火箭的射程有显著影响？表3燃料-推进器-射程数据表推进器1推进器2推进器3燃料158.256.265.352.641.260.8燃料42.850.548.4燃料360.170.939.258.373.240.7燃料475.858.248.771.551.041.4Matlab程序：disp2

7、=58.252.649.142.860.158.375.871.5;50.570.973.258.251.0;65.360.851.648.439.240.748.741.4;p=anova2(disp2,2)输出结果：方差分析表ANOVATableSourceSSdfMSFProbFColumns370.98082185.49049.39390.003506Rows261.675387.2254.41740.025969Interaction1768.69256294.782114.92886.1511e-005Error236.951219.7458Total263

8、8.298323曰些显著。方差分析上机练习为研究广告的效果考察4种广告方式：当地报纸（paper）、当地广播（radio）店内销售员（people）和店内展示（display）的效果。共设有144个销售点，每种广告随机抽取36个销售点记录销售额，分布在6个地区的144个销售点的销售情况生成的数据集ADS见下表。数据集ADS中有3个变量：AD表示广告的类型、AREA表示地区、SALES表示销售额（单位：千元）请完成以下练习：（1）概括下列数据：用箱形图、条形图直观地呈现四种广告方式下销售量的分布情况；计算四种广告方式下销售量的均值、方差、标准差、最大和最小值；进行单因素方差分析：检验四种广告方式

9、下销售量数据是否服从正态分布，方差是否相等；检验四种广告方式下的销售量是否有显著差异(a=0.01)；若四种广告方式下的销售量有显著差异，指出哪些类型的广告效果有显著的不同？在设计广告效果的试验时，虽然地区差异对销售量的影响并不是我们感兴趣的，但希望排除这一因素的影响。数据集ADS记录了各个销售点所在的地区AREA。试用双因素方差分析方法分析销售数据，并指出广告方式和地区对销售量是否有显著影响(a=0.01,01)?广告方式(AD)与地区(AREA)之间有无交互效应？表ADS数据集中的数据广告方式(变量:AD)销售额(单位：千元)(变量SALES)地区1地区2地区3地区4地区5地区6当地报纸7

10、557767775768194547087657962(paper)68758372638856866575666670868477806276627870当地广播6951907733791006168637673(radio)1006073688374547879837468757053737966815769656565店内销售员6367858087707564406751616450(people)58827862407670777562877077686155427165783783店内展示5261617657336961664165584445(display)41448652606

11、943515058757552616052524563435560参考答案(1)箱形图：boxplot(ads)结果：有异常值。(其它：略)(2)正态性检验Paper：Hist(X1,6)12108642*0556065707580859095频数直方图分布的正态性检验：normplot(X1)NormalProbabilityPlot均服从正态分布。单因素方差分析ANOVATableSourceSSdfMSFProbFColumns5866.083331955.361113.48318.8495e-008Error20303.2222140145.023Total26169.3056143P

12、=8.8495e-008FColumns1444.22225288.84441.95820.089763Rows5866.083331955.361113.25591.5637e-007Interaction11581577.20.523360.92341Erro5083Total26169.3056143从以上分析结果可知：0.05P1=0.0897630.1，地区对检验水平有一定影响，但不显著。P2=1.5637e-0070.010.1，地区和广告方式对销售量无交互效应。二Matlab作回归分析回归分析的相关数学理论可以参见概率论与数理统计教程，下面仅以示例说明如

13、何利用matlab处理回归分析。1.一元线性回归分析【例1】为了了解百货商店销售额x与流通费率（反映商业活动的一个质量指标，指每元商品流转额所分摊的流通费用）y之间的关系，收集了九个商店的有关数据，见下表1试建立流通费率y与销售额x的回归方程。表1销售额与流通费率数据样本点销售额x（万元）流通费率y11.57.024.54.837.53.6410.53.1513.52.7616.52.5719.52.4822.52.3925.52.2【分析】：首先绘制散点图以直观地选择拟合曲线，这项工作可结合相关专业领域的知识和经验进行，有时可能需要多种尝试。选定目标函数后进行线性化变换，针对变换后的线性目标

14、函数进行回归建模与评价，然后还原为非线性回归方程。【Matlab数据处理】：【Stepl】：绘制散点图以直观地选择拟合曲线x二10.513.516.519.522.525.5;y二7.02.32.2;plot（x,y,-o）输出图形见图1。图1销售额与流通费率数据散点图根据图1，初步判断应以幂函数曲线为拟合目标，即选择非线性回归模型，目标函数为：y二axb(bF)沁0；模型方差的估计值2=竺=0.0012。n-2【注】：严格来讲，模型评价工作应在逆线性化变换后进行；但是，若所建立的线性回归方程不理想，则相应的非线性回归方程必定不理想。【St

15、ep3】拟线性化变换求非线性回归方程(若选择为非线性模型)%逆线性化变换A=exp(b(1)B=b(2)运行结果为：A=8.5173；B=-0.4259。即非线性回归方程为：y=8.5173x-o.4259。2多元线性回归模型(pl)：y=卩+卩x+卩x+卩x+8,戸N(0Q2)01122ppU求得经验回归方程：y=P+Px+Px+Px01122pp统计量：总偏差平方和：SST=X(y-y)2，其自由度为f二n-1；iTi=1回归平方和：ssr工(y-y)2，其自由度为f=p；iRi=1残差平方和：SSE=丫(y-y)2，其自由度为f二n-p-1；iiEi=1它们之间有关系：SST=SSR+S

16、SE。多元回归分析的相关数学理论可以参见多元数据分析，下面仅以示例说明如何利用Matlab作多元回归分析。【例2】参见教材P294：10.1牙膏的销售量。【下面只描述运行程序的过程，应该按照规定格式书写报告】。符号说明：x：表示价格差；1x：广告费用；2y：销售量。【Stepl】：绘制散点图以直观地选择拟合曲线clearclcxl=-0.050.250.6000.250.200.l50.05-0.l50.l50.200.l00.400.450.350.300.500.500.40-0.05-0.05-0.l00.200.l00.500.60-0.0500.050.55;x2=5.506.757

17、.255.507.006.506.755.255.256.006.506.257.006.906.807.l07.006.806.506.256.006.507.006.806.806.505.756.80;y=7.388.5l9.527.509.338.288.757.877.l08.007.898.l59.l08.868.908.879.269.008.757.957.657.278.008.508.759.2l8.277.677.939.26;hl=figure;plot(xl,y,+);h2=figure;plot(x2,y,o);图1y对xl的散点图109.598.587.5755.

18、566.577.5图2y对x2的散点图分析图1,可以发现，随着xl的增加，y的值有比较明显的线性增长趋势；分析图2,当x增大时，y有向上弯曲的趋势，可用二次多项式进行逼近；因此可以选择如下方程作为初步的回归模型：y=卩+卩x+卩x+卩x2+8,N(u,02)0112232U【Step2】：模型求解(理论方法：最小二乘法)alpha=0.05;v=ones(length(x1),1)xlx2(x2.八2);b,bint,r,rint,stats=regress(y,v,alpha)计算结果：b=17.32441.3070-3.69560.3486bint=5.728228.92060.68291

19、.9311-7.49890.10770.03790.6594r=-0.0988-0.0795-0.1195-0.04410.4660-0.01330.29120.2735-0.23510.1031-0.40330.17470.0400-0.15040.12840.1637-0.0527-0.1907-0.0870-0.0165-0.1292-0.3002-0.2933-0.1679-0.21770.11160.30350.06930.24740.2270rint=-0.52700.3294;-0.53090.3718;-0.51060.2716;-0.47310.3848;0.08130.85

20、07;-0.46090.4343;-0.13740.7197;-0.08700.6340;-0.59600.1258;-0.32800.5341;-0.81900.0125;-0.26180.6112;-0.40320.4832;-0.59330.2925;-0.32070.5775;-0.28410.6116;-0.48300.3776;-0.62480.2434;-0.53480.3609;-0.44230.4092;-0.56090.3024;-0.71810.1177;-0.72430.1377;-0.55480.2190;-0.64490.2095;-0.29940.5226;-0.

21、10370.7106;-0.37140.5099;-0.18070.6755;-0.18900.6430stats=0.905482.94090.00000.0490【Step3】结果分析回归模型为：y=17.3244+1.3070 x3.6959x+0.3486x2122从结果数据来看，模型整体可用。但也有缺陷，可以改进。【Step4】销售量的预测设需要预测的点为：x=(x,x,00102,x0p)则预测值为=0+0 x+0 x+0101202(c)ijpxp=(XtX)-1,竝(xx)(xx匕,x0ii0jjiji=1j=1x,i=1,2,pkix,i=1,2,nk=1则在x处的区间预测为

22、：0(np1)d*A)(1)(yt(np1)d*A,y+10(la)0I新的回归模型为2【模型改进】：当两个因素是不独立时，引入交叉项xx,12y=0+0 x+0 x+0 x2+0 xx+80112232412alpha=0.05;v二ones(length(xl),l)xlx2(x2.2)(xl.*x2);b,bint,r,rint,stats=regress(y,v,alpha)输出结果：b=29.ll33ll.l342-7.60800.67l2-l.4777bint=l3.70l344.5252;l.977820.2906;-l2.6932-2.5228;0.2538l.0887;-2.

23、85l8-0.l037r=-0.044l;-0.l229;0.0299;-0.0745;0.384l;-0.0472;0.233l;0.0287;-0.066l;0.0297;-0.4372;0.l763;0.0356;-0.l382;0.l027;0.l270;0.0048;-0.l435;-0.l0l6;0.0050;-0.0389;-0.l334;-0.3272;-0.3274;-0.2l02;0.l4l2;0.3250;0.l096;0.2342;0.2455rint=-0.44250.3542;-0.54080.295l;-0.3l0l0.3698;-0.47360.3247;0.0

24、2450.7437;-0.46400.3695;-0.l6740.6337;-0.23690.2943;-0.375l0.2430;-0.369l0.4284;-0.8ll8-0.0627;-0.23060.5832;-0.37880.4499;-0.552l0.2757;-0.3l720.5226;-0.29l70.5456;-0.39440.4039;-0.54900.262l;-0.5l930.3l60;-0.39260.4026;-0.43600.3582;-0.50450.2378;-0.72l20.0667;-0.6326-0.022l;-0.60850.l88l;-0.23980

25、.5223;-0.04840.6984;-0.29880.5l8l;-0.l6500.6335;-0.l39l0.6302stats=0.920972.777l0.00000.0426结果分析：效果更好。3.逐步回归方法要点：【Stepl】根据问题所属专业领域的理论和经验提出对因变量可能有影响的所有自变量；【Step2】计算每一个自变量对因变量的相关系数，按其绝对值从大到小排序；【Step3】取相关系数绝对值最大的那个自变量建立一元线性回归模型，检验所得回归方程的显著性，若检验表明回归效果则转入【Step4】,若检验表明回归效果不显著则停止建模；【Step4】进行变量的追加、剔除和回归方程的更

26、新操作。Matlab命令：【命令1】：stepwisefit【调用格式】：b,se,pval,inmodel,stats,nextstep,history=stepwisefit（x,y,paraml,value1,param2,value2,）【参数说明】：X：p个自变量的n个观测值的nxp矩阵；Y:因变量的n个观测值的nx1矩阵；penter:设置回归方程显著性检验的显著性概率上限，缺省值为0.05；premove：设置回归方程显著性检验的显著性概率下限，缺省值为0.10；display：用来指明是否强制显示建模过程信息，取值为on（显示，缺省设置）和off（不显示）。【例3】某种水泥在凝

27、固时放出的热量（单位：卡/克）Y与水泥中的四种化学成分所占的百分比有关，现测得13组数据如下表：编号XXnXOX4Y172666078.52129155274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.4作回归分析。【Matlab程序】：clearclcloadhaldb,se,pval,inmodel,stats,nextstep,history=stepw

28、isefit(ingredients,heat,penter,0.10,display,off);%自变量的筛选和模型参数估计信息inmodel,b0=ercept,b%回归方程显著性整体检验信息Allp=stats.pval,rmse=stats.rmse%回归方程显著性分别检验信息P=stats.PVAL输出结果：inmodel=1100；b0=52.5773；b=1.46830.66230.2500-0.2365；Allp=4.4066e-009；rmse=2.4063；P=0.00000.00000.20890.2054。结果分析：最优回归方程为y二52.5773+1

29、.4683x+0.6623x，回归方程显著性整体检验12和分别检验均为高度显著，模型标准误差估计为2.4063。【命令2】：stepwise【调用格式】：stepwise（x,y,inmodel,penter,premove）【说明】：创建多元线性回归分析的逐步回归法建模的交互式图形环境。【图形界面说明】：窗口1：CoefficientswitherrorBars绘出各个解释变量回归系数的估计，圆点表示点估计值，横线表示置信区间（有色线段表示90%置信区间，黑色线段表示95%置信区间）。窗口的右侧给出回归系数的点估计值（Coeff）、显著性检验的t统计量的值（t-test）和显著性概率p值（p

30、-val）.窗口2：ModelHistory该窗口绘出的圆点表示历次建模的模型标准差a的估计。两个窗口中间输出的是当前模型的有关信息，包括：Intercept：模型截距（常数项）的估计；RMSE:模型标准差a的估计；R-square:可决系数；Adj-R-sq：校正可决系数；F：模型整体性检验的F统计量的值；p：模型整体性检验的显著性概率。窗口I右侧的三个按钮：NextStep:在回归方程中按相关系数绝对值大小逐次引入解释变量，如无解释变量可引入时，按钮不可用；AllSteps：直接给出“只进不出方式建模的最终结果（注意，此时的回归方程未必是最优回归方程）；Export：选择向Workspac

31、e传输的计算结果（有关变量名可由用户自定义）stepwise（ingredients,heat,1111,0.05,0.10）;三matlab作相关分析一、相关系数要初步研究变量之间的随机性关系，我们就要清楚，研究的对象是二元或多元的随机向量，利用的是成对观测数据。首先绘制一张散点图，直观上大致判断两两变量之间是否存在某种关系。MATLAB命令(散点图)：gscatter-两个变量的散点图.用法：gscatter(x,y)lsline-在散点图上增加最小二乘拟合线.用法：lslinegplotmatrix-矩阵散点图。用法：gplotmatrix(x,y)。其中x,y都是矩阵，行数相同。例如:

32、x=normrnd(0,1,100,3);y=normrnd(1,2,100,2);gplotmatrix(x,y)*;t撫.曦*養*+*:*金A，*+-22-202-202如果认为两个变量之间存在着某种直线关系，我们可以用相关系数来刻画这种关系。首先引入如下样本相关系数的概念：对二元总体(X，Y)的样本(x,yi),i二12,n，定义样本相关系数为苴中S2二1工(X-X)2,S2X0niY0二1工(Y-Y)2ni分别为X和Y的样本方差,i=1i=11,_-s2y=工（X.X）a-Y）叫x与y之间的样本协方差。这是一个重要统计XYniii-i量，与总体相关系数P（X,Y）相对应。那么，怎样充分

33、发挥这个统计量的作用呢？下面我们讲讲如何利用它对总体相关系数p：=p（X,Y）作假设检验和区间估计。原假设为H0：P二；对立假设为H1:P丰0-在原假设成立的情况下，可以证明下面的统计量服从自由度为n-2的t分布:t=r所以给定检验水平a，可得原假设的否定域lltlt（n-Ia/2MATLAB命令：corrcoef用法：r,p,rlo,rup=corrcoef（x）其中：x矩阵；r：相关矩阵；p：p-值;rlo：置信下限；rup:置信上限二、偏相关分析基本描述：控制其它变量的情况下研究两个变量之间的线性关系.因变量自变量J.自变量2*.原假设:两个变量之间的偏相关系数为，、MATLAB命令：p

34、artialcorr用法：RHO,PVAL=PARTIALCORR（X,Z）X是由多个变量的样本值构成的矩阵，Z是由控制变量构成的矩阵，RHO是偏相关系数矩阵，PVAL是对应的p-值。例：财政收入预测问题：财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试对各因素之间的关系进行初步分析。年份国民收工业总农业总总人口就业人固定资产财政入（亿元）产值（亿元）产值（亿元）口（万人）投资（亿元）收入（亿元）19525983494615748220729441841953586455475587962136489216195

35、4707520491602662183297248195573755852961465223289825419568257155566282823018150268195783779857564653237111392861958102812355986599426600256357195911141681509672072617333844419601079187044466207258803805061961757115643465859255901382711962677964461672952511066230196377910465146917226640852661964943125058470499