I’m receiving the following error:
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: validate transformed params: psi[8742] is -nan, but must be greater than or equal to 0 (in 'model3a6550c0a439_thresholds' at line 30)
Chain 4:
Chain 4: Initialization between (-2, 2) failed after 100 attempts.
Chain 4: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
error occurred during calling the sampler; sampling not done
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
I receive the error for all four chains. I did find this post "https://discourse.mc-stan.org/t/error-evaluating-the-log-probability-at-the-initial-value/5254," and I checked my priors and all have a lower bound at zero. It’s possible that I need to set the prior on fsA, but then it would be a matrix prior, and I’m uncertain what I should set that prior to be.
R version is 3.6.0. I’m running Stan through rStan 2.19.3. StanHeaders is 2.19.2.
This is my Stan 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;
}
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<lower=0> width[n_total];
vector[n_levels] factor_alpha;
vector[n_subjects] subject_alpha;
vector[n_subjects] subject_beta;
matrix[n_subjects,n_levels] interaction_alpha;
for (l in 1:(n_levels-1)) {
factor_alpha[l] = fA[l];
}
factor_alpha[n_levels] = -1*sum(fA)/(n_levels-1);
for (l in 1:(n_subjects-1)) {
subject_alpha[l] = sA[l];
subject_beta[l] = sB[l];
}
subject_alpha[n_subjects] = -1*sum(sA)/(n_subjects - 1);
subject_beta[n_subjects] = -1*sum(sB)/(n_subjects - 1);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_alpha[sj,l] = fsA[sj,l];
}
interaction_alpha[sj,n_levels] = -1 * sum(fsA[sj,]);
}
for (l in 1:n_levels-1) {
interaction_alpha[n_subjects,l] = -1 * sum(fsA[,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] = 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] = (1-lapse[ subject[tr] ])*( (1-chance_performance)*inv_logit(4.4/width[tr] * (intensity[tr]-threshold[tr])) + chance_performance) + chance_performance*lapse[ subject[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);
fA ~ normal(0,1);
sA ~ normal(0,1);
sB ~ normal(0,1);
muw ~ gamma(2,.5);
lapse ~ beta(lapse_alpha,lapse_beta);
correct ~ bernoulli(psi);
}
This is the code I use from R:
dat <- ace.threshold.t2
dat <- subset(dat, !is.na(norm))
dat$condition <- factor(dat$condition)
dat$pid <- factor(dat$pid)
nTotal <- dim(dat)[1]
nCond <- length(unique(dat$condition))
nSubj <- length(unique(dat$pid))
intensity <- dat$norm
condition <- as.numeric(dat$condition)
pid <- as.numeric(dat$pid)
correct <- dat$correct_button == "correct"
chancePerformance <- 1/2
stanData <- list(n_total=nTotal, n_levels=nCond, n_subjects = nSubj, subject = pid, intensity=intensity, level=condition, correct=correct, chance_performance=chancePerformance)
fit.norm.diff <- stan(file="thresholds.stan", data=stanData, cores = 4)
Thank you very much!
James
Here is a link to a snippet of the data:
