Next Article in Journal
A Note on Parametric Kinds of the Degenerate Poly-Bernoulli and Poly-Genocchi Polynomials
Next Article in Special Issue
An Asymptotic Test for Bimodality Using The Kullback–Leibler Divergence
Previous Article in Journal
Numerical Study of Bubble Behavior under Gradient Flows during Subcooled Flow Boiling in Vertical Flow Channel
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Likelihood-Based Inference for the Asymmetric Beta-Skew Alpha-Power Distribution

by
Guillermo Martínez-Flórez
1,
Roger Tovar-Falón
1,* and
Marvin Jimémez-Narváez
2
1
Departamento de Matemáticas y Estadística, Facultad de Ciencias Básicas, Universidad de Córdoba, Montería 230027, Colombia
2
Departamento de Ingeniería Industrial, Facultad de Ciencias e Ingenierías, Universidad del Sinú, Montería 230027, Colombia
*
Author to whom correspondence should be addressed.
Symmetry 2020, 12(4), 613; https://doi.org/10.3390/sym12040613
Submission received: 18 March 2020 / Revised: 25 March 2020 / Accepted: 26 March 2020 / Published: 13 April 2020
(This article belongs to the Special Issue Symmetric and Asymmetric Bimodal Distributions with Applications)

Abstract

:
This paper introduces a new family of asymmetric distributions that allows to fit unimodal as well as bimodal and trimodal data sets. The model extends the normal model by introducing two parameters that control the shape and the asymmetry of the distribution. Basic properties of this new distribution are studied in detail. The problem of estimating parameters is addressed by considering the maximum likelihood method and Fisher information matrix is derived. A small Monte Carlo simulation study is conducted to examine the performance of the obtained estimators. Finally, two data set are considered to illustrate the developed methodology.

1. Introduction

Probability distributions for modeling data with high asymmetry and bimodality have been proposed by several authors. They stand out those distributions that have been supported under the general structure of the skew-symmetric probability density function (pdf) given by
f S S ( z ; λ ) = 2 f ( z ) F ( λ z ) , z R ,
where f ( · ) is a symmetric pdf, F ( · ) is an absolutely continuous symmetric cumulative distribution function (cdf) and λ R is a shape parameter. A particular unimodal case of density function given in (1) is the skew-normal (SN) distribution [1], which has pdf given by
f S N ( z ; λ ) = 2 ϕ ( z ) Φ ( λ z ) , z R ,
where ϕ ( · ) and Φ ( · ) are the pdf and cdf of the standard normal distribution, respectively. Among the densities of bimodal type supported under general structure (1), it can be stood out the proposals of Arnold et al. [2], Gómez et al. [3] and Kim [4]. Based on the alpha-power family of distributions by Durrans [5] with pdf given by
f A P ( z ; α ) = α f ( z ) F ( z ) α , z R ,
where F ( · ) is a cdf absolutely continuous, f = d F and α R + is a shape parameter, Bolfarine et al. [6] studied a new bimodal distribution. Other proposals for fitting data with bimodal behavior have been considered by Elal-Olivero [7], Elal-Olivero et al. [8] and Ma and Genton [9]. For example, Ma and Genton [9] introduced the family of distributions with pdf
f ( z ) = 2 ϕ ( z ) Φ w ( z ) , z R ,
where ϕ ( · ) and Φ ( · ) are the pdf and cdf of the standard normal distribution, respectively; and w ( z ) is a polynomial of odd order. Elal-Olivero et al. [8] defined the bimodal elliptical skew-normal (BESN) model by multiplying the pdf of the SN distribution by a polynomial of even order resulting the pdf given by
f B E S N ( z ; α ) = 2 1 + α z 2 1 + α ϕ ( z ) Φ ( λ z ) , z R ,
where λ R , is an asymmetry parameter, and α 0 is a shape parameter. They showed that the model fits data of bimodal type for α > 0 . The case α = 0 is reduced to the SN distribution and for α = λ = 0 , the standard normal distribution is obtained. An asymmetric extension of the bimodal-normal model which has pdf given by
f B N ( z ) = z 2 ϕ ( z ) , z R ,
was proposed by Elal-Olivero [7]. This extension is called the alpha-skew-normal (ASN) distribution and its pdf is
f A S N ( z ; α ) = ( 1 α z ) 2 + 1 2 + α 2 ϕ ( z ) , z R ,
for α R . This model fits data with bimodal shape for α 0 , and for α = 0 the model is reduced to the normal distribution. Properties of the model and the statistical inference for the parameters can be seen in Elal-Olivero [7]. Recently Shafiei et al. [10] presented a generalization of the ASN model by adding a parameter β resulting a more general and flexible model. This addition allows to modeling data sets with the possibility of until four modes and model is denominated alpha-beta skew-normal (ABSN). The pdf of the ABSN model is given by
f A B S N ( z ; α ) = ( 1 α z β z 3 ) 2 + 1 α 2 + 15 β 2 + 6 α β + 2 ϕ ( z ) , x R ,
where α ,   β R . Note that, for β = 0 the ASN model of Elal-Olivero [7] is obtained, and for α = β = 0 the normal distribution is followed.
One of the great disadvantages of the widely known SN model by Azzalini [1], is that the information matrix is singular when λ = 0 . This same characteristic is presented by the BESN model for the case α = λ = 0 . The bimodal models by Arnold et al. [2] and Gómez et al. [3] also present the problem of singularity of the information matrix. Elal-Olivero [7] shows that ASN distribution has singular information matrix when the shape parameter α is equal to zero, and given that ABSN model of Shafiei et al. [10] contains the ASN model as a particular case, it can be shown that for α = β = 0 , the ABSN model also has singular information matrix, a situation also presented by the proposal of Ma and Genton [9]. In this way, models of bimodal type that result from including polynomials of even order in the density of normal or skew-normal model; or polynomials of odd order in the argument of skew-normal model, acquire the same characteristic of the skew-normal model in its information matrix; situation that also happens with those models of bimodal type that contain this model as a special case, as it happens with the proposals of Arnold et al. [2] and Gómez et al. [3].
On the other hand, Pewsey et al. [11] studied the alpha-power (AP) model given in (3) for the special case F ( · ) = Φ ( · ) , which is denominated the power-normal (PN) model and they found that α values near to one (greater and smaller than one), the information matrix is non-singular, being this an advantage in the inferential process and the asymptotic properties of the maximum likelihood estimator (MLE). Further, Martínez-Flórez et al. [12] extended the SN distribution to the alpha-power family of distributions, and they obtained a more flexible model in terms of asymmetry and kurtosis than SN and PN models. The bimodal model based on the PN distribution of Bolfarine et al. [6] also presents the same property of information matrix non-singular. Hence, given the characteristics of the information matrix of the alpha-power family of distributions, in this paper we propose an extension alpha-power of a multimodal model by using a polynomial of even degree. The introduced model is a more flexible distribution in terms of asymmetry and kurtosis than those existing in the literature and with information matrix non-singular.
The rest of this paper is organized as follows—Section 2 presents the beta-skew alpha-power model and discusses its main properties. In particular, we show how bimodality and trimodality shape are obtained. We consider a location-scale family and the inference process is carried out by using maximum likelihood method. In addition, some properties of the special case of the beta-skew-normal model are studied in details. In Section 3 a small Monte Carlo simulation is presented. In Section 4 two real data application are reported and compares it with several rival models.

