

# Code that generates the p-values in Table 1 (analysis of the sensorless drive diagonistic dataset).


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


dd<-read.table("SDD.txt", sep=" ")

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

#The test based on the halfspace depth and MST between the (1, 2)-th pair. 


a<-dd[dd[, 49]==1,][-49]
b<-dd[dd[, 49]==2,][-49]

pointsn<-as.matrix(rbind(a, b))


d<-c(1,2,3,48)

r1<-vector(length=4)

#The test based on the halfspace depth. 


for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-pointsn%*%e
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  a1<-mdepth.TD(q1, q1, scale=TRUE)$dep 
  b1<-mdepth.TD(q2, q1, scale=TRUE)$dep 
  t<-sum(outer(a1,b1,"<")*1)/(m*n)
  r1[j]<-abs(t-0.5)/(sqrt((1/m+1/n)/12))
}

halfspace<-1-pnorm(r1)


#The test based on the MST. 

r1<-vector(length=4)

VarD<-c(0, 0.453, 0.648, 1.192) #These are values from Table 1 of Henze-Penrose (1999), used for computing the null variance of the FR test statistic. 

for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-as.matrix(pointsn%*%e)
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  s<-2*(m*n)/(m+n)^2
  D<-s*(s+0.5*VarD[j]*(1-2*s))
  points<-as.matrix(rbind(q1, q2))
  
  x<-as.matrix(dist(points, upper = TRUE, method="euclidean"))
  child<-spantree(x)$kid
  edges<-t(rbind(c(2:(m+n)), child))
  
  r<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
  r1[j]<-(r-s*(m+n-1))/sqrt(D*(m+n))
  print(j)
}

mst<-pnorm(r1)

pvalues12<-rbind(halfspace, mst) 

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


#The test based on the halfspace depth and MST between the (1, 6)-th pair. 


a<-dd[dd[, 49]==1,][-49]
b<-dd[dd[, 49]==6,][-49]

pointsn<-as.matrix(rbind(a, b))


d<-c(1,2,3,48)

r1<-vector(length=4)

#The test based on the halfspace depth. 


for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-pointsn%*%e
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  a1<-mdepth.TD(q1, q1, scale=TRUE)$dep 
  b1<-mdepth.TD(q2, q1, scale=TRUE)$dep 
  t<-sum(outer(a1,b1,"<")*1)/(m*n)
  r1[j]<-abs(t-0.5)/(sqrt((1/m+1/n)/12))
}

halfspace<-1-pnorm(r1)


#The test based on the MST. 

r1<-vector(length=4)

for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-as.matrix(pointsn%*%e)
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  s<-2*(m*n)/(m+n)^2
  D<-s*(s+0.5*VarD[j]*(1-2*s))
  points<-as.matrix(rbind(q1, q2))
  
  x<-as.matrix(dist(points, upper = TRUE, method="euclidean"))
  child<-spantree(x)$kid
  edges<-t(rbind(c(2:(m+n)), child))
  
  r<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
  r1[j]<-(r-s*(m+n-1))/sqrt(D*(m+n))
  print(j)
}

mst<-pnorm(r1)

pvalues16<-rbind(halfspace, mst) 

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


#The test based on the halfspace depth and MST between the (2,6)-th pair. 


a<-dd[dd[, 49]==2,][-49]
b<-dd[dd[, 49]==6,][-49]

pointsn<-as.matrix(rbind(a, b))


d<-c(1,2,3,48)

r1<-vector(length=4)

#The test based on the halfspace depth. 


for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-pointsn%*%e
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  a1<-mdepth.TD(q1, q1, scale=TRUE)$dep 
  b1<-mdepth.TD(q2, q1, scale=TRUE)$dep 
  t<-sum(outer(a1,b1,"<")*1)/(m*n)
  r1[j]<-abs(t-0.5)/(sqrt((1/m+1/n)/12))
}

halfspace<-1-pnorm(r1)


#The test based on the MST. 

r1<-vector(length=4)

for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-as.matrix(pointsn%*%e)
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  s<-2*(m*n)/(m+n)^2
  D<-s*(s+0.5*VarD[j]*(1-2*s))
  points<-as.matrix(rbind(q1, q2))
  
  x<-as.matrix(dist(points, upper = TRUE, method="euclidean"))
  child<-spantree(x)$kid
  edges<-t(rbind(c(2:(m+n)), child))
  
  r<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
  r1[j]<-(r-s*(m+n-1))/sqrt(D*(m+n))
  print(j)
}

mst<-pnorm(r1)

pvalues26<-rbind(halfspace, mst) 

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


#The test based on the halfspace depth and MST between the (4, 9)-th pair. 


a<-dd[dd[, 49]==4,][-49]
b<-dd[dd[, 49]==9,][-49]

pointsn<-as.matrix(rbind(a, b))


d<-c(1,2,3,48)

r1<-vector(length=4)

#The test based on the halfspace depth. 


for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-pointsn%*%e
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  a1<-mdepth.TD(q1, q1, scale=TRUE)$dep 
  b1<-mdepth.TD(q2, q1, scale=TRUE)$dep 
  t<-sum(outer(a1,b1,"<")*1)/(m*n)
  r1[j]<-abs(t-0.5)/(sqrt((1/m+1/n)/12))
}

halfspace<-1-pnorm(r1)


#The test based on the MST. 

r1<-vector(length=4)

for(j in 1:4)
{
  e<-as.matrix(princomp(pointsn)$loadings[, c(1:d[j])])
  q<-as.matrix(pointsn%*%e)
  q1<-as.matrix(q[c(1:5319),])
  q2<-as.matrix(q[c(5320:10638),])
  
  m<-dim(q1)[1]
  n<-dim(q2)[1]
  s<-2*(m*n)/(m+n)^2
  D<-s*(s+0.5*VarD[j]*(1-2*s))
  points<-as.matrix(rbind(q1, q2))
  
  x<-as.matrix(dist(points, upper = TRUE, method="euclidean"))
  child<-spantree(x)$kid
  edges<-t(rbind(c(2:(m+n)), child))
  
  r<-length(which(edges[1:(m-1),][,2]>m))+length(which(edges[m:(m+n-1),][,2]<=m))
  r1[j]<-(r-s*(m+n-1))/sqrt(D*(m+n))
  print(j)
}

mst<-pnorm(r1)

pvalues49<-rbind(halfspace, mst) 

