A Novel One Parameter Unit Distribution: Median Based Unit Rayleigh: Properties and Estimations

2.1. Basic Functions (PDF, CDF, Sf, HR, rHR)

The following are the probability density function (PDF), cumulative distribution function (CDF), survival function, hazard function, and reversed hazard function as shown in equation (5-9) respectively:

$$f\left(y\right)={\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)},0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(5)

$$F\left(y\right)={3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}},0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(6)

$$S\left(y\right)=1-F\left(Y\right)=1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right),0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(7)

$$h\left(y\right)={\displaystyle \frac{f\left(y\right)}{S\left(y\right)}}={\displaystyle \frac{\frac{6}{{\alpha }^2}\left(1-y^{\frac{1}{{\alpha }^2}}\right)y^{\left(\frac{2}{{\alpha }^2}-1\right)}}{1-\left({3y}^{\frac{2}{{\alpha }^2}}-{2y}^{\frac{3}{{\alpha }^2}}\right)}},0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(8)

$$rh\left(y\right)={\displaystyle \frac{f\left(y\right)}{F\left(y\right)}}={\displaystyle \frac{\frac{6}{{\alpha }^2}\left(1-y^{\frac{1}{{\alpha }^2}}\right)y^{\left(\frac{2}{{\alpha }^2}-1\right)}}{{3y}^{\frac{2}{{\alpha }^2}}-{2y}^{\frac{3}{{\alpha }^2}}}},0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(9)

2.2. Quantile Function

The CDF of the MBUR can be written as a third degree polynomial

$$u=F\left(y\right)={3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}=-2{\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^3+3{\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^2$$

The inverse of this CDF is used to obtain y, the real root of this 3rd polynomial function is $y=F^{-1}\left(y\right)$ as shown in (equation 10):

$$y={\left\{-.5\left(cos\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]-\sqrt{3}sin\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]\right)+.5\right\}}^{{\alpha }^2}\dots \dots \dots \dots \dots \dots \dots \dots ..$$

(10)

Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 illustrate the specified functions for various values of alpha for a random variable X distributed as MBUR.

To generate random variable distributed as MBUR:

• Generate uniform random variable (0,1): $u~uniform\left(\mathrm{0,1}\right)$
• Choose the parameter alpha.
• Substitute the above values of $u\in \left(\mathrm{0,1}\right)$and the chosen alpha in the quantile function, to obtain x distributed as $x~MBUR\left(\alpha \right)$

2.3. Discussion and Analysis of the Above Functions

See supplementary materials (section 1).

2.4. rth Raw Moments for a Random Variable y Distributed as MBUR is Shown in Equation (11)

