matlab与统计回归分析

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

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

3、1 85 68 92 90;73 79 60 75 81;P=anova1(Score)输出结果：方差分析表和箱形图ANOVA TableSource SS df MS FProb>FColumns 604.93332302.46674.25610.040088Error852.81271.0667Total 1457.733314seulaV6012Column Number3908075p 值小于 0.05，拒绝零假设，认为3 种教学方法存在显著差异。例 2（双因素方差分析）为了考察4 种不同燃料与3 种不同型号的推进器对火箭射程（单位：海里）的影响，做了12 次试验，得数据如表2

4、所示。表 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

5、的所有样本（X 中的所有列样本）取自相同的总体；H0B ：因素B 的所有样本（X 中的所有行样本）取自相同的总体。当 reps>1 时， 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.7p=anova2(disp1,1)输出结果：方差分析表ANOVA TableSourceSSdfMSProb>FColumns157.5952.530.430

6、59 0.73875Rows223.8467111.92330.917430.44912Error731.986121.9967Total1113.4167 11p 值均大于0.05 ， 所以可以接受零假设H0A 和 H0B， 认为燃料和推进器对火箭的射程没有显著影响。例3（双因素方差分析） 设火箭的射程在其它条件基本相同时与燃料种类及推进器型号有关。现在考虑4 种不同的燃料及3 种不同型号的推进器，对于每种搭配个发射了火箭两次，得数据见表3。问各自变量和自变量的交互效应是否对火箭的射程有显著影响？表 3 燃料-推进器-射程数据表推进器1推进器2推进器3燃料158.256.265.352.64

7、1.260.8燃料42.850.548.4燃料360.170.939.258.373.240.7燃料475.858.248.771.551.041.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.4p=anova2(disp2,2)输出结果：方差分析表ANOVA TableSourceSSdf MSProb>FColumns370.9808 2

8、185.49049.39390.003506Rows261.67587.2254.41740.025969Interaction1768.6925294.782114.9288 6.1511e-005Error236.951219.7458Total2638.298323显著。方差分析上机练习为研究广告的效果，考察 4 种广告方式：当地报纸 （ paper） 、 当地广播 （ radio） 、店内销售员（people ）和店内展示（display）的效果。共设有144 个销售点，每种广告随机抽取36 个销售点记录销售额，分布在 6 个地区的144 个销售点的销售情况生成的数据集ADS 见下表。

9、数据集ADS 中有 3 个变量： AD 表示广告AREA 表示地区、SALES 表示销售额（单位： 千元） 。 请完成以下练习：（1） 概括下列数据：用箱形图、条形图直观地呈现四种广告方式下销售量的 分布情况；计算四种广告方式下销售量的均值、方差、标准差、最大和最小值；(2) 进行单因素方差分析：检验四种广告方式下销售量数据是否服从正态分布，方差是否相等；检验四种广告方式下的销售量是否有显著差异(0.01 ) ；若四种广告方式下的销售量有显著差异，指出哪些类型的广告效果有显著的不同？(3) 在设计广告效果的试验时，虽然地区差异对销售量的影响并不是我们感兴趣的，但希望排除这一因素的影响。数据集A

10、DS 记录了各个销售点所在的地区 AREA。 试用双因素方差分析方法分析销售数据，并指出广告方式和地区对销售量是否有显著影响(0.01,0.1 )？广告方式(AD)与地区 (AREA)之间有无交互效应？表 ADS 数据集中的数据广告方式(变量:AD)销售额(单位：千元)(变量SALES )地区1地区2地区 3地区4地区 5地区 6当地报纸75 57 7677 7576 8194 54 7087 6579 62(paper)68 75 83726388 56 86657566 6670 8684 7780 6276627870当地广播69 5190 7733 79100 6168 6376 73

11、(radio)100607368837454 78 7983 7468 7570 53 7379 6681 5769656565店内销售员63 67 8580 8770 7564 40 6751 6164 50( people )58 82 78624076 70 77756287 707768 615542 716578 3783店内展示52 61 6176 5733 6961 66 4165 5844 45( display )41 44 86526069 43 51505875 7552 6160 5252 4563435560参考答案1)箱形图：boxplot(ads) 结果：有异常值

12、。(其它：略)2)正态性检验Paper： Hist(X1,6)50556065707580859095频数直方图分布的正态性检验：normplot(X1)0.990.980.950.90Normal Probability Plot505 .7 .5 .2000 ytilbaborP0 5 21 .1 .0 .0 .00 0 00箱形图为：55606558805990均服从正态分布。ANOVA TableSourceSSdfMSProb>FColumns5866.0833 31955.361113.48318.8495e-008Error20303.2222140 145.023Tota

