北京工業(yè)大學(xué)數(shù)學(xué)建模選修第七次作業(yè)

1、數(shù)學(xué)建模第七次作業(yè)一、 回歸分析第一問：R程序> 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,10

2、1,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(y1+x)> summary(lm.sol, correlation = T, symbolic.cor = F)運(yùn)行結(jié)果：一元回歸模型：y=56.16+0.58x相關(guān)系數(shù)檢驗(yàn)結(jié)果為三顆星，p值均極小。第二問：R程序> new<-data.frame(x=40)> pp<-predict(lm.sol, new, interval="predictio

3、n")> pp fit lwr upr1 79.35816 62.86166 95.85466如果年齡為40歲，血壓預(yù)測(cè)值為79.36，血壓范圍為：（62.86,95.85）第三問：R程序pc<-predict(lm.sol, new, interval="confidence")pcmatplot(new$x, cbind(pp, pc,-1), type="l", xlab="X", ylab="Y", lty=c(1,5,5,2,2), col=c("blue", &

4、quot;red", "red", "brown", "brown"), lwd=2)points(x,y, cex=1.4, pch=21, col="red", bg="orange")legend(0.1, 63, c("Points", "Fitted", "Prediction", "Confidence"), pch=c(19, NA, NA, NA), lty=c(NA, 1,5,2), co

5、l=c("orange", "blue", "red", "brown")savePlot("predict", type="eps")二、 回歸分析和逐步回歸 第一問：將數(shù)據(jù)寫入.txt文件中，命名“data.txt”> health<-read.table("data.txt",header=T)> lm.sol<-lm(yx1+x2+x3+x4+x5+x6,data=health)> summary(lm.sol)運(yùn)行結(jié)

6、果：Call:lm(formula = y x1 + x2 + x3 + x4 + x5 + x6, data = health)Residuals: Min 1Q Median 3Q Max -5.3906 -0.9856 0.0749 1.0215 5.4075 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 104.86493 12.12725 8.647 7.75e-09 *x1 -0.24074 0.09459 -2.545 0.01778 * x2 -0.07456 0.05328 -1.399 0

7、.17448 x3 -2.62442 0.37249 -7.046 2.77e-07 *x4 -0.02532 0.06467 -0.391 0.69891 x5 -0.35992 0.11757 -3.061 0.00536 * x6 0.28766 0.13437 2.141 0.04266 * -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 2.267 on 24 degrees of freedomMultiple R-squared: 0.8552, Adjusted R-squared:

8、 0.819 F-statistic: 23.62 on 6 and 24 DF, p-value: 5.818e-09多元線性回歸方程：y=104.86-0.24x1-0.07x2-2.62x3-0.03x4-0.36x5-0.29x6第二問：由運(yùn)行結(jié)果知：方程F檢驗(yàn)P值小于0.05，否定原假設(shè)，即通過檢驗(yàn)。相關(guān)系數(shù)R2=0.8552,調(diào)整后R2=0.819，即因變量與自變量是相關(guān)的。系數(shù)中，常數(shù)項(xiàng)、x1、x3、x5、x6的P值小于0.05相關(guān)度較高，而x2,x4的相關(guān)度較低。第三問：R程序：> lm.step<-step(lm.sol,k=4)運(yùn)行結(jié)果：Start: AIC=7

9、0.8y x1 + x2 + x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC- x4 1 0.787 124.09 66.999- x2 1 10.062 133.37 69.233<none> 123.31 70.801- x6 1 23.546 146.85 72.219- x1 1 33.277 156.58 74.208- x5 1 48.153 171.46 77.021- x3 1 255.042 378.35 101.557Step: AIC=67y x1 + x2 + x3 + x5 + x6 Df Sum of Sq RSS AIC-

10、 x2 1 9.59 133.69 65.307<none> 124.09 66.999- x6 1 23.97 148.07 68.474- x1 1 32.69 156.78 70.247- x5 1 49.94 174.03 73.483- x3 1 327.46 451.56 103.040Step: AIC=65.31y x1 + x3 + x5 + x6 Df Sum of Sq RSS AIC<none> 133.69 65.307- x6 1 19.28 152.96 65.482- x1 1 28.09 161.77 67.218- x5 1 44.4

11、0 178.09 70.197- x3 1 358.81 492.50 101.730從程序的運(yùn)行結(jié)果可見，當(dāng)用全部變量作回歸方程時(shí)，AIC值為70.8，接下來(lái)的數(shù)據(jù)表顯示，如果去掉X2,X4,則相應(yīng)的ACI值為65.31，此時(shí)ACI值達(dá)到最小。對(duì)得到的回歸方程進(jìn)行顯著性分析,編寫R程序：> summary(lm.step)運(yùn)行結(jié)果：Call:lm(formula = y x1 + x3 + x5 + x6, data = health)Residuals: Min 1Q Median 3Q Max -4.9590 -1.2603 -0.1512 1.1796 4.8132 Coeffi

12、cients: Estimate Std. Error t value Pr(>|t|) (Intercept) 100.07910 11.57739 8.644 4.02e-09 *x1 -0.21266 0.09099 -2.337 0.02740 * x3 -2.76824 0.33138 -8.354 7.79e-09 *x5 -0.33957 0.11555 -2.939 0.00683 * x6 0.25535 0.13188 1.936 0.06378 . -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual sta

13、ndard error: 2.268 on 26 degrees of freedomMultiple R-squared: 0.843, Adjusted R-squared: 0.8188 F-statistic: 34.9 on 4 and 26 DF, p-value: 4.219e-10由運(yùn)行結(jié)果知：方程F檢驗(yàn)P值小于0.05，否定原假設(shè)，即通過檢驗(yàn)。相關(guān)系數(shù)R2=0.843,調(diào)整后R2=0.8188，即因變量與自變量是相關(guān)的。系數(shù)中，常數(shù)項(xiàng)、x1、x3、x5的P值小于0.05相關(guān)度較高，而x6稍微大于0.05，相關(guān)度大致認(rèn)為較高。軟件得到的回歸方程為：y=+100.07910-0.

14、21266 x1-2.76824x3-0.33957x5+0.25535x6三、 方差分析第一問：R程序：> coating<-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)+ )> coating.lm<-lm(xA,data=coating)> anova(coatin

15、g.lm)運(yùn)行結(jié)果：由運(yùn)行結(jié)果知，P值小于0.05，拒絕原假設(shè)，即不同涂料間存在顯著差異。第二問：R程序> attach(coating)> tapply(x,A,mean) 1 2 3 229.6 309.9 427.8由運(yùn)行結(jié)果知，第一組與二、三組有明顯差異，二三組間無(wú)明顯差異。四、 方差分析第一問R程序：> nd<-data.frame(+ y=c(14,10,11,11,13,9,10,12,9,7,10,8,7,11,6,10,5,11,13,14,12,13,14,10),+ B=gl(3,8,24,labels=paste("B",1:

16、3,sep=""),+ A=gl(3,2,24,labels=paste("A",1:3,sep="")+ )> nd.aov<-aov(yB+A+B:A,data=nd)> summary(nd.aov)結(jié)果：由運(yùn)行結(jié)果知，顯著性水平=0.05下，提高因素B對(duì)結(jié)果的影響是顯著的，提高因素A和AB交互作用對(duì)結(jié)果的影響不顯著。第二問：R程序> nd.lm<-lm(yB+A+B:A,data=nd> anova(nd.lm)Analysis of Variance TableResponse: y Df

17、 Sum Sq Mean Sq F value Pr(>F) B 2 44.333 22.1667 4.0549 0.0391 *A 2 8.867 4.4333 0.8110 0.4630 B:A 4 12.633 3.1583 0.5777 0.6832 Residuals 15 82.000 5.4667 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1> attach(nd)The following object is masked from coating: A> tapply(y,A,mean) A1 A2 A3 11.375 9.875 10.000 > tapply(y,