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.04408607
Example 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.