A competitive family to the Beta and Kumaraswamy generators: Properties, Regressions and Applications

Abstract

We define two new flexible families of continuous distributions to fit real data by compoun-ding the Marshall–Olkin class and the power series distribution. These families are very competitive to the popular beta and Kumaraswamy generators. Their densities have linear representations of exponentiated densities. In fact, as the main properties of thirty five exponentiated distributions are well-known, we can easily obtain several properties of about three hundred fifty distributions using the references of this article and five special cases of the power series distribution. We provide a package implemented in R software that shows numerically the precision of one of the linear representations. This package is useful to calculate numerical values for some statistical measurements of the generated distributions. We estimate the parameters by maximum likelihood. We define a regression based on one of the two families. The usefulness of a generated distribution and the associated regression is proved empirically.

Key words
generating function; Marshall–Olkin family; maximum likelihood; moment; distribution

Introduction

The Marshall–Olkin (“MO“) family (Marshall & Olkin 1997) adds one parameter to a parent distribution. Let $G\left(z\right)=G(z;𝛕)$ be the parent cumulative distribution function (cdf) of a random variable $Z$ with parameter vector $𝛕={({\tau }_1,\dots ,{\tau }_q)}^{\top }$. The survival function and probability density function (pdf) of $Z$ are $\bar{G}\left(z\right)=\bar{G}(z;𝛕)$ and $g\left(z\right)=g(z;𝛕)$, respectively.

The cdf $H\left(z\right)$ and survival function $\bar{H}\left(z\right)=1-H\left(z\right)$ of the MO class with baseline $G(z;𝛕)$ are

$$H\left(z\right)=H(z;\alpha ,𝛕)=\frac{G(z;𝛕)}{1-\bar{\alpha }\bar{G}(z;𝛕)},z\in ℝ,\alpha >0,$$(1)

and

$$\bar{H}\left(z\right)=\bar{H}(z;\alpha ,𝛕)=\frac{\alpha \bar{G}(z;𝛕)}{1-\bar{\alpha }\bar{G}(z;𝛕)},$$(2)

respectively, where $\bar{\alpha }=1-\alpha$.

Equation (1) can generate many continuous distributions from popular ones. The MO-G density function can be expressed as

$$h\left(z\right)=h(z;\alpha ,𝛕)=\frac{\alpha g(z;𝛕)}{{[1-\bar{\alpha }\bar{G}(z;𝛕)]}^2}.$$(3)

For $\alpha =1$, $h\left(z\right)=g(z;𝛕)$ is the simplest case of (3). Marshall & Olkin 1997 pioneered the MO-Weibull (MOW) distribution which is a useful extension of the Weibull.

Consider $N$ random variables $Z_1,\dots ,Z_N$ independent and identically distributed (i.i.d.) with cdf $H\left(z\right)$ and pdf $h\left(z\right)$ given by (1) and (3), respectively. Here, $N$ is a discrete random variable with support $\{1,2,\dots \}$. Henceforth, let $X=max\{Z_1,\dots ,Z_N\}$ and $Y=min\{Z_1,\dots ,Z_N\}$ be two random variables assuming that $N$ has the zero-truncated power series (PS) distribution with probability mass function (pmf)

$$p_n=\mathbb{P}(N=n;\theta )=\frac{a_n{\theta }^n}{C\left(\theta \right)},n=1,2,\dots ,$$(4)

where $a_n>0$ (for $n\ge 1$), $\theta$ is called the power parameter and $C\left(\theta \right)={\sum }_{n=1}^{\infty }{a}_n{\theta }^n>0$. The probability generating function (pgf) of $N$ is $P\left(z\right)=E\left(z^N\right)=C\left(z\theta \right)/C\left(\theta \right)$.

Five important distributions are special cases of (4): the zero-truncated Poisson, logarithmic, negative binomial, geometric and zero-truncated binomial distributions.

The cdf of $X=max\{Z_1,\dots ,Z_N\}$ conditional given $N=n$ is

$$\begin{array}{c}\hfill F_X(x\mid N=n)=\mathbb{P}[X\le x|N=n]=H{(x;\alpha ,𝛕)}^n,\end{array}$$

and then the unconditional cdf of $X$ follows from (4)

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }H{(x;\alpha ,𝛕)}^n\frac{a_n{\theta }^n}{C\left(\theta \right)}=\frac{C\left(\theta H(x;\alpha ,𝛕)\right)}{C\left(\theta \right)}.\end{array}$$(5)

The conditional cdf of $Y=min\{Z_1,\dots ,Z_N\}$ under $N=n$ is

$$\begin{array}{c}\hfill F_Y(y\mid N=n)=\mathbb{P}[Y\le y|N=n]=1-\bar{H}{(y;\alpha ,𝛕)}^n,\end{array}$$

and then the unconditional cdf of $Y$ follows from (4) as

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\sum _{n=1}^{\infty }\bar{H}{(y;\alpha ,𝛕)}^n\frac{a_n{\theta }^n}{C\left(\theta \right)}=1-\frac{C\left(\theta \bar{H}(y;\alpha ,𝛕)\right)}{C\left(\theta \right)}.\end{array}$$(6)

Equations (5) and (6) define two Marshall–Olkin Power Series-G (MOPS-G) families under baseline G. They provide a strong motivation for explaining the failure time of any mechanism formed by an unknown number $N$ of identical and independent (parallel or serial) components. The densities of $X$ and $Y$ are obtained by differentiating (5) and (6). We emphasize that these equations can generate many MOPS models. For each baseline G, we can generate ten ($2\times 5$) associated models from the five discrete distributions in Equation (4). For $\alpha =1$, we have the Power Series-G (PS-G) classes under baseline G.

