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

Parametric Estimation for a New Pareto-Type Model
Based on Constant-Stress Partially Accelerated

Censoring Data

Ahmed A. Soliman1, Gamal A. Abd-Elmougod2, Osama M. Taha1,∗,
Al-Wageh A. Farghal1

1 Department of Mathematics, Faculty of Science, Sohag University, Sohag 82524, Egypt
2 Department of Mathematics and Computer Science, Faculty of Science,
Damanhour University, Damanhour, Egypt

Abstract. In Low failure or High longevity, accelerated life tests (ALTs) are used to speed up
tests. This paper investigates precise estimation issues related to point and interval estimations
for a new Pareto-type (NPT) distribution. The data are collected under constant stress partially
(ALTs) concerning Type-I generalized hybrid censored samples (Type-I GHCS). The point max-
imum likelihood estimates (MLEs) of the model parameters and accelerated factor are obtained
with the help of the expectation–maximization (EM) algorithm. Also, Bayes estimates under
various loss functions with the help of the Metropolis-Hastings (MH) algorithm method are con-
structed. The asymptotic confidence intervals, bootstrap confidence interval, and highest posterior
density (HPD) credible intervals are derived. The performance of different estimators is compared
with the help of a Monte Carlo simulation study. Finally, we analyze a real data set to show the
applicability of the model considered.

2020 Mathematics Subject Classifications: 62F10, 62F12, 62F15, 62F40, 62N02

Key Words and Phrases: A new Pareto-type distribution, Accelerated life tests model, Type-I
generalized hybrid censoring scheme, EM algorithm, Maximum likelihood estimation, Bayesian
estimation

1. Introduction

There are many real-life scenarios where data demand a probability distribution with
decreasing and/or upside down bathtub (unimodal) shaped failure rate functions. The
probability that an event occurs at a given moment, given that it has not yet occurred,
decreases over time a phenomenon known as a decreasing failure rate. Patients undergoing

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5788

Email addresses: a a sol@hotmail.com (A. A. Soliman),
gam amin@yahoo.com (G. A. Abd-Elmougod), osama.taha@science.sohag.edu.eg (O. M. Taha),
alwageh ahmed@science.sohag.edu.eg (Al-Wageh A. Farghal)

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



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 2 of 34

heart transplants, for example, face an increasing hazard of death in the first few days
after the transplant while the body adjusts to the new organ. As the patient improves,
the hazard rate reduces Ref.[1]. The upside-bathtub-shaped failure rate function would
be appropriate in this case. To represent a model of failure rate with a relatively high
rate of failure in the middle of the predicted lifetime, an unimodal hazard rate function
is used. When product failures are driven by fatigue and corrosion, the related failure
rates have unimodal forms Ref.[2]. In addition, in some medical scenarios, such as breast
cancer and infection with novel viruses, the hazard rate is unimodal, as demonstrated by
Ref.[3] and Ref.[4]. The NPT distribution exhibits both a decreasing and an upside-down
bathtub (unimodal) shaped failure rate function Ref.[5]. NPT distribution is considered a
generalization of the well-known Pareto distribution. When the ”probability” or fraction
of the population that owns a small amount of wealth per person is relatively high, and
then steadily decreases as wealth increases, as in Ref.[6], experimenter seeks a probabil-
ity distribution describing this case, the NPT model is an excellent alternative to the
most common model Pareto distribution. Because it is used to analyze various income
and reliability data, this distribution might be a versatile model. It is used in actuar-
ial sciences, reliability, finance, and climatology to describe the occurrence of extreme
weather events. In light of these characteristics, multiple authors adopted this model in
data modeling, particularly censored observations. As a result, our motivation to validate
this distribution throughout this article stems from its practical usefulness in the various
areas outlined above, as well as the supporting references. Among these references, [7]
proposed the comparisons methods of estimation for the model, Ref.[8] studied the ML,
Bayes estimation and prediction of the model under progressive type-II censoring. Also,
Ref.[9] derived simpler expressions for many relevant economic inequality and risk indices
using the incomplete beta function, Ref.[10] presented a generalized of the NPT distribu-
tion. Recently, Ref.[11] developed inference based on NPT records with applications to
precipitation and Covid-19 data. The NPT distribution, denoted by NPT(θ, σ) is specified
by the following Probability Density Function (PDF), Cumulative Distribution Function
(CDF), Survival Function (SF), and Hazard Rate function (HRF) given by

f(x) =
2θσθxθ−1

(xθ + σθ)2
; x ≥ σ, θ, σ > 0

F (x) = 1− 2σθ

xθ + σθ

S(x) =
2σθ

xθ + σθ

h(x) =
θxθ−1

xθ + σθ
.

(1)

Different potential scenarios concerning the reliability and quality of items face the
challenge of lacking adequate knowledge about the failure of such products under normal
operating conditions. To address this, experimenters uses accelerated life testing (ALT)
or partially accelerated life testing (PALT) to quickly gather sufficient failure data and



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 3 of 34

understand the relationship between failures and external stress factors. These tests can
significantly save time, labor, resources, and costs. In ALTs, all units are subjected to
stress levels higher than usual to induce failures more quickly. In contrast, PALTs involve
testing units under both normal and elevated stress conditions in order to collect more
failure data within a constrained time frame without subjecting all units to high stress.
PALTs are particularly beneficial in scenarios where time and cost are pressing concerns.
The information obtained from such tests can be utilized to estimate the failure behavior
of the units under normal conditions. Ref.[12] performed a proper statistical modeling of
PALTs, where the authors considered the tempered random variable model for PALTs. In
contrast, with a step stress ALT model, the stress level changes at specific times or when
a designated number of units fail. As a result, the stress on all the products increases
gradually according to some type of rule. Many authors have extensively researched the
step stress models. For recent literature on this topic, Refs. [13],[14], and [15]. To bet-
ter extrapolate the lifetime under normal use conditions, partially accelerated life test
(PALT) was developed. It divides all products into groups and test each in normal and
accelerated environments. Two sorts of prevalent PALTs are known as the constant stress
partially accelerated life tests (CSPALTs) and the step stress partially accelerated life tests
(SSPALTs). CSPALTs set different groups of units under use and stress condition respec-
tively. By comparison, SSPALTs start with normal conditions and switch to accelerated
conditions if the unit does not cease to function before a prefixed point in time τ . The
specific illustration of the differences between the two PALTs are displayed in Fig.(1).

Figure 1: A schematic diagram of (a) CSPALTs and (b) SSPALTs

Our primary focus in this paper lies within the realm of CSPALTs. Numerous re-
search studies have also given attention to CSPALTs, as illustrated by, for instance,
Refs.[16],[17],[18],[19],[20], and [21]. Even though the major goal of partial ALT is to
reduce the testing period of the experiment, the experimenter spends a significant amount



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 4 of 34

of time waiting for all test units to reach the point of failure. To overcome this problem,
censoring schemes are employed to terminate life tests. Type-I and Type-II censoring are
commonly regarded as two main scheme for conducting lifetime experiments, in which we
terminate the experiments at a specific time point or when a given amount of failures
occur. Ref.[22] merged the two censoring techniques mentioned earlier and introduced the
initial hybrid censoring scheme, referred to as the Type-I hybrid censoring scheme (Type-I
HCS). Within the Type-I HCS , the life testing experiment concludes either at the des-
ignated time T or when a predetermined count of m items encounter failure. Therefore,
the cessation moment for Type-I HCS is given by T ∗ = min(Xm:n, T ), where Xm:n rep-
resents the duration until m out of n testing items experience failure. Many researchers
have worked on Type-I HCS, including Ref.[23]. One downside of Type-I HCS is that a
small number of failures may occur till after a set duration T ∗. Hence, drawing statisti-
cal conclusions within this type of scheme might prove to be unworkable. To address this
drawback and enhance the effectiveness of estimators within life-testing experiments, while
also ensuring a specific count of failures before the experiment concludes, and concurrently
reducing testing time and costs associated with unit failures, a generalized hybrid Type-I
censoring scheme was introduced by Ref.[24]. The Type-I Generalized Hybrid Censoring
scheme (Type-I GHCS) approach ensures a minimum threshold of failures, which helps
alleviate the shortfall present in Type-I HCS. We suppose, {X1:n, . . . , Xn:n} constitute a
sequence of ordered observations concerning the lifetimes of failures within a set of size n.
The values of l, m, and T are defined in advance, where l < m < n represents the small-
est acceptable number of failures fixed before the experiment, m represents the intended
quantity of failures, and T represents a certain time instant. These three censoring models
mentioned previously can be stated as follows;

• Type-I CS: terminate at T .

• Type-I HCS: terminate at T ∗ = min(Xm,n, T ).

• Type-I GHCS: terminate at T ∗ = max(Xl:n,min(Xm:n, T )).

This article focuses on Type-I GHCS, which is classified into three categories such as;

• {X1:n < · · · < Xl:n}, When, Xl:n > T .

• {X1:n < · · · < Xd∗:n}, When, Xl:n < T < Xm:n.

• {X1:n < · · · < Xm:n}, When, Xm:n < T .

Additional literature on Type-I GHCS has recently been offered, including Refs.[25],[26],[27],
[28] and, [29]. The purpose of this article is parameters estimation of the NPT distribu-
tion based on a partially constant-stress ALT under Type-I GHCS. Maximum likelihood
(ML) and Bayesian methods will be used to estimate the parameters of the NPT distri-
bution. The Newton-Raphson (NR) and EM algorithms are developed and described in
detail when constructing maximum likelihood estimates (MLEs) and their corresponding
confidence intervals. When the Bayes approach is applied, the posterior distribution and



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 5 of 34

the corresponding Bayes estimate often require computing integrals, which can be dif-
ficult, especially when using complex high- or low-dimensional models. We applied the
MCMC method to estimate the unknown parameters under squared error loss function
with dependent priors for the parameters. The Gibbs within Metropolis samples can be
used to produce samples from posterior distributions using various MCMC algorithms.
Extensive computer studies can be used to compare the performance of Bayes estimators
against that of classical MLEs. Using the NR, EM, and MCMC algorithms, approximate
95% confidence intervals for unknown parameters can be generated. We will also compare
them based on their average lengths and coverage probability. This article is structured
as follows: In Section (2), the proposed model and its fundamental assumptions are dis-
cussed. Section (2), we display the point estimate with, maximum likelihood estimation
by using Newton-Raphson and expectation maximization method and Bayesian estimates
with the help of MCMC method. Section (2) deals with interval estimates. Approximate
confidence intervals, bootstrap confidence intervals, and credible intervals for the unknown
parameters are also constructed. The findings of the simulation investigation are outlined
in Section (3). A real-life example and the simulated example are presented and discussed
in Section (3). Finally, some concluding remarks are given in Section (4).

