Flexible Log-Linear Birnbaum–Saunders Model

mathematics Article Flexible Log-Linear Birnbaum–Saunders Model Guillermo Martínez-Flórez 1,2 , Inmaculada Barranco-Chamorro 3,* and Héctor W. Gómez 4   Citation: Martínez-Flórez, G.; Barranco-Chamorro, I.; Gómez, H.W. Flexible Log-Linear Birnbaum–Saunders Model. Mathematics 2021,9, 1188. https:// doi.org/10.3390/math9111188 Academic Editors: Carlos Agra Coelho and Tatjana von Rosen Received: 19 April 2021 Accepted: 17 May 2021 Published: 24 May 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Departamento de Matemáticas y Estadística, Facultad de Ciencias Básicas, Universidad de Córdoba, Córdoba 230027, Colombia; [email protected] 2Programa de Pós-Graduação em Modelagem e Métodos Quantitativos, Universidade Federal do Ceará, Fortaleza 60020-181, Brazil 3Departamento de Estadística e I.O., Facultad de Matemáticas, Universidad de Sevilla, 41012 Sevilla, Spain 4Departamento de Matemáticas, Facultad de Ciencias Básicas, Universidad de Antofagasta, Antofagasta 1240000, Chile; hector[email protected] *Correspondence: chamorr[email protected] Abstract: Rieck and Nedelman (1991) introduced the sinh-normal distribution. This model was built as a transformation of a N(0,1) distribution. In this paper, a generalization based on a flexible skew normal distribution is introduced. In this way, a more general model is obtained that can describe a range of asymmetric, unimodal and bimodal situations. The paper is divided into two parts. First, the properties of this new model, called flexible sinh-normal distribution, are obtained. In the second part, the flexible sinh-normal distribution is related to flexible Birnbaum–Saunders, introduced by Martínez-Flórez et al. (2019), to propose a log-linear model for lifetime data. Applications to real datasets are included to illustrate our findings. Keywords: flexible Birnbaum–Saunders distribution; flexible Sinh-Normal distribution; lifetime regression model; log Birnbaum–Saunders regression model 1. Introduction In this paper, a generalization of the sinh-normal distribution introduced by Rieck and Nedelman [ 1 ] is proposed. Recall that the sinh-normal model is defined as a transformation, h(·), of a standard normal distribution Y=h(Z)where Z∼N(0, 1). (1) Rieck and Nedelman [ 1 ] showed that the probability density function (pdf) of Y is symmetric and can be unimodal or bimodal. The most relevant application of (1) is that the sinh-normal distribution can be used to propose a log-linear model for lifetime data distributed as a Birnbaum–Saunders (see [ 1 ] for details). This point is of crucial importance since the Birnbaum–Saunders distribution is used in a variety of situations, such as to model lifetime, economics and environmental data. In these applications, departures of the BS model are often found. In recent years, some improvements have been introduced to deal with this problem. In this sense, we highlight the flexible Birnbaum–Saunders model, proposed by Matínez-Flórez et al. [ 2 ], where two parameters, δ∈R , λ> 0, are added to the usual Birnbaum–Saunders model, in such a way that λ controls asymmetry (skewness) and δis a shape parameter related to unimodality/bimodality of our proposal. Throughout this paper, the next points are addressed. First, details about the precedents, in which our proposal is based on, are given in Section 2. Second, the flexible sinh-normal distribution is defined as result of applying (1) to Z distributed as a flexible skew-normal distribution, Z∼FSN(δ , λ) with δ∈R and λ> 0. Recall that the FSN(δ , λ) model was introduced by Gómez et al. [ 3 ], as a generalization of Azzalini skew-normal model [ 4 ], and it can deal with unimodal and bimodal situations. Mathematics 2021,9, 1188. https://doi.org/10.3390/math9111188 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 1188 2 of 23 The properties of the flexible sinh-normal distribution are given in Section 3. These are explicit expressions for the pdf and cumulative distribution function (cdf), shape and modes of the distribution, location and scale changes, approximations, moments and quantiles, among others. All of them are based on the fact that the flexible sinh-normal is obtained as a given function of Z∼FSN(δ , λ) . The obtained properties are compared to those in the Rieck and Nedelman model. In addition, we discuss in which sense our proposal is more general than others previously introduced in literature. The relationship between the flexible sinh-normal distribution and the flexible Birnbaum–Saunders via the exponential (or neperian logarithm), which allows us to build a regression model for lifetime data, is also proven there. Section 4is devoted to maximum likelihood estimation for the parameters in the flexible sinh-normal distribution. Third, the flexible Log-linear Birnbaum–Saunders regression model is studied in Section 5. This can be used for non-negative lifetime variables following the flexible Birnbaum–Saunders model given in [ 2 ]. The model is built, particular cases of interest are cited, scores and maximum likelihood equations are given. Applications to real datasets are included in Section 6. In Section 6.1, a positively skewed environmental dataset is considered. There, it is proven that the flexible sinhnormal model provides a better fit than other precedent skew models. In Section 6.2, a symmetric bimodal dataset is considered. Only a few models in the literature exhibit both properties. They are considered along with a mixture of normal distributions. It is shown that the flexible sinh-normal submodel with skewness parameter equal to zero provides a better fit than the other ones. Section 6.3 deals with failure time data in a fatigue test, where the values of stress are considered as known covariate. It is shown that the flexible log-linear Birnbaum–Saunders provides a superior fit to the log-Birnbaum–Saunders [ 1 ] and the log-skew-Birnbaum–Saunders model. A simulation study is presented in Section 7. A final discussion is given as Section 8. 2. Materials and Methods For completeness, details about the most relevant precedents in which our proposal is based are introduced next. Flexible skew-normal distribution. The flexible skew-normal (FSN) model was obtained and studied in detail by [ 5 ]. They proved that the FSN model can be bimodal for certain values of δ . A random variable (rv) Zfollows a FSN distribution, Z∼FSN(δ,λ), if its pdf is given by f(z;δ,λ) = cδφ(|z|+δ)Φ(λz),z∈R,λ,δ∈R, (2) where φ and Φ are the pdf and cdf of the N( 0, 1 ) distribution, respectively, and c−1 δ=1−Φ(δ). Particular cases of interest in the flexible skew-normal are the following: - Taking δ= 0 in (2), the skew-normal distribution proposed by Azzalini [ 4 ] is obtained. - If λ= 0, then (2) reduces to the flexible normal distribution introduced by MartínezFlórez et al. [6]. - If δ=λ=0, then Z∼N(0, 1). Sinh-normal or log-Birnbaum–Saunders distribution. Rieck and Nedelman [ 1 ] developed the sinh-normal (SHN) distribution, which is given as the following transformation of a standard normal distribution Y=arcsinhγ 2Zη+ξwhere Z∼N(0, 1), (3) γ> 0 is a shape parameter, ξ∈R is a location parameter and η> 0 is a scale parameter. (3) is denoted as Y∼SHN(γ,ξ,η). Mathematics 2021,9, 1188 3 of 23 Proposition 1 (Properties of sinh-normal distribution).Let Y ∼SHN(γ,ξ,η). Then, 1. Y is symmetric about the location parameter ξ. 2. The pdf of Y is unimodal for γ≤2and bimodal for γ>2. 3. The mean and variance of Y are E[Y] = ξ and Var[Y] = η2v(γ) , where v(γ) is the variance when the scale parameter is equal to one ( η= 1). There is no closed-expression for v(γ) , but Rieck and Nedelman [ 1 ] provided asymptotic approximations for small and large values of γ . 4. Let Yγ∼SHN(γ , ξ , η) . Then, Uγ=2(Yγ−ξ) γη converges in law to a N( 0, 1 ) distribution as γ−→ 0+. The next lemma considers the particular case of a sinh-normal distribution in which the scale parameter is equal to 2, η= 2. This model is related to the Birnbaum–Saunders (BS) distribution. This result is the basis of the use of this distribution as a regression model, as shown in Section 5. Lemma 1. 1. If Y ∼SHN(γ,ξ, 2), then T =exp(Y)∼BS(γ, exp(ξ)). 2. Reciprocally, if T∼BS(α , β) , then Y=log T follows a sinh-normal distribution with shape parameter α , location parameter ξ=log β and scale parameter η= 2, that is Y∼ SHN(α, log β, 2). Proof. It can be seen in [1] (Theorem 1.1). Due to the result given in Lemma 1, the sinh-normal distribution is also named log-Birnbaum–Saunders distribution. 3. Results Flexible sinh-normal distribution. This model is obtained by applying the transformation introduced in (3) to Z∼FSN(δ,λ). The new model is called the flexible sinh-normal (FSHN) distribution. It is denoted as Y∼ FSHN(γ,ξ,η,δ,λ). By applying the change of variable formula, the pdf of Y,fY, is fY(y) = fZ(aγ,ξ,η(y))a0 γ,ξ,η(y)(4) where fZ denotes the pdf of Z∼FSN(δ , λ) given in (2), aγ,ξ,η(y) and a0 γ,ξ,η(y) = ∂ ∂yaγ,ξ,η(y)are aγ,ξ,η(y) = 2 γsinhy−ξ η,a0 γ,ξ,η(y) = 2 ηγcoshy−ξ η. (5) Explicitly, (4) can be written as fY(y) = cδ 2 ηγcoshy−ξ ηφ2 γsinhy−ξ η+δΦλ2 γsinhy−ξ η, (6) with cδ= (1−Φ(δ))−1. The following are particular cases of interest in the FSHN model: - If δ=0, then the FSHN model reduces to the skew sinh-normal distribution studied by Leiva et al. [7]. - If λ= 0, then a symmetric model denoted by FSHNλ=0(γ , ξ , η , δ) is obtained, which allows us to model symmetric bimodal data. - If δ=λ= 0, then the FSHN model reduces to the sinh-normal distribution introduced by Rieck and Nedelman [1]. Next, some results dealing with the inverse transformation to that introduced in (3) are given. Mathematics 2021,9, 1188 4 of 23 Proposition 2. Let Y ∼FSHN(γ,ξ,η,δ,λ). Then, aγ,ξ,η(Y) = 2 γsinhY−ξ η∼FSN(δ,λ). Proof. It is obtained by considering Z=aγ,ξ,η(Y) and applying the inverse transformation technique to Y=h(Z) = arcsinhγ 2Zη+ξ, with Z∼FSN(δ,λ). The following are particular cases of interest in Proposition 2: 1. If λ= 0, then Z=aγ,ξ,η(Y) follows the f lexible −normal(δ) model studied by Castillo et al. [8]. 2. If δ=0, then Z=aγ,ξ,η(Y)∼Skew −Normal(λ). 3. If λ=δ=0, then Z=aγ,ξ,η(Y)∼N(0, 1). Next, an explicit expression for the cumulative distribution function (cdf) of a FSHN distribution is given. Theorem 1. Let Y ∼FSHN(γ,ξ,η,δ,λ). Then, the cdf of Y is FY(y) =                cδΦBNλλδ √1+λ2,a(y)−δ,if y <ξ cδhΦBNλλδ √1+λ2,−δ+ΦBNλ−λδ √1+λ2,a(y) + δ−ΦBNλ−λδ √1+λ2,δi, if y ≥ξ, (7) where a(y) = aγ,ξ,η(y) is defined in (5), ΦBNλ(· , ·) denotes the cdf of a bivariate normal distribution, with mean vector µ0= (0, 0)and covariance matrix Ωλ=1ρλ ρλ1where ρλ=−λ √1+λ2. (8) Proof. Since Y∼FSHN(γ,ξ,η,δ,λ) Y=arcsinhγ 2Zη+ξ, with Z∼FSN(δ,λ). Taking into account that arcsinh(x) = logx+√x2+1 is a monotonically increasing function of x, we have the following relationship between the cdf’s de Yand Z: FY(y) = FZ(aγ,ξ,η(y)) (9) with aγ,ξ,η(y) = 2 γsinhy−ξ η. By applying the expression for the cdf of Z∼FSN(δ , λ) given in [ 2 ] we have that, for aγ,ξ,η(y)<0, or, equivalently, for y<ξ FZ(aγ,ξ,η(y)) = cδΦBNλλδ √1+λ2,a(y)−δ, and for aγ,ξ,η(y)≥0, or, equivalently, for y≥ξ FZ(aγ,ξ,η(y)) = cδΦBNλλδ √1+λ2,−δ+ΦBNλ−λδ √1+λ2,a(y) + δ−ΦBNλ−λδ √1+λ2,δ Mathematics 2021,9, 1188 5 of 23 where ΦBNλ(· , ·) denotes the cdf of a bivariate normal distribution, with mean vector µ0= (0, 0)and covariance matrix Ωλgiven in (8). Thus, the proposed result holds Next, some particular cases of interest for the cdf in the FSHN model are discussed. Corollary 1. Let Y ∼FSHN(γ,ξ,η,δ,λ). 1. If λ=0, then the cdf of Y reduces to FY(y) =    cδ 2Φ(a(y)−δ),if y <ξ cδ 2{Φ(a(y) + δ) + 1−2Φ(δ)},if y ≥ξ. (10) 2. If δ=0, then the cdf of Y reduces to FY(y) = 2ΦBNλ(0, a(y)),for y ∈R. Proof. 1. If λ= 0, then ρλ= 0, where ρλ is given in (8). Recall that, in the bivariate normal distribution, uncorrelation implies independence, therefore ΦBNλ=0(x1,x2) = Φ(x1)Φ(x2),∀(x1,x2). Taking into account that Φ( 0 ) = 1 / 2 and Φ(−δ) = 1 −Φ(δ) , we have that (7) reduces to (10). 2. If δ=0, then we have the cdf of the skew-sinh-normal distribution. The result in Corollary 1for λ= 0 corresponds to the model studied by Olmos et al. [9] , whereas the result for δ= 0 is a particular case of models studied by Vilca-Labra and Leiva-Sanchez [10]. Shape of fY(·). Proposition 3. Let Y∼FSHN(γ , ξ , η , δ , λ) .The pdf given in (16) can be bimodal. The modes would be given as the solution of the following non-linear equations 1. y∗ 1<ξsolution of δ+λφ(λs(y)) Φ(λs(y)) +s(y) c2(y)=s(y), (11) 2. y∗ 2>ξsolution of −δ+λφ(λs(y)) Φ(λs(y)) +s(y) c2(y)=s(y), (12) with s(y) = 2 γsinhy−ξ ηand c(y) = 2 γcoshy−ξ η. Proof. To obtain the maximums of fY , we take the natural logarithm of fY , except additive terms not depending of y. Thus, let us consider g(y) = ln(c(y)) + ln φ(|s(y)|+δ)+ln Φ(λs(y)) where s(y) = 2 γsinhy−ξ η and c(y) = 2 γcoshy−ξ η . Taking the first derivative of g(y) with respect to y,g0(y) = 0, after straightforward calculations, (11) and (12) are obtained. Proposition 4. Let Y∼FSHN(γ , ξ , η , δ , λ) with λ= 0. Then, the pdf of Y is symmetrical with respect to ξand can be bimodal for δ<0. Proof. Symmetry follows from the fact that cosh(·) and φ(|z|+δ) are even functions. Bimodality for δ<0 follows from Proposition 3. Proposition 5. Let Y ∼FSHN(γ,ξ,η,δ,λ). Then, the pdf of Y is non-differentiable at y =ξ. Mathematics 2021,9, 1188 6 of 23 Proof. This result follows from the fact that, if y=ξ , then s(ξ) = 2 γsinh( 0 ) = 0 and the non-differentiability of the absolute value function at zero. Some illustrative plots for the pdfs in the FSHN model are given in Figures 1–3. In Figures 1and 2, we mainly focus on the study of the shape parameters, δ and γ , so in every figure a positive or negative value of δ is fixed, and a range of possible values is considered for γ (less and greater than two), since our aim is to compare with the sinhnormal distribution, which is strongly unimodal for γ≤ 2 and bimodal for γ> 2 [ 1 ]. In this way, the effect of δ and/or γ can be appreciated. Moreover, positive and negative values of λ are considered to assess the effect of the skewness parameter. In Figure 3, λ= 0 is fixed and therefore the symmetric submodel is plotted. Additional details are given next for every figure. (a) −4 −2 0 2 4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 delta=−2.75, lambda=−0.25, 0.25 y f(y) −4 −2 0 2 4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y −4 −2 0 2 4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y −4 −2 0 2 4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y gamma=5.5 gamma=3.5 gamma=1.5 gamma=0.75 (b) −2 −1 0 1 2 3 0.0 0.5 1.0 1.5 2.0 delta=1.5, lambda=1 y f(y) −2 −1 0 1 2 3 0.0 0.5 1.0 1.5 2.0 y −2 −1 0 1 2 3 0.0 0.5 1.0 1.5 2.0 y −2 −1 0 1 2 3 0.0 0.5 1.0 1.5 2.0 y gamma=5.5 gamma=3.5 gamma=1.5 gamma=0.75 Figure 1. FSHN pdfs: ( a ) FSHN (5.5, 0, 1, − 2.75, − 0.25) (black), FSHN (3.5, 0, 1, − 2.75, − 0.25) (blue), FSHN (1.5, 0, 1, − 2.75, 0.25) (orange) and FSHN (0.75, 0, 1, − 2.75, 0.25) (red); and (b)FSHN (5.5, 0, 1, 1.5,1) (black), FSHN (3.5, 0, 1, 1.5,1) (blue), FSHN (1.5, 0, 1, 1.5,1) (orange) and FSHN (0.75, 0, 1, 1.5,1) (red). In Figure 1, without loss of generality for our purposes, the location and scale parameters are taken as ξ= 0 and η= 1. We intend to highlight the effect of δ on the bimodality of FSHN model. In Figure 1a, small values are considered for the skewness parameter, λ=− 0.25 and 0.25, and a clearly negative value of δ=− 2.75. All pdfs are bimodal. Note that the parameter γ∈ { 0.75, 1.5, 3.5, 5.5 } and even for values of γ≤ 2 the pdf is bimodal. These plots suggest that a negative value of δ induces bimodality and the effect of γ is attenuated. It is also possible to appreciate the effect of the sign of skewness parameter λ in Figure 1a; note that black and blue plots correspond to negative λ , whereas orange and red ones correspond to positive λ . On the other hand, Figure 1b, a positive value of δ= 1.5 is considered ( λ= 1). We highlight that the distribution is unimodal even for values of γ>2. In Figure 2, a small value of δ is fixed ( δ= 0.2 in Figure 2a and δ=− 0.2 in Figure 2b); the plots suggest that, in this case, the bimodality depends on the value of γ . However, we must also be conscious that, in these plots λ= 0.4, for large values of λ unimodal distributions are nearly always obtained for any values of γand δ. Mathematics 2021,9, 1188 7 of 23 (a) −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 delta=0.2, lambda=0.4 y f(y) −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y gamma=9.5 gamma=5.5 gamma=3 gamma=1.5 (b) −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 delta=−0.2, lambda=0.4 y f(y) −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y gamma=9.5 gamma=5.5 gamma=3 gamma=1.5 Figure 2. FSHN pdfs: ( a ) FSHN (9.5, 0, 1, 0.2,0.4) (black), FSHN (5.5, 0, 1, 0.2,0.4) (blue), FSHN (3, 0, 1, 0.2,0.4) (orange) and FSHN (1.5, 0, 1, 0.2,0.4) (red); and In ( b ) FSHN (9.5, 0, 1, − 0.2,0.4) (black), FSHN (5.5, 0, 1, − 0.2,0.4) (blue), FSHN (3, 0, 1, − 0.2,0.4) (orange) and FSHN (1.5, 0, 1, −0.2,0.4) (red). In Figure 3, the skewness parameter λ= 0 is fixed, therefore symmetric pdfs are always obtained. In Figure 3a, a negative value of δ=− 0.75 is fixed and all plots are bimodal (even for values of the parameter γ< 2 in which the sinh-normal distribution is unimodal). In Figure 3b, δ= 0.75 is fixed and unimodal symmetric distributions are obtained. (a) −5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 lambda=0, delta=−0.75 y −5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y −5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y −5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y gamma=4.5 gamma=2.5 gamma=1.5 gamma=1.25 (b) −5 0 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 lambda=0, delta=0.75 y −5 0 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 y −5 0 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 y −5 0 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 y gamma=4.5 gamma=2.5 gamma=1.5 gamma=1.25 Figure 3. FSHN pdfs: ( a ) FSHN (4.5, 0, 2.5, − 0.75,0) (black), FSHN (2.5, 0, 2.5, − 0.75, 0) (blue), FSHN (1.5, 0, 2.5, − 0.75, 0) (orange) and FSHN (1.25, 0, 2.5, − 0.75, 0) (red); and ( b ) FSHN (4.5, 0, 2.5, 0.75,0) (black), FSHN (2.5, 0, 2.5, 0.75, 0) (blue), FSHN (1.5, 0, 2.5, 0.75, 0) (orange) and FSHN (1.25, 0, 2.5, 0.75, 0) (red). Mathematics 2021,9, 1188 8 of 23 Next, we show that the FSHN model is closed under location and scale changes. Proposition 6. Let Y ∼FSHN(γ,ξ,η,δ,λ). Then, 1. aY ∼FSHN(γ,aξ,aη,δ,λ), a >0. 2. Y +b∼FSHN(γ,ξ+b,η,δ,λ),∀b∈R. 3. −Y∼FSHN(γ,−ξ,η,δ,−λ). Proof. 1. Let Z=aY with a> 0. The result proposed follows from the fact that fZ(z) = 1 afYz a. 2. Let Z=Y+b. Then, fZ(z) = fY(z−b). 3. Let Z=−Y . Then, fZ(z) = fY(−z) . Taking into account that cosh(·) is an even function, cosh(−x) = cosh(x) 2 γcosh−z−ξ η=2 γcoshz+ξ η=2 γcoshz−(−ξ) η, and, since sinh(·)is an odd function, sinh(−x) = −sinh(x) 2 γsinh−z−ξ η=−2 γsinhz+ξ η=−2 γsinhz−(−ξ) η, and the result proposed is obtained. Corollary 2. Let Y ∼FSHN(γ,ξ,η,δ,λ). Then, Y−ξ η∼FSHN(γ, 0, 1, δ,λ). Next, the limit behavior of FSHN model when γapproaches to 0+is studied. Proposition 7. Let Y∼FSHN(γ , ξ , η , δ , λ) . Then, Uγ=2(Y−ξ) γη converges in distribution to the flexible skew-normal FSN(δ,λ)distribution as γ−→ 0+. Proof. Let Uγ=2(Y−ξ) γη . By applying Proposition 6, Uγ=2(Y−ξ) γη ∼FSHN(γ , 0, 2 γ , δ , λ) . Therefore, its pdf is fUγ(u) = cδcoshγu 2φ2 γsinhγu 2+δΦλ2 γsinhγu 2,u∈R. (13) Considering that lim γ→0+coshγu 2=1, by L’Hôpital–Bernoulli rule, lim γ→0+ 2 γsinhγu 2=ulim γ→0+coshγu 2=u. Since φ(·)and Φ(·)are continuous differentiable functions, fUγ(u)−→ cδφ(|u|+δ)Φ(λu),u∈Ras γ→0+, which is the pdf of a flexible skew-normal (δ , λ) distribution. That is, the result proposed. Mathematics 2021,9, 1188 9 of 23 Remark 1. Given Y∼FSHN(γ , ξ , η , δ , λ) with γ close to zero, Proposition 7allows us to approximate the distribution of Y , properly normalized, by a flexible skew-normal (δ , λ) distribution. In addition, note that this result is analogous to the property of the sinh-normal distribution given in Proposition 1. Next, it is shown that the pth quantile of Y can be obtained from the pth quantile of Z∼FSN(δ,λ). Theorem 2. Let Y∼FSHN(γ , ξ , η , δ , λ) . Then, the p th quantile of Y , yp , with 0 <p< 1is given by yp=ξ+ηarcsinhγ 2zp(14) where zpdenotes the pth quantile of Z ∼FSN(δ,λ). Proof. It follows from (3) and the fact that the arcsinh function is monotonically increasing. Submodel with λ=0. For λ=0, ypcan be obtained from the quantiles of the N(0, 1)distribution. Corollary 3. If λ=0, then yp=     ξ+ηarcsinhγ 2δ+Φ−1(2p/cδ),if p <0.5 ξ+ηarcsinhnγ 2h−δ+Φ−1Φ(δ) + (2p−1) cδio,if p ≥0.5. (15) As it can be seen in [ 1 ], recall that the log-Birnbaum–Saunders regression model is based on the relationship between the Birnbaum–Saunders and the sinh-normal distribution. In our proposal, a regression model is introduced next that is based on the relationship between the flexible Birnbaum–Saunders, studied in [ 2 ] and the FSHN distribution introduced in this paper. First, we recall the expression of the pdf for the flexible Birnbaum–Saunders. Flexible Birnbaum–Saunders. Based on the flexible skew-normal model defined in (2), Martínez-Flórez et al. [ 2 ] proposed the flexible Birnbaum–Saunders (FBS) distribution, T∼FBS(α , β , δ , λ) , whose pdf is given by f(t;α,β,δ,λ) = t−3/2(t+β) 2αβ1/2(1−Φ(δ))φ(|at|+δ)Φ(λat),t>0, (16) with at at=at(α,β) = 1 α st β−rβ t!, (17) α>0, β>0, δ∈R,λ∈Rand φ(·)and Φ(·)the pdf and cdf of a N(0, 1), respectively. In (16), δ and λ are shape parameters. We highlight that λ is a parameter that controls asymmetry (or skewness) and δ is a shape parameter related to bimodality of FBS model. Specifically, it can be seen in [ 2 ] that (16) can be bimodal for some combinations of δ and λ parameters. The following are particular cases of interest in (16): - If λ=0, then the model introduced by Olmos et al. [9] is obtained. - If δ=0, then the skew-Birnbaum–Saunders is obtained. - If λ=δ=0, then the Birnbaum–Saunders is obtained. Theorem 3. Let Y ∼FSHN(γ,ξ, 2, δ,λ). Then, T =exp(Y)∼FBS(γ,exp(ξ),δ,λ). Mathematics 2021,9, 1188 16 of 23 Since SHN and FSHNλ=0are nested models, we can test H0:δ=0 versus H1:δ6=0, which is equivalent to compare SHN versus FSHNλ=0 . Taking into account that the Fisher information matrix is non-singular at (δ , λ) = ( 0, 0 ) , we can consider the likelihood ratio statistic Λ1=LSHN(b γ,b ξ,b η) LFSHN(b γ,b ξ,b η,b δ). It is obtained that − 2 log(Λ1) = 15.868, which is greater than the χ2 1,0.05 = 3.84. Thus, H0 is rejected, leading to the conclusion that the proposed FSHNλ=0 model fits better than the SHN considered in this case. In Figure 5, the histogram of bebe weights is plotted along with the proposed bimodal models. There, the good fit provided by our proposal can be checked. 6.3. Illustration 3 Dataset 3 consists of 40 independent observations, which correspond to the failure time, T , for hardened steel specimens in a rolling contact fatigue test. The observations were taken at each of four values of contact stress, which is the covariate x in our proposal. This dataset can be seen in the work of Chan et al. [ 24 ]. Let Yi=log Ti and consider the regression model Yi=β0+β1log xi+ei,i=1, . . . , 40. (36) We fit the log-Birnbaum–Saunders, the log-skew-Birnbaum–Saunders proposed by Lemonte et al. [ 13 ] and the log-flexible Birnbaum–Saunders introduced in this paper to this dataset. The estimated parameters along with the AIC and BIC are given in Table 4. According to AIC and BIC, log-flexible Birnbaum–Saunders provides a better fit than the other ones, so this model would be preferred. Table 4. MLEs in LBS, LSBS and LFBS models. Parameters LBS LSBS LFBS γ1.279 2.011 1.988 (0.143) (0.313) (0.527) β00.097 −0.961 −2.053 (0.170) (0.166) (0.530) β∗ 00.165 −0.392 β1−14.116 −13.870 −13.812 (1.571) (1.602) (1.248) δ−1.565 (0.430) λ−0.932 2.073 (0.174) (0.910) AIC 129.23 125.36 121.885 BIC 134.30 132.11 130.32 7. Simulation Next, a simulation study is presented to illustrate the performance of our results. Two features are studied: (1) the global performance of MLEs of parameters if the sample size increases; and (2) the effect of varying every shape parameter, γ , δ and λ , on the performance of the rest of estimators under consideration. Both issues are studied in the use of the flexible sinh-normal distribution as a regression model. Mathematics 2021,9, 1188 17 of 23 7.1. Simulation I The flexible sinh-normal distribution as a regression model with one covariate is considered, that is Yi=β0+β1xi+ei,i=1, . . . , n, where ei∼FSHN(γ , 0, 2, δ , λ) or, equivalently, Yi∼FSHN(γ , β0+β1xi , 2, δ , λ) independent for i=1, . . . , n. The vector of unknown parameters is (γ,β0,β1,δ,λ). As sample sizes, we consider n∈ { 25, 50, 75, 100, 500 } , m= 5000 simulations and values of the parameters (γ , β0 , β1 , δ , λ)=( 1, 1.75, 1.25, − 1.5, 2 ) . As x , covariate values of a U( 0, 1 ) distribution are considered. Values of the FSHN were generated from values of the FSN distribution, which were obtained by using the stochastic representation of FSN model given in [ 5 ]. The results are listed in Tables 5–7. As for measures of performance, the mean of estimates, their standard errors (sd) in parentheses, the relative bias and the square root of the mean squared error (MSE) are given. Table 5. Simulations for b γ. nbγ(sd) bias (bγ)√MSE 25 0.9968 (0.2878) 0.2026 0.3833 50 1.0674 (0.2338) 0.1461 0.2966 75 1.0972 (0.2221) 0.1222 0.2696 100 1.1235 (0.2196) 0.1012 0.2534 500 1.0960 (0.2114) 0.0432 0.2181 Table 6. Simulations for b β0,b β1. nb β0(sd) bias ( b β0)√MSE b β1(sd) bias ( b β1)√MSE 25 1.8205 (0.4539) 0.0403 0.4593 1.2866 (0.5022) 0.0293 0.5034 50 1.7290 (0.4054) 0.0120 0.4059 1.2681 (0.3424) 0.0145 0.3428 75 1.6813 (0.4027) 0.0393 0.4085 1.2621 (0.2702) 0.0097 0.2704 100 1.6348 (0.4074) 0.0659 0.4234 1.2617 (0.2373) 0.0094 0.2376 500 1.7141 (0.4475) 0.0205 0.5058 1.2514 (0.1027) 0.0011 0.1027 Table 7. Simulations for b δand bλ. nb δ(sd) bias (b δ)√MSE bλ(sd) bias (bλ)√MSE 25 −1.6190 (0.4530) 0.0793 0.4683 1.8495 (2.4948) 0.0753 2.4991 50 −1.6021 (0.3325) 0.0680 0.3478 2.0366 (1.0632) 0.0183 1.0638 75 −1.6048 (0.3013) 0.0699 0.3190 2.1029 (1.0358) 0.0515 1.0408 100 −1.6239 (0.2934) 0.0826 0.3185 2.1537 (1.0253) 0.0769 1.0366 500 −1.6025 (0.2981) 0.1083 0.3395 2.1683 (1.3550) 0.2342 1.4335 From the results in Tables 5–7, we can conclude that the MLEs are biased; in general, the standard error, relative bias and the squared root of MSE decrease if the sample size n increases. We highlight that, in this simulation, b β1 exhibits a good behavior, whereas the greater variability corresponds to bλ. Simulations were carried out by using the optim function of software R [ 25 ] by applying the Nelder–Mead method. 7.2. Simulation II Next, a sensitivity study about the effect of varying the shape parameters (and the sample size) on the estimates of other parameters is presented. Summaries varying the γ parameter are given in Table 8. Similarly, results varying the parameters δand λare given in Tables 9and 10, respectively. Mathematics 2021,9, 1188 18 of 23 Table 8. Empirical sd, relative bias and √MSE for the FSHN(γ, 1.75, 1.25, −1.5, 2)model. ˆγˆ β0ˆ β1ˆ δˆ λ γnsd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE 25 0.1223 0.0841 0.1241 0.1108 0.0214 0.1169 0.1600 0.0012 0.1600 0.7379 0.0040 0.7378 1.7371 0.0332 1.7382 50 0.0362 0.048 0.0381 0.0851 0.0124 0.0878 0.1122 0.0003 0.1122 0.3840 0.0005 0.3839 0.9076 0.0463 0.9123 75 0.0285 0.0338 0.0297 0.0732 0.0091 0.0749 0.0913 0.0019 0.0914 0.3288 0.0023 0.3288 0.8327 0.0698 0.8443 0.25 100 0.0246 0.0246 0.0254 0.0647 0.0078 0.0661 0.0766 0.0000 0.0766 0.3012 0.0093 0.3015 0.7395 0.0607 0.7494 200 0.0174 0.0126 0.0177 0.0508 0.0037 0.0512 0.0548 0.0005 0.0548 0.2465 0.0028 0.2465 0.4737 0.0365 0.4793 25 0.1801 0.0819 0.1903 0.3016 0.0561 0.3171 0.4195 0.0164 0.4199 0.5691 0.0069 0.5692 1.0871 0.0487 1.0914 50 0.1285 0.0437 0.1326 0.2413 0.0327 0.2479 0.2896 0.0012 0.2895 0.4314 0.0075 0.4315 1.1060 0.0696 1.1146 75 0.0989 0.0284 0.1012 0.2055 0.0208 0.2087 0.2326 0.0077 0.2327 0.3555 0.0071 0.3557 1.2154 0.0833 1.2266 0.75 100 0.0853 0.0221 0.0869 0.1829 0.0129 0.1843 0.1987 0.0023 0.1987 0.3147 0.0017 0.3147 0.8197 0.0753 0.8334 200 0.0596 0.0100 0.0600 0.1444 0.0062 0.1448 0.1401 0.0013 0.1401 0.2387 0.0007 0.2387 0.5718 0.0478 0.5797 25 0.3304 0.1037 0.3549 0.4336 0.1015 0.4685 0.5576 0.0212 0.5582 0.5984 0.0069 0.5985 3.4611 0.0179 3.4610 50 0.2370 0.0551 0.2468 0.3606 0.0534 0.3725 0.3872 0.0046 0.3872 0.4080 0.0031 0.4079 1.0082 0.0505 1.0131 75 0.1951 0.0388 0.2010 0.3106 0.0353 0.3166 0.3065 0.0065 0.3066 0.3286 0.0030 0.3286 0.9589 0.0750 0.9705 1.25 100 0.1703 0.0279 0.1738 0.2803 0.0255 0.2839 0.2596 0.0011 0.2596 0.2942 0.0010 0.2942 1.0168 0.0834 1.0303 200 0.1278 0.0138 0.1290 0.2219 0.0110 0.2228 0.1858 0.0027 0.1858 0.2198 0.0025 0.2198 0.6109 0.0472 0.6181 Mathematics 2021,9, 1188 19 of 23 Table 9. Empirical sd, relative bias and √MSE for the FSHN(1, 1.75, 1.25, δ, 2)model. ˆγˆ β0ˆ β1ˆ δˆ λ δn sd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE 25 0.7149 0.0835 0.7197 0.3610 0.0488 0.3709 0.5022 0.0080 0.5023 1.2311 0.1081 1.2358 5.1219 0.0603 5.1228 50 0.5221 0.0493 0.5244 0.2857 0.0213 0.2881 0.3528 0.0022 0.3528 0.7337 0.0750 0.7374 1.3591 0.0676 1.3657 75 0.1679 0.0347 0.1714 0.2515 0.0071 0.2518 0.2918 0.0073 0.2919 0.4105 0.0572 0.4144 0.9445 0.0657 0.9535 −1 100 0.1450 0.0270 0.1475 0.2311 0.0101 0.2317 0.2482 0.0003 0.2482 0.3644 0.0357 0.3661 0.8753 0.0660 0.8851 200 0.1009 0.0138 0.1018 0.1842 0.0041 0.1843 0.1683 0.0016 0.1683 0.2750 0.0186 0.2756 0.6144 0.0355 0.6185 25 0.9140 0.0934 0.9187 0.3478 0.0137 0.3486 0.4973 0.0149 0.4976 1.388 0.5287 1.4128 6.2775 0.0923 6.2795 50 0.3620 0.0535 0.3659 0.2714 0.0018 0.2714 0.3426 0.0003 0.3426 0.7485 0.3018 0.7635 1.3659 0.0951 1.3790 75 0.2473 0.0403 0.2506 0.2482 0.0012 0.2482 0.2754 0.0056 0.2755 0.5598 0.2215 0.5706 1.0731 0.0816 1.0853 −0.5 100 0.1735 0.0348 0.1770 0.2217 0.0041 0.2217 0.2415 0.0005 0.2415 0.4555 0.1987 0.4662 1.6011 0.0805 1.6090 200 0.1195 0.0216 0.1214 0.1850 0.0025 0.1850 0.1643 0.0005 0.1643 0.3406 0.1183 0.3456 0.7319 0.0442 0.7371 25 1.3881 0.0596 1.3892 0.3364 0.0093 0.3368 0.4892 0.0094 0.4893 2.0955 1.0817 2.1126 7.9980 0.2026 8.0075 50 0.6369 0.0448 0.6384 0.2713 0.0036 0.2713 0.3298 0.0109 0.3301 1.1466 0.6245 1.1570 2.6481 0.1318 2.6609 75 0.2728 0.0372 0.2753 0.2455 0.0026 0.2455 0.2654 0.0004 0.2653 0.6709 0.4761 0.6813 1.4155 0.1159 1.4342 −0.25 100 0.2067 0.0266 0.2084 0.2292 0.0074 0.2296 0.2272 0.0006 0.2272 0.5568 0.3704 0.5644 1.2191 0.1269 1.2451 200 0.1518 0.0164 0.1527 0.1963 0.0050 0.1965 0.1603 0.0002 0.1603 0.4422 0.2235 0.4456 0.9483 0.0807 0.9619 25 2.0897 0.1139 2.0925 0.3216 0.0211 0.3237 0.4466 0.0009 0.4466 3.4203 0.4546 3.4218 11.1958 0.479 11.2356 50 0.8816 0.0242 0.8818 0.2650 0.0007 0.2650 0.3036 0.0053 0.3037 1.6859 0.0111 1.6857 2.8473 0.2296 2.8838 75 0.5832 0.0087 0.5832 0.2383 0.0067 0.2386 0.2372 0.0038 0.2372 1.1711 0.1668 1.1717 2.1403 0.2225 2.1859 0.25 100 0.3810 0.0155 0.3813 0.2211 0.0100 0.2218 0.2083 0.0008 0.2083 0.9596 0.0649 0.9597 1.4984 0.2060 1.5538 200 0.1888 0.0036 0.1888 0.1847 0.0116 0.1858 0.1459 0.0002 0.1459 0.5933 0.1700 0.5947 0.9973 0.1412 1.0364 25 1.9152 0.1240 1.9191 0.2973 0.0191 0.2991 0.4197 0.0051 0.4197 3.4184 0.4688 3.4261 15.9834 0.5792 16.0237 50 1.0250 0.0667 1.0271 0.2435 0.0022 0.2435 0.2802 0.0050 0.2802 2.0445 0.2053 2.0469 3.6007 0.3098 3.6533 75 0.7143 0.0455 0.7157 0.2228 0.0076 0.2232 0.2284 0.0015 0.2284 1.4870 0.0934 1.4876 2.7262 0.2506 2.7716 0.5 100 0.3324 0.0137 0.3326 0.2076 0.0102 0.2084 0.1904 0.0029 0.1904 0.9070 0.0385 0.9071 1.3857 0.1926 1.4381 200 0.2089 0.0110 0.2091 0.1698 0.0124 0.1711 0.1359 0.0002 0.1359 0.6743 0.0530 0.6747 0.9694 0.1442 1.0113 25 1.9723 0.2052 1.9828 0.2573 0.0170 0.2590 0.3659 0.0037 0.3659 3.8193 0.5343 3.8562 18.1767 0.7324 18.2338 50 1.3344 0.1345 1.3410 0.2116 0.0025 0.2116 0.2548 0.0034 0.2548 2.8532 0.3136 2.8701 4.7951 0.3444 4.8439 75 0.8267 0.0764 0.8302 0.1906 0.0063 0.1909 0.2034 0.0013 0.2034 1.9203 0.1296 1.9245 2.4815 0.2837 2.5453 1.0 100 0.7316 0.0690 0.7347 0.1743 0.0106 0.1753 0.1711 0.0014 0.1711 1.7036 0.1055 1.7067 2.0484 0.2545 2.1105 200 0.4334 0.0374 0.4350 0.1350 0.0105 0.1362 0.1184 0.0000 0.1183 1.1408 0.0382 1.1413 1.1312 0.1508 1.1706 Mathematics 2021,9, 1188 20 of 23 Table 10. Empirical sd, relative bias and √MSE for the FSHN(1, 1.75, 1.25, −1.5, λ)model. ˆγˆ β0ˆ β1ˆ δˆ λ λn sd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE sd RB √MSE 25 0.2404 0.1088 0.2639 0.3883 0.0407 0.3948 0.5086 0.0004 0.5085 0.5320 0.1424 0.5732 0.5807 0.0151 0.5809 50 0.1763 0.0508 0.1835 0.3073 0.0110 0.3079 0.3521 0.0068 0.3522 0.3709 0.0838 0.3916 0.4488 0.0282 0.4496 75 0.1510 0.0325 0.1544 0.2648 0.0011 0.2648 0.2838 0.0019 0.2838 0.2999 0.0606 0.3134 0.3266 0.0193 0.3272 −1 100 0.1343 0.0284 0.1372 0.2301 0.0053 0.2303 0.2445 0.0029 0.2445 0.2604 0.0468 0.2696 0.2926 0.0096 0.2927 200 0.0979 0.0099 0.0983 0.1628 0.0013 0.1628 0.1693 0.0020 0.1693 0.1828 0.0194 0.1851 0.1893 0.0066 0.1894 25 0.2206 0.0935 0.2396 0.3742 0.0102 0.3746 0.5369 0.0040 0.5369 0.4774 0.1842 0.5515 0.3220 0.0099 0.3220 50 0.1729 0.0363 0.1766 0.2690 0.0003 0.2690 0.3778 0.0045 0.3778 0.3421 0.0819 0.3635 0.2121 0.0134 0.2122 75 0.1446 0.0258 0.1468 0.2160 0.0008 0.2160 0.2956 0.0060 0.2957 0.2785 0.0565 0.2911 0.1656 0.0009 0.1655 −0.5 100 0.1235 0.0212 0.1253 0.1817 0.0016 0.1817 0.2532 0.0014 0.2532 0.2390 0.0438 0.2479 0.1413 0.0004 0.1413 200 0.0870 0.0102 0.0876 0.1277 0.0008 0.1277 0.1754 0.0028 0.1754 0.1702 0.0224 0.1734 0.0984 0.0008 0.0984 25 0.2086 0.0831 0.2245 0.3642 0.0052 0.3643 0.5740 0.0008 0.5739 0.4812 0.1713 0.5454 0.2217 0.0212 0.2217 50 0.1527 0.0389 0.1576 0.2436 0.0003 0.2436 0.3869 0.0011 0.3869 0.3409 0.0785 0.3606 0.1442 0.0008 0.1442 75 0.1245 0.0247 0.1269 0.1979 0.0008 0.1979 0.3083 0.0039 0.3083 0.2841 0.0510 0.2942 0.1142 0.0040 0.1142 −0.25 100 0.1068 0.0194 0.1086 0.1714 0.0002 0.1714 0.2675 0.0016 0.2675 0.2410 0.0398 0.2482 0.0962 0.0064 0.0962 200 0.0774 0.0097 0.0780 0.1168 0.0008 0.1168 0.1813 0.0002 0.1813 0.1729 0.0183 0.1750 0.0688 0.0051 0.0688 25 0.2099 0.0859 0.2267 0.3646 0.0075 0.3648 0.5633 0.0033 0.5633 0.4799 0.1791 0.5499 0.2034 0.0335 0.2035 50 0.1537 0.0407 0.1590 0.2517 0.0003 0.2517 0.3856 0.0064 0.3857 0.3423 0.0830 0.3642 0.1473 0.0040 0.1473 75 0.1263 0.0288 0.1296 0.1979 0.0011 0.1979 0.3079 0.0018 0.3079 0.2813 0.0566 0.2938 0.1142 0.0006 0.1142 0.25 100 0.1094 0.0194 0.1111 0.1696 0.0032 0.1697 0.2644 0.0031 0.2644 0.2449 0.0399 0.2521 0.0984 0.0127 0.0984 200 0.0757 0.0092 0.0762 0.1185 0.0001 0.1185 0.1858 0.0011 0.1858 0.1689 0.0194 0.1714 0.0683 0.0008 0.0683 25 0.2258 0.0986 0.2463 0.3932 0.0166 0.3942 0.5549 0.0022 0.5549 0.4783 0.1898 0.5566 0.3089 0.0263 0.3091 50 0.1761 0.0447 0.1817 0.2753 0.0032 0.2754 0.373 0.0040 0.3730 0.3343 0.0893 0.3601 0.2124 0.0042 0.2124 75 0.1420 0.0317 0.1455 0.2131 0.0031 0.2131 0.2950 0.0021 0.2950 0.2779 0.0630 0.2935 0.1695 0.0016 0.1694 0.5 100 0.1236 0.0216 0.1255 0.1885 0.0037 0.1886 0.2571 0.0060 0.2572 0.2379 0.0434 0.2466 0.1416 0.0017 0.1416 200 0.088 0.0087 0.0884 0.1294 0.0014 0.1294 0.1765 0.0002 0.1765 0.1695 0.0175 0.1715 0.0969 0.0007 0.0969 25 0.2453 0.1065 0.2674 0.3931 0.0453 0.401 0.5165 0.0009 0.5165 0.5063 0.1383 0.5471 0.5795 0.0012 0.5795 50 0.1809 0.0613 0.1910 0.3007 0.0230 0.3034 0.3507 0.0021 0.3507 0.3503 0.0803 0.3704 0.4679 0.0108 0.4680 75 0.1550 0.0380 0.1595 0.2602 0.0117 0.2609 0.2878 0.0018 0.2877 0.2877 0.0541 0.2989 0.3588 0.0176 0.3592 1.0 100 0.1345 0.0308 0.1380 0.2310 0.0087 0.2315 0.2468 0.0000 0.2467 0.2534 0.0421 0.2611 0.2835 0.0060 0.2835 200 0.0982 0.0140 0.0992 0.1662 0.0037 0.1663 0.1725 0.0030 0.1726 0.1814 0.0223 0.1845 0.1915 0.0028 0.1915 Mathematics 2021,9, 1188 21 of 23 Comments to Table 8(varying γ>0). In this table, we consider β0=1.75, β1=1.25, δ=−1.5, λ=2, γ∈ {0.25, 0.75, 1.25} and sample sizes n∈ { 25, 50, 75, 100, 200 } . For c β0 and c β1 , it can be seen that the relative biases are small. Moreover, we can appreciate that the relative bias, standard error and √MSE decrease when the sample size increases, especially for c β1.b δand bλexhibit greater variability, but their relative bias, sd and √MSE also decrease when the sample size increases. As for b γ , we can point out that this estimator is also well behaved when n increases. Its relative bias for small sample sizes, n= 25, is small and the √MSE is a moderate value. Comments to Table 9(varying δ∈R). In this table, we consider γ= 1, β0= 1.75, β1= 1.25, λ= 2, δ∈ {− 1, − 0.5, − 0.25, 0.25, 0.5, 1 } and n∈ { 25, 50, 75, 100, 200 } . We can see that, for c β0 and c β1 , the relative bias is small even for n= 25, and the relative bias, sd and √MSE decrease when n increases. Again, b δ and bλ exhibit greater variability than the other estimators, but the error measures under consideration are satisfactory. The relative bias, sd and √MSE of b γ decrease when n increases, and their values are moderate even for small sample sizes. Comments to Table 10 (varying λ∈R). In this case, we consider γ= 1, β0= 1.75, β1= 1.25, δ=− 1.5, λ∈ {− 1, − 0.5, − 0.25, 0.25, 0.5, 1 } and n∈ { 25, 50, 75, 100, 200 } . Again, a similar behavior to the previously explained is obtained. As final conclusion, we point out that these simulation studies cover a variety of situations as for the shapes of the FSHN distribution (unimodal and bimodal) and suggest that the estimators of the parameters are consistent when the sample size increases [26]. 8. Discussion The BS distribution is an asymmetric model used for survival time data and material lifetime subject to stress as it can be seen in [ 27 ]. Due to its practical and theoretical interest, a number of generalizations can be found in literature. We can cite the extensions provided in [ 28 ] to the family of elliptical distributions; in [ 10 ] based on the elliptical asymmetric distributions; and the extended Birnbaum–Saunders (EBS) distribution introduced in [7]. Other generalizations intend to solve specific deficiencies observed when this model is fitted to a dataset. In this sense, we can cite the epsilon-Birnbaum–Saunders model introduced in [ 8 ] to accommodate outliers; the extension based on the slash-elliptical family of distributions given in [ 3 ]; and the generalized modified slash Birnbaum–Saunders, which is based on the work in [ 29 ] and proposed in [ 30 ]. All these extensions are appropriate to fit data with greater or smaller asymmetry (or kurtosis) than that of the usual BS model, but they are not appropriate for fitting bimodal data. This issue is of great interest, as discussed by Olmos et al. [ 9 ], Bolfarine et al. [ 31 ], Martínez-Flórez et al. [ 2 ] and Elal-Olivero et al. [ 32 ]. The flexible BS distribution was proposed to model skewness and to fit data with and without bimodality. In this paper, this model is spread out to be used as a regression model in those situations in which regression models based on other generalizations of BS, such as those proposed in [7,13,33–37], do not provide a satisfactory fit. Author Contributions: Conceptualization, G.M.-F., I.B.-C. and H.W.G.; formal analysis, G.M.-F., I.B.-C. and H.W.G.; investigation, G.M.-F., I.B.-C. and H.W.G.; methodology, G.M.-F., I.B.-C. and H.W.G.; Software, G.M.-F. and I.B.-C.; supervision, H.W.G.; validation, G.M.-F. and I.B.-C.; and visualization, H.W.G. All authors have read and agreed to the published version of the manuscript. Funding: The research of G. Martínez-Flórez was supported by project: Distribuições de Probabilidade Mutivariadas Assimétricas e Flexíveis, Universidade Federal do Ceará, Fortaleza, Brazil. The research of H. W. Gómez was supported by Grant SEMILLERO UA-2021, Chile. Data Availability Statement: Details about data availability are given in Sections 6.1–6.3. Mathematics 2021,9, 1188 22 of 23 Acknowledgments: G. Martínez-Flórez acknowledges the support given by Universidad de Córdoba, Montería, Colombia, and Universidade Federal do Ceará, Fortaleza, Brazil. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: BS Birnbaum–Saunders cdf cumulative distribution function FBS Flexible Birnbaum-Saunders FSHN Flexible sinh-normal FSN Flexible skew Normal pdf probability density function SHN Sinh-normal References 1. Rieck, J.R.; Nedelman, J.R. A log-linear model for the Birnbaum-Saunders distribution. Technometrics 1991,33, 51–60. 2. Martínez-Flórez, G.; Barranco-Chamorro, I.; Bolfarine, H.; Gómez, H.W. Flexible Birnbaum-Saunders Distribution. Symmetry 2019,11, 1305. [CrossRef] 3. Gómez, H.W.; Olivares-Pacheco, J.F.; Bolfarine, H. An extension of the generalized BirnbaumSaunders distribution. Stat. Probab. Lett. 2009,79, 331–338. [CrossRef] 4. Azzalini A. A class of distributions which includes the normal ones. Scand. J. Stat. 1985,12, 171–178. 5. Gómez, H.W.; Elal-Olivero, D.; Salinas, H.S.; Bolfarine, H. Bimodal Extension Based on the Skew-normal Distribution with Application to Pollen Data. Environmetrics 2011,22, 50–62. [CrossRef] 6. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. Censored bimodal symmetric-asymmetric families. Stat. Interface 2018 ,11, 237–249. [CrossRef] 7. Leiva, V.; Vilca-Labra, F.; Balakrishnan, N.; Sanhueza, A. A skewed sinh-normal distribution and its properties and application to air pollution. Commun. Stat. Theory Methods 2010,39, 426–443. [CrossRef] 8. Castillo, N.; Gomez, H.W.; Bolfarine, H. Epsilon Birnbaum-Saunders distribution family: Properties and inference. Stat. Pap. 2011,52, 871–883. [CrossRef] 9. Olmos, N.M.; Martínez-Flórez, G.; Bolfarine, H. Bimodal Birnbaum-Saunders Distribution with Application to Corrosion Data. Commun. Stat. Theory Methods 2017,46, 6240–6257. [CrossRef] 10. Vilca-Labra, F.; Leiva-Sanchez, V. A new fatigue life model based on the family of skew-elliptical distributions. Commun. Stat. Theory Methods 2006,35, 229–244. [CrossRef] 11. Calderín Ojeda, E.; Gómez Déniz, E.; Barranco Chamorro, I. Modelling Zero-Inflated Count Data With a Special Case of the Generalised Poisson Distribution. Astin Bull. 2019,49, 689–707. [CrossRef] 12. Rivera, P.A.; Barranco-Chamorro, I.; Gallardo, D.I.; Gómez, H.W. Scale mixture of Rayleigh distribution. Mathematics 2020 ,8, 1842. [CrossRef] 13. Lemonte, A. A log-Birnbaum-Saunders regression model with asymmetric errors. Improved statistical inference for the twoparameter Birnbaum-Saunders distribution. J. Stat. Comput. Stat. Simul. 2012,82, 1775–1787. [CrossRef] 14. Chambers, J.M.; Cleveland, W.S.; Kleiner, B.; Tukey, P.A. Graphical Methods for Data Analysis; Wadsworth: Belmont, CA, USA, 1983. 15. Gokhale, S.; Khare, M. Statistical behavior of carbon monoxide from vehicular exhausts in urban environments. Environ. Model. Softw. 2007,22, 526–535. [CrossRef] 16. Nadarajah, S. A truncated inverted beta distribution with application to air pollution data. Stoch. Environ. Res. Risk Assess. 2008 , 22, 285–289. [CrossRef] 17. Akaike, H. A new look at statistical model identification. IEEE Trans. Autom. Control 1974,19, 716–723. [CrossRef] 18. Kullback, S.; Leibler, R.A. On information and sufficiency. Ann. Math. Stat. 1951,22, 79–86. [CrossRef] 19. Vuong, Q. Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica 1989,57, 307–333. [CrossRef] 20. Arnold, B.C.; Gómez, H.W.; Salinas, H.S. On multiple constraint skewed models. Statistics 2009,43, 279–293. [CrossRef] 21. Elal-Olivero, D. Alpha-skew-normal distribution. Proyecciones J. Math. 2010,29, 224–240. [CrossRef] 22. Kim, H.J. On a class of two-piece skew-normal distribution. Statistic 2005,39, 537–553. [CrossRef] 23. Marin, J.-M.; Mengersen, K.; Robert, C. Bayesian modelling and inference on mixtures of distributions. In Handbook of Statistics; Rao, C., Dey, D., Eds.; Springer: New York, NY, USA, 2005; Volume 25, pp. 1–56. 24. Chan, P.S.; Ng, H.K.T.; Balakrishnan, N.; Zhou, Q. Point and interval estimation for extreme-value regression model under Type-II censoring. Comput. Stat. Data Anal. 2008,52, 4040–4058. [CrossRef] 25. R Development Core Team. A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021. Mathematics 2021,9, 1188 23 of 23 26. Barranco Chamorro, I.; Jimenez Gamero, M.D. Asymptotic Results in Partially Non-Regular Log-Exponential Distributions. J. Stat. Comput. Simul. 2012,82, 445–461. [CrossRef] 27. Birnbaum, Z.W.; Saunders, S.C. A New Family of Life Distributions. J. Appl. Probab. 1969,6, 319–327. [CrossRef] 28. Díaz-García, J.A.; Leiva-Sánchez, V. A new family of life distributions based on the elliptically contoured distributions. J. Stat. Plann. Inference 2005,128, 445–457. [CrossRef] 29. Reyes, J.; Barranco-Chamorro, I.; Gómez, H.W. Generalized modified slash distribution with applications. Commun. Stat. Theory Methods 2020,49, 2025–2048. [CrossRef] 30. Reyes, J.; Barranco-Chamorro, I.; Gallardo, D.I.; Gómez, H.W. Generalized Modified Slash Birnbaum-Saunders Distribution. Symmetry 2018,10, 724. [CrossRef] 31. Bolfarine, H.; Martínez–Flórez, G.; Salinas, H.S. Bimodal symmetric-asymmetric families. Commun. Stat. Theory Methods 2018 ,47, 259–276. [CrossRef] 32. Elal-Olivero, D.; Olivares-Pacheco, J.F.; Venegas, O.; Bolfarine, H.; Gómez, H.W. On Properties of the Bimodal Skew-Normal Distribution and an Application. Mathematics 2020,8, 703. [CrossRef] 33. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. An alpha-power extension for the Birnbaum-Saunders distribution. Statistics 2014,48, 896–912. [CrossRef] 34. Santana, L.; Vilca, F.; Leiva, V. Influence analysis in skew-Birnbaum-Saunders regression models and applications. J. Appl. Stat. 2011,38, 1633–1649. [CrossRef] 35. 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, 316–332. [CrossRef] [PubMed] 36. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. The Log-Linear Birnbaum-Saunders Power Model. Methodol. Comput. Appl. Probab. 2017,19, 913–933. [CrossRef] 37. Ortega, E.M.; Bolfarine, H.; Paula, G.A. Influence diagnostics in generalized log-gamma regression models. Comput. Stat. Data Anal. 2003,42, 165–186. [CrossRef]