PythonDataBasics_Univariate_Statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Data Analytics \n",
    "\n",
    "### Basic Univariate Statistics in Python \n",
    "\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Analytics: Basic Univariate Statistics\n",
    "\n",
    "Here's a demonstration of calculation of univariate statistics in Python. This demonstration is part of the resources that I include for my courses in Spatial / Subsurface Data Analytics and Geostatistics at the Cockrell School of Engineering and Jackson School of Goesciences at the University of Texas at Austin.  \n",
    "\n",
    "We will cover the following statistics:\n",
    "\n",
    "#### Measures of Centrality\n",
    "* Arithmetic Average / Mean\n",
    "* Median\n",
    "* Mode (most frequent binned)\n",
    "* Geometric Mean\n",
    "* Harmonic Mean\n",
    "* Power Law Average\n",
    "\n",
    "#### Measures of Dispersion\n",
    "* Population Variance\n",
    "* Sample Variance\n",
    "* Population Standard Deviation\n",
    "* Sample Standard Deviation\n",
    "* Range\n",
    "* Percentile w. Tail Assumptions\n",
    "* Interquartile Range\n",
    "\n",
    "#### Tukey Outlier Test\n",
    "* Lower Quartile/P25\n",
    "* Upper Quartile/P75\n",
    "* Interquartile Range\n",
    "* Lower Fence\n",
    "* Upper Fence\n",
    "* Calculating Outliers\n",
    "\n",
    "#### Measures of Shape\n",
    "* Skew\n",
    "* Excess Kurtosis\n",
    "* Pearson' Mode Skewness\n",
    "* Quartile Skew Coefficient\n",
    "\n",
    "#### Nonparmetric Cumulative Distribution Functions (CDFs)\n",
    "* plotting a nonparametric CDF\n",
    "* fitting a parametric distribution and plotting\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).   \n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "You will need to copy the data file to your working directory.  The dataset is available on my GitHub account in my GeoDataSets repository at:\n",
    "\n",
    "* Tabular data - [2D_MV_200wells.csv](https://github.com/GeostatsGuy/GeoDataSets/blob/master/2D_MV_200wells.csv)\n",
    "\n",
    "#### Importing Packages\n",
    "\n",
    "We will need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "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' is available at https://github.com/GeostatsGuy/GeoDataSets.  It is a comma delimited file with X and Y coordinates,facies 1 and 2 (1 is sandstone and 2 interbedded sand and mudstone), porosity (fraction), permeability (mDarcy) and acoustic impedance (kg/m2s*10^6). We load it with the pandas 'read_csv' function into a data frame we called 'df' and then preview it by printing a slice and by utilizing the 'head' DataFrame member function (with a nice and clean format, see below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>facies_threshold_0.3</th>\n",
       "      <th>porosity</th>\n",
       "      <th>permeability</th>\n",
       "      <th>acoustic_impedance</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.1700</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.2750</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.2970</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.0480</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.1230</td>\n",
       "      <td>3.979</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3375</td>\n",
       "      <td>2525</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1239</td>\n",
       "      <td>1.4680</td>\n",
       "      <td>2.337</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2295</td>\n",
       "      <td>1325</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1790</td>\n",
       "      <td>31.9330</td>\n",
       "      <td>3.491</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3715</td>\n",
       "      <td>3045</td>\n",
       "      <td>2</td>\n",
       "      <td>0.1914</td>\n",
       "      <td>116.7810</td>\n",
       "      <td>2.187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2865</td>\n",
       "      <td>2215</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1763</td>\n",
       "      <td>3.0030</td>\n",
       "      <td>2.048</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>55</td>\n",
       "      <td>1545</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1674</td>\n",
       "      <td>5.2130</td>\n",
       "      <td>2.251</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>2125</td>\n",
       "      <td>1105</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1369</td>\n",
       "      <td>3.6930</td>\n",
       "      <td>3.627</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>565</td>\n",
       "      <td>195</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1095</td>\n",
       "      <td>0.2627</td>\n",
       "      <td>2.860</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>3475</td>\n",
       "      <td>1035</td>\n",
       "      <td>1</td>\n",
       "      <td>0.1677</td>\n",
       "      <td>9.9140</td>\n",
       "      <td>2.742</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       X     Y  facies_threshold_0.3  porosity  permeability  \\\n",
       "0    565  1485                     1    0.1184        6.1700   \n",
       "1   2585  1185                     1    0.1566        6.2750   \n",
       "2   2065  2865                     2    0.1920       92.2970   \n",
       "3   3575  2655                     1    0.1621        9.0480   \n",
       "4   1835    35                     1    0.1766        7.1230   \n",
       "5   3375  2525                     1    0.1239        1.4680   \n",
       "6   2295  1325                     1    0.1790       31.9330   \n",
       "7   3715  3045                     2    0.1914      116.7810   \n",
       "8   2865  2215                     1    0.1763        3.0030   \n",
       "9     55  1545                     1    0.1674        5.2130   \n",
       "10  2125  1105                     1    0.1369        3.6930   \n",
       "11   565   195                     1    0.1095        0.2627   \n",
       "12  3475  1035                     1    0.1677        9.9140   \n",
       "\n",
       "    acoustic_impedance  \n",
       "0                2.009  \n",
       "1                2.864  \n",
       "2                3.524  \n",
       "3                2.157  \n",
       "4                3.979  \n",
       "5                2.337  \n",
       "6                3.491  \n",
       "7                2.187  \n",
       "8                2.048  \n",
       "9                2.251  \n",
       "10               3.627  \n",
       "11               2.860  \n",
       "12               2.742  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#df = pd.read_csv(\"2D_MV_200wells.csv\")                      # read a .csv file in as a DataFrame\n",
    "df = pd.read_csv(\"https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/2D_MV_200wells.csv\") # read a .csv file in as a DataFrame\n",
    "#print(df.iloc[0:5,:])                                       # display first 4 samples in the table as a preview\n",
    "df.head(n=13)                                                   # we could also use this command for a table preview "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "por = df['porosity'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the output from the above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "por is a <class 'numpy.ndarray'>, with shape of (200,)\n"
     ]
    }
   ],
   "source": [
    "print('por is a ' + str(type(por)) + ', with shape of ' + str(por.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's go through all the univariate statistics listed above one-by-one.\n",
    "\n",
    "#### Measures of Central Tendency\n",
    "\n",
    "Let's start with 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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity average is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_average = np.average(por)\n",
    "print('Porosity average is ' + str(round(por_average,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Median\n",
    "\n",
    "\\begin{equation}\n",
    "P50_x = F^{-1}_{x}(0.50)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity median is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_median = np.median(por)\n",
    "print('Porosity median is ' + str(round(por_median,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity mode is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_mode = stats.mode(np.round(por,2))\n",
    "print('Porosity mode is ' + str(round(por_mode,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity geometric mean is 0.15.\n"
     ]
    }
   ],
   "source": [
    "por_geometric = scipy.stats.mstats.gmean(por)\n",
    "print('Porosity geometric mean is ' + str(round(por_geometric,2)) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity harmonic mean is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_hmean = scipy.stats.mstats.hmean(por)\n",
    "print('Porosity harmonic mean is ' + str(round(por_hmean,2)) + '.')"
   ]
  },
  {
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity law mean for p = 1.0 is 0.15.\n"
     ]
    }
   ],
   "source": [
    "p = 1.0\n",
    "por_power = np.average(np.power(por,p))**(1/p)\n",
    "print('Porosity law mean for p = ' + str(p) + ' is ' + str(round(por_power,2)) + '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(200,)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 1.0\n",
    "np.power(por,p).shape"
   ]
  },
  {
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity population variance is 0.0011.\n"
     ]
    }
   ],
   "source": [
    "por_varp = stats.pvariance(por)\n",
    "print('Porosity population variance is ' + str(round(por_varp,4)) + '.')"
   ]
  },
  {
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0011.\n"
     ]
    }
   ],
   "source": [
    "por_var = stats.variance(por)\n",
    "print('Porosity sample variance is ' + str(round(por_var,4)) + '.')"
   ]
  },
  {
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample standard deviation is 0.0329.\n"
     ]
    }
   ],
   "source": [
    "por_stdp = stats.pstdev(por)\n",
    "print('Porosity sample standard deviation is ' + str(round(por_stdp,4)) + '.')"
   ]
  },
  {
   "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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity sample variance is 0.0329.\n"
     ]
    }
   ],
   "source": [
    "por_std = stats.stdev(por)\n",
    "print('Porosity sample variance is ' + str(round(por_std,4)) + '.')"
   ]
  },
  {
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity range is 0.17.\n"
     ]
    }
   ],
   "source": [
    "por_range = np.max(por) - np.min(por)\n",
    "print('Porosity range is ' + str(round(por_range,2)) + '.')"
   ]
  },
  {
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 13th percentile is 0.11.\n"
     ]
    }
   ],
   "source": [
    "p_value = 13\n",
    "por_percentile = np.percentile(por,p_value)\n",
    "print('Porosity ' + str(int(p_value)) + 'th percentile is ' + str(round(por_percentile,2)) + '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(por)"
   ]
  },
  {
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity interquartile range is 0.04.\n"
     ]
    }
   ],
   "source": [
    "por_iqr = scipy.stats.iqr(por)\n",
    "print('Porosity interquartile range is ' + str(round(por_iqr,2)) + '.')"
   ]
  },
  {
   "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": 21,
   "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(por, [25, 75])\n",
    "lower_fence = p25 - por_iqr * 1.5\n",
    "upper_fence = p75 + por_iqr * 1.5\n",
    "por_outliers = por[np.where((por > upper_fence) | (por < lower_fence))[0]]\n",
    "print('Porosity outliers by Tukey test include ' + str(por_outliers) + '.')\n",
    "por_outliers_indices = np.where((por > upper_fence) | (por < lower_fence))[0]\n",
    "print('Porosity outlier indices by Tukey test are ' + str(por_outliers_indices) + '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.06726, 0.05   , 0.06092])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "por[np.where((por > upper_fence) | (por < lower_fence))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06092"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "por[198]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Measures of Shape\n",
    "\n",
    "##### Pearson's Mode Skewness\n",
    "\n",
    "\\begin{equation}\n",
    "skew = \\frac{3 (\\overline{x} - mode_x)}{\\sigma_x}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity skew is 0.28.\n"
     ]
    }
   ],
   "source": [
    "por_skew = (por_average - por_mode)/por_std\n",
    "print('Porosity skew is ' + str(round(por_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": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity 3rd central moment is -1.22e-05.\n"
     ]
    }
   ],
   "source": [
    "por_cm = scipy.stats.moment(por,moment=3)\n",
    "print('Porosity 3rd central moment is ' + str(round(por_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": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity quartile skew coefficient is 0.14.\n"
     ]
    }
   ],
   "source": [
    "por_qs = ((np.percentile(por,75)-np.percentile(por,50))\n",
    "          -(np.percentile(por,50)-np.percentile(por,25))) /((np.percentile(por,75))-np.percentile(por,25))\n",
    "print('Porosity quartile skew coefficient is ' + str(round(por_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": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7AUlEQVR4nO3de3yU9Zn38c9FyIAmAmqCiJCWgqYnRVChalGg1q2llLrwaLei0kfrYbFddu22drePXdrd7W673dV22aptXMRu16WNB6xude2KUq13lUFNtY3SoCRgGhIRmAjmdD1/zD04jDlMEiZz+r5fL17JzNz3nWt+DLm4fvfvYO6OiIiIvG1UtgMQERHJNUqOIiIiKZQcRUREUig5ioiIpFByFBERSaHkKCIikkLJUSTPmNktZvb/sh3HYJnZX5nZD7Mdh0g6lBxlxJnZK2b2BzMrS3ruSjPbmMWwss7MVpjZLwc6zt2vcfdvDOH6r5jZfjOLhe3/72ZWPrRoB8/d/97drwxjebeZuZmNHur1zOx4M6sxs9fMbJ+Z/c7MVic+V+H128P322ZmvzCzi1OusdHMDoTHJP6cObx3KoVAyVGyZTTwZ9kOoj/D+cWdKWZWMsxLLHb3cmA2cAbw1UH+/JxoEzM7BvgVcARwprsfBXwUmABMTzp0Zvh+q4G1wL+a2ddSLnedu5cn/flVxt+A5DwlR8mWbwNfNLMJvb1oZmeZ2dNmtif8elbSaxvN7Btm9kRYMTxsZhXha4mK5Coz2xlWFdcnnTvHzH5lZm+Er/2rmUWSXnczW2lmLwMvh8/dbGaNZrbXzDab2byk4//GzH5iZj8KY6kzs5PM7Ctm1hKed37S8eOTqp0dZva3ZlZiZu8DbgHODKuXN8Lj15rZ983sQTNrBxaEz/1t0jWXmNmzYXy/N7OPDdT47r4D+G/gg+E1PmlmL4TtsjGMJ3H9V8zsy2b2PNBuZqMHOP7L4XvbZ2b1ZvaRpLb6UXjY4+HXN8L3e66ZvW5mJyddZ2JY6Vb28hb+AtgHLHf3V8L31Ojuf+buz/fyflvd/U7gWuArZnbsQG0kxU3JUbLlGWAj8MXUF8Kq4AHgu8CxwD8DD6T8QvsM8FlgIhDp5ToLgBOB84EbzOy88Plu4M+BCuBM4CPAn6ac+ylgLvD+8PHTwKnAMcCPgZ+Y2dik4xcDdwJHA1uAh4j/2zoB+Dpwa9KxdwBdwAxgVhjfle7+W+Aa4Fdh9TIh5b3+HXAUcEi3q5nNAdYBf0m8ajoHeIUBmNlU4OPAFjM7CfhPYBVQCTwI3J/8nwbgT4BF4c94T1/Hm1k1cB1wRljN/VEf8ZwTfp0Qvt/HgLuA5Sk/8xF339XL+ecBd7t7z0DvNcV9xHst5gzyPCkySo6STTcCn++lMlgEvOzud7p7l7v/J/A74kko4d/d/SV33w+sJ568kq1293Z3rwP+nfgvWtx9s7s/FV73FeKJ69yUc7/p7q+H18bdf+TubeE53wHGEO+mS9jk7g+5exfwE+IJ4x/cvZP4L/x3m9kEMzsOuABYFcbWAvwL8OkB2uk+d3/C3Xvc/UDKa1cAt7v7/4Sv73D33/VzrXvDqvSXwGPA3wMXAw+E1+gE/ol4d+VZSed9N6zM9g9wfHfYPu83s1J3f8Xdfz/A+0u4A/iMmSV+L11K/D8dvTkWeC3N6x4UxttK/D86Cd8NK+A3zCw62GtKYVJylKxx998APwNuSHlpMvBqynOvEq/EEpqTvn8TSB1Y0phy7mSAsMvzZ2bWbGZ7iSeHin7OxcyuN7Pfhl28bwDjU875Q9L3+4FWd+9OekwY37uAUuC1xC9j4sl5Iv1r7Oe1qUC6yQfgU+4+wd3f5e5/Gia7Q9o7rMYaObS9k2Po83h330q8ovwboMXM7jKzyekE5u4B0A6ca2bvJV5db+jj8Dbg+HSum8zMSon/5+X1pKe/ELbJBHefPdhrSmFScpRs+xrwOQ79RbyTeCJJVgXsGMR1p6acuzP8/vvEq9AT3X0c8FeApZx7cKua8P7il4GLgKPD7s49vZyTjkbgLaAi6ZfxOHf/QOrP7SuePq45vZ/X03FIe5uZEW+/5Pb2dI939x+7+4fDYxz4x15+Zl/v6Q7iXauXAj/tpVJOeAS4MKnKTNcS4t3avx7keVJklBwlq8JK47+ALyQ9/SBwkpl9Jhz8cTHx+38/G8Sl/5+ZHWlmHyB+b/K/wuePAvYCsbA6uXaA6xxF/JfpLmC0md0IjBtEHAe5+2vAw8B3zGycmY0ys+lmlujW/QMwJeVe30BqgM+a2UfC650Qvq/BWA8sCq9RClxPPIk/OdjjzazazBaa2RjgAPHKubuXa+wCeojfv0x2J3Ah8QS5rp+Y/5n438MdZvYugPC9/7OZnZJ6sJkdY2aXAGuAf3T3tn6uLaLkKDnh68DBOY/hL65PEP+l2wZ8CfiEu7cO4pqPAVuBXwD/5O4Ph89/kfgAl33AD3g7afblIeKjOl8i3pV4gP67OQdyGfEBRC8Cu4Gf8nb34P8CLwDNZpbWe3X3XxNP/v9CvKJ9jHdW3QNdo554Mvoe8ftxi4lP+egYwvFjgH8In28m3mX8V71c403ig4yeCLuYPxQ+3wREiVeWm/qJ+XXi9zg7gcDM9hH/u95D/O894Tkzi4XPXQn8ubvfmEazSJEzbXYshcTM3g1sA0rDATKSZ8zsdmCnuw9qDqbI4ZQTE3pFRODgf27+mPg0F5GsyVi3qpndbvFJ0L/p43Uzs++a2VYze97MNEpMpIiZ2TeA3wDfdvdt2Y5HilvGulXN7BwgBqxz9w/28vrHgc8Tn4g8F7jZ3edmJBgREZFByFjl6O6Pc+hcolRLiCdOd/engAlmNuh5SyIiIodbNu85nsCho/6awufeseqFmV0FXAUwduzY06qqqkYkwELT09PDqFEaoDxYarehU9sNTabbraen5+DPGOjndHV1YV1dlJjR7Y6PHs3o0bk7XOWll15qdffe1uMdlGy+w94mUffax+vutwG3AVRXV3t9fX0m4ypYGzduZP78+dkOI++o3YZObTc0mWy3IAi4deVKFhFfwPjqNWuYO7fvO1qDPT7bzCx1da0hyWZybOLQVUym8PYqJiIikgF10SjnxmKcUVpKrLOTumi032Q3d+5cWLOGumiUq2fPzunEeDhlMzluAK4zs7uID8jZE64gIiIiGWKRCLe/+ipvdndzV0kJl0UGXpBp7ty5RZMUEzKWHM3sP4H5QIWZNRFfQ7MUwN1vIb5E2MeJr1zxJvFVPkREJIMa6uuZ6U57SQkz3WnQbapeZSw5uvufDPC6Aysz9fNFROSdnPgGm2cS33H6HfPsBNDaqiIiRWV6dTU7gIe7u9kRPpZ3yt3xuCIictg11Ndz3vjxTC4v56RIBO/odX35oqfKUUSkSARBwPP3389LsRidzc1sAU6erZU7e6PKUUQkDwVBQF00ysmDmF5RF42yoqyMKbNmccf27cxYvLjoRqGmS5WjiEieCYKA76xYQdc3v8l3VqwgCIK0zjt59mweAJoOHODAxIksWbo0s4HmMVWOIiJ5ZkNtLXO2b2dJaSmxtjY21NamVQEW64T+oVByFBHJUX11nTrwBDDNnSeA9w3imsU4oX8o1K0qIpKDEmuaHl1Tw60rVx7Sdbpk6VK6qqr4xfjxdFVVqXs0A1Q5iojkgOQqEeKDZxYBS48/Hl577ZA1UOfOnctX166lLhrlcnWPZoQqRxGRLEutEtvb2w8Onql97TUe4J1TLubOncuV116rxJghqhxFRLJsQ20tc1paWFhVBQcOsP/NNzV4JsuUHEVEsqimpobH1q3jmDfeINbWxq+rqrj8yCMBDZ7JJiVHEZEsCYKAO1avZsW+fYwZPZqNZWXMWLyYsrKybIdW9HTPUUQkS+qiUeab8WhPD292ddEydqxGnuYIJUcRkSyxSIRNra28y52be3o4ZflydaPmCHWriohkSUN9PXPKyzmuvJxLIxEqp0zJdkgSUuUoIpIF2iEjt6lyFBHJgg21tSxqb2f6jBncs3u3dsjIMUqOIiIZ1Nv6qImq8cW2Ns5ua6O1qorLNRAnp6hbVUQkQ/paHzWxr+INs2bRcOyxqhpzkJKjiEiG1EWjnBuLccZbb3FuLEZdNApoX8V8oG5VEZEMsUiEtY2NdAL/ASyPRADtq5gPlBxFRDLEOzqYV1FBW2cn80pL8Y6Og69pabjcpm5VEZEM2dbUxKPNzUx4/XU2tbZiYeUouU/JUUQkA4IgYNOddzLPjJ1mzKuoOKRylNym5CgikgGJdVNfNWOyuyb55xklRxGRDNC6qflNA3JERA6zIAj45SOPMK+igqpx47i0q0vrpuYZVY4iIodRYuL/Sc8+y6bWVko7O3m5rExdqnlGlaOIyGG0obaWOS0tXFxVBcDGU0/lmlWr1KWaZ5QcRUQOgyAIuK+2lmD9el5sayPW1saWqiquV2LMS0qOIiLDlOhKHdvSwgVtbczUTht5T/ccRUSGIQgCbrnpJs6Nxbi8qoongN/v3q01U/OcKkcRkSE6OPimvZ21jY2sALqqqmhcvJirly5V1ZjHlBxFRIYgUTHOj8W4/MQTgfjgm6/qHmNBUHIUERmk1IoR4OXyco1KLSBKjiIig6TpGoVPyVFEZABBEFAXjR6cyP/8/fdrukaBU3IUEelHogt1EXBrTQ3HL1zIirIypsyaxR3bt2u6RoHSVA4RkX4kulAXjh3LIsCBB4CmAwc0XaOAqXIUEelDEASHdKH+uqqK65cuhaVLqYtGuXr2bFWNBUrJUUSkD3XRKMvcmTBpEhtisUO6UJUUC5u6VUVE+mCRCGsbG9nV3ExDLMb06upshyQjRJWjiEgfvKODeRUVtHV2Mq+0FO/oyHZIMkJUOYqI9GFbUxOPNjcz4fXX2dTaikUi2Q5JRoiSo4hIL4IgYNOddzLPjJ1mzKuoUOVYRJQcRUR6UReNMt+MV82Y7M4WOLgIgBQ+JUcRkV5YJMKm1lbe5c7NPT2csny5RqgWEQ3IERHpRUN9PXPKyzmuvJxLIxEqp0zJdkgyglQ5ioikSEz+fykWo7O5WV2qRUiVo4hIig21tSxqb2f6jBncs3u31k8tQkqOIiJJampqeGzdOo554w3OaWujtaqKy7V+atFRchQRCQVBwB2rV7Ni3z7GjB7NxrIyVY1FSvccRURCiekbj/b08GZXFy1jx2rXjSKl5CgiEtL0DUlQt6qISCixlmpZZyfLSkup0vSNoqXKUUQkpLVUJUHJUUQEraUqh1JyFBFBa6nKoZQcRUTQYBw5VEaTo5l9zMzqzWyrmd3Qy+vjzex+M3vOzF4ws89mMh4Rkb4cXEv1hBO4dPp0pmkwTlHLWHI0sxJgDXAB8H7gT8zs/SmHrQRedPeZwHzgO2amO+AiMqK0lqqkyuRUjjnAVndvADCzu4AlwItJxzhwlJkZUA68DnRlMCYRkXeoi0ZZ5s6ESZPYEItpVRzB3D0zFzZbBnzM3a8MH18KzHX365KOOQrYALwXOAq42N0f6OVaVwFXAVRWVp62fv36jMRc6GKxGOXl5dkOI++o3YYuX9pu586d7Gtu5ljgdTOOmTqVioqKrMWTL+2WixYsWLDZ3U8f7nUyWTlaL8+lZuI/Ap4FFgLTgf8xs03uvveQk9xvA24DqK6u9vnz5x/2YIvBxo0bUdsNntpt6PKh7YIg4JYbbuDDLS3sNMMnTmTUDTewbNmyrMWUD+1W6DI5IKcJmJr0eAqwM+WYzwJ3e9xWYBvxKlJEZERoCof0JpPJ8WngRDObFg6y+TTxLtRk24GPAJjZcUA10JDBmEREDqEpHNKbjHWrunuXmV0HPASUALe7+wtmdk34+i3AN4C1ZlZHvBv2y+7emqmYRERSHZzCUV7OpZEIlZrCIWR44XF3fxB4MOW5W5K+3wmcn8kYRET6kpjCMToWozIWY0tVFderS1XQrhwiUsQ0hUP6ouXjRKRoWSTC7a++ys7GRrbu3cv06upshyQ5QslRRIpWQ309M91pLylhpjsN9fXZDklyhJKjiBQtB14BTgi/ZmZJFMlHSo4iUrSmV1ezA3i4u5sd4WMR0IAcESliDfX1nDd+PJPLyzkpEtHmxnKQKkcRKUraiUP6o8pRRIqSpnFIf5QcRaSoBEFAXTTKtqYmnmxs5BLia1aepfuNkkTJUUSKRhAE3LpyJYuAdS0tzK+ogHHjOL+rS/cb5RC65ygiRaMuGuXcWIwz3nqL+WZsAY4eM4aXy8p0v1EOocpRRIqGRSKsbWykE9gEnLVqFbunTOHq2bN1v1EOoeQoIkUhCAJ++cgjzEvqSq2cMoUrr70226FJDlK3qogUvMS9xpOefZZNra2UdnaqK1X6pcpRRApaEATcctNNzI/FuPzEEwHYeOqpXLNqlbpSpU9KjiJSsGpqarj3619nFrC2Nb6P+svl5UqMMiAlRxEpSEEQcMfq1VzW1sbppaVQUcHGWbOUGCUtSo4iUpDqolHmm/FoTw/dBw6wZfx4vqrEKGnSgBwRKUgWibCptZV3uXNzTw+nLF+uxChpU+UoIgUhsSzcyeGcRe/o4PyqKipHj+bScNqGSLqUHEUk7wVBwHdWrGBhezvfKSvj+rVrOXn2bG4tK+Mk4MkxY1ioaRsyCEqOIpL3NtTWMmf7dpaUlhJra2NDbS1/961vwZo11EWjWgFHBk3JUUTyngOPd3dzQk8Pj7vzwfD5uXPnKinKkGhAjojkvenV1ewAHu7uZkf4WGQ4VDmKSN7zjg4WHXccZZ2dTC0t1fZTMmyqHEUk7yWmbRy7Zw+bWluxSCTbIUmeU+UoInmvob6eOeXltJeXc34kospRhk3JUUTyVhAE3FdbS7B+PeNiMSpjMbZUVXG9pm3IMCk5ikheSmxDNbalhQva2pg5Ywb37N7NjMWLNUJVhk3JUUTyUl00yrmxGBNKS6np7qZ8924OTJzIkqVLsx2aFAANyBGRvGSRCGsbG9nV3MwOIDjnHK5es0ZVoxwWqhxFJC8lD8JZFolQefbZSoxy2KhyFJG8EwQBz99/Py/FYnQ2N7MFOFmDcOQwUuUoInlnQ20ti9rbma5BOJIhqhxFJK8EQcBT69fzs+ZmttTX01pWpkE4ctipchSRvLKhtpYLdu3ifZEI6zo6OPK001Q1ymGnylFE8kpiB4693d3sHzWKyZMnZzskKUBKjiKSV7QDh4wEdauKSF5pqK/nvPHjmVxezklaR1UyRJWjiOQNTeGQkaLKUUTyRl00yjJ3JkyaxIZYTFM4JGNUOYpI3rBIhNtffZWdjY1s3btX9xslY5QcRSRvNNTXM9Od9pISZrrTUF+f7ZCkQA2YHM3sOjM7eiSCERHpz46dO2no6WFqSQmNJSV4tgOSgpVO5TgJeNrM1pvZx8zMMh2UiEiqIAho27yZfaNGcXdHB3srK7UyjmTMgMnR3b8KnAjUACuAl83s781seoZjExE5KLGe6lerqzl20iTmXnSRBuNIxqR1z9HdHWgO/3QBRwM/NbNvZTA2ERHg7SkcD7W18dzWrVpPVTJuwKkcZvYF4HKgFfgh8Jfu3mlmo4CXgS9lNkQRKXbahUNGWjrzHCuAP3b3V5OfdPceM/tEZsISEYlXjPfV1hKsX8+LbW2c3dZGa1UVl6tqlAxLJzlOS02MZnanu1/q7r/NUFwiUuSCIODWlSsZ29LCBW1tzFTVKCMoneT4geQHZlYCnJaZcERE4uqiUc6NxZhQWkpNdzflu3dzYOJE3WuUEdHngBwz+4qZ7QNOMbO94Z99QAtw34hFKCJFySIR1jY2squ5mR1AcM45XL1mjapGGRF9Vo7u/k3gm2b2TXf/ygjGJCJCQ309c8rLaS8vZ1kkQuXZZysxyojpMzma2Xvd/XfAT8zsHcveu3s0o5GJSNGqqanhsXXrOGbPHipjMbZUVXG9dt+QEdTfPcfrgc8B3+nlNQcWZiQiESlqQRBwx+rVrNi3jzGjR7OxrEyDcGTE9det+rnw64KRC0dEil1dNMp8Mx7t6eEcd1rGjuUqDcKREdZft+of93eiu999+MMRkWJnkQibWls5052be3q4cPlyVY0y4vrrVl3cz2sOKDmKyGHnHR3Mq6igrLOTZaWlVE2Zku2QpAj116362ZEMREQE3q4cLwF+DiyPRLIdkhSh/rpVl7v7j8zsL3p73d3/OXNhiUixSp7CcX4kgnd0ZDskKUL97cpRFn49qo8/Awr3f6w3s61mdkMfx8w3s2fN7AUze2wQsYtIgQmCgKfWr+eF3bt5a8cOtgAnawqHZEF/3aq3hl9XD+XC4TJza4CPAk3EN0ze4O4vJh0zAfg34GPuvt3MJg7lZ4lIYdhQW8sFu3bxvkiEdR0dHHnaaRqMI1kx4H6OZvYeM7vfzHaZWYuZ3Wdm70nj2nOAre7e4O4dwF3AkpRjPgPc7e7bAdy9ZbBvQEQKx46dO9nY2UkM6CgtZfLkydkOSYqUxfcx7ucAs6eIV4D/GT71aeDz7t7vf+fMbBnxivDK8PGlwFx3vy7pmJuAUuKLmx8F3Ozu63q51lXAVQCVlZWnrV+/Pq03J4eKxWKUl5dnO4y8o3YbusG0XXt7O6+98grdb71FKdBVWsoJ73kPZWVlA55baPSZG7oFCxZsdvfTh3uddHblMHe/M+nxj8zsuj6PTjqvl+dSM/Fo4jt8fAQ4AviVmT3l7i8dcpL7bcBtANXV1T5//vw0fryk2rhxI2q7wVO7Dd1g2u6H3/8+pd/7HpWdnWyIxTju8su56vOfz2yAOUqfuezrb7TqMeG3j4aDae4intwuBh5I49pNwNSkx1OAnb0c0+ru7UC7mT0OzAReQkSKSmIXjkuABuCs6upshyRFrL/KcTPxZJioAK9Oes2Bbwxw7aeBE81sGrCDeHfsZ1KOuQ/4VzMbDUSAucC/pBe6iBQSTeGQXNLfaNVpw7mwu3eF3a8PASXA7e7+gpldE75+i7v/1sx+DjwP9AA/dPffDOfnikj+SUzhOGL3bo7ds4ct06ZpFw7JqnTuOWJmHwTeD4xNPNfbwJlU7v4g8GDKc7ekPP428O104hCRwqQpHJJrBkyOZvY1YD7x5PggcAHwS2DA5Cgikg4HHu/u5gQz9o8axXRN4ZAsG3CeI7CM+GjS5nC91ZnAmIxGJSJFZXp1NTuAh7u72RE+FsmmdLpV97t7j5l1mdk4oAVIZxEAEZG0NNTXc9748UwuL+ckDcaRHJBO5fhMuMzbD4iPYI0Cv85kUCJSPIIg4Pn77+elWIzO5matpyo5YcDK0d3/NPz2lnBk6Th3fz6zYYlIsdhQW8ui9namz5jBPbt3M2PxYg3GkaxLd7TqHwMfJn7f/JfEp16IiAxLYgrHc83NnNPSQuu0aVy+dGm2wxJJa7TqvwEzeHtt1avN7Dx3X5nRyESk4GkKh+SqdCrHc4EPerhCuZndAdRlNCoRKWhBEFAXjbJj505eBKaVlNBRWqopHJIz0kmO9UAV8Gr4eCrqVhWRIQqCgFtXrmQR0NbezpuVlfyip4eusjKWqEtVckR/C4/fT/we43jgt2aWGKE6B3hyBGITkQK0obaWOS0tLKyqAmDz4sW8Z9o0Lp89W12qkjP6qxz/acSiEJGikJi28WJbG7G2Nn5dVcX1S5cqKUrO6W/h8ccS35vZccAZ4cNfu3tLpgMTkcISBAG33HQTy9x576xZ3LF9u6ZtSM4acBEAM7uI+KT//wNcBARmtizTgYlI4aipqeFvL7qIqU8+ydrGRn7X2sqBiRN1j1FyVjoDcv4aOCNRLZpZJfAI8NNMBiYihSEIAu5YvZrL2to4vbQUKirYOGsW16xapapRclY6yXFUSjdqG+ktOyciQl00ynwzHu3pofvAAbaMH89XlRglx6WTHH9uZg/x9iIAF5OyR6OISF8sEmFTaytnunNzTw8XLl+uxCg5r9/kaGYGfJf4YJwPAwbc5u73jEBsIlIAvKODeRUVlHV2sqy0lKopU7IdksiA+k2O7u5mdq+7nwbcPUIxiUgBsUiEx1pa+HR3Nw+WlHBZJJLtkEQGlM69w6fM7IyBDxMReaeG+npmutNeUsJMdxrq67MdksiA0kmOC4gnyN+b2fNmVmdmWj5ORNKyY+dOGnp6mFpSQmNJCZ7tgETSkM6AnAsyHoWIFKT29nbaNm9m36hR3N3Rwd7JkzW3UfJCf2urTgT+ivh2VXXAN91970gFJiL5b88bb8Q3Mq6u5p7du5l+0UUaqSp5ob9u1XVAO/A9oJz4qFURkbQEQcCbb7zBQ21tPLd1K63adUPySH/dqpPc/a/D7x8ys+hIBCQihaEuGuVoYMWkSWyIxbSOquSV/pKjmdnRxOc2ApQkP3b31zMdnIjkL4tEaN2zh57GRraWlHBWdXW2QxJJW3/JcTywmbeTI0CienTgPZkKSkTyX0N9Pe+dNIk/aAqH5KH+tqx69wjGISIFZsfOnUybNImpJSU84c77sh2QyCBoAXEROeyCIKBt82a6IT6Fo7JSg3Ekr6Qzz1FEZFDqolGWuTMhEqH86KM1hUPyjipHETnsLBJhbWMjnZ2dNMRiTNdgHMkzaSVHM/uwmX02/L7SzKZlNiwRyWeJnTi6SkqYV1GBd3RkOySRQRkwOZrZ14AvA18JnyoFfpTJoEQkv21rauLR5mZGd3WxqbUV004ckmfSqRwvBD5JfLUc3H0ncFQmgxKR/BUEAZvuvJN5ZnSAKkfJS+kkxw53d+JzGzGzssyGJCL5bENtLWe+9RZNo0cTMWMLcPLs2dkOS2RQ0kmO683sVmCCmX0OeAT4QWbDEpF8FAQBz99/Py/FYkzo6qK1pIRP3XijRqpK3hlwKoe7/5OZfRTYC1QDN7r7/2Q8MhHJOwencITrqR517LFccskl2Q5LZNAGTI5m9ufAT5QQRWQgiSkclwANwAVjxmQ7JJEhSWcRgHHEd+V4HbgL+Km7/yGzYYlIPmqor2dOeTnt5eWcH4mAe7ZDEhmSAe85uvtqd/8AsBKYDDxmZo9kPDIRySs1NTU8tm4dv92zh87mZrYARxx5ZLbDEhmSwSwf1wI0A23AxMyEIyL5KAgC7li9mhX79jFm9Gg2lpUxY/Fiyso0uF3yUzqLAFxrZhuBXwAVwOfc/ZRMByYi+aMuGmW+GY/29PBmVxctY8dqoXHJa+lM5XgXsMrdP+DuX3P3FzMdlIjkF4tE2NTayrvcubmnh1OWL9f0DclrfXarmtk4d98LfCt8fEzy6+7+eoZjE5E8kVhLtayzk2WlpVRNmZLtkESGpb97jj8GPgFsJr46jiW95sB7MhiXiOSRbU1NPN7czHIzfl5SwnKtpSp5rs/k6O6fCL9qBw4R6VPyWqo7zbSWqhSEdAbk/CKd50SkOCUG47xqxmR3raUqBaG/e45jgSOBCjM7mre7VccRn+8oInJwMM6Z4WCcCzUYRwpAf/ccrwZWEU+Em3k7Oe4F1mQ2LBHJB0EQ8MtHHmFeRQVV48ZxaVcXlRqMIwWgv3uONwM3m9nn3f17IxiTiOSBIAi4deVKTmpv5+HWVlYccQQvl5ezUF2qUgDS2ZXje2b2QeD9wNik59dlMjARyW0bamuZ09LCxVVVAGw89VSuWbVKXapSENLZleNrwHziyfFB4ALgl4CSo0gRCoKA+2prCdav58W2NmJtbWypquJ6JUYpIOmsrboMmAlscffPmtlxwA8zG5aI5KJEV+rYlhYuaGtj5owZ3LN7NzMWL1ZilIKSTnLc7+49ZtZlZuOIL0CuBQBEilBdNMq5sRgTSkup6e6mfPduDkycqHVUpeCkkxyfMbMJwA+Ij1qNAb/OZFAikpuSNzPeAQTnnMM1n/+8qkYpOOkMyPnT8NtbzOznwDh3fz6zYYlILkqsodrW2cmi0lKqzj5biVEKUn+LAPQ5HtvMZrt7NDMhiUiuSkz4vwT4OWgNVSlY/VWO3+nnNQcWHuZYRCTHNdTXM6e8nPbycs6PRLSGqhSs/hYBWDCSgYhIbguCgKfWr+eI3bs5ds8etkybxvWa8C8FKp15jpf19rwWARApLhtqa7lg1y7eF4mwrqODI087TfcbpWClM1r1jKTvxwIfAaJoEQCRouLA493dnGDG/lGjmD5Z+w9I4UpntOrnkx+b2XjgzoxFJCI5aXp1NQ8BD3d3s2PUKJZUV2c7JJGMGXA/x168CZyYzoFm9jEzqzezrWZ2Qz/HnWFm3Wa2bAjxiEiGJXbfWHTccZxdXc2yadM0GEcKWjr3HO8n3qMC8WT6fmB9GueVEN/a6qNAE/C0mW1w9xd7Oe4fgYcGF7qIjATtviHFKJ17jv+U9H0X8Kq7N6Vx3hxgq7s3AJjZXcAS4MWU4z4P1HLovU0RyRGJJeMWlJayv6KCjbNmafcNKXjm7gMfBYTrqh5Mpu7++gDHLwM+5u5Xho8vBea6+3VJx5wA/Jj4nMka4Gfu/tNernUVcBVAZWXlaevXD1i4Si9isRjl5eXZDiPvFHu7tba20rZ9O8cArwPHVlVRUVGR1rnF3nZDpXYbugULFmx299OHe510ulWvAr4B7Ad6ACPezTrQ4uPWy3Opmfgm4Mvu3m3W2+HhSe63AbcBVFdX+/z58wcKW3qxceNG1HaDV+zt9tdf+hL7vvtdqszY7s5RX/gCf/etb6V1brG33VCp3bIvnW7VvwQ+4O6tg7x2EzA16fEUYGfKMacDd4WJsQL4uJl1ufu9g/xZIpIhO3bupLWnh7MjEZ5w533ZDkhkBKSTHH9PfITqYD0NnGhm04gv4P9p4DPJB7j7tMT3ZraWeLfqvUP4WSKSAUEQ0LZ5M/tGjeLujg72Tp6s7amkKKSTHL8CPGlmAfBW4kl3/0J/J7l7l5ldR3wUaglwu7u/YGbXhK/fMvSwRWQk1EWjLHNnwgknsCEWY/pFF2kgjhSFdJLjrcD/AnXE7zmmzd0fBB5Mea7XpOjuKwZzbRHJvOT9GxuAszTxX4pEOsmxy93/IuORiEjOSd6/cV5pqSb+S9FIJzk+Go5YvZ9Du1X7ncohIvkrCALqolG2NTXxpPZvlCKUTnJMDKL5StJz6UzlEJE8lFgRZxGwrqWF+RUVMG4c53d1qXKUopHOwuPTBjpGRApHYkWcM0pLmW/GFmDmmDE8OWaMloyToqH9HEXkEIlBOJ3AJuCsVavYPWUKV8+erZGqUjS0n6OIHCJ1EE7VlClcee212Q5LZERpP0cROcS2piYeb25muRk/LynRIBwpShndz1FE8ksQBGy6807mmbHTjHkVFRqEI0UpY/s5ikj+qYtGmW/G7804x52fAYs0CEeKUCb3cxSRPGORCJtaWznTnZt7erhw+XINwpGi1GdyNLMZwHHu/ljK8/PMbIy7/z7j0YnIiGqor2dOeTnHlZdzaSRC5ZQp2Q5JJCv6u+d4E7Cvl+f3h6+JSAEJgoDn77+fl2IxOpub2QKcrC5VKVL9dau+292fT33S3Z8xs3dnLiQRyYaDO3BMmsSGWIwZixerS1WKVn/JcWw/rx1xuAMRkeza1tTE49u2sdyMhpIS7cAhRa2/btWnzexzqU+a2RXA5syFJCIjTVM4RA7VX+W4CrjHzC7h7WR4OhABLsxwXCIygjSFQ+RQfSZHd/8DcJaZLQA+GD79gLv/74hEJiIjRlM4RA6VzvJxjwKPjkAsIpIlmsIhcqihLB8nIgVEUzhE3imdFXJEpIBpCofIO6lyFClyFolw+6uvsrOxka179zJdUzhElBxFil1DfT0z3WkvKWGmOw319dkOSSTrlBxFipwDrwAnhF+9v4NFioSSo0iRm15dzQ7g4e5udoSPRYqdBuSIFDnv6GDRccdR1tnJ1NJSrYwjgipHkaKXWADg2D172NTaikUi2Q5JJOtUOYoUOe/oYF5FBW2dncxT5SgCqHIUKXrbmpp4tLmZCa+/rspRJKTkKFLEtBuHSO+UHEWK2IbaWs586y2aRo9m6qhRWjpOJKTkKFKkktdUndDVxQ+OOopP3Xijlo4TQQNyRIpW6pqqH7nsMq644opshyWSE1Q5ihQpi0RY29jIruZmGmIxTf4XSaLKUaRIJfZwbC8v5/xIRANxRJKochQpQkEQ8NT69bywezdv7dihgTgiKVQ5ihShDbW1XLBrF++LRFjX0cGRp52mgTgiSZQcRYpMEAQ89+yzjOruZtqYMXSUljJ98uRshyWSU5QcRYpIEATcunIlZ7e381PgwbFj6TrmGJYsXZrt0ERyipKjSBGpi0Y5NxZjQWkp+487jsazzuKrq1apS1UkhZKjSBHZ1tTE49u2ccCMTSUlLD/vPCVGkV5otKpIkdA6qiLpU3IUKRJaR1UkfUqOIkVA66iKDI7uOYoUAa2jKjI4qhxFioDWURUZHFWOIkXAOzqYV1FBW2cn80pLNRBHZACqHEWKgEUiPNbSwrjWVh5racEikWyHJJLTlBxFikBDfT0z3WkvKWGmOw319dkOSSSnKTmKFAEHXgFOCL96NoMRyQNKjiJFYHp1NTuAh7u72RE+FpG+aUCOSBFoqK/nvPHjmVxezkna2FhkQKocRQpc8gIAnc3NWhlHJA1KjiIFLAgCbrnpJpa5c8OsWTQceywzFi/WyjgiA1C3qkiBqqmp4d6vf51ZwNrWVlYAByZO1N6NImlQchQpQEEQcMfq1VzW1sbppaVQUcHGWbO4Rns3iqRFyVGkACV24HjcDDo72QLa1FhkEHTPUaTAaAcOkeFTchQpIKkDcLomTtQOHCJDoG5VkQIRBAG3rlzJSe3trG1s1AAckWFQchTJc0EQUBeNsm3bNhYBS2fMAGDjqadqAI7IEGU0OZrZx4CbgRLgh+7+DymvXwJ8OXwYA6519+cyGZNIIUlM11h2xBE8b8bLAK+9xstlZUqMIsOQseRoZiXAGuCjQBPwtJltcPcXkw7bBpzr7rvN7ALgNkD/mkXSkDxdY2ZpKUyaRP0nP8nuadO4evZsJUaRYchk5TgH2OruDQBmdhewBDiYHN39yaTjnwKmZDAekYKRGHgz34wnS0uhs5P79u/nq0uXKimKHAbmnpnNa8xsGfAxd78yfHwpMNfdr+vj+C8C700cn/LaVcBVAJWVlaetX78+IzEXulgsRnl5ebbDyDu51m7t7e20bt/OmJ4e9nZ0UD56NPuAY48/noqKimyHd4hca7t8oXYbugULFmx299OHe51MVo7Wy3O9ZmIzWwBcAXy4t9fd/TbiXa5UV1f7/PnzD1OIxWXjxo2o7QYv19rtr7/0JabedRefrKrilrY2fpXDA29yre3yhdot+zKZHJuAqUmPpwA7Uw8ys1OAHwIXuHtbBuMRyXuJCf4vtrURa2tjS1UV1+doYhTJZ5lcBOBp4EQzm2ZmEeDTwIbkA8ysCrgbuNTdX8pgLCIFYUNtLYva2/nTGTO0w4ZIBmWscnT3LjO7DniI+FSO2939BTO7Jnz9FuBG4Fjg38wMoOtw9BWLFKKamhoeW7eOY954g3Pa2mitquJyTfAXyYiMznN09weBB1OeuyXp+yuBdwzAEZFDJaZtrNi3jzGjR7OxrExVo0gGaW1VkRyXOm3jLXdaxo7VsnAiGaTl40RyWPJ6qQ+3tjKvooJ17lyuXTZEMkrJUSQH9bVe6kunnsq3NTpVJOOUHEVyTBAEfGfFCha2t/PUqFG8fMQRWi9VZIQpOYrkmA21tczZvp0lpaXEOjv5zYUXsvvss7VeqsgIUnIUyTE7du7kuc5OpprxBPC+yZO58tprsx2WSFHRaFWRHBIEAW2bN7Nv1Cju7uhgb2WlRqWKZIEqR5EckDwAZ0VZGVNOP507tm9n+kUXqStVJAuUHEWyLDFdYxFwX3s7LwMXAwcmTlTVKJIlSo4iWfKO6RrHHw+vvcbmhQu1YbFIlik5imRBX9M1HgCu1obFIlmn5CiSBZquIZLblBxFskDTNURym6ZyiIwwTdcQyX2qHEUyLDHw5uSwy7QuGtV0DZEcp8pRJIMS0zSOrqnh1pUrCYKAk2fP5gGg6cABTdcQyVGqHEUyqC4a5dxYjDPCgTd10Wj83uKaNdRFoxqAI5KjlBxFMmhbUxOPb9vGATPuKilheSQCwNy5c5UURXKYulVFMiQIAjbdeSfzzNhpxryKCryjI9thiUgalBxFMmRDbS1nvvUWTaNHM3XUKLYAJ8+ene2wRCQNSo4iGRAEAc/ffz8vxWJM6OriB0cdxaduvFFdqSJ5QslR5DALgoBbbrqJZe7cMGsWXRMn8pHLLuOKK67IdmgikiYNyBE5jBJTN05qb2dtYyMr0O4aIvlIyVHkMElUjPNjMS4/8UQANp56KtesWqXuVJE8o+QochjU1NRw79e/zixgbWsrAC+XlysxiuQpJUeRYQqCgDtWr+aytjZOLy2Figo2zpqlxCiSx5QcRYapLhplvhmP9vTQfeAAW8aP56tKjCJ5TaNVRYbJIhE2tbbyLndu7unhlOXLlRhF8pwqR5Fh8o4O5lVUUNbZybLSUqqmTMl2SCIyTKocRYYpUTkeu2cPm1pbsXD9VBHJX6ocRQYpdX9G7+jg/KoqGD2a87u6tH6qSAFQchQZhMSUjWVHHMGt5eWwZg0nz57NrWVlnAQ8OWYMC7V+qkjeU3IUSVPylI2ZpaUwaZL2ZxQpUEqOImlKnbLxs/37+WpYJWp/RpHCogE5ImnSlA2R4qHKUaQPiYE3x1dVARwceFM5ejSXdnVRqSkbIgVLyVEklDwKFeDWlStZBLRecQVBEGjgjUgRUXIU4e2tphYBt9bUcPzChSwClh5/PPeCBt6IFBklRxHiye/cWIwzSkuJdXZSDzwA8NprxOBgNamBNyLFQclRhPhgm7WNjXQC/wEsr65mydKlB+85KiGKFBclRylayfcYe1vlJlElbty4MduhisgI01QOKUqJe4xH19Rw68qVWCTCy2VlHD1mDC+XlR3sRhWR4qTKUYpGcqVYF40eHHDDa6+xu6ODqzXYRkRCSo5SUFIXBU9+Pnk06tnXXntwwM0DcDAhKimKCCg5SgFJTYCsWXMw2alSFJHBUHKUvJeoFrdt23ZIAqyLRg8mvZNnz44nTFWKIpIGJUfJa8nV4n3t7bwMhyTAhLlz52oCv4ikTclR8k5/A2s2L1zI7mnTek2AqhRFJF1KjpJXBhxYs3SpEqCIDJuSo+SVDbW1zGlpYWFVFRw4oIE1IpIRSo6SN4Ig4Kn163muuZm9LS08M20a12tgjYhkgJKj5I0NtbVcsGsX74tEWNfRwZGnnaakKCIZoeQoOS8xAGfHzp28CEwrKaGjtJTpkydnOzQRKVBKjpLTkgfgtLW382ZlJb/o6aGrrIwlS5dmOzwRKVBKjpKT+prYv3nxYt4zbRqXa/CNiGSQkqPknH4n9muqhoiMACVHyYq+FgiHd66D2t/EfhGRTFBylMOuv8SXeL2vBcKhl3VQVS2KyAhTcpTDaqDEB++sDJMXCAetgyoi2Tcq2wFI/giCgB9+//sEQdDnMcmJb1H4ONXJs2fzAFAbVoYnJy0QnjB37lyuvPZaJUYRyYqMJkcz+5iZ1ZvZVjO7oZfXzcy+G77+vJm987ek5IRERXh0TQ23rlzZZ4JMN/FdvWYNu6+4gqt7qSxFRLItY92qZlYCrAE+CjQBT5vZBnd/MemwC4ATwz9zge+HX+UwC4KA1l27CIJgSMlooK7QhHS7RLXkm4jkskxWjnOAre7e4O4dwF3AkpRjlgDrPO4pYIKZHZ/BmIpSouob3drab9XXn3QqwgR1iYpIvsvkgJwTgMakx028syrs7ZgTgNeSDzKzq4CrwodvmdlvDm+ohc2gchJUPLB5s5WAr/3Qh1oddg3hUmVr4UiHN//9Qx9qP9xx5rAKoDXbQeQptd3QqN2GrvpwXCSTydF6ec6HcAzufhtwG4CZPePupw8/vOKjthsatdvQqe2GRu02dGb2zOG4Tia7VZuAqUmPpwA7h3CMiIjIiMpkcnwaONHMpplZBPg0sCHlmA3AZeGo1Q8Be9z9tdQLiYiIjKSMdau6e5eZXQc8BJQAt7v7C2Z2Tfj6LcCDwMeBrcCbwGfTuPRtGQq5GKjthkbtNnRqu6FRuw3dYWk7c3/HLT4REZGiphVyREREUig5ioiIpMip5Dic5eYGOreQDbPdXjGzOjN79nANgc4nabTde83sV2b2lpl9cTDnFrJhtps+c/233SXhv9PnzexJM5uZ7rmFbJjtNvjPnLvnxB/ig3Z+D7wHiADPAe9POebjwH8Tnx/5ISBI99xC/TOcdgtfewWoyPb7yOG2mwicAfwd8MXBnFuof4bTbuFr+sz133ZnAUeH31+g33PDa7fw8aA/c7lUOQ5nubl0zi1UWqZv6AZsO3dvcfengc7BnlvAhtNuxS6dtnvS3XeHD58iPv87rXML2HDabUhyKTn2tZRcOsekc26hGk67QXxFoofNbHO4TF8xGc7nRp+5tw32vesz97aB2u4K4r0+Qzm3kAyn3WAIn7lc2ux4OMvNpbUMXYEa7jJ9Z7v7TjObCPyPmf3O3R8/rBHmruF8bvSZO9Rg3rs+c4fqte3MbAHxX/IfHuy5BWg47QZD+MzlUuU4nOXminkZumEt0+fuia8twD3Euy+KxXA+N/rMvW1Q712fuYHbzsxOAX4ILHH3tsGcW6CG025D+szlUnIcznJz6ZxbqIbcbmZWZmZHAZhZGXA+UEw7ngznc6PP3BDeuz5zA7edmVUBdwOXuvtLgzm3gA253Yb6mcuZblUfxnJzfZ2bhbcx4obTbsBxwD1mBvHPwo/d/ecj/BayJp22M7NJwDPAOKDHzFYRHyW3V5+5wbcb8a2Y9Jnr/9/rjcCxwL+F7dTl7qfr99zQ2o0h/p7T8nEiIiIpcqlbVUREJCcoOYqIiKRQchQREUmh5CgiIpJCyVFERCSFkqNICjPrDlfv/42Z/cTMjszQzzndzL4bfj/fzM4awjVWmdll4ffvDePeYmbThxnbqWb28aTHnxzqLhBmVmlmRTNdQwqDkqPIO+1391Pd/YNAB3BNOieZ2aDmDbv7M+7+hfDhfOK7CqQt/Hn/F/hx+NSngPvcfZa7/z7pODOzwf5bP5X43NhErBvc/R8GeY3EubuA18zs7KGcL5INSo4i/dsEzDCzY8zs3nCvuKfCZaows78xs9vM7GFgnZm9y8x+ER73i3DVDszs/4SV6HNm9nj43Hwz+5mZvZt4Av7zsPKbZ2bbzKw0PG6cxfejK02JbSEQDSdIfxxYBVxpZo+a2bvN7Ldm9m9AFJhqZt83s2fM7AUzW524iJmdYfH9754zs1+b2Xjg68DFYTwXm9kKM/vX8Pi+3uNai+8b+qSZNZjZsqRY7wUuOYx/LyIZpeQo0oewMrsAqANWA1vc/RTgr4B1SYeeRnwtx88A/0p8e7BTgP8AvhsecyPwR+4+E/hk8s9x91eAW4B/CSvWTcBGYFF4yKeBWndP3f7pbGBzeI0Hk66xIHy9Ooxllru/Cvx1uGLIKcC5ZnZKuBTXfwF/FsZ2HtAexvtfYTz/lfJz+3qPAMcTX/D5E0BypfkMMA+RPKHkKPJOR5jZs8R/oW8Haoj/wr8TwN3/Fzg2rLAANrj7/vD7M3m7m/NO3t4Z4AlgrZl9jvjyVwP5IW8v8/dZ4N97OeZ4YFc/13g13L8z4SIziwJbgA8QX86tGngt3HsRd9/r7l0DxNbXewS419173P1F4st2JbQAkwe4rkjOyJm1VUVyyH53PzX5CQsXZkyRWHuxvZ9rOYC7X2Nmc4lXg8+a2an9nIO7PxF2jZ4LlLh7bwsl7wfG9nOZg3GZ2TTgi8AZ7r7bzNaG5xrD3/Yo+fy3kr5PbrOxYbwieUGVo0h6Hie8Z2Zm84FWd9/by3FPEu8GJTz+l+E50909cPcbgVYO3X4HYB9wVMpz64D/pPeqEeC3wIw04x9HPFnuMbPjiHcXA/wOmGxmZ4RxHhV2J/cWT0Kv73EAJ1Fcu29InlNyFEnP3wCnm9nzxO+lXd7HcV8APhsedynwZ+Hz3zazOjP7DfFE+1zKefcDFyYG5ITP/QdwNPEE2Zv/Bs5JJ3h3f454d+oLwO3Eu3lx9w7gYuB7ZvYc8D/Eq7xHgfcnBuSk+R77swB4IJ1YRXKBduUQyVHhaM8l7n5pP8fcA3zJ3V8eucgGLxyhu8Tdd2c7FpF0KDmK5CAz+x7xrs+Pp2x4m3pcNXCcuz8+YsENkplVAme7+73ZjkUkXUqOIiIiKXTPUUREJIWSo4iISAolRxERkRRKjiIiIimUHEVERFL8fzrVAE5IzjbZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sort the data:\n",
    "por_sort = np.sort(por)\n",
    "\n",
    "# calculate the cumulative probabilities assuming known tails\n",
    "p = np.arange(len(por)) / (len(por) - 1)\n",
    "\n",
    "# plot the cumulative probabilities vs. the sorted porosity values\n",
    "plt.subplot(122)\n",
    "plt.scatter(por_sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7)\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric Porosity 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": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAGWCAYAAAANLkjgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMVklEQVR4nO3deXjU1dn/8fedjSUBAiTsRFkUKyKCAooioIhbESyU+lSt+BNXqqXVp1Zb973UrRYFlL1uaERBqPioBLDKCIZNRfZ9CwlLSAjZ5vz+mAmGGJJJwjDJ5PO6rrky891yz2GYO+d8z2LOOUREROQnEaEOQEREpLpRchQRESlByVFERKQEJUcREZESlBxFRERKUHIUEREpQclRpIYxs3Fm9lCo46goM3vQzN4IdRwigVBylJPOzDab2R4ziy22baSZpYQwrJAzsxFm9mV5xznn7nDOPVGJ6282sxwzy/KX/2Qzi6tctBXnnHvaOTfSH8upZubMLKqy1zOzlmY20cx2mdkhM/vRzB4r+lz5r5/tf78ZZva5mf2mxDVSzOyI/5iixwVVe6cSDpQcJVSigD+EOoiyVOWLO1jMLLKKlxjknIsDugM9gL9V8PdXizIxsybA10A94ALnXAPgMiAe6FDs0K7+99sJmAL8y8weKXG53zvn4oo9vg76G5BqT8lRQmUMcJ+ZxZe208x6m9kSMzvo/9m72L4UM3vCzP7rrzF8amYJ/n1FNZLbzGynv1Zxb7Fze5rZ12Z2wL/vX2YWU2y/M7NRZrYOWOff9rKZbTOzTDP71sz6FDv+UTN7z8z+7Y9llZmdbmYPmFma/7yBxY5vVKy2s8PMnjSzSDP7BTAOuMBfezngP36Kmb1mZnPNLBvo79/2ZLFrDjaz5f74NpjZFeUVvnNuB/Af4Cz/Na4xs+/95ZLij6fo+pvN7H4zWwlkm1lUOcff739vh8xsjZldWqys/u0/bKH/5wH/++1rZvvMrEux6zTz13QTS3kLfwIOATc45zb739M259wfnHMrS3m/6c656cCdwANm1rS8MpLaTclRQmUpkALcV3KHv1YwB/gn0BR4AZhT4gvtt8DNQDMgppTr9AdOAwYCfzGzAf7thcAfgQTgAuBS4K4S5w4BegFn+l8vAc4BmgBvAe+ZWd1ixw8CpgONgWXAPHz/t1oDjwPjix07FSgAOgLd/PGNdM6tBu4AvvbXXuJLvNengAbAMc2uZtYTmAb8L75a08XAZsphZm2Bq4BlZnY68DYwGkgE5gKzi//RAPwPcLX/d7Q/3vFm1gn4PdDDX5u7/DjxXOz/Ge9/vwuAd4AbSvzOz5xze0s5fwDwgXPOW957LeEjfK0WPSt4ntQySo4SSg8Dd5dSM7gaWOecm+6cK3DOvQ38iC8JFZnsnFvrnMsBZuBLXsU95pzLds6tAibj+6LFOfetc26x/7qb8SWuviXOfcY5t89/bZxz/3bOZfjPeR6og6+Zrsgi59w851wB8B6+hPGscy4f3xf+qWYWb2bNgSuB0f7Y0oAXgevKKaePnHP/dc55nXNHSuy7BZjknPs///4dzrkfy7jWh/5a6ZfAAuBp4DfAHP818oF/4Guu7F3svH/6a2Y55Rxf6C+fM80s2jm32Tm3oZz3V2Qq8FszK/peuhHfHx2laQrsCvC6R/njTcf3h06Rf/prwAfMLLWi15TwpOQoIeOc+w74GPhLiV2tgC0ltm3BVxMrsrvY88NAyY4l20qc2wrA3+T5sZntNrNMfMkhoYxzMbN7zWy1v4n3ANCoxDl7ij3PAdKdc4XFXuOP7xQgGthV9GWMLzk3o2zbytjXFgg0+QAMcc7FO+dOcc7d5U92x5S3vza2jWPLu3gMxz3eObceX43yUSDNzN4xs1aBBOac8wDZQF8zOwNf7XrWcQ7PAFoGct3izCwa3x8v+4ptvsdfJvHOue4VvaaEJyVHCbVHgFs59ot4J75EUlwSsKMC121b4tyd/uev4auFnuacawg8CFiJc48uVeO/v3g/MBxo7G/uPFjKOYHYBuQCCcW+jBs65zqX/L3Hi+c41+xQxv5AHFPeZmb4yq94ebtAj3fOveWcu8h/jAOeK+V3Hu89TcXXtHoj8H4pNeUinwHXFqtlBmowvmbtbyp4ntQySo4SUv6axrvAPcU2zwVON7Pf+jt//Abf/b+PK3Dph8ysvpl1xndv8l3/9gZAJpDlr53cWc51GuD7Mt0LRJnZw0DDCsRxlHNuF/Ap8LyZNTSzCDPrYGZFzbp7gDYl7vWVZyJws5ld6r9ea//7qogZwNX+a0QD9+JL4l9V9Hgz62Rml5hZHeAIvppzYSnX2At48d2/LG46cC2+BDmtjJhfwPfvMNXMTgHwv/cXzOzskgebWRMzux4YCzznnMso49oiSo5SLTwOHB3z6P/i+iW+L90M4M/AL51z6RW45gJgPfA58A/n3Kf+7ffh6+ByCHidn5Lm8czD16tzLb6mxCOU3cxZnt/h60D0A7AfeJ+fmge/AL4HdptZQO/VOfcNvuT/Ir4a7QJ+Xusu7xpr8CWjV/DdjxuEb8hHXiWOrwM869++G1+T8YOlXOMwvk5G//U3MZ/v374dSMVXs1xURsz78N3jzAc8ZnYI37/1QXz/7kVWmFmWf9tI4I/OuYcDKBap5UyLHUs4MbNTgU1AtL+DjNQwZjYJ2Omcq9AYTJETqVoM6BURgaN/3PwK3zAXkZAJWrOqmU0y3yDo746z38zsn2a23sxWmpl6iYnUYmb2BPAdMMY5tynU8UjtFrRmVTO7GMgCpjnnzipl/1XA3fgGIvcCXnbO9QpKMCIiIhUQtJqjc24hx44lKmkwvsTpnHOLgXgzq/C4JRERkRMtlPccW3Nsr7/t/m0/m/XCzG4DbgOoW7fuuUlJSSclwHDj9XqJiFAH5YpSuVWeyq5ygl1uXq/36O8o7/cUFBRgBQVEmlHoHC4qiqioqqWOwsLCYx5er/eY56U9nHPHPC9DunOutPl4KySUybG0QdSlvmPn3ARgAkCnTp3cmjVrghlX2EpJSaFfv36hDqPGUblVnsqucoJZbh6Ph/GjRnE1vgmMbx87ll69jn9HqyLHFxYWsn37djZv3syWLVvYsmUL27dvZ+fOnezatYudO3eSlpZGYWFpQ1/BzGjYsCGNGjWiYcOGNGzYkLi4OGJjY48+6tevT7169Y4+6tatS506dY4+rrvuupKza1VKKJPjdo6dxaQNP81iIiIiQbAqNZW+WVn0iI4mKz+fVampZSbHXr16wdixrEpN5fbu3enVqxfZ2dl8//33rFq1itWrV7N27VrWrVvHhg0byM/PP+b8xMREWrVqRatWrTjnnHNo3rw5iYmJRx9NmzalcePGNGnShIYNG1a5xnzddeVNVRyYUCbHWcDvzewdfB1yDvpnEBERkSCxmBgmbdnC4cJC3omM5HcxZU/IlJWVRW5uLgeys3n++edJTU1l48aNR5s269atS8eOHfnFL37B4MGD6dChA6eeeiqnnHIKSUlJ1K1bt8zrV1dBS45m9jbQD0gws+345tCMBnDOjcM3RdhV+GauOIxvlg8REQmijWvW0NU5siMj6eocG0vcptq/fz+LFi1iwYIFpKSksHz5crxe38pgp5xyCueddx6/+93v6NKlC126dKF9+/ZheV85aMnROfc/5ex3wKhg/X4REfk5h2+BzQvwrTjd2TlSU1OZM2cOc+bM4ZtvvsE5R506dTj//PP561//Sq9evejRowfNmpW3gEz40Aw5IiK1SIdOnfjEOaYXFLDEORa+/jrP/uMfmBk9evTg4Ycfpn///vTq1avGNomeCEqOIiK1xHfffceEceNYV1DAMq+XKDPObteOe0aP5sorr6xVNcPyKDmKiISxw4cP89577zF+/Hi+/vprDGgeGck1MTEcPuUU/jJuXJm9VWsrJUcRkRrI4/GwKjWVLv7hFSXt3LmTl19+mfHjx3Pw4EE6derEr4cO5cr16zkzLo6pW7dyypAhSozHEX5djEREwpzH4+H5ESMoeOYZnh8xAo/Hc3Tfjz/+yMiRI2nXrh3/+Mc/uPzyy1mwYAGrV6/m3v/9XxZFRbH9yBGONGvG4KFDQ/guqjfVHEVEaphZycn03LqVwdHRZGVkMCs5mcTERB5++GHeeust6tSpw8iRI7n33ntp37790fNKG9AvpVNyFBGppo7XdOqA/wLtnONz50j77DP+/uKLREdH8+c//5l7772XxMTSpxft1auXkmIA1KwqIlINFc1p2njiRMaPGnVM0+ngoUPJa9uWf5jxf0eOsGLlSkaOHMn69et59tlnj5sYJXCqOYqIVAPFa4ngmwP1amBoy5awa9cxc6Dm5OTwQ24umw8epG/fvrz++uucdtppIYw+/KjmKCISYiVridnZ2XTp3p05QPKuXcwBunTvTnp6OjfeeCP9+/fHzPj4449JSUlRYgwC1RxFREJsVnIyPdPSuCQpCY4cIefw4Z91ntm3bx9DhgwhIyODv/3tbzzwwAPUr18/1KGHLSVHEZEQmjhxIgumTaPJgQNkZWTwTVISN/mTXq9evTjrrLO47777GDduHGeddRaffPIJXbt2DXHU4U/JUUQkRDweD1Mfe4wRhw5RJyqKlNhYOg4aRGxsLADLly9n+PDhrF+/nvvuu48nnniiVs93ejLpnqOISIisSk2lnxnzvV4OFxSQVrfu0YH506dP54ILLuDw4cN88cUXjBkzRonxJFLNUUQkRCwmhkXp6VzgHC97vVx7ww1069aN4cOH89FHH9GvXz/effddTQgeAkqOIiIhsnHNGnrGxdE8Lo4bY2Ko07Ah/fv356uvvuK+++7jmWeeISpKX9OhoFIXEQkBj8fDytmzicrKIjEri4UtWrDqlVfYv38/jzzyCI8++mioQ6zVlBxFREJgVnIyV2dn06FjR17bvZsFu3YRFxfH/PnzOXz4cKjDq/WUHEVEgqi0+VGLao0/ZGTQcM8ePsrPp1Xr1qSkpNChQwdSUlJCG7Sot6qISLAcb37UVampjIiNpUvr1kzLy6NFixakpqbSoUOHEEcsRZQcRUSCZFVqKn2zsuiRm0vfrCxWpaYCvqngXkxP58l162jVoAFvvf22JguvZtSsKiISJBYTw5Rt28gH3gRuiIkBYMWKFfx3yxY6n3kmY199lb59+4Y0Tvk5JUcRkSBxeXn0SUggIz+fPtHRuLw8XnvtNe666y6uuuoqkpOTNbC/mlJyFBEJkk3bt7Nw925uMOOTyEjaLVnC5MmTGTRoEO+99x516tQJdYhyHLrnKCISBB6Ph0XTp9PHjJ1mNI+NZerUqQwYMECJsQZQchQRCYKieVO3mHG4sJD30tPp1KkTH3zwgRJjDaDkKCISBEXzptYvLOSF/HziGzdmwYIFNGjQINShSQCUHEVETjCPx8OXn31G1/h43neOuKgo/vePf9RwjRpEHXJERE6gooH/px46xPO7d5PvHJf/4hf0v+yyUIcmFaCao4jICTQrOZkee/awct8+srxeLrzoIh6YNOno1HFSM6jmKCJyAng8Hj5KTsYzYwa79uxhdX4+XRITeWrMGCXGGkjJUUSkioqaUuumpdFmzx6+yM/n9Pr1ufqmm5QYayglRxGRKvB4PIx76SX6ZWVRJyGBG7Zt49SYGHp16sSQYcNCHZ5UkpKjiEglFdUYT8/OZuLWrWwsLCQyMpJBt9zC9ao11mhKjiIilVC8xnjTaaeRvGsXOw4e5OWXX+aee+4JdXhSRUqOIiIVVLzGOGXbNr48eJBvDx7k5ptvVmIME0qOIiIVNCs5mZ5pafwmKYk9ubmM3bqV7t278/rrr4c6NDlBlBxFRMrh8XhYlZpKl+7dAVg5ezY/ZGSwPz2dyYWFNGjYkI8//pjIyMgQRyonipKjiEgZippQrwbGT5xIy0suYURsLG26dePWFSvIzMlh3uzZtGzZMtShygmk5CgiUoaiJtRLkpLgyBG+BeYASWlprMrKYtiwYQwcODDUYcoJpuQoInIcHo/naBNqVkYG3yQlce/QoRwaOJDfDB9OUlISU6dODXWYEgRKjiIix7EqNZVhzhHfogWzsrLoOGgQvXr14qabbuJgZib/+eQT6tevH+owJQg08biIyHFYTAxTtm1j7+7dbMzKooN/seJp06bx4IMP0rNnz1CHKEGi5CgichwuL48+CQlkNGpEn4QEMvft4/bbb6d79+489NBDoQ5PgkjNqiIix7Fp+3YW7t7NDWZ8EhlJ7syZHDhwgC+++ILo6OhQhydBpJqjiEgpPB4Pi6ZPp48ZO81oHRvL119/zf3330+XLl1CHZ4EmZKjiEgpVqWm0s+MLWYkeL3M3L+ftm3b8re//S3UoclJoOQoIlIKi4lhUXo6pzjHYwUFHC4oYNq0adStWzfUoclJoOQoIlKKjWvW0DMujvyEBPZ7vVzUuzf9+vULdVhykqhDjohICUWD/yMPHeLb3FyiIyN56JFHQh2WnESqOYqIlDArOZmrs7Np16QJ271eBl5xhaaIq2VUcxQRKWbixIksmDaNuP37+TIvj4R69XhQnXBqHdUcRUT8PB4PUx97jBGHDlHo9XIYuGb4cM4///xQhyYnmZKjiIhf0fCNDwoK+LyggKS4OG67885QhyUhoGZVERE/i4lh4d69bMrPJwK4duRIevXqFeqwJARUcxQR8XN5eTStX5+tznF548Z0Pv30UIckIaLkKCLit3bTJuZkZNDWjOwjR7CYmFCHJCGiZlUREXydcd6bMIFcoH90NKckJuLy8kIdloSIao4iIsBXCxeyKzOT1hER9DZjGdCle/dQhyUhouQoIgLMmTePXOe4KjKSl71ezr7hBnXGqcWCmhzN7AozW2Nm683sL6Xsb2Rms81shZl9b2Y3BzMeEZHSrF+/nvnz59O1bl06t23LjR060K5Nm1CHJSEUtORoZpHAWOBK4Ezgf8zszBKHjQJ+cM51BfoBz5uZ7oCLyEl1xx13YM7R0oz83bvVpCpBrTn2BNY75zY65/KAd4DBJY5xQAMzMyAO2AcUBDEmEZFjeDwePv/8c37ZpAl3tGzJ6rg4Og4apCbVWi6YvVVbA9uKvd4OlPy0/QuYBewEGgC/cc55S17IzG4DbgNITEwkJSUlGPGGvaysLJVdJajcKq+6l51zjtGjR9OoUSP6/+lP5NetSy8zmrRtG9K4q3u51QbBTI5WyjZX4vXlwHLgEqAD8H9mtsg5l3nMSc5NACYAdOrUyWlNtcpJSUnRenSVoHKrvOpedvPmzWPlypV0iI9nz+OPs88M16wZEX/5C8OGDQtZXNW93GqDYDarbgfaFnvdBl8NsbibgQ+cz3pgE3BGEGMSEQF8tcaHHnqIpk2aMLxBA7aY0co53W8UILjJcQlwmpm183eyuQ5fE2pxW4FLAcysOdAJ2BjEmEREAJg9ezZLlixh0DXX8HVGBqc4pyEcclTQmlWdcwVm9ntgHhAJTHLOfW9md/j3jwOeAKaY2Sp8zbD3O+fSgxWTiAiA1+vloYceomPHjrRo2pRmcXE0j4vjxpgYEjWEQwjy9HHOubnA3BLbxhV7vhPQ8toiclIlJyezcuVKHn30UZa+8w5RWVkkZmWxLCmJe9WkKmhuVRGpZQoLC3nkkUc488wzaZmQwDDniG/RgllZWRrCIUcpOYpIrfL222+zevVq3n//fQ4cOMCkLVu4rrCQ9ZGR9O7UKdThSTWhuVVFpNYoKCjgscce45xzzuHaa69l45o1dHWO7MhIujrHxjVrQh2iVBOqOYpIrfHuu++yfv16Zs6cSUREBA7YDFwALATOCml0Up2o5igitYLX6+Wpp57irLPO4pprrgGgQ6dO7AA+LSxkh/+1CKjmKCK1xAcffMDq1at5++23iYjw1Qs2rlnDgEaNaBUXx+kxMVrcWI5SzVFEwp5zjieffJLTTz+dX//614BvwvGVs2ezNitLK3HIz6jmKCJhb86cOaxYsYIpU6YQGRkJwKrUVA3jkONSchSRsOac44knnuDUU0/lt7/9LR6Ph1WpqWzavp2vtm3jenxzVmoYhxSn5CgiYe2zzz7jm2++Yfz48aSmpjJ+1CiuBqalpdEvIQEaNmRgQYHuN8oxdM9RRMLak08+SevWrbnppptYlZpK36wseuTm0s+MZUDjOnVYFxur+41yDNUcRSRsLV68mIULF/LCCy9Qp04dLCaGKdu2kQ8sAnqPHs3+Nm24vXt33W+UYyg5ikjYGjNmDI0bN+bWW2/F4/Hw5Wef0adYU2pimzaMvPPOUIcp1ZCaVUUkLK1du5aZM2dy11138f333zN+1ChOX76cRenpROfnqylVyqSao4iEpeeff56YmBh69+7NuJdeol9WFjeddhoAKeecwx2jR6spVY5LyVFEws6ePXuYOnUq559/Pq/deSfdgCnpvnXU18XFKTFKuZQcRSTsvPLKK+Tl5ZG9Zg03ZGZyXnQ0JCSQ0q2bEqMERPccRSSsZGVlMXbsWLp17cqVMTHM93rxHDnCMlBilIApOYpIWHnjjTc4cOAAV1x1FYvS0znFOV72ejn7hhuUGCVgalYVkbDg8XhYvmQJzz33HH369KFdmzbEJSWRGBXFjf5hGyKBUnIUkRrP4/Hw/IgRNNq7l90ZGfzpT3+iS/fujI+N5XTgqzp1uETDNqQClBxFpMablZxMz61beTc3lyZm7Nuzx9eEOnYsq1JTNQOOVJiSo4jUeA74KD+fpYWFdI6MxPyLGffq1UtJUSpFHXJEpMbr0KkT3xcWEg1ERUbSQctPSRWp5igiNV7Gnj0c8Hq5MDaW/k2bavkpqTLVHEWkxpu/aBEOuLqwkEXp6VhMTKhDkhpONUcRqdGys7NJSUnh9Dp1qNeqFQNjYlRzlCpTzVFEaiyPx8PQIUM4cuQIzYH83btZBlptQ6pMyVFEaiSPx8O4u+7Ck5JCm4gI/nb66Wxs2pSOgwaph6pUmZKjiNRIq1JTSUxL40BBAQkREWw8cIAjzZoxeOjQUIcmYUDJUURqJIuJYdKuXTQAiIjAc/HF3D52rGqNckKoQ46I1Eip33zDvsJCBsTH079FCxIvvFCJUU4Y1RxFpMbxeDzMe/99HNDTvxyVOuHIiaSao4jUOMnvvMPuAwe4qGFD9jVqpE44csKp5igiNYrH4+Hj6dPJ9no5MyeH9NhYdcKRE07JUURqlI/ef5/MfftIMmO/c9Q/91zVGuWEU3IUkRplx65d7HCOAZGRHImMpFWrVqEOScKQkqOI1Chr1q8nAsgFdoBW4JCgUIccEakx9u7dy9KlSzmnXj16tGxJZ82jKkGimqOI1BiPP/44hYWFxDuneVQlqJQcRaRG8Hq9vPvuu3SqV497WrVidVychnBI0Cg5ikiN8Nlnn7F3716soICd27axPjNT9xslaJQcRaRGeO2116hfvz6XmJEdGUlX59i4Zk2ow5IwVW5yNLPfm1njkxGMiEhptm/fzuzZs2l/6qlscY62kZFsi4zEhTowCVuB1BxbAEvMbIaZXWFmFuygRESKe+ONN/B6vSQeOcKhiAg+yMsjMzFRM+NI0JSbHJ1zfwNOAyYCI4B1Zva0mXUIcmwiIhQUFPD666/T7tRTGZ6fz986daJpixb0Gj5cnXEkaAK65+icc8Bu/6MAaAy8b2Z/D2JsIiLMnj2bnTt30iQvj3kZGaxYv17zqUrQlTsJgJndA9wEpANvAP/rnMs3swhgHfDn4IYoIrXZ+PHjaRAXxwgzTu/YkZn792sIhwRdIDPkJAC/cs5tKb7ROec1s18GJywREZg5cybz5s3j1EaN+GzfPnL27SM9KYmbVGuUIAukWbVdycRoZtMBnHOrgxKViNR6Ho+Hx+66C4Df5uZyV8eObGzaVLVGOSkCqTl2Lv7CzCKBc4MTjoiIz/IlS9iSnk73evVY5fXSdv9+jjRrpnuNclIct+ZoZg+Y2SHgbDPL9D8OAWnARyctQhGplb5bs4YDBQWcX1jIDsBz8cXcPnasao1yUhy35uicewZ4xsyecc49cBJjEhHhk7lziY2IoEPr1rSpU4fECy9UYpSTpqya4xn+p++ZWfeSj5MUn4jUQmPGjGH9xo20NMO7Z49W35CTrqx7jvcCtwLPl7LPAZcEJSIRqdU8Hg//evJJAO6IieEHrb4hIVBWs+qt/p/9T144IlLbrVi6lAPZ2bSIiCDO6yWtbl1uUyccOcmOmxzN7Fdlneic++DEhyMitd0P69aRWVjIJVFRvOz1cu0NN6jWKCddWc2qg8rY5wAlRxE54RampFA/IoLzEhLoEhNDUps2oQ5JaqGymlVvPpmBiIjs2bOHld99R8uICJpnZvImcENMTKjDklqorGbVG5xz/zazP5W23zn3QvDCEpHaaNq0aRQWFnJ5kyZkx8czMCYGl5cX6rCkFipr+rhY/88Gx3mUy7/+4xozW29mfznOMf3MbLmZfW9mCyoQu4iEEeccY8eOpVGdOuw+dIjcHTs0hENCpqxm1fH+n49V5sL+aebGApcB2/EtmDzLOfdDsWPigVeBK5xzW82sWWV+l4jUfF999RVbtmzh1zEx3BQTw7S8POqfe64640hIlDvxuJm1N7PZZrbXzNLM7CMzax/AtXsC651zG51zecA7wOASx/wW+MA5txXAOZdW0TcgIuFh4sSJREVFken1kgXkRUfTqlWrUIcltZT51jEu4wCzxfhqgG/7N10H3O2cK/PPOTMbhq9GONL/+kagl3Pu98WOeQmIxje5eQPgZefctFKudRtwG0BiYuK5M2bMCOjNybGysrKIi4sLdRg1jsqt8gItu+zsbIYOHUqP7t25YdgwooGC6Ghat29PbGxsueeHG33mKq9///7fOufOq+p1AlmVw5xz04u9/reZ/f64Rxc7r5RtJTNxFL4VPi4F6gFfm9li59zaY05ybgIwAaBTp06uX79+Afx6KSklJQWVXcWp3Cov0LJ7/fXXyc3Npd+OHbQYO5ZZWVk0v+kmbrv77uAHWQ3pMxd6ZfVWbeJ/Ot/fmeYdfMntN8CcAK69HWhb7HUbYGcpx6Q757KBbDNbCHQF1iIitcbEiRNp1aoVi/bu5RQzNgK9O3UKdVhSi5VVc/wWXzIsqgHeXmyfA54o59pLgNPMrB2wA19z7G9LHPMR8C8ziwJigF7Ai4GFLiLh4Pvvv8fj8XBpv36c+8MPZMfFaQiHhFxZvVXbVeXCzrkCf/PrPCASmOSc+97M7vDvH+ecW21mnwArAS/whnPuu6r8XhGpWYo64uSuX8/3+/fT9OBBlrVrx70awiEhFMg9R8zsLOBMoG7RttI6zpTknJsLzC2xbVyJ12OAMYHEISLhJS8vj+nTp9OhXTsG79jBLzSEQ6qJQIZyPAK84n/0B/4OXBPkuESkFpgzZw7p6emc3aULCwsLySwsJCciQkM4JOTKTY7AMHy9SXf751vtCtQJalQiUitMmjSJli1bMuDyy9kBfFpYyA6ggzrjSIgF0qya45zzmlmBmTUE0oBAJgEQETmuXbt28Z///If77ruPLevXM6BRI1rFxXG6OuNINRBIzXGpf5q31/H1YE0FvglmUCIS/qZPn05hYSHdu3dn5ezZrM3KIn/3bs2nKtVCuTVH59xd/qfj/D1LGzrnVgY3LBEJZ845Jk+eTO/evVmxdClXZ2fToWNHZu7fT8dBg9QZR0IukJojZvYrM3sBuBvoENyQRCTcLV68mB9//JG+ffuyeMYMPt69m2Vr1pAeG8vgoUNDHZ5I+TVHM3sV6MhPc6vebmYDnHOjghqZiIStyZMnU79+ffKys7ly714N4ZBqJ5AOOX2Bs5x/hnIzmwqsCmpUIhK2srOzefPNN+napQvpGRlsANpFRpIXHU0HDeGQaiKQ5LgGSAK2+F+3xTejjYhIhY0ZM4bDhw9zxf79LPn2Ww4nJvK510uBmlSlGilr4vHZ+OZQbQSsNrOiHqo9ga9OQmwiEoamTJ5MQmQkoxISSMnN5dtBg2jfrh03de+uJlWpNsqqOf7jpEUhIrXCzJkz2bJ1K2dERzN5+XK+SUri3qFDlRSl2ilr4vEFRc/NrDnQw//yG+dcWrADE5Hw4vF4ePzhhzHgmbPP5tO0NA3bkGorkLlVh+Mb9P9rYDjgMbNhwQ5MRMLHxIkTeeLXv2bz6tU0iojgYFYWR5o10z1GqbYC6ZDzV6BHUW3RzBKBz4D3gxmYiIQHj8fD1Mceo1daGnMKCxnctCkp3bpxx+jRqjVKtRXIJAARJZpRMwI8T0SEVamp9DMjOT+fukBevXpKjFLtBVJz/MTM5vHTJAC/ocQajSIix2MxMczfu5ftXi+xZnT73e+UGKXaKzM5mpkB/8TXGeciwIAJzrmZJyE2EQkDLi+PJvXqkZ+Tw5BmzWjXpk2oQxIpV5nJ0TnnzOxD59y5wAcnKSYRCSMWE8OC/ftpBmw4cICLYmJCHZJIuQK5d7jYzHqUf5iIyM99u3gxB53j3KgozgE2rlkT6pBEyhVIcuyPL0FuMLOVZrbKzDR9nIgE5CuPB4BfRkezLTISF+J4RAIRSIecK4MehYiEpUOHDrF29Woam7GgoIDMVq00tlFqhLLmVm0GPIhvuapVwDPOucyTFZiI1HyLFy8mp6CA0UlJHHCODsOHq6eq1AhlNatOA7KBV4A4fL1WRUQC4vF4mP/ZZ0QBDfbu1ULGUqOU1azawjn3V//zeWaWejICEpHwsOSrr0hdvpxLGjRgfb16mkdVapSykqOZWWN8YxsBIou/ds7tC3ZwIlJzfbtyJbl5eZxdWMg3+fn07tQp1CGJBKys5NgI+JafkiNAUe3RAe2DFZSI1Hyff/YZLRITaZaZSVfnNIRDapSylqw69STGISJhZMOGDWzbvp1rr7qKpAUL+Mo5fhHqoEQqQBOIi8gJ98wzzwDQ89xz+SAvj8zERHXGkRpFyVFETiiv18uHM2dyVv36/KJZM+IaN6aXhnBIDaPkKCIn1IIFC8jYt4/CggLy8/PZmJVFB3XGkRomoORoZheZ2c3+54lm1i64YYlITTVlyhTq1a3L4GbNKIiMpE9CAi4vL9RhiVRIucnRzB4B7gce8G+KBv4dzKBEpGY6dOgQ77//Pp3OOIMv09KIKihgUXo6ppU4pIYJpOZ4LXANvtlycM7tBBoEMygRqZnef/99Dh8+jNu5kz5m5IFqjlIjBZIc85xzDt/YRswsNrghiUhNNWXKFJo0bsxA59geFUWMGcuALt27hzo0kQoJJDnOMLPxQLyZ3Qp8Brwe3LBEpKbZsGEDCxcupFlUFOuys4kvKCA9MpIhDz+snqpS45S7ZJVz7h9mdhmQCXQCHnbO/V/QIxORGmXatGmYGbc1aED7Bg2YlZVFg6ZNuf7660MdmkiFlZsczeyPwHtKiCJyPF6vl2nTpnHmmWcya9Mmrgc2AlfWqRPq0EQqJZBm1Yb4VuVYZGajzKx5sIMSkZpl4cKFbN68mY7t2tEzLo7sFi0YmJQEzoU6NJFKKTc5Oucec851BkYBrYAFZvZZ0CMTkRpjypQp1KtXjz3ffMPqgwfJ372bZUC9+vVDHZpIpZTbrFpMGrAbyACaBSccEalpsrKymDFjBo0iIrglK4s6UVGkxMbScdAgYmPVuV1qpkAmAbjTzFKAz4EE4Fbn3NnBDkxEaob33nuPnJwcroyLY77Xy+GCAtLq1tVE41KjBXLP8RRgtHOus3PuEefcD8EOSkRqjilTptC8eXM2ZmZyinO87PVy9g03aPiG1GjHbVY1s4bOuUzg7/7XTYrvd87tC3JsIlLNFY1tvPaaazhr+XJi8/MZFh1NUps2oQ5NpErKqjm+5f/5LbDU//PbYq9FpJabMmUKERERtDnlFObv3k38vn2aS1XCwnFrjs65X/p/agUOEfkZr9fL1KlT6dmzJ8s//JA+Zuw001yqEhYC6ZDzeSDbRKR2+eKLL9i2bRtnn3km/czYYkYr5zSXqoSF4yZHM6vrv8+YYGaNzayJ/3EqvvGOIlKLTZ48mfj4eLr16MGi9HR1xpGwUtY4x9uB0fgS4beA+bdnAmODG5aIVGcHDhzggw8+4KqrrsKzYAF9EhJIatiQGwsKSFRnHAkDZd1zfBl42czuds69chJjEpFqbsaMGRw5coS8H37gdODT9HRG1KvHurg4LlGTqoSBQFbleMXMzgLOBOoW2z4tmIGJSPU1efJkEhISuCori+tOOQWAlHPO4Y7Ro9WkKmEhkFU5HgH64UuOc4ErgS8BJUeRWuidd95h8eLFdIyP59N9+8jet49lSUncq8QoYSSQGXKGAZcCu51zNwNdAa1DI1ILeTwenho9GgP+JzeXuzp2ZGPTpnQcNEiJUcJKIBOP5zjnvGZWYGYN8U1A3j7IcYlINbR8yRK2pafTo359lhcW0mr/fo40a6Z5VCXsBJIcl5pZPPA6vl6rWcA3wQxKRKqn79eu5WBhIT0KCvjaDM/FF3PH3Xer1ihhJ5AOOXf5n44zs0+Ahs65lcENS0Sqoy8XLiQ2IoJWTZpwdUwMSRdeqMQoYamsiceP2x/bzLo751KDE5KIVEdpaWms/O47WkVE0CwzkzeBGzSHqoSpsmqOz5exzwGXnOBYRKQa+/e//01hYSGXNWlCdnw8A2NiNIeqhK2yJgHofzIDEZHqyznH2LFjaRgTw55Dh8jNzmZZu3bcqwH/EqYCGef4u9K2axIAkdpj6dKlbNy4kaExMdwcE8O0vDzqn3uu7jdK2Aqkt2qPYs/r4hvzmIomARCpNSZNmkRUVBTZQGZhITkREXRopfUHJHwF0lv17uKvzawRMD1oEYlItXL48GHeeustevToQVpqKp8WFrIjIoLBnTqFOjSRoAlkhpySDgOnBXKgmV1hZmvMbL2Z/aWM43qYWaGZDatEPCISRDNnziQzM5OGdetydfPmXNipE8PatVNnHAlrgdxznI2vdyr4kumZwIwAzovEt7TVZcB2YImZzXLO/VDKcc8B8yoWuoicDC+88AINYmLos3Mnn2n1DaklArnn+I9izwuALc657QGc1xNY75zbCGBm7wCDgR9KHHc3kMyx9zZFpBpYv349qampDE1I4MboaHITEkjp1k2rb0jYC+Se4wIA/7yqUf7nTZxz+8o5tTWwrdjr7cAx/5vMrDVwLb4xk8dNjmZ2G3AbQGJiIikpKeWFLaXIyspS2VVCbS63119/HTOjx1134WnUiNbA2UlJ5OTkBFQmtbnsqkLlFnqBNKveBjwB5ABewPA1s5Y3+biVss2VeP0ScL9zrtCstMP9Jzk3AZgA0KlTJ9evX7/ywpZSpKSkoLKruNpabgUFBfz2t7+lQ/v27HjuOSLN2OocDe65h6f+/veArlFby66qVG6hF0iHnP8FOjvnTnXOtXfOtXPOBbIqx3agbbHXbYCdJY45D3jHzDbjWxrrVTMbEsC1RSTI/vOf/7Br1y7atm7NRq+XtpGRbIuM/NlfuCLhKJB7jhvw9VCtqCXAaWbWDtgBXAf8tvgBzrl2Rc/NbArwsXPuw0r8LhE5wd544w2aNm1KvT17OBQRwQd5eWS2aqXlqaRWCCQ5PgB8ZWYeILdoo3PunrJOcs4VmNnv8fVCjQQmOee+N7M7/PvHVT5sEQmmXbt2MWfOHC679FKGb9lCfOvWzMrKosPw4eqII7VCIMlxPPAFsArfPceAOefmAnNLbCs1KTrnRlTk2iISPFOnTqWwsJA+ffsy5amnuB7YCPTWwH+pJQJJjgXOuT8FPRIRqRacc7zxxhv07duXZo0b0ychgYz8fPpER2vgv9QagXTImW9mt5lZSzNrUvQIemQiEhILFixgw4YNdDz1VDZt386i9HSaHjzIovR0TOs3Si0RSM2xqBPNA8W2BTKUQ0RqoGeeeYaYyEguWbmScV98Qb+EBGjYkIEFBao5Sq0RyCQA7co7RkTCQ0ZGBl988QX9GjTgovx8fjRjGdC1Th2+qlNHU8ZJraH1HEXkqKlTp1JQUMCBI0f4dONGFgG9R49mf5s23N69u3qqSq2h9RxFBPB1xJkwYQId2rXjysLCo51wktq0YeSdd4Y6PJGTqtwOOc65u4s9bgW6AborLxJmFi1axJo1a/hF587M372b+H371AlHaq2grucoIjXHhAkTiIuLY9+yZfQxY6cZfRIS1AlHaqWgrecoIjVHRkYG77//Pr3PP5+LNm1igxkXO8fHwNXqhCO1UDDXcxSRGmL69Onk5ubSt39/UpYs4QLneNnr5dobblAnHKmVjpsczawj0LxoPcdi2/uYWR3n3IagRyciQVfUEef8888nLzubnnFxNI+L48aYGBLbtAl1eCIhUdY9x5eAQ6Vsz/HvE5Ew8OWXX7J69WoGDBjAytmzWZuVRf7u3SwDuqhJVWqpsppVT3XOrSy50Tm31MxODV5IInIyvfbaazRs2JAWTZsyzDniW7RgVlYWHQcNUpOq1Fpl1RzrlrGv3okOREROvrS0NN5//31GjBjBzr17eWPTJnbv2MHGrCw6aAUOqcXKSo5LzOzWkhvN7Bbg2+CFJCIny6RJk8jPz+eCCy5g0fTpGsIh4ldWs+poYKaZXc9PyfA8fBMAXBvkuEQkyAoLCxk/fjz9+vUja/9++plpCIeI33GTo3NuD9DbzPoDZ/k3z3HOfXFSIhORoJo3bx6bN2/mueee49ChQyxKT9cQDhG/QFblmA/MPwmxiMhJ9Oqrr9K8eXOGDBnCY3/7m4ZwiBRTmenjRKSG27x5M3PnzuXWW29l2bJlGsIhUoKSo0gtNGHCBMyM2267jVWpqQxzjhEtWrA6Lk5DOERQchSpdfLy8pg4cSK//OUvadu2LRYTw6QtW9i5bRvrMzM1hEMEJUeRWic5OZm0tDTu9K/RuHHNGro6R3ZkJF2dY+OaNSGOUCT0lBxFapl//vOfnHbaaQwcOBDwLbmzGWjt/+mOe6ZI7aHkKFKLLF26lMWLFzNq1CgiInz//Tt06sQO4NPCQnb4X4vUdoEsWSUiYeKVV14hLi6OESNGHN3m8vK4unlzYvPzaRsdrZlxRFDNUaTWSEtL45133uGmm26iUaNGR7dbTAyL0tNpevAgi9LTsZiYEEYpUj2o5ihSS7z++uvk5eXx+9///pjtLi+PPgkJZOTn00c1RxFANUeRWiE/P5/XXnuNyy67jDPOOOOYfZu2b2f+7t3E79unmqOIn2qOIrXAzJkz2bFjB6+99tox2z0ej1bjECmFao4itcArr7xC+/btueqqq47ZPis5mQtyc9keFUXbiAhNHSfip+QoEuaWL1/Ol19+yahRo4iMjDy63ePxHJ1TNb6ggNcbNGDIww9r6jgR1KwqEvZefPFFYmNj+X//7/8ds71oTtX4Fi2YlZXFpb/7HbfcckuIohSpXlRzFAljO3fu5O233+aWW24hPj7+mH0WE8OUbdvYu3s3G7OyNPhfpBjVHEXC2NixYykoKOAPf/jDz/ZtXLOGnnFxZMfFMTAmRh1xRIpRzVEkTB0+fJhx48Zx7bXX0r59+2P2eTweFs+Ywff795O7Y4c64oiUoJqjSJiaNm0a+/bt449//OPP9s1KTubKvXv5RUwM0/LyqH/uueqII1KMkqNIGPJ6vbz44ov06NGDCy+88Jh9Ho+HFcuXE1FYSLs6dciLjqZDq1YhilSkelJyFAlDc+fOZe3atbz99tuY2dHtHo+H8aNGcWF2Nu8Dc+vWpaBJEwYPHRq6YEWqISVHkTD0wgsv0LZtW4aWSHqrUlPpm5VF/+hocpo3Z1vv3vxt9Gg1qYqUoOQoEmaWLVvG/PnzGTNmDNHR0cfs27R9Ows3beKIGYsiI7lhwAAlRpFSqLeqSJj5+9//TsOGDRk5cuQx2zWPqkjglBxFwsiGDRuYMWMGd9xxx88G/WseVZHAKTmKhJHnn3+eqKgoRo8efcx2zaMqUjG65ygSJvbs2cOkSZO46aabaNmy5TH7NI+qSMWo5igSJv75z3+Sl5fHfffd97N9mkdVpGJUcxQJA5mZmbz66qv86le/4vTTT//ZfpeXR5+EBDLy8+kTHa2OOCLlUM1RJAxMmDCBAwcOcP/995e632JiWJCWRsP0dBakpWExMSc5QpGaRclRpIbLzc3lxRdf5JJLLqFHjx6lHrNxzRq6Okd2ZCRdnWPjmjUnOUqRmkXJUaSGmzp1Kjt37uQvf/nLcY9xwGagtf+nOymRidRcSo4iNVh+fj5PP/00vXr1YsCAAcc9rkOnTuwAPi0sZIf/tYgcnzrkiNRg06dPZ8uWLbz66qvHTDBe0sY1axjQqBGt4uI4XQsbi5RLNUeRGqqgoICnnnqKc889lyuvvPK4xxWfACB/927NjCMSANUcRWqot956i40bN/LRRx8dt9bo8XgY99JLDHOOM7p1Y+rWrXQcNEgz44iUQ8lRpAYqLCzkySefpGvXrgwaNKjUYyZOnMiHjz9ON2BKejojgCPNmmntRpEAKDmK1EDvvvsu69atIzk5udRao8fjYepjj/G7jAzOi46GhARSunXjDq3dKBIQJUeRGqao1njWWWcxZMiQUo8pWoFjoRnk57MMtKixSAUoOYrUMO+++y6rV6/m3XffJSLi533qijrgRGVl0bawkNfj47lNK3CIVIh6q4rUIPn5+Tz88MOcffbZDBs27Gf7i3fA+Uu3bhQ0a6YVOEQqQTVHkRpk8uTJbNiwgdmzZ/+s1ujxeBg/ahSnZ2czZds2dcARqQIlR5Ea4siRIzz++ONccMEFXH311Ue3ezweVqWmsmnTJq4GhnbsCEDKOeeoA45IJQU1OZrZFcDLQCTwhnPu2RL7rweKlhHIAu50zq0IZkwiNdVrr73Gjh07mD59+tEeqkXDNYbVq8dKM9YB7NrFuthYJUaRKghacjSzSGAscBmwHVhiZrOccz8UO2wT0Nc5t9/MrgQmAPrfLFLCoUOHePrppxkwYAD9+/cHjh2u0TU6Glq0YM0117C/XTtu795diVGkCoJZc+wJrHfObQQws3eAwcDR5Oic+6rY8YuBNkGMR6TGeumll0hPT+epp54Cfup408+Mr6KjIT+fj3Jy+NvQoUqKIieAORecxWvMbBhwhXNupP/1jUAv59zvj3P8fcAZRceX2HcbcBtAYmLiuTNmzAhKzOEuKyuLuLi4UIdR44S63A4ePMj1119Pt27deOKJJ8jOziZ961bqeL1k5uURFxXFIaBpy5YkJCSELM7ShLrsaiqVW+X179//W+fceVW9TjBrjqVN9lhqJjaz/sAtwEWl7XfOTcDX5EqnTp1cv379TlCItUtKSgoqu4oLdbn94Q9/ICcnh3HjxtG5c2f++uc/0/add7gmKYlxGRl8XY073oS67GoqlVvoBTM5bgfaFnvdBthZ8iAzOxt4A7jSOZcRxHhEapy1a9fy6quvMnLkSDp37nx0gP8PGRlkZWSwLCmJe6tpYhSpyYI5CcAS4DQza2dmMcB1wKziB5hZEvABcKNzbm0QYxGpke6//37q1q3L448/Dvimhbs6O5u7OnZkY9OmWmFDJEiCVnN0zhWY2e+BefiGckxyzn1vZnf4948DHgaaAq/6u6YXnIi2YpFwsGDBAj788EOefPJJmjdvzsSJE1kwbRpNDhzg4owM0pOSuEkD/EWCIqjjHJ1zc4G5JbaNK/Z8JPCzDjgitZ3X6+W+++6jdevW/PGPfzw6bGPEoUPUiYoiJTZWtUaRINLcqiLV0Ntvv83SpUt5+umnWbVq1THDNnKdI61uXU0LJxJEmj5OpJrJycnhwQcfpHv37px22mlH50v9ND2dPgkJTHOOm7TKhkhQKTmKVDPPPvssW7duZfivfsXHM2ceM1/q2nPOYYx6p4oEnZKjSDWyYcMGnn32Wdo2bEiH5GTei4hgXb16mi9V5CRTchSpJpxz3HPPPeAcN+XlMTgzk6z8fL679lr2X3ih5ksVOYmUHEWqidmzZzN37lzO696dZatWsbCggP8Cv2jVipF33hnq8ERqFfVWFakGcnJy+MMf/kC7du1olp3NoYgIPsjLIzMxUb1SRUJAyVGkGrj77rvZvHkzvXv25P/FxfH3886jaYsW9Bo+XE2pIiGg5CgSYu+99x5TJk3iosaNObhiBe9mZ7P9yBGONGumWqNIiOieo0iIeDweVixdypPPPkvdiAhmnHsuX+3fz7eXXKIFi0VCTMlRJAQ8Hg/PjxhB1J49bNu/n/NatOCr/fuZA9yuBYtFQk7JUSQEZiUnc/rmzbx85AgdIyI485JL2H/RRaotilQTSo4iIbB9xw4+yc0lH2gXE0Pr1q01XEOkGlGHHJGTzOPxsGr+fNKc44yICAqbN1fHG5FqRjVHkSDzeDysSk2li7/J9Mv581mflsaZcXGcHx9PEw3XEKl2lBxFgsjj8TB+1CiuBsZPnIj7179InjWLw14vtyYlsbJePdUaRaohJUeRIFqVmkrfrCx6REeTlZ/PKy+9xNdff829995LXIcO6oAjUk0pOYoE0abt21m4aRNHzJhoxvKNG7niiisYM2YMZhbq8ETkONQhRyRIPB4Pi6ZPp48Z24BthYXUqVOHSZMmKTGKVHOqOYoEyazkZC7IzWV7VBTpubnsLijguSefpGXLlqEOTUTKoeQoEgQej4eVs2cTlZWFNz+f/xQU0KdPH/785z+HOjQRCYCaVUVOMI/Hw7iXXmKYc4w480w+LSwkMSGBTz75JNShiUiAVHMUOYGKhm6cnp3NpK1b2e4chcDYV1+lfv36oQ5PRAKkmqPICVJUY+yblcUDHTuSV6cOG48c4bEnn+TXv/51qMMTkQpQzVHkBJg4cSIfPv443YAp6ekszsxk8YEDXHfddTzwwAOhDk9EKkjJUaSKPB4PUx97jN9lZHBedDS7Gzbkjb17Ofvss5k2bVqowxORSlByFKmiVamp9DNjvtdLRk4O07OyaOrvgBMdHR3q8ESkEnTPUaSKLCaGRenptPR6eTQ/H29kJPPmzdN4RpEaTMlRpIpcXh4XNW3KJxER5DrHbbfcQteuXUMdlohUgZpVRarIYmKYumsX2woK6BAdzTnnnhvqkESkipQcRSqo5PqMn/7nP2wrKKBvkyYMTEzE5eWFOkQRqSIlR5EKKBqyMaxePcbHxTH7iiuYkZzMqfHx3HXqqXxixqXdu4c6TBGpIiVHkQAVH7LRNTqaJbGxPP3001x++eX89a9/Zc1332l9RpEwoeQoEqDiQzZWHT7MqwcP0uXss0lOTiY2NpY+ffqEOkQROUHUW1UkQEVDNgoKC3klP59mzZuzYMECYmNjQx2aiJxgqjmKHEdRx5uWSUmAb8hGm8aNeXvnTlrUqcP//ulPxMfHhzZIEQkKJUcRv+K9UAHGjxrF1UD6Lbfg8XjYnZHBmzt20DE2lh4dO9K7b9/QBiwiQaPkKMJPS01dDYyfOJGWl1zC1cDQli35EHjtX/9i2ptvck63bvy/G2+kR+/e6ngjEsaUHEXwdbbpm5VFj+hosvLzWQPMAdzOnfz744/5aO5cBg4cyMyZM7Uuo0gtoA45Ivg620zZto1PN25kyrZtdOjUiZtffJExkZF8NHcuN998M7Nnz1ZiFKklVHOUWqv4PUaXl8fApCSIimJgQQGZ+/bx4IMP8s2SJYwcOZIJEyZgZqEOWUROEiVHqZVK3mO88M47WRcby+nArMJCvn3pJTIyMnjnnXdo3ry5EqNILaNmVak1PB4Pb7z22tEaY1GHm6vxDdO47V//YuZpp/HJ+vVERESwYMECfvOb34Q6bBEJASVHCSvFE2DJ7eNHjaLxxImMHzUKi4lhDpC8axdzgA5nnMErr7zCm++8w6WXXsqyZcvUG1WkFlOzqoSNkk2ljB17NMEVrymyaxf78/K4fexYVqWm0js6mjvuuIP169fz5JNP8sADDxARob8bRWozfQNIjVdUW5yVnHxMU+mq1NSjx3Tp3v2YmmKX7t3p2rUrazdt4vbbb+fIkSN8/vnn/PWvf1ViFBHVHKVmK15b/Cg7m3UA/gR4e7Glo3r16gX+mmLR9u7du7N69Wpuu+02xowZQ8OGDUPxFkSkGlJylBqn+BCMks2l315yCfvbtSt16ahevXrRrl07HnjgASZNmkSbNm2YN28eAwcODM0bEZFqS8lRapTShmDMgZ9qi0OHltqRpqCggLFjx/LII4+QnZ3Nfffdx0MPPaTaooiUSslRapRZycn0TEvjkqQkOHLkmI41pdUWvV4v7733Ho8++ig//vgjAwcO5OWXX+aMM84I0TsQkZpAyVFqDI/Hw+IZM1ixezeZaWksbdeOe/0JsWRSdM4xa9YsHn74YVauXEnnzp358MMPueaaazSgX0TKpW55UmPMSk7myr17uTMmhm+9Xuqfe+7PkmJeXh7//ve/6d69O0OGDCEnJ4c333yTFStWMHjwYCVGEQmIao5S7RV1wNmxcyc/AO0iI8mLjqZDq1ZHj9m3bx8TJkzglVdeYefOnZx55plMnjyZG264gagofcxFpGL0rSHVWvEOOBnZ2RxOTORzr5eC2Fh+OWQIn376KZMmTWLmzJnk5eVx2WWXMXHiRC6//HLVEkWk0pQcpVoqqi1u2rTpmKEaS3/5S6Kjozm8cyf/8z//w9atW2ncuDG33347t956K126dAl16CISBpQcpdopObB/rXOsWbeOtw8cYP/bb7Njxw4iIyO59NJL+fvf/87gwYOpW7duqMMWkTCi5CghUXwgf8lONatSU+mZl0dmRAT79uwhNSeH948cITIykgEDBvDoo48yZMgQEhISQhS9iIQ7JUc54cpKfEX7iw/kL3j5ZWJjY1m8eDGLFi3i888/Z8+ePQDUi45mwOWXc+ONNzJgwACaNGlykt+NiNRGSo5yQpW1Mgb4epUmz5hB1N69fG7Gor17md6vHwUFBQC0aNGCiy++mLZt29Kgbl2uGDSI888/P0TvRkRqKyVHCVh5NUKAFUuX0js/n7YNGtB41y6ee/ppGicksH79etauXcvu3buPHlsvIoJG9eszfPhwBg0aRI8ePWjfvr16mYpIyAU1OZrZFcDLQCTwhnPu2RL7zb//KuAwMMI5l/qzC0nIOOfIycnhiy++YOL999O9oIC3CgroNXw4DRo0YNeuXezatYvdu3ezc+dOtm3bdrQWCMDGjbRo0YKOHTtyxRVX0LlzZzp37kxubi57d+7k7FIG8ouIhFrQkqOZRQJjgcuA7cASM5vlnPuh2GFXAqf5H72A1/w/w5Jz7piH1+st9bXX6z36KCwsPPqztEdBQcHRR35+/jGPvLw88vLyyM3N5ccff2TPnj18+OGHNG/enJycHI4cOUJ2djaHDx8mOzub7OxsDh06dPSRmZnJwYMHycvLO/oePvT/nP+s7++c+Ph4WrZsSYsWLbjgggsYPnw4hYWFHMnKoteFFzJkyBDi4uJOfmGLiFRBMGuOPYH1zrmNAGb2DjAYKJ4cBwPTnHMOWGxm8WbW0jm363gXXbduXZW+bH2/KrD9JY893uvSfhY9AvmdoVKvXj3q169PbGws9evXp379+jRo0IA2bdrQoEEDGjRoQOPGjYmPj+fAgQMsnjGDvlFRrIiK4pbnnuOyyy7TEAoRCUvBTI6tgW3FXm/n57XC0o5pDRyTHM3sNuA2/8vc7Ozs705sqOHNILEFJBSCRYLbDekO9ubk5JCTk0NGRkagl4pdCPUdHP7ommuygxlzNZMApIc6iBpKZVc5KrfK63QiLhLM5Fhar4qSVahAjsE5NwGYAGBmS51z51U9vNpHZVc5KrfKU9lVjsqt8sxs6Ym4TjBX5dgOtC32ug2wsxLHiIiInFTBTI5LgNPMrJ2ZxQDXAbNKHDML+J35nA8cLOt+o4iIyMkQtGZV51yBmf0emIdvKMck59z3ZnaHf/84YC6+YRzr8Q3luDmAS08IUsi1gcquclRulaeyqxyVW+WdkLKz6tqTUkREJFSC2awqIiJSIyk5ioiIlFCtkqOZXWFma8xsvZn9pZT9Zmb/9O9faWbdAz03nFWx3Dab2SozW36iukDXJAGU3Rlm9rWZ5ZrZfRU5N5xVsdz0mSu77K73/z9daWZfmVnXQM8NZ1Ust4p/5kpOaRaqB75OOxuA9kAMsAI4s8QxVwH/wTc+8nzAE+i54fqoSrn5920GEkL9Pqpx2TUDegBPAfdV5NxwfVSl3Pz79Jkru+x6A439z6/U91zVys3/usKfuepUczw63ZxzLg8omm6uuKPTzTnnFgPxZtYywHPDVVXKrbYrt+ycc2nOuSVAfkXPDWNVKbfaLpCy+8o5t9//cjG+8d8BnRvGqlJulVKdkuPxppIL5JhAzg1XVSk38M1I9KmZfeufpq82qcrnRp+5n1T0vesz95Pyyu4WfK0+lTk3nFSl3KASn7nqtJ5jVaabC2gaujBV1Wn6LnTO7TSzZsD/mdmPzrmFJzTC6qsqnxt95o5Vkfeuz9yxSi07M+uP70v+ooqeG4aqUm5Qic9cdao5VmW6udo8DV2VpulzzhX9TANm4mu+qC2q8rnRZ+4nFXrv+syVX3ZmdjbwBjDYOZdRkXPDVFXKrVKfueqUHKsy3Vwg54arSpebmcWaWQMAM4sFBgK1acWTqnxu9JmrxHvXZ678sjOzJOAD4Ebn3NqKnBvGKl1ulf3MVZtmVVeF6eaOd24I3sZJV5VyA5oDM80MfJ+Ft5xzn5zktxAygZSdmbUAlgINAa+ZjcbXSy5Tn7mKlxu+pZj0mSv7/+vDQFPgVX85FTjnztP3XOXKjUp+z2n6OBERkRKqU7OqiIhItaDkKCIiUoKSo4iISAlKjiIiIiUoOYqIiJSg5ChSgpkV+mfv/87M3jOz+kH6PeeZ2T/9z/uZWe9KXGO0mf3O//wMf9zLzKxDFWM7x8yuKvb6msquAmFmiWZWa4ZrSHhQchT5uRzn3DnOubOAPOCOQE4yswqNG3bOLXXO3eN/2Q/fqgIB8/++/we85d80BPjIOdfNObeh2HFmZhX9v34OvrGxRbHOcs49W8FrFJ27F9hlZhdW5nyRUFByFCnbIqCjmTUxsw/9a8Ut9k9ThZk9amYTzOxTYJqZnWJmn/uP+9w/awdm9mt/TXSFmS30b+tnZh+b2an4EvAf/TW/Pma2ycyi/cc1NN96dNElYrsESPUPkL4KGA2MNLP5Znaqma02s1eBVKCtmb1mZkvN7Hsze6zoImbWw3zr360ws2/MrBHwOPAbfzy/MbMRZvYv//HHe49TzLdu6FdmttHMhhWL9UPg+hP47yISVEqOIsfhr5ldCawCHgOWOefOBh4EphU79Fx8czn+FvgXvuXBzgbeBP7pP+Zh4HLnXFfgmuK/xzm3GRgHvOivsS4CUoCr/YdcByQ750ou/3Qh8K3/GnOLXaO/f38nfyzdnHNbgL/6Zww5G+hrZmf7p+J6F/iDP7YBQLY/3nf98bxb4vce7z0CtMQ34fMvgeI1zaVAH0RqCCVHkZ+rZ2bL8X2hbwUm4vvCnw7gnPsCaOqvYQHMcs7l+J9fwE/NnNP5aWWA/wJTzOxWfNNflecNfprm72ZgcinHtAT2lnGNLf71O4sMN7NUYBnQGd90bp2AXf61F3HOZTrnCsqJ7XjvEeBD55zXOfcDvmm7iqQBrcq5rki1UW3mVhWpRnKcc+cU32D+iRlLKJp7MbuMazkA59wdZtYLX21wuZmdU8Y5OOf+628a7QtEOudKmyg5B6hbxmWOxmVm7YD7gB7Ouf1mNsV/rlH1ZY+Kn59b7HnxMqvrj1ekRlDNUSQwC/HfMzOzfkC6cy6zlOO+wtcMiv/4L/3ndHDOeZxzDwPpHLv8DsAhoEGJbdOAtym91giwGugYYPwN8SXLg2bWHF9zMcCPQCsz6+GPs4G/Obm0eIqU+h7LcTq1a/UNqeGUHEUC8yhwnpmtxHcv7abjHHcPcLP/uBuBP/i3jzGzVWb2Hb5Eu6LEebOBa4s65Pi3vQk0xpcgS/Mf4OJAgnfOrcDXnPo9MAlfMy/OuTzgN8ArZrYC+D98tbz5wJlFHXICfI9l6Q/MCSRWkepAq3KIVFP+3p6DnXM3lnHMTODPzrl1Jy+yivP30B3snNsf6lhEAqHkKFINmdkr+Jo+ryqx4G3J4zoBzZ1zC09acBVkZonAhc65D0Mdi0iglBxFRERK0D1HERGREpQcRURESlByFBERKUHJUUREpAQlRxERkRL+P6wjVHhymiCGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "por_values = np.linspace(0.0,0.25,100) \n",
    "fit_mean, fit_stdev = norm.fit(por,loc = por_average, scale = por_std) # fit MLE of the distribution \n",
    "cumul_p = norm.cdf(por_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(por_sort, p, c = 'red', edgecolors = 'black', s = 10, alpha = 0.7)\n",
    "plt.plot(por_values,cumul_p, c = 'black')\n",
    "plt.xlabel('Porosity (fraction)'); plt.ylabel('Cumulative Probability'); plt.grid(); \n",
    "plt.title('Nonparametric Porosity 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": [
    "#### 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 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}