Robust Classification via Finite Mixtures of Matrix Variate Skew-t Distributions

Abstract

Analysis of matrix variate data is becoming increasingly common in the literature, particularly in the field of clustering and classification. It is well known that real data, including real matrix variate data, often exhibit high levels of asymmetry. To address this issue, one common approach is to introduce a tail or skewness parameter to a symmetric distribution. In this regard, we introduce here a new distribution called the matrix variate skew-t distribution (MVST), which provides flexibility, in terms of heavy tail and skewness. We then conduct a thorough investigation of various characterizations and probabilistic properties of the MVST distribution. We also explore extensions of this distribution to a finite mixture model. To estimate the parameters of the MVST distribution, we develop an EM-type algorithm that computes maximum likelihood (ML) estimates of the model parameters. To validate the effectiveness and usefulness of the developed models and associated methods, we performed empirical experiments, using simulated data as well as three real data examples, including an application in skin cancer detection. Our results demonstrate the efficacy of the developed approach in handling asymmetric matrix variate data.

1. Introduction

The advent of modern data-collection technologies, such as electronic sensors, cell phones and web browsers, has resulted in an abundance of multivariate data sources. Much of these data can be represented as matrix variate (three-way) data, with two ways associated to the row and column dimension of each matrix variate observation and the third one representing subjects (see [1]). Matrix data can occur in different application domains, such as spatial multivariate data, longitudinal data on multiple response variables or spatio-temporal data. For this reason, statistical methods that can effectively utilize three-way data have become increasingly popular. The matrix variate normal (MVN) distribution is one of the most commonly used matrix variate elliptical distributions. However, for many real phenomena, the tails of the MVN distribution are lighter than required, with a direct impact on the corresponding model. In particular, in robust statistical analysis, heavy-tailed distributions are essential, and these include slash and t distributions. Matrix variate t (MVT) distribution has been discussed by [2], and some distributional properties of it have also been studied by [2].

Flexibility and robustness are often lacking in symmetric models when dealing with highly asymmetric data. To address this issue, a recognized method is to add a tail or skewness parameter to a symmetric distribution. Several formulations have been discussed in the literature in the form of continuous mixtures of normal variables, where a mixing variable operates on the mean or on the variance, or on both the mean and the variance of a multivariate normal variable. A general formulation was presented in [3], which encompasses a large number of existing constructions involving continuous mixtures of normal variables. Given a real-valued function $r(u,w)$ and a positive-valued function $s(u,w)$, a generalized mixture of a p variate normal distribution is given by

$$\begin{array}{c}\hfill \mathit{Y}\stackrel{d}{=}\mathit{\xi }+r(U,W)\mathit{\gamma }+s(U,W)\mathit{X},\end{array}$$

(1)

where $\stackrel{d}{=}$ denotes equality in distribution, $\mathit{\xi }\in {\mathbb{R}}^p$, $\mathit{\gamma }\in {\mathbb{R}}^p$, $\mathit{X}\sim N_p(\mathbf{0},\mathbf{\Sigma })$ and U and W are univariate random variables, with $(\mathit{X},U,W)$ being mutually independent.

The representation in (1) can be extended to the matrix variate case as

$$\begin{array}{c}\hfill \mathit{Y}\stackrel{d}{=}\mathit{M}+r(U,W)\mathsf{\Lambda }+s(U,W)\mathit{Z},\end{array}$$

(2)

where $\mathit{M}$ and $\mathsf{\Lambda }$ are $n\times p$ matrices representing the location and skewness, respectively, $\mathit{Z}\sim N_{n\times p}(\mathbf{0},\mathbf{\Sigma },\mathsf{\Psi })$ and $(\mathit{Z},U,W)$ are mutually independent. It is worth noting that the univariate nature of functions $r(u,w)$ and $s(u,W)$ simplifies the stochastic representation in (2), leading to more suitable properties for $\mathit{Y}$, and it also facilitates easier parameter estimation. Furthermore, the representation in (1) can be considered, after rearranging into a vector (denoted by $\mathrm{Vec}\left(\mathit{Y}\right)$), as 

$$\begin{array}{c}\hfill \mathrm{Vec}\left(\mathit{Y}\right)\stackrel{d}{=}\mathrm{Vec}\left(\mathit{M}\right)+r(U,W)\mathrm{Vec}\left(\mathsf{\Lambda }\right)+s(U,W)\mathrm{Vec}\left(\mathit{Z}\right).\end{array}$$

(3)

In this work, we introduce and study in detail finite mixtures of a new simple matrix variate skew-t (FM-MVST) distribution, based on (2), for dealing with clustering and classification of asymmetric and heavy-tailed matrix variate data. The proposed model’s simplicity in both density function and stochastic representation leads to a convenient strategy for parameter estimation using the expectation–conditional maximization either (ECME; [4]) algorithm, which is a variant of the EM algorithm [5]. In addition, using simulated and real datasets, we show how the proposed EM algorithm can be implemented for determining the ML estimates of the model parameters for the finite mixture of the proposed model.

The rest of this paper is organized as follows. Section 2 discusses the relevant previous research. Section 3 presents the formulation of the MVST distribution and discusses how the ECME algorithm can be proposed for ML estimation of model parameters. In Section 4, the finite mixture of MVST distributions is defined, and then, the implementation of the EM algorithm for fitting the FM-MVST model is presented. The proposed methods are illustrated by two simulation studies in Section 5 and also by the analysis of three real data datasets in Section 6. Finally, some concluding remarks and possible avenues for future work are outlined in Section 7.

2. Related Studies

To address highly asymmetric data and to have more flexibility, the functions $r(u,w)$ and $s(u,w)$ in (2) can also be considered as p-dimensional random variables. This approach may introduce complexity in both density-function and parameter-estimation issues. For instance, ref. [6] extended the scale and shape mixtures of multivariate skew normal distributions to a matrix variate setting and studied special cases and their properties. Here, we concentrate on simplifying the proposed model and the associated estimation procedure by focusing on the univariate case of $r(u,w)$ and $s(u,w)$. However, it should be noted that the proposed model here is distinct from the model presented in [6], and it cannot be viewed as a special case of it.

Based on (2), various cases have been introduced. For example, ref. [7] introduced a matrix variate skew-t distribution using the following matrix variate normal variance–mean mixture representation:

$$\mathit{Y}\stackrel{d}{=}\mathit{M}+W\mathsf{\Lambda }+W^{1/2}\mathit{Z},$$

(4)

where $\mathit{M}$ and $\mathsf{\Lambda }$ are $n\times p$ matrices representing the location and skewness, respectively, $\mathit{Z}\sim N_{n\times p}(\mathbf{0},\mathbf{\Sigma },\mathsf{\Psi })$ and $W\sim IG(\nu /2,\nu /2)$ with $IG(·)$ denoting the inverse gamma distribution. Herein, we denote the random variable with representation (4) by $\mathit{Y}\sim MVSTIG(\mathit{M},\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda },\nu )$ and its resulting density function is given by

$$\begin{array}{cc}\hfill & f_{MVSTIG}(\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda },\nu )\hfill \\ & ={\displaystyle \frac{2{(\nu /2)}^{\nu /2}exp\{\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{M}){\mathsf{\Psi }}^{-1}{\mathsf{\Lambda }}^{\top }\right]}{{\left(2\pi \right)}^{np/2}{|\mathbf{\Sigma }|}^{p/2}{\left|\mathsf{\Psi }\right|}^{n/2}\mathsf{\Gamma }(\nu /2)}}{\left({\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu }{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })}}\right)}^{-{\displaystyle \frac{\nu +np}{4}}}\hfill \\ & \times {\kappa }_{-(\nu +np)/2}\left(\sqrt{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })\left(\delta \right(\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu )}\right),\mathit{Y}\in {\mathbb{R}}^{n\times p},\hfill \end{array}$$

