OLiHL Distribution

Introduction

This R markdown document is a brief study of the Odd Lindley Half Logistic Distribution, which was first proposed in the article “A new one-parameter lifetime distribution and its regression model with applications”, from the authors:

Eliwa MS, Altun E, Alhussain ZA, Ahmed EA, Salah MM, Ahmed HH, et al. (2021) A new one-parameter lifetime distribution and its regression model with applications. PLoS ONE 16(2): e0246969. https://doi.org/10.1371/journal.pone.0246969

The Odd Lindley Half Logistic Distribution probability density function

\[\begin{equation} f(x; \lambda) \boldsymbol{=} \frac{\lambda^{2}(1+e^x)}{4(1+\lambda)} e^{ \boldsymbol{-}\frac{\lambda}{2}(e^{x}-1) + x}; x > 0 \end{equation}\]

The distribution is implemented as it follows

Random Variable generation from quantile function

set.seed(42)

lambda0 = 2
n = 10000
u = runif(n,0,1)

quantile_func <- function(u,param){
  Q = log( -2/param *(1 + lambertWm1((1+param)*(u-1)*exp(-1-param))) - 1)
  
}

x = quantile_func(u,lambda0)

Empirical plotting

data_dist <- data.frame(u = rep(u, 3),lambda = factor(rep(1:3, each = length(u))),
                   xvalues = unlist(lapply(1:3, function(lambda) quantile_func(u, lambda))))


ggplot(data_dist, aes(x = xvalues, fill = lambda, color = lambda)) +
  geom_histogram(aes(y = ..density..), position = "identity", alpha = 0.2, bins = 30) +
  geom_density(aes(color = lambda), alpha = 0, size = 1) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1") +
  xlim(0, 4) +
  ylim(0, 1.2) +
  labs(color = "Lambda", fill = "Lambda", title = "Empirical Density and Histograms", x = "Values", y = "Density") +
  theme_light()

ggplot(data_dist, aes(x = xvalues, fill = lambda, color = lambda))+
  stat_ecdf(geom = "step", lwd = 1)+
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1") +
  xlim(0, 4) +
  ylim(0, 1) +
  labs(color = "Lambda", fill = "Lambda", title = "Empirical CDF", x = "Values", y = "CDF") +
  theme_light()

Theoretical pdf, cdf, hrf

pdf <- function(x,shape){
  fx = ((shape**2)*(1+exp(x)))/(4*(1+shape)) * exp( -(shape/2)*(exp(x) - 1) + x )
}

cdf <- function(x,shape){
  Fx = 1 - ( ( shape*exp(x) + shape +2 ) / ( 2*(shape + 1) ) ) * exp( -(shape/2)*(exp(x) - 1)) 
}

hrf <- function(x,shape){
  h = ((shape**2)*(1+exp(x))) / ( 2*( shape*(1+exp(-x)) + 2*exp(-x)) )
}

space <- seq(0, 4, length.out = 1000)
space_hrf <- seq(0, 5, length.out = 1000)

PDF

data_pdf <- data.frame(
  x = rep(space, 5),
  lambda = factor(rep(c(1, 2, 3, 1.5, 0.5), each = length(space))),
  y = c(pdf(space, 1), pdf(space, 2), pdf(space, 3), pdf(space, 1.5), pdf(space, 0.5))
)

ggplot(data_pdf, aes(x = x, y = y, color = lambda)) +
  geom_line(size = 1) +
  scale_color_brewer(palette = "Set1") +
  xlim(0, 4) +
  ylim(0, 1.2) +
  labs(color = expression(lambda), title = "pdf OLiHL", x = "Values", y = "Density") +
  theme_classic()

CDF

data_cdf <- data.frame(
  x = rep(space, 5),
  lambda = factor(rep(c(1, 2, 3, 1.5, 0.5), each = length(space))),
  y = c(cdf(space, 1), cdf(space, 2), cdf(space, 3), cdf(space, 1.5), cdf(space, 0.5))
)

ggplot(data_cdf, aes(x = x, y = y, color = lambda)) +
  geom_line(size = 1) +
  scale_color_brewer(palette = "Set1") +
  xlim(0, 4) +
  ylim(0, 1) +
  labs(color = expression(lambda), title = "cdf OLiHL", x = "Values", y = "Density")+
  theme_classic()

