**1. Introduction**

Modeling heavy-tailed data is one of the important aspects in many engineering and medical domains. Initial work on this topic was carried out by Pareto [1] to model income data. In later years, the applications of Pareto, particularly the type II Lomax distribution (see [2]), usually referred to as Lomax (Lx) distribution, branched into scientific fields such as engineering sciences, actuarial sciences, medicine, income, and many more. The distribution function (cdf) and probability density function (pdf) of the Lx distribution are given by

$$F\_{Lx}(\mathbf{x}; \boldsymbol{\eta}) = 1 - \left(\frac{\boldsymbol{\beta}}{\mathbf{x} + \boldsymbol{\beta}}\right)^{\alpha}$$

and

$$f\_{Lx}(x;\eta) = \frac{a}{\beta} \left(\frac{\beta}{x+\beta}\right)^{a+1}, \quad x \ge 0, \eta = (a,\beta) > 0,$$

respectively, where *α* is a shape parameter, and *β* is a scale parameter. We have *FLx*(*x*; *η*) = *fLx*(*x*; *η*) = 0 for *x* < 0. References [3,4] considered the Lx distribution to model income and wealth data. Reference [5] used the Lx distribution as an alternative to the exponential, gamma, and Weibull distributions for heavy-tailed data. Reference [6] derived various estimation techniques based on the Lx distribution. References [7,8] examined the various structural properties and record value moments of the Lx distribution. Reference [9] extensively studied and extended the family of distributions that were used in the Lx distribution. Reference [10] considered the Lx distribution as an important distribution to model lifetime data, since it belongs to the family of decreasing hazard rate.

**Citation:** Nagarjuna, V.B.V.; Vardhan, R.V.; Chesneau, C. Nadarajah–Haghighi Lomax Distribution and Its Applications. *Math. Comput. Appl.* **2022**, *27*, 30. https://doi.org/10.3390/ mca27020030

Received: 1 December 2021 Accepted: 18 March 2022 Published: 1 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

133

In continuation of this, many researchers have proposed several distributions that deal with heavy-tailed data by generalizing the functional forms of the Lx distribution. It mainly consists of adding scale/shape parameters accordingly. A few to mention are the exponentiated Lx (EL) distribution in [11], beta Lx (BL) distribution in [12], Poisson Lx distribution in [13], exponential Lx (EXL) distribution in [14], gamma Lx (GL) distribution in [15], Weibull Lx (WL) distribution in [16], beta exponentiated Lx distribution in [17], power Lx distribution in [18], exponentiated Weibull Lx distribution in [19], Marshall– Olkin exponential Lx distribution in [20], type II Topp–Leone power Lx distribution in [21], Marshall–Olkin length biased Lomax distribution in [22], Kumaraswamy generalized power Lx distribution in [23] and sine power Lx distribution in [24]. For the purpose of this study, a retrospective on the EXL distribution is required. To begin, it is defined by the following cdf and pdf:

and

$$f\_{\text{EXL}}(\mathbf{x}; \boldsymbol{\xi}) = \frac{\lambda a}{\beta} \left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a+1} e^{-\lambda \left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}}, \quad \mathbf{x} \ge -\beta, \boldsymbol{\xi} = (a, \beta, \lambda) > 0, \mathbf{x}$$

−*λ β x*+*β* −*α*

*FEXL*(*x*; *ξ*) = 1 − *e*

respectively, where *α* is a shape parameter, and *β* and *λ* are scale parameters. We have *FEXL*(*x*; *ξ*) = *fEXL*(*x*; *ξ*) = 0 for *x* < −*β*. Thus, the EXL distribution combines the functionalities of the exponential and Lx distributions through a specific composition scheme. This scheme may be called the extended Lx scheme (it will be discussed mathematically later). As immediate remarks, the EXL distribution has three parameters and is with a lower bounded support. It is shown in [14] that the pdf of the EXL distribution is unimodal and has an increasing hazard rate function (hrf). Moreover, its quantile and moment properties are manageable. On the statistical side, by considering the aircraft windshield data collected in [25], it is proven in [14] that the EXL model outperforms several three- or four-parameter extensions of the Lx model, including the EL, BL, and GL models. Thus, strong evidence is for the use of the extended Lx scheme for the construction of efficient distributions and models.

