数学建模课后作业第七章.docx
《数学建模课后作业第七章.docx》由会员分享，可在线阅读，更多相关《数学建模课后作业第七章.docx（41页珍藏版）》请在冰点文库上搜索。
数学建模课后作业第七章
第七章.多元分析实验
基本实验
1.线性回归；
解：
由题可以得出如下的R程序：
>X1<-c,,,,,,,,,,239）
>X2<-c,,,,,,,,,,
>X3<-c,,,,,,,,,,
>Y<-c,,19,,,,,,,,
>
><-lm（Y~X1+X2+X3）
>summary
运行后可以得知；
Call:
lm（formula=Y~X1+X2+X3）
Residuals:
Min1QMedian3QMax
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）***
X1
X2***
X3*
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
Residualstandarderror:
on7degreesoffreedom
MultipleR-squared:
AdjustedR-squared:
F-statistic:
on3and7DF,p-value:
则可以得出Y关于X1、X2、X3的线性回归方程；
Y=X2+
由上述的结果可以得知方程的常量与X2显著性为***表示十分的显著，X3显著性为*表示显著，而X2为不显著。
（2）由
（1）中的数据可以得知新的分析函数anovaR程序如下：
X1<-c,,,,,,,,,,239）
X2<-c,,,,,,,,,,
X3<-c,,,,,,,,,,
Y<-c,,19,,,,,,,,
<-lm（Y~X1+X2+X3,data=blood）
summary
anova
运行后可以得出：
Min1QMedian3QMax
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）***
X1
X2***
X3*
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
Residualstandarderror:
on7degreesoffreedom
MultipleR-squared:
AdjustedR-squared:
F-statistic:
on3and7DF,p-value:
>
>anova
AnalysisofVarianceTable
Response:
Y
DfSumSqMeanSqFvaluePr（>F）
X11***
X21***
X31*
Residuals7
---
Signif.codes:
0‘***’‘**’‘*’‘.’
由此结果可以看出X1、X2、X3均能通过显著性检验，所以选择全部变量作回归方程是十分合理的。
然后可以得出对变量作逐步回归的R程序如下：
<-step
运行后：
Start:
AIC=
Y~X1+X2+X3
DfSumofSqRSSAIC
-X11
-X31
-X21
Step:
AIC=
Y~X2+X3
DfSumofSqRSSAIC
-X21
-X31
接着继续显著性R程序分析：
summary
运行可得：
Call:
lm（formula=Y~X2+X3）
Residuals:
Min1QMedian3QMax
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）***
X2***
X3***
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
Residualstandarderror:
on8degreesoffreedom
MultipleR-squared:
AdjustedR-squared:
F-statistic:
on2and8DF,p-value:
可以得出更加恰当与精确Y与X1、X2、X3的线性关系如下如下：
Y=++
（3）由
（1）与
（2）中的线性方程可分别得出对应的Y值为
和
由题目已知条件可以得出如下R程序
new<
<-predict,new,interval=”prediction”,level=
运行后可以得出：
fitlwrupr
[1,]
可以得出预测区间为[,],预测值大致为
置信区间则为[,]
2．方差分析I（单因素方差分析）；
解:
作出如下假设命令；
H0：
三个厂产品的零件强度无差异,即二者方差相同；
H1：
三个厂产品的零件强度无有异,即二者方差不相同；
由题可以得出关于三个工厂产品检测数据差异的方差分析R程序如下：
products<-c（115,116,98,83,103,107,118,116,73,89,85,97）
A<-factor（rep（1:
3,c（4,4,4）））
<-aov（products~A）
summary
运行程序可以得出：
>products<-c（115,116,98,83,103,107,118,116,73,89,85,97）
>A<-factor（rep（1:
3,c（4,4,4）））
><-aov（products~A）
>summary
DfSumSqMeanSqFvaluePr（>F）
A21304*
Residuals91192
---
Signif.codes:
0‘***’‘**’‘*’‘.’
由于Pr=<,所以假设H1成立
（2）由题目条件可以得出如下关于强度均值的R程序分析：
products<-c（115,116,98,83,103,107,118,116,73,89,85,97）
A<-factor（rep（1:
3,c（4,4,4）））
运行后可以得出：
PairwisecomparisonsusingttestswithpooledSD
data:
productsandA
12
2-
3
Pvalueadjustmentmethod:
holm
由此可以得出三个工厂产品的产品强度均不相同，且存在一定的差异；
（3）由题目要求可以得出如下多重分析的R软件：
1.
mouse<（
X=c（115,116,98,83,103,107,118,116,73,89,85,97）,
A=factor（rep（1:
3,c（4,4,4）））
）
<-lm（X~A,data=mouse）
anova
attach（mouse）
tapply（X,A,mean）
A）
A,="none"）
plot（X~A,col=5:
7,
main="Box-and-WhiskerPlotofMouseData"）
detach（mouse）
savePlot（"box_plot2",type="eps"）
products<-c（115,116,98,83,103,107,118,116,73,89,85,97）
运行后可得出：
>mouse<（
+X=c（115,116,98,83,103,107,118,116,73,89,85,97）,
+A=factor（rep（1:
3,c（4,4,4）））
+）
><-lm（X~A,data=mouse）
>anova
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr（>F）
A21304*
Residuals91192
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
>
>
>attach（mouse）
Thefollowingobject（s）aremasked_by_'.GlobalEnv':
A
>tapply（X,A,mean）
123
10311186
>
>
>A）
PairwisecomparisonsusingttestswithpooledSD
data:
XandA
12
2-
3
Pvalueadjustmentmethod:
holm
>
>A,="none"）
PairwisecomparisonsusingttestswithpooledSD
data:
XandA
12
2-
3
Pvalueadjustmentmethod:
none
>
>plot（X~A,col=5:
7,
+main="Box-and-WhiskerPlotofMouseData"）
>detach（mouse）
>
>savePlot（"box_plot2",type="eps"）
由以上结果可以得知三者之间存在显著地差异性；
2.
mouse<（
X=c（115,116,98,83,103,107,118,116,73,89,85,97）,
A=factor（rep（1:
3,c（4,4,4）））
）
<-lm（X~A,data=mouse）
anova
<-aov（X~A,data=mouse）
anova
（X~A,data=mouse）
（X~A,data=mouse,=T）
运行程序可以得出：
>mouse<（
+X=c（115,116,98,83,103,107,118,116,73,89,85,97）,
+A=factor（rep（1:
3,c（4,4,4）））
+）
><-lm（X~A,data=mouse）
>anova
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr（>F）
A21304*
Residuals91192
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
>
><-aov（X~A,data=mouse）
>anova
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr（>F）
A21304*
Residuals91192
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
>
>（X~A,data=mouse）
One-wayanalysisofmeans（notassumingequalvariances）
data:
XandA
F=,numdf=,denomdf=,p-value=
>（X~A,data=mouse,=T）
One-wayanalysisofmeans
data:
XandA
F=,numdf=2,denomdf=9,p-value=
从上述结果可以得出p-value=，p-value=二者均小于说明三个工厂的产品之间有显著性差异；
3．方差分析II（双因素方差分析）；
解：
（1）由实验数据可以得出如下的方差分析R程序：
tree<（
Y=c,,,,,,,,,,,,,,,,,,
A=gl（3,6,18,labels=paste（"A",1:
3,sep=""））,
B=gl（3,2,18,labels=paste（"B",1:
3,sep=""））
）
<-aov（Y~A+B+A:
B,data=tree）
summary
运行后可以得出：
>tree<（
+Y=c,,,,,,,,,,,,,,,,,,
+A=gl（3,6,18,labels=paste（"A",1:
3,sep=""））,
+B=gl（3,2,18,labels=paste（"B",1:
3,sep=""））
+）
><-aov（Y~A+B+A:
B,data=tree）
>summary
DfSumSqMeanSqFvaluePr（>F）
A2***
B2***
A:
B4***
Residuals9
---
Signif.codes:
0‘***’‘**’‘*’‘.’
由程序运行结果可以得知反应温度与反应压力的影响均是高度显著地的，而且二者之间也有着高度的交互作用显著性。
（2）由题中已知条件可以得到关于反应温度与压力均值表格如下
A1（60º）
A2（70º）
A3（80º）
均值
B1（2公斤）
B2（公斤）
4
B3（3公斤）
7
均值
进一步可以得到反应温度与压力综合的表格
产量均值
A1（60º）
A2（70º）
A3（80º）
B1（2公斤）
B2（公斤）
B3（3公斤）
从表中数据可以看出当反应压力B3（3公斤），反应温度为A3（80º）为最优的生产方案。
（3）从
（2）中表格的数据分析可以看出产量随着压力的增加先降低后增加，随温度的增长也是先增后降，但是其规律存在一定的差异，从表中数据可以看出压力与温度都很高时产量较大，但是压力与温度过高就会对设备、人员、技术提出很高的要求，使生产成本过高，而适当的压力与温度条件下也可以得到较高的产量，所以综合各方面的因素而言适当的反应温度与反应压力获得相对较高的产量降低生产的各方面要求以获取较大的生产利益为最优策略。
4.判别分析
解
（1）对Af蠓虫与Apf蠓虫做判别分析，
一．距离判别
由题目要求可以得出如下的R程序判别分析
<-function
（TrnX1,TrnX2,TstX=NULL,=FALSE）{
if（TstX）==TRUE）TstX<-rbind（TrnX1,TrnX2）
if（TstX）==TRUE）TstX<-t（TstX））
elseif（TstX）!
=TRUE）
TstX<（TstX）
if（TrnX1）!
=TRUE）TrnX1<（TrnX1）
if（TrnX2）!
=TRUE）TrnX2<（TrnX2）
nx<-nrow（TstX）
blong<-matrix（rep（0,nx）,nrow=1,byrow=TRUE,
dimnames=list（"blong",1:
nx））
mu1<-colMeans（TrnX1）;mu2<-colMeans（TrnX2）
if==TRUE||==T）{
S<-var（rbind（TrnX1,TrnX2））
w<-mahalanobis（TstX,mu2,S）-mahalanobis（TstX,mu1,S）
}
else{
S1<-var（TrnX1）;S2<-var（TrnX2）
w<-mahalanobis（TstX,mu2,S2）-mahalanobis（TstX,mu1,S1）
}
for（iin1:
nx）{
if（w[i]>0）
blong[i]<-1
else
blong[i]<-2
}
blong
}
1.对于输入标本,
classX1<-matrix（
c,,,,,,,,,,,,,,,,,,
nrow=9,byrow=F）
classX2<-matrix（
c,,,,,,,,,,,,
nrow=6,byrow=F）
test<-matrix（c,,nrow=1,byrow=T）
（classX1,classX2,test,=TRUE）
（classX1,classX2,test）
运行后可得：
>（classX1,classX2,test,=TRUE）
1
blong1
>（classX1,classX2,test）
1
blong2
可以得知当方差相同时标本,属于Af蠓虫；
当方差不相同时标本,属于Apf蠓虫；
2.对于标本（，）
classX1<-matrix（
c,,,,,,,,,,,,,,,,,,
nrow=9,byrow=F）
classX2<-matrix（
c,,,,,,,,,,,,
nrow=6,byrow=F）
test<-matrix（c,,nrow=1,byrow=T）
（classX1,classX2,test,=TRUE）
（classX1,classX2,test）
运行后可以得出；
>（classX1,classX2,test,=TRUE）
1
blong2
>（classX1,classX2,test）
1
blong2
可以得知当方差相同时标本属于Apf蠓虫；
当方差不相同时标本属于Apf蠓虫；
3．对于标本（，）；
classX1<-matrix（
c,,,,,,,,,,,,,,,,,,
nrow=9,byrow=F）
classX2<-matrix（
c,,,,,,,,,,,,
nrow=6,byrow=F）
test<-matrix（c,,nrow=1,byrow=T）
运行后可以得出：
>（classX1,classX2,test,=TRUE）
1
blong2
>（classX1,classX2,test）
1
blong1
可以得知当方差相同时标本属于Apf蠓虫；
当方差不相同时标本属于Af蠓虫；
二．Fisher判别
可以得出如下的判别程序：
<-function（TrnX1,TrnX2,TstX=NULL）{
if（TstX）==TRUE）TstX<-rbind（TrnX1,TrnX2）
if（TstX）==TRUE）TstX<-t（TstX））
elseif（TstX）!
=TRUE）
TstX<（TstX）
if（TrnX1）!
=TRUE）TrnX1<（TrnX1）
if（TrnX2）!
=TRUE）TrnX2<（TrnX2）
nx<-nrow（TstX）
blong<-matrix（rep（0,nx）,nrow=1,byrow=TRUE,
dimnames=list（"blong",1:
nx））
n1<-nrow（TrnX1）;n2<-nrow（TrnX2）
mu1<-colMeans（TrnX1）;mu2<-colMeans（TrnX2）
S<-（n1-1）*var（TrnX1）+（n2-1）*var（TrnX2）
mu<-n1/（n1+n2）*mu1+n2/（n1+n2）*mu2
w<-（TstX-rep（1,nx）%o%mu）%*%solve（S,mu2-mu1）;
for（iin1:
nx）{
if（w[i]<=0）
blong[i]<-1
else
blong[i]<-2
}
blong
}
对于标本1（,）有如下的判别程序：
classX1<（
x1=c,,,,,,,,,
x2=c,,,,,,,,
）
classX2<（
x1=c,,,,,,
x2=c,,,,,
）
test<-matrix（c,,nrow=1,byrow=T）
（classX1,classX2,test）
运行后可以得出：
>test<-matrix（c,,nrow=1,byrow=T）
>（classX1,classX2,test）
1
blong2
可以得知无论方差是否相同，标本都属于Apf蠓虫；
对于标本2（,）有如下的判别程序：
classX1<（
x1=c,,,,,,,,,
x2=c,,,,,,,,
）
classX2<（
x1=c,,,,,,
x2=c,,,,,
）
test<-matrix（c,,nrow=1,byrow=T）
（classX1,classX2,test）
运行后可以得出：
>test<-matrix（c,,nrow=1,byrow=T）
>（classX1,classX2,test）
1
blong2
可以得知无论方差是否相同，标本都属于Apf蠓虫
对于标本3（,）有如下的判别程序：
classX1<（
x1=c,,,,,,,,,
x2=c,,,,,,,,
）
classX2<（
x1=c,,,,,,
x2=c,,,,,
）
test<-matrix（c,,nrow=1,byrow=T）
（classX1,classX2,test）
运行程序后可以得出：
>test<-matrix（c,,nrow=1,byrow=T）
>（classX1,classX2,test）
1
blong2
可以得知无论方差是否相同，标本都属于Apf蠓虫；
（2）由题中已知条件可以得出如下的lda（）和qda（）函数判别程序：
对于标本1（,）
classX1<-matrix（
c,,,,,,,,,,,,,,,,,,
ncol=2）
classX2<-matrix（
c,,,,,,,,,,,,
ncol=2）
test<-c,
train<-rbind（classX1,classX2）
test<-rbind（test）
f<-factor（c（rep（“1”,9）,rep（”2”,6）））
m<-lda（train,f）
predict（m,test）$class
n<-qda（train,f）
predict（n,test）$class
运行程序后可以得出：
>predict（m,test）$class
[1]
升级会员 