Interactive_Sivia_Coin_Toss.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Interactive Bayesian Coin Demonstration from Sivia (1996)\n",
    "\n",
    "### Sivia, D.S., 1996, Data Analysis: A Bayesian Tutorial\n",
    "\n",
    "* interactive plot demonstration with ipywidget package\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "#### The Bayesian Coin Example\n",
    "\n",
    "I have a coin and you need to figure out if it is a fair coin!\n",
    "\n",
    "* a fair coin would have a 50% probability of heads (and a 50% probability of tails)\n",
    "\n",
    "You start with your prior assessment of my coin, a prior distribution over the probability of heads $Prob(Coin)$\n",
    "\n",
    "* it could be based on how honest you think I am\n",
    "\n",
    "Then you perform a set of coin tosses to build a likelihood distribution, $P(Tosses | Coin)$\n",
    "\n",
    "* the more coin tosses, the narrower this distribution\n",
    "\n",
    "Then you update the prior distribution with the likelihood distribution to get the posterior distribution, $P(Coin | Tosses)$.\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "P( Coin | Tosses ) = \\frac{P( Tosses | Coin ) P( Coin )}{P( Tosses )}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "#### Objective \n",
    "\n",
    "Provide an example and demonstration for:\n",
    "\n",
    "1. interactive plotting in Jupyter Notebooks with Python packages matplotlib and ipywidgets\n",
    "2. provide an intuitive hands-on example of Bayesian updating   \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "#### Load the Required Libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interactive                        # widgets and interactivity\n",
    "from ipywidgets import widgets                            # widgets and interactivity\n",
    "import matplotlib; import matplotlib.pyplot as plt        # plotting\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "plt.rc('axes', axisbelow=True)                  # set axes and grids in the background for all plots\n",
    "import numpy as np                                        # working with arrays\n",
    "from scipy.stats import triang                            # parametric distributions\n",
    "from scipy.stats import binom\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make Our Interactive Plot \n",
    "\n",
    "For this demonstration we will: \n",
    "\n",
    "* declare a set of 4 widgets in a HBox (horizontal box of widgets). \n",
    "\n",
    "\n",
    "* define a function 'f' that will read the output from these widgets and make a plot\n",
    "\n",
    "You may have some flicker and lag.  I have not tried to optimize performance for this demonstration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 4 slider bars for the model input\n",
    "a = widgets.FloatSlider(min=0.0, max = 1.0, value = 0.5, description = 'coin bias',continuous_update=False)\n",
    "d = widgets.FloatSlider(min=0.01, max = 1.0, value = 0.1, step = 0.01, description = 'coin uncert.',continuous_update=False)\n",
    "b = widgets.FloatSlider(min = 0, max = 1.0, value = 0.5, description = 'prop. heads',continuous_update=False)\n",
    "c = widgets.IntSlider(min = 5, max = 1000, value = 100, description = 'coin tosses',continuous_update=False)\n",
    "\n",
    "ui = widgets.HBox([a,d,b,c],)\n",
    "\n",
    "def f(a, b, c, d):                                        # function to make the plot   \n",
    "    heads = int(c * b)\n",
    "    tails = c - heads\n",
    "    \n",
    "    x = np.linspace(0.0, 1.0, num=1000)\n",
    "    \n",
    "    prior = norm.pdf(x,loc = a, scale = d)\n",
    "    prior = prior / np.sum(prior)\n",
    "    \n",
    "    plt.subplot(221)\n",
    "    \n",
    "    plt.plot(x, prior,color='red',alpha=0.8,linewidth=3)  # prior distribution of coin fairness\n",
    "    plt.arrow(0.48,0.04,-0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Tails',[0.10,0.042],rotation=0.0)\n",
    "    plt.arrow(0.52,0.04,0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Heads',[0.55,0.042],rotation=0.0)\n",
    "    plt.annotate('Double Tailed Coin',[0.02,0.02],rotation=90.0)\n",
    "    plt.annotate('Double Headed Coin',[0.97,0.02],rotation=90.0)\n",
    "    plt.xlim(0.0,1.0)\n",
    "    plt.xlabel('Coin Bias = P(Coin Heads)'); plt.ylabel('Density'); plt.title(r'$Prior$ $Distribution$ - What do you think of the coin before the tosses?')\n",
    "    plt.ylim(0, 0.05)\n",
    "    plt.axvline(x=0.5,color='grey',linestyle='--');\n",
    "    plt.annotate('Fair Coin',[0.49,0.02],rotation=90.0)\n",
    "    plt.grid()\n",
    "    plt.gca().xaxis.grid(True, which='major',linewidth = 1.0); plt.gca().xaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "    \n",
    "    plt.subplot(222)                                      # results from the coin tosses \n",
    "    plt.pie([heads, tails],labels = ['heads','tails'],radius = 0.5*(c/1000)+0.5, autopct='%1.1f%%', colors = ['#ff9999','#66b3ff'], explode = [.02,.02], wedgeprops = {\"edgecolor\":\"k\",'linewidth': 1}  )\n",
    "    plt.title(str(c) + ' Coin Tosses')\n",
    "    \n",
    "    likelihood = binom.pmf(heads,c,x)\n",
    "    likelihood = likelihood/np.sum(likelihood)\n",
    "    \n",
    "    plt.subplot(223)                                      # likelihood distribution given the coin tosses\n",
    "    plt.plot(x, likelihood,color='red',alpha=0.8,linewidth=3)\n",
    "    plt.arrow(0.48,0.04,-0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Tails',[0.10,0.042],rotation=0.0)\n",
    "    plt.arrow(0.52,0.04,0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Heads',[0.55,0.042],rotation=0.0)\n",
    "    plt.annotate('Double Tailed Coin',[0.02,0.02],rotation=90.0)\n",
    "    plt.annotate('Double Headed Coin',[0.97,0.02],rotation=90.0)\n",
    "    plt.xlim(0.0,1.0); plt.ylim(0, 0.05)\n",
    "    plt.xlabel('P(Tosses | Coin Bias)'); plt.ylabel('Density'); plt.title(r'$Likelihood$ $Distribution$ - What do the coin tosses tell you about the coin?')\n",
    "    plt.axvline(x=0.5,color='grey',linestyle='--');\n",
    "    plt.annotate('Fair Coin',[0.49,0.02],rotation=90.0)\n",
    "    plt.grid()\n",
    "    plt.gca().xaxis.grid(True, which='major',linewidth = 1.0); plt.gca().xaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "    \n",
    "    post = prior * likelihood\n",
    "    post = post / np.sum(post)    \n",
    "    \n",
    "    plt.subplot(224)                                      # posterior distribution\n",
    "    plt.plot(x, post,color='red',alpha=0.8,linewidth=3)\n",
    "    plt.arrow(0.48,0.04,-0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Tails',[0.10,0.042],rotation=0.0)\n",
    "    plt.arrow(0.52,0.04,0.4,0,linewidth=1,facecolor='grey',edgecolor='black',head_length=.03)\n",
    "    plt.annotate('Increasingly Weighted Heads',[0.55,0.042],rotation=0.0)\n",
    "    plt.annotate('Double Tailed Coin',[0.02,0.02],rotation=90.0)\n",
    "    plt.annotate('Double Headed Coin',[0.97,0.02],rotation=90.0)\n",
    "    plt.xlim(0.0,1.0); plt.ylim(0, 0.05)\n",
    "    plt.xlabel('P(Coin Bias | Tosses)'); plt.ylabel('Density'); plt.title(r\"$Posterior$ $Distribution$ - What do you think of the coin after the tosses?\")    \n",
    "    plt.axvline(x=0.5,color='grey',linestyle='--');\n",
    "    plt.annotate('Fair Coin',[0.49,0.02],rotation=90.0)\n",
    "    plt.grid()\n",
    "    plt.gca().xaxis.grid(True, which='major',linewidth = 1.0); plt.gca().xaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.6, wspace=0.2, hspace=0.3)\n",
    "    plt.show()\n",
    "\n",
    "interactive_plot = widgets.interactive_output(f, {'a': a, 'd': d, 'b': b, 'c': c})\n",
    "interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Coin Example from Sivia, 1996, Data Analysis: A Bayesian Tutorial\n",
    "* interactive plot demonstration with ipywidget package\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Problem\n",
    "\n",
    "What is the PDF for the coin probability of heads, P(Coin Heads)?  Start with a prior model and update with coin tosses.\n",
    "\n",
    "* **coin bias**: expectation for your prior distribution for probability of heads\n",
    "\n",
    "* **coin uncert.**: standard deviation for your prior distribution for probability of heads\n",
    "\n",
    "* **prop. heads**: proportion of heads in the coin toss experiment\n",
    "\n",
    "* **coin tosses**: number of coin tosses in thecoin toss experiment\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1d49ebdc2101493abe1fb479d5afb41f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatSlider(value=0.5, continuous_update=False, description='coin bias', max=1.0), FloatSlider(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e5c0aac98ce54dc88df1b684b07dc0fe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a simple demonstration of interactive plots in Jupyter Notebook Python with the ipywidgets and matplotlib packages. \n",
    "\n",
    "I have many other demonstrations on data analytics and machine learning, e.g. on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}