
#Power of the CovTest for the table in Figure 3.


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



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

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





for(j in 1:5)
{
  
  
  
  delta1<-h/(m+n)^(1/2)	
  mu0=rep(0, d[j])
  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)
    points<-rbind(points1, points2)
    rp1[i]<-(m+n)*log(det((t(points)%*%points)/(m+n)))-m*log(det((t(points1)%*%points1)/m))-n*log(det((t(points2)%*%points2)/n))
    
  }
  
  
  for(i in 1:iter)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu0, sigma1)
    points<-rbind(points1, points2)
    rp[i]<-(m+n)*log(det((t(points)%*%points)/(m+n)))-m*log(det((t(points1)%*%points1)/m))-n*log(det((t(points2)%*%points2)/n))
    
    #rpn[i]<-(r-s*(m+n-1))/sqrt(D*(m+n))
    
  }
  
  power[j]<-length(which(rp> quantile(rp1, 0.95)))/iter
  #power2[j]<-length(which(rpn< qnorm(0.05)))/iterations
  
}


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




