Extended Generalized Sinh-Normal Distribution

: Positively skewed data sets are common in different areas, and data sets such as material fatigue, reaction time, neuronal reaction time, agricultural engineering, and spatial data, among others, need to be ﬁtted according to their features and maintain a good quality of ﬁt. Skewness and bimodality are two of the features that data sets like this could present simultaneously. So, ﬂexible statistical models should be proposed in this sense. In this paper, a general extended class of the sinh-normal distribution is presented. Additionally, the asymmetric distribution family is extended, and as a natural extension of this model, the extended Birnbaum–Saunders distribution is studied as well. The proposed model presents a better goodness of ﬁt compared to the other studied models.


Introduction
When materials are exposed to pressure or stress levels, material structural damage could occur. This is known as material fatigue, and a statistical model to fit random variables to model the failure time of fatigue for material was proposed by Birnbaum and Saunders (1969) [1], known in the literature as Birnbaum-Saunders distribution and generally denoted by BS(α, β), where α > 0 represents a shape parameter and β > 0 is a scale parameter and the median of the distribution. Later, Desmond (1985) [2] showed that the BS distribution describes the time failure that occurs when some kind of damage is accumulated after a given time.
A distribution associated with that of Birnbaum and Saunders (BS) is the sinh-normal (SHN) distribution. This distribution, introduced by Rieck and Nederman (1991) [3], is based on a nonlinear transformation of a normal distribution. Therefore, let Z = 2 α sinh Y−ξ σ ∼ N(0, 1), where ξ and σ are location and scale parameters, respectively, and α is a shape parameter. Then, random variable (r.v.) Y follows a sinh-normal distribution denoted by SHN(α, ξ, σ). A probability density function (pdf) of a random variable with a SHN distribution is given by where b y = b(y, α, ξ, σ) = 2 α sinh y−ξ σ , b y = 2 ασ cosh y−ξ σ is a derivative of b y with respect to y, and φ(.) is the pdf of the normal distribution. It can be shown that the s-th moment of a random variable with a SHN distribution is given by where From this result, we obtain that µ 1 = E(Y) = ξ, µ 2 = ξ 2 + σ 2 c 2 (α) and var(Y) = σ 2 c 2 (α). It can also be shown, from the two central moments, that E(Y − ξ) 3 = 0 is a symmetric distribution with respect to ξ. In general, we find that E(Y − ξ) s = σ s c s (α). (4) Asymmetric extensions of the SHN model have been considered, for instance, by Leiva [7], who presented the Proportional Hazard BS model (PHBS).
An important feature of this distribution is its robustness concerning the estimation of its parameters-an aspect that was analyzed by Barros, Paula, and Leiva (2010) [8]. Moreover, extensions of this distribution to an elliptical family and to a skew elliptical family (this latter is known as the double generalized BS distribution) have been studied by Díaz-García et al. (2005) [9] and by Vilca-Labra and Leiva-Sánchez (2006) [10], respectively. Martínez-Flórez et al. [11] present an extension to a power-skew-elliptical family. Other types of extensions have also been considered by authors such as Castillo et al. (2011) [12], Cordeiro and Lemonte. (2014) [13], and Reyes et al. (2018) [14].
All these extensions are particularly characterized by fitting skew unimodal data while not being appropriate to fit bimodal data. However, Martinez-Florez et al. (2017) [6] and Olmos et al. (2017) [15] recently presented BS models to fit positive bimodal data. Likewise, Bolfarine et al. (2011) [16] introduced another model to fit positive bimodal data generated by the log-skew-normal distribution. Cortés et al. (2018) [17] presented a class of extended distribution, generated by the pdf, g(x), of the bimodal-normal distribution. They specifically defined a general class of distributions with a pdf given by where g(x) is a pdf, ≥ 0 is a shape parameter, and h is a continuous positive function such that κ = E g (h(X)) < ∞. This distribution is called the "general class of distributions". Further, the authors study the normal, t-student, Laplace, and BS distributions as special cases. As special cases of this family, from Elal-Olivero et al. (2010) [18] and the bimodal logskew-normal of Bolfarine et al. (2011) [16], the class of bimodal skew-elliptical distributions can be found in the literature.
In this study, we analyze and study a general extended class of the SHN distributions, which has an extra parameter to the SHN model and introduces flexibility to the SHN distribution. Additionally, an extension of the BS model is presented. We highlight that this kind of model could be applied to data sets related to material fatigue [2], reaction and neuronal reaction time [19], agricultural engineering [20], and spatial data (see [21,22], among others). Note that this distribution can fit bimodal data sets, which could be present in problems where the population is divided by groups such as gender or different levels of HIV-RNA. This paper unfolds as follows: Section 2 presents the extended sinh-normal distribution (ESHN) model. Section 3 provides the statistical approximation of the moments of an ESHN random variable. Section 4 outlines the properties of the extended generalized sinh-normal (EGSHN) distribution. Section 5 shows the features of the extended sinhnormal regression model. Section 6 presents the results of a simulation study to analyze the properties of the EGSHN model. Section 7 develops two numerical illustrations to evaluate the relevance of the EGSHN. Section 8 discusses the statistical and practical implications of the proposed distribution.

