Generate from a normal distribution in model block

Modeling
1

Hi there,
I’m trying to fit a model that looks like this:
\hat V = Normal(V, \sigma_v)
\Delta EV = P_L {\hat V _L}^{\rho} - {\hat V _{SB}}^{\rho}
p(choose\ lottery) = ({1 + e^{-\Delta EV}})^{-1}

V is a vector of lottery values that I pass in, \sigma_v and \rho are the parameters. The problem is, there is no function in stan that allows to me to generate \hat V from V in the model block. normal_rng can only be used in generated quantities. Is there an easy solution to this? Many thanks for any help!

2

Hey! Did you try something like this

...
parameters{
  ...
  vector[N] v_hat;
  ...
}
model{
  v_hat ~ normal(v, sigma_v);
  ...
}

?

Edit: Or something like

...
parameters{
  ...
  vector[N] z;
  ...
}
transformed parameters{
  vector[N] v_hat = v + sigma_v * z;
model{
  z ~ std_normal();
  ...
}
1 Like
3

Hi Max, thanks for the great tip! The second option works for me.

1 Like
4

Actually, now I realized this approach has one problem: it will save iter * N parameter z, which is not practical when fitting on a large dataset. My final goal is to fit this model on a large dataset with more than 1M trials. Is there any other way to implement this without declaring extra parameters? @Max_Mantei

5

The problem is that \hat V hast to be a parameter (you can only use RNG functions inside the generated quantities and transformed data blocks in Stan). In the second approach you can move the line vector[N] v_hat = v + sigma_v * z; inside the model block, so it’ll not save v_hat. Also, you can specify which parameters to save in your specific Stan interface (you’ll still need the memory while fitting).

2 Likes
6

Alternatively, you can use the offset and multiplier syntax. Where this:

...
parameters{
  ...
  vector[N] z;
  ...
}
transformed parameters{
  vector[N] v_hat = v + sigma_v * z;
model{
  z ~ std_normal();
  ...
}

is equivalent to:

...
parameters{
  ...
  vector<offset=v, multiplier=sigma_v>[N] v_hat;
  ...
}
model{
  v_hat ~ normal(v, sigma_v);
  ...
}

That way you only save the values of v_hat, but the estimation is the same

4 Likes
7

I always forget about this… Good catch @andrjohns! Thank you!

8

It took me a deep dive into stan manual to find this, thanks for the insight

9

I’m not sure if I should start a new thread for this, but I’m gonna post it here anyway. Basically here is a hierarchical version of the simple model. The simple model worked fine on synthetic data. However, when running this model on synthetic population data, the fits end up with > 80% divergent transitions and large Rhats. Increasing adapt_delta and increasing iter from 2000 to 3000 didn’t seem to help much, so didn’t increasing dataset size.
My intuition tells that I should reparameterize the model, but honestly I don’t know how. I’m quite new to stan so any help is greatly appreciated!!

// hierarchical model fitting one species as rho-sigmav agents
functions {
  real choicex(real rho, real lott_value_hat, real lott_prob, real sb_value_hat, real rew_multi) {
    real y;
    real u1;	// Lottery utility
    real u2;	// Surebet utility
    
    u1 = lott_prob * (rew_multi * lott_value_hat) ^ rho;
    u2 = (rew_multi * sb_value_hat) ^ rho;
    y = u1 - u2; 
    return y;
  }
}
data {
  int<lower=0> N; // Number of trials we have
  int<lower=0> K; // number of subjects in each species
  int individual[N]; // vector of subjid indexes
  vector[N] lott_value; // Lottery value for that choice
  vector[N] lott_prob; // Lottery probabilities for that choice
  vector[N] sb_value; // Surebet values for that choice
  vector[N] total_rew_multi; // total reward multiplier = base_reward * rew_multi
  int<lower=0,upper=1> y[N]; // Choices we observe (1 if they pick lottery)
}
parameters {
  real<lower=0, upper=4> rho_s; // species-level rho
  real<lower=0> sigmav_s; // species-level sigma_v
  real<lower=0> sigma_rho; // standard deviation for species rho
  real<lower=0> sigma_sigmav; // standard deviation for species sigma_v
  vector<lower=0, upper=4>[K] rho_i; // individual-level rho
  vector<lower=0> [K] sigmav_i; // individual-level sigma_v
  vector[N] z_lott; // z score of lott_value_hat
  vector[N] z_sb; // z score of sb_value_hat
}
model {
  vector[N] thetas;
  real lott_value_hat; // placeholder for perceived lottery value
  real sb_value_hat; // placeholder for perceived surebet value
  
  // set weak priors 
  rho_s ~ normal(1, 0.5); 
  sigmav_s ~ normal(3, 1); 
  sigma_rho ~ normal(0.5, 0.3);
  sigma_sigmav ~ normal(0.5, 0.3);
  z_lott ~ std_normal(); // implies lott_value_hat ~ normal(lott_value, sigma_v * z_lott)
  z_sb ~ std_normal(); // implies sb_value_hat ~ normal(sb_value, sigma_v * z_sb)

  // draw individual parameters from the species & population distribution
  for(k in 1:K){
    rho_i[k] ~ normal(rho_s, sigma_rho);
    sigmav_i[k] ~ normal(sigmav_s, sigma_sigmav);
  }
  
  // fit the actual model to each trial
  for(n in 1:N){
    lott_value_hat = lott_value[n] + sigmav_i[individual[n]] * z_lott[n]; 
    sb_value_hat = sb_value[n] + sigmav_i[individual[n]] * z_sb[n];
    thetas[n] = choicex(rho_i[individual[n]], lott_value_hat, lott_prob[n], sb_value_hat, total_rew_multi[n]);
  }
  y ~ bernoulli_logit(thetas);
}