New wiener distribution estimating parameters strangely from RWiener package

[Joseph Burling]

I'm trying to figure out why the parameters estimated when using the new wiener_log function differ quite dramatically from simulated data using the RWiener package. I generated some random response time values using the rwiener() function from the RWiener package, then tried to estimate the parameters from the simulated data in Stan. The Stan manual refers me to this R package so I assume the parameterization from the wiener distribution is the same in both Stan and RWiener. I'm uploading some example code for testing.

For example, simulated data from the RWiener package using parameters alpha = 3.25, tau = 0.25, beta = 0.4, delta = 0.65 will give me the following in Stan, 7.25, 0.22, 0.71, 0.93. This occurs even with a larger N and setting different priors on these four parameters. Any ideas?

EDIT: in case attachment isn't found

library(rstan)

library(loo)

library(RWiener)

# if running chains in parallel

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

# wiener parameters

a <- 3.25

t <- 0.25

b <- 0.4

v <- 0.65

n <- 150

# random RT values from RWiener using parameters above

rt <- rwiener(

    n = n,

    alpha = a,

    tau = t,

    beta = b,

    delta = v

)

# make stan data list

stan_data <- list(rt = rt$q,

                  acc = ifelse(rt[, "resp"] == "upper", 1, 0),

                  n = nrow(rt))

# number of chains to run

n_chains <- 5

# use inits from above dist. parameters

stan_inits <- lapply(1:n_chains, function(i) {

    list(

        a_mu = max(c(0.0001, a)),

        t_mu = max(c(0.0001, t)),

        v_mu = max(c(0.0001, v)),

        a_raw = 0.01,

        t_raw = 0.01,

        v_raw = 0.01,

        a_sd = 1,

        t_sd = 1,

        v_sd = 1,

        b = b

    )

})

# STAN code

stan_model <-

"

data {

    int<lower=1>                n;

    vector[n]                   rt;

    vector<lower=0,upper=1>[n]  acc;

}

parameters {

    real<lower=0,upper=10> a_mu;

    real<lower=0,upper=3>  a_sd;

    real                   a_raw;

    real<lower=0,upper=3>  t_mu;

    real<lower=0,upper=2>  t_sd;

    real                   t_raw;

    real<lower=0,upper=5>  v_mu;

    real<lower=0,upper=2>  v_sd;

    real                   v_raw;

    real<lower=0,upper=1>  b;     // <- beta

}

transformed parameters {

    real a; // <- alpha

    real t; // <- tau

    real v; // <- delta

    a <- a_mu + a_sd * a_raw; // a ~ normal(a_mu, a_sd)

    t <- t_mu + t_sd * t_raw; // t ~ normal(t_mu, t_sd)

    v <- v_mu + v_sd * v_raw; // v ~ normal(v_mu, v_sd)

}

model {

    a_raw ~ normal(0,1);

    t_raw ~ normal(0,1);

    v_raw ~ normal(0,1);

    b ~ beta(2, 2);

    // likelihood

    for (i in 1:n) {

        if (acc[n] == 1) {

            rt[i] ~ wiener(a, t, b, v);    // upper bound

        } else {

            rt[i] ~ wiener(a, t, 1-b, -v); // lower bound

        }

    }

}

generated quantities {

    vector[n] log_lik;

    for (i in 1:n) {

        if (acc[n] == 1) { // upper

            log_lik[i] <- wiener_log(rt[i], a, t, b, v);

        } else { // lower

            log_lik[i] <- wiener_log(rt[i], a, t, 1-b, -v);

        }

    }

}

"

# run model

stan_fit <-

    stan(

        model_name = "diffusion_test",

        model_code = stan_model,

        data = stan_data,

        init = stan_inits,

        chains = n_chains

    )

# print stan model parameters

print(stan_fit)

# print elpd_loo

print(loo(extract_log_lik(stan_fit)), digits = 3)

# print RWiener likelihood

print(wiener_likelihood(c(a,t,b,v), rt))

# Plots -------------------------------------------------------------------

# stan traceplot

traceplot(stan_fit)

# plot distributions from RWiener

plot_rwiener <- function(rt) {

    rt_u <- rt$q[rt$resp == "upper"]

    rt_l <- rt$q[rt$resp == "lower"]

    xlim <- max(c(rt_u, rt_l))

    rt_seq <- seq(0, xlim, .005)

    d_rt_u <-

        dwiener(rt_seq, a, t, b, v, resp = rep("upper", length(rt_seq)))

    d_rt_l <-

        dwiener(rt_seq, a, t, b, v, resp = rep("lower", length(rt_seq)))

    ylim <- max(c(d_rt_u, d_rt_l))

    par(mfrow = c(3, 1))

    plot(

        x = rt_seq,

        y = d_rt_u,

        type = "l",

        xlim = c(0, xlim),

        ylim = c(-ylim, ylim),

        main = "Density of upper/lower responses",

        ylab = "density",

        xlab = "quantile",

        col = "green"

    )

    lines(x = rt_seq, y = -d_rt_l, col = "red")

    xbrks <- seq(0, xlim, length.out = 30)

    hist(

        rt_u,

        freq = FALSE,

        breaks = xbrks,

        xlim = c(0, xlim),

        ylim = c(0, ylim * 2),

        main = "upper",

        xlab = paste("N =", length(rt_u))

    )

    hist(

        rt_l,

        freq = FALSE,

        breaks = xbrks,

        xlim = c(0, xlim),

        ylim = c(0, ylim * 2),

        main = "lower",

        xlab = paste("N =", length(rt_l))

    )

}

plot_rwiener(rt)

[Guido Biele]

is there a specific reason you don't simply specify a prior for a (and v and t)?*

it looks as if your statistical model is more complex than the true generative model.

i.e. when you generate the data 

a <- 3.25 

but in your stan model 

a <- a_mu + a_sd * a_raw.

so  I think 

a ~ normal(0,2);

after specifying real<lower=0,upper=10> a; 

in the parameters block and removing a_mu, a_raw, and a_sd should work better.

(same for v and t).

maybe you can try this simpler model first

(I had no issues to recover the correct wiener parameter values 

when testing such a simple model a few weeks ago).

personally, I would try this simpler model first and see if I can recover parameters, 

before trying the more complex model. but if you do this and you use 

a <- a_mu + a_sd * a_raw;

in stan, then you might also want to generate the data in R accordingly

I also would use more informative priors than the (implicit) flat priors you used for 

a_mu, a_sd, ... you have in your model.

guido

*(maybe this is a preparation to run the full DDM with between trial variability for v,a,t?)

[Joseph Burling]

I had the priors on the wiener parameters setup for a hierarchical model later down the line, which is why you might see the mu, sd, raw, etc..

However, even when setting normal priors constrained on some range as suggested, I still get weird estimates, in this case alpha=5.27, tau=0.12, v=0.92, b=0.58. Next I tried just using uniform priors on a, t, v, b and I get 7.27, 0.11, 1.04, 0.66 respectively (see code below).

I guess I'm not sure what the proper constraints should be, but I've seen normal priors used before but with the JAGS version of the wiener distribution.

        a = a,

        t = t,

        v = v,

        b = b

    )

})

# STAN code

stan_model <-

"

data {

    int<lower=1>                n;

    vector[n]                   rt;

    vector<lower=0,upper=1>[n]  acc;

}

parameters {

    real<lower=0,upper=10> a;

    real<lower=0,upper=5>  t;

    real<lower=0,upper=5>  v;

    real<lower=0,upper=1>  b;

}

model {

    //a ~ normal(0, 2);

    //t ~ normal(0, 2);

    //v ~ normal(0, 2);

    //b ~ beta(2, 2);

[Guido Biele]

Hi,

I don't have time to look at the model in detail. 

but one thing that might be an issue i that the 150 trials 

you have are not a lot to allow a reliable estimation of parameters, I believe.

Here is some code for a simple parameter recovery test that shows that 

one can recover wfpt distribution parameters with acceptable accuracy with stan:

https://github.com/gbiele/stan_wiener_test

maybe you can compare this with your code to see what the problem might be.

One issue I have not been able to solve is that using the wiener in stan

requires lots of RAM. For the simple test in the linked code, 2500 trials

(500 per condition) is already too much for a machine with 8GB RAM.

(using rstan, I haven't tested that in cmdstan yet)

guido

[Joseph Burling]

Hi,

Thanks for providing me with your model to look at. The main differences I noticed was setting the max of tau to the minimum response time from the data, and using a half-cauchy prior with small scale. You also allowed delta to be negative. I made all these changes and still the parameters don't match up with the original I used to sample the data. The chains are mixing well, so that doesn't seem to be the issue. I then tried sampling 1000 trials instead of 150, still no luck. Below is my final attempt.

library(rstan)

library(loo)

library(RWiener)

# if running chains in parallel

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

# wiener parameters

a <- 3.25

t <- 0.25

b <- 0.4

v <- 0.65

n <- 1000

# random RT values from RWiener using parameters above

rt <- rwiener(

    n = n,

    alpha = a,

    tau = t,

    beta = b,

    delta = v

)

# make stan data list

stan_data <- list(rt = rt$q,

                  acc = ifelse(rt[, "resp"] == "upper", 1, 0),

                  n = nrow(rt),

                  min_t = min(rt$q))

# number of chains to run

n_chains <- 5

# use inits from above dist. parameters

stan_inits <- lapply(1:n_chains, function(i) {

    list(

        a = a,

        t = t,

        v = v,

        b = b

    )

})

# STAN code

stan_model <-

"

data {

    int<lower=1>                n;

    vector[n]                   rt;

    vector<lower=0,upper=1>[n]  acc;

    real<lower=0>               min_t;

}

parameters {

    real<lower=0,upper=6>  a;

    real<lower=.1,upper=min_t>  t;

    real<lower=-2,upper=2> v;

    real<lower=0,upper=1>  b;

}

model {

    a ~ normal(0, 2);

    t ~ cauchy(0, 0.1);

    v ~ normal(0, 2);

    b ~ beta(1.25, 1.25);

print(lapply(extract(stan_fit), mean))

# Plots -------------------------------------------------------------------

# stan traceplot

#traceplot(stan_fit)

#plot_rwiener(rt)

[Joseph Burling]

I was able to actually recover the parameters in STAN generated from the RWiener package. It turns out I made a dumb mistake when indexing accuracy for the wiener likelihood. Just uploading correct code here just in case.

# install these packages if not already installed, run:

#    install.packages(c("rstan", "loo", "RWiener"))

library(rstan) # need to load namespace due to Rcpp dependencies

# for running multiple chains in parallel, uncomment next two lines

#rstan::rstan_options(auto_write = TRUE)

#options(mc.cores = parallel::detectCores())

# wiener parameters to try and estimate in STAN

a <- 3

t <- 0.25

b <- 0.33

v <- 1

n <- 50

# random RT values from RWiener package using parameters from above

rt <- RWiener::rwiener(

    n = n,

    alpha = a,

    tau = t,

    beta = b,

    delta = v

)

# make stan data list

stan_data <- list(rt = rt$q,

                  acc = ifelse(rt[, "resp"] == "upper", 1, 0),

                  n = nrow(rt),

                  min_t = min(rt$q),

                  max_t = max(rt$q))

# number of chains to run

n_chains <- 3

# use initial values from above weiner parameters

stan_inits <- lapply(1:n_chains, function(i) {

    list(

        a = a,

        t = t,

        v = v,

        b = b

    )

})

# STAN model code

stan_model <-

"

data {

    int<lower=1>                n;

    vector[n]                   rt;

    vector<lower=0,upper=1>[n]  acc;

    real<lower=0>               min_t;

    real<lower=0>               max_t;

}

parameters {

    real<lower=0,upper=100>          a; // speed-acc tradeoff / boundary width

    real<lower=0.1,upper=min_t>      t; // motor latency / non-decision time

    real<lower=-max_t,upper=max_t>   v; // processing speed / drift rate

    real<lower=0,upper=1>            b; // response bias / boundary bias

}

model {

    a ~ gamma(0.75, 0.25);

    v ~ normal(0, max_t/3);

    b ~ beta(1.25, 1.25);

    // implied: t ~ uniform(0.1, mininum RT)

    // likelihood

    for (i in 1:n) {

        if (acc[i] == 1) {

            rt[i] ~ wiener(a, t, b, v);    // upper bound

        } else {

            rt[i] ~ wiener(a, t, 1 - b, -1 * v); // lower bound

        }

    }

}

generated quantities {

    vector[n] log_lik;

    for (i in 1:n) {

        if (acc[i] == 1) { // upper boundary

            log_lik[i] <- wiener_log(rt[i], a, t, b, v);

        } else {           // lower boundary

            log_lik[i] <- wiener_log(rt[i], a, t, 1 - b, -1 * v);

        }

    }

}

"

# Fit STAN model ----------------------------------------------------------

stan_fit <-

    rstan::stan(

        model_name = "diffusion_test",

        model_code = stan_model,

        data = stan_data,

        init = stan_inits,

        chains = n_chains,

        iter = n_chains * 1000

    )

# print stan model parameters

print(stan_fit)

# print elpd_loo, expected log-prob density of leave-one-out cross validation with importance sampling

print(loo::loo(loo::extract_log_lik(stan_fit)), digits = 3)

# print log likelihood from RWiener package

print(RWiener::wiener_likelihood(c(a,t,b,v), rt))

# print means of posterior distributions for each parameter

print(unlist(lapply(rstan::extract(stan_fit), mean)))

# Plots -------------------------------------------------------------------

# rstan traceplot of parameters

rstan::traceplot(stan_fit, pars = c("a", "t", "b", "v"))

# function to plot distributions from RWiener

plot_rwiener <- function(rt) {

    rt_u <- rt$q[rt$resp == "upper"]

    rt_l <- rt$q[rt$resp == "lower"]

    xlim <- max(c(rt_u, rt_l))

    rt_seq <- seq(0, xlim, .005)

    d_rt_u <-

        RWiener::dwiener(rt_seq, a, t, b, v, resp = rep("upper", length(rt_seq)))

    d_rt_l <-

        RWiener::dwiener(rt_seq, a, t, b, v, resp = rep("lower", length(rt_seq)))

[Bob Carpenter]

Thanks much for following up with an example.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.