The minimum ($Y$) and maximum ($X$) statistics can be applied in several series and parallel systems with identical components and have many industrial and biological applications. In parallel systems, the random variable $Y$ models the time of the first component to fail, while $X$ models the time for the breakout system. A dual interpretation can be given for systems with serial components. These random variables are also very useful in oncology. For example, suppose we are studying a recurrence of a certain type of cancerous tumor of an individual after undergoing any kind of treatment. So, the time for the first cell to activate to produce cancer cells can be modeled by the generated distribution of $Y$, while the disease manifestation (if it occurs only after an unknown number of factors have been active) can be modeled by the generated distribution of $X$.

Four new distributions based on the MOPS construction are introduced for illustrative purposes in Section Four special models. We derive linear representations for the densities of $X$ and $Y$ in Section Expansions. A package in R is presented in Section Numerical evaluation to calculate numerically several mathematical properties for the generated distributions based on the linear representations. General structural properties for the two families are addressed in Section Properties. In Section Estimation, we estimate the parameters for one of the families. We introduce in Section Regression the Marshall–Olkin Truncated Poisson Weibull regression defined from one of the families. In Section Two simulation studies, some simulations examine the accuracy of the maximum likelihood estimates (MLEs) and the quantile residuals (qrs). Two applications prove the utility of our finding in Section Applications. Finally, we offer concluding remarks in Section Conclusions.

Four special models

First, consider the zero-truncated Poisson in (4). The cdfs of the Marshall–Olkin Zero-Truncated Poisson-G (MOTP-G) distributions are determined from Equations (5) and (6) as

$$\begin{array}{c}\hfill F_X\left(x\right)={({\mathrm{\text{e}}}^{\theta }-1)}^{-1}[exp\{\theta H(x;\alpha ,𝛕)\}-1]\end{array}$$(7)

and

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-{({\mathrm{\text{e}}}^{\theta }-1)}^{-1}[exp\left\{\theta \bar{H}(y;\alpha ,𝛕)\right\}-1].\end{array}$$(8)

The Weibull cdf with scale parameter $\lambda >0$ and shape parameter $\beta >0$ is (for $x\ge 0$)

$$\begin{array}{c}\hfill G(z;\lambda ,\beta )=1-exp[-{\left(\lambda z\right)}^{\beta }].\end{array}$$

Then, the cdf and survival function of the MO-Weibull (MOW) distribution are

$$\begin{array}{c}\hfill H\left(z\right)=H(z;\alpha ,\lambda ,\beta )=\frac{1-exp[-{\left(\lambda z\right)}^{\beta }]}{1-\bar{\alpha }exp[-{\left(\lambda z\right)}^{\beta }]}\end{array}$$(9)

and

$$\begin{array}{c}\hfill \bar{H}\left(z\right)=\bar{H}(z;\alpha ,\lambda ,\beta )=\frac{\alpha exp[-{\left(\lambda z\right)}^{\beta }]}{1-\bar{\alpha }exp[-{\left(\lambda z\right)}^{\beta }]},\end{array}$$(10)

respectively.

By inserting the last two formulae in Equations (7) and (8) and differentiating the resulting expressions, we obtain the MOTP-Weibull (MOTPW) densities

\begin{aligned} \label{MOTPW1} f_X(x)=\frac{\alpha\,\theta\,\beta\,\lambda^{\beta}\,x^{\beta-1}\,\rm{e}^{-u}}{(\rm{e}^\theta-1)\,(1-\bar{\alpha}\rm{e}^{-u})^2}\, \exp\Bigg[\frac{\theta\,(1-\rm{e}^{-u})}{(1-\bar{\alpha}\rm{e}^{-u})}\Bigg]\end{aligned}(11)

and

\begin{aligned} \label{MOTPW2} f_Y(y)=\frac{\alpha\,\theta\,\beta\,\lambda^{\beta}\,x^{\beta-1}\,\rm{e}^{-u}}{(\rm{e}^\theta-1)\,(1-\bar{\alpha}\rm{e}^{-u})^2}\, \exp\Bigg[\frac{\alpha\,\theta\,\rm{e}^{-u}}{(1-\bar{\alpha}\rm{e}^{-u})}\Bigg],\end{aligned}(12)

respectively, where $u=u\left(x\right)={\left(\lambda x\right)}^{\beta }$ in $f_X\left(x\right)$ and $u=u\left(y\right)={\left(\lambda y\right)}^{\beta }$ in $f_Y\left(y\right)$.

Second, consider the geometric distribution in (4). The cdfs of the Marshall–Olkin Geometric-G (MOG-G) classes follow from Equations (5) and (6)

$$\begin{array}{c}\hfill F_X\left(x\right)=\frac{(1-\theta )}{\theta }[\frac{\theta H(x;\alpha ,𝛕)}{1-\theta H(x;\alpha ,𝛕)}]\end{array}$$(13)

and

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\frac{(1-\theta )}{\theta }[\frac{\theta \bar{H}(y;\alpha ,𝛕)}{1-\theta \bar{H}(y;\alpha ,𝛕)}].\end{array}$$(14)

The Burr XII (BXII) cdf is (for $x>0$)

$$\begin{array}{c}\hfill G(z;\beta ,\lambda )=1-{(1+z^{\beta })}^{-\lambda },\end{array}$$(15)

where $\beta >0$ and $\lambda >0$ are shape parameters. For $\lambda =1$ and $\beta =1$ in Equation (15), we have the log-logistic (LL) and Lomax distributions, respectively.

Hence, the cdf and survival function of the Marshall–Olkin Burr XII (MOBXII) distribution are

$$\begin{array}{c}\hfill H\left(z\right)=H(z;\alpha ,\lambda ,\beta )=\frac{1-{(1+z^{\beta })}^{-\lambda }}{1-\bar{\alpha }{(1+z^{\beta })}^{-\lambda }}\end{array}$$(16)

and

$$\begin{array}{c}\hfill \bar{H}\left(z\right)=\bar{H}(z;\alpha ,\lambda ,\beta )=\frac{\alpha {(1+z^{\beta })}^{-\lambda }}{1-\bar{\alpha }{(1+z^{\beta })}^{-\lambda }},\end{array}$$(17)

