数学建模实验八.docx
- 文档编号:3687750
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:33
- 大小:1.10MB
数学建模实验八.docx
《数学建模实验八.docx》由会员分享,可在线阅读,更多相关《数学建模实验八.docx(33页珍藏版)》请在冰豆网上搜索。
数学建模实验八
数学建模实验八
(多元分析实验)
基本实验
1.回归分析
解答:
由于要分析血压与年龄两个变量之间的回归关系,可做一元回归,设年龄为自变量x,血压为因变量y.
(1)设方程为:
,(其中
和
是未知常数,为回归系数,
服从正态分布
)
利用R软件中的lm()可以非常方便地求出回归参数
和
,并作相应的检验,程序如下:
>rm(list=ls(all=TRUE))
>X=c(27,21,22,24,25,23,20,20,29,24,25,28,26,38,32,33,31,34,
++37,38,33,35,30,31,37,39,46,49,40,42,43,46,43,44,46,47,
++45,49,48,40,42,55,54,57,52,53,56,52,50,59,50,52,58,57)
>Y=c(73,66,63,75,71,70,65,70,79,72,68,67,79,91,76,69,66,73,
++78,87,76,79,73,80,68,75,89,101,70,72,80,83,75,71,80,96,
++92,80,70,90,85,76,71,99,86,79,92,85,71,90,91,100,80,109)
>lm.sol=lm(Y~1+X)
>summary(lm.sol)
Call:
lm(formula=y~1+x)
Residuals:
Min1QMedian3QMax
-16.4786-5.7877-0.07845.611719.7813
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)56.156933.9936714.061<2e-16***
x0.580030.096955.9832.05e-07***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
8.146on52degreesoffreedom
MultipleR-squared:
0.4077,AdjustedR-squared:
0.3963
F-statistic:
35.79on1and52DF,p-value:
2.05e-07
由结果可见,Estimate项中给出了回归方程的系数估计,即
,
,这两个参数的估计还是比较准确的,都有“***”,即显著性比较好。
该方程总体也通过了F统计数的检验,其p-value:
2.05e-07<<0.05,由此得到的回归方程为:
Y=56.15693+0.58003x
做残差分析,画出散点图:
>y.res=residuals(lm.sol)
>y.rst=rstandard(lm.sol)
>op=par(mfrow=c(2,2),mar=.1+c(4,4,1,1),oma=c(0,0,2,0))
>plot(fitted(lm.sol),y.res,
+cex=1.2,pch=21,col="red",bg="orange",
+xlab="FittedValue",ylab="Residuals")
>plot(fitted(lm.sol),y.rst,
+cex=1.2,pch=21,col="red",bg="orange",
+xlab="FittedValue",ylab="standardResiduals")
残差图:
从图中可以看出,回归值的大小与残差的波动大小有关系,属于异方差情况,即等方差性的假设有问题。
(2)正态性&方差齐次性检验
>shapiro.test(y.res)
Shapiro-Wilknormalitytest
data:
y.res
W=0.9901,p-value=0.9342
因此,残差满足正态性假设
由上面残差图的结果可知,不满足方差齐次性的要求。
(3)Box-Cox变换
选取lambda为-2.0,作box-cox变换,然后重新计算残差散点图,如下图中的左下图,从图中可以看出残差散点图有很大改善,大部分点均匀落在标准化残差的2以内,满足齐次性要求。
拟合后的曲线见右下图曲线。
总体程序如下:
rm(list=ls(all=TRUE))
X=c(27,21,22,24,25,23,20,20,29,24,25,28,26,38,32,33,31,34,37,38,33,35,30,31,37,39,46,49,40,42,43,46,43,44,46,47,45,49,48,40,42,55,54,57,52,53,56,52,50,59,50,52,58,57)
Y=c(73,66,63,75,71,70,65,70,79,72,68,67,79,91,76,69,66,73,78,87,76,79,73,80,68,75,89,101,70,72,80,83,75,71,80,96,92,80,70,90,85,76,71,99,86,79,92,85,71,90,91,100,80,109)
lm.sol=lm(Y~X);summary(lm.sol)
Y.rst=rstandard(lm.sol)
op=par(mfrow=c(2,2),mar=.1+c(4,4,1,1),oma=c(0,0,2,0))
plot(fitted(lm.sol),Y.rst,
cex=1.2,pch=21,col="red",bg="orange",
xlab="FittedValue",ylab="standardResiduals")
library(MASS)
####确定参数lambda
boxcox(lm.sol,lambda=seq(-10,2,by=0.01))
####Box-Cox变换后,作回归分析
lambda=-2;Ylam=(Y^lambda-1)/lambda
lm.lam=lm(Ylam~X);summary(lm.lam)
####变换后残差与预测散点图
plot(fitted(lm.lam),rstandard(lm.lam),
cex=1.2,pch=21,col="red",bg="orange",
xlab="FittedValue",ylab="standardResiduals")
beta0=lm.lam$coefficients[1]
beta1=lm.lam$coefficients[2]
curve((1+lambda*(beta0+beta1*x))^(1/lambda),
from=min(X),to=max(X),col="blue",lwd=2,
xlab="X",ylab="Y")
points(X,Y,pch=21,cex=1.2,col="red",bg="orange")
mtext("Box-CoxTransformations",outer=TRUE,cex=1.5)
par(op)
2.回归分析和逐步回归
解答:
(1)设方程为:
,(其中
是未知常数,为回归系数,
服从正态分布)
利用R软件中的lm()可以非常方便地求出回归参数
,并作相应的检验,程序如下:
>rm(list=ls(all=TRUE))
>X1=c(44,40,44,42,38,47,40,43,44,38,44,45,45,47,54,49,51,51,48,49,57,54,56,50,51,54,51,57,49,48,52)
>X2=c(89.47,75.07,85.84,68.15,814.27,77.45,75.98,81.19,81.42,81.87,73.03,87.66,66.45,79.15,83.12,81.42,69.63,77.91,91.63,73.37,73.37,79.38,76.32,70.87,67.25,91.63,73.71,59.08,76.32,61.24,82.78)
>X3=c(11.37,10.07,8.65,8.17,9.22,11.63,11.95,10.85,13.08,8.63,10.13,14.03,11.12,10.6,10.33,8.95,10.95,10,10.25,10.08,12.63,11.17,9.63,8.92,11.08,12.88,10.47,9.93,9.4,11.5,10.5)
>X4=c(62,62,45,40,55,58,70,64,63,48,45,56,51,47,50,44,57,48,48,67,58,62,48,48,48,44,59,49,56,52,53)
>X5=c(178,185,156,166,178,176,176,162,174,170,168,186,176,162,166,180,168,162,162,168,174,156,164,146,172,168,186,148,186,170,170)
>X6=c(182,185,168,172,180,176,180,170,176,186,168,192,176,164,170,185,172,168,164,168,176,165,166,155,172,172,188,155,188,176,172)
>Y=c(44.609,45.313,54.297,59.571,49.874,44.811,45.681,49.091,39.442,60.055,50.541,37.388,44.754,47.273,51.855,49.156,40.836,46.672,46.774,50.388,39.407,46.08,45.441,54.625,45.118,39.203,45.79,50.545,48.673,47.92,47.467)
>lm.sol=lm(Y~X1+X2+X3+X4+X5+X6)
>summary(lm.sol)
Call:
lm(formula=Y~X1+X2+X3+X4+X5+X6)
Residuals:
Min1QMedian3QMax
-5.3906-0.98560.07491.02155.4075
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)104.8649312.127258.6477.75e-09***
X1-0.240740.09459-2.5450.01778*
X2-0.074560.05328-1.3990.17448
X3-2.624420.37249-7.0462.77e-07***
X4-0.025320.06467-0.3910.69891
X5-0.359920.11757-3.0610.00536**
X60.287660.134372.1410.04266*
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
2.267on24degreesoffreedom
MultipleR-squared:
0.8552,AdjustedR-squared:
0.819
F-statistic:
23.62on6and24DF,p-value:
5.818e-09
(2)由运行结果可得:
回归方程为:
Y=104.86493-0.24074X1-0.07456X2-2.62442X3-0.02532X4-0.35992X5+0.28766X6
p-value值为5.818e-09<0.05方程本身是通过检测的
常数项显著;X1项系数显著;X2项系数不显著;X3项系数极为显著;X4项系数不显著;X5项系数高度显著;X6项系数显著。
有两项系数没有通过检验,总体来说拟合并不理想。
(3)做逐步回归分析:
>lm.step=step(lm.sol)
Start:
AIC=56.8
Y~X1+X2+X3+X4+X5+X6
DfSumofSqRSSAIC
-X410.787124.0954.999
-X2110.062133.3757.233
-X6123.546146.8560.219
-X1133.277156.5862.208
-X5148.153171.4665.021
-X31255.042378.3589.557
Step:
AIC=55
Y~X1+X2+X3+X5+X6
DfSumofSqRSSAIC
-X219.59133.6955.307
-X6123.97148.0758.474
-X1132.69156.7860.247
-X5149.94174.0363.483
-X31327.46451.5693.040
由运行结果可得:
程序去掉X4项,AIC值降为55
>summary(lm.step)
Call:
lm(formula=Y~X1+X2+X3+X5+X6)
Residuals:
Min1QMedian3QMax
-5.4716-0.9252-0.01280.95895.4057
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)103.9969611.719198.8743.37e-09***
X1-0.232250.09050-2.5660.01666*
X2-0.072410.05209-1.3900.17676
X3-2.686910.33081-8.1221.78e-08***
X5-0.364610.11495-3.1720.00398**
X60.289970.131952.1980.03746*
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
2.228on25degreesoffreedom
MultipleR-squared:
0.8542,AdjustedR-squared:
0.8251
F-statistic:
29.3on5and25DF,p-value:
1.083e-09
由运行结果可得:
回归系数检验的水平已有显著提升,但X2项系数仍然不显著。
>drop1(lm.step)
Singletermdeletions
Model:
Y~X1+X2+X3+X5+X6
DfSumofSqRSSAIC
X1132.69156.7860.247
X219.59133.6955.307
X31327.46451.5693.040
X5149.94174.0363.483
X6123.97148.0758.474
由结果可得:
去掉X2项的时候,AIC值和残差平方值上升最小,因此在做回归分析时,要去掉X2,编写如下程序:
>lm.opt=lm(Y~X1+X3+X5+X6);
>summary(lm.opt)
Call:
lm(formula=Y~X1+X3+X5+X6)
Residuals:
Min1QMedian3QMax
-4.9590-1.2603-0.15121.17964.8132
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)100.0791011.577398.6444.02e-09***
X1-0.212660.09099-2.3370.02740*
X3-2.768240.33138-8.3547.79e-09***
X5-0.339570.11555-2.9390.00683**
X60.255350.131881.9360.06378.
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
2.268on26degreesoffreedom
MultipleR-squared:
0.843,AdjustedR-squared:
0.8188
F-statistic:
34.9on4and26DF,p-value:
4.219e-10
由结果可得:
除X6不太显著外,其余常数项的检测结果均显著。
p-value=4.219e-10<0.05,比之前的值都低。
由此得到了最优回归方程:
Y=11.57739-0.21266X1-2.76824X3-0.33957X5+0.25535X6
3.
方差分析I
解答:
(1)我们分别使用R软件中的aov()函数以及lm()函数进行方差分析:
>rm(list=ls(all=TRUE))
>time=data.frame(
X=c(148,76,393,520,236,134,55,166,415,153,513,264,433,94,535,327,214,135,280,304,335,643,216,536,128,723,258,380,594,465),
A=factor(rep(1:
3,c(10,10,10))))
>attach(time)
>plot(X~A,xlab="FactorA",ylab="Time-SpanDataX",
main="Box-and-WhiskerPlot",
col=c("yellow","orange","lightgreen","lightblue"))
>detach()
>time.aov=aov(X~A,data=time)
>summary(time.aov)
DfSumSqMeanSqFvaluePr(>F)
A2198772993863.4820.0452*
Residuals2777067128543
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
>time.lm=lm(X~A,data=time)
>anova(time.aov)
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr(>F)
A2198772993863.48190.04515*
Residuals2777067128543
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
由结果可知,两种函数计算出的p值为0.0452、0.04515,均小于0.05,拒绝原假设,即三种涂料的磨损平均时间存在明显差异。
磨损时间试验的箱线图:
(2)使用多重t检验分析差异,下面给出了调整p值前后的结果:
>####求数据在各水平下的均值
>attach(time)
Thefollowingobjectsaremaskedfromtime(pos=3):
A,X
>tapply(X,A,mean)
123
229.6309.9427.8
>####作多重t检验
>pairwise.t.test(X,A)
PairwisecomparisonsusingttestswithpooledSD
data:
XandA
12
20.297-
30.0420.261
Pvalueadjustmentmethod:
holm
>
>pairwise.t.test(X,A,p.adjust.method="none")
PairwisecomparisonsusingttestswithpooledSD
data:
XandA
12
20.297-
30.0140.130
Pvalueadjustmentmethod:
none
由结果可知,无论是否使用holm方法调整p值,结果都为:
1、2;2、3差异不显著
1、3差异显著
即三种涂料的磨损平均时间,A、B;B、C的差异不显著,A、C的差异较显著。
其实从箱线图也能观察到该结论。
4.
方差分析II
解答:
使用oneway.test()做方差分析:
>rm(list=ls(all=TRUE))
>time=data.frame(
X=c(148,76,393,520,236,134,55,166,415,153,513,264,433,94,535,327,214,135,280,304,335,643,216,536,128,723,258,380,594,465),
A=factor(rep(1:
3,c(10,10,10))))
oneway.test(X~A,data=time,var.equal=T)
One-wayanalysisofmeans
data:
XandA
F=3.4819,numdf=2,denomdf=27,p-value=0.04515
由结果可知p-value=0.04515,与问题3得到的结果相同,拒绝原假设,即三种涂料的磨损平均时间存在明显差异。
使用Kruskal-Wallis秩和检验:
>time=data.frame(
X=c(148,76,393,520,236,134,55,166,415,1
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 建模 实验
![提示](https://static.bdocx.com/images/bang_tan.gif)