(5)

where

$$\begin{array}{cc}\hfill \delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })& =\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{M}){\mathsf{\Psi }}^{-1}{(\mathit{Y}-\mathit{M})}^{\top }\right],\hfill \\ \hfill \rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })& =\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}\mathsf{\Lambda }{\mathsf{\Psi }}^{-1}{\mathsf{\Lambda }}^{\top }\right],\hfill \end{array}$$

and ${\kappa }_x$ is the modified Bessel function of the third kind with index x. From (4), some other matrix variate skew distributions can be introduced by assuming different distributions for W. For more details, one may refer to [8].

Ref. [9] introduced a new family of matrix variate distributions, based on the matrix variate mean mixture of normal (MVMMN) distributions, as 

$$\begin{array}{c}\hfill \mathit{Y}\stackrel{d}{=}\mathit{M}+W\mathsf{\Lambda }+\mathit{Z}.\end{array}$$

(6)

Based on (6), three special cases, including the restricted matrix variate skew-normal (RMVSN), exponentiated MVMMN (MVMMNE) and mixed-Weibull MVMMN (MVMMNW) distributions, have been studied by using half-normal, exponential and Weibull distributions for W, respectively. Several other skew matrix variate distributions have also been discussed in the literature; see [10,11,12].

One common statistical challenge faced by researchers is identifying sub-populations or clusters within multivariate data. Recently, researchers have explored the use of finite-mixture models for matrix variate data in applications such as image analysis, genetics and neuroscience. These models offer a flexible framework for capturing complex patterns in the data and can provide insights into the underlying sub-populations or clusters; see [8,13,14,15,16,17,18].

3. Methodology

3.1. The Model

An $n\times p$-variate random matrix $\mathit{Y}$ is said to have a matrix variate skew-t (MVST) distribution, with $n\times p$ location matrix $\mathit{M}$, $n\times n$ and $p\times p$ scale matrices $\mathbf{\Sigma }$ and $\mathsf{\Psi }$, $n\times p$ shape matrix $\mathsf{\Lambda }$ and flatness parameters $\nu$, if its probability density function (pdf) is   

$$\begin{array}{cc}\hfill & f_{MVST}(\mathit{Y};\mathit{\theta })\hfill \\ & ={\displaystyle \frac{2{(\nu /2)}^{\nu /2}\mathsf{\Gamma }\left(\frac{\nu +np}{2}\right){|\mathbf{\Sigma }|}^{-p/2}{\left|\mathsf{\Psi }\right|}^{-n/2}}{{\left(2\pi \right)}^{np/2}\mathsf{\Gamma }(\nu /2)\sqrt{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}}{\left({\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}{2}}\right)}^{-{\displaystyle \frac{\nu +np}{2}}}\hfill \\ & \times T_{(\nu +np)}\left(\Delta (\mathit{Y};\mathit{\theta })\sqrt{{\displaystyle \frac{\nu +np}{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}}}\right),\mathit{Y}\in {\mathbb{R}}^{n\times p},\hfill \end{array}$$

(7)

where $\mathit{\theta }=(\mathit{M},\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda },\nu )$ denotes all the model parameters, $\Delta (\mathit{Y};\mathit{\theta })=\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{M}){\mathsf{\Psi }}^{-1}{\mathsf{\Lambda }}^{\top }\right]/\sqrt{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}$ and $T_{\nu }(·)$ denotes the cumulative distribution function (cdf) of the student’s t distribution with $\nu$ degrees of freedom. The MVST distribution reduces to the RMVSN distribution [9] when $\nu \to \infty$.

Moreover, the MVST distribution possesses the stochastic representation

$$\begin{array}{c}\hfill \mathit{Y}\stackrel{d}{=}\mathit{M}+W^{-1/2}\left(U\mathsf{\Lambda }+\mathit{Z}\right),\end{array}$$

(8)

where $\mathit{Z}\sim N_{n\times p}(\mathbf{0},\mathbf{\Sigma },\mathsf{\Psi })$, $W\sim \mathsf{\Gamma }(\nu /2,\nu /2)$ and $U\sim TN(0,1)I_{(0,\infty )}$. Herein, $TN(\mu ,{\sigma }^2)I_{\mathbb{A}}$ represents a doubly truncated normal distribution defined in the interval $\mathbb{A}=\{a_1<x<a_2\}$, and $I_{\mathbb{A}}$ denotes the indicator function of set $\mathbb{A}$. It is important to highlight that the MVSTIG distribution, represented by (4), utilizes a single mixing random variable. In contrast, the proposed MVST distribution, represented by (8), incorporates two mixing random variables, W and U, which can significantly increase the model’s flexibility.

From (8), it is easy to show that

$$\begin{array}{cc}\hfill E\left(\mathit{Y}\right)& =\mathit{M}+{\displaystyle \frac{\mathsf{\Gamma }\left(\frac{\nu -1}{2}\right)}{\mathsf{\Gamma }\left(\frac{\nu }{2}\right)}}\left({\displaystyle \frac{\nu }{\sqrt{\pi }}}\right)\mathsf{\Lambda },\hfill \\ \hfill \mathrm{Vec}\left(\mathit{Y}\right)& \sim rS{T}_{np}\left(\mathrm{Vec}\left(\mathit{M}\right),\mathsf{\Psi }\otimes \mathbf{\Sigma },\mathrm{Vec}\left(\mathsf{\Lambda }\right)\right),\hfill \end{array}$$

where ⊗ is the Kronecker product and $rS{T}_p$ denotes the p-variate restricted skew-t distribution (see [19,20]).

The stochastic representation given in (8) not only facilitates random number generation, but also enables the implementation of the EM algorithm for determining the maximum likelihood (ML) estimates of the parameters of the MVST distribution. This leads to the hierarchical representation

$$\begin{array}{cc}\hfill \mathit{Y}\left|\right(\gamma ,w)& \sim N_{n\times p}(\mathit{M}+\gamma \mathsf{\Lambda },w^{-1}\mathbf{\Sigma },\mathsf{\Psi }),\hfill \\ \hfill \gamma |w& \sim TN(0,w^{-1})I_{(0,\infty )},\hfill \\ \hfill W& \sim \mathsf{\Gamma }(\nu /2,\nu /2),\hfill \end{array}$$

(9)

where $\gamma =W^{-1/2}U$ and W are treated as latent variables. Then, $(\mathit{Y},W,\gamma )$ has the joint pdf as

$$\begin{array}{c}\hfill f_{\mathit{Y},W,\gamma }(\mathit{Y},w,\gamma )=2w^{1/2}{\varphi }_{n\times p}(\mathit{Y};\mathit{M}+\gamma \mathsf{\Lambda },w^{-1}\mathbf{\Sigma },\mathsf{\Psi })\varphi \left(w^{1/2}\gamma \right)g(w;\nu /2,\nu /2),\end{array}$$

(10)

where $\varphi (·)$ and ${\varphi }_{n\times p}(·;\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })$ are the pdfs of $N(0,1)$ and $N_{n\times p}(\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })$, respectively, and $g(·;\alpha ,\beta )$ denotes the pdf of the gamma distribution with mean $\alpha /\beta$.

Integrating out W and $\gamma$, respectively, from (10), we obtain the joint pdfs

