Merge pull request #22 from kiante-fernandez/ssmDev

updated DDM to reflect tuple output

docs/src/DDM.md

@@ -0,0 +1,125 @@
+# Diffusion Decision Model
+
+The Diffusion Decision Model (DDM; Ratcliff et al., 2016) is a model of speeded decision-making in two-choice tasks. The DDM assumes that evidence accumulates over time, starting from a certain position, until it crosses one of two boundaries and triggers the corresponding response (Ratcliff & McKoon, 2008; Ratcliff & Rouder, 1998; Ratcliff & Smith, 2004). Like other Sequential Sampling Models, the DDM comprises psychologically interpretable parameters that collectively form a generative model for reaction time distributions of both responses.
+
+The drift rate (ν) determines the rate at which the accumulation process approaches a decision boundary, representing the relative evidence for or against a specific response. The distance between the two decision boundaries (referred to as the evidence threshold, α) influences the amount of evidence required before executing a response. Non-decision-related components, including perceptual encoding, movement initiation, and execution, are accounted for in the DDM and reflected in the τ parameter. Lastly, the model incorporates a bias in the evidence accumulation process through the parameter z, affecting the starting point of the drift process in relation to the two boundaries. The z parameter in DDM is relative to a (i.e. it ranges from 0 to 1).
+
+One last parameter is the within-trial variability in drift rate (σ), or the diffusion coefficient. The diffusion coefficient is the standard deviation of the evidence accumulation process within one trial. It is a scaling parameter and by convention it is kept fixed. Following Navarro & Fuss, (2009), we use the σ = 1 version.
+
+# Example
+In this example, we will demonstrate how to use the DDM in a generic two alternative forced choice task.
+
+## Load Packages
+The first step is to load the required packages.
+
+```@example DDM
+using SequentialSamplingModels
+using Plots
+using Random
+
+Random.seed!(8741)
+```
+
+## Create Model Object
+In the code below, we will define parameters for the DDM and create a model object to store the parameter values. 
+
+### Drift Rate
+
+The average slope of the information accumulation process. The drift gives information about the speed and direction of the accumulation of information. Typical range: -5 < ν < 5
+
+```@example DDM
+ν=1.0
+```
+
+### Boundary Separation
+
+The amount of information that is considered for a decision. Large values indicates response caution. Typical range: 0.5 < α < 2
+
+```@example DDM 
+α = 0.80
+```
+
+### Non-Decision Time
+
+The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 
+
+```@example DDM 
+τ = 0.30
+```
+
+### 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 DDM 
+z = 0.25
+```
+
+### DDM Constructor 
+
+Now that values have been assigned to the parameters, we will pass them to `DDM` to generate the model object.
+
+```@example DDM 
+dist = DDM(ν, α, τ, z)
+```
+
+## Simulate Model
+
+Now that the model is defined, we will generate $10,000$ choices and reaction times using `rand`. 
+
+ ```@example DDM 
+ choices,rts = rand(dist, 10_000)
+```
+
+## Compute PDF
+The PDF for each observation can be computed as follows:
+
+ ```@example DDM 
+pdf.(dist, choices, rts)
+```
+## Compute Log PDF
+Similarly, the log PDF for each observation can be computed as follows:
+
+ ```@example DDM 
+logpdf.(dist, choices, rts)
+```
+
+## Plot Simulation
+The code below overlays the PDF on reaction time histograms for each option.
+
+ ```@example DDM 
+# rts for option 1
+rts1 = rts[choices .== 1]
+# rts for option 2 
+rts2 = rts[choices .== 2]
+# probability of choosing 1
+p1 = length(rts1) / length(rts)
+t_range = range(.30, 2, length=100)
+# pdf for choice 1
+pdf1 = pdf.(dist, (1,), t_range)
+# pdf for choice 2
+pdf2 = pdf.(dist, (2,), t_range)
+# histogram of retrieval times
+hist = histogram(layout=(2,1), leg=false, grid=false,
+     xlabel="Reaction Time", ylabel="Density", xlims = (0,1.5))
+histogram!(rts1, subplot=1, color=:grey, bins = 200, norm=true, title="Choice 1")
+plot!(t_range, pdf1, subplot=1, color=:darkorange, linewidth=2)
+histogram!(rts2, subplot=2, color=:grey, bins = 150, norm=true, title="Choice 2")
+plot!(t_range, pdf2, subplot=2, color=:darkorange, linewidth=2)
+# weight histogram according to choice probability
+hist[1][1][:y] *= p1
+hist[2][1][:y] *= (1 - p1)
+hist
+```
+
+# References
+
+Navarro, D., & Fuss, I. (2009). Fast and accurate calculations for first-passage times in Wiener diffusion models. https://doi.org/10.1016/J.JMP.2009.02.003
+
+Ratcliff, R., & McKoon, G. (2008). The Diffusion Decision Model: Theory and Data for Two-Choice Decision Tasks. Neural Computation, 20(4), 873–922. https://doi.org/10.1162/neco.2008.12-06-420
+
+Ratcliff, R., & Rouder, J. N. (1998). Modeling Response Times for Two-Choice Decisions. Psychological Science, 9(5), 347–356. https://doi.org/10.1111/1467-9280.00067
+
+Ratcliff, R., & Smith, P. L. (2004). A comparison of sequential sampling models for two-choice reaction time. Psychological Review, 111 2, 333–367. https://doi.org/10.1037/0033-295X.111.2.333
+
+Ratcliff, R., Smith, P. L., Brown, S. D., & McKoon, G. (2016). Diffusion Decision Model: Current Issues and History. Trends in Cognitive Sciences, 20(4), 260–281. https://doi.org/10.1016/j.tics.2016.01.007

docs/src/lca.md

@@ -1,6 +1,6 @@
 # Leaky Competing Accumulator
 
-The Leaky Competing Accumulator (LCA; Brown & Heathcote, 2008) is a sequential sampling model in which evidence for options races independently. The LBA makes an additional simplification that evidence accumulates in a linear and ballistic fashion, meaning there is no intra-trial noise. Instead, evidence accumulates deterministically and linearly until it hits the threshold.
+The Leaky Competing Accumulator (LCA; Usher & McClelland, 2001) is a sequential sampling model in which evidence for options races independently. The LBA makes an additional simplification that evidence accumulates in a linear and ballistic fashion, meaning there is no intra-trial noise. Instead, evidence accumulates deterministically and linearly until it hits the threshold.
 
 # Example
 In this example, we will demonstrate how to use the LBA in a generic two alternative forced choice task. 
@@ -126,4 +126,4 @@ hist
 ```
 # References
 
