A weighted hybrid discrete probability model: Mathematical framework, statistical analysis, estimation techniques, simulation-based ranking, and goodness-of-fit evaluation for over-dispersed data

Abbreviations: HDWC: Hybrid discrete weighted combining; PMF: Probability mass function; CDF: Cumulative distribution function; HRF: Hazard rate function; RHRF: Reversed hazard rate function; IFR: Increasing failure rate; MGF: Moment generating function; PGF: Probability generating function; COV: Covariance; IOD: Index of dispersion; Sk and Ku: Skewness and kurtosis; MRL: Mean residual life; SF: Survival function; RCOV: Residual coefficient of variation; SSM: Stress-strength model; EVP: Expected value principle; EPP: Exponential premium principle; MLE: Maximum likelihood estimation; ME: Moment estimation; MPSE: Maximum product of spacings estimator; PCE: Percentile estimator; ADE: Anderson-Darling estimator; RADE: Right-tail Anderson-Darling estimator; LSE: Least-squares estimator; WLSE: Weighted least-squares estimator; CVME: Cramer-Von-Mises estimator; MCMC: Markov chain Monte Carlo; MSE: Mean squared errors; GOF: Goodness-of-fit; AIC: Akaike information criterion; CAIC: Corrected Akaike information criterion; BIC: Bayesian information criterion; HQIC: Hannan-Quinn information criterion; D.F: Degree of freedom; EDS: Empirical descriptive statistics; KS: Kolmogorov-Smirnov statistic; P-P: Probability-Probability

1.   Introduction

In recent years, there has been remarkable progress in developing discrete probability distributions, which are essential for modeling count data across a wide range of disciplines, including statistical quality control, engineering, epidemiology, biology, agriculture, medicine, and social sciences. These distributions are used to represent diverse scenarios, such as the number of species in an ecological community, customer purchases, regional accidents, or votes received by candidates in an election. Count data, which generally represent the occurrence of random events within a specific time frame, are often viewed as the outcomes of underlying continuous-time count processes. Discretizing continuous data is frequently preferred because it simplifies analysis and interpretation in many practical applications. Such data, consisting of nonnegative integers (including zero), typically lack an upper limit, making them ideal candidates for modeling with discrete probability distributions.

While classical models like the Poisson, geometric, and negative binomial distributions have long been used for this purpose, they come with notable limitations. To address these challenges, statisticians have developed various methods for generating discrete distributions as counterparts to continuous ones, enriching the tools available for analyzing count data. These include the survival discretization approach (see Nakagawa and Osaki ^[1]), infinite series discretization technique (see Kulasekera and Tonkyn ^[2] and Sato et al. ^[3]), round-off method (see Smith ^[4]), binning method (see Roederer et al. ^[5]), minimizing Cramer-von Mises-type distances method (see Barbiero and Hitaj ^[6]), reversed hazard function discretization method (see Ghosh et al. ^[7]), compound two-phase approach (see Chakraborty ^[8]), and others. For more details (see Kotsiantis and Kanellopoulos ^[9]). These methods offer different approaches to discretizing continuous distributions, each with its own advantages and implications. For a comprehensive understanding of these discretization methods and their differences, interested readers are encouraged to refer to Chakraborty ^[8], which provides a detailed survey on the topic.

Additionally, there is another approach for proposing a discrete probability model known as the combining approach. This method involves combining two probability mass functions (PMFs) to create a new PMF, a technique commonly used in probability theory and statistics. The combining approach is often applied in situations where the outcomes of interest can be modeled by multiple distributions, aiming to produce a combined distribution that captures the characteristics of both original distributions. Various methods can be employed for combining PMFs, including convolution, mixture models, and weighted averages. Convolution, a mathematical operation, merges two functions to generate a third function. In the context of probability distributions, convolution merges the PMFs of two random variables to derive the PMF of their sum, as described in Casella and Berger ^[10], which is particularly beneficial for dealing with independent random variables. Mixture models, on the other hand, are probabilistic models representing a variable's distribution as a mixture of several component distributions. Each component distribution is weighted by a mixing coefficient, resulting in an overall distribution that is a weighted sum of these components, as outlined in Bishop ^[11]. This approach is valuable for modeling data originating from different subpopulations, each with its distinct distribution. Furthermore, weighted averages can also be employed to combine PMFs by taking the weighted averages of their values, as discussed in Gelman et al. ^[12]. This method is frequently utilized in Bayesian statistics, where prior distributions are combined with likelihood functions to yield a posterior distribution. The weighted combining approach is vital for analyzing count-type data across diverse sectors, especially in situations where accurately measuring quality characteristics on a numerical scale is challenging. This challenge is particularly pronounced in statistical quality control scenarios involving defects in manufactured items or service error frequencies. By integrating dependencies or correlations between observations within integer-valued time series data models, this approach enhances the accuracy and reliability of performance evaluation, addressing a critical aspect in count-type data analysis where the assumption of observation independence in traditional control charts may not hold. As a result, the discrete weighted PMF stands out as a potent tool that boosts the flexibility, accuracy, and interpretability of probabilistic models, making it indispensable in various analytical and modeling tasks across different domains.

