Composing distributions

Introduction

ComposedDistributions.jl composes per-event delay distributions into one object that describes a whole record: named events linked by delays, wired into a tree.

What we do here

Each section is a small runnable example rather than a full analysis. We: 2. Compose a record from per-event delays with compose, starting from plain Distributions.jl leaves.

3. Build the five composers directly (Sequential, Parallel, Resolve, Compete, Choose) and see how they nest.

4. Score and simulate from one composed object.

5. Attach parameters and priors with params_table and build_priors.

6. Edit an assembled tree with update, prune and splice.

The operator map

Every operator here falls into one of three families: structural composers wire branches into a tree, combinators add or difference whole delays, and the introspection verbs read or edit an assembled object. See the Concepts verb map for the full list, grouped the same way.

Packages used

We use Distributions for the delay distributions and Random for reproducibility.

using ComposedDistributions
using ComposedDistributions: update
using ConvolvedDistributions
using Distributions
using Random

Composing a record

A record is a set of events linked by delays. compose is the front-end: it takes a friendly description and lowers it to a nested stack of the composers, without introducing a new tree type. A NamedTuple of bare distributions names each branch off one shared origin (a Parallel); a Vector value is a chain of steps (a Sequential).

onset_admit = LogNormal(1.5, 0.4);

admit_death = Gamma(2.0, 1.0);

Two branches off one onset: an onset-to-admission delay and an onset-to-notification delay.

parallel_stack = compose((onset_admit = onset_admit,
    onset_notif = Gamma(1.5, 1.0)));

event_names(parallel_stack)

(:onset, :admit, :notif)

A composed stack is a distribution like any leaf. It simulates a named event record with rand,

rand(Xoshiro(1), parallel_stack)

(onset_admit = 4.356925919792406, onset_notif = 1.8400307468613506)

and scores one with logpdf. The whole toolkit below is built on these two operations.

logpdf(parallel_stack, rand(Xoshiro(1), parallel_stack))

-2.89126290386135

A chain step is a Vector: onset to admission, then admission to death.

chain = compose((path = [onset_admit, admit_death],));

event_names(chain)

(:event_1, :event_2, :event_3)

The same stack also lowers from an explicit Tables.jl table, so a column-oriented data source builds a composer without a hand-written NamedTuple. A table has name and dist columns, one row per branch. A chain column folds rows sharing a non-zero id into a Sequential, and a compete/prob column pair folds rows into a Resolve node whose prob entries are the branch probabilities.

table = [
    (name = :death, dist = Gamma(1.5, 1.0), compete = 1, prob = 0.3),
    (name = :discharge, dist = Gamma(2.0, 1.5), compete = 1, prob = 0.7),
    (name = :onset_notif, dist = Gamma(1.5, 1.0), compete = 0,
        prob = missing)];

table_stack = compose(table);

event_names(table_stack)

(:event_1, :death, :discharge, :notif)

The five composers

Each front-end lowers to these composers, which can also be built directly. They differ in how the branches relate.

Sequential is a conjunctive chain: each step adds an independent delay onto the previous event. The lowercase sequential verb is the public constructor; name the steps with name => dist pairs.

seq = sequential(:onset_admit => onset_admit, :admit_death => admit_death);

Parallel places independent branches off one shared origin, built with the parallel verb.

par = parallel(:onset_admit => onset_admit, :onset_notif => admit_death);

Resolve is a disjunction where exactly one outcome occurs, governed by fixed branch probabilities that sum to one. A death-versus-discharge split makes the death probability the case-fatality ratio.

cfr = 0.3;

resolution = resolve(:death => (Gamma(1.5, 1.0), cfr),
    :discharge => (Gamma(2.0, 1.5), 1 - cfr));

The last outcome's probability may be omitted (a bare name => delay). It then takes the residual 1 - sum(of the others), so the discharge probability 1 - cfr need not be written out.

resolution_residual = resolve(:death => (Gamma(1.5, 1.0), cfr),
    :discharge => Gamma(2.0, 1.5));

Its marginal is the time to resolution regardless of which outcome occurs.

mean(resolution)

2.55

A Resolve node carries its own time-to-resolution event slot alongside the named per-outcome slots, so its flat event layout pairs that resolution slot (defaulting to :event_1) with the outcome names.

event_names(resolution)

(:event_1, :death, :discharge)

Compete is the racing-hazard sibling of Resolve, built with the compete verb from bare name => delay outcomes (no probabilities). The cause-specific delays race, the first to fire wins, and which outcome wins is coupled to when it fires. Where resolve takes the winning probability as a fixed parameter, compete derives it from the hazards: the marginal any-event time is min of the racing delays, with survival the product of the per-cause survivals.

racing = compete(:death => Gamma(1.5, 1.0), :discharge => Gamma(2.0, 1.5));

Competing outcomes works through resolve versus compete in full, including when to reach for each.

