Hello everyone,

I am new to modelling and Stan, and I am interested in implementing a simple Prospect Theory model to my data (initially to data from a single participant as a test, but I would like to make it hierarchical further down the line). Within the model, the components of an option (value * probability) are transformed via the following functions:

value\_corrected = value^{alpha}

probability\_corrected = probability^{gamma}/{(probability^{gamma}+(1-probability)^{gamma})^{(1/gamma)}}

And finally, the probability of accepting the risky option is given by

1/{1+exp(-mu*value\_difference)}

So, I implemented the following model in Stan:

```
data {
int<lower=0> N;//the number of trials for a single participant
int accept[N]; //a vector holding the sequence of responses
vector[N] RM; //a vector holding the risky magnitude for each trial
vector[N] SM; //a vector holding the safe magnitude for each trial
vector[N] RP; //a vector holding the risky probability for each trial
}
parameters {
real<lower=0.1,upper=5> alpha;
real<lower=0.1,upper=5> gamma;
real<lower=0,upper=2> mu;
}
model {
//priors
//alpha~normal(0.8,0.1);
//gamma~cauchy(0.8,1);
//mu~uniform(0.5,1.5); //reminder: no need to specify uniform prior here if the parameter is defined as bounded
//vectors for holding the transformed values/probs
vector[N] risky_value;
vector[N] safe_value;
vector[N] risky_prob;
vector[N] VD;
vector[N] VD_cor;
//transforming the raw values/probs
for (i in 1:N){
risky_value[i] <- pow(RM[i],alpha);
safe_value[i] <- pow(SM[i],alpha);
risky_prob[i] <- pow(RP[i], gamma)/pow((pow(RP[i],gamma)+pow(1-RP[i],gamma)),inv(gamma));
}
//computing the value difference and the probability of accepting the risky option
for (i in 1:N){
VD[i] <- risky_value[i]*risky_prob[i] - safe_value[i];
VD_cor[i] <- 1/1+exp(-mu * VD[i]);
}
for (i in 1:N){
accept[i] ~ bernoulli_logit (VD_cor[i]);
}
}
```

The model seems to converge as indicated by low RHat values, occasionally there are a couple of divergent transitions, seemingly nothing too severe. But my biggest issue is that the estimated parameter values are not what I would expect them to be. For instance, I had a condition where the participants were asked to be risk-averse, and another where they were supposed to be more risk-taking; in this model, this should translate into the alpha parameter being lower in the former condition and higher in the latter one. This was also confirmed by a maximum likelihood estimate of these parameters, which I am using as a sort of ‘ground truth’. Yet when I fit this model to data of a single participant from these two conditions (separately), this is very often not the case.

So, because I am new to all this, I was wondering if you have any ideas why that might be? I have tried playing around with the priors and the parameter limits in the Stan model, but they alone were not responsible for the estimated parameters to be ‘off’. Is there something obvious that I’m missing? For instance, do the parameters get transformed ‘under the hood’ by stan, meaning that the model is not sampling in the parameter space where I think it’s sampling from? Or have I simply specified something the wrong way, which is causing the model to misbehave?

Many thanks for any insights!