respectively.

By inserting the last two formulae in Equations (13) and (14) and differentiating the resulting expressions with respect to $x$ and $y$, respectively, we obtain the MOG-Burr XII (MOGBXII) densities

$$\begin{array}{c}\hfill f_X\left(x\right)=\frac{\alpha \beta \lambda (1-\theta )x^{\beta -1}{(1+x^{\beta })}^{-\lambda -1}}{{[1-\theta -(1-\alpha -\theta ){(1+x^{\beta })}^{-\lambda }]}^2}\end{array}$$(18)

and

$$\begin{array}{c}\hfill f_Y\left(y\right)=\frac{\alpha \beta \lambda (1-\theta )x^{\beta -1}{(1+x^{\beta })}^{-\lambda -1}}{{\{1-[1-(1-\theta )\alpha ]{(1+x^{\beta })}^{-\lambda }\}}^2}.\end{array}$$(19)

For the MOTPW and MOGBXII distributions (to the maximum $X$) referred to (11) and (18), some plots of the densities and cumulative functions are displayed in Figures 1 and 2, respectively. The various forms of the densities indicate more flexibility than the parent distributions.

Figure 1
Plots of the density and cumulative functions of the MOTPW distribution under four scenarios. (a) $𝛂=30$, $𝛌=2$, $𝛃=1.5$, and varying $𝛉$. (b) $𝛂=30$, $𝛌=2$, $𝛃=1.5$, and varying $𝛉$. (c) $𝛉=0.09$, $𝛌=2$, $𝛃=1.5$, and varying $𝛂$. (d) $𝛉=0.09$, $𝛌=2$, $𝛃=1.5$, and varying $𝛂$.

Figure 2
Plots of the density and cumulative functions of the MOGBXII distribution under four scenarios. (a) $𝛂=10$, $𝛌=2$, $𝛃=1.5$, and varying $𝛉$. (b) $𝛂=10$, $𝛌=2$, $𝛃=1.5$, and varying $𝛉$. (c) $𝛉=0.9$, $𝛌=2$, $𝛃=1.5$, and varying $𝛂$. (d) $𝛉=0.9$, $𝛌=2$, $𝛃=1.5$, and varying $𝛂$.

We can note increasing, decreasing, and unimodal shapes for the hrf of the MOTPW distribution in Figure 3. Also, we see a slightly different hrf with increasing, decreasing and increasing shape.

Figure 3
Plots of the hrf of the MOTPW model.

Graphics comparing the histograms from two simulated data sets and the MOTPW and MOGBXII densities of $X$ under specified parameters are reported in Figure 4. They show good agreement between the simulated values and these densities.

Figure 4
Plots of the MOTPW (a) and MOGBXII (b) densities and histograms of simulated data.

Expansions

We obtain useful linear representations for the density functions of $X$ and $Y$ for two separated cases $\alpha \in (0,1)$ and $\alpha >1$. For $\alpha =1$, we have $H(z;1,𝛕)=G(z;𝛕)$.

By inserting (1) in Equation (5) and letting $\bar{G}\left(x\right)=\bar{G}(x;𝛕)$, we can write

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }\frac{p_{n}G{\left(x\right)}^n}{{[1-\bar{\alpha }\bar{G}\left(x\right)]}^n}.\end{array}$$(20)

First, we consider the density of the maximum $X$ when $\alpha \in (0,1)$. For $\left|z\right|<1$ and $n=1,2,\dots$, the negative binomial expansion holds

$${(1-z)}^{-n}=\sum _{k=0}^{\infty }\left(\genfrac{}{}{0ex}{}{-n}{k}\right){(-z)}^k.$$(21)

Expanding ${[1-\bar{\alpha }\bar{G}\left(z\right)]}^{-n}$ as in Equation (21) since $\alpha \in (0,1)$, we have

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }\left(\genfrac{}{}{0ex}{}{-n}{k}\right){(-\bar{\alpha })}^k{p}_{n}G{\left(x\right)}^n{[1-G\left(x\right)]}^k.\end{array}$$

Henceforth, let $T_s\sim$exp-G$\left(s\right)$ be the exponentiated-G (exp-G) random variable with power parameter $s>0$. Its cdf and pdf are ${\Pi }_s\left(x\right)={\Pi }_s(x;𝛕)=G{(x;𝛕)}^s$ and ${\pi }_s\left(x\right)={\pi }_s(x;𝛕)=sG{(x;𝛕)}^{s-1}g(x;𝛕)$, respectively. Many exp-G properties have been studied exhaustively by several authors (Tahir & Nadarajah 2015). We can write

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }{w}_{n,0}{\Pi }_n\left(x\right)+\sum _{n=1}^{\infty }\sum _{k=1}^{\infty }{w}_{n,k}{\Pi }_n\left(x\right){[1-G\left(x\right)]}^k,\end{array}$$

where $w_{n,k}=w_{n,k}(\alpha ,\theta )=\left(\genfrac{}{}{0ex}{}{-n}{k}\right){(-\bar{\alpha })}^k{p}_n$ for $n=1,2,\dots$ and $k=0,1,\dots$ Further, using the binomial theorem, we obtain

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }{w}_{n,0}{\Pi }_n\left(x\right)+\sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\sum _{i=0}^k{w}_{n,k,i}{\Pi }_{n+i}\left(x\right),\end{array}$$

where $w_{n,k,i}={(-1)}^i\left(\genfrac{}{}{0ex}{}{k}{i}\right)w_{n,k}$ for $i=0,1,\dots ,k$.

By differentiating the last equation, we obtain

$$\begin{array}{c}\hfill f_X\left(x\right)=\sum _{n=1}^{\infty }{w}_{n,0}{\pi }_n\left(x\right)+\sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\sum _{i=0}^k{w}_{n,k,i}{\pi }_{n+i}\left(x\right).\end{array}$$(22)

