数学模型实验8.docx

数学模型实验8.docx

《数学模型实验8.docx》由会员分享，可在线阅读，更多相关《数学模型实验8.docx（17页珍藏版）》请在冰豆网上搜索。

数学模型实验8

数学建模实验八多元分析实验

1回归分析

解：

（1）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）

summary（lm.sol,correlation=T,symbolic.cor=F）

得到的结果如下：

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

由计算结果知：

系数

不显著，

极为显著。

相关系数和方程通过检验。

得到的一元线性回归模型为：

（2）R程序如下：

new<-data.frame（x=7）

predict（lm.sol,new,interval="prediction",level=0.95）

得到的结果如下：

fitlwrupr

12690.2272454.9712925.484

预测值为2690.227,区间为 （2454.971,2925.484）

（3）输入程序如下：

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","brown","brown"）,

lwd=2）

points（x,y,cex=1.4,pch=21,col="red",bg="orange"）

legend（4.1,3000,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"）

得到的结果如下：

分析：

从图中看出所有的点都落在预测区间内。

只有第7个点不在置信区间内，其余都在置信区间内。

可以看出拟合情况较好。

2回归分析

解：

（1）输入程序如下：

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）

influence.measures（lm.sol）

op<-par（mfrow=c（2,2）,mar=0.4+c（4,4,1,1）,oma=c（0,0,2,0））

plot（lm.sol,1:

4）

par（op）

savePlot（"diagnoses1",type="eps"）

运行结果如下：

Call:

lm（formula=y~1+x,data=intellect）

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

>

>influence.measures（lm.sol）

Influencemeasuresof

lm（formula=y~1+x,data=intellect）:

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

计算结果说明：

"*"表示强影响点，所以2,7点为强影响点。

再分析回归诊断图：

这里共四张图，第1张是残差图，可以看出残差的方差不满足齐性。

第二张是正态QQ图，除7号点外，残差满足正态性。

第三张是标准差的平方根与预测值的散点图，到小于1.5。

第四张给出了Cook距离值，从图上来看，2号点的Cook距离最大。

（2）这说明2号点可能是强影响点。

由于7号点是异常值点，将在后面的计算中剔除，2号点是强影响点，加权计算减少它的影响。

 （3）运行程序如下：

n<-length（intellect$x）

weights<-rep（1,n）;weights[2]<-c（0.5）

lm.correct<-lm（y~1+x,data=intellect,subset=-2,

weights=weights）

summary（lm.correct）

attach（intellect）

plot（x,y,cex=1.2,pch=21,col="red",bg="orange"）

abline（lm.sol,col="blue",lwd=2）

text（x[c（19,18）],y[c（19,18）],

label=c（"19","18"）,adj=c（1.5,0.3））

detach（）

abline（lm.correct,col="red",lwd=2,lty=5）

legend（30,120,c（"Points","Regression","CorrectReg"）,

pch=c（19,NA,NA）,lty=c（NA,1,5）,

col=c（"orange","blue","red"））

savePlot（"diagnoses2",type="eps"）

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）

运行结果如下：

Call:

lm（formula=y~1+x,data=intellect,subset=-7,weights=weights）

Residuals:

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

结果说明：

从图中看出，7号点的残差过大，2号点对整个回归直线有较大影响。

这说明前面的回归诊断是合理的。

图中实线是原始数据计算的结果，虚线是修正数据后的计算结果。

经过检验，得到的图形如下：

从图中看出，所得的结果均有所改善。

3回归分析和逐步回归

解:

（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）

lm.step<-step（lm.sol）

summary（lm.step）

得到的结果如下：

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

得到的结果表明：

显著，

高度显著，

、

不显著。

并且没有通过F检验。

（2）所以需要step（）作逐步回归。

计算结果如下：

>

>lm.step<-step（lm.sol）

Start:

AIC=111.27

Y~X1+X2+X3

DfSumofSqRSSAIC

-X2115.75599.4109.3

5583.7111.3

-X31830.66414.4111.8

-X114363.49947.2119.7

Step:

AIC=109.32

Y~X1+X3

DfSumofSqRSSAIC

5599.4109.3

-X31833.26432.6109.8

-X115169.510768.9119.1

>

>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

结果说明：

