第二章 計(jì)算機(jī)專業(yè)測試

弟一早

2.1

>x<-c(l,2,3)；y<-c(4,5,6)

>e<-c(l,l,l)

>z<-2*x+y+e;z

[1]71013

>zl<-crossprod(x,y);zl

Li]

[1J32

>z2<-outer(x,y);z2

Li]L2]13]

[1,]456

[2,]81012

[3,]121518

2.2

(1)>A<-matrix(l:20,nrow=4);B<-matrix(l:20,nrow=4,byrow=T)

>C<-A+B;C

⑵〉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(l,5),rep(2z3),rep(3,4),rep(4,2));x

2.4

>H<-matrix(nrow=5/ncol=5)

>for(iin1:5)

+for(jin1:5)

+H[iJ]<-l/(i+j-l)

(1)>det(H)

(2)>solve(H)

(3)>eigen(H)

2.5

>studentdata<?data.frame(姓名=c('張三李四，王五，趙六，丁一')

+,性別二女,「男》女二男)女)年齡二31%15匚,16,「14?15,),

+身高:《156曾65*577162,,,159)體重=4427497415「52,/455))

2.6

>write.table(studentdata,file='student.txt')

>write.csv(studentdata,file='student.csv')

2.7

count<-function(n)

if(n<=0)

prints要求輸入一個正整數(shù))

else(

repeat{

if(n%%2==0)

n<-n/2

else

n<-(3*n+l)

if(n==l)break

)

print(運(yùn)算成功，)}

)

弟二早

3.1

首先將數(shù)據(jù)錄入為x。利用data_outline函數(shù)。如下

>data_outline(x)

3.2

>hist(x,freq=F)

