GitHub

ForwardBackward.jl

ForwardBackward.jl is a Julia package for evolving discrete and continuous states under a variety of processes.

Overview

The package implements forward and backward methods for a number of useful processes. For times $s < t < u$:

  • Xt = forward(Xs, P, t-s) computes the distribution $\propto P(X_t | X_s)$
  • Xt = backward(Xu, P, u-t) computes the likelihood $\propto P(X_u | X_t)$ (ie. considered as a function of $X_t$)

Where P is a Process, and each of $X_s$, $X_t$, $X_u$ can be DiscreteState or ContinuousState, or scaled distributions (CategoricalLikelihood or GaussianLikelihood) over states, where the uncertainty (and normalizing constants) propogate. States and Likelihoods hold arrays of points/distributions, which are all acted upon by the process.

Since CategoricalLikelihood and GaussianLikelihood are closed under elementwise/pointwise products, to compute the (scaled) distribution at $t$, which is $P(X_t | X_s, X_u) ∝ P(X_t | X_s)P(X_u | X_t)$, we also provide ⊙:

Xt = forward(Xs, P, t-s) ⊙ backward(Xu, P, u-t)

One use-cases is drawing endpoint conditioned samples:

rand(forward(Xs, P, t-s) ⊙ backward(Xu, P, u-t))

or

endpoint_conditioned_sample(Xs, Xu, P, t-s, u-t)

For some processes where we don't support propagation of uncertainty, (eg. the ManifoldProcess), endpoint_conditioned_sample is implemented directly via approximate simulation.

Processes

  • Continuous State
    • BrownianMotion
    • OrnsteinUhlenbeck
    • Deterministic (where endpoint_conditioned_sample interpolates)
  • Discrete State:
    • GeneralDiscrete, with any Q matrix, where propogation is via matrix exponentials
    • UniformDiscrete, with all rates equal
    • PiQ, where any event is a switch to a draw from the stationary distribution
    • UniformUnmasking, where switches occur from a masked state to any other states

Installation

using Pkg
Pkg.add("ForwardBackward")

Quick Start

using ForwardBackward

# Create a Brownian motion process
process = BrownianMotion(0.0, 1.0)  # drift = 0.0, variance = 1.0

# Define start and end states
X0 = ContinuousState(zeros(10))     # start at origin
X1 = ContinuousState(ones(10))      # end at ones

# Sample a path at t = 0.3
sample = endpoint_conditioned_sample(X0, X1, process, 0.3)

Examples

Discrete State Process

# Create a process with uniform transition rates
process = UniformDiscrete()
X0 = DiscreteState(4, [1])    # 4 possible states, starting in state 1
X1 = DiscreteState(4, [4])    # ending in state 4

# Sample intermediate state
sample = endpoint_conditioned_sample(X0, X1, process, 0.5)

Manifold-Valued Process

using Manifolds

# Create a process on a sphere
M = Sphere(2)                  # 2-sphere
process = ManifoldProcess(0.1) # with some noise

# Define start and end points
p0 = [1.0, 0.0, 0.0]
p1 = [0.0, 0.0, 1.0]
X0 = ManifoldState(M, [p0])
X1 = ManifoldState(M, [p1])

# Sample a path
sample = endpoint_conditioned_sample(X0, X1, process, 0.5)

Endpoint-conditioned samples on a torus

using ForwardBackward, Manifolds, Plots

#Project Torus(2) into 3D (just for plotting)
function tor(p; R::Real=2, r::Real=0.5)
    u,v = p[1], p[2]
    x = (R + r*cos(u)) * cos(v)
    y = (R + r*cos(u)) * sin(v)
    z = r * sin(u)
    return [x, y, z]
end

#Define the manifold, and two endpoints, which are on opposite sides (in both dims) of the torus:
M = Torus(2)
p1 = [-pi, 0.0]
p0 = [0.0, -pi]