Building upon these principles, the proposed paper introduces a new PMF with a single parameter, known as the hybrid discrete weighted combining (HDWC) distribution. The HDWC model is distinct within the family of statistical models due to its capacity to handle complex data behaviors that other classical distributions may not effectively represent. While the geometric, Poisson, and negative binomial distributions each have their strengths, they struggle with data exhibiting heavy tails, overdispersion, or extreme values. In contrast, the HDWC model is specifically designed to address these challenges. The geometric distribution is useful for modeling the number of trials until the first success in Bernoulli trials, where the data follows a memoryless property. However, it cannot account for overdispersion or heavy tails, limitations that the HDWC model overcomes with its ability to model skewed and heavy-tailed data. Similarly, while the Poisson distribution is widely used to model count data with a constant mean and variance, it fails to handle overdispersion or data with greater variability than the mean. The HDWC model excels in this area, as it accommodates both overdispersion and outlier-sensitive, heavy-tailed data, offering a more robust fit for real-world datasets. The negative binomial distribution extends the Poisson model to address overdispersion by introducing an additional parameter. However, it assumes a specific variance structure, which may not always match observed data. The HDWC model offers greater flexibility, effectively capturing leptokurtic shapes and data with increasing failure rates, making it especially valuable for analyzing extreme data and phenomena with high variability.

The HDWC model presents distinct advantages compared to the geometric, Poisson, and negative binomial distributions, especially in its capacity to handle data with heavy tails, overdispersion, and extreme values. Its closed-form statistical properties and versatility in addressing complex datasets make it a valuable tool for researchers in fields like engineering sustainability, actuarial analysis, and medical data. Thus, the motivation for introducing the HDWC model can be summarized as follows:

● It allows for expressing statistical properties in closed forms with a simple formulation.

● It is suitable for modeling positively skewed data with a leptokurtic shape.

● It can be used to analyze overdispersed data.

● It enables the discussion of extreme and outlier data, fitting well under heavy-tailed shapes.

● It is considered a suitable choice for researchers analyzing the phenomenon of increasing failure rates.

The paper is organized as follows: Section 2 introduces the mathematical framework of the HDWC distribution, while Section 3 explores its key statistical properties. Section 4 derives various mathematical estimation methods for practical applications. In Section 5, these methods are evaluated and compared using a ranked-simulation approach to identify optimal model parameters across different sample sizes. Section 6 demonstrates the HDWC model's flexibility in fitting real-world datasets, highlighting its performance against well-established PMFs in the statistical literature, with applications in system reliability and treatment analysis. Finally, Section 7 concludes with a summary of key findings and suggestions for future research directions.

2.   Mathematical framework of the hybrid discrete weighted combining distribution

Let $f_{1}(x; \theta) = (1-\theta)\theta^{x}$ represent a geometric distribution where the probability of success in each trial is $\theta$, and the probability of failure is $1-\theta$. On the other hand, $f_{2}(x; \theta) = (1-\theta^{2})\theta^{2\, x}$ can be viewed as a geometric distribution applied specifically to even integers $x$, rather than the conventional geometric distribution that applies to all nonnegative integers. Let $f(x; \theta) = 2\, f_{1}(x; \theta)-f_{2}(x; \theta)$, where $\, f_{1}(x; \theta)$ and $f_{2}(x; \theta)$ represent two PMFs such that $2\, f_{1}(x; \theta)-f_{2}(x; \theta)\geq0$ for all $x = 0, 1, 2, \cdots$ and $0 < \theta < 1$. Then, the PMF of the HDWC distribution can be formlated as

The conditions for a valid PMF, including nonnegativity and normalization, have been verified. Figure 1 illustrates the potential PMF shapes of the HDWC distribution, revealing unimodal right-skewed patterns with extended right tails.

The corresponding cumulative distribution function (CDF), denoted as $F(x; \theta)$, can be written as

The hazard rate function (HRF) can be listed as

The reversed hazard rate function (RHRF) can be reported as

Figure 2 presents the HRF for the HDWC model across different parameter values. The analysis indicates that the failure rate of the HDWC distribution exhibits a monotonically increasing or increasing-constant shape. A log-concave distribution has discrete increasing failure rate (IFR) (see Barlow ^[13]). The relationship between the IFR and other discrete aging classes like new better than used (NBU) can be explored in Lai and Xie ^[14].

The finding that the failure rate of the HDWC distribution shows either a consistently increasing trend or a steady, constant increase has significant implications across various domains of data analysis and analytics. In reliability engineering, this behavior aids in predicting system failures and scheduling maintenance. Healthcare analytics can leverage it to assess medical treatment effectiveness over time. Financial risk management benefits from understanding investment failure rates, while manufacturing and quality control can optimize processes based on product failure trends. Additionally, in survival analysis and event history studies, this insight proves valuable for analyzing time-to-event data in epidemiology, social sciences, and marketing, enhancing decision-making and risk mitigation strategies across diverse sectors.

