Modified BAS iteration method for absolute value equation

1.   Introduction

Lifetime distributions are integral to reliability theory and survival analysis, offering critical tools for analyzing real-world data in fields such as finance, manufacturing, biological sciences, physics, and engineering ^[1,2]. Traditional distributions like the exponential and Weibull are valued for their simplicity and interpretability. The exponential distribution, with its assumption of a constant failure rate, is often used to describe system reliability at the component level ^[3,4,5]. However, real-world systems frequently exhibit non-constant hazard rates, requiring more flexible distributions to capture increasing or decreasing failure behaviors. Although the Weibull distribution can accommodate both increasing and decreasing failure rates, it often proves inadequate for complex datasets involving multiple failure causes or censored data, necessitating more advanced techniques to accurately represent intricate phenomena ^[6].

To address the limitations of classical lifetime distributions in capturing the diverse behaviors observed in real-world reliability and survival data, recent research has increasingly focused on constructing new families of distributions through the compounding of continuous lifetime distributions with truncated discrete distributions ^[7,8,9]. This innovative approach enhances the flexibility and applicability of standard techniques by introducing additional parameters that allow for better representation of skewness, kurtosis, and hazard rate shapes. Specifically, continuous distributions such as the exponential, Weibull, Rayleigh, and gamma have been compounded with discrete distributions like the Poisson, geometric, binomial, negative-binomial, logarithmic, and more generally, power series distributions. The foundational work by Adamidis & Loukas ^[10], which introduced the exponential-geometric distribution, laid the groundwork for a wide class of compound distributions. This approach was later expanded by Adamidis et al. ^[11] and Kuş ^[12], who explored further generalizations and alternative mixing distributions, leading to a surge in the development of new techniques. These compound distributions exhibit greater flexibility, enabling better fitting of datasets with complex features such as increasing, decreasing, bathtub-shaped, or upside-down bathtub hazard functions. Notable contributions in this domain include extended versions of well-known lifetime distributions that have demonstrated improved performance in applications spanning reliability engineering, biostatistics, and actuarial science ^{[13,14,15,16,17]}. These distributions not only offer improved statistical properties—such as the ability to accommodate over-dispersion and unobserved heterogeneity—but also maintain mathematical tractability, allowing for explicit expressions of the probability density, cumulative distribution, moments, and hazard functions. For a comprehensive overview and classification of these compounding models and their theoretical properties, interested readers may refer to the reviews by Goyal et al. ^[7], Maurya & Nadarajah ^[8], Rahmouni & Orabi ^[18,19], and Tahir & Cordeiro ^[9].

Building on these advancements, we focus on the need to analyze intermediate failure events, as prior studies often concentrated on the minimum lifetime, $X_{(1)} = \min\{Y_i\}_{i = 1}^{N}$, or the maximum lifetime, $X_{(N)} = \max\{Y_i\}_{i = 1}^{N}$, of a system's components. Many real-world applications, however, demand the analysis of the $k$-th failure. Order statistics, arranged in ascending order as $X_{(1)} \leq X_{(2)} \leq \cdots \leq X_{(N)}$, are essential in reliability analysis ^{[20,21,22,23]} and also apply to economics, where they provide distributional criteria for evaluating welfare and inequalities in income and wealth ^[24,25]. To address this need, we propose the Weibull-generalized shifted geometric (WGSG) distribution, a new family of compound distributions that generalizes the Weibull-geometric (WG) distribution introduced by Barreto-Souza et al. ^[26]. The WG distribution describes the time to the first failure in a system with a random number of components, $N$, where each component's lifetime, $Y_i$, follows a Weibull distribution, and $N$ follows a geometric distribution. Our approach extends this model to capture the $k$-th order statistic, $X_{(k)}$, enabling the analysis of any failure order, such as the second, third, or $k$-th lifetime. To illustrate this approach, consider a machine producing $N$ units (e.g., light bulbs or wire fuses), each with an independent Weibull-distributed lifetime $Y_i$ ($i = 1, 2, \ldots, N$). If $N$ follows a shifted geometric distribution—adjusting the starting point of the geometric progression to $n \geq k$ while preserving its properties—the $k$-th order statistic $X_{(k)}$ provides a versatile method for analyzing system reliability. This shifted geometric distribution is particularly valuable in scenarios where counting begins at a specific point $k$, such as in quality control, queueing theory, inventory management, and reliability systems with latent defects ^[27].

The WGSG distribution introduces additional shape parameters to capture a broad spectrum of failure behaviors, making it well-suited for analyzing the $k$-th lifetime in systems where failure occurs after $k$ units fail. This model aligns with the analysis of $k$-out-of-$n$ systems, where a system fails when $n-k+1$ or more components fail ^{[28,29,30,31]}. Unlike traditional approaches that assume binary states (fully functional or failed), the WGSG distribution can accommodate degraded states in multi-state systems, which are prevalent in practical applications ^{[32,33,34,35]}.

