#
# Chapter 4 -- Bayesian Computation With R
# Demonstrates Dirichlet on Electoral Votes
#
# Remove all objects just to be safe
#
rm(list=ls(all=TRUE))
#
library(LearnBayes)
#
data(election.2008)
attach(election.2008)
#
# Assume 500 voters for each of the 51 (D.C. is included) "states"
# Let Q_Oj and Q_Mj be the proportion of voters for Obama and McCain
# in the jth State
# Assume that the Prior distribution of 51 proprotions for Obama and McCain
# are uniformly distributed
# Then the Posterior distribution of voters is a Dirichelt distribution with
# parameters:
# (500*Q_Oj+1, 500*Q_Mj+1, 500*(1-Q_Oj - Q_Mj)+1)
#
# Function to Compute Obama Probability
#
prob.Obama <- function(j)
{
p <- rdirichlet(5000, 500*c(M.pct[j],O.pct[j],100-M.pct[j]-O.pct[j])/100+1)
mean(p[,2] > p[,1])
}
#
# Function To Simulate the "coin flips" for Obama wins in the 51 states -- then
# sum the corresponding electoral vote
# Note -- rbinom(n, size, prob)
#
# Has the arguments:
#
# n number of observations. If length(n) > 1, the length is taken to be the number required.
# size number of trials (zero or more).
# prob probability of success on each trial.
#
sim.election <- function()
{
winner <- rbinom(51,1,Obama.win.probs)
sum(EV*winner) # Note that EV is in election.2008
}
#
# Compute the Obama Win Probability for all states by using the sapply
# function in R
# &&&&&&&&&&&&& NOTE -- THIS IS DOCUMENTED UNDER LAPPLY !!! &&&&&&&&&&&&&&&&&&&&&&&
# sapply is a user-friendly version of lapply by default returning a vector or matrix if appropriate.
# Usage
#
# sapply(X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)
#
# Arguments
# X a vector (atomic or list) or an expressions vector. Other objects (including classed objects) will be coerced by as.list.
# FUN the function to be applied to each element of X: see ‘Details’. In the case of functions like +, %*%, etc., the function name must be backquoted or quoted.
# ... optional arguments to FUN.
# simplify logical; should the result be simplified to a vector or matrix if possible?
# USE.NAMES logical; if TRUE and if X is character, use X as names for the result unless it had names already.
# n number of replications.
# expr expression (language object, usually a call) to evaluate repeatedly.\
#
Obama.win.probs <- sapply(1:51,prob.Obama)
#
# Use the Obama win probabilities by flipping the 51 "coins" 1000 times
# to simulate the output of the election -- This is a lot like the
# parametric bootstrap
#
sim.EV <- replicate(1000,sim.election())
hist(sim.EV,min(sim.EV):max(sim.EV),xlab="",ylab="",main="",col="red",font=2)
# Main title
mtext("Histogram of 1000 Simulated Draws of the Total Electoral Vote\nFor Senator Obama in the 2008 Election",side=3,line=1.00,cex=1.2,font=2)
# x-axis title
mtext("Simulated Electoral Votes",side=1,line=2.75,font=2,cex=1.2)
# y-axis title
mtext("Frequency",side=2,line=2.75,font=2,cex=1.2)
#
abline(v=365,lwd=3) # Obama Received 365 electoral votes
text(375,30,"Actual \n Obama \n Total")
#
# this detaches the data -- pos=3 is were it is stored
#
detach(election.2008, pos=2)
#detach(election.2008, pos=3)
#detach(election.2008, pos=4)