Maryedith Uche Anozie^1*, Ben Ifeanyichukwu Oruh², Chukwuemeka Onwuzurike Omekara³, Samuel Ugochukwu Enogwe^3*
1Department of Mathematics Education, Federal College of Education (Technical), Umunze, Nigeria.
²Department of Mathematics, Michael Okpara University of Agriculture, Umudike, Nigeria.
³Department of Statistics, Michael Okpara University of Agriculture, Umudike, Nigeria.
DOI: 10.4236/ajor.2026.161001   PDF    HTML   XML   84 Downloads   257 Views   Citations

Abstract

The reliance on exponential assumptions in classical queueing models often leads to a misrepresentation of real-world systems where service and interarrival times exhibit more complex variability. This paper introduces the Nova distribution, a novel one-parameter lifetime distribution developed as a mixture of exponential and gamma components. We derive its fundamental statistical properties and demonstrate its applicability by modeling interarrival and service times in a G/G/1 queueing system. Using real-world data from a banking service facility, we show that the Nova distribution provides a superior fit compared to established one-parameter models like Lindley and Shanker, as measured by Akaike and Bayesian Information Criteria. By integrating the Nova distribution into a G/G/1 framework, we derive essential performance measures and an associated economic cost model. The results confirm that the Nova-based queueing model offers a more accurate and cost-effective tool for analyzing and optimizing service systems, particularly under the high-utilization conditions typical in many real-world scenarios.

Keywords

Queueing Theory, G/G/1 Queue, Nova Distribution, Mixture Distribution, Maximum Likelihood Estimation, Performance Measures

Share and Cite:

1. Introduction

Queueing theory is vital for analyzing and optimizing service systems. The accuracy of any queueing model hinges on how well its underlying probability distributions represent real-world interarrival and service times. While the M/M/1 queue is tractable, its exponential assumptions are often violated in practice [1] [2]. The more general G/G/1 queue accommodates arbitrary distributions but requires a judicious choice of model to balance flexibility and complexity. This has spurred interest in one-parameter mixture distributions, such as the Lindley and Shanker distributions, which offer more flexibility than the exponential distribution [3] [4].

This study proposes the Nova distribution, a new one-parameter mixture distribution designed to better capture the empirical characteristics of service and arrival processes. The objectives are to:

1) Derive the statistical properties of the Nova distribution.

2) Empirically validate its superiority over existing one-parameter distributions.

3) Formulate a G/G/1 queueing model with Nova-distributed interarrival and service times.

4) Derive key performance measures and a cost model for system optimization.

2. The Nova Distribution

2.1. Probability Density Function (PDF) and Derivation

The Nova distribution is a mixture of an exponential distribution and a gamma distribution with a shape parameter of 2.

Theorem 1: A random variable $X$ is said to have a Nova distribution with parameter $\lambda >0$ , denoted by $X~Nova\left(\lambda \right)$ , if its probability density function (PDF) is given by:

$f\left(x;\lambda \right)=\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x},x>0,\lambda >0$ (1)

Proof: The Nova distribution can be derived as a mixture of exponential and gamma distributions. Let $f_1\left(x\right)=\lambda {\text{e}}^{-\lambda x}$ be the PDF of an exponential distribution with parameter $\lambda$ , and $f_2\left(x\right)=\frac{{\lambda }^3}{2}{x}^2{\text{e}}^{-\lambda x}$ be the PDF of a gamma distribution with shape parameter 4 and rate parameter $\lambda$ .

Define the mixing proportion $p=\frac{{\lambda }^2+1}{{\lambda }^2+6}$ . Then, the mixture distribution is given by:

$\begin{array}{c}f\left(x;\lambda \right)=p{f}_1\left(x\right)+\left(1-p\right)f_2\left(x\right)\\ =\frac{{\lambda }^2+1}{{\lambda }^2+6}\left(\lambda {\text{e}}^{-\lambda x}\right)+\left(1-\frac{{\lambda }^2+1}{{\lambda }^2+6}\right)\left(\frac{{\lambda }^3}{2}{x}^2{\text{e}}^{-\lambda x}\right)\\ =\frac{{\lambda }^2+1}{{\lambda }^2+6}\left(\lambda {\text{e}}^{-\lambda x}\right)+\left(\frac{5}{{\lambda }^2+6}\right)\left(\frac{{\lambda }^3}{2}{x}^2{\text{e}}^{-\lambda x}\right)\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x},x>0,\lambda >0\end{array}$

which is the PDF given in (1) and completes the proof of Theorem 1.

Corollary 1: The PDF of the one-parameter Nova distribution is a proper density function.

Proof: To be a proper density function, the PDF in (1) must satisfy the following two conditions:

1) $f\left(x;\lambda \right)\ge 0$ for all $x>0$ .

2) $\begin{array}{c}{\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x;\lambda \right)\text{d}x}={\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\text{d}x}\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right){\displaystyle \underset{0}{\overset{\infty }{\int }}{\text{e}}^{-\lambda x}\text{d}x}+\frac{5{\lambda }^2}{2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^2{\text{e}}^{-\lambda x}\text{d}x}\right]=1\end{array}$ .

Hence, the proposed distribution is a proper density function.

2.2. Cumulative Distribution Function of the Nova Distribution

Theorem 2: The cumulative distribution function (CDF) of $X~Nova\left(\lambda \right)$ is given by:

$F\left(x;\lambda \right)=1-\left[1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda x},x>0,\lambda >0$ (2)

Proof:

$\begin{array}{c}F\left(x;\lambda \right)={\displaystyle \underset{0}{\overset{x}{\int }}f\left(t;\lambda \right)\text{d}t}=\frac{\lambda }{{\lambda }^2+6}{\displaystyle \underset{0}{\overset{x}{\int }}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{t}^2\right]{\text{e}}^{-\lambda t}\text{d}t}\\ =\frac{\lambda }{{\lambda }^2+6}{\displaystyle \underset{0}{\overset{x}{\int }}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{t}^2\right]{\text{e}}^{-\lambda t}\text{d}t}\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right){\displaystyle \underset{0}{\overset{x}{\int }}{\text{e}}^{-\lambda t}\text{d}t}+\frac{5{\lambda }^2}{2}{\displaystyle \underset{0}{\overset{x}{\int }}{t}^2{\text{e}}^{-\lambda t}\text{d}t}\right]\end{array}$ (3)

But

${\displaystyle \underset{0}{\overset{x}{\int }}{\text{e}}^{-\lambda t}\text{d}t}={\frac{{\text{e}}^{-\lambda t}}{-\lambda }|}_0^x={\frac{{\text{e}}^{-\lambda t}}{\lambda }|}_x^0=\frac{1}{\lambda }\left({\text{e}}^{-0}-{\text{e}}^{-\lambda x}\right)=\frac{1}{\lambda }\left(1-{\text{e}}^{-\lambda x}\right)$ (4)

and

