**1. Introduction**

The introduction of new statistical distributions, defined on both the whole real line and the positive real line, is required for the interpretation of real-world occurrences. In both of these lines, a large number of probability distributions for real data sets have been presented recently. Several distributions that fall within these two categories have demonstrated their importance from both a theoretical and practical standpoint. It is true, however, that in many practical disciplines such as economics, medicine, biology, and others, there is a problem of bounded phenomena, or uncertainty resulting from factors such as rates, indices, proportions, and test scores. It is also important to note that, in comparing distributions with unbounded support to those with bounded support, one can find notable scarceness. As a result, in the recent past, some authors have focused on developing distributions that are defined on the bounded interval using any one of the parent distribution transformation techniques. For better understanding, we refer to recent papers in [1,2].

The beta distribution is definitely the most famous distribution among distributions defined in the (0, 1) interval. While the beta distribution is useful for modeling data on the unit interval, many other distributions have been proposed and studied over time. The readers may refer to the Topp–Leone distribution (see [3]), Kumaraswamy distribution (see [4]) and transformed Leipnik distribution (see [5]). However, in recent times, there has been a growing interest by statisticians in proposing distributions defined by the unit interval corresponding to any continuous distribution. The readers may refer to the log-Lindley distribution (see [6]), exponentiated Topp–Leone distribution (see [7]), unit inverse Gaussian distribution (see [8]), unit Lindley distribution (see [9]), unit Weibull distribution

**Citation:** Krishna, A.; Maya, R.; Chesneau, C.; Irshad, M.R. The Unit Teissier Distribution and Its Applications. *Math. Comput. Appl.* **2022**, *27*, 12. https://doi.org/ 10.3390/mca27010012

Academic Editors: Sandra Ferreira and Nicholas Fantuzzi

Received: 29 December 2021 Accepted: 29 January 2022 Published: 1 February 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/).

(see [10]), unit Burr-III distribution (see [11]), and unit half normal distribution (see [12]), among others. As a result, a lot of contributions are made to the existing literature. This recent leap in the number of research papers devoted to proposing new distributions in the unit interval exhibits their growing relevance. Despite the fact that many distributions have been proposed and studied as alternatives, there is still no agreement on which distribution is preferred.

In this article, a new bounded distribution is introduced by considering the baseline distribution as the Teissier distribution (*TD*) which was first introduced by the French biologist Georges Teissier in modeling the mortality of domestic animal species as a result of pure aging (see [13]). A continuous random variable *Y* is said to have the *TD* with parameter *θ* > 0 if its probability density function (pdf) is given by

$$f\_Y(y; \theta) = \theta \left(\epsilon^{\theta y} - 1\right) \exp\left(\theta y - \epsilon^{\theta y} + 1\right), \qquad y > 0. \tag{1}$$

The cumulative distribution function (cdf) of *Y* is given by

$$F\_Y(y; \theta) = 1 - \exp\left(\theta y - e^{\theta y} + 1\right), \qquad y > 0.1$$

Despite its importance, the *TD* has been overlooked in the statistical literature. Ref. [14] examined the location version of this model and explored a characterization based on life expectancy as part of a demographic study, four decades after it was introduced. In a later study, Muth [15] found that a *TD*-derived distribution has a heavier tail than wellknown lifetime distributions such as the gamma, lognormal, and Weibull distributions. This distribution was later used by [16] to estimate the lifetime distribution for a German dataset based on used car prices. The *TD* and its location version have been forgotten since [16], and no further references appear in the literature until 2008. Ref. [17] showed that as the parameter decreases to zero, the model approaches the classical exponential distribution. The authors of [18] reintroduced this distribution and its scaled version, and studied its important properties using a generalized integro-exponential function [18,19] can be considered as the two-parameter extensions of the *TD*.

The study of the *TD* has gained momentum in both theoretical and applied perspectives after the work of [19–23] substantiate this claim. Meanwhile, the *TD* and its various variants proved their signature compared to other models for various real-life data sets. This is the chief motivation for looking into the *TD* defined on the unit interval, so-called the unit Teissier distribution (*UTD*), and its competency compared to other popular distributions defined on the unit interval. Finally, it can be concluded that it performs well compared to other distributions for the real-life data sets considered here.

