Getting started

ConvolvedDistributions builds the distribution of a sum (X + Y, a convolution), a signed gap (X - Y), or a product (X * Y) of independent random variables, for any pair of Distributions.jl univariate distributions. Closed forms are used where they exist and an AD-safe Gauss-Legendre quadrature everywhere else, so the results can be scored, truncated, and differentiated inside a fitting loop. This page walks through the main entry points; the Public API has the full interface.

See Installation for how to install the package.

Convolving distributions

convolved returns the distribution of the sum of independent components. For pairs with a known closed form (Normal + Normal, equal-scale Gamma, equal-rate Exponential) it delegates to the analytic result; any other pair uses the numeric quadrature fallback.

using ConvolvedDistributions, Distributions

# An incubation period plus a reporting delay, say.
d = convolved(Gamma(2.0, 1.0), LogNormal(0.5, 0.4))

cdf(d, 5.0)

0.8051810060082844

Densities and moments work the same way; the mean and variance are exact component sums, not quadrature results.

pdf(d, 5.0), mean(d), var(d)

(0.13671207573887373, 3.7860384307500734, 2.553488101144678)

More than two components can be passed as varargs, a tuple, or a vector.

d3 = convolved(Gamma(2.0, 1.0), LogNormal(0.5, 0.4),
    Exponential(2.0))
mean(d3)

5.786038430750073

Composing convolutions

The result of convolved is itself a UnivariateDistribution, so it can be a component of another call. Here d (the incubation plus reporting delay from above) is convolved with a further processing delay, and the nested distribution matches the flat three-component one.

total = convolved(d, Exponential(2.0))
mean(total), mean(d3), cdf(total, 8.0), cdf(d3, 8.0)

(5.786038430750073, 5.786038430750073, 0.8304533961544223, 0.8304533961544223)

Nesting is how a multi-stage delay built in one part of a model is reused in another without unpacking its components. The Visualising convolutions tutorial plots a nested convolution against its flat equivalent.

Evaluating cdf or pdf over a vector of points shares one quadrature window solve across the batch, which is much cheaper than mapping the scalar call.

cdf(d, [1.0, 2.5, 5.0, 10.0])

4-element Vector{Float64}:
 0.001816868167571861
 0.21402716540748604
 0.8051810060082844
 0.9966058316148315

Differences

difference is the dual of the sum: the distribution of Z = X - Y for independent X and Y, with support on both sides of zero. Normal - Normal uses the closed form; everything else uses the numeric cross-correlation path.

z = difference(Gamma(3.0, 1.0), LogNormal(0.5, 0.4))
cdf(z, 0.0)

0.2683441902952708

A symmetric difference is centred on zero:

zs = difference(Normal(1.0, 1.0), Normal(1.0, 1.0))
mean(zs), cdf(zs, 0.0)

(0.0, 0.5)

Either side can itself be a Convolved, so the signed gap between a multi-stage delay and a single delay is one call. Here cdf(gap, 0.0) is the probability that the two-stage delay d is shorter than the single delay.

gap = difference(d, Gamma(2.5, 1.0))
mean(gap), cdf(gap, 0.0)

(1.2860384307500734, 0.25861733505565326)

Products

product builds the multiplicative member: the distribution of Z = X * Y for independent X and Y, as when a delay is stretched by an independent multiplicative factor. LogNormal * LogNormal uses the closed form (the log-parameters add); everything else uses the numeric Mellin quadrature. Both components must have non-negative support; sign-crossing supports throw an error and are future work.

w = product(Gamma(3.0, 1.0), LogNormal(0.0, 0.3))
mean(w), cdf(w, 3.0)

(3.1380835797261506, 0.5754383333307876)

The mean and variance are the exact independent-product moments, and a Convolved (or another combination) can itself be a component, so a multi-stage delay scaled by a factor is one call.

scaled = product(d, LogNormal(0.0, 0.2))
mean(scaled), mean(d)

(3.8625214804440207, 3.7860384307500734)

