We present the R package fourierin (Basulto-Elias 2017) for evaluating functions defined as Fourier-type integrals over a collection of argument values. The integrals are finitely supported with integrands involving continuous functions of one or two variables. As an important application, such Fourier integrals arise in so-called “inversion formulas”, where one seeks to evaluate a probability density at a series of points from a given characteristic function (or vice versa) through Fourier transforms. This paper intends to fill a gap in current R software, where tools for repeated evaluation of functions as Fourier integrals are not directly available. We implement two approaches for such computations with numerical integration. In particular, if the argument collection for evaluation corresponds to a regular grid, then an algorithm from Inverarity (2002) may be employed based on a fast Fourier transform, which creates significant improvements in the speed over a second approach to numerical Fourier integration (where the latter also applies to cases where the points for evaluation are not on a grid). We illustrate the package with the computation of probability densities and characteristic functions through Fourier integrals/transforms, for both univariate and bivariate examples.
Continuous Fourier transforms commonly appear in several subject areas, such as physics and statistics. In probability theory, for example, continuous Fourier transforms are related to the characteristic function of a distribution and play an important role in evaluating probability densities from characteristic functions (and vice versa) through inversion formulas (cf. Athreya and Lahiri (2006)). Similar Fourier-type integrations are also commonly required in statistical methods for density estimation, such as kernel deconvolution (cf. Meister (2009)).
At issue, the Fourier integrals of interest often cannot be solved in elementary terms and typically require numerical approximations. As a compounding issue, the oscillating nature of the integrands involved can cause numerical integration recipes to fail without careful consideration. However, Bailey and Swarztrauber (1994) present a mid-point integration rule in terms of appropriate discrete Fourier transforms, which can be efficiently computed using the Fast Fourier Transform (FFT). Inverarity (2002) extended this characterization to the multivariate integral case. These works consequently offer targeted approaches for numerically approximating types of Fourier integrals of interest (e.g., in the context of characteristic or density functions).
Because R is one of the most popular programming languages among statisticians, it seems worthwhile to have general tools available for computing such Fourier integrals in this software platform. However, we have not found any R package that specifically performs this type of integral in general, though this integration does intrinsically occur in some statistical procedures. See Stirnemann, A. Samson, and F. C. C. from Claire Lacour (2012) for an application in kernel deconvolution where univariate Fourier integrals are required. Furthermore, beyond the integral form, the capacity to handle repeated calls for such integrals is another important consideration. This need arises when computing a function, that is itself defined by a Fourier integral, over a series of points. Note that this exact aspect occurs when determining a density function from characteristic function (or vice versa), so that the ability to efficiently compute Fourier integrals over a collection of arguments is crucial.
The intent of the package fourierin explained here is to help in computing such Fourier-type integrals within R. The main function of the package serves to calculate Fourier integrals over a range of potential arguments for evaluation and is also easily adaptable to several definitions of the continuous Fourier transform and its inverse (cf. Inverarity (2002)). (That is, the definition of a continuous Fourier transform may change slightly from one context to another, often up to normalizing constants, so that it becomes reasonable to provide a function that can accommodate any given definition through scaling adjustments.) If the points for evaluating Fourier integrals are specified on regular grid, then the package allows use of the FFT for particularly fast numerical integration. However, the package also allow the user to evaluate such integrals at arbitrary collections of points that need not constitute a regular grid (though, in this case, the FFT cannot be used and computations naturally become slower). The latter can be handy in some situations; for example, evaluations at zero can provide moments of a random variable when computing derivatives of a characteristic function from the probability density. The heavy computations in fourierin are performed in C++ via the RcppArmadillo package (cf. Eddelbuettel and Sanderson (2014)).
The rest of the paper has four sections. We first describe the Fourier integral for evaluation and its numerical approximation in “Fourier Integrals and Fast Fourier Transform (FFT).” We then illustrate how package fourierin may be used in univariate and bivariate cases of Fourier integration. In “Speed Comparison,” we demonstrate the two approaches (FFT-based or not) for computing Fourier integrals, both in one and two dimensions and at several grid sizes. We provide evidence that substantial time savings occur when using the FFT-based method for evaluation at points on a grid. Finally, in “Summary,” we present conclusions and upcoming extensions to the R package.
For
This package aims to compute Fourier integrals at several points
simultaneously, which namely involves computation of the integral
At given argument
To reiterate, the goal is then to evaluate the Fourier integral
((1)) over some set of argument points
The R package fourierin can take advantage of such FFT representations for the fast computation of Fourier integrals evaluated on a regular grid. The package can also be used to evaluate Fourier integrals at arbitrary discrete sets of points. The latter becomes important when one wishes to evaluate the a continuous Fourier transform at only a few specific points (that may not necessarily constitute a regular grid). We later compare evaluation time of Fourier integrals on a regular grid, both using the FFT and without using it in fourierin.
In this section we present examples to illustrate use of the fourierin package. We begin with a univariate example which considers how to compute continuous Fourier transforms evaluated on a regular grid (therefore using the FFT operation) as well as how the computations proceed at three specified points not on a regular grid (where the FFT is not be used). The second example considers a two dimensional, or bivariate, case of Fourier integration.
The code that follows shows how the package can be used in univariate
cases. The example we consider is to recover a
For illustration, the limits of integration are set from
## -------------------------------------------------------------------
## Univariate example
## -------------------------------------------------------------------
## Load packages
library(fourierin)
library(dplyr)
library(purrr)
library(ggplot2)
## Set functions
df <- 5
cf <- function(t) (1 - 2i*t)^(-df/2)
dens <- function(x) dchisq(x, df)
## Set resolutions
resolutions <- 2^(6:8)
## Compute integral given the resoltion
recover_f <- function(resol){
## Get grid and density values
out <- fourierin(f = cf, lower_int = -10, upper_int = 10,
lower_eval = -3, upper_eval = 20,
const_adj = -1, freq_adj = -1,
resolution = resol)
## Return in dataframe format
out %>%
as_data_frame() %>%
transmute(
x = w,
values = Re(values),
resolution = resol)
}
## Density approximations
vals <- map_df(resolutions, recover_f)
## True values
true <- data_frame(x = seq(min(vals$x), max(vals$x), length = 150),
values = dens(x))
univ_plot <-
vals %>%
mutate(resolution = as.character(resolution),
resolution = gsub("64", "064", resolution)) %>%
ggplot(aes(x, values)) +
geom_line(aes(color = resolution)) +
geom_line(data = true, aes(color = "true values"))
univ_plot
## Evaluate in a nonregular grid
set.seed(666)
new_grid <- rchisq(n = 3, df = df)
resolutions <- 2^(6:9)
fourierin(f = cf, lower_int = -10, upper_int = 10,
eval_grid = new_grid,
const_adj = -1, freq_adj = -1,
resolution = 128) %>%
c () %>% Re() %>%
data_frame(x = new_grid, fx = .)
## Function that evaluates the log-density on new_grid at different
## resolutions (i.e., number of points to approximate the integral in
## the Fourier inversion formula).
approximated_fx <- function (resol) {
fourierin(f = cf, lower_int = -10, upper_int = 12,
eval_grid = new_grid,
const_adj = -1, freq_adj = -1,
resolution = resol) %>%
c() %>% Re() %>%
{data_frame(x = new_grid,
fx = dens(new_grid),
diffs = abs(. - fx),
resolution = resol)}
}
## Generate table
tab <-
map_df(resolutions, approximated_fx) %>%
arrange(x) %>%
mutate(diffs = round(diffs, 7)) %>%
rename('f(x)' = fx,
'absolute difference' = diffs)
tab
absolute difference | resolution | |||
---|---|---|---|---|
3.0585883 | 0.1541368 | 0.0002072 | 64 | |
0.0000476 | 128 | |||
0.0000472 | 256 | |||
0.0000471 | 512 | |||
6.4144242 | 0.0874281 | 0.0000780 | 64 | |
0.0000155 | 128 | |||
0.0000148 | 256 | |||
0.0000147 | 512 | |||
11.7262677 | 0.0151776 | 0.0000097 | 64 | |
0.0000025 | 128 | |||
0.0000022 | 256 | |||
0.0000022 | 512 |
Observe that the first call of the fourierin
function above has the
default argument use_fft = TRUE
. Therefore, this computation uses the
the FFT representation described in Inverarity (2002) for regular
grids, which is substantially fast (Figure 5, as
described later, provides timing comparisons without the FFT for
contrast). Also note that, when a regular evaluation grid is used,
fourierin
returns a list with both the Fourier integral values and the
evaluation grid. Figure 1 shows the resulting plot
generated. A low resolution (64) for numerical integration has been
included in order to observe differences between the true density and
its recovered version using Fourier integrals.
fourierin()
function for univariate function at
resolution At the bottom of the code above, we also show how fourierin()
works
when a non-regular “evaluation grid” is provided. Observe that, in this
case, one directly specifies separate points for evaluation of the
integral in addition to separately specifying a resolution level for
integration. This aspect is unlike the evaluation case on a regular
grid. Consequently, only the Fourier integral values are returned, which
is also unlike the regular grid case (where the evaluation grid is
returned with corresponding integrals in a list). Note that the function
In a second example, to illustrate how the fourierin()
function works
for bivariate functions, we use a bivariate normal density
Below is the code for this bivariate case using a regular evaluation
grid, where the output is a complex matrix whose components are Fourier
integrals corresponding to the gridded set of bivariate arguments. As
illustration, the limits of integration are set from
## -------------------------------------------------------------------
## Bivariate example
## -------------------------------------------------------------------
## Load packages
library(fourierin)
library(tidyr)
library(dplyr)
library(purrr)
library(lattice)
library(ggplot2)
## Set functions to be tested with their corresponding parameters.
mu <- c(-1, 1)
sig <- matrix(c(3, -1, -1, 2), 2, 2)
## Multivariate normal density, x is n x d
f <- function(x) {
## Auxiliar values
d <- ncol(x)
z <- sweep(x, 2, mu, "-")
## Get numerator and denominator of normal density
num <- exp(-0.5*rowSums(z * (z %*% solve(sig))))
denom <- sqrt((2*pi)^d*det(sig))
return(num/denom)
}
## Characteristic function, s is n x d
phi <- function (s) {
complex(modulus = exp(-0.5*rowSums(s*(s %*% sig))),
argument = s %*% mu)
}
## Evaluate characteristic function for a given resolution.
eval <- fourierin(f,
lower_int = c(-8, -6), upper_int = c(6, 8),
lower_eval = c(-4, -4), upper_eval = c(4, 4),
const_adj = 1, freq_adj = 1,
resolution = 2*c(64, 64),
use_fft = T)
## Evaluate true and approximated values of Fourier integral
dat <- eval %>%
with(crossing(y = w2, x = w1) %>%
mutate(approximated = c(values))) %>%
mutate(true = phi(matrix(c(x, y), ncol = 2)),
difference = approximated - true) %>%
gather(value, z, -x, -y) %>%
mutate(real = Re(z), imaginary = Im(z)) %>%
select(-z) %>%
gather(part, z, -x, -y, -value)
## Surface plot
wireframe(z ~ x*y | value*part, data = dat,
scales =
list(arrows=FALSE, cex= 0.45,
col = "black", font = 3, tck = 1),
screen = list(z = 90, x = -74),
colorkey = FALSE,
shade=TRUE,
light.source= c(0,10,10),
shade.colors = function(irr, ref,
height, w = 0.4)
grey(w*irr + (1 - w)*(1 - (1 - ref)^0.4)),
aspect = c(1, 0.65))
## Contours of values
biv_example1 <-
dat %>%
filter(value != "difference") %>%
ggplot(aes(x, y, z = z)) +
geom_tile(aes(fill = z)) +
facet_grid(part ~ value) +
scale_fill_distiller(palette = "Reds")
biv_example1
## Contour of differences
biv_example2 <-
dat %>%
filter(value == "difference") %>%
ggplot(aes(x, y, z = z)) +
geom_tile(aes(fill = z)) +
facet_grid(part ~ value) +
scale_fill_distiller(palette = "Spectral")
biv_example2
The result of fourierin()
was stored above in eval
, which is a list
with three elements: two vectors with the corresponding evaluation grid
values in each coordinate direction and a complex matrix containing the
Fourier integral values. If we do not wish to evaluate the Fourier
integrals on a regular grid and instead wish to evaluate these at, say
w
and the function will return a vector of size
fourierin
function for univariate function at
resolution Corresponding to this bivariate example, we have generated three plots to compare the approximation from Fourier integrals to the underlying truth (i.e., compare the approximated and true characteristic functions of the bivariate normal distribution). In Figure 2, we present the surface plots of the approximated and the true values, as well as their differences for both the real and imaginary parts. One observes that differences are small, indicating the adequacy of the numerical integration.
fourierin
function for univariate function at
resolution For a different perspective of the resulting Fourier integration, Figure 3 presents a contour plot showing the approximated and true values of the bivariate normal characteristic function, for both real and imaginary parts. We show a tile plot of the differences in Figure 4. Observe that the range of differences in Figure 4 is relatively much smaller than the values in Figure 3.
fourierin
function for univariate function at
resolution Through a small numerical study, here we compare the differences in
execution times using fourierin()
for integration at points on a
regular grid, both with or without FFT steps, considering univariate and
bivariate Fourier integrals. Figure 5 shows timing
results for a univariate example of the integral in ((1))
evaluated on a grid, while Figure 6 presents timing
results for a bivariate example. Note that the reported time units
differ between these figures, as the bivariate case naturally requires
more time. These figures provide evidence that, for evaluating integrals
on a regular grid, the FFT option in fourierin()
creates large
advantages in time.
The code that was used to generate Figure 5 and 6 is below.
## -------------------------------------------------------------------
## Univariate speed test
## -------------------------------------------------------------------
library(fourierin)
library(dplyr)
library(purrr)
library(ggplot2)
library(microbenchmark)
## Test speed at several resolutions
resolution <- 2^(3:8)
## Function to be tested
myfnc <- function(t) exp(-t^2/2)
## Aux. function
compute_times <- function(resol){
out <-
microbenchmark(
fourierin_1d(f = myfnc, -5, 5, -3, 3, -1, -1, resol),
fourierin_1d(f = myfnc, -5, 5, -3, 3, -1, -1, resol,
use_fft = FALSE),
times = 5) %>%
as.data.frame()
## Rename levels
levels(out$expr) <- c("yes", "no")
## Obtain median of time.
out %>%
group_by(expr) %>%
summarize(time = median(time*1e-6),
resolution = resol) %>%
rename(FFT = expr)
}
speed1 <- resolution %>%
map_df(compute_times) %>%
mutate(resolution = as.factor(resolution)) %>%
ggplot(aes(resolution, log(time), color = FFT)) +
geom_point(size = 2, aes(shape = FFT)) +
geom_line(aes(linetype = FFT, group = FFT)) +
ylab("time (in log-milliseconds)")
speed1
## -------------------------------------------------------------------
## Bivariate test
## -------------------------------------------------------------------
## Load packages
library(fourierin)
library(dplyr)
library(purrr)
library(ggplot2)
library(microbenchmark)
## Test speed at several resolutions
resolution <- 2^(3:7)
## Bivariate function to be tested
myfnc <- function(x) dnorm(x[, 1])*dnorm(x[, 2])
## Aux. function
compute_times <- function(resol){
resol <- rep(resol, 2)
out <-
microbenchmark(
fourierin(myfnc,
lower_int = c(-8, -6), upper_int = c(6, 8),
lower_eval = c(-4, -4), upper_eval = c(4, 4),
const_adj = 1, freq_adj = 1,
resolution = resol),
fourierin(myfnc,
lower_int = c(-8, -6), upper_int = c(6, 8),
lower_eval = c(-4, -4), upper_eval = c(4, 4),
const_adj = 1, freq_adj = 1,
resolution = resol, use_fft = FALSE),
times = 3) %>%
as.data.frame()
## Rename levels
levels(out$expr) <- c("yes", "no")
## Obtain median of time.
out %>%
group_by(expr) %>%
summarize(time = median(time*1e-9),
resolution = resol[1]) %>%
rename(FFT = expr)
}
## Values
comparison <-
resolution %>%
map_df(compute_times)
fctr_order <-
unique(comparison$resolution) %>%
paste(., ., sep = "x")
## Plot
speed2 <- comparison %>%
mutate(resolution = paste(resolution, resolution, sep = "x"),
resolution = ordered(resolution, levels = fctr_order)) %>%
ggplot(aes(resolution, log(time), color = FFT)) +
geom_point(size = 2, aes(shape = FFT)) +
geom_line(aes(linetype = FFT, group = FFT)) +
ylab("time (in log-seconds)")
speed2
Continuous Fourier integrals/transforms are useful in statistics for computation of probability densities from characteristic functions, as well as the reverse, when describing probability structure; see the “Examples” section for some demonstrations. The usefulness and potential application of Fourier integrals, however, also extends to other contexts of physics and mathematics, as well as to statistical inference (e.g., types of density estimation). For this reason, we have developed the fourierin package as a tool for computing Fourier integrals over collections of evaluation points, where repeat evaluation steps and often complicated numerical integrations are involved. When evaluation points fall on a regular grid, fourierin allows use of a Fast Fourier Transform as a key ingredient for rapid numerical approximation of Fourier-type integrals.
In “Speed Comparison,” we presented evidence of the gain in time when
using this fast implementation of fourierin()
on regular grids, while
we also illustrated the versatility of fourierin in “Examples”
section. At present (version 0.2.1), the fourierin package performs
univariate and bivariate Fourier integration. An extension of the
package to address higher dimensional integration will be included in
future versions.
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Basulto-Elias, et al., "fourierin: An R package to compute Fourier integrals", The R Journal, 2017
BibTeX citation
@article{RJ-2017-044, author = {Basulto-Elias, Guillermo and Carriquiry, Alicia and Brabanter, Kris De and Nordman, Daniel J.}, title = {fourierin: An R package to compute Fourier integrals}, journal = {The R Journal}, year = {2017}, note = {https://rjournal.github.io/}, volume = {9}, issue = {2}, issn = {2073-4859}, pages = {72-83} }