Refactor DDM tests and debug/ID bug

src/DDM.jl

@@ -4,14 +4,14 @@
     Model object for the Ratcliff (Full) Diffusion Decision Model.
 
 # Parameters
-    - `ν`: 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).
-    - `η`:  across-trial-variability of drift rate. Typical range: 0 < η < 2. Default is 0.
-    - `sz`: across-trial-variability of starting point. Typical range: 0 < sz < 0.5. Default is 0.
-    - `st`: across-trial-variability of non-decision time. Typical range: 0 < st < 0.2. Default is 0.
-    - `σ`: diffusion noise constant. Default is 1.
+- `ν`: 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).
+- `η`:  across-trial-variability of drift rate. Typical range: 0 < η < 2. Default is 0.
+- `sz`: across-trial-variability of starting point. Typical range: 0 < sz < 0.5. Default is 0.
+- `st`: across-trial-variability of non-decision time. Typical range: 0 < st < 0.2. Default is 0.
+- `σ`: diffusion noise constant. Default is 1.
 
 # Constructors 
 
@@ -40,6 +40,7 @@ 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.
+
 Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85, 59–108. https://doi.org/10.1037/0033-295X.85.2.59
 """
 mutable struct DDM{T<:Real} <: AbstractDDM
@@ -150,7 +151,13 @@ function _pdf_sv(d::DDM, rt; ϵ::Real = 1.0e-12)
         return _pdf(DDM(ν, α, τ, z, 0, 0, 0, σ), rt; ϵ)
     end
 
-    return _pdf(DDM(ν, α, τ, z, 0, 0, 0, σ), rt; ϵ)  + (  ( (α*z*η)^2 - 2*ν*α*z - (ν^2)*(rt-τ) ) / (2*(η^2)*(rt-τ)+2)  ) - log(sqrt((η^2)*(rt-τ)+1)) + ν*α*z + (ν^2)*(rt-τ)*0.5
+    # α *= 2.0
+    # return _pdf(DDM(0, α, τ, z, 0, 0, 0, σ), rt; ϵ)  + (  ( (α*z*η)^2 - 2*ν*α*z - (ν^2)*(rt-τ) ) / (2*(η^2)*(rt-τ)+2)  ) - log(sqrt((η^2)*(rt-τ)+1)) + ν*α*z + (ν^2)*(rt-τ)*0.5
+    #based on:https://github.com/AGhaderi/spatial_attenNCM/blob/main/Neuro-Cognitive%20Models/Models/Hier_Model/lambda_modelt.stan
+    XX = (rt - τ) / α^2 #use normalized time
+
+    return _pdf(DDM(ν, α, τ, z, 0, 0, 0, σ), rt; ϵ)  + (  ( (α*(1-z)*η)^2 - 2*ν*α*(1-z) - (ν^2)*XX ) / (2*(η^2)*XX+2)  ) - log(sqrt((η^2)*XX+1)) + ν*α*(1-z) + (ν^2)*XX*0.5
+    # return _pdf(DDM(0, α, τ, z, 0, 0, 0, σ), rt; ϵ)  + ((α * z * η)^2 - 2 * α * ν * z - (ν^2) * (rt-τ) ) ./ (2 * (η^2) * (rt-τ) .+ 2) - 0.5 * sqrt(η^2 * (rt-τ) .+ 1) - 2 * α
 end
 
 """
@@ -251,7 +258,6 @@ Calculate the 1-dimensional Simpson's numerical integration for a drift diffusio
 https://en.wikipedia.org/wiki/Simpson%27s_rule 
 
 """
-# Simpson's Method one dimentional case
 function _simpson_1D(x::Real, ν::Real, η::Real, α::Real, z::Real, τ::Real, ϵ::Real, lb_z::Real, ub_z::Real, n_sz::Int, lb_t::Real, ub_t::Real, n_st::Int)
 
     n = max(n_st, n_sz)
@@ -291,7 +297,7 @@ function _simpson_1D(x::Real, ν::Real, η::Real, α::Real, z::Real, τ::Real, 
     S = S - y  # the last term should be f(b) and not 2*f(b) so we subtract y
     S = S / ((ub_t - lb_t) + (ub_z - lb_z))  # the right function if pdf_sv()/sz or pdf_sv()/st
 
-    return (ht + hz) * S / 3
+    return ((ht + hz) * S / 3)
 
 end
 
@@ -316,7 +322,6 @@ Calculate the 2-dimensional Simpson's numerical integration for a drift diffusio
 - The function calls `_simpson_1D` to perform the 1-dimensional integrations over each dimension.
 
 """
