Dear Stan/Bayesian model experts,
I am using a prospect theory + drift diffusion model in Stan. I fit the model with 4 MCMC chains. The weird thing is one draw’s likelihood computed in the generated quantities block is NaN value, however, if i extract the parameter value of that draw and manually compute the likelihood in R instead of in stan, I found the log-likelihood value is normal (not NaN value). And all the other draws do not have this error. Here is the stan code:
When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use
functions {
real ddm_(real rt, real response, real v, real A, real z, real ndt){
if (response == 2){
return(wiener_lpdf(rt|A,ndt,z,v));
}
else{
return(wiener_lpdf(rt|A,ndt,1-z,-v));
}
}
}
data {
int<lower=1> ns; // subject number
int<lower=1> nt; // trial number
int value[ns, nt];
int ratio[ns, nt];
int choice[ns,nt];
int sign[ns,nt];
real rt[ns, nt];
vector[ns] minrt;
}
parameters {
vector[7] mu_pr; // group-level mean hyper parameter
vector<lower=0>[7] sigma; // group-level standard deviation
vector[ns] alpha_raw;
vector[ns] tau_raw;
vector[ns] eta_raw;
vector[ns] A_raw;
vector[ns] z_g_raw;
vector[ns] z_l_raw;
vector[ns] t0_raw;
}
transformed parameters {
vector[ns] alpha;
vector[ns] tau;
vector[ns] eta;
vector[ns] A;
vector[ns] t0;
vector[ns] z_g;
vector[ns] z_l;
real drift[ns,nt];
real z_t[ns,nt];
real utility;
// matt-trick
alpha = exp(mu_pr[1] + sigma[1] * alpha_raw);
tau = exp(mu_pr[2] + sigma[2] * tau_raw);
eta = exp(mu_pr[3] + sigma[3] * eta_raw);
A = exp(mu_pr[4] + sigma[4] * A_raw);
z_g = inv_logit(mu_pr[5] + sigma[5] * z_g_raw);
z_l = inv_logit(mu_pr[6] + sigma[6] * z_l_raw);
t0 = inv_logit(mu_pr[7] + sigma[7] * t0_raw) .* minrt;// avoid non-decision time > rt
// subject loop
for (i in 1:ns) {
// trial loop
for (t in 1:nt) {
utility = .5 * sign[i,t] * pow((value[i,t] * ratio[i,t]), alpha[i]) -
sign[i,t] * pow(value[i,t],alpha[i]);
drift[i,t] = eta[i] * (2 * inv_logit(tau[i]*utility) - 1);
z_t[i,t] = (sign[i,t]+1)/2 * z_g[i] - (sign[i,t]-1)/2 *z_l[i];
} // end of trial loop
} // end of subject loop
}
model {
// hyper parameters
mu_pr ~ std_normal();
// Gelman recommend if the number of individual level is small, we should adopt a prior with a thinner tail like half-normal
sigma ~ normal(0,3);
//sigma ~ cauchy(0,3);
// random effect parameters
alpha_raw ~ std_normal();
tau_raw ~ std_normal();
eta_raw ~ std_normal();
A_raw ~ std_normal();
z_g_raw ~std_normal();
z_l_raw ~std_normal();
t0_raw ~ std_normal();
for(i in 1:ns){
for(t in 1:nt){
target += ddm_(rt[i,t],choice[i,t],drift[i,t],A[i],z_t[i,t],t0[i]);
}
}
}
generated quantities {
real log_lik[ns]; ## for model comparison
{ // local section, this saves time and space
for (i in 1:ns) {
log_lik[i] = 0;
for (t in 1:nt) {
log_lik[i] += ddm_(rt[i,t],choice[i,t],drift[i,t],A[i],z_t[i,t],t0[i]);
} // end of i loop
}
}
}
here is the log-likelihood returned by stan:
And the log-likelihood computed with rtdists:
