ConvolvedDistributions.jl

Appearance

Getting started ​

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

Installation ​

ConvolvedDistributions is not yet registered, so install it from GitHub:

julia
using Pkg
Pkg.add(url = "https://github.com/EpiAware/ConvolvedDistributions.jl")

Once registered this becomes Pkg.add("ConvolvedDistributions").

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

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.13671207573887384, 3.7860384307500734, 2.553488101144678)

More than two components can be passed as varargs, a tuple, or a vector, and Convolved distributions nest.

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

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.0018172341322499784
 0.21402687896669537
 0.8051810053750608
 0.9966058316148325

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

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)

Convolving a timeseries ​

convolve_series discretises the delay to a PMF over the unit grid and causally convolves a numeric series with it. 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. This is the renewal-style observation layer, and the PMF masses depend differentiably on the delay parameters, so it composes with gradient-based fitting.

julia
infections = [0.0, 1.0, 3.0, 6.0, 8.0, 5.0, 2.0]
convolve_series(d, infections)
7-element Vector{Float64}:
 0.0
 0.001816868167571864
 0.09890781607000303
 0.5527470757071438
 1.6267554651249296
 3.3081292284109938
 4.8105160043151765

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.6726395669907121)

Truncation and scoring ​

Convolved and Difference 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.9709174085776435

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.5051558447780256
julia
length(rand(truncated(d, 0.0, 8.0), 100))
100

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