We now move to the density of the maximum $X$ when $\alpha >1$. We modify the denominator in (20)

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }\frac{p_{n}G{\left(x\right)}^n}{{\alpha }^n{[1-(1-{\alpha }^{-1})G\left(x\right)]}^n}\end{array}$$

and then apply Equation (21) to find

$$\begin{array}{c}\hfill F_X\left(x\right)=\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }{v}_{n,k}{\Pi }_{n+k}\left(x\right),\end{array}$$

where $v_{n,k}=v_{n,k}(\alpha ,\theta )={(-1)}^k\left(\genfrac{}{}{0ex}{}{-n}{k}\right){\alpha }^{-n}{(1-{\alpha }^{-1})}^k{p}_n$ (for $n=1,2,\dots$ and $k=0,1,\dots$). By differentiating $F_X\left(x\right)$, the density of $X$ follows as

$$\begin{array}{c}\hfill f_X\left(x\right)=\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }{v}_{n,k}{\pi }_{n+k}\left(x\right).\end{array}$$(23)

Next, we consider the density of the minimum $Y$. By inserting (2) in Equation (6), we have

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\sum _{n=1}^{\infty }\frac{{\alpha }^n{p}_n\bar{G}{\left(y\right)}^n}{{[1-\bar{\alpha }\bar{G}\left(y\right)]}^n}.\end{array}$$(24)

For $\alpha \in (0,1)$, we apply expansion (21) in the last equation to

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }{q}_{n,k}\bar{G}{\left(y\right)}^{n+k},\end{array}$$

where $q_{n,k}=q_{n,k}(\alpha ,\theta )={(-1)}^k\left(\genfrac{}{}{0ex}{}{-n}{k}\right){\bar{\alpha }}^k{\alpha }^n{p}_n$ for $n=1,2,\dots$ and $k=0,1,\dots$

By using the binomial theorem in $\bar{G}{\left(y\right)}^{n+k}$, we have

$$\begin{array}{c}\hfill F_Y\left(y\right)=1+\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }\sum _{i=0}^{n+k}{q}_{n,k,i}{\Pi }_i\left(y\right),\end{array}$$

where $q_{n,k,i}={(-1)}^{i+1}\left(\genfrac{}{}{0ex}{}{n+k}{i}\right)q_{n,k}$ for $i=0,1,\dots ,n+k$.

By differentiating $F_Y\left(y\right)$, the density of $Y$ can be expressed as

$$\begin{array}{c}\hfill f_Y\left(y\right)=\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }\sum _{i=1}^{n+k}{q}_{n,k,i}{\pi }_i\left(y\right).\end{array}$$(25)

We now obtain the density of $Y$ when $\alpha >1$. By changing the denominator in Equation (24), we have

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\sum _{n=1}^{\infty }\frac{p_n\bar{G}{\left(y\right)}^n}{{[1-(1-{\alpha }^{-1})G\left(y\right)]}^n}.\end{array}$$

Applying expansion (21) in the last equation

$$\begin{array}{c}\hfill F_Y\left(y\right)=1-\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }{t}_{n,k}\bar{G}{\left(y\right)}^{n}G{\left(y\right)}^k,\end{array}$$

where $t_{n,k}=t_{n,k}(\alpha ,\theta )={(-1)}^k{(1-{\alpha }^{-1})}^k\left(\genfrac{}{}{0ex}{}{-n}{k}\right)p_n$ (for $n=1,2,\dots$ and $k=0,1,\dots$).

Using the binomial theorem, we can rewrite $F_Y\left(y\right)$ as

$$\begin{array}{c}\hfill F_Y\left(y\right)=1+\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }\sum _{i=0}^n{t}_{n,k,i}{\Pi }_{i+k}\left(y\right),\end{array}$$

where $t_{n,k,i}={(-1)}^{i+1}\left(\genfrac{}{}{0ex}{}{n}{i}\right)t_{n,k}$ for $i=0,1,\dots$ By simple differentiation

$$\begin{array}{c}\hfill f_Y\left(y\right)=\sum _{n=1}^{\infty }\sum _{k=0}^{\infty }\sum _{i=0}^n{t}_{n,k,i}{\pi }_{i+k\ge 1}\left(y\right),\end{array}$$(26)

where ${\pi }_{i+k\ge 1}\left(y\right)$ is the exp-G density with power parameter $i+k\ge 1$.

Equations (22), (23), (25) and (26) are the main results of this section. These linear representations have great utility for deriving structural properties of the maximum $X$ and minimum $Y$ from well-known exp-G properties. More than thirty five exp-G models have been studied so far and then it is possible to construct at least three hundred fifty ($70\times 5$) MOPS-G models with properties determined from those exp-G properties. We can use statistical platforms with ten terms to have precise results.

Numerical evaluation

In order to evaluate the analytical results presented in the previous sections, a package was implemented using the R programming language (R Core Team 2022). The MarshallOlkinPSG package was constructed in a generic way, that is, its most important functions allow generalizations for any baseline G distribution or even inform a zero-truncated PS distribution.

The library code can be obtained from GitHub at https://github.com/prdm0/MarshallOlkinPSG. On the library’s website (see https://prdm0.github.io/MarshallOlkinPSG) it is possible to have more information on the functions implemented through the documentation and usage examples.

To install the package hosted and maintained on GitHub, it is necessary to previously install the remotes library. With the prerequisite met, the package MarshallOlkinPSG can be installed as:

# Install the remotes package: # install.packages("remotes") remotes::install_github("prdm0/MarshallOlkinPSG", force = TRUE)

The function eq_19() implements Equation (23) and compares, for example, with the exact MOTPW density in Equation (11). To facilitate comparison, the function pdf_theorical() implements this density function. By doing help(eq_19) it is possible to access an example of comparison of the two equations. Note that Equation (23) approximates (11) very well when finite sums are taken in applied problems. In other words, the results achieved by the function eq_19() approximates very well those from pdf_theorical(). The function eq_19() will also allow any baseline cdf $G\left(x\right)$ as an argument of eq_19().

