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)}\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }0<y<1\mathrm{ },\mathrm{ }\mathrm{ }\alpha >0$$

(5)

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

(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)\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }0<y<1\mathrm{ }\mathrm{ },\alpha >0$$

(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)}}\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }0<y<1\mathrm{ }\mathrm{ },\mathrm{ }\alpha >0$$

(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}}}}\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }0<y<1\mathrm{ }\mathrm{ },\mathrm{ }\alpha >0\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 3^rd 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}\mathrm{ }sin\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]\right)+.5\right\}}^{{\alpha }^2}$$

(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. r^th Raw Moments for a Random Variable y Distributed as MBUR is Shown in Equation (11)

$$\begin{gathered} E\left(y^r\right)={\int }_0^1{y}^r\mathrm{ }{\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\mathrm{ }dy\mathrm{ }={\displaystyle \frac{6}{\left(2+r{\alpha }^2\right)\left(3+r{\alpha }^2\right)}}\dots \\ E\left(y\right)={\displaystyle \frac{6}{\left(2+{\alpha }^2\right)\left(3+{\alpha }^2\right)}}\mathrm{ }\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }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)}}\mathrm{ },\mathrm{ }\mathrm{ }E\left(y^4\right)={\displaystyle \frac{6}{\left(2+4{\alpha }^2\right)\left(3+4{\alpha }^2\right)}} \\ \mathrm{v}\mathrm{a}\mathrm{r}\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}} \end{gathered}$$

(11)

Coefficient of Skewness: is shown in equation (12)

$$\begin{gathered} 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{ }\mathrm{o}\mathrm{f}\mathrm{ }\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 .\mathrm{ } \end{gathered}$$

(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\mathrm{ }\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\mathrm{ }\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\mathrm{ }\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}} \\ \mathrm{coefficient} \mathrm{of} \mathrm{Kurtosis}={\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 \end{gathered}$$

(13)

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

$$CV={\displaystyle \frac{S}{\mu }}$$

(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. r^th 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)}}$$

(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)

$$\begin{gathered} 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 \end{gathered}$$

(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)

$$\begin{gathered} L\left(p\right)={\displaystyle \frac{{\int }_0^{t}y\mathrm{ }f\left(y\right)dy}{{\int }_0^{1}yf\left(y\right)dy}}=\mathrm{ }{y}^{{\displaystyle \frac{2}{{\alpha }^2}}+1}\left\{\left(3+{\alpha }^2\right)-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\left(2+{\alpha }^2\right)\mathrm{ }\right\}\dots \dots \\ B\left(p\right)={\displaystyle \frac{L\left(p\right)}{F\left(y\right)}}={\displaystyle \frac{\mathrm{ }{y}^{\frac{2}{{\alpha }^2}+1}\left\{\left(3+{\alpha }^2\right)-y^{\frac{1}{{\alpha }^2}}\left(2+{\alpha }^2\right)\mathrm{ }\right\}}{{3y}^{\frac{2}{{\alpha }^2}}-{2y}^{\frac{3}{{\alpha }^2}}}} \end{gathered}$$

(17)

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

(18)

$$\mathrm{G}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{ }\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=1-2\mathrm{ }{\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)\mathrm{ }\mathrm{ }$$

(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\}\mathrm{ }$$

(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 .$$

(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$$

(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)

$$\begin{gathered} \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}}}} \end{gathered}$$

(26)

$$\begin{gathered} \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 \end{gathered}$$

(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}Y hence : X{\le }_{st}Y$$

Theorem: let $X~MBUR\left({\alpha }_1\right)$ and $Y~\mathrm{ }MBUR\left({\alpha }_2\right)$  if  ${\alpha }_1<\mathrm{ }{\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)}}}\mathrm{ }$$

(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{d LRO}{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\}$$

(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).

$$\begin{gathered} \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 \end{gathered}$$

(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):

$$\begin{gathered} 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 \\ ={\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}$$

(32)

$$={\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$$

(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}^{\mathrm{\infty }}{\left(-1\right)}^m\mathrm{ }\left({\displaystyle \genfrac{}{}{0pt}{}{s}{m}}\right){\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^2}}}\right)}^m$$

(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 j^th 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}$$

(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}$$

(36)

The CDF of the j^th 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 j^th 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}$$

