Interactive_Confidence_Interval.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 Confidence Interval Demonstration\n",
    "\n",
    "### Boostrap and Analytical Confidence Intervals\n",
    "\n",
    "* we calculate the confidence interval in the mean with boostrap and compare to the analytical expression\n",
    "\n",
    "* with this workflow we all provide an interactive plot demonstration with matplotlib and ipywidget packages\n",
    "\n",
    "#### Michael Pyrcz, 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",
    "#### Confidence Intevals\n",
    "\n",
    "Confidence intervals are the uncertainty in a sample statistic or model parameter\n",
    "\n",
    "* for uncertainty in the sample mean we have: \n",
    "\n",
    "    * center on the sample proportion, $\\hat{p}$\n",
    "\n",
    "    * the standard error in the proportion for the dispersion (spread)\n",
    "\n",
    "    * Student's t distributed for small samples and Gaussian distributed for large sample sizes\n",
    "    \n",
    "The analytical form for small samples ($n \\lt 30$) is:\n",
    "\n",
    "\\begin{equation}\n",
    "CI: \\hat{p} \\pm t_{\\frac{\\alpha}{2},n-1} \\times \\frac {\\sqrt{p(1-p)}}{\\sqrt{n}}\n",
    "\\end{equation}\n",
    "\n",
    "where the sampling distribution of the proportion is student's t distributed with number of samples, $n$ - 1, degrees of freedom and $\\alpha$ is the signficance level divded by 2 for the two tails.\n",
    "\n",
    "When the number of samples is large ($n \\ge 30$) then the analytical form converges to Gaussian distributed:\n",
    "\n",
    "\\begin{equation}\n",
    "CI: \\hat{p} \\pm N_{\\frac{\\alpha}{2}} \\times \\frac {\\sqrt{p(1-p)}}{\\sqrt{n}}\n",
    "\\end{equation}\n",
    "\n",
    "#### Bootstrap\n",
    "\n",
    "Uncertainty in the sample statistics\n",
    "* one source of uncertainty is the paucity of data.\n",
    "* do 200 or even less wells provide a precise (and accurate estimate) of the mean? standard deviation? skew? P13?\n",
    "\n",
    "Would it be useful to know the uncertainty in these statistics due to limited sampling?\n",
    "* what is the impact of uncertainty in the mean porosity e.g. 20%+/-2%?\n",
    "\n",
    "**Bootstrap** is a method to assess the uncertainty in a sample statistic by repeated random sampling with replacement.\n",
    "\n",
    "Assumptions\n",
    "* sufficient, representative sampling, identical, idependent samples\n",
    "\n",
    "Limitations\n",
    "1. assumes the samples are representative \n",
    "2. assumes stationarity\n",
    "3. only accounts for uncertainty due to too few samples, e.g. no uncertainty due to changes away from data\n",
    "4. does not account for boundary of area of interest \n",
    "5. assumes the samples are independent\n",
    "6. does not account for other local information sources\n",
    "\n",
    "The Bootstrap Approach (Efron, 1982)\n",
    "\n",
    "Statistical resampling procedure to calculate uncertainty in a calculated statistic from the data itself.\n",
    "* Does this work?  Prove it to yourself, for uncertainty in the mean solution is standard error: \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_\\overline{x} = \\frac{\\sigma^2_s}{n}\n",
    "\\end{equation}\n",
    "\n",
    "Extremely powerful - could calculate uncertainty in any statistic!  e.g. P13, skew etc.\n",
    "* Would not be possible access general uncertainty in any statistic without bootstrap.\n",
    "* Advanced forms account for spatial information and sampling strategy (game theory and Journel’s spatial bootstrap (1993).\n",
    "\n",
    "Steps: \n",
    "\n",
    "1. assemble a sample set, must be representative, reasonable to assume independence between samples\n",
    "\n",
    "2. optional: build a cumulative distribution function (CDF)\n",
    "    * may account for declustering weights, tail extrapolation\n",
    "    * could use analogous data to support\n",
    "\n",
    "3. For $\\ell = 1, \\ldots, L$ realizations, do the following:\n",
    "\n",
    "    * For $i = \\alpha, \\ldots, n$ data, do the following:\n",
    "\n",
    "        * Draw a random sample with replacement from the sample set or Monte Carlo simulate from the CDF (if available). \n",
    "\n",
    "6. Calculate a realization of the sammary statistic of interest from the $n$ samples, e.g. $m^\\ell$, $\\sigma^2_{\\ell}$. Return to 3 for another realization.\n",
    "\n",
    "7. Compile and summarize the $L$ realizations of the statistic of interest.\n",
    "\n",
    "This is a very powerful method.  Let's try it out and compare the result to the analytical form of the confidence interval for the sample mean. \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 confidence intervals and compare to statistical boostrap   \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                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "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",
    "import pandas as pd                                       # working with DataFrames\n",
    "import seaborn as sns                                     # for matrix scatter plots\n",
    "from scipy.stats import triang                            # parametric distributions\n",
    "from scipy.stats import binom\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import uniform\n",
    "from scipy.stats import triang\n",
    "from scipy.stats import t\n",
    "from scipy import stats                                   # statistical calculations\n",
    "import random                                             # random drawing / bootstrap realizations of the data\n",
    "from matplotlib.gridspec import GridSpec                  # nonstandard subplots\n",
    "import math                                               # square root operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Declare functions\n",
    "\n",
    "These are some of the functions from GeostatsPy required by the new program.  We will declare them here and then in the future integrate the new indicator kriging program into the package properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_grid():\n",
    "    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().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   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make a Synthetic Dataset\n",
    "\n",
    "This is an interactive method to:\n",
    "\n",
    "* select a parametric distribution\n",
    "\n",
    "* select the distribution parameters\n",
    "\n",
    "* select the number of samples and visualize the synthetic dataset distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the synthetic dataset\n",
    "bins = np.linspace(0,1000,1000)\n",
    "\n",
    "# interactive calculation of the sample set (control of source parametric distribution and number of samples)\n",
    "l = widgets.Text(value='                                Simple Boostrap Demonstration, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "a = widgets.IntSlider(min=0, max = 100, value = 2, step = 1, description = '$n_{red}$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "a.style.handle_color = 'red'\n",
    "\n",
    "b = widgets.IntSlider(min=0, max = 100, value = 3, step = 1, description = '$n_{green}$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "b.style.handle_color = 'green'\n",
    "\n",
    "c = widgets.IntSlider(min=1, max = 16, value = 3, step = 1, description = '$L$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "c.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([a,b,c],)                               # basic widget formatting           \n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def f_make(a, b, c):                                      # function to take parameters, make sample and plot\n",
    "    red_freq = make_data(a, b, c)\n",
    "\n",
    "    labels = ['Red', 'Green']\n",
    "    nrows = int(np.sqrt(c)+0.9); ncols = int(c / nrows + 0.9)\n",
    "    print(nrows,ncols)\n",
    "    plt.clf()\n",
    "    \n",
    "    for i in range(0, c): \n",
    "        plt.subplot(ncols,nrows,i + 1) \n",
    "        draw = [red_freq[i],a + b - red_freq[i]]\n",
    "        plt.grid(zorder=0, color='black', axis = 'y', alpha = 0.2); plt.ylim(0,a + b); \n",
    "        plt.ylabel('Frequency'); plt.xlabel('Balls Drawn')\n",
    "        plt.yticks(np.arange(0,a + b + 1,max(1,round((a+b)/10))))\n",
    "        barlist = plt.bar(labels,draw,edgecolor = \"black\",linewidth = 1,alpha = 0.8); plt.title('Realization #' + str(i+1),zorder = 1) \n",
    "        barlist[0].set_color('r'); barlist[1].set_color('g')\n",
    "            \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=ncols, top=nrows + 0.2 * nrows, wspace=0.2, hspace=0.2)\n",
    "    plt.show()\n",
    "\n",
    "def make_data(a, b, c):                                   # function to check parameters and make sample   \n",
    "    prop_red = np.zeros(c)\n",
    "    for i in range(0, c): \n",
    "        prop_red[i] = np.random.multinomial(a+b,[a/(a+b),b/(a+b)], size = 1)[0][0]\n",
    "    return prop_red\n",
    "\n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(f_make, {'a': a, 'b': b, 'c': c})\n",
    "interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Bootstrap Demonstration - Drawing Red and Green Balls from a Virtual Cowboy Hat\n",
    "\n",
    "* drawing red and green balls from a hat with replacement to access uncertainty in the proportion \n",
    "\n",
    "* interactive plot demonstration with ipywidget, matplotlib packages\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",
    "Let's simulate bootstrap, resampling with replacement from a hat with $n_{red}$ and $n_{green}$ balls\n",
    "\n",
    "* **$n_{red}$**: number of red balls in the sample (placed in the hat)\n",
    "\n",
    "* **$n_{green}$**: number of green balls in the sample (placed in the hat)\n",
    "\n",
    "* **$L$**: number of bootstrap realizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b1e8a58ded314882964e8cc7ecb0a23c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                Simple Boostrap Demonstration, Michael Pyrcz, Assoc…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a353af7b01884706bcc7bb4f522ff074",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'stream', 'text': '2 2\\n', 'name': 'stdout'}, {'output_type': 'display_data', …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summarizing Bootstrap Uncertainty\n",
    "\n",
    "* Run more bootstrap realizations and evaluate the uncertianty model\n",
    "\n",
    "Now instead of looking at each bootstrap result, let's make many and summarize with:\n",
    "\n",
    "* **box and whisker plot** of the red and green ball frequencies\n",
    "\n",
    "* **histograms** of the red and green ball frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the synthetic dataset\n",
    "bins = np.linspace(0,1000,1000)\n",
    "\n",
    "# interactive calculation of the sample set (control of source parametric distribution and number of samples)\n",
    "l2 = widgets.Text(value='     Confidence Interval for Proportions, Analytical and Bootstrap Demonstration, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "a2 = widgets.IntSlider(min=0, max = 100, value = 20, step = 1, description = '$n_{red}$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "a2.style.handle_color = 'red'\n",
    "\n",
    "b2 = widgets.IntSlider(min=0, max = 100, value = 30, step = 1, description = '$n_{green}$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "b2.style.handle_color = 'green'\n",
    "\n",
    "c2 = widgets.IntSlider(min=5, max = 1000, value = 1000, step = 1, description = '$L$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "c2.style.handle_color = 'gray'\n",
    "\n",
    "alpha = widgets.FloatSlider(min=0.01, max = 0.40, value = 0.05, step = 0.01, description = r'$\\alpha$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "alpha.style.handle_color = 'gray'\n",
    "\n",
    "uib = widgets.HBox([a2,b2,c2,alpha],)                            # basic widget formatting           \n",
    "uib2 = widgets.VBox([l2,uib],)\n",
    "\n",
    "def s_make(a, b, c, alpha):                                      # function to take parameters, make sample and plot\n",
    "    \n",
    "    dof = a + b - 1\n",
    "    \n",
    "    red_freq = make_data(a, b, c)\n",
    "    pred = red_freq/(a+b)\n",
    "    red_prop = (a / (a+b))\n",
    "    red_SE = math.sqrt((red_prop * (1.0 - red_prop)) / (a+b))\n",
    "    \n",
    "    green_freq = (a + b) - red_freq\n",
    "    pgreen = green_freq/(a+b)\n",
    "    green_prop = (b / (a+b))\n",
    "    green_SE = math.sqrt((green_prop * (1.0 - green_prop)) / (a+b))\n",
    "    \n",
    "    prop_red = red_freq / (a + b)\n",
    "    prop_green = green_freq / (a + b)\n",
    "    labels = ['Red Balls', 'Green Balls']\n",
    "    bins = np.linspace(0,a + b, a + b)\n",
    "    \n",
    "    fig = plt.figure(constrained_layout=False)\n",
    "    gs = GridSpec(3, 2, figure=fig)\n",
    "    \n",
    "    ax1 = fig.add_subplot(gs[:, 0])\n",
    "    boxplot = ax1.boxplot([pred,pgreen],labels = labels, notch = True, sym = '+',patch_artist=True)\n",
    "    colors = ['red','green']\n",
    "    for patch, color in zip(boxplot['boxes'], colors):\n",
    "        patch.set_facecolor(color)\n",
    "    for patch, color in zip(boxplot['medians'], colors):\n",
    "        patch.set_color('black')\n",
    "    \n",
    "    ax1.set_ylim([0,1])\n",
    "    ax1.grid(zorder=0, color='black', axis = 'y', alpha = 0.2)\n",
    "    ax1.set_ylabel('Proportion of Balls'); ax1.set_xlabel('Ball Color');ax1.set_title('Bootstrap Uncertainty - Proportion Distributions')\n",
    "    ax1.grid(True, which='major',axis='y',linewidth = 1.0); ax1.grid(True, which='minor',axis='y',linewidth = 0.2) # add y grids\n",
    "    ax1.tick_params(which='major',length=7); ax1.tick_params(which='minor', length=4)\n",
    "    ax1.xaxis.set_minor_locator(AutoMinorLocator()); ax1.yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "\n",
    "    cumul_prob = np.linspace(0.0,1.0,100)\n",
    "    if a <= 30 or b <= 30:\n",
    "        red_prop_values = t.ppf(cumul_prob, dof)\n",
    "        red_lower = t.ppf(alpha/2, dof); red_upper = t.ppf(1-alpha/2, dof)\n",
    "    else:\n",
    "        red_prop_values = norm.ppf(cumul_prob)\n",
    "        red_lower = norm.ppf(alpha/2); red_upper = norm.ppf(1-alpha/2)\n",
    "    \n",
    "    red_prop_values = red_prop_values * red_SE + red_prop \n",
    "    red_lower = red_lower * red_SE + red_prop\n",
    "    red_upper = red_upper * red_SE + red_prop\n",
    "    \n",
    "    cumul_prob = np.linspace(0.01,0.99,100)\n",
    "    if a <= 30 or b <= 30:\n",
    "        green_prop_values = t.ppf(cumul_prob, dof)\n",
    "        green_lower = t.ppf(alpha/2, dof); green_upper = t.ppf(1-alpha/2, dof)\n",
    "    else:\n",
    "        green_prop_values = norm.ppf(cumul_prob)\n",
    "        green_lower = norm.ppf(alpha/2); green_upper = norm.ppf(1-alpha/2)\n",
    "        \n",
    "    green_prop_values = green_prop_values * green_SE + green_prop \n",
    "    green_lower = green_lower * green_SE + green_prop\n",
    "    green_upper = green_upper * green_SE + green_prop    \n",
    "\n",
    "    ax2 = fig.add_subplot(gs[0, 1])\n",
    "    ax2.hist(prop_red,cumulative = True, density = True, alpha=0.7,color=\"red\",edgecolor=\"black\",linewidth=2,bins = np.linspace(0,1,50), label = 'Bootstrap')\n",
    "    ax2.plot([red_lower,red_lower],[0,1],color='black',linewidth=2,linestyle='--',label='Lower/Upper')\n",
    "    ax2.plot([red_upper,red_upper],[0,1],color='black',linewidth=2,linestyle='--')\n",
    "    ax2.plot([red_prop,red_prop],[0,1],color='black',linewidth=3,label='Exp.')\n",
    "    ax2.set_title('Uncertainty in Proportion of Red Balls'); ax2.set_xlabel('Proportion of Red Balls'); ax2.set_ylabel('Cumulative Probability') \n",
    "    ax2.set_xlim([0,1]); ax2.set_ylim([0,1])   \n",
    "    \n",
    "    ax2.plot(red_prop_values, cumul_prob, color = 'black', linewidth = 2, label = 'Analytical')\n",
    "    ax2.legend()\n",
    "    \n",
    "    ax3 = fig.add_subplot(gs[1, 1])\n",
    "    ax3.hist(prop_green,cumulative = True, density = True, alpha=0.7,color=\"green\",edgecolor=\"black\",linewidth=2,bins = np.linspace(0,1,50), label = 'Bootstrap')\n",
    "    ax3.plot([green_lower,green_lower],[0,1],color='black',linewidth=2,linestyle='--',label='Lower/Upper')\n",
    "    ax3.plot([green_upper,green_upper],[0,1],color='black',linewidth=2,linestyle='--')\n",
    "    ax3.plot([green_prop,green_prop],[0,1],color='black',linewidth=3,label='Exp.')\n",
    "    ax3.set_title('Uncertainty in Proportion of Green Balls'); ax3.set_xlabel('Proportion of Green Balls'); ax3.set_ylabel('Cumulative Probability') \n",
    "    ax3.set_xlim([0,1]); ax3.set_ylim([0,1])         \n",
    "\n",
    "    ax3.plot(green_prop_values, cumul_prob, color = 'black', linewidth = 2, label = 'Analytical')\n",
    "    ax3.legend()\n",
    "    \n",
    "    ax4 = fig.add_subplot(gs[2, 1])\n",
    "    ax4.hist(prop_green,cumulative = False, density = True, alpha=0.7,color=\"green\",edgecolor=\"black\",linewidth=2, bins = np.linspace(0,1,50), label = 'Bootstrap Prop. Green')\n",
    "    ax4.hist(prop_red,cumulative = False, density = True, alpha=0.7,color=\"red\",edgecolor=\"black\",linewidth=2, bins = np.linspace(0,1,50), label = 'Bootstrap Prop. Red')\n",
    "        \n",
    "    ax4.set_title('Confidence Interval in Proportion of Red and Green Balls (Alpha = ' + str(alpha) + ')'); ax3.set_xlabel('Proportion of Green Balls')\n",
    "    ax4.set_xlabel('Proportion of Red and Green Balls'); ax4.set_ylabel('Frequency') \n",
    "    ax4.set_xlim([0,1])\n",
    "    \n",
    "    prop_values = np.linspace(0.0,1.0,100)\n",
    "    if a <= 30 and b <= 30:\n",
    "        green_density = t.pdf(prop_values,loc = green_prop, df = dof, scale = green_SE)\n",
    "    else:\n",
    "        green_density = norm.pdf(prop_values,loc = green_prop, scale = green_SE)\n",
    "    ax4.plot(prop_values, green_density, color = 'black', linewidth = 5,zorder=99)          \n",
    "    ax4.plot(prop_values, green_density, color = 'green', linewidth = 3, label = 'Analytical Prop. Green',zorder=100)\n",
    "    \n",
    "    if a <= 30 and b <= 30:\n",
    "        red_density = t.pdf(prop_values,loc = red_prop, df = dof, scale = red_SE)\n",
    "    else:\n",
    "        red_density = norm.pdf(prop_values,loc = red_prop, scale = red_SE)\n",
    "    ax4.plot(prop_values, red_density, color = 'black', linewidth = 5,zorder=99)           \n",
    "    ax4.plot(prop_values, red_density, color = 'red', linewidth = 3, label = 'Analytical Prop. Red',zorder=100)\n",
    "        \n",
    "    ax4.fill_between(prop_values, 0, green_density, where = prop_values <= green_lower, facecolor='green', interpolate=True, alpha = 0.9,zorder=101)\n",
    "    ax4.fill_between(prop_values, 0, green_density, where = prop_values >= green_upper, facecolor='green', interpolate=True, alpha = 0.9,zorder=101)\n",
    "    ax4.fill_between(prop_values, 0, red_density, where = prop_values <= red_lower, facecolor='darkred', interpolate=True, alpha = 0.9,zorder=101)\n",
    "    ax4.fill_between(prop_values, 0, red_density, where = prop_values >= red_upper, facecolor='darkred', interpolate=True, alpha = 0.9,zorder=101)\n",
    "    ax4.legend()\n",
    "\n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.5, top=3.0, wspace=0.2, hspace=0.3)\n",
    "    plt.show()\n",
    "\n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(s_make, {'a': a2, 'b': b2, 'c': c2, 'alpha': alpha})\n",
    "interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Bootstrap and Analytical Confidence Interval Demonstration for Sample Proportions\n",
    "\n",
    "* drawing red and green balls from a hat with replacement to access uncertainty in the proportion \n",
    "\n",
    "* run many bootstrap realizations and summarize the results and compare to the analytical sampling distribution for the proportion\n",
    "\n",
    "* interactive plot demonstration with ipywidget, matplotlib packages\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, The 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",
    "Let's simulate bootstrap, resampling with replacement from a hat with $n_{red}$ and $n_{green}$ balls\n",
    "\n",
    "* **$n_{red}$**: number of red balls in the sample (placed in the hat)\n",
    "\n",
    "* **$n_{green}$**: number of green balls in the sample (placed in the hat)\n",
    "\n",
    "* **$L$**: number of bootstrap realizations\n",
    "\n",
    "* **$\\alpha$**: alpha level for the confidence interval (significance level)\n",
    "\n",
    "and then compare the uncertainty in the proportion of balls to the analytical expression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence Interval Demonstration for Sample Proportions, Analytical and Bootstrap\n",
    "\n",
    "#### Michael Pyrcz, 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "90e585d00a8845abb1b53b70ba0b32ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='     Confidence Interval for Proportions, Analytical and Bootstrap Demonstration, M…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4bff03c928e542f9a5e7e891c9d382ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 432x288 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(uib2, interactive_plot)                           # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Confidence Interval Demonstration for Uncertainty in the Population Mean\n",
    "\n",
    "#### Michael Pyrcz, Professor, The 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",
    "Let's calculate the confidence interval, given:\n",
    "\n",
    "* **$n$**: number of samples\n",
    "\n",
    "* **$\\overline{x}$**: sample mean\n",
    "\n",
    "* **$s_x$**: sample standard deviation\n",
    "\n",
    "* **$\\alpha$**: alpha level for the confidence interval (significance level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interactive calculation of the sample set (control of source parametric distribution and number of samples)\n",
    "l2 = widgets.Text(value='     Confidence Interval for Mean, Analytical Explained, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "n = widgets.IntSlider(min=1, max = 20, value = 10, step = 1, description = '$n$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "n.style.handle_color = 'red'\n",
    "\n",
    "smean = widgets.FloatSlider(min=1.0, max = 4.0, value = 2.0, step = 0.11, description = '$\\overline{x}$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "smean.style.handle_color = 'green'\n",
    "\n",
    "sstdev = widgets.FloatSlider(min=0.1, max = 2.0, value = 1.0, step = 0.1, description = '$s_x$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "sstdev.style.handle_color = 'gray'\n",
    "\n",
    "alpha = widgets.FloatSlider(min=0.01, max = 0.40, value = 0.05, step = 0.01, description = r'$\\alpha$',orientation='horizontal',layout=Layout(width='400px', height='20px'),continuous_update=False)\n",
    "alpha.style.handle_color = 'gray'\n",
    "\n",
    "uib = widgets.HBox([n,smean,sstdev,alpha],)                            # basic widget formatting           \n",
    "uib3 = widgets.VBox([l2,uib],)\n",
    "\n",
    "def ci_make(n,smean,sstdev,alpha):                                      # function to take parameters, make sample and plot\n",
    "    dof = n-1 \n",
    "    tscore = np.round(-1*stats.t.ppf(alpha/2.0, df=dof),2)\n",
    "    se = np.round(sstdev/(np.sqrt(n)),2)\n",
    "    smean = np.round(smean,2)\n",
    "    sstdev = np.round(sstdev,2)\n",
    "    alpha = np.round(alpha,2)\n",
    "    xval = np.linspace(0,10,1000)\n",
    "    pdf = stats.t.pdf(xval,df = dof, loc = smean, scale = se)\n",
    "    lower_CI = np.round(smean - tscore*se,2)\n",
    "    upper_CI = np.round(smean + tscore*se,2)\n",
    "    \n",
    "    plt.plot(xval,pdf,color='black',lw=2,zorder=100)\n",
    "    plt.vlines(smean,0,stats.t.pdf(smean,df = dof, loc = smean, scale = sstdev/np.sqrt(n)),color='red',ls='--',zorder=10)\n",
    "    plt.vlines(lower_CI,0,stats.t.pdf(lower_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n)),color='darkred',zorder=10)\n",
    "    plt.vlines(upper_CI,0,stats.t.pdf(lower_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n)),color='darkred',zorder=10)\n",
    "    plt.fill_between(xval,pdf,where= xval< lower_CI, color='red',alpha=0.7)\n",
    "    plt.fill_between(xval,pdf,where= xval> upper_CI, color='red',alpha=0.7)\n",
    "    add_grid()\n",
    "    \n",
    "    vset = +0.4\n",
    "    plt.annotate(r'$CI_{\\mu} \\rightarrow \\overline{x} \\pm t_{\\left(\\frac{\\alpha}{2},dof\\right)} \\times \\frac{s_{x}}{\\sqrt{n}}$',[6.3,1.2+vset],size=12)\n",
    "    # plt.annotate(r'$CI_{\\mu} = $' + str(smean) + r'$ \\pm $' + str(tscore) + r' $\\times$ ',[6.3,1.1],size=12)\n",
    "    plt.annotate(r'$\\overline{x} = $' + str(smean),[6.4,1.35+vset],color='red'); plt.plot([6.8,7.0],[1.33+vset,1.27+vset],color='red')\n",
    "    plt.annotate(r'$\\alpha = $' + str(alpha),[7.4,1.355+vset],color='red'); plt.plot([7.8,7.75],[1.33+vset,1.26+vset],color='red')\n",
    "    plt.annotate(r'$s_{x} = $' + str(sstdev),[8.6,1.355+vset],color='red'); plt.plot([8.85,8.75],[1.33+vset,1.27+vset],color='red')\n",
    "    plt.annotate(r'$dof = $' + str(dof),[7.4,1.05+vset],color='red'); plt.plot([7.9,7.8],[1.15+vset,1.10+vset],color='red')\n",
    "    plt.annotate(r'$n = $' + str(n),[8.8,1.05+vset],color='red'); plt.plot([8.8,8.9],[1.15+vset,1.10+vset],color='red')\n",
    "    \n",
    "    vset = +0.2\n",
    "    plt.plot([8.7,9.1],[1.115+vset,1.115+vset],lw=0.5,color='black')\n",
    "    plt.annotate(r'$CI_{\\mu} \\rightarrow $' + str(smean) + r'$ \\pm $' + str(tscore) + r' $\\times$ ',[6.3,1.1+vset],size=12)\n",
    "    plt.annotate(sstdev,[8.70,1.13+vset]); plt.annotate(np.round(np.sqrt(n),2),[8.68,1.07+vset])\n",
    "    #plt.annotate(r'$CI_{\\mu} = $' + str(smean) + '$ \\pm t_{\\left(\\frac{alpha}{2},dof\\right)} \\times \\frac{\\sigma_{s}}{\\sqrt{n}}$',[6.3,1.1],size=12)\n",
    "    \n",
    "    vset = +0.1\n",
    "    plt.annotate(r'$CI_{\\mu} \\rightarrow $' + str(smean) + r'$ \\pm $' + str(tscore) + r' $\\times$ ' + str(se),[6.3,0.95+vset],size=12)\n",
    "    plt.annotate('statistic',[7.1,0.65+vset],rotation=270.0,ha='left',color='red'); plt.plot([7.23,7.23],[0.86+vset,0.93+vset],color='red')\n",
    "    plt.annotate('-score',[7.8,0.7+vset],rotation=270.0,ha='left',color='red'); plt.plot([7.93,7.93],[0.86+vset,0.93+vset],color='red')\n",
    "    plt.annotate('standard error',[8.8,0.48+vset],rotation=270.0,ha='left',color='red'); plt.plot([8.93,8.93],[0.86+vset,0.93+vset],color='red')\n",
    "    \n",
    "    plt.annotate(r'$CI_{\\mu} \\rightarrow [$' + str(lower_CI) + ' , ' + str(upper_CI) + r'$]$',[6.3,0.5+vset],size=12)\n",
    "    \n",
    "    plt.annotate('lower interval = ' + str(lower_CI),[lower_CI-0.1,stats.t.pdf(lower_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+0.05],rotation=90.0)\n",
    "    plt.annotate('upper interval = ' + str(upper_CI),[upper_CI-0.1,stats.t.pdf(upper_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+0.05],rotation=90.0)\n",
    "    plt.annotate('$\\overline{x} = $' + str(smean),[smean-0.1,stats.t.pdf(smean,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+0.05],rotation=90.0)\n",
    "    \n",
    "    \n",
    "    #plt.annotate(r'$F^{-1}_x($' + str(np.round(P,2)) + '$)$ = ' + str(np.round(x,2)),xy=[x+0.003,0.08],rotation=270.0)\n",
    "    \n",
    "    plt.xlim([0,10]); plt.ylim([0,2.0]); plt.ylabel('Density'); plt.xlabel('Value'); plt.title('Analytical Confidence Interval for Population Mean')\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.1, top=1.2, wspace=0.2, hspace=0.3); plt.show()\n",
    "    \n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot3 = widgets.interactive_output(ci_make, {'n':n,'smean':smean,'sstdev':sstdev,'alpha':alpha})\n",
    "interactive_plot3.clear_output(wait = True)                # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Confidence Interval Demonstration for Uncertainty in the Population Mean\n",
    "\n",
    "#### Michael Pyrcz, 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d14da64f46b485a8e9cdc8daffbe2ef",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='     Confidence Interval for Mean, Analytical Explained, Michael Pyrcz, Professor, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f73af8c617146eabee8ea7d272742fb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 432x288 with 1 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(uib3, interactive_plot3)                           # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observations\n",
    "\n",
    "Some observations:\n",
    "\n",
    "* sampling distribution for proportions become discrete with too few samples, as only $n$ cases are possible\n",
    "\n",
    "* enough bootstrap realizations are required for stable statistics\n",
    "\n",
    "* the analytical sampling distribution for the uncertainty in the sample proportion matches the results from bootstrap\n",
    "\n",
    "\n",
    "#### 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. Professor, Cockrell School of Engineering and 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 (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": 2
}