非參數(shù)統(tǒng)計(jì)(第二版)習(xí)題R程序

1、P37.例2.1build.price-c(36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35);build.pricehist(build.price,freq=FALSE)#直方圖lines(density(build.price),col=red)#連線#方法一：m-mean(build.price);m#均值D-var(build.price)#方差SD-sd(build.price)#標(biāo)準(zhǔn)差St=(m-37)/(SD/sqrt(length(build.price);t#t統(tǒng)計(jì)量計(jì)算檢驗(yàn)統(tǒng)計(jì)量t=1 -0.1412332#方法二：t.tes

2、t(build.price-37)#課本第38頁(yè)例2.2binom.test(sum(build.price37),length(build.price),0.5)#課本40頁(yè)例2.3P-2*(1-pnorm(1.96,0,1);P1 0.04999579P1 例2.4 p p1-2*(pnorm(-0.9487,0,1);p11 0.3427732例2.5（P45）scores-c(95,89,68,90,88,60,81,67,60,60,60,63,60,92,60,88,88,87,60,73,60,97,91,60,83,87,81,90);length(scores)#輸入向量求長(zhǎng)

3、度ss-c(scores-80);sst-0t1-0for(i in 1:length(ss)if (ssi0) t-t+1#求小于80的個(gè)數(shù)else t1 t;t11 131 15binom.test(sum(scores80),length(scores),0.75)p-value = 0.0014360.01Cox-Staut趨勢(shì)存在性檢驗(yàn)P47例2.6year-1971:2002;yearlength(year)rain-c(206,223,235,264,229,217,188,204,182,230,223,227,242,238,207,208,216,233,233,274,2

4、34,227,221,214,226,228,235,237,243,240,231,210)length(rain)#(1)該地區(qū)前10年降雨量是否變化？t1=0for (i in 1:5)if (rainiraini+5) t1-t1+1t1k-0:t1-1sum(dbinom(k,5,0.5)# =0.1875y-6/(25);y# =0.1875#(2)該地區(qū)前32年降雨量是否變化？t=0for (i in 1:16)if (rainiraini+16) t-t+1tk1-0:min(t,16-t)-1sum(dbinom(k1,16,0.5)# =0.0002593994pbinom

5、(max(k1),16,0.5)#= 0.0002593994y1-(1+16)/(216);y1#=0.0002593994plot(year,rain)abline(v=(1971+2002)/2,col=2)lines(year,rain)anova(lm(rain(year)隨機(jī)游程檢驗(yàn)（P50）例2.8client-c(F,M,M,M,M,M,F,M,M,F,M,M,M,M,F,M,F,M,M,M,F,F,F,M,M,M);clientn-length(client);nn1-sum(client=M);n1n0-n-n1;n0t1-0for (i in 1:(length(clie

6、nt)-1)if (clienti=clienti+1) t1-t1else t1-t1+1R-t1+1;R#=12#find rejection region（不寫）rl-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0);rlru-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0);ru#=15.33476（課本為ru=17）例2.9shuju39-data.frame(read.table(SHUJU39.txt,header=TRUE);shuju39attach(shuju39)sum.a=0sum.b=0sum.c=0for (i in

7、1:length(id)if (pinzhongi=A) sum.a-sum.a+chanliangielse if (pinzhongi=B) sum.b-sum.b+chanliangielse fuhao-sum.c-sum.c+chanliangisum.a;sum.b;sum.cma-sum.a/4mb-sum.b/4mc-sum.c/4ma;mb;mcfuhao0) fuhaoi0) fuhaoi0) fuhaoi-+else fuhaoi-fuhao#利用上題編程解決檢驗(yàn)的隨機(jī)性n-length(fuhao);nn1-sum(fuhao=+);n1n0-n-n1;n0t1-0fo

