

#Power of the Friedman-Rafsky Test for the tables in Figure 6 and 7.

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

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

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


d=c(10, 30, 50, 70, 100)

m=60
n=40
iter=500
iterations=100
h=2

power=vector(length=5)

for(j in 1:5)
{
  
rpnull=vector(length=iter)
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]))

#Calculating the null distribution

for(i in 1:iter)
{
  
  points1<-rmnorm(m, mu0, sigma0)
  points2<-rmnorm(n, mu0, 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))
  
  rpnull[i]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
  
}


#Calculating the power of the Friedman-Rafsky Test in dimension 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]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))


}

power[j]<-length(which(rp< quantile(rpnull, 0.05)))/iterations

}

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

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

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


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


d=c(10, 30, 50, 70, 100) 

m=60
n=40
iter=500
iterations=100
h=2


power=vector(length=5)

for(j in 1:5)
{
 
  
  rpnull=vector(length=iter)
  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]))
  
  
  #Calculating the null distribution
  
  for(i in 1:iter)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu0, 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))
    
    rpnull[i]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
    
  }
  
  
  #Calculating the power of the Friedman-Rafsky Test in dimension 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]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
    
    
  }
  
  power[j]<-length(which(rp< quantile(rpnull, 0.05)))/iterations
  
}

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


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


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

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


d=c(10, 30, 50, 70, 100)
m=60
n=40
iter=500
iterations=100
h=2
power=vector(length=5)

for(j in 1:5)
{
  
  rpnull=vector(length=iter)
  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]))
  
  #Calculating the null distribution
  
  for(i in 1:iter)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu0, 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))
    
    rpnull[i]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
    
  }
  
  
  #Calculating the power of the Friedman-Rafsky Test in dimension 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]<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
    
    
  }
  
  power[j]<-length(which(rp< quantile(rpnull, 0.05)))/iterations
  
}

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



