Median Based Unit Weibull Distribution (MBUW): Do the Higher Order Probability Weighted Moments (PWM) Add More Information over the Lower Order PWM in Parameter Estimation

1. Introduction

Iman Attia was the first to introduce Median Based Unit Weibull distribution (MBUW) [1]. Given a random variable y that is distributed as Median Based Unit Weibull distribution (MBUW), the PDF, CDF, and quantile functions are as follows:

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

(1)

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

(2)

$$y=F^{-1}\left(y\right)={\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 }^{\beta }}$$

(3)

Various methods are employed to estimate the parameters of a distribution. Maximum Likelihood Estimation (MLE) is popular because it leads to efficient, asymptotically minimum variance estimators, although these estimators may not necessarily be unbiased. The method of moments is straightforward to apply and can serve as starting values for the numerical procedures associated with MLE.

[2] advocated for Probability Weighted Moments (PWMs) and initiated their use in hydrology for small sample sizes, as MLE does not always perform well with limited data. PWMs are a leading alternative to method of moments (MM) and MLE for fitting statistical distributions to data, especially those expressed in inverse form. For instance, if y is a random variable and F represents the cumulative distribution function (CDF) for y, the value of y can be expressed as a function of F; $y=y\left(F\right).$ These distributions include the Gumbel, Weibull, and logistic distributions, as well as lesser-known types such as Tukey’s symmetric lambda, Thomas’ Wakeby, and Mielke’s kappa.

[3] advocated for PWMs, demonstrating their superior performance over other estimators through Monte Carlo simulations.For most of these distributions, relatively simple expressions for the parameters can be derived, as noted by [4], including several for which parameter estimates are not easily obtained using MLE or conventional moments. These distributions relate to various fields, including hydrology, resource management, and the forecasting of weather parameters such as temperature, precipitation, wind velocity, floods, droughts, and rainfall.

Many researchers have utilized PWMs, including [5] for estimating parameters of the type II extreme value distribution, and [6] who applied the generalized form of PWMs to estimate the two-parameter Weibull distribution. [7] introduced a new class of generalized PWMs for estimating parameters of the Pareto distribution, while [8] implemented partial PWMs for censored data from the generalized extreme value distribution (GEV). Many authors like [9,10,11,12,13,14,15,16] have applied PWMs to various distributions, mainly focusing on the GEV.

PWMs are more robust than conventional moments in the presence of outliers in the data and yield more efficient estimators from small samples regarding the underlying probability distribution. They are computed as expectations of certain functions of a random variable. A significant advantage over conventional moments is that PWMs are, by definition, linear functions of the ordered data, resulting in reduced sensitivity to sampling variability. Various methods are employed to estimate the parameters of a distribution.

This paper is organized as follows: Section 2 introduces the methodology and the definition of Probability Weighted Moments (PWMs). It explicates the classic PWMs method for parameter estimation. In this section the author derives the PWMs for the Median Biased Unit Weibull (MBUW) and clarifies how to utilize these moments for parameter estimation. The author also illustrates the derivation of the asymptotic variance for different estimators. Section 3 demonstrates Monte Carlo simulation study to validate the premise of the paper. Section 4 explains the results of the simulation study. Section 5 explores the application of this method in real data analysis. Section 6 comprehends the conclusions. Finally, Section 7 suggests the future work.

2. Materials and Methods:

2.1. Definitions of Probability Weighted Moments (PWMs)

The population PWM, as shown in equations (4-5), is defined as $M_{p,r,s}$, where $p,r,and s$ are real numbers.

$$M_{p,r,s}=E\left(y^p {\left[F\left(y\right)\right]}^{r } {\left[1-F\left(y\right)\right]}^s\right)=\underset{-\infty }{\overset{\infty }{\int }}{y}^p {\left[F\left(y\right)\right]}^{r } {\left[1-F\left(y\right)\right]}^s f\left(y\right)dy$$

(4)

$$M_{p,r,s}=E\left(y^p {\left[F\left(y\right)\right]}^{r } {\left[1-F\left(y\right)\right]}^s\right)=\underset{0}{\overset{1}{\int }}{\left[y\left(F\right)\right]}^p {\left[F\right]}^{r } {\left[1-F\right]}^s dF$$

(5)

where $y\left(F\right)$ is a solution to $F\left(y\right)=F$ , in other words, the distribution has a closed form quantile function, inverse CDF.

In practice, the value of $p$ is chosen to be 1, and so $M_{1,r,s}$ are used to estimate the parameter. Here, $M_{\mathrm{1,0},0}$ represents the mean. If this mean exists, then $M_{1,r,s}$ exist for any real positive values of r and s. These values are often restricted to small positive integers. Choosing p=1 has the dual advantage of not unduly overweighting sample values while also leading to a class of linear L-moments that exhibit asymptotic normality, [17,18]. While only small positive integers are needed to estimate the parameters of distributions, using real numbers, regardless of their size, can offer significant advantages. This concept was explored by [19], who extended PWMs into generalized PWMs to include those with $p\ne 1$. PWMs is a generalization of the classic method of moments when $r=s=0$ , $M_{\mathrm{1,0},0}$ , are the non-central moments of order p. [20] modified the method to accommodate the models without an analytic CDF and quantile function. [2] and [18] strongly advised employing $M_{1,r,s}$ due to the significantly clearer relationships that exist between parameters and moments. This approach simplifies analysis and enhances our understanding of these critical relationships. The empirical estimate of $M_{1,r,s}$ is typically less sensitive to outliers and exhibits good properties, especially when the sample size is small. For convenience, several authors chose to use p=1 and non-negative integer values for r and s. This approach is referred to as the classic PWM method. When p=1, r and s are non-negative, the following equations (6-7) are defined:

$${\alpha }_r=M_{\mathrm{1,0},s}=E\left(y{\left(1-F\left(y\right)\right)}^S\right) , s=\mathrm{0,1},\dots . \& r=0$$

(6)

$${\beta }_s=M_{1,r,0}=E\left(y{\left(F\left(y\right)\right)}^r\right) , r=\mathrm{0,1},\dots . \& s=0$$

(7)

Both ${\alpha }_r$ and ${\beta }_s$ are related by the following equation (8):

$${\alpha }_r=\sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right){\beta }_i and {\beta }_s=\sum _{i=0}^r{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right){\alpha }_i$$

(8)

For non-negative integers values of r and s that are as small as possible, both ${\alpha }_r$ and ${\beta }_s$ are equivalent. [21] defined the sample unbiased estimators for PWMs ${\alpha }_r$ and ${\beta }_s$ as the followings in equation (9):

$${\widehat{\alpha }}_r={\displaystyle \frac{1}{n}}\sum _{i=1}^{n-r}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n-i}{r}\right)}{\left(\genfrac{}{}{0pt}{}{n-1}{r}\right)}}{y}_{i:n } and {\widehat{\beta }}_s={\displaystyle \frac{1}{n}}\sum _{i=1+s}^n{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{i-1}{s}\right)}{\left(\genfrac{}{}{0pt}{}{n-1}{s}\right)}}{y}_{i:n }$$

(9)

The biased sample estimator for the PWM is defined as

${\widehat{M}}_{1,r,s}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^r{\left(1-p_{i:n}\right)}^s$, so it can be defined in equations (10-11) for:

$${\widehat{\alpha }}_r=M_{\mathrm{1,0},s}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^0{\left(1-p_{i:n}\right)}^s$$

(10)

$${\widehat{\beta }}_s=M_{1,r,0}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^r{\left(1-p_{i:n}\right)}^0$$

(11)

where $p_{i:n}={\displaystyle \frac{i-b}{n }} , 0\le b\le 1 or p_{i:n}={\displaystyle \frac{i-b}{n+1-2b }} ,-0.5\le b\le 0.5$

[22] empirically concluded that moderated biased estimates of the PWMs could produce more accurate estimates of upper quantiles. In this paper, $M_{\mathrm{1,0},1}$ and $M_{\mathrm{1,1},0}$ are used as a system of equations to estimate the parameters of the Median Based Unit Weibull (MBUW) distribution. Higher moments $M_{\mathrm{1,0},2}$ and $M_{\mathrm{1,2},0}$  are also used to estimate the parameters and compare the results with lower moments $M_{\mathrm{1,0},1}$  and $M_{\mathrm{1,1},0}$  . Parameter estimation using PWMs is carried out by equating the analytic expression of the population PWMs by the corresponding sample estimates of PWMs and solving the resulting systems of equations in terms of the parameters.

2.2. Calculating PWMs for MBUW

2.2.1. Calculating $M_{\mathrm{1,0},1}$: Equations (12-20)

 Substitute in equation 4 to get: $M_{p,r,s}=M_{\mathrm{1,0},1}=E\left(y {\left[1-F\left(y\right)\right]}^s\right)=\underset{0}{\overset{1}{\int }}y {\left[1-F\left(y\right)\right]}^s f\left(y\right)dy$

$$M_{\mathrm{1,0},1}=\underset{0}{\overset{1}{\int }}y {\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^s {\displaystyle \frac{6}{{\alpha }^{\beta }}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}-1\right)}dy$$

(12)

$$=\underset{0}{\overset{1}{\int }} {\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^s {\displaystyle \frac{6}{{\alpha }^{\beta }}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)}dy$$

$$=\underset{0}{\overset{1}{\int }} {\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)} {\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^{s}dy-\underset{0}{\overset{1}{\int }}{\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{3}{{\alpha }^{\beta }}}\right)}{\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^{s}dy$$

(13)

$$=\underset{0}{\overset{1}{\int }} A\left(y\right)dy$$

$$-\underset{0}{\overset{1}{\int }}B\left(y\right)dy$$

(14)

Using binomial expansion in equation (13):

$$\begin{gathered} {\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^s=\sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right){\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)}^i \\ \sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right) {\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}\right)}^i{\left(1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right)}^i= \\ \sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right) {\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}\right)}^i\sum _{m=0}^i{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{i}{m}}\right) {\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right)}^m= \\ \sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right) 3^i y^{{\displaystyle \frac{2i}{{\alpha }^{\beta }}}}\sum _{m=0}^i{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{i}{m}}\right){\left({\displaystyle \frac{2}{3}}\right)}^m y^{{\displaystyle \frac{m}{{\alpha }^{\beta }}}} \end{gathered}$$

(15)

Now exchange integration with summation in equations (13) & (14):

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

(16)

where the integral is:

$$\begin{gathered} \underset{0}{\overset{1}{\int }} y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)}{y}^{{\displaystyle \frac{2i}{{\alpha }^{\beta }}}} y^{{\displaystyle \frac{m}{{\alpha }^{\beta }}}} dy={\displaystyle \frac{{\alpha }^{\beta }}{2+2i+m+{\alpha }^{\beta }}} \\ so: \underset{0}{\overset{1}{\int }} A\left(y\right)dy=\underset{0}{\overset{1}{\int }} {\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)} {\left[1-\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right)\right]}^{s}dy= \\ {\displaystyle \frac{6}{{\alpha }^{\beta }}}\sum _{i=0}^s{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{s}{i}}\right) 3^i \sum _{m=0}^i{\left(-1\right)}^m\left({\displaystyle \genfrac{}{}{0pt}{}{i}{m}}\right){\left({\displaystyle \frac{2}{3}}\right)}^m \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2i+m+{\alpha }^{\beta }}}\right) \end{gathered}$$

(17)

The same steps follow for:

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

(18)

Substitute s=1 in equations (17) and (18):

$$\begin{gathered} \underset{0}{\overset{1}{\int }} A\left(y\right)dy= \\ {\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left(3\right)}^0 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{0}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(0\right)+0+{\alpha }^{\beta }}}\right)+{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{0}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(0\right)+1+{\alpha }^{\beta }}}\right)\right]\right. \\ \left.{+\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left(3\right)}^1 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+0+{\alpha }^{\beta }}}\right)+{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+1+{\alpha }^{\beta }}}\right)\right]\right\} \\ ={\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{{\alpha }^{\beta }}{2+{\alpha }^{\beta }}}-{\displaystyle \frac{3{\alpha }^{\beta }}{4+{\alpha }^{\beta }}}+{\displaystyle \frac{2{\alpha }^{\beta }}{5+{\alpha }^{\beta }}}\right\}=6\left\{{\displaystyle \frac{1}{2+{\alpha }^{\beta }}}-{\displaystyle \frac{3}{4+{\alpha }^{\beta }}}+{\displaystyle \frac{2}{5+{\alpha }^{\beta }}}\right\} \end{gathered}$$

