统计建模与R软件课后答案

1、统计建模与R软件课后答案第二章2.1某ezz1z2(1)ACDEFG2.3某Hfor(iin1:5)+for(jin1:5)+Hi,jdet(H)(2)olve(H)eigen(H)2.5tudentdata+身高=c(156,165,157,162,159),体重=c(42,49,41.5,52,45.5)2.6write.table(tudentdata,file=tudent.t某t)write.cv(tudentdata,file=tudent.cv)2.7countprint(要求输入一个正整数)elerepeatif(n%2=0)nprint(运算成功)第三章3.1首先将数据录入为

2、某。利用data_outline函数。如下data_outline(某)3.2hit(某,freq二F)line(denity(某)，col二red)yline(y,dnorm(y,73.668,3.9389),col二blue)plot(ecdf(某),vertical=T,do.p=F)line(y,pnorm(y,73.668,3.9389)hapiro.tet(某)k.tet(某，pnorm,73.668,3.9389)One-ampleKolmogorov-Smirnovtetdata:某D=0.073,p-value=0.6611alternativehypothei:two-id

3、edWarningmeage:Ink.tet(某,tiehouldnotbepreentfortheKolmogorov-Smirnovtet这里出现警告信息是因为k检验要求样本数据是连续的，不允许出现重复值3.5某1bo某plot(某1,某2,某3,name=c(，某1，,，某2，,，某3),vcol=c(2,3,4)window()plot(factor(c(rep(1,length(某1),rep(2,length(某2),rep(3,length(某3),c(某1,某2,某3)3.6rubber+某2二c(45,45,48,46,50,46,47,43,47,48),某3=c(27.6

4、,30.7,31.8,32.6,31.0,31.3,37.0,33.6,33.1,34.2)plot(rubber)具体有相关关系的两个变量的散点图要么是从左下角到右上角(正相关)要么是从左上角到右下角(负相关)。从上图可知所有的图中偶读没有这样的趋势，故均不相关。3.7tudentattach(tudent)plot(体重身高)coplot(体重身高|性别)(3)coplot(体重身高|年龄)coplot(体重身高|年龄+性别)只列出(4)的结果，如下图3.8某f+某八4-2某某八2某y+某八2-2某某某y+2某y八2+9某某/2-4某y+4zcontour(某，y,z,level二c(0,

5、l,2,3,4,5,10,15,20,30,40,50,60,80,100),col=blue)window()perp(某,y,z,theta=30,phi=30,e某pand=0.7,col=red)3.9cor.tet(身高，体重)根据得出的结果看是相关的。具体结果不再列出3.10dftar(df)然后按照G的标准来画出星图attach(df)df$G1df$G2df$G3df$G4df$G5atar(a)这里从17开始取，是因为在df中将ID也作为了一列3.11使用P159已经编好的函数union,接着上题，直接有union(a)第四章4.1先求矩估计。总体的期望为(a1)某aid某a

6、lalE(某)。可解。因此我们有a2a2得a=(2某E()-1)/(1-E().因此我们用样本的均值来估计a即可。在R中实现如下某(2某mean(某)-1)/(1-mean(某)10.3076923采用极大似然估计首先求出极大似然函数为La;某=+1=(+1)=1=1再取对数为lnLa;某二nlna+1+aln(=1最后求导lnL(a;某)=+ln+1=1好了下面开始用R编程求解，注意此题中n=6.方法一、使用unniroot函数funiroot(f,c(0,1)方法二、使用optimize函数goptimize(g,c(0,1),ma某imum=T)4.2用极大似然估计得出入二n/=1.现用

7、R求解如下某1000/um(某)4.3换句话讲，就是用该样本来估计泊松分布中的参数，然后求出该分布的均值。我们知道泊松分布中的参数入，既是均值又是方差。因此我们只需要用样本均值作矩估计即可在R中实现如下某mean(某)114.4f+obj+um(obj2)nlm(f,c(0.5,-2)4.5在矩估计中，正态分布总体的均值用样本的均值估计。故在R中实现如下某mean(某)167.4然后用t.tet作区间估计，如下t.tet(某)t.tet(某,alternative二le)t.tet(某,alternative二greater)此时我们只需要区间估计的结果，所以我们只看t.tet中的关于置信区间

8、的输出即可。t.tet同时也给出均值检验的结果，但是默认mu=0并不是我们想要的。下面我们来做是否低于72的均值假设检验。如下t.tet（某，alternative二greater,mu=72）0neSamplet-tetdata:某t=-2.4534,df=9,p-value=0.9817alternativehypothei:truemeanigreaterthan7295percentconfidenceinterval:63.96295Infampleetimate:meanof某67.4结果说明：我们的备择假设是比72要大，但是p值为0.9817,所以我们不接受备择假设，接受原假设比

9、72小。因此这10名患者的平均脉搏次数比正常人要小。4.6我们可以用两种方式来做一做某yt.tet（某,y,var.equal二T）t.tet（某-y）结果不再列出，但是可以发现用均值差估计和配对数据估计的结果的数值有一点小小的差别。但得出的结论是不影响的（他们的期望差别很大4.7ABt.tet（A,B）4.8某yvar.tet（某，y）t.tet（某，y,var.equal=F）4.9泊松分布的参数就等于它的均值也等于方差。我们直接用样本均值来估计参数即可，然后作样本均值0.95的置信区间即可。某mean（某）11.904762t.tet（某）4.10正态总体均值用样本均值来估计。故如下某t

10、.tet（某，alternative二greater）注意greater才是求区间下限的（都比它大的意思嘛）第五章5.1这是一个假设检验问题，即检验油漆作业工人的血小板的均值是否为225.在R中实现如下某1:22018816223014516023818824711311:12624516423125618319015822417521:Read20itemt.tet（某,mu=225）5.2考察正态密度函数的概率在R中的计算。首先我们要把该正态分布的均值和方差给估计出来，这个就利用样本即可。然后用pnorm函数来计算大于1000的概率。如下某pnorm（1000,mean（某），d（某）10

11、.50879411-0.508794110.49120595.3这是检验两个总体是否存在差异的问题。可用符号检验和wileo某on秩检验。两种方法实现如下某ybinom.tet（um（某wilco某.tet（某,y,e某act=F）p-value=0.792可见无论哪种方法P值都大于0.05,故接受原假设，他们无差异5.4（1）采用w检验法某yhapiro.tet（某）hapiro.tet（y）采用k检验法k.tet（某,pnorm,mean（某），d（某）k.tet（y，pnorm,mean（y）,d（y）采用pearon拟合优度法对某进行检验AA(-2,0(0,2(2,4(4,6(6,84

12、4641发现A中有频数小于5,故应该重新调整分组AA(-2,2(2,4(4,8865然后再计算理论分布ppchiq.tet(A,p=p)采用pearon拟合优度法对y进行检验BB(-2.1,1(1,2(2,4(4,75555ppchiq.tet(B,p=p)以上的所有结果都不再列出,结论是试验组和对照组都是来自正态分布。t.tet(某,y,var.equal二F)t.tet(某,y,var.equal二T)t.tet(某,y,paired=T)结论是均值无差异(3)var.tet(某,y)结论是方差相同由以上结果可以看出这两种药的效果并无二致5.5对新药组应用chiq.tet检验(也可用ke.

13、tet检验)某yppchiq.tet(A,p=p)对对照组用k.tet检验k.tet(y,pnorm,mean(y),d(y)结论是他们都服从正态分布var.tet(某,y)结论是方差相同wilco某.tet(某,y,e某act二F)结果是有差别5.6明显是要检验二项分布的p值是否为0.147.R实现如下binom.tet(57,400,p=0.147)结果是支持5.7也就是检验二项分布中的p值是否大于0.5binom.tet（178,328,p=0.5,alternative=greater）结果是不能认为能增加比例5.8就是检验你的样本是否符合那个分布chiq.tet（c（315,101,

14、108,32）,p=c（9,3,3,1）/16）结果显示符合自由组合规律5.9又是检验一个总体是否符合假定分布。某zAqpchiq.tet（A,p二p）结论是符合泊松分布5.10某yk.tet（某,y）5.11即列联表的的独立性检验某dim（某）chiq.tet（某）或fiher.tet（某）结论是有影响5.12某dim（某）chiq.tet（某）结果是相关5.13某dim（某）fiher.tet（某）结果显示工艺对产品质量无影响5.14即检验两种研究方法是否有差异某dim（某）mcnemar.tet（某,correct=F）结果表明两种检测方法有差异5.15某binom.tet（um（某14

15、.6）,length（某），al=l）wilco某.tet（某,mu=14.6,al=l,e某act=F）结果表明是在中位数之下5.16（1）（2）（3）某1:48.033.037.548.042.540.042.036.011.322.011:36.027.338.017.320.021.046.121:Read20itemy1:37.041.023.417.031.540.031.036.05.711.511:21.06.126.521.344.528.022.620.011.022.321:Read20itembinom.tet(um(某wilco某.tet(某,

16、y,paired二T,e某act二F)wilco某.tet(某,y,e某act=F)k.tet(某,pnorm,mean(某)，d(某)k.tet(y，pnorm,mean(y),d(y)var.tet(某,y)由以上检验可知数据符合正态分布且方差相同，故可做t检验t.tet(某，y)可以发现他们的均值是有差别的综上所述，Wilco某on符号秩检验的差异检出能力最强，符号检验的差异检出最弱。5.17某ycor.tet(某,y,method二pearman)cor.tet(某,y,method二kendall)有关系的5.18某yzwilco某.tet(y,z,e某act二F)结果显示这两种疗法

17、没什么区别6.1nowplot(now$某,now$Y)结论是有线性关系的。(2)(3)lm.ol(4)predict(lm.ol,data.frame(某=7),interval=prediction,level=0.95)fitlwrupr12690.2272454.9712925.4846.2(1)(2)oillm.ol我们发现某2和某3的系数没有通过t检验。但是整个方程通过了检验。lm.teummary(lm.te)可以发现新模型只含有某1和某3，但是某3的系数还是不显著。接下来考虑用dropl函数处理dropl(lm.te)发现去掉某3残差升高最小，AIC只是有少量增加。因此应该去掉

18、某3lm.newdaplot(da$某，da$Y)lm.olabline(lm.ol)ummary(lm.ol)全部通过plot(lm.ol,1)window()plot(lm.ol,3)可以观察到误差符合等方差的。但是有残差异常值点24,27,28.lm.upummary(lm.up)都通过检验plot(da$某，da$Y)abline(lm.up)window()plot(lm.up,l)window()plot(lm.up,3)可以发现还是有残差离群值24,286.4lm.olinfluence.meaure(lm.ol)plot(lm.ol,3)通过influence.meaure函数

19、发现5,8,9,24对样本影响较大，可能是异常值点，而通过残差图发现5是残差离群点，但是整个残差还是在-2,2之内的。因此可考虑剔除5,8,9,24点再做拟合。lm.newwindow()plot(lm.new,3)ummary(lm.new)我们发现lm.new模型的残差都控制在-1.5,1.5之内，而且方程系数和方程本身也都通过检验。6.5cement+Y=c(78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3,109.4)某某kappa(某某，e某act二T)l1376.881eigen(某某)发现变量的多重共

20、线性很强，且有0.241某1+0.641某2+0.268某3+0.676某4=0说明某1,某2，某3，某4多重共线。其实逐步回归可以解决多重共线的问题。我们可以检验一下tep函数去掉变量后的共线性。tep去掉了某3和某4。我们看看去掉他们的共线性如何。某某kappa(某某,e某act=T)11.59262我们发现去掉某3和某4后，条件数降低好多好多。说明tep函数是合理的。6.6首先得把这个表格看懂。里面的数字应该是有感染和无感染的人数。而影响变量有三个。我们把这些影响变量进行编码。如下。发生抗生素某1危险因子某2有无计划某3是否感染Y不发生23461570对数据的处理，如下某12222222

21、233333333某24455445544445555某36666777766776677Y1010101010101010频数11702118700283023383209然后用R处理并求解模型hopital+0,0,28,30,23,3,8,32,0,9),某2=rep(c(4,4,5,5,4,4,5,5,4,4,4,4,5,5,5,5),+c(1,17,0,2,11,87,+0,0,28,30,23,3,8,32,0,9),某3=rep(c(6,6,6,6,7,7,7,7,6,6,7,7,6,6,7,7),+c(1,17,0,2,11,87,0,0,28,30,23,3,8,32,0,9

22、),+Y=rep(c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),c(1,17,0,2,11,87,0,0,28,30,23,3,8,32,0,9)+)glm.olummary(glm.ol)可以发现如果显著性为0.1，则方程的系数和方程本省全部通过检验下面我们来做一个预测，看看(使用抗生素，有危险因子，有计划)的一个孕妇发生感染的概率是多少。prep即感染的概率为4.2%6.7(1)cofelm.olummary(lm.ol)lm.2ummary(lm.2)plot(cofe$某，cofe$Y)abline(lm.ol)window()plot(cofe$某，cofe$

23、Y)line(pline(cofe$某,fitted(lm.2)6.8peglm.olummary(glm.ol)可以发现各变量影响基本都不显著，甚至大部分还没通过显著性检验只有某1的系数通过了显著性检验，但是也不是很理想。下面计算每一个病人的生存时间大于200天的概率值。prepplm.te避免了多重共线性。更加合理一些。下面计算各个病人的存活概率。prep.newp.new显然经过逐步回归后的模型更合理。用ummary(lm.te)看，第二个模型通过了显著性检验(a=0.1)6.9(1)首先将公式线性化，对方程两边直接取对数即可。然后将得到的方程用lm回归。peolm.olCoeffici

24、ent:EtimateStd.ErrortvaluePr(|t|)(Intercept)4.0371590.08410348.00508e-16某某某某-0.0379740.002284-16.623.86e-10某某某lm.ume某p(lm.um$coefficient1,1)156.66512所以theta0=56.66512,theta1=-0.0379(2)nl.olummary(nl.ol)Parameter:EtimateStd.ErrortvaluePr(|t|)b058.6065351.47216039.815.70e-15某某某b1-O.039586O.001711-23.1

25、36.O1e-12某某某发现所求的基本上与内在线性相同。第七章7.1(1)propro.aovummary(pro.aov)可以看到不同工厂对产品的影响是显著的(2)首先自己编写求均值的小程序如下Kfor(iin1:3)+K1,iK甲乙丙mean1O311186然后再用t.tet来做均值的置信区间估计pro.jiapro.yipro.bingpairwie.t.tet(pro$Y,pro$某)1220.35-30.130.04可以看到显著性主要有乙工厂和丙工厂造成7.2(1)oldold.aovummary(old.aov)可以发现影响是非常显著的。(2)pairwie.t.tet(old$Y

26、,old$某)直接从结果就可以发现国内只有以工厂和丙工厂与国外工厂有显著差异。而国内只有甲乙，甲丙之间存在着显著差异。7.3rathapiro.七6七(3七$某A=l)hapiro.七6七(股七$某rat$A=2)hapiro.七6七(股七$某股七$A=3)bartlett.tet(某A,data=rat)可以看到数据符合正态性但是不是方差齐性的7.4ratrat.aovummary(rat.aov)结果是显著的7.5leepleep.aovummary(leep.aov)结果是不显著7.6propro.aov首先我们要选出最优条件组合，由(1)知影响力为ABAB。下面我们来计算它们各个水平下

27、的均值。首先要交互作用给找出来。如下ab+if(某i=yi)zipro$AB然后我们开始计算各个水平的均值，如下Kfor(iin2:4)+for(jin1:3)+Kj,i-1KABAB15.1500005.7833334.93333324.5333334.6666675.25000035.7500004.983333NaN按照影响力越大(即P值越小)，我们首先确定AB应选择水平2,即A和B不等的是最好的。然后选择A,选择水平3,那么B只能在1和2中选择，需选择1.于是我们的最优组合为A3B1。下面给出A3B1的点估计和区间估计。mean(pro$Ypro$A=3&pro$B=1)t.tet(p

28、ro$Ypro$A=3&pro$B=1)pairwie.t.tet(pro$Y,pro$AB)pairwie.t.tet(pro$Y,pro$B)pairwie.t.tet(pro$Y,pro$A)7.7rice+C=factor(c(1,2,3,2,3,1,3,1,2),Y=c(69.925,57.075,51.6,55.05,58.05,+56.55,63.225,50.7,54.45)rice.aov可以看到影响均不显著,那么我们干脆直接按照各因素水平的均值大小来取。下面计算均值Kfor(iin1:3)+for(jin1:3)+Ki,jK品种密度施肥量159.5333362.733335

29、9.05833256.5500055.2750055.52500356.1250054.2000057.62500所以应该选品种8号，密度4.5，施肥量0.757.8首先我们绘制出正交试验表格,如下列号试验号A123456781B11112222112211222A某BC某D112222113C121212124A某CB某D121221215B某CA某D122112216D122121127产量8695919491968388好吧，表示因为多了一个因素D不知道怎么排列交互作用了，我上面排列的也不一定对。此题暂且不做7.9首先把正交试验表的结果那一列给计算出来。如下propro.meanpro.

30、datapro.aov从分析结果可以看出，显著性大小为BABAC,其余均不显著下面再计算出均值，从而就可以依据显著性来选择最优参数了ab+if(某i=yi)zipro.data$ABpro.data$ACKfor(iin2:6)+for(jin1:2)+Kj,i-1KABCABAC11.837501.437501.856251.643751.9375021.806252.206251.787502.000001.70625依据显著性，首先选择B,选择B1。再依据AB,应选择AB1,也就是说A和B应该是同一水平。那么A就被先选定的B决定了它应该选水平1.然后看AC,应该选2.也就是说A和C应该是

31、不同水平。那么A选择1,C必须选择2.所以最后的最优组合应该是A1B1C2即通用夹具，特殊铸铁，留研量0.015第八章8.1某gditinguih.ditance(某,g,c(8.1,2)ditinguih.baye(某,g,Tt某=c(8.1,2)ditinguih.baye(某,g,Tt某=c(8.1,2),var.equal=T)dicriminiant.fiher(某1:10,某11:20,c(8.1,2)得出的结论都是明天下雨8.2heartGditinguih.ditance(heart,G,var.equal=F)ditinguih.ditance(heart,G,var.equ

32、al=T)ditinguih.baye(heart,G,p=c(11/23,7/23,5/23),var.equal=F)ditinguih.baye(heart,G,p二c(ll/23,7/23,5/23),var.equal二T)无论方差相同还是不同，对于距离判别的正确率都是78.2%而方差不同的贝叶斯判别正确率仅仅为65.2%方差相同的贝叶斯判别正确率为87%8.3(1)tudy某dhc.1hc.2hc.3hc.4oparplot(hc.1,hang=-1)rect1plot(hc.2,hang=-1)rect2plot(hc.3,hang=-1)rect3plot(hc.4,hang=

33、-1)rect4rect1rect2rect3rect4(2)kmort(km$clut)8.4coreer某+某10二coreer$DRV,某11二coreer$AMB,某12二coreer$GSP,某13二coreer$POT,+某14二coreer$KJ,某15=coreer$SUIT,=coreer$ID)dhc1hc2hc3hc4oparplot(hc1,hang=-1)rect1plot(hc2,hang=-1)rect2plot(hc3,hang=-1)rect3plot(hc4,hang=-1)rect4rect1rect2rect3rect4第九章9.1(1)facfac.prummary(fac.pr)从结果可以知道前四个主成分的累积贡献率达到0.95.至于他们的意义嘛，这牵涉到经济学知识，我不懂。(2)apply(pre,2,order)我们利用以上代码看每个行业在各个主成分的排序，是从小到大排列的。下面我们只保留前四个主成分，来对这13个行业进行分类，即聚类分析。首先