## A new unit-bimodal distribution based on correlated Birnbaum-Saunders random variables

Roberto Vila ∗ 1 , Helton Saulo † 1 , Felipe Quintino ‡ 1 , and Peter Z¨rnig o § 1

1 Department of Statistics, University of Bras´ ılia, 70910-900, Bras´ ılia, Brazil

November 5, 2024

## Abstract

In this paper, we propose a new distribution over the unit interval which can be characterized as a ratio of the type Z d = Y/ X ( + Y ) where X and Y are two correlated BirnbaumSaunders random variables. The density of Z may be unimodal or bimodal. Simple expressions for the cumulative distribution function, moment-generating function and moments are obtained. Moreover, the stress-strength probability between X and Y is calculated explicitly in the symmetric case, that is, when the respective scale parameters are equal. Two applications of the ratio distribution are discussed.

Keywords. Unit distribution, Type I ratio, Maximum likelihood method, Monte Carlo simulation, R+ software.

## 1 Introduction

Traditionally, the activities in statistical distribution theory have been concentrated on distributions with unbounded support. Only in recent years, attention has been devoted to the construction of distributions with bounded support. Bounded data have applications in economics, finance, biology and medicine among others, see e.g., the introductory remarks in [25] and [26]. Such data arise naturally in the form of rates and proportions. Some of the most prominent bounded distributions are the beta distribution and their generalizations and the distributions of Kumaraswamy and Topp-Leone.

A natural way to create a model for bounded data is to transform one or two unbounded random variables in a bounded variable Z . For example, Z = X/ X ( + 1), Z = exp( -X ), Z = tanh( X ), Z = min { X,Y } / max { X,Y } , Z = ( X/Y ) 1 { X&lt;Y } and Z = X/ X ( + Y ) have the support (0 1), ,

∗ rovig161@gmail.com

‡ felipesquintino2@gmail.com

† heltonsaulo@gmail.com

§ peter@unb.br

when X and Y are positive variables. In applications of the last model it is usually assumed that X and Y are independent. We consider some specific approaches using these transformations or more complicated ones. These techniques are simple transformations of one or two variables: [25] derive a bimodal beta distribution, applying a quadratic transformation technique due to Elal-Olivero to the beta distribution. [2] apply an approach of Marshall and Olkin to the reduced Kies distribution. [9] propose a general framework in which a monotone function is applied to a single variable. [28] applies a procedure due to [10] to a uniformly distributed variable, resulting in a three-parameter family of distributions over the unit interval. The authors in [27] (unit log) use the approach Z = X/ X ( + Y ) where X and Y are the components of a two-dimensional random vector, following a bivariate log-symmetric distribution. In this model, X and Y may be correlated.

The distribution of ratios of the latter type is of great importance in practice, for example in modeling COVID-19-related death rates, see [6].

On the other hand, the Birnbaum-Saunders (BS) fatigue-life model, introduced by [4], is based on renewal theory principles. This model focuses on the number of cycles needed to exceed a critical threshold and cause fatigue crack propagation. In the same year, [5] provided maximum likelihood estimates for the BS model parameters. Later, Desmond (1985) [11], in his work on stochastic models of failure in random environments, offered an alternative derivation of the distribution using a biological model and relaxed several of the assumptions made by [4]. The BS distribution has been extensively studied and applied across various fields, including business, industry, economics, finance, management, engineering, and medical sciences. For a review of recent work on the BS distribution, see [24] and the references therein. Bivariate (and multivariate) versions of the BS distribution have also been explored. [16] introduced the bivariate BS distribution with five parameters using the bivariate normal distribution function. [15] revisited the bivariate BS distribution and discussed several new properties. [24] introduced a mean-based bivariate BS distribution, demonstrating that the mean is one of its parameters.

In the present paper we propose the distribution of a random variable of the type Z = Y/ X ( + Y ), where X and Y are two (possibly correlated) Birnbaum-Saunders random variables. The article is structured as follows. In Section 2, we introduce the unit-bimodal Birnbaum-Saunders model, which we name by type 1 and denote by UBBS1 to differentiate from the unit-bimodal BirnbaumSaunders (UBBS) distribution studied in [18], and some of its forms are outlined in graphs. Note that bimodality in the model is observed as a function of its shape and correlation parameters. For a unit Birnbaum-Saunders model based on the transformation Z = exp( -X ), with X having a BS distribution, see [19, 20]. Since the variable Z in [19, 20] exhibits decreasing, upside down bathtub and bathtub shaped density, this model is not suitable to capture bimodal pattern in data. In Section 3, we show one of the main properties of the UBBS1 model, that is, the corresponding random variable originally comes from a type I ratio [3, 13] of two variables following a bivariate BS distribution.

Also in this section, we derive the cumulative distribution function (CDF), which is used to obtain simple expressions of the moment-generating function and moments, and to explicitly calculate the stress-strength probability in the case of the presence of correlation and equality of scale parameters. Finally, the maximum likelihood and maximum product of spacings methods for parameter estimation are presented. In Section 4, we carry out a Monte Carlo simulation study to evaluate the performance of the maximum likelihood and maximum product of spacings estimators by means of their bias and root mean square error. In this section, we also compare the performance

of the estimation methods when bimodality is present. In Section 5, we present two applications to income-consumption and body mass data. The first data set presents unimodal characteristics, and in this case, the new UBBS1 model proved to be more appropriate than the Beta and UBBS models (from [18]). The second data set has bimodal characteristics and was previously modeled by [18] using the UBBS model. Now, we show that the UBBS1 model also provides a good fit for this data. Finally, in Section 6, some concluding remarks are made.

## 2 The unit-bimodal Birnbaum-Saunders distribution of type 1

We say that a continuous random variable Z , with support (0 1), follows a unit-bimodal BS dis-, tribution of type 1 (UBBS1) with parameter vector θ ρ = ( α , α , β 1 2 1 , β 2 , ρ ) ⊤ , α 1 &gt; 0 , α 2 &gt; 0 , β 1 &gt; 0 , β 2 &gt; 0 and ρ ∈ -( 1 1), denoted by , Z ∼ UBBS1( θ ρ ), if its probability density function (PDF) is given by (for 0 &lt; z &lt; 1 and s = 1 /z -1)

$$\Lambda _ { T } ( z ; \theta _ { \rho } ) \ & = \frac { \exp \left \{ \frac { 1 } { 1 - \rho ^ { 2 } } \left ( \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } \right ) \right \} } { 4 \pi \alpha _ { 1 } \alpha _ { 2 } \sqrt { 1 - \rho ^ { 2 } } } \, \frac { ( s + 1 ) ^ { 2 } } { s } \, \exp \left \{ - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) } \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) \right \} \\ & \quad \times \left \{ \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) K _ { 0 } \left ( \frac { \sqrt { u _ { \rho } v _ { \rho } } } { 1 - \rho ^ { 2 } } \right ) + \left [ \sqrt { \frac { ( v _ { \rho } / u _ { \rho } ) s } { \beta _ { 1 } \beta _ { 2 } } } + \sqrt { \frac { 1 } { ( v _ { \rho } / u _ { \rho } ) s } } \right ] K _ { 1 } \left ( \frac { \sqrt { u _ { \rho } v _ { \rho } } } { 1 - \rho ^ { 2 } } \right ) \right \},$$

where

$$u _ { \rho } \equiv \frac { s } { \alpha _ { 1 } ^ { 2 } \beta _ { 1 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } \beta _ { 2 } } - \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { s } { \beta _ { 1 } \beta _ { 2 } } }, \quad v _ { \rho } \equiv \frac { \beta _ { 1 } } { \alpha _ { 1 } ^ { 2 } s } + \frac { \beta _ { 2 } } { \alpha _ { 2 } ^ { 2 } } - \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s } }.$$

Furthermore, in formula (1),

$$K _ { \lambda } ( x ) & = \frac { 1 } { 2 } \int _ { 0 } ^ { \infty } w ^ { \lambda - 1 } \exp \left \{ - \frac { 1 } { 2 } \, x \left ( w + \frac { 1 } { w } \right ) \right \} d w, \quad x > 0, \lambda \in \mathbb { R }, \\ & = \frac { \pi } { 2 \sin ( \lambda \pi ) } \left [ I _ { - \lambda } ( x ) - I _ { \lambda } ( x ) \right ]$$

is the modified Bessel function of the third kind (also known as modified Bessel function of the second kind) [1, 14], with

$$I _ { \lambda } ( x ) = \left ( \frac { x } { 2 } \right ) ^ { \lambda } \sum _ { k = 0 } ^ { \infty } \frac { 1 } { k! \Gamma ( \lambda + k + 1 ) } \left ( \frac { x } { 2 } \right ) ^ { 2 k }$$

