EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 4, Article Number 6653
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

The Mgamma Distribution: Statistical Properties and
Application to Real-Life Data

M. I. Khan1, Molay Kumar Ruidas2,∗

1 Department of Mathematics, Islamic University of Madinah, Madinah,
Kingdom of Saudi Arabia, 42351
2 Department of Statistics, Trivenidevi Bhalotia College, Kazi Nazrul University, Raniganj,
West Bengal, India, 713347

Abstract. We propose a new one-parameter continuous distribution, termed the Mgamma dis-
tribution, obtained as a finite mixture of two Gamma components having different shapes. Key
statistical properties are derived, including moments, generating functions, stochastic ordering, and
reliability measures such as hazard rate and mean residual life. Parameter estimation is addressed
using the method of moments and maximum likelihood, with simulation studies confirming the
consistency and efficiency of estimators. The applicability of the Mgamma distribution is demon-
strated on four real datasets from engineering, hydrology, medicine, and queueing systems. Model
selection criteria (AIC, BIC, CAIC) and goodness-of-fit tests (KS, Cramér–von Mises, Anderson–
Darling) consistently show that Mgamma provides a superior fit compared to classical models such
as Exponential, Lindley, Xgamma, Shanker, and Akash distributions. These findings establish
Mgamma as a flexible and robust alternative for modeling skewed and heavy-tailed lifetime data.

2020 Mathematics Subject Classifications: 60E05, 62E99

Key Words and Phrases: Finite mixture, Gamma distribution, maximum likelihood estimation,
survival properties, stochastic ordering

1. Introduction

In recent years, research has focused on developing continuous distributions that are
more flexible in fitting real lifetime data sets. Some distributions that are finite mixtures
of Gamma distributions introduced by Geoffrey et al. (1988) [1], mixture of Exponential
and Gamma distributions introduced by Ghitany et al. (2008) [2] and Maiti et al. (2024)
[3] used to model lifetime data. Recently, Sen et al. (2016) [4] introduced Xgamma, the
Lindley distribution by DV Lindley (1958) [5], Shanker (2015) [6] introduced the Shanker
distribution, Gamma-Lindley (Zeghdoudi and Nedjar, (2016) [7]) and several other mix-
ture distributions are derived using the literature of Johnson, Kotz and Balakrishnan

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6653

Email addresses: khanizhar@iu.edu.sa (M. I. Khan),
molayruidas@gmail.com (M. K. Ruidas)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 2 of 22

(1994) [8] to fit real-life data. However, these distributions are not as sophisticated in
modeling datasets more accurately. Various classical distributions are modified and ex-
tended to find new distributions like extended Gamma introduced by Bakouch (2012) [9],
two-parameter Xgamma Weibull by Yousof et al. (2021) [10], Exponentiated survival func-
tion of Beta introduced by Mahmoud et al. (2024)[11], new beta power flexible Weibull
(NBPF-Weibull) by Tang et al. (2024)[12], GE Weibull by kamal et al. (2023)[13] and
extended Xgamma by Saha et al. (2019)[14], which serve the purpose of flexible data
modelling. In this study, we introduce a new one-parameter continuous type of distribu-
tion using a similar approach used by Messaadia and Zeghdoudi (2018) [15] but taking a
different mixing probability.

The rest of the paper is organized as follows. In Section 3, we define the probability
density (pdf) and cumulative distribution (cdf) functions of the proposed distribution.
Section 4 derives its statistical properties, including raw moments, mean deviation, mo-
ment and cumulant generating functions, skewness, and kurtosis. The stochastic ordering
that characterizes the structural behavior of the distribution is discussed in Section 5.
Reliability properties, such as the hazard function and mean residual life, are presented in
Section 7. Section 8 provides results for the order statistics. Parameter estimation meth-
ods, including the method of moments and maximum likelihood estimation, are described
in Section 9. A Monte Carlo simulation study investigating bias and mean square error
of the estimators is given in Section 10. Applications of the proposed distribution to four
real-life datasets are discussed in Section 11. Finally, Section 12 concludes the paper with
a summary of findings and potential directions for future research.

2. Methodology

The new distribution is a finite mixture of two Gamma distributions having differ-
ent shape parameters with certain mixing probabilities. The form of the finite mixture
distribution is given by

g(x) =
k∑

i=1

πig(xi).

Let g1(x) ∼ Gamma(2, θ) and g2(x) ∼ Gamma(3, θ) and π1 =
θ

1 + θ
be the mixing