Choose is a data-selected disjunction: the alternatives are independent sub-models with different origins, and a data field picks which one applies to a record. Neither Parallel nor Resolve (both shared origin) expresses this.

selector = choose(:index => onset_admit, :sourced => admit_death);

Scoring names the active alternative through the kind keyword.

logpdf(selector, 3.0; kind = :index)

-1.6047353862744176

Nesting

The composers nest, so trees of arbitrary depth are built by composing on composers. A compose result drops into another compose as a branch.

early = compose((onset_admit = onset_admit, onset_notif = admit_death));

nested = compose((early = early, late = chain));

event_names(nested)

(:onset, :admit, :notif, :event_1, :event_2)

A pre-built composer is a valid Sequential step, so a chain can carry a Resolve resolution as its terminal step. Naming the chain steps gives the simulated record readable event names.

tree = compose((
    path = sequential(:onset_admit => onset_admit,
        :admit_resolve => resolution),
    onset_notif = admit_death));

The flat event layout of a tree is derived from the edge names.

event_names(tree)

(:onset, :admit, :death, :discharge, :notif)

Scoring and simulation from one object

The composer is dual-purpose: it scores observed records and simulates new ones. A rand of a tree returns a full named event record, and logpdf reads one straight back.

record = rand(Xoshiro(7), tree)

(path_onset_admit = 3.7632484408708633, path_admit_resolve = 0.477544094272877, onset_notif = 0.8838923479842206)

The labelled draw round-trips through logpdf, either as the record or as its flat vector of values.

(from_record = logpdf(tree, record),
    from_vector = logpdf(tree, collect(values(record))))

(from_record = -3.804342220775429, from_vector = -3.804342220775429)

A Resolve node scores and samples its marginal time to resolution directly.

logpdf(resolution, 3.0)

-1.8610823128872573

rand(node; outcome = true) draws which outcome occurs and its time as a compact (outcome, time) pair, so a standalone draw tells you which outcome won.

rand(Xoshiro(7), resolution; outcome = true)

(:death, 4.894059791798188)

Combining whole delays

Where the composers wire named branches into a tree, the combinators join two whole delays algebraically. convolved forms the sum X + Y (the total of two independent delays), and difference forms the dual X - Y. Both are ConvolvedDistributions' own verbs, extended here for composed tree operands (a bare distribution works unchanged too).

total = convolved(Gamma(2.0, 1.0), LogNormal(0.5, 0.4));

(convolved_mean = mean(total),
    summed_mean = mean(Gamma(2.0, 1.0)) + mean(LogNormal(0.5, 0.4)))

(convolved_mean = 3.7860384307500734, summed_mean = 3.7860384307500734)

Reading the composed marginal

A composed chain has a marginal delay from its origin to its final event. The moments are additive over the steps, so mean and var of a Sequential sum the per-step moments.

chain_moments = sequential(:onset_admit => Gamma(2.0, 1.0),
    :admit_death => LogNormal(0.5, 0.4));

mean(chain_moments), var(chain_moments)

(3.7860384307500734, 2.553488101144678)

observed_distribution collapses that chain to the single convolved distribution of its origin-to-final gap, integrating the intermediate event out. Its mean matches the chain's overall mean.

collapsed = observed_distribution(chain_moments);

(collapsed_mean = mean(collapsed), chain_mean = mean(chain_moments))

(collapsed_mean = 3.7860384307500734, chain_mean = 3.7860384307500734)

Parameters and priors

A composed distribution carries a flat inventory of its free parameters. params_table lists one row per scalar parameter, keyed by the edge path and the parameter name, with the support a prior must respect. It prints as a table and is a Tables.jl source, so tbl.edge / tbl.param read its columns.

template = compose((onset_admit = Gamma(2.0, 1.0),
    admit_death = LogNormal(0.5, 0.4)));

tbl = params_table(template)

params_table (4 rows)
  edge         param  value  support     prior
  ───────────  ─────  ─────  ──────────  ─────
  onset_admit  shape  2.0    (0.0, Inf)       
  onset_admit  scale  1.0    (0.0, Inf)       
  admit_death  mu     0.5    (0.0, Inf)       
  admit_death  sigma  0.4    (0.0, Inf)

Its columns are accessed by name.

tbl.edge, tbl.param

([:onset_admit, :onset_admit, :admit_death, :admit_death], [:shape, :scale, :mu, :sigma])

build_priors takes that table (any Tables.jl source with edge, param, value, support columns) and derives a default prior per row from that leaf's support: a positive scale parameter gets a positive-truncated prior, a location parameter an unbounded one, a [0, 1] probability a Uniform(0, 1). So build_priors(tbl) alone yields a complete set, defined against the table rather than by hand-matching the tree.

priors = build_priors(tbl);

priors.onset_admit.shape

Truncated(Distributions.Normal{Float64}(μ=2.0, σ=2.0); lower=0.0)

