function (KW,nvec=c(8,10,15),nsim=1000) 
{
# This function simulates the p-value of an observed
# Kruskal-Wallis statistic KW, computed from samples
# of sizes nvec (default c(8,10,15)). The p-value is 
# based on nsim simulations.
#---------------------------------------------------
N<-sum(nvec)
k <- length(nvec)
nvec2<-cumsum(nvec)
nvec1<-c(0,nvec2[1:(k-1)])+1
out<- NULL
x <-list()
for(i in 1:nsim){
xx <- sample(1:N,replace=F)
for(j in 1:k){
x[[j]]<-xx[nvec1[j]:nvec2[j]]
}
out[i]<-kruskal.test(x)$statistic
}
y<-mean(out>=KW)
names(y)<-"p-value"
y
}