Let’s imagine a scenario where a binary variable is drawn from a vector of probabilities that we can convert to the log-odds scale:

```
library(boot)
library(tidyverse)
library(rstan)
set.seed(0)
N <- 10000
logodds_true <- rnorm(N)
p_true <- inv.logit(logodds_true)
Y <- rbinom(N, 1, p_true)
```

We come up with a measured log-odds for each outcome that contains some error:

```
sigma_error <- 0.5
logodds_meas <- logodds_true + rnorm(N, sd = sigma_error)
```

I’ve tried to write a model that mimics this data generating process:

```
data {
int<lower = 0> N;
vector[N] logodds_meas;
int<lower = 0, upper = 1> Y[N];
}
parameters {
vector[N] logodds_true;
real<lower = 0> sigma_error;
}
model {
// prior not super important -- going with half normal
sigma_error ~ normal(0, 1);
// noisy measurements of true log-odds
(logodds_meas - logodds_true) ~ normal(0, sigma_error);
Y ~ bernoulli_logit(logodds_true);
}
```

But for some reason, I can’t recover `sigma_error`

– the posterior always gives back far too large of an estimate, even with excessively high `N`

. For example, the posterior mean in the above example is around 1.2, whereas the actual standard deviation of the error was 0.5.

Anyone know what I’m doing wrong here?