PythonDataBasics_ConfidenceInterval_HypothesisTesting.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",
    "## Data Analytics \n",
    "\n",
    "### Confidence Intervals and Hypothesis Testing in Python in Python \n",
    "\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, The University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a tutorial / demonstration of **Confidence Intervals and Hypothesis Testing in Python**.  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 and proportion\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",
    "#### Project 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",
    "\n",
    "#### Load the required libraries\n",
    "\n",
    "The following code loads the required libraries.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os                                                   # to set current working directory \n",
    "import numpy as np                                          # arrays and matrix math\n",
    "import scipy.stats as stats                                 # statistical methods\n",
    "import pandas as pd                                         # DataFrames\n",
    "import math                                                 # square root\n",
    "import statistics                                           # statistics\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",
    "#### Declare Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def welch_dof(x,y):                                          # DOF for Welch's test from https://pythonfordatascienceorg.wordpress.com/welch-t-test-python-pandas/\n",
    "        dof = (x.var()/x.size + y.var()/y.size)**2 / ((x.var()/x.size)**2 / (x.size-1) + (y.var()/y.size)**2 / (y.size-1))\n",
    "        return dof"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(\"C:\\PGE337\")                                  # 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": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     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\n"
     ]
    },
    {
     "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": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#df = pd.read_csv(\"PorositySample2Units.csv\")                # read a .csv file in as a DataFrame\n",
    "df = pd.read_csv(r\"https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/PorositySample2Units.csv\") # load data from Prof. Pyrcz's github\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": 5,
   "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": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe().transpose()   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we extract the X1 and X2 unit porosity samples from the DataFrame into separate arrays called 'X1' and 'X2' for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = df['X1'].values\n",
    "X2 = df['X2'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Interval for the Mean\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.  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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The confidence interval for the X1 mean is [0.151 0.178]\n"
     ]
    }
   ],
   "source": [
    "ci_mean_95_x1 = stats.t.interval(0.95, len(df)-1, loc=np.mean(X1), scale=stats.sem(X1))\n",
    "print('The confidence interval for the X1 mean is ' + str(np.round(ci_mean_95_x1,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Interval for the Mean By-Hand\n",
    "\n",
    "Let's not use the interval function and do each part by-hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistic +/- t-score x Standard Error\n",
      "   0.164  +/-   2.09  x     0.00606\n",
      "The confidence interval by-hand for the X1 mean is [0.15 0.18]\n"
     ]
    }
   ],
   "source": [
    "smean = np.mean(X1); sstdev = np.std(X1); n = len(X1)\n",
    "dof = len(df)-1\n",
    "tscore = stats.t.ppf([0.025,0.975], df = dof)\n",
    "SE = sstdev/math.sqrt(n)\n",
    "lower_CI,upper_CI = np.round(smean + tscore*SE,2)\n",
    "print('Statistic +/- t-score x Standard Error')\n",
    "print('   ' + str(np.round(smean,3)) + '  +/-   ' + str(np.round(tscore[1],2)) + '  x     ' + str(np.round(SE,5)))\n",
    "print('The confidence interval by-hand for the X1 mean is ' + str(np.round([lower_CI,upper_CI],3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's Plot the Result\n",
    "\n",
    "Plot the full distribution and the upper and lower confidence intervals for the uncertainty in the population mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAF6CAYAAAA5yXBHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABIHklEQVR4nO3dd5xU1f3/8dcHlt4Ell6tWAEVEbGASowair18Y0FRLIkl6verJkaNYvsmGs0vomKDaKIgNkT5WhA0KhZANEoRBaV3YYFFgd3z++PcgWHZMuXeKcv7+Xjs487ccu7nzszOZ865555rzjlEREQkd9XIdgAiIiJSOSVrERGRHKdkLSIikuOUrEVERHKckrWIiEiOU7IWERHJcUrWEiozq2dmr5nZOjN7wcx+bWZvVbL+ZDO7JJMxVheVvXbmPW1mP5rZpxHtf7CZfRBF2ZlgZn3NbFEa2z9qZn8MMyaRiihZ76LM7L/MbKqZbTCzpWY2wcyOCqHoM4BWQHPn3JnOuX86504IodyMMLORZjYswXVvN7Nno44pRUcBvwDaO+d6pluYmXU2M2dmBemHVm75t5vZluDzuNbMPjKzI6LYVyrK+2HinLvcOXdnBPu6PXitry4z/9pg/u1h71Nyn5L1LsjMrgMeBO7GJ9aOwHBgUAjFdwK+cc5tDaGsai2qxBfoBHzvnNuY7IYRx1WZ0c65hkAL4APgJTOzLMWSbd8AF5aZd0EwX3ZBSta7GDNrAtwB/MY595JzbqNzbotz7jXn3H8H69QxswfNbEnw96CZ1QmW9TWzRWZ2vZmtCGrlFwXL/gTcCpwd1JCGlK2RmNkvzGx20Ez+d8DKxHexmc0Kmm/fNLNOccucmV1uZnOD5Q/Hf5mb2aXBtuvNbKaZHRLMb2tmL5rZSjObX7bGUslrFatNXmhmC8xslZn9IVh2IvD7uGP9Ivb6mtmTweuy2MyGmVnNYNlgM/vQzP5qZmuAO4Na5IFx+2xhZpvMrKWZNTWz8UHcPwaP2ycQ9xDgCeCIILY/xb0+35rZGjMbZ2Zty7y2vzGzucDccop9P5iuDco8Im7bvwTxzTezk+LmV/haVMY5twUYBbQGmgfv37gg7m/N7NK4fdxuZmPNbHTwvk83s25ljmuvuOcVtpyY2U1m9l3c5+fUYP5+wKNxr+fa8spK4PWt8LNbjs+A+mZ2QLD9AUC9YH58zP3NbIZtb43oWtXxBMsGm9kHFb13knuUrHc9RwB1gZcrWecPQC+gO9AN6AncEre8NdAEaAcMAR42s6bOudvwtfXRzrmGzrkn4ws1s0LgxaCsQuA74Mi45afgE+Bp+NrVv4HnysTWHzgsiOss4JfBtmcCt+NrH42BgcBqM6sBvAZ8EcR7PHCtmf2ykuMv6yigS7DtrWa2n3Pu/8ocayxBjAK2AnsBBwMnAPHnlQ8H5gEt8T+aXgLOjVt+FvCec24F/v/zaXwtuSOwCfh7VcEGr/vlwJQgttvM7DjgnqD8NsAPwPNlNj0liG//coo9JpjuFpQ5Je545uDfz/8FnoxLQlW9FuUy/8NwMLDIObcK/xlYBLTFn2a528yOj9tkEPAC0Az4F/CKmdWqaj/l+A44Gv/Z/hPwrJm1cc7NYsfXc7dyYk7k9S33s1uJZ/CfZ/C17H+U2echwFPAZUBz4DFgXPD6VXg8cUVU9t5JrnHO6W8X+gN+DSyrYp3vgJPjnv8S36QK0BefNArilq8AegWPbweejVs2GPggeHwB8HHcMsN/CV8SPJ8ADIlbXgMoBjoFzx1wVNzyMcBNweM3gWvKOZbDgQVl5t0MPF3BsY8EhgWPOwf7bB+3/FPgnAqOtRXwM1Avbt65wKS416JsLP2AeXHPPwQuqCC27sCPcc8nx167ctbd9roHz58E/jfueUNgC9A57rU9rpLPROy1KCizj2/jntcP1mld1WtRTvm3A5uBtcHn6V3gUKADUAI0ilv3HmBk3Hbxn6kawFLg6Ljj2quC97cv/gdBRcc8AxhU3utZTlmJvL7lfnYreC2exf9AWwDUCqYdgvm3B+s9AtxZZts5QJ8Ej6fc966i10N/2f3L1rkpyZ7VQKGZFbiKzyu3xdcMYn4I5m0ro8y2xfgvp6q0BRbGnjjnnJktjFveCXjIzO6Pm2f4GnEsnmUV7LcD/kdGWZ2AtrGmy0BNfK09URXts7x91QKWxlVQahB3zGUeg09K9czs8GA/3QlaPcysPvBX4ESgabB+IzOr6ZwrSSJ+8K/99NgT59wGM1uNf22/ryC2RGx7bZxzxcFxN8TXcqt6Lcoa45w7L35G8Lqscc6tj5v9A9Aj7nn8Z6rUfA/v+M9rQszsAuA6/A+T2HEUJrh5Iq9vop+jWBkLzOxbfAvOXOfcwjIV307AhWZ2Vdy82kEsiRxPRe+d5CAl613PFOAnfJPn2ArWWYL/Ivg6eN4xmJeupfikCvjLi+Kf479073LO/TOFshcCe1Ywf75zbu8UyqxK2VvWLcTXJgsr+SG0wzZBchmDr3UuB8bHJabr8c3vhzvnlplZd+BzypznT1DsPQXAzBrgm04XV3I8FcadgERei0QsAZqZWaO416UjO8Yd/5mqAbRn++e1GF9rjGmNb83Zgfm+EY/jT3VMcc6VmNkMtr/WVR1/Iq9vKv6Bb+q+qJxlsf+Xu8ouSOB4JM/onPUuxjm3Dt8J7GEzO8XM6ptZLTM7ycz+N1jtOeAW852dCoP1w7hE6XXgADM7zXyP46vxX54xjwI3x3WqaRKci07EE8ANZnaoeXsFX1ifAkVmdqP5a8BrmtmBZnZYCMezHOgcJAicc0uBt4D7zayxmdUwsz3NrE8V5fwLOBt/iuJfcfMb4U85rDWzZsBtacT6L+AiM+senNO8G/jEOfd9gtuvBEqBPRJZOY3Xomw5C4GPgHvMrG7QgWoIEP+D7tC4z9S1+B8JHwfLZgD/FbzvJwIV7b8BPiGvBDDfafLAuOXLgfZmVruC7dN9fSsyGn+uf0w5yx4HLjezw4PPfAMz+5WZNUrgeCTPKFnvgpxzD+Cbx27B/zMvBH4LvBKsMgyYCnwJ/AffvJfQtcdV7HcVcCZwL745fm/8OdrY8peB+4DnzawI+ApIqIeqc+4F4C78l+b64FiaBc3FA/DNy/OBVfjE3iTd48F3agLfkS3WBHoBvilyJvAjvvWiTTnbxsf+CbAR33w5IW7Rg/gewKvwyef/Ug3UOTcR+CO+g99SfCvEOUlsX4x/fT8Meh73SmCzpF+LCpyLb8pdgj9FcJtz7u245a/if+z8CJwPnOZ8j3KAa/Dv/1r8j6FXytuBc24mcD++5Wk5cBBxn0386YqvgWVmtqqc7dN6fSvinNvknHvHObepnGVTgUvxnQ5/BL7Fn4tO5Hgkz5hzybZuiYjkBvMDhOxV9ly3SHWjmrWIiEiOizRZm9nvzOxrM/vKzJ4Lzjk1M7O3zQ8O8LaZNa26JBERkV1XZM3gZtYOP2Tg/s65TUGP1zfwAy6scc7da2Y3AU2dczdGEoSIiEg1EHUzeAH+GtIC/OUTS/CjDY0Klo/CX0IkIiIiFYgsWTvnFgN/wY+8sxRY55x7C2gVXNYRu7yjZVQxiIiIVAeRDYoSnIseBOyOv2ziBTNLuMemmQ0FhgI0aNDg0H333TeKMEVERCIzbdq0Vc65FumWE+UIZv3wI0fFLsp/CegNLA8Gx18aDCq/oryNnXMjgBEAPXr0cFOnTo0wVBERkfCZ2Q9Vr1W1KM9ZLwB6BSNkGX7Yu1nAOLbfp/VC/IAGIiIiUoHIatbOuU/MbCx+9Kut+DGNR+AHih9j/p67C/AjWomIiEgFIr2Rh/P3Ny47nvHP+Fq2iIiIJEAjmImIiOQ4JWsREZEcp2QtIiKS45SsRSJ2wQUXZDsEEclzkXYwE9nVDBw4cIfnzjkmTZrE2rVrARg3blwWohKRfKdkLRKiRYsWsf/++3PJJZdgZjjnmDp1Ktdff322QxORPKZmcJEQTZ06lUMPPZS77rqLJk2a0LdvX+rVq0efPn3o06dPtsMTkTylmrVIiGrUqMHvfvc7zjzzTH73u9/RsmVLtm7dmu2wRCTPKVmLRKB9+/a88MILvP766+y2227ZDkdE8pyawUVC9Mknn1BUVATApk2b+OSTT/jiiy+48cYbWbduXZajE5F8pWQtEqKLL76Y+vXrA3DNNddQVFTETTfdRP369bnooouyHJ2I5Cs1g4uEqLS0lIIC/281depUpk+fDsBRRx1F9+7dsxiZiOQz1axFQnTggQfy9NNPA9CtWzdi92H/5ptvqFWrVjZDE5E8pmQtEqInnniC9957jz333JOZM2dyxBFHsMcee3DppZfyxBNPZDs8EclTagYXCVGTJk0YOXIk69evZ968eWzdupX27dvTqlWrbIcmInlMNWuRCDRq1Ihu3bpx6KGHbkvUseZxEZFkKVmLZMhtt92W7RBEJE+pGVwkRF27di13vnOO5cuXZzgaEakulKxFQrR8+XLefPNNmjZtusN85xy9e/fOUlQiku+UrEVC1L9/fzZs2FDuNdV9+/bNeDwiUj2Ycy7bMVSpR48eLna9qoiISL4ws2nOuR7plqMOZiIiIjlOyVokQ5YtW5btEEQkTylZi2TIkCFDsh2CiOQpJWuRDHn99dezHYKI5Cn1BhcJ2ezZs3n11VdZvHgxZkbbtm0ZOHAg++23X7ZDE5E8pZq1SIjuu+8+zjnnHJxz9OzZk8MOOwznHOeeey733ntvtsMTkTylS7dEQrTPPvvw9ddf73Q7zM2bN3PAAQcwd+7cLEUmItmgS7dEclCNGjVYsmTJTvOXLl1KjRr6dxOR1OictUiIHnzwQY4//nj23ntvOnToAMCCBQv49ttv+fvf/57l6EQkXylZi4ToxBNP5JtvvuHTTz9l8eLFOOdo3749hx12GDVr1sx2eCKSpyJL1mbWBRgdN2sP4FbgH8H8zsD3wFnOuR+jikMk02rUqEGvXr2yHYaIVCORnURzzs1xznV3znUHDgWKgZeBm4CJzrm9gYnBcxEREalApnq8HA9855z7ARgEjArmjwJOyVAMIhn17rvv7jAVEUlVppL1OcBzweNWzrmlAMG0ZYZiEMmoG264YYepiEiqIk/WZlYbGAi8kOR2Q81sqplNXblyZTTBiWRAPoxlICK5LRM165OA6c655cHz5WbWBiCYrihvI+fcCOdcD+dcjxYtWmQgTBERkdyUiWR9LtubwAHGARcGjy8EXs1ADCIiInkr0mRtZvWBXwAvxc2+F/iFmc0NlmnAZBERkUpEOiiKc64YaF5m3mp873CRaq1hw4YANGrUKMuRiEi+02DFIhF5//33d5iKiKRKyVokRBpSVESioGQtEiJdpiUiUVCyFgmRmQFw7bXXKnGLSGiUrEUi0LBhQwYOHMjGjRsBeOuttzjyyCOzHJWI5CvdIlMkAsOGDeNf//oXffv2pU6dOjRo0IB779VViiKSGiVrkQhMnDiRxx9/nAYNGrB06VKefPJJunTpku2wRCRPqRlcJAJ33XUXd955J5MnT2bs2LGcffbZuvuWiKTM8qETTI8ePdzUqVOzHYZIlWrUqEFpaelO85cuXcrpp5/ORx99lIWoRCRbzGyac65HuuWoZi0SovISNUCbNm2YOHFihqMRkepCyVokQ+rVq5ftEEQkTylZi4iI5DglaxERkRynZC0StpISePVVeOYZCAZFERFJh5K1SJhKS+G00+h7yin0veAC6NkT1qzJdlQikueUrEXC9PjjMG4cFBRA3bowcybccEO2oxKRPKdkLRKWTZvgD3/wj/fZB7p29Ul71CiYNSu7sYlIXlOyFgnL2LGwejUceij9Bw+m/5lnwsUX+6bxRx/NdnQiksc0gplIWPr0gfffhyeegCFD/Lxp06BHDygshMWLoXbt7MYoIhmlEcxEcsnKlfDBBz4Zn3XW9vmHHAIHHACrVoHGBheRFClZi4ThjTd8c/exx0KjRvTt25e+ffuCGZxyil9nwoRsRigieUzJWiQMb7zhpwMG7LzspJP8VMlaRFKkZC2SLufg3//2j/v123n54YfDbrvB3Lnw/feZjExEqgkla5F0zZ8PS5f6TmT77LPz8oICOOoo//iDDzIbm4hUC0rWIumKJeAjj/TnqMsTS9YffpiZmESkWinIdgAieS+WrGMJGTgrvkc4+EQev66ISBKUrEXSVU6yvvLKK3dcp0cPqFULvv4aNmyAhg0zGKCI5Ds1g4uko6jIDyVau7a/pjpQXFxMcXHx9vXq1vXXWzsHX3yRhUBFJJ8pWYuk4z//8dMDDthhdLKTTz6Zk08+ecd1u3f3088/z0xsIlJtKFmLpGPGDD/t1q3qdQ8+eMdtREQSFGmyNrPdzGysmc02s1lmdoSZNTOzt81sbjBtGmUMIpGKJd5YrbkysWStmrWIJCnqmvVDwP855/YFugGzgJuAic65vYGJwXOR/BQ7/5xIso7Vvr/6CrZsiSwkEal+IkvWZtYYOAZ4EsA5t9k5txYYBIwKVhsFnBJVDCKR2rp1+znrrl2rXr9xY9hrL9i8GWbOjDY2EalWorx0aw9gJfC0mXUDpgHXAK2cc0sBnHNLzaxlhDGIRGfuXPjpJ+jUCZrueDZn8ODB5W/TrRt8+62vXSdynltEhGiTdQFwCHCVc+4TM3uIJJq8zWwoMBSgY8eO0UQoko6vvvLTgw7aaVGFyXq//fx01qxoYhKRainKc9aLgEXOuU+C52PxyXu5mbUBCKYrytvYOTfCOdfDOdejRYsWEYYpkqI5c/x03313WrRq1SpWrVq18zZK1iKSgsiStXNuGbDQzLoEs44HZgLjgAuDeRcCr0YVg0ikYsm6S5edFp1xxhmcccYZO28TS9Y6Zy0iSYh6uNGrgH+aWW1gHnAR/gfCGDMbAiwAzow4BpFozJ7tp+XUrCvUpYu/2ce33/oe4bVqRRObiFQrkSZr59wMoEc5i46Pcr8ikXOu0pp1herX9x3Svv/eJ+xYTVtEpBIawUwkFcuWwfr1vhd4YWFy2+q8tYgkSclaJBXxteqK7mFdESVrEUmSbpEpkooqmsCvuOKKireNneOOlSEiUgUla5FUVJGszz777Iq33WsvP/3uu5CDEpHqSs3gIqmo5BprgIULF7Jw4cLyt91zTz9VshaRBKlmLZKK2GVbFdSszz//fAAmT56888J27fy9r5cvhw0boGHDiIIUkepCNWuRZG3e7C+9qlFjey05GTVrwu67+8fz5oUamohUT0rWIslasABKS6F9e6hTJ7UyYuetv/02vLhEpNpSshZJVqw2vMceqZeh89YikgQla5FkKVmLSIapg5lIsubP99NKkvX1119feRlqBheRJChZiyQrVrOOdRIrx4ABAyovQzVrEUmCmsFFkpVAM/icOXOYU9kIZZ07+2FKFyzwvctFRCqhZC2SrASawS+77DIuu+yyisuoUwc6dvS9yn/4IeQARaS6UbIWScaPP/q/+vWhRYv0yoo1heu8tYhUQclaJBnxtepk77ZVVufOfqqatYhUQclaJBkJNIEnrFMnP1WyFpEqKFmLJCOBnuAJU7IWkQTp0i2RZCRYs77llluqLkvJWkQSpGQtkowEa9b9+vWruiwlaxFJkJrBRZKxaJGfduxY6WozZsxgxowZlZfVvr2/c9eSJbrWWkQqpZq1SDJiybp9+0pXu/baa4EK7mcdU6sWtG3ry1y0KJxOayJSLalmLZKo9eth3TqoVw+aNQunTDWFi0gClKxFEhVfq073GuuYWLJesCCc8kSkWlKyFklUgk3gSYmd+1bNWkQqoWQtkqgokrWawUUkAepgJpKoJJL13XffnViZStYikgAla5FEJZGse/funViZStYikgA1g4skKolk/dFHH/HRRx9VXWZ8B7PS0jSCE5HqTDVrkUQlkax///vfA1VcZw3QoAE0bw6rV8Py5dCmTZpBikh1pJq1SKKi6GAGagoXkSpFmqzN7Hsz+4+ZzTCzqcG8Zmb2tpnNDaZNo4xBJBTFxbBmDdSuDYWF4ZatZC0iVchEzfpY51x351yP4PlNwETn3N7AxOC5SG5bvNhPY+N5h0nJWkSqkI1m8EHAqODxKOCULMQgkpyFC/007CZwgA4d/DT2g0BEpIyoO5g54C0zc8BjzrkRQCvn3FIA59xSM2sZcQwi6UvyfPWDDz6YeNmxMmM/CEREyog6WR/pnFsSJOS3zWx2ohua2VBgKEDHKm5HKBK5JJN19+7dEy87VmZsHyIiZUTaDO6cWxJMVwAvAz2B5WbWBiCYrqhg2xHOuR7OuR4tWrSIMkyRqiWZrN955x3eeeedxMpWshaRKkSWrM2sgZk1ij0GTgC+AsYBFwarXQi8GlUMIqFJMlkPGzaMYcOGJVZ2mzb+Ll7LlsGWLSkGKCLVWZTN4K2Al83fSrAA+Jdz7v/M7DNgjJkNARYAZ0YYg0g4orrGGqBWLZ+wlyyBpUu334lLRCQQWbJ2zs0DupUzfzVwfFT7FYlElMk6Vu6SJX4/StYiUoZGMBOpyk8/wcqVUFAALSO6eEHnrUWkEkrWIlVZssRP27aFmjWj2YeStYhUQjfyEKlK7Prn2OAlCXjssceS24eStYhUQslapCopnK/u0qVLcvtQshaRSqgZXKQqKSTr1157jddeey3xfShZi0glVLMWqUoKyfr+++8HYMCAAYltoGQtIpVQzVqkKlFftgW+8xr4zmwlJdHtR0TykpK1SFUykazr1PGXhZWU+JHMRETiKFmLVCUTyTq+fDWFi0gZStYildm8GZYvhxo1oHXraPelZC0iFVAHM5HKLF0KzvlzygWJ/7s888wzye8rdh23krWIlKFkLVKZ2IAoSTaBd0hiAJVtVLMWkQqoGVykMrHEmWTyHT16NKNHj05uX0rWIlIB1axFKpNi57JHHnkEgLPPPjvxjZSsRaQCqlmLVCZTPcHj96FkLSJlKFmLVCaTybpdOz9dvBhKS6Pfn4jkDSVrkcpkMlnXqwfNm8OWLf7+2SIiASVrkcpkMlnH70dN4SISJ6FkbWYvmtmvzEzJXXYdW7f666zNoE2bpDYdO3YsY8eOTX6fStYiUo5Ek+8jwH8Bc83sXjPbN8KYRHLDsmX+3HGrVlC7dlKbFhYWUlhYmPw+Y8k6dn23iAgJJmvn3DvOuV8DhwDfA2+b2UdmdpGZ1YoyQJGsSaMJfOTIkYwcOTL5fapmLSLlSLhZ28yaA4OBS4DPgYfwyfvtSCITybYURy+DEJL14sXJbysi1VZCg6KY2UvAvsAzwADn3NJg0WgzmxpVcCJZleLoZWlRzVpEypHoCGZPOOfeiJ9hZnWccz8753pEEJdI9mW6Jzhsv9ZayVpE4iTaDD6snHlTwgxEJOdkI1nHN4M7l7n9ikhOq7RmbWatgXZAPTM7GLBgUWOgfsSxiWRXNpJ1o0bQuDEUFcGPP0KzZpnbt4jkrKqawX+J71TWHnggbv564PcRxSSSG9JI1m+88UbVK1WkXTufrBctUrIWEaCKZO2cGwWMMrPTnXMvZigmkewrKYElS/zjtm2T3rx+/TQantq3h1mzfFN4166plyMi1UZVzeDnOeeeBTqb2XVllzvnHihnM5H8t2KFH8GsRQuoWzfpzYcPHw7AlVdemfy+1SNcRMqoqhm8QTBtGHUgIjklzfPVY8aMAVJM1vF33xIRoepm8MeC6Z8yE45IjkhjQJS0qWYtImUkeiOP/zWzxmZWy8wmmtkqMzsvwW1rmtnnZjY+eN7MzN42s7nBtGk6ByASiWz0BI/RtdYiUkai11mf4JwrAvoDi4B9gP9OcNtrgFlxz28CJjrn9gYmBs9Fcks2Ri+L0ZCjIlJGosk6drOOk4HnnHNrEtnIzNoDvwKeiJs9CBgVPB4FnJJgDCKZk82atZrBRaSMRIcbfc3MZgObgCvNrAXwUwLbPQj8D9Aobl6r2NjizrmlZtayvA3NbCgwFKBjx44JhikSkjST9eTJk1Pfd/PmUKcOrF0LGzdCgwZVbiIi1Vuit8i8CTgC6OGc2wJsxNeQK2Rm/YEVzrlpqQTmnBvhnOvhnOvRokWLVIoQSV02a9Zm6hEuIjtItGYNsB/+euv4bf5RyfpHAgPN7GSgLtDYzJ4FlptZm6BW3QZYkXTUIlEqLd2eJGNJM0l/+ctfALjhhhtSi6F9e5g3z/9o2Gef1MoQkWoj0d7gzwB/AY4CDgv+Kr3blnPuZudce+dcZ+Ac4F3n3HnAOODCYLULgVdTC10kIqtWwebNfqjPFEciGz9+POPHj089BvUIF5E4idasewD7OxfKbYDuBcaY2RBgAXBmCGWKhCebTeAx6hEuInESTdZfAa2BpansxDk3GZgcPF4NHJ9KOSIZkUvJWjVrESHxZF0IzDSzT4GfYzOdcwMjiUokm7I5elmMmsFFJE6iyfr2KIMQySkh1Kzr1auXXgxqBheROAkla+fce2bWCdjbOfeOmdUHakYbmkiWhDB62YQJE9KLQc3gIhIn0d7glwJjgceCWe2AVyKKSSS7cuGcdatWUKMGLF/ue6aLyC4t0eFGf4O/broIwDk3Fyh35DGRvBdCsr7zzju58847U4+hoADatPGPl6bUr1NEqpFEk/XPzrltP++DgVHCuIxLJLc4tz1ZpzggCsDEiROZOHFierGoKVxEAokm6/fM7PdAPTP7BfAC8Fp0YYlkyZo18NNP0KQJNGpU9fpRUo9wEQkkmqxvAlYC/wEuA94AbokqKJGsyYXz1THqES4igUR7g5ea2SvAK865ldGGJJJFuZisVbMW2eVVWrM273YzWwXMBuaY2UozuzUz4YlkWEgDojRv3pzmzZunF4uawUUkUFXN+lp8L/DDnHPzAcxsD+ARM/udc+6vEccnklkh1axffPHF9GNRM7iIBKo6Z30BcG4sUQM45+YB5wXLRKoXNYOLSA6qKlnXcs6tKjszOG9dK5qQRLIohNHLAG6++WZuvvnm9GJp29ZPlyzx99gOwZo1a7jjjjt44okncM5x11130b9/f/77v/+bH3/8MZR9iEj4qkrWlQ2dpGGVJKO+/PLLbY+3bNnCsGHDGDhwIL///e8pLi4OZych1aynTJnClClT0oulbl0oLIStW2HFivTKCpx33nls3LiRadOmceyxx7Js2TJuvPFG6tWrx+DBg0PZh4iEr6pz1t3MrKic+QbUjSAekQoNHjyY6dOnA3DTTTexevVqrr/+el555RUuv/xy/vGPf6S3g/gBUXKhGRx8HKtW+bhat067uCVLlvDGG2/gnKN9+/ZMnjwZgKOPPpru3bunXb6IRKPSZO2c0806JGc4t33QvIkTJ/LZZ59Rq1YtjjnmGLp165b+Dtatg40boWFDaNw4/fLC0K4dzJjhk3WPHmkXV1payo8//sj69evZsGED33//PZ07d2b16tVs1hjkIjkr0VtkimTdunXrePnllyktLeXnn3+mVi3fbcLMMLP0dxBfqw6jvDCE3CP85ptvZt999wXgqaee4pJLLsHMmDlzJrfddlso+xCR8ClZS97o06cP48aNA6BXr14sX76cVq1asWzZMgoLC9PfQYhN4O3DakYP+Vrrc889l7POOgvnHAUFBQwaNIgZM2bQrl072sRuHCIiOUfJWvLG008/Xe781q1bp3/TDAhtQBSAZ599Nu0ygEgu36pZc/vZrYKCAnoEzeuzZ8/eVusWkdyS6NjgIjnt7bffTr+QWLJO87KtUGVwYJQTTjgh8n2ISGpUs5ZqYciQISxYsCC9QkJM1tdeey0ADz74YHoFhdwMfvXVV5c73znH2rVrQ9mHiIRPyVryxsCBA8ud75xj9erV6e8gxGQ9Y8aMtMsAdmwGdy7tjm9PP/00999/P3Xq1Nlp2XPPPZdW2SISHSVryRv//ve/efbZZ2nYsOEO851zfPrpp+nvIBebwRs39vfVXr8e1q6Fpk3TKu6www7jwAMPpHfv3jstu/3229MqW0Sio2QteaNXr17Ur1+fPn367LSsS5cu6RXuXG4ma/BN4bNn+9p1msl67Nix1K1b/nhG8+fPL3e+iGSfkrXkjQkTJlS47P3330+v8DVrYNMmX5PNlQFRYtq3356sDzooraKaNWsWUlAikknqDS4Codeq99lnH/bZZ59QyspUj3A1g4vkLiVrqRaGDh2aXgEhXmMNMGLECEaMGBFKWWH3CK/IoYceGmn5IpI6JWupFi677LL0CsjV89WQsftaDxgwINLyRSR1Omct1ULatcKQk3Wsph9K7TrkZvA333yTV155hcWLF2NmtG3blkGDBnHiiSeGUr6IhE/JWvLG1q1befLJJ3n55ZdZsmTJDolmyJAh227skZKQk/U333wTSjlAqM3g1157Ld988w0XXHDBtvHLFy1axN/+9jcmTJjAQw89lPY+RCR8kSVrM6sLvA/UCfYz1jl3m5k1A0YDnYHvgbOccz9GFYdUH+effz677bYbt99++w6JZtSoUZx33nmMHj069cJjiTCXm8FDqFm/8cYb5f6QOPvss9lnn32UrEVyVJQ165+B45xzG8ysFvCBmU0ATgMmOufuNbObgJuAGyOMQ6qJ6dOnM2fOnB3mtW/fnl69eqXf8zqXz1kXFkLt2vDjj/5+2w0apFxU3bp1+fTTT+nZs+cO8z/77LMKr78WkeyLLFk75xywIXhaK/hzwCCgbzB/FDAZJWtJQNOmTXnhhRc4/fTTqVHD940sLS3lhRdeoGk6g4WUluZ2zdrM167nzfM/KtK4M9bIkSO54oorWL9+/bbWiYULF9K4cWNGjhwZUsAiErZIz1mbWU1gGrAX8LBz7hMza+WcWwrgnFtqZi2jjEGqj+eff54bb7yRK6+8cltyXrt2LcceeyzPP/986gWvXAmbN0OzZlC/fiixdu/ePZRytunYMZRkfcghh/DJJ5+wbNkyFi9ejHOO9u3b07p16xCDFZGwRZqsnXMlQHcz2w142cwOTHRbMxsKDAXo2LFjNAFKXuncufO289KrV6/GOUdhYWH6BUfQBJ723bbK6tTJT3/4IZTiWrdurQQtkkcycp21c24tvrn7RGC5mbUBCKYrKthmhHOuh3OuR4sWLTIRpuSR5s2bh5OoIbfPV8fEfrCmextQEclLkSVrM2sR1Kgxs3pAP2A2MA64MFjtQuDVqGIQSUgEyfq8887jvPPOC628sGvWIpJfomwGbwOMCs5b1wDGOOfGm9kUYIyZDQEWAGdGGINI1SJI1ovCHm1MNWuRXVpkNWvn3JfOuYOdc12dcwc65+4I5q92zh3vnNs7mK6JKgapvoqKinaYpiUfmsEjqFkfddRRO0xFJHdpbHDJS3379t1hmpZ8SNax2BYtgpKSUIosLi4GYOPGjaGUJyLRUbKWvOYv509TPiTrevWgZUvYsgWWLct2NCKSYRobXHZtW7fCkiX+cWwM7hAcccQRoZW1TceOsGKFP28dYqwikvuUrGXXtnixb1Zu2xbq1Amt2HvuuSe0srbp1AmmTvXnraP4MSAiOUvN4JLXzCy9AmIdtmIduHJZyD3CQzmFICIZoWQteSmWaNJOON9/76edO6dXThmnn346p59+eqhlbkvWIfUI/+tf/7rDVERyl5rBJS/Fhh1N67aYEFnNevXq1aGWB2yPMaSadag96kUkUqpZS16K3RIz7VtjRlSzjkTINWsRyR9K1rJry6dz1iHXrEUkfyhZy64tn2rWzZv7663XrfN/IrLL0Dlr2XWVlm6vpYZ8G9bjjz8+1PIAMPO169mzfdwHHZR0EVdddVWlPej/9re/pROhiEREyVryxoABAypNNOPGjUuuwKVL/YhgLVtC/fppRrejP/7xj6GWt03Hjmkl6x49ekQQlIhETcla8sYNN9wQboH5dL46Js0belx44YVVryQiOUfJWvJGnz59wi0wwvPVJ510EgATJkwIt+CQBkZZuXIl9913HzNnzuSnn37aNv/dd99Nq1wRiYY6mEnemTt3LmeccQb7778/e+yxx7a/pEVYs960aRObNm0KvdywbpX561//mv3224/58+dz22230blzZw477LAQAhSRKChZS9656KKLuOKKKygoKGDSpElccMEFnH/++ckXlE89wWNiscZiT9Hq1asZMmQItWrVok+fPjz11FN8/PHHaYcnItFQspa8s2nTJo4//nicc3Tq1Inbb789tebbfDxnvfvufjpvXlrF1KpVC4A2bdrw+uuv8/nnn7No0aJ0oxORiOicteSdunXrUlpayt57783f//532rVrx4oVK5IvKJas86lm3bYt1K7tb5W5YQM0bJhSMbfccgvr1q3j/vvv56qrrqKoqEhjhIvkMCVryTsPPvggxcXF/O1vf+OPf/wjkyZNYtSoUckV4lykNev+/fuHXiYANWr42vWcOTB/fkqXbwEcfvjhNGnShCZNmjBp0qSQgxSRsClZS94pKCigYcOGNGzYkKeffjq1QlauhE2boFkzaNQo3ACJ4DKzeHvskXay7t27N7vvvjtnn302p512Gk2bNg05SBEJk85ZS9657rrr2HffffnjH//I119/nVohsXO+sXPA+STW8z2N89Zz585l2LBhfP311xx66KH079+fZ599NqQARSRsStaSdyZNmsTkyZNp0aIFQ4cO5aCDDmLYsGHJFfLdd366557hB4i/7WRkt54MqZNZz549eeCBB/j0009p1qyZBkwRyWFK1pKXWrduzdVXX82jjz5K9+7dueOOO5IrIOJkHakQatZFRUWMGjWKk046id69e9OmTRs+/fTTkAIUkbDpnLXknVmzZjF69GjGjh1L8+bNOeecc7j//vuTKySWrFMZTCXbQkjW3bp145RTTuHWW2/liCOOCCkwEYmKkrXknYsuuohzzz2Xt956i7Zt26ZWSD7XrGPN4PPn+17tldzcpDwlJSWceuqpPPDAAxEEJyJRULKWvFJSUsKee+7JNddck15B+ZysGzeGwkJYtQqWLYM2bZLavGbNmnzxxRcRBSciUVCylrxSs2ZNVq9ezebNm6ldu3ZqhWzc6JNc7drQrl24AQbOOuusSMrdZvfdfbKeNy/pZA3QvXt3Bg4cyJlnnkmDBg22zT/ttNPCjFJEQqJkLXmnU6dOHHnkkQwcOHCHRHPdddclVkD8ZVs1a0YQIVx55ZWRlLvNHnvAZ5/5YznyyKQ3X7NmDc2bN99hmFYzU7IWyVFK1pJ32rZtS9u2bSktLWX9+vXJF5CBJvDi4mIA6tevH80O0uxklvJgMiKSFUrWknduu+02ADZu3LhDzTphGUjWJ598MgCTJ0+OZgexZD1/fkqbf/PNN1xxxRUsX76cr776ii+//JJx48Zxyy23hBikiIRF11lL3pkyZQr7778/++23HwBffPFFcs3O+dy5LCbNmvWll17KPffcs+3uW127duX5558PKzoRCVlkydrMOpjZJDObZWZfm9k1wfxmZva2mc0NphqUWJJy7bXX8uabb9K8eXPAXzP8/vvvJ15ALMFVh2T97bcpbV5cXEzPnj13mFdQoIY2kVwVZc16K3C9c24/oBfwGzPbH7gJmOic2xuYGDwXSUqHDh12eF4zmY5i+TwgSkyHDr43+9Kl/laZSSosLOS7777Dgmu0x44dS5sUepWLSGZE9lPaObcUWBo8Xm9ms4B2wCCgb7DaKGAycGNUcUj106FDBz766CPMjM2bN/O3v/1tW5N4lbZuhe+/94/z8SYeMTVrwl57wcyZMHcuHHxwUps//PDDDB06lNmzZ9OuXTt23313/vnPf0YUrIikKyPtXmbWGTgY+ARoFSRynHNLzaxlJmKQ6uPRRx/lmmuuYfHixbRv354TTjiB4cOHJ7bxwoU+YbdrB/XqRRbj4MGDIyt7m3328cn6m2+STtZmxjvvvMPGjRspLS2lUaNGzE+xs5qIRC/yZG1mDYEXgWudc0WW4NCIZjYUGArQsWPH6AKUvDNnzpydaoEffvghRyZyvXGGOpdlLFmDT9ZJOv3005k+ffoOvenPOOMMpk2bFlZ0IhKiSJO1mdXCJ+p/OudeCmYvN7M2Qa26DbCivG2dcyOAEQA9evRwUcYp+eWqq65i+vTpVc4r15w5frr33hFEtt2qVasAf244Ml26+GkSyXr27Nl8/fXXrFu3jpdeemnb/KKiIn766aewIxSRkESWrM1XoZ8EZjnn4u8YMA64ELg3mL4aVQxSvUyZMoWPPvqIlStX7nATiqKiIkpKShIrJJasY4kuImeccQYQ4XXWkFLNes6cOYwfP561a9fy2muvbZvfqFEjHn/88bAjFJGQRFmzPhI4H/iPmc0I5v0en6THmNkQYAFwZoQxSDWyefNmNmzYwNatW3cYuaxx48aMHTs2sUIylKwzIj5ZJ3j3rUGDBjFo0CCmTJmiW2OK5JEoe4N/AFT07XF8VPuV6qtPnz706dOHwYMH06lTp9QKqU7JukULaNIE1q6FlSuhZeJ9Nffaay/uvvtuvv/+e7Zu3bpt/lNPPRVBoCKSLo2CIHnn559/ZujQoTslmvibUpRr0yZYsAAKCvL7GusYM1+7/uwzX7tOIlkPGjSIo48+mn79+iV3jbqIZIWSteSdM888k8svv5xLLrkkuUQzd65vLt5jDwiG2cx78cn6qKMS3qy4uJj77rsvwsBEJExK1pJ3CgoKuOKKK5LfMINN4CnFl4oUL9/q378/b7zxxrYbjohIblOylrwzYMAAhg8fzqmnnkqdOnW2zW/WrFnlG2YwWZ999tmR7wNIOVk/9NBD3H333dSpU4datWrhnMPMKCoqiiBIEUmXkrXknVGjRgHw5z//eds8M2NeVXegymCyXrhwIbDzGOahSzFZp3QfcBHJGiVryTspD4uZwWR9/vnnAxFfZw3bB3f59lsoKfFjhldi9uzZ7LvvvhUOIHPIIYeEHaGIhEDJWvLGu+++y3HHHbfDyFvxTjvttIo3dq56XbYV06iRH+d88WKYP9/f3KMSDzzwACNGjOD666/faZmZVd2jXkSyQsla8sZ7773Hcccdt8PIWzFmVnmyXrIEioqgaVN/fXJ1csABPll//XWVyXrEiBEATJo0KRORiUhIlKwlb/zpT38C4Omnn05+46+/9tMDD0xopK+8csAB8NZb/hgHDcp2NCISgRrZDkAkI776yk8PPDC7cUThgAP8NPaDRESqHdWsZdeQ4WRd3jnhyMSSdewYq+CcY9GiRdH3VBeR0KhmLXmltLSUjz76KPkNY4ksltgiNmDAAAYMGJCRfbH//n46ezbEDb9aETPjlFNOiTYmEQmVkrXklRo1aiRfay0thZkz/eMMJes5c+YwJ9b7PGqNG0PHjrB5M3z3XUKb9OrVi88++yziwEQkLErWkndOOOEEXnzxRZxziW3www+wcSO0bg2FhdEGF7jsssu47LLLMrIvIOnz1pMmTaJXr17sueeedO3alYMOOoiuXbtGGKCIpEPnrCXvPPDAA2zcuJGaNWtSr169qofKrM6dy2IOOAAmTPDJurJL2AITJkzIQFAiEhbVrCXvrF+/ntLSUrZs2UJRURHr16+vfEzrWG0zQ03gWZFkzbpTp04sXLiQd999l06dOlG/fn1KS0sjDFBE0qFkLXnHOcezzz7LnXfeCfhxuD/99NOKN9gVatYHHeSnX36Z0Op/+tOfuO+++7jnnnsA2LJlC+edd15U0YlImpSsJe9ceeWVTJkyhX/9618ANGzYkN/85jcVbxBLYNU5WR9wgB8XfM4cKC6ucvWXX36ZcePG0aBBAwDatm2rm3uI5DCds5a888knnzB9+nQOPvhgAJo2bcrmzZvLX3nTJt8TvEYNyGAHqltuuSVj+wKgbl1/Cdd//uP/Dj+80tVr166NmWHBaG4bN27MRJQikiLVrCXv1KpVi5KSkm2JZuXKldSoUcFH+auv/N2o9t0X6tfPWIz9+vWjX79+GdsfAN27++nnn1e56llnncVll13G2rVrefzxx+nXrx+XXnpptPGJSMpUs5a8c/XVV3PqqaeyYsUK/vCHPzB27FiGDRtW/sqxW0Fm+NaPM2bMAKB7LIFmwsEHwzPPJJSsb7jhBt5++20aN27MN998wx133MEvfvGLDAQpIqlQspa88+tf/5pDDz2UiRMn4pzjlVdeYb/99it/5VjiCprMM+Xaa68FMnA/63ixY0wgWQMcdNBBbNq0CTPjoFgHNRHJSWoGl7xz6623snDhQgYPHsxvf/vbihM1ZK1mnRWxWvx//lPlsKNPPPEEPXv25KWXXmLs2LH06tWLp556KvoYRSQlqllL3uncuTPPPfccV199NY0aNeLoo4/mmGOOYVDZ20Nu3eoTF2xPZNXZbrvB7rvD/Pm+V3gl15X/+c9/5vPPP6d58+YArF69mt69e3PxxRdnKFgRSYZq1pJ3Lr74Yp566ikmTZrEeeedxwsvvFD+NcKzZ8NPP8Eee/hEtitIsCm8ffv2NGrUaNvzRo0a6S5cIjlMNWvJO5dccgkzZ86kVatWHH300YwdO5ZDymvmjjWBZ/h8dVYdfDC89BJMnQqVDHLSrl07Dj/8cAYNGoSZ8eqrr9KzZ08eeOABAK677rpMRSwiCVCylryzevVqSkpK2G233WjWrBmFhYUUFJTzUY7dVSoL56vvvvvujO8TgF69/PSTTypdbc8992TPPffc9jx2CkEDo4jkJkv4zkVZ1KNHDzd16tRshyE5ZtasWbz55pv89a9/paSkhEWLFu24Qo8eMG0avPsuHHtsdoLMtHXroGlTqFULioqgTp1KVy8qKsLMdmgSF5HwmNk051yPdMtRzVryzvjx4/n3v//N+++/z48//shxxx3H0UcfveNKxcXwxRd+5LLDDst4jB999BEAvXv3zuyOmzSB/fbzo7bNmFHhSGZTp07loosu2laTbtKkCU899RSHHnpoBoMVkUQpWUvemTBhAscccwzXXHMNbdu2LX+ladN8b/Bu3aBhw8wGCPz+978HMnyddUyvXj5Zf/xxhcn64osvZvjw4dt+5HzwwQdcdNFFfJngjUBEJLMi6w1uZk+Z2Qoz+ypuXjMze9vM5gbTplHtX6qvhx9+mL59+zJ9+nTGjx/PihUrdl5pyhQ/PeKIzAaXC2LnrT/+uMJVYpe8xRx11FFqChfJYVFeujUSOLHMvJuAic65vYGJwXORpLzwwgv07NmTF154gTFjxnD44YczduzYHVdSsq40Wffs2ZPLLruMyZMn895773HllVdu+wE0PdaLXkRyRqQdzMysMzDeOXdg8HwO0Nc5t9TM2gCTnXNdqipHHcwkXrdu3Xj77bdp2bIl4G/k0a9fP7744gu/gnPQpg0sXw7ffAN7753xGPv27QtkqRm8pMRfV75hAyxdCq1b77TKsZV0uDMz3n333QgDFNl15GsHs1bOuaUAQcJumeH9SzVQWlq6LVEDNG/enNLS0u0rfP+9T9TNm8Nee2U+wGyrWdOfq544ET74AM44Y6dVJk2alIXARCRVOdvBzMyGAkMBOnbsmOVoJJeceOKJ/PKXv+Tcc88FYPTo0Zx88snbV3j/fT/t3RuC22hm2oMPPpiV/W7Tt69P1pMmlZus77jjjnI3u/XWWyMOTERSkelkvdzM2sQ1g5fTM8hzzo0ARoBvBs9UgJL7/vznP/Piiy/y4Ycf4pxj6NChnHrqqdtXiDXhHndcdgIkw7fGLE/QDE8FzfANGjTY9vinn35i/Pjxld8QRUSyKtPnrP8MrHbO3WtmNwHNnHP/U1U5OmctCXMOOnaERYv8ddZdu2YljHfeeQeAfv36ZWX/bN7sB0cpLoZly6BVq0pX//nnnxk4cCBvvvlmhgIU2TXk/DlrM3sO6AsUmtki4DbgXmCMmQ0BFgBnRrV/qX4aNWqEldOs7ZzDzCgqKoJvv/WJurAQDjwwC1F6w4YNA7KYrGvXhiOPhLff9rXrs8+udPXi4mLmzZuXmdhEJGmRJWvn3LkVLDo+qn1K9ZbQuNWxJvC+ff3oZbuyvn0rTNYHHXTQth8+JSUlrFy5UuerRXJYznYwE0lJDpyvzhmxy7PK6fk9fvz4bY8LCgpo1apV+TdDEZGcoP9OqT5KS7cnJiVrfyOTxo1hzhyYPx92333bok6dOmUxMBFJ1i7eTijVyqefwsqV0KkT7LNPtqPJvlq14IQT/OPXX89uLCKSFiVrqT5ee81PBwzI2vXVMY899hiPPfZYVmMAoH9/P1WyFslragaX6iN2HnbAgOzGAXTpUuUouplx0kn+h8ukSbBxI8RdXy0i+UM1a6kefvgBvvzS3w6zT59sR8Nrr73Ga7Gafja1bAk9e8LPP/sRzUQkLylZS/UQq1X/8pdQp052YwHuv/9+7r///myH4f3qV346blx24xCRlClZS97q27fvtrtb8fLLfpoDTeA555RT/PTll/3IZnF2eA1FJGcpWUv+W7LEX19duzYMGpTtaHLPgQfC/vvDmjV+kBQRyTtK1pL/nn/ejwn+q1/5+zjLjswguEMZzz2X3VhEJCVK1pL//vlPP/31r7MbRy475xw/feUVf3MPEckrkd51Kyy665aUZ+HChfDtt3Q47jg/Utfy5VC3brbDAoLYgA4dOmQ5kjg9e8Jnn/nadZC8czJOkWokrLtuqWYteatDhw50mDDBPzn99JxJ1BDElmsJ8IIL/DRusJacjFNEdqJkLXlr9DPPMHr4cP/k8suzG0wZo0ePZvTo0dkOY0fnn+8HRZk8GWbNAnI0ThHZiZK15K1H7rqLRzZuhMMO8028OeSRRx7hkUceyXYYO2rSBM47zz8OYsvJOEVkJ0rWkp+cg8WL/ePf/ja7seSTK67w01GjIJH7g4tITlCylvw0eTJs2AA1a8JZZ2U7mvzRrRscfTQUFcHDD2c7GhFJkJK15Kc77vDTgoKc6liWF265xU/vvx9KSrIbi4gkRMla8s977/maNfh7NktyfvEL6NULVq3yo7+JSM7TLTIlv5SWwvXXAzC2eXPo3j278VRg7Nix2Q6hYmZw661w8smMXbMGYpe/iUjOUs1a8suzz8K0aVC7NoV7701h7drZjqhchYWFFBYWZjuMip14Ihx3HIVr11L4179mOxoRqYKSteSPFSu21app25aRRUWMXLQouzFVYOTIkYwcOTLbYVTMDP7f/2NkjRqMfPRR0AiBIjlNyVryg3Nw5ZX+PGujRtClCyOXLGFk7PKtHJPzyRpg//0Z2bYtIwFOPRU2bcpyQCJSESVryQ/PPAMvvuh7f++1lzqWhaVzZ1/LXrQIfvObbEcjIhVQspbc99lnMHSof9y2rf+TcNSs6W+CYgZPP+0HSxGRnKNkLblt/nw45RT4+Wdo3hz2398nFglPQQHsvbd/PGQIvPNOduMRkZ0oWUvumjcP+vb11wI3bOhH3yrQ1YaR6NDBN4mXlMCAATBpUrYjEpE4up+15KZPPvGdnpYu9Ym6e3fYbbcdVineuBFq1aL+m29mJcTKFBcXA1C/fv0sR1K54uJiOPts6m/a5EeCmz7dv+YFBfDPf2ooV5E06X7WUj2Vlvoxq485xieNRo3KTdQA9WvWpH7NmhkPMRH169fP+UQNQZyx1gozOOQQX8veuhXOPhuuugp++im7QYqIkrXkkGnT4Nhj/V20Nm+GwkJ/+8tyEjXA8IULGf7DD5mNMUHDhw9neOxe2zls+PDhDJ83b/sMM+jaFfbd1z/++999P4EJE/zlcyKSFUrWkl3OwYcfwmmnQY8e8P77ULs27LGHT9SV1E7HLF/OmGXLMhhs4saMGcOYMWOyHUaVxowZw5iy16qb+cvjevf2TePz58PJJ8MRR8Arr/hat4hkVFZ665jZicBDQE3gCefcvdmIQ7KkpMRfjjVhAowZA7Nn+/k1akCLFj5RNG2a3RjFvwfHHQdz5/rOfrF+BC1bwrnnQv/+cOSRUK9etiMVqfYynqzNrCbwMPALYBHwmZmNc87NzHQsEjHnYOVKXzObNw9mzPBN3dOmwdq129erVctfltWxo0/WujQrd9SoAV26wJ57wnffwYIFftjXhx7yf3XqQM+evl9B166+ybxjR2jdWj33RUKUjf+mnsC3zrl5AGb2PDAIqDhZFxdvH7s4/rxZ2XNoubKsOsa1ZYu/1vmnn/y07OP162HNGli92k/XrIHly/17V546daBJEz8gR+fOuid1riso8El7n338e/v997BunX9///1v/xevZk1o08Yn7aZNfb+D2F/9+v79rlNn57+CAv8DIZm/RH/cVZf1ZJeUjWTdDlgY93wRcHilW8ya5c9fSv6pWdPXnAsK/Jdx/frQrJlP1LGe3CUlsHFj8mWXloYb667KOf/6J/p61q3rO6CB7wi4apVvKSku9j/qtmzx7+miRf5PZBdWD+qEUU42knV5Px936mZqZkOBoQAFUNoF1kcdWKathdq7weZsx1GWlftwhxkW93yndX6EWk1hSymUlpSUlFBS4t/fDRt8zXvhwrKbpBqnbYJNZrY0lAKrVgisSmYDy4PaUhPYc5+PP8ZBmL9+rAbUqOn7pZTL7fh/X1Ffc7fTgwStzdH/r3St1XHllRLoEkY52UjWi4AOcc/bA0vKruScGwGMADCzqXNCuKg815jZ1OXV9LiWVtPjCmNwg1xjZlPXVtPjqq7/Xzqu/GFmoYzolY1Ltz4D9jaz3c2sNnAOMC4LcYiIiOSFjNesnXNbzey3wJv4JrKnnHNfZzoOERGRfJGVayucc28AbySxyYioYskyHVd+0XHlFx1XftFxVSIvbuQhIiKyK9NwoyIiIjkuq8nazE40szlm9q2Z3VTO8n3NbIqZ/WxmNySzbTalelxm1sHMJpnZLDP72syuyWzklUvn/QqW1zSzz81sfGYiTkyan8PdzGysmc0O3rcjMhd55dI8rt8Fn8GvzOw5M8uZUWsSOK5fm9mXwd9HZtYt0W2zKdXjqgbfGxW+X8HyfP3eqOxzmPz3hnMuK3/4zmXfAXsAtYEvgP3LrNMSOAy4C7ghmW3z9LjaAIcEjxsB31SH44pbfh3wL2B8to8nrOMCRgGXBI9rA7tl+5hC+By2A+YD9YLnY4DB2T6mJI6rN9A0eHwS8Emi2+bpceX790a5xxW3PF+/Nyo8rlS+N7JZs9427KhzbjMQG3Z0G+fcCufcZ8CWZLfNopSPyzm31Dk3PXi8HpiF/+LMBem8X5hZe+BXwBOZCDYJKR+XmTUGjgGeDNbb7Jxbm5Goq5bW+4XvfFrPzAqA+pQzFkKWJHJcHznnfgyefowfyyGhbbMo5eOqBt8bFb1f+f69Ue5xpfq9kc1kXd6wo4l+wNLZNmqhxGZmnYGDgU/CCStt6R7Xg8D/EO4oWWFI57j2AFYCTwfNdE+YWYOwA0xRysflnFsM/AVYACwF1jnn3go9wtQke1xDgAkpbptJ6RzXNtXge6PscT1I9fjeiD+ulL43spmsExp2NIJto5Z2bGbWEHgRuNY5VxRKVOlL+bjMrD+wwjk3LdyQQpHO+1UAHAI84pw7GNgI5Mp50HTer6b4WsLuQFuggZmdF2Js6Uj4uMzsWPyX5I3JbpsF6RxXbH5ef2+UPa7q8r1RzvuV0vdGNpN1QsOORrBt1NKKzcxq4f/h/umceynk2NKRznEdCQw0s+/xzUXHmdmz4YaXsnQ/h4ucc7FazFj8P2EuSOe4+gHznXMrnXNbgJfw599yQULHZWZd8U2ng5xzq5PZNkvSOa68/96o4Ljy/nujks9h0t8b2UzW6Qw7mstDlqYcm5kZ/jzGLOfcAxHGmIqUj8s5d7Nzrr1zrnOw3bvOuVypqaVzXMuAhWYWG6j/eCq71WtmpfM/sgDoZWb1g8/k8fjzoLmgyuMys474HxjnO+e+SWbbLEr5uPL9e6Oi48r3741Kjiu1742oe81V9gecjO+5+B3wh2De5cDlwePW+F8hRcDa4HHjirbNlb9Ujws4Ct+U8iUwI/g7OdvHE8b7FVdGX3KoV2cIn8PuwNTgPXuFoPdnLvyleVx/AmYDXwHPAHWyfTxJHNcTwI9x/0NTK9s2V/5SPa5q8L1R4fsVV0Y+fm9U9jlM+ntDI5iJiIjkOI1gJiIikuOUrEVERHKckrWIiEiOU7IWERHJcUrWIiIiOU7JWiTPmdlkM/tlmXnXmtnwStbvkZnoRCQMStYi+e85/KAM8c4J5otINaBkLZL/xgL9zawObLuZQ1vgv8xsanCP4z+Vt6GZbYh7fIaZjQwetzCzF83ss+DvyGB+HzObEfx9bmaNIj42EcEPKC4iecw5t9rMPgVOBF7F16pHA/c459aYWU1gopl1dc59mWCxDwF/dc59EAyb+CawH3AD8Bvn3IfBjSN+Cv2ARGQnqlmLVA/xTeGxJvCzzGw68DlwALB/EuX1A/5uZjPwYx43DmrRHwIPmNnVwG7Oua0hxS8ilVCyFqkeXgGON7NDgHr4MYlvAI53znUFXgfqlrNd/HjD8ctrAEc457oHf+2cc+udc/cClwT7+NjM9o3gWESkDCVrkWrAObcBmAw8ha9VN8bfJ3edmbUCTqpg0+Vmtp+Z1QBOjZv/FvDb2BMz6x5M93TO/cc5dx/+RgRK1iIZoGQtUn08B3QDnnfOfYFv/v4an8A/rGCbm4DxwLvA0rj5VwM9zOxLM5uJv5sQwLVm9pWZfQFsAiaEfxgiUpbuuiUiIpLjVLMWERHJcUrWIiIiOU7JWkREJMcpWYuIiOQ4JWsREZEcp2QtIiKS45SsRUREcpyStYiISI77/0nMTzIKHjsPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xval = np.linspace(0.1,0.26,1000)\n",
    "tpdf = stats.t.pdf(xval,loc=smean,df=dof,scale=SE)\n",
    "plt.plot(xval,tpdf,color='red',lw=2); plt.xlabel('Values'); plt.ylabel('Density'); plt.title('Confidence Interval for the Population Mean')\n",
    "plt.xlim([0.1,0.26]); plt.ylim([0,80])\n",
    "\n",
    "plt.annotate('lower interval = ' + str(lower_CI),[lower_CI-0.0015,stats.t.pdf(lower_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+5],rotation=90.0)\n",
    "plt.annotate('upper interval = ' + str(upper_CI),[upper_CI-0.0015,stats.t.pdf(upper_CI,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+6],rotation=90.0)\n",
    "plt.annotate('$\\overline{x} = $' + str(np.round(smean,3)),[smean-0.0015,stats.t.pdf(smean,df = dof, loc = smean, scale = sstdev/np.sqrt(n))+2],rotation=90.0)\n",
    "    \n",
    "plt.fill_between(xval,tpdf,where= xval < lower_CI, color='red',alpha=0.7)\n",
    "plt.fill_between(xval,tpdf,where= xval > upper_CI, color='red',alpha=0.7)\n",
    "plt.vlines(smean,0,65,color='black',ls='--')\n",
    "plt.vlines(lower_CI,0,7,color='black',ls='--'); plt.vlines(upper_CI,0,7,color='black',ls='--')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.1, wspace=0.3, hspace=0.4); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Interval for the Proportion\n",
    "\n",
    "Now I demonstrate the calculation of the confidence interval for a proportion.\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{p} \\pm t_{\\frac{\\alpha}{2},n-1} \\times \\frac {\\hat{p}(1-\\hat{p})}{\\sqrt{n}} \n",
    "\\end{equation}\n",
    "\n",
    "First we need to make a new categorical feature to demonstrate confidence intervals for proportions.\n",
    "\n",
    "* we will assign samples as category 'high porosity' when their porosity is greater than 18%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New binary feature, 1 = high porosity, 0 = low porosity\n",
      "[1 0 0 1 1 1 0 0 0 0 0 0 1 0 0 0 0 1 1 0]\n",
      "\n",
      "Proportion of high porosity rock in well X1 is 0.35\n"
     ]
    }
   ],
   "source": [
    "CX1 = np.where(X1 < 0.18,0,1)\n",
    "prop_CX1 = np.sum(CX1 == 1)/CX1.shape[0]\n",
    "print('New binary feature, 1 = high porosity, 0 = low porosity')\n",
    "print(CX1)\n",
    "\n",
    "print('\\nProportion of high porosity rock in well X1 is ' + str(np.round(prop_CX1,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to calculate the confidence interval in the proportion with the SciPy stats t interval function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The confidence interval for the X1 proportion of high porosity is [0.244 0.456]\n"
     ]
    }
   ],
   "source": [
    "SE_prop = (prop_CX1*(1-prop_CX1))/math.sqrt(CX1.shape[0])\n",
    "ci_prop_95_x1 = stats.t.interval(0.95, CX1.shape[0]-1, loc=prop_CX1, scale=SE_prop)\n",
    "print('The confidence interval for the X1 proportion of high porosity is ' + str(np.round(ci_prop_95_x1,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence Interval for the Proportion By-Hand\n",
    "\n",
    "Let's repeat this calculation by-hand without the SciPy stats function to test our knowledge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistic +/- t-score x Standard Error\n",
      "   0.35   +/-   2.09  x     0.0509\n",
      "The confidence interval for the X1 proportion of high porosity is [0.244 0.456]\n"
     ]
    }
   ],
   "source": [
    "sample_prop = np.sum(CX1 == 1)/CX1.shape[0]\n",
    "tscore = stats.t.ppf([0.025,0.975], df = len(df)-1)\n",
    "SE_prop = (prop_CX1*(1-prop_CX1))/math.sqrt(CX1.shape[0])\n",
    "ci_prop_BH_95_x1 = sample_prop + tscore * SE_prop\n",
    "print('Statistic +/- t-score x Standard Error')\n",
    "print('   ' + str(np.round(sample_prop,3)) + '   +/-   ' + str(np.round(tscore[1],2)) + '  x     ' + str(np.round(SE_prop,4)))\n",
    "print('The confidence interval for the X1 proportion of high porosity is ' + str(np.round(ci_prop_BH_95_x1,3)))"
   ]
  },
  {
   "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",
    "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.\n",
    "\n",
    "##### Pooled Variance t-test Difference in Means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t statistic is -2.98 and the p-value is 0.00499\n"
     ]
    }
   ],
   "source": [
    "t_pooled, p_pooled = stats.ttest_ind(X1,X2) # assuminng equal variance\n",
    "print('The t statistic is ' + str(np.round(t_pooled,2)) + ' and the p-value is ' + str(np.round(p_pooled,5)))"
   ]
  },
  {
   "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 probability 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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t critical lower and upper values are [-2.02  2.02]\n"
     ]
    }
   ],
   "source": [
    "n1 = len(X1); n2 = len(X2)\n",
    "smean1 = np.mean(X1); smean2 = np.mean(X2)\n",
    "sstdev1 = np.std(X1); sstdev2 = np.std(X2,)\n",
    "dof = len(X1)+len(X2)-2\n",
    "t_critical = np.round(stats.t.ppf([0.025,0.975], df=dof),2)\n",
    "print('The t critical 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-critical 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 plot this result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAF6CAYAAADMN/v3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABOnklEQVR4nO3debzUY//H8denU2m1VUqLsrRZu7srwn3LnuyEECKyhXRn7+a+xc8aZUuElDXZU5Yi3OJWEm4qQirZSmnVev3+uGac6TinM+c033PNfOf9fDzmMd8z852Zz8yZ6/uZ6/peiznnEBERkexQKXQAIiIiUkiJWUREJIsoMYuIiGQRJWYREZEsosQsIiKSRZSYRUREsogSc44zs2Zm5syscoafd5yZnZHJ58xFZtbDzP4T1WPN7D0z+0v5ossMM+tkZvMq6LW2M7NlZlZQEa9XQgzPmVnnDDzPMjPbIbFd3cxeNrPfzOyZxG03mNkCM/txU18rBDP7m5nNDB1HPlJiLgczm21mBxW5rdwH8E197Sg45w5zzj2a7v6JQrwscVme+LGwLOWyXVljSDzHTmXYf7iZrU683q9m9oaZtSrr61YUMzsSWOqc+zjx97/MbE2Rz21x4BhnmNlZxdx+iZlNKevzOefmOOdqOefWZSbCcrkZuLGkOxM/VNan/A/mmdkoM2uful/ifXyT+LMrUB+o45w7wcyaAP8AdnbONYjqjUTJOfeuc65leR6bOB46M7ujyO3HJG4fnpEgY0qJWTIiUYhrOedqAbskbt4yeZtzbk4FhXJrIobGwM/A8Ap63fI4DxhZ5LanUz6zWs65LQPElepR4PRibj8tcV/aMt2qU17OuQ+Bzc2s3UZ2m5/4HtUG9gJmAO+a2YEl7N8U+NI5tzbl74XOuZ/LGp95cTg2fw2cVOT/fjrwZaB4ckYc/vlZx8wuM7Nni9x2t5kNSmxPNLObzOzDRNPXi2a2dcq+R5nZ52a2OLFv68TtI4HtgJcTv+QvT3mJU81sTqLp7JqU56pkZlea2ddmtjDxy3/rxH3VzOyxxO2LzWyymdVPifHsxPZOZvZ2ItYFZvZ0GT+PLczsITP7wcy+TzTxFWzsuc3sncTDP0m815PK8prOuRXAE8CuiedrnXhPixOf7VFF4hthZr+Y2Xdm1r+kA6OZtUrUxH81s5lmdmLKfXXM7CUzW2JmHwI7buQzqQocALyd7nsys4MTNdjfzOyexOeW/B/9y8weS9l3g1McZnammU03s6Vm9o2ZnZvmy44E9jWzpinP3RrYHXjSzA43s48T73mumf2rmBh6mtkc4M2yxGWJJnYz+4eZ/Zz4/pyZcn91MxuY+J/9Zmb/MbPqifv2MrNJif/3J2bWqcj7mggcXtqbd94859y1wDDglpTXd4nv77+Ba/FJaFniPbwBNEz8Pby0mBLfzRvN7D1gBbBDKd+14WZ2r5m9kvjs/mtmO6bcv0vKY38ys6sTt5d4PCjKipziMN9a18/MPk183k+bWbWNfHw/Ap8BhyYevzWwN/BSkdfZ2OeyKd+PLmb2ReKx35tZv43Eml2cc7qU8QLMBg4qclsP4D+J7W2B5fgaI0BlfO3tr4m/JwLf45NGTeBZ4LHEfS0Sjz0YqAJcDswCqhb32kAzwAEPAtWBPYBVQOvE/X2AD/A1yM2AocCTifvOBV4GagAFwF+BzVNiPDux/SRwDf6HXDVg31I+n2RMlRN/v5B43ZrANsCHwLmlPXfiOXYqw/9lOHBDYrsWPjG/m/gcZwFXA8mEuBRomdh3BPAivnbUDP+Lvmcx/9eawFzgzMT/tC2wANglcf9TwKjEfrsm/sf/KSHWXYDlRW77V/J7UMz+dYEl+CbTKsClwNqU/9EGjy3mf3A4/oeCAfvhD/5tE/d1AuZt5HN9A+if8vdNwAspj90t8f/bHfgJOKZIDCMSn0n1csS1Frg+8Z67JO7fKnH/vfjvaSP893dv/He8EbAwsX8lfFlaCNRLeQ99gedKeL/Ffh747816oGbR72cxn/8Gz1FaTIn3MSfxvagMbMHGv2vDgV+BDon7HweeStxXG/gB35ReLfH3nqUdD0r7HPDHng+BhsDWwHTgvBIe2wP4D3AKvhUI4ILE690ADE/zc9mU78cPwN8S21slH5cLl+AB5OIl8QVdBixOuawg5SAMjAPOSWwfAXyRct9E4OaUv3cGVuMPLv8ERqXcVwl/gO+U8trFJebGKbd9CHRLbE8HDky5b1tgTaIwnwVMAnYv5j1OpPCgPwJ4IPU1Svl8kjFVxp93WwVUT7n/ZOCt0p6b8iXm3xP/jx/xv8x3BP6W+LtSyr5P4g+mBYn4dk6571xgYmK7B4WJ+STg3SKvORS4LvE8a4BWKff9HyUn5n2AH4vc9q/E9yD1e5X8nE4HPkjZ14B5pJmYi3n9F4BLEtud2Hhi7g7MTPk+zgGOLWHfQcCdRWLYYRPiWpm6L/4H7l6JOFYCexTzHFcAI4vc9hpwRsrf5wBvlhBDsZ8H0CoRe6Oi389iPv8NnqO0mPDl7fqU+0r8rqV814el3NcFmJFSvj4u4b2VeDwo7XPAH3u6p/x9K3B/Ca/TA5+Yq+N/rG2B/0GwDxsm5lL/V+X5fiS25+DL8uYlfbez9aKm7PI7xjm3ZfKC/zWY6lH8AY3EddFziXNTtr/D/+Kri/81+l3yDufc+sS+jUqJJ7Xn5wp8jRH8ua7nE81Ei/EFcx0+YY7EF4KnzGy+md1qZlWKee7L8YngQ/PNwH/qDLQRTRPv7YeUGIbia86b+tzFuT3xP2ngnDvKOfc1/jOdm/gsk77Df6Z18bXo74q5r7j3smfyfSTey6lAA6Ae/odI0f9rSRbhazJFjUr9Xjnn9k/c3jD1uZ0/8swt5vHFMrPDzOyDRNPmYvyBvG6aD38O2NbM9sIfDGsArySed08ze8v8aYDf8OfNiz5viXGmEddCV3jeFgq/23XxtcGvi3napsAJRf5P++KTUFJt/A+fsmiET8ZlfVy6Mc0tsn9J37Wkksp8E4r/XJLPW9LxIB0lvWaxnHMr8d+V/kBd59x7xcRT4ueyCd8PgOMT+39n/rRPxzTfY3BZ0Rkjpl4AhpjZrvga8+VF7m+Ssr0d/lfrAmA+vmkQ8B1BEvt+n7jJlTGOucBZxRSIpH8D/zazZsBYYCbwUOoOzrkf8TUMzGxfYLyZveOcm5Xm66/CF8q1Re/cxOdO13ygiZlVSknO2+GbrBfgP/umwBcp933/p2fx7+Vt59zBRe8wf858Lf5/NSPleUrylX+YNXLOFfdaRf1Ayncm5XuRtByfMJMapOy7Gf50yenAi865NWb2Av4HUamccyvMbHTi8dXxTaarE3c/AdwDHOac+918P4qiibnY7+wmxrUA3zqyI/BJkfvm4mth52zk8a2LeVxpjgWmOueWl/Fx6caU+jmV+F1L87VO3sh9GzseRGEE8Cb+WFNcPMV+Lhn43k4Gjk5UNnrjTzM12fijsoNqzBFxzv0OjMYfuD50f+6V3N3MdjazGvhzJKOdH0IyCjjczA5MfKH+gU9skxKP+wnYoQyh3A/caInOO2ZWz8yOTmzvb2a7JZLKEnyC+tMwFjM7wcwaJ/5chD+ApDXcxTn3A/A6MNDMNk90PtnRzPZL47n/9F7Nd7jplN5b/8N/8YnrcjOrknj8kfgEk/zMbzSz2onPqS/wWDHPMwZoYWanJZ6nipm1N7PWied5DviXmdUws52BMzbyuawBxuPPm6XjFWAXMzvOfMepi9mw9jQN+Lv5ccJbAFel3FcVfz7xF2CtmR0GHJLm6yY9im9ePZ4Ne2PXBn5NJOUO+HOK6Sp3XIkfWA8Dd5hZQzMrMLOOiYP5Y8CRZnZo4vZq5jsKNU55iv3wp5s2yrxGZnYdcDa+n0J5pBNTqhK/a2m81higgZn1MbPNEt/rPRP3lXg8iNDb+HPHdxdz38Y+l3J/P8ysqpmdamZbJMraEtI8ZmUDJeZoPYqv/RZtxiZx23B801A1/IEW59xMfNP33fhawZHAkSk1lJuA/olmn3R6GQ7Gn2t93cyW4s/zJAtpA/yPhyX4Jq23KT4htQf+a2bLEs91iXPu2zReO+l0fCH7Ap98R1PYhLex5/4X8GjivZ6YKKzL8D0905b47I4CDsN/pvcBpzvnkjXbi/CJ+xv8ebEn8Af9os+zFH9g6Iavhf+I76W7WWKX3vhmtB/x/9tHSgltKH7YUapkz97UyzbOuQXACfgxuAuB5sAftR7n3BvA08CnwEf4g3Nq3Bfjf4AswifPDXrGpuEd4Dfg+0RNJOkC4PrEd+vaxGukJQNx9cN/FybjO0Ldgu9HMBc4Gp9Ef8HXyi4jcbwzPx55ufPDpkrSMPGdXJZ4/t3w/TxeL0N8fygtpmL2L+27trHXWopPhEcmHvcVkDwlsrHjQSScN8E592sx95X4uWTg+3EaMNvMluBPsXQvZf+sYYmT5BIB85NqzAAaOOeWpNw+Ed9RZFio2HKRmXXH90q9qtSdc4T5SWkucolJRsr42Inoe1Rm5ocyPuScGxs6FpHi6BxzRMyPg+2Lby5dUtr+UjrnXHG1+ZzmnNs3dAz5xjl3fOgYRDYm0qZsM+tsfmD8LDO7soR9OpnZNPM9ctOebCGbmVlNfPPwwfihNCIiImmJrCk70aHoS3xymoc/T3Oyc+6LlH22xHdq6uycm5M4l1bmKexERETiIsoacwdglnPum0Tnm6fwJ/lTnYKffWcOgJKyiIjkuygTcyM2HDA/jz9P2tAC2Mr8PLEfmVlxk+WLiIjkjSg7fxU3CLxou3ll/PzMB+InLnjfzD5wzm2w+oiZ9QJ6AdSsWfOvrVpl7Up+5Tczsexpy3KtsiYiIlnuo48+WuCcq1faflEm5nlsOMtKY/x4vKL7LEjMpLPc/IpCe1BkWTDn3AP4+ZRp166dmzKlzMvAZr9Onfz1xIkho8iY22+/HYB+/XJnQReRXKNyllvMbGPT9P4hyqbsyUBzM9ve/BJ33fjz4PAXgb+ZWWXzM2DtiZ/oQnLcmDFjGDNmTOk7iki5qZzFU2Q1ZufcWjPrjV8koQB42Dn3uZmdl7j/fufcdDN7FT9b0Xr8ain/iyomERGRbBfpBCOJmXXGFrnt/iJ/3wbcFmUcIiIiuUJzZYuIiGQRTckpkahevXroEERiT+UsnpSYJRLjxpW6op6IbCKVs3hSU7aIiEgWUWKWSAwYMIABAwaEDkMk1lTO4kmJWSIxYcIEJkyYEDoMkVhTOYsnJWYREZEsosQsIiKSRZSYRUREsoiGS0kk6tSpEzoEkdhTOYsnJWaJxLPPPhs6BJHYUzmLJzVli4iIZBElZonEVVddxVVXXRU6DJFYUzmLJzVlSyTef//90CGIxJ7KWTypxiwiIpJFlJhFRESyiBKziIhIFtE5ZolE48aNQ4cgEnsqZ/GkxCyReOyxx0KHIBJ7KmfxpKZsERGRLKLELJHo06cPffr0CR2GSKypnMWTmrIlEtOmTQsdgkjsqZzFk2rMIiIiWUSJWUREJIsoMYuIiGQRnWOWSLRo0SJ0CCKxp3IWT0rMEokHHnggdAgisadyFk9qyhYREckiSswSiV69etGrV6/QYYjEmspZPKkpWyLx5Zdfhg5BJPZUzuJJNWYREZEsosQsIiKSRZSYpdxmzJjBhAkTWLZs2Qa3v/rqq4EiEhHJfUrMUi533XUXRx99NHfffTe77rorL7744h/3XX311bRp04Y2bdqEC1AkD6icxZM6f2WDRx4JHUGZPfjgg3z00UfUqlWL2bNn07VrV2bPns0ll1yCc45BgwaFDlEk9lTO4kk15mxw3XWhIyizdevWUatWLQCaNWvGxIkTGTduHH379sU5Fzg6kfjQKaP8o8RcUXbfvfjLbrvBTz+Fjq7MGjRosMGSc7Vq1WLMmDEsWLCAzz77jO7du9O9e/dwAYrEQGmnjFTO4klN2RXlp5/gtddgq602vN052HvvMDFtghEjRlC58oZfn8qVKzNixAjOPfdcrrnmmkCRicRHaaeM5s2bFzpEiYASc0U54ghYtgyK66jRqRPMn1/REW2Sxo0bM2PGDB599FG+//57zIyGDRty1FFHsc8++4QOTyQWijtl1LVrV7777judMooxNWVXlIcegn33Lf6+J56o2Fgy4JZbbqFbt2445+jQoQPt27fHOcfJJ5/MzTffHDo8kVgo7ZSRxJNqzFIuDz30EJ9//jlVqlTZ4Pa+ffuyyy670KhRo0CRicSHThnlJyXmirJ6NTz1FDRsCAcd5GvJkyZB69aQg5PQV6pUifnz59O0adMNbv/hhx+oVKkSHTt2DBSZSHw0bty4xPv22WcflbOYslw7T9GuXTs3ZcqU0GGU3amnwtq1sGIFbLmlP9983HEwYYLvAPbdd36/iRNDRpm2V199ld69e9O8eXOaNGkCwJw5c5g1axb33HMPnTt3DhyhSDzMmDGDF1988U99OVq3bh06NCkjM/vIOdeu1P2UmCvI7rvDp5/65Nyoke/sVVDgk/Iee8DWW/v9ciQxA6xfv54PP/yQ77//HuccjRs3pn379hQUFIQOTSQWbrnlFp588km6dev2R+153rx5PPXUU3Tr1o0rr7wycIRSFukm5kibss2sMzAYKACGOeduLnJ/J+BF4NvETc85566PMqZg1q/3zdnLl/ta82+/+WS8ahWsWRM6unKpVKkSe+21159u//HHH7nwwgsBePbZZys6LJHYKK0vx+TJkwGVs7iJLDGbWQFwL3AwMA+YbGYvOee+KLLru865I6KKI2v07AmtWsG6dXDjjXDCCbDDDvDBB9CtG7z1VugIM6Znz54sX748dBgiOa+0vhwLFy4MFJlEKcoacwdglnPuGwAzewo4GiiamPPDpZfCSSf57YYN4fTTYfx4OOcc6NAhVon5lVdeoVOnTqHDEMl5gwYN4sADDyyxL4eGJsZTlIm5ETA35e95wJ7F7NfRzD4B5gP9nHOfRxhTWCtXwvPPw9y5ULkyNG8OBx8cOioRyVKdO3fmyy+/LLEvhxJzPEWZmK2Y24r2NJsKNHXOLTOzLsALQPM/PZFZL6AXwHbbbZfhMCvI4MEwZgzstx9MnuxnAJs7Fzp2hPvuCx1ducyYMYPvv/+ePffc84/ZiUCT64tk0sb6ckg8RZmY5wFNUv5ujK8V/8E5tyRle6yZ3WdmdZ1zC4rs9wDwAPhe2dGFHKFhw2DaNN8Tu29f6NLF98A+91w4+mjYYovQEZbJXXfdxb333kvr1q3p2bMngwcP5uijjwb85PrHHnts4AhF4q1nz54ceOCBocOQCESZmCcDzc1se+B7oBtwSuoOZtYA+Mk558ysA36K0Pj2Zli71ifmVatg6VJ/23bb5WSv7NIm1//nP/8ZOkSRWHvllVdChyARiSwxO+fWmllv4DX8cKmHnXOfm9l5ifvvB7oC55vZWmAl0M3l2sDqdJ19NrRvD3vtBe+8A1dc4W//5ZfCMcw5RJPri1SMjZ0y0kQ+8aQJRirS55/D9Omw665+6FSqZC/mHJlg5IADDuCOO+6gTcpqWWvXruWss87i8ccf55BDDgFg3LhxgSIUyX2pp4ymTZu2wSmjtm3bUr9+fUDlLFdkxQQjUsQuu/hLUTnYiaO4yfUXLFigyfVFMqi0U0YrV64MHaJEQIk5G/TsGTqCMitucv0uXbowdepUrccskiE6ZZSftB5zNohJJw4dKEQyS+sx5yclZsmYc845J3QIIrEyYsQIGjRosMFtyfWY33nnnUBRSdTUlF2RZsyA77+HPfeElN6VxGRCjgsuuOCP7SOOiP/05yJRK209ZpWzeFKv7Ipy111w773QurWfaGTwYD+xCEDbtrD55n47R3pli4hI2ahXdrZ58EH46CNfU549G7p29deXXOLXZBYREUHnmCvOunWFzdfNmvma8bhxfnrOHE/Mb7755gbXAJ06ddIKUyIZpHKWP5SYK0qDBr4JO6lWLb+oxYIFkOO9K/v167fBtYhknspZ/lBirigjRvjknGrBAn97THpX5lp/BZFcpHIWf0rMFaVx4z8n5i5d/LUm5BARkQQl5pD0y1dERIpQr+yQYjwhx4knnhg6BJHYUzmLJyXmkFIm5Mhlybl8a9eu/cdtF8TkvYlkC5Wz/KEJRrJFji37WJoVK1YAUKNGjcCRiMSXyllu0QQjElSXRMe2iTH5oSGSjVTO4kmdv0JIThCQMlGAiIgIqMYcRr9+MHVq4XUMnHnmmdSqVYu2bdvSvn17nHOYWeiwRGJF5Sw/KDGHlGPn9zfmkUceYeXKlUydOpXx48fz1Vdf0aJFi9BhicSKyll+UGKWjLnlllv473//S8OGDalZs2bocERiSeUs/nSOWTJm8eLF7LXXXvTv359mzZrRo0eP0CGJxI7KWfwpMUvGbL311qxbt45tttmGffbZRwcMkQionMWfmrJDSC7/mDJRQBxce+21zJ8/n4svvpjtt9+eBQsWULdu3dBhicSKyln8qcYcQnI1qZisKpV0zz330LBhQx566CFeffVVunbtGjokkdhROYs/JWbJmO++++6P7Tlz5gSMRCS+VM7iT4k5lMcfhy+/DB1FRq1fv553332X9evXs2bNmtDhiMSSyln8KTGHUq+eX8Ti73+H446Db74JHdEmu/XWW/nkk08455xzqFOnTuhwRGJJ5Sz+1PkrlEMOgbffhhtvhKVLoVWr0BFtkmuvvZZ169bRpk0brrzySs6J8ZKWIqGonOUHJeaQlizxU3LuuiusWxc6mk1y/fXX89NPP/Hxxx/z7LPPUlBQQK9evUKHJRIrKmf5Qcs+hrRmDQwZAtOmwQcfwDbb5PSyj3fccQd9+/YFYMaMGbTK8VYAkWykcpa7tOxjLqhSBS6+2G8n12POQYsXL+bSSy9l5syZVKtWjT322IPBgwczcOBAmjRpEjo8kVhQOcsf6vwV0oUXQnLWnl9/DRrKpthyyy1p3Lgxv//+O8uXL+eLL75g+vTpnHbaaaFDE4kNlbP8ocQcUtWqsMMOfnvx4qChbKoBAwYwduxYdtttNxYuXMiCBQtChyQSOypn+UGJOaQaNeC33/y55lWrQkezyUaMGEHnzp258soradmyZehwRGJJ5Sz+dI45pH//Gx54wDdpb7NN6Gg22YsvvkiNGjU4+OCDQ4ciElsqZ/GnGnMoZ54Jl14K1avDRRfB1luHjmiTPffcc+y00048//zzzJw5M3Q4IrGkchZ/qjGH8sgjsHIlfPwxvPkmfPUVtGgROqpN8swzz9C7d286d+5M06ZNqZVcRUtEMkblLP5UYw7pppvgzjt9cm7cOHQ0myx1cv1p06Zx5JFHBoxGJJ5UzuJPiTmkX3+FZ57x55m//z50NJssdXL9b775Rs1sIhFQOYs/JeaQNtvMT8lplrNTcqb+ek+dXH/69Omce+65ASMTiQ+Vs/yixBzSjTfChAnQq5dfbSoHHXvssX9sFxQU0Lt3b8455xzq1q0bMCqReFE5yy/q/BXS//0fTJ/uxzNXrx46mjIZNWoUU6dOZenSpUyfPp0WLVpQUFAAQK9evdg6Br3MRUJTOctPSswhJc8xr1kDzZpB8+ahI0rbPvvsw++//86wYcPo27cvM2fOZMstt6Rhw4ZUz7EfGSLZSuUsPykxh5Q8x7zbbjl3jrlRo0acfvrp7Ljjjuyzzz4A/Prrr3z77be0atWKww8/PHCEIrlP5Sw/Rbrso5l1BgYDBcAw59zNJezXHvgAOMk5N3pjzxmrZR9//x3uvhu++AKmTIE6dXJ22ceXX355g2Eb48ePB+Cggw4KFZJI7LzxxhsbzPilcpZb0l32MbLEbGYFwJfAwcA8YDJwsnPui2L2ewP4HXg4rxIzwMKFPiEnl33M0cS8++678+mnn4YOQyTW2rZty9SpU0OHIeWUDesxdwBmOee+SQT0FHA08EWR/S4CngXaRxhLdjr+eKhf3y9ksWgRbLVV6IjKregPvGnTpgHQpk2big9GJE+onMVTlIm5ETA35e95wJ6pO5hZI+BY4ADyMTG3auWHTAE0apTTidnMNvi7T58+AEzM0RYAkWxx5plnYmY455gzZw5nnXXWH/d98803gMpZ3ESZmK2Y24q2mw8CrnDOrSt6YN/gicx6Ab0Atttuu0zFF8Z330HTpn771Vf94hV77OEnGRGpaKtW+SF706f7UQKrVvkfiE2aQNu2sVhcJdf16NHjj+3//Oc/nHHGGX/8fd111wWISKIWZWKeBzRJ+bsxML/IPu2ApxJJuS7QxczWOudeSN3JOfcA8AD4c8xRBVwhjj3W98QGn5gnT4bbb/dDpmbMCBub5Ic1a+CFF+DJJ+H112H58pL33WMPOOEE6NHDt+pIhdtvv/3+2K5du/YGf0s8RZmYJwPNzWx74HugG3BK6g7Oue2T22Y2HBhTNCnHxqhRPiEvXeprJy1a+Nm+unSBK6+E1q1DR7hJ6tevHzoEKc3atTBsmD99Mm9e4e3Nm/she/XrQ9Wqvr/DV1/5xVU++cRf/v1vOPVUf53rrVY5rGrVqqFDkAoQWWJ2zq01s97Aa/jhUg875z43s/MS998f1WtnpX328cOjhg2Dvn1h5kzYckto2DDnZv0qzhtvvBE6BNmYDz6As87yPwoBWraE88+H447zzdbFWbUK3ngDHn0UnnsOhg+Hp5+Gq6/2PyYraxqEivbBBx+EDkEqQKTjmKOQ88Ol3nvPJ2nw5/S+/dZ3AktOFBCTThyTJk0CYO+99w4cSZ5bs8bXcm+6Cdavhx13hJtv9gm5Uhmmyv/6a7jmGp+YATp2hMcfh+233/jjJFIqZ7kl+DjmqOR8Yi5Jjo9jliy0aJE/Pzxhgu9ceNllcP31fsa58powAc44wy9TuvnmfkrZQw7JXMwiMZZuYtbqUhKJSZMm/fFrXgL4+mvYay+fSOvX9z/4brll05IywIEH+nPORx8NS5b4PhIPPZSRkKXsVM7iSSeJJBJXX301oPGVQcycCQccAPPn+17VL72U2Q5bder4c85XX+2T/dln+9Myl12WudeQtKicxZNqzCJxMn26Py0yfz7stx+8+240vagrVfLnqocM8X9ffjnccUfmX0ckDykxi8TF3Llw8MHw44++xvzKK1C7drSved558OCDfvsf/yhM1CJSbkrMInGwaBEcdpjvlLXvvvDyy1CzZsW89tlnw/2J0Y+9e8OYMRXzuiIxpcQskuvWrvULonz+Oey8sz+nXKNGxcZw7rlw3XV+SNZJJ8FHH1Xs64vEiDp/SSQGDRoUOoT8ccUV8NZb0KABjBsXbjGU666Db76BkSMLp56tWzdMLHlC5SyelJglElqGroI8/bTvdFW5MoweHXa6TDM/s91XX/mZxk49FcaOhYKCcDHFnMpZPKkpWyIxfvx4xo8fHzqMePv8cz/NJvjknJxRLqSqVf288HXr+gUyrr8+dESxpnIWT5r5K1vEbOavTon3o/GVEVm1Ctq3h88+8zXTkSOza+nQN96AQw8t3D7wwLDxxJTKWW7RzF8icda/v0/KO+7oe0RnU1IGP2zruuvAOb9k5OLFoSMSyRlKzCK55s03YeBAf+72scegVq3QERXvmmtgzz39EpMXXRQ6GpGcocQskkuWLPGLSDgH//ynnw87W1WuDCNG+GVNH3vMd04TkVIpMYvkkquv9jXQ9u19jTTbtWgBt93mt887DxYsCBuPSA7QcCmJxNChQ0OHED/vvw/33edrosOG+etccMEF8Oyzfqx1v34wfHjoiGJD5SyeVGOWSLRs2ZKWLVuGDiM+Vq+Gc87xTdiXXQa77x46ovSZ+Q5qm20Gjz7qz5FLRqicxZMSs0Ti5Zdf5uWXXw4dRnzcdpsft7zTTv7ccq5p0aKw6f288+D338PGExMqZ/GkxCyRGDhwIAMHDgwdRjx88w0MGOC3hw71naly0RVX+Lm8v/rKLxkpm0zlLJ6UmEWyXb9+fkKR7t39co65qmrVwmUhb70V5swJG49IllJiFslmEybA88/7JRzjUMv8+9/hxBNh5UpfgxaRP1FiFslWa9dCnz5+++qroVGjoOFkzK23QrVq8NRT8O67oaMRyTpKzCLZ6oEH4H//g+23h759Q0eTOU2bFtaWL7kE1q0LG49IlsmRgZCSa0aOHBk6hNy2aFFh7+uBA30NM04uvxweegg+/tiPa+7ZM3REOUnlLJ5UY5ZINGnShCZNmoQOI3fdeiv8+qtfdeyYY0JHk3k1asAtt/jt667z55ylzFTO4kmJWSLx9NNP8/TTT4cOIzfNnw+DB/vtW27JvpWjMqVbN/jLX+D77+Gee0JHk5NUzuJJiVkiMWTIEIYkh8ZI2Vx/va9BHn88dOgQOproVKoEN93kt2+6SUtDloPKWTwpMYtkky+/9PNgV6oEN9wQOproHXKIb65ftMg334uIErNIVunf3/dSPussaNUqdDTRMyscnz1okG/GF8lzSswi2WLKFHjmGd8D+7rrQkdTcfbcE447zjffJ6ceFcljSswi2eLaa/31RRdB48ZhY6loN97om++HDYNvvw0djUhQaY1jNrNngYeBcc659dGGJHEwevTo0CHklg8/hHHjoFYtP8Y337Rq5ecCHzEC/u//4MEHQ0eUE1TO4indGvMQ4BTgKzO72czy4OSXbIq6detSt27d0GHkjn//21/37g35+rn17+9rzcOHw+zZoaPJCSpn8ZRWYnbOjXfOnQq0BWYDb5jZJDM708yqRBmg5Kbhw4czfPjw0GHkhsmTYexYv1DFP/4ROppwmjeHU07xc4Qnh1HJRqmcxVPa55jNrA7QAzgb+BgYjE/Ub0QSmeQ0HTDK4Prr/fWFF+ZvbTkpWWt+5BH47rvQ0WQ9lbN4Sisxm9lzwLtADeBI59xRzrmnnXMXAbWiDFAk1j76CMaM8VNU9usXOprwWraEk0+GNWtUa5a8lW6NeZhzbmfn3E3OuR8AzGwzAOdcu8iiE4m7ZG35ggugXr2wsWSL/v39+OaHH4Y5c0JHI1Lh0k3MxU1B9H4mAxHJOx9/DC+9BNWrq7acqlUrP4/2mjWFk4+I5JGNJmYza2BmfwWqm9lfzKxt4tIJ36wtIuWVbKo9/3yoXz9sLNnmn//0teaHHoIffggdjUiFKm0c86H4Dl+NgTtSbl8KXB1RTBIDY8eODR1CdvvqKxg9GqpWze+e2CVp3RqOPRaee86vtKWac7FUzuJpozVm59yjzrn9gR7Ouf1TLkc5556roBglB9WoUYMaNdSoUqLbbwfn4LTToGHD0NFkpyuu8NdDhsBvv4WNJUupnMVTaU3Z3RObzcysb9FLBcQnOeq+++7jvvvuCx1GdvrhBz+JhhlcdlnoaLJXhw6w//6wZAncf3/oaLKSylk8ldb5q2biuhZQu5iLSLFGjRrFqFGjQoeRnQYPhtWr4Zhj/PAgKVmy1nznnfD772FjyUIqZ/G00XPMzrmhiet/V0w4IjH322++aRYKk46U7JBDoE0bmDbNz6Pdq1foiEQil+4EI7ea2eZmVsXMJpjZgpRm7o09rrOZzTSzWWZ2ZTH3H21mn5rZNDObYmb7ludNiOSMoUN902ynTn65Q9k4s8IfMLfd5teqFom5dMcxH+KcWwIcAcwDWgAbPTlmZgXAvcBhwM7AyWa2c5HdJgB7OOfaAGcBw9IPXSTH/P67b5IF1ZbLomtX2GEHmDXL99IWibl0E3NyoYouwJPOuV/TeEwHYJZz7hvn3GrgKeDo1B2cc8uccy7xZ03AIRJXI0fCjz/CHnvAoYeGjiZ3VK5cOAHLzTf73uwiMZZuYn7ZzGYA7YAJZlYPKK0nRiNgbsrf8xK3bcDMjk089yv4WvOfmFmvRFP3lF9++SXNkCWkiRMnMnHixNBhZI9163xTLPjaslnYeHJNjx6wzTYwdSq8+WboaLKGylk8pbvs45VAR6Cdc24NsJwitd9iFHfk+dNPXefc8865VsAxwIASXv8B51w751y7eppPWHLRyy/7SUWaNYMTTggdTe6pXh0uushvJ08HiMRU2ss+Aq2Bk8zsdKArcEgp+88DmqT83RiYX9LOzrl3gB3NLM/XvYuH22+/ndtvvz10GNlj0CB/3aePb5qVsjvvPKhWDV55BWbMCB1NVlA5i6d0e2WPBG4H9gXaJy6lrSo1GWhuZtubWVWgG/BSkefdycy36ZlZW6AqsLBM70Cy0pgxYxgzZkzoMLLDtGnw9ttQuzaceWboaHJX3bpwxhl+O/lDJ8+pnMVTuj/d2wE7p3TUKpVzbq2Z9QZeAwqAh51zn5vZeYn77weOB043szXASuCksryGSE646y5/feaZsPnmYWPJdX36+CFnjz4KN9zgk7VIzKTblP0/oEFZn9w5N9Y518I5t6Nz7sbEbfcnkjLOuVucc7s459o45zo65/5T1tcQyWq//AJPPOE7eyXPkUr5tWoFhx/uh55pmk6JqXQTc13gCzN7zcxeSl6iDEwkFoYOhVWrfDLZaafQ0cRD38Q0/ffc4z9bkZhJtyn7X1EGIfFTvXr10CGEt3o1JBcYuOSSsLHEyf77+7Hgn3wCTz7ph1LlKZWzeEorMTvn3jazpkBz59x4M6uBP28sUqxx48aFDiG80aP9SlK77AIHHhg6mvgw87XmM86AO+7w13k6LlzlLJ7S7ZV9DjAaGJq4qRHwQkQxicTD4MH++uKL8zZxRKZbN2jQAD77DCZMCB2NSEale475QmAfYAmAc+4rYJuogpLcN2DAAAYMKHa+mPzwwQfw4Yew9dbQvdT1XqSsqlaF3r399h13hI0loLwvZzGVbmJelZjvGgAzq4zmtZaNmDBhAhPyuSaTrC336gU1aoSNJa7OPdfPCDZuHHzxRehogsj7chZT6Sbmt83saqC6mR0MPAO8HF1YIjls3jx45hkoKIALLggdTXylTjhy991hYxHJoHQT85XAL8BnwLnAWKB/VEGJ5LT77vOLVhx/PDRpUvr+Un7JseEjRsCiRWFjEcmQdBexWI/v7HWBc66rc+5BzdAlUoyVK+GBB/y2hkhFb+ed4eCDYcUKeOih0NGIZMRGE7N5/zKzBcAMYKaZ/WJm11ZMeJKr6tSpQ506dUKHUfEefxwWLoR27aBjx9DR5IeLL/bX99zjWyrySN6Ws5grbRxzH3xv7PbOuW8BzGwHYIiZXeqc0/prUqxnn302dAgVz7nCTl+XXKIhUhWlSxfYcUf4+mu/vOYxx4SOqMLkZTnLA6U1ZZ8OnJxMygDOuW+A7on7RCTprbfgf//z42tPPDF0NPmjUqXCc83JBUNEclhpibmKc25B0Rudc78AVaIJSeLgqquu4qqrrgodRsVK1pbPP9+Ps5WK06MH1Krlfxx9+mnoaCpMXpazPFBaYl5dzvskz73//vu8//77ocOoOMlm1KpV/fhaqVhbbFE4Z3YeDZ3Ku3KWJ0pLzHuY2ZJiLkuB3SoiQJGccM89/hzzKadA/fqho8lPyebsxx7zHfBEctRGE7NzrsA5t3kxl9rOOTVliwAsXQoPP+y3NUQqnBYt4LDD/FrNw4aFjkak3NKdYERESjJ8OCxZAn//O7RpEzqa/JYcOnXvvbB2bdhYRMpJiVki0bhxYxo3bhw6jOitX1/YE1i15fAOOcTXnOfOhRdeCB1N5PKmnOUZy7UJvNq1a+emTJkSOozM69TJX0+cGDIKKasxY+DII6FpU98BrEDLlAd3zz3+fPPf/gbvvBM6GpE/mNlHzrl2pe2nGrPIpkgOkerdW0k5W5xxBmy+Obz7Lnz8cehoRMpMiVki0adPH/r06RM6jGh9/jmMH++XdezZM3Q0klS7Npx1lt+O+YQjeVHO8pASs0Ri2rRpTJs2LXQY0Uoe9M84A7baKmwssqELL/RToj7xBPz8c+hoIpMX5SwPKTGLlMevv8LIkX472RNYssdOO8Hhh8Pq1fDgg6GjESkTJWaR8njwQb/E46GHQqtWoaOR4iR7yd93H6xZEzYWkTJQYhYpq7Vr/ThZAJ3fy14HHgitW8P8+fDcc6GjEUmbErNEokWLFrRo0SJ0GNF4/nk/TrZlSz9uVrKTWeFphmTv+ZiJdTnLYxrHnC00jjl37LsvvPeerzVfcEHoaGRjli+Hxo1h8WL48ENo3z50RJLHNI5ZJApTpvikvMUWcLqWJM96NWvC2Wf77ZjWmiV+lJglEr169aJXr16hw8i85MH97LP9+r+S/Xr3hkqVYNQo+OGH0NFkVGzLWZ5TYpZIfPnll3z55Zehw8isH36Ap5/2B/nevUNHI+lq2hSOOcb3zB4yJHQ0GRXLciZKzCJpu/9+f3A/+mho1ix0NFIWyaFT99/vl4UUyWJKzCLpWLXKH9RBq0jlor/9zS/J+csv8NRToaMR2SglZpF0PPWUn9pxjz38usuSW8wKf1ANHgw5NhpF8osSs0SiTZs2tGnTJnQYmeFcYaevPn38QV5yT7duUK8eTJvmV56KgViVM/mDxjFnC41jzl7vvutryfXqwZw5UK1a6IikvK69FgYMgOOOg2efDR2N5BmNYxbJlEGD/PV55ykp57rzz4cqVeCFF2D27NDRiBRLiVki0b17d7p37x46jE337bf+IF6lij+oS27bdls48URYv75wvvMcFptyJhtQYpZIzJs3j3nz5oUOY9Pdc48/iHfr5g/qkvuSncCGDfNTduaw2JQz2YASs0hJli71B2/QKlJx0r49dOzo588eMSJ0NCJ/osQsUpJHHoElS3zHr7ZtQ0cjmZSsNd91l28REckiSswixVm3zh+0QbXlODruOGjUCGbMgDfeCB2NyAaUmCUSHTt2pGPHjqHDKL9XXoGvv/ZTbx51VOhoJNOqVClcsjOHV53K+XImxdI45myhcczZ5YAD4K234I474NJLQ0cjUViwAJo08XNnz5gBLVuGjkhiLivGMZtZZzObaWazzOzKYu4/1cw+TVwmmdkeUcYjkpZp03xSrlULzjordDQSlbp14dRT/fbdd4eNRSRFZInZzAqAe4HDgJ2Bk81s5yK7fQvs55zbHRgAPBBVPFKxjj/+eI4//vjQYZRPsmmzZ0/YYouwsUi0kp3Ahg/3vbRzTE6XMylRlDXmDsAs59w3zrnVwFPA0ak7OOcmOecWJf78AGgcYTxSgRYuXMjChQtDh1F2P/0ETzzh58O+6KLQ0UjUdtsN9t/fj2d++OHQ0ZRZzpYz2agoE3MjYG7K3/MSt5WkJzAuwnhESjdkCKxe7Tt87bhj6GikIiRrzXff7XvjiwQWZWIubgmeYnuamdn++MR8RQn39zKzKWY25ZdffslgiCIpfv/dJ2bQEKl8csQRsMMOfu7sF18MHY1IpIl5HtAk5e/GwPyiO5nZ7sAw4GjnXLFtMs65B5xz7Zxz7erVqxdJsCIbrLm8336ho5GKUlBQ2PP+9tvDxiICVI7wuScDzc1se+B7oBtwSuoOZrYd8BxwmnPuywhjkQp24IEHhg6hbJwrXEXq0ku15nK+OfNMvyTk++/DpEmw996hI0pLzpUzSUuk45jNrAswCCgAHnbO3Whm5wE45+43s2HA8cB3iYesLW2Ml8YxSyTGj4eDD4ZttvFrLm+2WeiIpKL17w833gjHHgvPPRc6GomhdMcxR1ljxjk3Fhhb5Lb7U7bPBs6OMgaRtNx2m7+++GIl5XzVu7f/HrzwAnz1FTRvHjoiyVOaklMicdhhh3HYYYeFDiM9n3wCr78ONWtqzeV81qABdO/uT2vceWfoaNKSU+VM0qbELJFYuXIlK1euDB1GepIdfnr2hK23DhuLhNW3r79+5BHIgREgOVXOJG1KzJLf5s71vbFTe+ZK/tplF+jSZcOhcyIVTIlZ8tugQbB2LZxwgl9JSqRfP399zz2g2qgEoMQs+WvxYnggMT37ZZcFDUWySKdO0Latb8oeOTJ0NJKHIu2VLfnriCOOCB1C6YYOhWXL/BKPbduGjkayhZmvNZ9yCgwcCGefDZWysw6TE+VMykzrMWcLjWOuWKtWwfbbww8/wLhx0Llz6Igkm6xZAzvt5Me0v/iinztdZBNlxXrMIlnriSd8Ut5tNzj00NDRSLapUqVwvvRbbw0aiuQfJWaJRKdOneiUbAXINuvXFw6R6tdP029K8c4+G7baCt57D959N3Q0xcrqciblpsQs+eeVV+CLL6BRI+jWLXQ0kq1q1y5ck/v//i9sLJJXlJglvzjn50MGX1uuWjVsPJLdLr7Yzwj36qvw8ceho5E8ocQs+eXNN+G//4W6deGcc0JHI9muTh3o1ctv33RT2FgkbygxS35J1pYvvdTXhERK849/+M5go0fDzJmho5E8oHHMEokTTzwxdAh/9v778NZbsMUWcOGFoaORXNGoEfToAQ8+CLfcAg8/HDqiP2RlOZNNpnHM2ULjmKN3xBG+49c118ANN4SORnLJrFnQsqWfaOTrr2G77UJHJDlI45glqBUrVrBixYrQYRSaNs0n5Ro1CseniqRrp53gpJP8vOoDB4aO5g9ZV84kI5SYJRJdunShS5cuocMolBzucu65vuOXSFldeaW/fvDBrFkSMuvKmWSEErPE34wZvuNO1aq+I49Ieey+uz8dsnIl3HFH6GgkxpSYJf5uvtmPX+7Rw3fkESmv/v399d13w4IFYWOR2FJilnibNQseewwKCuCKK0JHI7luzz3hsMNg+fLCaV1FMkyJWeJtwABYtw5OPx122CF0NBIH//qXv77nnqw51yzxonHMEokePXqEDsFPBvHYY1C5Mvzzn6Gjkbjo0AEOP9z38r/ttqCrT2VFOZOM0zjmbKFxzJl36ql+ecdevWDo0NDRSJxMnuwTdI0a8O23sM02oSOSHKBxzBLUggULWBCyc8wXX8CTT/qpFK+5JlwcEk/t2/se2itWBK0xBy9nEgklZolE165d6dq1a7gA/vUv3xP7nHM0S5NEI3mu+b774KefgoQQvJxJJJSYJX4+/RSeeQY22wyuvjp0NBJXf/0rHHWUH9d8yy2ho5EYUWKW+EnWZM47T+OWJVrXXeev77sP5s0LG4vEhhKzxMtHH8Hzz0P16oVTKIpEpW1bOOEEWLWq8AehyCZSYpZ4SSbjCy+EBg3CxiL54YYb/AQ2jzwC06eHjkZiQOOYJRLnn39+xb/oG2/A+PF+veWrrqr415f81KIFnH22H5J3zTXw3HMV9tJByplETuOYs4XGMW+a9euhXTv4+GM/N7am35SKNH++Xxpy5Up4/33Ya6/QEUkW0jhmCWru3LnMnTu34l7wqad8Um7UCC6+uOJeVwSgYcPCdb6vvNIP1asAFV7OpEIoMUskTjvtNE477bSKebFVqwonEbn+et/xS6SiXX45bLUVvP02vPZahbxkhZYzqTBKzJL7hg6F2bNh5539YhUiIWy5ZeG4+Suv9KdXRMpBiVly22+/+RWkAG66yS9YIRLKhRdC48bwyScwYkToaCRHKTFLbrvhBr9g/b77wpFHho5G8l316v4HIviRAcuWhY1HcpISs+Sur76CwYPBDAYN8tcioZ1yil956scfC5O0SBmo3U8i8Y9//KMiXgTWrIEzz/TzFotkg0qV/A/Gjh1h4EC/kEqzZpG8VIWUM6lwGsecLTSOuWxefx0OPRRq1fI1Z83yJdkmuR74CSfAqFGho5EsoHHMEtTMmTOZOXNmNE++di1ceqnf7t9fSVmy0803+3POzzwD774byUtEWs4kGCVmicS5557LueeeG82T338/fPEF7Lhj4aQOItmmSRM/thn893Tduoy/RKTlTIJRYpbc8ssvcO21fnvgQL/mski2uvxyP3xq6lT/g1IkDUrMklsuvxwWLYKDD/aL1Itksxo1fEcw8JOP/Phj2HgkJygxS+545x0YPtzXku+9V8OjJDcceyx06QJLlviRBCKliDQxm1lnM5tpZrPM7E+r1ptZKzN738xWmVm/KGORHLd6NSSXuLvqKmjePGw8Iukyg7vvhmrVfC/tCRNCRyRZLrJxzGZWANwLHAzMAyab2UvOuS9SdvsVuBg4Jqo4JIz+/ftn9gnvuMN3+GreXEs6Su7ZYQc/gqB/f7jgAvj004z0j8h4OZOsEGWNuQMwyzn3jXNuNfAUcHTqDs65n51zk4E1EcYhARx00EEcdNBBmXmy2bP9qlEA993nax4iuaZfP2jVCr78Em67LSNPmdFyJlkjysTcCEhdKHRe4jbJA9OmTWPatGmb/kTOQa9efgH6k08GHYQkV222mf9hCX6O9xkzNvkpM1bOJKtEmZiL65lTrmnGzKyXmU0xsym//PLLJoYlFaFPnz70ycQY44cfhjfegK23hjvv3PTnEwlp//39FLKrVvnrTRzbnLFyJlklysQ8D2iS8ndjYH55nsg594Bzrp1zrl29evUyEpzkgHnzoG9fv3333VC/fth4RDLhjjugYUP44IPCoVQiKaJMzJOB5ma2vZlVBboBL0X4ehInzsG55/ohJkcd5ZuxReJgyy1h6FC/fc01fq53kRSRJWbn3FqgN/AaMB0Y5Zz73MzOM7PzAMysgZnNA/oC/c1snpltHlVMkkMeewzGjvUHsfvv15hliZcjjoDu3eH336FnT1i/PnREkkUiXfbROTcWGFvktvtTtn/EN3GLFJo7Fy6+2G8PGgTbbhs0HJFIDB7s+0+8+67fTi7MInlP6zFLJP7v//6vfA9ctw5OPx0WL/a1itNPz2hcIllj663hgQfg6KPhyivhwANh993L9BTlLmeS1ZSYJRJ77713+R44cKBfk3qbbeChh9SELfF21FF+OOADD8App8DkyX6pyDSVu5xJVtNc2RKJSZMmMWnSpLI96KOP/MxI4OfE3mabjMclknXuuANatIDPP/c15zIoVzmTrKcas0Ti6quvBmDixInpPWD5cjj1VFizBi66CA47LLrgRLJJzZp+Du299oK77oLOndP+/pe5nElOUI1ZwnMOeveGmTNhl13glltCRyRSsf76VxgwwG+fcQZ8/33YeCQoJWYJ7+GHfdN19erw5JNlOscmEhuXXeannP3lFzjpJN96JHlJiVnC+vhjuPBCv33//bDbbmHjEQmloMA3aTdqBO+9V+bzzRIfSswSzuLF0LWrnze4Vy8NjRKpVw9GjYLKlX2nsGefDR2RBKDOXxKJQYMGbXyH9ev9ubRvvoG2bTVnsEjS3nv7ZSEvvdQvdLHLLn65yGKUWs4kJykxSyTatGmz8R3694eXXvJTbo4erTWWRVJdcglMmgTPPANHHgn//a+fkKSIUsuZ5CQ1ZUskxo8fz/jx44u/c+RIuOkmf05t9GjYfvuKDU4k25nBI49AmzYwa5Y/5VNMZ7CNljPJWaoxSyRuuOEGAA466KAN75g0Cc4+22/ffbefhlBE/qxmTd+q1L49vPWWH98/ZMgGs+GVWM4kp6nGLBVn9mw49lhYvdqPWz7//NARiWS3Jk3gxRdhs838UpF33RU6IqkASsxSMX7+GQ45xF8ffDDceWfoiERyw557+mZt8B3Cnn46bDwSOSVmid7SpdCli18Qvk0bf165ss6iiKTt5JN9vwzn4LTTQOeVY02JWaK1ahUcd5xfoGKHHWDcONh889BRieSeK66APn18J7Bjj4UpU0JHJBFRtUUiMXToUFi7tvDXff368Prr0KBB6NBEcpOZXxb1l1/g8cehSxeGPvII7LRT6Mgkw5SYJRItd9zRrxb1zDO+hjxuHOy4Y+iwRHJbpUr+fPOvv8K4cbQ86yzfY1tiRU3Zknlr1/Ly/vvz8qhRPim/9hr85S+hoxKJhypV/FSdBx/Myz//zMt77w1ffBE6KskgJWbJrDVroHt3Bv7nPww0g1df9evMikjmVK8OL77IwK22YuBvv8E++8Dnn4eOSjJEiVkyZ8UK3yklOZyjWjU/1ENEMq96ddh1Vz/GefFi6NjRT90pOU+JWTLj11/9WrKvvOKb2mrU8FNuikh0KlXyHSpr1/bDEvfbz7dSSU5TYpZNN28e/O1v8P77vpa8885alEKkolSq5Fef2nZbPzzxiCPgscdCRyWbQIlZNs2HH0KHDr7zSY0afom6Jk1CRyWSXypV8sunNmsG69b5YYrXXuuXV5Wco8Qs5ffEE/D3v8MPP/imtDZt/K92YORf/sLIpk3DxicScyNHjmRk167+DzN/zjm5dvOAAX5yn+XLwwUo5aLELGW3bh1cdZUfp7xqFdSp41fASVkvtkn16jSpWjVgkCLx16RJE5psueWGN+60ky+PBQV+AYwOHfwCMpIzlJilbObP9528br7Z/0Jv0gTatfPN2Cmenj+fpxctChSkSH54+umnefrTT/98R/36vt9HtWr+NNPuu8MLL1R4fFI+SsySvldfhT32gIkToWpVaN7cF/gqVf6065DZsxmyYEHFxyiSR4YMGcKQDz8s/s5atfyppjp1fI/tY4/1y62uWlWxQUqZKTFL6VasgL594bDDYMECP5tXmzbQosUGi7aLSJapWtVP8NOypS+r997rW7g++yx0ZLIRSsyyce+955PwnXf6gr3ttn7SkG22CR2ZiKTDzLdu7b23n4zkf//zPbgHDPAz9UnWUWKW4i1d6hdl/9vf/DrK1atD69a+QG+2WejoRKSsttoK9t8fGjb0K79de63vJDZ1aujIpAglZtmQc34YVMuWMGiQv61+fd8ctsMOaroWyWWVK/sf1x06+GbuTz7xTdvnnedn75OsoMQshaZO9VP6nXqqH5tcs6afxatdO79dBqPbtWP09ttHFKiIAIwePZrRJ59c9gdusw0ccIAfVeEcDB3qh1kNGaLm7SygxCwwYwaccAL89a/w7ru+l3WzZv6c1Pbbl6uWXLdqVepW1nLfIlGqW7cudcv4o/kPlSv7URZ//7vv0LloEVxwgZ+g5MknNWtYQErM+ezrr6FnTz+N5ujRflq/+vV9M1dy1ZpyGj53LsMXLsxgsCJS1PDhwxm+qeeIN9/c9yVp08aX+W++gVNO8X+/+KISdACq0uSjKVPg1lv9Yuvr1/sacd26vpZcv35GziMPnzsXVq+mxyY/k4iUZPjw4TBnDj02dZSEGTRu7DuGzZ7tO3x+9hkcc4wfFnnFFf4Ulzp+VgjVmPPFmjU+ER9wgO+J+cwz/vY6dfwv4z339MvHqXOXSP6qVMl38jzoIJ+Qq1aFL7/0LWtNm8INN/jZ/yRSSsxx9+23cPXVvpNH167w1lt+Dt369f055b32gkaNlJBFpFBBgU/MBx3kT2tVrw4//QT//Cdst52fRezVV/28+ZJxasqOox9/9DXip56CSZMKb69WrbDJuujE9yIiRVWq5I8XTZv6xDxrFixe7OfdfuEF38p20knQrZtvddMP/IxQYo6Lr7+GsWN9YZk4sbDDRqVKfmKBevV84SpmXmsRkY0y80m4QQNYudIfb374wVcCBg/2l6ZN4cQT4fDD/YgOHWvKzZxzoWMok3bt2rkpU6aEDiPzOnXy1xMnprf/ihW+Njx2rL/MnFl4nxlssYVPyNtt59dKrmAr1q2D+fOpMXWq/3EgIhm3YsUKuPVWarz6qv/xXZGcg4ULfWexhQs3HP+8xRZwyCHQpYufbUxrswNgZh8559qVtp9qzLli0SKfiN95x481njJlw4JQUOALw+ab+96VW2wRtFmpRkGBErJIxGrUqOE7aIWQHM1Rt65P0j//DHPn+qbu337zp9OSnUybNvXjpffbD/bd18/dreNDiZSYs9GCBfDxx/DRR342ro8+8mMLi6pRwyfirbbyHbiqVav4WEtw3+zZsHgxF4QORCTG7rvvPvjgg/DlzMx3KK1f3/+9bBnMmeNr0suWwXffwciR/gL+uPWXv/hZBf/6Vz9N6I47+klPRIk5mFWr/K/LmTP9zFszZ/rm6Xr1fGIuqlIln4hr1/brrG67rd/O0s4Wo+bPh9Wrwx8wRGJs1KhRMGcOF2Tbam+1avnpfMHXpn/91Z+PXrQIli+HJUvg7bf9JSm5xnvr1v6xrVv7v5s1g623ztpjXRQiTcxm1hkYDBQAw5xzNxe53xL3dwFWAD2cc7m91Ilz/kv388/+8tNPftzfnDn+V+N33/ntH3/0+xYnmYRr1vS14K239s1FGtwvIrnGzM+XUKdO4W0rV/rj48KFPlGvXAmrV8Pnn/tLUTVr+gSdvGy3na+c1K/vO6TVr++PkQUFFfSmohVZYjazAuBe4GBgHjDZzF5yzn2RstthQPPEZU9gSOI6nDVr/Jdk5Upfg12+3J8vKemyeLG/JJPwzz/72nA6qlb1ibd6df8cBQV+7tottlCTjojEV/Xq/rxzaqewNWt8jXrRIr/s7O+/+2S9apU/DpeUtJMqVfItjvXr+x8BW27pL1ttVbid/HuLLXyyT1aAktebbZYVNfMoj/4dgFnOuW8AzOwp4GggNTEfDYxwvmv4B2a2pZlt65z7ocRn/flnGDjQ/xPXrPHriia30/l7zRr/D1+xYsMEnNxeu3bT33mlSn6oQPJSubK/3mwz/8+vVcs3Q1etWvglSI43Tv1VKSKSL6pU8ateFdcsv2qVrwgtXVpYw163bsPj+9q1vnL000/ljyHZWpmasKtV88fqqlV9jMntjV0qV/aXgoINr9MUZWJuBMxN+Xsef64NF7dPI6DkxDx3LvTrl6EQS1Cpkk+YqdeVKvkPt+h1atKtXt3/EzfbrPQeh2vXbvgjIDnuePny6N5XRVq/vuSmehHJHOf8caNGjdCRRKtmzY0vP7tuna90rVzpE/nq1Rsm7fXr/WXdOn9xrvA4ldxev953Vlu2LJK3sA2kVfOKMjEX1x5Q9Eidzj6YWS+gF0ABrG0NK8lQg0Pixf741zhwzidJV5HTzf0G1baA33nrrQp7zaitg3VWULACKKY3W86ri95XLont+9oa1raAOg5iswzUH8fDimdWzGVTnzSZZ5ZAk3T2jzIxz2PDIBoDRWc/T2cfnHMPAA8AmNmU6WkM0M41Zjbl55i+r3QG1Ocava/cEuf3tTCm7yuux8N09otyhPdkoLmZbW9mVYFuwEtF9nkJON28vYDfNnp+WUREJOYiqzE759aaWW/gNfxwqYedc5+b2XmJ++8HxuKHSs3CD5c6M6p4REREckGkY3Kcc2PxyTf1tvtTth1wYRmf9oEMhJaN9L5yi95XbtH7yi15/b5ybhELERGRONMs4iIiIlkkpxOzmfUzM2dmdUPHkglmNsDMPjWzaWb2upk1DB1TJpjZbWY2I/HenjezLUPHlAlmdoKZfW5m680s53uQmllnM5tpZrPM7MrQ8WSCmT1sZj+b2f9Cx5JJZtbEzN4ys+mJ7+AloWPKBDOrZmYfmtkniff179AxZYqZFZjZx2Y2prR9czYxm1kT/HSfc0LHkkG3Oed2d861AcYA1waOJ1PeAHZ1zu0OfAlcFTieTPkfcBzwTuhANlXKFLqHATsDJ5vZzmGjyojhQOfQQURgLfAP51xrYC/gwpj8v1YBBzjn9gDaAJ0TI3bi4BJgejo75mxiBu4ELqeYCUlylXNuScqfNYnJe3POve6cS05z9gF+vHrOc85Nd87NDB1Hhvwxha5zbjWQnEI3pznn3gF+DR1Hpjnnfkgu+OOcW4o/4DcKG9Wmc15y2q0qiUvOHwfNrDFwODAsnf1zMjGb2VHA9865T0LHkmlmdqOZzQVOJT415lRnAeNCByF/UtL0uJLlzKwZ8Bfgv4FDyYhEk+804GfgDedcHN7XIHxFMq3Z2bJ2CSMzGw80KOaua4CrgUMqNqLM2Nj7cs696Jy7BrjGzK4CegPXVWiA5VTa+0rscw2+Ce7xioxtU6TzvmIirelxJbuYWS3gWaBPkRa3nOWcWwe0SfRFed7MdnXO5WwfATM7AvjZOfeRmXVK5zFZm5idcwcVd7uZ7QZsD3zil3OmMTDVzDo4536swBDLpaT3VYwngFfIkcRc2vsyszOAI4ADXQ6N0SvD/yvXpTU9rmQPM6uCT8qPO+eeCx1PpjnnFpvZRHwfgZxNzMA+wFFm1gWoBmxuZo8557qX9ICca8p2zn3mnNvGOdfMOdcMf0BpmwtJuTRm1jzlz6OAGaFiySQz6wxcARzlnFsROh4pVjpT6EqWMF8reQiY7py7I3Q8mWJm9ZKjNsysOnAQOX4cdM5d5ZxrnMhX3YA3N5aUIQcTc8zdbGb/M7NP8U31sRgCAdwD1AbeSAwFu7+0B+QCMzvWzOYBHYFXzOy10DGVV6JzXnIK3enAKOfcRlalzw1m9iTwPtDSzOaZWc/QMWXIPsBpwAGJMjUtUSPLddsCbyWOgZPx55hLHV4UN5r5S0REJIuoxiwiIpJFlJhFRESyiBKziIhIFlFiFhERySJKzCIiIllEiVkkRsxsopkdWuS2PmZ230b2z/mVsUTiRIlZJF6exE9ikKpb4nYRyQFKzCLxMho4wsw2gz8WOGgInGJmUza2xq2ZLUvZ7mpmwxPb9czsWTObnLjsk7h9v5TJLT42s9oRvzeRvJC1c2WLSNk55xaa2Yf4+YVfxNeWnwZucs79mlh3eYKZ7e6c+zTNpx0M3Omc+4+ZbYefHaw10A+40Dn3XmIxhd8z/oZE8pBqzCLxk9qcnWzGPtHMpgIfA7sAO5fh+Q4C7kksxfcSfhL+2sB7wB1mdjGwZcqa2yKyCZSYReLnBeBAM2sLVAcW4Wu3BzrndsevWlatmMelzs+ben8loKNzrk3i0sg5t9Q5dzNwduI1PjCzVhG8F5G8o8QsEjPOuWXAROBhfG15c2A58JuZ1QcOK+GhP5lZazOrBBybcvvr+AUuADCzNonrHROrvd0CTAGUmEUyQIlZJJ6eBPYAnnLOfYJvwv4cn6zfK+ExVwJjgDeBH1JuvxhoZ2afmtkXwHmJ2/skVkP7BFgJjMv82xDJP1pdSkREJIuoxiwiIpJFlJhFRESyiBKziIhIFlFiFhERySJKzCIiIllEiVlERCSLKDGLiIhkESVmERGRLPL/rDi4ZfFdn0kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xval = np.linspace(-4.0,4.0,1000)\n",
    "tpdf = stats.t.pdf(xval,loc=0,df=dof,scale=1.0)\n",
    "plt.plot(xval,tpdf,color='red',lw=2); plt.xlabel('Values'); plt.ylabel('Density'); plt.title('Hypothesis Test, Pooled (Equal Variance) Difference in Means')\n",
    "plt.xlim([-4.0,4.0]); plt.ylim([0,0.6])\n",
    "\n",
    "plt.annotate('$-t_{critical} = $' + str(t_critical[0]),[t_critical[0]-0.18,0.4],rotation=90.0)\n",
    "plt.annotate('$+t_{critical} = $' + str(t_critical[1]),[t_critical[1]-0.18,0.4],rotation=90.0)\n",
    "plt.annotate('$t_{pooled} = $' + str(np.round(t_pooled,3)),[t_pooled-0.18,0.4],rotation=90.0,color='red')\n",
    "    \n",
    "plt.fill_between(xval,tpdf,where= xval < t_critical[0], color='red',alpha=0.7)\n",
    "plt.fill_between(xval,tpdf,where= xval > t_critical[1], color='red',alpha=0.7)\n",
    "plt.vlines(t_critical,0,7,color='black',ls='--')\n",
    "plt.vlines(t_pooled,0,0.7,color='red')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.1, wspace=0.3, hspace=0.4); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Welch's t-test Difference in Means \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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The t statistic is -2.98 and the p-value is 0.0055\n"
     ]
    }
   ],
   "source": [
    "t_Welch,p_Welch = stats.ttest_ind(X1, X2, equal_var = False) # allowing for difference in variance\n",
    "print('The t statistic is ' + str(np.round(t_Welch,2)) + ' and the p-value is ' + str(np.round(p_Welch,5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we can see by $p$ that we will clearly reject the null hypothesis.  \n",
    "\n",
    "Now let's plot the result!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xval = np.linspace(-4.0,4.0,1000)\n",
    "tpdf = stats.t.pdf(xval,loc=0,df=welch_dof(X1,X2),scale=1.0)\n",
    "plt.plot(xval,tpdf,color='red',lw=2); plt.xlabel('Values'); plt.ylabel('Density'); plt.title('Hypothesis Test, Welch\\'s (Unequal Variance) Difference in Means')\n",
    "plt.xlim([-4.0,4.0]); plt.ylim([0,0.6])\n",
    "\n",
    "plt.annotate('$-t_{critical} = $' + str(t_critical[0]),[t_critical[0]-0.18,0.4],rotation=90.0)\n",
    "plt.annotate('$+t_{critical} = $' + str(t_critical[1]),[t_critical[1]-0.18,0.4],rotation=90.0)\n",
    "plt.annotate('$t_{Welch} = $' + str(np.round(t_Welch,3)),[t_Welch-0.18,0.4],rotation=90.0,color='red')\n",
    "    \n",
    "plt.fill_between(xval,tpdf,where= xval < t_critical[0], color='red',alpha=0.7)\n",
    "plt.fill_between(xval,tpdf,where= xval > t_critical[1], color='red',alpha=0.7)\n",
    "plt.vlines(t_critical,0,7,color='black',ls='--')\n",
    "plt.vlines(t_pooled,0,0.7,color='red')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.1, wspace=0.3, hspace=0.4); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### F-test Difference in Variances\n",
    "\n",
    "Let's now compare the variances with the F-test for difference in variances.  \n",
    "\n",
    "\\begin{equation}\n",
    "H_0: \\frac{\\sigma^{2}_{X_2}}{\\sigma^{2}_{X_1}} = 1.0\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "H_1: \\frac{\\sigma^{2}_{X_2}}{\\sigma^{2}_{X_1}} > 1.0\n",
    "\\end{equation}\n",
    "\n",
    "Note, by ordering the variances we eliminate the case of $\\sigma^{2}_{X_2} \\lt \\sigma^{2}_{X_1}$.\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value = 1 - stats.f.cdf(np.var(X2)/np.var(X1), dfn=len(X2)-1, dfd=len(X1)-1)\n",
    "p_value"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compact Code Solution - Hypothesis Test, Student's t-test for Difference in Means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "\n",
    "df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/PorositySample2Units.csv')\n",
    "X1 = df['X1'].values; X2 = df['X2'].values\n",
    "\n",
    "alpha = 0.05\n",
    "t_critical = stats.t.ppf(alpha/2,len(X1)+len(X2)-2)\n",
    "\n",
    "t_pooled, p_pooled = stats.ttest_ind(X1,X2) # assuminng equal variance\n",
    "if(t_pooled < t_critical or t_pooled > -1*t_critical):\n",
    "    print('Equal Variance: t-critical and t-statistic are ' + str(np.round(t_critical,2)) + ' ≤ ' + str(np.round(t_pooled,2)) + \n",
    "      ' ≤ ' + str(np.round(-1*t_critical,2)) + '; therefore, reject the null hypothesis')\n",
    "else:\n",
    "    print('Equal Variance: t-critical and t-statistic are ' + str(np.round(t_critical,2)) + ' ≤ ' + str(np.round(t_pooled,2)) + \n",
    "      ' ≤ ' + str(np.round(-1*t_critical,2)) + '; therefore, fail to reject the null hypothesis')\n",
    "    \n",
    "    alpha = 0.05\n",
    "t_critical = stats.t.ppf(alpha/2,len(X1)+len(X2)-2)\n",
    "\n",
    "t_welch, p_welch = stats.ttest_ind(X1,X2,equal_var = False) # assuming unequal variance, Welch's t-test\n",
    "if(t_welch < t_critical or welch > -1*t_critical):\n",
    "    print('Welch\\'s test: t-critical and t-statistic are ' + str(np.round(t_critical,2)) + ' ≤ ' + str(np.round(t_welch,2)) + \n",
    "      ' ≤ ' + str(np.round(-1*t_critical,2)) + '; therefore, reject the null hypothesis')\n",
    "else:\n",
    "    print('Welch\\'s test: t-critical and t-statistic are ' + str(np.round(t_critical,2)) + ' ≤ ' + str(np.round(t_welch,2)) + \n",
    "      ' ≤ ' + str(np.round(-1*t_critical,2)) + '; therefore, fail to reject the null hypothesis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compact Code Solution - Hypothesis Test, F-test for Difference in Variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f-statistic and f-critical are 2.67 > 2.17; therefore, reject the null hypothesis\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "\n",
    "df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/PorositySample2Units.csv')\n",
    "X1 = df['X1'].values; X2 = df['X2'].values\n",
    "\n",
    "alpha = 0.05\n",
    "\n",
    "var_X1 = np.var(X1,ddof=1); var_X2 = np.var(X2,ddof=1) # sample variance\n",
    "f_stat = np.max([var_X1,var_X2])/np.min([var_X1,var_X2])\n",
    "if var_X1 > var_X2:\n",
    "    f_critical = stats.f.ppf(1-alpha,len(X1)-1,len(X2)-1) \n",
    "else:\n",
    "    f_critical = stats.f.ppf(1-alpha,len(X2)-1,len(X1)-1) \n",
    "if f_stat > f_critical:\n",
    "    print('f-statistic and f-critical are ' + str(np.round(f_stat,2)) + ' > ' + str(np.round(f_critical,2)) \n",
    "          + '; therefore, reject the null hypothesis')\n",
    "else:\n",
    "    print('f-statistic and f-critical are ' + str(np.round(f_stat,2)) + ' ≤ ' + str(np.round(f_critical,2)) \n",
    "          + '; therefore, fail to reject the null hypothesis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "We are just scratching the surface for confidence intervals and hypothesis tests.  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 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",
    "*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",
    "On twitter I'm the @GeostatsGuy.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}