Does anyone know how to model an individual specific variance of the error-term in the simple Radon example? So, the model is

y_i = \alpha_{j[i]} + \beta x_{i} + \epsilon_i, with \epsilon_{i} \sim N(0,\sigma^2_{i}). The standard stan code with \epsilon_{i} \sim N(0,\sigma^2) equals:

```
data {
int<lower=0> J;
int<lower=0> N;
int<lower=1,upper=J> county[N];
vector[N] x;
vector[N] y;
}
parameters {
vector[J] a;
real b;
real mu_a;
real<lower=0,upper=100> sigma_a;
real<lower=0,upper=100> sigma;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- a[county[i]] + x[i] * b;
}
model {
sigma_a ~ uniform(0, 100);
a ~ normal (mu_a, sigma_a);
b ~ normal (0, 1);
sigma ~ uniform(0, 100);
y ~ normal(y_hat, sigma);
}
```

My next question would be how to do this in a three-level model, for example \epsilon would then follow \epsilon \sim N(0,\sigma^2_{ij}).