$$\begin{array}{cc}\hfill f_{\mathit{Y},\gamma }(\mathit{Y},\gamma )& ={\displaystyle \frac{2{(\nu /2)}^{\nu /2}\mathsf{\Gamma }\left(\frac{\nu +np+1}{2}\right){|\mathbf{\Sigma }|}^{-p/2}{\left|\mathsf{\Psi }\right|}^{-n/2}}{{\left(2\pi \right)}^{(np+1)/2}\mathsf{\Gamma }(\nu /2)}}\hfill \\ & {\left({\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right){\gamma }^2-2\eta (\mathit{Y};\mathit{\theta })\gamma +\nu }{2}}\right)}^{-{\displaystyle \frac{\nu +np+1}{2}}},\hfill \end{array}$$

(11)

where $\eta (\mathit{Y};\mathit{\theta })=\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{M}){\mathsf{\Psi }}^{-1}{\mathsf{\Lambda }}^{\top }\right]$, and

$$\begin{array}{cc}\hfill f_{\mathit{Y},W}(\mathit{Y},w)& =\frac{2{(\nu /2)}^{\nu /2}{|\mathbf{\Sigma }|}^{-p/2}{\left|\mathsf{\Psi }\right|}^{-n/2}}{{\left(2\pi \right)}^{np/2}\mathsf{\Gamma }(\nu /2)\sqrt{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}{w}^{\frac{\nu +np}{2}-1}\hfill \\ & \times exp\left\{-{\displaystyle \frac{w\left(\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })\right)}{2}}\right\}\mathsf{\Phi }\left(w^{1/2}\Delta (\mathit{Y};\mathit{\theta })\right),\hfill \end{array}$$

(12)

where $\mathsf{\Phi }(·)$ denotes the cdf of the standard normal distribution.

Dividing (10) by (11), we obtain

$$\begin{array}{cc}\hfill & W\left|\right(\mathit{Y},\gamma )\hfill \\ & \sim \mathsf{\Gamma }\left({\displaystyle \frac{\nu +np+1}{2}},{\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right){\gamma }^2-2\eta (\mathit{Y};\mathit{\theta })\gamma +\nu }{2}}\right).\hfill \end{array}$$

(13)

Additionally, dividing (10) by (12), we obtain

$$\begin{array}{c}\hfill \gamma |(\mathit{Y},w)\sim TN\left({\displaystyle \frac{\eta (\mathit{Y};\mathit{\theta })}{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}},{\displaystyle \frac{w^{-1}}{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}\right)I_{(0,\infty )}.\end{array}$$

(14)

From (7) and (12), it is easy to see that

$$\begin{array}{c}\hfill f\left(w\right|\mathit{Y})=C{w}^{\frac{\nu +np}{2}-1}exp\left\{-{\displaystyle \frac{w\left(\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })\right)}{2}}\right\}\mathsf{\Phi }\left(w^{1/2}\Delta (\mathit{Y};\mathit{\theta })\right),\end{array}$$

(15)

where

$$C={\displaystyle \frac{{\left(\frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}{2}\right)}^{\frac{\nu +np}{2}}}{\mathsf{\Gamma }\left(\frac{\nu +np}{2}\right)T_{(\nu +np)}\left(\Delta (\mathit{Y};\mathit{\theta })\sqrt{\frac{\nu +np}{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}}\right)}}.$$

(16)

Furthermore, by using (7) and (11), it can be shown that

$$\begin{array}{c}\hfill \gamma |\mathit{Y}\sim Tt\left({\displaystyle \frac{\eta (\mathit{Y};\mathit{\theta })}{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}},{\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}{\left(\rho \right(\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1)(\nu +np)}},\nu +np\right)I_{(0,\infty )},\end{array}$$

(17)

where $Tt(\mu ,{\sigma }^2,\nu )I_{\mathbb{A}}$ represents a doubly truncated t distribution with $\nu$ degrees of freedom defined in the interval $\mathbb{A}=\{a_1<x<a_2\}$. From the conditional density in (15), we find

$$\begin{array}{c}\hfill E\left(W\right|\mathit{Y})=C_0{\displaystyle \frac{T_{(\nu +np+2)}\left(\Delta (\mathit{Y};\mathit{\theta })\sqrt{C_2}\right)}{T_{(\nu +np)}\left(\Delta (\mathit{Y};\mathit{\theta })\sqrt{C_0}\right)}},\end{array}$$

(18)

where

$$\begin{array}{c}\hfill C_0={\displaystyle \frac{\nu +np}{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}},C_2={\displaystyle \frac{\nu +np+2}{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}}.\end{array}$$

Additionally, using the law of iterated expectations, we can obtain

$$\begin{array}{c}\hfill E\left(\gamma W\right|\mathit{Y})={\displaystyle \frac{\eta (\mathit{Y};\mathit{\theta })}{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}E\left(W\right|\mathit{Y})+{\displaystyle \frac{1}{\sqrt{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}}\zeta \left(\mathit{Y}\right)\end{array}$$

(19)

and

$$\begin{array}{cc}\hfill E\left({\gamma }^{2}W\right|\mathit{Y})& ={\displaystyle \frac{1}{\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1}}+{\displaystyle \frac{{\eta }^2(\mathit{Y};\mathit{\theta })}{{\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right)}^2}}E\left(W\right|\mathit{Y})\hfill \\ & +{\displaystyle \frac{\eta (\mathit{Y};\mathit{\theta })}{{\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right)}^{3/2}}}\zeta \left(\mathit{Y}\right),\hfill \end{array}$$

(20)

where

$$\begin{array}{cc}\hfill \zeta \left(\mathit{Y}\right)& ={\displaystyle \frac{\mathsf{\Gamma }\left(\frac{\nu +np+1}{2}\right)}{\sqrt{2\pi }{T}_{(\nu +np)}\left(\Delta (\mathit{Y};\mathit{\theta })\sqrt{C_0}\right)\mathsf{\Gamma }\left(\frac{\nu +np}{2}\right)}}{\left({\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu }{2}}\right)}^{-{\displaystyle \frac{\nu +np+1}{2}}}\hfill \\ & \times {\left({\displaystyle \frac{\delta (\mathit{Y};\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu -{\Delta }^2(\mathit{Y};\mathit{\theta })}{2}}\right)}^{{\displaystyle \frac{\nu +np}{2}}}.\hfill \end{array}$$

3.2. Parameter Estimation via the ECME Algorithm

Suppose $\mathit{Y}=({\mathit{Y}}_1,\dots ,{\mathit{Y}}_N)$ constitutes a set of $n\times p$-dimensional observed samples of size N arising from the MVST model. In the EM framework, the latent variables are $\mathit{w}=(w_1,\dots ,w_N)$ and $\mathit{\gamma }=({\gamma }_1,\dots ,{\gamma }_N)$. With these, the complete data are given by ${\mathit{Y}}_c=(\mathit{Y},\mathit{w},\mathit{\gamma })$.

According to (10), the log likelihood function of $\mathit{\theta }$ corresponding to the complete data ${\mathit{Y}}_c$, excluding additive constants and terms that do not involve parameters of the model, is given by

$$\begin{array}{cc}\hfill {\ell }_c(\mathit{\theta }\mid {\mathit{Y}}_c)& ={\displaystyle \frac{1}{2}}\sum _{i=1}^N\{\nu log\left({\displaystyle \frac{\nu }{2}}\right)-2log\mathsf{\Gamma }\left({\displaystyle \frac{\nu }{2}}\right)-plog|\mathbf{\Sigma }|-nlog\left|\mathsf{\Psi }\right|\hfill \\ & +2\eta ({\mathit{Y}}_i;\mathit{\theta }){\gamma }_i{w}_i-\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right){\gamma }_i^2{w}_i\hfill \\ & -\left(\delta ({\mathit{Y}}_i;\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu \right)w_i+\left(\nu +np-1\right)log{w}_i\}.\hfill \end{array}$$

(21)

In the kth iteration, the E step requires the calculation of the so-called Q function, which is the conditional expectation of (21), given the observed data $\mathit{Y}$ and the current estimate ${\widehat{\mathit{\theta }}}^{\left(k\right)}$, where the superscript ^(k) denotes the updated estimates at the kth step of the iterative process. To evaluate the Q function, we then need the following conditional expectations:

$$\begin{array}{c}\hfill {\widehat{w}}_i^{\left(k\right)}=E(W_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}),{\widehat{\kappa }}_{1i}^{\left(k\right)}=E({\gamma }_i{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}),{\widehat{\kappa }}_{2i}^{\left(k\right)}=E({\gamma }_i^2{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}),\end{array}$$

