(完整版)非参数统计(第二版)习题R程序

1、60,88,88,87,60,73,60,97,91,60,83,87,81,90);length(scores)# 输入向量求长度build.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)# 标准差 St=(m-37)

2、/(SD/sqrt(length(build.price);t#t 统计量计算检验统计量t=1 -0.1412332# 方法二： t.test(build.price-37)# 课本第 38 页例 2.2binom.test(sum(build.price37),length(build.price),0.5)# 课本 40 页例 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,

3、60,60,60,63,60,92,ss-c(scores-80);sst-0t1-0for(i in 1:length(ss)if (ssi0) t-t+1# 求小于 80 的个数else t1 t;t11 131 15 binom.test(sum(scores80),length(scores),0.75)p-value = 0.0014360.01Cox-Staut 趋势存在性检验 P47例 2.6year-1971:2002;year length(year)rain-c(206,223,235,264,229,217,188,204,182,230,223,227,242,238,

4、207,208,216,233,233,274,234,227,221 ,214,226,228,235,237,243,240,231,210) length(rain)#(1) 该地区前 10 年降雨量是否变化？t1=0for (i in 1:5)if (rainiraini+5) t1-t1+1t1k-0:t1-1sum(dbinom(k,5,0.5)# =0.1875y-6/(2A5);y# =0.1875P37. 例 2.1#(2) 该地区前 32 年降雨量是否变化？t=0for (i in 1:16)if (rainiraini+16) t-t+1tk1-0:min(t,16-t)

5、-1sum(dbinom(k1,16,0.5)# =0.0002593994pbinom(max(k1),16,0.5)#= 0.0002593994y1-(1 + 16)/(2X6);y1#=0.0002593994plot(year,rain)abline(v=(1971+2002)/2,col=2)lines(year,rain)anova(lm(rain(year)随机游程检验( P50)例 2.8client-c(F,M,M,M,M,M,F,M,n1-sum(client=M);n1n0-n-n1;n0 t1-0 for (i in1:(length(client)-1) if (c

6、lienti=clienti+1)t1-t1 else t1-t1+1R-t1+1;R#=12 #find rejection region (不写)例 2.9shuju39-data.frame(read.table(SHUJU39.txt,header=TRUE);shuju39attach(shuju39)sum.a=0sum.b=0sum.c=0for (i in 1:length(id)if (pinzhongi=A) sum.a-sum.a+chanliangielse if (pinzhongi=B) sum.b-sum.b+chanliangielse fuhao-sum.c-

7、sum.c+chanliangisum.a;sum.b;sum.c ma-sum.a/4 mb-sum.b/4mc-sum.c/4ma;mb;mc fuhao0) fuhaoi0)fuhaoi0)fuhaoi-+else fuhaoi-fuhao# 利用上题编程解决检验的随机性n-length(fuhao);nn1-sum(fuhao=+);n1 n0-n-n1;n0t1-0;clientrl-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=1

8、7 )n-length(client);nfor (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 这个名字的程辑包library(zoo)# 不存在叫 zoo 这个名字的程辑包library(tseri

9、es)# 不存在叫 tseries 这个名字的程 辑包run1=factor(c(1,1,1,0,rep(1,7),0,1,1,0,0,rep(1,6),0,r ep(1,4),0,rep(1,5),rep(0,4),rep(1,13)； run1y=factor(run1)runs.test(y)# 错误 : 没有 runs.test 这个函数Wilcoxon 符号秩检验W+ 在零假设下的精确分布# 下面的函数 dwilxonfun 用来计算 W+ 分布密度函数， 即P(W+=x) 的一个参考程序！dwilxonfun=function(N)a=c(1,1) #when n=1 freque

10、ncy of W+=1 or o n=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)+tp=a/(2Ai) #de nsity of Wilcox distribut whe n size=NpN=19 #sample size of expected distribution of W+y-dwilxonfun(N);y# 计算 P(W+=x) 中的 x 取值的 R 参考程序！d

11、wilxonfun=function(N)a=c(1,1) #when n=1 frequency of W+=1 or on=1pp=NULL #distribute of all size from 2 to Naa=NULL #frequency of all size from 2 to Nfor (i in 2:N)t=c(rep(0,i),a)a=c(a,rep(0,i)+tp=a/(2Ai) #density of Wilcox distribut when size=NaN=19 #sample size of expected distribution of W+y-dwil

