A New Generalization of the Generalized Inverse Rayleigh Distribution with Applications

Abstract

In this article, a new four-parameter lifetime model called the beta generalized inverse Rayleigh distribution (BGIRD) is defined and studied. Mixture representation of this model is derived. Curve’s behavior of probability density function, reliability function, and hazard function are studied. Next, we derived the quantile function, median, mode, moments, harmonic mean, skewness, and kurtosis. In addition, the order statistics and the mean deviations about the mean and median are found. Other important properties including entropy (Rényi and Shannon), which is a measure of the uncertainty for this distribution, are also investigated. Maximum likelihood estimation is adopted to the model. A simulation study is conducted to estimate the parameters. Four real-life data sets from difference fields were applied on this model. In addition, a comparison between the new model and some competitive models is done via information criteria. Our model shows the best fitting for the real data.

1. Introduction

Modeling and analysis of lifetime phenomena are important aspects of statistical work in a wide variety of scientific and technological fields. The field of lifetime data analysis has grown and expanded rapidly with respect to methodology, theory, and fields of applications. In the context of modeling the real-life phenomena, continuous probability distributions and many generalizations or transformation methods have been proposed. These generalizations, obtained either by adding one or more shape parameters or by changing the functional form of the distribution, make the models more sufficient for many applications.

The inverse Rayleigh distribution has many applications in lifetime studies. Sometimes it could not fit skewed data or near symmetrical data [1]. This paper proposes a solution by introducing a new two-parameter extension of this distribution through the use of beta generator distribution. Furthermore, the extended parameters control the skewness of the proposed model.

The generalized inverse Rayleigh distribution (GIRD) is a very useful life distribution. The GIRD is an important distribution in statistics and operations research. Ref. [2] estimated parameters of GIRD based on randomly censored samples. Ref. [3] discussed the inference of GIRD with real data applications. Moreover, ref. [4] estimated parameters of GIRD based on progressive type II censoring. In addition, ref. [5] studied the characterizations through the generalized inverted Rayleigh distribution.

The probability density function (pdf) of the generalized inverse Rayleigh distribution is given by

$$g\left(x\right)=\frac{2\eta }{{\rho }^2{x}^3}{e}^{-{\left(\rho x\right)}^{-2}}{[1-e^{-{\left(\rho x\right)}^{-2}}]}^{\eta -1},x>0,\eta ,\rho >0.$$

(1)

and the corresponding cumulative distribution function (cdf) is provided by

$$G\left(x\right)=1-{[1-e^{-{\left(\rho x\right)}^{-2}}]}^{\eta },x>0,\eta ,\rho >0.$$

(2)

where $\eta$ and $\rho$ are the shape and the scale parameters, respectively.

Recently, attempts have been made to define new families of probability distributions that extend well-known families of distributions and at the same time provide great flexibility in modeling data in practice. One such class of distributions is generated from the beta random variable. This class of generalized distributions has been receiving considerable attention over the last years. A general class of distribution was introduced by [6] as:

$$F\left(x\right)=I_{G\left(x\right)}(a,b),a>0,b>0,$$

(3)

where $G\left(x\right)$ is the cdf of a baseline distribution and $I_u(a,b)$ is the incomplete beta function ratio such that

$$I_u(a,b)=\frac{\beta (u;a,b)}{\beta (a,b)}=\frac{1}{\beta (a,b)}{\int }_0^u{\varpi }^{a-1}{(1-\varpi )}^{b-1}d\varpi ,a,b>0,$$

(4)

where $\beta (a,b)={\int }_0^1{\varpi }^{a-1}{(1-\varpi )}^{b-1}d\varpi$ is the beta function. The skewness of the distribution is controlled by the two shape parameters a and b.

The pdf of the beta-G distribution has the following form

$$f\left(x\right)=\frac{g\left(x\right)}{\beta (a,b)}G{\left(x\right)}^{a-1}{[1-G\left(x\right)]}^{b-1},$$

(5)

where $g\left(x\right)$ is the pdf of a baseline distribution.

We can obtain the beta generalized inverse Rayleigh distribution (BGIRD) by using the pdf and cdf of GIRD in (1) and (2) as a baseline distribution for (5) and (3), which is a generalization of GIRD. The BGIRD is expected to be more flexible in real applications.

Some classes of beta-generated distributions have received considerable attention in recent years. Ref. [7] studied the beta inverse Rayleigh distribution. They provided various properties, including the quantile function, moments, mean deviations, Bonferroni and Lorenz curves, Rényi and Shannon entropies, and order statistics, as well as the maximum likelihood estimates. In addition, Ref. [8] studied the beta Rayleigh distribution and discussed some properties of the distribution. Additionally, maximum likelihood estimation and the information matrix were obtained. Ref. [9] proposed the beta generalized inverted exponential distribution and derived various statistical properties. They obtained the maximum likelihood estimators, asymptotic Fisher information matrix, and confidence interval estimates of the parameters. Moreover, applications on real data sets were provided.

The aim of this study is to introduce a generalization of the generalized inverted Rayleigh distribution termed the beta generalized inverted Rayleigh distribution. We hope that this generalization shall attract wide applications. In this study, many important characteristics of the distribution are studied. As well as estimating the parameters of the distribution using the maximum likelihood estimation method and Bayes estimation method under complete samples. A simulation study is adopted to find the estimates of this distribution according to the methods mentioned above for different sample sizes. Finally, the distribution was applied to five real data in different fields. As well, the distribution is compared with other models using information criteria.

2. The Beta Generalized Inverse Rayleigh Distribution

In this section, we introduce the four parameter beta generalized inverse Rayleigh distribution. The cdf of the BGIRD could be written using the incomplete beta function defined in (3) as follows

$$\begin{array}{ccc}\hfill F\left(x\right)& =& I_{1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta }}(\nu ,\tau )=\frac{\beta (1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta };\nu ,\tau )}{\beta (\nu ,\tau )}\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& =& \frac{1}{\beta (\nu ,\tau )}{\int }_0^{1-{[1-e^{-{\left(\rho x\right)}^{-2}}]}^{\eta }}{\varpi }^{\nu -1}{(1-\varpi )}^{\tau -1}d\varpi ,x>0,\nu ,\tau ,\eta and\rho >0,\hfill \end{array}$$

(7)

where $\eta$, $\nu$, and $\tau$ are shape parameters and $\rho$ is a scale parameter.

2.1. Probability Density Function of BGIRD

The pdf of the BGIRD takes the form

$$\begin{array}{c}f\left(x\right)=\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2{x}^3}{e}^{-{\left(\rho x\right)}^{-2}}{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta \tau -1}{[1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta }]}^{\nu -1},\\ x>0,\nu ,\tau ,\eta and\rho >0.\end{array}$$

(8)

By using the power series expansion, for $b>0$ a real non-integer number and $\left|Z\right|<1$,

$${(1-Z)}^{b-1}={\displaystyle \sum _{i=0}^{\infty }}{(-1)}^i\left(\genfrac{}{}{0pt}{}{b-1}{i}\right)Z^i,-1<Z<1,$$

(9)

we can rewrite the pdf of the BGIRD in (8) as an infinite power series in the following forms

$$\begin{array}{c}f\left(x\right)=\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2{x}^3}{e}^{-{\left(\rho x\right)}^{-2}}{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k}{\nu \beta (\nu -k,k+1)}{[1-e^{-{\left(\rho x\right)}^{-2}}]}^{\eta (\tau +k)-1},\\ x>0,\nu ,\tau ,\eta and\rho >0,\end{array}$$