being the modified Bessel function of the first kind. When λ ∈ Z the following integral representation of I λ is valid [1, p. 376]: I λ ( x ) = (1 /π ) ∫ π 0 exp { x cos( θ ) } cos( λθ )d θ.

Remark 2.1. A simple algebraic manipulation shows that, when β 1 = β 2 , f Z (1 / 2 -w ; θ ρ ) = f Z (1 / 2 + w ; θ ρ ), for all 0 &lt; w &lt; 1 / 2. That is, The UBBS1 PDF in (1) is symmetric around z 0 = 1 2 provided / β 1 = β 2 (see Figure 4). Furthermore, in this case, the median and the mean of a UBBS1 distribution both occur at z 0 = 1 2. /

When ρ = 0, that is, X and Y are independent, the PDF (1) of Z is simplified as (for 0 &lt; z &lt; 1 and s = 1 /z -1)

$$f _ { Z } ( z ; \theta _ { 0 } ) \ = \ & \frac { \exp \left ( \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } \right ) } { 4 \pi \alpha _ { 1 } \alpha _ { 2 } } \ \frac { ( s + 1 ) ^ { 2 } } { s } \\ & \times \left \{ \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) K _ { 0 } \left ( \sqrt { u v } \right ) + \left [ \sqrt { \frac { ( v / u ) s } { \beta _ { 1 } \beta _ { 2 } } } + \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { ( v / u ) s } } \right ] K _ { 1 } \left ( \sqrt { u v } \right ) \right \},$$

where u = u 0 and v = v 0 , and θ 0 = θ ρ =0 = ( α , α , β 1 2 1 , β 2 , 0) ⊤ .

Remark 2.2. Let S = X/Y be the type II ratio, where X ∼ BS( α , β 1 1 ) and Y ∼ BS( α , β 2 2 ) are independent random variables (that is, ρ = 0), and let f S ( s ; θ 0 ) denote its respective PDF. By using the identity f S ( s ; θ 0 ) = ∫ ∞ 0 yf X ( sy f ) Y ( y )d y it is simple to check that the PDF of S is (for s &gt; 0)

$$\dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\ f _ { S } ( s ; \theta _ { 0 } ) \ & = \frac { \exp \left ( \frac { 1 } { \alpha _ { 1 } ^ { 2 } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } } \right ) } { 4 \pi \alpha _ { 1 } \alpha _ { 2 } } \frac { 1 } { s } \\ & \quad \times \left \{ \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) K _ { 0 } \left ( \sqrt { u v } \right ) + \left [ \sqrt { \frac { ( v / u ) s } { \beta _ { 1 } \beta _ { 2 } } } + \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { ( v / u ) s } } \right ] K _ { 1 } \left ( \sqrt { u v } \right ) \right \} \\ & = \frac { 1 } { ( s + 1 ) ^ { 2 } } \ f _ { Z } \left ( \frac { 1 } { s + 1 } ; \theta _ { 0 } \right ), \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdclots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots.$$

where f Z ( z ; θ 0 ) is as given in (4). The PDF (5) of S with β 1 = β 2 = 1 has appeared in reference [21].

Figures 1, 2, 3 and 4 show the behavior of f Z for some parameter choices. Note that the density modality changes according to the values of α 1 and α 2 parameters. Figure 4 shows the symmetry of the distribution in the case where the scale parameters are equal. Figure 2 indicates that, for θ ρ = ( α , 1 0 7 . , 1 1 . , 0 9 . , 0 6) . ⊤ , the density exhibits a bimodal behavior when α 1 ⩾ 1 5. . A similar pattern is observed when α 1 is fixed and α 2 ⩾ 2.

Remark 2.3. Note that the equation f ′ Z ( z ; θ ρ ) = 0, where f ′ Z ( z ; θ ρ ) is the derivative of f Z ( z ; θ ρ ) in (1) with respect to z , can equivalently be written as G z ( ) = H z ( ), 0 &lt; z &lt; 1, where G and H are non-polynomial functions (in terms of modified Bessel functions). In other words, the critical points of f Z are determined by the intersection points of (the graphs of the functions) G and H . The analysis of the intersection points of G and H can not be performed analytically, but graphical and numerical analyses are required. From Figures 1, 2 and 3 we come to the conclusion that the number of intersection points of G and H is bounded above by 3. That is, the function f Z has at most 3 critical points. Using this fact plus the fact that lim z → 0 + f Z ( z ; θ ρ ) = lim z → 1 -f Z ( z ; θ ρ ) = 0 we conclude that f Z is unimodal (in the case that it has a single critical point) or bimodal (in the case of exactly 3 critical points).

Figure 1: Plot of the PDF f Z with varying parameters α 1 (left) and α 2 (right).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-22c2ff00-0a75-464a-b81c-17692465e76e.jpg)

Figure 2: Plot of the PDF f Z with varying parameters β 1 (left) and β 2 (right).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-76f940e3-54e3-458b-ad83-cc55259529cb.jpg)

## a 1 = 1 ,   a 2 = 0.7 , b 1 = 1.1 ,   b 2 = 0.9

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-d0aa4db6-f9a5-4f2f-9f3c-2f8209806e68.jpg)

z

Figure 3: Plot of the PDF f Z with varying parameter ρ .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-18bbeb88-1145-431b-b292-f7046fb62e8e.jpg)

z

Figure 4: Plot of the PDF f Z with varying parameters β 1 = β 2 .

## 3 Characterization and estimation

In this section, we establish some mathematical results of the UBBS1 distribution. In particular, we define the UBBS1 model as a ratio, present the cumulative distribution function, the moment-

generating function and moments, and calculate the stress-strength probability. We also derive parameter estimators for the UBBS1 model based on maximum likelihood and maximum product of spacings methods.

## 3.1 UBBS1 model arising as a ratio

In this section, we establish one of the most important properties of the UBBS1 model. This property characterizes a UBBS1 random variable Z in (1) as a ratio of the form Y/ X ( + Y ), where the random vector ( X,Y ) ⊤ follows a bivariate Birnbaum-Saunders distribution as discussed in [15]. Indeed, let X,Y and Z be continuous and positive random variables such that

$$Z \stackrel { d } { = } \frac { Y } { X + Y },$$

being = equality in distribution, and the ratio in (6) has support in the unit interval (0 d , 1). The random variable Z = Y/ X ( + Y ) is known in the literature as type I ratio [3, 13].

From (6), it is simple to verify that the cumulative distribution function (CDF) of Z is given by

$$F _ { Z } ( z ) = \mathbb { P } \left ( \frac { X } { Y } \geq s \right ) = 1 - \int _ { 0 } ^ { \infty } \mathbb { P } \left ( X \leqslant s y | Y = y \right ) f _ { Y } ( y ) \mathrm d y, \quad s = \frac { 1 } { z } - 1.$$

Hence, the PDF of Z can be expressed as

$$f _ { Z } ( z ) = ( s + 1 ) ^ { 2 } \int _ { 0 } ^ { \infty } y f _ { X, Y } ( s y, y ) \text{d} y.$$

Remark 3.1. Equation in (8) can also be obtained using the Jacobian method by determining the marginal with respect to Z of the random vector ( Z, U ) ⊤ , where Z = Y/ X ( + Y ) and U = Y .

Now, let's suppose that ( X,Y ) ⊤ has a bivariate Birnbaum-Saunders distribution with parameter vector θ ρ = ( α , α , β 1 2 1 , β 2 , ρ ) ⊤ , denoted by ( X,Y ) ⊤ ∼ BS( θ ), that is, its joint density is given by

$$f _ { X, Y } ( x, y ) = \phi _ { 2 } ( a ( x ; \alpha _ { 1 }, \beta _ { 1 } ), a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) ; \rho ) a ^ { \prime } ( x ; \alpha _ { 1 }, \beta _ { 1 } ) a ^ { \prime } ( y ; \alpha _ { 2 }, \beta _ { 2 } ),$$

where ϕ 2 ( u, v ; ρ ) denotes the standard bivariate normal density function;

$$\phi _ { 2 } ( u, v ; \rho ) = \frac { 1 } { 2 \pi \sqrt { 1 - \rho ^ { 2 } } } \, \exp \left \{ - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \, ( u ^ { 2 } + v ^ { 2 } - 2 \rho u v ) \right \}, \quad ( u, v ) ^ { \top } \in \mathbb { R } ^ { 2 }$$

and

$$a ( t ; \alpha, \beta ) = \frac { 1 } { \alpha } \left ( \sqrt { \frac { t } { \beta } } - \sqrt { \frac { \beta } { t } } \right ), \ \ a ^ { \prime } ( t ; \alpha, \beta ) = \frac { 1 } { 2 \alpha t } \left ( \sqrt { \frac { t } { \beta } } + \sqrt { \frac { \beta } { t } } \right ), \ \ t > 0.$$

