
#Computes the power of the CovTest for the plot in Figure 3 (normal scale problem in dimension d=10). 


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


d=10
m=1000
n=500
iter=500
iterations=100

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

power=vector(length=20)


rp1=vector(length=iter)
rp=vector(length=iterations)

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


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(j in 1:20)
{
  
  
  delta1<-h[j]/(m+n)^(1/2)  
  sigma1=sigma0+delta1*diag(rep(1, d))
  

  
  for(i in 1:iterations)
  {
    
    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)))/iterations
  #power2[j]<-length(which(rpn< qnorm(0.05)))/iterations
  
}


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




