Very simple hierarchical model divergencies - am I missing something obvious?

Modeling
1

Hello Community,

I am fitting a test model, really simple. There are 2 hierarchies of normal distributions. I get divergencies and I don’t understand why. Is there anything obvious I am missing.

The only challenging aspect is that there is little data.

The model


data{
int N;
int K;
vector[N] y;
int idx[N];
}
parameters{
  real mu_mu;
  real<lower=0> mu_sigma;
  
  real<lower=0> sigma_mu;
  real<lower=0> sigma_sigma;
  
  vector<offset=mu_mu, multiplier=mu_sigma>[K] mu;
  vector<lower=0>[K] sigma;
  

}
model{
y ~ normal(mu[idx], sigma[idx]);

mu ~ normal(mu_mu, mu_sigma);
sigma~ normal(sigma_mu, sigma_sigma);

sigma_sigma ~ student_t(3, 0, 5);
mu_sigma ~ student_t(3, 0, 5);
}

Rscript to replicate

library(rstan)
library(tidyverse)

mod = "
data{
int N;
int K;
vector[N] y;
int idx[N];
}
parameters{
  real mu_mu;
  real<lower=0> mu_sigma;
  
  real<lower=0> sigma_mu;
  real<lower=0> sigma_sigma;
  
  vector<offset=mu_mu, multiplier=mu_sigma>[K] mu;
  vector<lower=0>[K] sigma;
  

}
model{
y ~ normal(mu[idx], sigma[idx]);

mu ~ normal(mu_mu, mu_sigma);
sigma~ normal(sigma_mu, sigma_sigma);

sigma_sigma ~ student_t(3, 0, 5);
mu_sigma ~ student_t(3, 0, 5);
}
"
my_mod = stan_model(model_code = mod)

.data = 
  tibble(idx = c(1,2), .mean = c(1, 2), .sigma = c(1, 1)) %>%
  mutate(y = map2(.mean, .sigma, ~ rnorm(100, .x, .y))) %>%
  unnest(y) 


fit = sampling(my_mod, data = list(N=nrow(.data), y = .data$y, idx = .data$idx, K=(.data$idx %>% unique %>% length)))

image
image1180×618 112 KB

Although you see a funnel, I am using non-centered parametrization (trough offset and multiplier)

1 Like
2

As a first step in debugging, would it help if you used std_normal() here first? (and comment out hyperparameters for sigma)

sigma ~ std_normal()