Hi all,
I’ve been working with the super computer; that sample of roughly 32 participants runs fine. No problem.
Then I move to all 1000 participants, 14 conditions (and four time points each of which is run separately). It takes about two days to run with each chain (4) using a separate core, and then I get the following error.
Error: cannot allocate vector of size 17.4 Gb
Error in FUN(X[[i]], ...) :
trying to get slot "mode" from an object (class "try-error") that is not an S4 object
Calls: stan ... sampling -> sampling -> .local -> sapply -> lapply -> FUN
In addition: Warning message:
In parallel::mclapply(1:chains, FUN = callFun, mc.preschedule = FALSE, :
1 function calls resulted in an error
Execution halted
Warning message:
system call failed: Cannot allocate memory
Unfortunately, I don’t know all of the details on the rstan version, etc., since it’s running from the super computer, but I can email the folks running the system and can ask them.
It seems like these comprise multiple errors: one of which is solved by reinstalling Rcpp (found here: Error: cannot allocate vector of size 13813.6 Gb)
The other errors seem more evasive. For instance, I received the same 17.4gb warning after the first attempt, and then I changed the output so that I was only pulling out the threshold and width parameters,
fit.w <- summary(fit.norm_subj.cond.t1, pars = c("w"))
fit.w <- fit.w$summary
write.csv(fit.w, "fit.w.1")
fit.m <- summary(fit.norm_subj.cond.t1, pars = c("m"))
fit.m <- fit.m$summary
write.csv(fit.m, "fit.m.1")
and I still got the same warning 17.4gb warning.
As a second question, I’m sure there are ways that I could speed the process along (faster than 2 days). Let me know if there are any resources you’d recommend.
Thanks!
Here is the code
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)
}
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<lower=0,upper=1> mean_beta; // Mean of beta prior for individual lapse rate
real<lower=0> betaEta; // Precision parameter of beta prior for individual lapse rate
vector<lower=0,upper=1>[n_total] psi;
real<lower=0> lapse_alpha;
real<lower=0> lapse_beta;
real m[n_subjects,n_levels];
real<lower=0> w[n_subjects,n_levels];
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);
}