Log-normal Distribution

Intuitively, the lognormal distribution states that observed data after log-transformation follows a normal distribution.

\[F(t) = \Phi_{nor}(\frac {\log t-\mu}{\sigma})\]

\[f(t) = \frac {1}{\sigma t}\phi_{nor}(\frac {\log t-\mu}{\sigma}) = \frac {1}{\sigma t\sqrt{2\pi}}\exp(-\frac{1}{2}(\frac {\log t - \mu}{\sigma})^2)\]

\(\mu = \rm {location\ parameter}\)

\(\sigma = \rm {shape\ parameter}\)

3.

Based on the dataset \(ShockAbsorber.txt\) on moodle, use the lognormal distribution \(\rm {LN}(\mu, \sigma)\) to fit the data by the maximum likelihood (ML) method and provide the following estimates.

shock.absorber <- read.table("ShockAbsorber.txt", header = T)
shock.absorber$Censor <- ifelse(shock.absorber$Censor == "Failure", 1, 0)
head(shock.absorber)
##   Distance     Mode Censor
## 1     6700    Mode1      1
## 2     6950 Censored      0
## 3     7820 Censored      0
## 4     8790 Censored      0
## 5     9120    Mode2      1
## 6     9660 Censored      0

Parameters estimation.

\[\hat{\mu} = 10.145\]

\[\hat{\sigma} = 0.113\]

log.likelihood.ln <- function(mu, sigma, data){
  lj <- c()
  for(j in 1:length(data[,1])){
    likelihood <- data[j, 3]*dlnorm(data[j, 1], mu, sigma, log = TRUE) + (1-data[j, 3])*log(1-plnorm(data[j, 1], mu, sigma))
    lj <- c(lj, likelihood)
  }
  return(sum(lj))
}

minus.log.likelihood.ln <- function(pars, data){
  # if loop: for parameters <=0, the dweibull returns NaNs
  if (pars[1]<=0 | pars[2]<=0){
    return(100000)
  } else {
    return(-log.likelihood.ln(pars[1], pars[2], data))
  }
}


mu0 <- 10
sig0 <- 0.5
initial <- c(mu0, sig0)
par.hat <- optim(initial, minus.log.likelihood.ln, data = shock.absorber, hessian = TRUE)
est.par <- c(par.hat$par[1], par.hat$par[2])
est.par
## [1] 10.1447348  0.5301075

\[\hat{t_{0.1}} = \exp(\hat {\mu} + \Phi^{-1}_{nor}(0.1)\hat {\sigma}) = 12905.06\]

\[\hat F(10000) = \Phi_{nor}(\frac {\log 10000-\hat {\mu}}{\hat {\sigma}}) = 0.0389797\]

# t_p
t_0.1<- exp(est.par[1]+ qnorm(0.1)*est.par[2])
t_0.1
## [1] 12905.06
# F(10000)
F_10000<- pnorm((log(10000)-est.par[1])/est.par[2])
F_10000
## [1] 0.0389797
# plot
plot(density(shock.absorber$Distance), main = "", xlab = "t")
curve(dlnorm(x, est.par[1], (est.par[2])^2), add = T, col =2)
lines(c(t_0.1, t_0.1), c(0, dlnorm(t_0.1, est.par[1], (est.par[2])^2)), col = 3)

Standard error.

\[\hat {se}_{\hat {\mu}} = 0.1441874\]

\[\hat {se}_{\hat {\sigma}} = 0.1126943\]

cov.var.M<- round(solve(par.hat$hessian), 5)
se_mu<- sqrt(cov.var.M[1,1])
se_sigma<- sqrt(cov.var.M[2,2])
se_mu ; se_sigma
## [1] 0.1441874
## [1] 0.1126943

\[c = \left[ \begin{matrix} \frac {\partial t_p}{\partial \mu} \\ \frac {\partial t_p}{\partial \sigma} \end{matrix} \right] = \left. \left[ \begin{matrix} t_p \\ t_p\Phi^{-1}_{nor}(p) \end{matrix} \right]\right| _{\begin{matrix}\mu = \hat{\mu} \\ \sigma = \hat {\sigma} \\ p = 0.1\end{matrix}}\]

\[\hat {Var}_{\hat {t_p}} = \hat {c}^T \hat{\sum}_{\hat {\mu},\hat {\sigma}} \hat {c}\]

\[\hat {se}_{\hat {t_p}} = \sqrt{\hat {Var}_{\hat {t_p}}} = 1666.879\]

c <- matrix(t_0.1*c(1, qnorm(0.1)), nrow = 2)
se_t0.1<- as.numeric(sqrt(t(c)%*%cov.var.M%*%c))
se_t0.1
## [1] 1666.879

\[c = \left[ \begin{matrix} \frac {\partial F(t)}{\partial \mu} \\ \frac {\partial F(t)}{\partial \sigma} \end{matrix} \right] = \left. \left[ \begin{matrix} \frac {\phi_{nor}(z)}{-\hat {\sigma}} \\ \frac {\phi_{nor}(z)z}{-\hat {\sigma}} \end{matrix} \right]\right| _{\begin{matrix}\mu = \hat{\mu} \\ \sigma = \hat {\sigma} \\ t = 10000\end{matrix}} ,\ \ \ \ z = \frac {\log(10000) - \hat{\mu}}{\hat {\sigma}}\]

\[\hat {Var}_{\hat {F(t)}} = \hat {c}^T \hat{\sum}_{\hat {\mu},\hat {\sigma}} \hat {c}\]

\[\hat {se}_{\hat {F(t)}} = \sqrt{\hat {Var}_{\hat {F(t)}}} = 0.02562305\]

z<- (log(10000) - est.par[1])/est.par[2]
c <- matrix(dnorm(z)/-est.par[2]*c(1, z), nrow = 2)
se_F10000<- as.numeric(sqrt(t(c)%*%cov.var.M%*%c))
se_F10000
## [1] 0.02562305

Confidence interval.

\(\mu\)

mu_l<- est.par[1] - qnorm(0.975)*se_mu
mu_u<- est.par[1] + qnorm(0.975)*se_mu
c(mu_l, mu_u)
## [1]  9.862133 10.427337

\(\sigma\)

w <- exp(qnorm(0.975)*se_sigma/est.par[2])
sigma_l<- est.par[2]/w
sigma_u<- est.par[2]*w
c(sigma_l, sigma_u)
## [1] 0.3494693 0.8041162

\(t_{0.1}\)

w <- exp(qnorm(0.975)*se_t0.1/t_0.1)
tp_l<- t_0.1/w
tp_u<- t_0.1*w
c(tp_l, tp_u)
## [1] 10018.78 16622.84

\(F(10000)\)

w <- exp(qnorm(0.975)*se_F10000/(F_10000*(1-F_10000)))
F10000_l<- F_10000/(F_10000 + (1-F_10000)*w)
F10000_u<- F_10000/(F_10000 + (1-F_10000)/w)
c(F10000_l, F10000_u)
## [1] 0.01050253 0.13419942

Table:

estimates \(\mu\) \(\sigma\) \(t_p\) \(F(10000)\)
ML estimate 10.145 0.113 12905.06 0.0389797
Standard error of ML estimate 0.1441874 0.1126943 1666.879 0.02562305
Confidence interval [9.862, 10.427] [0.349, 0.804] [10018.78, 16622.84] [0.011, 0.134]

4.

Proof.