第十章 统计方法建模
10.1 牙膏的销售量
10.2 软件开发人员的薪金
10.3 酶促反应
10.4 投资额与国民生产总值和物价指数回归模型是用统计分析方法建立的最常用的一类模型数学建模的基本方法 机理分析 测试分析通过对数据的统计分析,找出与数据拟合最好的模型
不涉及回归分析的数学原理和方法
通过实例讨论如何选择不同类型的模型
对软件得到的结果进行分析,对模型进行改进由于客观事物内部规律的复杂及人们认识程度的限制,
无法分析实际对象内在的因果关系,建立合乎机理规律的数学模型。
10.1 牙膏的销售量问题建立牙膏销售量与价格、广告投入之间的模型预测在不同价格和广告费用下的牙膏销售量收集了 30个销售周期本公司牙膏销售量、价格、
广告费用,及同期其它厂家同类牙膏的平均售价
9.260.556.804.253.7030
7.930.055.803.853.8029

8.510.256.754.003.752
7.38-0.055.503.803.851
销售量
(百万支 )
价格差
(元)
广告费用
(百万元 )
其它厂家价格 (元 )
本公司价格 (元 )
销售周期基本模型
y ~公司牙膏销售量
x1~其它厂家与本公司 价格差
x2~公司广告费用
110 xy
222210 xxy
5 5,5 6 6,5 7 7,5
7
7,5
8
8,5
9
9,5
10
x2
y
- 0,2 0 0,2 0,4 0,6
7
7,5
8
8,5
9
9,5
10
x1
y
22322110 xxxy
x1,x2~解释变量 (回归变量,自变量 )
y~被解释变量(因变量)
0,?1,?2,?3 ~回归系数
~随机 误差( 均值为零的正态分布随机变量)
MATLAB 统计工具箱模型求解
[b,bint,r,rint,stats]=regress(y,x,alpha)
输入
x= ~n?4数据矩阵,第 1列为全 1向量
]1[ 2221 xxx
alpha(置信 水平,0.05)
22322110 xxxy
b~?的 估计值
bint~b的置信区间
r~残差向量 y-xb
rint~r的置信区间
Stats~
检验统计量
R2,F,p
y~n维数据向量 输出由数据 y,x1,x2估计?
参数 参数估计值 置信区间
17.3244 [5.7282 28.9206]
1.3070 [0.6829 1.9311 ]
-3.6956 [-7.4989 0.1077 ]
0.3486 [0.0379 0.6594 ]
R2=0.9054 F=82.9409 p=0.0000
0
1
2
3
结果分析
y的 90.54%可由模型确定参数 参数估计值 置信区间
17.3244 [5.7282 28.9206]
1.3070 [0.6829 1.9311 ]
-3.6956 [-7.4989 0.1077 ]
0.3486 [0.0379 0.6594 ]
R2=0.9054 F=82.9409 p=0.0000
0
1
2
3
22322110 xxxy
F远超过 F检验的临界值
p远小于?=0.05
2的置信区间包含零点
(右端点距零点很近 )
x2对因变量 y 的影响不太显著
x22项显著 可将 x2保留在模型中模型从整体上看成立
22322110 xxxy销售量预测价格差 x1=其它厂家 价格 x3-本公司 价格 x4
估计 x3 调整 x4
控制价格差 x1=0.2元,投入广告费 x2=650万元销售量预测区间为 [7.8230,8.7636](置信度 95%)
上限用作库存管理的目标值 下限用来把握公司的现金流若估计 x3=3.9,设定 x4=3.7,则可以 95%的把握知道销售额在 7.8320?3.7? 29(百万元)以上控制 x1 通过 x1,x2预测 y
2 9 3 3.8 22322110 xxxy (百万支 )
模型改进
x1和 x2对 y
的 影响独立
22322110 xxxy
21422322110 xxxxxy
参数 参数估计值 置信区间
17.3244 [5.7282 28.9206]
1.3070 [0.6829 1.9311 ]
-3.6956 [-7.4989 0.1077 ]
0.3486 [0.0379 0.6594 ]
R2=0.9054 F=82.9409 p=0.0000
0
1
2
3
参数 参数估计值 置信区间
29.1133 [13.7013 44.5252]
11.1342 [1.9778 20.2906 ]
-7.6080 [-12.6932 -2.5228 ]
0.6712 [0.2538 1.0887 ]
-1.4777 [-2.8518 -0.1037 ]
R2=0.9209 F=72.7771 p=0.0000
3
0
1
2
4
x1和 x2对 y
的影响有交互作用两模型销售量预测 比较
21422322110 xxxxxy
22322110 xxxy
2933.8y (百万支 )
区间 [7.8230,8.7636]
区间 [7.8953,8.7592]
3272.8y (百万支 )
控制价格差 x1=0.2元,投入广告费 x2=6.5百万元预测区间长度更短略有增加y?
x2=6.5
x1=0.2
- 0,2 0 0,2 0,4 0,6
7,5
8
8,5
9
x1
y?
- 0,2 0 0,2 0,4 0,6
7,5
8
8,5
9
x1
y?
5 6 7 8
7,5
8
8,5
9
9,5
10
x2
y?
5 6 7 8
8
8,5
9
9,5
10
1 0,5
x2
y?
22322110 xxxy 21422322110 xxxxxy
两模型 与 x1,x2关系的 比较y?
交互作用影响的讨论
2221.0 6 7 1 2.07 5 5 8.72 2 6 7.30?
1 xxy x
价格差 x1=0.1
价格差 x1=0.3 2
223.0 6 7 1 2.00 5 1 3.84 5 3 5.32? 1 xxy x
21422322110 xxxxxy
5 3 5 7.72?x
加大广告投入使销售量增加
( x2大于 6百万元)
价格差较小时增加的速率更大
5 6 7 8
7,5
8
8,5
9
9,5
10
1 0,5
x 1 = 0,1x 1 = 0,3
x2
y?
1.03.0 11 xx yy
价格优势会使销售量增加价格差较小时更需要靠广告来吸引顾客的眼球完全二次多项式模型
22521421322110 xxxxxxy
MATLAB中有命令 rstool直接求解
0 0,2 0,4
7,5
8
8,5
9
9,5
10
5,5 6 6,5 7
x1 x2
y?
)?,?,?,?,?,?(? 543210从输出 Export 可得
10.2 软件开发人员的薪金资历 ~ 从事专业工作的年数;管理 ~ 1=管理人员,0=非管理人员;教育 ~ 1=中学,2=大学,3=更高程度建立模型研究薪金与资历、管理责任、教育程度的关系分析人事策略的合理性,作为新聘用人员薪金的参考编号薪金 资历管理教育
01 13876 1 1 1
02 11608 1 0 3
03 18701 1 1 3
04 11283 1 0 2

编号薪金 资历管理教育
42 27837 16 1 2
43 18838 16 0 2
44 17483 16 0 1
45 19207 17 0 2
46 19346 20 0 1
46名软件开发人员的档案资料
分析与假设 y~ 薪金,x1 ~资历(年)
x2 = 1~ 管理人员,x2 =0~ 非管理人员
1=中学
2=大学
3=更高

其它中学
,
,x
0
1
3

其它大学
,
,x
0
1
4
资历每加一年薪金的增长是常数;
管理、教育、资历之间无交互作用教育
443322110 xaxaxaxaay线性回归模型
a0,a1,…,a4是待估计的回归系数,?是随机误差中学,x3=1,x4=0 ;
大学,x3=0,x4=1;
更高,x3=0,x4=0
模型求解
443322110 xaxaxaxaay
参数 参数估计值 置信区间
a0 11032 [ 10258 11807 ]
a1 546 [ 484 608 ]
a2 6883 [ 6248 7517 ]
a3 -2994 [ -3826 -2162 ]
a4 148 [ -636 931 ]
R2=0.957 F=226 p=0.000
R2,F,p? 模型整体上可用资历增加 1年薪金增长 546
管理人员薪金多
6883
中学程度薪金比更高的少 2994
大学程度薪金比更高的多 148
a4置信区间包含零点,
解释不可靠 !
中学,x3=1,x4=0;大学,x3=0,x4=1; 更高:
x3=0,x4=0,
x2 =1~ 管理,x2
=0~ 非管理
x1~资历 (年 )
残差分析方法结果分析
443322110 xaxaxaxaay
残差 yye
e 与资历 x1的关系
0 5 10 15 20
- 2 0 0 0
- 1 0 0 0
0
1000
2000
e与管理 —教育组合的关系
1 2 3 4 5 6
- 2 0 0 0
- 1 0 0 0
0
1000
2000
残差全为正,或全为负,管理 —教育组合处理不当残差大概分成 3个水平,
6种管理 —教育组合混在一起,未正确反映 。 应在模型中增加管理 x2与教育x
3,x4的交互项组合 1 2 3 4 5 6
管理 0 1 0 1 0 1
教育 1 1 2 2 3 3
管理与教育的组合
426325443322110 xxaxxaxaxaxaxaay
进一步的模型 增加管理 x2与教育 x3,x4的交互项参数 参数估计值 置信区间
a0 11204 [11044 11363]
a1 497 [486 508]
a2 7048 [6841 7255]
a3 -1727 [-1939 -1514]
a4 -348 [-545 –152]
a5 -3071 [-3372 -2769]
a6 1836 [1571 2101]
R2=0.999 F=554 p=0.000
R2,F有改进,所有回归系数置信区间都不含零点,模型完全可用消除了不正常现象异常数据 (33号 )应去掉
0 5 10 15 20
- 1 0 0 0
- 5 0 0
0
500
e ~ x1
1 2 3 4 5 6
- 1 0 0 0
- 5 0 0
0
500
e ~组合去掉异常数据后 的结果参数 参数估计值 置信区间
a0 11200 [11139 11261]
a1 498 [494 503]
a2 7041 [6962 7120]
a3 -1737 [-1818 -1656]
a4 -356 [-431 –281]
a5 -3056 [-3171 –2942]
a6 1997 [1894 2100]
R2= 0.9998 F=36701 p=0.0000
0 5 10 15 20
- 2 0 0
- 1 0 0
0
100
200
e ~ x1
1 2 3 4 5 6
- 2 0 0
- 1 0 0
0
100
200
e ~组合
R2,0.957? 0.999? 0.9998
F,226? 554? 36701
置信区间长度更短残差 图十分正常最终模型的结果可以应用模型应用制订 6种管理 —教育组合人员的,基础,薪金 (资历为 0)
组合 管理 教育 系数 ―基础”薪金
1 0 1 a0+a3 9463
2 1 1 a0+a2+a3+a5 13448
3 0 2 a0+a4 10844
4 1 2 a0+a2+a4+a6 19882
5 0 3 a0 11200
6 1 3 a0+a2 18241
426325443322110 xxaxxaxaxaxaxaay
中学,x3=1,x4=0 ;大学,x3=0,x4=1; 更高,x3=0,x4=0
x1=0; x2 = 1~ 管理,x2 =0~ 非管理大学程度管理人员比更高程度管理人员的薪金高大学程度非管理人员比更高程度非管理人员的薪金略低对定性因素 (如管理、教育 ),可以 引入 0-1变量 处理,
0-1变量的个数应比定性因素的水平少 1
软件开发人员的薪金残差分析方法 可以发现模型的缺陷,引入交互作用项常常能够改善模型剔除异常数据,有助于得到更好的结果注:可以直接对 6种管理 —教育组合引入 5个 0-1变量
10.3 酶促反应问题研究酶促反应( 酶催化反应) 中嘌呤霉素对反应速度与底物 (反应物) 浓度之间关系的影响建立数学模型,反映该酶促反应的速度与底物浓度以及经嘌呤霉素处理与否之间的关系设计了两个实验,酶经过嘌呤霉素处理;酶未经嘌呤霉素处理。实验数据见下表,
方案底物浓度 (ppm) 0.02 0.06 0.11 0.22 0.56 1.10
反应速度处理 76 47 97 107 123 139 159 152 191 201 207 200
未处理 67 51 84 86 98 115 131 124 144 158 160 /
基本模型 Michaelis-Menten模型
y ~ 酶促反应的速度,x ~底物浓度
x
xxfy
2
1),(

1,?2 ~ 待定 系数底物浓度较小时,反应速度大致与浓度成正比;
底物浓度很大、渐进饱和时,反应速度趋于固定值。
酶促反应的基本性质
x
y
0
1
实验数据
0 0,5 1 1,5
0
50
100
150
200
250
经嘌呤霉素处理
x
y
0 0,5 1 1,5
0
50
100
150
200
250
未经嘌呤霉素处理
x
y
线性化模型经嘌呤霉素处理后实验数据的估计结果参数 参数估计值( × 10-3) 置信区间( × 10-3)
1 5.107 [3.539 6.676]
2 0.247 [0.176 0.319]
R2=0.8557 F=59.2975 p=0.0000
8 0 2 7.1 9 5?/1? 11 0 4 8 4 1.0?/ 122
x
xy
2
1
xy
111
1
2
1?

对?1,?2非线性 对?1,?2线性
x
1
21
线性化模型结果分析
x较大时,y有较大偏差1/x较小时有很好的线性趋势,1/x较大时出现很大的起落? 参数估计时,x较小( 1/x很大)的数据控制了回归参数的确定
0 10 20 30 40 50
0
0,0 0 5
0,0 1
0,0 1 5
0,0 2
0,0 2 5
1/y
1/x
xy
11
21
0 0,5 1 1,5
0
50
100
150
200
250
x
xy
2
1
x
y
[beta,R,J] = nlinfit (x,y,’model’,beta0)
beta的置信区间
MATLAB 统计工具箱输入 x~自变量 数据矩阵
y ~因变量数据向量
beta ~参数的估计值
R ~残差,J ~估计预测误差的 Jacobi矩阵
model ~模型的函数 M文件名
beta0 ~给定的参数初值输出
betaci =nlparci(beta,R,J)
非线性模型参数估计
function y=f1(beta,x)
y=beta(1)*x./(beta(2)+x);
x
xy
2
1
x= ; y= ;
beta0=[195.8027 0.04841];
[beta,R,J]=nlinfit(x,y,’f1’,beta0);
betaci=nlparci(beta,R,J);
beta,betaci
beta0~线性化模型估计结果非线性模型结果分析参数 参数估计值 置信区间
1 212.6819 [197.2029 228.1609]
2 0.0641 [0.0457 0.0826 ]
画面左下方的 Export
输出其它统计结果。
拖动画面的十字线,得
y的预测值和预测区间剩余标准差 s= 10.9337
x
xy
2
1
最终反应速度为半速度点 (达到最终速度一半时的 x值 )为
6 8 3 1.2 1 2?1
0641.0? 2
其它输出 命令 nlintool 给出交互画面
0 0,5 1 1,5
0
50
100
150
200
250
o ~原始数据
+ ~ 拟合结果
0 0,2 0,4 0,6 0,8 1
- 5 0
0
50
100
150
200
250
混合反应模型
x1为底物浓度,x2为一示性变量
x2=1表示经过处理,x2=0表示未经处理
β1是未经处理的最终反应速度
γ1是经处理后最终反应速度的增长值
β2是未经处理的反应的半速度点
γ2是经处理后反应的半速度点的增长值在同一模型中考虑嘌呤霉素处理的影响
x
xy
2
1
1222
1211
)( xx
xxy



)(
o ~原始数据
+ ~拟合结果混合模型求解 用 nlinfit 和 nlintool命令
,17001,6001,05.002 01.002
估计结果和预测剩余标准差 s= 10.4000
参数 参数估计值 置信区间
1 160.2802 [145.8466 174.7137]
2 0.0477 [0.0304 0.0650 ]
1 52.4035 [32.4130 72.3941 ]
2 0.0164 [-0.0075 0.0403]
2置信区间包含零点,表明?2对因变量 y的影响不显著
1222
1211
)( xx
xxy



)(
参数初值 (基于对数据的分析 )
经嘌呤霉素处理的作用不影响半速度点参数未经处理经处理
o ~原始数据
+ ~拟合结果未经处理经处理简化的混合模型简化的混合模型 形式简单,参数置信区间 不含零点剩余标准差 s = 10.5851,比一般混合模型略大
1222
1211
)( xx
xxy



)(
12
1211
x
xxy

)(
估计结果和预测参数参数估计值 置信区间
1 166.6025 [154.4886 178.7164]
2 0.0580 [0.0456 0.0703 ]
1 42.0252 [28.9419 55.1085]
一般混合模型与简化混合模型预测比较实际值 一般模型预测值 Δ(一般 模型 ) 简化模型预测值 Δ(简化 模型 )
67 47.3443 9.2078 42.7358 5.4446
51 47.3443 9.2078 42.7358 5.4446
84 89.2856 9.5710 84.7356 7.0478
… … … … …
191 190.8329 9.1484 189.0574 8.8438
201 190.8329 9.1484 189.0574 8.8438
207 200.9688 11.0447 198.1837 10.1812
200 200.9688 11.0447 198.1837 10.1812
简化混合模型的预测区间较短,更为实用、有效
1222
1211
)( xx
xxy



)(
12
1211
x
xxy

)(
预测区间为预测值?Δ
注:非线性模型拟合程度的评价无法直接利用线性模型的方法,但 R2 与 s仍然有效。
酶促反应反应速度与底物浓度的关系 非线性 关系求解 线性模型 求解非线性模型机理分析嘌呤霉素处理对反应速度与底物浓度关系的影响混合模型发现问题,
得参数初值引入 0-1变量简化模型检查 参数置信区间 是否包含零点
10.4 投资额与国民生产总值和物价指数问题建立投资额模型,研究 某地区 实际投资额与国民生产总值 ( GNP ) 及物价指数 ( PI ) 的关系
2.06883073.0424.5201.00001185.9195.010
1.95142954.7474.9190.96011077.6166.49
1.78422631.7401.9180.9145992.7144.28
1.63422417.8423.0170.8679944.0149.37
1.50422163.9386.6160.8254873.4133.36
1.40051918.3324.1150.7906799.0122.85
1.32341718.0257.9140.7676756.0125.74
1.25791549.2206.1130.7436691.1113.53
1.15081434.2228.7120.7277637.797.42
1.05751326.4229.8110.7167596.790.91
物价指数国民生产总值投资额年份序号物价指数国民生产总值投资额年份序号根据对未来 GNP及 PI的估计,预测未来投资额该地区 连续 20年的统计数据时间序列中同一变量的顺序观测值之间存在 自相关以时间为序的数据,称为 时间序列分析许多经济数据在时间上有一定的 滞后 性需要诊断并消除数据的自相关性,建立新的模型若采用普通回归模型直接处理,将会出现不良后果投资额与国民生产总值和物价指数
…………………… 1.32341718.0257.9140.7676756.0125.74
1.25791549.2206.1130.7436691.1113.53
1.15081434.2228.7120.7277637.797.42
1.05751326.4229.8110.7167596.790.91
物价指数国民生产总值投资额年份序号物价指数国民生产总值投资额年份序号基本回归模型投资额与 GNP及物价指数间均有很强的线性关系
tttt xxy 22110
t ~年份,yt ~ 投资额,x1t~ GNP,x2t ~ 物价指数
0,?1,?2 ~回归系数
x1t
yt
x2t
yt
t ~对 t相互 独立的零均值正态随机变量基本回归模型的结果与分析
ttt xxy 21 479.8596185.0725.322
MATLAB 统计工具箱参数 参数估计值 置信区间
0 322.7250 [224.3386 421.1114]
1 0.6185 [0.4773 0.7596]
2 -859.4790 [-1121.4757 -597.4823 ]
R2= 0.9908 F= 919.8529 p=0.0000
剩余标准差
s=12.7164
没有考虑时间序列数据的 滞后性影响
R2= 0.9908,拟合度高模型优点模型缺点可能忽视了随机误差存在 自相关 ;如果存在自相关性,用此模型会有不良后果自相关性的定性诊断 残差诊断法
ttt yye
模型残差作残差 et~et-1 散点图大部分点落在第 1,3象限?t存在正的自相关大部分点落在第 2,4象限自相关性直观判断在 MATLAB工作区中输出
et为随机误差?t 的估计值
- 3 0 - 2 0 - 1 0 0 10 20
- 3 0
- 2 0
- 1 0
0
10
20
et-1
et
t存在负的自相关基本回归 模型的随机误差项?t 存在正的自相关自回归 性 的 定量诊断自回归模型
ttttttt uxxy 122110,
ρ~自相关系数 1||0,?1,?2 ~回归系数
ρ= 0 无 自相关性
ρ> 0
ρ< 0
如何估计 ρ
如何消除自相关 性
D-W统计量
D-W检验
ut ~对 t相互 独立的零均值正态随机变量存在负 自相关性存在正 自相关性广义差分法
D-W统计量与 D-W检验

n
t
t
n
t
tt
e
ee
DW
2
2
2
2
1 )(
检验水平,样本容量,
回归变量数目
D-W分布 表
n
t
t
n
t
tt
e
ee
2
2
2
1
12 )(12
n较大



n
t
t
n
t
tt eee
2
2
2
1 /
401?1 DW?
DW4-dU 44-dLdUdL 20
正自相关负自相关不能确定不能确定无自相关
20 DW?
01 DW? 41 DW?
检验 临界值 dL和 dU 由 DW值的大小确定 自相关性广义差分变换
)1(0*0
以?*0,?1,?2 为 回归系数的普通回归模型原模型
DW值
D-W
检验无自相关有自相关广义差分 继续此过程原模型新模型新模型
tttt uxxy *22*11*0*
步骤原模型
ttttttt uxxy 122110,
,1* ttt yyy? 2,1,1,* ixxx tiitit?变换
)(12DW 21? DW
不能确定 增加数据量;选用其它方法投资额新模型的建立
DWold < dL
作变换原模型残差 et
样本容量 n=20,回归变量数目 k=3,?=0.05
查表临界值 dL=1.10,dU=1.54
DWold=
0.8754
原模型有正自相关
1* 5 6 2 3.0 ttt yyy
2,1,5623.0 1,* ixxx tiitit

n
t
t
n
t
tt
e
ee
DW
2
2
2
2
1 )(
5623.02/1 DW?
DW4-dU 44-dLdUdL 20
正自相关负自相关不能确定不能确定无自相关参数 参数估计值 置信区间
*0 163.4905 [1265.4592 2005.2178]
1 0.6990 [0.5751 0.8247]
2 -1009.0333 [-1235.9392 -782.1274]
R2= 0.9772 F=342.8988 p=0.0000
tttt uxxy *22*11*0*
21*0*2*1*,,,,估计系数由数据 ttt xxy
总体效果良好 剩余标准差
snew= 9.8277 < sold=12.7164
投资额新模型的建立
1* 5 6 2 3.0 ttt yyy 2,1,5623.0 1,* ixxx tiitit
新模型的自相关性检验
dU< DWnew < 4-dU
新模型残差 et
样本容量 n=19,回归变量数目 k=3,?=0.05
查表临界值 dL=1.08,dU=1.53
DWnew=
1.5751
新模型无自相关性
DW4-dU 44-dLdUdL 20
正自相关负自相关不能确定不能确定无自相关
1,2,2
1,1,11
3794.5670333.1009
3930.0699.05623.04905.163?



tt
tttt
xx
xxyy
*2*1* 033.1009699.04905.163? ttt xxy新模型还原为原始变量一阶自回归模型一阶自回归 模型残差 et比 基本回归 模型要小
0 5 10 15 20
- 3 0
- 2 0
- 1 0
0
10
20
新模型 et~ *,原模型 et~ +
残差图比较
0 5 10 15 20
0
100
200
300
400
500
新模型?t ~ *,新模型?t ~ +
拟合图比较模型结果比较
ttt xxy 21 4 7 9.8 5 96 1 8 5.07 2 5.3 2 2
基本回归模型一阶自回归模型
1,2,2
1,1,11
3794.5670333.1009
3930.0699.05623.04905.163?



tt
tttt
xx
xxyy
投资额预测对未来投资额 yt 作预测,需先 估计出未来的国民生产总值 x1t 和物价指数 x2t
设已知 t=21时,x1t =3312,x2t=2.1938
7 6 3 8.4 6 9ty一阶自回归模型
2.06883073.0424.520
1.95142954.7474.919
1.78422631.7401.918
0.7436691.1113.53
0.7277637.797.42
0.7167596.790.91
物价指数国民生产总值投资额年份序号物价指数国民生产总值投资额年份序号一阶自回归模型 7 6 3 8.4 6 9
ty
基本回归模型 6 7 2 0.4 8 5
ty
t 较小是由于 yt-1=424.5过小所致