(10)

and

$$\begin{array}{c}f\left(x\right)=\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2{x}^3}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}{e}^{-(i+1){\left(\rho x\right)}^{-2}}}{\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)},\\ x>0,\nu ,\tau ,\eta and\rho >0,\end{array}$$

(11)

where $\nu$ is a real number.

2.2. Cumulative Distribution Function of BGIRD

From (10) and (11), the corresponding cdf, respectively, can be written as follows

$$F\left(x\right)=\frac{1}{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k[1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta (\tau +k)}]}{\nu (\tau +k)\beta (\nu -k,k+1)},x>0,\nu ,\tau ,\eta and\rho >0,$$

(12)

and

$$\begin{array}{c}F\left(x\right)=\frac{\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}{e}^{-(i+1){\left(\rho x\right)}^{-2}}}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)},\\ x>0,\nu ,\tau ,\eta and\rho >0.\end{array}$$

(13)

Figure 1 shows the curves of the BGIR probability density function.

From Figure 1, we note that the distribution is unimodal and positively skewed.

2.3. Mixture Representation

In this section, we derive mixture representations for the pdf and cdf of X in order to obtain a simple form for the BGIRD pdf. By using Equation (10), the pdf of X can be written as

$$f\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}2(i+1){\upsilon }_i{\rho }^{-2}{x}^{-3}{e}^{-(i+1){\left(\rho x\right)}^{-2}},$$

(14)

$${\upsilon }_i=\frac{{(-1)}^{k+i}\eta \Gamma \left(\nu \right)\Gamma \left(\eta (\tau +k)\right)}{\beta (\nu ,\tau )k!\Gamma (\nu -k)(i+1)!\Gamma \left(\eta \right(\tau +k)-i)}.$$

Equation (14) can be rewritten as

$$f\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{\upsilon }_i{h}_{i+1}\left(x\right),$$

(15)

where $h_{i+1}\left(x\right)$ is the inverse Rayleigh (IR) pdf with scale parameter (${\rho }^{-2}(i+1)$).

Equation (15) reveals that the BGIRD density function can be expressed as a mixture of IR densities. Therefore, several of its structural properties can be derived from those of the inverse Rayleigh distribution (IRD).

By integrating (15), we obtain

$$F\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{\upsilon }_i{H}_{i+1}\left(x\right),$$

where $H_{i+1}\left(x\right)$ is the cdf of the IR model with scale parameter (${\rho }^{-2}(i+1)$).

2.4. The Reliability Function

Suppose X is a BGIR random variable which represents the lifetime of a unit, and t represents time, then the probability that a unit X survives beyond time t is called the reliability at time t. The reliability, $R\left(t\right)$, at time t is given as follows:

$$\begin{array}{cc}\hfill R\left(t\right)& =1-F\left(t\right)\hfill \\ \hfill & =1-\frac{1}{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k[1-{(1-e^{-{\left(\rho t\right)}^{-2}})}^{\eta (\tau +k)}]}{\nu (\tau +k)\beta (\nu -k,k+1)},t>0,\nu ,\tau ,\eta and\rho >0.\hfill \end{array}$$

Figure 2 shows the reliability function of the BGIRD for various parameter choices.

2.5. The Hazard Function

The rate at which failure of a unit occurs per unit of time relative to the proportion of the population which has not yet failed, is the hazard function, $h\left(t\right)$. The hazard function of the BGIRD is given as

$$\begin{array}{cccc}\hfill h\left(t\right)& =\frac{f\left(t\right)}{1-F\left(t\right)}\hfill & & \\ \hfill & =\frac{\frac{2\eta e^{-{\left(\rho t\right)}^{-2}}}{\beta \left(\nu ,\tau \right){\rho }^2{x}^3}{\sum }_{k=0}^{\infty }\frac{{\left(-1\right)}^k}{\nu \beta \left(\nu -k,k+1\right)}{\left[1-e^{-{\left(\rho t\right)}^{-2}}\right]}^{\eta \left(\tau +k\right)-1}}{1-\frac{1}{\beta \left(\nu ,\tau \right)}{\sum }_{k=0}^{\infty }\frac{{\left(-1\right)}^k\left[1-{\left(1-e^{-{\left(\rho t\right)}^{-2}}\right)}^{\eta \left(\tau +k\right)}\right]}{\nu \left(\tau +k\right)\beta \left(\nu -k,k+1\right)}},t>0,\nu ,\tau ,\eta and\rho >0.\hfill & \end{array}$$

Figure 3 shows the hazard function of the BGIRD for various parameter choices.

It is shown from Figure 3 that the hazard function increases and then decreases. That is, the upside down bathtub hazard function which indicates that the risk of failing decreases as soon as the item has passed a specific time, during which it may have experienced some type of stress. Thus, the BGIRD shows good statistical behavior based on this function and could be a flexible model for fitting data in many fields.

2.6. Special Sup-Models

The BGIRD in Equation (8) represents a generalization of several distributions that have been considered in the literature.

• In particular, BGIRD becomes GIRD ($\eta$, $\rho$) when $\nu$ and $\tau =1$.
• The beta inverse Rayleigh distribution BIRD ($\nu ,\tau ,\rho$) is clearly a special case of BGIRD when $\eta$ = 1.
• IRD ($\rho$) can be obtained from (8) by making $\nu ,\tau ,\eta$ = 1.
• In addition, the exponentiated Rayleigh distribution ERD ($\eta$, $\rho$) is a special case of BGIRD when $\nu$ = $\tau$ = 1 and the random variable Y = 1/X.
• If $\nu$ = $\tau$ = 1 and $\eta$ = 1 in Equation (8), the random variable Y = 1/X has the Rayleigh distribution ($\rho$).

3. Statistical Properties

3.1. Quantile Function

The quantile function of the BGIRD corresponding to (6) is

$$q\left(u\right)={\left[\frac{-1}{{\rho }^{2}log\left[1-{[1-I_u^{-1}(\nu ,\tau )]}^{\frac{1}{\eta }}\right]}\right]}^{\frac{1}{2}},0<u<1,$$

(16)

where $I_u^{-1}(\nu ,\tau )$ is the inverse of the incomplete beta function defined in (3)

$$I_u(\nu ,\tau )=\frac{1}{\beta (\nu ,\tau )}{\int }_0^u{\varpi }^{\nu -1}{(1-\varpi )}^{\tau -1}d\varpi .$$

3.2. Median

We can derive the median of the BGIRD as the following:

Set $F\left(m\right)=\frac{1}{2}$, then by using Equation (12), we have:

$$\frac{1}{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k[1-{(1-e^{-{\left(\rho m\right)}^{-2}})}^{\eta (\tau +k)}]}{\nu (\tau +k)\beta (\nu -k,k+1)}=\frac{1}{2}$$

when a is positive and $\left|{(1-e^{-{\left(\rho m\right)}^{-2}})}^{\eta }\right|<1$, we can take the first approximation by putting k = 0 and $\nu$ = 1, then we get:

$${(1-e^{-{\left(\rho m\right)}^{-2}})}^{\eta \tau }=\frac{1}{2}$$

After some simplifications, we can write the median of the BGIRD in the form:

$$m={\left[\frac{1}{-{\rho }^{2}log\left[1-{\left(\frac{1}{2}\right)}^{\frac{1}{\eta \tau }}\right]}\right]}^{\frac{1}{2}}.$$

(17)

3.3. Mode

The mode of the BGIRD can be found by solving the following equation with respect to x:

$$\frac{df\left(x\right)}{dx}=0.$$