3.   Statistical characteristics of the hybrid discrete weighted combining distribution

3.1.   Recurrence relation for probabilities and long-tailedness

The PMF described in Eq (2.1) exhibits log-concavity across all $\theta$ values, with the ratio of successive probabilities $\frac {f(x+1;\theta)}{f(x; \theta)}$ being a decreasing function in $x$ for all $\theta$ values. Given that for a discrete distribution, the value of $B = \lim_{x\rightarrow \infty}\frac{f(x+1;\theta)}{f(x; \theta)}$ indicates the relative long-tailedness of the distribution, and comparing it to the Poisson distribution (which has $B = 0$), we can conclude that the HDWC distribution exhibits longer tails than the Poisson distribution. Long-tailedness refers to the slow decay of probabilities in the tail of a distribution, where the probability of observing large values remains higher compared to short-tailed distributions. Unlike light-tailed distributions such as the Poisson, long-tailed distributions decay more slowly, with their tail probabilities taking longer to approach zero. This characteristic is particularly important in fields like risk assessment, finance, and insurance, where rare but extreme events can have significant consequences.

3.2.   Log-concavity

The HDWC distribution exhibits log-concavity. Specifically, since the ratio of successive probabilities is a strictly decreasing function for $x\geq0$, it follows that

This condition is sufficient for the log-concavity of a PMF. A log-concave distribution exhibits an increasing failure hazard rate, is strongly unimodal, has all moments that exist, retains its log-concavity when truncated, and its convolution with any other discrete distribution remains both unimodal and log-concave.

3.3.   The moments and associated metrics

The moment generating function (MGF) and probability generating function (PGF) are fundamental tools in probability theory and statistics with diverse practical applications. The MGF, denoted as $M_{X}(t)$, for a random variable $X$, is defined as the expected value of $e^{tX}$ where $t$ is a parameter. It is a powerful tool for finding moments of a distribution. Specifically, the $n$th moment of $X$ can be obtained by differentiating the MGF $n$ times and evaluating at $t = 0$, i.e., $E(X^{n}) = M_{X}^{(n)}(0)$. This makes the MGF useful for calculating moments without directly integrating over the distribution. The PGF, denoted as $\psi_{X}(z)$, for $X$ is defined as the expected value of $z^{X}$ where $z$ is a parameter typically restricted to the interval $[0, 1]$. It is particularly useful for discrete random variables and is employed to derive probabilities of various events related to the random variable $X$. For the HDWC distribution, the MGF and PGF can be formulated as

and

MGFs are extensively employed in statistical analysis to compute moments, crucial for determining central tendencies, variances, index of dispersion and other statistical properties of distributions, while PGFs are utilized in analyzing discrete distributions, particularly in problems involving counting and combinatorial aspects. In reliability engineering, MGFs aid in analyzing component and system lifetimes, predicting failure rates, and assessing system reliability over time. Queueing theory benefits from both MGFs and PGFs, assisting in the analysis of arrival and service times, queue lengths, and overall system performance metrics. Actuarial science utilizes MGFs for modeling insurance risks, mortality rates, and financial aspects related to insurance and pensions. In epidemiology and public health, probability generating functions play a key role in modeling disease spread, calculating infection probabilities, and evaluating the impact of interventions. In machine learning and data science, these functions are integral to probability distributions used in various models, such as the Poisson distribution for count data and the exponential distribution for survival analysis. For the HDWC distribution, the raw moments, say $E(X^{n})$, can be formulated as

Using Eq (3.3), the initial four moments can be represented in closed forms as

They are called "raw" because they are calculated directly from the data without any adjustments or transformations. The skewness, kurtosis, covariance (COV), and index of dispersion (IOD) can be calculated from the first four moments. Figure 3 illustrates several descriptive statistics pertaining to the HDWC distribution.

The data depicted in Figure 3 indicates that the HDWC model is adept at analyzing positively skewed data with a leptokurtic shape and can also be used to investigate phenomena associated with overdispersion "IOD($X$) $> 1$". Moreover, a positive covariance signifies that the variables tend to co-vary in the same direction.

3.4.   Order statistics

Order statistics refer to the statistical analysis of ranked data points within a dataset. They are particularly useful in understanding the distribution and characteristics of extreme values, such as the minimum and maximum values, as well as percentiles and quartiles. By arranging data points in ascending or descending order, order statistics allow researchers to analyze the spread and central tendencies of a dataset more effectively. For instance, in environmental studies, order statistics can help identify the most extreme weather events, such as the hottest day of the year or the coldest month, which is crucial for risk assessment, infrastructure planning, and climate change analysis. Suppose we have a random sample $X_{1}, X_{2}, \cdots, X_{n}$ of size $n$ from the CDF of the HDWC model. The CDF of the $i^{th}$ order statistic can be expressed as

