A novel two-parameter discrete exponentiated exponential distribution with its neutrosophic extension: mathematical characterization and real-world count data modeling

1 Introduction

The analysis of count data (non-negative integer-valued observations) arises across a broad spectrum of scientific and applied disciplines, including epidemiology, reliability engineering, actuarial science, ecology, and operational research. Typical examples include the number of equipment failures within a maintenance cycle, hospital visits per patient per year, disease outbreak counts in public health surveillance, and defective units in manufacturing inspection. Conventional discrete probability distributions such as the Poisson, binomial, geometric, and negative binomial serve as the classical foundation for modeling such data; however, they impose restrictive structural assumptions that limit their applicability in practice.

The discretization of continuous lifetime distributions has emerged as a principled and increasingly popular methodology for constructing new discrete distributions that inherit the theoretical richness of their continuous counterparts while remaining suited to count data. Among the available discretization techniques (including mixed Poisson, infinite series, and moment-matching approaches) the survival function discretization method, introduced by Nakagawa and Osaki (1975), has attracted the greatest attention due to its mathematical tractability, its preservation of essential reliability properties at integer support points, and the interpretability of the resulting PMF as the probability of discrete “time-to-event” outcomes.

On the basis of this framework, researchers have introduced numerous discrete distributions through the systematic discretization of established continuous models, including discrete Weibull (Nakagawa and Osaki, 1975), discrete normal (Roy, 2003), discrete Rayleigh (Roy, 2004), discrete half-normal (Kemp, 2008), discrete Burr XII (Krishna and Pundir, 2009), discrete Lindley (Gómez-Déniz and Caldeŕın-Ojeda, 2011), discrete log-logistic (Para and Jan, 2016), discrete Xgamma (Maiti et al., 2018), discrete Burr-Hatke (El-Morshedy et al., 2020), discrete natural Lindley (Al-Babtain et al., 2020), discrete Marshall-Olkin Weibull (Opone et al., 2021), discrete power-Ailamujia (Alghamdi et al., 2022), discrete Ramos-Louzada (Eldeeb et al., 2022), discrete Bilal (Altun et al., 2022), discrete XLindley (Eldeeb et al., 2024), discrete Hjorth (Haj Ahmad and Elshahhat, 2025), and discrete Muth (Elsayed and Hussein, 2025), reflecting ongoing innovations in discretization theory. Despite such advancements, a persistent challenge in count data modeling remains the joint accommodation of three empirically common features: (i) flexible dispersion structure, encompassing both over-dispersion and under-dispersion; (ii) flexible hazard rate shape, allowing both increasing and decreasing failure rates governed by a single interpretable parameter; and (iii) tractable closed-form expressions for reliability quantities such as the mean residual life and the reverse hazard rate. Most existing discrete distributions address only a subset of these requirements. The Poisson distribution is structurally equidispersed and exhibits a constant-like hazard behavior, rendering it unsuitable for over- or under-dispersed data. The geometric and discrete exponential distributions accommodate only monotone decreasing hazard rates. The negative binomial, while handling over-dispersion, cannot accommodate under-dispersion. Several two-parameter discrete distributions, such as the discrete Burr XII and the discrete power-Ailamujia, offer shape flexibility but do not provide explicit closed-form expressions for the mean residual life or reverse hazard rate, and the theoretical connection between the dispersion transition and the shape parameter in these models remains implicit rather than analytically transparent.

The continuous exponentiated exponential (EE) distribution, introduced by Gupta and Kundu (1999) and further studied by Gupta and Kundu (2001), addresses analogous modeling challenges in the continuous setting: the shape parameter α governs both the hazard rate shape and the tail behavior of the distribution. Its discrete analog has not been studied in the literature, which motivates the present work.

The broader family of exponential-type lifetime models has continued to expand in recent years, reflecting sustained interest in flexible reliability and lifetime distributions for engineering and medical applications. Recent contributions include the Marshall–Olkin extended exponential distribution for engineering and failure-time data (Ragab et al., 2025a), a further extension of the extended exponential model validated on medical and failure-time datasets (Ragab et al., 2025b), the generalized Dinesh–Umesh–Sanjay generalized exponential distribution for engineering data (Alrashidi and Ragab, 2023), a three-parameter exponential model applied to engineering and medical data (Almutiry et al., 2021), the DUS-transformed generalized linear failure rate distribution (Elgarhy et al., 2025), a discrete Burr III distribution developed for radiation and biological count data (Alomair et al., 2025b), and a unit extended exponential distribution for bounded proportion-type data (Ragab et al., 2024). These studies underscore the continued relevance of exponential-based lifetime models across application domains, both continuous and discrete, and motivate the present development of a discrete counterpart with comparable flexibility for count-data settings.

Rationale. Despite this active body of work, no existing discrete distribution combines, within a single interpretable shape parameter, the ability to switch between increasing and decreasing hazard rates and between over- and under-dispersion, while retaining fully closed-form expressions for every standard reliability quantity. This is precisely the gap the present paper addresses: by discretizing the exponentiated exponential distribution, we obtain a two-parameter count model in which the shape parameter α governs both properties simultaneously, without sacrificing analytical tractability. The remainder of this Introduction develops the discretization framework and the neutrosophic preliminaries needed for the extension in Section 8, before Section 2 formally introduces the DEE distribution.

The EE distribution was selected as the baseline for discretization, rather than alternative continuous lifetime models (e.g., gamma, Weibull, or generalized Rayleigh), for three specific reasons. First, the EE survival function admits a fully elementary closed-form expression, S(x; α, λ) = 1−(1−e^−λ*x*)α, involving no special functions, unlike the gamma or generalized Rayleigh families, whose survival functions require the incomplete gamma or incomplete beta function. Applying the survival discretization technique of Equation 3 to such an elementary form guarantees that every reliability quantity of the resulting discrete model, namely the PMF, CDF, hazard rate, reverse hazard rate, moment generating function, and mean residual life, remains expressible in closed algebraic form, as demonstrated throughout Sections 2 and 3. This is a genuine tractability advantage over discretized gamma- or Weibull-type models, whose discrete counterparts typically require numerical evaluation of the underlying continuous CDF at every integer support point rather than an explicit algebraic formula.

Second, the EE distribution nests the exponential distribution as the special case α = 1, so the single shape parameter α has a transparent interpretation as a departure from the memoryless (constant-hazard) baseline: values α < 1 and α>1 correspond to interpretable, and oppositely signed, deviations in hazard shape. This interpretability carries over directly to the discrete setting, where the fitted α^ can be read diagnostically from data, as illustrated in Section 6.

