A Matrix Variate Skew-t Distribution
arXiv:1703.01364v3 [stat.ME] 13 Apr 2017

Michael P.B. Gallaugher and Paul D. McNicholas
Dept. of Mathematics & Statistics, McMaster University, Hamilton, Ontario, Canada.

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.
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
2.1

Background
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
1

matrix-normal distribution. For notional clarity, we use X to denote a realization of a random
matrix X . An n × p random matrix X follows a matrix variate normal distribution with
location parameter M and scale matrices Σ and Ψ of dimensions n×n and p×p, respectively.
We write X ∼ Nn×p (M, Σ, Ψ) to denote such a random matrix and the density of X can
be written



1
1
−1
−1
0
. (1)
f (X | M, Σ, Ψ) =
− tr Σ (X − M)Ψ (X − M)
np
p
n exp
2
(2π) 2 |Σ| 2 |Ψ| 2
One well known property of the matrix variate normal distribution (Harrar and Gupta,
2008) is
X ∼ Nn×p (M, Σ, Ψ) ⇐⇒ vec(X ) ∼ Nnp (vec(M), Ψ ⊗ Σ)
(2)
where Nnp (·) is the multivariate normal density with dimension np, vec(M) is the vectorization of M, and ⊗ 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, Doğ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 X takes the form
Z ∞
f (x) =
φp (x|µ + wα, wΣ)h(w|θ)dw,
0

which is equivalent to the representation

√

X = µ + Wα +

W V,

(3)

where V ∼ Np (0, Σ) and W > 0 is a latent random variable with density h(w|θ). The
multivariate skew-t
distribution with ν degrees of freedom arises as a special case with

ν ν
W ∼ IG 2 , 2 , where IG(·) denotes the inverse Gamma distribution with density function


β α −α−1
β
f (x|α, β) =
x
exp −
.
Γ(α)
x
2

2.3

The Generalized Inverse Gaussian Distribution

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

(
)
ay + yb
y λ−1
√
f (y|a, b, λ) =
exp −
,
2
2Kλ ( ab)
a
b

where

1
Kλ (x) =
2

Z

∞

y
0

λ−1




x
1
exp −
y+
dy
2
y

is the modified Bessel function of the third kind with index λ. Several functions of GIG
random variables have tractable expected values, e.g.,
r
√
b Kλ+1 ( ab)
√
E(Y ) =
,
(4)
a Kλ ( ab)
√
  r
1
a Kλ+1 ( ab) 2λ
√
E
=
− ,
(5)
Y
b Kλ ( ab)
b
r !
√
1
∂
b
√
Kλ ( ab),
E(log Y ) = log
+
(6)
a
Kλ ( ab) ∂λ
where a, b ∈ R+ , λ ∈ R, and Kλ (·) is the modified Bessel function of the third kind with
index λ. These results will prove to be useful for parameter estimation for the matrix-variate
skew-t distribution.

3
3.1

Methodology
A Matrix Variate Skew-t Distribution

We will say that an n × p random matrix X has a matrix variate skew-t distribution,
MVSTn×p (M, A, Σ, Ψ, ν), if X can be written
√
(7)
X = M + WA + WV ,

where M and A are n × p matrices, V ∼ Nn×p (0, Σ, Ψ), and W ∼ IG ν2 , ν2 . Analogous to
its multivariate counterpart, M is a location matrix, A is a skewness matrix, Σ and Ψ are
scale matrices, and ν is the degrees of freedom. It then follows that
X |w ∼ Nn×p (M + wA, wΣ, Ψ)

3

and thus the joint density of X and W is
f (X, w | ϑ) = f (X | w)f (w)
=

ν
2
np
2

ν
2
p
2

n
2

Γ( ν2 )

w−

ν+np
−1
2

(2π) |Σ| |Ψ|




1 
−1
−1
0
tr Σ (X − M − wA)Ψ (X − M − wA) + ν ,
× exp −
2w

(8)

where ϑ = (M, A, Σ, Ψ, ν). We note that the exponential term in (8) can be written as




