数学建模之统计回归模型.doc

数学建模大作业

统计回归模型

摘要
某公司想用全行业的销售额作为自变量来预测公司的销售额，题目给出了1977—1981此公司的销售额和行业销售额的分季度数据表格。通过对所给数据的简单分析，我们可以看出：此公司的销售额有随着行业销售额的增加而增加的趋势，为了更加精确的分析题目所给的数据，得出科学的结论，从而达到合理预测的目的。我们使用时间序列分析法，参照课本统计回归模型例4，做出了如下的统计回归模型。
在问题一中，我们使用MATLB数学软件，画出了数据的散点图，通过观察散点图，发现公司的销售额和行业销售额之间有很强的线性关系，于是我们用线性回归模型去拟合，发现有很好的拟合性。但是这种情况下，并没有考虑到数据的自相关性，所以我们做了下面几个问题的分析来对这个数学模型进行优化。
在问题二中，通过建立了公司销售额对全行业销售额的回归模型，并使用DW检测诊断随机误差项的自相关性。通过计算和查DW表比较后发现随即误差存在正自相关，也就是说前面的模型有一定的局限性，预测结果存在一定的偏差，还有需要改进的地方。
在问题三中，因为在问题二中得出随即误差存在正自相关，为了消除随机误差的自相关性，我们建立了一个加入自相关后的回归模型。并对其作出了分析和验证，我们发现加入自相关后的回归模型更加合理。通过使用我们建立的模型对公司的销售额进行预测，发现和实际的销售额很接近，也就是说模型效果还不错。
关键词：销售额、回归模型、自相关性
一、问题提出
某公司想用全行业的销售额作为自变量来预测公司的销售额，下表给出了1977-1981年公司销售额和行业销售额的分季度数据（单位：百万元）.
（1） 画出数据的散点图，观察用线性回归模型拟合是否合适。
（2） 监理公司销售额对全行业销售额的回归模型，并用DW检验诊断随机误差项的自相关性。
（3） 建立消除了随机误差项自相关性后的回归模型。
年
季

公司销售额
行业销售额
年
季

公司销售额
行业销售额
1977

1978

1979
1
2
3
4
1
2
3
4
1
2
1
2
3
4
5
6
7
8
9
10
20.96
21.40
21.96
21.52
22.39
22.76
23.48
23.66
24.10
24.01
127.3
130.0
132.7
129.4
135.0
137.1
141.2
142.8
145.5
145.3
1979

1980

1981
3
4
1
2
3
4
1
2
3
4
11
12
13
14
15
16
17
18
19
20
24.54
24.30
25.00
25.64
26.36
26.98
27.52
27.78
28.24
28.78
148.3
146.6
150.2
153.1
157.3
160.7
164.2
165.6
168.7
171.7

二、基本假设
假设一：模型中ε（对时间t）相互独立。
三、符号说明
公司销售额：（百万）
行业销售额：（百万）
概念介绍：
1.自相关：自相关（auto correlation），又称序列相关（serial correlation）是指总体回归模型的随机误差项之间存在的相关关系。即不同观测点上的误差项彼此相关。
2．置信区间：如果P（）=1-α，α=0.1或0.05，则称区间[a,b]为的置信度为1-α的置信区间。
3.时间序列：时间序列法是一种定量预测方法，亦称简单外延方法。时间序列即按时间的推移或排布会对规律的变化有所影响。
四、问题分析
问题一：表中的数据是以时间为顺序的。由于前期的销售额对后期的投资一般有明显的影响，从而对后期的后期的销售额造成影响。因此在此模型中应考虑到存在自相关，我们可以先建立基本的回归模型，然后再进行自相关性诊断，并建立新的回归模型。
问题二：在问题一之后，就可以接着求出问题二，然后利用DW检验诊断随机误差项的自相关性。
问题三：进行了自相关诊断后，将自相关加入模型中，建立消除了随机误差项自相关性的回归模型。
五、模型的建立与求解
5.1 问题一
5.1.1 问题一的分析
表中数据是以时间为序的，建立基本的回归模型。
5.1.2 问题一模型的建立
基本回归模型：
设该公司第时间的公司销售额为，行业销售额为 。为了大致分析和的关系，首先利用表中的数据作出对关系作出散点图，如下（见图中的“+”）：做散点图：

