A New Exponentiated Power Distribution for Modeling Censored Data with Applications to Clinical and Reliability Studies

1. Introduction

The study of lifetime distributions plays a crucial role in reliability analysis and survival modeling. [1] pioneered a notable approach by developing a one-parameter lifetime distribution based on a two-component mixture: an exponential distribution and a gamma distribution (with a shape parameter of 2), both sharing a common scale parameter $\theta$, and mixed with a proportion of $p={\displaystyle \frac{{\theta }^2}{{\theta }^2+1}}$. This methodology was inspired by [2]’s foundational work on mixing proportions. Distributions derived from this innovative method often exhibit a desirable bathtub-shaped hazard rate, making them highly suitable for various real-life reliability applications. Their appeal lies in their simplicity, consolidating properties of both exponential and gamma distributions into a single-parameter framework. This simplicity has subsequently spurred numerous investigations, leading to the development of other mixture-based distributions, including the Chris-Jerry distribution by [3], the Odoma distribution by [4], the Pranav distribution by [5], the Aradhana distribution by [6], and the Fav-Jerry distribution by [7], among others. To enhance the flexibility of baseline distributions, particularly in capturing complex data characteristics such as skewness, kurtosis, and varied goodness-of-fit, power distributions have emerged as a significant advancement. This approach introduces an additional shape parameter, which invariably improves the adaptability of the resulting distribution. For example, the power Lindley distribution, proposed by [8], remarkably captures diverse hazard rate shapes (increasing, decreasing, and bathtub), a capability the original one-parameter Lindley distribution lacks. Similarly, [9] introduced the two-parameter power Shanker distribution, adept at modeling increasing, decreasing, and constant hazard rates. These power parameters critically modify the shape of the distribution by influencing skewness and tail behavior, enabling a better fit for skewed or heavy/light-tailed data, and facilitating the modeling of non-monotonic hazard rates often observed in survival and reliability data. Such enhancements lead to improved goodness-of-fit over base models, as widely supported by empirical applications using real-life datasets. The increasing complexity of real-world data, as summarized in Table (Table 1) along with standard limitations and advanced solutions, necessitates the development of sophisticated distributions capable of accurately capturing and explaining these intricate behaviors.

Beyond power transformations, another potent method for increasing distributional flexibility is through the exponentiation approach. This technique, often involving the addition of an extra shape parameter, is well-regarded for its capacity to modify the tail behavior and further enrich the variety of hazard rate shapes a distribution can model. Applying this exponentiation method to the Power Shanker distribution yields the Exponentiated Power Shanker (EPS) distribution. The EPS distribution consequently inherits enhanced flexibility from both the power and exponentiation processes, making it highly adaptable for modeling diverse data patterns and achieving superior goodness-of-fit. Furthermore, for practical applications in regression analysis, particularly with positive-valued data such as survival times, it is common to apply a logarithm transformation to the random variable. This transformation allows a linear predictor, dependent on covariates, to directly influence the scale or location parameters of the distribution, leading to a more interpretable and tractable regression framework. Numerous studies underscore the critical role of clinical and demographic factors, such as age and various comorbidities (e.g., cardiovascular disease, asthma, diabetes, neurological disorders, and obesity), as significant predictors of survival times [28,29,30,31,32]. These variables, commonly treated as covariates in statistical models, are essential for enhancing the interpretability of regression-based survival analyses in biomedical contexts. The inherent complexities of such data, including individual heterogeneity and the frequent presence of censoring, necessitate highly flexible and robust modeling approaches. To address these challenges, we propose the novel Exponentiated Power Shanker (EPS) distribution and its regression counterpart, the Log-Exponentiated Power Shanker (LEPS) regression model. By integrating the LEPS within a regression framework, our model offers an enhanced capability to analyze the lifetimes of COVID-19 patients in Brazil, particularly under the prevalent censoring conditions found in medical data. This robust framework allows for the accurate capture of covariate effects, such as age, heart disease, asthma, diabetic condition, neurological condition, and obesity, on patient lifetimes, thereby facilitating a deeper understanding of disease dynamics.

The primary objectives of this study encompass the introduction of the Exponentiated Power Shanker (EPS) distribution, a three-parameter extension of the Shanker distribution, designed to enhance modeling flexibility for diverse data behaviors, including right-skewed and heavy-tailed phenomena prevalent in reliability and survival analysis. This involves a systematic derivation and examination of its core theoretical properties, such as complete and incomplete moments, entropy, and moment-generating function. Furthermore, the study aims to apply Maximum Likelihood Estimation for parameter estimation of the EPS distribution, assessing its performance through a Monte Carlo simulation across various sample sizes and parameter values. The research also extends the EPS framework into a regression model, the Log-Exponentiated Power Shanker (LEPS) regression, to incorporate covariate data, with parameter estimation also performed via MLE, thereby enabling the modeling of covariate influence on the distribution’s parameters. Ultimately, the study showcases the applicability and superior goodness-of-fit of the proposed distributions using real-world datasets. Specifically, the EPS distribution is fitted to failure times of item data and mortality rates of HIV/AIDS patients in Germany. The LEPS regression model is then applied to analyze the lifetimes (in days) of COVID-19 patients in Brazil, correlating these lifetimes with clinical covariates including age, heart disease, asthma, diabetic condition, neurological condition, and obesity. A comprehensive evaluation, encompassing comparisons against existing distributions using various goodness-of-fit statistics and graphical diagnostics, confirms the models’ performance.

The novelty of this work resides in several key aspects. The proposed Exponentiated Power Shanker (EPS) distribution innovatively combines power transformation and exponentiation approaches, yielding a three-parameter distribution significantly more adaptable than its predecessors. This dual enhancement introduces an additional shape parameter, enabling it to capture complex data characteristics such as varying skewness, kurtosis, and a wider array of hazard rate shapes, including bathtub, L-shape, increasing, and decreasing, which simpler distributions lack. The EPS distribution is particularly designed to effectively model complex data behaviors like right-skewed and heavy-tailed phenomena that often challenge less flexible distributions. The development and specific application of the Log-Exponentiated Power Shanker (LEPS) regression model are novel, especially for analyzing positive-valued data like lifetimes in a regression context. By applying a logarithm transformation, the model facilitates linear predictors based on covariates to influence the distribution’s parameters, offering a more interpretable framework for biomedical and reliability data with censoring. The study’s practical utility is robustly demonstrated through diverse real-world applications: fitting the EPS distribution to failure times of item data and mortality rates of HIV/AIDS patients in Germany, and utilizing the LEPS regression model for analyzing COVID-19 patient lifetimes in Brazil against relevant clinical covariates. This showcases its broad applicability and effectiveness in challenging data environments. Finally, the work includes a rigorous evaluation of the proposed distributions, comprising formal derivations of structural properties, a Monte Carlo simulation for parameter estimation, and extensive comparisons with several existing distributions using multiple goodness-of-fit criteria and graphical diagnostics, collectively affirming its superior modeling capacity.

Following the introduction, Section 2 is dedicated to a detailed construction of the EPS distribution itself. With this as a foundation, Section 3 examines its structural properties. It subsequently focuses on parameter estimation in Section 4, where the Maximum Likelihood Estimation (MLE) methodology is outlined. In Section 5, the discussion extends to the application of the EPS distribution under a regression framework. For verifying the theoretical derivations and estimation techniques, a comprehensive simulation study is presented in Section 6. Then Section 7 addresses real-world applications, specifically demonstrating the application of the distribution with the fitting of Log-EPS Regression in Section 7.1 and its application on reliability data in Section 7.2. Finally, Section 8 summarises the paper with an overview of the most significant findings and potential future research directions.

2. The Exponentiated Power Shanker (EPS) Distribution

The cumulative distribution function (CDF) of Power Shanker distribution due to [9] is given as follows

$$G(x;\theta ,\alpha )=1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }};x>0,$$

(1)

where $\theta$, $\alpha >0$ are the scale and shape parameters respectively. The probability density function (PDF) corresponding to equation (1) is

$$g(x;\theta ,\alpha )={\displaystyle \frac{\alpha {\theta }^2}{({\theta }^2+1)}}{x}^{\alpha -1}(\theta +x^{\alpha })e^{-\theta x^{\alpha }}.$$

(2)

[33] and [34] developed the exponentiated family of distributions. By definition, a random variable X is said to assume exponentiated-G distribution if its CDF and PDF are respectively

$$F(x;c,\xi )=G{(x;\xi )}^c,$$

(3)

and

$$f(x;c,\xi )=cg(x;\xi )G{(x;\xi )}^{c-1},$$

(4)

where $\xi$ is the parameter vector of $G(.)$. By substituting Equation (1) into Equation (3) yields Equation (5). Similarly, substituting Equations (1) and (2) into Equation (4) yields Equation (6). Hence, the CDF and PDF of the $X\sim$ EPS $(c,\theta ,\alpha )$ are expressed as follows;

