Scaling up a hierarchical model


Thanks @stijn! I will try.


You could use lambda[i] ~ gamma(mu * inv_sigma, inv_sigma),
with inv_sigma = 1 / sigma.

mean(rgamma(10000,1 * 5, 5))
[1] 0.9961446
var(rgamma(10000,1 * 5, 5))
[1] 0.2008218


First, you need to use a non-centered parameterization of scale, which won’t speed up iterations but will greatly improve effective sample size.

parameters {
  vector[M] scale_std:

transformed parameters {
  vector<lower = 0>[M] scale = exp(mu + scale_std * sigma);
model {
  scale_std ~ normal(0, 1);

You can compound declare define, so that would be

vector<lower = 0>[M] rate = inv(scale);

We dont’ recommend diffuse priors like these:

shape ~ gamma(0.01, 0.01);

Everything should be much better behaved with something with thinner tails and more mass concentrated around where you think the values should be.

Also, the cauchy(0, 5) on sigma is very wide compared to what sigma probably is.


Many thanks @Bob_Carpenter! I have a few quick questions.

Q1. Is the following not good?

parameters {
  vector<lower=0>[M] scale_std:
transformed parameters {
  vector<lower=0>[M] scale = mu + scale_std * sigma;

Q2. When we do non-centered param, what hyperprior distribution should I use for mu and sigma? It seems that I may leave them undefined?

Q3. I cannot use matrix for discrete distributions, correct?


A1: that’s fine.

A2: if you don’t specify a prior for a variable, it’s implicitly uniform on the values that satisfy its constraints. We recommend at least weakly informative priors everywhere.

A3: Not sure what you mean by matrix for discrete distributions. It’s not fully vectorized, if that’s what you mean. The index of the manual has the legal instantiations including vectorizations (the ones that end in s).