where

The PMF of $i^{th}$ order statistics for $x = 1, 2, \cdots$ can be formulated as

and

By referring to Eq (3.4), it is possible to calculate the $r^{th}$ moments of order statistics, and subsequently, derive the L-moments. L-moments are statistical measures used to describe the shape, location, and scale of a probability distribution. They are particularly useful in situations where traditional moments like mean and variance may not adequately capture the distribution's characteristics, especially in cases of skewed or heavy-tailed data. L-moments are derived from linear combinations of order statistics, providing information about the distribution's central tendency, dispersion, and skewness. They are commonly used in hydrology, meteorology, and other fields where analyzing extreme values and fitting distributions to observed data is important for risk assessment and decision-making.

3.5.   Various entropies

Residual entropy and cumulative residual entropy are measures used in information theory to quantify the uncertainty or unpredictability of a random variable or system. Residual entropy captures the remaining uncertainty after observing a sequence of events, while cumulative residual entropy accounts for the accumulated uncertainty across multiple observations. These measures are particularly useful in assessing the information content or randomness in data sequences, such as in cryptography, where residual entropy helps evaluate the security of encryption algorithms based on the unpredictability of ciphertexts. The residual entropy of variable $X$ and its cumulative expression can be formulated, respectively, as follows:

and

By utilizing algebraic techniques with respect to the HDWC model, we can derive

and

On the other hand, Shannon entropy and Rényi entropy are also measures of uncertainty in information theory but differ in their sensitivity to outliers and the concentration of probability mass. Shannon entropy, the most widely used entropy measure, characterizes the average information content in a probability distribution, providing insights into the diversity and spread of data. In contrast, Rényi entropy generalizes Shannon entropy by introducing a parameter that controls the emphasis on rare events, making it more robust to outliers and suitable for analyzing distributions with heavy tails or extreme values. These entropy measures find applications in various fields such as data compression, pattern recognition, and statistical physics, where understanding the information content and structure of data is critical for modeling and decision-making. The Shanon entropy for the HDWC can be formulated as

where

The Rényi entropy of variable $X$, obtained through the use of generalized binomial expansion, is defined as

where $\gamma > 0$ and $\gamma\neq1$. Some numerical entropy of the HDWC model can be listed in Table 1.

Table 1 quantifies the uncertainty linked to its parameter combinations. Entropy increases with both $\gamma$ and $\theta$, indicating greater uncertainty. For fixed $\gamma$, entropy grows nonlinearly with $\theta$, while for fixed $\theta$, it decreases as $\gamma$ increases. Smaller $\gamma$ values, associated with heavier tails, lead to higher entropy and greater randomness, whereas larger $\gamma$ values result in a more concentrated distribution and lower entropy. Larger $\theta$ values amplify entropy, reflecting a broader and less predictable distribution, with entropy rising faster at higher $\theta$. The interplay between $\gamma$ and $\theta$ shapes the overall randomness. For example, a small $\gamma$ and large $\theta$ produce the highest entropy, while the reverse yields the lowest. These insights are valuable for applications in engineering, medicine, and actuarial sciences, where entropy helps assess the suitability of the HDWC distribution for modeling uncertainty. The relationship between entropy and parameters guides parameter selection and model fitting, leveraging the distribution's flexibility to match specific data or phenomena.

3.6.   Mean residual life and its residual coefficient of variation function

The mean residual life (MRL) and its variance function are important concepts in survival analysis, which deals with the time until an event of interest occurs. In the case of discrete random variables, such as in discrete survival analysis, these concepts are adapted to work with discrete time points. Considering $F(x; .)$ as the CDF of an element with a finite first moment, where $X$ denotes the random variable associated with $F(x; .)$, in the discrete context, the MRL, say $\varpi(k; \theta)$, is defined as follows

Consider $X$ to be an HDWC random variable, then the MRL can be formulated as

The MRL of the HDWC distribution is a decreasing function. The relationship between the HRF and MRL can be established as

Thus, the survival function (SF) of the HDWC model can be expressed as

The function for the variance of residual life (VRL), denoted as $\Upsilon_{VRL}(k)$, is defined as

The random variable $X$ exhibits increasing (decreasing) VRL if

The residual coefficient of variation (RCOV), denoted as $\Xi(k; \theta)$, can be explicitly derived as

The HRF, MRL, and VRL of the HDWC model are interconnected in the following manner

Further, the HRF, MRL, VRL, and RCOV of the HDWC distribution are interconnected as

3.7.   Stress-strength model

The stress-strength model (SSM) is a statistical framework used to assess reliability and risk in systems involving discrete random variables. It compares the strength of a system or material (represented by a random variable) against the stress it is subjected to, often in engineering, manufacturing, and quality control contexts. For example, in mechanical engineering, this model can evaluate the likelihood of a component failing under varying stress conditions. Similarly, in quality control, it helps determine if a product meets certain strength criteria given manufacturing variations. The model is also applicable in pharmaceuticals to assess the efficacy of drug formulations concerning patient-specific factors. Let $X_{1}\, \sim\ $HDWC$(\theta_{1})$ and $X_{2}\, \sim\ $HDWC$(\theta_{2})$, then

