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

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)