Composing distributions
Introduction
ComposedDistributions.jl composes per-event delay distributions into one object that describes a whole record. The same object scores observed records and simulates new ones, so a model is built once and used in both directions.
What we do here
Each section is a small runnable example rather than a full analysis. We:
Compose a record from per-event delays with
compose, starting from plain Distributions.jl leaves.Build the five composers directly (
Sequential,Parallel,Resolve,Compete,Choose) and see how they nest.Score and simulate from one composed object.
Attach parameters and priors with
params_tableandbuild_priors.
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.
Structural composition (wire branches into a tree)
├─ compose lower a NamedTuple / table / matrix to the stack
├─ sequential conjunctive chain (steps add up)
├─ parallel independent branches off one shared origin
├─ resolve one outcome occurs by fixed probability
├─ compete racing hazards, first to fire wins
├─ choose a data field picks the branch
├─ tie tie leaves at paths into one parameter group
└─ shared tag a leaf as a tied parameter group at build time
Combination (add or difference whole delays)
├─ convolve_distributions the sum X + Y (Convolved)
└─ difference the dual X - Y (Difference)
Reading and editing (inspect or edit an assembled object)
├─ mean / var the composed marginal moments
├─ observed_distribution collapse a chain to its convolved total
├─ params_table the flat free-parameter inventory
├─ event / event_names fetch a child / the per-record key names
└─ update / prune / splice / tie edit values, topology or tiesPackages used
We use Distributions for the delay distributions and Random for reproducibility.
using ComposedDistributions
using Distributions
using RandomComposing 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.89126290386135A 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.55A 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));The marginal any-event time is the min of the racing delays, so its survival is the product of the per-cause survivals.
(racing_ccdf = ccdf(racing, 3.0),
product_ccdf = ccdf(Gamma(1.5, 1.0), 3.0) * ccdf(Gamma(2.0, 1.5), 3.0))(racing_ccdf = 0.04531440427588506, product_ccdf = 0.045314404275885074)Reach for compete when the outcome split is driven by competing risks on a shared clock, so the probabilities follow from the delays rather than being set; reach for resolve when the split is a fixed probability independent of timing; reach for choose (below) when a data field selects which sub-model a record follows.
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.6047353862744176Nesting
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.8610823128872573rand_outcome draws which outcome occurs and its time as a compact (outcome, time) pair, so a standalone draw tells you which outcome won.
ComposedDistributions.rand_outcome(Xoshiro(7), resolution)(:death, 4.894059791798188)Combining whole delays
Where the composers wire named branches into a tree, the combinators join two whole delays algebraically. convolve_distributions forms the sum X + Y (the total of two independent delays), and difference forms the dual X - Y. Both are re-exported from ConvolvedDistributions, so they are reachable through ComposedDistributions alone.
total = convolve_distributions(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.shapeTruncated(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}:
:rateSyntax 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 |
convolve_distributions(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
composelowers a NamedTuple, table, or matrix to the same composer stack.Sequential,Parallel,Resolve,CompeteandChooseare 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:
logpdfreads a record,randgenerates one.convolve_distributionsanddifferencecombine two whole delays algebraically.meanandvarread the composed marginal moments, andobserved_distributioncollapses a chain to its convolved total.params_tableandbuild_priorsattach parameters and support-derived priors to the same object.updateedits the tree:path => new_nodereplaces nodes keeping the shape,pruneandspliceare the two topology edits, andtielinks leaves into one parameter group.
Where next
- The Public API lists every composer, combinator and introspection verb with its full docstring.