>lines(density(x),col='red,)

>y<-min(x):max(x)

,

>lines(y,dnorm(yz73.668,3.9389)zcol=blue')

>plot(ecdf(x),verticals=7；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「pnorm73.668,3.9389)

One-sampleKolmogorov-Smirnovtest

data:x

D=0.073,p-value=0.6611

alternativehypothesis:two-sided

Warningmessage:

Inks.test(x,"pnorm**,73.668,3.9389):

tiesshouldnotbepresentfortheKolmogorov-Smirnovtest

這里出現(xiàn)警告信息是因?yàn)閗s檢驗(yàn)要求樣本數(shù)據(jù)是連續(xù)的，不允許出

現(xiàn)重復(fù)值

3.5

>X1<-C(2,4,3,2A7,7,2/2,5,4);X2<-C(5,6,8,5,101.7,12,12,6,6);X3<-C(7/11,6/6,

7,9,5,5,10,6,3,10)

,

>boxplot(xl,x2,x3,names=c(xl7x27x3'),vcol=c(2,3,4))

>windows()

>plot(factor(c(rep(lJength(xl))/rep(2Jength(x2))/rep(3Jength(x3))))/c(xl/

x2,x3))

3.6

>rubber<-data.frame(xl=c(65/70/70,69766,67,68,72,66,68),

+x2=c(45,45,48,46,50,46,47,43,47,48),x3=c(27.6,30.7,31.8,32.6,31.0,313

,37.0,33.6,33.1,34.2))

>plot(rubber)

具體有相關(guān)關(guān)系的兩個變量的散點(diǎn)圖要么是從左下角到右上角（正相

關(guān)），要么是從左上角到右下角（負(fù)相關(guān)）。從上圖可知所有的圖中偶

讀沒有這樣的趨勢，故均不相關(guān)。

3.7

（1）>student<-read.csv（'3.7.csv'）

>attach（student）

>plot（體重~身高）

（2）>coplot（體重~身高|性別）

（3）>coplot（體重~身高|年齡）

（4）>coplot（體重~身高|年齡+性別）

只列出（4）的結(jié)果，如下圖

Given:

111213141516

身高

3.8

>x<-seq(-2,3,0.5);y<-seq(-l,7,0.5)

>f<-function(x,y)

+xA4-2*xA2*y+xA2-2*x*y+2*yA2+9*x/2-4*y+4

>z<-outer(x,y,f)

>contour(x,y,z,levels=c(0,l,2,3,4,5,10,15,20,30,40,50,60,80,100),col='blu

e')

>windows()

>persp(x/y/z,theta=30,phi=30,expand=0.7,col='red')

3.9

>contest(身高，體重)

根據(jù)得出的結(jié)果看是相關(guān)的。具體結(jié)果不再列出

3.10

>df<-read.csv(,48名求職者得分.csv')

>stars(df)

然后按照G的標(biāo)準(zhǔn)來畫出星圖

>attach(df)

>df$Gl<-(SC+LC+SMS+DRV+AMB+GSP+P0T)/7

>df$G2<-(FL+EXP+SUIT)/3

>df$G3<-(LA+HON+KJ)/3

>df$G4<-AA

>df$G5<-APP

>a<-scale(df[,17:21])

>stars(a)

這里從17開始取，是因?yàn)樵赿f中將ID也作為了一列

3.11

使用P159已經(jīng)編好的函數(shù)unison,接著上題，直接有

>unison(a)

第四章

4.1

(1)先求矩估計(jì)。

+O0|1_

總體的期望為J(a-1)y啰=怒。因此我們有亍言=E(x)o可解

-CO

得a=(2*E(x)-1)/(1-E(幻).因此我們用樣本的均值來估計(jì)a即可。在

R中實(shí)現(xiàn)如下

>x<-c(0.1,0.2,0.9,0.8,0.7,0.7)

>(2*mean(x)-l)/(l-mean(x))

[1]0.3076923

(2)采用極大似然估計(jì)

首先求出極大似然函數(shù)為

nn

L(a;x)=J-1(a+l)xf=(a+l)nxf

1=1i=l

再取對數(shù)為

InL(a;x)=nln(a+1)+aln(1]xt

i=l

最后求導(dǎo)

n

31nL(a;x)n+ln口陽

daQ+1

i=l

好了下血開始用R編程求解，注意此題中n=6.

方法一、

使用unniroot函數(shù)

>f<-function(a)6/(a+l)+sum(log(x))

>uniroot(f,c(0,l))

方法二、

使用optimize函數(shù)

>g<-function(a)6*log(a+l)+a*sum(log(x))

>optimize(g,c(0,l)/maximum=T)

4.2

用極大似然估計(jì)得出入=n/IXi孫現(xiàn)用R求解如下

>x<-c(rep(5/365),rep(15,245),rep(25,150),rep(35/100),rep(45,70),rep(55,

45),rep(65,25))

>1000/sum(x)

4.3

換句話講，就是用該樣本來估計(jì)泊松分布中的參數(shù)，然后求出該分布

的均值。我們知道泊松分布中的參數(shù)入，既是均值又是方差。因此我

們只需要用樣本均值作矩估計(jì)即可

在R中實(shí)現(xiàn)如下

>x<-c(rep(0,17),rep(l,20),rep(2,10),rep(3,2)jep(4,l))

>mean(x)

[1]1

4.4

>f<-function(x){

+obj<-c(-13+x[l]+((5-x[2])*x[2]-2)*x[2],(-29+x[l]+((x[2]+l)*x[2]-14)*x[2])

)

+sum(objA2)}

>nlm(f,c(0.5,-2))

4.5

在矩估計(jì)中，正態(tài)分布總體的均值用樣本的均值估計(jì)。故在R中實(shí)現(xiàn)

如下

>x<-c(54,67,68,78,70,66,67,70,65,69)

>mean(x)

[1]67.4

然后用t.test作區(qū)間估計(jì)，如下

>t.test(x)

>t.test(x,alternative='less')

>t.test(x,alternative='greater')

此時我們只需要區(qū)間估計(jì)的結(jié)果，所以我們只看"est中的關(guān)于置信

區(qū)間的輸出即可。t.test同時也給出均值檢驗(yàn)的結(jié)果，但是默認(rèn)mu=0

并不是我們想要的。下面我們來做是否低于72的均值假設(shè)檢驗(yàn)。如

下

>t.test(x,alternative:'greater',mu=72)

OneSamplet-test

data:x

t=-2.4534,df=9,p-value=0.9817

alternativehypothesis:truemeanisgreaterthan72

95percentconfidenceinterval:

63.96295Inf

sampleestimates:

meanofx

67.4

結(jié)果說明：我們的備擇假設(shè)是比72要大，但是p值為0.9817,所以

我們不接受備擇假設(shè)，接受原假設(shè)比72小。因此這10名患者的平均

脈搏次數(shù)比正常人要小。

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,var.equal=T)

>t.test(x-y)

結(jié)果不再列出，但是可以發(fā)現(xiàn)用均值差估計(jì)和配對數(shù)據(jù)估計(jì)的結(jié)果的

數(shù)值有一點(diǎn)小小的差別。但得出的結(jié)論是不影響的(他們的期望差別

很大)

4.7

>A<-c(0.143,0.142,0.143,0.137)

>B<-c(0.140,0.142,0.136,0.138,0.140)

>t.test(A,B)

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)

>t.test(x,y,var.equal二F)

4.9

泊松分布的參數(shù)就等于它的均值也等于方差。我們直接用樣本均值來

估計(jì)參數(shù)即可，然后作樣本均值0.95的置信區(qū)間即可。

>x<-c(rep(0/7),rep(lz10),rep(2z12),rep(3/8),rep(4z3)jep(5,2))

>mean(x)

[1]1.904762

>t.test(x)

4.10

正態(tài)總體均值用樣本均值來估計(jì)。故如下

>x<-c(1067,919,1196,785,1126,936,918,H56,920,948)

>t.test(x/alternative='greater')

注意greater才是求區(qū)間下限的(都比它大的意思嘛)

第五章

5.1

這是一個假設(shè)檢驗(yàn)問題，即檢驗(yàn)油漆作業(yè)工人的血小板的均值是否為

225.在R中實(shí)現(xiàn)如下

>x<-scan()

1:220188162230145160238188247113

11:126245164231256183190158224175

21:

Read20items

>t.test(x,mu=225)

5.2

考察正態(tài)密度函數(shù)的概率在R中的計(jì)算。首先我們要把該正態(tài)分布的

均值和方差給估計(jì)出來，這個就利用樣本即可。然后用pnorm函數(shù)

來計(jì)算大于1000的概率。如下

>x<-c(1067z919,1196,785,1126,936,918,1156,920,948)

>pnorm(1000zmean(x),sd(x))

[1]0,5087941

>1-0.5087941

[1]0,4912059

5.3

這是檢驗(yàn)兩個總體是否存在差異的問題。可用符號檢驗(yàn)和Wilcoxon

秩檢驗(yàn)。兩種方法實(shí)現(xiàn)如下

>x<-c(113,120,138,120,100,118,138,123)

>y<-c(138,116,125,136,110,132,130,HO)

>binom.test(sum(x<y),length(x))

p-value=1

>wilcox.test(x,y,exact=F)

p-value=0.792

可見無論哪種方法P值都大于0.05,故接受原假設(shè)，他們無差異

5.4

(1)采用w檢驗(yàn)法

>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<-c(3.7,6.5,5,520.8,0.2,0.6,3.4,6.6廠1.1,638,2,1.6,2,2.2,L2,3.1」.7廠2)

>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)度法對x進(jìn)行檢驗(yàn)

>A<-table(cut(x/br=c(-2/0/2/4,6,8)))

>A

(-2,0](0,2](2,4](4,6](6,8]

44641

發(fā)現(xiàn)A中有頻數(shù)小于5,故應(yīng)該重新調(diào)整分組

>A<-table(cut(xzbr=c(-2,2,4,8)))

>A

(-2,2](2,4](4,8]

865

然后再計(jì)算理論分布

>p<-pnorm(c(-2z2,4,8),mean(x),sd(x))

>p<-c(p[2],p[3]-p[2],l-p[3])

最后檢驗(yàn)

>chisq.test(A,p=p)

采用pearson擬合優(yōu)度法對y進(jìn)行檢驗(yàn)

>B<-table(cut(y,br=c(-2.1,l,2,4,7)))

>B

(-2.1,1](1,2](2,4](4,刀

5555

>p<-pnorm(c(l,2,4),mean(y),sd(y))

>p<-c(p[l],p[2]-p[l],p[3]-p[2],l-p[3])

>chisq.test(B,p=p)

以上的所有結(jié)果都不再列出，結(jié)論是試驗(yàn)組和對照組都是來自正態(tài)分

布。

(2)>t.test(x,y,var.equal=F)

>t.test(x,yzvar.equal=T)

>t.test(x,y,paired=T)

結(jié)論是均值無差異

(3)>var.test(x,y)

結(jié)論是方差相同

由以上結(jié)果可以看出這兩種藥的效果并無二致

5.5

(1)對新藥組應(yīng)用chisq.test檢驗(yàn)(也可用ke.test檢驗(yàn))

>x<-c(126,125,136,128,123,138,142,116,110,108,115,140)

>y<-c(162,172,177,170,175,152,157,159,160,162)

>p<-pnorm(c(105,125,145),mean(x),sd(x))

>p<-c(p[2],l-p[2])

>chisq.test(A,p=p)

對對照組用ks.test檢驗(yàn)

>ks.test(y,'pnorm',mean(y)/sd(y))

結(jié)論是他們都服從正態(tài)分布

(2)>var.test(xzy)

結(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)愕臉颖臼欠穹夏莻€分布

>chisq.test(c(315,101,108,32),p=c(9z3,3,l)/16)

結(jié)果顯示符合自由組合規(guī)律

5.9

又是檢驗(yàn)一個總體是否符合假定分布。

>x<-0:5;y<-c(92,68,28,11,1,0)

>z<-rep(x,y)

>A<-table(cut(z,br=c(-l,0,1,2,5)))

>q<-ppois(c(0,l,2,5),mean(z))

>p<-c(q[l],q[2]-q[l],q[3]-q[2],l-q[3])

>chisq.test(A,p=p)

結(jié)論是符合泊松分布

5.10

>x<-c(2.36,3.14,7.52,3.48,2.76,5.43,6.54,7.41)

>y<-c(4.38,4.25,6.53,3.28,7.21,6.55)

>ks.test(x,y)

5.11

即列聯(lián)表的的獨(dú)立性檢驗(yàn)

>x<-c(358,229,2492,2754)

>dim(x)<-c(2,2)

>chisq.test(x)或〉fisher.test(x)

結(jié)論是有影響

5.12

>x<-c(45,12,10,46,20,28,28,23,30,11,12,35)

>dim(x)<-c(4,3)

>chisq.test(x)

結(jié)果是相關(guān)

5.13

>X<-C(3,4,6,4)

>dim(x)<-c(2,2)

>fisher.test(x)

結(jié)果顯示工藝對產(chǎn)品質(zhì)量無影響

5.14

即檢驗(yàn)兩種研究方法是否有差異

>x<-c(58/2/34,42,7,8,9,17)

>dim(x)<-c(3,3)

>mcnemar.test(x,correct=F)

結(jié)果表明兩種檢測方法有差異

5.15

>x<-c(13.32,13.06,14.02,11.86,13.58,13.77,13.51,14.42,14.44,15.43)

>binom.test(sum(x>14.6),length(x)/al=T)

>wilcox.test(x,mu=14.6,al=T,exact=F)

結(jié)果表明是在中位數(shù)之下

5.16

(1)(2)(3)

>x<-scan()

1:48.033.037.548.042.540.042.036.011.322.0

11:36,027.314.232.152.038.017.320.021.046.1

21:

Read20items

>y<-scan()

1:37.041.023.417.031.540.031.036.05.711.5

11:21.06.126.521.344.528.022.620.011.022.3

21:

Read20items

>binom.test(sum(x<y),length(x))

>wilcox.test(x,y,paired=1exact=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(xzy)

可以發(fā)現(xiàn)他們的均值是有差別的

(5)綜上所述，Wilcoxon符號秩檢驗(yàn)的差異檢出能力最強(qiáng)，符號檢驗(yàn)的差異檢出最弱。

5.17

>x<-c(24,17,20,41,52,23,46,18,15,29)

>y<-c(8,lA7,9,5,10,3,2,6)

>cor.test(x,y,method='spearman,)

>cor.test(x,y,method='kendair)

有關(guān)系的

5.18

>x<-l:5

>y<-c(rep(x,c(0,l,9,7,3)))

>z<-c(rep(x,c(2,2,HAl)))

>wilcox.test(y,z,exact=F)

結(jié)果顯示這兩種療法沒什么區(qū)別

弟八早

6.1

(1)>snow<-data.frame(X=c(5.1,3.5,7.1,6.2,8.8/7.8/4.5,5,6,8.0,6.4),

+Y=c(1907,1287,2700,2373,3260,3000,1947,2273,3113,2493))

>plot(snow$X,snow$Y)

結(jié)論是有線性關(guān)系的。

(2)(3)

>lm.sol<-lm(Y/"l+X,data=snow);summary(lm.sol)

結(jié)果是方程是顯著的

(4)>predict(lm.sol/data.frame(X=7),interval='prediction'Jevel=0.95)

fitIwrupr

12690.2272454.9712925.484

6.2

(1)(2)

>soil<-data.frame(Xl=c(0.4z0.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=c(52,23,19,34,24,65,44,31,

+29,58,37,46,50,44,56,36,58,51),X3=c(158,163,37,157,59,123,46,117,

+173,112,111,114,134,73,168,143z202z124)zY=c(64,60,71,61,54,77,81,

+93,93,51,76,96,77,93,95,54,168,99))

>lm.sol<-lm(Y~l+X14-X2+X3,data=soil);summary(lm.sol)

我們發(fā)現(xiàn)X2和X3的系數(shù)沒有通過t檢驗(yàn)。但是整個方程通過了檢驗(yàn)。

(3)>lm.ste<-step(lm.sol)

>summary(lm.ste)

可以發(fā)現(xiàn)新模型只含有XI和X3,但是X3的系數(shù)還是不顯著。接下

來考慮用dropl函數(shù)處理

>dropl(lm.ste)

發(fā)現(xiàn)去掉X3殘差升高最小，AIC只是有少量增加。因此應(yīng)該去掉X3

>lm.new<-lm(Y~Xl,data=soil);summary(lm.new)

此時發(fā)現(xiàn)新模型Im.new系數(shù)顯著且方程顯著

6.3

(1)>da<-data.frame(X=c(l,1,1,1,2,2,233,3,4,4,4,5,6,6,6,7,7,7,8,8,8,

+9,11,12,12,12),Y=c(0.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))

>plot(da$X,da$Y)

>lm.sol<-lm(Y~X,data=da)

>abline(lm.sol)

(2)>summary(lm.sol)

全部通過

(3)>plotflm.sol,!)

>windows()

>plot(lm.sol,3)

可以觀察到誤差符合等方差的。但是有殘差異常值點(diǎn)24,27,28.

(4)>lm.up<-update(lm.sol,sqrt(.)~.)

>summary(lm.up)

都通過檢驗(yàn)

>plot(da$X,da$Y)

>abline(lm.up)

>windows()

>plot(lm.up,l)

>windows()

>plot(lm.up,3)

可以發(fā)現(xiàn)還是有殘差離群值24,28

6.4

>lm.sol<-lm(Y^l+Xl+X2,data=toothpaste);summary(lm.sol)

>influence.measures(lm.sol)

>plot(lm.sol,3)

通過influence.measures函數(shù)發(fā)現(xiàn)5,8,9,24對樣本影響較大，可能是

異常值點(diǎn)，而通過殘差圖發(fā)現(xiàn)5是殘差離群點(diǎn)，但是整個殘差還是在

卜2,2］之內(nèi)的。因此可考慮剔除5,8,9,24點(diǎn)再做擬合。

>Im.newv-lm(Y~1+Xl+X2,data二toothpaste,subset=c(-5,-8,-9,-24))

>windows()

>plot(lm.new,3)

>summary(lm.new)

我們發(fā)現(xiàn)Im.new模型的殘差都控制在之內(nèi)，而且方程系數(shù)

和方程本身也都通過檢驗(yàn)。

6.5

>cement<-data.frame(Xl=c(7,1,11,11,7,11,3^1,2,21,1,11,10),

+X2=c(26,29,56,31,52,55,71,31,54,47,40,66,68),

+X3=c(6,15,8,8,6,9,17,22,18,4,23,9,8),

+X4=c(60,52,20,47,33,22,6,44,22,26,34,12,12),

+￥=0(78.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<-cor(cement[l:4])

>kappa(XX,exact=T)

[1]1376.881

>eigen(XX)

發(fā)現(xiàn)變量的多重共線性很強(qiáng)，且有

0.241X1+0.641X2+0.268X3+0.676X4=0

說明X1,X2,X3,X4多重共線。其實(shí)逐步回歸可以解決多重共線的問

題。我們可以檢驗(yàn)一下step函數(shù)去掉變量后的共線性。step去掉了

X3和X4。我們看看去掉他們的共線性如何。

>XX<-cor(cement[l:2])

>kappa(XX,exact=T)

[1]1.59262

我們發(fā)現(xiàn)去掉X3和X4后，條件數(shù)降低好多好多。說明step函數(shù)是

合理的。

6.6

首先得把這個表格看懂。里面的數(shù)字應(yīng)該是有感染和無感染的人數(shù)。

而影響變量有三個。我們把這些影響變量進(jìn)行編碼。如下。

發(fā)生不發(fā)生

抗生素XI23

危險因子

X245

有無計(jì)劃

X367

是否感染Y10

對數(shù)據(jù)的處理，如下

XIX2X3Y頻數(shù)

24611

246017

25610

25602

247111

247087

25710

25700

346128

346030

347123

34703

35618

356032

35710

35709

然后用R處理并求解模型

>hospital<-data.frame(Xl=rep(c(2/2/2/2/2/2/2/2/3/3/3/3/3/3/3/3)/C(l/17/0/2

,11,87,

+0,0,28,30,23,3,8,32,0,9)),X2=rep(c(4,4,5,5,4,4,5,5,4,4,4,4,5,5,5,5),

+6(1,17,0,2,11,87,

+0,0,28,30,23,3,8,32/0z9))zX3=rep(c(6/6z6z6/7/7z7z7/6z6,7,7,6,67,7),

+c(l17,0,2,Cl,8700,28,30,23,3,8,32,0,9)),

+

Y=rep(c(l,0,l,04,0,1,0,1,0,1,0,1,0,1,0),6(1,17,0,2,11,87,0,0,28,30,23,3,8,

32,0,9))

+)

>glm.sol<-glm(Y~Xl+X2+X3,family=binomial’data二hospital)

>summary(glm.sol)

可以發(fā)現(xiàn)如果顯著性為0.1,則方程的系數(shù)和方程本省全部通過檢驗(yàn)。

下面我們來做一個預(yù)測，看看(使用抗生素，有危險因子，有計(jì)劃)

的一個孕婦發(fā)生感染的概率是多少。

>pre<-predict(glm.sol,data.frame(Xl=2,X2=4,X3=6))

>p<-exp(pre)/(l+exp(pre));p

1

0.04240619

即感染的概率為4.2%

6.7

(1)>cofe<-data.frame(X=c(0,0,l,l,2,2,3,344,5,5,6,6),Y=c(508.1,498.4,

+568.2,5773651.7,657,7134,697575537589787.6,792.1,841.4,831.

8))

>lm.sol<-lm(Y~X,data=cofe)

>summary(lm.sol)

(2)>lm.s2<-lm(Y,vX+l(XA2),data=cofe)

>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<-read.csv('6.8.csv'zheader=T)

>glm.sok-glm(Y~Xl+X2+X3+X4+X5,family二binomial,data二pe)

>summary(glm.sol)

可以發(fā)現(xiàn)各變量影響基本都不顯著，甚至大部分還沒通過顯著性檢驗(yàn)。

只有XI的系數(shù)通過了顯著性檢驗(yàn)，但是也不是很理想。下面計(jì)算每

一個病人的生存時間大于200天的概率值。

>pre<-predict(glm.solzdata.frame(Xl=pe$Xl,X2=pe$X2,X3=pe$X3,X4=pe

$X4,X5=pe$X5))

>p<-exp(pre)/(l+exp(pre))

>P

(2)>lm.ste<-step(glm.sol)

結(jié)果是只保留了變量XI和X4o

避免了多重共線性。更加合理一些。下面計(jì)算各個病人的存活概率。

>pre<-predict(lm.ste/data.frame(Xl=pe$Xl,X2=pe$X2,X3=pe$X3/X4=pe$

X4,X5=pe$X5))

>p.new<-exp(pre)/(l+exp(pre))

>p.new

顯然經(jīng)過逐步回歸后的模型更合理。用summary(lm.ste)<,第二個模

型通過了顯著性檢驗(yàn)(a=0.1)

6.9

(1)首先將公式線性化，對方程兩邊直接取對數(shù)即可。然后將得到的

方程用Im回歸。

>peo<-data.frame(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))

>lm.sol<-lm(log(Y)~l+X,data=peo);summary[lm.sol)

Coefficients:

EstimateStd.ErrortvaluePr(>|t|)

(Intercept)4.0371590.08410348.005.08e-16***

X-0.0379740.002284-16.623.86e-10***

>lm.sum<-summary(lm.sol)

>exp(lm.sum$coefficients[l,l])

[1]56.66512

所以theta0=56.66512,thetal=-0.0379

(2)>nls.sol<-nls(Y~b0*exp(bl*X),data=peo,start=list(b0=50/bl=0))

>summary(nls.sol)

Parameters:

EstimateStd.ErrortvaluePr(>|t|)

bO58.6065351.47216039.815.70e-15***

bl-0.0395860.001711-23.136.01e-12***

發(fā)現(xiàn)所求的基本上與內(nèi)在線性相同。

第七章

7.1

(l)>pro<-data.frame(Y=c(115/116,98,83,103,107,118,116,73,89,85,97),

+X=factor(rep(l:3/rep(4,3))))

>pro.aov<-aov(Y~X,data=pro)

>summary(pro.aov)

可以看到不同工廠對產(chǎn)品的影響是顯著的

(2)首先自己編寫求均值的小程序如下

EI

>K<-matrix(0,nrow=l,ncol=3,dimnames=list('mean'/c('R72i7W')))

>for(iin1:3)

+K[lJ]<-mean(pro$Y[pro$X==i])

>K

甲乙丙

mean10311186

然后再用t.test來做均值的置信區(qū)間估計(jì)

>pro.jia<-t.test(pro$Y[pro$X==l]);pro.jia

>pro.yi<-t.test(pro$Y[pro$X==2]);pro.yi

>pro.bing<-t.test(pro$Y[pro$X==3]);pro.bing

(3)>pairwise.t.test(pro$Y;pro$X)

12

20.35-

30.130.04

可以看到顯著性主要有乙工廠和內(nèi)工廠造成

7.2

(1)>old<-data.frame(Y=c(20/18,19,17,15,16,13,18,22,17,26,19,2678,

+23,25,24,25/18,22/27,24/12z14),X=factor(rep(l:4,c(10,6,6,2))))

>old.aov<-aov(Y~X,data=old)

>summary(old.aov)

可以發(fā)現(xiàn)影響是非常顯著的。

(2)>pairwise.t.test(old$\；old$X)

直接從結(jié)果就可以發(fā)現(xiàn)國內(nèi)只有以工廠和丙工廠與國外工廠有顯著

差異。而國內(nèi)只有甲乙，甲丙之間存在著顯著差異。

7.3

>rat<-data.frame(X=c(30/27/35,35,29,33,32,36,26,41,33,31,43,45,53,44,

+51,53,54,37,47,57,48,42,82,66,66,86,56,52,76,83,72,73,59,53),

+A=gl(3,12))

>shapiro.test(rat$X[A==l])

>shapiro.test(rat$X[rat$A==2])

>shapiro.test(rat$X[rat$A==3])

>bartlett.test(X^A,data=rat)

可以看到

數(shù)據(jù)符合正態(tài)性但是不是方差齊性的

7.4

>rat<-data.frame(Y=c(2.79,2.69,3.11,3.47,1.77,2.44,2.83,2.52,3.83,

+3.15,4.7,3.97,2.03,2.87,3.65,5.09,5.41,3.47,4.92,4.07,2.18,3.13,3.77,

+4.26),X=gl(3,8))

>rat.aov<-aov(Y~X,data=rat)

>summary(rat.aov)

結(jié)果是顯著的

7.5

>sleep<-data.frame(Y=c(23.1,57.6,10.5,23.6,H.9,54.6,21.0,20.3,22.7,

+53.2,9.7,19.6,13.8,47.1,13.6,23.6,22.5,53.7,10.8,21.1,13.7,39.2,

+13.7,16.3,22.6,53.1,8.3,21.6,13.3,37.0,14.8,14.8),X=gl(4,8))

>sleep.aov<-aov(Y~X,data=sleep)

>summary(sleep.aov)

結(jié)果是不顯著

7.6

(1)>pro<-data.frame(Y=C(4.6,4.3,6.L6.5,6.86.4,6.3,6.7,3.4,3.84.0,3.8,

+4.7,43,3.9,3.5,6.5,7.0),A=gl(3,2/18),B=gl(3,6,18))

,v

>pro.aov<-aov(YA+B+A:Bzdata=pro);summary(pro.aov)

結(jié)果是A和B及其交互作用都是十分顯著的

(2)首先我們要選出最優(yōu)條件組合，由(1)知影響力為

AB>A>B。下面我們來計(jì)算它們各個水平下的均值。首先要交互作用

給找出來。如下

>ab<-function(x,y){

+n<-length(x);z<-rep(O,n)

+for(iinl:n)

+if(x[i]==y[i]){z[i]<-l}else{z[i]<-2}

+factor(z)}

>pro$AB<-ab(pro$A,pro$B)

然后我們開始計(jì)算各個水平的均值，如下

,,,,

>K<-matrix(0,nrow=3,ncol=3,dimnames=list(l:3/c(A,B7AB)))

>for(iin2:4)

+for(jin1:3)

+K[j,i-l]<-mean(pro$Y[pro[i]==j])

>K

ABAB

15.1500005.7833334.933333

24.5333334.6666675.250000

35.7500004.983333NaN

按照影響力越大（即P值越小），我們首先確定AB應(yīng)選擇水平2,即

A和B不等的是最好的。然后選擇A,選擇水平3,那么B只能在1

和2中選擇，需選擇1.于是我們的最優(yōu)組合為A3B1。

下面給出A3B1的點(diǎn)估計(jì)和區(qū)間估計(jì)。

>mean(pro$Y[pro$A==3&pro$B==l])

>t.test(pro$Y[pro$A==3&pro$B==l])

(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<-data.frame(A=gl(3,3)/B=gl(3/l,9)/

+C=factor(c(l,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<-aov(Y~A+B+C,data=rice);summary(rice.aov)

可以看到影響均不顯著，那么我們干脆直接按照各因素水平的均值大

小來取。下面計(jì)算均值

>K<-matrix(0,nrow=3,ncol=3,dimnames=list(l:3,c('品種密度施肥量

')))

>for(iin1:3)

+for(jin1:3)

+K[iJ]<-mean(rice$Y[rice[j]==i])

>K

品種密度施肥量

159.5333362.7333359.05833

256.5500055.2750055.52500

356.1250054.2000057.62500

所以應(yīng)該選品種8號，密度4.5,施肥量Q75

7.8

首先我們繪制出正交試驗(yàn)表格，如下

列號1234567

A*BA*CB*C

試驗(yàn)號AEC*DCB*DA*DD產(chǎn)量

1111111186

2111222295

3122112291

4122221194

5212121291

6212212196

7221122183

S221211288

好吧，表示因?yàn)槎嗔艘粋€因素D不知道怎么排列交互作用了，我上面

排列的也不一定對。此題暫且不做

7.9

首先把正交試驗(yàn)表的結(jié)果那?列給計(jì)算出來。如下

>pro<-matrix(c(1.5,L7,1.3,1.5,l,1.2,l,l,2.5,2.2,3.2,2,2.5,2.5,1.5,2.8,

+1.5,1.8,1.7,1.5,1,2.5,131.5,1.8,151.8,2.2,192.6,232),ncol=4,

+byrow=T)

>pro.mean<-apply(pro,l,rnean)

現(xiàn)在可以輸入正交試驗(yàn)表了，如下

>pro.data<-data.frame(Y=pro.mean,A=gl(2,4),B=gl(2,2,8)/C=gl(2,l,8))

進(jìn)行分析

>pro.aov<-aov(Y/x,A+B+C+A:B+A:C+B:C,data=pro.data);summary(pro.aov)

從分析結(jié)果可以看出，顯著性大小為B>AB>AC,其余均不顯著

下面再計(jì)算出均值，從而就可以依據(jù)顯著性來選擇最優(yōu)參數(shù)了

>ab<-function(x,y){

+n<-length(x);z<-rep(O,n)

+for(iinl:n)

+if(x[i]==y[i]){z[i]<-l}else{z[i]<-2}

+factor(z)}

>pro.data$AB<-ab(pro.data$A,pro.data$B)

>pro.data$AC<-ab(pro.data$A,pro.data$C)

,,,,

>K<-matrix(0,nrow=2,ncol=5,dimnames=list(l:2/c(A7B7C7AB/'AC')))

>for(iin2:6)

+for(jin1:2)

+K[j,i-l]<-mean(pro.data$Y[pro.data[i]==j])

>K

ABCABAC

11.837501.437501.856251.643751.93750

21.806252.206251.787502.000001.70625

依據(jù)顯著性，首先選擇B,選擇Bl。再依據(jù)AB,應(yīng)選擇AB1,也就是

說A和B應(yīng)該是同一水平。那么A就被先選定的B決定了它應(yīng)該選

水平L然后看AC,應(yīng)該選2.也就是說A和C應(yīng)該是不同水平。那么

A選擇1,C必須選擇2.

所以最后的最優(yōu)組合應(yīng)該是A1B1C2

即通用夾具，特殊鑄鐵，留研量0.015

第八章

8.1

>x<-matrix(c(-l.9,3-2,-6.9,10.4,5.2,2,5,2.5,7.3,0,6.8,12.7,0.9,-15.4,

+-12.5廣2.5,1.5,1.3,3.8,6.8,0.2,0.2,-0.1,7.5,0.4,14.6,2.7,8.3,2.1,0.8,

+-4.6,4.3,-1.7,10.9,-2.6,13.1,2.642.8,-2.8,10),ncol=2,byrow=T)

>g<-gl(2,l,20)

>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.fishe「(x[l:10,Lx[ll:20,],c(8.1,2))

得出的結(jié)論都是明天下雨

8.2

>heart<-read.csv('8.2.csv',header二T)

>G<-factor(rep(l:3,c(11,7,5)))

>distinguish.distance(heart,G,var.equal=F)

>distinguish.distance(heart,G,var.equal=T)

>distinguish.bayes(heart,G,p=c(ll/23,7/23,5/23),var.equal=F)

>distinguish.bayes(heartzG,p=c(ll/23,7/23,5/23),var.equal=T)

無論方差相同還是不同，對于距離判別的正確率都是78.2%

而方差不同的貝葉斯判別正確率僅僅為65.2%

方差相同的貝葉斯判別正確率為87%

8.3

(1)>study<-read.csv('8.3.csv',header=T)

>X<-data.frame(xl=study$xl,x2=study$x2/x3=study$x3,s=stud

y$地區(qū))

>d<-dist(X)

>hc.l<-hclust(d,method=,complete,)

>hc.2<-hclust(d,method='average')

>hc.3<-hclust(d,method='centroid,)

>hc.4<-hclust(d,method='ward')

>opar<-par(mfrow=c(2,2))

>plot(hc.l,hang=-l)

>rectl<-rect.hclust(hc.l,k=4)

>plot(hc.2,hang=-l)

>rect2<-rect.hclust(hc.27k=4)

>plot(hc.3,hang=-l)

>rect3<-rect.hclust(hc.3,k=4)

>plot(hc.4,hang=-l)

>rect4<-rect.hclust(hc.4,k=4)

下面是各種方法分類的結(jié)果

>recti

>rect2

>rect3

>rect4

(2)>km<-kmeans(scale(X),4,nstart=20)

>sort(km$clust)

8.4

>coreer<-read.csv(,48名求職者得分.csv',header=T)

>X<-data.frame(xl=coreer$FLzx2=coreer$APBx3=coreer$AA,x4=coreer$L

A/+x5=coreer$SC,x6=coreer$LC,x7=coreer$HON,x8=coreer$SMS,x9=core

er$EXP,

+xlO=coreer$DRV/xll=coreer$AMB,xl2=coreer$GSP,xl3=coreer$POT,

+xl4=coreer$KJ,xl5=coreer$SUIT,s=coreer$ID)

>d<-as.dist(l-cor(X))

>hcl<-hclust(d,method='complete')

>hc2<-hclust(d,method='average')

>hc3<-hclust(d,method='centroid,)

>hc4<-hclust(d,method='ward')

>opar<-par(mfrow=c(2z2))

>plot(hcl,hang=-l)

>rectl<-rect.hclust(hcl,5)

>plot(hc2,hang=-l)

>rect2<-rect.hclust(hc2z5)

>plot(hc3,hang=-l)

>rect3<-rect.hclust(hc3,5)

>plot(hc4,hang=-l)

>rect4<-rect.hclust(hc4,5)

下面打印出分類的結(jié)果

>recti

>rect2

>rect3

>rect4

第九章

9.1

(1)>fac<-read.csv('9.1.csv',header=T)

>fac.pr<-princomp(fac,cor=T);

>summary(fac.pr)

從結(jié)果可以知道前四個主成分的累積貢獻(xiàn)率達(dá)到0.95.至于他們的意

義嘛，這牽涉到經(jīng)濟(jì)學(xué)知識，我不懂。

(2)>apply(pre,2,order)

我們利用以上代碼看每個行業(yè)在各個主成分的排序，是從小到大排列

的。

下面我們只保留前四個主成分，來對這13個行業(yè)進(jìn)