(22)

which have explicit expressions as given earlier, but also the expectation

$${\widehat{\kappa }}_{3i}^{\left(k\right)}=E(log{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}),$$

(23)

which is difficult to evaluate explicitly. So, we perform the ECME algorithm, which replaces the conditional maximization Q function (CMQ) step with the conditional maximization log likelihood (CML) step, to avoid computing the expectation in (23).

Substituting (22) and (23) into (21), we obtain the following expression for the Q function:   

$$\begin{array}{cc}\hfill Q(\mathit{\theta }\mid {\widehat{\mathit{\theta }}}^{\left(k\right)})& ={\displaystyle \frac{1}{2}}\sum _{i=1}^N\{\nu log\left({\displaystyle \frac{\nu }{2}}\right)-2log\mathsf{\Gamma }\left({\displaystyle \frac{\nu }{2}}\right)-plog|\mathbf{\Sigma }|-nlog\left|\mathsf{\Psi }\right|\hfill \\ & +2\eta ({\mathit{Y}}_i;\mathit{\theta }){\widehat{\kappa }}_{1i}^{\left(k\right)}-\left(\rho (\mathbf{\Sigma },\mathsf{\Psi },\mathsf{\Lambda })+1\right){\widehat{\kappa }}_{2i}^{\left(k\right)}\hfill \\ & -\left(\delta ({\mathit{Y}}_i;\mathit{M},\mathbf{\Sigma },\mathsf{\Psi })+\nu \right){\widehat{w}}_i^{\left(k\right)}+\left(\nu +np-1\right){\widehat{\kappa }}_{3i}^{\left(k\right)}\}.\hfill \end{array}$$

(24)

The CMQ steps are implemented, to update estimates of $\mathit{\theta }$ in the order of $\mathit{M}$, $\mathbf{\Sigma }$, $\mathsf{\Psi }$, $\mathsf{\Lambda }$ and $\nu$ by maximizing, one by one, the Q function obtained in the E step. After some algebraic manipulations, they are summarized in the following CMQ and CML steps:

• CMQ step 1: Fixing $\mathsf{\Lambda }={\widehat{\mathsf{\Lambda }}}^{\left(k\right)}$, we update ${\widehat{\mathit{M}}}^{\left(k\right)}$ by maximizing (24) with respect to $\mathit{M}$, leading to

  $${\widehat{\mathit{M}}}^{(k+1)}={\displaystyle \frac{{\sum }_{i=1}^N{\widehat{w}}_i^{\left(k\right)}{\mathit{Y}}_i-{\widehat{\mathsf{\Lambda }}}^{\left(k\right)}{\sum }_{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}}{{\sum }_{i=1}^N{\widehat{w}}_i^{\left(k\right)}}};$$
• CMQ step 2: Fixing $\mathit{M}={\widehat{\mathit{M}}}^{(k+1)}$, $\mathsf{\Psi }={\widehat{\mathsf{\Psi }}}^{\left(k\right)}$ and $\mathsf{\Lambda }={\widehat{\mathsf{\Lambda }}}^{\left(k\right)}$, we then update ${\widehat{\mathbf{\Sigma }}}^{\left(k\right)}$ by maximizing (24) over $\mathbf{\Sigma }$, yielding

  $$\begin{array}{cc}\hfill {\widehat{\mathbf{\Sigma }}}^{(k+1)}& ={\displaystyle \frac{1}{Np}}\{\sum _{i=1}^N{\widehat{w}}_i^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right){\widehat{\mathsf{\Psi }}}^{-1\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)}^{\top }\hfill \\ & +{\widehat{\mathsf{\Lambda }}}^{\left(k\right)}{\widehat{\mathsf{\Psi }}}^{-1\left(k\right)}{\widehat{\mathsf{\Lambda }}}^{\top \left(k\right)}\sum _{i=1}^N{\widehat{\kappa }}_{2i}^{\left(k\right)}-\sum _{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right){\widehat{\mathsf{\Psi }}}^{-1\left(k\right)}{\widehat{\mathsf{\Lambda }}}^{\top \left(k\right)}\hfill \\ & -{\widehat{\mathsf{\Lambda }}}^{\left(k\right)}{\widehat{\mathsf{\Psi }}}^{-1\left(k\right)}\sum _{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)}^{\top }\};\hfill \end{array}$$
• CMQ step 3: Fixing $\mathit{M}={\widehat{\mathit{M}}}^{(k+1)}$, $\mathbf{\Sigma }={\widehat{\mathbf{\Sigma }}}^{(k+1)}$ and $\mathsf{\Lambda }={\widehat{\mathsf{\Lambda }}}^{\left(k\right)}$, we update ${\widehat{\mathsf{\Psi }}}^{\left(k\right)}$ by maximizing (24) over $\mathsf{\Psi }$, yielding

  $$\begin{array}{cc}\hfill {\widehat{\mathsf{\Psi }}}^{(k+1)}& ={\displaystyle \frac{1}{Nn}}\{\sum _{i=1}^N{\widehat{w}}_i^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)}^{\top }{\widehat{\mathbf{\Sigma }}}^{-1(k+1)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)\hfill \\ & +{\widehat{\mathsf{\Lambda }}}^{\top \left(k\right)}{\widehat{\mathbf{\Sigma }}}^{-1(k+1)}{\widehat{\mathsf{\Lambda }}}^{\left(k\right)}\sum _{i=1}^N{\widehat{\kappa }}_{2i}^{\left(k\right)}-\sum _{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)}^{\top }{\widehat{\mathbf{\Sigma }}}^{-1(k+1)}{\widehat{\mathsf{\Lambda }}}^{\left(k\right)}\hfill \\ & -{\widehat{\mathsf{\Lambda }}}^{\top \left(k\right)}{\widehat{\mathbf{\Sigma }}}^{-1(k+1)}\sum _{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)\};\hfill \end{array}$$
• CMQ step 4: Fixing $\mathit{M}={\widehat{\mathit{M}}}^{(k+1)}$, we obtain ${\widehat{\mathsf{\Lambda }}}^{(k+1)}$ by maximizing (24) over $\mathsf{\Lambda }$, yielding

  $$\begin{array}{c}\hfill {\widehat{\mathsf{\Lambda }}}^{(k+1)}={\displaystyle \frac{{\sum }_{i=1}^N{\widehat{\kappa }}_{1i}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}^{(k+1)}\right)}{{\sum }_{i=1}^N{\widehat{\kappa }}_{2i}^{\left(k\right)}}}.\end{array}$$

An update of ${\widehat{\nu }}^{\left(k\right)}$ can be achieved by directly maximizing the constrained actual log likelihood function. This gives rise to the following CML step:

• CML step: Update ${\widehat{\nu }}^{\left(k\right)}$ by optimizing the following constrained log likelihood function:

  $${\widehat{\nu }}^{(k+1)}=arg\underset{\nu }{max}\sum _{i=1}^{N}log{f}_{MVST}\left({\mathit{Y}}_i;{\widehat{\mathit{M}}}^{(k+1)},{\widehat{\mathbf{\Sigma }}}^{(k+1)},{\widehat{\mathsf{\Psi }}}^{(k+1)},{\widehat{\mathsf{\Lambda }}}^{(k+1)},\nu \right).$$

4. Fitting Finite Mixtures of MVST Distributions

4.1. The Model

We consider N independent random variables ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_N$ observed from a G-component mixture of MVST distributions, whose pdf is given by

$$\begin{array}{c}\hfill f({\mathit{Y}}_i;\mathsf{\Theta })=\sum _{g=1}^G{\pi }_g{f}_{MVST}({\mathit{Y}}_i;{\mathit{\theta }}_g),\end{array}$$

where $0\le {\pi }_g\le 1$, ${\sum }_{i=1}^G{\pi }_g=1$ and $\mathsf{\Theta }$ is the set containing all the parameters of the considered mixture model. To pose this mixture model as an incomplete data problem, we introduce allocation variables ${\mathit{Z}}_i=(Z_{i1},\dots ,Z_{iG})$, where a particular element $Z_{ig}$ is equal to 1 if $Y_i$ belongs to group g and is equal to zero, otherwise. Observe that ${\mathit{Z}}_i$ follows a multinomial random vector with one trial and cell probabilities ${\pi }_1,\dots ,{\pi }_G$, denoted by ${\mathit{Z}}_i\sim \mathcal{M}(1;{\pi }_1,\dots ,{\pi }_G)$. The hierarchical representation in (9), originally designed for the single distribution, can be extended to the mixture modeling framework, as follows:

$$\begin{array}{cc}\hfill {\mathit{Y}}_i|({\gamma }_i,w_i,Z_{ig}=1)& \sim N_{n\times p}({\mathit{M}}_g+{\gamma }_i{\mathsf{\Lambda }}_g,w_i^{-1}{\mathbf{\Sigma }}_g,{\mathsf{\Psi }}_g),\hfill \\ \hfill {\gamma }_i|(w_i,Z_{ig}=1)& \sim TN(0,w_i^{-1})I_{(0,\infty )},\hfill \\ \hfill W_i|Z_{ig}=1& \sim \mathsf{\Gamma }({\nu }_g/2,{\nu }_g/2),\hfill \\ \hfill Z_{ig}& \sim \mathcal{M}(1;{\pi }_1,\dots ,{\pi }_G).\hfill \end{array}$$

(25)

It then follows from the hierarchical structure in (25), on the basis of the observed data $\mathit{Y}=({\mathit{Y}}_1,\dots ,{\mathit{Y}}_N)$ and the latent data $\mathit{w}=(w_1,\dots ,w_N)$, $\mathit{\gamma }=({\gamma }_1,\dots ,{\gamma }_N)$ and $\mathit{Z}=(Z_1,\dots ,Z_N)$, excluding additive constants, the complete data log likelihood function of $\mathsf{\Theta }$ based on the complete data ${\mathit{Y}}_c=(\mathit{Y},\mathit{w},\mathit{\gamma },\mathit{Z})$ is

$$\begin{array}{cc}\hfill {\ell }_c(\mathit{\theta }\mid {\mathit{Y}}_c)& ={\displaystyle \frac{1}{2}}\sum _{i=1}^N\sum _{g=1}^G{Z}_{ig}\{2log{\pi }_g+{\nu }_{g}log\left({\displaystyle \frac{{\nu }_g}{2}}\right)-2log\mathsf{\Gamma }\left({\displaystyle \frac{{\nu }_g}{2}}\right)-plog\left|{\mathbf{\Sigma }}_g\right|\hfill \\ & -nlog|{\mathsf{\Psi }}_g|+2\eta ({\mathit{Y}}_i;{\mathit{\theta }}_g){\gamma }_i{w}_i-\left(\rho ({\mathbf{\Sigma }}_g,{\mathsf{\Psi }}_g,{\mathsf{\Lambda }}_g)+1\right){\gamma }_i^2{w}_i\hfill \\ & -\left(\delta ({\mathit{Y}}_i;{\mathit{M}}_g,{\mathbf{\Sigma }}_g,{\mathsf{\Psi }}_g)+{\nu }_g\right)w_i+\left({\nu }_g+np-1\right)log{w}_i\}.\hfill \end{array}$$

(26)

The expected value of (26) to start the E step, given the current parameter ${\mathsf{\Theta }}^{\left(k\right)}$, requires some conditional expectations, including

$$\begin{array}{cc}\hfill {\widehat{z}}_{ig}^{\left(k\right)}& =E\left(Z_{ig}\right|{\mathit{Y}}_i,{\widehat{\mathsf{\Theta }}}^{\left(k\right)})={\displaystyle \frac{{\widehat{\pi }}_g^{\left(k\right)}{f}_{MVST}({\mathit{Y}}_i;{\widehat{\mathit{\theta }}}_g^{\left(k\right)})}{f({\mathit{Y}}_i;{\widehat{\mathsf{\Theta }}}^{\left(k\right)})}},\hfill \\ \hfill {\widehat{w}}_{ig}^{\left(k\right)}& =E(W_i\mid {\mathit{Y}}_i,Z_{ig}=1,{\widehat{\mathsf{\Theta }}}^{\left(k\right)})=E(W_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}_g^{\left(k\right)}),\hfill \\ \hfill {\widehat{\kappa }}_{1ig}^{\left(k\right)}& =E({\gamma }_i{W}_i\mid {\mathit{Y}}_i,Z_{ig}=1,{\widehat{\mathsf{\Theta }}}^{\left(k\right)})=E({\gamma }_i{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}_g^{\left(k\right)}),\hfill \\ \hfill {\widehat{\kappa }}_{2ig}^{\left(k\right)}& =E({\gamma }_i^2{W}_i\mid {\mathit{Y}}_i,Z_{ig}=1,{\widehat{\mathsf{\Theta }}}^{\left(k\right)})=E({\gamma }_i^2{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}_g^{\left(k\right)})\hfill \end{array}$$

(27)

and

$${\widehat{\kappa }}_{3ig}^{\left(k\right)}=E(log{W}_i\mid {\mathit{Y}}_i,Z_{ig}=1,{\widehat{\mathsf{\Theta }}}^{\left(k\right)})=E(log{W}_i\mid {\mathit{Y}}_i,{\widehat{\mathit{\theta }}}_g^{\left(k\right)})$$

for which we utilize the CML step, as mentioned in the preceding section. Consequently, the conditional expectation of the complete data log likelihood is obtained as   

$$\begin{array}{cc}\hfill Q(\mathit{\theta }\mid {\widehat{\mathit{\theta }}}^{\left(k\right)})& ={\displaystyle \frac{1}{2}}\sum _{i=1}^N\sum _{g=1}^G{\widehat{z}}_{ig}^{\left(k\right)}\{2log{\pi }_g+{\nu }_{g}log\left({\displaystyle \frac{{\nu }_g}{2}}\right)-2log\mathsf{\Gamma }\left({\displaystyle \frac{{\nu }_g}{2}}\right)-plog\left|{\mathbf{\Sigma }}_g\right|\hfill \\ & -nlog|{\mathsf{\Psi }}_g|+2\eta ({\mathit{Y}}_i;{\mathit{\theta }}_g){\widehat{\kappa }}_{1ig}^{\left(k\right)}-\left(\rho ({\mathbf{\Sigma }}_g,{\mathsf{\Psi }}_g,{\mathsf{\Lambda }}_g)+1\right){\widehat{\kappa }}_{2ig}^{\left(k\right)}\hfill \\ & -\left(\delta ({\mathit{Y}}_i;{\mathit{M}}_g,{\mathbf{\Sigma }}_g,{\mathsf{\Psi }}_g)+{\nu }_g\right){\widehat{w}}_{ig}^{\left(k\right)}+\left({\nu }_g+np-1\right){\widehat{\kappa }}_{3ig}^{\left(k\right)}\}.\hfill \end{array}$$

