MLE Normal

Here we fit the MLE for the normal

\(X_1, X_2, ...,X_n\) iid \(Normal(\mu, \sigma^2)\), find the MMEs.

p.53 \(f(x|\mu,\sigma)= \frac{1}{\sigma \sqrt{2\pi}}e^\frac{-(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

use the density estimator since this is a continuous random variable

hist(x.orig, probability=T)
lines(density(x.orig))

MLE fitted normal p.255

mu.hat.mle <- mean(x)
mu.hat.mle

[1] 70.01163

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

[1] 3.098744

library(fitdistrplus)

Loading required package: MASS

Loading required package: survival

fit.norm <- fitdist(x.orig, "norm", method = "mle")
summary(fit.norm)

Fitting of the distribution ' norm ' by maximum likelihood 
Parameters : 
      estimate Std. Error
mean 70.011634  0.5657506
sd    3.098744  0.4000459
Loglikelihood:  -76.49806   AIC:  156.9961   BIC:  159.7985 
Correlation matrix:
     mean sd
mean    1  0
sd      0  1

plot the fitted model to the simulated data

denscomp(list(fit.norm), legendtext = c("Normal"))

cdfcomp(list(fit.norm), legendtext = c("Normal"))

qqcomp(list(fit.norm), legendtext = c("Normal"))

ppcomp(list(fit.norm), legendtext = c("Normal"))

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

[1] 12.57640 13.91014

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

[1] 68.90276 71.12051

parametric bootstrap the MLE

fit.norm <- fitdist(x.orig, "norm", method = "mle")
summary(fit.norm)

Fitting of the distribution ' norm ' by maximum likelihood 
Parameters : 
      estimate Std. Error
mean 70.011634  0.5657506
sd    3.098744  0.4000459
Loglikelihood:  -76.49806   AIC:  156.9961   BIC:  159.7985 
Correlation matrix:
     mean sd
mean    1  0
sd      0  1

b1 <- bootdist(fit.norm, niter=101)
print(b1)

Parameter values obtained with parametric bootstrap 
      mean       sd
1 69.16186 2.961610
2 70.81909 2.343714
3 69.91526 3.790972
4 70.23152 3.080772
5 69.69590 2.965678
6 70.97888 2.990729

plot(b1)

plot(b1, enhance=TRUE)

summary(b1)

Parametric bootstrap medians and 95% percentile CI 
        Median     2.5%     97.5%
mean 70.061757 68.67938 70.978470
sd    2.965678  2.25645  3.803492

quantile(b1)

(original) estimated quantiles for each specified probability (non-censored data)
            p=0.1    p=0.2    p=0.3    p=0.4    p=0.5    p=0.6    p=0.7   p=0.8
estimate 66.04043 67.40367 68.38665 69.22658 70.01163 70.79669 71.63662 72.6196
            p=0.9
estimate 73.98283
Median of bootstrap estimates
            p=0.1    p=0.2    p=0.3    p=0.4    p=0.5   p=0.6    p=0.7    p=0.8
estimate 66.22283 67.53438 68.48182 69.28089 70.06176 70.7742 71.59122 72.48585
            p=0.9
estimate 73.68261

two-sided 95 % CI of each quantile
          p=0.1    p=0.2    p=0.3    p=0.4    p=0.5    p=0.6    p=0.7    p=0.8
2.5 %  64.54336 65.96744 66.94257 67.87320 68.67938 69.48556 70.26836 71.23212
97.5 % 67.55883 68.72942 69.57110 70.23741 70.97847 71.73471 72.50875 73.49470
          p=0.9
2.5 %  72.38655
97.5 % 75.04397

CIcdfplot(b1, CI.output = "quantile")