Simulating Matrix-T Variables By Scale Mixtures

Why Simulate a Matrix-T Variable as a Mixture?

A matrix-normal model has Gaussian tails conditional on its covariance parameters. A matrix-t model can instead represent uncertainty in a precision matrix: conditional on that precision, the draw is matrix-normal, but the precision itself varies from draw to draw. The resulting mixture has heavier tails. This representation is useful both for generating matrix-t samples and for checking how degrees of freedom affect extreme observations.

We use $X$ for a $p\times q$ random matrix, with location $M$, row scale matrix $\Sigma$, and column scale matrix $\Omega$. Both scale matrices are symmetric positive definite. Our convention is $X\mid W\sim MN_{p,q}(M,W^{-1},\Omega)$: the conditional covariance of $\operatorname{vec}(X)$ is $\Omega\otimes W^{-1}$. The mixing matrix $W$ is a row precision, not a row covariance.

Definition: Let $\nu>p-1$, and let $W\sim \operatorname{Wishart}_p(\nu,\Sigma^{-1})$, using the convention that a Wishart distribution with scale $V$ has mean $\nu V$. Conditional on $W$, draw $X\sim MN_{p,q}(M,W^{-1},\Omega)$. The marginal distribution of $X$ is the matrix-t distribution with degrees of freedom $\nu$, location $M$, row scale $\Sigma$, and column scale $\Omega$, under the scale convention made explicit below.

The Two Random Draws

A simulation has two stages. First draw a positive-definite Wishart precision matrix. Then, given that matrix, draw a matrix-normal observation whose row covariance is its inverse. If $L L'=W^{-1}$ and $C C'=\Omega$ are lower-triangular Cholesky factorizations, and $Z$ is a $p\times q$ matrix of independent standard normal variables, the second stage is

$$ X=M+LZC'. $$

