Non-centered parameterization of the multivariate student-t distribution

Modeling
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