On the other hand, recently, a generalized version of the exponential distribution was given by Nadarajah and Haghighi [26]. It can be presented as an alternative to the Weibull, gamma, and exponentiated exponential (EE) distributions. It is called the Nadarajah– Haghighi (NH) distribution. The cdf and pdf of the NH distribution are

$$F\_{NH}(\mathbf{x};\mathbf{\tau}) = 1 - \mathbf{e}^{1 - (1 + b\mathbf{x})^a}$$

and

$$f\_{NH}(\mathbf{x};\tau) = ab(1+b\mathbf{x})^{a-1}e^{1-(1+b\mathbf{x})^a}, \quad \mathbf{x}\geq 0, \ \tau = (a,b)>0,$$

respectively. We have *FNH*(*x*; *τ*) = *fNH*(*x*; *τ*) = 0 for *x* < 0. Among its main features, the pdf can have decreasing and uni-modal shapes, and the hrf exhibits increasing, decreasing, and constant shapes. According to [26], if the pdfs of the gamma, Weibull, and exponentiated exponential are monotonically decreasing, then it is not possible to allow increasing hrf. However, such a hrf property can be achieved by the NH distribution.

In light of the above research work, we present a new distribution based on the extended Lx scheme with the use of the NH distribution as the main generator. It is called the NH Lx (NHLx) distribution. In this sense, the NHLx distribution is to the NH distribution what the EXL distribution is to the exponential distribution. The NHLx distribution can also be presented as a generalization of the EXL distribution through the introduction of an additional shape parameter. We investigate the theoretical and practical facets of the NHLx distribution. Among its functional features, it has four parameters, it is lower-bounded (as with the EXL distribution, with a bound governed by a scale parameter), its pdf exhibits non-increasing and inverted J-shaped curves, and its hrf possesses increasing, decreasing, and upside-down bathtub shapes. This combination of qualitative characteristics is rare for a lower-bounded distribution and, in this way, it has better functionality to model

lifetime data than the EXL and Lx distributions, among others. We illustrate this aspect by considering four different data sets referenced in the literature.

The rest of the article covers the following aspects: Section 2 presents the most important functions of the NHLx distribution, namely the cdf, pdf, hrf, and quantile function (qf), along with a graphical analysis when necessary. Section 3 is devoted to moment analysis and related functions. Section 4 concerns the maximum likelihood estimates of the NHLx model parameters. The above section is completed by a simulation study in Section 5. Concrete applications of the NHLx model are developed in Section 6. A conclusion is formulated in Section 7.

#### **2. NHLx Distribution**

In order to understand the essence of the NHLx distribution, let us describe more precisely the extended Lx scheme on the basis of the EXL distribution. One can remark that *FEXL*(*x*; *ξ*) = *FE* <sup>1</sup> 1−*F*<sup>∗</sup> *Lx*(*x*;*η*); *<sup>λ</sup>* , where *FE*(*x*; *λ*) denotes the cdf of the exponential distribution with parameter *λ*, and *F*∗ *Lx*(*x*; *η*) = 1 − *β x*+*β α* for *x* ≥ −*β*, and *F*<sup>∗</sup> *Lx*(*x*; *η*) = 0 otherwise. Thus, *F*∗ *Lx*(*x*; *η*) can be thought of as a support-extended version of *FLx*(*x*; *η*) over the semi-finite interval [−*β*, ∞). It is worth noting that *F*<sup>∗</sup> *Lx*(*x*; *η*) is not a cdf anymore, but it is increasing and satisfies lim*x*→−*<sup>β</sup>* <sup>1</sup> 1−*F*<sup>∗</sup> *Lx*(*x*;*η*) <sup>=</sup> 0 and lim*x*→<sup>∞</sup> <sup>1</sup> 1−*F*<sup>∗</sup> *Lx*(*x*;*η*) <sup>=</sup> <sup>∞</sup>, which ensure that *FEXL*(*x*; *ξ*) as a cdf is mathematically correct. It is worth noting that it can be applied to any lifetime distribution in place of the generator exponential distribution.