-# Simpson's Method two dimentional case
 function _simpson_2D(x::Real, ν::Real, η::Real, α::Real, z::Real, τ::Real, ϵ::Real, lb_z::Real, ub_z::Real, n_sz::Int, lb_t::Real, ub_t::Real, n_st::Int)
 
     ht = (ub_t-lb_t)/n_st
@@ -338,33 +343,33 @@ function _simpson_2D(x::Real, ν::Real, η::Real, α::Real, z::Real, τ::Real, 
     S = S - y  # the last term should be f(b) and not 2*f(b) so we subtract y
     S = S / (ub_t - lb_t)
 
-    return ht * S / 3
+    return (ht * S / 3)
 
 end
 
 logpdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12) = log(pdf(d, choice, rt; ϵ))
 
 
-"""
-    cdf(d::DDM, choice, rt; ϵ::Real = 1e-7)
+# """
+#     cdf(d::DDM, choice, rt; ϵ::Real = 1e-7)
 
-Compute the Cumulative Distribution Function (CDF) for the Ratcliff Diffusion model. This function uses 6 Gaussian quadrature for numerical integration.
+# Compute the Cumulative Distribution Function (CDF) for the Ratcliff Diffusion model. This function uses 6 Gaussian quadrature for numerical integration.
 
-# Arguments
-- `d`: an instance of DDM Constructor
-- `choice`: an input representing the choice.
-- `rt`: response time.
-- `ϵ`: a small constant to avoid division by zero, defaults to 1e-7.
+# # Arguments
+# - `d`: an instance of DDM Constructor
+# - `choice`: an input representing the choice.
+# - `rt`: response time.
+# - `ϵ`: a small constant to avoid division by zero, defaults to 1e-7.
 
-# Returns
-- `y`: an array representing the CDF of the Diffusion Decision model.
+# # Returns
+# - `y`: an array representing the CDF of the Diffusion Decision model.
 
-# Reference
-Tuerlinckx, F. (2004). The efficient computation of the cumulative distribution and probability density functions in the diffusion model, Behavior Research Methods, Instruments, & Computers, 36 (4), 702-716.
+# # Reference
+# Tuerlinckx, F. (2004). The efficient computation of the cumulative distribution and probability density functions in the diffusion model, Behavior Research Methods, Instruments, & Computers, 36 (4), 702-716.
 
-# See also
-- Converted from cdfdif.c C script by Joachim Vandekerckhove: https://ppw.kuleuven.be/okp/software/dmat/
-"""
+# # See also
+# - Converted from cdfdif.c C script by Joachim Vandekerckhove: https://ppw.kuleuven.be/okp/software/dmat/
+# """
 # function cdf(d::DDM, choice, rt; ϵ::Real = 1e-7)
 
 #     (ν, α, τ, z, η, sz, st, σ) = params(d)
@@ -391,22 +396,22 @@ Tuerlinckx, F. (2004). The efficient computation of the cumulative distribution
 
 #     return y
 # end
-"""
-    _cdf(d::DDM{T}, choice, rt, prob; ϵ::Real = 1e-7) where {T<:Real}
+# """
+#     _cdf(d::DDM{T}, choice, rt, prob; ϵ::Real = 1e-7) where {T<:Real}
 
-A helper function to compute the Cumulative Distribution Function (CDF) for the Diffusion Decision model.
+# A helper function to compute the Cumulative Distribution Function (CDF) for the Diffusion Decision model.
 
-# Arguments
-- `d`: an instance of DDM Constructor
-- `choice`: an input representing the choice.
-- `rt`: response time.
-- `prob`: a probability value.
-- `ϵ`: a small constant to avoid division by zero, defaults to 1e-7.
+# # Arguments
+# - `d`: an instance of DDM Constructor
+# - `choice`: an input representing the choice.
+# - `rt`: response time.
+# - `prob`: a probability value.
+# - `ϵ`: a small constant to avoid division by zero, defaults to 1e-7.
 
-# Returns
-- `Fnew`: the computed CDF for the given parameters.
+# # Returns
+# - `Fnew`: the computed CDF for the given parameters.
 
-"""
+# """
 # function _cdf(d::DDM, choice, rt, prob; ϵ::Real = 1e-7)
 
 #     (ν, α, τ, z, η, sz, st, σ) = params(d)
@@ -584,36 +589,33 @@ function rand(rng::AbstractRNG, d::DDM)
 end
 
 function _rand_rejection(rng::AbstractRNG, d::DDM; N::Int = 1)
-    (ν, α, τ, z, η, sz, st, σ) = params(d)
+    (;ν, α, τ, z, η, sz, st, σ) = d
 
-    if η == 0
-        η = 1e-16
-    end
+    η = η == 0 ? eps() : η
 
     # Initialize output vectors
-    choice = fill(0, N)
-    rt = fill(0.0, N)
+    T = zeros(N)  
+    XX = fill(0, N)  
 
     # Called sigma in 2001 paper
     D = σ^2 / 2
 
     # Program specifications
     ϵ = eps()  # precision from 1.0 to next double-precision number
-    Δ = ϵ
 
     for n in 1:N
         r1 = randn(rng)
         μ = ν + r1 * η
         bb = z - sz / 2 + sz * rand(rng)
         zz = bb * α
-        finish = 0
-        totaltime = 0
-        startpos = 0
+        finish = false
+        totaltime = 0.0
+        startpos = 0.0
         Aupper = α - zz
         Alower = -zz
         radius = min(abs(Aupper), abs(Alower))
 
-        while finish == 0
+        while !finish
             λ = 0.25 * μ^2 / D + 0.25 * D * π^2 / radius^2
             # eq. formula (13) in 2001 paper with D = sigma^2/2 and radius = Alpha/2
             F = D * π / (radius * μ)
@@ -622,14 +624,14 @@ function _rand_rejection(rng::AbstractRNG, d::DDM; N::Int = 1)
             prob = exp(radius * μ / D)
             prob = prob / (1 + prob)
             dir_ = 2 * (rand(rng) < prob) - 1
-            l = -1
-            s2 = 0
-            s1 = 0
+            l = -1.0
+            s2 = 0.0
+            s1 = 0.0
             while s2 > l
                 s2 = rand(rng)
                 s1 = rand(rng)
-                tnew = 0
-                told = 0
+                tnew = 0.0
+                told = 0.0
                 uu = 0
                 while abs(tnew - told) > ϵ || uu == 0
                     told = tnew
@@ -645,21 +647,21 @@ function _rand_rejection(rng::AbstractRNG, d::DDM; N::Int = 1)
             totaltime += t
             dir_ = startpos + dir_ * radius
             ndt = τ - st / 2 + st * rand(rng)
-            if (dir_ + Δ) > Aupper
-                rt[n] = ndt + totaltime
-                choice[n] = 1
-                finish = 1
-            elseif (dir_ - Δ) < Alower
-                rt[n] = ndt + totaltime
-                choice[n] = 2
-                finish = 1
+            if (dir_ + ϵ) > Aupper
+                T[n] = ndt + totaltime
+                XX[n] = 1
+                finish = true
+            elseif (dir_ - ϵ) < Alower
+                T[n] = ndt + totaltime
+                XX[n] = 2
+                finish = true
             else
                 startpos = dir_
                 radius = minimum(abs.([Aupper, Alower] .- startpos))
             end
         end
     end
-    return (choice=choice,rt=rt)
+    return (choice=XX,rt=T)
 end
 
 function _rand_stochastic(rng::AbstractRNG, d::DDM; N::Int = 1, nsteps::Int=300, step_length::Int=0.01)

test/ddm_tests.jl

@@ -372,68 +372,78 @@
         using Test
         # tested against rtdists
         ## eta
-        test_val1 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 1, .5)
-        @test test_val1 ≈ 2.129834 atol = 1e-1
+        # test_val1 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 1, .5)
+        # @test test_val1 ≈ 2.129834 atol = 1e-1
 
-        test_val2 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 2, .5)
-        @test test_val2 ≈ 0.2889796 atol = 1e-1
+        # test_val2 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 2, .5)
+        # @test test_val2 ≈ 0.2889796 atol = 1e-1
 
-        test_val3 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 1, .35)
-        @test test_val3 ≈ 0.6633653 atol = 1e-1
+        # test_val3 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 1, .35)
+        # @test test_val3 ≈ 0.6633653 atol = 1e-1
 
-        test_val4 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 2, .35)
-        @test test_val4 ≈ 0.4449858 atol = 1e-1
+        # test_val4 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.08, sz = 0.0, st = 0.0, σ = 1.0), 2, .35)
+        # @test test_val4 ≈ 0.4449858 atol = 1e-1
         #try a bunch of combinations from values from rtdists
