
#Power of the Chen-Friedman Test for the tables in Figure 2, 3, and 4.


library(cluster)
library(vegan)
library(mnormt)


########################################################

#Power of the Chen-Friedman Test for the table in Figure 2 (normal location).


d=c(4, 10, 20, 30, 50)

m=1000
n=500
iterations=100
h=2

power=vector(length=5)

for(j in 1:5)
{
  
  rp=vector(length=iterations)
  
  delta1<-h*rep(1, d[j])/sqrt(m+n)  	
  mu0=rep(0, d[j])
  mu1=delta1
  sigma0=diag(rep(1, d[j]))
  

  for(i in 1:iterations)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu1, sigma0)
    points<-rbind(points1,points2)
    x<-as.matrix(vegdist(points, binary=FALSE, method="euclidean"))
    child<-spantree(x)$kid
    edges<-t(rbind(c(2:(m+n)), child))
    
    rp[i]<-g.tests(edges, c(1:m), c((m+1):(m+n)), test.type="g", perm=0)$generalized$test.statistic
      
  }
  
  power[j]<-length(which(rp> qchisq(0.95, 2)))/iterations
  
}

write(as.vector(power), file="FR_NewNormal.txt", ncolumns=1)


##################################################

#Power of the Chen-Friedman Test for the table in Figure 3 (normal scale).


library(cluster)
library(vegan)
library(mnormt)


d=c(4, 10, 20, 30, 50)

m=1000
n=500
iterations=100
h=2

power=vector(length=5)

for(j in 1:5)
{
  
  rp=vector(length=iterations)
  
  mu0=rep(0, d[j])
  delta1<-h/(m+n)^(1/2)  
  sigma0=diag(rep(1, d[j]))
  sigma1=sigma0+delta1*diag(rep(1, d[j]))
  
  
  for(i in 1:iterations)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu0, sigma1)
    points<-rbind(points1,points2)
    x<-as.matrix(vegdist(points, binary=FALSE, method="euclidean"))
    child<-spantree(x)$kid
    edges<-t(rbind(c(2:(m+n)), child))
    
    rp[i]<-g.tests(edges, c(1:m), c((m+1):(m+n)), test.type="g", perm=0)$generalized$test.statistic
    
  }
  
  power[j]<-length(which(rp> qchisq(0.95, 2)))/iterations
  
}

write(as.vector(power), file="FR_NewScale.txt", ncolumns=1)


##################################################

#Power of the Chen-Friedman Test for the table in Figure 4 (lognormal location).


library(cluster)
library(vegan)
library(mnormt)


d=c(4, 10, 20, 30, 50)

m=1000
n=500
iterations=100
h=2

power=vector(length=5)

for(j in 1:5)
{
  
  rp=vector(length=iterations)
  
  delta1<-h*rep(1, d[j])/sqrt(m+n)  	
  mu0=rep(0, d[j])
  mu1=delta1
  sigma0=diag(rep(1, d[j]))
  
  
  for(i in 1:iterations)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu1, sigma0))
    points<-rbind(points1,points2)
    x<-as.matrix(vegdist(points, binary=FALSE, method="euclidean"))
    child<-spantree(x)$kid
    edges<-t(rbind(c(2:(m+n)), child))
    
    rp[i]<-g.tests(edges, c(1:m), c((m+1):(m+n)), test.type="g", perm=0)$generalized$test.statistic
    
  }
  
  power[j]<-length(which(rp> qchisq(0.95, 2)))/iterations
  
}

write(as.vector(power), file="FR_NewLogNormal.txt", ncolumns=1)