By quantifying the relationship between stress and strength through probability distributions, the SSM aids in decision-making regarding design tolerances, product specifications, and risk mitigation strategies.

3.8.   Premium principles

Insurance premiums are determined using premium principles, which consider the risk levels associated with different events. Over time, various premium principles have been developed. In this section, we introduce some of these principles, assuming that the loss distribution follows the HDWC distribution. Within this context, let $\lambda\geq0$ represent the risk loading parameter.

3.8.1.   Expected value principle

The expected value principle (EVP) is a cornerstone in insurance and risk management, guiding the calculation of insurance premiums. It asserts that the premium assessed for insurance protection must equate to the anticipated losses' expected value, modified by a risk loading component. The EVP can be articulated as

where EVP$(\lambda; .)$ represents the insurance premium, $\lambda$ is the risk loading factor, and $E(Z)$ denotes the expected value of losses or the distribution's expected value. The term $(1+\lambda)$ "Risk Loading" signifies the risk loading factor, included to cover miscellaneous expenses and ensure a profit margin for the insurer. The value of $\lambda$ varies depending on factors like administrative costs, claims processing, underwriting, and the desired profit level of the insurer. The EVP of the HDWC distribution, denoted as EVP($\lambda; \theta$), is defined by the equation:

The EVP serves as a fundamental concept in insurance pricing, applicable across diverse insurance categories such as property, liability, health, and life insurance. Insurers use the EVP to set premiums based on expected losses, adjusting for risk loading to ensure fairness in pricing. This methodology aligns premiums with the expected costs of coverage, enabling insurers to manage expenses, generate profits, and ensure that policyholders pay appropriate premiums relative to the risks transferred. From an economic standpoint, the EVP promotes risk transfer, enhances efficiency through risk pooling, and minimizes the financial impacts of uncertainty. However, it operates under the assumption of known and accurately estimated loss distributions, overlooking factors such as moral hazard and market competition that can influence pricing and market dynamics.

3.8.2.   Exponential premium principle

The exponential premium principle (EPP) is a key aspect of insurance pricing strategies, employing an exponential function to calculate premiums. This principle acknowledges that the likelihood of encountering significant loss events escalates quickly as risk levels elevate. The EPP asserts that insurance premiums should escalate exponentially in proportion to the risk level linked to the insured event. Put differently, as the probability of severe losses increases, premiums should surge at an escalating pace to sufficiently address the potential expenses associated with such occurrences. The EPP is determined by solving for EPP in the equation:

Here, $v$ refers to an individual's wealth, and $\varphi(t) = -e^{-\lambda t}$ signifies the exponential utility function. As a result, the EPP is derived as:

where $M_{Z}(\lambda)$ is the MGF. The EPP of the HDWC distribution, say EPP$(\lambda; \theta)$, is characterized as

This principle finds widespread application across different insurance sectors, such as property insurance (e.g., protection against natural disasters), liability insurance (e.g., coverage for high-risk activities), and health insurance (e.g., coverage for pre-existing medical conditions). Within these sectors, premiums are frequently determined using sophisticated actuarial models that integrate the EPP to guarantee sufficient coverage for various risks.

4.   Various estimation approaches

4.1.   Maximum likelihood estimation (MLE)

Suppose we have a random sample $X_{1}$, $X_{2}$, $\ldots$, $X_{k}$ obtained from the HDWC model. The log-likelihood, say $L(\underline{x}|\theta)$, function for the HDWC distribution can be expressed as follows:

Taking the derivative of the $L(\underline{x}|\theta)$ with respect to $\theta$ and equating it to zero, we obtain

It is not feasible to find an analytical solution for this equation. Therefore, it necessitates the use of a numerical iterative method such as the Newton-Raphson method within R software or other optimization techniques.

4.2.   Moment estimation (ME)

If we have a random sample $X_{1}$, $X_{2}$, $\ldots$, $X_{k}$ derived from the HDWC model, subsequently, the moment estimator, referred to as $\hat{\theta}_{M}$, can be formulated as

where $\overline{X}\ $represents the sample mean of the observed values. Estimation techniques like MLE and ME offer powerful tools for parameter estimation, but they come with notable computational and methodological challenges, especially when implemented in R. MLE, while efficient and consistent for large datasets, often requires iterative numerical methods such as Newton-Raphson, which can face convergence issues, instability due to singular or near-singular Hessians, and challenges with non-convex likelihood surfaces. These issues are amplified in high-dimensional models or with poor initial parameter values. ME, on the other hand, provides a simpler alternative but is less efficient, more prone to bias and variance in small samples, and heavily dependent on the choice of moments, which may not capture the data's complexity.

