Creating a posterior distribution for classic coin flipping in python using grid search

Question

I'm reading the book "Bayesian Analysis with Python" and the author provides some python code designed to show the grid search method of obtaining an approximate posterior distribution for the classic coin flipping example. In this example we set a prior on the probability of obtaining a head, then given some data our likelihood is the binomial distribution.

So from my understanding using grid search we will break the interval $[0,1]$ into chunks, in our case 100. We have a prior probability $P(\theta=\theta_{0})$ for each of the points in the discretised grid we multiply this by $P(D|\theta=\theta_{0})$ and we multiply these together to get an un-normalised estimate at each location.

The code is the following

def posterior_grid_approx(grid_points=100, heads=6, tosses=9):
    grid = np.linspace(0, 1, grid_points)
    prior=np.repeat(5,grid_points)
    likelihood=stats.binom.pmf(heads, tosses, grid)
    unstd_posterior = likelihood * prior
    posterior = unstd_posterior / unstd_posterior.sum()    
    return grid, posterior

I don't understand why the prior is set to a list of $5$'s isn't this meant to be a probability value? From my understanding they've broken the grid into 100 different points obtained a binomial pmf for each of these locations and then multiplied the number by 5. How is this a valid prior?

I'd really appreciate if someone could clearly explain the reasoning behind this as i like to understand the code as clearly as the concept. Thanks!

Answer 1

In short

The problem assumes a uniform prior distribution function. All possible $P(H)$ are equally likely. Because they are standardizing the probability distribution function at the end it does not matter what value is placed in

prior=np.repeat(1,grid_points)

---

Likelihood function

The likelihood function answers the question, how probable is the prior which we assumes.

For example, when flipping a fair coin we can assume that the probability of getting a heads is $P(H) = 0.5$. Now, if we do an experiment and in 3 flips we observe 3 heads then $P(HHH|P(H)=0.5) = 0.5^3$. This tells us that their is only a 0.125 chance that we selected the correct prior.

$Posterior \propto Likelihood \times Prior$

The code

As you observed it is strange that they set their entire list with the value $5$. However this will have no effect on the results because the before last line standardizes the results. Try changing the $5$ to a different value and you will see that your results remain unchanged.

Let's go through the code

grid = np.linspace(0, 1, grid_points)

This creates an array between $0$ and $1$. This is the different probabilities we will assume to calculate our likelihood. For example if we use 5 grid_points we have $[0, 0.25, 0.5, 0.75, 1]$. So we will assume that $P(H)$ is equal to these values in turn and get the resulting likelihood for each value.

likelihood=stats.binom.pmf(heads, tosses, grid)

We can now set out prior however we wish. We can assume that the coin is fair, or we can have no clue and just assume that it's a uniform distribution. This is essentially what they are doing here. The value they chose makes no difference in the results. However, I think for clarity making it equal to 1 is more obvious.

prior=np.repeat(1,grid_points)

To calculate the posterior we multiply the prior by the likelihood

unstd_posterior = likelihood * prior

And then to get a probability distribution function we standardize the function.

posterior = unstd_posterior / unstd_posterior.sum()