Interactive_Bootstrap.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 Bootstrap Demonstration\n",
    "\n",
    "### Boostrap for Uncertainty in Sample Statistics 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",
    "#### 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.\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 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",
    "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"
   ]
  },
  {
   "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": 2,
   "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='                                      Boostrap 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 = 'Dataset: Mean / Mode',orientation='vertical',layout=Layout(width='170px', 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 = 'Dataset: St. Deviation',orientation='vertical',layout=Layout(width='130px', height='200px'))\n",
    "d.style.handle_color = 'green'\n",
    "b = widgets.FloatSlider(min = 0, max = 100.0, value = 0.5, description = 'Dataset: Minimum',orientation='vertical',layout=Layout(width='130px', height='200px'))\n",
    "b.style.handle_color = 'red'\n",
    "c = widgets.IntSlider(min = 0, max = 100, value = 100, description = 'Dataset: Maximum',orientation='vertical',layout=Layout(width='130px', height='200px'))\n",
    "c.style.handle_color = 'orange'\n",
    "n = widgets.IntSlider(min = 2, max = 1000, value = 100, description = 'Dataset: Number Samples',orientation='vertical',layout=Layout(width='180px', height='200px'))\n",
    "n.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([dist,a,d,b,c,n],)                      # basic widget formatting           \n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def f_make(dist,a, b, c, d, n):                           # function to take parameters, make sample and plot\n",
    "    dataset = make_data(dist,a, b, c, d, n)\n",
    "    \n",
    "    plt.subplot(111) \n",
    "    plt.hist(\n",
    "        dataset,\n",
    "        alpha=0.2,\n",
    "        color=\"red\",\n",
    "        edgecolor=\"black\",\n",
    "        bins=bins) \n",
    "    plt.xlim(0.0,100.0); plt.title('Synthetic Dataset'); plt.ylabel('Frequency'); plt.xlabel('Data Values')\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",
    "def make_data(dist,a, b, c, d, n):                        # function to check parameters and make sample   \n",
    "    if dist == 'Uniform':\n",
    "        if b >= c:\n",
    "            print('Invalid uniform distribution parameters')\n",
    "            return None\n",
    "        dataset = uniform.rvs(size=n, loc = b, scale = c, random_state = 73073).tolist()\n",
    "        return dataset\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 = triang.rvs(size=n, loc = b, c = (a-b)/interval, scale = interval, random_state = 73073).tolist()\n",
    "        return dataset\n",
    "    elif dist == 'Gaussian':\n",
    "        dataset = norm.rvs(size=n, loc = a, scale = d, random_state = 73073).tolist()\n",
    "        return dataset\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, 'n': n})\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": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd41971bb7294b1ba27547e38e42f77d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                      Boostrap Demonstration, Michael Pyrcz, Associ…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3a6295ccafcb4dca8913f9bdb6016157",
       "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": [
    "We now have a synthetic dataset to work with. Now we can:\n",
    "\n",
    "* assign the samples to a 1D ndarray \n",
    "\n",
    "* make a DataFrame with the samples\n",
    "\n",
    "* check the summary statistics\n",
    "\n",
    "This is our sample set that we will apply ot bootstrap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DataSet</th>\n",
       "      <td>27.0</td>\n",
       "      <td>57.671481</td>\n",
       "      <td>11.604119</td>\n",
       "      <td>22.400407</td>\n",
       "      <td>51.761746</td>\n",
       "      <td>58.371597</td>\n",
       "      <td>67.016427</td>\n",
       "      <td>76.082631</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         count       mean        std        min        25%        50%  \\\n",
       "DataSet   27.0  57.671481  11.604119  22.400407  51.761746  58.371597   \n",
       "\n",
       "               75%        max  \n",
       "DataSet  67.016427  76.082631  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = make_data(dist.value,a.value, b.value, c.value, d.value, n.value)\n",
    "df = pd.DataFrame({'DataSet':data})\n",
    "df.describe().transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform Bootstrap \n",
    "\n",
    "Now we take our synthetic dataset, sampled from the parametric distributioin above, and apply it to statistical bootstrap.\n",
    "\n",
    "* we calculate the sampling distributions / uncertainty in the mean and standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5wAAAGWCAYAAAAUtoyeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xkZX3v+8/XGa6KDobGcHXQoAlIRJngZHs5RHRLjAHdJxg40WAkQROMsmO2gvEo8WwSkuM126hBJYAXcCIqxMQoEtGYLeKMEG6NAQRhZIT2MgKjMgF++4+1JtY03dPVPb26uqo/79erXrXWsy71W09dnvrVetZTqSokSZIkSZpvDxt0AJIkSZKk0WTCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmntEQkuS7JEQvwOPsnuTfJsnna3+lJPjwf+5IkDQfbrMFLUkl+bh73d2uS58zX/mb52HN+PSV5ZpJvzHNIS4oJp2at/cD4cfsB/YMk/5Bkv3nY78r2w215n+vP6wdhH4/3siRfnqJ8wT9AZ1tXAFV1cFVd1uf+53xMVXVbVT2iqh6Y7bZJjkiyfi6P2+f+z2nr7ehJ5e9sy1/W1WNLGgzbrIeU22Zt/TiLuc3aN8mFSb6b5IdJrtnSTs2lTheztn3enOSe9nZtkj9P8qj52P8sX09bvVer6l+q6onzEcdSZcKpufr1qnoEsBdwJ/C/BhzPQ4zKh/Bko3pcC+jfgRO2zLT1eSxw88AiktQ126wBGdXjWiAfAm4HHgv8DPDbNK/fRW07nvO/rKrdgDHgd4DVwL8mefi8BaeBMOHUdqmqnwAfBw7aUpbkUUnOSzKR5FtJ3pjkYe2yh7Xz30pyV7vell+vvtTeb2x/if7lJD+X5IvtL3vfTfKxdj9b1v23dt3f3PJLY5LXJ/kO8LdJdk/y6TaWH7TT+/bEeln7C9oV7WNclOTRc62P9he6v25/Qb8nyVeTPL5n+cFJLkny/SR3JnlDT72cmuTmJN9LsmZLHD2/Yp6Y5Dbgn6epq8cn+ed2++8m+UiSFT2P/Z+/AKfp8rOmrf970nQ1WdUu+xCwP/D37b5f1x7PH0461quTvHCKOtjqV9e2jv+/JP/aPtbnkuwxxXYPBz4D7N0+7r1J9m4X7zhVrO12e7e/AE8kuSXJq2d4mv4eeHqS3dv5o4Crge9MiuflScbb181nkzy2Z9m7ktye5O4k65I8s2fZtHUrabBss7Zmm7Xo26xfAs6pqk1VdX9VXVlVn2mXzbVO/7itix8m+ViSnXuW/48kG5LckeTlk47315Jcmabduz3J6VPUYe9zTpKXtu+d7yX5k20c51aq6idV9TXgaJpE+3d6HmvKtjnJ+5K8dVLMFyX5o55j3/J6OjzJV5JsbI/33Ul2bJdN+17t2e8vtK+Tje3ze3TPsm2+p5asqvLmbVY34FbgOe30rsC5wHk9y88DLgJ2A1bSnFE6sV32cuAm4HHAI4BPAB9ql60ECljes6/zgT+h+XFkZ+AZPcsK+Lme+SOA+4G/AHYCdqH5oPq/2zh3A/4O+FTPNpcB3waeBDwcuBD48DTH/TLgyzPUxznA94HDgeXAR4AL2mW7ARuA17bHshvwtHbZKcDlwL5t7H8DnD+pXs5rY9xlmrr6OeC57fZjNI3RO6eJ83TgJ8DzgWXAnwOXT7VuO/9i4Ks9808GvgfsOEV9bBVbW8c3A09oY78MOHOaOj4CWD+pbNpY29fFOuBNwI40r6tvAs+bZv/nAP8TOAv4/bZsDXA88GXgZW3ZC2lep7/QPo9vBP53z35eQvPaWt4+n98Bdu6nbr1587awN2yztlUf52CbtVVsLK426/PAvwLHAftvK+5Z1OkVwN7Ao4Fx4JXtsqNozp5ueW19lJ7XbHush7TH8Ivtui/cxnN+EHAv8Kw2nrfTvN6fM82xngP8zynKzwM+1k5P2za3j3M7kHZ+d+DHwN5TvJ4Oozl7uryNfRw4ZYb36vp2eoc2hje0z+GzgXuAJ870nlrKt4EH4G34bu2b9l5gY/vhcQdwSLtsGXAfcFDP+q8ALmunLwX+oGfZE4H/6HnTT/7wPI8mOdh3ijim+kDYTPvFf5rYDwV+0DN/GT0NSfsBuRlYNsW2L6O/xvsDPcueD9zQTh8PXDlNXOPAkT3ze01RL4/rWf6Quppiny/sfTwe2nh/ftJx/3iqddv5ndoP0APb+bcC75nmcbeKra3jN/Ys/wPgn6bZ9gimbrynjBV4GnDbpPVPA/52mv2fQ5NwPgP4CvAomkZzF7ZOOD9D+4WznX8Y8CPgsdPs9wfAk/upW2/evC3sDdusqerDNmua2FhcbdbuwJnAdcADwFXAL21nnb6kZ/4vgfe102dPem09gUmv2Un7fifwjkmx9D7nb6In0aJJRDcz+4TzTOCSdnrathkIcBvwrHbZ7wH/PN1rZNJjnAJ8cob36paE85k0PzI/rGf5+cDpM72nlvLNLrWaqxdW1QqaD/VXAV9M8rPAHjS/+HyrZ91vAfu003tPsWw58JhpHud1NB8iV7TdFl4+zXpbTFTTZQqAJLsm+Zu2S8fdNL/2rcjWo9HdPimeHdrjmOz+dtlkO9A0tFv0ds38Ec2v4gD7Mf11go8FPtl2z9hI05g/wNb1cvuUW7aS7JnkgiTfbo/1w9Mcx3Rx7pxprruoqvtozgS+JE1Xs+Npri3p13R1Mtftt8T6WJruTBt76u4NTP96AqCqvkzz6+8bgU9X1Y8nrfJY4F09+/w+zetwH4Akr2279PywXf4otq7rvutW0oKwzfop26yZLYo2q6p+UFWnVtXB7TpXAZ9KkqnW77NOpzu2vXnoa6t3309L8oW2K/APgVdOse/e7bfaX1VtojnLPFv70LTBsI22uZrs7gKa5xrg/6E5u/gQSZ6Qprv6d9p6+rMpjmU6ewO3V9WDPWW9nxmw/a+fkWPCqe1SVQ9U1SdoGppnAN+lacge27Pa/jRdgKD5ZXnysvtpzjLVFPv/TlX9XlXtTfOr83uy7VH+Ju/jtTS/SD+tqh5J0+UCmg+oLXpHK9y/jf+7U+z7NmD/3g/6JLsCezLpg3katwPT9eO/HfjVqlrRc9u5qr7ds05NM73Fn7flv9ge60vY+jhnY6r9nwv8FnAk8KOq+soc9z3bx92W24FbJtXbblX1/D62/TDN6+O8afb7ikn73aWq/nea6zVfT9Nla/f2S+wPmXtdS1ogtlm2WfNswdqsqvouzZnaLd1h57tON/DQ11avjwIXA/tV1aOA902x796Yttpf+9r7mT5j2bLNI4DnAP/SFk3bNrfLzwd+o72u82k0Xc6n8l7gBpoz4I+kSfr7rac7gP3aHzK26P3M0BRMOLVd0jiGptvHeDXDiq8BzkiyW/um/yOaL/fQfBj89yQHtB8kf0bTN/9+YAJ4kOaahi37PzY/HTDhBzQfZluGLr+zd91p7EbTh39jmgEN3jzFOi9JclD7YfgW4OM19fDoX6W5LuPUJDunGTDgTGAt/TXenwZ+NskpSXZq6+dp7bL30dTZY9vjHmvrdToPqav2WO9tj3Uf4H/0EdN0HlK3bWP9IPA2ZvdL8Wwf92fS/zDoVwB3pxl0Y5cky5I8Kckv9bHtX9Fc6/KlKZa9DzgtycHwn4OKHNsu243mC+cEsDzJm4BH9hmvpAGyzbLNmmedtllJ/qJdvjzJbsDvAzdV1feY/zpdA7ys57U1+bW3G/D9qvpJksNpziBuy8eBFyR5RpoBed5Cn3lH+3o7DPgUzfvob9tF22qbqaoraerlA8Bnq2rjNA+xG3A3cG+Sn6ep117beq9+FdgEvC7JDmn+2/PXac6uahomnJqrv09yL80b9gzghKq6rl32hzRvxm/SXBf3UZprA2jvP0TzJf8WmsbwDwGq6kftvv617S6xmmaEtq+2j3Ux8JqquqXd1+nAue26L54mznfSXJ/3XZoBDv5pinU+RNPn/js0AyNMOWJc20Xn12j78rfHtzfw4rYrxzZV1T00Cc6vt491I/Ar7eJ3tcf3uST3tLE+bar9tPuaqq7+FHgqzdm2f6AZ3GKu/hx4Y7vvP+4pP49m0IBO/tS6qm6g+YL3zfax955h/Qdo6vNQmtfTd2kamhkb/6r6flVdOtVzV1WfpBnI44K2u821wK+2iz9Lcx3Jv9N8afsJM3QdkzRwtlm2WfNuAdqsXYFP0lx//E2as+1Ht/ua1zqtZvTbd9KMMHtTe9/rD4C3tM/3m2gS1G3t7zrgZJr30waaxHGm/yx9Xbv/79M8d+uA/9J2x52pbd7ifJqzoh/dxuP8MU3CfA/wfuBjk5afzjTv1araTPMc/CrN8/ce4Lfb14KmkT4+c6SRleQymhH+PjDoWIZBkt8GTqqqZww6FklaamyzZsc2S1ocPMMpqS9tF5s/oBmBUZKkRcs2S1o8TDglzSjJ82iui7iTbXdTkSRpoGyzpMXFLrWSJEmSpE54hlOSJEmS1Imh/iPyPfbYo1auXDnoMCRJi9C6deu+W1Vjg45jGNm+SpKmMpe2dagTzpUrV7J27dpBhyFJWoSS9PNfg5qC7askaSpzaVvtUitJkiRJ6oQJpyRJkiSpEyackiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpEyackiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpEyackiQNiST7JflCkvEk1yV5TVv+6CSXJLmxvd+9Z5vTktyU5BtJnje46CVJS5EJpyRJw+N+4LVV9QvAauDkJAcBpwKXVtWBwKXtPO2y44CDgaOA9yRZNpDIJUlLkgmnJElDoqo2VNXX2+l7gHFgH+AY4Nx2tXOBF7bTxwAXVNV9VXULcBNw+MJGLUlaykw4JUkaQklWAk8Bvgo8pqo2QJOUAnu2q+0D3N6z2fq2bKr9nZRkbZK1ExMTXYUtSVpiTDglSRoySR4BXAicUlV3b2vVKcpqqhWr6qyqWlVVq8bGxuYjTEmSWD7oACQNl2suv5zNGzfOuN6OK1ZwyOrVCxCRtLQk2YEm2fxIVX2iLb4zyV5VtSHJXsBdbfl6YL+ezfcF7li4aDUf/NyVNMxMOCXNyuaNGzmsj7Mf6+ySJ827JAE+CIxX1dt7Fl0MnACc2d5f1FP+0SRvB/YGDgSuWLiINR/83JU0zEw4JUkaHk8HXgpck+SqtuwNNInmmiQnArcBxwJU1XVJ1gDX04xwe3JVPbDwYUuSlioTTkmShkRVfZmpr8sEOHKabc4AzugsKEmStsFBgyRJkiRJnTDhlCRJkiR1wi61kiRJI+Dm8fG+1nM0W0kLyYRTkiRpBDy4aZOj2UpadOxSK0mSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqROcJZ5JlSa5M8ul2/tFJLklyY3u/e8+6pyW5Kck3kjyv69gkSZIkSd1ZiDOcrwHGe+ZPBS6tqgOBS9t5khwEHAccDBwFvCfJsgWIT5IkSZLUgU4TziT7Ar8GfKCn+Bjg3Hb6XOCFPeUXVNV9VXULcBNweJfxSZIkSZK60/UZzncCrwMe7Cl7TFVtAGjv92zL9wFu71lvfVu2lSQnJVmbZO3ExEQ3UUuSJEmStltnCWeSFwB3VdW6fjeZoqweUlB1VlWtqqpVY2Nj2xWjJEmSJKk7yzvc99OBo5M8H9gZeGSSDwN3JtmrqjYk2Qu4q11/PbBfz/b7And0GJ8kSZIkqUOdneGsqtOqat+qWkkzGNA/V9VLgIuBE9rVTgAuaqcvBo5LslOSA4ADgSu6ik+SJEmS1K0uz3BO50xgTZITgduAYwGq6roka4DrgfuBk6vqgQHEJ0mSJEmaBwuScFbVZcBl7fT3gCOnWe8M4IyFiEmSJEmS1K2F+B9OSZIkSdISZMIpSZIkSeqECackSUMkydlJ7kpybU/Zx5Jc1d5uTXJVW74yyY97lr1vcJFLkpaiQQwaJEmS5u4c4N3AeVsKquo3t0wneRvww571b66qQxcsOkmSephwSpI0RKrqS0lWTrUsSYAXA89eyJgkSZqOXWolSRodzwTurKobe8oOSHJlki8meeZ0GyY5KcnaJGsnJia6j1SStCSYcEqSNDqOB87vmd8A7F9VTwH+CPhokkdOtWFVnVVVq6pq1djY2AKEKklaCkw4JUkaAUmWA/8N+NiWsqq6r/3/a6pqHXAz8ITBRChJWopMOCVJGg3PAW6oqvVbCpKMJVnWTj8OOBD45oDikyQtQSackiQNkSTnA18BnphkfZIT20XHsXV3WoBnAVcn+Tfg48Arq+r7CxetJGmpc5RaSZKGSFUdP035y6YouxC4sOuYJEmajmc4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJGmIJDk7yV1Jru0pOz3Jt5Nc1d6e37PstCQ3JflGkucNJmpJ0lK1fNABSJKkWTkHeDdw3qTyd1TVW3sLkhwEHAccDOwNfD7JE6rqgYUIVIvTzePjfa+744oVHLJ6dYfRSBp1JpySJA2RqvpSkpV9rn4McEFV3QfckuQm4HDgKx2FpyHw4KZNHDY21te66yYmOo5G0qizS60kSaPhVUmubrvc7t6W7QPc3rPO+rbsIZKclGRtkrUTJhmSpHliwilJ0vB7L/B44FBgA/C2tjxTrFtT7aCqzqqqVVW1aqzPs1+SJM3EhFOSpCFXVXdW1QNV9SDwfppus9Cc0dyvZ9V9gTsWOj5J0tJlwilJ0pBLslfP7IuALSPYXgwcl2SnJAcABwJXLHR8kqSly0GDJEkaIknOB44A9kiyHngzcESSQ2m6y94KvAKgqq5Lsga4HrgfONkRaiVJC8mEU5KkIVJVx09R/MFtrH8GcEZ3EUmSND271EqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkT/g+nJEnSPLrm8svZvHHjjOvtuGIFh6xevQARSdLgmHBKkiTNo80bN3LY2NiM662bmFiAaCRpsOxSK0mSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOtFZwplk5yRXJPm3JNcl+dO2/NFJLklyY3u/e882pyW5Kck3kjyvq9gkSZIkSd3r8gznfcCzq+rJwKHAUUlWA6cCl1bVgcCl7TxJDgKOAw4GjgLek2RZh/FJkiRJkjrUWcJZjXvb2R3aWwHHAOe25ecCL2ynjwEuqKr7quoW4Cbg8K7ikyRJkiR1q9NrOJMsS3IVcBdwSVV9FXhMVW0AaO/3bFffB7i9Z/P1bdnkfZ6UZG2StRP+f5UkSZIkLVqdJpxV9UBVHQrsCxye5EnbWD1T7WKKfZ5VVauqatVYH3+qLEmSJEkajAUZpbaqNgKX0VybeWeSvQDa+7va1dYD+/Vsti9wx0LEJ0mSJEmaf12OUjuWZEU7vQvwHOAG4GLghHa1E4CL2umLgeOS7JTkAOBA4Iqu4pMkSZIkdWt5h/veCzi3HWn2YcCaqvp0kq8Aa5KcCNwGHAtQVdclWQNcD9wPnFxVD3QYnyRJkiSpQ50lnFV1NfCUKcq/Bxw5zTZnAGd0FZMkSZIkaeEsyDWckiRJkqSlx4RTkiRJktQJE05JkiRJUidMOCVJkiRJnTDhlCRJkiR1woRTkiRJktQJE05JkiRJUidMOCVJkiRJnTDhlCRJkiR1woRTkiRJktQJE05JkoZIkrOT3JXk2p6y/z/JDUmuTvLJJCva8pVJfpzkqvb2vsFFLklaipYPOgBJkjQr5wDvBs7rKbsEOK2q7k/yF8BpwOvbZTdX1aELG6L6cfP4eF/r3TI+zmFjYx1HI0ndMOGUJGmIVNWXkqycVPa5ntnLgd9YyJg0Nw9u2tRXInnj2rULEI0kdcMutZIkjZaXA5/pmT8gyZVJvpjkmdNtlOSkJGuTrJ2YmOg+SknSkmDCKUnSiEjyJ8D9wEfaog3A/lX1FOCPgI8meeRU21bVWVW1qqpWjdl9U5I0T0w4JUkaAUlOAF4A/FZVFUBV3VdV32un1wE3A08YXJSSpKXGhFOSpCGX5CiaQYKOrqof9ZSPJVnWTj8OOBD45mCilCQtRQ4aJEnSEElyPnAEsEeS9cCbaUal3Qm4JAnA5VX1SuBZwFuS3A88ALyyqr4/kMAlSUuSCackSUOkqo6foviD06x7IXBhtxFJkjQ9u9RKkiRJkjphwilJkiRJ6oQJpyRJkiSpEyackiRJkqRO9JVwJnlS14FIkrSU2LZKkpaCfs9wvi/JFUn+IMmKTiOSJGlpsG2VJI28vhLOqnoG8FvAfsDaJB9N8txOI5MkaYTZtkqSloK+r+GsqhuBNwKvB/4v4K+S3JDkv3UVnCRJo8y2VZI06vq9hvMXk7wDGAeeDfx6Vf1CO/2ODuOTJGkk2bZKkpaC5X2u927g/cAbqurHWwqr6o4kb+wkMkmSRpttqyRp5PWbcD4f+HFVPQCQ5GHAzlX1o6r6UGfRSZI0umxbJUkjr99rOD8P7NIzv2tbJkmS5sa2VZI08vpNOHeuqnu3zLTTu3YTkiRJS4JtqyRp5PWbcG5K8tQtM0kOA368jfUlSdK22bZKkkZev9dwngL8XZI72vm9gN/sJiRJkpYE21ZJ0sjrK+Gsqq8l+XngiUCAG6rqPzqNTJKkEWbbKklaCvo9wwnwS8DKdpunJKGqzuskKkmSlgbbVknSSOsr4UzyIeDxwFXAA21xATaKkiTNgW2rJGkp6PcM5yrgoKqqLoORJGkJsW2VJI28fkepvRb42S4DkSRpibFtlSSNvH7PcO4BXJ/kCuC+LYVVdXQnUUkaejePj/e13o4rVnDI6tUdRyMtSratkqSR12/CeXqXQUgavGsuv5zNGzfOuN4t4+McNjY243oPbtrU13rrJib6ik8aQacPOgBJkrrW79+ifDHJY4EDq+rzSXYFlnUbmqSFtHnjxr4SxBvXrl2AaKTRZ9sqSVoK+rqGM8nvAR8H/qYt2gf4VFdBSZI06mxbJUlLQb+DBp0MPB24G6CqbgT27CooSZKWANtWSdLI6zfhvK+qNm+ZSbKc5r/CJEnS3Ni2SpJGXr8J5xeTvAHYJclzgb8D/r67sCRJGnm2rZKkkddvwnkqMAFcA7wC+EfgjV0FJUnSEmDbKkkaef2OUvsg8P72JkmStpNtqyRpKegr4UxyC1NcV1JVj5v3iCRJWgJsWyVJS0FfCSewqmd6Z+BY4NHzH44kSUuGbaskaeT1dQ1nVX2v5/btqnon8OyOY5MkaWTZtkqSloJ+u9Q+tWf2YTS/yu7WSUSSJC0Bc21bk5wNvAC4q6qe1JY9GvgYsBK4FXhxVf2gXXYacCLwAPDqqvrs/B2FJEnb1m+X2rf1TN9P25jNezSSJC0dc21bzwHeDZzXU3YqcGlVnZnk1Hb+9UkOAo4DDgb2Bj6f5AlV9cD2hy9J0sz6HaX2V7oORJKkpWSubWtVfSnJyknFxwBHtNPnApcBr2/LL6iq+4BbktwEHA58ZS6PLUnSbPXbpfaPtrW8qt4+P+FIkrQ0zHPb+piq2tButyHJnm35PsDlPeutb8umiuck4CSA/ffffxYPLUnS9PoaNIjmupLfp2mk9gFeCRxEc62J13JKkjR7C9G2Zoqyh/wVC0BVnVVVq6pq1djY2Dw9vCRpqev3Gs49gKdW1T0ASU4H/q6qfrerwCRJGnHz2bbemWSv9uzmXsBdbfl6YL+e9fYF7tiOmCVJmpV+z3DuD2zumd9MMxKeJEmam/lsWy8GTminTwAu6ik/LslOSQ4ADgSumONjSJI0a/2e4fwQcEWST9J0xXkRW4+OJ0mSZmdObWuS82kGCNojyXrgzcCZwJokJwK3AccCVNV1SdYA19OMhHuyI9RKkhZSv6PUnpHkM8Az26LfqaoruwtLkqTRNte2taqOn2bRkdM9DnDG3KKUJGn79NulFmBX4O6qehewvu2aI0mS5s62VZI00vpKOJO8meb/vE5ri3YAPtxVUJIkjTrbVknSUtDvGc4XAUcDmwCq6g78OxRJkraHbaskaeT1m3Burqqi/e+uJA/vLiRJkpYE21ZJ0sjrN+Fck+RvgBVJfg/4PPD+7sKSJGnk2bZKkkZev6PUvjXJc4G7gScCb6qqSzqNTJKkEWbbKklaCmZMOJMsAz5bVc8BbAglSdpOtq2SpKVixi617R9E/yjJo2az4yT7JflCkvEk1yV5TVv+6CSXJLmxvd+9Z5vTktyU5BtJnjfro5EkaQjMtW2VJGnY9NWlFvgJcE2SS2hH0wOoqldvY5v7gddW1deT7Aasa7d/GXBpVZ2Z5FTgVOD1SQ4CjgMOBvYGPp/kCW2jLEnSqJlL2ypJ0lDpN+H8h/bWt6raAGxop+9JMg7sAxwDHNGudi5wGc3/kB0DXFBV9wG3JLkJOBz4ymweV5KkITHrtlWSpGGzzYQzyf5VdVtVnbs9D5JkJfAU4KvAY9pklKrakGTPdrV9gMt7Nlvflk3e10nASQD777//9oQlSdKCm6+2VZKkYTDTNZyf2jKR5MK5PECSRwAXAqdU1d3bWnWKsnpIQdVZVbWqqlaNjY3NJSRJkgZpu9tWSZKGxUwJZ0x85E4AABMDSURBVG8S+LjZ7jzJDjTJ5keq6hNt8Z1J9mqX7wXc1ZavB/br2Xxf4I7ZPqYkSYvcdrWtkiQNk5kSzppmekZJAnwQGK+qt/csuhg4oZ0+Abiop/y4JDslOQA4ELhiNo8pSdIQmHPbKknSsJlp0KAnJ7mb5tfYXdpp2vmqqkduY9unAy+lGYHvqrbsDcCZwJokJwK3AcfS7Oy6JGuA62lGuD3ZEWolSSNoe9pWSZKGyjYTzqpaNtcdV9WXmfq6TIAjp9nmDOCMuT6mJEmL3fa0rZIkDZuZutRKkiRJkjQnJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTywcdgCRJ2n5Jngh8rKfoccCbgBXA7wETbfkbquofFzg8SdISZcIpSdIIqKpvAIcCJFkGfBv4JPA7wDuq6q0DDE+StETZpVaSpNFzJHBzVX1r0IFIkpY2E05JkkbPccD5PfOvSnJ1krOT7D7VBklOSrI2ydqJiYmpVpEkadZMOCVJGiFJdgSOBv6uLXov8Hia7rYbgLdNtV1VnVVVq6pq1djY2ILEKkkafSackiSNll8Fvl5VdwJU1Z1V9UBVPQi8Hzh8oNFJkpYUE05JkkbL8fR0p02yV8+yFwHXLnhEkqQly1FqJUkaEUl2BZ4LvKKn+C+THAoUcOukZZIkdcqEU5KkEVFVPwJ+ZlLZSwcUjiRJdqmVJEmSJHXDhFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHVi+aADkCRJ0uJ08/h4X+vtuGIFh6xe3XE0koaRCackSZKm9OCmTRw2NjbjeusmJhYgGknDyC61kiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpEw4aJEmS1IdrLr+czRs3zrjeLePjfQ20I0lLgQmnJElSHzZv3NhXInnj2rULEI0kDQe71EqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTJpySJEmSpE6YcEqSJEmSOmHCKUmSJEnqhAmnJEmSJKkTywcdgCRJmh9JbgXuAR4A7q+qVUkeDXwMWAncCry4qn4wqBglSUuLZzglSRotv1JVh1bVqnb+VODSqjoQuLSdlyRpQZhwSpI02o4Bzm2nzwVeOMBYJElLjAmnJEmjo4DPJVmX5KS27DFVtQGgvd9zqg2TnJRkbZK1ExMTCxSuJGnUeQ2nJEmj4+lVdUeSPYFLktzQ74ZVdRZwFsCqVauqqwAlSUuLZzglSRoRVXVHe38X8EngcODOJHsBtPd3DS5CSdJSY8IpSdIISPLwJLttmQb+K3AtcDFwQrvaCcBFg4lQkrQU2aVWkqTR8Bjgk0mgad8/WlX/lORrwJokJwK3AccOMEZJ0hLTWcKZ5GzgBcBdVfWktmza/wJLchpwIs1/h726qj7bVWySJI2aqvom8OQpyr8HHLnwEUmS1G2X2nOAoyaVTflfYEkOAo4DDm63eU+SZR3GJkmSJEnqWGcJZ1V9Cfj+pOLp/gvsGOCCqrqvqm4BbqIZ6ECSJEmSNKQWetCg6f4LbB/g9p711rdlD+H/hEmSJEnScFgso9RmirIp/wOsqs6qqlVVtWpsbKzjsCRJkiRJc7XQCed0/wW2HtivZ719gTsWODZJkiRJ0jxa6IRzuv8Cuxg4LslOSQ4ADgSuWODYJEmSJEnzqMu/RTkfOALYI8l64M3AmUzxX2BVdV2SNcD1wP3AyVX1QFexSZIkSZK611nCWVXHT7Noyv8Cq6ozgDO6ikeSJEmStLAWy6BBkiRJkqQRY8IpSZIkSeqECackSZIkqRMmnJIkSZKkTphwSpIkSZI6YcIpSZIkSeqECackSZIkqRMmnJIkSZKkTphwSpIkSZI6YcIpSZIkSeqECackSZIkqRMmnJIkSZKkTphwSpIkSZI6YcIpSZIkSeqECackSZIkqRMmnJIkSZKkTphwSpIkSZI6YcIpSZIkSeqECackSZIkqRMmnJIkSZKkTphwSpIkSZI6YcIpSZIkSerE8kEHIEmSpOF28/h4X+vtuGIFh6xe3XE0khYTE05JkiRtlwc3beKwsbEZ11s3MbEA0UhaTOxSK0nSCEiyX5IvJBlPcl2S17Tlpyf5dpKr2tvzBx2rJGnp8AynJEmj4X7gtVX19SS7AeuSXNIue0dVvXWAsUmSligTTkmSRkBVbQA2tNP3JBkH9hlsVJKkpc4utZIkjZgkK4GnAF9ti16V5OokZyfZfZptTkqyNsnaCa+zkyTNExNOSZJGSJJHABcCp1TV3cB7gccDh9KcAX3bVNtV1VlVtaqqVo31MfiLJEn9MOGUJGlEJNmBJtn8SFV9AqCq7qyqB6rqQeD9wOGDjFGStLSYcEqSNAKSBPggMF5Vb+8p36tntRcB1y50bJKkpctBgyRJGg1PB14KXJPkqrbsDcDxSQ4FCrgVeMVgwpMkLUUmnJIkjYCq+jKQKRb940LHIknSFnaplSRJkiR1woRTkiRJktQJE05JkiRJUidMOCVJkiRJnTDhlCRJkiR1woRTkiRJktQJ/xZFGnHXXH45mzdunHG9W8bHOWxsbAEi2trN4+N9rbfjihUcsnp1x9FIkiRpPplwSiNu88aNfSWSN65duwDRPNSDmzb1Fd+6iYkFiEaSJEnzyS61kiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpEyackiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpE8sHHYAkSdIgXXP55WzeuHHG9W4ZH+ewsbEFiGh03Tw+3td6O65YwSGrV3ccjaSFYMIpSZKWtM0bN/aVSN64du0CRDPaHty0qa+6XjcxsQDRSFoIdqmVJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ0w4JUmSJEmdMOGUJEmSJHXChFOSJEmS1AkTTkmSJElSJ5YPOgBJkiSp183j432tt+OKFRyyenXH0UjaHiackiRJWlQe3LSJw8bGZlxv3cTEAkQjaXvYpVaSJEmS1AnPcEqSpJF0zeWXs3njxhnXu2V8vK+zaZKk2Vt0CWeSo4B3AcuAD1TVmQMOSZKkobZU29bNGzf2lUjeuHbtAkQjSUvToko4kywD/hp4LrAe+FqSi6vq+sFGJi0+/nIvqR+2rRplDi4kLX6LKuEEDgduqqpvAiS5ADgGsFGUJllqv9z7pUKaM9tWjSwHF5IWv1TVoGP4T0l+Aziqqn63nX8p8LSqelXPOicBJ7WzTwKuXfBAR8MewHcHHcQQst7mzrqbO+tubp5YVbsNOohB66dtbctHrX0dhfeNx7A4eAyLg8ewOMy6bV1sZzgzRdlWGXFVnQWcBZBkbVWtWojARo11NzfW29xZd3Nn3c1NktE4vb/9ZmxbYfTaV49hcfAYFgePYXEYlWOY7TaL7W9R1gP79czvC9wxoFgkSRoFtq2SpIFZbAnn14ADkxyQZEfgOODiAcckSdIws22VJA3MoupSW1X3J3kV8FmaodvPrqrrtrHJWQsT2Uiy7ubGeps7627urLu5sd6YU9sKo1F3HsPi4DEsDh7D4rAkj2FRDRokSZIkSRodi61LrSRJkiRpRJhwSpIkSZI6MbQJZ5KjknwjyU1JTh10PItVkv2SfCHJeJLrkrymLX90kkuS3Nje7z7oWBejJMuSXJnk0+289daHJCuSfDzJDe1r75etu/4k+e/te/XaJOcn2dm6m1qSs5PcleTanrJp6yrJaW2b8Y0kzxtM1IvbKLStSW5Nck2Sq4bpr3Fm+3pejKY5htOTfLt9Pq5K8vxBxrgto/CdaRvHMDTPA0Db9l2R5N/a4/jTtnyYnovpjmHYnovt/i48lAlnkmXAXwO/ChwEHJ/koMFGtWjdD7y2qn4BWA2c3NbVqcClVXUgcGk7r4d6DTDeM2+99eddwD9V1c8DT6apQ+tuBkn2AV4NrKqqJ9EM8HIc1t10zgGOmlQ2ZV21n3vHAQe327ynbUvUGrG29Veq6tAh+7+7c+jz9byIncNDjwHgHe3zcWhV/eMCxzQbo/CdabpjgOF5HgDuA55dVU8GDgWOSrKa4XoupjsGGK7nYru/Cw9lwgkcDtxUVd+sqs3ABcAxA45pUaqqDVX19Xb6HpoXzD409XVuu9q5wAsHE+HilWRf4NeAD/QUW28zSPJI4FnABwGqanNVbcS669dyYJcky4Fdaf4v0bqbQlV9Cfj+pOLp6uoY4IKquq+qbgFuomlL9FO2rQM0y9fzojTNMQyNUfjOtI1jGCrVuLed3aG9FcP1XEx3DENjvr4LD2vCuQ9we8/8eobwzbTQkqwEngJ8FXhMVW2A5sMJ2HNwkS1a7wReBzzYU2a9zexxwATwt20XjA8keTjW3Yyq6tvAW4HbgA3AD6vqc1h3szFdXdluzGxU6qiAzyVZl+SkQQeznUblvf+qJFe3XW4XbRfIXqPwnWnSMcCQPQ9tV86rgLuAS6pq6J6LaY4Bhue5mJfvwsOacGaKsqH6xWChJXkEcCFwSlXdPeh4FrskLwDuqqp1g45lCC0Hngq8t6qeAmxicXd5WTTaRucY4ABgb+DhSV4y2KhGhu3GzEaljp5eVU+l6Rp8cpJnDTqgJe69wONpuhRuAN422HBmNgrfmaY4hqF7Hqrqgao6FNgXODzJkwYd02xNcwxD8VzM53fhYU041wP79czvS9PtTFNIsgPNh85HquoTbfGdSfZql+9F88uLfurpwNFJbqXpVvbsJB/GeuvHemB9z694H6dJQK27mT0HuKWqJqrqP4BPAP8F6242pqsr242ZjUQdVdUd7f1dwCcZ7q7TQ//er6o72y/dDwLvZ5E/H6PwnWmqYxi256FXe1nOZTTXBw/Vc7FF7zEM0XMxb9+FhzXh/BpwYJIDkuxIMxDExQOOaVFKEppr6car6u09iy4GTminTwAuWujYFrOqOq2q9q2qlTSvr3+uqpdgvc2oqr4D3J7kiW3RkcD1WHf9uA1YnWTX9r17JM31N9Zd/6arq4uB45LslOQA4EDgigHEt5gNfdua5OFJdtsyDfxX4Nptb7WoDf17f8sX09aLWMTPxyh8Z5ruGIbpeQBIMpZkRTu9C80PsjcwXM/FlMcwLM/FfH4XTtUw9paBdgjhd9KM4nh2VZ0x4JAWpSTPAP4FuIaf9r9+A01//jXA/jRfco+tqqG90L9LSY4A/riqXpDkZ7DeZpTkUJoLzHcEvgn8Ds0PXNbdDNph03+TZqTBK4HfBR6BdfcQSc4HjgD2AO4E3gx8imnqKsmfAC+nqdtTquozAwh7URv2tjXJ42jOakLTvf+jw3IMs309L0bTHMMRNF0HC7gVeMWW678Wm1H4zrSNYzieIXkeAJL8Is2ANMtovz9U1VuG6XvYNo7hQwzRcwHb/114aBNOSZIkSdLiNqxdaiVJkiRJi5wJpyRJkiSpEyackiRJkqROmHBKkiRJkjphwilJkiRJ6oQJp7RIJLksyfMmlZ2S5D3b2Obe7iOTJGk42bZKg2fCKS0e59P8sW6v49pySZI0e7at0oCZcEqLx8eBFyTZCSDJSmBv4Koklyb5epJrkhwzecMkRyT5dM/8u5O8rJ0+LMkXk6xL8tkke7Xlr05yfZKrk1zQ/eFJkrTgbFulAVs+6AAkNarqe0muAI4CLqL5BfZjwI+BF1XV3Un2AC5PcnFV1Uz7TLID8L+AY6pqIslvAmcALwdOBQ6oqvuSrOjosCRJGhjbVmnwTDilxWVL158tjeLLgQB/luRZwIPAPsBjgO/0sb8nAk8CLkkCsAzY0C67GvhIkk8Bn5rHY5AkaTGxbZUGyIRTWlw+Bbw9yVOBXarq6233nTHgsKr6jyS3AjtP2u5+tu4iv2V5gOuq6peneKxfA54FHA38v0kOrqr75+9QJElaFGxbpQHyGk5pEamqe4HLgLP56YAGjwLuahvEXwEeO8Wm3wIOSrJTkkcBR7bl3wDGkvwyNN2Akhyc5GHAflX1BeB1wArgEV0dlyRJg2LbKg2WZzilxed84BP8dFS9jwB/n2QtcBVww+QNqur2JGtouvLcCFzZlm9O8hvAX7WN5XLgncC/Ax9uywK8o6o2dntYkiQNjG2rNCDp49poSZIkSZJmzS61kiRJkqROmHBKkiRJkjphwilJkiRJ6oQJpyRJkiSpEyackiRJkqROmHBKkiRJkjphwilJkiRJ6sT/ARReFU1DVOtLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Summary Statistics for Bootstrap for Uncertainty in the Mean:\n",
      "DescribeResult(nobs=1000, minmax=(48.84610178059382, 64.34314390256239), mean=57.68535131216555, variance=4.8788353846381165, skewness=-0.09255827111358368, kurtosis=-0.03699400493390792)\n",
      "P10: 54.831, P50: 57.732, P90: 60.489\n",
      "\n",
      "Summary Statistics for Bootstrap for Uncertainty in the Standard Deviation:\n",
      "DescribeResult(nobs=1000, minmax=(5.769922838747761, 17.623458895883527), mean=11.018357175839355, variance=3.92583002943624, skewness=0.12807364716183928, kurtosis=-0.29031128339252144)\n",
      "P10: 8.412, P50: 11.051, P90: 13.481\n"
     ]
    }
   ],
   "source": [
    "L = 1000                                                  # set the number of realizations\n",
    "mean = np.zeros(L); stdev = np.zeros(L)                   # declare arrays to hold the realizations of the statistics\n",
    "for l in range(0, L):                                     # loop over realizations\n",
    "    samples = random.choices(df['DataSet'].values, weights=None, cum_weights=None, k=len(df))\n",
    "    mean[l] = np.average(samples)\n",
    "    stdev[l] = np.std(samples)\n",
    "    \n",
    "plt.subplot(121)                                          # plot the distribution for uncertainty in the mean     \n",
    "plt.hist(mean,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=np.linspace(0,100,40)); plt.xlim(0,100); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Mean')\n",
    "\n",
    "plt.subplot(122)                                          # plot the distribution for uncertainty in the standard deviation\n",
    "plt.hist(stdev,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=np.linspace(0,40,40)); plt.xlim(0,40); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Standard Deviation')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "plt.show()   \n",
    "\n",
    "# provide summary statistics, P10, P50 and P90\n",
    "print('Summary Statistics for Bootstrap for Uncertainty in the Mean:') \n",
    "print(stats.describe(mean))\n",
    "print('P10: ' + str(round(np.percentile(mean,10),3)) + ', P50: ' + str(round(np.percentile(mean,50),3)) + ', P90: ' + str(round(np.percentile(mean,90),3))) \n",
    "\n",
    "print('\\nSummary Statistics for Bootstrap for Uncertainty in the Standard Deviation:')\n",
    "print(stats.describe(stdev))\n",
    "print('P10: ' + str(round(np.percentile(stdev,10),3)) + ', P50: ' + str(round(np.percentile(stdev,50),3)) + ', P90: ' + str(round(np.percentile(stdev,90),3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Change the Number of Data\n",
    "\n",
    "Let's change the number of data drawn to observe the change in uncertainty\n",
    "\n",
    "* we will assume the same dataset and not recalculate it each time\n",
    "\n",
    "* we will just sample with replacement with the new number of samples for each bootstrap realization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the synthetic dataset\n",
    "bins = np.linspace(0,1000,1000)\n",
    "\n",
    "l = widgets.Text(value='                           Boostrap Demonstration with Modified Number of Data, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "n = widgets.IntSlider(min = 2, max = 1000, value = 100, description = 'New Number Samples',orientation='horizontal',layout=Layout(width='800px', height='20px'))\n",
    "n.style.handle_color = 'gray'\n",
    "\n",
    "ui3 = widgets.VBox([l,n],)\n",
    "\n",
    "def f_rerun(n):\n",
    "    L = 1000                                 # set the number of realizations\n",
    "    mean2 = np.zeros(L); stdev2 = np.zeros(L)    # declare arrays to hold the realizations of the statistics\n",
    "    for l in range(0, L):                      # loop over realizations\n",
    "        samples = random.choices(df['DataSet'].values, weights=None, cum_weights=None, k=n)\n",
    "        mean2[l] = np.average(samples)\n",
    "        stdev2[l] = np.std(samples)\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.hist(mean,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=np.linspace(0,100,40),label='Original'); plt.xlim(0,100); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Mean')\n",
    "    plt.hist(mean2,alpha=0.2,color=\"blue\",edgecolor=\"black\",bins=np.linspace(0,100,40),label='New'); plt.xlim(0,100); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Mean')\n",
    "    plt.legend()\n",
    "\n",
    "    plt.subplot(122)\n",
    "    plt.hist(stdev,alpha=0.2,color=\"red\",edgecolor=\"black\",bins=np.linspace(0,40,40),label='Original'); plt.xlim(0,40); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Standard Deviation')\n",
    "    plt.hist(stdev2,alpha=0.2,color=\"blue\",edgecolor=\"black\",bins=np.linspace(0,40,40),label='New'); plt.xlim(0,40); plt.ylabel('Frequency'); plt.xlabel('Values'); plt.title('Bootstrap Uncertainty in the Standard Deviation')\n",
    "    plt.legend()\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)  \n",
    "    plt.show()  \n",
    "\n",
    "interactive_plot3 = widgets.interactive_output(f_rerun, {'n': n})\n",
    "interactive_plot3.clear_output(wait = True)                # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display the GUI for Modifying the Number of Data\n",
    "\n",
    "We display the GUI now. Select the desired number of data\n",
    "\n",
    "* observed the change (new) vs. the original (original) bootstrap uncertainty as the number of data ranges from less than to greater than the original number of data  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "129642d3855f47789e2c0b042114872c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                           Boostrap Demonstration with Modified Number of Data, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "310a2b73e3154ef6874c199dc7a3b74c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui3, interactive_plot3)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observations\n",
    "\n",
    "Some observations form changing the number of data:\n",
    "\n",
    "* with less data the uncertainty increases\n",
    "\n",
    "* with more data the uncertianty decreases \n",
    "\n",
    "* there is a bias low in standard deviation with too few samples, as we fail to observe the dispersion well\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. 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
}