A Julia package for simulating diffusion models and compute their first passage time densities.
No guarantee is provided for the correctness of the implementation.
The code is licensed under the MIT License.
The library provides Julia classes and functions to define diffusion models, sample their first-passage times and boundaries, and compute their first-passage time density functions. For now, only diffusion models with two absorbing boundaries are supported. The library supports time-changing boundaries and drifts.
The diffusion models assume a drifting and diffusing particle x(t) that starts at x(0) = 0and whose time-course follows
dx = mu(t) dt + sig(t) dW ,
where mu(t) is the current drift, sig(t) is the current diffusion standard deviation (for now assumed to be sig(t)=1, always), and dW is a Wiener process. Diffusion is terminated as soon as the particle reaches either the upper boundary theta_u(t) or lower boundary theta_t(t). The library requires theta_u(0) > 0 and theta_l(0) < 0. The time at which either boundary is reached is the first-passage time. The associated densities, g_u(t) and g_l(t), are the joint densities over bounds and first-passage times, such that
integral_0^infinity (g_u(t) + g_l(t)) dt = 1 .
The library provides specialised classes for time-invariant drifts, mu(t) = mu_0 for all t, time-invariant bounds, theta(t) = theta_0 for all t, symmetric bounds, theta_l(t) = - theta_u(t) (see below).
The easiest way to install DiffModels.jl is by using the Julia Package Manager at the Julia Pkg prompt:
pkg> add https://github.com/DrugowitschLab/DiffModels.jl
For installation with Julia 0.6.x, use the following command at the Julia prompt:
julia> Pkg.clone("git://github.com/DrugowitschLab/DiffModels.jl.git")
Then call Pkg.dir() to find the package installation directly, open a terminal and navigate to this directory. For there, use
# cd DiffModels
# git checkout 1cd0b4a418a9d78f6f19442526cb9de91507bbbe
in the terminal to roll back to the last version compatible with Julia 0.6.x.
The library is based on assembling diffusion models from drift and boundary specifications. If not constant, all drifts/boundaries are specified by vectors in time steps of dt, with the first vector element corresponding to t=0. In most cases, these vectors need to be sufficiently long to cover the whole time of relevance. The first-passage times cannot computed beyond this time, and samples cannot be drawn after it. This restriction does not apply to time-invariant drifts/bounds, which do not feature this limitation.
The available drifts are defined in src/drift.jl. All drifts are based on the abstract base class AbstractDrift. The following constructors are available:
ConstDrift(mu::Real, dt::Real)
VarDrift(mu::Vector{Float64}, dt::Real)
VarDrift(mu::Vector{Float64}, dt::Real, maxt::Real)ConstDrift is a constant drift of size mu. dt needs to be nonetheless specified, to evaluate the step size for diffusion model sampling. VarDrift is a drift that changes over time according to the mu vector that specifies this drift in steps of dt. If maxt is given, the mu vector is extended by repeating its last element until time maxt.
The available boundaries are defined in src/bound.jl. Single boundaries are based on the abstract base class 'AbstractBound'. For such boundaries, the following constructors are available:
ConstBound(b::Real, dt::Real)
LinearBound(b0::Real, bslope::Real, dt::Real)
VarBound(b::Vector{Float64}, bg::Vector{Float64}, dt::Real)
VarBound(b::Vector{Float64}, dt::Real)ConstBound is a constant boundary at b. LinearBound is a linearly changing boundary that, at time t is located at b0 + t * bslope. VarBound is a time-varying boundary that changes over time according to the vector b in steps of dt. bg is its time derivative, and needs to contain the same number of elements as b. If not specified (last constructor), it is estimated from b by finite differences.
Boundary pairs are based on the abstract base class AbstractBounds. For such pairs, the following constructors are available:
SymBounds(b::T) where T <: AbstractBound
typealias VarSymBounds SymBounds{VarBound}
typealias LinearSymBounds SymBounds{LinearBound}
typealias ConstSymBounds SymBounds{ConstBound}
AsymBounds(upper::T1, lower::T2) where {T1 <: AbstractBound, T2 <: AbstractBound}
VarAsymBounds AsymBounds{VarBound, VarBound}
ConstAsymBounds AsymBounds{ConstBound, ConstBound}
typealias ConstBounds Union(ConstSymBounds, ConstAsymBounds)SymBounds and AsymBounds specify symmetric boundaries (around zero) and asymmetric boundaries, respectively. SymBounds needs to be constructed with the upper boundary, and the lower boundaries is mirrored around zero. AsymBounds is constructed with two boundaries, where upper is the upper boundary, and lower is the negative lower boundary.
In the most general case, the first-passage time densities are computed by
pdf(d::AbstractDrift, b::AbstractBounds, tmax::Real)This function returns the vectors g1 and g2 that contain the densities for the upper and the lower boundary, respecively, in steps of dt (as specified by drift and bounds) up to time tmax.
For some combination of constant drifts and constant boundaries, significantly more efficient functions are available:
pdf(d::ConstDrift, b::ConstBounds, tmax::Real)
pdfu(d::ConstDrift, b::ConstBounds, t::Float64)
pdfl(d::ConstDrift, b::ConstBounds, t::Float64)
pdful(d::ConstDrift, b::ConstBounds, t::Float64)pdf returns two vectors, g1 and g2, as before. pdfu and pdfl compute the first-passage time density at the upper and lower boundary, respecively, only at time t. pdful returns both densities for this time.
Diffusion model sampling is based on the Sampler framework from Distributions.jl. The basic idea to is create a sampler s from a drift/boundary specification, on which rand(s) is called. All calls to rand(s) return a t, bound pair, where t is the first-passage time, and bound is true if the upper boundary was reached, and false otherwise. Currently, the following samplers exist:
sampler(d::AbstractDrift, b::AbstractBounds)
sampler(d::ConstDrift, b::ConstSymBounds)
sampler(d::ConstDrift, b::ConstAsymBounds)The first sampler is a generic diffusion model sampler that draws samples by simulating full trajectories in steps of dt. All samples beyond the time that the drift/bound are specified return t = Inf and a random bound. The other samplers use a specialised and significantly faster method, based on rejection sampling.
In general, the library computes the first-passage time densities by finding the solution to an integral equation, as described in
Smith PL (2000). Stochastic Dynamic Models of Response Time and Accuracy: A Foundational Primer. Journal of Mathematical Psychology, 44 (3). 408-463.
For constant drift and bounds, it instead uses a much faster method, based on an infinite series expansion of these densities, as described in
Cox DR and Miller HD (1965). The Theory of Stochastic Processes. John Wiley & Sons, Inc.
and
Navarro DJ and Fuss IG (2009). Fast and accurate calculations for first-passage times in Wiener diffusion models. Journal of Mathematical Psychology, 53, 222-230.
Samples are in the most general case drawn by simulating trajectories by the Euler–Maruyama method. For diffusion models with constant drift and (symmetric or asymmetric) boundaries, the following significantly faster method based on rejection sampling is used:
Drugowitsch J (2016). Fast and accurate Monte Carlo sampling of first-passage times from Wiener diffusion models. Scientific Reports 6, 20490; doi: 10.1038/srep20490.