$$F(x;c,\theta ,\alpha )={\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^c,x>0,c,\alpha ,\theta >0,$$

(5)

and

$$f(x;c,\theta ,\alpha )={\displaystyle \frac{c\alpha {\theta }^2}{{\theta }^2+1}}(\theta +x^{\alpha })x^{\alpha -1}{e}^{-\theta x^{\alpha }}{\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}.$$

(6)

The hazard rate function (HRF) is given as

$$h(x;c,\theta ,\alpha )={\displaystyle \frac{\frac{c\alpha {\theta }^2}{{\theta }^2+1}(\theta +x^{\alpha })x^{\alpha -1}{e}^{-\theta x^{\alpha }}{\left\{1-\left[1+\frac{\theta x^{\alpha }}{{\theta }^2+1}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}}{1-{\left\{1-\left[1+\frac{\theta x^{\alpha }}{{\theta }^2+1}\right]e^{-\theta x^{\alpha }}\right\}}^c}}.$$

To verify the normalization condition, we evaluate

$$I={\int }_0^{\infty }{\displaystyle \frac{c\alpha {\theta }^2}{{\theta }^2+1}}(\theta +x^{\alpha })x^{\alpha -1}{e}^{-\theta x^{\alpha }}{\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}dx.$$

Substituting $t=x^{\alpha }$, we get $dt=\alpha x^{\alpha -1}dx$, yielding:

$$I={\int }_0^{\infty }{\displaystyle \frac{c{\theta }^2}{{\theta }^2+1}}(\theta +t)e^{-\theta t}{\left\{1-\left[1+{\displaystyle \frac{\theta t}{{\theta }^2+1}}\right]e^{-\theta t}\right\}}^{c-1}dt.$$

Now, expand the bracketed term using the binomial expansion:

$${\left(1-\left[1+{\displaystyle \frac{\theta t}{{\theta }^2+1}}\right]e^{-\theta t}\right)}^{c-1}=\sum _{k=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{k}}\right){(-1)}^k{\left(1+{\displaystyle \frac{\theta t}{{\theta }^2+1}}\right)}^k{e}^{-k\theta t}.$$

Rewriting the integral using the summation, we have:

$$I={\displaystyle \frac{c{\theta }^2}{{\theta }^2+1}}\sum _{k=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{k}}\right){(-1)}^k{\int }_0^{\infty }(\theta +t)e^{-(k+1)\theta t}{\left(1+{\displaystyle \frac{\theta t}{{\theta }^2+1}}\right)}^{k}dt.$$

Each term inside the summation is a gamma-like integral. Using the gamma function identity,

$${\displaystyle \frac{1}{{\lambda }^{m+1}}}{\int }_0^{\infty }{t}^m{e}^{-\lambda t}dt={\displaystyle \frac{\mathsf{\Gamma }(m+1)}{{\lambda }^{m+1}}},$$

we find that the integral simplifies to 1. This confirms that $f(x;c,\theta ,\alpha )$ is properly normalized.

Figure 1 and Figure 2 collectively illustrate the substantial flexibility of the PDF for the EPS distribution, parameterized by c, $\theta$, and $\alpha$. Figure 1 showcases the capacity of the EPS distribution to exhibit various unimodal shapes, ranging from markedly right-skewed (as seen with the black, red, blue, and orange curves) to distinctly left-skewed (exemplified by the green dashed-dotted curve). These plots demonstrate that varying the parameters allows the mode of the distribution to shift across the range of x values, accommodating different central tendencies and degrees of asymmetry in the data. Figure 2 further highlights the versatility of the EPS distribution, particularly its ability to model shapes commonly encountered in survival and reliability analysis. The black and red dashed curves, with their sharp initial peaks near zero followed by a decline, indicate a decreasing PDF shape, which is characteristic of phenomena where events are more likely to occur early in the observation period. Most notably, the green dashed-dotted curve in Figure (Figure 2 demonstrates a complex, irregular shape, suggesting that the EPS distribution is capable of capturing more intricate data patterns beyond simple unimodal or monotonically decreasing forms. The orange dashed-dotted curve in this figure provides another example of a unimodal, right-skewed shape. Together, both figures affirm that by adjusting its three parameters, the EPS distribution can adeptly fit a wide array of empirical data distributions, including those with varying skewness, different peak locations, and even non-standard or decreasing probability densities, making it a robust tool for diverse modeling applications.

Figure 3 and Figure 4 graphically depict the versatility of the EPS distribution in modeling various shapes of the HRF, which is crucial for understanding instantaneous failure rates in survival and reliability analysis. Figure 3 demonstrates the EPS distribution’s capacity to exhibit the widely recognized bathtub shape, where the hazard initially decreases, remains relatively constant, and then increases. This figure also illustrates the ability to capture continuously increasing and decreasing hazard rate patterns depending on the chosen parameter combinations of c, $\theta$, and $\alpha$. Complementing Figure 3, Figure 4 further accentuates the diverse HRF shapes achievable by the EPS distribution. It reinforces the presence of prominent bathtub shapes, along with clear examples of increasing and decreasing hazard functions. Notably, some parameter settings in Figure 4 demonstrate scenarios where the hazard rate can be nearly constant or exhibit sharp, pronounced decreases, indicating early life failures. The broad range of HRF configurations across both figures underscores the inherent flexibility of the EPS distribution, making it highly adaptable for modeling complex time-to-event data that might exhibit different risk characteristics over time, from early life reliability issues to wear-out phases or simply monotonically changing risk profiles.

2.1. Linear Representation

Using the generalized binomial theorem on ${\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}$, the pdf in Equation (6) can be expressed as

$$\begin{array}{cc}\hfill f(x;c,\alpha ,\theta )& ={\displaystyle \frac{c\alpha }{{\theta }^2+1}}\left(\theta +x^{\alpha }\right)\sum _{j=0}^{\infty }\sum _{i=j}^{\infty }\sum _{k=0}^j{(-1)}^{i+k}\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{j+k-1}{k}}\right)\hfill \\ \hfill & \times {\theta }^{2-j-2k}{x}^{\alpha (j+1)-1}{e}^{-(i+1)\theta x^{\alpha }}.\hfill \end{array}$$

3. Structural Characteristics

The $r^{th}$ crude moment of $X\sim$ EPS $(c,\alpha ,\theta )$ is defined as

$${\mu }_r^{\prime }={\displaystyle \frac{c\alpha {\theta }^2}{{\theta }^2+1}}{\int }_0^{\infty }(\theta +x^{\alpha })x^{r+\alpha -1}{e}^{-\theta x^{\alpha }}{\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}dx.$$

Applying binomial theorem on ${\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}$, we obtain

$${\mu }_r^{\prime }={\displaystyle \frac{c\alpha {\theta }^2}{{\theta }^2+1}}\sum _{i=0}^{\infty }{(-1)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right){\int }_0^{\infty }(\theta +x^{\alpha })x^{r+\alpha -1}{e}^{-(1+i)\theta x^{\alpha }}{\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]}^{i}dx.$$

(7)

We further decompose ${\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]}^i$ using binomial expansion to have Equation (7) reduce to

$${\mu }_r^{\prime }=c\alpha \sum _{i=0}^{\infty }\sum _{j=0}^i{(-1)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right){\displaystyle \frac{{\theta }^{j+2}}{{({\theta }^2+1)}^{j+1}}}{\int }_0^{\infty }(\theta +x^{\alpha })x^{\alpha (1+j)+r-1}{e}^{-(i+1)\theta x^{\alpha }}dx.$$

Intuitively, we replace ${\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^i}$ by ${\displaystyle \sum _{j=0}^{\infty }}{\displaystyle \sum _{i=j}^{\infty }}$ and complete the integration as follows;

$$\begin{array}{cc}\hfill {\mu }_r^{\prime }& =c\sum _{j=0}^{\infty }\sum _{i=j}^{\infty }{(-1)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right){\displaystyle \frac{{\theta }^{j+2}}{{({\theta }^2+1)}^{j+1}}}\left\{{\displaystyle \frac{\theta \mathsf{\Gamma }\left[1+j+\frac{r}{\alpha }\right]}{{\left[(1+j)\theta \right]}^{1+j+\frac{r}{\alpha }}}}+{\displaystyle \frac{\mathsf{\Gamma }\left[2+j+\frac{r}{\alpha }\right]}{{\left[(1+j)\theta \right]}^{2+j+\frac{r}{\alpha }}}}\right\}.\hfill \end{array}$$

The incomplete $r^{th}$ moment of X defined as $m_r\left(z\right)=\underset{0}{\overset{z}{\int }}{x}^{r}f\left(x\right)dx$ can be written;