The function eval_plot_moptw() allows to validating numerically Equation (23) by means of plots. The true parameters for the MOTPW density are: $\alpha =1.20$, $\theta =1.50$, $\beta =1.33$, and $\lambda =2$. In addition, we require just a few terms in the sums to obtain a reasonable level of precision as shown in the plots in Figure 5, where six or eight terms provide very accurate approximations.

Figure 5
Numerical evaluation of (19) with finite sums, where N and K denote the upper limits of terms in the related sums with the running indices n and k, respectively.

Properties

We now provide some mathematical properties of $T_s$ that can be easily utilized in the linear representations of the previous section to find the corresponding properties of $X$ and $Y$.

The $n$th ordinary moment of $T_s$ has the form

$${\mu }_n^{\prime }=E\left(T_s^n\right)=s{\int }_{-\infty }^{\infty }{t}^{n}G{(t;𝛕)}^{s-1}g(t;𝛕)dt=s{\int }_0^1{Q}_G{\left(u\right)}^n{u}^{s-1}du,$$(27)

where $Q_G(u;𝛕)=G^{-1}(u;𝛕)$ is the qf of $G$.

Explicit expressions for several exp-G moments can be determined from (27).

The $n$th incomplete moment of $T_s$ follows the previous algebra

$$\begin{array}{c}\hfill m_n\left(y\right)=E(T_s^n|T_s<y)=s{\int }_0^{G(y;𝛕)}{Q}_G{\left(u\right)}^n{u}^{s-1}du,\end{array}$$(28)

where the integral can be calculated for the great majority of G distributions. The first incomplete moment $m_1\left(y\right)$ is the most important case of (28) to find mean deviations and Lorenz and Bonferroni curves.

The moment generating function (mgf) of $T_s$ follows as

\[\label{mgf} M(w)=E(\rm{e}^{w T_s})=s\,\int_{-\infty}^{\infty} \rm{e}^{w t}\,G(t;\pmb{\tau})^{s-1}\,g(t;\pmb{\tau}) dt=s\,\int_{0}^{1} \exp\left[w\, Q_G(u)\right]\,u^{s-1} du.\](29)

The mgfs of exp-G distributions con be determined from Equation (29).

Estimation

The MLEs are appropriate at least in large samples to determine confidence intervals for the parameters. We consider the random variable $X$ defined from Equations (3) and (5) for any baseline G with any unknown parameter vector $\psi =(\alpha ,\theta ,\tau \tau {)}^T$. By simple differentiation of (5), the density of $X$ takes the form

$$\begin{array}{c}\hfill f_X(x;\alpha ,\theta ,𝛕)=\frac{\alpha \theta g(x;𝛕)C^{\prime }\left(\theta H(x;\alpha ,𝛕)\right)}{C\left(\theta \right){[1-\bar{\alpha }\bar{G}(x;𝛕)]}^2},\end{array}$$(30)

where $C^{\prime }(\cdot )$ follows from (4) and $H(x;\alpha ,𝛕)=G(x;𝛕)/[1-\bar{\alpha }\bar{G}(x;𝛕)]$.

The log-likelihood function for $\psi$ from a random sample $x_1,\dots ,x_n$ of $X$ is

$$\begin{array}{rlr}\ell & =& \ell (\psi )=\log \left[\frac{\alpha \theta }{C(\theta )}\right]+\sum _{i=1}^n\log [g(x_i;\tau \tau )]+\sum _{i=1}^n\log \left[C^{\mathrm{\prime }}(\theta H(x_i;\alpha ,\tau \tau ))\right]\\ & -& 2\sum _{i=1}^n\log [1-\overline{\alpha }\overline{G}(x;\tau \tau )].\end{array}$$(31)

A similar development can be conducted for the random variable $Y$ defined from Equation (6) for any baseline G.

We can find the MLE $\hat{\psi }$ by maximizing Equation (31) using the MaxBFGS sub-routine (Ox program), optim function (R), and PROC NLMIXED (SAS). The AdequacyModel package can also maximize (31) using the PSO (particle swarm optimization) approach from the quasi-Newton BFGS, Nelder-Mead and simulated-annealing methods to maximize the log-likelihood function and it does not require initial values. Details are available at Marinho et al. 2019 and https://github.com/prdm0/AdequacyModel.

These scripts can be executed for a wide range of initial values and may lead to more than one maximum. However, in these cases, we consider the MLEs corresponding to the largest value of the maximum log-likelihood. There are sufficient conditions for the existence of these estimates such as compactness of the parameter space and the concavity of the log-likelihood function, but they can exist even when the conditions are not satisfied. In general, there is no explicit solution for the estimates from maximizing (31), but we can establish theoretical conditions on their existence and uniqueness for very special models by examining the ranges of the score components.

Regression

Consider that $X_1,\dots ,X_n$ are independent random variables from any distribution in (11) assuming that the parameters $\lambda$ and $\lambda$ vary through them. We propose a new regression based on the response variable in (11) with the systematic components

$$\begin{array}{r}{\lambda }_i=\exp ({\mathbf{v}}_i^T{\eta }_1)\text{and}{\beta }_i=\exp ({\mathbf{v}}_i^T{\eta }_2),i=1,\dots ,n,\end{array}$$(32)

respectively, where \({\bf v}_{i}^{T}=(v_{i1},\ldots,v_{ip})\), ${\eta }_1=({\eta }_{11},\dots ,{\eta }_{1p}{)}^T$ and ${\eta }_2=({\eta }_{21},\dots ,{\eta }_{2p}{)}^T$. Equations (11) and (32) define the MOTPW regression. For $\alpha =1$, it follows the truncated Poisson Weibull (TPW) regression.

In a similar manner, we can construct many other regressions based on other MOPS-G distributions defined from Equations (5) and (6).

The log-likelihood function for the vector $\psi =(\alpha ,\theta ,{\eta }_1^T,{\eta }_2^T{)}^T$ from the MOTPW regression can be reduced to

