# Generate from a normal distribution in model block

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!

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

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

1 Like

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

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

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

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

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);
}