probability, then g(x) is our desired probability distribution. Although the Mgamma
distribution may be a particular form of one-parameter polynomial Exponential (OPPE)
distribution introduced by Bouchahed and Zeghdoudi (2018) [16] taking r = 2 and a0 =
0, a1 = 1, a2 = 1/2, it is easy to use by non-statisticians due to its flexible model. Whereas
for OPPE, it is difficult to determine which values of ai yield a good model, i.e., choice
of ai’s are very difficult to estimate. Moreover, it is a general form, difficult for non-
statisticians to use for modeling. Also, Bouchahed and Zeghdoudi (2018) did not test the
bias and MSE, they only provided the formula.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 3 of 22

3. The Mgamma distribution

Most of the one-parameter continuous distributions, such as Lindley, Shanker and
Xgamma etc., are derived directly from simple transformations or mixtures of Exponential
and Gamma with a single mixing scheme, whereas Mgamma is a finite mixture of two
Gamma distributions with different shape parameters but a single scale parameter, which
makes it different and richer than other classical distributions. A random variable X is
said to have the Mgamma distribution with parameter θ if its pdf is defined as

f(x) =

θ3x

(
1 +

x

2

)
e−θx

1 + θ
; x > 0, θ > 0 (1)

0 2 4 6 8 10

0.0
0.2

0.4
0.6

0.8
1.0

Plot of pdf of Mgamma distribution

x

f(x)

θ = 0.5

θ = 0.8

θ = 1

θ = 3

θ = 5

Figure 1: p.d.f. plot of Mgamma distribution for different θ

denoted by x ∼ M(θ), its cdf is given by

F (x) = 1−
[
1 +

θx

1 + θ

(
1 + θ +

θx

2

)]
e−θx; x > 0, θ > 0 (2)

To find the mode of the distribution, we first differentiate equation (1), w. r. t θ

d

dx
f(x; θ) =

θ3

1 + θ
[1 + x− θx(1 + x/2)] e−θx (3)

We observe that for θ < 1 ,
(1− θ) +

√
θ2 + 1

θ
is the unique critical point for which f(x) is

maximized. And for θ > 1 , f(x) is decreasing.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 4 of 22

0 2 4 6 8 10

0.0
0.2

0.4
0.6

0.8
1.0

Plot of cdf of Mgamma distribution

x

F(x
)

θ = 0.5

θ = 0.8

θ = 1

θ = 3

θ = 5

Figure 2: c.d.f. plot of Mgamma distribution for different θ

So, the mode(x) of the new distribution is given by

mode(x) =
(1− θ) +

√
θ2 + 1

θ
; θ < 1 (4)

Remark 1. From the above it is seen that mode(x) < median(x) < mean(x) which is
same as in the Xgamma and Exponential distribution. We get value of mode for different
values of θ

4. Moments and related measures

Let x ∼ Mgamma(θ), then the kth raw moment is given by

E(Xk) =

∫ ∞

0
xk

θ3x
(
1 + x

2

)
e−θx

θ + 1
dx.

=
θ3

θ + 1

∫ ∞

0
xk+1

(
1 + x

2

)
e−θx dx

=
θ3

θ + 1

[∫ ∞

0
xk+1e−θx dx+ 1

2

∫ ∞

0
xk+2e−θx dx

]
=

θ3

θ + 1

[
Γ(k + 2)

θk+2
+ 1

2

Γ(k + 3)

θk+3

]
=

θ3

θ + 1

[
(k + 1)!

θk+2
+ 1

2

(k + 2)(k + 1)!

θk+3

]
=

θ3(k + 1)!

θk+3(θ + 1)

[
θ + k+2

2

]



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 5 of 22

=
(k + 1)!

θk(θ + 1)

(
θ +

k + 2

2

)
. (5)

Thus, the kth order raw moment is given by

E(Xk) =
(k + 1)!

θk(θ + 1)

(
θ +

k

2
+ 1

)
, θ > 0, k > −2.

Therefore putting k = 1, 2, 3 and 4, we get the first, second, third and fourth order
moments about the origin as:

µ1
′ =

3 + 2θ

θ(1 + θ)
= mean(x),

µ2
′ =

6(2 + θ)

θ2(1 + θ)
,

µ3
′ =

12(5 + 2θ)

θ3(1 + θ)

and

µ4
′ =

120(3 + θ)

θ4(1 + θ)
.

The variance is given by

µ2 = Var(X) = E(X2)−
(
E(X)

)2
=

6(θ + 2)

θ2(θ + 1)
−

(
2θ + 3

θ(θ + 1)

)2

=
6(θ + 2)

θ2(θ + 1)
− (2θ + 3)2

θ2(θ + 1)2

=
6(θ + 2)(θ + 1)

θ2(θ + 1)2
− (2θ + 3)2

θ2(θ + 1)2

=
6(θ + 2)(θ + 1)− (2θ + 3)2

θ2(θ + 1)2

=
6(θ2 + 3θ + 2)− (4θ2 + 12θ + 9)

θ2(θ + 1)2

=
6θ2 + 18θ + 12− 4θ2 − 12θ − 9

θ2(θ + 1)2

=
2θ2 + 6θ + 3

θ2(θ + 1)2
.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 6 of 22

µ3 = E(X3)− 3E(X)E(X2) + 2(E(X))3

=
12(2θ + 5)

θ3(θ + 1)
− 3 · 2θ + 3

θ(θ + 1)
· 6(θ + 2)

θ2(θ + 1)
+ 2 · (2θ + 3)3

θ3(θ + 1)3

=
12(2θ + 5)(θ + 1)2 − 18(2θ + 3)(θ + 2)(θ + 1) + 2(2θ + 3)3

θ3(θ + 1)3

=
2
(
2θ3 + 9θ2 + 9θ + 3

)
θ3(θ + 1)3

µ4 = E(X4)− 4E(X3)E(X) + 6(E(X2)E(X))− 3E(X)4

=
24θ4 + 144θ3 + 252θ2 + 180θ + 45

θ4(1 + θ)4

Again, putting the value of (µ2) and (µ4) we get kurtosis (γ2) as,

γ2 =
µ4

µ2
2

=

24θ4 + 144θ3 + 252θ2 + 180θ + 45

θ4(θ + 1)4

(2θ2 + 6θ + 3)2

θ4(θ + 1)4

=
3(8θ4 + 48θ3 + 84θ2 + 60θ + 15)

4θ4 + 24θ3 + 48θ2 + 36θ + 9

=
24θ4 + 144θ3 + 252θ2 + 180θ + 45

(2θ2 + 6θ + 3)2

Figure 3: The graph of kurtosis for different θ



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 7 of 22

Now, putting the value of (µ2) and (µ3) we get the value of skewness (γ1) as ,

γ1 =
µ3

µ
3/2
2

=

2
(
2θ3+9θ2+9θ+3

)
θ3(θ+1)3

(2θ2+6θ+3)3/2

θ3(θ+1)3

=
2
(
2θ3 + 9θ2 + 9θ + 3

)
(2θ2 + 6θ + 3)3/2

.

Figure 4: The graph of skewness for different θ



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 8 of 22

From the Fig. 3 and Fig. 4 we can conclude the following observations,

• As the numerator and denominator are positive for all θ > 0 in (γ1), so γ1 > 0 i.e.
the distribution is right-skewed.

• Again, for kurtosis, we found that (γ2 > 3), ∀θ > 0, so the distribution is leptokurtic
(heavy-tailed, more peaked than normal).

4.1. Mean Deviation

The mean deviation about mean of random variable X, having pdf in equation (1) is
obtained as

MDx̄ = E(|x− µ|)

=

∫ ∞

0
|x− µ|f(x)dx

=

∫ µ

0
|µ− x|f(x)dx+

∫ ∞

µ
|x− µ|f(x)dx

= 2µF (µ)− µ+ 2

∫ ∞

µ
xf(x)dx

= 2µF (µ)− µ+
1

θ(1 + θ)

(
2θ(θ2µ2 + 2θµ+ 2) + (θ3µ3 + 3θ2µ2 + 6θµ+ 6)

)
e−θµ

4.2. Moment and Cumulant Generating Function

It is also possible to express the moment generating function in terms of the moment

MX(t) = E(etX) =

∫ ∞

0
etx

θ3x(1 +
x

2
)e−θx

1 + θ

 dx

=
θ3

1 + θ

∫ ∞

0
x

(
1 +

x

2

)
e−(θ−t)xdx

=
θ3

1 + θ

(∫ ∞

0
xe−(θ−t)xdx+

∫ ∞

0

x

2
e−(θ−t)x

)
dx

=
θ3

1 + θ

(
Γ2

(θ − t)2
+

Γ2

2(θ − t)2

)
=

θ3(1 + θ − t)

(1 + θ)(θ − t)2

The cumulant generating function of the random variable X is given by

kX(t) = logMX(t)



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 9 of 22

KX(t) = log

(
θ3(1 + θ − t)

(1 + θ)(θ − t)2

)
= 3 log θ + log(1 + θ − t)− log(1 + θ)− 2 log(θ − t)

5. Stochastic interpretation of the parameter θ

Stochastic orders are important measures for comparing the behavior of random vari-
ables shown by Shaked and Santikumar (1995) [17]) Let X and Y are the two random
variables with cumulative distribution functions FX(.) and FY (.) respectively. Then X is
said to be smaller than Y in the

• Stochastic order (X ≺s Y ), ifFX(t) ≥ FY (t), ∀t

• Convex order (X ≺cx Y ) if for all convex functions ϕ and provided expectation exist,
E[ϕ(X) ≤ E[ϕ(Y )].

• Hazard rate order (X ≺hr Y ), ifhX(t) ≥ hY (t), ∀t.

• Likelihood Ratio order (X ≺lr Y ), if fX(t)
fY (t) is decreasing in t.

Remark 2. Likelihood ratio order ⇒ Hazard rate order ⇒ Stochastic order if E[X] =
E[Y] , then Convex order ⇐⇒ Stochastic order.

Theorem 1. Let Xi ∼ M(θi), i = 1, 2 be two random variables. If θ1 ≤ θ2, thenX1 ≺lr

X2, X1 ≺hr X2, X1 ≺s X2andX1 ≺cx X2

Proof. We have

fX1(t)

fX2(t)
=

θ31(1 + θ2)

θ32(1 + θ1)
e−(θ1−θ2)t

differentiating w.r.t t, we get

d

dt

fX1(t)

fX2(t)
= −(θ1 − θ2)e

−(θ1−θ2)t < 0, ∀θ1 ≥ θ2

Hence, we see that
fX1

(t)

fX2
(t) decreases in t, which holds the above theorem. (proved)

6. Incomplete moments and Inequality curves

The rth incomplete moment µr(t) for a random variable X with pdf f(x) and cdf F(x)
is given by Sen et al. (2018) [18])

µr(t) =

∫ t

0
xrf(x)dx



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 10 of 22

So, when x ∼ Mgamma(θ), the rth incomplete moment is obtained as

µr(t) =

∫ t

0
xr

θ3x(1 +
x

2
)e−θx

1 + θ

 dx

=
θ3

1 + θ

(∫ ∞

0
xr+1e−θxdx+

∫ ∞

0

1

2
xr+2e−θxdx

)
=

θ3

1 + θ

[
γ(r + 2, θt) +

1

2
γ(r + 3, θt)

]
where γ(α, x) =

∫ x
0 uα−1e−u du is lower incomplete Gamma function. When a non-negative

continuous random variable X has pdf f (x) and cdf F (x), the Lorenz and Bonferroni curves
are defined by

l(p) =
1

µ

∫ q

0
xf(x)dx

B(p) =
1

pµ

∫ q

0
xf(x)dx

respectively, where µ = E(X) and q = F 1(p) for 0 < p < 1. So, when x ∼ Mgamma(θ)
the Lorenz and Bonferroni curves are given by

l(p) =
θ4

3 + θ

[
γ(3, θq) +

1

2
γ(4, θq)

]

B(p) =
θ4

p(3 + θ)

[
γ(3, θq) +

1

2
γ(4, θq)

]
The Lorenz and Bonferroni curves confirm that Mgamma can capture income/wealth in-
equality patterns due to its flexible tail behavior, offering better interpretability compared
to simpler one-parameter models

7. Reliability Properties

In reliability engineering, a batch of products may contain both weaker items (e.g.,
Gamma(2, θ)) and stronger items (long-lasting Gamma(3, θ)). Similarly, in the medical
field, the patients’ survival times may be a mixture of fast progressors and slow progressors,
etc. A mixture allows non-monotone hazard rates, which are often observed in practice
but impossible with simple one-parameter models. So Mgamma distribution is very much
suitable in this field with single- parameter and mixture distributions.

The hazard rate of Mgamma distribution is

h(x) =
f(x)

1− F (x)



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 11 of 22

=
θ3x(1 +

x

2
)e−θx[

1 +
θx

1 + θ

(
1 + θ +

θx

2

)]
e−θx

=

θ3x

(
1 +

x

2

)
[
1 + θ + θx

(
1 + θ +

θx

2

)]
From Fig. 5, it is evident that unlike the Exponential and Lindley distributions, the

hazard rate of the Mgamma distribution is non-monotone.

(i) For small θ, h(x) is increasing.

(ii) For moderate θ, h(x) shows a bathtub shape.

(iii) For large θ, h(x) is monotone increasing like Lindley.

0 5 10 15 20

0
1

2
3

4
5

Hazard plot of Mgamma distribution

x

h(x
)

θ = 0.8

θ = 1

θ = 3

θ = 0.5

Figure 5: Hazard plot of Mgamma distribution for different values of θ

Mean remaining life

For a continuous random variable X with pdf f(x) and cdf F (x), the mean residual
life (mrl) function is defined as

m(x) = E[(X − x)|X > x]



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 12 of 22

=
1

1− F (x)

∫ ∞

x
[1− F (t)]dt

=

∫∞
x F̄ (t)dt

F̄ (x)

Where,

F̄ (t) =

(
1 +

θt

1 + θ
(1 + θ +

θt

2

)
e−θt

and, ∫ ∞

x
F̄ (t)dt =

1

1 + θ

(
(1 + θ)x+

2(1 + θ)

θ
+ x+

1

θ
+

θx2

2

)
(6)

For the Mgamma distribution the mrl function is given by

m(x) =
(1 + θ)x+

2(1 + θ)

θ
+ x+

1

θ
+

θx2

2
1 + θ + θx(1 + θ + θx/2)

Figure 6 shows mrl functions for different choices of θ.

0 2 4 6 8 10

0
1

2
3

4
5

Mean remaining life of Mgamma distribution

x

mr
l(x)

θ = 0.8

θ = 1

θ = 3

θ = 5

Figure 6: mrl plot of Mgamma distribution for different values of θ

From the plot we found that

(i) For x = 0, m(0) = µ =
2θ + 3

θ(1 + θ)
which is the mean of Mgamma.

(ii) The mrl(x) decreases with both x and θ, indicating that the expected remaining life
becomes shorter.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 13 of 22

8. Order Statistics

Let Y(1), Y(2), .., Y(n) be the ordered random sample drawn from the Mgamma Distri-

bution, then the c.d.f. and p.d.f. of the rth order statistics is given by

Gr(y) =

n∑
i=r

(
n

i

)
[G(y; θ;λ)]i [1−G(y; θ;λ)]n−i

gr(y) =
n!

(r − 1)!(n− r)!
[G(y; θ;λ)]r−1 [1−G(y; θ;λ)]n−r g(y; θ;λ)

So, the p.d.f. and c.d.f. of y(1) and y(n) is given by

gY1(y) =
nθ3x(1 + x/2)e−θx

1 + θ

[(
1 +

θx

1 + θ
(1 + θ +

θx

2
)

)
e−θx

]n−1

gYn(y) =
nθ3x(1 + x/2)e−θx

1 + θ

[
1−

(
1 +

θx

1 + θ
(1 + θ +

θx

2
)

)
e−θx

]n−1

Gy1(y) = 1−
[(

1 +
θx

1 + θ
(1 + θ +

θx

2
)

)
e−θx

]n
Gyn(y) =

[
1−

(
1 +

θx

1 + θ
(1 + θ +

θx

2
)

)
e−θx

]n
9. Estimation of the parameters

Here we discuss two estimation methods viz. (i) Method of moments and (ii) Methods
of maximum likelihood estimation.

9.1. Method of moment

Let x1, x2, x3, .., xn be a random sample of size n drawn from the equation (1). We

compare the sample mean x̄ =
1

n

∑n
i=1 xi with the first order raw moment µ1

′ to get the

mom estimator of θ as

θmom =
1

2x̄

(
−x̄+

√
x̄2 + 8x̄+ 4 + 2

)
, x > 0

The method of moment estimator θmom is positively biased, consistent and asymptot-
ically normal.

For large sample the 100(1-α)% confidence interval of θ is given by θ ± Zα

2

√
1

nσ2



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 14 of 22

9.2. Maximum likelihood estimation

Let X1, X2, ..., Xn be a random sample of size n drawn independently from the pdf 1.
The maximum likelihood estimator (mle) of θ is given by

L(θ|x1, x2, .., xn) =
n∏

i=1

θ3xi(1 + xi/2)e
−θxi

1 + θ

The log likelihood equation is as follows

L(θ|x) = 3nlog(θ)− nlog(1 + θ) +
n∑
i

log [xi(1 + xi/2)]− θ
n∑
i

xi (7)

To have the maximum likelihood estimator of θ, we have to solve I(θ|x) = ∂L(θ|x)
∂θ = 0,

Hence, we have

3n

θ
− n

1 + θ
−

n∑
i=1

xi = 0

Using Newton-Raphson iteration method, we find the θ̂mle for which the log-likelihood

equation is maximized. The initial solution for the iteration method is θ0 =
2

x̄
. We

continue the iteration process

θ(i) = θ(i−1) − I(θ(i−1)|x)
I ′(θ(i−1)|x)

we stop when θ(i) ≡ θ(i−1). All computations were performed in the updated version of R
software [19], following the procedures book [20]. The MLE is consistent and asymptot-
ically normal under standard regularity conditions (iid samples, differentiability of likeli-
hood). Convergence of Newton–Raphson is generally fast; in our simulation, it stabilized
within 4–6 iterations.

10. Simulation Study

To generate the random sample from Mgamma distribution, we proceed as follows:
Let F (x) = p, where p follows uniform(0, 1). We choose the solution that lies in (0, 1).
Here G(x) is the cdf of Mgamma random variable. Then, we follow the approach of
generating sample from Mgamma distribution. Hence the algorithm is

1. Generate pi ∼ uniform(0, 1), i = 1, 2, . . . , n.

2. Generate Vi ∼ Gamma(2, θ), i = 1, 2, . . . , n.

3. Generate Wi ∼ Gamma(3, θ), i = 1, 2, . . . , n.

4. If Ui ≤ θ
1+θ , then set Xi = Vi. Otherwise, set Xi = Wi.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 15 of 22

We perform Monte Carlo simulation in R programming considering N=10000 times for
different values of n=5, 10, 15 and 20 and θ=0.8, 1, 2 and 3 . The following two measures
were computed

1. Average bias of the simulated estimates θ̂i, i = 1, 2, . . . , N : 1
N

∑N
i=1(θ̂i − θ)

2. Average Mean Square Error (MSE) of the simulated estimates θ̂i, i = 1, 2, . . . , N :
1
N

∑N
i=1(θ̂i − θ)2

Figure 7: mse plot of Mgamma distribution for different value of θ

It is observed that both the MSE and bias decrease as the sample size increases and θ
decreases.

11. Data Analysis

In this section, we have compared our proposed model with the existing one-parameter
distributions by fitting some real-life datasets. The first dataset is obtained from Smith
and Naylor (1987) [21], which represents the strength of 1.5 cm glass fibres, measured at
the National Physical laboratory, England.

The second dataset is taken from Bakouch et. al (2021)[22] which concerns the ex-
ceedances of flood peaks (in m3/s) of the Wheaton River near Carcross in Yukon Territory,
Canada that consist of 72 exceedances for the years 1958–1984 and rounded to one decimal
place.

The third dataset is about the remission times (in weeks) for 20 leukaemia patients
randomly assigned to a certain treatment taken from Al-Babtain et al. (2020) [23].



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 16 of 22

Figure 8: bias plot of Mgamma distribution for different value of θ

The last dataset is obtained from Sharma et al. (2017) [24], which consists of the
waiting time (in mins) before service for 100 bank customers.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 17 of 22

Fitting to strength of 1.5 cm glass fibres

strength (in cms)

fr
e

q
u

e
n

c
y
 d

e
n

s
it
y

0 10 20 30 40 50 60 70

0
.0

0
0

.0
2

0
.0

4

Exponential
Xgamma
Lindley
Mgamma
Shanker
akash

Fitting exceedence of flood peaks 
 of the wheaton river

flood peaks (in m^3/s)

fr
e

q
u

e
n

c
y
 d

e
n

s
it
y

0 10 20 30 40 50 60 70

0
.0

0
0

.0
4

0
.0

8

Exponential
Xgamma
Lindley
Mgamma
Shanker
akash

Fitting of Remission times for 
 20 leukaemia patients

times (in weeks)

fr
e

q
u

e
n

c
y
 d

e
n

s
it
y

0 10 20 30 40 50

0
.0

0
0

.0
4

0
.0

8

Exponential
Xgamma
Lindley
Mgamma
Shanker
akash

Fitting waiting time before service 
  of 100 bank customers

time(in mins)

fr
e

q
u

e
n

c
y
 d

e
n

s
it
y

0 10 20 30 40

0
.0

0
0

.0
4

0
.0

8

Exponential
Xgamma
Lindley
Mgamma
Shanker
akash

Figure 9: Histogram and Fitted densities for all the four datasets.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 18 of 22

0

200

400

600

100 200 300 400

Theoretical Quantiles (Mgamma)

E
m

p
ir

ic
a

l 
Q

u
a

n
ti
le

s
 (

D
a

ta
)

QQ Plot for Mgamma Distribution

0

20

40

60

0 10 20 30

Theoretical Quantiles (Mgamma)

E
m

p
ir

ic
a

l 
Q

u
a

n
ti
le

s
 (

D
a

ta
)

QQ Plot for Mgamma Distribution

0

10

20

30

40

50

10 20 30 40 50

Theoretical Quantiles (Mgamma)

E
m

p
ir

ic
a

l 
Q

u
a

n
ti
le

s
 (

D
a

ta
)

QQ Plot for Mgamma Distribution

0

10

20

30

40

0 10 20 30

Theoretical Quantiles (Mgamma)

E
m

p
ir

ic
a

l 
Q

u
a

n
ti
le

s
 (

D
a

ta
)

QQ Plot for Mgamma Distribution

Figure 10: Q-Q plot for all the four datasets using Mgamma distribution accordingly.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 19 of 22

Table 1: Model comparison for Dataset I: Strength of 1.5 cm glass fibres

Criteria Exponential Lindley Xgamma Shanker Akash Mgamma

MLE 0.6636 0.9961 1.3376 0.9562 1.3554 1.8354
AIC 88.83 81.27 85.80 81.14 81.86 58.66
KS Statistic 0.4186 0.3862 0.4196 0.3661 0.3705 0.3361
Cramér–von Mises 3.862 3.332 3.881 3.113 3.192 2.496
Anderson–Darling 18.426 16.425 18.468 15.406 15.736 12.590

Table 2: Model comparison for Dataset II: Exceedances of flood peaks (Wheaton River)

Criteria Exponential Lindley Xgamma Shanker Akash Mgamma

MLE 0.8194 0.1530 0.2044 0.1646 0.2412 0.2422
AIC 252.120 264.228 266.470 279.790 291.070 251.648
KS Statistic 0.1421 0.2408 0.2562 0.2272 0.3030 0.1325
Cramér–von Mises 0.2319 0.8198 1.0660 0.7775 1.5190 0.1489
Anderson–Darling 1.4680 7.4400 8.8136 6.9887 19.100 1.4231

Table 3: Model comparison for Dataset III: Remission times of 20 leukaemia patients

Criteria Exponential Lindley Xgamma Shanker Akash Mgamma

MLE 0.0511 0.0977 0.1383 0.1019 0.1522 0.1521
AIC 79.459 78.872 79.398 79.244 81.790 78.746
KS Statistic 0.1143 0.1189 0.1589 0.1095 0.2007 0.0104
Cramér–von Mises 0.0660 0.0344 0.0725 0.0424 0.1366 0.0181
Anderson–Darling 0.410 0.2960 0.5644 0.3351 1.4800 0.2704

Table 4: Model comparison for Dataset IV: Waiting times of 100 bank customers

Criteria Exponential Lindley Xgamma Shanker Akash Mgamma

MLE 0.1012 0.1865 0.2634 0.1983 0.2952 0.2989
AIC 329.020 319.037 321.020 317.629 320.629 317.062
KS Statistic 0.1730 0.0676 0.0624 0.0881 0.1002 0.0601
Cramér–von Mises 0.7153 0.0581 0.0839 0.1280 0.2170 0.0535
Anderson–Darling 4.2280 0.4863 0.6447 0.9044 1.3797 0.0326



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 20 of 22

From the table 1 - 4, we find that Mgamma provides the best fit, capturing the skewness
and heavy tails of the data. Across all four datasets i.e. spanning engineering, hydrology,
medical survival, and queueing contexts the Mgamma distribution consistently outper-
formed the classical one-parameter models. The mixture structure of Mgamma allows it
to model skewness and heavy-tailed behaviour more effectively. Also the model accommo-
dates non-monotone hazard rates (e.g., bathtub-shaped or increasing-decreasing), which
are common in practice but not captured by simpler distributions. Even with small sam-
ples (as in the remission time data), Mgamma provides superior fit. In every dataset,
Mgamma yielded the lowest AIC and other criteria, reinforcing its statistical adequacy.
Thus, the Mgamma distribution can be considered a highly competitive and flexible alter-
native for real-life data modelling, outperforming several classical and recently proposed
lifetime distributions. It is well known that distributions with more parameters usually
provide greater flexibility for data modeling. However, the presence of multiple parameters
often makes estimation more complex and computationally demanding. In contrast, the
proposed Mgamma distribution has the advantage of involving only a single parameter,
which is straightforward to estimate, while still offering substantial flexibility. Notably, it
provides a superior fit compared to several existing multi-parameter models.

12. Conclusion

In this article, we introduced a new one-parameter continuous distribution, the so
called Mgamma distribution, by a finite mixture approach of two Gamma distributions
with different shape parameters. We studied several of its statistical properties including
moments, moment and cumulant generating function, order statistics, mean deviation,
incomplete moments and inequality curves. We generated random variables and studied
their mean square error and bias and found that both decrease as the sample size increases.
The estimation of model parameters is approached by the method of maximum likelihood
and method of moments estimation. Two real-life data sets are analyzed using various
model selection criteria. It is found that the new distribution provides a consistently better
fit than the other probability models available in the literature. For future studies, this
new distribution can be modified using a different technique to make it more flexible in
data modelling, such as developing its pseudo form, considering the ratio of two random
variables, introducing transmuted versions with an additional parameter, or exploring
Bayesian estimation approaches.

Data availability

The data that support the findings of this study are available from the corresponding
author upon reasonable request.

Conflict of interest

There is no conflict of interest.



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 21 of 22

Acknowledgements

The authors are grateful to the Deanship of Graduate Studies and Scientific Re-
search, Islamic University of Madinah, Saudi Arabia, for the support provided to the
Post-Publication Program (4). The authors also thank the unknown referees and editors
Eyup Cetin and Baris Kiremitci for their valuable comments and suggestions to improve
the quality of the article.

References

[1] Geoffrey J McLachlan and Kaye E Basford. Mixture models: Inference and applica-
tions to clustering, volume 38. M. Dekker New York, 1988.

[2] Mohamed E Ghitany, Barbra Atieh, and Saralees Nadarajah. Lindley distribution
and its application. Mathematics and computers in simulation, 78(4):493–506, 2008.

[3] Sudhansu S Maiti, Molay Kumar Ruidas, and Sumanta Adhya. Natural discrete one
parameter polynomial exponential family of distributions and the application. Annals
of Data Science, 11(3):1051–1076, 2024.

[4] Subhradev Sen, Sudhansu S Maiti, and N Chandra. The xgamma distribution: sta-
tistical properties and application. Journal of Modern Applied Statistical Methods,
15(1):38, 2016.

[5] Dennis V Lindley. Fiducial distributions and bayes’ theorem. Journal of the Royal
Statistical Society. Series B (Methodological), pages 102–107, 1958.

[6] R Shanker. Shanker distribution and its applications. International journal of statis-
tics and Applications, 5(6):338–348, 2015.

[7] Halim Zeghdoudi and Sihem Nedjar. Gamma lindley distribution and its application.
Journal of Applied Probability and Statistics, 11(1):129–138, 2016.

[8] Norman Lloyd Johnson, Samuel Kotz, and Narayanaswamy Balakrishnan. Continuous
univariate distributions. 1994.

[9] Hassan S Bakouch, Bander M Al-Zahrani, Ali A Al-Shomrani, Vitor AA Marchi,
and Francisco Louzada. An extended lindley distribution. Journal of the Korean
Statistical Society, 41(1):75–85, 2012.

[10] Haitham M Yousof, Mustafa Ç Korkmaz, and Subhradev Sen. A new two-parameter
lifetime model. Annals of Data Science, 8:91–106, 2021.

[11] AAM Mahmoud, Rana H Khashab, Zakiah I Kalantan, Sulafah M Saleh Binhimd,
Amani S Alghamdi, and Said G Nassr. On bivariate compound exponentiated survival
function of the beta distribution: Estimation and prediction. Journal of Radiation
Research and Applied Sciences, 17(2):100886, 2024.

[12] Xiangming Tang, Jin-Taek Seong, Randa Alharbi, Aned Al Mutairi, and Said G Nasr.
A new probabilistic model: Theory, simulation and applications to sports and failure
times data. Heliyon, 10(4), 2024.

[13] Mustafa Kamal, Meshayil M Alsolmi, Aned Al Mutairi, Eslam Hussam, Man-
ahil SidAhmed Mustafa, Said G Nassr, et al. A new distributional approach: estima-



M. I. Khan, M. K. Ruidas / Eur. J. Pure Appl. Math, 18 (4) (2025), 6653 22 of 22

tion, monte carlo simulation and applications to the biomedical data sets. Networks
& Heterogeneous Media, 18(4), 2023.

[14] Mahendra Saha, Abhimanyu Singh Yadav, Arvind Pandey, Shivanshi Shukla,
and Sudhansu S Maiti. The extended xgamma distribution. arXiv preprint
arXiv:1909.01103, 2019.

[15] Hamouda Messaadia and Halim Zeghdoudi. Zeghdoudi distribution and its appli-
cations. International Journal of Computing Science and Mathematics, 9(1):58–65,
2018.

[16] Lahsen Bouchahed and Halim Zeghdoudi. A new and unified approach in generalizing
the lindley’s distribution with applications. Statistics, 61, 2018.

[17] Moshe Shaked, J George Shanthikumar, and YL Tong. Stochastic orders and their
applications. SIAM Review, 37(3):477–478, 1995.

[18] Subhradev Sen, N Chandra, and Sudhansu S Maiti. On properties and applications of
a two-parameter xgamma distribution. Journal of Statistical Theory and Applications,
17(4):674–685, 2018.

[19] R Core Team. R: A Language and Environment for Statistical Computing. R Foun-
dation for Statistical Computing. Vienna, Austria, 2025.

[20] Ross Ihaka and Robert Gentleman. R: a language for data analysis and graphics.
Journal of computational and graphical statistics, 5(3):299–314, 1996.

[21] Richard L Smith and J. C. Naylor. A comparison of maximum likelihood and bayesian
estimators for the three-parameter weibull distribution. Journal of the Royal Statis-
tical Society: Series C (Applied Statistics), 36(3):358–369, 1987.

[22] Hassan S Bakouch, Abdus Saboor, and Muhammad Nauman Khan. Modified beta
linear exponential distribution with hydrologic applications. Annals of Data Science,
8(1):131–157, 2021.

[23] Abdulhakim A Al-Babtain, Abdul Hadi N Ahmed, and Ahmed Z Afify. A new discrete
analog of the continuous lindley distribution, with reliability applications. Entropy,
22(6):603, 2020.

[24] Vikas Kumar Sharm, Sanjay Kumar Singh, and Umesh Singh. Classical and bayesian
methods of estimation for power lindley distribution with application to waiting time
data. Communications for statistical applications and methods, 24(3):193–209, 2017.