## Chapter 3
##
## Code to support Video Vignette
##
## Videos and supporting code are not a complete portrayal
## of Chapter content.

## heights of women from class
y <- c(65, 65, 69, 67, 66, 66, 67, 64, 68, 68, 68, 69, 70, 65, 70, 
       67, 65, 62, 69, 68, 65, 68, 64, 64)
n <- length(y)

## what google says is typical of American women
mu <- 64.5          ## mu_0 from the notes
sigma2 <- 2.5^2     ## assume known

##
## a: Likelihood
##

## MLE mu
muhat = mean(y)
muhat

## MLE sigma2
sum((y - muhat)^2)/n
sum((y - muhat)^2)/(n-1)
s2 <- var(y)
s2

##
## b: Inference
##

## z-test samping distribution
ygrid <- seq(61.5, 67.5, length=1000)
plot(ygrid, dnorm(ygrid, mu, sqrt(sigma2/n)), type="l", 
     xlab="ybar = muhat", ylab="density N(64.5, 2.5^2/n)")
abline(v=muhat, col=2, lwd=2)
legend("topleft", "observed", lty=1, col=2, lwd=2, bty="n")

## z-test p-value
pval <- 2*integrate(dnorm, muhat, Inf, mean=mu, sd=sqrt(sigma2/n))$value
pval
2*pnorm(muhat, mu, sqrt(sigma2/n), lower.tail=FALSE)

## standardizing
se <- sqrt(sigma2/n)
z <- (muhat - mu)/se
z

## p-value from standardization
2*pnorm(-abs(z))

## t-test

## MC
N <- 100000
Ss <- rep(NA, N)                   
for(i in 1:N) {  
  sigma2s <- (n - 1)*s2/rchisq(1, n - 1)
  Ys <- rnorm(n, mu, sqrt(sigma2s))
  Ss[i] <- mean(Ys)                     
}

## p-value
pval <- 2*mean(Ss > muhat)
pval

## closed form via standarization
se <- sqrt(s2/n)           ## same as z, but with sigma2 = s2
t <- (muhat - mu)/se       ## same as z, different se above
t

## p-value
pval.exact <- 2*pt(-abs(t), df=n - 1)
pval.exact

## comparison
c(mc=pval, exact=pval.exact, lib=t.test(y, mu=mu)$p.value)

## 
## c: Confidence Intervals
##

## level = 1 - alpha
alpha <- 0.05

## first via MC
muhats <- rep(NA, N)                   
for(i in 1:N) {  
  sigma2s <- (n - 1)*s2/rchisq(1, n - 1)
  Ys <- rnorm(n, muhat, sqrt(sigma2s))        ## use muhat instead of mu
  muhats[i] <- mean(Ys)                     
}
CImc.mu <- quantile(muhats, c(alpha/2, 1 - alpha/2))
CImc.mu

## then analytically
q <- qt(c(alpha/2, 1 - alpha/2), n-1)
CI.mu <- muhat + q*se
CI.mu