Why `PEs` produces NaN?

Modeling
1 
data {
  int<lower=1> N; // participants
  int<lower=1> T; // trials (=180)
  int<lower=1, upper=T> Tsubj[N];
  int<lower=1, upper=15> rounds[N,T];  // 15 rounds per PGG, total 12 PGG
  int<lower=0, upper=1> decision[N,T]; // decision for contribution
  int<lower=0, upper=1> result[N,T];   // result
  int<lower=2, upper=4> threshold[N,T];// k=2 or k=4
  real<lower=2, upper=4> thres[N,T];   // k=2 or k=4
  real<lower=0,upper=1> ratioC[N,T];   // ratio of contributors
  real<lower=0,upper=1> prior2[N,T];     // prior belief of #cont in k=2
  real<lower=0,upper=1> prior4[N,T];     // prior belief of #cont in k=4
}

parameters {
  vector[4] mu_p;
  real mu_l;
  vector<lower=0>[4] sigma;
  //Matt-trick
  vector[N] w_pr;
  vector[N] pi_pr;
  vector[N] lambda_pr;
  vector[N] alpha_pr;
  //vector[N] theta_pr;
}

transformed parameters{
  vector<lower=0,upper=1>[N] w;
  vector<lower=-1,upper=1>[N] pi;
  vector<lower=-10,upper=0>[N] lambda;
  vector<lower=0,upper=1>[N] alpha;
  //vector<lower=0,upper=1>[N] theta;
  
  for (i in 1:N) {
    w[i]      = Phi_approx( mu_p[1] + sigma[1] * w_pr[i] );
    pi[i]     = Phi_approx( mu_p[2] + sigma[2] * pi_pr[i] ) * 2 - 1;
    alpha[i]  = Phi_approx( mu_p[3] + sigma[3] * alpha_pr[i] ); 
    lambda[i] = Phi_approx( mu_l    + sigma[4] * lambda_pr[i] ) * (-10);
  }
}

model {
  //hyperparameters
  mu_p  ~ normal(0, 1.0); 
  mu_l  ~ normal(-5, 1.0);
  sigma ~ normal(0, 1.0);
  //individual parameters
  w_pr ~ normal(0, 1.0);
  pi_pr ~ cauchy(0, 1.0);
  lambda_pr ~ cauchy(0, 1.0);
  alpha_pr ~ normal(0, 1.0);
  //theta_pr ~ normal(0, 1.0);
  
  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real pCont;  //overall probability for contribution
      real I; real G; real Q;  //individual, group, total utility
      real gamma1; //belief of other's choice
      real PEs;                //social prediction error
      real K; real eG;
      real eI;
      int free;
      free = 5-threshold[i, t];   //desired limit of free-riders   
      
      //initial belief of other's choice
      if (rounds[i, t] == 1){
        if (threshold[i, t] == 2) gamma1 = 1 - prior2[i, t];
        else                      gamma1 = 1 - prior4[i, t];
      }
      //update belief via social reinforcement learning
      else{
        PEs = (1 - ratioC[i, t-1]) - gamma1; //social prediction error
        gamma1 = gamma1 + alpha[i] * PEs;    //update belief by PEs
        print("alpha PEs: ", alpha[i], PEs);
      } 
      if (gamma1 > 1-1e-3) gamma1 = 1-1e-3;
      if (gamma1 < 1e-3)   gamma1 = 1e-3;
      
      //compute Individual Utility
      eI = choose(free, 4) * gamma1^(free) * (1-gamma1)^(5-free);
      I = lambda[i] + eI * 2 + pi[i] * eI * 2 * 4;
        
      //compute Group Utility
      K = thres[i,t] / 5;
      eG = binomial_cdf(free, 4, gamma1);
      G = ((1-K^(16-rounds[i,t]))/(1-K)) * 2 * eG;
        
      //compute Total Utility
      Q = w[i] * I + (1-w[i]) * G;
      pCont = inv_logit(Q);
      decision[i,t] ~ bernoulli(pCont);
    }
  }
}

From this stan code, I debugged with ‘print("alpha PEs: ", alpha[i], PEs);’
And PEs = (1 - ratioC[i, t-1]) - gamma1; yielded NaN.

ratioC[i, t-1] comes from the data, ranging from 0 to 1,
gamma1 was initiated by
if (rounds[i, t] == 1){
** if (threshold[i, t] == 2) gamma1 = 1 - prior2[i, t];**
** else gamma1 = 1 - prior4[i, t];**
** }**
When I printed ratioC[i, t-1], gamma1, alpha[i] they had appropriate value like
Chain 1: rat: 0.5
alpha PEs: 0.335983nan
What is the matter…? Frustrating

2

I think this is about programming, not really Stan. gamma1 is not defined before printing alphaPEs

1 Like
3

Oh I had to initiate gamma1 outside the 2nd for loop!
Thanks!

1 Like