I am trying to fit a model to a couple (19) datasets. For 16 datasets, the model runs smoothly (no divergences, reasonable Rhat, and very good ESSs). But it seems, three datasets are a bit hard to estimate. I determined the best possible initial values to prevent infinite log-likelihoods, however, I cannot get rid of the “Gradient evaluated at the initial value is not finite.” error.
This is my model (which is cumulative prospect theory):
functions {
real prob_weighting ( real p, real gamma) {
return exp(-(-log(p))^gamma);
}
}
data {
////////////////////////
// Context variables:
int<lower=0> N;
int<lower=0> Ntrials;
int<lower=0> S1; // number of states of first option
int<lower=0> S2; // number of states of second option
int<lower=1, upper=N> sbj[Ntrials]; // subject identifier
vector[4] prior_mu;
vector[4] prior_sd;
////////////////////////
// Predictor variables:
// IMPORTANT: Note, that outcomes and corresponding probabilities are expected
// to be sorted according to the outcomes within each option and trial!
// Positive outcomes should come first and be sorted decreasingly, losses
// should come second and be sorted increasingly!
// In addition, probabilities should already be included as CUMULATIVE Probabilities!
vector[S1] outcomes1[Ntrials]; //values of the outcomes of the first option in each trial
vector[S1] cumprobs1[Ntrials]; // cumulative outcome probabilities of first option in each trial
vector[S2] outcomes2[Ntrials]; //values of the outcomes of the second option in each trial
vector[S2] cumprobs2[Ntrials]; // cumulative outcome probabilities of second option in each trial
////////////////////////
// Outcome variables:
int<lower=0, upper=1> choice[Ntrials]; // Final choice: 1=chooseX; 2 = chooseY
}
parameters {
// population means
real mu_log_dscale;
real mu_log_alpha;
real mu_log_gamma;
real mu_log_lambda;
// population SDs
real<lower=0> sd_log_dscale;
real<lower=0> sd_log_alpha;
real<lower=0> sd_log_gamma;
real<lower=0> sd_log_lambda;
// non-centered subject effects
vector[N] z_log_dscale;
vector[N] z_log_alpha;
vector[N] z_log_gamma;
vector[N] z_log_lambda;
}
transformed parameters {
// subject parameters
vector<lower=0>[N] dscale;
vector<lower=0>[N] alphas;
vector<lower=0>[N] gammas;
vector<lower=0>[N] lambdas;
dscale = log1p_exp(mu_log_dscale + sd_log_dscale * z_log_dscale);
alphas = log1p_exp(mu_log_alpha + sd_log_alpha * z_log_alpha);
gammas = log1p_exp(mu_log_gamma + sd_log_gamma * z_log_gamma);
lambdas = log1p_exp(mu_log_lambda + sd_log_lambda * z_log_lambda);
}
model {
// population means
mu_log_dscale ~ normal(prior_mu[1] , prior_sd[1]);
mu_log_alpha ~ normal(prior_mu[2] , prior_sd[2]);
mu_log_gamma ~ normal(prior_mu[3] , prior_sd[3]);
mu_log_lambda ~ normal(prior_mu[4] , prior_sd[4]);
// population SDs
sd_log_dscale ~ normal(0, 1);
sd_log_alpha ~ normal(0, 1);
sd_log_gamma ~ normal(0, 1);
sd_log_lambda ~ normal(0, 1);
// non-centered latent effects
z_log_dscale ~ std_normal();
z_log_alpha ~ std_normal();
z_log_gamma ~ std_normal();
z_log_lambda ~ std_normal();
// Likelihood
real drift;
int curM;
vector[S1] cur_outcomes1;
vector[S2] cur_outcomes2;
vector[S1] weights1;
vector[S2] weights2;
real prev_weighted_cumprob, weighted_cumprob;
int switchsign;
for (j in 1:Ntrials){
cur_outcomes1 = outcomes1[j];
cur_outcomes2 = outcomes2[j];
switchsign = 0;
prev_weighted_cumprob = 0.0; // reset prop_weight
for (out in 1:S1) {
if (cur_outcomes1[out] > 0) {
weighted_cumprob = prob_weighting(cumprobs1[j][out], gammas[sbj[j]]);
weights1[out] = weighted_cumprob - prev_weighted_cumprob;
prev_weighted_cumprob = weighted_cumprob;
cur_outcomes1[out] = cur_outcomes1[out]^alphas[sbj[j]];
} else {
if (switchsign==0) {
prev_weighted_cumprob = 0.0;
switchsign=1;
}
weighted_cumprob = prob_weighting(cumprobs1[j][out], gammas[sbj[j]]);
weights1[out] = weighted_cumprob - prev_weighted_cumprob;
prev_weighted_cumprob = weighted_cumprob;
cur_outcomes1[out] = - lambdas[sbj[j]] * (-cur_outcomes1[out])^alphas[sbj[j]];
}
}
switchsign = 0;
prev_weighted_cumprob = 0.0;
for (out in 1:S2) {
if (cur_outcomes2[out] > 0) {
weighted_cumprob = prob_weighting(cumprobs2[j][out], gammas[sbj[j]]);
weights2[out] = weighted_cumprob - prev_weighted_cumprob;
prev_weighted_cumprob = weighted_cumprob;
cur_outcomes2[out] = cur_outcomes2[out]^alphas[sbj[j]];
} else {
if (switchsign==0) {
prev_weighted_cumprob = 0.0;
switchsign=1;
}
weighted_cumprob = prob_weighting(cumprobs2[j][out], gammas[sbj[j]]);
weights2[out] = weighted_cumprob - prev_weighted_cumprob;
prev_weighted_cumprob = weighted_cumprob;
cur_outcomes2[out] = - lambdas[sbj[j]] * (-cur_outcomes2[out])^alphas[sbj[j]];
}
}
drift = dscale[sbj[j]] *
(sum( cur_outcomes1 .* weights1)-
sum( cur_outcomes2 .* weights2));
// Choice probability
choice[j] ~ bernoulli_logit(drift);
}
}
There you see that I even started tuning my priors to be suitable for the data. With this, I could actually fit two of the three problematic datasets, but the last one still resists all my attempts.
It contains about 270 choices from 64 participants each (in total 17320 observations).
I used the mean of participant-level maximum likelihood estimates for the parameters for the initialization of the group-level priors for the means (mu_…) their standard deviation for the sd_…'s, and the scaled individual offsets for the random effects (z_…), so the initials should be fine actually.
Also, pooling all participants by ignoring the hierarchical structure and investigating the likelihood surface for the population means produces reasonably smooth curves.
I attach the problematic sub-data and a minimal R-script to reproduce the error.
Is there something problematic with my model or any strategy that I could use to find out what exactly is going wrong?
Fit_choice_data.R (9.0 KB)
choice_data.csv (780.6 KB)