《统计建模与R软件》书本课后习题答案
统计建模与R软件习题仅供参考
原假设:油漆工人的血小板计数与正常成年男子无差异。
备择假设:油漆工人的血小板计数与正常成年男子有差异。
因为p值小于,拒绝原假设,则认为油漆工人的血小板计数与正常成年男子有差异。
由程序结果可知,x<=1000的概率为,所以x大于1000的概率为.P值为,大于,接受原假设,则两种方法治疗后患者血红蛋白无差异。
(1)以上检验中p值均大于,接受原假设,数据来自正态分布。
(2)三种检验的结果都显示两组数据均值无差异。
(3)因为p值为,大于,所以接受原假设,两组数据方差相同。
(1)因为两组数据p值都大于,所以均接受原假设,两组数据都服从正态分布。
(2)P=>,接受原假设,可认为两组样本方差相同。
(3)P值小于,拒绝原假设,两组有差别。
P 值>,故接受原假设,表示调查结果支持该市老年人口的看法。
不能认为这种处理能增加母鸡的比例。
接受原假设,符合自由组合定律。
因为p值大于,接受原假设,可以认为X服从poission分布。
因为p大于,接受原假设,所以两样本来自同一总体。
因为p值小于,拒绝原假设,有影响。
因为P 值< ,所以拒绝原假设,B与C不独立。
由于p值大于,故两变量独立,两种工艺对产品的质量没有影响。
因为p值大于,接受原假设,可以认定两种方法测定结果相同。
(1)p<,拒绝原假设,中位数小于(2)p<,所以拒绝原假设,中位数小于(1)p值大于,接受原假设,无差别。
(2)p=<,拒绝原假设,有差别。
(3)p=<,拒绝原假设,有差别。
(4)可认为两组数据方差相同。
综上,该数据可做t检验,p值小于,拒绝原假设,有差别。
(5)综上所述,Wilcoxon符号秩检验的差异检出能力最强,符号检验的差异检出最弱。
(6)P值均小于,接受原假设,二者有关系,呈正相关。
因为p值大于,接受原假设,尚不能认为新方法的疗效显着优于原疗法。
(完整版)统计建模与R软件课后答案
第二章2.1> x<-c(1,2,3);y<-c(4,5,6)> e<-c(1,1,1)> z<-2*x+y+e;z[1] 7 10 13> z1<-crossprod(x,y);z1[,1][1,] 32> z2<-outer(x,y);z2[,1] [,2] [,3][1,] 4 5 6[2,] 8 10 12[3,] 12 15 182.2(1)> A<-matrix(1:20,nrow=4);B<-matrix(1:20,nrow=4,byrow=T) > C<-A+B;C(2)> D<-A%*%B;D(3)> E<-A*B;E(4)> F<-A[1:3,1:3](5)> G<-B[,-3]> x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x2.4> H<-matrix(nrow=5,ncol=5)> for (i in 1:5)+ for(j in 1:5)+ H[i,j]<-1/(i+j-1)(1)> det(H)(2)> solve(H)(3)> eigen(H)2.5> studentdata<-data.frame(姓名=c('张三','李四','王五','赵六','丁一')+ ,性别=c('女','男','女','男','女'),年龄=c('14','15','16','14','15'),+ 身高=c('156','165','157','162','159'),体重=c('42','49','41.5','52','45.5')) 2.6> write.table(studentdata,file='student.txt')> write.csv(studentdata,file='student.csv')2.7count<-function(n){if (n<=0)print('要求输入一个正整数')repeat{if (n%%2==0)n<-n/2elsen<-(3*n+1)if(n==1)break}print('运算成功')}}第三章3.1首先将数据录入为x。
统计建模与R软件课后参考答案(可编辑修改word版)
第二章2.1> x<-c(1,2,3);y<-c(4,5,6)> e<-c(1,1,1)> z<-2*x+y+e;z[1] 7 10 13>z1<-crossprod(x,y);z1[,1][1,] 32>z2<-outer(x,y);z2[,1] [,2] [,3][1,] 4 5 6[2,] 8 10 12[3,] 12 15 182.2(1) > A<-matrix(1:20,nrow=4);B<-matrix(1:20,nrow=4,byrow=T) >C<-A+B;C(2) > D<-A%*%B;D(3) > E<-A*B;E(4) > F<-A[1:3,1:3](5) > G<-B[,-3]2.3>x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x2.4>H<-matrix(nrow=5,ncol=5)>for (i in 1:5)+ for(j in 1:5)+ H[i,j]<-1/(i+j-1)(1)> det(H)(2)> solve(H)(3)> eigen(H)2.5>studentdata<-data.frame(姓名=c('张三','李四','王五','赵六','丁一') + ,性别=c('女','男','女','男','女'),年龄=c('14','15','16','14','15'),+ 身高=c('156','165','157','162','159'),体重=c('42','49','41.5','52','45.5')) 2.6>write.table(studentdata,file='student.txt')>write.csv(studentdata,file='student.csv')2.7count<-function(n){if (n<=0)print('要求输入一个正整数')else{ repeat{if (n%%2==0)n<-n/2elsen<-(3*n+1)if(n==1)break}print('运算成功')}}第三章3.1首先将数据录入为x。
R软件课后习题解第四章-
习题4.1计算代码如下(极大似然估计):x<-c(0.1,0.2,0.9,0.8,0.7,0.7)n<- length(x)f<-function(a){n/(a+1)+sum(log(x))} uniroot(f,c(0,1))$root[1] 0.211182$f.root[1] -3.844668e-05$iter[1] 5$estim.prec[1] 6.103516e-05习题4.2(指数分布)x<-c(5,15,25,35,45,55,65)v<-c(365,245,150,100,70,45,25)y<-x*vf<-function(k){1000/k-sum(y)} uniroot(f,c(0,100))$root[1] 0.05002618$f.root[1] -10.46586$iter[1] 14$estim.prec[1] 6.103516e-05估计值为0.05002618习题4.3(极大似然估计)(作业)解:设n x x x ⋯,,21为字样n ξξξ,...,21的一组观测值,所以似然函数为λλλλλλn n x n i i x n e x x x e x x x x L L n i i i -=-∑==⋯==∏!!...!!),,;()(211211取对数,所以 ∑∑+=-⋅+-=n i in i i x x n L 11)!l n (ln )(ln λλλ 求偏导,并令其等于0,01ln 1=+-=∂∂∑=ni i x n L λλ 解得,-=x ^λ,所以λ的极大似然估计量为-ξ.x<-rep(0:6,c(17,20,10,2,1,0,0))mean(x)[1] 1习题4.4y<-function(x){f<-c(-13+x[1]+((5-x[2])*x[2]-2)*x[2],-29+x[1]+((x[2]+1)*x[2]-14)*x[2]);sum(f^2)}x0<-c(0.5,-2)nlm(y,x0)$minimum[1] 48.98425$estimate[1] 11.4127791 -0.8968052$gradient[1] 1.415447e-08 -1.435296e-07$code[1] 1$iterations[1] 16习题4.5X<- c(54,67,68,78,70,66,67,70,65,69)t.test(X)One Sample t-testdata: Xt = 35.947, df = 9, p-value = 4.938e-11alternative hypothesis: true mean is not equal to 095 percent confidence interval:63.1585 71.6415sample estimates:mean of x67.4因此,10名患者平均脉搏在95%的置信区间为[63.16,71.64] 10个人的平均脉搏为67.4,所以这10名患者的平均脉搏属不低于正常人的平均脉搏习题4.6(作业)x<-c(140,137,136,140,145,148,140,135,144,141)y<-c(135,118,115,140,128,131,130,115,131,125)t.test(x,y)Welch Two Sample t-testdata: x and yt = 4.6287, df = 13.014, p-value = 0.0004712alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:7.359713 20.240287sample estimates:mean of x mean of y140.6 126.8所以u1-u2的置信区间为[7.53626 ,20.06374]习题4.7x<-c(0.143,0.142,0.143,0.137)y<-c(0.140,0.142,0.136,0.138,0.140)t.test(x,y,var.equal=TRUE) ## 注意:如果方差相同,需要声明var.equal=TRUE Two Sample t-testdata: x and yt = 1.198, df = 7, p-value = 0.2699alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:-0.001996351 0.006096351sample estimates:mean of x mean of y0.14125 0.13920习题4.8(作业)x<-c(140,137,136,140,145,148,140,135,144,141)y<-c(135,118,115,140,128,131,130,115,131,125)var.test(x,y)F test to compare two variancesdata: x and yF = 0.2353, num df = 9, denom df = 9, p-value = 0.04229alternative hypothesis: true ratio of variances is not equal to 195 percent confidence interval:0.05845276 0.94743902sample estimates:ratio of variances0.2353305t.test(x,y)Welch Two Sample t-testdata: x and yt = 4.6287, df = 13.014, p-value = 0.0004712alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:7.359713 20.240287sample estimates:mean of x mean of y140.6 126.8所以所求置信区间为[7.359713,20.240287]习题4.9X<-rep(0:6,c(7,10,12,8,3,2,0))t.test(X)One Sample t-testdata: Xt = 9.0895, df = 41, p-value = 2.238e-11alternative hypothesis: true mean is not equal to 095 percent confidence interval:1.4815562.327968sample estimates:mean of x1.904762所以估计值为1.904762,置信区间为[1.481556,2.327968]习题4.10(作业)X<-c(1067,919,1196,785,1126,936,918,1156,920,948)t.test(X,alternative="greater")One Sample t-testdata: Xt = 23.9693, df = 9, p-value = 9.148e-10alternative hypothesis: true mean is greater than 095 percent confidence interval:920.8443 Infsample estimates:mean of x997.1所以有95%的灯泡寿命在920.8443h以上.习题5.1x<-c(220, 188, 162, 230, 145, 160, 238, 188, 247, 113, 126, 245, 164, 231, 256, 183, 190, 158, 224, 175)t.test(x,mu=225)One Sample t-testdata: xt = -3.4783, df = 19, p-value = 0.002516alternative hypothesis: true mean is not equal to 22595 percent confidence interval:172.3827 211.9173sample estimates:mean of x192.15原假设:油漆工人的血小板计数与正常成年男子无差异。
统计建模与R语言习题6.6答案
统计建模与R语言习题6.6答案Homework6.6 P333解:设抗生素为x1(有用抗生素为1,没有为0)危险因子为x2(有危险因子为1,没有为0)计划为x3(事先有计划为1,临时决定为0)感染为Y整个R语言计算过程:doc<-data.frame(x1<-c(1,1,1,1,0,0,0,0),x2<-c(1,1,0,0,1,1,0,0),x3<-c(1,0,1,0,1,0,1,0) ,success<-c(1,11,0,0,28,23,8,0),fail<-c(17,87,2,0,30,3,32,9))doc$Ymat<-cbind(doc$success,doc$fail)glm.sol<-glm(Ymat~x1+x2+x3,family=binomial,data=doc) summary(glm.sol)pre<-predict(glm.sol,data.frame(x1=1,x2=1,x3=0))p<-exp(pre)/(1+exp(pre));ppre<-predict(glm.sol,data.frame(x1=1,x2=1,x3=1))p<-exp(pre)/(1+exp(pre));ppre<-predict(glm.sol,data.frame(x1=0,x2=1,x3=1))p<-exp(pre)/(1+exp(pre));ppre<-predict(glm.sol,data.frame(x1=1,x2=0,x3=1))p<-exp(pre)/(1+exp(pre));p回归模型输出结果如下:Call:glm(formula = Ymat ~ x1 + x2 + x3, family = binomial, data = doc)Deviance Residuals:1 2 3 4 5 6 7 8 0.26470 -0.07162 -0.15231 0.00000 -0.785201.49623 1.21563 -2.56229 Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -0.8207 0.4947 -1.659 0.0971 .x1 -3.2544 0.4813 -6.761 1.37e-11 ***x2 2.0299 0.4553 4.459 8.25e-06 ***x3 -1.0720 0.4254 -2.520 0.0117 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 83.491 on 6 degrees of freedomResidual deviance: 10.997 on 3 degrees of freedomAIC: 36.178Number of Fisher Scoring iterations: 5从上述结果即有常数项以及x项的系数β分别为-0.8207,-3.2544, 2.0299,-1.0720,并且回归方程通过检验得到回归模型:P=exp(-0.8207-3.2544x1+2.0299x2-1.0720x3)/(1+exp(-0.8207-3.2544x1+2.0299x2-1.0720x3))假设医生需要对一产妇进行剖腹手术,令x1=1,x2=1,x3=0即在手术中是用抗生素且有危险因子和临时决定进行手术,那么pre<-predict(glm.sol,data.frame(x1=1,x2=1,x3=0))p<-exp(pre)/(1+exp(pre));p10.1145421即11.4521%,也就是说在临时决定进行手术的感染率为11.4521% 那么当x3=1时pre<-predict(glm.sol,data.frame(x1=1,x2=1,x3=1))> p<-exp(pre)/(1+exp(pre));p10.04240619即4.240619%,也就是说在事先有计划进行手术的感染率为4.240619%。
统计建模与R软件(薛毅)第九章答案
第九章9.1(1)利用主成分确定了8个指标的主成分,有4个,如图(21)(2)用order()分别对4个主成分的预测值进行排序,结果是如下表(26),而利用kmeans()进行动态排序得到如下分类:第1类:建材(6),森工(7),食品(8),纺织(9),皮革(11);第2类:机械(5);第3类:电力(2),煤炭(3),缝纫(10)造纸(12);第4类:冶金(1)化学(4),文教艺术用品(13)。
成分13个行业排序结果第一主成分: 5 1 3 2 4 6 13 11 9 7 12 10 8 第二主成分: 5 8 4 9 10 1 13 12 7 11 6 2 3 第三主成分:8 1 5 3 9 12 7 10 2 6 11 4 13 第四主成分:11 6 5 7 10 13 12 9 1 8 3 2 4表(26)各行业按主成分得分进行排序结果图(21)主成分碎石图图(22)第一主成分与第二主成分下的散点图习题程序与结论:> industry<-data.frame(+X1=c(90342,4903,6735,49454,139190,12215,2372,11062,17111,1206,2150,5251,14341),+X2=c(52455,1973,21139,36241,203505,16219,6572,23078,23907,3930,5704,6155,13203),+X3=c(101091,2035,3767,81557,215898,10351,8103,54935,52108,6126,6200,10383,19396),+X4=c(19272,10313,1780,22504,10609,6382,12329,23804,21796,15586,10870,16875,14691),+ X5=c(82.0,34.2,36.1,98.1,93.2,62.5,184.4,370.4,221.5,330.4,184.2,146.4,94.6),+ X6=c(16.1,7.1,8.2,25.9,12.6,8.7,22.2,41.0,21.5,29.5,12.0,27.5,17.8),+X7=c(197435,592077,726396,348226,139572,145818,20921,65486,63806,1840,8913,78796,6354), +X8=c(0.172,0.003,0.003,0.985,0.628,0.066,0.152,0.263,0.276,0.437,0.274,0.151,1.574) )> industry.pr<-princomp(industry,cor=T)> summary(industry.pr) ####做主成分分析,得到4个主成分,累积贡献率达94.68% Importance of components:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 1.7620762 1.7021873 0.9644768 0.80132532 0.55143824Proportion of Variance 0.3881141 0.3621802 0.1162769 0.08026528 0.03801052Cumulative Proportion 0.3881141 0.7502943 0.8665712 0.94683649 0.98484701Comp.6 Comp.7 Comp.8Standard deviation 0.29427497 0.179400062 0.0494143207Proportion of Variance 0.01082472 0.004023048 0.0003052219Cumulative Proportion 0.99567173 0.999694778 1.0000000000> load<-loadings(industry.pr) ####求出载荷矩阵> loadLoadings:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8X1 -0.477 -0.296 -0.104 0.184 0.758 0.245X2 -0.473 -0.278 -0.163 -0.174 -0.305 -0.518 0.527X3 -0.424 -0.378 -0.156 -0.174 -0.781X4 0.213 -0.451 0.516 0.539 0.288 -0.249 0.220X5 0.388 -0.331 -0.321 -0.199 -0.450 0.582 0.233X6 0.352 -0.403 -0.145 0.279 -0.317 -0.714X7 -0.215 0.377 -0.140 0.758 -0.418 0.194X8 -0.273 0.891 -0.322 0.122Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125Cumulative Var 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000> plot(load[,1:2])> text(load[,1],load[,2],adj=c(-0.4,-0.3))> screeplot(industry.pr,npcs=4,type="lines") ####得出主成分的碎石图> biplot(industry.pr) ####得出在第一,第二主成分之下的散点图> p<-predict(industry.pr) ####预测数据,讲预测值放入p中> order(p[,1]);order(p[,2]);order(p[,3]);order(p[,4]);####将预测值分别以第一,第二,第三,第四主成分进行排序[1] 5 1 3 2 4 6 13 11 9 7 12 10 8[1] 5 8 4 9 10 1 13 12 7 11 6 2 3[1] 8 1 5 3 9 12 7 10 2 6 11 4 13[1] 11 6 5 7 10 13 12 9 1 8 3 2 4> kmeans(scale(p),4) ####将预测值进行标准化,并分为4类K-means clustering with 4 clusters of sizes 5, 1, 4, 3Cluster means:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.61 0.5132590 -0.03438438 -0.3405983 -0.5130031 0.2355151 0.224410402 -2.5699693 -1.32913757 -0.4848689 -0.9460127 -0.9000187 -0.064979503 0.2381581 0.72871986 -0.2995918 0.3126036 -0.4744091 -0.197097104 -0.3163193 -0.47127333 1.1287426 0.7535380 0.5400265 -0.08956137Comp.7 Comp.81 -0.38197798 -0.74748552 -0.67500209 0.45695483 0.09063069 0.98269154 0.74078975 -0.2167643Clustering vector:[1] 4 3 3 4 2 1 1 1 1 3 1 3 4Within cluster sum of squares by cluster:[1] 19.41137 0.00000 24.49504 16.61172(between_SS / total_SS = 37.0 %)Available components:[1] "cluster" "centers" "totss" "withinss" "tot.withinss"[6] "betweenss" "size"9.2####用数据框的形式输入数据####用数据框的形式输入数据sale<-data.frame(X1=c(82.9,88.0,99.9,105.3,117.7,131.0,148.2,161.8,174.2,184.7),X2=c(92,93,96,94,100,101,105,112,112,112),X3=c(17.1,21.3,25.1,29.0,34.0,40.0,44.0,49.0,51.0,53.0),X4=c(94,96,97,97,100,101,104,109,111,111),Y=c(8.4,9.6,10.4,11.4,12.2,14.2,15.8,17.9,19.6,20.8))####作线性回归lm.sol<-lm(Y~X1+X2+X3+X4,data=sale)summary(lm.sol)显示结果Call:lm(formula = Y ~ X1 + X2 + X3 + X4, data = sale)Residuals:1 2 3 4 5 6 70.024803 0.079476 0.012381 -0.007025 -0.288345 0.216090 -0.1420858 9 100.158360 -0.135964 0.082310Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) -17.66768 5.94360 -2.973 0.03107 *X1 0.09006 0.02095 4.298 0.00773 **X2 -0.23132 0.07132 -3.243 0.02287 *X3 0.01806 0.03907 0.462 0.66328X4 0.42075 0.11847 3.552 0.01636 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.2037 on 5 degrees of freedomMultiple R-squared: 0.9988, Adjusted R-squared: 0.9978F-statistic: 1021 on 4 and 5 DF, p-value: 1.827e-07模型通过t检验和F检验,因此回归方程为:Y=-17.66768+0.09006X1-0.23132X2+0.01806X3+0.42075X4 Y 是销售量,X1是居民可支配收入X2是该类消费品平均价格指数,X1和X2越高Y越高这与实际情况不符,原因是4个变量存在多重共线性,对变量作主成分回归,先作主成分分析。
统计建模与R软件第六章习题(回归分析)下篇
统计建模与R软件第六章习题(回归分析)下篇#6.6剖腹产后感染> data<-data.frame(+ x1=rep(c(1,1,1,1,1,0,0,0,0,0,0,0),c(1,17,2,11,87,28,30,23,3,9,8,32)),+ x2=rep(c(1,1,0,1,1,1,1,1,1,0,0,0),c(1,17,2,11,87,28,30,23,3,9,8,32)),+ x3=rep(c(1,1,1,0,0,1,1,0,0,0,1,1),c(1,17,2,11,87,28,30,23,3,9,8,32)),+ y=rep(c(1,0,0,1,0,1,0,1,0,0,1,0),c(1,17,2,11,87,28,30,23,3,9,8,32))+ )#x1表⽰是否使⽤了抗⽣素,x2表⽰是否有危险因⼦,x3表⽰是否有事先准备剖腹产,y表⽰是否有感染> glm.sol<-glm(y~x1+x2+x3,family=binomial,data=data)> summary(glm.sol)Call:glm(formula = y ~ x1 + x2 + x3, family = binomial, data = data)Deviance Residuals:Min 1Q Median 3Q Max-1.7149 -0.5298 -0.4933 0.7227 2.5141Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -0.8207 0.4947 -1.659 0.0971 .x1 -3.2544 0.4813 -6.761 1.37e-11 ***x2 2.0299 0.4553 4.459 8.25e-06 ***x3 -1.0720 0.4254 -2.520 0.0117 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 299.01 on 250 degrees of freedomResidual deviance: 226.52 on 247 degrees of freedomAIC: 234.52Number of Fisher Scoring iterations: 5#logistic回归模型为:P=exp(-0.8207-3.2544x1+2.0299x2-1.0720x3)/(1+exp(-0.8207-3.2544x1+2.0299x2-1.0720x3))#⽤上述回归模型作观测,使⽤抗⽣素,⽆危险因⼦,有剖腹产计划的感染概率为:> pre<-predict(glm.sol,data.frame(x1=1,x2=0,x3=1))> p<-exp(pre)/(1+exp(pre));p10.005783046 #感染概率为0.578%,⾮常⾮常低#不使⽤抗⽣素,有危险因⼦,⽆剖腹产计划的感染率为:> pre1<-predict(glm.sol,data.frame(x1=0,x2=1,x3=0))> p<-exp(pre1)/(1+exp(pre1));p10.7701639 #感染概率为77%,相当⾼了,回归⽅程中是否使⽤抗⽣素的回归系数绝对值最⼤,这个剖腹产还是要使⽤抗⽣素的=_=#6.7肺癌病⼈> cancer<-data.frame(+ x1=c(70,60,70,40,40,70,70,80,60,30,80,40,60,40,20,50,50,40,80,70,+ 60,90,50,70,20,80,60,50,70,40,30,30,40,60,80,70,30,60,80,70),+ x2=c(64,63,65,69,63,48,48,63,63,53,43,55,66,67,61,63,66,68,41,53,+ 37,54,52,50,65,52,70,40,26,44,54,59,69,50,62,68,39,49,64,67),+ x3=c(5,9,11,10,58,9,11,4,14,4,12,2,25,23,19,4,16,12,12,8,+ 13,12,8,7,21,28,13,13,22,36,9,87,5,22,4,15,4,11,10,18),+ x4=rep(c(1,2,3,0,1,2,3,0),c(7,7,2,4,8,4,3,5)),+ x5=rep(c(1,0),c(21,19)),+ y=rep(c(1,0,1,0,1,0,1,0,1,0,1),c(1,11,1,5,2,1,3,1,1,12,2))+ )#x1表⽰⽣活⾏动能⼒评分,x2表⽰病⼈年龄,x3表⽰由诊断到直⼊研究时间(⽉),x4表⽰肿瘤类型(0是磷癌,1是⼩型细胞癌,2是腺癌,3是⼤型细胞癌),x5表⽰两种化疗⽅法(1是常规,0是试验新法),y表⽰病⼈的⽣存时间(0是⽣存时间短,⼩于200天;1是⽣存时间长,⼤于等于200天)6.7.1logistic回归模型> glm.c<-glm(y~x1+x2+x3+x4+x5,family=binomial,data=cancer)> summary(glm.c)Call:glm(formula = y ~ x1 + x2 + x3 + x4 + x5, family = binomial,data = cancer)Deviance Residuals:Min 1Q Median 3Q Max-1.73402 -0.64522 -0.22773 0.09734 2.21941Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -7.22297 4.39926 -1.642 0.1006x1 0.10048 0.04322 2.325 0.0201 *x2 0.01716 0.04417 0.389 0.6976x3 0.01828 0.05405 0.338 0.7353x4 -1.07233 0.58653 -1.828 0.0675 .x5 -0.62473 0.96204 -0.649 0.5161---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 44.987 on 39 degrees of freedomResidual deviance: 28.329 on 34 degrees of freedomAIC: 40.329Number of Fisher Scoring iterations: 6#模型不显著,肿瘤类型x4与化疗⽅法x5应该是主要的影响因素#计算各病⼈⽣存时间⼤于等于200天的概率估计值> n=c(cancer$y==1)> n[1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE[14] FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE[27] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE[40] TRUE> pre=glm.c$fitted[n];pre1 13 19 20 22 23 240.33257413 0.08518899 0.75283057 0.56015746 0.86922215 0.09686322 0.4315437926 39 400.75903351 0.89058079 0.78407564> p=exp(pre)/(1+exp(pre));p1 13 19 20 22 23 24 260.5823856 0.5212844 0.6797952 0.6364890 0.7045838 0.5241969 0.6062422 0.681143939 400.7090100 0.68655786.7.2逐步回归法选取⾃变量> ew<-step(glm.c)Start: AIC=40.33y ~ x1 + x2 + x3 + x4 + x5Df Deviance AIC- x3 1 28.430 38.430- x2 1 28.484 38.484- x5 1 28.750 38.75028.329 40.329- x4 1 32.483 42.483- x1 1 38.297 48.297Step: AIC=38.43y ~ x1 + x2 + x4 + x5Df Deviance AIC- x2 1 28.564 36.564- x5 1 28.955 36.95528.430 38.430- x4 1 32.563 40.563- x1 1 38.474 46.474Step: AIC=36.56y ~ x1 + x4 + x5Df Deviance AIC- x5 1 29.073 35.07328.564 36.564- x4 1 32.892 38.892- x1 1 38.478 44.478Step: AIC=35.07y ~ x1 + x4Df Deviance AIC29.073 35.073- x4 1 33.535 37.535- x1 1 39.131 43.131#剩下x1与x4⾃变量> summary(ew)Call:glm(formula = y ~ x1 + x4, family = binomial, data = cancer) Deviance Residuals:Min 1Q Median 3Q Max-1.4825 -0.6617 -0.1877 0.1227 2.2844Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) -6.13755 2.73844 -2.241 0.0250 *x1 0.09759 0.04079 2.393 0.0167 *x4 -1.12524 0.60239 -1.868 0.0618 .---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1)Null deviance: 44.987 on 39 degrees of freedomResidual deviance: 29.073 on 37 degrees of freedomAIC: 35.073Number of Fisher Scoring iterations: 6#模型显著性提⾼#计算逐步回归后新模型各病⼈⽣存时间⼤于等于200天的概率估计值> pre1=ew$fitted[n];pre11 13 19 20 22 23 240.39371565 0.07358866 0.84149409 0.66674959 0.82054044 0.08444317 0.3937156526 39 400.63277649 0.84149409 0.66674959> p1=exp(pre1)/(1+exp(pre1));p11 13 19 20 22 23 24 260.5971768 0.5183889 0.6987798 0.6607750 0.6943510 0.5210983 0.5971768 0.653118839 400.6987798 0.6607750#6.8⾃动咖啡售货机与咖啡销售> x<-c(0,0,1,1,2,2,3,3,4,4,5,5,6,6) #咖啡机数量> y<-c(508.1,498.4,568.2,577.3,651.7,657,713.4,697.5,755.3,758.9,787.6,792.1,841.4,831.8) #咖啡销量6.8.1线性回归模型> lm.coffee<-lm(y~x)> summary(lm.coffee)Call:lm(formula = y ~ x)Residuals:Min 1Q Median 3Q Max-25.400 -11.484 -3.779 14.629 24.921Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 523.800 8.474 61.81 < 2e-16 ***x 54.893 2.350 23.36 2.26e-11 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 17.59 on 12 degrees of freedomMultiple R-squared: 0.9785, Adjusted R-squared: 0.9767F-statistic: 545.5 on 1 and 12 DF, p-value: 2.265e-11#线性回归⽅程为y=523.800+54.893x6.8.2多项式回归模型> plot(x,y)#从散点图上看,⽤⼆次曲线拟合更好> lm.coffee1<-lm(y~x+I(x^2))> summary(lm.coffee1)Call:lm(formula = y ~ x + I(x^2))Residuals:Min 1Q Median 3Q Max-10.6643 -5.6622 -0.4655 5.5000 10.6679Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 502.5560 4.8500 103.619 < 2e-16 ***x 80.3857 3.7861 21.232 2.81e-10 ***I(x^2) -4.2488 0.6063 -7.008 2.25e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 7.858 on 11 degrees of freedomMultiple R-squared: 0.9961, Adjusted R-squared: 0.9953F-statistic: 1391 on 2 and 11 DF, p-value: 5.948e-14#多项式回归⽅程为:y=502.5560+80.3857x-4.2488x^2 相⽐线性回归,残差标准误降低,R平⽅达到0.99616.8.3画出散点图与拟合曲线> plot(x,y)> xfit<-seq(0,6,0.01)> yfit<-predict(lm.coffee1,data.frame(x=xfit))> lines(xfit,yfit)#6.9重伤病⼈出院后恢复情况> x<-c(2,5,7,10,14,19,26,31,34,38,45,52,53,60,65) #住院天数> y<-c(54,50,45,37,35,25,20,16,18,13,8,11,8,4,6) #预后指数6.9.1内在线性模型⽅法> plot(x,y)> lm.patient<-lm(log(y)~x)> summary(lm.patient)Call:lm(formula = log(y) ~ x)Residuals:Min 1Q Median 3Q Max-0.37241 -0.07073 0.02777 0.05982 0.33539Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 4.037159 0.084103 48.00 5.08e-16 ***x -0.037974 0.002284 -16.62 3.86e-10 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.1794 on 13 degrees of freedomMultiple R-squared: 0.9551, Adjusted R-squared: 0.9516F-statistic: 276.4 on 1 and 13 DF, p-value: 3.858e-10#线性回归⽅程为:lny=4.037159-0.037974x#计算估计值> yev=exp(4.037159)*exp(-0.037974*x);yev[1] 52.520890 46.865838 43.438277 38.761171 33.298855 27.540372 21.111861 17.460916 [9] 15.580856 13.385165 10.260781 7.865696 7.572604 5.804996 4.8011196.9.2⽤⾮线性⽅法计算估计值> nls.patient<-nls(y~a*exp(b*x),start=list(a=10,b=-0.01))> summary(nls.patient)Formula: y ~ a * exp(b * x)Parameters:Estimate Std. Error t value Pr(>|t|)a 58.606568 1.472160 39.81 5.70e-15 ***b -0.039586 0.001711 -23.13 6.01e-12 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.951 on 13 degrees of freedomNumber of iterations to convergence: 6Achieved convergence tolerance: 5.107e-07#⾮线性模型为:y=58.606568 * exp(-0.039586 * x)#计算估计值> yev=58.606568 * exp(-0.039586 * x);yev[1] 54.145495 48.082427 44.422441 39.448135 33.671197 27.624768 20.938944 17.178881 [9] 15.255236 13.021200 9.869772 7.481062 7.190702 5.450388 4.471647#使⽤nlm( )函数> fn<-function(p,X,Y)+ {+ f<-Y-p[1]*exp(p[2]*X)+ res<-sum(f^2)+ f1<-exp(p[2]*X)+ f2<-p[1]*exp(p[2]*X)*X+ J<-cbind(f1,f2)+ res+ }> nlm(fn,p=c(50, -0.01),x,y,hessian=TRUE)$minimum[1] 49.4593$estimate[1] 58.60712251 -0.03958742$gradient[1] -1.774625e-05 -9.522651e-03$hessian[,1] [,2][1,] 7.022561 4285.646[2,] 4285.646407 5161077.388$code[1] 2$iterations[1] 33#计算估计值> yev=58.60712251*exp(-0.03958742*x);yev[1] 54.145853 48.082541 44.422420 39.447948 33.670846 27.624284 20.938369 17.178287 [9] 15.254644 13.020620 9.869235 7.480580 7.190229 5.449975 4.471276#来看下估计值与实际值的差> y-yev[1] -0.14585337 1.91745922 0.57758014 -2.44794812 1.32915418 -2.62428435[7] -0.93836901 -1.17828719 2.74535642 -0.02062025 -1.86923491 3.51941952[13] 0.80977134 -1.44997507 1.52872352。
统计模拟与R相关资料习题答案
第一篇:R介绍
– R是
• 一个开放(GPL)的统计编程环境 • 一种语言,是S语言(由AT&T Bell实验室的Rick Becker, John Chambers,Allan Wilks开发)的一 种方言(dialect) 之一,另一则为S-plus. • 一种软件,是集统计分析与图形直观显示于一体的统计 分析
• 参考资料 随软件所附pdf文档(help->manuals),随版 本更新: – W.N. Venables, D.M. Smith and the R DCT: Introduction to R -- Notes on R: A Programming Environment for Data Analysis and Graphics, 2003. /R web/Rnotes/R.html – R DCT, The R Environment for Statistical Computing and Graphics -- Reference Index,2003. – R DCT, R Data Import/Export, 2003. – R DCT, R Language Definition,2003 – R DCT, Writing R Extensions,2003
– R作为一个计划(project),最早(1995年)是由 Auckland大学统计系的Robert Gentleman 和Ross Ihaka开始编制,目前由R核心开发小 组(R Development Core Team – 以后用R DCT表示)维护,他们完全自愿、工作努力负责, 并将全球优秀的统计应用软件打包提供给我们。我 们可以通过R计划的网站()了解有关R的最新信息和使用说明, 得到最新版本的R软件和基于R的应用统计软件包 .
统计建模与R软件假设检验习题答案
04 习题答案解析
习题一答案解析
答案:D
解析:根据题目描述,A、B、C三个选项都是描述数据特征的,而D选项是描述数据来源的,与题目 要求的“数据特征”不符。
习题二答案解析
答案:B
解析:根据题目要求,需要选择一个假设检验的方法。A选项是参数检验,适用于总体分布已知的情况;B选项是非参数检验 ,适用于总体分布未知或不符合正态分布的情况;C选项是回归分析,用于研究变量之间的关系;D选项是聚类分析,用于数 据的分类。根据题目描述,由于总体分布未知且不符合正态分布,所以选择B选项。
模型评估
01
02
03
交叉验证
将数据集分成训练集和测 试集,使用训练集训练模 型,在测试集上评估模型 的性能。
均方误差
衡量预测值与实际值之间 的误差,越小越好。
准确率
衡量分类模型正确预测的 比例,越高越好。
02 R软件基础
R软件介绍
总结词
R软件是一种开源的统计计算和图形绘制软件,广泛应用于数据分析和统计建模 。
解析:根据题目要求,需要选择一个统计量 来描述数据的集中趋势。A选项是平均数, 是最常用的描述集中趋势的统计量;B选项 是中位数,描述数据的中位数位置;C选项 是众数,描述数据中出现次数最多的数;D 选项是标准差,描述数据的离散程度。根据
题目要求,选择A选项。
习题五答案解析
答案:C
解析:根据题目要求,需要选择一个统计量来检验两 个样本是否来自同一个总体。A选项是t检验,适用于 两个正态分布的总体;B选项是卡方检验,适用于分 类数据的比较;C选项是F检验,适用于两个总体方差 的比较;D选项是z检验,适用于总体比例的比较。根 据题目要求,选择C选项。
假设检验的优缺点
统计建模与R软件课后答案
第二章2.1> x<-c(1,2,3);y<-c(4,5,6)> e<-c(1,1,1)> z<-2*x+y+e;z[1] 7 10 13> z1<-crossprod(x,y);z1[,1][1,] 32> z2<-outer(x,y);z2[,1] [,2] [,3][1,] 4 5 6[2,] 8 10 12[3,] 12 15 182.2(1) >A<-matrix(1:20,nrow=4);B<-matrix(1:20,nrow=4,byrow=T) > C<-A+B;C(2)> D<-A%*%B;D(3)> E<-A*B;E(4)> F<-A[1:3,1:3](5)> G<-B[,-3]2.3> x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x2.4> H<-matrix(nrow=5,ncol=5)> for (i in 1:5)+ for(j in 1:5)+ H[i,j]<-1/(i+j-1)(1)> det(H)(2)> solve(H)(3)> eigen(H)2.5> studentdata<-data.frame(姓名=c('张三','李四','王五','赵六','丁一')+ ,性别=c('女','男','女','男','女'),年龄=c('14','15','16','14','15'),+ 身高=c('156','165','157','162','159'),体重=c('42','49','41.5','52','45.5'))2.6> write.table(studentdata,file='student.txt')> write.csv(studentdata,file='student.csv')2.7count<-function(n){if (n<=0)print('要求输入一个正整数')else{repeat{if (n%%2==0)n<-n/2elsen<-(3*n+1)if(n==1)break}print('运算成功')}}第三章3.1首先将数据录入为x。
统计建模与r软件 课后习题答案
统计建模与r软件课后习题答案统计建模与R软件课后习题答案在统计建模与R软件课程中,学生们经常需要完成一系列的习题来巩固所学知识。
这些习题涉及到统计建模的理论和实践,以及如何使用R软件来进行数据分析和建模。
在本文中,我们将给出一些常见的统计建模与R软件课后习题的答案,希望能够帮助学生更好地理解课程内容。
1. 线性回归模型习题:使用R软件对给定数据集进行线性回归分析,并给出回归方程和相关系数。
答案:在R软件中,可以使用lm()函数来进行线性回归分析。
例如,对于数据集data,可以使用以下代码进行线性回归分析:```model <- lm(y ~ x, data=data)summary(model)```其中,y和x分别表示因变量和自变量。
通过summary()函数可以得到回归方程和相关系数等信息。
2. 逻辑回归模型习题:使用R软件对给定数据集进行逻辑回归分析,并给出回归方程和模型拟合度。
答案:逻辑回归分析可以使用glm()函数来进行。
例如,对于数据集data,可以使用以下代码进行逻辑回归分析:```model <- glm(y ~ x, data=data, family=binomial)summary(model)```其中,y和x分别表示因变量和自变量,family=binomial表示使用二项分布进行逻辑回归分析。
通过summary()函数可以得到回归方程和模型拟合度等信息。
3. 方差分析习题:使用R软件对给定数据集进行方差分析,并给出各组之间的差异是否显著。
答案:在R软件中,可以使用aov()函数来进行方差分析。
例如,对于数据集data,可以使用以下代码进行方差分析:```model <- aov(y ~ group, data=data)summary(model)```其中,y和group分别表示因变量和自变量。
通过summary()函数可以得到各组之间的差异是否显著等信息。
统计建模与R软件第七章习题答案
统计建模与R软件第七章习题答案统计建模与R软件第七章习题答案(方差分析)Ex7.1(1)>lamp<-data.frame(X=c(115,116,98,83,103,107,118,116,73,89,85,97),A=f a ctor(rep(1:3,c(4,4,4))))> lamp.aov<-aov(X~A,data=lamp);summary(lamp.aov)Df Sum Sq Mean Sq F value Pr(>F)A 2 1304 652.0 4.923 0.0359 *Residuals 9 1192 132.4P值小于0.05,有显著差异。
(2)对甲的区间估计:> a<-c(115,116,98,83)> t.test(a)One Sample t-testdata: at = 13.1341, df = 3, p-value = 0.0009534alternative hypothesis: true mean is not equal to 095 percent confidence interval:78.04264 127.95736sample estimates:mean of x103或者用这个命令更简单:>attach(lamp)> t.test(X[A==1])乙的均值估计为111,95%置信区间为99.59932, 122.40068。
丙的均值估计为86,95%置信区间为70.08777, 101.91223。
(3)多重检验:> attach(lamp)P值不做调整:> pairwise.t.test(X,A,p.adjust.method = "none")Pairwise comparisons using t tests with pooled SDdata: X and A1 22 0.351 -3 0.066 0.013P值进行Holm调整:P value adjustment method: none> pairwise.t.test(X,A,p.adjust.method = "holm",data)Pairwise comparisons using t tests with pooled SDdata: X and A1 22 0.35 -3 0.13 0.04P value adjustment method: holm不论采取哪种方法,都可看出乙和丙有显著差异。
《统计建模与R软件》书本课后习题答案
第二章答案:Ex2.1x<-c(1,2,3)y<-c(4,5,6)e<-c(1,1,1)z=2*x+y+ez1=crossprod(x,y)#z1为x1与x2的内积或者x%*%yz2=tcrossprod(x,y)#z1为x1与x2的外积或者x%o%yz;z1;z2要点:基本的列表赋值方法,内积和外积概念。
内积为标量,外积为矩阵。
Ex2.2A<-matrix(1:20,c(4,5));AB<-matrix(1:20,nrow=4,byrow=TRUE);BC=A+B;C#不存在AB这种写法E=A*B;EF<-A[1:3,1:3];FH<-matrix(c(1,2,4,5),nrow=1);H#H起过渡作用,不规则的数组下标G<-B[,H];G要点:矩阵赋值方法。
默认是byrow=FALSE,数据按列放置。
取出部分数据的方法。
可以用数组作为数组的下标取出数组元素。
Ex2.3x<-c(rep(1,times=5),rep(2,times=3),rep(3,times=4),rep(4,times=2));x #或者省略times=,如下面的形式x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x要点:rep()的使用方法。
rep(a,b)即将a重复b次Ex2.4n <- 5; H<-array(0,dim=c(n,n))for (i in 1:n){for (j in 1:n){H[i,j]<-1/(i+j-1)}};HG <- solve(H);G #求H的逆矩阵ev <- eigen(H);ev #求H的特征值和特征向量要点:数组初始化;for循环的使用待解决:如何将很长的命令(如for循环)用几行打出来再执行?每次想换行的时候一按回车就执行了还没打完的命令...Ex2.5StudentData<-data.frame(name=c("zhangsan","lisi","wangwu","zhaoliu","dingyi"),sex=c("F","M", "F","M","F"),age=c("14","15","16","14","15"),height=c("156","165","157","162","159"),weight=c( "42","49","41.5","52","45.5"));StudentData要点:数据框的使用待解决:SSH登陆linux服务器中文显示乱码。
统计建模与R软件课后答案修订稿
统计建模与R软件课后答案文档编制序号:[KKIDT-LLE0828-LLETD298-POI08]第二章> x<-c(1,2,3);y<-c(4,5,6)> e<-c(1,1,1)> z<-2*x+y+e;z[1] 7 10 13> z1<-crossprod(x,y);z1[,1][1,] 32> z2<-outer(x,y);z2[,1] [,2] [,3][1,] 4 5 6[2,] 8 10 12[3,] 12 15 18(1) > A<-matrix(1:20,nrow=4);B<-matrix(1:20,nrow=4,byrow=T) > C<-A+B;C(2)> D<-A%*%B;D(3)> E<-A*B;E(4)> F<-A[1:3,1:3](5)> G<-B[,-3]> x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x> H<-matrix(nrow=5,ncol=5)> for (i in 1:5)+ for(j in 1:5)+ H[i,j]<-1/(i+j-1)(1)> det(H)(2)> solve(H)(3)> eigen(H)> studentdata<(姓名=c('张三','李四','王五','赵六','丁一')+ ,性别=c('女','男','女','男','女'),年龄=c('14','15','16','14','15'), + 身高=c('156','165','157','162','159'),体重=c('42','49','','52',''))> (studentdata,file='')> (studentdata,file='')count<-function(n){if (n<=0)print('要求输入一个正整数')else{repeat{if (n%%2==0)n<-n/2elsen<-(3*n+1)if(n==1)break}print('运算成功')}}第三章首先将数据录入为x。
R软件课后习题第五章
第五章5.1####写出求正态总体均值检验的R程序(程序名:mean.test1.R)mean.test1<-function(x, mu=0, sigma=-1, side=0){source("P_value.R")n<-length(x); xb<-mean(x)if (sigma>0){z<-(xb-mu)/(sigma/sqrt(n))P<-P_value(pnorm, z, side=side)data.frame(mean=xb, df=n, Z=z, P_value=P)}else{t<-(xb-mu)/(sd(x)/sqrt(n))P<-P_value(pt, t, paramet=n-1, side=side)data.frame(mean=xb, df=n-1, T=t, P_value=P)}}####写出求P值的R程序(程序名:P_value.R)P_value<-function(cdf, x, paramet=numeric(0), side=0){n<-length(paramet)P<-switch(n+1,cdf(x),cdf(x, paramet),cdf(x, paramet[1], paramet[2]),cdf(x, paramet[1], paramet[2], paramet[3]))if (side<0) Pelse if (side>0) 1-Pelseif (P<1/2) 2*Pelse 2*(1-P)}####输入数据,再调用函数mean.test1()>x<-c(220,188,162,230,145,160,238,188,247,113,126,245,164,231,256,183,190,158,224,175) > source("mean.test1.R")> a<-mean.test1(x, mu=225,side=0)> a得到:mean df T P_value1 192.15 19 -3.478262 0.002516436可知,P值小于0.05,故与正常值存在差异5.2####输入数据,再调用函数mean.test1()> x<-c(1067,919,1196,785,1126,936,918,1156,920,948)> source("mean.test1.R")> mean.test1(x, mu=1000,side=1)得到:mean df T P_value1 997.1 9 -0.06971322 0.5270268所以灯泡寿命为1000小时以上的概率是0.47297325.3####写出两总体均值检验的R程序(程序名:mean.test2.R)mean.test2<-function(x, y,sigma=c(-1, -1), var.equal=FALSE, side=0){source("P_value.R")n1<-length(x); n2<-length(y)xb<-mean(x); yb<-mean(y)if (all(sigma>0)){z<-(xb-yb)/sqrt(sigma[1]^2/n1+sigma[2]^2/n2)P<-P_value(pnorm, z, side=side)data.frame(mean=xb-yb, df=n1+n2, Z=z, P_value=P)}else{if (var.equal == TRUE){Sw<-sqrt(((n1-1)*var(x)+(n2-1)*var(y))/(n1+n2-2))t<-(xb-yb)/(Sw*sqrt(1/n1+1/n2))nu<-n1+n2-2}else{S1<-var(x); S2<-var(y)nu<-(S1/n1+S2/n2)^2/(S1^2/n1^2/(n1-1)+S2^2/n2^2/(n2-1))t<-(xb-yb)/sqrt(S1/n1+S2/n2)}P<-P_value(pt, t, paramet=nu, side=side)data.frame(mean=xb-yb, df=nu, T=t, P_value=P)}}####输入数据,再调用函数mean.test2()> x<-c(113,120,138,120,100,118,138,123)> y<-c(138,116,125,136,110,132,130,110)> source("mean.test2.R")> mean.test2(x, y, var.equal=TRUE, side=0)得到:mean df T P_value1 -3.375 14 -0.5659672 0.5803752P值大于0.05,故接受原假设5.4####写出均值已知和均值未知两种情况方差比检验的R程序(程序名:var.test2.R)var.test2<-function(x, y,mu=c(Inf,Inf),side=0){source("P_value.R")n1<-length(x); n2<-length(y)if (all(all(mu<Inf)){Sx2<-sum((x-mu[1])^2)/n1;Sy2<-sum((y-mu[2])^2)/n2df1=n1;df2=n2}else{Sx2<-var(x); Sy2<-var(y);df1=n1-1;df2=n2-1}r<-Sx2/Sy2P<-P_value(pf, r, paramet=c(df1,df2), side=side)data.frame(rate=r, df1=df1, df2=df2,F=f, P_value=P)}}####输入数据>x<-c(-0.70,-5.60,2.00,2.80,0.70,3.50,4.00,5.80,7.10,-0.50,2.50,-1.60,1.70,3.00,0.40,4.50,4.60,2.5 0,6.00,-1.40)> a<-shapiro.test(x)> aShapiro-Wilk normality testdata: xW = 0.9699, p-value = 0.7527>0.05>y<-c(3.70,6.50,5.00,5.20,0.80,0.20,0.60,3.40,6.60,-1.10,6.00,3.80,2.00,1.60,2.00,2.20,1.20,3.10, 1.70,-2.00)> b<-shapiro.test(y)> bShapiro-Wilk normality testdata: yW = 0.971, p-value = 0.7754>0.05由以上可知,两组数据均为正态分布####输入数据,再调用函数mean.test2()>x<-c(-0.70,-5.60,2.00,2.80,0.70,3.50,4.00,5.80,7.10,-0.50,2.50,-1.60,1.70,3.00,0.40,4.50,4.60,2.5 0,6.00,-1.40)>y<-c(3.70,6.50,5.00,5.20,0.80,0.20,0.60,3.40,6.60,-1.10,6.00,3.80,2.00,1.60,2.00,2.20,1.20,3.10, 1.70,-2.00)> source("mean.test2.R")> a<-mean.test2(x, y, var.equal=TRUE, side=0);amean df T P_value1 -0.56 38 -0.641872 0.5248097> b<-mean.test2(x, y, var.equal=FALSE, side=0);bmean df T P_value1 -0.56 36.08632 -0.641872 0.525013> c<-t.test(x-y, alternative = "two.sided");cOne Sample t-testdata: x - yt = -0.6464, df = 19, p-value = 0.5257alternative hypothesis: true mean is not equal to 095 percent confidence interval:-2.373146 1.253146sample estimates:mean of x-0.56以上P值均大于0.05,故均值无差异。
统计建模与R软件 习题6.4 答案
Homework6.4 P332解:这个过程中回归R语言计算过程如下:toothpaste<-data.frame(X1<-c(-0.05,0.25,0.60,0,0.25,0.2,0.15,0.05,-0.15,0.15,0.2,0.1,0.4,0.45,0.35,0.3,0.5,0.5,0.4,-0.05,-0 .05,-0.1,0.2,0.1,0.5,0.6,-0.05,0,0.05,0.55),X2<-c(5.5,6.75,7.25,5.5,7,6.5,6.75,5.25,5.25,6,6.5,6.25,7,6.9,6.8,6.8,7.1,7,6.8,6.5,6.25,6,6.5,7,6.8, 6.8,6.5,5.75,5.8,6.8),Y<-c(7.38,8.51,9.52,7.5,9.33,8.28,8.75,7.87,7.1,8,7.89,8.15,9.1,8.86,8.9,8.87,9.26,9,8.75,7.95,7.6 5,7.27,8,8.5,8.75,9.21,8.27,7.67,7.93,9.26))lm.sol<-lm(Y~X1+X2,data=toothpaste);summary(lm.sol)influence.measures(lm.sol)回归结果输出如下:Call:lm(formula = Y ~ X1 + X2, data = toothpaste)Residuals:Min 1Q Median 3Q Max-0.49779 -0.12031 -0.00867 0.11084 0.58106Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 4.4075 0.7223 6.102 1.62e-06 ***X1 1.5883 0.2994 5.304 1.35e-05 ***X2 0.5635 0.1191 4.733 6.25e-05 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.2383 on 27 degrees of freedomMultiple R-squared: 0.886, Adjusted R-squared: 0.8776F-statistic: 105 on 2 and 27 DF, p-value: 1.845e-13从上述输出的结果知道,模型通过t检验和F检验,因此可以得到回归模型:.4XXY++=407556352.015883.1回归诊断结果如下:Influence measures oflm(formula = Y ~ X1 + X2, data = toothpaste) :dfb.1_ dfb.X1 dfb.X2 dffit cov.r cook.d hat inf1 -0.05536 -6.96e-03 0.050306 -0.0807 1.281 2.25e-03 0.13012 0.04401 2.83e-02 -0.048123 -0.0922 1.152 2.92e-03 0.04703 -0.01269 6.49e-02 0.010195 0.1295 1.277 5.78e-03 0.13394 -0.00929 -2.98e-03 0.008653 -0.0118 1.299 4.81e-05 0.13805 -0.68922 -4.88e-01 0.719765 0.9115 0.536 2.18e-01 0.0909 *6 0.01121 1.59e-02 -0.016637 -0.0860 1.133 2.54e-03 0.03487 -0.26976 -2.67e-01 0.290196 0.3922 0.997 4.98e-02 0.07968 1.17309 6.39e-01 -1.129656 1.2670 0.895 4.69e-01 0.2493 *9 -0.03920 -5.27e-05 0.035418 -0.0601 1.373 1.25e-03 0.1860 *10 -0.02101 -1.08e-02 0.019442 -0.0295 1.194 3.02e-04 0.063611 0.05653 8.02e-02 -0.083875 -0.4338 0.670 5.42e-02 0.034812 0.00259 -1.74e-02 0.001766 0.0546 1.160 1.03e-03 0.041913 -0.04609 1.35e-02 0.047372 0.1279 1.167 5.61e-03 0.065614 -0.00181 -8.57e-02 0.001894 -0.1784 1.149 1.08e-02 0.070715 -0.01729 1.83e-02 0.019442 0.0994 1.149 3.39e-03 0.047616 -0.05590 -1.53e-02 0.060772 0.1449 1.118 7.15e-03 0.046517 -0.01349 3.00e-02 0.012894 0.0786 1.223 2.13e-03 0.090618 0.00123 -1.01e-01 0.000402 -0.1972 1.173 1.33e-02 0.088119 -0.00397 -5.62e-02 0.002754 -0.1299 1.149 5.78e-03 0.056620 0.04972 6.92e-02 -0.054542 -0.0782 1.320 2.11e-03 0.155121 0.11714 2.36e-01 -0.139343 -0.2970 1.142 2.97e-02 0.102022 0.08114 3.85e-01 -0.126203 -0.5601 0.929 9.84e-02 0.104123 0.04240 6.02e-02 -0.062902 -0.3253 0.841 3.29e-02 0.034824 0.01971 1.90e-02 -0.020667 -0.0235 1.378 1.90e-04 0.1873 *25 -0.13024 -3.05e-01 0.133977 -0.4177 1.038 5.69e-02 0.098226 0.01670 3.14e-02 -0.017393 0.0367 1.351 4.67e-04 0.1714 *27 -0.35156 -4.89e-01 0.385694 0.5528 1.100 9.94e-02 0.155128 0.01703 -1.01e-04 -0.014952 0.0296 1.223 3.03e-04 0.085529 0.14921 3.39e-02 -0.134903 0.2241 1.140 1.70e-02 0.080230 0.10055 2.05e-01 -0.104290 0.2545 1.227 2.21e-02 0.1309从上面诊断输出结果可以得到对回归模型存在影响的有样本点5、8、9、24、26,从dfb.1_角度看,点8对响应变量影响最大,从dffit cov.r 角度看,点24、26对自变量影响最大,那么考虑去掉样本点8、24、26除去影响点后整个R语言的计算过程如下:toothpaste<-data.frame( X1=c(-0.05,0.25,0.60,0,0.20,0.15,-0.15,0.15,0.10,0.40,0.45,0.35,0.30,0.5 0,0.50,0.40,-0.05,-0.05,-0.10,0.20,0.10,0.50,0.60,-0.05,0,0.05,0.55),X2=c(5.50,6.75,7.25,5.50,6.50,6.75,5.25,6.00,6.25,7.00,6.90,6.80,6.80,7.10,7.00,6.80,6.50,6.25,6. 00,6.50,7.00,6.80,6.80,6.50,5.75,5.80,6.80),Y=c(7.38,8.51,9.52,7.50,8.28,8.75,7.10,8.00,8.15,9.10,8.86,8.90,8.87,9.26,9.00,8.75,7.95,7.65,7.27,8.00,8.50,8.75,9.21,8.27,7.67,7.93,9.26))lm.sol<-lm(Y~X1+X2,data=toothpaste);summary(lm.sol)回归输出结果如下:Call:lm(formula = Y ~ X1 + X2, data = toothpaste)Residuals:Min 1Q Median 3Q Max-0.37130 -0.10114 0.03066 0.10016 0.30162Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 4.0759 0.6267 6.504 1.00e-06 ***X1 1.5276 0.2354 6.489 1.04e-06 ***X2 0.6138 0.1027 5.974 3.63e-06 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.1767 on 24 degrees of freedomMultiple R-squared: 0.9378, Adjusted R-squared: 0.9327F-statistic: 181 on 2 and 24 DF, p-value: 3.33e-15模型通过t 检验和F 检验,因此回归模型为:26138.015276.10759.4X X Y ++=∧。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
第二章答案:Ex2.1x<-c(1,2,3)y<-c(4,5,6)e<-c(1,1,1)z=2*x+y+ez1=crossprod(x,y)#z1为x1与x2的内积或者x%*%yz2=tcrossprod(x,y)#z1为x1与x2的外积或者x%o%yz;z1;z2要点:基本的列表赋值方法,内积和外积概念。
内积为标量,外积为矩阵。
Ex2.2A<-matrix(1:20,c(4,5));AB<-matrix(1:20,nrow=4,byrow=TRUE);BC=A+B;C#不存在AB这种写法E=A*B;EF<-A[1:3,1:3];FH<-matrix(c(1,2,4,5),nrow=1);H#H起过渡作用,不规则的数组下标G<-B[,H];G要点:矩阵赋值方法。
默认是byrow=FALSE,数据按列放置。
取出部分数据的方法。
可以用数组作为数组的下标取出数组元素。
Ex2.3x<-c(rep(1,times=5),rep(2,times=3),rep(3,times=4),rep(4,times=2));x #或者省略times=,如下面的形式x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x要点:rep()的使用方法。
rep(a,b)即将a重复b次Ex2.4n <- 5; H<-array(0,dim=c(n,n))for (i in 1:n){for (j in 1:n){H[i,j]<-1/(i+j-1)}};HG <- solve(H);G #求H的逆矩阵ev <- eigen(H);ev #求H的特征值和特征向量要点:数组初始化;for循环的使用待解决:如何将很长的命令(如for循环)用几行打出来再执行?每次想换行的时候一按回车就执行了还没打完的命令...Ex2.5StudentData<-data.frame(name=c("zhangsan","lisi","wangwu","zhaoliu","dingyi"),sex=c("F","M","F","M","F"),a ge=c("14","15","16","14","15"),height=c("156","165","157","162","159"),weight=c("42","49","41. 5","52","45.5"));StudentData要点:数据框的使用待解决:SSH登陆linux服务器中文显示乱码。
此处用英文代替。
Ex2.6write.table(StudentData,file="studentdata.txt")#把数据框StudentData在工作目录里输出,输出的文件名为studentdata.txt.StudentData_a<-read.table("studentdata.txt");StudentData_a#以数据框的形式读取文档studentdata.txt,存入数据框StudentData_a中。
write.csv(StudentData_a,"studentdata.csv")#把数据框StudentData_a在工作目录里输出,输出的文件名为studentdata.csv,可用Excel打开.要点:读写文件。
read.table("file")write.table(Rdata,"file")read.csv("file")write.csv(Rdata,"file")外部文件,不论是待读入或是要写出的,命令中都得加双引号。
Ex2.7Fun<-function(n){if(n <= 0)list(fail="please input a integer above 0!")else{repeat{if(n==1) breakelse if(n%%2==0){n<-n/2}else n<- 3*n+1}list("sucess!")}}在linux下新建一个R文件,输入上述代码,保存为"2.7.R"然后在当前目录下进入R环境,输入source("2.7.R"),即打开了这个程序脚本。
然后就可以执行函数了。
输入Fun(67),显示"sucess!"输入Fun(-1),显示$fail[1] "please input a integer above 0!"待解决:source("*.R")是可以理解为载入这个R文件吧?如何在R环境下关闭R文件呢?OK,自己写的第一个R程序~~Ex3.1新建txt文件如下:3.1.txt74.3 79.5 75.0 73.5 75.8 74.0 73.5 67.2 75.8 73.5 78.8 75.6 73.5 75.0 75.872.0 79.5 76.5 73.5 79.5 68.8 75.0 78.8 72.0 68.8 76.5 73.5 72.7 75.0 70.478.0 78.8 74.3 64.3 76.5 74.3 74.7 70.4 72.7 76.5 70.4 72.0 75.8 75.8 70.476.5 65.0 77.2 73.5 72.7 80.5 72.0 65.0 80.3 71.2 77.6 76.5 68.8 73.5 77.280.5 72.0 74.3 69.7 81.2 67.3 81.6 67.3 72.7 84.3 69.7 74.3 71.2 74.3 75.072.0 75.4 67.3 81.6 75.0 71.2 71.2 69.7 73.5 70.4 75.0 72.7 67.3 70.3 76.573.5 72.0 68.0 73.5 68.0 74.3 72.7 72.7 74.3 70.4编写一个函数(程序名为data_outline.R)描述样本的各种描述性统计量。
data_outline<-function(x){n<-length(x)m<-mean(x)v<-var(x)s<-sd(x)me<-median(x)cv<-100*s/mcss<-sum((x-m)^2)uss<-sum(x^2)R <- max(x)-min(x)R1 <-quantile(x,3/4)-quantile(x,1/4)sm <-s/sqrt(n)g1 <-n/((n-1)*(n-2))*sum((x-m)^3)/s^3g2 <-((n*(n+1))/((n-1)*(n-2)*(n-3))*sum((x-m)^4)/s^4-(3*(n-1)^2)/((n-2)*(n-3)))data.frame(N=n,Mean=m,Var=v,std_dev=s,Median=me,std_mean=sm,CV=cv,CSS=css,USS=uss,R= R,R1=R1,Skewness=g1,Kurtosis=g2,s=1)}进入R,source("data_outline.R") #将程序调入内存serumdata<-scan("3.1.txt");serumdata #将数据读入向量serumdata。
data_outline(serumdata)结果如下:N Mean Var std_dev Median std_mean CV CSS USS R 1 100 73.696 15.41675 3.926417 73.5 0.3926417 5.327857 1526.258 544636.3 20R1 Skewness Kurtosis1 4.6 0.03854249 0.07051809要点:read.table()用于读表格形式的文件。
上述形式的数据由于第七行缺几个数据,故用read.table()不能读入。
scan()可以直接读纯文本文件。
scan()和matrix()连用还可以将数据存放成矩阵形式。
X<-matrix(scan("3.1.txt",0),ncol=10,byrow=TRUE) #将上述数据放置成10*10的矩阵。
scan()还可以从屏幕上直接输入数据。
Y<-scan()然后按提示输入即可。
结束输入时按回车即可。
Ex3.2>hist(serumdata,freq=FALSE,col="purple",border="red",density=3,angle=60,main=paste("the histogram of serumdata"),xlab="age",ylab="frequency")#直方图。
col是填充颜色。
默认空白。
border是边框的颜色,默认前景色。
density是在图上画条纹阴影,默认不画。
angle是条纹阴影的倾斜角度(逆时针方向),默认45度。
main, xlab, ylab是标题,x和y坐标轴名称。
>lines(density(serumdata),col="blue")#密度估计曲线。
>x<-64:85> lines(x,dnorm(x,mean(serumdata),sd(serumdata)),col="green") #正态分布的概率密度曲线> plot(ecdf(serumdata),verticals=TRUE,do.p=FALSE) #绘制经验分布图> lines(x,pnorm(x,mean(serumdata),sd(serumdata)),col="blue") #正态经验分布> qqnorm(serumdata,col="purple") #绘制QQ图> qqline(serumdata,col="red") #绘制QQ直线Ex3.3> stem(serumdata,scale=1) #作茎叶图。