Hence, the density (8) of Z is written as

$$f _ { Z } ( z ) = ( s + 1 ) ^ { 2 } \int _ { 0 } ^ { \infty } y \phi _ { 2 } ( a ( s y ; \alpha _ { 1 }, \beta _ { 1 } ), a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) ; \rho ) a ^ { \prime } ( s y, \alpha _ { 1 }, \beta _ { 1 } ) a ^ { \prime } ( y, \alpha _ { 2 }, \beta _ { 2 } ) d y. \quad ( 1 0 )$$

Since

$$a ^ { 2 } ( s y ; \alpha _ { 1 }, \beta _ { 1 } ) + a ^ { 2 } ( y ; \alpha _ { 2 }, \beta _ { 2 } ) - 2 \rho a ( s y ; \alpha _ { 1 }, \beta _ { 1 } ) a ( y ; \alpha _ { 2 }, \beta _ { 2 } )$$

$$\begin{matrix} \end{matrix}$$

$$= \left ( u _ { \rho } y + \frac { v _ { \rho } } { y } \right ) - 2 \left [ \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) \right ],$$

with u ρ and v ρ as given in (2), note that the function ϕ 2 ( a sy ( ; α , β 1 1 ) , a ( y α , β ; 2 2 ); ρ ) in (10) can be expressed as

$$\phi _ { 2 } ( a ( s y ; \alpha _ { 1 }, \beta _ { 1 } ), & \, a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) ; \rho ) = \frac { 1 } { 2 } \, \exp \left \{ - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \left ( u _ { \rho } y + \frac { v _ { \rho } } { y } \right ) \right \} \\ & \times \frac { 1 } { \pi \sqrt { 1 - \rho ^ { 2 } } } \, \exp \left \{ \frac { 1 } { 1 - \rho ^ { 2 } } \left [ \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) \right ] \right \}.$$

Furthermore, notice that

$$a ^ { \prime } ( s y ; \alpha _ { 1 }, \beta _ { 1 } ) a ^ { \prime } ( y ; \alpha _ { 2 }, \beta _ { 2 } ) = \frac { 1 } { 4 \alpha _ { 1 } \alpha _ { 2 } s y } \left [ \frac { 1 } { y ^ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s } } + \frac { 1 } { y } \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) + \sqrt { \frac { s } { \beta _ { 1 } \beta _ { 2 } } } \right ].$$

By replacing (11) and (12) in formula (10), we have

$$\psi _ { \ } & \psi _ { \ } & \psi _ { \ } & \psi _ { \ } & \psi _ { \ } & \psi _ { \ } \\ f _ { Z } ( z ) & = \frac { 1 } { \pi \sqrt { 1 - \rho ^ { 2 } } } \, \frac { ( s + 1 ) ^ { 2 } } { 4 \alpha _ { 1 } \alpha _ { 2 } s } \, \exp \left \{ \frac { 1 } { 1 - \rho ^ { 2 } } \left [ \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) \right ] \right \} \\ & \quad \times \left [ \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s } } \, \mathbb { I } _ { 1 } + \left ( \sqrt { \frac { \beta _ { 2 } s } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s } } \right ) \mathbb { I } _ { 2 } + \sqrt { \frac { s } { \beta _ { 1 } \beta _ { 2 } } } \, \mathbb { I } _ { 3 } \right ],$$

where and

$$\mathbb { I } _ { 3 } = \frac { 1 } { 2 } \, \int _ { 0 } ^ { \infty } \exp \left \{ - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \left ( u _ { \rho } y + \frac { v _ { \rho } } { y } \right ) \right \} \mathrm d y = \sqrt { \frac { v _ { \rho } } { u _ { \rho } } } \, K _ { 1 } \left ( \frac { \sqrt { u _ { \rho } v _ { \rho } } } { 1 - \rho ^ { 2 } } \right ).$$

In the above, K λ ( x ), x &gt; 0, denotes the modified Bessel function of the third kind defined in (3).

Finally, by using that the Bessel function is symmetric with respect to the index λ : K λ ( x ) = K -λ ( x ), and by employing the above formulas of I 1 , I 2 and I 3 , from (13) the PDF of Z can be expressed as in (1).

$$\mathbb { I } _ { 1 } & = \frac { 1 } { 2 } \int _ { 0 } ^ { \infty } \frac { 1 } { y ^ { 2 } } \, \exp \left \{ - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \left ( u _ { \rho } y + \frac { v _ { \rho } } { y } \right ) \right \} \mathrm d y = \sqrt { \frac { u _ { \rho } } { v _ { \rho } } } \, K _ { - 1 } \left ( \frac { \sqrt { u _ { \rho } v _ { \rho } } } { 1 - \rho ^ { 2 } } \right ), \\ \mathbb { I } _ { 2 } & = \frac { 1 } { 2 } \int _ { 0 } ^ { \infty } \frac { 1 } { y } \, \exp \left \{ - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \left ( u _ { \rho } y + \frac { v _ { \rho } } { y } \right ) \right \} \mathrm d y = K _ { 0 } \left ( \frac { \sqrt { u _ { \rho } v _ { \rho } } } { 1 - \rho ^ { 2 } } \right ) \$$

## 3.2 The cumulative distribution function

For ( X,Y ) ⊤ ∼ BS( θ ρ ), it holds that X ∼ BS( α , β 1 1 ) and Y ∼ BS( α , β 2 2 ). Furthermore, the conditional CDF of X , given Y = y , is given by [16, Theorem 3.1]

$$\mathbb { P } \left ( X \leqslant x | Y = y \right ) = \Phi \left \{ \frac { a ( x ; \alpha _ { 1 }, \beta _ { 1 } ) - \rho a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) } { \sqrt { 1 - \rho ^ { 2 } } } \right \},$$

where a t ( ; α, β ) is as given in (9) and Φ is the CDF of the standard normal distribution. Hence, from (7) it follows that

$$& \L ^ { \prime } _ { Z } ( z ; \boldsymbol \theta _ { \rho } ) = 1 - \int _ { 0 } ^ { \infty } \mathbb { P } \left ( X \leqslant s y | Y = y \right ) f _ { Y } ( y ) \mathrm d y, \quad s = \frac { 1 } { z } - 1, \\ & \quad = 1 - \int _ { 0 } ^ { \infty } \Phi \left \{ \frac { a ( s y ; \alpha _ { 1 }, \beta _ { 1 } ) - \rho a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) } { \sqrt { 1 - \rho ^ { 2 } } } \right \} \phi ( a ( y ; \alpha _ { 2 }, \beta _ { 2 } ) ) a ^ { \prime } ( y ; \alpha _ { 2 }, \beta _ { 2 } ) \mathrm d y,$$

where ϕ denotes the PDF of the standard normal distribution. By using the change of variable w = a y ( ; α , β 2 2 ), we write the above expression as

$$1$$

$$= \frac { 1 } { 2 } - \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { a ( s a ^ { - 1 } ( w ; \alpha _ { 2 }, \beta _ { 2 } ) ; \alpha _ { 1 }, \beta _ { 1 } ) - \rho w } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \right \} \phi ( w ) \text{d} w,$$

$$& 1 - \int _ { - \infty } ^ { \infty } \Phi \left \{ \frac { a ( s a ^ { - 1 } ( w ; \alpha _ { 2 }, \beta _ { 2 } ) ; \alpha _ { 1 }, \beta _ { 1 } ) - \rho w } { \sqrt { 1 - \rho ^ { 2 } } } \right \} \phi ( w ) \mathrm d w \\ & = \frac { 1 } { 2 } - \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { a ( s a ^ { - 1 } ( w ; \alpha _ { 2 }, \beta _ { 2 } ) ; \alpha _ { 1 }, \beta _ { 1 } ) - \rho w } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \right \} \phi ( w ) \mathrm d w,$$

where we have used the well-known identity Φ( x ) = (1 / 2)[1 + erf( x/ √ 2)]. In the above, a -1 ( ; t α, β ) denotes the inverse function of a t ( ; α, β ). As a -1 ( ; t α, β ) = ( β/ 4)[ αt + √ ( αt ) 2 +4] 2 , the integral above can be written as

$$\frac { 1 } { 2 } - \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { \left ( \frac { \alpha _ { 2 } } { \alpha _ { 1 } } \sqrt { \frac { s \beta _ { 2 } } { \beta _ { 1 } } } - \rho \right ) } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \, w + \frac { \sqrt { 2 } ( s \beta _ { 2 } - \beta _ { 1 } ) } { \alpha _ { 1 } \sqrt { s \beta _ { 1 } \beta _ { 2 } ( 1 - \rho ^ { 2 } ) } } \frac { 1 } { \alpha _ { 2 } w + \sqrt { ( \alpha _ { 2 } w ) ^ { 2 } + 4 } } \right \} \phi ( w ) \text{d} w.$$

