**1. Introduction**

Burr devised a dynamic family of probability distributions based on the Pearson differential equations. The Burr XII (BXII) and Burr III (BIII) distributions are widely used models from the system of Burr distributions. On the contrary, according to [1], the Burr X (BX) model has also gained much attention from applied statisticians along with the BXII and BIII models. The prime reason is that these densities exists in simpler forms and can yield a range of shapes to model a variety of scenarios in diverse scientific fields. The authors in [2] are of the view that the most adaptable of these three is BIII, especially in environmental, reliability, and survival sciences. The BIII distribution is also called the Dagum distribution in studies of income, wage, and wealth distribution [3]. In the actuarial literature, it is known as the inverse Burr distribution [4] and the kappa distribution in the meteorological literature [5]. As per [4], it is a prime case of the four-parameter generalised Beta-II distribution. In order to follow the ambit regarding the scope of this provision, we now shift our attention to the BIII distribution. For a random variable *X* defined on a positive real line, the cumulative distribution function (cdf) and probability density function (pdf) of two-parameter BIII distribution, respectively, are given below:

$$F(\mathbf{x}; \mathbf{c}, k) = \left(\mathbf{1} + \mathbf{x}^{-\mathbf{c}}\right)^{-k} \tag{1}$$

**Citation:** Jamal, F.; Abuzaid, A.H.; Tahir, M.H.; Nasir, M.A.; Khan, S.; Mashwani, W.K. New Modified Burr III Distribution, Properties and Applications. *Math. Comput. Appl.* **2021**, *26*, 82. https://doi.org/ 10.3390/mca26040082

Academic Editors: Nicholas Fantuzzi and Fernando Casas

Received: 10 November 2021 Accepted: 15 December 2021 Published: 20 December 2021

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

**Copyright:** © 2021 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/).

and

$$f(\mathbf{x}; \mathbf{c}, k) = \mathbf{c} \, k \, \mathbf{x}^{-c-1} \left( 1 + \mathbf{x}^{-c} \right)^{-k-1},\tag{2}$$

where *c*, *k* > 0 are the shape parameters.

The shape parameter plays a significant role in yielding the hazard rate of BIII distribution, which can be decreasing or unimodal. Thus, it cannot be used to model lifetime data with a bathtub-shaped hazard function, such as human mortality and deterioration modelling. For the last few decades, statisticians have been developing various extensions and modifications in Weibull distribution due to its simple functional form. The two-parameter flexible Weibull extension of [6] has a hazard function that can be increasing, decreasing, or bathtub shaped. Zhang and Xie [7] studied the characteristics and application of the truncated Weibull distribution, which has a bathtub-shaped hazard function. A three-parameter model, called exponentiated Weibull distribution, was introduced by [8]. Another three-parameter model is referred to as the extended Weibull distribution by [9]. Xie et al. [10] proposed a three-parameter modified Weibull extension with a bathtubshaped hazard function. A new modified Weibull distribution by the authors in [11] has been presented with increasing and a bathtub-shaped hazard function.

Various extensions of BIII distribution have been studied in the literature. In reference [12], the authors studied low-flow frequency analysis in hydrology with threeparameter-modified BIII distribution with supreme interest in the lower tail of a distribution. Çankaya et al. [13] extended the BIII model by adding a skew parameter with an epsilon skew extension approach. Modi and Gill [14] introduced the unit BIII model. Haq et al. [15] introduced the unit-modified BIII model. Ali et al. [16] re-parameterized BIII distribution and proposed the modified BIII (MBIII) distribution with the following cdf:

$$F(\mathbf{x}) = \left(1 + \mu \mathbf{x}^{-\mathbf{c}}\right)^{\frac{-k}{\mu}} \qquad \mathbf{x} > \mathbf{0},\tag{3}$$

where *c*, *k*, and *μ* are the shape parameters. The authors claimed that the newly structured model is a limiting case of generalized inverse Weibull, BIII, and log-logistic distribution. Still, the density of the improved model can only model positively skewed data, which greatly dented the proposition of the model in the first place. Other extensions are mostly based on the generalized families of distributions that sare complex in nature. Some of them are mentioned as: Beta Dagum by [17], Modified BIII by [18], Marshall Olkin BIII by [19], Gamma BIII by [20], and Gamma BIII by [21]. However, we feel that a flexible model with computationally simpler functional forms is still presently needed. Motivated by a lack of availability of literature related to the modified BIII distribution, we present a much more flexible new modification of BIII distribution. The cdf of the new, modified BIII (NMBIII) distribution is defined as