1 δ(X; M, Σ, Ψ) + ν
−1
−1 0
+ wρ(A, Σ, Ψ) ,
exp tr(Σ (X − M)Ψ A ) × exp −
2
w
where
δ(X; M, Σ, Ψ) = tr(Σ−1 (X − M)Ψ−1 (X − M)0 ) and ρ(A, Σ, Ψ) = tr(Σ−1 AΨ−1 A0 ).
Therefore, the marginal density of X is
Z ∞
f (X) =
f (X, w)dw
0

=

ν
2

ν
2
p
2

np
2


exp tr(Σ−1 (X − M)Ψ−1 A0 )

n
2

(2π) |Σ| |Ψ| Γ( ν2 )



Z ∞
1 δ(X; M, Σ, Ψ) + ν
− ν+np
−1
w 2
exp −
×
+ wρ(A, Σ, Ψ) dw.
2
w
0

Making the change of variables given by
p
ρ(A, Σ, Ψ)

y=p
w
δ(X; M, Σ, Ψ) + ν
we can write
fMVST (X | ϑ) =

2

ν
2

 ν2

exp { tr(Σ−1 (X − M)Ψ−1 A0 )}
np

p

n



δ(X; M, Σ, Ψ) + ν
ρ(A, Σ, Ψ)

− ν+np
4

(2π) 2 |Σ| 2 |Ψ| 2 Γ( ν2 )
p

× K− ν+np
[ρ(A, Σ, Ψ)] [δ(X; M, Σ, Ψ) + ν] .
2

The density of 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 Σ and Ψ are
unique only up to a multiplicative constant. Indeed, if we let Σ̃ = vΣ and Ψ̃ = (1/v)Ψ,
v 6= 0, the likelihood is unchanged. This identifiability issue can be resolved, for example,
4

by setting tr(Σ) = n or tr(Ψ) = p. Note that Ψ̃ ⊗ Σ̃ = Ψ ⊗ Σ, so the estimate of the
Kronecker product is unique.
For the purposes of parameter estimation, note that the conditional density of W is
f (w | X) =

f (X | w)f (w)
f (X)

λ


[ρ(A, Σ, Ψ)/(δ(X; M, Σ, Ψ) + ν)] 2 wλ−1
ρ(A, Σ, Ψ)w + [δ(X; M, Σ, Ψ) + ν]/w
p
=
.
exp −
2
2Kλ ( ρ(A, Σ, Ψ)[δ(X; M, Σ, Ψ) + ν])

Therefore,
W | X ∼ GIG (ρ(A, Σ, Ψ), δ(X; M, Σ, Ψ) + ν, λ) ,
where λ = −(ν + np)/2.
Finally, we note that
X ∼ MVSTn×p (M, A, Σ, Ψ, ν) ⇐⇒ vec(X ) ∼ MSTnp (vec(M), vec(A), Ψ ⊗ Σ, ν), (9)
where MSTnp (·) denotes the multivariate skew-t distribution with location parameter vec(M),
skewness parameter vec(A), scale matrix Ψ⊗Σ, and ν 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 X1 , X2 , . . . XN from an n×p matrix variate skewt distribution. As with the multivariate skew-t distribution, we proceed as if the observed
data is incomplete, and introduce the latent variables W1 , . . . , Wn . The complete-data loglikelihood is then
N
ν  Np
  ν  ν X
