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(δ::T, v::T) where T <: 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.OrnsteinUhlenbeck — Type
OrnsteinUhlenbeck(μ::T, v::T, θ::T) where T <: 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.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.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(Xt::StateLikelihood, process::Process, t)
backward(Xt::State, process::Process, t)

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

Returns

The backward-propagated state or likelihood

source
ForwardBackward.endpoint_conditioned_sample — Method
endpoint_conditioned_sample(X0, X1, p, tF, tB)
endpoint_conditioned_sample(X0, X1, p, t)
endpoint_conditioned_sample(X0, X1, p::Deterministic, tF, tB)

Generate a sample from the endpoint-conditioned process.

Parameters

  • X0: Initial state
  • X1: Final state
  • p: The stochastic process
  • t, tF: Forward time
  • tB: Backward time (defaults to 1-t for single time parameter)

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(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(Xt::StateLikelihood, process::Process, t)
forward(Xt::State, process::Process, t)

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

Returns

The forward-propagated state or likelihood

source
ForwardBackward.interpolate — Method
interpolate(X0::ContinuousState, X1::ContinuousState, tF, tB)

Linearly interpolate between two continuous states.

Parameters

  • X0: Initial state
  • X1: Final state
  • tF: Forward time
  • tB: 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

```