(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.

$$sample\mathrm{ }mean=\overline{y}=\mathrm{ }{\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_i$$

(38)

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

(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}$$

(40)

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

$$\begin{gathered} 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-\mathrm{ }{\displaystyle \frac{6}{\overline{y}}}\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\dots \dots {\alpha }^2={\displaystyle \frac{-5+\sqrt{25-4\left(1\right)\left(6-\mathrm{ }\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},\mathrm{ }\mathrm{ }25-24\left({\displaystyle \frac{\overline{y}-1}{\overline{y}}}\right)>0\mathrm{ },\mathrm{ }\alpha >0\mathrm{ } \end{gathered}$$

(41)

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).

$$\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]\mathrm{ }\prod _{i=1}^n{y}_i^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}\mathrm{ }\mathrm{ }\mathrm{ },\mathrm{ }\mathrm{ }\mathrm{ }0<y<1\mathrm{ },\mathrm{ }\alpha >0\dots L$$

(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)$$

(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)\mathrm{ }$$

(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}}$$

(45)

Alpha can be estimated numerically using Newton Raphson algorithm.

$${\alpha }_1={\alpha }_{0\mathrm{ }}-{\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}\mathrm{ }\left(\alpha -\widehat{\alpha }\right)~N\left(0,\mathrm{ }{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}}\mathrm{ }\sum _{i=1}^{n}log\left[F\left(y_i\right)-F\left(y_{i-1}\right)\right]\mathrm{ }$$

(46)

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

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

(47)

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

(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(x_i\right)\right]+log\left[1-F\left(x_{n-i+1}\right)\right]\right\}$$

(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^{\text{'}\left(x_i\right)}}{F\left(x_i\right)}}+{\displaystyle \frac{-F^{\text{'}\left(x_{n-i+1}\right)}}{1-F\left(x_{n-i+1}\right)}}\right\}$$

(50)

$${\displaystyle \frac{d^{2}AD}{d{\alpha }^2}}=-\sum _{i=1}^n\left({\displaystyle \frac{2i-1}{n}}\right)\left\{{\displaystyle \frac{F^{\text{'}\text{'}}\left(x_i\right)}{F\left(x_i\right)}}-{\displaystyle \frac{ {\left[F^{\text{'}\left(x_i\right)}\right]}^2}{{\left[F\left(x_i\right)\right]}^2}}-{\displaystyle \frac{F^{\text{'}\text{'}}\left(x_{n-i+1}\right)}{1-F\left(x_{n-i+1}\right)}}-{\displaystyle \frac{{\left[F^{\text{'}}\left(x_{n-i+1}\right)\right]}^2}{{\left[1-F\left(x_{n-i+1}\right)\right]}^2}}\right\}$$