Then, by using (8), we get:

$$f\left(x\right)\left[\frac{2}{{\rho }^2{x}^3}-\frac{3}{x}-\frac{2(\eta \tau -1)}{{\rho }^2{x}^3}{\left[e^{{\left(\rho x\right)}^{-2}}-1\right]}^{-1}+\frac{2\eta (\nu -1)}{{\rho }^2{x}^3}{\left[e^{{\left(\rho x\right)}^{-2}}-1\right]}^{-1}{\left[(1-e^{-{\left(\rho x\right)}^{-2}})-1\right]}^{-1}\right]=0,$$

since $f\left(x\right)>0$, the mode is the solution of the following equation:

$$\frac{2}{{\rho }^2{x}^3}-\frac{3}{x}-\frac{2(\eta \tau -1)}{{\rho }^2{x}^3}{\left[e^{{\left(\rho x\right)}^{-2}}-1\right]}^{-1}+\frac{2\eta (\nu -1)}{{\rho }^2{x}^3}{\left[e^{{\left(\rho x\right)}^{-2}}-1\right]}^{-1}{\left[(1-e^{-{\left(\rho x\right)}^{-2}})-1\right]}^{-1}=0.$$

(18)

Equation (18) is a non-linear equation and cannot have an analytic solution in x. Therefore, it has to be solved numerically. The mode of the BGIRD can be obtained by solving (18) numerically using the Newton–Raphson method.

3.4. Moments

The rth moment of the BGIRD random variable X is given by:

$${\mu }_r^{\prime }={\int }_0^{\infty }{x}^{r}f\left(x\right)dx$$

using the form of the pdf in (14), we have:

$${\mu }_r^{\prime }={\int }_0^{\infty }{\displaystyle \sum _{i=0}^{\infty }}2(i+1){\upsilon }_i{\rho }^{-2}{x}^{-3+r}{e}^{-(i+1){\left(\rho x\right)}^{-2}}dx$$

By setting $w=(i+1){\left(\rho x\right)}^{-2}$, we get:

$$\begin{array}{cc}\hfill {\mu }_r^{\prime }& ={\displaystyle \sum _{i=0}^{\infty }}{\upsilon }_i{(i+1)}^{\frac{r}{2}}{\rho }^{-r}{\int }_0^{\infty }{e}^{-w}{w}^{-\frac{r}{2}}dw\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^{\infty }}{\upsilon }_i{(i+1)}^{\frac{r}{2}}{\rho }^{-r}\Gamma (1-\frac{r}{2}),\hfill \end{array}$$

(19)

where $\Gamma (.)$ is the gamma function, which is defined as $\Gamma \left(z\right)={\int }_0^{\infty }{e}^{-t}{t}^{z-1}dt.$

Setting $r=1$ in Equation (19), we obtain the mean of X.

$$\begin{array}{cc}\hfill {\mu }_1^{\prime }& ={\displaystyle \sum _{i=0}^{\infty }}{\upsilon }_i{(i+1)}^{\frac{1}{2}}{\rho }^{-1}\Gamma \left(\frac{1}{2}\right),\Gamma \left(\frac{1}{2}\right)=\sqrt{\pi }.\hfill \end{array}$$

(20)

We noticed that the even moments of (X) do not exist, as the second and fourth moments.

Using the relation between the central and non-central moments, we obtain the nth central moment of X, say ${\mu }_n$, as follows

$$\begin{array}{cc}\hfill {\mu }_n& ={\displaystyle \sum _{r=0}^n}{\displaystyle \sum _{i=0}^{\infty }}\left(\genfrac{}{}{0pt}{}{n}{r}\right){(-{{\mu }_1}^{\prime })}^{n-r}{\upsilon }_i{(i+1)}^{\frac{r}{2}}{\rho }^{-r}\Gamma (1-\frac{r}{2}).\hfill \end{array}$$

3.5. Harmonic Mean

The harmonic mean is given by (see [10]):

$$H_M\left(X\right)=\frac{1}{E\left(X^{-1}\right)}={\left[{\int }_0^{\infty }{x}^{-1}f\left(x\right)dx\right]}^{-1}.$$

Using Equation (11), the harmonic mean of the BGIRD can be derived as follows:

Let

$$\begin{array}{cc}\hfill I& ={\int }_0^{\infty }{x}^{-1}f\left(x\right)dx\hfill \\ & ={\int }_0^{\infty }\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2{x}^4}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}{e}^{-(i+1){\left(\rho x\right)}^{-2}}}{\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}dx\hfill \end{array}$$

By setting $u=(i+1){\left(\rho x\right)}^{-2}$, we get:

$$\begin{array}{cc}\hfill I& =\frac{\rho \eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{{(i+1)}^{\frac{3}{2}}\nu \left(\eta (\tau +k)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}{\int }_0^{\infty }{u}^{\frac{1}{2}}{e}^{-u}du\hfill \\ & =\frac{\rho \eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{{(i+1)}^{\frac{3}{2}}\nu \left(\eta (\tau +k)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}A\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill A& ={\int }_0^{\infty }{u}^{\frac{1}{2}}{e}^{-u}du=\Gamma \left(\frac{3}{2}\right)=\frac{1}{2}\Gamma \left(\frac{1}{2}\right)=\frac{\sqrt{\pi }}{2}\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill I& =\frac{\rho \eta \Gamma \left(\frac{3}{2}\right)}{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{{(i+1)}^{\frac{3}{2}}\nu \left(\eta (\tau +k)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}\hfill \\ & =\frac{\rho \eta \sqrt{\pi }}{2\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{{(i+1)}^{\frac{3}{2}}\nu \left(\eta (\tau +k)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}\hfill \end{array}$$

Therefore, the harmonic mean of BGIRD is given by:

$$\begin{array}{cc}\hfill H_M\left(X\right)& ={\left[I\right]}^{-1}\hfill \\ & ={\left[\frac{\rho \eta \sqrt{\pi }}{2\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{{(i+1)}^{\frac{3}{2}}\nu \left(\eta (\tau +k)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}\right]}^{-1}.\hfill \end{array}$$

(21)

3.6. Skewness and Kurtosis

Based on the fact that the second and fourth moments of the BGIRD are non-existent, usual skewness and kurtosis based on moments could not be found.

However, measures based on quantile, such as Bowley skewness (see [11]) and Moors kurtosis (see [12]), can quantify asymmetry and the peakedness of a given distribution. These measures exist even when moments are not available.

Bowley skewness and Moors kurtosis are defined, respectively, by

$$\begin{array}{cc}\hfill B& =\frac{q(3/4)-2q(1/2)+q(1/4)}{q(3/4)-q(1/4)}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill M& =\frac{q(7/8)-q(5/8)+q(3/8)-q(1/8)}{q(6/8)-q(2/8)},\hfill \end{array}$$

where $q(.)$ is defined in (16).

The Bowley skewness, Moors kurtosis, mean, median, mode, and harmonic mean of the BGIRD for various values of $\nu$, $\tau$, $\eta$, and $\rho$ are shown in

Table 1. The skewness, kurtosis, mean, median, mode, and harmonic mean of the BGIRD.

$\mathit{\nu }$	$\mathit{\tau }$	$\mathit{\eta }$	$\mathit{\rho }$	Skewness	Kurtosis	Mean	Median	Mode	Harmonic Mean
2	2	2	1.5	0.13731	1.29818	0.637261	0.601615	0.545538	0.592602
2	2	2	2	0.13731	1.29818	0.477946	0.451212	0.409154	0.444451
2	2	0.5	2	0.31697	1.58531	1.42599	0.93221	0.60726	0.862293
2	2	1.5	2	0.162677	1.32328	0.543969	0.50147	0.439831	0.492135
2	2	5	2	0.0792963	1.25709	0.357224	0.349688	0.336093	0.346181
1	2	2	2	0.136143	1.29595	0.392529	0.368785	0.331165	0.360099
3	2	2	2	0.137703	1.29892	0.536444	0.507568	0.462275	0.501744
2	0.75	2	2	0.240794	1.4393	0.770942	0.63456	0.498165	0.616217
2	1	2	2	0.203663	1.37855	0.64575	0.566454	0.470633	0.553656
2	4	2	2	0.0921975	1.26287	0.387276	0.376634	0.35771	0.372261
1	1	1.5	0.5	0.24225	1.43601	2.46538	2.00588	1.53797	1.92415
1	1	1	1.5	0.30686	1.57048	3.54491	0.800748	0.54433	0.75225

.

Table 1 shows the results of studying BGIRD’s behavior. Kurtosis and skewness values remain constant for fixed values of $\nu$, $\tau$, and $\eta$. We also note that, with increasing $\rho$, the mean, median, mode, and harmonic mean are decreasing. As the value of $\eta$ increases and the other parameters are fixed, the skewness, kurtosis, mean, median, mode and harmonic mean decrease. Furthermore, the skewness, kurtosis, mean, mode, and harmonic mean increase as we increase the value of $\nu$ and with the stability of the values of $\tau$, $\eta$, and $\rho$. Additionally, for different values of $\tau$ and fixed values of $\nu$, $\eta$, and $\rho$, the skewness, kurtosis, mean, median, mode, and harmonic mean decrease as $\tau$ increases. Moreover, we found that our results for $\nu$ = $\tau$ = 1 are exactly the same as the results in [5]. All graphs and computations presented were carried out by Mathematica 12.0.

3.7. Mean Deviations

The mean deviation is a measure of dispersion derived by computing the mean of the absolute values of the differences between the observed values of a variable and the mean or median of the variable. It is called average deviation (see [13]).

The mean deviation about the mean and the mean deviation about the median are, respectively, defined by:

$$D\left(\mu \right)=E\left(\left|X-\mu \right|\right)$$

and

$$D\left(m\right)=E\left(\left|X-m\right|\right)$$

where $\mu =E\left(X\right)$ and $m=q(1/2)$.

The mean deviation about the mean and the median is given by the following theorems.

3.7.1. The Mean Deviation about the Mean

Theorem 1.

If X follows BGIRD, then the mean deviation about the mean is in the form

$$D\left(\mu \right)=\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}\left[\mu e^{-(i+1){\left(\rho \mu \right)}^{-2}}-{\left(\frac{i+1}{{\rho }^2}\right)}^{\frac{1}{2}}\Gamma (\frac{1}{2},\frac{i+1}{{\rho }^2{\mu }^2})\right]}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}.$$

Proof. 

The mean deviation about the mean can be written as

$$\begin{array}{cccc}\hfill D\left(\mu \right)& =E\left(\left|X-\mu \right|\right)={\int }_0^{\infty }\left|x-\mu \right|f\left(x\right)dx\hfill & \hfill & \\ \hfill & ={\int }_0^{\mu }(\mu -x)f\left(x\right)dx+{\int }_{\mu }^{\infty }(x-\mu )f\left(x\right)dx\hfill & \hfill & \\ \hfill & =2{\int }_0^{\mu }(\mu -x)f\left(x\right)dx\hfill & \hfill & \\ \hfill & =2\mu F\left(\mu \right)-2{\int }_0^{\mu }xf\left(x\right)dx\hfill & \hfill & \\ \hfill & =2\mu F\left(\mu \right)-2{\int }_0^{\mu }xdF\left(x\right)\hfill & \hfill & \\ \hfill & =2{\int }_0^{\mu }F\left(x\right)dx\hfill & \hfill & \hfill \end{array}$$

(22)

by using Equation (13), the mean deviation of the BGIRD can be derived as:

$$\begin{array}{cccc}\hfill D\left(\mu \right)& =2{\int }_0^{\mu }\frac{\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}{e}^{-(i+1){\left(\rho x\right)}^{-2}}}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}dx\hfill & \hfill & \\ \hfill & =\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}\hfill & \hfill & \\ \hfill & \times {\int }_0^{\mu }{e}^{-(i+1){\left(\rho x\right)}^{-2}}dx\hfill & \hfill & \\ \hfill & =\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}B,\hfill & \hfill \end{array}$$

where

$$\begin{array}{cccc}\hfill B& ={\int }_0^{\mu }{e}^{-(i+1){\left(\rho x\right)}^{-2}}dx\hfill & \hfill & \\ \hfill & =\mu e^{-(i+1){\left(\rho \mu \right)}^{-2}}-{\left(\frac{i+1}{{\rho }^2}\right)}^{\frac{1}{2}}\Gamma (\frac{1}{2},\frac{i+1}{{\rho }^2{\mu }^2})\hfill & \hfill \end{array}$$

then, we get the mean deviation about the mean of the BGIRD as:

$$\begin{array}{cc}\hfill D\left(\mu \right)& =\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}\left[\mu e^{-(i+1){\left(\rho \mu \right)}^{-2}}-{\left(\frac{i+1}{{\rho }^2}\right)}^{\frac{1}{2}}\Gamma (\frac{1}{2},\frac{i+1}{{\rho }^2{\mu }^2})\right]}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}.\hfill \end{array}$$

(23)

□

3.7.2. The Mean Deviation about the Median

Theorem 2.

If X follows BGIRD, then the mean deviation about the median is in the form

$$D\left(m\right)=\mu -m+\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}\left[m{e}^{-(i+1){\left(\rho m\right)}^{-2}}-{\left(\frac{i+1}{{\rho }^2}\right)}^{\frac{1}{2}}\Gamma (\frac{1}{2},\frac{i+1}{{\rho }^2{m}^2})\right]}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}.$$

Proof. 

The mean deviation about the median can be written as

$$\begin{array}{cc}\hfill D\left(m\right)& =E\left(\left|X-m\right|\right)={\int }_0^{\infty }\left|x-m\right|f\left(x\right)dx\hfill \\ \hfill & =2{\int }_0^m(m-x)f\left(x\right)dx+{\int }_0^{\infty }(x-m)f\left(x\right)dx\hfill \\ \hfill & =2mF\left(m\right)-2{\int }_0^{m}xf\left(x\right)dx+E(X-m)\hfill \\ \hfill & =2mF\left(m\right)-2{\int }_0^{m}xf\left(x\right)dx+\mu -m\hfill \\ \hfill & =\mu -m+2{\int }_0^{m}F\left(x\right)dx\hfill \end{array}$$

(24)

The third term in Equation (24) is the same as Equation (22) where the upper limit of the integration is m instead of $\mu$. Hence, by substituting the result of (23) in (24), we get:

$$\begin{array}{cc}\hfill D\left(m\right)& =\mu -m+\frac{2\eta }{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^{k+i}\left[m{e}^{-(i+1){\left(\rho m\right)}^{-2}}-{\left(\frac{i+1}{{\rho }^2}\right)}^{\frac{1}{2}}\Gamma (\frac{1}{2},\frac{i+1}{{\rho }^2{m}^2})\right]}{(i+1)\nu \left(\eta \right(\tau +k\left)\right)\beta (i+1,\eta (\tau +k)-i)\beta (\nu -k,k+1)}.\hfill \end{array}$$