Therefore, the CDF of the UBBS1 random variable Z is given by (for s = 1 /z -1 and 0 &lt; z &lt; 1)

$$F _ { Z } ( z ; \theta _ { \rho } ) = \frac { 1 } { 2 } - \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { \left ( \frac { \alpha _ { 2 } } { \alpha _ { 1 } } \sqrt { \frac { s \beta _ { 2 } } { \beta _ { 1 } } } - \rho } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \, w + \frac { \sqrt { 2 } ( s \beta _ { 2 } - \beta _ { 1 } ) } { \alpha _ { 1 } \sqrt { s \beta _ { 1 } \beta _ { 2 } ( 1 - \rho ^ { 2 } ) } } \, \frac { 1 } { \alpha _ { 2 } w + \sqrt { ( \alpha _ { 2 } w ) ^ { 2 } + 4 } } \right \} \phi ( w ) \mathrm d w.$$

## 3.3 Moment-generating function and moments

Let Z ∼ UBBS1( θ ρ ). Since Z has unitary support it is clear that the moment-generating function (MGF), denoted by M Z ( ) = t E [exp( tZ )], t ∈ R , and the (raw) moments, denoted by µ n = E ( Z n ),

n ∈ N , are finite. In particular, Z has a finite MGF M Z ( ) in some open interval containing 0, then t it holds M ( n ) Z ( ) t ∣ ∣ t =0 = µ n , and M Z ( ) may be expanded in a power series about 0 as follows t

$$M _ { Z } ( t ) = 1 + \sum _ { n = 1 } ^ { \infty } \frac { \mu _ { n } } { n! } \, t ^ { n }.$$

On the other hand, by using the well-known formula E ( X p ) = p ∫ ∞ 0 x p -1 P ( X &gt; x )d x , X &gt; 0, p &gt; 0, we get

$$\mu _ { n } = n \int _ { 0 } ^ { 1 } z ^ { n - 1 } [ 1 - F _ { Z } ( z ; \boldsymbol \theta _ { \rho } ) ] \mathrm d z,$$

where F Z ( z ; θ ρ ) is as given in (14). Hence, from (14),

$$\mu _ { n } = \frac { 1 } { 2 } + \frac { n } { 2 } \int _ { 0 } ^ { 1 } z ^ { n - 1 } \left [ \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { \left ( \frac { \alpha _ { 1 } } { \alpha _ { 1 } } \sqrt { \frac { \pi } { \beta _ { 1 } } } - \rho \right ) } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \, w \\ + \, \frac { \sqrt { 2 } ( s \beta _ { 2 } - \beta _ { 1 } ) } { \alpha _ { 1 } \sqrt { s \beta _ { 1 } \beta _ { 2 } ( 1 - \rho ^ { 2 } ) } } \, \frac { 1 } { \alpha _ { 2 } w + \sqrt { ( \alpha _ { 2 } w ) ^ { 2 } + 4 } } \right \} \phi ( w ) d w \right ] d z.$$

In the special case that β 1 = β 2 , we have µ 1 = 1 2 (see Remark 2.1). / In general, the integral in (15) cannot be solved analytically, so the use of numerical method is essential for its respective calculation.

Table 1 shows some values of moments (with n = 1 , . . . , 10 50 100) in (15) for different choices , , of ρ in the special case when α 1 = α 2 = β 1 = β 2 = 1.

| ρ     | -0.900    | -0.750   | -0.500   | -0.250   | -0.100   | 0.000    | 0.100    | 0.250    | 0.500    | 0.750    | 0.900     |
|-------|-----------|----------|----------|----------|----------|----------|----------|----------|----------|----------|-----------|
| µ 1   | 0.500     | 0.500    | 0.500    | 0.500    | 0.500    | 0.500    | 0.500    | 0.500    | 0.500    | 0.500    | 0.500     |
| µ 2   | 0.259     | 0.269    | 0.279    | 0.284    | 0.285    | 0.286    | 0.285    | 0.284    | 0.279    | 0.269    | 0.259     |
| µ 3   | 0.139     | 0.153    | 0.168    | 0.176    | 0.178    | 0.178    | 0.178    | 0.176    | 0.168    | 0.153    | 0.139     |
| µ 4   | 0.076     | 0.092    | 0.108    | 0.116    | 0.118    | 0.119    | 0.118    | 0.116    | 0.108    | 0.092    | 0.076     |
| µ 5   | 0.043     | 0.057    | 0.072    | 0.080    | 0.082    | 0.083    | 0.082    | 0.080    | 0.072    | 0.057    | 0.043     |
| µ 6   | 0.025     | 0.037    | 0.050    | 0.057    | 0.059    | 0.060    | 0.059    | 0.057    | 0.050    | 0.037    | 0.025     |
| µ 7   | 0.015     | 0.024    | 0.036    | 0.042    | 0.044    | 0.044    | 0.044    | 0.042    | 0.036    | 0.024    | 0.015     |
| µ 8   | 0.009     | 0.016    | 0.026    | 0.032    | 0.033    | 0.034    | 0.033    | 0.032    | 0.026    | 0.016    | 0.009     |
| µ 9   | 0.005     | 0.011    | 0.019    | 0.024    | 0.026    | 0.026    | 0.026    | 0.024    | 0.019    | 0.011    | 0.005     |
| µ 10  | 0.003     | 0.008    | 0.015    | 0.019    | 0.020    | 0.020    | 0.020    | 0.019    | 0.015    | 0.008    | 0.003     |
| µ 50  | < 10 - 7  | < 10 - 5 | < 10 - 4 | < 10 - 4 | < 10 - 4 | < 10 - 4 | < 10 - 4 | < 10 - 4 | < 10 - 4 | < 10 - 5 | < 10 - 7  |
| µ 100 | < 10 - 11 | < 10 - 7 | < 10 - 6 | < 10 - 5 | < 10 - 5 | < 10 - 5 | < 10 - 5 | < 10 - 5 | < 10 - 6 | < 10 - 7 | < 10 - 11 |

Table 1: Moments values for different choices of ρ .

## 3.4 The stress-strength probability

In what follows we explicitly calculate the stress-strength probability R ≡ P ( X &lt; Y ) in the case where ( X,Y ) ⊤ ∼ BS( θ ρ ), that is, X and Y can admit correlation.

Indeed, by (7), we have

$$R = \mathbb { P } \left ( \frac { X } { Y } < 1 \right ) = 1 - \mathbb { P } \left ( \frac { X } { Y } \geqslant 1 \right ) = 1 - F _ { Z } ( z ; \boldsymbol \theta _ { \rho } ) \Big | _ { z = 1 / 2 }.$$

Hence, as z = 1 2 if and only if / s = 1, by using (14), we get

$$R = \frac { 1 } { 2 } + \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { \left ( \frac { \alpha _ { 2 } } { \alpha _ { 1 } } \sqrt { \frac { \beta _ { 2 } } { \beta _ { 1 } } } - \rho \right ) } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \, w + \frac { \sqrt { 2 } ( \beta _ { 2 } - \beta _ { 1 } ) } { \alpha _ { 1 } \sqrt { \beta _ { 1 }, \beta _ { 2 } ( 1 - \rho ^ { 2 } ) } } \, \frac { 1 } { \alpha _ { 2 } w + \sqrt { ( \alpha _ { 2 } w ) ^ { 2 } + 4 } } \right \} \phi ( w ) \text{d} w.$$

In the particular case β 1 = β 2 , from (16), the stress-strength probability R is simply written as

$$R = \frac { 1 } { 2 } + \frac { 1 } { 2 } \int _ { - \infty } ^ { \infty } \text{erf} \left \{ \frac { \left ( \frac { \alpha _ { 2 } } { \alpha _ { 1 } } - \rho \right ) } { \sqrt { 2 ( 1 - \rho ^ { 2 } ) } } \, w \right \} \phi ( w ) \text{d} w = \frac { 1 } { 2 }, \quad \forall \rho \in ( - 1, 1 ),$$

where in the last equality we have used the fact that erf( x ) and ϕ are odd and even functions, respectively.

To our knowledge, in this setting ( β 1 = β 2 ), identity R = 1 2 was known only for the special / case in which X and Y are independent (that is, ρ = 0); see, e.g. [24] and [22].

## 3.4.1 Maximum likelihood estimation

