Three Skewed Matrix Variate Distributions

TL;DR

This paper introduces three new matrix variate skew distributions, extending the normal distribution to better model three-way data with skewness and flexibility, including parameter estimation and simulation.

Contribution

It presents three novel matrix variate skew distributions, their properties, and methods for parameter estimation, filling a gap in the modeling of skewed three-way data.

Findings

Derived moment generating functions for the distributions.

Presented maximum likelihood estimation methods.

Illustrated with simulated data.

Abstract

Three-way data can be conveniently modelled by using matrix variate distributions. Although there has been a lot of work for the matrix variate normal distribution, there is little work in the area of matrix skew distributions. Three matrix variate distributions that incorporate skewness, as well as other flexible properties such as concentration, are discussed. Equivalences to multivariate analogues are presented, and moment generating functions are derived. Maximum likelihood parameter estimation is discussed, and simulated data is used for illustration.

Tables

Table 1. Table 1 : Component-wise averages and standard deviations for the estimated parameters for each of the three distributions.

Generalized Hyperbolic
$𝐌$ (sd)	$𝐀$ (sd)	$𝚺$ (sd)	$𝚿$ (sd)	$\lambda$ (sd)	$\omega$ (sd)
				1.63 (2.42)	4.08 (1.33)
Variance-Gamma
$𝐌$ (sd)	$𝐀$ (sd)	$𝚺$ (sd)	$𝚿$ (sd)	$\gamma$ (sd)
				4.20 (1.04)
Normal Inverse Gaussian
$𝐌$ (sd)	$𝐀$ (sd)	$𝚺$ (sd)	$𝚿$ (sd)	$\tilde{\gamma }$ (sd)
				2.12 (0.50)

Equations

$$f(\mathbf{X}|\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })=\frac{1}{(2\pi {)}^{\frac{np}{2}}|\mathbf{\Sigma }{|}^{\frac{p}{2}}|\mathbf{\Psi }{|}^{\frac{n}{2}}}\exp \left\{-\frac{1}{2}\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}(\mathbf{X}-\mathbf{M}){\mathbf{\Psi }}^{-1}(\mathbf{X}-\mathbf{M}{)}^{\prime }\right)\right\},$$

$$\mathcal{X}\sim {\mathcal{N}}_{n\times p}(\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })⟺\text{vec}(\mathcal{X})\sim {\mathcal{N}}_{np}(\text{vec}(\mathbf{M}),\mathbf{\Psi }\otimes \mathbf{\Sigma }),$$

$$f(y|\delta ,\gamma )=\frac{\delta }{\sqrt{2\pi }}\exp \{\delta \gamma \}{y}^{-\frac{3}{2}}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{y}+{\gamma }^{2}y\right)\right\}$$

$$f(y|a,b,\lambda )=\frac{{\left(a/b\right)}^{\frac{\lambda }{2}}{y}^{\lambda -1}}{2K_{\lambda }(\sqrt{ab})}\exp \left\{-\frac{ay+b/y}{2}\right\}$$

$$K_{\lambda }(u)=\frac{1}{2}{\int }_0^{\infty }{z}^{\lambda -1}\exp \left\{-\frac{u}{2}\left(z+\frac{1}{z}\right)\right\}dz$$

$$\begin{array}{rl}& \mathbb{E}(Y)=\sqrt{\frac{b}{a}}\frac{K_{\lambda +1}(\sqrt{ab})}{K_{\lambda }(\sqrt{ab})},\mathbb{E}\left(1/Y\right)=\sqrt{\frac{a}{b}}\frac{K_{\lambda +1}(\sqrt{ab})}{K_{\lambda }(\sqrt{ab})}-\frac{2\lambda }{b},\\ & \mathbb{E}(\log Y)=\log \left(\sqrt{\frac{b}{a}}\right)+\frac{1}{K_{\lambda }(\sqrt{ab})}\frac{\partial }{\partial \lambda }{K}_{\lambda }(\sqrt{ab}).\end{array}$$

$$g(y|\omega ,\eta ,\lambda )=\frac{{\left(w/\eta \right)}^{\lambda -1}}{2\eta K_{\lambda }(\omega )}\exp \left\{-\frac{\omega }{2}\left(\frac{w}{\eta }+\frac{\eta }{w}\right)\right\},$$

