
#Power of the Hoteling's T^2 Test for the tables in Figure 2, 3, and 4.


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

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

#Power of the Hotelling's T^2 Test for the table in Figure 2 (normal location).



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

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

for(j in 1:5)
{
  rp=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<-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[j])))/iter

}

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


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


#Power of the Hotelling's T^2 Test for the table in Figure 3 (normal scale).



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

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

for(j in 1:5)
{
  rp=rpnull=vector(length=iter)
  
  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, 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[j])))/iter

}


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


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

#Power of the Hotelling's T^2 Test for the table in Figure 4 (lognormal location).



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

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

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, 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[j])))/iter

}

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