$$F(\mathbf{x}; \mathbf{c}, k, \lambda) = \left(\mathbf{1} + \mathbf{x}^{-\mathbf{c}} \mathbf{e}^{-\lambda \mathbf{x}}\right)^{-k} \qquad \mathbf{x} > \mathbf{0},\tag{4}$$

where the e−*λ<sup>x</sup>* is the additional factor, with *λ* as the rate parameter and *c*, *k* are power parameters of the baseline model.

It is worth mentioning that when we use the additional term to add flexibility in the model, we specifically refer to the ability of the proposed model to fit a diverse range of real life phenomena. Additionally, flexibility may also be associated with the instantaneous failure rate or hazard rate, and is more commonly known as risk function. By selecting precise values for the shape parameters, the hazard rate function of the NMBIII distribution can take on a variety of appealing shapes. Generally speaking, the classical models deal with normal extreme observations. A new modification of BIII distribution will also enable us to observe the tail behaviour of the distribution, which is skewed in nature. Further, the BIII distribution has a monotonic decreasing and unimodal hazard rate function, but due to its modification, NMBIII has monotonic, decreasing, increasing, unimodal, bathtub, and approximately constant hazard-rate shapes. Moreover, many standard distributions are nested models or limiting cases of the Burr system of distributions, which include the

Weibull, exponential, logistic, generalised logistic, Gompertz, normal, extreme value, and uniform distributions. The NMBIII distribution outperforms most of these competitive existing models. When *λ* = 0, NMBIII distribution reduces to BIII distribution. When *λ* = 0 and *k* = 1, then NMBIII distribution gives us log-logistic distribution. When *k* = 1, then NMBIII distribution gives us modified log-logistic distribution (new). When *c* = 0 and *k* = 1, the NMBIII distribution reduces to logistic distribution. When *c* = 1, it reduces to modified skew logistic distribution (new). When *c* = 0 and *λ* = 1, it reduces to generalized logistic distribution type I or Burr type II, or this type has also been called the "skewlogistic" distribution (see [22]). In a nutshell, with the proposed NMBIII, we seek and hope to attract applied researchers from all scientific community to utilize it in the significant modelling of real-life scenarios.