$$f(\mathbf{x})=\int_{0}^{\infty}\phi_{p}(\mathbf{x}|\mbox{\boldmath$\mu$}+w\boldsymbol{\alpha},w\mathbf{\Sigma})h(w|{\boldsymbol{\theta}})dw,$$

$$\mathbf{X}=\mbox{\boldmath$\mu$}+W\boldsymbol{\alpha}+\sqrt{W}\mathbf{V},$$

$$f(w|a,b)=\frac{b^a}{\Gamma (a)}{w}^{a-1}\exp \{-bw\}$$

$$\mathcal{X}=\mathbf{M}+W\mathbf{A}+\sqrt{W}\mathcal{V},$$

$$\mathcal{X}|w\sim {\mathcal{N}}_{n\times p}\left(\mathbf{M}+w\mathbf{A},w\mathbf{\Sigma },\mathbf{\Psi }\right)$$

$$f$$

$$\times \exp \left\{-\frac{1}{2w}\left(\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}){\mathbf{\Psi }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}{)}^{\prime }\right)+\omega \right)-\omega w/2\right\},$$

$$\exp \left\{\text{mbox}tr({\mathbf{\Sigma }}^{-1}(\mathbf{X}-\mathbf{M}){\mathbf{\Psi }}^{-1}{\mathbf{A}}^{\prime })\right\}\times \exp \left\{-\frac{1}{2}\left[\frac{\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega }{w}+w\left(\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right)\right]\right\},$$

$$f(\mathbf{X})$$

$$\times \frac{1}{2}{\int }_0^{\infty }{w}^{\lambda -\frac{np}{2}-1}\exp \left\{-\frac{1}{2}\left[\frac{\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega }{w}+w\left(\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right)\right]\right\}dw.$$

$$y=\frac{\sqrt{\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega }}{\sqrt{\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega }}w,$$

$$f_{\text{MVGH}}(\mathbf{X}|\mathbf{\vartheta })=$$

$$\times K_{\left(\lambda -np/2\right)}\left(\sqrt{\left[\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right]\left[\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right]}\right),$$

$$f(w|$$

$$={\left(\frac{\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega }{\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega }\right)}^{\frac{\left(\lambda -np/2\right)}{2}}\frac{w^{\lambda -np/2-1}}{2K_{\left(\lambda -np/2\right)}\sqrt{\left[\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right]\left[\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega \right]}}$$

$$\times \exp \left\{-\frac{(\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+\omega )w+[\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\omega ]/w}{2}\right\}.$$

$$\begin{array}{rl}f(\mathbf{X},w|\mathbf{\vartheta })=& \frac{{\gamma }^{\gamma }}{(2\pi {)}^{\frac{np}{2}}|\mathbf{\Sigma }{|}^{\frac{p}{2}}|\mathbf{\Psi }{|}^{\frac{n}{2}}\Gamma (\gamma )}{w}^{\gamma -\frac{np}{2}-1}\\ & \times \exp \left\{-\frac{1}{2w}\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}){\mathbf{\Psi }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}{)}^{\prime }\right)-\gamma w\right\}.\end{array}$$

$$f_{\text{MVVG}}(\mathbf{X}|\mathbf{\vartheta })=$$

$$\times K_{\left(\gamma -\frac{np}{2}\right)}\left(\sqrt{\left[\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+2\gamma \right]\left[\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })\right]}\right),$$

$$f(\mathbf{X},w|\mathbf{\vartheta })$$

$$\times \exp \left\{-\frac{1}{2w}\left(\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}){\mathbf{\Psi }}^{-1}(\mathbf{X}-\mathbf{M}-w\mathbf{A}{)}^{\prime }\right)+1\right)-\frac{w{\tilde{\gamma }}^2}{2}+\tilde{\gamma }\right\},$$

$$f_{\text{MVNIG}}(\mathbf{X}|\mathbf{\vartheta })$$

$$\times K_{-(1+np)/2}\left(\sqrt{\left[\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })+{\tilde{\gamma }}^2\right]\left[\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+1\right]}\right),$$