$\begin{array}{c}{\displaystyle \underset{0}{\overset{x}{\int }}{t}^2{\text{e}}^{-\lambda t}\text{d}t}=uv-{\displaystyle \underset{0}{\overset{x}{\int }}v\text{d}u}=t^2{\frac{{\text{e}}^{-\lambda t}}{-\lambda }|}_0^x-{\displaystyle \underset{0}{\overset{x}{\int }}\frac{{\text{e}}^{-\lambda t}}{-\lambda }\left(2t\text{d}t\right)}\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{2}{\lambda }{\displaystyle \underset{0}{\overset{x}{\int }}t{\text{e}}^{-\lambda t}\text{d}t}\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{2}{\lambda }\left[t{\frac{{\text{e}}^{-\lambda t}}{-\lambda }|}_0^x-{\displaystyle \underset{0}{\overset{x}{\int }}\frac{{\text{e}}^{-\lambda t}}{-\lambda }\text{d}t}\right]\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{2}{\lambda }\left[\frac{-x{\text{e}}^{-\lambda x}}{\lambda }+\frac{1}{\lambda }{\displaystyle \underset{0}{\overset{x}{\int }}{\text{e}}^{-\lambda t}\text{d}t}\right]\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{2}{\lambda }\left[\frac{-x{\text{e}}^{-\lambda x}}{\lambda }+\frac{1}{\lambda }{\cdot \frac{{\text{e}}^{-\lambda t}}{-\lambda }|}_0^x\right]\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{2}{\lambda }\left[\frac{-x{\text{e}}^{-\lambda x}}{\lambda }+\frac{1}{{\lambda }^2}\left(1-{\text{e}}^{-\lambda x}\right)\right]\\ =\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }+\frac{-2x{\text{e}}^{-\lambda x}}{{\lambda }^2}+\frac{1}{{\lambda }^3}\left(1-{\text{e}}^{-\lambda x}\right)\end{array}$ (5)

Putting (4) and (5) into (3) gives

$\begin{array}{c}F\left(x;\lambda \right)=\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)\cdot \frac{1}{\lambda }\left(1-{\text{e}}^{-\lambda x}\right)+\frac{5{\lambda }^2}{2}\left(\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }-\frac{2x{\text{e}}^{-\lambda x}}{{\lambda }^2}+\frac{2}{{\lambda }^3}\left(1-{\text{e}}^{-\lambda x}\right)\right)\right]\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)\cdot \frac{1}{\lambda }\left(1-{\text{e}}^{-\lambda x}\right)+\frac{5{\lambda }^2}{2}\left(\frac{-x^2{\text{e}}^{-\lambda x}}{\lambda }\frac{-2x{\text{e}}^{-\lambda x}}{{\lambda }^2}-\frac{2}{{\lambda }^3}{\text{e}}^{-\lambda x}+\frac{2}{{\lambda }^3}\right)\right]\\ =\frac{1}{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)\left(1-{\text{e}}^{-\lambda x}\right)+\frac{5\lambda }{2}\left(-x^2{\text{e}}^{-\lambda x}-\frac{2x{\text{e}}^{-\lambda x}}{\lambda }-\frac{2}{{\lambda }^2}{\text{e}}^{-\lambda x}+\frac{2}{{\lambda }^2}\right)\right]\\ =\frac{1}{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)\left(1-{\text{e}}^{-\lambda x}\right)-\frac{5\lambda x^2}{2}{\text{e}}^{-\lambda x}-5x{\text{e}}^{-\lambda x}-\frac{5}{\lambda }{\text{e}}^{-\lambda x}+\frac{5}{\lambda }\right]\\ =\frac{1}{{\lambda }^2+6}\left[\left({\lambda }^2+1+\frac{5}{\lambda }\right)-\left({\lambda }^2+1+\frac{5}{\lambda }\right){\text{e}}^{-\lambda x}-\frac{5\lambda x^2}{2}{\text{e}}^{-\lambda x}-5x{\text{e}}^{-\lambda x}\right]\end{array}$

which simplifies to Theorem 2.

2.3. Survival and Hazard Functions of the Nova Distribution

The survival function of $X~Nova\left(\lambda \right)$ is

$S\left(x;\lambda \right)=1-F\left(x;\lambda \right)=\left[1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda x},x>0,\lambda >0$ (6)

The hazard rate function of $X~Nova\left(\lambda \right)$ is

$\begin{array}{c}h\left(x;\lambda \right)=\frac{f\left(x;\lambda \right)}{S\left(x;\lambda \right)}\\ =\frac{\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}}{\left[1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda x}}\\ =\frac{2\lambda \left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]}{2\left({\lambda }^2+6\right)+5{\lambda }^2{x}^2+10\lambda x}\end{array}$ (7)

2.4. Moments of the Proposed Nova Distribution

Theorem 3: The r-th raw moment of $X~Nova\left(\lambda \right)$ is

${{\mu }^{\prime }}_r=E\left(X^r\right)=\frac{r\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(r+1\right)\left(r+1\right)\right]}{{\lambda }^r\left({\lambda }^2+6\right)},r=1,2,3,\cdots$ (8)

Proof: By definition

$\begin{array}{c}{{\mu }^{\prime }}_r=E\left(X^r\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^r\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\text{d}x}\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right){\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^r{\text{e}}^{-\lambda x}}\text{d}x+\frac{5{\lambda }^2}{2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^{r+2}{\text{e}}^{-\lambda x}\text{d}x}\right]\end{array}$ (9)

Using the gamma function

${\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^n{\text{e}}^{-\lambda x}\text{d}x}=\frac{n!}{{\lambda }^{n+1}}$ ,

we obtain

$\begin{array}{c}{{\mu }^{\prime }}_r=E\left(X^r\right)\\ =\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)\frac{r!}{{\lambda }^{r+1}}+\frac{5{\lambda }^2}{2}\frac{\left(r+2\right)!}{{\lambda }^{r+3}}\right]\\ =\frac{r!\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(r+1\right)\left(r+2\right)\right]}{{\lambda }^r\left({\lambda }^2+6\right)}\end{array}$ (10)

This completes the proof.

Corollary 3: The first four raw moments of $X~Nova\left(\lambda \right)$ are:

1) First raw moment (Mean) is:

$\begin{array}{c}{{\mu }^{\prime }}_1=\mu =E\left(X\right)=\frac{1!\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(2\right)\left(3\right)\right]}{{\lambda }^1\left({\lambda }^2+6\right)}\\ =\frac{{\lambda }^2+1+15}{\lambda \left({\lambda }^2+6\right)}=\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\end{array}$ (11)

2) Second raw moment is:

$\begin{array}{c}{{\mu }^{\prime }}_2=E\left(X^2\right)=\frac{2!\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(3\right)\left(4\right)\right]}{{\lambda }^2\left({\lambda }^2+6\right)}\\ =\frac{2\left({\lambda }^2+1+30\right)}{{\lambda }^2\left({\lambda }^2+6\right)}=\frac{2\left({\lambda }^2+31\right)}{{\lambda }^2\left({\lambda }^2+6\right)}\end{array}$ (12)

3) Third raw moment is:

$\begin{array}{c}{{\mu }^{\prime }}_3=E\left(X^3\right)=\frac{3!\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(4\right)\left(5\right)\right]}{{\lambda }^3\left({\lambda }^2+6\right)}\\ =\frac{6\left({\lambda }^2+1+50\right)}{{\lambda }^3\left({\lambda }^2+6\right)}=\frac{6\left({\lambda }^2+51\right)}{{\lambda }^3\left({\lambda }^2+6\right)}\end{array}$ (13)

4) Fourth raw moment is:

$\begin{array}{c}{{\mu }^{\prime }}_4=E\left(X^4\right)=\frac{4!\left[\left({\lambda }^2+1\right)+\frac{5}{2}\left(5\right)\left(6\right)\right]}{{\lambda }^4\left({\lambda }^2+6\right)}\\ =\frac{24\left({\lambda }^2+1+75\right)}{{\lambda }^4\left({\lambda }^2+6\right)}=\frac{24\left({\lambda }^2+76\right)}{{\lambda }^4\left({\lambda }^2+6\right)}\end{array}$ (14)

Corollary 4: The first four central moments of $X~Nova\left(\lambda \right)$ are:

1) First central moment is given by

$\begin{array}{c}{\mu }_1={{\mu }^{\prime }}_1-{{\mu }^{\prime }}_1\\ =\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}-\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}=0\end{array}$ (15)

2) Second raw moment (variance) is given by

$\begin{array}{c}{\mu }_2={{\mu }^{\prime }}_2-{\left(\mu \right)}^2\\ =\frac{2\left({\lambda }^2+31\right)}{{\lambda }^2\left({\lambda }^2+6\right)}-{\left[\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\right]}^2\\ =\frac{2\left({\lambda }^2+31\right)\left({\lambda }^2+6\right)-{\left({\lambda }^2+16\right)}^2}{{\lambda }^2{\left({\lambda }^2+6\right)}^2}\\ =\frac{2\left({\lambda }^4+37{\lambda }^2+186\right)-\left({\lambda }^4+32{\lambda }^2+256\right)}{{\lambda }^2{\left({\lambda }^2+6\right)}^2}\\ =\frac{{\lambda }^4+42{\lambda }^2+116}{{\lambda }^2{\left({\lambda }^2+6\right)}^2}\end{array}$ (16)

3) Third raw moment is given by

$\begin{array}{c}{\mu }_3={{\mu }^{\prime }}_3-3\mu {{\mu }^{\prime }}_2+2{\left(\mu \right)}^3\\ =\frac{6\left({\lambda }^2+51\right)}{{\lambda }^3\left({\lambda }^2+6\right)}-3\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\cdot \frac{2\left({\lambda }^2+31\right)}{{\lambda }^2\left({\lambda }^2+6\right)}+2{\left[\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\right]}^3\\ =\frac{2\left({\lambda }^6+78{\lambda }^4+378{\lambda }^2+676\right)}{{\lambda }^3{\left({\lambda }^2+6\right)}^3}\end{array}$ (17)

4) Fourth raw moment is given by

$\begin{array}{c}{\mu }_4={{\mu }^{\prime }}_4-4\mu {{\mu }^{\prime }}_3+6{\mu }^2{{\mu }^{\prime }}_2-3{\left(\mu \right)}^4\\ =\frac{24\left({\lambda }^2+76\right)}{{\lambda }^4\left({\lambda }^2+6\right)}-4\cdot \frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\cdot \frac{6\left({\lambda }^2+51\right)}{{\lambda }^3\left({\lambda }^2+6\right)}\\ +6{\left(\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\right)}^2\cdot \frac{2\left({\lambda }^2+31\right)}{{\lambda }^2\left({\lambda }^2+6\right)}-3{\left[\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+6\right)}\right]}^4\\ =\frac{9{\lambda }^8+996{\lambda }^6+10584{\lambda }^4+45216{\lambda }^2+63744}{{\lambda }^4{\left({\lambda }^2+6\right)}^4}\end{array}$ (18)

2.5. Skewness and Kurtosis of the Nova Distribution

The skewness of $X~Nova\left(\lambda \right)$ is

$\begin{array}{c}{\gamma }_1=\frac{{\mu }_3}{{\mu }_2^{3/2}}=\frac{2\left({\lambda }^6+78{\lambda }^4+378{\lambda }^2+676\right)/{\lambda }^3{\left({\lambda }^2+6\right)}^3}{{\left[\left({\lambda }^4+42{\lambda }^2+116\right)/{\lambda }^2{\left({\lambda }^2+6\right)}^2\right]}^{3/2}}\\ =\frac{2\left({\lambda }^6+78{\lambda }^4+378{\lambda }^2+676\right)}{{\left({\lambda }^4+42{\lambda }^2+116\right)}^{3/2}}\end{array}$ (19)

The kurtosis of $X~Nova\left(\lambda \right)$ is

$\begin{array}{c}{\gamma }_2=\frac{{\mu }_4}{{\mu }_2^2}=\frac{9{\lambda }^8+996{\lambda }^6+10584{\lambda }^4+45216{\lambda }^2+63744/{\lambda }^4{\left({\lambda }^2+6\right)}^4}{{\left[\left({\lambda }^4+42{\lambda }^2+116\right)/{\lambda }^2{\left({\lambda }^2+6\right)}^2\right]}^2}\\ =\frac{9{\lambda }^8+996{\lambda }^6+10584{\lambda }^4+45216{\lambda }^2+63744}{{\left({\lambda }^4+42{\lambda }^2+116\right)}^{3/2}}\end{array}$ (20)

2.6. Quantile Function of the Nova Distribution

$F\left(Q\left(p;\lambda \right);\lambda \right)=p$

$1-\left[1+\frac{5{\lambda }^2{Q}^2+10\lambda Q}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda Q}=p$

$\left[1+\frac{5{\lambda }^2{Q}^2+10\lambda Q}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda Q}=1-p$ (21)

This is a transcendental equation that cannot be solved analytically in closed form for $Q$ . Hence, we use the Newton-Raphson method to obtain the approximate solution. For a given probability $p$ and parameter $\lambda$ , we solve

$g\left(Q\right)=\left[1+\frac{5{\lambda }^2{Q}^2+10\lambda Q}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda Q}-\left(1-p\right)=0$ (22)

The iteration equation of the Newton-Raphson method for finding an approximate root of the equation is

$Q_{n+1}=Q_n-\frac{g\left(Q\right)}{g^{\prime }\left(Q\right)},n=0,1,2,\cdots$ (23)

where

$g^{\prime }\left(Q\right)=\left[\frac{10{\lambda }^{2}Q+10\lambda }{2\left({\lambda }^2+6\right)}-\lambda \left(1+\frac{5{\lambda }^2{Q}^2+10\lambda Q}{2\left({\lambda }^2+6\right)}\right)\right]{\text{e}}^{-\lambda Q}$ (24)

2.7. Median of the Proposed Distribution

The median of $X~Nova\left(\lambda \right)$ satisfies

$\left[1+\frac{5{\lambda }^2{Q}_{0.5}^2+10\lambda Q_{0.5}}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda Q_{0.5}}=0.5$ (25)

2.8. Order Statistics for the Nova Distribution