Let { Z i : i = 1 , . . . , n } be a univariate random sample of size n from the UBS( θ ρ ) distribution with PDF as in (1), and let z i be the sample observation of Z i . Then, the log-likelihood function for θ ρ = ( α , α , β 1 2 1 , β 2 , ρ ) ⊤ , α 1 &gt; 0 , α 2 &gt; 0 , β 1 &gt; 0 , β 2 &gt; 0 and ρ ∈ ( -1 1), is given by (without , additive constant)

$$\ell ( \boldsymbol \theta _ { \rho } ) \ & = \frac { n } { 1 - \rho ^ { \rho } } \left ( \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } \right ) - n \log \left ( 4 \alpha _ { 1 } \alpha _ { 2 } \pi \sqrt { 1 - \rho ^ { 2 } } \right ) - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) } \sum _ { i = 1 } ^ { n } \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) \\ & + \sum _ { i = 1 } ^ { n } \log \left \{ \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) K _ { 0 } \left ( \frac { \sqrt { \frac { \pi _ { i }, i \rho, i } { 1 - \rho ^ { 2 } } } } { 1 - \rho ^ { 2 } } \right ) + \left [ \sqrt { \frac { ( v _ { i }, i \rho, i ) s _ { i } } { \beta _ { 1 }, \beta _ { 2 } } } + \sqrt { \frac { \beta _ { 1 }, \beta _ { 2 } } { ( v _ { i }, i \rho, i ) s _ { i } } } \right ] K _ { 1 } \left ( \frac { \sqrt { \frac { \pi _ { i }, i \rho, i } { 1 - \rho ^ { 2 } } } } { 1 - \rho ^ { 2 } } \right ) \right \},$$

where

$$s _ { i } = \frac { 1 } { z _ { i } } - 1, \ \ i = 1, \dots, n,$$

$$& s _ { i } - \frac { - } { z _ { i } } - \frac { 1 } { }, \quad \imath - \frac { 1 } { }, \dots, \imath, \\ & u _ { \rho, i } \equiv \frac { s _ { i } } { \alpha _ { 1 } ^ { 2 } \beta _ { 1 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } \beta _ { 2 } } - \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { s _ { i } } { \beta _ { 1 } \beta _ { 2 } } }, \quad v _ { \rho, i } \equiv \frac { \beta _ { 1 } } { \alpha _ { 1 } ^ { 2 } s _ { i } } + \frac { \beta _ { 2 } } { \alpha _ { 2 } ^ { 2 } } - \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s _ { i } } },$$

and K λ ( x ), x &gt; 0, is the modified Bessel function of the third kind defined in (3).

The necessary conditions for the occurrence of a maximum (or a minimum) are

$$\frac { \partial \ell ( \boldsymbol \theta _ { \rho } ) } { \partial \alpha _ { 1 } } = 0, \quad \frac { \partial \ell ( \boldsymbol \theta _ { \rho } ) } { \partial \alpha _ { 2 } } = 0, \quad \frac { \partial \ell ( \boldsymbol \theta _ { \rho } ) } { \partial \beta _ { 1 } } = 0, \quad \frac { \partial \ell ( \boldsymbol \theta _ { \rho } ) } { \partial \beta _ { 2 } } = 0, \quad \frac { \partial \ell ( \boldsymbol \theta _ { \rho } ) } { \partial \rho } = 0,$$

known as the likelihood equations, with

$$\frac { \partial \ell ( \theta _ { \rho } ) } { \partial \alpha _ { 1 } } = - \frac { 2 n } { ( 1 - \rho ^ { 2 } ) \alpha _ { 1 } ^ { 3 } } - \frac { n } { \alpha _ { 1 } } + \frac { \rho } { \alpha _ { 1 } ^ { 2 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) } \sum _ { i = 1 } ^ { n } \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) + \sum _ { i = 1 } ^ { n } \frac { \frac { \partial G ( \theta _ { \rho }, s _ { i } ) } { \partial \alpha _ { 1 } } } { G ( \theta _ { \rho }, s _ { i } ) },$$

$$\frac { \partial \ell ( \theta _ { \rho } ) } { \partial \alpha _ { 2 } } = - \frac { 2 n } { ( 1 - \rho ^ { 2 } ) \alpha _ { 2 } ^ { 3 } } - \frac { n } { \alpha _ { 2 } } + \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } ^ { 2 } ( 1 - \rho ^ { 2 } ) } \sum _ { i = 1 } ^ { n } \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) + \sum _ { i = 1 } ^ { n } \frac { \frac { \partial G ( \theta _ { \rho }, s _ { i } ) } { \partial \alpha _ { 2 } } } { G ( \theta _ { \rho }, s _ { i } ) },$$

$$\frac { \partial \ell ( \theta _ { \rho } ) } { \partial \beta _ { 1 } } = - \frac { \rho } { 2 \alpha _ { 1 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) } \sum _ { i = 1 } ^ { n } \left ( \frac { 1 } { \sqrt { \beta _ { 1 } \beta _ { 2 } s _ { i } } } - \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } ^ { 3 } } } \right ) + \sum _ { i = 1 } ^ { n } \frac { \frac { \partial G ( \theta _ { \rho }, s _ { i } ) } { \partial \beta _ { 1 } } } { G ( \theta _ { \rho }, s _ { i } ) },$$

$$\frac { \partial \ell ( \theta _ { \rho } ) } { \partial \beta _ { 2 } } = - \frac { \rho } { 2 \alpha _ { 1 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) } \sum _ { i = 1 } ^ { n } \left ( \frac { 1 } { \sqrt { \beta _ { 1 } \beta _ { 2 } s _ { i } } } - \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } ^ { 3 } s _ { i } ^ { 3 } } } \right ) s _ { i } + \sum _ { i = 1 } ^ { n } \frac { \frac { \partial \log G ( \theta _ { \rho }, s _ { i } ) } { \partial \beta _ { 2 } } } { G ( \theta _ { \rho }, s _ { i } ) },$$

$$\frac { \partial \ell ( \theta _ { \rho } ) } { \partial \rho } = \frac { 2 n \rho } { ( 1 - \rho ^ { 2 } ) ^ { 2 } } \left ( \frac { 1 } { \alpha _ { 1 } ^ { 2 } } + \frac { 1 } { \alpha _ { 2 } ^ { 2 } } \right ) - \frac { n \rho } { 1 - \rho ^ { 2 } } - \frac { 1 + \rho ^ { 2 } } { \alpha _ { 1 } \alpha _ { 2 } ( 1 - \rho ^ { 2 } ) ^ { 2 } } \sum _ { i = 1 } ^ { n } \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) + \sum _ { i = 1 } ^ { n } \frac { \frac { \partial G ( \theta _ { \rho }, s _ { i } ) } { \partial \rho } } { G ( \theta _ { \rho }, s _ { i } ) },$$

where we denote

$$G ( \boldsymbol \theta _ { \rho }, s _ { i } ) \ = \left ( \sqrt { \frac { \beta _ { 2 } s _ { i } } { \beta _ { 1 } } } + \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } } \right ) K _ { 0 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) + \left [ \sqrt { \frac { ( v _ { \rho, i } / u _ { \rho, i } ) s _ { i } } { \beta _ { 1 } \beta _ { 2 } } } + \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { ( v _ { \rho, i } / u _ { \rho, i } ) s _ { i } } } \right ] K _ { 1 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ),$$

and its partial derivatives are obtained using the known fact that ∂K ν ( x /∂x ) = -K ν -1 ( x ) -( ν/x K ) ν ( x ).

Note that

$$G \left ( \theta _ { \rho }, s _ { i } \right ) \ = \beta _ { 2 } s _ { i } \, a ^ { \prime } \left ( \beta _ { 2 } s _ { i } ; \frac { 1 } { 2 }, \beta _ { 1 } \right ) K _ { 0 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) + \frac { v _ { \rho, i } } { u _ { \rho, i } } \, s _ { i } \, a ^ { \prime } \left ( \frac { v _ { \rho, i } } { u _ { \rho, i } } \, s _ { i } ; \frac { 1 } { 2 }, \beta _ { 1 } \beta _ { 2 } \right ) K _ { 1 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ),$$

where a ′ ( t, α, β ) is given in (9).

We require additional partial derivatives of G ( θ ρ , s i ). Essentially, the product and chain rules are employed for deriving the subsequent functions:

## (a) Partial derivatives with respect to α 1 :

$$\frac { \partial v _ { \rho, i } } { \partial \alpha _ { 1 } } = - \frac { 2 \beta _ { 1 } } { \alpha _ { 1 } ^ { 3 } s _ { i } } + \frac { 2 \rho } { \alpha _ { 1 } ^ { 2 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s _ { i } } }, \quad \frac { \partial u _ { \rho, i } } { \partial \alpha _ { 1 } } = - \frac { 2 s _ { i } } { \alpha _ { 1 } ^ { 3 } \beta _ { 1 } } + \frac { 2 \rho } { \alpha _ { 1 } ^ { 2 } \alpha _ { 2 } } \sqrt { \frac { s _ { i } } { \beta _ { 1 } \beta _ { 2 } } } ;$$