(19)

For

$$\begin{gathered} \underset{0}{\overset{1}{\int }} B\left(y\right)dy= \\ {\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{ \left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left(3\right)}^0\left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{0}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(0\right)+0+{\alpha }^{\beta }}}\right)+{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{0}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(0\right)+1+{\alpha }^{\beta }}}\right)\right]\right. \\ \left.{+\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left(3\right)}^1 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+0+{\alpha }^{\beta }}}\right)+{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+1+{\alpha }^{\beta }}}\right)\right]\right\} \\ ={\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{{\alpha }^{\beta }}{3+{\alpha }^{\beta }}}-{\displaystyle \frac{{3\alpha }^{\beta }}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{{2\alpha }^{\beta }}{6+{\alpha }^{\beta }}}\right\}=6\left\{{\displaystyle \frac{1}{3+{\alpha }^{\beta }}}-{\displaystyle \frac{3}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{2}{6+{\alpha }^{\beta }}}\right\} \end{gathered}$$

(20)

Substitute equations (19) and (20) into equation (14)

$$M_{\mathrm{1,0},1}={\displaystyle \frac{360+108{\alpha }^{\beta }}{1044{\alpha }^{\beta }+580{\alpha }^{2\beta }+155{\alpha }^{3\beta }+20{\alpha }^{4\beta }+{\alpha }^{5\beta }+720}}$$

2.2.2. Calculating $M_{\mathrm{1,1},0}$ Equations (21-26)

$$M_{\mathrm{1,1},0}=E\left(y^p {\left[F\left(y\right)\right]}^{r } \right)=\underset{0}{\overset{1}{\int }}{y}^p {\left[F\left(y\right)\right]}^{r } f\left(y\right)dy$$

$$\underset{0}{\overset{1}{\int }}y {\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right]}^r {\displaystyle \frac{6}{{\alpha }^{\beta }}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}-1\right)}dy$$

(21)

$$=\underset{0}{\overset{1}{\int }} {\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)} {\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right]}^{r}dy-\underset{0}{\overset{1}{\int }}{\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{3}{{\alpha }^{\beta }}}\right)}{\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right]}^{r}dy=\underset{0}{\overset{1}{\int }} C\left(y\right)dy-\underset{0}{\overset{1}{\int }}D\left(y\right)dy$$

(22)

$$\underset{0}{\overset{1}{\int }} C\left(y\right)dy=\underset{0}{\overset{1}{\int }} {\displaystyle \frac{6}{{\alpha }^{\beta }}}{y}^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)} {\left[{3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^{\beta }}}}\right]}^{r}dy={\displaystyle \frac{6}{{\alpha }^{\beta }}} \underset{0}{\overset{1}{\int }} y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)}{\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}\right)}^r {\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right]}^{r}dy$$

Using binomial expansion: ${\left[1-{{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right]}^r=\sum _{i=0}^r{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right)}^i$

$$So \underset{0}{\overset{1}{\int }} C\left(y\right)dy={\displaystyle \frac{6}{{\alpha }^{\beta }}} \underset{0}{\overset{1}{\int }} y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)}{\left({3y}^{{\displaystyle \frac{2}{{\alpha }^{\beta }}}}\right)}^r \sum _{i=0}^r{\left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({{\displaystyle \frac{2}{3}}y}^{{\displaystyle \frac{1}{{\alpha }^{\beta }}}}\right)}^{i}dy$$

(23)

Exchange the integral and the sum ${\displaystyle \frac{6}{{\alpha }^{\beta }}} \sum _{i=0}^r{3^r \left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^i\underset{0}{\overset{1}{\int }} y^{\left({\displaystyle \frac{2}{{\alpha }^{\beta }}}\right)}{y}^{{\displaystyle \frac{2r}{{\alpha }^{\beta }}}} y^{{\displaystyle \frac{i}{{\alpha }^{\beta }}}}dy=$

$$\underset{0}{\overset{1}{\int }} C\left(y\right)dy={\displaystyle \frac{6}{{\alpha }^{\beta }}} \sum _{i=0}^r{3^r \left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^i\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2r+i+{\alpha }^{\beta }}}\right)$$

(24)

The same steps are followed and give:

$$\underset{0}{\overset{1}{\int }} D\left(y\right)dy={\displaystyle \frac{6}{{\alpha }^{\beta }}} \sum _{i=0}^r{3^r \left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^i\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2r+i+{\alpha }^{\beta }}}\right)$$

(25)

Substitute equation (24) and equation (25) into equation (22) and let r=1

$$\begin{gathered} M_{\mathrm{1,1},0}={\displaystyle \frac{6}{{\alpha }^{\beta }}} \sum _{i=0}^r{3^r \left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^i\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2r+i+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2r+i+{\alpha }^{\beta }}}\right)\right] \\ ={\displaystyle \frac{6}{{\alpha }^{\beta }}} \{{3^1 \left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^0\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+\left(0\right)+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+\left(0\right)+{\alpha }^{\beta }}}\right)\right] \\ {{+3}^1 \left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^1\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+\left(1\right)+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+\left(1\right)+{\alpha }^{\beta }}}\right)\right]\} \\ ={\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{3{\alpha }^{\beta }}{4+{\alpha }^{\beta }}}-{\displaystyle \frac{{2\alpha }^{\beta }}{5+{\alpha }^{\beta }}}-{\displaystyle \frac{{3\alpha }^{\beta }}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{{2\alpha }^{\beta }}{6+{\alpha }^{\beta }}}\right\}={\displaystyle \frac{18}{4+{\alpha }^{\beta }}}-{\displaystyle \frac{30}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{12}{6+{\alpha }^{\beta }}} \\ M_{\mathrm{1,1},0}={\displaystyle \frac{6\left(10+{\alpha }^{\beta }\right)}{\left(4+{\alpha }^{\beta }\right)\left(5+{\alpha }^{\beta }\right)\left(6+{\alpha }^{\beta }\right)}}={\displaystyle \frac{60+6{\alpha }^{\beta }}{74{\alpha }^{\beta }+15{\alpha }^{2\beta }+{\alpha }^{3\beta }+120}} \end{gathered}$$

(26)

To sum up, PWMs method for estimating the parameters using $M_{\mathrm{1,0},1}$  and $M_{\mathrm{1,1},0}$  :

Step 1: Calculate the population PWMs for the order p=1

$$M_{\mathrm{1,0},1}={\displaystyle \frac{360+108{\alpha }^{\beta }}{1044{\alpha }^{\beta }+580{\alpha }^{2\beta }+155{\alpha }^{3\beta }+20{\alpha }^{4\beta }+{\alpha }^{5\beta }+720}}$$

(27)

$$M_{\mathrm{1,1},0}={\displaystyle \frac{60+6{\alpha }^{\beta }}{74{\alpha }^{\beta }+15{\alpha }^{2\beta }+{\alpha }^{3\beta }+120}}$$

(28)

Step 2: Calculate the estimated sample PWMs whether the unbiased or biased estimators and equate these estimators with the corresponding population PWMs.

$${\widehat{\alpha }}_{r=1}={\displaystyle \frac{1}{n}}\sum _{i=1}^{n-1}{\displaystyle \frac{\left(n-i\right)}{\left(n-1\right)}} y_{i:n}=M_{\mathrm{1,0},1}, or {\widehat{\alpha }}_r={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^0{\left(1-p_{i:n}\right)}^s=M_{\mathrm{1,0},1}$$

$${\widehat{\beta }}_{s=1}={\displaystyle \frac{1}{n}}\sum _{i=1+s}^n{\displaystyle \frac{\left(i-1\right)}{\left(n-1\right)}} y_{i:n}=M_{\mathrm{1,1},0} , or {\widehat{\beta }}_s={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^1{\left(1-p_{i:n}\right)}^0=M_{\mathrm{1,1},0}$$

Step 3: The above equations construct a system of equations to be solved numerically. In this paper, the author used the Levenberg-Marquardt (LM) algorithm. The objective functions to be minimized are equations (27-28). Differentiate the previous equations (27-28) with respect to alpha and beta.

The Jacobian matrix is $\left[\begin{array}{cc}{\displaystyle \frac{\partial }{\partial \alpha }}M\_101 & {\displaystyle \frac{\partial }{\partial \beta }}M\_101\\ {\displaystyle \frac{\partial }{\partial \alpha }}M\_110& {\displaystyle \frac{\partial }{\partial \beta }}M\_110\end{array}\right]$

$$\left(y-f\left({\mathit{\theta }}^{\left(\mathit{a}\right)}\right)\right)=\left[\begin{array}{c}{\widehat{\alpha }}_r-{\displaystyle \frac{360+108{\alpha }^{\beta }}{1044{\alpha }^{\beta }+580{\alpha }^{2\beta }+155{\alpha }^{3\beta }+20{\alpha }^{4\beta }+{\alpha }^{5\beta }+720}}\\ {\widehat{\beta }}_s-{\displaystyle \frac{60+6{\alpha }^{\beta }}{74{\alpha }^{\beta }+15{\alpha }^{2\beta }+{\alpha }^{3\beta }+120}}\end{array}\right], \& {\mathit{ }\mathit{\theta }}^{\left(\mathit{a}\right)}=\left[\begin{array}{c}{\alpha }^{\left(a\right)}\\ {\beta }^{\left(a\right)}\end{array}\right]$$

$$\mathrm{Apply} \mathrm{LM} \mathrm{algorithm}: {\theta }^{\left(a+1\right)}={\theta }^{\left(a\right)}+{\left[J^{\text{'}}\left({\theta }^{\left(a\right)}\right)J\left({\theta }^{\left(a\right)}\right)+{\lambda }^{\left(a\right)}{I}^{\left(a\right)}\right]}^{-1}\left[J\text{'}\left({\theta }^{\left(a\right)}\right)\right]\left(y-f\left({\theta }^{\left(a\right)}\right)\right).$$

where the parameters used in the first iteration are the initial guess, then they are updated according to the sum of squares of errors. The LM algorithm is an iterative algorithm. $J^{\text{'}}\left({\theta }^{\left(a\right)}\right)$  is the Jacobian function which is the first derivative of the objective function evaluated at the initial guess ${\theta }^{\left(a\right)}$  and ${\lambda }^{\left(a\right)}$  is a damping factor that adjusts the step size in each iteration direction. The starting value of this factor is usually 0.001 and according to the sum square of errors (SSE) in each iteration this damping factor is adjusted:

$${SSE}^{\left(a+1\right)}\ge {SSE}^{\left(a\right)} so update: {\lambda }_{updated}=10\ast {\lambda }_{old}$$

${SSE}^{\left(a+1\right)}< {SSE}^{\left(a\right)} so update: {\lambda }_{updated}={\displaystyle \frac{1}{10}}\ast {\lambda }_{old}$ , $f\left({\theta }^{\left(a\right)}\right)$ is the objective function (population PWMs, $M_{\mathrm{1,0},1}$ & $M_{\mathrm{1,1},0 }$) evaluated at the initial guess, and $y$ is the sample estimates of population PWMs , the $M_{\mathrm{1,0},1}$ and $M_{\mathrm{1,1},0 }$. Steps of the LM algorithm:

• Start with the initial guess of parameters (alpha and beta).
• Substitute these values in the objective function and the Jacobian.
• Choose the damping factor, say lambda=0.001
• Substitute in the equation (LM equation) to get the new parameters.
• Calculate the SSE at these parameters and compare this SSE value with the previous one when using the initial parameters to adjust for the damping factor.
• Update the damping factor accordingly as previously explained.
• Start new iteration with the new parameters and the new updated damping factor, i.e, apply the previous steps many times till convergence is achieved or a pre-specified number of iterations is accomplished.

The value of this quantity: ${\left[J^{\text{'}}\left({\theta }^{\left(a\right)}\right)J\left({\theta }^{\left(a\right)}\right)+ {\lambda }^{\left(a\right)}{I}^{\left(a\right)}\right]}^{-1}$ can be considered a good approximation to the variance–covariance matrix of the estimated parameters. Standard errors for the estimated parameters are the square root of the diagonal of the elements in this matrix.

2.2.3. Calculating $M_{\mathrm{1,0},2}$ Equation (29)

Substitute s=2 in equation (13) and solve to get equation (29):

$$\begin{gathered} {{\displaystyle \frac{6}{{\alpha }^{\beta }}}\{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{2}{0}}\right){\left(3\right)}^0 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{0}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(0\right)+0+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{0}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(0\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{0}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(0\right)+2+{\alpha }^{\beta }}}\right)\right] \\ {+\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{2}{1}}\right){\left(3\right)}^1 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+0+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{1}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(1\right)+2+{\alpha }^{\beta }}}\right)\right] \\ {+\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{2}{2}}\right){\left(3\right)}^2 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{2}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+0+{\alpha }^{\beta }}}\right)+{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{2}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{2}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+2+{\alpha }^{\beta }}}\right)\right]\} \\ {-{\displaystyle \frac{6}{{\alpha }^{\beta }}}\{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{2}{0}}\right){\left(3\right)}^0 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{0}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(0\right)+0+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{0}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(0\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{0}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(0\right)+2+{\alpha }^{\beta }}}\right)\right] \\ {+\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{2}{1}}\right){\left(3\right)}^1 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{1}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+0+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{1}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(1\right)+2+{\alpha }^{\beta }}}\right)\right] \\ {+\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{2}{2}}\right){\left(3\right)}^2 \left[{\left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{2}{0}}\right){\left({\displaystyle \frac{2}{3}}\right)}^0\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+0+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{2}{1}}\right){\left({\displaystyle \frac{2}{3}}\right)}^1 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+1+{\alpha }^{\beta }}}\right) \\ +{\left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{2}{2}}\right){\left({\displaystyle \frac{2}{3}}\right)}^2 \left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+2+{\alpha }^{\beta }}}\right)\right]\} \\ {M}_{\mathrm{1,0},2}={\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{{\alpha }^{\beta }}{2+{\alpha }^{\beta }}}-{\displaystyle \frac{6{\alpha }^{\beta }}{4+{\alpha }^{\beta }}}+{\displaystyle \frac{{4\alpha }^{\beta }}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{9{\alpha }^{\beta }}{6+{\alpha }^{\beta }}}-{\displaystyle \frac{{12 \alpha }^{\beta }}{7+{\alpha }^{\beta }}}+{\displaystyle \frac{{4\alpha }^{\beta }}{8+{\alpha }^{\beta }}}\right\} \\ -{\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{{\alpha }^{\beta }}{3+{\alpha }^{\beta }}}-{\displaystyle \frac{6{\alpha }^{\beta }}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{{4\alpha }^{\beta }}{6+{\alpha }^{\beta }}}+{\displaystyle \frac{9{\alpha }^{\beta }}{7+{\alpha }^{\beta }}}-{\displaystyle \frac{{12 \alpha }^{\beta }}{8+{\alpha }^{\beta }}}+{\displaystyle \frac{{4\alpha }^{\beta }}{9+{\alpha }^{\beta }}}\right\} \\ {M}_{\mathrm{1,0},2}=6\left\{{\displaystyle \frac{1}{2+{\alpha }^{\beta }}}-{\displaystyle \frac{6}{4+{\alpha }^{\beta }}}+{\displaystyle \frac{4}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{9}{6+{\alpha }^{\beta }}}-{\displaystyle \frac{12}{7+{\alpha }^{\beta }}}+{\displaystyle \frac{4}{8+{\alpha }^{\beta }}}\right\} \\ -6\left\{{\displaystyle \frac{1}{3+{\alpha }^{\beta }}}-{\displaystyle \frac{6}{5+{\alpha }^{\beta }}}+{\displaystyle \frac{4}{6+{\alpha }^{\beta }}}+{\displaystyle \frac{9}{7+{\alpha }^{\beta }}}-{\displaystyle \frac{12}{8+{\alpha }^{\beta }}}+{\displaystyle \frac{4}{9+{\alpha }^{\beta }}}\right\} \\ {M}_{\mathrm{1,0},2}={\displaystyle \frac{K}{L}} where : \\ K=58320 { \alpha }^{\beta } +6480 {\alpha }^{2\beta }+120960 \end{gathered}$$

