DOI: 10.5281/zenodo.17521050
title: Simulation and dataset-application code for the distribution "The Half-Logistic Generalized Power Lindley Distribution: Theory and Applications"

description: Simulation and dataset-application code for the distribution “The Half-Logistic Generalized Power Lindley Distribution: Theory and Applications” are provided:

library(ggplot2)library(goftest)library(stats4)

data <- c(  0.014, 0.034, 0.059, 0.061, 0.069, 0.080, 0.123, 0.142, 0.165, 0.210,  0.381, 0.464, 0.479, 0.556, 0.574, 0.839, 0.917, 0.969, 0.991, 1.064,  1.088, 1.091, 1.174, 1.270, 1.275, 1.355, 1.397, 1.477, 1.578, 1.649,  1.702, 1.893, 1.932, 2.001, 2.161, 2.292, 2.326, 2.337, 2.628, 2.785,  2.811, 2.886, 2.993, 3.122, 3.248, 3.715, 3.790, 3.857, 3.912, 4.100,  4.106, 4.116, 4.315, 4.510, 4.580, 5.267, 5.299, 5.583, 6.065, 9.701)stopifnot(all(is.finite(data)), all(data > 0))

## ---- Model: CDF and PDF (Lindley) ----## F(x;θ) = 1 - (1 + θ x/(θ+1)) exp(-θ x)## f(x;θ) = (θ^2/(θ+1)) (1+x) exp(-θ x)F_lindley <- function(x, theta){  if(!is.finite(theta) || theta <= 0) return(rep(NA_real_, length(x)))  x <- pmax(x, .Machine$double.xmin)  1 - (1 + theta*x/(theta+1)) * exp(-theta*x)}logpdf_lindley <- function(x, theta){  if(!is.finite(theta) || theta <= 0) return(rep(-Inf, length(x)))  if(any(!is.finite(x)) || any(x <= 0)) return(rep(-Inf, length(x)))  x <- pmax(x, .Machine$double.xmin)  2*log(theta) - log(theta+1) + log1p(x) - theta*x}f_lindley <- function(x, theta) exp(logpdf_lindley(x, theta))

Negative log-likelihood on log-parameter ## param = a with theta = exp(a)nll <- function(a){  theta <- exp(a)  llvec <- logpdf_lindley(data, theta)  if(any(!is.finite(llvec))) return(1e300)  -sum(llvec)}

MLE via optim (BFGS with Hessian) fit <- optim(log(1), nll, method = "BFGS", hessian = TRUE)if(fit$convergence != 0) warning("optim did not fully converge; consider different starts.")

theta_hat <- exp(fit$par[1])

H <- fit$hessian            # 1x1Vlog <- tryCatch(1/H[1,1], error=function(e) NA_real_)se_theta <- if(is.finite(Vlog) && Vlog > 0) sqrt(Vlog) * theta_hat else NA_real_

logL <- -fit$valuen <- length(data); k <- 1AIC  <- -2*logL + 2*kBIC  <- -2*logL + log(n)*kHQIC <- -2*logL + 2*k*log(log(n))CAIC <- -2*logL + k*(log(n) + 1)

ks_res <- ks.test(data, function(x) F_lindley(x, theta_hat))ad_res <- ad.test(data, null = F_lindley, theta = theta_hat)

cat("=== Estimated Parameters (Lindley) ===\n")cat(sprintf("theta_hat = %.6f  (SE = %s)\n",            theta_hat, ifelse(is.finite(se_theta), sprintf("%.6f", se_theta), "NA")))cat("\n=== Information Criteria ===\n")cat(sprintf("Log-Likelihood: %.6f\n", logL))cat(sprintf("AIC  : %.6f\n", AIC))cat(sprintf("BIC  : %.6f\n", BIC))cat(sprintf("HQIC : %.6f\n", HQIC))cat(sprintf("CAIC : %.6f\n", CAIC))cat("\n=== KS Test (vs fitted CDF) ===\n")cat(sprintf("KS Statistic: %.6f\np-value     : %.6f\n", ks_res$statistic, ks_res$p.value))cat("\n=== Anderson–Darling Test (goftest::ad.test) ===\n")cat(sprintf("AD Statistic: %.6f\np-value     : %.6f\n", ad_res$statistic, ad_res$p.value))