In R, implementing these methods demands careful handling of numerical precision, algorithmic stability, and computational efficiency. Challenges include the iterative nature of MLE, handling divergence or oscillation in numerical methods, and optimizing performance for large datasets. Users must choose starting values judiciously, apply robust optimization algorithms, and leverage advanced features in R packages like optim and bbmle to mitigate these issues. By addressing these computational hurdles and understanding the trade-offs of each method, practitioners can effectively apply MLE and ME to a wide range of statistical modeling problems while ensuring reliable and accurate results.

4.3.   Maximum product of spacings estimator (MPSE)

In this section, we explore the HDWC parameter estimation using the MPSE method with a complete sample. Suppose we have a random sample $X_{1}$, $X_{2}$, $\ldots$, $X_{k}$ obtained from the HDWC distribution. For $h = 1, 2, \ldots, m+1$, consider the following:

to be the uniform spacings of a random sample from the HDWC model, where $F\left(x_{(0)}|\theta\right) = 0$, $F\left(x_{(m+1)}|\theta\right) = 1$, and $\sum_{h = 1}^{m+1}W_{h}(\theta) = 1$. The MPSE of $\theta$, say $\widehat{\theta}_{MPS}$, can be derived by maximizing the geometric mean of the spacings

with respect to the parameter $\theta$.

4.4.   Percentile estimator (PCE)

Consider $s_{h} = h/\left(k+1\right)$ to be an unbiased estimator of $F\left(x_{(h)}|\theta\right)$. Hence, the PCE of the parameter $\theta$, denoted by $\widehat{\theta}_{PC}$, can be reported by minimizing

with respect to the parameter $\theta$, where $\varpi\left(s_{h}\right) = F^{-1}\left(x_{(h)}|\theta\right)$ is the quantile function of the HDWC model.

4.5.   Anderson–Darling and right-tail Anderson–Darling estimators

Consider a random sample $X_{1}$, $X_{2}$, $\ldots$, $X_{k}$ drawn from the HDWC model. The Anderson-Darling estimator (ADE) serves as another form of minimum distance estimator. The ADE for the HDWC parameter, denoted as $\widehat{\theta}_{AD}$, is obtained by minimizing

Regarding the parameter $\theta$, the model undergoes optimization, whereas the right-tail Anderson-Darling estimator (RADE) for the model parameter is obtained by minimizing

with respect to the parameter $\theta$.

4.6.   Weighted (least-squares) estimators

Let's consider a random sample from the HDWC model, denoted by order statistics $X_{(1)}, X_{(2)}, \cdots, X_{(k)}$. The least-squares estimator (LSE) of the HDWC parameter, represented as $\widehat{\theta}_{LS}$, is derived by solving the nonlinear equation defined as follows:

with respect to the parameter $\theta$. The weighted LSE (WLSE), say $\widehat{\theta}_{WLS}$, can be derived by solving the nonlinear equation defined by

with respect to the parameter $\theta$.

4.7.   Cramer-Von-Mises estimator (CVME)

The CVME occurs due to the difference between the estimated CDF and the empirical CDF. Determining the CVME of the HDWC parameter entails solving the nonlinear equation specified as follows:

with respect to the parameter $\theta$.

5.   Simulation ranking techniques: Different estimators

In this section, we evaluate the performance of MPSE, ME, ADE, MLE, LSE, RADE, PCE, CVME, and WLSE with respect to the sample size '$n$' using the R software and HDWC parameters. The process of generating a random variable $X$ from the HDWC distribution starts with generating a value $Y$ from the continuous distribution. Next, the obtained $Y$ value is discretized to generate $X$, where $X$ is defined as the largest integer less than or equal to $Y$. To replicate this procedure, we conduct Markov Chain Monte Carlo (MCMC) simulations using various schemes. The evaluation is conducted through a simulation study:

1) Generate $N = 10000$ samples of various sizes "$n_{i}; i = 1, 2, 3, 4, 5$" from the HDWC model as follows:

● Scheme Ⅰ: $\theta = 0.05$ $|$ $n_{1} = 20, \ n_{2} = 100, \ n_{3} = 250, \ n_{4} = 500, \ n_{5} = 750.$

● Scheme Ⅱ: $\theta = 0.2$ $|$ $n_{1} = 20, \ n_{2} = 100, \ n_{3} = 250, \ n_{4} = 500, \ n_{5} = 750.$

● Scheme Ⅲ: $\theta = 0.5$ $|$ $n_{1} = 20, \ n_{2} = 100, \ n_{3} = 250, \ n_{4} = 500, \ n_{5} = 750.$

● Scheme Ⅳ: $\theta = 0.8$ $|$ $n_{1} = 20, \ n_{2} = 100, \ n_{3} = 250, \ n_{4} = 500, \ n_{5} = 750.$

2) Compute the MPSE, ME, ADE, MLE, LSE, RADE, PCE, CVME, and WLSE for the 10,000 samples, say $\widehat{\theta}_{k}$ for $k = 1, 2, ...,$ 10,000.