(25)

□

3.8. Rényi and Shannon Entropies

There are two important entropy measures: the Shannon entropy and its generalization, which is known as the Rényi entropy. The entropy of a random variable quantifies its associated uncertainty (see [14]).

The Rényi and Shannon entropies are defined, respectively, by:

$$\begin{array}{cc}\hfill R_{\theta }\left(X\right)& =\frac{1}{1-\theta }log\left[J\left(\theta \right)\right],\hfill \end{array}$$

where $J\left(\theta \right)={\int }_0^{\infty }{f}^{\theta }\left(x\right)dx,\theta >0and\theta \ne 1.$

$$\begin{array}{cc}\hfill W\left(X\right)& =-E\left[log\right(f\left(x\right)\left)\right].\hfill \end{array}$$

Now, the Shannon entropy can be obtained by taking the limit of the Rényi entropy when $\theta \to 1$, as follows:

$$W\left(X\right)=\underset{\theta \to 1}{lim}\left(R_{\theta }\left(X\right)\right)=\underset{\theta \to 1}{lim}\frac{1}{1-\theta }log\left({\int }_0^{\infty }{f}^{\theta }\left(x\right)dx\right).$$

The following theorems presented the forms of Rényi and Shannon entropies for the BGIRD.

3.8.1. The Rényi Entropy for the BGIRD

Theorem 3.

If X follows BGIRD, then the Rényi entropy is in the form

$$\begin{array}{cc}\hfill R_{\theta }\left(X\right)& =\frac{\theta log\eta }{1-\theta }-log\left(2\rho \right)-\frac{\theta log\beta (\nu ,\tau )}{1-\theta }+\frac{1}{1-\theta }log\hfill \\ & \times \left[\Gamma \left(\frac{3\theta -1}{2}\right){\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right){(\theta +i)}^{-\left(\frac{3\theta -1}{2}\right)}\right].\hfill \end{array}$$

Proof. 

By using Equation (8), the Rényi entropy can be written as

$$\begin{array}{cc}\hfill J\left(\theta \right)& ={\int }_0^{\infty }{\left(\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2{x}^3}\right)}^{\theta }{e}^{-\theta {\left(\rho x\right)}^{-2}}{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\theta (\eta \tau -1)}{[1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta }]}^{\theta (\nu -1)}dx\hfill \\ & ={\left(\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2}\right)}^{\theta }{\int }_0^{\infty }{x}^{-3\theta }{e}^{-\theta {\left(\rho x\right)}^{-2}}{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\theta (\eta \tau -1)}{[1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta }]}^{\theta (\nu -1)}dx\hfill \end{array}$$

By setting $u={\left(\rho x\right)}^{-2}$, we get:

$$\begin{array}{cc}\hfill J\left(\theta \right)& ={\left(\frac{2\eta \rho }{\beta (\nu ,\tau )}\right)}^{\theta }{\left(2\rho \right)}^{-1}{\int }_0^{\infty }{u}^{\frac{3}{2}(\theta -1)}{e}^{-\theta u}{(1-e^{-u})}^{\theta (\eta \tau -1)}{[1-{(1-e^{-u})}^{\eta }]}^{\theta (\nu -1)}du,\hfill \end{array}$$

by applying the series defined in Equation (9), we get:

$$\begin{array}{cc}\hfill J\left(\theta \right)& ={\left(\frac{2\eta \rho }{\beta (\nu ,\tau )}\right)}^{\theta }{\left(2\rho \right)}^{-1}{\int }_0^{\infty }{u}^{\frac{3}{2}(\theta -1)}{e}^{-\theta u}{\displaystyle \sum _{k=0}^{\infty }}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right){(1-e^{-u})}^{\eta k+\theta (\eta \tau -1)}du\hfill \\ & ={\left(\frac{2\eta \rho }{\beta (\nu ,\tau )}\right)}^{\theta }{\left(2\rho \right)}^{-1}{\int }_0^{\infty }{u}^{\frac{3}{2}(\theta -1)}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right)\hfill \\ & \times e^{-u(\theta +i)}du\hfill \end{array}$$

$$\begin{array}{cc}\hfill J\left(\theta \right)& ={\left(\frac{2\eta \rho }{\beta (\nu ,\tau )}\right)}^{\theta }{\left(2\rho \right)}^{-1}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right)\hfill \\ \hfill & \times {\int }_0^{\infty }{u}^{\frac{3}{2}(\theta -1)}{e}^{-u(\theta +i)}du\hfill \\ \hfill & ={\left(\frac{2\eta \rho }{\beta (\nu ,\tau )}\right)}^{\theta }{\left(2\rho \right)}^{-1}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right)\frac{\Gamma (\frac{3\theta }{2}-\frac{1}{2})}{{(\theta +i)}^{\left(\frac{3\theta -1}{2}\right)}}.\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill log\left[J\right(\theta \left)\right]& =(\theta -1)log\left(2\rho \right)+\theta log\eta -\theta log\beta (\nu ,\tau )\hfill \\ & +\left[\Gamma \left(\frac{3\theta -1}{2}\right){\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right){(\theta +i)}^{-\left(\frac{3\theta -1}{2}\right)}\right].\hfill \end{array}$$

Hence,

$$\begin{array}{cc}\hfill R_{\theta }\left(X\right)& =\frac{\theta log\eta }{1-\theta }-log\left(2\rho \right)-\frac{\theta log\beta (\nu ,\tau )}{1-\theta }+\frac{1}{1-\theta }log\hfill \\ & \times \left[\Gamma \left(\frac{3\theta -1}{2}\right){\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i=0}^{\infty }}{(-1)}^{k+i}\left(\genfrac{}{}{0pt}{}{\theta (\nu -1)}{k}\right)\left(\genfrac{}{}{0pt}{}{\eta k+\theta (\eta \tau -1)}{i}\right){(\theta +i)}^{-\left(\frac{3\theta -1}{2}\right)}\right].\hfill \end{array}$$

(26)

□

Figure 4 shows the curve of the BGIR Rényi entropy.

3.8.2. The Shannon Entropy for the BGIRD

Theorem 4.

If X follows BGIRD, then the Shannon entropy is in the form

$$\begin{array}{cc}\hfill W\left(X\right)& =-log\left(\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2}\right)+\frac{1}{\beta (\nu ,\tau )}{\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k}{\nu \left(\eta \right(\tau +k\left)\right)\beta (\nu -k,k+1)}\hfill \\ \hfill & \times (\frac{3\eta }{2}{\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^i\left[log\left(\frac{i+1}{{\rho }^2}\right)\right]}{(i+1)\beta (i+1,\eta (\tau +k)-i)}\hfill \\ \hfill & +\eta {\displaystyle \sum _{i=0}^{\infty }}\frac{{(-1)}^i}{{(i+1)}^2\beta (i+1,\eta (\tau +k)-i)}\hfill \\ \hfill & +\eta (\eta \tau -1){\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=1}^{\infty }}\frac{{(-1)}^i}{j(i+j+1)\beta (i+1,\eta (\tau +k)-i)}\hfill \\ \hfill & +(\nu -1){\displaystyle \sum _{j=1}^{\infty }}{\displaystyle \sum _{l=0}^{\infty }}\frac{{(-1)}^l\left(\eta (\tau +k)\right)}{(j+\tau +k)j\beta (l+1,\eta (j+\tau +k)-l)(l+1)}).\hfill \end{array}$$

The proof is presented in Appendix A.

4. Order Statistics

