Drawing from the multivariate Student's t-distribution So I am trying to figure out if there is a nice decomposition for sampling from the multivariate Student's t-distribution like there is for sampling from the multivariate normal distribution:  http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
There are quite a few R packages out there for doing this but none of them seem intuitive to me to implement.  What I am trying to accomplish it to sample from the multivariate t-distribution given that I have a specified number of degrees of freedom, a known mean and covariance matrix.
I need to be able to sample from the multivariate t-distribution because I have analytically integrated out parameters from my posterior distribution.  I.e., 
$$p(\psi|y)=\int_{R^p}\int_{R^+}p(\beta,\sigma,\psi|y)d\sigma^2d\beta$$
which after the integration results in $\psi\sim\text{Multivariate Student-t}(\mu,\Sigma,\nu)$. Since I have a closed analytical form for $\psi$, I would like to be able to sample directly from the multivariate Student t-distribution.
 A: The Wikipedia article on multivariate t explains that when $Y$ has a $N(0, \Sigma)$ distribution and independently $U$ has a $\chi_\nu^2$ distribution, then
$$X = \mu + \frac{Y}{\sqrt{\frac{U}{\nu}}} = \mu + Y\sqrt{\frac{\nu}{U}}$$
has a $t_\nu(\mu, \Sigma)$ distribution.
This generalizes the univariate case of Student's t, which is the ratio of a standard Normal and a (scaled) $\chi$ variate: we recognize in this formula the classic ratio of a deviation (expected to be zero) to a standard error, which is predicted by the null hypothesis of the Student t-test.
Whence realizations of $X$ ($p$-vectors) can be had from the same algebraic combination of independent realizations of $Y$ (also $p$-vectors) and $U$ (numbers), at a cost therefore of generating $p+1$ random variates for every $p$ coefficients of $X$: that's reasonably efficient.
This R code illustrates the case $p=2$.
library(MASS) # mvrnorm()
mu <- c(7, 5); sigma <- matrix(c(1, 1/2, 1/2, 1), 2); nu <- 4
n <- 4e3 # Number of draws
y <- t(t(mvrnorm(n, rep(0, length(mu)), sigma) * sqrt(nu / rchisq(n, nu))) + mu)

You can check that colMeans(y) approximates mu, var(y) approximates sigma multiplied by $\nu/(\nu-2)$ (an observation kindly offered by Stéphane Laurent, whose claim is justified because the variance of a Student t distribution is $\nu/(\nu-2)$), and you can plot(y) to see the distribution.