-Usher, M., & McClelland, J. L. (2001). The time course of perceptual choice: the leaky, competing accumulator model. Psychological Review, 108(3), 550.
+Usher, M., & McClelland, J. L. (2001). The time course of perceptual choice: The leaky, competing accumulator model. Psychological Review, 108 3, 550–592. https://doi.org/10.1037/0033-295X.108.3.550

src/DDM.jl

@@ -4,10 +4,25 @@
     Model object for the Standard Diffusion Decision Model.
 
 # Fields
-    - `ν`: drift rate
-    - `α`: evidence threshold 
-    - `τ`: non-decision time 
-    - `z`: mean starting point
+    - `ν`: drift rate. Average slope of the information accumulation process. The drift gives information about the speed and direction of the accumulation of information. Typical range: -5 < ν < 5
+    - `α`: boundary threshold separation. The amount of information that is considered for a decision. Typical range: 0.5 < α < 2
+    - `τ`: non-decision time. The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 
+    - `z`: starting point. 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 
+
+```julia
+using SequentialSamplingModels
+dist = DDM(ν = 1.0, α = 0.8, τ = 0.3 z = 0.25) 
+choice,rt = rand(dist, 10)
+like = pdf.(dist, choice, rt)
+loglike = logpdf.(dist, choice, rt)
+```
+
+# References
+    
+Ratcliff, R., & McKoon, G. (2008). The Diffusion Decision Model: Theory and Data for Two-Choice Decision Tasks. Neural Computation, 20(4), 873–922.
+
 """
 @concrete mutable struct DDM{T1,T2,T3,T4} <: SequentialSamplingModel
     ν::T1
@@ -16,42 +31,41 @@
     z::T4
 end
 
-# 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)
+Base.broadcastable(x::DDM) = Ref(x)
 
+function params(d::DDM)
+    (d.ν, d.α, d.τ, d.z)    
+end
 
-Base.broadcastable(x::DDM) = Ref(x)
+loglikelihood(d::DDM, data) = sum(logpdf.(d, data...))
 
 """
