{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Demonstration: Bayesian Coin Tossing"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Bayesian updating\n",
"\n",
"$ \\newcommand{\\thetavec}{\\boldsymbol{\\theta}} \\newcommand{\\pr}{\\textrm{p}}$\n",
"Recall Bayes' theorem with $\\thetavec$ the vector of parameters we seek and information $I$ is kept implicit.\n",
"\n",
"$$\n",
" \\overbrace{\\pr(\\thetavec \\mid \\textrm{data},I)}^{\\textrm{posterior}} =\n",
" \\frac{\\color{red}{\\overbrace{\\pr(\\textrm{data} \\mid \\thetavec,I)}^{\\textrm{likelihood}}} \\times\n",
" \\color{blue}{\\overbrace{\\pr(\\thetavec \\mid I)}^{\\textrm{prior}}}}\n",
" {\\color{darkgreen}{\\underbrace{\\pr(\\textrm{data} \\mid I)}_{\\textrm{evidence}}}}\n",
"$$\n",
"\n",
"If we view the prior as the initial information we have about $\\thetavec$, summarized as a probability density function, then Bayes' theorem tells us how to update that information after observing some data: this is the posterior pdf. Here we will give some examples of how this plays out when tossing coins.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Determining the bias of a coin"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The idea here is that we are observing successive flips of a coin, which is a proxy for any process that has a binary outcome. There is a definite true probability for getting heads, which we'll label $p_h$, but we don't know what it is. We start with a preconceived notion of the probability expressed in terms of a prior pdf for $p_h$, i.e., $\\pr(p_h)$. With each flip of the coin, we have more information, so our goal is to update our expectation of $p_h$, meaning we want the posterior $p(p_h \\vert \\mbox{\\# tosses, \\# heads})$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"rng = np.random.default_rng(seed=999) # for reproducibility\n",
"pHtrue=0.6 # biased coin\n",
"flips=rng.random(2**12) # simulates 4096 coin flips\n",
"heads=flips"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pH=np.linspace(0,1,1000)\n",
"fig, axs = plt.subplots(nrows=4,ncols=3,sharex=True,sharey='row',figsize=(14,14))\n",
"axs_vec=np.reshape(axs,-1)\n",
"axs_vec[0].plot(pH,prior(pH))\n",
"for ndouble in range(11):\n",
" ax=axs_vec[1+ndouble]\n",
" ax.plot(pH,posterior(pH,heads[:2**ndouble]))\n",
" ax.text(0.1, 0.8, '$N={0}$'.format(2**ndouble), transform=ax.transAxes)\n",
"for row in range(4): axs[row,0].set_ylabel('$p(p_H|D_\\mathrm{obs},I)$')\n",
"for col in range(3): axs[-1,col].set_xlabel('$p_H$')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Please, feel free to play with this setup yourself by modifying the code above. Things to explore:\n",
"* The choice of `pH`.\n",
"* The number of data.\n",
"* The choice of prior (i.e. the function `prior`).\n",
"* Is there a difference between the results obtained for many sequential updates (i.e. updating the posterior after each toss), and a single update with all data? Can you show this result analytically?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Degree of belief intervals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are going to compute some Bayesian confidence intervals, aka DoB intervals, aka credibility intervals....\n",
"You shoud go through this section once using the numbers provided, and then come back and run it again. You can also employ priors other than those provided if you wish."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import scipy.stats as stats"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"def dist_stuff(dist):\n",
" \"\"\"\n",
" Finds the median, mean, and 68%/95% credible intervals for the given \n",
" 1-d distribution (which is an object from scipy.stats). \n",
" \"\"\"\n",
" # For x = median, mean: return x and the value of the pdf at x as a list\n",
" median = [dist.median(), dist.pdf(dist.median())] \n",
" mean = [dist.mean(), dist.pdf(dist.mean())]\n",
" # The left and right limits of the credibility interval are returned\n",
" cred68 = dist.interval(0.68)\n",
" cred95 = dist.interval(0.95)\n",
" return median, mean, cred68, cred95"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"def dist_mode(dist, x):\n",
" \"\"\"\n",
" Return the mode (maximum) of the 1-d distribution for array x.\n",
" \"\"\"\n",
" x_max_index = dist.pdf(x).argmax()\n",
" # Return x of the maximum and the value of the pdf at that x \n",
" mode = [x[x_max_index], dist.pdf(x[x_max_index])]\n",
" return mode"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Then we use this to write an object that will give us back the mean, 68%, and 95% intervals for a uniform prior. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def print_uniform_prior_measures(N,heads):\n",
" \"\"\"\n",
" Prints out the mean, and 68 and 95 CIs for a uniform prior. Note that this means alpha=beta=1.\n",
" \"\"\"\n",
" x = np.linspace(0, 1, 301) # mesh for posterior plots (enough so smooth)\n",
" median, mean, cred68, cred95 = dist_stuff(stats.beta(1+heads,1+N-heads))\n",
" mode=dist_mode(stats.beta(1+heads,1+N-heads),x)\n",
" print(f'For a uniform prior, and {heads} heads out of {N} tosses:')\n",
" print(f' mean={mean[0]:.3f}; Mode={mode[0]:.3f}')\n",
" print(f' 68% DoB interval=[{cred68[0]:.3f},{cred68[1]:.3f}]')\n",
" print(f' 95% DoB interval=[{cred95[0]:.3f},{cred95[1]:.3f}]')\n",
" return "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Now we fill in the values for N and heads from running the widget. "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"For a uniform prior, and 4 heads out of 10 tosses:\n",
" mean=0.417; Mode=0.400\n",
" 68% DoB interval=[0.276,0.558]\n",
" 95% DoB interval=[0.167,0.692]\n"
]
}
],
"source": [
"num_coin_tosses=10\n",
"print_uniform_prior_measures(num_coin_tosses,np.sum(heads[:num_coin_tosses]))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The actual value of p_H is 0.6\n"
]
}
],
"source": [
"print('The actual value of p_H is',pHtrue)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"So, at least for the data I generated, the 68% DoB for a uniform prior does not contain the true value. Should you be surprised? "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"def print_frequentist_estimators(N,heads):\n",
" \"\"\"\n",
" Finds the mean, and 68 and 95 CIs. Note that this means alpha_1=beta_1=1.\n",
" \"\"\"\n",
" mean=heads/N\n",
" sigma=np.sqrt(mean*(1-mean)/N)\n",
" print(f'For {heads} heads out of {N} tosses, the frequentist 68% confidence interval=',\\\n",
" f'({mean-sigma:.3f},{mean+sigma:.3f})')\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"For 4 heads out of 10 tosses, the frequentist 68% confidence interval= (0.245,0.555)\n"
]
}
],
"source": [
"print_frequentist_estimators(num_coin_tosses,np.sum(heads[:num_coin_tosses]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Which Bayesian estimator is the frequentist mean closest to?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Is the frequentist 68% confidence interval the same as the Bayesian 68% DoB interval? If so, should they be? If not, why are they different?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Now we will also generate the summary statistics for the other priors. (What is coded is for the default values. After running through the exercise you can come back and try and change it; indeed, you should do that if you comparing to results where you altered the prior above.)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def print_likely_fair_prior_measures(N,heads):\n",
" \"\"\"\n",
" Prints out the mean, and 68 and 95 CIs for the \"coin is likely fair\" prior. This means alpha=beta=30.\n",
" \"\"\"\n",
" x = np.linspace(0, 1, 301) # mesh for posterior plots (enough so smooth)\n",
" median, mean, cred68, cred95 = dist_stuff(stats.beta(30+heads,30+N-heads))\n",
" mode=dist_mode(stats.beta(1+heads,1+N-heads),x)\n",
" print('For the trusting-person`s prior, and',heads,'heads out of',N,'tosses:')\n",
" print (f'Mean={mean[0]:.3f}; Mode={mode[0]:.3f}')\n",
" print (f'68% DoB interval=[{cred68[0]:.3f},{cred68[1]:.3f}]')\n",
" print (f'95% DoB interval=[{cred95[0]:.3f},{cred95[1]:.3f}]')\n",
" return "
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"For the trusting-person`s prior, and 4 heads out of 10 tosses:\n",
"Mean=0.486; Mode=0.400\n",
"68% DoB interval=[0.426,0.545]\n",
"95% DoB interval=[0.370,0.602]\n"
]
}
],
"source": [
"print_likely_fair_prior_measures(num_coin_tosses,np.sum(heads[:num_coin_tosses]))"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def print_likely_unfair_prior_measures(N,heads):\n",
" \"\"\"\n",
" Prints out the mean, and 68 and 95 CIs for the \"coin is likely unfair\" prior. This means alpha=beta=0.2.\n",
" \"\"\"\n",
" x = np.linspace(0, 1, 301) # mesh for posterior plots (enough so smooth)\n",
" median, mean, cred68, cred95 = dist_stuff(stats.beta(0.2+heads,0.2+N-heads))\n",
" mode=dist_mode(stats.beta(1+heads,1+N-heads),x)\n",
" print('For the nasty-suspicious-mind prior, and',heads,'heads out of',N,'tosses:')\n",
" print (f'Mean={mean[0]:.3f}; Mode={mode[0]:.3f}')\n",
" print (f'68% DoB interval=[{cred68[0]:.3f},{cred68[1]:.3f}]')\n",
" print (f'95% DoB interval=[{cred95[0]:.3f},{cred95[1]:.3f}]')\n",
" return "
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"For the nasty-suspicious-mind prior, and 4 heads out of 10 tosses:\n",
"Mean=0.404; Mode=0.400\n",
"68% DoB interval=[0.253,0.555]\n",
"95% DoB interval=[0.144,0.699]\n"
]
}
],
"source": [
"print_likely_unfair_prior_measures(num_coin_tosses,np.sum(heads[:num_coin_tosses]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"So what is the best approach in this case? Objectivity? Trust? Suspicion?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Now having printed out the results for my particular coin-tossing experiment you should play with things and see what the different summary statistics give for other \"experimental runs\".\n",
"\n",
"Don't forget to change the seed for the random number generator, otherwise you will get exactly the same sequence of coin tosses."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}