---
title: "Random effect on standard deviation"
output: html_notebook
---

```{r}
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
```

Stan model

```{stan output.var=model1}
data {
  int<lower = 0> N;    // Number of records
  int<lower = 0> M;    // Number of blocks
  vector[N] Y;         // Data
  int<lower = 1, upper = M> B[N]; // Index of blocks
}
parameters {
  real mu;
  real<lower = 0> sigma_bar;
  vector[M] log_epsilon;
  real<lower = 0> sigma_e;
}
transformed parameters {
  vector<lower = 0>[M] sigma = sigma_bar * exp(log_epsilon);
}
model {
  Y ~ normal(mu, sigma[B]);
  log_epsilon ~ normal(0, sigma_e);
}
```

Generate data

```{r}
set.seed(123)
n_block <- 20
n_block_obs <- 10
mu <- 3
sd <- 1
sd_block <- 0.5
ranef <- rnorm(n_block, 0, sd_block)
stdev <- sd * exp(ranef)
y <- sapply(1:n_block, function(i)
  rnorm(n_block_obs, mu, stdev[i]))
y <- c(y)
block <- rep(1:n_block, each = n_block_obs)
df <- data.frame(block = as.factor(block), y = y)
ggplot(df) +
  geom_point(aes(x = block, y = y))
```

Fit to the simulated data

```{r fit, cache=TRUE}
data <- list(N = length(y), M = n_block, Y = y, B = block)
fit <- sampling(model1, data = data, iter = 2000)
print(fit)
```