This is our model that does seem to converge:
functions {
//Custom quantile function for lognormal distribution
matrix log_normal_qf(vector x, vector m, vector s) {
//takes a vector sampling the x-axis in terms of cumulative density, plus vectors of mu and sigma parameters
return exp(rep_matrix(s' * sqrt(2), size(x)) .* rep_matrix(inv_erfc(-2 * x + 2), size(m)) + rep_matrix(m', size(x)));
}
//Custom copula function
real normal_copula(real u, real v, real rho) {
real rho_sq = square(rho);
return (0.5 * rho * (-2. * u * v + square(u) * rho + square(v) * rho)) / (-1. + rho_sq) - 0.5 * log1m(rho_sq);
}
//Likelihood model for CASANDRE
real casandre_log(array[] matrix resps, vector guess, matrix sds, matrix sm, vector nconf, vector pconf) {
//declare and initialize log likelihood variable
real llhC;
llhC = 0;
//loop over participants to calculate log likelihood
for (n in 1:size(resps)) {
//declare the index-sizing variable
int m[rows(resps[n])];
//build the logical vector of non-zero trials
for (tt in 1:size(m)) {
m[tt] = sum(resps[n][tt]) == 1;
}
//declare the likelihood variable
matrix[rows(sm),4] lhC;
//check if subject has any non-zero data before running
if (sum(resps[n]) != 0) {
//declare incrementing and index variables
int t;
int ind[sum(m)];
//initialize incrementing variable
t = 1;
//declare and calculate mu parameters for normal cdfs
matrix[rows(sm),rows(sds)] avgs;
avgs = rep_matrix(sm[,n],rows(sds)).*rep_matrix(sds[,n]',rows(sm));
//loop over trials
for (tr in 1:rows(sm)){
//check if trial has any non-zero data before running
if (sum(resps[n][tr]) != 0) {
//declare sample vector
matrix[3,rows(sds)] raws;
//sample the cumulative density of normals along the transformed x-axis for each confidence bin
for (rws in 1:rows(sds)) {
raws[1,rws] = normal_cdf(-nconf[n],avgs[tr,rws],sds[rws,n]);
}
for (rws in 1:rows(sds)) {
raws[2,rws] = normal_cdf(0,avgs[tr,rws],sds[rws,n]);
}
for (rws in 1:rows(sds)) {
raws[3,rws] = normal_cdf(pconf[n],avgs[tr,rws],sds[rws,n]);
}
//declare the cumulative likelihood variable
vector[3] ratiodist;
//take the mean of the sampled densities
ratiodist[1] = mean(raws[1]);
ratiodist[2] = mean(raws[2]);
ratiodist[3] = mean(raws[3]);
//implement a lapse rate parameter to absorb unmodeled noise, and calculate likelihoods of each choice possibility in the trial
lhC[tr,1] = ((guess[n]/4) + ((1-guess[n])*ratiodist[1]));
lhC[tr,2] = ((guess[n]/4) + ((1-guess[n])*(ratiodist[2]-ratiodist[1])));
lhC[tr,3] = ((guess[n]/4) + ((1-guess[n])*(ratiodist[3]-ratiodist[2])));
lhC[tr,4] = ((guess[n]/4) + ((1-guess[n])*(1-ratiodist[3])));
//add this trial to the index
ind[t] = tr;
//increment the index
t = t + 1;
}
}
//multiply the choice matrix by the log of the likelihood matrix
llhC += sum(columns_dot_product(resps[n][ind], log(lhC[ind])));
}
}
//return the log likelihood for all data given the sampled parameters
return llhC;
}
//Partial sum function
real partial_sum_casandre_log(array[] matrix slice_n_resps, int start, int end, vector guess, matrix sds, matrix sm, vector nconf, vector pconf) {
//dynamically slice data into the log likelihood function
return casandre_log(slice_n_resps, guess[start:end], sds[:,start:end], sm[:,start:end], nconf[start:end], pconf[start:end]);
}
}
data {
//declare the number of subjects
int ns;
//declare the number of trials in the largest unpadded data set for each experiment
int nt_p;
int nt_v;
//declare the number of conditions and number of contrast levels in the perceptual domain
int ncond;
int ncont;
//declare the response matrix (one subject per array) for perceptual domain
array[ncond,ncont,ns] matrix[nt_p,4] resps_p;
//declare the response matrix (one per subject) for the value-based domain
array[ns] matrix[nt_v,4] resps_v;
//declare the experimental values matrices (one subject per matrix in each array) for perceptual domain
array[ncond,ncont] matrix[nt_p,ns] vals_p;
//declare the experimental values matrices (one subject per matrix) for value-based domain
matrix[nt_v,ns] vals_v;
//declare the number of x-axis samples and the vector sampling the x-axis in terms of cumulative density
int sampn;
vector[sampn] sampx;
}
parameters {
//Perceptual parameters
vector<lower=0,upper=1>[ns] guess_p; //lapse rate
matrix<lower=0>[ncont,ns] sens_p; //stimulus sensitivities
matrix[ncont,ns] crit_p; //decision criteria
matrix<lower=0>[ncond,ns] meta_p; //meta-uncertainties
matrix<lower=0>[ncond,ns] nconf_p; //negative confidence criteria
matrix<lower=0>[ncond,ns] pconf_p; //positive confidence criteria
//Value-based parameters
vector<lower=0,upper=1>[ns] guess_v; //lapse rate
vector<lower=0>[ns] sens_v; //stimulus sensitivities
vector<lower=0>[ns] meta_v; //meta-uncertainties
vector<lower=0>[ns] nconf_v; //negative confidence criteria
vector<lower=0>[ns] pconf_v; //positive confidence criteria
//Perceptual hyperparameters
vector[ncont] ssm_p; //mu hyperparameters of each contrast level's stimulus sensitivity lognormal prior
vector<lower=0>[ncont] sss_p; //sigma hyperparameters of each contrast level's stimulus sensitivity lognormal prior
vector[ncont] scm_p; //mu hyperparameters of each contrast level's decision criterion normal prior
vector<lower=0>[ncont] scs_p; //sigma hyperparameters of each contrast level's s decision criterion normal prior
vector[ncond] mum_p; //mu hyperparameters of each condition's meta-uncertainty lognormal prior
vector<lower=0>[ncond] mus_p; //sigma hyperparameters of each condition's meta-uncertainty lognormal prior
vector[ncond] nccm_p; //mu hyperparameters of each condition's negative confidence criterion lognormal prior
vector<lower=0>[ncond] nccs_p; //sigma hyperparameters of each condition's negative confidence criterion lognormal prior
vector[ncond] pccm_p; //mu hyperparameters of each condition's positive confidence criterion lognormal prior
vector<lower=0>[ncond] pccs_p; //sigma hyperparameters of each condition's positive confidence criterion lognormal prior
//Value-based hyperparameters
real ssm_v;
real<lower=1>sss_v;
real mum_v;
real<lower=1>mus_v;
real nccm_v;
real<lower=1> nccs_v;
real pccm_v;
real<lower=1> pccs_v;
//Correlation parameter
vector<lower=-1,upper=1>[ncond] rho;
}
model {
//Perceptual hyperpriors
ssm_p ~ normal(0,1);
sss_p ~ lognormal(0,1);
scm_p ~ normal(0,1);
scs_p ~ lognormal(0,1);
mum_p ~ normal(0,1);
mus_p ~ lognormal(0,1);
nccm_p ~ normal(0,1);
nccs_p ~ lognormal(0,1);
pccm_p ~ normal(0,1);
pccs_p ~ lognormal(0,1);
//Value-based hyperpriors
ssm_v ~ normal(0,1);
sss_v ~ lognormal(0,1);
mum_v ~ normal(0,1);
mus_v ~ lognormal(0,1);
nccm_v ~ normal(0,1);
nccs_v ~ lognormal(0,1);
pccm_v ~ normal(0,1);
pccs_v ~ lognormal(0,1);
//Perceptual priors
guess_p ~ beta(1,197.0/3.0);
for (cont in 1:ncont) {
sens_p[cont] ~ lognormal(ssm_p[cont],sss_p[cont]);
crit_p[cont] ~ normal(scm_p[cont],scs_p[cont]);
}
for (cond in 1:ncond) {
meta_p[cond] ~ lognormal(mum_p[cond],mus_p[cond]);
nconf_p[cond] ~ lognormal(nccm_p[cond],nccs_p[cond]);
pconf_p[cond] ~ lognormal(pccm_p[cond],pccs_p[cond]);
}
//Value-based priors
guess_v ~ beta(1,197.0/3.0);
sens_v ~ lognormal(ssm_v,sss_v);
meta_v ~ lognormal(mum_v,mus_v);
nconf_v ~ lognormal(nccm_v,nccs_v);
pconf_v ~ lognormal(pccm_v,pccs_v);
//declare the transformed x-axis matrices
array[ncond] matrix[sampn,ns] xtrans_p;
array[ncond] matrix[sampn,ns] sds_p;
//Perceptual likelihood model (with local variable calculations)
for (cond in 1:ncond) {
//calculate the transformed x-axis matrices
xtrans_p[cond] = log_normal_qf(sampx,-0.5*log1p(meta_p[cond]'.^2),sqrt(log1p(meta_p[cond]'.^2)));
sds_p[cond] = 1./xtrans_p[cond];
for (cont in 1:ncont) {
//check if condition and contrast level combination has non-zero data
if (sum(abs(vals_p[cond][cont])) != 0) {
//declare matrices for calculating each contrast level's transformed data
array[ncont] matrix[nt_p,ns] om_p;
array[ncont] row_vector[ns] cm_p;
array[ncont] matrix[nt_p,ns] sm_p;
//orientations and decision criteria in terms of standard deviations
om_p[cont] = vals_p[cond][cont].*rep_matrix(sens_p[cont],nt_p);
cm_p[cont] = crit_p[cont].*sens_p[cont];
//transform the data
sm_p[cont] = om_p[cont]-rep_matrix(cm_p[cont],nt_p);
//hand the data and parameters to the reduce sum wrap function for dynamic within-chain parallelization
target += reduce_sum(partial_sum_casandre_log, resps_p[cond][cont], 1, guess_p, sds_p[cond], sm_p[cont], nconf_p[cond]', pconf_p[cond]');
}
}
}
//Value-based likelihood model (with local variable calculations)
//declare the transformed x-axis matrices
matrix[sampn,ns] xtrans_v;
matrix[sampn,ns] sds_v;
//calculate the transformed x-axis matrices
xtrans_v = log_normal_qf(sampx,-0.5*log1p(meta_v.^2),sqrt(log1p(meta_v.^2)));
sds_v = 1./xtrans_v;
//declare matrices for calculating transformed data
matrix[nt_v,ns] om_v;
matrix[nt_v,ns] sm_v;
//orientations in terms of standard deviations
om_v = vals_v.*rep_matrix(sens_v',nt_v);
//transform the data
sm_v = om_v;
//hand the data and parameter to the reduce sum wrap function for dynamic within-chain parallelization
target += reduce_sum(partial_sum_casandre_log, resps_v, 1, guess_v, sds_v, sm_v, nconf_v, pconf_v);
//Correlation copula
vector[ncond] y_p;
real y_v;
for (n in 1:ns) {
for (cond in 1:ncond) {
y_p[cond] = inv_Phi(lognormal_cdf(meta_p[cond,n],mum_p[cond],mus_p[cond]));
y_v = inv_Phi(lognormal_cdf(meta_v[n],mum_v,mus_v));
target += normal_copula(y_p[cond],y_v,rho[cond]);
}
}
}