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(要求輸入一個正整數)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這里出現警告信息是因為ks檢驗要求樣本數據是連續的，不允許出現重復值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)具體有相關關系的兩個變量的散點圖要么是從左下角到右上角（正相關），要么是從左上角到右下角（負相關）。從上圖可知所有的圖中偶讀沒有這樣的趨勢，故均不相關。3.7（1） student attach(student) plot(體重身高)（2） coplot(體重身高|性別)（3） coplot(體重身高|年齡)（

5、4） coplot(體重身高|年齡+性別)只列出（4）的結果，如下圖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(身高,體重)根據得出的結果看是相關的。具體結果不再列出3.10 df stars(df)然后按照G的標準來畫出星圖 attach(df) df$G1 df$G2 df$G3 df$G4 df$G

6、5 a stars(a)這里從17開始取，是因為在df中將ID也作為了一列3.11使用P159已經編好的函數unison，接著上題，直接有 unison(a)第四章4.1（1）先求矩估計。總體的期望為。因此我們有。可解得a=(2*E（x）-1）/(1-E(x).因此我們用樣本的均值來估計a即可。在R中實現如下 x (2*mean(x)-1)/(1-mean(x)1 0.3076923（2）采用極大似然估計首先求出極大似然函數為La;x=i=1na+1xia=(a+1)ni=1nxia再取對數為lnLa;x=nlna+1+aln(i=1nxi最后求導lnL(a;x)a=na+1+lni=1nxi

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

8、樣本的均值估計。故在R中實現如下 x mean(x)1 67.4然后用t.test作區間估計，如下 t.test(x) t.test(x,alternative=less) t.test(x,alternative=greater)此時我們只需要區間估計的結果，所以我們只看t.test中的關于置信區間的輸出即可。t.test同時也給出均值檢驗的結果，但是默認mu=0并不是我們想要的。下面我們來做是否低于72的均值假設檢驗。如下 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結果說明：我們的備擇假設是比72要大，但是p值為0.9817，所以我們不接受備擇假設，接受原假設比72小。因此這10名患者的平均脈搏次數比正常人要小。4.6我們可以用兩種方式來做一做 x y t.test(x,y,var.equal=T) t.test(x-y)結果不再列出，但是可以發現用均值差估計和配對數據估計的結果

10、的數值有一點小小的差別。但得出的結論是不影響的（他們的期望差別很大）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泊松分布的參數就等于它的均值也等于方差。我們直接用樣本均值來估計參數即可，然后作樣本均值0.95的置信區間即可。 x mean(x)1 1.904762 t.test(x)4.10正態總體均值用樣本均值來估計。故如下 x t.test(x,alternative=greater)注意greater才是求區間下限的（都比它大的意思嘛）第五章5.1這是一個假設檢驗問題，即檢驗油漆作業工人的血小板的均

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

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

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

14、)結論是他們都服從正態分布（2） var.test(x,y)結論是方差相同（3） wilcox.test(x,y,exact=F)結果是有差別5.6明顯是要檢驗二項分布的p值是否為0.147.R實現如下 binom.test(57,400,p=0.147)結果是支持5.7也就是檢驗二項分布中的p值是否大于0.5 binom.test(178,328,p=0.5,alternative=greater)結果是不能認為能增加比例5.8就是檢驗你的樣本是否符合那個分布 chisq.test(c(315,101,108,32),p=c(9,3,3,1)/16)結果顯示符合自由組合規律5.9又是檢驗一個

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

16、.6),length(x),al=l) wilcox.test(x,mu=14.6,al=l,exact=F)結果表明是在中位數之下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)由以上檢驗可知數據符合正態分布且方差相同，故可做t檢驗 t.test(x,y)可以發現他們的均值是有差別的（5）綜上所述，Wil

17、coxon符號秩檢驗的差異檢出能力最強，符號檢驗的差異檢出最弱。5.17 x y cor.test(x,y,method=spearman) cor.test(x,y,method=kendall)有關系的5.18 x y z wilcox.test(y,z,exact=F)結果顯示這兩種療法沒什么區別第六章6.1（1） snow plot(snow$X,snow$Y)結論是有線性關系的。（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)可以發現新模型只含有X1和X3，但是X3的系數還是不顯著。接下來考慮用drop1函數處理 drop1(lm.ste)發現去掉X3殘差升高最小，AIC只是有少量增加。因此應該去掉X3 lm.new da plot(da$X,da$Y) lm.sol abline(lm.sol)（2） summary(lm.sol)全部通過（3） plot(lm.sol,1) windows() plot(lm.sol,3)可以觀察到誤差符合等方差的。但是有殘差異常值點24,27,28.（4） lm.u

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

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

21、低好多好多。說明step函數是合理的。6.6首先得把這個表格看懂。里面的數字應該是有感染和無感染的人數。而影響變量有三個。我們把這些影響變量進行編碼。如下。發生不發生抗生素X123危險因子X245有無計劃X367是否感染Y10對數據的處理，如下X1X2X3Y頻數246112460172561025602247111247087257102570034612834603034712334703356183560323571035709然后用R處理并求解模型hospital glm.sol summary(glm.sol)可以發現如果顯著性為0.1，則方程的系數和方程本省全部通過檢驗。下面我們來做

