Why my sampling not done

Hi all,

I’m trying to use bayesian hierarchial modeling method to fit my data into a simple Expected Utility model. However I’m stuck with keeping getting an error saying that “error occurred during calling the sampler: sampling not done”. Is there any way to track what’s wrong with my code (see below)? Any help is much apprecated! Thank you!


data {
  int N; // number of trials total (across participants)
  int nsubj;  // number of subjects
  int subjs[N];
  int choices[N]; // the binary choice vector
  real x1a[N];
  real x2a[N];
  real x3a[N];
  real x4a[N];
  real x1b[N];
  real x2b[N];
  real x3b[N];
  real x4b[N];
  real p1a[N];
  real p2a[N];
  real p3a[N];
  real p4a[N];
  real p1b[N];
  real p2b[N];
  real p3b[N];
  real p4b[N];
  int  sid[N]; // the subject index
}

parameters {
  // parameters for group-level distributions of individuals' baseline parameters
  real meanalpha; // alpha in EU
  real<lower=0> sdalpha;
  real meannoise; //noise
  real<lower=0> sdnoise;
  
  // individual-level parameters
  real alpha[nsubj];
  real noise[nsubj];
}

model {
  real p[N]; // trialwise likelihood
  real ua; // utility of the option a
  real ub; // utility of the option b
  
  // Group-level distribution paramter priors
  meanalpha ~ normal(0,100);
  sdalpha ~ inv_gamma(0.001,0.001);
  meannoise ~ normal(0,100);
  sdnoise ~ inv_gamma(0.001,0.001);
  

  
  // Hierarchy
  alpha ~ normal(meanalpha, sdalpha);
  noise ~ normal(meannoise, sdnoise);


  for (i in 1:N) {
    for (t in 1:subjs[i]) {
       // calcualte the utility of the gamble and the alternative using alpha and m
      // then use the softmax function to estiamte the likelihood of gambling in each trial
      ua=p1a[t]*x1a[t]^alpha[t]+p2a[t]*x2a[t]^alpha[t]+p3a[t]*x3a[t]^alpha[t]+p4a[t]*x4a[t]^alpha[t];
      ub=p1b[t]*x1b[t]^alpha[t]+p2b[t]*x2b[t]^alpha[t]+p3b[t]*x3b[t]^alpha[t]+p4b[t]*x4b[t]^alpha[t];
      if(choices[t]==1)
      p[t]=inv_logit((ua-ub)/noise[t]);
      else
      p[t]=inv_logit((ub-ua)/noise[t]);
      choices ~ bernoulli(p);
    }
  }
}

When I execute my code, I see message below:

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: bernoulli_lpmf: Probability parameter[1] is nan, but must be finite! (in 'modelc468ceada7_EU' at line 65)

Chain 1:
Chain 1: Initialization between(-1,1) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]): Initialization failed."
[1]"error occurred during calling the sampler: sampling not done"

The error points to line 65

      choices ~ bernoulli(p);

My first observation is that this line is inside a loop so it shouldn’t be vectorized. You probably meant

      choices[i] ~ bernoulli(p[i]);

(Actually, the error says parameter[1] is nan so there’s a problem calculating p regardless of vectorization.)
The lines right before are also a bit weird

      if(choices[t]==1)
      p[t]=inv_logit((ua-ub)/noise[t]);
      else
      p[t]=inv_logit((ub-ua)/noise[t]);

These are complementary probabilities chosen depending on the value of choices[t]. That’s basically a Bernoulli distribution. But there’s already a Bernoulli distribution on the next line and these two end up sort of canceling each other out. Maybe the model you wanted is

      choices[i] ~ bernoulli_logit((ub-ua)/noise[t]);

(note the use of bernoulli_logit(..) as a shorthand for bernoulli(inv_logit(..)))
…except that’s not quite right because of p[i] vs p[t] indexing confusion. I don’t quite understand what the inner loop is supposed to do. You calculate several probabilities for each trial but there’s only one choice per trial so something does not add up. The model also does not make any use of int sid[N]; // the subject index.
Assuming each trial has only one subject (numbered by sid) the model that computes the utility for that subject only and predicts the subject’s choice is

  for (i in 1:N) {
    // calcualte the utility of the gamble and the alternative using alpha and m
    // then use the softmax function to estiamte the likelihood of gambling in each trial
    int t = sid[i];
    ua=p1a[i]*x1a[i]^alpha[t]+p2a[i]*x2a[t]^alpha[t]+p3a[i]*x3a[i]^alpha[t]+p4a[i]*x4a[i]^alpha[t];
    ub=p1b[i]*x1b[i]^alpha[t]+p2b[i]*x2b[t]^alpha[t]+p3b[i]*x3b[i]^alpha[t]+p4b[i]*x4b[i]^alpha[t];
    choices[i] ~ bernoulli_logit((ua-ub)/noise[t]);
  }

The exponentiation x1a[i]^alpha[t] only works if the x1a is positve so make sure there are no zeros or negative values in the data.

  real<lower=0> x1a[N];
  real<lower=0> x2a[N];
  real<lower=0> x3a[N];
  real<lower=0> x4a[N];
  real<lower=0> x1b[N];
  real<lower=0> x2b[N];
  real<lower=0> x3b[N];
  real<lower=0> x4b[N];
2 Likes

Dear Niko,

Thank you very much for your reply. As you pointed out, I think there are big problems in the for loops. You are right that there’s only one choice per trial. But each subject has multiple trials. What I tend to achieve is the model that computes the utility of each trial for that subject and predicts the subject’s choice for that trial. Would you mind show me how to achieve this?

In addition, in your last note, you mention that the exponentiation xia[i]^alpha[t] only works if the x1a is positive. However, in my data, there are indeed some cases where x1a takes the value of zero. This is how the trials were created. Is there any way to get around this? Should I leave those x1a as blank in the dataset instead of “0”?

Thank you very much!

The subject in the ith trial is subject number sid[i], right? Then the above code should do that.

Well, the question is, what is the subject to do when one of xs is zero? The model allows a subject to have a negative alpha but if a subect with alpha=-1 sees x1a=0.0 their utility is 0^{-1} which is not going to give a meaningful prediction. I assume you want the subject to behave as if the utility was zero (i.e. ignore that predictor). Then one way to write the model is

for (i in 1:N) {
    int t = sid[i];
    ua = 0.0;
    if (x1a[i] != 0.0)
        ua += p1a[i]*x1a[i]^alpha[t];
    if (x2a[i] != 0.0)
        ua += p2a[i]*x2a[i]^alpha[t];
    if (x3a[i] != 0.0)
        ua += p3a[i]*x3a[i]^alpha[t];
    if (x4a[i] != 0.0)
        ua += p4a[i]*x4a[i]^alpha[t];
    ub = 0.0;
    if (x1b[i] != 0.0)
        ub += p1b[i]*x1b[i]^alpha[t];
    if (x2b[i] != 0.0)
        ub += p2b[i]*x2b[i]^alpha[t];
    if (x3b[i] != 0.0)
        ub += p3b[i]*x3b[i]^alpha[t];
    if (x4b[i] != 0.0)
        ub += p4b[i]*x4b[i]^alpha[t];
    choices[i] ~ bernoulli_logit((ua-ub)/noise[t]);
}
1 Like

This worked perfectly! Thank you very much!