(29)

$$L=663696{ \alpha }^{\beta }+509004 {\alpha }^{2\beta }+214676 {\alpha }^{3\beta }+54649{\alpha }^{4\beta }+8624 {\alpha }^{5\beta }+826 {\alpha }^{6\beta }+44 {\alpha }^{7\beta }+{\alpha }^{8\beta }+362880$$

2.2.4. Calculating $M_{\mathrm{1,2},0}$ Equation (30)

Substitute r=2 in equation (21) and solve to get equation (30):

$$\begin{gathered} M_{\mathrm{1,2},0} ={\displaystyle \frac{6}{{\alpha }^{\beta }}} \sum _{i=0}^2{3^r \left(-1\right)}^i\left({\displaystyle \genfrac{}{}{0pt}{}{r}{i}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^i\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2r+i+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2r+i+{\alpha }^{\beta }}}\right)\right] \\ ={\displaystyle \frac{6}{{\alpha }^{\beta }}} \{{3^2 \left(-1\right)}^0\left({\displaystyle \genfrac{}{}{0pt}{}{2}{0}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^0\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+\left(0\right)+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+\left(0\right)+{\alpha }^{\beta }}}\right)\right] \\ {{+ 3}^2 \left(-1\right)}^1\left({\displaystyle \genfrac{}{}{0pt}{}{2}{1}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^1\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+\left(1\right)+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+\left(1\right)+{\alpha }^{\beta }}}\right)\right] \\ {{+ 3}^2 \left(-1\right)}^2\left({\displaystyle \genfrac{}{}{0pt}{}{2}{2}}\right) {\left({\displaystyle \frac{2}{3}}\right)}^2\left[\left({\displaystyle \frac{{\alpha }^{\beta }}{2+2\left(2\right)+\left(2\right)+{\alpha }^{\beta }}}\right)-\left({\displaystyle \frac{{\alpha }^{\beta }}{3+2\left(2\right)+\left(2\right)+{\alpha }^{\beta }}}\right)\right]\} \\ {M}_{\mathrm{1,2},0}={\displaystyle \frac{6}{{\alpha }^{\beta }}}\left\{{\displaystyle \frac{9{\alpha }^{\beta }}{6+{\alpha }^{\beta }}}-{\displaystyle \frac{{12\alpha }^{\beta }}{7+{\alpha }^{\beta }}}+{\displaystyle \frac{{4\alpha }^{\beta }}{8+{\alpha }^{\beta }}}-{\displaystyle \frac{{9\alpha }^{\beta }}{7+{\alpha }^{\beta }}}+ {\displaystyle \frac{{12\alpha }^{\beta }}{8+{\alpha }^{\beta }}}-{\displaystyle \frac{{4\alpha }^{\beta }}{9+{\alpha }^{\beta }}}\right\} \\ {M}_{\mathrm{1,2},0}=6\left\{{\displaystyle \frac{9}{6+{\alpha }^{\beta }}}-{\displaystyle \frac{21}{7+{\alpha }^{\beta }}}+{\displaystyle \frac{16}{8+{\alpha }^{\beta }}}-{\displaystyle \frac{{4\alpha }^{\beta }}{9+{\alpha }^{\beta }}}\right\} \\ {M}_{\mathrm{1,2},0}={\displaystyle \frac{ 150 {\alpha }^{\beta }+6 {\alpha }^{2\beta } + 1008 }{1650 {\alpha }^{\beta }+335 {\alpha }^{2\beta } + 30 {\alpha }^{3\beta } + {\alpha }^{4\beta }+ 3024}} \end{gathered}$$

(30)

In summary, PWMs method for estimating the parameters using $M_{\mathrm{1,0},2}$ and $M_{\mathrm{1,2},0}$ can be explained in the following steps: the first step is to calculate the population PWMs for the order p=1, using equations (29-30). The second step is to calculate the estimated sample PWMs, whether the unbiased or biased estimators, and equate these estimators with the corresponding population PWMs.

