bayes_anim.R

library(tidyverse)
library(gganimate)
library(magick)

# Number of Bern trials to run
n_trials <- 100

# Grid for plotting posterior
x <- seq(0,1, length.out = 500)

# Simulate cumulative number of successes in 100 Bern(0.7) trials
success <- map_dbl(1:n_trials, ~rbinom(1,1,0.7)) %>% cumsum()

df <-
  tibble(trial=1:n_trials, success=success, failure=trial-success) %>%
  mutate(p1 = success+1) %>%
  mutate(p2 = failure+1)

# Create a data frame where each row is a point on the PDF of one of the
# posterior distributions. So 500 grid points x 100 posteriors = 50,000 rows
df <-
  map2_dfr(df$p1, df$p2, ~tibble(x =x, density=dbeta(x, .x, .y), p1=.x, p2=.y), .id = "trial") %>%
  mutate(trial = as.numeric(trial))

# Compute 95% credible interval endpoints
df$lower <- map2_dbl(df$p1, df$p2, ~qbeta(0.025, shape1=.x, shape2=.y))
df$upper <- map2_dbl(df$p1, df$p2, ~qbeta(0.975, shape1=.x, shape2=.y))

# Create posterior PDF animation
p <- 
  df %>% 
  ggplot(aes(x=x, y=density)) +
  geom_line() +
  geom_line(data = tibble(x=c(0.7, 0.7), y=c(0, max(df$density))),
            mapping = aes(x=x, y=y),
            linetype=2) +
  geom_vline(aes(xintercept=lower), linetype=3, color="gray", alpha=0.7) +
  geom_vline(aes(xintercept=upper), linetype=3, color="gray", alpha=0.7) +
  labs(title = "Binomial Proportion Inference",
       subtitle="(Uniform prior) True p = 0.7",
       x = "p",
       y= "Probability Density",
       caption = "Posterior distribution with 95% Credible Interval") +
  cowplot::theme_cowplot(font_family = "Source Sans Pro",
                         font_size = 10) +
  theme(legend.position = "none") +
  transition_states(trial, state_length = 1, transition_length = 2, wrap=FALSE) +
  ease_aes('quartic-in-out')

# Specify animation parameters
p_gif <- animate(p,
                 nframes=600,
                 fps=50,
                 end_pause = 50,
                 height = 3,
                 width = 2.7,
                 units = "in",
                 res = 300)

# Put frames into magick format
p_mgif <- image_read(p_gif)

##################
# Create animation for bar chart

df <-
  tibble(trial=1:n_trials, Success=success, Failure=trial-success) %>%
  pivot_longer(cols=Success:Failure)

b <-
df %>%
  ggplot(aes(x=name, y=value, fill=name)) +
  geom_bar(stat = "identity", position = 'dodge') +
  labs(title = "Observations: {closest_state}",
       x="",
       y= "Frequency",
       caption = " ") +
  cowplot::theme_cowplot(font_family = "Source Sans Pro",
                         font_size = 10) +
  theme(legend.position = "none")+
  transition_states(trial, state_length = 1, transition_length = 2, wrap=FALSE) +
  ease_aes('quartic-in-out')

# Specify animation parameters
b_gif <- animate(b,
                 nframes=600,
                 fps=50,
                 end_pause = 50,
                 height = 3,
                 width = 2,
                 units = "in",
                 res = 300)

# Put frames into magick format
b_mgif <- image_read(b_gif)

##################################

# Put plots side by side

# Merge the top two panels
# Create the first frame
new_gif <- image_append(c(p_mgif[1], b_mgif[1]))

# Append the other frames
for(i in 2:600){
  combined <- image_append(c(p_mgif[i], b_mgif[i]))
  new_gif <- c(new_gif, combined)
}

# Save the combined animation to a GIF
anim_save(animation = new_gif, filename = "out.gif")

# Convert GIF to MP4 with ffmpeg using the command line
system('ffmpeg -i out.gif -movflags faststart -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" beta_bayes.mp4')