8、r (i in 1:(length(fuhao)-1)if (fuhaoi=fuhaoi+1) t1-t1else t1-t1+1R-t1+1;R#find rejection regionrl-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0);rlru-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0);ru例2.10（P52）library(quadprog)# 不存在叫quadprog這個(gè)名字的程輯包library(zoo)# 不存在叫zoo這個(gè)名字的程輯包library(tseries)# 不存在叫tseries這個(gè)名字的程輯包run1=f

9、actor(c(1,1,1,0,rep(1,7),0,1,1,0,0,rep(1,6),0,rep(1,4),0,rep(1,5),rep(0,4),rep(1,13)；run1y=factor(run1)runs.test(y)# 錯(cuò)誤: 沒(méi)有runs.test這個(gè)函數(shù)Wilcoxon符號(hào)秩檢驗(yàn)W+在零假設(shè)下的精確分布#下面的函數(shù)dwilxonfun用來(lái)計(jì)算W+分布密度函數(shù)，即P(W+=x)的一個(gè)參考程序！dwilxonfun=function(N)a=c(1,1) #when n=1 frequency of W+=1 or on=1pp=NULL #distribute of all s

10、ize from 2 to Naa=NULL #frequency of all size from 2 to N for (i in 2:N) t=c(rep(0,i),a) a=c(a,rep(0,i)+t p=a/(2i) #density of Wilcox distribut when size=N pN=19 #sample size of expected distribution of W+y-dwilxonfun(N);y#計(jì)算P(W+=x)中的x取值的R參考程序！dwilxonfun=function(N)a=c(1,1) #when n=1 frequency of W+

11、=1 or on=1pp=NULL #distribute of all size from 2 to Naa=NULL #frequency of all size from 2 to N for (i in 2:N) t=c(rep(0,i),a) a=c(a,rep(0,i)+t p=a/(2i) #density of Wilcox distribut when size=N aN=19 #sample size of expected distribution of W+y-dwilxonfun(N);length(y)-1hist(y,freq=FALSE)lines(densit

12、y(y),col=red)例2.12（P59）ceo-c(310,350,370,377,389,400,415,425,440,295,325,296,250,340,298,365,375,360,385);length(ceo)#方法一wilcox.test(ceo-320)#方法二ceo.num320);ceo.numn=length(ceo)binom.test(ceo.num,n,0.5)例2.13(P61)a-c(62,70,74,75,77,80,83,85,88)walsh-NULLfor (i in 1:(length(a)-1) for (j in (i+1):lengt

13、h(a) walsh-c(walsh,(ai+aj)/2) walsh=c(walsh,a)NW=length(walsh);NWmedian(walsh)2.5單組數(shù)據(jù)的位置參數(shù)置信區(qū)間估計(jì)(P61)例2.14stu-c(82,53,70,73,103,71,69,80,54,38,87,91,62,75,65,77);stualpha=0.05rstu-sort(stu);rstuconff1-alpha)conff-c(conff,i,j,conf) confflength(conff)min-103-38;minc-seq(1,(length(conff)-1),3);cfor(i i

14、n c)col-c(rstuconffi,rstuconffi+1,conffi+2)min1-rstuconffi+1-rstuconffiif (min1min)min-min1;l-iprint(col)col1-c(rstuconffl,rstuconffl+1,conffl+2);col1min例2.14“stu-c(82,53,70,73,103,71,69,80,54,38,87,91,62,75,65,77);stualpha=0.05n=length(stu);nconf=pbinom(n,n,0.5)-pbinom(0,n,0.5);conffor(k in 1:n) co

15、nf=pbinom(n-k,n,0.5)-pbinom(k,n,0.5) if (conf1-alpha)loc=k-1;break print(loc)（剩余的例題參考程序在課本）3.6正態(tài)記分檢驗(yàn)例2.18baby1-c(4,6,9,15,31,33,36,65,77,88)baby=(baby1-34);babybaby.mean=mean(baby);baby.mean例2.18qiuzhi-function(x) n=length(x) a=rep(2,n) for (i in 1:n) ai=sum(x=xi) afuhaoy) sgni=1 else if (xi=y) sgni

16、=0 else sgni=-1 sgnn1-length(baby)babyzhi=qiuzhi(baby)q=(n1+1+babyzhi)/(2*n1+2)babysgn-fuhao(baby,34)babysgn=sign(baby1-34);babysgns=qnorm(q,0,1)W-t(s)%*%babysgn;Wsd-sum(s*babysgn)2);sdT=W/sd;T2.7分布的一致性檢驗(yàn)例2.19shuju1-data.frame(month=c(1:6),customers=c(27,18,15,24,36,30);shuju1attach(shuju1)n-sum(cus

17、tomers);nexpect-rep(1,6)*(1/6)*n;expectx.squ=sum(customers-expect)2)/25;x.squ#方法一value-qchisq(1-0.05,length(customers)-1);value#方法二pvalue-1-pchisq(x.squ,length(customers)-1);pvalue例2.20shuju2-data.frame(chongshu=c(0:6),zhushu=c(10,24,10,4,1,0,1);shuju2attach(shuju2)n=sum(zhushu);nlamda-sum(chongshu*

