数学建模之统计回归模型

1、数学建模大作业统计回归模型摘要某公司想用全行业的销售额作为自变量来预测公司的销售额，题目给出了19771981此公司的销售额和行业销售额的分季度数据表格。通过对所给数据的简单分析，我们可以看出：此公司的销售额有随着行业销售额的增加而增加的趋势，为了更加精确的分析题目所给的数据，得出科学的结论，从而达到合理预测的目的。我们使用时间序列分析法，参照课本统计回归模型例4，做出了如下的统计回归模型。在问题一中，我们使用MATLB数学软件，画出了数据的散点图，通过观察散点图，发现公司的销售额和行业销售额之间有很强的线性关系，于是我们用线性回归模型去拟合，发现有很好的拟合性。但是这种情况下，并没有考虑到数

2、据的自相关性，所以我们做了下面几个问题的分析来对这个数学模型进行优化。在问题二中，通过建立了公司销售额对全行业销售额的回归模型，并使用DW检测诊断随机误差项的自相关性。通过计算和查DW表比较后发现随即误差存在正自相关，也就是说前面的模型有一定的局限性，预测结果存在一定的偏差，还有需要改进的地方。在问题三中，因为在问题二中得出随即误差存在正自相关，为了消除随机误差的自相关性，我们建立了一个加入自相关后的回归模型。并对其作出了分析和验证，我们发现加入自相关后的回归模型更加合理。通过使用我们建立的模型对公司的销售额进行预测，发现和实际的销售额很接近，也就是说模型效果还不错。关键词：销售额、回归模型、

3、自相关性一、问题提出某公司想用全行业的销售额作为自变量来预测公司的销售额，下表给出了1977-1981年公司销售额和行业销售额的分季度数据（单位：百万元）.（1） 画出数据的散点图，观察用线性回归模型拟合是否合适。（2） 监理公司销售额对全行业销售额的回归模型，并用DW检验诊断随机误差项的自相关性。（3） 建立消除了随机误差项自相关性后的回归模型。年季公司销售额行业销售额年季公司销售额行业销售额19771978197912341234121234567891020.9621.4021.9621.5222.3922.7623.4823.6624.1024.01127.3130.0132.7129

4、.4135.0137.1141.2142.8145.5145.319791980198134123412341112131415161718192024.5424.3025.0025.6426.3626.9827.5227.7828.2428.78148.3146.6150.2153.1157.3160.7164.2165.6168.7171.7二、基本假设假设一：模型中（对时间t）相互独立。三、符号说明公司销售额：（百万）行业销售额：（百万）概念介绍：1.自相关：自相关（auto correlation），又称序列相关（serial correlation）是指总体回归模型的随机误差项之间存

5、在的相关关系。即不同观测点上的误差项彼此相关。2置信区间：如果P（）=1-，=0.1或0.05，则称区间a,b为的置信度为1-的置信区间。3.时间序列：时间序列法是一种定量预测方法，亦称简单外延方法。时间序列即按时间的推移或排布会对规律的变化有所影响。 四、问题分析问题一：表中的数据是以时间为顺序的。由于前期的销售额对后期的投资一般有明显的影响，从而对后期的后期的销售额造成影响。因此在此模型中应考虑到存在自相关，我们可以先建立基本的回归模型，然后再进行自相关性诊断，并建立新的回归模型。问题二：在问题一之后，就可以接着求出问题二，然后利用DW检验诊断随机误差项的自相关性。问题三：进行了自相关诊断

6、后，将自相关加入模型中，建立消除了随机误差项自相关性的回归模型。五、模型的建立与求解5.1 问题一5.1.1 问题一的分析表中数据是以时间为序的，建立基本的回归模型。5.1.2 问题一模型的建立基本回归模型：设该公司第时间的公司销售额为，行业销售额为 。为了大致分析和的关系，首先利用表中的数据作出对关系作出散点图，如下（见图中的“+”）：做散点图：可以看出，随着行业销售额的增加，公司销售额增大，而且两者有很强的线性关系，图中的直线说明两者呈线性模型，因此本题用线性回归模型拟合非常合适。5.2 问题二5.2.1 问题二的分析从问题一中的图形可以看出，随着行业销售额的增加，公司销售额增大，而且两者