Let $X_1,X_2,\cdots ,X_n$ be a random sample from $X~Nova\left(\lambda \right)$ , and let $X_{\left(1\right)}\le X_{\left(2\right)}\le \cdots \le X_{\left(n\right)}$ be the order statistics. The PDF of the kth order statistic is:

$\begin{array}{c}{f}_{X_{\left(k\right)}}\left(x\right)=\frac{n!}{\left(k-1\right)!\left(n-k\right)!}{\left[F\left(x\right)\right]}^{k-1}{\left[1-F\left(x\right)\right]}^{n-k}f\left(x\right)\\ =\frac{n!}{\left(k-1\right)!\left(n-k\right)!}{\left[1-\left(1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right){\text{e}}^{-\lambda x}\right]}^{k-1}\\ \times {\left[\left(1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right){\text{e}}^{-\lambda x}\right]}^{n-k}\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\\ =\frac{n!\lambda }{\left(k-1\right)!\left(n-k\right)!\left({\lambda }^2+6\right)}{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{k-1}\\ \times {\left[A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-k}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x\left(n-k+1\right)}\end{array}$ (26)

where

$A\left(x\right)=\left(1+\frac{5{\lambda }^2{x}^2+10\lambda x}{2\left({\lambda }^2+6\right)}\right)$ (27)

In particular, PDF of the minimum $X_{\left(1\right)}$ and PDF of the maximum $X_{\left(2\right)}$

$\begin{array}{c}{f}_{X_{\left(1\right)}}\left(x\right)=n{\left[1-F\left(x\right)\right]}^{n-1}f\left(x\right)\\ =n{\left[A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-1}\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\\ =\frac{n\lambda }{{\lambda }^2+6}{\left[A\left(x\right)\right]}^{n-1}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-n\lambda x}\end{array}$ (28)

and

$\begin{array}{c}{f}_{X_{\left(n\right)}}\left(x\right)=n{\left[F\left(x\right)\right]}^{n-1}f\left(x\right)\\ =n{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-1}\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\\ =\frac{n\lambda }{{\lambda }^2+6}{\left[1-A\left(x\right)\right]}^{n-1}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\end{array}$ (29)

The corresponding CDF for the minimum $X_{\left(1\right)}$ and CDF maximum $X_{\left(2\right)}$ are:

$F_{X_{\left(1\right)}}\left(x\right)=1-{\left[1-F\left(x\right)\right]}^n=1-{\left[A\left(x\right)\right]}^n{\text{e}}^{-n\lambda x}$ (30)

and

$F_{X_{\left(n\right)}}\left(x\right)=n\left[F\left(x\right)\right]{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^n$ (31)

The joint distribution of the order statistics is

$\begin{array}{l}{f}_{X_{\left(i\right)},X_{\left(j\right)}}\left(x,y\right)\\ =\{\begin{array}{l}\frac{n!}{\left(i-1\right)!\left(j-i-1\right)!\left(n-j\right)!}{\left[F\left(x\right)\right]}^{i-1}{\left[F\left(y\right)-F\left(x\right)\right]}^{j-i-1}\\ \times {\left[1-F\left(y\right)\right]}^{n-j}f\left(x\right)f\left(y\right),\text{for }x<y\\ \frac{n!}{\left(i-1\right)!\left(j-i-1\right)!\left(n-j\right)!}{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{i-1}{\left[A\left(y\right){\text{e}}^{-\lambda y}-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{j-i-1}\\ \times {\left[1-A\left(x\right){\text{e}}^{-\lambda y}\right]}^{n-j}\times {\left(\frac{\lambda }{{\lambda }^2+6}\right)}^2\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]\\ \times \left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{y}^2\right]{\text{e}}^{-\lambda \left(x+y\right)},\text{for }x<y\end{array}\end{array}$ (32)

The range $R=X_{\left(n\right)}-X_{\left(1\right)}$ has the PDF

$\begin{array}{c}{f}_R\left(r\right)=n\left(n-1\right){\displaystyle \underset{0}{\overset{\infty }{\int }}{\left[F\left(x+r\right)-F\left(x\right)\right]}^{n-2}f\left(x\right)f\left(x+r\right)\text{d}x}\\ =\frac{n\left(n-1\right){\lambda }^2}{{\left({\lambda }^2+6\right)}^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left[A\left(x\right){\text{e}}^{-\lambda x}-A\left(x+r\right){\text{e}}^{-\lambda \left(x+r\right)}\right]}^{n-2}}\\ \times \left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{\left(x+r\right)}^2\right]{\text{e}}^{-\lambda \left(2x+r\right)}\text{d}x\end{array}$ (33)

The r-th moment of the kth order statistic is

$\begin{array}{c}E\left[X_{\left(k\right)}^r\right]=\frac{n!}{\left(k-1\right)!\left(n-k\right)!}{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^r{\left[F\left(x\right)\right]}^{k-1}{\left[1-F\left(x\right)\right]}^{n-k}f\left(x\right)\text{d}x}\\ =\frac{n!\lambda }{\left(k-1\right)!\left(n-k\right)!\left({\lambda }^2+6\right)}\\ \times {\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^r{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{k-1}{\left[A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-k}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x\left(n-k+1\right)}\text{d}x}\end{array}$ (34)

In particular, the mean of first-order statistic $X_{\left(1\right)}$ is

$E\left[X_{\left(1\right)}\right]=\frac{n\lambda }{{\lambda }^2+6}{\displaystyle \underset{0}{\overset{\infty }{\int }}x{\left[A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-1}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-n\lambda x}\text{d}x}$ (35)

Also, the mean of n-th order statistic $X_{\left(n\right)}$ is

$E\left[X_{\left(n\right)}\right]=\frac{n\lambda }{{\lambda }^2+6}{\displaystyle \underset{0}{\overset{\infty }{\int }}x{\left[1-A\left(x\right){\text{e}}^{-\lambda x}\right]}^{n-1}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-n\lambda x}\text{d}x}$ (36)

2.9. Mean Deviation of the Nova Distribution

The mean deviation of the Nova distribution is

$\begin{array}{l}MD=E\left[\left|X-\mu \right|\right]={\displaystyle \underset{0}{\overset{\mu }{\int }}\left|x-\mu \right|f\left(x\right)\text{d}x}\\ ={\displaystyle \underset{0}{\overset{\mu }{\int }}\left(\mu -x\right)f\left(x\right)\text{d}x}+{\displaystyle \underset{\mu }{\overset{\infty }{\int }}\left(x-\mu \right)f\left(x\right)\text{d}x}=2\mu F\left(\mu \right)-2{\displaystyle \underset{0}{\overset{\mu }{\int }}xf\left(x\right)\text{d}x}\\ =2\mu \left[1-\left(1+\frac{5{\lambda }^2{\mu }^2+10\lambda \mu }{2\left({\lambda }^2+6\right)}\right){\text{e}}^{-\lambda \mu }\right]-2{\displaystyle \underset{0}{\overset{\mu }{\int }}x\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\text{d}x}\\ =2\mu \left[1-\left(1+\frac{5{\lambda }^2{\mu }^2+10\lambda \mu }{2\left({\lambda }^2+6\right)}\right){\text{e}}^{-\lambda \mu }\right]-\frac{2\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right){\displaystyle \underset{0}{\overset{\mu }{\int }}x{\text{e}}^{-\lambda x}\text{d}x}+\frac{5{\lambda }^2}{2}{\displaystyle \underset{0}{\overset{\mu }{\int }}{x}^2{\text{e}}^{-\lambda x}\text{d}x}\right]\\ =2\mu \left[1-A\left(\mu \right){\text{e}}^{-\lambda x}\right]-\frac{2\lambda }{{\lambda }^2+6}I\end{array}$ (37)

where

$I=\frac{{\lambda }^2+1}{{\lambda }^2}\left[1-\left(1+\lambda \mu \right){\text{e}}^{-\lambda \mu }\right]+\frac{5{\lambda }^2}{2{\lambda }^4}\left[6-6\left(6+6\lambda \mu +3{\lambda }^2{\mu }^2+{\lambda }^3{\mu }^3\right){\text{e}}^{-\lambda \mu }\right]$

$A\left(\mu \right)=1+\frac{5{\lambda }^2{\mu }^2+10\lambda \mu }{2\left({\lambda }^2+6\right)}$

and

$\mu =\frac{{\lambda }^2+16}{\lambda \left({\lambda }^2+1\right)}$

2.10. Bonferroni and Lorenz Curves for the Nova Distribution

The Bonferroni curve measures inequality and is defined for a given probability level $p$ as

$\begin{array}{c}B\left(p\right)=\frac{1}{p\mu }{\displaystyle \underset{0}{\overset{x_p}{\int }}xf\left(x\right)\text{d}x}\\ =\frac{1}{p\mu }{\displaystyle \underset{0}{\overset{x_p}{\int }}x\frac{\lambda }{{\lambda }^2+6}\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}^2\right]{\text{e}}^{-\lambda x}\text{d}x}\\ =\frac{1}{p\mu }\frac{\lambda }{{\lambda }^2+6}I\left(x_p\right)\end{array}$ (38)