This paper investigates the properties of the WGSG distribution, explores parameter estimation techniques, and demonstrates its practical applications using real datasets. Building on the compounding methodology of Adamidis & Loukas ^[10] and Kuş ^[12], we aim to provide a versatile tool for reliability and survival analysis that addresses the limitations of traditional approaches in capturing complex failure dynamics. The paper is structured as follows: Section 2 introduces the WGSG distribution and its probability density function (PDF). Section 3 discusses its properties, including reliability, the failure rate function, the moment generating function, $r$-th moments, and Shannon entropy. Section 4 presents parameter estimation approaches, such as maximum likelihood estimation (MLE), the expectation-maximization (EM) algorithm, and method of moments. Section 5 presents the simulation study and Section 6 illustrates the WGSG distribution's practicality with real datasets. Finally, Section 7 concludes the paper and suggests directions for future research.

2.   The WGSG distribution

Let $X$ be a random variable following the WGSG distribution, parameterized by $\varphi = (\alpha, \theta, p)$, where $\alpha > 0$, $\theta = \frac{1}{\beta} > 0$, and $p \in (0, 1)$. Then, the PDF of the WGSG distribution is:

The cumulative distribution function (CDF) is given by:

This formulation characterizes the mixture between the shifted geometric distribution and the $k^{\text{th}}$-order statistic of Weibull-distributed random variables. A notable special case occurs when $k = 1$, where the WGSG density simplifies to the Weibull-geometric (WG) distribution, which models the time to the first failure, $X_{(1)} = \min\{Y_i\}_{i = 1}^{N}$, as introduced in Barreto-Souza et al. ^[26]. Figures 1 and 2 display the PDF and CDF of the WGSG distribution for parameter values $\alpha \in \{0.3, 0.5, 1, 1.5, 2, 2.5\}$, $\theta = 1$, $p = 0.3$, and $k \in \{1, 2, 3, 4, 5\}$. These plots illustrate the distribution's flexibility in modeling a wide range of reliability scenarios. When $\alpha < 1$ (e.g., 0.3, 0.5), the PDF exhibits a sharp decline and a singularity at $x = 0$, capturing early failure behavior. In contrast, for $\alpha > 1$ (e.g., 1.5, 2, 2.5), the PDF becomes unimodal, with its peak shifting leftward and becoming sharper as $\alpha$ increases. Higher values of $k$ yield flatter peaks and more rapid tail decay in the PDF, reflecting the compounding effect of higher-order statistics. Similarly, the CDF becomes steeper as $\alpha$ and $k$ increase, indicating lighter tails and a greater concentration of probability mass near the origin. These characteristics underline the WGSG distribution's utility in reliability engineering as well as in survival analysis involving extreme event times and in risk assessment for rare failures. The parameters $k$ and $p$ offer fine control over the distribution's shape, hazard rate, and tail behavior, enhancing its applicability across various practical contexts.

The Weibull-like term $\alpha \theta^\alpha x^{\alpha-1} e^{-(\theta x)^\alpha}$ governs the overall shape of the PDF, typically exhibiting unimodal behavior—initially increasing, reaching a peak, and then decreasing. The geometric component, involving the parameters $p$ and $k$, adjusts the probability structure and influences the tail behavior potentially making it heavier or lighter depending on $p$ and $k$. For small values of $x$, the Weibull-like term dominates, causing the PDF to increase as $x$ grows, with the exponential term $e^{-(\theta x)^\alpha}$ having minimal effect. As $x$ continues to increase, the exponential term starts to decay and the PDF reaches a maximum at some intermediate value, where the Weibull-like component and the geometric factors balance. For large $x$, the exponential term takes over, causing the PDF to decline. The geometric component further ensures that the distribution tail decays smoothly without forming multiple modes.

To establish that the WGSG distribution is unimodal, we examine the first derivative of the PDF with respect to $x$. The critical points of the PDF, where the derivative equals zero, indicate the mode. The first derivative is positive for small $x$, showing an increasing density, and negative for larger $x$, indicating a decreasing density. The presence of a single critical point where the first derivative changes sign from positive to negative suggests unimodality. Further verification via the second derivative or concavity analysis strengthens this conclusion. Therefore, the WGSG distribution is unimodal when $\alpha > 1$, characterized by a single peak followed by a steady decline, with no secondary modes present.

As $p \to 1$, the WGSG distribution becomes a Weibull distribution with shape parameter $\alpha$ and scale parameter $\theta$ because the shifted geometric component effectively becomes deterministic, i.e., $P(N = k) \to 1$. The mixture model collapses to a single Weibull component, and the number of components (determined by $N$) becomes one. The limiting behavior is given by:

Proof of the PDF and CDF expressions in Eqs (2.1) and (2.2): Let $Y_1, Y_2, \dots, Y_n$ be a random sample of size $n$ from a Weibull distribution. The $k$-th order statistic, $X_{(k)}$, is the $k$-th smallest value in the sample. The PDF of $X_{(k)}$ is:

where $f(x; \alpha, \theta) = \alpha \theta^{\alpha} x^{\alpha-1} e^{-(\theta x)^\alpha}$ is the PDF of a Weibull random variable and $F(x) = 1 - e^{-(\theta x)^\alpha}$ is the CDF with scale parameter $\theta > 0$ and shape parameter $\alpha > 0$. Thus, the PDF of the $k$-th order Weibull distribution is:

The probability mass function (PMF) of the distribution of $N$ that begins counting from $k$ is:

where $p$ is the probability of success on each trial. The marginal PDF is obtained by summing over all possible values of $N$:

Substituting $f_{X_{(k)}}(x; \alpha, \theta, n)$ and $P(N = n)$:

Let $n = m + k$, where $m = n - k$. The summation becomes:

Using the rising factorial (Pochhammer symbol) $(k+1)_m = \frac{\Gamma(m+k+1)}{\Gamma(k+1)}$, the summation can be expressed as:

This series is the generalized hypergeometric function ${}_1F_0(-; k; z)$, where $z = p e^{-(\theta x)^\alpha}$. Thus:

Substituting back:

Then,

where:

Let $u = 1 - e^{-(\theta t)^\alpha}$. Then, $du = \alpha \theta^\alpha t^{\alpha-1} e^{-(\theta t)^\alpha} \, dt$.

The PDF in terms of $u$ is:

Thus, the CDF becomes:

3.   Proprieties

This section characterizes the WGSG distribution by examining its survival and hazard functions (including asymptotic tails), moment generating function, $r$-th moments, limiting behavior, relationships with standard distributions in reliability and survival analysis (e.g., Weibull, geometric, and exponential), and entropy.

3.1.   Survival and hazard functions

The survival function $S_X(x)$ of the WGSG distribution represents the probability that a random variable $X$ exceeds a given value $x$:

The corresponding hazard function $h_X(x)$, which quantifies the instantaneous failure rate, is given by:

The tail behavior of the WGSG distribution governs the likelihood of extreme values, critical for risk assessment. As $x \to \infty$, $e^{-(\theta x)^\alpha} \to 0$. Letting $z = e^{-(\theta x)^\alpha}$, $F_X(x) = \left(\frac{1 - z}{1 - p z} \right)^k$, so $F_X(x) \to 1$ and $S_X(x) \to 0$, indicating a light-tailed distribution.

For large $x$, $1 - e^{-(\theta x)^\alpha} \to 1$ and $1 - p e^{-(\theta x)^\alpha} \to 1$, yielding $f_X(x) \sim k (1 - p) \alpha \theta^\alpha x^{\alpha - 1} e^{-(\theta x)^\alpha}$. The exponential term drives rapid decay, modulated by $x^{\alpha - 1}$, with the tail exhibiting exponential decay for $\alpha = 1$, subexponential decay for $\alpha > 1$, or stretched exponential decay for $\alpha < 1$, accommodating diverse extreme event scenarios. To characterize the tail and hazard behavior, we derive asymptotic expansions for $S_X(x)$ and $h_X(x)$.

Theorem 1 (Right-tail and hazard asymptotics). For $X \sim \mathit{\text{WGSG}}(\alpha, \theta, p, k)$, the following hold:

(1) As $x \to \infty$, the survival function satisfies:

(2) The hazard function behaves as:

Proof. (Right-tail asymptotics) Let $z = e^{-(\theta x)^\alpha}$, so $z \to 0^+$ as $x \to \infty$. The survival function is:

Expand:

Using $(1 - p z)^{-1} = 1 + p z + p^2 z^2 + \mathcal{O}(z^3)$:

Compute:

Thus:

Factorize:

Simplify the coefficient:

With $z = e^{-(\theta x)^\alpha}$:

□

Proof. (Hazard rate asymptotics) As $x \to \infty$, using $f_X(x) \sim k (1 - p) \alpha \theta^\alpha x^{\alpha - 1} e^{-(\theta x)^\alpha}$ and $S_X(x) \sim k (1 - p) e^{-(\theta x)^\alpha}$:

As $x \to 0^+$, approximate $e^{-(\theta x)^\alpha} \approx 1 - (\theta x)^\alpha$. Then:

For the PDF:

Since $S_X(x) \approx 1$:

□

The hazard rate, illustrated in Figure 3, transitions from $\frac{k \alpha \theta^{\alpha k}}{(1 - p)^k} x^{\alpha k - 1}$ near the origin to $\alpha \theta^\alpha x^{\alpha - 1}$ for large $x$. For $\alpha < 1$, it decreases, capturing early failures; for $\alpha > 1$, it increases, reflecting wear-out phases; and for $\alpha = 1$, it remains constant.

The log-survival function, depicted in Figure 4, confirms the WGSG's light-tailed nature, with $\log(S_X(x)) \approx \log(k (1 - p)) - (\theta x)^\alpha$ for large $x$. Smaller $\alpha$ values yield slower decay, suitable for prolonged survival, while larger $\alpha$ values produce steeper declines, indicating rapid failure onset.

The WGSG distribution's hazard rate supports complex patterns, such as bathtub-shaped hazard functions, common in reliability studies. The parameter $p$ enhances adaptability to varied tail behaviors, and expansions like $(1 - (1 - p) z)^k$ may arise in series representations, making the WGSG distribution a robust tool for precise tail and hazard estimates in reliability engineering and survival analysis.

3.1.1.   Limiting properties

The limiting behavior of the WGSG distribution reveals its flexibility and parameter sensitivity, with the survival function defined in Equation (3.1). As the shape parameter $\alpha \to 0^+$, the approximation $e^{-(\theta x)^\alpha} \to e^{-1}$ yields:

Such a constant survival function implies a degenerate, uniform-like distribution over $x$, rather than a point mass at $x = 0$. Conversely, as $\alpha \to \infty$:

The resulting step-function concentrates probability mass around $x = 1/\theta$. For the scale parameter, as $\theta \to 0^+$:

The distribution becomes improper, with probability mass dispersed toward infinity. As $\theta \to \infty$:

The distribution collapses to a point mass at $x = 0$. As the mixing parameter $p \to 0^+$:

The result is the exponentiated Weibull distribution, commonly applied in reliability modeling for independent failure modes. As $p \to 1^-$:

All mass concentrates at $x = 0$. For large $k$:

The distribution approaches a point mass at $x = 0$. As $k \to 0^+$:

Probability mass shifts toward infinity, producing an improper distribution. These limits often lead to degenerate distributions (e.g., point masses at $x = 0$ or $x = 1/\theta$, or improper distributions), highlighting parameter roles but offering limited practical utility. The exponentiated Weibull, obtained as $p \to 0^+$, is particularly valuable for reliability engineering. For practical applications, moderate parameter values, such as $\alpha \in \{0.3, 0.5, 1, 1.5, 2, 2.5\}$, $k \in \{1, 2, 3, 4, 5\}$, and $p = 0.3$, ensure flexible modeling without degeneracy, allowing precise control over tail behavior, as illustrated in Figures 3 and 4.

3.1.2.   Associations with other distributions

The WGSG distribution, parameterized by shape $\alpha > 0$, scale $\theta > 0$, weight $0 \leq p < 1$, and order $k \in \mathbb{Z}^+$, reduces to classical distributions used in reliability and survival analysis under specific parameter limits.

Exponentiated Weibull distribution: As $p \to 0^+$, the WGSG becomes the exponentiated Weibull distribution:

For $k = 1$, this is the standard Weibull distribution with shape $\alpha$ and scale $\theta^{-1}$:

For $k > 1$, the exponentiated Weibull models transformed failure times, extending the Weibull's flexibility. The parameter $p$ enriches hazard and tail behaviors.

Exponential distribution: Setting $\alpha = 1$, $p \to 0$, and $k = 1$, the WGSG simplifies to the exponential distribution:

This describes memoryless processes, such as inter-event times in a Poisson process, with constant hazard rate $\theta$. For $k > 1$, the survival function $1 - \left(1 - e^{-\theta x} \right)^k$ is non-standard, lacking memorylessness.

Rayleigh distribution: When $\alpha = 2$, $p \to 0$, and $k = 1$, the WGSG reduces to the Rayleigh distribution:

This models life data under accumulated stress, such as wear-related failures, aligning with the Weibull for $\alpha = 2$.

3.2.   Moment-generating function (MGF)

The moment-generating function (MGF) of a random variable $X$ is defined as:

Substituting this into the definition of $M_X(t)$ gives:

Using the binomial series expansion,

Substituting into $M_X(t)$:

Using the known integral formula:

where $E_{1/\alpha, 1}(z)$ is the Mittag-Leffler function defined as:

We set $c = i+j+1$ and obtain:

3.3.   $r$-th Moment

The $r$-th moment is given by:

Substituting the PDF:

Expanding the numerator and the denominator using the binomial series:

Thus,

Using the substitution $y = (i + j + 1)^{\frac{1}{\alpha}} \theta x$ and the gamma function identity:

we obtain:

Mean:

The first moment, or mean, is given by $E(X) = \mu_1$:

For $k = 1$, the formula simplifies:

This sum involves a generalized Hurwitz zeta function:

Variance:

The variance is given by:

Setting $r = 2$, we get:

For $k = 1$, this simplifies to:

This sum involves a generalized Hurwitz zeta function:

Thus, the variance $\text{Var}(X)$ is:

3.4.   Entropy

Entropy is a fundamental concept in information and communication theory, introduced by ^[36]. It quantifies the average missing information or uncertainty associated with a random source. Shannon defined the entropy of a continuous random variable $X$ as:

where $\log\left(\frac{1}{f_X(x)}\right)$ is referred to as the uncertainty associated with the outcome $x$. The probability distribution $f_X(x)$ fully describes the probabilistic characteristics of a random variable. However, when comparing two or more probability distributions, entropy provides a quantitative measure to compare the randomness of different distributions.

For the new family of distributions, the entropy is given by:

Expanding the logarithm term using its properties:

Taking expectations, the entropy simplifies to:

where $\mathbb{E}[\log X]$ captures the spread of the distribution and $\mathbb{E}[(\theta X)^\alpha]$ relates to the moment of $X$. The last two expectation terms capture the influence of the geometric component and tail behavior.

The partial derivative with respect to the parameters are given by:

The residual entropy, also known as conditional entropy in this context, quantifies the uncertainty in a random variable $X_{(k)}$ given that it exceeds a threshold $t \geq 0$. For the WGSG distribution, this corresponds to the differential entropy of the residual lifetime distribution $X_{(k)} - t \mid X_{(k)} > t$. The residual entropy is defined as:

where

Applying the logarithmic property $\log \left(\frac{a}{b} \right) = \log a - \log b$, the expression simplifies as follows:

Since $\log S_{X_{(k)}}(t)$ is constant with respect to $x$, the second integral becomes:

Noting that $\int_{t}^{\infty} f_{X_{(k)}}(x) \, dx = S_{X_{(k)}}(t)$, we have:

Thus:

This form highlights two components: the entropy integral over the tail of the distribution (from $t$ to $\infty$) and a correction term involving the survival probability at time $t$.

For $t = 0$, the survival function is:

Thus:

which is the differential entropy of $X_{(k)}$, often referred to as the Shannon entropy in this context:

This measure is valuable in reliability engineering and information theory, as it quantifies the uncertainty in a system's remaining lifetime after surviving up to time $t$, aiding in reliability analysis and decision-making under uncertainty.

4.   Parameter estimation

4.1.   Maximum likelihood estimation

Let $x = (x_1, x_2, \ldots, x_n)$ be a random sample drawn from the WGSG distribution with the unknown parameter vector $\varphi = (\alpha, \theta, p)$. Using the PDF of the WGSG distribution, the log-likelihood function is given by:

Substituting $f_X(x_i; \varphi)$ into the log-likelihood function, we obtain:

Thus,

The partial derivatives with respect to $\alpha$, $\theta$, and $p$ are:

These nonlinear equations are typically solved numerically using iterative methods such as the Newton-Raphson method, the expectation-maximization (EM) algorithm, and gradient-based optimization (e.g., BFGS).

4.2.   EM algorithm

Let the observed data be $\{x_i\}_{i = 1}^n$ and the latent variable $N$, which follows a truncated geometric distribution. The joint PDF is:

Thus, the complete data log-likelihood function is:

E-step (Expectation step): In the E-step, we calculate the expected value of the latent variable $N$ given the observed data $x_i$. The posterior distribution of $N$ given $x_i$ is:

Thus, the expected value of $N$ given $x_i$ is:

Substituting for the posterior:

Simplifying, we get the following expression:

M-step (Maximization step): In the M-step, we maximize the expected complete-data log-likelihood with respect to the parameters $\alpha$, $\theta$, and $p$. We calculate the updates for each parameter.

where

and

4.3.   Method of moments estimation

The method of moments (MM) provides consistent estimators for the parameters $\alpha > 0$, $\theta > 0$, and $0 \leq p < 1$ of the WGSG distribution by equating sample moments to their theoretical counterparts. For a fixed-order $k \in \mathbb{Z}^+$, the estimation leverages the first three moments to solve a nonlinear system.

The $r$-th raw moment of the WGSG distribution is:

where:

which equals $\mathbb{E}\left[(I + J + 1)^{-\left(\frac{r}{\alpha} + 1\right)}\right]$, with $I \sim \text{NB}(k, 1 - p)$ (negative binomial, failures before $k$ successes, success probability $1 - p$) and $J \sim \text{Bin}(k - 1, 1/2)$ independent. The series converges absolutely due to the rapid decay of $(i + j + 1)^{-\left(\frac{r}{\alpha} + 1\right)}$.

For a sample of size $n$, the $r$-th sample moment is $\hat{\mu}_r = n^{-1} \sum_{i = 1}^n x_i^r$. Matching moments for $r = 1, 2, 3$ yields:

where $M_r(\alpha, p, k) = \Gamma\left(\frac{r}{\alpha} + 1\right) S_r(\alpha, p, k)$. We solve for $\theta$:

Substituting into Eq (4.2):

From Eqs (4.4) and (4.3):

and equating with (4.5):

This nonlinear equation can be solved numerically using a root-finding algorithm, such as the Newton-Raphson method or a bracketing method (e.g., bisection), implemented in software like $\texttt{R}$ or $\texttt{Python}$. For instance, in $\texttt{R}$, the $\texttt{nleqslv}$ package can be used by defining the following system of equations:

and iteratively adjusting $\alpha$ and $p$ to minimize the functions $h_1$ and $h_2$. Appropriate initial guesses (e.g., $\alpha = 1$, $p = 0.5$) and parameter bounds ($\alpha > 0$, $0 \leq p < 1$) help ensure convergence. The function $S_r(\alpha, p, k)$ is approximated by truncating the infinite series at a sufficiently large upper index $i_{\text{max}} \approx 1000$, where additional terms become negligible. Once $\alpha$ is estimated, the remaining parameters $\theta$ and $p$ are computed from Eqs (4.5) and (4.6), respectively.

4.3.1.   Bayesian estimation

Bayesian estimation provides a robust framework for estimating the parameters $\alpha > 0$, $\theta > 0$, and $0 \leq p < 1$ of the WGSG distribution, leveraging prior knowledge to yield full posterior distributions, while treating the fixed parameter $k \in \mathbb{N}$ as a known value that defines a distinct distribution for each integer. For a random sample $\mathbf{x} = (x_1, x_2, \ldots, x_n)$ from the WGSG distribution with the PDF given in Eq (2.1), the likelihood function is:

Independent priors are assigned to reflect prior beliefs: $\alpha \sim \text{Gamma}(a_\alpha, b_\alpha)$, $\theta \sim \text{Gamma}(a_\theta, b_\theta)$, and $p \sim \text{Beta}(a_p, b_p)$, chosen for their domain compatibility and computational tractability, with hyperparameters $(a_\alpha, b_\alpha)$, $(a_\theta, b_\theta)$, and $(a_p, b_p)$ set to weakly informative values (e.g., $a_\alpha = b_\alpha = 0.001$) to approximate noninformative priors when prior knowledge is limited. The posterior distribution is defined as:

The posterior expected loss is:

The joint posterior density, proportional to the likelihood times the priors, is:

but its complexity, driven by the WGSG's nested exponential and polynomial terms, prevents a closed-form solution, necessitating Markov Chain Monte Carlo (MCMC) methods like the Metropolis-Hastings algorithm or Hamiltonian Monte Carlo, implemented via tools such as $\texttt{Stan}$ or $\texttt{JAGS}$, to approximate posterior summaries. Under mean square error (MSE) loss, the Bayes estimator for a parameter $\varphi \in \{\alpha, \theta, p\}$, denoted $\hat{\varphi}_{\text{Bayes}}$, is the posterior mean:

which minimizes the expected posterior loss and is approximated using MCMC samples, with analogous expressions for $\hat{\alpha}_{\text{Bayes}}$, $\hat{\theta}_{\text{Bayes}}$, and $\hat{p}_{\text{Bayes}}$.

5.   Simulation study

To evaluate the performance of the maximum likelihood (ML), method of moments (MM), and Bayesian estimation methods for the WGSG distribution, we conducted a simulation study. We generated $1000$ random samples of sizes $n = 20$, 50, and 100 from the WGSG distribution, considering three parameter settings for $\varphi = (\alpha, \theta, p)$: $(0.5, 1, 0.5)$, $(1, 1, 0.5)$, and $(2, 1, 0.5)$, combined with $k \in \{1, 2, 3\}$. The simulations were implemented using the R programming language.

For each sample, we estimated the parameters $\widehat{\varphi} = (\widehat{\alpha}, \widehat{\theta}, \widehat{p})$ and computed the bias and mean squared error (MSE) of the estimators. The bias and MSE are defined as:

where $\text{Var}(\widehat{\varphi})$ is the variance of the estimator. An estimator $\widehat{\varphi}_1$ is considered more efficient than $\widehat{\varphi}_2$ if $\text{MSE}(\widehat{\varphi}_1) < \text{MSE}(\widehat{\varphi}_2)$.

The results are summarized in Table 1, which reports the bias and mean squared error (MSE) of the parameter estimates for the ML, MM, and Bayesian methods across all scenarios. The findings reveal that for each order statistic $k \in \{1, 2, 3\}$, the MSEs decrease as the sample size $n$ increases, indicating improved estimation precision with larger samples. Moreover, the Bayesian method consistently exhibits lower MSEs compared to the ML and MM approaches, suggesting greater efficiency under these simulation settings. As the sample size increases to $n = 100$, the bias of the ML estimators tends to decrease, aligning with their known property of asymptotic unbiasedness. However, the MSE results suggest that ML estimation may struggle in small samples when applied to the WGSG distribution, primarily due to the sensitivity of numerical optimization routines to initial values and the complexity of the likelihood surface. In contrast, the Bayesian approach employs Metropolis-Hastings MCMC sampling to generate posterior samples, with parameters estimated using posterior summaries such as the mean or median. While the Bayesian estimators may be influenced by the choice of prior and exhibit some bias, they often achieve lower MSE by effectively regularizing parameter estimates, especially in small-sample scenarios where likelihood-based estimates may be unstable. The Bayesian MSE for $\widehat{p}$ is notably low (e.g., $0.004$ at $n = 20, k = 2, \alpha = 2$) compared to ML ($0.142$) and MM ($0.038$), underlining its superior efficiency. Similarly, for $\widehat{\alpha}$, the Bayesian bias is substantially lower (e.g., $-0.592$ at $n = 20, k = 3, \alpha = 2$) than the corresponding ML bias (1.957), demonstrating the stabilizing effect of the prior. Furthermore, Bayesian estimation shows greater robustness to changes in $k$ and $\alpha$: for example, the MSE for $\widehat{\alpha}$ increases only modestly from $0.021$ at $k = 1, \alpha = 0.5$ to $0.414$ at $k = 3, \alpha = 2$, whereas ML and MM estimators exhibit more pronounced increases.

The MM approach estimates parameters by equating the theoretical moments of the WGSG distribution to their sample counterparts and solving for $\alpha$, $\theta$, and $p$. This method is computationally simpler than ML and Bayesian approaches, as it avoids complex likelihood evaluations and MCMC sampling. However, it tends to be less reliable for distributions with intricate moment structures such as the WGSG, which are sensitive to sample variability and outliers. The simulation results reflect this limitation: the MM approach consistently yields the highest bias and MSE across all parameters. For example, at $n = 20, k = 3, \alpha = 0.5$, the MM bias for $\widehat{\alpha}$ is $1.490$, and the MSE is $2.863$—substantially worse than ML (bias: $0.487$, MSE: $0.323$) and Bayesian (bias: $-0.067$, MSE: $0.014$) results in the same setting. Even at larger sample sizes ($n = 100$), the MM estimator remains biased (e.g., $1.264$ at $k = 3, \alpha = 0.5$), suggesting limited improvement with increased data. This underperformance likely stems from the difficulty of accurately estimating higher-order moments in the WGSG distribution. Although ML estimation offers favorable large-sample properties, its finite-sample performance can be inferior to the Bayesian approach. The Bayesian method provides more stable and accurate estimates across various parameter settings and sample sizes. The MM approach, while computationally attractive, proves to be the least reliable among the three, particularly in small-sample contexts and in the presence of complex moment behavior.

6.   Applications

In this section, we compare the adequacy of the WGSG distribution for 5 special cases ($k = 1, 2, 3, 4, 5$) compared with the Weibull, Weibull-geometric (WG), and Weibull-Poisson (WP) distributions using two real datasets. The first data, sourced from ^[26], represent the fatigue life (rounded to the nearest thousand cycles) for $67$ specimens of Alloy T7987 that failed before accumulating $300$ thousand cycles. The dataset is as follows:

(94,118,139,159,171,189,227, 96,121,140,159,172,190,256, 99,121,141,159,173,196,257, 99,123,141,159,176,197,269,104,129,143,162,177,203,271,108,131,144,168,180,205,274,112,133,149,168,180,211,291,114,135,149,169,184,213,117,136,152,170,187,224,117,139,153,170,188,226).

The second dataset represents the time intervals (in hours) between successive failures of the air conditioning systems in a fleet of seven Boeing 720 airplanes, specifically aircraft numbers 7910, 7911, 7912, 7913, 7914, 7915, and 7916. This subset, consisting of 125 observations, was analyzed by ^[37]. It is part of a larger dataset of 213 observations, originally studied by ^[38] and further explored by ^{[10,12,39,40]}. The 125 observations for the specified aircraft are as follows:

(74, 57, 48, 29,502, 12, 70, 21, 29,386, 59, 27,153, 26,326, 55,320, 56,104,220,239, 47,246,176,182, 33, 15,104, 35, 23,261, 87, 7,120, 14, 62, 47,225, 71,246, 21, 42, 20, 5, 12,120, 11, 3, 14, 71, 11, 14, 11, 16, 90, 1, 16, 52, 95, 97, 51, 11, 4,141, 18,142, 68, 77, 80, 1, 16,106,206, 82, 54, 31,216, 46,111, 39, 63, 18,191, 18,163, 24, 50, 44,102, 72, 22, 39, 3, 15,197,188, 79, 88, 46, 5, 5, 36, 22,139,210, 97, 30, 23, 13, 14,359, 9, 12,270,603, 3,104, 2,438, 50,254, 5,283, 35, 12).

Weibull distribution:

Weibull-geometric (WG) distribution:

Weibull-Poisson (WP) distribution:

Table 2 presents the estimated parameters for various lifetime models fitted to the fatigue life data and air conditioning system failure data. The WGSG distribution, evaluated across different order statistics $k$, demonstrates flexibility in capturing varying data characteristics through changes in its parameters. For the numerical computation of maximum likelihood estimation, we use the L-BFGS-B method, as proposed by ^[41]. This method is a limited-memory variation of the BFGS quasi-Newton approach, applied when solving a system of simultaneous equations. The maximum likelihood estimates (MLEs) of the parameters, along with the $p$-values of the K-S statistics for these models, are summarized in Table 3. The results also include the AIC, BIC, and Shannon's entropy as model selection criteria. It is important to note that the Kolmogorov-Smirnov (K-S) test compares an empirical distribution with a known (non-estimated) distribution. This test helps determine if a sample originates from a population that follows a specific distribution (H$_{0}$: the data follow a specified distribution).

The analysis of the fatigue life data reveals that the WGSG distribution with $k = 4$ provides the best fit among all tested models, outperforming the Weibull, Weibull-geometric (WG), and Weibull-Poisson (WP) distributions across all key metrics. Specifically, the WGSG ($k = 4$) model achieves the highest log-likelihood ($-347.367$), indicating the best overall fit to the data. It also records the lowest Akaike Information Criterion (AIC = $700.735$) and Bayesian Information Criterion (BIC = $707.349$), which balance goodness-of-fit with model complexity, reinforcing its superiority. Furthermore, it exhibits the lowest entropy ($5.185$), suggesting it offers the most precise representation of the data's underlying distribution. The Kolmogorov-Smirnov (K-S) test supports the adequacy of all models, with $p$-values exceeding $0.05$ (WGSG $k = 4$: $0.941$), indicating that none can be rejected based on fit to the empirical distribution.

Notably, the WG and WGSG ($k = 1$) models, which are mathematically equivalent, achieve the highest K-S $p$-value ($0.988$) and perform exceptionally well, with a log-likelihood of $-348.552$, AIC of $703.104$, and BIC of $709.718$. This excellent fit reflects the WGSG's reduction to the WG distribution when $k = 1$. However, the WGSG ($k = 4$) edges out these models by further optimizing the AIC, BIC, and entropy, suggesting it captures additional nuances in the data. In contrast, the standard Weibull model performs the worst, with the lowest log-likelihood ($-353.292$), highest AIC ($710.584$), highest BIC ($714.993$), and highest entropy ($5.273$). This indicates that the Weibull distribution is too simplistic to fully model the complexity of fatigue life data.