The article is structured as follows: In Section 2, we focus our attention on the idea behind the new modification. {In Section 3, we acquaint the readers with some of the structural properties including the linear expansion, moments, mode, moment-generating functions, order statistics, and stochastic ordering of NMBIII distribution. In Section 4, model parameters are estimated by maximum likelihood method, and the Fisher information matrix is derived. Section 5 gives the simulation method based on complete and incomplete samples (middle censored). In Section 6, three data sets on complete and middlecensored data sets have been employed to established the authenticity of the proposed model to the readers. Section 7 consists of the concluding remarks and discussions.

#### **2. The New Modified BIII Model**

The modified Weibull (MW) distribution (see [23] has the cumulative survival function that is the product of the Weibull cumulative hazard function *αx<sup>β</sup>* and e*λx*. Hence, the distribution function was found to be

$$F(\mathbf{x}) = \left(1 - \mathbf{e}^{-\alpha x^{\beta}} \mathbf{e}^{\lambda x}\right) \mathbf{y}$$

which was later generalized to exponentiated form by [24] using Lehmann alternative-I.

In the same vein, Equation (4) has been modified. The pdf corresponding to (4) is given as:

$$f(\mathbf{x}; \mathbf{c}, k, \lambda) = \frac{k\left(\lambda + \frac{\mathbf{c}}{\mathbf{x}}\right)}{\mathbf{x}^{\mathbf{c}} \mathbf{e}^{\lambda \cdot \mathbf{x}}} \left(1 + \mathbf{x}^{-\mathbf{c}} \mathbf{e}^{-\lambda \cdot \mathbf{x}}\right)^{-k-1}.\tag{5}$$

The corresponding survival and hazard functions of NMBIII are, respectively, given by:

$$S(\mathbf{x}; \mathbf{c}, k, \lambda) = 1 - \left(\mathbf{1} + \mathbf{x}^{-\mathbf{c}} \mathbf{e}^{-\lambda \cdot \mathbf{x}}\right)^{-k} \tag{6}$$

and

$$h(\mathbf{x}; \mathbf{c}, k, \lambda) = \frac{k\left(\lambda + \frac{\mathbf{c}}{\mathbf{x}}\right)}{\mathbf{x}^{\epsilon} \mathbf{e}^{\lambda \cdot \mathbf{x}}} \frac{\left(1 + \mathbf{x}^{-\epsilon} \mathbf{e}^{-\lambda \cdot \mathbf{x}}\right)^{-k-1}}{1 - \left(1 + \mathbf{x}^{-\epsilon} \mathbf{e}^{-\lambda \cdot \mathbf{x}}\right)^{-k}}.\tag{7}$$

If a new random variable *y* is defined as *y* = <sup>1</sup> *<sup>x</sup>* in Equation (4), then we obtain the following model, referred to as modified Burr XII distribution, with cdf and pdf, respectively, as under

$$G(y) = 1 - \left(1 + \frac{y^c}{\mathbf{e}^{\frac{1}{y}}}\right)^{-k} \tag{8}$$

and

$$g(y) = \frac{k\left(\varepsilon + \frac{1}{y}\right)}{\mathbf{e}^{\frac{1}{y}}} y^{\varepsilon - 1} \left(1 + \frac{y^{\varepsilon}}{\mathbf{e}^{\frac{1}{y}}}\right)^{-k - 1}.\tag{9}$$

As far as we can tell, Equations (4) and (8) are first modifications of BIII distribution and BXII distributions, respectively. Thus, the proposed distribution in (4) is more

flexible and has tractable tail properties than its parent BIII distribution as well as MBIII distributions. The shapes of pdf and hrf are presented in Figures 1 and 2, respectively.

**Figure 1.** Density function of NMBIII distribution.

**Figure 2.** Hazard function for NMBIII distribution.

Figure 1 represents the different shapes of the proposed model, i.e., bimodal, reversed-J, right skewed, approximate left-skewed, and symmetrical shapes for different parameter values. Figure 2 reflects the different shapes of hazard function, which are increasing, decreasing, bathtub, upside-down bathtub, and nearly constant for different parameter values. The proposed distribution is more flexible and tractable than its parent BIII distribution, as well as MBIII distributions (see in Table 1).


**Table 1.** Sub models of NMBIII distributions.

#### **3. Some Properties of NMBIII**

In this section, we will provide some significant properties of the NMBIII distribution such as *r*th moment, *s*th incomplete moment, moment generating function, skewness, kurtosis, mode, and order statistics.

#### *3.1. Useful Expansion*

The generalized binomial theorem or power series is given by:

$$(1+z)^{-b-1} = \sum\_{i=0}^{\infty} \binom{b+i}{i} (-1)^i z^i. \tag{10}$$

Using series expansion in (10), Equation (4) becomes

$$f(\mathbf{x}; \mathbf{c}, k, \lambda) = \sum\_{i=0}^{\infty} \binom{k+i}{i} (-1)^i \frac{k\left(\lambda + \frac{\mathbf{c}}{\mathbf{x}}\right)}{\mathbf{x}^{\epsilon(i+1)} \mathbf{e}^{\lambda \cdot \mathbf{x}(i+1)}}.\tag{11}$$

This expression can be used to obtain the following properties of the NMBIII distribution.

#### *3.2. Moments*

The *r*th moment of NMBIII distribution is given by:

$$\begin{split} m\_r' &= E(X') = \int\_0^\infty x^r f(x) \, dx \\ &= \sum\_{i=0}^\infty \binom{k+i}{i} \left(-1\right)^i \int\_0^\infty x^{r-\zeta(i+1)} \left(\lambda + \frac{c}{x}\right) \mathbf{e}^{-\lambda \left(i+1\right)x} \, dx \\ &= \lambda \sum\_{i=0}^\infty a\_i \int\_0^\infty x^{r-\zeta(i+1)} \mathbf{e}^{-\lambda \left(i+1\right)x} \, dx + c \sum\_{i=0}^\infty a\_i \int\_0^\infty x^{r-\zeta(i+1)-1} \mathbf{e}^{-\lambda \left(i+1\right)x} \, dx \\ &= \lambda \sum\_{i=0}^\infty a\_i \Gamma(r - c\zeta(i+1) - 1) \left[\frac{1}{\lambda(i+1)}\right]^{r-\zeta(i+1)-1} \\ &+ c \sum\_{i=0}^\infty a\_i \Gamma(r - c\zeta(i+1)) \left(\frac{1}{\lambda(i+1)}\right)^{r-\zeta(i+1)} \\ &= \lambda \sum\_{i=0}^\infty a\_i \frac{\Gamma(r - c\zeta(i+1) - 1)}{\left(\lambda(i+1)\right)^{r-\zeta(i+1)}} \left(\frac{1}{i+1} + c(r - c\zeta(i+1) - 1)\right), \end{split} \tag{12}$$
 
$$\text{where } a\_i = \left(\begin{array}{c} k + i \\ & i \end{array}\right) (-1)^i \text{ and } \Gamma(a) \, b^a = \int\_0^\infty x^{a-1} e^{-bx} \, dx \text{ is gamma function.}$$

**Remark 1.** *By submitting r* = 1 *in Equation (13), one can find mean of the NMBIII distribution.*

The *s*th incomplete moment of NMBIII distribution is

$$\begin{split} T\_s'(\mathbf{x}) &= \lambda \sum\_{i=0}^{\infty} a\_i \gamma \left( r - c(i+1) - 1, \frac{\mathbf{x}}{\lambda(i+1)} \right) \left( \frac{1}{\lambda(i+1)} \right)^{r - c(i+1) - 1} \\ &+ c \sum\_{i=0}^{\infty} a\_i \gamma \left( r - c(i+1), \frac{\mathbf{x}}{\lambda(i+1)} \right) \left( \frac{1}{\lambda(i+1)} \right)^{r - c(i+1)}. \end{split} \tag{13}$$

The application of incomplete moment refers to the mean deviations and Bonferroni and Lorenz curves. These curves are useful in economics reliability, demography, insurance, and medicine, to mention few.

## *3.3. Moment-Generating Function*

The moment-generating function of NMBIII distribution is given by:

$$\begin{split} \mathcal{M}u(t) &= E(\boldsymbol{d}^{T}) = \int\_{i=0}^{\infty} \boldsymbol{\varepsilon}^{t \cdot x} f(\boldsymbol{x}) \, d\boldsymbol{x} \\ &= \sum\_{i=0}^{\infty} \binom{k+i}{i} \left(-1\right)^{i} \int\_{i=0}^{\infty} \mathbf{x}^{-\varepsilon(i+1)} \left(\lambda + \frac{c}{\lambda}\right) \mathbf{e}^{\left(t - \lambda\left(i+1\right)\right) \cdot \mathbf{x}} \, d\boldsymbol{x} \\ &= \sum\_{i=0}^{\infty} a\_{i} \int\_{i=0}^{\infty} \mathbf{x}^{-\varepsilon(i+1)} \left(\lambda + \frac{c}{\lambda}\right) \mathbf{e}^{\left(t - \lambda\left(i+1\right)\right) \cdot \mathbf{x}} \, d\mathbf{x} \\ &= \sum\_{i=0}^{\infty} a\_{i} \left(\lambda \int\_{i}^{\infty} \mathbf{x}^{-\varepsilon(i+1)} \mathbf{e}^{\left(t - \lambda\left(i+1\right)\right) \cdot \mathbf{x}} \, d\mathbf{x} + c \int\_{0}^{\infty} \mathbf{x}^{-\varepsilon(i+1) - 1} \, \mathbf{e}^{\left(t - \lambda\left(i+1\right)\right) \cdot \mathbf{x}} \, d\mathbf{x}\right) \\ &= \sum\_{i=0}^{\infty} a\_{i} \left(\lambda \frac{\Gamma(1 - \varepsilon(i+1))}{\left(\lambda(i+1) - t\right)^{1 - \varepsilon(i+1)}} + c \frac{\Gamma(-\varepsilon(i+1))}{\left(\lambda(i+1) - t\right)^{-\varepsilon(i+1)}}\right) . \end{split} \tag{14}$$

The skewness and kurtosis of the NMBIII distribution can be obtained numerically by the following expression.

$$\alpha = \frac{m\_3' - 3\, m\_2'\, m\_1' + 2\, m\_1'}{\left\{m\_2' - (m\_2')^2\right\}^{3/2}} \tag{15}$$

