Greetings, I am trying to put together an errors-in-variables linear regression model in stan (as a starting point of its censored version). The model I have, even though it estimates the slope and intercept correctly, always understimates the error sigma (I have checked that it is not a confusion between variance and sd…). The initial model I tried is below:

```
data{
int n;
real xobs[n];
real yobs[n];
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
real x[n];
}
model {
// priors
alpha ~ normal(0, 5);
beta ~ normal(0, 5);
sigma ~ normal(0, 5);
// model structure
for (i in 1:n){
xobs[i] ~ normal(x[i], sigma);
}
for (i in 1:n){
yobs[i] ~ normal(alpha + beta*x[i], sigma);
}
}
```

And I have run this in R and there I generate the data via (I set intercept to 0 for testing purposes but still underestimates sigma)

```
n <- 1000
slope <- 1
intercept <- 0
sigma <- 1
xobs <- runif(n, -3, 3) + rnorm(n,sigma)
yobs <- slope*xobs + intercept + rnorm(n,sigma)
```

I get fitted sigma values of around 0.7 all the time. When I increase the real sigma to 2 then it is around 1.4.

I have also tried a more complicated model (inspired by stan users guide measurement error model which is

```
data {
int n;
real xobs[n];
real yobs[n];
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
vector[n] x;
real mu_x;
real<lower=0> sigma_x;
}
model {
x ~ normal(mu_x, sigma_x); // prior
xobs ~ normal(x, sigma); // measurement model
yobs ~ normal(alpha + beta * x, sigma);
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ cauchy(0, 5);
sigma_x ~ cauchy(0, 5);
mu_x ~ normal(0, 10);
}
```

This one also produces the similar underestimated sigmas. I have tried it for slope = 2, intercept = -1.5 and sigma= 1 and the mean value of the posterior samples (which look nice as a histogram) come out around 2, -0.49 and 0.46.

For people who are familiar with pymc3, I had previously done the second model above in pymc3 (should be more or less identical) which does give the correct sd. And I am super confused because to me the second model above and the model below are identical (apart from using HalfNormal instead of Cauchy for sigmas). And help would be much appreciated.

```
with pm.Model() as model:
alpha = pm.Normal("intercept", 0., 5.)
beta = pm.Normal("slope", 0., 5.)
sigma_x = 0.01 + pm.HalfNormal("sigma_x", 5.)
sigma = 0.01 + pm.HalfNormal("sigma", 5.)
mu_x = pm.Normal("x_mu", 0., 10.)
x = pm.Normal("x", mu_x, sigma_x, shape=x_obs.shape)
y = pm.Deterministic("y_pred", alpha + beta * x)
x_likelihood = pm.Normal("x_obs", mu=x, sigma=sigma, observed=x_obs,
shape=x_obs.shape)
y_likelihood = pm.Normal("y", mu=y, sd=sigma, observed=y_obs,
shape=y_obs.shape)
trace = pm.sample(400, tune=1000, chains=6, cores=6, target_accept=0.95,
return_inferencedata=False, init='jitter+adapt_diag')
```