Convolving a timeseries

convolve_series causally convolves a numeric series with a delay PMF on the unit lag grid. If the series holds the expected events at times 0, 1, ..., t (say infections), the result is the expected downstream counts at the same times, the renewal-style observation layer.

A discrete delay is read straight off its own PMF (the lag-k mass is pdf(delay, k)):

infections = [0.0, 1.0, 3.0, 6.0, 8.0, 5.0, 2.0]
convolve_series(Poisson(2.0), infections)

7-element Vector{Float64}:
 0.0
 0.1353352832366127
 0.6766764161830635
 1.8946939653125778
 3.6991644084674142
 5.097629001912411
 5.178830171854379

A continuous delay carries no mass on the integer grid until it is discretised, and discretising it is an explicit modelling choice, so convolve_series will not do it silently. Discretise it first with discretise_pmf (raw CDF-difference masses: interval-censored secondary event, exact primary), then convolve. For a primary event known only to the day use CensoredDistributions.jl's double-interval-censored masses instead. The masses depend differentiably on the delay parameters, so this composes with gradient-based fitting.

pmf = discretise_pmf(d, length(infections) - 1)
convolve_series(pmf, infections)

7-element Vector{Float64}:
 0.0
 0.0018168681675718615
 0.09890781607000299
 0.5527470757071435
 1.626755465124929
 3.308129228410993
 4.810516004315176

Choosing the solver

Both constructors take a method keyword. The default AnalyticalSolver() uses the closed form when one exists and falls back to quadrature; NumericSolver() forces the quadrature path, which is mainly useful for testing and for comparing the two.

da = convolved(Normal(0.0, 1.0), Normal(1.0, 2.0))
dn = convolved(Normal(0.0, 1.0), Normal(1.0, 2.0);
    method = NumericSolver())
cdf(da, 2.0), cdf(dn, 2.0)

(0.6726395769907114, 0.6726395669907115)

Truncation and scoring

Convolved, Difference, and Product compose with Distributions.truncated, so right-truncated (or doubly truncated) scoring works out of the box. This is the usual pattern for fitting delay data observed up to a cut-off.

td = truncated(d, 0.0, 8.0)
logpdf(td, 5.0)

-1.970917408577644

Quantiles and sampling truncated distributions

There is no closed-form inverse CDF for a generic convolution, so quantile lives in an extension that is loaded when both Optimization.jl and OptimizationOptimJL.jl are present. It finds the quantile by a Nelder-Mead inversion of cdf. Loading it also enables rand on truncated wrappers, which routes through the base quantile.

using Optimization, OptimizationOptimJL

quantile(d, 0.5)     # median by inverse-CDF root-find

3.5052589047225053

length(rand(truncated(d, 0.0, 8.0), 100))

100

Nothing else on this page needs the extension. rand on a bare Convolved, Difference, or Product samples the components directly.

Gradients

The cdf, pdf, and logpdf paths are AD-safe by construction. The quadrature uses fixed nodes, the window clamp is shielded from the tape, and the gamma CDF carries analytic derivative rules (supplied by EpiAwareADTools.jl). Gradients with respect to the component parameters are tested on ForwardDiff, ReverseDiff, Enzyme (forward and reverse), and Mooncake (forward and reverse) on every CI run; the per-backend badges in the README track their status.

Learning more

• New to the package? Installation covers installing by URL and the optional quantile extension; for setting up Julia itself, see the EpiAware site.

• Common questions (solver choice, the timeseries form, the extension, AD support) are answered in the FAQ.

• Want the full interface? See the Public API.

• Curious how the numeric layer is put together? The internal quadrature (integrate, gl_integrate, GaussLegendre) is documented in the Internal API, and an Integrals.jl backend is available as an extension.

• Contributing, or adding a new member of the combination family? Start from the developer documentation, the Contributing guide, and Adding a new combination.

• Want to report a problem or ask a question? Open an issue or start a discussion on the GitHub repository.