Non-centered parameterization of the multivariate student-t distribution

Modeling
1

In the multivariate normal case, we have

mu + L * z

where mu is the mean vector, L is the cholesky factor of the covariance matrix and z are independent standard normal variables. My question is how to correctly incorporate the degrees of freedom parameter nu into the above formula to create multivariate-t random variables? I remember having seens this somewhere already, but I can’t find it right now…

2

https://groups.google.com/d/msg/stan-users/nk03n8hX4Us/-pTN5IGjCQAJ

One may optimize it using the inv_chi_square and turn the division in a multiplication.

3 Likes
3

Thanks, that’s what I was looking for! :-)

4

Hello @paul.buerkner @andre.pfeuffer,

The parameterization: mu+L*z you mentioned for the multivariate normal is valid when the covariance matrix is known (defined in the data block).

Is there any difference when the covariance matrix is unknown (defined in the parameters block)?

Thanks a lot.

5

It’s still valid.

6

I think it’s a bit more complicated in my case. I’ll create a new topic.

7

Could somebody re-post the solution at the end of this thread?

The link to the google group seems to no longer be valid.

8 
data {
  int<lower=1> p;
  vector[p] mu;
  cholesky_factor_cov[p,p] L;
  real<lower=0> nu;
}
parameters {
  vector[p] z;
  real<lower=0> u;
}
transformed parameters {
  vector[p] x;
  x = mu + sqrt(nu / u) * (L * z); // distributed multi_student_t
}
model {
  target += normal_lpdf(z | 0, 1);
  target += chi_square_lpdf(u | nu);
}
3 Likes
9

To consider nu as a parameter rather than data, you can also use this function (derivation below):

real multivariate_t_scale_lpdf(vector x,
                               real nu) {
        int N = size(x);
        vector[N] xx = nu/(x^2);
        real lp = 0.0;
        for (n in 1:N) {
            lp += log(2) + log(nu) + chi_square_lpdf(xx[n] | nu) - 3*log(x[n]);
        }
        return lp;
    }

From some tinkering, I found that the quantile function for the outcome of interest is

y = \sqrt{\frac{\nu}{\text{Q-Chi-Square}(1-p, \nu)}} 
library(tidyverse)

# simulate outcome w/df = 7
nu <- 7
raw <- sqrt(nu/rchisq(1e4, nu))

# quantiles for checking simulated results against empirical
sim <- 
  tibble(x = seq(from = 0, to = 1, length.out = 100),
         q = quantile(raw, probs = seq(from = 0, to = 1, length.out = 100)))

# check the quantile function
tibble(x = seq(from = 0, to = 1, length.out = 1000)) %>%
  mutate(q = sqrt(nu/qchisq(1 - x, nu))) %>%
  ggplot(aes(x = x,
             y = q)) + 
  geom_line() +
  geom_point(data = sim,
             color = "royalblue") +
  theme_minimal()

Created on 2024-05-06 with reprex v2.1.0

Some rearrangement to get p on its own

\begin{align*} \text{Q-Chi-Square}(1-p, \nu) &= \frac{\nu}{y^2} \end{align*}

Set a=\frac{\nu}{y^2} and recalling that the quantile and cumulative functions are inverses:

\begin{align*} \text{CDF-Chi-Square}(a, \nu) &= 1-p \\ p &= 1 - \text{CDF-Chi-Square}(a, \nu) \end{align*}

The pdf is the derivative of the cdf, so by the chain rule, you get the following:

\begin{align*} d &= -\text{PDF-Chi-Square}(a, \nu) \times a' \\ d &= -\text{PDF-Chi-Square}\left(\frac{\nu}{y^2}, \nu \right) \times \frac{-2\nu}{y^3} \\ d &= 2 \times \text{PDF-Chi-Square}\left(\frac{\nu}{y^2}, \nu \right) \times \frac{\nu}{y^3} \end{align*}

The derived density function checks out against the simulated density:

library(tidyverse)

# simulate outcome w/df = 7
nu <- 7
raw <- sqrt(nu/rchisq(1e4, nu))

tibble(x = seq(from = 0, to = 6, length.out = 1000)) %>%
  mutate(d = dchisq(nu/(x^2), nu) * (2*nu)/(x^3)) %>%
  ggplot(aes(x = x)) + 
  geom_line(aes(y = d)) +
  geom_density(data = tibble(x = raw),
               color = "royalblue") +
  theme_minimal()

Created on 2024-05-06 with reprex v2.1.0