12、xonfun(N);length(y)-1 hist(y,freq=FALSE)lines(density(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.num n=length(ceo)binom.test(ceo.num,n,0.5)例 2.13(P61) a-c(62,70,74,75,77,80,83,85,

13、88)walsh-NULLfor (i in 1:(length(a)-1)for (j in(i+1):length(a) walsh-c(walsh,(ai+aj)/2) walsh=c(walsh,a) NW=length(walsh);NWmedian(walsh)2.5 单组数据的位置参数置信区间估计 (P61)例 2.14 stu-c(82,53,70,73,103,71,69,80,54,38,87,91,62,75,65,77);stualpha=0.05 rstu-sort(stu);rstuconff1-alpha)conff-c(conff,i,j,conf)conff

14、length(conff) min-103-38;minc-seq(1,(length(conff)-1),3);c for(i in c)col-c(rstuconffi,rstuconffi+1,conffi+2)min1-rstuconffi+1-rstuconffi if(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=le

15、ngth(stu);nconf=pbinom(n,n,0.5)-pbinom(0,n,0.5);conffor(k in 1:n)conf=pbinom(n-k,n,0.5)-pbinom(k,n,0.5)if (conf1-alpha)loc=k-1;breakprint(loc)(剩余的例题参考程序在课本)3.6 正态记分检验例 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,

16、n)for (i in 1:n)ai=sum(x=xi)afuhaoy)sgni=1else if (xi=y)sgni=0elsesgni=-1sgnn1-length(baby)babyzhi=qiuzhi(baby) q=(n1+1+babyzhi)/(2*n1+2)babysgn-fuhao(baby,34)babysgn=sign(baby1-34);babysgn s=qnorm(q,0,1)W-t(s)%*%babysgn;Wsd-sum(s*babysg n)A2);sdT=W/sd;T2.7 分布的一致性检验例 2.19shuju1-data.frame(month=c(1:

17、6),customers=c(27,18,15,24,36,30);shuju1attach(shuju1)n-sum(customers);nexpect-rep(1,6)*(1/6)*n;expectx.squ=sum(customers-expect)A2)/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

18、,10,4,1,0,1);shuju2attach(shuju2)n=sum(zhushu);nlamda-sum(chongshu*zhushu)/n;lamdap-dpois(chongshu,lamda);p n*px.squ=sum(zhushuA2)/(n *p)-n; x.squ# 方法一value-qchisq(1-0.05,length(zhushu)-1);value# 方法二pvalue-1-pchisq(x.squ,length(zhushu)-1);pvalue例 2.21 shuju3-c(36,36,37,38,40,42,43,43,44,45,48,48,50,

19、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,66,68,68,70,73,73,75);shuju3n=length(shuju3)n0=sum(shuju330 & shuju340 & shuju350 & shuju360 & shuju370 & shuju380);n6 nn-c(n0,n1,n2,n3,n4,n5,n6);nn #计算 45 位学生 体重分类的频数！shuju3.mean=mean(shuju3);shuju3.meanshuju3.var=var(shuju3

20、);shuju3.varshuju3.sd=sd(shuju3);shuju3.sde0=pnorm(30,shuju3.mean,shuju3.sd)e1=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,sh

21、uju3.sd)-pnorm(60,shuju3.mean,shuju3.sd)e5=pnorm(80,shuju3.mean,shuju3.sd)-pnorm(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 (nnA2)/(ee)-n; x.squ# 方法一value-qchisq(1-0.05,length(ee)-1);value# 方法二pvalue-1-pchisq(x.squ

22、,length(ee)-1);pvalue例 2.22healthy-c(87,77,92,68,80,78,84,77,81,80,80,77,92,86,76,80,81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,75,80,78);healthy#md.xy-quantile(xy,0.25) # 利用 p 分位数的检 验tmd.xy)lx-length(x)ly-length(y)lxy-lx+lyAmd.xy)if (alt=greater)w-1-phyper(A,lx,ly,t)else if (alt=less)w-phyper(

23、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,620)Y-c(780,754,740,712,693,680,621)# 方法一：BM.tes

24、t(X,Y,less)# 方法二：XY-c(X,Y)md.xy-median(xy) # 利用中位数的检验ks.test(healthy,pnorm,80,6)第三章#Brown_Mood 中位数#Brown-Mood 中位数检验程序BM.test-function(x,y,alt)xy-c(x,y)tmd.xy)lx-length(X)ly-length(Y)lxy-lx+lyAmd.xy)# 没有修正时的情形pvalue1-pnorm(A,lx*t/(lx+ly),sqrt(lx*ly*t*(lx+ly-t)/(lx+lyF3);pvalue1# 修正时的情形pvalue2-pnorm(A

25、,lx*t/(lx+ly)-0.5,sqrt(lx*ly*t*(lx+ly-t)/(lx+ly)A3);pvalue23.2 、Wilcoxon-Mann-Whitney秩和检验# 求两样本分别的秩和的程序 .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-Q

26、iuzhi(weight.low,weight.high)#wy=50wxy-wy-n*(n+1)/2;wxy#=22mean-m*n/2var-m*n*(m+n+1)/12pvalue-1-2*pnorm(wxy,mean-0.5,var);pvalue# 方法二wilcox.test(weight.high,weight.low)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=10000T

