Animation of Gaussian and Poissonian functions and their confidence intervals

gauss_v_poisson_CI.R

library(tidyverse)
library(ggplot2)
library(survival)

lambdas <- seq(300, 2)
pic_dir <- "pics"
fps <- 24
output_video <- "output.mp4"

# construct dataframes
df <- tibble()
lines <- tibble()
errors <- tibble()
for(lambda in lambdas){
  range <- seq(0, (2*lambda)+3) # +3 to show at small lambda values
  
  df_to_append <- tibble(x = c(range, range),
                         y = c(dpois(x = range, lambda = lambda), dnorm(x = range, mean = lambda, sd = sqrt(lambda))),
                         type = c(rep("poisson", length(range)), rep("gauss", length(range))),
                         lambda = rep(lambda, 2*length(range)))
  df <- df %>% bind_rows(df_to_append)
  
  ci_95_poisson <- cipoisson(k = lambda, p = .95, method="exact")
  ci_95_gauss_u <- qnorm(.975, sd = sqrt(lambda), mean = lambda)  # 95% CI on RHS
  ci_95_gauss_l <- qnorm(.025, sd = sqrt(lambda), mean = lambda)  # 95% CI on LHS
  
  
  lines_to_append <- tibble(x = c(ci_95_poisson["upper"], ci_95_poisson["lower"], ci_95_gauss_u, ci_95_gauss_l),
                            type = c("poisson", "poisson", "gauss", "gauss"),
                            lambda = rep(lambda, 4),
                            u_or_l = rep(c("u", "l"), 2))
  lines <- lines %>% bind_rows(lines_to_append)
  
  upper_rel_error <- 100*abs((ci_95_gauss_u - ci_95_poisson["upper"]) / ci_95_gauss_u)
  lower_rel_error <- 100*abs((ci_95_gauss_l - ci_95_poisson["lower"]) / ci_95_gauss_l)
  
  errors_to_append <- tibble(upper = upper_rel_error, lower = lower_rel_error, lambda = lambda)
  errors <- errors %>% bind_rows(errors_to_append)
}

for(i in 1:length(lambdas)){
  lambda_sample <- lambdas[i]
  
  filename <- file.path(pic_dir, paste0("pic_", i, ".png"))
  
  df_to_plot <- df %>% filter(lambda == lambda_sample)
  
  lines_to_plot <- lines %>% 
    filter(lambda == lambda_sample)
  
  poisson_lines <- lines_to_plot %>% filter(type == "poisson")
  
  gauss_lines <- lines_to_plot %>% filter(type == "gauss")
  
  upper_error <- str_pad(signif(filter(errors, lambda == lambda_sample)$upper, 3), width = 5, side = "left")
  lower_error <- str_pad(signif(filter(errors, lambda == lambda_sample)$lower, 3), width = 5, side = "left")
  
  caption_text <- paste0("upper CI difference: ", upper_error, "%    \nlower CI difference: ", lower_error, "%    ")
  
  p <- df_to_plot %>% 
    ggplot(aes(x = x, y = y, color = type)) +
    geom_line(size = 2) +
    ylim(0, NA) +
    labs(title = "Gaussian & Poissonian 95% CI comparison", x = "", y = "", caption = caption_text) + 
    geom_point(data = filter(df_to_plot, type == "poisson"), size = 4, show.legend = FALSE) +
    geom_vline(data = lines_to_plot, aes(xintercept = x, color = type), show.legend = FALSE) +
    geom_vline(aes(xintercept = lambda_sample), color = "grey", show.legend = FALSE) +
    annotate("rect", xmin = filter(gauss_lines, u_or_l == "l")$x, 
             xmax = filter(poisson_lines, u_or_l == "l")$x, 
             ymin = 0, ymax = max(df_to_plot$y),
             fill = "black", alpha = .1) +
    annotate("rect", xmin = filter(gauss_lines, u_or_l == "u")$x, 
             xmax = filter(poisson_lines, u_or_l == "u")$x, 
             ymin = 0, ymax = max(df_to_plot$y),
             fill = "black", alpha = .1) +
    xlim(0, 2*lambda_sample+3) +
    scale_x_continuous(breaks = c(lambda_sample), labels = c(lambda_sample)) +
    theme_minimal() +
    theme(legend.title = element_blank(),
          legend.position = "top",
          legend.justification = "left",
          legend.text = element_text(size = 20, family = "sans", face = "bold", margin = margin(0, 40, 0, 0)),
          plot.title = element_text(family = "sans", face = "bold", size = 32, margin = margin(20, 20, 40, 20)),
          plot.caption = element_text(family = "mono", size = 20, margin = margin(0, 0, 15, 0)),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          axis.line.x = element_line(colour = "black"),
          axis.text.x = element_text(size = 20),
          axis.text.y = element_blank())
  
  scale <- 1
  ggsave(plot = p, filename = filename, width = scale*19.2, height = scale*10.8, dpi = 100)
}

#### FFMPEG ####
command <- paste0("ffmpeg -y -r ", fps, " -f image2 -s 1920x1080 -i ", pic_dir,
                  "/pic_%d.png -vcodec libx264 -crf 20 -pix_fmt yuv420p ", output_video)
system(command = command)

gauss_v_poisson_CI.R

Creates an animation showing similar Gaussian and Poissonian functions, along with their 95% confidence intervals. Animation can be found here.

---