$${\widehat{\alpha }}_{r=2}={\displaystyle \frac{1}{n}}\sum _{i=1}^{n-r}{\displaystyle \frac{\left(n-i\right)\left(n-i-1\right)}{\left(n-1\right)\left(n-2\right)}} y_{i:n}=M_{\mathrm{1,0},2} or {\widehat{\alpha }}_r={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^0{\left(1-p_{i:n}\right)}^2=M_{\mathrm{1,0},2}$$

$${\widehat{\beta }}_{s=2}={\displaystyle \frac{1}{n}}\sum _{i=1+s}^n{\displaystyle \frac{\left(i-1\right)\left(i-2\right)}{\left(n-1\right)\left(n-2\right)}} y_{i:n}=M_{\mathrm{1,2},0} or {\widehat{\beta }}_s={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_{i:n}{\left(p_{i:n}\right)}^2{\left(1-p_{i:n}\right)}^0=M_{\mathrm{1,2},0}$$

The last step is to construct a system of the above equations to solve it numerically. Using the Levenberg-Marquardt (LM) algorithm to minimize equations (29) and (30). Differentiate the previous equations (29-30) with respect to alpha and beta

$$\mathrm{The} \mathrm{Jacobian} \mathrm{matrix} \mathrm{is} \left[\begin{array}{cc}{\displaystyle \frac{\partial }{\partial \alpha }}M\_102 & {\displaystyle \frac{\partial }{\partial \beta }}M\_102\\ {\displaystyle \frac{\partial }{\partial \alpha }}M\_120& {\displaystyle \frac{\partial }{\partial \beta }}M\_120\end{array}\right]$$

$$\left(y-f\left({\mathit{\theta }}^{\left(\mathit{a}\right)}\right)\right)=\left[\begin{array}{c}{\widehat{\alpha }}_r-{\displaystyle \frac{K}{L}}\\ {\widehat{\beta }}_s-{\displaystyle \frac{ 150 {\alpha }^{\beta }+6 {\alpha }^{2\beta }+1008 }{1650 {\alpha }^{\beta }+335 {\alpha }^{2\beta }+30 {\alpha }^{3\beta }+{\alpha }^{4\beta }+3024}}\end{array}\right], \& {\mathit{\theta }}^{\left(\mathit{a}\right)}=\left[\begin{array}{c}{\alpha }^{\left(a\right)}\\ {\beta }^{\left(a\right)}\end{array}\right]$$

Apply the LM algorithm : ${\theta }^{\left(a+1\right)}= {\theta }^{\left(a\right)}+ {\left[J^{\text{'}}\left({\theta }^{\left(a\right)}\right)J\left({\theta }^{\left(a\right)}\right)+ {\lambda }^{\left(a\right)}{I}^{\left(a\right)}\right]}^{-1}\left[J\text{'}\left({\theta }^{\left(a\right)}\right)\right]\left(y-f\left({\theta }^{\left(a\right)}\right)\right)$.

See Appendix A for the Jacobian matrices for both $M_{\mathrm{1,0},1} , M_{\mathrm{1,1},0}$ and $M_{\mathrm{1,0},2} , M_{\mathrm{1,2},0}$

2.3. Asymptotic Distribution of PWM Estimators

For optimal evaluations of estimators derived from various Probability Weighted Moments (PWMs) method, it is essential to leverage analytical expressions for their variances, as established by asymptotic theory. This approach not only streamlines the process but also eliminates the need for extensive and potentially flawed computer simulations. Unfortunately, in numerous cases, deriving straightforward expressions for the asymptotic variance of moments and parameter estimators proves to be a significant challenge. Moreover, as highlighted by Hosking1986 and further explored by [4], the reliability of asymptotic variance expressions diminishes, particularly when working with small sample sizes.

In his 1986 study, Hosking offered critical insights into the asymptotic variances of generalized Pareto parameters estimated through classical PWMs. These asymptotic variances can be misleading, significantly misrepresenting true sampling variances when the sample size (n) is below 50. The asymptotic distribution of sample PWMs manifests as a linear combination of order statistics. [23] proved that the vector of $\widehat{b}=({\widehat{b}}_1,{\widehat{b}}_2)$ has asymptotically a multivariate Normal Distribution with mean $\widehat{b}=({\widehat{b}}_1,{\widehat{b}}_2)$ and covariance matrix $n^{-1}V$ . Supporting this, [24] employed the delta method to show that the covariance matrix can be articulated in terms of the variance of the parameter being estimated, hence, using the delta method, the covariance matrix has the expression of $n^{-1} GV{G}^T,$ where V is the variance of the parameter. The quantity ${\left[J^{\text{'}}\left({\theta }^{\left(a\right)}\right)J\left({\theta }^{\left(a\right)}\right)+ {\lambda }^{\left(a\right)}{I}^{\left(a\right)}\right]}^{-1}$ derived from applying the LM algorithm for parameter estimation serves as a reliable approximation of the variance–covariance matrix for the estimated parameters, known as the V matrix. Moreover, the G matrix, integral to the LM algorithm, functions as the Jacobian matrix, reflecting the relationships between the parameters being estimated. Harnessing these methods can vastly improve the accuracy and reliability of estimators, ultimately leading to more credible and impactful results in statistical analysis.

Using Taylor series expansion of the $g_1\left(\theta \right)=M_{101}$ and $g_2\left(\theta \right)=M_{110}$ around the ${\theta }_0$ gives

$g\left(\widehat{\theta }\right)=g\left({\theta }_0\right)+\left(\widehat{\theta }-{\theta }_0\right)g^{\text{'}}\left({\theta }_0\right)$ . Therefore, the delta method can be applied. Taking the variance on both sides yields $var\left[g\left(\widehat{\theta }\right)\right]=\left[g^{\text{'}}\left({\theta }_0\right)\right]var\left(\widehat{\theta }\right){\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T$ where $g^{\text{'}}\left({\theta }_0\right)={\left.{\displaystyle \frac{\partial g\left({\theta }_0\right)}{\partial {\theta }_0}}\right|}_{{\theta }_0=\widehat{\theta }}$ which is the Jacobian matrix in LM algorithm. Hence, $var\left(\widehat{\theta }\right)={\left[{\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]\right]}^{-1}$ and $var\left[g\left(\widehat{\theta }\right)\right]=\left[\begin{array}{cc}var\left(M_{101}\right)& cov\left(M_{101},M_{110}\right)\\ cov\left(M_{101},M_{110}\right)& var\left(M_{110}\right)\end{array}\right]$. This matrix, ${\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]$ may have a high condition number that inflates the $var\left(\widehat{\theta }\right)$ .

To ameliorate this, an appropriate regularization factor$\left(\xi \right)$ can be added to the matrix like this ${\left[{\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]+ \xi I\right]}^{-1}$, where $\xi ={\displaystyle \frac{{\lambda }_{max}-C^{\ast }{\lambda }_{min}}{C^{\ast }-1}}$ such that ${\lambda }_{max}$ is the maximum eigenvalue of the matrix ${\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]$ and ${\lambda }_{min}$ is the minimum eigenvalue for the same matrix. While the $C^{\ast }$ is the desired condition number. The desired value of this number is usually from 5 up to 10. Other types of variance is the value of this quantity ${\left[J^{\text{'}}\left({\theta }^{\left(a\right)}\right)J\left({\theta }^{\left(a\right)}\right)+ {\lambda }^{\left(a\right)}{I}^{\left(a\right)}\right]}^{-1}$ which can be considered a good approximation to the variance-covariance matrix of the estimated parameters $\widehat{\theta }$ . The delta method is applied to derive variance of the function of the estimated parameters.

The same discussion is applied to $M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}$. See Appendix B for the variance and covariance of the lower and higher order PWMs.

3. Monte Carlo Simulation Study

3.1. Simulation Study for the Parameters

A simulation was conducted by generating 1000 replicates (N) of the MBUW random variable $y$, each replicate with different sample sizes (n). The author chosen suitable pairs of $\alpha$ and $\beta$ parameters respectively as follows (0.5, 0.6), (0.5, 1.3), (1.1, 0.6), and (1.1, 1.6) . The sample sizes were 20, 50, 100, and 500. For each sample, the parameters were estimated using the PWMs method (the low and high orders), maximum likelihood (MLE) and the Method of Moment (MOM, the first and the second raw moments). For the MOM, $m_1^{\text{'}}\left(y\right)={\mu }_1^{\text{'}}$ and , $m_2^{\text{'}}\left(y\right)={\mu }_2^{\text{'}}$ where $m_1^{\text{'}}\left(y\right)={\displaystyle \frac{1}{n}}\sum _i^n{y}_i=\overline{y}$ and $m_2^{\text{'}}\left(y\right)={\displaystyle \frac{1}{n}}\sum _i^n{y}_i^2$. ${\mu }_1^{\text{'}}={\displaystyle \frac{6}{\left(2+{\alpha }^{\beta }\right)\left(3+{\alpha }^{\beta }\right)}}$ and ${\mu }_2^{\text{'}}={\displaystyle \frac{6}{\left(2+{2\alpha }^{\beta }\right)\left(3+2{\alpha }^{\beta }\right)}}$. For each method, where $\widehat{\theta }=\left[\begin{array}{c}\widehat{\alpha }\\ \widehat{\beta }\end{array}\right]$, the statistical indices used for the comparisons are:

• $Average Absolute Bias \left(AAB\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N\left|\widehat{\mathit{\theta }}-\mathit{\theta }\right|$
• $Mean Square Error \left(MSE\right)= {\displaystyle \frac{1}{N}}\sum _{j=1}^N{\left(\widehat{\mathit{\theta }}-\mathit{\theta }\right)}^2$
• $Mean Relative Error \left(MRE\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N{\displaystyle \frac{\left|\widehat{\mathit{\theta }}-\mathit{\theta }\right|}{\theta }}$

$The empirical variance for each parameter= MSE-{bias}^2$.

The 2.5^th and the 97.5^th quantiles of the sampling distribution for each parameter are also recorded. The number of the valid samples, which shows GoF, is recorded. These results are illuminated in Tables (1-4). Each table shows the results of the comparisons between the methods for each sample size. The method yielding the estimator with the least bias is colored. Tables (5-8) show the empirical variance-covariance matrices obtained from the sampling distribution of the parameters. These last four tables demonstrate the Pearson correlation coefficient between the parameters, the condition number and eigenvalues of the empirical covariance matrix. Table (9) expounds the variance decay rate for each simulation validating the asymptotic consistency of the estimator $\left({\widehat{\theta }}_n\right)$as n increases. In asymptotic theory, for a consistent estimator ${\widehat{\theta }}_n$ of a parameter $\theta$, the variance typically behaves like $var\left({\widehat{\theta }}_n\right)\approx {\displaystyle \frac{C}{n}}$ for large n, where C is some constant. This means that the variance decays inversely with sample size, it decays as 1/n. Taking log of both sides of $var\left({\widehat{\theta }}_n\right)\approx {\displaystyle \frac{C}{n}}$ gives $\log \left[var\left({\widehat{\theta }}_n\right)\right]\approx \log \left(C\right)-\log n$ so the slope in the log-log plot of the variance vs. sample size should be approximately (-1). For each simulation the last two values of the variance of each parameter at n=500 and n=100 will be substituted in this log equation, in other words, the slope is calculated as

${\displaystyle \frac{\log \left[var\left({\widehat{\theta }}_{n=500}\right)\right]-\log \left[var\left({\widehat{\theta }}_{n=100}\right)\right]}{\log 500-\log 100}}\approx -1$ . This is calculated for $\left(\widehat{\alpha }\right)\&\left(\widehat{\beta }\right)$ . Tables (10-13) illustrate the computational time for each simulation and the relative efficiency of each method vs. other methods. The most efficient method (the one with the least MSE among other methods in the row, for $\left(\widehat{\alpha }\right)\&\left(\widehat{\beta }\right)$) are colored. The computational time is recorded in seconds per replicate and also the total time of all replicates is recorded. The relative efficiency of estimator A obtained from method A to estimator B gained from method B is defined as $RE( A vs. B) ={\displaystyle \frac{{\widehat{MSE}}_B}{{\widehat{MSE}}_A}}$, this means that the estimator A is more efficient than the estimator B.

3.2. Simulation Study for the Function of the Parameters (Validation of the Delta Method)

Another simulation was conducted to validate the asymptotic normality of the delta method as follows: for sample size, n=500, alpha=0.5 , beta=1.3, and N=1000

Step 1: for each replicate the $g_1\left(\widehat{\theta }\right)={\widehat{M}}_{101}$ and $g_2\left(\widehat{\theta }\right)={\widehat{M}}_{110}$ .

Step 2: $g_1\left({\theta }_0\right)=M_{101}$ and $g_2\left({\theta }_0\right)=M_{110}$ where ${\theta }_0=\left[\begin{array}{cc}\alpha =0.5,& \beta =1.3\end{array}\right]$ .

Step 3: ${error}_{101}={\widehat{M}}_{101}-M_{101}$and ${error}_{110}={\widehat{M}}_{110}-M_{110}$Step 4: calculate the variance-covariance between ${error}_{101}$ and ${error}_{110}$, say it is $\left(co{v}_{emp}\right)$. This is the empirical covariance matrix (2by2 matrix).

Step 5: obtain the theoretical covariance matrix $\left(co{v}_{theo}\right)$ using the delta method applied on each replicate then take the mean, let’s call this $mean\_co{v}_{theo}$Step 6: apply the singular value decomposition (SVD) on the error matrix $error=\left[\begin{array}{cc}{error}_{101}& {error}_{110}\end{array}\right]$ . The error matrix is N by 2 matrix.

$\left[\begin{array}{cc}{error}_{101}& {error}_{110}\end{array}\right]=U_1\ast S_1\ast V_1^T$, where $S_1$ is the large non-zero singular value, and $V_1$ is the corresponding principal vector, ${\lambda }_{empirical}=S_1^2/\left(N-1\right)$. Let’s call $s_{emp}=\left[\begin{array}{cc}{error}_{101}& {error}_{110}\end{array}\right]\ast V_1$, is the projection of the error matrix in the direction of the principle vector. Let’s call $Z_{empirical}=s_{emp}/\sqrt{{\lambda }_{empirical}}$ , this standardized empirical Z should be distributed as standard normal. These Zs’ are the standardized errors after projection along the direction of maximum variability.

Step 7: get the trace of the $co{v}_{emp}$ and the trace of the $Mean\_co{v}_{theo}$ . Let’s call $scalar=trace\left(co{v}_{emp.}\right)/trace\left(mean\_co{v}_{theo}\right)$. Let’s call ${cov}_{theo\_scaled}=scalar\ast mean\_co{v}_{theo} .$Step 8: apply the singular value decomposition (SVD) on $mean\_co{v}_{theo}$, so the $mea{n}_{co{v}_{theo}}=U_2\ast S_2\ast V_2^T$ where $S_2$ is the large non-zero singular value, and $V_2$ is the corresponding principal vector, Let’s call $s_{theo}=\left[\begin{array}{cc}{error}_{101}& {error}_{110}\end{array}\right]\ast V_2$, is the projection of the error matrix in the direction of the maximum variability of the theoretical covariance matrix.

Step 9: calculate ${\lambda }_{theo\_scaled}=scalar\ast V_2\ast mea{n\_cov}_{theo}\ast V_2^T$ . Then calculate $Z_{theo\_scaled}=s_{theo}/\sqrt{{\lambda }_{theo\_scaled}}$. This standardized theoretical Z should be distributed as standard normal.

The intuition behind this algorithm is that the the theoretical covariance matrix (obtained from the delta method) should approximate the empirical covariance matrix from the replicates. The empirical and the theoretical standardized error can be defined respectively as

$Z_{emp}={\mathsf{\Sigma }}_{emp}^{-1/2}\left(g\left(\widehat{\theta }\right)-g\left({\theta }_0\right)\right)$ and $Z_{theo\_scaled}={\mathsf{\Sigma }}_{theo\_scaled}^{-1/2}\left(g\left(\widehat{\theta }\right)-g\left({\theta }_0\right)\right)$ and to hold for asymptotic normality these Zs’ should be distributed as standard normal. A consistent observation across simulations (for all pairs of parameter and at different sample sizes) was that both covariance matrices (empirical and theoretical) shared nearly identical eigenvectors, as revealed by the singular value decomposition(SVD). This indicates that the empirical and theoretical covariance matrices capture a similar dependence structure between the two parameters. However, notable differences were observed in their singular values, particularly for large sample size. The theoretical matrix was a full rank for small sample size (n=20) but it is rank-one as n increases (n=500). This degeneration reflects the increasing linear dependency between the two parameter estimators as the sample size grows, implying that one parameter becomes an almost deterministic function of the other in large samples. To adjust for the overall discrepancy between the empirical and the theoretical covariance, a scalar rescaling factor was applied to the estimated theoretical matrix. This factor is $trace\left(co{v}_{emp.}\right)/trace\left(mean\_co{v}_{theo}\right)$ , where the trace of covariance matrix represents the total marginal variance. Multiplying the theoretical covariance matrix by this scalar aligns the overall variance magnitude of the theoretical covariance with that of the empirical one while preserving its correlation structure. In other words, the rescaling equalizes the total variability but does not alter the orientation (the dependence structure, the eigenvectors) between the parameters. After this adjustment,; the standardized vector, the theoretical Z, exhibited approximately zero mean and unit standard deviation, confirming that the scaling effectively normalized the parameter estimators (function of parameters) .Hence, the asymptotic normality of the delta method was thus validated: as n increased, the distribution of the standardized estimators (PWMs are functions of the parameters) approached the standard normal, and the theoretical covariance structure aligned closely with the empirical results, despite the asymptotic rank deficiency inherent in the parameter dependence. Figure 3 illustrates the histogram, the fitted standard normal curves, and the QQ plot for the Z empirical and the Z theoretical (scaled).

4. Results and Discussion of the Simulation Study

Table 1, Table 2, Table 3 and Table 4 show that at various sample sizes where $\left(\alpha =0.5 \& \beta =0.6\right)$ the estimated $\left(\widehat{\alpha }\right)$obtained by MLE has the least Average Absolute Bias (AAB) and the least Mean Square Error (MSE) among the methods while the estimated $\left(\widehat{\beta }\right)$obtained by $M_{\mathrm{1,0},2}{, M}_{\mathrm{1,2},0}$ has the least AAB and MSE among the methods. This is also true at simulation runs using the pair $\left(\alpha =0.5 \& \beta =1.3\right)$ and the pair $\left(\alpha =1.1 \& \beta =1.6\right)$ . However, at true values of the pair $\left(\alpha =1.1 \& \beta =0.6 \right)$, the higher order PWMs ($M_{\mathrm{1,0},2}{, M}_{\mathrm{1,2},0}$) produce $\left(\widehat{\alpha }\right)$and $\left(\widehat{\beta }\right)$ with the least AAB and MSE among the methods. The method yielding an estimate of each parameter with the least bias and least MSE is colored in each table. The tables record the number of valid samples which pass the GoF in each run of the simulation. Method of Moments (MOM) give non-identifiable estimates of the parameters at $\left(\alpha =1.1 \& \beta =0.6 \right)$ and $\left(\alpha =1.1 \& \beta =1.6\right)$ especially at small sample size; n=20 and n=50. Moreover, at large sample sizes; n=100 and n=500 the number of valid samples gained by (MOM) are less than the number of valid samples obtained by other methods at the same values of the parameters pair. MOMs almost give the highest AAB among other methods at different sample sizes and at different pairs of parameters. MOMs almost yield the highest MSE among other methods across different sample sizes and across different pairs of the simulated parameters. The AAB of the higher order PWMs is approximately equal to that obtained from the lower order PWMs especially at $\left(\alpha =0.5 \& \beta =1.3\right)$ and at almost all sample sizes. The standard error of each of the estimated parameter is recorded beside the estimated value of the parameter. It is also obvious that the AAB and the MSE decrease as the sample size increases. This is almost true for all methods and for both estimated parameters. And this supports that the methods yield consistent estimators for the parameters.

Table 5 shows the variance-covariance matrix for each estimated parameter $\left(\widehat{\alpha }\right)$ and $\left(\widehat{\beta }\right)$obtained by the lower order PWMs ($M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}$). These estimated parameters are perfectly correlated as seen in all simulation across different sample sizes and across different pairs. These matrices show high condition number and one eigenvalue that is nearly zero. This may lead to identifiability problem and inflate the estimated variance which impairs the confidence interval construction. The lower order PWMs are derived from the order statistics so correlation is inherent in the method. The question is does the distribution contribute to this high dependency? To answer this question; the correlation between the estimated parameters obtained from deploying other methods were investigated using the sampling distribution gained from the simulation. Table 6 depicts the sampling distribution obtained from applying the higher order PWMs ($M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}$) for parameters estimation. The same outcomes were observed. But this may be due to the PWMs inheriting the high dependency from the order statistics. But as the same results were also returned from MLE and MOMs as shown in Table 7 and Table 8, this is a strong evidence that the distribution is the main root cause. Although, the method can contribute, but when every reasonable method shows the same correlation patterns, so the distribution itself is likely to cause this coupling. The Weibull distribution and any derived distribution from the Weibull, like the one between our hands (MBUW), the parameters control related features of the sample data (mean, spread, and tail) by giving information about some combinations of the parameters more strongly than about each parameter separately. Whether the dependency relation is linear, nonlinear, monotonic, or non- onotonic, the method can influence this. Figure 1 shows the relationship between the estimated $\left(\widehat{\alpha }\right)$ and estimated $\left(\widehat{\beta }\right)$ obtained from the simulation using the lower order PWMs at n=20. Table 6, Table 7 and Table 8show the sampling distribution of the estimated parameters obtained by the simulation study. The empirical variance-covariance matrices and the Pearson correlation coefficient are recorded. The condition number and the eigenvalues of the covariance matrix for each simulation are also shown. Other figures for the different shapes of the dependency are seen in Appendix C.

The importance to recognize that the cause of the correlation is the distribution itself is the application of the practical remedies to lessen this correlation and hence to deflate the variance and to construct a valid confidence interval. Some of these remedies is to reparameterize the parameter and use the log of the parameter. These remedies can be a blueprint for future studies.

Table 9. The variance decay of the estimated $\left(\widehat{\alpha }\right) \& \left(\widehat{\beta }\right)$ in different simulations.

	PWM: M101,M110 $\mathit{\alpha }=0.5$$\text{\&} \mathit{\beta }=0.6$	PWM: M101,M110 $\mathit{\alpha }=0.5$$\text{\&} \mathit{\beta }=1.3$	PWM: M101,M110 $\mathit{\alpha }=1.1$$\text{\&} \mathit{\beta }=0.6$	PWM: M101,M110 $\mathit{\alpha }=1.1$$\text{\&} \mathit{\beta }=1.6$
Slope For$\mathit{\alpha }$	$-0.9846$	$-1.0206$	$-1.0283$	$-1.0048$
Slope For$\mathit{\beta }$	$-0.9882$	$-1.0206$	$-1.0277$	$-0.9933$
	PWM: M102,M120 $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	PWM: M102,M120 $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	PWM: M102,M120 $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	PWM: M102,M120 $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For$\mathit{\alpha }$	$-0.9876$	$-1.006$	$-1.0204$	$-0.9932$
Slope For$\mathit{\beta }$	$-0.9820$	$-1.0061$	$-1.0244$	$-0.9990$
	MLE $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	MLE $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	MLE $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	MLE $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For α	$-0.9429$	$-1.0035$	$-0.9755$	$-09830$
Slope For$\mathit{\beta }$	$-0.9607$	$-0.7299$	$-0.6679$	$-0.1561$
	MOM $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	MOM $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	MOM $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	MOM $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For α	$-1.0047$	$-1.0627$	$-1.9940$	$-1.6376$
Slope For β	$-0.9940$	$-1.9704$	$-0.8719$	$-0.6409$

The slope that is not approximately equal to (-1) is colored

Table 9. The variance decay of the estimated $\left(\widehat{\alpha }\right) \& \left(\widehat{\beta }\right)$ in different simulations.

	PWM: M101,M110 $\mathit{\alpha }=0.5$$\text{\&} \mathit{\beta }=0.6$	PWM: M101,M110 $\mathit{\alpha }=0.5$$\text{\&} \mathit{\beta }=1.3$	PWM: M101,M110 $\mathit{\alpha }=1.1$$\text{\&} \mathit{\beta }=0.6$	PWM: M101,M110 $\mathit{\alpha }=1.1$$\text{\&} \mathit{\beta }=1.6$
Slope For$\mathit{\alpha }$	$-0.9846$	$-1.0206$	$-1.0283$	$-1.0048$
Slope For$\mathit{\beta }$	$-0.9882$	$-1.0206$	$-1.0277$	$-0.9933$
	PWM: M102,M120 $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	PWM: M102,M120 $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	PWM: M102,M120 $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	PWM: M102,M120 $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For$\mathit{\alpha }$	$-0.9876$	$-1.006$	$-1.0204$	$-0.9932$
Slope For$\mathit{\beta }$	$-0.9820$	$-1.0061$	$-1.0244$	$-0.9990$
	MLE $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	MLE $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	MLE $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	MLE $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For α	$-0.9429$	$-1.0035$	$-0.9755$	$-09830$
Slope For$\mathit{\beta }$	$-0.9607$	$-0.7299$	$-0.6679$	$-0.1561$
	MOM $\mathit{\alpha }=0.5$&$\mathit{\beta }=0.6$	MOM $\mathit{\alpha }=0.5$&$\mathit{\beta }=1.3$	MOM $\mathit{\alpha }=1.1$&$\mathit{\beta }=0.6$	MOM $\mathit{\alpha }=1.1$&$\mathit{\beta }=1.6$
Slope For α	$-1.0047$	$-1.0627$	$-1.9940$	$-1.6376$
Slope For β	$-0.9940$	$-1.9704$	$-0.8719$	$-0.6409$

The slope that is not approximately equal to (-1) is colored

Table (9) expounds the pattern of the variance decay. The slopes in this table are obtained from the empirical variances for the estimated parameters recorded in Tables (5-8). Using the log (variance) at n=500, log (variance) at n=100, log (n=500) and log (n=100), the slope can be calculated as previously explained in Section 3, (Monte Carlo Simulation Study). Figure 2 also supports this asymptotic variance decay at a rate (1/n) as n enlarges, where the estimated parameter $\left({\widehat{\theta }}_n\right)$ is asymptotically consistent at large sample size. Table (9) elucidates that the lower and the higher PWMs show better validity for this asymptotic consistency rather than the MLE and MOMs.

Table 10, Table 11, Table 12 and Table 13 show the relative efficiency of each method against other methods. The most efficient method among the others is colored. The higher order PWMs ($M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}$) are slightly more efficient than the lower order PWMs ($M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}$). This can be confirmed by calculating the relative efficiency of the lower order to the higher order PWMs. Across all the sample sizes, this RE at the pairs $\left(\alpha =0.5 \& \beta =0.6 \right)$ is in the range from 0.94 to 0.95 for both estimated parameters. While at the pairs $\left(\alpha =0.5 \& \beta =1.3 \right)$ , it is in the range from 0.998 to1. Furthermore, the RE at the pair $\left(\alpha =1.1 \& \beta =0.6 \right)$ is around 0.8. Additionally, the ER at the pairs $\left(\alpha =1.1 \& \beta =1.6 \right)$ is in the range from 0.94 to 0.95. So they are more or less close to each other if taking into account the computational time reported in seconds per replicate. At sample size 20, the time taken by the lower PWMs are slightly less than that taken by the higher order PWMs to accomplish the estimation. This is also true for sample size 50 except at the pair $\left(\alpha =0.5 \& \beta =1.3 \right)$ where the time taken by the lower order PWMs is 0.1289 seconds per replicate while this time is 0.1264 seconds per replicate for the higher order PWMs to complete the task of estimation. At sample size 100, the lower order PWMs take less time than the time of the higher order PWMs. At sample size 500, and at pairs $\left(\alpha =0.5 \& \beta =0.6 \right)$ and $\left(\alpha =0.5 \& \beta =1.3 \right)$, the lower order PWMs are faster than the higher order PWMs , in contrast to the pair $\left(\alpha =1.1 \& \beta =0.6 \right)$ and the pair $\left(\alpha =1.1 \& \beta =1.6 \right)$, where the higher order PWMs are faster than lower order PWMs. The time taken by both the lower and higher PWMs is linearly increasing across different sample sizes as n increases. The time taken by the MLE is less than that by PWMs of both order. The time taken by MOMs is the least among other methods. The time taken by MLE and MOMs does not vary markedly across sample sizes as seen in PWMs.

Table 10, Table 11, Table 12 and Table 13 almost depict that the descending order of the efficient methods from the most efficient to the least efficient in estimating $\left(\alpha \right)$ parameter across different sample sizes and across different pairs is MLE, $\left(M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}\right)$, $\left(M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}\right)$ then MOM. The exception for this is the pair $\left(\alpha =1.1 \& \beta =0.6 \right)$, which shows the descending order is $\left(M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}\right)$, MLE, $\left(M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}\right)$ and finally MOM across different sample sizes excluding n=500 & for this pair $\left(\alpha =1.1 \& \beta =0.6 \right)$ the order is $\left(M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}\right)$, $\left(M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}\right)$ , MLE then MOM. While the descending order of the efficient methods in estimating $\left(\beta \right)$ across different sample sizes and different pairs is $\left(M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}\right)$, $\left(M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}\right)$, MLE then MOM. The exception for this is the pair $\left(\alpha =0.5 \& \beta =0.6 \right)$, which shows the descending order is $\left(M_{\mathrm{1,0},2},M_{\mathrm{1,2},0}\right)$, $\left(M_{\mathrm{1,0},1},M_{\mathrm{1,1},0}\right)$, MOM and lastly MLE across different sample sizes. This last order is also true when n=500 and the pair is $\left(\alpha =0.5 \& \beta =1.3 \right)$.

Figure 3. Shows the histogram, fitted standard normal curve and QQ plot for the error derived from the difference between the function of the parameter at the estimated parameters and the true parameters; $g\left(\widehat{\theta }\right)-g\left({\theta }_0\right)$. The error standardized by the empirical covariance matrix obtained from the empirical distribution is a standard normal variable (left upper and lower panels). The same error standardized by the covariance matrix obtained from rescaling the theoretical covariance matrix obtained via applying the delta method on each replicate is also a standard normal variable. (right upper and lower panels). This is the results of the simulation using n=500 , N=1000 , $\alpha =0.5$ and$\beta =1.3$ and the lower order PWMs (M101,M110).

Figure 3. Shows the histogram, fitted standard normal curve and QQ plot for the error derived from the difference between the function of the parameter at the estimated parameters and the true parameters; $g\left(\widehat{\theta }\right)-g\left({\theta }_0\right)$. The error standardized by the empirical covariance matrix obtained from the empirical distribution is a standard normal variable (left upper and lower panels). The same error standardized by the covariance matrix obtained from rescaling the theoretical covariance matrix obtained via applying the delta method on each replicate is also a standard normal variable. (right upper and lower panels). This is the results of the simulation using n=500 , N=1000 , $\alpha =0.5$ and$\beta =1.3$ and the lower order PWMs (M101,M110).

Figure 3 shows the validation of asymptotic normality of the delta method. The error obtained from the difference between $g\left(\widehat{\theta }\right)$ and the $g\left({\theta }_0\right)$ where ${\theta }_0$ are the pair of parameter $\left(\alpha =0.5 \& \beta =1.3 \right)$ , the sample size is 500 and the function is the lower order PWMs , is standardized by the empirical covariance matrix . The same error is also standardized by the theoretical covarinace matrix obtained via applying the delta method. However, both matrices the empirical and the theoretical matrices are nearly singular whose one of their eigenvalues are nearly zero. So both matrices are subjected to singular value decomposition (SVD) so as to project the error along the principle vector with high variability (the vector corresponding to the large non-zero singular value). And to match the empirical with the theoretical matrix, the last one is multiplied by a small scalar that is the ratio between the trace of the empirical covariance and the trace of the theoretical covariance. Both standardized error should be a standard normal variable as shown in figure (3). This procedure was repeated for the higher order PWMs with different sample sizes and different pairs of the parameters. The same results were obtained which confirms the asymptotic normality of the delta method. The regularity conditions essential to apply the delta method is that the function should be a continuously differentiable function which also holds for the PWMs of both orders.

5. Real Data Analysis

The OECD platform is available at https://stats.oecd.org/index.aspx?DataSetCode=BLI

In OECD platform, many indicators are recorded. The author transformed these indicators into unit ratios and a conducted distributional fit for those indicators. In this paper, the author discussed three indicators; the water quality, the educational attainment and the self-reported health. Table 14 shows the descriptive analysis of the indicators. Figure 4 shows the boxplot of each indicator.

Water quality: is the percentage of the people who report being satisfied with the quality of their water. This reflects the proportion of people who have access to clean water supply. The values of this indicator are 0.92, 0.92, 0.79, 0.90, 0.62, 0.82, 0.87, 0.89, 0.93, 0.86, 0.97, 0.78, 0.91, 0.67, 0.81, 0.97, 0.8, 0.77, 0.77, 0.87, 0.82, 0.83, 0.83, 0.85, 0.75, 0.91, 0.85, 0.98, 0.82, 0.89, 0.81, 0.93, 0.76, 0.97 ,0.96, 0.62, 0.82, 0.88, 0.7, 0.62, 0.72.

Educational attainment: is presented as the percentage of a given population who has completed a specific level of education mainly the tertiary education. The values of this indicator are 0.84, 0.86, 0.8, 0.92, 0.67, 0.59, 0.43, 0.94, 82, 0.91, 0.91, 0.81, 0.86, 0.76, 0.86, 0.76, 0.85, 0.88, 0.63, 0.89, 0.89, 0.94, 0.74, 0.42, 0.81, 0.81, 0.82, 0.93, 0.55, 0.92, 0.90, 0.63, 0.84, 0.89, 0.42, 0.82, 0.92, 0.57, 0.95, 0.48

Self-reported health is presented as the percentage of the population that rates their health as good or very good and this reflects access to health services, healthy living conditions, and better life-style choices (diets and exercises). The values of this indicator are 0.85, 0.71, 0.74, 0.89, 0.60, 0.8, 0.73, 0.62, 0.7, 0.57, 0.68, 0.67, 0.66, 0.79, 0.58, 0.77, 0.84, 0.74, 0.73, 0.37, 0.34, 0.47, 0.46, 0.72, 0.66, 0.75, 0.86, 0.75, 0.6, 0.5, 0.65, 0.67, 0.75, 0.76, 0.81, 0.67, 0.73, 0.88, 0.43

The variables exhibit left skewness with different degrees. Educational attainment shows significant left skewness then the self-reported health and lastly the water quality indicator. Water quality shows mesokurtic shape while the educational attainment and the self-reported health show slightly more than excess kurtosis (leptokurtic, the kurtosis coefficient is more than 3). So the data exhibit different shapes.

The different unit distributions tested was Beta, Kumaraswamy, unit Lindley, Topp Leone and the new MBUW distribution. MLE method was applied using the Nelder Mead optimizer in MATLAB. For each distribution, the estimated parameters, variance, Log-Likelihood, AIC, CAIC, BIC, and HQIC are recorded. Also the KS-test, AD, CVM test are reported. Figures for the fitted CDFs and fitted PDFs are shown.

Beta Distribution:$f\left(y;\alpha ,\beta \right)={\displaystyle \frac{\mathsf{\Gamma }\left(\alpha +\beta \right)}{\mathsf{\Gamma }\left(\alpha \right)\mathsf{\Gamma }\left(\beta \right)}}{y}^{\alpha -1}{\left(1-y\right)}^{\beta -1} , 0<y<1 , \alpha >0, \beta >0$

Kumaraswamy Distribution:$f\left(y;\alpha ,\beta \right)=\alpha \beta y^{\alpha -1}{\left(1-y^{\alpha }\right)}^{\beta -1} , 0<y<1 , \alpha >0, \beta >0$

Topp-Leone Distribution:$f\left(y;\theta \right)=\theta \left(2-2y\right) {\left(2y-y^2\right)}^{\theta -1} , 0<y<1 , \theta >0$

Unit-Lindley: $f\left(y;\theta \right)={\displaystyle \frac{{\theta }^2}{1+\theta }}{\left(1-y\right)}^3\exp \left({\displaystyle \frac{-\theta y}{1-y}}\right) , 0<y<1 , \theta >0$

Table 15 shows that all the unit distributions fit the water quality data well. The variance-covariance obtained after fitting the MBUW distribution is not identified, although the other statistical indices like AIC, CAIC, BIC, HQIC, LL are comparable to other distributions. Figure 5 and Figure 6 show that the fitted CDF and the fitted PDF of the MBUW align with the fitted CDFs and fitted PDFs of both the Beta distribution and the Kumaraswamy distribution. But the problem is mainly the variance obtained from employing the MLE method utilizing Nelder Mead optimizer in MATLAB. Using the values of the estimated parameter as initial guesses to substitute in the LM algorithm and PWMs (lower and higher orders) gives the results shown in Table 18 and in Figure 11 and Figure 12.

Table 16. Estimators and validation indices for the educational attainment dataset.

	Beta		Kumaraswamy		MBUW		Topp-Leone	Unit-Lindley
theta	$\mathit{\alpha }=6.0373$		$\mathit{\alpha }=5.4735$		$\mathit{\alpha }=0.4333$		12.3585	0.2891
	$\mathit{\beta }=1.7313$		$\mathit{\beta }=1.9272$		$\mathit{\beta }=1.3477$
Var-cov	2.9116	0.9141	1.0268	0.3651			3.8183	0.0011
	0.9141	0.3410	0.3651	0.2227
SE	1.7063		1.0133				1.9540	0.0332
	0.5839		0.4719
AIC	-48.1413		-48.8385		-48.0574		-40.9642	-57.3738
CAIC	-47.817		-48.5142		-47.7331		-40.8589	-57.2685
BIC	-44.7636		-45.4607		-44.6796		-39.2753	-55.6849
HQIC	-46.92		-47.6172		-46.8361		-40.3535	-56.7631
LL	26.0707		26.4193		26.0287		21.4821	29.6869
K-S Value	0.1481		0.1434		0.1577		0.2598	0.0697
H0	Fail to reject		Fail to reject		Fail to reject		Reject	Fail to reject
P-value	0.1613		0.1841		0.1219		0.0023	0.9442
AD	1.1861		1.1575		1.4386		4.5709	0.2596
CVM	0.2068		0.1985		0.2546		0.8485	0.0344

Table 16. Estimators and validation indices for the educational attainment dataset.

	Beta		Kumaraswamy		MBUW		Topp-Leone	Unit-Lindley
theta	$\mathit{\alpha }=6.0373$		$\mathit{\alpha }=5.4735$		$\mathit{\alpha }=0.4333$		12.3585	0.2891
	$\mathit{\beta }=1.7313$		$\mathit{\beta }=1.9272$		$\mathit{\beta }=1.3477$
Var-cov	2.9116	0.9141	1.0268	0.3651			3.8183	0.0011
	0.9141	0.3410	0.3651	0.2227
SE	1.7063		1.0133				1.9540	0.0332
	0.5839		0.4719
AIC	-48.1413		-48.8385		-48.0574		-40.9642	-57.3738
CAIC	-47.817		-48.5142		-47.7331		-40.8589	-57.2685
BIC	-44.7636		-45.4607		-44.6796		-39.2753	-55.6849
HQIC	-46.92		-47.6172		-46.8361		-40.3535	-56.7631
LL	26.0707		26.4193		26.0287		21.4821	29.6869
K-S Value	0.1481		0.1434		0.1577		0.2598	0.0697
H0	Fail to reject		Fail to reject		Fail to reject		Reject	Fail to reject
P-value	0.1613		0.1841		0.1219		0.0023	0.9442
AD	1.1861		1.1575		1.4386		4.5709	0.2596
CVM	0.2068		0.1985		0.2546		0.8485	0.0344

Table 16 shows that all the unit distributions fit the educational attainment data well except the Topp-Leone distribution. The same trouble of the variance-covariance obtained after fitting the MBUW distribution is shown, even though the other statistical indices like AIC, CAIC, BIC, HQIC, LL are more or less analogous to other fitting distributions. Figure 7 and Figure 8 show that the fitted CDF and the fitted PDF of the MBUW is side by side with the fitted CDFs and fitted PDFs of both the Beta distribution and the Kumaraswamy distribution. But then again; the trouble is essentially the variance gained from deploying the MLE method consuming Nelder Mead optimizer in MATLAB. Treating the values of the estimated parameters as initial guesses to replace in the LM algorithm and PWMs (lower and higher orders) produces the outcomes shown in Table 19 and in Figure 13 and Figure 14.

Table 17. Estimators and validation indices for the self-reported health dataset.

	Beta		Kumaraswamy		MBUW		Topp-Leone	Unit-Lindley
theta	$\mathit{\alpha }=8.0698$		$\mathit{\alpha }=5.7073$		$\mathit{\alpha }=0.5771$		7.3278	0.5929
	$\mathit{\beta }=3.8357$		$\mathit{\beta }=5.2250$		$\mathit{\beta }=1.2824$
Var	3.3978	1.6515	0.7936	1.2015			1.3768	0.0058
	1.6515	0.9425	1.2015	2.5192
SE	1.8433		0.8908				1.1734	0.0693
	0.9708		1.5872
AIC	-46.2684		-47.3865		-38.5536		-44.4298	-40.09617
CAIC	-45.9351		-47.0531		-38.2202		-44.3217	-39.9617
BIC	-42.9413		-44.0594		-35.2264		-42.7662	-38.4062
HQIC	-45.0746		-46.1927		-37.3598		-43.8329	-39.4729
LL	25.1342		25.6932		21.2768		23.2149	21.0349
K-S Value	0.0867		0.0684		0.1625		0.1234	0.1400
H0	Fail to reject		Fail to reject		Fail to reject		Fail to reject	Fail to reject
P-value	0.6674		0.8577		0.2284		0.5512	0.2102
AD	0.4674		0.3414		1.4935		0.8744	1.6862
CVM	0.0807		0.0547		0.2492		0.1397	0.2886

Table 17. Estimators and validation indices for the self-reported health dataset.

	Beta		Kumaraswamy		MBUW		Topp-Leone	Unit-Lindley
theta	$\mathit{\alpha }=8.0698$		$\mathit{\alpha }=5.7073$		$\mathit{\alpha }=0.5771$		7.3278	0.5929
	$\mathit{\beta }=3.8357$		$\mathit{\beta }=5.2250$		$\mathit{\beta }=1.2824$
Var	3.3978	1.6515	0.7936	1.2015			1.3768	0.0058
	1.6515	0.9425	1.2015	2.5192
SE	1.8433		0.8908				1.1734	0.0693
	0.9708		1.5872
AIC	-46.2684		-47.3865		-38.5536		-44.4298	-40.09617
CAIC	-45.9351		-47.0531		-38.2202		-44.3217	-39.9617
BIC	-42.9413		-44.0594		-35.2264		-42.7662	-38.4062
HQIC	-45.0746		-46.1927		-37.3598		-43.8329	-39.4729
LL	25.1342		25.6932		21.2768		23.2149	21.0349
K-S Value	0.0867		0.0684		0.1625		0.1234	0.1400
H0	Fail to reject		Fail to reject		Fail to reject		Fail to reject	Fail to reject
P-value	0.6674		0.8577		0.2284		0.5512	0.2102
AD	0.4674		0.3414		1.4935		0.8744	1.6862
CVM	0.0807		0.0547		0.2492		0.1397	0.2886

Table 17 shows that all the unit distributions fit the self-reported health data well. The same dilemma of the variance-covariance obtained after fitting the MBUW distribution is displayed, albeit the other statistical indices like AIC, CAIC, BIC, HQIC, LL are slightly less than that of the Beta, Kumaraswamy, and the Topp –Leone distributions. Figure 9 and Figure 10 show that the fitted CDF and the fitted PDF of the MBUW is lined up with the Beta and Kumaraswamy at the lower tail and is aligned with the Topp-Leone and the Unit Lindley distribution at the upper tail. Nevertheless, the misfortune is fundamentally the variance returned from implementing the MLE method executing Nelder Mead optimizer in MATLAB. Employing the values of the estimated parameters as initial guesses to insert in the LM algorithm and PWMs (lower and higher orders) generates the sequels shown in Table (20) and in Figure 15 and Figure 16.

Table 18. comparisons of the results of the PWMs method for parameter estimation of the water quality dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.3564			0.3529
	$\beta$	1.4307			1.4316
Var-cov matrix of parameter		10.7370		22.9418	627.8954	2.4064e+3
		22.9418		94.1037	2.4064e+3	9.3821e+3
		Eigenvalues=4.8406, 100			Eigenvalue=10.0408 , 1e+4
AD		0.3418 ( p=0.8810)			0.3820 (p=0.865)
CVM		0.0531 (p=0.8420)			0.0637 (p=0.785)
KS & p-value		0.0980 (p=0.7901)			0.1061 (p=0.7051
H0		Fail to reject			Fail to reject
SSE		1.8169e-19			1.7110e-19
${\widehat{\alpha }}_r$, unbiased estimator		0.3891			0.2491
${\widehat{\beta }}_s$, unbiased estimator		0.4441			0.3041
Sig.of $\alpha$parameter		1.648(p=0.1072)			0.5950(p=0.5552)
Sig.of $\beta$parameter		3.3898(p=0.0016)			1.2365(p=0.2235)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0181	0.0096		0.0214	0.0080
		0.0096	0.0051		0.0080	0.0030
		Eigenvalues=0.0232, 8.6736e-19 Determinant=0 Trace=0.0232			Eigenvalues=0.0244 , -4.3368e-19 Determinant=0 Trace=0.0244
Var-cov between M101&M110		8.8987e-4	0.0042
		0.0042	0.0020
Var-cov between M102 & M120					9.5798e-4	0.0197
					0.0197	0.002
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0468	0.0361		0.3518	0.2717
		0.0361	0.1781		0.2717	1.3403
		$\xi = 5.3356$to achieved condition number=5			$\xi = 0.7092$to achieved condition number=5
Jacobian matrix		-0.3796	0.0977		-0.2863	0.0735
		-1.4718	0.3786		-0.1067	0.0274

Table 18. comparisons of the results of the PWMs method for parameter estimation of the water quality dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.3564			0.3529
	$\beta$	1.4307			1.4316
Var-cov matrix of parameter		10.7370		22.9418	627.8954	2.4064e+3
		22.9418		94.1037	2.4064e+3	9.3821e+3
		Eigenvalues=4.8406, 100			Eigenvalue=10.0408 , 1e+4
AD		0.3418 ( p=0.8810)			0.3820 (p=0.865)
CVM		0.0531 (p=0.8420)			0.0637 (p=0.785)
KS & p-value		0.0980 (p=0.7901)			0.1061 (p=0.7051
H0		Fail to reject			Fail to reject
SSE		1.8169e-19			1.7110e-19
${\widehat{\alpha }}_r$, unbiased estimator		0.3891			0.2491
${\widehat{\beta }}_s$, unbiased estimator		0.4441			0.3041
Sig.of $\alpha$parameter		1.648(p=0.1072)			0.5950(p=0.5552)
Sig.of $\beta$parameter		3.3898(p=0.0016)			1.2365(p=0.2235)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0181	0.0096		0.0214	0.0080
		0.0096	0.0051		0.0080	0.0030
		Eigenvalues=0.0232, 8.6736e-19 Determinant=0 Trace=0.0232			Eigenvalues=0.0244 , -4.3368e-19 Determinant=0 Trace=0.0244
Var-cov between M101&M110		8.8987e-4	0.0042
		0.0042	0.0020
Var-cov between M102 & M120					9.5798e-4	0.0197
					0.0197	0.002
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0468	0.0361		0.3518	0.2717
		0.0361	0.1781		0.2717	1.3403
		$\xi = 5.3356$to achieved condition number=5			$\xi = 0.7092$to achieved condition number=5
Jacobian matrix		-0.3796	0.0977		-0.2863	0.0735
		-1.4718	0.3786		-0.1067	0.0274

Figure 11. Shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for water quality dataset.

Figure 11. Shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for water quality dataset.

Figure 12. Shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for water quality dataset.

Figure 12. Shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for water quality dataset.

Table 19. Comparisons of the results of the PWMs method for parameter estimation of the educational attainment dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.4310			0.4462
	$\beta$	1.3483			1.3442
Var-cov matrix of parameter		679.8845		2.5075e+3	679.7206	2.4966e+3
		2.5075e+3		9.3254e+3	2.4966e+3	9.3313e+3
		Eigenvalues=5.2672, 10000			Eigenvalue=10.9762 , 1e+4
AD		1.4027 ( p=0.2090)			1.6993( p=0.1330)
CVM		0.2420 (p=0.1990)			0.3333(p=0.1110)
KS & p-value		0.1787 (p=0.1371)			0.2041(p=0.0615)
H0		Fail to reject			Fail to reject
SSE		2.2362e-19			2.4395e-19
${\widehat{\alpha }}_r$, unbiased estimator		0.3485			0.2139
${\widehat{\beta }}_s$, unbiased estimator		0.4325			0.2979
Sig.of $\alpha$parameter		2.0345(p=0.0487)			0.6592(p=0.5136)
Sig.of $\beta$parameter		3.2983(p=0.0021)			1.0279(p=0.3104)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0190	0.0107		0.0213	0.0089
		0.0107	0.0060		0.0089	0.0037
		Eigenvalues=0.0250, 2.6021e-18 Determinant=0 Trace=0.0250			Eigenvalues=0.0250 , 4.3368e-19 Determinant=0 Trace=0.0250
Var-cov between M101&M110		6.7399e-4	0.0038
		0.0038	0.0020
Var-cov between M102 & M120					6.7302e-4	0.0264
					0.0264	0.0020
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0449	0.0355		0.4582	0.3613
		0.0355	0.1671		0.3613	1.7103
		$\xi = 5.6607$to achieve a condition number=5			$\xi = 0.5534$to achieve a condition number=5
Jacobian matrix		-0.3668	0.0987		-0.2691	0.0721
		-0.2059	0.0554		-0.1119	0.0300

Table 19. Comparisons of the results of the PWMs method for parameter estimation of the educational attainment dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.4310			0.4462
	$\beta$	1.3483			1.3442
Var-cov matrix of parameter		679.8845		2.5075e+3	679.7206	2.4966e+3
		2.5075e+3		9.3254e+3	2.4966e+3	9.3313e+3
		Eigenvalues=5.2672, 10000			Eigenvalue=10.9762 , 1e+4
AD		1.4027 ( p=0.2090)			1.6993( p=0.1330)
CVM		0.2420 (p=0.1990)			0.3333(p=0.1110)
KS & p-value		0.1787 (p=0.1371)			0.2041(p=0.0615)
H0		Fail to reject			Fail to reject
SSE		2.2362e-19			2.4395e-19
${\widehat{\alpha }}_r$, unbiased estimator		0.3485			0.2139
${\widehat{\beta }}_s$, unbiased estimator		0.4325			0.2979
Sig.of $\alpha$parameter		2.0345(p=0.0487)			0.6592(p=0.5136)
Sig.of $\beta$parameter		3.2983(p=0.0021)			1.0279(p=0.3104)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0190	0.0107		0.0213	0.0089
		0.0107	0.0060		0.0089	0.0037
		Eigenvalues=0.0250, 2.6021e-18 Determinant=0 Trace=0.0250			Eigenvalues=0.0250 , 4.3368e-19 Determinant=0 Trace=0.0250
Var-cov between M101&M110		6.7399e-4	0.0038
		0.0038	0.0020
Var-cov between M102 & M120					6.7302e-4	0.0264
					0.0264	0.0020
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0449	0.0355		0.4582	0.3613
		0.0355	0.1671		0.3613	1.7103
		$\xi = 5.6607$to achieve a condition number=5			$\xi = 0.5534$to achieve a condition number=5
Jacobian matrix		-0.3668	0.0987		-0.2691	0.0721
		-0.2059	0.0554		-0.1119	0.0300

Figure 13. shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for educational attainment dataset.

Figure 13. shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for educational attainment dataset.

Figure 14. shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for educational attainment dataset.

Figure 14. shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for educational attainment dataset.

Table 20. comparisons of the results of the PWMs method for parameter estimation of the self-reported health dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.5843			0.5752
	$\beta$	1.2806			1.2829
Var-cov matrix of parameter		11.610		21.6717	71.1767	230.3148
		21.6717		94.6865	230.3148	942.8902
		Eigenvalues=6.2964 , 100			Eigenvalue=14.0669 , 1000
AD		1.4443 ( p=0.2010)			1.5100 ( p=0.1780)
CVM		0.2382 (p=0.2100)			0.2529 (p=0.1900)
KS & p-value		0.1549 (p=0.2768)			0.1645 (p=0.2167)
H0		Fail to reject			Fail to reject
SSE		5.9631e-19			0
${\widehat{\alpha }}_r$, unbiased estimator		0.3019			0.1865
${\widehat{\beta }}_s$, unbiased estimator		0.3775			0.2621
Sig.of $\alpha$parameter		2.7038(p=0.0102)			0.7048 (p=0.4852)
Sig.of $\beta$parameter		3.0045(p=0.0047)			0.7988 (p=0.4293)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0172	0.0108		0.0205	0.0099
		0.0108	0.0068		0.0099	0.0047
		Eigenvalues=0.0240, 1.7347e-18 Determinant=0 Trace=0.0240			Eigenvalues=0.0253 , 0 Determinant=0 Trace=0.0253
Var-cov between M101& M110		4.0448e-4	0.0029
		0.0029	0.0020
Var-cov between M102 & M120					4.2498e-4	0.0347
					0.0347	0.0020
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0467	0.0352		0.6661	0.5053
		0.0352	0.1817		0.5053	2.5787
		$\xi = 5.2546$to achieve a condition number=5			$\xi = 0.3698$to achieve a condition number=5
Jacobian matrix		-0.3062	0.0717		-0.2317	0.0574
		-1.5296	0.3579		-0.1112	0.0276

Table 20. comparisons of the results of the PWMs method for parameter estimation of the self-reported health dataset using $M_{\mathrm{1,0},1} \& M_{\mathrm{1,1},0}$

		Using unbiased Sample estimator for $M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0}$			Using unbiased Sample estimator for $M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0}$
thetas	$\alpha$	0.5843			0.5752
	$\beta$	1.2806			1.2829
Var-cov matrix of parameter		11.610		21.6717	71.1767	230.3148
		21.6717		94.6865	230.3148	942.8902
		Eigenvalues=6.2964 , 100			Eigenvalue=14.0669 , 1000
AD		1.4443 ( p=0.2010)			1.5100 ( p=0.1780)
CVM		0.2382 (p=0.2100)			0.2529 (p=0.1900)
KS & p-value		0.1549 (p=0.2768)			0.1645 (p=0.2167)
H0		Fail to reject			Fail to reject
SSE		5.9631e-19			0
${\widehat{\alpha }}_r$, unbiased estimator		0.3019			0.1865
${\widehat{\beta }}_s$, unbiased estimator		0.3775			0.2621
Sig.of $\alpha$parameter		2.7038(p=0.0102)			0.7048 (p=0.4852)
Sig.of $\beta$parameter		3.0045(p=0.0047)			0.7988 (p=0.4293)
Variance of the function of the parmaeter, after the delta method application. Determinant and trace of this matrix		0.0172	0.0108		0.0205	0.0099
		0.0108	0.0068		0.0099	0.0047
		Eigenvalues=0.0240, 1.7347e-18 Determinant=0 Trace=0.0240			Eigenvalues=0.0253 , 0 Determinant=0 Trace=0.0253
Var-cov between M101& M110		4.0448e-4	0.0029
		0.0029	0.0020
Var-cov between M102 & M120					4.2498e-4	0.0347
					0.0347	0.0020
$var\left(\widehat{\theta }\right)$and associated $\xi$		0.0467	0.0352		0.6661	0.5053
		0.0352	0.1817		0.5053	2.5787
		$\xi = 5.2546$to achieve a condition number=5			$\xi = 0.3698$to achieve a condition number=5
Jacobian matrix		-0.3062	0.0717		-0.2317	0.0574
		-1.5296	0.3579		-0.1112	0.0276

Figure 15. shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for self-reported health dataset.

Figure 15. shows the eCDF vs. fitted CDF for MLE & low order PWMs, M101, M110 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for self-reported health dataset.

Figure 16. shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for self-reported health dataset.

Figure 16. shows the eCDF vs. fitted CDF for MLE & high order PWMs, M102, M120 (upper left subplot), histogram & fitted PDFs ( upper right subplot), QQ plot (lower left) and PP plot (lower right) for self-reported health dataset.

Table 18 shows the results of the water quality dataset. The values of the estimated parameters using the lower PMWs and higher PWMs are nearly equal. The approximated variance-covariance obtained from the LM algorithm using the lower order PWMs is less than that obtained from using the higher order PWMs. Applying the delta method on both lower and higher PWMs gives variances comparable to each other. The determinant for the matrix of the lower order (0.0232) is less than that of the higher order (0.0244). The Jacobian matrix is an ill condition matrix for both PWMs. Using a regularization factor, which is lesser for the higher order PWMs, enhances the estimated variance of the estimated parameters. This regularized factor is chosen so the condition number of this matrix $\left({\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]\right)$ is 5. This is because the estimated $var\left(\widehat{\theta }\right)={\left[{\left[g^{\text{'}}\left({\theta }_0\right)\right]}^T {\left(var\left[g\left(\widehat{\theta }\right)\right]\right)}^{-1}\left[g^{\text{'}}\left({\theta }_0\right)\right]+ \xi I\right]}^{-1}$ requires this factor $\left(\xi \right)$. The standard errors of the estimated parameters are obtained by taking the square root of the diagonal of this last matrix. The estimated variance; $var\left(\widehat{\theta }\right) ;$ obtained from the lower order PWMs is less than that from the higher PWMS. The covariance between the samples PWMs is also reported; $var\left[g\left(\widehat{\theta }\right)\right]=\left[\begin{array}{cc}var\left(M_{101}\right)& cov\left(M_{101},M_{110}\right)\\ cov\left(M_{101},M_{110}\right)& var\left(M_{110}\right)\end{array}\right]$ and it is less for the lower order PWMs than the covariance for the higher order PWMs. The AD, CVM, and KS tests are statistically significant for both order PWMs. The estimated $\left(\widehat{\beta }\right)$ obtained from lower PWMs is statistically significant. Figure 11 and Figure 12 show the fitted CDFs and the fitted PDFs for the lower and higher order PWMs respectively. The curves of the fitted CDF and the fitted PDF obtained from MLE are perfectly aligned with the curves obtained from using both orders of PWMs. The QQ plot and PP plot are almost perfectly aligned with the diagonal.

Table 19 clarifies the same type of the results for the educational attainment dataset. The values of the estimated parameters are nearly alike. The approximated variance yielded from the LM algorithm is nearly comparable between the two orders of the PWMs. After applying the delta method, the determinant of the covariance matrix obtained from both orders of the PWMs is identical (0.025). The variance of the estimated parameter, $var\left(\widehat{\theta }\right)$ , is less for lower order PWMs than that for the higher order PWMs. The AD, CVM, KS tests are statistically significant. The estimated $\left(\widehat{\alpha }\right) \&\left(\widehat{\beta }\right)$ obtained from lower PWMs are statistically significant. Figure 13 and Figure 14 disclose the fitted CDFs and the fitted PDFs for the lower and higher order PWMs respectively. The fitted CDF and the fitted PDF exhibit curves gained from MLE that are seamlessly aligned with those curves attained from using both orders of PWMs. The QQ plot and PP plot display fair alignment with the diagonal at the upper tail and at the lower tail respectively.

Table 20 expresses similar type of the outcomes for the self-reported health dataset. The values of the estimated parameters are nearly indistinguishable. The approximated variance gained from the LM algorithm is nearly comparable between the two orders of the PWMs. After applying the delta method, the determinant of the covariance matrix obtained from lower order PWMs (0.0240) is less than that produced from the higher order PWMs (0.0253). The variance of the estimated parameter, $var\left(\widehat{\theta }\right)$ , is less for lower order PWMs than that for the higher order PWMs. The AD, CVM, KS tests are statistically significant. The estimated $\left(\widehat{\alpha }\right) \&\left(\widehat{\beta }\right)$obtained from lower PWMs are statistically significant. Figure 15 and Figure 16 unveil the fitted CDFs and the fitted PDFs for the lower and higher order PWMs respectively. The curves of the fitted CDF and the fitted PDF obtained from MLE are perfectly aligned with the curves obtained from using both orders of PWMs. The QQ plot and PP plot illustrate reasonable alignment with the diagonal.

6. Conclusions

The classic (PWMs) method can be effectively utilized for the parameter estimation of the new Median Based Unit Weibull (MBUW) distribution. This method is robust against outliers and is more straightforward to obtain than the maximum likelihood estimator. In the context of the fitting datasets employed in this study, the unbiased estimators for both order of PWMs yielded comparable results. The PWMs method offers several advantages over alternative estimation techniques. It is both rapid and easy to compute, consistently producing feasible values for the estimated parameters and the estimated variances. Although the higher order PWMs take much longer time to execute at lager sample sizes than MLE, they are more efficient in estimating $\left(\widehat{\beta }\right).$ Furthermore, PWMs estimators exhibit asymptotic normal distributions. It was observed that higher-order moments like $(M_{\mathrm{1,0},2} and M_{\mathrm{1,2},0})$ did not significantly contribute more information than the commonly used moments $(M_{\mathrm{1,0},1} and M_{\mathrm{1,1},0})$. The relative efficiency of the lower order PWMs to the higher order PWMs is in the range from 0.8 to 1 depending on the values of the parameters. The Average Absolute Biases for both orders are more or less comparable and nearly equal so there is no big difference between them. The computational time is less for lower order PWMs than the time for the higher PWMs. The estimated parameters from both orders of the PWMs are asymptotically consistent. The delta method used to estimate the variance of the function of the parameters is asymptotically normal. Using regularization factor to estimate the variance of the estimated parameters, after accounting for the covariance between the sample PWMs, can be a solution to mitigate the dependency between the parameters. The main root cause of this dependency between the parameters is the distribution itself. The method may configure the shape of the dependency, being linear or non-linear and monotonic or non-monotonic.

In this paper, the formula employed for defining PWMs is primarily based on the binomial expansion of the cumulative distribution function and the survival function, thus integrating with respect to the variable (dy) rather than integrating with respect to the cumulative distribution function (dF). The quantile function of the MBUW is difficult to integrate in a sense to obtain a system of equations that can be solved numerically.

7. Future Work

PWMs serve as a powerful foundation for L-moments, paving the way for a wealth of analytical possibilities. In future research, we can leverage the estimation of various L-moments, including L-skewness and L-kurtosis, alongside the L-moment method for precise parameter estimation. Furthermore, by extending the $\left(\beta \right)$parameter of the Median Based unit Weibull (MBUW) distribution to accommodate negative values, we can broaden our analytical horizons. In such scenarios, the Generalized Probability-Weighted Moments (GPWMs) method stands out as an essential tool. Additionally, we can effectively apply Partial GPWMs to handle censored data, ensuring a comprehensive approach to our analyses. The Bayesian methods with different loss functions can be consumed in the future to estimate the parameters. The correlation between the two parameters is a proposal for future studies to alleviate this dependency and hence to enhance the estimation and inference procedures. The tail indices can be further assessed to see if they can benefit from higher order PWMs. These directions not only enhance our statistical methodologies but also significantly enrich our understanding of the underlying data distributions.