Order statistics deal with the properties and applications of ordered random variables and their functions. In the study of many natural problems related to flood, longevity, breaking strength, atmospheric temperature, atmospheric pressure, wind, etc., the future possibilities in the recurrence of extreme situations are of much importance and accordingly the problem of interest in these cases reduces to that of the extreme observations.

Let $X_1,X_2,....,X_n$ be a random sample from the BGIRD with pdf and cdf as in (14) and (6), respectively, and let $X_{1:n}\le X_{2:n}\le ....\le X_{n:n}$ denote the order statistics obtained from this sample. The pdf of the rth order statistic $X_{r:n}$ is given by (see [15]).

$$\begin{array}{cc}\hfill {f_X}_{r:n}\left(x\right)& =\frac{1}{\beta (r,n-r+1)}{\left(F\left(x\right)\right)}^{r-1}{(1-F\left(x\right))}^{n-r}f\left(x\right).\hfill \end{array}$$

Here, we present an explicit expression for the density function ${f_X}_{r:n}\left(x\right)$ of the rth order statistic $X_{r:n}$ in a random sample of size n from the BGIRD.

Using (6), the pdf of $X_{r:n}$ for the BGIRD is given by

$$\begin{array}{cc}\hfill {f_X}_{r:n}\left(x\right)& =\frac{1}{\beta (r,n-r+1)}f\left(x\right){\left(\frac{\beta (1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta };\nu ,\tau )}{\beta (\nu ,\tau )}\right)}^{r-1}\hfill \\ & \times {\left(1-\frac{\beta (1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta };\nu ,\tau )}{\beta (\nu ,\tau )}\right)}^{n-r}\hfill \end{array}$$

by applying the series defined in Equation (9), we get:

$$\begin{array}{c}\hfill {f_X}_{r:n}\left(x\right)=\frac{f\left(x\right)}{\beta (r,n-r+1)}{\displaystyle \sum _{j=0}^{n-r}}{(-1)}^j\left(\genfrac{}{}{0pt}{}{n-r}{j}\right){\left(\frac{\beta (1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta };\nu ,\tau )}{\beta (\nu ,\tau )}\right)}^{r-1+j}\end{array}$$

Using the series expression for the incomplete beta function:

$$\begin{array}{cc}\hfill I_x(a,b)& =\frac{\beta (x;a,b)}{\beta (a,b)}={\displaystyle \sum _{k=a}^{a+b-1}}\left(\genfrac{}{}{0pt}{}{a+b-1}{k}\right)x^k{(1-x)}^{a+b-1-k},\hfill \end{array}$$

the pdf of $X_{r:n}$ can be written as

$$\begin{array}{cc}\hfill {f_X}_{r:n}\left(x\right)& =\frac{f\left(x\right)}{\beta (r,n-r+1)}{\displaystyle \sum _{j=0}^{n-r}}{(-1)}^j\left(\genfrac{}{}{0pt}{}{n-r}{j}\right)\hfill \\ & \times {\left\{{\displaystyle \sum _{k=\nu }^{\nu +\tau -1}}\left(\genfrac{}{}{0pt}{}{\nu +\tau -1}{k}\right){(1-{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta })}^k{(1-e^{-{\left(\rho x\right)}^{-2}})}^{\eta (\nu +\tau -1-k)}\right\}}^{r-1+j}.\hfill \end{array}$$

5. Maximum Likelihood Estimation Method

In this section, we determine the maximum likelihood estimates (MLEs) of the parameters of the BGIRD from complete samples. Let $X_1,X_2,\cdots ,X_n$ be a random sample of size n from the BGIRD $\underline{\theta }$= ($\nu$, $\tau$, $\eta$, $\rho$).

The likelihood function in this case can be written as:

$$L(\underline{\theta },\underline{x})={\displaystyle \prod _{i=1}^n}f\left(x_i\right)$$

(27)

The likelihood function for the BGIRD is given by:

$$\begin{array}{cc}\hfill L(\underline{\theta },\underline{x})& ={\left(\frac{2\eta }{\beta (\nu ,\tau ){\rho }^2}\right)}^n{\displaystyle \prod _{i=1}^n}{\left(x_i\right)}^{-3}{e}^{-{\sum }_{i=1}^n{\left(\rho x_i\right)}^{-2}}{\displaystyle \prod _{i=1}^n}{(1-e^{-{\left(\rho x_i\right)}^{-2}})}^{\eta \tau -1}\hfill \\ \hfill & \times {\displaystyle \prod _{i=1}^n}{\left[1-{\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right)}^{\eta }\right]}^{\nu -1}\hfill \end{array}$$

(28)

and the log-likelihood function is obtained as:

$$\ell =logL(\underline{\theta },\underline{x})={\displaystyle \sum _{i=1}^n}logf\left(x_i\right).$$

(29)

For the BGIRD, we have

$$\begin{array}{cc}\hfill logL(\underline{\theta },\underline{x})& =nlog\left(\eta \right)+nlog\left(2\right)-nlog\left(\beta (\nu ,\tau )\right)-2nlog\left(\rho \right)-3{\displaystyle \sum _{i=1}^n}log\left(x_i\right)-{\displaystyle \sum _{i=1}^n}{\left(\rho x_i\right)}^{-2}\hfill \\ \hfill & +(\eta \tau -1){\displaystyle \sum _{i=1}^n}log(1-e^{-{\left(\rho x_i\right)}^{-2}})+(\nu -1){\displaystyle \sum _{i=1}^n}log[1-{(1-e^{-{\left(\rho x_i\right)}^{-2}})}^{\eta }].\hfill \end{array}$$

(30)

The first derivatives of the log-likelihood function with respect to the components of $\underline{\theta }$ are given by

$$\frac{\partial logL}{\partial \nu }=\frac{-n}{\beta (\nu ,\tau )}{\varphi }_1+{\displaystyle \sum _{i=1}^n}log\left[1-{\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right)}^{\eta }\right],$$

$$\begin{array}{cc}\hfill {\varphi }_1& =\frac{\partial \beta (\nu ,\tau )}{\partial \nu }=\frac{\Gamma \left(\tau \right)\left[\Gamma (\nu +\tau ){\Gamma }^{\prime }\left(\nu \right)-\Gamma \nu \frac{\partial \Gamma (\nu +\tau )}{\partial \nu }\right]}{{\left[\Gamma (\nu +\tau )\right]}^2}\hfill \\ & =\beta (\nu ,\tau )\left[\psi \left(\nu \right)-\psi (\nu +\tau )\right],\hfill \end{array}$$

where $\psi \left(x\right)=\frac{d[ln\Gamma (x\left)\right]}{dx}=\frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}$ is called the Psi (Digamma) Function [16].

$$\begin{array}{c}\hfill \frac{\partial logL}{\partial \nu }=-n\left[\psi \left(\nu \right)-\psi (\nu +\tau )\right]+{\displaystyle \sum _{i=1}^n}log\left[1-{\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right)}^{\eta }\right].\\ \hfill \frac{\partial logL}{\partial \tau }=\frac{-n}{\beta (\nu ,\tau )}{\varphi }_2+\eta {\displaystyle \sum _{i=1}^n}log\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right),\end{array}$$

(31)

where

