matlab与统计回归分析.doc

一 Matlab作方差分析方差分析是分析试验（或观测）数据的一种统计方法。在工农业生产和科学研究中，经常要分析各种因素及因素之间的交互作用对研究对象某些指标值的影响。在方差分析中，把试验数据的总波动（总变差或总方差）分解为由所考虑因素引起的波动（各因素的变差）和随机因素引起的波动（误差的变差），然后通过分析比较这些变差来推断哪些因素对所考察指标的影响是显著的，哪些是不显著的。【例1】（单因素方差分析）一位教师想要检查3种不同的教学方法的效果，为此随机地选取水平相当的15位学生。把他们分为3组，每组5人，每一组用一种方法教学，一段时间以后，这位教师给15位学生进行统考，成绩见下表1。问这3种教学方法的效果有没有显著差异。表1 学生统考成绩表方法成绩甲7562715873乙7185689290丙7379607581Matlab中可用函数anova1()函数进行单因子方差分析。调用格式：p=anova1(X)含义：比较样本 mn的矩阵X中两列或多列数据的均值。其中，每一列表示一个具有m个相互独立测量的独立样本。返回：它返回X中所有样本取自同一总体（或者取自均值相等的不同总体）的零假设成立的概率p。解释：若p值接近0（接近程度有解释这自己设定），则认为零假设可疑并认为至少有一个样本均值与其它样本均值存在显著差异。Matlab程序：Score=75 62 71 58 73;81 85 68 92 90;73 79 60 75 81;P=anova1(Score)输出结果：方差分析表和箱形图ANOVA TableSource SS df MS F ProbFColumns 604.93332302.46674.25610.040088Error 852.8 1271.0667Total 1457.7333 14由于p值小于0.05，拒绝零假设，认为3种教学方法存在显著差异。例2（双因素方差分析）为了考察4种不同燃料与3种不同型号的推进器对火箭射程（单位：海里）的影响，做了12次试验，得数据如表2所示。表2 燃料-推进器-射程数据表推进器1推进器2推进器3燃料158.256.265.3燃料燃料360.170.939.2燃料475.858.248.7在Matlab中利用函数 anova2函数进行双因素方差分析。调用格式：p=anova2(X,reps)含义：比较样本X中两列或两列以上和两行或两行以上数据的均值。不同列的数据代表因素A的变化，不同行的数据代表因素B的变化。若在每个行-列匹配点上有一个以上的观测量，则参数reps指示每个单元中观测量的个数。返回：当 reps=1（默认值）时，anova2将两个p值返回到向量p中。 H0A：因素A的所有样本（X中的所有列样本）取自相同的总体； H0B：因素B的所有样本（X中的所有行样本）取自相同的总体。 当reps1时，anova2还返回第三个p值： H0AB：因素A与因素B没有交互效应。解释：如果任意一个p值接近于0，则认为相关的零假设不成立。Matlab程序：disp1=58.2 56.2 65.3;49.1 54.1 51.6;60.1 70.9 39.2;75.8 58.2 48.7;p=anova2(disp1,1)输出结果：方差分析表ANOVA TableSource SS df MS F ProbFColumns 157.59 3 52.53 0.430590.73875Rows 223.84672 111.9233 0.917430.44912Error 731.98612 1.9967Total 1113.416711由于燃料和推进器对应的p值均大于0.05，所以可以接受零假设H0A和H0B，认为燃料和推进器对火箭的射程没有显著影响。例3（双因素方差分析）设火箭的射程在其它条件基本相同时与燃料种类及推进器型号有关。现在考虑4种不同的燃料及3种不同型号的推进器，对于每种搭配个发射了火箭两次，得数据见表3。问各自变量和自变量的交互效应是否对火箭的射程有显著影响？表3 燃料-推进器-射程数据表推进器1推进器2推进器3燃料158.252.656.241.265.360.8燃料249.142.854.150.551.648.4燃料360.158.370.9燃料475.871.558.251.048.741.4Matlab程序：disp2=58.2 52.6 49.1 42.8 60.1 58.3 75.8 71.5;56.2 41.2 54.1 50.5 70.9 73.2 58.2 51.0;65.3 60.8 51.6 48.4 39.2 40.7 48.7 41.4;p=anova2(disp2,2)输出结果：方差分析表ANOVA TableSource SS df MS F ProbFColumns 370.98082185.49049.39390.003506Rows 261.675 387.225 4.41740.025969Interaction1768.69256294.782114.92886.1511e-005Error 236.95 1219.7458Total 2638.2983 23显著。方差分析上机练习为研究广告的效果，考察4种广告方式：当地报纸（paper)、当地广播（radio）、店内销售员（people）和店内展示（display）的效果。共设有144个销售点，每种广告随机抽取36个销售点记录销售额，分布在6个地区的144个销售点的销售情况生成的数据集ADS见下表。数据集ADS中有3个变量：AD表示广告的类型、AREA表示地区、SALES表示销售额（单位：千元）。请完成以下练习：(1) 概括下列数据：用箱形图、条形图直观地呈现四种广告方式下销售量的分布情况；计算四种广告方式下销售量的均值、方差、标准差、最大和最小值；(2) 进行单因素方差分析：检验四种广告方式下销售量数据是否服从正态分布，方差是否相等；检验四种广告方式下的销售量是否有显著差异（）；若四种广告方式下的销售量有显著差异，指出哪些类型的广告效果有显著的不同？(3) 在设计广告效果的试验时，虽然地区差异对销售量的影响并不是我们感兴趣的，但希望排除这一因素的影响。数据集ADS记录了各个销售点所在的地区AREA。试用双因素方差分析方法分析销售数据，并指出广告方式和地区对销售量是否有显著影响（）？广告方式(AD)与地区(AREA)之间有无交互效应？表 ADS数据集中的数据广告方式（变量:AD）销售额（单位：千元）（变量SALES）地区1地区2地区3地区4地区5地区6当地报纸(paper)75 57 7668 75 8377 75 7266 66 7676 81 6370 86 6294 54 7088 56 8687 65 6584 77 7879 62 7580 62 70当地广播(radio)69 51 10054 78 7990 77 6083 74 6933 79 7368 75 65100 61 6870 53 7368 63 8379 66 6576 73 7481 57 65店内销售员（people）63 67 8558 82 7880 87 6287 70 7770 75 4068 61 5564 40 6776 70 7751 61 7542 71 6564 50 6278 37 83店内展示（display）52 61 6141 44 8676 57 5275 75 6333 69 6052 61 4361 66 4169 43 5165 58 5060 52 5544 45 5852 45 60参考答案（1）箱形图：boxplot(ads) 结果：有异常值。（其它：略）（2）正态性检验 Paper： Hist(X1,6)频数直方图分布的正态性检验： normplot(X1) 均服从正态分布。单因素方差分析ANOVA TableSource SS df MS F ProbFColumns 5866.08333 1955.3611 13.48318.8495e-008Error 20303.2222140145.023Total 26169.3056143P=8.8495e-008FColumns1444.2222 5288.8444 1.9582 0.089763Rows 5866.0833 31955.3611 13.2559 1.5637e-007Interaction1158 1577.2 0.52336 0.92341Error 17701 120147.5083Total26169.3056 143从以上分析结果可知：0.05P1=0.0897630.1，地区对检验水平有一定影响，但不显著。P2=1.5637e-0070.010.1，地区和广告方式对销售量无交互效应。二 Matlab作回归分析回归分析的相关数学理论可以参见概率论与数理统计教程，下面仅以示例说明如何利用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【分析】：首先绘制散点图以直观地选择拟合曲线，这项工作可结合相关专业领域的知识和经验进行，有时可能需要多种尝试。选定目标函数后进行线性化变换，针对变换后的线性目标函数进行回归建模与评价，然后还原为非线性回归方程。【Matlab数据处理】：【Step1】：绘制散点图以直观地选择拟合曲线x=1.5 4.5 7.5 10.5 13.5 16.5 19.5 22.5 25.5;y=7.0 4.8 3.6 3.1 2.7 2.5 2.4 2.3 2.2;plot(x,y,-o)输出图形见图1。图1 销售额与流通费率数据散点图根据图1，初步判断应以幂函数曲线为拟合目标，即选择非线性回归模型，目标函数为：其线性化变换公式为：线性函数为：【Step2】：线性化变换即线性回归建模（若选择为非线性模型）与模型评价% 线性化变换u=log(x);v=log(y);% 构造资本论观测值矩阵mu=ones(length(u),1) u;alpha=0.05;% 线性回归计算b,bint,r,rint,states=regress(v,mu,alpha)输出结果：b = 2.1421; -0.4259表示线性回归模型 中：lna=2.1421，b=-0.4259；即拟合的线性回归模型为；bint = 2.0614 2.2228; -0.4583 -0.3934表示拟合系数lna和b的100(1-alpha)%的置信区间分别为：2.0614 2.2228和-0.4583 -0.3934；r = -0.0235 0.0671 -0.0030 -0.0093 -0.0404 -0.0319 -0.0016 0.0168 0.0257表示模型拟合残差向量；rint = -0.0700 0.0230 0.0202 0.1140 -0.0873 0.0813 -0.0939 0.0754 -0.1154 0.0347 -0.1095 0.0457 -0.0837 0.0805 -0.0621 0.0958 -0.0493 0.1007表示模型拟合残差的100(1-alpha)%的置信区间；states =0.9928 963.5572 0.0000 0.0012表示包含、方差分析的F统计量、方差分析的显著性概率；模型方差的估计值。【注】：严格来讲，模型评价工作应在逆线性化变换后进行；但是，若所建立的线性回归方程不理想，则相应的非线性回归方程必定不理想。【Step3】：拟线性化变换求非线性回归方程（若选择为非线性模型）% 逆线性化变换A=exp(b(1)B=b(2)运行结果为：A = 8.5173；B = -0.4259。即非线性回归方程为：。2.多元线性回归模型（p1）：求得经验回归方程：统计量：总偏差平方和：，其自由度为；回归平方和：，其自由度为；残差平方和：，其自由度为；它们之间有关系：SST=SSR+SSE。多元回归分析的相关数学理论可以参见多元数据分析，下面仅以示例说明如何利用Matlab作多元回归分析。【例2】参见教材P294：10.1 牙膏的销售量。【下面只描述运行程序的过程，应该按照规定格式书写报告】。符号说明：表示价格差；：广告费用；：销售量。【Step1】：绘制散点图以直观地选择拟合曲线clearclcx1=-0.05 0.25 0.60 0 0.25 0.20 0.15 0.05 -0.15 0.15 0.20 0.10 0.40 0.45 0.35 0.30 0.50 0.50 0.40 -0.05 -0.05 -0.10 0.20 0.10 0.50 0.60 -0.05 0 0.05 0.55;x2=5.50 6.75 7.25 5.50 7.00 6.50 6.75 5.25 5.25 6.00 6.50 6.25 7.00 6.90 6.80 6.80 7.10 7.00 6.80 6.50 6.25 6.00 6.50 7.00 6.80 6.80 6.50 5.75 5.80 6.80;y=7.38 8.51 9.52 7.50 9.33 8.28 8.75 7.87 7.10 8.00 7.89 8.15 9.10 8.86 8.90 8.87 9.26 9.00 8.75 7.95 7.65 7.27 8.00 8.50 8.75 9.21 8.27 7.67 7.93 9.26;h1=figure;plot(x1,y,+);h2=figure;plot(x2,y,o);图1 y对x1的散点图图2 y 对 x2的散点图分析图1，可以发现，随着x1的增加，y的值有比较明显的线性增长趋势；分析图2，当x增大时，y有向上弯曲的趋势，可用二次多项式进行逼近；因此可以选择如下方程作为初步的回归模型：【Step2】：模型求解(理论方法：最小二乘法)alpha=0.05;v=ones(length(x1),1) x1 x2 (x2.2);b,bint,r,rint,stats=regress(y,v,alpha)计算结果：b = 17.3244 1.3070 -3.6956 0.3486bint = 5.7282 28.9206 0.6829 1.9311 -7.4989 0.1077 0.0379 0.6594r = -0.0988 -0.0795 -0.1195 -0.0441 0.4660 -0.0133 0.2912 0.2735 -0.2351 0.1031 -0.4033 0.1747 0.0400 -0.1504 0.1284 0.1637 -0.0527 -0.1907 -0.0870 -0.0165 -0.1292 -0.3002 -0.2933 -0.1679 -0.2177 0.1116 0.3035 0.0693 0.2474 0.2270rint = -0.5270 0.3294; -0.5309 0.3718; -0.5106 0.2716; -0.4731 0.3848; 0.0813 0.8507; -0.4609 0.4343; -0.1374 0.7197; -0.0870 0.6340; -0.5960 0.1258; -0.3280 0.5341; -0.8190 0.0125; -0.2618 0.6112; -0.4032 0.4832; -0.5933 0.2925; -0.3207 0.5775; -0.2841 0.6116; -0.4830 0.3776; -0.6248 0.2434; -0.5348 0.3609; -0.4423 0.4092; -0.5609 0.3024; -0.7181 0.1177; -0.7243 0.1377; -0.5548 0.2190; -0.6449 0.2095; -0.2994 0.5226; -0.1037 0.7106; -0.3714 0.5099; -0.1807 0.6755; -0.1890 0.6430stats = 0.9054 82.9409 0.0000 0.0490【Step3】结果分析回归模型为：从结果数据来看，模型整体可用。但也有缺陷，可以改进。【Step4】销售量的预测设需要预测的点为：，则预测值为记，则在处的区间预测为：【模型改进】：当两个因素是不独立时，引入交叉项，新的回归模型为alpha=0.05;v=ones(length(x1),1) x1 x2 (x2.2) (x1.*x2);b,bint,r,rint,stats=regress(y,v,alpha)输出结果：b = 29.1133 11.1342 -7.6080 0.6712 -1.4777bint = 13.7013 44.5252; 1.9778 20.2906; -12.6932 -2.5228; 0.2538 1.0887; -2.8518 -0.1037r = -0.0441; -0.1229; 0.0299; -0.0745; 0.3841; -0.0472; 0.2331; 0.0287; -0.0661; 0.0297; -0.4372; 0.1763; 0.0356; -0.1382; 0.1027; 0.1270; 0.0048; -0.1435; -0.1016; 0.0050; -0.0389; -0.1334; -0.3272; -0.3274; -0.2102; 0.1412; 0.3250; 0.1096; 0.2342; 0.2455rint = -0.4425 0.3542; -0.5408 0.2951; -0.3101 0.3698; -0.4736 0.3247; 0.0245 0.7437; -0.4640 0.3695; -0.1674 0.6337; -0.2369 0.2943; -0.3751 0.2430; -0.3691 0.4284; -0.8118 -0.0627; -0.2306 0.5832; -0.3788 0.4499; -0.5521 0.2757; -0.3172 0.5226; -0.2917 0.5456; -0.3944 0.4039; -0.5490 0.2621; -0.5193 0.3160; -0.3926 0.4026; -0.4360 0.3582; -0.5045 0.2378; -0.7212 0.0667; -0.6326 -0.0221; -0.6085 0.1881; -0.2398 0.5223; -0.0484 0.6984; -0.2988 0.5181; -0.1650 0.6335; -0.1391 0.6302stats = 0.9209 72.7771 0.0000 0.0426结果分析：效果更好。3.逐步回归方法要点：【Step1】 根据问题所属专业领域的理论和经验提出对因变量可能有影响的所有自变量；【Step2】计算每一个自变量对因变量的相关系数，按其绝对值从大到小排序；【Step3】取相关系数绝对值最大的那个自变量建立一元线性回归模型，检验所得回归方程的显著性，若检验表明回归效果则转入【Step4】,若检验表明回归效果不显著则停止建模；【Step4】进行变量的追加、剔除和回归方程的更新操作。Matlab命令：【命令1】：stepwisefit【调用格式】：b,se,pval,inmodel,stats,nextstep,history=stepwisefit(x,y,param1,value1,param2,value2,)【参数说明】：X：p个自变量的n个观测值的矩阵；Y：因变量的n个观测值的矩阵；penter：设置回归方程显著性检验的显著性概率上限，缺省值为0.05；premove：设置回归方程显著性检验的显著性概率下限，缺省值为0.10；display：用来指明是否强制显示建模过程信息，取值为on（显示，缺省设置）和off（不显示）。【例3】某种水泥在凝固时放出的热量(单位：卡/克)Y与水泥中的四种化学成分所占的百分比有关，现测得13组数据如下表：编号X1X2X3X4Y172666078.52129155274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.4作回归分析。【Matlab程序】：clearclcload haldb,se,pval,inmodel,stats,nextstep,history=stepwisefit(ingredients,heat,penter,0.10,display,off);% 自变量的筛选和模型参数估计信息inmodel,b0=ercept,b% 回归方程显著性整体检验信息Allp=stats.pval,rmse=stats.rmse% 回归方程显著性分别检验信息P=stats.PVAL输出结果：inmodel = 1 1 0 0；b0 = 52.5773；b = 1.4683 0.6623 0.2500 -0.2365；Allp = 4.4066e-009；rmse = 2.4063；P = 0.0000 0.0000 0.2089 0.2054。结果分析：最优回归方程为，回归方程显著性整体检验和分别检验均为高度显著，模型标准误差估计为2.4063。【命令2】：stepwise【调用格式】：stepwise(x,y,inmodel,penter,premove)【说明】：创建多元线性回归分析的逐步回归法建模的交互式图形环境。【图形界面说明】：窗口1：Coefficients with error Bars绘出各个解释变量回归系数的估计，圆点表示点估计值，横线表示置信区间（有色线段表示90%置信区间，黑色线段表示95%置信区间）。窗口的右侧给出回归系数的点估计值（Coeff）、显著性检验的t统计量的值（t-test）和显著性概率p值（p-val）.窗口2：Model History该窗口绘出的圆点表示历次建模的模型标准差的估计。两个窗口中间输出的是当前模型的有关信息，包括：Intercept：模型截距（常数项）的估计；RMSE：模型标准差的估计；R-square：可决系数；Adj-R-sq：校正可决系数；F：模型整体性检验的 F 统计量的值；p：模型整体性检验的显著性概率。窗口I右侧的三个按钮：Next Step：在回归方程中按相关系数绝对值大小逐次引入解释变量，如无解释变量可引入时，按钮不可用；All Steps：直接给出“只进不出”方式建模的最终结果（注意，此时的回归方程未必是最优回归方程）；Export：选择向Workspace传输的计算结果（有关变量名可由用户自定义）stepwise(ingredients,heat,1 1 1 1,0.05,0.10);三 matlab作相关分析一、相关系数 要初步研究变量之间的随机性关系，我们就要清楚，研究的对象是二元或多元的随机向量，利用的是成对观测数据。首先绘制一张散点图，直观上大致判断两两变量之间是否存在某种关系。MATLAB命令（散点图）：gscatter - 两个变量的散点图. 用法：gscatter(x,y)lsline - 在散点图上增加最小二乘拟合线. 用法：lslinegplotmatrix 矩阵散点图。用法：gplotmatrix(x,y)。 其中x，y都是矩阵，行数相同。例如：x=normrnd(0,1,100,3);y=normrnd(1,2,100,2);gplotmatrix(x,y) 如果认为两个变量之间存在着某种直线关系，我们可以用相关系数来刻画这种关系。首先引入如下样本相关系数的概念：对二元总体(X，Y)的样本，定义样本相关系数为其中分别为X和Y的样本方差，叫X与Y之间的样本协方差。这是一个重要统计量，与总体相关系数相对应。那么，怎样充分发挥这个统计量的作用呢？下面我们讲讲如何利用它对总体相关系数作假设检验和区间估计。 原假设为对立假设为 在原假设成立的情况下，可以证明下面的统计量服从自由度为n-2的t分布： 所以给定检验水平，可得原假设的否定域。 MATLAB命令：corrcoef用法：r,p,rlo,rup=corrcoef(x)其中：x 矩阵；r：相关矩阵；p：p-值；rlo：置信下限；rup：置信上限二、偏相关分析基本描述：控制其它变量的情况下研究两个变量之间的线性关系.因变量 自变量1 自变量2原假设:两个变量之间的偏相关系数为0MATLAB命令：partialcorr用法：RHO,PVAL = PARTIALCORR(X,Z) X是由多个变量的样本值构成的矩阵，Z是由控制变量构成的矩阵，RHO是偏相关系数矩阵，PVAL是对应的p-值。例：财政收入预测问题：财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试对各因素之间的关系进行初步分析。年份国民收入（亿元）工业总产值(亿元)农业总产值（亿元）总人口（万人）就业人口（万人）固定资产投资（亿元）财政收入(亿元)19525983494615748220729441841953586455475587962136489216195470752049160266218329724819557375585296146522328982541956825715556628282301815026819578377985756465323711139286195810281235598659942660025635719591114168150967207261733384441960107918704446620725880380506196175711564346585925590138271196267796446167295251106623019637791046514691722664085266196494312505847049927736129323196511521581632725382867