18、zhushu)/n;lamdap-dpois(chongshu,lamda);pn*px.squ=sum(zhushu2)/(n*p)-n;x.squ#方法一value-qchisq(1-0.05,length(zhushu)-1);value#方法二pvalue-1-pchisq(x.squ,length(zhushu)-1);pvalue例2.21shuju3-c(36,36,37,38,40,42,43,43,44,45,48,48,50,50,51,52,53,54,54,56,57,57,57,58,58,58,58,58,59,60,61,61,61,62,62,63,63,65,

19、66,68,68,70,73,73,75);shuju3n=length(shuju3)n0=sum(shuju330 & shuju340 & shuju350 & shuju360 & shuju370 & shuju380);n6nn-c(n0,n1,n2,n3,n4,n5,n6);nn #計(jì)算45位學(xué)生體重分類的頻數(shù)！shuju3.mean=mean(shuju3);shuju3.meanshuju3.var=var(shuju3);shuju3.varshuju3.sd=sd(shuju3);shuju3.sde0=pnorm(30,shuju3.mean,shuju3.sd)e1=

20、pnorm(40,shuju3.mean,shuju3.sd)-pnorm(30,shuju3.mean,shuju3.sd)e2=pnorm(50,shuju3.mean,shuju3.sd)-pnorm(40,shuju3.mean,shuju3.sd)e3=pnorm(60,shuju3.mean,shuju3.sd)-pnorm(50,shuju3.mean,shuju3.sd)e4=pnorm(70,shuju3.mean,shuju3.sd)-pnorm(60,shuju3.mean,shuju3.sd)e5=pnorm(80,shuju3.mean,shuju3.sd)-pnor

21、m(70,shuju3.mean,shuju3.sd)e6=1-pnorm(80,shuju3.mean,shuju3.sd)e=c(e0,e1,e2,e3,e4,e5,e6);eee=n*c(e0,e1,e2,e3,e4,e5,e6);eex.squ=sum(nn2)/(ee)-n;x.squ#方法一value-qchisq(1-0.05,length(ee)-1);value#方法二pvalue-1-pchisq(x.squ,length(ee)-1);pvalue例2.22healthy-c(87,77,92,68,80,78,84,77,81,80,80,77,92,86,76,80,

22、81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,75,80,78);healthyks.test(healthy,pnorm,80,6)第三章#Brown_Mood中位數(shù)#Brown-Mood中位數(shù)檢驗(yàn)程序BM.test-function(x,y,alt) xy-c(x,y) md.xy-median(xy) #利用中位數(shù)的檢驗(yàn) #md.xy-quantile(xy,0.25) #利用p分位數(shù)的檢驗(yàn) tmd.xy) lx-length(x) ly-length(y) lxy-lx+ly Amd.xy) if (alt=greater) w-1-p

23、hyper(A,lx,ly,t) else if (alt=less) w-phyper(A,lx,ly,t) conting.table=matrix(c(A,lx-A,lx,t-A,ly-(t-A),ly,t,lxy-t,lxy),3,3) -c(X,Y,X+Y) MXY,MXY,TOTAL) dimnames(conting.table)-list(, ) list(contingency.table=conting.table,p.vlue=w)例3.2X-c(698,688,675,656,655,648,640,639

24、,620)Y-c(780,754,740,712,693,680,621)#方法一：BM.test(X,Y,less)#方法二：XY-c(X,Y)md.xy-median(XY)tmd.xy)lx-length(X)ly-length(Y)lxy-lx+lyAmd.xy)#沒(méi)有修正時(shí)的情形pvalue1-pnorm(A,lx*t/(lx+ly),sqrt(lx*ly*t*(lx+ly-t)/(lx+ly)3);pvalue1#修正時(shí)的情形pvalue2-pnorm(A,lx*t/(lx+ly)-0.5,sqrt(lx*ly*t*(lx+ly-t)/(lx+ly)3);pvalue23.2、Wi

25、lcoxon-Mann-Whitney秩和檢驗(yàn)#求兩樣本分別的秩和的程序.Qiuzhi-function(x,y)n1-length(y)yyyy,1) wm例3.3weight.low=c(134,146,104,119,124,161,107,83,113,129,97,123)m=length(weight.low)weight.high=c(70,118,101,85,112,132,94)n=length(weight.high)#方法一：wy-Qiuzhi(weight.low,weight.high)#wy=50wxy-wy-n*(n+1)/2;wxy#=22mean-m*n/2

26、var-m*n*(m+n+1)/12pvalue-1-2*pnorm(wxy,mean-0.5,var);pvalue#方法二wilcox.test(weight.high,weight.low)例3.4 Mx-My的R參考程序：x1-c(140,147,153,160,165,170,171,193)x2-c(130,135,138,144,148,155,168)n1-length(x1)n2-length(x2)th.hat-median(x2)-median(x1)B=10000Tboot=c(rep(0,1000) #vector of length Bootstrapfor (i