Third, and most important for the present application, the EE shape parameter α simultaneously governs both the hazard rate shape and the dispersion regime of the discretized model (see Section 2 and Remark 5, a dual interpretive role that is not shared by discretizations of the gamma or Weibull families, where the shape and scale parameters typically influence hazard shape and dispersion jointly and non-monotonically rather than through a single interpretable parameter. These three considerations, closed-form tractability, transparent parameter interpretation, and dual hazard/dispersion control, jointly motivate the choice of the EE distribution as the discretization baseline for this work.

The specific contributions of this paper are as follows:

• We introduce the Discrete Exponentiated Exponential (DEE) distribution and establish closed-form expressions for all fundamental distributional quantities, including the PMF, CDF, survival and hazard rate functions, reverse hazard rate function, moment generating and characteristic functions, raw and central moments up to fourth order, and the mean residual life function.
• We demonstrate analytically, via the Glaser technique, that the DEE hazard rate function is monotonically increasing for 0 < α < 1 and monotonically decreasing for α>1, and show numerically that the dispersion index transitions continuously from over-dispersion to under-dispersion as α increases, a dual flexibility governed by a single parameter that is not shared by most competing models.
• We conduct a Monte Carlo simulation study that not only assesses finite-sample MLE performance but also benchmarks the DEE estimator behavior against those of four competing discrete distributions under identical data-generating conditions.
• We validate the DEE distribution through three real-world count datasets from oncology, education, and nephrology.
• We develop a neutrosophic extension of the DEE distribution (DNEE) for modeling count data subject to epistemic uncertainty and indeterminacy, with accompanying numerical illustrations of the neutrosophic mean, variance, and PMF under varying indeterminacy.

Gupta and Kundu (2001) proposed a flexible two-parameter lifetime distribution as a generalization of the exponential distribution, subsequently recognized as the exponentiated exponential (EE) distribution (see also Gupta and Kundu, 1999). A random variable X is said to follow the EE distribution if its probability density function (PDF) and cumulative distribution function (CDF) are expressed as

f(x;α,λ)=αλe-λx(1-e-λx)α-1, x>0,(1)

and

F(x;α,λ)=(1-e-λx)α, x≥0,λ>0,α>0.(2)

The survival discretization approach is employed when a continuous random variable (r.v.) Y is defined on R+. In accordance with the methodology outlined by Nakagawa and Osaki (1975), if the discrete random variable X has a survival function Y, then the probability mass function (PMF) of the discrete random variable X is:

P(X=x)=S(x)-S(x+1), x=0,1,2,…(3)

This formulation preserves the survival probabilities at integer points, ensuring that the discrete distribution retains essential reliability properties of its continuous analog.

1.1 Neutrosophic data analysis

The neutrosophic statistical framework was first proposed by Smarandache (1998) in 1998, being a generalized philosophical approach that encompasses both fuzzy logic and intuitionistic fuzzy logic in its theory construct. The neutrosophic probability measure was conceptualized by Smarandache (2013) as a function NP:Y → [0, 1]³ operating on a neutrosophic sample space Y, with its mapping defined by NP(S)=(ch(S),ch(neut(S)),ch(anti(S)))=(δ,ρ,S), where δ, ρ, and S denote the degrees of truth, indeterminacy, and falsity, respectively, satisfying 0≤δ,ρ,S≤1 and 0≤δ+ρ+S≤3. Scholars have developed neutrosophic variants for a wide range of distributions, such as exponential, binomial, Poisson, uniform, Lognormal, Kumaraswamy, Weibull, normal, generalized Pareto, and gamma see (Ahsan-ul Haq¹, 2022; Alhabib et al., 2018; Khan et al., 2021b; Alhasan and Smarandache, 2019; Khan et al., 2022, 2021a; Patro and Smarandache, 2016). The concept of neutrosophic random variables has been formally defined by Zeina and Hatip (2021), who have formulated basic principles, while Granados (Granados and Sanabria, 2021; Granados et al., 2021; Granados, 2022) has further analyzed the mathematical properties, independence, and discrete distributions with neutrosophic parameters.

Recent studies have shown the effectiveness of the framework in goodness-of-fit tests, neutrosophic time series models for predictions, and other modeling techniques such as logarithmic models, moving average models, and linear models (Alhabib and Salama, 2020; Alhabib¹ and Salama, 2020; Cruzaty et al., 2020). Ahsan-ul Haq and Zafar (2023) generalized a one-parameter discrete Ramos-Louzada distribution to its neutrosophic counterpart and showed its applicability in modeling count data with indeterminacy. Further applications in survival analysis include the neutrosophic extension of the Topp-Leone extended exponential distribution for characterizing uncertain lifetime data in medical contexts (Hammood et al., 2025). The above-mentioned studies reflect the increasing need for neutrosophic models in dealing with real-world data that is prone to vagueness, incompleteness, and imprecision.

1.2 Preliminaries

This subsection establishes fundamental notations and concepts that form the basis for our proposed distribution. We denote by S the sample space, ℝ the set of real numbers, and Φ a sample space event. The symbols X_*N* and Y_*N* represent neutrosophic random variables. We present key definitions and properties from neutrosophic probability theory that are essential for constructing the neutrosophic extension of our discrete model.

Definition 1. Let X be a real-valued classical random variable defined on a probability space (S,F,P) with the mapping:

X:S→ℝ

where S represents the event space. A neutrosophic random variable X_*N* is then constructed as:

XN:S→ℝ(I)

and expressed as

XN=X+I

*where I denotes the indeterminacy component satisfying 0 ≤ I ≤ 1.*

Theorem 1. *(See Granados, 2022) Assume a neutrosophic random variable X_*N* = X+I, where F_*X*(x) = P(X ≤ x) represents the CDF of the classical random variable X. The following relationships are established:*

FXN(x)=FX(x-I),

fXN(x)=fX(x-I)

*where F_{X_*N*} and f_{X_*N*} denote the CDF and PDF of the neutrosophic random variable X_*N*, respectively.*

Theorem 2. *(See Granados, 2022) Assume X_*N* = X+I represents a neutrosophic random variable. The corresponding expected value and variance satisfy:*

E(XN)=E(X)+I and Var(XN)=Var(X).

Corollary 1. *If X has moment generating function MX(t)=E(etX), then the neutrosophic random variable X_*N* = X+I has moment generating function*

MXN(t)=etIMX(t).

Remark 1. An important property of the theoremsAn important property of the theorems is that, even as the neutrosophic extension shifts the distribution by indeterminacy parameter I, it retains the dispersion characteristics of the original distribution. This property proves to be advantageous in the analysis of reliability, as the variability structure of the distribution must be retained.

The main aim of this work is to develop new theoretical results and a neutrosophic extension for the discrete probability distribution obtained by applying the survival discretization approach to the exponentiated exponential distribution, hereafter referred to as the “Discrete Exponentiated Exponential Distribution-DEE” for expository convenience [see the remark following Equation 3 in Section 2 for its correspondence with existing literature]. Building on this distributional form, the model admits closed-form expressions for its probability mass function (PMF), cumulative distribution function (CDF), survival function, hazard rate function, moments, and reliability functions. The parameters are estimated using the maximum likelihood estimation (MLE) method, and the performance of the method is evaluated using a simulation study that includes competitor benchmarking. The usefulness and competitiveness of the DEE-based framework developed here are shown using analysis of real-world datasets from diverse application domains. Moreover, to model count data with uncertainty and indeterminacy, a neutrosophic version of the model is developed, accompanied by numerical illustrations.

2 Discrete exponentiated exponential distribution

Employing the survival discretization technique specified in Equation 3, the probability mass function of the proposed model is derived by applying Equation 3 to the EE density and CDF (Equations 1–2) as

P(x;γ,α)=(1-γ(x+1))α-(1-γx)α, x∈ℕ0,γ>0,α>0(4)

where γ = e^−λ and ℕ₀ = {0, 1, 2, …, w} for 0 < w < ∞.

To explore the shape behavior of the DEE distribution, visual depictions of its probability mass function are rendered in Figure 1 (Equation 4) for distinct choices of the parameters γ and α.

The associated cumulative distribution function and survival function are given by:

F(x;γ,α)=(1-γ(x+1))α, x∈ℕ0(5)

and

S(x;γ,α)=1-(1-γ(x+1))α, x∈ℕ0(6)

Remark 2. We note that the PMF and CDF above (Equations 4 and 5) coincide, under the parameter correspondence γ≡p, with the discrete generalized exponential distribution of the second type, DGE2(α, p), introduced by Nekoukhou et al. (2013). We derived this model independently via the survival-function discretization of the continuous exponentiated exponential distribution of Gupta and Kundu (1999); on comparison with the existing literature, this discretization is found to correspond to DGE2(α, p), and we cite Nekoukhou et al. (2013) accordingly. The present work extends the study of this distribution through several new theoretical results, including a strict log-concavity theorem for the PMF (Section 3.2), an explicit uniqueness proof for the maximum likelihood estimator of α (Section 4), closed-form expressions for the minimum and maximum order statistics (Section 3.5), a reparameterized and more interpretable form of the reverse hazard rate function (Section 2), a comparative simulation study benchmarking estimator performance against several competing discrete models, validation on three real-world datasets spanning education, nephrology, and oncology, and a novel neutrosophic extension (Section 8) that, to the best of our knowledge, has not previously been developed for this distribution.

Using the PMF (Equation 4) and survival function (Equation 6) above, the hazard rate function of the DEE is

h(x;γ,α)=P(x)S(x)=(1-γ(x+1))α-(1-γx)α1-(1-γ(x+1))α, x∈ℕ0(7)

The reverse hazard rate is also known as the past failure rate. This is an important reliability measure that calculates the conditional probability of a failure occurring just before a specified time. This reliability measure is particularly useful when the interest is in retrospective failure analysis. This is particularly useful in post-failure analysis for industrial systems as well as retrospective risk analysis for clinical studies. The reverse hazard rate is defined as the ratio of the probability mass function to the cumulative distribution function:

rh(x;γ,α)=1-(1-γx1-γ(x+1))α, x∈ℕ0(8)

Equation 8 is already fully closed-form and elementary: it involves a single ratio raised to the power α,with no special function or unevaluated series. Its interpretation is made more transparent by writing rh(x; γ, α) = 1−[1−ψ(x)]^α, where

ψ(x)=γx-γx+11-γx+1=γx(1-γ)1-γx+1, x∈ℕ0,(9)

is the elementary proportional decrement between the consecutive discretization terms 1−γ^*x* and 1−γ^*x*+1. In this reparameterization, ψ(x)∈(0, 1) decreases monotonically in x toward 0, and the reverse hazard rate is simply 1 minus (1−ψ(x))^α: the shape parameter α acts as a single amplifier of the elementary decrement ψ(x), exactly mirroring the role α plays in the forward hazard rate function h(x; γ, α). This reparameterization (Equation 9) shows that rh(x; γ, α) is governed by the same interpretable discretization mechanism, γ↔γ^*x*+1, as the forward hazard rate, without introducing any additional notation beyond the model's own parameters γ and α.

The hazard rate behavior is analyzed using the Glaser technique (Glaser, 1980) in conjunction with the PMF of the DEE distribution, applied to the hazard rate function (Equation 7) defined above.

ρ(x)=-p′(x)p(x)=α(-ln γ)γx(1-γx)α-1-γ(x+1)(1-γ(x+1))α-1(1-γ(x+1))α-(1-γx)α,

hence

ρ′(x)=-α(α-1)(ln γ)2γx(1-γx)α-2(1-γx+1)α-(1-γx)α

As ρ′(x)>0 for 0 < α < 1 and ρ′(x) < 0 for α>1, the hazard rate function of the DEE distribution is an increasing function of x when 0 < α < 1 and a decreasing function of x when α>1. Figure 2 presents a graphical portrayal of the hazard rate function of the DEE distribution, highlighting its various shape characteristics.

Remark 3. The parameter α thus serves as a shape-regime switch: values 0 < α < 1 correspond to distributions with increasing hazard rates, appropriate for modeling components that experience wear-out or deterioration over time (e.g., aging electronics, biological tissues); values α>1 correspond to decreasing hazard rates, suitable for systems that exhibit early-life failures or infant mortality (e.g., burn-in periods in semiconductor manufacturing). This interpretable role of α as a hazard-shape controller is a key practical advantage of the DEE distribution over one-parameter discrete models, which fix the hazard behavior a priori.

3 Statistical properties of the DEE distribution

This section obtains and discusses the key mathematical properties of the DEE distribution. The closed-form expressions for the moment generating function, characteristic function, and reliability measures are derived. The properties are further investigated numerically for specific values of the parameters.

3.1 Moments and generating function

Let X be a random variable following the DEE distribution. The moment generating function, denoted by M_*X*(t), is derived as follows:

MX(t)=E(etX)=∑x=0∞etxP(X=x) =∑x=0∞etx(1-γ(x+1))α-∑x=0∞etx(1-γx)α =∑x=0∞etx(1-γ(x+1))α-et∑x=0∞etx(1-γ(x+1))α =(1-et)∑x=0∞etx(1-γ(x+1))αMX(t)=(1-et)∑k=0∞(-1)k(αk)γk1-etγk, |etγk|<1.

The r^th moment about the origin is obtained by:

μr′=drdtrMX(t)|t=0

Thus, the first four raw moments of the DEE distribution are expressed as:

μ1′=∑k=0∞(-1)k(αk)γk(1-γk)2,μ2′=∑k=0∞(-1)k(αk)γk(1+γk)(1-γk)3,μ3′=∑k=0∞(-1)k(αk)γk(1+4γk+γ2k)(1-γk)4,μ4′=∑k=0∞(-1)k(αk)γk(1+11γk+11γ2k+γ3k)(1-γk)5.

The variance of the distribution is computed as

σ2=∑k=0∞(-1)k(αk)γk(1+γk)(1-γk)3-(∑k=0∞(-1)k(αk)γk(1-γk)2)2.

Remark 4. Although the raw moments μ1′,μ2′ above are given as infinite series in k rather than in fully reduced closed form, their combination in σ2=μ2′-(μ1′)2 is guaranteed to be positive for every admissible (γ, α), 0 < γ < 1, α>0. This follows directly from the general definition σ2=Var(X)=E[(X-μ1′)2]≥0, which holds for every non-degenerate random variable regardless of the specific closed-form expression used to evaluate μ1′ and μ2′; equality is attained only in the degenerate limit and does not occur for any finite (γ, α) in the admissible parameter space, since the DEE distribution has support on at least two integers for all such (γ, α). This structural guarantee is confirmed numerically in Table 1, where the Variance column is strictly positive across the entire range of (α, γ) values considered.

Table 1. Numerical characteristics of the DEE distribution for selected (α, γ) values.

The moments about the mean of the DEE distribution are given by:

μ=∑k=0∞(-1)k(αk)γk(1-γk)2,μ2=∑k=0∞(-1)k(αk)γk(1-γk+γ2k)(1-γk)4,μ3=∑k=0∞(-1)k(αk)γk(1-3γk+3γ2k-γ3k)(1-γk)5,μ4=∑k=0∞(-1)k(αk)γk(1-4γk+6γ2k-4γ3k+γ4k)(1-γk)6.

Table 1 reports key distributional measures such as the mean, variance (Var), dispersion index (DI), skewness (CS), kurtosis (CK), and coefficient of variation (CV).

Table 1 presents descriptive statistics of the DEE distribution for selected values of (α, γ). It is observed that as α increases for a fixed γ, the mean rises gradually, while the variance grows more slowly, causing the dispersion index (DI) to decline. For small values of α, DI exceeds unity, indicating overdispersion commonly observed in high-variability count data, such as accident frequencies or rare disease occurrences. For larger values of α, DI falls below unity, reflecting underdispersion typical of controlled processes.

Meanwhile, both skewness (CS) and kurtosis (CK) decrease substantially with increasing α. The declining skewness suggests greater symmetry with shorter right tails, whereas the decreasing kurtosis indicates lighter tails and a flatter peak, characteristics of more evenly distributed counts rather than data concentrated at low values with occasional extreme observations.

Remark 5. A key insight from Table 1 is the dispersion transition induced by α: as α increases for a fixed γ, the dispersion index crosses unity from above, transitioning the distribution from over-dispersed to under-dispersed. This transition occurs continuously and is governed solely by α, enabling the analyst to select an appropriate regime from the data, a property not available in the Poisson (equidispersion by construction), the negative binomial (over-dispersion only), or most one-parameter discrete lifetime models. The resulting analytical and practical flexibility makes the DEE distribution particularly well-suited for count data arising in contexts where the dispersion regime is unknown a priori or varies across subpopulations.

3.2 Log-concavity of the probability mass function

A discrete probability mass function P(x) on ℕ₀ is said to be log-concave if P(x)²≥P(x−1)P(x+1) for every x≥1, equivalently, if the ratio P(x+1)/P(x) is non-increasing in x. Log-concavity is a strong structural property, and the following theorem establishes it for the DEE distribution.

Theorem 3. *For α≥1, the probability mass function P(x; γ, α) of the DEE distribution is strictly log-concave in x on ℕ₀.*

Proof. By construction, P(x;γ,α)=S(x)-S(x+1)=∫xx+1f(t;α,λ)dt, where f(t; α, λ) = αλe^−λ*t*(1−e^−λ*t*)^α−1 is the continuous EE density underlying the discretization. Differentiating ln f(t; α, λ) twice with respect to t shows that f is log-concave in t whenever α≥1, the well-established log-concavity regime of the EE density (Gupta and Kundu, 2001). It is a standard preservation result for interval probabilities of a log-concave density (An, 1998) that if f is log-concave on ℝ₊, then the sequence P(x)=∫xx+1f(t)dt, x ∈ ℕ₀, is itself log-concave. Since the DEE PMF is exactly of this interval-probability form, the result follows.

Remark 6. Strict log-concavity of P(x; γ, α) for α≥1 carries several practically useful implications. First, it guarantees that P(x; γ, α) is strongly unimodal, ruling out multimodal or irregular PMF shapes and ensuring a single, well-defined mode, which aids interpretability in applied count-data modeling. Second, log-concave count distributions have tails that decay at least geometrically, consistent with the closed-form existence of the moment generating function established below. Third, the class of log-concave discrete distributions is closed under convolution: if X₁, X₂ are independent log-concave count variables, then X₁+X₂ is also log-concave. This is directly relevant to the integer-valued autoregressive (INAR) time-series extension of the DEE distribution proposed as future work in Section 9, since it guarantees that aggregated or superposed DEE counts (for α≥1) retain a well-behaved, unimodal structure. Fourth, log-concavity of the underlying likelihood contributions is known to promote numerical stability of Newton-Raphson-type optimization (Section 4), consistent with the fast empirical convergence of the MLEs reported in Section 5. For 0 < α < 1, the underlying continuous EE density is not globally log-concave, so Theorem 3 is stated for the α≥1 regime; the increasing-hazard regime 0 < α < 1 is instead characterized analytically via the Glaser technique above.

3.3 Characteristic function

The characteristic function is a very important concept in probability theory that uniquely represents a probability distribution and helps to find moments by differentiating the function. The characteristic function is widely used in proving convergence in distribution, investigating the sum of independent random variables, and finding limiting distributions in statistical inference. The characteristic function of the DEE distribution, denoted by ϕ_*X*(t) is obtained as follows:

ϕX(t)=E(eitX)=∑x=0∞eitxP(X=x) =∑x=0∞eitx(1-γx+1)α-∑x=0∞eitx(1-γx)α =∑x=0∞eitx(1-γx+1)α-eit∑x=0∞eitx(1-γx+1)α =(1-eit)∑x=0∞eitx(1-γx+1)αϕX(t)=(1-eit)∑k=0∞(-1)k(αk)γk1-eitγk, |eitγk|<1.

3.4 Mean residual life function

The mean residual life (MRL) function is a significant reliability measure that captures the expected remaining lifetime of a system given that it has survived up to time t, and is instrumental in evaluating maintenance strategies, burn-in testing procedures, and warranty optimization. For the DEE distribution, the MRL function simplifies to

MRL=ε(i)=E(X-i∣X≥i)=11-F(i-1;γ,α)∑j=i+1S[1-F(j-1;γ,α)], i∈ℕ0,

where ℕ0={0,1,2,…,S} and 0<S<∞. The explicit form of the MRL function for the DEE distribution is

MRL=∑j=i+1S[1-(1-γj)α]1-(1-γi)α, i∈ℕ0,0<γ<1,α>0.

Remark 7. The closed-form tractability of the MRL function for the DEE distribution is practically significant. In warranty analysis and maintenance planning, the MRL quantifies the expected remaining service life of a component given that it has already operated without failure for i time units. For the DEE distribution, the MRL can be computed directly from the fitted parameters γ^ and α^ without numerical integration, unlike many competing distributions. This makes the DEE model directly applicable to post-installation component monitoring and maintenance scheduling in reliability-critical systems.

3.5 Order statistics

Let X₁, …, X_*n* be a random sample of size n from the DEE distribution with CDF F(x; γ, α) and survival function S(x; γ, α), and let X₍₁₎ ≤ ⋯ ≤ X_(*n*) denote the corresponding order statistics. The minimum and maximum order statistics admit substantially simplified closed-form expressions relative to a general r-th order statistic, because the sample size n enters only as an integer power, in contrast with the real-valued shape parameter α.

Maximum. The CDF of X_(*n*) is F(n)(x)=[F(x;γ,α)]n=(1-γx+1)nα, which is already fully closed-form (it is simply the DEE CDF with shape parameter nα in place of α). Consequently, the PMF of the maximum is

P(X(n)=x)=(1-γx+1)nα-(1-γx)nα, x∈ℕ0,(10)

requiring no series expansion at all.

Minimum. The survival function of X₍₁₎ is S(1)(x)=[S(x;γ,α)]n=[1-(1-γx+1)α]n. Since n is a positive integer, the binomial theorem gives the finite expansion

S(1)(x)=∑j=0n(-1)j(nj)(1-γx+1)jα, x∈ℕ0,(11)

a sum of exactly n+1 elementary terms rather than an infinite series, and the PMF of the minimum follows as P(X₍₁₎ = x) = S₍₁₎(x−1)−S₍₁₎(x), likewise a finite sum of n+1 terms.

General order statistic. For a general r-th order statistic (1 < r<n), the PMF involves (nr)-weighted combinations of F(x)^*k*[1−F(x)]^*n*−*k* over k≥r; because F(x; γ, α) = (1−γ^*x*+1)^α is itself already a non-integer power of an elementary term, expanding [1−F(x)]^*n*−*r* via the binomial theorem, as in Equation 11, again produces a finite sum in the integer exponent (n−r), but each resulting term still contains F(x)^*k* = (1−γ^*x*+1)^*kα* raised to the real power α, which does not reduce further. Expressions of this type for the general order statistic are therefore most compactly left in this finite-binomial-sum form rather than an infinite series; no further closed-form reduction is available because α is real-valued.

Remark 8. Unlike the raw and central moments of Section 3.1, which genuinely require infinite hypergeometric-type series in the binomial coefficient (αk) because α enters as a non-integer exponent applied to a variable base, the minimum and maximum order statistics in Equations 10–11 require no infinite series: the sample size n is a positive integer, so all binomial expansions in n terminate after n+1 terms. Consequently, there is no truncation error to control for the extreme order statistics; Equation 11 is exact for any finite n, and its evaluation cost grows linearly, O(n), in the sample size. For the general r-th order statistic, the only infinite-series component that can arise is inherited from the underlying moment or generating-function expressions of Section 3.1 if higher moments of X_(*r*) are required; in that case, the geometric decay of the terms γ^*k* in the binomial series (αk)γk (with 0 < γ < 1) ensures fast convergence, and truncating the series at K terms incurs a remainder that decays geometrically in K, so that K on the order of a few dozen terms is sufficient for numerical precision to machine tolerance in typical applications.

4 Parameter estimation

In this section, we describe the estimation of the parameters of the DEE distribution using the maximum likelihood estimation (MLE) method. A simulation study is also conducted to compare the efficiency of the MLE for different parameter values and sample sizes. Let x₁, x₂, …, x_*n* be a random sample of size n from the DEE distribution with parameters γ and α. The log-likelihood function is given by

ℓ(γ,α;x)=∑i=1nlog[(1-γxi+1)α-(1-γxi)α](12)

The likelihood equations are obtained by differentiating Equation 12 with respect to γ and α, yielding the following nonlinear equations.

∂ℓ∂γ=∑i=1n(13)

αγxi-1[xi(1-γxi)α-1-(xi+1)γ(1-γxi+1)α-1](1-γxi+1)α-(1-γxi)α=0

and

∂ℓ∂α=∑i=1n(14)

(1-γxi+1)αlog(1-γxi+1)-(1-γxi)αlog(1-γxi)(1-γxi+1)α-(1-γxi)α=0

Theorem 4. *For fixed γ∈(0, 1) and fixed data x₁, …, x_*n*, the score (Equation 14) has at most one solution α>0; that is, whenever a maximum likelihood estimate α^ exists for fixed γ, it is unique.*

Proof. Write ui=1-γxi+1 and vi=1-γxi, so that 0 < v_*i*<u_*i* < 1 for every i (since γ∈(0, 1) and x_*i*≥0). Equation 14 can be written as g(α)=∑i=1nφ(α;ui,vi)=0, where, for fixed 0 < v<u < 1,

φ(α;u,v)=uαln u-vαln vuα-vα.

It suffices to show that φ(α; u, v) is strictly monotonic in α for every fixed pair u≠v, since g is then a sum of functions sharing the same direction of monotonicity, hence itself strictly monotonic, and a strictly monotonic continuous function admits at most one zero. Let a = u^α>0, b = v^α>0, p = ln u < 0, q = ln v < 0, with p≠q since u≠v. Using da/dα = ap and db/dα = bq, direct differentiation gives

∂φ∂α=(ap2-bq2)(a-b)-(ap-bq)2(a-b)2.

Expanding the numerator,

(ap2-bq2)(a-b)-(ap-bq)2=-ab(p-q)2,

so that

∂φ∂α(α;u,v)=-uαvα(ln u-ln v)2(uα-vα)2<0(15)

for every α>0, since u^α, v^α>0 and ln u≠ln v whenever u≠v. Inequality (Equation 15) shows φ(α; u_*i*, v_*i*) is strictly decreasing in α for every observation i, so g(α)=∑i=1nφ(α;ui,vi) is strictly decreasing in α on (0, ∞). A strictly decreasing continuous function crosses zero at most once, establishing the uniqueness of the root α^ whenever one exists.

Remark 9. Inequality (Equation 15) shows that the strict monotonicity of the score function stems entirely from the term (ln u−ln v)²>0, which is strictly positive whenever the consecutive discretization terms u = 1−γ^*x*+1 and v = 1−γ^*x* differ, a condition satisfied for every count observation x≥0 given γ∈(0, 1). This confirms that the fast, stable convergence of the Newton-Raphson iterations reported in Section 3 is not incidental but is guaranteed by the strict monotonicity of the profile score function in α implied by Theorem 4.

In the absence of an analytical solution to the equation under consideration, numerical solutions are derived using an iterative method. The maximum likelihood estimates (MLEs) of the parameters γ and α are obtained by solving Equations 13, 14 iteratively, wherein the Newton-Raphson method is employed as the numerical optimization technique. All related computations are performed using R software (R Core Team, 2020).

5 Simulation analysis

In this section, a simulation study is performed to assess the finite-sample behavior of the MLEs of the DEE parameters γ and α, and crucially, to benchmark the DEE estimator performance against those of four established competing discrete distributions: the Discrete Burr XII (DBXII), the Discrete Natural Lindley (DNL), the Discrete Power-Ailamujia (DPA), and the Discrete Ramos-Louzada (DRL). The competitor models were selected to represent both one-parameter and two-parameter families and a range of hazard-rate behaviors, providing a comprehensive baseline for comparison.

Using the inverse transformation method, random samples are generated from the DEE distribution for n = 25, 50, 100, 200, 300 with N = 10, 000 replications. The following four parameter combinations are examined: (γ = 0.3, α = 0.5), (γ = 0.3, α = 2.5), (γ = 0.5, α = 0.5), and (γ = 0.5, α = 1.5). Let ϑ∈{γ, α} denote either parameter generically. The performance of the MLEs is evaluated using the average estimate (AVE), average bias (AB), mean square error (MSE), and mean relative error (MRE), defined as

AVE(ϑ)=1N∑i=1Nϑ^i, AB(ϑ)=1N∑i=1N|ϑ^i-ϑ|,MSE(ϑ)=1N∑i=1N(ϑ^i-ϑ)2, MRE(ϑ)=1N∑i=1N|ϑ^i-ϑ|ϑ.

The DEE simulation results are reported in Table 2 and the corresponding heatmaps are presented in Figure 3. From Table 2, the following conclusions are drawn:

• The AVE values of γ^ and α^ approach their respective true values as n increases across all cases, confirming the consistency of the MLEs.
• The AB, MSE, and MRE of γ^ and α^ decrease monotonically as n increases in all cases, establishing the asymptotic efficiency of the MLEs.
• The parameter γ is estimated with greater precision than α across all cases and sample sizes, as evidenced by consistently smaller AB and MSE values for γ^.

Table 2. Simulated performance measures of the MLEs for the parameters γ and α.

Competitor benchmarking. To assess whether the DEE MLE converges faster and with lower MSE than analogous estimators for competing models, we additionally generated 10,000 samples of sizes n = 25, 50, 100, 200, 300 from the DBXII, DNL, DPA, and DRL distributions (at their true parameter values), and computed the MLE MSE for the respective scale/shape parameters under each competitor model. Table 3 reports the MSE of the primary scale parameter estimator for each model at n = 25 and n = 100, representative of small and moderate samples. The DEE MLEs exhibit competitive convergence rates, with the MSE of γ^ declining substantially faster than those of the one-parameter competitors (DNL and DRL) and comparably to the two-parameter competitors (DBXII and DPA). Notably, for the over-dispersed parameter setting (γ = 0.5, α = 0.5), the DEE MLE achieves smaller MSE than the DRL scale estimator even at n = 25, demonstrating practical estimation efficiency despite the two-parameter structure of the DEE model.

Table 3. Comparative MSE of the primary scale parameter estimator across DEE and four competing models at representative sample sizes (N = 10, 000 replications).

Parameters for each competitor set to values matching similar mean behavior as DEE Case III (γ = 0.5, α = 0.5).

Remark 10. The larger absolute MSE of γ^ in the DEE model (relative to one-parameter competitors) is expected: the DEE simultaneously estimates two parameters from the same sample. Normalizing by the number of parameters and the parameter range, the DEE estimation efficiency is broadly comparable to two-parameter competitors such as DBXII and DPA. The key advantage of the DEE is that the two parameters serve orthogonal interpretive roles (scale and dispersion/hazard regime), providing diagnostically richer estimates without a proportionate loss of precision.

6 Application

This section evaluates the practical performance of the DEE distribution on three real-world count datasets from distinct application domains: education, nephrology, and oncology. In each case, the DEE model is compared against seven established discrete distributions: Poisson (Poi), discrete Burr XII (DBXII), discrete Natural Lindley (DNL), discrete Burr-Hatke (DBH), discrete power-Ailamujia (DPA), discrete Ramos-Louzada (DRL), and discrete Poisson-Haq (DPH) (Alomair et al., 2025a). The probability mass functions of these competing models are summarized in Table 4.

Table 4. Probability mass functions of competing models.

Model adequacy is assessed using four complementary criteria: the maximum log-likelihood value (ℓ), the Akaike information criterion (AIC), the Bayesian information criterion (BIC), and a χ² goodness-of-fit statistic with its associated p-value. The χ² statistic is computed by grouping each dataset's support into categories chosen to keep the expected frequency in every category numerically stable (following standard grouping practice for count data, Cochran, 1954), comparing the observed count in each category against the count implied by each fitted model, and evaluating the statistic against a χ² distribution with degrees of freedom equal to the number of categories minus one minus the number of estimated parameters. Lower AIC and BIC, together with higher log-likelihood and χ² p-values, indicate superior model fit. For each dataset, we additionally present a table comparing the observed frequency against the frequency expected under every fitted model, giving a direct, transparent view of fit quality across the full range of the data.

Dataset I (mathematics exam scores): Final mathematics examination scores for n = 48 slow-learner students at the Indian Institute of Technology at Kanpur (Bakouch et al., 2014), a genuine non-negative integer count dataset. Descriptive statistics: mean = 25.8958, variance = 346.1379, skewness = 1.3317, kurtosis = 4.3233, indicating a positively skewed, heavy-tailed distribution with substantial over-dispersion (DI = 13.3665>1). Exploratory plots appear in Figure 4; MLEs and model comparison measures are in Table 5; fitted PMFs are shown in Figure 5.

Table 5. Parameter estimates with standard errors and goodness-of-fit statistics for Dataset I (mathematics exam scores, n = 48).

Bold values indicate the DEE model, which achieves the lowest χ² statistic and highest p-value among all eight fitted models for this dataset.

Dataset II (kidney cyst counts): Counts of kidney cysts in patients treated with corticosteroids, reported by Chan et al. (2010) and later analyzed by Al-Bossly et al. (2023). The dataset comprises n = 110 observations ranging from 0 to 11; the observed frequency distribution is given in Table 6. Descriptive statistics: mean = 1.3909, variance = 6.1118, skewness = 2.2612, kurtosis = 7.8066, confirming notable over-dispersion (DI = 4.3941>1). Exploratory plots appear in Figure 6; MLEs and model comparison measures are in Table 7; fitted PMFs are shown in Figure 7.

Table 6. Observed frequency distribution of kidney cyst counts (n = 110).

Table 7. Parameter estimates with standard errors and goodness-of-fit statistics for Dataset II (kidney cyst counts, n = 110).

Bold values indicate the DEE model, which achieves the lowest χ² statistic and highest p-value among all eight fitted models for this dataset.

Dataset III (breast cancer tumor size): Tumor size measurements (in mm) of female breast cancer patients diagnosed in the United States, obtained from a breast cancer registry (Namdari, 2023). A random sample of n = 50 observations was drawn using R software; the variable of interest is the recorded primary tumor dimension at diagnosis, transformed by (x−1) to conform to the model's non-negative integer support. Descriptive statistics: mean = 26.0200, variance = 240.9996, skewness = 1.6773, kurtosis = 6.3769, indicating a positively skewed, heavy-tailed distribution with substantial over-dispersion (DI = 9.2621>1). Exploratory plots appear in Figure 8; MLEs and model comparison measures are in Table 8; fitted PMFs are shown in Figure 9.

Table 8. Parameter estimates with standard errors and goodness-of-fit statistics for Dataset III (breast cancer tumor size, n = 50).

Bold values indicate the DEE model, which achieves the lowest χ² statistic and highest p-value among all eight fitted models for this dataset.

Remark 11. The three datasets were chosen to span distinct application domains and over-dispersion levels: Dataset I (education) has a known over-dispersed structure and has been used previously to benchmark discrete distributions (Bakouch et al., 2014); Dataset II (nephrology) is a canonical over-dispersed event-count dataset; and Dataset III (oncology) contributes a different measurement type, an anatomical dimension recorded at diagnosis rather than a naturally occurring event count, together with a different dispersion level. In every case the variable of interest is a genuine non-negative integer count (Dataset III after a location shift to conform to the model's support), ensuring that the DEE distribution is validated on data types appropriate to its intended use.

Remark 12. For Dataset I, the χ² test groups the support into class intervals of width 10 (0–9, 10–19, 20–29, 30–39, ≥40). Owing to the wide range and heavy right tail of the exam-score data (from 4 to 86), the Poisson and DBH models cannot accommodate the substantial over-dispersion present and yield very large χ² values, correctly reflecting their inadequacy. The DEE distribution achieves the lowest χ² statistic and highest p-value among all eight models (χ² = 1.423, p = 0.4910), confirming its superior fit, and Table 9 presents the corresponding observed-vs.-expected frequency comparison.

Table 9. Observed vs. expected frequencies for all fitted models, Dataset I (mathematics exam scores, n = 48).

Categories are grouped in class intervals of width 10, with the final category pooling x≥40. Bold values indicate the expected frequencies under the DEE model, shown for direct comparison with the observed counts.

Remark 13. For Dataset II, the χ² test groups the support into the categories x = 0, 1, …, 7, ≥8, chosen so that the observed count in every category is at least 2. The DEE distribution provides the best fit among all eight models (χ² = 0.778, p = 0.9926), and Table 10 presents the corresponding observed-vs.-expected frequency comparison.

Table 10. Observed vs. expected frequencies for all fitted models, Dataset II (kidney cyst counts, n = 110).

The last category pools x≥8 to keep expected counts numerically stable, following standard practice for discrete goodness-of-fit testing (Cochran, 1954). Bold values indicate the expected frequencies under the DEE model, shown for direct comparison with the observed counts.

Remark 14. The near-unity χ² p-value for Dataset II (p = 0.9926) arises because the DEE's estimated α^=0.3828<1 correctly captures the increasing-hazard structure of kidney cyst counts: patients with more cysts are at progressively higher conditional risk of additional cysts, a clinically plausible pattern that most competing models cannot represent without structural modification. The DBXII (p = 0.6905) and DPA (p = 0.9450) achieve competitive but lower p-values due to their inability to parameterically encode an increasing-hazard shape. The DBH achieves p = 0.8892 but this model is structurally restricted to decreasing-hazard behavior; its fit here reflects a broadly adequate envelope rather than a mechanistically correct model. The DEE's advantage in this dataset is therefore structural, not merely numerical.

Remark 15. As with Dataset I, the χ² test for Dataset III groups the support into class intervals of width 10 (0–9, 10–19, 20–29, 30–39, ≥40). The DEE distribution again achieves the lowest χ² statistic and highest p-value among all eight fitted models (χ² = 7.834, p = 0.0199), together with the lowest AIC and BIC, confirming its comparative superiority. The absolute p-value is more modest here than for Datasets I and II, reflecting the smaller sample size (n = 50) and the wide, heavy-tailed spread of clinical tumor size measurements; nonetheless, DEE remains the clear best-fitting model by a substantial margin over every competitor, as Table 11 illustrates.

Table 11. Observed vs. expected frequencies for all fitted models, Dataset III (breast cancer tumor size, n = 50).

Categories are grouped in class intervals of width 10, with the final category pooling x≥40. Bold values indicate the expected frequencies under the DEE model, shown for direct comparison with the observed counts.

7 Discussion

The empirical results across all three datasets offer insights that go beyond the goodness-of-fit statistics and merit further discussion. In each application, the fitted parameters carry interpretable diagnostic content. For Dataset I (mathematics examination scores), γ^=0.0663 and α^=2.687; α^>1, pointing to a decreasing hazard structure, which reflects the fact that the probability of a student scoring extremely high decreases as the score threshold increases within a population of slow learners. For Dataset II (kidney cyst counts), γ^=0.2883 and α^=0.3828; here α^<1, indicating an increasing hazard rate, meaning patients with more cysts are at progressively higher conditional risk of exhibiting additional cysts, a clinically plausible pattern in corticosteroid-induced kidney dysmorphogenesis. For Dataset III (breast cancer tumor size), γ^=0.9208 and α^=4.4763; α^>1 indicates a decreasing hazard rate, consistent with the clinical expectation that very large primary tumor dimensions are progressively rarer at the time of diagnosis. The ability of the same distributional family to accommodate both hazard regimes across different datasets without model modification is a key practical advantage of the DEE framework.

The competitive fit of the DEE distribution across all three datasets can be attributed to its parametric efficiency: two parameters that jointly control both the location/scale (γ) and the shape/dispersion regime (α). Several competing models, namely the Poisson, DNL, DRL, and DPH, are one-parameter families and therefore cannot adjust the dispersion index to match the data, leading to systematic misfit in over-dispersed samples. The two-parameter competitors (DBXII, DPA) offer shape flexibility but do not exhibit the same analytically transparent separation between dispersion-increasing and dispersion-decreasing regimes governed by a single parameter. In Dataset II, where the DEE's estimated α^<1 correctly identifies the increasing-hazard structure of the cyst count data, this distinction is particularly apparent: competing models achieve reasonable but inferior fit because they cannot represent this structure without modification.

From the theoretical and empirical evidence presented, the DEE distribution is most suitable when the data exhibit over-dispersion (DI>1), which is the case in all three datasets analyzed here, or when the context implies a monotone hazard structure that can be diagnosed empirically by inspecting whether α^≷1. The DEE distribution is less appropriate for data with non-monotone hazard rate shapes, such as bathtub or upside-down bathtub forms, as the current formulation does not accommodate such behavior; this limitation identifies a natural direction for future work.

8 Neutrosophic version of DEE distribution

Let Y_*N* denote a neutrosophic random variable. We say that Y_*N* follows a neutrosophic discrete Exponentiated Exponential distribution, denoted by Y_*N*~DEE(γ_*N*, α_*N*), where the neutrosophic parameters γ_*N* and α_*N* are characterized as sets containing one or multiple elements, typically represented as intervals to capture uncertainty or imprecision in the data. The neutrosophic probability mass function (pmf) of Y_*N* is defined as

fY(y−I)={(1−γN y−I+1)αN−(1−γN y−I)αN,y=I,1+I,2+I,…,0,otherwise,

where I represents the indeterminacy component inherent in neutrosophic theory, allowing the distribution to model datasets with vague, incomplete, or uncertain observations.

Proof

We need to verify that the neutrosophic PMF sums to unity:

∑y=I∞P(y)=∑y=I∞[(1-γNy-I+1)αN-(1-γNy-I)αN]

Let z = y−I. When y = I, we obtain z = 0, and as y → ∞, it follows that z → ∞. Therefore,

∑y=I∞P(y)=∑z=0∞[(1-γNz+1)αN-(1-γNz)αN].

This is a telescoping series. By the telescoping property and taking the limit as n → ∞, we obtain:

∑y=I∞P(y)=limn→∞[(1-γNn+1)αN-(1-γN0)αN].

Since 0 < γ_*N* < 1, we have limn→∞γNn+1=0. Also, (1-γN0)αN=(1-1)αN=0. Therefore,

∑y=I∞P(y)=(1-0)αN-0=1.

By Theorem 2, the expressions for the mean and variance are given by:

E(YN)={∑k=0∞(-1)k(αNk)γNk(1-γNk)2}+I,Var(YN)=∑k=0∞(-1)k(αNk)γNk(1+γNk)(1-γNk)3- (∑k=0∞(-1)k(αNk)γNk(1-γNk)2)2.

These series converge rapidly for typical parameter values with 0 < γ_*N* < 1 and α_*N*>0.

8.1 Neutrosophic numerical illustration

To demonstrate the practical behavior of the DNEE distribution under varying degrees of indeterminacy, Table 12 reports the neutrosophic mean E(Y_*N*), variance Var(Y_*N*), and dispersion index DI(Y_*N*) = Var(Y_*N*)/E(Y_*N*) for selected parameter combinations (γ_*N*, α_*N*) and indeterminacy levels I∈{0, 0.25, 0.50, 0.75, 1.00}.

Table 12. Neutrosophic mean, variance, and dispersion index of the DNEE distribution for selected (γ_*N*, α_*N*) and indeterminacy level I.

Table 12 reveals several practically important patterns. First, as predicted by Theorem 1.2.3, Var(Y_*N*) remains unchanged across all levels of I: the indeterminacy component shifts the mean upward without inflating the spread. Second, because E(Y_*N*) increases with I while Var(Y_*N*) is fixed, the dispersion index DI(Y_*N*) decreases monotonically as I grows. This means that higher indeterminacy effectively moves the DNEE distribution toward under-dispersion. Third, for the parameter setting (γ_*N* = 0.3, α_*N* = 0.5), the DNEE transitions from over-dispersion (DI = 2.000 at I = 0) to under-dispersion (DI = 0.250 at I = 1) purely through the indeterminacy component, providing a mechanism for modeling ambiguous dispersion regimes in datasets with indeterminate observations. These properties confirm that the DNEE framework extends the DEE's flexibility to uncertain environments while preserving the distributional integrity guaranteed by the classical model.

9 Conclusion

This paper introduces and fully characterizes the Discrete Exponentiated Exponential (DEE) distribution, a new two-parameter discrete lifetime model derived via the survival discretization of the continuous exponentiated exponential distribution. The principal theoretical contribution is the analytical demonstration that the single shape parameter α simultaneously controls two fundamental modeling features: the hazard rate shape (increasing for 0 < α < 1, decreasing for α>1) and the dispersion regime (over-dispersion for small α, under-dispersion for large α). This dual parametric role is established through the Glaser monotonicity analysis and the numerical investigation of the dispersion index, providing a theoretically grounded explanation for the distributional flexibility of the DEE model, a flexibility that competing one-parameter discrete distributions such as the Poisson, geometric, and discrete natural Lindley cannot match without structural modification.

Closed-form expressions are derived for all standard distributional quantities, including the PMF, CDF, survival and hazard rate functions, reverse hazard rate, moment generating and characteristic functions, raw and central moments up to fourth order, and the mean residual life function. The tractability of these expressions, particularly the closed-form MRL, facilitates direct application in reliability and maintenance analysis without numerical integration.

Model parameters are estimated via maximum likelihood, with numerical optimization using the Newton-Raphson method. A comprehensive simulation study (N = 10, 000 replications across five sample sizes) confirms consistency and asymptotic efficiency of the MLEs, and includes an explicit benchmarking of DEE MLE convergence rates against those of four competing discrete models (DBXII, DNL, DPA, DRL) under identical data-generating conditions, demonstrating that the DEE achieves competitive estimation efficiency despite its two-parameter structure.

The practical utility of the DEE distribution is demonstrated through three real-world count datasets spanning mathematics exam scores, kidney cyst counts, and breast cancer tumor size measurements. In all three applications, the DEE model performs competitively or superiorly against seven discrete distributions on all goodness-of-fit criteria. Crucially, the fitted parameter values carry interpretable content: datasets with α^>1 exhibit decreasing hazard rates consistent with their subject-matter context, while the kidney cyst dataset with α^<1 exhibits an increasing hazard rate reflecting cumulative count risk, demonstrating that the DEE framework provides not merely statistical fit, but diagnostically meaningful parameter estimates.

Finally, the neutrosophic generalization of the DEE distribution (DNEE) extends the model's applicability to count data characterized by vagueness, incompleteness, or epistemic indeterminacy. A new numerical illustration of the DNEE demonstrates how the indeterminacy component I affects the mean and dispersion index while preserving the variance, providing practitioners with concrete guidance on interpreting DNEE parameter estimates in uncertain data environments.

Despite these contributions, the present study has several limitations that should be acknowledged. First, the DEE distribution is a two-parameter model, and while the shape parameter α jointly governs hazard shape and dispersion, the model does not accommodate non-monotone hazard shapes (e.g., bathtub or upside-down bathtub), which arise in some reliability and survival applications. Second, the real-data validation, while spanning three distinct domains, is limited to three datasets of moderate size (n = 48, 110, and 50); the smaller Dataset III yielded a more modest (though still best-in-class) goodness-of-fit p-value, suggesting that estimation precision and fit quality may be sensitive to sample size in small-to-moderate samples. Third, parameter estimation in this study relies exclusively on maximum likelihood; the finite-sample behavior of alternative estimators (e.g., method of moments or Bayesian estimators with informative priors) has not been investigated. Fourth, the neutrosophic extension currently treats the indeterminacy parameter I as fixed rather than estimated from data, which limits its immediate applicability to datasets where the degree of indeterminacy is unknown a priori.

Future research directions include: count regression models incorporating the DEE as the response distribution; integer-valued autoregressive (INAR) time series formulations based on the DEE; zero-inflated and hurdle variants for excess-zero count data; Bayesian estimation frameworks with informative priors for reliability applications; bivariate and multivariate extensions; extensions to accommodate non-monotone (bathtub or upside-down bathtub) hazard rate shapes; and estimation procedures that treat the neutrosophic indeterminacy parameter I as a data-driven quantity rather than a fixed input.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Vikas Chaurasiya, National Forensic Sciences University, India

Reviewed by: Ibrahim E. Ragab, Alexandria University, Egypt

Fatma Zohra Seghier, Higher Normal School of Technological Education, Algeria