- (b) Partial derivatives with respect to α 2 :

$$\frac { \partial v _ { \rho, i } } { \partial \alpha _ { 2 } } = - \frac { 2 \beta _ { 2 } } { \alpha _ { 2 } ^ { 3 } } + \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } ^ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s _ { i } } }, \quad \frac { \partial u _ { \rho, i } } { \partial \alpha _ { 2 } } = - \frac { 2 } { \alpha _ { 2 } ^ { 3 } \beta _ { 2 } } + \frac { 2 \rho } { \alpha _ { 1 } \alpha _ { 2 } ^ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s _ { i } } } ;$$

- (c) Partial derivatives with respect to β 1 :

$$\frac { \partial v _ { \rho, i } } { \partial \beta _ { 1 } } = \frac { 1 } { \alpha _ { 1 } ^ { 2 } s _ { i } } - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 2 } } { \beta _ { 1 } s _ { i } } }, \quad \frac { \partial u _ { \rho, i } } { \partial \beta _ { 1 } } = - \frac { s _ { i } } { \alpha _ { 1 } ^ { 2 } \beta _ { 1 } ^ { 2 } } + \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { s _ { i } } { \beta _ { 1 } ^ { 3 } \beta _ { 2 } } } ;$$

- (d) Partial derivatives with respect to β 2 :

$$\frac { \partial v _ { \rho, i } } { \partial \beta _ { 2 } } = \frac { 1 } { \alpha _ { 2 } ^ { 2 } } - \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 1 } } { \beta _ { 2 } s _ { i } } }, \quad \frac { \partial u _ { \rho, i } } { \partial \beta _ { 2 } } = - \frac { 1 } { \alpha _ { 2 } ^ { 2 } \beta _ { 2 } ^ { 2 } } + \frac { \rho } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { s _ { i } } { \beta _ { 1 } \beta _ { 2 } ^ { 3 } } } ;$$

## (e) Partial derivatives with respect to ρ :

$$\frac { \partial v _ { \rho, i } } { \partial \rho } = - \frac { 2 } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { \beta _ { 1 } \beta _ { 2 } } { s _ { i } } }, \quad \frac { \partial u _ { \rho, i } } { \partial \rho } = - \frac { 2 } { \alpha _ { 1 } \alpha _ { 2 } } \sqrt { \frac { s _ { i } } { \beta _ { 1 } \beta _ { 2 } } }.$$

Finally, for θ k ∈ { α , α , β 1 2 1 , β 2 , ρ } , we have

$$\frac { \partial } { \partial \theta _ { k } } K _ { 0 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) = - K _ { - 1 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) \frac { \partial } { \partial \theta _ { k } } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ),$$

and

$$\frac { \partial } { \partial \theta _ { k } } K _ { 1 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) = \left [ - K _ { 0 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) - \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) ^ { - 1 } K _ { 1 } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ) \right ] \frac { \partial } { \partial \theta _ { k } } \left ( \frac { \sqrt { u _ { \rho, i } v _ { \rho, i } } } { 1 - \rho ^ { 2 } } \right ).$$

No closed-form solution to the maximization problem is available. So, the maximum likelihood estimator of θ ρ can only be obtained via numerical optimization.

## 3.5 Maximum product of spacings method

The maximum product of spacings (MPS) method was introduced by [7, 8] and [23] as a different approach to estimating parameters of continuous univariate distributions compared to maximum likelihood estimation. Let { Z i : i = 1 , . . . , n } be a univariate random sample of size n from the UBBS1( θ ρ ) distribution with PDF as in (1), and let z i be the sample observation of Z i . Then, the uniform spacings of { z i : i = 1 , . . . , n } are given by

$$\Delta _ { i } ( \theta _ { \rho } ) = F _ { Z } ( z _ { i \colon n } ; \theta _ { \rho } ) - F _ { Z } ( z _ { i - 1 \colon n } ; \theta _ { \rho } ), \ \ 1, \dots, n,$$

with F Z ( z 0: n ; θ ρ ) = 0, F Z ( z n +1: n ; θ ρ ) = 1, and ∆ ( 0 θ ρ ) + ∆ ( 1 θ ρ ) + · · · +∆ n +1 ( θ ρ ) = 1.

The maximum product of spacings estimator, ̂ θ ρ , is obtained by maximizing the geometric mean of the spacings

$$G ( \theta _ { \rho } ) = \left ( \prod _ { i = 1 } ^ { n + 1 } \Delta _ { i } ( \theta _ { \rho } ) \right ) ^ { \frac { 1 } { n + 1 } },$$

or, equivalently, maximizing the function

$$H ( \theta _ { \rho } ) = \frac { 1 } { 1 + n } \sum _ { i = 1 } ^ { n + 1 } \log ( \Delta _ { i } ( \theta _ { \rho } ) ).$$

## 4 Simulation study

In this section, we carry out a Monte Carlo simulation study for evaluating the performance of the above described maximum likelihood and maximum product of spacings methods for the UBBS1 distribution (1). The simulation scenario considers the following setup: 300 Monte Carlo replications, sample size n ∈ { 100 200 400 800 , , , } , and vector of true parameters ( α , α , β 1 2 1 , β 2 ) ⊤ =

## Algorithm 1 Generator of random numbers from a UBBS1 distribution

- 1: Generate n random numbers from the bivariate normal distribution, X = ( X , X 1 2 ) ⊤ ∼

$$\Delta _ { 2 } \left ( \mu = \left ( 0, 0 \right ) ^ { \top }, \Sigma = \left ( \begin{matrix} 1 & \rho \\ \rho & 1 \end{matrix} \right ) \right ), \, \text{as}$$

- · Generate a matrix W of size n × 2 containing 2 n random numbers from N(0 1). ,
- · Compute the decomposition Σ = Q Q ⊤ .
- · Apply the transformation X = WQ + Jµ ⊤ , where J is a column vector of ones.
- · Obtain the matrix X = ( X , X 1 2 ) ⊤ of size n × 2.
- 2: Compute n random numbers from T i ∼ BS( α , β i i ) by using

$$T _ { i } = \beta _ { i } \left ( \frac { \alpha _ { i } X _ { i } } { 2 } + \sqrt { \frac { \alpha _ { i } ^ { 2 } X _ { i } ^ { 2 } } { 4 } + 1 } \right ) ^ { 2 }, \ \ i = 1, 2.$$

- 3: Compute n random numbers from Z ∼ UBBS1( θ ρ = ( α , α , β 1 2 1 , β 2 , ρ ) ⊤ ) by using Z = T 1 T 1 + T 2 .

(0 5 . , 0 5 1 0 1 0) , , . , . ⊤ , ρ ∈ { 0 10 . , 0 25 . , 0 5 . , 0 75 . } . A generator of random numbers for the UBBS1 distribution is summarized in Algorithm 1.

The performance and recovery of the estimators are evaluated by means of relative bias (RB), and root mean square error (RMSE), which are calculated from the Monte Carlo replicates, as

$$\widehat { R B } ( \widehat { \gamma } ) \ = \ \frac { 1 } { N } \sum _ { i = 1 } ^ { N } \left | \frac { ( \widehat { \gamma } ^ { ( i ) } - \gamma ) } { \gamma } \right |, \ \widehat { R M S E } ( \widehat { \gamma } ) = \sqrt { \frac { 1 } { N } \sum _ { i = 1 } ^ { N } ( \widehat { \gamma } ^ { ( i ) } - \gamma ) ^ { 2 } },$$

where γ and ̂ γ ( ) i are the true parameter value and its i -th estimate, and N is the number of Monte Carlo replications. The steps for the Monte Carlo simulation study are as follows: i. set the values of the parameters of the UBBS1 distribution; ii. generate 300 samples of size n from the chosen model; iii. estimate the model parameters using the maximum product of spacings method for each sample; and iv. compute the empirical RB and RMSE.

The maximum likelihood and maximum product of spacings estimation results so obtained are shown in Figures 5-8, wherein the empirical RB and RMSE are both plotted. The following results are observed: i) In general, when the sample size increases, both RB and RMSE decrease for all estimators and cases; ii) in the maximum product of spacings case, we observe that the RB of ̂ α 1 , ̂ α 2 , ̂ β 1 , and ̂ β 2 tends to increase as ρ increases. In the case of ̂ ρ , we observe the opposite behavior; and iii) in the maximum likelihood case, the impact of ρ on the RB of ̂ α 1 , ̂ α 2 , ̂ β 1 , ̂ β 2 are somewhat undefined; iv) the RB for ̂ α 1 , ̂ α 2 , and ̂ ρ is smaller for the estimators based on the maximum product of spacings method compared to maximum likelihood, while for ̂ β 1 and ̂ β 2 , the RB is smaller for the estimators based on the maximum likelihood method. We observe a similar behavior in the RMSE case, except for ̂ ρ , where the maximum product of spacings method performs better.