Based on the extended Lx scheme with the NH distribution as a generator, the cdf and pdf of the NHLx distribution are specified by

$$F\_{NHLx}(\mathbf{x}; \boldsymbol{\zeta}) = 1 - e^{1 - \left(1 + b\left(\frac{\beta}{x+\beta}\right)^{-\alpha}\right)^{\alpha}}$$

and

$$\begin{aligned} f\_{NHLx}(\mathbf{x};\boldsymbol{\zeta}) &= \frac{ab\boldsymbol{\alpha}}{\beta} \left(\frac{\boldsymbol{\beta}}{\mathbf{x}+\boldsymbol{\beta}}\right)^{-a+1} \left(1+b\left(\frac{\boldsymbol{\beta}}{\mathbf{x}+\boldsymbol{\beta}}\right)^{-a}\right)^{a-1} \boldsymbol{\varepsilon}^{1-\left(1+b\left(\frac{\boldsymbol{\beta}}{\mathbf{x}+\boldsymbol{\beta}}\right)^{-a}\right)^{a}} \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \boldsymbol{\varepsilon} \ge -\beta \; \; \boldsymbol{\zeta} = (\boldsymbol{a},\boldsymbol{b},\boldsymbol{a},\boldsymbol{\beta}) > 0, \end{aligned}$$

respectively, where *a* and *α* are shape parameters, and *b* and *β* are scale parameters. We have *FNHLx*(*x*; *ζ*) = *fNHLx*(*x*; *ζ*) = 0 for *x* < −*β*. Thus, the cdf has been derived from the following formula: *FEXL*(*x*; *<sup>ξ</sup>*) = *FNH* <sup>1</sup> 1−*F*<sup>∗</sup> *Lx*(*x*;*η*); *<sup>τ</sup>* , *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>. By taking *<sup>a</sup>* <sup>=</sup> 1, we remark that *FNHLx*(*x*; *ζ*) = *FEXL*(*x*; *ξ*); the NHLx distribution is reduced to the EXL distribution with *λ* = *b*. The asymptotic properties of the pdf depend on the values of *α* mainly; with the use of standard asymptotic techniques, we establish that

$$\lim\_{\substack{\mathbf{x}\to-\boldsymbol{\beta}}} f\_{NHLx}(\mathbf{x};\boldsymbol{\zeta}) = \begin{cases} \begin{array}{cc} \infty & \text{if } \boldsymbol{\alpha} < 1 \\ \frac{\mathrm{Id}}{\boldsymbol{\beta}} & \text{if } \boldsymbol{\alpha} = 1 \\ 0 & \text{if } \boldsymbol{\alpha} > 1 \end{array} , \quad \lim\_{\substack{\mathbf{x}\to\boldsymbol{\infty}}} f\_{NHLx}(\mathbf{x};\boldsymbol{\zeta}) = 0. \end{cases}$$

Figure 1 completes these asymptotic results by showing some curves of the pdf for several parameter values.

In Figure 1, we see that the pdf can be inverted J decreasing or have uni-modal shapes. It is very flexible to skewness, peakedness, and platness curves at a small value of *β* (at least), and different selected parameter values of *a*, *b*, and *α*. Such flexibility is not observed for the pdf of the EXL distribution, as visually shown in the figures in [14].

The analysis of the corresponding hrf is now examined. By applying the definition *hNHLx*(*x*; *ζ*) = *fNHLx*(*x*; *ζ*)/[1 − *FNHLx*(*x*; *ζ*)], it is given by

$$h\_{NHLx}(\mathbf{x}; \boldsymbol{\xi}) = \frac{ab\mathbf{a}}{\beta} \left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a+1} \left(1 + b\left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}\right)^{a-1}, \quad \mathbf{x} \ge -\beta\sqrt{\beta}$$