and

$$\beta = \frac{m\_4' - 4\, m\_3' m\_1' + 6\, m\_2' \, (m\_1')^2 - 3\, (m\_1')^4}{\left\{m\_2' - (m\_2')^2\right\}^2},\tag{16}$$

where *m <sup>r</sup>* is the *r*th moment can be obtained form Equation (13).

**Remark 2.** *The mode of the NMBIII distribution can be obtained as follows: taking the* log *of Equation* (5)*, one obtains*

$$\log f(\mathbf{x}) = \log k + \log \left(\lambda + \frac{c}{\mathbf{x}}\right) - c \log \mathbf{x} - \lambda \mathbf{x} - (k+1) \log \left(1 + \mathbf{x}^{-c} e^{-\lambda \cdot \mathbf{x}}\right), \tag{17}$$

*Taking derivative with respect to x, we get*

$$\frac{d}{dx}\log f(x) = \frac{-\frac{1}{x^2}}{\lambda + \frac{c}{x}} - \frac{c}{x} - \lambda + (k+1)\frac{x^{-\varepsilon}e^{-\lambda x}\left(\lambda + \frac{c}{x}\right)}{1 + x^{-\varepsilon}e^{-\lambda x}},\tag{18}$$

*by setting the above expression equal to zero and solving for x, one can find the mode. The numerical values of the first four moments are given in Table 2.*