op <- par(mfrow = c(1, 2))

hist(data, probability = TRUE, breaks = 30,     main = "Lindley PDF",     col = "lightgray", xlab = "Data", border = "white")curve(f_lindley(x, theta_hat), add = TRUE, lwd = 2, col = "blue")

plot(ecdf(data), main = "Lindley CDF",     xlab = "Data", ylab = "Cumulative Probability")curve(F_lindley(x, theta_hat), add = TRUE, lwd = 2, col = "red")legend("bottomright", legend = c("Empirical CDF", "Fitted CDF"),       col = c("black", "red"), lwd = 2, bty = "n")

par(op)

***************************************

# G(x;θ,β) = [1 - A] / [1 + A],  A = (1 + θ x^β/(θ+1)) * exp(-θ x^β)G_cdf <- function(x, theta, beta){  z <- theta * x^beta  A <- (1 + z/(theta+1)) * exp(-z)  (1 - A) / (1 + A)}

logpdf_g <- function(x, theta, beta){  z    <- theta * x^beta  log_num <- log(2) + 2*log(theta) + log(beta) +             log1p(x^beta) - log(theta+1) + (beta-1)*log(x) - z  A <- (1 + z/(theta+1)) * exp(-z)  log_den <- 2*log1p(A)                 # log( (1+A)^2 )  log_num - log_den}g_pdf <- function(x, theta, beta){ exp(logpdf_g(x, theta, beta)) }# G(x) = (1 - A)/(1 + A) ⇒ A = (1 - u)/(1 + u)# A = ((θ+1)+t)/(θ+1) * e^{-t}, t = θ x^β# ((θ+1)+t) e^{-t} = A(θ+1) ⇒ s=t+(θ+1):  s e^{-s} = A(θ+1)e^{-(θ+1)}# s e^{-s} = c  ⇒  -s = W(-c). Gerçek çözüm için W_{-1} dalı gerekir.

lambertW_m1 <- function(z, tol=1e-12, maxit=50){  # z in (-1/e, 0). Çıktı W_{-1}(z) <= -1  if (any(z >= 0 | z <= -1/exp(1))) stop("lambertW_m1: z ∈ (-1/e, 0) olmalı.")  w <- log(-z) - log(-log(-z))   # iyi başlangıç (z ~ 0- için)  for(it in 1:maxit){    ew <- exp(w); f <- w*ew - z    w1 <- w + 1    denom <- ew*w1 - (w + 2)*f/(2*w1)    dw <- f/denom    w <- w - dw    if (max(abs(dw)) < tol) break  }  w}

rg_thlgpl <- function(n, theta, beta){  u <- runif(n)  u <- pmin(pmax(u, 1e-12), 1-1e-12)  A <- (1 - u)/(1 + u)  c <- A * (theta+1) * exp(-(theta+1))  z <- -c  # W_{-1} dalı:  w <- lambertW_m1(z)  s <- -w  t <- s - (theta+1)  y <- t / theta          # y = x^β  x <- y^(1/beta)  x}

nll_g <- function(par, dat){  a <- par[1]; b <- par[2]                theta <- exp(a); beta <- exp(b)  ll <- sum(logpdf_g(dat, theta, beta))  if(!is.finite(ll)) ll <- -1e300  -ll}

mle_once <- function(x, init=log(c(theta=1, beta=1))){  fit <- optim(init, nll_g, dat=x, method="L-BFGS-B",               control=list(maxit=200))  a <- fit$par[1]; b <- fit$par[2]  th_hat <- exp(a); be_hat <- exp(b)  H <- tryCatch(optimHess(fit$par, nll_g, dat=x),                error=function(e) matrix(NA_real_, 2, 2))  V <- tryCatch(solve(H), error=function(e) matrix(NA_real_, 2, 2))  se <- rep(NA_real_, 2)  if(all(is.finite(diag(V))) && all(diag(V) > 0)){    se_log <- sqrt(diag(V))    se <- se_log * c(th_hat, be_hat)   # delta yöntemi  }  list(est=c(theta=th_hat, beta=be_hat), se=c(theta=se[1], beta=se[2]))}

make_table_MLE <- function(theta, beta, n_vals=c(100,200,300,400), R=40,                           init=log(c(theta=1, beta=1))){  set.seed(123)  labs <- c("θ","β"); out <- list(); true <- c(theta, beta)  for(n in n_vals){    ests <- low <- up <- matrix(NA_real_, R, 2)    for(r in 1:R){      x <- rg_thlgpl(n, theta, beta)      m <- mle_once(x, init=init)      ests[r,] <- m$est      low[r,]  <- m$est - 1.96*m$se      up[r,]   <- m$est + 1.96*m$se    }    est_bar <- colMeans(ests, na.rm=TRUE)    bias    <- est_bar - true    mse     <- colMeans((ests - matrix(true, R, 2, byrow=TRUE))^2, na.rm=TRUE)    low_m   <- colMeans(low, na.rm=TRUE)    up_m    <- colMeans(up,  na.rm=TRUE)    len_m   <- up_m - low_m    out[[length(out)+1]] <- data.frame(      n=n, Parameter=labs,      ML_Estimate=est_bar, Bias=bias, MSE=mse,      Lower_Limit=low_m, Upper_Limit=up_m, Length=len_m,      row.names=NULL    )  }  do.call(rbind, out)}

gibbs_once <- function(x, niter=3000, burn=1000, thin=10,                       start=log(c(theta=1, beta=1)),                       prior=list(kt=1, lt=0.1, kb=1, lb=0.1), # Gam(k, l) oran: lt                       prop_sd=c(0.25, 0.25)){  a <- start[1]; b <- start[2]  keep <- seq(burn+1, niter, by=thin)  out  <- matrix(NA_real_, length(keep), 2); k <- 0  logpost <- function(a,b){    theta<-exp(a); beta<-exp(b)    lp <- sum(logpdf_g(x, theta, beta))    with(prior, { lp <- lp + kt*a - lt*theta + kb*b - lb*beta })    if(!is.finite(lp)) lp <- -1e300    lp  }  lcur <- logpost(a,b)  for(it in 1:niter){    ap <- rnorm(1,a,prop_sd[1]); lp <- logpost(ap,b); if(log(runif(1)) < lp - lcur){ a<-ap; lcur<-lp }    bp <- rnorm(1,b,prop_sd[2]); lp <- logpost(a,bp); if(log(runif(1)) < lp - lcur){ b<-bp; lcur<-lp }    if(it %in% keep){ k <- k+1; out[k,] <- c(exp(a),exp(b)) }  }  colnames(out) <- c("theta","beta")  list(mean=colMeans(out),       q=apply(out, 2, quantile, probs=c(.025,.5,.975)))}

make_table_GIBBS <- function(theta, beta,                             n_vals=c(150,200,300,400), R=30,                             niter=3000, burn=1000, thin=10){  set.seed(456)  labs <- c("θ","β"); out <- list(); true <- c(theta, beta)  for(n in n_vals){    ests <- low <- up <- matrix(NA_real_, R, 2)    for(r in 1:R){      x  <- rg_thlgpl(n, theta, beta)      bz <- gibbs_once(x, niter=niter, burn=burn, thin=thin)      ests[r,] <- bz$mean      low[r,]  <- bz$q[1,]; up[r,] <- bz$q[3,]    }    est_bar <- colMeans(ests, na.rm=TRUE)    bias    <- est_bar - true    mse     <- colMeans((ests - matrix(true, R, 2, byrow=TRUE))^2, na.rm=TRUE)    low_m   <- colMeans(low, na.rm=TRUE)    up_m    <- colMeans(up,  na.rm=TRUE)    len_m   <- up_m - low_m    out[[length(out)+1]] <- data.frame(      n=n, Parameter=labs,      Gibbs_Estimate=est_bar, Bias=bias, MSE=mse,      Lower_Limit=low_m, Upper_Limit=up_m, Length=len_m,      row.names=NULL    )  }  do.call(rbind, out)}set.seed(7)

tab_gibbs <- make_table_GIBBS(theta=1.5, beta=0.9,                              n_vals=c(100,200,300,400), R=10,                              niter=3000, burn=1000, thin=10)

tab_mle   <- make_table_MLE(theta=1.5, beta=0.9,                            n_vals=c(100,200,300,400), R=10,                            init=log(c(theta=1.2, beta=0.8)))  # iyi başlangıç

print(tab_gibbs)print(tab_mle)

****************************

## MLE & Bayes (RW-MH)make_summary_table <- function(estimates, true_vals, lows, ups, method_name){  # estimates: R x p (R tekrar, p parametre)  # lows, ups: R x p  (alt/üst CI)  p <- length(true_vals)  est_bar <- colMeans(estimates, na.rm=TRUE)  bias    <- est_bar - true_vals  mse     <- colMeans((estimates - matrix(true_vals, nrow(estimates), p, byrow=TRUE))^2, na.rm=TRUE)  low_m   <- colMeans(lows, na.rm=TRUE)  up_m    <- colMeans(ups, na.rm=TRUE)  len_m   <- up_m - low_m  data.frame(    Parameter = colnames(estimates),    Estimate  = est_bar,    Bias      = bias,    MSE       = mse,    Lower     = low_m,    Upper     = up_m,    Length    = len_m,    Method    = method_name,    row.names = NULL  )}

make_table_MLE_fmt <- function(theta, beta, n_vals=c(100,200,300,400), R=50){  true <- c(theta, beta)  labs <- c("θ","β")  out <- list()  for(n in n_vals){    ests <- lows <- ups <- matrix(NA_real_, R, 2, dimnames=list(NULL,labs))    for(r in 1:R){      x <- rg_newmodel(n, theta, beta)      fit <- mle_once(x)      ests[r,] <- fit$est      if(all(is.finite(fit$se))){        lows[r,] <- fit$est - 1.96*fit$se        ups[r,]  <- fit$est + 1.96*fit$se      }    }    tab <- make_summary_table(ests, true, lows, ups, "MLE")    tab$n <- n    out[[length(out)+1]] <- tab  }  do.call(rbind, out)}

make_table_BAYES_fmt <- function(theta, beta, n_vals=c(100,200,300,400), R=50,                                 niter=3000, burn=1000, thin=10){  true <- c(theta, beta)  labs <- c("θ","β")  out <- list()  for(n in n_vals){    ests <- lows <- ups <- matrix(NA_real_, R, 2, dimnames=list(NULL,labs))    for(r in 1:R){      x  <- rg_newmodel(n, theta, beta)      bz <- rw_mh_once(x, niter=niter, burn=burn, thin=thin)      ests[r,] <- bz$mean      lows[r,] <- bz$q[1,]      ups[r,]  <- bz$q[3,]    }    tab <- make_summary_table(ests, true, lows, ups, "Bayes")    tab$n <- n    out[[length(out)+1]] <- tab  }  do.call(rbind, out)}

set.seed(77)tab_mle   <- make_table_MLE_fmt(theta=4.5, beta=0.8,                                n_vals=c(100,200,300,400), R=30)tab_bayes <- make_table_BAYES_fmt(theta=4.5, beta=0.8,                                  n_vals=c(100,200,300,400), R=30)

print(tab_mle)print(tab_bayes)

********************

## HL–G Power Lindley — MLE + Simple Bayesian data <- c(  0.1405, 0.1566, 0.1577, 0.1604, 0.1608, 0.2215, 0.2994, 0.3131, 0.3246, 0.3247,  0.3295, 0.3300, 0.3379, 0.3397, 0.3523, 0.3589, 0.3933, 0.4176, 0.4258, 0.4356,  0.4421, 0.4444, 0.4505, 0.4558, 0.4683, 0.4733, 0.4846, 0.4889, 0.5096, 0.5177,  0.5278, 0.5347, 0.5433, 0.5442, 0.5508, 0.5527, 0.5606, 0.5607, 0.5671, 0.5753,  0.5828, 0.6030, 0.6050, 0.6136, 0.6261, 0.6395, 0.6469, 0.6512, 0.6816, 0.6994,  0.7048, 0.7292, 0.7430, 0.7455, 0.7798, 0.7984, 0.8147, 0.8230, 0.8302, 0.8342,  0.9794)stopifnot(all(is.finite(data)), all(data > 0))## G(x;θ,β)= (1 - A)/(1 + A),  A = (1 + θ x^β/(θ+1)) exp(-θ x^β)G_cdf <- function(x, theta, beta){  if (!is.finite(theta) || !is.finite(beta) || theta <= 0 || beta <= 0)    return(rep(NA_real_, length(x)))  x <- pmax(x, .Machine$double.xmin)  z <- theta * x^beta  A <- (1 + z/(theta+1)) * exp(-z)  (1 - A)/(1 + A)}

## log g(x;θ,β) = log( 2 * θ^2 * β/(θ+1) * x^{β-1} * (1+x^β) * e^{-θ x^β} ) - 2*log(1+A)logpdf_HLPL <- function(x, theta, beta){  if (!is.finite(theta) || !is.finite(beta) || theta <= 0 || beta <= 0)    return(rep(-Inf, length(x)))  if (any(!is.finite(x)) || any(x <= 0)) return(rep(-Inf, length(x)))  x <- pmax(x, .Machine$double.xmin)  z <- theta * x^beta  log_num <- log(2) + 2*log(theta) + log(beta) - log(theta+1) +    (beta-1)*log(x) + log1p(pmin(x^beta, 1e300)) - z  A <- (1 + z/(theta+1)) * exp(-z)  log_den <- 2*log1p(pmin(A, 1e300))  log_num - log_den}g_pdf <- function(x, theta, beta) exp(logpdf_HLPL(x, theta, beta))

nll <- function(par_log){  th <- exp(par_log[1]); be <- exp(par_log[2])  lv <- logpdf_HLPL(data, th, be)  if (any(!is.finite(lv))) return(1e300)  -sum(lv)}start_log <- log(c(theta=1, beta=1))fit_mle <- optim(start_log, nll, method="BFGS", hessian=TRUE)if(fit_mle$convergence != 0) warning("optim did not fully converge; consider different starts.")

theta_hat <- exp(fit_mle$par[1]); beta_hat <- exp(fit_mle$par[2])

H <- fit_mle$hessianV <- tryCatch(solve(H), error=function(e) matrix(NA_real_, 2, 2))se_theta <- se_beta <- NA_real_if (all(is.finite(V)) && all(diag(V) > 0)){  se_log <- sqrt(diag(V))  se_theta <- se_log[1] * theta_hat  se_beta  <- se_log[2] * beta_hat}

logL <- -fit_mle$valuen <- length(data); k <- 2AIC  <- -2*logL + 2*kBIC  <- -2*logL + log(n)*kHQIC <- -2*logL + 2*k*log(log(n))CAIC <- -2*logL + k*(log(n) + 1)

ks_res <- ks.test(data, function(xx) G_cdf(xx, theta_hat, beta_hat))

ad_stat <- NA_real_; ad_p <- NA_real_if (isTRUE(has_goftest)) {  adr <- goftest::ad.test(data, null=G_cdf, theta=theta_hat, beta=beta_hat)  ad_stat <- as.numeric(adr$statistic); ad_p <- as.numeric(adr$p.value)} else {  u <- sort(pmin(pmax(G_cdf(data, theta_hat, beta_hat), 1e-12), 1-1e-12))  i <- seq_along(u); n <- length(u)  ad_stat <- -n - (1/n) * sum((2*i - 1) * (log(u) + log(1 - rev(u))))  ad_p <- NA_real_}

bayes_hlgpl <- function(x,                        niter=15000, burn=5000, thin=10,                        start=c(theta=1, beta=1),                        prior=list(kt=1, lt=0.1, kb=1, lb=0.1), # Gamma(shape, rate)                        prop_sd=c(0.25, 0.25), seed=123){

  set.seed(seed)  logpost <- function(a,b){    th <- exp(a); be <- exp(b)    lp <- sum(logpdf_HLPL(x, th, be))    if(!is.finite(lp)) return(-Inf)    with(prior, lp + (kt-1)*a - lt*th + (kb-1)*b - lb*be + (a+b)) # + Jacobian  }

  a <- log(start["theta"]); b <- log(start["beta"])  cur <- logpost(a,b)  if(!is.finite(cur)){    a <- log(1); b <- log(1); cur <- logpost(a,b)    if(!is.finite(cur)) stop("Start/prior invalid.")  }

  keep <- seq(burn+1, niter, by=thin)  out  <- matrix(NA_real_, nrow=length(keep), ncol=2, dimnames=list(NULL,c("theta","beta")))  kacc_a <- 0L; kacc_b <- 0L; jj <- 0L

  for(it in 1:niter){    ap <- rnorm(1, a, prop_sd[1]); lp <- logpost(ap,b)    if (is.finite(lp) && log(runif(1)) < lp - cur){ a <- ap; cur <- lp; kacc_a <- kacc_a+1L }    bp <- rnorm(1, b, prop_sd[2]); lp <- logpost(a,bp)    if (is.finite(lp) && log(runif(1)) < lp - cur){ b <- bp; cur <- lp; kacc_b <- kacc_b+1L }    if (it %in% keep){ jj <- jj+1L; out[jj,] <- c(exp(a), exp(b)) }  }

  list(draws=out,       accept=c(a=kacc_a/niter, b=kacc_b/niter),       mean=colMeans(out),       sd=apply(out,2,sd),       ci95=apply(out,2,quantile, probs=c(.025,.975)))}

fit_bayes <- bayes_hlgpl(  x=data,  niter=15000, burn=5000, thin=10,  start=c(theta=theta_hat, beta=beta_hat),   # MLE başlangıç  prior=list(kt=1, lt=0.10, kb=1, lb=0.10),  prop_sd=c(0.25, 0.25),  seed=7)

cat("\n=== MLE RESULTS — HL–G PL ===\n")cat(sprintf("theta_hat = %.6f (SE= %s)\n", theta_hat, ifelse(is.finite(se_theta), sprintf("%.6f",se_theta),"NA")))cat(sprintf("beta_hat  = %.6f (SE= %s)\n", beta_hat,  ifelse(is.finite(se_beta),  sprintf("%.6f",se_beta), "NA")))cat(sprintf("LogL= %.6f | AIC= %.6f | BIC= %.6f | HQIC= %.6f | CAIC= %.6f\n", logL, AIC, BIC, HQIC, CAIC))cat(sprintf("KS: D= %.6f, p= %.6f\n", ks_res$statistic, ks_res$p.value))cat(sprintf("AD: A^2= %.6f, p= %s\n", ad_stat, ifelse(is.na(ad_p),"NA",sprintf("%.6f",ad_p))))

cat("\n=== Bayesian (RW–MH) RESULTS — HL–G PL ===\n")cat(sprintf("Accept rates  : a=%.3f, b=%.3f\n", fit_bayes$accept[1], fit_bayes$accept[2]))cat(sprintf("Posterior mean: theta=%.6f, beta=%.6f\n", fit_bayes$mean[1], fit_bayes$mean[2]))cat(sprintf("Posterior sd  : theta=%.6f, beta=%.6f\n", fit_bayes$sd[1], fit_bayes$sd[2]))cat(sprintf("95%% CrI theta : [%.6f, %.6f]\n", fit_bayes$ci95[1,1], fit_bayes$ci95[2,1]))cat(sprintf("95%% CrI beta  : [%.6f, %.6f]\n", fit_bayes$ci95[1,2], fit_bayes$ci95[2,2]))

tab_compare <- data.frame(  Parameter         = c("θ","β"),  MLE_Estimate      = c(theta_hat, beta_hat),  MLE_SE            = c(se_theta,  se_beta),  Bayesian_Estimate = c(fit_bayes$mean[1], fit_bayes$mean[2]),  Bayesian_SE       = c(fit_bayes$sd[1],   fit_bayes$sd[2]),  check.names = FALSE)num_cols <- sapply(tab_compare, is.numeric)tab_rounded <- tab_comparetab_rounded[ , num_cols] <- lapply(tab_rounded[ , num_cols], function(v) round(v, 6))

cat("\n=== Table: MLE vs Bayesian (Estimates & SE) ===\n")print(tab_rounded, row.names = FALSE)