2. Methodology

In this particular instance, unit lifetime has NPT distribution in the model that is
being developed. The point estimates of Parameters are created using the Bayesian, EM
algorithm and MLE approaches. Additionally, interval estimators are developed using
bootstrap methods, HPD credible intervals, and the asymptotic property of MLEs.

2.1. Modeling

Suppose in a partially constant-stress ALT model under Type-I GHCS, n indicates
that items are separated into two groups with sizes n1 and n2 = n − n1. The normal
condition is assigned to items n1, and the stress condition is assigned to items n2. The
hazard failure rate of an item under the stress condition is given by h2(X) = λh1(X),
where h1(X) is the FR function under normal conditions, and λ > 1 is the acceleration
factor. Suppose that the lifetime of an item follows the two-parameter NPT distribution.
As a result of Eq.(1), the PDF, CDF, and FR function of the lifetime X under accelerated
conditions are given respectively by;

f2(x) =
2λθλσθλxθ−1

(xθ + σθ)1+λ
; x ≥ σ

F2(x) = 1− 2λσθλ

(xθ + σθ)λ

h2(x) =
θλxθ−1

xθ + σθ
.

(2)



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 6 of 34

Let, Xj = {Xj1, . . . , Xjdj}, j = 1, 2, and X = {X1, X2}, be the set of data under normal
and accelerate conditions, respectively. Also, suppose fixed integers 1 < lj < mj < nj and
fixed time T ∈ (0,∞). If the lthj failure occurs before time T , terminate the experiment at

min{Xmj :nj , T}. If the lthj failure occurs after time T , terminate the experiment at Xlj :nj
.

Under this setting, the experimenter would ideally like to observe mj failures but is willing
to accept a bare minimum of lj failures. The realization of the failure time sample can
be presented by X(j1) ≤ X(j2) ≤ · · · ≤ X(jDj), Dj represents the number of failures before
time point T , where;

Dj =


lj if Xlj > T

d∗j if Xlj < T < Xmj

mj if Xmj < T

Moreover, let Cj denote the recorded lifetimes for all survived components (terminated
time), then Cj can also be expressed as;

Cj =


Xlj if Xlj > T

T if Xlj < T < Xmj

Xmj if Xmj < T

For given type-I GHC sample under partially constant-stress ALT,Xj = {Xj1, . . . , Xjdj}, j =
1, 2 the joint likelihood function can be written as

L(Φ|X) =

2∏
j=1

nj !

(nj −Dj)!
(1− Fj(Cj |Φ))−(Dj−nj)e

∑Dj
i=1 log(fj(xi|Φ)), (3)

where Φ = {θ, σ, λ}.

2.2. Point Estimation

In this subsection, we discuss ML and Bayes point estimations under accelerated type-
I GHC sample from NTP distribution.In the ML approach, we estimate the parameters
using the Newton-Raphson and EM algorithms. In the Bayesian approach, we use the
MCMC method.

2.2.1. ML estimation

The purpose of ML estimation is to identify the values of model parameters that optimize
the likelihood function across the whole parameter space. The ML technique is a com-
monly employed statistical inference technique that may be applied to a diverse range of
distributions and models. The ML estimator possesses several key attributes, including
consistency, asymptotic efficiency, asymptotic unbiasedness, and asymptotic normality.
These quantiles have contributed to the widespread adoption and significance of the MLE



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 7 of 34

as a prominent technique for statistical model fitting. The related PDF and CDF are
given in Eqs.(1) and (2). Therefore, the likelihood function is derived as shown below;

L(θ, σ, λ|X) ∝ 2λn2θD1+D2λD2σθ(D1+λD2)

[
1 + (

C1

σ
)
θ ](D1−n1)[

1 + (
C2

σ
)
θ ]λ(D2−n2)

×
2∏

j=1

Dj∏
i=1

x
−(Sj(λ) θ+1)
ji[

1 + ( σ
xji

)θ
]Sj(λ)+1

, σ ≤ min{x11, x21}, (4)

and

Sj(λ) =

