## How to solve the Fourier Series problem

install.packages("fourierin", dependencies = FALSE)
> library("dplyr", lib.loc="~/R/win-library/3.6")
> library("ggplot2", lib.loc="~/R/win-library/3.6")
install.packages("purrrlyr", dependencies = FALSE)
install.packages("purrr", dependencies = FALSE)
library("dplyr", lib.loc="~/R/win-library/3.6")
> library("lattice", lib.loc="C:/Program Files/R/R-3.6.1/library")
> 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,
+                      resolution = resol)
+     ## Return in dataframe format
+     out %>%
+         as_data_frame() %>%
+         transmute(
+             x = w,
+             values = Re(values),
+             resolution = resol)
+ }
> vals <- map_df(resolutions, recover_f)
>
> true <- data_frame(x = seq(min(vals\$x), max(vals\$x), length = 150),
+                    values = dens(x))
Warning message:
`data_frame()` is deprecated, use `tibble()`.
This warning is displayed once per session.
>
> 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
> 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,
+           resolution = 128) %>%
+     c () %>% Re() %>%
+     data_frame(x = new_grid, fx = .)
# A tibble: 3 x 2
x     fx
<dbl>  <dbl>
1  6.41 0.0874
2 11.7  0.0152
3  3.06 0.154
> approximated_fx <- function (resol) {
+     fourierin(f = cf, lower_int = -10, upper_int = 12,
+               eval_grid = new_grid,
+               resolution = resol) %>%
+         c() %>% Re() %>%
+         {data_frame(x = new_grid,
+                     fx = dens(new_grid),
+                     diffs = abs(. - fx),
+                     resolution = resol)}
+ }
>
> tab <-
+     map_df(resolutions, approximated_fx) %>%
+     arrange(x) %>%
+     mutate(diffs = round(diffs, 7)) %>%
+     rename('f(x)' = fx,
+            'absolute difference' = diffs)
> tab
# A tibble: 12 x 4
x `f(x)` `absolute differe~ resolution
<dbl>  <dbl>              <dbl>      <dbl>
1  3.06 0.154           0.000207          64
2  3.06 0.154           0.0000476        128
3  3.06 0.154           0.0000472        256
4  3.06 0.154           0.0000471        512
5  6.41 0.0874          0.000078          64
6  6.41 0.0874          0.0000155        128
7  6.41 0.0874          0.0000148        256
8  6.41 0.0874          0.0000147        512
9 11.7  0.0152          0.0000097         64
10 11.7  0.0152          0.0000025        128
11 11.7  0.0152          0.0000022        256
12 11.7  0.0152          0.0000022        512
> library("dplyr", lib.loc="~/R/win-library/3.6")
> library("lattice", lib.loc="C:/Program Files/R/R-3.6.1/library")
mu <- c(-1, 1)
> sig <- matrix(c(3, -1, -1, 2), 2, 2)
> 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),
+                   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)
Error in crossing(y = w2, x = w1) : could not find function "crossing"
> 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,
+           light.source= c(0,10,10),
+                                   height, w = 0.4)
+               grey(w*irr + (1 - w)*(1 - (1 - ref)^0.4)),
+           aspect = c(1, 0.65))
Error in wireframe.formula(z ~ x * y | value * part, data = dat, scales = list(arrows = FALSE,  :
>
> mu <- c(-1, 1)
> sig <- matrix(c(3, -1, -1, 2), 2, 2)
> 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)
+ }
> phi <- function(s) {
+     complex(modulus = exp(- 0.5*rowSums(s*(s %*% sig))),
+             argument = s %*% mu)
+ }
>
> eval <- fourierin_2d(f, lower_int = c(-8, -6), upper_int = c(6, 8),
+                      lower_eval = c(-4, -4), upper_eval = c(4, 4),
+                      resolution = c(128, 128))
>
> 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
>
> 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)
> 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)
> > myfnc <- function(t) exp(-t^2/2)
> out <- fourierin(f = myfnc, lower_int = -5, upper_int = 5,
+                  lower_eval= -3, upper_eval = 3, const_adj = -1,
+                  freq_adj = -1, resolution = 64)
> grid <- out\$w
> values <- Re(out\$values)
> plot(grid, values, type = "l", col = 3)
> lines(grid, dnorm(grid), col = 4)
> myfnc <- function(t) dgamma(t, shape, rate)
> shape <- 5
> rate <- 3
> out <- fourierin(f = myfnc, lower_int = 0, upper_int = 6,
+                  lower_eval = -4, upper_eval = 4,
> grid <- out\$w
> re_values <- Re(out\$values)
> im_values <- Im(out\$values)
> true_cf <- function(t, shape, rate) (1 - 1i*t/rate)^-shape
> true_re <- Re(true_cf(grid, shape, rate))
> true_im <- Im(true_cf(grid, shape, rate))
> plot(grid, re_values, type = "l", col = 3)
> lines(grid, true_re, col = 4)
> plot(grid, im_values, type = "l", col = 3)
> lines(grid, true_im, col = 4)
> mu <- c(-1, 1)
> sig <- matrix(c(3, -1, -1, 2), 2, 2)
> 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)
+ }
> phi <- function(s) {
+     complex(modulus = exp(- 0.5*rowSums(s*(s %*% sig))),
+             argument = s %*% mu)
+ }
> 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 = c(128, 128))
> t1 <- eval\$w1
> t2 <- eval\$w2
> t <- as.matrix(expand.grid(t1 = t1, t2 = t2))
> approx <- eval\$values
> true <- matrix(phi(t), 128, 128)
> 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)
> 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)