Interactive_Central_Limit_Theorem.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 of the Central Limit Theorem Demonstration\n",
    "\n",
    "### Summation or Averaging of Random Variables Tutorial\n",
    "\n",
    "* demonstrate the practical impact of summation or averaging on distributions\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",
    "#### Central Limit Theorem\n",
    "\n",
    "The summation of independent random variables will approach Gaussian distributed as the number of random variables become large\n",
    "\n",
    "* for any distribution shape for the individual random variables\n",
    "\n",
    "* same for averaging since it is just the application of a scalar to the summation ($\\frac{1}{m}$)\n",
    "\n",
    "We can state this as:\n",
    "\n",
    "\\begin{equation}\n",
    "Y \\sim N\\left[\\right]\n",
    "\\end{equation}\n",
    "\n",
    "where $Y$ is the summation of a large number of independent random variables, $X_i$, where $i = 1,\\ldots,m$.\n",
    "\n",
    "\\begin{equation}\n",
    "Y = \\sum_{i=1}^{m \\to \\infty} X_i\n",
    "\\end{equation}\n",
    "\n",
    "from statistical expectation we can calculate the central tedency as:\n",
    "\n",
    "\\begin{equation}\n",
    "E[Y] = E \\left[\\sum_{i=1}^{m \\to \\infty} X_i \\right] = \\sum_{i=1}^{m \\to \\infty} E \\left[X_i \\right]  \n",
    "\\end{equation}\n",
    "\n",
    "and under the assumption of independence the dispersion as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_Y^2 = \\sum_{i=1}^{m \\to \\infty} \\sigma_X^2  \n",
    "\\end{equation}\n",
    "\n",
    "therefore, we have our distribution as:\n",
    "\n",
    "\\begin{equation}\n",
    "Y \\sim N \\left[ \\sum_{i=1}^{m \\to \\infty} E[X_i], \\sum_{i=1}^{m \\to \\infty} \\sigma_X^2 \\right]\n",
    "\\end{equation}\n",
    "\n",
    "we can simplify by assuming the same central tendency and dispersion for all random variables, $X_i$, with central tendency and dispersion:\n",
    "\n",
    "\\begin{equation}\n",
    "E[X_i] = \\mu\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "Var[X_i] = \\sigma^2\n",
    "\\end{equation}\n",
    "\n",
    "now we have:\n",
    "\n",
    "\\begin{equation}\n",
    "Y \\sim N \\left[ m \\mu, m \\sigma_X^2 \\right]\n",
    "\\end{equation}\n",
    "\n",
    "for the case of the average instead of the summation of random variables, $X_i$:\n",
    "\n",
    "\\begin{equation}\n",
    "Y = \\frac{1}{m} \\sum_{i=1}^{m \\to \\infty} X_i\n",
    "\\end{equation}\n",
    "\n",
    "and all with the same central tendency and dispersion the distribution of $Y$ is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "Y \\sim N \\left[ \\mu, \\frac{\\sigma^2}{n} \\right]\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\n",
    "**Monte Carlo Simulation** is a method to sample from a single or set of random variables to assess the emergent distribution. \n",
    "\n",
    "* also known as Monte Carlo methods and Monte Carlo experiments\n",
    "\n",
    "* a method in statistics for generating draws from probability distributions empirically\n",
    "\n",
    "* powerful due to it's broad applicability\n",
    "\n",
    "We can apply Monte Carlo methods ot sample from the sum of our random variables.  We proceed as follows:\n",
    "\n",
    "* for all $X_i$ random variables draw a random value (random realization), $x_i$\n",
    "\n",
    "* calculate the average of the random realizations, $y_i = \\frac{1}{n} \\sum_{i=1}^{n} x_i$\n",
    "\n",
    "* repeat over a large number of samples, $n$, and observe the resulting distribution\n",
    "\n",
    "* compare the experimental CDF to the Gaussian distribution predicted from the Central Limit Theorem\n",
    "\n",
    "assess the uncertainty in a sample statistic by repeated random sampling with replacement.\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 demonstration of the central limit theorem   \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": 48,
   "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",
    "import numpy as np                                        # working with arrays\n",
    "import pandas as pd                                       # working with DataFrames\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 import stats                                   # statistical calculations\n",
    "import random                                             # random drawing / bootstrap realizations of the data\n",
    "import math                                               # square root operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Specify the Distributions of $X_i$ and the Number of Random Variables, $m$, and Realizations, $L$\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 random variables, $m$ \n",
    "\n",
    "* select the number of Monte Carlo realizations, $L$, to sample the distribution of $Y = \\sum_{i=1}^m X_i$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the synthetic dataset\n",
    "bins = np.linspace(0,100,50)\n",
    "\n",
    "# interactive calculation of the sample set (control of source parametric distribution and number of samples)\n",
    "l = widgets.Text(value='                              Central Limit Theorem Demonstration, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "dist = widgets.Dropdown(\n",
    "    options=['Triangular', 'Uniform', 'Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='Dataset Distribution:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "a = widgets.FloatSlider(min=0.0, max = 100.0, value = 0.5, description = '$X_i$: Mean / Mode',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "a.style.handle_color = 'blue'\n",
    "d = widgets.FloatSlider(min=0.01, max = 30.0, value = 5.0, step = 1.0, description = '$X_i$: St.Dev.',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "d.style.handle_color = 'green'\n",
    "b = widgets.FloatSlider(min = 0, max = 100.0, value = 0.5, description = '$X_i$: Min.',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "b.style.handle_color = 'red'\n",
    "c = widgets.IntSlider(min = 0, max = 100, value = 100, description = '$X_i$: Max.',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "c.style.handle_color = 'orange'\n",
    "m = widgets.IntSlider(min = 1, max = 20, value = 1, description = '$m$ ',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "m.style.handle_color = 'gray'\n",
    "L = widgets.IntSlider(min = 1, max = 1000, value = 100, description = '$L$ ',orientation='vertical',layout=Layout(width='110px', height='200px'))\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([dist,a,d,b,c,m,L],)                   # basic widget formatting   \n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def f_make(dist,a, b, c, d, m, L):                       # function to take parameters, make sample and plot\n",
    "    dataset, average, stdev = make_average_data(dist,a, b, c, d, m, L)\n",
    "    stdev = stdev / math.sqrt(m)\n",
    "\n",
    "    plt.subplot(221)\n",
    "    sample, null, null = make_average_data(dist,a, b, c, d, m=1, L=1000)\n",
    "    plt.hist(sample,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=bins) \n",
    "    plt.xlim(0.0,100.0); plt.title('Distribution, $X_i$'); plt.ylabel('Frequency'); plt.xlabel('Values')\n",
    "    \n",
    "    plt.subplot(222) \n",
    "    plt.hist(dataset,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=bins) \n",
    "    plt.xlim(0.0,100.0); plt.title('Distribution, $Y$'); plt.ylabel('Frequency'); plt.xlabel('Values')\n",
    "    \n",
    "    plt.subplot(223)    \n",
    "    plt.hist(dataset,cumulative = True, density = True, alpha=0.2,color=\"red\",edgecolor=\"black\", bins = bins, label = 'Bootstrap')\n",
    "    plt.xlim(0.0,100.0); plt.title('Comparison to Gaussian CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability') \n",
    "    \n",
    "    cumul_prob = np.linspace(0.0,1.0,100)\n",
    "    prop_values = norm.ppf(cumul_prob) \n",
    "    prop_values = prop_values * stdev + average           \n",
    "    plt.plot(prop_values, cumul_prob, color = 'black', linewidth = 2, dashes = (5,2,1,2), label = 'Analytical')\n",
    "\n",
    "    plt.legend()\n",
    "             \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=2.6, wspace=0.2, hspace=0.2)\n",
    "    plt.show()\n",
    "\n",
    "def make_average_data(dist,a, b, c, d, m, L):            # function to check parameters and make samples   \n",
    "    average = 0.0; stdev = 0.0\n",
    "    if dist == 'Uniform':\n",
    "        if b >= c:\n",
    "            print('Invalid uniform distribution parameters')\n",
    "            return None\n",
    "        dataset = np.mean(uniform.rvs(size=[m,L], loc = b, scale = c, random_state = 73073),axis = 0)\n",
    "        average = uniform.mean(loc = b, scale = c)\n",
    "        stdev = uniform.std(loc = b, scale = c)\n",
    "        return dataset, average, stdev\n",
    "    elif dist == 'Triangular':\n",
    "        interval = c - b\n",
    "        if b >= a or a >= c or interval <= 0:\n",
    "            print('Invalid triangular distribution parameters')\n",
    "            return None        \n",
    "        dataset = np.mean(triang.rvs(size=[m,L], loc = b, c = (a-b)/interval, scale = interval, random_state = 73073), axis = 0)\n",
    "        average = triang.mean(loc = b, c = (a-b)/interval, scale = interval)\n",
    "        stdev = triang.std(loc = b, c = (a-b)/interval, scale = interval)\n",
    "        return dataset, average, stdev\n",
    "    elif dist == 'Gaussian':\n",
    "        dataset = np.mean(norm.rvs(size=[m,L], loc = a, scale = d, random_state = 73073), axis = 0)\n",
    "        average = norm.mean(loc = a, scale = d)\n",
    "        stdev = norm.std(loc = a, scale = d)\n",
    "        return dataset, average, stdev\n",
    "\n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(f_make, {'dist': dist,'a': a, 'd': d, 'b': b, 'c': c, 'm': m, 'L':L})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display the GUI for Building the Synthetic Dataset\n",
    "\n",
    "We display the GUI now.  Select the desired parametric distribution and associated parameters.\n",
    "\n",
    "* if the parameters are invalid (e.g. traingular mode > max) an error message should display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bf528fbfe0dd4b6e989bf649662bdcc2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                              Central Limit Theorem Demonstration, Michael Pyrcz, A…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6a3c25b5ab5f483bad559c05134dbf90",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observations\n",
    "\n",
    "Some observations:\n",
    "\n",
    "* if $X_i$ are Gaussian, $Y$ is Gaussian distributed\n",
    "\n",
    "* if $X_i$ are uniform then convergence occurs over 5 or more RVs averaged\n",
    "\n",
    "* triangular distribution converges faster than uniform\n",
    "\n",
    "#### Comments\n",
    "\n",
    "This was a simple demonstration central limit theorem with 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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}