A convolved distribution

Creates the distribution of the sum of two or more independent random variables using numerical convolution of their densities. The density of the result is approximated on a grid using a Fast Fourier Transform (FFT), and then interpolated with stats::approxfun().

The cumulative distribution function and quantile function are derived from the approximated density via numerical integration and inversion respectively.

This is primarily intended to be used via arithmetic on distributions: dist1 + dist2 or dist1 + dist2 + dist3 + .... Distributions with known closed-form sums (e.g. two dist_normal()) will use the exact result rather than this approximation. Chaining + automatically performs a single k-way FFT convolution rather than nested binary convolutions, avoiding compounding approximation errors.

dist_convolved(...)

Arguments

...

Two or more distribution vectors to sum. All will be recycled to a common length.

Details

Let \(Z = X_1 + X_2 + \cdots + X_k\) where the \(X_i\) are independent random variables.

Mean: $$E(Z) = \sum_{i=1}^{k} E(X_i)$$

Variance: $$\mathrm{Var}(Z) = \sum_{i=1}^{k} \mathrm{Var}(X_i)$$

Probability density function (p.d.f): $$f_Z(z) = (f_{X_1} * f_{X_2} * \cdots * f_{X_k})(z)$$

The k-way convolution is computed in a single FFT pass: each component density is evaluated on a common grid, transformed via FFT, all transforms are multiplied element-wise, and a single inverse FFT yields the result. This avoids compounding approximation errors from nested binary convolutions.

The accuracy of the FFT approximation can be controlled by passing n (number of grid points, default 2^12) and tail_p (tail probability used to find finite grid bounds for distributions with infinite support, default 1e-6) to density(), cdf(), quantile(), or generate().

See also

stats::convolve(), stats::approxfun()

Examples

# Sum of a lognormal and an exponential (no closed-form result)
d <- dist_convolved(dist_lognormal(0, 1), dist_exponential(1))
d
#> <distribution[1]>
#> [1] lN(0, 1) + Exp(1)

density(d, 2)
#> [1] 0.2696637
cdf(d, 2)
#> [1] 0.4921358
quantile(d, 0.5)
#> [1] 2.029754
generate(d, 5)
#> [[1]]
#> [1] 1.511044 9.565109 1.704500 1.546822 2.707052
#> 

# Three distributions from different families
d3 <- dist_convolved(dist_lognormal(0, 0.5), dist_gamma(2, 1), dist_exponential(2))
d3
#> <distribution[1]>
#> [1] lN(0, 0.25) + Γ(2, 1) + Exp(2)

# Via arithmetic
d2 <- dist_lognormal(0, 1) + dist_exponential(1)
density(d2, 2)
#> [1] 0.2696637

# Mean and variance are computed exactly from components
mean(d)
#> [1] 2.648721
variance(d)
#> [1] 5.670774