Cool. I didn’t know covariance matrices were flexible enough for that kind of constraint, but I guess I should learn how to think about matrices as bundles of equations.
Your implementation does give me what I need, the scaling. When I apply that to the constrained implementation in Stan,
transformed data {
int<lower = 2> N = 4;
}
parameters {
vector[N - 1] mu_prefix;
}
transformed parameters {
vector[N] mu = append_row(mu_prefix, -sum(mu_prefix));
}
model {
mu ~ normal(0, inv(sqrt(1 - inv(N))));
}
I get the correct marginals:
mean se_mean sd
mu[1] -0.04 0.02 1.03
mu[2] 0.02 0.02 1.02
mu[3] 0.04 0.02 1.01
mu[4] -0.02 0.02 1.01
With more iterations the mean and sd are correct to more decimal places, as expected.
The direct constrained version might be more efficient, because it doesn’t require the matrix multiply A_qr * x_raw
, which is a quadratic operation with lots of autodiff to do for the N x N
Jacobian. On the other hand, it’s less of a normal sampling and reconstruct approach, which should mix very well.