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)
```