Robust finite mixture modeling of multivariate unrestricted skew-normal generalized hyperbolic distributions

Abstract

In this paper, we introduce an unrestricted skew-normal generalized hyperbolic (SUNGH) distribution for use in finite mixture modeling or clustering problems. The SUNGH is a broad class of flexible distributions that includes various other well-known asymmetric and symmetric families such as the scale mixtures of skew-normal, the skew-normal generalized hyperbolic and its corresponding symmetric versions. The class of distributions provides a much needed unified framework where the choice of the best fitting distribution can proceed quite naturally through either parameter estimation or by placing constraints on specific parameters and assessing through model choice criteria. The class has several desirable properties, including an analytically tractable density and ease of computation for simulation and estimation of parameters. We illustrate the flexibility of the proposed class of distributions in a mixture modeling context using a Bayesian framework and assess the performance using simulated and real data.

Appendix

1.1 A.1. Proof of Propositions 1 to 6

In this appendix, we prove Propositions 1 to 6.

Proof of Proposition 1

By considering (7),

• (a):
• (b):

$\square$

Proof of Proposition 2

By considering the stochastic representation (7) and the fact that $\varvec{W}_{0}$ (and so ${\varvec{W}}$) are uncorrelated, this subject proved. In the case of $\varvec{\Lambda }^{*}=\left( {{\begin{array}{cc} {\varvec{\Lambda }_{p\times q} }&{} {\mathbf{0}_{p\times m} } \\ \end{array} }} \right)$, relation (7) for ${\varvec{Y}}\sim \mathrm{SUNGH}_{p,q+m} \left( {{\varvec{\mu }} ,\varvec{\Sigma },\varvec{\Lambda }^{*},\varpi } \right)$ is equivalent to ${\varvec{Y}}={\varvec{\mu }} +\varvec{\Lambda }^{*}{\varvec{W}}+\kappa \left( U \right) ^{1/2}\varvec{\Sigma }^{1/2}{\varvec{W}}_1 ={\varvec{\mu }} +\varvec{\Lambda }{\varvec{W}}^{\left( 1 \right) }+\kappa \left( U \right) ^{1/2}\varvec{\Sigma }^{1/2}{\varvec{W}}_1$, where ${\varvec{W}}^{\left( 1 \right) }$ is the first q components of W, and in the case of $\varvec{\Lambda }^{*}=\left( {{\begin{array}{cc} {\mathbf{0}_{p\times m} }&{} {\varvec{\Lambda }_{p\times q} } \\ \end{array} }} \right)$, relation (7) for ${\varvec{Y}}\sim \mathrm{SUNGH}_{p,q+m} \left( {{\varvec{\mu }} ,\varvec{\Sigma },\varvec{\Lambda }^{*},{\varvec{\varpi }} } \right)$ is equivalent to ${\varvec{Y}}={\varvec{\mu }} +\varvec{\Lambda }^{*}{\varvec{W}}+\kappa \left( U \right) ^{1/2}\varvec{\Sigma }^{1/2}{\varvec{W}}_1 ={\varvec{\mu }} +\varvec{\Lambda }{\varvec{W}}^{\left( 2 \right) }+\kappa \left( U \right) ^{1/2}\varvec{\Sigma }^{1/2}{\varvec{W}}_1$, where ${\varvec{W}}^{\left( 2 \right) }$ is the last q components of ${\varvec{W}}$$\square$

Proof of Proposition 3

By considering the stochastic representation (7), we have that ${\varvec{b}}+{\varvec{BY}}={\varvec{b}}+{\varvec{B}}{\varvec{\mu }} +{\varvec{B}}\varvec{\Lambda }{\varvec{W}}+\kappa \left( U \right) ^{1/2}\left( {{\varvec{B}}\varvec{\Sigma }{\varvec{B}}^{\top }} \right) ^{1/2}{\varvec{W}}_1$$\square$

Proof of Proposition 4

By considering Proposition 3, with ${\varvec{b}}=\mathbf{0}$ and the matrix ${\varvec{B}}$ in the form of $\left( {{\begin{array}{cc} {{\varvec{I}}_{p_1 } }&{} {{\varvec{0}}_{p_1 \times p_2 } } \\ \end{array} }} \right)$ or $\left( {{\begin{array}{cc} {{\varvec{0}}_{p_2 \times p_1 } }&{} {{\varvec{I}}_{p_2 } } \\ \end{array} }} \right)$, respectively, this subject proved $\square$

Proof of Proposition 5