As in the case of real-life data sets that vary over positive real lines and possess a bathtub-shaped hazard rate function (hrf), some real-life data sets that are defined on the unit interval also possess a bathtub-shaped hrf. In the former case, there is a tremendous amount of work in the literature on various models, such as [24–26]. As far as the latter case is considered, the papers in [8,27] showed that the unit inverse Gaussian and logit slash distributions, respectively, provide a bathtub-shaped hrf. These two distributions are more intricate because of the presence of the log function in their pdfs, and so their cdfs are not obtained in closed form. Hence, simulation experiments from these two distributions are very difficult. For some real data sets, the logit slash distribution does not yield a numerical estimation of parameters. Moreover, the hrfs of both of these distributions are not analytically tractable, even though they sketch the graph of the hrf and provide its various shapes. In addition, Korkmaz [27] proved that the logit slash distribution possesses a N-shaped (modified bathtub shape) and w-shaped bathtub shaped hrf. Another distribution defined on the unit interval, which also possesses a non-monotone shaped hrf, is the unit Modified Burr III distribution (see [28]), which has three parameters. Models with a lower number of parameters are preferable for various reasons, such as ease of simulation drawing and less difficulty in inferential aspects.

The *UTD* solves these two problems by providing closed form expressions for pdf, cdf, moments, inequality measures, and entropy measures. Secondly, it possesses monotone (increasing and decreasing) as well as non-monotone functions, such as bathtub-shaped and N-shaped hrf (modified bathtub shaped) with a smaller number of parameters. Based on [29], N-shaped hrfs appear in mortality among breast cancer patients.

The rest of the paper is presented as follows. In Section 2, the *UTD* is introduced and the statistical properties such as moments, incomplete moments, probability-weighted moments, and quantile function are studied. The Shannon entropy and extropy are derived in Section 3. Discussion on different methods of estimation, such as maximum likelihood (ML), ordinary, weighted least square and Bayesian estimation are carried out in Section 4. In Section 5, a simulation study is performed to evaluate the performance of the model parameter estimates. In Section 6, the proposed distribution is elucidated with two real data sets. Finally, the conclusions are presented in Section 7.

## **2. The Unit Teissier Distribution**

In this Section, the *UTD* is presented and we determine some of its statistical properties.

## *2.1. Presentation*

Mathematically, the *UTD* is derived from the transformation *X* = *e*−*<sup>Y</sup>* by taking into account the definition of the pdf given by (1). Its definition is formalized below.

**Definition 1.** *A random variable X is said to follow the UTD with parameter θ* > 0*, if its pdf is of the following form:*

$$f(\mathbf{x}; \theta) = \theta \left(\mathbf{x}^{-\theta} - 1\right) \mathbf{x}^{-(\theta+1)} \mathbf{e}^{-\mathbf{x}^{-\theta} + 1}, \qquad \mathbf{x} \in (0, 1). \tag{2}$$

The corresponding cdf is given by

$$F(\mathbf{x}; \theta) = \mathbf{x}^{-\theta} e^{-\mathbf{x}^{-\theta} + 1}. \tag{3}$$

From (2) and (3), the survival function *S*(*x*; *θ*), hrf *h*(*x*; *θ*) and reversed hrf *τ*(*x*; *θ*) of the *UTD* are obtained as follows (for *x* ∈ (0, 1)):

$$S(\mathbf{x};\theta) = 1 - \mathbf{x}^{-\theta}e^{-\mathbf{x}^{-\theta} + 1},$$

$$h(\mathbf{x};\theta) = \frac{f(\mathbf{x};\theta)}{S(\mathbf{x};\theta)} = \frac{\theta\left(\mathbf{x}^{-\theta} - 1\right)\mathbf{x}^{-(\theta+1)}e^{-\mathbf{x}^{-\theta} + 1}}{1 - \mathbf{x}^{-\theta}e^{-\mathbf{x}^{-\theta} + 1}}\tag{4}$$

and

$$\pi(\mathfrak{x};\theta) = \frac{f(\mathfrak{x};\theta)}{F(\mathfrak{x};\theta)} = \frac{\theta\left(\mathfrak{x}^{-\theta} - 1\right)}{\mathfrak{x}}.$$

### *2.2. Shapes of the pdf and hrf*

As immediate shape properties for the pdf and hrf, we can say that

$$\lim\_{x \to 0} f(x; \theta) = \lim\_{x \to 0} h(x; \theta) = 0,$$

and that the *UTD* is unimodal; the mode is the value of *x* maximizing *f*(*x*; *θ*), which is given by

$$x\_{modc} = 2^{-\frac{1}{\theta}} \left( \frac{1 + 3\theta - \sqrt{5\theta^2 + 2\theta + 1}}{1 + \theta} \right)^{\frac{1}{\theta}}.1$$

In order to obtain an overview of the shapes of the pdf (2) and hrf (4), the corresponding plots for different choices of the parameter *θ* are given in Figures 1 and 2, respectively.

**Figure 1.** Plots of various shapes of the pdf of the *UTD* for varying values of the parameter *θ*.

**Figure 2.** *Cont*.

**Figure 2.** Plots of various shapes of the hrf of the *UTD* for varying values of the parameter *θ*.

The graphs of the hrf for various combinations of parameters show various shapes, including increasing, decreasing, bathtub, and N-shaped (modified bathtub shape). According to [29], mortality among breast cancer patients is characterized by an N-shaped hrf. This is one of the prominent properties of the *UTD*.

