Reproducing "Circular Inference" model in Stan - Divergent iterations issue

[Riccardo Fusaroli]

I am currently trying to reproduce the circular inference model sketched in this paper: https://www.nature.com/articles/ncomms14218 .

When I try to fit the model I get lots of divergent iterations (and Rhat goes wild) and I could use some help in figuring out whether I'm making a trivial mistake somewhere, or as Michael Betancourt mentioned on Twitter "Your model is multimodal or has an awkward geometry and needs reparameterization". In that case, i could use some help in figuring out the best way of tacking that.

Sample Data (5 out of ca 80 participants): https://www.dropbox.com/s/g6jcj85a3w6gscg/d2.csv?dl=0

The task is a gumball task, where participants have to guess whether the next gumball extracted is red (or green). The dataset includes: 

• Participant Number, 
• ChoiceRed: whether the participant indicated the next gumball is red or not
  
• Self: an index of how many red gumballs the participant has seen in the urn (expressed on a scale from 0 to 1)
• lSelf: the logit of Self
• Others: an index of how many other people the participant sees guessing red (expressed on a scale from 0 to 1)
• lOthers: the logit of Others

First as a sanity check I implemented a weighted bayesian model (using the rethinking package as a crutch), which gives a few divergences, but nice traceplots and Rhats=1 on both simulated and real data. The model follows the paper previously linked, only adding the possibility for biases towards red (no matter the info) with a[Participant] and pooling information across participants.

data{

    int<lower=1> N;

    int<lower=1> N_Participant;

    int ChoiceRed[N];

    int Participant[N];

    real lSelf[N];

    real lOthers[N];

}

parameters{

    vector[N_Participant] a;

    real muA;

    vector[N_Participant] wOthers;

    real muO;

    vector[N_Participant] wSelf;

    real muS;

}

model{

    vector[N] mu;

    muS ~ normal( 0.5 , 0.1 );

    wSelf ~ normal( muS , 0.3 );

    muO ~ normal( 0.5 , 0.1 );

    wOthers ~ normal( muO , 0.3 );

    muA ~ normal( 0.5 , 0.1 );

    a ~ normal( muA , 0.1 );

    for ( i in 1:N ) {

        mu[i] = a[Participant[i]] + log((wSelf[Participant[i]] * exp(lSelf[i]) + 1 - wSelf[Participant[i]])/((1 -      wSelf[Participant[i]]) * exp(lSelf[i]) + wSelf[Participant[i]])) +      log((wOthers[Participant[i]] * exp(lOthers[i]) + 1 - wOthers[Participant[i]])/((1 -          wOthers[Participant[i]]) * exp(lOthers[i]) + wOthers[Participant[i]]));

    }

    ChoiceRed ~ binomial_logit( 1 , mu );

}

generated quantities{

    vector[N] mu;

    real dev;

    dev = 0;

    for ( i in 1:N ) {

        mu[i] = a[Participant[i]] + log((wSelf[Participant[i]] * exp(lSelf[i]) + 1 - wSelf[Participant[i]])/((1 -      wSelf[Participant[i]]) * exp(lSelf[i]) + wSelf[Participant[i]])) +      log((wOthers[Participant[i]] * exp(lOthers[i]) + 1 - wOthers[Participant[i]])/((1 -          wOthers[Participant[i]]) * exp(lOthers[i]) + wOthers[Participant[i]]));

    }

    dev = dev + (-2)*binomial_logit_lpmf( ChoiceRed | 1 , mu );

}

Next I tried to implement the circular inference model. Sorry for the ugliness of the code, I'm still figuring out how to pass functions to the model definition.

data{

    int<lower=1> N;

    int<lower=1> N_Participant;

    int ChoiceRed[N];

    int Participant[N];

    real lSelf[N];

    real lOthers[N];

}