and *hNHLx*(*x*; *ζ*) = 0 for *x* < −*β*. Contrary to the pdf, the asymptotic properties of the hrf mainly depend on the values of *a* and *α*; we have

$$\lim\_{\mathbf{x}\to-\mathbf{\tilde{\beta}}} h\_{NHLx}(\mathbf{x};\mathbf{\tilde{\zeta}}) = \begin{cases} \infty & \text{if } a < 1\\ \frac{ba}{\mathbf{\tilde{\beta}}} & \text{if } a = 1\\ 0 & \text{if } a > 1 \end{cases}, \quad \lim\_{\mathbf{x}\to\mathbf{\infty}} h\_{NHLx}(\mathbf{x};\mathbf{\tilde{\zeta}}) = \begin{cases} \infty & \text{if } aa > 1\\ \frac{b^a}{\mathbf{\tilde{\beta}}} & \text{if } aa = 1\\ 0 & \text{if } aa < 1 \end{cases}$$

.

In full generality, the possible shapes of the hrf are determinant for modeling purposes: the more different shapes it has, the more the associated model is applicable to a wide panel of data sets.

**Figure 1.** Curves of the pdf of the NHLx distribution for various parameter values, but with the fixed value: *β* = 0.005.

Figure 2 presents the identified shapes for the hrf of the NHLx distribution. From Figure 2, we see that the hrf can be increasing, decreasing, or upside-down bathtub-shaped, with flexible convex–concave properties. In particular, these curve modulations are possible thanks to the variation of the new additional parameters *a*. We are far beyond the curve possibilities of the hrf of the EXL distribution, which is only increasing according to [14]. Thus, from one perspective, the NHLx distribution adds a new shape parameter *a* to the EXL distribution in a thorough fashion, considerably improving its modeling properties.

The qf of the NHLx distribution is now studied. To begin, it is defined in function of *FNHLx*(*u*; *ζ*) by *QNHLx*(*u*; *ζ*) = *F*−<sup>1</sup> *NHLx*(*u*; *ζ*), *u* ∈ (0, 1). After some mathematical development, we establish that

$$Q\_{NHLx}(u; \mathcal{J}) = \beta \left\{ \left[ \frac{1}{b} \left( (1 - \log(1 - u))^{\frac{1}{\sigma}} - 1 \right) \right]^{\frac{1}{\sigma}} - 1 \right\}, \quad u \in (0, 1).$$

Based on this qf, the main quartiles of the NHLx distribution can be explicated: by taking *u* = 1/4, *u* = 1/2, and *u* = 3/4 into *QNHLx*(*u*; *ζ*), we get the first, second, and third quartiles. In addition, several quantile-based functions, and skewness and kurtosis measures,

can be listed and analyzed (see [27]). In addition, various quantile regression models can be constructed (see [28]).

**Figure 2.** Curves of the hrf of the NHLx distribution for various parameter values.

#### **3. Moment Properties of the NHLx Distribution**

The moment properties of the NHLx distribution are now under investigation. First, for a random variable *X* with the NHLx distribution and any integer *r*, the *r*th moment of *X* is defined by

$$\mu\_r^\* = E(X^r) = \int\_{-\infty}^{\infty} x^r f\_{NHLx}(x;\zeta) d\infty.$$

which can be explicated as

$$\mu\_r^\* = \int\_{-\beta}^{\infty} \mathbf{x}^r \frac{ab\mathbf{a}}{\beta} \left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a+1} \left(1 + b\left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}\right)^{a-1} e^{1 - \left(1 + b\left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}\right)^a} d\mathbf{x} \dots$$

For given distribution parameters, this integral can be computed numerically with the help of scientific software. An analytical expression involving sums is given in the next proposition.

**Proposition 1.** *Let X be a random variable with the NHLx distribution. Then, its rth moment can be expressed as*