**Table 2.** The numerical values of the first four moments (*m <sup>r</sup>*,*r* = 1, 2, 3, 4), skewness (*α*) and kurtosis (*β*) of the NMBIII for some parameter values.


## *3.4. Order Statistics*

The density function *fi*:*n*(*x*) of the *i*-th order statistic, for *i* = 1, ... , *n*, from i.i.d. random variables *X*1,..., *X*<sup>2</sup> following MBIII distribution is simply given by:

$$F\_{i:n}(\mathbf{x}) = \frac{n!}{(i-1)!(n-i)!} \sum\_{j=0}^{n-i} \binom{n-i}{j} \frac{(-1)^j}{j+i} F(\mathbf{x})^{j+i}.\tag{19}$$

The corresponding pdf is

$$f\_{\rm i:n}(\mathbf{x}) = \frac{n!}{(i-1)!(n-i)!} \sum\_{j=0}^{n-i} \binom{n-i}{j} (-1)^j f(\mathbf{x}) F(\mathbf{x})^{j+i-1}.\tag{20}$$

Using the pdf and cdf of NMBIII in Equations (4) and (5), we obtain

$$F\_{\rm i:n}(\mathbf{x}) = \frac{n!}{(i-1)!(n-i)!} \sum\_{j=0}^{n-i} \binom{n-i}{j} \frac{(-1)^j}{j+i} \left[1 + \mathbf{x}^{-\mathbf{c}} \mathbf{e}^{-\mathbf{A} \cdot \mathbf{x}}\right]^{-k(j+i)}.\tag{21}$$

Using series expansion in (10), we obtain

$$F\_{\rm i:n}(\mathbf{x}) = \sum\_{j=0}^{n-i} b\_j \sum\_{l=0}^{\infty} \binom{k(j+i)+l}{l} (-1)^l \mathbf{x}^{-\varepsilon l} \mathbf{e}^{-\lambda l \ge \varepsilon},\tag{22}$$

where *bj* = *<sup>n</sup>*! (*i*−1)!(*n*−*i*)! *<sup>n</sup>* <sup>−</sup> *<sup>i</sup> j* (−1)*<sup>j</sup> <sup>j</sup>*+*<sup>i</sup>* . Similarly, following the above algebra, we have

$$f\_{\hat{i}:n}(\mathbf{x}) = \sum\_{j=0}^{n-i} a\_j \sum\_{l=0}^{\infty} \binom{j+i+l}{l} (-1)^l \mathbf{x}^{-\mathbf{c}\left(l+1\right)} \left(\lambda + \frac{\mathbf{c}}{\mathbf{x}}\right) \mathbf{e}^{-\lambda\left(l+1\right)\cdot\mathbf{x}},\tag{23}$$

where *aj* = *k <sup>n</sup>*! (*i*−1)!(*n*−*i*)! *<sup>n</sup>* <sup>−</sup> *<sup>i</sup> j* (−1)*<sup>j</sup>* .

#### *3.5. Stochastic Ordering*

The concept of stochastic ordering is frequently used to show the ordering mechanism in life-time distributions. For more details about stochastic ordering, see [26]. A random variable is said to be stochastically greater (*X* ≤*st Y*) than *Y* if *FX*(*x*) ≤ *FY*(*x*) for all *x*. In the similar way, *X* is said to be stochastically lower (*X* ≤*st Y*) than *Y* in the


The stochastic orders defined above are related to each other, as the following implications.

$$X \leq\_{\text{rhr}} Y \Leftarrow X \leq\_{\text{llr}} Y \Rightarrow X \leq\_{\text{hr}} Y \Rightarrow X \leq\_{\text{stl}} Y \Rightarrow X \leq\_{\text{nvrl}} Y. \tag{24}$$

Let *X*<sup>1</sup> ∼ *NMBIII*(*c*1, *k*1, *λ*1) and *X*<sup>2</sup> ∼ *NMBIII*(*c*2, *k*2, *λ*2). Then, according to the definition of likelihood ratio ordering *<sup>f</sup>*(*x*) *g*(*x*) ,

$$f(\mathbf{x}) = \frac{k\_1 \left(\lambda\_1 + \frac{c\_1}{\mathbf{x}}\right)}{\mathbf{x}^{c\_1} \mathbf{e}^{\lambda\_1 \cdot \mathbf{x}}} \left(1 + \mathbf{x}^{-c\_1} \mathbf{e}^{-\lambda\_1 \cdot \mathbf{x}}\right)^{-k\_1 - 1} \text{.}\tag{25}$$

$$g(\mathbf{x}) = \frac{k\_2 \left(\lambda\_2 + \frac{c\_2}{x}\right)}{x^{c\_2} \mathbf{e}^{\lambda\_2 x}} \left(1 + x^{-c\_2} \mathbf{e}^{-\lambda\_2 x}\right)^{-k\_2 - 1},\tag{26}$$

and

$$\frac{f(\mathbf{x})}{g(\mathbf{x})} = \frac{k\_1}{k\_2} \frac{\left(\lambda\_1 + \frac{c\_1}{x}\right)}{\left(\lambda\_2 + \frac{c\_2}{x}\right)} \frac{\mathbf{x}^{c\_1} \mathbf{e}^{\lambda\_1 x}}{\mathbf{x}^{c\_1} \mathbf{e}^{\lambda\_1 x}} \frac{\left(1 + \mathbf{x}^{-c\_1} \mathbf{e}^{-\lambda\_1 x}\right)^{-k\_1 - 1}}{\left(1 + \mathbf{x}^{-c\_2} \mathbf{e}^{-\lambda\_2 x}\right)^{-k\_2 - 1}}.\tag{27}$$

Taking log on both sides and taking the derivative with respect to *x*, we obtain