where

$I\left(x_p\right)=\frac{{\lambda }^2+1}{{\lambda }^2}\left[1-\left(1+\lambda x_p\right){\text{e}}^{-\lambda x_p}\right]+\frac{5{\lambda }^2}{2{\lambda }^4}\left[6-6\left(6+6\lambda x_p+3{\lambda }^2{x}_p^2+{\lambda }^3{x}_p^3\right){\text{e}}^{-\lambda x_p}\right]$

and the p-th quantile $x_p$ is obtained by solving:

$\left[1+\frac{5{\lambda }^2{x}_p^2+10\lambda x_p}{2\left({\lambda }^2+6\right)}\right]{\text{e}}^{-\lambda x_p}=1-p$ (39)

Similarly, the Lorenz curve measures the cumulative share of total income or wealth held by the bottom proportion of the population and is defined as

$B\left(p\right)=\frac{1}{\mu }{\displaystyle \underset{0}{\overset{x_p}{\int }}xf\left(x\right)\text{d}x}=\frac{\lambda }{\mu \left({\lambda }^2+6\right)}I\left(x_p\right)$ (40)

2.11. Inequality Measures for the Nova Distribution

The Gini coefficient and Bonferroni index are the two measures of inequality derived for the Nova distribution, and they are given respectively as

$G=1-2{\displaystyle \underset{0}{\overset{1}{\int }}L\left(p\right)\text{d}p}$ (41)

and

$B=1-{\displaystyle \underset{0}{\overset{1}{\int }}B\left(p\right)\text{d}p}$ (42)

2.12. Stress-Strength Reliability for the Nova Distribution

The Stress-Strength Reliability is given by

$\begin{array}{c}R=P\left(X>Y\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \underset{y}{\overset{\infty }{\int }}{f}_X\left(x\right)f_Y\left(y\right)\text{d}x\text{d}y}}={\displaystyle \underset{0}{\overset{\infty }{\int }}{F}_Y\left(x\right)f_x\left(x\right)\text{d}y}\\ ={\displaystyle \underset{0}{\overset{\infty }{\int }}\left[1-\left(1+\frac{5{\lambda }_2^2{x}^2+10{\lambda }_{2}x}{2\left({\lambda }_2^2+6\right)}{\text{e}}^{-{\lambda }_{2}x}\right)\right]\frac{{\lambda }_1}{{\lambda }_1^2+6}\left[\left({\lambda }_1^2+1\right)+\frac{5{\lambda }_1^2}{2}{x}^2\right]{\text{e}}^{-{\lambda }_{1}x}\text{d}x}\\ =\frac{{\lambda }_1}{{\lambda }_1^2+6}{\displaystyle \underset{0}{\overset{\infty }{\int }}\left[\left({\lambda }_1^2+1\right)+\frac{5{\lambda }_1^2}{2}{x}^2\right]{\text{e}}^{-{\lambda }_{1}x}\text{d}x}\\ -\frac{{\lambda }_1}{{\lambda }_1^2+6}{\displaystyle \underset{0}{\overset{\infty }{\int }}\left(1+\frac{5{\lambda }_2^2{x}^2+10{\lambda }_{2}x}{2\left({\lambda }_2^2+6\right)}\right)\left[\left({\lambda }_1^2+1\right)+\frac{5{\lambda }_1^2}{2}{x}^2\right]{\text{e}}^{-\left({\lambda }_1+{\lambda }_2\right)x}\text{d}x}\\ =1-\frac{{\lambda }_1}{{\lambda }_1^2+6}{\displaystyle \underset{0}{\overset{\infty }{\int }}\left(1+\frac{5{\lambda }_2^2{x}^2+10{\lambda }_{2}x}{2\left({\lambda }_2^2+6\right)}\right)\left[\left({\lambda }_1^2+1\right)+\frac{5{\lambda }_1^2}{2}{x}^2\right]{\text{e}}^{-\left({\lambda }_1+{\lambda }_2\right)x}\text{d}x}\\ =1-\left[\left({\lambda }_1^2+1\right)+\frac{5{\lambda }_1^2}{2}{x}^2+\frac{5{\lambda }_2^2\left({\lambda }_1^2+1\right)x^2}{2\left({\lambda }_2^2+6\right)}+\frac{25{\lambda }_1^2{\lambda }_2^2}{4\left({\lambda }_2^2+6\right)}{x}^4+\frac{10{\lambda }_2\left({\lambda }_1^2+1\right)}{2\left({\lambda }_2^2+6\right)}x+\frac{25{\lambda }_1^2{\lambda }_2}{4\left({\lambda }_2^2+6\right)}{x}^3\right]\\ =1-\frac{{\lambda }_2}{{\lambda }_1^2+6}\left[\frac{{\lambda }_1^2+1}{\tau }+\frac{5{\lambda }_1^2}{{\tau }^3}+\frac{5{\lambda }_2^2\left({\lambda }_1^2+1\right)}{2\left({\lambda }_2^2+6\right){\tau }^3}+\frac{25{\lambda }_1^2{\lambda }_2^2}{2\left({\lambda }_2^2+6\right){\tau }^5}+\frac{5{\lambda }_2\left({\lambda }_1^2+1\right)}{\left({\lambda }_2^2+6\right){\tau }^2}+\frac{15{\lambda }_1^2{\lambda }_2}{\left({\lambda }_2^2+6\right){\tau }^4}\right]\end{array}$ (43)

