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



library(cluster)
library(mnormt)
library(localdepth)
library(fda.usc)


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

#Power of the test based on the halfspace depth for the table in Figure 2 (normal location).


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

m=1000
n=500
iter=100
h=2
power=vector(length=5)

for(j in 1:5)
{
  rp=rpnull=vector(length=iter)
  
  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:iter)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu1, sigma0)
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    t<-sum(outer(a,b,"<")*1)/(m*n)
    rp[i]<-(t-0.5)/(sqrt((1/m+1/n)/12))
    
  }
  power[j]<-mean((abs(rp)> qnorm(0.975))*1)
  
}

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

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

#Power of the test based on the halfspace depth for the table in Figure 3 (normal scale).


library(cluster)
library(mnormt)
library(localdepth)
library(fda.usc)

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

m=1000
n=500
iter=100
h=2
power=vector(length=5)

for(j in 1:5)
{
  rp=rpnull=vector(length=iter)
  
  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:iter)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu0, sigma1)
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    t<-sum(outer(a,b,"<")*1)/(m*n)
    rp[i]<-(t-0.5)/(sqrt((1/m+1/n)/12))
    
  }
  power[j]<-mean((abs(rp)> qnorm(0.975))*1)
  
}

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

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

#Power of the test based on the Mahalanobis depth for the table in Figure 3 (lognormal location).
#Can be changed to obtain the power of the test based on halfspace depth by changing mdepth.MhD to mdepth.TD below.


library(cluster)
library(mnormt)
library(localdepth)
library(fda.usc)

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

m=1000
n=500
iter=100
h=2
power=vector(length=5)

for(j in 1:5)
{
  rp=rpnull=vector(length=iter)
  
  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:iter)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu1, sigma0))
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    t<-sum(outer(a,b,"<")*1)/(m*n)
    rp[i]<-(t-0.5)/(sqrt((1/m+1/n)/12))
    
  }
  power[j]<-mean((abs(rp)> qnorm(0.975))*1)
  
}

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



