Chain 1 Rejecting initial value: Chain 1 Error evaluating the log probability at the initial value. Chain 1 Exception: Exception: wiener_lpdf: Random variable = 0.560084, but must be greater than nondecision time = 2.48104

Interfaces CmdStan
1

Hello everyone, I encountered such a problem when running the cmdstanr package, I have been thinking about it for a long time and don’t know how to solve it, I hope someone can help me take a look at it if you have time, I appreciate it.Here is my code.Thank you!

load("H:/L2dat.RData")
data <- Level2dat
summary(data)
# Resp: Right (1), Left (2)
# Cond: Easy (1), Medium (2), Hard (3)
# Stim: Right (1), Left (2)
# -> Use response coding
library(cmdstanr)
#set_cmdstan_path("~/Psychology/DDM/cmdstan/")  # set correspondingly
Ncnds = 3
chains = 1
# initial values
init.stan = function(){
  L = list()
  for (i in 1:chains) {
    L[[i]] = list(
      a = 1 + runif(1, -0.2, 0.2),
      v = 2.5 + runif(Ncnds, -0.5, 0.5),
      w = 0.5 + runif(1, -0.05, 0.05),
      t0 = 0.2 + runif(1, -0.1, 0.1),   
      sv = 1 + runif(1, -0.1, 0.1),
      sw = 0.1 + runif(1, -0.05, 0.05),
      st0 = 0.2 +runif(1, -0.05, 0.05)
    )
  }
  return(L)
}
# Compiling model
mc <- cmdstan_model("H:/level2/level2.stan", cpp_options = list(stan_threads = TRUE))
stan.data = list(
  N = nrow(data),
  Ncnds = 3,
  rt = data$RT,
  stim = as.numeric(data$Stim),
  resp = as.numeric(data$Resp),
  cnd = as.numeric(data$Cond)
)
outdir = "fits"
dir.create(outdir, recursive = T, showWarnings = FALSE)
m <- mc$sample(data=stan.data,
               init = init.stan(),
               max_treedepth = 5,  # reduces number of function calls/iteration
               refresh = 25,
               iter_sampling = 250,
               iter_warmup = 250,
               chains = chains,
               parallel_chains = 1,
               threads_per_chain = 8,  # set accordingly
               output_dir = outdir,
               output_basename = "fit",
               save_warmup = TRUE
)

Here’s the STAN code.

functions {
real partial_sum_fullddm(array[] real rt_slice, int start, int end,
  real a, array[] real v, real w, real t0, real sv, real sw, real st0,
  array[] int stim, array[] int resp, array[] int cnd) {
    real ans = 0;
    for (i in start:end) {
      if (stim[i] == 2) {  // stimulus "Left": v_L=v, w_L=w
        if (resp[i] == 2) {  // response "Left" -> v'=v_L=v, w'=w_L=w
          ans += wiener_lpdf(rt_slice[i+1-start] | a, v[cnd[i]], w, t0, sv,  sw, st0);
        } else {  // response "Right" -> switch boundary: v'=-v_L=-v, w'=1-w_L=1-w
          ans += wiener_lpdf(rt_slice[i+1-start] | a, -v[cnd[i]], 1-w, t0, sv,  sw, st0);
        }
      } else {  // stimulus "Right": v_R=-v, w_R=w
        if (resp[i] == 1) {  // response "Right" -> switch boundary v'=-v_R=v, w'=1-w_R=1-w
          ans += wiener_lpdf(rt_slice[i+1-start] | a, v[cnd[i]], 1-w, t0, sv, sw, st0);
        } else {  // response "Left" -> v'=v_R=-v, w'=w_R=w
          ans += wiener_lpdf(rt_slice[i+1-start] | a, -v[cnd[i]], w, t0, sv,  sw, st0);
        }
      }
    }
    return ans;
  }
}

data {
int<lower=0> N;  // number of trials
int<lower=1> Ncnds;  // number of conditions: 3
array[N] real rt;  // response times (seconds)
array[N] int<lower=1, upper=Ncnds> cnd;  // condition (1, 2, 3)
array[N] int<lower=1, upper=2> stim;  // stimuli (1, 2)
array[N] int<lower=1, upper=2> resp;  // responses (1, 2)
}

parameters {
// parameter artificially constrained following Table B1 (Heathcote's priors)
real<lower=0, upper=2> a;  // boundary separation (a>0)
real<lower=0.001, upper=0.999> w;  // relative starting point
array[Ncnds] real v;  // mean drift for each condition type
real<lower=0> t0; // mean non-decision-time

real<lower=0, upper=1> sw;  // variability of relative starting point
real<lower=0> sv;  // variability of drift-rate
real<lower=0> st0; // variability of non-decision-time
}

model {
// priors
a ~ normal(1, 1)T[0,2];
v ~ normal(2.5, 3);
t0 ~ normal(0.3, 1)T[0,];
w ~ normal(0.5, 1)T[0,1];

sv ~ normal(1, 3)T[0,];
st0 ~ normal(0.5, 1)T[0,1];
sw ~ normal(0.5, 1)T[0,1];

target += reduce_sum(partial_sum_fullddm, rt, 1,
  a, v, w, t0, sv, sw, st0, stim, resp, cnd);
}