13、l 26169.3056143P=8.8495e-008<0.01 ，四种广告方式下的销售量是否有显著差异。对应的selaS1 paper2 radio34peopledisplay00090087006530当地媒体和电台广播与点内展示有显著不同。3)ANOVA TableError 17701120 147.5083SourceSSdfMSFProb>FColumns1444.22225 288.84441.95820.089763Rows5866.083331955.3611 13.25591.5637e-007Interaction 115815 77.20.523360.

14、92341Total 26169.3056143从以上分析结果可知：0.05<P1=0.089763<0.1 ，地区对检验水平有一定影响，但不显著。P2=1.5637e-007<0.01<0.1 ，无论哪种检验水平，广告方式对销售量都有显著影响；P3=0.92341>0.1 ，地区和广告方式对销售量无交互效应。二 Matlab 作回归分析回归分析的相关数学理论可以参见概率论与数理统计教程，下面仅以示例说明如何利用matlab 处理回归分析。1. 一元线性回归分析【例1】 为了了解百货商店销售额x 与流通费率（反映商业活动的一个质量指标，指每元商品流转额所分摊的流通

15、费用）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 】 ：绘制散点图以直观地选择拟合曲线

16、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。105010152025图 1 销售额与流通费率数据散点图根据图 1，初步判断应以幂函数曲线为拟合目标，即选择非线性回归模型，目标函数为：y axb (b 0)其线性化变换公式为：v ln y,u lnx线性函数为：v ln a bu【 Step2】 ：线性化变换即线性回归建模(若选择为非线性模型)与模型评价% 线性化变换u=log(x)'v=log(y)

17、9;% 构造资本论观测值矩阵mu=ones(length(u),1) u;alpha=0.05;% 线性回归计算b,bint,r,rint,states=regress(v,mu,alpha)输出结果：b = 2.1421; -0.4259表示线性回归模型v ln a bu 中： lna=2.1421 ， b=-0.4259 ；即拟合的线性回归模型为y 2.1421 0.4259x；bint = 2.0614 2.2228; -0.4583 -0.3934表示拟合系数lna 和 b 的 100(1-alpha)% 的置信区间分别为：2.0614 2.2228 和 -0.4583 -0.3934

18、；r = -0.0235 0.0671 -0.0030 -0.0093 -0.0404 -0.0319 -0.00160.0168 0.0257表示模型拟合残差向量；rint = -0.0700 0.02300.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.0012SSR表示包含R2 SSR 0.992

19、8、SSR963.5572、SSE/(n 2)SST方差分析的F 统计量 F SSR/ fRSSE/ fE方差分析的显著性概率p P(F(1,n 2) F) 0；模型方差的估计值?2 SSE 0.0012。n2【注】 ：严格来讲，模型评价工作应在逆线性化变换后进行；但是，若所建立的线性回归方程不理想，则相应的非线性回归方程必定不理想。【 Step3】 ：拟线性化变换求非线性回归方程(若选择为非线性模型)% 逆线性化变换A=exp(b(1)B=b(2)运行结果为：A = 8.5173 ； B = -0.4259 。即非线性回归方程为：0.4259y 8.5173x。2. 多元线性回归模型(p&g

20、t;1) ：y01x12x2pxp,N(0, 2)求得经验回归方程：y?0?1x1?2x2?pxp统计量：n总偏差平方和：SST(yi y)2，其自由度为fT n 1；i1n回归平方和：SSR(y?i y)2，其自由度为fR p；i1n残差平方和：SSE(yi y?i )2 ，其自由度为fE n p 1；i1它们之间有关系：SST=SSR+SS。 E多元回归分析的相关数学理论可以参见多元数据分析，下面仅以示例说明如何利用 Matlab 作多元回归分析。【例2】参见教材P294： 10.1 牙膏的销售量。【下面只描述运行程序的过程，应该按照规定格式书写报告】。符号说明：x1 ：表示价格差；x2：

21、广告费用；y ：销售量。【 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.450.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 00.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.906.80 6.80 7.10 7.00 6.80 6.50 6

22、.25 6.00 6.50 7.00 6.80 6.80 6.50 5.755.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.868.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.939.26;h1=figure;plot(x1,y,'+' );h2=figure;plot(x2,y,'o' );109.598.587.57-0.2-0.30.

23、40.50.6109.598.587.5755.566.577.51 y 对 x1 的散点图图 2 y 对 x2 的散点图分析图1，可以发现，随着x1 的增加， y 的值有比较明显的线性增长趋势；分析图2，当x 增大时， y 有向上弯曲的趋势，可用二次多项式进行逼近；因此可以选择如下方程作为初步的回归模型：y 01x12x23x22, N(u, 2)【 Step2】 ：模型求解( 理论方法：最小二乘法)alpha=0.05;v=ones(length(x1),1) x1' x2' (x2.2)'b,bint,r,rint,stats=regress(y',v,a

24、lpha)计算结果：b = 17.3244 1.3070 -3.6956 0.3486bint = 5.7282 28.92060.6829 1.9311-7.4989 0.10770.0379 0.6594r = -0.0988-0.07950.2912 0.2735-0.23510.10310.1284 0.1637-0.0527-0.1907-0.2933 -0.1679-0.1195-0.0441-0.40330.1747-0.0870-0.01650.4660-0.01330.0400-0.1504-0.1292-0.3002-0.2177 0.1116 0.3035 0.0693

