The Bivariate Unit-Sinh-Normal Distribution and Its Related Regression Model

Abstract

In this paper, a new bivariate absolutely continuous probability distribution is introduced. The new distribution, which is called the bivariate unit-sinh-normal (BVUSHN) distribution, arises by applying a transformation to the bivariate Birnbaum–Saunders distribution (BVBS). The main properties of the new proposal are studied in detail. In addition, from the new distribution, the BVUSHN regression model is also introduced. For both the bivariate probability distribution and the respective associated regression model, parameter estimation is conducted from a classical approach by using the maximum likelihood method together with the two-step estimation method. A small Monte Carlo simulation study is carried out to evaluate the behavior of the used estimation method and the properties of the estimators. Finally, for illustrative purposes, two applications with real data are presented in which the usefulness of the proposals is evidenced.

1. Introduction

Birnbaum and Saunders [1,2] developed an important two-parameter lifetime distribution for fatigue failure caused under dynamic load. The probability density function (pdf) of the Birnbaum–Saunders (BS) distribution is given by

$f_T\left(t\right)=\varphi \left(a_t\right)A_t,t>0,$

where $a_t={\alpha }^{-1}(\sqrt{t/\beta }-\sqrt{\beta /t})$, $A_t=t^{-3/2}(t+\beta )/\left[2\alpha \sqrt{\beta }\right]$, $\varphi (·)$ is the pdf of the standard normal distribution, $\alpha >0$ is the shape parameter, and $\beta >0$ is a scale parameter. We use the conventional notation to denote this distribution, that is, $T\sim \mathrm{BS}(\alpha ,\beta )$. Extensions of the BS distribution have been studied for many authors, for example, [3,4,5,6,7,8,9].

In addition to the work mentioned above, recent proposals were made by Balakrishnan and Kundu [10], who made a detailed review of the developments that have taken place in relation to the BS distribution. In particular, the authors addressed different aspects, including interpretations, constructions, generalizations, inference processes, and extensions to multivariate cases, among others. The authors also presented some open problems that could be considered in future research. Furthermore, Athayde et al. [11] studied the failure rate shape of the Birnbaum–Saunders-logistic distribution. In addition, the authors discussed and compared some robustness issues related to the Birnbaum–Saunders, Birnbaum–Saunders-logistic, and Birnbaum–Saunders-t distributions. On the other hand, Arrué et al. [12] built an extension of the BS distribution from the skew-normal model and gained flexibility in skewness and kurtosis. The characterization of this new distribution was presented and they illustrated the usefulness of their proposal with a motivating example of fatigue life data.

A more recent extension of the BS distribution to the case of bimodal data was presented by Reyes et al. [13]. This proposal has proven to be useful for the analysis of data from the environment and medical sciences, which demonstrates the wide field of applicability for this model in the different areas of knowledge by considering not only data from engineering, for example.

The BS distribution has also been extended to models suitable for fitting data in the interval $(0,1)$ as rates, proportions, or indices. Mazucheli et al. [14], for example, presented a type of BS distribution with support in the unit interval $(0,1)$, which is a new alternative to the beta and Kumaraswamy distributions. These authors have called it as unit-Birnbaum–Saunders (UBS) distribution. The UBS distribution is obtained from the transformation $X=exp(-T)$, where $T\sim \mathrm{BS}(\alpha ,\beta )$. The pdf of the UBS distribution is given by

$f\left(x\right)=\frac{1}{2x\alpha \beta \sqrt{2\pi }}\left[{\left(-\frac{\beta }{logx}\right)}^{1/2}+{\left(-\frac{\beta }{logx}\right)}^{3/2}\right]exp\left\{\frac{1}{2{\alpha }^2}\left[\frac{logx}{\beta }+\frac{\beta }{logx}+2\right]\right\},$

where $x\in (0,1)$, $\alpha >0$ is the shape parameter, and $\beta >0$ is a scale parameter. The corresponding cumulative distribution function (cdf) is given by

$F\left(x\right)=1-\mathsf{\Phi }\left[\frac{1}{\alpha }\left\{{\left(-\frac{logx}{\beta }\right)}^{1/2}-{\left(-\frac{\beta }{logx}\right)}^{1/2}\right\}\right],x\in (0,1),$

where $\mathsf{\Phi }(·)$ is the cdf of the standard normal distribution. We use the notation $X\sim \mathrm{UBS}(\alpha ,\beta )$.

The BS distribution has also been used to define linear regression models, known as log-BS models. In this framework, it is assumed that $Y=log\left(T\right)$, where $T\sim \mathrm{BS}(\alpha ,\beta )$, and the errors of the log-linear model follow the sinh-normal (SHN) distribution of Rieck and Nedelman [15]. The pdf of the SHN distribution is given by

$\phi \left(y\right)=B_y\varphi \left(b_y\right),y\in \mathbb{R},$

where $b_y=b(y;\alpha ,\mu ,\sigma )=(2/\alpha )sinh((y-\mu )/\sigma )$, $B_y=B(y;\alpha ,\mu ,\sigma )={\sigma }^{-1}(2/\alpha )cosh((y-\mu )/\sigma )$, $\alpha >0$ is the shape parameter, $\mu$ is a location parameter, and $\sigma >0$ is a scale parameter. We use the notation $Y\sim \mathrm{SHN}(\alpha ,\mu ,\sigma )$. Extensions of the SHN distribution are reported in Martínez-Flórez et al. [8], Barros et al. [16], Leiva et al. [17], Lemonte [18], and and Santana et al. [19].

Extensions of the BS distribution to the bivariate setup were given by Kundu et al. [20] and, for the multivariate case, Lemonte et al. [21] introduced the multivariate skew-normal distribution. Likewise, for the log-BS regression model, Lemonte [22] proposed the multivariate BS regression model, Kundu [23] studied the bivariate SHN distribution, and Martínez-Flórez et al. [24] considered the asymmetric multivariate BS regression model. In addition, Martínez-Flórez et al. [25] proposed an exponentiated multivariate extension of the BS log-linear model, that is, an extension to the alpha-power family of distribution.

To study the relationship between variables when the dependent variable is observed in the unit interval, such as indices, rates, or proportions, the most widely known and used model is the beta regression model. Other authors developed works in this same direction using the beta distribution; see, for example, [26,27,28,29,30,31,32]. On the other hand, situations in which the limits of the interval $(0,1)$ are included in the support of the response variable have also been considered; see, for example, [33,34,35,36,37].

The chief goal of this paper is to introduce a bivariate regression model to deal with response variables in the unit interval. To do so, we first introduce a bivariate probability distribution supported on the unit interval which arises from the bivariate SHN distribution. The article is organized as follows. Section 2 presents the bivariate BS and SHN distributions. The bivariate UBS distribution is introduced in Section 3, and some of its structural properties are also derived. The results of a simulation study and its respective discussion are presented in Section 6. In Section 4, we introduce the bivariate unit SHN distribution, and we derive some of its properties. The bivariate unit SHN regression model is introduced in Section 5. Applications to real data are provided in Section 7.

2. The Bivariate BS and SHN Distributions

Kundu et al. [20] extended the univariate BS distribution to the bivariate case. The bivariate random vector $(T_1,T_2)$ follows a bivariate BS distribution if it can be expressed as

$(T_1,T_2)=\left({\beta }_1{\left[\frac{{\alpha }_1}{2}{Z}_1+\sqrt{{\left(\frac{{\alpha }_1}{2}{Z}_1\right)}^2+1}\right]}^2,{\beta }_2{\left[\frac{{\alpha }_2}{2}{Z}_2+\sqrt{{\left(\frac{{\alpha }_2}{2}{Z}_2\right)}^2+1}\right]}^2\right),$

where the bivariate random vector $(Z_1,Z_2)$ follows a bivariate normal distribution with parameters $(0,0,1,1,\rho )$ and $-1<\rho <1.$ It follows that the cdf of the bivariate BS distribution is given by

$F_{T_1,T_2}(t_1,t_2)={\mathsf{\Phi }}_2(a_{t_1},a_{t_2};\rho ),t_1>0,t_2>0,$

where ${\mathsf{\Phi }}_2(·)$ is the cdf of the bivariate normal distribution and 

$a_{t_j}=a_{t_j}({\alpha }_j,{\beta }_j)=\frac{1}{{\alpha }_j}\left(\sqrt{\frac{t_j}{{\beta }_j}}-\sqrt{\frac{{\beta }_j}{t_j}}\right),j=1,2.$

Hence, the joint pdf of $(T_1,T_2)$ can be written as

$f_{T_1,T_2}(t_1,t_2)={\varphi }_2(a_{t_1},a_{t_2};\rho )A_{t_1}{A}_{t_2},$

where ${\varphi }_2(x,y;\rho )$ denotes the pdf of the bivariate standard normal distribution and

$A_{t_j}=\frac{t_j^{-3/2}[t_j+{\beta }_j]}{2{\alpha }_j{\beta }_j^{1/2}},j=1,2.$

