Hi,
sorry for taking quite long to respond.
One thing I’ve found useful in a different context where I also had two sources of variability that were hard to distinguish was to have a parameter for total variability and then a proportion of the variability to be accounted for by one of the components, e.g. (for just one sd parameter):
parameters {
real<lower=0> total_sd;
real<lower=0,upper=1> varianceG_frac;
}
transformed parameters {
//simplified from sqrt((total_sd^ 2) * varianceG_frac)
real<lower=0> sdG = total_sd * sqrt(varianceG_frac);
//simplified sqrt((total_sd^ 2) * (1 - varianceG_frac))
real<lower=0> sdE = total_sd * sqrt(1 - varianceG_frac);
}
Note that we need to transform from sd to variance, as variance is a linear operator while sd is not.
This way total_sd becomes usually well identified and independent of varianceG_frac (in the model I’ve used the posterior for this parameter was basically uniform over [0,1] indicating that we couldn’t learn the decomposition, so the model is effectively integrating it out - and the sampling worked quite well).
Does that make sense? And do you think this answers you inquiry?
Not sure if you could decompose the whole covariance matrix this way, but it might be possible… (my linear algebra is too rusty for that). Or maybe just doing this for the standard deviations would be enough.
Best of luck with your model!