I was mostly referencing this post from that thread, which described the singular covariance matrix that fulfills that constraint.

Following @andre.pfeuffer’s suggestion, which is much more stable than mine, here is a worked example. I also stand corrected about the scaling on the variances, it depends if you care about the marginal variance or the variance of the parameters.

So in your example, where the last parameter is just the negative sum of the previous parameters, it does not have the right variance. But if you pick the right transformation, you can have N-1 free parameters, put the prior on those, and end up with N parameters with the right distribution.

```
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
x[1] 0.00 0.02 1.01 -2.02 -0.68 0.00 0.68 1.99 4000 1
x[2] -0.02 0.02 1.00 -1.96 -0.70 0.00 0.66 1.93 4000 1
x[3] 0.00 0.02 0.98 -1.90 -0.66 -0.01 0.65 1.91 4000 1
x[4] 0.00 0.02 0.99 -1.97 -0.66 0.00 0.66 1.96 4000 1
x[5] 0.01 0.02 0.99 -1.91 -0.65 0.02 0.66 1.96 4000 1
x[6] 0.00 0.02 1.01 -1.99 -0.66 -0.01 0.66 2.03 4000 1
x[7] -0.01 0.02 1.01 -2.01 -0.69 -0.02 0.67 1.98 4000 1
x[8] 0.02 0.02 1.00 -1.95 -0.63 0.05 0.70 1.92 4000 1
x[9] 0.00 0.02 1.00 -1.96 -0.66 -0.01 0.66 1.94 4000 1
sum_x 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4000 1
var_x 1.12 0.01 0.56 0.32 0.70 1.03 1.42 2.44 1987 1
```

```
transformed data{
int N = 9;
matrix [N,N] A = diag_matrix(rep_vector(1,N));
matrix [N,N-1] A_qr;
for(i in 1:N-1)
A[N,i] = -1;
A[N,N] = 0;
A_qr = qr_Q(A)[,1:(N-1)];
}
parameters {
vector [N-1] x_raw;
}
transformed parameters{
vector [N] x = A_qr * x_raw;
}
model {
x_raw ~ normal(0,1/sqrt(1-inv(N)));
}
generated quantities{
real sum_x = sum(x);
real var_x = variance(x);
}
```