-        sz_values = [0.00,0.05, 0.1, 0.2, 0.3]
-        st_values = [0.00,0.05, 0.1, 0.2, 0.3]
-        η_values  =  [0.00,0.08, 0.9, 0.13, 0.16]
+        sz_values = [0.00, 0.05, 0.10, 0.30]
+        st_values = [0.00, 0.05, 0.10, 0.30]
+        η_values  = [0.00, 0.08, 0.16]
+
         combinations = [(sz, st, η) for sz in sz_values for st in st_values for η in η_values]
 
-        test_vals_upper = [2.131129, 2.129384, 2.124152, 2.101896, 2.065028, 2.540625, 
-        2.538263, 2.531186, 2.501185, 2.45186, 3.025843, 3.022579, 3.01281, 
-        2.971536, 2.90416, 2.866027, 2.865954, 2.865732, 2.864744, 2.862934, 
-        1.910684, 1.910636, 1.910488, 1.909829, 1.908623, 2.129834, 2.128092, 
-        2.122868, 2.100647, 2.063835, 2.539411, 2.537052, 2.529984, 2.500024, 
-        2.450762, 3.024903, 3.021643, 3.011883, 2.970646, 2.903329, 2.865857, 
-        2.865785, 2.865565, 2.864583, 2.862783, 1.910572, 1.910523, 1.910377, 
-        1.909722, 1.908522, 1.983902, 1.98245, 1.978094, 1.95952, 1.9286, 
-        2.399602, 2.397588, 2.391552, 2.365916, 2.323594, 2.913938, 2.911015, 
-        2.902262, 2.86524, 2.804655, 2.846349, 2.846329, 2.846268, 2.845959, 
-        2.845273, 1.897566, 1.897553, 1.897512, 1.897306, 1.896849, 2.127715, 
-        2.125978, 2.120768, 2.098603, 2.061881, 2.537423, 2.535069, 2.528017, 
-        2.498121, 2.448964, 3.023364, 3.020108, 3.010363, 2.969186, 2.901967, 
-        2.86558, 2.865508, 2.865291, 2.864319, 2.862536, 1.910387, 1.910339, 
-        1.910194, 1.909546, 1.908358, 2.125965, 2.124231, 2.119032, 2.096913, 
-        2.060266, 2.53578, 2.53343, 2.526391, 2.496549, 2.447477, 3.02209, 
-        3.018839, 3.009105, 2.967979, 2.900839, 2.865351, 2.86528, 2.865064, 
-        2.864101, 2.862332, 1.910234, 1.910187, 1.910043, 1.909401, 1.908221]
-
-        test_vals_lower = [0.288417, 0.288327, 0.288055, 0.286843, 0.284638, 0.343836, 
-        0.343781, 0.343611, 0.342803, 0.341161, 0.409503, 0.409585, 0.409821, 
-        0.410666, 0.411506, 0.387875, 0.389054, 0.3926, 0.407855, 0.43373, 
-        0.258583, 0.25937, 0.261734, 0.271903, 0.289153, 0.28898, 0.288889, 
-        0.288615, 0.287391, 0.285169, 0.344436, 0.34438, 0.344207, 0.343387, 
-        0.341724, 0.410134, 0.410215, 0.410447, 0.411278, 0.412097, 0.388387, 
-        0.389566, 0.393111, 0.408359, 0.434222, 0.258925, 0.259711, 0.262074, 
-        0.272239, 0.289481, 0.354834, 0.354647, 0.354081, 0.351634, 0.34743, 
-        0.415484, 0.415315, 0.414806, 0.412568, 0.408615, 0.486035, 0.485984, 
-        0.485823, 0.484998, 0.483146, 0.451208, 0.452322, 0.455669, 0.470067, 
-        0.494477, 0.300805, 0.301548, 0.303779, 0.313378, 0.329651, 0.289902, 
-        0.28981, 0.289531, 0.28829, 0.286039, 0.345418, 0.34536, 0.345183, 
-        0.344343, 0.342648, 0.411169, 0.411247, 0.411475, 0.412283, 0.413065, 
-        0.389227, 0.390406, 0.393948, 0.409185, 0.435029, 0.259485, 0.26027, 
-        0.262632, 0.27279, 0.290019, 0.290664, 0.290571, 0.290289, 0.289034, 
-        0.286759, 0.346231, 0.346172, 0.34599, 0.345134, 0.343412, 0.412025, 
-        0.412102, 0.412325, 0.413114, 0.413866, 0.389923, 0.391101, 0.394641, 
-        0.409869, 0.435697, 0.259949, 0.260734, 0.263094, 0.273246, 0.290465]
+        # for (sz, st, η) in combinations println("$sz, $st, $η") end
+
+        test_vals_upper = [1.81597, 1.814461, 1.809956, 1.974649, 1.973321, 1.969353, 
+        2.123612, 2.122583, 2.119502, 1.624452, 1.624509, 1.62468, 1.815153, 
+        1.813646, 1.809147, 1.973663, 1.972337, 1.968375, 2.122466, 2.121439, 
+        2.118364, 1.624535, 1.624593, 1.624764, 1.812704, 1.811203, 1.806723, 
+        1.970707, 1.969387, 1.965443, 2.119034, 2.118012, 2.114953, 1.624786, 
+        1.624843, 1.625014, 1.785586, 1.784153, 1.779875, 1.938027, 1.936776, 
+        1.933035, 2.081122, 2.080158, 2.077274, 1.627546, 1.627602, 1.627769]
+
+        test_vals_lower = [0.074023, 0.074416, 0.075599, 0.080491, 0.080889, 0.082087, 
+        0.086563, 0.08696, 0.088155, 0.066216, 0.066472, 0.067242, 0.074072, 
+        0.074465, 0.075648, 0.080556, 0.080954, 0.082151, 0.086653, 0.087051, 
+        0.088245, 0.066413, 0.066669, 0.067439, 0.074218, 0.074611, 0.075792, 
+        0.080751, 0.081149, 0.082345, 0.086923, 0.08732, 0.088513, 0.067004, 
+        0.06726, 0.068031, 0.075803, 0.076192, 0.077358, 0.082871, 0.083265, 
+        0.084447, 0.089867, 0.09026, 0.091442, 0.073773, 0.074034, 0.074818]
 
         for (i, (sz, st, η)) in enumerate(combinations)