#When non-zero, the process will diffuse. When 0, the process is deterministic:
for P in [ManifoldProcess(0), ManifoldProcess(0.05)]
    #When non-zero, the endpoints will be slightly noised:
    for perturb_var in [0.0, 0.0001] 
     
        #We'll generate endpoint-conditioned samples evenly spaced over time:
        t_vec = 0:0.001:1

        #Set up the X0 and X1 states, just repeating the endpoints over and over:
        X0 = ManifoldState(M, [perturb(M, p0, perturb_var) for _ in t_vec])
        X1 = ManifoldState(M, [perturb(M, p1, perturb_var) for _ in t_vec])

        #Independently draw endpoint-conditioned samples at times t_vec:
        Xt = endpoint_conditioned_sample(X0, X1, P, t_vec)

        #Plot the torus:
        R = 2
        r = 0.5
        u = range(0, 2π; length=100)  # angle around the tube
        v = range(0, 2π; length=100)  # angle around the torus center
        pl = plot([(R + r*cos(θ))*cos(φ) for θ in u, φ in v], [(R + r*cos(θ))*sin(φ) for θ in u, φ in v], [r*sin(θ) for θ in u, φ in v],
            color = "grey", alpha = 0.3, label = :none, camera = (30,30))

        #Map the points to 3D and plot them:
        endpts = stack(tor.([p0,p1]))
        smppts = stack(tor.(eachcol(tensor(Xt))))
        scatter!(smppts[1,:], smppts[2,:], smppts[3,:], label = :none, msw = 0, ms = 1.5, color = "blue", alpha = 0.5)
        scatter!(endpts[1,:], endpts[2,:], endpts[3,:], label = :none, msw = 0, ms = 2.5, color = "red")
        savefig("torus_$(perturb_var)_$(P.v).svg")
    end
end

torus_0.0_0.svg:

Image

torus_0.0001_0.svg:

Image

torus_0.0_0.05.svg:

Image

torus_0.0001_0.05.svg:

Image

API Reference

ForwardBackward.BrownianMotion — Type
BrownianMotion(δ::T1, v::T2) where T1 <: Real where T2 <: Real
BrownianMotion(v::Real)
BrownianMotion()

Brownian motion process with drift δ and variance v.

Parameters

  • δ: Drift parameter (default: 0)
  • v: Variance parameter (default: 1)

Tip

The rate parameters must match the type of the state (eg. Float32 both both process and state), but you can avoid this if you use integer process parameters.

Examples

# Standard Brownian motion
process = BrownianMotion()

# Brownian motion with drift 0.5 and variance 2.0
process = BrownianMotion(0.5, 2.0)
source
ForwardBackward.CategoricalLikelihood — Type
CategoricalLikelihood(dist::AbstractArray, log_norm_const::AbstractArray)
CategoricalLikelihood(K::Int, dims...; T=Float64)
CategoricalLikelihood(dist::AbstractArray)

Probability distribution over discrete states.

Parameters

  • dist: Probability masses for each state
  • log_norm_const: Log normalization constants
  • K: Number of categories (for initialization)
  • dims: Additional dimensions for initialization
  • T: Numeric type (default: Float64)
source
ForwardBackward.ContinuousState — Type
ContinuousState(state::AbstractArray{<:Real})

Representation of continuous states.

Parameters

  • state: Array of current state values

Examples

# Create a continuous state
state = ContinuousState(randn(100))
source
ForwardBackward.GaussianLikelihood — Type
GaussianLikelihood(mu::AbstractArray, var::AbstractArray, log_norm_const::AbstractArray)

Gaussian probability distribution over continuous states.

Parameters

  • mu: Mean values
  • var: Variances
  • log_norm_const: Log normalization constants
source
ForwardBackward.ManifoldProcess — Type
ManifoldProcess(v::T)

A stochastic process on a Riemannian manifold with drift variance v.

Parameters

  • v: Drift variance parameter (default: 0)
source
ForwardBackward.ManifoldState — Type
ManifoldState(M::AbstractManifold, state::AbstractArray{<:AbstractArray})

Represents a state on a Riemannian manifold.

Parameters

  • state: Array of points on the manifold
  • M: The manifold object
ManifoldState(M::ProbabilitySimplex, x::AbstractArray{<:Integer}; softner! = ForwardBackward.soften!)

Convert a discrete array to points on a probability simplex. maximum(x) must be <= manifold_dimension(M)+1. By default this moves the points slightly away from the corners of the simplex (see soften!).

source
ForwardBackward.OUBridgeDistVarSched — Type
struct OUBridgeDistVarSched <: ContinuousProcess

OU-style endpoint-conditioned bridge with drift θ * (xc - Xt) and a variance-to-go schedule defined by a distribution in absolute time.

Let S(t) = 1 - cdf(vdist, t) (survival). Define R(t) = Rscale * S(t). For a < b < c, with Xa = xa and Xc = xc: phi(t2,t1) = exp(-θ * (t2 - t1)) r = S(b) / S(a) μb = xc + phi(b,a) * r * (xa - xc) Var_b = (Rscale * S(b) * (S(a) - S(b))) / (S(a) * phi(c,b)^2)

This formulation uses only survival ratios (no clamping), so it’s partition-consistent. Require S(a) > 0 (i.e., cdf(vdist, a) < 1).

source
ForwardBackward.OUBridgeExpVar — Type
struct OUBridgeExpVar <: ContinuousProcess

Endpoint-conditioned OU bridge with time-varying noise.

  • Drift of the conditional SDE is Λ(t)*(xc - Xt), with Λ including θ and the bridge term.
  • Noise schedule uses the same mixture-of-exponentials parameterization as your OU: σ²(t) ≈ a0 + sum_k w[k] * exp(β[k] * t).
  • the endpoint x_c provides the target, acting like the mean in an OU process.
  • Only endpoint-conditioned single-time sampling is provided.

Parameters

  • θ :: Real # base endpoint-reversion knob (mixing-rate)
  • a0 :: Real # baseline noise scale
  • w :: AbstractVector{<:Real} # weights for exp components
  • β :: AbstractVector{<:Real} # exponents (can be negative/positive)

Notes

  • Sampling at time b uses closed-form Gaussian conditioning with: φ(t,s) = exp(-θ(t - s)) Vb = _ounoiseQ(a, b, θ, a0, w, β) Vc = ounoiseQ(a, c, θ, a0, w, β) C{b,c} = exp(-θ(c - b)) * Vb mb = xc + φ(b,a)*(xa - xc), mc = xc + φ(c,a)*(xa - x_c)
  • No numerical integration; uses your _ou_noise_Q.
source
ForwardBackward.OrnsteinUhlenbeck — Type
OrnsteinUhlenbeck(μ::T1, v::T2, θ::T3) where T1 <: Real where T2 <: Real where T3 <: Real
OrnsteinUhlenbeck()

Ornstein-Uhlenbeck process with mean μ, variance v, and mean reversion rate θ.

Parameters

  • μ: Long-term mean (default: 0)
  • v: Variance parameter (default: 1)
  • θ: Mean reversion rate (default: 1)

Tip

The rate parameters must match the type of the state (eg. Float32 both both process and state), but you can avoid this if you use integer process parameters.

Examples

# Standard OU process
process = OrnsteinUhlenbeck()

# OU process with custom parameters
process = OrnsteinUhlenbeck(1.0, 0.5, 2.0)
source
ForwardBackward.OrnsteinUhlenbeckExpVar — Type
OrnsteinUhlenbeckExpVar(μ, θ::Real, a0::Real, w::AbstractVector{<:Real}, β::AbstractVector{<:Real})
OrnsteinUhlenbeckExpVar()
OrnsteinUhlenbeckExpVar(μ, θ, v)
OrnsteinUhlenbeckExpVar(μ, θ, v_at_0, v_at_1; dec = -0.1)

Ornstein–Uhlenbeck process with time-varying instantaneous variance v(t) = a0 + sum(w[k] * exp(β[k] * t) for k).

Parameters

  • μ : Long-run mean (default 0)
  • θ : Mean-reversion rate (default 1)
  • a0 : Baseline variance level (default 1)
  • w : Weights of exponential components (default empty)
  • β : Exponents (same length as w) (default empty)

Notes

  • Keep w ≥ 0 and a0 ≥ 0 if you want v(t) ≥ 0.
  • Handles the limits θ → 0 and β[k] + 2θ → 0 in a numerically stable way.
  • OrnsteinUhlenbeckExpVar(μ, θ, vat0, vat1; dec = -0.1) sets up a process where the variance decays nearly linearly from vat0 to vat1 over the interval [0, 1].
source
ForwardBackward.PiQ — Type
PiQ(r::Real, π::Vector{<:Real}; normalize=true)
PiQ(π::Vector{<:Real}; normalize=true)

Discrete process that switches to states proportionally to the stationary distribution π with rate r.

Parameters

  • r: Overall switching rate (default: 1.0)
  • π: Target stationary distribution (will always be normalized to sum to 1)
  • normalize: Whether to normalize the expected substitutions per unit time to be 1 when r = 1 (default: true)

Examples

# Process with uniform stationary distribution
process = PiQ(ones(4) ./ 4)

# Process with custom stationary distribution and rate
process = PiQ(2.0, [0.1, 0.2, 0.3, 0.4])
source
ForwardBackward.ScheduledGaussianMixture — Type
struct ScheduledGaussianMixture <: ContinuousProcess

Gaussian bridge for single-time endpoint-conditioned sampling.

Design:

  • Drift of the conditional SDE is Λ(t) * (xc - Xt).
  • Two schedules, specified in absolute time (no renormalization): • Mean approach via an integrable rate k(t) with integral ∫{t1}^{t2} k = k0(t2 - t1) + k1(Fk(t2) - Fk(t1)), where Fk is the CDF of kdist. • Variance-to-go R(t) via a survival function R(t) = Rscale * Sv(t), Sv(t) = 1 - Fv(t), where Fv is the CDF of vdist. Then R(c)=0.

Bridge at time b (a < b < c), with Xa=xa and Xc=xc:

  • φ(t2,t1) = exp(-∫_{t1}^{t2} k).
  • mean: μb = xc + φ(b,a) * (R(b)/R(a)) * (xa - xc)
  • variance: Var_b = R(b) * (R(a) - R(b)) / ( R(a) * φ(c,b)^2 )

Absolute-time CDFs ensure consistency across partitions: sampling (a→b1→b2) matches sampling (a→b2) marginally.

Fields:

  • k0, k1 : Real
  • kdist : UnivariateDistribution (absolute-time CDF for k)
  • Rscale : Real (scales R)
  • vdist : UnivariateDistribution (absolute-time CDF for R)
source
ForwardBackward.UniformDiscrete — Type
UniformDiscrete(μ::Real)
UniformDiscrete()

Discrete process with uniform transition rates between states, scaled by μ.

Parameters

  • μ: Rate scaling parameter (default: 1)

Tip

The rate parameters must match the type of the state (eg. Float32 both both process and state), but you can avoid this if you use integer process parameters.

source
ForwardBackward.UniformUnmasking — Type
UniformUnmasking(μ::Real)
UniformUnmasking()

Mutates only a mask (the last state index) to any other state (with equal rates). When everything is masked, μ=1 corresponds to one substitution per unit time.

Parameters

  • μ: Rate parameter (default: 1)

Tip

The rate parameters must match the type of the state (eg. Float32 both both process and state), but you can avoid this if you use integer process parameters.

source
ForwardBackward.:⊙ — Method
⊙(a::CategoricalLikelihood, b::CategoricalLikelihood; norm=true)
⊙(a::GaussianLikelihood, b::GaussianLikelihood)

Compute the pointwise product of two likelihood distributions. For Gaussian likelihoods, this results in another Gaussian. For categorical likelihoods, this results in another categorical distribution.

Parameters

  • a, b: Input likelihood distributions
  • norm: Whether to normalize the result (categorical only, default: true)

Returns

A new likelihood distribution of the same type as the inputs.

source
ForwardBackward.backward! — Method
backward!(Xdest::StateLikelihood, Xt::State, process::Process, t)
backward!(Xdest, Xt, process::Process, t1, t2) = backward!(Xdest, Xt, process, t2 - t1) #For time-homogeneous processes
backward(Xt::StateLikelihood, process::Process, t)
backward(Xt::State, process::Process, t)
backward(Xt, process::Process, t1, t2) = backward!(Xt, process, t2 - t1) #For time-homogeneous processes
backward(Xt::StateLikelihood, process::Process, t1, t2)
backward(Xt::State, process::Process, t1, t2)

Propagate a state or likelihood backward in time according to the process dynamics.

Parameters

  • Xdest: Destination for in-place operation
  • Xt: Final state or likelihood
  • process: The stochastic process
  • t: Time to propagate backward, for the single-time call
  • t1, t2: Start and end times, for the two-time call

Returns

The backward-propagated state or likelihood

source
ForwardBackward.endpoint_conditioned_sample — Method
endpoint_conditioned_sample(Xa, Xc, P::Process, t)
endpoint_conditioned_sample(Xa, Xc, P::Process, tF, tB)
endpoint_conditioned_sample(Xa, Xc, P::Process, t_a, t_b, t_c)
endpoint_conditioned_sample(Xa, Xc, P::Deterministic, tF, tB)

Generate a sample from the endpoint-conditioned process.

Parameters

  • Xa: Initial state
  • Xc: Final state
  • P: The stochastic process

Time argumenrs

  • t: For single-time call, this samples at time=t assuming endpoints at time=0 and time=1.
  • tF, tB: For two-time call, this assumes tF is the forward time and tB is the backward time (allowed for time-homogeneous processes)
  • t_a, t_b, t_c: If 3-time call, this samples at time=tb assuming endpoints at time=ta and time=t_c.

