I’m trying to wrap my head around a rather simple coin flip experiment.
Say I have studied a coin and obtained a posterior for the probability of getting a head when flipping it. How do I then use this posterior to get a prediction for how many heads I’ll see given n flips?
You want something called a Posterior Predictive Distribution. Your Bayesian model is a generative model. In essence, this means that if you knew the probability of the coin landing on heads (call it $\theta$), then you could generate predictions from your likelihood. Using the posterior to draw multiple $\theta$ and then passing them to the likelihood allows us to make predictions.
Let's see this in action. Let's assume we put a beta prior on our coin, observed some data, and now our posterior is something like $\operatorname{Beta}(10,12)$. Here is some python code to generate predictions or how many heads we might see in 100 flips of the same coin.
from scipy.stats import beta, binom
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Define the posterior to be a Beta(10,12)
posterior = beta(a=10, b=12)
def likelihood(theta, flips):
# Binomial Likelhood
return binom(n=flips, p=theta).rvs()
# Pass the posterior draws to the likelihood
predictions = likelihood(theta=posterior.rvs(1000), flips = 100)
plt.hist(predictions, edgecolor='white', bins = np.arange(10,100,10));
The result is the following:
We can then summarize our simulations using the mean (or some other summary function). So for what we've drawn here, we estimate that the number of heads in 100 flips would be approximately 45 or 46, and our prediction interval (depending on how you define it) is 23--46
Let's say you updated your prior with incoming data as you said, and obtained a posterior. Let's call it as $f_P(p)$. And call the number of heads in the next $n$ flips as $X$. Then, $P(X=k|p)$ is in binomial form. In order to get a final expression for $P(X=k)$, you need to integrate wrt posterior: $$P(X=k)=\int_{p\in\mathcal{P}}P(X=k|p)f_P(p)dp={n\choose k}\int_0^1 p^k(1-p)^{n-k}f_P(p)dp$$
If you've $f_P(p)$ in Beta form, this integral is going to have a trivial solution.
Let $\Theta$ be the RV for the probability of heads. Let`s denote the posterior with $f(\theta|X)$.
You now want $P(X_1=x_1, ..., X_n=x_n)$. where $X_i$ is the RV for outcome of the $i-th$ trial and $x_i \in \{0,1\}$.
Then $P(X_1=x_1, ..., X_n=x_n) = \int_\Theta \Pi_{i=1}^n\theta^{x_i}(1-\theta)^{1-x_i} f(\theta|X)d\theta$
This means that you "take a weighted average" over all possible values for $\theta$ and the weighting happends according to the posterior probability.
