数学建模
第十三讲 统计模型
数学实验
与
模型一 牙膏的销售量
确定关系:
牙膏销售量——价格、广告投入
内部规律复杂数据统计分析
常用模型回归模型×数学原理软件
30个销售周期数据:
销售量、价格、广告费用、同类产品均价
30
29
…
…
…
…
…
…
2
1
销售量(百万支)
价差(元)
广告(百万元)
它厂价 (元)
公司价 (元)
销售周期
分析
y ~公司牙膏销售量
x1~其它厂家与本公司价格差
x2~公司广告费用
解释变量
(回归变量, 自变量)
被解释变量(因变量)
y与x1
y与 x2
Matlab 统计分析
版本 s增加一个统计量: 剩余方差s2
[b , bint , r , rint , stats] = regress( y , X , alpha )
3、回归分析
statistics toolbox
解释变量:矩阵
显著性水平:
系数估计值
置信区间
残差
置信区间
被解释变量:列
检验统计量:R2,F,p
随机误差:正态分布均值为零
回归系数
rcoplot(r,rint)
残差及其置信区间作图
于是
考察
y与x1
y与 x2
有
结果
模型
多元回归模型
MATLAB
t1
x=[ones(size(x1)),x1,x2,x2.^2];
[b,bint,r,rint,stats]=regress(y,x)
程序
结果分析
即:
R2= F= p=
[ ]
3
[ ]
2
[ ]
1
[ ]
0
置信区间
参数估计值
参数
显著性 :整体显著
y的%可由模型确定、 F远超过F检验的临界值、 p远小于=
x2 :2置信区间包含零点——不显著;3 显著
销售量预测
控制价格差x1=元,投入广告费x2=650万元
得:
价差x1=它厂价x3-公司价x4
估计x3,调整x4
销售量预测区间为 [,](置信度95%)?
上限用作库存管理的目标值
若估计x3=,设定x4=
可以95%的把握知道销售额在 29(百万元)以上
(百万支)
控制x1
预测y
模型改进
x1和x2对y的影响有交互作用
交互项
3
置信区间
R2= F= p=
估计值
置信区间
4
R2= F= p=
2
1
0
估计值
参数
比较:置信区间、 R2
比较:销售量预测
控制价格差x1=元,投入广告费x2=百万元
(百万支)
区间 [,]
区间 [,]
(百万支)
预测区间长度更短
略有增加
比较:两模型y与x1,x2关系
x2=
x1=
x1
x1
x2
x2
讨论:交互作用影响
价格差 x1=
价格差 x1=
广告投入y
( x2大于6百万元)
价格差较小时
增加的速率更大
x2
价格优势y
价格差较小广告作用大
x1
x2
完全二次多项式模型
MATLAB
相应面分析
Rstool(x,y,’model’,alpha,’xname’,’yname’)
linear
interaction
quadratic
purequadratic
线性项
常数项、线性项、交叉项
交叉项、平方项
常数项、线性项、平方项
相应面分析 rstool(x(:,2:3),y)
all: beta rmse residuals
系数 均方差 残差
x1
x2
Export :
评注
回归模型
数据、经验、图形
回归变量、函数形式
Matlab求解
统计分析: R2、F、p
回归系数置信区间包含0点改进
添加二次项、交叉项……
模型二 软件开发人员的薪金
薪金——资历、岗位、学历
建立模型:分析人事策略的合理性,作为新聘用人员薪金的参考
资历~ 从事专业工作的年数;管理~ 1=管理人员,0=非管理人员;教育~ 1=中学,2=大学,3=更高程度
46名软件开发人员的档案资料
1
0
20
19346
46
2
0
17
19207
45
…
…
…
…
…
3
0
1
11608
02
1
1
1
13876
01
教育
管理
资历
薪金
编号
模型假设
假设: y~ 薪金,x1 ~资历(年)
x2 = 1~ 管理人员,0~ 非管理人员
1 ~中学 2 ~大学 3 ~更高?
假设: 资历每加一年薪金的增长是常数;
管理、教育、资历之间无交互作用
教育=
模型:线性回归
回归系数 随机误差
中学:x3=1, x4=0 ;大学:x3=0, x4=1; 更高:x3=0, x4=0
模型求解
x1~资历(年)
x2 = 1~ 管理,0~ 非管理
中学:x3=1, x4=0;大学:x3=0, x4=1; 更高:x3=0, x4=0
Matlab程序:
:
序号、工资y、资历x1、管理x2、学历、x3、x4、xx
:
M=dlmread('');
x1=M(:,3);x2=M(:,4);x3=M(:,6);x4=M(:,7);y=M(:,2);
x=[ones(size(x1)) x1 x2 x3 x4 ]
[b,bi,r,ri,s]=regress(y,x)
结果
R2,F, p 模型整体上可用
资历增加1年薪金增长546
管理人员薪金多6883
中学程度薪金比更高的少2994
大学程度薪金比更高的多148
a4置信区间包含零点 解释不可靠!
置信区间
估计值
参数
R2= F=226 p=
[ -636 931 ]
148
a4
[-3826 -2162]
-2994
a3
[6248 7517]
6883
a2
[ 484 608 ]
546
a1
[10258 11807]
11032
a0
结果分析
残差分析法
残差
与资历x1的关系
残差大概分成3个水平
6种管理—教育组合混在一起,未正确反映
Matlab:
残差分析
与管理x2—教育x3 、x4的关系
残差全为正,或全为负,管理—教育组合处理不当
应在模型中增加管理x2与教育x3, x4的交互项
组合
1
2
3
4
5
6
管理
0
1
0
1
0
1
教育
1
1
2
2
3
3
管理与教育的组合
模型改进
增加管理x2与教育x3, x4的交互项
R2,F有改进
回归系数置信区间不含零点
模型可用
[-545 –152]
-348
a4
[-3372 -2769]
-3071
a5
置信区间
估计值
参数
R2= F=554 p=
[1571 2101]
1836
a6
[-1939 -1514]
-1727
a3
[6841 7255]
7048
a2
[486 508]
497
a1
[11044 11363]
11204
a0
Matlab:
残差分析
消除了不正常现象
异常数据(33号)去掉
e ~ x1
e ~组合
Matlab:
模型改进
去掉异常数据后的结果
R2:
F: 226
554
36701
置信区间长度更短
[-431 –281]
-356
a4
[-3171 –2942]
-3056
a5
置信区间
估计值
参数
R2= F=36701 p=
[1894 2100]
1997
a6
[-1818 -1656]
-1737
a3
[6962 7120]
7041
a2
[494 503]
498
a1
[11139 11261]
11200
a0
残差分析
残差图正常
模型的结果可以应用
~ x1
~组合
模型应用
制订基础薪金
资历为0 : x1= 0
管理—教育组合:6种
大学程度管理人员比更高程度管理人员的薪金高
大学程度非管理人员比更高程度非管理人员的薪金略低
18241
11200
19882
10844
13448
9463
基础薪金
1
0
1
0
1
0
管理
管理+更高
非管理+更高
管理+大学
非管理+大学
管理+中学
非管理+中学
a0+a2
3
系数
教育
组合
6
a0
3
5
a0+a2+a4+a6
2
4
a0+a4
2
3
a0+a2+a3+a5
1
2
a0+a3
1
1
教育
1中学:x3=1, x4=0
2大学:x3=0, x4=1
3更高:x3=0, x4=0
评注
对定性因素:如管理、教育
可以引入0-1变量处理
0-1变量的个数应比定性因素的水平少1
残差分析:可以发现模型的缺陷
引入交互作用项常常能够改善模型
剔除:异常数据
有助于得到更好的结果
另:
可以直接对6种管理—教育组合引入5个0-1变量
模型三 酶促反应
酶促反应
由酶作为催化剂催化进行的化学反应
生物体内的化学反应绝大多数属于酶促反应
酶促反应中酶作为高效催化剂使得反应以极快的速度(103~1017倍)或在一般情况下无法反应的条件下进行
酶是生物体内进行各种化学反应最重要的因素
酶促反应动力学
问题:
研究酶促反应中——嘌呤霉素——对反应速度与底物(反应物)浓度之间关系的影响
方案
建立数学模型,反映该酶促反应的速度与底物浓度以及经嘌呤霉素处理与否之间的关系
设计了两个实验
酶经过嘌呤霉素处理
酶未经嘌呤霉素处理
实验数据:
反应速度
115
98
86
84
51
67
未处理
反应速度
139
123
107
97
47
76
处理
200
207
201
191
152
159
处理
底物浓度(ppm)
/
160
158
144
124
131
未处理
底物浓度(ppm)
分析
酶促反应的基本性质
底物浓度较小时,反应速度大致与浓度成正比;
底物浓度很大、渐进饱和时,反应速度趋于固定值
Michaelis-Menten模型
待定系数 =(1 , 2)
基本模型
酶促反应的速度
底物浓度
数据分析
实验数据:散点图
y ~酶促反应速度
x ~底物浓度
x
y
0
1
经嘌呤霉素处理
x
y
未经嘌呤霉素处理
x
y
模型
线性化模型
经嘌呤霉素处理后实验数据的估计结果
估计值
β2=2 /1
β1=1/1
参数
R2= F= p=
[ ]
2
[ ]
1
置信区间
(×10-3)
估计值
(×10-3)
参数
对1 , 2非线性
对1, 2线性
结果分析
x较大时,y有较大偏差
1/x较小时有很好的线性趋势,1/x较大时出现很大的起落
线性化:参数估计时
x较小(1/x很大)的数据控制了回归参数的确定
改进:非线性模型
1/y
1/x
x
y
Matlab 统计分析
beta的置信区间
[beta , R , J] = nlinfit ( x,y, ’model’ , beta0 )
3、回归分析:非线性
statistics toolbox
解释变量:矩阵
模型函数
参数估计值
残差
参数初值
被解释变量:列
预测误差的Jacobi矩阵
betaci =nlparci(beta,R,J)
Matlab 程序
function y=f1(beta, x)
y=beta(1)*x./(beta(2)+x);
x=………… ; y=…………;
beta0=[ ];
[beta,R,J]=nlinfit(x,y,’f1’,beta0);
betaci=nlparci(beta,R,J);
beta, betaci
[ ]
[ ]
结果分析
最终反应速度为
半速度点(达到最终速度一半时的x值 )为
效果
[ ]
2
[ ]
1
置信区间
估计值
参数
o -原始数据
+ -拟合结果
其它输出
给出交互画面
拖动画面的十字线,得y的预测值和预测区间
画面左下方的Export 输出其它统计结果。
剩余标准差s=
nlintool (x,y,’model’,beta)
混合反应模型
在同一模型中考虑嘌呤霉素处理的影响
底物浓度
示性变量
x2示性变量:x2=1表示经过处理,x2=0表示未经处理
未经处理的最终反应速度
经处理后最终反应速度增长值
未经处理的反应的半速度点
经处理后反应的半速度点增长值
Matlab 程序
用nlinfit 和 nlintool命令
参数初值:基于对数据的分析
o ~原始数据
+ ~拟合结果
[ ]
2
[ ]
1
[ ]
2
[ ]
1
置信区间
估计值
参数
经嘌呤霉素处理的作用不影响半速度点参数
未经处理
经处理
估计结果和预测
置信区间包含零点,对因变量影响不显著
简化混合模型
估计结果和预测
o ~原始数据
+ ~拟合结果
未经处理
经处理
简化的混合模型形式简单
参数置信区间不含零点
剩余标准差 s = ,比一般混合模型略大
[ ]
1
[ ]
2
[ ]
1
置信区间
估计值
参数
结果分析
一般混合模型与简化混合模型预测比较
200
207
201
191
…
…
…
…
…
84
51
67
Δ
简化
Δ
一般
实际值
简化混合模型的预测区间较短,更为实用、有效
评注
酶促反应
注:非线性模型拟合程度的评价无法直接利用线性模型的方法,但R2 与s仍然有效。
反应速度与底物浓度的关系
非线性关系
求解线性模型
求解非线性模型
机理分析
嘌呤霉素处理对反应速度与底物浓度关系的影响
混合模型
发现问题,得参数初值
引入0-1变量
简化模型
检查参数置信区间是否包含零点
END