Instantly share code, notes, and snippets.

# ctufts/index.html

Last active July 9, 2017 23:23
Star You must be signed in to star a gist
coin flip example
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 Coin Flip Problem

I’m not the best at calculus (or at least it’s not the most natural way for me to understand things), so I went about the problem using Gibbs sampling (MCMC). The uniform prior was assumed in the other responses, which is the same as a beta distribution with parameters a = 1 and b = 1. The plot below is targeted to provide a sense of what a beta distribution looks like as you adjust parameters. In simple terms, a beta distribution is a probability distribution of probabilities. Think of it like this, if I’ve flipped a coin 100 times and 51 were heads and 49 were tails, I’d be pretty confident the coin is fair. If you look at the plots below, as the alpha and beta values increase the distribution becomes more narrow, i.e. the probability of theta being 0.5 goes up.

library(ggplot2)
a <- c(1, 3, 10)
b <- c(1, 3, 10)
params <- cbind(a,b)
ds <- NULL
n <- seq(0,1,0.01)
for(i in 1:nrow(params)){
ds <- rbind(data.frame(x = n, y = dbeta(n, params[i,1], params[i,2]),
parameters = paste0("\U03B1 = ",params[i,1],
", \U03B2 = ", params[i,2])), ds)
}

ggplot(ds, aes(x = x, y = y, color=parameters)) + geom_line() +
labs(x = '\U03B8', y = 'Probability Density') +
scale_color_discrete(name=NULL)

So, why do you need this beta distribution? Primarily due to the method I’m going to use - Gibbs sampling. Gibbs sampling allows for an end run around that ‘evidence’ term (i.e. marginal probability). In this case it is an easily solveable integral (for those good at calculus), but I’m more of a addition/subtraction type person. I’m not going to go into detail about how conjugate priors (beta to bernoulli) and gibbs sampling work - if you are interested go here for more info.

Let’s simulate 10k iterations using a uniform prior (i.e. beta distribution with parameters 1 and 1).

# uniform prior (i.e. beta of 1, 1)
a = 1
b = 1
N1 = 10 # total flips

niters = 50000
burnin = 1000

thetas = rep(0, nrow = (niters-burnin))
for (i in 1:niters){
theta = c(rbeta(n = 1, shape1 = a + z1, shape2 = b + N1 - z1))
if (i >= burnin){
thetas[(i-burnin)] = theta
}
}

p_1_1 <- mean(thetas > 0.45 & thetas < 0.55)

ggplot(data.frame(thetas), aes(x = thetas)) +
geom_histogram(fill="#E4002B", color = "#7A99AC") +
scale_x_continuous(limits = c(0,1))

We get a probability of 12.88% that the coin meets the ‘relatively fair’ contraints of 0.45 to 0.55. Now lets’s adjust those priors. You might ask why adjust the priors? Unless someone constructed this coin to specifically be loaded to the point that it never goes heads, it’s probably not the best assumption that the coin has the same probability of being fair as it does being loaded to never be heads. (In a practical sense: If I’m trying to cheat I’m not going to be greedy, have to lose a few flips here an there, otherwise I’ll get caught)

Let’s change up to a beta distribution of 3 and 3.

a = 3
b = 3
N1 = 10 # total flips

niters = 50000
burnin = 1000

thetas = rep(0, nrow = (niters-burnin))
for (i in 1:niters){
theta = c(rbeta(n = 1, shape1 = a + z1, shape2 = b + N1 - z1))
if (i >= burnin){
thetas[(i-burnin)] = theta
}
}

p_3_3 <- mean(thetas > 0.45 & thetas < 0.55)

ggplot(data.frame(thetas), aes(x = thetas)) +
geom_histogram(fill="#E4002B", color = "#7A99AC") +
scale_x_continuous(limits = c(0,1))