ConvolvedDistributions.jl

Appearance

1

Getting started ​

ConvolvedDistributions builds the distribution of a sum (X + Y, a convolution), a signed gap (X - Y), a product (X * Y), or a ratio (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.

julia
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.

julia
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.

julia
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.

julia
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 Convolving distributions 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.

julia
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.

julia
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:

julia
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.

julia
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.

julia
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.

julia
scaled = product(d, LogNormal(0.0, 0.2))
mean(scaled), mean(d)
(3.8625214804440207, 3.7860384307500734)

Ratios ​

ratio builds the quotient member: the distribution of Z = X / Y for independent X and Y, as when a rate, a proportion, or a normalised measurement is formed from two independent uncertain quantities. Normal(0, σx) / Normal(0, σy) uses the closed form (a Cauchy); Gamma / Gamma and Chisq / Chisq also have closed forms; everything else uses the numeric branch-split quadrature. Unlike product, either component may have two-sided support — the denominator only needs to avoid probability mass at zero.

julia
r = ratio(Gamma(3.0, 1.0), Gamma(2.0, 1.0))
mean(r), cdf(r, 2.0)
(3.0, 0.5925925925925924)

The sign-crossing headline pair — a Ratio of two zero-mean Normals — matches its Cauchy closed form exactly, including when the numerator or denominator itself straddles zero.

julia
rc = ratio(Normal(0.0, 2.0), Normal(0.0, 0.5))
cdf(rc, 0.0), cdf(Cauchy(0.0, 4.0), 0.0)
(0.5, 0.5)

mean/var/std throw unless an analytic pair applies: E[X / Y] = E[X] E[1/Y] needs an inverse moment of the denominator, which this package does not compute and which need not exist.

A Ratio can itself be a component of another combination only when both its numerator and its denominator are non-negative: outside that regime the ratio's tails are Cauchy-like, so nesting it throws rather than silently narrowing the outer window (see the FAQ).

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)):

julia
t = 0:30
infections = 100 .* exp.(-((t .- 10.0) .^ 2) ./ 40.0)
convolve_series(Poisson(2.0), infections)
31-element Vector{Float64}:
  1.1108996538242308
  4.008141722979866
  8.526856383040526
 14.494207016175757
 22.040612466668186
 31.330073874163926
 42.26184773370586
 54.34073062309471
 66.68501861759444
 78.13217137485263
  ⋮
  9.919612593388146
  6.208280252442309
  3.7216254767162114
  2.1374312845880192
  1.1764298210822512
  0.6206879051664294
  0.31400439684842807
  0.1523624271040959
  0.07092950565089116

See CensoredDistributions.jl for tools that turn a continuous delay into the masses this takes.

julia
maxlag = length(infections) - 1
pmf = pdf.(NegativeBinomial(5, 0.5), 0:maxlag)
convolve_series(pmf, infections)
31-element Vector{Float64}:
  0.25651562069968376
  1.0537698117111787
  2.6240620962620347
  5.164368050762789
  8.858346206111577
 13.866428019808481
 20.283325643091214
 28.078676924114166
 37.04218347085685
 46.75433145856008
  ⋮
 33.807455928780485
 26.539865161903737
 20.325968118376395
 15.207455414959044
 11.130079097937406
  7.979069966780508
  5.610250749929932
  3.8737836505860237
  2.6298948485084055
julia
sum(pmf)
0.9999982672743509

A delay that changes over the window is passed as one delay per time point; the Convolving a timeseries tutorial works through it.

julia
delays = [Poisson(λ) for λ in range(3.0, 1.0; length = length(infections))]
convolve_series(delays, infections)
31-element Vector{Float64}:
  0.40867714384640674
  1.928494010662295
  5.048157769762642
  9.940060355409067
 16.698109411262777
 25.44938016082706
 36.22138116457638
 48.72961625807533
 62.242871888035765
 75.58954753940728
  ⋮
  9.368099242258317
  5.54931436317456
  3.1156879489771225
  1.6572189336439593
  0.8346002454646066
  0.3977269735514724
  0.17922764642080175
  0.07631500112318403
  0.03067836891767108

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.

julia
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, Product, and Ratio 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.

julia
td = truncated(d, 0.0, 8.0)
logpdf(td, 5.0)
-1.970917408577644

Distributions.censored composes the same way, clamping rather than renormalising: the kept region's density is unchanged and the trimmed tails become point masses at the bounds. This package adds no censoring machinery of its own — that stays CensoredDistributions.jl's job — but the generic wrapper works out of the box because Convolved, Difference, Product, and Ratio implement the standard UnivariateDistribution interface.

julia
cd = censored(d, 0.0, 8.0)
logpdf(cd, 5.0)
-1.9898782016345686

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.

julia
using Optimization, OptimizationOptimJL

quantile(d, 0.5)     # median by inverse-CDF root-find
3.5052589047225053
julia
length(rand(truncated(d, 0.0, 8.0), 100))
100

Nothing else on this page needs the extension. rand on a bare Convolved, Difference, Product, or Ratio 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 ​