$$m_r\left(z\right)={\displaystyle \frac{c\alpha {\theta }^2}{{\theta }^2+1}}{\int }_0^z(\theta +x^{\alpha })x^{r+\alpha -1}{e}^{-\theta x^{\alpha }}{\left\{1-\left[1+{\displaystyle \frac{\theta x^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x^{\alpha }}\right\}}^{c-1}dx.$$

The same steps used in the complete moment yield the incomplete moment, hence

$$\begin{array}{c}\hfill m_r\left(z\right)=c\sum _{j=0}^{\infty }\sum _{i=j}^{\infty }{(-1)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right){\displaystyle \frac{{\theta }^{j+2}}{{({\theta }^2+1)}^{j+1}}}\left\{{\displaystyle \frac{\theta \gamma \left[1+j+\frac{r}{\alpha },(1+j)\theta z^{\alpha }\right]}{{\left[(1+j)\theta \right]}^{1+j+\frac{r}{\alpha }}}}+{\displaystyle \frac{\gamma \left[2+j+\frac{r}{\alpha },(1+j)\theta z^{\alpha }\right]}{{\left[(1+j)\theta \right]}^{2+j+\frac{r}{\alpha }}}}\right\}.\end{array}$$

The Rény entropy of $X\sim$ EPS $(c,\alpha ,\theta )$ is

$$R_E={\displaystyle \frac{1}{1-\omega }}log\left\{\underset{0}{\overset{\infty }{\int }}{f}^{\omega }(x;c,\alpha ,\theta )dx\right\},$$

where $f(x;c,\alpha ,\theta )$ is the pdf of EPS distribution and $\omega \ne 1$. Using the general binomial expansion for the decomposition, the Rény entropy is given as

$$\begin{array}{cc}\hfill R_E& ={\displaystyle \frac{1}{1-\omega }}log\{{\left({\displaystyle \frac{c}{{\theta }^2+1}}\right)}^{\omega }\sum _{k=0}^{\infty }\sum _{i=j}^{\infty }\sum _{j=k}^{\infty }\sum _{h=0}^k{(-1)}^{i+k}\left({\displaystyle \genfrac{}{}{0pt}{}{\omega (c-1)}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{j+k-1}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\omega }{h}}\right)\hfill \\ \hfill & \times {\theta }^{2k+j+h-2\omega }{\displaystyle \frac{\mathsf{\Gamma }\left(\alpha (j-h+2)-\omega +1\right)}{{\left[(\omega +i)\theta \right]}^{\alpha (j-h+2)-\omega +1}}}\}.\hfill \end{array}$$

The moment generating function $M_X\left(t\right)$ and the characteristic function ${\Phi }_X\left(it\right)$ are respectively

$$\begin{array}{cc}\hfill M_X\left(t\right)& =\left({\displaystyle \frac{c}{{\theta }^2+1}}\right)\sum _{j=0}^{\infty }\sum _{i=j}^{\infty }\sum _{r=0}^{\infty }\sum _{k=0}^j{(-1)}^{i+k}\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{j+k-1}{k}}\right){\displaystyle \frac{t^r}{r!}}{\theta }^{2-j-2k}\hfill \\ \hfill & \times \left\{{\displaystyle \frac{\theta \mathsf{\Gamma }\left[r+\alpha (j+1)\right]}{{\left[(i+1)\theta \right]}^{r+\alpha (j+1)}}}{\displaystyle \frac{\mathsf{\Gamma }\left[r+\alpha (j+2)\right]}{{\left[(i+1)\theta \right]}^{r+\alpha (j+2)}}}\right\},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\Phi }_X\left(it\right)& =\left({\displaystyle \frac{c}{{\theta }^2+1}}\right)\sum _{j=0}^{\infty }\sum _{i=j}^{\infty }\sum _{r=0}^{\infty }\sum _{k=0}^j{(-1)}^{i+k}\left({\displaystyle \genfrac{}{}{0pt}{}{c-1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{j+k-1}{k}}\right){\displaystyle \frac{{\left(it\right)}^r}{r!}}{\theta }^{2-j-2k}\hfill \\ \hfill & \times \left\{{\displaystyle \frac{\theta \mathsf{\Gamma }\left[r+\alpha (j+1)\right]}{{\left[(i+1)\theta \right]}^{r+\alpha (j+1)}}}{\displaystyle \frac{\mathsf{\Gamma }\left[r+\alpha (j+2)\right]}{{\left[(i+1)\theta \right]}^{r+\alpha (j+2)}}}\right\}.\hfill \end{array}$$

Figure (Figure 5) evidently illustrates the influence of parameters c and $\alpha$ on the most significant statistical moments of the EPS distribution: mean, variance, skewness, and kurtosis. Mean and variance plots exhibit a typically homogeneous tendency. Both moments increase as a function of c and $\alpha$. Exactly, for a fixed $\alpha$, there is an initial steep increase with c (by some 2.0-2.5), after which the rate of increase slows down, suggesting saturation or asymptotic behavior. Similar positive correlation with increasing $\alpha$ is observed, particularly when c is small. This close similarity in behaviour between the mean and variance indicates that while the central tendency of the distribution itself varies as a result of changes in c and $\alpha$, its overall variability increases proportionally.

In contrast, plots for skewness and kurtosis portray a more intricate and non-linear association with the parameters. For skewness, although high positive values are seen when both c and $\alpha$ are low, the distribution becomes more symmetric (skewness gets closer to zero) as they increase. Interestingly, there are areas where a rise in c may reduce skewness at first, even to negative values in some instances, before a possible turnaround. The kurtosis plot also exhibits dynamic behavior. Extremely high kurtosis, indicating heavy tails and spiked peakness (leptokurtic distributions) occurs when c and $\alpha$ are low. However, when c and $\alpha$ increase in size, the kurtosis decreases, suggesting a move towards distributions with lighter tails and less spiky peaks, perhaps near to mesokurtic or even platykurtic ones. These plots illustrate that while c and $\alpha$ always influence the location and scale (mean and variance) of the EPS distribution, their influence on its shape (skewness and kurtosis) is considerably more complex and non-monotonic. The discussed sensitivity of skewness and kurtosis to the smaller values of c and $\alpha$ illustrates the way small changes in these parameters can significantly alter the asymmetry and tail behavior of the distribution. Conversely, for higher values of the parameters, the distribution will converge to a less extreme and more symmetric shape. Such intimate knowledge is central to correct interpretation of model fits as well as parameter selection in practical applications of the EPS distribution.

4. Maximum Likelihood Estimation

Let $x_1,x_2,\dots ,x_n$ be n independent and identically distributed random variables from the EPS distribution, and let $\Im =(c,\theta ,\alpha )$ be the vector of unknown parameters. Then, the log-likelihood function for ℑ is given by:

$$\begin{array}{cc}\hfill \ell (\Im )& =n\left[logc+log\alpha +2log\theta -log({\theta }^2+1)\right]-\theta \sum _{i=1}^n{x}_i^{\alpha }+\sum _{i=1}^{n}log(\theta +x_i^{\alpha })\hfill \\ \hfill & +(\alpha -1)\sum _{i=1}^{n}log\left(x_i\right)+(c-1)\sum _{i=1}^{n}log\left\{1-\left[1+{\displaystyle \frac{\theta x_i^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x_i^{\alpha }}\right\}.\hfill \end{array}$$

Taking partial derivatives of the log-likelihood with respect to the parameters yields:

$${\displaystyle \frac{\partial \ell }{\partial c}}={\displaystyle \frac{n}{c}}+\sum _{i=1}^{n}log\left\{1-\left[1+{\displaystyle \frac{\theta x_i^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x_i^{\alpha }}\right\}.$$

Hence, setting ${\displaystyle \frac{\partial \ell }{\partial c}}=0$, we obtain the MLE for c:

$$\widehat{c}=-{\displaystyle \frac{n}{\sum _{i=1}^{n}log\left\{1-\left[1+\frac{\theta x_i^{\alpha }}{{\theta }^2+1}\right]e^{-\theta x_i^{\alpha }}\right\}}}.$$

(8)

The partial derivative with respect to $\alpha$ is:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \ell }{\partial \alpha }}& ={\displaystyle \frac{n}{\alpha }}-\theta \sum _{i=1}^n{x}_i^{\alpha }log{x}_i+\sum _{i=1}^n{\displaystyle \frac{x_i^{\alpha }log{x}_i}{\theta +x_i^{\alpha }}}+\sum _{i=1}^{n}log\left(x_i\right)\hfill \\ \hfill & +(c-1)\sum _{i=1}^n{\displaystyle \frac{log{x}_i{e}^{-\theta x_i^{\alpha }}\left[\theta ({\theta }^2+1)x_i^{\alpha }-\theta x_i^{\alpha }+{\theta }^2{x}_i^{2\alpha }\right]}{({\theta }^2+1)\left\{1-\left[1+\frac{\theta x_i^{\alpha }}{({\theta }^2+1)}\right]e^{-\theta x_i^{\alpha }}\right\}}}.\hfill \end{array}$$

(9)

The partial derivative with respect to $\theta$ is:

$${\displaystyle \frac{\partial \ell }{\partial \theta }}={\displaystyle \frac{2n}{\theta }}-{\displaystyle \frac{2n\theta }{{\theta }^2+1}}-\sum _{i=1}^n{x}_i^{\alpha }+\sum _{i=1}^n{\displaystyle \frac{1}{\theta +x_i^{\alpha }}}-(c-1)\sum _{i=1}^n{\displaystyle \frac{e^{-\theta x_i^{\alpha }}\left[\frac{x_i^{\alpha }(1-{\theta }^2)}{{({\theta }^2+1)}^2}-x_i^{\alpha }\left(1+\frac{\theta x_i^{\alpha }}{{\theta }^2+1}\right)\right]}{1-\left[1+\frac{\theta x_i^{\alpha }}{({\theta }^2+1)}\right]e^{-\theta x_i^{\alpha }}}}.$$

(10)

Setting ${\displaystyle \frac{\partial \ell }{\partial \alpha }}=0$ and ${\displaystyle \frac{\partial \ell }{\partial \theta }}=0$ and solving numerically, the MLEs for $\theta$ and $\alpha$ are obtained. These can be computed in R using the optim() function. Let

$$\begin{array}{cc}\hfill A_i& =1+{\displaystyle \frac{\theta x_i^{\alpha }}{{\theta }^2+1}},B_i=e^{-\theta x_i^{\alpha }}{D}_i=1-A_i{B}_i=1-\left[1+{\displaystyle \frac{\theta x_i^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x_i^{\alpha }}\hfill \\ \hfill M_i& ={\displaystyle \frac{x_i^{\alpha }(1-{\theta }^2)}{{({\theta }^2+1)}^2}}-x_i^{\alpha }{A}_i,N_i=B_i{M}_i,E_i=\theta ({\theta }^2+1)x_i^{\alpha }-\theta x_i^{\alpha }+{\theta }^2{x}_i^{2\alpha }.\hfill \end{array}$$

The second partial derivatives are then given by:

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2\ell }{\partial c^2}}& =-{\displaystyle \frac{n}{c^2}};{\displaystyle \frac{{\partial }^2\ell }{\partial c\partial \theta }}={\displaystyle \frac{{\partial }^2\ell }{\partial \theta \partial c}}=-\sum _{i=1}^n{\displaystyle \frac{B_i{M}_i}{D_i}};{\displaystyle \frac{{\partial }^2\ell }{\partial c\partial \alpha }}={\displaystyle \frac{{\partial }^2\ell }{\partial \alpha \partial c}}=\sum _{i=1}^n{\displaystyle \frac{log{x}_i{B}_i{E}_i}{({\theta }^2+1)D_i}}\hfill \\ \hfill {\displaystyle \frac{{\partial }^2\ell }{\partial {\theta }^2}}& =-{\displaystyle \frac{2n}{{\theta }^2}}-{\displaystyle \frac{2n(1-{\theta }^2)}{{({\theta }^2+1)}^2}}-\sum _{i=1}^n{\displaystyle \frac{1}{{(\theta +x_i^{\alpha })}^2}}-(c-1)\sum _{i=1}^n{\displaystyle \frac{\frac{\partial N_i}{\partial \theta }{D}_i-N_i\frac{\partial D_i}{\partial \theta }}{D_i^2}}\hfill \\ \hfill & \mathrm{where}{\displaystyle \frac{\partial D_i}{\partial \theta }}=-B_i{x}_i^{\alpha }\left({\displaystyle \frac{1-{\theta }^2}{{({\theta }^2+1)}^2}}-A_i\right)\mathrm{and}{\displaystyle \frac{\partial N_i}{\partial \theta }}={\displaystyle \frac{\partial \left(B_i{M}_i\right)}{\partial \theta }}\hfill \\ \hfill {\displaystyle \frac{{\partial }^2\ell }{\partial {\alpha }^2}}& =-{\displaystyle \frac{n}{{\alpha }^2}}-\theta \sum _{i=1}^n{x}_i^{\alpha }{(log{x}_i)}^2+\sum _{i=1}^n{\displaystyle \frac{\theta x_i^{\alpha }{(log{x}_i)}^2}{{(\theta +x_i^{\alpha })}^2}}+(c-1)\sum _{i=1}^n{\displaystyle \frac{log{x}_i}{({\theta }^2+1)D_i^2}}\left[\left({\displaystyle \frac{\partial \left(B_i{E}_i\right)}{\partial \alpha }}\right)D_i-\left(B_i{E}_i\right)\left({\displaystyle \frac{\partial D_i}{\partial \alpha }}\right)\right]\hfill \\ \hfill & \mathrm{where}{\displaystyle \frac{\partial D_i}{\partial \alpha }}={\displaystyle \frac{\theta x_i^{\alpha }log{x}_i{B}_i({\theta }^2+\theta x_i^{\alpha })}{{\theta }^2+1}}\mathrm{and}{\displaystyle \frac{\partial \left(B_i{E}_i\right)}{\partial \alpha }}=B_i{x}_i^{\alpha }log{x}_i\left(-\theta E_i+\theta ({\theta }^2+1)-\theta +2{\theta }^2{x}_i^{\alpha }\right)\hfill \\ \hfill {\displaystyle \frac{{\partial }^2\ell }{\partial \alpha \partial \theta }}& ={\displaystyle \frac{{\partial }^2\ell }{\partial \theta \partial \alpha }}=-\sum _{i=1}^n{x}_i^{\alpha }log{x}_i-\sum _{i=1}^n{\displaystyle \frac{x_i^{\alpha }log{x}_i}{{(\theta +x_i^{\alpha })}^2}}-(c-1)\sum _{i=1}^n\left({\displaystyle \frac{\frac{\partial N_i}{\partial \alpha }{D}_i-N_i\frac{\partial D_i}{\partial \alpha }}{D_i^2}}\right)\hfill \\ \hfill & \mathrm{where}{\displaystyle \frac{\partial D_i}{\partial \alpha }}={\displaystyle \frac{\theta x_i^{\alpha }log{x}_i{B}_i({\theta }^2+\theta x_i^{\alpha })}{{\theta }^2+1}}\mathrm{and}{\displaystyle \frac{\partial N_i}{\partial \alpha }}=B_i{x}_i^{\alpha }log{x}_i\left(-\theta M_i+\left[{\displaystyle \frac{1-{\theta }^2}{{({\theta }^2+1)}^2}}-1-{\displaystyle \frac{2\theta x_i^{\alpha }}{{\theta }^2+1}}\right]\right)\hfill \end{array}$$

The Hessian matrix $H(\Im )$ for the parameter vector $\Im =(c,\theta ,\alpha )$ is a $3\times 3$ symmetric matrix:

$$H(\Im )=\left(\begin{array}{ccc}{\displaystyle \frac{{\partial }^2\ell }{\partial c^2}}& {\displaystyle \frac{{\partial }^2\ell }{\partial c\partial \theta }}& {\displaystyle \frac{{\partial }^2\ell }{\partial c\partial \alpha }}\\ {\displaystyle \frac{{\partial }^2\ell }{\partial \theta \partial c}}& {\displaystyle \frac{{\partial }^2\ell }{\partial {\theta }^2}}& {\displaystyle \frac{{\partial }^2\ell }{\partial \theta \partial \alpha }}\\ {\displaystyle \frac{{\partial }^2\ell }{\partial \alpha \partial c}}& {\displaystyle \frac{{\partial }^2\ell }{\partial \alpha \partial \theta }}& {\displaystyle \frac{{\partial }^2\ell }{\partial {\alpha }^2}}\end{array}\right).$$

The concept of MLE relies on finding the parameter values that maximize the likelihood function (or equivalently, the log-likelihood function) for a given set of observed data. While the method is widely used, it is essential to understand the conditions under which the MLE is guaranteed to exist and be unique. Without such guarantees, numerical optimization procedures like the optim() function in R might converge to a local maximum, fail to converge, or yield multiple possible solutions, leading to ambiguous or incorrect parameter estimates. For a general parametric statistical model, let $f\left(x\right|\Im )$ be the PDF of the random variable X, where ℑ is the vector of unknown parameters belonging to a parameter space $\Theta$. Given n independent and identically distributed (i.i.d.) observations $x_1,x_2,\dots ,x_n$, the log-likelihood function is given by:

$\ell (\Im )={\displaystyle \sum _{i=1}^n}logf\left(x_i\right|\Im ).$ The MLE $\widehat{\Im }$ is defined as, $\widehat{\Im }=arg\underset{\Im \in \Theta }{max}\ell (\Im ).$ The existence and uniqueness of $\widehat{\Im }$ depend on the properties of the log-likelihood function $\ell (\Im )$ and the parameter space $\Theta$.

4.1. Conditions for Existence

A maximum likelihood estimator is guaranteed to exist under the following general conditions:

• The log-likelihood function $\ell (\Im )$ must be continuous over the parameter space $\Theta$. This ensures that the function does not have "jumps" that would prevent a maximum from being attained. For the EPS distribution, the log-likelihood function $\ell (\Im )$ is a composite of continuous functions (logarithms, exponentials, sums, products), so it is generally continuous wherever it is defined. However, careful consideration must be given to the domains of logarithmic terms, e.g., $logc$, $log\alpha$, $log\theta$, $log({\theta }^2+1)$, $log{x}_i$, and the term $log\left\{1-\left[1+{\displaystyle \frac{\theta x_i^{\alpha }}{{\theta }^2+1}}\right]e^{-\theta x_i^{\alpha }}\right\}$. These terms imply that $c>0$, $\alpha >0$, $\theta >0$, $x_i>0$, and the argument of the last logarithm must be strictly positive.
• If the parameter space $\Theta$ is a compact set (i.e., closed and bounded), and the log-likelihood function is continuous on $\Theta$, then by the Extreme Value Theorem, a maximum of $\ell (\Im )$ is guaranteed to exist within $\Theta$. In practice, parameter spaces for distributions are often open (e.g., $(0,\infty )$ for scale or shape parameters), which are not compact. In such cases, alternative conditions are needed, such as:

  (a)

  $\ell (\Im )\to -\infty$ as ℑ approaches the boundary of $\Theta$ or goes to infinity: This condition ensures that the maximum does not lie on the boundary or at infinity, forcing it into the interior of the parameter space. For the EPS distribution, for example, as $c\to 0^{+}$ or $\alpha \to 0^{+}$ or $\theta \to 0^{+}$, certain terms in $\ell (\Im )$ tend to $-\infty$, which is a desirable property for ensuring an interior maximum. Similarly, if parameters become excessively large, the likelihood should tend to zero.

4.2. Conditions for Uniqueness

Even if an MLE exists, it may not be unique. Uniqueness is typically ensured by:

• If the log-likelihood function $\ell (\Im )$ is strictly concave over the parameter space $\Theta$, and if a maximum exists, then that maximum is unique. A strictly concave function has at most one global maximum.

  (a)

  Mathematically, strict concavity can be checked by examining the Hessian matrix of the log-likelihood function, $H(\Im )=\left[{\displaystyle \frac{{\partial }^2\ell }{\partial {\Im }_j\partial {\Im }_k}}\right]$. If $H(\Im )$ is negative definite for all $\Im \in \Theta$, then $\ell (\Im )$ is strictly concave. This is often the most challenging condition to verify analytically for complex likelihood functions.

  (b)

  For the EPS distribution, one would need to compute the second partial derivatives with respect to c, $\alpha$, and $\theta$, and then form the $3\times 3$ Hessian matrix. Demonstrating that this matrix is negative definite across the entire parameter space can be very difficult.

4.3. Implications for the EPS Distribution

For the EPS distribution and its MLEs for $(c,\theta ,\alpha )$:

• The analytical solution for $\widehat{c}$ (Equation 8) implicitly depends on $\theta$ and $\alpha$. This means that $\widehat{c}$ exists uniquely given specific values of $\theta$ and $\alpha$. However, the global MLE for c depends on the numerically optimized values of $\theta$ and $\alpha$.
• Since the likelihood equations for $\theta$ and $\alpha$ (Equations 9 and 10 set to zero) do not have closed-form solutions, numerical methods like optim() are necessary. The success and reliability of these numerical methods heavily depend on the underlying properties of the log-likelihood function related to existence and uniqueness:

  (a)

  Local Maxima: If the log-likelihood function for the EPS distribution is not strictly concave, optim() might converge to a local maximum instead of the global maximum. This highlights the importance of choosing good initial values for the optimization (e.g., from moment estimators or a grid search) and potentially running the optimization from multiple starting points to increase confidence in finding the global optimum.

  (b)

  Boundary Issues: The parameters $c,\theta ,\alpha$ must be positive. The numerical optimization routine must be constrained to respect these bounds (e.g., using method="L-BFGS-B" with lower arguments in R’s optim() function). If the true MLE lies on the boundary of the parameter space, standard derivative-based methods may not apply directly, and the maximum might not be at a point where the gradient is zero.

  (c)

  Identifiability: A fundamental prerequisite for uniqueness is that the model must be identifiable, meaning that different parameter vectors must correspond to different probability distributions. If the EPS distribution is not identifiable (i.e., different combinations of $c,\theta ,\alpha$ can produce the exact same distribution), then unique MLEs cannot exist. While most well-defined distributions are identifiable, it’s a theoretical point to consider.

5. Regression

Suppose one makes the transformation $Y=log\left(X\right)$ where $X\sim$ EPS $(c,\theta ,\alpha )$ with pdf in Equation (6). Set $\alpha ={\displaystyle \frac{1}{\sigma }}$ and $\theta =e^{-{\displaystyle \frac{\mu }{\sigma }}}$, the LEPS density, cdf, and survival function for $y\in \mathbb{R}$ are given by

$$\begin{array}{cc}\hfill f(y;c,\sigma ,\mu )& ={\displaystyle \frac{c{e}^{\left(\frac{y-2\mu }{\sigma }\right)}}{\sigma \left(e^{-\frac{2\mu }{\sigma }}+1\right)}}\left(e^{-{\displaystyle \frac{\mu }{\sigma }}}+e^{{\displaystyle \frac{y}{\sigma }}}\right)e^{-e^{\left({\displaystyle \frac{y-\mu }{\sigma }}\right)}}\hfill \\ \hfill & \times {\left\{1-\left[1+{\displaystyle \frac{e^{\left(\frac{y-\mu }{\sigma }\right)}}{e^{-\frac{2\mu }{\sigma }}+1}}\right]e^{-e^{\left({\displaystyle \frac{y-\mu }{\sigma }}\right)}}\right\}}^{c-1},\hfill \end{array}$$

(11)

$$F(y;c,\sigma ,\mu )={\left\{1-\left[1+{\displaystyle \frac{e^{\left(\frac{y-\mu }{\sigma }\right)}}{e^{-\frac{2\mu }{\sigma }}+1}}\right]e^{-e^{\left({\displaystyle \frac{y-\mu }{\sigma }}\right)}}\right\}}^c,$$

and

$$S(y;c,\sigma ,\mu )=1-{\left\{1-\left[1+{\displaystyle \frac{e^{\left(\frac{y-\mu }{\sigma }\right)}}{e^{-\frac{2\mu }{\sigma }}+1}}\right]e^{-e^{\left({\displaystyle \frac{y-\mu }{\sigma }}\right)}}\right\}}^c,$$

where $c,\theta ,\sigma >0$ and $\mu \in \mathbb{R}$. Hence, for $X\sim$ EPS $(c,\theta ,\alpha )$, $Y=log\left(X\right)\sim$ LEPS $(c,\sigma ,\mu )$. Define $Z={\displaystyle \frac{(Y-\mu )}{\sigma }}$, the survival and density functions are respectively

$$S(z;c,\sigma ,\mu )=1-{\left\{1-\left[1+{\displaystyle \frac{e^z}{e^{-\frac{2\mu }{\sigma }}+1}}\right]e^{-e^z}\right\}}^c,$$

and

$$f(z;c,\sigma ,\mu )={\displaystyle \frac{cw\left(z\right)e^{-e^z}}{\sigma \left(e^{-\frac{2\mu }{\sigma }}+1\right)}}{\left\{1-\left[1+{\displaystyle \frac{e^z}{e^{-\frac{2\mu }{\sigma }}+1}}\right]e^{-e^z}\right\}}^{c-1},$$

(12)

where $w\left(z\right)=e^{z-{\displaystyle \frac{2\mu }{\sigma }}}+e^{2z}$. Using equation (11), we construct a parametric regression model for the response variable $Y_i$ and a vector of explanatory variables ${\mathbf{v}}_{\mathbf{i}}^{\prime }=(v_{i1},v_{i2},\dots ,v_{ip})$ as

$$y_i={\mathbf{v}}^{\prime }\beta +\sigma z_i,i=1,2,\dots ,n,$$

(13)

where ${\mu }_i={\mathbf{v}}^{\prime }\beta$, $\beta ={({\beta }_1,\dots ,{\beta }_p)}^{\prime }$ is the vector of unknown regression coefficients and z is the random error with density in Equation (12). Define the survival and density functions of $Y_i|{\mathbf{v}}^{\prime }$ are

$$S\left(y\right|{\mathbf{v}}^{\prime })=1-{\left\{1-\left[1+{\displaystyle \frac{e^{z_i}}{e^{-\frac{2{\mu }_i}{\sigma }}+1}}\right]e^{-e^{z_i}}\right\}}^c,$$

and

$$f\left(y\right|{\mathbf{v}}^{\prime })={\displaystyle \frac{c\left(e^{z_i-\frac{2{\mu }_i}{\sigma }}+e^{2z_i}\right)e^{-e^{z_i}}}{\sigma \left(e^{-\frac{2{\mu }_i}{\sigma }}+1\right)}}{\left\{1-\left[1+{\displaystyle \frac{e^{z_i}}{e^{-\frac{2{\mu }_i}{\sigma }}+1}}\right]e^{-e^{z_i}}\right\}}^{c-1},$$

(14)

where $w\left(z_i\right)=e^{z_i-{\displaystyle \frac{2{\mu }_i}{\sigma }}}+e^{2z_i}$ and $z_i={\displaystyle \frac{y_i-{\mu }_i}{\sigma }}.$

5.1. MLE of For Right-Censored Sample

Suppose the lifetime X of n individuals diagnosed with COVID-19 virus is EPS $(\wp )$ distributed. Let $y|{\mathbf{v}}^{\prime }$ be the response variable of a parametric regression model from the EPS $(\wp )$ distribution with pdf in Equation (14). To estimate the parameters $\beta$ of the transformed model, the n individuals are quarantined and subjected to routine treatment at the same time. After time $\left(t\right)$, $(n-m)$ individuals recovers. If the lifetimes of the other $m(>0)$ individuals are denoted by $y_1,y_2,...,y_m$. Then, the likelihood of $\beta$ can be expressed as

$$\ell \left(y\right|{\mathbf{v}}^{\prime },\beta )=\left(\prod _{i=0}^{m}f\left(y\right|{\mathbf{v}}^{\prime })\right)\prod _{i=m}^{n}P(Y>t),$$

where $f\left(y\right|{\mathbf{v}}^{\prime })$ is the pdf in equation (14). It is witty to write

$$P(Y>t)=P(Z>t)=\underset{t}{\overset{\infty }{\int }}{\displaystyle \frac{cw\left(z_i\right)e^{-e^{z_i}}}{\sigma \left(e^{-\frac{2{\mu }_i}{\sigma }}+1\right)}}{\left\{1-\left[1+{\displaystyle \frac{e^{z_i}}{e^{-\frac{2{\mu }_i}{\sigma }}+1}}\right]e^{-e^{z_i}}\right\}}^{c-1}dz.$$

The estimates of the unknown parametric regression coefficients $\beta$ are obtained using numerical iteration implemented in R.

6. Simulation Study

The EPS distribution is simulated over three scenarios to examine the validity of the MLEs. The acceptance-rejection technique is employed to simulate random samples of size $n=50,100,200,\mathrm{and}400$ from the distribution. The aforementioned procedure is replicated 1000 times, and the average estimates (AEs), the biases, and the mean squared errors (MSEs) are calculated.

The algorithm for simulating random samples is based on the acceptance-rejection method:

• Generate t from the density

  $$h\left(t\right)=\alpha \theta t^{\alpha -1}{e}^{-\theta t^{\alpha }}.$$
• Generate $u\sim \mathrm{Uniform}(0,1)$.
• If

  $$u\le {\displaystyle \frac{f\left(t\right)}{Nh\left(t\right)}},$$

  set $x=t$, where $f(·)$ is the PDF given in equation (6) and

  $$N=max\left[{\displaystyle \frac{f\left(t\right)}{h\left(t\right)}}\right].$$
  Otherwise, return to step 1.

The function $h\left(t\right)$ serves as the proposal density within the acceptance-rejection method. This particular form, $h\left(t\right)=\alpha \theta t^{\alpha -1}{e}^{-\theta t^{\alpha }}$, corresponds to the probability density function of a two-parameter Weibull distribution. Its selection is motivated by two key factors. Firstly, random samples can be efficiently generated from a Weibull distribution using standard techniques (e.g., inverse transform sampling), making it a practical choice for the proposal step. Secondly, for the acceptance-rejection method to be efficient, the proposal density $h\left(t\right)$ should ideally envelope the target distribution $f\left(t\right)$ (the PDF of the EPS distribution) when scaled by a constant N. By choosing a proposal density that shares a similar functional form or general shape with the target distribution, the acceptance rate of the algorithm is maximized, thereby reducing the computational time required to generate the desired number of samples from the EPS distribution. This efficiency is crucial for conducting a large number of replications in the simulation study.

Maximum likelihood estimates of the parameters are computed for every replication. Performance of the estimators is assessed by computing the following measures over the 1000 replications:

$$\begin{array}{cc}\hfill \mathrm{Average}\mathrm{Estimate}\left(\mathrm{AE}\right)& ={\displaystyle \frac{1}{1000}}\sum _{i=1}^{1000}{\widehat{\theta }}_i,\hfill \\ \hfill \mathrm{Bias}& =\mathrm{AE}-\theta ,\hfill \\ \hfill \mathrm{Mean}\mathrm{Squared}\mathrm{Error}\left(\mathrm{MSE}\right)& ={\displaystyle \frac{1}{1000}}\sum _{i=1}^{1000}{({\widehat{\theta }}_i-\theta )}^2,\hfill \end{array}$$

(15)

where ${\widehat{\theta }}_i$ is the estimate of the parameter $\theta$ in the $i^{\mathrm{th}}$ replication.

These estimates indicate the bias and accuracy of the MLEs under different sample sizes and parameter values. Low MSE and bias are indicators of high estimator quality. The results in Tables (Table 2) and (Table 3) indicate that with increased sample size, estimators gain precision and stability, supporting the asymptotic behavior of the MLE.

Table 2 and Table 3 present simulation results for EPS model for two population (true) parameter values: $(c=1.25,\theta =3,\alpha =2.25),(c=1.20,\theta =2.5,\alpha =1.75)$ and $(c=2.0,\theta =1.5,\alpha =0.5),(c=1.75,\theta =2.5,\alpha =0.75)$. Data was generated for sample sizes $n=50,100,300,$ and $600,$ and MLEs were calculated for each case with 1000 Monte Carlo replications. The results indicate that, for both sets of true parameters, the estimators of c, $\theta$, and $\alpha$ have positive biases that substantially decrease with increasing sample size (n). For larger sample sizes ($n=300$ and $n=600$), the biases are very small, and it can be observed that the estimators are increasingly getting more accurate (unbiased) as n increases. This is how the asymptotic properties of MLEs must be. In addition, the MSEs of all parameters have a sharp and significant decrease with a increase in sample size, showing a striking increase in the accuracy of estimates for larger samples. The sudden decline in failure rates. For $n\ge 100$, the maximum likelihood method converges perfectly for almost all replications (no failures), indicating considerably improved numerical stability of the estimation algorithm. While a few failures in the smallest sample size ($n=50$) were found for the initial set of parameters for Table 2, the values are very small compared to normal convergence failure rates.

In summary, the findings suggest that the MLEs of the EPS model behave well with low bias and MSE at large sample sizes and exhibit good convergence properties for varying sample sizes and true parameter values. These are indications that the log-likelihood function is better or the optimization process is more stable, which leads to reliable estimates of the EPS distribution.

Figure 6. (a) Bias (b) MSE plots for the LEPS Regression Simulation

Figure 6. (a) Bias (b) MSE plots for the LEPS Regression Simulation

From Table 4, the simulation results for the MLEs of the LEPS regression model, with true parameters set as $c=0.7$, $\sigma =1.0$, ${\beta }_0=0.5$, ${\beta }_1=-0.4$, ${\beta }_2=0.3$, ${\beta }_3=-0.2$, ${\beta }_4=0.1$, ${\beta }_5=-0.1$, and ${\beta }_6=0.05$, demonstrate distinct behaviors across varying sample sizes. For parameter c, its bias exhibits some fluctuation and does not consistently diminish towards zero, suggesting that larger sample sizes might be necessary for its full asymptotic properties to manifest, despite its MSE reducing significantly with increased sample size. By contrast, while the MSE for parameter $\sigma$ consistently decreases, a persistent positive bias indicates the presence of a finite-sample bias even as precision improves. The intercept ${\beta }_0$ consistently shows a negative bias; although its MSE decreases, the magnitude of this bias does not decrease monotonically, possibly reflecting complex finite-sample dynamics. In contrast, the coefficients for the covariates (${\beta }_1$ to ${\beta }_6$) exhibit strong convergence of their mean estimates to the true values. Crucially, both the magnitude of their bias and their Mean Squared Error (MSE) consistently decrease with larger sample sizes, thereby affirming the expected asymptotic properties of MLEs for these parameters and confirming enhanced accuracy and precision. Overall, the findings align with the desirable asymptotic properties of MLEs for most parameters, indicating improved accuracy and precision as sample size increases. The very low incidence of simulation failures further attests to the excellent computational stability of the MLE procedure, an observation visually supported by Figure 7.

Figure (Figure 6) presents the simulation results for the MLEs of the LEPS regression model, illustrating the bias and Mean Squared Error (MSE) of the parameter estimates across varying sample sizes. Panel (a) displays the bias of the LEPS MLEs. It is evident that as the sample size increases from 50 to 1000, the bias for most parameters, including $\sigma$ and the regression coefficients ${\beta }_0$ through ${\beta }_6$, generally converges towards zero. This diminishing bias with increasing sample size suggests that the MLEs are asymptotically unbiased, a desirable property for these estimators. However, parameter c shows a more fluctuating bias that remains somewhat distant from zero, indicating that larger sample sizes than those explored might be required for its bias to fully dissipate, or that it exhibits a higher finite-sample bias. Panel (b) illustrates the MSE of the LEPS MLEs. Consistent with theoretical expectations for well-behaved estimators, the MSE for all parameters systematically decreases as the sample size increases. This trend confirms the asymptotic efficiency of the MLEs, indicating that the precision of the parameter estimates improves with larger datasets. The rate of decrease is notably steep for smaller sample sizes (e.g., from 50 to 300), after which the decline becomes more gradual, suggesting that while improvements in precision continue, the marginal gains diminish beyond a certain sample size within the simulated range. The overall pattern across both plots supports the statistical consistency of the LEPS MLEs, as both bias and variance (implicitly captured by MSE) generally decrease with increasing sample size.

7. Applications

In this section, we demonstrate the utility of the proposed LEPS regression model in explaining the impact of clinical factors on the lifetimes of COVID-19 patients. We also show its effectiveness in fitting traditional datasets, specifically times between failures for repairable items and the death rate of HIV/AIDS patients in Germany.

7.1. Fitting LEPS Regression Model to Lifetime of COVID-19 Patients

The dataset comprises the lifetime (in days) of 322 individuals diagnosed with COVID-19 through RT-PCR screening in Campinas, Brazil. These data were previously studied by [35]. A notable characteristic of this sample, is that approximately $66.45\%$ of the individuals have censored observations. This implies that for these participants, the event of death due to COVID-19 was not observed within the study’s follow-up period. For censored cases, we know only that their survival time extends beyond the recorded observation duration, meaning they were alive at the last point of contact or at the study’s conclusion. The variable `ind` serves as a unique, non-informative identifier for each participant in the dataset. Its role is solely for record management and to differentiate individual entries, carrying no inherent statistical value or order. Based on Equation (13), the variable ’censor’ is ${\delta }_i$ which is the crucial censoring indicator. It is a binary variable that takes a value of 0 if the observation for individual i is censored, indicating that the individual was alive at the end of the follow-up period or was lost to follow-up, without experiencing the event. Conversely, ${\delta }_i$ takes a value of 1 if the observed lifetime event, specifically death due to COVID-19, occurred for individual i. Age is $v_{i1}$ representing the age of individual i, measured in full years. Age is a continuous covariate widely recognized as a significant demographic factor influencing the severity and outcomes of COVID-19 infection. The variable ’heart’ is $v_{i2}$. This binary variable would indicate the presence or absence of heart disease (e.g., cardiovascular disease, ischemic heart disease, heart failure) for individual i. A value of 1 would denote the presence of such a condition, and 0 its absence. Cardiac comorbidities are well-established risk factors for severe COVID-19. The variable ’asthma’ is $v_{i3}$, a binary variable which captures the presence or absence of asthma for individual i. A value of 1 would indicate a diagnosis of asthma, and 0 its absence. Respiratory conditions, including asthma, are frequently investigated as potential modifiers of COVID-19 severity. The variable ’diab’ is $v_{i4}$, a binary variable indicating the presence or absence of diabetes mellitus for individual i. A value of 1 signifies that the individual has a diagnosis of diabetes mellitus (or that this condition was reported), while 0 indicates the absence of diabetes or that this information was not available or not applicable. Diabetes is a common comorbidity consistently linked with increased risk of severe COVID-19 and adverse outcomes. The variable ’neuro’ $v_{i5}$, a binary variable which denote the presence or absence of a neurological condition for individual i. A value of 1 would signify the existence of such a condition, and 0 its absence. Neurological involvement and pre-existing neurological disorders have been identified in relation to COVID-19 outcomes. The variable ’obesity’ $v_{i6}$, a binary variable which indicates whether individual i is classified as obese. A value of 1 would signify obesity, and 0 its absence. Obesity is widely recognized as a significant risk factor for severe COVID-19 and increased mortality.

Figure (Figure 7a) represents the normal probability graph that shows that the lifetimes of the COVID-19 patients follow the normal distribution. Figure (Figure 7b) is the quantile residual plot against their index (the order of observations) in the dataset. From the plot, there are potential outliers as points lie below and above the lower and upper implicit boundaries (-3,3). Figure (Figure 7c) is the histogram of the age of the patients. It is evident from the histogram that the data is right-skewed.

For the fitting of the LEPS regression model to the lifetimes (in days) of COVID-19 patients, several established parametric regression models were selected for comparative analysis. These include the log-exponentiated power Chris-Jerry (LEPCJ) by [36], log-exponentiated power Akash (LEPA) by [37], log-power Prakaamy (LPP), log-exponentiated power Lindley (LEPL), log-power Rama (LPR), log-power Lomax (LPLo), log-power Zeghdoudi (LPZ), log-power Suja (LPSU), log-exponentiated power Ishita (LEPI), log-power Ishita (LPI), and the log-exponentiated Weibull (LEW) model.

Table (Table 5) shows comparisons of the model selection metrics of various distributions to the COVID-19 data. LEPS achieves the lowest Akaike Information Criterion (AIC = 426.4954), Consistent AIC (CAIC = 427.3471), and Hannan-Quinn Information Criterion (HQIC = 440.0577), and thus ranks first by consequence. LEPI and LEW models follow with slightly higher values in these metrics. In contrast, models such as LPP and LEPR have very much larger information criteria values, which refer to poorer fit. Bayesian Information Criterion (BIC) ranks LPZ last among all the best-fitting models, even though its AIC and CAIC values are higher than those of LEPS. Overall, the results reveal that the LEPS distribution presents the most consistent and best-fitting model for the data among the competitors under consideration.

For each distributional model in Table (Table 6), the regression model is

$$y_i={\beta }_0+{\beta }_1·\mathrm{Age}+{\beta }_2·\mathrm{Heart}+{\beta }_3·\mathrm{Asthma}+{\beta }_4·\mathrm{Diabetes}+{\beta }_5·\mathrm{Neurological}\mathrm{Condition}+{\beta }_6·\mathrm{Obesity}.$$

The statistically significant predictors (meaning $p<0.05$) within models tend to consist of age, heart disease, and neurological condition, and sometimes diabetic condition, which are constant negative or positive effects depending on the model. The signs of coefficients define the effect direction: age (usually negative and significant), heart disease (negative and significant), neurological condition (negative and significant in most models), Others (asthma, diabetic condition and obesity) vary less in significance. Parameter c distinguishes between models and is linked with the shape or flexibility of fit of the distribution, and $\sigma$ is representative of scale with varying precision.

For the LEPS model, the regression model for the lifetimes (in days) of COVID-19 patients ($y_i$) is given by:

$$\begin{array}{cc}\hfill y_i& =4.52354-0.01763\mathrm{Age}-0.30804\mathrm{Heart}+0.08840\mathrm{Asthma}-0.53958\mathrm{Diabetes}\hfill \\ \hfill & -0.43946\mathrm{Neurological}\mathrm{Condition}-0.05539\mathrm{Obesity}.\hfill \end{array}$$

Standard errors are given in parentheses, and p-values are shown in brackets below:

$$\begin{array}{cc}\hfill {\beta }_0& =4.52354\left(0.35328\right)\hfill \\ \hfill {\beta }_1& =-0.01763\left(0.00445\right),p=9.20958\times {10}^{-5}\hfill \\ \hfill {\beta }_2& =-0.30804\left(0.12678\right),p=1.56573\times {10}^{-2}\hfill \\ \hfill {\beta }_3& =0.08840\left(0.11187\right),p=0.43002\hfill \\ \hfill {\beta }_4& =-0.53958\left(0.31905\right),p=0.09176\hfill \\ \hfill {\beta }_5& =-0.43946\left(0.15208\right),p=0.00412\hfill \\ \hfill {\beta }_6& =-0.05539\left(0.18159\right),p=0.76054.\hfill \end{array}$$

The scale and shape parameters of the LEPS distribution are $c=0.32055,\sigma =0.39854(\mathrm{with}\mathrm{standard}\mathrm{error}0.17345).$ These are specific to the LEPS distribution; where

(i)

Intercept (${\beta }_0$): Estimated lifetime for the baseline patient with all predictors at zero (i.e. youngest age, no heart disease, no asthma, no diabetes, no neurological condition, no obesity) is about 4 years and 6 months.

(ii)

Age (${\beta }_1=-0.01763$): For every additional year of age, the lifetime is expected to decrease by roughly 0.018 years, holding all else equal. This suggests a negative correlation between survival time and age.

(iii)

Heart disease (${\beta }_2=-0.30804$): Patients with heart disease suffers a lifetime reduction by about 0.31 years.

(iv)

Asthma (${\beta }_3=0.08840$): Asthma appears to be positively correlated with lifetime, increasing it by about 0.088 years but it is not significant statistically since (p-value $>0.05$).

(v)

Diabetes (${\beta }_4=-0.53958$): lifetime for diabetic patients is significantly low at about 0.54 years, suggesting diabetes seriously reduces the life expectancy of a COVID-19 patient.

(vi)

Neurological condition (${\beta }_5=-0.43946$): Neurological conditions also decrease lifetime by about 0.44 years.

(vii)

Obesity (${\beta }_6=-0.05539$): Obesity has a very weak negative effect on lifetime, but the effect is small and not statistically significant.

This model quantifies to what degree clinical parameters influence COVID-19 patients’ lifetimes through differences in survivals, which identifies high-risk groups. Across all the models considered, a consistent pattern emerges regarding the significance of covariates. Specifically, the parameters ${\beta }_1$, ${\beta }_2$, and ${\beta }_5$ were found to be statistically significant in every model. These parameters correspond to age, heart disease, and neurological condition, respectively, suggesting their strong association with the response variable under study. For models such as LEPS, LEPCJ, LEPA, LEPL, LEPR, LPL, LEPI, and LEW, the parameters c, $\sigma$, ${\beta }_0$, ${\beta }_3$, ${\beta }_4$, and ${\beta }_6$ were consistently not significant. This implies that these covariates do not have a meaningful contribution to the outcome in these models. In contrast, for the LPP, LPZ, LEPSU, and LPI models, the parameter c is not applicable, yet the insignificance of $\sigma$, ${\beta }_0$, ${\beta }_3$, ${\beta }_4$, and ${\beta }_6$ remains consistent. It is noteworthy that despite the structural differences in the models (such as the presence or absence of parameter c), the core set of significant parameters (${\beta }_1$, ${\beta }_2$, and ${\beta }_5$) remains unchanged. This reinforces the robustness of their influence and suggests that interventions or predictions based on these variables are likely to be reliable across various modeling frameworks. The repeated insignificance of parameters like ${\beta }_0$, ${\beta }_3$, ${\beta }_4$, and ${\beta }_6$ across models indicates that these variables may be redundant or have weak explanatory power in the context of this analysis. It may be beneficial in future studies to consider removing or replacing these covariates to streamline model complexity without sacrificing predictive accuracy.

The 95% confidence intervals for the parameters are computed with the following formula, $\mathrm{CI}={\widehat{\beta }}_j\pm 1.96\times \mathrm{SE}\left({\widehat{\beta }}_j\right).$ For instance, in the LEPS model, the coefficient for age (${\beta }_1$) had a $[95\%\mathrm{CI}=-0.01763\pm (1.96\times 0.00445)=[-0.02635,-0.00891]$, indicating a statistically significant and negative association with lifetime. Confidence intervals for other parameters provide similar insights and can be found in Tables (Table 7) and (Table 8). Above all, this consistency in parameter significance across different link functions and error structures provides strong evidence for the reliability of key predictors and offers direction for model refinement. Confidence intervals were computed for each parameter estimate to assess the precision of the estimates.

7.2. Fitting EPS to Reliability Data

The first data set represent the times between failures for repairable items initially studied by [38] and reported in Table (Table 9). The second data set is Death rate due to HIV/AIDS in Germany from year 2000 to 2020 reported in https://platform.who.int/mortality/themes/theme-details/topics/indicator-groups/indicator-group-details/MDB/hiv-aids and presented in Table (Table 10). In Table (Table 11), we present a summary of basic statistics to gain insight from the two data sets. From Table (Table 11), the repairs times dataset $(n=30)$ is smaller on its mean (0.0154) and median (0.0124), and has a steep right skew (skewness = 1.2955) that indicates a longer right tail. It is also smaller with an IQR of 0.0123 and a range of 0.0462. Two outliers $(0.0436,0.0473)$ indicate extreme values. A high kurtosis of 4.3192 indicates a more pointed distribution with heavier tails. In contrast, the HIV/AIDS mortality rate data set $(n=21)$ has a mean of 0.5203, median of 0.5395, and is very slightly left-skewed with skewness -0.3335, showing slight asymmetry toward smaller values. The inter-quartile range (IQR) is 0.1590, with a total range of 0.3786, and a standard deviation of 0.1090, showing moderate variability. Its kurtosis of 2.1101 shows a very slightly flatter distribution than normal.

The proposed EPS distribution is compared with new generalized logistic-x transformed exponential (NGLXTE) distribution introduced by [39], Logistic-Weibull distribution by [40], power Zeghdoudi distribution by [41], power Ishita distribution by [42], power Prakaamy distribution by [43], power Rama distribution by [44] and Power Lomax distribution introduced by [45]. The model selection criteria are the Akaike information criterion (AIC), consistent Akaike information criterion (CAIC), Bayesian information criterion (BIC) and Hannan-Quinn information criterion (HQIC) which are computed using;

$$\begin{array}{cc}\hfill \mathrm{AIC}& =-2\ell \left(\widehat{\Im }\right)+2k\hfill \\ \hfill \mathrm{CAIC}& =-2\ell \left(\widehat{\Im }\right)+{\displaystyle \frac{2kn}{(n-k-1)}}\hfill \\ \hfill \mathrm{BIC}& =-2\ell \left(\widehat{\Im }\right)+klog\left(n\right)\hfill \\ \hfill \mathrm{HQIC}& =-2\ell \left(\widehat{\Im }\right)+2klog\left(n\right),\hfill \end{array}$$

(16)

where $\ell \left(\widehat{\Im }\right)$ is the log-likelihood function in Equation (), k is the number if parameters, n represents the sample size and $\widehat{\Im }$ is the vector of parameters; for the EPS model, $\widehat{\Im }=(\widehat{c},\widehat{\theta },\widehat{\alpha }).$ For goodness of fit test, let an $X_{\left(1\right)},X_{\left(2\right)},\dots ,X_{\left(n\right)}$ be an ordered random sample from the EPS $(c,\theta ,\alpha )$, where c, $\theta$ and $\alpha$ are unknown; the Cramér-von Mises $\left(W_n^2\right)$, Anderson-Darling $A_n^2$ and Kolmogorov-Smirnov $\left(D_n\right)$ statistics are computed using the following;

$$\begin{array}{cc}\hfill W_n^2& ={\displaystyle \frac{1}{12n}}+\sum _{j=1}^n{\left[F\left(x_j,\widehat{c},\widehat{\theta },\widehat{\alpha }\right)-{\displaystyle \frac{2j-1}{n}}\right]}^2\hfill \\ \hfill A_n^2& =-n-\sum _{j=1}^n{\displaystyle \frac{2j-1}{n}}\left[ln\left\{F\left(x_j,\widehat{c},\widehat{\theta },\widehat{\alpha }\right)\right\}+ln\left\{1-F\left(x_j,\widehat{c},\widehat{\theta },\widehat{\alpha }\right)\right\}\right]\hfill \\ \hfill D_n& =\underset{j}{max}\left\{{\displaystyle \frac{j}{n}}-F\left(x_j,\widehat{c},\widehat{\theta },\widehat{\alpha }\right),F\left(x_j,\widehat{c},\widehat{\theta },\widehat{\alpha }\right)-{\displaystyle \frac{j-1}{n}}\right\}\hfill \end{array}$$

(17)

Table 12. Model Assessment and Goodness of Fit Criteria for the Data sets

Data	Distribution	LL	AIC	CAIC	BIC	HQIC	${\mathit{W}}_{\mathit{n}}^2$	${\mathit{A}}_{\mathit{n}}^2$	${\mathit{D}}_{\mathit{n}}$	P-value
Data-I	EPS	98.54	-191.0839	-190.1608	-186.8803	-189.7391	0.0178	0.1346	0.0660	0.9995
	NGLXTE	96.82	-189.6349	-189.1904	-186.8325	-188.7383	0.0712	0.5074	0.1219	0.7643
	Logistic-Weibull	97.74	-189.4824	-188.5593	-185.2788	-188.1376	0.0303	0.1958	0.0738	0.9968
	Power Zeghdoudi	231.89	-459.7825	-459.3380	-456.9801	-458.8860	0.0189	0.1435	0.0661	0.9994
	Power Ishita	98.24	-192.4895	-192.0450	-189.6871	-191.5929	0.0280	0.2115	0.0749	0.9960
	Power Prakaamy	98.24	-192.4895	-192.045	-189.6871	-191.5929	0.0280	0.2115	0.0749	0.9960
	Power Rama	-196.4895	-192.4895	-192.045	-189.6871	-191.5929	0.0220	0.2116	0.0748	0.9960
	Power Lomax	98.46	-190.9299	-190.0068	-186.7263	-189.5851	-	-	367.17	$2.2\times {10}^{-16}$
Data-II	EPS	17.95	-29.8916	-28.4798	-26.7580	-29.2115	0.0225	0.1958	0.0868	0.9932
	NGLXTE	17.84	-31.6816	-31.0149	-29.5925	-31.2282	0.0275	0.2202	0.0945	0.9827
	Logistic-Weibull	15.77	-25.5328	-24.1211	-22.3993	-24.8528	0.0923	0.6032	0.1170	0.9044
	Power Zeghdoudi	31.70	-59.4087	-58.7420	-57.3196	-58.9553	0.0466	0.3275	0.1054	0.9546
	Power Ishita	17.79	-31.5737	-30.9070	-29.4846	-31.1203	0.0326	0.2473	0.0999	0.9710
	Power Prakaamy	17.79	-31.5738	-30.9072	-29.4848	-31.1205	0.0326	0.2473	0.0999	0.9710
	Power Rama	-35.57	-31.5714	-30.9047	-29.4823	-31.1180	0.0326	0.2474	0.1000	0.9707
	Power Lomax	17.78	-29.5634	-28.1516	-26.4298	-28.8833	0.0756	0.5523	0.9074	2.2$\times {10}^{-16}$

Table 12. Model Assessment and Goodness of Fit Criteria for the Data sets

Data	Distribution	LL	AIC	CAIC	BIC	HQIC	${\mathit{W}}_{\mathit{n}}^2$	${\mathit{A}}_{\mathit{n}}^2$	${\mathit{D}}_{\mathit{n}}$	P-value
Data-I	EPS	98.54	-191.0839	-190.1608	-186.8803	-189.7391	0.0178	0.1346	0.0660	0.9995
	NGLXTE	96.82	-189.6349	-189.1904	-186.8325	-188.7383	0.0712	0.5074	0.1219	0.7643
	Logistic-Weibull	97.74	-189.4824	-188.5593	-185.2788	-188.1376	0.0303	0.1958	0.0738	0.9968
	Power Zeghdoudi	231.89	-459.7825	-459.3380	-456.9801	-458.8860	0.0189	0.1435	0.0661	0.9994
	Power Ishita	98.24	-192.4895	-192.0450	-189.6871	-191.5929	0.0280	0.2115	0.0749	0.9960
	Power Prakaamy	98.24	-192.4895	-192.045	-189.6871	-191.5929	0.0280	0.2115	0.0749	0.9960
	Power Rama	-196.4895	-192.4895	-192.045	-189.6871	-191.5929	0.0220	0.2116	0.0748	0.9960
	Power Lomax	98.46	-190.9299	-190.0068	-186.7263	-189.5851	-	-	367.17	$2.2\times {10}^{-16}$
Data-II	EPS	17.95	-29.8916	-28.4798	-26.7580	-29.2115	0.0225	0.1958	0.0868	0.9932
	NGLXTE	17.84	-31.6816	-31.0149	-29.5925	-31.2282	0.0275	0.2202	0.0945	0.9827
	Logistic-Weibull	15.77	-25.5328	-24.1211	-22.3993	-24.8528	0.0923	0.6032	0.1170	0.9044
	Power Zeghdoudi	31.70	-59.4087	-58.7420	-57.3196	-58.9553	0.0466	0.3275	0.1054	0.9546
	Power Ishita	17.79	-31.5737	-30.9070	-29.4846	-31.1203	0.0326	0.2473	0.0999	0.9710
	Power Prakaamy	17.79	-31.5738	-30.9072	-29.4848	-31.1205	0.0326	0.2473	0.0999	0.9710
	Power Rama	-35.57	-31.5714	-30.9047	-29.4823	-31.1180	0.0326	0.2474	0.1000	0.9707
	Power Lomax	17.78	-29.5634	-28.1516	-26.4298	-28.8833	0.0756	0.5523	0.9074	2.2$\times {10}^{-16}$

Model fit was checked using the log-likelihood (LL), Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Corrected AIC (CAIC), and Hannan-Quinn Information Criterion (HQIC). Models with higher log-likelihood and lower information criteria are considered better for parameter estimation. In Data I, Power Zeghdoudi distribution possessed the highest log-likelihood (LL = 231.89) and lowest AIC, BIC, CAIC, and HQIC, thereby indicating the best parameter estimation fit. However, the log-likelihood value is significantly higher than the rest, and the disparity may reflect an inconsistency in scaling or computation error that needs to be verified. Omitting Power Zeghdoudi, the Power Ishita and Power Prakaamy distributions were found to perform competitively with similar low information criteria. The EPS and Logistic-Weibull models also showed good performance with good values of log-likelihood and acceptable values of information criteria. The Power Lomax model was far worse based on very poor goodness-of-fit statistics and low log-likelihood. In Data II, the Power Zeghdoudi model again produced the best estimation performance, with maximum LL (31.70) and minimum information criteria. The next best were NGLXTE, Power Ishita, and Power Prakaamy, all of whom had closely competing values. The EPS and Logistic-Weibull models ranked slightly below in parameter estimation accuracy. Goodness-of-fit was assessed according to the Cramér–von Mises (W), Anderson–Darling (A), and Kolmogorov–Smirnov (KS) statistics and their corresponding p-values. Small W, A, KS statistics and large p-values for models indicate closer fit to empirical data. For Data I, all of EPS, Power Zeghdoudi, Power Ishita, Power Prakaamy, and Logistic-Weibull distributions provided excellent fit with very small KS values (0.066 to 0.075) and very large p-values (larger than 0.996). Among them, the EPS model provided the largest p-value of 0.9995, indicating the best fit to empirical distribution. On the contrary, the Power Lomax distribution provided very poor fit with very large KS statistics and a p-value near zero. Even in Data II, the same performance pattern was observed. The distributions of EPS, Power Zeghdoudi, Power Ishita, Power Prakaamy, and NGLXTE all had low KS values (0.0868 to 0.1000) and high p-values (0.9707 to 0.9932), which suggest very good reliability. The Power Lomax distribution again had very high KS statistics and essentially zero p-values, which suggest poor fit.

Table 13 presents the maximum likelihood estimates (MLEs) and standard errors of the parameters of various statistical distributions fitted to Data-I and Data-II. For the EPS distribution, the parameter estimates vary significantly across datasets, with $\widehat{c}$, $\widehat{\theta }$, and $\widehat{\alpha }$ taking higher values and standard errors in Data-I, reflecting higher variability. The NGLXTE model yields more stable estimates with relatively smaller standard errors, particularly for Data-II, reflecting higher model accuracy. The Logistic-Weibull distribution contains large standard errors in both data, particularly in $\widehat{\lambda }$ and $\widehat{\beta }$, which could indicate issues with parameter identifiability or fit to the data. The Power Zeghdoudi, Power Ishita, Power Prakaamy, and Power Rama distributions yield stable estimates between Data-I and Data-II, particularly in $\widehat{\theta }$ and $\widehat{\alpha }$, with Data-II providing lower standard errors, indicating more precise inference. Lastly, the Power Lomax distribution displays extreme differences between datasets in parameter size and standard error, most notably in $\widehat{\alpha }$ and $\widehat{\lambda }$, with Data-II estimates being considerably larger. Overall, the results indicate that parameter estimates and their errors vary quite noticeably between models and datasets, with certain distributions exhibiting more stability and precision for specific data.

In Figure (Figure 8), the data plot for Data-I overlays a violin plot with an embedded boxplot, and a histogram with a superimposed density curve. The histogram indicates right-skewness with dense low values and an elongated tail. The red density curve verifies the skewness and identifies the peak at a low value. The violin plot also visually supports the shape of the distribution. The box plot shows a small median and narrowed interquartile range (IQR), which is characteristic of right-skewness. For Data-II, the histogram shows right skewness, though perhaps less pronounced than for Data-I. The density curve supports this skewness. The violin and boxplot show slightly higher median and IQR than for Data-I. In Figure (Figure 9), plot for Data-I; Empirical CDF (red line) closely matches the EPS CDF fitted (black line / dots), indicating a good fit. The same assertion for Data-II on the basis of a similar trend as Data-I since empirical and EPS CDFs closely match. bThe empirical survival function (red line) closely matches the theoretical survival function of the EPS. The TTT plot of Data-I in Figure (Figure 10) is concave up, reflecting increasing hazard rate. Similar to Data-I, the TTT plot of Data-II reflects increasing hazard rate behavior. Data-I PP Plot’s red dots are close to the diagonal, which provides a good fit of the EPS distribution. For Data-II PP Plot, dots are also close to the diagonal and provide a strong fit for Data-II. In Figure (Figure 11), the quantiles of Data-I coincide closely with quantiles of the EPS distribution with slight discrepancies. Similarly, Data-II quantiles follow theoretical EPS quantiles closely, confirming a close fit. Both data sets are right-skewed and appear overall more concentrated towards small values for Data-I. Both exhibit increasing hazard rates as confirmed by TTT plots. The EPS distribution gives a good and consistent fit to both data sets with respect to: Empirical CDFs and survival functions, PP plots (probability correspondence), and QQ plots (quantile correspondence). Robustness of the overall fit with all the diagnostics provides the evidence of the appropriateness of the EPS distribution as the model for both the data sets.

8. Conclusion

The EPS distribution, proposed here, has strong modeling flexibility and empirical relevance for lifetime and reliability data. Through thorough structural analysis and simulation, the distribution has been proven to model a broad spectrum of density and hazard function shapes—varying from unimodal and skewed to bathtub and monotonic behaviors. Simulation results show that while the MLEs can become biased in small samples, especially for specific parameter settings, the estimators become more stable and accurate as the sample size increases. Application of the EPS model on two real data sets—repair time and HIV/AIDS mortality rates—documents its pragmatic solidness. In both cases, the EPS distribution shows an excellent fit, performing competitively or better than other proven distributions such as Power Ishita, Logistic-Weibull, and Power Zeghdoudi according to likelihood-based criteria and goodness-of-fit tests. Diagnostic plots (QQ, PP, CDF, and TTT plots) graphically authenticate the better or competitive fit of the model. Hence, the EPS distribution is a powerful member of the family of flexible life-time models, highly suitable in the broad range of applied contexts in survival analysis, reliability analysis, and public health research. Its extension to regression further broadens its suitability for covariate-dependent response modeling.