7、有很强的线性关系，图中的直线说明两者呈线性模型，因此可建立一元线性回归模型。5.2.2 问题二模型的建立由题意建立一元线性回归模型 （1）模型（1）中除了行业销售额和公司销售额的影响外，影响的其他因素都包含在随机误差内，这里假设（对t相互独立）且服从均值为零的正态分布N(0, )。5.2.3 问题二模型的求解根据表中的数据。对模型（1）直接利用MATLAB统计工具箱求解（具体算法见附录），得到的回归系数估计值及置信区间（置信水平=0.05）、检验统计量,的结果见下表：参数参数估计值参数置信区间-1.4548 【-1.9047 -1.0048】0.1763【0.1732 0.1793】R=1.0

8、e+004 *0.0001F=1.0e+004 *1.4888P=1.0e+004 *0.0000将参数估计值代入（1）得到： (2)用MATLAB中rstool命令得到的交互式画面见图 （1） ，由此可以得出不同水平下的预测值及其置信区间。通过左下方的Export下拉式菜单。可以输出模型的统计结果。图1自相关性诊断与处理方法 从表面上来看得到的基本模型（2）拟合度（R）非常之高，接近你100%，应该很满意了，但是，这个模型并没有考虑到我们的数据是一个时间序列（将原表中的数据打乱不影响模型（2）的结果）。实际上对于时间序列数据做回归分析时，模型的随机误差有可能存在相关性，违背模型关于（对时间t

9、）相互独立的基本假设，其他相关因素对公司销售额的影响肯能也有时间上的延续，包含在随机误差中，即随机误差会出现自相关性。残差可以作为随机误差的估计值，画出的散点图，能够从直观上判断的自相关性。模型（2）的残差可在计算过程中得到表1，以及数据的图见图 2t12345e-0.0282-0.06420.01980.16160.0443t678910e0.04410.0412-0.0608-0.0968-0.1516t1112131415e-0.1505-0.0555-0.02550.10330.0828t1617181920e0.10340.02630.0395-0.047-0.0359表 1图 2为

10、了对的字相关性做定量的诊断，并在确诊后得到新的结果，我们考虑如下模型 ， （3）其中是自相关系数，|1,相互独立且服从均值为0的正态分布。若=0，则退化为普通的回归模型；若>0，则随机误差存在正的自相关；若<0，则随机误差存在负的自相关。利用D-W检验诊断自相关现象如下：利用MATLAB算出:=0.0980 =0.1326DW=0.7388 =0.6306(具体程序见附录)因为DW2（1-），所以 0DW4,若的估计值在0附近，则DW的值在2附近，的自相关行很弱，若在正负1附近，则DW接近0或4，的自相关性很强。5.2.4 问题二结果的分析及验证要根据DW的具体数值确定是否存在自相

11、关，查D-W分布表，可以得到检验的临界值和，然后根据区间来确定。利用表1给出的残差,根据以上式子可得出DW=0.7388，对于显著性水平=0.05，n=20,k=2,查D-W分布表，得到检验的临界值=1.2和=1.4 .现在DW<,因此可以认为随即误差存在正自相关，而且可得出=0.6306。 5.3 问题三5.3.1 问题三的分析题目要求建立消除了随机误差项目自相关性后的回归模型，即是加入了自相关后的回归模型，下面我们将自相关性加入问题中。5.3.2 问题三模型的求解加入自相关后的回归模型 = 做变换 ， （4）则模型（3）转化为 ， （5）其中相互独立且服从均值为零的正态分布，所以（5

12、）是普通回归模型。以的估计值带入（3）和（4）做变换，利用变换后的数据 ，估计模型（5）的参数，得到的表见表2，还可以得出剩余标准差rmse=0.08828.最后将模型（5）的变量还原为原始变量。得到的结果如下 （6）表 2参数参数估计值参数置信区间-0.3951 -0.7481 -0.04220.17380.1675 0.1800R=1.0e+003 *0.0010F=1.0e+003 *3.4621P=1.0e+003 *0.00005.3.4 问题三结果的分析及验证当然应该对模型（6）也作一次自相关检验，即诊断随机误差是否还存在自相关，从模型（6）的残差可以计算出DW=1.65，对于显著