We now compare the performance of the maximum likelihood and maximum product of spacings methods when bimodality is present. We consider θ ρ = (1 6 0 7 1 1 . , . , . , 0 9 . , 0 6) . ⊤ with n = 100, that

Figure 5: Monte Carlo simulation results for the UBBS1 distribution based on the maximum product of spacings method.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-1c5c4781-1c25-4585-888f-2ce102d46dec.jpg)

Figure 6: Monte Carlo simulation results for the UBBS1 distribution based on the maximum product of spacings method.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-1444afea-978a-4304-8874-68a62e58e27c.jpg)

is, the density exhibits a bimodal behavior; see Table 2. The results show that, with the exception of the ̂ β 2 case, the performance is better for estimates based on the MPS method; we considered other sample sizes and found similar results.

Figure 7: Monte Carlo simulation results for the UBBS1 distribution based on the maximum likelihood method.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-1e069487-5768-494c-ae5a-cf35129bec13.jpg)

Figure 8: Monte Carlo simulation results for the UBBS1 distribution based on the maximum likelihood method.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-2252738f-bf0e-45ee-b69a-b9e0646ad841.jpg)

## 5 Application

We now give two applications using real data previously analyzed in the literature. We study the validity of the new model for the data sets and show that UBBS1 distribution presents good fits in both cases.

Table 2: RB and RMSE (between brackets) from simulated data for the UBBS1 parameter estimates based on the maximum likelihood (ML) and maximum product of spacings (MPS) methods.

| Method   | ̂ α 1     | ̂ α 2     | ̂ β 1     | ̂ β 2     | ̂ ρ       |
|----------|----------|----------|----------|----------|----------|
| MPS      | 0.1523   | 0.1608   | 0.0767   | 0.3426   | 0.3074   |
|          | (0.2437) | (0.1125) | (0.0843) | (0.3084) | (0.1844) |
| ML       | 0.2721   | 0.8957   | 0.1714   | 0.0111   | 1.6970   |
|          | (0.4353) | (0.6270) | (0.1885) | (0.0100) | (1.0182) |

## 5.1 Income-consumption data

In order to evaluate the new probabilistic model, a public data set of income and consumption data was used. The data (available at https://www.bancaditalia.it/statistiche/ tematiche/indagini-famiglie-imprese/bilanci-famiglie/documentazione/ricerca/ricerca.html? min anno pubblicazione=2008&amp;max anno pubblicazione=2008, accessed on July 22, 2024)

comes from the Bank of Italy's Survey on Household Income and Wealth from the year 2008.

The income comprises payroll income, pensions and net transfers, net self-employment income, and property income. Consumption consists of expenditures on both durable and non-durable goods. Durable expenditure represents the net balance of purchased and sold goods. Non-durable expenditures include monthly expenses such as rent, food, and other essentials, as well as annual expenses.

We used the data from the survey, available as the data sets RFAM08 (income as Y variable) and RISFAM08 (expenditures as C variable), according to the Survey on Household Income and Wealth 2008 data description file. We removed from the data set lines whose income or consumption was negative or unavailable. Thus, the analyzed data set included information from 7,957 families, after removing 20 lines of the data.

The same data set was also used in [12] and [17], however, they used a copula approach to model the data and compared income and consumption through stress-strength measures of the type P ( X &lt; Y ). Our framework involves constructing a new variable U = X/ X ( + Y ). In summary, when X and Y represent expenditures and income, respectively, U &lt; 1 / 2 indicates that the family ended the year with more expenditures than income. If U &gt; 1 / 2, the opposite is true. The case U = 1 2 means expenditures and income are equal. /

Table 3 presents descriptive statistics for U = X/ X ( + Y ). The estimated parameters of UBBS1 model was ̂ θ ρ = (0 275 0 274 1 041 1 331 0 149). Figure 9 shows the fit of the UBBS1 distribution . , . , . , . , . to the density and ECDF of the data. The histogram of U shows that with high (81.5%) probability U &gt; 1 2. That means, in general, the families have income greater than expenditures. /

Table 3: Descriptive statistics for income-consumption data.

| n     |   Min. |   1st Qu. |   Median |   Mean |   3rd Qu. |   Max. |   Std. dv. |     CS |    CK |
|-------|--------|-----------|----------|--------|-----------|--------|------------|--------|-------|
| 7,957 |  0.004 |     0.514 |    0.553 |  0.557 |     0.608 |  0.935 |      0.086 | -0.419 | 3.343 |

Given that the U data has unit support, we can make a comparison between the utilization of the UBBS1 model and other families of unit distributions, such as Beta and the unit bimodal BS

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-f8c48c5f-ebba-4584-8113-9189e298bb3d.jpg)

0.0

ECDF

fitted UBBS1

0.2

0.4

0.6

0.8

1.0

x

Figure 9: Fitted UBBS1 PDFs (left) and empirical CDF (right) for income-consumption data.

(UBBS ) models, as illustrated in Table 4. ∗ To obtain the estimation results for Beta and UBBS models in Table 4, we used the package AdequacyModel of the R software and we used the function ' goodness.fit ' to obtain maximum likelihood estimates. Deciding the better distribution to model the data was oriented by AIC (Akaike information criterion) and BIC (Bayesian information criterion) criteria. In Table 4, these criteria suggest UBBS1 fits the data better than Beta and UBBS models.

Table 4: Estimated parameters for income-consumption data and model selection with AIC and BIC.

| Model   | ˆ θ                                 | ll max    | AIC        | BIC        |
|---------|-------------------------------------|-----------|------------|------------|
| UBBS1   | (0.275, 0.274, 1.041, 1.331, 0.149) | -7,062.74 | -14,135.49 | -14,170.40 |
| Beta    | (7.492, 6.390)                      | -6,877.40 | -13,750.79 | -13,736.83 |
| UBBS    | (1.726, 0.589, 9.545)               | -3,614.60 | -7,223.20  | -7,202.25  |

## 5.2 Body mass data

Biomedical measurements of 202 athletes from different sports are presented in the Australian Institute of Sport (AIS) data set. The authors in [18] proposed the UBBS model and showed its fit on the distribution of the percentage of the body fat complement of the athletes. In this subsection, we show that the new UBBS1 model is an alternative distribution to fit that data. The data were imported directly through the software R by the command:

∗ The UBBS model was proposed in [18].

CDF(x)

1.0

0.8

0.6

0.4

0.2

0.0

l

## library(sn) data(ais)

## u=1-(ais$Bfat/100)