Indeed, conditional on $W$, the mean is $M$, and $\operatorname{vec}(LZC')=(C\otimes L)\operatorname{vec}(Z)$. Since $\operatorname{vec}(Z)$ has identity covariance, the resulting covariance is $(CC')\otimes(LL')=\Omega\otimes W^{-1}$, as required.

For real degrees of freedom $\nu>p-1$, the Wishart draw can be generated by Bartlett’s decomposition. Let $A$ be lower triangular, with independent entries: $A_{ii}^2\sim\chi^2_{\nu-i+1}$ for $i=1,\ldots,p$, and $A_{ij}\sim N(0,1)$ for $i>j$. If $C_\Sigma C_\Sigma'=\Sigma^{-1}$, then

W=C_\Sigma A A'C_\Sigma'.

The chi-square degrees of freedom are positive because $\nu>p-1$. This construction uses the Wishart scale $\Sigma^{-1}$, not $\Sigma$. When $\nu$ is an integer, one alternative is to sum $\nu$ independent Gaussian outer products with covariance $\Sigma^{-1}$; Bartlett’s construction also covers noninteger $\nu$.

Worked Example: A Scalar Draw From the Mixture

Take $p=q=1$, location $M=2$, row scale $\Sigma=4$, column scale $\Omega=9$, and $\nu=5$. In one possible draw, suppose the chi-square variate used for the Wishart precision is $4$. Since the scale of $W$ is $\Sigma^{-1}=1/4$, the resulting precision is

W=\frac{1}{4}\times 4=1.

The conditional variance is $\Omega/W=9$, so the conditional standard deviation is $3$. If the standard-normal draw is $Z=1.2$, then

X=M+\frac{\sqrt{\Omega}}{\sqrt W}Z =2+\frac{3}{1}(1.2)=5.6.

This is one possible matrix-t draw. The random scale comes from the Wishart stage; the normal draw supplies the conditional Gaussian fluctuation.

Why the Mixture Has a Matrix-T Density

The Wishart density, for $W$ positive definite, is proportional to \[ |W|^{(\nu-p-1)/2}\exp\{-\operatorname{tr}(\Sigma W)/2\}. \] The conditional matrix-normal density contributes a factor $|W|^{q/2}$. Define $D=X-M$ and $B=D\Omega^{-1}D'$. Its exponential term is $\exp\{-\operatorname{tr}(WB)/2\}$. Thus the product of the conditional and mixing densities has a $W$-dependent part

|W|^{(\nu+q-p-1)/2} \exp\{-\operatorname{tr}[W(\Sigma+B)]/2\}.

For $a>(p-1)/2$, the matrix gamma integral is \[ \int_{W>0}|W|^{a-(p+1)/2} \exp\{-\operatorname{tr}(HW)/2\}\,dW =2^{ap}\Gamma_p(a)|H|^{-a}, \] where $H$ is positive definite and $\Gamma_p$ is the multivariate gamma function. Here $a=(\nu+q)/2$, which satisfies the condition because $\nu>p-1$ and $q\geq1$. Integrating out $W$ therefore gives a density proportional to $|\Sigma+B|^{-(\nu+q)/2}$.

Derived result: Under the stated mixture, the density of $X$ is $$ f(X)= \frac{\Gamma_p((\nu+q)/2)} {\Gamma_p(\nu/2)\pi^{pq/2}|\Sigma|^{q/2}|\Omega|^{p/2}} \left|I_p+\Sigma^{-1}(X-M)\Omega^{-1}(X-M)'\right|^{-(\nu+q)/2}. $$ The determinant in the density is well-defined and positive: its matrix is similar to $I_p+\Sigma^{-1/2}B\Sigma^{-1/2}$, which is positive definite. This formula fixes the scale convention: $\Sigma$ and $\Omega$ are scale matrices in the density, not generally the marginal covariance matrices.

The determinant power describes how the mixture changes the Gaussian decay. The degrees of freedom also control which moments exist. For example, in the scalar case, $W\sim\operatorname{Gamma}(\nu/2,\text{ scale }2/\Sigma)$. Its inverse has finite expectation only when $\nu>2$, and then $E[W^{-1}]=\Sigma/(\nu-2)$. Thus for $p=q=1$, the marginal variance is $\Sigma\Omega/(\nu-2)$ when $\nu>2$, rather than $\Sigma\Omega$.

Worked Example: A Two-by-Two Matrix Draw

Let $p=q=2$, $\nu=5$, $M=0$, and $\Sigma=\Omega=I_2$. For illustration, suppose the Bartlett variables take the values $A_{11}=2$, $A_{21}=1$, and $A_{22}=\sqrt{3}$. These correspond to $A_{11}^2=4$, $A_{22}^2=3$, and one standard-normal off-diagonal draw of $1$. Then

A=\begin{pmatrix}2&0\\1&\sqrt3\end{pmatrix}, \qquad W=AA'=\begin{pmatrix}4&2\\2&4\end{pmatrix}.

The determinant of $W$ is $12$, so \[ W^{-1}=\frac{1}{12}\begin{pmatrix}4&-2\\-2&4\end{pmatrix} =\begin{pmatrix}1/3&-1/6\\-1/6&1/3\end{pmatrix}. \] A lower-triangular Cholesky factor is \[ L=\begin{pmatrix}1/\sqrt3&0\\-1/(2\sqrt3)&1/2\end{pmatrix}, \] as direct multiplication confirms: $LL'=W^{-1}$. With the standard-normal matrix draw \[ Z=\begin{pmatrix}1&0\\0&-1\end{pmatrix}, \] and $C=I_2$, the matrix-t draw is

X=LZ= \begin{pmatrix}1/\sqrt3&0\\-1/(2\sqrt3)&-1/2\end{pmatrix}.

The nonzero off-diagonal entries in $W^{-1}$ induce dependence between the two rows conditional on $W$. The displayed $X$ is one draw, not a summary of the marginal distribution.

Checking the Effect of Degrees of Freedom

A useful simulation check is to compare scalar contrasts with the same univariate t scale and different degrees of freedom. In the scalar mixture, set $W\sim\operatorname{Wishart}_1(\nu,1/\nu)$, so $W=U/\nu$ with $U\sim\chi^2_\nu$. Draw $X\mid W\sim N(0,1/W)$. Then $X=Z/\sqrt{U/\nu}$, the standard Student t variable with $\nu$ degrees of freedom. Holding this scale fixed matters: using $\Sigma=1$ instead would give $Z/\sqrt U$, whose scale changes with $\nu$.

Worked Example: Comparing Tail Probabilities

Compare standard t variables with $\nu=3$ and $\nu=15$, using the threshold $|X|>2$. For a standard t variable, \[ P(|X|>x)=I_{\nu/(\nu+x^2)}(\nu/2,1/2), \] where $I_y(a,b)$ is the regularized incomplete beta function. At $x=2$, this gives $I_{3/7}(3/2,1/2)\approx0.139$ for $\nu=3$, and $I_{15/19}(15/2,1/2)\approx0.064$ for $\nu=15$. The smaller degrees of freedom produce the larger probability of exceeding this threshold.

In $N=10{,}000$ independent simulated draws, the corresponding expected exceedance counts are $10{,}000(0.139)\approx1{,}390$ and $10{,}000(0.064)\approx640$. These are expectations, not guaranteed simulation outputs. For a simulated count $K$, compare $K/N$ with the theoretical probability; the Monte Carlo standard error is approximately $\sqrt{p(1-p)/N}$. At these probabilities and sample size, it is about $0.0035$ for $p=0.139$ and $0.0024$ for $p=0.064$.

Simulation Procedure and Common Pitfall

For each independent matrix-t replicate, use a fresh Wishart draw and a fresh standard-normal matrix. The following sequence makes the roles of the parameters explicit.

Set valid parameters.
Choose $\nu>p-1$, a $p\times p$ positive-definite $\Sigma$, a $q\times q$ positive-definite $\Omega$, and a $p\times q$ location $M$.

Draw the precision.
Generate $W\sim\operatorname{Wishart}_p(\nu,\Sigma^{-1})$, for example by Bartlett’s decomposition.

Draw conditionally.
Factor $W^{-1}=LL'$ and $\Omega=CC'$, generate independent standard-normal entries in $Z$, and set $X=M+LZC'$.

Check the output.
Repeat independently, then compare summaries or contrast exceedance frequencies with theoretical quantities whose moments or probabilities exist.

Common mistake: drawing a covariance when the construction calls for a precision. If $W\sim\operatorname{Wishart}_p(\nu,\Sigma^{-1})$, using $W$ as the conditional row covariance instead of $W^{-1}$ changes the model. The conditional covariance must be $\Omega\otimes W^{-1}$, so the matrix-normal draw uses a Cholesky factor of $W^{-1}$. Another frequent error is to compare different degrees of freedom while silently changing the univariate t scale; state the mixing scale before interpreting a tail comparison.

This mixture is useful because it turns a heavy-tailed matrix model into familiar conditional Gaussian and Wishart draws. It also makes simulation checks interpretable: the Gaussian randomness controls variation within a precision draw, while the random precision controls how often unusually large conditional scales occur.

Check Your Understanding

Use the mixture definition and scale convention to answer the following.

For $p=3$, what condition on $\nu$ ensures that the Wishart precision draw is well-defined?
In the conditional law, which matrix is the row covariance: $W$ or $W^{-1}$?
Explain why $M+LZC'$ has conditional covariance $\Omega\otimes W^{-1}$ when $LL'=W^{-1}$ and $CC'=\Omega$.
For the scalar standard-t mixture, why must the Wishart scale be $1/\nu$ rather than $1$?
What should be compared between simulated exceedance frequencies and theoretical tail probabilities, and why will they not usually be identical?