27、boot=c(rep(0,1000)Bootstrapfor (i in 1:B)xx1=sample(x1,5,T)replacement from x1xx2=sample(x2,5,T)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(Tboot,0.975),2*th-qu

28、antile (Tboot,0.025);Percentile.confProvotal.conf=c(quantile(Tboot,0.025),quantile(Tboot,0.975);Provotal.confth.hat3.3 、Mood 方差检验qiuzhi-function(x,y)xy-c(x,y)zhi-NULLfor (i in 1:length(x)zhi=xy)例 3.4Mx-My 的 R 参考程序：#vector of length#sample of size n1 with#sample of size n2 withzhi引例：x1-c(48,56,59,61,

29、84,87,91,95)x2-c(2,22,49,78,85,89,93,97)zhi_x1=qiuzhi(x1,x2);zhi_x1#zhi_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;var1M 仁 sum(zhi_x1-mean_R)A2);M1p_value=2

30、*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,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)A2);M2# 方法一：查附表 9# 方法二：p_value=2*(1-pnorm(M2,mean2-0.5,

31、sqrt(var2)p_value# 方法三x14=sum(A4,-mean(x1)A2);x143.4 、 Moses 方差检验qiuzhi-function(x,y)xy-c(x,y)zhi-NULLfor (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# 随机分组

32、a1=sample(x1,3,F)xx2=NULLfor(i in 1:m1)if(sum(a1=x1i)=0) xx2=c(xx2,x1i)a2=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)A2);x11x12=sum(A2,-mean(x1)A2);x12x13=sum(A3,-mean(x1)A2);x13Z=1/(sqrt(var2)*(M2-mean2+0.5);Zx-factor(lever);x xy-d

33、ata.frame(y,x) attach(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)A2);ssT # 计算总的平方和x1.mean-mean(x1)x2.mean-mean(x2)x3.mean-mean(x3)sse-sum(sum(x1-x1.m

34、ean)A2),sum(x2-x2.mean)A2),sum(x3-x3.mean)A2);sse# 计算误差平方和sst-ssT-sse;sst # 计算组间平方和F-(sst/2)/(sse/(length(y)-3);F # 计算方差分析的 F 检验统计量# 临界值的计算value-qf(0.95,2,length(y)-3);value# 计算 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.

35、9,12.0,8.4,9.2,8.5,9.5);xueyesst1-sum(xueye-mean(xueye)A2);sst1a=matrix(xueye,ncol=5);aSSA-c(x11,x12,x13,x14);SSAB-matrix(x21:9,ncol=3);By1 仁 sum(B1,-mean( x2)A2);y11y12=sum(B2,-mean( x2)A2);y12SSB-c(y11,y12,y13);SSBzhi_SSA=qiuzhi(SSA,SSB);zhi_SSAzhi_SSB=qiuzhi(SSB,SSA);zhi_SSBS=sum(zhi_SSA);STM=S-4

36、*(4+1)/2;TM# 方法一 (查附表 4)拒绝域 C=(TMW(0.975,m1,m2)其中 W(0.975,m1,m2)=m1*m2-W(0.025,m1,m2).# 方法二( Wilcoxon 秩和检验)wilcox.test(SSA,SSB)# 方法二( Mann-Whitney 秩和检验)m=length(SSA);mn=length(SSB);nmean_AB=m*n/2;mean_ABvar_AB=m*n*(m+n+1)/12;var_AB p_value=1-pnorm(S,mean_AB,sqrt(var_AB);p_value第四章4.1 、试验设计和方差分析的基本概念

37、回顾#R 软件中单因素方差分析的函数例 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)quzu-apply(a,2,sum);quzuchuli-apply(a,1,sum);chuli k=5b=4ssb=1/4*sum(quzuA2)-sum(quzu)A2/(k*b);ssbsst=1/5*sum(chulL2)-sum(chuli)A2/(k*b);sstsse=sst1-ssb-sst

38、;ssemssb=ssb/(k-1);mssbmsst=sst/(b-1);msstmsse=sse/(k*b-k-b+1);msse F1=mssb/msse;F1F2=msst/msse;F2 value1=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-NULLfor (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

39、(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,c,d);bzhiczhi=qiuzhi(c,a,b,d);czhidzhi=qiuzhi(d,a,b,c);dzhiH=12/( n*(n +1)*(sum(azhi)A2/length(a)+sum(bzhi )A2/le ngth(b)+sum(czhi)A2/length(c)+sum(dzhi)A2/length(d)-(3*(n+1)方法一： value=qchisq(1-0.05,3);

40、value 方法二：pvalue=1-pchisq(H,3);pvaluemean=c(mean(a),mean(b),mean(c),mean(d)# 两两比较的程序bjiao=function(azhi,bzhi,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/sehuizong=c(d,se,dab,qnorm(1-0

41、.05,0,1)huizongbjiao(azhi,bzhi,czhi,dzhi) bjiao(czhi,dzhi,azhi,bzhi)4.3 、Jonckheere-Terpstra 检验例 4.5x=c(125,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);xm ym=mean(y);ym zm=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)zhiij-0for(i in 1:n1)zhiij=zhiij+sum(xyi)+sum(x=yi)/2zhiijw12=Wij(x,y);w12w