where $\tau ={\lambda }_1+{\lambda }_2$ .

2.13. Maximum Likelihood Estimation of Parameter of the Nova Distribution

Let $X_1,X_2,\cdots ,X_n$ be a random sample from $Nova\left(\lambda \right)$ . The likelihood function is:

$L\left(\lambda \right)={\left(\frac{\lambda }{{\lambda }^2+6}\right)}^n{\displaystyle \prod _{i=1}^n\left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}_i^2\right]{\text{e}}^{-\lambda x_i}}$ (44)

The log-likelihood function is:

$\ln L\left(\lambda \right)=n\ln \lambda -n\ln \left({\lambda }^2+6\right)+{\displaystyle \sum _{i=1}^n\ln \left[\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}_i^2\right]}-\lambda {\displaystyle \sum _{i=1}^n{x}_i}$ (45)

The maximum likelihood estimator is obtained by solving numerically the equation:

$\frac{\text{d}}{\text{d}\lambda }\ln L\left(\lambda \right)=0$

$\frac{n}{\lambda }-\frac{2n\lambda }{{\lambda }^2+6}+{\displaystyle \sum _{i=1}^n\frac{2\lambda +5\lambda x_i^2}{\left({\lambda }^2+1\right)+\frac{5{\lambda }^2}{2}{x}_i^2}}-{\displaystyle \sum _{i=1}^n{x}_i}=0$ (46)

3. Application of the Nova Distribution to G/G/1 Queueing Theory

3.1. Formulation of the G/G/1 Queueing Model Based on the Nova Distribution

To formulate the single-server, single-queue system, denoted as G_N/G_N/1, we assume that:

1) The system consists of a single service channel and a single waiting line, i.e., a single-server, single-queue system.

2) Customers are served on a first-come, first-served (FCFS) basis.

3) The time between consecutive customer arrivals is modeled using the Nova distribution with parameter ${\lambda }_a$ , which captures the variability in arrival patterns. That is,

$f_A\left(t\right)=\frac{{\lambda }_a}{{\lambda }_a^2+6}\left[\left({\lambda }_a^2+1\right)+\frac{5{\lambda }_a^2}{2}{t}^2\right]{\text{e}}^{-{\lambda }_{a}t}$ (47)

$E\left(A\right)=\frac{{\lambda }_a^2+16}{{\lambda }_a\left({\lambda }_a^2+6\right)}$ (48)

${\eta }_a=\frac{1}{E\left(A\right)}=\frac{1}{\frac{{\lambda }_a^2+16}{{\lambda }_a\left({\lambda }_a^2+6\right)}}=\frac{{\lambda }_a\left({\lambda }_a^2+6\right)}{{\lambda }_a^2+16}$ (49)

$\begin{array}{c}{C}_a^2=\frac{Var\left(A\right)}{{\left[E\left(A\right)\right]}^2}=\frac{\left({\lambda }_a^4+42{\lambda }_a^2+116\right)/{\lambda }_a^2\left({\lambda }_a^2+6\right)}{{\left[\left({\lambda }_a^2+16\right)/{\lambda }_a\left({\lambda }_a^2+6\right)\right]}^2}\\ =\frac{\left({\lambda }_a^4+42{\lambda }_a^2+116\right){\left({\lambda }_a^2+6\right)}^2}{{\lambda }_a^2+16}\end{array}$ (50)

4) The time required to serve each customer is also modeled using a Nova distribution with parameter ${\mu }_s$ , representing the stochastic nature of service durations. Thus,

$f_S\left(t\right)=\frac{{\mu }_s}{{\mu }_s^2+6}\left[\left({\mu }_s^2+1\right)+\frac{5{\mu }_s^2}{2}{t}^2\right]{\text{e}}^{-{\mu }_{s}t}$ (51)

$E\left(S\right)=\frac{{\mu }_s^2+16}{{\mu }_s\left({\mu }_s^2+6\right)}$ (52)

${\eta }_s=\frac{1}{E\left(S\right)}=\frac{1}{\frac{{\mu }_s^2+16}{{\mu }_s\left({\mu }_s^2+6\right)}}=\frac{{\mu }_s\left({\mu }_s^2+6\right)}{{\mu }_s^2+16}$ (53)

$\begin{array}{c}{C}_s^2=\frac{Var\left(S\right)}{{\left[E\left(S\right)\right]}^2}\\ =\frac{\left({\mu }_s^4+42{\mu }_s^2+116\right)/{\mu }_s^2\left({\mu }_s^2+6\right)}{{\left[\left({\mu }_s^2+16\right)/{\mu }_s\left({\mu }_s^2+6\right)\right]}^2}\\ =\frac{\left({\mu }_s^4+42{\mu }_s^2+116\right){\left({\mu }_s^2+6\right)}^2}{{\mu }_s^2+16}\end{array}$ (54)

5) The system is assumed to be stable, meaning that the traffic intensity satisfies $\rho <1$ , ensuring that the queue does not grow indefinitely over time.

3.2. Derivation of Performance Measures of the Lindley-Based G/G/1 Queueing Model

Then, the traffic intensity ρ is defined as the ratio of

$\rho =\lambda E\left(S\right)=\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\mu }_a^2+16\right)}{{\mu }_s\left({\mu }_a^2+6\right)\left({\lambda }_a^2+16\right)}$ (55)

3.2.1. Average Waiting Time in Queue ($W_q$ )

The average waiting time in the queue is derived by applying the Pollaczek-Khinchine formula for the general G/G/1 case, which incorporates the variance of both inter-arrival and service times:

$W_q=\frac{\lambda \left(E\left(S^2\right)+E\left(A^2\right)-2E\left(A\right)E\left(S\right)\right)}{2\left(1-\rho \right)}$ (56)

However, a more direct and commonly used form of the G/G/1 approximation, based on the Kingman’s formula or Kiefer-Wolfowitz formula, is more suitable for this context s more suitable for this context because the Pollaczek-Khinchine formula is exact only for M/G/1 queues (requiring Poisson arrivals), whereas Kingman’s formula provides a robust heavy-traffic approximation for the general G/G/1 case with arbitrary arrival and service distributions:

$W_q=\frac{\rho E\left(S^2\right)+{\eta }_{a}Var\left(A\right)}{2\left(1-\rho \right)}$ (57)

Substituting the first and second raw moments of the Nova distribution for inter-arrival and service times, we derive the expression for the average waiting time in the queue as:

$\begin{array}{c}{W}_q=\frac{\rho \left[\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right]+\left(\frac{{\lambda }_a\left({\lambda }_a^2+6\right)}{{\lambda }_a^2+16}\right)\left(\frac{{\lambda }_a^4+42{\lambda }_a^2+116}{{\lambda }_a^2\left({\lambda }_a^2+6\right)}\right)}{2\left(1-\rho \right)}\\ =\frac{\rho \left[\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right]+\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\lambda }_a^4+42{\lambda }_a^2+116\right)}{{\lambda }_a^2\left({\lambda }_a^2+16\right)\left({\lambda }_a^2+6\right)}}{2\left(1-\rho \right)}\end{array}$ (58)

3.2.2. Average Time in System ($W_s$ )

By Little’s Law, the average time a patient spends in the system is the sum of the average waiting time and the average service time.

$W_s=W_q+E\left(S\right)$ (59)

Substituting the derived formula for $W_q$ and the mean of the proposed Nova service time:

$W_s=\frac{\rho \left[\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right]+\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\lambda }_a^4+42{\lambda }_a^2+116\right)}{{\lambda }_a^2\left({\lambda }_a^2+16\right)\left({\lambda }_a^2+6\right)}}{2\left(1-\rho \right)}+\frac{{\mu }_s^2+16}{{\mu }_s\left({\mu }_s^2+6\right)}$ (60)

3.2.3. Average Number in Queue ($L_q$ )

The average number of patients in the queue is given by Little’s Law:

$L_q={\eta }_a{W}_q$

$L_q=\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left[\rho \left(\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right)+\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\lambda }_a^4+42{\lambda }_a^2+116\right)}{{\lambda }_a^2\left({\lambda }_a^2+16\right)\left({\lambda }_a^2+6\right)}\right]}{2\left(1-\rho \right)\left({\lambda }_a^2+16\right)}$ (61)

3.2.4. Average Number in System ($L_s$ )

The average number of patients in the system is also given by Little’s Law:

$L_s={\eta }_a{W}_s$

$L_s={\eta }_a\cdot \left(\frac{\rho \left(\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right)+\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\lambda }_a^4+42{\lambda }_a^2+116\right)}{{\lambda }_a^2\left({\lambda }_a^2+16\right)\left({\lambda }_a^2+6\right)}}{2\left(1-\rho \right)}+\frac{{\mu }_s^2+16}{{\mu }_s\left({\mu }_s^2+6\right)}\right)$ (62)

3.3. G/G/1 Lindley Queueing Cost Model

To derive the cost function for the G/G/1 Nova distributed interarrival and service time distribution, we observe that the total cost function is a sum of the cost of waiting (per unit of time) and cost of service (per unit of time). Notably, the cost of waiting is given by

$T_W=C_w{W}_q{\eta }_a$ (63)

where

$W_q$ is the average waiting time in the queue per customer;

$C_w$ is the cost incurred per unit of customer waiting time;

${\eta }_a$ is the number of customers arriving per unit of time.

Also, the cost of service per unit of time is

$T_S=C_s$ (64)

where

$C_s$ is the fixed cost per unit of time for having a single server.

Consequently, the total cost is the sum of the cost of waiting and the cost of service, both measured per unit of time, given by

$T_C=C_w{W}_q{\eta }_a+C_s$ (65)

Substituting the parameters gives

$T_C=C_w\frac{\rho \left[\frac{2\left({\mu }_s^2+31\right)}{{\mu }_s^2\left({\mu }_s^2+6\right)}\right]+\frac{{\lambda }_a\left({\lambda }_a^2+6\right)\left({\lambda }_a^4+42{\lambda }_a^2+116\right)}{{\lambda }_a^2\left({\lambda }_a^2+16\right)\left({\lambda }_a^2+6\right)}}{2\left(1-\rho \right)}{\eta }_a+C_s$ (66)

3.4. Data Collection

The data for this study were collected from First Bank Plc, Douglas Road, Owerri Branch, Imo State, Nigeria, specifically from the teller section where customers perform over-the-counter transactions such as deposits, withdrawals, and account inquiries. Although commercial banks generally operate as multi-server queueing systems, with multiple tellers serving customers simultaneously. The choice of a single-channel queue was deliberate and analytically justified for three reasons. First, it provides a controlled abstraction to validate the Nova distribution’s efficacy for modeling service times without the confounding complexity of multi-server interactions. Second, in the observed environment, tellers often operate with dedicated queues; therefore, analyzing one queue in isolation offers direct insights into per-server performance metrics like individual utilization. Finally, establishing a robust single-server model is a necessary prerequisite for future extensions to multi-server (G/G/c) systems. This approach reflects realistic queueing behavior, as customers often form individual lines before their preferred tellers, resulting in distinct, nearly independent service queues that align with single-server assumptions.

However, this single-server simplification introduces limitations regarding the generalizability of absolute performance metrics to the entire branch. A G/G/1 model cannot capture key multi-server dynamics, such as the reduced average wait times from customer pooling in a single, shared queue or the effects of load balancing between tellers. Consequently, while the conclusions on the statistical fit of the Nova distribution remain valid, the precise waiting time or cost estimates are most accurate for the per-server context modeled. This acknowledged limitation underscores the value of extending this framework to a G/G/c model in future research.

The study utilized primary data obtained through direct observation of teller-customer service operations. The data were quantitative, time-based observations of queueing events and service durations, measured in minutes. The variables recorded include customer identification number (assigned sequentially to each observed customer), arrival time, interarrival time, service start time, service time, service end time, waiting time, queue length at arrival, system time, waiting cost (₦), service cost (₦), and total cost (₦). These variables were selected because they represent the fundamental parameters required for deriving queue performance measures, including average waiting time, mean queue length, server utilization rate, and total cost of service operation.

Data collection was carried out using a digital stopwatch, transaction timestamps, and manual recording sheets. The stopwatch was used to measure interarrival and service times with precision, while teller system timestamps were cross-checked to validate accuracy. A structured observation sheet was designed to record all relevant time points for each customer. Two trained research assistants worked alongside the researcher to ensure continuous observation and minimize human error. The raw data were later transferred into Microsoft Excel for organization, cleaning, and subsequent statistical analysis. The use of direct observation was justified by the need for accuracy and objectivity in measuring time-dependent events in the queueing system. Unlike secondary records, which may not provide precise timing information, direct

4. Numerical Results and Discussion

4.1. Goodness-of-Fit Analysis

Table 1. Goodness-of-fit for interarrival times.

Distribution	Estimated Parameter	Log-Lik	AIC	BIC
Nova	3.5241	−85.234	172.469	174.876
Shanker	3.6123	−86.789	175.578	177.986
Lindley	3.8345	−89.901	181.802	184.210

Table 2. Goodness-of-fit for service times.

Distribution	Estimated Parameter	Log Likelihood	AIC	BIC
Nova	1.2189	−192.346	386.691	389.099
Shanker	1.2456	−194.789	391.578	393.986
Lindley	1.2890	−201.901	405.802	408.210

The Nova distribution consistently achieved the lowest AIC and BIC values, indicating a statistically superior fit.

4.2. Queueing Performance and Cost Analysis

Table 3. Queueing performance measures.

Distribution	$\rho$	$W_q$	$L_q$	$W$
Nova	0.9345	3.2345	2.8765	4.1234
Lindley	0.9456	4.1234	3.4567	5.2345
Shanker	0.9289	3.4567	3.0123	4.5678
Exponential	0.9123	2.8765	2.3456	3.7890

