#HOTELLING'S T2 TEST FOR SINGLE/PAIRED POPULATIONS:

#LOAD DATA & HYPOTHESIS:
X=read.table("c:/DATA/Multivariate\\/T5-1.DAT")
X              #DATA TABLE

#=============================================
#FUNCTION FOR ONE SAMPLE HOTELLING'S T2 TEST:
#X     = dataset
#mu0   = hypothesis vector
#alpha = alpha of the test

OS.Hotelling.T2 <- function(X,mu0,alpha)
{

  n=nrow(X)
  p=ncol(X)
  Xbar=colMeans(X)
  S=cov(X)

  Tsq=n%*%((Xbar-mu0)%*%solve(S)%*%(Xbar-mu0))
  Tsq

  Sigmahat=((n-1)/n)*S
  XX=as.matrix(sweep(X,MARGIN=2,mu0,FUN="-"))
  Sigma0hat=t(XX)%*%XX
  Wilks.lambda=(det(n*(Sigmahat))/(det(Sigma0hat)))^(n/2)
  Wilks.lambda

  C=(((n-1)*p)/(n-p))*qf(1-alpha,p,n-p)
  C

  F=Tsq*(n-p)/((n-1)*p)
  df1=p
  df2=n-p

  P=1-pf(F,p,n-p)
  P

  lambda=eigen(S)$values
  eigenvec=eigen(S)$vectors
  FF=(p*(n-1))/(n*(n-p))
  CT=sqrt(FF*qf(1-alpha,p,n-p))
  CT
  ct=qt(1-alpha/2,n-1)
  ct
  cb=qt(1-alpha/(2*p),n-1)
  cb

  sigma=as.numeric(vector(length=p))
  T2.CI.lower=as.numeric(vector(length=p))
  T2.CI.upper=as.numeric(vector(length=p))
  Univar.CI.lower=as.numeric(vector(length=p))
  Univar.CI.upper=as.numeric(vector(length=p))
  Bonf.CI.lower=as.numeric(vector(length=p))
  Bonf.CI.upper=as.numeric(vector(length=p))
  L=vector(length=p)

  for (i in 1:p) 
  {
    L[i]=sqrt((1/n)*C*lambda[i])
    sigma[i]=S[i,i]

    T2.CI.lower[i]=Xbar[i]-CT*sqrt(sigma[i])
    T2.CI.upper[i]=Xbar[i]+CT*sqrt(sigma[i])

    Univar.CI.lower[i]=Xbar[i]-ct*sqrt(sigma[i]/n)
    Univar.CI.upper[i]=Xbar[i]+ct*sqrt(sigma[i]/n)

    Bonf.CI.lower[i]=Xbar[i]-cb*sqrt(sigma[i]/n)
    Bonf.CI.upper[i]=Xbar[i]+cb*sqrt(sigma[i]/n)

  }

  CI=cbind(T2.lower=T2.CI.lower,B.lower=Bonf.CI.lower,t.lower=Univar.CI.lower,Mean=Xbar,
         t.upper=Univar.CI.upper,B.upper=Bonf.CI.upper,T2.upper=T2.CI.upper)

  result=list(hypothesis.vector=mu0,HotellingT2=Tsq,WilksLambda=Wilks.lambda,CV=C,Prob=P,
           CI=CI,eigenvectors=eigenvec,eigenvalues=lambda,mean.vector=Xbar,CT=CT)
  class(result) = 'OS.Hotelling.T2'
  result
  cat("\n","One Sample Hotelling's T2","\n","\n",    
     "Hypothesis Vector:       (",mu0,")","\n",
     "T2 Ellipsoid half lengths: (",L,")","\n",
     "Hotelling's T2 Statistic: ",Tsq,"\n",
     "Equivalent F Statistic :  ",F,"\n",
     "F degrees of freedom:      (",df1,df2,")","\n",
     "Wilks's Lambda:           ",Wilks.lambda,"\n",
     "alpha:                    ",alpha,"\n",
     "Critical Value:           ",C,"\n",
     "Probability:              ",P,"\n","\n", 
  "Confidence Intervals: T2 - Bonferroni - t - Mean - t - Bonferroni - T2","\n","\n")
  print(CI)
  cat("\n")

  return(result)

}
#=============================================

RES=OS.Hotelling.T2(X,mu0=c(4,50,10),alpha=0.05)

#=============================================

library(ICSNP)
HotellingsT2(X,mu=c(4,50,10))


#=============================================
library(rrcov)
T2.test(X,mu=c(4,50,10),conf.level=0.95,method="c")


#=============================================

#FUNCTION PLOT CI:
#M = data
#R = return from function OS.Hotelling.T2 as a list
#X = variable to plot along the X axis
#Y = variable to plot along the Y axis

plot.CI <- function(M,R,X,Y)
{
#  lambda=R$eigenvalues
#  eigenvec=R$eigenvectors
  mu0=R$hypothesis.vector
  Xbar=R$mean.vector
#  CT=R$CT

#  L=CT*sqrt(lambda)
  

  plot(M[,X],M[,Y],type='n',xlab=X,ylab=Y)

  points(M[,X],M[,Y],pch=19,cex=0.5,col='blue')
  points(mu0[X],mu0[Y],pch=2,col='red')
  points(Xbar[X],Xbar[Y],pch=3,col='blue')

  abline(v=R$CI[X,1],col='red')
  abline(v=R$CI[X,7],col='red')
  abline(v=R$CI[X,2],col='brown',lty=3)
  abline(v=R$CI[X,6],col='brown',lty=3)
  abline(v=R$CI[X,3],col='brown',lty=2)
  abline(v=R$CI[X,5],col='brown',lty=2)

  abline(h=R$CI[Y,1],col='red')
  abline(h=R$CI[Y,7],col='red')
  abline(h=R$CI[Y,2],col='brown',lty=3)
  abline(h=R$CI[Y,6],col='brown',lty=3)
  abline(h=R$CI[Y,3],col='brown',lty=2)
  abline(h=R$CI[Y,5],col='brown',lty=2)

}

plot.CI(X,RES,1,2)

plot.CI(X,RES,1,3)

plot.CI(X,RES,2,3)