sisfft
An R implementation of sisFFT and aFFT-C.
See also: a Python implementation https://github.com/huonw/sisfft-py
Example of using sisFFT to compute a pvalue, along with doing a raw convolution power:
require(sisfft)
log_pmf <- log(c(0.1, 0.3, 0.2, 0.4))
L <- 10 # 10-fold convolution
s0 <- 26 # threshold
error_limit <- 1e-3
# compute the pvalue of s0 directly
print(exp(sisfft::log_pvalue(log_pmf, s0, L, error_limit)))
# => [1] 0.04408607
# now compute the full 10-fold convolution (delta limit of zero means no truncation)
log_full_conv <- sisfft::log_convolve_power(log_pmf, L, 1.0/error_limit, 0.0)
# use this to compute the pvalue
pvalue <- sum(exp(log_full_conv)[s0:length(log_full_conv)])
print(pvalue)
# => [1] 0.04408607Example of aFFT-C:
require(sisfft)
pmf1 <- c(0.1, 0.3, 0.2, 0.4)
pmf2 <- c(0.25, 0.25, 0.25, 0.25)
log_pmf1 <- log(pmf1)
log_pmf2 <- log(pmf2)
alpha <- 1e3 # accuracy parameter
# Compute pmf1 * pmf2 in two ways:
# first, with R's built in convolve (note that this uses FFT, but this
# particular convolution is guaranteed to be accurate)
direct <- convolve(pmf1, rev(pmf2), type = "o")
# and the with aFFT-C
via_afftc <- sisfft::log_convolve(log_pmf1, log_pmf2, alpha)
print(direct)
# => [1] 0.025 0.100 0.150 0.250 0.225 0.150 0.100
print(exp(via_afftc))
# => [1] 0.025 0.100 0.150 0.250 0.225 0.150 0.100
# They're not exactly the same, but aFFT-C is well within the error bounds:
rel_error <- (exp(via_afftc) - direct) / direct
print(rel_error)
# => [1] 1.387779e-15 -1.387779e-16 0.000000e+00 -4.440892e-16 -1.233581e-16
# [6] 3.700743e-16 1.387779e-16
