GeostatsPy_Confidence_Hypothesis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### GeostatsPy: Confidence Intervals and Hypothesis Testing for Subsurface Data Analytics in Python \n",
    "\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)\n",
    "\n",
    "##### Reporting Uncertainty and Significance\n",
    "\n",
    "With confidence intervals and hypothesis testing we have the opportunity to report uncertainty and to report significance in our statistics. These methods are based on standard methods with their associated limitations and assumptions.  See the lecture and workflow on Bootstrap https://git.io/fhgUW for a general, empirical approach to assess uncertianty in statistics.\n",
    "\n",
    "This is a tutorial / demonstration of **Confidence Intervals and Hypothesis Testing in Python** for Subsurface Modeling.  In Python, the SciPy package, specifically the Stats functions (https://docs.scipy.org/doc/scipy/reference/stats.html) provide excellent tools for efficient use of statistics. \n",
    "\n",
    "I have previously provided these examples worked out by-hand in Excel (https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg_R.xlsx) and also in R (https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg.R).  In all cases, I use the same dataset available as a comma delimited file (https://git.io/fxLAt).    \n",
    "\n",
    "This tutorial includes basic, typical confidence interval and hypothesis testing methods that would commonly be required for Engineers and Geoscientists including:\n",
    "\n",
    "1. Student-t confidence interval for the mean\n",
    "2. Student-t hypothesis test for difference in means (pooled variance)\n",
    "3. Student-t hypothesis test for difference in means (difference variances), Welch's t Test\n",
    "3. F-distribution hypothesis test for difference in variances \n",
    "\n",
    "##### Caveats\n",
    "\n",
    "I have not included all the details, specifically the test assumptions in this document.  These are included in the accompanying course notes, Lec08_hypothesis.pdf.\n",
    "\n",
    "#### Workflow Goal\n",
    "\n",
    "0. Introduction to Python in Jupyter including setting a working directory, loading data into a Pandas DataFrame.\n",
    "1. Learn the basics for working with confidence intervals and hypothesis testing in Python.  \n",
    "2. Demonstrate the efficiency of using Python and SciPy package for statistical analysis.\n",
    "3. Learn how to quantify uncertainty and significance in samples.\n",
    "\n",
    "#### Objective \n",
    "\n",
    "In the PGE 383: Stochastic Subsurface Modeling class I want to provide hands-on experience with building subsurface modeling workflows. Python provides an excellent vehicle to accomplish this. I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. \n",
    "\n",
    "The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    \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. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. 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",
    "You will need to copy the data file to your working directory.  They are available here:\n",
    "\n",
    "* Tabular data - PorositySample2Units.csv at https://git.io/fhrM8\n",
    "\n",
    "There are exampled below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code. \n",
    "\n",
    "#### Load the required libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geostatspy.GSLIB as GSLIB                       # GSLIB utilies, visualization and wrapper\n",
    "import geostatspy.geostats as geostats                 # GSLIB methods convert to Python     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os                                              # to set current working directory \n",
    "import numpy as np                                     # arrays and matrix math\n",
    "import scipy.stats as st                               # statistical methods\n",
    "import pandas as pd                                    # DataFrames\n",
    "import matplotlib.pyplot as plt                        # plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the working directory\n",
    "\n",
    "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time).  Also, in this case make sure to place the required (see below) data file in this directory.  When we are done with this tutorial we will write our new dataset back to this directory.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"C:\\PGE383\")                                  # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading Data\n",
    "\n",
    "Let's load the provided dataset. 'PorositySamples2Units.csv' is available at https://github.com/GeostatsGuy/GeoDataSets. It is a comma delimited file with 20 porosity measures from 2 rock units from the subsurface, porosity (as a fraction). We load it with the pandas 'read_csv' function into a data frame we called 'df' and then preview it by printing a slice and by utilizing the 'head' DataFrame member function (with a nice and clean format, see below).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "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>X1</th>\n",
       "      <th>X2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.21</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.17</td>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.15</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.13</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     X1    X2\n",
       "0  0.21  0.20\n",
       "1  0.17  0.26\n",
       "2  0.15  0.20\n",
       "3  0.20  0.19\n",
       "4  0.19  0.13"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"PorositySample2Units.csv\")           # read a .csv file in as a DataFrame\n",
    "#print(df.iloc[0:5,:])                                 # display first 4 samples in the table as a preview\n",
    "df.head()                                              # we could also use this command for a table preview "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is useful to review the summary statistics of our loaded DataFrame.  That can be accomplished with the 'describe' DataFrame member function.  We transpose to switch the axes for ease of visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>X1</th>\n",
       "      <td>20.0</td>\n",
       "      <td>0.1645</td>\n",
       "      <td>0.027810</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.1500</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X2</th>\n",
       "      <td>20.0</td>\n",
       "      <td>0.2000</td>\n",
       "      <td>0.045422</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.1675</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.23</td>\n",
       "      <td>0.30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    count    mean       std   min     25%   50%   75%   max\n",
       "X1   20.0  0.1645  0.027810  0.11  0.1500  0.17  0.19  0.21\n",
       "X2   20.0  0.2000  0.045422  0.11  0.1675  0.20  0.23  0.30"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe().transpose()                              # visualize summary statistics "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we extract the X1 and X2 unit porosity samples from the DataFrame into separate arrays called 'Por1' and 'Por2' for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "por1 = df['X1'].values                                 # extract well 1 porosity to a ndarray\n",
    "por2 = df['X2'].values                                 # extract well 2 porosity to a ndarray          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should first visualize the distributions.  It is convenient to plot the distributions on top of each other.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5QAAAGWCAYAAAAHYcxqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X2cnVV56P3flYQQQ0IIMAQhkUSlsRiDGLQQEEdBXhRR+/gGVauI9LFaQainoB5R69P2VG2xHnskRURFRaC+BAoxVNygTwhIUAkYbQgSTGKHAI5hCDCZ5Dp/7J2wM8xkdvbse78Mv+/nM5+5X9a97mvNDntx3S9rRWYiSZIkSdLuGtfqACRJkiRJncmEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkSZIk1cWEUpIkSZJUFxNKaQQRcU9EdLc6jsEi4hMRcUVleXZEZERMaHVcg1Xien6r45AktQ/71tGxb1U7MaFUW4uI+yPi8Yjoi4ieiPhKRExpZgyZ+cLMLFXi2dHR7K6IOD0ifjlo243DbLug7oCHPvcHIuKOiHgyIi5vZN2jERGfjYjVEfFoRPwqIt7Z6pgkaayzb20M+1apzIRSneB1mTkFeAnwUuBju1tBm1xdvBn444jogh0xHQ5MHrTtaOCWBp97A/Bp4LIG1ztajwGvA6YBfw58PiIWtjYkSXpGsG8dPftWCRNKdZDMXA/cAMwDiIiDImJxRDwSEfdGxHu3l61c7bwmIq6IiE3AuyJiz4i4OCI2VH4ujog9K+X3j4jrIqK3Ut+PI2JcZd/9EXFCRJwMfAR4a+Wq7i8i4s0RsaI6zog4PyK+N0T8G4D7gOMqm14C3EO5M6zeNg64o6qN/x4RGyPiNxHxwTr/dt/JzO8BD49UNiKeFxE3RcTDEfFQRHwjIvap2n9/RPx1RNwVEX+IiG9HxKSq/R+OiN9V/sZnjhDXRZn5q8zclpm3AT+m3OlLkprAvtW+VRotE0p1jIiYBbwG+Fll07eAdcBBwJuAv4uI46sOeT1wDbAP8A3go8BRwIspX718GU9dkT2/UlcXMINy55bV58/MJcDfAd/OzCmZeTiwGJgTEX9cVfTtwNeHacYtPNXBHUf5S/4ng7Ytz8z+Sqd7LfAL4GDgeODciDhp+L9SQwTw95T/rn8MzAI+MajMW4CTgTnAfOBdAJX/Mfhr4NXAocAJNZ804lmUr5LfM5rgJUm1s2+1b5VGy4RSneB7EdFLuXO4mXLnNgs4FvibzHwiM38OXAq8o+q4WzPze5UrdI8DfwZ8KjMfzMyNwCerym8Bng0ckplbMvPHmblTpzeUzHwS+Dbljo6IeCEwG7humEOqr5i+nHKn9+NB226uLL8U6MrMT2Vmf2beB/wb8LaR4hqNzLw3M2/MzCcrf6d/Al4xqNi/ZOaGzHyEcsf84sr2twBfycy7M/Mxnt5Z7sqXKHfwPxhdCyRJNbBvtW+VGsKEUp3gDZm5T2Yekpl/WenADgIeycxHq8qtpXy1cbvfDqrnoEqZ6vIHVZY/A9wLLI2I+2L3Xtz/KnBGRATlTvSqSmc4lFuA+RExnfIV3Vsz81fAsyvbjuWpdzwOAQ6qPCrUW+n4P0L5Km9hIuKAiLgyItZXHmm6Ath/ULH/rlreDGwfzOEgdv67V/+9d3XOz1B+3OottfzPhiRp1Oxb7VulhjChVKfaAOwbEVOrtj0HWF+1PvjLcwPljqS6/AaAzHw0M8/PzOdSfpH9vEGP+AxXJ5m5HOinfAX0DIZ/JIfKldANwNnAA5nZV9l1a2XbFGB5Zdtvgd9UOvztP1Mz8zXD1d8gf0+5nfMzc2/KV4ijxmN/R/kxnu2eM9IBEfFJ4BTgxMzctJuxSpIax761OPatGrNMKNWRMvO3wDLg7yNiUkTMB95D+X2O4XwL+FhEdEXE/sDHKV8hJCJOjYjnV66EbgK2Vn4G6wFmbx9UoMrXgP8NDGTmT0YI/8fAeZXf2/2ksu2OylVigNuBTRHxNxHxrIgYHxHzIuKlI9T/NBExofJy/3hgfOVvNtzofFOBPqA3Ig4GPrwbp7qK8iANh0XEZOCiEeK6kPL/KLw6M0cc1ECSVBz71t1j3yqVmVCqk51O+Z2KDcB3gYsy88ZdlP805RHe7gJWAndWtkH5Jff/pPxlfyvwr1mZH2uQqyu/H46IO6u2f53yYyXDXkGtcjNwAOWObrsfV7btGNI8M7dSvqL7YuA3wEOU32WZVsM5BvsY8DhwAeWroo8z/BDxn6Q8It4fgP8AvlPrSTLzBuBi4CbKjzndNMIhf0f5SuvqKI/u1xcRH6n1fJKkhrNvrZ19qwSEj1RLo1cZRe1B4CWZubrV8UiS1OnsW6XO4B1KqTHeB/zUDk+SpIaxb5U6wHDPeTdERHwIOIvyS8grgXdn5hNFnlNqtoi4n/KL9W9ocSiSOlhEXAacCjyYmfOG2B/A5ynPGbgZeFdm3jm4nDQW2LdKnaOwO5SVF44/CBxZ6RjHU/AcP1IrZObsyrDrPxu5tCQN63LKk5oP5xTK76QdSnnkyv/ThJiklrBvlTpH0Y+8TgCeVRnxajKVYaQlSdLOMvMW4JFdFHk98LUsWw7sExHPbk50kiQNrbBHXjNzfUR8FniA8qhXSzNz6eByEXE25SutTJo0acFznjPi1Dptbdu2bYwb19mvptqG9mAb2oNtaA//9V//9VBmdrU6jhY7mJ0nN19X2fa7wQXtW9uPbWgPg9uwbds26LABKpPaJ7BstW2ZQ8wyCjEuyG2d9XcfbCy0AeC+39w36v61sIQyIqZTvpo6B+gFro6It2fmFdXlMnMRsAhg7ty5+etf/7qokJqiVCrR3d3d6jBGxTa0B9vQHmxDe4iIta2OoQ0M9f+QQ/7fjH1r+7EN7WFwG1YsWcKCrs66VlXq6aF7xoxWh1GTJctW0zX98Kdt75m8lhmbD2lBRI0zFtoAcOQ7Dht1/1rkZaYTgN9k5sbM3EJ5vp2FBZ5PkqSxbB0wq2p9Jr5KIklqsSITygeAoyJicmVkuuOBVQWeT5KksWwx8M4oOwr4Q2Y+7XFXSZKaqch3KG+LiGuAO4EB4GdUHr+RJEk7i4hvAd3A/hGxDrgI2AMgM78EXE95ypB7KU8b8u7WRCpJ0lMKnYcyMy+i3CHWbcuWLaxbt44nnuiM6SunTZvGqlWNvxE7adIkZs6cyR577NHwuiVJrZeZp4+wP4H3N+Jc9q1l9q2SNHqFJpSNsG7dOqZOncrs2bMpPznb3h599FGmTp3a0Dozk4cffph169YxZ86chtYtSXrmsW+1b5WkRmn7sZ+feOIJ9ttvv47o8IoSEey3334dcyVZktTe7FvtWyWpUdo+oQSe0R3edv4NJEmNZL/i30CSGqEjEkpJkiRJUvtp+3coB1u5fDn9vb0Nq2/iPvvwoqOOGnb/hz70IQ455BDOPfdcAE466SRmzZrFpZdeCsD555/PwQcfzHnnnTdsHVOmTKGvr4/777+fU089lbvvvvtpZU4++WSWL1/Osccey3XXXTfKVkmSVDv7VklSvTouoezv7WVBV1fD6luxceMu9y9cuJCrr76ac889l23btvHQQw+xadOmHfuXLVvGxRdfPOo4PvzhD7N582YuueSSUdclSdLusG+VJNXLR15HcMwxx7Bs2TIA7rnnHubNm8fUqVP5/e9/z5NPPsmqVas44ogjAPjMZz7DK17xCubPn89FF+3ebCnHH398w0ewkySpHdm3StLY0XF3KJvtoIMOYsKECTzwwAMsW7aMo48+mvXr13Prrbcybdo05s+fz8SJE1m6dCmrV6+mVCoxZcoUTjvtNG655RaOO+64VjdBkqS2Yt8qSWOHCWUNtl9JXbZsGeeddx7r169n2bJlTJs2jYULFwKwdOlSli5dyrHHHsu4cePo6+tj9erVdnqSJA3BvlWSxgYTyhosXLiQZcuWsXLlSubNm8esWbP43Oc+x957782ZZ54JlCdIvvDCCznjjDN8vEaSpBHYt0rS2OA7lDU45phjuO6669h3330ZP348++67L729vdx6660cffTRQHmEussuu4y+vj4A1q9fz4MPPtjKsCVJalv2rZI0NnTcHcqJ++wz4uhxu1vfSF70ohfx0EMPccYZZ+y0ra+vj/333x+AE088kVWrVnHCCScwbtw4pkyZwhVXXMEBBxxQUxwvf/nL+dWvfkVfXx8zZ87ky1/+MieddFJ9jZIkaTfYt2qsGjwlzub+flYsWbJjfeKWLa0Iq2GW37WG3r5trQ5jWFsGJrY6BDVBxyWUu5rXqijjx4/faThzgMsvv/xp5c455xzOPPPMpz2Ws/3K6uzZs4ecJwvgxz/+cWOClSRpN9m3aqwaPCVOqaenoVPktFpv3za6ph/e6jD0DOcjr5IkSZKkuphQSpIkSZLqYkIpSZIkSaqLCaUkSZIkqS4mlJIkSZKkuphQSpIkSZLq0nHThixfvpLe3v6G1bfPPhM56qgXDbv/Qx/6EIcccgjnnnsuUJ5kedasWVx66aUAnH/++Rx88MGcd955w9YxZcoU+vr6uP/++zn11FOfNrz5z3/+c973vvexadMmxo8fz0c/+lHe+ta3NqB1kiSNzL5VklSvjksoe3v76epa0LD6Nm5cscv9Cxcu5Oqrr+bcc89l27ZtPPTQQzvNm7Vs2TIuvvjiUcUwefJkvva1r3HooYeyYcMGFixYwEknncQ+NUwMLUnSaNm3SpLq5SOvIzjmmGNYtmwZAPfccw/z5s1j6tSp/P73v+fJJ59k1apVHHHEEQB85jOf4RWveAXz58/noosuqvkcf/RHf8Shhx4KwEEHHcQBBxzAxo0bG98YSZLagH2rJI0dHXeHstkOOuggJkyYwAMPPMCyZcs4+uijWb9+PbfeeivTpk1j/vz5TJw4kaVLl7J69WpKpRJTpkzhtNNO45ZbbuG4447brfPdfvvt9Pf387znPa+gFkmS1Fr2rZI0dphQ1mD7ldRly5Zx3nnnsX79epYtW8a0adNYuHAhAEuXLmXp0qUce+yxjBs3jr6+PlavXr1bnd7vfvc73vGOd/DVr36VceO8eSxJGrvsWyVpbDChrMHChQtZtmwZK1euZN68ecyaNYvPfe5z7L333px55pkAZCYXXnghZ5xxBlOnTt3tc2zatInXvva1fPrTn+aoo45qdBMkSWor9q2SNDZ4qa4GxxxzDNdddx377rsv48ePZ99996W3t5dbb72Vo48+GiiPUHfZZZfR19cHwPr163nwwQdrqr+/v583vvGNvPOd7+TNb35zYe2QJKld2LdK0tjQcXco99ln4oijx+1ufSN50YtexEMPPcQZZ5yx07a+vj72339/AE488URWrVrFCSecwLhx45gyZQpXXHEFBxxwwIj1X3XVVdxyyy08/PDDXH755QBcfvnlvPjFL66vUZIk7Qb7VklSvTouodzVvFZFGT9+/E7DmQM7Oqdq55xzDmeeeebTHsvZfmV19uzZT5snC+Dtb387b3/72xsXsCRJu8G+VZJULx95lSRJkiTVpePuUEqSJEkAK5cvp7+3d9j9E7dsaWI0jbf8rjX09m3baVv/jD1ZsmY1AFsGRn68XCpaRySUmUlEtDqMlsrMVocgSRpD7FvtW8eC/t5eFnR1tTqMwvT2baNr+uE7beuZsJau6S9oUUTS07X9I6+TJk3i4YcffkZ/6WcmDz/8MJMmTWp1KJKkMcC+1b5Vkhql7e9Qzpw5k3Xr1rFx48ZWh1KTJ554opDOadKkScycObPh9UqSnnnsW8vsWyVp9No+odxjjz2YM2dOq8OoWalU4ogjjmh1GJIkDcu+VZLUKIU98hoRcyPi51U/myLi3KLOJ0mSJElqrsLuUGbmr4EXA0TEeGA98N2izidJkiRJaq5mDcpzPLAmM9c26XySJEmSpII16x3KtwHfGmpHRJwNnA3Q1dVFqVRqUkjF6Ovrsw1toJ3b8Nhjj7N168gjK2YOsHjx9SOWGz8+2GuvZzUitIZr58+hVrZBkiRpeIUnlBExETgNuHCo/Zm5CFgEMHfu3Ozu7i46pEKVSiVsQ+u1cxuWLFnBgQcuGLFcT0+JGTO6Ryy3ceMKurtHrq8V2vlzqJVtkCRJGl4zHnk9BbgzM3uacC5JkiRJUpM0I6E8nWEed5UkSZIkda5CE8qImAy8GvhOkeeRJEmSJDVfoe9QZuZmYL8izyFJkiRJao1mTRsiSZIkSRpjTCglSZIkSXUxoZQkSZIk1aXweSglSZKkZ6Lld62ht29b3cdvGZjYwGikYphQSpIkSQXo7dtG1/TDWx2GVCgfeZUkSZIk1cWEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkSZIk1cWEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkSZIk1cWEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkqU1ExMkR8euIuDciLhhi/3Mi4kcR8bOIuCsiXtOKOCVJ2s6EUpKkNhAR44EvAqcAhwGnR8Rhg4p9DLgqM48A3gb8a3OjlCRpZyaUkiS1h5cB92bmfZnZD1wJvH5QmQT2rixPAzY0MT5Jkp5mQqsDkCRJABwM/LZqfR3wJ4PKfAJYGhF/BewFnDBURRFxNnA2QFdXF6VSqdGxNlVfX59taAPt2IbN/f2UenpqLt83MLBb5Uerf8ae9ExY29A6B8b10zO5sXU2m20YW0woJUlqDzHEthy0fjpweWZ+LiKOBr4eEfMyc9tOB2UuAhYBzJ07N7u7u4uIt2lKpRK2ofXasQ0rlixhQVdXzeVLPT10z5hRYEQ7W7JmNV3TX9DQOnsmr2XG5kMaWmez2YaxxUdeJUlqD+uAWVXrM3n6I63vAa4CyMxbgUnA/k2JTpKkIZhQSpLUHn4KHBoRcyJiIuVBdxYPKvMAcDxARPwx5YRyY1OjlCSpigmlJEltIDMHgA8APwBWUR7N9Z6I+FREnFYpdj7w3oj4BfAt4F2ZOfixWEmSmsZ3KCVJahOZeT1w/aBtH69a/iVwTLPjkiRpON6hlCRJkiTVxYRSkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVxYRSkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVpdCEMiL2iYhrIuJXEbEqIo4u8nySJEmSpOaZUHD9nweWZOabImIiMLng80mSJEmSmqSwhDIi9gaOA94FkJn9QH9R55MkSZIkNVeRdyifC2wEvhIRhwMrgHMy87HqQhFxNnA2QFdXF6VSqcCQitfX12cb2kA7t+HRTY+wedN/j1hu3B4DrF3zHyOW28qTlEqPNiK0hmvnz6FWtkGSJGl4RSaUE4CXAH+VmbdFxOeBC4D/WV0oMxcBiwDmzp2b3d3dBYZUvFKphG1ovXZuwxc+fQkLZx83YrmeCWuZ0X/IiOWW3X8Lb37LnzYitIZr58+hVrZBkiRpeEUOyrMOWJeZt1XWr6GcYEqSJEmSxoDCEsrM/G/gtxExt7LpeOCXRZ1PkiRJktRcRY/y+lfANyojvN4HvLvg80mSJEmSmqTQhDIzfw4cWeQ5JEmSJEmtUeQ7lJIkSZKkMcyEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkSZIk1cWEUpIkSZJUFxNKSZIkSVJdTCglSZIkSXUxoZQkSZIk1cWEUpIkSZJUFxNKSZIkSVJdJrQ6AEmSJGkoK5cvp7+3d9j9E7dsKfT8y+9aQ2/ftrqP3zIwsYHRSO3JhFKSJEltqb+3lwVdXS07f2/fNrqmH96y80udwEdeJUmSJEl1MaGUJEmSJNXFhFKSJEmSVBcTSkmSJElSXUwoJUmSJEl1MaGUJEmSJNXFhFKSJEmSVBcTSkmSJElSXUwoJUmSJEl1MaGUJEmSJNXFhFKSJEmSVBcTSkmSJElSXUwoJUmSJEl1MaGUJKmBIuKOiHh/RExvdSySJBXNhFKSpMZ6G3AQ8NOIuDIiToqIaHVQkiQVYUKrA5AkaSzJzHuBj0bE/wROBS4DtkXEZcDnM/ORlgYoNdHK5cvp7+2t+/iJW7Y0MJqnW37XGnr7tg27f8vAxELPL40FJpSSJDVYRMwH3g28Bvh34BvAscBNwItbGJrUVP29vSzo6mp1GMPq7dtG1/TDWx2G1NFMKCVJaqCIWAH0Al8GLsjMJyu7bouIY1oXmSRJjWdCKUlSY705M++r3hARczLzN5n5p60KSpKkIjgojyRJjXVNjdskSep43qGUJKkBIuIFwAuBaRFRfSdyb2BSa6KSJKlYhSaUEXE/8CiwFRjIzCOLPJ8kSS00l/KorvsAr6va/ijw3pZEJElSwZpxh/KVmflQE84jSVLLZOb3ge9HxNGZeWur45EkqRl85FWSpAaIiP+Rmf8InBERpw/en5kfbEFYkiQVquiEMoGlEZHAJZm5aHCBiDgbOBugq6uLUqlUcEjF6uvrsw0t8vhjj5FbtwIwkMn1ixePus7+/n4mTmzspMb7zZpKz+S1I5YbGNdfU7kZz9uvbT+vTv23VM02aDesqvy+o6VRSJLUREUnlMdk5oaIOAC4MSJ+lZm3VBeoJJmLAObOnZvd3d0Fh1SsUqmEbWiNFUuWsODAAwEo9fTQPWPGqOu88oYbeMMpp4y6nmoXXPNN3nziWSOW65m8lhmbDxmx3Jr7b+Etp7+pEaE1XKf+W6pmG1SrzLy28vurrY5FkqRmKTShzMwNld8PRsR3gZcBt+z6KEmSOk9EXEv5yZwhZeZpTQxHkqSmKCyhjIi9gHGZ+Whl+UTgU0WdT5KkFvtsqwOQJKnZirxDOQP4bkRsP883M3NJgeeTJKllMvPmVscgSVKzFZZQZuZ9wOFF1S9JUjuJiKsy8y0RsZKdH30NIDNzfg11nAx8HhgPXJqZ/zBEmbcAn6ic4xeZeUYj4pckqR5OGyJJUmOcU/l9aj0HR8R44IvAq4F1wE8jYnFm/rKqzKHAhZQHvft9ZdA7SZJaZlyrA5AkaSzIzN9Vfq8FnqT8lM584MnKtpG8DLg3M+/LzH7gSuD1g8q8F/hiZv6+cq4HGxW/JEn18A6lJEkNFBFnAR8HbqL8uOsXIuJTmXnZCIceDPy2an0d8CeDyvxR5Rz/P+XHYj8x1PgEzvHcfp6pbdjc30+pp6eYgOrQNzCwUzz9M/akZ0It13vaR63zVLcz2zC2mFBKktRYHwaOyMyHASJiP2AZMFJCGUNsGzwNyQTgUKAbmAn8OCLmZWbvTgc5x3Pbeaa2YcWSJSzo6iomoDoMnqd6yZrVdE1/QQsj2n21zlPdzmzD2OIjr5IkNdY64NGq9UfZ+c7jro6bVbU+E9gwRJnvZ+aWzPwN8GvKCaYkSS3hHUpJkhogIs6rLK4HbouI71O+w/h64PYaqvgpcGhEzKnU8TZg8Aiu3wNOBy6PiP0pPwJ7XwPClySpLiaUkiQ1xtTK7zWVn+2+X8vBmTkQER8AfkD5/cjLMvOeiPgUcEdmLq7sOzEifglsBT68/dFaSZJawYRSkqQGyMxPNqCO64HrB237eNVyAudVfiRJajkTSkmSGigiuoD/AbwQmLR9e2a+qmVBSZJUEAflkSSpsb4B/AqYA3wSuJ/y+5GSJI05JpSSJDXWfpn5ZWBLZt6cmWcCR7U6KEmSiuAjr5IkNdaWyu/fRcRrKU/9MbOF8UiSVBgTSkmSGuvTETENOB/4ArA38KHWhiRJUjFqSigjYl5m3l10MJIkdbrMvK6y+Afgla2MRZKkotX6DuWXIuL2iPjLiNin0IgkSepgEfHciLg2Ih6KiAcj4vsR8dxWxyVJUhFqSigz81jgz4BZwB0R8c2IeHWhkUmS1Jm+CVwFHAgcBFwNfKulEUmSVJCaR3nNzNXAx4C/AV4B/EtE/Coi/rSo4CRJ6kCRmV/PzIHKzxVAtjooSZKKUOs7lPOBdwOvBW4EXpeZd0bEQcCtwHeKC1GSpPYXEftWFn8UERcAV1JOJN8K/EfLApMkqUC1jvL6v4F/Az6SmY9v35iZGyLiY4VEJklSZ1lBOYGMyvpfVO1L4G+bHpEkSQWrNaF8DfB4Zm4FiIhxwKTM3JyZXy8sOkmSOkRmzml1DJIkNVutCeV/AicAfZX1ycBSYGERQUmS1KkiYg/gfcBxlU0l4JLM3NKyoKRnqOV3raG3b9uO9f4Ze7Jkzeod61sGJrYiLGlMqTWhnJSZ25NJMrMvIiYXFJMkSZ3s/wB7AP9aWX9HZdtZLYtIeobq7dtG1/TDd6z3TFhL1/QXtDAiaeypNaF8LCJekpl3AkTEAuDxEY6RJOmZ6KWZeXjV+k0R8YuWRSNJUoFqTSjPBa6OiA2V9WdTHrVOkiTtbGtEPC8z1wBExHOBrS2OSZKkQtSUUGbmTyPiBcBcyqPX/cp3QSRJGtKHKU8dch/lPvMQylNvSZI05tR6hxLgpcDsyjFHRASZ+bVCopIkqQNVRkF/HDiUnS/CPtnSwCRJKkhNCWVEfB14HvBznnpsJwETSkmSKjJzW0R8LjOPBu5qdTySJBWt1juURwKHZWYWGYwkSWPA0oj4f4Dv2G9Kksa6WhPKu4EDgd8VGIskSWPBecBewEBEPEH5sdfMzL1bG5YkSY1Xa0K5P/DLiLgd2PEeSGaeVkhUkiR1qMyc2uoYJElqlloTyk8UGYQkSZ0uIg4APgI8n/L7k/+QmZtaG5UkScUaV0uhzLwZuB/Yo7L8U+DOAuOSJKnTfA14DPgCMBX4l9aGI0lS8Wod5fW9wNnAvpRHez0Y+BJwfHGhSZLUUQ7MzI9Wln8QEV54lSSNebU+8vp+4GXAbQCZubryaI8kSSqLiJhOeRAegPHV65n5SMsikySpILUmlE9mZn9EuY+MiAmU56EcUUSMB+4A1mfmqXVFKUlS+5sGrOCphBKeej0kgec2PSJJkgpWa0J5c0R8BHhWRLwa+Evg2hqPPQdYBThcuiRpzMrM2a2OQZKkZqtpUB7gAmAjsBL4C+B64GMjHRQRM4HXApfWG6AkSZIkqT3VdIcyM7cB/1b52R0XA/+D8mh3Q4qIsykP+ENXVxelUmk3T9Fe+vryeuGXAAAeEUlEQVT6bEOLbO7vp9TTA0DfwMCO5dHY83nPa0g91V70qiPombx2xHID4/prKjf9OXuzePH1jQhth/Hjg732etao6+nUf0vVbIMkSdLwah3l9TcM8c5kZg77PkhEnAo8mJkrIqJ7uHKZuQhYBDB37tzs7h62aEcolUrYhtZYsWQJC7q6ACj19NA9Y8ao67zyzjvpPuWUUddT7YIrf8ibTzxrxHI9k9cyY/MhI5a79j9v56xzPt6I0HbYuHEF3d0LRl1Pp/5bqmYbJEmShlfrO5RHVi1PAt5MeQqRXTkGOC0iXlM5Zu+IuCIz3777YUqS1Dki4ljg0Mz8SkR0AVMy8zetjkuSpEar6R3KzHy46md9Zl4MvGqEYy7MzJmVQQreBtxkMilJGusi4iLgb4ALK5v2AK5oXUSSJBWn1kdeX1K1Oo7yHcth34uUJOkZ7I3AEVSmDMnMDRFhnylJGpNqfeT1c1XLA8D9wFtqPUlmloBSreUlSepg/ZmZEZEAEbFXqwOSJKkotY7y+sqiA5EkaYy4KiIuAfaJiPcCZ7L7o6RLktQRan3k9bxd7c/Mf2pMOJIkdbbM/GxEvBrYBMwFPp6ZN7Y4LEmSCrE7o7y+FFhcWX8dcAvw2yKCkiSpU0XEh4CrTSIlSc8EtSaU+wMvycxHASLiE5Q7y5En05Mk6Zllb+AHEfEIcCVwTWb2tDgmSZIKUdO0IcBzgP6q9X5gdsOjkSSpw2XmJzPzhcD7gYOAmyPiP1scliRJhaj1DuXXgdsj4rtAUh4S/WuFRSVJUud7EPhv4GHggBbHIklSIWod5fX/i4gbgJdXNr07M39WXFiSJHWmiHgf8FagC7gGeG9m/rK1UUmSVIxa71ACTAY2ZeZXIqIrIuZk5m+KCkySpA51CHBuZv681YFIklS0WqcNuYjySK9zga8AewBXAMcUF5okSZ0jIvbOzE3AP1bW963en5mPtCQwSZIKVOsdyjcCRwB3AmTmhoiYWlhUkiR1nm8CpwIrKI83EFX7EnhuK4KSJKlItSaU/ZmZEZEAEbFXgTFJktRxMvPUyu85rY5FkqRmqTWhvCoiLgH2iYj3AmcC/1ZcWJIkdaaI+GFmHj/SNmksWLl8Of29vcPun7hlSxOjkdQKtY7y+tmIeDWwifJ7lB/PzBsLjUySpA4SEZMoD2C3f0RM56lHXvemPB+lNOb09/ayoKur1WFIaqERE8qIGA/8IDNPAEwiJUka2l8A51JOHlfwVEK5Cfhiq4KSJKlIIyaUmbk1IjZHxLTM/EMzgpIkqdNk5ueBz0fEX2XmF1odjyRJzVDrO5RPACsj4kbgse0bM/ODhUQlSVKHyswvRMQ84DBgUtX2r7UuKkmSilFrQvkflR9JkrQLlbmbuyknlNcDpwA/AUwoJUljzi4Tyoh4TmY+kJlfbVZAkiR1uDcBhwM/y8x3R8QM4NIWxyRJUiHGjbD/e9sXIuLfC45FkqSx4PHM3AYMRMTewIPAc1sckyRJhRjpkdeoWrYzlCRpZHdExD6U52teAfQBt7c2JEmSijFSQpnDLEuSpCFk5l9WFr8UEUuAvTPzrlbGJElSUUZKKA+PiE2U71Q+q7JMZT0zc+9Co5MkqUNExEt2tS8z72xmPJIkNcMuE8rMHN+sQCRJ6nCf28W+BF7VrEAkSWqWWqcNkSRJu5CZr2x1DJIkNZsJpSRJDRQR7xxqe2Y6D6UkacwxoZQkqbFeWrU8CTgeuBMwoZQkjTkmlJIkNVBm/lX1ekRMA75ey7ERcTLweWA8cGlm/sMw5d4EXA28NDPvGF3EkiTVb1yrA5AkaYzbDBw6UqGIGA98ETgFOAw4PSIOG6LcVOCDwG0NjlOSpN3mHUpJkhooIq7lqbmbx1FODq+q4dCXAfdm5n2Veq4EXg/8clC5vwX+EfjrhgQsSdIomFBKktRYn61aHgDWZua6Go47GPht1fo64E+qC0TEEcCszLwuIoZNKCPibOBsgK6uLkqlUo2ht6e+vj7b0AaGasPm/n5KPT2tCQh47PEn2bp1+P05YxI9E9buWB8Y10/P5LXDH9ABbEN7GAttaBQTSkmSGigzbwaIiL2p9LMRsW9mPjLCoTFUdTt2RowD/hl4Vw0xLAIWAcydOze7u7trCb1tlUolbEPrDdWGFUuWsKCrqzUBAUuWrebA6YfvulD/U4s9k9cyY/MhxQZVMNvQHsZCGxrFhFKSpAaq3B38W+BxYBvlRDGB545w6DpgVtX6TGBD1fpUYB5QigiAA4HFEXGaA/NIklrFhFKSpMb6MPDCzHxoN4/7KXBoRMwB1gNvA87YvjMz/wDsv309IkrAX5tMSpJayVFeJUlqrDWUR3bdLZk5AHwA+AGwCrgqM++JiE9FxGkNjlGSpIYo7A5lREwCbgH2rJznmsy8qKjzSZLUJi4ElkXEbcCT2zdm5gdHOjAzrweuH7Tt48OU7R5dmJIkjV6Rj7w+CbwqM/siYg/gJxFxQ2YuL/CckiS12iXATcBKyu9QSpI0ZhWWUGZmAn2V1T0qPzn8EZIkjQkDmXleq4OQJKkZCh2UJyLGAyuA5wNfzMzbhijjXFltplPb8PCmPtZvehyA3CNYvOaBUdeZz5nd8Pm1XvSqI2qat6jW+Y2OftWh9PSUGhBZ1bkHNjfk30Cn/luqZhtUhx9V+rZr2fmR15GmDZEkqeMUmlBm5lbgxRGxD/DdiJiXmXcPKuNcWW2mU9vwhU9fwsLZxwHQM2EtM/pHPzfQ1f95Kaefc8bIBXfDBVf+kDefeNaI5Wqd3+jam27nrHNOb0RoO2zcuILu7gWjrqdT/y1Vsw2qw/YvjQurttUybYgkSR2nKdOGZGZvZXjzk4G7RyguSVLHysw5rY5BkqRmKXKU1y5gSyWZfBZwAvC/ijqfJEntICLeOdT2zPxas2ORJKloRd6hfDbw1cp7lOMoz6d1XYHnkySpHby0ankScDxwJ2BCKUkac4oc5fUu4Iii6pckqR1l5l9Vr0fENODrLQpHkqRCjWt1AJIkjXGbgUNbHYQkSUVoyqA8kiQ9U0TEtTw17/I44DDgqtZFJElScUwoJUlqrM9WLQ8AazNzXauCkSSpSCaUkiQ1QEQ8H5iRmTcP2v7yiNgzM9e0KDSpbiuXL6e/txeAzf39rFiyZKf9E7dsaUVYktqICaUkSY1xMfCRIbY/Xtn3uuaGI41ef28vC7q6ACj19OxYlqTtHJRHkqTGmF0Z4XwnmXkHMLv54UiSVDwTSkmSGmPSLvY9q2lRSJLURCaUkiQ1xk8j4r2DN0bEe4AVLYhHkqTC+Q6lJEmNcS7w3Yj4M55KII8EJgJvbFlUkiQVyIRSkqQGyMweYGFEvBKYV9n8H5l5UwvDkiSpUCaUkiQ1UGb+CPhRq+OQJKkZfIdSkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVxYRSkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVxYRSkiRJklSXCa0OQJIkSWPT8rvW0Nu3rbD6twxMLKxuSbUxoZQkSVIhevu20TX98FaHIalAPvIqSZIkSaqLCaUkSZIkqS4mlJIkSZKkuphQSpIkSZLqYkIpSZIkSaqLCaUkSZIkqS4mlJIkSZKkuphQSpIkSZLqYkIpSZIkSaqLCaUkSZIkqS4mlJIkSZKkuphQSpIkSZLqUlhCGRGzIuJHEbEqIu6JiHOKOpckSZIkqfkmFFj3AHB+Zt4ZEVOBFRFxY2b+ssBzSpIkSZKapLA7lJn5u8y8s7L8KLAKOLio80mSJEmSmqvIO5Q7RMRs4AjgtiH2nQ2cDdDV1UWpVGpGSIXp6+sbk214/LHHyK1bG3qe/v5+Jk6c2LD69ps1lZ7JawEYGNe/Y3k0XvSqIyj19Iy6nsF11hJbrW04+lWH0tNTakBkVece2NyQf8dj9b+HTjMW2iCpPS2/aw29fduG3b9loHH9vKT2VHhCGRFTgH8Hzs3MTYP3Z+YiYBHA3Llzs7u7u+iQClUqlRiLbVixZAkLDjywoee58oYbeMMppzSsvguu+SZvPvEsAHomr2XG5kNGXectN93In51zxqjrqXbBlT/cEeeu1NqGa2+6nbPOOb0Roe2wceMKursXjLqesfrfQ6cZC22Q1J56+7bRNf3wVochqYUKHeU1IvagnEx+IzO/U+S5JEmSJEnNVeQorwF8GViVmf9U1HkkSZIkSa1R5B3KY4B3AK+KiJ9Xfl5T4PkkSZIkSU1U2DuUmfkTIIqqX5IkSZLUWoW+QylJkiRJGrtMKCVJkiRJdTGhlCRJkiTVxYRSkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVxYRSkiRJklQXE0pJkiRJUl1MKCVJahMRcXJE/Doi7o2IC4bYf15E/DIi7oqIH0bEIa2IU5Kk7UwoJUlqAxExHvgicApwGHB6RBw2qNjPgCMzcz5wDfCPzY1SkqSdmVBKktQeXgbcm5n3ZWY/cCXw+uoCmfmjzNxcWV0OzGxyjJIk7WRCqwOQJEkAHAz8tmp9HfAnuyj/HuCGoXZExNnA2QBdXV2USqUGhdgafX19tqEgjz/2GLl167D7I5NSTw8AfQMDO5a365+xJz0T1hYaYyMNjOunZ3LnxDsU29AexkIbGsWEUpKk9hBDbMshC0a8HTgSeMVQ+zNzEbAIYO7cudnd3d2gEFujVCphG4qxYskSFhx4YE1lSz09dM+YsdO2JWtW0zX9BUWEVoieyWuZsbmzXz22De1hLLShUUwoJUlqD+uAWVXrM4ENgwtFxAnAR4FXZOaTTYpNkqQh+Q6lJEnt4afAoRExJyImAm8DFlcXiIgjgEuA0zLzwRbEKEnSTkwoJUlqA5k5AHwA+AGwCrgqM++JiE9FxGmVYp8BpgBXR8TPI2LxMNVJktQUPvIqSVKbyMzrgesHbft41fIJTQ9KkqRd8A6lJEmSJKkuJpSSJEmSpLqYUEqSJEmS6mJCKUmSJEmqiwmlJEmSJKkuJpSSJEmSpLqYUEqSJEmS6mJCKUmSJEmqiwmlJEmSJKkuJpSSJEmSpLqYUEqSJEmS6jKh1QFIkiSpPiuXL6e/t7fu4ydu2bLL/cvvWkNv3zYA+mfsyZI1q3fav2VgYt3nljQ2mFBKkiR1qP7eXhZ0dRVWf2/fNrqmHw5Az4S1dE1/QWHnktSZfORVkiRJklQXE0pJkiRJUl1MKCVJkiRJdTGhlCRJkiTVpbCEMiIui4gHI+Luos4hSZIkSWqdIu9QXg6cXGD9kiRJkqQWKiyhzMxbgEeKql+SJEmS1Fotn4cyIs4Gzgbo6uqiVCq1NqBR6uvrG5NteHhTH+s3Pd7Q8zz57FksXvNAw+qb+/L59ExeC8DAuP4dy6Ots5Exbq+zlthqbcPRrzqUnp5SAyJ7Sl/f71m8uGfU9WQOsHjx9QD09/czcWLjJsAePz7Ya69nNaw+gMcee5ytW3OnbdVtqEcRce6usfC9JEmS2lPLE8rMXAQsApg7d252d3e3NqBRKpVKjMU2fOHTl7Bw9nENPc/Xv3kxF7z73MbVd+VT9fVMXsuMzYeMus6vXPndhsYIO8e5K7W24dqbbuesc05vRGg73HnnlZxyyhtHXU9PT4kZM7oBuOGGKznllDeMus7tNm5cQXf3gobVB7BkyQoOPHDnOqvbUI8i4txdY+F7SZIktSdHeZUkSZIk1cWEUpIkSZJUlyKnDfkWcCswNyLWRcR7ijqXJEmSJKn5CnuHMjMb+1KXJEmSJKmt+MirJEmSJKkuJpSSJEmSpLq0fNoQSZKkZ6qVy5fT39tb9/GrVt7Hxkn1Hz+SLQONmz9Y0thkQilJktQi/b29LOjqqvv4jZN66Zp+eAMjkqTd4yOvkiRJkqS6mFBKkiRJkupiQilJkiRJqosJpSRJkiSpLiaUkiRJkqS6mFBKkiRJkupiQilJkiRJqosJpSRJkiSpLiaUkiRJkqS6mFBKkiRJkupiQilJkiRJqosJpSRJkiSpLiaUkiRJkqS6mFBKkiRJkupiQilJkiRJqosJpSRJkiSpLhNaHYAkSdIz1V2rHmDj6t66j98yMLGB0UjS7jOhlCRJapG+x5Ku2Ye3OgxJqpuPvEqSJEmS6mJCKUmSJEmqiwmlJEmSJKkuJpSSJEmSpLqYUEqSJEmS6mJCKUmSJEmqiwmlJEmSJKkuJpSSJEmSpLqYUEqSJEmS6mJCKUmSJEmqiwmlJEmSJKkuJpSSJEmSpLoUmlBGxMkR8euIuDciLijyXJIkdbqR+s2I2DMivl3Zf1tEzG5+lJIkPaWwhDIixgNfBE4BDgNOj4jDijqfJEmdrMZ+8z3A7zPz+cA/A/+ruVFKkrSzIu9Qvgy4NzPvy8x+4Erg9QWeT5KkTlZLv/l64KuV5WuA4yMimhijJEk7icwspuKINwEnZ+ZZlfV3AH+SmR8YVO5s4OzK6jzg7kICap79gYdaHcQo2Yb2YBvag21oD3Mzc2qrgyhSLf1mRNxdKbOusr6mUuahQXXZt7Yf29AebEN7sA3tY9T964RGRTKEoa6YPi17zcxFwCKAiLgjM48sMKbC2Yb2YBvag21oD2OlDa2OoQlq6TftWzuUbWgPtqE92Ib20Yj+tchHXtcBs6rWZwIbCjyfJEmdrJZ+c0eZiJgATAMeaUp0kiQNociE8qfAoRExJyImAm8DFhd4PkmSOlkt/eZi4M8ry28Cbsqi3l2RJKkGhT3ympkDEfEB4AfAeOCyzLxnhMMWFRVPE9mG9mAb2oNtaA+2oQMM129GxKeAOzJzMfBl4OsRcS/lO5Nvq6HqsfC3sw3twTa0B9vQHsZCG6AB7ShsUB5JkiRJ0thW5COvkiRJkqQxzIRSkiRJklSXpiSUEXFyRPw6Iu6NiAuG2L9nRHy7sv+2iJhdte/CyvZfR8RJzYh3KPW2ISJmR8TjEfHzys+Xmh17VYwjteG4iLgzIgYq86FV7/vziFhd+fnzwcc2yyjbsLXqc2jpAFE1tOO8iPhlRNwVET+MiEOq9nXKZ7GrNrTFZ1FDG/7fiFhZifMnEXFY1b5O+W4asg2d9N1UVe5NEZERcWTVtrb4HFphLPStlVjsXzvjO73t+1f71o75HOxbm6CpfWtmFvpDeWCBNcBzgYnAL4DDBpX5S+BLleW3Ad+uLB9WKb8nMKdSz/iiY25wG2YDdzc75jrbMBuYD3wNeFPV9n2B+yq/p1eWp3dSGyr7+lr9OexGO14JTK4sv6/q31MnfRZDtqFdPosa27B31fJpwJLKcid9Nw3Xho75bqqUmwrcAiwHjmynz6Fd/260ed/agHZ0zL9h7F/bpQ32re3RBvvWNmhDpVxD+tZm3KF8GXBvZt6Xmf3AlcDrB5V5PfDVyvI1wPEREZXtV2bmk5n5G+DeSn3NNpo2tIsR25CZ92fmXcC2QceeBNyYmY9k5u+BG4GTmxH0IKNpQzuppR0/yszNldXllOejg876LIZrQ7uopQ2bqlb34qkJ5Dvmu2kXbWgXtXy/Avwt8I/AE1Xb2uVzaIWx0LeC/Wsnfae3e/9q39oe7FvbQ1P71mYklAcDv61aX1fZNmSZzBwA/gDsV+OxzTCaNgDMiYifRcTNEfHyooMdxmj+lp30OezKpIi4IyKWR8QbGhvabtnddrwHuKHOY4symjZAe3wWNbUhIt4fEWsof+F+cHeObYLRtAE65LspIo4AZmXmdbt77Bg2FvpWaozF/rV4Y6F/tW/toM/BvrVwTe1bC5uHsspQVxEHZ/HDlanl2GYYTRt+BzwnMx+OiAXA9yLihYOubDTDaP6WnfQ57MpzMnNDRDwXuCkiVmbmmgbFtjtqbkdEvB04EnjF7h5bsNG0Adrjs6ipDZn5ReCLEXEG8DHKk8p31OcwTBs64rspIsYB/wy8a3ePHePGQt8K9q/t8lmMhf7VvrWDPgf71sI1tW9txh3KdcCsqvWZwIbhykTEBGAa5Qmbazm2GepuQ+V28cMAmbmC8nPIf1R4xE83mr9lJ30Ow8rMDZXf9wEl4IhGBrcbampHRJwAfBQ4LTOf3J1jm2A0bWiXz2J3/5ZXAtuv+HbU51BlRxs66LtpKjAPKEXE/cBRwOLK4AHt8jm0wljoW6kxFvvX4o2F/tW+tYM+hyr2rcVobt+axb8UOoHyy81zeOql0BcOKvN+dn7h/qrK8gvZ+aXQ+2jNy7mjaUPX9pgpvxi7/v+2d78hdlxlHMe/vybFZNNqGyJipFCCERGkRluR+vdFRasGKgahStGm+EKU0tQiaBRKQVBfWEQopY10KQTBoNJQUWptVouSliSmbIJttTb1Tan/CjYxom0fX8zZeN3sbu7ejbv35n4/cLhnJufMnGfn3nkyM3fuAOuHMYaetpOc/qMBT9PdqH5xq49aDBcDr2j1DcDvmOPm5GGJgy4JPAVsnjV/ZLbFAjEMxbboM4bNPfWtwIFWH6V903wxjNy+qbWf4r8/HDAU22ElSp/bfqhz61mIY+Tew5hfV/q9ZG4djhjMrUMQw6z2Uywhty5XUB8CnmwfgJ1t3m10Z1YA1gB76G76fBTY1NN3Z+v3BHD1cm+QpcYAfAw42jbMIWDrEMdwBd1ZiRPAX4GjPX23t9h+D1w/ajEAVwLTbTtMAzesVAx9xvEg8BxwuJW9I7gt5oxhmLZFHzF8u31+DwP76NkZj9C+ac4YRmnfNKvtFC3pDdN2GMa/GyOQW5cSxyi9hzG/DksM5tbhiMHcOgQxzGo7xRJya1onSZIkSZIWZTnuoZQkSZIknYM8oJQkSZIkDcQDSkmSJEnSQDyglCRJkiQNxANKSZIkSdJAPKDU2EjyUpLDSY4k2ZNk4v+4rtvag4dJctNi15XOQ0le2aZvTPLbJLvPwtg+nWRjz/SuJG8acFmfT3L9UsckSRpd5tdTyza/aiz52BCNjSTHq+qCVt8NHKyqb/XRL3SflZcHXO8xumf7/GURfT4MXFVVO9r043TPAXp6VrvVVfXiIsczBdxSVQcW02+eZU0Av6qqLUtdliRpNJlfT/WZwvyqMeQVSo2rh4HXAyS5uZ1VPZLkpjbv0nbG8g66B9NekuTaJNOt3Tdau1VJJtu86SQzCWoyybYkNwIbgX1J9iW5IcntM4NI8pkkcyXdTwL3tTZ3ApuAvUl2JLk1yV1JHgDubWN9OMmhVq7sWf4X27geS/L1JNuAy4Hd7Wzy2iRTSS5v7U+Lsc0/nuRrbTn7k7wGoKr+ARxL8vazslUkSaPO/Gp+1bipKotlLApwvL2upksmnwXeBkwD64ALgKPAFuBS4GXgHa3PRuCPwKtb/4eAa1r/n/Ws46L2Oglsa/VjwIZWXwc8BZzfpn8NvHmOsT4DXNgz3buMW4GDwNo2PQGsafXNwIFWv7otf6JNr2+vU3RndOmdni/G1qaAra3+TeArPf13Al9Y6e1rsVgslpUp5lfzq2W8i1coNU7WJjkMHKDbsX8XeBfwo6o6UVXHgR8C727tn6mq/a1+BTBVVX+u7iswu4H3AH8ANiX5TpIPAn9faABVdYIukXwkyRvpEt/0HE3XV9ULCyxqb1WdbPXzgbuTTAN7gJn7Na4C7qnuLCdV9beFxrZAjAD/Au5v9YN0/yGY8Se6ZClJGk/m14WZX3VOW73SA5CW0cmqekvvjHb/xnxO9Dadq0FVPZ/kMuADwOeAjwPbzzCOXcCXgceBe+Zp82KS82r++0p6x7YDeA64jO5r7P/sGfNibpJe6G/x76qaWdZL/O++Yw1w8vQukqQxYX5dmPlV5zSvUGrc/RK4JslEknXAR+nu/5jtEeC9STYkWQVcC/wiyQbgvKr6AfBV4K1z9H0BuHBmoqoeAS4BPgF8b55xPUF3X0c/XgU825LjdcCqNv8BYHvaL+AlWT/XeM4UYx/rfwNwpM+xSpLGg/n1DDH2sX7zq0aCB5Qaa1V1iO5+jEfpdvi7quo3c7R7FvgSsA94DDhUVfcBrwOm2ld9Jlub2e4CfpJkX8+879P9etvz8wztx8D7+gzjDuBTSfbTJZ8Tbcw/BfYCB9r4bmntJ4E7Z340oI8Yz+SdwIN9jlWSNAbMr+ZXjQ8fGyKtgCT3A7dX1c/n+ffXAvdW1fuXd2SLk2QLcHNVXbfSY5EkyfwqLT+vUErLKMlFSZ6ku99kzmQHp85m3p324OUhtoHuq0iSJK0Y86u0crxCKUmSJEkaiFcoJUmSJEkD8YBSkiRJkjQQDyglSZIkSQPxgFKSJEmSNBAPKCVJkiRJA/kPC87JWmvPvUMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(121)\n",
    "plt.hist(por1, facecolor='red',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label='Well 1')\n",
    "plt.hist(por2, facecolor='blue',bins=np.linspace(0.0,0.4,20),alpha=0.2,density=False,edgecolor='black',label = 'Well 2')\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,8.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Frequency'); plt.title('Porosity Well 1 and 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.hist(por1, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.2,density=True,cumulative=True,edgecolor='black',label='Well 1')\n",
    "plt.hist(por2, facecolor='blue',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.2,density=True,cumulative=True,edgecolor='black',label='Well 2')\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Well 1 and 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram and cumulative distribution functions indicate the porosity distributions in both wells 1 and 2. Ocular inspection does indicate: \n",
    "\n",
    "1. not a lot of samples are available for each (n=20) \n",
    "2. the distributions look different\n",
    "\n",
    "Let's use statistics to learn more about the subsurface from these two wells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Intervals\n",
    "\n",
    "Let's first demonstrate the calculation of the confidence interval for the sample mean at a 95% confidence level.  This could be interpreted as the interval over which there is a 95% confidence that it contains the true population mean.  We use the student's t distribution as we assume we do not know the variance and the sample size is small. \n",
    "\n",
    "\\begin{equation}\n",
    "x̅ \\pm t_{\\frac{\\alpha}{2},n-1} \\times \\frac {s}{\\sqrt{n}} \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The confidence interval for the Well 1 mean porosity is 0.164 +/- 0.013, with a range of 0.151, 0.178\n",
      "The confidence interval for the Well 2 mean porosity is 0.2 +/- 0.021, with a range of 0.179, 0.221\n"
     ]
    }
   ],
   "source": [
    "ci_95_por1 = st.t.interval(0.95, len(df)-1, loc=np.mean(por1), scale=st.sem(por1))\n",
    "print('The confidence interval for the Well 1 mean porosity is ' + str(round(np.mean(por1),3)) + ' +/- ' + str(round(ci_95_por1[1]-np.mean(por1),3))\n",
    "      + ', with a range of ' + str(round(ci_95_por1[0],3)) + ', ' + str(round(ci_95_por1[1],3)))\n",
    "\n",
    "ci_95_por2 = st.t.interval(0.95, len(df)-1, loc=np.mean(por2), scale=st.sem(por2))\n",
    "print('The confidence interval for the Well 2 mean porosity is ' + str(round(np.mean(por2),3)) + ' +/- ' + str(round(ci_95_por2[1]-np.mean(por2),3))\n",
    "      + ', with a range of ' + str(round(ci_95_por2[0],3)) + ', ' + str(round(ci_95_por2[1],3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using Confidence Intervals to Build an Uncertainty Model\n",
    "\n",
    "The confidence interval for the mean could be applied directly as an uncertainty model! The uncertainty in mean porosity at each well could be applied to calculate the uncertainty in the mean porosity applied to model porosity around the two wells. This method would account for the limited number of porosity samples available. If more porosity samples we available then this uncertainty would decrease (note the $\\sqrt{n}$ in the denominator of the confidence interval in the mean calculation above.\n",
    "\n",
    "One possible work flow would be to apply the affine correction from the *GeostatsPy* package to formulate a the P5 and P95 porosity distributions for porosity in each well. Note, these distribution scenarios are based only on uncertainty in the mean porosity and assumes the distribution tails are not anchored (are free to shift)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAGWCAYAAAAQQ1yjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X+cHXV56PHPQwIEDCwKEYPhRyjc3AABhMUfVUIsWiMarC1W8EejRGgRAZXVovT6C4qtsVdRqCjFK2oRwVsFelGxxRW9tsKCEYqBS0ADCRuIBCTRBDbkuX+cSXKy7Nk92T2ze3583q/XvvbMzHe+83zP2Z3nPHNm5kRmIkmSJElSWXaY6AAkSZIkSe3NwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT5UmIu6OiHnjuL0PR8Q/NbC/dRFxYPH4KxFxYQP7viwi/kej+msHEbFf8ZxPmuhYGi0ifh0Rryoefywivj7RMUlqXebXYfs2vw5iflWzsPBsU8U/4vpiR/NIRPyviJg6njFk5qGZ2VvEM6adQUT0RsSGiFgbEU9GxO0RcV5E7Fy1vYsy81119jViu8ycmpkPjDbmqu29IyJ+Mqjvv8rMC8ba9xDbelYCj4gDIiIjYnKjtzdCLM8a93Ay88HiOX+mjr5LHVNEfD8iPlg1/cJie0PNe0EDt7tTRHyr+P/N8XxjKak+5tcR+zK/lsz8OqrtvjQifhARayJidURcGxHTG9W/6mPh2d4WZOZU4CjgGOBvtreD8d6ZjuA9mbkbMB04FzgZuDEiopEbabIxt6Q2eA5vAY6rmp4L3DPEvPsyc1WDt/0T4G1Ao/uV1Djm11FosjG3pDZ4Dicqvz4X+BJwALA/sBb4Xw3sX3Ww8OwAmbkS+C5wGEBE7BMR1xdHfZZFxGmb2xZHTr8VEV+PiCeBd0TEzhHx2Yh4uPj57OYjoRGxV0T8a0Q8UfT344jYoVj264h4VUTMBz4MvLk4QvyLiHhTRNxeHWdEnBsR36ljPL8rjvSeCLwMeF1V7F8vHk8pxvBYEdttEbF3RPwtcCxwSRHLJUX7jIgzI+I+4L6qeQdVbXqv4mjZ2oj4UUTsX7R71tHBzUd9I2I2cBnwsmJ7TxTLtzlyGhGnFa/FmuK12adqWUbEX0XEfRHxeERcOpY3A8Xr0hMRd0bEbyPimxExpWr5GyJiSXHk+/7i9SMiuiLiiojoj4iVEXFhFKftFEdf/29EfCYi1gDfrDHu10XEz4u+H4qIj1Vtd5vnsXgOLyj6XRsRN0XEXkXzW4rfTxT9H1c8d3Oq+nt+VD6VmDaKp+kW4OWb/5ap/M18FugeNG9zHETE64vn7YmI+GlEHL69G83MpzPzs5n5E2DEI9OSJpb51fxazfxal4nKr9/NzGsz88nM/D1wCfDyUcSvMbDw7AARsS9wAvDzYtY3gBXAPsBJwEURcXzVKm8AvgXsAfwzcD7wUuBI4AjgxWw9untu0dc0YG8qCTCrt5+Z3wMuAr5ZnOpxBHA9MLNIHJu9DfhavePKzAeBPio7qMEWAl3AvsCewF8B6zPzfODHVI7uTs3M91St8yfAS4BDamzyrcAFwF7AEirPzUgxLi22/R/F9vYY3CYi/gj4JPDnVI42LweuHtTs9VSOqh9RtHvNSNsewZ8D84GZwOHAO4pYXgx8FfgAldd/LvDrYp0rgY3AQcCLgD8Gqk+pegnwAPB8Kq/lUOP+HfAXRd+vA86IiD8ZJs63AO8s+twJ6Cnmzy1+71H0/yMqz9nbqtY9Bfi3zFw94rPxbLcCO1N5vjdv7wfAskHzbgGIiKOALwN/SeXv7YvA9VF1qpqk9mN+Nb8Owfw6vGbJr3OBu8fYh7aThWd7+05xJOwnwI+oJMB9gVcAf52ZGzJzCfBPwNur1vuPzPxOZm7KzPVUEsInMvPRYifz8ar2A1R25vtn5kBm/jgzt0mMQ8nMp6gctXsbQEQcSuX0h3/dzjE+DDxviPkDVHZQB2XmM5l5e2Y+OUJfn8zMNcWYh/J/MvOWIvbzqRxt3Hc74x3KW4EvZ+YdRd8fKvo+oKrN32XmE8WbgR9SeZMyFp/LzIczcw1wQ1V/i4pYflC8/isz856I2Bt4LfDe4oj4o8BnqJyOtdnDmfn5zNxY6znMzN7MvKvo+04qb9KOG6pt4X9l5v8r+ruG4cd9JfCWqiOmb2c73mgNivMp4GfA3Ih4HpUE/ACVN1Wb5x1C5f8K4DTgi5n5s+Lv7UrgKSpvKCW1H/Or+bUW8+swmiG/Fp+YfoTKQQCNIwvP9vYnmblHZu6fme8udi77AGsyc21Vu+XAC6umHxrUzz5Fm+r2m09VWUzlKNVNEfFARJy3HfFt3pEFlZ3YNcUOaXu8EFgzxPyvAd8Hro7K6UufiogdR+hr8LhrLs/MdcV296ndvG7bPL9F34+x7WtSfZ3D74FaN7LYCAwe547ApuJnpP72Be4fot/9i376i1NdnqBy1PH5VW1Gev6IiJdExA+jcmH/b6kctd1rmFXqHTeZ+TMqR3yPi4j/TuXI8fU14ri7OIVoXUQMdUQfKkdb51I54r/5Jg4/qZr3UGZuft32B87d/NwUz8++NObvQ1LzMb+aX8H8OlQcTZ1fo3KK93eBczLzx6PpQ6Nn4dl5HgaeFxG7Vc3bD1hZNT34iOrDVP7xq9s/DJCZazPz3Mw8EFgAvH/QaUW1+iQz/xN4mspO5i1s59Gz4mjo0VSOkg3ueyAzP56ZhwB/SOVUmr+oFcsI8zfbcvQ1KncwfB6V5+F3xexdq9pW34ltpH63eX4j4jlUjiavrLlGbQ9SObJdbSaVnfimZzd/loeAP6gx/ylgr+LN1h6ZuXtmHlrVZvA4hxr3VVSS1b6Z2UXlOpXRXE9T6zm9kspR/rcD38rMDUOuXLkj5NTip1biuYXK3+Zctv6N/V8q14RsOQ2o8BDwt1XPzR6ZuWtmfmP7hiWphZlfza/DMb9uNSH5NSrXDv8bcEFmjuoTW42NhWeHycyHgJ8Cn4zKDQIOp3L6x3DXU3wD+JuImFZcfP4RYPNNBl4fEQcVR1WfpHJDlKFuivIIcEDVaRqbfZXKBd4bs3JDlRFFxK4RcRxwHZVrBW4cos0rI2JOVC7Of5LKqUGb43oEOLCebQ1yQkS8IiJ2onItys8y86Hi9KiVwNsiYlJEnMq2yeURYEax3lCuAt4ZEUcW1yxcVPT961HE+L+B10XEHxex7EPleqHB17TUckURy/ERsUNUbmn+3zOzH7gJ+IeI2L1Y9gfF61DLUOPejconAhuK613esv1DBGA1lSPMg1/HrwFvpJIcvzrKvjf7KZVrZd5GkRgz8/Fi229j28R4OfBXxRHniIjnROVGD7sN7nQkUbnZyOabUexU/J829M6SkhrP/LolFvPr0MyvW417fo2IFwI3A5dm5mVjjF+jZOHZmU6hctTuYeDbwEcz8wfDtL+Qyk0G7gTuAu4o5gEcTOXo0TrgP4B/zOK7xQa5tvj9WETcUTX/a1TuBljPkadLImItlR3uZ6kkgfk1jjS+gMoNHJ4EllK5VmDz95xdDJwUlTvYfa6O7W52FfBRKqcAHU3l2pHNTqNyrcBjwKFUdqqb3UzlAvZVEfGbwZ1m5r8D/6MYTz+VpHry4Hb1yMy7qby+nyzi/A8q11J8vM71b6Vys4HPAL+l8rxtPlr8F1RuQPBL4HEqz+9w34E11LjfDXyieB0/QuW6ku2WlTvS/S3wf4tTb15azF9B5e8zGeJI/Si2cTuVmyD8V9WiH1M5BeqWqrZ9VP4GLqHy3CyjuKHEKNwLrKdyKtj3i8f7D7uGpGZhfjW/1lrf/LrtNsY7v76LSjH90dh6KvC6UQ1AoxY58nXqUmkiYhfgUeCozLxvouNR64uIL1O5EcN2f6+eJLUL86sazfyqsWr1L6FV6zsDuM2kqEaIyp0K/5TK7eglqZOZX9Uw5lc1Qmmn2kbElyPi0Yj4rxrLIyI+F5Uv9b0zKt/Tow4SEb8GzqHyXWXSmETEBVRO2Vmcmb+a6HiksphfNRLzqxrJ/KpGKe1U24iYS+W6hK9m5mFDLD8BOIvKFy+/BLg4M19SSjCSJLUJ86skqRWV9olnZt7C0N//tNkbqCTNLG77vUdEDHchtSRJHc/8KklqRRN5jecL2fYLcVcU8/oHN4yI04HTAaZMmXL0fvvtNy4BlmXTpk3ssENr31DYMUyMp57aQPVZCpMnT2bjxo1bpsd8BsME3Wts8o47snFgYGI23iCOoTk8tGLFbzJz2kTHMcHqyq/m1ubUDuNwDM1hPMcw+P3J9qq17uTJO7JxYwPy0gTeS9XcutVEFp5DfS/dkH8Wmfkl4EsAs2bNynvvvbfMuErX29vLvHnzJjqMMXEME6O7e3/6+p6/Zbq3dxHz5l2xdfmcu+hbPPp/6+6F6+nrGv+vjOw9+1PM+9wHx327jeQYmkPA8omOoQnUlV/Nrc2pHcbhGJrDeI5h8PuT7V6/xvuX3qcvZN5OY7+J7kS9vwFza7WJPJSzAti3anoGle+9kiRJo2d+lSQ1nYksPK8H/qK4+95Lgd9m5rNOs5UkSdvF/CpJajqlnWobEd8A5gF7RcQK4KPAjgCZeRlwI5U77i0Dfg+8s6xYJElqF+ZXSVIrKq3wzMxTRliewJmN2NbAwAArVqxgw4YNjeiudF1dXSxdunSiw9jGlClTmDFjBjvuuONEhyJJGsZ45Vdza2OYXyWpYiJvLtQwK1asYLfdduOAAw4gYmIuHN4ea9euZbfddpvoMLbITB577DFWrFjBzJkzJzocSVITMLeOnflVkrZq7ftEFzZs2MCee+7ZEomxGUUEe+65Z8sc1ZYklc/cOnbmV0naqi0KT8DEOEY+f5KkwcwNY+dzKEkVbVN4SpIkSZKaU1tc4znYggVz6e9v3HeIT5++PzfccMuwbSZNmsScOXPYuHEjs2fP5sorr2TXXXfle9/7Hueccw7PPPMM73rXuzjvvPMAeMc73sGPfvQjurq6APjKV77CkUce2bCYJUlqJHOrJGks2rLw7O9fTl/ftIb11909cqLdZZddWLJkCQBvfetbueyyyzjnnHM488wz+cEPfsCMGTM45phjOPHEE9l338r3ei9evJiTTjqpYXFKklQWc6skaSw81bYExx57LMuWLePWW2/loIMO4sADD2SnnXbi5JNP5rrrrpvo8CRJajnmVklqbRaeDbZx40a++93vMmfOHFauXLnlCCzAjBkzWLly5Zbp888/n8MPP5z3ve99PPXUUxMRriRJTc/cKkmtz8KzQdavX8+RRx5Jd3c3++23H4sWLaLyHd7b2nx3u09+8pPcc8893HbbbaxZs4a///u/H++QJUlqauZWSWofbXmN50Sovg5lsxkzZvDQQw9tmV6xYgX77LMPANOnTwdg55135p3vfCef/vSnxy9YSZJagLlVktqHn3iW6JhjjuG+++7jV7/6FU8//TRXX301J554IgD9/f0AZCbf+c53OOywwyYyVEmSWoK5VZJaU1t+4jl9+v513S1ve/objcmTJ3PJJZfwmte8hmeeeYZTTz2VQw89lLVr1/LWt76V1atXk5kceeSRXHbZZQ2LV81rrF9HsGbVI3TPWbVletEZ6+k5664t09O7nn0K2jbb79lA/+rabaZnAn7ZuaRnM7dKE2c07x8WLeqhp2chAKsefITctLGM0ADY+OQzvHCnB0e9/qQdoHvhwLPmL/rwJnouWj+W0ADf3zSLtiw8R/pesDKsW7duyPknnHACJ5xwwrPm33zzzWWHpCY01q8j6J6zir7FW/9te5+ObaZH3P7qpK9ruB2vO2VJQzO3ShNnNO8fensnb1mn8v5h5zJCq/S/cP0I7y9Gp3cSDerX9zfNwFNtJUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJWqLb9OZcH8ufSvbOB3jb1wf2743vC3kZ80aRJz5szZMn3yySdz3nnnNSyGoTzxxBNcddVVvPvd796u9T72sY8xdepUenp6SopMktRuzK3DM7dK0vDasvDsX7mcvsWj/67Ewbo/MHKi3WWXXViyZEnDtlmPJ554gn/8x3/c7uQoSdL2MrdKksbCU21L9Nvf/pZZs2Zx7733AnDKKadw+eWXAzB16lTOPfdcjjrqKI4//nhWr14NwP3338/8+fM5+uijOfbYY7nnnnsAeOSRR3jjG9/IEUccwRFHHMFPf/pTzjvvPO6//36OPPJIPvCBDwCwePFijjnmGA4//HA++tGPbonlb//2b5k1axavetWrtsQjSVKrMbdKUmuy8GyQ9evXc+SRR275+eY3v0lXVxeXXHIJ73jHO7j66qt5/PHHOe200wD43e9+x1FHHcUdd9zBcccdx8c//nEATj/9dD7/+c9z++238+lPf3rLEdezzz6b4447jl/84hfccccdHHroofzd3/0df/AHf8CSJUtYvHgxN910E/fddx+33norS5Ys4fbbb+eWW27h9ttv5+qrr+bnP/85//Iv/8Jtt902Yc+TJEn1MrdKUvtoy1NtJ0Kt04Fe/epXc+2113LmmWfyi1/8Ysv8HXbYgTe/+c0AvO1tb+NP//RPWbduHT/96U9505vetKXdU089BcDNN9/MV7/6VaByzUtXVxePP/74Ntu66aabuOmmm3jRi14EwLp167jvvvtYu3Ytb3zjG9l1110BOPHEExs4ckmSymFulaT2YeFZsk2bNrF06VJ22WUX1qxZw4wZM4ZsFxFs2rSJPfbYY9TXs2QmH/rQh/jLv/zLbeZ/9rOfJSJG1ackSc3G3CpJrcdTbUv2mc98htmzZ/ONb3yDU089lYGBAaCSNL/1rW8BcNVVV/GKV7yC3XffnZkzZ3LttdcClWS3+Uju8ccfzxe+8AUAnnnmGZ588kl222031q5du2Vbr3nNa/jyl7/MunXrAFi5ciWPPvooc+fO5dvf/jbr169n7dq13HDDDeM2fkmSGs3cKkmtpy0/8Zz+wv3rulve9vQ3ks3XoWw2f/58Tj31VP7pn/6JW2+9ld122425c+dy4YUX0tPTw3Oe8xzuvvtujj76aLq6uvjmN78JwD//8z9zxhlncOGFFzIwMMDJJ5/MEUccwcUXX8zpp5/OFVdcwaRJk/jCF77Ay172Ml7+8pdz2GGH8drXvpbFixezdOlSXvaylwGVmyx8/etf56ijjuLNb34zRx55JPvvvz/HHntsw54bSVJnMLeaWyVpLNqy8Bzpe8HK8Mwzzww5f+nSpVse/8//+T8BthxJveCCC7jgggu2aT9z5ky+973vPaufvffem+uuu+5Z86+66qptps855xzOOeecZ7U7//zzOf/880cYhSRJQzO3mlslaSw81VaSJEmSVKq2/MSzFWy+VkTtZcGCufT31z4Vbfr0Edaffw/9KwdqLl/z4ADdC7cuX/ThTfRctL7u+KZnAt4MQ1J7MreqVZXx/mHRGevpOeuuyvpdObb4ejbQv7p2H76/UD0sPKUG6u9fTl/ftNGvv3KAvsW1/y27Fw7Q17V1x947iW2mR2ZSkCSp2ZTx/qH36Rj2PcV29b86R3i/4fsLjcxTbSVJkiRJpbLwlCRJkiSVysJTkiRJklSqtrzGc8HcufQvb+B3je2/PzfcMvxt5CdNmsScOXPYuHEjs2fP5sorr2TXXXfl4osv5vLLLyczOe2003jve98LwMc+9jEuv/xypk2rnM9/0UUXccIJJzQsZkmSGsncKkkai7YsPPuXL6dv2ugv0B6su45Eu8suu7BkyRIA3vrWt3LZZZfxx3/8x1x++eXceuut7LTTTsyfP5/Xve51vOAFLwDgfe97Hz09PQ2LU5KksphbJUlj4am2JTj22GNZtmwZS5cu5aUvfSm77rorkydP5rjjjuPb3/72RIcnSVLLMbdKUmuz8GywjRs38t3vfpc5c+Zw2GGHccstt/DYY4/x+9//nhtvvJGHHnpoS9tLLrmEww8/nFNPPZXHH398AqOWJKl5mVslqfVZeDbI+vXrOfLII+nu7ma//fZj0aJFzJ49m7/+67/m1a9+NfPnz+eII45g8uTK2c1nnHEG999/P0uWLGH69Omce+65EzwCSZKai7lVktpHW17jORGqr0OptmjRIhYtWgTAhz/8YWbMmAHA3nvvvaXNaaedxutf//rxCVSSpBZhbpWk9uEnniV79NFHAXjwwQf5l3/5F0455RQA+vv7t7T59re/zWGHHTYh8UmS1GrMrZLUetryE8/p++9f193ytqe/0fqzP/szHnvsMXbccUcuvfRSnvvc57J27Vo++MEPsmTJEiKCAw44gC9+8YsNi1eSpEYzt0qSxqItC8+RvhesDOvWrRty/o9//OMh53/ta18rMxxJkhrK3CpJGgtPtZUkSZIklcrCU5IkSZJUqrYpPDNzokNoaT5/kqTBzA1j53MoSRVtcY3nlClTeOyxx9hzzz2JiIkOp+VkJo899hhTpkyZ6FAkSU3C3Dp25tf2sWDBXPr7G3dzrcHWrHqE7jmrRr/+gwN0LxzYZt6iD2+i56L1Yw0NgOmZgPsBjU1bFJ4zZsxgxYoVrF69eqJDqcuGDRuaLglNmTJly/egSZJkbm0M82t76O9fTl/ftNL6756zir7Fo39b3r1wgL6ubQvD3kk8a97oWXRq7Nqi8Nxxxx2ZOXPmRIdRt97eXl70ohdNdBiSJNVkbpUkNVLbXOMpSZIkSWpOFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSmXhKUmSJEkqlYWnJEmSJKlUFp6SJEmSpFJZeEqSJEmSSlVq4RkR8yPi3ohYFhHnDbF8v4j4YUT8PCLujIgTyoxHkqRWZ26VJLWi0grPiJgEXAq8FjgEOCUiDhnU7G+AazLzRcDJwD+WFY8kSa3O3CpJalVlfuL5YmBZZj6QmU8DVwNvGNQmgd2Lx13AwyXGI0lSqzO3SpJaUmRmOR1HnATMz8x3FdNvB16Sme+pajMduAl4LvAc4FWZefsQfZ0OnA4wbdq0o6+55ppSYh4v69atY+rUqRMdxpg4hqEtXXoXs2dPHv36v1zP7BlRe/mvNjF70tbpdXvPYOojK0a9vWbRDuNwDM3hlWf13J6Z3RMdR1nMrbW1Q16C9hhHO45hrPl9JCPl/xHXH/T+ANpjn+4YmkOjcmuZheebgNcMSo4vzsyzqtq8v4jhHyLiZcAVwGGZualWv7Nmzcp77723lJjHS29vL/PmzZvoMMbEMQytu3t/+vqmjX79OXfRt7h2YuteuJ6+rq2JqffsTzHvcx8c9faaRTuMwzE0h7hvU7sXnubWGtohL0F7jKMdxzDW/D6SkfL/iOsPen8A7bFPdwzNoVG5tcxTbVcA+1ZNz+DZp/ssAq4ByMz/AKYAe5UYkyRJrczcKklqSWUWnrcBB0fEzIjYicoNDq4f1OZB4HiAiJhNJTmuLjEmSZJamblVktSSSis8M3Mj8B7g+8BSKnfYuzsiPhERJxbNzgVOi4hfAN8A3pFlnfsrSVKLM7dKklpVeVdJA5l5I3DjoHkfqXr8S+DlZcYgSVI7MbdKklpRmafaSpIkSZJk4SlJkiRJKpeFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkko1eaIDkFrJggVz6e9fXnP59OkjrD//HvpXDtRcvubBAboX1l4+PROIkcKUJEnbYXB+X7Soh56ehVumx5rfRzK9K4fvv2cD/atrt/H9gVqBhae0Hfr7l9PXN230668coG9x7X+77oUD9HUNlzhMKpIkNdrg/N7bO3m78v1I+X2s+len7w/U8jzVVpIkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUqhELz4joi4gzI+K54xGQJEntztwqSeo09XzieTKwD3BbRFwdEa+JiCg5LkmS2pm5VZLUUSaP1CAzlwHnR8T/AF4PfBnYFBFfBi7OzDUlxyiNmwUL5tLfv7zm8jWrHqF7zqpR97/mwQG6Fw7UXD49E/C9p9TuzK3S+Bopv0+fPsL68++hf+Uw+bsrRxtapf+eDfSvrt2H7w/UDkYsPAEi4nDgncAJwP8G/hl4BXAzcGRp0UnjrL9/OX1902ou756zir7Fdf3bDL3+wgH6uoZLHCYVqVOYW6XxM1J+H3H9lQNjyv8j9r86fX+gtjfif1BE3A48AVwBnJeZTxWLfhYRLy8zOEmS2pG5VZLUaeo5dPOmzHygekZEzMzMX2Xmn5YUlyRJ7czcKknqKPXcXOhbdc6TJEn1MbdKkjpKzU88I+K/A4cCXRFRffR1d2BK2YFJktRuzK2SpE413Km2s6jcaW8PYEHV/LXAaWUGJUlSmzK3SpI6Us3CMzOvA66LiJdl5n+MY0ySJLUlc6skqVMNd6rtBzPzU8BbIuKUwcsz8+xSI5Mkqc2YWyVJnWq4U22XFr/7xiMQSZI6gLlVktSRhjvV9obi95XjF44kSe3L3CpJ6lTDnWp7A5C1lmfmiaVEJElSmzK3SpI61XCn2n563KKQJKkzmFslSR1puFNtfzSegUiS1O7MrZKkTrVDrQURcU3x+66IuLPq566IuLOeziNifkTcGxHLIuK8Gm3+PCJ+GRF3R8RVoxuGJEnNz9wqSepUw51qe07x+/Wj6TgiJgGXAq8aN77SAAAgAElEQVQGVgC3RcT1mfnLqjYHAx8CXp6Zj0fE80ezLUmSWoS5VZLUkWp+4pmZ/cXv5cBTwBHA4cBTxbyRvBhYlpkPZObTwNXAGwa1OQ24NDMfL7b16PYPQZKk1mBulSR1qsiseXO9SoOIdwEfAW4GAjgO+ERmfnmE9U4C5mfmu4rptwMvycz3VLX5DvD/gJcDk4CPZeb3hujrdOB0gGnTph19zTXX1D3AZrRu3TqmTp060WGMSbuOYenSu5g9u/aJAEt/uZ7ZM2LU21z6q03MnjTq1Z9l3d4zmPrIisZ1OEHaYRyOoTm88qye2zOze6LjGIm5tfHaIS9Be4yjGccwUn4fbN26vZg69Tdb1x9j/h9Jo98fQHvs0x1Dc2hUbq2n8LwX+MPMfKyY3hP4aWbOGmG9NwGvGZQcX5yZZ1W1+VdgAPhzYAbwY+CwzHyiVr+zZs3Ke++9t56xNa3e3l7mzZs30WGMSbuOobt7f/r6ptVcp3vOXfQtrj9xPWv9hevp62pc4uo9+1PM+9wHG9bfRGmHcTiG5hD3bWqVwtPc2mDtkJegPcbRjGMYKb8P1tu7iHnzrti6/hjz/0ga/f4A2mOf7hiaQ6Nya81TbausANZWTa8FHqpzvX2rpmcADw/R5rrMHMjMXwH3AgfX0bckSa3M3CpJ6ig1D91ExPuLhyuBn0XEdVS+9PoNwK119H0bcHBEzCz6OBl4y6A23wFOAb4SEXsB/w14YLtGIElSizC3SpI61XDnDOxW/L6/+Nnsuno6zsyNEfEe4PtUrjH5cmbeHRGfAPoy8/pi2R9HxC+BZ4APbD7tSJKkNmRulSR1pJqFZ2Z+fKydZ+aNwI2D5n2k6nEC7y9+JElqa+ZWSVKnGvEq6YiYBnwQOBSYsnl+Zv5RiXFJktS2zK2SpE5Tz82F/hm4B5gJfBz4NZVrTCRJ0uiYWyVJHaWewnPPzLwCGMjMH2XmqcBLS45LkqR2Zm6VJHWUer6QaKD43R8Rr6Ny2/YZ5YUkSVLbM7dKkjpKPYXnhRHRBZwLfB7YHXhfqVFJktTezK2SpI4yYuGZmf9aPPwt8Mpyw5Ekqf2ZWyVJnWbEazwj4sCIuCEifhMRj0bEdRFx4HgEJ0lSOzK3SpI6TT03F7oKuAZ4AbAPcC3wjTKDkiSpzZlbJUkdpZ7CMzLza5m5sfj5OpBlByZJUhszt0qSOkrNazwj4nnFwx9GxHnA1VSS4puB/zMOsUmS1FbMrZKkTjXczYVup5IMo5j+y6plCVxQVlCSJLUpc6skqSPVLDwzc+Z4BiJJUrszt0qSOtWIX6cSETsCZwBzi1m9wBczc6DmSpIkqSZzq7R9FiyYS3//8lGvP336CP3Pv4f+lVv//RadsZ6es+7aun7X2C7BXtCzgf7VtfuYntUnQkjtacTCE/gCsCPwj8X024t57yorKEmS2py5VdoO/f3L6eubVl7/KwfoW7z1bXHv07HN9Jj7X530dQ1XWFp0qv3V8x91TGYeUTV9c0T8oqyAJEnqAOZWSVJHqefrVJ6JiD/YPFF8wfUz5YUkSVLbM7dKkjpKPZ94foDKbd8foHIewP7AO0uNSpKk9mZulSR1lGELz4jYAVgPHAzMopIc78nMp8YhNkmS2o65VZLUiYYtPDNzU0T8Q2a+DLhznGKSJKltmVslSZ2onms8b4qIP4sIb7clSVJjmFslSR2lnms83w88B9gYERuonBKUmbl7qZFJktS+zK2SpI4yYuGZmbuNRyCSJHUKc6skqdPUPNU2Ip4fEZ+NiH+NiIsiwqOwkiSNgblVktSphrvG86vA74DPA7sBnxuXiCRJal/mVklSRxruVNsXZOb5xePvR8Qd4xGQJEltzNwqSepIwxWeERHPpXLDA4BJ1dOZuabs4CRJajPmVklSRxqu8OwCbmdrcgTYfGQ2gQPLCkqSpDZlbpUkdaSahWdmHjCOcUiS1PbMrZKkTjXczYUkSZIkSRozC09JkiRJUqksPCVJkiRJpaqr8IyIV0TEO4vH0yJiZrlhSZLU3sytkqROMmLhGREfBf4a+FAxa0fg62UGJUlSOzO3SpI6TT2feL4ROBH4HUBmPgzsVmZQkiS1OXOrJKmj1FN4Pp2ZSeX7xYiI55QbkiRJbc/cKknqKPUUntdExBeBPSLiNODfgMvLDUuSpLZmbpUkdZTJIzXIzE9HxKuBJ4FZwEcy8welRyZJUpsyt0qSOs2IhWdEvA+41oQoSVJjmFslSZ2mnlNtdwe+HxE/jogzI2LvsoOSJKnNmVslSR1lxMIzMz+emYcCZwL7AD+KiH8rPTJJktqUuVWS1Gnq+cRzs0eBVcBjwPPLCUeSpI5ibpUkdYQRC8+IOCMieoF/B/YCTsvMw8sOTJKkdmVulSR1mhFvLgTsD7w3M5eUHYwkSR3C3CpJ6ig1C8+I2D0znwQ+VUw/r3p5Zq4pOTZJktqKuVWS1KmG+8TzKuD1wO1AAlG1LIEDS4xLkqR2ZG6VJHWkmoVnZr6++D1z/MKRJKl9mVslSZ1qxGs8I+LfM/P4keZJrWDBgrn09y8HYNGiHnp6Fm6zfM2qR+ies6rm+tO7cvj+ezbQv7p2m+k5+AMOSZ3I3CpJ6jTDXeM5BdgV2CsinsvWd8u7U/nOManl9Pcvp69vGgC9vZO3PN6se84q+hbXc8+tGv2vTvq6hissLTqlTmZulSR1quHeYf8l8F4qifB2tibHJ4FLS45LkqR2ZG6VJHWk4a7xvBi4OCLOyszPj2NMkiS1JXOrJKlTjXhOYWZ+PiIOAw4BplTN/2qZgUmS1K7MrZKkTlPPzYU+CsyjkhxvBF4L/AQwOUqSNArmVklSp9mhjjYnAccDqzLzncARwM6lRiVJUnszt0qSOko9hef6zNwEbIyI3YFH8QuuJUkaC3OrJKmj1PO9EX0RsQdwOZU78K0Dbi01KkmS2pu5VZLUUeq5udC7i4eXRcT3gN0z885yw5IkqX2ZWyVJnaZm4RkRRw23LDPvKCckSZLak7lVktSphvvE8x+GWZbAHzU4FkmS2p25VZLUkWoWnpn5yvEMRJKkdmdulSR1qnq+x/Mvhprvl1xLkjQ65lZJUqep5662x1Q9nkLle8fuwC+5liRptMytkqSOUs9dbc+qno6ILuBr9XQeEfOBi4FJwD9l5t/VaHcScC1wTGb21dO3JEmtytwqSeo0O4xind8DB4/UKCImAZcCrwUOAU6JiEOGaLcbcDbws1HEIklSOzC3SpLaWj3XeN5A5U57UClUDwGuqaPvFwPLMvOBop+rgTcAvxzU7gLgU0BPnTFLktTSzK2SpE4TmTl8g4jjqiY3Asszc8WIHVdO8Zmfme8qpt8OvCQz31PV5kXA32Tmn0VEL9Az1OlAEXE6cDrAtGnTjr7mmnpyc/Nat24dU6dOnegwxqRVx7B06V3Mnl053rJu3V5MnfqbbZf/cj2zZ8To+//VJmZPGlOI22Xd3jOY+siI/45Nrx3G4RiawyvP6rk9M7snOo6RmFsbr1Xz0mDtMI7RjGHZsnsZGHi65vIdd4SDDqrn1iQ1+r9vAwMDtd/z7jgJDpq+Nf+vyxcyNVZuXX/FJgYGRr15dgQOGsf3B9Ae+3TH0BwalVvrucbzRwARsfvm9hHxvMxcM8KqQ7173/IfHxE7AJ8B3lFHDF8CvgQwa9asnDdv3kirNLXe3l4cw8To6VlIX980AHp7FzFv3hXbLj/rLvoWjz6x9Vy0nr6u0Reu26v37E8x73MfHLftlaUdxuEYtD3MrY3XqnlpsHYYx2jGUJ2fy7C9+b336QuZt9PfbF3/4+Ob3xuhHfbpjqG91HOq7elUTtlZD2yikvQSOHCEVVcA+1ZNzwAerpreDTgM6I0IgBcA10fEid4EQZLUzsytkqROU8+hnw8Ah2bmb0Zsua3bgIMjYiawEjgZeMvmhZn5W2CvzdPDnQ4kSVKbMbdKkjpKPXe1vZ/K3fa2S2ZuBN4DfB9YClyTmXdHxCci4sTt7U+SpDZibpUkdZR6PvH8EPDTiPgZ8NTmmZl59kgrZuaNwI2D5n2kRtt5dcQiSVI7MLdKkjpKPYXnF4GbgbuoXIciSZLGxtwqSeoo9RSeGzPz/aVHIklS5zC3SpI6Sj3XeP4wIk6PiOkR8bzNP6VHJklS+zK3SpI6Sj2feG6+W96HqubVc8t3SZI0NHOrJKmjjFh4ZubM8QhEkqROYW6VJHWaEQvPiPiLoeZn5lcbH44kSe3P3CpJ6jT1nGp7TNXjKcDxwB2AyVGSpNExt0qSOko9p9qeVT0dEV3A10qLSJKkNmdulSR1mnruajvY74GDGx2IJEkdzNwqSWpr9VzjeQOVO+1BpVA9BLimzKAkSWpn5lZJUqep5xrPT1c93ggsz8wVJcUjSVInMLdKkjpKzcIzIg4C9s7MHw2af2xE7JyZ95cenTTIggVz6e9fPur116x6hO45qwBYdMZ6es66a9vlDw7QvXBg1P1PzwRi1OtLam/mVklSpxruE8/PAh8eYv76YtmCUiKShtHfv5y+vmmjXr97zir6Flf+7Hufji2PtyxfOEBf11gKR4tOScMyt0qSOtJwNxc6IDPvHDwzM/uAA0qLSJKk9mVulSR1pOEKzynDLNul0YFIktQBzK2SpI40XOF5W0ScNnhmRCwCbi8vJEmS2pa5VZLUkYa7xvO9wLcj4q1sTYbdwE7AG8sOTJKkNmRulSR1pJqFZ2Y+AvxhRLwSOKyY/X8y8+ZxiUySpDZjbpUkdaoRv8czM38I/HAcYpEkqSOYWyVJnWa4azwlSZIkSRozC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklSqyRMdgCRJklrbggVz6e9fDsCiRT309CzcrvXXrHqE7jmrygit0v+DA3QvHKi7/aIPb6LnovVbpqdnAlFCZFLnsPCUJEnSmPT3L6evbxoAvb2TtzyuV/ecVfQtLu9taffCAfq66i8ceycxqL1FpzRWnmorSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEpSZIkSSqVhackSZIkqVQWnpIkSZKkUk2e6AAkSZI0sRYsmEt///JRrz99+gj9z7+H/pUDtdfvylFvG2BBzwb6V9fuY3omEGPahqSxsfCUJEnqcP39y+nrm1Ze/ysH6Ftc3tvO/tVJX9dwhaVFpzTRPNVWkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklQqC09JkiRJUqksPCVJkiRJpbLwlCRJkiSVysJTkiRJklSqUgvPiJgfEfdGxLKIOG+I5e+PiF9GxJ0R8e8RsX+Z8UiS1OrMrZKkVlRa4RkRk4BLgdcChwCnRMQhg5r9HOjOzMOBbwGfKiseSZJanblVktSqyvzE88XAssx8IDOfBq4G3lDdIDN/mJm/Lyb/E5hRYjySJLU6c6skqSVFZpbTccRJwPzMfFcx/XbgJZn5nhrtLwFWZeaFQyw7HTgdYNq0aUdfc801pcQ8XtatW8fUqVMnOowxKWsMy5bdy8DA0zWX77gjHHTQ5Nrr37eBgYHaf9PPPJ1MisrjvabP4Df9K7btHzho0naFPKHW7T2DqY+sGLlhk2uHcTiG5vDKs3puz8zuiY6jLObW2toht8LEjWPp0ruYPbt2ft0e69btxdSpv9m2/1+uZ/aMaEj/Q1n6q03MbmD+bof9oWNoDu0whkbl1sbsYYY21N5lyIogIt4GdAPHDbU8M78EfAlg1qxZOW/evAaFODF6e3txDEPr6VlIX9+00a9/1l30La79Z929cD19XZU/zd6zP8WbLv3gqLfVDHrP/hTzPtfaY4D2GIdj0Dgxt9bQDrkVJm4cY82/1Xp7FzFv3hXb9j9Cfh6rnou25vdGaIf9oWNoDu0whkYps/BcAexbNT0DeHhwo4h4FXA+cFxmPlViPJIktTpzqySpJZV5jedtwMERMTMidgJOBq6vbhARLwK+CJyYmY+WGIskSe3A3CpJakmlFZ6ZuRF4D/B9YClwTWbeHRGfiIgTi2aLganAtRGxJCKur9GdJEkdz9wqSWpVZZ5qS2beCNw4aN5Hqh6/qsztS5LUbsytkqRWVOaptpIkSZIkWXhKkiRJkspl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKpWFpyRJkiSpVBaekiRJkqRSWXhKkiRJkkpl4SlJkiRJKtXkiQ5AnWXBgrn09y+vuXzNqkfonrNq1P2veXCA7oUDNZdPzwRi1P1LktSKRsq/06ePsf/599C/spJ/F52xnp6z7tq2/64cfv2eDfSvHr7NcMzvUvOz8NS46u9fTl/ftJrLu+esom/x6P8suxcO0Nc1XOIxKUmSOs9I+XfM/a8c2JK/e5+O7c7l/atzhPw9EvO71Ow81VaSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlmjzRAai5LFgwl/7+5SO2W7Soh56ehc+av+rBR8hNG2uuN2kSdM9ZVXP5mgcH6F44UF+wQ5ieCcSo15ckqRnVm59rWbPqkWHz71hV5+9FH95Ez0Xrt2t987fU/iw8tY3+/uX09U0bsV1v7+Qh23XPWUXf4p1Hvf3uhQP0dY0l8Zi0JEntp978XEslP5f3tq86f/dOYhS53PwttTtPtZUkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWy8JQkSZIklcrCU5IkSZJUKgtPSZIkSVKpLDwlSZIkSaWaPNEBSJIktbIFC+bS37+8IX0tWtRDT8/CZ81fs+oRuuesGnW/07ty2OULejbQv3r4NsP2nwnEqNeX1P4sPCVJksagv385fX3TGtJXb+/kIfvqnrOKvsXlvW3rX530dY2lcLTolDQ8T7WVJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlsvCUJEmSJJXKwlOSJEmSVCoLT0mSJElSqSw8JUmSJEmlKrXwjIj5EXFvRCyLiPOGWL5zRHyzWP6ziDigzHgkSWp15lZJUisqrfCMiEnApcBrgUOAUyLikEHNFgGPZ+ZBwGeAvy8rHkmSWp25VZLUqsr8xPPFwLLMfCAznwauBt4wqM0bgCuLx98Cjo+IKDEmSZJamblVktSSIjPL6TjiJGB+Zr6rmH478JLMfE9Vm/8q2qwopu8v2vxmUF+nA6cXk4cB/1VK0ONnL+A3I7Zqbo6hObTDGKA9xuEYmsOszNxtooMoi7l1WO3w9wvtMQ7H0BwcQ3NohzE0JLdObkQkNQx1dHVwlVtPGzLzS8CXACKiLzO7xx7exHEMzcExNI92GIdjaA4R0TfRMZTM3FpDO4wB2mMcjqE5OIbm0C5jaEQ/ZZ5quwLYt2p6BvBwrTYRMRnoAtaUGJMkSa3M3CpJakllFp63AQdHxMyI2Ak4Gbh+UJvrgYXF45OAm7Osc38lSWp95lZJUksq7VTbzNwYEe8Bvg9MAr6cmXdHxCeAvsy8HrgC+FpELKNyNPbkOrr+UlkxjyPH0BwcQ/Noh3E4hubQDmOoydw6rHYYA7THOBxDc3AMzcExFEq7uZAkSZIkSVDuqbaSJEmSJFl4SpIkSZLK1VSFZ0TMj4h7I2JZRJw3xPKdI+KbxfKfRcQBVcs+VMy/NyJeM55xD4pxVGOIiAMiYn1ELCl+Lhvv2KtiHGkMcyPijojYWHynXPWyhRFxX/GzcPC642WMY3im6nUYfNOOcVPHGN4fEb+MiDsj4t8jYv+qZa3yOgw3hlZ5Hf4qIu4q4vxJRBxStaxV9ktDjqGV9ktV7U6KiIyI7qp5TfE6TBRza2v8DZtbx4e5tTlehyIW82sL7Juq2o09v2ZmU/xQuUnC/cCBwE7AL4BDBrV5N3BZ8fhk4JvF40OK9jsDM4t+JrXYGA4A/qtFXocDgMOBrwInVc1/HvBA8fu5xePnttIYimXrWuR1eCWwa/H4jKq/pVZ6HYYcQ4u9DrtXPT4R+F7xuJX2S7XG0DL7paLdbsAtwH8C3c30OjTzc4e5tVnGcADm1mYYg7m1ecZhfm2CMRTtGpJfm+kTzxcDyzLzgcx8GrgaeMOgNm8Ariwefws4PiKimH91Zj6Vmb8ClhX9jbexjKFZjDiGzPx1Zt4JbBq07muAH2Tmmsx8HPgBMH88gh5kLGNoFvWM4YeZ+fti8j+pfJ8ftNbrUGsMzaKeMTxZNfkcYPMd21pmvzTMGJpFPftWgAuATwEbquY1y+swUcytzcHc2hzMrc3D/NocxjW/NlPh+ULgoarpFcW8Idtk5kbgt8Ceda47HsYyBoCZEfHziPhRRBxbdrA1jOW5bKXXYThTIqIvIv4zIv6ksaHVbXvHsAj47ijXLctYxgAt9DpExJkRcT+VnfLZ27PuOBjLGKBF9ksR8SJg38z81+1dt82ZW1vkb7ikdRvJ3Nqar0Mz5lYwv0KL7JsamV9L+x7PURjqyOTgowK12tSz7ngYyxj6gf0y87GIOBr4TkQcOuhIyXgYy3PZSq/DcPbLzIcj4kDg5oi4KzPvb1Bs9ap7DBHxNqAbOG571y3ZWMYALfQ6ZOalwKUR8Rbgb4CF9a47DsYyhpbYL0XEDsBngHds77odwNzaAn/DJa7bSObWOtctWTvkVjC/tsS+qdH5tZk+8VwB7Fs1PQN4uFabiJgMdFH5cux61h0Pox5D8TH1YwCZeTuV86T/W+kRP9tYnstWeh1qysyHi98PAL3AixoZXJ3qGkNEvAo4HzgxM5/annXHwVjG0FKvQ5Wrgc1HkFvqdaiyZQwttF/aDTgM6I2IXwMvBa4vboDQLK/DRDG3tsbfcFnrNpK5tYVehybPrWB+bZV9U2PzazbBBcZZuUB1MpULtWey9eLWQwe1OZNtbx5wTfH4ULa9uPUBJuYi47GMYdrmmKlc4LsSeF4zjqGq7Vd49g0QfkXlovvnFo9bbQzPBXYuHu8F3McQF1k3wxioJIv7gYMHzW+Z12GYMbTS63Bw1eMFQF/xuJX2S7XG0HL7paJ9L1tvftAUr8NE/dT5+ptbm2AMVW2/grl1Iv+WzK3NMw7zaxOMYVD7XsaQX8f9D22EwZ8A/L/in+X8Yt4nqBytAZgCXEvl4tVbgQOr1j2/WO9e4LWtNgbgz4C7ixfwDmBBE4/hGCpHOX4HPAbcXbXuqcXYlgHvbLUxAH8I3FW8DncBi5p4DP8GPAIsKX6ub8HXYcgxtNjrcHHxv7sE+CFVO+wW2i8NOYZW2i8NattLkRib6XVo1ucOc2uzjMHc2hxjMLc2zzjMr00whkFtexlDfo1iJUmSJEmSStFM13hKkiRJktqQhackSZIkqVQWnpIkSZKkUll4SpIkSZJKZeEp/f/27j+0qjKO4/j74xTcpiVjES0EkYoIwiyL6PcfRkUFBiOwkMrojyjEmQT9AgmC6o/8IxBRwyGMIKlwGIWVW0kxY1vKjCyoZv9IP0jIrUWp3/44z5XTdu/dnXrdj/t5wcN9nrPnec73bGzfnXPuc4+ZmZmZmVWVTzzNRpF0UtIBSYck7ZTUUMV9vZwe8oyktRPdlzJ7JV2Q2mskfSup4xzE9qikllx7m6SrznCupyU9drYxmZnZ9OTcenpu51arWX6citkokoYiYl6qdwB9EfFGBeNE9jt16gz3O0j2bKTfJzDmXmB5RLSl9mGy5yj9NKrf7Ig4McF4uoH1EdE7kXEl5moAvoiIpWc7l5mZTT/OrafHdOPcajXKdzzNytsHXAYgaV26UntI0tq0bVG6CrqJ7AHACyWtlDSQ+r2W+tVJak/bBiQVklm7pFZJa4AWoEtSl6THJW0sBCHpCUnFEvTDwK7UZzOwGOiU1CZpg6QtkvYAO1Ks+yT1p3JTbv5nU1wHJb0qqRVYBnSkK9T1krolLUv9xxxj2j4k6ZU0T4+kiwEi4i9gUNIN5+SnYmZm05lzq3Or1aKIcHFxyRVgKL3OJks8TwLXAQNAIzAP+AZYCiwCTgE3pjEtwM/ARWn8XmBFGv9xbh8L0ms70Jrqg0BzqjcCPwBzUvtL4OoisR4B5ufa+Tk2AH1AfWo3AHNT/XKgN9XvSfM3pHZTeu0mu0pMvl3qGFOfAO5P9deBF3PjXwCemeyfr4uLi4vL+S/Orc6tLi6+42k2Vr2kA0AvWRJ4C7gFeD8ihiNiCHgPuDX1PxIRPal+PdAdEb9F9vabDuA24EdgsaQ3Jd0N/FkugIgYJks690m6kixJDhTp2hQRx8tM1RkRI6k+B9gqaQDYCRTWlCwHtkd25ZSI+KNcbGWOEeAfYHeq95H981DwK1liNTOz2uPcWp5zq814syc7ALMpaCQirslvSGtMShnOdy3WISKOSVoC3AU8BTwIrB4njm3A88BhYHuJPickzYrSa1/ysbUBvwBLyN5m/3cu5oks9i73vfg3IgpzneT/f2PmAiNjh5iZWQ1wbi3PudVmPN/xNKvM58AKSQ2SGoEHyNaojLYfuF1Ss6Q6YCXwmaRmYATS7G0AAAFUSURBVFZEvAu8BFxbZOxxYH6hERH7gYXAQ8DbJeL6jmztSSUuBI6mRLoKqEvb9wCrlT71T1JTsXjGO8YK9n8FcKjCWM3MbOZzbh3nGCvYv3OrTRs+8TSrQET0k60Z+YosOWyLiK+L9DsKPAd0AQeB/ojYBVwKdKe3GbWnPqNtAT6U1JXb9g7ZJ9YdKxHaB8AdFR7GJuARST1kiWo4xfwR0An0pvjWp/7twObCByBUcIzjuRn4pMJYzcxshnNudW612uLHqZhNYZJ2Axsj4tMSX78E2BERd57fyCZG0lJgXUSsmuxYzMystjm3mk0O3/E0m4IkLZD0PdmamKKJEU5fId2q9JDrKayZ7G1QZmZmk8K51Wxy+Y6nmZmZmZmZVZXveJqZmZmZmVlV+cTTzMzMzMzMqsonnmZmZmZmZlZVPvE0MzMzMzOzqvKJp5mZmZmZmVXVfwrm5RZQaImRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stdev1 = np.std(por1); stdev2 = np.std(por2)\n",
    "por1r05 = GSLIB.affine(por1,tmean=ci_95_por1[0],tstdev=stdev1)\n",
    "por1r95 = GSLIB.affine(por1,tmean=ci_95_por1[1],tstdev=stdev1)\n",
    "por2r05 = GSLIB.affine(por2,tmean=ci_95_por2[0],tstdev=stdev2)\n",
    "por2r95 = GSLIB.affine(por2,tmean=ci_95_por2[1],tstdev=stdev2)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.hist(por1r05, facecolor='yellow',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P05')\n",
    "plt.hist(por1, facecolor='orange',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='Expected')\n",
    "plt.hist(por1r95, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P95')\n",
    "\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Distribution Uncertainty - Well 1')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.hist(por2r05, facecolor='yellow',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P05')\n",
    "plt.hist(por2, facecolor='orange',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='Expected')\n",
    "plt.hist(por2r95, facecolor='red',bins=np.linspace(0.0,0.4,50),histtype=\"stepfilled\",alpha=0.8,density=True,cumulative=True,edgecolor='black',label='P95')\n",
    "\n",
    "plt.xlim([0.0,0.4]); plt.ylim([0,1.0])\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.title('Porosity Distribution Uncertainty - Well 2')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just calculated a scenario-based uncertainty model for the porosity distribution around wells 1 and 2. Of course, we could actually sample porosity means continuously from our confidence calculation as we have access to the complete distribution for uncertainty in the mean porosity of both wells.\n",
    "\n",
    "Communicating uncertainty is powerful, but always remember to state the assumptions. For example, here we assumed:\n",
    "\n",
    "1. The population distribution of porosity for each well is Gaussian distributed\n",
    "2. That the samples are independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can check the Excel file linked above with the confidence interval calculated by hand and confirm that this result is correct.\n",
    "\n",
    "#### Hypothesis Testing\n",
    "\n",
    "The confidence intervals help with uncertainty in the distributions of porosity.  Now let's try to figure out if:\n",
    "\n",
    "1. wells 1 and 2 drilled into the same type of rock? \n",
    "2. did something change between the 2 wells? \n",
    "3. different units are being compared between the 2 wells (issues with stratal correlation)?\n",
    "\n",
    "Now, let's try the t test, hypothesis test for difference in means.  This test assumes that the variances are similar along with the data being Gaussian distributed (see the course notes for more on this).  This is our test:\n",
    "\n",
    "\\begin{equation}\n",
    "H_0: \\mu_{X1} = \\mu_{X2}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "H_1: \\mu_{X1} \\ne \\mu_{X2}\n",
    "\\end{equation}\n",
    "\n",
    "For the resulting t-statistic and p-value we run this command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t statistic is -2.981 and the p-value is 0.005\n"
     ]
    }
   ],
   "source": [
    "t_pooled, p_pooled = st.ttest_ind(por1,por2)\n",
    "print('The t statistic is ' + str(round(t_pooled,3)) + ' and the p-value is ' + str(round(p_pooled,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value, $p$, is the symmetric interval probaiblity our outside.  In other words the $p$ reported is 2 x cumulative probaiblity of the t statistic applied to the sampling t distribution.  Another way to look at it, if one used the $\\pm t_{t_{statistic},.d.f}$ statistic as thresholds, $p$ is the probability being outside this symmetric interval. So we will reject the null hypothesis if $p \\lt \\alpha$.  From the p-value alone it is clear that we would reject the null hypothesis and accept the alternative hypothesis that the means are not equal.  \n",
    "\n",
    "In case you want to compare the t-statistic to t-critical, we can apply the inverse of the student's t distribution at $\\frac{\\alpha}{2}$ and $1-\\frac{\\alpha}{2}$ to get the upper and lower critcal values.       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t crical lower and upper values are [-2.02  2.02]\n"
     ]
    }
   ],
   "source": [
    "t_critical = st.t.ppf([0.025,0.975], df=len(por1)+len(por2)-2)\n",
    "print('The t crical lower and upper values are ' + str(np.round(t_critical,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe that, as expected, the t-statistic is outside the t-critcal interval.  These results are exactly what we got when we worked out the problem by hand in Excel, but so much more efficient!  \n",
    "\n",
    "Now let's try the t-test, hypothesis test for difference in means allowing for unequal variances, this is also known as the Welch's t test.  All we have to do is set the parameter 'equal_var' to false, note it defaults to true (e.g. the command above). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-2.9808897468855644, pvalue=0.005502572350112333)"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "st.ttest_ind(por1, por2, equal_var = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we can see by $p$ that we will clearly reject the null hypothesis.  \n",
    "\n",
    "Let's now compare the variances with the F-test for difference in variances.  \n",
    "\n",
    "\\begin{equation}\n",
    "H_0: \\sigma^{2}_{X2} \\le \\sigma^{2}_{X1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "H_1: \\sigma^{2}_{X2} \\gt \\sigma^{2}_{X1}\n",
    "\\end{equation}\n",
    "\n",
    "Note, by ordering the variances we eliminate the case of $\\sigma^{2}_{X2} \\lt \\sigma^{2}_{X1}$.\n",
    "\n",
    "Details about the test are available in the course notes (along with assumptions such as Gaussian distributed) and this example is also worked out by hand in the linked Excel workbook.  We can accomplish the F-test in with SciPy.Stats the function with one line of code if we calculate the ratio of the sample variances ensuring that the larger variance is in the numerator and get the degrees of freedom using the len() command, ensuring that we are consistent with the numerator degrees of freedom set as 'dfn' and the denominator degrees of freedom set as 'dfd'.  We take a p-value of $1-p$ since the test is configured to be a single, right tailed test.    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p-value for the F-test for difference is variances is 0.019\n"
     ]
    }
   ],
   "source": [
    "p_value = 1 - st.f.cdf(np.var(por2)/np.var(por1), dfn=len(por2)-1, dfd=len(por1)-1)\n",
    "print('The p-value for the F-test for difference is variances is ' + str(round(p_value,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we would clearly reject the null hypothesis since $p \\lt alpha$ and assume that the variances are not equal.  Note this is a single tail test so our p-critical is 0.05 for a 95% confidence level.\n",
    "\n",
    "#### Reporting the Hypothesis Test Results\n",
    "\n",
    "The difference in the means and variances were tested for statistical significance for porosity measured in wells 1 and 2.  In all cases we rejected the null hypotheses and adopted the alternative hypotheses that the porosity distributions' means and variances for well 1 and 2 are statistically significant in their differences. It is possible that Well 1 and 2 have drilled into different rock units or have encountered nonstationary behavoirs in a single rock unit. This result may be applied to update subsurface maps and consider a difference in the subsurface at each well for development decision making."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "We are just scratching the surface for the application of confidence intervals and hypothesis tests to subsurface modeling.  Once again, there are a lot of details left out of the problem formulation and assumptions, see the course notes for more coverage.  By running the same confidence interval and hypothesis tests: 1) by hand in Excel and with 2) R and now here in 3) Python code, I hope this demonstration will enable and encourage more engineers and scientists to make these R and Python tools part of their common practice. I'm always happy to discuss,\n",
    "\n",
    "I hope this was helpful,\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)"
   ]
  }
 ],
 "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
}