parameters{

    vector[N_Participant] a;

    vector[N_Participant] aSelf;

    vector[N_Participant] aOthers;

    vector[N_Participant] wSelf;

    vector[N_Participant] wOthers;

    real muA;

    real muS;

    real muO;

    real<lower=0> scaleS;

    real<lower=0> scaleO;

}

model{

    vector[N] mu;

    scaleO ~ normal( 0 , 3 );

    scaleS ~ normal ( 0 , 3 );

    muO ~ normal( 0 , 1 );

    muS ~ normal( 0 , 1 );

    muA ~ normal( 0 , 1 );

    wOthers ~ normal( muO , 3 );

    wSelf ~ normal( muS , 3 );

    for ( j in 1:N_Participant ) aOthers[j] ~ cauchy( 0 , scaleO );

    for ( j in 1:N_Participant ) aSelf[j] ~ cauchy( 0 , scaleS );

    a ~ normal( muA , 3 );

    for ( i in 1:N ) {

        mu[i] = a[Participant[i]]+

    #F(lSelf + fOthers + fSelf, wSelf) +

    log((wSelf[Participant[i]]*exp(lSelf + 

                                  log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]]))  + 

                                  log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]])) )+

           1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(lSelf + 

                                                               log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]]))  + 

                                                               log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]])) )+

                                    wSelf[Participant[i]])) +

    #F(lOthers + fSelf + fOthers,wOthers)

    log((wOthers[Participant[i]]*exp(lOthers + 

                                    log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]]))  + 

                                    log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]])) )+

           1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(lOthers + 

                                                                   log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]]))  + 

                                                                   log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]])) )+wOthers[Participant[i]]))

    }

    ChoiceRed ~ binomial_logit( 1 , mu );

}

generated quantities{

    vector[N] mu;

    real dev;

    dev = 0;

    for ( i in 1:N ) {

        mu[i] = a[Participant[i]]+

    #F(lSelf + fOthers + fSelf, wSelf) +

    log((wSelf[Participant[i]]*exp(lSelf + 

                                  log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]]))  + 

                                  log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]])) )+

           1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(lSelf + 

                                                               log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]]))  + 

                                                               log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]])) )+

                                    wSelf[Participant[i]])) +

    #F(lOthers + fSelf + fOthers,wOthers)

    log((wOthers[Participant[i]]*exp(lOthers + 

                                    log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]]))  + 

                                    log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]])) )+

           1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(lOthers + 

                                                                   log((wSelf[Participant[i]]*exp(aSelf[Participant[i]]*lSelf)+1-wSelf[Participant[i]])/((1-wSelf[Participant[i]])*exp(aSelf[Participant[i]]*lSelf)+wSelf[Participant[i]]))  + 

                                                                   log((wOthers[Participant[i]]*exp(aOthers[Participant[i]]*lOthers)+1-wOthers[Participant[i]])/((1-wOthers[Participant[i]])*exp(aOthers[Participant[i]]*lOthers)+wOthers[Participant[i]])) )+wOthers[Participant[i]]))

    }

    dev = dev + (-2)*binomial_logit_lpmf( ChoiceRed | 1 , mu );

}

Many thanks for the attention.

[Ben B]

I don't know about the model you're showing, but I think your idea of using functions to clean up your code is good. Maybe it's just a bug that's making the inference hard. I removed the generated quantities section for brevity (but you should go ahead and add it back). I only tested that this compiles. I did *not* actually do a sampling, so double check me, but rewriting your harder problem with Fs might make it easier to see what's going on.

functions {

  real F(real L, real w) {

    real eL = exp(L);

    return log((w * eL + 1 - w)/((1 - w) * eL + w));

  }

}

data{

  int<lower=1> N;

  int<lower=1> N_Participant;

  int ChoiceRed[N];

  int Participant[N];

  real lSelf[N];

  real lOthers[N];

}

parameters{

  vector[N_Participant] a;

  real muA;

  vector[N_Participant] wOthers;

  real muO;

  vector[N_Participant] wSelf;

  real muS;

}