2. The Asymmetric Beta-Skew Alpha-Power Distribution

In this section, we introduce a new multimodal asymmetric distribution by considering the asymmetric beta-skew-normal (BSN) model and by incorporating an additional parameter. The new family of distributions extend the usual normal model and others distributions are also particular cases from this model.
Definition 1.
The random variable Z is said to have a beta-skew alpha-power distribution, which we will denote as Z BSAP ( β , α ) , if Z has the following pdf
f B S A P ( z ; β , α ) = α f B S N ( z ; β ) F B S N ( z ; β ) α 1 , z R ,
for α > 0 and β R . The functions f B S N ( z ; β ) and F B S N ( z ; β ) are the pdf and cdf of the BSN distribution given by
f B S N ( z ; β ) = 1 β z 3 2 + 1 2 + 15 β 2 ϕ ( z )
and
F B S N ( z ; β ) = Φ ( z ) + 4 β 15 β 2 z + 2 β z 2 5 β 2 z 3 β 2 z 5 2 + 15 β 2 ϕ ( z ) ,
respectively. Here, ϕ ( · ) and Φ ( · ) are de pdf and cdf of the standard normal distribution, respectively.
Figure 1 depicts the shape of the beta-skew alpha-power (BSAP) distribution for some selected values of α and β parameters. It can see from the graph that BSAP distribution has unimodal, bimodal and multimodal (three modes) behavior. Notice that, when β increases, the density function takes a bimodal shape. The graphs in the figure also show a bimodal shape for great values of α .
If Z BSAP ( β , α ) , the following properties are deduced immediately from the definition
(i)
If α = 1 , then Z BSN ( β ) .
(ii)
If β = 0 , then Z PN ( α ) .
(iii)
If α = 1 and β = 0 , then Z N ( 0 , 1 ) .
Proposition 1.
The density function (9) has at most three modes.
Proof. 
Given that α is an asymmetry parameter, without loss of generality we take α = 1 (the case of the BSN distribution) and by differentiating Equation (10) we obtain
f B S N ( z ; β ) = z ϕ ( z ) 2 + 15 β 2 ( β 2 z 6 6 β 2 z 4 2 β z 3 + 6 β z + 2 ) ,
thus, f B S N ( z ; β ) has at most seven zeros. It can be seen that for 0 < | β | < + , the value z = 0 is a zero of the f B S N ( z ; β ) function and it can be shown that f B S N ( 0 ; β ) = 2 / π / ( 2 + 15 β 2 ) < 0 , therefore, at z = 0 occurs a maximum when 0 < | β | < + . By analyzing the polynomial of the BSN model, we making g ( z ) = ( 1 β z 3 ) 2 + 1 = β 2 z 6 2 β z 3 + 2 , hence g ( z ) = β 2 z 6 + 2 β z 3 + 2 . For β > 0 the number of positive roots of the polynomial g ( z ) is two, while g ( z ) does not change of sign, hence, the polynomial does not have negative roots. It is concluded that the polynomial has at most four complex roots. For β < 0 , the polynomial would have two or no negative roots and four complex roots. Using computational methods, it can be shown that for values of β [ 0 , 1000 ] the resulting polynomial from the derivative β 2 z 6 6 β 2 z 4 2 β z 3 + 6 β z + 2 has at least two non-real complex roots, therefore, it would have at most four real roots, and therefore it is concluded that at most they would have three maximums in the model. □
Proposition 2.
Let Z BSAP ( β , α ) ,
(i) 
The cdf of Z, which we denote by F B S A P ( t , β , α ) , is given by
F B S A P ( z , β , α ) = Φ ( z ) + 4 β 15 β 2 z + 2 β z 2 5 β 2 z 3 β 2 z 5 2 + 15 β 2 ϕ ( z ) α .
(ii) 
The survival function, denoted by S B S A P ( z ) , is
S B S A P ( z , β , α ) = 1 Φ ( z ) + 4 β 15 β 2 z + 2 β z 2 5 β 2 z 3 β 2 z 5 2 + 15 β 2 ϕ ( z ) α .
(iii) 
The Hazard function h B S A P ( z , β , α ) , is
h B S A P ( z , β , α ) = f B S A P ( z , β , α ) S B S A P ( z , β , α ) .
Figure 2 shows the shape of the survival function for some selected values of α and β . In these graphs, it can be seen that the curve becomes increasingly horizontal in the interval ( 0 , 1 ) as β increases, and the probability of survival is greater for larger α values when β is constant. It is also important to note that, regardless of the values of β and α , the survival approaches to zero when z tends to infinite.
Proposition 3.
Let Z BSAP ( β , α ) , then the k-th moment of Z is given by
E Z k = α 2 + 15 β 2 + 2 z r 2 β z r + 3 + β 2 z r + 6 ϕ ( z ) F B S N ( z ; β , α ) α 1 d z .
Proof. 
The proof is direct from definition of expected value. □
The expected value, the variance and the indices of skewness and kurtosis of the BSAP model, which are denoted by μ , σ 2 , γ 1 and γ 2 respectively, can be obtained from (16) by using the following expressions
(i)
μ = μ 1 .
(ii)
σ 2 = μ 2 μ 1 2 .
(iii)
γ 1 = μ 3 3 μ 1 μ 2 + 2 μ 1 3 μ 2 μ 1 2 3 / 2 .
(iv)
γ 2 = μ 4 4 μ 1 μ 3 + 6 μ 1 2 μ 2 3 μ 1 4 μ 2 μ 1 2 2 3 .
Here μ k = E [ Z k ] for k = 1 , 2 , 3 , 4 .
Remark 1.
We calculated γ 1 and γ 2 by using (16) and numerical integration for the BSAP ( β , α ) model for α ( 0.01 , 1000 ) and β ( 50 , 50 ) . We obtain
3.17 < γ 1 < 1.15 and 0.00014 < γ 2 < 31.18 .
Notice that, the length of the admissible intervals for the skewness and the kurtosis parameters of the BSAP distribution are larger than the corresponding intervals of the ASN, the SN, the PN and the skew-normal alpha-power (SNAP) distributions, which are ( 0.811 , 0.811 ) , ( 1.300 , 0.7489 ) ; ( 0.995 , 0.995 ) , ( 3.000 , 3.869 ) ; ( 0.6115 , 0.9007 ) , ( 1.283 , 1.3556 ) and ( 0.14676 , 0.9953 ) , ( 1.5328 , 2.4386 ) , respectively. See References [1,7,11,12] for more details.
Some additional properties for the special case when α = 1 are presented to follow.
Lemma 1.
Let Z BSN ( β ) , then
E Z 2 k = 2 + β 2 ( 2 k + 1 ) ( 2 k + 3 ) ( 2 k + 5 ) 2 + 15 β 2 j = 1 k ( 2 j 1 ) , for k = 1 , 2 ,
E Z 2 k 1 = 2 β ( 2 k + 1 ) 2 + 15 β 2 j = 1 k ( 2 j 1 ) , for k = 1 , 2 ,
Proof. 
Knowing that E Y 2 k = j = 1 k 2 j 1 and E Y 2 k 1 = 0 , for k = 1 , 2 , when Y N 0 , 1 , we obtain
E Z 2 k = 1 2 + 15 β 2 2 E Y 2 k + β 2 E Y 2 ( k + 3 ) , with Y N ( 0 , 1 ) = 1 2 + 15 β 2 2 j = 1 k ( 2 j 1 ) + β 2 j = 1 k + 3 ( 2 j 1 ) = 1 2 + 15 β 2 2 j = 1 k ( 2 j 1 ) + β 2 j = 1 k ( 2 j 1 ) ( 2 k + 1 ) ( 2 k + 3 ) ( 2 k + 5 ) = 2 + β 2 ( 2 k + 1 ) ( 2 k + 3 ) ( 2 k + 5 ) 2 + 15 β 2 j = 1 k ( 2 j 1 ) ,
E Z 2 k 1 = 1 2 + 15 β 2 2 E Y 2 k 1 2 β E Y 2 ( k + 1 ) + β 2 E Y 2 ( k + 3 ) 1 = 1 2 + 15 β 2 2 β E Y 2 ( k + 1 ) = 1 2 + 15 β 2 2 β j = 1 k + 1 ( 2 j 1 ) = 2 β ( 2 k + 1 ) 2 + 15 β 2 j = 1 k ( 2 j 1 ) .
 □
