Bayesian update with new information Based on http://www.databozo.com/2013/09/15/Bayesian_updating_of_probability_distributions.html

bayesian_updating_distributions.R

#------------------------------------------------------------------------------#
# Bayesian update with new information
# Based on http://www.databozo.com/2013/09/15/Bayesian_updating_of_probability_distributions.html
#------------------------------------------------------------------------------#
# Load packages
library('ggplot2')

# Define functions 
normalize <- function(pos) {
  # Normalizes the possibilities so that sum = 1
  pos_sum <- sum(pos)
  pos <- pos/pos_sum
  return(pos)
}

# Update function - heads
update_success <- function(pos) {
  # Run when heads is flipped. Updates the possibilities by multiplying each 
  # current possibility by its hypothesis value, then normalizing.
  pos <- pos * hyp
  pos <- normalize(pos)
  return(pos)
}

# Update function - tails
update_failure <- function(pos) {
  # Run when tails is flipped. Updates the possibilities by multiplying each
  # current possibility by the probability the hypothesis is false (1-hyp).
  pos <- pos * (1 - hyp)
  pos <- normalize(pos)
  return(pos)
}

# Plotting function - Plot Distribution
pd <- function(pos) {
  # Function creates a qplot to display likelihood of hypothesis vs.
  # hypothesis for chance of heads.
  # y-axis limits set to .05 for consistency across plots, but may need
  # to adjust if different prior or really unfair coin used.
  qplot(x=hyp, y=pos, geom="bar", stat="identity", fill = I("blue")) +
    geom_vline(x=.50, color="red") +
    xlab("Hypothesis for chance of heads") + 
    ylab("Likelihood of hypothesis") + 
    ylim(0,0.05)
}

# Initial state
pos <- rep(1,101) # Possibilities (starting with uniform distribution)
hyp <- 0:100/100  # Hypotheses for chance of heads (0% to 100%)
pos <- normalize(pos)
pd(pos)

# Start flipping
# Heads (h=1, t=0)
pos <- update_success(pos)
pd(pos)

# Tails (h=1, t=1)
pos <- update_failure(pos)
pd(pos)

# Heads (h=2, t=1)
pos <- update_success(pos)
pd(pos)

# Head (h=3, t=1)
pos <- update_success(pos)
pd(pos)

# Tails (h=3, t=2)
pos <- update_failure(pos)
pd(pos)