Since ${\varvec{Y}}=\left( {{\varvec{Y}}_1^\top ,{\varvec{Y}}_2^\top } \right) ^{\top }$, from part b) of the Proposition 1, we have $\hbox {Var}\left[ {\varvec{Y}} \right] =\left( {\hbox {Cov}\left( {{\varvec{Y}}_i ,{\varvec{Y}}_j } \right) } \right) _{i,j=1,2} =\left( {\varvec{\Sigma }_{ij} +\varvec{\Lambda }_i \left[ {\left( {k_2 -k_1^2 } \right) \frac{2}{\pi }{} \mathbf{1}_q \mathbf{1}_q^\top -\frac{2}{\pi }k_2 I_q } \right] \varvec{\Lambda }_j^\top } \right)$. Thus, if $\varvec{\Sigma }_{12} =\mathbf{0}$, then $\hbox {Cov}\left( {{\varvec{Y}}_1 ,{\varvec{Y}}_2 } \right) =\varvec{\Lambda }_1 \big [ \left( {k_2 -k_1^2 } \right) \frac{2}{\pi }{} \mathbf{1}_q \mathbf{1}_q^\top -\frac{2}{\pi }k_2 I_q \big ]\varvec{\Lambda }_2^\top$, thus following that each of the conditions $\varvec{\Lambda }_1 =\mathbf{0}$ or $\varvec{\Lambda }_2 =\mathbf{0}$ leads to $\hbox {Cov}\left( {{\varvec{Y}}_1 ,{\varvec{Y}}_2 } \right) =\mathbf{0}$$\square$

Proof of Proposition 6

The first part follows by applying Proposition 2 in the Proposition 4. For the proof of the second result, note from the proof of Proposition 5 that

$$\begin{aligned}&\hbox {Cov}\left( {{\varvec{Y}}_1 ,{\varvec{Y}}_2 } \right) =\left( {{\begin{array}{cc} {\varvec{\Lambda }_{11} }&{} {{\varvec{0}}_{p_1 \times q_2 } } \\ \end{array} }} \right) _{p_1 \times q}\\&\quad \left[ {\left( {k_2 -k_1^2 } \right) \frac{2}{\pi }{} \mathbf{1}_q \mathbf{1}_q^\top -\frac{2}{\pi }k_2 {\varvec{I}}_q } \right] _{q\times q} \left( {{\begin{array}{cc} {\mathbf{0}_{p_2 \times q_1 } }&{} {\varvec{\Lambda }_{22} } \\ \end{array} }} \right) _{q\times p_2 }^\top . \end{aligned}$$

Thus, using the partitions ${\varvec{I}}_q =\hbox {diag}\left( {{\varvec{I}}_{q_1 } ,{\varvec{I}}_{q_2 } } \right)$ and $\mathbf{1}_q =\left( {\mathbf{1}_{q_1 }^\top ,\mathbf{1}_{q_2 }^\top } \right) ^{\top }$ we obtain the proof $\square$

1.2 A.2. Matrix variate priors for skewness matrix

Considering the matrix variate priors in the form of $\varvec{\Lambda }_k \sim MN_{p,q} \left( {{\varvec{N}}_k ,{\varvec{S}}_k ,{\varvec{F}}_k } \right) ,k=1,\ldots ,K$, where MN denotes the matrix normal distributions, this leads to the following posteriors instead of (19) as follows:

$\left. {\hbox {vec}(\varvec{\Lambda }_k )} \right| \varvec{\Theta }_{\left( {-\varvec{\Lambda }_k } \right) } ,{\varvec{y}},{\varvec{u}},{\varvec{w}},z_i =k\sim N_{pq} \left( {{\varvec{\mu }} ,\varvec{\Sigma }} \right) ;k=1,\ldots ,K$, where

$$\begin{aligned} {\varvec{\mu }}&=\varvec{\Sigma }\left[ {\left( {\mathbf{S}_k \otimes {\varvec{F}}_k } \right) ^{-1}\hbox {vec}\left( {{\varvec{N}}_k } \right) +\mathop \sum \limits _{B_k } \kappa \left( {u_{ik} } \right) ^{-1}\left( {{\varvec{M}}_{ik}^\top \otimes \varvec{\Sigma }_k^{-1} } \right) } \right] ,\\&\varvec{\Sigma }=\left[ {\left( {\mathbf{S}_k \otimes {\varvec{F}}_k } \right) ^{-1}+\mathop \sum \limits _{B_k } \kappa \left( {u_{ik} } \right) ^{-1}\left( {\varvec{\Sigma }_k^{-1}\otimes {\varvec{L}}_{ik} } \right) } \right] ^{-1}, \end{aligned}$$

(19a)

where ${\varvec{L}}_{ik} ={\varvec{w}}_{ik} {\varvec{w}}_{ik}^\top$ and ${\varvec{M}}_{ik} =\left( {{\varvec{y}}_i -{\varvec{\mu }} _k } \right) {\varvec{w}}_{ik}^\top$, for which $\otimes$ denotes the Kronecker product and $\hbox {vec}$ denotes the vectorization of a matrix (a linear transformation which converts the matrix into a column vector).

Using these forms for the Gibbs updates may improve mixing and convergence to a stationary distribution. However, they involve the use of matrix variate distributions for which users may not be familiar; hence, a simpler (computational) update is provided in the main text.