3) Caculate the bias and the mean squared errors (MSE) for $N =$ 10,000 samples as

The MSE quantifies the average squared deviation between predicted and actual values, where a lower MSE indicates closer predictions to the actual values.

4) The empirical results of simulation are listed in the Tables 2–6.

As the sample size '$n$' increases, the trend observed in Tables 2 to 6 shows a clear reduction in the bias of the parameter $\theta$, approaching zero, while the MSE of the HDWC parameter also decreases toward zero. The findings indicate that the derived estimators consistently improve in performance as the sample size increases. Across varying sample sizes, all estimation methods perform well, with Table 6 demonstrating that the ME method outperforms MLE in simulations based on a ranking approach. This ranking process evaluates the methods against criteria such as bias, variance, MSE, efficiency, consistency, and robustness. Simulated data are generated, the estimation methods are applied, and their metrics are compared to determine their effectiveness. The final ranking balances these criteria, providing insights into the most suitable method for the problem at hand.

The choice between MLE and ME ultimately depends on the specific context, including factors like data size, computational complexity, and model requirements. Each method offers distinct advantages: MLE is celebrated for its efficiency, consistency, and precision, particularly in large datasets, making it a powerful choice when computational resources and well-defined likelihood functions are available. On the other hand, ME is a simpler, computationally lighter alternative, particularly advantageous when the likelihood function is difficult to derive or when computational simplicity is a priority. Selecting the appropriate method requires a careful balance between accuracy, simplicity, and the resources available for analysis.

6.   Data analysis

In this section, we showcase the versatility of the HDWC distribution in accurately modeling a variety of datasets from different domains. We assess the fitted distributions using the Chi-square ($\chi^{2}$) test and its corresponding P-value in dataset Ⅰ. In datasets Ⅱ and Ⅲ, the tested models are compared using the Kolmogorov-Smirnov (KS) test with its corresponding P-value and some criteria, namely, the maximized log-likelihood ($l$), Akaike information criterion (AIC), corrected AIC (CAIC), Bayesian IC (BIC), and Hannan-Quinn IC (HQIC). Since there is a limited number of frequencies for each observation in datasets Ⅱ and Ⅲ, the Pearson's Chi-square statistic cannot be employed for an inference test. Therefore, the KS measure is adequate in this case. The codes can be listed in the following link: https://docs.google.com/document/d/13XpIsgrCGryY9Ma1777NQ46fDpz32Ufc/edit. We will compare the performance of the HDWC distribution against other competing distributions listed in Table 7.

6.1.   Dataset Ⅰ

The dataset provided represents the count of computer breakdowns occurring over a period of 128 consecutive weeks of operation (Hand et al. ^[21]). Nonparametric plots are suitable for the initial visualization of this dataset, and the corresponding results can be shown in Figure 4. It is observed that there are some extreme values, displaying a right-skewed distribution. Table 8 presents the MLEs for the parameter(s), $-l$, and Chi-square test values, accompanied by their corresponding P-values, for all competing distributions in dataset Ⅰ.

Of all the distributions tested, the HDWC distribution shows the best performance as the most suitable distribution for fitting dataset Ⅰ. The estimated PMFs of all the distributions tested for dataset Ⅰ are illustrated in Figure 5. From this figure, the same results are those derived from the Table 8.

Table 9 lists different estimators for dataset Ⅰ, and it was noted that the MLE and ME are the best estimators for the dataset Ⅰ. In addition, the CVME, LSE, and WLSE techniques work quite well for modeling this data.

Table 10 presents numerical summaries of the empirical descriptive statistics (EDS).

The characteristics of dataset Ⅰ clearly show a right-skewed distribution with leptokurtic features and overdispersion.

6.2.   Dataset Ⅱ

The second dataset comprises the failure times for a sample of 15 electronic components subjected to an accelerated life test (refer to Lawless ^[22]). Nonparametric plots for dataset Ⅱ are illustrated in Figure 6. It is noted that there are no extreme values, yet the distribution remains right-skewed.

The MLEs and goodness of fit (GOF) measures for this data are reported in Table 11.

As observed, at a significance level of $0.05$, the HDWC, NBI, Geo, DINH, and DR models (for full names, refer to Table 7) perform effectively in modeling dataset Ⅱ, with the HDWC distribution being the most flexible. Figure 7 displays the probability-probability (P-P) plots for dataset Ⅱ, which support the empirical findings presented in Table 11.

Table 12 lists various estimators for dataset Ⅱ, and it was noted that all techniques work quite well for modeling this data except the RADE.

The empirical CDFs plots for dataset Ⅱ for all competing distributions are given in Figure 8 (left panel) and the empirical CDFs plots for dataset Ⅱ for the HDWC using the different estimation methods are reported in Figure 8 (right panel) which prove the empirical results mentioned in Tables 11 and 12.

Table 13 presents numerical summaries of the EDS.

Dataset Ⅱ clearly exhibits a right-skewed distribution with leptokurtic characteristics and overdispersion.

