## DiffModels.jl

Diffusion Model simulation and first-passage time densities in Julia
Author DrugowitschLab
Popularity
6 Stars
Updated Last
4 Months Ago
Started In
July 2014

# DiffModels.jl   A Julia package for simulating diffusion models and compute their first passage time densities.

No guarantee is provided for the correctness of the implementation.

## Content

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).

## Installation

### Julia 0.7.x and above

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
``````

### Julia 0.6.x

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.

## Usage

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.

### Drift

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`.

### Boundaries

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.

### First-passage time densities

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.

### Drawing first-passage time and boundary samples

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.

## References

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.