I need to add this to the manual. An N-dimensional tensor with variance tau^2 where each dimension sums to 0 through an isotropic sum-to-zero transform will have variance of tau^2 * prod( (d_k - 1) / d_k) so to correct you need tau^2 * prod(d_k / (d_k - 1)) as the correction.
I deliberately coded this to take in the dimensions to be suitable for a sum-to-zero N-d array but it works fine for matrix
functions {
real s2z_unit_normal_correction (array[] int s2z_dims) {
int N = num_elements(s2z_dims);
vector[N] d = to_vector(s2z_dims);
return sqrt(prod(d ./ (d - 1)));
}
real s2z_unit_normal_lpdf (vector z, array[] int s2z_dims) {
int N = num_elements(s2z_dims);
vector[N] d = to_vector(s2z_dims);
return -dot_self(z) + N * 2 * prod(d ./ (d - 1));
}
}
data {
int<lower=1> N;
int<lower=1> M;
int<lower=0, upper=1> flag;
}
parameters {
sum_to_zero_matrix[N, M] X;
}
model {
if (flag == 0) {
to_vector(X) ~ normal(0, s2z_unit_normal_correction(dims(X)));
} else {
to_vector(X) ~ s2z_unit_normal(dims(X));
}
}
Mostly a reminder for @spinkney who said he would add this on slack
Message
I need to add this to the manual. An N-dimensional tensor with variance tau^2 where each dimension sums to 0 through an isotropic sum-to-zero transform will have variance of tau^2 * prod( (d_k - 1) / d_k) so to correct you need tau^2 * prod(d_k / (d_k - 1)) as the correction.
I deliberately coded this to take in the dimensions to be suitable for a sum-to-zero N-d array but it works fine for matrix