## *2.3. Moments*

Analytically, the *UTD*'s non-central moments can be expressed in terms of the upper incomplete gamma function. This is specified in the result below.

**Result 1.** *For any non-negative integer r, the rth non-central moment of a random variable X with the UTD is*

$$\mu\_r' = E(X^r) = e\{\Gamma\left(-\frac{r}{\theta} + 2, 1\right) - \Gamma\left(-\frac{r}{\theta} + 1, 1\right)\},$$

*where E denotes the expectation operator, e* ≈ 2.718282 *(the exponential function taken at the value* 1*) and* Γ(*a*, *b*) = <sup>∞</sup> *<sup>b</sup> <sup>t</sup>a*−1*e*−*<sup>t</sup> dt is the upper incomplete gamma function.*

**Proof.** Based on (2) and the definition of non-central moment, as well as the change of variable *u* = *x*−*θ*, and considering the definition of the upper incomplete gamma function, we have

$$\begin{split} \mu'\_r &= \int\_0^1 \mathbf{x}^r f(\mathbf{x}; \boldsymbol{\theta}) d\mathbf{x} = \int\_1^\infty \boldsymbol{u}^{-\frac{r}{\theta}} (\boldsymbol{u} - 1) \boldsymbol{e}^{-\boldsymbol{u} + 1} d\boldsymbol{u} \\ &= \varepsilon \{ \boldsymbol{\Gamma} \left( -\frac{r}{\theta} + 2, 1 \right) - \boldsymbol{\Gamma} \left( -\frac{r}{\theta} + 1, 1 \right) \}. \end{split}$$

This ends the proof.

The statistical measures such as the mean (*μ*) and variance (*σ*2) can be calculated numerically by using R software. Figure 3a,b illustrate the behavior of the mean and variance for varying values of the parameter *θ*, respectively.

**Figure 3.** Plots of the (**a**) mean and (**b**) variance of the *UTD* for varying values of the parameter *θ*.

#### *2.4. Incomplete Moments*

As a generalized version of the non-central moments, the *UTD*'s incomplete moments can be also expressed in terms of the upper incomplete gamma function.

**Result 2.** *For any non-negative integer <sup>k</sup> and <sup>x</sup>* ∈ (0, 1)*, the <sup>k</sup>th incomplete moment at <sup>x</sup> of a random variable X with the UTD is*

$$m\_k(\mathbf{x}; \theta) = E(\mathbf{X}^k \mathbf{1}\_{X \le \mathbf{x}}) = \mathfrak{e}\left\{ \Gamma\left(-\frac{k}{\theta} + 2, \mathbf{x}^{-\theta}\right) - \Gamma\left(-\frac{k}{\theta} + 1, \mathbf{x}^{-\theta}\right) \right\}.$$

**Proof.** By the definition of the *k*th incomplete moment, we have

$$m\_k(x; \theta) = \int\_0^\chi t^k f(t; \theta) dt.$$

In the above integral expression, after making the change of variable *u* = *t* <sup>−</sup>*θ*, the proof is similar to the one of Result 1, so the details are excluded here.

There is an interesting aspect of the fact that the first incomplete moment *m*1(*x*; *θ*) can be used to compute the mean deviation from the mean (*μ*) given by *E*[|*X* − *μ*|] = 2*μF*(*μ*; *θ*) − 2*m*1(*μ*; *θ*).

#### *2.5. Probability-Weighted Moments*

The probability-weighted moments of the *UTD* are now under investigation.

**Result 3.** *For any non-negative integers r and s, the* (*r*,*s*)*th probability-weighted moment of a random variable X with the UTD is given as*

$$\begin{aligned} m\_{r,s}(\mathbf{x};\theta) &= E(X^r F(X;\theta)^s) \\ &= \frac{\mathfrak{c}^{s+1}}{(s+1)^{-\frac{s}{3}+s+2}} \Big\{ \Gamma\left(-\frac{r}{\theta}+s+2, s+1\right) - (s+1)\Gamma\left(-\frac{r}{\theta}+s+1, s+1\right) \Big\}. \end{aligned}$$

**Proof.** By the definition of the (*r*,*s*)*th* probability-weighted moment and by making the change of variable *u* = (*s* + 1)*x*−*θ*, we obtain

