北京工业大学数学建模8多元分析实验12.docx
《北京工业大学数学建模8多元分析实验12.docx》由会员分享，可在线阅读，更多相关《北京工业大学数学建模8多元分析实验12.docx（22页珍藏版）》请在冰点文库上搜索。
北京工业大学数学建模8多元分析实验12
数学建模作业8—多元分析实验
一、基本实验
1.回归分析
为估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X（米）与当年灌溉面积Y（公顷），测得连续10年的数据如表8.1所示。
表8.110年中最大积雪深度与当年灌溉面积的数据
X
Y
X
Y
1
5.1
1907
6
7.8
3000
2
3.5
1287
7
4.5
1947
3
7.1
2700
8
5.6
2273
4
6.2
2373
9
8.0
3113
5
8.8
3260
10
6.4
2493
（1）建立一元线性回归模型，求解，并验证系数、方程或相关系数是否通过检验；
（2）观测得今年的数据是X=7米，给出今年灌溉面积的预测值、预测区间和置信区间（α=0.05）；
（3）将数据散点、回归预测值、回归的预测区间和置信区间均画在一张图上，分析线性回归的拟合情况。
解：
（1）利用R软件中的lm（）求回归参数并作相应的检验。
写出相应的R程序，如下：
x<-c（5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4）
y<-c（1907,1287,2700,2373,3260,3000,1947,2273,3113,2493）
lm.sol<-lm（y~1+x）
summary（lm.sol）
运行结果如下：
Call:
lm（formula=y~1+x）
Residuals:
Min1QMedian3QMax
-128.591-70.978-3.72749.263167.228
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）140.95125.111.1270.293
x364.1819.2618.9086.33e-08***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
96.42on8degreesoffreedom
MultipleR-squared:
0.9781,AdjustedR-squared:
0.9754
F-statistic:
357.5on1and8DF,p-value:
6.33e-08
系数检验：
常数项截距Intercept的p值>0.05,接受0假设，说明不显著。
x的系数的p值<0.05,拒绝0假设，x的系数极为显著，有***标示。
方程检验：
F-statistic的p值为6.33e-08<0.05,拒绝0假设，通过回归方程的检验。
相关系数检验：
R2=0.9781,调整的R2=0.9754,说明因变量与自变量是相关的。
除了常数项截距，回归方程通过了回归参数的检验与回归方程的检验，得到的回归方程为：
（2）利用predict（）函数进行预测，编写R程序：
new<-data.frame（x=7）
predict（lm.sol,new,interval="prediction",level=0.95）
predict（lm.sol,new,interval="confidence"）
运行结果分别为：
predict（lm.sol,new,interval="prediction",level=0.95）
fitlwrupr
12690.2272454.9712925.484
>predict（lm.sol,new,interval="confidence"）
fitlwrupr
12690.2272613.352767.105
有分析结果可知：
今年灌溉渠的预测值为2690.227公顷；
预测区间是[2613.35,2767.105]（单位：
公顷）；
置信区间（α=0.05）为[2454.971,2925.484]（单位：
公顷）。
（3）参考exam0806.R,编写R程序：
x<-c（5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4）
y<-c（1907,1287,2700,2373,3260,3000,1947,2273,3113,2493）
lm.sol<-lm（y~1+x）
new<-data.frame（x=seq（3.5,8.8,by=0.1））
pp<-predict（lm.sol,new,interval="prediction"）
pc<-predict（lm.sol,new,interval="confidence"）
matplot（new$x,cbind（pp,pc[,-1]）,type="l",
xlab="X",ylab="Y",lty=c（1,5,5,2,2）,
col=c（"blue","red","red","green","green"）,
lwd=2）
points（x,y,cex=1.4,pch=21,col="red",bg="orange"）
legend（4,3500,
c（"Points","Fitted","Prediction","Confidence"）,
pch=c（19,NA,NA,NA）,lty=c（NA,1,5,2）,
col=c（"orange","blue","red","brown"））
savePlot（"predict",type="eps"）
运行结果如图8.1：
图中显示散点、回归预测值、回归的预测区间和置信区间。
散点紧密的散布在回归直线两侧，只有一点落在置信区间外侧。
散点分布与回归直线趋势一致，总体效果还可以。
图8.1数据的散点、回归预测值、预测区间与置信区间
2.回归诊断
对1题得到的回归模型做回归诊断：
（1）残差是否满足齐性、正态性的条件；
（2）那些点可能是异常值点；
（3）如果有异常值点，则去掉异常值点，再作回归分析。
解：
利用函数plot（）绘制各种残差的散点图，编写程序如下：
调用influence.measures（）函数并作回归诊断图，R程序如下：
influence.measures（lm.sol）
op<-par（mfrow=c（2,2）,mar=0.4+c（2,2,1,1）,
oma=c（2,2,0,0））
plot（lm.sol,1:
4）
par（op）
Influencemeasuresof
lm（formula=y~1+x）:
dfb.1_dfb.xdffitcov.rcook.dhatinf
1-0.34940.2706-0.44791.1650.099410.157
2-1.67001.5077-1.73200.8591.065030.413*
30.0232-0.0475-0.10531.4610.006270.126
4-0.02710.0056-0.08891.4230.004470.100
50.5656-0.6935-0.82081.4440.326490.349
6-0.04730.06630.09641.5940.005280.190
71.2525-1.05771.40850.4440.580590.229*
80.2329-0.15310.37861.1200.071170.120
9-0.18840.25360.34651.4740.064570.215
100.01330.00460.07301.4320.003020.100
图8.2回归诊断图
先进性回归诊断结果分析，得到的回归诊断结果共7列，第1,2列为dfbetas值（对应于常数和变量x）;第3列为DFFITS准则值；第4列为COVRATIO准则值；第5列为Cook距离；第6列为帽子值（高杠杆值）；第7列为影像店记号。
有诊断结果可知2号和7号点是强影响点。
下面分析回归诊断图8.2：
第一张图是残差图，散点基本在联系周围，可以认为残差的方差满足齐性。
第二张图是正态QQ图，除7号点外，基本都在一条直线上，也就是说除了7号点外，残差满足正态性。
第三张图是标准差的平方根与预测值的散点图，7号值最大接近1.4，可能是异常值点。
第四张图是Cook距离值，从图上来看2号距离最大，说明2号点可能是强影响点（高杠杆点）。
综上：
（1）残差满足齐性和正态性条件。
（2）7号点可能是异常值点，2号点是强影响点
（3）将2号点和7号点去掉，再作回归分析。
处理强影响点：
将7号点去掉，2号点赋值0.5的权重，加权计算减少它的影响。
编写R程序如下：
intellect<-data.frame（
x=c（5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4）,
y=c（1907,1287,2700,2373,3260,3000,1947,2273,3113,2493）
）
lm.sol<-lm（y~1+x,data=intellect）
summary（lm.sol）
n<-length（intellect$x）
weights<-rep（1,n）;weights[2]<-0.5
lm.correct<-lm（y~1+x,data=intellect,subset=-7,
weights=weights）
summary（lm.correct）
运行结果为：
Call:
lm（formula=y~1+x,data=intellect,subset=-7,weights=weights）
WeightedResiduals:
Min1QMedian3QMax
-93.962-60.875-9.21336.037115.036
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）64.97118.670.5470.601
x373.7517.4021.4851.19e-07***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
71.08on7degreesoffreedom
MultipleR-squared:
0.9851,AdjustedR-squared:
0.9829
F-statistic:
461.6on1and7DF,p-value:
1.192e-07
与修正之前相比，R2与调整的R2有所提高。
回归直线为：
检验：
修正后是否有效，看一下回归诊断的结果，编写R程序：
op<-par（mfrow=c（2,2）,mar=0.4+c（4,4,1,1）,oma=c（0,0,2,0））
plot（lm.correct,1:
4）
par（op）
运行结果见图8.3修正后的回归诊断图，从四张图中可见残差齐次性、正态性，标准差的平方根与预测值，以及Cook距离等结果均有所改善。
图8.3修正后的回归诊断图
3.回归分析和逐步回归
研究同一地区土壤所含可给态磷（Y）的情况，得到18组数据如表8.2所示。
表中X1为土壤内所含无机磷浓度，X2为土壤内溶于K2CO3溶液并受溴化物水解的有机磷，X3为土壤内溶于K2CO3溶液但不溶于溴化物水解的有机磷。
表8.2某地区土壤所含可给态磷的情况
X1
X2
X3
Y
X1
X2
X3
Y
1
0.4
52
158
64
10
12.6
58
112
51
2
0.4
23
163
60
11
10.9
37
111
76
3
3.1
19
37
71
12
23.1
46
114
96
4
0.6
34
157
61
13
23.1
50
134
77
5
4.7
24
59
54
14
21.6
44
73
93
6
1.7
65
123
77
15
23.1
56
168
95
7
9.4
44
46
81
16
1.9
36
143
54
8
10.1
31
117
93
17
26.8
58
202
168
9
11.6
29
173
93
18
29.9
51
124
99
（1）建立多元线性回归方程模型，求解，并验证系数、方程或相关系数是否通过检验；
（2）作逐步回归分析。
解：
（1）输入数据，做多元线性回归：
cement<-data.frame（
X1=c（0.4,0.4,3.1,0.6,4.7,1.7,9.4,10.1,11.6,12.6,10.9,23.1,23.1,21.6,23.1,1.9,26.8,29.9）,
X2=c（52,23,19,34,24,65,44,31,29,58,37,46,50,44,56,36,58,51）,
X3=c（158,163,37,157,59,123,46,117,173,112,111,114,134,73,168,143,202,124）,
Y=c（64,60,71,61,54,77,81,93,93,51,76,96,77,93,95,54,168,99）
）
lm.sol<-lm（Y~X1+X2+X3,data=cement）
summary（lm.sol）
运行结果为：
Call:
lm（formula=Y~X1+X2+X3,data=cement）
Residuals:
Min1QMedian3QMax
-28.349-11.383-2.65912.09548.807
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）43.6500718.054422.4180.02984*
X11.785340.539773.3080.00518**
X2-0.083290.42037-0.1980.84579
X30.161020.111581.4430.17098
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
19.97on14degreesoffreedom
MultipleR-squared:
0.5493,AdjustedR-squared:
0.4527
F-statistic:
5.688on3and14DF,p-value:
0.009227
由运行结果可知：
常数项一般相关，X1系数高度相关，但是X2,X3系数的P值>0.05,不相关，没有通过检验。
方程F检验P值<0.05,通过检验。
相关系数R2,及调整R2在0.5左右，相关性一般。
（2）下面用函数step（）做逐步回归：
lm.step<-step（lm.sol）
运行结果为：
Start:
AIC=111.27
Y~X1+X2+X3
DfSumofSqRSSAIC
-X2115.75599.4109.32
5583.7111.27
-X31830.66414.4111.77
-X114363.49947.2119.66
Step:
AIC=109.32
Y~X1+X3
DfSumofSqRSSAIC
5599.4109.32
-X31833.26432.6109.82
-X115169.510768.9119.09
从程序的运行结果可见，当用全部变量作回归方程时，AIC值为111.27，接下来的数据表显示，如果去掉X2,则相应的ACI值为109.32，此时ACI值达到最小，因此R软件自动去掉变量X2,进行下一轮计算。
R软件认为此时得到“最优”回归方程。
下面我们分析一下回归结果，用函数summary（）提取相关信息，看看是否是最优的回归方程:
summary（lm.step）
Call:
lm（formula=Y~X1+X3,data=cement）
Residuals:
Min1QMedian3QMax
-29.713-11.324-2.95311.28648.679
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）41.479413.88342.9880.00920**
X11.73740.46693.7210.00205**
X30.15480.10361.4940.15592
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
19.32on15degreesoffreedom
MultipleR-squared:
0.5481,AdjustedR-squared:
0.4878
F-statistic:
9.095on2and15DF,p-value:
0.002589
有结果可见，回归系数检验的显著性水平有提高，但是X3系数的检验水平仍不理想。
下面用drop1来做逐步回归：
drop1（lm.step）
Singletermdeletions
Model:
Y~X1+X3
DfSumofSqRSSAIC
5599.4109.32
X115169.510768.9119.09
X31833.26432.6109.82
从运算结果可见如果去掉X3,则AIC值从109,32增加到109.82，增加的最少。
去掉X3,残差的平方和上升的也最少，拟合越好的方程，残差的平方和应越小，从这两项指标，应该再去掉变量X3.
lm.opt<-lm（Y~1+X1,data=cement）;
summary（lm.opt）
运行结果为：
Call:
lm（formula=Y~1+X1,data=cement）
Residuals:
Min1QMedian3QMax
-31.486-8.282-1.6745.62359.337
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）59.25907.42007.9865.67e-07***
X11.84340.47893.8490.00142**
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
20.05on16degreesoffreedom
MultipleR-squared:
0.4808,AdjustedR-squared:
0.4484
F-statistic:
14.82on1and16DF,p-value:
0.001417
最终结果还是满意的，尽管R2还是0.5左右。
但是系数检验的显著性水平已经很好了。
最后得到的“最优”回归方程为：
4.方差分析I（单因素方差分析）
24只小鼠随机地分成3组，每组8只，每组喂养不同的饲料组，喂养一定时间后，测得小鼠肝中含铁量，结果如表8.3所示。
表8.3不同饲料组小鼠肝脏中铁含量（单位：
ug/g）
饲料
1
2
3
4
5
6
7
8
A
1.00
1.01
1.13
1.14
1.70
2.01
2.23
2.63
B
0.96
1.23
1.54
1.96
2.94
3.68
5.59
6.96
C
2.07
3.72
4.50
4.90
6.00
6.84
8.23
10.33
（1）试用单因素方差分析方法分析不同饲料的小鼠肝中的铁含量是否有显著差异？
（2）如果有显著差异，那些水平之间有显著差异？
（3）数据是否满足正态性和方差齐性的要求，如果不满足请选择合适的分析方法，并给出你的最终结论。
解：
（1）利用方差分析表的方法进行分析：
设小白鼠喂养的饲料为因素，三种不同的饲料为三个水平，喂养饲料后的小白鼠肝脏中铁含量视为来自三个正态分布总体
（i=1,2,3）的样本观测值。
问题归结为：
。
参考exam0813.R编写R程序：
mouse<-data.frame（
X=c（1.00,1.01,1.13,1.14,1.70,2.01,2.23,2.63,0.96,1.23,1.54,1.96,2.94,3.68,5.59,6.96,2.07,3.72,4.50,4.90,6.00,6.84,8.23,10.33）,
A=factor（rep（1:
3,c（8,8,8）））
）
mouse.lm<-lm（X~A,data=mouse）
anova（mouse.lm）
运行结果为：
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr（>F）
A273.11836.5599.1040.001422**
Residuals2184.3294.016
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
由
升级会员 