北京工业大学-数学建模8-多元分析实验201312

数学建模作业8—多元分析实验一、基本实验1.回归分析为估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X(米)与当年灌溉面积Y(公顷)，测得连续10年的数据如表8.1所示。表8.110年中最大积雪深度与当年灌溉面积的数据XYXY15.1190767.8300023.5128774.5194737.1270085.6227346.2237398.0311358.83260106.42493（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.228Coefficients:EstimateStd.ErrortvaluePr(|t|)(Intercept)140.95125.111.1270.293x364.1819.2618.9086.33e-08***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:96.42on8degreesoffreedomMultipleR-squared:0.9781,AdjustedR-squared:0.9754F-statistic:357.5on1and8DF,p-value:6.33e-08系数检验：常数项截距Intercept的p值0.05,接受0假设，说明不显著。x的系数的p值0.05,拒绝0假设，x的系数极为显著，有***标示。方程检验：F-statistic的p值为6.33e-080.05,拒绝0假设，通过回归方程的检验。相关系数检验：R2=0.9781,调整的R2=0.9754,说明因变量与自变量是相关的。除了常数项截距，回归方程通过了回归参数的检验与回归方程的检验，得到的回归方程为：^140.95364.18YX（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)fitlwrupr12690.2272454.9712925.484predict(lm.sol,new,interval=confidence)fitlwrupr12690.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)Influencemeasuresoflm(formula=y~1+x):dfb.1_dfb.xdffitcov.rcook.dhatinf1-0.34940.2706-0.44791.1650.099410.1572-1.67001.5077-1.73200.8591.065030.413*30.0232-0.0475-0.10531.4610.006270.1264-0.02710.0056-0.08891.4230.004470.10050.5656-0.6935-0.82081.4440.326490.3496-0.04730.06630.09641.5940.005280.19071.2525-1.05771.40850.4440.580590.229*80.2329-0.15310.37861.1200.071170.1209-0.18840.25360.34651.4740.064570.215100.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.5lm.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.036Coefficients:EstimateStd.ErrortvaluePr(|t|)(Intercept)64.97118.670.5470.601x373.7517.4021.4851.19e-07***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:71.08on7degreesoffreedomMultipleR-squared:0.9851,AdjustedR-squared:0.9829F-statistic:461.6on1and7DF,p-value:1.192e-07与修正之前相比，R2与调整的R2有所提高。回归直线为：^64.97373.75YX检验：修正后是否有效，看一下回归诊断的结果，编写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某地区土壤所含可给态磷的情况X1X2X3YX1X2X3Y10.452158641012.6581125120.423163601110.9371117633.11937711223.1461149640.634157611323.1501347754.72459541421.644739361.765123771523.1561689579.4444681161.93614354810.131117931726.858202168911.629173931829.95112499（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