(28)

Thus, the implementation of the ECME algorithm proceeds as follows:

• E step: Given $\mathsf{\Theta }$$={\widehat{\mathsf{\Theta }}}^{\left(k\right)}$, compute ${\widehat{z}}_{ig}^{\left(k\right)}$, ${\widehat{w}}_{ig}^{\left(k\right)}$, ${\widehat{\kappa }}_{1ig}^{\left(k\right)}$ and ${\widehat{\kappa }}_{2ig}^{\left(k\right)}$ given in (27), for $i=1,\dots ,N$ and $g=1,\dots ,G$;
• CM step 1: Calculate

  $${\widehat{\pi }}_g^{(k+1)}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)};$$
• CM step 2: Update ${\widehat{\mathit{M}}}_g^{\left(k\right)}$ as

  $$\begin{array}{c}\hfill {\widehat{\mathit{M}}}_g^{(k+1)}={\displaystyle \frac{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{w}}_{ig}^{\left(k\right)}{\mathit{Y}}_i-{\widehat{{\mathsf{\Lambda }}_g}}^{\left(k\right)}{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}}{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{w}}_{ig}^{\left(k\right)}}};\end{array}$$
• CM step 3: Update ${\widehat{\mathbf{\Sigma }}}_g^{\left(k\right)}$ as

  $$\begin{array}{cc}\hfill {\widehat{\mathbf{\Sigma }}}_g^{(k+1)}& ={\displaystyle \frac{1}{p{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}}}\{\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{w}}_{ig}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right){\widehat{\mathsf{\Psi }}}_g^{-1\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)}^{\top }\hfill \\ & +{\widehat{\mathsf{\Lambda }}}_g^{\left(k\right)}{\widehat{\mathsf{\Psi }}}_g^{-1\left(k\right)}{\widehat{\mathsf{\Lambda }}}_g^{\top \left(k\right)}\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{2ig}^{\left(k\right)}-\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right){\widehat{\mathsf{\Psi }}}_g^{-1\left(k\right)}{\widehat{\mathsf{\Lambda }}}_g^{\top \left(k\right)}\hfill \\ & -{\widehat{\mathsf{\Lambda }}}_g^{\left(k\right)}{\widehat{\mathsf{\Psi }}}_g^{-1\left(k\right)}\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)}^{\top }\};\hfill \end{array}$$
• CM step 4: Update ${\widehat{\mathsf{\Psi }}}_g^{\left(k\right)}$ as

  $$\begin{array}{cc}\hfill {\widehat{\mathsf{\Psi }}}_g^{(k+1)}& ={\displaystyle \frac{1}{n{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}}}\{\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{w}}_{ig}^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)}^{\top }{\widehat{\mathbf{\Sigma }}}_g^{-1(k+1)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)\hfill \\ & +{\widehat{\mathsf{\Lambda }}}_g^{\top \left(k\right)}{\widehat{\mathbf{\Sigma }}}_g^{-1(k+1)}{\widehat{\mathsf{\Lambda }}}_g^{\left(k\right)}\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{2ig}^{\left(k\right)}-\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}{\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)}^{\top }{\widehat{\mathbf{\Sigma }}}_g^{-1(k+1)}{\widehat{\mathsf{\Lambda }}}_g^{\left(k\right)}\hfill \\ & -{\widehat{\mathsf{\Lambda }}}_g^{\top \left(k\right)}{\widehat{\mathbf{\Sigma }}}_g^{-1(k+1)}\sum _{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)\};\hfill \end{array}$$
• CM step 5: Update ${\widehat{\mathsf{\Lambda }}}_g^{\left(k\right)}$ as

  $$\begin{array}{c}\hfill {\widehat{\mathsf{\Lambda }}}_g^{(k+1)}={\displaystyle \frac{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{1ig}^{\left(k\right)}\left({\mathit{Y}}_i-{\widehat{\mathit{M}}}_g^{(k+1)}\right)}{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(k\right)}{\widehat{\kappa }}_{2ig}^{\left(k\right)}}};\end{array}$$
• CML step: Update ${\widehat{\mathbf{\nu }}}^{\left(k\right)}=({\widehat{\nu }}_1,\dots ,{\widehat{\nu }}_G)$ by optimizing the constrained log likelihood function as 

  $${\widehat{\mathbf{\nu }}}^{(k+1)}=arg\underset{\mathbf{\nu }}{max}\sum _{i=1}^{N}log\left(\sum _{g=1}^G{\widehat{\pi }}_g^{(k+1)}{f}_{MVST}\left({\mathit{Y}}_i;{\widehat{\mathit{M}}}_g^{(k+1)},{\widehat{\mathbf{\Sigma }}}_g^{(k+1)},{\widehat{\mathsf{\Psi }}}_g^{(k+1)},{\widehat{\mathsf{\Lambda }}}_g^{(k+1)},{\nu }_g\right)\right).$$

4.2. Initialization

In order to speed up the convergence process, it is important to establish a set of reasonable starting values. To start the ECME algorithm for fitting the FM-MVST model, an intuitive scheme to partition data into G components ${\left\{{\mathit{Y}}_g^{\left(0\right)}\right\}}_{g=1}^G$ is to create an initial partition of data ${\left\{\mathrm{Vec}\left({\mathit{Y}}_i\right)\right\}}_{i=1}^N$, using the K-means algorithm [21,22]. This yields a validate estimate of ${\widehat{z}}_{ig}^{\left(0\right)}$, which, in turn, yields ${\widehat{\pi }}_g^{\left(0\right)}={\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(0\right)}/N$. Then, we compute the sample mean, covariance matrix of rows and covariance matrix of columns of ${\mathit{Y}}_g^{\left(0\right)}$ as good initial estimates for ${\widehat{\mathit{M}}}_g^{\left(0\right)}$, ${\widehat{\mathbf{\Sigma }}}_g^{\left(0\right)}$ and ${\widehat{\mathsf{\Psi }}}_g^{\left(0\right)}$, as follows:

$$\begin{array}{cc}\hfill {\widehat{\mathit{M}}}_g^{\left(0\right)}& ={\displaystyle \frac{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(0\right)}{\mathit{Y}}_i}{{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(0\right)}}},\hfill \\ \hfill {\widehat{\mathbf{\Sigma }}}_g^{\left(0\right)}& ={\displaystyle \frac{{\sum }_{i=1}^N{\sum }_{j=1}^p{\widehat{z}}_{ig}^{\left(0\right)}\left({\mathit{y}}_{ij}-{\mathit{m}}_{gj}^{\left(0\right)}\right){\left({\mathit{y}}_{ij}-{\mathit{m}}_{gj}^{\left(0\right)}\right)}^{\top }}{p{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(0\right)}}},\hfill \\ \hfill {\widehat{\mathsf{\Psi }}}_g^{\left(0\right)}& ={\displaystyle \frac{{\sum }_{i=1}^N{\sum }_{r=1}^n{\widehat{z}}_{ig}^{\left(0\right)}{\left({\mathit{y}}_{i.r}-{\mathit{m}}_{g.r}^{\left(0\right)}\right)}^{\top }\left({\mathit{y}}_{i.r}-{\mathit{m}}_{g.r}^{\left(0\right)}\right)}{n{\sum }_{i=1}^N{\widehat{z}}_{ig}^{\left(0\right)}}},\hfill \end{array}$$