Lemma 2.
Let μ, σ 2 , γ 1 and γ 2 the mean, variance and the indices of skewness and kurtosis, respectively of the BSN ( β ) model,
(i) 
μ = 6 β 2 + 15 β 2 .
(ii) 
σ 2 = 1575 β 4 + 204 β 2 + 4 2 + 15 β 2 2 .
(iii) 
γ 1 = 21600 β 5 + 2088 β 3 48 β 1575 β 4 + 204 β 2 + 4 3 / 2 .
(iv) 
γ 2 = 3189375 β 8 + 1474200 β 6 + 182952 β 4 + 6624 β 2 + 48 4 + 204 β 2 + 1575 β 4 2 .

2.1. Location and Scale Extension for BSAP Model

We can also consider a generalization of the BSAP ( β , α ) distribution by adding location and scale parameters. Next definition gives the generalization of the BSAP model.
Definition 2.
Let Z BSAP ( β , α ) . The beta-skew alpha-power density of location and scale is defined as the distribution of X = ξ + η Z , for ξ R and η > 0 . The corresponding density function is given by
f B S A P ( x ; θ ) = α ( 1 β z 3 ) 2 + 1 ϕ ( z ) σ ( 2 + 15 β 2 ) × Φ ( z ) + 4 β 15 β 2 z + 2 β z 2 5 β 2 z 3 β 2 z 5 2 + 15 β 2 ϕ ( z ) α 1 ,
where z = ( x ξ ) / η and θ = ( ξ , η , β , α ) . We denote this as X BSAP ( θ ) .
The distribution function associated to the density (20) is given by
F B S A P ( x ; θ ) = Φ ( z ) + 4 β 15 β 2 z + 2 β z 2 5 β 2 z 3 β 2 z 5 2 + 15 β 2 ϕ ( z ) α ,
and the k-th moment of the random variable X is
E X k = j = 0 k k j ξ k j η j E Z j ,
where Z BSAP ( β , α ) .

2.2. Maximum Likelihood Estimation for BSAP Model

