### Here we fit the MLE for the normal

# X_1, X_2, ...,X_n iid Normal(mu, sigma2), find the MLEs.
# p.54 f(x|mu,sigma) = (1/(sigma*sqrt(2*pi)))*exp(-(x-mu)^2/2*sigma^2)

# Consider a random sample from the normal distribution

n <- 30
mu <- 70
sigma <- 3
sigma2 <- 3^2

x <- rnorm(n, mu, sigma)

x.orig <- x

x.mean <- mean(x)
x.sd <- sd(x)
x.sigma <- sqrt(((n-1)/n)*var(x))

# make a histogram

hist(x,probability=T)

x.min <- min(hist(x)$breaks)
x.max <- max(hist(x)$breaks)
y.max <- max(hist(x)$density) + 0.05

# use the density estimator since this is a continuous random variable

plot(density(x, n=50, window="g"),type="l",xlab="x",ylab="density")

# MLE fitted normal p.255

mu.hat.mle <- mean(x)
sigma.hat.mle <- sqrt((1/length(x))*sum((x-mean(x))^2))

# plot the fitted model to the simulated data

xx <- seq(x.min,x.max,0.1)

yy <- dnorm(xx, mu.hat.mle, sigma.hat.mle)

plot(xx,yy, main="MLE fitted normal",type="l")

# make plot of histogram with the fitted model

hist(x,probability=T,xlim=c(x.min,x.max),ylim=c(0,y.max),main="MLE fitted normal",
  xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(x.min,x.max),ylim=c(0,y.max),xlab="", ylab="")

# make a plot of the density estimator with the fitted model.

plot(density(x, n=50, window="g"),type="l",xlim=c(x.min,x.max),ylim=c(0,y.max),
   main="MLE fitted normal",xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(x.min,x.max),ylim=c(0,y.max),
   xlab="x",ylab="density")

# exact confidence interval for mu

mu.ci.exact <- c(x.mean - 1.96*x.sd/sqrt(n), x.mean + 1.96*x.sd/sqrt(n))
mu.ci.exact

# large sample confidence interval for mu

mu.ci.large <- c(mu.hat.mle - 1.96*sigma.hat.mle/sqrt(n), mu.hat.mle +
	1.96*sigma.hat.mle/sqrt(n))
mu.ci.large

# parametric bootstrap the MLE

B <- 1000

mu.hat.mle.star <- numeric(B)
sigma.hat.mle.star <- numeric(B)

for(i in 1:B){
	x <- rnorm(n,mu.hat.mle,sigma.hat.mle)
	mu.hat.mle.star[i] <- mean(x)
	sigma.hat.mle.star[i] <- sqrt(((n-1)/n)*var(x))
}

plot(density(mu.hat.mle.star, n=50, window="g"),type="l",xlab="x",ylab="density")
plot(density(sigma.hat.mle.star, n=50, window="g"),type="l",xlab="x",ylab="density")

# bootstap mu

mu.mle.boot.m <- mean(mu.hat.mle.star)
mu.mle.boot.s <- sd(mu.hat.mle.star)

alpha <- 0.05
alpha.b <- c((alpha/2),(1-alpha/2))

# percentile method

mu.mle.boot.ci <- quantile(mu.hat.mle.star,alpha.b)
mu.mle.boot.ci

# bootstrap confidence interval - presented in the Rice book

mu.hat.mle.delta <- mu.hat.mle.star - mu.hat.mle
delta.up <- quantile(mu.hat.mle.delta,(1-alpha/2))
delta.low <- quantile(mu.hat.mle.delta,(alpha/2))
boot.LB <- mu.hat.mle - delta.up
boot.LB
boot.UB <- mu.hat.mle - delta.low
boot.UB           


# bootstap sigma

sigma.mle.boot.m <- mean(sigma.hat.mle.star)
sigma.mle.boot.s <- sd(sigma.hat.mle.star)

# percentile method

sigma.mle.boot.ci <- quantile(sigma.hat.mle.star,c(0.025,0.975))
sigma.mle.boot.ci

# bootstrap confidence interval - presented in the Rice book

sigma.hat.mle.delta <- sigma.hat.mle.star - sigma.hat.mle
delta.up <- quantile(sigma.hat.mle.delta,(1-alpha/2))
delta.low <- quantile(sigma.hat.mle.delta,(alpha/2))
boot.LB <- sigma.hat.mle - delta.up
boot.LB
boot.UB <- sigma.hat.mle - delta.low
boot.UB