title: The polar-generalized normal distribution
authors: Masoud Faridi Majid Jafari Khaledi
text: 
  
  
  
   The polar-generalized normal distribution
      
   Masoud  Faridi,
  Majid Jafari Khaledi^a
  
  masoudfaridi@modares.ac.ir
  
  ^a Department of Statistics, Tarbiat Modares University, P.O. Box 14115-134, Tehran, Iran
  
  jafari-m@modares.ac.ir
  
      
  ============================================================================================================================================================================================
  
  
  This paper introduces  an extension to the normal distribution  through the polar method to capture bimodality and asymmetry, which are often observed characteristics  of empirical data. The later two features are entirely controlled by a separate scalar parameter. Explicit expressions for the cumulative distribution function, the density function and the moments were derived. The stochastic representation of the distribution facilitates implementing Bayesian estimation via the Markov chain Monte Carlo methods. Some real-life  data as well as simulated data are analyzed  to illustrate the flexibility of the distribution for modeling  asymmetric bimodality.
  
  
   Keywords and phrases: Skewness; Bimodality; Polar method; Beta distribution; Bayesian estimation.
  
  
  § INTRODUCTION
  
  
  The normal distribution is often used  to model data with symmetric distributions.  However, the distribution of data were collected from a wide range of applications reveals asymmetry. Examples arise, for instance, in wind speed data (Zhu and Genton, 2012), temperature data (North et al., 2011), soil and air pollution data (Jafari Khaledi and Rivaz, 2009; Zareifard and Jafari Khaledi; 2013), precipitation data (Xu and Genton, 2017) and diabetes data (Kristensen et al., 2010).
  The skew-normal (SN) distribution, introduced by Azzalini (1985), has been used successfully to model data sets that have asymmetric behavior. The normal distribution is a special case of this distribution. Despite attractive features of this distribution, it is not suitable when the distribution of the data are biomodal. Since many data sets from applied fields are bimodal (for example gene expression patterns in breast cancer (Wang et al., 2009; Ertel, 2010), material characterization in engineering (Dierickx et al.,  2000), water vapor and temperature in meteorology (Zhang et al. , 2003; Dahdouh and Jafari Khaledi, 2020),  there is extensive research on modeling bimodality. The mixture of normal distributions is often used in many different  areas to model data with biomodal distributions. This strategy is often criticized for its non-identifiability problem (McLachlan and  Peel, 2000; Marin, Mengersen and  Robert, 2005). Gomez et al. (2011) introduced a bimodal extension of the skew-normal distribution with application to a pollen dissemination data set, but it suffers complexity and estimation problems inherited from the skewed normal model.
  
  The main objective of this paper is to provide a parsimonious probability distribution that could be used as an alternative to the mixture models and can  account for both asymmetry and bimodality. The basic idea relies on extending the standard polar method.
  In this method, one generate a random pair (X,Y)=(R cos(2π U) , R sin (2π U)) where
  R and U are random polar coordinates. The random variables R and U are independent.  The random variable U is uniformly distributed on the unit interval  (0,1),
   and R has a given distribution that is easy to sample from.  In the context of one-liners, we thus have
  
   
      X    =    R cos( 2π U).
  
   In this case, if R ∼√(χ^2_2), χ^2_2 has a chi-squared distribution with 2 degrees of freedom,
     then the distribution of X is normal. The well-known Box-Muller formula for normal random variables, X = √(-2 log U_2)cos( 2π U_1)
   is thus a simple one-linear, where U_1 and U_2 are independent random variables from the uniform distribution on the unit interval (0, 1).
  
  In our construction, instead of the uniform distribution, we propose to consider the beta distribution for the random variable U in (<ref>). We call the resulting
  construction the polar-generalized normal
  (PGN) distribution. We refer to Alexander et al. (2012),  Mameli and  Musio (2013 and 2016), Alleva and  Giommi (2016) and Eugene, Lee and  Famoye (2002) for information on other generalizations of normal distribution that present bimodality and involve the beta distribution. The PGN distribution can model asymmetry and bimodality and include as special cases the normal distribution. In particular, it contains easily interpretable and estimable parameters. Our approach allows for testing normality and symmetry within the PGN family via Beysian framework.
  
  The remainder of this paper is organized as follows: Section 2 is devoted to  the extension of the polar method. Moreover, basic properties and other characteristics of the PGN distribution are presented.  In particular, some expressions are derived which simplify moments derivation.
  In Section 3, the stochastic representation for a random variable following the PGN distribution is used  to implement Bayesian estimation. The results of a simulation study reported in Section 4 reveal satisfactory behaviour of the Bayesian estimates.  In Section 5, we provide some empirical applications. Finally, this paper concludes with a brief review of the main results.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § THE PGN MODEL  
  
  Below we consider an  extension of the polar method given in equation (<ref>).
  Let V, U  be independent random variables.
  Assume that V has a chi-squared distribution with  k
   degrees of freedom, that is, the density is given by
   f_V(v) =1/2^k/2Γ(k/2) v^k/2-1 e^-v/2 for 0<v< +∞ where Γ(· ) is the gamma function and k>0. Suppose U has a beta distribution with parameters 2μ and 2(1-μ), that is, the density is given by
   f_U(u)=1/B(2μ,2(1-μ)) u^2μ-1 (1-u)^2(1-μ)-1 for 0<u<1 where 0<μ<1 and
   B(α ,β) =Γ(α)Γ(β)/Γ(α +β) is the beta function. We define the polar-generalized normal distribution as the distribution of the product
  
      X=√(V)cos(π U).
  
   We will use the short notation PGN(μ , k) distribution.
    The standard normal distribution is a special case of the PGN(μ , k)  distribution when μ is 0.5 and k  is 2.
  In the following subsections, we study distributional properties and
  present some expansions for
  moments of  PGN(μ , k).
  
  
  
  
  
  
   §.§ The distributional properties
  
  
  
  
  
  
  The cumulative distribution function (cdf) of PGN(μ , k) is given by
  
      F_X(x)   =   {[ I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))
      -1/Γ(k/2)∫_0.5^1γ(k/2 , x^2/2cos^2(π u) ) f_U(u) du                                                                              if   x<0,;   I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))+1/Γ(k/2)∫_0^0.5γ(k/2 , x^2/2cos^2(π u) )f_U(u) du                                                                            if    x≥ 0, ].
  
  where I_c(α , β ) is the regularized incomplete beta function,
  γ(s ,t) is the
  regularized incomplete gamma function and f_U(u) is the density of
  beta distribution with parameters 2μ and 2(1-μ). Specifically,  I_c(α , β )=∫_0^c u^α -1(1-u)^β -1 du and
  γ(s ,t)=∫_0^t y^s-1 e^-y dy.
  
  
  
  
  
      F_X(x)   =   P(X≤ x)=P(√(V)cos(π U)≤ x)=
      ∫_0^1 P(√(V)cos(π u)≤ x)f_U(u) du
       
         =   ∫_0^0.5 P(√(V)cos(π u)≤ x)f_U(u) du +
      ∫_0.5^1 P(√(V)cos(π u)≤ x)f_U(u) du.
  
  
  
  The first term of the Equation (<ref>) is equal to
  
  
      ∫_0^0.5 P(√(V)cos(π u)≤ x)f_U(u) du    =   {[                                    0                            if    x<0; ∫_0^0.5 P(√(V)≤x/cos(π u) )f_U(u) du                           if    x≥ 0 ].
      
         =   {[                                                 0                                        if    x<0,; 1/Γ(k/2)∫_0^0.5γ(k/2 , x^2/2cos^2(π u) )f_U(u) du                                       if    x≥ 0, ].
  
  
  and the second term of the Equation (<ref>) is equal to
  
  
      ∫_0.5^1 P(√(V)cos(π u)≤ x)f_U(u) du    =   {[ ∫_0.5^1 P(√(V)≥x/cos(π u) )f_U(u) du                            if    x<0;                   ∫_0.5^1  f_U(u) du                           if    x≥ 0 ].
      
         =   {[ ∫_0.5^1 f(u) du
      -∫_0.5^1 P(√(V) < x/cos(π u) )f_U(u) du                                               if    x<0;                                      ∫_0.5^1  f_U(u) du                                              if    x≥ 0 ].
      
         =   {[ I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))
      -1/Γ(k/2)∫_0.5^1γ(k/2 , x^2/2cos^2(π u) ) f_U(u) du                                                                             if    x<0,;                                                     I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))                                                                            if    x≥ 0. ].
  
  
  
  Then, the cdf can be computed as follows:
  
  
      F_X(x)   =   {[ I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))
      -1/Γ(k/2)∫_0.5^1γ(k/2 , x^2/2cos^2(π u) ) f_U(u) du                                                                              if   x<0,;   I_0.5(2(1-μ) , 2μ )/B(2μ , 2(1-μ))+1/Γ(k/2)∫_0^0.5γ(k/2 , x^2/2cos^2(π u) )f_U(u) du                                                                            if    x≥ 0. ].    ▪
  
  
  
  
  
  
  
  The probability density function (pdf) of PGN(μ , k) has a two-peice form as
  
      f_X(x)= (x^2)^k-1/2/Γ(k/2) 2^k/2-1×{[ ∫_0.5^11/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du                                               if    x<0,; ∫_0^0.51/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du                                              if    x≥ 0. ].
  
  
  
  
  The proof of this result follows directly by the derivation.
  By the derivation of F_X(x), f_X(x) can be computed as follows:
  
  
      f_X(x)   =   {[ -1/Γ(k/2)∫_0.5^1d/dxγ(k/2 , x^2/2cos^2(π u) ) f_U(u) du                                                     x<0;   1/Γ(k/2)∫_0^0.5d/dxγ(k/2 , x^2/2cos^2(π u) )f_U(u) du                                                    x≥ 0 ].  
         =   {[ |x| (x^2)^k/2-1/Γ(k/2) 2^k/2-1∫_0.5^11/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du                                                                                   x<0,; |x| (x^2)^k/2-1/Γ(k/2) 2^k/2-1∫_0^0.51/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du                                                                                  x≥ 0. ].   ▪
  
  
  
  
  
  
  
  
  
  
  
  
  We could equivalently write the density of PGN(μ , k) as follows:
  
  
  
      f_X(x)    =    (x^2)^k-1/2/Γ(k/2) 2^k/2-1×∫_0.5^11/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u))(f_U(u)1_{ x<0 } +f_U(1-u)1_{ x≥ 0 }) du
       
         =    (x^2)^k-1/2/Γ(k/2) 2^k/2-1×∫_0^0.51/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u))(f_U(1-u)1_{ x<0 } +f_U(u)1_{ x≥ 0 }) du.
  
  
  
  
  
  
  
  
  
  
  
  
  
  By changing  variable v to 1-u and knowing that cos(π (1-v))=-cos(π v), we obtain
  
  
      ∫_0.5^11/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du
         =   ∫_0^0.51/((cos(π (1-v)))^2)^k/2
      e^-(x^2/2(cos(π (1-v)))^2) f_U(1-v) dv
      
         =   ∫_0^0.51/(cos^2(π v))^k/2
      e^-(x^2/2cos^2(π v))  f_U(1-v) dv.
  
  In a similar way we have
  
      ∫_0^0.51/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u)) f_U(u) du
         =   ∫_0.5^11/(cos^2(π v))^k/2
      e^-(x^2/2cos^2(π v))  f_U(1-v) dv.   ▪
  
  
  
  
  
  
  
  
  
  
  
  The probability density function of PGN(μ , k) distribution for different choices of the parameters μ and k are plotted in Figures 1 and 2.
  
  Figure <ref> shows some specific densities of the distribution  PGN(μ , k) and illustrates  the flexibility of this family.
  The distribution may be
  symmetric, skewed to the left, or skewed to the right. At the same time, it can be unimodal or
  bimodal.
  It is noteworthy that the densities can display quite
  different shapes depending on the values of the two parameters. In particular, it can
  be symmetric when μ=0.5   or asymmetric when μ≠ 0.5.
  For μ<0.5, it is easy that to see P(U<1/2)>1/2, or equivalently,
   P(0<U<1/2)>P(1/2<U<1). Accordingly,   P(cos (π U)>0)>P(cos (π U)<0)
  and hence  P(X>0)>P(X<0). So, μ controls the symmetry of the distribution. It confirmed
  by the plots of PGN(μ , k) in Figures <ref>, <ref>, <ref>, <ref>, <ref> and <ref>. As we expected,
  it is observed that for μ <0.5 the distribution is left skewed and for μ >0.5 is right skewed. The distribution also
  approaches symmetric  when μ tends to 0.5.  Asymmetry increases as μ tends to zero or one. So we call μ  the asymmetric  parameter because it regulates the asymmetry of the density function.
  
  We call k the peak parameter because by increasing k , the distribution  tends to have two peaks. These
  peaks are distinguished for μ close to 0.5, but as μ tends to 0 or 1 ,the small peak tends to smooth out.
  From the plots of PGN(μ , k) in Figures <ref>, <ref> and <ref>
  it is observed that for k <2 and μ =0.5  the distribution has a symmetric shape and the distribution tails  are shorter than the tails of a Normal distribution. It is also observed that by reducing k the distribution  tends to has a sharp peak. For k<2 and μ≠ 0.5 the distribution is skewed and by reducing k the distribution  tends to has a sharp peak.
  
  
  Figure <ref> shows under μ=0.5, the distribution is symmetric and for k = 2,  it is also a unimodal distribution (it has one peak). However,
  by increasing k, the distribution
  bifurcates, becomes bimodal and the distance between peaks increase.
  This happens since as k  increases, the mode of the distribution of the random variable √(V) increases, and consequently,  the   distance between
  the mode and the zero point increases. Accordingly, by multiplying cos (π U), which is between -1 and 1, to √(V),  two distinct modes are created (at both sides of zero point) and the distance between two modes increases.
  
  The PGN distribution becomes
  bimodal for certain values of the parameters μ and k, and the analytical solution of μ and k ,
  where the distribution becomes bimodal,cannot
  be solved algebraically. However,  by observing density function of PGN(μ,k) for different values of the parameters μ and k,  we obtain Table (<ref>). The table
   shows a grid of values where the distribution
  is bimodal.
  In the table, “1” indicates that PGN distribution is unimodal and “2” indicates that PGN distribution is bimodal.
  
  
  
  
  
  
  
  
  If μ=0.5, the cdf of the  PGN(μ,k) model  is  given by
  
      f_X(x)=
       (x^2)^k-1/2e^-x^2/2/√(π)Γ(k/2) 2^k/2
        U(1/2 , k+1/2,x^2/2),
  
  where U(a,b,z) is Tricomi's (confluent hypergeometric) function, that is,
  
  
      U(a,b,z)=
      1/Γ (a)∫_0^∞ e^-zu u^a-1 (1+u)^b-a-1 du.
  
  
  
  
  
  
  By Equation <ref> and knowing that for μ=0.5, f_U(u)=1, we have
  
  
  
      f_X(x)   =   |x| (x^2)^k/2-1/Γ(k/2) 2^k/2-1∫_0^0.51/(cos^2(π u))^k/2
      e^-(x^2/2cos^2(π u))  du
      
         =    (x^2)^k-1/2e^-x^2/2/Γ(k/2) 2^k/2-1∫_0^0.5 (1+ tan^2(π u))^k/2
      e^-x^2/2( tan^2(π u))  du
      
         =    (x^2)^k-1/2e^-x^2/2/Γ(k/2) 2^k/2-1∫_0^+∞ (1+ t)^k/2
      e^-x^2/2 t1/π 2(1+t)√(t)dt
      
         =    (x^2)^k-1/2e^-x^2/2/Γ(k/2) 2^k/2-11/2π∫_0^+∞ e^-x^2/2 t t^1/2-1(1+ t)^k/2-1
        dt
        
         =    (x^2)^k-1/2e^-x^2/2/Γ(k/2) 2^k/21/π∫_0^+∞ e^-x^2/2 t t^1/2-1(1+ t)^k+1/2-1/2-1
        dt
        
         =    (x^2)^k-1/2e^-x^2/2/Γ(k/2) 2^k/21/√(π)  U(1/2 , k+1/2,x^2/2).   ▪
  
  
  
  
  In the rest of this section, we provide some
  distributional properties of the PGN(μ,k)
  distribution.
  
  
  
  If μ=0.5, the density of the PGN(μ,k) model is symmetric about zero.
  
  
   The proof is simple since the cdf <ref> satisfy  f_X(-x)=f_X(x). ▪
  
  
  
  
  
  
  Table <ref> shows some special cases of the  Tricomi confluent hypergeometric function (Wolfram Research, 2020).
  
  
  
   The standard normal distribution is a special case of the PGN(μ , k) distribution.
  
  
  
   For the special case μ=0.5 and  k=2,
  U(1/2,3/2,z)=1/√(z), we obtain the standard normal distribution since
  f_X(x)=(x^2)^2-1/2e^-x^2/2/√(π)Γ(2/2) 2^2/21/√(x^2/2)=e^-x^2/2/√(2π). ▪
  
  
  
  
  
  
  
  
  
  For the special case μ=0.5 and  k=4, the PGN(μ , k) distribution  is an equal mixture of two different families, normal distribution and Chi distribution with 3 degree of freedom.
  
  
  
  
  
  
  According to U(1/2,5/2,z)=2z+1/2z^3/2,
  we have
  
      f_X(x)    =   (x^2)^4-1/2e^-x^2/2/√(π)Γ(4/2) 2^4/2
        U(1/2 , 5/2,x^2/2) 
         =   (x^2)^3/2e^-x^2/2/√(π)Γ(2) 2^22x^2/2+1/2(x^2/2)^3/2
         =   (x^2+1)e^-x^2/2/2√(2π)
         =1/2x^2e^-x^2/2/√(2π)
         +1/2e^-x^2/2/√(2π),
  
  which is an equal mixture of two different families, normal distribution and Chi distribution with 3 degree of freedom. ▪
  
  
  With regard to Corollary <ref>,  the distribution is bimodal, as you will see in Figure(<ref>). From the table
  we can see for k≥ 3 the pdf  can be written as a mixture of two pdf and this result could be seen in Figure(<ref>).
  Thus, the PGN distribution provides great flexibility in modeling symmetric,
  able to accommodate both unimodal and bimodal cases.
  
  
  For μ→ 0, the PGN(μ,K) distribution becomes χ_k, the Chi distribution with
  k degrees of freedom, that is, the density is given by
  f_X(x)=x^k-1e^-x^2/2/2^k/2-1Γ(k/2) for x>0. For μ→ 1 the PGN(μ,K) distribution becomes -χ_k.
  
  
  
  For μ→ 0, the beta distribution becomes a  Degenerate distribution at point 0 and for μ→ 1, the beta distribution becomes a  Degenerate distribution at point 1. ▪
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Let X∼ PGN(μ , k), the fraction P(X>0)/P(X<0) only depends on μ, we  denote
   φ (μ)= P(X>0)/P(X<0). We have
  
  
      φ (μ )>1     ⟺   μ < 0.5, 
      φ (μ )< 1     ⟺   μ > 0.5, 
      φ (μ )=1     ⟺   μ = 0.5.
  
  
  
  
  
      φ (μ )   =   P(X>0)/P(X<0)=P(√(V)cos(π U)>0)/P(√(V)cos(π U)<0)
         =   P(U<0.5)/P(U>0.5)=1/B(2μ , 2(1-μ))∫_0^0.5 u^2μ-1(1-u)^1-2μ du/1/B(2μ , 2(1-μ))∫_0.5^1 u^2μ-1(1-u)^1-2μ du
         =   ∫_0^0.5 (u/1-u)^2μ-1 du/∫_0.5^1 (u/1-u)^2μ-1 du = ∫_0^0.5 (u/1-u)^2μ-1 du/∫_0^0.5 (u/1-u)^1-2μ du
         =   ∫_0^0.5 h(u)^2μ-1 du/∫_0^0.5 h(u)^1-2μ du,
  
  with h(u)=u/1-u it is easy to check that 0<h(u)<1 for every u ∈ (0,0.5). Hence, if 2μ -1<0, i.e. if μ< 0.5,
  we have that h(u)^2μ -1 > h(u)^1-2μ, ∀ u ∈ (0,0.5), so that φ(μ) >1. On the other hand, the inequality is reversed if μ > 0.5 and, finally, φ(μ) =1
  if μ=0.5. ▪
  
  
  
  
  
  
  Figure <ref> shows φ(μ) for μ∈ (0,1). Accordingly  we came to the conclusion that μ controls the
  asymmetric of the distribution .
  On the other hand, μ controls the allocation of mass to each side of the point 0.
  To be more specific, the asymmetry of the distribution  is only controlled by μ.
  
  
  
   The  PGN(μ , k) distribution can not be symmetric around any a≠ 0.
   
  
   
   To see this, let the  PGN(μ , k) distribution be symmetric around a. Let a>0, we have  ∀ x>0    P(X<a-x) = P(X>a+x).
   For x=a we obtain
  
   
      P(X<0) = P(X>2a).
  
  
   The left hand side of <ref> only depend of μ, so by changing k it does not change, But as changing
   k, the  right hand side of <ref> will change. So the only way to the  PGN(μ , k) distribution
   be symmetric is for a=0. The proof of this for a<0 is similar to the proof for
   a>0. ▪
   
  
  
  From the above Proposition and Proposition (<ref>), we deduce that the PGN(μ, k) is symmetric around 0, if and only if μ=0.5
  
  
  
  
  
  
  
  
  
  
  
  
  Let X∼ PGN(μ , k), then -X∼ PGN((1-μ) , k).
  
  
  Knowing that
  1-U ∼ Beta(2(1-μ),2μ),
  the proof of this result follows directly from definition (<ref>) as follows:
  
      -X=-√(V)cos(π U)=√(V)cos(π -π U)=√(V)cos(π(1- U))∼  PGN(1-μ,k).   ▪
  
  
  
  From the above proof we can deduce that if the distribution of X is bimodal, the distribution of -X is bimodal too. We can
  see this in Table (<ref>).
  By the table the region of bimodality of the PGN distribution for μ and 1-μ is same. Since by the above lemma,
  if X∼ PGN(μ , k), -X∼ PGN((1-μ) , k) and k,that controls  bimodality of the distribution, does not change.
  
  
  
  Let V∼χ_k^2 and U ∼ Beta (2μ , 2(1-μ)) the following properties are deduced immediately from
  the definition:
  
  
  
   
      √(V)sin(π ( U-1/2))    ∼     PGN(μ ,k). 
      √(V)sin(π (1/2- U))    ∼     PGN(1-μ , k).
      √(V)( 2cos^2(π U/2) -1   )		   ∼     PGN(μ ,k).
      √(V)(1- 2sin^2(π U/2)   )		   ∼     PGN(μ ,k).
      √(V)( cos^2(π U/2) -sin^2(π U/2)   )		   ∼      PGN(μ ,k).
      √(V)1/2( e^iπ U+e^-iπ U)		   ∼     PGN(μ ,k).
  
  
  
  
  
  
  
  
  The proof is obvious from both Definition (<ref>)  and
   
      cos(π U)   =   1/2( e^iπ U+e^-iπ U)
      
         =   
       2cos^2(π U/2) -1
       
         =   
      1- 2sin^2(π U/2)
       
         =   cos^2(π U/2) -sin^2(π U/2).   ▪
  
  
  
  
  
  
  
  
  
  
  
   §.§ Moments derivation
  
  
  
  In this section, some results on the
  moment properties of PGN(μ ,k) distribution
  are obtained.
  
  
  
  Let U has a beta distribution with parameters
  a and b. Define  I_(m)(a,b)=E(cos^m( π U)), so
  
  
  
  
      I_(1)(a , b )   =   1/2( _1F_1(a;a+b;iπ) +  _1F_1(a;a+b;-iπ) ), 
      
      I_(2)(a , b )   =   1/4(2+
       _1F_1(a;a+b;i2π) +  _1F_1(a;a+b;-i2π)
       ),
       
      
      I_(3)(a , b )   =   3/8(
        _1F_1(a;a+b;iπ) +  _1F_1(a;a+b;-iπ)  ) +
       1/8(
       _1F_1(a;a+b;i3π) +  _1F_1(a;a+b;-i3π) ),
        
      
      I_(4)(a , b )   =   1/16(6+4
       _1F_1(a;a+b;i2π) + 4 _1F_1(a;a+b;-i2π)
       +
      _1F_1(a;a+b;i4π) +  _1F_1(a;a+b;-i4π)
      ),
  
  
   where
    
      (a)_(n)=Γ (a+n)/Γ(a)=∏_k=0^n-1 (a+k),
  
   and
   _1F_1(a;b;z) is the Kummer's (confluent hypergeometric) functions, that is,
  
   
      _1F_1(a;b;z)=∑_n=0^∞(a)_(n)/(b)_(n)z^n/n!
       = ∫_0^11/B(a,b) e^zu u^a-1 (1-u)^b-a-1 du.
  
  
  
  
  
  
  
  
  
  Define   J_(m)(a,b)=E(cos(m π U)) so
   
      J_(m)(a,b)   =   1/B(a,b)∫_0^1cos (mπ u) u^a-1 (1-u)^b-1  du
      
         =   1/B(a,b)∫_0^11/2(e^i m π u+e^-i m π u) u^a-1 (1-u)^b-1  du
      
         =   1/2 ( CF(i m π) +CF(i m π))
      
         =   1/2 ( _1F_1(a;a+b;i m π)+_1F_1(a;a+b;-i m π)),
  
   where CF(· ) is the characteristic function of the beta distribution. According to
  
  
  
  
  
  
      cos^2(x)   =   1/2(1+cos(2x)),
      cos^3(x)   =   1/4 (3cos(x)+cos(3x)), 
      cos^4(x)   =   1/8 (3+4cos(2x)+cos(4x)),
  
  
  we conclude that
  
  
      I_(1)(a , b )   =   J_(1)(a , b )
      = 1/2( _1F_1(a;a+b;iπ) +  _1F_1(a;a+b;-iπ) ), 
      
      I_(2)(a , b )   =   1/2(1+ J_(2)(a , b ) )
       
         =   1/2(1+
       1/2_1F_1(a;a+b;i2π) +  1/2_1F_1(a;a+b;-i2π)
       )
        
         =   1/4(2+
       _1F_1(a;a+b;i2π) +  _1F_1(a;a+b;-i2π)
       ),
       
      
      I_(3)(a , b )   =   3/4 J_(1)(a , b )+1/4
       J_(3)(a , b )
       
         =   3/8(
        _1F_1(a;a+b;iπ) +  _1F_1(a;a+b;-iπ)  ) +
       1/8(
       _1F_1(a;a+b;i3π) +  _1F_1(a;a+b;-i3π) ),
        
      
      I_(4)(a , b )   =   1/8(3+4 J_(2)(a , b )+ J_(4)(a , b)  )  
         =   1/8(3+4
       1/2{_1F_1(a;a+b;i2π) +  _1F_1(a;a+b;-i2π) }
      +
       1/2{_1F_1(a;a+b;i4π) +  _1F_1(a;a+b;-i4π) })
       
         =   1/16(6+4
       _1F_1(a;a+b;i2π) + 4 _1F_1(a;a+b;-i2π)
       +
      _1F_1(a;a+b;i4π) +  _1F_1(a;a+b;-i4π)
      ).   ▪
  
  
  
  
  We will now determine the moments of PGN(μ,k) using Lemma <ref>.
  
  
  Let X∼ PGN(μ,k) , we can deduce the following properties:
  
  
      E(X)   =   Γ(k+1/2)/√(2)Γ(k/2)( _1F_1(2μ ;2;iπ) +  _1F_1(2μ ;2; -iπ) ),  
      
         E(X^2)   =   k/4(2+ _1F_1(2μ ;2;i2π) +  _1F_1(2μ ;2;-i2π) ), 
      
        E(X^3)    =    (k +1)   √(2)Γ(k+1/2)/8Γ(k/2){3   _1F_1(2μ;2;iπ) + 3 _1F_1(2μ;2;-iπ)  )
       +   _1F_1(2μ ;2;i3π) +  _1F_1(2μ ;2;-i3π)  }, 
      
        E(X^4)    =    k (k +2)/16(6+4
       _1F_1(2μ ;2;i2π) + 4 _1F_1(2μ ;2;-i2π)
      +
      _1F_1(2μ ;2;i4π) +  _1F_1(2μ ;2;-i4π)
      ).
  
  
  
  
  
  
  Since
   V and U are independent, we have
  
      E(X^m)=E(√(V)^m) E(cos^m(U))=Γ(k+m/2)/Γ(k/2)2^m/2 E(cos^m(π U))=Γ(k+m/2)/Γ(k/2)2^m/2 I_m(2μ,2(1-μ)).
  
  
  
  
  
  
  
      E(X)   =   Γ(k+1/2)/Γ(k/2)√(2)1/2( _1F_1(2μ ;2;iπ) +  _1F_1(2μ ;2; -iπ) )
      
         =   Γ(k+1/2)/√(2)Γ(k/2)( _1F_1(2μ ;2;iπ) +  _1F_1(2μ ;2; -iπ) ).
  
  
  
      E(X^2)   =   E(√(V)^2) E(cos^2(π U))= Γ(k+2/2)/Γ(k/2)2^2/2 E(cos^2(π U))
       
         =   
      k I_(2)(2μ,2(1-μ))
       
         =   
      k  1/4(2+ _1F_1(μ ;1;i2π) +  _1F_1(μ;2;-i2π) )
       
         =   k/4(2+ _1F_1(2μ ;2;i2π) +  _1F_1(2μ ;2;-i2π) ).
  
  
  
  
      E(X^3)   =E(√(V)^3) E(cos^3(π U))=  Γ(k+3/2)/Γ(k/2)2^3/2 E(cos^k(π U))
      =
      (k +1)   Γ(k+1/2)/Γ(k/2)√(2)  I_(3)(2μ,2(1-μ))
      
         =
      (k +1)   Γ(k+1/2)/Γ(k/2)√(2){3/8(
        _1F_1(2μ ; 2;iπ) +  _1F_1(2μ ; 2;-iπ)  )
       +  1/8(
       _1F_1(2μ ; 2;i3π) +  _1F_1(2μ ; 2;-i3π) ) }
         =
      (k +1)   Γ(k+1/2)/Γ(k/2)√(2)/8{3   _1F_1(2μ ; 2;iπ) + 3 _1F_1(2μ ; 2;-iπ)  )
       +   _1F_1(2μ ; 2;i3π) +  _1F_1(2μ ; 2;-i3π)  }
         =
      (k +1)   Γ(k+1/2)/Γ(k/2)√(2)/8{3   _1F_1(2μ ;2;iπ) + 3 _1F_1(2μ ;2;-iπ)  )
       +   _1F_1(2μ ;2;i3π) +  _1F_1(2μ ;2;-i3π)  }
         =
       (k +1)   √(2)Γ(k+1/2)/8Γ(k/2){3   _1F_1(2μ ;2;iπ) + 3 _1F_1(2μ ;2;-iπ)  )
      +    _1F_1(2μ ;2;i3π) +  _1F_1(2μ ;2;-i3π)  }.
  
    
  
  
  
      E(X^4)   =   E(√(V)^4) E(cos^4(π U))=
      Γ(k+4/2)/Γ(k/2)2^4/2
       E(cos^4(π U))
      =
      k (k +2)  I_(4)(2μ,2(1-μ))
       
         =   
       k (k +2)
      1/16(6+4
       _1F_1(2 μ ; 2;i2π) + 4 _1F_1(2 μ ; 2;-i2π)
      +
      _1F_1(2 μ ; 2;i4π) +  _1F_1(2 μ ; 2;-i4π)
      )
      
         =    k (k +2)/16(6+4
       _1F_1(2 μ ; 2;i2π) + 4 _1F_1(2 μ ; 2;-i2π)
      +
      _1F_1(2 μ ; 2;i4π) +  _1F_1(2 μ ; 2;-i4π)
      ).   ▪
  
  
  
  
  
  
  
   The variance of the PGN(μ, k) is
  
  
      Var(X)   =   k/4(2+ _1F_1(2μ ;2;i2π) +  _1F_1(2μ ;2;-i2π) )
      -(
      Γ(k+1/2)/√(2)Γ(k/2)( _1F_1(2μ ;2;iπ) +  _1F_1(2μ ;2; -iπ) )
       )^2.
  
  
  
  
  If X∼ PGN(μ , k), the moment of order m is given by
  
  
      E(X^m)   =   Γ(k+m/2)/Γ(k/2)2^m/2 I_m(2μ,2(1-μ))
      
         =   Γ(k+m/2)/Γ(k/2)2^k-2m/2×{[ (  ∑_r=0^2r<m m  2r 
       ( _1F_1(2 μ ; 2;i  (m-2r)π)+_1F_1(2 μ ; 2;-i  (m-2r)π))
        )		+		 m m/2                                                                                    m=2,4,⋯ ,;        (  ∑_r=0^2r<m m  2r 
       ( _1F_1(2 μ ; 2;i  (m-2r)π)+_1F_1(2 μ ; 2;-i  (m-2r)π))
            )                                                                                    m=1,3,⋯ .;                                                                                              ].
  
  
  
  
  
  
  
  
  With regard to
  
  
      cos^m(x)=
      {[ 1/2^m-1(  ∑_r=0^2r<m m  2r cos (m-2r)x        )		+		1/2^mm m/2                                                      m=2,4,⋯ .;                1/2^m-1(  ∑_r=0^2r<m m  2r cos (m-2r)x        )                                                      m=1,3,⋯ ,;                                                                ].
  
  we obtain
  
      I_(m)(a,b)   =   ∫_0^11/B(a,b)cos^m(π u) u^a-1 (1-u)^b-1 du
      
         =   {[ 1/2^m-1(  ∑_r=0^2r<m m  2r ∫_0^11/B(a,b)cos( (m-2r)π u)  u^a-1 (1-u)^b-1 du
        )		+		1/2^mm m/2                                                                                        m=2,4,⋯;            1/2^m-1(  ∑_r=0^2r<m m  2r ∫_0^11/B(a,b)cos( (m-2r)π u)  u^a-1 (1-u)^b-1 du
            )                                                                                        m=1,3,⋯;                                                                                                ].
      
         =   {[ 1/2^m-1(  ∑_r=0^2r<m m  2r 1/2 ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
        )		+		1/2^mk m/2                                                                                               m=2,4,⋯;            1/2^m-1(  ∑_r=0^2r<m m  2r 1/2 ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
            )                                                                                               m=1,3,⋯;                                                                                                       ].
      
         =   {[ 1/2^m(  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
        )		+		1/2^mm m/2                                                                                         m=2,4,⋯ ,;            1/2^m(  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
            )                                                                                         m=1,3,⋯ .;                                                                                                   ].
  
  Thus,
  
  
      E(X^m)   =    E(√(V)^m)  E(cos^m(π U)) = E(√(V)^m)  E(cos^m(π U))
      
         =   Γ(k+m/2)/Γ(k/2)2^k/2 I_(m)(a,b)
      
         =   Γ(k+m/2)/Γ(k/2)2^k/2×{[ 1/2^m(  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
        )		+		1/2^mm m/2                                                                                           m=2,4,⋯;            1/2^m(  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
            )                                                                                           m=1,3,⋯;                                                                                                   ].
      
         =   Γ(k+m/2)/Γ(k/2)2^k-2m/2×{[ (  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
        )		+		 m m/2                                                                                m=2,4,⋯ ,;        (  ∑_r=0^2r<m m  2r 
       ( _1F_1(a;a+b;i  (m-2r)π)+_1F_1(a;a+b;-i  (m-2r)π))
            )                                                                                m=1,3,⋯ ,;                                                                                          ].
  
  
  the result is proved by replacing a=2μ and b=2(1-μ). ▪
  
  
  
  From the Equation (<ref>), E(X^m)  can be easily calculated.
   Table <ref> in Appendix  gives the value of
   I_m(2μ,2(1-μ)) under different values of μ and m.
   The hypergeo package (Hankin, 2015) of R  routines  is used for numerical calculation of
  hypergeometric functions.  For example, under μ=0.3
  and k=5,
  
      E(X)    =   Γ(k+1/2)/Γ(k/2)√(2) I_1(2μ,2(1-μ))= 2.1277 ×  0.46=0.9787, 
      
        E(X^2)   =   Γ(k+2/2)/Γ(k/2)2^2/2 EI_2(2μ,2(1-μ)) = 5× 0.5995=2.9975,
      
        E(X^3)   =   Γ(k+3/2)/Γ(k/2)2^3/2 I_3(2μ,2(1-μ))= 12.7661×  0.3942=5.0324,  
      
        E(X^4)   =   Γ(k+4/2)/Γ(k/2)2^4/2
       I_4(2μ,2(1-μ))=35× 0.4916=17.2060.
  
  
  
  
  
  
  
  § BAYESIAN ANALYSIS
  
  
  Here we specify a Bayesian model that accounts for asymmetry and bimodality. In what follows,
  the posterior inference from this model is discussed.
  Let X_1, X_2, ⋯, X_n be n independent and identically distributed random variables from PGN(μ , k).
  The Bayesian model specification requires a prior distributions for all the unknown
  parameters, i.e., θ=(μ , k).
   In the absence of prior information
  and in order to guarantee posterior property, we adopt proper but diffuse priors. For
  convenience but not always optimal, we suppose that elements of θ are independent
   so the joint prior distribution
  π (θ)=π (μ , k)=π (μ ) π (k ).
  we adopt the following prior distributions
  
      μ   ∼    Uniform(0,1), 
      
         k    ∼   Γ (k_0 , k_1).
  
  Here, the hyper-parameters k_0 and k_1 are known (to be specified) positive scalers. We choose the non-informative uniform prior for parameter μ to allow the data to select appropriate value for μ.
  Then the posterior distribution is proportional
  to
  
      π(θ|𝐱)   ∝   π(𝐱 |θ ) π(θ)
         ∝     (∏_i=1^np(x_i |μ , k))π (μ) π (k),
  
  where 𝐱=(x_1,x_2, ⋯, x_n)  and  p denotes  the density function of PGN(μ , k) given in  (<ref>). The joint posterior distribution of all unknown quantities involved is given by
  
      π(θ, 𝐮|𝐱)   ∝   ∏_i=1^np(x_i|u_i , k)p(u_i|μ)π (μ) π (k),
  
  where 𝐮=(u_1,u_2, ⋯, u_n) and 
      p(x_i|u_i , k)     =    {[ 1/|cos (π u_i)|2^k/2-1Γ(k/2)( x_i/cos (π u_i))^k/2-1  e^- x_i^2/2cos^2 (π u_i)                                                                        u_i>1/2                                                                          x_i<0; 1/|cos (π u_i)|2^k/2-1Γ(k/2)( x_i/cos (π u_i))^k/2-1  e^- x_i^2/2cos^2 (π u_i)                                                                        u_i<1/2                                                                            x>0;                                                                              0                                                                           O.W. ].
       
         =    1/|cos (π u_i)|2^k/2-1Γ(k/2)(| x_i/cos (π u_i)|)^k/2-1  e^- x_i^2/2cos^2 (π u_i) 1_{x_i×(1/2-u_i)}.
  
  The posterior distribution is in a complicated form and we will use a Markov chain Monte Carlo (MCMC) method to generate samples from posterior distributions. To generate samples from the posterior distribution, we exploit the full conditional distributions in a Gibbs sampling framework.
  The full conditional p(u_i|x_i, μ, k) does not defines a standard probability distribution, so sampling from this distribution is not simply practicable. Some methods like rejection sampling, importance sampling or Metropolis-Hastings sampler can be used for sampling from this full conditional. We  use Metropolis-Hastings sampler with generate proposals of beta distribution that the mean of it, centered on the current value. That is, in the Gibbs  sampling algorithm, at step  t, u_i^(t+1) is generated from Beta(2u_i^(t),2(1-u_i^(t))).
  
  The full conditional posterior distribution  of the parameter k, i. e. π (k |  𝐱, 𝐮, μ),  can be computed as follows:
  
      π (k | 𝐱, 𝐮 , μ)    ∝   
      p(x_1|u_1, k)
      p(x_2|u_2, k) ⋯ p(x_n|u_n, k)  π (k) 
         ∝   π (k) ∏_i=1^n1/|cos (π u_i)|^k/22^k/2-1Γ(k/2)  e^- x_i^2/2cos^2 (π u_i)
       (u_i/1-u_i)^2μ-1
        1_{x_i×(1/2-u_i)}
         ∝   π (k  ) ∏_i=1^n1/|2cos (π u_i)|^k/2Γ(k/2)
        1_{x_i×(1/2-u_i)}
         ∝   π (k  )/2^n k/2|∏_i=1^ncos (π u_i)|^k/2Γ^n(k/2)∏_i=1^n1_{x_i×(1/2-u_i)}.
  
  The full conditional of k does not define a standard probability distribution too.
  The Metropolis–Hastings sampler can be used for sampling from this full conditional with generate proposals
  of an exponential distribution, centered on the current value. That is, in the Gibbs  sampling algorithm, at step t, k^(t+1) is generated from Exp(k^(t)).
  
  
  The full conditional posterior distribution  of the parameter k, π (μ | 𝐱, 𝐮, k),  is given by
  
      π (μ | 𝐱, 𝐮 , k)
          ∝   
       p(u_1| μ )
      p(u_2| μ ) ⋯ p(u_n| μ )  π (μ  )
      
         ∝   π (μ  )/Γ^n(2μ) Γ^n(2(1-μ))(∏_i=1^nu_i/1-u_i)^2μ-1.
  
  
  The full conditional
  does not defines a standard probability distribution too. The Metropolis–Hastings sampler can be used for sampling from this full conditional with generate proposals
  of beta distribution, centered on the current value.That is, in the Gibbs  sampling algorithm, at step t, μ^(t+1) is generated from Beta(2μ^(t),2(1-μ^(t))).
  
  
  Based on these full conditional distributions, we use the following algorithm for sampling from the
  joint posterior distribution:
  The main steps of the Gibbs  sampling algorithm
  at step (t + 1) as
  
  
    
    * Draw u^(t+1)_i from p(u_i|x_i,μ^(t),k^(t)) for i=1,2,⋯, n.
    
    * Draw k^(t+1) from π (k|𝐱,𝐮^(t+1),μ^(t)).
    
    * Draw μ^(t+1) from π (μ|𝐱,𝐮^(t+1),k^(t+1)).
    
    * Iterate above steps until we get the appropriate number of MCMC samples.
  
  
  
  
  
  § SIMULATION STUDY
  
  
  
  
  
  
  
  Here we assess the finite sample behaviour of the Bayesian estimators of the parameters in a PGN model using a simulation study. We perform 100 Monte Carlo replications.  In each replication a random sample of size n=50 is drawn from PGN(μ , k). In order to evaluate the accuracy of posterior estimates under different scenarios, we consider various different values for the parameters μ and k.  The true parameter values used in the data generating process are μ=0.1,0.2,0.25,0.3,0.4,0.5,0.6,0.7,0.75,0.8,0.9 and k=1,2,5,10,15. The parameters are estimated from the sample by the Bayesian approach.  For a comparison, we
  computed two criteria: the average of bias (Bias) and the root of mean-square error
  (RMSE). The OpenBUGS is used  for fitting the model. We generate 2,000 MCMC samples after a burn-in 1,000 iterations in one chain with different starting values. We resampled the components every 10 iterations to improve mixing. Visual inspection of the trace plots does not
  reveal any obvious mixing or convergence problems, and autocorrelations seem within reasonable levels.
  
  
  
  The posterior results are summarized in Tables <ref> and  <ref>.
  As seen,   the results indicate overwhelming support for the proposed model.
  Overall, it is observed that the Bias and RMSE   of the parameter estimates are relatively small.
  This illustrates the effectiveness of the estimation procedure outlined in Section 3.
  As expected, the smallest Bias's correspond to   μ=0.5. This is a
  direct consequence of the fact that we have maximum information when the density is symmetric. There is a small bias in the estimation of the parameters when μ is close to 0 or 1.   It must be noted that  the results have been reported based
  on small sample sizes.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § ILLUSTRATIVE EXAMPLE
  
  In this section, we apply the location and scale extension of the standard PGN model to five illustrative data sets. This extension, denoted by PGN(β , σ^2 , μ ,k), is given by the distribution of
  Y = β  + σ X, where X ∼ PGN(μ ,k) and  β∈ R , σ  > 0 are the location and scale parameters, respectively.
  
  The data sets  contain clinical annotations and the log expression of 150 genes for a set of 444 lung
  cancer patients. The  data (lung.dataset) are  available in oompaData package from https://cran.r-project.org/package=oompaDatahttps://cran.r-project.org/.
  The original data are also available in FTP from ftp://caftpd.nci.
  nih.gov/pub/caARRAY/experiments/caArray_beer-00153/here and in the Gene Expression Omnibus at http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68571here. The data were
  log transformed by mapping the expression value x to log 2(1+x). The PGN model is fitted to  the gene ♯ 37145 (row 2), the gene ♯ 202459 (row 12), the gene ♯ 208288 (row 63),  the gene ♯ 215456 (row 106) and the gene ♯ 216437 (row 113). These data sets display bimodality and asymmetry (see Figure 6).
  The parameters are estimated via the Bayesian approach as described in Section 3.
  The priors of the parameters  β and σ^2
  are specified as:
  
      β | σ^2     ∼   Normal(β_0 , cσ^2) ,
      τ =1/σ^2   ∼   Γ(τ_0 ,τ_1 ).
  
  Here, the hyper-parameters c, τ_0 and τ_1   are known (to be specified) positive scalers and β_0  is a real-valued scaler. The hyper-parameter β_0 are set to be close to mean of the data. The hyper-parameters τ_0 and τ_1 are set to be some values such that
  to ensure  large prior variance for the parameter τ.
  If the histogram of the data be unimodal the hyper-parameters k_0 and k_1 are chosen in the sense  that E(k) ≈ 2. For the bimodality cases, these hyperparameters are considered such that E(K)>2.
  
  
  
  
  We generate 2,000 MCMC samples after
  a burn-in 1,000 iterations in one chain with different starting values. We resampled the components every 10 iterations.
  The posterior means and the 95%
   highest probability density (HPD) credible intervals of the model parameters are given in Table <ref>.
    The HPD interval is one of the most useful tools to measure the posterior uncertainty.
  The short HPD intervals for the parameters, in all data set,
  indicate the accuracy of the estimates. Our approach allows for testing symmetry and normality  within the PGN  family via the HPD intervals. The central position μ=0.5 allows
   for testing symmetry because under
  μ=0.5 the PGN family is symmetric.
   In a similar way,  a normality test
   is H_0: μ=0.5 &    k=2. The HPD intervals are used  to make inferences about hypothesis testing of parameters. Specifically, in the first four data sets, the HPD reject the hypothesis
  μ=0.5, so these data sets are not symmetric. However, for the data set 5,
  although we reject the normality assumption but  we  can not reject that the distribution is symmetric.
  
  
  The observed
  and estimated densities plotted in Figure 6. As seen, the PGN model is able to capture different distribution shapes which are observed
  in the  data sets. This result confirms the great flexibility of the model in accommodating bimodality and asymmetry.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSIONS
  
  This article proposed a new class of generalized normal distributions which includes as a special case the normal distribution. Two well known distributions,  i.e.
  the Chi distribution and the half normal distribution, are limiting cases of the general construction. The proposed model provides a flexible approach to account
  for  asymmetry and biomodality which are two pervasive features of the complex data.
  The peak parameter and the asymmetric parameter  control the biamodality and asymmetry of the model, respectively. A Gibbs algorithm was developed for
  the Bayesian inference. A simulation study is performed. The simulation results established that
  the proposed model has an appropriate ability to estimate
  model parameters. Several applications to real data have shown the usefulness of the PGN model for applied statistical research. There is a
  bias in the estimation of the parameters when μ is close 0 or 1, so further research could be conducted to obtain bias corrections
  for these estimators, thus reducing their systematic errors in the estimation of the parameters.
  It is also provide a test for symmetry and a test for normality in the context of Bayesian framework.
  The class was in the context of univariate, so further research could be conducted to extend this class to multivariate contexts. It would be interesting if the properties and the flexibility of
  this extension  were to be explored and compared in a further work.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   Alleva, G., and  Giommi, A. (2016). Topics in theoretical and applied statistics. Cham: Springer
  International Publishing.
  
   Alexander, C., Cordeiro, G.M., Ortega, E.M.M. and Sarabia, J.M. (2012). Generalized beta-generated
  distributions. Computational Statistics and Data Analysis, 56, 1880-189.
  
   Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian Journal of Statistics, 12, 171–178
  
   McLachan, G., and Peel, D. (2000). Mixture Models: Inference and Applications to Clustering. Marcel Dekker, N.Y.
  
   Dahdouh, O., and Jafari Khaledi, M. (2020). Generalized Spatial Stick-Breaking Processes. Communications in Statistics - Simulation and Computation, 1-20.
  
   Dierickx, D., Basu, B., Vleugels, J., and Van der Biest, O. (2000). Statistical extreme value modeling of particle size distributions: experimental grain size distribution type estimation and parameterization of
  sintered zirconia. Mater. Characteriz.,45(1): 61–70.
  
   Ertel, A. (2010). Bimodal gene expression and biomarker discovery.Cancer Inform, 9, 11–14.
  
   Eugene, N., Lee, C., and Famoye, F. (2002). Beta-normal distribution and its application. Commun Stat Theory Methods, 31(4), 497–512.
  
   Gómez, HW., Elal-Olivero, D., Salinas, HS., et al. (2011). Bimodal extension based on the skew-normal
  distribution with application to pollen data. Environmetrics, 22(1), 50–62.
  
  
  Hankin, R. K. (2015). Numerical Evaluation of the Gauss Hypergeometric Function with the hypergeo Package. R J., 7(2), 81.
  
   Khaledi, M.J., and Rivaz, F. (2009). Empirical Bayes spatial prediction using a Monte Carlo EM algorithm. Stat. Methods Appl., 18: 35–47.
  
   Kristensen, P.L., Pedersen-Bjergaard, U., Schalkwijk, C., Olsen, N.V., and Thorsteinsson, B. (2010). Erythropoietin and vascular endothelial growth factor as risk markers for severe hypoglycaemia in type 1
  diabetes.Eur. J. Endocrinol. 163(3): 391–398.
  
  
   Mameli, V., and Musio, M. (2013). A generalization of the skew-normal distribution: the beta skew-normal.
  Commun Stat Theory Methods, 42(12), 2229–2244.
  
   Mameli, V., and Musio, M. (2016). Some new results on the beta skew-normal distribution. In: Alleva G,
  Giommi A, editors. Topics in theoretical and applied statistics. Cham: Springer International
  Publishing;  p. 25–35.
  
   Marin, JM.,  Mengersen, K., and  Robert, C. (2005). Bayesian modeling and inference on mixtures of distributions. Handbook of Statistics, 25, 459–503.
  
   North, G. R., Wang, J., and Genton, M. G. (2011) Correlation models for temperature fields
  Journal of Climate, 24, 5850-5862.
  
  
  
   Wang, J., Wen, S., Symmans, W.F., Pusztai, L., and Coombes, K.R. (2009). The bimodality index: a criterion
  for discovering and ranking bimodal signatures from cancer gene expression profiling data. Cancer
  Inform, 7, 199–216.
  
   Xu, G., and Genton, M.G. (2017). Tukey g-and-h random fields. Journal of the American Statistical Association, 112(519), 1236-1249
  
   Zareifard, H., and Jafari Khaledi, M. (2013). Non-Gaussian modelling of  spatial data using scale mixing of
  a unified skew Gaussian process. Journal of Multivariate Analysis, 114, 16-28.
  
   Zhang, C., Mapes, B.E., and Soden, B.J. (2003). Bimodality in tropical water vapor.Quart. J. Roy. Meteorol. Soc, 129(594), 2847–2866.
  
  
   Zhu, X., and Genton, MG.(2012). Short-term wind speed forecasting for power system operation. International Statistical Review , 38, 2-23. Available: http://dx.doi.org/10.1111/j.1751-5823.2011.00168.x.
  
   Wolfram Research. (2020). Wolfram alpha. Wolfram Research, Inc. Champaign, IL.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § APPENDIX 1
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Zhang, H., and El-Shaarawi, A. (2010). On spatial skew-Gaussian processes and applications.
  Environmetrics. 21(1), 33-47.
  
  
  
  
  
  