We consider a random sample X = ( X 1 , , X n ) of size n, from the BSAP ( θ ) distribution, where  θ = ξ , η , β , α . The log-likelihood function is given by
( θ ; X ) = n log α n log η n log 2 + 15 β 2 n log ( 2 π ) + i = 1 n log 1 β z i 3 2 + 1 1 2 i = 1 n z i 2 + α 1 i = 1 n log F B S N z i ,
which is a continuos function in each parameter. Thus, by differentiating the log-likelihood function we obtain the following likelihood equations
( θ ; X ) ξ = i = 1 n 6 β 1 β z i 3 z i 2 1 β z i 3 2 + 1 + z i ( α 1 ) i = 1 n f B S N z i F B S N z i = 0 ,
( θ ; X ) η = i = 1 n 6 β 1 β z i 3 z i 3 1 β z i 3 2 + 1 1 + z i 2 ( α 1 ) i = 1 n z i f B S N z i F B S N z i = 0 ,
( θ ; X ) β = 2 i = 1 n z i 3 1 β z i 3 1 β z i 3 2 + 1 + 15 β 2 + 15 β 2 + α 1 i = 1 n W i f B S N z i F B S N z i = 0 ,
( θ ; X ) α = n α + i = 1 n log F B S N ( z i ) = 0 .
where W i = ( 8 60 β 2 60 β z i + 4 z i 2 30 β 2 z i 2 20 β z i 3 4 β z i 5 ) / ( 2 + 15 β 2 ) 2 . The solutions of likelihood Equations (21)–(24) provide the maximum likelihood estimators (MLEs) of ξ , η , α and β , which can be obtained by numerical method such as Newton-Rapshon type procedure. Under certain regularity conditions, the elements of the Fisher information may be calculated as
i θ p θ q = E 2 ( θ ; X ) θ p θ q , for p , q = 1 , 2 , 3 , 4 .
The Cramér–Rao bound states that the inverse of the Fisher information is a lower bound on the variance of any unbiased estimator. Thus, we can find a lower bound for the standard errors (SE) of the MLEs as the square root of the diagonal elements of the observed Fisher information matrix. The observed information matrix can be obtained by taking the second partial derivatives of the log-likelihoiod function and multiplying by −1, that is,
j θ p θ q = 2 ( θ ; X ) θ p θ q , for p , q = 1 , 2 , 3 , 4 ,
where θ 1 = ξ , θ 2 = η , θ 3 = β , and θ 4 = α . The elements of the expected and observed Fisher information | I ( θ ) | and | J ( θ ) | , respectively, are given in the Appendix A. It can be shown that, for  β = 0 and α = 1 (the normal distribution case), we have that | I ( θ ) | 0 , therefore, I ( θ ) is non-singular. Hence, covariance matrix of the parameter vector is the inverse of the Fisher information matrix, that is, Var ( θ ^ ) = I 1 ( θ ) . Thus, the MLE θ ^ converges to a normal distribution, that is
θ ^ d N 4 ( θ , I 1 ( θ ) ) .
Remark 2.
The Fisher information matrix for the BSN model is given by
I ( θ ) = 36 β 2 d 4 75 β 2 + 2 σ 2 ( 2 + 15 β 2 ) 36 β 2 d 5 + 6 β σ 2 ( 2 + 15 β 2 ) 12 β d 5 6 σ ( 2 + 15 β 2 ) 36 β 2 d 5 + 6 β σ 2 ( 2 + 15 β 2 ) 36 β 2 d 6 330 β 2 + 4 σ 2 ( 2 + 15 β 2 ) 12 β ( 15 d 6 ) σ 2 ( 2 + 15 β 2 ) 12 β d 5 6 σ ( 2 + 15 β 2 ) 12 β ( 15 d 6 ) σ 2 ( 2 + 15 β 2 ) 60 β 2 d 6 + 8 d 6 900 β 2 ( 2 + 15 β 2 ) 2
where d k = E Z k ( 1 β Z 3 ) 2 / ( 1 β Z 3 ) 2 + 1 , for k = 4 , 5 , 6 , with Z N 0 , 1 , In particular, if  β = 0  then
I ( θ ) = 1 σ 2 0 3 σ 0 2 σ 2 0 3 σ 0 15
whose determinant is I ( θ ) = 12 / σ 4 0 , hence, I ( θ ) is non-singular and Var ( θ ^ ) = I 1 ( θ ) . Hence, it follows that the MLE θ ^ converges to the normal distribution
θ ^ d N 3 ( θ , I 1 ( θ ) ) .

3. Simulation Study

In order to study the performance of the MLEs of the parameters in BSAP model, we conducted a Monte Carlo simulation study with samples sizes n = 40 , 80, 120, 160 and 320. The true values of the parameters were taken as ξ = 0 , η = 1 , α = 0.5, 2.0 and 5.0; and β = 1.0, 1.75 and 3.0. For each combination of parameters and sample sizes, we generated 5000 samples from the BSAP model. To evaluate estimators performance were considered the absolute value of the bias (B), and the squared root of the mean squared error (RMSE). They are given by
B ( θ ^ i ) = 1 5000 j = 1 5000 θ ^ i ( j ) θ i ,   RMSE ( θ ^ i ) = 1 5000 j = 1 5000 θ ^ i ( j ) θ i 2 ,
respectively, where θ ^ i ( j ) is the estimate of θ i for the j-th sample, for θ i θ = ( ξ , η , β , α ) . MLEs were computed by using optim function in R Development Core Team [13].
From Table 1 we can see that, as the sample sizes increase, the bias (in absolute value) and the squared root of the mean squared error decrease, indicating a good behavior of the MLEs of the parameters in BSAP model. Then, it follows that for large sample sizes, MLEs are asymptotically consistent.

4. Real Data Applications

In this section, we illustrate the proposed model by considering two real data sets. In the first application we consider the data on the otis IQ scores of 52 non-white males hired by a large insurance company in 1971. In the second application we use the geyser data set, available in R Development Core Team [13].

4.1. Application 1: The Otis IQ Scores Data