We use the notation $(T_1,T_2)\sim \mathrm{BVBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$. Similarly, the bivariate SHN distribution of the random vector $(Y_1,Y_2)$, studied by Kundu [23], is defined by the joint distribution function

$G_{Y_1,Y_2}(y_1,y_2)={\mathsf{\Phi }}_2(b_{y_1},b_{y_2};\rho ),y_1\in \mathbb{R},y_2\in \mathbb{R},$

where $b_{y_j}=b(y_j;{\alpha }_j,{\mu }_j,{\sigma }_j)$ for $j=1,2$. Hence, the joint pdf of $(Y_1,Y_2)$ can be written as

$g_{Y_1,Y_2}(y_1,y_2)={\varphi }_2(b_{y_1},b_{y_2};\rho )B_{y_1}{B}_{y_2},y_1\in \mathbb{R},y_2\in \mathbb{R},$

where $B_{y_j}=B(y_j;{\alpha }_j,{\mu }_j,{\sigma }_j)$. The notation $(Y_1,Y_2)\sim \mathrm{BVSHN}({\alpha }_1,{\alpha }_2,{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,\rho )$ is used to denote this distribution. The reader is referred to Kundu et al. [20] and Kundu [23] for more details about the BVBS and BVSHN distributions, respectively.

3. The Bivariate UBS Distribution

We now introduce the bivariate UBS distribution and study some of its properties and the process of estimating its parameters. The bivariate random vector $(X_1,X_2)$ is said to have a bivariate UBS distribution with parameters ${\alpha }_1$, ${\beta }_1$, ${\alpha }_2$, ${\beta }_2$, and $\rho$ if the joint cdf of $(X_1,X_2)$ is given by

(1)$F_{X_1,X_2}(x_1,x_2)=1-\mathsf{\Phi }\left({\lambda }_{x_1}\right)-\mathsf{\Phi }\left({\lambda }_{x_2}\right)+{\mathsf{\Phi }}_2({\lambda }_{x_1},{\lambda }_{x_2};\rho ),$

where $x_1\in (0,1)$, $x_2\in (0,1)$, and 

${\lambda }_{x_j}={\lambda }_{x_j}({\alpha }_j,{\beta }_j)=\frac{1}{{\alpha }_j}\left(\sqrt{-\frac{log{x}_j}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{x}_j}}\right).$

Hence, the joint pdf of $(X_1,X_2)$ which is obtained by deriving (1) can be expressed as

(2)$f_{X_1,X_2}(x_1,x_2)={\varphi }_2({\lambda }_{x_1},{\lambda }_{x_2};\rho ){\mathsf{\Lambda }}_{x_1}{\mathsf{\Lambda }}_{x_2},$

where

${\mathsf{\Lambda }}_{x_j}=\frac{{(-log{x}_j)}^{-3/2}[{\beta }_j-log{x}_j]}{2x_j{\alpha }_j{\beta }_j^{1/2}},j=1,2.$

For the joint pdf given in (2), we use the notation $\mathrm{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$, where $\rho$ is the correlation coefficient. The bivariate distribution $(X_1,X_2)\sim \mathrm{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$ is obtained from the vector transformation $(X_1,X_2)=(exp(-T_1),exp(-T_2))$, where $(T_1,T_2)\sim \mathrm{BVBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho ).$

From the previous definition, it follows that $(Z_1,Z_2)\sim \mathrm{BN}(0,0,1,1,\rho )$ with

$(Z_1,Z_2)=\left(\frac{1}{{\alpha }_1}\left(\sqrt{-\frac{log{X}_1}{{\beta }_1}}-\sqrt{-\frac{{\beta }_1}{log{X}_1}}\right),\frac{1}{{\alpha }_2}\left(\sqrt{-\frac{log{X}_2}{{\beta }_2}}-\sqrt{-\frac{{\beta }_2}{log{X}_2}}\right)\right),$

where $\mathrm{BN}(0,0,1,1,\rho )$ denotes the bivariate standard normal distribution with correlation coefficient $\rho$. Therefore, for the marginal distributions, it follows that

$Z_j=\frac{1}{{\alpha }_j}\left(\sqrt{-\frac{log{X}_j}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{X}_j}}\right)\sim \mathrm{N}(0,1)\mathrm{and}{Z}_j^2\sim {\chi }_1^2,\mathrm{for}j=1,2,$

and, hence, if $(X_1,X_2)\sim \mathrm{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$, then

$(T_1,T_2)=(-log{X}_1,-log{X}_2)\sim \mathrm{BVBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho ).$

We have the following theorem.

Theorem 1.

Let $(X_1,X_2)\sim \mathit{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$. Hence:

• (i) 

  $X_j\sim \mathit{UBS}({\alpha }_j,{\beta }_j),$ for $j=1,2.$

• (ii) 

  The cdf of $X_1$, given $X_2=x_2$, is

  $F_{X_1\mid X_2=x_2}(x_1\mid X_2=x_2)=1-\mathsf{\Phi }\left(\frac{{\lambda }_{x_1}-\rho {\lambda }_{x_2}}{\sqrt{1-{\rho }^2}}\right).$

• (iii) 

  The pdf of $X_1$, given $X_2=x_2$, is

  $f_{X_1\mid X_2=x_2}(x_1\mid X_2=x_2)=\frac{1}{\sqrt{2\pi }\sqrt{1-{\rho }^2}}exp\left(-\frac{1}{2}{\left(\frac{{\lambda }_{x_1}-\rho {\lambda }_{x_2}}{\sqrt{1-{\rho }^2}}\right)}^2\right).$

• (iv) 

  The joint survival function of $(X_1,X_2)$ is

  $S_{X_1,X_2}(x_1,x_2)=\mathsf{\Phi }\left({\lambda }_{x_1}\right)+\mathsf{\Phi }\left({\lambda }_{x_2}\right)-{\mathsf{\Phi }}_2({\lambda }_{x_1},{\lambda }_{x_2};\rho ).$

Remark 1.

We can generate random samples from the BVUBS $({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$ distribution according to the algorithm described below:

• (i) 

  For $j=1,2$, obtain $T_j$ following the described algorithm in Kundu et al. [20].

• (ii) 

  For each value of $T_j$, compute

  $X_j=exp(-T_j),j=1,2.$

  Hence, $(X_1,X_2)\sim \mathit{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$.

Since the marginal $X_j\sim \mathrm{UBS}({\alpha }_j,{\beta }_j)$ for $j=1,2$, the mean and variance of $X_j$ are

$\mathbb{E}\left(X_j\right)=\frac{2{\alpha }_j^2{\beta }_j+\sqrt{2{\alpha }_j^2{\beta }_j+1}+1}{2(2{\alpha }_j^2{\beta }_j+1)}exp\left[-\frac{\sqrt{2{\alpha }_j^2{\beta }_j+1}-1}{{\alpha }_j^2}\right],$

$Var\left(X_j\right)=\frac{4{\alpha }_j^2{\beta }_j+\sqrt{4{\alpha }_j^2{\beta }_j+1}+1}{2(4{\alpha }_j^2{\beta }_j+1)}exp\left[-\frac{\sqrt{4{\alpha }_j^2{\beta }_j+1}-1}{{\alpha }_j^2}\right]-{\mathbb{E}}^2\left(X_j\right).$

The product moments of $(X_1,X_2)$, say, $m_{rs}=\mathbb{E}\left(X_1^r{X}_2^s\right)$, are complicated to obtain in closed-form expression and, hence, they have to be obtained from numerical calculations. However, for $\rho =0$, it follows that

$\mathbb{E}\left(X_1{X}_2\right)=exp\left[-\sum _{j=1}^2\frac{\sqrt{2{\alpha }_j^2{\beta }_j+1}-1}{{\alpha }_j^2}\right]\prod _{j=1}^2\frac{2{\alpha }_j^2{\beta }_j+\sqrt{2{\alpha }_j^2{\beta }_j+1}+1}{2(2{\alpha }_j^2{\beta }_j+1)}.$

We say that a function defined in ${\mathbb{R}}^2$ is totally positive of order 2 (abbreviated $T{P}_2$) for $x_1<x_2$ and $y_1<y_2,$ with $x,y\in \mathbb{R}$ if

$g(x_1,y_1)g(x_2,y_2)\ge g(x_2,y_1)g(x_1,y_2).$

We have the following theorem.

Theorem 2.

Let $(X_1,X_2)\sim \mathit{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$. For $0<\rho <1$, we have that $(X_1,X_2)$ is $T{P}_2$.

Proof. 

The proof follows from the fact that for $(X_1,X_2)\sim \mathrm{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$ it follows that $(T_1,T_2)=(-log{X}_1,-log{X}_2)\sim \mathrm{BVBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$. Furthermore, let $f_{\mathrm{BVBS}}(·)$ and $f_{\mathrm{BVUBS}}(·)$ be the joint pdfs of the BVBS and BVUBS distributions, respectively. Then, $f_{\mathrm{BVUBS}}(-log{x}_1,-log{x}_2)=\frac{1}{x_1{x}_2}{f}_{\mathrm{BVBS}}(-log{x}_1,-log{x}_2)$. Since $(T_1,T_2)$ is TP2, then for $0<\rho <1$ and for $x_{11}<x_{12}$ and $x_{21}<x_{22}$, we have that $0<-log{x}_{12}<-log{x}_{11}$ and $0<-log{x}_{22}<-log{x}_{21}$. Note that

$\begin{array}{ccc}& & f_{\mathrm{BVBS}}(-log{x}_{12},-log{x}_{22})f_{\mathrm{BVBS}}(-log{x}_{11},-log{x}_{21})\hfill \\ & & -f_{\mathrm{BVBS}}(-log{x}_{11},-log{x}_{22})f_{\mathrm{BVBS}}(-log{x}_{12},-log{x}_{21})\ge 0,\hfill \end{array}$

and, hence,

$f_{\mathrm{BVUBS}}(x_{11},x_{21})f_{\mathrm{BVUBS}}(x_{12},x_{22})-f_{\mathrm{BVUBS}}(x_{12},x_{21})f_{\mathrm{BVUBS}}(x_{11},x_{22})\ge 0.$

Therefore, for $0<\rho <1$, $(X_1,X_2)$ is TP2.    □

Parameter Estimation

Estimating the unknown parameters of a parametric model lies at the heart of distribution theory and its applications. In this section, the estimation of the BVUBS model parameters is investigated by the maximum likelihood (ML) method.

Let $\{(x_{11},x_{12}),\cdots ,(x_{n1},x_{n2})\}$ be a sample of size n from the $\mathrm{BVUBS}({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )$ distribution. The log-likelihood function for the parameter vector $\mathit{\theta }={({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )}^{\top }$ is given by

$\begin{array}{cc}\hfill \ell \left(\mathit{\theta }\right)& =-nlog\left({\alpha }_1{\alpha }_2{\beta }_1{\beta }_2\right)-\frac{n}{2}log(1-{\rho }^2)\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^{2}log\left[{\left(-\frac{{\beta }_j}{log{x}_{ij}}\right)}^{1/2}+{\left(-\frac{{\beta }_j}{log{x}_{ij}}\right)}^{3/2}\right]\hfill \\ \hfill & -\frac{1}{2(1-{\rho }^2)}\sum _{i=1}^n[\sum _{j=1}^2\frac{1}{{\alpha }_j^2}{\left(\sqrt{-\frac{log{x}_{ij}}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{x}_{ij}}}\right)}^2\hfill \\ \hfill & -2\rho \prod _{j=1}^2\frac{1}{{\alpha }_j}\left(\sqrt{-\frac{log{x}_{ij}}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{x}_{ij}}}\right)].\hfill \end{array}$

We have that

$\sqrt{-\frac{log{X}_j}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{X}_j}}\sim \mathrm{N}(0,{\alpha }_j^2),\mathrm{for}j=1,2.$

Since $(Z_1,Z_2)\sim \mathrm{BN}(0,0,1,1,\rho )$, it follows that

$\left(\sqrt{-\frac{log{X}_1}{{\beta }_1}}-\sqrt{-\frac{{\beta }_1}{log{X}_1}},\sqrt{-\frac{log{X}_2}{{\beta }_2}}-\sqrt{-\frac{{\beta }_2}{log{X}_2}}\right)\sim {\mathrm{N}}_2\left(\left(\genfrac{}{}{0pt}{}{0}{0}\right),\mathsf{\Sigma }\right),$

where

$\mathsf{\Sigma }=\left[\begin{array}{cc}{\alpha }_1^2& {\alpha }_1{\alpha }_2\rho \\ {\alpha }_1{\alpha }_2\rho & {\alpha }_2^2\end{array}\right].$

We have that, given ${\beta }_1$ and ${\beta }_2$, the ML estimators of ${\alpha }_1$, ${\alpha }_2$, and $\rho$ become

${\widehat{\alpha }}_j\left({\beta }_j\right)={\left(\frac{s_j}{{\beta }_j}+\frac{{\beta }_j}{r_j}-2\right)}^{1/2},j=1,2,$

where $s_j=(1/n){\sum }_{i=1}^n(-log{x}_j)$, $r_j={\left[(1/n){\sum }_{i=1}^n{(-log{x}_j)}^{-1}\right]}^{-1}$, and

$\widehat{\rho }({\beta }_1,{\beta }_2)=\frac{{\sum }_{i=1}^n{\prod }_{j=1}^2\left(\sqrt{-\frac{log{x}_{ij}}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{x}_{ij}}}\right)}{\sqrt{{\prod }_{j=1}^2\left[{\sum }_{i=1}^n{\left(\sqrt{-\frac{log{x}_{ij}}{{\beta }_j}}-\sqrt{-\frac{{\beta }_j}{log{x}_{ij}}}\right)}^2\right]}}.$

To estimate ${\beta }_1$ and ${\beta }_2$, we consider the following profiled log-likelihood function

$\begin{array}{ccc}\hfill {\ell }_p({\beta }_1,{\beta }_2)=& -& nlog\left({\beta }_1{\beta }_2{\widehat{\alpha }}_1\left({\beta }_1\right){\widehat{\alpha }}_2\left({\beta }_2\right)\right)-\frac{n}{2}log(1-\widehat{\rho }({\beta }_1,{\beta }_2))\hfill \\ & +& \sum _{i=1}^n\sum _{j=1}^{2}log\left[{\left(-\frac{{\beta }_j}{log{x}_{ij}}\right)}^{1/2}+{\left(-\frac{{\beta }_j}{log{x}_{ij}}\right)}^{3/2}\right].\hfill \end{array}$

The ML estimates for ${\beta }_1$ and ${\beta }_2$ need to be obtained through a numerical maximization using optimization algorithms such as, for example, the Newton–Raphson iterative technique. The optimization algorithms require the specification of initial values to be used in the iterative scheme. Our suggestion is to use as initial guesses the estimates obtained from the modified moments of the marginal distributions of $X_j$ for $j=1,2$, given that $X_j\sim \mathrm{UBS}({\alpha }_j,{\beta }_j)$. More details about these values, which are denoted by ${\overline{\alpha }}_j,{\overline{\beta }}_j$, and $\overline{\rho }$, can be found in Kundu et al. [20].

Let $\widehat{\theta }={({\widehat{\alpha }}_1,{\widehat{\alpha }}_2,{\widehat{\beta }}_1,{\widehat{\beta }}_2,\widehat{\rho })}^{\top }$ denote the ML estimator of $\theta ={({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,\rho )}^{\top }$. For $n\to \infty$, we have that $\sqrt{n}(\widehat{\theta }-\theta )\to {\mathrm{N}}_5({\mathbf{0}}_5,I^{-1}(\theta ))$, where $I\left(\theta \right)$ is the Fisher information matrix, whose elements are the same as those given for the Fisher information matrix of the BVBS distribution provided in Kundu et al. [20], just replacing $T_j$ by $-log{X}_j$.

4. The Bivariate Unit SHN Distribution

The SHN distribution studied by Rieck and Nedelman [15] is defined on the set of real numbers. Similarly, the bivariate SHN distributions of Kundu [23] and Martínez-Flórez et al. [24] are defined on ${\mathbb{R}}^2$. Therefore, its application to sets of variables that are distributed between zero and one would not always be a convenient option.

We now introduce a bivariate distribution for rates or proportions as an alternative to bivariate datasets supported by the unit rectangle $(0,1)\times (0,1)$. Since the SHN distribution is defined on the set of real numbers, to obtain an adequate transformation that allows us to define a random variable with a distribution in the unit interval $(0,1)$, we use the log-log complement transformation. We say that the variable Y follows the unit SHN (USHN) distribution with parameters $(\alpha ,\mu ,\sigma )$ if 

$Z=\frac{2}{\alpha }sinh\left(\frac{log(-log(1-Y\left)\right)-\mu }{\sigma }\right)\sim \mathrm{N}(0,1).$

Thus, the pdf of a random variable Y with a USHN distribution is given by

$\phi \left(y\right)={\mathsf{\Omega }}_y\varphi \left({\omega }_y\right),y\in (0,1)$

with

$\begin{array}{cc}\hfill {\mathsf{\Omega }}_y& =\mathsf{\Omega }(y;\alpha ,\mu ,\sigma )=\frac{\frac{2}{\alpha }cosh\left(\frac{log(-log(1-y\left)\right)-\mu }{\sigma }\right)}{\sigma (1-y)log{(1-y)}^{-1}},\hfill \\ \\ \hfill {\omega }_y& =\omega (y;\alpha ,\mu ,\sigma )=\frac{2}{\alpha }sinh\left(\frac{log(-log(1-y\left)\right)-\mu }{\sigma }\right),\hfill \end{array}$

where $\alpha >0$ is the shape parameter, $\mu$ is a location parameter, and $\sigma >0$ is a scale parameter. We use the notation $Y\sim \mathrm{USHN}(\alpha ,\mu ,\sigma )$.

The cdf of the random variable $Y\sim \mathrm{USHN}(\alpha ,\mu ,\sigma )$ takes the form

$\begin{array}{c}\hfill F\left(y\right)=F_{\mathrm{SHN}}(log(-log(1-y)))=\mathsf{\Phi }\left(\frac{2}{\alpha }sinh\left(\frac{log(-log(1-y\left)\right)-\mu }{\sigma }\right)\right),\end{array}$

where $F_{\mathrm{SHN}}(·)$ is the cdf of the $\mathrm{SHN}(\alpha ,\mu ,\sigma )$ distribution. In addition, if $Z\sim N(0,1)$, then

$Y=1-e^{-e^{\mu +\sigma \ast {sinh}^{-1}\left(\frac{\alpha Z}{2}\right)}}\sim \mathrm{USHN}(\alpha ,\mu ,\sigma ).$

Finally, the generation of a random variable with a $\mathrm{USHN}(\alpha ,\gamma ,\sigma )$ distribution can be obtained from the previous expression.

Bivariate Extension

The bivariate USHN distribution can be directly obtained from the transformation

$Y_j=1-e^{-e^{{\mu }_j+{\sigma }_j\ast {sinh}^{-1}\left(\frac{{\alpha }_j{Z}_j}{2}\right)}}\mathrm{with}{Z}_j\sim \mathrm{N}(0,1),j=1,2.$

The bivariate vector $(Y_1,Y_2)$ is said to have a bivariate USHN distribution with parameters ${\alpha }_1$, ${\mu }_1$, ${\sigma }_1$, ${\alpha }_2$, ${\mu }_2$, ${\sigma }_2$, and $\rho$ if the joint cdf of $(Y_1,Y_2)$ is given by

(3)$F_{Y_1,Y_2}(y_1,y_2)={\mathsf{\Phi }}_2({\omega }_{y_1},{\omega }_{y_2};\rho ),$

where

${\omega }_{y_j}=\omega (y_j;{\alpha }_j,{\mu }_j,{\sigma }_j)=\frac{2}{{\alpha }_j}sinh\left(\frac{log(-log(1-y_j))-{\mu }_j}{{\sigma }_j}\right).$

The joint pdf of $(Y_1,Y_2)$ which follows from (3) can be expressed as

(4)$f_{Y_1,Y_2}(y_1,y_2)={\varphi }_2({\omega }_{y_1},{\omega }_{y_2};\rho ){\mathsf{\Omega }}_{y_1}{\mathsf{\Omega }}_{y_2},$

where

${\mathsf{\Omega }}_{y_j}=\mathsf{\Omega }(y_j;{\alpha }_j,{\mu }_j,{\sigma }_j)={\displaystyle \frac{\frac{2}{{\alpha }_j}cosh\left(\frac{log(-log(1-y_j))-{\mu }_j}{{\sigma }_j}\right)}{{\sigma }_j(1-y_j)log{(1-y_j)}^{-1}}},j=1,2.$

The notation used is $(Y_1,Y_2)\sim \mathrm{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$, where $\rho$ is the correlation coefficient. From the previous definition, it follows that

$(Z_1,Z_2)\sim \mathrm{BN}(0,0,1,1,\rho ),$

where

$Z_j=\frac{2}{{\alpha }_j}sinh\left(\frac{log(-log(1-Y_j))-{\mu }_j}{{\sigma }_j}\right),j=1,2.$

We have the following theorems.

Theorem 3.

Let $(Y_1,Y_2)\sim \mathit{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$. Hence:

• (i) 

  $Y_j\sim USHN({\alpha }_j,{\mu }_j,{\sigma }_j),$ for $j=1,2.$

• (ii) 

  The cdf of $Y_1$, given $Y_2=y_2$, is

  $F_{Y_1\mid Y_2=y_2}(y_1\mid Y_2=y_2)=\mathsf{\Phi }\left(\frac{{\omega }_{y_1}-\rho {\omega }_{y_2}}{\sqrt{1-{\rho }^2}}\right).$

• (iii) 

  The pdf of $Y_1$, given $Y_2=y_2$, is

  $f_{Y_1\mid Y_2=y_2}(y_1\mid Y_2=y_2)=\frac{{\mathsf{\Omega }}_{y_1}}{\sqrt{2\pi }\sqrt{1-{\rho }^2}}exp\left(-\frac{1}{2}{\left(\frac{{\omega }_{y_1}-\rho {\omega }_{y_2}}{\sqrt{1-{\rho }^2}}\right)}^2\right).$

• (iv) 

  The joint survival function of $(Y_1,Y_2)$ is

  $S_{Y_1,Y_2}(y_1,y_2)=1-\mathsf{\Phi }\left({\omega }_{y_1}\right)-\mathsf{\Phi }\left({\omega }_{y_2}\right)+{\mathsf{\Phi }}_2({\omega }_{y_1},{\omega }_{y_2};\rho ).$

Theorem 4.

If $(Y_1,Y_2)\sim \mathit{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$, then the conditional hazard function of $Y_1$, given $Y_2=y_2$, is a non-decreasing function of $y_1,$ for all values of $y_2$ and ρ when ${\mu }_1>0$ and ${\sigma }_1<1.$

Marshall [38] defined the bivariate hazard rate of $Y_1$ and $Y_2$ as

$h(y_1,y_2)=\left(-\frac{\partial }{\partial y_1},-\frac{\partial }{\partial y_2}\right)logP(Y_1>y_1,Y_2>y_2)=(h_1(y_1,y_2),h_2(y_1,y_2)).$

Theorem 5.

If $(Y_1,Y_2)\sim \mathit{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$, then:

• (i) 

  For $y_2$ fixed and ${\sigma }_1<1$, $h_1(y_1,y_2)$ is a non-decreasing function of $y_1.$

• (ii) 

  For $y_1$ fixed and ${\sigma }_2<1$, $h_2(y_1,y_2)$ is a non-decreasing function of $y_2.$

Remark 2.

We can generate samples from the $\mathit{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$ distribution according to the algorithm described below:

• (i) 

  For $j=1,2$, generate $Z_j$ according to the algorithm given by Kundu [23];

• (ii) 

  For each value of $Z_j$, compute

  $Y_j=1-e^{-e^{{\mu }_j+{\sigma }_j\ast {sinh}^{-1}\left(\frac{{\alpha }_j{Z}_j}{2}\right)}},j=1,2.$

  Hence, $(Y_1,Y_2)\sim \mathit{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$.

5. Bivariate Unit SHN Regression Model

Suppose that ${\mathbf{Y}}_1=(Y_{11},Y_{12}),\cdots ,{\mathbf{Y}}_n=(Y_{n1},Y_{n2})$ are n independent bivariate vectors of proportions, where for each bivariate vector we assume that

${\mathbf{Y}}_i={(Y_{i1},Y_{i2})}^{\top }\sim \mathrm{BVUSHN}({\alpha }_1,{\mathbf{X}}_i{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{X}}_i{\mathbf{\beta }}_2,{\sigma }_2,\rho ),i=1,\dots ,n,$

where ${\mathbf{X}}_i=(x_{i1},x_{i2},\dots ,x_{ip})$ is a p-dimensional vector with values of the explanatory variables associated with the j-th observable response ${\mathbf{y}}_j$, such that ${\mathbf{Y}}_{ij}\sim \mathrm{USHN}({\alpha }_j,{\mathbf{X}}_i{\mathbf{\beta }}_j,{\sigma }_j)$. In addition, ${\mathbf{\beta }}_j$ is a vector of size p of unknown parameters for $j=1,2$. In other words, we assume that, for each $i=1,2,\dots ,n$ and $j=1,2$, the location parameter of the response variable $Y_{ij}$ satisfies the functional relationship

$log(-log(1-{\mu }_{ij}))={\mathbf{X}}_{ij}{\mathbf{\beta }}_j,$

which corresponds to the log-log complement link function.

Let $\mathbf{Y}=vec({\mathbf{Y}}_1,{\mathbf{Y}}_2)$ be of size $2n,$ $\mathbf{X}={({\mathbf{X}}_1^{\top },{\mathbf{X}}_2^{\top },\cdots ,{\mathbf{X}}_n^{\top })}^{\top }$ and ${\mathbf{X}}^{*}=diag(\mathbf{X},\mathbf{X})$ be of size $2n\times 2p$, where $vec(·)$ is the vec operator, and $\mathbf{\beta }={({\mathbf{\beta }}_1^{\top },{\mathbf{\beta }}_2^{\top })}^{\top }=({\beta }_1,{\beta }_2,\dots ,{\beta }_{2p})$ be a $2p$-vector of unknown parameters.

5.1. Parameter Estimation

The estimation of the parameter vector $\mathbf{\theta }={({\alpha }_1,{\alpha }_2,{\sigma }_1,{\sigma }_2,\mathbf{\beta },\rho )}^{\top }$ of size $5+2p$ can be carried out by maximizing the log-likelihood function

$\ell \left(\mathbf{\theta }\right)=\sum _{i=1}^n{\ell }_i\left(\mathbf{\theta }\right),$

where

$\begin{array}{cc}\hfill {\ell }_i\left(\mathbf{\theta }\right)& =log\left({\varphi }_2(\omega (y_{i1};{\alpha }_1,{\mathbf{X}}_{i1}\mathbf{\beta },{\sigma }_1),\omega (y_{i2};{\alpha }_2,{\mathbf{X}}_{i2}\mathbf{\beta },{\sigma }_2);\rho )\right)\hfill \\ \hfill & +\sum _{j=1}^{2}log\left(\mathsf{\Omega }(y_{ij};{\alpha }_j,{\mathbf{X}}_{ij}\mathbf{\beta },{\sigma }_j)\right),\hfill \end{array}$

with

$\omega (Y_{ij};{\alpha }_j,{\mathbf{X}}_{ij}\mathbf{\beta },{\sigma }_j)=\frac{2}{{\alpha }_j}sinh\left({\displaystyle \frac{log(-log(1-Y_{ij}))-{\mathbf{X}}_{ij}\mathbf{\beta }}{{\sigma }_j}}\right)\sim \mathrm{N}(0,1).$

We have that

(5)$V_{ij}=\frac{log(-log(1-Y_{ij}))-{\mathbf{X}}_{ij}\mathbf{\beta }}{{\sigma }_j}\sim \mathrm{SHN}({\alpha }_j,0,1),$

for $i=1,2,\dots ,n$ and $j=1,2$. Thus,

$\left(e^{V_{i1}}-e^{-V_{i1}},e^{V_{i2}}-e^{-V_{i2}}\right)\sim {\mathrm{N}}_2\left(\left(\genfrac{}{}{0pt}{}{0}{0}\right),\left(\begin{array}{cc}{\alpha }_1^2& {\alpha }_1{\alpha }_2\rho \\ {\alpha }_1{\alpha }_2\rho & {\alpha }_2^2\end{array}\right)\right).$

For ${\sigma }_1$, ${\sigma }_2$, and $\mathbf{\beta }$ fixed, with 

$g(y_{ij};{\mathbf{X}}_{ij}\mathbf{\beta },{\sigma }_j)=e^{v_{ij}}-e^{-v_{ij}},$

it follows that

${\widehat{\alpha }}_j(\mathbf{\beta },{\sigma }_j)={\left(\frac{1}{n}\sum _{i=1}^n{g}^2(y_{ij};{\mathbf{X}}_{ij}\mathbf{\beta },{\sigma }_j)\right)}^{1/2},$

and

$\widehat{\rho }(\mathbf{\beta },{\sigma }_1,{\sigma }_2)=\frac{{\sum }_{i=1}^{n}g(y_{i1};{\mathbf{X}}_{i1}\mathbf{\beta },{\sigma }_1)g(y_{i2};{\mathbf{X}}_{i2}\mathbf{\beta },{\sigma }_2)}{\sqrt{{\sum }_{i=1}^n{g}^2(y_{i1};{\mathbf{X}}_{i1}\mathbf{\beta },{\sigma }_1)}\sqrt{{\sum }_{i=1}^n{g}^2(y_{i2};{\mathbf{X}}_{i2}\mathbf{\beta },{\sigma }_2)}}.$

Thus, to estimate ${\sigma }_1$, ${\sigma }_2$, and $\mathbf{\beta }$, we use the profiled likelihood approach; that is, the estimates of these parameters are obtained by maximizing

${\ell }_p(\mathbf{\beta },{\sigma }_1,{\sigma }_2)=\ell ({\widehat{\alpha }}_1(\mathbf{\beta },{\sigma }_1),{\sigma }_1,\mathbf{\beta },{\widehat{\alpha }}_2(\mathbf{\beta },{\sigma }_2),{\sigma }_2,\widehat{\rho }(\mathbf{\beta },{\sigma }_1,{\sigma }_2)),$

regarding ${\sigma }_1$, ${\sigma }_2$, and $\mathbf{\beta }$. This maximization procedure requires iterative numerical methods such as Newton–Raphson or quasi-Newton.

Denoting the estimators of $\mathbf{\beta },$${\sigma }_1$, and ${\sigma }_2$ by $\widehat{\mathbf{\beta }},$${\widehat{\sigma }}_1$, and ${\widehat{\sigma }}_2$, respectively, then the estimators of ${\alpha }_1,$${\alpha }_2$,and $\rho$ are given by ${\widehat{\alpha }}_1={\widehat{\alpha }}_1(\widehat{\mathbf{\beta }},{\widehat{\sigma }}_1),$${\widehat{\alpha }}_2={\widehat{\alpha }}_2(\widehat{\mathbf{\beta }},{\widehat{\sigma }}_2)$, and $\widehat{\rho }=\widehat{\rho }(\widehat{\mathbf{\beta }},{\widehat{\sigma }}_1,{\widehat{\sigma }}_2)$, respectively. The elements of the Fisher information matrix can be obtained by calculating the second derivatives of the log-likelihood function regarding the parameters of the model and thus finding the expectation numerically.

5.2. Two-Step Estimation

From (5), it follows that

$log(-log(1-Y_{ij}))={\mathbf{X}}_{ij}\mathbf{\beta }+{\sigma }_j{V}_{ij}\sim \mathrm{SHN}({\alpha }_j,{\mathbf{X}}_{ij}\mathbf{\beta },{\sigma }_j),$

for $i=1,2,\dots ,n$ and $j=1,2$; that is, the bivariate USHN regression model can be seen as

$W_{ij}=log(-log(1-Y_{ij}))={\mathbf{X}}_{ij}\mathbf{\beta }+{ϵ}_{ij}={\mu }_{ij}+{ϵ}_{ij},$

where ${ϵ}_{ij}\sim \mathrm{SHN}({\alpha }_j,0,{\sigma }_j)$ for $i=1,2,\dots ,n$ and $j=1,2$. Since $W_{ij}\sim \mathrm{SHN}({\alpha }_j,{\mu }_{ij},{\sigma }_j)$ for $i=1,2,\dots ,n$ and $j=1,2$, we can use the two-step method presented in Kundu [23] to estimate the parameters of the bivariate USHN regression model.

The two-step estimation method was provided by Joe [39] and it is an interesting procedure because it reduces the dimension of the parameter space in each of its stages. In our case, the estimation of the parameters will be determined in two independent stages, where in the first stage, the estimates of the marginals of each variable under study will be obtained, that is, $2+p$ will be estimated in the marginal distribution and, in the second stage, the association parameter $\rho$ will be estimated. The ML estimators obtained by the two-stage procedure are consistent and the variance-covariance matrix is obtained using the methodology presented in Joe [39].

Expressing the joint cdf of $W_1$ and $W_2$ in terms of the normal copula, we have that

$\begin{array}{ccc}\hfill F_{W_1{W}_2}(w_1,w_2)& =& {\mathsf{\Phi }}_2(a(w_1;{\alpha }_1,{\mu }_{W_1},{\sigma }_1),a(w_2;{\alpha }_2,{\mu }_{W_2},{\sigma }_2);\rho )\hfill \\ & =& C_G(F_{W_1}\left(w_1\right),F_{W_2}\left(w_2\right);\rho )=C_G(u,v;\rho ),\hfill \end{array}$

where $C_G(·,·;\rho )$ denotes the Gaussian copula. In addition, $u,v\sim U(0,1)$ and $F_{W_j}(·)$ for $j=1,2$ is the cdf of $W_j$. Then, let $c_G(u,v;\rho )={\partial }^2{C}_G(u,v;\rho )/\partial u\partial v$ be the pdf of the Gaussian copula. It follows that the log-likelihood function can be expressed as

$\begin{array}{ccc}\hfill \ell \left(\theta \right)& =& \sum _{i=1}^{n}log{c}_G\left(F(w_{i1};{\alpha }_1,{\mu }_{i1},{\sigma }_1),F(w_{i2};{\alpha }_2,{\mu }_{i2},{\sigma }_2);\rho \right)\hfill \\ & +& \sum _{i=1}^{n}logf(w_{i1};{\alpha }_1,{\mu }_{i1},{\sigma }_1)+\sum _{i=1}^{n}logf(w_{i2};{\alpha }_2,{\mu }_{i2},{\sigma }_2),\hfill \end{array}$

where $f(w_{ij};{\alpha }_j,{\mu }_{ij},{\sigma }_j)$ is the pdf of $W_{ij}$ for $i=1,2,\dots ,n$ and $j=1,2$ with ${\mu }_{ij}={\mathbf{X}}_{ij}{\mathbf{\beta }}_j.$ Then, to implement the algorithm in two stages, we initially estimate ${\alpha }_1,$${\sigma }_1$, and ${\mathbf{\beta }}_1$ by maximizing

$\sum _{i=1}^{n}logf(w_{i1};{\alpha }_1,{\mu }_{i1},{\sigma }_1)=\sum _{i=1}^{n}log\varphi \left(b(w_{i1};{\alpha }_1,{\mathbf{\beta }}_1,{\sigma }_1)\right)+\sum _{i=1}^{n}logB(w_{i1};{\alpha }_1,{\mathbf{\beta }}_1,{\sigma }_1).$

Setting ${\sigma }_1$ and ${\mathbf{\beta }}_1$, the maximization of the previous function with respect to ${\alpha }_1$ can be obtained as

${\overline{\alpha }}_1\left({\mathbf{\beta }}_1{\sigma }_1\right)={\left(\frac{1}{n}\sum _{i=1}^n{g}^2(w_{i1};{\mathbf{\beta }}_1,{\sigma }_1)\right)}^{1/2},$

where $g(·)$ was defined before. Then, the estimates of ${\mathbf{\beta }}_1$ and ${\sigma }_1$, say, ${\overline{\mathbf{\beta }}}_1$ and ${\overline{\sigma }}_1$, can be obtained by maximizing

$\sum _{i=1}^{n}log\varphi \left(b(w_{i1};{\overline{\alpha }}_1({\mathbf{\beta }}_1,{\sigma }_1),{\mathbf{\beta }}_1,{\sigma }_1)\right)+\sum _{i=1}^{n}logB(w_{i1};{\overline{\alpha }}_1({\mathbf{\beta }}_1,{\sigma }_1),{\mathbf{\beta }}_1,{\sigma }_1).$

In a similar way, we estimate ${\alpha }_2,$${\sigma }_2$, and ${\mathbf{\beta }}_2$ by maximizing

$\sum _{i=1}^{n}logf(w_{i2};{\alpha }_2,{\mu }_{i2},{\sigma }_2)=\sum _{i=1}^{n}log\varphi \left(b(w_{i2};{\alpha }_2,{\mathbf{\beta }}_2,{\sigma }_2)\right)+\sum _{i=1}^{n}logB(w_{i2};{\alpha }_2,{\mathbf{\beta }}_2,{\sigma }_2).$

Finally, to estimate $\rho$, we maximize

$h\left(\rho \right)=\sum _{i=1}^{n}log{c}_G\left(u_{i1},u_{i2};\rho \right),$

where $u_{i1}=F_{W_1}(w_{i1},{\overline{\alpha }}_1,{\overline{\mu }}_{i1},{\overline{\sigma }}_1)$ and $u_{i2}=F_{W_2}(w_{i2},{\overline{\alpha }}_2,{\overline{\mu }}_{i2},{\overline{\sigma }}_2)$, where ${\overline{\alpha }}_1,{\overline{\mathbf{\beta }}}_1,{\overline{\sigma }}_1,$${\overline{\alpha }}_2,{\overline{\mathbf{\beta }}}_2,$ and ${\overline{\sigma }}_2$ are the ML estimates of ${\alpha }_1,$${\mathbf{\beta }}_1,$${\sigma }_1$, ${\alpha }_2,$${\mathbf{\beta }}_2$, and ${\sigma }_2$ obtained in step 1. After intensive calculations, it is found that the ML estimator of $\rho$ is given by

$\overline{\rho }=\frac{{\sum }_{i=1}^n{z}_{i1}{z}_{i2}}{\sqrt{{\sum }_{i=1}^n{z}_{i1}^2}\sqrt{{\sum }_{i=1}^n{z}_{i2}^2}},$

where $z_{i1}={\mathsf{\Phi }}^{-1}\left(u_{i1}\right)=b(w_{i1},{\overline{\alpha }}_1,{\overline{\mu }}_{i1},{\overline{\sigma }}_1)$ and $z_{i2}={\mathsf{\Phi }}^{-1}\left(u_{i2}\right)=b(w_{i2},{\overline{\alpha }}_2,{\overline{\mu }}_{i2},{\overline{\sigma }}_2).$

6. Simulation Study

This section presents the results of a simulation study carried out with the aim of studying the asymptotic behavior of the parameter estimators in the BVSHN regression model as well as in the proposed method. We considered the regression model given by

(6)${\mathbf{Y}}_i={(Y_{i1},Y_{i2})}^{\top }\sim \mathrm{BVUSHN}({\alpha }_1,{\mathbf{X}}_i{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{X}}_i{\mathbf{\beta }}_2,{\sigma }_2,\rho ),\mathrm{for}i=1,\dots ,n,$

where ${\mathbf{\beta }}_1^{\top }={({\beta }_{10},{\beta }_{11})}^{\top }$, ${\mathbf{\beta }}_2^{\top }={({\beta }_{20},{\beta }_{21})}^{\top }$, and ${\mathbf{X}}_i=(1,x_i)$. In this study, there were 5000 Monte Carlo samples of sizes $n=30,50,100,200,500$, and 1000 from the model given in (6). The values of the covariate x were generated independently from a uniform distribution $\mathrm{U}(0,1)$, and the following three scenarios for the parameter vector $\mathbf{\theta }={({\alpha }_1,{\mathbf{\beta }}_1^{\top },{\sigma }_1,{\alpha }_2,{\mathbf{\beta }}_2^{\top },{\sigma }_2,\rho )}^{\top }$ were considered:

• Scenario 1: ${\mathbf{\beta }}_1={(0.25,0.75)}^{\top }$, ${\mathbf{\beta }}_2={(-0.50,0.25)}^{\top }$, ${\sigma }_1={\sigma }_2=1.0$, ${\alpha }_1=1.50$, ${\alpha }_2=2.50$, and $\rho =0.25,0.50,0.75$.

• Scenario 2: ${\mathbf{\beta }}_1={(0.25,0.75)}^{\top }$, ${\mathbf{\beta }}_2={(-0.50,0.25)}^{\top }$, ${\sigma }_1={\sigma }_2=1.0$, ${\alpha }_1=2.00$, ${\alpha }_2=0.50$, and $\rho =-0.25,-0.75$.

• Scenario 3: ${\mathbf{\beta }}_1={(1.50,-0.50)}^{\top }$, ${\mathbf{\beta }}_2={(-1.25,0.75)}^{\top }$, ${\alpha }_1=0.75$, ${\alpha }_2=1.0$, ${\sigma }_1=0.50$, ${\sigma }_2=1.50$, and $\rho =-0.50$; and ${\mathbf{\beta }}_1={(1.50,-0.50)}^{\top }$, ${\mathbf{\beta }}_2={(-1.25,0.25)}^{\top }$, ${\alpha }_1=1.0$, ${\alpha }_2=3.0$, ${\sigma }_1=1.75$, ${\sigma }_2=1.50$, and $\rho =0.50$.

Maximum likelihood estimates were obtained using the maxLik function of the statistical software R Development Core Team [40]. For the optimization of the log-likelihood function, iterative methods based on the Newton–Raphson algorithm were considered. To evaluate the performance of the estimators, the absolute value of the bias (AVB) and the root of the mean square error (RMSE) were considered, which are given by:

$\mathrm{AVB}\left({\theta }_i\right)=\frac{1}{5000}\left|\sum _{k=1}^{5000}\left({\widehat{\theta }}_j^{\left(k\right)}-{\theta }_j\right)\right|;\mathrm{RMSE}\left({\theta }_i\right)=\sqrt{\frac{1}{5000}\sum _{k=1}^{5000}{\left({\widehat{\theta }}_j^{\left(k\right)}-{\theta }_j\right)}^2}$

respectively, where ${\widehat{\theta }}_j^{\left(k\right)}$ is the estimator of ${\theta }_j$ for the jth sample, for ${\theta }_i\in {({\alpha }_1,{\mathbf{\beta }}_1^{\top },{\sigma }_1,{\alpha }_2,{\mathbf{\beta }}_2^{\top },{\sigma }_2,\rho )}^{\top }$.

To generate the random sample within the simulation process, we implemented Algorithm 1 based on the distribution function method. We initially defined:

• n: the number of rows of the $n\times 3\_\mathrm{matrix}$.

• $(Y_1,Y_2):$ a random vector with a $\mathrm{BVSHN}({\alpha }_1,{\mathbf{X}}_i{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{X}}_i{\mathbf{\beta }}_2,{\sigma }_2,\rho )$ distribution.

• $U_1:$ a random variable with a standard normal distribution, $\mathrm{N}(0,1)$ in short.

• $U_2:$ a random variable with a standard normal distribution, $\mathrm{N}(0,1)$ in short.

• $X:$ a random variable with a uniform distribution, $\mathrm{U}(0,1)$ in short.

• $Z_1:$ an auxiliary random variable.

• $Z_2:$ an auxiliary random variable.

Algorithm 1: Algorithm to generate a random sample from the $\mathrm{BVSHN}({\alpha }_1,{\mathbf{X}}_i{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{X}}_i{\mathbf{\beta }}_2,{\sigma }_2,\rho )$ regression model

For $j=1,2$, we generated $Z_j$ according to the algorithm given by Kundu [23]. For each value of $Z_j$, we computed

$Y_j=1-e^{-e^{{\mu }_j+{\sigma }_j\ast {sinh}^{-1}\left(\frac{{\alpha }_j{Z}_j}{2}\right)}},j=1,2.$

Hence, $(Y_1,Y_2)\sim \mathrm{BVUSHN}({\alpha }_1,{\mu }_1,{\sigma }_1,{\alpha }_2,{\mu }_2,{\sigma }_2,\rho )$.

The results of the simulation study are presented in Table A1, Table A2 and Table A3 given in the Appendix A. The Figure 1 shows the behavior of the AVB and the RMSE for each of the parameter estimators for scenario 1, while the Figure 2 does the same for scenario 2. It should be clarified that, for scenario 3, the respective graphs are not presented for reasons of space; however, these graphs can be constructed from Table A3.

In the Figure 1 (and in Table A1), it can be seen that, in general, when the sample size increases, both the AVB and the RMSE of the MLE decrease and tend to zero. This behavior can be seen for all estimators of the components of the parameter vector $\theta$. In particular, for small sample sizes ($n=30$ or 60), the estimators ${\alpha }_1$ and ${\alpha }_2$ present relatively high AVB and RMSE, compared to the other estimators; however, these quantities decrease considerably when $n\ge 100$. In relation to the estimators of the coefficients, ${\beta }_{0j}$ and ${\beta }_{1j}$ for $j=1,2$; these present low AVB for small sample sizes and approach zero as n increases. Something similar occurs with the RMSE of these estimators. Finally, the $\widehat{\rho }$ estimator of $\rho$ shows a stable behavior, with considerably low AVB and RMSE.

The above comments are generally equally valid for scenarios 2 and 3 and can be seen in Table A2 and Table A3, respectively. In particular, for simulation scenario 2, it can be better appreciated in the Figure 2.

7. Real Data Applications

Proportion, rate, or indices data arise in different areas of knowledge such as economics, agriculture, medicine, engineering, and social sciences, among others. Some practical situations in these fields include, for example: the Human Development Index or the illiteracy rate in a given country, the proportion of diseased teeth (periodontal disease) in patients, the mortality rate due to traffic accidents, the proportion of deaths caused by smoking or other exposure factors, the proportion of non-compliant items on a production line, and the percentage of a family’s salary spent on entertainment, among others. In this section, the usefulness of the proposed models is presented. The BVUBS distribution and BVUBS regression models are fitted to real datasets.

7.1. First Application

In the first illustration, we consider a dataset obtained from [41]. The data are related to the time in minutes $X_1$ elapsed until a soccer team scores a first kick goal, while $X_2$ records the time of the first goal scored by either of the two teams in dispute, either through a penalty kick, a free kick, or any other direct kick. To confirm that the observed data belong to the interval $\left(0.1\right)$, all the measurements were divided by the total time that a soccer game lasts, that is, 90 min. The interest in this case is therefore to model the proportion of time that any team and the home team score the first kick goal.

To fit these bivariate data, we consider the bivariate distributions: the BVBS, bivariate log-normal (BVLN), and BVUBS distributions. To compare these bivariate distributions, we consider the AIC and BIC, given by

$AIC=-2\ell (\widehat{\mathit{\theta }})+2k,BIC=-2\ell (\widehat{\mathit{\theta }})+klog(n),$

where k is the number of parameters and $\ell (\widehat{\mathit{\theta }})$ means the maximum value of the log-likelihood function. For fitting the bivariate distributions, we use the maxLik function of R Development Core Team [40] like in Section 6.

The ML estimates (standard errors between parentheses) are presented in Table 1.

Therefore, according to the AIC and BIC values, the BVUBS distribution outperforms the other bivariate distributions and so should be preferable to fit these data. Contour plots for the BVLN, BVBS, and BVUBS distributions are presented in Figure 3.

7.2. Second Application

The interest in this application is to model the Human Development Index, denoted by $Y_1$, and illiteracy rate, $Y_2$, as functions of the natural logarithm of life expectancy measured in years ($log\left(LE\right)$, which is denoted by $X_1$) and the percentage of people with high poverty level (HPL, in this case denoted by $X_2$) by using the BVUSHN regression model. This dataset was taken from 195 districts in Peru and more details about this dataset can be found at http://www.pe.undp.org.

The ML estimates (standard errors between parentheses) are given in Table 2.

From the previous estimates, we compute

$e_{ij}=\frac{2}{{\widehat{\alpha }}_j}sinh\left(\frac{log(-log(1-y_{ij}))-({\widehat{\beta }}_{j0}+{\widehat{\beta }}_{j1}{X}_{i1}+{\widehat{\beta }}_{j2}{X}_{i2})}{{\widehat{\sigma }}_j}\right).$

Therefore, we have that $(e_{i1},e_{i2})\sim \mathrm{BN}(0,0,1,1,\widehat{\rho }).$ The Shapiro–Wilk and Doornik–Hansen bivariate normality tests for the vector $(e_{i1},e_{i2})$ return the values 0.97669 (p-value $<0.001$) and $21.41961$ ( p-value $<0.001$), respectively; that is, the assumption of multivariate normality is rejected. The contour plot for the BN distribution is presented in Figure 4a. Calculating the Mahalanobis distance, we find that the observations {11, 13, 14, 27, 142, 143, 192, 193, 195} are considered outliers. Eliminating this set of observations, we estimate the parameters of the model, obtaining now the ML estimates (standard errors between parentheses) presented in Table 3.

The Shapiro–Wilk and Doornik–Hansen bivariate normality tests for the vector $(e_{i1},e_{i2})$ return the values 0.9914 (p-$\mathrm{value}>0.05$) and 3.9346 (p-$\mathrm{value}>0.05$), respectively; that is, the assumption of multivariate normality is not rejected.

8. Concluding Remarks

The statistical modeling of proportions, rates, or indices data is a very common activity in many areas of knowledge, such as medicine, economy, engineering, or social sciences. Although there are some proposals for modeling this type of data, in the statistical literature there are few alternatives capable of jointly modeling two or more variables of this type, even less so with the possibility of considering covariates to explain the variability of the responses.

In this work, from a transformation applied to a vector with a bivariate unit-Birnbaum–Saunders distribution, we have obtained a new bivariate distribution which is called the bivariate unit-sinh-normal (BVSHN) distribution. The new distribution is an absolutely continuous distribution and is useful for modeling data as described previously, that is, responses over the region $(0,1)\times (0,1)$ jointly. Based on the introduced proposal, the extension to the case of the regression model for bounded responses in $(0,1)\times (0,1)$ was proposed, which was called the BVUSHN regression model. In particular, for the introduced distribution, the joint probability density function and joint cumulative distribution were specified. Conditional distributions and the joint survival function were also presented. The estimation of the model parameters was carried out from a classical approach by using the maximum likelihood method together with the two-step method proposed by Joe [39]. A Monte Carlo simulation study was also carried out to evaluate the benefits and limitations of the introduced proposals, yielding good asymptotic results for the parameter estimators.

The BVUSHN distribution and the BVUSHN regression model showed great flexibility in fitting data to $(0,1)$, which is frequently encountered in many practical scenarios. The results obtained in the two applications showed better results than other existing methodologies in the literature, and therefore we conclude that our proposals are viable alternatives in the field of distribution theory and regression models.

Future works based on the results of this proposal could contemplate the extension to the case of more than two response variables and make inferences in the models considering a Bayesian perspective.

Author Contributions

Conceptualization, G.M.-F.; Data curation, A.J.L. and G.M.-A.; Formal analysis, G.M.-F., A.J.L. and G.M.-A.; Funding acquisition, R.T.-F.; Investigation, G.M.-F., A.J.L. and G.M.-A.; Methodology, G.M.-F., A.J.L. and G.M.-A.; Project administration, G.M.-F.; Resources, G.M.-A. and R.T.-F.; Software, G.M.-F. and A.J.L.; Supervision, G.M.-F., A.J.L. and G.M.-A.; Validation, A.J.L. and G.M.-A.; Visualization, G.M.-F. and A.J.L.; Writing—original draft, G.M.-F., A.J.L., G.M.-A. and R.T.-F.; Writing—review & editing, G.M.-F., A.J.L. and R.T.-F. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Details about the available data are given in Section 7.

Acknowledgments

G. Martínez-Flórez and R. Tovar-Falón acknowledge the support given by Universidad de Córdoba, Montería, Colombia.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 1.

Figure 1. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 1.

Figure 1. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 1.

Figure 2. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 2.

Figure 2. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 2.

Figure 2. Asymptotic behavior of the MLE for the BVSHN regression model for scenario 2.

Figure 3. Contour plots of bivariate distributions: (a) BVLN, (b) BVBS, and (c) BVUBS.

Figure 4. Contour plot of the bivariate vector [Forumla omitted. See PDF.]: (a) with outliers and (b) without outliers.

Table 1

Estimation of the parameters with their standard errors for the BVBS, BVLN, and BVUBS distributions.

Estimator	BVBS	BVLN	BVUBS
${\widehat{\alpha }}_1$	$0.8305 \left(0.0971\right)$		$0.8236 \left(0.0962\right)$
${\widehat{\alpha }}_2$	$1.0659 \left(0.1253\right)$		$1.1160 \left(0.1301\right)$
${\widehat{\beta }}_1$	$0.3204 \left(0.0402\right)$		$0.6840 \left(0.0846\right)$
${\widehat{\beta }}_2$	$0.2209 \left(0.0339\right)$		$0.7786 \left(0.1218\right)$
${\widehat{\mu }}_1$		$-0.9637 \left(0.1154\right)$
${\widehat{\mu }}_2$		$-1.3233 \left(0.1556\right)$
${\widehat{\sigma }}_1$		$0.7024 \left(0.0816\right)$
${\widehat{\sigma }}_2$		$0.9469 \left(0.1100\right)$
$\widehat{\rho }$	$0.5558 \left(0.1152\right)$	$0.5066 \left(0.1221\right)$	$0.5966 \left(0.1065\right)$
AIC	$15.7486$	$9.6148$	$-11.8823$
BIC	$29.4732$	$25.3395$	$-3.8277$

Table 2

Estimation of the parameters with their standard errors for the BVUSHN regression model.

Estimator	BVUBS
${\widehat{\alpha }}_1$	$0.0173 \left(0.0051\right)$
${\widehat{\alpha }}_2$	$0.1786 \left(0.0952\right)$
${\widehat{\beta }}_{10}$	$-3.9374 \left(0.8521\right)$
${\widehat{\beta }}_{11}$	$-0.5342 \left(0.0420\right)$
${\widehat{\beta }}_{12}$	$0.9140 \left(0.1979\right)$
${\widehat{\beta }}_{20}$	$6.9501 \left(6.4588\right)$
${\widehat{\beta }}_{21}$	$3.5038 \left(0.3266\right)$
${\widehat{\beta }}_{22}$	$-2.3048 \left(1.5009\right)$
${\widehat{\sigma }}_1$	$6.7323 \left(2.0491\right)$
${\widehat{\sigma }}_2$	$5.6380 \left(2.9657\right)$
$\widehat{\rho }$	$-0.8435 \left(0.0251\right)$

Table 3

Estimation of the parameters with their standard errors for the BVUSHN regression model without outliers.

Estimator	BVUBS
${\widehat{\alpha }}_1$	$0.0145 \left(0.0038\right)$
${\widehat{\alpha }}_2$	$0.1412 \left(0.0554\right)$
${\widehat{\beta }}_{10}$	$-3.8892 \left(0.6937\right)$
${\widehat{\beta }}_{11}$	$-0.5408 \left(0.0339\right)$
${\widehat{\beta }}_{12}$	$0.9046 \left(0.1611\right)$
${\widehat{\beta }}_{20}$	$6.9470 \left(5.5195\right)$
${\widehat{\beta }}_{21}$	$3.6462 \left(0.2752\right)$
${\widehat{\beta }}_{22}$	$-2.3220 \left(1.2821\right)$
${\widehat{\sigma }}_1$	$6.8814 \left(1.9256\right)$
${\widehat{\sigma }}_2$	$6.0481 \left(2.3610\right)$
$\widehat{\rho }$	$-0.8172 \left(0.0274\right)$

Appendix A. Tables of Simulations

This section presents the results of the simulation study described in Section 6.

Table A1

Asymptotic behavior of the MLE for the BVSHN regression model with parameter vector $\theta ={({\alpha }_1,{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{\beta }}_2,{\sigma }_2,\rho )}^{\top }$. Simulation scenario 1.

		$\mathit{\rho }=0.25$			$\mathit{\rho }=0.50$			$\mathit{\rho }=0.75$
$\mathbf{n}$	Estimates	Mean	AVB	RMSE	Mean	AVB	RMSE	Mean	AVB	RMSE
	${\widehat{\alpha }}_1$	2.8547	1.3547	5.0699	2.9825	1.4825	6.5773	2.7954	1.2954	4.9636
	${\widehat{\beta }}_{10}$	0.2643	0.0143	0.0588	0.2651	0.0151	0.0531	0.2631	0.0131	0.0538
	${\widehat{\beta }}_{11}$	0.7113	0.0387	0.1753	0.7072	0.0428	0.1610	0.7146	0.0354	0.1597
	${\widehat{\sigma }}_1$	0.8288	0.1712	0.3139	0.8469	0.1531	0.3678	0.8306	0.1694	0.2754
30	${\widehat{\alpha }}_2$	4.1517	1.6517	9.2441	4.1888	1.6888	9.2486	4.2104	1.7104	10.1232
	${\widehat{\beta }}_{20}$	−0.4963	0.0037	0.0900	−0.4829	0.0171	0.0854	−0.4901	0.0099	0.0924
	${\widehat{\beta }}_{21}$	0.2350	0.0150	0.2729	0.2082	0.0418	0.2628	0.2169	0.0331	0.2786
	${\widehat{\sigma }}_2$	0.8754	0.1246	0.1746	0.8841	0.1159	0.2115	0.8754	0.1246	0.2247
	$\widehat{\rho }$	0.2457	0.0043	0.0334	0.4944	0.0056	0.0222	0.7461	0.0039	0.0088
	${\widehat{\alpha }}_1$	2.1775	0.6775	1.7276	2.2095	0.7095	1.7357	2.0779	0.5779	1.4032
	${\widehat{\beta }}_{10}$	0.2607	0.0107	0.0288	0.2644	0.0144	0.0314	0.2578	0.0078	0.0264
	${\widehat{\beta }}_{11}$	0.7120	0.0380	0.0833	0.7082	0.0418	0.0948	0.7198	0.0302	0.0793
	${\widehat{\sigma }}_1$	0.9144	0.0856	0.2469	0.9037	0.0963	0.3308	0.9426	0.0574	0.2417
50	${\widehat{\alpha }}_2$	3.3182	0.8182	3.0038	3.3641	0.8641	3.1734	3.2990	0.7990	2.3721
	${\widehat{\beta }}_{20}$	−0.4971	0.0029	0.0467	−0.4901	0.0099	0.0461	−0.4935	0.0065	0.0476
	${\widehat{\beta }}_{21}$	0.2626	0.0126	0.1341	0.2168	0.0332	0.1476	0.2236	0.0264	0.1482
	${\widehat{\sigma }}_2$	0.9412	0.0588	0.1204	0.9359	0.0641	0.1263	0.9097	0.0903	0.0834
	$\widehat{\rho }$	0.2458	0.0042	0.0204	0.4946	0.0054	0.0125	0.7464	0.0036	0.0047
	${\widehat{\alpha }}_1$	1.8309	0.3309	0.5297	1.7845	0.2845	0.4954	1.7743	0.2743	0.4677
	${\widehat{\beta }}_{10}$	0.2592	0.0092	0.0139	0.2642	0.0142	0.0140	0.2574	0.0074	0.0127
	${\widehat{\beta }}_{11}$	0.7136	0.0364	0.0397	0.7106	0.0394	0.0413	0.7209	0.0291	0.0383
	${\widehat{\sigma }}_1$	0.9412	0.0588	0.1147	0.9619	0.0381	0.2267	0.9613	0.0387	0.1071
100	${\widehat{\alpha }}_2$	2.8044	0.3044	0.8265	2.8406	0.3406	0.8715	2.8515	0.3515	0.8446
	${\widehat{\beta }}_{20}$	−0.4972	0.0028	0.0223	−0.4961	0.0039	0.0203	−0.4936	0.0064	0.0193
	${\widehat{\beta }}_{21}$	0.2391	0.0109	0.0657	0.2315	0.0185	0.0605	0.2243	0.0257	0.0612
	${\widehat{\sigma }}_2$	0.9753	0.0247	0.0405	0.9648	0.0352	0.0378	0.9597	0.0403	0.0408
	$\widehat{\rho }$	0.2462	0.0038	0.0091	0.4947	0.0053	0.0058	0.7465	0.0035	0.0021
	${\widehat{\alpha }}_1$	1.6703	0.1703	0.2198	1.6547	0.1547	0.1925	1.6470	0.1470	0.1822
	${\widehat{\beta }}_{10}$	0.2564	0.0064	0.0064	0.2635	0.0135	0.0067	0.2568	0.0068	0.0063
	${\widehat{\beta }}_{11}$	0.7166	0.0334	0.0202	0.7161	0.0339	0.0195	0.7236	0.0264	0.0191
	${\widehat{\sigma }}_1$	0.9568	0.0432	0.0454	0.9628	0.0372	0.0593	0.9616	0.0384	0.0397
200	${\widehat{\alpha }}_2$	2.6581	0.1581	0.3260	2.6704	0.1704	0.3513	2.6478	0.1478	0.2894
	${\widehat{\beta }}_{20}$	−0.5016	0.0016	0.0092	−0.4962	0.0038	0.0097	−0.4946	0.0054	0.0092
	${\widehat{\beta }}_{21}$	0.2409	0.0091	0.0297	0.2320	0.0180	0.0293	0.2253	0.0247	0.0279
	${\widehat{\sigma }}_2$	0.9806	0.0194	0.0169	0.9790	0.0210	0.0167	0.9779	0.0221	0.0145
	$\widehat{\rho }$	0.2469	0.0031	0.0044	0.4948	0.0052	0.0028	0.7476	0.0024	0.0010
	${\widehat{\alpha }}_1$	1.5992	0.0992	0.0727	1.5899	0.0899	0.0677	1.5875	0.0875	0.0602
	${\widehat{\beta }}_{10}$	0.2558	0.0058	0.0030	0.2581	0.0081	0.0027	0.2563	0.0063	0.0023
	${\widehat{\beta }}_{11}$	0.7215	0.0285	0.0090	0.7165	0.0335	0.0085	0.7262	0.0238	0.0074
	${\widehat{\sigma }}_1$	0.9578	0.0422	0.0152	0.9645	0.0355	0.0144	0.9627	0.0373	0.0126
500	${\widehat{\alpha }}_2$	2.5654	0.0654	0.1294	2.5518	0.0518	0.1170	2.5547	0.0547	0.0966
	${\widehat{\beta }}_{20}$	−0.4992	0.0008	0.0039	−0.4962	0.0038	0.0036	−0.4958	0.0042	0.0035
	${\widehat{\beta }}_{21}$	0.2433	0.0067	0.0117	0.2326	0.0174	0.0114	0.2282	0.0218	0.0109
	${\widehat{\sigma }}_2$	0.9935	0.0065	0.0072	0.9936	0.0064	0.0066	0.9888	0.0112	0.0053
	$\widehat{\rho }$	0.2523	0.0023	0.0018	0.5048	0.0048	0.0012	0.7489	0.0011	0.0004
	${\widehat{\alpha }}_1$	1.5900	0.0900	0.0392	1.5749	0.0749	0.0353	1.5674	0.0674	0.0318
	${\widehat{\beta }}_{10}$	0.2542	0.0042	0.0014	0.2579	0.0079	0.0014	0.2537	0.0037	0.0011
	${\widehat{\beta }}_{11}$	0.7265	0.0235	0.0050	0.7168	0.0332	0.0046	0.7285	0.0215	0.0038
	${\widehat{\sigma }}_1$	0.9604	0.0396	0.0083	0.9818	0.0182	0.0077	0.9663	0.0337	0.0070
1000	${\widehat{\alpha }}_2$	2.5390	0.0390	0.0527	2.5270	0.0270	0.0558	2.5330	0.0330	0.0435
	${\widehat{\beta }}_{20}$	−0.5005	0.0005	0.0019	−0.4996	0.0004	0.0019	−0.5034	0.0034	0.0018
	${\widehat{\beta }}_{21}$	0.2536	0.0036	0.0058	0.2387	0.0113	0.0057	0.2342	0.0158	0.0058
	${\widehat{\sigma }}_2$	0.9943	0.0057	0.0030	0.9956	0.0044	0.0033	0.9908	0.0092	0.0024
	$\widehat{\rho }$	0.2486	0.0014	0.0009	0.4984	0.0016	0.0006	0.7500	0.0000	0.0002

Table A2

Asymptotic behavior of the MLE for the BVSHN regression model with parameter vector $\theta ={({\alpha }_1,{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{\beta }}_2,{\sigma }_2,\rho )}^{\top }$. Simulation scenario 2.

		$\mathit{\rho }=-0.25$			$\mathit{\rho }=-0.75$
$\mathbf{n}$	Estimates	Mean	AVB	RMSE	Mean	AVB	RMSE
	${\widehat{\alpha }}_1$	3.6337	1.6337	7.8127	3.6073	1.6073	6.9422
	${\widehat{\beta }}_{10}$	0.2711	0.0211	0.0731	0.2740	0.0240	0.0706
	${\widehat{\beta }}_{11}$	0.6263	0.1237	0.2298	0.6258	0.1242	0.2177
	${\widehat{\sigma }}_1$	0.8203	0.1797	0.2296	0.8134	0.1866	0.2227
30	${\widehat{\alpha }}_2$	2.0779	1.5779	4.8352	2.0682	1.5682	4.8138
	${\widehat{\beta }}_{20}$	−0.5049	0.0049	0.0103	−0.5060	0.0060	0.0093
	${\widehat{\beta }}_{21}$	0.2635	0.0135	0.0331	0.2809	0.0309	0.0281
	${\widehat{\sigma }}_2$	0.4800	0.5200	0.5077	0.4941	0.5059	0.5365
	$\widehat{\rho }$	−0.2414	0.0086	0.0353	−0.7372	0.0128	0.0097
	${\widehat{\alpha }}_1$	2.9569	0.9569	2.7829	2.8585	0.8585	2.3274
	${\widehat{\beta }}_{10}$	0.2704	0.0204	0.0386	0.2727	0.0227	0.0381
	${\widehat{\beta }}_{11}$	0.6294	0.1206	0.1230	0.6352	0.1148	0.1233
	${\widehat{\sigma }}_1$	0.8561	0.1439	0.1528	0.8541	0.1459	0.1058
50	${\widehat{\alpha }}_2$	1.3617	0.8617	1.4747	1.2661	0.7661	1.1609
	${\widehat{\beta }}_{20}$	−0.5044	0.0044	0.0048	−0.5057	0.0057	0.0046
	${\widehat{\beta }}_{21}$	0.2622	0.0122	0.0149	0.2768	0.0268	0.0139
	${\widehat{\sigma }}_2$	0.6078	0.3922	0.4427	0.6216	0.3784	0.4356
	$\widehat{\rho }$	−0.2418	0.0082	0.0200	−0.7372	0.0128	0.0044
	${\widehat{\alpha }}_1$	2.4804	0.4804	0.8374	2.4197	0.4197	0.6970
	${\widehat{\beta }}_{10}$	0.2695	0.0195	0.0190	0.2678	0.0178	0.0182
	${\widehat{\beta }}_{11}$	0.6297	0.1203	0.0643	0.6356	0.1144	0.0656
	${\widehat{\sigma }}_1$	0.8863	0.1137	0.0502	0.8933	0.1067	0.0416
100	${\widehat{\alpha }}_2$	0.9343	0.4343	0.4176	0.9180	0.4180	0.3931
	${\widehat{\beta }}_{20}$	−0.5033	0.0033	0.0022	−0.5040	0.0040	0.0021
	${\widehat{\beta }}_{21}$	0.2621	0.0121	0.0068	0.2757	0.0257	0.0069
	${\widehat{\sigma }}_2$	0.7464	0.2536	0.3241	0.7576	0.2424	0.3472
	$\widehat{\rho }$	−0.2424	0.0076	0.0105	−0.7394	0.0106	0.0022
	${\widehat{\alpha }}_1$	2.2802	0.2802	0.3519	2.2894	0.2894	0.3037
	${\widehat{\beta }}_{10}$	0.2677	0.0177	0.0087	0.2665	0.0165	0.0086
	${\widehat{\beta }}_{11}$	0.6345	0.1155	0.0404	0.6381	0.1119	0.0377
	${\widehat{\sigma }}_1$	0.9140	0.0860	0.0289	0.9101	0.0899	0.0239
200	${\widehat{\alpha }}_2$	0.7295	0.2295	0.1567	0.7187	0.2187	0.1455
	${\widehat{\beta }}_{20}$	−0.5018	0.0018	0.0012	−0.5033	0.0033	0.0011
	${\widehat{\beta }}_{21}$	0.2606	0.0106	0.0035	0.2756	0.0256	0.0039
	${\widehat{\sigma }}_2$	0.8844	0.1156	0.3166	0.8616	0.1384	0.2346
	$\widehat{\rho }$	−0.2439	0.0061	0.0044	−0.7396	0.0104	0.0011
	${\widehat{\alpha }}_1$	2.2205	0.2205	0.1428	2.1741	0.1741	0.1050
	${\widehat{\beta }}_{10}$	0.2666	0.0166	0.0038	0.2653	0.0153	0.0034
	${\widehat{\beta }}_{11}$	0.6392	0.1108	0.0238	0.6527	0.0973	0.0216
	${\widehat{\sigma }}_1$	0.9155	0.0845	0.0142	0.9255	0.0745	0.0113
500	${\widehat{\alpha }}_2$	0.6064	0.1064	0.0619	0.5909	0.0909	0.0546
	${\widehat{\beta }}_{20}$	−0.5012	0.0012	0.0005	−0.5024	0.0024	0.0004
	${\widehat{\beta }}_{21}$	0.2589	0.0089	0.0015	0.2751	0.0251	0.0019
	${\widehat{\sigma }}_2$	1.0595	0.0595	0.2508	0.9863	0.0137	0.2383
	$\widehat{\rho }$	−0.2450	0.0050	0.0018	−0.7428	0.0072	0.0005
	${\widehat{\alpha }}_1$	2.1977	0.1977	0.0846	2.1432	0.1432	0.0556
	${\widehat{\beta }}_{10}$	0.2646	0.0146	0.0021	0.2583	0.0083	0.0018
	${\widehat{\beta }}_{11}$	0.6483	0.1017	0.0198	0.6528	0.0972	0.0175
	${\widehat{\sigma }}_1$	0.9163	0.0837	0.0104	0.9304	0.0696	0.0076
1000	${\widehat{\alpha }}_2$	0.5417	0.0417	0.0330	0.5428	0.0428	0.0288
	${\widehat{\beta }}_{20}$	−0.5001	0.0001	0.0002	−0.5011	0.0011	0.0002
	${\widehat{\beta }}_{21}$	0.2574	0.0074	0.0008	0.2698	0.0198	0.0013
	${\widehat{\sigma }}_2$	0.9870	0.0130	0.2335	1.0131	0.0131	0.1658
	$\widehat{\rho }$	−0.2453	0.0047	0.0010	−0.7450	0.0050	0.0004

Table A3

Asymptotic behavior of the MLE for the BVSHN regression model with parameter vector $\theta ={({\alpha }_1,{\mathbf{\beta }}_1,{\sigma }_1,{\alpha }_2,{\mathbf{\beta }}_2,{\sigma }_2,\rho )}^{\top }$. Simulation scenario 3.

		$\mathit{\rho }=-0.50$			$\mathit{\rho }=0.50$
$\mathbf{n}$	Estimates	Mean	AVB	RMSE	Mean	AVB	RMSE
	${\widehat{\alpha }}_1$	2.2025	1.4525	4.8426	3.1497	2.1497	8.8825
	${\widehat{\beta }}_{10}$	1.5016	0.0016	0.0051	1.1611	0.3389	0.1740
	${\widehat{\beta }}_{11}$	−0.5083	0.0083	0.0157	−0.3015	0.1985	0.2205
	${\widehat{\sigma }}_1$	0.3621	0.1379	0.2777	0.8701	0.8799	1.1743
30	${\widehat{\alpha }}_2$	2.4597	1.4597	5.1021	4.6733	1.6733	14.0574
	${\widehat{\beta }}_{20}$	−1.2397	0.0103	0.0747	−1.4625	0.2125	0.2703
	${\widehat{\beta }}_{21}$	0.7346	0.0154	0.2165	0.3738	0.1238	0.6797
	${\widehat{\sigma }}_2$	1.0721	0.4279	0.7349	1.3157	0.1843	0.3368
	$\widehat{\rho }$	−0.5034	0.0034	0.0222	0.4833	0.0167	0.0305
	${\widehat{\alpha }}_1$	1.4589	0.7089	1.3758	2.2209	1.2209	3.1101
	${\widehat{\beta }}_{10}$	1.5016	0.0016	0.0032	1.1633	0.3367	0.1539
	${\widehat{\beta }}_{11}$	−0.4959	0.0041	0.0087	−0.3019	0.1981	0.1472
	${\widehat{\sigma }}_1$	0.4527	0.0473	0.2213	1.0710	0.6790	1.0229
50	${\widehat{\alpha }}_2$	1.7864	0.7864	1.6633	3.7823	0.7823	3.5003
	${\widehat{\beta }}_{20}$	−1.2541	0.0041	0.0391	−1.4490	0.1990	0.1574
	${\widehat{\beta }}_{21}$	0.7399	0.0101	0.1247	0.3655	0.1155	0.3518
	${\widehat{\sigma }}_2$	1.2342	0.2658	0.6093	1.3969	0.1031	0.1848
	$\widehat{\rho }$	−0.5020	0.0020	0.0128	0.4940	0.0060	0.0197
	${\widehat{\alpha }}_1$	1.0700	0.3200	0.4133	1.6780	0.6780	1.1032
	${\widehat{\beta }}_{10}$	1.5013	0.0013	0.0030	1.1741	0.3259	0.1244
	${\widehat{\beta }}_{11}$	−0.5028	0.0028	0.0036	−0.3036	0.1964	0.0917
	${\widehat{\sigma }}_1$	0.5405	0.0405	0.2218	1.2143	0.5357	0.7571
100	${\widehat{\alpha }}_2$	1.3041	0.3041	0.4207	3.3002	0.3002	1.0830
	${\widehat{\beta }}_{20}$	−1.2460	0.0040	0.0209	−1.4391	0.1891	0.0881
	${\widehat{\beta }}_{21}$	0.7555	0.0055	0.0592	0.3616	0.1116	0.1575
	${\widehat{\sigma }}_2$	1.4473	0.0527	0.4827	1.4376	0.0624	0.0812
	$\widehat{\rho }$	−0.5005	0.0005	0.0057	0.4944	0.0056	0.0127
	${\widehat{\alpha }}_1$	0.8753	0.1253	0.1390	1.4408	0.4408	0.5224
	${\widehat{\beta }}_{10}$	1.4989	0.0011	0.0027	1.1834	0.3166	0.1095
	${\widehat{\beta }}_{11}$	−0.5022	0.0022	0.0020	−0.3199	0.1801	0.0580
	${\widehat{\sigma }}_1$	0.5360	0.0360	0.1221	1.2220	0.5280	0.5589
200	${\widehat{\alpha }}_2$	1.1227	0.1227	0.1624	2.9238	0.0762	0.3892
	${\widehat{\beta }}_{20}$	−1.2485	0.0015	0.0088	−1.4318	0.1818	0.0560
	${\widehat{\beta }}_{21}$	0.7536	0.0036	0.0269	0.3467	0.0967	0.0770
	${\widehat{\sigma }}_2$	1.5244	0.0244	0.3129	1.4715	0.0285	0.0321
	$\widehat{\rho }$	−0.5005	0.0005	0.0029	0.4947	0.0053	0.0089
	${\widehat{\alpha }}_1$	0.7872	0.0372	0.0587	1.3271	0.3271	0.2423
	${\widehat{\beta }}_{10}$	1.5003	0.0003	0.0016	1.1914	0.3086	0.0991
	${\widehat{\beta }}_{11}$	−0.5005	0.0005	0.0010	−0.3287	0.1713	0.0393
	${\widehat{\sigma }}_1$	0.5240	0.0240	0.0806	1.2461	0.5039	0.3871
500	${\widehat{\alpha }}_2$	1.0394	0.0394	0.0544	3.0449	0.0449	0.1462
	${\widehat{\beta }}_{20}$	−1.2495	0.0005	0.0048	−1.4272	0.1772	0.0426
	${\widehat{\beta }}_{21}$	0.7475	0.0025	0.0103	0.3431	0.0931	0.0346
	${\widehat{\sigma }}_2$	1.5178	0.0178	0.1112	1.4739	0.0261	0.0118
	$\widehat{\rho }$	−0.5004	0.0004	0.0012	0.4950	0.0050	0.0070
	${\widehat{\alpha }}_1$	0.7559	0.0059	0.0275	1.2969	0.2969	0.1549
	${\widehat{\beta }}_{10}$	1.5002	0.0002	0.0009	1.1919	0.3081	0.0967
	${\widehat{\beta }}_{11}$	−0.5003	0.0003	0.0007	−0.3292	0.1708	0.0340
	${\widehat{\sigma }}_1$	0.5173	0.0173	0.0226	1.2678	0.4822	0.3275
1000	${\widehat{\alpha }}_2$	1.0108	0.0108	0.0247	2.9811	0.0189	0.0670
	${\widehat{\beta }}_{20}$	−1.2501	0.0001	0.0017	−1.4259	0.1759	0.0363
	${\widehat{\beta }}_{21}$	0.7507	0.0007	0.0035	0.3365	0.0865	0.0220
	${\widehat{\sigma }}_2$	1.5163	0.0163	0.0487	1.4860	0.0140	0.0052
	$\widehat{\rho }$	−0.4999	0.0001	0.0005	0.4981	0.0019	0.0064

References

1. Birnbaum, Z.W.; Saunders, S.C. A new family of life distributions. J. Appl. Probab.; 1969; 6, pp. 319-327. [DOI: https://dx.doi.org/10.2307/3212003]

2. Birnbaum, Z.W.; Saunders, S.C. Estimation for a family of life distributions with applications to fatigue. J. Appl. Probab.; 1969; 6, pp. 328-347. [DOI: https://dx.doi.org/10.2307/3212004]

3. Díaz-García, J.A.; Leiva-Sánchez, V. A new family of life distributions based on the elliptically contoured distributions. J. Statist. Plann. Inference; 2005; 128, pp. 445-457. [DOI: https://dx.doi.org/10.1016/j.jspi.2003.11.007]

4. Vilca-Labra, F.; Leiva-Sánchez, V. A new fatigue life model based on the family of skew-elliptical distributions. Commun. Stat.-Theory Methods; 2006; 35, pp. 229-244. [DOI: https://dx.doi.org/10.1080/03610920500440065]

5. Gómez, H.W.; Olivares-Pacheco, J.F.; Bolfarine, H. An extension of the generalized Birnbaun-Saunders distribution. Stat. Probab. Lett.; 2009; 79, pp. 331-338. [DOI: https://dx.doi.org/10.1016/j.spl.2008.08.014]

6. Castillo, N.O.; Gómez, H.W.; Bolfarine, H. Epsilon Birnbaum–Saunders distribution family: Properties and inference. Stat. Pap.; 2011; 52, pp. 871-883. [DOI: https://dx.doi.org/10.1007/s00362-009-0293-x]

7. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. An alpha-power extension for the Birnbaum–Saunders distribution. Stat. J. Theor. Appl. Stat.; 2014; 48, pp. 896-912. [DOI: https://dx.doi.org/10.1080/02331888.2013.846910]

8. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. The Log-Linear Birnbaum-Saunders Power Model. Methodol. Comput. Appl. Probab.; 2017; 19, pp. 913-933. [DOI: https://dx.doi.org/10.1007/s11009-016-9526-3]

9. Moreno-Arenas, G.; Martínez-Flórez, G.; Bolfarine, H. Power Birnbaum-Saunders Student t distribution. Rev. Integr.; 2017; 35, pp. 51-70.

10. Balakrishnan, N.; Kundu, D. Birnbaum-Saunders distribution: A review of models, analysis, and applications. Mathematics; 2018; 35, pp. 1-41. [DOI: https://dx.doi.org/10.1002/asmb.2348]

11. Athayde, E.; Azevedo, A.; Barros, M.; Leiva, V. Failure rate of Birnbaum–Saunders distributions: Shape, change-point, estimation and robustness. Braz. J. Probab. Stat.; 2019; 33, pp. 301-328. [DOI: https://dx.doi.org/10.1214/17-BJPS389]

12. Arrué, J.; Arellano, R.; Gomez, H.W.; Leiva, V. On a new type of Birnbaum-Saunders models and its inference and application to fatigue data. J. Appl. Stat.; 2020; 47, pp. 2690-2710. [DOI: https://dx.doi.org/10.1080/02664763.2019.1668365]

13. Reyes, R.; Arruúe, J.; Víctor Leiva, V.; Martin-Barreiro, C. A new Birnbaum–Saunders distribution and its mathematical features applied to bimodal real-world data from environment and medicine. Mathematics; 2021; 9, 1891. [DOI: https://dx.doi.org/10.3390/math9161891]

14. Mazucheli, J.; Menezes, A.; Dey, S. The unit-Birnbaum-Saunders distribution with applications. Chil. J. Stat.; 2018; 9, pp. 47-57.

15. Rieck, J.R.; Nedelman, J.R. A log-linear model for the Birnbaum-Saunders distribution. Technometrics; 1991; 33, pp. 51-60.

16. Barros, M.; Paula, G.A.; Leiva, V. A new class of survival regression models with heavy-tailed errors: Robustness and diagnostics. Lifetime Data Anal.; 2008; 14, pp. 316-332. [DOI: https://dx.doi.org/10.1007/s10985-008-9085-1]

17. Leiva, V.; Vilca, F.; Balakrishnan, N.; Sanhueza, A. A skewed sinh-normal distribution and its properties and application to air pollution. Commun. Stat. Theory Methods; 2010; 39, pp. 426-443. [DOI: https://dx.doi.org/10.1080/03610920903140171]

18. Lemonte, A. A log-Birnbaum-Saunders Regression Model with Asymmetric Errors. J. Stat. Comput. Simul.; 2012; 82, pp. 1775-1787. [DOI: https://dx.doi.org/10.1080/00949655.2011.595715]

19. Santana, L.; Vilca, F.; Leiva, V. Influence analysis in skew-Birnbaum-Saunders regression models and applications. J. Appl. Stat; 2011; 38, pp. 1633-1649. [DOI: https://dx.doi.org/10.1080/02664763.2010.515679]

20. Kundu, D.; Balakrishnan, N.; Jamalizadeh, A. Bivariate Birnbaum-Saunders distribution and associated inference. J. Multivar. Anal.; 2010; 101, pp. 113-125. [DOI: https://dx.doi.org/10.1016/j.jmva.2009.05.005]

21. Lemonte, A.J.; Martínez-Flórez, G.; Moreno-Arenas, G. Multivariate Birnbaum–Saunders distribution: Properties and associated inference. J. Stat. Comput. Simul.; 2015; 85, pp. 374-392. [DOI: https://dx.doi.org/10.1080/00949655.2013.823964]

22. Lemonte, A.J. Multivariate Birnbaum–Saunders regression model. J. Stat. Comput. Simul.; 2013; 83, pp. 2244-2257. [DOI: https://dx.doi.org/10.1080/00949655.2012.688054]

23. Kundu, D. Bivariate sinh-normal distribution and a related model. Braz. J. Probab. Stat.; 2015; 29, pp. 590-607. [DOI: https://dx.doi.org/10.1214/13-BJPS235]

24. Martínez-Flórez, G.; Azevedo-Farias, R.B.; Moreno-Arenas, G. Multivariate log-Birnbaum–Saunders regression models. Commun. Stat.-Theory Methods; 2017; 46, pp. 10166-10178. [DOI: https://dx.doi.org/10.1080/03610926.2016.1231818]

25. Martínez-Flórez, G.; Azevedo-Farias, R.B.; Tovar-Falón, R. An Exponentiated Multivariate Extension for the Birnbaum-Saunders Log-Linear Model. Mathematics; 2022; 10, 1299. [DOI: https://dx.doi.org/10.3390/math10081299]

26. Paolino, P. Maximum likelihood estimation of models with beta-distributed dependent variables. Political Anal.; 2001; 9, pp. 325-346. [DOI: https://dx.doi.org/10.1093/oxfordjournals.pan.a004873]

27. Cribari-Neto, F.; Vasconcellos, K.L.P. Nearly unbiased maximum likelihood estimation for the beta distribution. J. Stat. Comput. Simul.; 2002; 72, pp. 107-118. [DOI: https://dx.doi.org/10.1080/00949650212144]

28. Kieschnick, R.; McCullough, B.D. Regression analysis of variates observed on (0, 1): Percentages, proportions and fractions. Stat. Model.; 2003; 3, pp. 193-213. [DOI: https://dx.doi.org/10.1191/1471082X03st053oa]

29. Ferrari, S.; Cribari-Neto, F. Beta regression for modelling rates and proportions. J. Appl. Stat.; 2004; 31, pp. 799-815. [DOI: https://dx.doi.org/10.1080/0266476042000214501]

30. Vasconcellos, K.L.P.; Cribari-Neto, F. Improved maximum likelihood estimation in a new class of beta regression models. Braz. J. Probab. Stat.; 2005; 19, pp. 13-31.

31. Branscum, A.J.; Johnson, W.O.; Thurmond, M.C. Bayesian beta regression: Applications to household expenditure data and genetic distance between foot-and-mouth diseases viruses. Aust. N. Z. J. Stat.; 2007; 49, pp. 287-301. [DOI: https://dx.doi.org/10.1111/j.1467-842X.2007.00481.x]

32. Bayes, C.L.; Bazan, J.L.; García, C. A new robust regression model for proportions. Bayesian Anal.; 2012; 7, pp. 841-866. [DOI: https://dx.doi.org/10.1214/12-BA728]

33. Ospina, R.; Ferrari, S.L.P. Inflated beta distribution. Stat. Pap.; 2010; 51, pp. 111-126. [DOI: https://dx.doi.org/10.1007/s00362-008-0125-4]

34. Ospina, R.; Ferrari, S.L.P. A general class of zero-or-one inflated beta regression models. Comput. Stat. Data Anal.; 2012; 56, pp. 1609-1623. [DOI: https://dx.doi.org/10.1016/j.csda.2011.10.005]

35. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. Doubly censored power-normal regression models with inflation. Test; 2014; 24, pp. 265-286. [DOI: https://dx.doi.org/10.1007/s11749-014-0406-2]

36. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. Power-models for proportions with zero/one excess. Appl. Math. Inf. Sci.; 2018; 24, pp. 293-303. [DOI: https://dx.doi.org/10.18576/amis/120203]

37. Martínez-Flórez, G.; Tovar-Falón, R.; Barrera-Causil, C. Inflated Unit-Birnbaum-Saunders Distribution. Mathematics; 2022; 10, 667. [DOI: https://dx.doi.org/10.3390/math10040667]

38. Marshall, A.W. Some comments on the hazard gradient. Stoch. Process. Their Appl.; 1975; 3, pp. 293-300. [DOI: https://dx.doi.org/10.1016/0304-4149(75)90028-9]

39. Joe, H. Asymptotic efficiency of the two-stage estimation method for copula-based models. J. Multivar. Anal.; 2005; 94, pp. 401-419. [DOI: https://dx.doi.org/10.1016/j.jmva.2004.06.003]

40. R Development Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021; Available online: http://www.R-project.org (accessed on 28 February 2022).

41. Meintanis, S.G. Test of fit for Marshall-Olkin distributions with applications. J. Stat. Plan. Inference; 2007; 137, pp. 3954-3963. [DOI: https://dx.doi.org/10.1016/j.jspi.2007.04.013]