13、水平=0.05，n=19,k=2,查D-W分布表，得到检验的临界值=1.2和=1.40 .现在,可以认为随机误差不存在自相关。一次经变换得到的回归模型（6）是适用的。结果及其预测从机理上看，对于带滞后性的经济规律作用下的时间序列数据，加入自相关的模型（6）更为合理，我们将模型（1）与模型（6）的计算值与实际数据的比较，以及两个模型的残差，表示在表 3 中 表 3ty(实际数据)yy(模型1)yyy（模型2）eee1221.421.46421.464-0.06424.00E-06321.9621.9421.9150.019790.02521421.5221.35821.3990.16158-0.

14、04026522.3922.34622.4560.0443-0.11047622.7622.71622.7560.04407-0.04008723.4823.43923.4720.04124-0.033823.6623.72123.755-0.06084-0.03367924.124.19724.162-0.096850.0349341024.0124.16224.109-0.151590.052891124.5424.6924.595-0.150490.0952241224.324.35624.27-0.055520.085056132525.02524.988-0.025460.03766

15、1425.6425.53725.5170.103270.019971526.3626.27726.3320.08281-0.055271626.9826.87726.9170.10339-0.040491727.5227.49427.5440.02634-0.050071827.7827.7427.7440.03952-0.003491928.2428.28728.293-0.04701-0.006262028.7828.81628.765-0.035910.050926六、模型的评价与推广模型的评价与推广：此模型从最初的线性回归模型到DW检验诊断随机误差的自相关性，再到最后的消除了随机误差项

16、自相关性后的回归模型，模型逐步得到了优化。最后的结果预测可以看出，我们建立的这个模型的可靠性是非常高的。预测公司的销售额可以为公司的制定相应的生产计划或者购货数量提供依据，鉴于该模型的可靠性非常稳定，我们可以把此模型推广到公司其他产品的销售额或者某产品的市场销售额的预测。同时，通过这学期的建模课程的学习嘛，我们发现团队精神是数学建模是否取得好成绩的最重要的因素，一队三个人要相互支持，相互鼓励。切勿自己只管自己的一部分（数学好的只管建模，计算机好的只管编程，写作好的只管论文写作），很多时候，一个人的思考是不全面的，只有大家一起讨论才有可能把问题搞清楚，因此无论做任何板块，三个人要一起齐心才行，只

17、靠一个人的力量，要写出一篇高水平的文章几乎是不可能的。其实建模的过程就是大家互相鼓励，共同勉励的一个阶段，我们组从最开始的模拟训练时就十分注重团队的分工协作，在作业中都会总结教训，改进方法。另外我们还根据每个人的特长来进行分工，做到发挥优势，长短互补的效果。七、参考文献数学模型（第三版） 姜启源 谢金星 叶俊 高等教育出版社，2003年8月 数学建模上课课件：统计回归模型 应用回归分析 何晓群，刘文清 中国人民大学出版社，2001网上资源 DW表八、附录MATLAB运行程序附录1： >> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8

18、,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'yy=-1.4548+0.1763*xyy =20.9882 21.4642 21.9402 21.3584 22.3457 22.7159 23.4388 23.7208 24.1969 24.1616

19、24.6905 24.3555 25.0255 25.5367 26.2772 26.8766 27.4937 27.7405 28.2870 28.8159附录2：>> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'yy=-1.4548+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.

20、54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'e=y-yy;e1=e(2:20,:);e2=e(1:19,:);y0=sum(e1-e2).2);y1=sum(e1.2);DW=y0/y1;p=1-0.5*DWp = 0.6306>> >> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'yy=-1.454

21、8+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'e=y-yy;e1=e(2:20,:);e2=e(1:19,:);y0=sum(e1-e2).2)y0 = 0.0980>> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164

22、.2,165.6,168.7,171.7'yy=-1.4548+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'e=y-yy;e1=e(2:20,:);e2=e(1:19,:);y0=sum(e1-e2).2);y1=sum(e1.2)y1 = 0.1326>> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.

23、5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'yy=-1.4548+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'e=y-yy;e1=e(2:20,:);e2=e(1:19,:);y0=sum(e1-e2).2);y1=sum(e1.2);DW=y0/y1DW = 0.7388>&