where ${\mathit{m}}_{gj}^{\left(0\right)}$ and ${\mathit{m}}_{g.r}^{\left(0\right)}$ denote the j-th column and r-th row of ${\mathit{M}}_g^{\left(0\right)}$, respectively, and ${\mathit{y}}_{ij}$ and ${\mathit{y}}_{i.r}$ are the j-th column and r-th row of ${\mathit{Y}}_i$, respectively. The initial component skewness matrices, ${\mathsf{\Lambda }}_g^{\left(0\right)}$, are taken as the values randomly selected in the interval $(-1,1)$. Finally, we initialize ${\widehat{\mathbf{\nu }}}_g^{\left(0\right)}$ by taking it to be as small as 5 or 10.

4.3. Identifiability

Model identifiability is key to securing unique and consistent estimates of model parameters. With regard to the mixtures of MVST distributions, the estimates of ${\mathbf{\Sigma }}_g$ and ${\mathsf{\Psi }}_g$ are only unique up to a strictly positive constant. To resolve this issue, a constraint needs to be placed, such as setting the trace of ${\mathbf{\Sigma }}_g$ equal to n [13] or fixing $|{\mathbf{\Sigma }}_g|=1$ [23]. Herein, we set the first diagonal element of ${\mathbf{\Sigma }}_g$ as 1 [8]. This scaling procedure can be implemented at each iteration or at convergence, and either method has minimal impact on the final estimates and classifications achieved. To obtain the final parameter estimates, the resulting ${\mathbf{\Sigma }}_g$ is divided by the first diagonal element of ${\mathbf{\Sigma }}_g$, and then ${\mathsf{\Psi }}_g$ is multiplied by the first diagonal element of ${\mathbf{\Sigma }}_g$.

5. Empirical Study

5.1. Finite-Sample Properties of ML Estimators

Here, we conducted a simulation study for examining the accuracy of the parameter estimates obtained by using the proposed ECME algorithm in Section 3.2. We generated 500 Monte Carlo samples of sizes $N=250$, 500, 1000 and 2000 from the two-component FM-MVST model, in the two scenarios (low- and moderate-dimensional) described in Appendix A. Scenario I was characterized by matrices of size $3\times 4$, indicating that each matrix consisted of 3 rows and 4 columns. This structure allowed for a total of 12 elements within each matrix. On the other hand, Scenario II generated matrices of size $10\times 2$. In this case, each matrix was composed of 10 rows and 2 columns, resulting in a total of 20 elements per matrix. This increased number of elements was beneficial for examining the performance of the proposed algorithm in accurately recovering true parameters in moderate dimensional scenarios.

The accuracy of the obtained parameter estimates was assessed by the average of the root mean squared error (RMSE) of the elements of each estimated parameter. The results shown in

Table 1. Average RMSE based on 500 replications for the evaluation of ML estimates.

Scenario	N	${\mathit{\pi }}_1$	${\mathit{M}}_1$	${\mathit{M}}_2$	${\mathbf{\Sigma }}_1$	${\mathbf{\Sigma }}_2$	${\mathbf{\Psi }}_1$	${\mathbf{\Psi }}_2$	${\mathbf{\Lambda }}_1$	${\mathbf{\Lambda }}_2$	${\mathit{\nu }}_1$	${\mathit{\nu }}_2$
I	250	0.031	0.169	0.119	0.058	0.045	0.128	0.078	0.174	0.118	1.638	3.629
	500	0.020	0.129	0.087	0.044	0.041	0.093	0.064	0.126	0.088	1.215	3.128
	1000	0.015	0.091	0.059	0.032	0.031	0.066	0.052	0.088	0.056	0.803	2.712
	2000	0.011	0.062	0.044	0.027	0.028	0.051	0.043	0.062	0.044	0.702	2.230
II	250	0.035	0.260	0.193	1.348	0.376	0.785	0.440	0.257	0.187	2.544	2.426
	500	0.022	0.176	0.130	1.257	0.338	0.771	0.425	0.179	0.129	2.297	2.219
	1000	0.014	0.120	0.097	1.206	0.312	0.764	0.423	0.124	0.097	1.962	1.897
	2000	0.011	0.088	0.068	1.184	0.302	0.753	0.401	0.090	0.070	1.460	1.388

indicate the good performance of the proposed estimation method. Regardless of the considered scenario, it can be seen that the RMSE values all tended to zero with increasing sample size, indicating the satisfactory asymptotic properties of the ML estimates obtained by the proposed ECME algorithm.

5.2. Comparison of Classification Accuracy

To examine the classification accuracy of the FM-MVST model, we generated 1000 samples from each of the scenarios given in Appendix A. In each scenario, we compared the FM-MVST model described in Section 4 with finite mixtures of matrix variate normal (FM-MVN) and matrix variate t (FM-MVT) distributions, which were readily available in the R package MixMatrix. We also implemented the EM algorithm described in [8], to fit finite mixtures of MVSTIG (FM-MVSTIG) distributions. Additionally, the finite mixtures of the reduced RMVSN (FM-RMVSN) distributions were fitted as a sub-model of the FM-MVST model.

Model performance was assessed by comparing the classification accuracy and model selection criteria for all the fitted models. For classification accuracy, we report the adjusted rand index (ARI; [24]), which took the value of 1 when the two partitions perfectly matched, and the misclassification rate (MCR) of the map clustering for each model. ARI serves as a measure of the similarity between two data clusters, which can provide insights into the robustness of clustering results. By contrast, MCR focuses on the accuracy of classification, emphasizing the proportion of misclassified instances, thus offering a more direct assessment of predictive performance. Furthermore, the Bayesian information criterion (BIC; [25]) value was also reported as a model selection criterion. BIC incorporates both the goodness-of-fit and the complexity of the model, penalizing for overfitting, which can be particularly relevant when evaluating clustering algorithms.

We ran 100 simulations for each scenario and computed the classification accuracy and model selection criteria for each simulation.

Table 2. Simulation results, based on 100 replications, for performance comparison of four mixture models in two scenarios.

Scenario	Model	BIC	Std	ARI	Std	MCR	Std
	FM-MVN	55,419.78	831.20	0.82	0.20	0.05	0.06
	FM-MVT	45,004.95	787.39	0.90	0.17	0.03	0.05
I	FM-RMVSN	40,103.90	962.33	0.95	0.18	0.02	0.05
	FM-MVSTIG	38,215.01	859.47	0.97	0.08	0.01	0.02
	FM-MVST	38,170.98	804.07	0.98	0.05	0.01	0.01
	FM-MVN	85,450.81	705.80	0.91	0.15	0.07	0.06
	FM-MVT	76,011.44	720.48	0.93	0.14	0.08	0.05
II	FM-RMVSN	72,917.04	703.48	0.94	0.12	0.05	0.02
	FM-MVSTIG	69,892.52	694.08	0.95	0.10	0.04	0.02
	FM-MVST	69,839.90	673.31	0.97	0.09	0.02	0.01

presents the average BIC, ARI and MCR values along with their standard errors (Std), and the results are illustrated via the box plots shown in Figure 1. As one would expect, the model selection criteria selected the true model from which the data were generated. This outcome highlights the effectiveness of the selected metrics in distinguishing between models based on their ability to capture the underlying data structure. The consistency observed across the simulations further strengthens the case for the reliability of these model selection criteria in practical applications.

