I have a problem with a model which is reasonably complicated. However, I've determined that the underlying problem arises in a very simple situation. Say I have data which I create simply by sampling from a normal with mean/sd of mu and sigma. Then I "measure" each variate with some (normal) uncertainty, say u, so that my "observed" sample is
X ~ N(mu,sigma) + N(0,u)
To generate the posterior distributions of mu and sigma, I use a simple hierarchical model as follows:
int<lower=0> N; // Number of samples
vector[N] x_meas; // Measured values
vector<lower=0>[N] uncertainty; // Uncertainty in values
xtrue ~ normal(mu,sigma);
x_meas ~ normal(xtrue,uncertainty);
This works perfectly fine if u is a single constant. However, let's say u is dependent on the intrinsic value, so for instance in the simplest case, we could say
u = u_0 for x < mu,
u = u_1 for x > mu
I would expect the same model to work, since it incorporates individual values for the uncertainty. However, it doesn't seem to. Using mu=4, sigma=2, u_0 = 2 and u_1=0.2 for a sample of 4500, I get the following:
The variety of red lines/regions represent the pdf of the estimated hyperparameters. Clearly, they gravitate strongly towards a Gaussianized version of the observed distribution, rather than the underlying one. Furthermore, the likelihood of the steps in the chain are higher than the likelihood gotten from re-inputing the actual underlying distribution.
So, is this problem fundamentally ill-posed, unsolvable, or am I just doing something stupid?