

#Computes the power of the Chen-Friedman test for the plots in Figure 2, 3, and 4. 


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

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


#Computes the power of the Chen-Friedman test for the plot in Figure 2 (normal location problem in dimension d=10)


d=10
m=1000
n=500
iterations=100

h=seq(0, 3, length.out=20)

power=vector(length=20)

rp=vector(length=iterations)


for(j in 1:20)
{
  
  
  delta1<-h[j]*rep(1, d)/sqrt(m+n)	
  mu0=rep(0, d)
  mu1=delta1
  sigma0=diag(rep(1, d))
  
  
  
  
  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="FRNewNormal.txt", ncolumns=1)

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


#Computes the power of the Chen-Friedman test for the plot in Figure 3 (normal scale problem in dimension d=10)



d=10
m=1000
n=500
iterations=100

h=seq(0, 3, length.out=20)

power=vector(length=20)

rp=vector(length=iterations)


for(j in 1:20)
{
  
  
  mu0=rep(0, d)
  delta1<-h[j]/(m+n)^(1/2)  
  sigma0=diag(rep(1, d))
  sigma1=sigma0+delta1*diag(rep(1, d))
  
  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="FRNewScale.txt", ncolumns=1)

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

#Computes the power of the Chen-Friedman test for the plot in Figure 4 (lognormal location problem in dimension d=10)


d=10
m=1000
n=500
iter=500
iterations=100

h=seq(0, 3, length.out=20)

power=vector(length=20)


rpnull=vector(length=iter)
rp=vector(length=iterations)


for(j in 1:20)
{

  
  delta1<-h[j]*rep(1, d)/sqrt(m+n)	
  mu0=rep(0, d)
  mu1=delta1
  sigma0=diag(rep(1, d))
  
  
  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="FRNewLognormal.txt", ncolumns=1)

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

