
#Computes the power of the Hotelling's T^2 test for the plots in Figure 5 and 8. 



library(mnormt)

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

#Computes the power of the Hotelling's T^2 test 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)


rpnull=vector(length=iter)
rp=vector(length=iter)


for(j in 1:20)
{
  
  
  delta1<-h[j]*rep(1, d)/sqrt(m+n)	
  mu0=rep(0, d)
  mu1=delta1
  sigma0=diag(rep(1, d))
  
  
  for(i in 1:iter)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu1, sigma0)
    xmean<-apply(points1, 2, mean)
    ymean<-apply(points2, 2, mean)
    diff<-xmean-ymean
    vc<-(1/m+1/n)*(((m-1)*cov(points1)+(n-1)*cov(points2))/(m+n-2))
    rp[i]<-t(diff)%*%solve(vc)%*%diff

  }
  
  power[j]<-length(which(rp> qchisq(0.95, d)))/iter

}

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



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

#Computes the power of the Hotelling's T^2 test for the plot in Figure 8(b) (normal scale problem in dimension d=10)


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

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

power=vector(length=20)


rpnull=vector(length=iter)
rp=vector(length=iter)


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)
    xmean<-apply(points1, 2, mean)
    ymean<-apply(points2, 2, mean)
    diff<-xmean-ymean
    vc<-(1/m+1/n)*(((m-1)*cov(points1)+(n-1)*cov(points2))/(m+n-2))
    rp[i]<-t(diff)%*%solve(vc)%*%diff
    
  }
  
  power[j]<-length(which(rp> qchisq(0.95, d)))/iter
  
}

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



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

#Computes the power of the Hotelling's T^2 test for the plot in Figure 5(a) (lognormal with independent covaraince in dimension d=10)


library(mnormt)

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

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

power=vector(length=20)


rpnull=vector(length=iter)
rp=vector(length=iter)


for(j in 1:20)
{
  
  
  delta1<-h[j]*rep(1, d)/sqrt(m+n)	
  mu0=rep(0, d)
  mu1=delta1
  sigma0=diag(rep(1, d))
  
  
  for(i in 1:iter)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu1, sigma0))
    xmean<-apply(points1, 2, mean)
    ymean<-apply(points2, 2, mean)
    diff<-xmean-ymean
    vc<-(1/m+1/n)*(((m-1)*cov(points1)+(n-1)*cov(points2))/(m+n-2))
    rp[i]<-t(diff)%*%solve(vc)%*%diff

  }
  
  power[j]<-length(which(rp> qchisq(0.95, d)))/iter

}

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



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


#Computes the power of the Hotelling's T^2 test for the plot in Figure 5(b) (lognormal with dependent covaraince in dimension d=10)



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

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

power=vector(length=20)


rpnull=vector(length=iter)
rp=vector(length=iter)

mu0=rep(0, d)
sigma0=rep(1, d)%*%t(rep(1, d))+diag(rep(1,d))
#sigma0=as.matrix(read.table(file="Lognormal_Dependent_Sigma10.txt"))


for(j in 1:20)
{
  
  
  delta1<-h[j]*c(rep(1, d/2), rep(1, d/2))/sqrt(m+n)	
  mu1=delta1
 
  
  
  for(i in 1:iter)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu1, sigma0))
    xmean<-apply(points1, 2, mean)
    ymean<-apply(points2, 2, mean)
    diff<-xmean-ymean
    vc<-(1/m+1/n)*(((m-1)*cov(points1)+(n-1)*cov(points2))/(m+n-2))
    rp[i]<-t(diff)%*%solve(vc)%*%diff
    
  }
  
  power[j]<-length(which(rp> qchisq(0.95, d)))/iter
  
}



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

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