Dear Stan experts,
I am trying to fit a time-varying drift diffusion model(tDDM) adopted from Maier et al., 2020 NHB with cmdstanpy. However, there is a very weird thing. All parameters in generated quantities have no value. Is something wrong with my Stan code? Thanks.
data {
int<lower=1> ns;
int<lower=1> nt; // number of trials
vector [nt] choice[ns];
vector [nt] ss_money[ns];
vector [nt] ss_time[ns];
vector [nt] ll_money[ns];
vector [nt] ll_time[ns];
real<lower=0> rt[ns,nt]; // response time
vector<lower=0>[ns] minRT; //minimal RT per subject
}
parameters {
// Trial level parameters
real w_mu; // weight for the reward atttribute
real inter_mu;
real s_mu; // inverse temperature in softmax
real b_mu; //boundary
real z_mu; //starting point
real ndt_mu; // non-decision time
real d_mu; // scaling parameter for drift rate
real rst_mu; // relative starting time
real<lower=0> w_sigma;
real<lower=0> inter_sigma;
real<lower=0> s_sigma;
real<lower=0> b_sigma;
real<lower=0> z_sigma;
real<lower=0> ndt_sigma;
real<lower=0> d_sigma;
real<lower=0> rst_sigma;
// Indiviual level raw parameters
vector[ns] w_raw;
vector[ns] inter_raw;
vector[ns] s_raw;
vector[ns] b_raw;
vector[ns] z_raw;
vector[ns] ndt_raw;
vector[ns] d_raw;
vector[ns] rst_raw;
}
transformed parameters{
real w[ns]; //dicount rate in hyperbolic function
real inter[ns];
real rst[ns];
real s[ns];
real d[ns];
real b[ns];
real z[ns];//parameter that govern the maxmium value of drifte rate //bias parameter for after stress trial
real ndt[ns]; //none decision time parameter
for (i in 1:ns){
w[i] = exp(w_mu + w_sigma * w_raw[i]);
inter[i] = inter_mu + inter_sigma * inter_raw[i];
rst[i] = 2*Phi_approx(rst_mu + rst_sigma * rst_raw[i])*(minRT[i]) - minRT[i]; // range of rst is [-minRT ,minRT]
s[i] = exp(s_mu + s_sigma * s_raw[i]);
d[i] = exp(d_mu + d_sigma * d_raw[i]);
b[i] = exp(b_mu + b_sigma * b_raw[i]);
ndt[i] = Phi_approx(ndt_mu + ndt_sigma * ndt_raw[i])*(minRT[i]-0.1)+0.1;
z[i] = Phi_approx(z_mu + z_sigma * z_raw[i]);
}
}
model {
// group level hyerparameters for prior distribution
w_mu ~ normal(0,1);
rst_mu ~ normal(0,1);
s_mu ~ normal(0,1);
d_mu ~ normal(0,1);
b_mu ~ normal(0,1);
ndt_mu ~ normal(0,1);
z_mu ~ std_normal();
inter_mu ~ std_normal();
w_sigma ~ normal(0,3);
rst_sigma ~ normal(0,3);
s_sigma ~ normal(0,3);
d_sigma ~ normal(0,3);
b_sigma ~ normal(0,3);
ndt_sigma ~ normal(0,3);
z_sigma ~ normal(0,3);
inter_sigma ~ normal(0,3);
// individual parameters
w_raw ~ normal(0,1);
rst_raw ~std_normal();
s_raw ~ normal(0,1);
d_raw ~ normal(0,1);
b_raw ~ normal(0,1);
ndt_raw ~ normal(0,1);
z_raw ~ std_normal();
inter_raw ~ std_normal();
// subject loop
for (i in 1:ns){
real sign;
real w_sv; // weighted subjective value difference
real value_diff;
real time_diff;
real drift;
real choice_map;
real z_map;
// sign of the rst parameter
if (rst[i]>0){
sign = 1;
}
else{
sign = 0; // 1 for later -larger enter evidence accumulation first
}
for(t in 1:nt){
if (rt[i,t] == 100) continue;
choice_map = (choice[i,t] - 1) * 2 - 1; // 1 for ll, -1 for ss
value_diff = (ll_money[i,t] - ss_money[i,t]);
time_diff = (ll_time[i,t] - ss_time[i,t]);
w_sv = (inter[i] * rt[i,t] + w[i] * value_diff * (rt[i,t] - abs(rst[i]) * ( 1 - sign)) - time_diff * (rt[i,t] - abs(rst[i]) * sign))/rt[i,t] * s[i];
// calculate1 averaged drift rate
drift = d[i] * choice_map * (2 *inv_logit(w_sv) - 1);
z_map = z[i] * (choice[i,t] - 1) + (1 - z[i]) * abs(choice[i,t] - 2);
target += wiener_lpdf(rt[i,t] | b[i],ndt[i],z_map,drift) + machine_precision();
} // end trial loop
} // end subject loop
}
generated quantities{
matrix[ns,nt] log_lik;
// subject loop
for (i in 1:ns){
real sign;
real w_sv; // weighted subjective value difference
real value_diff;
real time_diff;
real drift;
real choice_map;
real z_map;
// sign of the rst parameter
if (rst[i]>0){
sign = 1;
}
else{
sign = 0; // 1 for later -larger enter evidence accumulation first
}
for(t in 1:nt){
log_lik[i,t] = 0;
if (rt[i,t] == 100) continue;
choice_map = (choice[i,t] - 1) * 2 - 1; // 1 for ll, -1 for ss
value_diff = (ll_money[i,t] - ss_money[i,t]);
time_diff = (ll_time[i,t] - ss_time[i,t]);
w_sv = (inter[i] * rt[i,t] + w[i] * value_diff * (rt[i,t] - abs(rst[i]) * ( 1 - sign)) - time_diff * (rt[i,t] - abs(rst[i]) * sign))/rt[i,t] * s[i];
z_map = z[i] * (choice[i,t] - 1) + (1 - z[i]) * abs(choice[i,t] - 2);
log_lik[i,t] = wiener_lpdf(rt[i,t] | b[i],ndt[i],z_map,drift);
} // end trial loop
} // end subject loop
}