22、一個預測，看看（使用抗生素，有危險因子，有計劃）的一個孕婦發生感染的概率是多少。 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)可以發現各變量影響基本都不顯著，甚至大部分還沒通過顯著性檢驗。只有X1的系數通過了顯著性檢驗，但是也不是很理想。下面計算每一個病人的

23、生存時間大于200天的概率值。pre p p（2） lm.stepre p.new p.new顯然經過逐步回歸后的模型更合理。用summary(lm.ste)看，第二個模型通過了顯著性檢驗（a=0.1）6.9(1) 首先將公式線性化，對方程兩邊直接取對數即可。然后將得到的方程用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 *發現所求的基本上與內在線性相同。第七章7.1（1）pro pro.aov summary(pro.aov)可以看到不同工廠對產品的影響是顯著的（2）首先自己編寫求均值的小程序如下 K for(i in 1:3)+ K1,i K 甲 乙 丙mean 10

25、3 111 86然后再用t.test來做均值的置信區間估計 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)可以發現影響是非常顯著的。（2） pairwise.t.test(old$Y,old$X)直接從結果就可以發現國內只有以工廠和丙工廠與國外工廠有顯著差異。而國內只有甲乙，甲丙之間存在著顯著差異。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)可以看到數據符合正態性但是不是方差齊性的7.4 rat rat.aov summary(rat.aov)結果是顯著的7.5 sleep sleep.aov summary(sleep.aov)結果是不顯著7.6（1） pro pro.aovAB。下面我們來計算它們各個水平下的均值。首先要交互作用給找出來。如下 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應選擇水平2，即A和B 不等的是最好的。然后選擇A，選擇水平3，那么B只能在1和2中選擇，需選擇1.于是我們的最優組合為A3B1。下面給出A3B1的點估計和區間估計。 mean(pro$Ypro$A=3&pro$B=1) t.test(pro$Ypro$A=3&pro$B=1)（

28、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.27500 55.525003 56.12500 54.20000 57.62500所以應該選品種8號，密度4.5，施肥量0.757.8首先我們繪制出正交試驗表格，如下列號1234567產量試驗號ABA*

29、BCA*CB*CDC*DB*DA*D11111111862111222295312211229141222211945212121291621221219672211221838221211288好吧，表示因為多了一個因素D不知道怎么排列交互作用了，我上面排列的也不一定對。此題暫且不做7.9首先把正交試驗表的結果那一列給計算出來。如下 pro pro.mean pro.datapro.aovABAC，其余均不顯著下面再計算出均值，從而就可以依據顯著性來選擇最優參數了 ab-function(x,y)+ n-length(x);z-rep(0,n)+ for( i in 1:n)+ if(xi=

30、yi)zi-1elsezi pro.data$AB pro.data$AC K for(i in 2:6)+ for(j in 1:2)+ Kj,i-1 K A B C AB AC1 1.83750 1.43750 1.85625 1.64375 1.937502 1.80625 2.20625 1.78750 2.00000 1.70625依據顯著性，首先選擇B，選擇B1。再依據AB,應選擇AB1，也就是說A和B應該是同一水平。那么A就被先選定的B決定了它應該選水平1.然后看AC，應該選2.也就是說A和C應該是不同水平。那么A選擇1，C必須選擇2.所以最后的最優組合應該是A1B1C2即通用夾

31、具，特殊鑄鐵，留研量0.015第八章8.1 x g distinguish.distance(x,g,c(8.1,2) distinguish.bayes(x,g,TstX=c(8.1,2) distinguish.bayes(x,g,TstX=c(8.1,2),var.equal=T) discriminiant.fisher(x1:10,x11:20,c(8.1,2)得出的結論都是明天下雨8.2 heart G distinguish.distance(heart,G,var.equal=F) distinguish.distance(heart,G,var.equal=T) distin

32、guish.bayes(heart,G,p=c(11/23,7/23,5/23),var.equal=F) distinguish.bayes(heart,G,p=c(11/23,7/23,5/23),var.equal=T)無論方差相同還是不同，對于距離判別的正確率都是78.2%而方差不同的貝葉斯判別正確率僅僅為65.2%方差相同的貝葉斯判別正確率為87%8.3（1） studyX d hc.1 hc.2 hc.3 hc.4 opar plot(hc.1,hang=-1) rect1 plot(hc.2,hang=-1) rect2 plot(hc.3,hang=-1) rect3 plot(hc.4,hang=-1) rect4 rect1 rect2 rect3 rect4（2） km sort(km$clust)8.4 coreerX d hc1 hc2 hc3 hc4 opar plot(hc1,hang=-1) rect1 plot(hc2,hang=-1) rect2 plot(hc3,hang=-1) rect3 plot(