$$\begin{array}{rl}\mathcal{X}\sim {\text{MVGH}}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi },\omega ,\lambda )& ⟺\text{vec}(\mathcal{X})\sim {\text{GH}}_{np}(\text{vec}(\mathbf{M}),\text{vec}(\mathbf{A}),\mathbf{\Psi }\otimes \mathbf{\Sigma },\omega ,\lambda ),\\ \mathcal{X}\sim {\text{MVVG}}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi },\gamma )& ⟺\text{vec}(\mathcal{X})\sim {\text{VG}}_{np}(\text{vec}(\mathbf{M}),\text{vec}(\mathbf{A}),\mathbf{\Psi }\otimes \mathbf{\Sigma },\gamma ),\\ \mathcal{X}\sim {\text{MVNIG}}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi },\tilde{\gamma })& ⟺\text{vec}(\mathcal{X})\sim {\text{NIG}}_{np}(\text{vec}(\mathbf{M}),\text{vec}(\mathbf{A}),\mathbf{\Psi }\otimes \mathbf{\Sigma },\tilde{\gamma }).\end{array}$$

$$M_{\mathcal{X}}(\mathbf{T})$$

$$=\exp \{\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{M})\}\mathbb{E}[\exp \{W\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{A}+\mathbf{T\Sigma }{\mathbf{T}}^{\prime }\mathbf{\Psi })\}]$$

$$=\exp \{\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{M})\}{M}_W(\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{A}+\mathbf{T\Sigma }{\mathbf{T}}^{\prime }\mathbf{\Psi })),$$

$$\exp \{\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{M})\}{\left[1-2\frac{\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{A}+\mathbf{T\Sigma }{\mathbf{T}}^{\prime }\mathbf{\Psi })}{\omega }\right]}^{-\frac{\lambda }{2}}\frac{K_{\lambda }\left(\sqrt{\omega (\omega -2\text{mbox}tr({\mathbf{T}}^{\prime }\mathbf{A}+\mathbf{T\Sigma }{\mathbf{T}}^{\prime }\mathbf{\Psi }))}\right)}{K_{\lambda }(\omega )}.$$

Full text

1 Introduction

Matrix variate distributions are useful in modelling three way data, e.g., multivariate longitudinal data. Although the matrix normal distribution is widely used, there is relative paucity in the area of matrix skewed distributions. Herein, we discuss matrix variate extensions of three already well established multivariate distributions using matrix normal variance-mean mixtures. Specifically, we consider a matrix variate generalized hyperbolic distribution, a matrix variate variance-gamma distribution, and a matrix variate normal inverse Gaussian (NIG) distribution. Along with the matrix variate skew-$t$ distribution, mixtures of these respective distributions have been used for clustering (Gallaugher and McNicholas, 2018); however, unlike the matrix variate skew-$t$ distribution (Gallaugher and McNicholas, 2017), their properties have yet to be discussed and this letter aims to fill that gap.

2 Background

2.1 The Matrix Variate Normal and Related Distributions

One of the most mathematically tractable examples of a matrix variate distribution is the matrix variate normal distribution. An $n\times p$ random matrix $\mathscr{X}$ follows a matrix variate normal distribution if its probability density function can be written as

[TABLE]

where $\mathbf{M}$ is an $n\times p$ location matrix, $\mathbf{\Sigma}$ is an $n\times n$ scale matrix for the rows of $\mathscr{X}$ and $\mathbf{\Psi}$ is a $p\times p$ scale matrix for the columns of $\mathscr{X}$. We denote this distribution by $\mathcal{N}_{n\times p}(\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi})$ and, for notational clarity, we will denote the random matrix by $\mathscr{X}$ and its realization by $\mathbf{X}$. One useful property of the matrix variate normal distribution, as given in Harrar and Gupta (2008), is

[TABLE]

where $\mathcal{N}_{np}(\cdot)$ denotes the multivariate normal distribution with dimension $np$, $\text{vec}(\cdot)$ denotes the vectorization operator, and $\otimes$ denotes the Kronecker product.

Although the matrix variate normal is probably the most well known matrix variate distribution, there are other examples. For example, the Wishart distribution (Wishart, 1928) was shown to be the distribution of the sample covariance matrix for a random sample from a multivariate normal distribution. There are also a few examples of a matrix variate skew normal distribution such as Chen and Gupta (2005), Domínguez-Molina et al. (2007) and Harrar and Gupta (2008). Most recently, Gallaugher and McNicholas (2017), considered a matrix variate skew-$t$ distribution using a matrix normal variance-mean mixture.

There are also a few examples of mixtures of matrix variate distributions. Anderlucci et al. (2015) considered a mixture of matrix variate normal distributions for clustering multivariate longitudinal data and Doğru et al. (2016) considered a mixture of matrix variate $t$ distributions.

2.2 The Inverse and Generalized Inverse Gaussian Distributions

