Skewed Multivariate Student-t distribution

GNU R package miscF has a function rmvst to draw samples from a skewed multivariate Student-t distribution. I translated it to Stan:


data {
  real<lower=1> nu;
  int p;
  vector[p] D;
  cov_matrix[p] Sigma;
}
transformed data {
  real sqrt_nu = sqrt(nu);
  matrix[p,p] L_Sigma = cholesky_decompose(Sigma);
}
parameters {
  vector<lower=0>[p] z;
  real<lower=0> w;
  vector[p] x;
}
transformed parameters {
  vector[p] Y = z .* D + sqrt(w * nu) * L_Sigma * x;
}
model {
  w ~ inv_chi_square(nu);
  z ~ normal(0, 1);
  x ~ normal(0, 1);
}

The main parts of R code are:

p <- nrow(D)
Y <- matrix(0, n, p)  

for(i in 1:n){
    u <- runif(p, 1/2, 1)
    z <- qnorm(u, 0, 1)
    w <- rgamma(1, nu/2, nu/2)
    Y[i,] <- mvrnorm(1, Mu+D%*%z, Sigma/w)
}
Y

I replaced qnorm(runif(10000,0.5,1),0,1) with abs(rnorm(10000)). Changed
gamma to inv_chi_squared.

Can somebody confirm, I did it right? I looks “too easy” for my 2 cents. (I didn’t add the mu)

multivariate_skew_t2a.R (833 Bytes)

Just going by the little blurb in Wikipedia (https://en.wikipedia.org/wiki/Multivariate_t-distribution#Definition), what you’ve written should imply:

z \sim N(0, \text{diag}(D))
Y \sim \text{multivariate_t}(\nu, z, \Sigma)

So it looks right enough to me :P?

I’d check the covariance of your samples against the analytical form on Wikipedia for various sigma and nu first. Just use \nu > 2 or whatever so the covariance is finite. Then you can check your student t against the implementation in that package to make sure they’re in line.

edit: Made Latex prettier