非参数统计(第二版)习题R程序P37.例2.1build.price<-c(36,32,31 ,25,28,36,40,32,41,26,3 5,35,32,87,33,35);build .pricehist(build.price,freq=FA LSE)#直方图lines(density(build.pric e),col="red")#连线#方法一：m<-mean(build.price); m#均值D<-var(build.price)#方差SD<-sd(build.price)#标准差St=(m-37)/(SD/sqrt(len gth(build.price)));t#t统计量计算检验统计量t=[1] -0.1412332 #方法二：t.test(build.price-37)#课本第38页例2.2binom.test(sum(build.p rice<37),length(build.p rice),0.5)#课本40页例2.3P<-2*(1-pnorm(1.96,0, 1));P[1] 0.04999579P1<-2*(1-pnorm(0.790 6,0,1));P1[1] 0.4291774> 例2.4>p<-2*(pnorm(-1.96,0,1 ));p[1] 0.04999579>>p1<-2*(pnorm(-0.9487 ,0,1));p1[1] 0.3427732例2.5（P45）scores<-c(95,89,68,90, 88,60,81,67,60,60,60,6 3,60,92,60,88,88,87,60,73,60,9 7,91,60,83,87,81,90);le ngth(scores)#输入向量求长度ss<-c(scores-80);sst<-0t1<-0for(i in 1:length(ss)){ if (ss[i]<0) t<-t+1#求小于80的个数else t1<-t1+1求大于80的个数}t;t1> t;t1 [1] 13[1] 15binom.test(sum(scores <80),length(scores),0.7 5)p-value =0.001436<0.01Cox-Staut趋势存在性检验P47例2.6year<-1971:2002;year length(year)rain<-c(206,223,235,26 4,229,217,188,204,182, 230,223,227,242,238,207,208,2 16,233,233,274,234,22 7,221,214,226,228,235,237,243,2 40,231,210)length(rain)#(1)该地区前10年降雨量是否变化？t1=0if(client[i]==client[i+1]) t1<-t1else t1<-t1+1}R<-t1+1;R#=12#find rejection region （不写）rl<-1+2*n1*n0/(n1+n 0)*(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(r ead.table("SHUJU39.txt",header =TRUE));shuju39 attach(shuju39) sum.a=0sum.b=0sum.c=0 for (i in 1:length(id)){ if (pinzhong[i]=="A") sum.a<-sum.a+chanlia ng[i]else if(pinzhong[i]=="B") sum.b<-sum.b+chanlia ng[i]elsefuhao<-sum.c<-sum.c+ chanliang[i]}sum.a;sum.b;sum.cma<-sum.a/4mb<-sum.b/4mc<-sum.c/4ma;mb;mcfuhao<-rep("a",12);fuh aofor (i in 1:length(id)){ if (pinzhong[i]=="A" & ((chanliang[i]-ma)>0)) fuhao[i]<-"+"else if(pinzhong[i]=="B" &((chanliang[i]-mb)>0)) fuhao[i]<-"+"else if(pinzhong[i]=="C" & ((chanliang[i]-mc)>0)) fuhao[i]<-"+"else fuhao[i]<-"-"}fuhao#利用上题编程解决检验的随机性n<-length(fuhao);nn1<-sum(fuhao=="+") ;n1n0<-n-n1;n0t1<-0for (i in1:(length(fuhao)-1)){ if(fuhao[i]==fuhao[i+1]) t1<-t1else t1<-t1+1} R<-t1+1;R#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例2.10（P52）library(quadprog)# 不存在叫‘quadprog’这个名字的程辑包library(zoo)# 不存在叫‘zoo’这个名字的程辑包library(tseries)# 不存在叫‘tseries’这个名字的程辑包run1=factor(c(1,1,1,0,r ep(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)# 错误: 没有"runs.test"这个函数Wilcoxon符号秩检验W+在零假设下的精确分布#下面的函数dwilxonfun 用来计算W+分布密度函数，即P(W+=x)的一个参考程序！dwilxonfun=function(N ){a=c(1,1) #when n=1 frequency of W+=1 or o n=1pp=NULL #distribute of all size from 2 to N aa=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/(2^i)#density of Wilcox distribut when size=N }p}N=19 #sample size of expected distribution of W+y<-dwilxonfun(N);y#计算P(W+=x)中的x取值的R参考程序！！dwilxonfun=function(N ){a=c(1,1) #when n=1 frequency of W+=1 or o n=1pp=NULL #distribute of all size from 2 to N aa=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/(2^i)#density of Wilcox distribut when size=N }a}N=19 #sample size of expected distribution of W+y<-dwilxonfun(N);lengt h(y)-1hist(y,freq=FALSE) lines(density(y),col="re d")例2.12（P59）ceo<-c(310,350,370,37 7,389,400,415,425,440, 295,325,296,250,340,298,3 65,375,360,385);length (ceo)#方法一wilcox.test(ceo-320)#方法二ceo.num<-sum(ceo>32 0);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 in1:(length(a)-1)){for (j in(i+1):length(a)){ walsh<-c(walsh,(a[i]+a [j])/2)}}walsh=c(walsh,a)NW=length(walsh);NW median(walsh)2.5单组数据的位置参数置信区间估计(P61)例2.14‘stu<-c(82,53,70,73,103 ,71,69,80,54,38,87,91,62,75,6 5,77);stualpha=0.05rstu<-sort(stu);rstu conff<-NULL;conffn=length(stu);nfor(i in 1:(n-1)){for (j in (i+1):n){ conf=pbinom(j,n,0.5)-p binom(i,n,0.5)if(conf>1-alpha){conff<-c(conff,i,j,conf)}}}confflength(conff) min<-103-38;minc<-seq(1,(length(conff) -1),3);cfor(i in c){col<-c(rstu[conff[i]],rst u[conff[i+1]],conff[i+2 ])min1<-rstu[conff[i+1]] -rstu[conff[i]]if(min1<min){min<-min 1;l<-i}print(col)}col1<-c(rstu[conff[l]],r stu[conff[l+1]],conff[l+ 2]);col1min例2.14“stu<-c(82,53,70,73,103 ,71,69,80,54,38,87,91,62,75,6 5,77);stualpha=0.05n=length(stu);nconf=pbinom(n,n,0.5)-pbinom(0,n,0.5);conf for(k in 1:n){conf=pbinom(n-k,n,0.5) -pbinom(k,n,0.5)if(conf<1-alpha){loc=k-1;break}}print(loc)（剩余的例题参考程序在课本）3.6正态记分检验例2.18baby1<-c(4,6,9,15,31,3 3,36,65,77,88)baby=(baby1-34);baby baby.mean=mean(baby );baby.mean例2.18qiuzhi<-function(x){n=length(x)a=rep(2,n)for (i in 1:n){a[i]=sum(x<=x[i]) }a}fuhao<-function(x,y){ n=length(x)sgn=rep(2,n)for(i in 1:n){if (x[i]>y)sgn[i]=1else if (x[i]==y)sgn[i]=0elsesgn[i]=-1}sgn}n1<-length(baby) babyzhi=qiuzhi(baby)q=(n1+1+babyzhi)/(2* n1+2)babysgn<-fuhao(baby,3 4)babysgn=sign(baby1-3 4);babysgns=qnorm(q,0,1)W<-t(s)%*%babysgn; Wsd<-sum((s*babysgn)^ 2);sdT=W/sd;T2.7分布的一致性检验例2.19shuju1<-data.frame(m onth=c(1:6), customers=c(27,18,15, 24,36,30));shuju1 attach(shuju1)n<-sum(customers);n expect<-rep(1,6)*(1/6) *n;expect x.squ=sum((customers-expect)^2)/25;x.squ#方法一value<-qchisq(1-0.05,le ngth(customers)-1);val ue#方法二pvalue<-1-pchisq(x.squ ,length(customers)-1);p value例2.20shuju2<-data.frame(ch ongshu=c(0:6), zhushu=c(10,24,10,4,1, 0,1));shuju2attach(shuju2)n=sum(zhushu);n lamda<-sum(chongshu *zhushu)/n;lamdap<-dpois(chongshu,lam da);pn*px.squ=sum((zhushu^2) /(n*p))-n;x.squ#方法一value<-qchisq(1-0.05,le ngth(zhushu)-1);value#方法二pvalue<-1-pchisq(x.squ ,length(zhushu)-1);pval ue例2.21shuju3<-c(36,36,37,38, 40,42,43,43,44,45,48,4 8,50,50,51,52,53,54,54,5 6,57,57,57,58,58,58,58, 58,59,60,61,61,61,62,6 2,63,63,65,66,68,68,70, 73,73,75);shuju3n=length(shuju3)n0=sum(shuju3<30);n 0n1=sum(shuju3>30 & shuju3<=40);n1 n2=sum(shuju3>40 & shuju3<=50);n2n3=sum(shuju3>50 & shuju3<=60);n3n4=sum(shuju3>60 & shuju3<=70);n4n5=sum(shuju3>70 & shuju3<=80);n5n6=sum(shuju3>80);n 6nn<-c(n0,n1,n2,n3,n4,n 5,n6);nn #计算45位学生体重分类的频数！shuju3.mean=mean(sh uju3);shuju3.mean shuju3.var=var(shuju3) ;shuju3.varshuju3.sd=sd(shuju3);s huju3.sde0=pnorm(30,shuju3.m ean,shuju3.sd)e1=pnorm(40,shuju3.m ean,shuju3.sd)-pnorm( 30,shuju3.mean,shuju3. sd)e2=pnorm(50,shuju3.m ean,shuju3.sd)-pnorm( 40,shuju3.mean,shuju3. sd)e3=pnorm(60,shuju3.m ean,shuju3.sd)-pnorm( 50,shuju3.mean,shuju3. sd)e4=pnorm(70,shuju3.m ean,shuju3.sd)-pnorm( 60,shuju3.mean,shuju3. sd)e5=pnorm(80,shuju3.m ean,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((nn^2)/(ee ))-n;x.squ#方法一value<-qchisq(1-0.05,le ngth(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,7 7,92,86,76,80,81,75,77,72,81,9 0,84,86,80,68,77,87,76, 77,78,92,75,80,78);healthyks.test(healthy,pnorm,8 0,6)第三章#Brown_Mood中位数#Brown-Mood中位数检验程序BM.test<-function(x,y,a lt){xy<-c(x,y)md.xy<-median(xy) #利用中位数的检验#md.xy<-quantile(xy,0.25) #利用p分位数的检验t<-sum(xy>md.xy)lx<-length(x)ly<-length(y)lxy<-lx+lyA<-sum(x>md.xy) if (alt=="greater"){w<-1-phyper(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") <-c(">MXY"," <MXY","TOTAL")dimnames(conting.tabl e)<-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.test(X,Y,"less")#方法二：XY<-c(X,Y)md.xy<-median(XY)t<-sum(XY>md.xy)lx<-length(X)ly<-length(Y)lxy<-lx+lyA<-sum(X>md.xy)#没有修正时的情形pvalue1<-pnorm(A,lx*t /(lx+ly),sqrt(lx*ly*t*(lx+ly-t)/( lx+ly)^3));pvalue1#修正时的情形pvalue2<-pnorm(A,lx*t /(lx+ly)-0.5,sqrt(lx*ly*t*(lx+ly-t)/( lx+ly)^3));pvalue2 3.2、Wilcoxon-Mann-Whitne y秩和检验#求两样本分别的秩和的程序.Qiuzhi<-function(x,y){ n1<-length(y)yy<-c(x,y)wm=0for(i in 1:n1){ wm=wm+sum(y[i]>yy, 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=50 wxy<-wy-n*(n+1)/2;w xy#=22mean<-m*n/2var<-m*n*(m+n+1)/1 2pvalue<-1-2*pnorm(wx y,mean-0.5,var);pvalue #方法二wilcox.test(weight.high ,weight.low)例3.4Mx-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)-me dian(x1)B=10000Tboot=c(rep(0,1000)) #vector of length Bootstrapfor (i 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 x2 Tboot[i]=median(xx2)-median(xx1)}th<-median(Tboot);th se=sd(Tboot) Normal.conf=c(th+qnor m(0.025,0,1)*se,th-qno rm(0.025,0,1)*se);Nor mal.confPercentile.conf=c(2*th-quantile(Tboot,0.975),2 *th-quantile(Tboot,0.025));Percenti le.confProvotal.conf=c(quantil e(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<-c(zhi,sum(x[i]>=x y))}zhi}引例：x1<-c(48,56,59,61,84,8 7,91,95)x2<-c(2,22,49,78,85,89 ,93,97)zhi_x1=qiuzhi(x1,x2);z hi_x1#zhi_x2=qiuzhi(x2,x1); zhi_x2#var_x1=var(x1);var_x 1#var_x2=var(x2);var_x 2 m=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,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(var 2)))p_value#方法三Z=1/(sqrt(var2))*(M2-mean2+0.5);Z3.4、Moses方差检验qiuzhi<-function(x,y){ xy<-c(x,y)zhi<-NULLfor (i in1:length(x)){ zhi<-c(zhi,sum(x[i]>=x y))}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#随机分组a1=sample(x1,3,F)xx2=NULLfor(i in 1:m1){if(sum(a1==x1[i])==0) xx2=c(xx2,x1[i])}a2=sample(xx2,3,F)xx3=NULLfor(i in 1:(m1-3)){if(sum(a2==xx2[i])== 0) xx3=c(xx3,x1[i])}a3=sample(xx3,3,F)x11=sum((A[1,]-mean( x1))^2);x11x12=sum((A[2,]-mean( x1))^2);x12x13=sum((A[3,]-mean( x1))^2);x13x14=sum((A[4,]-mean( x1))^2);x14SSA<-c(x11,x12,x13,x1 4);SSAB<-matrix(x2[1:9],ncol =3);By11=sum((B[1,]-mean( x2))^2);y11y12=sum((B[2,]-mean( x2))^2);y12 y13=sum((B[3,]-mean( x2))^2);y13SSB<-c(y11,y12,y13);S SBzhi_SSA=qiuzhi(SSA,SS B);zhi_SSAzhi_SSB=qiuzhi(SSB,SS A);zhi_SSBS=sum(zhi_SSA);STM=S-4*(4+1)/2;TM #方法一(查附表4)拒绝域C=(TM<W(0.025,m1,m 2)或者TM>W(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);n mean_AB=m*n/2;mea n_ABvar_AB=m*n*(m+n+1) /12;var_ABp_value=1-pnorm(S,me an_AB,sqrt(var_AB));p_ value第四章4.1、试验设计和方差分析的基本概念回顾#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);y lever<-c("B","A","C","C ","B","A","B","C","A","A ","C","B")x<-factor(lever);xxy<-data.frame(y,x) attach(xy)aov(formula=y~x,data =xy)aov.xy<-aov(formula=y ~x,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.me anssT<-sum((y-y.mean)^ 2);ssT #计算总的平方和x1.mean<-mean(x1)x2.mean<-mean(x2)x3.mean<-mean(x3) sse<-sum(sum((x1-x1. mean)^2),sum((x2-x2.mean)^2), sum((x3-x3.mean)^2)) ;sse #计算误差平方和sst<-ssT-sse;sst #计算组间平方和F<-(sst/2)/(sse/(lengt h(y)-3));F #计算方差分析的F检验统计量#临界值的计算value<-qf(0.95,2,length (y)-3);value#计算p-value值p.value<-1-pf(8,2,lengt h(y)-3);p.value表4.5xueye<-c(8.4,9.4,9.8,1 2.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((xueye-mea n(xueye))^2);sst1a=matrix(xueye,ncol=5 );aquzu<-apply(a,2,sum); quzuchuli<-apply(a,1,sum);c hulik=5b=4ssb=1/4*sum(quzu^2) -sum(quzu)^2/(k*b);ss bsst=1/5*sum(chuli^2)-sum(chuli)^2/(k*b);sst sse=sst1-ssb-sst;sse mssb=ssb/(k-1);mssb msst=sst/(b-1);msst msse=sse/(k*b-k-b+1) ;msseF1=mssb/msse;F1F2=msst/msse;F2value1=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 in1:length(w)){zhi<-c(zhi,sum(w[i]>= xy))}zhi}a<-c(80,203,236,252,2 84,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);az hibzhi=qiuzhi(b,a,c,d);bz hiczhi=qiuzhi(c,a,b,d);czh idzhi=qiuzhi(d,a,b,c);dz hiH=12/(n*(n+1))*(sum (azhi)^2/length(a)+su m(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);p valuemean=c(mean(a),mean (b),mean(c),mean(d))#两两比较的程序bjiao=function(azhi,bzh i,czhi,dzhi){{n=length(c(azhi,bzhi,c zhi,dzhi))av=sum(azhi)/length(a zhi)bv=sum(bzhi)/length(b zhi)se=sqrt(n*(n+1)/12*( 1/length(azhi)+1/lengt h(bzhi)))d=abs(av-bv)dab=d/sehuizong=c(d,se,dab,qn orm(1-0.05,0,1))}huizong}bjiao(azhi,bzhi,czhi,dzhi ) bjiao(czhi,dzhi,azhi,bzhi )4.3、Jonckheere-Terpstra检验例4.5x=c(125,136,116,101,1 05,109)y=c(122,114,131,120,1 19,127)z=c(128,142,128,134,1 35,131,140,129)xm=mean(x);xmym=mean(y);ymzm=mean(z);zmg=c(rep(1,6),rep(2,6),r ep(3,8))tapply(c(x,y,z),g,media n)JT.test(data=t(c(x,y,z)) ,class=g)Wij<-function(x,y){n1=length(y)zhiij<-0for(i in 1:n1){zhiij=zhiij+sum(x<y[i]) +sum(x==y[i])/2}zhiij}w12=Wij(x,y);w12w13=Wij(x,z);w13w23=Wij(y,z);w23#方法一：通过查表决策！# W=sum(w12,w13,w23) ;W#方法二：通过中心极限定律决策！#N=length(c(x,y,z))n1=length(x)n2=length(y)n3=length(z) E=(N^2-sum(n1^2,n2 ^2,n3^2))/4;E #计算J的数学期望#f1=function(n){f=n^2* (2*n+3)}Var=(f1(N)-sum(f1(n1), f1(n2),f1(n3)))/72;Var sd.Var=sqrt(Var);sd.Var #计算J的方差和标准差# z1=(W-E)/sd.Var;z1#可以通过拒绝域来决策# pvalue=2*(1-pnorm(z,0 ,1));pvalue #可以通过pvalue值来决策#例4.6jie=function(x,y){{ jiedian<-NULLshuju<-NULLxy<-c(x,y)y1<-unique(y)for (i in1:length(y1)){if(sum(xy==y1[i])>1){jiedian<-c(jiedian,sum (xy==y1[i]))shuju<-c(shuju,y1[i])}}j=c(jiedian,shuju)#shuju 输出是那些数据打结#}j}a=c(40,35,38,43,44,41) b=c(38,40,47,44,40,42) c=c(48,40,45,43,46,44) jie12=jie(a,b);jie12jie13=jie(a,c);jie13jie23=jie(b,c);jie23#例4.7（P128)qiuzhi=function(x){ n1=length(x)n2=length(unique(x)) zhi=NULLif (n1==n2){for (i in 1:n1){ zhi=c(zhi,sum(x<=x[i]) )}}else{for (i in 1:n1){ zhi=c(zhi,mean(sum(x <x[i],1):sum(x<=x[i])) )}}zhi}jiedian=function(x1){ n1=length(x1)x2=unique(x1)n2=length(x2)jie=NULLfor(i in 1:n2){n=sum(x1==x2[i])if (n>1)jie=c(jie,n)}jie}a1=c(85,82,82,79)a2=c(87,75,86,82)a3=c(90,81,80,76)a4=c(80,75,81,75)zhi1=qiuzhi(a1);zhi1 zhi2=qiuzhi(a2);zhi2 zhi3=qiuzhi(a3);zhi3 zhi4=qiuzhi(a4);zhi4a1=t(matrix(c(zhi1,zhi2 ,zhi3,zhi4),ncol=4));a1 b1=apply(a1,2,sum);b1 b=4 k=4Q=12/(b*k*(k+1))*su m(b1^2)-3*b*(k+1);Q jie1=jiedian(a1);jie1 jie2=jiedian(a2);jie2 jie3=jiedian(a3);jie3 jie4=jiedian(a4);jie4 jien=c(jie1,jie2,jie3,jie 4);jienjienn=sum(jien^3-jien) ;jiennt1=b*k*(k^2-1);t1Qc=Q/(1-jienn/t1);Qc5.3、Fisher精确性检验setwd("")getwd()例5.3medicine<-matrix(c(8,2 ,7,23),,2,byrow=T) medicinefisher.test(medicine) chisq.test(medicine)。