$$\begin{split} \frac{d}{dx} \left( \frac{f(\mathbf{x})}{g(\mathbf{x})} \right) &= \frac{c\_1}{\mathbf{x}\_i^2 \left( \lambda\_1 + \frac{c\_1}{\mathbf{x}} \right)} - \frac{c\_2}{\mathbf{x}\_i^2 \left( \lambda\_2 + \frac{c\_2}{\mathbf{x}} \right)} + \frac{c\_2 - c\_1}{\mathbf{x}\_i} + (\lambda\_2 - \lambda\_1) \\ &+ (k\_2 + 1) \frac{\mathbf{x}^{-\varepsilon\_2} \mathbf{e}^{\lambda\_2 \ge} \left( \lambda\_2 + \frac{c\_2}{\mathbf{x}} \right)}{1 + \mathbf{x}^{-\varepsilon\_2} \mathbf{e}^{\lambda\_2 \ge}} - (k\_1 + 1) \frac{\mathbf{x}^{-\varepsilon\_1} \mathbf{e}^{\lambda\_1 \ge} \left( \lambda\_1 + \frac{c\_1}{\mathbf{x}} \right)}{1 + \mathbf{x}^{-\varepsilon\_1} \mathbf{e}^{\lambda\_1 \ge}}, \end{split} \tag{28}$$

if *c*<sup>1</sup> = *c*<sup>2</sup> = *c* and *λ*<sup>1</sup> = *λ*<sup>2</sup> = *λ*, then *<sup>d</sup> d x f*(*x*) *<sup>g</sup>*(*x*) < 0 if (*k*<sup>2</sup> < *k*1) and then *X* <*lr Y*.

#### **4. Maximum Likelihood Estimation**

In this section, we will use the maximum-likelihood method to estimate the unknown parameters of the proposed model from complete samples only. Let *x*1, *x*2, ... , *xn* be a random sample of size n from the NMBIII family given in Equation (4) distribution. The log-likelihood function for the vector of parameter Θ = (*c*, *k*, *λ*)*<sup>T</sup>* can be expressed as

$$\begin{aligned} l(\Theta) &= n \log k - c \sum\_{i=1}^{n} \log x\_i - \lambda \sum\_{i=1}^{n} x\_i + \sum\_{i=1}^{n} \log \left( \lambda + \frac{c}{x} \right) \\ &- (k+1) \sum\_{i=1}^{n} \log \left[ 1 + x^{-c} e^{-\lambda \cdot x} \right] \end{aligned}$$

Taking the derivative with respect to *λ*, *c*, *k*, respectively, we get

$$\begin{split} \mathcal{U}\_{k} &= \frac{\partial \operatorname{l}(\boldsymbol{\Theta})}{\partial \boldsymbol{k}} = \frac{n}{k} - \sum\_{i=1}^{n} \log \Big( 1 + \mathbf{x}^{-\boldsymbol{c}} \mathbf{e}^{-\boldsymbol{\lambda} \cdot \mathbf{x}} \Big) \\ \mathcal{U}\_{\boldsymbol{k}} &= \frac{\partial \operatorname{l}(\boldsymbol{\Theta})}{\partial \boldsymbol{\lambda}} = -\sum\_{i=1}^{n} \boldsymbol{x}\_{i} + \sum\_{i=1}^{n} \left( \boldsymbol{\lambda} + \frac{\boldsymbol{c}}{\boldsymbol{x}} \right)^{-1} + (k+1) \sum\_{i=1}^{n} \left( \frac{\mathbf{x}^{-\boldsymbol{c}} \mathbf{e}^{-\boldsymbol{\lambda} \cdot \mathbf{x}} \boldsymbol{x}\_{i}}{1 + \mathbf{x}^{-\boldsymbol{c}} \mathbf{e}^{-\boldsymbol{\lambda} \cdot \mathbf{x}}} \right) \\ \mathcal{U}\_{\boldsymbol{c}} &= \frac{\partial \operatorname{l}(\boldsymbol{\Theta})}{\partial \boldsymbol{c}} = -\sum\_{i=1}^{n} \log \operatorname{x}\_{i} + \sum\_{i=1}^{n} \left( \frac{1}{\mathbf{x}\_{i} \left( \boldsymbol{\lambda} + \frac{\boldsymbol{c}}{\boldsymbol{x}} \right)} \right) + (k+1) \sum\_{i=1}^{n} \left( \frac{\mathbf{x}^{-\boldsymbol{c}} \mathbf{e}^{-\boldsymbol{\lambda} \cdot \mathbf{x}} \log \operatorname{x}\_{i}}{1 + \mathbf{x}^{-\boldsymbol{c}} \mathbf{e}^{-\boldsymbol{\lambda} \cdot \mathbf{x}}} \right) \end{split}$$

Setting *Uk*, *Uλ*, and *Uk* equal zero and solving these equations simultaneously yields the maximum likelihood estimates.

The observed information matrix for the parameter vector is given by

$$
\begin{pmatrix}
\mathcal{U}\_{kk} & \mathcal{U}\_{K\lambda} & \mathcal{U}\_{k\mathcal{L}} \\
\end{pmatrix},
$$

whose elements are given below

*Uk k* <sup>=</sup> <sup>−</sup> *<sup>n</sup> k*2 *Uk <sup>λ</sup>* = *n* ∑ *i*=1 *<sup>x</sup>*−*c*−<sup>1</sup> *<sup>i</sup>* <sup>e</sup>*<sup>λ</sup> xi* 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi Uk c* = − *n* ∑ *i*=1 *x*−*<sup>c</sup> <sup>i</sup>* <sup>e</sup>*<sup>λ</sup> xi* log *xi* 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi U<sup>λ</sup> <sup>c</sup>* = − *n* ∑ *i*=1 1 *xi λ* + *<sup>c</sup> xi* <sup>2</sup> + (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) *n* ∑ *i*=1 *x*1−2*<sup>c</sup>* e−<sup>2</sup> *<sup>λ</sup> xi* log *xi* 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi* <sup>2</sup> <sup>+</sup> <sup>e</sup>−*<sup>λ</sup> xi <sup>x</sup>*1−*<sup>c</sup> <sup>i</sup>* log *xi* 1 + *x*−*<sup>c</sup> e*−*<sup>λ</sup> xi Uλ λ* = − *n* ∑ *i*=1 1 *λ* + *<sup>c</sup> xi* <sup>2</sup> <sup>−</sup> (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) *n* ∑ *i*=1 *x*2−2*<sup>c</sup>* e−2˘xi 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi* <sup>2</sup> <sup>+</sup> *<sup>e</sup>*−*<sup>λ</sup> xi <sup>x</sup>*2−*<sup>c</sup> i* 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi Uc c* = − *n* ∑ *i*=1 1 *x*2 *i λ* + *<sup>c</sup> xi* <sup>2</sup> + (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) *n* ∑ *i*=1 *x*−2*<sup>c</sup>* e−<sup>2</sup> *<sup>λ</sup> xi*(log *xi*)<sup>2</sup> 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi* <sup>2</sup> <sup>+</sup> <sup>e</sup>−*<sup>λ</sup> xi <sup>x</sup>*−*<sup>c</sup> <sup>i</sup>* log *xi* 1 + *x*−*<sup>c</sup>* e−*<sup>λ</sup> xi* 

#### **5. Middle-Censoring**

The middle-censoring scheme is a non-parametric general censoring mechanism proposed by [27], where other censoring schemes can be obtained as special cases of this middle-censoring scheme (see [28]).

For *n* identical lifetimes *T*1, ... , *Tn* with a random censoring interval (*Li* ≤ *Ri*) at the *i*th item with some unknown bivariate distribution. Then, the exact value of *Ti* is observable only if *Ti* ∈/ [*Li* ≤ *Ri*]; otherwise, the interval (*Li* ≤ *Ri*) is observed.

Middle-censoring had previously been applied to exponential and Burr XII lifetime distributions (see [28,29]). Furthermore, it was extended to parametric models with covariates [30], and its robustness was investigated by [31].

In this section, we analyse the NMBIII lifetime data when they are middle-censored. Assume that *T*1, ... , *Tn* are *i.i.d.* NMBIII (*c*, *λ*, *k*) random variable and let *Zi* = *Ri* − *Li*, *i* = 1, ... , *n* be another random variable that defines the length of the censoring interval with exponential distribution with mean *γ*−1, where the left-censoring point for each individual *Li* is assumed to also be an exponential random variable with mean *θ*−1. Moreover, the *T <sup>i</sup>s*, *L i s*, and *Z i s* are all independent of each other and the observed data, and *X i s* are given by *Xi* = *Ti if Ti* <sup>∈</sup>/ (*Li* <sup>≤</sup> *Ri*), (*Li* ≤ *Ri*) *otherwise*.