6.3.   Dataset Ⅲ

This dataset includes leukemia remission times (in weeks) for 20 patients, as discussed by Damien and Walker ^[23] in relation to the concept of discretization. Nonparametric plots for dataset Ⅲ are displayed in Figure 9.

The MLEs and the GOF test and measures for this data are reported in Table 14.

As can be seen, based on significance level $0.05$, all tested distributions work quite well for analyzing dataset Ⅲ except DIR, DBH, DPa, and Poi, but the HDWC distribution is the best (For full names, refer to Table 7). Figure 10 presents the P-P plots for dataset Ⅲ, reinforcing the empirical results shown in Table 14.

Table 15 provides a list of various estimators for dataset Ⅲ, noting that all methods perform flexibility in modeling this data.

Figure 11 (left panel) displays the empirical CDF plots for dataset Ⅲ across all competing distributions, while Figure 11 (right panel) presents the empirical CDF plots for the HDWC distribution using various estimation methods. Together, these figures support the empirical findings discussed in Tables 14 and 15.

Table 16 provides numerical summaries of the EDS for dataset Ⅲ.

Dataset Ⅲ distinctly displays a right-skewed distribution with leptokurtic traits and overdispersion.

6.4.   Reliability analysis of systems and treatments: Insights from HRF, MRL, RHRF, and VRL trends

The analysis of the three datasets provides insights into their reliability characteristics from different perspectives. Dataset Ⅰ, representing the count of computer breakdowns over 128 consecutive weeks, reveals an increasing HRF, indicating a growing likelihood of breakdowns as time progresses, likely due to aging or wear-and-tear effects. The MRL decreases as more breakdowns occur, suggesting shorter expected operational times before subsequent failures. Similarly, the RHRF and VRL also decrease, reinforcing the notion of declining system reliability over time and highlighting the importance of proactive maintenance strategies. Dataset Ⅱ, involving failure times for 15 electronic components subjected to accelerated life tests, also shows an increasing HRF, suggesting higher failure probabilities as time advances under stress conditions. The decreasing MRL highlights the reduced lifespan of components as failures occur, while the declining RHRF and VRL confirm the ongoing degradation. These observations emphasize the necessity of robust design, quality control, and testing to enhance component reliability in challenging environments. Dataset Ⅲ, which contains leukemia remission times for 20 patients, demonstrates a similar reliability trend. The increasing HRF suggests that the risk of remission ending grows over time, potentially due to disease progression or diminishing treatment effectiveness. The decreasing MRL indicates shorter expected remission durations, while the declining RHRF and VRL highlight the reduced stability of remission periods. These findings underscore the importance of tailored treatment plans and continuous monitoring to improve patient outcomes. Figure 12 visually summarizes these findings, demonstrating consistent patterns across the datasets, all of which reflect progressively worsening reliability over time, calling for targeted strategies like preventive maintenance, rigorous testing, and personalized interventions to enhance performance and outcomes.

The reliability analysis of these datasets provides valuable insights into the progression of failures and the potential for improving performance and resilience. By understanding these patterns, researchers and practitioners can develop effective approaches tailored to the unique challenges presented by engineering systems, electronic components, and medical treatment plans.

7.   Conclusions and future directions

The article introduced a discrete distribution with a single parameter, developed using a weighted combining discretization approach, referred to as the HDWC distribution. This model offers a comprehensive set of statistical properties, all derived and expressed in closed forms. It has proven to be particularly effective for modeling datasets that exhibit right-skewness and leptokurtic shapes, making it well-suited for real-world data with heavy tails and extreme values. The HDWC distribution is particularly valuable in situations where the hazard rate function increases, especially when dealing with outliers, providing a more robust framework compared to traditional models. The article explored multiple estimation methods for the HDWC model, including MPSE, ME, ADE, MLE, LSE, RADE, PCE, CVME, and WLSE, and found that all methods were successful in estimating the model's parameter. Among these, the ME approach demonstrated the best performance. To further validate the effectiveness of the HDWC distribution, three real-world datasets were analyzed, showing that the HDWC distribution outperformed other competing distributions, such as the geometric, Poisson, and negative binomial models, in all aspects of the analysis. The article also outlined several promising directions for future research. These include extending the model to handle time-dependent stochastic processes for applications such as survival analysis and event prediction, adapting the model for spatially correlated data in fields like epidemiology and environmental studies, and incorporating Bayesian methods to improve parameter estimation, especially for small datasets. Additionally, the potential for combining the HDWC model with machine learning techniques or exploring hybrid models could further enhance its flexibility and predictive power, opening up new avenues for its application in complex data modeling.

Use of AI tools declaration

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

Acknowledgments

The authors extend their appreciation to Prince Sattam bin Abdulaziz University for funding this research work through the project number (PSAU/2024/01/31935).

Data Availability

The data that supports the findings of this study are available within the article.

Conflict of interest

There are no conflicts of interest.