$$\mu\_r^\* = \beta^r \varepsilon \sum\_{j=0}^r \sum\_{k=0}^\infty \binom{r}{j} \binom{\frac{j}{a}}{k} \frac{(-1)^{r-j+k}}{b^{\frac{j}{a}}} \Gamma\left(\frac{1}{a}\left(\frac{j}{a} - k\right) + 1, 1\right),$$

*where e* = exp(1) *and* Γ(*x*, *y*) = <sup>∞</sup> *<sup>x</sup> t <sup>y</sup>*−1*e*−*<sup>t</sup> dt with <sup>y</sup>* <sup>∈</sup> <sup>R</sup> *and <sup>x</sup>* <sup>&</sup>gt; <sup>0</sup>*, which defines the incomplete gamma function.*

**Proof.** Let us apply the following change of variable:

$$u = \left(1 + b\left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}\right)^a, \quad du = \frac{aba}{\beta} \left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a+1} \left(1 + b\left(\frac{\beta}{\mathbf{x} + \beta}\right)^{-a}\right)^{a-1} dx, \quad \mathbf{x} = \frac{\mathbf{a} \cdot \mathbf{x} + \mathbf{a} \cdot \mathbf{x}}{|\mathbf{a} + \mathbf{a}|}$$

which satisfies lim*x*→−*<sup>β</sup> u* = 1 and lim*x*→<sup>∞</sup> *u* = ∞. Then, we have

$$
\mu\_r^\* = \beta^r e \int\_1^\infty \left\{ \left[ \frac{1}{b} \left( u^{\frac{1}{a}} - 1 \right) \right]^{\frac{1}{a}} - 1 \right\}^r e^{-u} du.
$$

By applying the standard and generalized binomial formulas, we get

$$\begin{split} \mu\_r^\* &= \beta^r \varepsilon \sum\_{j=0}^r \binom{r}{j} \frac{(-1)^{r-j}}{b^{\frac{j}{\kappa}}} \int\_1^\infty \left( u^{\frac{1}{\kappa}} - 1 \right)^{\frac{j}{\kappa}} e^{-u} du \\ &= \beta^r \varepsilon \sum\_{j=0}^r \binom{r}{j} \frac{(-1)^{r-j}}{b^{\frac{j}{\kappa}}} \sum\_{k=0}^\infty \binom{\frac{j}{\kappa}}{k} (-1)^k \int\_1^\infty u^{\frac{1}{\kappa} \left( \frac{j}{\kappa} - k \right)} e^{-u} du \\ &= \beta^r \varepsilon \sum\_{j=0}^r \sum\_{k=0}^\infty \binom{r}{j} \binom{\frac{j}{\kappa}}{k} \frac{(-1)^{r-j+k}}{b^{\frac{j}{\kappa}}} \Gamma \left( \frac{1}{a} \left( \frac{j}{\kappa} - k \right) + 1, 1 \right). \end{split}$$

This ends the proof of Proposition 1.

Based on Proposition 1, the mean of *X* can be expanded as

$$\mu\_1^\* = \beta \varepsilon \sum\_{j=0}^1 \sum\_{k=0}^\infty \binom{1}{j} \binom{\frac{j}{a}}{k} \frac{(-1)^{1-j+k}}{b^{\frac{j}{a}}} \Gamma\left(\frac{1}{a}\left(\frac{j}{a} - k\right) + 1, 1\right).$$

and the moment of order 2 of *X* can be expressed as

$$
\mu\_2^\* = \beta^2 \varepsilon \sum\_{j=0}^2 \sum\_{k=0}^\infty \binom{2}{j} \binom{\frac{j}{\alpha}}{k} \frac{(-1)^{k-j}}{b^{\frac{j}{\alpha}}} \Gamma\left(\frac{1}{a}\left(\frac{j}{\alpha} - k\right) + 1, 1\right).
$$

From the above moments, we derive the variance of *X* by *V* = *μ*∗ <sup>2</sup> − (*μ*<sup>∗</sup> <sup>1</sup> )2. Several other moment measures can be expressed in a similar manner, including the dispersion index, coefficient of variation, moment skewness, and moment kurtosis. More details on the moment skewness and moment kurtosis will be provided later.

