
#Computes the power of the test based on halfspace-depth for the plots in Figure 5 and 8.
#Can be changed to obtain the power of the test based on Mahalonobis depth by changing mdepth.TD to mdepth.MhD in the code below. 



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

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

#Computes the power of the test based on halfspace depth for the plot in Figure 8(a) (normal location problem in dimension d=10)

d=10
m=1000
n=800
iter=100

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

power=vector(length=20)


rp=vector(length=iter)

mu0=rep(0, d)
sigma0=diag(rep(1, d))



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



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

#Computes the power of the test based on halfspace depth for the plot in Figure 8(b) (normal scale problem in dimension d=10)


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


d=10
m=1000
n=800
iter=100

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

power=vector(length=20)


rp=vector(length=iter)

mu0=rep(0, d)
sigma0=diag(rep(1, d))



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



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

#Computes the power of the test based on halfspace depth for the plot in Figure 5(a) (lognormal normal with independent coordinates in dimension d=10)


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


d=10
m=1000
n=800
iter=100

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

power=vector(length=20)


rp=vector(length=iter)

mu0=rep(0, d)
sigma0=diag(rep(1, d))



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



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

#Computes the power of the test based on halfspace depth for the plot in Figure 5(b) (lognormal normal with dependent coordinates in dimension d=10)


d=10
m=1000
n=800
iter=100

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

power=vector(length=20)


rp=vector(length=iter)

mu0=rep(0, d)
sigma0=rep(1, d)%*%t(rep(1, d))+diag(rep(1,d)) 



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