model{

  vector[N] mu;

  muS ~ normal( 0.5 , 0.1 );

  wSelf ~ normal( muS , 0.3 );

  muO ~ normal( 0.5 , 0.1 );

  wOthers ~ normal( muO , 0.3 );

  muA ~ normal( 0.5 , 0.1 );

  a ~ normal( muA , 0.1 );

  for ( i in 1:N ) {

    mu[i] = a[Participant[i]] + F(lSelf[i], wSelf[Participant[i]]) + F(lOthers[i], wOthers[Participant[i]]);

  }

  ChoiceRed ~ binomial_logit( 1 , mu );

}

Ben

[Bob Carpenter]

You probably want a non-centered parameterization for a
and wOthers.

> for ( i in 1:N ) {
> mu[i] = a[Participant[i]] + F(lSelf[i], wSelf[Participant[i]]) + F(lOthers[i], wOthers[Participant[i]]);
> }


>
> ChoiceRed ~ binomial_logit( 1 , mu );

This can be simplified to:

ChoiceRed ~ bernoulli_logit(mu);

- Bob

[Riccardo Fusaroli]

Thanks for the suggestions. I have now rewritten the models, but I still have issues in getting them to work.

When I run the weighted bayes model I get a bunch of "Logit transformed probability parameter is nan" and the model has a hard time reaching Rhat's of 1 (though it stays close, 1.01s and 1.02s).

When I run the circular inference model I get a "Error in stanc(file = file, model_code = model_code, model_name = model_name,  : "

Here is the data: https://www.dropbox.com/s/g6jcj85a3w6gscg/d2.csv?dl=0

And here the models

Weighted Bayes

functions {

  real F(real L, real w) {

    real eL = exp(L);

    return log((w * eL + 1 - w)/((1 - w) * eL + w));

  }

}

data{

    int<lower=1> N;

    int<lower=0, upper=1> ChoiceRed[N];

    real lSelf[N];

    real lOthers[N];

    int<lower=1> P;

    int<lower=1, upper=P> Participant[N];

}

parameters{

    vector[P] a;

    vector[P] wOthers;

    vector[P] wSelf;

}

model{

    vector[N] mu;

    wSelf ~ normal( 0.5 , 0.1 );

    wOthers ~ normal( 0.5 , 0.1);

    a ~ normal( 0.5 , 0.1 );

  for ( i in 1:N ) {

    mu[i] = a[Participant[i]] + F(lSelf[i], wSelf[Participant[i]]) + F(lOthers[i], wOthers[Participant[i]]);

  }

  ChoiceRed ~ bernoulli_logit(mu);

}

generated quantities{

    vector[N] mu;

    real dev;

    dev = 0;

    for ( i in 1:N ) {

    mu[i] = a[Participant[i]] + F(lSelf[i], wSelf[Participant[i]]) + F(lOthers[i], wOthers[Participant[i]]);

  }

    dev = dev + (-2)*bernoulli_logit_lpmf( ChoiceRed | mu );

}

Circular inference

functions {

  real F(real L, real w) {

    real eL = exp(L);

    return log((w * eL + 1 - w)/((1 - w) * eL + w));

  }

}

data{

  int<lower=1> N;

  int<lower=0, upper=1> ChoiceRed[N];

  real lSelf[N];

  real lOthers[N];

  int<lower=1> P;

  int<lower=1, upper=P> Participant[N];

}

parameters{

  vector[P] a;

  vector[P] wOthers;

  vector[P] wSelf;

}

