@JuliaRegistrator register()

I have been trying to use the WaldMixture in this SSM package. I am a long-time RT modeler, so I am familiar with both modeling and applications of these types of models. Here is the issue.

The "pdf" for the WaldMixture model does not actually produce a probability density function. That is, it does not produce a distribution that integrates to one (not even close for many parameters). For example, using parameters "v, sig, a, tnd = 3.0, 1.0, 2.0, 0.13;", the resulting density function integrates to ~1.31, and thus is not a pdf. The issue seems to be something with the drift variance since when sig=0.0, this problem goes away and it gets worse the larger sig becomes. I looked at the pdf code in src for this, and it seems to match the referenced paper, so not sure what is going on.

Here is some code that integrates and plots the appropriate pdf for different parameters using both the pdf built in here as well as a hand coded stochastic simulation pdf. Any thoughts? Am I doing something silly or is there an actual issue here?

using SequentialSamplingModels
using QuadGK
using Plots
using KernelDensity
using Distributions


#####
function Single_sim(v,a)

    dt = 0.001 # 1 ms
    sig_dt = 1.0 * sqrt(dt)

    Tfinal = 5.0;
    N_max = Int64(ceil(1.0/dt * Tfinal))

    x = 0.0
    thit = 1.0/eps(); # largest machine number
    for i = 1:N_max
        x = x + v*dt + randn() * sig_dt
        if x >= a
            thit = dt*i
            break
        end
    end

    return thit

end

#####
function WaldMixture_Nsim(v,σ,a,τ,Ns)

    vdist = Normal(v,σ);
    vvals = rand(vdist,Ns);

    rts = zeros(Ns)
    for i = 1:Ns
        rts[i] = Single_sim(vvals[i],a) + τ
    end

    return rts

end

v, sig, a, tnd = 3.0, 1.0, 2.0, 0.13; # Set parameters
dist = WaldMixture(v,sig,a,tnd); # Construct distribution
fun(x) = pdf(dist,x); # Generate pdf function to integrate
integral , error = quadgk(fun,tnd,5.0)  #Integrate

plotgrid = Vector(LinRange(tnd,5.0,1000)); # Set plot grid.
plotvals = pdf.(dist,plotgrid); # Evaluate WaldMixture on plot grid


#### Construct, integrate, and plot the pdf using stochastic simulation and density estimation.
Ns = 1_000_000; # Set number of simulations
rts = WaldMixture_Nsim(v,sig,a,tnd,Ns); # Simulate Ns wald accumulators
rts = rts[rts .< 1.0 / eps()]; # Remove any non-terminating rts
kde_rts = kde(rts,bandwidth = 0.0001); # Compute Kernel Density Estimate from data.
fun_kde(x) = pdf(kde_rts,x); # Construct a function to integrate and plot from kde.
integral , error = quadgk(fun_kde, tnd,5.0) # Integrate. 
plotvalskde = fun_kde.(plotgrid); # Evaluate kde distribution on plot grid.

# Plot "exact" and KDE approximation
plot(plotgrid,plotvals);
plot!(plotgrid,plotvalskde)

I'm having some trouble parametrizing the RDM, despite using priors that I'd think would be alright (but here's probably my mistake).
Very often, the sampling errors with the following:

ERROR: DomainError with -1.4187210312401369e-50:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
... [about log]
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math .\math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math .\special\log.jl:301
  [3] log
    @ .\special\log.jl:267 [inlined]
  [4] log
    @ C:\Users\domma\.julia\packages\ForwardDiff\PcZ48\src\dual.jl:240 [inlined]
  [5] logpdf
    @ C:\Users\domma\.julia\packages\SequentialSamplingModels\T1Z4T\src\RDM.jl:55 [inlined]
  [6] logpdf(d::RDM{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 4}}, r::Int64, rt::Float64)
    @ SequentialSamplingModels C:\Users\domma\.julia\packages\SequentialSamplingModels\T1Z4T\src\RDM.jl:144

Code:

using CSV
using DataFrames
using Turing
using Pigeons
using SequentialSamplingModels
using StatsModels
using StatsPlots
using GLMakie
using RCall

# Read data
df = CSV.read(download("https://raw.githubusercontent.com/RealityBending/DoggoNogo/main/study1/data/data_game.csv"), DataFrame)

# RDM model
@model function model_rdm(data; min_rt=minimum(data.rt), isi=nothing)

    # Priors for coefficients
    drift_intercept ~ filldist(truncated(Normal(6, 1), 0.0, Inf), 1)

    A ~ truncated(Normal(1, 0.3), 0.0, Inf)
    k ~ truncated(Normal(0.5, 0.1), 0.0, Inf)
    τ ~ truncated(Normal(0.2, 0.05), 0.0, min_rt)

    for i in 1:length(data)
        drift = drift_intercept
        data[i] ~ RDM(drift, k, A, τ)
    end
end
dat = [(choice=1, rt=df.RT[i]) for i in 1:length(df.RT)]
chain_rdm = sample(model_rdm(dat, min_rt=minimum(df.RT), isi=df.ISI), NUTS(), 500)
StatsPlots.plot(chain_rdm; size=(600, 2000))

Any tips?

Which makes me think that it would be amazing to add the docs, when listing all parameters, some info (where available, it can be a work-in-progress goal) to help people parametrize, e.g.: for sigma, "Cannot be negative", for drift in DDM: "typical range: [..., ....]"

@DominiqueMakowski, and @kiante-fernandez, for your situation awareness, predict is now working. Here is an example: https://itsdfish.github.io/SequentialSamplingModels.jl/dev/turing/

Code to test extensions. Will it actually work?

cd(@__DIR__)
using Pkg 
Pkg.activate("")
#using Revise
using SequentialSamplingModels
using Turing 

@model function wald_model(rts)
    ν ~ truncated(Normal(1.5, 1), 0, Inf)
    α ~ truncated(Normal(1, 1), 0, Inf)
    τ = 0.3
    rts ~ Wald(ν, α, τ)
    return (;ν, α, τ)
end

# generate simulated data
n_samples = 20
rts = rand(Wald(ν=1.5, α=.8, τ=.3), n_samples)
# make your model object 
model = wald_model(rts)
# estimate the parameters 
chain = sample(model, NUTS(1000, .85), MCMCThreads(), 1000, 1)
# make a submodel for generating the posterior predictive distribuiton
post_model = predict_posterior(Wald; model, func=mean, n_samples=20)
# generate posterior predictive distribution
post_pred_samples = generated_quantities(post_model, chain)

@JuliaRegistrator register() SequentialSamplingModels v0.1.4

@JuliaRegistrator register()

@JuliaRegistrator register()

@JuliaRegistrator register() SequentialSamplingModels v0.1.3

Few questions:

• Similar to LBA, can we enable dual SDs in LNR?
• Minor: Consider renaming the https://github.com/itsdfish/SequentialSamplingModels.jl/blob/master/src/LogNormalRace.jl file by LNR.jl for consistency
• In general, how much do you think there is room for optimizing (i.e., making faster) the code for the models (for instance for the logpdf calculation)? Would that impact sampling time or would that be marginal compared to the impact of AD?

Hey @itsdfish. @DominiqueMakowski and I noticed that the codebase creating custom distributions suitable for use with Turing.jl may be missing some methods. This omission can result in errors when users attempt to estimate parameters using Turing. see this thread.

To address this issue, I suggest adding the following methods to each distribution, following the guidelines outlined in the Distributions.jl documentation.

1. minimum(d::UnivariateDistribution)
2. maximum(d::UnivariateDistribution)
3. insupport(d::UnivariateDistribution, x::Real)
4. params(d::UnivariateDistribution)

While adding these methods, we can improve the usability, functionality, and standardization of using sequential sampling model distributions for use with Turing.

I am able to get started on this and opening a pull request with the necessary changes. Please let me know if there is any way I can assist or contribute further.

Best,

Kianté

I am encountering some problems with predict() despite using some code that previously worked (?), do you manage to reproduce the error:

using Turing
using SequentialSamplingModels
using Random
using LinearAlgebra
using StatsPlots
using Random
using Pigeons

Random.seed!(45461)


# Generate some data with known parameters
dist = LBA(ν=[3.0, 2.0], A=0.8, k=0.2, τ=0.3)
data = rand(dist, 100)

# Specify LBA model
@model function model_lba(data; min_rt=0.2)
    # Priors
    ν ~ MvNormal(zeros(2), I * 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    τ ~ Uniform(0.0, min_rt)

    # Likelihood
    data ~ LBA(; ν, A, k, τ)
end

chain = sample(model_lba(data; min_rt=minimum(data.rt)), NUTS(), 1000)

dat0 = [(missing) for i in 1:length(data.rt)]
pred = predict(model_lba(dat0; min_rt=minimum(data.rt)), chain)

2nd issue, with different error:

using DataFrames
using Random
using SequentialSamplingModels
using StatsModels
using Turing

# Generate data with different drifts for two conditions A vs. B
function make_data(; difference=1, n_groups=5, n_obs=100)
    n_obs_pergroup = n_obs ÷ n_groups

    # Create the first group (baseline / intercept)
    drift = [1.5, 0.5]  # Baseline
    df = DataFrame(rand(LBA(ν=drift, A=0.5, k=0.2, τ=0.3), n_obs_pergroup))

    # Compute the parameter change for each group
    change_drift = [
        difference / (n_groups - 1),
        difference / (n_groups - 1) / 2
    ]

    # Add new groups
    for g in 2:n_groups
        drift += change_drift  # new drift
        df = vcat(df, DataFrame(rand(LBA(ν=drift, A=0.5, k=0.2, τ=0.3), n_obs_pergroup)))
    end

    # Add condition column (if less than 5 groups, use letters for categorical groups)
    # Otherwise, assume continuous
    if n_groups <= 4
        df.x = repeat(["A", "B", "C", "D"][1:n_groups], inner=n_obs_pergroup)
    else
        df.x = repeat(range(0, 1, length=n_groups), inner=n_obs_pergroup)
    end

    return df
end



# Define models
@model function model_lba(data; min_rt=0.2, x=nothing)
    # Priors for auxiliary parameters
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    tau ~ Uniform(0.0, min_rt)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)

    # Priors for coefficients
    drift_intercept ~ filldist(Normal(0, 1), 2)
    drift_x ~ filldist(Normal(0, 1), 2)

    for i in 1:length(data)
        drifts = drift_intercept .+ drift_x * x[i]
        data[i] ~ LBA(; ν=drifts, τ=tau, A=A, k=k)
    end
end


# Generate data
Random.seed!(6)
df = make_data(difference=1, n_groups=2, n_obs=100)

# Format input data
f = @formula(rt ~ 1 + x)
f = apply_schema(f, schema(f, df))
_, predictors = coefnames(f)
X = modelmatrix(f, df)


# Sample model
data = [(choice=df.choice[i], rt=df.rt[i]) for i in 1:nrow(df)]  # Format data
chain = sample(model_lba(data; min_rt=minimum(df.rt), x=X[:, 2]), NUTS(), 1000)

# Predict
dat0 = [(missing) for i in 1:nrow(df)]
pred = predict(model_lba(dat0; min_rt=minimum(df.rt), x=X[:, 2]), chain)

@JuliaRegistrator register()

@JuliaRegistrator register() SequentialSamplingModels v0.1.1

@JuliaRegistrator register()

The team behind the paper showcasing the usage of LANs to SSMs is actively working on a Python package (https://github.com/lnccbrown/HSSM) that implements this "fast" approach to compute SSM parameters.

And so I was curious about what do you think of this approach? Is this something that could be implemented in Julia? Are other alternatives interesting as well? Is it going to be less relevant once Turing works with efficient ADs like Enzyme?

I feel like this question lies somewhere between SequentialSamplingModels and Turing itself, but let's assume the following basic LBA model:

using Turing
using SequentialSamplingModels
using Random
using LinearAlgebra
using DataFrames
using StatsPlots

# Generate data with different drifts for two conditions A vs. B
Random.seed!(254)
data = DataFrame(rand(LBA(ν=[3.0, 2.0], A=0.8, k=0.2, τ=0.3), 500))



@model function model_lba(data)
    data = (choice=data.choice, rt=data.rt)
    min_rt = minimum(data[2])

    # Priors
    drift ~ filldist(MvNormal(zeros(2), I * 2), 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    tau ~ Uniform(0.0, min_rt)

    # Likelihood
    data ~ LBA(; τ=tau, A=A, k=k, ν=drift)
    return (; data, drift, A, k, tau)
end


chain = sample(model_lba(data), Prior(), 50000)

(note that it is made to run with a dataframe input for convenience)

I can easily sample from priors. What I like to do next is to make generate predictions from these priors (to visualize a prior predictive check).

My first instinct was to simply run:

predict(model_lba(data), chain)

But it returns nothing - which I thought was a bug. The Turing team clarified the solution: it was to set the outcome variable to missing, which is straightforward in a linear model, but not here.

Attempt 1: Re-instantiating the model on a data with missing for choice:

predict(model_lba(DataFrames.transform(data, choice -> (choice=missing,))), chain)

ERROR: TypeError: non-boolean (Missing) used in boolean context
Stacktrace:
  [1] pdf
    @ C:\Users\domma\.julia\packages\SequentialSamplingModels\hMsCP\src\LBA.jl:108 [inlined]
  ...

Attempt 2: Re-write the model to have choice and rt as separate inputs, which would also make it more flexible (for instance to add predictor variables in the future).

@model function model_lba(choice, rt)
    min_rt = minimum(rt)

    # Priors
    drift ~ filldist(MvNormal(zeros(2), I * 2), 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    tau ~ Uniform(0.0, min_rt)

    # Likelihood
    (choice, rt) ~ LBA(; τ=tau, A=A, k=k, ν=drift)
    return (; choice, rt, drift, A, k, tau)
end

ERROR: LoadError: Malformed variable name (choice, rt)!

Turing doesn't seem to like the tuple output. So maybe we can workaround by creating the data tuple inside the model?

@model function model_lba(choice, rt)
    data = (choice=choice, rt=rt)
    min_rt = minimum(rt)

    # Priors
    drift ~ filldist(MvNormal(zeros(2), I * 2), 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    tau ~ Uniform(0.0, min_rt)

    # Likelihood
    data ~ LBA(; τ=tau, A=A, k=k, ν=drift)
    return (; data.choice, data.rt, drift, A, k, tau)
end

chain = sample(model_lba(data.choice, data.rt), Prior(), 50000)

ERROR: MethodError: no method matching iterate(::LBA{Matrix{Float64}, Float64, Float64, Float64})
  ...

Unfortunately no. I am not sure what else to try, and any pointers and thoughts are more than welcome ☺️

Side little question, does this package implements (or plans to add) tools for simple DDMs based on Wiener distributions? From my googling it doesn't seem to exist in Julia yet

Thanks again @itsdfish for your work on this, it's truly inspiring ☺️

Following up on a use-case I've mentioned somewhere else, in which once wants to estimate the drift rate(s) of an LBA as an intercept + a parameter. In a traditional linear regression, mu is set via the linear regression formula:

@model function linear_regression(x, y)
    σ² ~ truncated(Normal(0, 100); lower=0)
    intercept ~ Normal(0, sqrt(3))
    nfeatures = size(x, 2)
    coefficients ~ MvNormal(Zeros(nfeatures), 10.0 * I)

    # Linear formula here
    mu = intercept .+ x * coefficients
    return y ~ MvNormal(mu, σ² * I)
end

And I am wondering how to apply that to an LBA where the drift rate is a tuple.

Generate data

using Turing
using SequentialSamplingModels
using Random
using LinearAlgebra
using Distributions
using DataFrames
using StatsPlots
using StatsModels


# Generate data with different drifts for two conditions A vs. B
Random.seed!(254)
data1 = DataFrame(rand(LBA(ν=[3.0, 2.0], A=0.8, k=0.2, τ=0.3), 500))
data1[!, :condition] = repeat(["A"], nrow(data1))
data2 = DataFrame(rand(LBA(ν=[1.0, 1.5], A=0.7, k=0.2, τ=0.3), 500))
data2[!, :condition] = repeat(["B"], nrow(data2))

data = vcat(data1, data2)


# Format input data
f = @formula(rt ~ 1 + condition)
f = apply_schema(f, schema(f, data))

_, predictors = coefnames(f)
X = modelmatrix(f, data)

Here I am using the formula/modelmatrix API from StatsModels just to obtain the predictor data, which is simply a boolean matrix for the intercept and the condition:

1000×2 Matrix{Float64}:
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 1.0  0.0
 ⋮
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0

Here's my attempt at modelling this:

Model

@model function model_lba(data; intercept=nothing, condition=nothing, min_rt=minimum(data.rt))

    # Priors for auxiliary parameters
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    tau ~ Uniform(0.0, min_rt)

    # Priors over coefficients
    beta_intercept ~ Normal(0, 1)
    beta_condition ~ Normal(0, 0.5)

    # drift rate
    drift = beta_intercept + beta_condition * condition
    # drift ~ filldist(MvNormal(zeros(2), I * 2), 2)

    # Likelihood
    data ~ LBA(; τ=tau, A=A, k=k, ν=drift)
    return (; data, drift, A, k, tau)
end

Any suggestion are welcome

@JuliaRegistrator register()

We finally started piloting a project (a cognitive control game) that I am excited to analyze using SSMs. We gathered a few participants for the first "level" of the game, which is a simple reaction time task where people just have to answer as fast as they can when something appears (no choice).

Out of curiosity, I fitted the 3 models to the whole data (not caring for now about the presence of multiple participants) that I think are compatible with 1 choice data; Wald, LNR and LBA.

However, the LNR model gives me way-off predictions, and I'm not sure if I did something wrong or what's going on (see PPcheck plot). Do you see anything obvious that I might have overlooked?

Here's the reproducible code:

using CSV
using DataFrames
using Turing
using SequentialSamplingModels
using StatsPlots
using Pigeons
using GLMakie

# Read data
df = CSV.read(download("https://raw.githubusercontent.com/RealityBending/DoggoNogo/main/study1/data/data_game.csv"), DataFrame)

# Wald model
@model function model_wald(data; min_rt=minimum(data.rt))
    # Priors
    ν ~ truncated(Normal(5, 2), 0.0, Inf)
    α ~ truncated(Normal(.5, .4), 0.0, Inf)
    τ  ~ truncated(Normal(.2, .05), 0.0, min_rt)

    for i in 1:length(data)
        data[i] ~ Wald(; ν, α, τ)
    end
end

chain_wald = sample(model_wald(df.RT, min_rt=minimum(df.RT)), NUTS(), 1000)
# StatsPlots.plot(chain_wald; size=(600, 2000))
# summarystats(chain_wald)


# LBA model
@model function model_lba(data; min_rt=minimum(data.rt))
    ν ~ filldist(truncated(Normal(3, 5), 0.0, Inf), 1)
    σ ~ filldist(truncated(Normal(0, 1), 0.0, Inf), 1)
    A ~ truncated(Normal(0.4, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    τ ~ truncated(Normal(.2, .05), 0.0, min_rt)

    for i in 1:length(data)
        data[i] ~ LBA(; ν, σ, A, k, τ)
    end
end

dat = [(choice=1, rt=df.RT[i]) for i in 1:length(df.RT)]
chain_lba = sample(model_lba(dat, min_rt=minimum(df.RT)), NUTS(), 1000)


# LNR model
@model function model_lnr(data; min_rt=minimum(data.rt))
    ν ~ filldist(truncated(Normal(3, 5), 0.0, Inf), 1)
    σ ~ filldist(truncated(Normal(0, 3), 0.0, Inf), 1)
    τ ~ truncated(Normal(.2, .05), 0.0, min_rt)

    for i in 1:length(data)
        data[i] ~ LNR(; ν, σ, τ)
    end
end

chain_lnr = sample(model_lnr(dat, min_rt=minimum(df.RT)), NUTS(), 1000)

# PP Check
pred = predict(model_wald([missing for i in 1:100]; min_rt=minimum(df.RT)), chain_wald)
pred_rt_wald = Array(pred)[:, 2:2:end]

pred = predict(model_lba([(missing) for i in 1:100]; min_rt=minimum(df.RT)), chain_lba)
pred_rt_lba = Array(pred)[:, 2:2:end]

pred = predict(model_lnr([(missing) for i in 1:100]; min_rt=minimum(df.RT)), chain_lba)
pred_rt_lnr = Array(pred)[:, 2:2:end]

# Plot density
f = Figure()
ax1 = Axis(f[1, 1], title = "Wald")
for i in 1:size(pred_rt_wald)[2]
    lines!(ax1, Makie.KernelDensity.kde(pred_rt_wald[:, i]), color="green", alpha=0.05)
end
ax2 = Axis(f[2, 1], title = "LBA")
for i in 1:size(pred_rt_lba)[2]
    lines!(ax2, Makie.KernelDensity.kde(pred_rt_lba[:, i]), color="blue", alpha=0.05)
end
ax3 = Axis(f[3, 1], title = "LNR")
for i in 1:size(pred_rt_lnr)[2]
    lines!(ax3, Makie.KernelDensity.kde(pred_rt_lnr[:, i]), color="red", alpha=0.05)
end
GLMakie.xlims!(ax3, (0, 50))
for ax in [ax1, ax2]
    lines!(ax, Makie.KernelDensity.kde(df.RT), color="black")
    GLMakie.xlims!(ax, (0.15, 0.8))
end
f

I have been thinking about more general use cases, and one common thing I like to do is simulate data at the trial level or for multiple subjects. As an enhancement, I suggest we showcase these examples, perhaps under a new tab in the documentation titled "Simulating Data." Here is a quick working example of the trial level idea.

using SequentialSamplingModels
using Plots
using Random
using Distributions 
Random.seed!(2024)

#simulating trial level effects via regressions models

# Set up trial by trial parameters
intercepts = [0.3, 0.4, 0.5, 0.6]
x = rand(Uniform(1, 10), 1000) #generate sets of items values

v = (intercepts' .+ β .* x) #trail level drifts

A = 1.5
k = 0.5
τ = 0.0

Θ = hcat(v, repeat([A, k, τ]', 1000, 1))

num_rows, num_cols = size(Θ)

choices = Vector{Int64}(undef, num_rows)
rts = Vector{Float64}(undef, num_rows)

for i in 1:num_rows
    ν = Θ[i, 1:4]
    k = Θ[i, 5]
    A = Θ[i, 6]
    τ = Θ[i, 7]

    dist = RDM(;ν, k, A, τ)
    choice, rt = rand(dist, 1)

    choices[i] = choice[1]
    rts[i] = rt[1]
end

I can write this out along with a hierarchical example.

Minor: the σ parameters does not seem to be documented: https://itsdfish.github.io/SequentialSamplingModels.jl/dev/lba/ in the constructor, is this expected?

Also, side question, is it possible to fit an LBA with only one accumulator? Say to compare with a Wald model?

Hi @kiante-fernandez. I am moving a discussion of the PR here. Unfortunately, I couldn't figure out how to view the PR locally before merging it (it's possible that it would have worked, but I was in the wrong folder).

I would like to bring up some issues. First, the constructors were not working. I ended up commenting out the following to get it to work:

# function DDM(ν::T, α::T, τ::T, z::T; check_args::Bool=true) where {T <: Real}
#     check_args && Distributions.@check_args DDM (α, α > zero(α)) (τ, τ > zero(τ)) (z, z ≥ 0 && z ≤ 1)
#     return DDM{T}(ν, α, τ, z)
# end

# function DDM(ν::T, α::T, τ::T; check_args::Bool=true) where {T <: Real}
#     return DDM(ν, α, τ, 0.5; check_args=check_args)
# end

# DDM(ν::Real, α::Real, τ::Real, z::Real; check_args::Bool=true) = DDM(promote(ν, α, τ, z)...; check_args=check_args)
# DDM(ν::Real, α::Real, τ::Real; check_args::Bool=true) = DDM(promote(ν, α, τ)...; check_args=check_args)

I don't think you will be able to use DDM with a keyword argument because Julia unfortunately does not dispatch keyword arguments. A different approach will be needed if you would like to include checks on the correctness of the arguments. It might not be possible to opt in/out.

I'm not much of a fan about the negative RT scheme. It might be a clever way to make the model faster but I think it has two drawbacks worth considering. First, it has the potential to lead to confusion or errors. What if a person computes a mean over rts, or some other operation unaware of the presence of negative rts? Perhaps it might be easier to code correctly this way, but we only have to do that once. By contrast, users will have to use it appropriately everytime. A second drawback is that it is an exception to the general interface, which is to return an index vector of choices, and a vector of RTs. I think having a standardized interface is important because it sets up expectations for how the package works, and allows the user to potentially use the same code without having to modify depending on whether it uses LBA or DDM.

Another issue I noticed was that pdf returns NaN for negative RTs. I'm not sure why it does that. As a user, I would expect that I could pass simulated data to pdf without a problem. A similar problem exists for logpdf

Example

Random.seed!(5484)
dist = DDM(.5, .8, .3, .4)

rts = rand(dist, 10)

logpdf.(dist, rts)

Output

10-element Vector{Float64}:
 NaN
 NaN
  -2.0485034773320545
 NaN
 NaN
 NaN
  -2.6704531377602465
  -3.9876172218519876
 NaN
   0.34653192776453606

My guess is that there is a bug due to method ambiguity:

pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T<:Real}
pdf(d::DDM, t::Real; ϵ::Real = 1.0e-12)

I don't believe these methods can be distinguished based on their signatures, which causes the last one specified to be used. The last method has

    if τ ≥ t
        return T(NaN)
    end

which presumably is producing the NaN above.

But those are small details which should be easy to fix.

I am opposed to breaking the interface, but I wonder whether a solution would be to transform the data coming in and out? For example, in the pdf function, you might have:

rts = copy(_rts_in)
rts[choice .== 1] .= -rts[choice .==1]

Its less than ideal because it is copying data, but I would be amenable to that solution because it preserves the interface. If that is not a viable solution for you, I think the best approach moving forward would be keep the packages separate. In my opinion, the benefit of having models within the same package is that they follow the same interface. Eventually, I or someone else will add the full DDM (where the negative RT scheme possibly does not apply).

I am finishing up the DDM with across-trial variability in drift, starting point, and non-decision time.

For drift, you can use analytic integration of the likelihood function. See HDDM for a use case example. That works fine.

For the starting point and non-decision time, numerical integration is used (Ratcliff and Tuerlinckx, 2002). This part tends to be the most computationally expensive.

I am going to start with just coding up Simpson's Method. But I figured I would ask if you knew of some alternatives that are notoriously fast in Julia. I was looking at something like HCubature.jl or QuadGK.jl.

In this example, Bijectors throws an error when the data are transformed inside the main loop of the Turing model. The last example in #31 shows a workaround.

I am going through Myers (2022) and on the model goodness-of-fit section.

They mention different methods to compare models together, and I am to figure out how to do that using SSM & Turing. Apologies if that's slightly out of your depth, but I thought I'd start a discussion to help users achieve this critical step of their analysis.

Maximum Log Likelihood (LLE)

One of the most basic methods seems to be the max LLE.

For a target distribution

If I'm not mistaken, computing the max LLE between "observed" data and a distribution with known parameter is showcased in the example in rtdists:

objective_fun <- function(par, rt, response, distribution = "norm") {
  ...
  # get summed log-likelihood:
  d <- do.call(dLBA, args = c(rt=list(rt), response=list(response), spar, dist_par, 
                               distribution=distribution, silent=TRUE))
  if (any(d == 0)) return(1e6)
  else return(-sum(log(d)))
}

Am I correct in assuming that the equivalent can be obtained in SSM as follows:

m = LBA(ν=[2.4, 1.6], A=0.5, k=0.5, τ=0.5)
choice, rt = rand(m, 500)

-sum(logpdf.(m, choice, rt))

For a Turing model

A Turing chain has 2 (internal) parameters; lp and log_density

DataFrame(chain)[!, :log_density]

Which are the same. Do these correspond to the above? And do their mean/median would give us a reasonable goodness-of-fit index?

BIC / DIC

BIC (BIC = k*ln(n) − 2*LLE) is mentioned as a better index than the raw LLE. It is trivial to compute once we have LLE, but Turing has a method to compute DIC.

Thei method requires however to pass a function, that in the example they define as:

lpfun = function f(chain::Chains) # function to compute the logpdf values
    niter, nparams, nchains = size(chain)
    lp = zeros(niter + nchains) # resulting logpdf values
    for i = 1:nparams
        lp += map(p -> logpdf( ... , x), Array(chain[:,i,:]))
    end
    return lp
end

DIC, pD = dic(chain, lpfun)

That said I am not sure what to input in the ... of logpdf 🤔

Any thoughts or pointers are welcome!

It would be nice to have plot() methods for the various distributions to be able to quickly plot the density distributions, something like PyDDM has:

@JuliaRegistrator register()

The LBA docs lists sigma before tau in the list but the constructor expects sigma as last
(it's minor but it just took me a few minutes to debug that 😅)

@JuliaRegistrator register()

@kiante-fernandez, I noticed that cdf returns a large negative value for values near tau. Do you know if this is an expected edge case or perhaps a bug? I suspect it might be an edge case, but I don't know for certain.

using SequentialSamplingModels
dist = DDM(ν=2.0, α = 1.5, z = .5, τ = .30) 
cdf(dist, 1, 0.305)

Was just wondering about the possibility of implementing models with time-varying boundaries (or collapsing boundaries). That said, it might make sense to not have that as a priority, given some of their critiques.

• https://www.frontiersin.org/articles/10.3389/fpsyg.2014.01364/full (implementation)
• https://www.jneurosci.org/content/35/6/2476 (a critique of these models)
• https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5450920/ (another critique)

@JuliaRegistrator register()

@JuliaRegistrator register() SequentialSamplingModels v0.1.2

@JuliaRegistrator register()

Hi and thanks for this awesome package!

I was wondering if there was a roadmap for adding wrapper to conveniently fit models based on these distributions (using a GLM.jl-like interface) or a guide for writing Bayesian models with Turing.

Thanks a lot for any pointers :)

I would like to explore the possibility of defining a type for multialternative SSMs in hopes that this package plays well with Turing and other packages. The package currently works well in most cases. However, it does not work well with predict from Turing, as discussed here.

Here is one idea:

using Distributions
using Turing 
import Distributions: logpdf
import Distributions: loglikelihood
import Base: length

abstract type Mixed <: ValueSupport end 

const MixedMultivariateDistribution = Distribution{Multivariate, Mixed}

abstract type SSM1D <: ContinuousUnivariateDistribution end 

abstract type SSM2D <: MixedMultivariateDistribution end 

This defines MixedMultivariateDistribution which could potentially be used outside of SSMs. The type system is then split into 1D and 2D SSMs, which are abstract types.

The code below shows that this works for basic MCMC sampling. The question is how to get it to work with predict and friends.

# not really 2D, just for illustration
struct MyType{T<:Real} <: SSM2D
    n::Int
    x::T
end

logpdf(d::MyType,data::Int) = logpdf(Binomial(d.n, d.x), data)

loglikelihood(d::MyType,data::Int) = loglikelihood(Binomial(d.n, d.x), data)


@model function my_model(n, k)
    θ ~ Beta(1, 1)
    return k ~ MyType(n, θ)
end

chain = sample(my_model(10, 5), NUTS(), 3_000)

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!