Extended Sinh-Normal Distribution
We now propose an extension of the SHN case introduced by Rieck and Nederman (1991) [3] to a general class of distributions. Then, taking h(y) = b 2 y and g(y) = b y φ(b y ) in equation (7), where b y and b y are defined as in (1), the extended class of SHN distribution is defined through the pdf given by with ξ ∈ R and α, σ > 0 defined as the SHN distribution of Rieck and Nederman (1991) [3], and γ ∈ R+ as a shape parameter. This distribution is denoted by ESHN(α, ξ, σ, γ), which, for some values of the γ parameter, could be bimodal. In this context, γ could be considered as a bimodality parameter. It can be easily deduced that for γ = 0, the SHN distribution is obtained, and if γ → ∞, then f (y) → b 2 y b y φ(b y ), which is a new family of distributions. It can be also shown that a random variable where EN is the extended normal distribution, as studied by Cortés et al. (2018) [17]. Rieck et al. [3]). Next, it is possible to conclude that when α → 0, then f (y) ∼ EN. Figure 1 shows the pdf for an extended sinh-normal (ESHN) distribution with two parameters, respectively, where it can be seen that the distribution, for γ = 0, is symmetric and unimodal in cases (b) and (c), and for the other cases, it is strongly bimodal.

Extended Generalized Sinh-Normal Distribution
The ESHN distribution studied in the previous section has the main feature of bimodal symmetric data fitting, for which it is necessary to extend this model to asymmetric data; this asymmetric extension comes from the results found by Azzalini and Capitanio (2003) [23] and Azzalini (2005) [24], who showed that if g is a pdf symmetric around zero, and H is a cdf so that its density h is symmetric around zero as well, then for any odd function w(x), we find that Then, since the pdf ( f (y)) of the standard distribution, ESHN(α, 0, 1, γ), is continuous and symmetric around zero; H(·) = Φ(·) is an absolutely continuous distribution function that is symmetric around zero and whose density, φ(·), is also symmetric around zero; and w(y) = b y = b y (α, 0, 1, γ) = γ 2 α sinh(y), for constant values α, γ, is an odd function, then 2 f (y)Φ(λb y ) is a pdf for any λ ∈ R. Thus, the location-scale-extended generalized sinh-normal (EGSN) distribution is defined through the pdf given by where λ is a skewness parameter. So, in this work, this model is denoted by EGSHN (α, ξ, σ, γ, λ). Figure 2 shows the pdf of the EGSHN distribution for different values of the parameters. As observed, the distribution could be unimodal or bimodal depending on such values.
Likewise, for γ = 0, we obtain the asymmetric SHN distribution based on the model considered by Leiva et al. (2010) [4]. In addition, if γ → ∞, analogously to the ESHN case, we have the new family of distributions g(y) → 2b 2 y b y φ(b y )Φ(λb y ). Now, for γ = λ = 0, an SHN model is followed. It can also be shown that a random variable where ESN is the extended skew-normal distribution, which was introduced by Elal-Olivero et al. (2009) [25]. As in the ESHN model, it can be shown that when α → 0, then g(y) ∼ ESN.
For p, as in the case of the ESHN model, the pdf of the EGSHN(α, ξ, σ, γ, λ) distribution could be written as From (18), it follows that the cdf of EGSHN(α, ξ, σ, γ, λ) is given by where Φ(·; λ) and φ(·; λ) denote the cdf and pdf of the skew-normal distribution, respectively, with a location parameter of 0 and a scale parameter of 1. Then, the survival, risk (or hazard), and inverse risk functions of this distribution are, respectively, given by where S SSHN , h SSHN (t), and R SSHN (t) are the survival, hazard, and inverse risk functions of the skew-SHN function, SSHN(α, ξ, σ, λ), respectively.

