Hi,

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:

```
data {
int<lower=0> N; // Number of samples
vector[N] x_meas; // Measured values
vector<lower=0>[N] uncertainty; // Uncertainty in values
}
parameters {
real<lower=0,upper=8> mu;
real<lower=0> sigma;
vector[N] xtrue;
}
model {
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?