27、in 1:B)xx1=sample(x1,5,T) #sample of size n1 with replacement from x1xx2=sample(x2,5,T) #sample of size n2 with replacement from x2Tbooti=median(xx2)-median(xx1)th-median(Tboot);thse=sd(Tboot)Normal.conf=c(th+qnorm(0.025,0,1)*se,th-qnorm(0.025,0,1)*se);Normal.confPercentile.conf=c(2*th-quantile(Tboo

28、t,0.975),2*th-quantile(Tboot,0.025);Percentile.confProvotal.conf=c(quantile(Tboot,0.025),quantile(Tboot,0.975);Provotal.confth.hat3.3、Mood方差檢驗(yàn)qiuzhi-function(x,y)xy-c(x,y)zhi-NULLfor (i in 1:length(x)zhi=xy)zhi引例：x1-c(48,56,59,61,84,87,91,95)x2-c(2,22,49,78,85,89,93,97)zhi_x1=qiuzhi(x1,x2);zhi_x1#zh

29、i_x2=qiuzhi(x2,x1);zhi_x2#var_x1=var(x1);var_x1#var_x2=var(x2);var_x2m=length(x1);mn=length(x2);nmean_R=(m+n+1)/2;mean_Rmean1=m*(m+n+1)*(m+n-1)/12;mean1var1=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var1M1=sum(zhi_x1-mean_R)2);M1p_value=2*pnorm(M1,mean1-0.5,sqrt(var1)p_value例3.5X-c(4.5,6.5,7,10,12)Y-c(6,7.2,8

30、,9,9.8)zhi_X=qiuzhi(X,Y);zhi_Xm=length(X);mn=length(Y);nmean_R=(m+n+1)/2;mean_Rmean2=m*(m+n+1)*(m+n-1)/12;mean2var2=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var2M2=sum(zhi_X-mean_R)2);M2#方法一：查附表9#方法二：p_value=2*(1-pnorm(M2,mean2-0.5,sqrt(var2)p_value#方法三Z=1/(sqrt(var2)*(M2-mean2+0.5);Z3.4、Moses方差檢驗(yàn)qiuzhi-func

31、tion(x,y) xy-c(x,y) zhi-NULL for (i in 1:length(x) zhi=xy) zhi例3.6x1-c(8.2,10.7,7.5,14.6,6.3,9.2,11.9,5.6,12.8,5.2,4.9,13.5)m1=length(x1);m1x2-c(4.7,6.3,5.2,6.8,5.6,4.2,6.0,7.4,8.1,6.5)m2=length(x2);m2A-matrix(x1,ncol=3);A#隨機(jī)分組a1=sample(x1,3,F)xx2=NULLfor(i in 1:m1)if(sum(a1=x1i)=0) xx2=c(xx2,x1i)a2

