From 733cf793015f13bd00791f0e8ce89ea46649e82c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?= <61021880+kiante-fernandez@users.noreply.github.com> Date: Fri, 12 Jul 2024 18:23:32 +0200 Subject: [PATCH 1/3] new interface ratcliff DDM --- src/DDM.jl | 517 ++++++++++++++++++++++++++++------------------ test/ddm_tests.jl | 194 ++++++++--------- 2 files changed, 417 insertions(+), 294 deletions(-) diff --git a/src/DDM.jl b/src/DDM.jl index b9e09f50..9260654f 100644 --- a/src/DDM.jl +++ b/src/DDM.jl @@ -9,22 +9,26 @@ Model object for the standard Drift Diffusion Model. - `α`: boundary threshold separation. The amount of information that is considered for a decision. Typical range: 0.5 < α < 2 - `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). - `τ`: non-decision time. The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 +- `η`: 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 Two constructors are defined below. The first constructor uses positional arguments, and is therefore order dependent: - DDM(ν, α, z, τ) + DDM(ν, α, z, τ, η, sz, st, σ) The second constructor uses keywords with default values, and is not order dependent: - DDM(; ν = 1.00, α = 0.80, τ = 0.30, z = 0.50) + DDM(; ν = 1.00, α = 0.80, τ = 0.30, z = 0.50, η = 0.16, sz = 0.05, st = 0.10, σ = 1.0) # Example ```julia using SequentialSamplingModels -dist = DDM(ν = 1.0, α = 0.8, τ = 0.3, z = 0.25) +dist = DDM(ν = 1.00, α = 0.80, τ = 0.30, z = 0.50, η = 0.16, sz = 0.05, st = 0.10, σ = 1.0) choice,rt = rand(dist, 10) like = pdf.(dist, choice, rt) loglike = logpdf.(dist, choice, rt) @@ -33,24 +37,126 @@ 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. """ mutable struct DDM{T <: Real} <: AbstractDDM ν::T α::T z::T τ::T + η::T + sz::T + st::T + σ::T end -function DDM(ν, α, z, τ) - return DDM(promote(ν, α, z, τ)...) +function DDM(ν, α, z, τ, η, sz, st, σ) + return DDM(promote(ν, α, z, τ, η, sz, st, σ)...) end function params(d::DDM) - return (d.ν, d.α, d.z, d.τ) + (d.ν, d.α, d.z, d.τ, d.η, d.sz, d.st, d.σ) end -function DDM(; ν = 1.00, α = 0.80, τ = 0.30, z = 0.50) - return DDM(ν, α, z, τ) +function DDM(; ν = 1.00, α = 0.80, τ = 0.30, z = 0.50, η = 0.16, sz = 0.05, st = 0.10, σ = 1.0) + return DDM(ν, α, z, τ, η, sz, st, σ) +end + +function pdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12) + if choice == 1 + (ν, α, z, τ, η, sz, st, σ) = params(d) + # Transform ν, z if rt is upper bound response + return _pdf_Full(DDM(-ν, α, 1 - z, τ, η, sz, st, σ), rt; ϵ) + end + return _pdf_Full(d, rt; ϵ) +end + +""" + _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, simps_err::Real = 1e-6) where {T<:Real} + +Calculate the probability density function (PDF) for a Diffusion Decision Model (DDM) object. This +function applies numerical integration to account for variability in non-decision time and bias, as suggested +by Ratcliff and Tuerlinckx (2002). + +# Arguments +- `d::DDM{T}`: a DDM distribution object +- `rt`: reaction time. + +# Optional arguments +- `ϵ::Real`: a small constant to prevent divide by zero errors, default is 1.0e-12. +- `simps_err::Real = 1e-3`: Error tolerance for the numerical integration. + +""" +function _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, simps_err::Real = 1e-6) where {T<:Real} + (ν, α, z, τ, η, sz, st, σ) = params(d) + + # Check if parameters are valid + if (z < 0 || z > 1 || α <= 0 || τ < 0 || st < 0 || η < 0 || sz < 0 || sz > 1 || + (rt - (τ - st/2)) < 0 || (z + sz/2) > 1 || (z - sz/2) < 0 || (τ - st/2) < 0) + return zero(T) + end + + rt = abs(rt) + + if st < 1e-6 + st = 0 + end + if sz < 1e-6 + sz = 0 + end + + if sz == 0 + if st == 0 # η=0, sz=0, st=0 + return _pdf_sv(DDM(ν, α, z, τ, η, 0, 0, σ), rt - τ; ϵ) + else # η=0, sz=0, st≠0 + f = y -> _pdf_sv(DDM(ν, α, z, y[1], η, 0, 0, σ), rt - y[1]; ϵ) + result, _ = hcubature(f, (τ-st/2,), (τ+st/2,); rtol=simps_err, atol=simps_err) + return result / st # Normalize by the integration range + end + else # sz≠0 + if st == 0 # η=0, sz≠0, st=0 + f = y -> _pdf_sv(DDM(ν, α, y[1], τ, η, 0, 0, σ), rt - τ; ϵ) + result, _ = hcubature(f, (z-sz/2,), (z+sz/2,); rtol=simps_err, atol=simps_err) + return result / sz + else # η=0, sz≠0, st≠0 + f = y -> _pdf_sv(DDM(ν, α, y[1], y[2], η, 0, 0, σ), rt - y[2]; ϵ) + result, _ = hcubature(f, (z-sz/2, τ-st/2), (z+sz/2, τ+st/2); rtol=simps_err, atol=simps_err) + return result / (sz * st) + end + end +end + +""" + _pdf_sv(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} + +Computes the Probability Density Function (PDF) for a given Diffusion Decision Model +with across-trial variability in drift-rate. This function uses analytic integration of the likelihood function +for variability in drift-rate. + +# Arguments +- `d::DDM`: a DDM distribution constructor object +- `rt`: Reaction time for which the PDF is to be computed. + +# Returns +- Returns the computed PDF value. + +""" +function _pdf_sv(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} + (ν, α, z, τ, η, σ) = params(d) + + if t <= 0 + return zero(T) + end + + if η == 0 + return _pdf(DDM(ν, α, z, τ, 0, 0, 0, σ), t; ϵ) + end + + XX = t / (α^2) #use normalized time + p = _pdf(DDM(ν, α, z, τ, 0, 0, 0, σ), XX; ϵ) #get f(t|0,1,w) + # convert to f(t|v,a,w) + return exp(log(p) + ((α*z*η)^2 - 2*α*ν*z - (ν^2)*XX)/(2*(η^2)*XX+2)) / sqrt((η^2)*XX+1) / (α^2) + end ################################################################################ @@ -64,21 +170,33 @@ end # Wabersich & Vandekerckhove (2014) # ##################################### -function pdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12) - if choice == 1 - (ν, α, z, τ) = params(d) - return _pdf(DDM(-ν, α, 1 - z, τ), rt; ϵ) - end - return _pdf(d, rt; ϵ) -end +""" + _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} + +Computes the Probability Density Function (PDF) for a given Drift Diffusion Model (DDM). +The function uses normalized time and applies an infinite sum algorithm. The implementation +is based on the work of Navarro & Fuss (2009) and Wabersich & Vandekerckhove (2014). + +# Arguments +- `d::DDM{T}`: A DDM object containing the parameters of the Diffusion Model. +- `t::Real`: Time for which the PDF is to be computed. + +# Optional Arguments +- `ϵ::Real = 1.0e-12`: A very small number representing machine epsilon. + +# Returns +- Returns the computed PDF value. -# probability density function over the lower boundary +# See also: + - Converted from WienerDiffusionModel.jl repository orginally by Tobias Alfers: https://github.com/t-alfers/WienerDiffusionModel.jl +""" function _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} - (ν, α, z, τ) = params(d) - if τ ≥ t - return T(NaN) + (ν, α, z) = params(d) + + if t <= 0 + return zero(T) end - u = (t - τ) / α^2 #use normalized time + u = t K_s = 2.0 K_l = 1 / (π * sqrt(u)) @@ -91,7 +209,7 @@ function _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} K_s = max(2 + sqrt(-2u * log(2ϵ * sqrt(2 * π * u))), sqrt(u) + 1) end - p = exp((-α * z * ν) - (0.5 * (ν^2) * (t - τ))) / (α^2) + p = exp((-α * z * ν) - (0.5 * (ν^2) * (t))) / (α^2) # decision rule for infinite sum algorithm if K_s < K_l @@ -146,166 +264,166 @@ logpdf(d::DDM, data::Tuple) = logpdf(d, data...) # Blurton, Kesselmeier, & Gondan (2012) # ######################################### -function cdf(d::DDM, choice::Int, rt::Real = 10; ϵ::Real = 1.0e-12) - if choice == 1 - (ν, α, z, τ) = params(d) - return _cdf(DDM(-ν, α, 1 - z, τ), rt; ϵ) - end - - return _cdf(d, rt; ϵ) -end - -# cumulative density function over the lower boundary -function _cdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} - if d.τ ≥ t - return T(NaN) - end - - K_l = _K_large(d, t; ϵ) - K_s = _K_small(d, t; ϵ) - - if K_l < 10 * K_s - return _Fl_lower(d, K_l, t) - end - return _Fs_lower(d, K_s, t) -end - -# Large time representation of lower subdistribution -function _Fl_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} - (ν, α, z, τ) = params(d) - F = zero(T) - K_series = K:-1:1 - for k in K_series - F -= ( - k / (ν^2 + k^2 * π^2 / (α^2)) * - exp(-ν * α * z - 0.5 * ν^2 * (t - τ) - 0.5 * k^2 * π^2 / (α^2) * (t - τ)) * - sin(π * k * z) - ) - end - return _P_upper(ν, α, z) + 2 * π / (α^2) * F -end - -# Small time representation of the upper subdistribution -function _Fs_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} - (ν, α, z, τ) = params(d) - if abs(ν) < sqrt(eps(T)) - return _Fs0_lower(d, K, t) - end - - sqt = sqrt(t - τ) - - S1 = zero(T) - S2 = zero(T) - K_series = K:-1:1 - - for k in K_series - S1 += ( - _exp_pnorm(2 * ν * α * k, -sign(ν) * (2 * α * k + α * z + ν * (t - τ)) / sqt) - _exp_pnorm( - -2 * ν * α * k - 2 * ν * α * z, - sign(ν) * (2 * α * k + α * z - ν * (t - τ)) / sqt - ) - ) - - S2 += ( - _exp_pnorm(-2 * ν * α * k, sign(ν) * (2 * α * k - α * z - ν * (t - τ)) / sqt) - _exp_pnorm( - 2 * ν * α * k - 2 * ν * α * z, - -sign(ν) * (2 * α * k - α * z + ν * (t - τ)) / sqt - ) - ) - end - - return _P_upper(ν, α, z) + - sign(ν) * ( - ( - cdf(Normal(), -sign(ν) * (α * z + ν * (t - τ)) / sqt) - - _exp_pnorm(-2 * ν * α * z, sign(ν) * (α * z - ν * (t - τ)) / sqt) - ) + - S1 + - S2 - ) -end - -# Zero drift version -function _Fs0_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} - (_, α, z, τ) = params(d) - F = zero(T) - K_series = K:-1:0 - for k in K_series - F -= ( - cdf(Distributions.Normal(), (-2 * k - 2 + z) * α / sqrt(t - τ)) + - cdf(Distributions.Normal(), (-2 * k - z) * α / sqrt(t - τ)) - ) - end - return 2 * F -end -# Number of terms required for large time representation -function _K_large(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} - (ν, α, z, τ) = params(d) - x = t - τ - sqrtL1 = sqrt(1 / x) * α / π - sqrtL2 = sqrt( - max( - 1, - -2 / x * α * α / π / π * - (log(ϵ * π * x / 2 * (ν * ν + π * π / α / α)) + ν * α * z + ν * ν * x / 2) - ), - ) - return ceil(Int, max(sqrtL1, sqrtL2)) -end - -# Number of terms required for small time representation -function _K_small(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} - (ν, α, z, τ) = params(d) - if abs(ν) < sqrt(eps(T)) - return ceil( - Int, - max( - 0, - z / 2 - - sqrt(t - τ) / (2 * α) * quantile(Normal(), max(0, min(1, ϵ / (2 - 2 * z)))) - ) - ) - end - if ν > 0 - return _K_small(DDM(-ν, α, z, τ), t; ϵ = exp(-2 * α * z * ν) * ϵ) - end - S2 = z - 1 + 1 / (2 * ν * α) * log(ϵ / 2 * (1 - exp(2 * ν * α))) - S3 = (0.535 * sqrt(2 * (t - τ)) + ν * (t - τ) + α * z) / (2 * α) - S4 = - z / 2 - - sqrt(t - τ) / (2 * α) * quantile( - Normal(), - max( - 0, - min( - 1, - ϵ * α / (0.3 * sqrt(2 * π * (t - τ))) * - exp(ν^2 * (t - τ) / 2 + ν * α * z) - ) - ) - ) - return ceil(Int, max(0, S2, S3, S4)) -end -# Probability for absorption at upper barrier -function _P_upper(ν::T, α::T, z::T) where {T <: Real} - e = exp(-2 * ν * α * (1 - z)) - if isinf(e) - return 1.0 - end - if abs(e - 1) < sqrt(eps(T)) - return 1 - z - end - return (1 - e) / (exp(2 * ν * α * z) - e) -end - -# Calculates exp(a) * pnorm(b) using an approximation by Kiani et al. (2008) -function _exp_pnorm(a::T, b::T) where {T <: Real} - r = exp(a) * cdf(Distributions.Normal(), b) - if isnan(r) && b < -5.5 - r = (1 / sqrt(2)) * exp(a - b^2 / 2) * (0.5641882 / (b^3) - 1 / (b * sqrt(π))) - end - return r -end +# function cdf(d::DDM, choice::Int, rt::Real = 10; ϵ::Real = 1.0e-12) +# if choice == 1 +# (ν, α, z, τ) = params(d) +# return _cdf(DDM(-ν, α, 1 - z, τ), rt; ϵ) +# end + +# return _cdf(d, rt; ϵ) +# end + +# # cumulative density function over the lower boundary +# function _cdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} +# if d.τ ≥ t +# return T(NaN) +# end + +# K_l = _K_large(d, t; ϵ) +# K_s = _K_small(d, t; ϵ) + +# if K_l < 10 * K_s +# return _Fl_lower(d, K_l, t) +# end +# return _Fs_lower(d, K_s, t) +# end + +# # Large time representation of lower subdistribution +# function _Fl_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} +# (ν, α, z, τ) = params(d) +# F = zero(T) +# K_series = K:-1:1 +# for k in K_series +# F -= ( +# k / (ν^2 + k^2 * π^2 / (α^2)) * +# exp(-ν * α * z - 0.5 * ν^2 * (t - τ) - 0.5 * k^2 * π^2 / (α^2) * (t - τ)) * +# sin(π * k * z) +# ) +# end +# return _P_upper(ν, α, z) + 2 * π / (α^2) * F +# end + +# # Small time representation of the upper subdistribution +# function _Fs_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} +# (ν, α, z, τ) = params(d) +# if abs(ν) < sqrt(eps(T)) +# return _Fs0_lower(d, K, t) +# end + +# sqt = sqrt(t - τ) + +# S1 = zero(T) +# S2 = zero(T) +# K_series = K:-1:1 + +# for k in K_series +# S1 += ( +# _exp_pnorm(2 * ν * α * k, -sign(ν) * (2 * α * k + α * z + ν * (t - τ)) / sqt) - _exp_pnorm( +# -2 * ν * α * k - 2 * ν * α * z, +# sign(ν) * (2 * α * k + α * z - ν * (t - τ)) / sqt +# ) +# ) + +# S2 += ( +# _exp_pnorm(-2 * ν * α * k, sign(ν) * (2 * α * k - α * z - ν * (t - τ)) / sqt) - _exp_pnorm( +# 2 * ν * α * k - 2 * ν * α * z, +# -sign(ν) * (2 * α * k - α * z + ν * (t - τ)) / sqt +# ) +# ) +# end + +# return _P_upper(ν, α, z) + +# sign(ν) * ( +# ( +# cdf(Normal(), -sign(ν) * (α * z + ν * (t - τ)) / sqt) - +# _exp_pnorm(-2 * ν * α * z, sign(ν) * (α * z - ν * (t - τ)) / sqt) +# ) + +# S1 + +# S2 +# ) +# end + +# # Zero drift version +# function _Fs0_lower(d::DDM{T}, K::Int, t::Real) where {T <: Real} +# (_, α, z, τ) = params(d) +# F = zero(T) +# K_series = K:-1:0 +# for k in K_series +# F -= ( +# cdf(Distributions.Normal(), (-2 * k - 2 + z) * α / sqrt(t - τ)) + +# cdf(Distributions.Normal(), (-2 * k - z) * α / sqrt(t - τ)) +# ) +# end +# return 2 * F +# end +# # Number of terms required for large time representation +# function _K_large(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} +# (ν, α, z, τ) = params(d) +# x = t - τ +# sqrtL1 = sqrt(1 / x) * α / π +# sqrtL2 = sqrt( +# max( +# 1, +# -2 / x * α * α / π / π * +# (log(ϵ * π * x / 2 * (ν * ν + π * π / α / α)) + ν * α * z + ν * ν * x / 2) +# ), +# ) +# return ceil(Int, max(sqrtL1, sqrtL2)) +# end + +# # Number of terms required for small time representation +# function _K_small(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} +# (ν, α, z, τ) = params(d) +# if abs(ν) < sqrt(eps(T)) +# return ceil( +# Int, +# max( +# 0, +# z / 2 - +# sqrt(t - τ) / (2 * α) * quantile(Normal(), max(0, min(1, ϵ / (2 - 2 * z)))) +# ) +# ) +# end +# if ν > 0 +# return _K_small(DDM(-ν, α, z, τ), t; ϵ = exp(-2 * α * z * ν) * ϵ) +# end +# S2 = z - 1 + 1 / (2 * ν * α) * log(ϵ / 2 * (1 - exp(2 * ν * α))) +# S3 = (0.535 * sqrt(2 * (t - τ)) + ν * (t - τ) + α * z) / (2 * α) +# S4 = +# z / 2 - +# sqrt(t - τ) / (2 * α) * quantile( +# Normal(), +# max( +# 0, +# min( +# 1, +# ϵ * α / (0.3 * sqrt(2 * π * (t - τ))) * +# exp(ν^2 * (t - τ) / 2 + ν * α * z) +# ) +# ) +# ) +# return ceil(Int, max(0, S2, S3, S4)) +# end +# # Probability for absorption at upper barrier +# function _P_upper(ν::T, α::T, z::T) where {T <: Real} +# e = exp(-2 * ν * α * (1 - z)) +# if isinf(e) +# return 1.0 +# end +# if abs(e - 1) < sqrt(eps(T)) +# return 1 - z +# end +# return (1 - e) / (exp(2 * ν * α * z) - e) +# end + +# # Calculates exp(a) * pnorm(b) using an approximation by Kiani et al. (2008) +# function _exp_pnorm(a::T, b::T) where {T <: Real} +# r = exp(a) * cdf(Distributions.Normal(), b) +# if isnan(r) && b < -5.5 +# r = (1 / sqrt(2)) * exp(a - b^2 / 2) * (0.5641882 / (b^3) - 1 / (b * sqrt(π))) +# end +# return r +# end """ rand(dist::DDM) @@ -322,34 +440,35 @@ end # 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) + (ν, α, z, τ, η, sz, st, σ) = params(d) ϵ = 1.0e-15 - D = 0.005 # = 2*σ^2 => 1/200 - zn = (z * α) / 10 # absolute bias! - αn = α / 10 - νn = ν / 10 + D = σ^2 / 2 + + bb = z - sz / 2 + sz * rand(rng) + zz = (bb * α) + νν = ν + randn(rng) * η total_time = 0.0 start_pos = 0.0 - Aupper = αn - zn - Alower = -zn + Aupper = α - zz + Alower = -zz radius = min(abs(Aupper), abs(Alower)) λ = 0.0 F = 0.0 prob = 0.0 while true - if νn == 0 + if νν == 0 λ = (0.25D * π^2) / (radius^2) F = 1.0 prob = 0.5 else - λ = ((0.25 * νn^2) / D) + ((0.25 * D * π^2) / (radius^2)) - F = (D * π) / (radius * νn) + λ = ((0.25 * νν^2) / D) + ((0.25 * D * π^2) / (radius^2)) + F = (D * π) / (radius * νν) F = F^2 / (1 + F^2) - prob = exp((radius * νn) / D) + prob = exp((radius * νν) / D) prob = prob / (1 + prob) end @@ -379,13 +498,14 @@ function _rand_rejection(rng::AbstractRNG, d::DDM) total_time += abs(log(s1)) / λ dir = start_pos + dir * radius + ττ = τ - st / 2 + st * rand(rng) if (dir + ϵ) > Aupper - rt = total_time + τ + rt = total_time + ττ choice = 1 return (; choice, rt) elseif (dir - ϵ) < Alower - rt = total_time + τ + rt = total_time + ττ choice = 2 return (; choice, rt) else @@ -407,14 +527,17 @@ Returns 2 for the number of choice options n_options(d::DDM) = 2 function simulate(rng::AbstractRNG, model::DDM; Δt = 0.001) - (; ν, α, z) = model - x = α * z + (;ν, α, z, η, sz) = model + + zz = (z - sz/2) + ((z + sz/2) - (z - sz/2)) * rand(rng) + νν = rand(rng, Distributions.Normal(ν, η)) + x = α * zz t = 0.0 evidence = [x] time_steps = [t] while (x < α) && (x > 0) t += Δt - x += ν * Δt + rand(rng, Normal(0.0, 1.0)) * √(Δt) + x += νν * Δt + rand(rng, Normal(0.0, 1.0)) * √(Δt) push!(evidence, x) push!(time_steps, t) end diff --git a/test/ddm_tests.jl b/test/ddm_tests.jl index 8654af2e..b043297d 100644 --- a/test/ddm_tests.jl +++ b/test/ddm_tests.jl @@ -6,7 +6,7 @@ Random.seed!(654) include("KDE.jl") - dist = DDM(; ν = 1.0, α = 0.8, z = 0.5, τ = 0.3) + dist = DDM(; ν=1.0, α = 0.8, z = 0.5, τ = 0.3, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0) choice, rt = rand(dist, 10^6) rt1 = rt[choice .== 1] p1 = mean(choice .== 1) @@ -31,7 +31,7 @@ include("KDE.jl") Random.seed!(750) - dist = DDM(; ν = 2.0, α = 1.5, z = 0.5, τ = 0.30) + dist = DDM(; ν = 2.0, α = 1.5, z = 0.5, τ = 0.30, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0) choice, rt = rand(dist, 10^6) rt1 = rt[choice .== 1] p1 = mean(choice .== 1) @@ -56,7 +56,7 @@ using Random Random.seed!(7540) - dist = DDM(; ν = 1.0, α = 0.8, z = 0.5, τ = 0.3) + dist = DDM(; ν = 1.0, α = 0.8, z = 0.5, τ = 0.3, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0) choice, rt = rand(dist, 10^5) rt1 = rt[choice .== 1] p1 = mean(choice .== 1) @@ -81,7 +81,7 @@ using Random Random.seed!(2200) - dist = DDM(; ν = 2.0, α = 1.5, z = 0.5, τ = 0.30) + dist = DDM(; ν = 2.0, α = 1.5, z = 0.5, τ = 0.30, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0) choice, rt = rand(dist, 10^5) rt1 = rt[choice .== 1] p1 = mean(choice .== 1) @@ -99,42 +99,42 @@ @test y′ ≈ y rtol = 0.01 end - @safetestset "_exp_pnorm" begin - using SequentialSamplingModels: _exp_pnorm - using Test - true_values = [ - 0.003078896, - 0.021471654, - 0.067667642, - 0.113863630, - 0.132256388, - 0.008369306, - 0.058366006, - 0.183939721, - 0.309513435, - 0.359510135, - 0.022750132, - 0.158655254, - 0.500000000, - 0.841344746, - 0.977249868, - 0.061841270, - 0.431269694, - 1.359140914, - 2.287012135, - 2.656440558, - 0.168102001, - 1.172312572, - 3.694528049, - 6.216743527, - 7.220954098 - ] - cnt = 0 - for v1 ∈ -2:2, v2 ∈ -2:2 - cnt += 1 - @test true_values[cnt] ≈ _exp_pnorm(v1, v2) atol = 1e-5 - end - end + # @safetestset "_exp_pnorm" begin + # using SequentialSamplingModels: _exp_pnorm + # using Test + # true_values = [ + # 0.003078896, + # 0.021471654, + # 0.067667642, + # 0.113863630, + # 0.132256388, + # 0.008369306, + # 0.058366006, + # 0.183939721, + # 0.309513435, + # 0.359510135, + # 0.022750132, + # 0.158655254, + # 0.500000000, + # 0.841344746, + # 0.977249868, + # 0.061841270, + # 0.431269694, + # 1.359140914, + # 2.287012135, + # 2.656440558, + # 0.168102001, + # 1.172312572, + # 3.694528049, + # 6.216743527, + # 7.220954098 + # ] + # cnt = 0 + # for v1 ∈ -2:2, v2 ∈ -2:2 + # cnt += 1 + # @test true_values[cnt] ≈ _exp_pnorm(v1, v2) atol = 1e-5 + # end + # end # exp_pnorm = function(a, b) # { # r = exp(a) * pnorm(b) @@ -151,17 +151,17 @@ # } # } - @safetestset "K_large" begin - using SequentialSamplingModels - using SequentialSamplingModels: _K_large - using Test + # @safetestset "K_large" begin + # using SequentialSamplingModels + # using SequentialSamplingModels: _K_large + # using Test - test_val1 = _K_large(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 0.75; ϵ = 1e-10) - @test test_val1 ≈ 3.0 + # test_val1 = _K_large(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 0.75; ϵ = 1e-10) + # @test test_val1 ≈ 3.0 - test_val2 = _K_large(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 0.501; ϵ = 1e-10) - @test test_val2 ≈ 36 - end + # test_val2 = _K_large(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 0.501; ϵ = 1e-10) + # @test test_val2 ≈ 36 + # end # K_large = function(t, v, a, w, epsilon) # { # sqrtL1 = sqrt(1/t) * a/pi @@ -171,22 +171,22 @@ # K_large(.25, 1.5, .5, .25, 1e-10) # K_large(.001, 1.5, .5, .25, 1e-10) - @safetestset "K_small" begin - using SequentialSamplingModels - using SequentialSamplingModels: _K_small - using Test + # @safetestset "K_small" begin + # using SequentialSamplingModels + # using SequentialSamplingModels: _K_small + # using Test - # function(t, v, a, w, epsilon) - test_val1 = _K_small(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 0.75; ϵ = 1e-10) - @test test_val1 ≈ 16 + # # function(t, v, a, w, epsilon) + # test_val1 = _K_small(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 0.75; ϵ = 1e-10) + # @test test_val1 ≈ 16 - test_val2 = _K_small(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 0.501; ϵ = 1e-10) - @test test_val2 ≈ 16 + # test_val2 = _K_small(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 0.501; ϵ = 1e-10) + # @test test_val2 ≈ 16 - test_val3 = - _K_small(DDM(; ν = 0.50, α = 2.50, z = 0.50, τ = 0.30), 0.400; ϵ = 1e-10) - @test test_val3 ≈ 10 - end + # test_val3 = + # _K_small(DDM(; ν = 0.50, α = 2.50, z = 0.50, τ = 0.30, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 0.400; ϵ = 1e-10) + # @test test_val3 ≈ 10 + # end # K_small = function(t, v, a, w, epsilon) # { # if(abs(v) < sqrt(.Machine$double.eps)) # zero drift case @@ -203,35 +203,35 @@ # K_small(.001, 1.5, .5, .25, 1e-10) # K_small(.10, .50, 2.5, .50, 1e-10) - @safetestset "_P_upper" begin - using SequentialSamplingModels - using SequentialSamplingModels: _P_upper - using Test + # @safetestset "_P_upper" begin + # using SequentialSamplingModels + # using SequentialSamplingModels: _P_upper + # using Test - test_val1 = _P_upper(1.0, 0.5, 0.5) - @test test_val1 ≈ 0.3775407 atol = 1e-5 + # test_val1 = _P_upper(1.0, 0.5, 0.5) + # @test test_val1 ≈ 0.3775407 atol = 1e-5 - test_val2 = _P_upper(-1.0, 0.75, 0.25) - @test test_val2 ≈ 0.8693188 atol = 1e-5 + # test_val2 = _P_upper(-1.0, 0.75, 0.25) + # @test test_val2 ≈ 0.8693188 atol = 1e-5 - test_val3 = _P_upper(-1e10, 0.75, 0.25) - @test test_val3 ≈ 1 + # test_val3 = _P_upper(-1e10, 0.75, 0.25) + # @test test_val3 ≈ 1 - test_val4 = _P_upper(eps(), 0.75, 0.25) - @test test_val4 ≈ 0.75 - end + # test_val4 = _P_upper(eps(), 0.75, 0.25) + # @test test_val4 ≈ 0.75 + # end - @safetestset "_Fs_lower" begin - using SequentialSamplingModels - using SequentialSamplingModels: _Fs_lower - using Test + # @safetestset "_Fs_lower" begin + # using SequentialSamplingModels + # using SequentialSamplingModels: _Fs_lower + # using Test - test_val1 = _Fs_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 10, 0.75) - @test test_val1 ≈ 0.5955567 atol = 1e-5 + # test_val1 = _Fs_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 10, 0.75) + # @test test_val1 ≈ 0.5955567 atol = 1e-5 - test_val2 = _Fs_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 10, 0.501) - @test test_val2 ≈ 6.393096e-05 atol = 1e-8 - end + # test_val2 = _Fs_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 10, 0.501) + # @test test_val2 ≈ 6.393096e-05 atol = 1e-8 + # end # Fs_lower = function(t, v, a, w, K) # { @@ -253,17 +253,17 @@ # Fs_lower(.25, 1.5, .5, .25, 10) # Fs_lower(.001, 1.5, .5, .25, 10) - @safetestset "_Fl_lower" begin - using SequentialSamplingModels - using SequentialSamplingModels: _Fl_lower - using Test + # @safetestset "_Fl_lower" begin + # using SequentialSamplingModels + # using SequentialSamplingModels: _Fl_lower + # using Test - test_val1 = _Fl_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 10, 0.75) - @test test_val1 ≈ 0.5955567 atol = 1e-5 + # test_val1 = _Fl_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 10, 0.75) + # @test test_val1 ≈ 0.5955567 atol = 1e-5 - test_val2 = _Fl_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50), 10, 0.501) - @test test_val2 ≈ 0.001206451 atol = 1e-8 - end + # test_val2 = _Fl_lower(DDM(; ν = 1.5, α = 0.50, z = 0.25, τ = 0.50, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 10, 0.501) + # @test test_val2 ≈ 0.001206451 atol = 1e-8 + # end # Fl_lower = function(t, v, a, w, K) # { # F = numeric(length(t)) @@ -279,16 +279,16 @@ using SequentialSamplingModels using Test # tested against rtdists - test_val1 = pdf(DDM(; ν = 2.0, α = 1.0, z = 0.5, τ = 0.3), 1, 0.5) + test_val1 = pdf(DDM(; ν = 2.0, α = 1.0, z = 0.5, τ = 0.3, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 1, 0.5) @test test_val1 ≈ 2.131129 atol = 1e-5 - test_val2 = pdf(DDM(; ν = 2.0, α = 1.0, z = 0.5, τ = 0.3), 2, 0.5) + test_val2 = pdf(DDM(; ν = 2.0, α = 1.0, z = 0.5, τ = 0.3, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 2, 0.5) @test test_val2 ≈ 0.2884169 atol = 1e-5 - test_val3 = pdf(DDM(; ν = 0.8, α = 0.5, z = 0.3, τ = 0.2), 1, 0.35) + test_val3 = pdf(DDM(; ν = 0.8, α = 0.5, z = 0.3, τ = 0.2, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 1, 0.35) @test test_val3 ≈ 0.6635714 atol = 1e-5 - test_val4 = pdf(DDM(; ν = 0.8, α = 0.5, z = 0.3, τ = 0.2), 2, 0.35) + test_val4 = pdf(DDM(; ν = 0.8, α = 0.5, z = 0.3, τ = 0.2, η = 0.00, sz = 0.00, st = 0.00, σ = 1.0), 2, 0.35) @test test_val4 ≈ 0.4450956 atol = 1e-5 end From 5bc967fdfd7fc9713b2f18e5d88fefd4251ce08a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?= <61021880+kiante-fernandez@users.noreply.github.com> Date: Sat, 13 Jul 2024 13:06:28 +0200 Subject: [PATCH 2/3] tracking the use normalized time for debugging _pdf. --- src/DDM.jl | 49 +++++++++++++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/src/DDM.jl b/src/DDM.jl index 9260654f..443a079b 100644 --- a/src/DDM.jl +++ b/src/DDM.jl @@ -91,12 +91,12 @@ function _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, simps_err::Real = 1e-6) wh (ν, α, z, τ, η, sz, st, σ) = params(d) # Check if parameters are valid - if (z < 0 || z > 1 || α <= 0 || τ < 0 || st < 0 || η < 0 || sz < 0 || sz > 1 || - (rt - (τ - st/2)) < 0 || (z + sz/2) > 1 || (z - sz/2) < 0 || (τ - st/2) < 0) - return zero(T) - end + # if (z < 0 || z > 1 || α <= 0 || τ < 0 || st < 0 || η < 0 || sz < 0 || sz > 1 || + # (rt - (τ - st/2)) < 0 || (z + sz/2) > 1 || (z - sz/2) < 0 || (τ - st/2) < 0) + # return zero(T) + # end - rt = abs(rt) + rt = abs.(rt) if st < 1e-6 st = 0 @@ -109,7 +109,7 @@ function _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, simps_err::Real = 1e-6) wh if st == 0 # η=0, sz=0, st=0 return _pdf_sv(DDM(ν, α, z, τ, η, 0, 0, σ), rt - τ; ϵ) else # η=0, sz=0, st≠0 - f = y -> _pdf_sv(DDM(ν, α, z, y[1], η, 0, 0, σ), rt - y[1]; ϵ) + f = y -> _pdf_sv(DDM(ν, α, z, y[1], η, 0, 0, σ), rt .- y[1]; ϵ) result, _ = hcubature(f, (τ-st/2,), (τ+st/2,); rtol=simps_err, atol=simps_err) return result / st # Normalize by the integration range end @@ -119,7 +119,7 @@ function _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, simps_err::Real = 1e-6) wh result, _ = hcubature(f, (z-sz/2,), (z+sz/2,); rtol=simps_err, atol=simps_err) return result / sz else # η=0, sz≠0, st≠0 - f = y -> _pdf_sv(DDM(ν, α, y[1], y[2], η, 0, 0, σ), rt - y[2]; ϵ) + f = y -> _pdf_sv(DDM(ν, α, y[1], y[2], η, 0, 0, σ), rt .- y[2]; ϵ) result, _ = hcubature(f, (z-sz/2, τ-st/2), (z+sz/2, τ+st/2); rtol=simps_err, atol=simps_err) return result / (sz * st) end @@ -141,24 +141,30 @@ for variability in drift-rate. - Returns the computed PDF value. """ -function _pdf_sv(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} +function _pdf_sv(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T<:Real} (ν, α, z, τ, η, σ) = params(d) - if t <= 0 - return zero(T) + if !(t > zero(t)) + return zero(promote_type(T, typeof(t))) end - if η == 0 - return _pdf(DDM(ν, α, z, τ, 0, 0, 0, σ), t; ϵ) + if iszero(η) + # XX = t / (α^2) #use normalized time + return _pdf(DDM(ν, α, z, τ, zero(η), zero(η), zero(η), σ), t; ϵ) end + p = _pdf(DDM(ν, α, z, τ, zero(η), zero(η), zero(η), σ), t; ϵ) #get f(t|0,1,w) + XX = t / (α^2) #use normalized time - p = _pdf(DDM(ν, α, z, τ, 0, 0, 0, σ), XX; ϵ) #get f(t|0,1,w) # convert to f(t|v,a,w) return exp(log(p) + ((α*z*η)^2 - 2*α*ν*z - (ν^2)*XX)/(2*(η^2)*XX+2)) / sqrt((η^2)*XX+1) / (α^2) end +# function _pdf_sv(d::DDM{T}, t::AbstractVector; ϵ::Real = 1.0e-12) where {T<:Real} +# return [_pdf_sv(d, ti; ϵ=ϵ) for ti in t] +# end + ################################################################################ # Converted from WienerDiffusionModel.jl repository orginally by Tobias Alfers# # See https://github.com/t-alfers/WienerDiffusionModel.jl # @@ -196,17 +202,19 @@ function _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} if t <= 0 return zero(T) end - u = t +# u = t + u = t / (α^2) #use normalized time + - K_s = 2.0 - K_l = 1 / (π * sqrt(u)) + K_s = 2.0 #minimal kappa + K_l = 1 / (π * sqrt(u)) #boundary condition # number of terms needed for large-time expansion - if (π * u * ϵ) < 1 + if (π * u * ϵ) < 1 # if error threshold is set low enough K_l = max(sqrt((-2 * log(π * u * ϵ)) / (π^2 * u)), K_l) end # number of terms needed for small-time expansion if (2 * sqrt(2 * π * u) * ϵ) < 1 - K_s = max(2 + sqrt(-2u * log(2ϵ * sqrt(2 * π * u))), sqrt(u) + 1) + K_s = max(2 + sqrt(-2 * u * log(2 * sqrt(2 * π * u) * ϵ)), sqrt(u) + 1) end p = exp((-α * z * ν) - (0.5 * (ν^2) * (t))) / (α^2) @@ -225,7 +233,7 @@ function _small_time_pdf(u::T, z::T, K::Int) where {T <: Real} k_series = (-floor(Int, 0.5 * (K - 1))):ceil(Int, 0.5 * (K - 1)) for k in k_series - inf_sum += ((2k + z) * exp(-((2k + z)^2 / (2u)))) + inf_sum += ((z + 2 * k) * exp(-((z + 2 * k)^2/2/u))) end return inf_sum / sqrt(2π * u^3) @@ -236,7 +244,8 @@ function _large_time_pdf(u::T, z::T, K::Int) where {T <: Real} inf_sum = zero(T) for k = 1:K - inf_sum += (k * exp(-0.5 * (k^2 * π^2 * u)) * sin(k * π * z)) + # inf_sum += (k * exp(-0.5 * (k^2 * π^2 * u)) * sin(k * π * z)) + inf_sum += (k * exp(-(k^2)) * (π^2) * u / 2) * sin(k * π * z) end return π * inf_sum From 80791c22871b99d0322ae0bef8e066c9f4142236 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kiant=C3=A9=20Fernandez?= <61021880+kiante-fernandez@users.noreply.github.com> Date: Sun, 14 Jul 2024 22:17:22 +0200 Subject: [PATCH 3/3] validated REwritten wfpt --- src/DDM.jl | 117 ++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 81 insertions(+), 36 deletions(-) diff --git a/src/DDM.jl b/src/DDM.jl index 443a079b..c476b03a 100644 --- a/src/DDM.jl +++ b/src/DDM.jl @@ -202,54 +202,99 @@ function _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} if t <= 0 return zero(T) end -# u = t - u = t / (α^2) #use normalized time - + u = t / (α^2) # use normalized time + w = z / α # convert to relative start point + + # calculate number of terms needed for large t + if π * u * ϵ < 1 # if error threshold is set low enough + kl = sqrt(-2 * log(π * u * ϵ) / (π^2 * u)) # bound + kl = max(kl, 1 / (π * sqrt(u))) # ensure boundary conditions met + else # if error threshold set too high + kl = 1 / (π * sqrt(u)) # set to boundary condition + end - K_s = 2.0 #minimal kappa - K_l = 1 / (π * sqrt(u)) #boundary condition - # number of terms needed for large-time expansion - if (π * u * ϵ) < 1 # if error threshold is set low enough - K_l = max(sqrt((-2 * log(π * u * ϵ)) / (π^2 * u)), K_l) + # calculate number of terms needed for small t + if 2 * sqrt(2π * u) * ϵ < 1 # if error threshold is set low enough + ks = 2 + sqrt(-2 * u * log(2 * sqrt(2π * u) * ϵ)) # bound + ks = max(ks, sqrt(u) + 1) # ensure boundary conditions are met + else # if error threshold was set too high + ks = 2 # minimal kappa for that case end - # number of terms needed for small-time expansion - if (2 * sqrt(2 * π * u) * ϵ) < 1 - K_s = max(2 + sqrt(-2 * u * log(2 * sqrt(2 * π * u) * ϵ)), sqrt(u) + 1) + + # compute f(u|0,1,w) + p = 0.0 # initialize density + if ks < kl # if small t is better + K = ceil(Int, ks) # round to smallest integer meeting error + for k in -floor(Int, (K-1)/2):ceil(Int, (K-1)/2) # loop over k + p += (w + 2*k) * exp(-((w + 2*k)^2)/2/u) # increment sum + end + p = p / sqrt(2*π*u^3) # add constant term + else # if large t is better + K = ceil(Int, kl) + for k in 1:K + p += k * exp(-(k^2)*(π^2)*u/2) * sin(k*π*w) # increment sum + end + p *= pi # add constant term end + # convert to f(t|ν,α,w) + p *= exp(-ν*α*w - (ν^2)*t/2) / (α^2) + return p +end +# function _pdf(d::DDM{T}, t::Real; ϵ::Real = 1.0e-12) where {T <: Real} +# (ν, α, z) = params(d) - p = exp((-α * z * ν) - (0.5 * (ν^2) * (t))) / (α^2) +# if t <= 0 +# return zero(T) +# end +# # u = t +# u = t / (α^2) #use normalized time +# w = z / α # convert to relative start point - # decision rule for infinite sum algorithm - if K_s < K_l - return p * _small_time_pdf(u, z, ceil(Int, K_s)) - end - return p * _large_time_pdf(u, z, ceil(Int, K_l)) -end +# K_s = 2.0 #minimal kappa +# K_l = 1 / (π * sqrt(u)) #boundary condition +# # number of terms needed for large-time expansion +# if (π * u * ϵ) < 1 # if error threshold is set low enough +# K_l = max(sqrt((-2 * log(π * u * ϵ)) / (π^2 * u)), K_l) +# end +# # number of terms needed for small-time expansion +# if (2 * sqrt(2 * π * u) * ϵ) < 1 +# K_s = max(2 + sqrt(-2 * u * log(2 * sqrt(2 * π * u) * ϵ)), sqrt(u) + 1) +# end -# small-time expansion -function _small_time_pdf(u::T, z::T, K::Int) where {T <: Real} - inf_sum = zero(T) +# p = exp((-α * w * ν) - (0.5 * (ν^2) * (t))) / (α^2) - k_series = (-floor(Int, 0.5 * (K - 1))):ceil(Int, 0.5 * (K - 1)) - for k in k_series - inf_sum += ((z + 2 * k) * exp(-((z + 2 * k)^2/2/u))) - end +# # decision rule for infinite sum algorithm +# if K_s < K_l +# return p * _small_time_pdf(u, w, ceil(Int, K_s)) +# end - return inf_sum / sqrt(2π * u^3) -end +# return p * _large_time_pdf(u, w, ceil(Int, K_l)) +# end -# large-time expansion -function _large_time_pdf(u::T, z::T, K::Int) where {T <: Real} - inf_sum = zero(T) +# # small-time expansion +# function _small_time_pdf(u::T, z::T, K::Int) where {T <: Real} +# inf_sum = zero(T) - for k = 1:K - # inf_sum += (k * exp(-0.5 * (k^2 * π^2 * u)) * sin(k * π * z)) - inf_sum += (k * exp(-(k^2)) * (π^2) * u / 2) * sin(k * π * z) - end +# k_series = (-floor(Int, 0.5 * (K - 1))):ceil(Int, 0.5 * (K - 1)) +# for k in k_series +# inf_sum += ((z + 2 * k) * exp(-((z + 2 * k)^2/2/u))) +# end - return π * inf_sum -end +# return inf_sum / sqrt(2π * u^3) +# end + +# # large-time expansion +# function _large_time_pdf(u::T, z::T, K::Int) where {T <: Real} +# inf_sum = zero(T) + +# for k = 1:K +# # inf_sum += (k * exp(-0.5 * (k^2 * π^2 * u)) * sin(k * π * z)) +# inf_sum += (k * exp(-(k^2)) * (π^2) * u / 2) * sin(k * π * z) +# end + +# return π * inf_sum +# end 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; ϵ))