$$\begin{array}{cc}\hfill {\varphi }_2& =\frac{\partial \beta (\nu ,\tau )}{\partial \tau }=\frac{\Gamma \left(\nu \right)\left[\Gamma (\nu +\tau ){\Gamma }^{\prime }\left(\tau \right)-\Gamma \tau \frac{\partial \Gamma (\nu +\tau )}{\partial \tau }\right]}{{\left[\Gamma (\nu +\tau )\right]}^2}\hfill \\ & =\beta (\nu ,\tau )\left[\psi \left(\tau \right)-\psi (\nu +\tau )\right],\hfill \end{array}$$

$$\frac{\partial logL}{\partial \tau }=-n\left[\psi \left(\tau \right)-\psi (\nu +\tau )\right]+\eta {\displaystyle \sum _{i=1}^n}log\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right).$$

(32)

$$\frac{\partial logL}{\partial \eta }=\frac{n}{\eta }+{\displaystyle \sum _{i=1}^n}log(1-e^{-{\left(\rho x_i\right)}^{-2}})\left\{\tau -(\nu -1){\left[{(1-e^{-{\left(\rho x_i\right)}^{-2}})}^{-\eta }-1\right]}^{-1}\right\}.$$

(33)

$$\begin{array}{cc}\hfill \frac{\partial logL}{\partial \rho }& =\frac{-2n}{\rho }+\frac{2{\sum }_{i=1}^n{\left(x_i\right)}^{-2}}{{\rho }^3}-\frac{2{\sum }_{i=1}^n{\left(x_i\right)}^{-2}}{{\rho }^3}{\left[e^{{\left(\rho x_i\right)}^{-2}}-1\right]}^{-1}\hfill \\ & \times \left\{\left(\eta \tau -1\right)-\eta \left(\nu -1\right){\left[{\left(1-e^{-{\left(\rho x_i\right)}^{-2}}\right)}^{-\eta }\right]}^{-1}\right\}.\hfill \end{array}$$

(34)

Setting the four non-linear Equations (31)–(34) to zero and solving the resulting system of non-linear equations, we obtain the maximum likelihood estimators of the unknown parameters $\eta$, $\rho$, $\nu$, and $\tau$ of BGIRD. These equations are in implicit form, so they may be solved using numerical iteration, such as the Newton–Raphson method via Mathematica 12.0.

6. Simulation Study

In this part, a Monte Carlo simulation is investigated for estimating unknown parameters and the reliability function and hazard rate function of BGIRD. The simulation is conducted by using Mathematica 12.0, 1000 random samples of BGIRD were generated with values of n = (10, 20, 30, and 50) while choosing ($\eta$, $\rho$, $\nu$, $\tau$) = (0.5, 2, 2, 2). Average of the absolute relative bias (ARBias) and mean square error (MSE), where

$$ARBias\left(\widehat{\theta }\right)=\left|\frac{\widehat{\theta }-\theta }{\theta }\right|$$

$$MSE\left(\widehat{\theta }\right)=E{(\widehat{\theta }-\theta )}^2$$

The results are summarized and tabulated in

Table 2. Maximum likelihood estimates, ARBias, and MSE of the parameters $\underline{\theta }=(\nu ,\tau ,\eta ,\rho )$ = (2, 2, 0.5, 2) and R($t_0$) = 0.1519, h($t_0$) = 0.8740 at t${}_0$ = 2.

n		$\widehat{\mathit{\nu }}$	$\widehat{\mathit{\tau }}$	$\widehat{\mathit{\eta }}$	$\widehat{\mathit{\rho }}$	$\widehat{\mathit{R}}\left({\mathit{t}}_0\right)$	$\widehat{\mathit{h}}\left({\mathit{t}}_0\right)$
	MLEs	1.8771	2.1164	0.5575	1.9080	0.1219	1.0487
10	ARBias	0.0614	0.0582	0.1149	0.0460	0.1977	0.1999
	MSE	0.1384	0.1120	0.0204	0.0963	0.0056	0.1160
	MLEs	1.9092	2.0567	0.5353	1.9394	0.1348	0.9684
20	ARBias	0.0454	0.0283	0.0706	0.0303	0.1130	0.1080
	MSE	0.1217	0.0846	0.0133	0.0812	0.0032	0.0576
	MLEs	1.9511	2.0456	0.5190	1.9713	0.1412	0.9281
30	ARBias	0.0244	0.0228	0.0381	0.0143	0.0708	0.0619
	MSE	0.1191	0.0830	0.0064	0.0761	0.0019	0.0283
	MLEs	1.9430	2.0373	0.5131	1.9681	0.1447	0.9104
50	ARBias	0.0285	0.0187	0.0261	0.0159	0.0477	0.0417
	MSE	0.1160	0.0771	0.0057	0.0730	0.0014	0.0182
	MLEs	1.9714	2.0082	0.5049	1.9879	0.1490	0.8857
100	ARBias	0.0143	0.0041	0.0098	0.0060	0.0191	0.0135
	MSE	0.0373	0.0245	0.0019	0.0247	0.0005	0.0052

which contain the values of the ARBias and MSEs for estimating the parameters and the reliability function and hazard rate function of BGIRD.

Results of ML estimates for the four parameters, $R\left(t_0\right)$ and $h\left(t_0\right)$:

• From Table 2, we note that the MSEs of the ML estimates for BGIR($\nu$, $\tau$, $\eta$, $\rho$), $R\left(t_0\right)$, and $h\left(t_0\right)$ decrease as the sample size increases which show consistency of the estimated parameters.
• According to the simulation results given in Table 2, as the sample size n increases, the ARBias is close to zero, the mean estimates tend to be closer to the true parameter values.

7. Application

In this section, we provide four real data sets to illustrate the importance and flexibility of the BGIRD. The data sets that we applied on our model have been analyzed by [17] to assess the exbility of exponential, Lindley, and Akash distributions.

We shall compare the fit of the proposed BGIRD (and its sub-model namely: BIRD) with several other competitive models namely: the beta generalized inverse Weibull distribution (BGIWD) [18], and exponential, Lindley, and Akash distributions [17] with corresponding cumulative distribution functions (for x > 0) as follows

$$BGIWD:F\left(x\right)=\frac{1}{\beta (\nu ,\tau )}{\int }_0^{e^{-\theta {\left(\frac{\rho }{x}\right)}^{\eta }}}{Z}^{\nu -1}{(1-Z)}^{\tau -1}dZ,\nu ,\tau ,\rho ,\theta and\eta >0.$$

$$exponential:F\left(x\right)=1-e^{-\rho t}$$

$$Lindley:F\left(x\right)=1-\frac{\rho +1+\rho x}{\rho +1}{e}^{-\rho x}$$

$$Akash:F\left(x\right)=1-\left[1+\frac{\rho x(\rho x+2)}{{\rho }^2+2}\right]e^{-\rho t}$$

The previous models are chosen to compare fitting the selected data sets with the BGIRD. These models are distinguished by two distinct groups: beta-G distributions (BIRD and BGIWD) and one parameter distributions (exponential, Lindley, and Akash). Theoretically, our model has complex derivations but in the lifetime context it has a flexible failure rate function which has an upside down bathtub shape. On the other hand the exponential distribution has a constant failure rate function. Therefore, our model is better than the exponential distribution at this point.

In order to compare the models, we consider some goodness-of-fit measures including $-2\widehat{\ell }$, Akaike information criterion (AIC), and Bayesian information criterion (BIC),

$$AIC=2k-2\widehat{\ell }$$

and

$$BIC=kLog\left(n\right)-2\widehat{\ell }$$

where $\widehat{\ell }$ is the maximized log-likelihood, k is the number of parameters, and n is the sample size. For more discussion on these criteria, see [19]. In general, the model with minimum values for these statistics could be chosen as the best model to fit the data. The required numerical evaluations are carried out using Mathematica 12.0.

Data Set 1: This data represents the survival times of a group of patients suffering from Head and Neck cancer disease and treated using a combination of radiotherapy and chemotherapy. These data are: 12.2, 23.56, 23.74, 25.87, 31.98, 37, 41.35, 47.38, 55.46, 58.36, 63.47, 68.46, 78.26, 74.47, 81.43, 84, 92, 94, 110, 112, 119, 127, 130, 133, 140, 146, 155, 159, 173, 179, 194, 195, 209, 249, 281, 319, 339, 432, 469, 519, 633, 725, 817, 1776.

Data Set 2: The data set represents the failure times of the air conditioning system of an airplane. It has 30 observations as follows: 23, 261, 87, 7, 120, 14, 62, 47, 225, 71, 246, 21, 42, 20, 5, 12, 120, 11, 3, 14, 71, 11, 14, 11, 16, 90, 1, 16, 52, 95.

Data Set 3: The data set represents the lifetime data relating to relief times (in minutes) of patients receiving an analgesic. The data set consists of 20 observations and it is as follows: 1.1, 1.4, 1.3, 1.7, 1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3.0, 1.7, 2.3, 1.6, 2.0.

Data Set 4: These data are the strengths of 1.5 cm glass fibers, measured at the National Physical Laboratory, England. The data set is: 0.55, 0.93, 1.25, 1.36, 1.49, 1.52, 1.58, 1.61, 1.64, 1.68, 1.73, 1.81, 2, 0.74, 1.04, 1.27, 1.39, 1.49, 1.53, 1.59, 1.61, 1.66, 1.68, 1.76, 1.82, 2.01, 0.77, 1.11, 1.28, 1.42, 1.5, 1.54, 1.6, 1.62, 1.66, 1.69, 1.76, 1.84, 2.24, 0.81, 1.13, 1.29, 1.48, 1.5, 1.55, 1.61, 1.62, 1.66, 1.7, 1.77, 1.84, 0.84, 1.24, 1.3, 1.48, 1.51, 1.55, 1.61, 1.63, 1.67, 1.7, 1.78, 1.89.

Data Set 5: These data arose in testing on the cycle at which the Yarn failed. The data are the number of cycles until failure of the yarn and they are: 86, 146, 251, 653, 98, 249, 400, 292, 131, 169, 175, 176, 76, 264, 15, 364, 195, 262, 88, 264, 157, 220, 42, 321, 180, 198, 38, 20, 61, 121, 282, 224, 149, 180, 325, 250, 196, 90, 229, 166, 38, 337, 65, 151, 341, 40, 40, 135, 597, 246, 211, 180, 93, 315, 353, 571, 124, 279, 81, 186, 497, 182, 423, 185, 229, 400, 338, 290, 398, 71, 246, 185, 188, 568, 55, 55, 61, 244,20, 284, 393, 396, 203, 829, 239, 236, 286, 194, 277, 143, 198, 264, 105, 203, 124, 137, 135, 350, 193, 188.

The applicability of the BGIRD was demonstrated with four real lifetime data sets. The data sets were also fitted to other distributions,

Table 3. Parameter estimates, goodness-of-fit measures of the fitted distributions of Data Sets 1–4.

Data	Model	MLEs					Statistics
		$\widehat{\mathit{\rho }}$	$\widehat{\mathit{\nu }}$	$\widehat{\mathit{\tau }}$	$\widehat{\mathit{\eta }}$	$\widehat{\mathit{\theta }}$	AIC	BIC	−2 $\mathit{logL}$
Data 1	BGIRD	14.0184	74.0009	8.71669	0.15366		563.039	570.175	555.039
	BIRD	6899.62	4900.34	0.0543026			715.539	720.892	709.539
	BGIWD	4938.59	0.001181	0.0862678	0.750759	39.904	576.569	585.49	566.569
	Exponential	0.004475					566.02	567.80	564.02
	Lindley	0.00891					581.16	582.95	579.16
	Akash	0.013423					611.93	613.71	609.93
Data 2	BGIRD	5.0833	14.4244	2.9408	0.1921		316.199	321.80	308.199
	BIRD	0.5395	9.6731	0.2074			330.41	333.21	326.41
	BGIWD	4227.66	0.0013	0.0586	0.534319	40.7341	326.17	333.18	316.17
	Exponential	0.016779					307.26	308.66	305.26
	Lindley	0.033021					325.2	326.6	323.2
	Akash	0.050293					356.88	358.2	354.8
Data 3	BGIRD	17.5039	58.042	6.3644	0.3428		51.067	55.05	43.067
	BIRD	7569.7	5201.94	0.098989			133.896	136.883	127.896
	BGIWD	2464.11	0.0015	0.0396	0.4256	29.545	109.587	114.566	99.5871
	Exponential	0.526316					67.67	68.67	65.67
	Lindley	0.816118					60.50	62.50	63.49
	Akash	1.156923					61.52	62.51	59.52
Data 4	BGIRD	21.0277	69.3754	9.45318	0.3144		102.438	111.01	94.4377
	BIRD	0.65579	0.808613	1.00416			115.251	121.681	109.251
	BGIWD	37.7781	0.0052	0.174197	1.28585	2.63091	164.39	175.106	154.39
	Exponential	0.663647					179.66	181.80	177.66
	Lindley	0.996116					164.56	166.70	162.56
	Akash	1.355445					165.73	169.93	163.73
Data 5	BGIRD	12.0129	155.906	7.49802	0.205185		1285.79	1296.21	1277.79
	BIRD	19758	33641.7	0.0520308			1690.43	1698.24	1684.43
	BGIWD	5.54157	1.93615	3.77834	0.629758	9.71891	1289.55	1302.57	1279.55
	Exponential	0.004505					1282.52	1285.12	1280.52
	Lindley	0.00897					1253.34	1255.95	1251.34
	Akash	0.013514					1257.83	1260.43	1255.83

. The method of maximum likelihood estimation was used to estimate the parameters of the distributions. It is evident from Data Sets 3 and 4, the BGIRD has the lowest AIC, BIC, and −2 $logL$ values among all fitted models. Hence, this new distribution can be chosen as the best model for fitting these data sets. From Data Set 1, there is a strong competition between BGIRD and the exponential distribution. Exponential distribution is the best according to BIC but when compared to AIC value and −2 $logL$ value BGIRD is the best. In addition, as for the second data set, the exponential distribution gives better fit and the next distribution is BGIRD. As for the fifth and final data set, the most appropriate distribution was the Lindley distribution and our model is considered the third competitor distribution for fitting these data.

Figure 5 shows the curves of the empirical cdf of the data sets and the fitted cdfs.

8. Conclusions

We studied the beta generalized inverse Rayleigh distribution (BGIRD). This new unimodal distribution has a positively skewed curve for all values of the parameters. A comprehensive study of the statistical properties of the proposed distribution has been provided. The quantile function, median, mode, moments, harmonic mean, skewness, kurtosis, and the mean deviation from the mean and from the median have been obtained. Moreover, we have derived the Rényi and Shannon entropies. In addition, we have obtained the order statistics. We hope that the proposed model may be interesting for a wider range of statistical research. All results in this article generalize the generalized inverted Rayleigh distribution (GIRD) discussed by [5]. Maximum likelihood estimators of the BGIRD parameters are obtained. Simulation studies of Monte Carlo are conducted under various sample sizes to study the theoretical performance of the MLE of the parameters. Four real data sets are analyzed and a good fit for the data sets has been provided by the BGIRD.