-DDM(; ν = 0.50,
-    α = 0.08,
+DDM(; ν = 1.0,
+    α = 0.8,
     τ = 0.3
-    z = 0.4,
+    z = 0.25
     )
 
 Constructor for Diffusion Decision Model. 
     
 # Keywords 
-- `ν=.50`: drift rates 
-- `α=0.08`: evidence threshold 
-- `τ=0.3`: non-decision time 
-- `z=0.4`: mean starting point
+    - `ν`: drift rate. Average slope of the information accumulation process. The drift gives information about the speed and direction of the accumulation of information. Typical range: -5 < ν < 5
+    - `α`: boundary threshold separation. The amount of information that is considered for a decision. Typical range: 0.5 < α < 2
+    - `τ`: non-decision time. The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 
+    - `z`: starting point. 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 
 
+```julia
+using SequentialSamplingModels
+dist = DDM(ν = 1.0, α = 0.8, τ = 0.3 z = 0.25) 
+```
 """
-function DDM(; ν = 0.50,
-    α = 0.08,
+function DDM(; ν = 1.00,
+    α = 0.80,
     τ = 0.30,
-    z = 0.4)
-
+    z = 0.25
+    )
     return DDM(ν, α, τ, z)
 end
 
@@ -60,21 +74,18 @@ end
 #  See https://github.com/t-alfers/WienerDiffusionModel.jl                     #
 ################################################################################
 
-function params(d::DDM)
-    (d.ν, d.α, d.τ, d.z)    
-end
 #####################################
 # Probability density function      #
 # Navarro & Fuss (2009)             #
 # Wabersich & Vandekerckhove (2014) #
 #####################################
 
-function pdf(d::DDM, t::Real; ϵ::Real = 1.0e-12)
-    if t >= zero(t)
+function pdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12)
+    if choice == 1
         (ν, α, τ, z) = params(d)
-        return pdf(DDM(-ν, α, τ, 1-z), t; ϵ=ϵ)
+        return pdf(DDM(-ν, α, τ, 1-z), rt; ϵ=ϵ)
     end
-    return pdf(d, t; ϵ=ϵ)
+    return pdf(d, rt; ϵ=ϵ)
 end
 
 # probability density function over the lower boundary
@@ -83,7 +94,7 @@ function pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T<:Real}
     if τ ≥ t
         return T(NaN)
     end
-    u = (t - τ) / α^2
+    u = (t - τ) / α^2 #use normalized time
 
     K_s = 2.0
     K_l = 1 / (π * sqrt(u))
@@ -129,23 +140,35 @@ function _large_time_pdf(u::T, z::T, K::Int) where {T<:Real}
     return π * inf_sum
 end
 
-#logpdf(d::DDM, x::Int, t::Real; ϵ::Real = 1.0e-12) = log(pdf(d, x, t; ϵ=ϵ))
-logpdf(d::DDM, t::Real; ϵ::Real = 1.0e-12) = log(pdf(d, t; ϵ=ϵ))
+logpdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12) = log(pdf(d, choice, rt; ϵ=ϵ))
+#logpdf(d::DDM, t::Real; ϵ::Real = 1.0e-12) = log(pdf(d, t; ϵ=ϵ))
+
+function logpdf(d::DDM, data::T) where {T<:NamedTuple}
+    return sum(logpdf.(d, data...))
+end
+
+function logpdf(dist::DDM, data::Array{<:Tuple,1})
+    LL = 0.0
+    for d in data
+        LL += logpdf(dist, d...)
+    end
+    return LL
+end
 
-loglikelihood(d::DDM, t::Real) = sum(logpdf.(d, t...))
+logpdf(d::DDM, data::Tuple) = logpdf(d, data...)
 
 #########################################
 # Cumulative density function           #
 # Blurton, Kesselmeier, & Gondan (2012) #
 #########################################
 
-function cdf(d::DDM, t::Real; ϵ::Real = 1.0e-12)
-    if  t >= zero(t)
+function cdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12)
+    if choice == 1
         (ν, α, τ, z) = params(d)
-        return cdf(DDM(-ν, α, τ, 1-z), t; ϵ=ϵ)
+        return cdf(DDM(-ν, α, τ, 1-z), rt; ϵ=ϵ)
     end
 
-    return cdf(d, t; ϵ=ϵ)
+    return cdf(d, rt; ϵ=ϵ)
 end
 
 # cumulative density function over the lower boundary
@@ -260,7 +283,7 @@ function rand(rng::AbstractRNG, d::DDM)
     return _rand_rejection(rng, d)
 end
 
-# Rejection-based sampling method (Tuerlinckx et al., 2001 based on Lichters et al., 1995)
+# Rejection-based Method for the Symmetric Wiener Process(Tuerlinckx et al., 2001 based on Lichters et al., 1995)
 # adapted from the RWiener R package, note, here σ = 0.1
 function _rand_rejection(rng::AbstractRNG, d::DDM)
     (ν, α, τ, z) = params(d)
@@ -322,9 +345,13 @@ function _rand_rejection(rng::AbstractRNG, d::DDM)
         dir = start_pos + dir * radius
 
         if (dir + ϵ) > Aupper
-            return total_time + τ
+            rt = total_time + τ
+            choice = 1
+            return choice,rt
         elseif (dir - ϵ) < Alower
-            return -(total_time + τ)
+            rt = total_time + τ
+            choice = 2
+            return choice,rt
         else
             start_pos = dir
             radius = min(abs(Aupper - start_pos), (abs(Alower - start_pos)))
@@ -343,11 +370,12 @@ Generate `n_sim` random choice-rt pairs for the Diffusion Decision Model.
 """
 
 function rand(rng::AbstractRNG, d::DDM, n_sim::Int)
-    rts = fill(0.0, n_sim)
+    choice = fill(0, n_sim)
+    rt = fill(0.0, n_sim)
     for i in 1:n_sim
-        rts[i] =  _rand_rejection(rng, d)
+        choice[i],rt[i] = rand(d)
     end
-    return rts
+    return (choice=choice,rt=rt)
 end
 
 sampler(rng::AbstractRNG, d::DDM) = rand(rng::AbstractRNG, d::DDM)
@@ -418,5 +446,4 @@ sampler(rng::AbstractRNG, d::DDM) = rand(rng::AbstractRNG, d::DDM)
 
 # function params(d::RatcliffDDM)
 #     (d.ν, d.η, d.α, d.z, d.sz, d.τ, d.st, d.σ, d.Δt)    
-# end
-
+# end