Descriptive statistics for u (the athlete's body mass in the data set) are presented in Table 5. As the data set has unit support, the UBBS1 distribution is a candidate to model such data. Figure 10 shows the fit of UBBS1 to the data. Maximum likelihood estimates and maximum product of spacings estimates can be compared graphically. To obtain the ML estimates, we minimized the minus log-likelihood function using the optim function in the R software.

| Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   | Table 5: Descriptive statistics of athlete's body mass in the data set.   |
|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|---------------------------------------------------------------------------|
| Min.                                                                      | 1st Qu.                                                                   | Median                                                                    | Mean                                                                      | 3rd Qu.                                                                   | Max.                                                                      | Std. dv.                                                                  | CS                                                                        | CK                                                                        |
| 0.645                                                                     | 0.819                                                                     | 0.884                                                                     | 0.865                                                                     | 0.915                                                                     | 0.944                                                                     | 0.062                                                                     | -0.754                                                                    | -0.201                                                                    |

In Table 6, when comparing the UBBS1 model (with parameter estimates via product of spacings) with the Beta, UBBS, and Bimodal Beta (BBeta) models using information criteria, the BIC indicates the new model is better to fit the athlete's body mass data. However, when considering the AIC, the UBBS model should be chosen. The ML estimates for the Beta and UBBS models, presented in Table 6, were previously reported in Table 6 of [18]. For the BBeta model [25], the ML estimates were obtained using the ' optim ' function.

Table 6: Estimated parameters for the athlete's body mass data and model selection with AIC and BIC.

| Model       | ˆ θ                                          |      AIC |     BIC |
|-------------|----------------------------------------------|----------|---------|
| UBBS1 (MPS) | (0.149, 0.626, 0.296, 2.003, 0.586)          | -613.1   | -629.64 |
| UBBS1 (MLE) | (1.763, 1.106, 0.163, 1.097, 0.981)          | -619.3   | -602.75 |
| Beta        | (0.865, 0.173)                               | -585.32  | -578.7  |
| UBBS        | (0.284, 0.137, -1.353)                       | -626.64  | -616.72 |
| BBeta       | (26.1007122 4.3747393 0.9712073 -17.9150416) | -581.333 | -568.1  |

## 6 Concluding remarks

In this paper, we have introduced a new unit distribution. Its generation as a ratio of possibly dependent Birnbaum-Saunders random variables has been discussed in detail. As a more theoretical part, we have derived the cumulative distribution function to then obtain simple expressions of the moment-generating function, moments and the stress-strength probability. As practical issues, we presented maximum product of spacings estimators and applications to income-consumption and body mass data sets.

The research could be continued in several directions. One could investigate the entropy of the distribution and use the moments formulas in (15) to try to express shape characteristics like skewness and kurtosis. Instead of using the ratio Y/ X ( + Y ) one could explore other construction principles like those mentioned in the introduction to obtain a unit distribution. Finally, one could use a bivariate distribution other than the Birnbaum-Saunders to construct a unit distribution.

Figure 10: Plots for the athlete's body mass. On top, histogram and fitted UBBS1 PDFs (ML estimates on left and Maximum product of spacings estimates on right); on bottom, Empirical CDF and fitted UBBS1 CDF.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00100v2-figure-64555cb8-5b28-41ce-8d1d-3a82f1338537.jpg)

Acknowledgements This study was financed in part by the Coordena¸ c˜o de Aperfei¸ coamento a de Pessoal de N´ ıvel Superior - Brasil (CAPES) - Finance Code 001.

Disclosure statement There are no conflicts of interest to disclose.

## References

| [1]   | Abramowitz, M., Stegun, I. A.: (eds.): Handbook of Mathematical Functions, 9th edn. Dover, New York (1972).                                                                                           |
|-------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [2]   | Afifi, A. Z., Nassar, M., Kumar, D., Cordeiro, G.M.: A New Unit Distribution: Properties, Inference, and Applications. Electronic Journal of Applied Statistical Analysis. 15 (02), 460-484 (2022).   |
| [3]   | Bekker, A., Roux, J., Pham-Gia, T.: The type I distribution of the ratio of independent 'Weibullized' generalized beta-prime variables. Statistical Papers., 50 , 323-338 (2009).                     |
| [4]   | Birnbaum, Z. W., Saunders, S. C.: A new family of life distributions. Journal of Applied Probability. 6 (2), 319-327 (1969).                                                                          |
| [5]   | Birnbaum, Z. W., Saunders, S. C.: Estimation for a family of life distributions with applications to fatigue. Journal of Applied Probability. 6 (2), 328?347 (1969).                                  |
| [6]   | Bourguignon, M., Saulo, H., Gallardo, D.: Parametric quantile beta regression model. Inter- national Statistical Review. 92 , 106-129 (2024).                                                         |
| [7]   | Cheng, R. C. H., Amin, N. A. K.: Maximum product-of-spacings estimation with applications to the lognormal distribution, Technical Report, Department of Mathematics, University of Wales (1979).     |
| [8]   | Cheng, R. C. H., Amin, N. A. K.: Estimating parameters in continuous univariate distributions with a shifted origin. Journal of the Royal Statistical Society. Series B, 3 , 394-403 (1983).          |
| [9]   | Condino, F., Domma, F.: Unit Distributions: A General Framework, Some Special Cases, and the Regression Unit-Dagum Models. Mathematics. 11 , 2888 (2023).                                             |
| [10]  | Cort´s, e M. A., Elal-Olivero, D., Olivares-Pacheco, J. F.: A New Class of Distributions Gen- erated by the Extended Bimodal-Normal Distribution. Journal of Probability and Statistics. 1-10 (2018). |
| [11]  | Desmond, A.: Stochastic models of failure in random environments. Canadian Journal of Statistics. 13 (3), 171-183 (1985).                                                                             |
| [12]  | Domma, F., Giordano, S.: A stress-strength model with dependent variables to measure house- hold financial fragility. Statistical Methods & Applications. 21 , 375-389 (2012).                        |
| [13]  | Johnson, N. L., Kotz, S., Balakrishnan, N.: Continuous Univariate Distributions - Vol. 2. John Wiley & Sons, New York (1995).                                                                         |
| [14]  | Jorgensen, B.: Statistical Properties of the Generalized Inverse Gaussian Distribution, Lecture Notes in Statistics, vol. 9. Springer, New York (1982).                                               |
| [15]  | Kundu, D., Gupta, R. C.: On Bivariate Birnbaum-Saunders Distribution. American Journal of Mathematical and Management Sciences. 36 (1), 21-33 (2017).                                                 |

| [16]   | Kundu, D., Balakrishnan, N., Jamalizadeh, A.: Bivariate Birnbaum-Saunders distribution and associated inference. Journal of Multivariate Analysis. 101 , 113-125 (2010).                                                                                      |
|--------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [17]   | Lima, R., Quintino, F., da Fonseca, T., Ozelim, L., Rathie, P. N., Saulo, H.: Assessing the Impact of Copula Selection on Reliability Measures of Type P ( X < Y ) with Generalized Extreme Value Marginals. Modelling. 5 (1), 180-200 (2024).                |
| [18]   | Mart´nez-Fl´rez, ı o G., Olmos, N., Venegas, O.: Unit-bimodal Birnbaum-Saunders distribution with applications. Communications in Statistics-Simulation and Computation. 53 (5), 2173- 2192 (2024).                                                           |
| [19]   | Mazucheli, J., Menezes, A.F.B., Dey, S.: The unit-Birnbaum- Saunders distribution with ap- plications. Chilean Journal of Statistics, 9 (1), 47-57 (2018).                                                                                                    |
| [20]   | Mazucheli, J., Leiva, V., Alves, B., Menezes, A.F.B.: A New Quantile Regression for Modeling Bounded Data under a Unit Birnbaum- Saunders Distribution with Applications in Medicine and Politics. Symmetry, 13 , 682 (2021).                                 |
| [21]   | Teimouri, M., Hosseini, S. M., Nadarajah, S.: Ratios of Birnbaum-Saunders Random Variables. Quality Technology & Quantitative Management. 10 , 457-481 (2013).                                                                                                |
| [22]   | Quintino, F. S., Ozelim, L. C. d. S. M., Fonseca, T. A. d., Rathie, P. N.: Stress-Strength Reliability of the Type P ( X < Y ) for Birnbaum-Saunders Components: A General Result, Simulations and Real Data Set Applications. Modelling. 5 , 223-237 (2024). |
| [23]   | Ranneby, B.: The maximum spacing method. An estimation method related to the maximum likelihood method. Scandinavian Journal of Statistics. 11 , 93-112 (1984).                                                                                               |
| [24]   | Saulo, H., Le˜o, a J., Vila, R., Leiva, V., Tomazella, V.: On mean-based bivariate Birnbaum- Saunders distributions: Properties, inference and application. Communications in Statistics - Theory and Methods. 49 (24), 6032-6056 (2019).                     |
| [25]   | Vila, R., Alfaia, L., Menezes, A. F. B., Cankaya, ¸ M. N., Bourguignon, M.: AModel for Bimodal Rates and Proportions. Journal of Applied Statistics. 51 (4), 664-681 (2024).                                                                                  |
| [26]   | Vila, R., Balakrishnan, N, Bourguignon, M.: On the distribution of a random variable involved in an independent ratio. Communications in Statistics-Theory and Methods. (2024). https: //doi.org/10.1080/03610926.2024.2340608.                               |
| [27]   | Vila, R., Balakrishnan, N., Saulo, H., Z¨rnig, o P.: Unit-log-symmetric models: characterization, statistical properties and their applications to analyzing an internet access data. Quality and Quantity. 58 , 4779-4806 (2024).                            |
| [28]   | Z¨rnig, o P.: A new type of skew uniform distribution, Asian Journal of Statistical Sciences. 3 , 27-39 (2023).                                                                                                                                               |
Roberto Vila Helton Saulo Felipe Quintino Peter Zörnig
A new unit-bimodal distribution based on correlated Birnbaum-Saunders random variables
In this paper, we propose a new distribution over the unit interval which can
be characterized as a ratio of the type $Z=Y/(X+Y)$ where $X$ and $Y$ are two
correlated Birnbaum-Saunders random variables. The density of $Z$ may be
unimodal or bimodal. Simple expressions for the cumulative distribution
function, moment-generating function and moments are obtained. Moreover, the
stress-strength probability between $X$ and $Y$ is calculated explicitly in the
symmetric case, that is, when the respective scale parameters are equal. Two
applications of the ratio distribution are discussed.