fourierin_2d {fourierin}

Bivariate Fourier integrals

Description

It computes Fourier integrals for functions of one and two variables.

Usage

fourierin_2d(
  f,
  lower_int,
  upper_int,
  lower_eval = NULL,
  upper_eval = NULL,
  const_adj,
  freq_adj,
  resolution = NULL,
  eval_grid = NULL,
  use_fft = TRUE
)

Arguments

f

function or a vector of size m. If a function is provided, it must be able to be evaluated at vectors. If a vector of values is provided, such evaluations must have been obtained on a regular grid and the Fourier integral is faster is m is a power of 2.

lower_int

Lower integration limit(s).

upper_int

Upper integration limit(s).

lower_eval

Lower evaluation limit(s). It can be NULL if an evaluation grid is provided.

upper_eval

Upper evaluation limit(s). It can be NULL if an evaluation grid is provided.

const_adj

Factor related to adjust definition of Fourier transform. It is usually equal to 0, -1 or 1.

freq_adj

Constant to adjust the exponent on the definition of the Fourier transform. It is usually equal to 1, -1, 2pi or -2pi.

resolution

A vector of integers (faster if powers of two) determining the resolution of the evaluation grid. Not required if f is a vector.

eval_grid

Optional matrix with d columns with the points where the Fourier integral will be evaluated. If it is provided, the FFT will not be used.

use_fft

Logical value specifying whether the FFT will be used.

Value

If w is given, only the values of the Fourier integral are returned, otherwise, a list with three elements

w1

Evaluation grid for first entry

w2

Evaluation grid for second entry

values

m1 x m2 matrix of complex numbers, corresponding to the evaluations of the integral

Examples

##--- Recovering std. normal from its characteristic function -----
library(fourierin)

##-Parameters of bivariate normal distribution
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)
}

##-Approximate cf using Fourier integrals
eval <- fourierin_2d(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 = c(128, 128))

## Extract values
t1 <- eval$w1
t2 <- eval$w2
t <- as.matrix(expand.grid(t1 = t1, t2 = t2))
approx <- eval$values
true <- matrix(phi(t), 128, 128)        # Compute true values

##-This is a section of the characteristic functions
i <- 65
plot(t2, Re(approx[i, ]), type = "l", col = 2,
     ylab = "",
     xlab = expression(t[2]),
     main = expression(paste("Real part section at ",
                             t[1], "= 0")))
lines(t2, Re(true[i, ]), col = 3)
legend("topleft", legend = c("true", "approximation"),
       col = 3:2, lwd = 1)

##-Another section, now of the imaginary part
plot(t1, Im(approx[, i]), type = "l", col = 2,
     ylab = "",
     xlab = expression(t[1]),
     main = expression(paste("Imaginary part section at ",
                             t[2], "= 0")))
lines(t1, Im(true[, i]), col = 3)
legend("topleft", legend = c("true", "approximation"),
       col = 3:2, lwd = 1)