2022年统计建模与R软件课后答案.docx
精选学习资料 - - - - - - - - - 其次章> x<-c1,2,3;y<-c4,5,6 > e<-c1,1,1 > z<-2*x+y+e;z 1 7 10 13 > z1<-crossprodx,y;z1 ,1 1, 32 > z2<-outerx,y;z2 ,1 ,2 ,31, 4 5 62, 8 10 123, 12 15 181 > A<-matrix1:20,nrow=4;B<-matrix1:20,nrow=4,byrow=T > C<-A+B;C 2> D<-A%*%B;D 3> E<-A*B;E 4> F<-A1:3,1:3 5> G<-B,-3名师归纳总结 - - - - - - -第 1 页,共 43 页精选学习资料 - - - - - - - - - > x<-crep1,5,rep2,3,rep3,4,rep4,2;x> H<-matrixnrow=5,ncol=5 > for i in 1:5 + forj in 1:5 + Hi,j<-1/i+j-11> detH2> solveH3> eigenH> 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'> write.tablestudentdata,file='student.txt' > write.csvstudentdata,file='student.csv'count<-functionn if n<=0 print' 要求输入一个正整数 '名师归纳总结 - - - - - - -第 2 页,共 43 页精选学习资料 - - - - - - - - - else repeat if n%2=0 n<-n/2 else n<-3*n+1 ifn=1break print' 运算胜利 ' 第三章第一将数据录入为 x；利用 data_outline 函数；如下 > data_outlinex> histx,freq=F > linesdensityx,col='red' > y<-minx:maxx > linesy,dnormy,73.668,3.9389,col='blue' > plotecdfx,verticals=T,do.p=F > linesy,pnormy,73.668,3.9389名师归纳总结 - - - - - - -第 3 页,共 43 页精选学习资料 - - - - - - - - - > qqnormx> qqlinex> stemx> boxplotx> fivenumx> shapiro.testx> ks.testx,'pnorm',73.668,3.9389One-sample Kolmogorov-Smirnov testdata: xalternative hypothesis: two-sidedWarning message:In ks.testx, "pnorm", 73.668, 3.9389 :ties should not be present for the Kolmogorov-Smirnov test这里显现警告信息是由于 现重复值ks检验要求样本数据是连续的,不答应出名师归纳总结 - - - - - - -第 4 页,共 43 页精选学习资料 - - - - - - - - - >x1<-c2,4,3,2,4,7,7,2,2,5,4;x2<-c5,6,8,5,10,7,12,12,6,6;x3<-c7,11,6,6, 7,9,5,5,10,6,3,10 > boxplotx1,x2,x3,names=c'x1','x2','x3',vcol=c2,3,4 >windows >plotfactorcrep1,lengthx1,rep2,lengthx2,rep3,lengthx3,cx1, x2,x3> rubber<-data.framex1=c65,70,70,69,66,67,68,72,66,68, +x2=c45,45,48,46,50,46,47,43,47,48,x3=c27.6,30.7,31.8,32.6,31.0,31.3 ,37.0,33.6,33.1,34.2 > plotrubber名师归纳总结 - - - - - - -第 5 页,共 43 页精选学习资料 - - - - - - - - - 详细有相关关系的两个变量的散点图要么是从左下角到右上角正相关,要么是从左上角到右下角负相关 ；从上图可知全部的图中偶 读没有这样的趋势,故均不相关；1> student<-read.csv'3.7.csv' > attachstudent > plot体重身高2> coplot体重身高| 性别3> coplot体重身高| 年龄4> coplot体重身高| 年龄+性别 只列出 4的结果,如以下图名师归纳总结 - - - - - - -第 6 页,共 43 页精选学习资料 - - - - - - - - - > x<-seq-2,3,0.5;y<-seq-1,7,0.5 > f<-functionx,y + x4-2*x2*y+x2 -2*x*y+2*y2+9*x/2 -4*y+4 > z<-outerx,y,f >contourx,y,z,levels=c0,1,2,3,4,5,10,15,20,30,40,50,60,80,100,col='blu e' > windows > perspx,y,z,theta=30,phi=30,expand=0.7,col='red'> cor.test身高,体重 依据得出的结果看是相关的；详细结果不再列出> df<-read.csv'48名求职者得分 .csv' > starsdf 然后依据 G的标准来画出星图> attachdf > df$G1<-SC+LC+SMS+DRV+AMB+GSP+POT/7 > df$G2<-FL+EXP+SUIT/3 > df$G3<-LA+HON+KJ/3 > df$G4<-AA > df$G5<-APP名师归纳总结 - - - - - - -第 7 页,共 43 页精选学习资料 - - - - - - - - - > a<-scaledf,17:21> starsa这里从 17 开头取,是由于在df 中将 ID 也作为了一列使用 P159已经编好的函数 unison,接着上题,直接有 > unisona第四章1先求矩估量；总体的期望为a1 xa1dxa1；因此我们有a1E x ；可解a2a2得 a=2*E.-1/1-E. .因此我们用样本的均值来估量a 即可；在R中实现如下 > x<-c0.1,0.2,0.9,0.8,0.7,0.7 > 2*meanx-1/1-meanx2采纳极大似然估量 第一求出极大似然函数为 . .La;x = .+ 1. . .= .+ 1 . . .=1 .=1 再取对数为名师归纳总结 - - - - - - -第 8 页,共 43 页精选学习资料 - - - - - - - - - .ln La;x = nln a + 1 + aln . .=1 最终求导.lnLa;x=.+ 1 + ln .=1. .好了下面开头用 R编程求解,留意此题中n=6.方法一、使用 unniroot 函数> f<-functiona 6/a+1+sumlogx > unirootf,c0,1方法二、使用 optimize 函数> g<-functiona 6*loga+1+a*sumlogx > optimizeg,c0,1,maximum=T用极大似然估量得出 = n/ . .=1 . . .现用 R求解如下 >x<-crep5,365,rep15,245,rep25,150,rep35,100,rep45,70,rep55,45,rep65,25 > 1000/sumx换句话讲, 就是用该样原来估量泊松分布中的参数,然后求出该分布的均值；我们知道泊松分布中的参数 们只需要用样本均值作矩估量即可,既是均值又是方差；因此我名师归纳总结 - - - - - - -第 9 页,共 43 页精选学习资料 - - - - - - - - - 在 R中实现如下> x<-crep0,17,rep1,20,rep2,10,rep3,2,rep4,1 > meanx 1 1> f<-functionx +obj<-c-13+x1+5-x2*x2 -2*x2, -29+x1+x2+1*x2-14*x2 + sumobj2> nlmf,c0.5,-2在矩估量中, 正态分布总体的均值用样本的均值估量；故在 R中实现如下> x<-c54,67,68,78,70,66,67,70,65,69 > meanx然后用作区间估量,如下> t.testx> t.testx,alternative='less' > t.testx,alternative='greater' 此时我们只需要区间估量的结果, 所以我们只看中的关于置信区间的名师归纳总结 输出即可；同时也给出均值检验的结果,但是默认mu=0第 10 页,共 43 页- - - - - - -精选学习资料 - - - - - - - - - 并不是我们想要的；下面我们来做是否低于 下 > t.testx,alternative='greater',mu=72 One Sample t-test data: x72 的均值假设检验；如alternative hypothesis: true mean is greater than 72 95 percent confidence interval: 63.96295 Inf sample estimates: mean of x 结果说明：我们的备择假设是比72 要大,但是 p 值为,所以我们不接受备择假设, 接受原假设比 72 小；因此这 10 名患者的平均脉搏次 数比正常人要小；我们可以用两种方式来做一做> x<-c140,137,136,140,145,148,140,135,144,141 > y<-c135,118,115,140,128,131,130,115,131,125 > t.testx,y,var.equal=T > t.testx-y 结果不再列出, 但是可以发觉用均值差估量和配对数据估量的结果的名师归纳总结 - - - - - - -第 11 页,共 43 页精选学习资料 - - - - - - - - - 数值有一点小小的差异； 但得出的结论是不影响的 他们的期望差异很大> A<-c0.143,0.142,0.143,0.137 > B<-c0.140,0.142,0.136,0.138,0.140 > t.testA,B> x<-c140,137,136,140,145,148,140,135,144,141 > y<-c135,118,115,140,128,131,130,115,131,125 > var.testx,y > t.testx,y,var.equal=F泊松分布的参数就等于它的均值也等于方差；我们直接用样本均值来估量参数即可,然后作样本均值的置信区间即可；> x<-crep0,7,rep1,10,rep2,12,rep3,8,rep4,3,rep5,2 > meanx> t.testx正态总体均值用样本均值来估量；故如下 > x<-c1067,919,1196,785,1126,936,918,1156,920,948名师归纳总结 - - - - - - -第 12 页,共 43 页精选学习资料 - - - - - - - - - > t.testx,alternative='greater' 留意 greater 才是求区间下限的都比它大的意思嘛第五章这是一个假设检验问题, 即检验油漆作业工人的血小板的均值是否为225.在 R中实现如下 > x<-scan 1: 220 188 162 230 145 160 238 188 247 113 11: 126 245 164 231 256 183 190 158 224 175 21: Read 20 items > t.testx,mu=225考察正态密度函数的概率在R中的运算；第一我们要把该正态分布的均值和方差给估量出来,这个就利用样本即可；然后用 pnorm 函数 来运算大于 1000 的概率；如下 > x<-c1067,919,1196,785,1126,936,918,1156,920,948 > pnorm1000,meanx,sdx名师归纳总结 - - - - - - -第 13 页,共 43 页精选学习资料 - - - - - - - - - 这是检验两个总体是否存在差异的问题；可用符号检验和 wilcoxon 秩检验；两种方法实现如下 > x<-c113,120,138,120,100,118,138,123 > y<-c138,116,125,136,110,132,130,110 > binom.testsumx<y,lengthx p-value = 1 > wilcox.testx,y,exact=F可见无论哪种方法 P值都大于,故接受原假设,他们无差异1采纳 w 检验法 >x<-c-0.7,-5.6,2,2.8,0.7,3.5,4,5.8,7.1,-0.5,2.5,-1.6,1.7,3,0.4,4.5,4.6,2.5,6,-1.4 >y<-c3.7,6.5,5,5.2,0.8,0.2,0.6,3.4,6.6,-1.1,6,3.8,2,1.6,2,2.2,1.2,3.1,1.7,-2 > shapiro.testx > shapiro.testy 采纳 ks检验法> ks.testx,'pnorm',meanx,sdx > ks.testy,'pnorm',meany,sdy 采纳 pearson拟合优度法对 x 进行检验 > A<-tablecutx,br=c-2,0,2,4,6,8名师归纳总结 - - - - - - -第 14 页,共 43 页精选学习资料 - - - - - - - - - > A-2,0 0,2 2,4 4,6 6,8 14 4 6 4 发觉 A 中有频数小于 5,故应当重新调整分组> A<-tablecutx,br=c-2,2,4,8> A-2,2 2,4 4,8 8 6 5然后再运算理论分布> p<-pnormc-2,2,4,8,meanx,sdx> p<-cp2,p3-p2,1-p3最终检验> chisq.testA,p=p 采纳 pearson拟合优度法对 y 进行检验> B<-tablecuty,br=c-2.1,1,2,4,7> B-2.1,1 1,2 2,4 4,7 5 5 5 5 > p<-pnormc1,2,4,meany,sdy > p<-cp1,p2-p1,p3 -p2,1-p3> chisq.testB,p=p以上的全部结果都不再列出, 结论是试验组和对比组都是来自正态分名师归纳总结 - - - - - - -第 15 页,共 43 页精选学习资料 - - - - - - - - - 布；2> t.testx,y,var.equal=F> t.testx,y,var.equal=T > t.testx,y,paired=T 结论是均值无差异3> var.testx,y 结论是方差相同 由以上结果可以看出这两种药的成效并无二致1对新药组应用检验也可用检验> x<-c126,125,136,128,123,138,142,116,110,108,115,140 > y<-c162,172,177,170,175,152,157,159,160,162 > p<-pnormc105,125,145,meanx,sdx > p<-cp2,1-p2 > chisq.testA,p=p 对对比组用检验> ks.testy,'pnorm',meany,sdy 结论是他们都听从正态分布2> var.testx,y 结论是方差相同3> wilcox.testx,y,exact=F名师归纳总结 - - - - - - -第 16 页,共 43 页精选学习资料 - - - - - - - - - 结果是有差异明显是要检验二项分布的p 值是否为实现如下> binom.test57,400,p=0.147 结果是支持也就是检验二项分布中的 p 值是否大于> binom.test178,328,p=0.5,alternative='greater' 结果是不能认为能增加比例就是检验你的样本是否符合那个分布> chisq.testc315,101,108,32,p=c9,3,3,1/16 结果显示符合自由组合规律又是检验一个总体是否符合假定分布；> x<-0:5;y<-c92,68,28,11,1,0> z<-repx,y > A<-tablecutz,br=c-1,0,1,2,5> q<-ppoisc0,1,2,5,meanz > p<-cq1,q2-q1,q3 -q2,1-q3 > chisq.testA,p=p 结论是符合泊松分布名师归纳总结 - - - - - - -第 17 页,共 43 页精选学习资料 - - - - - - - - - > x<-c2.36,3.14,7.52,3.48,2.76,5.43,6.54,7.41 > y<-c4.38,4.25,6.53,3.28,7.21,6.55 > ks.testx,y即列联表的的独立性检验> x<-c358,229,2492,2754 > dimx<-c2,2 > chisq.testx或> fisher.testx 结论是有影响> x<-c45,12,10,46,20,28,28,23,30,11,12,35 > dimx<-c4,3 > chisq.testx 结果是相关> x<-c3,4,6,4 > dimx<-c2,2 > fisher.testx 结果显示工艺对产品质量无影响即检验两种讨论方法是否有差异名师归纳总结 - - - - - - -第 18 页,共 43 页精选学习资料 - - - - - - - - - > x<-c58,2,3,1,42,7,8,9,17 > dimx<-c3,3 > mcnemar.testx,correct=F 结果说明两种检测方法有差异> x<-c13.32,13.06,14.02,11.86,13.58,13.77,13.51,14.42,14.44,15.43 > binom.testsumx>14.6,lengthx,al='l' > wilcox.testx,mu=14.6,al='l',exact=F 结果说明是在中位数之下123> x<-scan21: Read 20 items > y<-scan21: Read 20 items > binom.testsumx<y,lengthx名师归纳总结 - - - - - - -第 19 页,共 43 页精选学习资料 - - - - - - - - - > wilcox.testx,y,paired=T,exact=F > wilcox.testx,y,exact=F4> ks.testx,'pnorm',meanx,sdx > ks.testy,'pnorm',meany,sdy > var.testx,y 由以上检验可知数据符合正态分布且方差相同,故可做 t 检验> t.testx,y 可以发觉他们的均值是有差异的5综上所述, Wilcoxon符号秩检验的差异检出才能最强,符号检验的差异检出最弱；> x<-c24,17,20,41,52,23,46,18,15,29 > y<-c8,1,4,7,9,5,10,3,2,6 > cor.testx,y,method='spearman' > cor.testx,y,method='kendall' 有关系的> x<-1:5 > y<-crepx,c0,1,9,7,3 > z<-crepx,c2,2,11,4,1 > wilcox.testy,z,exact=F 结果显示这两种疗法没什么区分名师归纳总结 - - - - - - -第 20 页,共 43 页精选学习资料 - - - - - - - - - 第六章1> snow<-data.frameX=c5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4, + Y=c1907,1287,2700,2373,3260,3000,1947,2273,3113,2493 > plotsnow$X,snow$Y 结论是有线性关系的；23> lm.sol<-lmY1+X,data=snow;summarylm.sol 结果是方程是显著的4> predictlm.sol,data.frameX=7,interval='prediction',level=0.95fit lwr upr2454.971 12> soil<-data.frameX1=c0.4,0.4,3.1,0.6,4.7,1.7,9.4,10.1,11.6,12.6, + 10.9,23.1,23.1,21.6,23.1,1.9,26.8,29.9,X2=c52,23,19,34,24,65,44,31, + 29,58,37,46,50,44,56,36,58,51,X3=c158,163,37,157,59,123,46,117, + 173,112,111,114,134,73,168,143,202,124,Y=c64,60,71,61,54,77,81, + 93,93,51,76,96,77,93,95,54,168,99 > lm.sol<-lmY1+X1+X2+X3,data=soil;summarylm.sol 我们发觉 X2和 X3的系数没有通过 t 检验；但是整个方程通过了检验；名师归纳总结 - - - - - - -第 21 页,共 43 页精选学习资料 - - - - - - - - - 3> lm.ste<-steplm.sol> summarylm.ste可以发觉新模型只含有X1 和 X3,但是 X3的系数仍是不显著；接下来考虑用 drop1 函数处理 > drop1lm.ste 发觉去掉 X3 残差上升最小, AIC只是有少量增加；因此应当去掉 X3> lm.new<-lmYX1,data=soil;summarylm.new 此时发觉新模型系数显著且方程显著1> da<-data.frameX=c1,1,1,1,2,2,2,3,3,3,4,4,4,5,6,6,6,7,7,7,8,8,8, + 9,11,12,12,12,Y=c0.6,1.6,0.5,1.2,2.0,1.3,2.5,2.2,2.4,1.2,3.5,4.1, + 5.1,5.7,3.4,9.7,8.6,4.0,5.5,10.5,17.5,13.4,4.5,30.4,12.4,13.4, + 26.2,7.4 > plotda$X,da$Y > lm.sol<-lmYX,data=da > ablinelm.sol2> summarylm.sol 全部通过3> plotlm.sol,1 > windows > plotlm.sol,3名师归纳总结 可以观看到误差符合等方差的；但是有残差反常值点24,27,28.第 22 页,共 43 页- - - - - - -精选学习资料 - - - - - - - - - 4> lm.up<-updatelm.sol,sqrt.> summarylm.up 都通过检验 > plotda$X,da$Y > ablinelm.up > windows > plotlm.up,1 > windows > plotlm.up,3 可以发觉仍是有残差离群值 24,28> lm.sol<-lmY1+X1+X2,data=toothpaste;summarylm.sol > influence.measureslm.sol > plotlm.sol,3 通过函数发觉 5,8,9,24对样本影响较大,可能是反常值点,而通过残差图发觉 5 是残差离群点,但是整个残差仍是在 可考虑剔除 5,8,9,24 点再做拟合；-2,2之内的；因此> lm.new<-lmY1+X1+X2,data=toothpaste,subset=c-5,-8,-9,-24 > windows > plotlm.new,3 > summarylm.new我们发觉模型的残差都掌握在-1.5,1.5之内,而且方程系数和方程本名师归纳总结 - - - - - - -第 23 页,共 43 页精选学习资料 - - - - - - - - - 身也都通过检验；> cement<-data.frameX1=c7,1,11,11,7,11,3,1,2,21,1,11,10, + X2=c26,29,56,31,52,55,71,31,54,47,40,66,68, + X3=c6,15,8,8,6,9,17,22,18,4,23,9,8, + X4=c60,52,20,47,33,22,6,44,22,26,34,12,12, +Y=c78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3,1 09.4 > XX<-corcement1:4 > kappaXX,exact=T> eigenXX 发觉变量的多重共线性很强,且有0.241X1+0.641X2+0.268X3+0.676X4=0 说明 X1,X2,X3,X4 多重共线；其实逐步回来可以解决多重共线的问题；我们可以检验一下step 函数去掉变量后的共线性；step 去掉了X3和 X4；我们看看去掉他们的共线性如何；> XX<-corcement1:2 > kappaXX,exact=T我们发觉去掉 X3 和 X4 后,条件数降低好多好多；说明 step 函数是 合理的；名师归纳总结 - - - - - - -第 24 页,共 43 页精选学习资料 - - - - - - - - - 第一得把这个表格看懂；里面的数字应当是有感染和无感染的人数；而影响变量有三个；我们把这些影响变量进行编码；如下；抗生素 X1 发生2 不发生3 危急因子X2 4 5 有无方案X3 6 7 是否感染 Y 1 0 对数据的处理,如下X1 2 X2 4 X3 6 Y 1 频数1 2 4 6 0 17 2 5 6 1 0 2 5 6 0 2 2 4 7 1 11 2 4 7 0 87 2 5 7 1 0 2 5 7 0 0 3 4 6 1 28 3 4 6 0 30 3 4 7 1 23 3 4 7 0 3 3 5 6 1 8 3 5 6 0 32 3 5 7 1 0 3 5 7 0 9 然后用 R处理并求解模型 >hospital<-data.frameX1=repc2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,c1,17,0,2 ,11,87, + 0,0,28,30,23,3,8,32,0,9,X2=repc4,4,5,5,4,4,5,5,4,4,4,4,5,5,5,5, + c1,17,0,2,11,87, + 0,0,28,30,23,3,8,32,0,9,X3=repc6,6,6,6,7,7,7,7,6,6,7,7,6,6,7,7,名师归纳总结 - - - - - - -第 25 页,共 43 页精选学习资料 - - - - - - - - - + c1,17,0,2,11,87,0,0,28,30,23,3,8,32,0,9, + Y=repc1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,c1,17,0,2,11,87,0,0,28,30,23,3,8, 32,0,9 + > glm.sol<-glmYX1+X2+X3,family=binomial,data=hospital > summaryglm.sol 可以发觉假如显著性为,就方程的系数和方程本省全部通过检验；下面我们来做一个猜测,看看使用抗生素,有危急因子,有方案的一个孕妇发生感染的概率是多少；> pre<-predictglm.s