$$\begin{split} m\_{r,s}(\mathbf{x};\theta) &= \int\_0^1 \mathbf{x}^r (F(\mathbf{x};\theta))^s f(\mathbf{x};\theta) d\mathbf{x} \\ &= \epsilon^{s+1} \theta \int\_0^1 \mathbf{x}^{r-\theta(s+1)} \mathbf{x}^{-(\theta+1)} \mathbf{e}^{-(s+1)} \mathbf{x}^{-\theta} d\mathbf{x} - \epsilon^{s+1} \theta \int\_0^1 \mathbf{x}^{r-\theta s} \mathbf{x}^{-(\theta+1)} \boldsymbol{\epsilon}^{-(s+1)} \mathbf{x}^{-\theta} d\mathbf{x} \\ &= \frac{\epsilon^{s+1}}{(s+1)^{-\frac{s}{\theta}+s+2}} \int\_{s+1}^\infty u^{-\frac{r}{\theta}+s+1} \mathbf{e}^{-u} du - \frac{\epsilon^{s+1}}{(s+1)^{-\frac{s}{\theta}+s+1}} \int\_{s+1}^\infty u^{-\frac{r}{\theta}+s} \boldsymbol{\epsilon}^{-u} d\mathbf{u} \\ &= \frac{\epsilon^{s+1}}{(s+1)^{-\frac{s}{\theta}+s+2}} \Big\{ \Gamma\left(-\frac{r}{\theta}+s+2, s+1\right) - (s+1)\Gamma\left(-\frac{r}{\theta}+s+1, s+1\right) \Big\}. \end{split}$$

The desired result is obtained.

#### *2.6. Mean Residual Life Function*

**Result 4.** *For any t* ∈ (0, 1)*, the mean residual life function at t of a random variable X with the UTD is*

$$\begin{split} r(t; \theta) &= E(X - t | X > t) \\ &= \frac{1}{1 - t^{-\theta} e^{-t^{-\theta} + 1}} \left\{ 1 - t - \frac{c}{\theta} \left[ \Gamma \left( 1 - \frac{1}{\theta}, 1 \right) - \Gamma \left( 1 - \frac{1}{\theta}, t^{-\theta} \right) \right] \right\}. \end{split}$$

**Proof.** By the definition of the mean residual life function and by making the change of variable *u* = *x*−*θ*, we obtain

$$\begin{split} r(t;\theta) &= \frac{1}{S(t;\theta)} \int\_{t}^{1} S(x;\theta) dx \\ &= \frac{1}{S(t;\theta)} \Big\{ 1 - t - \int\_{t}^{1} x^{-\theta} e^{-x^{-\theta} + 1} dx \Big\} \\ &= \frac{1}{S(t;\theta)} \Big\{ 1 - t - \frac{\varepsilon}{\theta} \int\_{1}^{t^{-\theta}} u^{-\frac{1}{\theta}} e^{-u} du \Big\} \\ &= \frac{1}{1 - t^{-\theta} e^{-t^{-\theta} + 1}} \Big\{ 1 - t - \frac{\varepsilon}{\theta} \Big[ \Gamma \Big( 1 - \frac{1}{\theta}, 1 \Big) - \Gamma \Big( 1 - \frac{1}{\theta}, t^{-\theta} \Big) \Big] \Big\}. \end{split}$$

The claimed result is achieved.

## *2.7. Quantile Function*