The derivation of the matrix distributions and parameter estimation discussed in Section 3, will rely heavily on the generalized inverse Gaussian distribution, and to a lesser extent the inverse Gaussian distribution. A random variable $Y$ follows an inverse Gaussian distribution if its probability density function is of the form

[TABLE]

for $y>0$, where $\delta,\gamma>0$. For notational purposes, we will denote this distribution by $\text{IG}(\delta,\gamma)$.

The generalized inverse Gaussian distribution has two different parameterizations, both of which will be useful. A random variable $Y$ has a generalized inverse Gaussian distribution parameterized by $a,b>0$ and $\lambda\in\mathbb{R}$, denoted by $\text{GIG}(a,b,\lambda)$, if its probability density function can be written as

[TABLE]

for $y>0$, where

[TABLE]

is the modified Bessel function of the third kind with index $\lambda$. Some expectations of functions of a GIG random variable with this parameterization have a mathematically tractable form, e.g.,

[TABLE]

Although this parameterization of the GIG distribution will be useful for parameter estimation, for the purposes of deriving the density of the matrix variate generalized hyperbolic distribution, it is more useful to take the parameterization

[TABLE]

where $\omega=\sqrt{ab}$ and $\eta=\sqrt{a/b}$. For notational clarity, we will denote the parameterization given in (2) by $\text{I}(\omega,\eta,\lambda)$.

2.3 Variance-Mean Mixtures

A $p$-variate random vector $\mathbf{X}$ defined in terms of a variance-mean mixture, has a probability density function of the form

[TABLE]

where the random variable $W>0$ has density function $h(w|{\boldsymbol{\theta}})$, and $\phi_{p}(\cdot)$ represents the density function of the $p$-variate Gaussian distribution. This representation is equivalent to writing

[TABLE]

where $\mu$ is a location parameter, $\boldsymbol{\alpha}$ is the skewness, $\mathbf{V}\sim\mathcal{N}_{p}({\bf 0},\mathbf{\Sigma})$ with $\mathbf{\Sigma}$ as the scale matrix, and $W$ has density function $h(w|{\boldsymbol{\theta}})$. Note that $W$ and $\mathbf{V}$ are independent. Many multivariate distributions can be obtained through a variance mean mixture by changing the distribution of $W$. For example, the $p$-dimensional generalized hyperbolic distribution, \text{GH}_{p}(\mbox{\boldmath\mu},\boldsymbol{\alpha},\mathbf{\Sigma},\psi,\chi,\lambda), as given in McNeil et al. (2005), was shown to arise as a special case of (3) by taking $W\sim\text{GIG}(\psi,\chi,\lambda)$. However, there was a restriction that $|\mathbf{\Sigma}|=1$. Simply relaxing this constraint results in an identifiability problem. In Browne and McNicholas (2015), this was discussed, and the authors proposed the reparameterization $\omega=\sqrt{\psi\chi},\eta=\sqrt{{\chi}/{\psi}}$. The representation of $\mathbf{X}$ is then as in (3), with $W\sim\text{I}(\omega,1,\lambda)$.

The $p$-dimensional variance-gamma distribution, \text{VG}_{p}(\mbox{\boldmath\mu},\boldsymbol{\alpha},\mathbf{\Sigma},\lambda,\psi), results as a limiting case of the generalized hyperbolic by taking $\lambda>0$, and $\chi\rightarrow 0$. The precise details can be found in McNicholas et al. (2017); in essence, the variance-gamma distribution also arises as a special case of (3), with $W\sim\text{gamma}(\lambda,\psi/2)$, where $\text{gamma}(a,b)$ denotes the gamma distribution with density

[TABLE]

for $w>0$, where $a,b>0$. However, we again have an identifiability issue using this representation if we remove the constraint $|\mathbf{\Sigma}|=1$. In McNicholas et al. (2017), the authors propose setting $\mathbb{E}(W)=1$, resulting in the reparameterization $\gamma\coloneqq\lambda=\psi/2$.

Finally, we have the $p$-dimensional Gaussian distribution, \text{NIG}_{p}(\mbox{\boldmath\mu},\boldsymbol{\alpha},\mathbf{\Sigma},\delta,\gamma). In Karlis and Santourian (2009), the authors derived the $p$-dimensional NIG distribution using a variance-mean mixture with $W\sim\text{IG}(\delta,\gamma)$. However, there was once again a restriction on the determinant of $\mathbf{\Sigma}$. To remove this restriction and maintain identifiability, Karlis and Santourian (2009) set $\delta=1$ and $\tilde{\gamma}=\gamma$.

