实验三-回归分析

实验三 回归分析
1.为了分析X 射线的杀菌作用，用200千伏的X 射线来照射细菌，每次照射6分钟用平板计数法估计尚存活的细菌数，照射次数记为 t ，照射后的细菌数y 如表1所示。

表1 X 射线照射次数与残留细菌数
试求：①给出y 与t 的二次函数回归模型；②在同一坐标系内做出原始数据与拟合结果的散点图；③预测t=16时残留的细菌数；④根据问题实际意义，你认为选择多项式函数是否合适？⑤给出非线性回归模型，并预测照射16次后细菌残留数目。

解：(1)实验程序： t=1:15;
y=[352 211 197 160 142 106 104 60 56 38 36 32 21 19 15];
rstool(t',y','purequadratic')
结果如图1所示：
4
6
8
10
12
50100150200250300350
图1
在Matlab 工作区中输入命令：beta ，rmse
beta =347.8967 -51.1394 1.9897 rmse =22.2649
所以y 与t 的二次回归模型函数：29897.11394.518967.347t t y +-= （2）画出同一坐标散点图，如图2所示，程序如下： [p,s]=polyfit(t,y,2); Y=polyconf(p,t,y); plot(t,y,'k+',t,Y,'r')
051015
50100150200250300350
400
图2 散点图
(3)当t=16时，计算程序如下： [p,s]=polyfit(t,y,2); Y=polyconf(p,16); 结果是：Y =39.0396
即说明预测残留的细菌数y=39.0396个；
(4)用二次函数计算出细菌残留数为39.0396，显然与实际不相符合。

根据实际问题的意义可知：尽管二次多项式拟合效果较好，但是用于预测并不理想。

因此，如何根据原始数据散点图的规律，选择适当的回归曲线是非常重要的，因此有必要研究非线性回归分析。

(5)由(2)散点图可知，可以假设将要拟合的的非线性模型为t /b ae y =对将要拟合的非线性模型t /b ae y =，建立的M-文件volum.m 如下： function yhat=volum(beta,t) yhat=beta(1)*exp(beta(2).*t);
%输入数据 t=1:15;
y=[352 211 197 160 142 106 104 60 56 38 36 32 21 19 15]; beta0=[150,0]';
%求回归系数
[beta,r,J]=nlinfit(t',y','volum',beta0); beta
y=nlpredci('volum',16,beta,r,J)
得结果：beta =400.0905 -0.2240，y =11.1014，即回归模型为：t e y 2240.0-0905.400=，那么根据此模型我们可以知道：当t=16时，残留的细菌数y=11.1014，很显然这样的结果会更令人满意！
2.某销售公司将库存占用资金情况、广告投入的费用、员工薪酬以与销售额等方面的数据作了汇总（表 2）,该公司试图根据这些数据找到销售额与其他变量之间的关系，以便进行销售额预测并为工作决策提供参考依据。

(1)建立销售额的回归模型；(2)如果未来某月库存资金额为150万元，广告投入预算为45万元，员工薪酬总额为27万元，试根据建立的回归模型预测该月的销售额。