24、gt; >> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'yy=-1.4548+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'y2=y(2:20,:)y3=y(1:19,

25、:)x2=x(2:20,:);x3=x(1:19,:);y4=y2-y3*p; x4=x2-x3*py2 = 21.4000 21.9600 21.5200 22.3900 22.7600 23.4800 23.6600 24.1000 24.0100 24.5400 24.3000 25.0000 25.6400 26.3600 26.9800 27.5200 27.7800 28.2400 28.7800y3 = 20.9600 21.4000 21.9600 21.5200 22.3900 22.7600 23.4800 23.6600 24.1000 24.0100 24.5400 2

26、4.3000 25.0000 25.6400 26.3600 26.9800 27.5200 27.7800 28.2400x4 = 49.7254 50.7227 45.7201 53.4011 51.9698 54.7455 53.7601 55.4511 53.5485 56.6747 52.8829 57.8810 58.3847 60.7560 61.5075 62.8635 62.0564 64.2736 65.3187>> x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,1

27、50.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'yy=-1.4548+0.1763*x ;y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'y2=y(2:20,:)y3=y(1:19,:)x2=x(2:20,:);x3=x(1:19,:);y4=y2-y3*p; x4=x2-x3*p;b1,bint1,r1,rint1,stats1=regress(y4

28、,ones(19,1),x4)y2 = 21.4000 21.9600 21.5200 22.3900 22.7600 23.4800 23.6600 24.1000 24.0100 24.5400 24.3000 25.0000 25.6400 26.3600 26.9800 27.5200 27.7800 28.2400 28.7800y3 = 20.9600 21.4000 21.9600 21.5200 22.3900 22.7600 23.4800 23.6600 24.1000 24.0100 24.5400 24.3000 25.0000 25.6400 26.3600 26.9

29、800 27.5200 27.7800 28.2400b1 = -0.3951 0.1738bint1 = -0.7481 -0.0422 0.1675 0.1800r1 = -0.0627 0.0466 0.1227 -0.0645 0.0056 0.0099 -0.0929 -0.0602 -0.0971 -0.0535 0.0311 0.0140 0.1250 0.0294 0.0648 -0.0218 0.0379 -0.0513 0.0170rint1 = -0.1941 0.0688 -0.0886 0.1817 0.0163 0.2291 -0.2012 0.0721 -0.13

30、37 0.1448 -0.1317 0.1516 -0.2252 0.0395 -0.1986 0.0782 -0.2284 0.0343 -0.1928 0.0858 -0.1083 0.1705 -0.1277 0.1556 -0.0003 0.2503 -0.1091 0.1679 -0.0693 0.1988 -0.1573 0.1137 -0.0981 0.1739 -0.1814 0.0788 -0.1128 0.1468stats1 = 1.0e+003 * 0.0010 3.4621 0 0.0000>> >> >> y=20.96,21.4

31、,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'yy=-1.4548+0.1763*x;x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'x2=x(2:20,:);x3=x(1:19,:);y3=y(1:19,:);yyy=-0.39

32、51+0.6306*y3+0.1738*x2-0.1096*x3;yyyy=yy(2:20,:);e1=yyyy-yyye1 = 0.0000 0.0252 -0.0403 -0.1105 -0.0401 -0.0330 -0.0337 0.0349 0.0529 0.0952 0.0851 0.0377 0.0200 -0.0553 -0.0405 -0.0501 -0.0035 -0.0063 0.0509>> y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.3

33、6,26.98,27.52,27.78,28.24,28.78'yy=-1.4548+0.1763*x;x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2,165.6,168.7,171.7'x2=x(2:20,:);x3=x(1:19,:);y3=y(1:19,:);yyy=-0.3951+0.6306*y3+0.1738*x2-0.1096*x3;yyyy=yy(2:20,:);e1=yyyy-yyy;e=y-yy;ee=e

34、(2:20,:)ee = -0.0642 0.0198 0.1616 0.0443 0.0441 0.0412 -0.0608 -0.0968 -0.1516 -0.1505 -0.0555 -0.0255 0.1033 0.0828 0.1034 0.0263 0.0395 -0.0470 -0.0359>> y=20.96,21.4,21.96,21.52,22.39,22.76,23.48,23.66,24.1,24.01,24.54,24.3,25,25.64,26.36,26.98,27.52,27.78,28.24,28.78'yy=-1.4548+0.1763*x;x=127.3,130.0,132.7,129.4,135,137.1,141.2,142.8,145.5,145.3,148.3,146.4,150.2,153.1,157.3,160.7,164.2