$$\begin{array}{rlr}l(\psi )& =& n\log \left[\frac{\alpha \theta }{\exp (\theta )-1}\right]+\sum _{i=1}^n\log ({\beta }_i)+\sum _{i=1}^n{\beta }_i\log ({\lambda }_i)+\sum _{i=1}^n({\beta }_i-1)\log (x_i)-\\ & & \sum _{i=1}^n({\lambda }_i{x}_i{)}^{{\beta }_i}-\sum _{i=1}^n\log \{1-\overline{\alpha }\exp [-({\lambda }_i{x}_i{)}^{{\beta }_i}]\}+\theta \sum _{i=1}^n\frac{\{1-\exp [-({\lambda }_i{x}_i{)}^{{\beta }_i}]\}}{\{1-\overline{\alpha }\exp [-({\lambda }_i{x}_i{)}^{{\beta }_i}]\}}.\end{array}$$(33)

We obtain the MOTPW distribution for ${\lambda }_i=\lambda$ and ${\beta }_i=\beta$.

Let $\hat{\psi }$ be the MLE of $\psi$. Equation (33) can also be maximized using the gamlss regression framework (Stasinopoulos & Rigby 2008) in R.

Two simulation studies

We perform two simulation studies. The first one examines the accuracy of the MLEs of the parameter estimates in the MOTPW distribution. The second one does the same for the MOTPW regression.

The MOTPW distribution

First, we evaluate the precision of the estimates in the MOTPW distribution based on 1,000 Monte Carlo simulations using the R software. The simulation procedure follows as:

• The inverse function $Q\left(u\right)=F^{-1}\left(u\right)$ comes from (7)

  $$Q\left(u\right)={\lambda }^{-1}{\{-log\left[\frac{\theta -log[uexp\left(\theta \right)-u+1]}{\theta +\alpha log[uexp\left(\theta \right)-u+1]-log[uexp\left(\theta \right)-u+1]}\right]\}}^{\frac{1}{\beta }}.$$(34)

• Generate $u\sim$ $U(0,1)$ and obtain the values $x=Q\left(u\right)$ of the MOTPW distribution.

The true parameters are $\lambda =3$, $\beta =1$, $\theta =1.5$ and $\alpha =0.7$. The average estimates (AEs), biases, and mean squared errors (MSEs) are listed in Table I. The three measures decrease steadily when $n$ becomes large.

Table I
Simulation results for the MOTPW distribution.

	$n=100$				$n=250$
Parameter	AE	Bias	MSE		AE	Bias	MSE
$\lambda$	3.001	0.001	0.005		2.996	-0.00	0.005
$\beta$	0.998	-0.002	0.013		1.004	0.004	0.012
$\theta$	1.585	0.085	0.081		1.567	0.0667	0.064
$\alpha$	0.569	-0.130	0.090		0.563	-0.137	0.089
	$n=500$				$n=1,000$
Parameter	AE	Bias	MSE		AE	Bias	MSE
$\lambda$	2.994	-0.006	0.004		2.995	-0.006	0.003
$\beta$	1.006	0.006	0.008		1.007	0.007	0.005
$\theta$	1.546	0.046	0.035		1.526	0.026	0.016
$\alpha$	0.597	-0.103	0.077		0.632	-0.068	0.063

The MOTPW regression

We perform some Monte Carlo simulations for some values of $n$ to investigate the accuracy of the MLEs in the MOTPW regression under four scenarios: Scenario 1: $\theta =0.6$ and $\alpha =0.4$; Scenario 2: $\theta =0.6$ and $\alpha =1.4$; Scenario 3: $\theta =1.7$ and $\alpha =0.4$; Scenario 4: $\theta =1.7$ and $\alpha =1.4$. We take values greater than and less than one for $\theta$ and $\alpha$.

The explanatory variables $v_1,\dots ,v_n$ are generated in the regression by taking ${\lambda }_i=0.5+0.8v_i$, ${\beta }_i=0.3+0.1v_i$, and $v_i\sim Bernoulli\left(0.5\right)$.

For each scenario and value of $n$, one thousand samples are generated from the MOTPW regression fitted to each generated data set. The quantities reported in Table II are in good agreement with the asymptotic results for the MLEs.

Table II
Simulation results for the MOTPW regression.

scenario 1
	$n=100$				$n=500$				$n=1,000$
Parameter	AE	Bias	MSE		AE	Bias	MSE		AE	Bias	MSE
${\gamma }_{10}$	0.614	0.114	0.074		0.557	0.057	0.029		0.527	0.027	0.017
${\gamma }_{11}$	0.785	-0.015	0.031		0.792	-0.008	0.006		0.798	-0.002	0.002
${\gamma }_{20}$	0.256	-0.044	0.031		0.271	-0.029	0.012		0.285	-0.015	0.007
${\gamma }_{21}$	0.101	0.001	0.031		0.101	0.001	0.006		0.102	0.002	0.003
$\theta$	0.734	0.134	0.164		0.651	0.051	0.093		0.621	0.021	0.070
$\alpha$	0.477	0.077	0.089		0.440	0.040	0.072		0.413	0.013	0.067
scenario 2
	$n=100$				$n=500$				$n=1,000$
Parameter	AE	Bias	MSE		AE	Bias	MSE		AE	Bias	MSE
${\gamma }_{10}$	0.684	0.184	0.149		0.567	0.067	0.045		0.528	0.028	0.025
${\gamma }_{11}$	0.779	-0.021	0.044		0.790	-0.009	0.007		0.798	-0.002	0.004
${\gamma }_{20}$	0.235	-0.065	0.045		0.272	-0.028	0.015		0.289	-0.011	0.009
${\gamma }_{21}$	0.096	-0.004	0.035		0.103	0.003	0.007		0.100	0.000	0.003
$\theta$	0.637	0.037	0.157		0.578	-0.023	0.086		0.558	-0.042	0.077
$\alpha$	1.722	0.322	0.382		1.530	0.130	0.180		1.467	0.067	0.123
scenario 3
	$n=100$				$n=500$				$n=1,000$