For the air conditioning system failures dataset, the WGSG distribution with $k = 2$ emerges as the best-fitting model, surpassing the Weibull, WG, and WP distributions across most metrics. It achieves the highest log-likelihood ($-685.279$), lowest AIC ($1376.557$), and lowest BIC ($1385.042$), indicating an optimal balance of fit and complexity. Additionally, it has the lowest entropy ($5.482$), suggesting a precise representation of the failure data. The K-S test yields a high p-value ($0.956$), confirming the model's adequacy and inability to be rejected.

The WG and WGSG ($k = 1$) models, again identical, perform strongly with a log-likelihood of $-686.117$, AIC of $1378.233$, BIC of $1386.718$, and a K-S $p$-value of $0.946$, closely trailing the WGSG ($k = 2$). Meanwhile, the WGSG ($k = 4$) achieves the highest K-S $p$-value ($0.976$), indicating an excellent fit to the empirical distribution, though its AIC ($1379.131$) and BIC ($1387.615$) are slightly higher than those of WGSG ($k = 2$). The standard Weibull model again performs the worst, with a log-likelihood of $-687.632$, AIC of $1379.265$, BIC of $1384.921$, and entropy of $5.501$, underscoring its inadequacy for capturing the complexity of air conditioning failure times. These results suggest that the WGSG model is particularly well-suited for modeling air conditioning system failures, offering practical value for maintenance scheduling and system reliability predictions. As with the fatigue life data, the increased complexity of the WGSG distribution requires careful validation to ensure generalizability.

To explore the empirical behaviors that the failure rate function may exhibit, we use the graphical method based on the total time on test (TTT) plot, introduced by ^[42]. In its empirical form, the TTT plot is constructed using the values $r|n$ and $G(r|n)$, where:

Here, $r = 1, \dots, n$, and $X_{(i:n)}$ denotes the order statistics of the sample. The TTT plot can take different forms, resembling the Gini index, and serves as a rough indicator of the shape of the failure rate function. Specifically, when the curve approaches a straight diagonal, a constant failure rate is appropriate, suggesting the data follow an exponential distribution. A convex curve indicates the data are likely from a decreasing failure rate (DFR) distribution.

The empirical and theoretical distribution fits are visually evaluated using Figures 5 and 6. Figure 5 compares the empirical PDFs and CDFs with their theoretical counterparts for both the fatigue life and air conditioning system failures datasets. For the fatigue life data, the WGSG distribution with $k = 4$ closely matches the empirical PDF and CDF, reinforcing the numerical evidence of its superior performance. Similarly, for the air conditioning system failures data, the WGSG ($k = 2$) provides a precise fit to the empirical distributions, particularly in the PDF, where it captures the failure time distribution's shape more accurately than the Weibull, WG, or WP distributions.

Figure 6 presents probability-probability (PP) and TTT plots for all distributions. The PP plots for both datasets show that the WGSG distribution lies closest to the diagonal line, indicating strong agreement between the theoretical and empirical distributions. The TTT plots offer additional insight into the failure rate behavior, which is explored further in the next section. For the fatigue life data, the TTT plot exhibits a convex shape, suggesting a decreasing failure rate (DFR), which the WGSG ($k = 4$) effectively captures, as evidenced by its superior numerical and visual fit. Likewise, the TTT plot for the air conditioning system failures data displays a convex shape, indicating a DFR, which the WGSG ($k = 2$) accurately reflects through its strong fit metrics and close alignment in the PP plot. These findings underscore the WGSG distribution's suitability for both datasets, particularly its ability to accommodate decreasing failure rate behaviors.

Across both datasets, the WGSG distribution consistently outperforms simpler alternatives like the Weibull, WG, and WP distributions, with the optimal $k$ value varying by dataset ($k = 4$ for fatigue life data and $k = 2$ for air conditioning system failures). This adaptability highlights the WGSG's capacity to balance fit and complexity, making it more effective than traditional distributions for diverse data structures. While the Weibull distribution remains statistically adequate (K-S $p$-values $> 0.05$), it consistently underperforms, suggesting it may be too rigid for reliability and survival data with heterogeneous failure patterns. The WGSG's superior performance, supported by both numerical metrics and graphical analysis, has significant implications for lifetime and reliability engineering, where precise failure time analysis can enhance design, maintenance, and safety practices.

7.   Conclusions

The WGSG distribution offers a versatile framework for lifetime data analysis, integrating Weibull and shifted geometric properties to model diverse failure patterns with precision. Its unique focus on the $k$-th order statistic enhances flexibility beyond models like the Weibull-geometric, capturing varying hazard rates and tail behaviors, as detailed in theoretical analyses of entropy and asymptotic properties. Simulations validate robust parameter estimation via maximum likelihood estimation, Expectation-Maximization, Bayesian methods, and the method of moments, while applications to real datasets demonstrate superior fit compared to Weibull, Weibull-Poisson, and Weibull-geometric distributions. This adaptability positions the WGSG as a valuable tool for reliability and risk assessment. Future work will extend it to censored data, multivariate structures, and alternative methods like minimum distance estimation and generalized method of moments, broadening its statistical impact.

Author contributions

M. Rahmouni: Conceptualization, visualization, methodology, writing—original draft, writing—review and editing, formal analysis, validation, investigation, software; D. Ziedan: visualization, data curation, writing—review, investigation. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Funding

This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Grant No. KFU251549].

Conflict of interest

The authors declare no conflicts of interest in this paper.