There are ~1000 participants and fourteen conditions for which I am estimating thresholds via psychometric functions. Intensity is the RW for the tasks, which lengthens as participants answer wrong and decreases (as participants answer right) via an adaptive paradigm.
I have access to a super computer, but no job can last longer than 48 hours; hence, I’m taking on trying to run using multi-threads, as one person on this forum suggested.
I’ve looked through the documentation and think that I am pretty close to having everything ready to go, but I’m wondering whether in the space marked HERE, I should use merely psi (as in the “model” section of the code), or use the entire psi function (as listed in “transformed parameters”).
I might also be way off. Hopefully not.
Thanks much!
functions {
real partial_sum_lpmf(int[] slice_n_subjects,
int start, int end,
int[] n_levels,
vector intensity,
vector lapse,
real mum,
real muw,
vector fA,
vector sA,
vector sB,
vector fsA) {
return binomial_logit_lupmf(slice_n_subjects |
n_levels[start:end],
**HERE** );
}
}
data{ // Everything that must be input by the user
int<lower=4> n_total; // Total number of trials to analyze
int<lower=2> n_subjects; // Number of unique observers
int<lower=2> n_levels; // Number of levels of Factor
real intensity[n_total]; // Intensity of the stimulus on each trial
int<lower=0> subject[n_total]; // Observer on each trial
int<lower=0> level[n_total]; // Level of Factor on each trial
int<lower=0,upper=1> correct[n_total]; // Whether the response was correct (1) on each trial
real<lower=0,upper=1> chance_performance; // Chance performance for experiment (e.g., 1/n for n-AFC)
int<lower=1> grainsize;
}
transformed data{
int df_levels_sA; // Number of degrees of freedom in the interaction
int n_levels_sA; // Number of levels in the interaction
df_levels_sA = (n_levels - 1) * (n_subjects - 1);
n_levels_sA = n_levels * n_subjects;
}
parameters{
vector<lower=0,upper=1>[n_subjects] lapse; // Observer's lapse rate
real mum;
real muw;
vector[n_levels-1] fA;
vector[n_subjects-1] sA;
vector[n_subjects-1] sB;
//matrix[n_subjects-1,n_levels-1] fsA;
vector[(n_subjects - 1) * (n_levels-1)] fsA;
}
transformed parameters {
real m[n_subjects,n_levels];
real w[n_subjects,n_levels];
real lapse_alpha;
real lapse_beta;
vector [n_total] psi;
{
real mean_beta; // Mean of beta prior for individual lapse rate
real betaEta; // Precision parameter of beta prior for individual lapse rate
real threshold[n_total];
real width[n_total];
vector[n_levels] factor_alpha = append_row(fA, -sum(fA));
vector[n_subjects] subject_alpha = append_row(sA, -sum(sA));
vector[n_subjects] subject_beta = append_row(sB, -sum(sB));
matrix[n_subjects,n_levels] interaction_alpha;
interaction_alpha = rep_matrix(0, n_subjects, n_levels);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_alpha[sj, l] = fsA[(sj - 1) * (n_levels - 1) + l];
}
interaction_alpha[sj, n_levels] = -1 * sum(interaction_alpha[sj,]);
}
for (l in 1:n_levels) {
interaction_alpha[n_subjects,l] = -1 * sum(interaction_alpha[,l]);
}
for (sj in 1:n_subjects) {
for (l in 1:n_levels) {
m[sj,l] = mum + factor_alpha[l] + subject_alpha[sj] + interaction_alpha[sj, l];
w[sj,l] = exp(muw + subject_beta[sj]);
}
}
for (tr in 1:n_total) {
threshold[tr] = m[subject[tr],level[tr]];
width[tr] = w[subject[tr],level[tr]];
psi[tr] = chance_performance
+ (1 - lapse[subject[tr]] - chance_performance)
* inv_logit((intensity[tr]-threshold[tr])/width[tr]);
if (is_nan(psi[tr])) {
print("lapse = ", lapse[subject[tr]]);
print("width = ", width[tr]);
print("threshold = ", threshold[tr]);
}
}
mean_beta = 0.01;
betaEta = 100;
lapse_alpha = mean_beta * betaEta;
lapse_beta = (1-mean_beta) * betaEta ;
}
}
model {
//mean_beta ~ beta(1,60);
//betaEta ~ gamma(1,0.01);
fsA ~ normal(0, inv(sqrt(1 - inv(n_levels_sA)))); // Produces a standard normal on on interaction_alpha
fA ~ normal(0, inv(sqrt(1 - inv(n_levels))));
sA ~ normal(0, inv(sqrt(1 - inv(n_subjects)))); // Produces a standard normal on on interaction_alpha
sB ~ normal(0, inv(sqrt(1 - inv(n_subjects)))); // Produces a standard normal on on interaction_alpha
muw ~ gamma(2,.5);
lapse ~ beta(lapse_alpha,lapse_beta);
correct ~ bernoulli(psi);
target += reduce_sum(partial_sum_lupmf, n_subjects, grainsize,
n_levels, intensity, lapse, mum, muw, fA, sA, sB, fsA);
}