32、=sample(xx2,3,F)xx3=NULLfor(i in 1:(m1-3)if(sum(a2=xx2i)=0) xx3=c(xx3,x1i)a3=sample(xx3,3,F)x11=sum(A1,-mean(x1)2);x11x12=sum(A2,-mean(x1)2);x12x13=sum(A3,-mean(x1)2);x13x14=sum(A4,-mean(x1)2);x14SSA-c(x11,x12,x13,x14);SSAB-matrix(x21:9,ncol=3);By11=sum(B1,-mean(x2)2);y11y12=sum(B2,-mean(x2)2);y12y1

33、3=sum(B3,-mean(x2)2);y13SSB-c(y11,y12,y13);SSBzhi_SSA=qiuzhi(SSA,SSB);zhi_SSAzhi_SSB=qiuzhi(SSB,SSA);zhi_SSBS=sum(zhi_SSA);STM=S-4*(4+1)/2;TM#方法一(查附表4)拒絕域C=(TMW(0.975,m1,m2)其中W(0.975,m1,m2)=m1*m2-W(0.025,m1,m2).#方法二（Wilcoxon秩和檢驗(yàn)）wilcox.test(SSA,SSB)#方法二（Mann-Whitney秩和檢驗(yàn)）m=length(SSA);mn=length(SSB);

34、nmean_AB=m*n/2;mean_ABvar_AB=m*n*(m+n+1)/12;var_ABp_value=1-pnorm(S,mean_AB,sqrt(var_AB);p_value第四章4.1、試驗(yàn)設(shè)計(jì)和方差分析的基本概念回顧#R軟件中單因素方差分析的函數(shù)例4.1#方法一：*Analysis of Variance Model *y-c(2.0,1.4,2.0,2.8,2.4,1.9,1.8,2.5,2.0,1.5,2.1,2.2);ylever-c(B,A,C,C,B,A,B,C,A,A,C,B)x-factor(lever);xxy-data.frame(y,x)attach(

35、xy)aov(formula=yx,data=xy)aov.xy-aov(formula=yx,data=xy)summary(aov.xy)#方法二：x1-c(1.4,1.9,2.0,1.5)x2-c(2.0,2.4,1.8,2.2)x3-c(2.6,2.8,2.5,2.1)y-c(x1,x2,x3);yy.mean-mean(y);y.meanssT-sum(y-y.mean)2);ssT #計(jì)算總的平方和x1.mean-mean(x1)x2.mean-mean(x2)x3.mean-mean(x3)sse-sum(sum(x1-x1.mean)2),sum(x2-x2.mean)2),s

36、um(x3-x3.mean)2);sse #計(jì)算誤差平方和sst-ssT-sse;sst #計(jì)算組間平方和F-(sst/2)/(sse/(length(y)-3);F #計(jì)算方差分析的F檢驗(yàn)統(tǒng)計(jì)量#臨界值的計(jì)算value-qf(0.95,2,length(y)-3);value#計(jì)算p-value值p.value-1-pf(8,2,length(y)-3);p.value表4.5xueye-c(8.4,9.4,9.8,12.2,10.8,15.2,9.8,14.4,8.6,9.8,10.2,9.8,8.8,9.8,8.9,12.0,8.4,9.2,8.5,9.5);xueyesst1-sum(

37、xueye-mean(xueye)2);sst1a=matrix(xueye,ncol=5);aquzu-apply(a,2,sum);quzuchuli-apply(a,1,sum);chulik=5b=4ssb=1/4*sum(quzu2)-sum(quzu)2/(k*b);ssbsst=1/5*sum(chuli2)-sum(chuli)2/(k*b);sstsse=sst1-ssb-sst;ssemssb=ssb/(k-1);mssbmsst=sst/(b-1);msstmsse=sse/(k*b-k-b+1);msseF1=mssb/msse;F1F2=msst/msse;F2val

38、ue1=qf(1-0.05,k-1,k*b-k-b+1)value2=qf(1-0.05,b-1,k*b-k-b+1)例4.3qiuzhi-function(w,x,y,z) xy-c(w,x,y,z) zhi-NULL for (i in 1:length(w) zhi=xy) zhia-c(80,203,236,252,284,368,457,393)b-c(133,180,100,160)c-c(156,295,320,448,465,481,279)d-c(194,214,272,330,386,475)azhi=qiuzhi(a,b,c,d);azhibzhi=qiuzhi(b,a,

39、c,d);bzhiczhi=qiuzhi(c,a,b,d);czhidzhi=qiuzhi(d,a,b,c);dzhiH=12/(n*(n+1)*(sum(azhi)2/length(a)+sum(bzhi)2/length(b)+sum(czhi)2/length(c)+sum(dzhi)2/length(d)-(3*(n+1)方法一：value=qchisq(1-0.05,3);value方法二：pvalue=1-pchisq(H,3);pvaluemean=c(mean(a),mean(b),mean(c),mean(d)#兩兩比較的程序bjiao=function(azhi,bzhi,

40、czhi,dzhi) n=length(c(azhi,bzhi,czhi,dzhi) av=sum(azhi)/length(azhi) bv=sum(bzhi)/length(bzhi) se=sqrt(n*(n+1)/12*(1/length(azhi)+1/length(bzhi) d=abs(av-bv) dab=d/se huizong=c(d,se,dab,qnorm(1-0.05,0,1) huizongbjiao(azhi,bzhi,czhi,dzhi)bjiao(czhi,dzhi,azhi,bzhi)4.3、Jonckheere-Terpstra檢驗(yàn)例4.5x=c(125,

41、136,116,101,105,109)y=c(122,114,131,120,119,127)z=c(128,142,128,134,135,131,140,129)xm=mean(x);xmym=mean(y);ymzm=mean(z);zmg=c(rep(1,6),rep(2,6),rep(3,8)tapply(c(x,y,z),g,median)JT.test(data=t(c(x,y,z),class=g)Wij-function(x,y)n1=length(y)zhiij-0 for(i in 1:n1) zhiij=zhiij+sum(xyi)+sum(x=yi)/2 zhiijw12=Wij(x,y);w12w13