These data set have been analyzed previously by Gupta and Gupta [14] and Sharafi and Behboodian [15]. See, Roberts [16] for more details. Table 2 shows the statistic summary for the otis scores data. The results shows that the data present positives skewness and kurtosis lower than the normal model, likewise, the histogram in Figure 3 shows that the data have more than one mode.
We fit the BSN and BSAP models to analyze this data set. To compare the proposed model, we also fit the multimodal alpha-beta skew-normal (ABSN) model by Shafiei et al. [10], and the asymmetric bimodal model (ETN) of Arnold et al. [2]. The fit of these models is carried out by using the maximum likelihood method and optim function of [13]. Table 3 shows the parameter estimates, together with their corresponding standard errors (SE). To compare fitted models, we use the AIC Akaike [17], corrected CAIC, BIC by Hastie and Tibshirani [18] and the HQIC or information criterion of Hannan and Quinn [19], namely
AIC = 2 ^ ( · ) + 2 p , CAIC = AIC + 2 p ( p + 1 ) n ( p + 1 ) , BIC = 2 ^ ( · ) + p log ( n )
and
HQIC = AIC 2 p + 2 p log ( log ( n ) ) ,
where p is the number of parameters in the considered model. According to any of these criteria, the BSAP and BSN seem to provide better fit to the otis IQ scores data than the ABSN and ETN models.
We can use the likelihood ratio (LR) test statistic to confirm the use the BSAP model instead of the BSN model, so we consider the following hypotheses,
H 0 : α = 1 ( BSN ( ξ , η , β ) ) vs H 1 : α 1 ( BSAP ( ξ , η , β , α ) ) ,
with LR test statistic,
Λ = L B S N ( θ ^ ; X ) L B S A P ( θ ^ ; X )
where L B S N ( · ) and L B S P N ( · ) denote the likelihood functions of the BSN and BSAP models, respectively. We obtain the value 2 log ( Λ ) = 2 ( B S N ( θ ^ ; X ) B S A P ( θ ^ ; X ) ) = 4.0166 , and comparing this quantity with χ 1 2 = 3.841 , the null hypotheses is rejected, that is, the BSAP model is more flexible than the asymmetric BSN model, taking into account the test results of Hartigan and Hartigann [20] and Hartigan [21].
Figure 3a,b show the behavior of the fitted models and the empirical cdf for the ETN, BSN and BSAP models. It can be seen from the figure that the BSAP model has the best fit against the ABSN, BSN and ETN models, while the BSN model has a better fit than the model ETN. Also, the graphs in Figure 3c,d show the QQplots for the BSAP and BSN models.

4.2. Application 2: Old Faithful Geyser Data

For the second application, we consider a data set consisting of 272 observations about the wait times between the eruptions (in minutes) of the old faithful geyser in Yellowstone National Park, Wyoming, U.S. Data set are available in the libraries stats and MASS of R Development Core Team [13]. More information of these data can be seen in Azzalini and Bowman [22] who takes a look at this set of data. Table 4 shows the summary statistic for data set. The results shows that the data set present negative asymmetry and kurtosis below the normal model.
In addition to the BSAP and BSN models, we also fit the ETN and the alpha-skew-normal (ASN) models, see Elal-Olivero [7]. The results of the fit of these models can be seen in Table 5. The standard errors of the estimators were obtained by using the observed information matrix, and again, to compare the fited models, we use the AIC, CAIC, BIC and CAIC criteria. According to any of these criteria, BSAP model seems to provide better fit to the old faithful gayser data than BSN, ASN and ETN models.
We tested the hypotheses
H 0 : α = 1 ( BSN ( ξ , η , β ) ) vs H 1 : α 1 ( BSAP ( ξ , η , β , α ) ) ,
with LR test statistic,
Λ = L B S N ( θ ^ ; X ) L B S A P ( θ ^ ; X )
where L B S N ( · ) and L B S P N ( · ) denote the likelihood functions of the BSN and BSAP models, respectively. We obtain the value 2 log ( Λ ) = 2 ( B S N ( θ ^ ; X ) B S A P ( θ ^ ; X ) ) = 8.884 , and comparing this quantity with χ 1 2 = 3.841 , the null hypotheses is rejected, that is, the BSAP model is more flexible than the asymmetric BSN model.
Figure 4 shows the behavior of the fitted models and the empirical cdf for the ETN, BSN and BSAP models. It can be seen that the BSAP model has the best fit against the BSN and ETN models, while the BSN model has a better fit than the ETN model. The graphs in Figure 4 also show the QQplot for the BSAP model.

5. Conclusions

In this paper, a new class of unimodal, as well as bimodal and trimodal, skew distribution was proposed. The main statistical properties of the model and the problem of the parameters estimation are studied in details by using maximum likelihood method. The model extends the usual normal distribution to trimodal asymmetric case and the BSN model is also a special case. Furthermore, we have shown that such distribution is more flexible than certain rival models and it fits better to some real data sets.

Author Contributions

All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

The research of G. Martínez-Flórez and R. Tovar-Falón was supported by Grant Proyecto Universidad de Córdoba: Familia bimodal de distribuciones de probabilidad skew-normal alpha-potencia, Code FCA-18-15 (Colombia).

Acknowledgments

The authors wish to acknowledge the Associate Editor and two anonymous referees for the constructive comments that helped to improve the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Information Matrix for the BSAP Model

