A Matrix Variate Skew-t Distribution

TL;DR

This paper introduces a matrix variate skew-t distribution derived from a mean-variance matrix normal mixture, along with an EM algorithm for parameter estimation, addressing a gap in modeling three-way data.

Contribution

It proposes the first matrix variate skew-t distribution and an associated estimation algorithm, expanding skewness modeling in three-way data analysis.

Findings

Demonstrates the distribution's applicability through simulations

Provides an EM algorithm for efficient parameter estimation

Addresses the lack of skew-t models in matrix variate settings

Abstract

Although there is ample work in the literature dealing with skewness in the multivariate setting, there is a relative paucity of work in the matrix variate paradigm. Such work is, for example, useful for modelling three-way data. A matrix variate skew-t distribution is derived based on a mean-variance matrix normal mixture. An expectation-conditional maximization algorithm is developed for parameter estimation. Simulated data are used for illustration.

Figures

Tables

Table 1. Table 1: Component wise averages and standard deviations for the estimated parameters for simulations 1 and 2.

Simulation	$𝐌\left(sd\right)$	$𝐀\left(sd\right)$	$𝚺\left(sd\right)$	$𝚿\left(sd\right)$	$\nu \left(sd\right)$
1					4.22 (0.63)
2					4.22 (0.92)

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(\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(x|\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}{x}^{-\alpha -1}\exp \left\{-\frac{\beta }{x}\right\}.$$

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

$$K_{\lambda }(x)=\frac{1}{2}{\int }_0^{\infty }{y}^{\lambda -1}\exp \left\{-\frac{x}{2}\left(y+\frac{1}{y}\right)\right\}dy$$

$$\mathbb{E}(Y)=\sqrt{\frac{b}{a}}\frac{K_{\lambda +1}(\sqrt{ab})}{K_{\lambda }(\sqrt{ab})},$$

$$\mathbb{E}\left(\frac{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}),$$

$$\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(\mathbf{X},w|\mathbf{\vartheta })$$

$$=\frac{{\frac{\nu }{2}}^{\frac{\nu }{2}}}{(2\pi {)}^{\frac{np}{2}}|\mathbf{\Sigma }{|}^{\frac{p}{2}}|\mathbf{\Psi }{|}^{\frac{n}{2}}\Gamma (\frac{\nu }{2})}{w}^{-\frac{\nu +np}{2}-1}$$

$$\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)+\nu \right]\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 })+\nu }{w}+w\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi })\right]\right\},$$

$$\begin{array}[]{ccc}\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})&\text{and}&\rho(\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi})=\,\mbox{tr}(\mathbf{\Sigma}^{-1}\mathbf{A}\mathbf{\Psi}^{-1}\mathbf{A}^{\prime}).\end{array}$$

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

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

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

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

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

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

$$f$$

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

$$W|\mathbf{X}\sim GIG\left(\rho (\mathbf{A},\mathbf{\Sigma },\mathbf{\Psi }),\delta (\mathbf{X};\mathbf{M},\mathbf{\Sigma },\mathbf{\Psi })+\nu ,\lambda \right),$$

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

$${\ell }_c(\mathbf{\vartheta })=$$

$$+\frac{1}{2}\sum _{i=1}^N\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}({\mathbf{X}}_i-\mathbf{M}){\mathbf{\Psi }}^{-1}{\mathbf{A}}^{\prime }\right)+\frac{1}{2}\sum _{i=1}^N\text{mbox}tr\left({\mathbf{\Sigma }}^{-1}\mathbf{A}{\mathbf{\Psi }}^{-1}({\mathbf{X}}_i-\mathbf{M}{)}^{\prime }\right)$$

$$-\frac{1}{2}\sum _{i=1}^N\frac{1}{w_i}\left[\text{mbox}tr({\mathbf{\Sigma }}^{-1}({\mathbf{X}}_i-\mathbf{M}){\mathbf{\Psi }}^{-1}({\mathbf{X}}_i-\mathbf{M}{)}^{\prime })+\nu \right]-\frac{1}{2}\sum _{i=1}^N{w}_i\text{mbox}tr({\mathbf{\Sigma }}^{-1}\mathbf{A}{\mathbf{\Psi }}^{-1}{\mathbf{A}}^{\prime }),$$

$$a_i^{(t+1)}$$

$$=\sqrt{\frac{\delta ({\mathbf{X}}_i;{\hat{\mathbf{M}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})+{\hat{\nu }}^{(t)}}{\rho ({\hat{\mathbf{A}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})}}\frac{K_{{\lambda }^{(t)}+1}({\kappa }^{(t)})}{K_{{\lambda }^{(t)}}({\kappa }^{(t)})}$$

$$b_i^{(t+1)}$$

$$=\sqrt{\frac{\rho ({\hat{\mathbf{A}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})}{\delta ({\mathbf{X}}_i;{\hat{\mathbf{M}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})+{\hat{\nu }}^{(t)}}}\frac{K_{{\lambda }^{(t)}+1}({\kappa }^{(t)})}{K_{{\lambda }^{(t)}}({\kappa }^{(t)})}+\frac{{\hat{\nu }}^{(t)}+np}{\delta ({\mathbf{X}}_i;{\hat{\mathbf{M}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})+{\nu }^{(t)}}$$

$$c_i^{(t+1)}$$

$$=\log \left(\sqrt{\frac{\delta ({\mathbf{X}}_i;{\hat{\mathbf{M}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})+{\hat{\nu }}^{(t)}}{\rho ({\hat{\mathbf{A}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})}}\right)+\frac{1}{K_{{\lambda }^{(t)}}({\kappa }^{(t)})}{\frac{\partial }{\partial \lambda }{K}_{\lambda }({\kappa }^{(t)})|}_{\lambda ={\lambda }^{(t)}}$$

$${\kappa }^{(t)}=\sqrt{[\rho ({\hat{\mathbf{A}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})][\delta ({\mathbf{X}}_i;{\hat{\mathbf{M}}}^{(t)},{\hat{\mathbf{\Sigma }}}^{(t)},{\hat{\mathbf{\Psi }}}^{(t)})+{\hat{\nu }}^{(t)}]},$$

Full text

Keywords: Matrix variate distribution; skew-$t$ distribution

1 Introduction

Matrix variate distributions have proven to be useful for modelling three-way data, such as multivariate longitudinal data. However, in most cases, the underlying distribution has been elliptical such as the matrix variate normal and the matrix variate $t$ distributions. However, there has been relatively little work done on matrix variate data that can account for skewness present in the data. The work that has been carried out in the area of matrix variate skew distributions is mostly limited to the matrix variate skew-normal distribution. Herein, we derive a matrix variate skew-$t$ distribution. The remainder of this paper is laid out as follows. In Section 2, some background is presented. In Section 3, the density of the matrix variate skew-$t$ distribution is derived and a parameter estimation procedure is given. Section 4 looks at some simulations, and we conclude with a summary and some future work (Section 5).

2 Background

2.1 Matrix Variate Distributions

One natural method to model three-way data is to use a matrix-variate distribution. There are many examples in the literature of such distributions, the most well-known being the matrix-normal distribution. For notional clarity, we use $\mathbf{X}$ to denote a realization of a random matrix $\mathscr{X}$. An $n\times p$ random matrix $\mathscr{X}$ follows a matrix variate normal distribution with location parameter $\mathbf{M}$ and scale matrices $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ of dimensions $n\times n$ and $p\times p$, respectively. We write $\mathscr{X}\sim\mathcal{N}_{n\times p}(\mathbf{M},\mathbf{\Sigma},\mathbf{\Psi})$ to denote such a random matrix and the density of $\mathscr{X}$ can be written

[TABLE]

One well known property of the matrix variate normal distribution (Harrar and Gupta, 2008) is

[TABLE]

where $\mathcal{N}_{np}(\cdot)$ is the multivariate normal density with dimension $np$, $\text{vec}(\mathbf{M})$ is the vectorization of $\mathbf{M}$, and $\otimes$ is the Kronecker product.

Although the matrix variate normal is arguably the most mathematically tractable, there are examples of non-normal cases. One famous example is the Wishart distribution (Wishart, 1928) arising as the distribution of the sample covariance matrix of a multivariate normal sample. More recently, however, there has been some work done in the area of matrix skew distributions such as the matrix-variate skew normal distribution, e.g., Chen and Gupta (2005), Domínguez-Molina et al. (2007), and Harrar and Gupta (2008). More information on matrix variate distributions can be found in Gupta and Nagar (1999). Very recently, there has also been work done in the area of finite mixtures. Specifically, Anderlucci et al. (2015) looked at clustering and classification of multivariate longitudinal data using a mixture of matrix variate normal distributions. Also, Doğru et al. (2016), looked at mixtures of matrix variate $t$ distributions.

2.2 Normal Variance-Mean Mixtures

Various multivariate distributions such as the multivariate $t$, and skew-$t$, the shifted asymmetric Laplace distribution, and the generalized hyperbolic distributions arise as special cases of a normal variance-mean mixture (cf. McNicholas, 2016, Ch. 6). In this formulation, the density of a $p$-dimensional random vector $\mathbf{X}$ takes the form

[TABLE]

which is equivalent to the representation

[TABLE]

where $\mathbf{V}\sim\mathcal{N}_{p}({\bf 0},\mathbf{\Sigma})$ and $W>0$ is a latent random variable with density $h(w|{\boldsymbol{\theta}})$. The multivariate skew-$t$ distribution with $\nu$ degrees of freedom arises as a special case with $W\sim\text{IG}\left(\frac{\nu}{2},\frac{\nu}{2}\right)$, where $\text{IG}(\cdot)$ denotes the inverse Gamma distribution with density function

[TABLE]

2.3 The Generalized Inverse Gaussian Distribution

A random variable $Y$ has a generalized inverse Gaussian (GIG) distribution with parameters $a,b$ and $\lambda$ if its density function can be written as

[TABLE]

where

[TABLE]

is the modified Bessel function of the third kind with index $\lambda$. Several functions of GIG random variables have tractable expected values, e.g.,

[TABLE]

[TABLE]

[TABLE]

where $a,b\in\mathbb{R}^{+}$, $\lambda\in\mathbb{R}$, and $K_{\lambda}(\cdot)$ is the modified Bessel function of the third kind with index $\lambda$. These results will prove to be useful for parameter estimation for the matrix-variate skew-$t$ distribution.

3 Methodology

3.1 A Matrix Variate Skew-$t$ Distribution

We will say that an $n\times p$ random matrix $\mathscr{X}$ has a matrix variate skew-$t$ distribution, $\text{MVST}_{n\times p}(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\nu)$, if $\mathscr{X}$ can be written

[TABLE]

where $\mathbf{M}$ and $\mathbf{A}$ are $n\times p$ matrices, $\mathscr{V}\sim\mathcal{N}_{n\times p}\left(\bf{0},\mathbf{\Sigma},\mathbf{\Psi}\right)$, and $W\sim\text{IG}\left(\frac{\nu}{2},\frac{\nu}{2}\right)$. Analogous to its multivariate counterpart, $\mathbf{M}$ is a location matrix, $\mathbf{A}$ is a skewness matrix, $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are scale matrices, and $\nu$ is the degrees of freedom. It then follows that

[TABLE]

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

[TABLE]

where ${\boldsymbol{\vartheta}}=(\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\nu)$. We note that the exponential term in (8) can be written as

[TABLE]

where

[TABLE]

Therefore, the marginal density of $\mathscr{X}$ is

[TABLE]

Making the change of variables given by

[TABLE]

we can write

[TABLE]

The density of $\mathscr{X}$, as derived here, is considered a matrix variate extension of the multivariate skew-$t$ density used by Murray et al. (2014a, b). As discussed by Dutilleul (1999) and Anderlucci et al. (2015) in the matrix variate normal case, the estimates of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are unique only up to a multiplicative constant. Indeed, if we let $\tilde{\mathbf{\Sigma}}=v\mathbf{\Sigma}$ and $\tilde{\mathbf{\Psi}}=(1/v)\mathbf{\Psi}$, $v\neq 0$, the likelihood is unchanged. This identifiability issue can be resolved, for example, by setting $\,\mbox{tr}(\mathbf{\Sigma})=n$ or $\,\mbox{tr}(\mathbf{\Psi})=p$. Note that $\tilde{\mathbf{\Psi}}\otimes\tilde{\mathbf{\Sigma}}=\mathbf{\Psi}\otimes\mathbf{\Sigma}$, so the estimate of the Kronecker product is unique.

For the purposes of parameter estimation, note that the conditional density of $W$ is

[TABLE]

Therefore,

[TABLE]

where $\lambda=-(\nu+np)/2$.

Finally, we note that

[TABLE]

where $\text{MST}_{np}(\cdot)$ denotes the multivariate skew-$t$ distribution with location parameter $\text{vec}(\mathbf{M})$, skewness parameter $\text{vec}(\mathbf{A})$, scale matrix $\mathbf{\Psi}\otimes\mathbf{\Sigma}$, and $\nu$ degrees of freedom. This can be easily seen from the representation given in (3) and the property of the matrix normal distribution given in (2). Note that the normal variance-mean mixture representation (7) as well as the relationship with the multivariate skew-$t$ distribution (9) present two convenient methods to generate random matrices from the matrix variate skew $t$ distribution. The former is used in Section 4.

3.2 Parameter Estimation

Suppose we observe a sample of $N$ matrices $\mathbf{X}_{1},\mathbf{X}_{2},\ldots\mathbf{X}_{N}$ from an $n\times p$ matrix variate skew-$t$ distribution. As with the multivariate skew-$t$ distribution, we proceed as if the observed data is incomplete, and introduce the latent variables $W_{1},\ldots,W_{n}$. The complete-data log-likelihood is then

[TABLE]

where $C$ is constant with respect to the parameters. We proceed by using an expectation-conditional maximization (ECM) algorithm (Meng and Rubin, 1993) described overleaf.

1) Initialization: Initialize the parameters $\mathbf{M},\mathbf{A},\mathbf{\Sigma},\mathbf{\Psi},\nu$. Set $t=0$

2) E Step: Update $a_{i},b_{i},c_{i}$, where

[TABLE]

where

[TABLE]

and

[TABLE]

3) First CM Step: Update the parameters $\mathbf{M},\mathbf{A},\nu$.

[TABLE]

where $\overline{a}^{(t+1)}=({1}/{N})\sum_{i=1}^{N}a_{i}^{(t+1)}$ and $\overline{b}^{(t+1)}=({1}/{N})\sum_{i=1}^{N}b_{i}^{(t+1)}$.

The update for the degrees of freedom cannot be obtained in closed form. Instead we solve (15) for $\nu$ to obtain $\hat{\nu}^{(t+1)}$.

[TABLE]

where $\varphi(\cdot)$ is the digamma function.

4) Second CM Step: Update $\mathbf{\Sigma}$

[TABLE]

5) Third CM Step: Update $\mathbf{\Psi}$

[TABLE]

6) Check Convergence: If not converged, set $t=t+1$ and return to step 2.

Note that there are several options for determining convergence of this ECM algorithm. In the simulations in Section 4, a criterion based on the Aitken acceleration (Aitken, 1926) is used. The Aitken acceleration at iteration $t$ is

[TABLE]

where $l^{(t)}$ is the (observed) log-likelihood at iteration $t$. The quantity in (18) can be used to derive an asymptotic estimate (i.e., an estimate of the value after very many iterations) of the log-likelihood at iteration $t+1$, i.e.,

[TABLE]

(cf. Böhning et al., 1994; Lindsay, 1995). As in McNicholas et al. (2010), we stop our EM algorithms when $l_{\infty}^{(t+1)}-l^{(t)}<\epsilon$, provided this difference is positive.

As discussed in Dutilleul (1999) and Anderlucci et al. (2015) for parameter estimation in the matrix variate normal case, the estimates of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are unique only up to a multiplicative constant. Indeed, if we let $\tilde{\mathbf{\Sigma}}=v\mathbf{\Sigma}$ and $\tilde{\mathbf{\Psi}}=(1/v)\mathbf{\Psi}$, $v\neq 0$, the likelihood is unchanged. However, we notice that, $\tilde{\mathbf{\Psi}}\otimes\tilde{\mathbf{\Sigma}}=\mathbf{\Psi}\otimes\mathbf{\Sigma}$, so the estimate of the Kronecker product is unique.

4 Simulations

We conduct two simulations to illustrate the estimation of the parameters. In both simulations, we take 50 different datasets of size 100, from a $3\times 4$ matrix skew-$t$ distribution. Also, in both simulations,

[TABLE]

and $\nu=4$. In simulation 1, we took the location and skewness matrix to be $\mathbf{M}_{1}$ and $\mathbf{A}_{1}$, respectively, and $\mathbf{M}_{2}$ and $\mathbf{A}_{2}$ in simulation 2, where

[TABLE]

In Figures 1 and 2, we show line plots of the marginals for each column (labelled V1, V2, V3, V4) of a typical dataset from simulations 1 and 2, respectively. The dashed red lines denote the means.

In Figure 1, the skewness in columns 1, 2, and 4, for simulation 1, is very prominent when visually compared to column 3, which has zero skewness. The skewness is also apparent in the lineplots for simulation 2, however, because the values of the skewness are generally less than those for simulation 1, it is not as prominent. The component-wise means of the parameters as well as the component wise standard deviations are given in Table 1. We see that the estimates of the mean matrix and skewness matrix are very close to the true value for both simulations. Moreover, we see that the estimates of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ also correspond approximately to the their true values, and thus so would the Kronecker product, which is not shown.

5 Discussion

The density of a matrix variate skew-$t$ distribution was derived. This distribution can be considered as a three-way extension of the multivariate skew-$t$ distribution. Parameter estimation was carried out using an ECM algorithm. Because the formulation of multivariate skew-$t$ distribution this work is based on is a special case of the generalized hyperbolic distribution, it is reasonable to postulate an extension to a broader class of matrix variate distributions. Ongoing work considers a finite mixture of matrix variate skew-$t$ distributions for clustering and classification of three-way data.

Acknowledgements

The authors gratefully acknowledge the financial support provided by the Vanier Canada Graduate Scholarships (Gallaugher) and the Canada Research Chairs program (McNicholas).