{
1, When, j = 1

λ, When, j = 2.

The ML estimators of the parameters θ, σ, and λ are the values that maximize the
likelihood function, as shown in Eq.(4).
Maximizing the likelihood function presents challenges, so it is more advantageous to
optimize its natural logarithm. The natural logarithm associated with Eq.(4), represented
as ℓ(θ, σ, λ|X), can be expressed as;

ℓ(θ, σ, λ|X) ∝ n2λ log 2 + (D1 +D2) log θ +D2 log λ+ θ(D1 + λD2) log σ

+ (D1 − n1) log (1 + (
C1

σ
)
θ

) + λ(D2 − n2) log (1 + (
C2

σ
)
θ

)

−
2∑

j=1

(Sj(λ)θ + 1)

Dj∑
i=1

log xji −
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

log(1 + (
σ

xji
)
θ
). (5)

The ML estimate of σ is σ̂ML = x1. By calculating the first derivatives of Eq.(5) with
respect to θ and λ, and then setting them equal to zero, we obtain the following system
of simultaneous equations.

∂ℓ∗(θ, λ|X)

∂θ
=

D1 +D2

θ
+ (D1 + λD2) log x1 +

(D1 − n1)(
C1
x1
)
θ
log (C1

x1
)

1 + (C1
x1
)
θ

−
2∑

j=1

Sj(λ)

Dj∑
i=1

log xji

+
λ(D2 − n2)(

C2
x1
)
θ
log (C2

x1
)

1 + (C2
x1
)
θ

−
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

( x1
xji

)θ log ( x1
xji

)

1 + ( x1
xji

)θ
= 0, (6)

and

∂ℓ∗(θ, λ|X)

∂λ
= n2 log 2 +

D2

λ
+D2 θ log x1 + (D2 − n2) log (1 + (

C2

x1
)
θ

)− θ

D2∑
i=1

log x2i

−
D2∑
i=1

log(1 + (
x1
x2i

)
θ
) = 0. (7)



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 8 of 34

Using Eq.(7) and for fixed θ, one can obtain the MLE of the unknown parameter λ as a
function of the parameter θ as shown below;

λ̂ML(θ) =
D2

−n2 log 2−D2θ log x1−(D2−n2) log (1+(
C2
x1

)
θ
)+θ

∑D2
i=1 log x2i+

∑D2
i=1 log(1+(

x1
x2i

)
θ
)

(8)

By substituting λ̂ML(θ) in the normal equation given by Eq.(6), the ML estimate of θ can
be obtained by solving the following non-linear equation;

∂ℓ∗(θ, λ|X)

∂θ
=

D1 +D2

θ
+ (D1 + λ̂ML(θ)D2) log x1 + (D1 − n1)

(C1
x1
)
θ
log (C1

x1
)

1 + (C1
x1
)
θ

+ λ̂ML(θ)(D2 − n2)
(C2
x1
)
θ
log (C2

x1
)

1 + (C2
x1
)
θ

−
D1∑
i=1

log x1i − λ̂ML(θ)

D2∑
i=1

log x2i

− 2

D1∑
i=1

( x1
x1i

)θ log ( x1
x1i

)

1 + ( x1
x1i

)θ
− (λ̂ML(θ) + 1)

D2∑
i=1

( x1
x2i

)θ log ( x1
x2i

)

1 + ( x1
x2i

)θ
, (9)

It is noteworthy to emphasis that the solution to Eq.(9) is analytically infeasible. Hence,
the explicit derivation of the maximum likelihood estimator θ̂ML poses a significant chal-
lenge. The required estimates can be obtained using numerical techniques, such as the
Newton-Raphson method. It can be noted from Eq.(8) that the maximum likelihood esti-
mator λ̂ML can be explicitly expressed as a function of the maximum likelihood estimator
of the parameter θ.

2.2.2. EM algorithm

The EM algorithm in this section is employed as a viable alternative method to approx-
imate the ML estimate of the unknown model parameters of the NPT distribution. This
is done by taking into account the incomplete nature of the available sample, which is
characterized as a partially constant-stress ALT under type-I GHCS. The algorithm in
the discussion was initially put forth by Ref.[30] and subsequently received comprehensive
analysis, along with its various adaptations, in the publication authored by Ref.[31]. The
EM algorithm involves a cyclic process that consists of an expectation step (E-step) and
a maximization step (M-step). The E-step involves the computation of an expectation
function for the log-likelihood of the current estimation evaluation, utilizing the given
parameters. On the other hand, the M-step entails the determination of the parameters
that optimize the expected log-likelihood obtained from the E-step. For j = 1, 2, let
Xj = {xj1, . . . , xjDj} be the observed data and Uj = {uj1, . . . , uj(nj−Dj)} be the cen-
sored data under the normal use and accelerated conditions, respectively. For a given
Dj , uj1, . . . , uj(nj−Dj) are not observable. We treated the censored observations as miss-
ing data. Thus, the combination of (Xj ,Uj) forms as Wj = (Xj ,Uj) represents the
complete partially constant-stress ALT failure data set, for which the likelihood function



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 9 of 34

Lc(w; Φ)can be expressed as;

Lc(w; Φ) ∝
2∏

j=1

{ Dj∏
i=1

fj(xji;Φ)

nj−Dj∏
i=1

fj(uji;Φ)

}
. (10)

After disregarding the constants in the above equation, the log-likelihood function denoted
as Lc(w; Φ) can be represented as ℓc(w; Φ);

ℓc(w; Φ) ∝ λn2 log 2 + (n1 + n2) log θ + n2 log λ+ θ(n1 + λn2) log σ

+ (θ − 1)
2∑

j=1

Dj∑
i=1

log xji −
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

log(xθji + σθ)

+ (θ − 1)
2∑

j=1

nj−Dj∑
i=1

log uji −
2∑

j=1

(Sj(λ) + 1)

nj−Dj∑
i=1

log(uθji + σθ). (11)

It is evident that the function ℓc(w; Φ) exhibits a monotonic increase with respected to
σ. Therefore, within the framework of the EM method, the maximum likelihood estimate
(MLE) of σ is given by σ̂EM = x1. By substituting the symbol σ into Eq.(11), we can
derive the profile log-likelihood function for the parameters θ and λ;

ℓc(w; Φ) ∝ λn2 log 2 + (n1 + n2) log θ + n2 log λ+ θ(n1 + λn2) log x1

+ (θ − 1)
2∑

j=1

Dj∑
i=1

log xji −
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

log(xθji + xθ1)

+ (θ − 1)
2∑

j=1

nj−Dj∑
i=1

log uji −
2∑

j=1

(Sj(λ) + 1)

nj−Dj∑
i=1

log(uθji + xθ1). (12)

The EM algorithm consists of two primary steps: The initial phase, known as the expec-
tation step (E-step), is succeeded by the maximization step (M-step), and this iterative
process continues until meeting the specified convergence criteria. During each iteration,
the missing data are imputed with expected values, resulting in the subsequent update of
parameter estimations.

• Expectation Step (E-step):
This phase entails calculating the conditional expectation of the log-likelihood, con-
sidering the incomplete data given the observed data. In order to meet the re-
quirements of the E-stage, it is necessary to calculate the pseudo-log-likelihood
function. This may be derived from the function ℓc(w; Φ) by substituting any
function of uji, denoted as g(uji), with the corresponding conditional expectation

”E

[
g(uji)|uji > Cj

]
”. consequently;

ℓs(θ, λ) = E

[
ℓc(w; θ, λ)|Uj

]
= λn2 log 2 + (n1 + n2) log θ + θ(n1 + λn2) log x1



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 10 of 34

+ (θ − 1)
2∑

j=1

nj−Dj∑
i=1

E

[
log uji

]
−

2∑
j=1

(Sj(λ) + 1)

nj−Dj∑
i=1

E

[
log(uθji + xθ1)

]

+ n2 log λ+ (θ − 1)
2∑

j=1

Dj∑
i=1

log xji −
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

log(xθji + xθ1) (13)

Following that, the expectation step (E-step) requires the calculation of the expected
value of the logarithm of uji, given that uji is greater than Cj . The expression may

be rewritten as E

[
log(uθji + xθ1)

∣∣∣∣uji > Cj

]
, where uji and x1 are variables, θ is a

parameter. The conditional probability function of the censored data, given the
observed data, can be computed according to Ref.[32].

fUj |Xj
(uji|xji; Φ) =

f(uji,Φ)

1− F (xji,Φ)

=
θSj(λ)(C

θ
j + σθ)Sj(λ)uθ−1

ji

(uθji + σθ)Sj(λ)+1
, uji > Cj , i = 1, . . . , Dj , j = 1, 2.

(14)

Hence, the required expected values of a NPT distribution from the left at Cj are,
respectively, given by;

E

[
log(uji)|uji > Cj

]
= θSj(λ)(C

θ
j + xθ1)

Sj(λ)

∫ ∞

Cj

yθ−1 log(y)

(yθ + xθ1)
Sj(λ)+1

dy

= log(Cj) + (Cθ
j + xθ1)

Sj(λ)

∫ ∞

Cj

1

y(yθ + xθ1)
Sj(λ)

dy

= ℑ1(Cj , θ, λ). (15)

and

E

[
log(uθji + xθ1)|uji > Cj

]
= θSj(λ)(C

θ
j + xθ1)

Sj(λ)

∫ ∞

Cj

yθ−1 log(yθ + xθ1)

(yθ + xθ1)
Sj(λ)+1

dy

= log(Cθ
j + xθ1) +

1

Sj(λ)

= ℑ2(Cj , θ, λ). (16)

• Maximization Step (M-step):
This step involves the maximization of the pseudo log-likelihood function Eq.(13)
with respect to (θ, λ). In this regard suppose (θ(s), λ(s)) denotes the sth stage esti-
mate of (θ, λ) then the next stage updated estimate (θ(s+1), λ(s+1)) is obtained by
maximizing the following function;

Q(θ, λ) = λn2 log 2 + (n1 + n2) log θ + θ(n1 + λn2) log x1 + (θ − 1)

2∑
j=1

Dj∑
i=1

log xji



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 11 of 34

−
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

log(xθji + xθ1) + (θ − 1)
2∑

j=1

(nj −Dj)ℑ1(Cj , θ(s), λ(s))

+ n2 log λ−
2∑

j=1

(Sj(λ) + 1)(nj −Dj)ℑ2(Cj , θ(s), λ(s)). (17)

By taking the derivatives of Eq.(17) with respect to θ and λ, respectively, and
equating them to zero, as follows;

∂Q(θ, λ)

∂θ
=

n1 + n2

θ
+ (n1 + λn2) log x(1) +

2∑
j=1

(nj −Dj)ℑ1(Cj , θ(s), λ(s))

+

2∑
j=1

Dj∑
i=1

log xji −
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

xθji log xji + xθ1 log x1

xθji + xθ1
= 0, (18)

and

∂Q(θ,λ)
∂λ = n2 log 2 +

n2
λ + θ n2 log x(1) −

∑D2
i=1 log(x

θ
2i + xθ1)− (n2 −D2)ℑ2(C2, θ(s), λ(s)) = 0.

(19)

For this purpose, we first find θ(s+1) by using the method of fixed point iteration as
Ref. [33]. Where we have;

h(θ) = n1+n2

−(n1+λ̂(θ)n2) log x(1)−
∑2

j=1

∑Dj
i=1 log xji−

∑2
j=1(nj−Dj)ℑ1(Cj ,θ(s),λ(s))+2A1+(λ̂(θ)+1)A2

,

(20)
with

Aj =

{ Dj∑
i=1

xθji log(xji) + xθ1 log(x1)

xθji + xθ1

}
and

λ̂(θ) =
n2

−n2 log(2)− θn2 log(x(1)) +
∑D2

i=1 log(x
θ
2i + xθ1) + (n2 −D2)ℑ2(C2, θ(s), λ(s))

.

(21)

Here, θ(s) and λ(s) are the estimate of θ and λ in the sth step, respectively. Finally

after finding θ(s+1) , the estimate λ(s+1) is derived as λ(s+1) = λ̂(θ(s+1)).

The MLEs of (θ, λ) can be obtained by repeating the E-step and M-step until convergence
is achieved. Now, an iterative process can be employed to obtain the necessary maximum
likelihood estimates of (θ, λ). This iterative procedure continues until.

| θ(s+1) − θ(s) |+ | λ(s+1) − λ(s) | < ϵ

for a predetermined small value of ϵ and some s. This approach converges to the local
maximum likelihood as the log-likelihood increases with each iteration. In the expectation-
maximization technique, we initialize the parameters using their maximum likelihood es-
timates based on the entire sample.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 12 of 34

2.2.3. Bayes estimation

Unlike traditional statistics, Bayesian estimation incorporates prior information about life
parameters. This approach combines the data provided with prior probabilities to infer
the parameters of interest, making the inference process more objective and reasonable.
A loss function is used to evaluate the difference between an estimated value and the true
value of the parameter. In this section, we focus on objective Bayesian estimation of the
unknown parameters θ, σ, and λ with respect to squared error (SE) and linear exponential
(LINEX) loss functions. For the selection of prior distribution of unknown parameters is
discussed and its simulation results are satisfactory in Ref.[8]. We adopt the same prior
distribution as this Ref.[8]. The joint prior probability density function (PDF) of θ and σ
is considered to be given by;

π∗
1(θ, σ) ∝ θaσθb−1c−θ, θ > 0, 0 < σ < d,

for θ and σ, where a, b, c, d are positive constants and db < c. This joint prior has been
used in [34, 35] for Bayesian inference of the two-parameter Pareto distribution.
Such a prior specifies π∗

1(θ) as a gamma distribution Ga(a, log c− b log d) and π∗
1(σ|θ) as

a power function distribution of the form;

π∗
1(σ|θ)) = b θ σbθ−1 d−bθ, 0 < σ < d.

By letting a = −1, b = 0, c = 1 and d −→ ∞ , this prior reduces to the non informative
prior;

π∗
1(σ, θ) ∝ 1

θσ
, θ, σ > 0

On the other hand, the prior for the acceleration factor λ is assumed to be the non-
informative prior, i.e.

π∗
2(λ) ∝ λ−1 , λ > 1.

Therefore, the joint prior PDF of θ, σ and λ can be written as,

π∗(θ, σ, λ) ∝ θaσθb−1c−θλ−1, θ > 0, 0 < σ < d, λ > 1. (22)

In combining the prior information from Eq.(22) with the likelihood function presented in
Eq.(4), the joint posterior distribution can be represented as follows;

π(θ, σ, λ|X) = A−1 L(θ, σ, λ|X) π∗(θ, σ, λ)

= A−1 2λn2 θ(D1+D2+a) c−θ λD2−1 σθ(D1+λD2+b)−1

[
1 + (

C1

σ
)
θ ](D1−n1)

×
[
1 + (

C2

σ
)
θ ]λ (D2−n2) D1∏

i=1

x
−(θ+1)
1i[

1 + ( σ
x1i

)θ
]2 D2∏

i=1

x
−(λθ+1)
2i[

1 + ( σ
x2i

)θ
]λ+1



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 13 of 34

= A−1 2λn2 θ(D1+D2+a+1)−1 c−θ λD2−1 σθ(D1+λD2+b)−1

[
1 + (

C1

σ
)
θ ](D1−n1)

×
[
1 + (

C2

σ
)
θ ]λ(D2−n2)

e−(θ+1)
∑D1

i=1 log x1i e
−2
∑D1

i=1 log(1+( σ
x1i

)θ)

× e−(λθ+1)
∑D2

i=1 log x2i e
−(λ+1)

∑D2
i=1 log(1+( σ

x2i
)θ)

. (23)

Where A is a normalizing constant, provided by;

A =

∫ ∞

1

∫ d

0

∫ ∞

0
L(θ, σ, λ|X) π∗(θ, σ, λ) dθdσdλ.

It is clear that Eq.(23) is analytically tricky. Furthermore, the Bayesian estimation of a
function with θ, σ and λ is also intractable because it is related to a ratio of three integrals.
For solving the corresponding ratio of three integrals, some approximate approaches have
been presented in the literature. Among them, the MH algorithm is a simulation method
with wide applications in sampling from posterior density functions.
In this article, we use the MH algorithm to derive approximate explicit forms for the
Bayesian estimates. In Bayesian statistics, the selection of the loss function is a funda-
mental step. There are many symmetric loss functions, among which is SE loss function
is well known for its good mathematical properties. Let Υ∗ be a Bayesian estimate of Υ.
The form of SE loss function is;

LSE(Υ
∗,Υ) = ( Υ−Υ∗ )2

Now, we compute the Bayes estimate of Υ(θ, σ, λ) under the SE loss function.

Υ∗(θ, σ, λ)SE = E[Υ(θ, σ, λ) |X]

=

∫ ∞

1

∫ d

0

∫ ∞

0
Υ(θ, σ, λ)π(θ, σ, λ |X) dθdσdλ. (24)

However, in many practical situations, overestimation and underestimation result in dif-
ferent losses, and the consequence is likely to be quite serious if one uses symmetric loss
function indiscriminately. In the literature, many different asymmetric loss functions were
used. Among them, LINEX is dominant, and this loss function can be expressed as;

LLL(Υ
∗,Υ) = eh

∗(Υ∗−Υ) − h∗(Υ∗ −Υ)− 1, h∗ ̸= 0.

The constant h∗ represented the weight of error on different decisions. Under the above
loss function, the Bayesian estimate of the function Υ(θ, σ, λ) can be calculated by;

Υ∗(θ, σ, λ)LL = − 1

h∗
ln[E(e−h∗Υ(θ,σ,λ) |X)]

= − 1

h∗
ln

[∫ ∞

1

∫ D∗

0

∫ ∞

0
e−h∗Υ(θ,σ,λ) π(θ, σ, λ |X) dθdσdλ

]
. (25)



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 14 of 34

For θ, σ and λ , the full conditionals can be expressed as;

π1(θ|σ, λ) ∝ θ(D1+D2+a)c−θσθ(D1+λD2+b)

[
1 + (

C1

σ
)
θ](D1−n1)[

1 + (
C2

σ
)
θ]λ(D2−n2)

× e
−θ
∑D1

i=1 log x1i−2
∑d1

i=1 log(1+( σ
x1i

)θ)−λθ
∑D2

i=1 log x2i−(λ+1)
∑D2

i=1 log(1+( σ
x2i

)θ)
(26)

π2(σ|θ, λ) ∝ σθ(D1+λD2+b−1)

[
1 + (

C1

σ
)
θ](D1−n1)[

1 + (
C2

σ
)
θ]λ(D2−n2)

× e
−2
∑D1

i=1 log(1+( σ
x1i

)θ)−(λ+1)
∑D2

i=1 log(1+( σ
x2i

)θ)
(27)

π3(λ|θ, σ) ∝ λD2−1 σθλD2

[
1 + (

C2

σ
)
θ]λ(D2−n2)

e
λn2 log 2−λθ

∑D2
i=1 log x2i−λ

∑D2
i=1 log(1+( σ

x2i
)θ)

(28)

The conditional posterior functions of the parameters θ, σ and λ in Eqs.(26–28) do not
present standard form. Therefore, we generate random samples from these distributions
using the Metropolis–Hastings (M–H) algorithm with proposal distribution.
Since Gibbs sampling is not a clear-cut alternative and the conditional posteriors of θ, σ
and λ in Eqs.(26–28) do not give conventional forms, the application of the M–H sampler
is necessary for the implementation of the Markov chain Monte Carlo (MCMC) approach.
A detailed discussion about MCMC and M-H algorithm can be found in Refs.[36],[37],
and [38]. The procedures outlined in Algo.(1) are utilized to derive Bayesian estimation
for parameters θ, σ and λ .



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 15 of 34

Algorithm 1. MH Sampling for Bayesian estimation.

Step 1: Set l = 1 and the initial value of the parameters as (θ(0), σ(0), λ(0)) =
(θ̂ML, σ̂ML, λ̂ML).

Step 2: Generate θ(l), σ(l) and λ(l) such that
θ(l) ∼ N(θ(l−1),Var(θ̂ML)), σ(l) ∼ N(σ(l−1),Var(σ̂ML)), and λ(l) ∼
N(λ(l−1),Var(λ̂ML)).

Step 3: Obtain; 

Pθ = min

(
π1(θ

l|σl−1, λl−1)

π1(θl−1|σl−1, λl−1)
, 1

)
,

Pσ = min

(
π2(σ

l|θl, λl−1)

π2(σl−1|θl, λl−1)
, 1

)
,

Pλ = min

(
π3(λ

l|θl, σl)

π3(λl−1|θl, σl)
, 1

)
.

Step 4: From a Uniform U(0, 1) distribution, we produce u1, u2 and u3.

Step 5: If u1 ≤ Pθ, set θ
(l) = θ(l), else θ(l) = θ(l−1). Similarly, If u2 ≤ Pσ, set σ

(l) = σ(l), else
σ(l) = σ(l−1) and If u3 ≤ Pλ, set λ

(l) = λ(l), else λ(l) = λ(l−1).

Step 5: Set l = l + 1. Then repeat steps 2 − 5 for B times to generate (θ(l), σ(l), λ(l)), for
l = 1, . . . , B.

Step 6: After discarding the first M number of burn-in samples, the remaining B−M samples
are used to obtain the Bayesian estimates of θ, σ, and λ, where the Bayes estimate
of Φ = Φ(θ, σ, λ) under the SE and LINEX loss functions can now be computed as;

Φ̂SE =
1

B −M

B∑
i=M+1

Φ(θ(i), σ(i), λ(i)) , Φ̂LL = − 1

h∗
log

[
1

B −M

B∑
i=M+1

e−h∗Φ(θ(i),σ(i),λ(i))

]
.

2.3. Interval estimation

In this section, the asymptotic confidence interval (CI), normal bootstrap confidence
interval (N-boot) and Bayesian CI (credibility interval) have been studied for parameters
θ, σ and λ.

2.3.1. Asymptotic Confidence Interval

Rather of acquiring specific estimates for unknown parameters, there may be a desire to
obtain a range of values that could potentially encompass these parameters with a desig-



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 16 of 34

nated level of probability. The aforementioned ranges are sometimes referred to as interval
estimations.
In this context, we employ the asymptotic properties of the maximum likelihood es-
timators (MLEs) to construct the asymptotic confidence intervals (ACIs) for the un-
known parameters Φ = (θ, σ, λ)T . Using the large sample theory reveals that the MLEs,
Φ̂ML = (θ̂ML, σ̂ML, λ̂ML) follow an asymptotic distribution that may be approximated
by a normal distribution. This normal distribution has a mean of Φ and a variance-
covariance matrix of I−1(Φ). In this study, we employ the asymptotic variance-covariance
matrix (AVCM) denoted as I−1(Φ̂) to estimate I−1(Φ). The estimation is achieved by
inverting the observed Fisher information matrix. In this particular instance, the AVCM
assumes the structure.

I−1(Φ̂) = −


∂2ℓ(Φ)
∂θ2

∂2ℓ(Φ)
∂θ∂σ

∂2ℓ(Φ)
∂θ∂λ

∂2ℓ(Φ)
∂σ2

∂2ℓ(Φ)
∂σ∂λ
∂2ℓ(Φ)
∂λ2


(Φ̂)

=

ϖ̂2
11 ϖ̂12 ϖ̂13

ϖ̂2
22 ϖ̂23

ϖ̂2
33

 , (29)

The notation Φ̂ indicates that the derivatives are being evaluated at the estimated values
of θ, σ, and λ. After a little calculation, one has;

∂2ℓ(Φ)

∂θ2
= −D1 +D2

θ2
+ (D1 − n1)

(
C1
σ

)θ
(log(C1

σ ))2(1 + (C1
σ )θ)−

(
(C1
σ )θ log(C1

σ )
)2

(1 + (C1
σ )θ)2

+ λ(D2 − n2)

(
C2
σ

)θ
(log(C2

σ ))2(1 + (C2
σ )θ)−

(
(C2
σ )θ log(C2

σ )
)2

(1 + (C2
σ )θ)2

−
2∑

j=1

(Sj(λ) + 1)

Dj∑
i=1

(
σ
xji

)θ
(log( σ

xji
))2(1 + ( σ

xji
)θ)−

(
( σ
xji

)θ log( σ
xji

)
)2

(1 + ( σ
xji

)θ)2
. (30)

∂2ℓ(Φ)

∂σ2
= (D1 − n1)

θ
σ2 (

C1
σ )θ(θ(C1

σ )θ + (C1
σ )θ + 1)

(1 + (C1
σ )θ)2

+ λ(D2 − n2)
θ
σ2 (

C2
σ )θ(θ(C2

σ )θ + (C2
σ )θ + 1)

(1 + (C2
σ )θ)2

− θ(D1 + λD2)

σ2
−

2∑
j=1

(Sj(λ) + 1)

Dj∑
i=1

θ
x2
ji
( σ
xji

)(θ−2)(θ − ( σ
xji

)θ − 1)

(1 + ( σ
xji

)θ)2
, (31)

∂2ℓ(Φ)

∂λ2
= −D2

λ2
, (32)

∂2ℓ(Φ)

∂θ∂σ
=

∂2ℓ(Φ)

∂σ∂θ
=

n1 + λn2

σ
+

(D1 − n1)σ
θ−1(1 + θ log σ)

σθ + Cθ
1

+
(D2 − n2)σ

θ−1λ(1 + θ log σ)

σθ + Cθ
2



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 17 of 34

+
(n1 −D1)θσ

θ−1(σθ log σ + Cθ
1 logC1)

(σθ + Cθ
1)

2
+

(n2 −D2)θσ
θ−1λ(σθ log σ + Cθ

2 logC2)

(σθ + Cθ
2)

2

− (1 + λ)

D2∑
i=1

(
σθ−1 + θσθ−1 log σ

σθ + xθi
− θσθ−1(σθ log σ + xθi log xi)

(σθ + xθi )
2

)

− 2

D1∑
i=1

(
σθ−1 + θσθ−1 log σ

σθ + xθi
− θσθ−1(σθ log θ + xθi log xi)

(σθ + xθi )
2

)
, (33)

∂2ℓ(Φ)

∂θ∂λ
=

∂2ℓ(Φ)

∂λ∂θ
=

(D2 − n2)(σ
θ log σ + Cθ

2 logD2)

σθ + Cθ
2

−
D2∑
i=1

σθ log σ + xαi log xi

σθ + xθi
, (34)

and

∂2ℓ(Φ)

∂σ∂λ
=

∂2ℓ(Φ)

∂λ∂σ
=

(D2 − n2)θσ
θ−1

σθ + Cθ
2

−
D2∑
i=1

θσθ−1

σθ + xθi
. (35)

Using the asymptotic distribution of the maximum likelihood estimators (MLE), specifi-
cally the multivariate normal distribution with a mean vector of zero and a variance−covariance
matrix as described in Eq.(29), it can be observed that the expression

√
n(Φ̂−Φ) follows

this distribution. Hence, for any arbitrary value of 0 < υ < 1, the 100(1−υ)% Asymptotic
Confidence Intervals (ACIs) for the unknown parameters can be expressed in the following
manner; 

θ̂ ± Zυ/2ϖ̂11,

σ̂ ± Zυ/2ϖ̂22,

λ̂± Zυ/2ϖ̂33.

(36)

The conventional tabular normal values Zυ/2 with tail (υ/2) are denoted correspond-
ingly.

2.3.2. Bootstrap Methods

The bootstrap method is a widely recognized resampling technique used in statistics to
estimate the sampling distribution of a statistic. It achieves this by repeatedly sampling
with replacement from the observed data. The main objective of this method is to offer
an empirical approximation of the sampling distribution of a statistic, which proves use-
ful for making inferences and constructing confidence intervals. The bootstrap method,
originally introduced by Ref.[39], has been rigorously examined and extensively debated
in the academic literature. In this subsection , we use the parametric bootstrap method
to construct normal bootstrap (N-boot) confidence interval (CI) for the unknown model
parameter (θ, σ, λ) , see Ref.[40]. The procedures outlined in Algo.(2) are utilized to derive
100(1− υ)% Normal bootstrap confidence intervals for these parameters.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 18 of 34

Algorithm 2. Normal bootstrap CI for (θ, σ, λ).

Step 1: Based on the original type-I GHCS Xj = {Xj1, . . . , Xjdj}, j = 1, 2, and X =

{X1, X2}, obtain the estimates of (θ, σ, λ), say (θ̂ML, σ̂ML, λ̂ML) .

Step 2: Employ the censoring plan nj , lj ,mj , T and (θ̂ML, σ̂ML, λ̂ML) to generate a type-I
GHCS bootstrap sample X∗

j = {X∗
j1, . . . , X

∗
jdj}, j = 1, 2.

Step 3: From the ordered observations obtained in Step 2, the bootstrap sample estimates
of (θ, σ, λ) are computed, namely (θ̂∗, σ̂∗, λ̂∗).

Step 4: To get B∗ bootstrap samples, repeat Steps 2-3 B several times.

Step 5: Arrange all (θ̂∗, σ̂∗, λ̂∗) in ascending order and denote,(
θ̂∗[1], θ̂∗[2], . . . θ̂∗[B

∗]

)
,

(
σ̂∗[1], σ̂∗[2], . . . σ̂∗[B∗]

)
, and

(
λ̂∗[1], λ̂∗[2], . . . λ̂∗[B∗]

)
Step 6: To rigorously construct the Normal bootstrap (N-boot) confidence intervals for

κj , j = 1, 2, 3, where κ1 = θ, κ2 = σ and κ3 = λ by the following, formulae:(
2κ̂∗jML − κ∗j − Z1−υ

2

√
S(κ̂∗j ), 2κ̂∗jML + κ∗j − Z1−υ

2

√
S(κ̂∗j )

)
where ZP is P th the quantile of the standard normal distribution,

κ∗j = (B∗)−1
B∗∑
i=1

κ̂
∗[i]
j , and S(κ̂∗j ) = (B∗ − 1)−1

B∗∑
i=1

(
κ̂
∗[i]
j − κ∗j

)2

2.3.3. Highest Posterior Density Credible Interval

In subs. (2.2.3), for Bayesian point estimation, B −M samples of (θ, σ, λ) were;(
θ̂(1), θ̂(2), . . . , θ̂(B−M)

)
,

(
σ̂(1), σ̂(2), . . . , σ̂(B−M)

)
, and

(
λ̂(1), λ̂(2), . . . , λ̂(B−M)

)
.

For arbitrary 0 < υ < 1, a 100(1− υ)% confidence interval for parameter (θ, σ, λ) can be
composed in the following form,(
θ̂[(B−M)υ

2
], θ̂[(B−M)( 1−υ

2
)]

)
,

(
σ̂[(B−M)υ

2
], σ̂[(B−M)( 1−υ

2
)]

)
,and

(
λ̂[(B−M)υ

2
], λ̂[(B−M)( 1−υ

2
)]

)
Where [(B −M)υ2 ] represents the maximum integer less than (B −M)υ2 ; similarly, [(B −
M)(1−υ

2 )] represents the maximum integer less than [(B − M)(1−υ
2 )]. Repeat the above

steps M∗ times to obtain M∗ interval estimates with the form above, from which the in-
terval with the smallest length is selected as the highest posterior density credible interval.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 19 of 34

3. Numerical Application

3.1. Monte Carlo Simulations

To investigate the effectiveness of frequentest and Bayesian estimations produced in
the preceding sections, Monte Carlo simulations were carried out. We compare the ef-
ficacy of different estimators and confidence intervals of the NPT parameters (θ, σ) and
the acceleration factor λ, based on the true parameter values of (θ, σ, λ), namely Set 1:
(0.5, 0.5, 1.3) and Set 2: (1.0, 1.0, 2.0). We generate 1000 type-I GHCS based on various
choices of (nj , (lj ,mj), T ), j = 1, 2. The sample size (nj , (lj ,mj), T ) is fixed from the data
under type-I GHCS and is set to (n1, n2) = (30, 40), (50, 60), (70, 80) with three sets of
fixed numbers (l1,m1) = (15, 25), (35, 45),
(55, 65), (l2,m2) = (20, 35), (40, 55), (60, 75) and T = 5.0, 7.0 presented respectively for
each size. The process of generating the type-I GHCS data is shown in Algo.(3).

Algorithm 3. Steps for generating type-I GHCS datasets.

Step 1: Create nj independent variables ϵj = (ϵj1, ϵj2, . . . , ϵjnj ) from Uniform(0, 1) distribu-
tion, j = 1, 2.

Step 2: For given lj ,mj , T , set Dj , j = 1, 2 as;

Dj =


lj if T < lj

d∗j if lj ≤ T < mj

mj if T ≥ mj

Step 3: Record the type-I GHCS data as Xj = (Xj1, Xj2, . . . , XjDj ), j = 1, 2, are derived
by the inverse function method;X1 = σ

(
1+ϵ1
1−ϵ1

)1/θ

X2 = σ
(

2
(1−ϵ2)1/λ

− 1
)1/θ

Once the samples have been generated, using R (4.1.3) software with the ”maxLik”
package introduced by Ref.[41], we obtain the MLEs of parameters (θ, σ, λ) based on both
NR and EM algorithm methods. We set the real values of the parameters (θ, σ, λ) as
the initial guesses of EM algorithm and the convergence is assumed when the absolute
differences between the successive estimates are less than 10−5. We computed the 100(1−
υ)% confidence intervals (CIs) for (θ, σ, λ) based on the asymptotic normal distribution
of the MLEs. We have also computed normal bootstrap (N-boot) confidence intervals as
well. We have reported the ALs and the CPs in each case. What’s more, for the numbers
of bootstrap re-sampling, we set it as B∗ = 1000. For comparing the performance of the



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 20 of 34

Bayesian estimates, the associated values of the hyper-parameters (a, b, c, d) were taken
as (0.02, 1.0, 4.0, 2.0) and, (0.05, 2.0, 6.0, 3.0) for the given Sets 1 and 2, respectively.
The M-H within the Gibbs algorithm given in Subs. (2.2.3) is used to obtain Bayesian
estimates. For this algorithm, we took a normal proposal density to generate θ and λ
from Eqs.(26) and (28), a power function distribution to generate σ from Eq.(27), using
the ‘coda’ package proposed by Ref.[42], and we take into account the MLEs as initial
guess values. We generate B = 11, 000 MCMC samples and discard the first M = 1000
values as burn-in period as described in Subs. (2.2.3). The last 10, 000 Markov chains
are used to establish the empirical distribution to estimate the posterior distribution of
(θ, σ, λ). Consequently, the Bayes estimates under the specified squared error (SE) and
linear exponential (LINEX) loss functions must be rigorously derived from the empirical
posterior distribution of (θ, σ, λ). Furthermore, the 100(1 − υ)% credible intervals for
(θ, σ, λ) must be precisely computed by extracting the two symmetric quantiles from the
same empirical posterior distribution of (θ, σ, λ), respectively. For evaluating the Bayes
estimators under LINEX loss function we take h∗ = (−0.8, 0.8). Also, ALs and the CPs
of Bayesian credible intervals are calculated. The interval estimates are computed using
the nominal υ = 0.05 significance level in each case and the performances of different
estimation methods are compared based on 1000 replications. A comparison of several
point estimations of Φj , j = 1, 2, 3, where Φ1 = θ,Φ2 = σ, and Φ3 = λ, is then made based
on two different criteria, namely average values (Avg) and mean-square errors (MSEs), by
the following formulae;

Avg = (N∗)−1
N∗∑
i=1

Φ̂
(i)
j and MSE = (N∗)−1

N∗∑
i=1

(
Φ̂
(i)
j − Φj

)2
.

Respectively, where N∗ is the number of generated sequence data, and Φ̂
(i)
j denotes the

calculated estimate at the ith sample of Φj , j = 1, 2, 3. Additionally, the evaluation of
various interval estimates of Φj , j = 1, 2, 3 is determined by two other standards, namely
the average confidence lengths (ACLs) and coverage percentages (CPs), by the following
formulae;

ACL(1−υ)%(Φp) = (N∗)−1
∑N∗

i=1

(
U∗
Φ̂

(i)
p

− L∗
Φ̂

(i)
p

)
and CP(1−υ)%(Φp) = (N∗)−1

∑N∗

i=11(
L∗
Φ̂
(i)
p

; U∗
Φ̂
(i)
p

)(Φp).

Respectively, where 1(·) is the indicator function, and (L∗
(·) , U

∗
(·)) are the lower and upper

bounds of asymptotic (or credible) interval estimates. The results of the Avg, MSEs, ALs
and CPs for the Bayes estimates for (θ, σ, λ) are strictly presented in Tabs.(1)-(4).
On the basis of the results reported in Tabs.(1)-(4), some points can be drawn which are
stated as follows;

• The proposed point (or interval) estimates of θ, σ and λ have shown good perfor-
mance based on given true parameter set.

• In every instance, as expected, the estimation results are satisfactory based on the
Avg. The MSEs of all estimates decrease as the sample size increases, confirming
the consistency of each estimation method.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 21 of 34

• The MLEs of θ, σ and λ using the EM algorithm have smaller MSEs than the MLEs
using the NR algorithm. Hence, the MLEs via the EM algorithm perform better
than those obtained by the NR method.

• In terms of Avg and MSEs, Bayes estimation using MCMC performs better than
the other approaches (ML, N-boot).

• Due to having the smallest MSE and narrowest width, MCMC CRIs are, overall,
the most satisfactory.

• The improvement of CPs and CIs is signifcant with increased total and efective
sample sizes.

• The estimates produced by the ML, bootstrap, and Bayesian approaches are highly
similar and have high CPs (around 0.95).

• In most simulations, the Bayes estimates outperform the MLEs for the estimation
of θ, σ, and λ. So, in general, we would recommend using the Bayes estimate of
the unknown parameters of NPT distribution based on type-I GHCS under partially
constant-stress ALT.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 22 of 34

blue

Table 1: Avg and MSEs of for the parameters with (θ, σ, λ) = (0.5, 0.5, 1.3).

(n1, n2)
(l1,m1)
(l2,m2)

T ParCriteria
MLE Bayesian

NR EM SELF
LLF

h∗ = −0.8 h∗ = 0.8

(30, 40)

5.0

θ̂
Avg 0.392 0.432 0.476 0.482 0.472
MSE 0.015 0.006 0.008 0.009 0.008

σ̂
Avg 0.521 0.521 0.514 0.517 0.512
MSE 0.001 0.001 0.003 0.004 0.003

λ̂
Avg 1.135 1.394 1.432 1.567 1.342

(15, 25) MSE 0.077 0.015 0.012 0.021 0.076
(20, 35)

10.0

θ̂
Avg 0.405 0.445 0.502 0.506 0.497
MSE 0.013 0.027 0.011 0.012 0.010

σ̂
Avg 0.528 0.528 0.595 0.598 0.592
MSE 0.001 0.001 0.013 0.024 0.014

λ̂
Avg 1.173 1.406 1.506 1.632 1.418
MSE 0.077 0.074 0.137 0.237 0.117

(50, 60)

5.0

θ̂
Avg 0.423 0.521 0.518 0.520 0.516
MSE 0.008 0.006 0.007 0.008 0.007

σ̂
Avg 0.515 0.515 0.547 0.548 0.546
MSE 0.001 0.001 0.006 0.006 0.005

λ̂
Avg 1.132 1.232 1.392 1.439 1.349

(35, 45) MSE 0.064 0.034 0.090 0.113 0.075
(40, 55)

10.0

θ̂
Avg 0.421 0.521 0.519 0.521 0.517
MSE 0.008 0.004 0.006 0.007 0.006

σ̂
Avg 0.516 0.516 0.569 0.571 0.568
MSE 0.001 0.001 0.007 0.008 0.007

λ̂
Avg 1.141 1.234 1.392 1.439 1.351
MSE 0.066 0.045 0.093 0.115 0.078

(70, 80)

5.0

θ̂
Avg 0.420 0.532 0.499 0.500 0.497
MSE 0.008 0.007 0.003 0.004 0.004

σ̂
Avg 0.511 0.511 0.520 0.529 0.521
MSE 0.002 0.002 0.003 0.004 0.003

λ̂
Avg 1.181 1.380 1.405 1.436 1.377

(55, 65) MSE 0.051 0.045 0.094 0.109 0.082
(60, 75)

10.0

θ̂
Avg 0.431 0.532 0.514 0.516 0.513
MSE 0.006 0.007 0.005 0.004 0.004

σ̂
Avg 0.509 0.509 0.520 0.521 0.519
MSE 0.001 0.001 0.002 0.003 0.002

λ̂
Avg 1.148 1.345 1.351 1.378 1.325
MSE 0.047 0.048 0.048 0.055 0.043



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 23 of 34

blue

Table 2: Avg and MSEs of for the parameters with (θ, σ, λ) = (1.0, 1.0, 2.0).

(n1, n2)
(l1,m1)
(l2,m2)

T ParCriteria
MLE Bayesian

NR EM SELF
LLF

h∗ = −0.8 h∗ = 0.8

(30, 40)

7.0

θ̂
Avg 1.098 1.072 0.934 0.936 0.928
MSE 0.086 0.079 0.066 0.062 0.073

σ̂
Avg 1.319 1.319 1.291 1.423 1.397
MSE 0.097 0.097 0.076 0.078 0.066

λ̂
Avg 1.519 1.572 1.373 1.468 1.386

(15, 25) MSE 0.244 0.278 0.052 0.078 0.053
(20, 35)

10.0

θ̂
Avg 1.154 1.131 0.879 0.913 0.869
MSE 0.083 0.080 0.044 0.052 0.043

σ̂
Avg 1.134 1.134 1.096 1.125 0.983
MSE 0.017 0.017 0.050 0.066 0.050

λ̂
Avg 1.370 1.412 1.531 1.612 1.435
MSE 0.171 0.188 0.103 0.126 0.105

(50, 60)

7.0

θ̂
Avg 1.014 0.983 0.869 0.923 0.895
MSE 0.033 0.034 0.030 0.036 0.029

σ̂
Avg 1.082 1.082 0.884 0.896 0.865
MSE 0.085 0.085 0.074 0.082 0.061

λ̂
Avg 1.714 1.787 1.731 1.752 1.693

(35, 45) MSE 0.203 0.269 0.090 0.099 0.089
(40, 55)

10.0

θ̂
Avg 1.023 1.005 0.867 0.902 0.849
MSE 0.030 0.031 0.023 0.032 0.019

σ̂
Avg 1.072 1.072 0.857 0.879 0.847
MSE 0.009 0.009 0.004 0.025 0.001

λ̂
Avg 1.543 1.581 1.725 1.865 1.658
MSE 0.125 0.141 0.122 0.156 0.069

(70, 80)

7.0

θ̂
Avg 1.021 0.991 0.962 0.976 0.936
MSE 0.017 0.018 0.015 0.026 0.011

σ̂
Avg 1.082 1.082 0.913 0.963 0.896
MSE 0.012 0.012 0.003 0.019 0.001

λ̂
Avg 1.485 1.537 1.822 1.836 1.756

(55, 65) MSE 0.037 0.042 0.022 0.056 0.015
(60, 75)

10.0

θ̂
Avg 0.944 0.928 0.918 0.968 0.897
MSE 0.020 0.023 0.018 0.032 0.009

σ̂
Avg 1.096 1.096 0.919 0.936 0.895
MSE 0.010 0.010 0.009 0.023 0.002

λ̂
Avg 1.597 1.634 1.716 1.768 1.659
MSE 0.142 0.169 0.139 0.189 0.106



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 24 of 34

blue

Table 3: ALs and CPs for the parameters with (θ, σ, λ) = (0.5, 0.5, 1.3).

(n1, n2)
(l1,m1) T Par Criteria MLE N-boot Bayesian
(l2,m2)

(30, 40)

5.0

θ̂
AL 0.295 0.529 0.324

95% CP 0.852 0.872 0.868

σ̂
AL 0.406 0.599 0.288

95% CP 0.897 0.886 0.854

λ̂
AL 1.933 1.789 1.324

(15, 25) 95% CP 0.912 93.8 0.867

(20, 35)

10.0

θ̂
AL 0.201 0.517 0.311

95% CP 0.789 0.897 0.967

σ̂
AL 0.451 0.599 0.315

95% CP 0.852 0.916 0.921

λ̂
AL 1.463 1.765 1.311

95% CP 0.890 0.923 0.873

(50, 60)

5.0

θ̂
AL 0.221 0.567 0.293

95% CP 0.893 0.936 0.827

σ̂
AL 0.348 0.577 0.173

95% CP 0.917 0.944 0.878

λ̂
AL 1.597 1.679 1.275

(35, 45) 95% CP 0.908 0.962 0.925

(40, 55)

10.0

θ̂
AL 0.167 0.568 0.292

95% CP 0.908 0.962 0.925

σ̂
AL 0.359 0.577 0.186

95% CP 0.923 0.938 0.932

λ̂
AL 1.271 1.672 1.266

95% CP 0.946 0.950 0.889

(70, 80)

5.0

θ̂
AL 0.150 0.608 0.221

95% CP 0.902 0.917 0.871

σ̂
AL 0.311 0.567 0.112

95% CP 0.918 0.915 0.885

λ̂
AL 1.245 1.563 1.035

(55, 65) 95% CP 0.902 0.956 0.938

(60, 75)

10.0

θ̂
AL 0.140 0.612 0.229

95% CP 0.919 0.923 0.909

σ̂
AL 0.310 0.566 0.118

95% CP 0.927 0.936 0.910

λ̂
AL 1.112 1.556 0.996

95% CP 0.884 0.922 0.918



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 25 of 34

blue

Table 4: ALs and CPs for the parameters with (θ, σ, λ) = (1.0, 1.0, 2.0).

(n1, n2)
(l1,m1) T Par Criteria MLE N-boot Bayesian
(l2,m2)

(30, 40)

5.0

θ̂
AL 1.164 0.534 0.775

95% CP 0.975 0.923 0.965

σ̂
AL 1.927 1.324 1.275

95% CP 0.960 0.882 97.2

λ̂
AL 1.444 1.653 1.022

(15, 25) 95% CP 0.935 0.893 0.965

(20, 35)

10.0

θ̂
AL 0.978 0.563 0.797

95% CP 0.938 0.913 0.923

σ̂
AL 1.083 0.698 0.828

95% CP 0.966 0.946 0.977

λ̂
AL 1.710 1.653 1.030

95% CP 0.975 0.943 96.3

(50, 65)

5.0

θ̂
AL 0.942 0.698 0.707

95% CP 0.965 0.930 0.973

σ̂
AL 1.363 0.963 0.673

95% CP 0.956 0.933 0.968

λ̂
AL 1.457 1.036 1.008

(35, 45) 95% CP 0.954 0.936 0.976

(40, 55)

10.0

θ̂
AL 0.721 0.563 0.633

95% CP 0.965 0.963 0.955

σ̂
AL 0.947 0.678 0.585

95% CP 0.955 0.906 0.965

λ̂
AL 1.671 1.324 1.029

95% CP 0.923 0.903 0.935

(70, 80)

5.0

θ̂
AL 0.701 0.456 0.653

95% CP 0.975 0.943 0.985

σ̂
AL 0.895 0.678 0.490

95% CP 0.925 0.930 0.965

λ̂
AL 1.440 0.896 0.984

(55, 65) 95% CP 0.955 0.939 0.965

(60, 75)

10.0

θ̂
AL 0.610 0.568 0.537

95% CP 0.985 0.936 0.966

σ̂
AL 0.891 0.783 0.483

95% CP 0.968 0.936 0.973

λ̂
AL 1.576 0.963 0.973

95% CP 0.966 0.936 0.956



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 26 of 34

3.2. Illustrative Example and real data analysis

In this section, we introduce a numerical investigation of the estimation methods dis-
cussed in previous sections for the NPT distribution using simulated data and a real data
set.

3.2.1. Type-I GHCS data

Here, the estimation methods described in the previous sections is applied to the set of
simulated type-I GHCS data under the constant-stress PALT. A data set of system lifetime
is generated from the NPT model with θ = 0.5, σ = 1.0 and λ = 1.5, respectively. Based
on (n1 = n2 = 60, l1 = l2 = 40,m1 = m2 = 50, T = 10), using the algorithm described in
Algo.(3). We simulate two samples of size D1 = 40 and D2 = 42 from the NPT(θ, σ) and
NPT(θ, σ, λ). The simulated data are given in Tab.(5).

Table 5: Simulated type-I GHC samples with constant-stress PALT.

blue

Normal condition 1.1471.1491.2061.2891.3101.4852.0212.1302.1492.2332.267

2.3242.3942.6322.6623.1663.4364.0724.1044.9185.2675.522

5.7786.4507.2198.4138.5279.2889.60010.3210.3411.4812.97

18.7220.1423.3426.7033.3437.0938.81

Accelerated condi-
tion

1.0391.1031.2351.2601.2801.3171.9662.0392.0822.1462.333

2.4102.4432.5652.6672.7592.8373.0283.0673.0903.5493.755

3.8073.9454.2344.3694.6164.8604.9455.1315.5195.9566.403

6.6086.8627.0307.0867.1487.2257.5798.3278.404

Based on the observed censored data in Tab.(5) and the methods adopted in Subs. (2.2.3),
we calculated the MLEs of θ, σ and λ. The values are shown in Tab.(6). In order to obtain
the Bayes estimates of θ, σ and λ, 10, 000 samples were generated based on the Metropo-
lis–Hastings algorithm. These Bayes estimates were obtained after discarding the initial
1000 burn-in iterative values. Bayes estimates were calculated with the informative priors,
and the hyper parameters as (a = 0.002, b = 2.0, c = 5.0, d = 2.0). For LLF, h∗ = −0.8
and h∗ = 0.8 were used as two options for the constant h∗. These choices give more weight
to underestimation and exaggeration, respectively. The MLEs of the unknown parame-
ters were considered as their initial values, while the diagonal elements of the reciprocal
of the observed Fisher information matrix were considered as the variance of the MLEs.
The estimates are displayed in Tab.(7). The 95% N-boot CIs, together with the HPD
credible intervals, are presented in Tab.(7). The numerical illustration demonstrates that
the means of the Bayes estimates and maximum likelihood estimates (MLEs) obtained by



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 27 of 34

Table 6: Estimated values of θ, σ, and λ.

Method−→ MLEs Bayesian

Parameter↓ NR EM SELF
LLF

h∗ = −0.8 h∗ = 0.8

θ = 0.5 Estimate 0.4544 0.4643 0.4706 0.4722 0.4691

σ = 1.0 Estimate 1.0390 1.0648 1.1499 1.1528 1.1469

λ = 1.5 Estimate 1.5732 1.4408 1.7870 1.8489 1.7324

Table 7: 95% confdence interval (CI) estimates of θ, σ and λ.

Method−→ MLEs N-boot Bayesian

θ = 0.5 95% CI (0.4079, 0.5010) (0.3704, 1.3103) (0.3567, 0.5448)

σ = 1.0 95% CI (0.6944, 1.3836) (0.9629, 2.0927) (0.9812, 1.2213)

λ = 1.5 95% CI (1.3264, 2.0545) (1.4541, 2.6006) (1.1472, 2.2459)

the Newton-Raphson (NR) or expectation-maximization (EM) algorithm for the unknown
parameters θ, σ, and λ exhibit a high degree of proximity to the true values.

3.2.2. Insulating Fluid

This application provides an analysis of the time-to-breakdown (measured in seconds) of
an insulating fluid during a voltage endurance test conducted under varying stress levels.
The longevity of various industrial components is contingent upon the durability of their
electrical insulation. Enhancing the efficiency of the life test can be achieved by imple-
menting accelerating circumstances, such as increasing the voltage. Ref.[43] conducted a
study on accelerated tests, wherein varying magnitudes of voltage were applied to distinct
experimental groups. According to Ref.[43], there are two stress levels that have been
identified; 30 Kilovolt, which represents typical use, and 32 Kilovolt, which represents
rapid stress. The data is shown in Tab.(8).

Table 8: oil breakdown times of insulating fluid.

Constant Stress Condition

194.90, 175.88, 144.12, 139.07, 47.30, 43.40, 22.66, 21.02, 20.46, 17.05, 7.74

Accelerated Stress Condition

215.10, 100.58, 89.29, 82.85, 53.24, 27.80, 15.93, 13.95, 9.88, 3.91, 2.75, 0.79, 0.69, 0.40, 0.27

Before performing parameter estimation, the adequacy of the dataset’s fit to the NPT dis-
tribution was assessed. The Kolmogorov-Smirnov test (KS test) was employed as a means
to quantify the deviation of the data from the underlying distribution. The KS distance
and p-value of the data obtained from the constant stress condition were determined
to be 0.2146 and 0.6178, respectively. Similarly, the KS distance and p-value of the
data obtained from the accelerated stress condition were found to be 0.2179 and 0.4149,



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 28 of 34

respectively. Thus, we ensured that the NPT distribution is considered an appropriate dis-
tribution for this data. For further illustration, the empirical cumulative distribution plot
with fitted theoretical NPT distribution CDF where the model parameters is estimated
from the complete real life data based on maximum likelihood method. In this goodness-
of-fit approach, probability-probability (P−P) and Quantile-Quantile (Q−Q) plots are
also provided in Figs.(2) and (3) which indicates that the NPT distribution provides a
good fit to the data as well.

50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

x1

F
1
(x

1
)

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

5
0

1
0
0

1
5
0

2
0
0

Theoretical Quantiles

S
a
m

p
le

 Q
u
a
n
ti
le

s

Estimated NP distribution

E
m

p
ir

ic
a
l 
c
u
m

u
la

ti
v
e
 d

is
tr

ib
u
ti
o
n

0.0 0.2 0.4 0.6 0.8 1.0

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

Figure 2: ECD, P-P, and Q-Q plots for normal conditions.

0 50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

0 50 100 150 200

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

x2

F
2
(x

2
)

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

0
5
0

1
0
0

1
5
0

2
0
0

Theoretical Quantiles

S
a
m

p
le

 Q
u
a
n
ti
le

s

Estimated NP distribution

E
m

p
ir

ic
a
l 
c
u
m

u
la

ti
v
e
 d

is
tr

ib
u
ti
o
n

0.0 0.2 0.4 0.6 0.8 1.0

0
.0

0
.2

0
.4

0
.6

0
.8

1
.0

Figure 3: ECD, P-P, and Q-Q plots for accelerated conditions.

The maximum likelihood estimates (MLEs) are presented in Tab.(9).In the absence of
prior knowledge about the unknown population parameters, adopting Bayes estimates is
justified. In this scenario, it is postulated that the prior distributions of θ and σ are
improper, specifically with parameters a = −1, b = 0, c = 1, and d = ∞, as there is a
lack of prior knowledge available. As previously mentioned, the Gibbs method utilized
Metropolis within it to produce a total of 11, 000 Markov Chain Monte Carlo (MCMC)
samples. These samples were generated using the maximum likelihood estimators (MLEs)
of the parameters θ, σ, and λ as the initial values at the onset of the procedure. Fig.(5)
displays trace plots for the initial 1000 Markov Chain Monte Carlo (MCMC) iterations of
the parameters θ, σ, and λ. The MCMC technique shows strong convergence. Further-
more, Fig.(4) displays the histogram plots of the generated samples for the parameters θ,
σ, and λ. The histograms of the obtained samples exhibit a strong resemblance to the
theoretical posterior density functions. For Bayesian estimations under asymmetric LLF,



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 29 of 34

it is known in the literature that, h∗ < 0 implies that underestimation results in more
penalty than overestimation and the reverse is true for h∗ > 0. When h∗ close to zero, the
LLF becomes symmetric and behaves roughly like the SELF. The MLEs relative to both
NR, EM techniques and Bayesian MCMC method with chosen h∗ = −3 and 3 of unknown
parameters are computed and listed in Tab.(9). Moreover, the results of 95% approximate
CI, N-boot and HPD intervals of unknown parameters are given in Tab.(10).

Table 9: Estimated values of θ, σ and λ from insulating fluid.

Method−→ MLEs Bayesian

Parameter↓ NR EM SELF
LLF

h∗ = −0.3 h∗ = 0.3

θ Estimate 0.2279 0.2111 0.2241 0.2233 0.2249
σ Estimate 0.2700 0.2734 0.3336 0.3316 0.3362
λ Estimate 1.4243 1.6554 1.7704 1.6526 1.9369

Table 10: 95% confdence interval (CI) estimates of θ, σ and λ from insulating fluid.

Method MLEs N-boot Bayesian

θ 95% CI (0.1544, 0.3015) (0.1676, 0.5900) (0.1008, 0.3144)

σ 95% CI (0.0252, 0.5652) (0.0037, 0.7070) (0.1014, 0.5016)

λ 95% CI (0.0777, 2.7709) (0.1113, 3.6775) (0.5968, 2.9076)



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 30 of 34

In conclusion, it can be inferred that the estimated distribution of the non parametric
NPT model has a strong alignment with the provided data.

0 4000 8000

0
.1

0
.3

0
.5

Iterations

θ

θ
D

e
n

s
it
y

0.1 0.3 0.5

0
1

2
3

4
5

6
Figure 4: Density (right) and Trace (left) plots of θ from insulating fluid data.

0 4000 8000

0
.0

0
.2

0
.4

0
.6

Iterations

σ

σ

D
e

n
s
it
y

0.0 0.2 0.4 0.6 0.8

0
.0

1
.0

2
.0

3
.0

Figure 5: Density (right) and Trace (left) plots of σ from insulating fluid data.

4. Conclusions

This study discusses multiple parameter estimation approaches to estimate the two
parameters of the NPT distribution. The estimation methods are applied to a constant
stress partially accelerated life test, utilising type-I GHCS data. The maximum likelihood
estimators (MLEs) and Bayes estimates of the parameters were obtained, along with the
calculation of the acceleration factor and the related confidence intervals. Additionally,
the EM technique has been utilised to derive the maximum likelihood estimators (MLEs)
for the undetermined parameters. The associated Hessian matrix is also presented in



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 31 of 34

0 4000 8000

0
1

2
3

4
5

Iterations

λ

λ

D
e

n
s
it
y

0 1 2 3 4 5

0
.0

0
.2

0
.4

0
.6

Figure 6: Density (right) and Trace (left) plots of λ from insulating fluid data.