表2 库存资金额、广告投入、员工薪酬、销售额汇总表（单位：万元）
月份 库存资金额(x1) 广告投入(x2) 员工薪酬总额(x3) 销售额(y)
1 75.
2 2 77.6
3 80.7
4 76.0
5 79.5
6 81.8
7 67.7
8 98.3
9 74.0 10 151.0 11 90.8 12 102.3 13 115.6 14 125.0 15 137.8 16 175.6 17 155.2 18 174.3
解：首先，作出因变量与各自变量的样本散点图，如图3所示，程序如下：
x1=[75.2 77.6 80.7 76.0 79.5 81.8 67.7 98.3 74.0 151.0 90.8 102.3 115.6 125.0 137.8 175.6 155.2 174.3];
x2=[30.6 31.3 33.9 29.6 32.5 27.9 24.8 23.6 33.9 27.7 45.5 42.6 40.0 45.8 51.7 67.2 65.0 65.4];
x3=[21.1 21.4 22.9 21.4 21.5 21.7 21.5 21.0 22.4 24.7 23.2 24.3 23.1 29.1 24.6 27.5 26.5 26.8];
y=[1090.4 1133.0 1242.1 1003.2 1283.2 1012.2 1098.8 826.3 1003.3 1554.6 1199.0 1483.1 1407.1 1551.3 1601.2 2311.7 2126.7 2256.5]; subplot(1,3,1),plot(x1,y,'g*'); subplot(1,3,2),plot(x2,y,'k+'); subplot(1,3,3),plot(x3,y,'ro');
21.1 21.4 22.9 21.4 21.5 21.7 21.5 21.0 22.4 24.7 23.2 24.3 23.1 29.1 24.6 27.5 26.5 26.8 1090.4 1133.0 1242.1 1003.2 1283.2 1012.2 1098.8 826.3 1003.3 1554.6 1199.0 1483.1 1407.1 1551.3 1601.2 2311.7 2126.7 2256.5
30.6 31.3 33.9 29.6 32.5 27.9 24.8 23.6 33.9 27.7 45.5 42.6 40.0 45.8 51.7 67.2 65.0 65.4
010*******
1000120014001600180020002200
2400050100800
10001200
14001600180020002200240020
2530
80010001200140016001800200022002400
图3 因变量y 与各自变量的样本散点图
从图上可以看出这些点大致分布在一条直线旁边，因此有较好的线性关系，
可以采用线性回归。

设回归方程为：3322110ˆˆˆˆˆx x x y ββββ+++=，建立M-文件输入如下程序：
x1=[75.2 77.6 80.7 76.0 79.5 81.8 67.7 98.3 74.0 151.0 90.8 102.3 115.6 125.0 137.8 175.6 155.2 174.3];
x2=[30.6 31.3 33.9 29.6 32.5 27.9 24.8 23.6 33.9 27.7 45.5 42.6 40.0 45.8 51.7 67.2 65.0 65.4];
x3=[21.1 21.4 22.9 21.4 21.5 21.7 21.5 21.0 22.4 24.7 23.2 24.3 23.1 29.1 24.6 27.5 26.5 26.8];
y=[1090.4 1133.0 1242.1 1003.2 1283.2 1012.2 1098.8 826.3 1003.3 1554.6 1199.0 1483.1 1407.1 1551.3 1601.2 2311.7 2126.7 2256.5]; n=18;m=3;
x=[ones(n,1),x1',x2',x3'];
[b,bint,r,rint,s]=regress(y',x,0.05); b,bint,r,rint,s
运行后即得到结果如表3所示
表3 对初步回归模型的计算结果
回归系数
回归系数的估计值
回归系数的置信区间 0β -53.9075 [-1011.2，903.4] 1β
5.7252 [2.0，9.4] 2β
15.2879 [6.4，24.2] 3β
9.5698
[-44.9，64.1]
12=R ,59=F ,0001.0<p ，172712=s
残差列向量r =[44.0394 59.3264 96.5750 -35.3239 179.3461 -36.4116 180.2202
-244.3408 -99.0812 84.1521 -184.5600 67.5082 -33.4048 -89.1104 -159.6274 69.7451 44.7425 56.2050]T
对应残差的)1(α-置信区间rint 如下：
[-228.8318，316.9105]、 [-214.8092，333.4620] [-173.5015，366.6515]、 [-311.2066，240.5589] [-75.9312，434.6233]、 [-313.1813， 240.3581] [-69.6357，430.0762]、 [-449.7576，-38.9240] [-365.7729，167.6105]、 [-69.8815，238.1857] [-428.0384，58.9185]、 [-208.3399，343.3563] [-312.3682， 245.5587]、[-199.0870，20.8662] [-415.9094，96.6547]、 [-172.6973，312.1875] [-207.2697，296.7547]、 [-186.8695，299.2794]。

因此得到初步的回归方程为：3215698.92879.157252.59075.53-ˆx x x y
+++=，当未来某月库存资金额为150万元，广告投入预算为45万元，员工薪酬总额为
27万元，那么根据所建立的回归模型可以预测出该月的销售额为1751.2万元。

3.葛洲坝机组发电耗水率的主要影响因素为库水位、出库流量。

现从数据库中将 2005年10月某天15时-16时06分范围内的出库流量、库水位对应的耗水率读取处理，数据如表4所示，试利用多元线性回归分析方法建立耗水率与出库流量、库水位的模型。

（表4数据来源：余波，多元线性回归分析在机组发电耗水率中的应用，计算机与现代化，2008（2））
表4 耗水率与出库流量、库水位数据
时间 库水位（米) 出库流量 机组发电耗水率 （年-月-天-时） （立方米） （立方米/万千瓦)
解：首先，作出耗水率y 与各自变量的样本散点图，如图4所示，程序如下：
x1=[65.08 65.10 65.12 65.17 65.21 65.37 65.38 65.39 65.40 65.43 65.47 65.53 65.62 65.58 65.70 65.84];
x2=[15607 15565 15540 15507 15432 15619 15536 15514 15519 15510 15489 15437 16355 14708 14393 14296];
y=[60.46 60.28 60.10 59.78 59.44 59.25 58.91 58.76 58.73 58.63 58.48 58.31 57.96 57.06 56.43 55.83]; subplot(1,2,1),plot(x1,y,'g*'); subplot(1,2,2),plot(x2,y,'k+');
65
65.5
66
55.5
5656.55757.55858.55959.56060.5 1.4
1.5
1.6
1.7x 10
4
55.5
5656.55757.55858.55959.56060.5
图4 耗水率y 与库水位、出库量关系散点图
65.08 65.10 65.12 65.17 65.21 65.37 65.38 65.39 65.40 65.43 65.47 65.53 65.62 65.58 65.70 65.84
2005-10-15:00 2005-10-15:02 2005-10-15:04 2005-10-15:06 2005-10-15:08 2005-10-15:10 2005-10-15:12 2005-10-15:14 2005-10-15:16 2005-10-15:18 2005-10-15:20 2005-10-15:22 2005-10-16:00 2005-10-16:02 2005-10-16:04 2005-10-16:06
15607 15565 15540 15507 15432 15619 15536 15514 15519 15510 15489 15437 16355 14708 14393 14296
60.46 60.28 60.10 59.78 59.44 59.25 58.91 58.76 58.73 58.63 58.48 58.31 57.96 57.06 56.43 55.83
从散点图中可以看出机组发电耗水率y 与库水位1x 有较好的线性关系，而与出库流量2x 的关系难以确定，可以采用建立二次函数的回归模型。

一般的多元二项式回归模型可表示为：εβ
βββ++
+++=∑≤≤m
k j k j jk
m m x x x x y ,1110...。

程序如下：
x=[x1',x2'];
rstool(x,y,'interaction'，0.05) % 包含线性项和完全二次项(交叉)
结果得到交互式化面，如图5所示
图5 耗水率y 与库水位、出库流量的一个交互界面
在左边图形下方的方框中输入65.40，右边图形下方的方框中输入15519，则画面左边的‘Predicted Y ’下方的数据变为58.7291,其置信区间在16162.07291.58±，即表明预测出库水位为65.40米，出库流量为15529立方米时的机组发电耗水率为58.7291立方米/万千瓦。

在画面左下方的下拉式菜单Export 里选“all ”，则beta 、rmse 和residuals 的值都发送到MATLAB 界面，只需在MATLAB 界面中输入：beta ，rmse ，residuals 即可得到结果如下：
beta =-559.9056 9.2385 0.0611 -0.0009 rmse= 0.1643
residuals =0.0079 -0.0189 -0.0674 -0.0963 -0.1520 0.2737 0.0638 -0.0149 0.0009 0.0605 0.1302 0.3036 -0.2672 -0.1403 -0.0481 -0.0356
故该多元线性回归模型：21210009.0-0611.02385.99056.559-x x x x y ++=，剩余标准差为0.1643，说明此回归模型的显著效果很好。

65.2
65.3
65.4
65.5
65.6
65.7
56
56.55757.55858.55959.56060.5 1.46
1.48
1.5
1.52
1.54
1.56
1.58
1.6
x 10
4。