Failed initialisation and parameter recovery - Initialization between (-2, 2) failed after 100 attempts

Modeling
, ,
1

Hi, I am new to Stan and have a question about an error message I keep getting regarding the initialisation of the chains. I’ve read though the user manual of Stan, and know that there are a couple of posts discussing the matter already, but none of the suggestions seem to work for me.

I am trying to model choice behaviour using prospect theory (PT). I simulate choice data using PT, which I then feed into my Stan model to recover my parameters.

In summary, the simulated participants have to choose between too prospects, which in turn differ in outcome x_{i} and probability p_{i}. All outcomes are gains (positive values). I calculate the utility of an outcome using a power utility function: u(x_{i}) = x^{a}. I calculate the probability weights using the Prelec weighting function w(p_{i}) = e^{-delta*(-log(p_{i}))^{gamma}}. The value of each prospect is then calculated by multiplying weight and utility. Choice probability is calculated using an inverse logit function. I want to recover the parameters for each individual separately.

When running the model I keep getting following error message when using simulated data of more than 10 participants. Additionally I am unable to recover my parameters properly. I would be extremely grateful for all input. I think that the issue may also be that my model is not set up correctly.

Chain 1: Initialization between (-2, 2) 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."
error occurred during calling the sampler; sampling not done


// Model: PT
// Takes choice data and gamble attributes as input
data {
  int<lower=0> N;
  int<lower=0> Nsubs;
  vector[N] magA; // outcome prospect A
  vector[N] magB; // outcome prospect B
  vector[N] probA; // probability prospect A
  vector[N] probB; // probability prospect B
  int <lower=0,upper=1> choice[N, Nsubs]; // choice outcome, 0 or 1
}

parameters {
  
  real<lower=0, upper=2> alpha[Nsubs]; 
  real <lower=0, upper=2> delta[Nsubs];
  real<lower=0, upper=2>  gamma[Nsubs];
  real <lower=0, upper=10> eta[Nsubs];
}

model {
     
  real uA[N, Nsubs];
  real uB[N, Nsubs];
  real pA[N, Nsubs];
  real pB[N, Nsubs];
  real vA[N, Nsubs];
  real vB[N, Nsubs];
  real p[N,Nsubs] ;
   
  alpha[Nsubs] ~ uniform(0, 2);
  delta[Nsubs] ~ uniform(0, 2);
  gamma[Nsubs] ~ uniform(0, 2);
  eta[Nsubs]  ~ uniform(0, 5);
 
for (i in 1:Nsubs) {  
  for (j in 1:N) {  
    
 uA[j,i] = magA[j]^alpha[i];
 uB[j,i] = magB[j]^alpha[i];
 pA[j,i] = exp(-delta[i]*(-log(probA[j]))^gamma[i]);
 pB[j,i] = exp(-delta[i]*(-log(probB[j]))^gamma[i]);
 vA[j,i] = uA[j,i]*pA[j,i];
 vB[j,i] = uB[j,i]*pB[j,i];
 p[j,i] = exp(eta[i]*(vA[j,i]))/(exp(eta[i]*(vA[j,i]))+exp(eta[i]*(vB[j,i])));
}

 choice[N ,i] ~ bernoulli(p[N,i]);
}
}

For the simulation I run:

rm(list = ls())

#-------------------------------------------------------------------------------

# load required libraries
library(rstan)
library(permute)
library(ggplot2)
library(GGally)
#library(rstanarm)
#-------------------------------------------------------------------------------

# set up choice task

# experimental set-up -----------------------------------------------------
Nblock = 1
Nchoice = 128 #
Nsubs = 10
Ntotal = Nblock * Nsubs

choiceMat = read.csv("gambles_risk_pref.csv")#, header = TRUE)

magA = as.numeric(rep(choiceMat[c(2:33),3],4))
magB = as.numeric(rep(choiceMat[c(2:33),6],4))
probA = as.numeric(rep(choiceMat[c(2:33),2],4))
probB = as.numeric(rep(choiceMat[c(2:33),5],4))

# simulation --------------------------------------------------------------------
uA = matrix(0, Nchoice, Nsubs)
uB = matrix(0, Nchoice, Nsubs)
pA = matrix(0, Nchoice, Nsubs)
pB = matrix(0, Nchoice, Nsubs)
vA = matrix(0, Nchoice, Nsubs)
vB = matrix(0, Nchoice, Nsubs)
p = matrix(0, Nchoice, Nsubs)
choice = matrix(0, Nchoice, Nsubs)
choice1 = matrix(0, Nchoice, Nsubs)

alpha = numeric(Nsubs)
delta = numeric(Nsubs)
gamma = numeric(Nsubs)
eta = numeric(Nsubs)

eta = numeric(Nsubs)#runif(1, 0, 10)
#set.seed()
for (j in 1:Nsubs) {

  alpha[j] = runif(1, 0, 2) #create range of values instead eg seq(0,2,0.1)
  delta[j] = 1#runif(1, 0, 2)#create range of values
  gamma[j] = 1#runif(1, 0, 1)#create range of values
  eta[j] = runif(1, 0, 3)#create range of values

for (i in 1:Nchoice){

  uA[i,j] = magA[i]^alpha[j];
  uB[i,j] = magB[i]^alpha[j];
  pA[i,j] = exp(-delta[j]*(-log(probA[i]))^gamma[j]);
  pB[i,j] = exp(-delta[j]*(-log(probB[i]))^gamma[j]);
  vA[i,j] = uA[i,j]*pA[i,j]
  vB[i,j] = uB[i,j]*pB[i,j]
  p[i,j] = exp(eta[j]*(vA[i,j]))/(exp(eta[j]*(vA[i,j]))+exp(eta[j]*(vB[i,j])))
  pp <- runif(1, 0, 1 ) 
  choice[i,j] = as.numeric(p[i,j] >= pp) 
}
}


data <- list(N = Nchoice, Nsubs=Nsubs, magA = magA, magB = magB, probB = probB, probA = probA, choice=choice)
fit <- stan(file = 'cpt1.stan', data = data)
2

Thanks for posting all the code. Can you post a snippet of the data or some simulated data? You can specify inits in the model call but typically I’d dive into the model with some print statements to see what values are causing problems.