Variance reparametrization

variance_reparametrization.Rmd

```{r}
library(cmdstanr)
```

```{r}
intercept_uncentered <- cmdstan_model(write_stan_file("data {
  int<lower=1> N;
  vector[N] signal;
  int<lower = 1> N_subj;
  vector[N] c_cloze;
  array[N] int<lower = 1, upper = N_subj> subj;
}
parameters {
  real<lower = 0> sigma;
  real<lower = 0>  tau;
  real alpha;
  real beta;
  vector[N_subj] z_u;
}
transformed parameters {
  vector[N_subj] u;
  u = z_u * tau;
 }
model {
  target += normal_lpdf(alpha| 0,10);
  target += normal_lpdf(beta | 0,10);
  target += normal_lpdf(sigma | 0, 50)  -
    normal_lccdf(0 | 0, 50);
  target += normal_lpdf(tau | 0, 20)  -
    normal_lccdf(0 | 0, 20);
  target += std_normal_lpdf(z_u);
  target += normal_lpdf(signal | alpha + u[subj] +
                        c_cloze .* beta, sigma);
}"))
```

```{r}
set.seed(455665)
dd <- list(N = 20, signal = rnorm(20),
           N_subj = 18,
           c_cloze = runif(20) - 0.5,
           subj = c(1,1,2,2,3:18))

res_orig <- intercept_uncentered$sample(dd)
res_orig
```

```{r}
bayesplot::mcmc_pairs(res_orig$draws(), pars = c("sigma", "tau", "alpha", "beta"))
bayesplot::mcmc_pairs(res_orig$draws(), pars = c("sigma", "tau"), transformations = "log")
```


```{r}
intercept_uncentered_split <- cmdstan_model(write_stan_file("data {
  int<lower=1> N;
  vector[N] signal;
  int<lower = 1> N_subj;
  vector[N] c_cloze;
  array[N] int<lower = 1, upper = N_subj> subj;
}
parameters {
  real<lower=0> sigma_total;
  real<lower=0,upper=1> sigma_split;
  real alpha;
  real beta;
  vector[N_subj] z_u;
}
transformed parameters {
  real<lower = 0> sigma;
  real<lower = 0>  tau;
  vector[N_subj] u;
  sigma = sqrt(sigma_total^2 * sigma_split);
  tau = sqrt(sigma_total^2 * (1 - sigma_split));
  u = z_u * tau;
 }
model {
  target += normal_lpdf(alpha| 0,10);
  target += normal_lpdf(beta | 0,10);
  target += normal_lpdf(sigma_total | 0, sqrt(50 ^ 2 + 20^2))  -
    normal_lccdf(0 | 0, sqrt(50 ^ 2 + 20^2));
  target += std_normal_lpdf(z_u);
  target += normal_lpdf(signal | alpha + u[subj] +
                        c_cloze .* beta, sigma);
}"))
```

```{r}
res_split <- intercept_uncentered_split$sample(dd)
res_split
```

```{r}
bayesplot::mcmc_pairs(res_split$draws(), c("sigma_total", "sigma_split", "alpha", "beta"))
bayesplot::mcmc_pairs(res_split$draws(), c("sigma_total", "sigma_split", "alpha", "beta"),
transformations = list("sigma_total" = "log", "sigma_split" = "qlogis", "alpha" = "identity", "beta" =  "identity"))
```