log_convolve: Logarithm of the convolution of 'exp(log_pmf1)' and...

In huonw/sisfft: Fast and accurate convolutions and tail probabilities

Description

log_convolve returns the logarithm of an approximation of the convolution of pmf1 = exp(log_pmf1) with pmf2 = exp(log_pmf2). The approximation has relative error bounded by beta, after accounting for the lower bound delta. Denote p the exact convolution of pmf1 and pmf2, q the exact convolution of those vectors after setting any value less than delta to zero, and finally v the approximation. The approximation satisfies (1 - beta) * q[i] <= v[i] <= (1 + beta) * p[i] for all i.

Usage

log_convolve(log_pmf1, log_pmf2, beta, delta = 0)

Arguments

log_pmf1	a vector that is the logarithm of a pmf.
log_pmf2	a vector that is the logarithm of a pmf.
beta	a positive number less than 0.5, to control the relative error. Note that this differs to Wilson, Keich 2016 which describes aFFT-C with alpha := 1 / beta.
delta	a non-negative number, to act as a lower bound on elements needed.

Details

log_convolve is a log-space implementation of trim-aFFT-C, which notably falls back to aFFT-C when delta = 0. This implementation uses the more natural beta = 1/alpha parameter in place of the alpha used in the published algorithm.

Value

A vector that is the logarithm of the approximation to the convolution.

References

Huon Wilson, Uri Keich, Accurate pairwise convolutions of non-negative vectors via FFT, Computational Statistics & Data Analysis, Volume 101, September 2016, Pages 300-315, ISSN 0167-9473, http://dx.doi.org/10.1016/j.csda.2016.03.010.

Huon Wilson, Uri Keich, Accurate small tail probabilities of sums of iid lattice-valued random variables via FFT, Submitted.

See Also

log_convolve_power log_pvalue

Examples

pmf1 <- c(0.5, 0.5)
pmf2 <- c(0.1, 0.2, 0.3, 0.4)

log_conv <- log_convolve(log(pmf1), log(pmf2), 1e-3, 0)
conv <- exp(log_conv)

exact <- c(0.05, 0.15, 0.25, 0.35, 0.20)
(exact * (1 - 1e-3) < conv) & (conv < exact * (1 + 1e-3))

# the same convolution with a bound that allows dropping the 0.1
log_conv_bound <- log_convolve(log(pmf1), log(pmf2), 1e-3, 0.15)
conv_bound <- exp(log_conv_bound)

exact_bound <- c(0.00, 0.10, 0.25, 0.35, 0.20)
(exact_bound * (1 - 1e-3) <= conv_bound) & (conv_bound <= exact * (1 + 1e-3))

# FFT-based convolution can have large relative error:
# relative error that behaves nicely with zeros:
rel <- function(x, true) {
    ifelse(true == 0, ifelse(x == 0, 0, Inf), abs((x - true) / true))
}

# example pmf from Wilson, Keich (2016).
b <- 1e-5
c <- 1e-20
a <- 1 - b - c
p <- c(0, a, b, c)

exact <- c(0, 0, a^2, 2*a*b, 2*a*c + b^2, 2*b*c, c^2)
fftc <- convolve(p, rev(p), type = "o")
afftc <- exp(log_convolve(log(p), log(p), 1e-3, 0))

rel(fftc, exact)
rel(afftc, exact)

huonw/sisfft documentation built on May 17, 2019, 9:13 p.m.