如果有结出现,Uii可稍作变更为U* = #(Xik < Xjl k = l,..., ni, l = l,..., ni)+#(Xik = Xjt k = 1, .,ni, l = 1,., nj)2J* =Zi<;Ut统计量其精确分布类似地可求出样本大的时候J-(N2-h=1n2)/4Z=→N(0, 1)V[N2(2N+3) -Zh=1 n2(2ni +3)]/72Jonckheere-Terpstra检验比Kruskal-Wallis检验有更强的势
其精确分布类似地可求出
例中,123生活方式3.77.39.0一个月后3.75.24.97.1减少的重量3.05.38.7(单位500g)3.95.72.76.5554ni=Ho : 01 = 02 = 03: H1 : 01 ≤ 02 ≤ 03(至少有一个严格不等式成立)容易得出U12=25,U23=14,U13=20,及J=59精确分布方法
例中, (至少有一个严格不等式成立) 精确分布方法
ex.Jonckhexact=function(m1=5,m2=5,m3=4,JTvalue=59)m<-m1+m2+m3;ex5.Jonckhexactch5=function(m)a<-rep(0,5);for (i in 1:(m-4)(for (j in (i+1):(m-3)for (kin(j+1):(m-2))(for ( in (k+1):(m-1)for (f in (I+1):m)(ind<-min((j-i),(k-i),(k-j),(1-i),(1-j),(1-k),(f-i),(f-j),(f-k),(f-I)if (ind>0)a<-rbind(a,c(i,j,k,1,f)))a[2:nrow(a),j)JTid1<-ex5.Jonckhexactch5(m1+m2+m3);n1<-nrow(JTid1)JTid2<-ex5.Jonckhexactch5(m2+m3);n2<-nrow(JTid2)const<-1:m;JT<-matrix(0,nrow=n1*n2,nco/=m+1)for(iin1:n1)(for(jin1:n2)temp1<-c(JTid1[i,J);JTj+(i-1)*n2,1:m1]<-temp1temp2<-(const[-temp1])[c(JTid2j,J)];JT[j+(i-1)*n2,(m1+1):(m1+m2)]<-temp2temp3<-const[-JTj+(i-1)*n2,1:(m2+m1)]]JT[j+(i-1)*n2,(m1+m2+1):(m1+m2+m3))<-temp3JTj+(i-1)*n2,m1+m2+m3+1<-sum(outer(temp2,temp1,">"))+sum(outer(temp3,temp1,">"))+sum(outer(temp3,temp2,">")))y<-JT[,m+1);pvalue<-(sum(y>=JTvalue))/(nrow(JT);hist(y,breaks=min(y):max(y))z<-c(o,hist(y,breaks=min(y):max(y)$counts)list(cbind("m1"=m1,"m2"=m2,"m3"=m3,"JTvalue"=JTvalue),cbind("n1"=n1,"n2"=n2,"n1timesn2"=nrow(JT)c("min"=min(y),"max"=max(y),"JTvalue"=JTvalue,pvalue),cbind(min(y):max(y),z,cumsum(z),cumsum(z)/(n1*n2),1-cumsum(z)/(n1*n2)
ex.Jonckhexact=function(m1=5,m2=5,m3=4,JTvalue=59){ m<-m1+m2+m3;ex5.Jonckhexactch5=function(m){ a<-rep(0,5);for (i in 1:(m-4)){for (j in (i+1):(m-3)){ for (k in (j+1):(m-2)){for (l in (k+1):(m-1)){ for (f in (l+1):m){ind<-min((j-i),(k-i),(k-j),(l-i),(l-j),(l-k),(f-i),(f-j),(f-k),(f-l)) if (ind>0) a<-rbind(a,c(i,j,k,l,f))}}}}}a[2:nrow(a),]} JTid1<-ex5.Jonckhexactch5(m1+m2+m3);n1<-nrow(JTid1) JTid2<-ex5.Jonckhexactch5(m2+m3);n2<-nrow(JTid2) const<-1:m;JT<-matrix(0,nrow=n1*n2,ncol=m+1) for (i in 1:n1){for (j in 1:n2){temp1<-c(JTid1[i,]);JT[j+(i-1)*n2,1:m1]<-temp1 temp2<-(const[-temp1])[c(JTid2[j,])];JT[j+(i-1)*n2,(m1+1):(m1+m2)]<-temp2 temp3<-const[-JT[j+(i-1)*n2,1:(m2+m1)]] JT[j+(i-1)*n2,(m1+m2+1):(m1+m2+m3)]<-temp3 JT[j+(i-1)*n2,m1+m2+m3+1]<-sum(outer(temp2,temp1,">"))+ sum(outer(temp3,temp1,">"))+sum(outer(temp3,temp2,">"))}} y<-JT[,m+1];pvalue<-(sum(y>=JTvalue))/(nrow(JT));hist(y,breaks=min(y):max(y)) z<-c(0,hist(y,breaks=min(y):max(y))$counts) list(cbind("m1"=m1,"m2"=m2,"m3"=m3,"JTvalue"=JTvalue), cbind("n1"=n1,"n2"=n2,"n1 times n2"=nrow(JT)), c("min"=min(y),"max"=max(y),"JTvalue"=JTvalue,pvalue), cbind(min(y):max(y),z,cumsum(z),cumsum(z)/(n1*n2),1-cumsum(z)/(n1*n2)))
Histogram ofy000000008o0009000000000102030405060y
Histogram of y y Frequency 0 10 20 30 40 50 60 0 2000 4000 6000 8000 10000
m1m2m3JTvalue4559[1,] 5n1n2n1 times n2126[1,]2002252252minJTvaluep-valuemax065595.272505e-04J=Ui;的精确分布表(m1=5,m2=5,m3=4)i<j
m1 m2 m3 JTvalue [1,] 5 5 4 59 n1 n2 n1 times n2 [1,] 2002 126 252252 min max JTvalue p-value 0 65 59 5.272505e-04 的精确分布表(m1=5,m2=5,m3=4)