第十章 统计回归模型 10.1 牙膏的销售量 10.2 软件开发人员的薪金 10.3 酶促反应 10.4 投资额与国民生产总值和 物价指数
数学建模的基本方法 机理分析 测试分析 由于客观事物内部规律的复杂及人们认识程度的限制,无法分析实际对象内在的因果关系,建立合乎机理规律的数学模型。 通过对数据的统计分析,找出与数据拟合最好的模型 回归模型是用统计分析方法建立的最常用的一类模型 不涉及回归分析的数学原理和方法 通过实例讨论如何选择不同类型的模型 对软件得到的结果进行分析,对模型进行改进
10.1 牙膏的销售量 问题 建立牙膏销售量与价格、广告投入之间的模型 预测在不同价格和广告费用下的牙膏销售量 10.1 牙膏的销售量 建立牙膏销售量与价格、广告投入之间的模型 问题 预测在不同价格和广告费用下的牙膏销售量 收集了30个销售周期本公司牙膏销售量、价格、广告费用,及同期其它厂家同类牙膏的平均售价 9.26 0.55 6.80 4.25 3.70 30 7.93 0.05 5.80 3.85 3.80 29 8.51 0.25 6.75 4.00 3.75 2 7.38 -0.05 5.50 1 销售量 (百万支) 价格差 (元) 广告费用 (百万元) 其它厂家价格(元) 本公司价格(元) 销售周期
~随机误差(均值为零的正态分布随机变量) 基本模型 x1 y y ~公司牙膏销售量 x1~其它厂家与本公司价格差 x2~公司广告费用 y~被解释变量(因变量) x2 y x1, x2~解释变量(回归变量, 自变量) 0, 1 , 2 , 3 ~回归系数 ~随机误差(均值为零的正态分布随机变量)
模型求解 由数据 y,x1,x2估计 [b,bint,r,rint,stats]=regress(y,x,alpha) 输入 输出 MATLAB 统计工具箱 由数据 y,x1,x2估计 [b,bint,r,rint,stats]=regress(y,x,alpha) 输入 输出 y~n维数据向量 b~的估计值 x= ~n4数据矩阵, 第1列为全1向量 bint~b的置信区间 r ~残差向量y-xb alpha(置信水平,0.05) rint~r的置信区间 参数 参数估计值 置信区间 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 Stats~ 检验统计量 R2,F, p
结果分析 y的90.54%可由模型确定 F远超过F检验的临界值 p远小于=0.05 模型从整体上看成立 x2对因变量y 的影响不太显著 参数 参数估计值 置信区间 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%可由模型确定 F远超过F检验的临界值 p远小于=0.05 模型从整体上看成立 2的置信区间包含零点(右端点距零点很近) x2对因变量y 的影响不太显著 x22项显著 可将x2保留在模型中
若估计x3=3.9,设定x4=3.7,则可以95%的把握知道销售额在 7.83203.7 29(百万元)以上 销售量预测 价格差x1=其它厂家价格x3-本公司价格x4 估计x3 调整x4 控制x1 通过x1, x2预测y 控制价格差x1=0.2元,投入广告费x2=650万元 (百万支) 销售量预测区间为 [7.8230,8.7636](置信度95%) 上限用作库存管理的目标值 下限用来把握公司的现金流 若估计x3=3.9,设定x4=3.7,则可以95%的把握知道销售额在 7.83203.7 29(百万元)以上
模型改进 x1和x2对y的影响独立 x1和x2对y的影响有交互作用 参数 参数估计值 置信区间 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 x1和x2对y的影响独立 x1和x2对y的影响有交互作用 参数 参数估计值 置信区间 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=0.2元,投入广告费x2=6.5百万元 (百万支) 区间 [7.8230,8.7636] (百万支) 区间 [7.8953,8.7592] 略有增加 预测区间长度更短
两模型 与x1,x2关系的比较 x1 x1 x2=6.5 x2 x2 x1=0.2
交互作用影响的讨论 价格差 x1=0.1 价格差 x1=0.3 价格优势会使销售量增加 x2 加大广告投入使销售量增加 价格差较小时增加的速率更大 价格差较小时更需要靠广告来吸引顾客的眼球
MATLAB中有命令rstool直接求解 完全二次多项式模型 MATLAB中有命令rstool直接求解 x1 x2 从输出 Export 可得
10.2 软件开发人员的薪金 建立模型研究薪金与资历、管理责任、教育程度的关系 分析人事策略的合理性,作为新聘用人员薪金的参考 10.2 软件开发人员的薪金 建立模型研究薪金与资历、管理责任、教育程度的关系 分析人事策略的合理性,作为新聘用人员薪金的参考 编号 薪金 资历 管理 教育 01 13876 1 02 11608 3 03 18701 04 11283 2 42 27837 16 43 18838 44 17483 45 19207 17 46 19346 20 46名软件开发人员的档案资料 资历~ 从事专业工作的年数;管理~ 1=管理人员,0=非管理人员;教育~ 1=中学,2=大学,3=更高程度
分析与假设 y~ 薪金,x1 ~资历(年) x2 = 1~ 管理人员,x2 = 0~ 非管理人员 1=中学2=大学3=更高 中学:x3=1, x4=0 ;大学:x3=0, x4=1; 更高:x3=0, x4=0 教育 资历每加一年薪金的增长是常数; 管理、教育、资历之间无交互作用 线性回归模型 a0, a1, …, a4是待估计的回归系数,是随机误差
模型求解 参数 参数估计值 置信区间 a0 11032 [ 10258 11807 ] a1 546 [ 484 608 ] a2 6883 [ 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 资历增加1年薪金增长546 管理人员薪金多6883 中学程度薪金比更高的少2994 大学程度薪金比更高的多148 R2,F, p 模型整体上可用 中学:x3=1, x4=0;大学:x3=0, x4=1; 更高:x3=0, x4=0. x2 = 1~ 管理,x2 = 0~ 非管理 x1~资历(年) a4置信区间包含零点,解释不可靠!
结果分析 残差分析方法 残差 组合 1 2 3 4 5 6 管理 教育 管理与教育的组合 e 与资历x1的关系 e与管理—教育组合的关系 教育 管理与教育的组合 残差 e 与资历x1的关系 e与管理—教育组合的关系 残差全为正,或全为负,管理—教育组合处理不当 残差大概分成3个水平, 6种管理—教育组合混在一起,未正确反映 。 应在模型中增加管理x2与教育x3, x4的交互项
进一步的模型 增加管理x2与教育x3, x4的交互项 参数 参数估计值 置信区间 a0 11204 [11044 11363] a1 497 [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 e ~ x1 e ~组合 R2,F有改进,所有回归系数置信区间都不含零点,模型完全可用 消除了不正常现象 异常数据(33号)应去掉
去掉异常数据后的结果 R2: 0.957 0.999 0.9998 残差图十分正常 F: 226 554 36701 e ~ x1 参数 参数估计值 置信区间 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 e ~组合 R2: 0.957 0.999 0.9998 F: 226 554 36701 置信区间长度更短 残差图十分正常 最终模型的结果可以应用
模型应用 制订6种管理—教育组合人员的“基础”薪金(资历为0) 大学程度管理人员比更高程度管理人员的薪金高 x1= 0; x2 = 1~ 管理,x2 = 0~ 非管理 中学:x3=1, x4=0 ;大学:x3=0, x4=1; 更高:x3=0, x4=0 组合 管理 教育 系数 “基础”薪金 1 a0+a3 9463 2 a0+a2+a3+a5 13448 3 a0+a4 10844 4 a0+a2+a4+a6 19882 5 a0 11200 6 a0+a2 18241 大学程度管理人员比更高程度管理人员的薪金高 大学程度非管理人员比更高程度非管理人员的薪金略低
软件开发人员的薪金 对定性因素(如管理、教育),可以引入0-1变量处理,0-1变量的个数应比定性因素的水平少1 残差分析方法可以发现模型的缺陷,引入交互作用项常常能够改善模型 剔除异常数据,有助于得到更好的结果 注:可以直接对6种管理—教育组合引入5个0-1变量
10.3 酶促反应 问题 方案 研究酶促反应(酶催化反应)中嘌呤霉素对反应速度与底物(反应物)浓度之间关系的影响 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 /
线性化模型 对1 , 2非线性 对1, 2线性 经嘌呤霉素处理后实验数据的估计结果 参数 参数估计值(×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
参数估计时,x较小(1/x很大)的数据控制了回归参数的确定 线性化模型结果分析 1/y 1/x x y 1/x较小时有很好的线性趋势,1/x较大时出现很大的起落 x较大时,y有较大偏差 参数估计时,x较小(1/x很大)的数据控制了回归参数的确定
[beta,R,J] = nlinfit (x,y,’model’,beta0) 非线性模型参数估计 MATLAB 统计工具箱 [beta,R,J] = nlinfit (x,y,’model’,beta0) 输入 x~自变量数据矩阵 y ~因变量数据向量 beta0~线性化模型估计结果 x= ; y= ; beta0=[195.8027 0.04841]; [beta,R,J]=nlinfit(x,y,’f1’,beta0); betaci=nlparci(beta,R,J); beta, betaci model ~模型的函数M文件名 beta0 ~给定的参数初值 输出 beta ~参数的估计值R ~残差,J ~估计预测误差的Jacobi矩阵 beta的置信区间 function y=f1(beta, x) y=beta(1)*x./(beta(2)+x); betaci =nlparci(beta,R,J)
非线性模型结果分析 其它输出 最终反应速度为 半速度点(达到最终速度一半时的x值 )为 命令nlintool 给出交互画面 参数 参数估计值 置信区间 1 212.6819 [197.2029 228.1609] 2 0.0641 [0.0457 0.0826 ] o ~原始数据 + ~ 拟合结果 最终反应速度为 半速度点(达到最终速度一半时的x值 )为 其它输出 命令nlintool 给出交互画面 拖动画面的十字线,得 y的预测值和预测区间 画面左下方的Export 输出其它统计结果。 剩余标准差s= 10.9337
混合反应模型 在同一模型中考虑嘌呤霉素处理的影响 x1为底物浓度, x2为一示性变量 x2=1表示经过处理,x2=0表示未经处理 β1是未经处理的最终反应速度 γ1是经处理后最终反应速度的增长值 β2是未经处理的反应的半速度点 γ2是经处理后反应的半速度点的增长值
混合模型求解 参数初值 估计结果和预测 2置信区间包含零点,表明2对因变量y的影响不显著 经嘌呤霉素处理的作用不影响半速度点参数 用nlinfit 和 nlintool命令 参数初值 (基于对数据的分析) o ~原始数据 + ~拟合结果 估计结果和预测 未经处理 经处理 参数 参数估计值 置信区间 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] 剩余标准差s= 10.4000 2置信区间包含零点,表明2对因变量y的影响不显著 经嘌呤霉素处理的作用不影响半速度点参数
简化的混合模型 估计结果和预测 简化的混合模型形式简单,参数置信区间不含零点 剩余标准差 s = 10.5851,比一般混合模型略大 参数 o ~原始数据 + ~拟合结果 未经处理 经处理 参数 参数估计值 置信区间 1 166.6025 [154.4886 178.7164] 2 0.0580 [0.0456 0.0703 ] 1 42.0252 [28.9419 55.1085] 简化的混合模型形式简单,参数置信区间不含零点 剩余标准差 s = 10.5851,比一般混合模型略大
一般混合模型与简化混合模型预测比较 简化混合模型的预测区间较短,更为实用、有效 预测区间为预测值 Δ 实际值 一般模型预测值 Δ(一般模型) 简化模型预测值 Δ(简化模型) 67 47.3443 9.2078 42.7358 5.4446 51 84 89.2856 9.5710 84.7356 7.0478 … 191 190.8329 9.1484 189.0574 8.8438 201 207 200.9688 11.0447 198.1837 10.1812 200 简化混合模型的预测区间较短,更为实用、有效
嘌呤霉素处理对反应速度与底物浓度关系的影响 酶促反应 机理分析 反应速度与底物浓度的关系 非线性关系 求解线性模型 发现问题,得参数初值 求解非线性模型 嘌呤霉素处理对反应速度与底物浓度关系的影响 引入0-1变量 混合模型 检查参数置信区间是否包含零点 简化模型 注:非线性模型拟合程度的评价无法直接利用线性模型的方法,但R2 与s仍然有效。
根据对未来GNP及PI的估计,预测未来投资额 10.4 投资额与国民生产总值和物价指数 问题 建立投资额模型,研究某地区实际投资额与国民生产总值 ( GNP ) 及物价指数 ( PI ) 的关系 根据对未来GNP及PI的估计,预测未来投资额 该地区连续20年的统计数据 年份序号 投资额 国民生产总值 物价 指数 年份 序号 投资额 国民生产总值 物价 指数 1 90.9 596.7 0.7167 11 229.8 1326.4 1.0575 2 97.4 637.7 0.7277 12 228.7 1434.2 1.1508 3 113.5 691.1 0.7436 13 206.1 1549.2 1.2579 4 125.7 756.0 0.7676 14 257.9 1718.0 1.3234 5 122.8 799.0 0.7906 15 324.1 1918.3 1.4005 6 133.3 873.4 0.8254 16 386.6 2163.9 1.5042 7 149.3 944.0 0.8679 17 423.0 2417.8 1.6342 8 144.2 992.7 0.9145 18 401.9 2631.7 1.7842 9 166.4 1077.6 0.9601 19 474.9 2954.7 1.9514 10 195.0 1185.9 1.0000 20 424.5 3073.0 2.0688
若采用普通回归模型直接处理,将会出现不良后果 投资额与国民生产总值和物价指数 分析 许多经济数据在时间上有一定的滞后性 以时间为序的数据,称为时间序列 时间序列中同一变量的顺序观测值之间存在自相关 若采用普通回归模型直接处理,将会出现不良后果 需要诊断并消除数据的自相关性,建立新的模型 … 1.3234 1718.0 257.9 14 0.7676 756.0 125.7 4 1.2579 1549.2 206.1 13 0.7436 691.1 113.5 3 1.1508 1434.2 228.7 12 0.7277 637.7 97.4 2 1.0575 1326.4 229.8 11 0.7167 596.7 90.9 1 物价 指数 国民生产总值 投资额 年份 序号 年份序号
基本回归模型 t ~年份, yt ~ 投资额,x1t~ GNP, x2t ~ 物价指数 投资额与 GNP及物价指数间均有很强的线性关系 0, 1, 2 ~回归系数 t ~对t相互独立的零均值正态随机变量
基本回归模型的结果与分析 剩余标准差 s=12.7164 模型优点 R2=0.9908,拟合度高 模型缺点 没有考虑时间序列数据的滞后性影响 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,拟合度高 模型缺点 没有考虑时间序列数据的滞后性影响 可能忽视了随机误差存在自相关;如果存在自相关性,用此模型会有不良后果
自相关性的定性诊断 模型残差 et为随机误差t 的估计值 在MATLAB工作区中输出 作残差 et~et-1 散点图 残差诊断法 模型残差 et-1 et et为随机误差t 的估计值 在MATLAB工作区中输出 作残差 et~et-1 散点图 大部分点落在第1, 3象限 t 存在正的自相关 大部分点落在第2, 4象限 t 存在负的自相关 基本回归模型的随机误差项t 存在正的自相关 自相关性直观判断
自回归性的定量诊断 D-W检验 自回归模型 0, 1, 2 ~回归系数 ρ~自相关系数 ut ~对t相互独立的零均值正态随机变量 ρ= 0 无自相关性 ρ> 0 存在正自相关性 存在负自相关性 ρ< 0 如何估计ρ D-W统计量 如何消除自相关性 广义差分法
D-W统计量与D-W检验 检验水平,样本容量,回归变量数目 由DW值的大小确定自相关性 检验临界值dL和dU n较大 DW 4-dU 4 2 正自 相关 负自 不能确定 无自相关 检验水平,样本容量,回归变量数目 D-W分布表 由DW值的大小确定自相关性 检验临界值dL和dU
广义差分变换 原模型 变换 新模型 以*0, 1 , 2 为回归系数的普通回归模型 步骤 原模型 DW值 D-W检验 无自相关 原模型 继续此过程 有自相关 新模型 增加数据量;选用其它方法 不能确定
投资额新模型的建立 原模型残差et DWold=0.8754 样本容量n=20,回归变量数目k=3,=0.05 原模型有正自相关 DWold < dL 查表 临界值dL=1.10, dU=1.54 作变换 DW 4-dU 4 4-dL dU dL 2 正自 相关 负自 不能确定 无自相关
投资额新模型的建立 总体效果良好 剩余标准差 snew= 9.8277 < sold=12.7164 参数 参数估计值 置信区间 *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 总体效果良好 剩余标准差 snew= 9.8277 < sold=12.7164
新模型的自相关性检验 新模型残差et DWnew=1.5751 样本容量n=19,回归变量数目k=3,=0.05 4-dU 4 4-dL dU dL 2 正自 相关 负自 不能确定 无自相关 新模型残差et DWnew=1.5751 样本容量n=19,回归变量数目k=3,=0.05 dU< DWnew < 4-dU 新模型无自相关性 查表 临界值dL=1.08, dU=1.53 新模型 还原为 原始变量 一阶自回归模型
模型结果比较 基本回归模型 一阶自回归模型 一阶自回归模型残差et比基本回归模型要小 新模型 et~ *,原模型 et~ + 残差图比较 新模型 ŷt ~ *,新模型 ŷt ~ + 拟合图比较 一阶自回归模型残差et比基本回归模型要小
投资额预测 对未来投资额yt 作预测,需先估计出未来的国民生产总值x1t 和物价指数 x2t 一阶自回归模型 2.0688 3073.0 424.5 20 1.9514 2954.7 474.9 19 1.7842 2631.7 401.9 18 0.7436 691.1 113.5 3 0.7277 637.7 97.4 2 0.7167 596.7 90.9 1 物价 指数 国民生产总值 投资额 年份 序号 年份序号 设已知 t=21时, x1t =3312,x2t=2.1938 基本回归模型 一阶自回归模型 ŷt 较小是由于yt-1=424.5过小所致