Table 4. Cost analysis.

Distribution	Wait Cost	Service Cost	Total Cost	Cost per Customer
Nova	N32345.00	N23456.78	N55801.78	N 55.80
Lindley	N41234.00	N24567.89	N65801.89	N 65.80
Shanker	N34567.00	N23678.90	N58245.90	N 58.25
Exponential	N28765.00	N22345.67	N51110.67	N 51.11

4.3. Discussion of Results

The results of this study demonstrate the practical and statistical value of the proposed Nova distribution within the G/G/1 queueing framework. The discussion is structured around three key insights: the distribution’s superior modeling capability, its impact on queue performance prediction, and the resulting economic implications for system management.

The goodness-of-fit analysis (Table 1, Table 2) provides unequivocal evidence that the Nova distribution outperforms the one-parameter Lindley and Shanker distributions for both interarrival and service time data. By achieving the lowest Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values, the Nova distribution is statistically identified as the best model among the candidates, balancing fit and parsimony. This superiority can be attributed to its genesis as a mixture of exponential and gamma components, which grants it an inherent flexibility to capture a wider range of distributional shapes—particularly the right-skewness and moderate tails commonly observed in real-world service processes. The Lindley and Shanker distributions, while useful, lack this specific mixture structure and are thus less adept at modeling the empirical nuances present in the banking data.

The derived performance measures (Table 3) reveal a critical finding: the Nova-based model predicts system congestion that is more severe yet more accurate than the estimates from fitted Lindley or Shanker models. For instance, the average number in the queue under the Nova model is 3.2345, which lies between the underestimated value from the Exponential fit (2.8765) and the overestimated values from the Lindley model (4.1234). This “balanced” prediction stems directly from the Nova distribution’s accurate capture of the second moments (variance) of the empirical data. Kingman’s approximation (Equation (57)), which incorporates these variances, shows that even slight misestimation of variability can significantly distort waiting time and queue length calculations. The Nova distribution, by providing a superior fit, inputs more realistic moments into the queueing formulas, yielding performance metrics that better reflect the true system behavior under high utilization (ρ ≈ 0.93).

The cost analysis (Table 4) translates the statistical and performance advantages into tangible economic terms. While a naive model like the exponential yields the lowest total cost (₦51110.67), this is an artifact of its underestimation of delays and its failure to fit the data. Among the models that actually fit the data (Nova, Lindley, Shanker), the Nova-based model results in the lowest total expected cost (₦55801.78). This is because it optimally balances the cost of customer waiting against the cost of service provision. The Lindley model, by overpredicting congestion, suggests excessively high waiting costs, which could lead managers to over-invest in capacity. The Nova model provides a more reliable baseline for decision-making.

For a bank branch manager, adopting the Nova-based model means using a tool that more faithfully represents observed customer flow. The accurate prediction of average queue length and waiting time enables precise staffing schedules, targeted process improvements, and reliable customer wait-time announcements. In capacity planning, the model’s accurate cost function allows for the identification of a service rate that minimizes total operational cost without compromising service quality, thereby enhancing both efficiency and customer satisfaction.

The choice to apply Kingman’s heavy-traffic approximation (over the exact Pollaczek-Khinchine formula) is justified by the non-Poisson, Nova-distributed arrival process. Kingman’s formula is a robust approximation for general G/G/1 queues, especially under high traffic intensity, making it the appropriate analytical tool for this study.

It is important to acknowledge the limitations that contextualize these findings. The data originate from a single banking facility, and the model assumes steady-state conditions and homogeneous servers. The single-server (G/G/1) formulation, while a valid simplification for analyzing per-teller performance as noted in Section 3.4, does not capture the full dynamics of a multi-teller queue with potential server cooperation or heterogeneous service rates. These factors present clear avenues for extending the Nova-based framework to G/G/c and more complex queueing networks in future work.

The integration of the Nova distribution into a G/G/1 queueing model creates a powerful, parsimonious, and cost-effective analytical tool. It advances queueing theory by offering a statistically superior one-parameter alternative for modeling service systems and provides practitioners with a more reliable foundation for operational optimization.

5. Conclusions and Recommendations

5.1. Conclusions

This study successfully introduced, characterized, and applied the novel one-parameter Nova distribution within the framework of queueing theory. Derived as a mixture of exponential and gamma components, the Nova distribution was shown to possess tractable statistical properties, making it a practical and flexible tool for stochastic modeling.

The primary contribution of this work is the empirical and analytical demonstration of the Nova distribution’s superiority. Applied to real-world data from a banking service facility, it provided a statistically superior fit for both interarrival and service times compared to established one-parameter alternatives like the Lindley and Shanker distributions, as rigorously confirmed by Akaike and Bayesian Information Criteria.

By integrating the Nova distribution into a G/G/1 queueing model, we derived key performance measures—such as average waiting time and queue length—and an associated economic cost model. The results confirm that the Nova-based model yields more realistic and balanced predictions of system congestion, particularly under the high-utilization conditions prevalent in many service environments. Crucially, this enhanced accuracy translates directly into economic value: among the well-fitting models, the Nova-based framework produced the most cost-effective operational profile, minimizing the total expected cost of service provision and customer waiting.

Therefore, this research establishes the Nova distribution not only as a valuable addition to the family of lifetime distributions but also as the foundation for a more accurate and reliable G/G/1 queueing model for analyzing, optimizing, and managing real-world service systems where data exhibit non-exponential patterns.

5.2. Recommendations

Based on the findings and conclusions of this study, the following actionable recommendations are offered for practitioners, researchers, and software developers:

1) Service industries—including banking, healthcare, and telecommunications—should adopt the Nova-based G/G/1 queueing model for capacity planning, staffing, and congestion management. The model’s proven accuracy in predicting realistic queue lengths and its associated cost-minimizing properties (Section 4.2) provide a reliable, data-driven tool for optimizing service levels and operational efficiency, leading to improved customer satisfaction and resource utilization.

2) Future work should focus on extending the Nova-based framework to multi-server (G/G/c) queueing systems. This is a logical and necessary progression, given that the single-server model used here serves as a foundational simplification of multi-channel service environments. Developing approximations or simulation protocols for the G/G/c-Nova queue would significantly broaden the model’s practical applicability. Further theoretical extensions to incorporate customer impatience (balking and reneging), server vacations, retrial behavior, or batch arrivals would further test and establish the Nova distribution’s versatility in advanced queueing theory.

3) To facilitate widespread adoption and validation, researchers and developers should create and disseminate open-source software libraries implementing the Nova distribution. Integrating its probability density, cumulative distribution, quantile functions, and random variate generation into popular statistical and simulation platforms (e.g., R, Python, MATLAB, Arena, AnyLogic) will lower the barrier to entry for both academic research and industrial application.

4) The promising results from this single-branch case study warrant broader validation across diverse service sectors and geographical contexts. Researchers are encouraged to apply the Nova queueing model to data from hospitals, call centers, retail, and transportation hubs to assess its generalizability and robustness under different operational conditions.