Thats great - thank you so much. Yes, the two models are dependent (I was trying to make my question shorter - sorry for that). What Iâ€™m trying to do is fit choices and reaction-times simultaneity using RL Q-learning.

```
transformed parameters {
//declare variables and parameters
vector<lower=0, upper=1>[Nsubjects] alpha;
vector<lower=0, upper=10>[Nsubjects] beta;
for (subject in 1:Nsubjects) {
alpha[subject] = Phi_approx(mu_aux[1] + sigma_aux[1] * alpha_individual_aux[subject]);
beta[subject] = Phi_approx(mu_aux[2] + sigma_aux[2] * beta_individual_aux[subject]) * 10;
mu[subject] = Phi_approx(mu_aux[3] + sigma_aux[3] * mu_individual_aux[subject]) ;
sigma[subject] = Phi_approx(mu_aux[3] + sigma_aux[3] * sigma_individual_aux [subject]) ;
lambda[subject] = Phi_approx(mu_aux[3] + sigma_aux[3] * lambda_individual_aux[subject]) ;
}
}
model {
// population level priors (hyper-parameters)
mu_aux ~ normal(0, 1);
sigma_aux ~ normal(0, 0.2);
// individual level priors (subjects' parameters)
alpha_individual_aux~ normal(0, 1);
beta_individual_aux ~ normal(0, 1);
mu_individual_aux ~ normal(0, 1);
sigma_individual_aux ~ normal(0, 1);
lambda_individual_aux ~ normal(0, 1);
//Likelihood function per subject per trial
for (subject in 1:Nsubjects){
vector[Narms] Qcard;
Qcard=Qvalue_initial;
for (trial in 1:Ntrials_per_subject[subject]){
//likelihood function
log_lik += categorical_logit_lpmf(action[subject, trial] | beta[subject] * Qcard);
log_lik += exp_mod_normal_lpdf(rt | mu+beta2*(Qcard[1]-Qcard[2)], sigma, lambda)
//Qvalues update
Qcard[action[subject,trial]] += alpha[subject] * (reward[subject,trial] - Qcard[action[subject,trial]]);
}
}
}
'''
```