diff --git a/.gitignore b/.gitignore index b7bdae0b..474a45ce 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,3 @@ Manifest.toml /temp -/benchmark/*.json -/benchmark/*.md \ No newline at end of file +/.vscode \ No newline at end of file diff --git a/docs/.gitignore b/docs/.gitignore index 93584992..1fa74a1a 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -1,4 +1,5 @@ build/ site/ docs/build/ -Manifest.toml \ No newline at end of file +Manifest.toml +settings.json \ No newline at end of file diff --git a/docs/src/DDM.md b/docs/src/DDM.md index 2e29adb0..effa2455 100644 --- a/docs/src/DDM.md +++ b/docs/src/DDM.md @@ -1,9 +1,13 @@ # 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 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). +However, in the Ratcliff Diffusion Decision Model, we also include across-trial variability parameters. These parameters were developed to explain specific discrepancies between the DDM and experimental data (Anderson, 1960; Laming, 1968; Blurton et al., 2017). The data exhibited a difference in mean RT between correct and error responses that could not be captured by the DDM. As a result, two parameters for across-trial variability were introduced to explain this difference: across-trial variability in the starting point (sz) to explain fast errors (Laming, 1968), and across-trial variability in drift rate (η) to explain slow errors (Ratcliff, 1978; Ratcliff and Rouder, 1998). Additionally, the DDM also showed a sharper rise in the leading edge of the response time distribution than observed in the data. To capture this leading edge effect, across-trial variability in non-decision time (st) was introduced. + +Previous work has validated predictions of these across-trial variability parameters (Wagenmakers et al., 2009). When compared to the DDM, the Ratcliff DDM improves the fit to the data. Researchers now often assume that the core parameters of sequential sampling models, such as drift rates, non-decision times, and starting points vary between trials. + 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 @@ -33,8 +37,11 @@ In the code below, we will define parameters for the DDM and create a model obje 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 +Across-trial-variability of drift rate. Standard deviation of a normal distribution with mean v describing the distribution of actual drift rates from specific trials. Values different from 0 can predict slow errors. Typical range: 0 < η < 2. Default is 0. + ```@example DDM ν=1.0 +η = 0.16 ``` ### Boundary Separation @@ -49,16 +56,22 @@ The amount of information that is considered for a decision. Large values indica The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 +Across-trial-variability of non-decisional components. Range of a uniform distribution with mean τ + st/2 describing the distribution of actual τ values across trials. Accounts for response times below t0. Reduces skew of predicted RT distributions. Typical range: 0 < τ < 0.2. Default is 0. + ```@example DDM τ = 0.30 +st = 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). +Across-trial-variability of starting point. Range of a uniform distribution with mean z describing the distribution of actual starting points from specific trials. Values different from 0 can predict fast errors. Typical range: 0 < sz < 0.5. Default is 0. + ```@example DDM z = 0.50 +sz = 0.05 ``` ### DDM Constructor @@ -66,7 +79,7 @@ z = 0.50 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) +dist = DDM(ν, α, τ, z, η, sz, st, σ) ``` ## Simulate Model @@ -108,6 +121,10 @@ plot!(dist; t_range=range(.301, 1, length=100)) # References +Blurton, S. P., Kesselmeier, M., & Gondan, M. (2017). The first-passage time distribution for the diffusion model with variable drift. Journal of Mathematical Psychology, 76, 7–12. https://doi.org/10.1016/j.jmp.2016.11.003 + +Laming, D. R. J. (1968). Information theory of choice-reaction times. Academic Press. + 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 @@ -117,3 +134,5 @@ Ratcliff, R., & Rouder, J. N. (1998). Modeling Response Times for Two-Choice Dec 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 + +Wagenmakers, E.-J. (2009). Methodological and empirical developments for the Ratcliff diffusion model of response times and accuracy. European Journal of Cognitive Psychology, 21(5), 641-671. diff --git a/docs/src/Ratcliff_DDM.md b/docs/src/Ratcliff_DDM.md deleted file mode 100644 index f6cd3ff4..00000000 --- a/docs/src/Ratcliff_DDM.md +++ /dev/null @@ -1,141 +0,0 @@ -# Ratcliff Diffusion Model - -The Ratcliff Diffusion Model (Ratcliff DDM; Ratcliff et al., 2016) is similar to the DDM. Like the DDM, the model 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). 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). - -However, the model differs in the inclusion of across-trial variability parameters. These parameters were developed to explain specific discrepancies between the DDM and experimental data (Anderson, 1960; Laming, 1968; Blurton et al., 2017). The data exhibited a difference in mean RT between correct and error responses that could not be captured by the DDM. As a result, two parameters for across-trial variability were introduced to explain this difference: across-trial variability in the starting point to explain fast errors (Laming, 1968), and across-trial variability in drift rate to explain slow errors (Ratcliff, 1978; Ratcliff and Rouder, 1998). Additionally, the DDM also showed a sharper rise in the leading edge of the response time distribution than observed in the data. To capture this leading edge effect, across-trial variability in non-decision time was introduced. - -Previous work has validated predictions of these across-trial variability parameters (Wagenmakers et al., 2009). When compared to the DDM, the Ratcliff DDM improves the fit to the data. Researchers now often assume that the core parameters of sequential sampling models, such as drift rates, non-decision times, and starting points vary between trials. - -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 RatcliffDDM -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 - -Across-trial-variability of drift rate. Standard deviation of a normal distribution with mean v describing the distribution of actual drift rates from specific trials. Values different from 0 can predict slow errors. Typical range: 0 < η < 2. Default is 0. - -```@example RatcliffDDM -ν=1.0 -η = 0.16 -``` - -### Boundary Separation - -The amount of information that is considered for a decision. Large values indicates response caution. Typical range: 0.5 < α < 2 - - -```@example RatcliffDDM -α = 0.80 -``` - -### Non-Decision Time - -The duration for a non-decisional processes (encoding and response execution). Typical range: 0.1 < τ < 0.5 - -Across-trial-variability of non-decisional components. Range of a uniform distribution with mean τ + st/2 describing the distribution of actual τ values across trials. Accounts for response times below t0. Reduces skew of predicted RT distributions. Typical range: 0 < τ < 0.2. Default is 0. - -```@example RatcliffDDM -τ = 0.30 -st = 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). - -Across-trial-variability of starting point. Range of a uniform distribution with mean z describing the distribution of actual starting points from specific trials. Values different from 0 can predict fast errors. Typical range: 0 < sz < 0.5. Default is 0. - -```@example RatcliffDDM -z = 0.25 -sz = 0.05 -``` - -### Ratcliff Diffusion Model Constructor - -Now that values have been assigned to the parameters, we will pass them to `RatcliffDDM` to generate the model object. - -```@example RatcliffDDM -dist = DDM(ν, α, τ, z) -``` - -## Simulate Model - -Now that the model is defined, we will generate $10,000$ choices and reaction times using `rand`. - - ```@example RatcliffDDM - choices,rts = rand(dist, 10_000) -``` - -## Compute PDF -The PDF for each observation can be computed as follows: - - ```@example RatcliffDDM -pdf.(dist, choices, rts) -``` -## Compute Log PDF -Similarly, the log PDF for each observation can be computed as follows: - - ```@example RatcliffDDM -logpdf.(dist, choices, rts) -``` - -## Plot Simulation -The code below overlays the PDF on reaction time histograms for each option. - - ```@example RatcliffDDM -# 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 - -Wagenmakers, E.-J. (2009). Methodological and empirical developments for the Ratcliff diffusion model of response times and accuracy. European Journal of Cognitive Psychology, 21(5), 641-671. - - diff --git a/src/DDM.jl b/src/DDM.jl index 973ea80b..7648f9e0 100644 --- a/src/DDM.jl +++ b/src/DDM.jl @@ -1,81 +1,189 @@ """ - DDM{T<:Real} <: SSM2D + DDM{T<:Real} <: AbstractDDM -Model object for the standard Drift Diffusion Model. + 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. # Constructors - DDM(ν, α, τ, z) + DDM(ν, α, τ, z, η, sz, st, σ) - DDM(; ν = 1.0, - α = 0.8, - τ = 0.3 - z = 0.25) - + DDM(; ν = 1.00, + α = 0.80, + τ = 0.30, + z = 0.25, + η = 0.16, + sz = 0.05, + st = 0.10, + σ = 1.0 + ) + # Example -```julia +````julia using SequentialSamplingModels -dist = DDM(ν = 1.0, α = 0.8, τ = 0.3, z = 0.25) +dist = DDM(ν = 1.0,α = 0.80,τ = 0.30,z = 0.25,η = 0.16,sz = 0.05,st = .10,σ = 1) choice,rt = rand(dist, 10) like = pdf.(dist, choice, rt) loglike = logpdf.(dist, choice, rt) -``` - -# References +```` +# 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} <: SSM2D +mutable struct DDM{T<:Real} <: AbstractDDM ν::T α::T τ::T z::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) - (d.ν, d.α, d.τ, d.z) + (d.ν, d.α, d.τ, d.z,d.η, d.sz, d.st, d.σ) end function DDM(; ν = 1.00, α = 0.80, τ = 0.30, - z = 0.50) - return DDM(ν, α, τ, z) + z = 0.25, + η = 0.16, + sz = 0.05, + st = 0.10, + σ = 1.0) + return DDM(ν, α, τ, z, η, sz, st, σ) end -################################################################################ -# Converted from WienerDiffusionModel.jl repository orginally by Tobias Alfers# -# See https://github.com/t-alfers/WienerDiffusionModel.jl # -################################################################################ - -##################################### -# Probability density function # -# Navarro & Fuss (2009) # -# 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; ϵ) + (ν, α, τ, z, η, sz, st, σ) = params(d) + return _pdf_Full(DDM(-ν, α, τ, 1-z, η, sz, st, σ), rt; ϵ) end - return _pdf(d, rt; ϵ) + return _pdf_Full(d, rt; ϵ) end -# probability density function over the lower boundary + +""" + _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, n_st::Int=2, n_sz::Int=2) 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. +- `n_st::Int`: specifies the number of subintervals in the Simpson's rule for the integration associated with non-decision time variability. Default is 2. +- `n_sz::Int`: specifies the number of subintervals in the Simpson's rule for the integration associated with starting point variability. Default is 2. + +""" +function _pdf_Full(d::DDM{T}, rt; ϵ::Real = 1.0e-12, n_st::Int=2, n_sz::Int=2) where {T<:Real} + + (ν, α, τ, z, η, sz, st, σ) = params(d) + + if τ ≥ rt + return T(NaN) + end + + if st < 1.0e-3 + st = 0 + end + if sz < 1.0e-3 + sz = 0 + end + + if sz==0 + if st==0 #sv=0,sz=0,st=0 + return _pdf_sv(d, rt; ϵ) + else #sv=0,sz=0,st=$ + return _simpson_1D(rt, ν, η, α, z, τ, ϵ, z, z, 0, τ-st/2., τ+st/2., n_st) + end + else #sz=$ + if st==0 #sv=0,sz=$,st=0 + return _simpson_1D(rt, ν, η, α, z, τ, ϵ, z-sz/2., z+sz/2., n_sz, τ, τ , 0) + else #sv=0,sz=$,st=$ + return _simpson_2D(rt, ν, η, α, z, τ, ϵ, z-sz/2., z+sz/2., n_sz, τ-st/2., τ+st/2., n_st) + end + end +end + +""" + _pdf_sv(d::DDM, rt; ϵ::Real = 1.0e-12) + +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, rt; ϵ::Real = 1.0e-12) + (ν, α, τ, z, η, sz, st, σ) = params(d) + + if η == 0 + return _pdf(DDM(ν, α, τ, z, 0, 0, 0, σ), rt; ϵ) + end + + # α *= 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 + +""" + _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. + +# 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) + (ν, α, τ, z, η, sz, st, σ) = params(d) + if τ ≥ t return T(NaN) end @@ -125,230 +233,488 @@ function _large_time_pdf(u::T, z::T, K::Int) where {T<:Real} 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; ϵ)) +""" + _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) -function logpdf(d::DDM, data::T) where {T<:NamedTuple} - return sum(logpdf.(d, data...)) -end +Calculate the 1-dimensional Simpson's numerical integration for a drift diffusion model with given parameters. This function is used for integrating over either the starting point or the non-decision time. -function logpdf(dist::DDM, data::Array{<:Tuple,1}) - LL = 0.0 - for d in data - LL += logpdf(dist, d...) - end - return LL -end +# Arguments +- `x::Real`: Reaction time for which the probability density function is being computed. +- `ν::Real`, `η::Real`, `α::Real`, `z::Real`, `τ::Real`: Parameters of the Ratcliff Drift Diffusion Model. +- `ϵ::Real`: A small constant to prevent divide by zero errors. +- `lb_z::Real`, `ub_z::Real`: Lower and upper bounds for z (starting point). +- `n_sz::Int`: Specifies the number of subintervals for Simpson's rule in the z dimension. +- `lb_t::Real`, `ub_t::Real`: Lower and upper bounds for t (non-decision time). +- `n_st::Int`: Specifies the number of subintervals for Simpson's rule in the t dimension. -logpdf(d::DDM, data::Tuple) = logpdf(d, data...) +# Returns +- Returns the Simpson's numerical integration of the PDF of the Ratcliff Drift Diffusion Model at the given reaction time `x`, over the specified bounds and subintervals. -######################################### -# Cumulative density function # -# 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 +# Note +- If `n_st` is 0, the function integrates over z (starting point). If `n_sz` is 0, the function integrates over t (non-decision time). - return _cdf(d, rt; ϵ) -end +# References -# 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) +https://en.wikipedia.org/wiki/Simpson%27s_rule + +""" +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) + + if n_st == 0 #integration over z + hz = (ub_z-lb_z)/n + ht = 0 + lb_t = τ + ub_t = τ + else #integration over t + hz = 0 + ht = (ub_t-lb_t)/n + lb_z = z + ub_z = z end - K_l = _K_large(d, t; ϵ) - K_s = _K_small(d, t; ϵ) + S = _pdf_sv(DDM(ν, α, lb_t, lb_z, η, 0, 0, 1), x; ϵ) + + y = 0 + z_tag = 0 + t_tag = 0 + + for i in 1:n + z_tag = lb_z + hz * i + t_tag = lb_t + ht * i + + y = _pdf_sv(DDM(ν, α, t_tag, z_tag, η, 0, 0, 1), x; ϵ) + + if isodd(i) + S += 4 * y + else + S += 2 * y + end - if K_l < 10*K_s - return _Fl_lower(d, K_l, t) end - return _Fs_lower(d, K_s, t) -end + + 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) -# 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 +""" + _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) - sqt = sqrt(t-τ) +Calculate the 2-dimensional Simpson's numerical integration for a drift diffusion model with given parameters. This function is used for integrating over both the starting point and the non-decision time. - S1 = zero(T) - S2 = zero(T) - K_series = K:-1:1 +# Arguments +- `x::Real`: Reaction time for which the probability density function is being computed. +- `ν::Real`, `η::Real`, `α::Real`, `z::Real`, `τ::Real`: Parameters of the Ratcliff Drift Diffusion Model. +- `ϵ::Real`: A small constant to prevent divide by zero errors. +- `lb_z::Real`, `ub_z::Real`: Lower and upper bounds for z (starting point). +- `n_sz::Int`: Specifies the number of subintervals for Simpson's rule in the z dimension. +- `lb_t::Real`, `ub_t::Real`: Lower and upper bounds for t (non-decision time). +- `n_st::Int`: Specifies the number of subintervals for Simpson's rule in the t dimension. - 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)) +# Returns +- Returns the Simpson's numerical integration of the PDF of the Ratcliff Drift Diffusion Model at the given reaction time `x`, over the specified bounds and subintervals in both dimensions. - S2 += (_exp_pnorm(-2*ν*α*k, sign(ν)*(2*α*k-α*z-ν*(t-τ))/sqt) - - _exp_pnorm(2*ν*α*k-2*ν*α*z, -sign(ν)*(2*α*k-α*z+ν*(t-τ))/sqt)) - end +# Note +- The function calls `_simpson_1D` to perform the 1-dimensional integrations over each dimension. - return _P_upper(ν, α, z) + sign(ν) * ((cdf(Normal(), -sign(ν) * (α*z+ν*(t-τ))/sqt) - - _exp_pnorm(-2*ν*α*z, sign(ν) * (α*z-ν*(t-τ)) / sqt)) + S1 + S2) -end +""" +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) -# 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 + ht = (ub_t-lb_t)/n_st + S = _simpson_1D(x, ν, η, α, z, τ, ϵ, lb_z, ub_z, n_sz, 0, 0, 0) -# 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 + t_tag = 0 + y = 0 + for i_t in 1:n_st + t_tag = lb_t + ht * i_t + y = _simpson_1D(x, ν, η, α, z, t_tag, ϵ, lb_z, ub_z, n_sz, 0, 0, 0) -# 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(π))) + if isodd(i_t) + S += 4 * y + else + S += 2 * y + end end - return r + + 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) + end +logpdf(d::DDM, choice, rt; ϵ::Real = 1.0e-12) = log(pdf(d, choice, rt; ϵ)) + + +# """ +# 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. + +# # 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. + +# # 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/ +# """ +# function cdf(d::DDM, choice, rt; ϵ::Real = 1e-7) + +# (ν, α, τ, z, η, sz, st, σ) = params(d) + +# size = length(rt) +# y = zeros(Float64, size) +# epsi = 1e-10 + +# #transform parameters +# α = α/10. +# τ = τ +# η = η/10. + epsi +# z = z*(α/10.) +# sz = sz*(α/10.) + epsi +# st = st + epsi +# ν = ν/10. + +# p_boundary = 0.0 + +# for i in 1:size +# y[i] = _cdf(d,rt[i],choice[i],p_boundary; ϵ) +# y[i] = (1 - p_boundary) + rt[i]*y[i] +# end + +# return y +# end +# """ +# _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. + +# # 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. + +# """ +# function _cdf(d::DDM, choice, rt, prob; ϵ::Real = 1e-7) + +# (ν, α, τ, z, η, sz, st, σ) = params(d) +# #Explcit recode of the choice from 2(lower) & 1(upper) to 0(lower) and 1(upper) +# #note we need to make sure this is consistent in the all the relative bound models +# if choice == 2 #lower +# choice = 0 +# elseif choice == 1 #upper +# choice = 1 +# end + +# # Initializing variables +# a2 = α*α +# Z_U = (1-choice)*z+choice*(α-z)+sz/2 +# Z_L = (1-choice)*z+choice*(α-z)-sz/2 +# lower_t = τ-st/2 +# upper_t = 0.0 +# Δ = 1e-29 +# min_rt=0.001 +# v_max = 5000 # maximum number of terms in a partial sum approximating infinite series + +# Fnew = 0.0 +# sum_z=0.0 +# sum_ν=0.0 +# p1 = 0.0 +# p0 = 0.0 +# sum_hist = zeros(3) +# denom = 0.0 +# sifa = 0.0 +# upp = 0.0 +# low = 0.0 +# fact = 0.0 +# exdif = 0.0 +# su = 0.0 +# sl = 0.0 +# zzz = 0.0 +# ser = 0.0 +# nr_ν = 6 +# nr_z = 6 + +# # Defining Gauss-Hermite abscissae and weights for numerical integration +# gk = [-2.3506049736744922818,-1.3358490740136970132,-.43607741192761650950,.43607741192761650950,1.3358490740136970132,2.3506049736744922818] +# w_gh = [.45300099055088421593e-2,.15706732032114842368,.72462959522439207571,.72462959522439207571,.15706732032114842368,.45300099055088421593e-2] +# gz = [-.93246951420315193904,-.66120938646626381541,-.23861918608319693247,.23861918608319712676,.66120938646626459256,.93246951420315160597] +# w_g = [.17132449237917049545,.36076157304813916138,.46791393457269092604,.46791393457269092604,.36076157304813843973,.17132449237917132812] + +# # Adjusting Gauss-Hermite abscissae and weights +# for i=1:nr_ν +# gk[i] = 1.41421356237309505*gk[i]*η+ν +# w_gh[i] = w_gh[i]/1.772453850905515882 +# end +# for i=1:nr_z +# gz[i] = (.5*sz*gz[i])+z +# end + +# # numerical integration +# for i=1:nr_z +# sum_ν=0.0 +# # numerical integration +# for m=1:nr_ν +# if abs(gk[m])>ϵ +# sum_ν+=(exp(-200*gz[i]*gk[m])-1)/(exp(-200*α*gk[m])-1)*w_gh[m] +# else +# sum_ν+=gz[i]/α*w_gh[m] +# end +# end +# sum_z+=sum_ν*w_g[i]/2 +# end +# prob = sum_z + +# if (rt-τ+st/2 > min_rt) # is t larger than lower boundary τ distribution? +# upper_t = min(rt, τ+st/2) +# p1 = prob*(upper_t-lower_t)/st # integrate probability with respect to t +# p0 = (1-prob)*(upper_t-lower_t)/st +# if rt > τ+st/2 # is t larger than upper boundary Ter distribution? +# sum_hist = zeros(3) +# for v in 1:v_max # infinite series +# sum_hist = circshift(sum_hist, 1) +# sum_ν = 0 +# sifa = π*v/α +# for m in 1:nr_ν # numerical integration with respect to xi +# denom = (100*gk[m]*gk[m] + (π*π)*(v*v)/(100*a2)) +# upp = exp((2*choice-1)*Z_U*gk[m]*100 - 3*log(denom) + log(w_gh[m]) - 2*log(100)) +# low = exp((2*choice-1)*Z_L*gk[m]*100 - 3*log(denom) + log(w_gh[m]) - 2*log(100)) +# fact = upp*((2*choice-1)*gk[m]*sin(sifa*Z_U)*100 - sifa*cos(sifa*Z_U)) - +# low*((2*choice-1)*gk[m]*sin(sifa*Z_L)*100 - sifa*cos(sifa*Z_L)) +# exdif = exp((-.5*denom*(rt-upper_t)) + log(1-exp(-.5*denom*(upper_t-lower_t)))) +# sum_ν += fact*exdif +# end +# sum_hist[3] = sum_hist[2] + v*sum_ν +# if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 +# break +# end +# end + +# Fnew = (p0*(1-choice) + p1*choice) - sum_hist[3]*4*π/(a2*sz*st) +# # cumulative distribution function for t and x +# elseif t <= τ+st/2 # is t lower than upper boundary Ter distribution? +# sum_ν = 0 +# for m in 1:nr_ν +# if abs(gk[m]) > ϵ +# sum_z = 0 +# for i in 1:nr_z +# zzz = (α - gz[i])*choice + gz[i]*(1 - choice) +# ser = -((α*a2)/((1 - 2*choice)*gk[m]*π*.01))*sinh(zzz*(1 - 2*x)*gk[m]/.01)/ +# (sinh((1 - 2*choice)*gk[m]*α/.01)^2) + +# (zzz*a2)/((1 - 2*choice)*gk[m]*π*.01)*cosh((α - zzz)*(1 - 2*choice)*gk[m]/.01)/ +# sinh((1 - 2*choice)*gk[m]*α/.01) +# sum_hist = zeros(3) +# for v in 1:v_max +# sum_hist = circshift(sum_hist, 1) +# sifa = π*v/α +# denom = (gk[m]*gk[m]*100 + (π*v)*(π*v)/(a2*100)) +# sum_hist[3] = sum_hist[2] + v*sin(sifa*zzz)*exp(-.5*denom*(rt - lower_t) - 2*log(denom)) +# if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 +# break +# end +# end +# sum_z += .5*w_g[i]*(ser - 4*sum_hist[3])*(π/100)/(a2*st)*exp((2*choice - 1)*zzz*gk[m]*100) +# end +# else +# sum_hist = zeros(3) +# su = -(Z_U*Z_U)/(12*a2) + (Z_U*Z_U*Z_U)/(12*α*a2) - (Z_U*Z_U*Z_U*Z_U)/(48*a2*a2) +# sl = -(Z_L*Z_L)/(12*a2) + (Z_L*Z_L*Z_L)/(12*α*a2) - (Z_L*Z_L*Z_L*Z_L)/(48*a2*a2) +# for v in 1:v_max +# sum_hist = circshift(sum_hist, 1) +# sifa = π*v/α +# denom = (π*v)*(π*v)/(a2*100) +# sum_hist[3] = sum_hist[2] + 1/(π*π*π*π*v*v*v*v)*(cos(sifa*Z_L) - cos(sifa*Z_U))* +# exp(-.5*denom*(rt - lower_t)) +# if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 +# break +# end +# end +# sum_z = 400*a2*α*(sl - su - sum_hist[3])/(st*sz) +# end +# sum_ν += sum_z*w_gh[m] +# end +# Fnew = (p0*(1 - choice) + p1*choice) - sum_ν +# end +# elseif rt - τ + st/2 <= min_rt # is t lower than lower boundary Ter distr? +# Fnew = 0 +# end + +# Fnew = Fnew > Δ ? Fnew : 0 + +# return Fnew + +# end + """ rand(dist::DDM) -Generate a random rt for the Diffusion Decision Model (negative coding) +Generate a random choice and rt for the Diffusion Decision Model # Arguments -- `dist`: model object for the Diffusion Decision Model. +- `dist`: model object for Diffusion Decision Model. + +method simulating the diffusion process: +_rand_rejection uses Tuerlinckx et al., 2001 rejection-based method for the general wiener process +_rand_stochastic uses the stochastic Euler method to simulate the stochastic differential equation + +# References + + Tuerlinckx, F., Maris, E., Ratcliff, R., & De Boeck, P. (2001). + A comparison of four methods for simulating the diffusion process. + Behavior Research Methods, Instruments, & Computers, 33, 443-456. + + Converted from Rhddmjagsutils.R R script by Kianté Fernandez + See also https://github.com/kiante-fernandez/Rhddmjags. """ function rand(rng::AbstractRNG, d::DDM) return _rand_rejection(rng, d) + # return _rand_stochastic(rng, d) 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) - - ϵ = 1.0e-15 - - D = .005 # = 2*σ^2 => 1/200 - zn = (z*α) / 10 # absolute bias! - αn = α / 10 - νn = ν / 10 - - total_time = 0.0 - start_pos = 0.0 - Aupper = αn - zn - Alower = -zn - radius = min(abs(Aupper), abs(Alower)) - λ = 0.0 - F = 0.0 - prob = 0.0 - - while true - if νn == 0 - λ = (0.25D*π^2) / (radius^2) - F = 1.0 - prob = .5 - else - λ = ((0.25*νn^2)/D) + ((0.25*D*π^2) / (radius^2)) - F = (D*π) / (radius*νn) - F = F^2 / (1+F^2) - prob = exp((radius*νn)/D) - prob = prob / (1+prob) +function _rand_rejection(rng::AbstractRNG, d::DDM; N::Int = 1) + (;ν, α, τ, z, η, sz, st, σ) = d + + η = η == 0 ? eps() : η + + # Initialize output vectors + 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 = false + totaltime = 0.0 + startpos = 0.0 + Aupper = α - zz + Alower = -zz + radius = min(abs(Aupper), abs(Alower)) + + 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 * μ) + F = F^2 / (1 + F^2) + # formula p447 in 2001 paper + prob = exp(radius * μ / D) + prob = prob / (1 + prob) + dir_ = 2 * (rand(rng) < prob) - 1 + l = -1.0 + s2 = 0.0 + s1 = 0.0 + while s2 > l + s2 = rand(rng) + s1 = rand(rng) + tnew = 0.0 + told = 0.0 + uu = 0 + while abs(tnew - told) > ϵ || uu == 0 + told = tnew + uu += 1 + tnew = told + (2 * uu + 1) * (-1)^uu * s1^(F * (2 * uu + 1)^2) + # infinite sum in formula (16) in BRMIC,2001 + end + l = 1 + s1^(-F) * tnew + end + # rest of formula (16) + t = abs(log(s1)) / λ + # is the negative of t* in (14) in BRMIC,2001 + totaltime += t + dir_ = startpos + dir_ * radius + ndt = τ - st / 2 + st * rand(rng) + 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=XX,rt=T) +end - r = rand(rng) - dir = r < prob ? 1 : -1 - l = -1.0 - s1 = 0.0 - s2 = 0.0 - - # Tuerlinckx et al. (2001; eq. 16) - while s2 > l - s1 = rand(rng) - s2 = rand(rng) - tnew = 0.0 - tδ = 0.0 - uu = zero(Int) - - while (abs(tδ) > ϵ) || (uu == 0) - uu += 1 - tt = 2*uu + 1 - tδ = tt * (uu % 2 == 0 ? 1 : -1) * (s1 ^ (F * tt^2)) - tnew += tδ - end +function _rand_stochastic(rng::AbstractRNG, d::DDM; N::Int = 1, nsteps::Int=300, step_length::Int=0.01) + (ν, α, τ, z, η, sz, st, σ) = params(d) - l = 1 + (s1^(-F)) * tnew - end + if η == 0 + η = 1e-16 + end - total_time += abs(log(s1)) / λ - dir = start_pos + dir * radius - - if (dir + ϵ) > Aupper - rt = total_time + τ - choice = 1 - return (;choice,rt) - elseif (dir - ϵ) < Alower - rt = total_time + τ - choice = 2 - return (;choice,rt) - else - start_pos = dir - radius = min(abs(Aupper - start_pos), (abs(Alower - start_pos))) + # Initialize output vectors + choice = fill(0, N) + rt = fill(0.0, N) + + for n in 1:N + random_walk = Array{Float64}(undef, nsteps) + start_point = (z - sz/2) + ((z + sz/2) - (z - sz/2)) * rand(rng) + ndt = (τ - st/2) + ((τ + st/2) - (τ - st/2)) * rand(rng) + drift = rand(rng, Distributions.Normal(ν, η)) + random_walk[1] = start_point * α + for s in 2:nsteps + random_walk[s] = random_walk[s-1] + rand(Distributions.Normal(drift * step_length, σ * sqrt(step_length))) + if random_walk[s] >= α + random_walk[s:end] .= α + rts[n] = s * step_length + ndt + choice[n] = 1 + break + elseif random_walk[s] <= 0 + random_walk[s:end] .= 0 + rts[n] = s * step_length + ndt + choice[n] = 2 + break + elseif s == nsteps + rts[n] = NaN + choice[n] = NaN + break + end end - end + end + return (choice=choice,rt=rts) +end + +""" + rand(dist::DDM, n_sim::Int) + +Generate `n_sim` random choice-rt pairs for the Diffusion Decision Model. + +# Arguments +- `dist`: model object for the Ratcliff DDM. +- `n_sim::Int`: the number of simulated rts +""" + +function rand(rng::AbstractRNG, d::DDM, n_sim::Int) + return _rand_rejection(rng, d, N = n_sim) end """ @@ -377,14 +743,19 @@ represent samples of evidence per time step and columns represent different accu - `Δt=.001`: size of time step of decision process in seconds """ function simulate(rng::AbstractRNG, model::DDM; Δt=.001) - (;ν,α,z) = model - x = α * z + (;ν, α, τ, z, η, sz, st, σ) = model + + start_point = (z - sz/2) + ((z + sz/2) - (z - sz/2)) * rand() + drift = rand(rng, Distributions.Normal(ν, η)) + + x = α * start_point 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 += drift * Δt + rand(rng, Distributions.Normal(0.0, 1.0)) * √(Δt) push!(evidence, x) push!(time_steps, t) end @@ -405,4 +776,127 @@ end # choice = (x < α) + 1 # t += rand(Uniform(τ - st / 2, τ + st / 2)) # return choice,t -# end \ No newline at end of file +# end + +####### +# old cdf code for debug checking +######################################### +# Cumulative density function # +# Blurton, Kesselmeier, & Gondan (2012) # +######################################### + +# function cdf(d::DDM, choice, rt; ϵ::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 diff --git a/src/RatcliffDDM.jl b/src/RatcliffDDM.jl deleted file mode 100644 index c8a57b23..00000000 --- a/src/RatcliffDDM.jl +++ /dev/null @@ -1,599 +0,0 @@ -""" - RatcliffDDM{T<:Real} <: SSM2D - - Model object for the Ratcliff Diffusion 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. - -# Constructors - - RatcliffDDM(ν, α, τ, z, η, sz, st, σ) - - RatcliffDDM(; ν = 1.00, - α = 0.80, - τ = 0.30, - z = 0.25, - η = 0.16, - sz = 0.05, - st = 0.10, - σ = 1.0 - ) - -# Example - -````julia -using SequentialSamplingModels -dist = RatcliffDDM(ν = 1.0,α = 0.80,τ = 0.30,z = 0.25,η = 0.16,sz = 0.05,st = .10,σ = 1) -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. -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 RatcliffDDM{T<:Real} <: SSM2D - ν::T - α::T - τ::T - z::T - η::T - sz::T - st::T - σ::T -end - -function RatcliffDDM(ν, α, τ, z, η, sz, st, σ) - return RatcliffDDM(promote(ν, α, τ, z, η, sz, st, σ)...) -end - -function params(d::RatcliffDDM) - (d.ν, d.α, d.τ, d.z,d.η, d.sz, d.st, d.σ) -end - -function RatcliffDDM(; ν = 1.00, - α = 0.80, - τ = 0.30, - z = 0.25, - η = 0.16, - sz = 0.05, - st = 0.10, - σ = 1.0) - return RatcliffDDM(ν, α, τ, z, η, sz, st, σ) -end - -# probability density function over the lower boundary - -# uses analytic integration of the likelihood function for variability in drift-rate -# function pdf_sv(d::RatcliffDDM, choice, rt; ϵ::Real = 1.0e-12) -# if choice == 1 -# (ν, α, τ, z, η, sz, st, σ) = params(d) -# ν = -ν -# z = 1 - z -# return pdf_sv(RatcliffDDM(ν, α, τ, z, η, sz, st, σ), rt; ϵ) -# end -# return pdf_sv(d, rt; ϵ) -# end - -function _pdf_sv(d::RatcliffDDM{T}, rt::Real; ϵ::Real = 1.0e-12) where {T<:Real} - (ν, α, τ, z, η, sz, st, σ) = params(d) - - if η == 0 - return _pdf(SequentialSamplingModels.DDM(ν, α, τ, z), rt; ϵ) - end - # if isless(ν,0) - # return pdf(SequentialSamplingModels.DDM(ν, α, τ, z), t; ϵ) + ( ( (α*z*η)^2 - 2*ν*α*z - (ν^2)*t ) / (2*(η^2)*t+2) ) - log(sqrt((η^2)*t+1)) + ν*α*z + (ν^2)*t*0.5 - # end - # return pdf(SequentialSamplingModels.DDM(ν, α, τ, z), t; ϵ) + ( ( (α*(1-z)*η)^2 + 2*ν*α*(1-z) - (ν^2)*t ) / (2*(η^2)*t+2) ) - log(sqrt((η^2)*t+1)) - ν*α*(1-z) + (ν^2)*t*0.5 - return _pdf(SequentialSamplingModels.DDM(ν, α, τ, z), rt; ϵ) + ( ( (α*z*η)^2 - 2*ν*α*z - (ν^2)*(rt-τ) ) / (2*(η^2)*(rt-τ)+2) ) - log(sqrt((η^2)*(rt-τ)+1)) + ν*α*z + (ν^2)*(rt-τ)*0.5 -end - -function pdf(d::RatcliffDDM, choice, rt; ϵ::Real = 1.0e-12) - if choice == 1 - (ν, α, τ, z, η, sz, st, σ) = params(d) - return _pdf(RatcliffDDM(-ν, α, τ, 1-z, η, sz, st, σ), rt; ϵ) - end - return _pdf(d, rt; ϵ) -end - -#use numerical integration for variability in non-decision time and bias (Ratcliff and Tuerlinckx, 2002) -function _pdf(d::RatcliffDDM{T}, rt; ϵ::Real = 1.0e-12, n_st::Int=2, n_sz::Int=2) where {T<:Real} - (ν, α, τ, z, η, sz, st, σ) = params(d) - - if τ ≥ rt - return T(NaN) - end - - if st < 1.0e-3 - st = 0 - end - if sz < 1.0e-3 - sz = 0 - end - - if sz==0 - if st==0 #sv=0,sz=0,st=0 - return _pdf_sv(d, rt; ϵ) - else #sv=0,sz=0,st=$ - return _simpson_1D(rt, ν, η, α, z, τ, ϵ, z, z, 0, τ-st/2., τ+st/2., n_st) - end - else #sz=$ - if st==0 #sv=0,sz=$,st=0 - return _simpson_1D(rt, ν, η, α, z, τ, ϵ, z-sz/2., z+sz/2., n_sz, τ, τ , 0) - else #sv=0,sz=$,st=$ - return _simpson_2D(rt, ν, η, α, z, τ, ϵ, z-sz/2., z+sz/2., n_sz, τ-st/2., τ+st/2., n_st) - end - end -end - -##################################################### -# Numerical Integration with Simpson's Method # -# 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) - - # @assert (n_sz & 1) == 0 && (n_st & 1) == 0 # n_st and n_sz have to be even - - n = max(n_st, n_sz) - - if n_st == 0 #integration over z - hz = (ub_z-lb_z)/n - ht = 0 - lb_t = τ - ub_t = τ - else #integration over t - hz = 0 - ht = (ub_t-lb_t)/n - lb_z = z - ub_z = z - end - - S = _pdf_sv(RatcliffDDM(ν, α, lb_t, lb_z, η, 0, 0, 1), x; ϵ) - - y = 0 - z_tag = 0 - t_tag = 0 - - for i in 1:n - z_tag = lb_z + hz * i - t_tag = lb_t + ht * i - - y = _pdf_sv(RatcliffDDM(ν, α, t_tag, z_tag, η, 0, 0, 1), x; ϵ) - - if isodd(i) - S += 4 * y - else - S += 2 * y - end - - end - - 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 - -end - -# 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) - # @assert (n_sz & 1) == 0 && (n_st & 1) == 0 # n_st and n_sz have to be even - # @assert (ub_t-lb_t)*(ub_z-lb_z)>0 && (n_sz*n_st)>0 # 2D-integration only - - ht = (ub_t-lb_t)/n_st - S = _simpson_1D(x, ν, η, α, z, τ, ϵ, lb_z, ub_z, n_sz, 0, 0, 0) - - t_tag = 0 - y = 0 - for i_t in 1:n_st - t_tag = lb_t + ht * i_t - y = _simpson_1D(x, ν, η, α, z, t_tag, ϵ, lb_z, ub_z, n_sz, 0, 0, 0) - - if isodd(i_t) - S += 4 * y - else - S += 2 * y - end - end - - 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 - -end - -logpdf(d::RatcliffDDM, choice, rt; ϵ::Real = 1.0e-12) = log(pdf(d, choice, rt; ϵ)) - -function logpdf(d::RatcliffDDM, data::T) where {T<:NamedTuple} - return sum(logpdf.(d, data...)) -end - -function logpdf(dist::RatcliffDDM, data::Array{<:Tuple,1}) - LL = 0.0 - for d in data - LL += logpdf(dist, d...) - end - return LL -end - -logpdf(d::RatcliffDDM, data::Tuple) = logpdf(d, data...) - -################################################################################ -# Calculate Cumulative Distribution Function # -# # -# Computes Cumulative Distribution Function for the Ratcliff Diffusion model # -# using 6 Gaussian quadrature for numerical integration # -# # -# References # -# 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. # -# # -# Converted from cdfdif.c C script by Joachim Vandekerckhove # -# See also https://ppw.kuleuven.be/okp/software/dmat/ # -################################################################################ - -function cdf(d::RatcliffDDM, choice, rt, p_outlier; w_outlier::Real = 0.1, ϵ::Real = 1e-7) - - (ν, α, τ, z, η, sz, st, σ) = params(d) - - DT = x - τ #make into decision time - - sum = 0.0 - for i = -N:N÷2 - d = 2*i + z - sum += exp(-d*d / (2*DT)) * d - end - return sum / sqrt(2π*DT*DT*DT) -end - -function _g_minus_large_time(x::Real, d::RatcliffDDM, N::Int) - """ - calculate the densities g- for the first exit time for large time values - """ - DT = x - τ #make into decision time - - sum = 0.0 - for i = 1:N - d = i * π - sum += exp(-0.5 * d*d * DT) * sin(d*z) * i - end - return sum * π -end - -function _cdf(d::RatcliffDDM{T}, choice, rt, prob; ϵ::Real = 1e-7) where {T<:Real} - - (ν, α, τ, z, η, sz, st, σ) = params(d) - #Explcit recode of the choice from 2(lower) & 1(upper) to 0(lower) and 1(upper) - #note we need to make sure this is consistent in the all the relative bound models - if choice == 2 #lower - choice = 0 - elseif choice == 1 #upper - choice = 1 - end - - # Initializing variables - a2 = α*α - Z_U = (1-choice)*z+choice*(α-z)+sz/2 - Z_L = (1-choice)*z+choice*(α-z)-sz/2 - lower_t = τ-st/2 - upper_t = 0.0 - Δ = 1e-29 - min_rt=0.001 - v_max = 5000 # maximum number of terms in a partial sum approximating infinite series - - Fnew = 0.0 - sum_z=0.0 - sum_ν=0.0 - p1 = 0.0 - p0 = 0.0 - sum_hist = zeros(3) - denom = 0.0 - sifa = 0.0 - upp = 0.0 - low = 0.0 - fact = 0.0 - exdif = 0.0 - su = 0.0 - sl = 0.0 - zzz = 0.0 - ser = 0.0 - nr_ν = 6 - nr_z = 6 - - # Defining Gauss-Hermite abscissae and weights for numerical integration - gk = [-2.3506049736744922818,-1.3358490740136970132,-.43607741192761650950,.43607741192761650950,1.3358490740136970132,2.3506049736744922818] - w_gh = [.45300099055088421593e-2,.15706732032114842368,.72462959522439207571,.72462959522439207571,.15706732032114842368,.45300099055088421593e-2] - gz = [-.93246951420315193904,-.66120938646626381541,-.23861918608319693247,.23861918608319712676,.66120938646626459256,.93246951420315160597] - w_g = [.17132449237917049545,.36076157304813916138,.46791393457269092604,.46791393457269092604,.36076157304813843973,.17132449237917132812] - - # Adjusting Gauss-Hermite abscissae and weights - for i=1:nr_ν - gk[i] = 1.41421356237309505*gk[i]*η+ν - w_gh[i] = w_gh[i]/1.772453850905515882 - end - for i=1:nr_z - gz[i] = (.5*sz*gz[i])+z - end - - # numerical integration - for i=1:nr_z - sum_ν=0.0 - # numerical integration - for m=1:nr_ν - if abs(gk[m])>ϵ - sum_ν+=(exp(-200*gz[i]*gk[m])-1)/(exp(-200*α*gk[m])-1)*w_gh[m] - else - sum_ν+=gz[i]/α*w_gh[m] - end - end - sum_z+=sum_ν*w_g[i]/2 - end - prob = sum_z - - if (rt-τ+st/2 > min_rt) # is t larger than lower boundary τ distribution? - upper_t = min(rt, τ+st/2) - p1 = prob*(upper_t-lower_t)/st # integrate probability with respect to t - p0 = (1-prob)*(upper_t-lower_t)/st - if rt > τ+st/2 # is t larger than upper boundary Ter distribution? - sum_hist = zeros(3) - for v in 1:v_max # infinite series - sum_hist = circshift(sum_hist, 1) - sum_ν = 0 - sifa = π*v/α - for m in 1:nr_ν # numerical integration with respect to xi - denom = (100*gk[m]*gk[m] + (π*π)*(v*v)/(100*a2)) - upp = exp((2*choice-1)*Z_U*gk[m]*100 - 3*log(denom) + log(w_gh[m]) - 2*log(100)) - low = exp((2*choice-1)*Z_L*gk[m]*100 - 3*log(denom) + log(w_gh[m]) - 2*log(100)) - fact = upp*((2*choice-1)*gk[m]*sin(sifa*Z_U)*100 - sifa*cos(sifa*Z_U)) - - low*((2*choice-1)*gk[m]*sin(sifa*Z_L)*100 - sifa*cos(sifa*Z_L)) - exdif = exp((-.5*denom*(rt-upper_t)) + log(1-exp(-.5*denom*(upper_t-lower_t)))) - sum_ν += fact*exdif - end - sum_hist[3] = sum_hist[2] + v*sum_ν - if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 - break - end - end - - Fnew = (p0*(1-choice) + p1*choice) - sum_hist[3]*4*π/(a2*sz*st) - # cumulative distribution function for t and x - elseif t <= τ+st/2 # is t lower than upper boundary Ter distribution? - sum_ν = 0 - for m in 1:nr_ν - if abs(gk[m]) > ϵ - sum_z = 0 - for i in 1:nr_z - zzz = (α - gz[i])*choice + gz[i]*(1 - choice) - ser = -((α*a2)/((1 - 2*choice)*gk[m]*π*.01))*sinh(zzz*(1 - 2*x)*gk[m]/.01)/ - (sinh((1 - 2*choice)*gk[m]*α/.01)^2) + - (zzz*a2)/((1 - 2*choice)*gk[m]*π*.01)*cosh((α - zzz)*(1 - 2*choice)*gk[m]/.01)/ - sinh((1 - 2*choice)*gk[m]*α/.01) - sum_hist = zeros(3) - for v in 1:v_max - sum_hist = circshift(sum_hist, 1) - sifa = π*v/α - denom = (gk[m]*gk[m]*100 + (π*v)*(π*v)/(a2*100)) - sum_hist[3] = sum_hist[2] + v*sin(sifa*zzz)*exp(-.5*denom*(rt - lower_t) - 2*log(denom)) - if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 - break - end - end - sum_z += .5*w_g[i]*(ser - 4*sum_hist[3])*(π/100)/(a2*st)*exp((2*choice - 1)*zzz*gk[m]*100) - end - else - sum_hist = zeros(3) - su = -(Z_U*Z_U)/(12*a2) + (Z_U*Z_U*Z_U)/(12*α*a2) - (Z_U*Z_U*Z_U*Z_U)/(48*a2*a2) - sl = -(Z_L*Z_L)/(12*a2) + (Z_L*Z_L*Z_L)/(12*α*a2) - (Z_L*Z_L*Z_L*Z_L)/(48*a2*a2) - for v in 1:v_max - sum_hist = circshift(sum_hist, 1) - sifa = π*v/α - denom = (π*v)*(π*v)/(a2*100) - sum_hist[3] = sum_hist[2] + 1/(π*π*π*π*v*v*v*v)*(cos(sifa*Z_L) - cos(sifa*Z_U))* - exp(-.5*denom*(rt - lower_t)) - if abs(sum_hist[1] - sum_hist[2]) < Δ && abs(sum_hist[2] - sum_hist[3]) < Δ && sum_hist[3] > 0 - break - end - end - sum_z = 400*a2*α*(sl - su - sum_hist[3])/(st*sz) - end - sum_ν += sum_z*w_gh[m] - end - Fnew = (p0*(1 - choice) + p1*choice) - sum_ν - end - elseif rt - τ + st/2 <= min_rt # is t lower than lower boundary Ter distr? - Fnew = 0 - end - - Fnew = Fnew > Δ ? Fnew : 0 - - return Fnew - -end - -function _add_outlier_cdf(y::Real, x::Real, p_outlier::Real; w_outlier::Real = 0.1) - #Ratcliff and Tuerlinckx, 2002 containment process - return y * (1 - p_outlier) + (x + (1. / (2 * w_outlier))) * w_outlier * p_outlier -end - -function _p_outlier_in_range(p_outlier) - return (p_outlier >= 0) & (p_outlier <= 1) -end - -""" - rand(dist::RatcliffDDM) - -Generate a random choice and rt for the Ratcliff Diffusion Model - -# Arguments -- `dist`: model object for Ratcliff Diffusion Model. -- `method`: method simulating the diffusion process. - "rejection" uses Tuerlinckx et al., 2001 rejection-based method for the general wiener process - "stochastic" uses the stochastic Euler method to directly simulate the stochastic differential equation - -# References - - Tuerlinckx, F., Maris, E., Ratcliff, R., & De Boeck, P. (2001). - A comparison of four methods for simulating the diffusion process. - Behavior Research Methods, Instruments, & Computers, 33, 443-456. - - Converted from Rhddmjagsutils.R R script by Kianté Fernandez - - See also https://github.com/kiante-fernandez/Rhddmjags. -""" -function rand(rng::AbstractRNG, d::RatcliffDDM) - # method::Char = "rejection" - return _rand_rejection(rng, d) - # method::Char = "stochastic" -# return _rand_stochastic(rng, d) -end - -function _rand_rejection(rng::AbstractRNG, d::RatcliffDDM; N::Int = 1) - (ν, α, τ, z, η, sz, st, σ) = params(d) - - if η == 0 - η = 1e-16 - end - - # Initialize output vectors - result = zeros(N) - T = zeros(N) - XX = zeros(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() - μ = ν + r1 * η - bb = z - sz / 2 + sz * rand() - zz = bb * α - finish = 0 - totaltime = 0 - startpos = 0 - Aupper = α - zz - Alower = -zz - radius = min(abs(Aupper), abs(Alower)) - - while finish == 0 - λ = 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 * μ) - F = F^2 / (1 + F^2) - # formula p447 in 2001 paper - prob = exp(radius * μ / D) - prob = prob / (1 + prob) - dir_ = 2 * (rand() < prob) - 1 - l = -1 - s2 = 0 - s1 = 0 - while s2 > l - s2 = rand() - s1 = rand() - tnew = 0 - told = 0 - uu = 0 - while abs(tnew - told) > ϵ || uu == 0 - told = tnew - uu += 1 - tnew = told + (2 * uu + 1) * (-1)^uu * s1^(F * (2 * uu + 1)^2) - # infinite sum in formula (16) in BRMIC,2001 - end - l = 1 + s1^(-F) * tnew - end - # rest of formula (16) - t = abs(log(s1)) / λ - # is the negative of t* in (14) in BRMIC,2001 - totaltime += t - dir_ = startpos + dir_ * radius - ndt = τ - st / 2 + st * rand() - if (dir_ + Δ) > Aupper - T[n] = ndt + totaltime - XX[n] = 1 - finish = 1 - elseif (dir_ - Δ) < Alower - T[n] = ndt + totaltime - XX[n] = 2 - finish = 1 - else - startpos = dir_ - radius = minimum(abs.([Aupper, Alower] .- startpos)) - end - end - end - return (choice=XX,rt=T) -end - -function _rand_stochastic(rng::AbstractRNG, d::RatcliffDDM; N::Int = 1, nsteps::Int=300, step_length::Int=0.01) - (ν, α, τ, z, η, sz, st, σ) = params(d) - - if η == 0 - η = 1e-16 - end - - # Initialize output vectors - choice = fill(0, N) - rt = fill(0.0, N) - - for n in 1:N - random_walk = Array{Float64}(undef, nsteps) - start_point = (z - sz/2) + ((z + sz/2) - (z - sz/2)) * rand() - ndt = (τ - st/2) + ((τ + st/2) - (τ - st/2)) * rand() - drift = rand(Distributions.Normal(ν, η)) - random_walk[1] = start_point * α - for s in 2:nsteps - random_walk[s] = random_walk[s-1] + rand(Distributions.Normal(drift * step_length, σ * sqrt(step_length))) - if random_walk[s] >= α - random_walk[s:end] .= α - rts[n] = s * step_length + ndt - choice[n] = 1 - break - elseif random_walk[s] <= 0 - random_walk[s:end] .= 0 - rts[n] = s * step_length + ndt - choice[n] = 2 - break - elseif s == nsteps - rts[n] = NaN - choice[n] = NaN - break - end - end - end - return (choice=choice,rt=rts) -end - -""" - rand(dist::RatcliffDDM, n_sim::Int) - -Generate `n_sim` random choice-rt pairs for the Ratcliff Diffusion Decision Model. - -# Arguments -- `dist`: model object for the Ratcliff DDM. -- `n_sim::Int`: the number of simulated rts -""" - -function rand(rng::AbstractRNG, d::RatcliffDDM, n_sim::Int) - return _rand_rejection(rng, d, N = n_sim) -end - -sampler(rng::AbstractRNG, d::RatcliffDDM) = rand(rng::AbstractRNG, d::RatcliffDDM) \ No newline at end of file diff --git a/src/SequentialSamplingModels.jl b/src/SequentialSamplingModels.jl index 050ec352..357f941c 100644 --- a/src/SequentialSamplingModels.jl +++ b/src/SequentialSamplingModels.jl @@ -36,6 +36,7 @@ module SequentialSamplingModels export AbstractLNR export AbstractPoissonRace export AbstractRDM + export AbstractDDM export AbstractWald export aDDM export CDDM @@ -50,6 +51,7 @@ module SequentialSamplingModels export SSM1D export SSM2D export ContinuousMultivariateSSM + # export RatcliffDDM export Wald export WaldMixture @@ -91,5 +93,5 @@ module SequentialSamplingModels include("ext_functions.jl") include("ex_gaussian.jl") include("poisson_race.jl") - include("RatcliffDDM.jl") + # include("RatcliffDDM.jl") end \ No newline at end of file diff --git a/src/type_system.jl b/src/type_system.jl index f2fa1dc1..06f2a694 100644 --- a/src/type_system.jl +++ b/src/type_system.jl @@ -72,6 +72,13 @@ An abstract type for the racing diffusion model. """ abstract type AbstractRDM <: SSM2D end +""" + AbstractDDM <: SSM2D + +An abstract type for the diffusion decision model. +""" +abstract type AbstractDDM <: SSM2D end + abstract type PDFType end """ diff --git a/src/utilities.jl b/src/utilities.jl index 38ff18a7..59ef1cb4 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,3 +1,281 @@ +abstract type Mixed <: ValueSupport end + +""" + SSM2D = Distribution{Multivariate, Mixed} + +An abstract type for sequential sampling models characterized by a multivariate choice-reaction time distribution. +Sub-types of `SSM2D` output a `NamedTuple` consisting of a vector of choices and reaction times. +""" +const SSM2D = Distribution{Multivariate, Mixed} + +""" + SSM1D <: ContinuousUnivariateDistribution + +An abstract type for sequential sampling models characterized by a single choice reaction time distribution. +Sub-types of `SSM1D` output a vector of reaction times. +""" +abstract type SSM1D <: ContinuousUnivariateDistribution end + +abstract type AbstractWald <: SSM1D end + +abstract type AbstractaDDM <: SSM2D end + +abstract type PDFType end + +""" + Exact <: PDFType + +Has closed-form PDF. +""" +struct Exact <: PDFType end + +""" + Approximate <: PDFType + +Has approximate PDF based on kernel density estimator. +""" +struct Approximate <: PDFType end + +get_pdf_type(d::SSM1D) = Exact +get_pdf_type(d::SSM2D) = Exact + +minimum(d::SSM1D) = 0.0 +maximum(d::SSM1D) = Inf + +minimum(d::SSM2D) = 0.0 +maximum(d::SSM2D) = Inf + +insupport(d::SSM2D, data) = data.rt ≥ minimum(d) && data.rt ≤ maximum(d) + +Base.broadcastable(x::SSM1D) = Ref(x) +Base.broadcastable(x::SSM2D) = Ref(x) + +vectorize(d::SSM2D, r::NamedTuple) = [r...] + +Base.length(d::SSM2D) = 2 + +rand(d::SSM2D) = rand(Random.default_rng(), d) +rand(d::SSM2D, n::Int) = rand(Random.default_rng(), d, n) + +""" + logpdf(d::SSM2D, data::NamedTuple) + +Computes the likelihood for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +logpdf(d::SSM2D, data::NamedTuple; kwargs...) = logpdf.(d, data.choice, data.rt; kwargs...) +logpdf(d::SSM2D, data::Vector{Real}; kwargs...) = logpdf(d, Int(data[1]), data[2]; kwargs...) + +""" + loglikelihood(d::SSM2D, data::NamedTuple) + +Computes the summed log likelihood for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +loglikelihood(d::SSM2D, data::NamedTuple) = sum(logpdf.(d, data...)) + +""" + pdf(d::SSM2D, data::NamedTuple) + +Computes the probability density for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +pdf(d::SSM2D, data::NamedTuple) = pdf.(d, data.choice, data.rt) + + +""" + 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_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) +end + +""" + n_options(dist::SSM2D) + +Returns the number of choice options based on the length of the drift rate vector `ν`. + +# Arguments + +- `d::SSM2D`: a sub-type of `SSM2D` +""" +n_options(d::SSM2D) = length(d.ν) + +""" + n_options(dist::SSM1D) + +Returns 1 for the number of choice options + +# Arguments + +- `d::SSM1D`: a sub-type of `SSM1D` +""" +n_options(d::SSM1D) = 1 + +abstract type Mixed <: ValueSupport end + +""" + SSM2D = Distribution{Multivariate, Mixed} + +An abstract type for sequential sampling models characterized by a multivariate choice-reaction time distribution. +Sub-types of `SSM2D` output a `NamedTuple` consisting of a vector of choices and reaction times. +""" +const SSM2D = Distribution{Multivariate, Mixed} + +""" + SSM1D <: ContinuousUnivariateDistribution + +An abstract type for sequential sampling models characterized by a single choice reaction time distribution. +Sub-types of `SSM1D` output a vector of reaction times. +""" +abstract type SSM1D <: ContinuousUnivariateDistribution end + +abstract type AbstractaDDM <: SSM2D end + +abstract type AbstractLBA <: SSM2D end + +abstract type AbstractWald <: SSM1D end + +abstract type PDFType end + +""" + Exact <: PDFType + +Has closed-form PDF. +""" +struct Exact <: PDFType end + +""" + Approximate <: PDFType + +Has approximate PDF based on kernel density estimator. +""" +struct Approximate <: PDFType end + +get_pdf_type(d::SSM1D) = Exact +get_pdf_type(d::SSM2D) = Exact + +minimum(d::SSM1D) = 0.0 +maximum(d::SSM1D) = Inf + +minimum(d::SSM2D) = 0.0 +maximum(d::SSM2D) = Inf + +insupport(d::SSM2D, data) = data.rt ≥ minimum(d) && data.rt ≤ maximum(d) + +Base.broadcastable(x::SSM1D) = Ref(x) +Base.broadcastable(x::SSM2D) = Ref(x) + +vectorize(d::SSM2D, r::NamedTuple) = [r...] + +Base.length(d::SSM2D) = 2 + +rand(d::SSM2D) = rand(Random.default_rng(), d) +rand(d::SSM2D, n::Int) = rand(Random.default_rng(), d, n) + +""" + logpdf(d::SSM2D, data::NamedTuple) + +Computes the likelihood for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +logpdf(d::SSM2D, data::NamedTuple; kwargs...) = logpdf.(d, data.choice, data.rt; kwargs...) +logpdf(d::SSM2D, data::Vector{Real}; kwargs...) = logpdf(d, Int(data[1]), data[2]; kwargs...) + +""" + loglikelihood(d::SSM2D, data::NamedTuple) + +Computes the summed log likelihood for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +loglikelihood(d::SSM2D, data::NamedTuple) = sum(logpdf.(d, data...)) + +""" + pdf(d::SSM2D, data::NamedTuple) + +Computes the probability density for a 2D sequential sampling model. + +# Arguments + +- `d::SSM2D`: an object for a 2D sequential sampling model +- `data::NamedTuple`: a NamedTuple of data containing choice and reaction time +""" +pdf(d::SSM2D, data::NamedTuple) = pdf.(d, data.choice, data.rt) + + +""" + 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_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) +end + +""" + n_options(dist::SSM2D) + +Returns the number of choice options based on the length of the drift rate vector `ν`. + +# Arguments + +- `d::SSM2D`: a sub-type of `SSM2D` +""" +n_options(d::SSM2D) = length(d.ν) + +""" + n_options(dist::SSM1D) + +Returns 1 for the number of choice options + +# Arguments + +- `d::SSM1D`: a sub-type of `SSM1D` +""" +n_options(d::SSM1D) = 1 + function Base.show(io::IO, ::MIME"text/plain", model::SSM1D) return _show(io::IO, model) end diff --git a/test/ddm_tests.jl b/test/ddm_tests.jl index aae27db0..3d26a152 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, α = .8, z = .5, τ = .3) + dist = DDM(ν=1.0, α = .8, z = .5, τ = .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 = .5, τ = .30) + dist = DDM(ν=2.0, α = 1.5, z = .5, τ = .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) @@ -49,6 +49,106 @@ @test y′ ≈ y rtol = .05 end + @safetestset "Full DDM pdf 1 (η)" begin + using SequentialSamplingModels + using Test + using Random + Random.seed!(654) + include("KDE.jl") + + dist = DDM(ν=1.0, α = .8, z = .5, τ = .3, η = 0.08, sz = 0.00, st = 0.00, σ = 1.0) + choice,rt = rand(dist, 10^6) + rt1 = rt[choice .== 1] + p1 = mean(choice .== 1) + p2 = 1 - p1 + approx_pdf = kernel(rt1) + x = range(.301, 1.5, length=100) + y′ = pdf(approx_pdf, x) * p1 + y = pdf.(dist, (1,), x) + @test y′ ≈ y rtol = .05 + + rt2 = rt[choice .== 2] + approx_pdf = kde(rt2) + y′ = pdf(approx_pdf, x) * p2 + y = pdf.(dist, (2,), x) + @test y′ ≈ y rtol = .05 + end + + @safetestset "Full DDM pdf 2 (sz)" begin + using SequentialSamplingModels + using Test + using Random + Random.seed!(654) + include("KDE.jl") + + dist = DDM(ν=1.0, α = .8, z = .5, τ = .3, η = 0.00, sz = 0.10, st = 0.00, σ = 1.0) + choice,rt = rand(dist, 10^6) + rt1 = rt[choice .== 1] + p1 = mean(choice .== 1) + p2 = 1 - p1 + approx_pdf = kernel(rt1) + x = range(.301, 1.5, length=100) + y′ = pdf(approx_pdf, x) * p1 + y = pdf.(dist, (1,), x) + @test y′ ≈ y rtol = .05 + + rt2 = rt[choice .== 2] + approx_pdf = kde(rt2) + y′ = pdf(approx_pdf, x) * p2 + y = pdf.(dist, (2,), x) + @test y′ ≈ y rtol = .05 + end + + @safetestset "Full DDM pdf 3 (st)" begin + using SequentialSamplingModels + using Test + using Random + Random.seed!(654) + include("KDE.jl") + + dist = DDM(ν=1.0, α = .8, z = .5, τ = .3, η = 0.00, sz = 0.00, st = 0.10, σ = 1.0) + choice,rt = rand(dist, 10^6) + rt1 = rt[choice .== 1] + p1 = mean(choice .== 1) + p2 = 1 - p1 + approx_pdf = kernel(rt1) + x = range(.301 + .10, 1.5, length=100) # τ and st start the lower bound + y′ = pdf(approx_pdf, x) * p1 + y = pdf.(dist, (1,), x) + @test y′ ≈ y rtol = .05 + + rt2 = rt[choice .== 2] + approx_pdf = kde(rt2) + y′ = pdf(approx_pdf, x) * p2 + y = pdf.(dist, (2,), x) + @test y′ ≈ y rtol = .05 + end + + @safetestset "Full DDM pdf 4 (η,sz,st)" begin + using SequentialSamplingModels + using Test + using Random + Random.seed!(654) + include("KDE.jl") + + dist = DDM(ν=1.0, α = .8, z = .5, τ = .3, η = 0.08, sz = 0.10, st = 0.10, σ = 1.0) + choice,rt = rand(dist, 10^6) + rt1 = rt[choice .== 1] + p1 = mean(choice .== 1) + p2 = 1 - p1 + approx_pdf = kernel(rt1) + x = range(.301 + .10, 1.5, length=100) # τ and st start the lower bound + y′ = pdf(approx_pdf, x) * p1 + y = pdf.(dist, (1,), x) + @test y′ ≈ y rtol = .05 + + rt2 = rt[choice .== 2] + approx_pdf = kde(rt2) + y′ = pdf(approx_pdf, x) * p2 + y = pdf.(dist, (2,), x) + @test y′ ≈ y rtol = .05 + end + @safetestset "DDM cdf 1" begin using SequentialSamplingModels using Test @@ -56,7 +156,7 @@ using Random Random.seed!(7540) - dist = DDM(ν=1.0, α = .8, z = .5, τ = .3) + dist = DDM(ν=1.0, α = .8, z = .5, τ = .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 +181,7 @@ using Random Random.seed!(2200) - dist = DDM(ν=2.0, α = 1.5, z = .5, τ = .30) + dist = DDM(ν=2.0, α = 1.5, z = .5, τ = .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) @@ -96,21 +196,21 @@ ecdf2 = ecdf(rt2) y′ = ecdf1.(x) * p2 y = cdf.(dist, (2,), x) - @test y′ ≈ y rtol = .01 + @test y′ ≈ y rtol = .02 #note the change in tol 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) @@ -127,17 +227,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, α=.50, z=.25, τ=.50), .75; ϵ = 1e-10) - @test test_val1 ≈ 3.0 + # test_val1 = _K_large(DDM(;ν=1.5, α=.50, z=.25, τ=.50), .75; ϵ = 1e-10) + # @test test_val1 ≈ 3.0 - test_val2 = _K_large(DDM(;ν=1.5, α=.50, z=.25, τ=.50), .501; ϵ = 1e-10) - @test test_val2 ≈ 36 - end + # test_val2 = _K_large(DDM(;ν=1.5, α=.50, z=.25, τ=.50), .501; ϵ = 1e-10) + # @test test_val2 ≈ 36 + # end # K_large = function(t, v, a, w, epsilon) # { # sqrtL1 = sqrt(1/t) * a/pi @@ -147,21 +247,21 @@ # 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, α=.50, z=.25, τ=.50), .75; ϵ = 1e-10) - @test test_val1 ≈ 16 + # # function(t, v, a, w, epsilon) + # test_val1 = _K_small(DDM(;ν=1.5, α=.50, z=.25, τ=.50, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), .75; ϵ = 1e-10) + # @test test_val1 ≈ 16 - test_val2 = _K_small(DDM(;ν=1.5, α=.50, z=.25, τ=.50), .501; ϵ = 1e-10) - @test test_val2 ≈ 16 + # test_val2 = _K_small(DDM(;ν=1.5, α=.50, z=.25, τ=.50, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), .501; ϵ = 1e-10) + # @test test_val2 ≈ 16 - test_val3 = _K_small(DDM(;ν=0.50, α=2.50, z=.50, τ=.30), .400; ϵ = 1e-10) - @test test_val3 ≈ 10 - end + # test_val3 = _K_small(DDM(;ν=0.50, α=2.50, z=.50, τ=.30, η = 0.00, sz = 0.0, st = 0.0, σ = 1.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 @@ -178,35 +278,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., .5, .5) - @test test_val1 ≈ 0.3775407 atol = 1e-5 + # test_val1 = _P_upper(1., .5, .5) + # @test test_val1 ≈ 0.3775407 atol = 1e-5 - test_val2 = _P_upper(-1., .75, .25) - @test test_val2 ≈ .8693188 atol = 1e-5 + # test_val2 = _P_upper(-1., .75, .25) + # @test test_val2 ≈ .8693188 atol = 1e-5 - test_val3 = _P_upper(-1e10, .75, .25) - @test test_val3 ≈ 1 + # test_val3 = _P_upper(-1e10, .75, .25) + # @test test_val3 ≈ 1 - test_val4 = _P_upper(eps(), .75, .25) - @test test_val4 ≈ .75 - end + # test_val4 = _P_upper(eps(), .75, .25) + # @test test_val4 ≈ .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, α=.50, z=.25, τ=.50), 10, .75) - @test test_val1 ≈ 0.5955567 atol = 1e-5 + # test_val1 = _Fs_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .75) + # @test test_val1 ≈ 0.5955567 atol = 1e-5 - test_val2 = _Fs_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .501) - @test test_val2 ≈ 6.393096e-05 atol = 1e-8 - end + # test_val2 = _Fs_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .501) + # @test test_val2 ≈ 6.393096e-05 atol = 1e-8 + # end # Fs_lower = function(t, v, a, w, K) # { @@ -228,17 +328,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, α=.50, z=.25, τ=.50), 10, .75) - @test test_val1 ≈ 0.5955567 atol = 1e-5 + # test_val1 = _Fl_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .75) + # @test test_val1 ≈ 0.5955567 atol = 1e-5 - test_val2 = _Fl_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .501) - @test test_val2 ≈ 0.001206451 atol = 1e-8 - end + # test_val2 = _Fl_lower(DDM(;ν=1.5, α=.50, z=.25, τ=.50), 10, .501) + # @test test_val2 ≈ 0.001206451 atol = 1e-8 + # end # Fl_lower = function(t, v, a, w, K) # { # F = numeric(length(t)) @@ -254,19 +354,103 @@ using SequentialSamplingModels using Test # tested against rtdists - test_val1 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3), 1, .5) + test_val1 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), 1, .5) @test test_val1 ≈ 2.131129 atol = 1e-5 - test_val2 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3), 2, .5) + test_val2 = pdf(DDM(;ν=2.0, α=1.0, z=.5, τ=.3, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), 2, .5) @test test_val2 ≈ 0.2884169 atol = 1e-5 - test_val3 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2), 1, .35) + test_val3 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), 1, .35) @test test_val3 ≈ 0.6635714 atol = 1e-5 - test_val4 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2), 2, .35) + test_val4 = pdf(DDM(;ν=.8, α=.5, z=.3, τ=.2, η = 0.00, sz = 0.0, st = 0.0, σ = 1.0), 2, .35) @test test_val4 ≈ 0.4450956 atol = 1e-5 end + @safetestset "full pdf" begin + using SequentialSamplingModels + 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_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_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.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] + + # 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) + # 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) + # @test 2.537052 ≈ pdf_value atol = 1e-1 + + end + @safetestset "simulate" begin using SequentialSamplingModels using Test