Nn
Nν
log
−
log |Σ| −
log |Ψ| − N log Γ
−
log(wi )
`c (ϑ) =C +
2
2
2
2
2
2 i=1
N

N

 1X

1X
+
tr Σ−1 (Xi − M)Ψ−1 A0 +
tr Σ−1 AΨ−1 (Xi − M)0
2 i=1
2 i=1
−

N
N
 1X
1X 1 
wi tr(Σ−1 AΨ−1 A0 ),
tr(Σ−1 (Xi − M)Ψ−1 (Xi − M)0 ) + ν −
2 i=1 wi
2 i=1

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

1) Initialization: Initialize the parameters M, A, Σ, Ψ, ν. Set t = 0
2) E Step: Update ai , bi , ci , where
(t+1)

ai

(t+1)

bi

(t+1)

ci

(t)

= E(Wi |Xi , ϑ̂ )
s
δ(Xi ; M̂(t) , Σ̂(t) , Ψ̂(t) ) + ν̂ (t) Kλ(t) +1 (κ(t) )
=
Kλ(t) (κ(t) )
ρ(Â(t) , Σ̂(t) , Ψ̂(t) )


(t)
1
=E
| Xi , ϑ̂
Wi
s
Kλ(t) +1 (κ(t) )
ρ(Â(t) , Σ̂(t) , Ψ̂(t) )
ν̂ (t) + np
=
+
δ(Xi ; M̂(t) , Σ̂(t) , Ψ̂(t) ) + ν̂ (t) Kλ(t) (κ(t) )
δ(Xi ; M̂(t) , Σ̂(t) , Ψ̂(t) ) + ν (t)

(10)

(11)

(t)

= E(log(Wi )|Xi , ϑ̂ )
s
!
δ(Xi ; M̂(t) , Σ̂(t) , Ψ̂(t) ) + ν̂ (t)
1
∂
Kλ (κ(t) )
+
= log
(t)
(t)
(t)
(t)
Kλ(t) (κ ) ∂λ
ρ(Â , Σ̂ , Ψ̂ )

where
(t)

κ

q
=

(12)
λ=λ(t)

[ρ(Â(t) , Σ̂(t) , Ψ̂(t) )][δ(Xi ; M̂(t) , Σ̂(t) , Ψ̂(t) ) + ν̂ (t) ],

and

ν (t) + np
2
3) First CM Step: Update the parameters M, A, ν.


PN
(t+1) (t+1)
bi
−1
i=1 Xi a
,
M̂(t+1) = PN
(t+1) (t+1)
−N
bi
i=1 a
λ(t) = −

 (t+1)

(t+1)
X
b
−
b
i
i
i=1
,
PN (t+1) (t+1)
a
b
−
N
i
i=1

(13)

PN
Â(t+1) =

(14)

P
P
(t+1)
(t+1)
(t+1)
where a(t+1) = (1/N ) N
and b
= (1/N ) N
.
i=1 ai
i=1 bi
The update for the degrees of freedom cannot be obtained in closed form. Instead we
solve (15) for ν to obtain ν̂ (t+1) .
N
1 X (t+1)
(t+1)
+1−ϕ
−
(b
log
+ ci
) = 0,
2
2
N i=1 i

ν 

ν 

where ϕ(·) is the digamma function.
6

(15)

4) Second CM Step: Update Σ
" N 


0
X (t+1) 
−1
1
Σ̂(t+1) =
bi
Xi − M̂(t+1) Ψ̂(t)
Xi − M̂(t+1)
N p i=1

0 

−1
−1
0
− Â(t+1) Ψ̂(t)
Xi − M̂(t+1) − Xi − M̂(t+1) Ψ̂(t) Â(t+1)
i
−1
0
(t+1)
Â(t+1) Ψ̂(t) Â(t+1)
+ ai

(16)

5) Third CM Step: Update Ψ
" N 
0


X (t+1) 
1
(t+1)
(t+1)
(t+1) −1
(t+1)
Xi − M̂
Xi − M̂
b
Σ̂
Ψ̂
=
N n i=1 i

 
0
0
−1
−1
− Â(t+1) Σ̂(t+1)
Xi − M̂(t+1) − Xi − M̂(t+1) Σ̂(t+1) Â(t+1)
i
−1
0
(t+1)
Â(t+1) Σ̂(t+1) Â(t+1)
+ ai

(17)

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
a(t) =

l(t+1) − l(t)
,
l(t) − l(t−1)

(18)

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.,
(t+1)
l∞
= l(t) +

1
(l(t+1) − l(t) )
(t)
1−a

(cf. Böhning et al., 1994; Lindsay, 1995). As in McNicholas et al. (2010), we stop our EM
(t+1)
algorithms when l∞ − l(t) < , 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 Σ and Ψ are unique only up to a
multiplicative constant. Indeed, if we let Σ̃ = vΣ and Ψ̃ = (1/v)Ψ, v 6= 0, the likelihood
is unchanged. However, we notice that, Ψ̃ ⊗ Σ̃ = Ψ ⊗ Σ, so the estimate of the Kronecker
product is unique.

7

4

Simulations

We conduct two simulations to illustrate the estimation of the
lations, we take 50 different datasets of size 100, from a 3 × 4
Also, in both simulations,



1 −0.5
1
0.5 0.1
 −0.5
1
0.5  ,
Σ =  0.5 1
Ψ=
 0.5 −0.5
0.1 0.5 1
0.1
0.6

parameters. In both simumatrix skew-t distribution.

0.5
0.1
−0.5
0.6 
,
1 −0.4 
−0.4
1

and ν = 4. In simulation 1, we took the location and skewness
respectively, and M2 and A2 in simulation 2, where



0 1 −1
0
1 −1



1 0
0 −1 ,
1 −1
A1 =
M1 =
0 1 −1
0
1 −1



1 −6 −1 −1
5 −4
1 ,
M2 =  −3
1 −4 −1
5

matrix to be M1 and A1 ,

0 1
0 1 ,
0 1




1
−1 0.5
0
A2 =  0.5 −0.5 0.5 0.5  .
0
0 0.5
0

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 Σ and Ψ
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.
8

a

b

0

20

V2

V1

−10

10

−20

0

1.0

1.5

2.0

2.5

3.0

1.0

1.5

2.0

X

2.5

3.0

2.5

3.0

X

c

d

5.0

2.5

V4

V3

20

0.0

10

−2.5

0

−5.0
1.0

1.5

2.0

2.5

3.0

1.0

1.5

2.0

X

X

Figure 1: Typical Marginals for Simulation 1 for (a) V1, (b) V2, (c) V3 and (d) V4. The
red dashed lines denote the means.
Table 1: Component wise averages and standard deviations for the estimated parameters for
simulations 1 and 2.
Simulation

M(sd)

A(sd)

Σ(sd)

Ψ(sd)





1

2

−0.04
 1.01
0.01

0.212
 0.181
0.185

1.04
0.03
1.04
0.176
0.216
0.206

−1.01
0.03
−0.97
0.175
0.158
0.137


−0.02
−1.01 
−0.01

0.176
0.151 
0.146

1.07
 1.01
1.02

0.197
 0.177
0.182



−6.01
4.98
−3.99
0.183
0.180
0.195

−0.99
−3.97
−0.98
0.166
0.202
0.190


−1.02
0.96 
4.99

0.153
0.159 
0.147



0.99
 −2.98
1.00

0.170
 0.218
0.177



1.03
 0.50
0.01

0.165
 0.147
0.121

−1.06
−1.04
−1.03
0.174
0.192
0.201

0.03
−0.01
−0.01
0.120
0.113
0.088



1.04
0.91

 0.45
1.03
1.04
0.09
 
0.180
0.111
0.167   0.074
0.169
0.055

−1.02
−0.51
−0.02
0.183
0.134
0.127

0.51
0.49
0.50
0.125
0.140
0.132



0.01
0.92 0.45
 0.45 0.91
0.52 
0.00
0.08 0.45
 
0.099
0.104 0.066
0.133   0.066 0.110
0.112
0.037 0.067

0.45
0.92
0.47
0.074
0.108
0.069

1.08 −0.53
 −0.53
1.06

 0.53 −0.53
0.12
0.64

0.083 0.074
 0.074 0.067

 0.061 0.058
0.072 0.064

1.09 −0.54
 −0.54

1.07

0.08
 0.54 −0.53
0.45 
0.11
0.63
0.92


0.087 0.068
0.037
 0.068 0.083

0.067  
 0.064 0.055
0.103
0.068 0.079


0.09
0.47 
0.92

0.055
0.069 
0.096

0.53
−0.53
1.08
−0.44
0.061
0.058
0.074
0.058
0.54
−0.53
1.06
−0.42
0.064
0.055
0.066
0.055

ν(sd)



0.12
0.64 

−0.44 
1.09

0.072
 4.22
0.064 
 (0.63)
0.058 
0.089

0.11
0.63 

−0.42 
1.05

0.068
 4.22
0.079 
 (0.92)
0.055 
0.090

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

9

a

b

15

10
0

V2

V1

5

−10
0

−5
−20

1.0

1.5

2.0

2.5

3.0

1.0

1.5

2.0

X

X

c

d

2.5

3.0

2.5

3.0

10

4

V4

V3

5

0

0

−4

1.0

1.5

2.0

2.5

3.0

X

1.0

1.5

2.0

X

Figure 2: Typical Marginals for Simulation 2 for (a) V1, (b) V2, (c) V3 and (d) V4. The
red dashed lines denote the means.

References
Aitken, A. C. (1926). A series formula for the roots of algebraic and transcendental equations.
Proceedings of the Royal Society of Edinburgh 45, 14–22.
Anderlucci, L., C. Viroli, et al. (2015). Covariance pattern mixture models for the analysis
of multivariate heterogeneous longitudinal data. The Annals of Applied Statistics 9 (2),
777–800.
Böhning, D., E. Dietz, R. Schaub, P. Schlattmann, and B. Lindsay (1994). The distribution
of the likelihood ratio for mixtures of densities from the one-parameter exponential family.
Annals of the Institute of Statistical Mathematics 46, 373–388.
Chen, J. T. and A. K. Gupta (2005). Matrix variate skew normal distributions. Statistics 39 (3), 247–253.
Doğru, F. Z., Y. M. Bulut, and O. Arslan (2016). Finite mixtures of matrix variate t
distributions. Gazi University Journal of Science 29 (2), 335–341.

10

Domı́nguez-Molina, J. A., G. González-Farı́as, R. Ramos-Quiroga, and A. K. Gupta (2007).
A matrix variate closed skew-normal distribution with applications to stochastic frontier
analysis. Communications in Statistics–Theory and Methods 36 (9), 1691–1703.
Dutilleul, P. (1999). The MLE algorithm for the matrix normal distribution. Journal of
Statistical Computation and Simulation 64 (2), 105–123.
Gupta, A. K. and D. K. Nagar (1999). Matrix Variate Distributions. Boca Raton: Chapman
& Hall/CRC Press.
Harrar, S. W. and A. K. Gupta (2008). On matrix variate skew-normal distributions. Statistics 42 (2), 179–194.
Lindsay, B. G. (1995). Mixture models: Theory, geometry and applications. In NSF-CBMS
Regional Conference Series in Probability and Statistics, Volume 5. California: Institute
of Mathematical Statistics: Hayward.
McNicholas, P. D. (2016). Mixture Model-Based Classification. Boca Raton: Chapman &
Hall/CRC Press.
McNicholas, P. D., T. B. Murphy, A. F. McDaid, and D. Frost (2010). Serial and parallel
implementations of model-based clustering via parsimonious Gaussian mixture models.
Computational Statistics and Data Analysis 54 (3), 711–723.
Meng, X.-L. and D. B. Rubin (1993). Maximum likelihood estimation via the ECM algorithm: a general framework. Biometrika 80, 267–278.
Murray, P. M., R. B. Browne, and P. D. McNicholas (2014a). Mixtures of skew-t factor
analyzers. Computational Statistics and Data Analysis 77, 326–335.
Murray, P. M., P. D. McNicholas, and R. B. Browne (2014b). A mixture of common skew-t
factor analyzers. Stat 3 (1), 68–82.
Wishart, J. (1928). The generalised product moment distribution in samples from a normal
multivariate population. Biometrika, 32–52.

11