3 Three Matrix Variate Skew Distributions

3.1 Matrix Normal Variance-Mean Mixture

We now present densities for matrix variate versions of the generalized hyperbolic, variance-gamma and NIG distributions. For all three of these distributions, we consider a matrix normal variance-mean mixture, where we can take the representation

[TABLE]

where $\mathscr{V}\sim\mathcal{N}_{n\times p}({\bf 0}_{n\times p},\mathbf{\Sigma},\mathbf{\Psi})$ with ${\bf 0}_{n\times p}$ representing the $n\times p$ zero matrix, $\mathbf{M}$ is an $n\times p$ location matrix, $\mathbf{A}$ is an $n\times p$ skewness matrix, and $W>0$ is a random variable with density $h({\boldsymbol{\theta}})$. We now derive three matrix variate distributions using this representation with different distributions for $W$.

3.2 A Matrix Variate Generalized Hyperbolic Distribution

We now derive the density of a matrix variate generalized hyperbolic distribution. In this case, to avoid the indentifiability issue discussed in Browne and McNicholas (2015), we take $W\sim\text{I}\left(\omega,1,\lambda\right)$, where $\omega$ is a concentration parameter and $\lambda$ is the index parameter. It then follows that

[TABLE]

and thus the joint density of $\mathbf{X}$ and $W$ is

[TABLE]

where ${\boldsymbol{\vartheta}}=(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\omega,\lambda)$.

We note that the exponential term in (5) can be written as

[TABLE]

where $\delta(\mathbf{X};\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi})=\,\mbox{tr}(\mathbf{\Sigma}^{-1}(\mathbf{X}-\mathbf{M})\mathbf{\Psi}^{-1}(\mathbf{X}-\mathbf{M})^{\prime})$ and $\rho(\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi})=\,\mbox{tr}(\mathbf{\Sigma}^{-1}\mathbf{A}\mathbf{\Psi}^{-1}\mathbf{A}^{\prime})$. Therefore, the marginal density of $\mathbf{X}$ is

[TABLE]

Making the change of variables given by

[TABLE]

(6) becomes

[TABLE]

where $\omega>0$ is a concentration parameter, and $\lambda\in\mathbb{R}$ is an index parameter.

We note that the density of $\mathscr{X}$, as derived here, resembles that of Browne and McNicholas (2015), and we denote this distribution by $\text{MVGH}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\lambda,\omega)$. For the purposes of parameter estimation, note that the conditional density of $W$ is

[TABLE]

Therefore, $W|\mathbf{X}\sim\text{GIG}\left(\rho(\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi})+\omega,\delta(\mathbf{X};\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi})+\omega,\lambda-{np}/{2}\right)$.

Note that a multiple scaled matrix variate generalized hyperbolic distribution was derived by Thabane and Safiul Haq (2004). While the distribution they derive is sometimes referred to as a matrix variate generalized hyperbolic distribution, the model of Thabane and Safiul Haq (2004) is in fact multiple scaled — a fact that may be confirmed by observing that they use a matrix variate distribution for the mixing variable $\mathbf{W}$. Not only does this mean that the distribution presented by Thabane and Safiul Haq (2004) is different to the matrix variate generalized hyperbolic distribution presented herein, but it also means that neither one of these distributions is a special case of the other. Some useful details about the multiple scaled generalized hyperbolic distribution are given by McNicholas (2016, Chp. 7).

3.3 A Matrix Variate Variance-Gamma Distribution

We now derive the density of a matrix variate variance-gamma distribution in much the same way as the generalized hyperbolic case. However, we now take $W\sim\text{gamma}(\gamma,\gamma)$, resulting in the joint distribution

[TABLE]

Following the same procedure as before, the density of $\mathscr{X}$ is then

[TABLE]

where $\gamma>0$. We will denote this distribution by $\text{MVVG}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\gamma)$. Note that $W|\mathbf{X}\sim\text{GIG}\left(\rho(\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi})+2\gamma,\delta(\mathbf{X};\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi}),\gamma-{np}/{2}\right)$.

3.4 A Matrix Variate NIG Distribution

Finally, we consider a matrix variate NIG distribution. Derived in much the same way as the previous distributions, we take $W\sim\text{IG}(1,\tilde{\gamma})$. The joint density of $\mathscr{X}$ and $W$ is

