


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


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

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

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

m=60
n=40
iter=500
h=2
power=vector(length=5)
rp=rp1=vector(length=iter)

for(j in 1:5)
{
  
  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, mu0, sigma0)
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    rp1[i]<-sum(outer(a,b,"<")*1)/(m*n)

  }
  
  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 
    
    rp[i]<-sum(outer(a,b,"<")*1)/(m*n)

  }

  power[j]<-(length(which(rp> quantile(rp1, 0.975)))+length(which(rp < quantile(rp1, 0.025))))/iter
  print(power[j])
}


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

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

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

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

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

m=60
n=40
iter=500
h=2
power=vector(length=5)
rp=rp1=vector(length=iter)

for(j in 1:5)
{
  
  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, sigma0)
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    rp1[i]<-sum(outer(a,b,"<")*1)/(m*n)
    
  }
  
  
  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 
    

    rp[i]<-sum(outer(a,b,"<")*1)/(m*n)
    
  }

  power[j]<-(length(which(rp> quantile(rp1, 0.975)))+length(which(rp < quantile(rp1, 0.025))))/iter
  print(power[j])
  
}

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

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

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


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

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

m=60
n=40
iter=500
h=2
power=vector(length=5)
rp=rp1=vector(length=iter)


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, mu0, sigma0))
    
    a<-mdepth.TD(points1, points1, scale=TRUE)$dep 
    b<-mdepth.TD(points2, points1, scale=TRUE)$dep 
    
    rp1[i]<-sum(outer(a,b,"<")*1)/(m*n)
    
  }
  
  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 
    
    rp[i]<-sum(outer(a,b,"<")*1)/(m*n)
    
  }

  power[j]<-(length(which(rp> quantile(rp1, 0.975)))+length(which(rp < quantile(rp1, 0.025))))/iter
  print(power[j])
  
}

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