25、0.2474 0.2270rint =-0.52700.3294;-0.53090.3718;-0.51060.2716;-0.4731 0.3848-0.1374 0.7197-0.3280 0.5341-0.4032 0.4832-0.2841 0.6116-0.5348 0.3609-0.7181 0.1177-0.6449 0.2095-0.3714 0.50990.0813 0.8507;-0.0870 0.6340;-0.8190 0.0125;-0.5933 0.2925;-0.4830 0.3776;-0.4423 0.4092;-0.7243 0.1377;-0.2994 0

26、.5226;-0.1807 0.6755;-0.4609 0.4343;-0.5960 0.1258;-0.2618 0.6112;-0.3207 0.5775;-0.6248 0.2434;-0.5609 0.3024;-0.5548 0.2190;-0.1037 0.7106;-0.1890 0.6430stats = 0.9054 82.9409 0.0000 0.0490Step3】结果分析回归模型为：2y? 17.3244 1.3070x1 3.6959x2 0.3486x22从结果数据来看，模型整体可用。但也有缺陷，可以改进。Step4】销售量的预测设需要预测的点为：x0 (x01

27、 ,x02 ,x0 p )' ，则预测值为y?0?0?1x01?2x02?px0SSE , n p1pp11(x0ixi )(x0jxj )cij , xin i1 j11xki ,i nk11,2, , pT1(cij)p p (XTX) 1，则在x0处的区间预测为：(y?0 t(1 )(n p 1)?* , y?0 t(1 )(n p 1)?* )：当两个因素是不独立时，引入交叉项x1x2，新的回归模型为y 01x12x23x224x1x2alpha=0.05;v=ones(length(x1),1) x1' x2' (x2.2)' (x1.*x2)'

28、;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;

29、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.230

30、6 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.6302sta

31、ts = 0.9209 72.7771 0.0000 0.0426结果分析：效果更好。3. 逐步回归方法要点：【 Step1 】 根据问题所属专业领域的理论和经验提出对因变量可能有影响的所有自变量；【 Step2】计算每一个自变量对因变量的相关系数，按其绝对值从大到小排序；【 Step3】取相关系数绝对值最大的那个自变量建立一元线性回归模型，检验所得回归方程的显著性，若检验表明回归效果则转入【 Step4】 , 若检验表明回归效果不显著则停止建模；【 Step4】进行变量的追加、剔除和回归方程的更新操作。Matlab 命令：【命令1】 ： stepwisefit【调用格式】：b,se,pval

32、,inmodel,stats,nextstep,history=stepwisefit（x,y, param1 ,value1, param2 ,value2, ）【参数说明】：X： p个自变量的n个观测值的n p 矩阵；Y：因变量的n个观测值的n 1 矩阵； penter ：设置回归方程显著性检验的显著性概率上限，缺省值为 0.05； premove：设置回归方程显著性检验的显著性概率下限，缺省值为0.10； display ：用来指明是否强制显示建模过程信息，取值为on（显示，缺省设置）和off （不显示）。【例3】 某种水泥在凝固时放出的热量（ 单位：卡 / 克 ）Y 与水泥中的四种化学

33、成分所占的百分比有关，现测得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,

34、he at, '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.4683 0.6623 0.2500 -0.2365；Allp = 4.4066e-009；rmse = 2.4063 ；P = 0.0000 0.0