$$E\left(y^r\right)={\int }_0^1{y}^r{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy={\displaystyle \frac{6}{\left(2+r{\alpha }^2\right)\left(3+r{\alpha }^2\right)}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(11)

$$E\left(y\right)={\displaystyle \frac{6}{\left(2+{\alpha }^2\right)\left(3+{\alpha }^2\right)}},E\left(y^2\right)={\displaystyle \frac{6}{\left(2+2{\alpha }^2\right)\left(3+2{\alpha }^2\right)}}$$

$$E\left(y^3\right)={\displaystyle \frac{6}{\left(2+3{\alpha }^2\right)\left(3+3{\alpha }^2\right)}},E\left(y^4\right)={\displaystyle \frac{6}{\left(2+4{\alpha }^2\right)\left(3+4{\alpha }^2\right)}}$$

$$\mathrm{var}\left(y\right)=E\left(y^2\right)-{\left[E\left(y\right)\right]}^2={\displaystyle \frac{3{\alpha }^4\left(13+10{\alpha }^2+{\alpha }^4\right)}{\left(3+5{\alpha }^2+2{\alpha }^4\right){\left(6+5{\alpha }^2+{\alpha }^4\right)}^2}}$$

Coefficient of Skewness: is shown in equation (12)

$$E{\displaystyle \frac{{\left(y-\mu \right)}^3}{{\sigma }^3}}={\displaystyle \frac{E\left(y^3\right)-3\mu E\left(y^2\right)+3{\mu }^{2}E\left(y\right)-{\mu }^3}{{\sigma }^3}}$$

$$={\displaystyle \frac{E\left(y^3\right)-3\mu \left[E\left(y^2\right)-\mu E\left(y\right)\right]-{\mu }^3}{{\sigma }^3}}={\displaystyle \frac{E\left(y^3\right)-3\mu \left[E\left(y^2\right)-\mu \mu \right]-{\mu }^3}{{\sigma }^3}}$$

$$\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{f}\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{w}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}={\displaystyle \frac{E\left(y^3\right)-3\mu {\sigma }^2-{\mu }^3}{{\sigma }^3}}={\displaystyle \frac{E\left(y^3\right)-\mu \left({3\sigma }^2+{\mu }^2\right)}{{\sigma }^3}}\dots \dots \dots \dots \dots \dots .$$

(12)

Coefficient of Kurtosis: is shown in equation (13)

$$\begin{gathered} E{\displaystyle \frac{{\left(y-\mu \right)}^4}{{\sigma }^4}}={\displaystyle \frac{E\left(y^4\right)-4\mu E\left(y^3\right)+6{\mu }^{2}E\left(y^2\right)-{3\mu }^4}{{\sigma }^4}} \\ ={\displaystyle \frac{E\left(y^4\right)-4\mu E\left(y^3\right)+6{\mu }^2\left[{\sigma }^2+{\mu }^2\right]-{3\mu }^4}{{\sigma }^4}} \\ ={\displaystyle \frac{E\left(y^4\right)-4\mu E\left(y^3\right)+6{\mu }^2{\sigma }^2+6{\mu }^4-{3\mu }^4}{{\sigma }^4}} \\ ={\displaystyle \frac{E\left(y^4\right)-4\mu E\left(y^3\right)+6{\mu }^2{\sigma }^2+3{\mu }^4}{{\sigma }^4}} \end{gathered}$$

$$\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{f}\mathrm{K}\mathrm{u}\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{s}={\displaystyle \frac{E\left(y^4\right)-4\mu E\left(y^3\right)+3{\mu }^2\left[{2\sigma }^2+{\mu }^2\right]}{{\sigma }^4}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(13)

Coefficient of Variation: is shown in see equation (14)

$$\mathrm{C}\mathrm{V}={\displaystyle \frac{S}{\mu }}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(14)

where S is standard deviation and mu is the mean

Figure 10 illustrates the graph for the specified coefficients. In subplot (1), the variance increases as the alpha level rises, reaching a maximum value between 0.05 and 0.06 at a parameter level between 1 and 1.5. After that point, the variance begins to decrease with further increases in the parameter values.

In subplot (2), the Fisher coefficient of skewness shows negative skewness (left skewness) at low parameter levels, reaching a zero value when the parameter equals 1. At this point, the probability density function (PDF) of the distribution becomes symmetrical around 0.5. Beyond this parameter value of 1, the coefficient of skewness increases as the alpha level rises, indicating right skewness in the distribution.

Subplot (3) reveals that the coefficient of kurtosis starts at approximately 6 when the alpha level is 0.1, decreasing to around 2 at an alpha level of 1. Following this, the coefficient rises again as the alpha level increases, demonstrating a wide variety of kurtosis shapes. The distribution exhibits a mesokurtic shape (kurtosis equals 3) at approximately parameter level of 0.7 and 1.5. It displays a leptokurtic shape, characterized by fatter tails, when there is positive excess kurtosis (kurtosis greater than 3), which occurs at parameter levels below 0.7 and above 1.5. Conversely, it shows a platykurtic shape, with thinner tails, when there is negative excess kurtosis (kurtosis less than 3), which happens approximately at parameter levels between 0.7 and 1.5. These coefficients, which reflect the shape of the PDF, are fundamental to the flexibility of the new distribution, allowing it to accommodate a wide variety of data shapes. They play a crucial role in making this distribution outperform other unit distributions in certain data analyses.

2.5. rth Incomplete Moments: for a Random Variable y Distributed as MBUR is Defined in Equation (15): (See Supplementary Materials Section1 for Derivation)

$$E\left(y^r|y<t\right)={\int }_0^t{y}^r{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy={\displaystyle \frac{6y^{\frac{2}{{\alpha }^2}+r}}{\left(2+r{\alpha }^2\right)}}-{\displaystyle \frac{6y^{\frac{3}{{\alpha }^2}+r}}{\left(3+r{\alpha }^2\right)}}\dots \dots \dots \dots \dots .$$

(15)

2.6. Stress- Strength Reliability

From a reliability prospective, if an element in a system has a random strength X that is strained with a random stress Y, this element will immediately break down if the stress overrides the strength. Conversely, it will adequately operate if the strength surpasses the stress. If X and Y are independent random variables denoting strength and stress respectively, and both follow MBUR distribution with parameters ${\alpha }_1$and ${\alpha }_2$ respectively, then the reliability measure of this element can be deduced from appropriate equation (16) .(see suppl. Mat. Section1)

$$R=Pr\left(Y<X\right)={\int }_0^{1}Pr\left(Y<X|X=x\right)f_X\left(x\right)dy$$

$$R=Pr\left(Y<X\right)={\int }_0^{1}F\left(x;{\alpha }_2\right)f\left(x;{\alpha }_1\right)dy$$

$$R={\alpha }_2^2\left\{{\displaystyle \frac{13}{\left({\alpha }_1^2+{\alpha }_2^2\right)}}-{\displaystyle \frac{18}{\left(2{\alpha }_1^2+{3\alpha }_2^2\right)}}-{\displaystyle \frac{12}{\left(3{\alpha }_1^2+{2\alpha }_2^2\right)}}\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(16)

2.7. Lorenz, Bonferroni curves and Gini Index

These indices have many applications in medicine, insurance, demography, and economics for studying wealth and poverty. They can be applied to variables defined as proportions where y is a random variable distributed as MBUR. Lorenz curve, Bonferroni curves and Gini index are defined in equation (17), (18), (19) respectively (see supplementary materials section 1 for derivation)

$$L\left(p\right)={\displaystyle \frac{{\int }_0^{t}yf\left(y\right)dy}{{\int }_0^{1}yf\left(y\right)dy}}=y^{{\displaystyle \frac{2}{{\alpha }^2}}+1}\left\{\left(3+{\alpha }^2\right)-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\left(2+{\alpha }^2\right)\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(17)

$$B\left(p\right)={\displaystyle \frac{L\left(p\right)}{F\left(y\right)}}={\displaystyle \frac{y^{\frac{2}{{\alpha }^2}+1}\left\{\left(3+{\alpha }^2\right)-y^{\frac{1}{{\alpha }^2}}\left(2+{\alpha }^2\right)\right\}}{{3y}^{\frac{2}{{\alpha }^2}}-{2y}^{\frac{3}{{\alpha }^2}}}}$$

$$B\left(p\right)={\displaystyle \frac{L\left(p\right)}{F\left(y\right)}}={\displaystyle \frac{y\left\{\left(3+{\alpha }^2\right)-y^{\frac{1}{{\alpha }^2}}\left(2+{\alpha }^2\right)\right\}}{3-{2y}^{\frac{1}{{\alpha }^2}}}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(18)

$$\mathrm{G}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=1-2{\int }_0^{1}L\left(p\right)dp=1-\left({\displaystyle \frac{{\alpha }^2\left(5+3{\alpha }^2\right)}{\left(1+{\alpha }^2\right)\left(3+2{\alpha }^2\right)}}\right)\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(19)

2.8. Renyi Entropy

Entropy quantifies the variation in uncertainty within a random variable, in the paper context, y is distributed as MBUR. Renyi entropy $T_R\left(\gamma \right)$is a well-known measure defined as follows in equation (20). For MBUR it is defined in equation (21):

$$T_R\left(\gamma \right)={\displaystyle \frac{1}{1-\gamma }}\log \left\{{\int }_0^1{\left[f\left(y\right)\right]}^{\gamma }dy\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(20)

$$T_R\left(\gamma \right)={\displaystyle \frac{1}{1-\gamma }}\log \left\{{\int }_0^1{\left[\left({\displaystyle \frac{6}{{\alpha }^2}}\right)\left(1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\right]}^{\gamma }dy\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(21)

Expanding the following term using the binomial expansion as in equation (22):

$${\left(1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^{\gamma }={\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{\gamma }{m}}\right){\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(22)

substitute equation (22) into equation (21) gives equations (23), (24) & (25).

$$={\displaystyle \frac{1}{1-\gamma }}\log \left\{{\int }_0^1\left[{\left({\displaystyle \frac{6}{{\alpha }^2}}\right)}^{\gamma }{\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{\gamma }{m}}\right){\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m{y}^{\left({\displaystyle \frac{2\gamma -\gamma {\alpha }^2}{{\alpha }^2}}\right)}\right]dy\right\}$$

The integral will pass to the variable for integration:

$$={\displaystyle \frac{1}{1-\gamma }}\log \left\{{\left({\displaystyle \frac{6}{{\alpha }^2}}\right)}^{\gamma }{\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{\gamma }{m}}\right){\int }_0^1\left[{\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m{y}^{\left({\displaystyle \frac{2\gamma -\gamma {\alpha }^2}{{\alpha }^2}}\right)}\right]dy\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(23)

$$={\displaystyle \frac{1}{1-\gamma }}\log \left\{{\left({\displaystyle \frac{6}{{\alpha }^2}}\right)}^{\gamma }{\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{\gamma }{m}}\right)\left({\displaystyle \frac{{\alpha }^2}{m+2\gamma -\gamma {\alpha }^2+{\alpha }^2}}\right)\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(24)

$$={\displaystyle \frac{1}{1-\gamma }}\log \left\{{6^{\gamma }\alpha }^{-2\left(\gamma -1\right)}{\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{\gamma }{m}}\right){\left(m+2\gamma -\gamma {\alpha }^2+{\alpha }^2\right)}^{-1}\right\}\dots \dots \dots \dots \dots \dots \dots \dots ..$$

(25)

2.9. Mean Residual Life Function: See Figure 11

It is defined for the MBUR random variable as shown in equation (26). Taking the limits at the ends of the unit interval is shown in equation (27)

$$\mu \left(y|\alpha \right)=E\left(Y-y|Y>y\right)={\displaystyle \frac{1}{S\left(y\right)}}{\int }_y^{1}S\left(y\right)dy$$

$$={\displaystyle \frac{\left(1-\frac{3{\alpha }^2}{2+{\alpha }^2}+\frac{2{\alpha }^2}{3+{\alpha }^2}\right)-y\left\{1-\frac{3{\alpha }^2{y}^{\frac{2}{{\alpha }^2}}}{2+{\alpha }^2}+\frac{3{\alpha }^2{y}^{\frac{2}{{\alpha }^2}}}{3+{\alpha }^2}\right\}}{{1-3y}^{\frac{2}{{\alpha }^2}}+{2y}^{\frac{3}{{\alpha }^2}}}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(26)

$$\underset{y\to 0}{\lim }\mu \left(y|\alpha \right)=1-{\displaystyle \frac{3{\alpha }^2}{2+{\alpha }^2}}+{\displaystyle \frac{2{\alpha }^2}{3+{\alpha }^2}}\mathrm{and}\underset{y\to 1}{\lim }\mu \left(y|\alpha \right)=0$$

(27)

2.10. Stochastic Ordering: See Figure 12

This ordering judges the comparative conduct of a variable. A random variable X is considered smaller than the random variable Y in the following orders:

• Stochastic order $\left(X{\le }_{st}Y\right)$ if $F_X\left(x\right)\ge F_Y\left(x\right)$ for all x.
• Hazard rate order $\left(X{\le }_{hr}Y\right)$ if $h_X\left(x\right)\ge h_Y\left(x\right)$ for all x.
• Mean residual life order $\left(X{\le }_{mrl}Y\right)$ if $m_X\left(x\right)\le m_Y\left(x\right)$ for all x.
• Likelihood ratio order $\left(X{\le }_{lr}Y\right)$ if ${\displaystyle \frac{f_X\left(x\right)}{f_Y\left(x\right)}}$ decreases in x.

The following results are due to Shaked and Shanthikumar (23). They used the results to evaluate the stochastic ordering of a distribution:

$X{\le }_{lr}Y\Rightarrow X{\le }_{hr}Y\Rightarrow X{\le }_{mrl}Yhence:X{\le }_{st}Y$

Theorem: let $X~MBUR\left({\alpha }_1\right)$ and $Y~MBUR\left({\alpha }_2\right)$ if ${\alpha }_1<{\alpha }_2$ , then

$X{\le }_{lr}Y$ , hence $X{\le }_{hr}Y$ , $X{\le }_{mrl}Y$ and $X{\le }_{st}Y$ .

Proof: see equation (28)

$${\displaystyle \frac{f_X\left(y;{\alpha }_1\right)}{f_Y\left(y;{\alpha }_2\right)}}={\displaystyle \frac{\frac{6}{{\alpha }_1^2}\left[1-y^{\frac{1}{{\alpha }_1^2}}\right]y^{\left(\frac{2}{{\alpha }_1^2}-1\right)}}{\frac{6}{{\alpha }_2^2}\left[1-y^{\frac{1}{{\alpha }_2^2}}\right]y^{\left(\frac{2}{{\alpha }_2^2}-1\right)}}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(28)

Taking the log of equation (28) gives equation (29)

$$\mathrm{l}\mathrm{n}\left({\displaystyle \frac{6}{{\alpha }_1^2}}\right)+ln\left(1-y^{{\displaystyle \frac{1}{{\alpha }_1^2}}}\right)+ln\left[y^{\left({\displaystyle \frac{2}{{\alpha }_1^2}}-1\right)}\right]-ln\left({\displaystyle \frac{6}{{\alpha }_2^2}}\right)-ln\left[1-y^{{\displaystyle \frac{1}{{\alpha }_2^2}}}\right]-ln\left[y^{\left({\displaystyle \frac{2}{{\alpha }_2^2}}-1\right)}\right]\dots \dots .$$

(29)

Taking the first derivative of the likelihood ratio order with respect to the variable y as shown in equation (30):

$${\displaystyle \frac{dLRO}{dy}}={\displaystyle \frac{2}{y}}\left\{{\displaystyle \frac{1}{{\alpha }_1^2}}-{\displaystyle \frac{1}{{\alpha }_2^2}}\right\}+{\displaystyle \frac{1}{y}}\left\{{\displaystyle \frac{\frac{-1}{{\alpha }_1^2}{y}^{\frac{1}{{\alpha }_1^2}}}{1-y^{\frac{1}{{\alpha }_1^2}}}}+{\displaystyle \frac{\frac{1}{{\alpha }_2^2}{y}^{\frac{1}{{\alpha }_2^2}}}{1-y^{\frac{1}{{\alpha }_2^2}}}}\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(30)

Figure 12 describes the behavior of the previous function defined in equation (30). Taking the limit at the ends of the unit interval gives equation (31).

$$\underset{y\to 1}{\lim }LRO=3\left\{{\displaystyle \frac{1}{{\alpha }_1^2}}-{\displaystyle \frac{1}{{\alpha }_2^2}}\right\}\mathrm{and}\underset{y\to 0}{\lim }LRO=+\infty$$

(31)

See supplementary materials (section1) for derivation of equation 31.( the Likelihood Ratio Order (LRO))

2.11. Probability Weighted Moments

Probability weighted moments are less vulnerable to extreme values and are tractable to obtain when ML estimators struggle for estimation. This can be deduced from the following equations (32) & (33):

$$E\left(y^r{\left[F\left(y\right)\right]}^s\right)={\int }_0^1{y}^r{\left[F\left(y\right)\right]}^{s}f\left(y\right)dy\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(32)

$$\begin{gathered} ={\int }_0^1{y}^r{\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^s{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy \\ ={\int }_0^1{y}^r{\left({3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}\right)}^s{\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]}^s{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy \end{gathered}$$

$$={\int }_0^1{y}^r\left({3^{s}y}^{{\displaystyle \frac{2s}{{\alpha }^2}}}\right){\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]}^s{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(33)

Binomial expansion of ${\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]}^s$

is shown in equation (34) :

$${\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]}^s={\sum }_{m=0}^{\infty }{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(34)

substitute equation (34) into equation (33).

$$\begin{gathered} ={\displaystyle \frac{3^{s}6}{{\alpha }^2}}{\int }_0^1{y}^r\left(y^{{\displaystyle \frac{2s}{{\alpha }^2}}}\right){\sum }_{m=0}^s{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy \\ ={\displaystyle \frac{3^{s}6}{{\alpha }^2}}{\sum }_{m=0}^s{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({\displaystyle \frac{2}{3}}\right)}^m{\int }_0^1{y}^r\left(y^{{\displaystyle \frac{2s}{{\alpha }^2}}}\right){\left(y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}dy \\ ={\displaystyle \frac{3^{s}6}{{\alpha }^2}}{\sum }_{m=0}^s{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({\displaystyle \frac{2}{3}}\right)}^m\left\{{\displaystyle \frac{{\alpha }^2}{\left(r{\alpha }^2+2s+m+2\right)}}-{\displaystyle \frac{{\alpha }^2}{\left(r{\alpha }^2+2s+m+3\right)}}\right\} \\ =3^{s}6{\sum }_{m=0}^s{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({\displaystyle \frac{2}{3}}\right)}^m\left\{{\displaystyle \frac{1}{\left(r{\alpha }^2+2s+m+2\right)}}-{\displaystyle \frac{1}{\left(r{\alpha }^2+2s+m+3\right)}}\right\} \end{gathered}$$

2.12. PDF and CDF of Order Statistics

Let $Y_1,Y_2,\dots ,Y_n$ be a sample randomly drawn from a MBUR distribution with sample size n and the corresponding order statistic distribution is $Y_{\left(1\right)},Y_{\left(2\right)},\dots ,Y_{\left(n\right)}$ so $Y_{\left(1\right)}\le Y_{\left(2\right)}\le \dots ,{\le Y}_{\left(n\right)}$. The PDF of the jth order statistics $Y_{\left(j\right)}$ is defined as:

$$f_{Y_{\left(j\right)}}\left(y;\alpha \right)={\displaystyle \frac{n!}{\left(j-1\right)!\left(n-j\right)!}}f\left(y;\alpha \right){\left[F\left(y;\alpha \right)\right]}^{j-1}{\left[1-F\left(y;\alpha \right)\right]}^{n-j}$$

where $c={\displaystyle \frac{n!}{\left(j-1\right)!\left(n-j\right)!}}$

$$f_{Y_{\left(j\right)}}\left(y;\alpha \right)=c\left[{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\right]{\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^{j-1}{\left[1-{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}+{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^{n-j}$$

PDF of the first or the smallest order statistics is given in equation (35):

$$f_{Y_{\left(1\right)}}\left(y;\alpha \right)={\displaystyle \frac{n!}{\left(n-1\right)!}}\left[{\displaystyle \frac{6}{{\alpha }^2}}\left(1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\right]{\left[1-{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}+{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^{n-1}\dots \dots \dots \dots \dots \dots \dots \dots .$$

(35)

PDF of the largest order statistics is given in equation (36):

$$f_{Y_{\left(n\right)}}\left(y;\alpha \right)={\displaystyle \frac{n!}{\left(n-1\right)!}}\left[{\displaystyle \frac{6}{{\alpha }^2}}\left(1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\right]{\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^{n-1}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(36)

The CDF of the jth order statistics is defined as follows

$$F_{j:n}\left(y\right)=\sum _j^n{C}_{n:j}{\left[F\left(y\right)\right]}^j{\left[1-F\left(y\right)\right]}^{n-j}$$

So CDF of the jth order statistics from a MBUR random sample is defined in equation (37)

$$F_{j:n}\left(y\right)=\sum _j^n\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right){\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^j{\left[1-{3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}+{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}\right]}^{n-j}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(37)

3.1. Methods of Estimations

3.1.1. Method of Moments (MOM)

Equating the sample's first moment equation (38), which is the mean, with the population's first moment equation (39) can provide an estimate for the parameter. This estimate can then be used as an initial guess in other methods that require numerical techniques to evaluate the parameter.

$$samplemean=\overline{y}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_i\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(38)

$$populationmean=E\left(y\right)={\displaystyle \frac{6}{\left(2+{\alpha }^2\right)\left(3+{\alpha }^2\right)}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(39)

Equate the population’s mean with the sample’s mean to estimate the parameter equation (40).

$$E\left(y\right)={\displaystyle \frac{6}{\left(2+{\alpha }^2\right)\left(3+{\alpha }^2\right)}}=\overline{y}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(40)

To find the estimator for the alpha parameter, find the root of equation (41):

$$c={\displaystyle \frac{6}{\overline{y}}}=\left(2+{\alpha }^2\right)\left(3+{\alpha }^2\right)=6+5{\alpha }^2+{\alpha }^4$$

$$0={\left({\alpha }^2\right)}^2+5{\alpha }^2+6-{\displaystyle \frac{6}{\overline{y}}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(41)

$${\alpha }^2={\displaystyle \frac{-5+\sqrt{25-4\left(1\right)\left(6-\frac{6}{\overline{y}}\right)}}{2\left(1\right)}}$$

$$\widehat{\alpha }={\left\{{\displaystyle \frac{-5+\sqrt{25-24\left(\frac{\overline{y}-1}{\overline{y}}\right)}}{2}}\right\}}^{.5},25-24\left({\displaystyle \frac{\overline{y}-1}{\overline{y}}}\right)>0,\alpha >0$$

3.1.2. Maximum Likelihood Estimation (MLE)

Taking the first derivative of the log likelihood of the PDF of the MBUR distribution with respect to the alpha parameter as shown below:

The objective function to be maximized is the log-likelihood function as shown in equation (42) & (43).

$$L\left(\alpha \right)={\left({\displaystyle \frac{6}{{\alpha }^2}}\right)}^n\prod _{i=1}^n\left[1-y_i^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]\prod _{i=1}^n{y}_i^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)},0<y<1,\alpha >0\dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(42)

Maximize the log-likelihood function in equation (43) and this is equivalent to maximizing the likelihood function in equation (42) under certain regularity conditions. This is carried out by differentiating it as in equation (44) & (45) and using non-linear optimization algorithm like Newton-Raphson algorithm or quasi-Newton algorithm whichever is suitable.

$$l\left(\alpha \right)=nln\left(6\right)-2nln\left(\alpha \right)+\sum _{i=1}^{n}ln\left[1-y_i^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]+\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)\sum _{i=1}^{n}ln\left(y_i\right)\dots \dots \dots \dots \dots \dots \dots ..$$

(43)

$${\displaystyle \frac{dl\left(\alpha \right)}{d\alpha }}={\displaystyle \frac{-2n}{\alpha }}+\sum _{i=1}^n{\displaystyle \frac{-y_i^{\frac{1}{{\alpha }^2}}\left(ln\left[y_i\right]\right)\left(-2{\alpha }^{-3}\right)}{1-y_i^{\frac{1}{{\alpha }^2}}}}-4{\alpha }^{-3}\sum _{i=1}^{n}ln\left(y_i\right)\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(44)

$${\displaystyle \frac{d^{2}l\left(\alpha \right)}{d{\alpha }^2}}={\displaystyle \frac{2n}{{\alpha }^2}}-{\displaystyle \frac{6}{{\alpha }^4}}\sum _{i=1}^n{\displaystyle \frac{y_i^{\frac{1}{{\alpha }^2}}\left(ln\left[y_i\right]\right)}{1-y_i^{\frac{1}{{\alpha }^2}}}}+{\displaystyle \frac{12}{{\alpha }^4}}\sum _{i=1}^{n}ln\left(y_i\right)-{\displaystyle \frac{4}{{\alpha }^6}}\sum _{i=1}^n{\displaystyle \frac{y_i^{\frac{1}{{\alpha }^2}}{\left(ln\left[y_i\right]\right)}^2}{{\left(1-y_i^{\frac{1}{{\alpha }^2}}\right)}^2}}\dots \dots \dots \dots \dots \dots \dots$$

(45)

Alpha can be estimated numerically using Newton Raphson algorithm.

$${\alpha }_1={\alpha }_0-{\left[{\displaystyle \frac{d^{2}l\left(\alpha \right)}{d{\alpha }^2}}\right]}^{-1}{\displaystyle \frac{dl\left(\alpha \right)}{d\alpha }}$$

The expected information matrix$J_y$ is defined as the negative expected value of the second derivative of the log-likelihood, evaluated at the estimated parameter. Under certain regularity conditions and with a large sample size, the inverse of this information matrix, ${\left[{\displaystyle \frac{d^{2}l\left(\alpha \right)}{d{\alpha }^2}}\right]}^{-1},$represents the variance of the estimated parameter. Consequently, this estimator is approximately normally distributed.

$$\sqrt{n}\left(\alpha -\widehat{\alpha }\right)~N\left(0,J_y^{-1}\right)$$

This is used to construct an approximate confidence interval for the parameter in the form of $\widehat{\alpha }\pm z_{1-{\displaystyle \frac{\gamma }{2}}}se\left(\widehat{\alpha }\right)$. SE is the square root of the inverse of the expected information matrix. $z_{1-{\displaystyle \frac{\gamma }{2}}}$ is the quantile of standard normal distribution.

3.1.3. Maximum Product of Spacing (MPS)

Maximize the following objective function in equation (46):

$$MPS={\displaystyle \frac{1}{n+1}}\sum _{i=1}^{n}log\left[F\left(y_i\right)-F\left(y_{i-1}\right)\right]\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(46)

Take the derivatives of this function as in equations (47) & (48)

$${\displaystyle \frac{dMPS}{d\alpha }}={\displaystyle \frac{1}{n+1}}\sum _{i=1}^n{\displaystyle \frac{\left[F\prime \left(y_i\right)-F\prime \left(y_{i-1}\right)\right]}{\left[F\left(y_i\right)-F\left(y_{i-1}\right)\right]}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(47)

$${\displaystyle \frac{d^{2}MPS}{d{\alpha }^2}}={\displaystyle \frac{1}{n+1}}\sum _{i=1}^n{\displaystyle \frac{\left[F″\left(y_i\right)F″\left(y_{i-1}\right)\right]}{\left[F\left(y_i\right)-F\left(y_{i-1}\right)\right]}}-{\displaystyle \frac{{\left[F^{\prime \left(y_i\right)}-F^{\prime }\left(y_{i-1}\right)\right]}^2}{{\left[F\left(y_i\right)-F\left(y_{i-1}\right)\right]}^2}}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(48)

Alpha can be estimated numerically using Newton Raphson method.

3.1.4. Anderson Darling Estimator (AD)

Minimize the following objective function in equation (49):

$$AD=-n-\sum _{i=1}^n\left({\displaystyle \frac{2i-1}{n}}\right)\left\{log\left[F\left(y_i\right)\right]+log\left[1-F\left(y_{n-i+1}\right)\right]\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(49)

Take the derivatives of this function as in equations (50) & (51)

$${\displaystyle \frac{dAD}{d\alpha }}=-\sum _{i=1}^n\left({\displaystyle \frac{2i-1}{n}}\right)\left\{{\displaystyle \frac{F^{\prime }\left(y_i\right)}{F\left(y_i\right)}}+{\displaystyle \frac{-F^{\prime }\left(y_{n-i+1}\right)}{1-F\left(y_{n-i+1}\right)}}\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(50)

$${\displaystyle \frac{d^{2}AD}{d{\alpha }^2}}=-\sum _{i=1}^n\left({\displaystyle \frac{2i-1}{n}}\right)\left\{{\displaystyle \frac{F^{″}\left(y_i\right)}{F\left(y_i\right)}}-{\displaystyle \frac{{\left[F^{\prime }\left(y_i\right)\right]}^2}{{\left[F\left(y_i\right)\right]}^2}}-{\displaystyle \frac{F^{″}\left(y_{n-i+1}\right)}{1-F\left(y_{n-i+1}\right)}}-{\displaystyle \frac{{\left[F^{\prime }\left(y_{n-i+1}\right)\right]}^2}{{\left[1-F\left(y_{n-i+1}\right)\right]}^2}}\right\}\dots \dots .$$

(51)

Alpha can be estimated numerically using Newton Raphson method.

3.1.5. Percentile method (PERC)

Minimize the following objective function in equation (52):

$$Percentile=\sum _{i=1}^n{\left\{y_i-F^{-1}\left(\alpha ,ecdf={\displaystyle \frac{i}{n+1}}\right)\right\}}^2\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(52)

Take the derivatives of this function as in equations (53) & (54)

$${\displaystyle \frac{dPercentile}{d\alpha }}=2\sum _{i=1}^n\left\{y_i-F^{-1}\right\}{\left\{-F^{-1}\right\}}^{\prime }\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(53)

$${\displaystyle \frac{d^{2}Percentile}{d{\alpha }^2}}=2\sum _{i=1}^n\left\{{\left(-F^{-1}\right)}^{\prime }{\left(-F^{-1}\right)}^{\prime }+\left(y_i-F^{-1}\right){\left(-F^{-1}\right)}^{″}\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(54)

Alpha can be estimated numerically using Newton Raphson method.

3.1.6. Cramer Von Mises(CVM)

Minimize the following objective function in equation (55):

$$CVM={\displaystyle \frac{1}{12n}}+\sum _{i=1}^n{\left\{F\left(y_i\right)-{\displaystyle \frac{2i-1}{2n}}\right\}}^2\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(55)

Take the derivatives of this function as in equations (56) & (57)

$${\displaystyle \frac{dCVM}{d\alpha }}=2\sum _{i=1}^n\left\{F\left(y_i\right)-{\displaystyle \frac{2i-1}{2n}}\right\}\left\{F\prime \left(x_i\right)\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(56)

$${\displaystyle \frac{d^{2}CVM}{d{\alpha }^2}}=2\sum _{i=1}^n\left\{F^{\prime }\left(y_i\right)F^{\prime }\left(y_i\right)+\left(F\left(y_i\right)-{\displaystyle \frac{2i-1}{2n}}\right)\left(F″\left(y_i\right)\right)\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(57)

Alpha can be estimated numerically using Newton Raphson method.

3.1.7. Least Squares Method (LS)

Minimize the following objective function in equation (58):

$$LS=\sum _{i=1}^n{\left\{F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right\}}^2\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(58)

Take the derivatives of this function as in equations (59) & (60)

$${\displaystyle \frac{\partial LS}{\partial \alpha }}=2\sum _{i=1}^n\left\{F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right\}{F}^{\prime }\left(y_i\right)\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(59)

$${\displaystyle \frac{{\partial }^{2}LS}{\partial {\alpha }^2}}=2\sum _{i=1}^n\left\{F\prime \left(y_i\right)F\prime \left(y_i\right)+\left(F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right)F″\left(y_i\right)\right\}\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ..$$

(60)

Alpha can be estimated numerically using Newton Raphson method.

3.1.8. Weighted Least Squares Method (WLS)

Minimize the following objective function in equation (61):

$$WLS=\sum _{i=1}^n{\left[{\displaystyle \frac{{\left(n+1\right)}^2\left(n+2\right)}{i\left(n+1-i\right)}}\right]\left\{F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right\}}^2\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

(61)

Take the derivatives of this function as in equations (62) & (63)

$${\displaystyle \frac{dLS}{d\alpha }}=2\sum _{i=1}^n\left[{\displaystyle \frac{{\left(n+1\right)}^2\left(n+2\right)}{i\left(n+1-i\right)}}\right]\left\{F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right\}{F}^{\prime }\left(y_i\right)\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots .$$

(62)

$${\displaystyle \frac{d^{2}LS}{d{\alpha }^2}}=2\sum _{i=1}^n\left[{\displaystyle \frac{{\left(n+1\right)}^2\left(n+2\right)}{i\left(n+1-i\right)}}\right]\left\{F\prime \left(y_i\right)F\prime \left(y_i\right)+\left(F\left(y_i\right)-{\displaystyle \frac{i}{n+1}}\right)F″\left(y_i\right)\right\}\dots \dots \dots \dots .$$

(63)

Alpha can be estimated numerically using Newton Raphson method.

3.2. Simulation

A simulation study is conducted utilizing the following sample sizes $n=\left(20,80,160,260,500\right)$ , and replicate N=1000 times. Various methods of estimation are utilized and compared with one another. The parameter alpha value chosen is $\alpha =\left(2.5\right)$. For alpha value $=\left(0.5\right)$ , see the results in supplementary material section 2.

Steps:

1-

Generate random variable from the MBUR Distribution with specified alpha

2-

Choose the sample size n.

3-

Replicate the method of estimation N times.

4-

Calculate various metrics to compare methods and assess the impact of increasing sample size on estimators.

a)

$AverageAbsoluteBias(AAB)={\displaystyle \frac{1}{N}}\sum _{i=1}^N\left|\widehat{\alpha }-\alpha \right|$

b)

$MeanSquareError\left(MSE\right)={\displaystyle \frac{1}{N}}\sum _{I=1}^N{\left(\widehat{\alpha }-\alpha \right)}^2$

c)

$RootofMeanSquareError\left(MSE\right)=\sqrt{MSE}$

d)

$MeanRelativeError(MRE)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\displaystyle \frac{\left|\widehat{\alpha }-\alpha \right|}{\alpha }}$

Also the mean of estimated alpha from the 1000 replicates is evaluated with the standard error. For the chosen alpha level 2.5, the following results are obtained in the successive Table 2, Table 3, Table 4, Table 5 and Table 6, mean in Table 2, SE in Table 3, AAB in Table 4, MSE in Table 5, MRE in Table 6. For better visualization of the results each table is represented by a Heat-map graph as shown in Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17

The tables indicate that increasing the sample size leads to a decreases in the SE, AAB, MSE and MRE indices. For each method used, the indices decrease as the sample size increases. The values obtained from the MLE and MPS methods are nearly equal, especially with larger sample sizes (n=260 and n=500). Similarly, the AD and WLS methods yield comparable results. Additionally, the CVM and LS methods show approximately equal indices as the sample size increases. Overall, the methods demonstrate consistent results regarding the estimation values. Figure 16 and Figure 17 display these findings

4.1. Analysis of the first data set. See Supplementary Materials (Section 3)

4.2. Analysis of the Second Data Set. See Table 12 & 13. See Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23

Table 12. Descriptive statistics of the second data set.

min	mean	std	skewness	kurtosis	25percentile	50perc	75perc	max
0.77	0.9005	0.064	-0.9147	2.6716	0.865	0.92	0.95	0.98

The data demonstrates a left skewness and a negative excess kurtosis, indicating a platykurtic shape. This is supported by the histogram and box plot shown in Figure 18.

Table 12. Descriptive statistics of the second data set.

min	mean	std	skewness	kurtosis	25percentile	50perc	75perc	max
0.77	0.9005	0.064	-0.9147	2.6716	0.865	0.92	0.95	0.98

The data demonstrates a left skewness and a negative excess kurtosis, indicating a platykurtic shape. This is supported by the histogram and box plot shown in Figure 18.

After fitting the MBUR distribution to the data, the scaled TTT plot reveals a concave shape, indicating an increased failure rate. This pattern is evident in the shape of the hazard rate function. Figure 18 illustrates this concavity through the theoretical scaled TTT plot, which supports the increase in the hazard rate. Similarly, Figure 19 displays the empirical scaled TTT plot, also confirming this increase in the hazard rate.

Out of the five distributions analyzed, all successfully fitted the data. However, the MBUR distribution emerged as the most effective model. This was evidenced by its superior validation indices compared to the other distributions, as illustrated in in Table 13.

The KS test effectively measures the maximum distance between the empirical cumulative distribution function (eCDF) and the theoretical cumulative distribution function (CDF). Its straightforward nature and broad applicability make it a valuable tool, as it imposes no assumptions on the distribution parameters. However, it is less sensitive to deviations in the tails of the distribution, as it primarily focuses on the center. In contrast, the AD test excels in detecting deviations in the tails and is particularly suited for distributions with extreme values. Despite this advantage, the necessity for calculating critical values for newly emerging distributions can hinder its application. The CVM test takes a different approach by measuring the overall distance between the eCDF and the theoretical CDF, treating all parts of the distribution equally. This means it effectively balances sensitivity to deviations across the tail and the center, making it a compelling choice in many scenarios. Given the complexity of skewed data, it is crucial to utilize more than one test. Each test highlights specific characteristics of the data, offering a more comprehensive understanding of the fitting distribution. When combined with visual aids, such as QQ plots and PP plots, this methodology significantly enhances the analysis, driving more informed decisions. Therefore, when assessing the goodness of fit of a distribution, it is important to consider the results of the three tests mentioned above, along with the information obtained from the QQ plot and PP plot. Key aspects to observe include how closely the points align with the diagonal, the degree of deviation from the diagonal, and the percentage of observations that deviate from it.

The MBUR model fits the second dataset, as evidenced by its failure to reject the Kolmogorov-Smirnov (KS) test. The QQ plot shows almost perfect alignment with the diagonal, with only slight deviations at the lower end of the distribution. Since the KS test is less sensitive to deviations in the tail of the distribution, the author also conducted the Anderson-Darling (AD) test and the Crámer-von Mises (CVM) test using Monte Carlo simulations.

The observed value of the AD test statistic was 0.3184, while the critical values obtained from the simulations were 2.4433 (95th quantile) and 3.0146 (97.5th quantile). Since the observed AD value is less than the critical values from the simulation, the author fails to reject the null hypothesis, indicating that MBUR could be a generating process for the data. The approximate p-value for this test was 0.929, which is greater than 0.025, further confirming that MBUR fits the second dataset.Additionally, the CVM test from the observed data revealed a value of 0.0407. The CVM test conducted using Monte Carlo simulations yielded critical values of 0.4578 (95th quantile) and 0.5781 (97.5th quantile). Again, since the observed CVM value is less than the critical values, the author fails to reject the null hypothesis that the data was generated by MBUR. The approximate p-value for this test is 0.936, which is also greater than 0.025, supporting the conclusion that MBUR fits the second dataset. Overall, combining various goodness-of-fit statistics with visualizations enhances the results of the analysis. This is shown in the Table 13 and Figure 20, Figure 21, Figure 22 and Figure 23

Table 13. Estimators and validation indices for the Second data set.

	Beta		Kumaraswamy		MBUR	Topp-Leone	Unit-Lindley
theta	$\alpha =21.7353$		$\alpha =16.5447$		0.3591	71.2975	0.1334
	$\beta =2.4061$		$\beta =2.772$
Var	86.461	9.0379	15.7459	3.2005	0.000837	254.1667	0.00045
	9.0379	1.0646	3.2005	1.0347
SE	9.2984		3.9681		0.0289	15.9426	0.0212
	1.0318		1.0172
AIC	-56.5056		-56.7274		-58.079	-56.6796	-57.3746
CAIC	-55.7997		-56.0215		-57.8567	-56.4574	-57.1523
BIC	-54.5141		-54.7359		-57.0832	-55.6839	-56.3788
HQIC	-56.1168		-56.3386		-57.8846	-56.4852	-57.1802
LL	30.2528		30.3637		30.0395	29.3398	29.6873
K-S	0.0974		0.0995		0.1309	0.1327	0.1057
H0	Fail to reject		Fail to reject		Fail to reject	Fail to reject	Reject to reject
P-value	0.9416		0.9513		0.8399	0.4627	0.954
AD	0.3828		0.3527		0.3184	0.9751	0.2749
CVM	0.0566		0.0498		0.0407	0.1719	0.0261

Table 13. Estimators and validation indices for the Second data set.

	Beta		Kumaraswamy		MBUR	Topp-Leone	Unit-Lindley
theta	$\alpha =21.7353$		$\alpha =16.5447$		0.3591	71.2975	0.1334
	$\beta =2.4061$		$\beta =2.772$
Var	86.461	9.0379	15.7459	3.2005	0.000837	254.1667	0.00045
	9.0379	1.0646	3.2005	1.0347
SE	9.2984		3.9681		0.0289	15.9426	0.0212
	1.0318		1.0172
AIC	-56.5056		-56.7274		-58.079	-56.6796	-57.3746
CAIC	-55.7997		-56.0215		-57.8567	-56.4574	-57.1523
BIC	-54.5141		-54.7359		-57.0832	-55.6839	-56.3788
HQIC	-56.1168		-56.3386		-57.8846	-56.4852	-57.1802
LL	30.2528		30.3637		30.0395	29.3398	29.6873
K-S	0.0974		0.0995		0.1309	0.1327	0.1057
H0	Fail to reject		Fail to reject		Fail to reject	Fail to reject	Reject to reject
P-value	0.9416		0.9513		0.8399	0.4627	0.954
AD	0.3828		0.3527		0.3184	0.9751	0.2749
CVM	0.0566		0.0498		0.0407	0.1719	0.0261

According to the AIC, AIC corrected, BIC, and Hannan–Quinn Information Criterion (HQIC), the MBUR distribution is the best fit for the data, followed by the Unit Lindley, Topp-Leone, Kumaraswamy, and finally the Beta distribution. This conclusion is based on the MBUR having the lowest values of these indices (or the largest negative values). However, it is worth noting that the MBUR has the second lowest value for the Anderson-Darling (AD) test, Cramer-von Mises (CVM) test, and Kolmogorov-Smirnov (KS) test, coming in just after the Unit Lindley. Figure 20 illustrates that the theoretical cumulative distribution functions (CDFs) for the various distributions closely follow the empirical CDF. Meanwhile, Figure 21 presents the fitted probability density functions (PDFs). An important observation from this analysis is that the metric values for the MBUR distribution are comparable to those of the Topp-Leone and Unit Lindley distributions, indicating that the new MBUR distribution has performed well in fitting the data. Figure 22 shows the quantile-quantile (QQ) plot for the fitted MBUR distribution, which exhibits nearly perfect alignment along the diagonal, with only slight deviations at the lower tail. The log-likelihood function is maximized at an alpha level of 0.3519. Finally, Figure 23 provides the QQ plot and probability-probability (PP) plot for the other distributions, which also demonstrate near-perfect alignment along the diagonal.

4.3. Analysis of the Third Data Set. See Supplementary Materials (Section 3)

4.4. Analysis of the Fourth Data Set. See Supplementary Materials (Section 3)

4.5. Analysis of the Fifth Data Set. See Table 14 & 15. See Figure 24, Figure 25, Figure 26, Figure 27, Figure 28 and Figure 29

The data shows a right skewness and a positive excess kurtosis (leptokurtic shape), which is supported by the histogram and box plot in Figure 24. The second approach illustrated in Figure 24 for calculating and graphing the TTT plot does not exhibit the typical convexity followed by concavity that characterizes the bathtub shape seen in the hazard rate function. In contrast, Figure 25, which employs the first approach for calculation and graphing, more accurately represents this relationship.

The MBUR distribution is the best fit for the time between failures data among the five distributions evaluated, followed by Kumaraswamy, Beta, and Topp-Leone. The Unit Lindley distribution, however, did not fit the data well. The MBUR has the most significant negative values for AIC, AIC corrected, BIC, and HQIC. Despite this, it is the second distribution to have the smallest values for the AD test and the CVM test. Figure 26 illustrates that the eCDF closely follows the theoretical CDF for the fitted distributions, particularly at the tails, though there is a slight deviation at the center. Figure 27 displays the PDFs for the various competing distributions. In Figure 28, the QQ plot demonstrates good alignment with the diagonal after fitting the MBUR, with the maximum likelihood estimate achieved at an alpha level of 1.7886. Figure 29 shows a generally close alignment with the diagonal line for the other fitted distributions, especially at the tails, with slight deviations at the center, indicating that these distributions capture the characteristics of the data well. The PP plot further illustrates that the Unit Lindley distribution does not align closely with the diagonal.

Table 15. Estimators and validation indices for the Fifth data set.

	Beta		Kumaraswamy		MBUR	Topp-Leone	Unit-Lindley
theta	$\alpha =0.6307$		$\alpha =0.6766$		1.7886	0.4891	4.1495
	$\beta =3.2318$		$\beta =2.936$
Var	0.071	0.2801	0.0198	0.1033	0.018	0.0104	0.5543
	0.2801	1.647	0.1033	0.9135
SE	0.2665		0.1407		0.1342	0.0213	0.7445
	1.2834		0.9558
AIC	-36.0571		-36.6592		-37.862	-35.5653	-27.007
CAIC	-35.4571		-36.0592		-37.6712	-35.3749	-26.8165
BIC	-33.7861		-34.3882		-36.7262	-34.4298	-25.8715
HQIC	-35.4859		-36.0881		-37.5764	-35.2798	-26.7214
LL	20.0285		20.3296		19.9310	18.7827	14.5035
K-S	0.1541		0.1393		0.1584	0.1962	0.3274
H0	Fail to reject		Fail to reject		Fail to reject	Fail to Reject	Reject
P-value	0.5918		0.7123		0.5575	0.2982	0.0107
AD	0.6886		0.5755		0.6703	1.1022	4.7907
CVM	0.1264		0.0989		0.1253	0.2149	0.8115

Table 15. Estimators and validation indices for the Fifth data set.

	Beta		Kumaraswamy		MBUR	Topp-Leone	Unit-Lindley
theta	$\alpha =0.6307$		$\alpha =0.6766$		1.7886	0.4891	4.1495
	$\beta =3.2318$		$\beta =2.936$
Var	0.071	0.2801	0.0198	0.1033	0.018	0.0104	0.5543
	0.2801	1.647	0.1033	0.9135
SE	0.2665		0.1407		0.1342	0.0213	0.7445
	1.2834		0.9558
AIC	-36.0571		-36.6592		-37.862	-35.5653	-27.007
CAIC	-35.4571		-36.0592		-37.6712	-35.3749	-26.8165
BIC	-33.7861		-34.3882		-36.7262	-34.4298	-25.8715
HQIC	-35.4859		-36.0881		-37.5764	-35.2798	-26.7214
LL	20.0285		20.3296		19.9310	18.7827	14.5035
K-S	0.1541		0.1393		0.1584	0.1962	0.3274
H0	Fail to reject		Fail to reject		Fail to reject	Fail to Reject	Reject
P-value	0.5918		0.7123		0.5575	0.2982	0.0107
AD	0.6886		0.5755		0.6703	1.1022	4.7907
CVM	0.1264		0.0989		0.1253	0.2149	0.8115

The MBUR model unequivocally fits the fifth dataset, as demonstrated by the results of the Kolmogorov-Smirnov (KS) test, which successfully failed to reject the null hypothesis that the data adheres to the MBUR distribution. This conclusion is visually substantiated by the QQ plot, which shows a strong alignment along the diagonal, indicating a robust correspondence between the theoretical distribution and the empirical data, though minor deviations can be observed at the lower end. To address the potential limitations of the KS test—particularly its insensitivity to deviations in the distribution tails—the author conducted additional analyses using the AD test and the CVM test, both of which utilized Monte Carlo simulations for enhanced accuracy. The AD test produced a statistic of 0.6703. The critical values from the simulations for this test were clear: 2.6428 for the 95th quantile and 3.3935 for the 97.5th quantile. The value for 2.5th quantile is 0.2309. The observed AD statistic is significantly less than these critical values, compellingly leading the author to fail to reject the null hypothesis. This strongly indicates that the MBUR model can indeed act as a generating process for the observed data. The p-value corresponding to this test was a definitive 0.594, exceeding the conventional significance threshold of 0.025, thus firmly establishing MBUR as an appropriate fit for the fifth dataset. Likewise, the CVM test reinforced this conclusion with an observed statistic of 0.1253. Critical values derived from Monte Carlo simulations were found to be 0.4858 (95th quantile) and 0.6099 (97.5th quantile). The value for the 2.5th quantile is 0.03. As with the AD test, the observed CVM statistic fell below these critical thresholds, leading to a resolute failure to reject the null hypothesis that the data originated from the MBUR model. The approximate p-value for this test, calculated at 0.485, also exceeds the 0.025 significance level, decisively affirming that MBUR is an excellent fit for the dataset in question.

Integrating various goodness-of-fit statistics with effective visualizations significantly enhances the analysis results, leading to clearer insights and more informed decisions. This is shown in the Table15 and Figure 26, Figure 27, Figure 28 and Figure 29

When using AIC and BIC to compare distributions that fit specific data, both metrics aim to balance maximizing model fit, reflected in the highest negative values of log-likelihood, with minimizing model complexity, which is represented by the number of parameters in the model. This balance helps avoid overfitting, particularly in cases where a model may be too complex and have too many parameters. Such complex models can capture not only the true underlying structure of the data but also random noise, leading to poor generalization when new data is introduced.

The log-likelihood (LL) measures how well a model fits the data; higher LL values indicate a better fit. However, simply adding more parameters (k) tends to increase LL, even if those additional parameters are not meaningful. Therefore, AIC and BIC serve as trade-offs between model fit and complexity, addressing the challenge of balancing overfitting (too complex a model) and underfitting (too simple a model). To mitigate overfitting, AIC and BIC introduce penalties for complexity that are proportional to the number of parameters (k). The AIC penalty is a linear penalty of 2k, whereas the BIC penalty is k * ln(n), which increases with sample size.

AIC and BIC are used to select the best-fitting distribution among candidates. They depend on LL, meaning that the model which better captures the structure of the data will display more negative values for AIC and BIC. In cases where the data exhibits complex features such as skewness and heavy tails, a model with more parameters may yield more negative AIC and BIC values. Conversely, if the data is simpler (e.g., symmetric with a small sample size), more straightforward models are often preferred.

The more negative the values of AIC and BIC, the better the model. By themselves, these values are meaningless, but they are useful for comparing models. A difference greater than 10 between two models suggests that the model with the more negative value is significantly better. AIC is typically used when the goal is prediction and the dataset is small, while BIC is preferred when identifying the true model is critical and the dataset is large, due to differing penalty structures.

Regarding the datasets discussed in this paper, they exhibit skewness and kurtosis, indicating their complexity. The sizes of these datasets range between 20 and 36, making them small to moderate in size. The new MBUR distribution can effectively fit all the data using just one parameter, resulting in a relatively small penalty from AIC and BIC. This represents an advantage of the MBUR distribution over other distributions, such as the beta and Kumaraswamy distributions, which require multiple parameters.