Trouble fitting a Q-Diffusion IRT model in Stan

Modeling
1

Dear stan forum users,

I ran into some trouble when fitting an IRT-Diffusion model via Stan. When I run the following code it first displays several messages of type:
“Rejecting initial value: Log probability evaluates to log(0), i.e. negative infinity. Stan can’t start sampling from this initial value.”
and then shows the error message "Initialization between (-2, 2) failed after 100 attempts. ". (I do not understand how an initialization between -2 and 2 appears considering the usage of invers-gamma priors for all model parameters - see model code below)

Therefore i fitted the same model but now with completely specified initial values (the true model parameter values used for generation). But then i get another error: It displays “Gradient evaluation took 0 seconds.
1000 transitions using 10 leapfrog steps would take 0 seconds…”
and then the program shuts down.

I would be grateful for any advices on how to fit the model properly.

Thank you for any help!

stan model

data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=0,upper=1> dec[N];  // decisions
  int<lower=1> k;  // number of items
  int<lower=1> J_1[N]; // vector which indicates subject id
  int<lower=1> I_1[N]; // vector which indicates item id
  int<lower=1> n; //number of subjects
  int<lower=1> p; //dimension of person parameter
}

transformed data {
  real min_Y = min(Y);
}

parameters {
  vector<lower=0>[n] theta;  // accuracy parameters; 1 per subject
  vector<lower=0>[n] g;  // speed_boundary parameter; 1 per subject
  vector<lower=0>[k] v;  // item (invers)drift parameter
  vector<lower=0>[k] a;  // item (invers)boundary parameter
  vector<lower=0, upper=min_Y>[k] ndt; // item non decision time
}

model {
  for (l in 1:k) {
    ndt[l] ~ inv_gamma(4, 3)T[0,min_Y];
    v[l] ~ inv_gamma(4, 3);
    a[l] ~ inv_gamma(4, 3);
  }
  for (l in 1:n) {
    theta[l] ~ inv_gamma(4, 3);
    g[l] ~ inv_gamma(4, 3);
  }
 
  for (j in 1:N) {
    if(dec[j]==1) {
     Y[j] ~ wiener(g[J_1[j]]/a[I_1[j]],ndt[I_1[j]],0.5,theta[J_1[j]]/v[I_1[j]]);
    }
    else {
     Y[j] ~ wiener(g[J_1[j]]/a[I_1[j]],ndt[I_1[j]],0.5,-theta[J_1[j]]/v[I_1[j]]);
   }
  }
}
2

The -2 to 2 is on the untransformed space (Stan samples on unconstrained reals and then does a transform to get the lower/upper bound constraints you specify).

The error, as you probably surmised, is that, for some reason, Stan wasn’t able to evaluate the log density or compute the gradients of the log density of the model.

You can get a code block by surrounding things with 3 backticks on each side (this one: `)

I’d assume the issue is that a constraint of the Wiener distribution is getting violated in here:

for (j in 1:N) {
  if(dec[j]==1) {
    Y[j] ~ wiener(g[J_1[j]]/a[I_1[j]],ndt[I_1[j]],0.5,theta[J_1[j]]/v[I_1[j]]);
  }  else{
    Y[j] ~ wiener(g[J_1[j]]/a[I_1[j]],ndt[I_1[j]],0.5,-theta[J_1[j]]/v[I_1[j]]);
  }
}

To confirm this, just comment out that bit of code and see if the model runs without it.

You can use prints to figure out exactly what went wrong:

real tau = ndt[I_1[j]];

if(Y[j] < 0 || Y[j] > tau) {
  print("Y[j] out of bounds!");
}

// etc...

Y[j] ~ wiener(g[J_1[j]]/a[I_1[j]], tau,0.5,theta[J_1[j]]/v[I_1[j]]);

I thought if there were constraint problems like that a more informative error should pop out… Ofc. sometimes these get swallowed up. Try commenting stuff out first.

Another thing you might try removing is the truncation here:
ndt[l] ~ inv_gamma(4, 3)T[0,min_Y];

You don’t need it since the parameters of the gamma are constant anyway (the normalizing term will be a constant), and I know the gamma stuff can be tricky numerically.

3

Thank you very much - I was able to fix it with your help!

4

Out of curiosity, what was the problem?