35、000 0.2089 0.2054。结果分析：最优回归方程为y? 52.5773 1.4683x1 0.6623x2，回归方程显著性整体检验和分别检验均为高度显著，模型标准误差估计为2.4063。【命令2】 ： stepwise【调用格式】：stepwise（x,y,inmodel,penter,premove）【说明】 ：创建多元线性回归分析的逐步回归法建模的交互式图形环境。【图形界面说明】：窗口 1： Coefficients with error Bars绘出各个解释变量回归系数的估计，圆点表示点估计值，横线表示置信区间（有色线段表示90%置信区间，黑色线段表示95%置信区间）。窗口的右

36、侧给出回归系数的点估计值（ Coeff ） 、 显著性检验的t 统计量的值（ t-test ）和显著性概率p值（p-val ） .窗口2： Model History该窗口绘出的圆点表示历次建模的模型标准差的估计。两个窗口中间输出的是当前模型的有关信息，包括：Intercept ：模型截距（常数项）的估计；RMS：模型标准差E的估计；R-square：可决系数；Adj-R-sq ：校正可决系数；F：模型整体性检验的F 统计量的值；p：模型整体性检验的显著性概率。窗口 I 右侧的三个按钮：Next Step ：在回归方程中按相关系数绝对值大小逐次引入解释变量，如无解释变量可引入时，按钮不可用；A

37、ll Steps ：直接给出“只进不出”方式建模的最终结果（注意，此时的回归方程未必是最优回归方程）；Export ：选择向Workspace传输的计算结果（有关变量名可由用户自定义）stepwise（ingredients,heat,1 1 1 1,0.05,0.10）;三 matlab 作相关分析一、相关系数要初步研究变量之间的随机性关系，我们就要清楚，研究的对象是二元或多元的随机向量，利用的是成对观测数据。首先绘制一张散点图 ，直观上大致判断两两变量之间是否存在某种关系。MATLAB 命令(散点图)：gscatter - 两个变量的散点图.用法： gscatter(x,y)lsline

38、- 在散点图上增加最小二乘拟合线.用法：lslinegplotmatrix 矩阵散点图。用法： gplotmatrix(x,y)x， y 都是矩阵，行数相同。例如：x=normrnd(0,1,100,3);y=normrnd(1,2,100,2);gplotmatrix(x,y)如果认为两个变量之间存在着某种直线关系，我们可以用相关系数来刻画这种关系。首 先 引 入 如 下 样 本 相 关 系 数 的 概 念 ： 对 二 元 总 体 (X， Y) 的 样 本 (xi, yi ), i 1,2, ,n ，定义样本相关系数为SX 0 SY0SX2Y1n其中SX 0ni1221n2(Xi X ) ,

39、SY0(Yi Y) 分别为 X和 Y的样本方差，ni1ni1SX2Y 1(Xi X)(Yi Y) 叫 X 与 Y之间的样本协方差。这是一个 重要统计量，与总体相关系数(X,Y) 相对应。那么， 怎样充分发挥这个统计量的作用呢？下面我们讲讲如何利用它对总体相关系数: (X ,Y) 作假设检验和区间估计。原假设为 H 0 :0; 对立假设为H 1 :0.在原假设成立的情况下，可以证明下面的统计量服从自由度为n-2 的 t分布：1r所以给定检验水平，可得原假设的否定域t t /2(n 2)MATLAB 命令： corrcoef用法： r,p,rlo,rup=corrcoef(x)其中： x 矩阵；r

40、：相关矩阵；p： p-值；rlo：置信下限；rup：置信上限二、偏相关分析基本描述：控制其它变量的情况下研究两个变量之间的线性关系.因变量自变量 1 自变量 2原假设:两个变量之间的偏相关系数为0MATLAB 命令： partialcorr用法： RHO,PVAL = PARTIALCORR(X,Z)X 是由多个变量的样本值构成的矩阵，Z 是由控制变量构成的矩阵，RHO 是偏相关系数矩阵，PVAL 是对应的p-值。例： 财政收入预测问题：财政收入与国民收入、工业总产值、农业总产值、总人口、 就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试对各因素之间的关系进行初步分析。年份 国民收工业总农业总总人口就业人固定资产财政入（亿元）产值（亿元）产值 （亿 元）（万 人）口（万人）投资（亿元）收入（亿元）19525983494615748220729441841953586455475587962136489216195470752049160266218329724819557375585296146522328982541956825715556628282301815026819578377985756465323711139286195810281235598