Parameter	AE	Bias	MSE		AE	Bias	MSE		AE	Bias	MSE
${\gamma }_{10}$	0.337	-0.161	0.079		0.483	-0.017	0.028		0.494	-0.007	0.019
${\gamma }_{11}$	0.819	0.019	0.019		0.802	0.002	0.005		0.799	-0.001	0.002
${\gamma }_{20}$	0.465	0.165	0.069		0.323	0.023	0.014		0.311	0.012	0.009
${\gamma }_{21}$	0.094	-0.006	0.033		0.101	0.001	0.005		0.101	0.001	0.003
$\theta$	1.349	-0.350	0.258		1.643	-0.057	0.035		1.679	-0.022	0.015
$\alpha$	0.460	0.060	0.083		0.429	0.029	0.079		0.407	0.007	0.069
scenario 4
	$n=100$				$n=500$				$n=1,000$
Parameter	AE	Bias	MSE		AE	Bias	MSE		AE	Bias	MSE
${\gamma }_{10}$	0.549	0.049	0.132		0.551	0.051	0.036		0.495	-0.005	0.015
${\gamma }_{11}$	0.796	-0.004	0.038		0.795	-0.006	0.006		0.798	-0.002	0.003
${\gamma }_{20}$	0.332	0.032	0.054		0.286	-0.014	0.012		0.307	0.007	0.006
${\gamma }_{21}$	0.096	-0.004	0.032		0.100	0.000	0.005		0.103	0.003	0.003
$\theta$	1.406	-0.294	0.222		1.643	-0.057	0.029		1.684	-0.016	0.013
$\alpha$	1.913	0.513	0.739		1.604	0.204	0.240		1.408	0.008	0.090

Residual analysis

We investigate the quantile residuals (qrs) to verity the adequacy of the response distribution to determine outliers in the MOTPW regression. The same approach can be adopted to many other regressions defined from the distributions in (5) and (6). The qrs are given by (Dunn & Smyth 1996)

$$\begin{array}{c}\hfill {qr}_i={\Phi }^{-1}\left\{{[exp\left(\theta \right)-1]}^{-1}exp\left\{\theta \frac{1-exp[-{\left({\lambda }_i{x}_i\right)}^{{\beta }_i}]}{1-\bar{\alpha }exp[-{\left({\lambda }_i{x}_i\right)}^{{\beta }_i}]}\right\}-1\right\},\end{array}$$(35)

where $\Phi (\cdot )$ is the normal cdf and ${\lambda }_i$ and ${\beta }_i$ are defined in Equation (32).

We consider the same scenarios for the simulations in Section Two Simulation Studies. For each fitted regression, the qrs are calculated from Equation (35). Figures 6, 7, 8, and 9 display QQ plots which show that the empirical distribution of these residuals is close to the standard normal distribution.

Figure 6
QQ plots for scenario 1 ($𝛉=0.6$ and $𝛂=0.4$). (a) $𝐧=100$. (b) $𝐧=500$. (c) $𝐧=1,000$.

Figure 7
QQ plots for scenario 2 ($𝛉=0.6$ and $𝛂=1.4$). (a) $𝐧=100$. (b) $𝐧=500$. (c) $𝐧=1,000$.

Figure 8
QQ plots for scenario 3 ($𝛉=1.7$ and $𝛂=0.4$). (a) $𝐧=100$. (b) $𝐧=500$. (c) $𝐧=1,000$.

Figure 9
QQ plots for scenario 4 ($𝛉=1.7$ and $𝛂=1.4$). (a) $𝐧=100$. (b) $𝐧=500$. (c) $𝐧=1,000$.

Applications

The beta Weibull (BW) and Kumaraswamy Weibull (KwW) distributions have been widely used to fit real data in the last ten years or so. We compare the MOTPW distribution with the BW and KwW distributions since all of them have four parameters. The BW density pioneered by Lee et al. 2007 is

\begin{aligned} f(x)=\frac{c\lambda^c}{B(a,b)}x^{c-1} {\rm exp}\{-b(\lambda x)^c\}[1-{\rm exp}\{-(\lambda x)^c\}]^{a-1},\,\,x>0,\end{aligned}

where all parameters are positive.

The KwW density introduced by Cordeiro & Castro 2011 has the form

$$\begin{array}{c}\hfill f\left(x\right)=abc{\lambda }^c{x}^{c-1}exp\{-{\left(\lambda x\right)}^c\}{[1-exp\{-{\left(\lambda x\right)}^c\}]}^{a-1}\{1-{[1-exp\{-{\left(\lambda x\right)}^c\}]}^a{\}}^{b-1},x>0,\end{array}$$

where all parameters are positive.

Application 1: Hourly dollar wage data

The first application refers to hourly dollar wages for $n=534$ US workers. These data are obtained from the SemiPar package (Wand et al. 2005). Table III lists the estimates, standard errors (SEs) in parentheses, and three classical statistics. The lowest values of these measures reveal that the MOTPW is the best model. Next, the likelihood ratio (LR) statistic for comparing the MOTPW and TPW models is $6.159(\mathrm{\text{p-value}}<0.013)$ which supports the wider distribution.

Figure 10a shows the histogram and the estimated MOTPW density. Figure 10b provides the empirical function and estimated MOTPW cdf, thus revealing that this distribution is appropriate for these data.

Figure 10
(a) Estimated MOTPW pdf. (b) Estimated MOTPW cdf and the empirical cdf.

Table III
Results for hourly dollar wage data.

Model	$log\left(\lambda \right)$	$log\left(\beta \right)$	$\theta$	$\alpha$	AIC	BIC	GD
MOTPW	-2.720	0.694	11.210	0.019	3031.288	3048.410	3023.288
	(0.141)	(0.085)	(3.020)	(0.007)
TPW	0.248	-0.541	31.100	(-)	3035.448	3048.289	3029.448
	(0.444)	(0.112)	(12.540)	(-)