Stochastic Representation
The stochastic representation of the EGSHN model is based on Elal-Olivero (2010) [18] and Elal-Olivero, et al. (2009) [25]. This is presented below. Definition 1. If the random variable X has a pdf given by then we say that X follows a bimodal-normal distribution, and it is denoted as X ∼ NB (see [18]).

Remark 1. Let W and U independent random variables with W ∼ χ 2
(3) , a chi-square distribution with three degrees of freedom, and U a pdf, such that P(U = −1) Remark 2. Let Y ∼ NB and considers the random variable X 2 which is defined as Then, , and say that U ∼ U(0, 1) is a uniform random variable, independent of X 1 and X 2 . If Proof.
Proof. Note that W(X) = λ 2 α sinh(X) is an odd function and the density of both the random variables −X and the distribution function of Z are symmetric around zero; thus, applying Lemma 1 as in [24], the result follows.

Location-Scale Extension
Let W ∼ EGSHN(α, γ, λ) and let V = µ + σW, then Then, the pdf of V is given by

Moments of an Extended Generalized Sinh-Normal Random Variable
For a random variable Y ∼ EGSHN(α, ξ, σ, γ, λ), the r-th moment is given by where
As in the EGSHN model, it can be shown that the cdf of the EGBS(α, β, γ, λ) distribution is given by Then, the survival, risk, or hazard, and inverse risk functions of this distribution are given, respectively, by and The moments of an EGBS random variable with parameters α, β, γ, and λ can be obtained by means of the following expression: where E BSSN (·) is the expectation operator of the Birnbaum-Saunders skew-normal distribution (BSSN). To calculate the mean, variance, and the skewness and kurtosis coefficients of the EGBS model, the expressions of the corollaries 2.1 and 2.2 of Vilca et al. [10] could be used. For (α, β, γ, λ), the log-likelihood function corresponding to the random sample t 1 , t 2 , . . . , t n is The score function is composed of the following elements: Φ(λa t i ) , with i = 1, . . . , n. Iterative numerical methods must be used to solve this system of nonlinear equations.

Extended Sinh-Normal Regression Model
One of our main goals is to develop a log-BS regression model based on the ESHN(α, ξ,  σ, γ, λ) model. This regression model will be an optimal alternative to the log-BS model introduced by Rieck et al. (1991) [3] in order to fit bimodal or survival asymmetric data. Now, the extended generalized log-linear BS regression model is defined following the considerations given by Rieck et al. (1991) [3] and considering that Y i is a dependent variable; that a set (p) of explanatory variables, denoted by x i = (x i1 , x i2 , · · · , x ip ) , is given; and that θ = (θ 1 , θ 2 , · · · , θ p ) is a p-dimensional vector of unknown parameters, where a linear predictor, ξ i = x i θ, is obtained for i = 1, 2, . . . , n.
Then, let us suppose that T 1 , T 2 , ...T n are independent and identically distributed random variables (i.i.d. r.v.) such that T i ∼ EGBS(α i , β i , γ i , λ i ). Now, let us suppose that the distribution of T i is independent of the set of explanatory variables, .., n, with θ = (θ 1 , θ 2 , ..., θ p ), being a p-dimensional vector of unknown parameters.