6. Real Data Analysis

In this section, we illustrate the results of applying the proposed methodology to three well-known real datasets.

6.1. Landsat Data

The first application concerned the Landsat data (LSD), originally obtained by NASA, and available at Irvine machine learning repository (http://archive.ics.uci.edu/ml, accessed on 1 September 2024). Multi-spectral satellite imagery allows for multiple observations over a spatial grid, resulting in matrix-valued observations. The LSD comprises lines that consist of four spectral values representing nine pixel neighborhoods in a satellite image. Essentially, each line corresponds to a $4\times 9$ observation matrix. Additionally, every observation matrix in the LSD is classified into one of six distinct categories: red soil, cotton crop, gray soil, damp gray soil, soil with vegetation stubble and very damp gray soil. For our analysis, we concentrated on three specific categories: red soil, gray soil and soil with vegetation stubble, which had sizes of 461, 397 and 237, respectively.

Table 3. Summary results from fitting various models to the LSD data.

Model	G	Log Likelihood	BIC	ARI	MCR
FM-MVN		−114,954.90	231,799.40	0.67	0.14
FM-MVT		−113,169.30	228,228.10	0.69	0.13
FM-RMVSN	3	−111,213.50	225,107.40	0.76	0.09
FM-MVSTIG		−110,920.90	224,543.20	0.79	0.07
FM-MVST		−110,836.60	224,374.60	0.82	0.06

presents a summary of the ML fitting results, including the maximized log likelihood values, BIC, ARI and MCR of the four fitted models. The results reveal that the log likelihood value for the FM-MVN distribution was lower than that for the FM-MVT distribution, indicating a poorer fit. In contrast, the skewed distributions (FM-MVST, FM-MVSTIG and FM-RMVSN) outperformed their respective models. Particularly noteworthy was the superior performance of the FM-MVST model. The estimated tailedness parameters were ${\widehat{\nu }}_1=0.47$, ${\widehat{\nu }}_2=0.44$ and ${\widehat{\nu }}_3=0.58$, indicating a distribution of matrix observations characterized by long-tailed behavior.

6.2. Apes Data

The second application considered the apes dataset included in the shapes R package [26]. The description of the dataset, taken from [27], is as follows. In an investigation to assess the cranial differences between the sexes of apes, 29 male and 30 female adult gorillas (Gorilla), 28 male and 26 female adult chimpanzees (Pan) and 30 male and 24 female adult orangutans (Pongo) were studied. Eight landmarks were chosen in the midline plane of each skull. These landmarks were anatomical landmarks and were located by an expert biologist. The dataset was stored as a list with two components: an array of coordinates in eight landmarks in two dimensions for each skull ($8\times 2$ observation matrix and $N=167$), and a vector of group labels ($G=6$).

All the competing models were fitted for $G=6$, and their fitting results are reported in

Table 4. Summary results from fitting various models to the apes data.

Model	G	Log Likelihood	BIC	ARI	MCR
FM-MVN		−7773.14	17,204.51	0.51	0.41
FM-MVT		−7609.66	16,877.54	0.56	0.32
FM-RMVSN	6	−6158.09	14,522.03	0.60	0.28
FM-MVSTIG		−6097.66	14,431.87	0.63	0.27
FM-MVST		−5970.42	14,177.41	0.67	0.25

. It is clear that the FM-MNV model provided the worst fitting performance, whereas FM-MVST was the best model. Similarly to the analyses in the previous section, this may be an indication that the components of FM-MVN were not skewed and heavy-tailed enough to adequately model the data. On a related note, the estimated tailedness parameters were ${\widehat{\nu }}_1=0.45$, ${\widehat{\nu }}_2=2.20$, ${\widehat{\nu }}_3=1.61$, ${\widehat{\nu }}_4=2.80$, ${\widehat{\nu }}_5=0.52$ and ${\widehat{\nu }}_6=1.61$, highlighting the presence of clusters with high levels of tailed behavior.

6.3. Melanoma Data

The performance of the FM-MVST model in skin cancer detection was demonstrated in the third and final application. The objective of the skin cancer detection project was to develop a framework for analyzing and assessing the risk of melanoma, using dermatological photographs taken with a standard consumer-grade camera. Segmentation of the lesion is a crucial step for developing a skin cancer detection framework. The objective, then, was to find the border of the skin lesion. It was important that this step was performed accurately, because many features used to assess the risk of melanoma are derived based on the lesion border. The set of images included images extracted from the public databases DermIS and DermQuest, along with manual segmentations (ground truth) of the lesions, available at https://uwaterloo.ca/vision-image-processing-lab/research-demos/skin-cancer-detection, accessed on 1 September 2024.

A skin image in $100\times 70$ pixels is displayed in Figure 2a. The next objective was to segment the image into two labels. We considered all the pixels of three numerical RGB components denoting red, green and blue intensities and a grayscale intensity, such as $y_i\in {[0,255]}^4$, which could be transformed into ${[0,1]}^4$. Upon considering each pixel as a $2\times 2$ matrix, each pixel was then be grouped into $G=2$ clusters, where every cluster was assumed to have a different distribution.

It follows from

Table 5. Summary results from fitting various models to the melanoma data.

Model	G	Log Likelihood	BIC	ARI	MCR
FM-MVN		63,941.64	−127,723.90	0.76	0.14
FM-MVT		65,509.39	−130,859.40	0.82	0.13
FM-RMVSN	2	65,601.15	−130,945.50	0.91	0.12
FM-MVSTIG		65,788.76	−131,303.10	0.93	0.11
FM-MVST		67,241.98	−134,209.50	0.95	0.09

that the FM-MVST model provided the best fit, in terms of BIC, as well as the lowest misclassification error for the binary classification of each pixel. The estimates of the tailedness parameters were ${\widehat{\nu }}_1=1.96$ and ${\widehat{\nu }}_2=0.92$, signifying the appropriateness of the use of heavy-tailed t distributions. Furthermore, the superiority of the FM-MVST model is reflected visually in Figure 2c–g, which depict the comparative segmentation performance of the fitted models in grayscale. The figures illustrate differences in identifying the lesion area, and they indicate that the proposed model exhibits a clearer boundary and a more consistent region of the lesion.

It is noteworthy that all the model selection criteria applied to the three datasets strongly favored the proposed FM-MVST model. The datasets varied significantly, in terms of dimensional and structural characteristics, which underscores the flexibility and effectiveness of the FM-MVST model. This model demonstrated a superior ability to accurately capture the skewness and leptokurtic features present in the data, outperforming the alternative models. The adaptability of the FM-MVST model across diverse datasets not only showcases its robustness but also reinforces its potential as a valuable tool for data analysis in various applications.

7. Concluding Remarks

We have introduced here a new family of matrix variate distributions that can capture both skewness and heavy-tailedness simultaneously. This MVST model was based on a stochastic representation that facilitated our developing an ECME algorithm for the maximum likelihood estimation of the model parameters. We evaluated the effectiveness and efficiency of the proposed algorithm through two simulation studies. Additionally, we used the proposed approach to analyze three real data datasets, demonstrating its capability in modeling asymmetric matrix variate data. Future developments of this approach could include accommodating censored data, determining the optimal number of mixing distributions for the clustering problem and including different distributions for the variables W and U in the stochastic representation. Another interesting extension could involve incorporating the FM-MVST distribution into a mixture-of-regression framework. We are currently looking into these problems, and we hope to report our findings in a future paper.