Editing a composed tree

update applies a set of parameter values back to a composed object, returning a distribution of the same structure.

updated = update(template, (onset_admit = (shape = 3.0, scale = 1.5),
    admit_death = (mu = 0.7, sigma = 0.5)));

mean(updated)

(onset_admit = 4.5, admit_death = 2.2818807653293036)

update also replaces whole nodes, not just their values. Passing path => new_node swaps the node at an address for a new distribution, keeping the tree shape. The address is the same one event reads: a bare name, a dotted Symbol, or a tuple of edge names.

replaced = update(template, :admit_death => Gamma(3.0, 1.5));

event(replaced, :admit_death)

Distributions.Gamma{Float64}(α=3.0, θ=1.5)

Two edits that change the tree shape are kept separate. prune drops a branch from a node (renormalising a Resolve arm's remaining probabilities), and splice inserts a step around a node.

three_way = resolve(:death => (Gamma(1.5, 1.0), 0.3),
    :discharge => (Gamma(2.0, 1.5), 0.4),
    :transfer => (Gamma(1.0, 1.0), 0.3));

resolution_tree = compose((resolution = three_way, onset = Gamma(1.0, 1.0)));

pruned = prune(resolution_tree, :resolution, :transfer);

The :transfer outcome is gone from the pruned node's event layout.

event_names(event(pruned, :resolution))

(:event_1, :death, :discharge)

splice wraps a node in a chain, here adding a reporting delay after the death branch.

spliced = splice(template, :admit_death;
    after = :death_report => Gamma(1.0, 2.0));

event_names(event(spliced, :admit_death))

(:admit, :death, :report)

tie links leaves at several paths into one free parameter group, so params_table inventories the tied occurrences once under the shared tag rather than as separate parameters.

shared_rate = compose((a = Gamma(2.0, 1.0), b = Gamma(2.0, 1.0)));

tied = tie(shared_rate, :a, :b; name = :rate);

unique(params_table(tied).edge)

1-element Vector{Symbol}:
 :rate

Syntax reference

Every public composition form on one object, with whether it preserves the tree shape.

Syntax	What it does	Shape
compose((a = d1, b = d2))	NamedTuple front-end; a Vector value is a chain	builds
compose(table)	Tables.jl name/dist source; a chain column folds rows into a Sequential, a compete/prob pair into a Resolve	builds
sequential(:a => d1, :b => d2)	a Sequential chain (steps add up)	builds
parallel(:a => d1, :b => d2)	a Parallel branch set (shared origin)	builds
resolve(:a => (d1, p1), :b => (d2, p2))	a Resolve node; the last prob may be omitted as the residual	builds
compete(:a => d1, :b => d2)	a Compete racing-hazard node (the winning probability derived from the hazards)	builds
choose(:a => d1, :b => d2)	a Choose disjunction (data picks the branch)	builds
convolved(d1, d2)	a Convolved sum X + Y	builds
difference(d1, d2)	a Difference X - Y	builds
shared(:tag, d)	tag a leaf as a tied parameter group	leaf wrap
tie(d, paths...; name)	tie leaves at paths into one group	yes
update(d, (a = (shape = 3,),))	replace free parameter values	yes
update(d, path => new_node)	replace a whole node	yes
prune(d, path...)	drop a branch (renormalise a Resolve arm)	no (topology)
splice(d, path; before, after)	insert a step at a node	no (topology)
event(d, path...)	fetch a child or descend a name path	read
event_tree(d)	the nested tree of event names	read
event_names(d)	the per-record key names	read
observed_distribution(d)	collapse a chain to its convolved total	read
params_table(d)	the flat free-parameter inventory	read

The address path in event / update / prune / splice / tie is the same in all: a bare Symbol, a dotted Symbol (:a.b), or a tuple of edge names. shared(:tag, d) and tie(d, paths...; name = :tag) are two spellings of the same tie: shared tags a leaf where it is built, tie walks the tree to the named leaves and wraps each in the same tie. Both make the tagged occurrences one free parameter.

Summary

• compose lowers a NamedTuple, table, or matrix to the same composer stack.

• Sequential, Parallel, Resolve, Compete and Choose are conjunctive chains, shared-origin branches, fixed-probability disjunctions, racing-hazard disjunctions and data-selected disjunctions.

• The composers nest, including a composer as a chain step.

• One object scores records and simulates them: logpdf reads a record, rand generates one.

• convolved and difference combine two whole delays algebraically.

• mean and var read the composed marginal moments, and observed_distribution collapses a chain to its convolved total.

• params_table and build_priors attach parameters and support-derived priors to the same object.

• update edits the tree: path => new_node replaces nodes keeping the shape, prune and splice are the two topology edits, and tie links leaves into one parameter group.

Where next

• The Public API lists every composer, combinator and introspection verb with its full docstring.