In this section, expressions for the elements of the observed and expected information matrix of the BSAP model are provided. The observed elements are denoted by j ξ ξ , j ξ η , i ξ β , i ξ α , j η η , j η β , j η α , j β β , j β α and j α α , and they can be calculated by using
j θ p θ q = 2 ( θ ; X ) θ p θ q , for p , q = 1 , 2 , 3 , 4 .
After some algebraic manipulations we obtain
j ξ ξ = 1 η 2 i = 1 n 6 β z i 5 β z i 3 2 1 β z i 3 2 + 1 36 β 2 z i 4 1 β z i 3 2 1 β z i 3 2 + 1 2 1 + α 1 η i = 1 n 6 β z i 2 1 β z i 3 2 + 15 β 2 ϕ z i F B S N z i + z i f B S N z i F B S N z i + f B S N 2 z i F B S N 2 z i , j ξ η = 2 η 2 i = 1 n 9 β z i 2 2 β z i 3 1 1 β z i 3 2 + 1 18 β 2 z i 5 1 β z i 3 2 1 β z i 3 2 + 1 2 z i + α 1 η 2 i = 1 n 6 β z i 3 1 β z i 3 2 + 15 β 2 ϕ z i F B S N z i + z i 2 1 f B S N z i F B S N z i + z i f B S N 2 z i F B S N 2 z i ,
j ξ β = 6 η i = 1 n z i 2 1 2 β z i 3 1 β z i 3 2 + 1 + 2 β z i 5 1 β z i 3 2 1 β z i 3 2 + 1 2 α 1 σ i = 1 n 2 z i 3 1 β z i 3 2 + 15 β 2 ϕ z i F B S N z i + 30 β 2 + 15 β 2 f B S N z i F B S N z i + W i ϕ z i f B S N z i F B S N 2 z i ,
j ξ α = 1 η i = 1 n f B S N z i F B S N z i , j η η = 1 η 2 i = 1 n 2 β z i 3 21 β z i 3 12 1 β z i 3 2 + 1 36 β 2 z i 6 1 β z i 3 2 1 β z i 3 2 + 1 2 + 1 3 z i 2 + α 1 η 2 i = 1 n 6 β z i 4 1 β z i 3 2 + 15 β 2 ϕ z i F B S N z i + z i 3 2 z i f B S N z i F B S N z i + z i 2 f B S N 2 z i F B S N 2 z i , j η β = 6 η i = 1 n z i 3 1 2 β z i 3 1 β z i 3 2 + 1 + 2 β z i 6 1 β z i 3 2 1 β z i 3 2 + 1 2 α 1 η i = 1 n [ 2 z i 4 1 β z i 3 2 + 15 β 2 ϕ z i F B S N z i + 30 β z i 2 + 15 β 2 f B S N z i F B S N z i + W i z i ϕ z i f B S N z i F B S N 2 z i ] ,
j η α = 1 η i = 1 n z i f B S N z i F B S N z i , j β β = 2 i = 1 n z i 6 1 β z i 3 2 + 1 2 z i 6 1 β z i 3 2 1 β z i 3 2 + 1 2 15 2 15 β 2 2 + 15 β 2 2 α 1 i = 1 n ϕ z i U i F B S N z i W i 2 ϕ z i F B S N 2 z i , j β α = i = 1 n W i ϕ z i F B S N z i , j α α = n α 2 ,
where
U i = 720 β + 1800 β 3 120 z i + 2700 β 2 z i 360 β z i 2 + 900 β 2 z i 2 40 z i 3 + 900 β 2 z i 3 8 z i 5 + 180 β 2 z i 5 , W i = 8 60 β 2 60 β z i + 4 z i 2 30 β 2 z i 2 20 β z i 3 4 β z i 5 2 + 15 β 2 2 .
The elements of the expected Fisher information matrix can be obtained by using
i θ p θ q = E 2 ( θ ; X ) θ p θ q , for p , q = 1 , 2 , 3 , 4 .
and letting
h j = E Z j ϕ z F B S N z , v j k = E Z j f B S N z F B S N z k , g j k = E Z j 1 β z 3 2 + 1 k , u j k = E Z j W k ϕ z f B S N z F B S N 2 z , a = E U ϕ z F v B S N z , and b k = E W k ϕ k z F B S N k z ,
where Z BSAP ( β , α ) ; we have
i ξ ξ = α 1 η 6 β 2 + 15 β 2 h 2 β h 5 + v 11 + v 02 1 η 2 6 β 5 β g 4 , 1 2 g 1 , 1 36 β 2 g 4 , 2 2 β g 7 , 2 + β 2 g 10 , 2 1 , i ξ η = α 1 η 2 6 β 2 + 15 β 2 h 3 β h 6 + v 21 v 01 + v 12 2 η 2 9 β 2 β g 5 , 1 g 2 , 1 18 β 2 g 5 , 2 2 β g 8 , 2 + β 2 g 11 , 2 g 1 , 0 , i ξ β = α 1 η 2 2 + 15 β 2 h 2 β h 6 + 30 β 2 + 15 β 2 v 01 + u 01 6 η g 2 , 1 2 β g 3 , 1 + 2 β g 5 , 2 2 β g 8 , 2 + β 2 g 11 , 2 , i ξ α = 1 η v 01 , i η η = α 1 η 6 β 2 + 15 β 2 h 4 β h 7 + v 31 2 v 11 + v 22 1 η 2 β 21 β g 6 , 1 12 g 3 , 1 36 β 2 g 6 , 2 2 β g 9 , 2 + β 2 g 12 , 2 + 1 3 g 1 , 0 , i η β = α 1 η 2 2 + 15 β 2 h 4 β h 7 + 30 β 2 + 15 β 2 v 11 + u 11 6 η g 3 , 1 2 β g 6 , 1 + 2 β g 6 , 2 2 β g 9 , 2 + β 2 g 12 , 2 ,
i β β = α 1 a b 1 2 g 6 , 1 2 g 6 , 2 2 β g 9 , 2 + β 2 g 12 , 2 15 2 15 β 2 2 + 15 β 2 2 , i β α = b 1 , i α α = 1 α 2 , i η α = 1 η v 11 .

