### 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

X11()
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

X11()
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)

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

# make plot of histogram with the fitted model

X11()
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.

X11()
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)
}

X11()
plot(density(mu.hat.mle.star, n=50, window="g"),type="l",xlab="x",ylab="density")
X11()
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