## October 30 Sweat Data
rm(list = ls())

## for the example we are doing an alpha=0.10 test

###################################################################################
## example from the book (Sweat data)
y.bar <- matrix(c(4.640, 45.400, 9.965), 3, 1)
S <- matrix(c(2.879,10.010,-1.810,10.010,199.788,-5.640,-1.810,-5.640,3.628),3,3)
n <- 20
p <- 3

## hypothesized mean value
mu.0 <- matrix(c(4,50,10),3,1)

## calculate the test statistic
T.sq.obs <- n*t(y.bar-mu.0)%*%solve(S)%*%(y.bar-mu.0)

## calculate the distribution under the null hypothesis
T.sq.null <- (((n-1)*p)/(n-p))*qf(0.90, df1=3, df2=n-p)


## So, dow we reject the null hypothesis?
T.sq.obs
T.sq.null

## let's see this visually
x <- seq(0,10, by=0.1)
plot(x, (((n-1)*p)/(n-p))*df(x, df1=3, df2=n-p), type="l", lwd=2)
abline(v=T.sq.null, col="red", lwd=2)
abline(v=T.sq.obs, col="green", lwd=2)



#############################################################################
## continue Sweat data example
## use Likelihood ratio test

Lambda.2.n.obs <- (1 + T.sq.obs/(n-1))^(-1)
Lambda.2.n.null <- (1 + T.sq.null/(n-1))^(-1)

Lambda.2.n.obs
Lambda.2.n.null

#############################################################################
## continue Sweat data example
## now let's look at the assymptotic test.  Here n is really to small for the test but let's procede.

Lambda.obs <- -2*log((1 + T.sq.obs/(n-1))^(-n/2))
Lambda.null <- qchisq(0.90,  df=3)

Lambda.obs
Lambda.null

## let's see this visually
x <- seq(0,10, by=0.1)
plot(x, dchisq(x, df=3), type="l", lwd=2)
abline(v=Lambda.null, col="red", lwd=2)
abline(v=Lambda.obs, col="green", lwd=2)

##############################################################################
## Do the likelihood ratio test again from first principles

## get mle's for unconstrained model
y <- as.matrix(read.table("sweat.txt", header=T))
mu.mle <- apply(y, 2, mean)
Sigma.mle <- ((n-1)/n)*cov(y)

## get mle for Sigma for the constrained model
Sum.Sq <- matrix(0,3,3)
for(i in 1:n){
Sum.Sq <- Sum.Sq + (y[i,]-mu.0)%*%t(y[i,]-mu.0)
}
Sigma.mle.null <- Sum.Sq/n

Lambda.obs.2 <- -2*log((det(Sigma.mle)/det(Sigma.mle.null))^(n/2))

Lambda.obs.2
Lambda.null