完整版)統(tǒng)計(jì)建模與R軟件課后答案44頁(yè)

1、第二章2.1 x-c(1,2,3);y e z z1 z2 A-matrix(1:20,nrow=4);B C D E F G x H for (i in 1:5)+ for(j in 1:5)+ Hi,j det(H)（2） solve(H)（3） eigen(H) 2.5 studentdata write.table(studentdata,file=student.txt) write.csv(studentdata,file=student.csv)2.7count-function(n)if (n=0)print(要求輸入一個(gè)正整數(shù))elserepeatif (n%2=0)n-n/

2、2elsen data_outline(x)3.2 hist(x,freq=F) lines(density(x),col=red) y lines(y,dnorm(y,73.668,3.9389),col=blue) plot(ecdf(x),verticals=T,do.p=F) lines(y,pnorm(y,73.668,3.9389) qqnorm(x) qqline(x)3.3 stem(x) boxplot(x) fivenum(x)3.4 shapiro.test(x) ks.test(x,pnorm,73.668,3.9389) One-sample Kolmogorov-S

3、mirnov testdata: xD = 0.073, p-value = 0.6611alternative hypothesis: two-sidedWarning message:In ks.test(x, pnorm, 73.668, 3.9389) : ties should not be present for the Kolmogorov-Smirnov test這里出現(xiàn)警告信息是因?yàn)閗s檢驗(yàn)要求樣本數(shù)據(jù)是連續(xù)的，不允許出現(xiàn)重復(fù)值3.5x1-c(2,4,3,2,4,7,7,2,2,5,4);x2-c(5,6,8,5,10,7,12,12,6,6);x3 boxplot(x1,x

4、2,x3,names=c(x1,x2,x3),vcol=c(2,3,4)windows()plot(factor(c(rep(1,length(x1),rep(2,length(x2),rep(3,length(x3),c(x1,x2,x3)3.6 rubber plot(rubber)具體有相關(guān)關(guān)系的兩個(gè)變量的散點(diǎn)圖要么是從左下角到右上角（正相關(guān)），要么是從左上角到右下角（負(fù)相關(guān)）。從上圖可知所有的圖中偶讀沒(méi)有這樣的趨勢(shì)，故均不相關(guān)。3.7（1） student attach(student) plot(體重身高)（2） coplot(體重身高|性別)（3） coplot(體重身高|年齡)（

5、4） coplot(體重身高|年齡+性別)只列出（4）的結(jié)果，如下圖3.8 x-seq(-2,3,0.5);y f zcontour(x,y,z,levels=c(0,1,2,3,4,5,10,15,20,30,40,50,60,80,100),col=blue) windows() persp(x,y,z,theta=30,phi=30,expand=0.7,col=red)3.9 cor.test(身高,體重)根據(jù)得出的結(jié)果看是相關(guān)的。具體結(jié)果不再列出3.10 df stars(df)然后按照G的標(biāo)準(zhǔn)來(lái)畫出星圖 attach(df) df$G1 df$G2 df$G3 df$G4 df$G

6、5 a stars(a)這里從17開始取，是因?yàn)樵赿f中將ID也作為了一列3.11使用P159已經(jīng)編好的函數(shù)unison，接著上題，直接有 unison(a)第四章4.1（1）先求矩估計(jì)?？傮w的期望為。因此我們有?？山獾胊=(2*E（x）-1）/(1-E(x).因此我們用樣本的均值來(lái)估計(jì)a即可。在R中實(shí)現(xiàn)如下 x (2*mean(x)-1)/(1-mean(x)1 0.3076923（2）采用極大似然估計(jì)首先求出極大似然函數(shù)為L(zhǎng)a;x=i=1na+1xia=(a+1)ni=1nxia再取對(duì)數(shù)為lnLa;x=nlna+1+aln(i=1nxi最后求導(dǎo)lnL(a;x)a=na+1+lni=1nxi

7、好了下面開始用R編程求解，注意此題中n=6.方法一、使用unniroot函數(shù) f uniroot(f,c(0,1)方法二、使用optimize函數(shù) g optimize(g,c(0,1),maximum=T)4.2用極大似然估計(jì)得出=n/i=1nxi.現(xiàn)用R求解如下x 1000/sum(x)4.3換句話講，就是用該樣本來(lái)估計(jì)泊松分布中的參數(shù)，然后求出該分布的均值。我們知道泊松分布中的參數(shù)，既是均值又是方差。因此我們只需要用樣本均值作矩估計(jì)即可在R中實(shí)現(xiàn)如下 x mean(x)1 14.4 f-function(x) +obj nlm(f,c(0.5,-2)4.5在矩估計(jì)中，正態(tài)分布總體的均值用

8、樣本的均值估計(jì)。故在R中實(shí)現(xiàn)如下 x mean(x)1 67.4然后用t.test作區(qū)間估計(jì)，如下 t.test(x) t.test(x,alternative=less) t.test(x,alternative=greater)此時(shí)我們只需要區(qū)間估計(jì)的結(jié)果，所以我們只看t.test中的關(guān)于置信區(qū)間的輸出即可。t.test同時(shí)也給出均值檢驗(yàn)的結(jié)果，但是默認(rèn)mu=0并不是我們想要的。下面我們來(lái)做是否低于72的均值假設(shè)檢驗(yàn)。如下 t.test(x,alternative=greater,mu=72) One Sample t-testdata: xt = -2.4534, df = 9, p-v

9、alue = 0.9817alternative hypothesis: true mean is greater than 7295 percent confidence interval: 63.96295 Infsample estimates:mean of x 67.4結(jié)果說(shuō)明：我們的備擇假設(shè)是比72要大，但是p值為0.9817，所以我們不接受備擇假設(shè)，接受原假設(shè)比72小。因此這10名患者的平均脈搏次數(shù)比正常人要小。4.6我們可以用兩種方式來(lái)做一做 x y t.test(x,y,var.equal=T) t.test(x-y)結(jié)果不再列出，但是可以發(fā)現(xiàn)用均值差估計(jì)和配對(duì)數(shù)據(jù)估計(jì)的結(jié)果

10、的數(shù)值有一點(diǎn)小小的差別。但得出的結(jié)論是不影響的（他們的期望差別很大）4.7 A B t.test(A,B)4.8 x y var.test(x,y) t.test(x,y,var.equal=F)4.9泊松分布的參數(shù)就等于它的均值也等于方差。我們直接用樣本均值來(lái)估計(jì)參數(shù)即可，然后作樣本均值0.95的置信區(qū)間即可。 x mean(x)1 1.904762 t.test(x)4.10正態(tài)總體均值用樣本均值來(lái)估計(jì)。故如下 x t.test(x,alternative=greater)注意greater才是求區(qū)間下限的（都比它大的意思嘛）第五章5.1這是一個(gè)假設(shè)檢驗(yàn)問(wèn)題，即檢驗(yàn)油漆作業(yè)工人的血小板的均

11、值是否為225.在R中實(shí)現(xiàn)如下 x t.test(x,mu=225)5.2考察正態(tài)密度函數(shù)的概率在R中的計(jì)算。首先我們要把該正態(tài)分布的均值和方差給估計(jì)出來(lái)，這個(gè)就利用樣本即可。然后用pnorm函數(shù)來(lái)計(jì)算大于1000的概率。如下 x pnorm(1000,mean(x),sd(x)1 0.5087941 1-0.50879411 0.49120595.3這是檢驗(yàn)兩個(gè)總體是否存在差異的問(wèn)題。可用符號(hào)檢驗(yàn)和wilcoxon秩檢驗(yàn)。兩種方法實(shí)現(xiàn)如下 x y binom.test(sum(x wilcox.test(x,y,exact=F)p-value = 0.792可見(jiàn)無(wú)論哪種方法P值都大于0.05

12、，故接受原假設(shè)，他們無(wú)差異5.4（1）采用w檢驗(yàn)法xy shapiro.test(x) shapiro.test(y)采用ks檢驗(yàn)法 ks.test(x,pnorm,mean(x),sd(x) ks.test(y,pnorm,mean(y),sd(y)采用pearson擬合優(yōu)度法對(duì)x進(jìn)行檢驗(yàn) A A(-2,0 (0,2 (2,4 (4,6 (6,8 4 4 6 4 1發(fā)現(xiàn)A中有頻數(shù)小于5，故應(yīng)該重新調(diào)整分組 A A(-2,2 (2,4 (4,8 8 6 5然后再計(jì)算理論分布 p p chisq.test(A,p=p)采用pearson擬合優(yōu)度法對(duì)y進(jìn)行檢驗(yàn) B B(-2.1,1 (1,2 (2

13、,4 (4,7 5 5 5 5 p p chisq.test(B,p=p)以上的所有結(jié)果都不再列出，結(jié)論是試驗(yàn)組和對(duì)照組都是來(lái)自正態(tài)分布。（2） t.test(x,y,var.equal=F) t.test(x,y,var.equal=T) t.test(x,y,paired=T)結(jié)論是均值無(wú)差異（3） var.test(x,y)結(jié)論是方差相同由以上結(jié)果可以看出這兩種藥的效果并無(wú)二致5.5（1）對(duì)新藥組應(yīng)用chisq.test檢驗(yàn)（也可用ke.test檢驗(yàn)） x y p p chisq.test(A,p=p)對(duì)對(duì)照組用ks.test檢驗(yàn) ks.test(y,pnorm,mean(y),sd(y

14、)結(jié)論是他們都服從正態(tài)分布（2） var.test(x,y)結(jié)論是方差相同（3） wilcox.test(x,y,exact=F)結(jié)果是有差別5.6明顯是要檢驗(yàn)二項(xiàng)分布的p值是否為0.147.R實(shí)現(xiàn)如下 binom.test(57,400,p=0.147)結(jié)果是支持5.7也就是檢驗(yàn)二項(xiàng)分布中的p值是否大于0.5 binom.test(178,328,p=0.5,alternative=greater)結(jié)果是不能認(rèn)為能增加比例5.8就是檢驗(yàn)?zāi)愕臉颖臼欠穹夏莻€(gè)分布 chisq.test(c(315,101,108,32),p=c(9,3,3,1)/16)結(jié)果顯示符合自由組合規(guī)律5.9又是檢驗(yàn)一個(gè)

15、總體是否符合假定分布。 x-0:5;y z A q p chisq.test(A,p=p)結(jié)論是符合泊松分布5.10 x y ks.test(x,y)5.11即列聯(lián)表的的獨(dú)立性檢驗(yàn) x dim(x) chisq.test(x)或 fisher.test(x)結(jié)論是有影響5.12 x dim(x) chisq.test(x)結(jié)果是相關(guān)5.13 x dim(x) fisher.test(x)結(jié)果顯示工藝對(duì)產(chǎn)品質(zhì)量無(wú)影響5.14即檢驗(yàn)兩種研究方法是否有差異 x dim(x) mcnemar.test(x,correct=F)結(jié)果表明兩種檢測(cè)方法有差異5.15 x binom.test(sum(x14

16、.6),length(x),al=l) wilcox.test(x,mu=14.6,al=l,exact=F)結(jié)果表明是在中位數(shù)之下5.16（1）（2）（3） x y binom.test(sum(x wilcox.test(x,y,paired=T,exact=F) wilcox.test(x,y,exact=F)（4） ks.test(x,pnorm,mean(x),sd(x) ks.test(y,pnorm,mean(y),sd(y) var.test(x,y)由以上檢驗(yàn)可知數(shù)據(jù)符合正態(tài)分布且方差相同，故可做t檢驗(yàn) t.test(x,y)可以發(fā)現(xiàn)他們的均值是有差別的（5）綜上所述，Wil

17、coxon符號(hào)秩檢驗(yàn)的差異檢出能力最強(qiáng)，符號(hào)檢驗(yàn)的差異檢出最弱。5.17 x y cor.test(x,y,method=spearman) cor.test(x,y,method=kendall)有關(guān)系的5.18 x y z wilcox.test(y,z,exact=F)結(jié)果顯示這兩種療法沒(méi)什么區(qū)別第六章6.1（1） snow plot(snow$X,snow$Y)結(jié)論是有線性關(guān)系的。（2）（3） lm.sol predict(lm.sol,data.frame(X=7),interval=prediction,level=0.95) fit lwr upr1 2690.227 2454.

18、971 2925.4846.2（1）（2） soil lm.sol lm.ste summary(lm.ste)可以發(fā)現(xiàn)新模型只含有X1和X3，但是X3的系數(shù)還是不顯著。接下來(lái)考慮用drop1函數(shù)處理 drop1(lm.ste)發(fā)現(xiàn)去掉X3殘差升高最小，AIC只是有少量增加。因此應(yīng)該去掉X3 lm.new da plot(da$X,da$Y) lm.sol abline(lm.sol)（2） summary(lm.sol)全部通過(guò)（3） plot(lm.sol,1) windows() plot(lm.sol,3)可以觀察到誤差符合等方差的。但是有殘差異常值點(diǎn)24,27,28.（4） lm.u

19、p summary(lm.up)都通過(guò)檢驗(yàn) plot(da$X,da$Y) abline(lm.up) windows() plot(lm.up,1) windows() plot(lm.up,3)可以發(fā)現(xiàn)還是有殘差離群值24,286.4 lm.sol influence.measures(lm.sol) plot(lm.sol,3)通過(guò)influence.measures函數(shù)發(fā)現(xiàn)5,8,9,24對(duì)樣本影響較大，可能是異常值點(diǎn)，而通過(guò)殘差圖發(fā)現(xiàn)5是殘差離群點(diǎn)，但是整個(gè)殘差還是在-2,2之內(nèi)的。因此可考慮剔除5,8,9,24點(diǎn)再做擬合。 lm.new windows() plot(lm.new,

20、3) summary(lm.new)我們發(fā)現(xiàn)lm.new模型的殘差都控制在-1.5,1.5之內(nèi)，而且方程系數(shù)和方程本身也都通過(guò)檢驗(yàn)。6.5 cement XX kappa(XX,exact=T)1 1376.881 eigen(XX)發(fā)現(xiàn)變量的多重共線性很強(qiáng)，且有0.241X1+0.641X2+0.268X3+0.676X4=0說(shuō)明X1,X2，X3，X4多重共線。其實(shí)逐步回歸可以解決多重共線的問(wèn)題。我們可以檢驗(yàn)一下step函數(shù)去掉變量后的共線性。step去掉了X3和X4。我們看看去掉他們的共線性如何。 XX kappa(XX,exact=T)1 1.59262我們發(fā)現(xiàn)去掉X3和X4后，條件數(shù)降

21、低好多好多。說(shuō)明step函數(shù)是合理的。6.6首先得把這個(gè)表格看懂。里面的數(shù)字應(yīng)該是有感染和無(wú)感染的人數(shù)。而影響變量有三個(gè)。我們把這些影響變量進(jìn)行編碼。如下。發(fā)生不發(fā)生抗生素X123危險(xiǎn)因子X(jué)245有無(wú)計(jì)劃X367是否感染Y10對(duì)數(shù)據(jù)的處理，如下X1X2X3Y頻數(shù)246112460172561025602247111247087257102570034612834603034712334703356183560323571035709然后用R處理并求解模型hospital glm.sol summary(glm.sol)可以發(fā)現(xiàn)如果顯著性為0.1，則方程的系數(shù)和方程本省全部通過(guò)檢驗(yàn)。下面我們來(lái)做

22、一個(gè)預(yù)測(cè)，看看（使用抗生素，有危險(xiǎn)因子，有計(jì)劃）的一個(gè)孕婦發(fā)生感染的概率是多少。 pre p cofe lm.sol summary(lm.sol)（2） lm.s2 summary(lm.s2)（3） plot(cofe$X,cofe$Y) abline(lm.sol) windows() plot(cofe$X,cofe$Y) lines(spline(cofe$X,fitted(lm.s2)6.8（1） pe glm.sol summary(glm.sol)可以發(fā)現(xiàn)各變量影響基本都不顯著，甚至大部分還沒(méi)通過(guò)顯著性檢驗(yàn)。只有X1的系數(shù)通過(guò)了顯著性檢驗(yàn)，但是也不是很理想。下面計(jì)算每一個(gè)病人的

23、生存時(shí)間大于200天的概率值。pre p p（2） lm.stepre p.new p.new顯然經(jīng)過(guò)逐步回歸后的模型更合理。用summary(lm.ste)看，第二個(gè)模型通過(guò)了顯著性檢驗(yàn)（a=0.1）6.9(1) 首先將公式線性化，對(duì)方程兩邊直接取對(duì)數(shù)即可。然后將得到的方程用lm回歸。 peo lm.sol|t|) (Intercept) 4.037159 0.084103 48.00 5.08e-16 *X -0.037974 0.002284 -16.62 3.86e-10 * lm.sum exp(lm.sum$coefficients1,1)1 56.66512所以theta0=56

24、.66512,theta1=-0.0379（2） nls.sol summary(nls.sol)Parameters: Estimate Std. Error t value Pr(|t|) b0 58.606535 1.472160 39.81 5.70e-15 *b1 -0.039586 0.001711 -23.13 6.01e-12 *發(fā)現(xiàn)所求的基本上與內(nèi)在線性相同。第七章7.1（1）pro pro.aov summary(pro.aov)可以看到不同工廠對(duì)產(chǎn)品的影響是顯著的（2）首先自己編寫求均值的小程序如下 K for(i in 1:3)+ K1,i K 甲 乙 丙mean 10

25、3 111 86然后再用t.test來(lái)做均值的置信區(qū)間估計(jì) pro.jia pro.yi pro.bing pairwise.t.test(pro$Y,pro$X) 1 2 2 0.35 - 3 0.13 0.04可以看到顯著性主要有乙工廠和丙工廠造成7.2（1） old old.aov summary(old.aov)可以發(fā)現(xiàn)影響是非常顯著的。（2） pairwise.t.test(old$Y,old$X)直接從結(jié)果就可以發(fā)現(xiàn)國(guó)內(nèi)只有以工廠和丙工廠與國(guó)外工廠有顯著差異。而國(guó)內(nèi)只有甲乙，甲丙之間存在著顯著差異。7.3 rat shapiro.test(rat$XA=1) shapiro.tes

26、t(rat$Xrat$A=2) shapiro.test(rat$Xrat$A=3) bartlett.test(XA,data=rat)可以看到數(shù)據(jù)符合正態(tài)性但是不是方差齊性的7.4 rat rat.aov summary(rat.aov)結(jié)果是顯著的7.5 sleep sleep.aov summary(sleep.aov)結(jié)果是不顯著7.6（1） pro pro.aovAB。下面我們來(lái)計(jì)算它們各個(gè)水平下的均值。首先要交互作用給找出來(lái)。如下 ab-function(x,y)+ n-length(x);z-rep(0,n)+ for( i in 1:n)+ if(xi=yi)zi-1else

27、zi pro$AB K for( i in 2:4)+ for(j in 1:3)+ Kj,i-1 K A B AB1 5.150000 5.783333 4.9333332 4.533333 4.666667 5.2500003 5.750000 4.983333 NaN按照影響力越大（即P值越?。覀兪紫却_定AB應(yīng)選擇水平2，即A和B 不等的是最好的。然后選擇A，選擇水平3，那么B只能在1和2中選擇，需選擇1.于是我們的最優(yōu)組合為A3B1。下面給出A3B1的點(diǎn)估計(jì)和區(qū)間估計(jì)。 mean(pro$Ypro$A=3&pro$B=1) t.test(pro$Ypro$A=3&pro$B=1)（3） pairwise.t.test(pro$Y,pro$AB) pairwise.t.test(pro$Y,pro$B) pairwise.t.test(pro$Y,pro$A)7.7 rice rice.aov K for(i in 1:3)+ for(j in 1:3)+ Ki,j K 品種 密度 施肥量1 59.53333 62.73333 59.058332 56.55000 55.