-            # Compute the pdf for the current values of sz, st, and η
-            pdf_value = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η=η, sz=sz, st=st, σ = 1.0), 1, .5)
-            @test test_vals_upper[i] ≈ pdf_value atol = 1e-1
-            pdf_value = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η=η, sz=sz, st=st, σ = 1.0), 2, .5)
-            @test test_vals_lower[i] ≈ pdf_value atol = 1e-1
-        end
+            # Print out the current parameter set
+           println("Parameter set [$i]: sz = $sz, st = $st, η = $η")
+           # Compute the pdf for the current values of sz, st, and η
+           pdf_value = pdf(DDM(;ν=2.0, α=1.6, z=.5, τ=.2, η=η, sz=sz, st=st, σ = 1.0), 1, .5)
+           println("Test value upper[$i]: ", test_vals_upper[i])
+           println("PDF value upper: ", pdf_value)
+           @test test_vals_upper[i] ≈ pdf_value atol = 1e-1
+           pdf_value = pdf(DDM(;ν=2.0, α=1.6, z=.5, τ=.2, η=η, sz=sz, st=st, σ = 1.0), 2, .5)
+           println("Test value lower[$i]: ", test_vals_lower[i])
+           println("PDF value lower: ", pdf_value)
+           @test test_vals_lower[i] ≈ pdf_value atol = 1e-1
+       end
+
+#     sz_values = [0.00]
+#     st_values = [0.00]
+#     η_values  = [0.00, 0.08, 0.16]
+
+#     combinations = [(sz, st, η) for sz in sz_values for st in st_values for η in η_values]
+
+#     test_vals_upper = [1.81597, 1.814461, 1.809956]
+#     test_vals_lower = [0.074023, 0.074416, 0.075599]
+
+#     for (i, (sz, st, η)) in enumerate(combinations)
+#         # Print out the current parameter set
+#        println("Parameter set [$i]: sz = $sz, st = $st, η = $η")
+#        # Compute the pdf for the current values of sz, st, and η
+#        pdf_value = pdf(DDM(;ν=2.0, α=1.6, z=.5, τ=.2, η=η, sz=sz, st=st, σ = 1.0), 1, .5)
+#        println("Test value upper[$i]: ", test_vals_upper[i])
+#        println("PDF value upper: ", pdf_value)
+#        @test test_vals_upper[i] ≈ pdf_value atol = 1e-3
+#        pdf_value = pdf(DDM(;ν=2.0, α=1.6, z=.5, τ=.2, η=η, sz=sz, st=st, σ = 1.0), 2, .5)
+#        println("Test value lower[$i]: ", test_vals_lower[i])
+#        println("PDF value lower: ", pdf_value)
+#        @test test_vals_lower[i] ≈ pdf_value atol = 1e-3
+#    end
 
         # 0.05, 0.05, 0.08
         # pdf_value = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η=.08, sz=.05, st=.05, σ = 1.0), 1, .5)