# Hierarchical Model - Hyper priors for a Gamma Distribution

I am working on building a hierarchical RL model for a parameter drawn from a gamma distribution (beta_1 variable). I am attempting to re-parameterize the gamma distribution in terms of its mean and variance and then set the hyperpriors from a uniform (0,20) but am getting an insane number of divergences so am definitely just doing something wrong. Could anyone point me in the right direction?

I am new to heirarchical modeling so any help will be greatly appreciated!

data {
int Trial;
int outcome[Trial,5];
int choice[Trial,5];
int stim1[Trial,5];
int stim2[Trial,5];
int feed[Trial,5];
int  times[Trial,5];
int<lower=0, upper=1> current_choice[Trial,5] ;

}

transformed data {

matrix[Trial,5] time_passed;
for(j in 1:5){
for(i in 2:Trial){
time_passed[i,j]=(times[i,j]-times[i-1,j])/(24*3600*1000);
}}
}

parameters {
// hyperparameters for learning rate
real<lower=0,upper=1>  w;
real<lower=0> kap;
// and beta
real<lower=0> m_b;
real<lower=0> s_b;

//per subject parameter
vector<lower=0,upper=1>[5] alpha_1;
vector<lower=0>[5] beta_1;
real  lambda_1;

}

transformed parameters {
real v_b= square(s_b);
real shape_b=square(m_b)/v_b;
real rate_b= m_b/v_b;
}

model {
vector[2] current;
real a;
real b;
real k;
// hyper prior for the learning rate - i.e. the overall learning rate
w ~ beta(1,9); // highest level priors
kap ~ gamma(.01,.01); //highest level priors - alpha didnt have a problem
k= kap+2;
a=w*(k-2)+1;
b=(1-w)*(k-2)+1;
//uniform higher level priors for beta_1 variable
m_b ~ uniform(1,20);
s_b ~ uniform(1,20);

for(n in 1:5){
vector [Trial] PE;
vector[105] ev;
alpha_1[n] ~ beta(a,b);
PE=rep_vector(0,Trial);
ev=rep_vector(0,105);
beta_1[n] ~ gamma(shape_b,rate_b);
lambda_1[n] ~ gamma(3,1);

for (t in 1:Trial){
if(t>1){

ev = (.5-ev)* exp(-lambda_1[n]*time_passed[t,n])+.5;

}
current[1] = ev[stim1[t,n]+1];
current[2] = ev[stim2[t,n]+1];

current_choice[t,n] ~ bernoulli_logit(beta_1[n]*current[2] - beta_1[n]*current[1]);

if(feed[t,n]==1){
PE[t] = outcome[t,n] - ev[choice[t,n]+1];
ev[choice[t,n]+1] +=  alpha_1[n] * PE[t];
}

}
}
}

1 Like

Before we start debugging the whole model, I’d like to make two requests:

i) Can you write a simpler version of your model that still displays the pathological behaviour?
ii) What do Gamma random variables drawn from this prior specification look like?
In other words, if I write
\mu_b \sim \operatorname{Uniform}(1, 20),
v_b \sim \operatorname{Uniform}(1, 20),
Y \sim \operatorname{Gamma}(\alpha(\mu_b, v_b), \beta(\mu_b, v_b)),
what does Y look like?
Finally, it’s important to note that the same way reparametrisation can make geometry nicer, it can also make it tougher. If you’re particularly unlucky, the reparametrisation with the most interpretable parameters will also be the one with worse geometry.

2 Likes

One more thing, maybe @maxbiostat can help to clear out.

The gamma distribution is not a location scale model, what I mean is that alpha hyper-parameter doesn’t behave like a mean.

So maybe re-parametrize a gamma in terms of a mean and variance is not a good idea

1 Like

Hi all, would try reparameterising a gamma in terms of its mean and shape (not variance) as the variance scales with the mean (given shape).

1 Like

Hi, have you made any progress with this issue? I’m curious! And good for posterity to report back ;)