Failure to recover delta

Modeling
,
1

I have generated some data using the rwiener function and the following parameters:
[s1] Alpha: 0.940 Beta: 0.354 Delta: 0.104 Tau: 0.371
[s2]Alpha: 1.01 Beta:0.183 Delta: 0.105 Tau: 0.373

However, when I put the data through my basic drift diffusion model of 4 parameters, which uses the function wiener_lpdf, it fails to recover the true delta value, and gives me the following output:

Mean: [s1]0.5925061
[s2]0.4015847
sd: [s1]0.305563750
[s2]0.222346652
2.5%: [s1]0.13145667
[s2]0.09094771
25%: [s1]0.3604417
[s2]0.2319065
50%: [s1]0.5512207
[s2]0.3632976
75%: [s1]0.7792915
[s2]0.5300637
97.5%: [s1]1.2819739
[s2]0.9295888
n_eff: [s1]7046.331
[s2]8127.745
R-hat: [s1]0.9998943
[s2]0.9998650

The parameters are specified like this

parameters {
    real raw_alpha[N_SUBJECTS];
    real raw_beta[N_SUBJECTS];
    real raw_delta[N_SUBJECTS];
    real raw_tau[N_SUBJECTS];
}
transformed parameters {
    vector<lower=0>[N_SUBJECTS] alpha_boundary_separation;
    vector<lower=0, upper=1>[N_SUBJECTS] beta_initial_bias;  // towards upper boundary
    vector<lower=0>[N_SUBJECTS] delta_drift_rate;  // towards upper boundary
    vector<lower=RT_LOWER_BOUND, upper=max(min_rt)>[N_SUBJECTS] tau_nondecision_time;

    for (s in 1:N_SUBJECTS) {
        
        alpha_boundary_separation[s] = exp(
            raw_alpha[s]
        );
        
        beta_initial_bias[s] = Phi_approx(          
            raw_beta[s]
        );
        
        delta_drift_rate[s] = exp(  
            raw_delta[s]
        );
        
        tau_nondecision_time[s] = Phi_approx(
            raw_tau[s]
        ) * (min_rt[s] - RT_LOWER_BOUND) + RT_LOWER_BOUND; 
   } 
}

And these are the priors (paremetized like this to allow bridgesampling)

model {
    target += normal_lpdf(raw_alpha | 0, 1);
    target += normal_lpdf(raw_beta | 0, 1);
    target += normal_lpdf(raw_delta | 0, 1);
    target += normal_lpdf(raw_tau | 0, 1);
}

I have no idea of what could be causing the problem, because diagnostics look fine. I’ve tried fixing all parameters but delta, and giving the parameters specific informative priors, but the estimation was still not right.

Thank you for reading. Any help will be much appreciated!

2

I have no idea what’s going on with this density. I thought it might be your parameterization, since I understand there are many of them for this density. So I tried just generating data with Stan and then fitting that so the parameterizations matched. First, it has a really hard time sampling without priors, so I threw some exponential priors on, which are the right scale. That sped up sampling, but we don’t get anything like the coverage we’d expect. Here’s the code I used.

Simulate data: wiener.stan

transformed data {
  real alpha = 0.94;
  real beta = 0.354;
  real delta = 0.104;
  real tau = 0.371;
}
parameters {
  real<lower=tau> y;
}
model {
  y ~ wiener(alpha, tau, beta, delta);
}

Fit: wiener-fit.stan

data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real<lower=0> alpha;
  real<lower=0, upper=max(y)> tau;
  real<lower=0, upper=1> beta;
  real<lower=0> delta;
}
model {
  y ~ wiener(alpha, tau, beta, delta);
  alpha ~ exponential(1);
  tau ~ exponential(1);
  delta ~ exponential(1);
}

R simulation (cmdstanr)

library('cmdstanr')

modelsim <- cmdstan_model('wiener.stan')
fitsim <- modelsim$sample(thin=20, parallel_chains=4)
d <- fitsim$draws('y')
ysim <- c(d[ , 1, 1], d[ , 2, 1], d[, 3, 1], d[, 4, 1])
N <- length(ysim)

model <- cmdstan_model('wiener-fit.stan')
fit <- model$sample(data = list(N = N, y = ysim), parallel_chains=4)

Sampling is fast and gives good ESS estimates, but the actual posterior is nowhere near calibrated:

> fit
 variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
    lp__  90.38  90.71 1.39 1.20 87.68 91.99 1.00     1352     1751
    alpha  2.51   2.30 0.80 0.54  1.71  3.97 1.00     1286     1135
    tau    0.36   0.36 0.02 0.02  0.33  0.38 1.00     1311     1054
    beta   0.68   0.68 0.08 0.09  0.55  0.82 1.00     1195     1111
    delta  2.61   2.60 0.25 0.25  2.21  3.04 1.00     1468     1743

So maybe someone can see a mistake in what I did or explain why this doesn’t work.