In addition to its remarkable property of being expressed in closed form, the *UTD* is also capable of representing the quantile function in terms of the negative branch of the Lambert W function. We recall that the Lambert W function is a multivalued complex function defined as the solution of the equation *W*(*z*)*eW*(*z*) = *z*, where *z* is a complex number. For any real numbers *z* ≥ −1/*e*, *W*(*z*) has two real branches. The real branch taking on values in [−1, ∞) is called the principal branch and is denoted by *W*0, and the one taking on values in (−∞, −1] is called the negative branch and is denoted by *W*−1. For a comprehensive review of this special function, readers are referred to [30]. From a computational perspective, the Lambert *W* function is available in computer algebra systems such as Maple (function LambertW), Mathematica (function ProductLog), and Matlab (function lambertW), and also in programming languages such as R [31] (functions lambert\_W0 and lambert\_Wm1 for *W*<sup>0</sup> and *W*−1, respectively, in the package gsl).

**Result 5.** *The quantile function of the UTD is given by*

$$Q(p; \theta) = \left[ -\mathcal{W}\_{-1} \left( -\frac{p}{\mathcal{e}} \right) \right]^{-\frac{1}{\theta}},\tag{5}$$

*where p* ∈ (0, 1) *and W*−<sup>1</sup> *denotes the negative branch of the Lambert W function.*

**Proof.** The cdf *F*(*x*; *θ*) should be inverted to obtain the quantile function of *UTD*. Thus, we need to solve *F*(*x*; *θ*) = *p* according to *x*, so

$$-\mathbf{x}^{-\theta}\mathbf{e}^{-\mathbf{x}^{-\theta}} = -\frac{p}{\varepsilon} \iff -\mathbf{x}^{-\theta} = \mathcal{W}\_{-1}\left(-\frac{p}{\varepsilon}\right) \Leftrightarrow \mathbf{x} = \left[-\mathcal{W}\_{-1}\left(-\frac{p}{\varepsilon}\right)\right]^{-\frac{1}{\theta}}.$$

The obtained *x* corresponds to the quantile function with respect to *p*. The desired result is achieved.

The quantile function in (5) can be used to obtain the following quantities:


$$S = \frac{Q(0.25; \theta) + Q(0.75; \theta) - 2M}{Q(0.75; \theta) - Q(0.25; \theta)};$$

• the Moors coefficient of kurtosis (with correction) defined by

$$T = \frac{Q(0.875; \theta) - Q(0.625; \theta) + Q(0.375; \theta) - Q(0.125; \theta)}{Q(0.75; \theta) - Q(0.25; \theta)} - 1.23.$$

#### **3. Shannon Entropy and Extropy**

This section discusses the Shannon entropy and extropy of the *UTD*.

#### *3.1. Shannon Entropy*

Conceptually, the Shannon entropy is the amount of information into a random variable. For a random variable *T* with pdf *f*(*t*) with support [0, 1], its Shannon entropy is defined as

$$H(T) = -\int\_0^1 f(t) \log f(t) dt.\tag{6}$$

The Shannon entropy for the *UTD* is obtained in terms of the following quantities:


The next result is about the expression of the Shannon entropy.

**Result 6.** *The Shannon entropy of a random variable T with the UTD has the following form*

$$\begin{split} H(T) &= - (\log \theta) (\epsilon \Gamma(2, 1) - 1) - \Gamma^1(2) - \frac{\epsilon(\theta + 1)}{\theta} \left( \Gamma^1(2, 1) - E\_1(1) \right) \\ &- 2\epsilon \Gamma(2, 1) + \epsilon \Gamma(3, 1) + 1. \end{split} \tag{7}$$

**Proof.** By considering the change of variables, *u* = *t* <sup>−</sup>*<sup>θ</sup>* and *v* = *t* <sup>−</sup>*<sup>θ</sup>* <sup>−</sup> 1 and taking into account (2), we obtain

$$\begin{split} \int\_{0}^{1} f(t) \log f(t) dt &= \epsilon (\log \theta) \left( \int\_{1}^{\infty} u e^{-u} du - \int\_{1}^{\infty} e^{-u} du \right) + \int\_{0}^{\infty} v e^{-v} (\log v) dv \\ &+ \frac{\epsilon (\theta + 1)}{\theta} \left( \int\_{1}^{\infty} u \log u e^{-u} du - \int\_{1}^{\infty} e^{-u} (\log u) du \right) \\ &+ \epsilon \int\_{1}^{\infty} (2u - u^{2} - 1) e^{-u} du. \end{split} \tag{8}$$

Substituting (8) in (6) and by using the quantities listed, (7) is obtained.

#### *3.2. Extropy*

Extropy is a complementary dual of Shannon entropy that has a variety of interesting applications, including appropriate scoring of forecasting distributions, comparing the uncertainty of two random variables, and astronomical measurements of heat distributions in galaxies, among others. For a random variable *T* with pdf *f*(*t*) with support [0, 1], its extropy is defined as

$$J(T) = -\frac{1}{2} \int\_0^1 f^2(t)dt. \tag{9}$$

**Result 7.** *The extropy of a random variable T with the UTD has the following form*

$$J(T) = \frac{\theta \epsilon^2}{2^{\frac{1}{\theta} + 3}} \left\{ \Gamma \left( \frac{1}{\theta} + 3, 2 \right) - \frac{1}{4} \Gamma \left( \frac{1}{\theta} + 4, 2 \right) - \Gamma \left( \frac{1}{\theta} + 2, 2 \right) \right\}. \tag{10}$$

**Proof.** By considering the change of variable *u* = 2*t* <sup>−</sup>*<sup>θ</sup>* and taking into account (2), the numerator of (9) can be obtained as

$$\begin{split} \int\_{0}^{1} f^{2}(t;\theta)dt &= \theta^{2}\varepsilon^{2} \int\_{0}^{1} t^{-2(\theta+1)}t^{-2\theta}\varepsilon^{-2t^{-\theta}}dt + \theta^{2}\varepsilon^{2} \int\_{0}^{1} t^{-2(\theta+1)}\varepsilon^{-2t^{-\theta}}dt \\ &- \theta^{2}\varepsilon^{2} \int\_{0}^{1} t^{-2(\theta+1)}2t^{-\theta}\varepsilon^{-2t^{-\theta}}dt \\ &= \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+4}} \int\_{2}^{\infty} u^{\frac{1}{2}+3}\varepsilon^{-u}du + \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+2}} \int\_{2}^{\infty} u^{\frac{1}{2}+1}\varepsilon^{-u}du - \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+2}} \int\_{2}^{\infty} u^{\frac{1}{2}+2}\varepsilon^{-u}du \\ &= \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+4}}\Gamma\left(\frac{1}{\theta}+4,2\right) + \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+2}}\Gamma\left(\frac{1}{\theta}+2,2\right) - \frac{\theta\varepsilon^{2}}{2^{\frac{1}{2}+2}}\Gamma\left(\frac{1}{\theta}+3,2\right). \end{split} \tag{11}$$

Substituting (11) in (9), (10) is obtained.

#### **4. Estimation and Inference**

*4.1. Maximum Likelihood Estimation*

Let *x*1, *x*2, ... , *xn* be a random sample of values of size *n* from a random variable *X* with the *UTD* of unknown parameter *θ*. Then the likelihood function of *θ* is given by

$$I(\theta) = \theta^n \prod\_{i=1}^n \mathbf{x}\_i^{-(\theta+1)} \prod\_{i=1}^n \left(\mathbf{x}\_i^{-\theta} - 1\right) e^{-\sum\_{i=1}^n \left(\mathbf{x}\_i^{-\theta} - 1\right)}.\tag{12}$$

Then the ML estimate (MLE) of *θ*, say ˆ *θ*, is obtained by maximizing *l*(*θ*) with respect to *<sup>θ</sup>*. Thus, for any *<sup>θ</sup>* <sup>&</sup>gt; 0, we have *<sup>l</sup>*(*θ*) <sup>≤</sup> *<sup>l</sup>*(<sup>ˆ</sup> *θ*). Practically, one can use the derivative technique; the partial derivative of log *l*(*θ*) with respect to the parameter is

$$\frac{\partial \log l(\theta)}{\partial \theta} = \frac{n}{\theta} - \sum\_{l=1}^{n} \frac{\boldsymbol{\chi}\_{l}^{-\theta} \log \boldsymbol{\chi}\_{l}}{\boldsymbol{\chi}\_{l}^{-\theta} - 1} - \sum\_{l=1}^{n} \log \boldsymbol{\chi}\_{l} + \sum\_{l=1}^{n} \boldsymbol{\chi}\_{l}^{-\theta} \log \boldsymbol{\chi}\_{l}.$$

and ˆ *θ* is obtained by solving the equation *<sup>∂</sup>* log *<sup>l</sup>*(*θ*) *∂θ* = 0 according to *θ*. Numerical optimization approaches employing mathematical tools such as R and Mathematica are the only way to achieve this.

#### Fisher Information Matrix and Asymptotic Confidence Interval

In order to carry out statistical inference on the parameters of the *UTD*, the 1 × 1 expected Fisher information matrix is needed. It can, however, be efficiently approximated by the observed Fisher information matrix *J*(ˆ *θ*) given by

$$J(\theta) = -\left. \frac{\partial^2 \log l(\theta)}{\partial \theta^2} \right|\_{\theta = \theta^\*}$$

and we can approximate variance of ˆ *θ* as

$$\operatorname{Var}(\hat{\theta}) \approx J^{-1}(\hat{\theta}).$$

Then the approximate 100(1 − *φ*)% two-sided normal confidence interval for *θ* is given by ˆ *θ* ± *Z<sup>φ</sup>* 2 & *Var*(ˆ *θ*), where *Z<sup>φ</sup>* 2 is the upper *<sup>φ</sup>* 2 *th* percentile of a standard normal distribution and *φ* is the significance level.

#### *4.2. Ordinary and Weighted Least-Squares Estimation*

The ordinary least square (LS) estimation and the weighted least square (WLS) estimation were proposed by [32] to estimate the parameters of the beta distribution. Let *x*1:*n*, *x*2:*n*,..., *xn*:*<sup>n</sup>* be the ordered values of *x*1, *x*2,..., *xn*. Let us set

$$R(\theta) = \sum\_{i=1}^{n} \left[ F(x\_{i:n}; \theta) - \frac{i}{n+1} \right]^2.$$

Then the LS estimate (LSE) of *θ*, say ˜ *θ*, is obtained by minimizing *R*(*θ*) with respect to *θ*. Thus, for any *θ* > 0, we have *R*(ˆ *θ*) ≤ *R*(*θ*). Practically, the LSE can also be obtained by solving the following equation:

$$\frac{\partial R(\theta)}{\partial \theta} = 2 \sum\_{i=1}^{n} \left[ F(\mathbf{x}\_{i:n}; \theta) - \frac{i}{n+1} \right] D(\mathbf{x}\_{i:n}; \theta) = 0\_{\prime \prime}$$

where

$$D(\mathbf{x}; \theta) = \frac{\partial F(\mathbf{x}; \theta)}{\partial \theta} = -\mathfrak{e}^{-\mathbf{x}^{-\theta} + 1} \mathfrak{x}^{-2\theta} (\mathfrak{x}^{\theta} - 1) \log(\mathfrak{x}),$$

according to *θ*. Similarly, the WLS estimate (WLSE) of *θ*, say ˘ *θ*, is obtained by minimizing the non-linear function

$$\mathcal{W}(\theta) = \sum\_{i=1}^{n} \frac{(n+1)^2 (n+2)}{i(n-i+1)} \left[ F(x\_{i:n}; \theta) - \frac{i}{n+1} \right]^2,$$

and it can also be obtained by solving the following equation:

$$\frac{\partial \mathcal{W}(\theta)}{\partial \theta} = 2 \sum\_{i=1}^{n} \frac{(n+1)^2 (n+2)}{i(n-i+1)} \left[ F(\mathbf{x}\_{i:n}; \theta) - \frac{i}{n+1} \right] D(\mathbf{x}\_{i:n}; \theta) = 0.1$$

#### *4.3. Bayesian Estimation*

In this subsection, the estimate of the *UTD* parameter is calculated by using Bayesian analysis. With this approach, prior knowledge about the problem can be incorporated. Here, the parameter should be given a prior density, and two types of priors are used as random priors: the half-Cauchy (HC) and the normal (N) priors. In the numerical integration, prior distributions that are not completely flat provide enough information to allow the numerical approximation algorithm to continue exploring the target posterior density. The HC distribution features such shapes, and its mean and variance do not exist, but its mode is equal to zero. Regarding the HC distribution, when the parameter becomes 25, the pdf is flat, but not entirely. According to [33], depending upon the information necessary, uniform or HC is a better choice of prior. Thus, for the parameter *θ*, *HC*(25) and *N*(0, 1000) are used as prior distributions.

• **Case 1:** *θ* ∼ *HC*(25), then the posterior pdf is given by

$$
\pi(\theta/\mathbf{x}) \propto l(\theta) \times \frac{2 \times 25}{\pi(\theta^2 + 25^2)}. \tag{13}
$$

• **Case 2:** *θ* ∼ *N*(0, 1000), then the posterior pdf is given by

$$
\pi(\theta/\mathbf{x}) \propto l(\theta) \times \frac{1}{\sqrt{2\pi \times 10^3}} e^{-\frac{\theta^2}{2 \times 10^3}}\,\tag{14}
$$

where *l*(*θ*) is derived from (12).

Since (13) and (14) are not in closed form, one may use numerical integration or MCMC methods.

## **5. Simulation Study**

*5.1. Simulation for ML, LS, and WLS Estimates*

In the simulation study, the Monte Carlo simulation was done in order to prove the efficiency of the model by using different estimation methods such as ML, LS, and WLS. The estimates were calculated for true values of parameters for *N* = 1000 samples of sizes 25, 50, 75, 100, 150, 200, and 500 and the following quantities were computed:


where *<sup>υ</sup>* ∈ { <sup>ˆ</sup> *θ*, ˜ *θ*, ˘ *θ*}, and the index *i* indicates the considered number of the sample. Tables 1–3 show the simulation results corresponding to the ML, LS, and WLS estimation methods. A graphical comparison of the MSEs obtained from the three methods is presented in Figure 4. Among these methods, the ML estimation provides accurate estimates of the parameter with the least MSE. As in all the estimation methods, the MSE decreases as sample size *n* increases, as expected.

**Table 1.** Simulation results: MLEs, bias, and MSE.




**Table 2.** Simulation results: LS estimates, bias, and MSE.


**Table 3.** Simulation results: WLS estimates, bias, and MSE.


**Figure 4.** Graphical comparison of the MSEs obtained from ML, LS, and WLS estimation methods for (**a**) *θ* = 0.26, (**b**) *θ* = 0.45, (**c**) *θ* = 0.60, (**d**) *θ* = 1.0, (**e**) *θ* = 1.5 and (**f**) *θ* = 2.0.

#### *5.2. Simulation for Bayesian Estimates*

For the *UTD* parameter *θ*, in Section 4.3, the HC distribution (Case 1) and N distribution (Case 2) were motivated as prior distributions. As a result, the posterior summary results, such as means, standard deviations (SDs), Monte Carlo errors (MCEs), lower bounds (LBs), and upper bounds (UBs) of the 95% confidence intervals and medians are summarized in Tables 4 and 5 for Case 1 and Case 2, respectively. In both cases, increasing sample size leads to a decrease in SD and MCE, which predicts a consistency in the Bayesian estimates.



**Table 4.** *Cont.*


**Table 5.** Posterior summary results (Case 2: *θ* ∼ *N*(0, 1000)).


**Table 5.** *Cont.*


## **6. Application**

*6.1. Methodology*

With the help of two real-life data sets, the superiority of the *UTD* is illustrated. The first data set is from [34]. The data are the maximum flood level (in millions of cubic feet per second) for the Susquehanna River at Harrisburg, Pennsylvania. The second data set represents times between failures of secondary reactor pumps (see [35]).

An analysis of the total time on test (TTT) plot is used to identify the shape of the underlying hrf of the data. We specify that the hrf is decreasing, increasing, bathtub-shaped, and upside-down bathtub-shaped if the empirical TTT transform is convex, concave, convex then concave, and concave then convex, respectively. Thus, it will be shown that the first data set has an increasing hrf, while the second data set has a bathtub-shaped hrf. For the sake of comparison, the following lifetime distributions were considered: Log-Lindley distribution (*LLD*) (see [6]), unit Lindley distribution (*ULD*) (see [9]), Topp– Leone distribution (*TLD*) (see [3]), one-parameter Kumaraswamy distribution (*OKD*), power distribution (*PD*), transmuted distribution (*TD*), two-parameter Kumaraswamy distribution (*KD*), beta distribution (*BD*), exponentiated Topp–Leone distribution (*ETLD*), and unit Burr-III distribution (*UBD*). Using R software, the MLEs of all these distributions' parameters are computed along with information criteria, which are listed below.


Here, log *l* denotes the estimated value of the maximum log-likelihood, *p* denotes the number of parameters and *n* denotes the number of observations. The Kolmogorov– Smirnov (KS) test is also used to test the goodness of fit for all the data sets of the *UTD* and other distributions. Nonparametrically, this test measures how close the empirical distribution and the fitted distribution are. The AIC, AICc, CAIC, and BIC measure the adequacy, while the KS test measures the fit of each distribution. All the computations were performed using R software.

#### *6.2. Flood Level Data*

The data set is from [34]. The data represent the maximum flood level (in millions of cubic feet per second) for the Susquehanna River at Harrisburg, Pennsylvania, over 20 four-year periods from 1890 to 1969. Table 6 provides the measurements of the data set.

**Table 6.** Flood level data.


In Table 7, the MLEs of the parameters of every distribution used and the observed KS test statistics of each distribution are given. We conclude that the *UTD* gives the best fit because it has the smallest KS value and the largest *p*-value. Table 8 shows that *UTD* has the largest maximum log-likelihood value and the smallest AIC, AICc, CAIC, and BIC values among other models. As a result, it can be concluded that the *UTD* model is more successful than the other models for the flood level data.

**Table 7.** The MLEs of the parameters with their standard errors, KS values, and the associated *p*-values.



**Table 8.** Log-likelihood, AIC, AICc, CAIC, and BIC values for flood level data.

Figure 5a,b represent the TTT plot and histogram of flood level data, respectively. The TTT plot indicates an increasing hrf, case covered by the *UTD*, and histogram depicts how well the proposed model fits the data, compared to other models.

**Figure 5.** (**a**) TTT plot and (**b**) histogram for the flood level data.

The estimated variance of the MLE ˆ *θ* of the *UTD* parameter *θ* for the flood level data is given by

$$I^{-1}(\delta) = 0.0112\dots$$

Therefore, an approximate 95% confidence interval for *θ* is [1.0114, 1.4254]. Table 9 gives the median of bootstrap estimates and bootstrap confidence intervals.

**Table 9.** The median of bootstrap estimators and bootstrap confidence intervals.


*6.3. Times between Failures of Secondary Reactor Pumps Data*

The data represent times between failures of secondary reactor pumps (see [35]). Here, a normalization operation is carried out by dividing the original data by 10, in order to obtain data between 0 and 1. Table 10 provides the measurements of the transformed data.

**Table 10.** Secondary reactor pumps data.


Listed in Table 11 are the MLEs of the parameters and the observed KS test statistic for each distribution. In terms of KS value and *p*-value, *UTD* is the best model. On the basis of Table 12, it can be seen that the *UTD* has the largest log-likelihood value, while the AIC, AICc, CAIC, and BIC values are the smallest. Therefore, the *UTD* model fits the secondary reactor pump failure data better than the other models.

**Table 11.** The MLEs of the parameters and their standard errors, KS values and the associated *p*-values.




Figure 6a,b represent the empirical TTT plot and histogram of the secondary reactor pumps data. A bathtub hrf is indicated by the TTT plot for the secondary reactor pumps data. From the histogram, it can be seen that the empirical line is closer to the fitted line of the *UTD* model than other models. The estimated variance of the MLE ˆ *θ* of the *UTD* parameter *θ* for the times between failures of secondary reactor pumps data data is calculated by

$$\int^{-1}(\hat{\theta}) = 0.0008$$

Therefore, an approximate 95% confidence interval for *θ* is given by [0.3057, 0.4193]. Table 13 provides the median of bootstrap estimates and bootstrap confidence intervals.

**Table 13.** The median of bootstrap estimators and bootstrap confidence intervals.


**Figure 6.** (**a**) TTT plot and (**b**) histogram for secondary reactor pumps data.