(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$$

(52)

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

$${\displaystyle \frac{d Percentile}{d\alpha }}=2\sum _{i=1}^n\left\{y_i-F^{-1}\right\}{\left\{-F^{-1}\right\}}^{\text{'}}$$

(53)

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

(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(x_i\right)-{\displaystyle \frac{2i-1}{2n}}\right\}}^2$$

(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(x_i\right)-{\displaystyle \frac{2i-1}{2n}}\right\}\left\{F\text{'}\left(x_i\right)\right\}$$

(56)

$${\displaystyle \frac{d^{2}CVM}{d{\alpha }^2}}=2\sum _{i=1}^n\left\{F^{\text{'}}\left(x_i\right)F^{\text{'}}\left(x_i\right)+\left(F\left(x_i\right)-{\displaystyle \frac{2i-1}{2n}}\right)\left(F\text{'}\text{'}\left(x_i\right)\right)\right\}$$

(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(x_i\right)-{\displaystyle \frac{i}{n+1}}\right\}}^2$$

(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(x_i\right)-{\displaystyle \frac{i}{n+1}}\right\}{F}^{\text{'}\left(x_i\right)}$$

(59)

$${\displaystyle \frac{{\partial }^{2}LS}{\partial {\alpha }^2}}=2\sum _{i=1}^n\left\{F\text{'}\left(x_i\right)F\text{'}\left(x_i\right)+\left(F\left(x_i\right)-{\displaystyle \frac{i}{n+1}}\right)F\text{'}\text{'}\left(x_i\right)\right\}$$

(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(x_i\right)-{\displaystyle \frac{i}{n+1}}\right\}}^2$$

(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(x_i\right)-{\displaystyle \frac{i}{n+1}}\right\}{F}^{\text{'} }\left(x_i\right)$$

(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\text{'}\left(x_i\right)F\text{'}\left(x_i\right)+\left(F\left(x_i\right)-{\displaystyle \frac{i}{n+1}}\right)F\text{'}\text{'}\left(x_i\right)\right\}$$

(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)

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

b)

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

c)

$Root of Mean Square Error\left(MSE\right)=\sqrt{MSE}$

d)

$Mean Relative Error\left(MRE\right)=\mathrm{ }{\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

Table 2 shows that as sample size increases, the estimated parameter approaches the true value. The percentile method is the least efficient to approach the true regardless the sample size. MOM and MPS have nearly comparable results at different sample sizes. AD and CVM have nearly similar results. And LS and WLS have similar results at large sample sizes. Figure 13 visually illustrates these results.

Table 3 shows that as sample size increases the standard error (SE) decreases. Percentile method has the highest SE followed by the MOM at all different sample sizes. The MLE & MPS have the lowest SE at n=500. MLE and MPS have nearly equal results at different sample sizes. This is also true as regards the pair of AD and WLS methods and the pair of CVM and LS methods. Figure 14 visually demonstrates these results.

Table 4 shows that as the sample size increases, the average absolute bias (AAB) decreases. The percentile method has the highest value of AAB followed by MOM at all different sample sizes. MLE and MPS yield near identical results at different sample sizes. AD and WLS have approximately equal results. CVM and LS have nearly similar results. Figure 15 visually depicts these findings.

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-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, Figure 23, Figure 24 and Figure 25

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. Plotting the empirical survival function reveals a rapid decay, suggesting a light tail, as illustrated in Figure 19. This observation is further reinforced by the log-log plot in Figure 20. When the author plots the logarithm of the observations against the logarithm of the empirical survival function, the author sees a concave curve, which indicates faster decay and supports the concept of a light tail. Additionally, a quantile analysis was conducted by comparing the empirical 1st, 5th, and 10th quantiles with the corresponding theoretical quantiles from the standard uniform distribution. The findings show a light tail: specifically, (0.7700 < 0.7721), (0.7750 < 0.7805), and (0.7900 < 0.7910). Fitting the MBUR model to the data produced an estimated alpha value of 0.3591, which is less than 1. This is consistent with the empirical survival function depicted in Figure 19, which bears similarity to the survival function in Figure 6, where alpha is also less than 1.

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 21 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 P-values for the estimators of alpha and beta parameters of the Beta distribution and Kumaraswamy distributions are significant $p<0.001$.P-values for the estimators of alpha of the MBUR distribution is significant $p<0.001$ , for the estimators of theta of the Topp-Leone distribution is significant $p<0.001$, and for the estimators of theta of the Unit Lindley distribution is significant $p<0.001$.

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 Table13 and Figure 22, Figure 23, Figure 24 and Figure 25

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 22 illustrates that the theoretical cumulative distribution functions (CDFs) for the various distributions closely follow the empirical CDF. Meanwhile, Figure 23 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 24 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 25 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 and Table 15. See Figure 26, Figure 27, Figure 28, Figure 29, Figure 30, Figure 31, Figure 32 and Figure 33

The data shows a right skewness and a positive excess kurtosis (leptokurtic shape), which is supported by the histogram and box plot in Figure 26. When plotting the empirical survival function, the author observes a slower decay, indicating heavier tails, as illustrated in Figure 27. This observation is further supported by the log-log plot in Figure 28. When the author plots the logarithm of the observations against the logarithm of the empirical survival function, the resulting straight line indicates a slower decay, which is characteristic of heavier tails. In quantile analysis, the author compared the empirical 99th quantile (0.6560) with the 99th theoretical quantile of a standard uniform distribution (0.6494). This comparison reveals that the empirical value is greater than the theoretical one (0.6560 > 0.6494), suggesting a heavier tail. Additionally, when fitting the MBUR model to the data, the author obtained an estimated alpha value of 1.7886, which is greater than 1. This is consistent with the empirical survival function shown in Figure 27, which is similar to the survival function depicted in Figure 6, where alpha is also greater than 1.

The second approach illustrated in Figure 26 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 29, 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 30 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 31 displays the PDFs for the various competing distributions. In Figure 32, 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 33 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. The P-values for the estimators of alpha and beta parameters of the Beta distribution and Kumaraswamy distributions are significant $p<0.001$. P-value for the estimators of alpha of the MBUR distribution is significant $p<0.001$. P-value for the estimators of theta of the Topp-Leone distribution is significant $p<0.001$.

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 30, Figure 31, Figure 32 and Figure 33

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.