2.
The shape, bimodality, and skew parameters do not consider x i ; i.e., α i = α, γ i = γ, and λ i = λ for i = 1, 2, ..., n. Let us suppose that Y i = log(T i ). Then, the extended generalized log-linear BS regression model is defined as where i ∼ EGSHN(α, 0, 2, γ, λ), for i = 1, . . . , n and y i is the log-survival for the i-th individual. This model is denoted by MRESHN(α, θ, 2, γ, λ). When γ = λ = 0 it follows the log-BS regression model, LBS(α, θ, 2), of Rieck et al. [3]; i.e., the MRESHN model is more flexible than the log-BS model in terms of skewness and kurtosis. When λ = 0, it follows that i ∼ ESHN(α, 0, 2, γ), for i = 1, . . . , n; then, important results are obtained. Thus, we find that E( i ) = 0 and Var( i ) = 4c 20 (α); additionally, as the errors are independent random variables, then for i = i , it follows that Cov( i , i ) = 0. Furthermore, considering that the explanatory variables are independent of the shape parameter, from the above results, it is therefore possible to conclude that Y i ∼ ESHN(α, x i θ, 2, γ) for i = 1, . . . , n; furthermore, like = 0, it can be shown that E(Y i ) = ξ i = x i θ, so the linear estimators for θ can be derived from the ordinary least squares method, whose solution is given byθ So, a biased estimator for c 2 (α) could bê (26) is a linear regression model similar to the models of the theory of linear models, with the characteristic that the random component follows a EGSHN(α, 0, 2, γ, λ) distribution; that is, assume that these parameters are located around zero and have a scale parameter equal to 2. Then, the interpretation of their parameters, with relation to the observed variable Y, is given in the same manner as the linear regression model.
The least squares estimator (θ * ) may be used to initialize the iterative numerical process for θ, and with these initial values, we can calculateα The elements of the observed information matrix, defined as minus the second derivative of the log-likelihood function, are denoted by I αα , I αθ , I αγ , I αλ I θ j θ k , I θ j γ , I θ j λ , I γγ , and I λλ , and are, respectively, given by (1 + γξ 2 i2 ) 2 , x ij x ik w i −ξ i2 + λξ 2 i1 (λξ i2 + w i ) , The information matrix (Ψ) could be estimated as the expected value of the elements of the observed information matrix, which must be calculated using numerical approximation methods. For λ = 0, the sub-matrix of vector (α, θ, γ) of the information matrix has the following elements: For γ > 0, the determinant of the information sub-matrix is not equal to zero; that is, |Ψ(α, θ, γ)| = 0. Thus, it is possible to conclude that the information matrix of the ESHN model is nonsingular. Likewise, for λ = 0, we find that x ij , ψ λγ = 0, and ψ λλ = n 2 where E(·) denotes the expected value function, and G(a, b) is the gamma distribution with parameters a and b, respectively. The rows or columns of the information matrix of the parameters vector (α, θ , γ, λ) of the regression model MRESHN in the case of λ = 0 are linearly independent, Therefore, the information matrix of this model is nonsingular, and its inverse is the variance-covariance matrix, Σ α,θ,γ,λ , of the estimator vector of maximum likelihood of the parameters vector; then, the estimated standard errors of the estimators are the square root of the matrixΣ α,θ,γ,λ .