[TABLE]

and the density of $\mathscr{X}$ is then

[TABLE]

where $\tilde{\gamma}>0$. We denote this distribution by $\text{MVNIG}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\tilde{\gamma})$, and note that $W|\mathbf{X}\sim\text{GIG}\left(\rho(\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi})+\tilde{\gamma}^{2},\delta(\mathbf{X};\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi})+1,-{(1+np)}/{2}\right)$.

3.5 Some Properties

One interesting element that we see for all three of these distributions is that there is a relationship between each of these matrix variate distributions and their multivariate counterparts. Specifically,

[TABLE]

These properties can be easily seen by using the representation of $\mathscr{X}$ given in (4) as well as the property of the matrix variate normal distribution given in (1).

We can also easily derive the moment generating functions for each of these three distributions. Using the representation for a random matrix $\mathscr{X}$ given in (4) and the moment generating function for the matrix variate normal distribution given in Dutilleul (1999), we have that the moment generating function in the general case of a matrix normal variance-mean mixture is

[TABLE]

where $M_{W}(\cdot)$ is the moment generating function of $W$. Therefore, in the case of the generalized inverse Gaussian distribution, we have that the moment generating function is

[TABLE]

For the variance-gamma distribution, the moment generating function is

[TABLE]

for $\,\mbox{tr}(\mathbf{T}^{\prime}\mathbf{A}+\mathbf{T}\mathbf{\Sigma}\mathbf{T}^{\prime}\mathbf{\Psi})<\gamma$ and, in the case of the NIG distribution, the moment generating function is

[TABLE]

Parameter estimation can be performed using expectation-conditional maximization (ECM) algorithms (Meng and Rubin, 1993) by treating the data as incomplete. Details are not given here but the algorithms are equivalent to one-component versions of the ECM algorithms described by Gallaugher and McNicholas (2018).

4 Example

We now consider a simple example for each of the three different distributions. Common elements between the distributions are as follows. We take 50 datasets each with 100 observations. For each distribution, we take

[TABLE]

and the scale matrices $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are

[TABLE]

We take the additional parameters to be $\lambda_{2}=-2$, $\omega=2$ for the matrix variate generalized hyperbolic, $\gamma_{2}=4$ for the matrix variate variance-gamma and $\tilde{\gamma}_{2}=2$ for the matrix variate NIG distribution. In Figure 1 (Appendix A), we show the marginal distributions of the columns for each distribution of a typical dataset. We label the columns V1, V2, V3, and V4. The marginal location (mode) is shown by the red dashed line.

The component-wise means and standard deviations (in brackets) of the parameter estimates are given in Table 1. For all three of the distributions, we get good average estimates in general. However, one unexpected outcome is the estimate for $\lambda$ for the matrix variate generalized hyperbolic distribution. The estimate is very different from the true value, and there is a very large amount of variation. We also notice a deflation, in absolute value, of the estimated entries of the skewness matrix $\mathbf{A}$ as well as a fair amount of variation. One possible explanation is that the generalized hyperbolic distribution is over-parameterized; in which case, the deflation in the estimates for the entries of $\mathbf{A}$ could be compensation for the increased value of $\lambda$.

5 Discussion

We derived the densities and described parameter estimation for three matrix variate skew distributions using a matrix normal variance-mean mixture. The three distributions were the matrix variate generalized hyperbolic, variance-gamma and NIG distributions, respectively. When looking at the estimates in the simulations, we obtained fairly good results. One exception was the average estimates of $\lambda$ and the skewness matrix $\mathbf{A}$ for the matrix variate generalized hyperbolic distribution. However, this could be due to over-parameterization.

One possible extension of the work herein is to consider multiple-scaled analogues of the matrix variate variance-gamma, NIG and skew-$t$ distributions. The resulting multiple scaled distributions would be arrived at in an analogous fashion to the multiple-scaled matrix variate generalized hyperbolic distribution of Thabane and Safiul Haq (2004). Finally, it would be interesting to consider placing a constrained covariance structure on $\mathbf{\Sigma}$ for possible use with multivariate longitudinal data, i.e., data where multiple quantities are measured over time.

Acknowledgements

The authors are grateful for the support of the Vanier Canada Graduate Scholarships (Gallaugher) and the Canada Research Chairs (McNicholas) programs.

Appendix A Figure