New to Stan/Bayesian modelling, apologies if this is too naive or not appropriate for these forums.

I have a somewhat novel model (I think), for which I obtain two measurements per object, `nmeas`

and `zeta`

. Across the population of objects there is an inherent `ninh`

that is a (tight) Gaussian, with mean `n_mu`

which I wish to estimate. I know that `ninh`

is strictly between 2 and 7. The model is that

`nmeas`

~ normal(`ninh`

,`zeta`

)

I generate some synthetic data in python via

```
import numpy as np
rng = np.random.default_rng(seed=0) #for reproducibility
ninh = rng.normal(3.5, 0.1, size=500) #assume n_mu=3.5
zeta = rng.lognormal(1, 0.5, size=500)
nmeas = rng.normal(ninh, zeta, size=500)
data_dic = {"NN": 500,
"nmeas": nmeas,
"zeta": zeta,
}
```

I can then recover `ninh`

(`n_mu`

) using the very simple Stan model

```
data {
int NN;
vector[NN] nmeas;
vector[NN] zeta;
}
parameters {
real<lower=2, upper=7> n_mu;
real<lower=2, upper=7> ninh;
}
model {
ninh ~ normal(n_mu, 0.1);
nmeas ~ normal(ninh, zeta);
}
```

While this is all well-and-good, in reality the expected distribution of `zeta`

across the population is a lognormal covering 12 orders of magnitude. That is, I need to make the above work with

```
zeta = rng.lognormal(7, 7, size=500)
```

At this point Stan rightfully has some nasty convergence complaints for me, and fails to recover the true value of `n_mu`

much of the time.

I believe the answer will lie in a smart re-parameterization of the model, but my initial attempts haven’t come up with much success. Any help would be appreciated!