Simulation Study
A simulation study is presented to analyze the performance of the maximum likelihood estimation of the parameters of the EGSHN(α, β 0 , β 1 , γ, λ) regression model. In general, m = 5000 simulations with n = 30, 60, 90, 120, and 500 were generated for different scenarios. So, the following model, with β 0 = 2.25 and β 1 = 0.75, is studied.
Note that, for i = 1, 2, · · · , n, x i takes values of a uniform random variable (0, 1), X ∼ U(0, 1), and ∼ EGSHN(α, 0, 2, γ, λ). The statistics of the empirical standard deviation (sd), relative bias (RB) and √ MSE for the EGSHN model are calculated. The results for each studied scenario are described as follows. Scenario 1 (Varying α): In this scenario, the used alpha values are 0.75, 1.75, and 2.75, keeping γ = 1.5 and λ = 1 fixed. For each studied case of the EGSHN(α, 2.25, 0.75, 1.5, 1) model (see Figure 3 and Table 1), it can be seen that for the parameters β 0 and β 1 , the relative biases are small, especially for the parameter β 1 . Note that the statistics RB, sd, and √ MSE decrease when the sample size is increased. Additionally, it can also be observed that the relative bias of the parameter α is not very important (not very large). Something similar is observed for the parameter λ. On the other hand, the relative bias and the root of the mean square error of the parameter γ are not very small, especially when the sample size is small (30 or 60).
The statistics RB, sd, and √ MSE for the parameters α, γ, λ decrease when the sample size increases. These results guarantee a lack of bias and the asymptotic consistency of the estimates of the parameters α, β 0 , β 1 , γ, λ.
The parameters α and λ show very small relative biases, but the parameter γ presents large relative biases and √ MSE. This behavior of the gamma parameter is more striking for small sample sizes.
The asymptotic consistency of the estimates of the parameters α, β 0 , β 1 , γ, λ can be guaranteed because the calculated statistics for α, γ, λ decrease when sample size increases.
Scenario 3 (Varying λ): In scenario 3, where the EGSHN(1.75, 2.25, 0.75, 2, λ) model is considered, the used lambda values are 0.5, 1.5, and 3, keeping α = 1.75 and γ = 2 fixed (see Figure 5 and Table 3). As in the previous simulation scenarios, the relative biases for the parameters β 0 and β 1 are small, especially for the parameter β 1 . It is possible to see that the statistics RB, sd, and √ MSE tend to decrease when the sample size increase. The RB and √ MSE of the parameter γ are large for small sample sizes. The parameter α presents small relative biases, but the parameter λ shows a large sd and √ MSE when λ = 3. In general, we can see that the asymptotic consistency of the estimates of the parameters α, β 0 , β 1 , γ, λ is satisfied.   In this paper, we employed the dataset studied by Hirose (1993) [28] to show the relevance of the EGSHN model. This dataset contains the results of an accelerated life test on polyethylene terephthalate (PET) (used in electrical isolation) in SF6 gas-insulated transformers. Such an accelerated life test was performed at four voltage levels (5, 7, 10, and 15), with 7, 15, 10, and 9 observations per level, respectively. The main purpose of the study was to analyze the resistance times (t) of the insulating films at different voltages (v). Hence, we here consider the following regression model: In this work, the SHN and ESHN models are fitted (see Ortega et al. [29]). To compare these models, the Akaike (AIC) information criterion (Akaike, [30]) and the corrected Akaike (AICC) criterion (Cavanaugh,[31]) are implemented. These measures are defined as follows: where p is the number of parameters andˆ (·) is the log-likelihood function evaluated for the MLEs of parameters. The best model is that with the smallest AIC or AICC value. To fit the bivariate model, we used the optim function of the statistical package R Core Team. The estimated parameters of these models, accompanied by their standard errors in parentheses, are obtained using the maximum likelihood method. Table 4, shows the results. Note that according to the AIC and AICC, the EGSHN and ESHN models present the best fits. In order to identify atypical observations and/or the misspecifications of models, we analyzed the transformation of the martingale residual, rMTi , proposed by Barros et al. [8]. These residuals are defined by where rM i = α i + log(S(e i ,θ)) is the martingale residual proposed by Ortega et al. [32], where α i = 0, 1 indicates whether the i-th observation is censored or not, respectively; sgn(rM i ) denotes the sign of rM i ; and S(e i ;θ) represents the survival function evaluated in e i , whereθ are the MLE for θ.
The plots of rMT i with generated confidence envelopes is presented in Figure 6. From this figure, we can see clearly that the EGSHN model fits better to the data than the SHN and ESHN models, since, in these cases, there are no observations that lie outside the envelopes. In Figure 6b, we observe two points: one in the border and the other outside it. There is also a point far from the set of observations. Since these points could be values influencing the estimates of the parameters, we calculated the generalized Cook's distance (GCD) and showed the residual components of the deviation plot. Figure 7 illustrates the behavior of these two statistics. As can be seen in this figure, observations 8, 23, and 34 could be possible influential values. Thus, to calculate their impact on the estimates, we computed the estimates by eliminating each of these observations or groups of them.
The relative change (RC), in percentage, of each parameter estimate is used to evaluate the effect of the potentially influential case. The RC is given by RC(θ (i) ) = 100 × |(θ − θ (i) )/θ|, whereθ (i) denotes the MLE of θ after removing the i-th observation. Table 5 lists the obtained RC values. According to this table, the relative changes of the MLE of parameter λ are excessively pronounced in all the models, mainly for observation 8

