PythonDataBasics_Statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "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 Science Basics in Python Series\n",
    "\n",
    "## Chapter VI: Basic Statistical Analysis in Python \n",
    "\n",
    "### Michael Pyrcz, Associate Professor, The University of Texas at Austin \n",
    "\n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Basic Univariate Statistics\n",
    "\n",
    "Here's a demonstration of the calculation and visualization for basic statistical analysis in Python. \n",
    "\n",
    "We will use the following Python packages:\n",
    "\n",
    "* [statistics](https://docs.python.org/3/library/statistics.html)\n",
    "\n",
    "* [SciPy](https://www.scipy.org/)\n",
    "\n",
    "* [MatPlotLib](https://matplotlib.org/)\n",
    "\n",
    "We will cover a variety of common, basic statistical analyses and displays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Basic Univariate Statistics\n",
    "\n",
    "This tutorial includes the methods and operations that would commonly be required for Engineers and Geoscientists working with Regularly Gridded Data Structures for the purpose of:\n",
    "\n",
    "1. Data Checking and Cleaning\n",
    "2. Data Mining / Inferential Data Analysis\n",
    "3. Predictive Modeling\n",
    "\n",
    "for Data Analytics, Geostatistics and Machine Learning."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General Definitions\n",
    "\n",
    "**Statistics**\n",
    "\n",
    "collecting, organizing, and interpreting data, as well as drawing conclusions and making decisions.\n",
    "\n",
    "**Geostatistics**\n",
    "\n",
    "a branch of applied statistics that integrates: \n",
    "\n",
    "1. the spatial (geological) context,\n",
    "\n",
    "2. the spatial relationships, \n",
    "\n",
    "3. volumetric support / scale, and \n",
    "\n",
    "4.  uncertainty.  \n",
    "\n",
    "**Data Analytics**\n",
    "\n",
    "use of statistics [with visualization] to support decision making.\n",
    "\n",
    "**Big Data Analytics** \n",
    "\n",
    "process of examining large and varied data sets (big data) to discover patterns and make decisions."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General Definitions\n",
    "\n",
    "**Variable** or **Feature** \n",
    "\n",
    "* any property measured / observed in a study (e.g., porosity, permeability, mineral concentrations, saturations, contaminant concentration)\n",
    "* the measure often requires significant analysis, interpretation and uncertainty, 'data softness'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General Definitions\n",
    "\n",
    "**Population** \n",
    "\n",
    "exhaustive, finite list of property of interest over area of interest. Generally the entire population is not accessible.\n",
    "\n",
    "* exhaustive set of porosity at each location within a reservoir\n",
    "\n",
    "**Sample** \n",
    "\n",
    "set of values, locations that have been measured\n",
    "\n",
    "* porosity data from well-logs within a reservoir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General Definitions\n",
    "\n",
    "**Parameters**\n",
    "\n",
    "summary measure of a population\n",
    "\n",
    "* population mean, population standard deviation, we rarely have access to this\n",
    "\n",
    "* model parameters is different and we will cover later.\n",
    "\n",
    "**Statistics**\n",
    "\n",
    "summary measure of a sample\n",
    "\n",
    "* sample mean, sample standard deviation, we use statistics as estimates of the parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Covered Parameters / Statistics\n",
    "\n",
    "We cover the following parameters and statistics.\n",
    "\n",
    "| Central Tendency | Dispersion | Outliers | Distributions Shape |\n",
    "| :--------------: | :--------: | :------: | :-----------------: |\n",
    "| Arithmetic Average / Mean | Variance | Tukey Outlier Test| Skew            |\n",
    "| Median            | Standard Deviation |                 | Excess Kurtosis |\n",
    "| Mode              | Range |                              | Person's Mode Skewness |\n",
    "| Geometric Mean    | Percentile |                         | Quartile Skew Coefficient |\n",
    "| Harmonic Mean     | Interquartile Range |                |                           |\n",
    "| Power Law Average |                     |                |                           |\n",
    "\n",
    "\n",
    "I have a lecture on these univariate statistics available on [YouTube](https://www.youtube.com/watch?v=wAcbA2cIqec&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=11&t=0s).   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Nonparmetric Cumulative Distribution Functions (CDFs)\n",
    "\n",
    "**nonparametric CDF** \n",
    "* plotting nonparametric distributions \n",
    "\n",
    "**fitting CDFs**\n",
    "* fitting a parametric distribution and plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup to run this demonstration:\n",
    "\n",
    "1. **Install Anaconda 3** on your machine from https://www.anaconda.com/download/. \n",
    "2. **Open Jupyter Notebook**, look for the Jupyter app on your system after installing Anaconda 3.\n",
    "3. **Load this Workflow** found here [PythonDataBasics_Statistics.ipynb](https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/PythonDataBasics_PedictiveMachineLearning.ipynb). \n",
    "4. **Load the data**, this workflow retreives the data from my GitHub [GeoDataSets Repository](https://github.com/GeostatsGuy/GeoDataSets). If you want to work locally, you will need to first download the data file to your working directory. The data file is found here, [2D_MV_200wells.csv](https://github.com/GeostatsGuy/GeoDataSets/blob/master/2D_MV_200wells.csv). Code is provided below to set the working directory and to load the data locally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load the required libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                              # ndarrys for gridded data\n",
    "import pandas as pd                             # DataFrames for tabular data\n",
    "import os                                       # set working directory, run executables\n",
    "import matplotlib.pyplot as plt                 # plotting\n",
    "import scipy                                    # statistics\n",
    "import statistics as stats                      # statistics like the mode\n",
    "from scipy.stats import norm                    # fitting a Gaussian distribution"
   ]
  },
  {
   "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). Set this to your working directory, with the above mentioned data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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 multivariate, spatial dataset [2D_MV_200wells.csv](https://github.com/GeostatsGuy/GeoDataSets). These are the features:  \n",
    "\n",
    "| Feature  | Units | Descriptions |\n",
    "| :------: | :--------------: | :--------- |\n",
    "| X, Y  | $meters$     | Location |\n",
    "| porosity | $fraction$   | rock void fraction |\n",
    "| permeability | $mDarcy$ | capability of a porous rock to permit the flow of fluids through its pore spaces |\n",
    "| acoustic impedance | $\\frac{kg}{m^2s} 10^6$ | rock bulk density times rock acoustic velocity |\n",
    "\n",
    "* load use the Pandas 'read_csv' function and rename the features for readable code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>facies</th>\n",
       "      <th>porosity</th>\n",
       "      <th>perm</th>\n",
       "      <th>ai</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>565</td>\n",
       "      <td>1485</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1184</td>\n",
       "      <td>6.170</td>\n",
       "      <td>2.009</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2585</td>\n",
       "      <td>1185</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1566</td>\n",
       "      <td>6.275</td>\n",
       "      <td>2.864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2065</td>\n",
       "      <td>2865</td>\n",
       "      <td>2</td>\n",
       "      <td>0.1920</td>\n",
       "      <td>92.297</td>\n",
       "      <td>3.524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3575</td>\n",
       "      <td>2655</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1621</td>\n",
       "      <td>9.048</td>\n",
       "      <td>2.157</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1835</td>\n",
       "      <td>35</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1766</td>\n",
       "      <td>7.123</td>\n",
       "      <td>3.979</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      X     Y  facies  porosity    perm     ai\n",
       "0   565  1485       1    0.1184   6.170  2.009\n",
       "1  2585  1185       1    0.1566   6.275  2.864\n",
       "2  2065  2865       2    0.1920  92.297  3.524\n",
       "3  3575  2655       1    0.1621   9.048  2.157\n",
       "4  1835    35       1    0.1766   7.123  3.979"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"2D_MV_200wells.csv\")          # read a .csv file in as a DataFrame\n",
    "df = df.rename(columns={'facies_threshold_0.3': 'facies','permeability':'perm','acoustic_impedance':'ai'}) # rename columns of the \n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract a Feature \n",
    "\n",
    "Let's extract one of the features, porosity, into a 1D ndarray and do our statistics on porosity.\n",
    "\n",
    "* then we can use NumPy's statistics methods\n",
    "\n",
    "* note this is a **shallow copy**, any changes to the array will change the feature in the DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are working with Porosity (fraction) from column porosity .\n"
     ]
    }
   ],
   "source": [
    "feature = 'Porosity'; fcol = 'porosity'; funits = '(fraction)'; f2units = '(fraction^2)'; fmin = 0.0; fmax = 0.25\n",
    "X = df[fcol].values\n",
    "print('We are working with ' + feature + ' ' + funits + ' from column ' + fcol +  ' .')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize the Feature Histogram\n",
    "\n",
    "To improve our understanding of the feature, let's visualize the featre distribuiton as a histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAW6UlEQVR4nO3dfbRddX3n8feHAAKKBppoUYiBlupgrTxcmazaVkWcolbFGceRcVzUUtOZ6rI+dOrjTGnXuJbOqtJ5aEujKOjgAz4gDNXWqETrTAO9wUiCwSIPCpIlUckAKSUFvvPH3tdeknuTk3vuPufeu9+vtc66e//O73f293d38j37/vbev52qQpLUHweNOwBJ0miZ+CWpZ0z8ktQzJn5J6hkTvyT1zMHjDmAQK1asqNWrV487DElaVDZt2vTDqlq5Z/miSPyrV69mcnJy3GFI0qKS5LszlTvUI0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9syju3JUWqy0bN7J75845tz90+XKevmbNPEYkmfilTu3euZPTVu41VcrANu3YMY/RSI3Oh3qSLEvyjSRXtevHJ7kmyU1JPpnk0K5jkCT9k1GM8f8OsG3a+nuBC6rqROBu4LwRxCBJanWa+JMcC7wI+GC7HuAM4NNtlUuAs7uMQZL0SF0f8f8x8HvAw+36TwE7q+rBdv0O4EkzNUyyNslkkskdjnNK0rzpLPEn+TXgrqraNL14hqo1U/uqWldVE1U1sXKIk2OSpEfq8qqeZwEvSfJC4DDgsTR/ASxPcnB71H8scGeHMUiS9tDZEX9Vvb2qjq2q1cArga9U1auAq4GXt9XOBa7oKgZJ0t7GcefuW4E3J/kOzZj/RWOIQZJ6ayQ3cFXVBmBDu3wLcPootitJ2ptz9UhSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9Y+KXpJ7p8mHrhyW5Nsk3k9yQ5A/a8ouT3Jpkc/s6uasYJEl76/IJXA8AZ1TVfUkOAb6e5Avte/+xqj7d4bYlSbPoLPFXVQH3tauHtK/qanuSpMF0OsafZFmSzcBdwPqquqZ9691Jrk9yQZJHzdJ2bZLJJJM7duzoMkxJ6pVOE39VPVRVJwPHAqcn+Xng7cBTgWcCRwNvnaXtuqqaqKqJlStXdhmmJPXKSK7qqaqdwAbgrKraXo0HgA8Dp48iBklSo8urelYmWd4uHw6cCdyY5Ji2LMDZwNauYpAk7a3Lq3qOAS5JsozmC+ayqroqyVeSrAQCbAb+fYcxSJL20OVVPdcDp8xQfkZX25Qk7Z937kpSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9Y+KXpJ7p8tGLhyW5Nsk3k9yQ5A/a8uOTXJPkpiSfTHJoVzFIkvbW5RH/A8AZVfUM4GTgrCRrgPcCF1TVicDdwHkdxiBJ2kNnib8a97Wrh7SvAs4APt2WX0LzwHVJ0oh0+bB12getbwJ+FvgT4GZgZ1U92Fa5A3jSLG3XAmsBVq1a1WWY0oJ187ZtQ7U/dPlynr5mzTxFo6Wi08RfVQ8BJydZDlwO/LOZqs3Sdh2wDmBiYmLGOtJS9/CuXZy2cuWc22/asWMeo9FSMZKreqpqJ7ABWAMsTzL1hXMscOcoYpAkNbq8qmdle6RPksOBM4FtwNXAy9tq5wJXdBWDJGlvXQ71HANc0o7zHwRcVlVXJfkW8Ikk/wX4BnBRhzFIkvbQWeKvquuBU2YovwU4vavtSpo/WzZuZPfOnXNu78nlhanTk7uSFrfdO3d6cnkJcsoGSeoZE78k9YyJX5J6xjF+aQkb9s7fW7dtG2qMXwuTiV9awoa98/emycl5jEYLhUM9ktQzJn5J6hkTvyT1jIlfknrGxC9JPWPil6SeMfFLUs8MlPiT/HzXgUiSRmPQI/4Lk1yb5LenHq4iSVqcBkr8VfVLwKuA44DJJB9L8vxOI5MkdWLgMf6qugl4F/BW4NnAf09yY5J/OVP9JMcluTrJtiQ3JPmdtvz8JN9Psrl9vXA+OiJJGsxAc/Uk+QXgNcCLgPXAi6vquiRPBP4G+OwMzR4E3tLWOxLYlGR9+94FVfVHw4cvSTpQg07S9j+BDwDvqKr7pwqr6s4k75qpQVVtB7a3y/cm2QY8ach4JUlDGnSo54XAx6aSfpKDkhwBUFUf3V/jJKtpnr97TVv0+iTXJ/lQkqMOOGpJ0pwNmvi/BBw+bf2Itmy/kjwG+Azwxqq6B/gz4GeAk2n+InjfLO3WJplMMrnD53ZK0rwZNPEfVlX3Ta20y0fsr1GSQ2iS/qVV9dm27Q+q6qGqephm+Oj0mdpW1bqqmqiqiZU+CEKS5s2giX9XklOnVpKcBty/j/okCXARsK2q3j+t/Jhp1V4GbB08XEnSsAY9uftG4FNJ7mzXjwH+zX7aPAt4NbAlyea27B3AOUlOBgq4DfitA4pYkjSUgRJ/Vf1tkqcCTwEC3FhV/7ifNl9v6+7p8wccpSRp3hzIM3efCaxu25yShKr6SCdRSZI6M+gNXB+luRJnM/BQW1yAiV+SFplBj/gngJOqqroMRpLUvUGv6tkK/HSXgUiSRmPQI/4VwLeSXAs8MFVYVS/pJCpJUmcGTfzndxmEJGl0Br2c86tJngycWFVfaufpWdZtaJKkLgz66MXXAp8G/rwtehLwua6CkiR1Z9CTu6+juRP3HvjJQ1ke31VQkqTuDJr4H6iq3VMrSQ6muY5fkrTIDJr4v5rkHcDh7bN2PwX87+7CkiR1ZdDE/zZgB7CFZlK1z9M8f1eStMgMelXP1Nz5H+g2HElS1wadq+dWZhjTr6oT5j0iSVKnDmSunimHAf8aOHr+w5EkdW2gMf6q+tG01/er6o+BMzqOTZLUgUGHek6dtnoQzV8AR3YSkSSpU4MO9bxv2vKDNI9MfMW+GiQ5jma+/p8GHgbWVdV/S3I08Emah7rcBryiqu4+oKglSXM26FU9z53DZz8IvKWqrktyJLApyXrg14EvV9V7kryN5lLRt87h8yVJczDoUM+b9/V+Vb1/hrLtwPZ2+d4k22jm+Hkp8Jy22iXABkz8kjQyB3JVzzOBK9v1FwNfA24fpHGS1cApwDXAE9ovBapqe5IZ5/xJshZYC7Bq1aoBw5Qk7c+BPIjl1Kq6FyDJ+cCnquo399cwyWOAzwBvrKp7kgy0wapaB6wDmJiYcF4gSZong07ZsArYPW19N83J2X1KcghN0r+0qj7bFv8gyTHt+8cAdw0crSRpaIMe8X8UuDbJ5TR38L6M5oqdWaU5tL8I2LbHOYArgXOB97Q/rzjQoCVJczfoVT3vTvIF4JfbotdU1Tf20+xZwKuBLUk2t2XvoEn4lyU5D/gezV3AkqQRGfSIH+AI4J6q+nCSlUmOr6pbZ6tcVV8HZhvQf96BBClJmj+DPnrx92kuuXx7W3QI8L+6CkqS1J1BT+6+DHgJsAugqu7EKRskaVEaNPHvrqqinZo5yaO7C0mS1KVBE/9lSf4cWJ7ktcCX8KEskrQoDXpVzx+1z9q9B3gK8J+ran2nkUmSOrHfxJ9kGfBXVXUmYLKXpEVuv0M9VfUQ8PdJHjeCeCRJHRv0Ov5/oLkRaz3tlT0AVfWGTqKSJHVm0MT/F+1LkrTI7TPxJ1lVVd+rqktGFZAkqVv7G+P/3NRCks90HIskaQT2l/inz7VzQpeBSJJGY3+Jv2ZZliQtUvs7ufuMJPfQHPkf3i7TrldVPbbT6CRJ826fib+qlo0qEEnSaAw6V48kaYnoLPEn+VCSu5JsnVZ2fpLvJ9ncvl7Y1fYlSTPr8oj/YuCsGcovqKqT29fnO9y+JGkGnSX+qvoa8OOuPl+SNDfjGON/fZLr26Ggo2arlGRtkskkkzt27BhlfJK0pI068f8Z8DPAycB24H2zVayqdVU1UVUTK1euHFV8krTkjTTxV9UPquqhqnqY5glep49y+5KkESf+JMdMW30ZsHW2upKkbgw6LfMBS/Jx4DnAiiR3AL8PPCfJyTTTP9wG/FZX25ckzayzxF9V58xQfFFX25NmsmXjRnbv3Dnn9ocuX87T16yZx4j65eZt24Zq7++/G50lfmkh2L1zJ6cNcXHAJq8oG8rDu3b5+1+AnLJBknrGxC9JPWPil6SecYxf2odhT07eum3bUGPcUhdM/NI+DHty8qbJyXmMRpofDvVIUs+Y+CWpZ0z8ktQzJn5J6hkTvyT1jIlfknrGxC9JPWPil6SeMfFLUs+Y+CWpZzpL/Ek+lOSuJFunlR2dZH2Sm9qfR3W1fUnSzLo84r8YOGuPsrcBX66qE4Evt+uSpBHqLPFX1deAH+9R/FLgknb5EuDsrrYvSZrZqMf4n1BV2wHan4+frWKStUkmk0zu8PFrkjRvFuzJ3apaV1UTVTWx0vnMJWnejDrx/yDJMQDtz7tGvH1J6r1RJ/4rgXPb5XOBK0a8fUnqvS4v5/w48DfAU5LckeQ84D3A85PcBDy/XZckjVBnj16sqnNmeet5XW1TkrR/C/bkriSpGyZ+SeoZE78k9YyJX5J6xsQvST1j4peknjHxS1LPmPglqWdM/JLUM53duSsBbNm4kd07d865/fduv51Vxx035/a3btvGac7uKj2CiV+d2r1z51CJ96bJSU479dSh2kt6JId6JKlnTPyS1DMmfknqGcf4JS1YN2/bNlT7YS8OOHT5cp6+Zs1QMSxEJn5JC9bDu3aN9eKATTt2zLntQjaWxJ/kNuBe4CHgwaqaGEccktRH4zzif25V/XCM25ekXvLkriT1zLgSfwFfTLIpydoxxSBJvTSuoZ5nVdWdSR4PrE9yY1V9bXqF9gthLcCqVavGEaMkLUljOeKvqjvbn3cBlwOnz1BnXVVNVNXESudakaR5M/LEn+TRSY6cWgb+BbB11HFIUl+NY6jnCcDlSaa2/7Gq+ssxxCFJvTTyxF9VtwDPGPV2JUkNL+eUpJ4x8UtSz5j4JalnTPyS1DMmfknqGadl1j4N+7B0H3YuLTwmfu3TfDwsXdLC4lCPJPWMiV+SesbEL0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqGW/gWuCGvXP2e7ffzqrjjptze++8lZYeE/8CNx93zp526qlDtZe0tIxlqCfJWUm+neQ7Sd42jhgkqa/G8bD1ZcCfAC8ATgLOSXLSqOOQpL4axxH/6cB3quqWqtoNfAJ46RjikKReSlWNdoPJy4Gzquo32/VXA/+8ql6/R721wNp29SnAt+e4yRXAD+fYdrGyz/1gn/thmD4/uar2Okk4jpO7maFsr2+fqloHrBt6Y8lkVU0M+zmLiX3uB/vcD130eRxDPXcA068vPBa4cwxxSFIvjSPx/y1wYpLjkxwKvBK4cgxxSFIvjXyop6oeTPJ64K+AZcCHquqGDjc59HDRImSf+8E+98O893nkJ3clSePlXD2S1DMmfknqmUWd+Pc39UOSRyX5ZPv+NUlWT3vv7W35t5P86ijjHsZc+5xkdZL7k2xuXxeOOva5GKC/v5LkuiQPtveITH/v3CQ3ta9zRxf1cIbs80PT9vGiuWhigD6/Ocm3klyf5MtJnjztvaW6n/fV5+H2c1UtyhfNieGbgROAQ4FvAiftUee3gQvb5VcCn2yXT2rrPwo4vv2cZePuU8d9Xg1sHXcfOujvauAXgI8AL59WfjRwS/vzqHb5qHH3qcs+t+/dN+4+dNTn5wJHtMv/Ydq/66W8n2fs83zs58V8xD/I1A8vBS5plz8NPC9J2vJPVNUDVXUr8J328xa6Yfq8GO23v1V1W1VdDzy8R9tfBdZX1Y+r6m5gPXDWKIIe0jB9XqwG6fPVVfX37epGmvt/YGnv59n6PLTFnPifBNw+bf2OtmzGOlX1IPD/gJ8asO1CNEyfAY5P8o0kX03yy10HOw+G2U9LeR/vy2FJJpNsTHL2/IbWmQPt83nAF+bYdqEYps8w5H5ezPPxDzL1w2x1Bpo2YgEaps/bgVVV9aMkpwGfS/K0qrpnvoOcR8Psp6W8j/dlVVXdmeQE4CtJtlTVzfMUW1cG7nOSfwdMAM8+0LYLzDB9hiH382I+4h9k6oef1ElyMPA44McDtl2I5tzndljrRwBVtYlmfPHnOo94OMPsp6W8j2dVVXe2P28BNgCnzGdwHRmoz0nOBN4JvKSqHjiQtgvQMH0efj+P+yTHECdHDqY5kXM8/3Ry5Gl71HkdjzzReVm7/DQeeXL3FhbHyd1h+rxyqo80J5S+Dxw97j4N299pdS9m75O7t9Kc8DuqXV7Q/Z2HPh8FPKpdXgHcxB4nDBfia8B/16fQHKycuEf5kt3P++jz0Pt57L+AIX95LwT+rv3lvLMt+0Oab0eAw4BP0Zy8vRY4YVrbd7btvg28YNx96brPwL8Cbmj/gV0HvHjcfZmn/j6T5uhpF/Aj4IZpbX+j/T18B3jNuPvSdZ+BXwS2tPt4C3DeuPsyj33+EvADYHP7urIH+3nGPs/HfnbKBknqmcU8xi9JmgMTvyT1jIlfknrGxC9JPWPil6SeMfFrwZk28+DWJJ9KckSH2/rD9iYZkrzxQLeVxleSPLZdf0OSbUkunYfYfj3JE6etfzDJSXP8rNcnec2wMWlp8HJOLThJ7quqx7TLlwKbqur9A7QLzb/pOU1eluQ2YKKqfngAbV4EnFlVb2rXb6S5L+TWPeodXM3cSQcSzwbgd6tq8kDazfJZRwD/p6oWw5286phH/Fro/hr4WfjJ/ORb29cb27LV7RH2n9LcmHZcknOSbGnrvbettyzJxW3ZliRTifriJC9P8gbgicDVSa5Ocl6SC6aCSPLaJDN9+bwKuKKtcyHNXdFXJnlTkvOTrEvyReAjbax/3c6lf12SX5z2+b/XxvXNJO9JM8/+BHBp+9fP4Uk2JJlo6+/Vx7b8viTvbj9nY5InAFQzy+NtSRbDLLTq2rjvXvPla88X7VzjNLe1X0EzF/lpNHcpPhp4DM1dyKfQzE3/MLCmbfNE4Hs0U1QcDHwFOLttv37aNpa3Py+mnfYAuA1Y0S4/muaOykPa9f8LPH2GWL8LHDltffpnnA9sAg5v148ADmuXTwQm2+UXtJ8/Nff60e3PDTR/gTB9fbY+tnWK9q5s4L8C75rW/p3AW8a9f32N/+URvxaiw5NsBiZpEtxFwC8Bl1fVrqq6D/gsMDW19HeramO7/ExgQ1XtqGZo5VLgV2jmRTkhyf9Ichawz1lJq2oXTUL9tSRPpfkC2DJD1aOr6t59fNSVVXV/u3wI8IEkW2im1Zgarz8T+HC1c69X1Y/3Fds++giwG7iqXd5E88U45S6aLw313GKelllL1/1VdfL0gv08TGbX9KozVaiqu5M8g+bBHa8DXkEzx8u+fBB4B3Aj8OFZ6jyY5KCa/bzC9NjeRDP3yjNohln/YVrMB3KybV+/i3+sqqnPeohH/h8/DLh/7ybqG4/4tVh8DTg7yRFJHg28jGb8f0/XAM9OsiLJMuAc4KtJVgAHVdVngP8EnDpD23uBI6dWquoamqlz/y3w8Vni+jbNuP4gHgdsb78kXk3z+D2ALwK/MXVFUZKjZ4pnf30cYPs/B2wdMFYtYSZ+LQpVdR3NePy1NInvg1X1jRnqbQfeDlxNOxNpVV1B83SjDe0Q0sVtnT2tA76Q5OppZZfRXA1z9yyh/QXwnAG78afAuUk20iThXW3MfwlcCUy28f1uW/9i4MKpk7sD9HF/nkUz46N6zss5pX1IchVwQVV9eZb3jwE+UlXPH21kBybJKcCbq+rV445F4+cRvzSDJMuT/B3N+YYZkz785Oj7A1M3cC1gK2iGuCSP+CWpbzzil6SeMfFLUs+Y+CWpZ0z8ktQzJn5J6pn/D4dZSwf7tikmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(X,color='red',alpha=0.2,edgecolor='black',bins=np.linspace(fmin,fmax,20))\n",
    "plt.xlabel(feature + ' ' + funits); plt.ylabel('Frequency');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### The Arithmetic Average / Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x} = \\frac{1}{n}\\sum^n_{i=1} x_i\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity average is 0.15 (fraction).\n"
     ]
    }
   ],
   "source": [
    "average = np.average(X)\n",
    "print(feature + ' average is ' + str(round(average,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### The Weighted Arithmetic Average / Mean\n",
    "\n",
    "Many of the following methods accept data weights, e.g. declustering\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x} = \\frac{\\sum^n_{i=1} \\lambda_i x_i}{\\sum^n_{i=1} \\lambda_i}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity average is 0.15 (fraction).\n"
     ]
    }
   ],
   "source": [
    "weights = np.ones(X.shape)\n",
    "wt_average = np.average(X,weights = weights)\n",
    "print(feature + ' average is ' + str(round(wt_average,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### Median\n",
    "\n",
    "\\begin{equation}\n",
    "P50_x = F^{-1}_{x}(0.50)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity median is 0.15 (fraction).\n"
     ]
    }
   ],
   "source": [
    "median = np.median(X)\n",
    "print(feature + ' median is ' + str(round(median,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### Mode\n",
    "\n",
    "The most common value. To do this we should bin the data, like into histogram bins/bars.  To do this we will round the data to the 2nd decimal place.  We are assume bin boundaries, $0.01, 0.02,\\ldots, 0.30$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity mode is 0.14 (fraction).\n"
     ]
    }
   ],
   "source": [
    "mode = stats.mode(np.round(X,2))\n",
    "print(feature + ' mode is ' + str(round(mode,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### Geometric Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_G = ( \\prod^n_{i=1} x_i )^{\\frac{1}{n}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity geometric mean is 0.15 (fraction).\n"
     ]
    }
   ],
   "source": [
    "geometric = scipy.stats.mstats.gmean(X)\n",
    "print(feature + ' geometric mean is ' + str(round(geometric,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Central Tendency\n",
    "\n",
    "##### Harmonic Mean\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_H = \\frac{n}{\\sum^n_{i=1} \\frac{1}{x_i}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity harmonic mean is 0.14 (fraction).\n"
     ]
    }
   ],
   "source": [
    "hmean = scipy.stats.mstats.hmean(X)\n",
    "print(feature + ' harmonic mean is ' + str(round(hmean,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Power Law Average\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{x}_p = (\\frac{1}{n}\\sum^n_{i=1}{x_i^{p}})^\\frac{1}{p}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity power law average for p = -0.5 is 0.14 (fraction).\n"
     ]
    }
   ],
   "source": [
    "power = -0.5\n",
    "power_avg = np.average(np.power(X,power))**(1/power)\n",
    "print(feature + ' power law average for p = ' + str(power) + ' is ' + str(round(power_avg,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's Visualize Some of the Measures of Central Tendency\n",
    "\n",
    "To visualize ocmpare these statistics, parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXQUZdb48e8dCCQBZAsKyhL4uYwsYQuIKJhBBIYRB0Yc4H1FER0Y0aO4juKGuCMq4gLiAo74AgqCouiwaARE1AQBUUQR2RRkEVB2Avf3R1d6GujudNKprk73/ZzTJ9WVW9230nBTeeqpW6KqGGOMSR5/8DoBY4wxsWWF3xhjkowVfmOMSTJW+I0xJslY4TfGmCRT3usEIpGRkaGZmZlep2FMXFi9ejX79+8nLS2Ns846y+t0TBzLz8/frqq1jl9fJgp/ZmYmeXl5XqdhjGsWL14MQPv27dm4cSMA9erVCxqbk5PDsmXLaNGiBbm5ubFK0ZRBIrI+2PoyUfiNSXTDhg0DIDc3l/79+/uXjXGDFX5j4szdd99dZEyDBg0iijMmGCv8xsSZzp07FxlTvXr1iOKMCabMFv7Dhw+zadMmDhw44HUqZUZqaip169YlJSXF61RMGGvXrgWgUaNGIWP27NnjH+c3prjKbOHftGkTVapUITMzExHxOp24p6rs2LGDTZs20bBhQ6/TMWEMHDgQCD/Gv2bNGoYOHWrnAUyJlNnCf+DAASv6xSAi1KxZk23btnmdiinC/fff73UKJsGV2cIPWNEvJvt5lQ0XXHCB1ymYBGdX7hoTZ1avXs3q1au9TsMkMCv8UZoxYwYiwrfffut1KiZBDB48mMGDB3udhklgZXqoJx5MnjyZ888/nylTpjB8+PCoXuvIkSOUK1eudBIzceGrJUs4tGtXkXFXXHIJAPkffHDMcoVq1WjWrt0J8Q0bNuThhx8u3WRN0rDCH4U9e/bwySef8NFHH3HJJZcwfPhw+vTpw5VXXkn37t0BGDBgAD169KBnz57ccccd5ObmcvDgQa677joGDx5Mbm4u999/P3Xq1GHZsmV888039OzZk40bN3LgwAFuvPFGBg0aBMDLL7/MY489xqmnnsoZZ5xBxYoVefbZZ9m2bRv//Oc/2bBhAwCjR4/mvPPO8+znYv7r0K5dtK51QquUE7QOGNcPXM4PcTK+atWqtG/fPvoETVJyvfCLSDkgD/hJVS8WkYbAFKAGsBTor6qHonmPoUOHsmzZsuiTDdCiRQtGjx4dNmbmzJl069aNM888kxo1arB06VL69u3L1KlT6d69O4cOHWL+/PmMHTuWl19+mapVq/LFF19w8OBBzjvvPLp06QLA559/zsqVK/3TLF955RVq1KjB/v37adOmDZdeeikHDx7kgQceYOnSpVSpUoVOnTrRvHlzAG688UZuuukmzj//fDZs2EDXrl1ZtWpVqf48jLsWL18OQPvmzVm5Zg0ATU8/PWT87t27Wbx4sRV/UyKxOOK/EVgFnOQ8fwx4SlWniMg44GpgbAzyKHWTJ09m6NChAPTt25fJkyfzwAMPcMMNN3Dw4EE++OADOnbsSFpaGnPmzGHFihVMmzYN8P3H/f7776lQoQJt27Y9Zm79mDFjmDFjBgAbN27k+++/Z8uWLVxwwQXUqFEDgMsuu4zvvvsOgHnz5vHNN9/4t//tt9/4/fffqVKlSkx+DiZ6w557DoDc8eO5fuRI/3IoP/74I8OGDbN5/KZEXC38IlIX+AvwEHCz+OYTdgL+xwl5FRhOlIW/qCNzN+zYsYMPP/yQlStXIiIcOXIEEWHkyJHk5OTwn//8h6lTp9KvXz/AdwHVM888Q9euXY95ndzcXCpVqnTM83nz5vHpp5+Snp5OTk4OBw4cQFVD5nL06FE+/fRT0tLS3NlZ47oXnCZtAI/feKOHmZhk4PasntHA7cBR53lNYJeqFjjPNwGnBdtQRAaJSJ6I5MXjRUfTpk3jiiuuYP369axbt46NGzfSsGFDFi1aRN++fZkwYQILFy70F/quXbsyduxYDh8+DMB3333H3r17T3jd3bt3U716ddLT0/n2229ZsmQJAG3btuXjjz9m586dFBQUMH36dP82Xbp04dlnn/U/L+1hL+O+szIzOcu550SbJk1o06SJtwmZhOZa4ReRi4GtqpofuDpIaNBDWVUdr6rZqppdK4KTY7E2efJkevXqdcy6Sy+9lP/7v/+jS5cuLFiwgM6dO1OhQgUArrnmGho3bkyrVq1o2rQpgwcPpqCg4ITX7datGwUFBWRlZXHPPffQzpnRcdpppzFs2DDOOeccOnfuTOPGjalatSrgGxrKy8sjKyuLxo0bM27cOJf33pS2WQsWMGvBAgCWrV7NMpvHb1wk4YYQonphkUeA/kABkIpvjH8G0BWoraoFInIuMFxVu4Z+JcjOztbjb8SyatUqzj77bFdyj1d79uyhcuXKFBQU0KtXLwYOHHjCL5+iJOPPzUv5H3wQ0ayeHGfmVu748ccs52/bRutu3Y6NtRuxmAiJSL6qZh+/3rUxflW9E7jTefMc4FZV/V8ReRPojW9mz5XA227lkGiGDx/OvHnzOHDgAF26dKFnz55ep2RcMPqWW4qMOf300z05t2USgxfz+P8FTBGRB4EvgZc9yKFMGjVqlNcpmBhoEcF9dCtXrmwtmU2JxaRlg6rmqurFzvJaVW2rqqer6mWqejAWORhTVnzx9dd88fXXYWN27tzJvHnzYpSRSTTWq8eYOHPb009z29NPh41Zv349Dz74YIwyMonGWjYYE2eevf12r1MwCc4KvzFxJlyrBmNKgw31REFE6N+/v/95QUEBtWrV4uKLLy7W6+Tk5FA4XbV79+7siqCbo0lci5cv9/fuMcYNdsQfhUqVKrFy5Ur2799PWloac+fO5bTTgl6IHLHZs2eXUnamrArs22OMG+yIP0p//vOfee+99wDf1byFvXkA9u7dy8CBA2nTpg0tW7bk7bd9lyzs37+fvn37kpWVRZ8+fdi/f79/m8zMTLZv3w5Az549ad26NU2aNGF8QBGoXLkyd911F82bN6ddu3b88ssvsdhV46IXhg3z9+sJXA7lzDPP5IUXXohFaiYBJcwRf05Ozgnr/v73vzNkyBD27dvn748faMCAAQwYMIDt27fTu3fvY74X6RWRffv2ZcSIEVx88cWsWLGCgQMHsnDhQgAeeughOnXqxCuvvMKuXbto27YtnTt35oUXXiA9PZ0VK1awYsUKWrVqFfS1g7VnrlmzJnv37qVdu3Y89NBD3H777bz44ovcfffdEeVr4lNhn57jl0NJT0/nrAjm+xsTjB3xRykrK4t169YxefLkE365zJkzh0cffZQWLVr4u2xu2LCBBQsWcPnll/u3z8rKCvraY8aM8R/VF7ZnBqhQoYL/PELr1q1Zt26deztoYiKwV8/H+fl8nJ8fNn779u3MmjUrFqmZBJQwR/zhjtDT09PDfj8jIyOqnieXXHIJt956K7m5uezYscO/XlWZPn160CMzX4fq0EK1ZwZISUnxb1+uXLmgzd5M2fLEpEkA9OjYkfucIZxwY/ybNm3iiSeeoEePHjHJzySWhCn8Xho4cCBVq1alWbNmx/wC6dq1K8888wzPPPMMIsKXX35Jy5Yt6dixI6+//jp/+tOfWLlyJStWrDjhNUO1ZzaJaZpz8xWAV+6918NMTDKwoZ5SULduXW4McvOMe+65h8OHD5OVlUXTpk255557ALj22mvZs2cPWVlZjBw5krZt256wbaj2zCYxZVSrRka1agA0qluXRnXrepyRSWR2xB+FPXv2nLAuJyfHf6I5LS0t6MyLtLQ0pkyZEvQ1A8fr33///SLft3fv3iecmDZlz0RnvH5Ajx7M++wzADqfc46XKZkEZoXfmDgQWPgffNnXsNYKv3GLFX5j4sxrI0YUGfPHP/6R1157LQbZmERkhd+YOFOvdu0iY1JTU6lXr14MsjGJyE7uGhNnPli8mA8WLw4bs3XrVqZOnRqjjEyicfNm66ki8rmILBeRr0Xkfmf9RBH5UUSWOQ+7jZAxAR6dOJFHJ04MG/Pzzz8zduzY2CRkEo6bQz0HgU6qukdEUoBFIlI4TeU2VZ3m4nsbU2ZNefhhr1MwCc61I371KZx3mOI81K3380JptWUObMxmTO2MDGpnZHidhklgro7xi0g5EVkGbAXmqupnzrceEpEVIvKUiFQMse0gEckTkbxt27a5mWaJBbZlBkqlLbMxgX17jHGDq4VfVY+oagugLtBWRJoCdwJ/BNoANYB/hdh2vKpmq2p2rVq13EwzKuHaMv/666/07NmTrKws2rVr52/NsGPHDrp06ULLli0ZPHgwqv/9Q2jSpEm0bduWFi1aMHjwYI4cORLbHTKee2LSJH/vHmPcEJPpnKq6S0RygW6qOspZfVBEJgC3lsZ7BGvLfLyLL76YW2+91R/vdlvm++67j5YtWzJz5kw+/PBDrrjiCpYtW8b999/P+eefz7333st7773n77W/atUqpk6dyieffEJKSgpDhgzh9ddf54orroj8B2HKpMBePYHLoTRp0oRp0+w0mSkZ1wq/iNQCDjtFPw3oDDwmInVUdbP42kv2BFa6lUMshGvLvGjRIqZPnw5Ap06d2LFjB7t372bBggW89dZbAPzlL3+hevXqAMyfP5/8/HzatGkD+G7YcvLJJ8dwb4xXCvv0HL8cSkpKChl2HsCUkJtH/HWAV0WkHL4hpTdU9V0R+dD5pSDAMuCfpfFmxW2rHBjvZlvm4xW2Uw7WlllVufLKK3nkkUdKnIspmwJbNrz14YcA/K1Tp5DxW7ZsYeLEiQwYMCAW6ZkE4+asnhWq2lJVs1S1qaqOcNZ3UtVmzrrLA2b+lFkDBw7k3nvvpVmzZsesL2y/DL5fNBkZGZx00knHrH///ffZuXMnABdeeCHTpk1j69atgO8cwfr162O4J8YrE2fN8hf/MVOmMCZEE79ChYXfmJKwlg2lIFRb5uHDh3PVVVeRlZVFeno6r776KuAb++/Xrx+tWrXiggsuoH79+gA0btyYBx98kC5dunD06FFSUlJ47rnnaNCgQUz3x8Re4E1X3n7ySQ8zMcnACn8UimrLXKNGDf8N1gPVrFmTOXPm+J8/9dRT/uU+ffrQp0+f0k/WlBlVK1f2OgWT4KxXjzFxYNRrrzHK6bY5dc4cpgYcGBhT2qzwGxMH3l24kHedacBjp01jrE3VNC6yoR5j4szsMWOKjGnWrBmzZ8+OQTYmEVnhNybOpKemFhlTrlw50tPTY5CNSUQ21GNMnJk0ezaTijia/+mnn3j++edjlJFJNFb4jYkzL82cyUszZ4aN2bZtG2+88UaMMjKJxoZ6olCuXDmaNWtGQUEBZ599Nq+++qr9+W2iNteO5I3L7Ig/CmlpaSxbtoyVK1dSoUIFxo0b5/p7WrfOxJdSvjwp5e2YzLjHCn8p6dChA2vWrAHgySefpGnTpjRt2pTRo0cDMHLkSMY4szVuuukmOjl9WObPn8/ll18OwJw5czj33HNp1aoVl112mf8CsczMTEaMGMH555/Pm2++GetdMzEW2L7BGDckTOHPycnx9y45fPgwOTk5THJ6mu/bt4+cnBz/zal3795NTk6Ov0Pm9u3bycnJYZbzn23Lli3Feu+CggLef/99mjVrRn5+PhMmTOCzzz5jyZIlvPjii3z55Zd07NjR3645Ly+PPXv2cPjwYRYtWkSHDh3Yvn07Dz74IPPmzWPp0qVkZ2fzZMCl+6mpqSxatIi+fftG9XMy8c8Kv3Gb/T0Zhf3799Oihe9e8R06dODqq69m7Nix9OrVi0qVKgHwt7/9jYULF3LttdeSn5/P77//TsWKFWnVqhV5eXksXLiQMWPGsGTJEr755hvOO+88AA4dOsS5557rfy9r45DYAnv1BC6H0qJFi6g6yprkljCFP/A/QUpKyjHP09PTj3letWrVsG2Za9euHdF7Fo7xBwrWirkwp8zMTCZMmED79u3Jysrio48+4ocffuDss8/mhx9+4KKLLmLy5MlBty/8RWKMMdFKmKGeeNGxY0dmzpzJvn372Lt3LzNmzKBDhw7+740aNYqOHTvSoUMHxo0bR4sWLRAR2rVrxyeffOI/T7Bv3z6+++47L3fFxFBgr54XZ8zgxRkzwsZv3LiRUaNGhY0xJhQr/KWsVatWDBgwgLZt23LOOedwzTXX0LJlS8A3HLR582bOPfdcTjnlFFJTU/2/FGrVqsXEiRPp16+f/x693377rZe7YmLo0xUr+NS5J/PUuXOZOndu2PgdO3bw7rvvxiI1k4DcvPViKrAAqOi8zzRVvU9EGgJT8N1ofSnQX1UPuZWHm4K1ZQa4+eabufnmm09Yf+GFF3L48GH/8+OP6Dt16sQXX3xxwnbr1q2LLlET96Y//rh/eZ7N4zcuc/OI/yDQSVWbAy2AbiLSDngMeEpVzwB2Ale7mIMxxpjjuHnrRQ24rWKK81CgE1DYc/ZVfDdcNyap3fnss9z57LMAPP/mmzxv12sYF7k6q8e50Xo+cDrwHPADsEtVC5yQTcBpIbYdBAwC/LcmPJ6qBr1puQku1Iwj473C8X2AWc71HkMuu4wfVq06Ifb3X39Fjxzh8N695H/wQdjXrVCtGs3atSvdZE2Z52rhV9UjQAsRqQbMAM4OFhZi2/HAeIDs7OwTYlJTU9mxYwc1a9a04h8BVWXHjh2kRtDy13jr/YB+/Ef37qV1rVrHfL9KSgotzzyT3AhahORv21bq+ZmyLybz+FV1l4jkAu2AaiJS3jnqrwv8XJLXrFu3Lps2bWKb/cOOWGpqKnXr1vU6DWOMx9yc1VMLOOwU/TSgM74Tux8BvfHN7LkSOPFu5BFISUmhYcOGpZWuMXHjaecivhv79QsZs37zZh546SXuueaaWKVlEoibs3rqAB+JyArgC2Cuqr4L/Au4WUTWADWBl13MwZgyZ/7nnzP/88/Dxuz8/fciY4wJxbUjflVdAbQMsn4t0Nat9zWmrHvnqae8TsHvqyVLOLRrV4m3t5PL8SlhevUYY0rfoV27Tji5XBx2cjk+WcsGY+JMYN8eY9xgR/zGxIGaVav6lwPn9IeSUr78MdsYUxxW+I2JA4G9egKXQ2nSqFFEccYEY4XfmAQW7Mrf4vhx1aqoxvhNfLLCb0wcKOzT88j11/OocwvROwYMCBm/9qefuPPZZ3nk+uvDvm6wK3+L4/u8vBJva+KXFX5j4sCO3bv9y8tWry4y/re9eyM6F2BMMFb4jYkD4++6y7885ZFHPMzEJAObzmmMMUnGCr8xcWDQQw8x6KGHAHjgpZd44KWXPM7IJDIb6jEmDny3fr1/eXXAcigVU1Koe8opbqZkEpgVfmPizKQHHigy5uyGDSOKMyYYG+oxxpgkE1HhF5GmbidijPG5d9w47i3i7lprNm5k6BNPxCgjk2giHeoZJyIVgInA/6lqyfu0GmPC2vjLL0XG7Nm/P6L5/sYEE1HhV9XzReQMYCCQJyKfAxNUda6r2RmThCbcd5/XKZgEF/EYv6p+D9yN7w5aFwBjRORbEflbsHgRqSciH4nIKhH5WkRudNYPF5GfRGSZ8+heGjtijDEmMhEd8YtIFnAV8BdgLtBDVZeKyKnAp8BbQTYrAG5x4qoA+SJS+BfCU6o6Kvr0jUk8gX17jHFDpGP8zwIvAsNUdX/hSlX9WUTuDraBqm4GNjvLv4vIKuC0KPM1JiGd2aCBfzmwb08o6RUrHrONMcURaeHvDuxX1SMAIvIHIFVV96lqkbcKEpFMfPff/Qw4D7heRK4A8vD9VbCzBLkbkzACe/UELodyZoMGEcUZE0ykY/zzgLSA5+nOuiKJSGVgOjBUVX8DxgL/D2iB7y+CoHPSRGSQiOSJSN42u2+nMcaUmkgLf6qq7il84iynF7WRiKTgK/qvq+pbzra/qOoRVT2Kb/iobbBtVXW8qmaranYtuxGESXCBvXpuHT2aW0ePDhv/3fr1/nhjiivSwr9XRFoVPhGR1sD+MPGIiAAvA6tU9cmA9XUCwnoBKyNP15jEVLNqVf89dPcfPMj+gwfDxu87ePCY/j7GFEekY/xDgTdF5GfneR2gTxHbnAf0B74SkWXOumFAPxFpASiwDhhcrIyNSUCBM3ie+9e/PMzEJINIL+D6QkT+CJwFCPCtqh4uYptFTuzxZhc7S2OMMaWmON052wCZzjYtRQRV/bcrWRmTZC697TYApj/+uL8Hz+hbbvEyJZPAIr2A6zV8M3GWAUec1QpY4TemFEQydz9Q5bQ0Wpx1lkvZmEQX6RF/NtBYVdXNZIwxkR3pn16vnv1FYEos0lk9K4HabiZijDEmNiIt/BnANyLyHxF5p/DhZmLGJKvrHnuM6x57LGzMqh9/5PJ77olRRibRRDrUM9zNJIwx/5VWsWKRMQcPH2ZTBH37jQkm0umcH4tIA+AMVZ0nIulAOXdTMyY5jRo61OsUTIKL9NaL/wCmAS84q04DZrqVlDHGGPdEOsZ/Hb4rcX8D/01ZTnYrKWOSWWDfHmPcEOkY/0FVPeRrvwMiUh7fPH5jTCk4NyvLv1zYsyeckypVOmYbY4oj0sL/sYgMA9JE5CJgCDDLvbSMSS6BvXoiufNWo9NOszt0mRKLdKjnDmAb8BW+pmqz8d1/1xhjTBkTUeFX1aOq+qKqXqaqvZ1lG+oxppRcettt/n49V91/P1fdf3/Y+K/XrvXHG1Nckfbq+ZEgY/qq2qjUMzImCQWO19c75ZQi4w8XFBS7v48xhYrTq6dQKnAZUKP00zEmOd3av79/ecQ//+lhJiYZRDrUsyPg8ZOqjgY6uZybMcYYF0Q61NMq4Okf8P0FUMWVjIxJQjmDBgGQO368vwfPpAce8DIlk8AiHep5ImC5AN8tE/8ebgMRqYevX39t4CgwXlWfFpEawFR8N3VZB/xdVXcWK2tjEthZDRoUGVO9ShUubNs2BtmYRBRpr54/leC1C4BbVHWpiFQB8kVkLjAAmK+qj4rIHfimitpNRo1x3HPNNUXGNKhTJ6I4Y4KJdKjn5nDfV9Ung6zbDGx2ln8XkVX4evz8Fchxwl4FcrHCb4wxMRPpBVzZwLX4CvdpwD+BxvjG+Ysc6xeRTKAl8BlwivNLofCXQ9CePyIySETyRCRv27ZtEaZpTNnX98476XvnnWFjVqxZw59vuCFGGZlEE+kYfwbQSlV/BxCR4cCbqlrk35oiUhmYDgxV1d8K+/0URVXHA+MBsrOz7WIxkzQiuZfu0aNH2X/gQAyyMYko0sJfHzgU8PwQvpOzYYlICr6i/7qqvuWs/kVE6qjqZhGpA2wtRr7GJLw7BgzwOgWT4CIt/K8Bn4vIDHxX8PbCN2MnJPEd2r8MrDruHMA7wJXAo87Xt4ubtDHGmJKLdFbPQyLyPtDBWXWVqn5ZxGbnAf2Br0RkmbNuGL6C/4aIXA1swHcVsDHGUdiDZ/rjj3uciUlUkR7xA6QDv6nqBBGpJSINVfXHUMGquggINaB/YXGSNCbRXdyhg385kj77NatWPWYbY4oj0umc9+Gb2XMWMAFIASbhO6o3xkQpsFdP4HIo9U45JaI4Y4KJdDpnL+ASYC+Aqv6MtWwwxpgyKdLCf8jpv68AIlLJvZSMST45gwb5+/VcctNNXHLTTWHjl333nT/emOKKdIz/DRF5AagmIv8ABgIvupeWMcllQI8e/mXrwWPcFumsnlHOvXZ/wzfOf6+qznU1M2OSSGDhv7FfPw8zMcmgyMIvIuWA/6hqZ8CKvTEu2L5rFwAZ1ap5nIlJBkUWflU9IiL7RKSqqtq93oxxQe/bbwd8/fgLe/C8P2aMlymZBBbpGP8BfBdizcWZ2QOgqtYlyphS1iOC+fm1qlfn7xddFINsTCKKtPC/5zyMMS4bclnRF7OfVqtWRHHGBBO28ItIfVXdoKqvxiohY0zRjhw9yr4DB0hPTfU6FVMGFTWPf2bhgohMdzkXYwzQecgQOg8ZEjbmqzVr6G79+E0JFTXUE9hrp5GbiRhjfPrY2L1xWVGFX0MsG2Nc8o9evbxOwSS4ogp/cxH5Dd+Rf5qzjPNcVfUkV7MzxhhT6sIWflUtF6tEjDE+hT14cseP9zgTk6iK04/fGOOSwJYNgcuh1K5ZM6I4Y4JxrfCLyCvAxcBWVW3qrBsO/APY5oQNU9XZbuVgTFlhhd/EUqRtmUtiItAtyPqnVLWF87Cibwy+Xj2F/XoOFxRwuKAgbPzhggJ/vDHF5VrhV9UFwK9uvb4xiaT37bf7+/VcNGQIFxUxj//rtWv98cYUlxdj/NeLyBVAHnCLqu4MFiQig4BBAPXr149hesbE3i2XX+5fvqZnTw8zMckg1oV/LPAAvmsCHgCewHdTlxOo6nhgPEB2drZdQ2ASWo+OHf3Ll3fv7mEmJhm4OcZ/AlX9RVWPqOpRfHfwslsNGQOsXreO1evWAbDvwAH2HTjgbUImocX0iF9E6qjqZudpL2BlLN/fmHg1+OGHAd/c/cIePDaP37jFzemck4EcIENENgH3ATki0gLfUM86YLBb729MWXVt795FxpyakRFRnDHBuFb4VTXYjUNfduv9jAnmqyVLOBTFtMcK1arRrF27UsyoaH26dCky5uQaNSKK89oPq1ZFtb0XP/9kYFfumoR2aNcuWteqVeLt87dtKzqolO3esweAqpUrh4w5cOgQG7dsoV7t2rFKq0SO7t1b5n7+ySCmJ3eNMUX7680389ebbw4b8+26dfS/994YZWQSjR3xGxNnbujb1+sUTIKzwm9MnPlbp05ep2ASnBV+Y8KI9uTkj6tWFXuMu7AHT0a1alG9tzGhWOE3JoxoT05+n5dX7G0Ke/DYPH7jFiv8xsSBwF49gcuh1D355IjijAnGCr8xcSCwV0/gcigZ1apFFGdMMDad05g4ENirZ8v27WzZvj1s/L4DB/zxxhSXFX5j4sDghx/29+vpO2wYfYcNCxv/3YYN/nhjisuGeoyJAw9fd51/+Y4BA7xLxCQFK/zGxIH2zZv7l7u1b+9hJiYZ2FCPMXFg8fLlLF6+HICNW7awccsWjzMyicyO+FTg+WQAAA4BSURBVI2JA8Oeew7wzd0v7MFj8/iNW6zwGxNn7r766iJjGtSuHVGcMcFY4TcmznQ+55wiY6qfdFJEccYE49oYv4i8IiJbRWRlwLoaIjJXRL53vlZ36/2NKavWbtrE2k2bwsbs2bePZatXxygjk2jcPLk7Eeh23Lo7gPmqegYw33lujAkwcMQIBo4YETZmzaZNDH3iiRhlZBKNm7deXCAimcet/iu++/ACvArkAv9yKwdjyqL7B9utqI27Yj3Gf4qqbgZQ1c0icnKoQBEZBAwCqF+/fozSM8Z7F7Ru7XUKJsHF7Tx+VR2vqtmqml0rira4xpQ1gX17jHFDrI/4fxGROs7Rfh1ga4zf35i4V9iDx+bxG7fEuvC/A1wJPOp8fTvG729MXArs1RO4HErDU0+NKM6YYFwr/CIyGd+J3AwR2QTch6/gvyEiVwMbgMvcen9jypLAXj2By6FUrVw5ojhjgnFzVk+/EN+60K33NKasKuzT0755c1auWQNA09NPDxm/e88eFi9fbsXflEjcntw1JpkMe+45f7+e60eO5PqRI8PG//jzz/54Y4rLWjYYEwdeCLjxyuM33uhhJiYZWOE3Jg6clZnpX27TpIl3iZikYEM9xsSBWQsWMGvBAgCWrV5tfXiMq+yI35g48MSkSQD06NjR34PH5vEbt1jhNybOjL7lliJjTq9bN6I4Y4Kxwm9MnGlx1llFxlROT48ozphgbIzfmDjzxddf88XXX4eN2fnbb8z77LMYZWQSjRV+Y+LMbU8/zW1PPx02Zv2WLTz48ssxysgkGhvqMa76askSDu3aVeLtN2zcSP169Uq8/Y+rVtG6jHV3ffb2271OwSQ4K/zGVYd27Yqq8H6fl0frVq2i2r6sCdeqwZjSYEM9xsSZxcuX+3v3GOMGO+I3Js4U9uCxefzGLVb4jYkDgb16ApdDObN+/YjijAnGCr8xcSCwV0/gcijpqakRxZV1P6xaFdX20U4OqFCtGs3atYsqh3hkhd+YOFDYp6dHx458nJ8PhL/p+vZdu5i1YAE9OnaMSX5eObp3r6eTA/K3bSvxtvHMk8IvIuuA34EjQIGqZnuRhzHxIrBXz30vvACEH+PftHUrT0yalPCF37jDyyP+P6nqdg/f35i4MS3gxiuv3Huvh5mYZGBDPcbEgYxq1fzLjerW9TATkwy8msevwBwRyReRQR7lYEzcmDhrFhNnzQJg3mefWR8e4yqvjvjPU9WfReRkYK6IfKuqCwIDnF8IgwDq16/vRY7GxExh0R/Qo4e/B0/nc87xMiWTwDwp/Kr6s/N1q4jMANoCC46LGQ+MB8jOztaYJ2mMR14bMaLImD9mZkYUZ0wwMR/qEZFKIlKlcBnoAqyMdR7GxKt6tWtTr3btsDGpFSoUGWNMKF6M8Z8CLBKR5cDnwHuq+oEHeRgTlz5YvJgPFi8OG7P111+ZOmdOjDIyiSbmQz2quhZoHuv3NaaseHTiRAC6tW8fMubn7dsZO20afbp0iVFWJpHYdE5j4syUhx/2OgWT4KzwGxNnamdkeJ2CSXDWj9+YODNrwQJ/7x5j3GBH/MbEmcC+Pca4wQq/MXEgsFdP4HIoTRo1iijOmGCs8Juwor1Zelm82bkXAnv1BC6HklK+fERxxgRjhd+EVRo3SzdFC2zZ8NaHHwLwt06dQsZv2bGDibNmMaBHj5jkZxKLndw1Jg4ENmkbM2UKY6ZMCRtfWPiNKQk74jcmDgTedOXtJ5/0MBOTDKzwGxNnqlau7HUKJsHZUI8xcWDUa68x6rXXAJg6Z4714TGussJvTBx4d+FC3l24EICx06Yxdto0jzMyicyGeoyJM7PHjCkyptnpp0cUZ0wwVviNiTPpqalFxpT7wx8iijMmGBvqMSbOTJo9m0mzZ4eN+WnbNp5/880YZWQSjR3xx7lor5zdsHEj9evVK/H2duVt7L00cyYAl3fvHjJm286dvDF3LkMuuyxWaZkEYoU/zpXGlbOtW7WKansTW3Off97rFEyC82SoR0S6ichqEVkjInd4kYMx8SqlfHlSytsxmXGPFzdbLwc8B/wZaAz0E5HGsc7DmHgV2L7BGDd4ccTfFlijqmtV9RAwBfirB3kYE5es8Bu3iarG9g1FegPdVPUa53l/4BxVvf64uEHAIOfpWcDqEr5lBrC9hNuWVbbPycH2OTlEs88NVPWEk4ReDCRKkHUn/PZR1fHA+CCxxXszkTxVzY72dcoS2+fkYPucHNzYZy+GejYBgfML6wI/e5CHMcYkJS8K/xfAGSLSUEQqAH2BdzzIwxhjklLMh3pUtUBErgf+A5QDXlHVr118y6iHi8og2+fkYPucHEp9n2N+ctcYY4y3rFePMcYkGSv8xhiTZMp04S+q9YOIVBSRqc73PxORzIDv3emsXy0iXWOZdzRKus8ikiki+0VkmfMYF+vcSyKC/e0oIktFpMC5RiTwe1eKyPfO48rYZR2dKPf5SMBnXGYmTUSwzzeLyDciskJE5otIg4DvJernHG6fo/ucVbVMPvCdGP4BaARUAJYDjY+LGQKMc5b7AlOd5cZOfEWgofM65bzeJ5f3ORNY6fU+uLC/mUAW8G+gd8D6GsBa52t1Z7m61/vk5j4739vj9T64tM9/AtKd5WsD/l0n8uccdJ9L43Muy0f8kbR++CvwqrM8DbhQRMRZP0VVD6rqj8Aa5/XiXTT7XBYVub+quk5VVwBHj9u2KzBXVX9V1Z3AXKBbLJKOUjT7XFZFss8fqeo+5+kSfNf/QGJ/zqH2OWplufCfBmwMeL7JWRc0RlULgN1AzQi3jUfR7DNAQxH5UkQ+FpEObidbCqL5nBL5Mw4nVUTyRGSJiPQs3dRcU9x9vhp4v4Tbxoto9hmi/JzLcu/XSFo/hIqJqG1EHIpmnzcD9VV1h4i0BmaKSBNV/a20kyxF0XxOifwZh1NfVX8WkUbAhyLylar+UEq5uSXifRaRy4Fs4ILibhtnotlniPJzLstH/JG0fvDHiEh5oCrwa4TbxqMS77MzrLUDQFXz8Y0vnul6xtGJ5nNK5M84JFX92fm6FsgFWpZmci6JaJ9FpDNwF3CJqh4szrZxKJp9jv5z9vokRxQnR8rjO5HTkP+eHGlyXMx1HHui8w1nuQnHntxdS9k4uRvNPtcq3Ed8J5R+Amp4vU/R7m9A7EROPLn7I74TftWd5bje31LY5+pARWc5A/ie404YxuMjwn/XLfEdrJxx3PqE/ZzD7HPUn7PnP4Aof3jdge+cH85dzroR+H47AqQCb+I7efs50Chg27uc7VYDf/Z6X9zeZ+BS4GvnH9hSoIfX+1JK+9sG39HTXmAH8HXAtgOdn8Ma4Cqv98XtfQbaA185n/FXwNVe70sp7vM84BdgmfN4Jwk+56D7XBqfs7VsMMaYJFOWx/iNMcaUgBV+Y4xJMlb4jTEmyVjhN8aYJGOF3xhjkowVfhN3AjoPrhSRN0Uk3cX3GuFcJIOIDC3ue4nPhyJykvP8BhFZJSKvl0JuA0Tk1IDnL4lI4xK+1vUiclW0OZnEYNM5TdwRkT2qWtlZfh3IV9UnI9hO8P2bLlHzMhFZB2Sr6vZibPMXoLOq3uQ8/xbfdSE/HhdXXn29k4qTTy5wq6rmFWe7EK+VDnyiqmXhSl7jMjviN/FuIXA6+PuTr3QeQ511mc4R9vP4LkyrJyL9ROQrJ+4xJ66ciEx01n0lIoWFeqKI9BaRG4BTgY9E5CMRuVpEnipMQkT+ISLBfvn8L/C2EzMO31XR74jITSIyXETGi8gc4N9OrgudXvpLRaR9wOvf7uS1XEQeFV+f/WzgdeevnzQRyRWRbCf+hH101u8RkYec11kiIqcAqK/L4zoRKQtdaI3bvL56zR72OP6B02sc32Xtb+PrRd4a31WKlYDK+K5CbomvN/1RoJ2zzanABnwtKsoDHwI9ne3nBrxHNefrRJy2B8A6IMNZroTvisoU5/lioFmQXNcDVQKeB77GcCAfSHOepwOpzvIZQJ6z/Gfn9Qt7r9dwvubi+wuEwOeh9tGJUZyrsoGRwN0B298F3OL152sP7x92xG/iUZqILAPy8BW4l4HzgRmquldV9wBvAYWtpder6hJnuQ2Qq6rb1De08jrQEV9flEYi8oyIdAPCdiVV1b34CurFIvJHfL8AvgoSWkNVfw/zUu+o6n5nOQV4UUS+wtdWo3C8vjMwQZ3e66r6a7jcwuwjwCHgXWc5H98vxkJb8f3SMEmuLLdlNolrv6q2CFxRxM1k9gaGBgtQ1Z0i0hzfjTuuA/6Or8dLOC8Bw4BvgQkhYgpE5A8a+rxCYG434eu90hzfMOuBgJyLc7It3M/isKoWvtYRjv0/ngrsP3ETk2zsiN+UFQuAniKSLiKVgF74xv+P9xlwgYhkiEg5oB/wsYhkAH9Q1enAPUCrINv+DlQpfKKqn+Frnfs/wOQQea3GN64fiarAZueXRH98t98DmAMMLJxRJCI1guVT1D5G8P5nAisjzNUkMCv8pkxQ1aX4xuM/x1f4XlLVL4PEbQbuBD7C6USqqm/ju7tRrjOENNGJOd544H0R+Shg3Rv4ZsPsDJHae0BOhLvxPHCliCzBV4T3Ojl/ALwD5Dn53erETwTGFZ7cjWAfi3Ievo6PJsnZdE5jwhCRd4GnVHV+iO/XAf6tqhfFNrPiEZGWwM2q2t/rXIz37IjfmCBEpJqIfIfvfEPQog/+o+8XCy/gimMZ+Ia4jLEjfmOMSTZ2xG+MMUnGCr8xxiQZK/zGGJNkrPAbY0ySscJvjDFJ5v8Dk/Y3Vpebnb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(X,color='red',alpha=0.2,edgecolor='black',bins=np.linspace(fmin,fmax,20))\n",
    "plt.xlabel(feature + ' ' + funits); plt.ylabel('Frequency')\n",
    "plt.axvline(x=average, ymin=0, ymax=1,color='black',label='Average')\n",
    "plt.axvline(x=median, ymin=0, ymax=1,color='black',label='Median',linestyle='--')\n",
    "plt.axvline(x=mode, ymin=0, ymax=1,color='black',label='Mode',linestyle='dashdot')\n",
    "plt.axvline(x=power_avg, ymin=0, ymax=1,color='black',label='Power',linestyle='dotted');\n",
    "plt.legend(loc='upper left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Dispersion\n",
    "\n",
    "##### Population Variance\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_{x} = \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity population variance is 0.0011 (fraction^2).\n"
     ]
    }
   ],
   "source": [
    "varp = stats.pvariance(X)\n",
    "print(feature + ' population variance is ' + str(round(varp,4)) + ' ' + f2units + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sample Variance\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_{x} = \\frac{1}{n-1}\\sum^n_{i=1}(x_i - \\overline{x})^2\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0011 (fraction^2).\n"
     ]
    }
   ],
   "source": [
    "var = stats.variance(X)\n",
    "print(feature + ' sample variance is ' + str(round(var,4)) + ' ' + f2units + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Population Standard Deviation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{x} = \\sqrt{ \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)^2 }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0329 (fraction^2).\n"
     ]
    }
   ],
   "source": [
    "stdp = stats.pstdev(X)\n",
    "print(feature + ' sample variance is ' + str(round(stdp,4)) + ' ' + f2units + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sample Standard Deviation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{x} = \\sqrt{ \\frac{1}{n-1}\\sum^n_{i=1}(x_i - \\mu)^2 }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0329 (fraction^2).\n"
     ]
    }
   ],
   "source": [
    "std = stats.stdev(X)\n",
    "print(feature + ' sample variance is ' + str(round(std,4)) + ' ' + f2units + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Range\n",
    "\n",
    "\\begin{equation}\n",
    "range_x = P100_x - P00_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity range is 0.17 (fraction).\n"
     ]
    }
   ],
   "source": [
    "range = np.max(X) - np.min(X)\n",
    "print(feature + ' range is ' + str(round(range,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Percentile\n",
    "\n",
    "\\begin{equation}\n",
    "P(p)_x = F^{-1}_{x}(p)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 13th percentile is 0.11 (fraction).\n"
     ]
    }
   ],
   "source": [
    "p_value = 13\n",
    "percentile = np.percentile(X,p_value)\n",
    "print(feature + ' ' + str(int(p_value)) + 'th percentile is ' + str(round(percentile,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Inter Quartile Range\n",
    "\n",
    "\\begin{equation}\n",
    "IQR = P(0.75)_x - P(0.25)_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity interquartile range is 0.04 (fraction).\n"
     ]
    }
   ],
   "source": [
    "iqr = scipy.stats.iqr(X)\n",
    "print(feature + ' interquartile range is ' + str(round(iqr,2)) + ' ' + funits + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tukey Test for Outliers\n",
    "\n",
    "Let's demonstrate the Tukey test for outliers based on the lower and upper fences.\n",
    "\n",
    "\\begin{equation}\n",
    "fence_{lower} = P_x(0.25) - 1.5 \\times [P_x(0.75) - P_x(0.25)]\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "fence_{upper} = P_x(0.75) + 1.5 \\times [P_x(0.75) - P_x(0.25)]\n",
    "\\end{equation}\n",
    "\n",
    "Then we declare samples values above the upper fence or below the lower fence as outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity outliers by Tukey test include [0.06726 0.05    0.06092].\n",
      "Porosity outlier indices by Tukey test are [110 152 198].\n"
     ]
    }
   ],
   "source": [
    "p25, p75 = np.percentile(X, [25, 75])\n",
    "lower_fence = p25 - iqr * 1.5\n",
    "upper_fence = p75 + iqr * 1.5\n",
    "outliers = X[np.where((X > upper_fence) | (X < lower_fence))[0]]\n",
    "print(feature + ' outliers by Tukey test include ' + str(outliers) + '.')\n",
    "outliers_indices = np.where((X > upper_fence) | (X < lower_fence))[0]\n",
    "print(feature + ' outlier indices by Tukey test are ' + str(outliers_indices) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's Visualize Outliers with a Box Plot (Box and Wisker Plot)\n",
    "\n",
    "The median is the orange line, P25 and P75 are the box and lower and upper fences are the wiskers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAD4CAYAAAD7CAEUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdPElEQVR4nO3deZhddZ3n8feHYhVkCZQzkIAJY3oggobmEplHpIkLBh5ZeqSFuACS6YgzIDaPPoDYg8bOCG60Ii5pwyKy2GIjGZSOCMFuFDAViJBAIyGyRFDiJMoSWRI+88c5hZebW1X3QJ1aUp/X85ynzvmd8/ud76lK7vf+zu8ssk1ERESnNhvuACIiYnRJ4oiIiEqSOCIiopIkjoiIqCSJIyIiKtl8uAMYCrvssosnTpw43GFERIwqS5Ys+b3t7tbyMZE4Jk6cSE9Pz3CHERExqkh6qF15TlVFREQlSRwREVFJEkdERFRSa+KQNEPSfZJWSDqzzfrTJd0j6S5JN0p6bVk+VdKtkpaX645tqnOJpF9LWlpOU+s8hoiIeKnaEoekLuBC4DBgCjBT0pSWze4EGrbfAFwNfK4sXwccb/v1wAzgHyXt2FTv47anltPSuo4hIiI2VmePYxqwwvZK288BVwFHNW9ge5HtdeXibcCEsvxXtu8v5x8FHgc2uiQsIiKGXp2JYzzwSNPyqrKsL7OA61sLJU0DtgQeaCqeW57COl/SVu0akzRbUo+kntWrV1ePPiIi2qozcahNWdtnuEt6P9AAPt9SvitwGfBB2y+UxWcBewEHAOOAM9q1aXue7YbtRnd3OisREYOlzhsAVwG7Ny1PAB5t3UjS24Gzgb+y/WxT+fbAD4FP2r6tt9z2Y+Xss5IuBj5WQ+wRL4vU7vvS4Mt7dGI41Zk4FgOTJU0CfgMcB7y3eQNJ+wHfBGbYfrypfEvgGuDbtr/XUmdX24+p+B96NLCsxmOIqKTqB7qkJIEYdWpLHLbXSzoFWAh0ARfZXi5pDtBjewHFqantgO+V39Qetn0k8B7gYGBnSSeWTZ5YXkF1uaRuilNhS4GT6zqGiIjYmMbCt51Go+E8qypGovQ4YiSTtMR2o7U8d45HREQlSRwREVFJEkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVJLEERERlSRxREREJUkcERFRSRJHRERUksQRERGVJHFEREQlSRwREVFJEkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVJLEERERldSaOCTNkHSfpBWSzmyz/nRJ90i6S9KNkl7btO4ESfeX0wlN5ftLurts8yuSVOcxRETES9WWOCR1ARcChwFTgJmSprRsdifQsP0G4Grgc2XdccA5wJuAacA5knYq63wdmA1MLqcZdR1DRERsrM4exzRghe2Vtp8DrgKOat7A9iLb68rF24AJ5fw7gRtsr7G9FrgBmCFpV2B727faNvBt4OgajyEiIlrUmTjGA480La8qy/oyC7h+gLrjy/kB25Q0W1KPpJ7Vq1dXDD0iIvpSZ+JoN/bgthtK7wcawOcHqNtxm7bn2W7YbnR3d3cQbkREdKLOxLEK2L1peQLwaOtGkt4OnA0cafvZAequ4s+ns/psMyIi6lNn4lgMTJY0SdKWwHHAguYNJO0HfJMiaTzetGohcKikncpB8UOBhbYfA56UdGB5NdXxwLU1HkNERLTYvK6Gba+XdApFEugCLrK9XNIcoMf2AopTU9sB3yuvqn3Y9pG210j6DEXyAZhje005/2HgEmAbijGR64mIiCGj4uKkTVuj0XBPT89whxGxEUmMhf+DMTpJWmK70VqeO8cjIqKSJI6IiKgkiSMiIipJ4oiIiEqSOCIiopIkjoiIqCSJIyIiKkniiIiISpI4IiKikiSOiIioJIkjIiIqSeKIiIhKkjgiIqKSJI6IiKgkiSMiIipJ4oiIiEqSOCIiopIkjoiIqCSJIyIiKqk1cUiaIek+SSskndlm/cGS7pC0XtIxTeXTJS1tmp6RdHS57hJJv25aN7XOY4iIiJfavK6GJXUBFwLvAFYBiyUtsH1P02YPAycCH2uua3sRMLVsZxywAvhx0yYft311XbFHRETfakscwDRghe2VAJKuAo4CXkwcth8s173QTzvHANfbXldfqBER0ak6T1WNBx5pWl5VllV1HHBlS9lcSXdJOl/SVu0qSZotqUdSz+rVq1/GbiMiop06E4falLlSA9KuwL7Awqbis4C9gAOAccAZ7eranme7YbvR3d1dZbcRAIwbNw5JtU5A7fsYN27cMP8mY1Mz4KkqSQ3gLcBuwJ+AZcBPbK8ZoOoqYPem5QnAoxXjew9wje3newtsP1bOPivpYlrGRyIGy9q1a7ErfdcZkXoTVMRg6bPHIelESXdQfMPfBrgPeBw4CLhB0qWS9uin7cXAZEmTJG1JccppQcX4ZtJymqrshaDif8PRFIksIiKGSH89jm2BN9v+U7uV5WWwkymujNqI7fWSTqE4zdQFXGR7uaQ5QI/tBZIOAK4BdgKOkPRp268v259I0WP5aUvTl0vqpjgVthQ4uaMjjYiIQaFNoSs+kEaj4Z6enuEOI0YZSZvMqapN4Thi6ElaYrvRWt7JGEc38LfAxObtbZ80mAFGRMTo0Ml9HNcC/w78BNhQbzgRETHSdZI4XmW77SWvEREx9nRyH8d1kg6vPZKIiBgVOkkcp1Ekj2ckPVlOT9QdWEREjEwDnqqy/eqhCCQiIkaHjh5yKOlI4OBy8Wbb19UXUkREjGQDnqqSdC7F6ap7yum0siwiIsagTnochwNTbb8AIOlS4E5goxczRUTEpq/Tp+Pu2DS/Qx2BRETE6NBJj+OzwJ2SFlE8H+pgigcfRkTEGNTJVVVXSrqZ4v0XAs6w/du6A4uIiJGpv8eq71X+/EtgV4r3azwC7FaWRUTEGNRfj+N0YDbwxTbrDLy1logiImJE6zNx2J5dzh5m+5nmdZK2rjWqiIgYsTq5qurnHZZFRMQY0GePQ9J/BsYD20jaj2JgHGB74FVDEFtERIxA/Y1xvBM4EZhAMc7RmzieAD5Rb1gRETFS9TfGcSlwqaR32/7+EMYUEREjWCdjHPtLevHOcUk7SfqHGmOKiIgRrJPEcZjtP/Qu2F5L8fyqAUmaIek+SSskbfRsK0kHS7pD0npJx7Ss2yBpaTktaCqfJOl2SfdL+q6kLTuJJSIiBkcniaNL0la9C5K2AbbqZ/ve7bqAC4HDgCnATElTWjZ7mGIc5Yo2TfzJ9tRyOrKp/DzgfNuTgbXArA6OISIiBkknieM7wI2SZkk6CbgBuLSDetOAFbZX2n4OuAo4qnkD2w/avgt4oZNgJYnixsOry6JLgaM7qRsREYOjk2dVfU7S3cDbKK6s+ozthR20PZ7iESW9VgFvqhDb1pJ6gPXAubZ/AOwM/MH2+qY2x7erLGk2xZ3v7LHHHhV2G1HwOdvDp0b/w6B9zvbDHUJsYjp6A6Dt64HrK7atNmWuUH8P249K2hO4qUxe7d513rZN2/OAeQCNRqPKfiMA0KefwB79/3Qk4U8NdxSxKenkDYAHSlos6SlJz5WD1u0+wFutAnZvWp4APNppYLYfLX+uBG4G9gN+D+woqTfhVWozIiJeuU7GOL4KzATuB7YB/gdwQQf1FgOTy6ugtgSOAxYMUAd48ZLfrcr5XYA3A/e4+Pq3COi9AusE4NpO2oyIiMHR0RsAba8AumxvsH0xML2DOuuBU4CFwL3AP9teLmmOpCMBJB0gaRXwN8A3JS0vq+8N9Ej6JUWiONf2PeW6M4DTJa2gGPOY3+nBRkTEK9fJGMe6ssewVNLngMeAbTtp3PaPgB+1lP3vpvnFFKebWuv9HNi3jzZXUlyxFRERw6CTHscHyu1OAZ6mGLd4d51BRUTEyNVvj6O8iW+u7fcDzwCfHpKoIiJixOq3x2F7A9Cdx3pERESvTsY4HgR+Vj4v6uneQttfqiuoiIgYuTpJHI+W02bAq+sNJyIiRrr+3gB4me0PUDzi48tDGFNERIxg/Y1x7C/ptcBJ5Q1545qnoQowIiJGlv5OVX0D+FdgT2AJL332lMvyiIgYY/rscdj+iu29gYts72l7UtOUpBERMUb1mTgkbQdg+8MDbRMREWNHf2Mc10r6Yvl61xcfMSJpz/KlTguBGfWHGBERI0mfYxy23ybpcOBDwJsl7UTxUqX7gB8CJ9j+7dCEGRERI0W/93G0e0hhRESMbR09Vj0iIqJXR6+OjRirpHZvQB5ddtppp+EOITYxSRwRfRiK941L2iTeax5jSyfvHP+CpNcPRTARETHydTLG8R/APEm3SzpZ0g51BxURESPXgInD9rdsvxk4HpgI3CXpCkkDvnc8IiI2PR1dVVW+CXCvcvo98EvgdElXDVBvhqT7JK2QdGab9QdLukPSeknHNJVPlXSrpOWS7pJ0bNO6SyT9WtLScpra4bFGRMQgGHBwXNKXgCOAm4D/Y/sX5arzJN3XT70u4ELgHcAqYLGkBbbvadrsYeBE4GMt1dcBx9u+X9JuwBJJC23/oVz/cdtXD3x4EREx2Dq5qmoZ8Enb69qsm9ZPvWnACtsrAcreyVHAi4nD9oPluheaK9r+VdP8o5IeB7qBPxAREcOqk1NV72tNGpJuBLD9x37qjQceaVpeVZZVImkasCXwQFPx3PIU1vmStuqj3mxJPZJ6Vq9eXXW3ERHRh/6ejrt1+cKmXVpe5DQR2K2DttvdOVXpgnVJuwKXAR+03dsrOYtirOUAYBxwRru6tufZbthudHd3V9ltRET0o79TVR8CPkqRJO5oKn+CYuxiIKuA3ZuWJ1C8u7wjkraneJjiJ23f1ltu+7Fy9llJF7Px+EhERNSov6fjfhn4sqRTbV/wMtpeDEyWNAn4DXAc8N5OKkraErgG+Lbt77Ws29X2YyqeBXE0xRhMREQMkT4Th6S32r4J+I2k/9663va/9New7fWSTgEWAl0UbxJcLmkO0GN7gaQDKBLETsARkj5t+/XAe4CDgZ0lnVg2eaLtpcDlkropToUtBU6ueMwREfEKqK/n5JQf4ueUp4Na2fZJ9YY2eBqNhnt6eoY7jIiN5FlVMZJJWmK70Vre36mqc8qfH6wzsIiIGF06ecjhaZK2V+Fb5Z3ehw5FcBERMfJ0ch/HSbafAA4FXgN8EDi31qgiImLE6iRx9N6PcThwse1f0v4ejYiIGAM6SRxLJP2YInEslPRq4IUB6kRExCaqk2dVzQKmAittr5O0M8XpqoiIGIMGTBy2X5A0AXhv+f7ln9r+v7VHFhERI1InV1WdC5xG8VTbe4CPSPps3YFFRMTI1MmpqsOBqb0PGZR0KXAnxcMGIyJijOnoDYDAjk3zeed4RMQY1kmP47PAnZIWUVyGezDpbUREjFn9Jo7yCbS3AAdSvP9CwBm2fzsEsUVExAjUb+KwbUk/sL0/sGCIYoqIiBGskzGO28rHn0dERHQ0xjEd+JCkh4CnKU5X2fYbao0sIiJGpE4Sx2G1RxEREaPGgKeqbD9EcTnuEeW0Y1kWERFjUEfv4wAup3ik+muA70g6te7AIiJiZOr0IYdvsv00gKTzgFuBC+oMLCIiRqZO38exoWl5A3kfR0TEmNVJ4rgYuF3SpyR9CrgNmN9J45JmSLpP0gpJZ7ZZf3D5Ktr1ko5pWXeCpPvL6YSm8v0l3V22+ZXyJsWIiBginQyOf4ni/RtrgLXAB23/40D1JHUBF1JclTUFmClpSstmDwMnAle01B0HnAO8CZgGnCNpp3L114HZwORymjFQLBERMXj6HOOQtDVwMvA64G7ga7bXV2h7GrDC9sqyvauAoygezQ6A7QfLda1vFHwncIPtNeX6G4AZkm4Gtrd9a1n+beBo4PoKcUVExCvQX4/jUqBBkTQOA75Qse3xwCNNy6vKsldSd3w5P2CbkmZL6pHUs3r16o6DjoiI/vV3VdUU2/sCSJoP/KJi2+3GHvwK63bcpu15wDyARqPR6X4jImIA/fU4nu+dqXiKqtcqYPem5QnAo6+w7qpy/uW0GRERg6C/xPFGSU+U05PAG3rnJT3RQduLgcmSJknaEjiOzp+wuxA4VNJO5aD4ocBC248BT0o6sLya6njg2g7bjIiIQdDnqSrbXa+kYdvrJZ1CkQS6gItsL5c0B+ixvaB86u41wE7AEZI+bfv1ttdI+gxF8gGY0ztQDnwYuATYhmJQPAPjERFDSPamf/q/0Wi4p6dnuMOI2IgkxsL/wRidJC2x3Wgt7/Sd4xEREUASR0REVJTEERERlSRxREREJUkcERFRSRJHRERUksQRERGVJHFEREQlSRwREVFJEkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVJLEERERlSRxREREJUkcERFRSRJHRERUksQRERGV1Jo4JM2QdJ+kFZLObLN+K0nfLdffLmliWf4+SUubphckTS3X3Vy22bvuNXUeQ0REvFRtiUNSF3AhcBgwBZgpaUrLZrOAtbZfB5wPnAdg+3LbU21PBT4APGh7aVO99/Wut/14XccQEREbq7PHMQ1YYXul7eeAq4CjWrY5Cri0nL8aeJsktWwzE7iyxjgjIqKCOhPHeOCRpuVVZVnbbWyvB/4I7NyyzbFsnDguLk9T/X2bRAOApNmSeiT1rF69+uUeQ0REtKgzcbT7QHeVbSS9CVhne1nT+vfZ3hd4Szl9oN3Obc+z3bDd6O7urhZ5RET0qc7EsQrYvWl5AvBoX9tI2hzYAVjTtP44Wnobtn9T/nwSuILilFjEiCCp0vRy6vTRyY4YMnUmjsXAZEmTJG1JkQQWtGyzADihnD8GuMm2ASRtBvwNxdgIZdnmknYp57cA3gUsI2KEsD0kU8Rw2ryuhm2vl3QKsBDoAi6yvVzSHKDH9gJgPnCZpBUUPY3jmpo4GFhle2VT2VbAwjJpdAE/Af6prmOIiIiNaSx8e2k0Gu7p6RnuMCIiRhVJS2w3Wstz53jEMLjyyivZZ5996OrqYp999uHKK3PFeYwetZ2qioj2rrzySs4++2zmz5/PQQcdxC233MKsWbMAmDlz5jBHFzGwnKqKGGL77LMPF1xwAdOnT3+xbNGiRZx66qksW5ZrPWLk6OtUVRJHxBDr6urimWeeYYsttnix7Pnnn2frrbdmw4YNwxhZxEtljCNihNh777255ZZbXlJ2yy23sPfeew9TRBHVJHFEDLGzzz6bWbNmsWjRIp5//nkWLVrErFmzOPvss4c7tIiOZHA8Yoj1DoCfeuqp3Hvvvey9997MnTs3A+MxaqTHERERlaTHETHEcjlujHa5qipiiOVy3BgtcjluEkeMELkcN0aLXI4bMULkctwY7ZI4IoZYLseN0S6D4xFDLJfjxmiXMY6IiGgrYxwRETEokjgiIqKSJI6IiKgkiSMiIipJ4oiIiEpqTRySZki6T9IKSWe2Wb+VpO+W62+XNLEsnyjpT5KWltM3mursL+nuss5XJKnOY4ioQ945HqNZbfdxSOoCLgTeAawCFktaYPueps1mAWttv07SccB5wLHlugdsT23T9NeB2cBtwI+AGcD1NR1GxKDLQw5jtKuzxzENWGF7pe3ngKuAo1q2OQq4tJy/Gnhbfz0ISbsC29u+1cUNKN8Gjh780CPqM3fuXObPn8/06dPZYostmD59OvPnz2fu3LnDHVpER+pMHOOBR5qWV5VlbbexvR74I7BzuW6SpDsl/VTSW5q2XzVAmwBImi2pR1LP6tWrX9mRRAyie++9l4MOOuglZQcddBD33nvvMEUUUU2diaNdz6H1NvW+tnkM2MP2fsDpwBWStu+wzaLQnme7YbvR3d1dIeyIeuUhhzHa1Zk4VgG7Ny1PAB7taxtJmwM7AGtsP2v7/wHYXgI8APxFuf2EAdqMGNHykMMY7ep8yOFiYLKkScBvgOOA97ZsswA4AbgVOAa4ybYldVMkkA2S9gQmAyttr5H0pKQDgduB44ELajyGiEGXhxzGaFdb4rC9XtIpwEKgC7jI9nJJc4Ae2wuA+cBlklYAayiSC8DBwBxJ64ENwMm215TrPgxcAmxDcTVVrqiKUWfmzJlJFDFq5em4ERHRVp6OGxERgyKJIyIiKkniiIiISpI4IiKikjExOC5pNfDQcMcR0cYuwO+HO4iIPrzW9kZ3UI+JxBExUknqaXfVSsRIllNVERFRSRJHRERUksQRMbzmDXcAEVVljCMiIipJjyMiIipJ4oiIiEqSOKJ2kjZIWippmaTvSXpVjfuaI+nt5fxHq+5L0oOS7i7jvVtS6+uOByPGByXt0s/6HSX9z8Heb1P7T/VRPkHStZLul/SApC9L2rKD9j7Rrn1Ju0m6enCijpEkYxxRO0lP2d6unL8cWGL7Sx3UE8W/0Rde5n4fBBq2O77BrrmOpP8K/Nj2a1/O/l9uXJImAtfZ3mcw99vU/ot/j6YyUbzj5uu2L5bURTFwv8b2x6u01679DuPqsr2har0YeulxxFD7d+B1AJJOL3shyyR9tCybKOleSV8D7gB2lzSz/Pa/TNJ55XZdki4py+6W9Hdl+SWSjpH0EWA3YJGkRZJmSTq/NwhJfytpoOS1PbC2qU67eA+QdJekrSVtK2m5pH0kHSLp3yRdI+keSd+QtNH/t3ZtAucC/6Xs9Xy+TZ0fSFpS7mt2U/lTkuZK+qWk2yT9p7J8kqRbJS2W9Jk+jvWtwDO2LwYoP8D/DjhJ0qsknSjpq037uq48xnOBbcpYL2+Jc6KkZeV8l6TPlzHcJelDZfkh5d/nCuDu8nf4w/IYlkk6doC/UQwH25ky1ToBT5U/NweupXgZ1/7A3cC2wHbAcmA/YCLwAnBgWWc34GGgu6x/E3B0Wf+Gpn3sWP68BDimnH8Q2KWc35biFcRblMs/B/ZtE+uDZVzLgHXAu8rytvGW6/4B+AJwIXBWWXYI8AywJ8WLzG5ojWuA38Gyfn6f48qf25Rx7lwuGziinP8c8MlyfgFwfDn/v3r/Hi1tfgQ4v035ncAbgBOBrzaVXwcc0vz3bfP3fvE4gNlN8WwF9ACTyt/T08Ckct27gX9qamuH4f73m2njKT2OGArbSFpK8WHxMMWbHw8CrrH9tO2ngH8B3lJu/5Dt28r5A4Cbba+2vR64nOINkSuBPSVdIGkG8ER/Adh+miLpvEvSXhQJ5O4+Np/u4jTRvsBXJW03QLxzgHcADYoP7F6/sL3Sxbf3K8s2mvXXZn8+IumXwG3A7hSvVgZ4juIDHWAJxQc3wJvL/QNc1kebokg8nZZXdShwfPnv4HZgZ/4c9y9s/7qcvxt4u6TzJL3F9h8HYd8xyOp853hErz/ZntpcUJ5T78vTzZu228D2WklvBN5J8S36PcBJA8TxLeATwH8AFw8UtO0HJP0OmNJXHKVxFD2GLYCtm+Jv/cBtXe6vzbYkHQK8HfhvttdJurncJ8DzLr+mU7xyufn/90Af/sspvu0372t7isT0APBGXnpqe2uqEXCq7YUt+ziEpr+37V9J2h84HPispB/bnlNxX1Gz9DhiuPwbcHR5/nxb4K8pxj9a3Q78laRdygHbmcBPVVyVtJnt7wN/D/xlm7pPAq/uXbB9O8UH4Xv58zfwPkl6DcXplIcGiHdeGcPlwHlNTUwrxxc2A44Fbunwd/CSuFvsAKwtk8ZewIEDHQfwM+C4cv59fWxzI/AqScdDMSYBfBG4xPY6itNrUyVtJml3YFpT3eclbTFADAuBD/duJ+kvymN+CUm7Aetsf4fi9F+7v2sMs/Q4YljYvkPSJcAvyqJv2b5TxRVFzds9JuksYBHFt9Yf2b627G1c3DTgfFab3cwDrpf0mO3pZdk/A1Ntr22zfa9FkjZQ9CDOtP074Hd9xHs8sN72FeWH7c8lvZVinOZWioHufSmSxDWd/A4AJP2sHFi+3i+9qulfgZMl3QXcR3G6aiCnAVdIOg34frsNbFvSXwNfk/T3FF8qf0TRQ4Mi+fyaP4//3NFUfR5wl6Q7bPeVmL5FcersjrK3uZpirKrVvsDnJb0APE8xHhYjTC7HjTFF0nUUg8A31ryfQ4CP2X5XnfuJGA45VRVjgoqb6n5FMd5Sa9KI2NSlxxEREZWkxxEREZUkcURERCVJHBERUUkSR0REVJLEERERlfx/DtQ56TacFU4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot(X); plt.ylabel(feature + ' ' + funits)\n",
    "plt.xticks([1], [feature + ' Boxplot and Outliers']); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Shape\n",
    "\n",
    "##### Pearson's Mode Skewness\n",
    "\n",
    "\\begin{equation}\n",
    "skew = \\frac{3 (\\overline{x} - P50_x)}{\\sigma_x}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity skew is -0.03.\n"
     ]
    }
   ],
   "source": [
    "skew = (average - median)/std\n",
    "print(feature + ' skew is ' + str(round(skew,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Population Skew, 3rd Central Moment\n",
    "\n",
    "\\begin{equation}\n",
    "\\gamma_{x} = \\frac{1}{n}\\sum^n_{i=1}(x_i - \\mu)^3\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 3rd cenral moment is -1.22e-05.\n"
     ]
    }
   ],
   "source": [
    "cm = scipy.stats.moment(X,moment=3)\n",
    "print(feature + ' 3rd cenral moment is ' + str(round(cm,7)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Quartile Skew Coefficient\n",
    "\n",
    "\\begin{equation}\n",
    "QS = \\frac{(P75_x - P50_x) - (P50_x - P25_x)}{(P75_x - P25_x)}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity quartile skew coefficient is 0.14.\n"
     ]
    }
   ],
   "source": [
    "qs = ((np.percentile(X,75)-np.percentile(X,50))\n",
    "          -(np.percentile(X,50)-np.percentile(X,25))) /((np.percentile(X,75))-np.percentile(X,25))\n",
    "print(feature + ' quartile skew coefficient is ' + str(round(qs,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the Nonparametric CDF\n",
    "\n",
    "Let's demonstrate plotting a nonparametric cumulative distribution function (CDF) in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df3ycZZnv8c+VH9NCQlsgaUt/REoLWVQsLUgPdAsFWYStWDntQdQidYsgW3SrqAePLiuse3aPu+zCul0pWrbAustWoxKEFRYlUEUeoVMkgobWFNo0jWlCt3QCbZLmOn/MM2U65Mck6WR+fd+vV16ZeeaZJ9fcDl697+e+r9vcHREREXlLSbYDEBERyTVKjiIiIimUHEVERFIoOYqIiKRQchQREUmh5CgiIpJCyVFE3sbM7jKzP892HCLZouQoBcnMXjGz35tZRdKxa82sIYthZZ2ZrTSznw11nrt/yt3/coR/46Nm9pyZxcxst5n9p5n9YfjaV82sx8z2hz8vm9k/mdlJSe9fbGZ94fsTPw+NJBaRkVJylEJWBvxZtoMYjJmVZTuGVGZWOor3fg64A/i/wBSgBvhnYGnSaf/h7scBJwBXAFOBzckJEmh198qkn8tHGpPISCg5SiH7W+DzZjapvxfN7Dwze9bM9oW/z0t6rcHM/tLMfh72cB4zs6rwtZPNzM3sOjNrDXtHNyW99xwz+4WZ/Xf42j+ZWSTpdTez1Wa2FdgaHrvTzHaa2etmttnMFiWd/1Uz+66Z/WsYS6OZnWZmXzKz9vB9lySdP9HM1od/e5eZfc3MSs3sdOAu4NywN/bf4fkbzOybZvaImXUBF4bHvpZ0zaVm9nwY3+/M7NJ+2nMicBuw2t2/7+5d7t7j7g+5+xdSzw9fexH4MLAHuCn1HJFsUXKUQvYc0AB8PvUFMzsBeBj4R+BE4O+Bh83sxKTTPgp8ApgMRPq5zoXAqcAlwM1mdnF4/BDwWaAKOBd4H/CnKe/9ELAAeGf4/FngTOK9qX8Dvmtm45POvxy4Hzge2AI8Svy/3+nEE9K6pHPvBXqBOcC8ML5r3f03wKeAX4S9seR/NHwU+CvgOOCIYVczOwe4D/gCMAk4H3iFtzsXGA/8oJ/XBuTuh4AHgUVDnSsyVpQcpdDdAnzazKpTji8Btrr7/e7e6+7/DvyWeBJK+Bd3f9nd3wQ2Ek9eyW4Ne0eNwL8AHwFw983u/kx43VeIJ64LUt771+7+Wnht3P1f3b0zfM/twDigNun8Te7+qLv3At8FqoG/cfce4AHgZDObZGZTgMuANWFs7cA/AFcN0U4PuvvP3b3P3Q+kvLYKuMfd/yt8fZe7/7afa5wIdIQxDlcr8X8YJEwLe96JnytHcE2REcu5+x0iR5O7/9rMfgTcDPwm6aVpwKspp79KvCeW0Jb0+A2gMuX8nSnvPQPAzE4j3hM9GziW+H9nmwd5L+Gw7LVhXA5MIN7zTPh90uM3iSehQ0nPCeObBpQDu80scX5J6t/rx2CvzwQeGeL9AJ1AlZmVjSBBTgdeS3re6u4zhnkNkaNGPUcpBn8BfJIjE18r8I6U82qAXcO47syU97aGj79JvBd6qrtPAP4PYCnvPbwdTnh/8X8DVwLHh8Od+/p5Tzp2AgeBKnefFP5McPd3pf7dgeIZ4Jqz0/jbvwAOEB8yTpuZlRDvsW8azvtEMknJUQqeu28D/gP4TNLhR4DTwmUHZWb2YeL3/340jEv/uZkda2bvIn5v8j/C48cBrwMxM/sD4IYhrnMc8XuEe4AyM7uFeM9x2Nx9N/AYcLuZTTCzEjObbWaJYd3fAzOSJwilYT3wCTN7X3i96eHnSv3b+4gPY681sw+FbVNuZpeZ2ddTzw9fOx34d+IzVv9+mB9XJGOUHKVY3AYcXvPo7p3AB4jPkOwEvgh8wN07hnHNJ4FtwE+Av3P3x8Ljnyc+wWU/8C3eSpoDeRT4T+Bl4sOzBxh6GHQwHyc+geglYC/wPSCxTOKnwItAm5ml9Vnd/ZfEk/8/EO/RPsnbe92Jc/8e+BzwFeLJfidwI/DDpNM+bGYx4L+BeuLtf5a7tyKSI0ybHYsMj5mdDGwHykc4+UREcpx6jiIiIikylhzN7J5wgfKvB3jdzOwfzWybmb1gZvMzFYuIiMhwZLLnuAF4WxWNJJcRX0B9KnAd8Rl+IjnP3V9xd9OQqkjhylhydPenOHLdUqqlwH0e9wwwKaW2ooiISFZkswjAdI6ckdcSHtudeqKZXUe8d8n48ePPqqmpGZMAC01fXx8lJbrNPFxqt5FT241Mptutr6/v8N8Y6u/09vZivb2UmnHIHS8ro6wsd+vHvPzyyx3unloRa9iy+Qn7W+Dc79RZd78buBugtrbWm5qaMhlXwWpoaGDx4sXZDiPvqN1GTm03MplstyAIWLd6NUuIFxe+fu1aFixYcNTOzzYzS618NSLZTI4tHFlhZAZvVRgREZEMaIxGuSAW473l5cR6emiMRgdNdgsWLIC1a2mMRrl+/vycToxHUzaTYz1wo5k9QHx3gn1hdQ8REckQi0S459VXeePQIR4oLeXjkaGLJS1YsKBokmJCxpKjmf07sJh4IeIW4vUtywHc/S7i5bv+mHiFkTeIV+AQEZEMam5qYq47XaWlzHWnWbep+pWx5OjuHxnidQdWZ+rvi4jI2+1qbaWjr4+FkQg/d+f0bAeUozSNTESkSARBQOfmzewvKeH73d28Xl3N0mXLsh1WTsrd+bgiInJUNUajLHdn0vTp1MdizL7yyqK7l5gu9RxFRIqERSJs2LmTPW1tNMdizK6tzXZIOUs9RxGRPBQEAY3RKGcMY3mFd3dzSU0NlJVxSW8v3t2d4Sjzl3qOIiJ5JggCbl+5kt6//mtuX7mSIAjSet8Z8+eztaKC48eNY2tFBWfM134PA1HPUUQkz9TX1XHOjh0sLS8n1tlJfV1dWr3HYl3QPxJKjiIiOWqgoVMHfg7McufnMKzlGMW4oH8kNKwqIpKDEjVNj1+/nnWrVx8xdLp02TJ6a2r4ycSJ9NbUaDlGBqjnKCKSZf31EBujUZYAy046CXbvPqIG6oIFC/jKhg00RqNco+HRjFByFBHJouRdL9atXw9r1wLxyTPr1q+H3bvju2GkTJ7R8GhmKTmKiGRRf7tkzDn9dE2eyTIlRxGRLNre0sJT27dzwIwHSktZkbRLhnqH2aMJOSIiWRIEAZvuv59FZrSasaiqSgvzc4SSo4hIljRGoyw241UzprmzBbQwP0coOYqIZIlFImzq6OAd7tzZ18d7VqzQMGqO0D1HEZEsaW5q4pzKSqZUVnJ1JEL1jBnZDklC6jmKiGRBEAS88NBDvByL0dPWpiHVHKOeo4hIFtTX1bGkq4vZc+bwg717mXP55RpSzSFKjiIiGdRf9ZtEr/Glzk4WdnbSUVPDNSoBl1M0rCoikiED1UdtjEZZWVHBzfPm0Xziieo15iAlRxGRDDlc/ebgQS6IxWiMRoH4vcWHgZYDBzgwebIKh+cgDauKiGSIRSJs2LmTHuA7cLj6jUrD5T4lRxGRDPHubhZVVdHZ08Oi8vIjqt+oNFxu07CqiEiGbG9p4Ym2Nia99hqbOjqwpLqpktuUHEVEMkB1U/ObkqOISAaobmp+U3IUEckA1U3Nb5qQIyJylAVBwM8ef5xFVVXUTJjA1b29qpuaZ9RzFBE5ihIL/097/nk2dXRQ3tPD1ooKDanmGfUcRUSOovq6Os5pb+fDNTUANJx5Jp9as0ZDqnlGyVFEZJQS9VMtEjlcMzXW2cmWmhpuUmLMS0qOIiKjkBhGXQLc2d7OqmOP5Q/mzePeHTtUMzWP6Z6jiMgoJIZRLxo/nssqKvjem2+qZmoBUM9RRGSEkreeSgyjfuiWW9jb3a2aqXlOyVFEZASCIOCuO+5gufsRw6irVq3KdmhyFCg5iogM0+HlGl1dbNi5k5WgYdQCo+QoIjJMWq5R+JQcRUSGkFiqkVjIr+UahU/JUURkEMlLNdatX89JF13EyooKZmi5RkHTUg4RkUEkL9VYAjjwMGi5RoFTz1FEZACpSzV+WVPDTcuWwbJlNEajWq5RwJQcRUQG0BiNstydSVOnUh+LHTGEqqRY2DSsKiIyAItE2LBzJ3va2miOxZhdW5vtkGSMqOcoIjIA7+5mUVUVnT09LCovx7u7sx2SjBH1HEVEBrC9pYUn2tqY9NprbOrowCKRbIckY0TJUUSkH0EQsOn++1lkRqsZi6qq1HMsIkqOIiL9qK+r49yDB2kpK2NmSQlb4HARACl8So4iIikSSzhejsWY1NvLt447jg/dcotmqBYRTcgREUmRuoTjfR//uHbbKDLqOYqIpNASDlHPUUQkRXNTE+dUVtJVWcklkYgm4hQh9RxFRJKsX7+eJ++7j9/s20dPW5sm4hQp9RxFREJBEHDvrbeycv9+xpWV0VBRoV03ipR6jiIiocZolMVmPNHXxxu9vbSPH69dN4qUkqOISMgiETZ1dPAOd+7s6+M9K1ao11ikNKwqIhJK1FKt6OlheXk5NTNmZDskyRL1HEVEQqqlKglKjiIiqJaqHEnJUUQE1VKVIyk5ikjRUy1VSZXR5Ghml5pZk5ltM7Ob+3m9xsyeMLMtZvaCmf1xJuMREelPopbqyqlT6Zo4UbVUJXPJ0cxKgbXAZcA7gY+Y2TtTTvsKsNHd5wFXAf+cqXhERAaiWqqSKpNLOc4Btrl7M4CZPQAsBV5KOseBCeHjiUBrBuMREelXYglHZ08Pi8rLNRFHMHfPzIXNlgOXuvu14fOrgQXufmPSOScBjwHHAxXAxe6+uZ9rXQdcB1BdXX3Wxo0bMxJzoYvFYlRWVmY7jLyjdhu5fGm71tZW9re1cSLwmhknzJxJVVVV1uLJl3bLRRdeeOFmdz97tNfJZM/R+jmWmok/Amxw99vN7FzgfjN7t7v3HfEm97uBuwFqa2t98eLFmYi34DU0NKC2Gz6128jlQ9sFQcBdN9/MH7a302qGT55Myc03s3z58qzFlA/tVugyOSGnBZiZ9HwGbx82XQVsBHD3XwDjgez9c01Eio6WcEh/MpkcnwVONbNZZhYhPuGmPuWcHcD7AMzsdOLJcU8GYxIROUxLOGQgGRtWdfdeM7sReBQoBe5x9xfN7DbgOXevB24CvmVmnyU+5LrSM3UTVEQkRWIJx6SpU6mPxbSEQw7LaOFxd38EeCTl2C1Jj18CFmYyBhGRgSSWcHwMaAbO0xIOCWlXDhEpWlrCIQNR+TgRKVoWifBkezsTOjp4sr1du3DIYUqOIlK0mpuamOtOV2kpc91pbmrKdkiSI5QcRaRo7Wptpbmvj5mlpewsLX3bQmwpXkqOIlKUgiCgc/Nm9peU8P3ubl6vrmbpsmXZDktyhCbkiEhROryMY/p06mMxZl95pdY3ymHqOYpIUdJOHDIY9RxFpChpGYcMRslRRIpKEAQ0RqNsb2nh6Y4OPgb8GFihZRySRMlRRIpGEASsW72aJcB97e0srqqCCRO4pLdXPUc5gu45ikjRaIxGuSAW470HD7LYjC3A8ePGsbWiQjtxyBHUcxSRopGYhNMDbALOW7OGvTNmcP38+ZqpKkdQchSRohAEAT97/HEWJQ2lVs+YwbU33JDt0CQHaVhVRApe4l7jac8/z6aODsp7ejSUKoNSz1FECloQBNx1xx0sjsW45tRTAWg480w+tWaNhlJlQEqOIlKw1q9fzw9vu415wIaODgC2VlYqMcqQlBxFpCAFQcC9t97Kxzs7Obu8HKqqaJg3T4lR0qLkKCIFqTEaZbEZT/T1cejAAbZMnMhXlBglTZqQIyIFySIRNnV08A537uzr4z0rVigxStrUcxSRgpAoC3dGuGbRu7u5pKaG6rIyrg6XbYikS8lRRPJeEATcvnIlF3V1cXtFBTdt2MAZ8+ezrqKC04Cnx43jIi3bkGFQchSRvFdfV8c5O3awtLycWGcn9XV1/NXXvw5r19IYjaoCjgybkqOI5D0Hnjp0iOl9fTzlzrvD4wsWLFBSlBHRhBwRyXuza2vZBTx26BC7wucio6Geo4jkPe/uZsmUKVT09DBTGxfLUaCeo4jkvcSyjRP37WNTRwemjYtllNRzFJG819zUxDmVlXRVVnJJJKKeo4yakqOI5K0gCHiwro5g40YmxGJUx2JsqanhJi3bkFFSchSRvJTYhmp8ezuXdXYyd84cfrB3L3Muv1wzVGXUlBxFJC81RqNcEIsxqbyc9YcOUbl3LwcmT2bpsmXZDk0KgCbkiEheskiEDTt3sqetjV1AcP75XL92rXqNclSo5ygieSl5Es7ySITqhQuVGOWoUc9RRPJOEAS88NBDvByL0dPWxhbgDE3CkaNIPUcRyTv1dXUs6epitibhSIao5ygieSUIAp7ZuJEftbWxpamJjooKTcKRo049RxHJK/V1dVy2Zw+nRyLc193NsWedpV6jHHXqOYpIXtnV2kpDTw8xoLu8nGnTpmU7JClASo4ikjeCIKBz82b2l5Tw/e5uXq+u1pCqZISGVUUkbzRGoyx3Z9L06dTHYsy+8koNqUpGqOcoInkjeeF/cyymfRslY9RzFJG84d3dLKqqorOnh0Xat1EySD1HEckbFonwZHs7Ezo6eLK9Xfs2SsYoOYpI3mhuamKuO12lpcx1p7mpKdshSYFSchSRvLGrtZXmvj5mlpays7QUz3ZAUrCGTI5m9pyZrTaz48ciIBGR/mgZh4yldHqOVwHTgGfN7AEze7+ZWYbjEhE5QqKe6ldqazlx6lQWaBmHZNCQydHdt7n7l4HTgH8D7gF2mNmtZnZCpgMUEUnswvFoZye/2rZN9VQl49K652hm7wFuB/4WqAOWA68DP81caCIicYle45/OmUPziSdqFw7JuCHXOZrZZuC/gfXAze5+MHwpMLOFmQxORIpbEAQ8WFdHsHEjL3V2srCzk46aGq5Rr1EyLJ0iAP/L3ZuTD5jZLHff7u7/M0NxiUiRC4KAdatXM769ncs6O5mrvRtlDKUzrPq9NI+JiBw1jdEoF8RivL+8nKcOHeJ3e/dyYPJk3WuUMTFgz9HM/gB4FzDRzJJ7iBOA8ZkOTESKW6KO6seAXUBw/vl86tOfVq9RxsRgw6q1wAeAScDlScf3A5/MZFAiIs1NTZxTWUlXZSXLIxGqFy5UYpQxM2BydPcHgQfN7Fx3/8UYxiQiRW79+vU8ed99nLBvH9WxGFtqarhp/vxshyVFZLBh1S+6+9eBj5rZR1Jfd/fPZDQyESlKQRBw7623snL/fsaVldFQUaFJODLmBhtW/U34+7mxCEREBOITcRab8URfH+e70z5+PNdpEo6MscGGVR8Kf987duGISLGzSIRNHR2c686dfX1csWKFeo0y5gYbVn0IBi567+4fzEhEIlLUEhNxplRWcnUkQvWMGdkOSYrQYMOqfzdmUYiI8FYN1bJYTBNxJKsGG1Z9ciwDERFpjEZZ7s6kqVOpj8U0EUeyZsAKOWa2MfzdaGYvJP00mtkL6VzczC41syYz22ZmNw9wzpVm9pKZvWhm/zayjyEihcAiEe559VVad+5k2+uvM7u2NtshSZEabFj1z8LfHxjJhc2sFFgL/BHQQnw/yHp3fynpnFOBLwEL3X2vmU0eyd8SkcLQ3NTEXHe6SkuZ605zU1O2Q5IiNWDP0d13h79fBQ4Cc4H3AAfDY0M5B9jm7s3u3g08ACxNOeeTwFp33xv+rfbhfwQRKRS7Wltp7utjZmkpO0tLB54RKJJh5j7418/MrgVuIb53owEXALe5+z1DvG85cKm7Xxs+vxpY4O43Jp3zQ+BlYCFQCnzV3X/cz7WuA64DqK6uPmvjxo1pf0B5SywWo7KyMtth5B2128gNp+26urrY/corHDp4kHKgt7yc6aecQkVFRWaDzEH6zo3chRdeuNndzx7tddLZsuoLwDx37wQwsxOBp4FBkyPxRJoqNROXAacCi4EZwCYze7e7//cRb3K/G7gboLa21hcvXpxG2JKqoaEBtd3wqd1Gbjht9+1vfpPyb3yD6p4e6mMxplxzDdd9+tOZDTBH6TuXfelsWdVCvNh4wn5gZ5rvm5n0fAbQ2s85D7p7j7tvB5qIJ0sRKTKJXTj2tLXRHItpMo5k1WBFAD4XPtwFBGb2IPGe31Lgl2lc+1ngVDObFV7jKuCjKef8EPgIsMHMqoDTgGZEpOgk78JxSSSCd3dnOyQpYoP1HI8Lf35HPIklhkQfBHYPdWF37wVuBB4lXqd1o7u/aGa3mVmius6jQKeZvQQ8AXwhMXwrIsUjCAKe2biRF/fu5eCuXWwBztDif8miwYoA3Drai7v7I8AjKcduSXrswOfCHxEpUvV1dVy2Zw+nRyLc193NsWedpcX/klVDTsgxs2rgi8C7gPGJ4+5+UQbjEpEi4sBThw4x3Yw3S0qYPW1atkOSIpfOhJzvAL8FZgG3Aq8Qv58oInJUzK6tZRfw2KFD7Aqfi2RTOks5TnT39Wb2Z2G91SfNTHVXReSoaW5q4uKJE5lWWclpmowjOSCdnmNP+Hu3mS0xs3nEl2WIiIxaYieOl2MxetraNBlHckI6PcevmdlE4CbgG8AE4LMZjUpEikZ9XR1LurqYPWcOP9i7VztxSE4YMjm6+4/Ch/uACzMbjogUk8QSjl+1tXF+ezsds2ZxzbJl2Q5LZOhhVTM7xcweMrMOM2s3swfN7JSxCE5ECltiCccNkQib+/q0hENyRjrDqv9GfOupK8LnVwH/DugbLCIjEgQBjdEou1pbeQmYVVpKd3m5lnBIzkgnOZq735/0/F/N7MYBzxYRGUQQBKxbvZolQGdXF29UV/OTvj56KypYqiFVyRGD1VY9IXz4hJndTHw/Rgc+DDw8BrGJSAGqr6vjnPZ2LqqpAWDz5ZdzyqxZXDN/voZUJWcM1nPcTDwZJraeuj7pNQf+MlNBiUhhSizbeKmzk1hnJ7+sqeGmZcuUFCXnDFZbddZYBiIihStxj3H79u2srKhgxrx53Ltjh5ZtSM5Kp7ZqOXADcH54qAFY5+49A75JRCQUBAG3r1zJRV1dPFNSwtZjjuHDwIHJk3WPUXJWOhNyvgmUA/8cPr86PHZtpoISkcJRX1fHOTt2sLS8nFhPD7++4gr2LlzI9brHKDksneT4Xnefm/T8p2b2q0wFJCKF5fCOG319POXOu6dN49obbsh2WCKDSqe26iEzm514EhYAOJS5kESkkGjHDclH6fQcv0B8OUcz8Zmr7wA+kdGoRKRgeHc3S6ZMoaKnh5nl5dpxQ/LCoMnRzEqAN4FTgVriyfG37n5wDGITkQKwvaWFp9raWGHGj0tLWRGJZDskkSENmhzdvc/Mbnf3c4EXxigmESkQXV1dbLr/fhaZ0WrGoqoq9RwlL6Rzz/ExM1tmZjb0qSIib3nzjTdYbMarZkxz116NkjfSuef4OaAC6DWzA8SHVt3dJ2Q0MhHJf2Zs6ujgXHfu7OvjihUrtHxD8kI6+zkeNxaBiEhhCYKA2P79LKqqombCBK7u7aV6xoxshyWSlgGHVc1sspndYWY/MrP/a2bqKYpIWhI7b4x/4w02dXRQ3tPD1ooKDalK3hjsnuN9QBfwDeA44B/HJCIRyXuN0SgXxGKcEE7CaZg3j+vXrtWQquSNwZLjVHf/srs/6u6fBt4zVkGJSH6zSIR7Xn2VfQcO8GR7O3948cVKjJJXBrvnaGZ2PG9tWVWa/NzdX8t0cCKSn5qbmpjrTh8w153mpqZshyQyLIMlx4nE93RMXsIRDX87cEqmghKR/ObAK8Di8Pe7sxiLyEgMtp/jyWMYh4gUkNm1tTwKvA7sApaqnqrkmXSKAIiIDEuinur4sjKWTJmiqjiSd5QcReSos0iETR0dlB06xKaODkz1VCXPKDmKyFHX3NTEOZWV9JWXc0lNjXqOknfSSo5m9odm9onwcbWZzcpsWCKSr4Ig4JmNG3lx7168u1v1VCUvDVk+zsz+Ajib+JZV/wKUA/8KLMxsaCKSj+rr6rhszx5Oj0T4PXDsWWdpjaPknXR6jlcAHyReLQd3byVeMUdE5G12tbbS0NNDDHAzpk2blu2QRIYtneTY7e5OfOkSZlaR2ZBEJF8FQUDn5s3sLynh+93dHCorY+myZdkOS2TY0tmyaqOZrQMmmdkngT8BvpXZsEQkHzVGoyx3Z9L06dTHYlSccIKGVCUvpbNl1d+Z2R8RX89bC9zi7v+V8chEJO9YJMKGnTv5GNAMXDZuXLZDEhmRdCbkfBb4rhKiiAwlsYSjq7KSSyIRcM92SCIjks49xwnAo2a2ycxWm9mUTAclIvln/fr1PHnfffxm3z562trYAhxz7LHZDktkRIZMju5+q7u/C1gNTAOeNLPHMx6ZiOSNIAi499ZbWbl/P/+rrIzfVFYy5/LLqajQ/D3JT8OpkNMOtAGdwOTMhCMi+agxGmWxGU/09fFGby/t48drlqrktSGTo5ndYGYNwE+AKuCT7q6Nj0XksEQt1Xe4c2dfH+9ZsUKzVCWvpbOU4x3AGnd/PtPBiEh+SkzEmVJZydWRCNUzZmQ7JJFRGbDnaGYTwodfB3aY2QnJP2MTnojkuuRaqgd37VItVSkIg/Uc/w34ALCZeHUcS3rNgVMyGJeI5InkWqr3dXerlqoUhAGTo7t/IPytHThEZEAOPHXoENPNeLOkhNmqpSoFIJ0JOT9J55iIFKfZtbXsAh47dIhd4XORfDdgz9HMxgPHAlVmdjxvDatOIL7eUUSKXBAE/Ozxx1kyZQo1EyZwWm+vNjaWgjDYPcfrgTXEE+Fm3kqOrwNrMxyXiOS4IAhYt3o1p3V18VhHByuPOYatlZVcpMk4UgAGu+d4J3CnmX3a3b8xhjGJSB6or6vjnPZ2PlxTA0DDmWfyqTVrNBlHCkI6u3J8w8zeDbwTGJ90/L5MBiYiuSkIAh6sqyPYuJGXOjuJdXaypaaGm5QYpYCksyvHXwCLiSfHR4DLgJ8BSo4iRSYxlDq+vZ3LOjuZO2cOP9i7lzmXX67EKAUlndqqy4H3AW3u/glgLqBN2kSKUGM0ygWxGO8vL+epQ4f43d69HJg8WXVUpeCkk0NUtSEAABomSURBVBzfdPc+oDesmtOOCgCIFKXEZsZ72trYBQTnn8/1a9eq1ygFJ53aqs+Z2STgW8RnrcaAX2Y0KhHJSd7dzaKqKjp7elhSXk7NwoVKjFKQ0pmQ86fhw7vM7MfABHd/IbNhiUguSuy+8THgx8CKSCTbIYlkxGBFAAZcrGRm8909mpmQRCQXJRb8L6qqggkTuEQL/qWADdZzvH2Q1xy46CjHIiI5KggCbl+5klmvvcaTr7/On2jBvxS4wYoAXDiWgYhI7qqvq+OcHTv4WHk5f+3O92bO5Ctf+5ruN0rBSmed48f7O64iACLF4/DOG319vAKcMW+eEqMUtHRmq7436fF44mseo6gIgEjRmF1by6OEO2+UlLBUO29IgUtntuqnk5+b2UTg/nQubmaXAncCpcC33f1vBjhvOfBd4L3u/lw61xaRsdPc1MTFEycyrbKS0yIRTcSRgpdOEYBUbwCnDnWSmZUS373jMuKl5z5iZu/s57zjgM8AwQhiEZEMC4KAFx56iJdjMXra2tgCnKGJOFLg0rnn+BDxWw4QT6bvBDamce1zgG3u3hxe5wFgKfBSynl/CXwd+HyaMYvIGGqMRlnuzqSpU6mPxVRHVYpCOvcc/y7pcS/wqru3pPG+6cDOpOctwBH/RZnZPGCmu//IzAZMjmZ2HXAdQHV1NQ0NDWn8eUkVi8XUdiNQ7O02qbqazuuuo5v4v3hPrKlJuz2Kve1GSu2Wfencc3wSIKyrWhY+PsHdXxvirdbPMT/8olkJ8A/AyjRiuBu4G6C2ttYXL1481FukHw0NDajthq/Y2+3LX/wi+//xH6kxY4c7x33mM/zV17+e1nuLve1GSu2WfUPeczSz68zs98ALwHPE66umM2mmBZiZ9HwG0Jr0/Djg3UCDmb0C/A+g3szOTi90ERkLu1pbae7rY2ZpKTtLS9/6F65IAUtnWPULwLvcvWOY134WONXMZgG7gKuAjyZedPd9QFXiuZk1AJ/XbFWR3BEEAZ2bN7O/pITvd3fz+rRp2p5KikI6yfF3xGeoDou795rZjcCjxJdy3OPuL5rZbcBz7l4/3GuKyNg6PBln+nTqYzFmX3mlJuNIUUgnOX4JeNrMAuBg4qC7f2aoN7r7I8AjKcduGeDcxWnEIiJjKLF/48eAZuA8Lf6XIpFOclwH/BRoBPoyG46I5JLk/RsXlZdr8b8UjXSSY6+7fy7jkYhIzgiCgMZolO0tLTyt/RulCKWTHJ8I1xk+xJHDqkMt5RCRPBQEAetWr2YJcF97O4u1f6MUoXSSY2KG6ZeSjjlwytEPR0SyrTEa5YJYjPeWl7PYjC3A3HHjeHrcOO3fKEUjnSIAs8YiEBHJDYlJOD3AJuC8NWvYO2MG18+fr5mqUjS0n6OIHCF1Ek7NjBlce8MN2Q5LZExpP0cROcL2lhaeamtjhRk/Li3VJBwpShndz1FE8ksQBGy6/34WmdFqxqKqKk3CkaKUsf0cRST/1NfVce7Bg7SUlTGzpER7N0rRyuR+jiKSRxKbGpfFYsw8dIhvTZrEdbfcokk4UpQyuZ+jiOSR+ro6lnR1MXvOHH6wdy/vu+oqVq1ale2wRLJiwORoZnOAKYn9HJOOLzKzce7+u4xHJyJjItFrfKmzk4WdnXTU1HCNdt+QIjbYPcc7gP39HH8zfE1ECkRi942VU6fym8pK5lx+uYZTpagNlhxPdvcXUg+G+y2enLGIRGTMWSTCPa++SuvOnWx7/XVma/cNKXKDJcfxg7x2zNEORESyp7mpibnudJWWMted5qambIckklWDJcdnzeyTqQfNbBWwOXMhichYc+AVYHr42wc7WaQIDDZbdQ3wAzP7GG8lw7OBCHBFpgMTkbEzu7aWR4HHDh1iV0kJSzWsKkVuwOTo7r8HzjOzC4F3h4cfdvefjklkIjJmmpuauHjiRKZVVnJaJKKqOFL00ikf9wTwxBjEIiJZkLz4vzoWY0tNDTepKo4UuXSKAIhIAUss45g0dSr1sZiWcYgwstqqIlJAtIxD5O2UHEWKnJZxiLydkqNIkdMyDpG3U3IUKXKza2vZRbiMI3wuUuw0IUekyHl3N0umTKGip4eZ5eVaxiGCeo4iRc8iETZ1dHDivn1s6ujAIpFshySSdeo5ihQ57+5mUVUVnT09LFLPUQRQz1Gk6G1vaeGJtjYmvfaaeo4iISVHkSIWBAGb7r+fRWa0mrGoqko9RxGUHEWKWn1dHecePEhLWRkzS0rYApyh0nEiSo4ixSpRU/XlWIxJvb1867jj+NAtt6h0nAiakCNStOrr6ljS1cXsOXP4wd69vO+qq1i1alW2wxLJCeo5ihShRK/x0c5OfrVtGx0VFSxdtizbYYnkDPUcRYqQduIQGZx6jiJFSDtxiAxOyVGkCGknDpHBKTmKFJkgCPjV88/zCjArEmFnaal24hBJoXuOIkUkCALWrV7Nwq4uvgc8Mn48vSecoMk4IimUHEWKSGM0ygWxGBeWl/PmlCnsPO88vrJmjSbjiKRQchQpIttbWnhq+3YOmLGptJQVF1+sxCjSD91zFCkSqqMqkj4lR5EioTqqIulTchQpAqqjKjI8uucoUgRSK+K87+MfVx1VkUGo5yhSBCwSYcPOnexpa6M5FlNFHJEhqOcoUgS8u5tFVVV09vSwqLxcE3FEhqCeo0gR2N7SwhNtbUx67TU2dXRgkUi2QxLJaUqOIgVOSzhEhk/JUaTANUajLDbjVTOmuWsJh0galBxFCpxFImzq6OAd7tzZ18d7VqzQEg6RIWhCjkiBS0zGqejpYXl5OTUzZmQ7JJGcp56jSIFL9BxP3LdPk3FE0qSeo0gBC4KAnz3+OIuqqmDCBC7p7dVkHJE0qOcoUqDWr1/P1668kplPP82mjg7Ke3rYWlGhyTgiaVDPUaQABUHAvbfeysc7Ozm7vByqqmiYN49Pae9GkbQoOYoUoMQOHE+ZQU8PW0CbGosMg4ZVRQqMduAQGT0lR5ECEgQBd91xB8vduXnePHonT9YOHCIjoGFVkQIRBAHrVq/mtK4uNuzcyUrgwOTJLF22LNuhieQdJUeRPBcEAY3RKNu3b2cJsGzOHAAazjxTE3BERiijydHMLgXuBEqBb7v736S8/jngWqAX2AP8ibu/msmYRArJ+vXr+eFtt7H8mGN4wYytALt3s7WiQolRZBQylhzNrBRYC/wR0AI8a2b17v5S0mlbgLPd/Q0zuwH4OvDhTMUkUkiSl2vMLS+HqVNp+uAH2TtrFtfPn6/EKDIKmew5ngNsc/dmADN7AFgKHE6O7v5E0vnPACsyGI9IwUhMvFlsxtPl5dDTw4NvvslXli1TUhQ5CszdM3Nhs+XApe5+bfj8amCBu984wPn/BLS5+9f6ee064DqA6urqszZu3JiRmAtdLBajsrIy22HknVxrt66uLjp27GBcXx+vd3dTWVbGfuDEk06iqqoq2+EdIdfaLl+o3Ubuwgsv3OzuZ4/2OpnsOVo/x/rNxGa2AjgbuKC/1939buBugNraWl+8ePFRCrG4NDQ0oLYbvlxqtyAIuHfdOhZv2cJVp57KX2/bxi9yeOJNLrVdPlG7ZV8mk2MLMDPp+QygNfUkM7sY+DJwgbsfzGA8InktdakGwNbKypxNjCL5LJPJ8VngVDObBewCrgI+mnyCmc0D1hEffm3PYCwiea8xGuWCWIwLy8t5U7VSRTIqY8nR3XvN7EbgUeJLOe5x9xfN7DbgOXevB/4WqAS+a2YAO9z9g5mKSSSfbW9p4ant2zlgxqbSUlZcfLESo0iGZHSdo7s/AjyScuyWpMcXZ/LvixSKIAjYdP/9LDKj1YxFVVXal1Ekg1RbVSTHJS/baBs/npklJWwB7csokkEqHyeSw5In4TzW0cGiqiruc+ca7bIhklFKjiI5aKB6qS+feSZ/q0k4Ihmn5CiSY4Ig4PaVK7moq4tnSkrYeswxqpcqMsaUHEVyTH1dHefs2MHS8nJiPT38+oor2LtwoeqliowhJUeRHLOrtZVf9fQw04yfA6dPm8a1N9yQ7bBEiopmq4rkkCAI6Ny8mf0lJXy/u5vXq6u1WbFIFqjnKJIDkifgrKyoYMbZZ3Pvjh3MvvJKDaWKZIGSo0iWJZZrLAEe7OpiK/FNTQ9Mnqxeo0iWKDmKZMnblmucdBLs3s3miy7ShsUiWabkKJIFAy3XeBi4XhsWi2SdkqNIFmi5hkhuU3IUyQIt1xDJbVrKITLGtFxDJPep5yiSYYmJN2eEQ6aN0aiWa4jkOPUcRTIosUzj+PXrWbd6NUEQcMb8+TwMtBw4oOUaIjlKPUeRDGqMRrkgFuO94cSbxmg0fm9x7Voao1FNwBHJUUqOIhkSBAGbfv5zXnn1VXpKSvgOsCISAWDBggVKiiI5TMlRJAMSw6kV7e3Mdadr6lQuiUTw7u5shyYiadA9R5EMSAynvr+8nFeAcT09bK2o4Iz587MdmoikQclRJAMsEmHDzp3saWtjFxCcfz7Xr12roVSRPKFhVZGjLAgCfvb44yyqqoIJE1je20v1woVKjCJ5RD1HkaMoca/xtOefZ1NHB+UaThXJS+o5ihwFb9thY84cABrOPJNPrVmjXqNInlFyFBmlgXbY2FpRocQokqeUHEVGSTtsiBQeJUeRUXLgqUOHmN7Xx1PuvFs7bIjkPU3IERml2bW17AIeO3SIXeFzEclv6jmKjJJ3d7NkyhQqenqYWV6uKjgiBUA9R5FRskiETR0dnLhvH5s6OrCwfqqI5C/1HEWGKXV/Ru/u5pKaGigr45LeXvUcRQqAkqPIMKxfv54f3nYby485hnWVlbB2LWfMn8+6igpOA54eN46LtOBfJO8pOYqkKQgC7r31Vj7e2cnc8nKYOlX7M4oUKCVHkTTV19Vx7sGDPGUGPT08+OabfCXsJWp/RpHCogk5ImkIgoAXHnqIl2MxJvX28q3jjuNDt9yihChSoNRzFBlAYuLNSTU1bPvNb1hZUcGMefO4d8cO3nfVVaxatSrbIYpIhig5ioSSZ6ECrFu9miVAx6pVWCTCw8CSAwc4MHkyS5cty2qsIpJZSo4ivLXV1BJg3fr1nHTRRfHdNU46iR8CHd3dXK9JNyJFQ8lRBGiMRrkgFuO9YfHwJuBhgN27icHhNY1KiiLFQclRhHiVmw07d9IDfAdYUVvL0mXLDt9zVFIUKS5KjlKU0qlyk+gpNjQ0ZDtcERljWsohRSdxf/H49etZt3o1QRBwxvz5bK2o4Phx49haUXF4Uo6IFCf1HKVoJHqL27dvPzzZht27VeVGRN5GyVEKRupQaepridmoD3Z1sRVg924eBq5XlRsRSaHkKAUhdSkGa9cekegao9EjeoubL7qIvbNmqZcoIv3SPUfJe0EQcNcdd3BBLMayk05iCfFkmOyM+fN5GKgLe4tLly3j2htuUGIUkX6p5yh5LdFjPK2riw07dwLwZGXl4aHShAULFuieooikTclR8k7yvcXDw6Vz5gDQcOaZfGrNmn6Tn+4piki6NKwqeSV1GUai5mnd7t1sragYMDGKiAyHeo6SV+rr6jinvZ2LamrgwAH2quapiGSAkqPkjSAIeGbjRn7V1sbr7e08N2sWN6nmqYhkgJKj5I36ujou27OH0yMR7uvu5tizzlJSFJGMUHKUnJeYgLOrtZWXgFmlpXSXlzN72rRshyYiBUrJUXJa8uL+zq4u3qiu5id9ffRWVGjDYRHJGCVHyUkD1UHdfPnlnDJrFtdo8o2IZJCSo+ScQeugLlumpCgiGafkKFkxWJFw1UEVkWxTcpSjbrDEl3h9sCLhZ8yfHz+u3qKIZImSoxxVQyU+eHvPsDEaPeIc1UEVkWxT+ThJWxAEfPub3yQIggHPSU58/e2OAW/fIeOMlCLhEE+Q2jVDRLIlo8nRzC41syYz22ZmN/fz+jgz+4/w9cDMTs5kPDJyqTVNB0qQ6Sa+69euZe+qVVzfT89SRCTbMjasamalwFrgj4AW4Fkzq3f3l5JOWwXsdfc5ZnYV8P+AD2cqpmIWBAEde/YQBMGIktFQQ6EJ6Q6JquSbiOSyTPYczwG2uXuzu3cDDwBLU85ZCtwbPv4e8D4zswzGVJQSvb6yjo5Be32DSadHmKAhURHJd+bumbmw2XLgUne/Nnx+NbDA3W9MOufX4Tkt4fPfhed0pFzrOuC68Om7gV9nJOgCZVA9FaoOgZWCt0GHw54RXKrC4FiHN4Cuox1nDqsCOoY8S/qjthsZtdvI1br7caO9SCZnq/bXA0zNxOmcg7vfDdwNYGbPufvZow+v+KjtRkbtNnJqu5FRu42cmT13NK6TyWHVFmBm0vMZQOtA55hZGTAReC2DMYmIiAwpk8nxWeBUM5tlZhHgKqA+5Zx64Jrw8XLgp56pcV4REZE0ZWxY1d17zexG4FGgFLjH3V80s9uA59y9HlgP3G9m24j3GK9K49J3ZyrmIqC2Gxm128ip7UZG7TZyR6XtMjYhR0REJF+pQo6IiEgKJUcREZEUOZUcR1Nuzsy+FB5vMrP3j2Xc2TbSdjOzk83sTTN7Pvy5a6xjz7Y02u58M4uaWW+4djf5tWvMbGv4c03qewvZKNvtUNJ3LnWSXsFLo+0+Z2YvmdkLZvYTM3tH0mv6zo2s3Yb/nXP3nPghPmnnd8ApQAT4FfDOlHP+FLgrfHwV8B/h43eG548DZoXXKc32Z8qDdjsZ+HW2P0OOt93JwHuA+4DlScdPAJrD38eHj4/P9mfK9XYLX4tl+zPkeNtdCBwbPr4h6b9XfedG0G7h82F/53Kp5ziacnNLgQfc/aC7bwe2hdcrBirTN3JDtp27v+LuLwB9Ke99P/Bf7v6au+8F/gu4dCyCzgGjabdil07bPeHub4RPnyG+Rhz0nRtpu41ILiXH6cDOpOct4bF+z3H3XmAfcGKa7y1Uo2k3gFlmtsXMnjSzRZkONseM5nuj79xbhvvZx5vZc2b2jJl96OiGlvOG23argP8c4XsLyWjaDUbwnculzY5HU24urTJ0BWo07bYbqHH3TjM7C/ihmb3L3V8/2kHmqNF8b/SdO9JwPnuNu7ea2SnAT82s0d1/d5Riy3Vpt52ZrQDOBi4Y7nsL0GjaDUbwnculnuNoys2l895CNeJ2C4ehOwHcfTPxMf3TMh5x7hjN90bfubcM67O7e2v4uxloAOYdzeByXFptZ2YXA18GPujuB4fz3gI1mnYb2Xcu2zdak26YlhG/wTyLt264vivlnNUcObFkY/j4XRw5IaeZ4pmQM5p2q060E/Eb3buAE7L9mXKp7ZLO3cDbJ+RsJz4x4vjwcVG03Sjb7XhgXPi4CthKysSKQv5J87/XecT/oXpqynF950bWbiP6zmX9Q6d8iD8GXg4/4JfDY7cR/1cAwHjgu8Qn3PwSOCXpvV8O39cEXJbtz5IP7QYsA14Mv2hR4PJsf5YcbLv3Ev9XaxfQCbyY9N4/Cdt0G/CJbH+WfGg34DygMfzONQKrsv1ZcrDtHgd+Dzwf/tTrOzfydhvpd07l40RERFLk0j1HERGRnKDkKCIikkLJUUREJIWSo4iISAolRxERkRRKjiIpkir4/9rMvmtmx2bwb90WLlzGzNYM929Z3E/NbEL4/DNm9hsz+85RiG2lmU1Lev5tM3vnCK91o5l9YrQxiYwVLeUQSWFmMXevDB9/B9js7n+fxvuM+H9TIyq2bWavAGe7e8cw3rMEuNjdPxs+/y3xdb7bU84r83hd3eHE0wB83t2fG877BrjWscDP3b2YquFIHlPPUWRwm4A5cHi/uF+HP2vCYyeHPbV/Jl5IYaaZfcTMGsPz/l94XqmZbQiPNZpZIpltMLPlZvYZYBrwhJk9YWarzOwfEkGY2SfNrL8E/THgwfCcu4hXOqo3s8+a2VfN7G4zewy4L4x1U7jPYtTMzku6/hfDuH5lZn9j8T0Yzwa+E/aijzGzBjM7Ozz/bZ8xPB4zs78Kr/OMmU0B8PhuCa+YWbHsliP5LttVD/Sjn1z7Idz7jXjJqgeJ7w13FvHqGhVAJfHKQvOI71vYB/yP8D3TgB3ES/OVAT8FPhS+/7+S/sak8PcGwvJqwCtAVfi4gnglkPLw+dPAGf3E+ipwXNLz5Gt8FdgMHBM+PxYYHz4+FXgufHxZeP3EXngnhL8biPdkSX4+0GcMz3HCSkvA14GvJL3/y8BN2f7fVz/6SedHPUeRtzvGzJ4HniOeBNYDfwj8wN273D0GfB9IbPH1qrs/Ez5+L9Dg7ns8Poz5HeB84nUhTzGzb5jZpcCgO5+4exfxpPMBM/sD4kmysZ9TT3D3/YNcqt7d3wwflwPfMrNG4uUEE/cPLwb+xcO98Nz9tcFiG+QzAnQDPwofbyb+j4eEduKJVSTn5dKWVSK54k13PzP5wBCbQ3cln9rfCe6+18zmEt+wdjVwJfE6mYP5NvB/gN8C/zLAOb1mVuID3+dMju2zxGtPziV+S+VAUszDmXwwWFv0uHviWoc48v9jxgNvvv0tIrlHPUeR9DwFfMjMjjWzCuAK4vcjUwXABWZWZWalwEeAJ82sCihx9zrgz4H5/bx3P3Bc4om7B8S36fko8O8DxNVE/D5jOiYCu8NEejVQGh5/DPiTxExZMzuhv3iG+oxp/P3TgF+nGatIVik5iqTB3aPE7w/+knhy+La7b+nnvN3Al4AnCHc7cfcHie9a3hAO124Iz0l1N/CfZvZE0rGNxGd57h0gtIeBxWl+jH8GrjGzZ4gnqq4w5h8D9cBzYXyfD8/fANyVmJCTxmccykLiOyeI5Dwt5RDJYWb2I+Af3P0nA7x+EnCfu//R2EY2PGY2D/icu1+d7VhE0qGeo0gOMrNJZvYy8fuf/SZGONyL+1aiCEAOqyI+nCySF9RzFBERSaGeo4iISAolRxERkRRKjiIiIimUHEVERFIoOYqIiKT4/5/a2auimD7EAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sort the data:\n",
    "sort = np.sort(X)\n",
    "\n",
    "# calculate the cumulative probabilities assuming known tails\n",
    "p = np.arange(len(X)) / (len(X) - 1)\n",
    "\n",
    "# plot the cumulative probabilities vs. the sorted porosity values\n",
    "plt.subplot(122)\n",
    "plt.scatter(sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7)\n",
    "plt.xlabel(feature + ' ' + funits); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric CDF')\n",
    "plt.ylim([0,1]); plt.xlim([0,0.25])\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit a Gaussian Distribution\n",
    "\n",
    "Let's fit a Gaussian distribution\n",
    "\n",
    "* we get fancy with Maximuum Likelihood Estimation (MLE) for the Gaussian parametric distribution fit mean and standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXhU5fn/8fedjWBC2BJACJssKVZkUURFLK6tXwWqWEXFiq0LlGqp1v5stVqXVkX5ggsi0CjiUkWwCmqVr9UALh2RoMYNxbAFCJAQlgkJ2Z7fHzPBELMMy2QyM5/XdeVi5pwzJ/echNxzP+dZzDmHiIiIfC8m1AGIiIg0N0qOIiIitSg5ioiI1KLkKCIiUouSo4iISC1KjiIiIrUoOYocIjN7wsz+EuIY/mpmzx7C6/5sZv8IRkxNxcyuMLMloY5DIpOSYxQzs3VmttXMkmpsu8bMskIYVsiZ2Xgze6+x45xzE5xz9zRFTIfCzEaYWZWZeWt8LQZwzv3dOXeN/7geZubMLK6R8/UxsxfMbLuZ7Tazb83sUTNLb4r3U5tz7jnn3LnBOLf53Ghmn5tZsZnlmdlLZtbfv3+umZWZ2R7/1+dmdp+Zta5xjvFmVlnr+j8WjHjlyFNylDjgd6EOoiGN/dEOBTOLDXUMAdrsnEuu8TXyUE5iZr0BD7AZGOScSwGGAd8Bpx25cJuNh/H9v7gRaAf0BV4Bzq9xzBTnXCsgDbgaOBl4v+aHTeDDWtf/t00TvhwuJUd5EPiDmbWpa6eZnWpmK8xsl//fU2vsyzKze8zsff+n5yVmlurfV12NXGdmm81si5ndXOO1J5nZh2a207/vMTNLqLHfmdkkM/sW+Na/7WEz2+ivWlaa2fAax//V/8n+WX8sOWbW18z+ZGbb/K87t8bxrc0s0/+9N5nZvWYWa2b9gCeAU/yf9Hf6j59rZjPN7A0zKwbO8G+7t8Y5R5vZJ/74vjOzn9VzTW/1799jZl+a2YU19o03s/fM7CEzKzKztWZ2Xo39Pc1sqf+1/wekNv4jrjOGms2xy/z/7vS/51PqeMlfgfedczc55/IAnHPbnHPTnXMv+M/Z1sxe81eWRf7H+6tK87VUnF1XDGaW6P/ZFfp/J1aYWcca1yTX/57XmtkVNa9VjfM19vsx38zm+c/zhZmdWM+16QNMAi5zzr3jnNvnnNvrr1Tvr328c67UObcCGAW0x5coG2Rm/+P/2e/x//79obHXSNNScpSPgSzgB/85zawd8DrwCL7/9P8LvG5m7Wscdjm+PwYdgIQ6znMG0Ac4F7i1xh/HSuD3+P64nwKcBfym1mt/DgwFjvU/XwEMxPdJ/nngJTNLrHH8SOAZoC2wCngL3+94F+BuYFaNY58GKoDewCB/fNc4574CJvD9J/6aHxouB/4GtAIOaHY1s5OAecAtQBvgdGAddfsOGA60Bu4CnjWzo2vsHwqs9l+bKUCmmZl/3/PASv++e4Cr6vkeB+N0/79t/O/5wzqOORtY2Mh5YoCngO5AN6AECLQZ8Sp816Mrvt+1CUCJ+aqwR4Dz/FXaqcAn9Zyjsd+PUcAL+H4+ixqI7Swgzzn3UYCxA+Cc2wP8H76fbWMygev97+k44J2D+V4SfEqOAnAHcIOZpdXafj7wrXPuGedchXPun8DX+JJQtaecc98450qA+fj+ONV0l3Ou2DmXg+8P52UAzrmVzrn/+s+7Dl/i+kmt197nnNvhPzfOuWedc4X+10wFWgAZNY5f7px7yzlXAbyEr7nrfudcOb4/ij3MrI2/IjkPmOyPbRswDRjbyHV61Tn3vnOuyjlXWmvfr4EnnXP/59+/yTn3dV0ncc695Jzb7D/uRXyV8Uk1DlnvnJvjnKvEl8SPBjqaWTdgCPAXfzWzDFjcSMyd/ZVY9dcljRxfn1Qgv/qJmf3Wfz6vmc3xv69C59xCf5W1B98Hido/0/qU40uKvZ1zlf7fj93+fVXAcWbW0jm3xTn3RV0nCOD34z3n3Bv+6/oMMKCeWNoDWwKMu7bN+JJztZNrXf+T/dvLgWPNLMU5V+Scyz7E7ydBouQoOOc+B14Dbq21qzOwvta29fgqsWr5NR7vBZJrHb+x1ms7A/ibPF8zs3wz2w38nR82EdZ8LWZ2s5l9Zb4m3p34Ko2ar9la43EJUOD/Q1j9HH983YF4YEv1Hy18ybkDDdvYwL6u+CrCRpnZL/3Nr9Xf+7ha72P/NXXO7a0Rd2egyDlXXOPY2j+f2jY759rU+JofSIx1KMSXpKvjesxfVU/Hdy0xs6PMbJaZrff/TJcBbSyw+7PP4Kv0XzBfM/wUM4v3v9dL8VWSW8zsdTP7UV0nCOD3o/bvaqLVfT/7gPd6kLoAO2o8/2+t6/9f//YxwP8A6/3N5HU1ZUsIKTlKtTuBazkw8W3Gl0hq6gZsOojzdq312s3+xzPxVaF9/J07/gxYrdfuXzLGf//o/wGXAG39f5h31fGaQGwE9gGpNf5opTjnflz7+9YXTz3n7NXYNzaz7sAc4LdAe//7+JzA3scWoK0d2OGjWwCva0wgS/P8B7iokWNuxlepDfX/TKuba6vfWzFwVI3jO+0PwLly59xdzrlj8TWdXgD80r/vLefcOfgS1tf4rt8BjvDvx3+A9PruSdbHzJLxNT8vb+xY59wK59xofB/IXsHX6iLNiJKjAOCcWwO8iK93XrU3gL5mdrmZxZnZpfju/712EKf+i7+i+DG+e5Mv+re3AnYDXn8lMLGR87TCd49wOxBnZncAKQcRx37OuS3AEmCqmaWYWYyZ9TKz6ibArfj+OCbUf5YfyASuNrOz/OfrUk+Fk4QvGW0HMLOr8VWOgcS9Ht894rvMLMHMTuPAJu5DtR1f0+UxDRzzV2C4mf2vmXUBMF/nq341jmmFr0Lf6b9ffWetc3wCjDWzeH/iubh6h5mdYWb9/VXmbnzNjpVm1tHMRvk/EOwDvPjuV9d2JH8/vgUeB/5pvuEwCebrMDTWzGq3rmBmLczsBHxJrgjf7YN6+c93hZm19jf5767nPUkIKTlKTXfj++MN+O4h4fsEfzO+pqY/Ahc45woO4pxLgTX4Po0/5JyrHrT9B3wdXPbgqwRerPvl+70F/Bv4Bl9TYikNN3M25pf4OhB9ie8P2gK+b0p7B/gCyDezgN6rv/PG1fjuXe7C975rV904574EpgIf4kvC/YH3DyLuy/F12NmBL/nMO4jX1snfdPs3fMMQat4Xq3nMN/iGKqQDn5rZHn/cm4HqiRCmAy2BAuC/wJu1TvMXfNV1Eb6OSM/X2NcJ389gN/AVvuv3LL6/UTf7v88OfPcwa3fcgiP/+3Ejvg47M4Cd+JrML+TAe7x/9F+HHfh+DiuBU2s1e9fnSmCdv/l5AjDuMGKVIDCnxY4lCMysB7AWiPd3kBERCRuqHEVERGoJWnI0syfNN/j683r2m5k9YmZrzOwzMxscrFhEREQORjArx7lAnTOE+J2Hb3B4H+A6fL0XJUI459Y550xNqiISjoKWHP0DlHc0cMhoYJ7z+S++8VCHOrZIRETkiAnlhM5dOLA3WZ5/2w9mpjCz6/BVlyQmJp7QrduRGNoVfaqqqoiJ0W3mg6Xrduh07Q5NsK9bVVXV/u/R2PepqKjAKiqINaPSOVxcHHFxh5c6Kisrf/BVVVW1/9/qr8rKSpxz+59XP26kI2mBc672bF8HLZTJsa7BuXW+Y+fcbGA2QEZGhlu9enUw44pYWVlZjBgxItRhhB1dt0Ona3dognndPB4PsyZN4nx8EydfP2MGQ4cOPSLHl5eXs2HDBtatW8f69etZv349eXl5bNmyhS1btrB582YKCgqoqqqq8/UxMTG0bt2atm3bkpKSQqtWrUhOTiYpKWn/vy1btjzgKzExkRYtWuz/uvTSSxubNSogoUyOeRw4e0o638+eIiIiQZCTnc1PvF6GxMfjLS8nJzu7weQ4dOhQmDGDnOxsrh88mKFDh7Jr1y4+//xzcnJy+Oqrr/jmm2/49ttvWbduHZWV389nYGZ07NiRLl26kJ6ezkknnUTHjh1JS0sjLS2N1NRU2rdvT9u2bWnXrh2tWrXi+zn2D82ll156WK+vFsrkuAj4rZm9gG9Q8y7/zCUiIhIklpDAk+vXs7eykhdiY/llQsMTQRUVFbF792627tjBa/fdR3Z2Nhs3fn9HLDk5mT59+nDCCScwduxYevXqRY8ePejevTvp6ekkNHL+5ipoydHM/gmMAFLNLA/fbB7xAM65J/BNTfY/+GZP2UsAa6CJiMjhyV29mgHOURwbywDnyK11m2rr1q0sXbqUpUuXkpWVxZdffrl/X9++fTnttNPo378/xx9/PP3796dr166HXe01R0FLjs65yxrZ7/AtKHrYysvLycvLo7S09ipCkSUxMZH09HTi4+NDHYqIhKlNmzdTUFXFsIQE3neOjKoq3nvvPV5//XXeeOMNPvvsM8BXEQ4bNozLL7+coUOHcuKJJ9KmTZ1rokekUDarHjF5eXm0atWKHj16ROQnGADnHIWFheTl5dGzZ89QhyMiYcjj8VC4ciW7zJhRWsq3iYm8M2cOD0ydSmxsLKeddhr3338/I0aMYPDgwVH9QTwikmNpaWlEJ0bw3dhu374927dvD3UoIhKGnHO8smAB3i1b+LSigj1VVbQoL+cXF17IyJEjOffcc6OqMmxMRCRHIKITY7VoeI8icmTt3LmTZ599ltmzZ5OTk0MMcEJsLCUJCUx6+GEmTJgQ6hCbpYhJjiIi0cTj8ZCTnU1///CK2tasWcPUqVN5+umnKSkp4cQTT+TKyy/nmBUrSG/Rgu0VFcRpVaZ6aeqKI+SRRx6hX79+tG3blvvvvx+AV1555YCeXiIiR4LH42Hq+PFU3HcfU8ePx+Px7N+3YsUKfvGLX9C3b1+eeuoprrjiClauXMmKFSuYdOONbEhJoW2LFnyblET/wVrvoT6qHI+Qxx9/nH//+98HdJZ55ZVXuOCCCzj22GNDGJmIRJpFCxdy0oYNjI6Px1tYyKKFC0lMTOTPf/4zb7zxBq1bt+bWW2/lxhtvpFOnTvtfV9eAfqmbkuMRMGHCBHJzcxk1ahS/+tWv+O6777j88stZtGgRS5cu5d5772XhwoX06tUr1KGKSBipr+nUAe8DPZ1jSVUVm157jb8/+CBt2rThvvvu4ze/+Q0pKSl1nnPo0KFKigGIuOQ4efJkPvnkkyN6zoEDBzJ9+vR69z/xxBO8+eabvPvuu7z22msAnHrqqYwaNYoLLriAiy+++IjGIyKRr+acprMyM6HGnKajx4xh1auvcs+WLXxWWkpCbi5/+tOfuOWWW2jbtm1oA48QEZccRUTCTV0VYk52NucDY44+GrZs2T8HqnOOvLw8snfuZNuePYwcOZJZs2Zx9NFa8e9Iirjk2FCFJyLS3NRVIQL0HzzY93zLFt9qGIMHs379eiZMmMCbb77JgAEDeOWVVzjllFNCGn+kirjk2Jy0atWKPXv2hDoMEWnG6lolo3e/fgd0nrlu0CC+/vprzjnnHJxzTJ8+nUmTJh32uopSP13ZIBo7dizXXnstjzzyCAsWLFCHHBH5gbV5eSxbu5ZSM16IjWVcjVUshg4dSq9evbj++ut5+eWXGT58OPPmzaNHjx6hCzhKKDkeIevWrQNg/PjxjB8/HoBhw4ZpnKOI1Mvj8bD8mWcYbsZmM4anpuLKyvbvf+edd7j88sspKipiypQp3HTTTcTGxoYw4uihSQBEREIkJzubEWasN6Ozc6zCd6/ROceUKVM455xzaNeuHStWrOCWW25RYmxCqhxFRELEEhJYXlDAKc7xcFUVF44bR79+/bjgggtYvnw5l1xyCZmZmSQnJ4c61Kij5CgiEiK5q1dzUnIyHZOTuTIhgaqEBIYOHco333zD1KlT+f3vf68FB0JEyVFEJAQ8Hg+fLV5MnNdLmtfL22lpZE+fTlxcHFOnTmXy5MmhDjGq6Z6jiEgILFq4kPOLi/lN797856ijWLZpE+3bt+fDDz9k4MCBoQ4v6qlyFBEJorpmv6muGr8sLIT8fN4sL+fYY48lKyuLtLQ08vLyQhy1qHIMkr/+9a889NBD9e7XclYika969pu2mZnMmjRp/9JSOdnZjE9KokNqKq+Ul9O7Vy9WrFhBWlpaiCOWakqOIaLkKBL59s9+s28fP/F6ycnOBuC4QYO4d/NmZm/YQK927Zj79NMcddRRIY5WalJyPIL+9re/kZGRwdlnn83q1asBmDNnDkOGDGHAgAGMGTOGvXv38sEHH7Bo0SJuueUWBg4cyHfffVfncSIS3iwhgbkbN7IkN5e5GzdiCQk453j99df5ZMsWhp1yCvMWL2bYsGGhDlVqidrk6PF4+MfMmQesoH04Vq5cyQsvvMCqVat4+eWXWbFiBQAXXXQRK1as4NNPP6Vfv35kZmbuX87qwQcf5JNPPqFXr151Hici4c2VlTE8NZXC1q0ZnppK1b59/PnPf+bee+/lmmuuYdl773HqqaeGOkypQ1R2yGlonbRDtXz5ci688ML9TSOjRo0C4PPPP+f2229n586deL1efvrTn9b5+kCPE5HwsTYvj2X5+Ywz483YWFovWcKrr77KhAkTmDFjBjExUVufNHtR+ZOpuU7a+f7nR0Jdg3XHjx/PY489Rk5ODnfeeSelpaV1vjbQ40QkPNSeN/WoFi149dVXueqqq5QYw0BU/nT6Dx7M68BC/zpp/QcPPuxznn766fzrX/+ipKSEPXv2sHjxYgD27NnD0UcfTXl5Oc8999z+42svZ1XfcSISnmrOm7qlooI3iooYNmwYc+bMUWIMA1HZrFpznbTra4w9OhyDBw/m0ksvZeDAgXTv3p3hw4cDcM899zB06FC6d+9O//799yfE2stZ1XeciISn6nlTj66s5B8VFaSnp7NkyRLi4+NDHZoEICqTI/gS5JFIijXddttt3HbbbT/YPnHixB9sq72c1cSJE+s8TkTCj8fj4b2336ZvSgrztm2jQ4sW/GHyZA3XCCNRmxxFRIKhusNf5127mLp1Ky1jY/lJ376cfNppoQ5NDoIavkVEjqBFCxcyeOtW3t22jUrgtDPP5KY5c454S5UEV8RUjs65iF/axTkX6hBEpA7V86daQgKfLlrEc1u2sL6ykpM7d+b2u+9WYgxDEZEcExMTKSwspH379hGbIJ1zFBYWkpiYGOpQRKSGmuOmH962ja779rG+spITUlI444orlBjDVEQkx/T0dPLy8ti+fXuoQwmqxMRE0tPTQx2GiNSwaOFCTtq2jTO7deNN4B/btnFSmzYce8wxjB4zJtThySGKiOQYHx9Pz549Qx2GiESZmktPbSko4NmyMjp16sSVf/wjQ049VVVjGIuI5Cgi0tQ8Hg9PTJ/Oxc7RZ+BALlm5kgozlixZQv/+/UMdnhwmJUcRkYNUfZ+xb3ExczdupG1BAZv8k4orMUYGDeUQETlI1fcZJ7RvT6/27flXYSHnnXce9957b6hDkyNElaOISCOqh2pUz8NcfZ8x33+fsUePHrz00ksR21s+Gik5iog0oPYSd0efeSbjk5LoMnAgl61aRQWwaNEikpKSQh2qHEFKjiIiDag5VIPSUlYCrwNHbdrEupISbrjhBt1njEBKjiIi9ag5VMNbWMhH3bpx85gxbD7lFC6//HIGDRrE9OnTQx2mBIGSo4hIPXKys7nYOdp06sQir5feI0cyZMgQzjzzTOLi4nj55Ze1NmOEUnIUEamHJSQwd+NGrgBygVMzMpg+fTpLly4lMzOTHj16hDhCCRYlRxGReriyMoanplJYXs7w+Hi2bNzIvfffz6hRo7j66qtDHZ4EkZKjiEg91ublsSw/n3Fm/Dsmhi3//CctW7Zk1qxZGrYR4ZQcRUTq4PF4WP7MMww3Y7MZrVq2ZOk33zBnzhw6deoU6vAkyHQnWUSkDosWLuSUffvIi4sjxYy3iooYNGgQv/71r0MdmjQBVY4iIrVUD+GI83rpWlnJ3wEXG8s///lPNadGCVWOIiK1VA/hGN+pE18lJrJj3z7uvPNOMjIyQh2aNBFVjiIitVQP4RjjHMtLS+ncuTN//OMfQx2WNCFVjiIiteSuXs1JycksatGCfc4xbuxYEhISQh2WNCFVjiIiNWRmZrJ03jxa7NzJ0rIyuqWkcNEll4Q6LGliqhxFRPw8Hg9P33UX4/fsYUdlJWbG6CuuYOjQoaEOTZqYkqOIiF9OdjYjzHi2vJxPKivp07o1V1x1VajDkhBQs6qIiJ8lJLB0+3ZyKipIAEZdf72qxiilylFExM+VlRHXogVFzjEmLY3e3buHOiQJESVHERG/nK+/JmvnTo6PiWFLcTGmHqpRS82qIiL4OuPMz8wEYHhcHO1SU3FlZSGOSkJFlaOICPDk7NlsLS6mT1wcA+PiWAX0Hzw41GFJiCg5ikjU83g8vDp/PrHA8JgY5rRqxc/vuEOdcaJYUJOjmf3MzFab2Rozu7WO/d3M7F0zW2Vmn5nZ/wQzHhGRuix88UW2er2MbduWqjZtOOuXv9TqG1EuaMnRzGKBGcB5wLHAZWZ2bK3DbgfmO+cGAWOBx4MVj4hIXZxzLHjlFRLMOLm0lFyvl16aYDzqBbNDzknAGudcLoCZvQCMBr6scYwDUvyPWwObgxiPiMgPvPTSS6xdu5ZR7drhbdGC4fHx6ogjmHMuOCc2uxj4mXPuGv/zK4Ghzrnf1jjmaGAJ0BZIAs52zq2s41zXAdcBpKWlnTB//vygxBzpvF4vycnJoQ4j7Oi6Hbrmfu3Ky8sZP348sbGx/OnGG0mLiWGHGe26diU1NTVkcTX369acnXHGGSudcyce7nmCWTnWtSJo7Ux8GTDXOTfVzE4BnjGz45xzVQe8yLnZwGyAjIwMN2LEiGDEG/GysrLQtTt4um6Hrrlfu1mzZrF582Z+nJrKxttvJ98M16EDMbfeysUXXxyyuJr7dYsGweyQkwd0rfE8nR82m/4amA/gnPsQSARC93FNRKJGaWkp9957L106d+Z/zMiLi6NrTIyGcAgQ3OS4AuhjZj3NLAFfh5tFtY7ZAJwFYGb98CXH7UGMSUQEgNmzZ5OXl0fX2Fi+LS6mTUWFhnDIfkFrVnXOVZjZb4G3gFjgSefcF2Z2N/Cxc24RcDMwx8x+j6/JdbwL1k1QERG/vXv38ve//52MPn2YEBNDm/h4Fnm9GsIh+wV1+jjn3BvAG7W23VHj8ZfAsGDGICJS24wZM9i6dSvjx49n7qOPcgWQC5yqIRzip7lVRSSq7NmzhwceeICf/exn9O7enYTUVArLyzWEQw6g6eNEJKo8/PDDFBYWcs899/jWb9y2jZSCApZu26ZVOGQ/JUcRiRq7du3ioYceYvTo0Zx44onkrl7NAOcojo1lgHPkrl4d6hClmVByFJGo8dhjj7Fr1y7uvPNOADZt3kxuVRVdY2PZGBv7g4HYEr2UHEUkKni9XqZNm8b555/PoEGD8Hg8FK5cyZ6YGF4uK2N3Whqjx4wJdZjSTKhDjohEhSeeeILCwkJuv/12AHKys7nYOdp06cIir5del1yi8Y2ynypHEYl4JSUlPPTQQ5x99tmcfPLJAFhCAnM3bmR7fr5W4pAfUOUoIhEvMzOTrVu38uKLL+7f5srKGK5hHFIPJUcRiWj79u3jgQce4LTTTuP000/H4/GQk53N2rw8Pigo4ArgTWCchnFIDUqOIhLR5s2bR15eHpmZmXz00UfMmjSJ84F527YxIjUVUlI4t6JClaMcQPccRSRiVVRUcP/99zNkyBDOOecccrKz+YnXy5B9+xhhxiqgbYsWfJuUpJU45ACqHEUkYi1cuJDc3FweeughzGx/J5xyYDlw6uTJFKWnc/3gweqpKgdQchSRiOScY8qUKfTt25fRo0fj8Xh47+23GV6jKTUtPZ1rJk4MdajSDKlZVUQi0jvvvEN2dja33HILK1asYNakSfT95BOWFxQQX16uplRpkCpHEYlIU6ZMoVOnTvTt25cnpk9nhNfLVX36AJA1cCATJk9WU6rUS8lRRCLOJ598wpIlS7jooot48MorGQTMLSgA4NvkZCVGaZSSo4hEnIceeoijjjqKLR4Pvyoq4sT4eEhNJWvQICVGCYjuOYpIRFm/fj0vvPACw045hbNjY3m3qgpPaSmrQIlRAqbkKCIRZdq0aZgZ5/z0pywvKKC7czxcVcXx48YpMUrA1KwqIhHB4/Hgee89Zs+ezWWXXUbb5GTO7daNtLg4rvQP2xAJlJKjiIQ9j8fD1PHjqdiyhZKSEs455xz69u3LrKQk+gIftGjBmRq2IQdByVFEwt6ihQs5Yf16HiktpVdMDF/n5HDllVfCjBnkZGdrBhw5aEqOIhL2HPBSeTmbnWNIbCzOv33o0KFKinJI1CFHRMJer4wMvq6spBVQEROjhYvlsKlyFJGwt2b1aoqdY2RKCgPbtNHyU3LYVDmKSNj7T1YWscC55eUsLyjAtHCxHCZVjiIS1jZu3MjHH3/MiS1bUnn00ZybkKDKUQ6bKkcRCVsej4dLLr4Y5xxtnaM8P59VoNU25LApOYpIWPJ4PMyYOJFVH39M/9hY/tC3L7nt29N75Ej1UJXDpuQoImEpJzubuM2b2VdVRauYGL4rKqK0QwdGjxkT6tAkAig5ikh4io9n/rZtpJtRYobn9NO5fsYMVY1yRKhDjoiEpWXvvEOxc1yQmsqA1FTShg1TYpQjRpWjiIQdj8fDf157jXjg+OJidcKRI06Vo4iEneeefpr8PXu4MDWVvJYt1QlHjjhVjiISVjweD4uff54qoNfu3RQkJakTjhxxSo4iElb+NX8+hbt2MSQ2lnXOcdQJJ6hqlCNOyVFEwsqHHg97gLPi4iiLj6dz586hDkkikJKjiIQNj8dDTnY2iWZ8V1nJ7rQ0NalKUKhDjoiEjX8vXkxRSQlXt2uHxcXR+5JL1KQqQaHKUUTCxrIPPsCAASUl5Hq9WrdRgkaVo4iEBa/Xy4cffkj/pCRKU1IYHh+v1TckaFQ5ikhYeP755yktLSWmvJyUggKWbtumdRslaJQcRaTZc84xc+ZMOqSlcRpQHBvLAOfIXVz5yS4AACAASURBVL061KFJhFJyFJFm76OPPuKTTz7hmB49WOscXWNj2Rgbiwt1YBKxGk2OZvaxmU0ys7ZNEZCISG0zZ87kqKOOos3OneyJieHlsjIN45CgCqRyHAt0BlaY2Qtm9lMzsyDHJSICwI4dO3jxxRfJ6NOH0aWl3J6RQftOnRiqYRwSRI0mR+fcGufcbUBf4HngSWCDmd1lZu2CHaCIRLenn36a0tJSWu3axVuFhXy6Zo3mU5WgC+ieo5kdD0wFHgQWAhcDu4F3gheaiEQ75xyzZs2iy9FHc1llJb/p3Zvc9u21CocEXaPjHM1sJbATyARudc7t8+/ymNmwYAYnItFt1qxZrF69mn7t2/PWrl0MKyykoFs3rlLVKEEWyCQAv3DO5dbcYGY9nXNrnXMXBSkuEYlyHo+Hh267jQQzxu3dy5A+ffhXUZGqRmkSgTSrLghwm4jIEfPf5ctZX1TEma1a8d+qKr4rKqK0Qwfda5QmUW/laGY/An4MtDazmhViCpAY7MBEJLr9d+VKKpxjSFkZrwOe009nwg03qGqUJtFQs2oGcAHQBhhZY/se4NpgBiUi0c05x9tvv02nuDjad+7MxQkJpA0bpsQoTabe5OicexV41cxOcc592IQxiUiUu+222ygoKOD4uDjK8/NZ1a0bNw8eHOqwJIo01Kz6R+fcFOByM7us9n7n3I1BjUxEopLH4+HJRx6hBTApIQFPcrI64UiTa6hZ9Sv/vx83RSAiIgCe99+ncO9eusfG4ior2ZaYyHXqhCNNrKFm1cX+f59uunBEJNp9lJ1NhXOcFhPDw1VVXDhunKpGaXINNasuhvonvXfOjQpKRCIS1d7+v/+jQ1wcA9PTyWjRgrT09FCHJFGooWbVh5osChERfPOobt22jePj46nYulUdcSRkGmpWXdqUgYiIzH3ySeLM+GPnzrxTUqKOOBIy9c6QY2bz/f/mmNlnNb5yzOyzQE5uZj8zs9VmtsbMbq3nmEvM7Esz+8LMnj+0tyEi4a6kpATPihW0NWP35s2s2b2bXhkZoQ5LolRDzaq/8/97waGc2MxigRnAOUAevvUgFznnvqxxTB/gT8Aw51yRmXU4lO8lIuHv5ZdfpqSkhJ/Gx1McG8sA58hdvTrUYUmUqrdydM5t8f+7HtgHDACOB/b5tzXmJGCNcy7XOVcGvACMrnXMtcAM51yR/3ttO/i3ICKRIDMzk+TkZMqco2tsLBtjY+vvESgSZIEsWXUNcAe+tRsNeNTM7nbOPdnIS7sAG2s8zwNq3zzo6/8e7wOxwF+dc2/WEcN1wHUAaWlpZGVlNRa21MHr9eraHQJdt0MX6LXbtGkT7777LqMvuIDzzjiDYuC8+Hi6HHNMVF57/c6FXiBLVt0CDHLOFQKYWXvgA6Cx5Gh1bKv9QTAO6AOMANKB5WZ2nHNu5wEvcm42MBsgIyPDjRgxIoCwpbasrCx07Q6ertuhC/Ta3XbbbZgZP/nySzp9+SWLvF46XnUV191wQ/CDbIb0Oxd6gSxZlYdvsvFqeziwImzodV1rPE8HNtdxzKvOuXLn3FpgNb5kKSJRoqKigrlz53LcccexKD+f7fn55Hq96owjIdXQJAA3+R9uAjxm9iq+ym808FEA514B9DGznv5zjAUur3XMK8BlwFwzS8XXzJqLiESNt956i82bN3PykCH03rqV4uRkzk1IwJWVhTo0iWINVY6t/F/f4Uti1U2irwJbGjuxc64C+C3wFr55Wuc7574ws7vNrHp2nbeAQjP7EngXuKW6+VZEokNmZiZt27Zlx6pVfFFUxL5Nm1gF9NfgfwmhhiYBuOtwT+6cewN4o9a2O2o8dsBN/i8RiTLbtm1j8eLFDB44kPO//JJ+CQnMKyvjqBNO0OB/CalAequmAX8EfgwkVm93zp0ZxLhEJAo899xzVFRUcHz//izLyaGLGSUxMfTq3DnUoUmUC6RDznPA10BP4C5gHb77iSIih8w5x5NPPsmQIUM4edgwNgFLKivZBOqMIyEXyFCO9s65TDP7nX++1aVmpnlXReSwrFy5ks8//5yZM2eSu3o1Z7duTefkZPqqM440A4FUjuX+f7eY2flmNgjfsAwRkUP25JNPkpiYSO/evfls8WK+8Xopz89XZxxpFgKpHO81s9bAzcCjQArw+6BGJSIRraSkhH/+859cdNFFvLtkCecXF9Ord2/+VVSklTikWWg0OTrnXvM/3AWcEdxwRCQavPLKK+zcuZOTTz6ZV6ZO5dP8fE7fto2Cnj25asyYUIcn0nizqpkdY2aLzazAzLaZ2atmdkxTBCcikempp56iW7dubNmwgfO2b2diQgIrq6o0hEOajUDuOT4PzAc6AZ2Bl4B/BjMoEYlcGzZs4O233+b4445j85YtvA+UxsZSFh9PZw3hkGYikORozrlnnHMV/q9n+eEE4iIiAfnb3/6Gc47zNmygcOVKdqel8Z/Wrano1o3RalKVZqKhuVXb+R++a2a34luP0QGXAq83QWwiEmGqqqqY/+KL9G3Rgstat6ZjaSkrR47kmJ49uWrwYDWpSrPRUIeclfiSYfXSU9fX2OeAe4IVlIhEptmzZ7Nz1y6OSUggc9UqPurWjZvHjFFSlGanoblVezZlICISuTweDznZ2cx8/HFaxsQwdfBg5m/apGEb0mwFMrdqPDARON2/KQuY5Zwrr/dFIiJ+Ho+HqePHc+qePXyxeTPdUlIoLC+ntEMH3WOUZiuQSQBmAvHA4/7nV/q3XROsoEQkcixauJCTNmzAnKPSOTJOOomiCy/ket1jlGYskOQ4xDk3oMbzd8zs02AFJCKRxQHLKitZU1ZGkhkDBwzgmokTQx2WSIMCGcpRaWa9qp/4JwCoDF5IIhJJemVkkOscXzlHSkwMvX/0o1CHJNKoQCrHW/AN58jF13O1O3B1UKMSkYjhyspon5iIlZXxi44dteKGhIUGk6OZxQAlQB8gA19y/No5t68JYhORCJC7cSOe3bv5kRmf7dzJ8QkJoQ5JpFENJkfnXJWZTXXOnQJ81kQxiUiEKC4u5rV//IN9QL+4OH6cmqrKUcJCIPccl5jZGDOzxg8VEfleyd692N69xAMjzLRWo4SNQO453gQkARVmVoqvadU551KCGpmIhL2S0lI+Ly5mQEwMM53jwnHjNHxDwkIg6zm2aopARCSyeDwest59lyrgpz16kJKQQFp6eqjDEglIQxOPdwD+DPTGd7/xfufc7qYKTETCl8fjYdakSXy4axctzciIi2NZUhJnqklVwkRD9xznAcXAo0Ar4JEmiUhEwl5Odjb9duzg6zVrODklhaWDB3P9jBlqUpWw0VCzaifn3G3+x2+ZWXZTBCQi4c8SEnhswwbMjL2lpZx29tlKjBJWGqoczczamlk7/9qOsbWei4jU6buvv6aospL+ffpwEpC7enWoQxI5KA1Vjq3xrelYcwhHdfXogGOCFZSIhLf1GzeyBzhzyBC+W7+e40IdkMhBamg9xx5NGIeIRJC1GzYQA/Tr359lCxYwOiMj1CGJHJRAJgEQEQnYnj17WLlyJQOTk0lp2ZLzNZ+qhCElRxE5ohYsWEBZWRmuvJy4ykqWFxRgmk9VwkwgM+SIiARs7ty5tGvblrPj4qiKj+fcbt1UOUrYCahyNLPTzOxq/+M0M+sZ3LBEJBx99913LFu2jDbO8eXOnbiyMs2nKmGp0crRzO4ETsS3ZNVTQDzwLDAsuKGJSLiZN28eAJeVlnJKQgJbgaNOOEFjHCXsBFI5XgiMwjdbDs65zfhmzBER2a+qqoqnn36azkcfzSeVlXgBZ0bnzp1DHZrIQQskOZY55xy+sY2YWVJwQxKRcLR06VLWr19PRzP2xMTwclkZlXFxjB4zJtShiRy0QDrkzDezWUAbM7sW+BUwJ7hhiUi4mTt3Li0TE5mYlESHxEQWeb0ktWunJlUJS4EsWfWQmZ0D7MZ33/EO59z/BT0yEQkbXq+XhQsXMuSkk3j+44+5AsgFzmvRItShiRySQDrk/B54SQlRROqzYMECiouL6ZGeTqevv6Y4OZlzExLAuVCHJnJIArnnmIJvVY7lZjbJzDoGOygRCS9z586lY8eOrHn7bb7atYvy/HxWAS2POirUoYkckkaTo3PuLufcj4FJQGdgqZm9HfTIRCQs5ObmsnTpUlqWlnK118sv4uL4KjmZ3iNHkpSk/nsSng5m+rhtQD5QCHQITjgiEm6efvppzIzRSUm8W1XF3ooKtiUmqpeqhLVGk6OZTTSzLOA/QCpwrXPu+GAHJiLNX/XYxn79+vHpzp10d46Hq6o4ftw49VKVsBbIUI7uwGTn3CfBDkZEwktWVhbr169n9AUXkFFQQMfkZK5MSCAtPT3UoYkclnorRzNL8T+cAmwws3Y1v5omPBFpzp566imSk5PZ+emnfFFUxL5NmzSXqkSEhirH54ELgJX4ZsexGvsccEwQ4xKRZm737t0sXLiQH/XtywXffku/hATmlZVpLlWJCPUmR+fcBf5/tQKHiPzA/PnzKSkpof9xx7Hs66/pYkZJTAy9NJeqRIBAOuT8J5BtIhJdnnrqKfr168fwESPYBCyprGQT0CsjI9ShiRy2eitHM0sEjgJSzawt3zerpuAb7ygiUWr16tV88MEHTJo0iff/8x/O79iRbikp9K2o0MLGEhEauud4PTAZXyJcyffJcTcwI8hxiUgz9vTTTxMbG0vB8uUMKCtjSUEB41u25NvkZM5UZxyJAA3dc3wYeNjMbnDOPdqEMYlIM1ZZWcm8efPo2b07I4qKuLRbNwCyBg5kwuTJ6owjESGQVTkeNbPjgGOBxBrb5wUzMBFpnh599FE2bdrEcampvLV3L97CQlZ168bNSowSQQJZleNOYAS+5PgGcB7wHqDkKBJlPB4PD999Ny1iYriiuJgT+vThX0VF9B45UolRIkogc6teDJwF5DvnrgYGAFqkTSQK/Xf5cvJ27uSc5GQ+qKriu6IiSjt00DyqEnECSY4lzrkqoMI/a842NAGASFT6aNUqKpxjcFkZmwDP6adz/YwZqhol4gQyt+rHZtYGmIOv16oX+CioUYlIs/Te8uV0io/nqLZtOT8+nm7DhikxSkQKpEPOb/wPnzCzN4EU59xnwQ1LRJqbVatWsWHjRvrEx9N+1y7eBMYlJIQ6LJGgaGgSgHoHK5nZYOdcdnBCEpHm6L777iMmJobRHTpA27acqwH/EsEaqhynNrDPAWce4VhEpJlavnw5ry5cSN/4eD4qLOS45GQN+JeI1tAkAGc0ZSAi0nxNe+ghyqqquCM2lg8rK1nQtSu333uv7jdKxApknOMv69quSQBEosenOTkkAlWVlawD+g8apMQoES2Q3qpDajxOxDfmMRtNAiASFTZu3MjadevoGBvL21VVbIqJYbRW3pAIF0hv1RtqPjez1sAzgZzczH4GPAzEAv9wzt1fz3EXAy8BQ5xzHwdybhFpGvPmzcM5x6g2bTi2dWv6JiSoI45EvEAmAahtL9CnsYPMLBbf6h3n4Zt67jIzO7aO41oBNwKeQ4hFRIKoqqqKxx9/nPZHHcW2khLK8/NZBfRXRxyJcIHcc1yMr3cq+JLpscD8AM59ErDGOZfrP88LwGjgy1rH3QNMAf4QYMwi0kSysrLYvHkz1x99NOe1bMkir1fzqEpUCOSe40M1HlcA651zeQG8rguwscbzPOCA/1FmNgjo6px7zczqTY5mdh1wHUBaWhpZWVkBfHupzev16todgmi+bn/7299ISkriuN/9jrL4eE4C2nfrFvD1iOZrdzh03UIvkHuOSwH886rG+R+3c87taOSlVsc2t3+nWQwwDRgfQAyzgdkAGRkZbsSIEY29ROqQlZWFrt3Bi9brtmPHDt5//32OO/ZYvrnzTkrN2OAcrW68kb9NmRLQOaL12h0uXbfQa/Seo5ldZ2Zbgc+Aj/HNrxpIp5k8oGuN5+nA5hrPWwHHAVlmtg44GVhkZicGFrqIBNNzzz3Hvn376Hz00eRWVdE1NpaNsbHff8IViWCBNKveAvzYOVdwkOdeAfQxs57AJmAscHn1TufcLiC1+rmZZQF/UG9VkdBzzjFnzhz69euHy81lT0wML5eVsbtzZy1PJVEhkN6q3+HroXpQnHMVwG+Bt4CvgPnOuS/M7G4zG3Ww5xORpvPxxx+Tk5PD4AEDuNg5burSheS2bRl6ySXqjCNRIZDK8U/AB2bmAfZVb3TO3djYC51zbwBv1Np2Rz3HjgggFhFpApmZmbRs2ZKhp57K3EWLuALIBU7V4H+JEoEkx1nAO0AOUBXccEQk1IqLi3n++ee55JJLaBkXx/DUVArLyxkeH6/B/xI1AkmOFc65m4IeiYg0Cy+99BJ79uyhY/v2rM3L44OCAq4Ard8oUSWQ5Piuf5zhYg5sVm1sKIeIhKH//d//pXWLFgzJyuKR7dsZkZoKKSlav1GiSiDJsbqH6Z9qbHPAMUc+HBEJpc8//5ycnBzGpqVxUlkZI8xYBQxo0YIPWrTQ+o0SNQKZBKBnUwQiIqE3e/Zs4uLi2OD1sqS4mOXAqZMnU5SezvWDB6unqkQNrecoIgDs3buXZ555hsEDB3LWtm37O+F0S0/nmokTQx2eSJMKZJzjkBpfw4G/AhqnKBJhFixYwM6dO+mTkcG7+fm02bGD5QUFmDrhSBQK6nqOIhI+Zs+eTdeuXVm/dCnDzdhsxvDUVHXCkagUtPUcRSR8fPHFF7z//vv0OeYYTi0rIy8ujq4xMVq7UaJWMNdzFJEwMWfOHOLj44nbtIlvvF66VlYyp00brrvjDnXCkagUzPUcRSQMlJSUMG/ePPr06sWFxcX06t2bfxUVcdbYsfz6178OdXgiIVFvcjSz3kDH6vUca2wfbmYtnHPfBT06EQm6BQsWUFRUREZyMm8VFjKssJCCbt24SqtvSBRr6J7jdGBPHdtL/PtEJALMnDmTDh06cH3Llozv1ImvkpPpPXKkmlMlqjWUHHs45z6rvdG/3mKPoEUkIk3m008/5cMPP2TEiBE8tWEDmzduZM3u3fTS6hsS5Rq655jYwL6WRzoQEWl6M2fOJDExkS4dOtDROYpjYxngHLmrV4c6NJGQaqhyXGFm19beaGa/BlYGLyQRaQq7d+/m2WefZezYsSS2bMk6oAuwju+7p4tEq4Yqx8nAv8zsCr5PhicCCcCFwQ5MRILr2Wefpbi4mIkTJ5KTk8NbwJLKSjbFxDBazaoS5epNjs65rcCpZnYGcJx/8+vOuXeaJDIRCRrnHI8//jiDBw9myJAhvLpgAWe3bk3n5GT6JiRoVhyJeoFMH/cu8G4TxCIiTeS9997jiy++4B//+AcfffQRny1eTJzXS5rXy6pu3bhZs+JIlDuU6eNEJMzNnDmT1q1bM3bsWHKys7nYOQ3jEKlByVEkymzbto0FCxZw1VVXkZSUhCUk8OT69RrGIVKDkqNIlMnMzKS8vJwJEyYAkLt6NQM0jEPkAEqOIlGkoqKCxx9/nLPOOot+/foBvmEb69AwDpGalBxFosgrr7xCXl4eN9zw/TKtvTIy2IR/GIf/uUi0C2RVDhGJEI8++ig9evTgggsu2L/NlZVxfseOJJWX0zU+XsM4RFDlKBI1PvvsM5YtW8akSZOIjY3dv90SElheUED7XbtYXlCAJSSEMEqR5kGVo0iUePTRR2nZsiW/+tWvDtjuysoYnppKYXk5w1U5igCqHEWiQmFhIc899xzjxo2jXbt2B+xbm5fHu/n5tNmxQ5WjiJ+So0gUyMzMpKSk5ICOOAAej4flzzzDcDM2mzE8NVWVowhKjiIRr7Kykscff5wRI0bQv3//A/YtWriQU/btIy8ujq4xMawC+mvqOBElR5FIt3jxYtavX19n1fjZ4sV84/XSpqKCOa1a8fM77tDUcSKoQ45IxJs2bRrdu3dn1KhRB2xftHAh5xcX06t3b/5VVMRZY8fy61//OkRRijQvqhxFItjHH3/MsmXL+N3vfkdc3PefhaurxrcKC/l0zRoKkpIYPWZMCCMVaV5UOYpEsGnTptGqVasfVITVK3G06dSJRV6vVuIQqUWVo0iEysvLY/78+Vx77bWkpKQcsE8rcYg0TMlRJEI99thjVFVV/aAjDmglDpHGKDmKRCCv18usWbMYM2YMPXr0OGCfx+Ph008+YR3QMyGBjbGxWolDpBbdcxSJQHPnzmXnzp3cdNNNB2z3eDzMmjSJYcXFLADeSEykol07dcYRqUXJUSTCVFZWMn36dE455RROPvnkA/blZGfzE6+XM+LjKenYkY2nnsrtkyerM45ILUqOIhFm0aJFfPfdd9x///0/2Lc2L49la9dSasby2FjGnX22EqNIHXTPUSSCOOd44IEHOOaYY/j5z39+wD7NoyoSOCVHkQiybNkyPB4Pf/jDHw4Y9A+aR1XkYCg5ikSQBx54gA4dOjB+/PgDtmseVZGDo+QoEiE+/fRT/v3vf/O73/2Oli1bHrCvekac8Z06Udy6NWf98peaR1WkAUqOIhFiypQpJCcnM3HixB/ss4QE5m7cyPb8fHK9Xs2II9II9VYViQBr167lxRdfZPLkybRt2/YH+11ZGcNTUyksL2d4fLw64og0QpWjSASYOnUqMTEx/P73v69z/9q8PN7Nz6fNjh0sLyjAEhKaOEKR8KLkKBLmtm/fzpNPPsmVV15Jly5dfrBfQzhEDp6So0iYmzZtGqWlpdxyyy117s/JzmaEGevN6OychnCIBEDJUSSM7dixg0cffZRLLrmEH/3oR3UeYwkJLC8ooLtzPFxVxfHjxmkIh0gj1CFHJIxNnz4dr9fL7bffXu8x1Z1xksrLuTg+nm7p6U0YoUh4UuUoEqZ27tzJww8/zJgxYzjuuOPqPa66cmy/a5c644gESJWjSJh65JFH2L17d4NVo8fj4b2332Z4aiqkpHBuRYU644gEQJWjSBjavXs306ZNY/To0QwcOLDOYzIzM7n3kkvo+sEHLC8oIL68nG+TktQZRyQAqhxFwtBjjz3Gzp07+ctf/lLnfo/Hw9N33cUvCws5MT4eUlPJGjSICVq7USQgSo4iYWbPnj1MnTqV888/nxNOOKHOY6pX4FhmBuXlrAItaixyENSsKhJmHn30UXbs2NFg1agVOEQOj5KjSBgpKipiypQpjBw5ss5k5/F4eGL6dC52jlsHDaKiQwetwCFyCNSsKhJGHnzwQXbt2sW99977g30ej4dZkybRt7iYuRs3Mh4o7dCB0WPGNHmcIuFOyVEkTOTn5/Pwww9z2WWXcfzxx+/f7vF4yMnOZu3atZwPjOndG4CsgQPVAUfkEAU1OZrZz4CHgVjgH865+2vtvwm4BqgAtgO/cs6tD2ZMIuHqvvvuY9++fdx11137t2VmZvLK3XdzccuWfGbGtwBbtvBtUpISo8hhCFpyNLNYYAZwDpAHrDCzRc65L2sctgo40Tm318wmAlOAS4MVk0i4Wr9+PU888QRXX301ffr0AQ4crjEgPh46dWL1qFEU9ezJ9YMHKzGKHIZgVo4nAWucc7kAZvYCMBrYnxydc+/WOP6/wLggxiMStu6++24A7rjjDuD7jjcjzPggPh7Ky3m1pITbx4xRUhQ5Asw5F5wTm10M/Mw5d43/+ZXAUOfcb+s5/jEg3zn3g54GZnYdcB1AWlraCfPnzw9KzJHO6/WSnJwc6jDCTqiv24YNG7j66qu56KKLmDRpEsXFxRRs2ECLqip2l5WRHBfHHqD90UeTmpoasjjrEuprF6503Q7dGWecsdI5d+LhnieYlaPVsa3OTGxm44ATgZ/Utd85NxuYDZCRkeFGjBhxhEKMLllZWejaHbxQX7fRo0eTlJTEjBkzWLt2LU/PmsWIVasY26cP961Zw4fNuONNqK9duNJ1C71gJsc8oGuN5+nA5toHmdnZwG3AT5xz+4IYj0jYycrKYtGiRfz9739n7dq1BwzVAPg2ObnZJkaRcBbMSQBWAH3MrKeZJQBjgUU1DzCzQcAsYJRzblsQYxEJO1VVVdx8881069aNyZMnk5OdzU+8Xq6Ii2O4f67U62fMUGIUCYKgVY7OuQoz+y3wFr6hHE86574ws7uBj51zi4AHgWTgJTMD2OCcGxWsmETCybPPPkt2djbPPvssLVu2ZG1eHsvWrqXUjOWxsYw7+2wlRpEgCeo4R+fcG8AbtbbdUePx2cH8/iLhau/evdx2222ceOKJXHbZZXg8HpY/8wzDzdhsxvDUVK3LKBJEmltVpBmaNm0aeXl5TJ06lRUrVuwftpGfmEjXmBhWgdZlFAkiTR8n0szk5+dz//338/Of/5wWLVrs74SzpKCA4ampzHOOq7TKhkhQKTmKNDP/7//9P0pLSzlx4EAWLVx4wHyp3wwcyIPqnSoSdEqOIs3Ie++9x7x58/hR+/a0z8zknZgYvm3ZUvOlijQxJUeRZqKiooJJkyaR0qoVV+7dy+iKCrzl5Xx+4YUUDRum+VJFmpCSo0gz8fjjj/PZZ58xYvhwPvjvf+kVE8P7/P/27j88qurO4/j7S0IgmURF+Q2yoFQX0YoRCWqLgPVnFESplW67UKJlXdw+Yn22dpW2aPuouKjrFouUKNoHWrXIDxWsoqAuNWND/IGtIhSwkB8kQBQyISSTnP1jbnAYA0wmhDtJPq/nuc/cmTn35nsON3xz7o9zYEjfvtx8661+hyfSoehuVZEksHPnTmbOnElOTg6B8nL2derEC7W17O3RQ5MVi/hAyVEkCUydOpVQKMSws8/mB5mZzB4+nFN69ybnxht1KlXEB0qOIj6bN28eK1euZFyPHhSvW8ezoRA7amqo6dlTvUYRn+iao4hPgsEg6wsKmDlrFr3S0vjdsGG8UlHB+rFjNWGxiM+UHEV8EAwGmTNlCl8UF7Nn3z5GnXoqr1RU8DIwTRMWi/hO5jo7gAAAEltJREFUyVHEByuWLKHvtm28UFPDiJQUBo0apcc1RJKIkqOID7bv2MHymhpOMKN7Whp99biGSFLRDTkix1kwGOTd115jLzCkUyfdeCOShNRzFGllwWCQDUVFnOOdMn3lpZf4++7djD3lFM7MyOAkPa4hknSUHEVaUTAY5Inp08kFnsjPp/aRR1j07LOkpqTw3X79WNe5s3qNIklIyVGkFW0oKuKSqiou6NyZqro6fnXvvWzatImHHnoIFwjoBhyRJKXkKNJKgsEgb69bx7bPPqOuUyd+XV/Ph59+yrRp07jzzjv9Dk9EjkDJUaQVNJ5ODZSXc65zVPTowbbiYnr26MGcOXP8Dk9EjkLJUaQVNJ5OPalzZxY4xzsVFewNh3n0wQcJBAJ+hyciR6FHOURagaWlsXD7dirKyviwoYHC6mp+OG0aU6ZM8Ts0EYmDkqPIMRYMBvm/1av5ZvfufNanDyX19Zw1ZAhz5871OzQRiZNOq4ocQ43XGs8IhXi5ooJPw2HSUlOZ8/DDpKSk+B2eiMRJyVHkGGh80H/r1q3kAhNOP51nduxgV00N8+fP58orr/Q7RBFpBiVHkRZqnGFjbChEQadObEpPZ+n77/NJVRUzZszg5ptv9jtEEWkmXXMUaaEVS5Yw4h//YPzevVxVUUFlnz4sLinh0ksv1WMbIm2Ueo4iLeSAt+rr6dfQwEv19byzbh3nnnsuy5Ytw8z8Dk9EEqCeo0gLnX7mmRQDy8Nh3qmrIz09nRdffJHMzEy/QxORBCk5irSQq63l8h49eMeMBjP+49Zb6d+/v99hiUgL6LSqSAu51FSeLC2lvL6eoWlpDBw82O+QRKSFlBxFmil6fsYRI0bw7OLFlNfXc32vXgzv1g1XW+t3iCLSQkqOIs2Qn5/PsnvvZWJ6Ok9kZvLbYcN4fe1ahvbsyaR+/VhpxtjsbL/DFJEWUnIUiVMwGOTpWbP41927ObdzZ/6Ulsbv169n6tSp3HLLLXz03nuan1GknVByFInTiiVLuPDAAd4yY93+/fz+iy8YM2YM8+fPJyUlhZEjR/odoogcI0qOInEIBoN8+OKLpFZVsa+2ltfDYYYOHcqqVas0ZqpIO6TkKHIYjTfe9BkwgM0ff8yUQICC3r357y1bGDRwIO+++y5dunTxO0wRaQVKjiKe6LtQAZ6YPp1cYFdeHpaWxgNlZRQWF3PqiSfy9DPPkJGR4W/AItJqlBxF+HKqqVzgifx8+owdSy5wQ58+LHWOp5YupbC4mOHZ2Tz62GNcfPHFfocsIq1IyVEE2FBUxCVVVVzQuTNVdXVsBF4GwiUlPL5oEW/9+c9MnjyZ/Px8XWMU6QCUHEUAS0tj4fbt1AGLgO+deSZjL7+cH912Gx8XFXHPPfcwa9YsOnXSiIsiHYGSo3RI0dcXc3JyIuOjDhgAqalcHg5TXlLCnDlz2LxlCz/5yU+47777/A5ZRI4jJUfpcGKvLzJ3LudkZ/NEIMAZwPMHDvDuww9jZrz66qt+hysiPtA5IukwgsEgC37zG1YsWXLwZptcItcbc3JyuPmxx3i6Vy9e27yZfv36UVBQwOjRo32OWkT8oJ6jtBuxp0pjv2vsLS4PhdgEUFrKy8C07GzKy8v5+c9/zurVq5k8eTKPP/64HtUQ6cCUHKVdaOpUaXSC3FBUdLC3SGkp68eOpXLQIKZlZ1NRUcGECROorKxkwYIFTJ06FTPzrS4i4j+dVpU2LxgMMu/RR7mkquqQU6XRzsnO5mVgiddbHH/DDUycNIl58+Zx7bXX0r17dwoKCsjLy1NiFBH1HKVta+wxnhEKsXD7dgDezMxkWsy0UTk5OTB3LhuKipiWnc2ePXs4++yzKSsr45577mHmzJmkpaX5UQURSUJKjtLmRF9bPHi6dPBgANYOG8a/3X57k9NG5eTk0LNnT2bMmMHy5cs566yzWLZsGcOHDz/ONRCRZKfTqtKmNPYUu+Xn88T06Vha2sHTpZsCgcMmxurqan72s58xZMgQVq9ezf33309RUZESo4g0ST1HaVNWLFnCiPJyxg4YADU1VNbWMi3qdGlsYqytreWpp57ivvvuo7i4mEmTJjF79mz69+/vUw1EpC1QcpQ2IxgMUvDcc3xQVsbe8nIKBw3ix15CjE2K4XCYRYsWMWvWLLZu3cqFF17I4sWLGTVqlE/Ri0hbotOq0masWLKEqyoquDUtjfUNDWScf/5XkmIoFGLu3LkMGTKEKVOm0K1bN1auXMm6deuUGEUkbuo5StJrvAGnuKSEvwGDUlKo7dyZ0/v2PVimuLiYuXPnMm/ePCorK8nJyeGhhx5i/PjxejRDRJpNyVGSWvTD/btDIap79OD1hgbCgQBXjxvHCy+8wJNPPsmqVasAmDBhAnfccQcXXXSRv4GLSJum5ChJqbG3uHXr1kNGtim85hoaGhrI2LGD66+/noqKCvr27ctdd91FXl4ep512mt+hi0g7oOQoSSd2HNRPGhr4cO9env/8cyoWLmTXrl106dKF3Nxc8vLyuPzyy0lN1aEsIseO/kcRXxxpkPAP16/nvJoadppRVlrKB9XVvFBXR5cuXbj66qv59re/TW5uLieccIJP0YtIe6fkKMfckRJf4/fRg4QfePhhunbtyjvvvMPbb7/NmjVr2LNnDwCZaWmMGzeOyZMnM2bMGDIzM49zbUSkI1JylGPqSLNjOOfYuXMnzy1ejJWXsxJYs2sXz4wZQ31DAwADBw4kNzeXfv36kdm1K5decQUjR470r0Ii0iEpOUrcjtYjBPigsJCRtbX0ycrixNJS7v/lLznx5JPZvHkzGzduZPfu3QfLZqWkkBUI8L0bb2TcuHGMGDFCI9eISFJo1eRoZlcC/wOkAAuccw/EfN8FeAY4H9gNfMc5t601Y5Lmcc4RCoVYs2YN+XfdxXnhMIvDYS6YOJFAIEBZWRmlpaWUlZVRXFxMSUkJ9fX1X+5g61b69+/P4MGDuf766xk6dChDhw7lwIEDlHz2GV9v4kF+ERG/tVpyNLMUYC5wGbAD+IuZrXDO/S2qWB5Q6ZwbbGY3AQ8C32mtmPzmnDtkaWhoaPJ9Q0MDDQ0N1NfXf+W1cQmHwwdfo5e6ujpqa2upq6ujrq6OAwcOUFtby8cff0z5zp0sXbqUXr16sX//fmpqagiFQoRCIaqrqwmFQuzdu5d9+/YdXD7//HPC4fDBOiz3XtfMng3AKaecQp8+fejduzejR49mwIABhMNh9u/bx8hvfIPrrruO9PR0H1pbRCRxrdlzHAFsds5tATCzPwDjgejkOB74hbf+R+DXZmbOOXe4nW7atCnhmzKOsNsmv4/3fVOvjUs8P9cPZkZ6ejqBQICMjAwyMjIIBAJkZWUxcOBAsrKyyMrKolu3bpx00klUVlZS8NxzjE5N5b3UVG6ZPZvLLrtMcyCKSLvUmsmxH7A96v0OIPb82cEyzrmwmX0BnALsii5kZj8Efui9PRAKhT5qlYjbKYMevaF7PVgKuDLY5ZyrqK6uprq6ujm7CrwJGQ6ql19zTai14k1C3Yk5JiVuarvEqN0Sd+ax2ElrJsemBrSM7ULFUwbn3HxgPoCZFTrnNAlfAtR2iVG7JU5tlxi1W+LMrPBY7Kc1Z+XYAZwa9b4/UHK4MmaWCpwI7GnFmERERI6qNZPjX4CvmdkgM0sDbgJWxJRZAUz21icCbxzpeqOIiMjx0GqnVb1riLcBfyLyKMeTzrm/mtm9QKFzbgWQD/zOzDYT6THeFMeu57dWzB2A2i4xarfEqe0So3ZL3DFpO1NHTURE5FCteVpVRESkTVJyFBERiZFUydHMrjSzjWa22czuauL7Lmb2rPd90MwGRn33U+/zjWZ2xfGM22+JtpuZDTSz/Wb2vrfMO96x+y2OthtlZkVmFjaziTHfTTazTd4yOXbb9qyF7VYfdczF3qTX7sXRdneY2d/M7EMze93M/inqOx1zibVb84+52CHN/FqI3LTzd+A0IA34ADgrpsy/A/O89ZuAZ731s7zyXYBB3n5S/K5TG2i3gcBHftchydtuIPB1ImMAT4z6/GRgi/fazVvv5nedkr3dvO+q/K5DkrfdGCDDW7816vdVx1wC7ea9b/Yxl0w9x4PDzTnnaoHG4eaijQee9tb/CFxqZuZ9/gfn3AHn3FZgs7e/jqAl7dbRHbXtnHPbnHMfAg0x214BvOac2+OcqwReA648HkEngZa0W0cXT9utcc41Dl1VQOQZcdAxl2i7JSSZkmNTw831O1wZ51wYaBxuLp5t26uWtBvAIDN7z8zeNLNvtnawSaYlx42OuS81t+5dzazQzArM7LpjG1rSa27b5QGrEty2PWlJu0ECx1wyzefYkuHm4hqGrp1qSbuVAgOcc7vN7HxgmZkNdc7tPdZBJqmWHDc65g7VnLoPcM6VmNlpwBtmtsE59/djFFuyi7vtzOx7wHDgkuZu2w61pN0ggWMumXqOLRluLp5t26uE2807Db0bwDm3nsg5/TNaPeLk0ZLjRsfcl5pVd+dcife6BVgLnHcsg0tycbWdmX0LuBsY55w70Jxt26mWtFtix5zfF1qjLpimErnAPIgvL7gOjSkznUNvLHnOWx/KoTfkbKHj3JDTknbr0dhORC50FwMn+12nZGq7qLIL+eoNOVuJ3BjRzVvvEG3XwnbrBnTx1rsDm4i5saI9L3H+vp5H5A/Vr8V8rmMusXZL6JjzvdIxlbga+NSr4N3eZ/cS+SsAoCvwPJEbbt4FTova9m5vu43AVX7XpS20G3AD8FfvQCsCrvW7LknYdhcQ+as1BOwG/hq17VSvTTcDP/C7Lm2h3YCLgA3eMbcByPO7LknYdquBncD73rJCx1zi7ZboMafh40RERGIk0zVHERGRpKDkKCIiEkPJUUREJIaSo4iISAwlRxERkRhKjiIxokbw/8jMnjezjFb8Wfd6Dy5jZrc392dZxBtmdoL3/kdm9rGZLToGsU0xs75R7xeY2VkJ7us2M/tBS2MSOV70KIdIDDOrcs5leuuLgPXOuYfj2M6I/E4lNNi2mW0DhjvndjVjm1zgW865Gd77T4g857s1plyqi4yr25x41gJ3OucKm7PdYfaVAaxzznWk0XCkDVPPUeTI3gYGw8H54j7yltu9zwZ6PbXHiQykcKqZTTKzDV65B71yKWa20Ptsg5k1JrOFZjbRzH4E9AXWmNkaM8szs0cagzCzW8ysqQT9L8Byr8w8IiMdrTCzGWb2CzObb2avAs94sb7tzbNYZGYXRe3/P724PjCzBywyB+NwYJHXi043s7VmNtwr/5U6ep9XmdmvvP0UmFkvABeZLWGbmXWU2XKkrfN71AMtWpJtwZv7jciQVcuJzA13PpHRNQJAJpGRhc4jMm9hAzDS26Yv8A8iQ/OlAm8A13nbvxb1M07yXhfiDa8GbAO6e+sBIiOBdPbe/xk4p4lYPwOyot5H7+MXwHog3XufAXT11r8GFHrrV3n7b5wL72TvdS2RnizR7w9XR6+MwxtpCZgN3BO1/d3Aj/3+99WiJZ5FPUeRr0o3s/eBQiJJIB/4BrDUORdyzlUBLwCNU3x95pwr8NYvANY65ypc5DTmImAUkXEhTzOz/zWzK4EjznzinAsRSTrXmNk/E0mSG5ooerJzbt8RdrXCObffW+8M/NbMNhAZTrDx+uG3gKecNxeec27PkWI7Qh0BaoGXvPX1RP54aFROJLGKJL1kmrJKJFnsd84Ni/7gKJNDh6KLNlXAOVdpZucSmbB2OnAjkXEyj2QB8F/AJ8BThykTNrNO7vDXOaNjm0Fk7MlziVxSqYmKuTk3HxypLeqcc437qufQ/2O6Avu/uolI8lHPUSQ+bwHXmVmGmQWACUSuR8YKApeYWXczSwEmAW+aWXegk3NuCTATyG5i231AVuMb51yQyDQ93wV+f5i4NhK5zhiPE4FSL5F+H0jxPn8VmNp4p6yZndxUPEerYxw//wzgozhjFfGVkqNIHJxzRUSuD75LJDkscM6910S5UuCnwBq82U6cc8uJzFq+1jtdu9ArE2s+sMrM1kR99hyRuzwrDxPay8DoOKvxODDZzAqIJKqQF/MrwAqg0IvvTq/8QmBe4w05cdTxaC4mMnOCSNLToxwiSczMXgIecc69fpjv+wDPOOcuO76RNY+ZnQfc4Zz7vt+xiMRDPUeRJGRmJ5nZp0SufzaZGOFgL+63jYMAJLHuRE4ni7QJ6jmKiIjEUM9RREQkhpKjiIhIDCVHERGRGEqOIiIiMZQcRUREYvw/NuPayhq8+dsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "values = np.linspace(fmin,fmax,100) \n",
    "fit_mean, fit_stdev = norm.fit(X,loc = average, scale = std) # fit MLE of the distribution \n",
    "cumul_p = norm.cdf(values, loc = fit_mean, scale = fit_stdev)\n",
    "\n",
    "# plot the cumulative probabilities vs. the sorted porosity values\n",
    "plt.subplot(122)\n",
    "plt.scatter(sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7,label='data')\n",
    "plt.plot(values,cumul_p, c = 'black',label='fit'); plt.legend(loc='upper left')\n",
    "plt.xlabel(feature + ' ' + funits); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric and Fit Gaussian CDFs')\n",
    "plt.ylim([0,1]); plt.xlim([0,0.25])\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of univariate statistics in Python.\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at [Python Demos](https://github.com/GeostatsGuy/PythonNumericalDemos) and a Python package for data analytics and geostatistics at [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy). \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "rise": {
   "theme": "white",
   "transition": "fade"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}