References

  1. Azzalini, A. A class of distributions which includes the normal ones. Scand. J. Stat. 1985, 12, 171–178. [Google Scholar]
  2. Arnold, B.C.; Gómez, H.W.; Salinas, H.S. On multiple constraint skewed models. Statistics 2009, 43, 279–293. [Google Scholar] [CrossRef]
  3. 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. [Google Scholar] [CrossRef]
  4. Kim, H.J. On a class of two–piece skew-normal distributions. Statistics 1985, 39, 537–553. [Google Scholar] [CrossRef]
  5. Durrans, S.R. Distributions of fractional order statistics in hydrology. Water Resour. Res. 1992, 28, 1649–1655. [Google Scholar] [CrossRef]
  6. Bolfarine, H.; Martínez-Flórez, G.; Salinas, H.S. Bimodal symmetric–asymmetric power–normal families. Commun. Stat. Theory Methods 2018, 47, 259–276. [Google Scholar] [CrossRef]
  7. Elal-Olivero, D. Alpha-skew-normal distribution. Proyecc. J. Math. 2010, 29, 224–240. [Google Scholar] [CrossRef] [Green Version]
  8. Elal-Olivero, D.; Gómez, H.W.; Quintana, F.A. Bayesian modeling using a class of bimodal skew–elliptical distributions. J. Stat. Plan. Inference 2010, 139, 1484–1492. [Google Scholar] [CrossRef]
  9. Ma, Y.; Genton, M.G. Flexible Class of Skew–Symmetric Distributions. Scand. J. Stat. 2004, 31, 459–468. [Google Scholar] [CrossRef] [Green Version]
  10. Shafiei, S.; Doostparast, M.; Jamalizadeh, A. The alpha–beta skew normal distribution: Properties and applications. Statistics 2016, 50, 338–349. [Google Scholar] [CrossRef]
  11. Pewsey, A.; Gómez, H.W.; Bolfarine, H. Likelihood–based inference for power distributions. Test 2012, 21, 775–789. [Google Scholar] [CrossRef]
  12. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. Skew–normal alpha-power model. Statistics 2014, 48, 1414–1428. [Google Scholar] [CrossRef]
  13. R Development Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019; Available online: http://www.R-project.org (accessed on 10 December 2019).
  14. Gupta, R.C.; Gupta, R.D. Generalized skew normal model. Test 2004, 13, 501–524. [Google Scholar] [CrossRef]
  15. Sharafi, M.; Behboodian, J. The Balakrishnan skew–normal density. Stat. Papers 2008, 49, 769–778. [Google Scholar] [CrossRef]
  16. Roberts, H. Data Analysis for Managers with MINITAB, 1st ed.; Scientific Press: Redwood City, CA, USA, 1988. [Google Scholar]
  17. Akaike, H. A new look at statistical model identification. IEEE Trans. Autom. Contr. 1974, 19, 716–722. [Google Scholar] [CrossRef]
  18. Hastie, T.J.; Tibshirani, R.J. Generalized Additive Models, 1st ed.; Chapman and Hall/CRC: New York, NY, USA, 1990. [Google Scholar]
  19. Hannan, E.J.; Quinn, B.G. The determination of the order of an autoregression. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 1979, 41, 190–195. [Google Scholar] [CrossRef]
  20. Hartigan, J.A.; Hartigan, P.M. The dip test of unimodality. Ann. Stat. 1985, 13, 70–84. [Google Scholar] [CrossRef]
  21. Hartigan, P.M. Algorithm AS 217: Computation of the dip statistic to test for unimodality. J. R. Stat. Soc. Ser. C (Appl. Stat.) 1985, 34, 320–325. [Google Scholar] [CrossRef]
  22. Azzalini, A.; Bowman, A.W. A look at some data on the old faithful geyser. J. R. Stat. Soc. Ser. C (Appl. Stat.) 1990, 39, 357–365. [Google Scholar] [CrossRef]