从程序的运行结果可以看到，用全部变量做回归方程时，AIC是为111.27，接下来显示的数据表告诉我们，如果去掉X2，则相应的AIC值为109.3；如果去掉变量X3，则相应的AIC值为111.8。

后面的类推。

由于去掉变量X2可以使AIC达到最小，因此，R软件自动去掉变量，进行下一轮计算。

得到了最优回归方程。

再用summary（）;提取相关信息。

由显示结果看到：

回归系数检验的显著性水平有很大提高，但变量X3系数检验的显著性水平仍不理想。

下面用drop（）函数计算。

drop1（lm.step）

得到的结果如下：

Singletermdeletions

Model:

Y~X1+X3

DfSumofSqRSSAIC

5599.4109.3

X115169.510768.9119.1

X31833.26432.6109.8

计算结果表明：

如果去掉X3，则AIC值会从109.3增加到109.8，是增加最少的。

另外，无掉X3，残差的平方和上升833.2，也是最少的。

因此，从这两项标准来看，应该再去掉变量X3

lm.opt<-lm（Y~X1,data=cement）;

summary（lm.opt）

得到的结果如下：

Call:

lm（formula=Y~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

计算结果表明：

这个结果是满意的，因为所有的检验均是显著的，最后得到的“最优”的回归方程为：

4方差分析

解：

输入程序如下：

lamp<-data.frame（

X=c（1,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.5,4.9,6.0,6.84,8.23,10.33）,

A=factor（rep（1:

3,c（8,8,8）））

）

lamp.aov<-aov（X~A,data=lamp）

summary（lamp.aov）

lamp.lm<-lm（X~A,data=lamp）

anova（lamp.aov）

attach（lamp）

plot（X~A,xlab="FactorA",ylab="Life-SpanDataX",

main="Box-and-WhiskerPlot",

col=c（"yellow","orange","lightgreen","lightblue"））

detach（）

savePlot（"box_plot",type="eps"）

得到的结果如下：

DfSumSqMeanSqFvaluePr（>F）

A273.11836.5599.1040.001422**

Residuals2184.3294.016

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

>

结果说明：

从P值0.0014<0.05看出，拒绝H0,说明不同私聊的小鼠肝中的铁含量有显著差异。

从上面的箱线图可以看出，A、B、C互相有显著差异。

用

>attach（lamp）

>tapply（X,A,mean）

得到的均值分别为

123

1.606253.107505.82375

（3）输入程序：

>tapply（X,A,shapiro.test）

得到结果：

$`1`

Shapiro-Wilknormalitytest

data:

X[[1L]]

W=0.8742,p-value=0.1656

$`2`

Shapiro-Wilknormalitytest

data:

X[[2L]]

W=0.8893,p-value=0.2306

$`3`

Shapiro-Wilknormalitytest

data:

X[[3L]]

W=0.985,p-value=0.9833

结果说明：

数据是满足正态分布。

检验方差：

bartlett.test（X,A）

得到结果：

Bartletttestofhomogeneityofvariances

data:

XandA

Bartlett'sK-squared=10.5677,df=2,p-value=0.005073

结果说明：

P值>0.05,方差齐性不满足。

当数据只满足正态性，但不满足齐性要求，可用oneway.text（）作方差分析。

oneway.test（X~A,data=lamp）

得到的结论如下：

One-wayanalysisofmeans（notassumingequalvariances）

data:

XandA

F=10.3592,numdf=2.00,denomdf=10.51,p-value=0.003271

最终结论：

与原方法相同，满足正态性但不满足方差齐性。

4方差分析

解：

输入程序：

tree<-data.frame（

Y=c（4.6,4.3,6.1,6.5,6.8,6.4,6.3,6.7,3.4,3.8,

4.0,3.8,4.7,4.3,3.9,3.5,6.5,7.0）,

A=gl（3,6,18,labels=paste（"A",1:

3,sep=""））,

B=gl（3,2,18,labels=paste（"B",1:

3,sep=""））

）

tree.aov<-aov（Y~A+B+A:

B,data=tree）

summary（tree.aov）

tree.lm<-lm（Y~A+B+A:

B,data=tree）

anova（tree.lm）

attach（tree）

tapply（Y,A,mean）

tapply（Y,B,mean）

matrix（tapply（Y,A:

B,mean）,nr=3,nc=3,byrow=T,

dimnames=list（levels（A）,levels（B）））

interaction.plot（A,B,Y,l