the paper. In the absence of explicit formulations for the maximum likelihood estima-
tors (MLEs) of certain parameters, we rely on the Bayesian approach for assistance. The
Bayes estimates are derived under the assumption of dependent priors, using the two loss
functions(squared error (SE) and linear exponential (LINEX)). The posterior distribu-
tions of the unknown parameters suggest that certain parameters do not conform to a
widely recognised distribution. As a result, we employed Metropolis-Hastings sampling
in the Gibbs sampling procedure to compute the Bayes estimates together with their
corresponding credible intervals. Moreover, the empirical findings indicate that the uti-
lization of informative priors in the Bayes approach yielded superior outcomes compared
to the maximum likelihood (ML) technique, regardless of whether the Newton-Raphson
(NR) or expectation-maximization (EM) algorithm was employed. Furthermore, it is well
acknowledged that the estimates obtained using the Expectation-Maximization (EM) ap-
proach tend to exhibit superior performance when compared to those acquired using the
Newton-Raphson (NR) method. This is seen in the significantly lower values of mean
squared errors (MSEs) and average widths of the interval estimators. In this study, our
primary focus has been on type-I GHCS and NPT distribution. However, it is worth
noting that the methodology discussed can also be applied to other distribution and cen-
soring schemes. There are several additional tasks that can be pursued in this area. These
topics present potential avenues for future investigation. In conclusion, we propose the
utilization of Markov Chain Monte Carlo (MCMC) and Expectation-Maximization (EM)
methodologies in conjunction with partially accelerated life testing techniques, specifically
employing type-I GHCS data, for the purpose of life testing and reliability modeling.

References

[1] David Collett. Modelling survival data in medical research. CRC Press, Boca Raton,
FL, 2023.

[2] Chin Diew Lai and Min Xie. Stochastic ageing and dependence for reliability. Springer,
New York, 2006.



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 32 of 34

[3] Mousa Abdi, Akbar Asgharzadeh, Hassan S. Bakouch, and Zahra Alipour. A new
compound gamma and Lindley distribution with application to failure data. Austrian
Journal of Statistics, 48(3):54–75, 2019.

[4] Romano Demicheli, Gianni Bonadonna, William J. M. Hrushesky, Michael W. Retsky,
and Pinuccia Valagussa. Menopausal status dependence of the timing of breast cancer
recurrence after surgical removal of the primary tumour. Breast Cancer Research,
6(6):689–696, 2004.

[5] Marcelo Bourguignon, Helton Saulo, and Rodrigo Nobre Fernandez. A new Pareto-
type distribution with applications in reliability and income data. Physica A: Statis-
tical Mechanics and its Applications, 457:166–175, 2016.

[6] Paduthol Godan Sankaran, N. Unnikrishnan Nair, and Preethi John. A family of
bivariate Pareto distributions. Statistica, 74(2):199–215, 2014.

[7] Ali Saadati Nik, Akbar Asgharzadeh, and Saralees Nadarajah. Comparisons of meth-
ods of estimation for a new Pareto-type distribution. Statistica, 79(3):291–319, 2019.

[8] A. Saadati Nik, Akbar Asgharzadeh, and Mohammad Z. Raqab. Estimation and
prediction for a new Pareto-type distribution under progressive type-II censoring.
Mathematics and Computers in Simulation, 190:508–530, 2021.

[9] José Maŕıa Sarabia, Vanesa Jorda, and Faustino Prieto. On a new Pareto-type distri-
bution with applications in the study of income inequality and risk analysis. Physica
A: Statistical Mechanics and its Applications, 527:121277, 2019.

[10] Kadir Karakaya, Yunus Akdoğan, A. Saadati Nik, Coşkun Kuş, and Akbar As-
gharzadeh. A generalization of new Pareto-type distribution. Annals of Data Science,
9:1–15, 2022.

[11] A. Saadati Nik, Akbar Asgharzadeh, and Ayman Baklizi. Inference based on new
Pareto-type records with applications to precipitation and COVID-19 data. Statistics,
Optimization & Information Computing, 11(2):243–257, 2023.

[12] Morris H. DeGroot and Prem K. Goel. Bayesian estimation and optimal designs in
partially accelerated life testing. Naval Research Logistics Quarterly, 26(2):223–235,
1979.

[13] Vilijandas Bagdonavicius and Mikhail Nikulin. Accelerated life models: modeling and
statistical analysis. Chapman and Hall/CRC, Boca Raton, FL, 2001.

[14] Ayon Ganguly and Debasis Kundu. Analysis of simple step-stress model in presence
of competing risks. Journal of Statistical Computation and Simulation, 86(10):1989–
2006, 2016.

[15] Man Ho Ling and X. W. Hu. Optimal design of simple step-stress accelerated life tests
for one-shot devices under Weibull distributions. Reliability Engineering & System
Safety, 193:106630, 2020.

[16] N. A. Abou-Elheggag, Al-Wageh A. Farghal, G. A. Abd-Elmougod, and Osama M.
Taha. Progressive first-failure censored samples in estimation and prediction of NH
distribution. Journal of Statistics Applications & Probability, 10(3):717–731, 2021.

[17] Abdullah Ali H. Ahmadini, Wali Khan Mashwani, Rehman Ahmad Khan Sherwani,
Shokrya S. Alshqaq, Farrukh Jamal, Miftahuddin Miftahuddin, Kamran Abbas, Faiza
Razaq, Mohammed Elgarhy, and Sanaa Al-Marzouki. Estimation of constant stress



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 33 of 34

partially accelerated life test for Fréchet distribution with type-I censoring. Mathe-
matical Problems in Engineering, 2021:5590406, 2021.

[18] Al-Wageh A. Farghal, Souha K. Badr, Hanaa Abu-Zinadah, and Gamal A. Abd-
Elmougod. Analysis of generalized inverted exponential competing risks model in
presence of partially observed failure modes. Alexandria Engineering Journal, 78:74–
87, 2023.

[19] Seunggeun Hyun and Jimin Lee. Constant-stress partially accelerated life testing
for log-logistic distribution with censored data. Journal of Statistics Applications &
Probability, 4(2):193–201, 2015.

[20] Nagwa M. Mohamed. Estimation on Kumaraswamy-inverse Weibull distribution with
constant stress partially accelerated life tests. Applied Mathematics & Information
Sciences, 15(4):503–510, 2021.

[21] Mazen Nassar and Farouq Mohammad A. Alam. Analysis of modified Kies expo-
nential distribution with constant stress partially accelerated life tests under type-II
censoring. Mathematics, 10(5):819, 2022.

[22] Benjamin Epstein. Truncated life tests in the exponential case. The Annals of Math-
ematical Statistics, 25(3):555–564, 1954.

[23] Nader Ebrahimi. Estimating the parameters of an exponential distribution from a
hybrid life test. Journal of Statistical Planning and Inference, 14(2-3):255–261, 1986.

[24] B. Chandrasekar, A. Childs, and N. Balakrishnan. Exact likelihood inference for the
exponential distribution under generalized type-I and type-II hybrid censoring. Naval
Research Logistics, 51(7):994–1004, 2004.

[25] Abdulaziz S. Alghamdi. Statistical inferences of competing risks generalized half-
logistic lifetime populations in presence of generalized type-I hybrid censoring scheme.
Alexandria Engineering Journal, 65:699–708, 2023.

[26] Baria A. Helmy, Amal S. Hassan, Ahmed K. El-Kholy, Rashad A. R. Bantan, and
Mohammed Elgarhy. Analysis of information measures using generalized type-I hybrid
censored data. Journal of Statistical Theory and Applications, 21(4):229–249, 2022.

[27] Laila A. Al-Essa, Ahmed A. Soliman, Gamal A. Abd-Elmougod, and Huda M. Al-
shanbari. Comparative study with applications for Gompertz models under competing
risks and generalized hybrid censoring schemes. Axioms, 12(10):973, 2023.

[28] Abdalla Rabie and Junping Li. E-Bayesian estimation for Burr-X distribution based
on generalized type-I hybrid censoring scheme. American Journal of Mathematical
and Management Sciences, 39(1):41–55, 2020.

[29] Ahmed A. Soliman, Gamal A. Abd-Elmougod, Alwageh Ahmed, and Osama Mo-
hamed Taha. Statistical inference of a new Pareto-type model under generalized
hybrid type-I censored samples. Sohag Journal of Sciences, 10(1):10–23, 2025.

[30] Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin. Maximum likelihood
from incomplete data via the EM algorithm. Journal of the Royal Statistical Society:
Series B (Methodological), 39(1):1–22, 1977.

[31] Geoffrey J. McLachlan and Thriyambakam Krishnan. The EM algorithm and exten-
sions. John Wiley & Sons, Hoboken, NJ, 2007.

[32] Richard A. Askey and Adri B. Olde Daalhuis. Generalized hypergeometric func-



A. A. Soliman et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5788 34 of 34

tions and Meijer G-function. NIST Digital Library of Mathematical Functions, 2010.
Available at https://dlmf.nist.gov.

[33] Debasis Kundu and Biswabrata Pradhan. Estimating the parameters of the gener-
alized exponential distribution in presence of hybrid censoring. Communications in
Statistics—Theory and Methods, 38(12):2030–2041, 2009.

[34] Thaung Lwin. Estimation of the tail of the Paretian law. Scandinavian Actuarial
Journal, 1972(2):170–178, 1972.

[35] Barry C. Arnold and S. James Press. Bayesian estimation and prediction for Pareto
data. Journal of the American Statistical Association, 84(408):1079–1084, 1989.

[36] Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin. Bayesian data
analysis. Chapman and Hall/CRC, Boca Raton, FL, 1995.

[37] W. Keith Hastings. Monte Carlo sampling methods using Markov chains and their
applications. Biometrika, 57(1):97–109, 1970.

[38] Nicholas Metropolis and Stanislaw Ulam. The Monte Carlo method. Journal of the
American Statistical Association, 44(247):335–341, 1949.

[39] Michael R. Chernick. Bootstrap methods: a practitioner’s guide. Wiley, New York,
1999.

[40] Yunus Akdoğan. On the confidence intervals of process capability index Cpm based
on a progressive type-II censored sample. Quality and Reliability Engineering Inter-
national, 38(5):2845–2861, 2022.

[41] Arne Henningsen and Ott Toomet. maxLik: a package for maximum likelihood esti-
mation in R. Computational Statistics, 26(3):443–458, 2011.

[42] Martyn Plummer, Nicky Best, Kate Cowles, Karen Vines, et al. CODA: convergence
diagnosis and output analysis for MCMC. R News, 6(1):7–11, 2006.

[43] Wayne B. Nelson. Accelerated testing: statistical models, test plans, and data analysis.
John Wiley & Sons, Hoboken, NJ, 2009.