Multivariate Likelihoods

Modeling
7

Nevermind! I got it to work! Code is below. I wasn’t postmultiplying where I should’ve been.

I was able to estimate a multivariate normal distribution without actually using the multivariate normal distribution.
It uses the logic of Y = Mu + SLZ. Essentially, it is finding a mu and cholesky-covariance matrix such that when applied to the Y matrix, yields uncorrelated standardized variates. This is the inversion of the process of sampling from a multivariate normal without actually using the multivariate density function, so to speak.

Given Y, find Mu and SL such that (Y - Mu)/(SL)’ ~ Normal(0,1). This should let people use multivariate forms fairly easily, so long as one can express a multivariate distribution as M + SLZ, which /is/ the case for a skew-normal according to Azzalini (and skew-T for that matter).
Beautiful.

 data{
  int N;
  int J;
  matrix[N,J] y;
}
parameters {
  vector[J] xi;
  vector<lower=0>[J] S;
  cholesky_factor_corr[J] L;
}
model {
  matrix[N,J] x;
  for(n in 1:N){
    x[n] = (y[n] - xi') / diag_pre_multiply(S,L)';
  }
  for(n in 1:N){
    x[n] ~ normal(0,1);
    target += log_determinant(inverse(diag_pre_multiply(S,L)));
  }
}

generated quantities{
  corr_matrix[J] R = multiply_lower_tri_self_transpose(L);
  cov_matrix[J] V = quad_form_diag(R,S);
}