I am attempting to conduct a robust correlation estimation by adapting This Post to 4 outcome variables (my real data have 9 outcomes).
I’m a little stuck on how to write the likelihood to estimate the additional scales and degrees of freedom parameters.
Some fake data:
sigma = c(10, 15, 10, 20)
rho = c(-0.6, 0.3, 0.6, -0.3, 0.6, -0.3)
cov.mat = matrix(c(sigma[1] ^ 2,
sigma[1] * sigma[2] * rho[1],
sigma[1] * sigma[3] * rho[2],
sigma[1] * sigma[4] * rho[3],
sigma[1] * sigma[2] * rho[1],
sigma[2] ^ 2,
sigma[2] * sigma[3] * rho[4],
sigma[2] * sigma[4] * rho[5],
sigma[1] * sigma[3] * rho[2],
sigma[2] * sigma[3] * rho[4],
sigma[3] ^ 2,
sigma[3] * sigma[4] * rho[6],
sigma[1] * sigma[4] * rho[3],
sigma[2] * sigma[4] * rho[5],
sigma[3] * sigma[4] * rho[6],
sigma[3] ^ 2),
nrow=4, byrow=T)
My current, failing attempt:
robust_code = "
data {
int<lower=1> N; // number of rows
//int<lower=1> K; // number of columns
vector[4] x[N]; // input data: rows (N) are observations, columns (4) are the variables
}
parameters {
vector[4] mu; // locations of the marginal t distributions
real<lower=0> sigma[4]; // scales of the marginal t distributions
real<lower=1> nu[4]; // degrees of freedom of the marginal t distributions
real<lower=-1, upper=1> rho[6]; // correlation coefficients
}
transformed parameters {
// Covariance matrix
cov_matrix[4] cov = [[sigma[1] ^ 2, sigma[1] * sigma[2] * rho[1], sigma[1] * sigma[3] * rho[2], sigma[1] * sigma[4] * rho[3]],
[sigma[1] * sigma[2] * rho[1], sigma[2] ^ 2, sigma[2] * sigma[3] * rho[4], sigma[2] * sigma[4] * rho[5]],
[sigma[1] * sigma[3] * rho[2], sigma[2] * sigma[3] * rho[4], sigma[3] ^ 2, sigma[3] * sigma[4] * rho[6]],
[sigma[1] * sigma[4] * rho[3], sigma[2] * sigma[4] * rho[5],sigma[3] * sigma[4] * rho[6],sigma[3] ^ 2]];
}
model {
// Likelihood
// Bivariate Student's t-distribution instead of normal for robustness
x ~ multi_student_t(nu, mu, cov);
// Noninformative priors on all parameters
sigma ~ normal(0,100);
mu ~ normal(0, 100);
nu ~ gamma(2, 0.1);
}
generated quantities {
// Random samples from the estimated bivariate t-distribution (for assessment of fit)
vector[2] x_rand;
x_rand = multi_student_t_rng(nu, mu, cov);
}"
On compile I get an error: Ill-typed arguments to '~' statement. No distribution 'multi_student_t' was found with the correct signature.
Do I need a loop or some vectorization for the likelihood?