Illustration 2
We used the dataset reported by the Center for Applied Statistics of the Institute of Mathematics and Statistics at the University of São Paulo to illustrate the relevance of the EGBSSN model. This dataset contains data on the amount of DNA within the nucleus (ploidy) of mammary cells (250 samples) from women with breast cancer. The ploidy variable exhibits a bimodal asymmetric performance, with D = 0.0399 and p-value = 0.0059 after performing the Hartigan and Hartigan (1985) [33] bimodality test. In addition, according to its descriptive statistics (presented in Table 6), this variable has a considerable positive skewness and high kurtosis. Additionally, Figure 8a shows a bimodal distribution of the data, which is why the EGBSSN model becomes an alternative to fit this kind of data. In this study, the Birnbaum-Saunders (BS) and Birnbaum-Saunders skew-normal models were fitted. We also considered the following bimodal models: the log-skew-normal model developed by Bolfarine et al. (2011) [16], the extended class (EBN) introduced by Cortez et al. (2018) [17], and the EGBS model. Table 7 presents the maximum likelihood estimates and AIC and AICC values of five models. As observed in this table, the BS and BSSN models provide a poor fit to the DNA dataset. Conversely, the EGBS model shows the best fit among the fitted models, which is explained by its flexibility to fit asymmetric bimodal data. Figures 8a,b shows the estimated densities of the models with the best fit (BSSN, BLSN, EBS, and EGBS); the empirical cumulative distribution functions of the BLSN, EBS, and EGBS models; and the parameter estimates. Note that the EGBS model provides the best fit when compared to the BSSN, BLSN, EBN, and BS models.

Discussion
Bimodality and skewness are two common features that may be present in data from engineering, geo-spatial, medicine, and other areas. The natural complexity of data from these areas needs to be fitted using models that offer great flexibility and goodness of fit. In some cases, the data only present positive skewness, or in other cases, the distribution is only bimodal, but these two features could be present simultaneously. Thus, in this paper, a distribution capable of fitting bimodal and positively skewed data sets was proposed. In addition, the extended Birnbaum-Saunders distribution was studied as well.
Although there are proposals in the literature such as those cited in the previous sections, which allow the fitting of asymmetric or bimodal data, our model has the characteristic of modeling data that simultaneously have these two characteristics and has great flexibility and goodness of fit for data with these conditions. It is a new, promising, and user-friendly option to consider in statistical analysis.
To conclude, the EGBS model was proposed as a new statistical distribution suitable to fit real data sets with positive skewness and bimodality that presents a great performance compared with models available in the literature and that aim to achieve the same objective. Funding: This research was funded by the Instituto Tecnológico Metropolitano (ITM) through the project (P20244) Priorización de zonas para la restauración ecológica y de uso público mediante la armonización de técnicas de mapeo participativo y modelación espacial multicriterio en el municipio de Belmira, Antioquia, Colombia.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.