Returns

A sample from the endpoint-conditioned distribution

Notes

For continuous processes, uses the forward-backward algorithm. For deterministic processes, uses linear interpolation.

source
ForwardBackward.endpoint_conditioned_sample — Method
endpoint_conditioned_sample(Xa::ContinuousState, Xc::ContinuousState,
                            P::OUBridgeDistVarSched, t_a, t_b, t_c) :: ContinuousState

Draw Xb | (Xa = Xa.state, X_c = Xc.state), broadcasting over array states and scalar/array times via expand. Assumes Distributions is available (uses cdf). Throws if S(a) = 1 - cdf(vdist, a) ≤ 0 anywhere.

source
ForwardBackward.endpoint_conditioned_sample — Method
endpoint_conditioned_sample(X0::ManifoldState, X1::ManifoldState, p::ManifoldProcess, tF, tB; Δt = 0.05)

Generate a sample from the endpoint-conditioned process on a manifold.

Parameters

  • X0: Initial state
  • X1: Final state
  • p: The manifold process
  • tF: Forward time
  • tB: Backward time
  • Δt: Discretized step size (default: 0.05)

Returns

A sample state at the specified time

Notes

Uses a numerical stepping procedure to approximate the endpoint-conditioned distribution.

source
ForwardBackward.forward! — Method
forward!(Xdest::StateLikelihood, Xt::State, process::Process, t)
forward!(Xdest, Xt, process::Process, t1, t2) = forward!(Xdest, Xt, process::Process, t2-t1) #For time-homogeneous processes
forward(Xt::StateLikelihood, process::Process, t)
forward(Xt::State, process::Process, t)
forward(Xt, process::Process, t1, t2) = forward!(Xt, process, t2 - t1) #For time-homogeneous processes
forward(Xt::StateLikelihood, process::Process, t1, t2)
forward(Xt::State, process::Process, t1, t2)

Propagate a state or likelihood forward in time according to the process dynamics.

Parameters

  • Xdest: Destination for in-place operation
  • Xt: Initial state or likelihood
  • process: The stochastic process
  • t: Time to propagate forward, for the single-time call
  • t1, t2: Start and end times, for the two-time call

Returns

The forward-propagated state or likelihood.

source
ForwardBackward.interpolate — Method
interpolate(X0::ContinuousState, X1::ContinuousState, tF, tB)
interpolate(X_a::ContinuousState, X_c::ContinuousState, t_a, t_b, t_c) = interpolate(X_a, X_c, t_b .- t_a, t_c .- t_b)

Linearly interpolate between two continuous states.

Parameters

  • X0: Initial state
  • X1: Final state
  • tF: Forward time
  • tB: Backward time
  • t_a, t_b, t_c: If 3-time call, this assumes t_b-t_a is the forward time and t_c-t_b is the backward time.

Returns

The interpolated state

source
ForwardBackward.interpolate — Method
interpolate(X0::ManifoldState, Xt::ManifoldState, tF, tB)

Interpolate between two states on a manifold using geodesics.

Parameters

  • dest: Destination state for in-place operation
  • X0: Initial state
  • Xt: Final state
  • tF: Time difference from initial state
  • tB: Time difference to final state

Returns

The interpolated state

source
ForwardBackward.perturb! — Method
perturb!(M::AbstractManifold, q, p, v)

Perturb a point p on manifold M by sampling from a normal distribution in the tangent space with variance v and exponentiating back to the manifold.

Parameters

  • M: The manifold
  • q: The point that is overwritten (for perturb!)
  • p: Original point
  • v: Variance of perturbation
source
ForwardBackward.soften! — Method
soften!(x::AbstractArray{T}, a = T(1e-5)) where T

In-place regularizes values of x slightly while preserving its sum along the first dimension.

(x .+ a) ./ sum(x .+ a, dims = 1)
source
ForwardBackward.stochastic — Method
stochastic(o::State)
stochastic(T::Type, o::State)

Convert a state to its corresponding likelihood distribution: A zero-variance (ie. delta function) Gaussian for the continuous case, and a one-hot categorical distribution for the discrete case.

Parameters

  • o: Input state
  • T: Numeric type for the resulting distribution (default: Float64)

Returns

A likelihood distribution corresponding to the input state.

source
ForwardBackward.tensor — Method
tensor(d::Union{State, StateLikelihood})

Convert a state or likelihood to its tensor (ie. multidimensional array) representation.

Returns

The underlying array representation of the state or likelihood.

source

```