HRF

data_hrf <- data.frame(
  x = rep(space_hrf, 5),
  lambda = factor(rep(c(1, 2, 3, 1.5, 0.5), each = length(space))),
  y = c(hrf(space_hrf, 1), hrf(space_hrf, 2), hrf(space_hrf, 3), hrf(space_hrf, 1.5), hrf(space_hrf, 0.5))
)

ggplot(data_hrf, aes(x = x, y = y, color = lambda)) +
  geom_line(size = 1) +
  scale_color_brewer(palette = "Set1") +
  ylim(0,25) +
  labs(color = expression(lambda), title = "hrf OLiHL", x = "Values", y = "Density")+
  theme_classic()

ESTIMATION

MLE for λ : 2.000518

Building empirical Standard Error and Confidence Interval

If you are maximizing a likelihood, then the covariance matrix of the estimates (observed information) is (asymptotically normal) the inverse of the negative of the Hessian. The standard errors are the square roots of the diagonal elements of the covariance matrix.

inv_fish <- solve(hessian)
se <- sqrt(inv_fish[1,1])
ub <- mle_lambda + (qnorm(0.025,0,1,lower.tail = FALSE)*se)
lb <- mle_lambda - (qnorm(0.025,0,1,lower.tail = FALSE)*se)

cat("CI at 95% confidence level = ", "[",lb,ub,"]")

CI at 95% confidence level =  [ 1.96908 2.031957 ]

Monte Carlo Simulation for sample size n = 25

Set a sample size \(n_{mc}\);
Set a number of simulations (sampling)
Initialize parameters \(\lambda_{i}\)
Generate random variables from distribution
Perform calculations for mle, bias, mse, approximate confidence intervals, and proportion of estimates contained in CI

set.seed(42)

n_mc = 25
n_sims = 1000

parameters <- c(0.1, 0.5, 1, 1.5, 3, 5)

results_mc <- list()
bias_mc <- list()
mse_mc <- list()
mle_mc <- list()
ci_mc <- list()
hessian_mc <- list()
proportion_ci_mc <- list()

nloglike_mc <- function(l, x) {
  n <- length(x)
  -1*(2*n*log(l)-n*log(4+4*l)+sum(log(1+exp(x)))-(l/2)*sum(exp(x)-1)+sum(x))
}

for (param in parameters) {
  
  u_mc <- matrix(runif(n_sims * n_mc, 0, 1), nrow = n_sims, ncol = n_mc)
  sims_mc <- apply(u_mc, c(1, 2), quantile_func, param = param)
  results_mc[[as.character(param)]] <- sims_mc
  
  mle_estimates_mc <- numeric(nrow(sims_mc))
  hessian_estimates_mc <- list()
  ci_param <- matrix(NA, nrow = nrow(sims_mc), ncol = 2)
  
  for ( i in 1:nrow(sims_mc) ) {
    x <- sims_mc[i, ]
    result_mle_mc <- optim(1, nloglike_mc, x=x, method="L-BFGS-B", lower=0.01, upper=100, hessian=TRUE)
    mle_estimates_mc[i] <- result_mle_mc$par
    hessian_estimates_mc[[i]] <- result_mle_mc$hessian
    
    se <- (diag(solve(result_mle_mc$hessian)))
    ci_param[i, ] <- mle_estimates_mc[i] + qnorm(c(0.025, 0.975)) * sqrt(se)
  }

  
  mle_mc[[as.character(param)]] <- mle_estimates_mc
  
  hessian_mc[[as.character(param)]] <- hessian_estimates_mc
  
  ci_mc[[as.character(param)]] <- ci_param
  
  within_ci <- apply(ci_param, 1, function(ci) param >= ci[1] && param <= ci[2])
  proportion_ci_mc[[as.character(param)]] <- mean(within_ci)
  
  bias_mc[[as.character(param)]] <- mean(mle_estimates_mc) - param
    
  mse_mc[[as.character(param)]] <- mean((mle_estimates_mc - param)^2)
  
}