### Here we fit the MLE for the normal

# X_1, X_2, ...,X_n iid Normal(mu, sigma2), find the MMEs.
# p.53 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 <- ((n-1)/n)*x.sd

# 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="MME 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] <- ((n-1)/n)*sd(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")

mu.mle.boot.m <- mean(mu.hat.mle.star)
mu.mle.boot.s <- sd(mu.hat.mle.star)
mu.mle.boot.ci <- quantile(mu.hat.mle.star,c(0.025,0.975))
mu.mle.boot.ci

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