ML estimation of the Matrix Variate t Distribution
Estimation of the Matrix Variate t Distribution
The parameters of the multivariate t distribution can be estimated using the EM algorithm. An EM algorithm for the multivariate t-distribution with no missing data was provided by Rubin (1983). This was variously extended and refined, as the EM algorithm can be quite slow to converge. One set of refinements split the M-step into a series of conditional maximization (CM) steps (Meng and Rubin 1993), and then a later refinement (called ECME - Expectation/Conditional Maximization Either) allowed for either conditional maximization steps or maximizing the constrained actual likelihood (Liu and Rubin 1994), which leads to dramatically faster convergence. The matrix variate t is an extension of the multivariate version (Dickey 1967); however, unlike the normal distribution, it cannot be treated as a rearrangement of the multivariate case, so the addition of the extra dimension poses a non-trivial problem for extending the results.
An n × p random matrix X is distributed as a matrix variate t random variable if it has probability density function as follows: $$f({\mathbf {X}};\nu,{\mathbf {M}}, {\boldsymbol \Omega }, {\boldsymbol \Sigma }) = {\frac {\Gamma _{p}\left({\frac {\nu +n+p-1}{2}}\right)}{(\pi )^{{\frac {np}{2}}}\Gamma _{p}\left({\frac {\nu +p-1}{2}}\right)}}|{\boldsymbol \Omega }|^{{-{\frac {n}{2}}}}|{\boldsymbol \Sigma }|^{{-{\frac {p}{2}}}}\times \left|{\mathbf {I}}_{n}+{\boldsymbol \Sigma }^{{-1}}({\mathbf {X}}-{\mathbf {M}}){\boldsymbol \Omega }^{{-1}}({\mathbf {X}}-{\mathbf {M}})^{{{\rm {T}}}}\right|^{{-{\frac {\nu +n+p-1}{2}}}}$$
With Ω and Σ covariance matrices of appropriate
dimension, M a mean
matrix, ν the degrees of
freedom parameter, and Γp() the
multivariate gamma function, which is implemented in the CholWishart
package available on CRAN.
For p = 1 and Σ = ν (or q = 1 and Ω = ν), this reduces to the familiar multivariate t distribution. If the p × q dimensional random variable X follows the matrix variate t distribution t(ν, M, Σ, Ω) with center M, positive definite spread matrices Σ (p × p) and Ω (q × q), and degrees of freedom ν, it can be shown that, given a certain weight matrix S, X has a matrix variate normal distribution and that S is Wishart-distributed. Specifically, Additionally, it can be shown the Wishart is the conjugate prior distribution for the parameter Σ, and thus the conditional posterior distribution of S, i.e., its distribution given (M, Ω, Σ, ν, X) is S|X, M, Ω, Σ, ν ∼ Wp(ν + p + q − 1, [(X − M)Ω−1(X − M)T + Σ]−1)
For an observed set of Xi, i = 1, 2, …, N, then, we can construct an ECM-style algorithm by augmenting our data with a set of weights, Si, and use this to estimate the parameters of the matrix variate t distribution. Briefly, we define a set of complete data sufficient statistics based on the Si and update them with their expectations in the E-step. Then, the other parameters are maximized in the CM steps.
The complete data sufficient statistics for the quantities to be estimated are: $$S_{SX} =\sum_{i = 1}^N S_i X_i; \quad S_S = \sum_{i = 1}^N S_i; \quad S_{XSX} = \sum_{i = 1}^N X_i^T S_i X_i \quad S_{|S|} = \sum_{i = 1}^N \log |Si|$$
E-step
Define at step t the set Θ(t) = (Xobs, ν(t), M(t), Σ(t), Ω(t)). Then, given these values, we have, based on the properties of the Wishart distribution, Si(t + 1) = E(Si|Θ(t)) = (ν(t) + p + q − 1)[(Xi − M(t))Ω(t)−1(Xi − M(t))T + Σ(t)]−1 From this we can derive the other expected values: $$ \begin{align*} S_{S}^{(t+1)} &= E(S_S | \Theta^{(t)}) = \sum_{i=1}^N S_{i}^{(t+1)} \\ S_{SX}^{(t+1)} &= E(S_{SX} | \Theta^{(t)}) = \sum_{i=1}^NE(S_iX_i|\Theta^{(t)}) = \sum_{i=1}^N S_{i}^{(t+1)}X_i \\ S_{XSX}^{(t+1)} &= E(S_{XSX} | \Theta^{(t)}) = \sum_{i=1}^NE(X_i^TS_iX_i|\Theta^{(t)}) = \sum_{i=1}^N X_i^TS_{i}^{(t+1)}X_i \end{align*} $$ This last value only needs to be computed if ν is unknown: $$ \begin{align*} S_{|S|}^{(t+1)} &= E(S_{|S|} | \Theta^{(t)}) = E \left[\sum_{i = 1}^N \log| S_i| \big| \Theta^{(t)} \right]\\ &=N\psi_{p}\left({\frac {\nu^{(t)} + p + q -1}{2}}\right) + Np \log 2 +\sum_{i=1}^N \log \left|\frac{S_{i,obs}^{(t+1)}}{\nu^{(t)} + p + q -1}\right| \end{align*} $$
CM Steps
Similarly to the case of the multivariate t, it is more efficient to partition the maximization step into multiple steps. This is an ECME algorithm that first maximizes the expected log-likelihood for (M, Σ, Ω) and then maximizes the actual log-likelihood over ν given (M, Σ, Ω), similar to .
Based on the updated weight matrices Si(t + 1) and statistics based on Θ(t) and Xobs,
$$ \begin{align*} \widehat{M} &= \left(\sum_{i = 1}^NS_i^{(t+1)}\right)^{-1}\sum_{i=1}^NS_iX_i = {S_{S}^{(t+1)}}^{-1} S_{SX}^{(t+1)} \\ \widehat{\Omega} &= \frac{1}{Np} \sum_{i = 1}^N (X_i - {M}^{(t)})^T S_i^{(t+1)} (X_i - {M}^{(t)}) =\frac{1}{Np} \left(S_{XSX}^{(t+1)} - {S_{SX}^{(t+1)}}^T {S_{S}^{(t+1)}}^{-1}{S_{SX}^{(t+1)}} \right) \\ \widehat{\Sigma}^{-1} &= \frac{1}{N(\nu^{(t)}+p-1)}\sum_{i = 1}^N S_i^{(t+1)} = \frac{S_S^{(t+1)}}{N(\nu^{(t)}+p-1)} \end{align*} $$ And, again, with the set of Si observed, the MLE of ν can be obtained: $$ N \frac{d}{d\nu}\log \Gamma_p ((\nu+p-1)/2) - \frac{1}{2} (S_{|S|} - Np\log2 + N \log |\widehat{\Sigma}|) = 0$$ Note that ν > p − 1. Specifically, we have: $$ \begin{align*} 0 &= N \psi_p((\nu+p-1)/2) - \left(N\psi_{p}\left({\frac {\nu + p + q -1}{2}}\right) +\sum_{i=1}^N \log \left|\frac{S_{i,obs}^{(t+1)}}{\nu + p + q -1}\right| - N \log \left|\frac{S_S^{(t+1)}}{N(\nu+p-1)}\right|\right) \\* &= \psi_p((\nu+p-1)/2) - \left(\psi_{p}\left({\frac {\nu + p + q -1}{2}}\right) +\frac{1}{N}\sum_{i=1}^N \log \left|Z_{i,obs}^{(t+1)}\right| + p \log \frac{N(\nu + p -1)}{\nu + p +q -1} - \log \left|Z_S^{(t+1)}\right|\right) \end{align*}$$ where Z* is the appropriate S* statistic with (ν + p + q − 1) factored out and ψp is the p-dimensional digamma function. This can be solved for ν using a 1-dimensional search.
Applications and results
With ν known
If the degrees of freedom parameter, ν known, the estimation is fairly
straightforward. The procedure is similar to the multivariate t or the matrix variate normal. In
this case, the interface is just like the interface for the
MLmatrixnorm() function:
set.seed(20190622)
sigma = (1/7) * rWishart(1, 7, 1*diag(3:1))[,,1]
A = rmatrixt(n=100,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
V = sigma, df = 7)
results=MLmatrixt(A, df = 7)
print(results)
#> $mean
#> [,1] [,2] [,3]
#> [1,] 99.8291911 -1.000708e+02 25.01042
#> [2,] -0.0425028 -8.261356e-03 -1000.02357
#>
#> $U
#> [,1] [,2]
#> [1,] 1.00000000 -0.08039981
#> [2,] -0.08039981 0.97069319
#>
#> $V
#> [,1] [,2] [,3]
#> [1,] 1.0000000 0.21375613 0.11094767
#> [2,] 0.2137561 0.14872070 0.04098332
#> [3,] 0.1109477 0.04098332 0.12860660
#>
#> $var
#> [1] 5.148569
#>
#> $nu
#> [1] 7
#>
#> $iter
#> [1] 44
#>
#> $tol
#> [1] 1.155394e-08
#>
#> $logLik
#> [1] -386.7623
#>
#> $convergence
#> [1] TRUE
#>
#> $call
#> MLmatrixt(data = A, df = 7)There are two restrictions possible for the mean matrices:
row.mean = TRUE will force a common mean within a row and
col.mean = TRUE will force a common mean within a column.
Setting both will ensure a constant mean for the entire system.
Restrictions on U and
V, the row-wise
variance and column-wise variance, are possible with
row.variance and col.variance commands.
The options for variance restrictions are the same as for the
MLmatrixnorm() function. Currently the options for variance
restrictions are to impose an AR(1) structure by providing the
AR(1) option, a compound symmetry structure by providing
the CS option, to impose a correlation matrix structure by
specifying correlation or corr, or to impose
an identical and independent structure by specifying
Independent or I. This works by using
uniroot to find the appropriate ρ which sets the derivative of the
log-likelihood to zero for the AR and CS
options - it is not fast but if this is the true structure it will be
better in some sense than an unstructured variance matrix. The ρ parameter should be > 0 and is forced to be non-negative. If
the data behaves incompatibly with those restrictions, the function will
provide a warning and exit with the current model fit.
With ν unknown
Estimation of ν, the
degrees of freedom parameter, is slow and the principal mathematical
difficulty of the matrix-variate t distribution. It is performed
using ECME. Generally, a fair amount more data are needed in order to
have good convergence properties for the estimator, but they have not
been derived analytically. Here you can see the recovery of the
parameter for a few sample sizes. Because of the relative slowness of
running a longer simulation, this only includes one set of examples. I
give code for a longer and larger simulation than what is plotted if
you’re really interested below. What is plotted below is only the
df = 10 example with 75 trials and a maximum number of
iterations max.iter = 20. The full simulation may take
several minutes.
| samplesize | min | max | median | mean | sd |
|---|---|---|---|---|---|
| 25 | 7.059192 | 15.76929 | 10.316514 | 10.512701 | 2.0603677 |
| 50 | 7.220806 | 14.47148 | 10.424858 | 10.415426 | 1.6347640 |
| 100 | 8.281614 | 11.85454 | 9.826301 | 9.859941 | 0.8019631 |
As expected, increased sample size leads to better results in recovering the parameter. The results for smaller sample sizes would be more divergent if left to run until convergence.
Use for classification
Using the t distribution
works in both matrixlda() and matrixqda() as
expected by specifying method = "t" and providing either a
single parameter (for lda or qda) for the
degrees of freedom or a vector as long as the number of classes (for
qda). Additional parameters for fitting can be passed
through the ... to MLmatrixt() just as for the
normal case, including estimating the degrees of freedom parameter. The
qda will only estimate nu with it varying
between groups, it will not estimate a common nu.
A <- rmatrixt(30, mean = matrix(0, nrow=2, ncol=3), df = 10)
B <- rmatrixt(30, mean = matrix(c(1,0), nrow=2, ncol=3), df = 10)
C <- rmatrixt(30, mean = matrix(c(0,1), nrow=2, ncol=3), df = 10)
ABC <- array(c(A,B,C), dim = c(2,3,90))
groups <- factor(c(rep("A",30),rep("B",30),rep("C",30)))
prior = c(30,30,30)/90
matlda <- matrixlda(x = ABC,grouping = groups, prior = prior,
method = 't', nu = 10, fixed = TRUE)
predict(matlda, newdata = ABC[,,c(1,31,61)])
#> $class
#> [1] A B C
#> Levels: A B C
#>
#> $posterior
#> [,1] [,2] [,3]
#> [1,] 9.997064e-01 2.849389e-04 8.693417e-06
#> [2,] 4.174839e-05 9.999508e-01 7.439619e-06
#> [3,] 8.908508e-06 4.857627e-06 9.999862e-01Session info
This vignette was built using rmarkdown.
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] magrittr_2.0.3 dplyr_1.1.4 ggplot2_3.5.1 MixMatrix_0.2.8
#> [5] rmarkdown_2.28
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_1.8.9 highr_0.11 compiler_4.4.1
#> [5] tidyselect_1.2.1 Rcpp_1.0.13 jquerylib_0.1.4 scales_1.3.0
#> [9] yaml_2.3.10 fastmap_1.2.0 R6_2.5.1 labeling_0.4.3
#> [13] generics_0.1.3 knitr_1.48 tibble_3.2.1 maketools_1.3.1
#> [17] munsell_0.5.1 bslib_0.8.0 pillar_1.9.0 rlang_1.1.4
#> [21] utf8_1.2.4 cachem_1.1.0 xfun_0.48 sass_0.4.9
#> [25] sys_3.4.3 cli_3.6.3 withr_3.0.2 digest_0.6.37
#> [29] grid_4.4.1 CholWishart_1.1.4 lifecycle_1.0.4 vctrs_0.6.5
#> [33] evaluate_1.0.1 glue_1.8.0 farver_2.1.2 buildtools_1.0.0
#> [37] fansi_1.0.6 colorspace_2.1-1 tools_4.4.1 pkgconfig_2.0.3
#> [41] htmltools_0.5.8.1All the code for easy copying
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
set.seed(20190622)
sigma = (1/7) * rWishart(1, 7, 1*diag(3:1))[,,1]
A = rmatrixt(n=100,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
V = sigma, df = 7)
results=MLmatrixt(A, df = 7)
print(results)
### Here is the long simulation
library(ggplot2)
set.seed(20181102)
df = c(5, 10, 20)
df5 <- rep(0,200)
df10 <- rep(0,200)
df100 <- rep(0,200)
df550 <- rep(0,200)
df1050 <- rep(0,200)
df2050 <- rep(0,200)
df5100 <- rep(0,200)
df10100 <- rep(0,200)
df20100 <- rep(0,200)
meanmat = matrix(0,5,3)
U = diag(5)
V = diag(3)
for(i in 1:200){
df5[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 5, n = 35, U =U, V =V), fixed = FALSE)$nu
df10[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 10, n = 35, U =U, V =V), fixed = FALSE)$nu
df100[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 20, n = 35, U =U, V =V), fixed = FALSE)$nu
df550[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 5, n = 50, U =U, V =V), fixed = FALSE)$nu
df1050[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 10, n = 50, U =U, V =V), fixed = FALSE)$nu
df2050[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 20, n = 50, U =U, V =V), fixed = FALSE)$nu
df5100[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 5, n = 100, U =U, V =V), fixed = FALSE)$nu
df10100[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 10, n = 100, U =U, V =V), fixed = FALSE)$nu
df20100[i] = MLmatrixt(rmatrixt(mean = meanmat,
df = 20, n = 100, U =U, V =V), fixed = FALSE)$nu
}
truedataframe = data.frame(truedf = factor(c(5,10,20),
label = c('5 df', '10 df', '20 df')),
estdf = c(5,10,20))
dfdataframe = data.frame(truedf = factor(rep(rep(c(5,10,20), each = 200),3),
label = c('5 df', '10 df', '20 df')),
estdf = c(df5, df10, df100, df550, df1050, df2050, df5100, df10100, df20100),
samplesize = factor(rep(c(35,50,100), each = 600)))
library(tidyverse)
denseplot <- ggplot(data = subset(dfdataframe, estdf < 200),
aes(x=estdf, fill=samplesize)) +
geom_density(alpha = .5) +
geom_vline(data = truedataframe,
mapping = aes(xintercept = estdf),
size = .5) +
theme_bw() +
theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(),
strip.text = element_text(size = 8),
legend.justification=c(1,0), legend.position=c(.95,.4),
legend.background = element_blank(),
legend.text =element_text(size = 8), legend.title = element_text(size = 8)) +
ggtitle("Density plot of estimated degrees of freedom compared to actual") +
xlab(NULL) +
ylab(NULL) +
scale_fill_manual(values = c("#050505", "#E69F00", "#56B4E9"),
name = "Sample Size") +
facet_wrap(factor(truedf)~., scales="free") +
NULL
denseplot
knitr::kable(dfdataframe %>% group_by(truedf, samplesize) %>%
summarize(min = min(estdf), max = max(estdf),
median = median(estdf),
mean=mean(estdf),
sd = sd(estdf)))
### Here ends the long simulation
##### Here is what is really run
set.seed(20190621)
df10 <- rep(0,50)
df1050 <- rep(0,50)
df10100 <- rep(0,50)
for(i in 1:50){
df10[i] = suppressWarnings(MLmatrixt(rmatrixt(mean = matrix(0,5,3),df = 10, n = 25), fixed = FALSE, df = 5, max.iter = 20)$nu)
df1050[i] = suppressWarnings(MLmatrixt(rmatrixt(mean = matrix(0,5,3),df = 10, n = 50), fixed = FALSE, df = 5, max.iter = 20)$nu)
df10100[i] = suppressWarnings(MLmatrixt(rmatrixt(mean = matrix(0,5,3),df = 10, n = 100), fixed = FALSE, df = 5, max.iter = 20)$nu)
}
dfdataframe = data.frame(label = c('10 df'),
estdf = c(df10, df1050, df10100),
samplesize = factor(rep(c(25,50,100), each = 50)))
library(ggplot2)
library(dplyr)
library(magrittr)
denseplot <- ggplot(data = subset(dfdataframe, estdf < 200),aes(x=estdf, fill=samplesize)) +
geom_density(alpha = .5) +
geom_vline(mapping = aes(xintercept = 10),
size = .5) +
theme_bw() +
theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), strip.text = element_text(size = 8),
legend.justification=c(1,0), legend.position=c(.95,.4),
legend.background = element_blank(),
legend.text =element_text(size = 8), legend.title = element_text(size = 8)) +
ggtitle("Density plot of estimated degrees of freedom compared to actual") +
xlab(NULL) +
ylab(NULL) +
scale_fill_manual(values = c("#050505", "#E69F00", "#56B4E9"), name = "Sample Size") +
# facet_wrap(factor(truedf)~., scales="free") +
NULL
denseplot
knitr::kable(dfdataframe %>% group_by(samplesize) %>%
summarize(min = min(estdf), max = max(estdf),
median = median(estdf),
mean=mean(estdf),
sd = sd(estdf)))
#### Here ends what is really run
A <- rmatrixt(30, mean = matrix(0, nrow=2, ncol=3), df = 10)
B <- rmatrixt(30, mean = matrix(c(1,0), nrow=2, ncol=3), df = 10)
C <- rmatrixt(30, mean = matrix(c(0,1), nrow=2, ncol=3), df = 10)
ABC <- array(c(A,B,C), dim = c(2,3,90))
groups <- factor(c(rep("A",30),rep("B",30),rep("C",30)))
prior = c(30,30,30)/90
matlda <- matrixlda(x = ABC,grouping = groups, prior = prior,
method = 't', nu = 10, fixed = TRUE)
predict(matlda, newdata = ABC[,,c(1,31,61)])
sessionInfo()