R code to check the posterior prob of a rigged coin-tossing experiment

cointoss.R

# calculate the prob of P(fake|data)
calc_prob <- function(heads)
{
    P_fake = .1 # my prior belief of the experiment being rigged
    P_real = 1 - P_fake

    P_data_given_fake = dnorm(heads, mean=10000, sd=10) # P(data|fake)
    P_data_given_real = dbinom(heads, size=20000, p=.5) # P(data|real)

    # bayes theorem
    P_data = P_data_given_fake * P_fake + P_data_given_real * P_real
    P_fake_given_data = P_data_given_fake * P_fake / P_data

    P_fake_given_data
}

calc_prob(10000)
calc_prob(10093)

library(ggplot2)

x <- 9900:10100
fake_dist <- dnorm(x, mean=10000, sd=10)
real_dist <- dbinom(x, size=20000, p=.5)
posterior_dist <- calc_prob(x)

# make the data
data <- data.frame(x=rep(x, 2), freq=c(fake_dist, real_dist),
distribution=rep(c("Fake Normal(10000, 10)", "Real Binomial(20000, .5)"), each=length(x)))

# plot the data
g <- ggplot(data, aes(x=x, y=freq, group=distribution))
g + geom_line(aes(colour=distribution), size=1) + xlab("# of heads") + ylab("Density") +
opts(title="Distribution for fake and real experiment")
ggsave("cointoss.pdf")