P37.例2.1
build.price<-c(36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35);build.price
hist(build.price,freq=FALSE)#直方图 lines(density(build.price),col=\连线 #方法一：m<-mean(build.price);m#均值 D<-var(build.price)#方差 SD<-sd(build.price)#标准差S
t=(m-37)/(SD/sqrt(length(build.price)));t#t统计量计算检验统计量 t=
[1] -0.1412332
#方法二：t.test(build.price-37)#课本第38页 例2.2
binom.test(sum(build.price<37),length(build.price),0.5)#课本40页 例2.3
P<-2*(1-pnorm(1.96,0,1));P [1] 0.04999579
P1<-2*(1-pnorm(0.7906,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,63,60,92,
60,88,88,87,60,73,60,97,91,60,83,87,81,90);length(scores)#输入向量求长度 ss<-c(scores-80);ss t<-0 t1<-0
for(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] 15
binom.test(sum(scores<80),length(scores),0.75) p-value = 0.001436<0.01 Cox-Staut趋势存在性检验P47 例2.6
year<-1971:2002;year length(year)
rain<-c(206,223,235,264,229,217,188,204,182,230,223, 227,242,238,207,208,216,233,233,274,234,227,221,214,
226,228,235,237,243,240,231,210) length(rain)
#(1)该地区前10年降雨量是否变化？ t1=0 for (i in 1:5){
if (rain[i] sum(dbinom(k,5,0.5))#=0.1875 y<-6/(2^5);y#=0.1875 #(2)该地区前32年降雨量是否变化？ t=0 for (i in 1:16){ if (rain[i] k1<-0:min(t,16-t)-1 sum(dbinom(k1,16,0.5))#=0.0002593994 pbinom(max(k1),16,0.5)#= 0.0002593994 y1<-(1+16)/(2^16);y1#=0.0002593994 plot(year,rain) abline(v=(1971+2002)/2,col=2) lines(year,rain) anova(lm(rain~(year))) 随机游程检验（P50） 例2.8 client<-c(\\\n<-length(client);n n1<-sum(client==\n0<-n-n1;n0 t1<-0 for (i in 1:(length(client)-1)){ if (client[i]==client[i+1]) t1<-t1 else t1<-t1+1 } R<-t1+1;R#=12 #find rejection region（不写） rl<-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0));rl ru<-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0));ru#=15.33476（课本为ru=17） 例2.9 shuju39<-data.frame(read.table (\attach(shuju39) sum.a=0 sum.b=0 sum.c=0 for (i in 1:length(id)){ if (pinzhong[i]==\else if (pinzhong[i]==\sum.b+chanliang[i] elsefuhao<-sum.c<-sum.c+chanliang[i] } sum.a;sum.b;sum.c ma<-sum.a/4 mb<-sum.b/4 mc<-sum.c/4 ma;mb;mc fuhao<-rep(\for (i in 1:length(id)){ if (pinzhong[i]==\fuhao[i]<-\ else if (pinzhong[i]==\fuhao[i]<-\ else if (pinzhong[i]==\fuhao[i]<-\elsefuhao[i]<-\} fuhao #利用上题编程解决检验的随机性 n<-length(fuhao);n n1<-sum(fuhao==\n0<-n-n1;n0 t1<-0 for (i in 1:(length(fuhao)-1)){ if (fuhao[i]==fuhao[i+1]) t1<-t1 else t1<-t1+1 } R<-t1+1;R #find rejection region rl<-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0));rl ru<-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,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)))；run1 y=factor(run1) runs.test(y)# 错误: 没有\这个函数 Wilcoxon符号秩检验 W+在零假设下的精确分布 #下面的函数dwilxonfun用来计算W+分布密度函数，即P(W+=x)的一个参考程序！ dwilxonfun=function(N){ a=c(1,1) #when n=1 frequency of W+=1 or o n=1 pp=NULL #distribute of all size from 2 to N aa=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/(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=1 pp=NULL #distribute of all size from 2 to N aa=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/(2^i) #density of Wilcox distribut when size=N } a } N=19 #sample size of expected distribution of W+ y<-dwilxonfun(N);length(y)-1 hist(y,freq=FALSE) lines(density(y),col=\例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.num<-sum(ceo>320);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,88) walsh<-NULL for (i in 1:(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,65,77);stu alpha=0.05 rstu<-sort(stu);rstu conff<-NULL;conff n=length(stu);n for(i in 1:(n-1)){ for (j in (i+1):n){ conf=pbinom(j,n,0.5)-pbinom(i,n,0.5) if (conf>1-alpha){conff<-c(conff,i,j,conf)} } } conff length(conff) min<-103-38;min c<-seq(1,(length(conff)-1),3);c for(i in c){ col<-c(rstu[conff[i]],rstu[conff[i+1]],conff[i+2]) min1<-rstu[conff[i+1]]-rstu[conff[i]] if (min1 print(col) } col1<-c(rstu[conff[l]],rstu[conff[l+1]],conff[l+2]);col1 min 例2.14“ stu<-c(82,53,70,73,103,71,69, 80,54,38,87,91,62,75,65,77);stu alpha=0.05 n=length(stu);n conf=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.18 baby1<-c(4,6,9,15,31,33,36,65,77,88) baby=(baby1-34);baby baby.mean=mean(baby);baby.mean 例2.18 qiuzhi<-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]=1 else if (x[i]==y) sgn[i]=0 else sgn[i]=-1 } sgn } n1<-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)%*obysgn;W sd<-sum((s*babysgn)^2);sd T=W/sd;T 2.7分布的一致性检验 例2.19 shuju1<-data.frame(month=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,length(customers)-1);value #方法二 pvalue<-1-pchisq(x.squ,length(customers)-1);pvalue 例2.20 shuju2<-data.frame(chongshu=c(0:6), zhushu=c(10,24,10,4,1,0,1));shuju2 attach(shuju2) n=sum(zhushu);n lamda<-sum(chongshu*zhushu)/n;lamda p<-dpois(chongshu,lamda);p n*p x.squ=sum((zhushu^2)/(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,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);shuju3 n=length(shuju3) n0=sum(shuju3<30);n0 n1=sum(shuju3>30 & shuju3<=40);n1 n2=sum(shuju3>40 & shuju3<=50);n2 n3=sum(shuju3>50 & shuju3<=60);n3 n4=sum(shuju3>60 & shuju3<=70);n4 n5=sum(shuju3>70 & shuju3<=80);n5 n6=sum(shuju3>80);n6 nn<-c(n0,n1,n2,n3,n4,n5,n6);nn #计算45位学生体重分类的频数！ shuju3.mean=mean(shuju3);shuju3.mean shuju3.var=var(shuju3);shuju3.var shuju3.sd=sd(shuju3);shuju3.sd e0=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,shuju3.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);e ee=n*c(e0,e1,e2,e3,e4,e5,e6);ee x.squ=sum((nn^2)/(ee))-n;x.squ #方法一 value<-qchisq(1-0.05,length(ee)-1);value #方法二 pvalue<-1-pchisq(x.squ,length(ee)-1);pvalue 例2.22 healthy<-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 ks.test(healthy,pnorm,80,6) 第三章 #Brown_Mood中位数 #Brown-Mood中位数检验程序 BM.test<-function(x,y,alt){ 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+ly A<-sum(x>md.xy) if (alt==\ {w<-1-phyper(A,lx,ly,t)} else if (alt==\ {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) col.name<-c(\ row.name<-c(\ dimnames(conting.table)<-list(row.name, col.name) list(contingency.table=conting.table,p.vlue=w) } 例3.2 X<-c(698,688,675,656,655,648,640,639,620) Y<-c(780,754,740,712,693,680,621) #方法一： BM.test(X,Y,\#方法二： XY<-c(X,Y) md.xy<-median(XY) t<-sum(XY>md.xy) lx<-length(X) ly<-length(Y) lxy<-lx+ly A<-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-Whitney秩和检验 #求两样本分别的秩和的程序. Qiuzhi<-function(x,y){ n1<-length(y) yy<-c(x,y) wm=0 for(i in 1:n1){ wm=wm+sum(y[i]>yy,1) } wm } 例3.3 weight.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;wxy#=22 mean<-m*n/2 var<-m*n*(m+n+1)/12 pvalue<-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=10000 Tboot=c(rep(0,1000)) #vector of length Bootstrap for (i in 1:B) { xx1=sample(x1,5,T) #sample of size n1 with replacement from x1 xx2=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+qnorm(0.025,0,1)*se,th-qnorm(0.025,0,1)*se);Normal.conf Percentile.conf=c(2*th-quantile(Tboot,0.975),2*th-quantile (Tboot,0.025));Percentile.conf Provotal.conf=c(quantile(Tboot,0.025),quantile(Tboot,0.975));Provotal.conf th.hat 3.3、Mood方差检验 qiuzhi<-function(x,y){ xy<-c(x,y) zhi<-NULL for (i in 1:length(x)){ zhi<-c(zhi,sum(x[i]>=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 #zhi_x2=qiuzhi(x2,x1);zhi_x2 #var_x1=var(x1);var_x1 #var_x2=var(x2);var_x2 m=length(x1);m n=length(x2);n mean_R=(m+n+1)/2;mean_R mean1=m*(m+n+1)*(m+n-1)/12;mean1 var1=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var1 M1=sum((zhi_x1-mean_R)^2);M1 p_value=2*pnorm(M1,mean1-0.5,sqrt(var1)) p_value 例3.5 X<-c(4.5,6.5,7,10,12) Y<-c(6,7.2,8,9,9.8) zhi_X=qiuzhi(X,Y);zhi_X m=length(X);m n=length(Y);n mean_R=(m+n+1)/2;mean_R mean2=m*(m+n+1)*(m+n-1)/12;mean2 var2=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var2 M2=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);Z 3.4、Moses方差检验 qiuzhi<-function(x,y){ xy<-c(x,y) zhi<-NULL for (i in 1:length(x)){ zhi<-c(zhi,sum(x[i]>=xy)) } zhi } 例3.6 x1<-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);m1 x2<-c(4.7,6.3,5.2,6.8,5.6,4.2, 6.0,7.4,8.1,6.5) m2=length(x2);m2 A<-matrix(x1,ncol=3);A#随机分组 a1=sample(x1,3,F) xx2=NULL for(i in 1:m1){ if(sum(a1==x1[i])==0) xx2=c(xx2,x1[i]) } a2=sample(xx2,3,F) xx3=NULL for(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);x11 x12=sum((A[2,]-mean(x1))^2);x12 x13=sum((A[3,]-mean(x1))^2);x13 x14=sum((A[4,]-mean(x1))^2);x14 SSA<-c(x11,x12,x13,x14);SSA B<-matrix(x2[1:9],ncol=3);B y11=sum((B[1,]-mean(x2))^2);y11 y12=sum((B[2,]-mean(x2))^2);y12 y13=sum((B[3,]-mean(x2))^2);y13 SSB<-c(y11,y12,y13);SSB zhi_SSA=qiuzhi(SSA,SSB);zhi_SSA zhi_SSB=qiuzhi(SSB,SSA);zhi_SSB S=sum(zhi_SSA);S TM=S-4*(4+1)/2;TM #方法一(查附表4) 拒绝域C=(TM 其中W(0.975,m1,m2)=m1*m2-W(0.025,m1,m2). #方法二（Wilcoxon秩和检验） wilcox.test(SSA,SSB) #方法二（Mann-Whitney秩和检验） m=length(SSA);m n=length(SSB);n mean_AB=m*n/2;mean_AB var_AB=m*n*(m+n+1)/12;var_AB p_value=1-pnorm(S,mean_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(\ x<-factor(lever);x xy<-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);y y.mean<-mean(y);y.mean ssT<-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/(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.5 xueye<-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);xueye sst1<-sum((xueye-mean(xueye))^2);sst1 a=matrix(xueye,ncol=5);a quzu<-apply(a,2,sum);quzu chuli<-apply(a,1,sum);chuli k=5 b=4 ssb=1/4*sum(quzu^2)-sum(quzu)^2/(k*b);ssb sst=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);msse F1=mssb/msse;F1 F2=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.3 qiuzhi<-function(w,x,y,z){ xy<-c(w,x,y,z) zhi<-NULL for (i in 1:length(w)){ zhi<-c(zhi,sum(w[i]>=xy)) } zhi } a<-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);azhi bzhi=qiuzhi(b,a,c,d);bzhi czhi=qiuzhi(c,a,b,d);czhi dzhi=qiuzhi(d,a,b,c);dzhi H=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);pvalue mean=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/se huizong=c(d,se,dab,qnorm(1-0.05,0,1))} huizong } bjiao(azhi,bzhi,czhi,dzhi) bjiao(czhi,dzhi,azhi,bzhi) 4.3、Jonckheere-Terpstra检验 例4.5 x=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);zm g=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(x w12=Wij(x,y);w12 w13=Wij(x,z);w13 w23=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.6 jie=function(x,y){{ jiedian<-NULL shuju<-NULL xy<-c(x,y) y1<-unique(y) for (i in 1: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);jie12 jie13=jie(a,c);jie13 jie23=jie(b,c);jie23 #例4.7（P128) qiuzhi=function(x){ n1=length(x) n2=length(unique(x)) zhi=NULL if (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 jiedian=function(x1){ n1=length(x1) x2=unique(x1) n2=length(x2) jie=NULL for(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);zhi4 a1=t(matrix(c(zhi1,zhi2,zhi3,zhi4),ncol=4));a1 b1=apply(a1,2,sum);b1 b=4 k=4 Q=12/(b*k*(k+1))*sum(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,jie4);jien jienn=sum(jien^3-jien);jienn t1=b*k*(k^2-1);t1 Qc=Q/(1-jienn/t1);Qc 5.3、Fisher精确性检验 setwd(\getwd() 例5.3 medicine<-matrix(c(8,2,7,23),,2,byrow=T) medicine fisher.test(medicine) chisq.test(medicine)