Model	$log\left(\lambda \right)$	$a$	$b$	$log\left(c\right)$	AIC	BIC	GD
KwW	-0.601	12.124	0.317	0.060	3034.039	3051.160	3026.039
	(0.023)	(0.802)	(0.013)	(0.014)
BW	-5.453	2.327	126.000	0.216	3084.086	3101.208	3076.086
	(0.020)	(0.067)	(0.009)	(0.007)

Application 2: Diabetes data

We consider two variables from the data reported by Reaven & Miller 1979: the response $x_i$ is the relative weight defined by the ratio between the actual weight and the expected weight (given the person’s height), and the explanatory variable $v_{i1}$ indicates the diagnostic group (0 =normal, 1= chemical diabetes, 2 = overt diabetes). The diagnostic group has three levels and then we have two dummy variables $\left(d_{ij}\right)$ (for $i=1,\dots ,145$ and $j=1,2$). The objective is to know what are the relations among the relative weight and the levels of the diagnostic group.

The systematic components for the MOTPW regression are

$${\lambda }_i=exp({\eta }_{10}+{\eta }_{11}{d}_{i1}+{\eta }_{12}{d}_{i2})\mathrm{\text{and}}{\beta }_i=exp({\eta }_{20}+{\eta }_{21}{d}_{i1}+{\eta }_{22}{d}_{i2}),i=1,\dots ,145.$$

The measures for the fitted regressions are reported in Table IV. Clearly, the MOTPW is the best regression for these data.

Table IV
Measures for diabetes data.

Model	AIC	BIC	GD
MOTPW	-194.316	-170.502	-210.316
TPW	-191.726	-170.889	-205.726
KwW	-188.769	-164.955	-204.769
BW	-185.607	-161.793	-201.607

Table V provides the estimates, SEs and $p$-values for the best regression.

Table V
Results for diabetes data.

Parameter	Estimate	SE	p-Value
${\eta }_{10}$	0.065	0.093	0.489
${\eta }_{11}$	-0.119	0.036	0.001
${\eta }_{12}$	-0.049	0.028	0.082
${\eta }_{20}$	1.719	0.245	$<$0.001
${\eta }_{21}$	0.373	0.140	0.009
${\eta }_{22}$	0.131	0.141	0.355
$\theta$	12.401	8.866
$\alpha$	0.095	0.079

We note that the co-variable $d_{i1}$ is significant and $d_{i2}$ is not. So, there is a real difference between normal and chemical diabetes groups in relation to relative weight and no difference between normal and overt diabetes groups to relative weight. The same findings can be seen in Figure 12.

The LR statistic to compare the MOTPW and TPW regressions is $w=4.590$ ($p$-value=0.032) that indicates that the fist regression is superior to the second regression to these data in terms of model fitting.

The plot of the residuals reported in Figure 11a does not detect outliers and departures from the general assumptions. The worm plot (Buuren & Fredriks 2001) of the residuals in Figure 11b and the QQ plot displayed in Figure 11c show the adequacy of the MOTPW regression for the current data.

Figure 11
(a) Residual plot. (b) Worm plots. (c) QQ plot.

A graphical comparison from the estimated cdfs in Figure 12 also supports the regression analysis.

Figure 12
Estimated cdf and the empirical cdf.

Conclusions

We define two flexible Marshall–Olkin–Power-Series (MOPS) families of continuous distributions which can be very useful to fit real data. They are obtained by combining the Marshall–Olkin class (Marshall & Olkin 1997) and the power series distribution. Hundreds of continuous distributions can be easily formulated from the two families. We discuss some special distributions and maximum likelihood estimation. We introduce the Marshall–Olkin Truncated Poisson Weibull regression associated with one of the families. Some mathematical properties of these families are presented. We provide a package implemented in R software which can be used to determine numerically some mathematical properties for any distribution in the new families. The utility of the proposed models is proved empirically in two applications.

ACKNOWLEDGMENTS

We gratefully acknowledge from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil.

• BUUREN S VAN & FREDRIKS M. 2001. Worm plot: a simple diagnostic device for modelling growth reference curves. Stat Med 20(8): 1259–1277.
• CORDEIRO GM & CASTRO M DE. 2011. A new family of generalized distributions. J Stat Comput Simul 81(7): 883–898.
• DUNN PK & SMYTH GK. 1996. Randomized quantile residuals. J Comput Graph Stat 5(3): 236–244.
• LEE C, FAMOYE F & OLUMOLADE O. 2007. Beta-Weibull distribution: some properties and applications to censored data. J Mod Appl Stat Meth 6(1): 173–186.
• MARINHO PRD, SILVA RB, BOURGUIGNON M, CORDEIRO GM & NADARAJAH S. 2019. AdequacyModel: An R package for probability distributions and general purpose optimization. PloS One 14(8): e0221487.
• MARSHALL AW & OLKIN I. 1997. A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families. Biometrika 84(3): 641–652.
• R CORE TEAM. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. URL https://www.R-project.org/
  » https://www.R-project.org/
• REAVEN G & MILLER R. 1979. An attempt to define the nature of chemical diabetes using a multidimensional analysis. Diabetologia 16(1): 17–24.
• STASINOPOULOS DM & RIGBY RA. 2008. Generalized additive models for location scale and shape (GAMLSS) in R. J Stat Soft 23: 1–46.
• TAHIR MH & NADARAJAH S. 2015. Parameter induction in continuous univariate distributions: Well-established G families. An Acad Bras Cienc 87: 539–568.
• WAND M, COULL B, FRENCH J, GANGULI B, KAMMANN E, STAUDENMAYER J & ZANOBETTI A. 2005. SemiPar 1.0. R package. URL http://cran.r-project.org
  » http://cran.r-project.org

Publication Dates

• Publication in this collection
  18 July 2022
• Date of issue
  2022

History

• Received
  28 Dec 2020
• Accepted
  7 June 2021