# ----------------------------------------------------
# Create animation for visualising likelihood function
# ----------------------------------------------------

# Create a working directory for this file
# Add an images folder
# Run script
# Create GIF using Photoshop

# Set binomial distribution parameters
n       = 10  # Sample size
k       = 0:n # Discrete probability space
k.heads = 2   # Number of heads thrown

# Set theta range and interval
thetas = seq(0, 1, .01)
thetas = c(thetas, rev(thetas[c(-1,-length(thetas))]))

# Set bar colors
col = rep("beige", n+1) # all bars beige
col[k.heads+1] = 'cyan' # Data to cyan

# First manually create images folder. Write files to images folder.
png(file="images/binomial_likelihood_function%02d.png", width=800, height=300)

for(theta in thetas) {

  # Set layout for two plots
  layout(matrix(1:2, 1, 2))

  # Plot binomial distribution
  title = bquote('Binomial distribution for' ~ theta ~ '=' ~ .(theta))

  barplot( dbinom(k, n, theta),
           main      = title,
           names.arg = 0:10,
           xlab      = "number of head",
           ylab      = "P(head)",
           beside    = TRUE,
           col       = col,
           las       = 1,
           ylim      = c(0,1)
         )

  # Plot Likelihood function
  plot( thetas, dbinom(k.heads, n,thetas),
        type = 'l',
        lwd  = 3,
        xlim = c(0,1),
        main = bquote("Likelihood function for" ~ .(k.heads)~ "/" ~.(n) ~ "heads"),
        ylab = expression(paste('Likelihood - ', "P(D|", theta, ")" )),
        xlab = expression(theta),
        ylim = c(0,1),
        las  = 1
      )

  # Plot likelihood dot on likelihood function
  points(theta, dbinom(k.heads, n, theta),
         cex = 2,
         pch = 21,
         bg  = "cyan")

  # Add tracking line
  lines(c(theta,theta), c(0,dbinom(k.heads, n, theta)), lty='dashed')

}

dev.off()