############### Earwig Example ####################

gen<-matrix(0,8,8)
dimnames(gen)[[2]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
dimnames(gen)[[1]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
gen[2,1]<-.3
gen[3,1:2]<-c(.14,.5)
gen[4,1:3]<-c(.23,.5,.54)
gen[5,1:4]<-c(.3,.4,.5,.61)
gen[6,1:5]<-c(-.04,.04,.11,.03,.15)
gen[7,1:6]<-c(.02,.09,.14,-.16,.11,.14)
gen[8,1:7]<-c(-.09,-.06,.05,-.16,.03,-.06,.36)
gen<-diag(1,8)+gen+t(gen)

dis1<-matrix(0,8,8)
dimnames(dis1)[[2]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
dimnames(dis1)[[1]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
dis1[2,1]<-1
dis1[3,1:2]<-c(2,1)
dis1[4,1:3]<-c(1,2,3)
dis1[5,1:4]<-c(2,3,4,1)
dis1[6,1:5]<-c(3,4,5,2,1)
dis1[7,1:6]<-c(2,3,4,3,4,5)
dis1[8,1:7]<-c(1,2,3,2,3,4,1)
dis1<-diag(0,8)+dis1+t(dis1)

dis2<-matrix(0,8,8)
dimnames(dis2)[[2]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
dimnames(dis2)[[1]]<-c("EA","AF","Ma","OR","AU","NE","SA","NA")
dis2[2,1]<-1
dis2[3,1:2]<-c(2,1)
dis2[4,1:3]<-c(1,1,1)
dis2[5,1:4]<-c(2,1,1,1)
dis2[6,1:5]<-c(3,2,2,2,1)
dis2[7,1:6]<-c(2,1,2,2,2,3)
dis2[8,1:7]<-c(1,2,3,2,3,4,1)
dis2<-diag(0,8)+dis2+t(dis2)

iter <- 10000
Z1   <-  Z2 <- rep(0,iter)

for (i in 1:iter){
  set.seed(i)
  tmp1 <- sample(1:8,size=8,replace=F)
  tmp2 <- gen[tmp1,tmp1]
  Z1[i] <- sum(tmp2[lower.tri(tmp2)]*dis1[lower.tri(dis1)]) 
  Z2[i] <- sum(tmp2[lower.tri(tmp2)]*dis2[lower.tri(dis2)]) 
}

Z1.obs <- sum(gen[lower.tri(gen)]*dis1[lower.tri(dis1)])
Z2.obs <- sum(gen[lower.tri(gen)]*dis2[lower.tri(dis2)])

pval1<-(length(which(Z1<=Z1.obs))/iter) #one-sided p-value 
pval2<-(length(which(Z2<=Z2.obs))/iter) #one-sided p-value 
c(pval1,pval2)

############### Swedish pine saplings example ####################

library(spatstat)
data(swedishpines) #original data is in the unit 1/10 meter

#Plot
plot(swedishpines,chars=19,main="")

#Mead's test
cnt <- c(6,2,4,6,5,4,4,6,3,4,2,3,5,6,4,7)
mead <- function(x){
  t.bar <- mean(x)
  tmp1  <- 16*(t.bar^2)
  tss   <- sum(x^2)-tmp1
  bss   <- 0.25*(sum(x[1:4])^2+sum(x[5:8])^2+
           sum(x[9:12])^2+sum(x[13:16])^2)-tmp1
  q     <- bss/tss
  q
}

q.obs <- mead(cnt)
q.obs #observed value

iter <- 5000
q.rt <- rep(0,iter)
for (i in 1:iter){
  set.seed(i)
  cnt.rand <- sample(cnt,size=16,replace=F)
  q.rt[i]  <- mead(cnt.rand)
}

length(which(q.rt>=q.obs))/iter #p-value

#Test based on distance

dis <- matrix(0,71,71)
x <- swedishpines$x/10
y <- swedishpines$y/10
for (i in 2:71){
  for (j in 1:(i-1)){
    dis[i,j]<-sqrt(((x[i]-x[j])^2+(y[i]-y[j])^2))
    dis[j,i]<-dis[i,j]
  }
}

FUN1<-function(dis,k){
  FUN <- function(x){
    sort(x)[(k+1)]
  }
  apply(dis,1,FUN)
}

mean(FUN1(dis,1)) #q1

q.obs <- rep(0,10)
for (i in 1:10){
  q.obs[i] <- mean(FUN1(dis,i))
}

iter   <- 500
q.rand <- matrix(0,iter,10)

for (it in 1:iter){
  set.seed(it)
  x.r   <- runif(1:71,min=0,max=10)
  y.r   <- runif(1:71,min=0,max=10)
  dis.r <- matrix(0,71,71)
  for (i in 2:71){
    for (j in 1:(i-1)){
      dis.r[i,j] <- sqrt(((x.r[i]-x.r[j])^2+(y.r[i]-y.r[j])^2))
      dis.r[j,i] <- dis.r[i,j]
    }
  }
  for (i in 1:10){
    q.rand[it,i] <- mean(FUN1(dis.r,i))
  }
}

FUN2 <- function(x){
  as.vector(quantile(x,prob=c(0.01,0.99)))[1:2]
}

ci.rand <- apply(q.rand,2,FUN2) #CI band from null