model{

  vector[N] mu;

  wSelf ~ normal( 0.5 , 0.3 );

  wOthers ~ normal( 0.5 , 0.3 );

  a ~ normal( 0.5, 0.1 );

  for ( i in 1:N ) {

    mu[i] = a[Participant[i]] + F(lSelf[i] + F(aOthers[Participant[i]]*lOthers[i],wOthers[Participant[i]]) + F(aSelf[Participant[i]]*lSelf[i],wSelf[Participant[i]]), wSelf[Participant[i]]) + F(lOthers[i] + F(aSelf[Participant[i]]*lSelf[i],wSelf[Participant[i]]) + F(aOthers[Participant[i]]*lOthers[i],wOthers[Participant[i]]), wOthers[Participant[i]])

   }

  ChoiceRed ~ bernoulli_logit(mu);

}

generated quantities{

    vector[N] mu;

    real dev;

    dev = 0;

    for ( i in 1:N ) {

    mu[i] = a[Participant[i]] + F(lSelf[i] + F(aOthers[Participant[i]]*lOthers[i],wOthers[Participant[i]]) + F(aSelf[Participant[i]]*lSelf[i],wSelf[Participant[i]]), wSelf[Participant[i]]) + F(lOthers[i] + F(aSelf[Participant[i]]*lSelf[i],wSelf[Participant[i]]) + F(aOthers[Participant[i]]*lOthers[i],wOthers[Participant[i]]), wOthers[Participant[i]])

   }

    dev = dev + (-2)*bernoulli_logit_lpmf( ChoiceRed | mu );

}

Thanks again for the help

[Ben B]

I think the "Logit transformed probability parameter is nan" errors make sense. That would happen whenever (w * eL + 1 - w)/((1 - w) * eL + w) < 0 (so you can't take the log of it).

Do you think the results of the weighted Bayes model is correct?

The error I got when I tried to compile the second model was "aOthers" doesn't exist. Which seems true. Are these missing parameters?

Ben

[Riccardo Fusaroli]

Thanks again for the response. 

Weighted bayes looks correct. But the problem with the circular inference is still there. The aOthers slipped out as I was copying the model here, sorry. I attach the model file, so to avoid that issue again.

When I run it I get:

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  c++ exception (unknown reason)

[Ben B]

It looks like you're missing semi-colons at the end of the two long lines (mu[i] = ...). If you add those, the model should compile (I'm assuming this is what the error was). I didn't try to run it.

Ben

[Riccardo Fusaroli]

Oh darn, sorry for the rookie mistake, I really appreciate you putting up with them. That's corrected. Now I have a clean working model, but I'm back to my initial problem: 

- the model has a really hard time estimating the different parameters, even when using a humongous amount of iterations. See the very low n_eff, high Rhat in the attached image (and there's lot of brown noise in the traceplots) .

Here the code I used to run the model (and a link to the data: https://www.dropbox.com/s/i4um11kfpt8ycem/d1.csv?dl=0):

data <- list(

  Participant = as.integer(d1$Participant),

  ChoiceRed = d1$ChoiceRed,

  lSelf= d1$lSelf,

  lOthers= d1$lOthers,

  N=nrow(d1),

  P = length(unique(d1$Participant))

)

fit3 <- stan(

  file = "CircularInference.stan",  

  data = data,    

  chains = 2,             

  warmup = 3e3,          

  iter = 1e4,            

  cores = 2,               

  refresh = 1000,          

  init=0

)

[Riccardo Fusaroli]

small addendum: if I run the model on a subset of participants (2-5), the model seems to have fewer problems in converging, with n_eff being decent and Rhat <1.1. 

[Ben B]

I don't know about the divergent iterations thing, but if you have a low neff and a high Rhat, it'd be a good idea to go look at the posterior distributions/traceplots.

You might look at the posterior distributions of each chain to make sure everything is converging to the same solution. Multimodal distributions could give these the Rhat/neff metrics trouble I think.

You might also look at pairplots for your parameters too see if things are really tightly correlated. Maybe there are identifiability issues that are giving the sampler problems (that could be addressed with some priors).

Do you have simulated data for this problem? Might also scale around the size of the data and see if anything changes if so.