可以看出，随着行业销售额的增加，公司销售额增大，而且两者有很强的线性关系，图中的直线说明两者呈线性模型，因此本题用线性回归模型拟合非常合适。
5.2 问题二
5.2.1 问题二的分析
从问题一中的图形可以看出，随着行业销售额的增加，公司销售额增大，而且两者有很强的线性关系，图中的直线说明两者呈线性模型，因此可建立一元线性回归模型。
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.0e+004 *0.0001
F=1.0e+004 *1.4888
P=1.0e+004 *0.0000

将参数估计值代入（1）得到：
(2)
用MATLAB中rstool命令得到的交互式画面见图 （1） ，由此可以得出不同水平下的预测值及其置信区间。通过左下方的Export下拉式菜单。可以输出模型的统计结果。

图1
自相关性诊断与处理方法 从表面上来看得到的基本模型（2）拟合度（R）非常之高，接近你100%，应该很满意了，但是，这个模型并没有考虑到我们的数据是一个时间序列（将原表中的数据打乱不影响模型（2）的结果）。实际上对于时间序列数据做回归分析时，模型的随机误差有可能存在相关性，违背模型关于（对时间t）相互独立的基本假设，其他相关因素对公司销售额的影响肯能也有时间上的延续，包含在随机误差中，即随机误差会出现自相关性。
残差可以作为随机误差的估计值，画出的散点图，能够从直观上判断的自相关性。模型（2）的残差可在计算过程中得到表1，以及数据的图见图 2
t
1
2
3
4
5
e
-0.0282
-0.0642
0.0198
0.1616
0.0443
t
6
7
8
9
10
e
0.0441
0.0412
-0.0608
-0.0968
-0.1516
t
11
12
13
14
15
e
-0.1505
-0.0555
-0.0255
0.1033
0.0828
t
16
17
18
19
20
e
0.1034
0.0263
0.0395
-0.047
-0.0359
表 1

图 2
为了对ε的字相关性做定量的诊断，并在确诊后得到新的结果，我们考虑如下模型
， （3）
其中是自相关系数，||1,相互独立且服从均值为0的正态分布。
若=0，则退化为普通的回归模型；若>0，则随机误差存在正的自相关；若<0，则随机误差存在负的自相关。
利用D-W检验诊断自相关现象如下：
利用MATLAB算出:
=0.0980 =0.1326
DW=0.7388 =0.6306
(具体程序见附录)
因为DW≈2（1-），所以 0≤DW≤4,若的估计值在0附近，则DW的值在2附近，的自相关行很弱，若在正负1附近，则DW接近0或4，的自相关性很强。
5.2.4 问题二结果的分析及验证
要根据DW的具体数值确定是否存在自相关，查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）是普通回归模型。
以的估计值带入（3）和（4）做变换，利用变换后的数据 ，估计模型（5）的参数，得到的表见表2，还可以得出剩余标准差rmse=0.08828.
最后将模型（5）的变量还原为原始变量。得到的结果如下
（6）
表 2
参数
参数估计值
参数置信区间

-0.3951
[ -0.7481 -0.0422]

0.1738
[0.1675 0.1800]
R=1.0e+003 *0.0010
F=1.0e+003 *3.4621
P=1.0e+003 *0.0000

5.3.4 问题三结果的分析及验证
当然应该对模型（6）也作一次自相关检验，即诊断随机误差是否还存在自相关，从模型（6）的残差可以计算出DW=1.65，对于显著水平α=0.05，n=19,k=2,查D-W分布表，得到检验的临界值=1.2和=1.40 .现在,可以认为随机误差不存在自相关。一次经变换得到的回归模型（6）是适用的。
结果及其预测
从机理上看，对于带滞后性的经济规律作用下的时间序列数据，加入自相关的模型（6）更为合理，我们将模型（1）与模型（6）的计算值与实际数据的比较，以及两个模型的残差，表示在表 3 中
表 3
t
y(实际数据)
yy(模型1)
yyy（模型2）
ee
e1
2
21.4
21.464
21.464
-0.0642
4.00E-06
3
21.96
21.94
21.915
0.01979
0.02521
4
21.52
21.358
21.399
0.16158
-0.04026
5
22.39
22.346
22.456
0.0443
-0.11047
6
22.76
22.716
22.756
0.04407
-0.04008
7
23.48
23.439
23.472
0.04124
-0.033
8
23.66
23.721
23.755
-0.06084
-0.03367
9
24.1
24.197
24.162
-0.09685
0.034934
10
24.01
24.162
24.109
-0.15159
0.05289
11
24.54
24.69
24.595
-0.15049
0.095224
12
24.3
24.356
24.27
-0.05552
0.085056
13
25
25.025
24.988
-0.02546
0.03766
14
25.64
25.537
25.517
0.10327
0.01997
15
26.36
26.277
26.332
0.08281
-0.05527
16
26.98
26.877
26.917
0.10339
-0.04049
17
27.52
27.494
27.544
0.02634
-0.05007
18
27.78
27.74
27.744
0.03952
-0.00349
19
28.24
28.287
28.293
-0.04701
-0.00626
20
28.78
28.816
28.765
-0.03591
0.050926

六、模型的评价与推广
模型的评价与推广：此模型从最初的线性回归模型到DW检验诊断随机误差的自相关性，再到最后的消除了随机误差项自相关性后的回归模型，模型逐步得到了优化。最后的结果预测可以看出，我们建立的这个模型的可靠性是非常高的。预测公司的销售额可以为公司的制定相应的生产计划或者购货数量提供依据，鉴于该模型的可靠性非常稳定，我们可以把此模型推广到公司其他产品的销售额或者某产品的市场销售额的预测。
同时，通过这学期的建模课程的学习嘛，我们发现团队精神是数学建模是否取得好成绩的最重要的因素，一队三个人要相互支持，相互鼓励。切勿自己只管自己的一部分（数学好的只管建模，计算机好的只管编程，写作好的只管论文写作），很多时候，一个人的思考是不全面的，只有大家一起讨论才有可能把问题搞清楚，因此无论做任何板块，三个人要一起齐心才行，只靠一个人的力量，要写出一篇高水平的文章几乎是不可能的。其实建模的过程就是大家互相鼓励，共同勉励的一个阶段，我们组从最开始的模拟训练时就十分注重团队的分工协作，在作业中都会总结教训，改进方法。另外我们还根据每个人的特长来进行分工，做到发挥优势，长短互补的效果。
七、参考文献
《数学模型》（第三版） 姜启源 谢金星 叶俊 高等教育出版社，2003年8月
数学建模上课课件：统计回归模型
《应用回归分析》 何晓群，刘文清 中国人民大学出版社，2001
网上资源 DW表
八、附录
MATLAB运行程序
附录1：
>> 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]';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
yy =20.9882
21.4642
21.9402
21.3584
22.3457
22.7159
23.4388
23.7208
24.1969
24.1616
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.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*DW
p =
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.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)
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.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.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/y1
DW =
0.7388
>>
>> 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,:)
x2=x(2:20,:);x3=x(1:19,:);y4=y2-y3*p;
x4=x2-x3*p
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.7800
y3 =
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.9800
27.5200
27.7800
28.2400
x4 =
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,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,:)
x2=x(2:20,:);x3=x(1:19,:);y4=y2-y3*p;
x4=x2-x3*p;[b1,bint1,r1,rint1,stats1]=regress(y4,[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.7800
y3 =
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.9800
27.5200
27.7800
28.2400
b1 =
-0.3951
0.1738
bint1 =
-0.7481 -0.0422
0.1675 0.1800
r1 =
-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.0170
rint1 =
-0.1941 0.0688
-0.0886 0.1817
0.0163 0.2291
-0.2012 0.0721
-0.1337 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.1468
stats1 =
1.0e+003 *
0.0010 3.4621 0 0.0000
>>
>>
>> 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,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
e1 =
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.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.3951+0.6306*y3+0.1738*x2-0.1096*x3;yyyy=yy(2:20,:);e1=yyyy-yyy;e=y-yy;ee=e(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,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
yyy =
21.4642
21.9150
21.3987
22.4562
22.7560
23.4718
23.7545
24.1619
24.1087
24.5953
24.2705
24.9878
25.5168
26.3325
26.9171
27.5437
27.7440
28.2933
28.7650
>>
>> 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]';yy=-1.4548+0.1763*x
yy =
20.9882
21.4642
21.9402
21.3584
22.3457
22.7159
23.4388
23.7208
24.1969
24.1616
24.6905
24.3555
25.0255
25.5367
26.2772
26.8766
27.4937
27.7405
28.2870
28.8159
>> 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,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;e1=yyyy-yyy;t=2:20;subplot(1,2,1);plot(y2,yyy,'+')
hold on;yyyy=yy(2:20,:);plot(y2,yyyy,'o')
hold off;subplot(1,2,2);plot(t,ee,'+')
hold on;plot(t,e1,'o')

17