The two following points can be proven by following the lines of the proof of Proposition 1.

• The *r*th moment of *X* about the mean can be expressed as

$$\begin{aligned} \mu\_{\mathcal{I}} &= E\left[ (X - \mu\_1^\*)^r \right] \\ &= \beta^r e \sum\_{j=0}^r \sum\_{k=0}^\infty \binom{r}{j} \binom{\frac{j}{\alpha}}{k} \frac{(-1)^{r-j+k}}{b^{\frac{j}{\alpha}}} \left( 1 + \frac{\mu\_1^\*}{\beta} \right)^{r-j} \Gamma \left( \frac{1}{a} \left( \frac{j}{\alpha} - k \right) + 1, 1 \right). \end{aligned}$$

Based on it, the standard moment skewness measure is defined by *SK* <sup>=</sup> *<sup>μ</sup>*3/*<sup>V</sup>* <sup>3</sup> 2 , and the standard moment kurtosis measure is defined by *KU* = *μ*4/*V*2, among other moment measures.

• The *r*th unconditional moment of *X* at a certain *t* > 0 can be expanded as

$$\begin{split} \mu\_{r}(t) &= E[X^{r} \mid X \le t] = \frac{\beta^{r} e}{1 - e^{1 - \left(1 + b\left(\frac{\beta}{t + \beta}\right)^{-a}\right)^{a}}} \sum\_{j=0}^{r} \sum\_{k=0}^{\infty} \binom{r}{j} \left(\frac{\frac{j}{a}}{k}\right) \frac{(-1)^{r - j + k}}{b^{\frac{j}{a}}} \times \frac{1}{t} \\ & \left[ \Gamma\left(\frac{1}{a}\left(\frac{j}{a} - k\right) + 1, 1\right) - \Gamma\left(\frac{1}{a}\left(\frac{j}{a} - k\right) + 1, \left(1 + b\left(\frac{\beta}{t + \beta}\right)^{-a}\right)^{a}\right) \right]. \end{split}$$

It is immediate that lim*t*→<sup>∞</sup> *μr*(*t*) = *μr*. The unconditional moments are useful in the expression of various important functions, such as the mean residual life and reversed mean residual life functions. For more information on these functions, see [29].

#### **4. Maximum Likelihood Estimates of the Parameters**

We now consder the NHLx distribution as a statistical model, and we assume that the parameters *a*, *b*, *α*, and *β* are unknown. We aim to give some details on the maximum likelihood estimates (MLEs) of the parameters. First, let *n* be a positive integer, *X*1, *X*2, ... *Xn* be independent and identically distributed random variables drawn from the NHLx distribution, and *x*1, *x*2, ... , *xn* be corresponding observations. Then, provided that inf(*x*1, *x*2, ... , *xn*) ≥ −*β*, the likelihood function and log-likelihood functions are defined by

$$\begin{aligned} L(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n; \boldsymbol{\xi}) &= \prod\_{i=1}^n f\_{\text{NHLx}}(\mathbf{x}\_i; \boldsymbol{\xi}) \\ &= \left(\frac{ab\boldsymbol{\alpha}}{\beta}\right)^n \prod\_{i=1}^n \left(\frac{\beta}{\mathbf{x}\_i + \beta}\right)^{-a+1} \prod\_{i=1}^n \left(1 + b\left(\frac{\beta}{\mathbf{x}\_i + \beta}\right)^{-a}\right)^{a-1} e^{-\sum\_{i=1}^n \left(1 + b\left(\frac{\beta}{\mathbf{x}\_i + \beta}\right)^{-a}\right)^{a-1}} \end{aligned}$$

and

$$\begin{aligned} \ell(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n; \zeta) &= \log L(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n; \zeta) \\ \ell &= n \log \left( \frac{aba}{\beta} \right) + (1 - a) \sum\_{i=1}^n \log \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right) + (a - 1) \sum\_{i=1}^n \log \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right) \\ &+ n - \sum\_{i=1}^n \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right)^a, \end{aligned}$$

respectively. Then, the MLEs of the parameters *a*, *b*, *α*, and *β*, say *a*ˆ, ˆ *b*, *α*ˆ, and *β*ˆ, respectively, are defined by

$$\mathcal{J} = (\mathfrak{a}, \mathfrak{b}, \mathfrak{a}, \mathfrak{b}) = \mathfrak{argmax}\_{\mathbb{Q}} \ell(\mathfrak{x}\_1, \mathfrak{x}\_2, \dots, \mathfrak{x}\_n; \mathbb{Q}).$$

In the case where *β* is known and we have surely inf(*x*1, *x*2, ... , *xn*) ≥ −*β*, the MLEs of *a*, *b*, and *α* are the solution of the following equations: *<sup>∂</sup> ∂a* -(*x*1, *x*2, ... , *xn*; *ζ*) = 0, *∂ ∂b* -(*x*1, *x*2,..., *xn*; *ζ*) = 0 and *<sup>∂</sup> ∂α* -(*x*1, *x*2,..., *xn*; *ζ*) = 0, where

$$\begin{aligned} \frac{\partial}{\partial a} \ell(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n; \zeta) &= \frac{n}{a} + \sum\_{i=1}^n \log \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right) \\ &- \sum\_{i=1}^n \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right)^a \log \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right), \end{aligned}$$

$$\begin{aligned} \frac{\partial}{\partial b} \ell(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n; \xi) &= \frac{n}{b} + (a-1) \sum\_{i=0}^n \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right)^{-1} \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \\ &- a \sum\_{i=1}^n \left( 1 + b \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \right)^{a-1} \left( \frac{\beta}{\mathbf{x}\_i + \beta} \right)^{-a} \end{aligned}$$

and

$$\begin{split} &\frac{\partial}{\partial\alpha}\ell(\mathbf{x}\_{1},\mathbf{x}\_{2},\ldots,\mathbf{x}\_{n};\xi) = \frac{n}{\alpha} - \sum\_{i=1}^{n}\log\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right) \\ &+ (1-a)b\sum\_{i=1}^{n}\left(1+b\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right)^{-a}\right)^{-1}\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right)^{-a}\log\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right) \\ &+ ab\sum\_{i=1}^{n}\left(1+b\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right)^{-a}\right)^{a-1}\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right)^{-a}\log\left(\frac{\beta}{\mathbf{x}\_{i}+\beta}\right). \end{split}$$

The above expressions do not have closed-form solutions; hence, they are to be solved numerically by iterative methods. These numerical values can be easily obtained using specific tools in statistical software such as the R software, and the MLE of *β* is obtained by taking its first-order statistics, as in [14]. It is also possible to determine the values of the standard errors (SEs) of the MLEs. For more information, see [30].

Based on the MLEs, we define the estimated pdf of the NHLx distribution by *fNHLx*(*x*; ˆ *ζ*). Conceptually, the curve of this estimated function must be close to the shape of the histogram of the data, among other visual criteria.

#### **5. Simulation Study**

In this section, we perform 1000 Monte Carlo simulation studies for three different sets of parameters and each of the sample sizes of *n* ∈ {50, 100, 250, 500, 750, 1000}. By considering the order (*a*, *b*, *α*, *β*), these sets of parameters are Set I = (0.5, 1.5, 5, 0.5), Set II = (0.5, 1.5, 4, 0.75), and Set III = (1.5, 0.5, 4, 0.5). Table 1 shows the mean MLEs (MMLEs), biases and mean squared errors (MSEs) of the studies.

From Table 1, it can be observed that as the sample size increases, the biases and MSEs of the MLEs decrease, and with the increase in the sample sizes, the MMLEs are closer to the true parameter values. These results prove the accuracy of the considered parameter strategy estimation.

**Table 1.** Simulation results related to the MLEs of the NHLx model parameters.

