### HW 6 solutions
### AHW3
### December 2, 2007

rm(list = ls())
############################################
## 1
y <- as.matrix(read.table("bird.txt", header=T))
n <- nrow(y)
p <- ncol(y)

## null hypothesis
mu.0 <- matrix(c(200, 250),2,1)

## get the unconstrained mle for Sigma
Sigma.mle.uncon <- ((n-1)/n)*cov(y)
Sigma.mle.uncon

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

## likelihood ratio test statistic
lambda.obs <- (det(Sigma.mle.uncon)/det(Sigma.mle.con))^(n/2)
L.test.obs <- -2*log(lambda.obs)
L.test.obs

## critical value
alpha <- 0.05
crit.val <- qchisq(1-alpha, p)
crit.val

## decision
L.test.obs > crit.val


##############################################
## 2
y.g <- as.matrix(read.table("milkGas.txt", header=T))
y.d <- as.matrix(read.table("milkDiesel.txt", header=T))

n.g <- nrow(y.g)
n.d <- nrow(y.d)
p <- ncol(y.g)


# sample means also mles
y.bar.g <- apply(y.g, 2, mean)
y.bar.d <- apply(y.d, 2, mean)

# sample covariances
S.g <- cov(y.g)
S.d <- cov(y.d)

# pooled covariance
S.pooled <- ((n.g-1)*S.g+(n.d-1)*S.d)/(n.g+n.d-2)

## calculate mles from the sample estimators
S.g.mle <- ((n.g-1)/n.g)*S.g
S.g.mle

S.d.mle <- ((n.d-1)/n.d)*S.d
S.d.mle

S.pooled.mle <- ((n.g+n.d-2)/(n.g+n.d))*S.pooled
S.pooled.mle

# conduct the likelhood ratio test
L.test.obs <- -2*((n.g/2)*log(det(S.g.mle)) + (n.d/2)*log(det(S.d.mle)) - ((n.g+n.d)/2)*log(det(S.pooled.mle)) )
L.test.obs

## critical value
alpha <- 0.05
crit.val <- qchisq(1-alpha, p*(p+1)/2)
crit.val

## decision
L.test.obs > crit.val

###################################################
## 3
FT <- read.table("turtlesFemale.txt", header=T)
MT <- read.table("turtlesMale.txt", header=T)

p <- ncol(FT)
n.f <- nrow(FT)
n.m <- nrow(MT)

# sample means 
y.bar.f <-apply(FT, 2, mean)
y.bar.m <- apply(MT, 2, mean)

# sample covariances
S.f <- cov(FT)
S.f

S.m <- cov(MT)
S.m

# pooled covariance
S.pooled <- ((n.f-1)*S.f+(n.m-1)*S.m)/(n.f+n.m-2)
S.pooled

##
C <- matrix(c(1,1,1),1,3)
C

## Hotellings T^2
T.sq.obs <- t(y.bar.f-y.bar.m)%*%t(C)%*%solve( ((1/n.f)+(1/n.m))*C%*%S.pooled%*%t(C))%*%C%*%(y.bar.f-y.bar.m)
T.sq.obs

## critical value
alpha <- 0.05
crit.val <- qf(1-alpha,1, n.f+n.m-2)
crit.val

## decision
T.sq.obs > crit.val

###################################################
## 4 
y.t <- read.table("flyCarteri.txt", header=T)
y.c <- read.table("flyTorrens.txt", header=T)

p <- ncol(y.t)
n.t <- nrow(y.t)
n.c <- nrow(y.c)

# sample means 
y.bar.t <-apply(y.t, 2, mean)
y.bar.t

y.bar.c <- apply(y.c, 2, mean)
y.bar.c

# sample covariances
S.t <- cov(y.t)
S.t

S.c <- cov(y.c)
S.c

####(a)

## large sample statistic
LS.obs <- t(y.bar.t-y.bar.c)%*%solve( (1/n.t)*S.t +(1/n.c)*S.c)%*%(y.bar.t-y.bar.c)
LS.obs

## critical value
alpha <- 0.05
crit.val <- qchisq(1-alpha,p)
crit.val

## decision
LS.obs > crit.val


######## (b)
LS.int <- NULL
for(k in 1:p){
x.k <- sqrt( ((1/n.t)*S.t +(1/n.c)*S.c)[k,k])*sqrt(qchisq(1-alpha,p))
LS.k <- c(y.bar.t[k]-y.bar.c[k] - x.k, y.bar.t[k]-y.bar.c[k] + x.k)
LS.k <- round(LS.k,3)
LS.int <- rbind(LS.int, LS.k)
}
colnames(LS.int) <- c("lower", "upper")
rownames(LS.int) <- names(y.t)
LS.int