Figure 1. Density function of BSAP ( β , α ) distribution. (a) BSAP ( 0.25 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (b) BSAP ( 0.75 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (c) BSAP ( 1.5 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line).
Figure 1. Density function of BSAP ( β , α ) distribution. (a) BSAP ( 0.25 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (b) BSAP ( 0.75 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (c) BSAP ( 1.5 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line).
Symmetry 12 00613 g001
Figure 2. Survival function of BSAP ( β , α ) distribution. (a) BSAP ( 0.25 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (b) BSAP ( 0.75 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (c) BSAP ( 1.5 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line).
Figure 2. Survival function of BSAP ( β , α ) distribution. (a) BSAP ( 0.25 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (b) BSAP ( 0.75 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line). (c) BSAP ( 1.5 , α ) for α = 0.75 (dotted-dashed line), α = 1.0 (dotted line), α = 2.0 (dashed line) and α = 3.0 (solid line).
Symmetry 12 00613 g002
Figure 3. (a) Histogram for the otis data set. BSAP model (solid line), BSN model (dashed line), ETN model (dotted line). (b) Empiric cdf for BSAP model (dashed line), BSN (dotted lines) and ETN model (dotted-dashed line). (c) QQplot for fitted BSAP model. (d) QQplot fitted BSN model.
Figure 3. (a) Histogram for the otis data set. BSAP model (solid line), BSN model (dashed line), ETN model (dotted line). (b) Empiric cdf for BSAP model (dashed line), BSN (dotted lines) and ETN model (dotted-dashed line). (c) QQplot for fitted BSAP model. (d) QQplot fitted BSN model.
Symmetry 12 00613 g003
Figure 4. (a) Histogram for the old faithful geyser data. BSAP model (solid line), BSN model (dotted line), ETN model dashed( line). (b) Empiric cdf for BSAP model (dashed line), BSN (dotted lines) and ETN model (dotted-dashed line). (c) QQplot for fitted BSAP model.
Figure 4. (a) Histogram for the old faithful geyser data. BSAP model (solid line), BSN model (dotted line), ETN model dashed( line). (b) Empiric cdf for BSAP model (dashed line), BSN (dotted lines) and ETN model (dotted-dashed line). (c) QQplot for fitted BSAP model.
Symmetry 12 00613 g004aSymmetry 12 00613 g004b
Table 1. Asymptotic behavior of the MLEsfor the BSAP ( ξ , η , β , α ) model.
Table 1. Asymptotic behavior of the MLEsfor the BSAP ( ξ , η , β , α ) model.
ξ ^ η ^ α ^ β ^
α β n | B | RMSE | B | RMSE | B | RMSE | B | RMSE
0.51.0400.01400.26470.01200.08082.28316.19430.00780.1752
800.01050.17910.00530.05800.63563.41200.00610.1259
1200.00650.14330.00400.04610.25662.29420.00170.0940
1600.00360.12120.00180.03980.09960.73870.00240.0803
3200.00100.08350.00110.02700.03320.16500.00140.0557
1.75400.01970.21700.00930.07286.46259.35730.00470.1543
800.01650.15180.00640.05153.96088.25400.00230.1078
1200.01130.12700.00460.04162.67147.64730.00180.0884
1600.00970.10740.00260.03551.65306.17340.00010.0751
3200.00440.07480.00140.02540.37832.30640.00010.0532
3.0400.01520.19380.00850.06859.120211.77130.00630.1449
800.01460.13710.00480.04798.216711.02280.00230.1007
1200.01250.11270.00340.03897.617210.90910.00240.0819
1600.01110.09860.00290.03346.490710.21520.00150.0718
3200.00940.07090.00140.02403.41558.92610.00130.0504
2.02.0400.01920.21100.00830.06151.13335.46920.01950.4858
800.00790.14140.00290.04280.14691.34800.01480.3384
1200.00700.11270.00270.03400.06320.61120.00440.2674
1600.00360.09330.00120.02860.03390.16860.00650.2219
3200.00210.06670.00080.02050.01710.11220.00740.1594
1.75400.02390.18010.00660.05236.447011.3890.01020.4745
800.01440.12760.00350.03682.48067.47040.00010.3290
1200.00840.10120.00180.02881.34705.69210.00770.2633
1600.00390.08750.00110.02560.70144.08470.00790.2326
3200.00260.05890.00110.01800.11610.64280.00410.1546
3.0400.01710.16970.00420.050611.852212.43020.02240.4768
800.01510.11730.00280.03417.855912.67500.00850.3171
1200.01230.09480.00310.02766.753411.93300.01550.2591
1600.00920.08600.00150.02465.301411.22550.00520.2308
3200.00530.05980.00120.01722.28747.90280.00730.1593
5.01.0400.01350.47980.00230.12243.35128.26190.85764.6909
800.00700.27750.00290.07430.98874.80770.21701.7650
1200.00350.20310.00360.05560.34422.80670.11570.9897
1600.00270.16460.00290.04590.11341.18070.07590.7488
3200.00270.11030.00140.03080.02700.16280.03980.4915
1.75400.07400.60270.01080.12979.372610.968012.9476.5324
800.01180.30070.00030.07126.869410.67970.32102.4019
1200.00610.20770.00180.05225.013210.31830.12111.3300
1600.00510.16090.00210.04163.51248.95490.05460.8401
3200.00240.10440.00210.02780.96934.85860.01510.5193
3.0400.16270.68330.02800.140010.861312.68771.98768.3799
800.02190.32550.00190.072010.808111.42700.41403.1350
1200.01120.20300.00040.04959.001210.96690.10591.2182
1600.00600.16180.00210.04038.51419.85830.05730.9396
3200.00090.10390.00230.02666.14168.09370.01970.5351
Table 2. Statistic summary for the otis IQ scores data.
Table 2. Statistic summary for the otis IQ scores data.
MeanVarianceSkewnessKurtosis
106.6538.3090.3642.337
Table 3. Parameter estimates (SE) for the fitted models to the otis IQ scores data.
Table 3. Parameter estimates (SE) for the fitted models to the otis IQ scores data.
EstimateBSAPBSNABSNETN
ξ ^ 116.98 (1.010)108.210 (0.946)110.442 (0.459)110.812 (1.336)
η ^ 2.764 (0.338)4.160 (0.280)3.625 (0.187)8.48 (0.985)
β ^ −0.326 (0.107)0.367 (0.076)1.225 (0.611)1.482 (1.025)
α ^ 0.120 (0.040) −0.828 (0.177)−0.592 (0.280)
AIC366.97368.99400.74370.92
CAIC367.82369.49401.59371.77
BIC374.78374.84408.55378.72
HQIC369.96371.23403.73373.91
Table 4. Statistic summary for the the old faithful geyser data.
Table 4. Statistic summary for the the old faithful geyser data.
MeanVarianceSkewnessKurtosis
70.89713.594−0.4141.843
Table 5. Parameters estimates(SE) for the fitted models to the old faithful gayser data.
Table 5. Parameters estimates(SE) for the fitted models to the old faithful gayser data.
EstimateBSAPBSNASNETN
ξ ^ 61.969 (0.850)68.052 (0.376)67.228 (0.242)66.899 (0.916)
η ^ 6.778 (0.231)5.836 (0.122)8.109 (0.204)13.036 (0.589)
β ^ 0.672 (0.064)−0.679 (0.066) 1.547 (0.336)
α ^ 2.637 (0.267) 25.111 (8.807)0.337 (0.107)
AIC2092.422099.312127.042156.88
CAIC2092.562099.392127.122157.02
BIC2106.852110.122137.862171.31
HQIC2098.212103.652131.382162.67

Share and Cite

MDPI and ACS Style

Martínez-Flórez, G.; Tovar-Falón, R.; Jimémez-Narváez, M. Likelihood-Based Inference for the Asymmetric Beta-Skew Alpha-Power Distribution. Symmetry 2020, 12, 613. https://doi.org/10.3390/sym12040613

AMA Style

Martínez-Flórez G, Tovar-Falón R, Jimémez-Narváez M. Likelihood-Based Inference for the Asymmetric Beta-Skew Alpha-Power Distribution. Symmetry. 2020; 12(4):613. https://doi.org/10.3390/sym12040613

Chicago/Turabian Style

Martínez-Flórez, Guillermo, Roger Tovar-Falón, and Marvin Jimémez-Narváez. 2020. "Likelihood-Based Inference for the Asymmetric Beta-Skew Alpha-Power Distribution" Symmetry 12, no. 4: 613. https://doi.org/10.3390/sym12040613

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop