From 2e962927ac6ae471acb452516d453ae86738ed8c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Mon, 26 Feb 2024 21:58:29 -0800
Subject: [PATCH 01/10] Add stDDM to SequentialSamplingModels

---
 src/SequentialSamplingModels.jl |   3 +
 src/stDDM.jl                    | 207 ++++++++++++++++++++++++++++++++
 src/type_system.jl              |   7 ++
 3 files changed, 217 insertions(+)
 create mode 100644 src/stDDM.jl

diff --git a/src/SequentialSamplingModels.jl b/src/SequentialSamplingModels.jl
index 050ec352..b814785e 100644
--- a/src/SequentialSamplingModels.jl
+++ b/src/SequentialSamplingModels.jl
@@ -36,6 +36,7 @@ module SequentialSamplingModels
     export AbstractLNR
     export AbstractPoissonRace
     export AbstractRDM 
+    export AbstractstDDM 
     export AbstractWald
     export aDDM
     export CDDM
@@ -49,6 +50,7 @@ module SequentialSamplingModels
     export PoissonRace
     export SSM1D
     export SSM2D
+    export stDDM
     export ContinuousMultivariateSSM
     export Wald
     export WaldMixture 
@@ -92,4 +94,5 @@ module SequentialSamplingModels
     include("ex_gaussian.jl")
     include("poisson_race.jl")
     include("RatcliffDDM.jl")
+    include("stDDM.jl")
 end
\ No newline at end of file
diff --git a/src/stDDM.jl b/src/stDDM.jl
new file mode 100644
index 00000000..2409ab75
--- /dev/null
+++ b/src/stDDM.jl
@@ -0,0 +1,207 @@
+"""
+    stDDM{T<:Real} <: AbstractstDDM
+
+An object for the starting-time diffusion decision model.
+
+# Parameters 
+
+- `ν`:  vector of drift rates for attribute one and two
+- `α`:  evidence threshold
+- `τ`:  non-decision time 
+- `s`:  initial latency bias (positive for attribute two, negative for attribute one)
+- `z`:  initial evidence 
+- `η`:  vector of variability in drift rate for attribute one and two
+- `σ`:  diffusion noise 
+- `Δt`: time step 
+
+# Constructors 
+
+    stDDM(ν, α, τ, s, z, η, σ, Δt)
+
+    stDDM(;ν = [0.5,0.6], 
+        α = 1.0, 
+        τ = .300, 
+        s = 0.50, 
+        z = 0.50, 
+        η = [1.0,1.0], 
+        σ = 1,
+        Δt = .001)
+        
+# Example 
+
+```julia 
+using SequentialSamplingModels
+
+ν = [0.5, 0.6]
+α = 1.0
+τ = 0.300
+s = 0.50
+z = 0.50
+η = [1.0, 1.0]
+σ = 1
+Δt = 0.001
+
+# Create stDDM model instance
+dist = stDDM(;ν, α, τ, s, z, η, σ, Δt)
+
+choices,rts = rand(dist, 500)
+```
+
+# References
+
+Amasino, D.R., Sullivan, N.J., Kranton, R.E. et al. Amount and time exert independent influences on intertemporal choice. Nat Hum Behav 3, 383–392 (2019). https://doi.org/10.1038/s41562-019-0537-2
+Barakchian, Z., Beharelle, A.R. & Hare, T.A. Healthy decisions in the cued-attribute food choice paradigm have high test-retest reliability. Sci Rep, (2021). https://doi.org/10.1038/s41598-021-91933-6
+Chen, HY., Lombardi, G., Li, SC. et al. Older adults process the probability of winning sooner but weigh it less during lottery decisions. Sci Rep, (2022). https://doi.org/10.1038/s41598-022-15432-y
+Lombardi, G., & Hare, T. Piecewise constant averaging methods allow for fast and accurate hierarchical Bayesian estimation of drift diffusion models with time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
+Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate of information accumulation. Nat Hum Behav 5, 1698–1706 (2021). https://doi.org/10.1038/s41562-021-01154-0
+"""
+
+mutable struct stDDM{T<:Real} <: AbstractstDDM
+    ν::Vector{T}
+    α::T
+    τ::T
+    s::T
+    z::T
+    η::Vector{T}
+    σ::T
+    Δt::T
+end
+
+function stDDM(ν, α, τ, s, z, η, σ, Δt)
+    _, α, τ, s, z,_ , σ, Δt = promote(ν[1], α, τ, s, z, η[1], σ, Δt)
+    ν = convert(Vector{typeof(τ)}, ν)
+    η = convert(Vector{typeof(τ)}, η)
+    return stDDM(ν, α, τ, s, z, η, σ, Δt)
+end
+
+function stDDM(;ν = [0.5,0.6], 
+    α = 1.0, 
+    τ = .300, 
+    s = 0.50, 
+    z = 0.50, 
+    η = fill(1.0, length(ν)), 
+    σ = 1,
+    Δt = .001)
+
+    return stDDM(ν, α, τ, s, z, η, σ, Δt)
+end
+
+function params(d::AbstractstDDM)
+    (d.ν, d.α, d.τ, d.s, d.z, d.η, d.σ, d.Δt)    
+end
+
+get_pdf_type(d::AbstractstDDM) = Approximate
+
+"""
+    rand(dist::AbstractstDDM)
+
+Generate a random choice-rt pair for starting-time diffusion decision model.
+
+# Arguments
+- `dist`: model object for the starting-time diffusion decision model. 
+"""
+function rand(rng::AbstractRNG, dist::AbstractstDDM)
+    return simulate_trial(rng, dist)
+end
+
+"""
+    rand(dist::AbstractstDDM, n_sim::Int)
+
+Generate `n_sim` random choice-rt pairs for the starting-time diffusion decision model.
+
+# Arguments
+- `dist`: model object for the starting-time diffusion decision model.
+- `n_sim::Int`: the number of simulated choice-rt pairs  
+"""
+function rand(rng::AbstractRNG, dist::AbstractstDDM, n_sim::Int)
+    choices = fill(0, n_sim)
+    rts = fill(0.0, n_sim)
+    for i in 1:n_sim
+        choices[i],rts[i] = simulate_trial(rng, dist)
+    end
+    return (;choices,rts) 
+end
+
+# noise(rng, σ) = rand(rng, Normal(0, σ))
+
+"""
+simulate_trial(rng::AbstractRNG, dist::AbstractstDDM, TMax)
+
+Generate a single simulated trial from the starting-time diffusion decision model.
+
+# Arguments
+
+- `rng`: a random number generator
+- `model::AbstractstDDM`: a starting-time diffusion decision model object
+- `TMax`: total/max time for simulation
+"""
+
+function simulate_trial(rng::AbstractRNG, dist; TMax=6)
+    (;ν, α, τ, s, z, η, σ, Δt) = dist
+
+    lt = Int(TMax / Δt)
+    vec_tν1 = ones(Int, lt)
+    vec_tν2 = ones(Int, lt)
+    aux = abs(Int(s / Δt))
+    
+    if s > 0
+        vec_tν1[1:aux] .= 0
+    elseif s < 0
+        vec_tν2[1:aux] .= 0
+    end
+    
+    t = TMax
+    choice = 0  # Initialize choice with a default value
+
+    X = z * α
+    flag = false
+    cont = 1
+    while !flag && cont <= lt
+        # noise = noise(rng, σ) * sqrt(Δt)
+        noise = rand(rng, Normal(0, σ)) * sqrt(Δt)
+
+        X += (ν[1] * η[1] * vec_tν1[cont] + ν[2] * η[2] * vec_tν2[cont]) * Δt + noise
+        
+        if X > α
+            # t = τ + cont * Δt
+            choice = 1
+            flag = true
+        elseif X < 0
+            # t = -τ - cont * Δt
+            choice = 2
+            flag = true
+        end
+        t = τ + cont * Δt
+
+        cont += 1
+    end
+    
+    return (;choice,rt=t)
+
+end
+
+
+# """
+#     simulate(model::AbstractstDDM; _...)
+
+# Returns a matrix containing evidence samples of the stDDM decision process. In the matrix, rows 
+# represent samples of evidence per time step and columns represent different accumulators.
+
+# # Arguments
+
+# - `model::AbstractstDDM`: a starting-time diffusion decision model diffusion model object
+# """
+# function simulate(rng::AbstractRNG, model::AbstractstDDM, _...)
+#     (;ν,α,z,Δt) = model
+#     x = α * z
+#     t = 0.0
+#     evidence = [x]
+#     time_steps = [t]
+#     while (x < α) && (x > 0)
+#         t += Δt
+#         x += ν * Δt + rand(rng, Normal(0.0, 1.0)) * √(Δt)
+#         push!(evidence, x)
+#         push!(time_steps, t)
+#     end
+#     return time_steps,evidence
+# end
diff --git a/src/type_system.jl b/src/type_system.jl
index f2fa1dc1..b0e7c374 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -65,6 +65,13 @@ An abstract type for the Poisson race model.
 """
 abstract type AbstractPoissonRace <:SSM2D end
 
+"""
+    AbstractstDDM <: SSM2D
+
+An abstract type for the starting-time diffusion decision model.
+"""
+abstract type AbstractstDDM <: SSM2D end
+
 """
     AbstractRDM <: SSM2D
 

From 95564b88363e67e95115522ce39db816a102b468 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sat, 9 Mar 2024 23:22:44 -0800
Subject: [PATCH 02/10] revisions and update stDDM model constructor to include
 correlation between drift rates

---
 src/stDDM.jl       | 158 +++++++++++++++++++++++----------------------
 src/type_system.jl |  20 +++---
 2 files changed, 92 insertions(+), 86 deletions(-)

diff --git a/src/stDDM.jl b/src/stDDM.jl
index 2409ab75..bd4fa3d3 100644
--- a/src/stDDM.jl
+++ b/src/stDDM.jl
@@ -11,12 +11,13 @@ An object for the starting-time diffusion decision model.
 - `s`:  initial latency bias (positive for attribute two, negative for attribute one)
 - `z`:  initial evidence 
 - `η`:  vector of variability in drift rate for attribute one and two
-- `σ`:  diffusion noise 
+- `ρ`:  correlation between drift rate for attributes
+- `σ`:  diffusion noise
 - `Δt`: time step 
 
 # Constructors 
 
-    stDDM(ν, α, τ, s, z, η, σ, Δt)
+    stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
 
     stDDM(;ν = [0.5,0.6], 
         α = 1.0, 
@@ -24,6 +25,7 @@ An object for the starting-time diffusion decision model.
         s = 0.50, 
         z = 0.50, 
         η = [1.0,1.0], 
+        ρ = 0.00,
         σ = 1,
         Δt = .001)
         
@@ -38,11 +40,12 @@ using SequentialSamplingModels
 s = 0.50
 z = 0.50
 η = [1.0, 1.0]
+ρ = 0.00
 σ = 1
 Δt = 0.001
 
 # Create stDDM model instance
-dist = stDDM(;ν, α, τ, s, z, η, σ, Δt)
+dist = stDDM(;ν, α, τ, s, z, η, ρ, σ, Δt)
 
 choices,rts = rand(dist, 500)
 ```
@@ -63,15 +66,16 @@ mutable struct stDDM{T<:Real} <: AbstractstDDM
     s::T
     z::T
     η::Vector{T}
+    ρ::T
     σ::T
     Δt::T
 end
 
-function stDDM(ν, α, τ, s, z, η, σ, Δt)
-    _, α, τ, s, z,_ , σ, Δt = promote(ν[1], α, τ, s, z, η[1], σ, Δt)
+function stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
+    _, α, τ, s, z,_ , ρ, σ, Δt = promote(ν[1], α, τ, s, z, η[1], ρ, σ, Δt)
     ν = convert(Vector{typeof(τ)}, ν)
     η = convert(Vector{typeof(τ)}, η)
-    return stDDM(ν, α, τ, s, z, η, σ, Δt)
+    return stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
 end
 
 function stDDM(;ν = [0.5,0.6], 
@@ -80,14 +84,15 @@ function stDDM(;ν = [0.5,0.6],
     s = 0.50, 
     z = 0.50, 
     η = fill(1.0, length(ν)), 
+    ρ = 0.0,
     σ = 1,
     Δt = .001)
 
-    return stDDM(ν, α, τ, s, z, η, σ, Δt)
+    return stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
 end
 
 function params(d::AbstractstDDM)
-    (d.ν, d.α, d.τ, d.s, d.z, d.η, d.σ, d.Δt)    
+    (d.ν, d.α, d.τ, d.s, d.z, d.η, d.ρ, d.σ, d.Δt)    
 end
 
 get_pdf_type(d::AbstractstDDM) = Approximate
@@ -104,26 +109,6 @@ function rand(rng::AbstractRNG, dist::AbstractstDDM)
     return simulate_trial(rng, dist)
 end
 
-"""
-    rand(dist::AbstractstDDM, n_sim::Int)
-
-Generate `n_sim` random choice-rt pairs for the starting-time diffusion decision model.
-
-# Arguments
-- `dist`: model object for the starting-time diffusion decision model.
-- `n_sim::Int`: the number of simulated choice-rt pairs  
-"""
-function rand(rng::AbstractRNG, dist::AbstractstDDM, n_sim::Int)
-    choices = fill(0, n_sim)
-    rts = fill(0.0, n_sim)
-    for i in 1:n_sim
-        choices[i],rts[i] = simulate_trial(rng, dist)
-    end
-    return (;choices,rts) 
-end
-
-# noise(rng, σ) = rand(rng, Normal(0, σ))
-
 """
 simulate_trial(rng::AbstractRNG, dist::AbstractstDDM, TMax)
 
@@ -133,45 +118,46 @@ Generate a single simulated trial from the starting-time diffusion decision mode
 
 - `rng`: a random number generator
 - `model::AbstractstDDM`: a starting-time diffusion decision model object
-- `TMax`: total/max time for simulation
+- `max_steps`: total/max time for simulation
 """
 
-function simulate_trial(rng::AbstractRNG, dist; TMax=6)
-    (;ν, α, τ, s, z, η, σ, Δt) = dist
-
-    lt = Int(TMax / Δt)
-    vec_tν1 = ones(Int, lt)
-    vec_tν2 = ones(Int, lt)
-    aux = abs(Int(s / Δt))
+function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; max_steps=6)
+    (;ν, α, τ, s, z, η, ρ, σ, Δt) = d
     
-    if s > 0
-        vec_tν1[1:aux] .= 0
-    elseif s < 0
-        vec_tν2[1:aux] .= 0
-    end
+    lt = Int(max_steps / Δt)
+    start_step = abs(Int(s / Δt))
     
-    t = TMax
+    t = τ
     choice = 0  # Initialize choice with a default value
 
     X = z * α
-    flag = false
+    deciding = true
     cont = 1
-    while !flag && cont <= lt
-        # noise = noise(rng, σ) * sqrt(Δt)
-        noise = rand(rng, Normal(0, σ)) * sqrt(Δt)
 
-        X += (ν[1] * η[1] * vec_tν1[cont] + ν[2] * η[2] * vec_tν2[cont]) * Δt + noise
+    Ρ = [1.0 ρ; ρ 1.0]    
+    # Convert the standard deviations to a diagonal matrix
+    Σ_diag = [η[1] 0; 0 η[2]]
+    # Calculate the covariance matrix from correlation
+    Σ = Σ_diag * Ρ * Σ_diag
+
+    evidence = rand(MvNormal([ν[1], ν[2]], Σ))
+
+    while deciding && cont <= lt
+        δ1 = cont ≤ start_step && s > 0 ? 0.0 : 1.0  
+        δ2 = cont ≤ start_step && s < 0 ? 0.0 : 1.0 
+
+        noise = rand(rng, Normal(0, σ)) * √(Δt)
+
+        X += (evidence[1] * δ1 + evidence[2] * δ2) * Δt + noise
         
         if X > α
-            # t = τ + cont * Δt
             choice = 1
-            flag = true
+            deciding = false
         elseif X < 0
-            # t = -τ - cont * Δt
             choice = 2
-            flag = true
+            deciding = false
         end
-        t = τ + cont * Δt
+        t += Δt
 
         cont += 1
     end
@@ -181,27 +167,47 @@ function simulate_trial(rng::AbstractRNG, dist; TMax=6)
 end
 
 
-# """
-#     simulate(model::AbstractstDDM; _...)
-
-# Returns a matrix containing evidence samples of the stDDM decision process. In the matrix, rows 
-# represent samples of evidence per time step and columns represent different accumulators.
-
-# # Arguments
-
-# - `model::AbstractstDDM`: a starting-time diffusion decision model diffusion model object
-# """
-# function simulate(rng::AbstractRNG, model::AbstractstDDM, _...)
-#     (;ν,α,z,Δt) = model
-#     x = α * z
-#     t = 0.0
-#     evidence = [x]
-#     time_steps = [t]
-#     while (x < α) && (x > 0)
-#         t += Δt
-#         x += ν * Δt + rand(rng, Normal(0.0, 1.0)) * √(Δt)
-#         push!(evidence, x)
-#         push!(time_steps, t)
-#     end
-#     return time_steps,evidence
-# end
+"""
+    simulate(model::AbstractstDDM)
+
+Returns a matrix containing evidence samples of the stDDM decision process. In the matrix, rows 
+represent samples of evidence per time step and columns represent different accumulators.
+
+# Arguments
+
+- `model::AbstractstDDM`: a starting-time diffusion decision model diffusion model object
+"""
+function simulate(rng::AbstractRNG, model::AbstractstDDM)
+    (;ν, α, s, z, η, ρ, σ, Δt) = model
+
+    x = α * z
+    t = 0.0
+    evidence = [x]
+    time_steps = [t]
+    cont = 1
+    start_step = abs(Int(s / Δt))
+
+    Ρ = [1.0 ρ; ρ 1.0]    
+    # Convert the standard deviations to a diagonal matrix
+    Σ_diag = [η[1] 0; 0 η[2]]
+    # Calculate the covariance matrix from correlation
+    Σ = Σ_diag * Ρ * Σ_diag
+    
+    while (x < α) && (x > 0)
+        t += Δt
+
+        δ1 = cont ≤ start_step && s > 0 ? 0.0 : 1.0  
+        δ2 = cont ≤ start_step && s < 0 ? 0.0 : 1.0 
+
+        increment = rand(rng, MvNormal([ν[1], ν[2]], Σ))
+        noise = rand(rng, Normal(0, σ)) * √(Δt)
+
+        x += (increment[1] * δ1 + increment[2] * δ2) * Δt + noise
+
+        push!(evidence, x)
+        push!(time_steps, t)
+        cont += 1
+    end
+
+    return time_steps,evidence
+end
\ No newline at end of file
diff --git a/src/type_system.jl b/src/type_system.jl
index b0e7c374..e7e867aa 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -113,7 +113,7 @@ Base.broadcastable(x::ContinuousMultivariateSSM) = Ref(x)
 
 Base.length(d::SSM2D) = 2
 
-rand(d::SSM2D) = rand(Random.default_rng(), d)
+rand(d::SSM2D; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM, n::Int; kwargs...) = rand(Random.default_rng(), d, n; kwargs...)
 
@@ -166,15 +166,15 @@ with more than one choice option.
 
 # Arguments
 - `d::SSM2D`: a 2D sequential sampling model.
-- `n_sim::Int`: the number of simulated choices and rts  
-"""
-function rand(rng::AbstractRNG, d::SSM2D, N::Int)
-    choice = fill(0, N)
-    rt = fill(0.0, N)
-    for i in 1:N
-        choice[i],rt[i] = rand(rng, d)
-    end
-    return (;choice,rt)
+- `N::Int`: the number of simulated choices and rts  
+"""
+function rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)  
+    choice = fill(0, N)  
+    rt = fill(0.0, N)  
+    for i in 1:N  
+        choice[i],rt[i] = rand(rng, d; kwargs...)  
+    end  
+    return (;choice,rt)  
 end
 
 """

From 2c9f1052bb0cf81030b5fc5ea0c1f0798b8fa616 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sun, 10 Mar 2024 16:30:38 -0700
Subject: [PATCH 03/10] Add StatsBase package and import cor2cov function added
 markdown for stDDM

---
 Project.toml                    |   1 +
 docs/src/stDDM.md               | 106 ++++++++++++++++++++++++++++++++
 src/SequentialSamplingModels.jl |   1 +
 src/stDDM.jl                    |  12 ++--
 4 files changed, 114 insertions(+), 6 deletions(-)
 create mode 100644 docs/src/stDDM.md

diff --git a/Project.toml b/Project.toml
index 109eb5fe..78ec5c09 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/stDDM.md b/docs/src/stDDM.md
new file mode 100644
index 00000000..3a7bee9c
--- /dev/null
+++ b/docs/src/stDDM.md
@@ -0,0 +1,106 @@
+# Starting-time Diffusion Decision Model (stDDM)
+
+The relative starting time drift diffusion model (stDDM) characterizes the contributions of multiple unique attributes to the rate of evidence accumulation. Compared to the DDM, which assumes a constant evidence accumulation rate within each trial, the stDDM allows different attributes to enter the evidence accumulation process at various time points relative to one another. By doing so, the stDDM quantifies both the weights given to each attribute and their onset times (Amasino et al., 2019; Barakchian et al., 2021; Chen et al., 2022; Maier et al., 2020; Sullivan and Huettel, 2021).
+
+# Example
+In this example, we will demonstrate how to use the stDDM in a generic two-alternative forced-choice task with two arbitrary attributes.
+
+## Load Packages
+The first step is to load the required packages.
+
+```@example stDDM
+using SequentialSamplingModels
+using Plots
+using Random
+
+Random.seed!(8741)
+```
+
+## Create Model Object
+In the code below, we will define parameters for the stDDM and create a model object to store the parameter values. 
+
+### Drift Rates
+The drift rates control the speed and direction with which information accumulates, with one drift rate per attribute (e.g., taste and health, payoff and delay, self and other).
+```@example stDDM
+ν = [2.5,2.0]
+```
+
+### Threshold
+The threshold α represents the amount of evidence required to make a decision.
+```@example stDDM 
+α = 1.5
+```
+
+### Non-Decision Time
+Non-decision time is an additive constant representing encoding and motor response time. 
+```@example stDDM 
+τ = 0.30
+```
+
+### starting time
+The starting time parameter \(s\) denotes how much earlier one attribute begins to affect the evidence accumulation process relative to the other(s). If \(s\) is negative, attribute 1 evidence is accumulated before attribute 2 evidence; if \(s\) is positive, attribute 1 evidence is accumulated after attribute 2 evidence. The absolute value of \(s\) indicates the difference in starting times for the two attributes.
+```@example stDDM 
+s = 0.10 
+```
+
+### Starting Point
+An indicator of an an initial bias towards a decision. The z parameter is relative to a (i.e. it ranges from 0 to 1).
+```@example stDDM 
+z = 0.50
+```
+
+### Drift Rates Dispersion
+Dispersion parameters of the drift rate are drawn from a multivariate normal distribution, with the mean vector ν describing the distribution of actual drift rates from specific trials. The standard deviation or across-trial variability is captured by the η vector, and the corresponding correlation between the two attributes is denoted by ρ.
+```@example stDDM
+η = [1.0,1.0]
+ρ = 0.3
+```
+
+### Diffusion Noise
+Diffusion noise is the amount of within trial noise in the evidence accumulation process. 
+```@example stDDM 
+σ = 1.0
+```
+### Time Step
+The time step parameter $\Delta t$ is the precision of the discrete time approxmation. 
+```@example stDDM 
+Δt = .001
+```
+
+### stDDM Constructor 
+Now that values have been asigned to the parameters, we will pass them to `stDDM` to generate the model object.
+```@example stDDM 
+dist = stDDM(;ν, α, τ, s, z, η, ρ, σ, Δt)
+```
+
+## Simulate Model
+Now that the model is defined, we will generate $10,000$ choices and reaction times using `rand`. 
+ ```@example stDDM 
+ choices,rts = rand(dist, 10_000)
+```
+
+## Compute Choice Probability
+The choice probability $\Pr(C=c)$ is computed by passing the model and choice index to `cdf`.
+ ```@example stDDM 
+cdf(dist, 1)
+```
+To compute the joint probability of choosing $c$ within $t$ seconds, i.e., $\Pr(T \leq t \wedge C=c)$, pass a third argument for $t$.
+
+## Plot Simulation
+The code below overlays the PDF on reaction time histograms for each option.
+ ```@example stDDM 
+histogram(dist)
+plot!(dist; t_range=range(.301, 3.0, length=100))
+```
+
+# References
+
+Amasino, D.R., Sullivan, N.J., Kranton, R.E. et al. Amount and time exert independent influences on intertemporal choice. Nat Hum Behav 3, 383–392 (2019). https://doi.org/10.1038/s41562-019-0537-2
+
+Barakchian, Z., Beharelle, A.R. & Hare, T.A. Healthy decisions in the cued-attribute food choice paradigm have high test-retest reliability. Sci Rep, (2021). https://doi.org/10.1038/s41598-021-91933-6
+
+Chen, HY., Lombardi, G., Li, SC. et al. Older adults process the probability of winning sooner but weigh it less during lottery decisions. Sci Rep, (2022). https://doi.org/10.1038/s41598-022-15432-y
+
+Maier, S.U., Raja Beharelle, A., Polanía, R. et al. Dissociable mechanisms govern when and how strongly reward attributes affect decisions. Nat Hum Behav 4, 949–963 (2020). https://doi.org/10.1038/s41562-020-0893-y
+
+Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate of information accumulation. Nat Hum Behav 5, 1698–1706 (2021). https://doi.org/10.1038/s41562-021-01154-0
\ No newline at end of file
diff --git a/src/SequentialSamplingModels.jl b/src/SequentialSamplingModels.jl
index b814785e..a280f70e 100644
--- a/src/SequentialSamplingModels.jl
+++ b/src/SequentialSamplingModels.jl
@@ -28,6 +28,7 @@ module SequentialSamplingModels
     import Distributions: rand
     import Distributions: std
     import StatsAPI: params
+    import StatsBase:cor2cov
 
     export AbstractaDDM
     export AbstractCDDM
diff --git a/src/stDDM.jl b/src/stDDM.jl
index bd4fa3d3..b278ac8e 100644
--- a/src/stDDM.jl
+++ b/src/stDDM.jl
@@ -53,9 +53,13 @@ choices,rts = rand(dist, 500)
 # References
 
 Amasino, D.R., Sullivan, N.J., Kranton, R.E. et al. Amount and time exert independent influences on intertemporal choice. Nat Hum Behav 3, 383–392 (2019). https://doi.org/10.1038/s41562-019-0537-2
+
 Barakchian, Z., Beharelle, A.R. & Hare, T.A. Healthy decisions in the cued-attribute food choice paradigm have high test-retest reliability. Sci Rep, (2021). https://doi.org/10.1038/s41598-021-91933-6
+
 Chen, HY., Lombardi, G., Li, SC. et al. Older adults process the probability of winning sooner but weigh it less during lottery decisions. Sci Rep, (2022). https://doi.org/10.1038/s41598-022-15432-y
+
 Lombardi, G., & Hare, T. Piecewise constant averaging methods allow for fast and accurate hierarchical Bayesian estimation of drift diffusion models with time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
+
 Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate of information accumulation. Nat Hum Behav 5, 1698–1706 (2021). https://doi.org/10.1038/s41562-021-01154-0
 """
 
@@ -133,13 +137,9 @@ function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; max_steps=6)
     X = z * α
     deciding = true
     cont = 1
-
     Ρ = [1.0 ρ; ρ 1.0]    
-    # Convert the standard deviations to a diagonal matrix
-    Σ_diag = [η[1] 0; 0 η[2]]
-    # Calculate the covariance matrix from correlation
-    Σ = Σ_diag * Ρ * Σ_diag
-
+    Σ = cor2cov(Ρ,η)
+   
     evidence = rand(MvNormal([ν[1], ν[2]], Σ))
 
     while deciding && cont <= lt

From 20a51cea7d7577212a16ca9d9c294bcf97539f8d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sun, 10 Mar 2024 16:51:54 -0700
Subject: [PATCH 04/10] fix method overwrite

---
 src/type_system.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/type_system.jl b/src/type_system.jl
index bdbead10..bcaacda8 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -120,11 +120,10 @@ Base.broadcastable(x::ContinuousMultivariateSSM) = Ref(x)
 
 Base.length(d::SSM2D) = 2
 
-rand(d::SSM2D; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::SSM2D; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM, n::Int; kwargs...) =
-    rand(Random.default_rng(), d, n; kwargs...)
+rand(Random.default_rng(), d, n; kwargs...)
 
 """
     rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)

From 29e7976a33ddd662c35cdecf42db743a493c46c2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sun, 10 Mar 2024 16:55:02 -0700
Subject: [PATCH 05/10] method overwrite bug

---
 src/type_system.jl | 31 ++++++-------------------------
 1 file changed, 6 insertions(+), 25 deletions(-)

diff --git a/src/type_system.jl b/src/type_system.jl
index bcaacda8..b81a2ba5 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -134,17 +134,17 @@ with more than one choice option.
 # Arguments
 
 - `d::SSM2D`: a 2D sequential sampling model.
-- `n_sim::Int`: the number of simulated choices and rts  
+- `N::Int`: the number of simulated choices and rts  
 
 # Keywords
 
 - `kwargs...`: optional keyword arguments 
 """
-function rand(rng::AbstractRNG, d::SSM2D, n_sim::Int)
-    choice = fill(0, n_sim)
-    rt = fill(0.0, n_sim)
-    for i = 1:n_sim
-        choice[i], rt[i] = rand(rng, d)
+function rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)
+    choice = fill(0, N)
+    rt = fill(0.0, N)
+    for i = 1:N
+        choice[i], rt[i] = rand(rng, d; kwargs...)
     end
     return (; choice, rt)
 end
@@ -192,25 +192,6 @@ pdf(d::SSM2D, data::NamedTuple, args...; kwargs...) =
 
 pdf(d::SSM2D, data::AbstractArray{Real,2}) = pdf(d, Int(data[1]), data[2])
 
-"""
-    rand(rng::AbstractRNG, d::SSM2D, N::Int)
-
-Default method for Generating `n_sim` random choice-rt pairs from a sequential sampling model 
-with more than one choice option.
-
-# Arguments
-- `d::SSM2D`: a 2D sequential sampling model.
-- `N::Int`: the number of simulated choices and rts  
-"""
-function rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)  
-    choice = fill(0, N)  
-    rt = fill(0.0, N)  
-    for i in 1:N  
-        choice[i],rt[i] = rand(rng, d; kwargs...)  
-    end  
-    return (;choice,rt)  
-end
-
 """
     cdf(d::SSM2D, choice::Int, ub=10)
 

From 41c192b36f3ff6f63bbfe41400107764d0c69f05 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sun, 10 Mar 2024 16:58:33 -0700
Subject: [PATCH 06/10] Rm RatcliffDDM.jl

---
 src/SequentialSamplingModels.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/SequentialSamplingModels.jl b/src/SequentialSamplingModels.jl
index cfc7209a..e109c7a8 100644
--- a/src/SequentialSamplingModels.jl
+++ b/src/SequentialSamplingModels.jl
@@ -95,6 +95,5 @@ module SequentialSamplingModels
     include("ext_functions.jl")
     include("ex_gaussian.jl")
     include("poisson_race.jl")
-    include("RatcliffDDM.jl")
     include("stDDM.jl")
 end

From a4f95ed2461780fb6f4c6597df2bfedc15478892 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sun, 10 Mar 2024 17:04:35 -0700
Subject: [PATCH 07/10] Update stDDM documentation

---
 docs/make.jl      | 1 +
 docs/src/stDDM.md | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/docs/make.jl b/docs/make.jl
index 6c5d8c16..52adc5ac 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -33,6 +33,7 @@ makedocs(
             "Muti-attribute attentional drift diffusion Model" => "maaDDM.md",
             "Poisson Race" => "poisson_race.md",
             "Racing Diffusion Model (RDM)" => "rdm.md",
+            "Starting-time Drift Diffusion Model (stDDM)" => "stDDM.md",
             "Wald Model" => "wald.md",
             "Wald Mixture Model" => "wald_mixture.md",
         ],
diff --git a/docs/src/stDDM.md b/docs/src/stDDM.md
index 3a7bee9c..13590a28 100644
--- a/docs/src/stDDM.md
+++ b/docs/src/stDDM.md
@@ -1,4 +1,4 @@
-# Starting-time Diffusion Decision Model (stDDM)
+# Starting-time Drift Diffusion Model (stDDM)
 
 The relative starting time drift diffusion model (stDDM) characterizes the contributions of multiple unique attributes to the rate of evidence accumulation. Compared to the DDM, which assumes a constant evidence accumulation rate within each trial, the stDDM allows different attributes to enter the evidence accumulation process at various time points relative to one another. By doing so, the stDDM quantifies both the weights given to each attribute and their onset times (Amasino et al., 2019; Barakchian et al., 2021; Chen et al., 2022; Maier et al., 2020; Sullivan and Huettel, 2021).
 

From 11ced8e3314d561c973076968a83ba904d9ef62d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Mon, 11 Mar 2024 12:49:33 -0700
Subject: [PATCH 08/10] added intial stddm tests

---
 test/stddm_tests.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 test/stddm_tests.jl

diff --git a/test/stddm_tests.jl b/test/stddm_tests.jl
new file mode 100644
index 00000000..7f47956c
--- /dev/null
+++ b/test/stddm_tests.jl
@@ -0,0 +1,45 @@
+@safetestset "stDDM Tests" begin
+    @safetestset "stDDM predictions" begin
+        using SequentialSamplingModels
+        using Test
+        using Random
+        Random.seed!(8414)
+
+        parms = (ν = [1.0, 1.2], α = 0.8, τ = 0.30 , s = 0.0, z = 0.50)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+        # use validated simulator
+        #Lombardi, G., & Hare, T. Piecewise constant averaging methods allow for fast and 
+        # accurate hierarchical Bayesian estimation of drift diffusion models with 
+        # time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
+        #code to validate
+        #https://github.com/galombardi/method_HtSSM_aDDM/tree/master/RecoveryFitting/stDDM
+        @test mean(choice .== 1) ≈  0.8192000  atol = 1e-2
+        @test mean(rt[choice.==1]) ≈ 0.4332338 atol = 1e-2
+        @test mean(rt[choice.==2]) ≈ 0.4555747 atol = 1e-2
+        @test std(rt[choice.==1]) ≈ 0.1106120 atol = 1e-2
+        @test std(rt[choice.==2]) ≈ 0.1318918 atol = 1e-2
+    end
+
+    @safetestset "simulate" begin
+        using SequentialSamplingModels
+        using Test
+        using Random
+
+        Random.seed!(8477)
+        α = 0.80
+        Δt = 0.0001
+        dist = stDDM(; α, ν = [3, 3], Δt)
+
+        time_steps, evidence = simulate(dist)
+
+        @test time_steps[1] ≈ 0
+        @test length(time_steps) == length(evidence)
+        @test evidence[end] ≈ α atol = 0.05
+
+        dist = stDDM(; α, ν = [-3, -3], Δt)
+        time_steps, evidence = simulate(dist)
+        @test evidence[end] ≈ 0.0 atol = 0.05
+    end
+
+end

From 53ced0aebfbf121d3856c10df50811d3ae985d60 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sat, 16 Mar 2024 22:16:28 -0700
Subject: [PATCH 09/10] Refactor stDDM model and update documentation added
 tests

---
 docs/src/stDDM.md   |   6 +-
 src/stDDM.jl        |  88 ++++++++++-------------
 src/type_system.jl  |   1 +
 test/stddm_tests.jl | 165 +++++++++++++++++++++++++++++++++++++++++---
 4 files changed, 198 insertions(+), 62 deletions(-)

diff --git a/docs/src/stDDM.md b/docs/src/stDDM.md
index 13590a28..b6ff0f01 100644
--- a/docs/src/stDDM.md
+++ b/docs/src/stDDM.md
@@ -19,8 +19,8 @@ Random.seed!(8741)
 ## Create Model Object
 In the code below, we will define parameters for the stDDM and create a model object to store the parameter values. 
 
-### Drift Rates
-The drift rates control the speed and direction with which information accumulates, with one drift rate per attribute (e.g., taste and health, payoff and delay, self and other).
+### Drift Rate
+The drift rate controls the speed and direction in which information accumulates. Here, each drift coefficient indicates the weighting strengths given to the first and second attributes (e.g., taste and health, payoff and delay, self and other), respectively, to the total drift rate in a given trial, where the drift rate accumulates relative evidence in favor of an option.
 ```@example stDDM
 ν = [2.5,2.0]
 ```
@@ -38,7 +38,7 @@ Non-decision time is an additive constant representing encoding and motor respon
 ```
 
 ### starting time
-The starting time parameter \(s\) denotes how much earlier one attribute begins to affect the evidence accumulation process relative to the other(s). If \(s\) is negative, attribute 1 evidence is accumulated before attribute 2 evidence; if \(s\) is positive, attribute 1 evidence is accumulated after attribute 2 evidence. The absolute value of \(s\) indicates the difference in starting times for the two attributes.
+The starting time parameter $s$ denotes how much earlier one attribute begins to affect the evidence accumulation process relative to the other(s). If $s$ is negative, attribute 1 evidence is accumulated before attribute 2 evidence; if $s$ is positive, attribute 1 evidence is accumulated after attribute 2 evidence. The absolute value of $s$ indicates the difference in starting times for the two attributes.
 ```@example stDDM 
 s = 0.10 
 ```
diff --git a/src/stDDM.jl b/src/stDDM.jl
index b278ac8e..9c596db5 100644
--- a/src/stDDM.jl
+++ b/src/stDDM.jl
@@ -5,29 +5,20 @@ An object for the starting-time diffusion decision model.
 
 # Parameters 
 
-- `ν`:  vector of drift rates for attribute one and two
-- `α`:  evidence threshold
-- `τ`:  non-decision time 
+- `ν`:  vector of drift rate weights for attribute one and two
+- `σ`:  diffusion noise
 - `s`:  initial latency bias (positive for attribute two, negative for attribute one)
 - `z`:  initial evidence 
 - `η`:  vector of variability in drift rate for attribute one and two
 - `ρ`:  correlation between drift rate for attributes
-- `σ`:  diffusion noise
-- `Δt`: time step 
+- `α`:  evidence threshold
+- `τ`:  non-decision time 
 
 # Constructors 
 
-    stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
-
-    stDDM(;ν = [0.5,0.6], 
-        α = 1.0, 
-        τ = .300, 
-        s = 0.50, 
-        z = 0.50, 
-        η = [1.0,1.0], 
-        ρ = 0.00,
-        σ = 1,
-        Δt = .001)
+    stDDM(ν, σ, s, z, η, ρ, α, τ)
+
+    stDDM(;ν = [0.5,0.6],σ = 1,s = 0.50, z = 0.50, η = [1.0,1.0], ρ = 0.00, α = 1.0, τ = .300)
         
 # Example 
 
@@ -35,17 +26,16 @@ An object for the starting-time diffusion decision model.
 using SequentialSamplingModels
 
 ν = [0.5, 0.6]
-α = 1.0
-τ = 0.300
+σ = 1
 s = 0.50
 z = 0.50
 η = [1.0, 1.0]
 ρ = 0.00
-σ = 1
-Δt = 0.001
+α = 1.0
+τ = 0.300
 
 # Create stDDM model instance
-dist = stDDM(;ν, α, τ, s, z, η, ρ, σ, Δt)
+dist = stDDM(;ν, σ, s, z, η, ρ, α, τ)
 
 choices,rts = rand(dist, 500)
 ```
@@ -65,38 +55,36 @@ Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate o
 
 mutable struct stDDM{T<:Real} <: AbstractstDDM
     ν::Vector{T}
-    α::T
-    τ::T
+    σ::T
     s::T
     z::T
     η::Vector{T}
     ρ::T
-    σ::T
-    Δt::T
+    α::T
+    τ::T
 end
 
-function stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
-    _, α, τ, s, z,_ , ρ, σ, Δt = promote(ν[1], α, τ, s, z, η[1], ρ, σ, Δt)
+function stDDM(ν, σ, s, z, η, ρ, α, τ)
+    _, σ ,s ,z ,_ ,ρ , α, τ = promote(ν[1],σ, s, z, η[1], ρ, α, τ)
     ν = convert(Vector{typeof(τ)}, ν)
     η = convert(Vector{typeof(τ)}, η)
-    return stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
+    return stDDM(ν, σ, s, z, η, ρ, α, τ)
 end
 
 function stDDM(;ν = [0.5,0.6], 
-    α = 1.0, 
-    τ = .300, 
+    σ = 1,
     s = 0.50, 
     z = 0.50, 
     η = fill(1.0, length(ν)), 
     ρ = 0.0,
-    σ = 1,
-    Δt = .001)
-
-    return stDDM(ν, α, τ, s, z, η, ρ, σ, Δt)
+    α = 1.0, 
+    τ = .300
+    )
+    return stDDM(ν, σ, s, z, η, ρ, α, τ)
 end
 
 function params(d::AbstractstDDM)
-    (d.ν, d.α, d.τ, d.s, d.z, d.η, d.ρ, d.σ, d.Δt)    
+    (d.ν, d.σ, d.s, d.z, d.η, d.ρ, d.α, d.τ)    
 end
 
 get_pdf_type(d::AbstractstDDM) = Approximate
@@ -107,14 +95,16 @@ get_pdf_type(d::AbstractstDDM) = Approximate
 Generate a random choice-rt pair for starting-time diffusion decision model.
 
 # Arguments
+- `rng`: a random number generator
 - `dist`: model object for the starting-time diffusion decision model. 
+- `Δt`: time-step for simulation
 """
-function rand(rng::AbstractRNG, dist::AbstractstDDM)
-    return simulate_trial(rng, dist)
+function rand(rng::AbstractRNG, dist::AbstractstDDM; kwargs...)
+    return simulate_trial(rng, dist; kwargs...)
 end
 
 """
-simulate_trial(rng::AbstractRNG, dist::AbstractstDDM, TMax)
+    simulate_trial(rng::AbstractRNG, dist::AbstractstDDM;  Δt, max_steps)
 
 Generate a single simulated trial from the starting-time diffusion decision model.
 
@@ -122,12 +112,12 @@ Generate a single simulated trial from the starting-time diffusion decision mode
 
 - `rng`: a random number generator
 - `model::AbstractstDDM`: a starting-time diffusion decision model object
+- `Δt`: time-step for simulation
 - `max_steps`: total/max time for simulation
 """
+function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; Δt = .001, max_steps=6)
+    (;ν, σ, s, z, η, ρ, α, τ) = d
 
-function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; max_steps=6)
-    (;ν, α, τ, s, z, η, ρ, σ, Δt) = d
-    
     lt = Int(max_steps / Δt)
     start_step = abs(Int(s / Δt))
     
@@ -149,7 +139,7 @@ function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; max_steps=6)
         noise = rand(rng, Normal(0, σ)) * √(Δt)
 
         X += (evidence[1] * δ1 + evidence[2] * δ2) * Δt + noise
-        
+        # X += (evidence[1] * a1Δ * δ1 + evidence[2] * a1Δ * δ2) * Δt + noise #note something to consider and fix is attribute difference
         if X > α
             choice = 1
             deciding = false
@@ -168,17 +158,19 @@ end
 
 
 """
-    simulate(model::AbstractstDDM)
+    simulate(rng::AbstractRNG, model::AbstractstDDM; Δt)
 
 Returns a matrix containing evidence samples of the stDDM decision process. In the matrix, rows 
 represent samples of evidence per time step and columns represent different accumulators.
 
 # Arguments
 
+- `rng`: a random number generator
 - `model::AbstractstDDM`: a starting-time diffusion decision model diffusion model object
+- `Δt`: time-step for simulation
 """
-function simulate(rng::AbstractRNG, model::AbstractstDDM)
-    (;ν, α, s, z, η, ρ, σ, Δt) = model
+function simulate(rng::AbstractRNG, model::AbstractstDDM; Δt = .001)
+    (;ν, σ, s, z, η, ρ, α, τ) = model
 
     x = α * z
     t = 0.0
@@ -188,11 +180,7 @@ function simulate(rng::AbstractRNG, model::AbstractstDDM)
     start_step = abs(Int(s / Δt))
 
     Ρ = [1.0 ρ; ρ 1.0]    
-    # Convert the standard deviations to a diagonal matrix
-    Σ_diag = [η[1] 0; 0 η[2]]
-    # Calculate the covariance matrix from correlation
-    Σ = Σ_diag * Ρ * Σ_diag
-    
+    Σ = cor2cov(Ρ,η) 
     while (x < α) && (x > 0)
         t += Δt
 
diff --git a/src/type_system.jl b/src/type_system.jl
index b81a2ba5..553a237c 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -148,6 +148,7 @@ function rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)
     end
     return (; choice, rt)
 end
+rand(d::SSM2D, N::Int; kwargs...) = rand(Random.default_rng(), d, N; kwargs...)
 
 """
     logpdf(d::SSM2D, data::NamedTuple) 
diff --git a/test/stddm_tests.jl b/test/stddm_tests.jl
index 7f47956c..8e7ad2fb 100644
--- a/test/stddm_tests.jl
+++ b/test/stddm_tests.jl
@@ -14,12 +14,159 @@
         # time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
         #code to validate
         #https://github.com/galombardi/method_HtSSM_aDDM/tree/master/RecoveryFitting/stDDM
-        @test mean(choice .== 1) ≈  0.8192000  atol = 1e-2
-        @test mean(rt[choice.==1]) ≈ 0.4332338 atol = 1e-2
-        @test mean(rt[choice.==2]) ≈ 0.4555747 atol = 1e-2
-        @test std(rt[choice.==1]) ≈ 0.1106120 atol = 1e-2
-        @test std(rt[choice.==2]) ≈ 0.1318918 atol = 1e-2
+
+        @test mean(choice .== 1) ≈  0.8192  atol = 1e-1
+        @test mean(rt[choice.==1]) ≈ 0.4332 atol = 1e-1
+        @test mean(rt[choice.==2]) ≈ 0.4555 atol = 1e-1
+        @test std(rt[choice.==1]) ≈ 0.1106 atol = 1e-1
+        @test std(rt[choice.==2]) ≈ 0.1318 atol = 1e-1
+
+        parms = (ν = [2, 1], α =1.5, τ = 0.30 , s = 0.20, z = 0.50)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+
+        #R:      0.8883 0.6639 0.6303 0.2377 0.3146
+        #Python: 0.8918 0.6694 0.6285 0.2478 0.3156
+        #Julia:  0.8897 0.6635 0.6298 0.2464 0.3037
+
+        @test mean(choice .== 1)   ≈  0.8883 atol = 1e-1
+        @test mean(rt[choice.==1]) ≈  0.6639  atol = 1e-1
+        @test mean(rt[choice.==2]) ≈  0.6303  atol = 1e-1
+        @test std(rt[choice.==1])  ≈  0.2377 atol = 1e-1
+        @test std(rt[choice.==2])  ≈  0.3146 atol = 1e-1
+
+        parms = (ν = [.7, .9], α =2, τ = 0.40 , s = 0.60, z = 0.40)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+
+        #R:      0.7417 1.1632 1.0788 0.5329 0.6273
+        #Python: 0.7345 1.1686 1.0444 0.5515 0.6125
+        #Julia:  0.7335 1.1642 1.0546 0.533 0.6285
+
+        @test mean(choice .== 1)   ≈  0.7417 atol = 1e-1
+        @test mean(rt[choice.==1]) ≈  1.1632  atol = 1e-1
+        @test mean(rt[choice.==2]) ≈  1.0788  atol = 1e-1
+        @test std(rt[choice.==1])  ≈  0.5329 atol = 1e-1
+        @test std(rt[choice.==2])  ≈  0.6273 atol = 1e-1
+
     end
+    
+#   R:
+#   Rstddm <- function(d_v, d_h, thres, nDT, tIn, bias, sd_n, N, seed = NULL) {
+#   if (!is.null(seed)) set.seed(seed)
+  
+#   T <- 6 # Total time
+#   dt <- 0.001 # Time step
+#   lt <- as.integer(T / dt) # Number of time steps
+#   vec_tHealth <- rep(1, lt)
+#   vec_tVal <- rep(1, lt)
+#   aux <- abs(as.integer(tIn / dt))
+  
+#   # Adjusting health and value vectors based on tIn
+#   if (tIn > 0) {
+#     vec_tVal[1:aux] <- 0
+#   } else if (tIn < 0) {
+#     vec_tHealth[1:aux] <- 0
+#   }
+  
+#   simulate_ddm <- function() {
+#     X <- bias * thres # Initial accumulation value
+#     flag <- FALSE
+#     cont <- 0
+#     vecOut <- T
+    
+#     Sigma = matrix(c(1,0,0,1),2,2)
+#     d_vht <- MASS::mvrnorm(1, c(d_v, d_h), Sigma)
+#     while (!flag && cont < lt) {
+#       noise <- rnorm(1, mean = 0, sd = sd_n) * sqrt(dt)
+#       X <- X + (d_vht[[1]] *1* vec_tVal[cont + 1] + d_vht[[2]] *1* vec_tHealth[cont + 1]) * dt + noise #note 1's are VD and HD
+      
+#       if (X > thres) {
+#         flag <- TRUE
+#         vecOut <- nDT + cont * dt
+#       } else if (X < 0) {
+#         flag <- TRUE
+#         vecOut <- -nDT - cont * dt
+#       }
+#       cont <- cont + 1
+#     }
+    
+#     return(vecOut)
+#   }
+  
+#   results <- numeric(N)
+#   for (i in 1:N) {
+#     results[i] <- simulate_ddm()
+#   }
+  
+#   return(results)
+# }
+#res = Rstddm(d_v = 1, d_h = 1.2, thres = 0.8, nDT = 0.30, tIn = 0.0, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = Rstddm(d_v = 2, d_h = 1, thres = 1.5, nDT = 0.30, tIn = 0.2, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = Rstddm(d_v = .7, d_h = .9, thres = 2, nDT = 0.40, tIn = 0.6, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#
+#   Python:
+# import numpy as np
+# def pystddm(d_v, d_h, thres, nDT, tIn, bias, sd_n, N, seed=None):
+#     # Set a seed for reproducibility if provided
+#     if seed is not None:
+#         np.random.seed(seed)
+    
+#     # Constants and Preparations
+#     T = 6  # Total time
+#     dt = 0.001  # Time step
+#     lt = int(T / dt)  # Number of time steps
+#     vec_tHealth = np.ones(lt)
+#     vec_tVal = np.ones(lt)
+#     aux = abs(int(tIn / dt))
+    
+#     # Adjusting health and value vectors based on tIn
+#     if tIn > 0:
+#         vec_tVal[:aux] = 0
+#     elif tIn < 0:
+#         vec_tHealth[:aux] = 0
+    
+#     def simulate_ddm():
+#         X = bias * thres  # Initial accumulation value
+#         flag = False
+#         cont = 0
+#         vecOut = T
+        
+#         Sigma = np.array([[1, 0], [0, 1]])
+#         d_vht = np.random.multivariate_normal([d_v, d_h], Sigma)
+        
+#         while not flag and cont < lt:
+#             noise = np.random.normal(0, sd_n) * np.sqrt(dt)
+#             X += (d_vht[0] * vec_tVal[cont] + d_vht[1] * vec_tHealth[cont]) * dt + noise
+            
+#             if X > thres:
+#                 flag = True
+#                 vecOut = nDT + cont * dt
+#             elif X < 0:
+#                 flag = True
+#                 vecOut = -nDT - cont * dt
+#             cont += 1
+        
+#         return vecOut
+    
+#     # Apply the simulation for each output element
+#     results = [simulate_ddm() for _ in range(N)]
+    
+#     return results
+
+#res = pystddm(d_v = 1, d_h = 1.2, thres = 0.8, nDT = 0.30, tIn = 0.0, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = pystddm(d_v = 2, d_h = 1, thres = 1.5, nDT = 0.30, tIn = 0.2, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = pystddm(
+#     d_v=.7,    # Weight of the drift for the first attribute
+#     d_h=.9,    # Weight of the drift for the second attribute
+#     thres=2,  # Decision threshold
+#     nDT=0.4,  # Additional decision time
+#     tIn=0.60,  # Time influence of the second attribute
+#     bias=.4,   # Initial value (starting point) of the accumulation
+#     sd_n=1,   # Standard deviation of the noise
+#     N=10000,    # Number of simulations
+#     seed=8414  # Seed for reproducibility
+# )
 
     @safetestset "simulate" begin
         using SequentialSamplingModels
@@ -29,16 +176,16 @@
         Random.seed!(8477)
         α = 0.80
         Δt = 0.0001
-        dist = stDDM(; α, ν = [3, 3], Δt)
+        dist = stDDM(; α, ν = [3, 3])
 
-        time_steps, evidence = simulate(dist)
+        time_steps, evidence = simulate(dist; Δt)
 
         @test time_steps[1] ≈ 0
         @test length(time_steps) == length(evidence)
         @test evidence[end] ≈ α atol = 0.05
 
-        dist = stDDM(; α, ν = [-3, -3], Δt)
-        time_steps, evidence = simulate(dist)
+        dist = stDDM(; α, ν = [-3, -3])
+        time_steps, evidence = simulate(dist; Δt)
         @test evidence[end] ≈ 0.0 atol = 0.05
     end
 

From 13742ba01ae4f749bb6a61e7496c84e7f08f32e2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?=
 <61021880+kiante-fernandez@users.noreply.github.com>
Date: Sat, 16 Mar 2024 22:29:37 -0700
Subject: [PATCH 10/10] Update stDDM.md with parameter changes

---
 docs/src/stDDM.md | 28 +++++++++++-----------------
 1 file changed, 11 insertions(+), 17 deletions(-)

diff --git a/docs/src/stDDM.md b/docs/src/stDDM.md
index b6ff0f01..6ca7dd81 100644
--- a/docs/src/stDDM.md
+++ b/docs/src/stDDM.md
@@ -24,17 +24,10 @@ The drift rate controls the speed and direction in which information accumulates
 ```@example stDDM
 ν = [2.5,2.0]
 ```
-
-### Threshold
-The threshold α represents the amount of evidence required to make a decision.
-```@example stDDM 
-α = 1.5
-```
-
-### Non-Decision Time
-Non-decision time is an additive constant representing encoding and motor response time. 
+### Diffusion Noise
+Diffusion noise is the amount of within trial noise in the evidence accumulation process. 
 ```@example stDDM 
-τ = 0.30
+σ = 1.0
 ```
 
 ### starting time
@@ -56,21 +49,22 @@ Dispersion parameters of the drift rate are drawn from a multivariate normal dis
 ρ = 0.3
 ```
 
-### Diffusion Noise
-Diffusion noise is the amount of within trial noise in the evidence accumulation process. 
+### Threshold
+The threshold α represents the amount of evidence required to make a decision.
 ```@example stDDM 
-σ = 1.0
+α = 1.5
 ```
-### Time Step
-The time step parameter $\Delta t$ is the precision of the discrete time approxmation. 
+
+### Non-Decision Time
+Non-decision time is an additive constant representing encoding and motor response time. 
 ```@example stDDM 
-Δt = .001
+τ = 0.30
 ```
 
 ### stDDM Constructor 
 Now that values have been asigned to the parameters, we will pass them to `stDDM` to generate the model object.
 ```@example stDDM 
-dist = stDDM(;ν, α, τ, s, z, η, ρ, σ, Δt)
+dist = stDDM(;ν, σ, s, z, η, ρ, α, τ,)
 ```
 
 ## Simulate Model