I am wondering how Stan handles constraints on correlations. For example, with a 3 by 3 matrix, the correlations cannot be uniform between -1 and 1. Yet, when I fit this size correlation matrix, I do not receive any error unless I include generated quantities which states that a value has been rejected. I take this to mean that a whether the matrix is positive definite is evaluated and then rejected if not before proceeding. And hence, the sampler does not throw an error and stop.
Note this model is for a simple example, which I hope to scale to a 4 by 4 random effects co-variance matrix. Also, while I do know that the LKJ is currently recommended by the Stan team, in my application I need independent priors for each correlation.
model <- '
// Pearson Correlation
data {
int<lower=0> n;
vector[3] x[n];
}
parameters {
vector[3] mu;
vector<lower=0>[3] sigma;
vector<lower=-1,upper=1>[3] r;
}
transformed parameters {
cov_matrix[3] S;
// Reparameterization
S[1,1] = square(sigma[1]);
S[1,2] = r[1] * sigma[1] * sigma[2];
S[2,1] = r[1] * sigma[1] * sigma[2];
S[1,3] = r[2] * sigma[1] * sigma[3];
S[3,1] = r[2] * sigma[3] * sigma[1];
S[2,3] = r[3] * sigma[2] * sigma[3];
S[3,2] = r[3] * sigma[3] * sigma[2];
S[2,2] = square(sigma[2]);
S[3,3] = square(sigma[3]);
}
model {
// Priors
mu ~ normal(0, 100);
sigma[1] ~ cauchy(0,3);
sigma[2] ~ cauchy(0,3);
sigma[3] ~ cauchy(0,3);
r[1] ~ uniform(-1, 1);
r[2] ~ uniform(-1, 1);
r[3] ~ uniform(-1, 1);
// Data
x ~ multi_normal(mu, S);
}
generated quantities{
vector[3] ppc;
ppc = multi_normal_rng(mu, S);
}'
Here is the R code:
x <- scale(mtcars[,1:3])
mod <- stan_model(model_code = model, verbose = TRUE)
fit <- sampling(mod, data=list(x = x, n = nrow(x)))