GEO4310-Lecture04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lecture 4: Normal distribution and other continuous distributions\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outline\n",
    "* Continuous probability distributions: the basics\n",
    "    * Moments\n",
    "* Normal distribution\n",
    "    * Standard normal distribution\n",
    "    * 68-95-99 rule\n",
    "    * Reproducible properties\n",
    "    * Constructing pdfs to data\n",
    "    * Central limit theorem\n",
    "* Normal approximations to other distributions\n",
    "    * Binomial distribution\n",
    "    * Poisson distribution\n",
    "* Lognormal distribution\n",
    "* Chi square distribution\n",
    "* Student t distribution\n",
    "* F distribution\n",
    "* Extreme value distributions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pretty printing has been turned ON\n"
     ]
    }
   ],
   "source": [
    "## Import packages needed to do calculations and plotting\n",
    "# required to keep the plots in the notebook, not in an external window\n",
    "%matplotlib inline\n",
    "%pprint\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "from scipy import stats                 # required to use norm.pdf, norm.cdf\n",
    "import matplotlib.pyplot as plt         # required to plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous probability distributions: the basics\n",
    "###### From discrete to continuous distributions\n",
    "You've talked about discrete probability distributions. Some variables may be measured infinitely detailed. In that case, you may reduce the bin width and get more and more bins. This shape approximates a curve. Some of the curves are given by mathematical functions, which we'll talk about today.\n",
    "\n",
    "###### Properties of probability distributions\n",
    "* The area under the curve = 1\n",
    "* $P(a \\leq X \\leq b)$ is the area under the curve from $a$ to $b$\n",
    "* The function never evaluates negative, that is $f(x) \\geq 0$\n",
    "\n",
    "The first point, can also be stated as \"the integral of the probability density function equals one\", or \n",
    "\n",
    "$F(x) = \\int_{-\\infty}^{\\infty} f(x) dx = 1$. \n",
    "\n",
    "**note** Use this to solve question 3a in exercise 1.\n",
    "\n",
    "###### How do we characterize large datasets? Let's look at a dataset from all possible angles.\n",
    "\n",
    "###### Measures of central tendency\n",
    "* Mean (also called average, or expected value)\n",
    "* Median\n",
    "* Mode\n",
    "\n",
    "How should we define these measures mathematically to solve question 3d, 3e and 3f in exercise 1?\n",
    "\n",
    "###### Measures of spread\n",
    "* Variance\n",
    "* Standard deviation\n",
    "* Coefficient of variation\n",
    "\n",
    "###### Measure of skewness\n",
    "* Skewness\n",
    "\n",
    "__Learning goal 1__: understand what these measures are and where they are marked off the curve.\n",
    "\n",
    "__Learning goal 2__: understand what happens if the mean, standard deviation or skewness is changed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "Assume a normal distribution (precipitation) as an example of a continuous distribution. Where would these numbers be marked off the curve? (e.g. outside expected bounds)\n",
    "\n",
    "Next, find the numbers for the daily Uppsala precipitation (in EX1 - it is done for monthly data), and try to sketch the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean =  2.063778939507092\n",
      "Median =  0.0\n",
      "Mode =  0    0.0\n",
      "dtype: float64\n",
      "Variance =  25.120308724028334\n",
      "Stdev =  5.012016432936781\n",
      "Cv =  2.428562641565593\n",
      "Skewness =  8.15401600952\n"
     ]
    }
   ],
   "source": [
    "precip_daily=pd.read_table(\"../data/P_Uppsala.dat\",header=None,names=['year','month','day','prec'])\n",
    "P_daily = precip_daily.as_matrix()\n",
    "type(P_daily)\n",
    "\n",
    "# Measures of central tendency\n",
    "print('Mean = '  , precip_daily['prec'].mean())    # The balance point of the dataset. Influenced by outliers. \n",
    "print('Median = ', precip_daily['prec'].median())  # The middle value of the dataset, _sorted_ from smallest to largest. 50%.\n",
    "print('Mode = '  , precip_daily['prec'].mode())    # The most common values in the dataset.\n",
    "\n",
    "# Measures of spread\n",
    "print('Variance = ', precip_daily['prec'].var())     # Variance is the average squared deviation in the dataset.\n",
    "print('Stdev = '   , precip_daily['prec'].std())     # Standard deviation is the square root of the variance. The units of std and mean are the same.\n",
    "print('Cv = '      , precip_daily['prec'].std()/precip_daily['prec'].mean())  # The coefficient of variation is std relative to the mean.  \n",
    "\n",
    "# Measure of skewness\n",
    "print('Skewness = ', precip_daily['prec'].skew())    # Skewness is the degree of asymmetry in the distribution.\n",
    "\n",
    "# Let's sketch the distribution of daily precipitation according to these numbers. (on the blackboard)\n",
    "\n",
    "# Daily observations of rainfall are known to have heavy tails, because it is possible to get a rainfall value that is \n",
    "# very much larger than the median/mode, but values that are much less than the median/mode do not occur. (The median\n",
    "# and mode of a rainfall station is 0 if more than half of the days do not experience rain. 1. sept 2015, an extreme\n",
    "# event resulted in 118 mm rainfall recorded in 24 hours in Postmyr and 104 mm in Lier.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DescribeResult(nobs=4017, minmax=(0.0, 136.59999999999999), mean=2.0637789395070949, variance=25.120308724028455, skewness=8.150970878096675, kurtosis=149.7806267084826)\n"
     ]
    }
   ],
   "source": [
    "# MATRIX FORM, instead of dataframe. Consider converting into matrix if time.\n",
    "# Measures of central tendency\n",
    "from scipy.stats import describe\n",
    "\n",
    "print(describe(P_daily[:,3]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The normal distribution\n",
    "\n",
    "How do we know whether a dataset follows the normal distribution, and if it is skewed, how can we find out what distribution it follows? Let's try the normal distribution and see if it fits our daily precipitation. First, some basics about the normal distribution.\n",
    "\n",
    "We'll start by generating some artificial data to demonstrate some characteristics of the Normal distribution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 500 normally distributed data having mean = 16 and standard deviation = 3 \n",
    "# (which is not unrealistic summer temperatures in Oslo).\n",
    "generated_data = np.random.normal(loc=16, scale=3, size=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([   3.,   11.,   20.,   41.,   49.,  100.,   92.,   89.,   54.,\n",
       "          25.,    9.,    7.]),\n",
       " array([  6.61616892,   8.15844482,   9.70072073,  11.24299663,\n",
       "         12.78527253,  14.32754844,  15.86982434,  17.41210024,\n",
       "         18.95437615,  20.49665205,  22.03892795,  23.58120386,  25.12347976]),\n",
       " <a list of 12 Patch objects>)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADnFJREFUeJzt3X/sXXddx/Hny5WpDAKt/aaWbaUzaTDDKCzfEARcFjsV\nN7JO/1i2iKm6pCEZOPwRLJIA/5B0/iAao5jqJo0uG3MM27ChjAoh/sG0G2Nb12EHdNDZX0BkoAkw\nePvHPTM3X7732+6e7/3RT5+P5Jtzzud8zve8d76fvnbuufecm6pCktSuH5p1AZKkyTLoJalxBr0k\nNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY1bM+sCANavX1+bN2+edRmSdFZ58MEHv1pVC6fr\nNxdBv3nzZg4cODDrMiTprJLkqTPp56UbSWqcQS9JjTPoJalxBr0kNc6gl6TGnTbok9yW5GSSx4ba\n1iW5P8nhbrp2aN07kzyZ5PNJfmlShUuSzsyZnNF/EHjjkradwP6q2gLs75ZJcilwPfDKbpu/SnLe\nqlUrSXreThv0VfVp4OtLmrcBe7r5PcC1Q+13VtW3q+pLwJPAa1apVknSGMa9Rr+hqo5188eBDd38\nhcBXhvod7dokSTPS+87Yqqokz/sbxpPsAHYAbNq0qW8ZOgdt3nnvVPZzZNfVU9mPNCnjntGfSLIR\noJue7NqfBi4e6ndR1/YDqmp3VS1W1eLCwmkf1SBJGtO4Qb8P2N7Nbwf2DrVfn+SHk1wCbAH+vV+J\nkqQ+TnvpJskdwBXA+iRHgfcAu4C7ktwIPAVcB1BVB5PcBTwOPAvcVFXfm1DtkqQzcNqgr6obRqza\nOqL/+4D39SlKkrR6vDNWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BL\nUuMMeklqXO/n0Uutm8Zz733mvSbJM3pJapxBL0mNM+glqXEGvSQ1zqCXpMb5qRtpDvjJHk2SZ/SS\n1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN\nM+glqXG9gj7J7yQ5mOSxJHck+ZEk65Lcn+RwN127WsVKkp6/sYM+yYXAbwOLVfVTwHnA9cBOYH9V\nbQH2d8uSpBnpe+lmDfCjSdYALwT+C9gG7OnW7wGu7bkPSVIPYwd9VT0N/AnwZeAY8I2q+jiwoaqO\ndd2OAxt6VylJGlufSzdrGZy9XwK8DLggyZuH+1RVATVi+x1JDiQ5cOrUqXHLkCSdRp9LN1cCX6qq\nU1X1XeAe4HXAiSQbAbrpyeU2rqrdVbVYVYsLCws9ypAkraRP0H8ZeG2SFyYJsBU4BOwDtnd9tgN7\n+5UoSepjzbgbVtUDSe4GHgKeBT4L7AZeBNyV5EbgKeC61ShUkjSesYMeoKreA7xnSfO3GZzdS5Lm\ngHfGSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalx\nBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQ\nS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS43oFfZKXJrk7yRNJDiX52STr\nktyf5HA3XbtaxUqSnr++Z/R/DvxzVf0k8DPAIWAnsL+qtgD7u2VJ0oyMHfRJXgJcDtwKUFXfqar/\nBrYBe7pue4Br+xYpSRpfnzP6S4BTwN8l+WySv01yAbChqo51fY4DG/oWKUkaX5+gXwNcBnygql4N\n/A9LLtNUVQG13MZJdiQ5kOTAqVOnepQhSVpJn6A/Chytqge65bsZBP+JJBsBuunJ5Tauqt1VtVhV\niwsLCz3KkCStZOygr6rjwFeSvKJr2go8DuwDtndt24G9vSqUJPWypuf2bwNuT3I+8EXgNxn8z+Ou\nJDcCTwHX9dyHJKmHXkFfVQ8Di8us2trn90qSVo93xkpS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG\nGfSS1Li+N0xJOkts3nnvVPZzZNfVU9mPzpxn9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0k\nNc6gl6TGGfSS1DjvjNVETOsuTEmn5xm9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+gl\nqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWpc76BPcl6Szyb5aLe8Lsn9SQ5307X9y5QkjWs1\nzuhvBg4NLe8E9lfVFmB/tyxJmpFeXzyS5CLgauB9wO92zduAK7r5PcCngD/osx+tLr8URDq39D2j\n/zPgHcD3h9o2VNWxbv44sKHnPiRJPYwd9EneBJysqgdH9amqAmrE9juSHEhy4NSpU+OWIUk6jT5n\n9K8HrklyBLgT+Pkk/wCcSLIRoJueXG7jqtpdVYtVtbiwsNCjDEnSSsYO+qp6Z1VdVFWbgeuBf62q\nNwP7gO1dt+3A3t5VSpLGNonP0e8CfiHJYeDKblmSNCO9PnXznKr6FINP11BVXwO2rsbvlST1552x\nktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9J\njTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4\ng16SGrdm1gVIasvmnfdOfB9Hdl098X20xDN6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Lixgz7J\nxUk+meTxJAeT3Ny1r0tyf5LD3XTt6pUrSXq++pzRPwv8XlVdCrwWuCnJpcBOYH9VbQH2d8uSpBkZ\nO+ir6lhVPdTNfxM4BFwIbAP2dN32ANf2LVKSNL5VuTM2yWbg1cADwIaqOtatOg5sGLHNDmAHwKZN\nm1ajjLPeNO4olHTu6f1mbJIXAR8G3l5Vzwyvq6oCarntqmp3VS1W1eLCwkLfMiRJI/QK+iQvYBDy\nt1fVPV3ziSQbu/UbgZP9SpQk9dHnUzcBbgUOVdX7h1btA7Z389uBveOXJ0nqq881+tcDvw48muTh\nru0PgV3AXUluBJ4CrutXoiSpj7GDvqr+DciI1VvH/b2SpNXlnbGS1DiDXpIaZ9BLUuMMeklqnEEv\nSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNW5WvEjwX+DV/ks5WntFLUuMM\neklqnEEvSY3zGr2ks8603jM7suvqqexn0jyjl6TGGfSS1DiDXpIa5zV6SRphGu8FTON9AM/oJalx\nBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY1r4oYpvxREkkbzjF6SGjexoE/yxiSfT/Jk\nkp2T2o8kaWUTCfok5wF/CfwycClwQ5JLJ7EvSdLKJnVG/xrgyar6YlV9B7gT2DahfUmSVjCpoL8Q\n+MrQ8tGuTZI0ZTP71E2SHcCObvFbST4/gzLWA1+dwX6fD2tcHWdDjXB21GmNq2M98NXc0ut3vPxM\nOk0q6J8GLh5avqhr+39VtRvYPaH9n5EkB6pqcZY1nI41ro6zoUY4O+q0xtUxzRondenmP4AtSS5J\ncj5wPbBvQvuSJK1gImf0VfVskrcC/wKcB9xWVQcnsS9J0somdo2+qu4D7pvU718lM710dIascXWc\nDTXC2VGnNa6OqdWYqprWviRJM+AjECSpcc0HfZJXJHl46OeZJG9f0ueKJN8Y6vPuKdR1W5KTSR4b\naluX5P4kh7vp2hHbTuXxEiNq/OMkTyR5JMlHkrx0xLZHkjzaHc8DU67xvUmeHvp7XjVi21kexw8N\n1XckycMjtp3Wcbw4ySeTPJ7kYJKbu/a5GZMr1Dg3Y3KFGmc7JqvqnPlh8MbwceDlS9qvAD465Vou\nBy4DHhtq+yNgZze/E7hlxH/DF4CfAM4HPgdcOsUafxFY083fslyN3bojwPoZHcf3Ar9/BmNhZsdx\nyfo/Bd494+O4Ebism38x8J8MHl8yN2NyhRrnZkyuUONMx2TzZ/RLbAW+UFVPzbqQqvo08PUlzduA\nPd38HuDaZTad2uMllquxqj5eVc92i59hcI/EzIw4jmdipsfxOUkCXAfcMYl9n6mqOlZVD3Xz3wQO\nMbibfW7G5Kga52lMrnAcz8TEjuO5FvTXM/of1Ou6l34fS/LKaRY1ZENVHevmjwMblukzT4+X+C3g\nYyPWFfCJJA92d0FP29u6v+dtIy43zMtx/DngRFUdHrF+6scxyWbg1cADzOmYXFLjsLkZk8vUOLMx\nec4EfXfj1jXAPy6z+iFgU1X9NPAXwD9Ns7bl1OC13Nx+JCrJu4BngdtHdHlDVb2KwRNMb0py+dSK\ngw8wePn7KuAYg0sj8+oGVj6bn+pxTPIi4MPA26vqmeF18zImR9U4T2NymRpnOibPmaBn8Md9qKpO\nLF1RVc9U1be6+fuAFyRZP+0CgRNJNgJ005PL9Dnt4yUmLclvAG8Cfq37x/8DqurpbnoS+AiDl6VT\nUVUnqup7VfV94G9G7HsejuMa4FeBD43qM83jmOQFDMLp9qq6p2ueqzE5osa5GpPL1TjrMXkuBf3I\nM6ckP95dKyXJaxgcl69Nsbbn7AO2d/Pbgb3L9Jnp4yWSvBF4B3BNVf3viD4XJHnxc/MM3ix7bLm+\nE6px49Dir4zY9zw8puNK4ImqOrrcymkex2783wocqqr3D62amzE5qsZ5GpMr1DjbMTnJd6Dn5Qe4\ngEFwv2So7S3AW7r5twIHGbzL/RngdVOo6Q4GL+G+y+Ba3I3AjwH7gcPAJ4B1Xd+XAfcNbXsVg3fz\nvwC8a8o1PsngOuLD3c9fL62RwUvUz3U/B2dQ498DjwKPMPiHsnHejmPX/sHnxuBQ31kdxzcwuCzz\nyNDf9qp5GpMr1Dg3Y3KFGmc6Jr0zVpIady5dupGkc5JBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINe\nkhpn0EtS4/4PLkC8sShobnEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7acaabf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot a histogram of generated data.\n",
    "plt.hist(generated_data, bins = 12, normed = False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What have we plotted above? Describe the shape. <br/ >\n",
    "What names would you give to the axes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Example for the students to try at home:\n",
    "# Change the mean to 18\n",
    "generated_data_changed_mean = np.random.normal(loc=18, scale=3, size=500)\n",
    "# plt.hist(generated_data,              bins = 12, normed = False)               # plot the original dataset again\n",
    "# plt.hist(generated_data_changed_mean, bins = 12, normed = False, alpha = 0.7)  # then plot the data with changed mean, transparent\n",
    "# What is the effect of changing the mean? It is easier to see when a larger dataset is generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Example for the students to try at home, continued:\n",
    "# Change the standard deviation to 5 \n",
    "generated_data_changed_std = np.random.normal(loc=16, scale=5, size=500)\n",
    "# plt.hist(generated_data,             bins = 12, normed = False)                # plot the original dataset again\n",
    "# plt.hist(generated_data_changed_std, bins = 12, normed = False, alpha = 0.7)   # then plot the data with changed standard deviation\n",
    "# What is the effect of increasing the standard deviation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Before we plot normal distribution curve, let's generate a larger dataset\n",
    "generated_data = np.random.normal(loc=16, scale=3, size=10000)    # generate a larger sample of normally distributed data\n",
    "# plt.hist(generated_data,             bins = 50, normed = True)    # plot the generated data. 'normed' must equal True.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX59/HPNZN9AwJBkYCgBDDsNrIoLSpuiIorQuvS\n+muRWqu1/dWl0qfW/tQ+rb+2+hS1uKMgKKIgorHgUtyQsMgqgogQNgPREMieuZ4/ZjhMYiADTHJm\nud6vV16e+8w5mW/G4Zoz97nPfURVMcYYEz88bgcwxhjTuqzwG2NMnLHCb4wxccYKvzHGxBkr/MYY\nE2es8BtjTJyxwm+MMXHGCr8xxsQZK/zGGBNnEkLZSEQuAB4CvMATqvrnRo/3Bp4GTgXuVtUHGz3u\nBYqAbap6UXPP16FDB+3WrVtIf4AxxhhYunTpblXNCWXbZgt/oGhPBs4FioElIjJXVdcGbVYK3AJc\neohfcyuwDsgKJVS3bt0oKioKZVNjjDGAiHwV6rahdPUMBjaq6iZVrQFmAGOCN1DVr1V1CVDbRJhc\nYDTwRKihjDHGtJxQCn9nYGtQuziwLlT/AG4HfEewjzHGmBbSoid3ReQi4GtVXRrCthNEpEhEikpK\nSloyljHGxLVQCv82oEtQOzewLhRnAJeIyGb8XURni8jzTW2oqlNUtUBVC3JyQjo/YYwx5iiEUviX\nAHki0l1EkoBxwNxQfrmq3qWquaraLbDf26p6zVGnNcYYc8yaHdWjqnUicjNQiH8451OqukZEJgYe\nf0xEjsc/XDML8InIr4B8Vd3bgtmNMcYcBYnEO3AVFBSoDec0xpjQichSVS0IZVu7ctcYY+KMFX5j\njIkzIU3ZYEy86Xbn60e0/eY/j26hJMaEnx3xG2NMnLEjfmNCUFf2NRXr36diw2Kqt38GvnoAvBnt\nST35NN4Y4eG8887D6/W6nNSY5lnhN+YwfDWVlH3wAnuL5jjFPlj9vj3s+/RNLrzwTQYOHMijjz7K\n0KFDXUhqTOisq8eYQ6jaspLtT9zE3k9mN1n0G1uxYgXDhg3jpptuoqamphUSGnN07IjfmCbsX7eI\n3a//L9TXOeuSO59CRv9zST35NLzp7VBfPdXb1lHx+UfUrXmLyspKAB599FE2bNjA7NmzyczMdOtP\nMOaQ7IjfmEYeeeQRds/9i1P0PalZtB99G8f96C9k9D8Pb3o7AMTjJaVLX7JH/oy1a9dy8cUXO79j\nwYIFnH322diEgyYS2RG/MUGmT5/OL37xC6edkJ3LcVffS0JWx8Pu161bN+bMmcN9993H73//ewCK\niorIHTCc48Y/gCcxucH2NvzTuMmO+I0JWLx4MTfccIPTTurUi+Ov+UuzRf8AEWHSpEn861//wuPx\n/9Oq2fE5e954iEicGsXELyv8xgBbt25lzJgxVFdXA5DYvgvHXX0v3tSQ7hbawIQJE/jHP/7htCvW\n/YeyD2eELasxx8q6ekzc8/l8/OhHP2LXrl0AZGdnk3LF/8GTnB7y72h8pa9qdzIGjmLfijcAKHt/\nOild+5HSpW/4ghtzlOyI38S9hx56iEWLFgHg9XqZPXs2ie06HdPvFBGyz7mR5K79A2uUPa//HV9N\n5TGmNebYWeE3cW3dunXcddddTvvuu+9mxIgRYfnd4k2gw0W/dr451JXt4pt3ngzL7zbmWFjhN3Gr\nvr6e66+/3unXHzRoEHfffXdYnyMhswPtzp3otPeteJPKL5eH9TmMOVJW+E3cevrpp1myZAkASUlJ\nPPvssyQlJYX9edLzzySt5+lOu3TBv6itrQ378xgTKiv8Ji6VlZXxu9/9zmnfdddd9OvXr0WeS0TI\nPu/nSFIaAHWlxUyePLlFnsuYUFjhN3HpT3/6k3NVbZcuXbj99ttb9Pm86e1oe8Y4p33PPffYVb3G\nNTac08SF4OGWtaXb2P7kwXH2lYPGk3/vOy2eIfN7F1P+aSF1pdsoKytzLvYyprXZEb+JO98uet6Z\nbTM5tw9pvb/fKs8r3kSyz/6p037iiSfYsGFDqzy3McGs8Ju4UlOymYrPFjntdmf/FBFptedPPfk0\nUk4cCPgvHLv33ntb7bmNOSCkwi8iF4jIehHZKCJ3NvF4bxH5SESqReS/g9Z3EZF3RGStiKwRkVvD\nGd6YI1X2/nRnOTVvKMmd8lo9Q9vv/8hZnj59OuvWrWv1DCa+NVv4RcQLTAZGAfnAeBHJb7RZKXAL\n8GCj9XXAb1Q1HxgK/KKJfY1pFTW7vqDi8w+ddtszfuhKjuTOpzBq1CjAjvqNO0I54h8MbFTVTapa\nA8wAxgRvoKpfq+oSoLbR+h2quiywXA6sAzqHJbkxR+jboKP9tJ6nk3TcSa5l+eMf/+gsz5w5k9Wr\nV7uWxcSfUAp/Z2BrULuYoyjeItINGAQsPsTjE0SkSESKbJibCbea3Vuo3HjwrddmuDtH+wecdtpp\nXHTRRQCoKn/9619dzWPiS6uc3BWRDOBl4FequrepbVR1iqoWqGpBTk5Oa8QycWTvJ7Od5dS8oSTl\ndHMvTMCkSZOc5enTp1NcXOxiGhNPQin824AuQe3cwLqQiEgi/qI/TVVnN7e9MeG2fft29q9512ln\nDb7CvTBBhgwZwve/7x9KWldXx0MPPeRyIhMvQrmAawmQJyLd8Rf8cUBI35PFP07uSWCdqv7tqFMa\nE6LG8+IDfPPuM+Dz3z83ufMppOSe0sqpvutAzoqOZwL+4aX/+/BkJk2aRJs2bdwLZuJCs0f8qloH\n3AwU4j85+6KqrhGRiSIyEUBEjheRYuDXwCQRKRaRLOAM4FrgbBFZEfi5sMX+GmMa8VVXUB64GQpA\n1uDLXEzzXaknn0ZCdi4AWlPJlClTXE5k4kFIUzao6nxgfqN1jwUt78TfBdTY+0DrXR1jTCP7Vi1A\nq/cDkNDuBFJ7DHE5UUMiHrIGX0bpm/8PgIcffpjbbruNhASbTcW0HLty18QsVaV8+cHjlayCMYjH\n62KipmX0OQtPmr97p7i4mNdee83lRCbWWeE3Mavqq0+pK/WPlJGkNNL7nu1yoqZJQhIZA8532jZl\ns2lpVvhNzCpfNs9Zzug3Ek9SqotpDi9z4AUg/n+OCxcu5LPPPnM5kYllVvhNTKrb+zWVGz9x2pmD\nIntMQUJWR1J7DHbajzzyiItpTKyzwm9iUvnyN0B9AKScOIDE9l2a2cN9made5Cw/88wzlJeXu5jG\nxDIr/CbmaH0d+1b922kHF9RIlnLiAHr16gVAeXk5M2fOdDmRiVVW+E3MqfxiCb793wLgzchu0IUS\nyUSEiRMnOu0nnnjCxTQmllnhNzFn38q3nOX0fudE5BDOQ7nmmmtITEwEYPHixaxatcrlRCYWWeE3\nMaWufDeVm5Y67Yx+57qY5sh16NCByy47eHXxk08+6WIaE6us8JuYsn/VwqCTuv1JbNfJ5URH7qc/\nPXhf3ueee46qqioX05hYZIXfxAyfz9egmyej//mH2TpyjRw5khNPPBGA0tJSXn31VZcTmVhjhd/E\njPfee4+6sl0AeFIySOs5zOVER8fj8XDDDTc4bevuMeFmM0GZmDF16lRnOT1/BJKQ5GKao3Nguua6\nvV3wz2+oLFiwkG3bttG5s9211ISHHfGbmLB//35mzZrltNP7jnQxzbFLyOpIyon9Ai1l2rRpruYx\nscUKv4kJr7zyCvv27QMgITuXpOPzXE507NL7HPzwevbZZ1FVF9OYWGKF38SEZ5991lnO6DcS/83f\noltar9ORxGQA1q5dy7Jly1xOZGKFFX4T9YqLi1m4cGGgJaTnn+VqnnDxJKWS1vN0px18DsOYY2GF\n30S9adOmOd0gKScOICGrg8uJwif4XMX06dOpra11MY2JFVb4TdR7/vnnneX0vrFxtH9AStd+eDP9\nH2S7d+/mrbfeamYPY5pnhd9EtVWrVrF69WoAUlMbdo3EAvF4Sc8f4bSnT5/uYhoTK2wcv4lqwYVw\nzJgxfBTBd9k6Wun5I9i7+GUAXnjpZRYddzmepBQANv95tJvRTJQK6YhfRC4QkfUislFE7mzi8d4i\n8pGIVIvIfx/JvsYcLZ/P16Dw//CHP3QxTctJzOlOYvuuAGhtNZUbP3Y5kYl2zRZ+EfECk4FRQD4w\nXkTyG21WCtwCPHgU+xpzVD788EO2bNkCQHZ2NuefH51z8zRHRBp09+xf+56LaUwsCOWIfzCwUVU3\nqWoNMAMYE7yBqn6tqkuAxkMOmt3XmKMVfLR/1VVXkZQUfVM0hCotqPBXfrmM+ooyF9OYaBdK4e8M\nbA1qFwfWheJY9jXmkGpra3nxxReddqx28xyQ2PZ4kk/o7W/46qlY/4G7gUxUi5hRPSIyQUSKRKSo\npKTE7Tgmwi1cuJA9e/YAkJuby/Dhw11O1PKCj/r3r/uPi0lMtAul8G8DugS1cwPrQhHyvqo6RVUL\nVLUgJycnxF9v4tWMGTOc5bFjx+LxRMwxTItJ7zUcxP93Vm9dQ135HpcTmWgVyr+WJUCeiHQXkSRg\nHDA3xN9/LPsa06Tq6uoGNycZN26ci2lajzejHSld+wZaSsX6913NY6JXs+P4VbVORG4GCgEv8JSq\nrhGRiYHHHxOR44EiIAvwicivgHxV3dvUvi31x5j4UFhYSFmZ/+TmSSedREFBgcuJWk9a7x9Q9dVK\nACrWLXI5jYlWIV3AparzgfmN1j0WtLwTfzdOSPsacyxmzpzpLI8dOzYmZuIMVVrPYZS+9Qioj+rt\nn/HVV185t2k0JlSx3zFqYkpFRQVz5sxx2ldffbWLaVqfN60NKd0GOu3gkU3GhMoKv4kq8+fPZ//+\n/QD06tWLAQMGuJyo9aWf8gNnOfjbjzGhssJvokrwEe7VV18dV908B6TlDQWvv5d26dKlbNq0yeVE\nJtpY4TdRo6Kigtdff91pjx071sU07vGkZJDabZDTfumll1xMY6KRFX4TNebPn09FRQUAp5xyCn36\n9HE5kXvSeh+8YM0KvzlSVvhN1AgucFdddZWLSdyX1mMIeKy7xxwdK/wmKlRUVDBv3jynHe+F35OS\nQWr3g909s2bNcjGNiTZ2IxYTFd544w2nm6d3796Mfm4zIl+5nMpdab2GU/nFEsD/bej22293OZGJ\nFnbEb6JC426eeBzN01ha3hASExMBKCoq4ssvv3Q5kYkWVvhNxKusrLRuniZ4UjJI6HLwOoZTb/gT\n3e583fkx5lCs8JuIV1hY6Fy01bNnT/r27dvMHvEjeHRPxfoPXUxiookVfhPxXn75ZWfZunkaSu0x\nBDxeAGp2rKdur93LwjTPCr+JaNXV1cyde3Am7yuuuMLFNJHHm5pJStf+TrviczvqN82zwm8i2oIF\nC9i7dy/gn4J54MCBzewRf9J6neEs2y0ZTSis8JuIFjw+/corr7Runiak9Rx28M5cxevszlymWTaO\n30SsE387h+IZBwv/1B3HMdNGq3yHN60NyV36Ur1lJaBUbviIzFMvcjuWiWB2xG8iVtWWlfiq9gHg\nzcwhqVNPlxNFrvSg7p791t1jmmGF30Ss4P7qtF6nWzfPYaT2HAb4X5/qrWuoryhzN5CJaFb4TUSq\nr6+nYsPHTjv4BKb5roSMbJJzT/E31EfF5x+5G8hENCv8JiItWrQIX+Co1ZuRTXLn3i4ninxpPYNG\n99iwTnMYVvhNRAq+aCs1bxgi9lZtTlqvYc5y1Vef8s0337iYxkQy+9dkIo7P52P27NlOO63X6S6m\niR4JWR1J6pTnb/jqee2119wNZCJWSIVfRC4QkfUislFE7mzicRGRhwOPrxSRU4Meu01E1ojIahF5\nQURSwvkHmNizePFitm/fDoAnNYuULjY3T6iCz4UEf2syJlizhV9EvMBkYBSQD4wXkfxGm40C8gI/\nE4BHA/t2Bm4BClS1L+AFxoUtvYlJwRdtpeUNRQJz0ZjmBffzFxYWUl5e7mIaE6lCOeIfDGxU1U2q\nWgPMAMY02mYMMFX9PgbaikinwGMJQKqIJABpwPYwZTcxSFUbHKnaaJ4jk9iuE4kduwP+eY6Cb05v\nzAGhFP7OwNagdnFgXbPbqOo24EFgC7ADKFPVt5p6EhGZICJFIlJUUmIzDMarZcuW8dVX/jtrSXI6\nKSf2b2YP01haz4PnRKy7xzSlRU/uikg7/N8GugMnAOkick1T26rqFFUtUNWCnJycloxlIliDo/0e\ngxFvootpolPwt6T58+c7t6w05oBQCv82oEtQOzewLpRtzgG+VNUSVa0FZgM2RMM0ybp5wiOpQ1cS\nsnMB/03qCwsLXU5kIk0ohX8JkCci3UUkCf/J2bmNtpkLXBcY3TMUf5fODvxdPENFJE3819uPBNaF\nMb+JIWvWrOHzzz8HID09nZRug1xOFL1sdI85nGYLv6rWATcDhfiL9ouqukZEJorIxMBm84FNwEbg\nceCmwL6LgVnAMmBV4PmmhPuPMLEheDTP6NGj8SQmu5gmuqUHXfvw2muvUV1d7WIaE2lEVd3O8B0F\nBQVaVFTkdgzTyvr168fq1asBmDlzJrcvS3c5UfRSVRJm3cqmTZsAmDdvHqNHj3Y5lWlJIrJUVQtC\n2dau3DUR4fPPP3eKfkpKChdeeKHLiaKbiDS4TaV195hgVviN67rd+TpDfnaf05YuA+n7P++5mCg2\nBBf+OXPmUFtb62IaE0ms8JuI0HDufRvNEw6DBw+mSxf/YLvS0lLeffdddwOZiGGF37iu9tud1Oz6\nwt/wJJDWY7C7gWJE4+6e4JPnJr5Z4Teuq1h/cO741O6D8CTbSd1wCS78r7zyCvX19S6mMZHCCr9x\nnXXztJzTTz+dTp3802aVlJSwaNEilxOZSGCF37hq69at1OxY7294vKT2GOJuoBjj8Xi4/PLLnbZ1\n9xiwwm9cFnzDlZSu/fGmZrqYJrZ0u/N1ut35OrP2HJxT8bFnX8Dn87mYykQCK/zGVS+99JKzbN08\nLSO5Sx88aW0AqN9Xyocf2v14450VfuOabdu28cEHgf598ZDWc9jhdzBHRTxe0vIOvrbBH7YmPlnh\nN64Jvpo0pWs/vIGjUhN+ab2HO8uzZs2y7p44Z4XfuKZBN09QYTLhl9K1H57ULAC2b9/ORx995HIi\n4yYr/MYV27dvb9jNk2fdPC1JPN4GXWk2uie+WeE3rpg9ezYHZoZN6doXb3pblxPFvrRe1t1j/Kzw\nG1c0HM1j3TytIbi7p7i4mMWLF7ucyLjFCr9pdTt27HCuIPV4bDRPaxFvAml5Q532iy++6GIa4yYr\n/KbVzZo1y+nmGTFiBN70di4nih/BJ9Ffeukl6+6JU1b4TaubOXOms3z11Ve7mCT+pJw4gPbt2wP+\n6yjsYq74ZIXftKqtW7c6o3m8Xm+DeWRMyxOPt8GMncEfwiZ+WOE3rSp4GOHZZ59NTk6Oi2niU/C3\nrFmzZtlUzXHICr9pVdbN474RI0bQsWNHAHbu3GlTNcehkAq/iFwgIutFZKOI3NnE4yIiDwceXyki\npwY91lZEZonIZyKyTkRsCEec2rx5szOEMCEhgcsuu8zlRPHJ6/Vy5ZVXOm0b3RN/mi38IuIFJgOj\ngHxgvIjkN9psFJAX+JkAPBr02EPAm6raGxgArAtDbhOFggvMueeeS3Z2totp4lvj7p66ujoX05jW\nFsoR/2Bgo6puUtUaYAYwptE2Y4Cp6vcx0FZEOolIG+AHwJMAqlqjqt+GMb+JIi+88IKzbN087ho+\nfDgnnHAC4L8z18KFC11OZFpTKIW/M7A1qF0cWBfKNt2BEuBpEVkuIk+IiN1QNQ599tlnrFixAoDk\n5GTr5nGZx+Np8OEb/KFsYl9Ln9xNAE4FHlXVQcB+4DvnCABEZIKIFIlIUUlJSQvHMq0tuLCMHj2a\nrKwsF9MYgPHjxzvLr7zyClVVVS6mMa0plMK/DegS1M4NrAtlm2KgWFUPTAoyC/8HwXeo6hRVLVDV\nAhviF1tUtUHhDy44xj0FBQWcfPLJAOzdu5f58+e7nMi0llAK/xIgT0S6i0gSMA6Y22ibucB1gdE9\nQ4EyVd2hqjuBrSLSK7DdSGBtuMKb6LBs2TI2bNgAQGZmJqNHj3Y5kQEQkQYfwtbdEz+aLfyqWgfc\nDBTiH5HzoqquEZGJIjIxsNl8YBOwEXgcuCnoV/wSmCYiK4GBwP1hzG+iQHBBufTSS0lNTXUxjQkW\nXPjnzZvH3r17XUxjWktCKBup6nz8xT143WNBywr84hD7rgAKjiGjiWI+n6/BRVvWzRNZ8vPz6d+/\nPytXrqSqqopXX32V6667zu1YpoWFVPiNOVrvvfcexcXFAHTo0IFzzjnH5USm252vN2iXtTsVWAnA\ntGnTrPDHASv8pkU9//zzzvK4ceNITEx0MY1pSnr+CL597xkA3vr3AnJvfo6EjIMX123+s52TiTU2\nV49pMVVVVQ0mZbvmmmtcTGMOJSErh+Su/fwN9VGx7j/uBjItzgq/aTHBJwt79OjB4MGDXU5kDiU9\n/yxnef/ad90LYlqFFX7TYoK7ea655hpExMU05nDSe50OXn83XM3OjdTu3trMHiaaWR+/aRF79uxp\ncEHQj370o++cVDSRw5OSQdrJp1Hxuf+OXPvWvku7H1zrcirTUuyI37SImTNnUltbC8CQIUPo0aOH\ny4lMc9L7BHX3rHkbVbsfb6yywm9axLPPPussX3utHTlGg9STC/CkZAJQv7eE6i2rXU5kWooVfhN2\n69at45NPPgEgKSmJcePGuZzIhEK8iaTnj3Da+1YvcDGNaUlW+E3YBR/tX3zxxbRv397FNOZIpPcd\n6SxXrP8QX02li2lMS7GTuyas6uvreXDy4077P/Sxk7pRJOn4HiR26Ert7i1obRUV6z8Armx2PxNd\n7IjfhNWCBQuo31cKgCetLandm5yF20QoEWlw1L9vtd2ZKxZZ4Tdh9cwzzzjL6fkjEK99qYw26X3O\nAvGXhuotq/jyyy9dTmTCzQq/CZs9e/Ywe/Zsp53RzyZki0YJGdmkdB/ktJ9++mkX05iWYIXfhM20\nadOoqakBIOn4PJI6dnc5kTlaGf3Pc5afeuop6uvrXUxjws0KvwkLVeWJJ55w2hkDzjvM1ibSpfUY\njCetDQDbtm2jsLDQ5UQmnKzwm7AoKipi1apVAEhiMumnjGhmDxPJxJtIRtBJ3uAPdRP9rPCbsAgu\nDGm9huNJTnMxjQmHjP7nOsuvvfYau3btcjGNCScr/OaY7d+/v8F9da2bJzYktu9Ccud8AOrq6hpc\nmGeimxV+c8ymT59OeXk5AL169XKKhYl+wR/iU6ZMweezidtigRV+c0xUlUcffdRp33jjjTbvfgxJ\n6z2ctm3bAvDFF1+wYIHN3xMLrPCbY7JkyRKWL18OQEpKCtdff73LiUw4eRIb/j8N/pA30Sukwi8i\nF4jIehHZKCJ3NvG4iMjDgcdXisipjR73ishyEZkXruAmMgQXgnHjxpGdnX2YrU00mjhxorM8d+5c\niouLXUxjwqHZwi8iXmAyMArIB8aLSONO3FFAXuBnAtD4sOBWYN0xpzURpbS0lBkzZjjt4AJhYscF\nz3zh3Izd5/PRd9ztNvFelAvliH8wsFFVN6lqDTADGNNomzHAVPX7GGgrIp0ARCQXGA3YQOAY0+PK\n26mqqgIg6biTGTu7xApCjMoceKGzvO/TQrS+zsU05liFUvg7A8F3Xi4OrAt1m38AtwM2HCCG1NfX\nU778YJHPGDjKTurGsLSeQ/GmtwOgfl+pc29eE51a9OSuiFwEfK2qS0PYdoKIFIlIUUlJSUvGMmEw\nf/586r7dCfhv1J3e50x3A5kWJd5EMgZe4LTLl77mYhpzrEIp/NuALkHt3MC6ULY5A7hERDbj7yI6\nW0Seb+pJVHWKqhaoakFOTk6I8Y1bHnroIWc5Y8D5eBJTXExjWkPmwAvB459mu3rbOoqKilxOZI5W\nKIV/CZAnIt1FJAkYB8xttM1c4LrA6J6hQJmq7lDVu1Q1V1W7BfZ7W1WvCecfYFrfmjVrWLgwcIMO\n8ZB56mh3A5lW4c1oR/op33faDz/8sItpzLFotvCrah1wM1CIf2TOi6q6RkQmisiBYRzzgU3ARuBx\n4KYWymsiQPA/+LS8oSRkdXQxjWlNmd+7xFmeMWMGO3fudDGNOVoh9fGr6nxV7amqJ6vqfYF1j6nq\nY4FlVdVfBB7vp6rf+Q6oqu+q6kXhjW9a2+7du3nuueecdmbBJYfZ2sSa5E55JHc+BYDa2loeeeQR\nlxOZo2FX7pojMnnyZCorKwH/EM7k3D4uJzKtLbPg4GjuyZMns3//fhfTmKNhhd+ErKKign/+859O\nO2vw5TaEMw6l9RxGQpvjAP9FfHZrxuhjd8I2zTpwUVb5stcp3b0bAG9WR9J6D3czlnGJeLxkDb6M\n0n8/BsCvJv0Pf93SFfF4Adj8ZzvZH+nsiN+ERH317F3yitPOOu1S5x+6iT/p/c7Bk5oFQH3ZLirW\nf+ByInMkrPCbkFSs/7DBBVvBd2cy8ceTmELmoINH9nsXv4yqupjIHAkr/KZZqj7KPprptDMGjcaT\nlOpiIhMJMr93EZKQDEDNri+o3GQXdEULK/ymWZUbFlNbshkASUwhy4ZwGsCb1oaMAec77bIPZthR\nf5Swwm8OS1Up+/Dg1MuZgy7Em9bGxUQmkmQNuRy8/jEiNTvWU7V5hcuJTCis8JvDeuONN6jZ9QUA\nkpBE1mmXuZzIRJKEzA5k9D94X96yD1+wo/4oYIXfHJKqcu+99zrtjAEX4M1o52IiE4naDL3y4ORt\nxWt55513XE5kmmOF3xzS3LlzWbx4sb/hTfB/rTemkYSsjmT0Pdtp33333XbUH+Gs8Jsm1dfXM2nS\nJKedOWg0CZkdXExkIlmb08c5ff0ff/wx8+bZ7bUjmV25axzBt03ct+Yd9qxeDfhH8rQZepVbsUwU\nSGjTkcyBo5wbtNx9992MHj0aj8eOLSOR/V8x36H1dZS9P81pZ512Kd70ti4mMtGgzbCxSOCGPKtW\nreKFF15wOZE5FCv85jvKlze8rWLWYBvJY5rnTW/XYObOSZMmUVVV5WIicyhW+E0D9ZXllH0w3Wln\nDR2LJzmizoXVAAAN60lEQVTdxUQmmrQZfBnt27cHYPPmzQ1u0WkihxV+00DZhzPwVe0DIKHt8WR9\n72KXE5lo4knJ4J577nHa9913H7t27XIvkGmSFX7jqC3dRvmyg6Mx2p75EyQh0cVEJhrdeOON9O7d\nG4Dy8nL+8Ic/uJzINGaF3wD+i7W+efsJ8NUDkJzbh7Sep7ucykSjvN+/xZ6+Vzvtf015nBN+Yjdm\njyRW+A0Ac+bMofKLJYGW0O7sn9rdtcxRSz2pgJRug/wN9bHnrUfw+XzuhjIOK/yGffv2ccsttzjt\njIHnk9wpz8VEJtqJCNnn3HhwArft63nyySddTmUOsMJvuPfee9m6dSsAnrQ2tB3xY3cDmZiQ2D6X\nNoOvcNp33HEHJSUlLiYyB4RU+EXkAhFZLyIbReTOJh4XEXk48PhKETk1sL6LiLwjImtFZI2I3Bru\nP8Acm+XLl/P3v//dabc78wa8KRkuJjKxJGvYWOfG7N988w233Xaby4kMhDBlg4h4gcnAuUAxsERE\n5qrq2qDNRgF5gZ8hwKOB/9YBv1HVZSKSCSwVkX832te4pKamhp/85CfU1dUBkNylL+lBk20Zc6w8\niclknzuRr2f9EYBp06bxVuVJpOUNcbaxm7O3vlCO+AcDG1V1k6rWADOAMY22GQNMVb+PgbYi0klV\nd6jqMgBVLQfWAZ3DmN8cg/vvv59PP/0UgJSUFNpf8Es7oWvCLvXk00jvc5bTLi38J/WV5S4mMqEU\n/s7A1qB2Md8t3s1uIyLdgEHA4qaeREQmiEiRiBRZP2DLW7FiBffdd5/Tvv/++0nMts9k0zLajZyA\nN91/L4f6/d/wzcIpLieKb61ycldEMoCXgV+p6t6mtlHVKapaoKoFOTk5rRErbu3fv5/x48c7XTxn\nnHFGg1E9xoSbNzWT7PNvdtr717zD/rXvuZgovoVS+LcBXYLauYF1IW0jIon4i/40VZ199FFNuNx2\n22189tlnAKSlpfHUU0/h9XpdTmViXVrekAZdPnsKJ1MbmAzQtK5QCv8SIE9EuotIEjAOmNtom7nA\ndYHRPUOBMlXdIf4O4yeBdar6t7AmN0flpZde4vHHH3fa//znP+nZs6eLiUw8yT735yS07QSA1lSw\ne+5fqa2tdTlV/Gm28KtqHXAzUIj/5OyLqrpGRCaKyMTAZvOBTcBG4HHgpsD6M4BrgbNFZEXg58Jw\n/xEmNOvXr+enP/2p0x43bhw//vGP3Qtk4o4nOY0Ol/wWPP5vmDU71vPb3/7W5VTxRyLx3pgFBQVa\nVFTkdoyYsnfvXoYMGeJ08XTv3p3ly5fTpk0bZ5vgO3AZ05LKFr/Mt+8+7bSnTp3Ktdde62Ki6Cci\nS1W1IJRt7daLMexAIVf1UfLK/VRu8Bd9SUimasRtDHjgfTfjmTiWNfhyaravp+LzDwGYMGECffr0\n4dRTT3U5WXywKRviwLfvTaVyw8dOu/2oX5J03EkuJjLxTkRof+GvSGzvHxNSVVXFJZdcQnFxscvJ\n4oMV/hhXvnw+exfPctqZBWNIzz/TvUDGBHiS08i5fBJZWVkAbNu2jVGjRlFWVuZysthnhT+GVWxY\nTOm/H3PaqSefRruzbnAxkTENJWZ3Zvbs2SQk+HudV69ezeWXX051dbXLyWKbFf4YVVhYSMmcB0D9\nc6AndcqjwyV3IB4br28iy8iRIxtM2fz2229z1VVXUVNT42Kq2GYnd2PQ22+/zaWXXgr1/itzE9oe\nT8cr/g+epBSXkxnTtOuuu46tW7cyadIkAF577TXa9TuLDmO+e7Bik7odOzvijzFvvvkmF198MVVV\nVQB4s3I4btz9zjwpxkSq3/3ud9xxxx1Ou+LzDymZ82e0zo78w80KfwyZMWMGF198MRUVFQB4M9pz\n3Lj7SWjT0eVkxjRPRHjggQcazNlf+flHfD3rHnzVFS4miz3W1RMDVJV//OMf/OY3v+HABXldu3al\n7vy7SWzXyeV0xhxe4wsHNfFssk7bzN4lrwBQ9dVKdr1wFzlX/J6EzA5uRIw5dsQf5aqrq/mv//ov\nfv3rXztFPz8/nw8++MCmWTZRSURoe9YNtB1xvbOuZtcX7Jz6a6q3r3cxWeywwh/FtmzZwllnncXT\nTx+89H3YsGH85z//ITc318VkxhwbEaHN0KvIvuCXIP4yVb+vlJ3T7+Txxx8nEqeaiSZW+KPUK6+8\nwoABA/joo4+cdddffz3vvPMO7du3dzGZMeGTOeB8Ol79JzwH7gNdX8uECRMYP368Xeh1DKzwR5nS\n0lKuv/56Lr/8cr799lsAvF4vDz74IE8//TTJyckuJzQmvFJPHMDx1/2dxA4nOutmzpxJv379ePPN\nN11MFr2s8EcJVWXGjBnk5+czdepUZ33Xrl1ZtGgRv/nNb+x+uSZmJbbrxPHX/Y2MARc467Zu3cqo\nUaO47rrr2LnTbuhyJGxa5iiwdOlSbr31Vj744IMG69NO+QHZ592E98DXYGPiwIMFlfz85z9n9+7d\nzrrMzEwmTZrErbfeGrffeo9kWmY74o9gq1at4sorr6SgoKBB0fdmZJNz+e/JueR2K/om7lx55ZWs\nXbuW8ePHO+vKy8u544476NGjB4899phN99AMO+KPMCfeMY/qravYu+RVKjcuAYL+/3gSyPrexbQ5\nYxye5HTXMhoTKSo3LeWbt5+gds/WButzc3O55ZZb+NnPfkbbtm1dSte6juSI3wp/hNizZw/Tpk3j\nv//n79SWbP7O46k9h9FuxI9tbL4xjWh9HeUr3qDsw5n4Kr5t8JgkppDW+/v8+19/ZOjQoTF9HswK\nf5T49ttvmTdvHi+++CKFhYVNfj1N7TGEtsN/SNJxJ7uQ0Jjo4aup8t9/4pPZ3/kAAMjLy2Ps2LFc\nddVV9O/fP+Y+BKzwR6i6ujqWLVvG22+/zT2PTKO6eK0zbXIwSUgmve9ZZBWMce5QZIwJjdbVsH/t\nu+xdMofa3V81uU1ubi4XXngh55xzDmeeeSY5OTmtnDL8rPBHAJ/Px+bNm/n0008pKipi8eLFfPLJ\nJ5SXlx9yn6Tj88jofy7p+SOsD9+YY6Sq1Gz/jH0r/83+zxahNZWH3PaUU05hyJAhDBkyhIEDB9Kv\nXz/S06Pr32DYC7+IXAA8BHiBJ1T1z40el8DjFwIVwI9VdVko+zYlmgu/z+dj5MiRLF269LBF/oCk\n43uQ1vN00noPJ7HdCa2Q0Jj446utpurLZez/7H2qNhXhq95/2O1FhB49evDAAw9wxRVXtFLKY3Mk\nhb/Z2TlFxAtMBs4FioElIjJXVdcGbTYKyAv8DAEeBYaEuG9M8Xg8lJSUHLLoezM7kNK1HyknDiC1\n+/fwZtg8+ca0NE9iMmk9h5HWcxjqq6d62zoqv1xO9ZZVVO/4HHx1DbZXVTZs2BCz1wSEMi3zYGCj\nqm4CEJEZwBgguHiPAaaq/+vDxyLSVkQ6Ad1C2DcmBE8tW6L+qWM9qVkkdexOUseTSOrUk+QTeuHN\nyom5k0rGRBPxeEnp0peULn0B8NVWUbPrC2q2r6d65xfUfv0ltaXFoD4mzt9NwvuvN/Mbo++uYKEU\n/s5A8CDZYvxH9c1t0znEfcOq8dzebmj7/Wtpd+YNeDOyrcgbE+E8iSmk5PYhJbePs07raqjZvQVv\niPP/h6vutNYHSMTciEVEJgATAs19ItJSE293AHY3u1V8sNfCz16Hg+y1OKjVXwv5v8e0+4nNb+IX\nSuHfBgSPKcwNrAtlm8QQ9gVAVacAU0LIc0xEpCjUEyCxzl4LP3sdDrLX4qBYfi1CmatnCZAnIt1F\nJAkYB8xttM1c4DrxGwqUqeqOEPc1xhjTipo94lfVOhG5GSjEPyTzKVVdIyITA48/BszHP5RzI/7h\nnD853L4t8pcYY4wJSUh9/Ko6H39xD173WNCyAr8IdV+XtXh3UhSx18LPXoeD7LU4KGZfi4i8ctcY\nY0zLsfn4jTEmzsRV4ReRzSKySkRWiEh0zglxFETkKRH5WkRWB63LFpF/i8iGwH/j4hLiQ7wW94jI\ntsD7YoWIXOhmxtYgIl1E5B0RWSsia0Tk1sD6uHtfHOa1iNn3RVx19YjIZqBAVeNqnLKI/ADYh//q\n6r6BdX8BSlX1zyJyJ9BOVe9wM2drOMRrcQ+wT1UfdDNbawpcWd9JVZeJSCawFLgU+DFx9r44zGsx\nlhh9X8TVEX+8UtX/AKWNVo8Bng0sP4v/jR7zDvFaxB1V3XFgIkVVLQfW4b/SPu7eF4d5LWJWvBV+\nBRaIyNLAlcLx7LjAtRYAO4Hj3AwTAX4pIisDXUEx370RTES6AYOAxcT5+6LRawEx+r6It8I/XFUH\n4p9N9BeBr/1xLzAcN376/L7rUeAkYCCwA/hfd+O0HhHJAF4GfqWqe4Mfi7f3RROvRcy+L+Kq8Kvq\ntsB/vwZewT/zaLzaFejbPNDH+bXLeVyjqrtUtV5VfcDjxMn7QkQS8Re6aao6O7A6Lt8XTb0Wsfy+\niJvCLyLpgRM3iEg6cB6w+vB7xbS5wPWB5euBOS5mcdWBQhdwGXHwvgjcPOlJYJ2q/i3oobh7Xxzq\ntYjl90XcjOoRkZPwH+WD/4rl6ap6n4uRWo2IvACciX+2wV3AH4BXgReBrsBXwFhVjfmTnod4Lc7E\n/3Vegc3AjUH93DFJRIYDi4BVwIEbP/8Of992XL0vDvNajCdG3xdxU/iNMcb4xU1XjzHGGD8r/MYY\nE2es8BtjTJyxwm+MMXHGCr8xxsQZK/zGGBNnrPAbY0ycscJvjDFx5v8DStggvpTPxmIAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7a9394be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the normal distribution curve. \n",
    "# To do this, we need to import matplotlib.mlab as mlab and use the function normpdf\n",
    "my_mean = 16\n",
    "my_std = 3\n",
    "\n",
    "generated_data = np.random.normal(loc=my_mean, scale=my_std, size=10000)# generate a larger sample of normally distributed data\n",
    "plt.hist(generated_data, bins = 50, normed = True)                      # plot the generated data. 'normed' must equal True.\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "plt.plot(x,stats.norm.pdf(x,my_mean,my_std), color=\"black\", linewidth=3)# plot the normal distribution curve\n",
    "\n",
    "# Does the curve fit the data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal distribution curve is given as:\n",
    "\n",
    "$$f(x) = \\frac{1}{\\sigma \\sqrt{2 \\pi}}\\exp(\\frac{-(x-\\mu)^2}{2 \\sigma^2})$$\n",
    "\n",
    "Which parameters does it require? How do they change the curve?\n",
    "\n",
    "$f(x)$ is called the probabiliby density function, or pdf.\n",
    "\n",
    "How do we read off any percentiles from this curve, for instance the 90th percentile?\n",
    "The answer is: we can't. The 90th percentile is the x value up to which the area under the curve is 0.9.\n",
    "\n",
    "We may plot the distribution in a different way to make it easier to read off percentiles: <br/ >\n",
    "The cumulative distribution function, or CDF.\n",
    "(The normal cdf is not given as a formula because the integral of $f(x)$ is not analytically solvable.)\n",
    "\n",
    "Learning goal 3: understand the difference between the pdf and the CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc7a8fe0470>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FeX9/vH3hxBAkLDIYlkEqhSkihvFWmx/VVFBKpRK\nBUqptZVIMWpAkcUqNIiiICIVFSwRY2WzBQqioFaxP0QpgaKCiqRAIYCApBggBLI83z8Sp0nYDnBO\n5pw59+u6uDjPzCRzz3WSm2HOLOacQ0REgqWK3wFERCT8VO4iIgGkchcRCSCVu4hIAKncRUQCSOUu\nIhJAKncRkQBSuYuIBJDKXUQkgKr6teIGDRq4li1b+rV6EZGYtHr16q+ccw1Ptpxv5d6yZUsyMzP9\nWr2ISEwys/+EspwOy4iIBJDKXUQkgFTuIiIBpHIXEQkglbuISACdtNzNLN3MdpvZuuPMNzObbGZZ\nZvaxmV0e/pgiInIqQtlznwF0OcH8rkDr0j/JwHNnHktERM7ESc9zd879w8xanmCRHkCGK3le34dm\nVtfMvuWc2xmmjCIScM45Dhw4wN69e8nJyWHfvn3s27eP3NxccnNz2b9/PwcPHuTgwYPk5eWRl5dH\nfn6+9+fw4cMcPnyY+vXrs3TpUr83JyqE4yKmpsC2MuPs0mlHlbuZJVOyd895550XhlWLSLRzzvHl\nl1+yadMmtmzZwtatW9m6dSs7duxg586d7Nq1i927d5Ofn3/G62rUqFEYEgdDpV6h6pybBkwD6NCh\ng57MLRIwO3fuZM2aNXz88cesW7eOzz//nA0bNnDw4MFKWf+RI0cqZT2xIBzlvh1oXmbcrHSaiARY\nYWEhmZmZLF++nPfff5+VK1eyc+fpH40966yzOOecc6hfvz716tWjbt26JCUlUadOHc4++2zOPvts\natWqRc2aNTnrrLOoWbMmNWrUoHr16t6fGjVqhHELY1s4yn0hkGJms4Erga91vF0kmLZu3crixYt5\n4403WLZsGfv37w/p6+rWrcv5559Pq1ataNGiBc2bN6dZs2Y0adKExo0b07hxY2rVqhXh9PHlpOVu\nZrOAHwMNzCwbGAUkAjjnngdeB24CsoA84PZIhRWRypeVlcXcuXN59dVXWbt27QmXrVmzJpdddhmX\nXnopF198Me3ataNt27Y0aNAAM6ukxAKhnS3T9yTzHXBX2BKJiO/279/P7NmzmTFjBitWrDjucs2b\nN+eaa66hU6dO/OAHP+DCCy8kISGhEpPK8fh2y18RiT5ffPEFf/zjH3nppZeOecilWrVqXHPNNfzk\nJz/hxhtv5IILLtAeeZRSuYsIq1atYty4ccyfP5+S/4z/T9WqVenatSu9e/fm5ptvJikpyaeUcipU\n7iJxbO3atTz88MMsWrToqHlt27blzjvvpF+/fjRseNIH/0iUUbmLxKHt27czcuRIMjIyjprXrVs3\nBg8ezLXXXqtDLjFM5S4SRwoKCnjyyScZM2YMeXl53nQzo3fv3owcOZKLL77Yx4QSLip3kTjx4Ycf\nMmDAANatK3+D15tvvpmxY8eq1ANG93MXCbgjR44wYsQIOnXqVK7YL774Yt5++20WLlyoYg8g7bmL\nBNjnn39O3759y118VLNmTcaMGcM999xD1aqqgKDSOysSUDNnziQ5ObncTbuuvfZapk+fTsuWLf0L\nJpVCh2VEAqagoICUlBT69evnFXuNGjWYNGkSb731loo9TmjPXSRAcnJy+PnPf84777zjTWvdujV/\n+ctfaN++vY/JpLJpz10kIDZt2sT3v//9csXeq1cvMjMzVexxSOUuEgCrV6/mqquuYuPGjd60tLQ0\n5s6dq9sFxCkdlhGJccuWLePmm2/mwIEDQMnx9T//+c/ccsstPicTP6ncRWLYm2++SY8ePbznj9ar\nV49FixbRqVMnn5OJ31TuIjHq9ddfp2fPnt5zQ5s0acJbb71Fu3btfE4m0UDlLhKD3nnnHX72s595\nxX7eeefxzjvvcP755/ucTKKFyl0kxqxYsYLu3btz+PBhAFq2bMmyZcto0aKFz8kkmuhsGZEYsn79\nerp16+ZdnNSsWTPeeecdFbscReUuEiOys7Pp0qUL+/btA6BRo0a8/fbbtGrVyudkEo1U7iIx4Ouv\nv6Zr165kZ2cDULt2bZYsWUKbNm18TibRSuUuEuUKCwvp06ePd7veqlWrMm/ePC677DKfk0k0U7mL\nRLmhQ4eyZMkSb5yenk7nzp19TCSxQOUuEsXS09OZNGmSN37wwQfp37+/j4kkVqjcRaLU6tWrGTRo\nkDfu2bMnaWlpPiaSWKJyF4lCe/fupVevXt657BdddBEZGRlUqaJfWQmNflJEokxxcTH9+/dny5Yt\nACQlJTFv3jzOPvtsf4NJTFG5i0SZiRMn8sYbb3jjjIwMWrdu7WMiiUUqd5EosmrVKkaMGOGNhw4d\nSo8ePXxMJLFK5S4SJXJzc+nTpw+FhYUAXHnllYwdO9bnVBKrVO4iUSI1NZVNmzYBJcfZZ82aRWJi\nos+pJFap3EWiwIIFC3jxxRe98dSpU3XPGDkjIZW7mXUxsw1mlmVmw48xv46ZLTKzj8xsvZndHv6o\nIsG0a9cuBgwY4I379u1Lnz59fEwkQXDScjezBGAK0BVoB/Q1s4qPerkL+NQ5dwnwY+BJM6sW5qwi\ngeOcY+DAgXz11VcANG3alClTpvicSoIglD33jkCWc26Tc+4IMBuo+PG9A2qbmQFnAzlAYViTigTQ\nq6++yoIFC7zxiy++SL169XxMJEERSrk3BbaVGWeXTivrGeBCYAfwCXCvc644LAlFAmrPnj2kpKR4\n4+TkZK6//nofE0mQhOsD1RuBtUAT4FLgGTNLqriQmSWbWaaZZe7ZsydMqxaJTampqXzze9CsWTOe\neOIJnxNJkIRS7tuB5mXGzUqnlXU7MM+VyAI2A20rfiPn3DTnXAfnXIeGDRuebmaRmLd06VJmzpzp\njadNm0adOnV8TCRBE0q5rwJam1mr0g9J+wALKyyzFbgOwMwaA22ATeEMKhIUeXl5/O53v/PG/fr1\no2vXrj4mkiCqerIFnHOFZpYCLAUSgHTn3HozG1g6/3lgDDDDzD4BDBjmnPsqgrlFYtYjjzzC5s2b\nAahXrx4TJ070OZEE0UnLHcA59zrweoVpz5d5vQO4IbzRRILn008/Zfz48d54/PjxNGrUyMdEElS6\nQlWkkjjnSElJ8e4d88Mf/pDbb9f1fhIZKneRSjJ37lzeffddABISEnj22Wf18A2JGP1kiVSCAwcO\ncN9993nju+++m4suusjHRBJ0KneRSvDoo4+yfXvJGcSNGzdm9OjR/gaSwFO5i0TY5s2by50R88QT\nT+icdok4lbtIhD3wwAPeg647duzIL3/5S58TSTxQuYtE0Hvvvcdf/vIXbzxp0iR9iCqVQj9lIhFS\nXFzMkCFDvPEvfvELrrrqKh8TSTxRuYtEyKxZs1izZg0ANWrUYNy4cT4nkniicheJgPz8fEaOHOmN\nBw8eTPPmzU/wFSLhpXIXiYA//vGPbN26FYCGDRsyfPhRT6cUiSiVu0iY5eTk8Oijj3rjUaNGkZR0\n1OMNRCJK5S4SZo8//jj79u0DoHXr1iQnJ/ucSOKRyl0kjLKzs5k8ebI3Hjt2LImJiT4mknilchcJ\no7S0NPLz8wG44oor6NWrl8+JJF6p3EXCZMOGDaSnp3vjcePGYWY+JpJ4pnIXCZPRo0dTVFQEwHXX\nXUfnzp19TiTxTOUuEgYff/wxs2fP9sZlz5YR8YPKXSQMRo0a5b3u3r07HTt29DGNiMpd5IxlZmay\nYMECb5yWluZjGpESKneRM/Twww97r2+99VYuueQSH9OIlFC5i5yBlStX8sYbbwBgZnrCkkQNlbvI\nGShb5n379uXCCy/0L4xIGSp3kdP0wQcfsGTJEgCqVKnCQw895HMikf9RuYucpj/84Q/e6759+9K2\nbVsf04iUp3IXOQ0rV65k6dKlgPbaJTqp3EVOw5gxY7zXffr0oU2bNj6mETmayl3kFK1Zs4bFixcD\nJWfIPPjggz4nEjmayl3kFD3yyCPe6169etGuXTsf04gcm8pd5BSsW7eO+fPne+Pf//73PqYROT6V\nu8gpKHtDsJ/+9Ke0b9/exzQix6dyFwlRVlYWc+bM8cY61i7RTOUuEqJx48ZRXFwMwA033ECHDh18\nTiRyfCGVu5l1MbMNZpZlZsOPs8yPzWytma03s/fCG1PEX9u2bSMjI8Mba69dol3Vky1gZgnAFOB6\nIBtYZWYLnXOfllmmLvAs0MU5t9XMGkUqsIgfxo8fT0FBAQBXX301P/rRj3xOJHJioey5dwSynHOb\nnHNHgNlAjwrL/AKY55zbCuCc2x3emCL+2bNnD3/605+88ciRI31MIxKaUMq9KbCtzDi7dFpZ3wHq\nmdkyM1ttZr861jcys2QzyzSzzD179pxeYpFKNnnyZA4dOgTAZZddRpcuXXxOJHJy4fpAtSpwBdAN\nuBF4yMy+U3Eh59w051wH51yHhg0bhmnVIpGTm5vLM888442HDx+OmfmYSCQ0Jz3mDmwHmpcZNyud\nVlY2sNc5dxA4aGb/AC4BvghLShGfTJ06lX379gHQunVrbrnlFp8TiYQmlD33VUBrM2tlZtWAPsDC\nCsv8DbjazKqaWU3gSuCz8EYVqVz5+fk89dRT3viBBx4gISHBx0QioTvpnrtzrtDMUoClQAKQ7pxb\nb2YDS+c/75z7zMyWAB8DxcCfnHPrIhlcJNJefvlldu7cCUCTJk3o37+/z4lEQhfKYRmcc68Dr1eY\n9nyF8XhgfPiiifinqKiI8eP/9+M8ePBgqlev7mMikVOjK1RFjmH+/Pls3LgRgLp165KcnOxzIpFT\no3IXqcA5x+OPP+6NBw0aRFJSko+JRE6dyl2kgnfffZfMzEwAqlevzj333ONzIpFTp3IXqeCJJ57w\nXt9+++00btzYxzQip0flLlLGRx99VO7B1/fdd5/PiUROj8pdpIyye+233HILF1xwgY9pRE6fyl2k\n1JYtW8o9jOOBBx7wMY3ImVG5i5R66qmnKCoqAuDaa6/VwzgkpqncRYC9e/eWu63v0KFDfUwjcuZU\n7iLAc889R15eHgDt27fnxhtv9DmRyJlRuUvcO3ToEJMnT/bGQ4cO1W19Jeap3CXuvfTSS3zz8Jjm\nzZvTu3dvnxOJnDmVu8S1oqIinnzySW88ePBgEhMTfUwkEh4qd4lrCxYsICsrCyi5Qdgdd9zhcyKR\n8FC5S9xyzpW7re+gQYOoXbu2j4lEwkflLnFr+fLlrFy5EoBq1apx9913+5xIJHxU7hK3yt5q4Lbb\nbuPcc8/1MY1IeKncJS59+umnvPbaawCYmW4QJoGjcpe4NGHCBO919+7dadOmjY9pRMJP5S5xZ8eO\nHfz5z3/2xrpBmASRyl3iztNPP01BQQEAnTp14gc/+IHPiUTCT+UucSU3N5fnn3/eG+sGYRJUKneJ\nK1OnTiU3NxeANm3acPPNN/ucSCQyVO4SNw4fPsykSZO88f3330+VKvoVkGDST7bEjZkzZ7Jjxw4A\nzj33XPr37+9zIpHIUblLXCguLi53q4HU1FSqV6/uYyKRyFK5S1x47bXX+OyzzwCoXbs2AwcO9DmR\nSGSp3CUuPP74497rgQMHUqdOHR/TiESeyl0Cb/ny5axYsQIouUFYamqqz4lEIk/lLoE3btw473X/\n/v1p0qSJj2lEKofKXQJt3bp1LF68GCi5QZguWpJ4oXKXQCt7W9+ePXvqBmESN0IqdzPrYmYbzCzL\nzIafYLnvmVmhmfUKX0SR07NlyxZmzpzpjYcNG+ZjGpHKddJyN7MEYArQFWgH9DWzdsdZ7nHgzXCH\nFDkdEyZMoKioCIBrrrmGjh07+pxIpPKEsufeEchyzm1yzh0BZgM9jrHc3cBfgd1hzCdyWnbt2sX0\n6dO98YgRI3xMI1L5Qin3psC2MuPs0mkeM2sK9ASeC180kdP39NNPk5+fD8Dll19O586dfU4kUrnC\n9YHqJGCYc674RAuZWbKZZZpZ5p49e8K0apHyvv76a6ZMmeKNR4wYgZn5mEik8lUNYZntQPMy42al\n08rqAMwu/QVqANxkZoXOuQVlF3LOTQOmAXTo0MGdbmiRE3n22We92/p+5zvfoWfPnj4nEql8oZT7\nKqC1mbWipNT7AL8ou4BzrtU3r81sBvBaxWIXqQx5eXk89dRT3nj48OEkJCT4mEjEHyctd+dcoZml\nAEuBBCDdObfezAaWzn/+hN9ApBK98MILfHPI77zzzuOXv/ylz4lE/BHKnjvOudeB1ytMO2apO+d+\nfeaxRE7d4cOHy93W94EHHiAxMdHHRCL+0RWqEhgZGRls317ycVDjxo35zW9+43MiEf+o3CUQCgoK\neOyxx7zxfffdx1lnneVjIhF/qdwlEGbOnMnmzZsBqF+/vh7GIXFP5S4xr6ioiLFjx3rjIUOGULt2\nbR8TifhP5S4xb86cOWzcuBGAunXrkpKS4nMiEf+p3CWmFRUV8cgjj3jje+65R4/QE0HlLjFu7ty5\n5R58fe+99/qcSCQ6qNwlZhUVFZGWluaN7733XurXr+9jIpHooXKXmDVnzhw+//xzAJKSkhgyZIjP\niUSih8pdYlLFvfbU1FTq1avnYyKR6KJyl5g0c+ZMNmzYAECdOnVITU31OZFIdFG5S8wpKChg9OjR\n3njw4MHaaxepQOUuMWfGjBls2rQJKLkadfDgwT4nEok+KneJKfn5+eWOtQ8bNoykpCQfE4lEJ5W7\nxJSpU6eSnZ0NlNz58a677vI5kUh0UrlLzNi/f3+5q1FHjhxJrVq1fEwkEr1U7hIzJk6cyFdffQVA\nixYtuPPOO31OJBK9VO4SE3bv3s2ECRO8cVpaGtWrV/cxkUh0U7lLTBg7diwHDhwA4Lvf/S79+vXz\nOZFIdFO5S9T797//zXPPPeeNx44dS0JCgo+JRKKfyl2i3ogRIygoKADg6quvpnv37j4nEol+KneJ\nah988AGvvvqqN54wYQJm5mMikdigcpeo5Zxj6NCh3vjWW2/lyiuv9DGRSOxQuUvUevXVV3n//fcB\nSExM5LHHHvM5kUjsULlLVDp06FC5vfaUlBS+/e1v+5hIJLao3CUqTZw4ka1btwLQoEEDHn74YZ8T\nicQWlbtEnR07dpQ7BDNmzBjq1q3rYyKR2KNyl6hz//33c/DgQQAuuugi7rjjDp8TicQelbtElWXL\nljFr1ixv/PTTT1O1alUfE4nEJpW7RI2CggJSUlK8ce/evbn22mt9TCQSu1TuEjWeeeYZ1q9fD0Ct\nWrXK3ShMRE6Nyl2iwrZt23jooYe88ahRo2jWrJmPiURim8pdosLdd9/tfYj63e9+l3vvvdfnRCKx\nLaRyN7MuZrbBzLLMbPgx5vczs4/N7BMzW2Fml4Q/qgTVggUL+Nvf/uaNp06dSrVq1XxMJBL7Tlru\nZpYATAG6Au2AvmbWrsJim4H/55y7GBgDTAt3UAmmffv2lXsOanJyMp06dfIxkUgwhLLn3hHIcs5t\ncs4dAWYDPcou4Jxb4Zz7b+nwQ0AHSyUk999/Pzt27ABKHng9btw4nxOJBEMo5d4U2FZmnF067Xh+\nC7xxrBlmlmxmmWaWuWfPntBTSiC9+eabTJ8+3Rs/++yz1KtXz8dEIsER1g9UzewaSsp92LHmO+em\nOec6OOc6NGzYMJyrlhiTm5vLgAEDvPHPf/5zfvazn/mYSCRYQrn0bzvQvMy4Wem0csysPfAnoKtz\nbm944klQpaamejcGO+ecc3jmmWd8TiQSLKHsua8CWptZKzOrBvQBFpZdwMzOA+YB/Z1zX4Q/pgTJ\nggULePHFF73xlClTaNSokY+JRILnpHvuzrlCM0sBlgIJQLpzbr2ZDSyd/zzwMHAO8GzpI9AKnXMd\nIhdbYtWuXbvKHY7p27cvvXv39jGRSDCZc86XFXfo0MFlZmb6sm7xR3FxMd26dWPJkiUANG3alE8+\n+UQfooqcAjNbHcrOs65QlUozceJEr9gBZsyYoWIXiRCVu1SKf/7zn4wYMcIbDxs2jM6dO/uYSCTY\nVO4ScXv37qV3794UFhYCcOWVVzJmzBifU4kEm8pdIqqoqIh+/fqxZcsWAOrUqcOsWbNITEz0N5hI\nwKncJaLS0tJYunSpN87IyKBVq1Y+JhKJDyp3iZi//vWvpKWleeORI0fSvXt3HxOJxA+Vu0TEv/71\nL371q195486dO5crehGJLD15WMJu586ddO/enby8PAAuuOAC5syZQ0JCgs/JjtZy+OJjTt8yrlsl\nJxEJL+25S1jl5uZy0003kZ2dDZR8gLpo0SLq16/vczKR+KI9dwmbI0eOcMstt7B27VoAEhISmDNn\nDm3btq2U9R9vL7wy1qE9fYk2KncJi6KiIm677Tbefvttb9q0adO48cYbfUx1+irjHwqRSFK5yxkr\nLi4mOTmZ2bNne9PS0tL4zW9+E5H1qXhFTk7lLmfEOcc999xDenq6N+2uu+7i97//vY+poocO44hf\nVO5y2oqLixk0aBBTp071pt1+++1MnjyZ0ls/xw39b0KijcpdTkthYSEDBgxgxowZ3rQ+ffrwwgsv\nUKVKeE7CUmGKnD6Vu5yyvLw8+vbty8KF/3sgV//+/UlPTz+tc9lV4iLhp3KXU/LVV1/RvXt3Pvjg\nA2/ab3/7W6ZOnRqVFylFKx2Ll0jTRUwSsvXr19OxY8dyxT58+HBeeOEFFbtIlFG5S0gWLlzIVVdd\nxebNmwEwMyZPnsxjjz0Wdx+eisQCHZaREyooKODBBx9k/Pjx3rRatWrxyiuv0KNHj1P6Xjq2LlJ5\nVO5yXJs2baJ///6sWLHCm9aiRQsWLlxI+/btj/t1KvHTp2PxEi4qdzmKc4709HRSU1M5cOCAN/2m\nm24iIyODc845x8d08UmlL6dKx9ylnKysLK6//nruuOMOr9gTEhJ49NFHWbRokYpdJEZoz12AknPX\nx48fz7hx48jPz/emt2nThpdffpnvfe97R32NDr+IRC/tuce5oqIiMjIyaNOmDaNHj/aKvUqVKgwe\nPJg1a9Ycs9hFJLppzz1OFRUVMW/ePEaNGsVnn31Wbt7ll1/OtGnTuOKKK3xKJ6HSsXg5HpV7nDl0\n6BCvvPIK48eP54svvig3r1GjRjz66KP8+te/1kVJIjFO5R4nsrKymDZtGtOnTycnJ6fcvNq1azNk\nyBCGDBlCUlLSUV+rY+uxR3v0onIPsD179jB//nxefvllli9fftT8OnXqMGjQIO677z6dBRMnVPrx\nQ+UeMFlZWSxevJiFCxfy3nvvUVRUdNQyCXUaU/uybtS+tAszi2syc/yHPiQVkUhSuccw5xxbtmxh\nxYoVLFu2jL///e/evV8qSkhIoFu3bgwYMIC7/lGMVdExdfmfEx160159bFK5x4iioiI2b97MJ598\nwkcffcTq1avJzMzkyy+/POHXVW/ajpptO1Hrwh/xUa16pCxHxS4SB0IqdzPrAjwNJAB/cs6NqzDf\nSuffBOQBv3bOrQlz1kBzzpGTk8P27dvZtm0b//nPf9i8eTObNm3iiy++YOPGjRw+fPik36dmzZpc\nd911vH+kBWed34GqtRtUQnoJMh2nj00nLXczSwCmANcD2cAqM1vonPu0zGJdgdalf64Eniv9O/Cc\ncxw5coT8/Hzy8/M5dOgQhw4dIi8vj4MHD3Lw4EH2799Pbm4uubm57Nu3j//+97/k5OSQk5PDnj17\n2L17N7t27eLIkSOnvH6rdhbVm7SletMLqdHyEqp/6zt8nJBI7Qhsq0hZKv3oFsqee0cgyzm3CcDM\nZgM9gLLl3gPIcM454EMzq2tm33LO7Qx74ihxww03sGzZMgoKCiptnVVq1aVagxYkNmxJtcbnU+3c\nC0is31SHWSSqnOqps/rHIDJCKfemwLYy42yO3is/1jJNgcCV+zc/uLu+2BX2YrfqtUioVY+qSQ1J\nqN2AqnUbk1j3XKrW/RaJ9ZtSpcbZYV2fSDQI13UU+keivEr9QNXMkoHk0uEBM9sQwdU1AL6K4PcP\nO3f4IIWHD1KYkx2Obxdz2x9G2vY4ZI/Hzba3CGWhUMp9O9C8zLhZ6bRTXQbn3DRgWijBzpSZZTrn\nOlTGuqJRPG+/tl3bLqHdFXIV0NrMWplZNaAPsLDCMguBX1mJ7wNfB/l4u4hItDvpnrtzrtDMUoCl\nlJwKme6cW29mA0vnPw+8TslpkFmUnAp5e+Qii4jIyYR0zN059zolBV522vNlXjvgrvBGO2OVcvgn\nisXz9mvb41M8b/tRrKSXRUQkSPQkJhGRAApkuZvZFjP7xMzWmlmm33kiyczSzWy3ma0rM62+mb1l\nZhtL/67nZ8ZIOs72jzaz7aXv/1ozu8nPjJFiZs3N7F0z+9TM1pvZvaXTA//+n2Db4+K9D0UgD8uY\n2Ragg3Mu8Oe8mtmPgAOUXCF8Uem0J4Ac59w4MxsO1HPODfMzZ6QcZ/tHAweccxP8zBZpZvYt4FvO\nuTVmVhtYDfwU+DUBf/9PsO23EgfvfSgCueceT5xz/wByKkzuAbxU+volSn7oA+k42x8XnHM7v7lB\nn3NuP/AZJVeGB/79P8G2S6mglrsD3jaz1aVXxcabxmWuM/gSaOxnGJ/cbWYflx62CdxhiYrMrCVw\nGbCSOHv/K2w7xNl7fzxBLfernXOXUnK3yrtK/+sel0pPUw3esbcTew74NnApJfc3etLfOJFlZmcD\nfwVSnXO5ZecF/f0/xrbH1Xt/IoEsd+fc9tK/dwPzKbmzZTzZVXpM8ptjk7t9zlOpnHO7nHNFzrli\n4AUC/P6bWSIl5faKc25e6eS4eP+Pte3x9N6fTODK3cxqlX7AgpnVAm4A1p34qwJnIXBb6evbgL/5\nmKXSfVNspXoS0Pe/9CE504HPnHMTy8wK/Pt/vG2Pl/c+FIE7W8bMvk3J3jqUXIE70zk31sdIEWVm\ns4AfU3I3wF3AKGABMBc4D/gPcKtzLpAfOh5n+39MyX/LHbAFuDOI9zoys6uB/w98AhSXTh5JybHn\nQL//J9iG/AYuAAAAPklEQVT2vsTBex+KwJW7iIgE8LCMiIio3EVEAknlLiISQCp3EZEAUrmLiASQ\nyl1EJIBU7iIiAaRyFxEJoP8D2HQPQrG3HlgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7a9274f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To match the data with the curve, we first need to find the cumulative values of our histogram.\n",
    "generated_data = np.random.normal(loc=my_mean, scale=my_std, size=10000)# generate a larger sample of normally distributed data\n",
    "plt.hist(generated_data,             bins = 50, normed = True, cumulative=False)          # plot the generated data. 'normed' must equal True.\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "plt.plot(x,stats.norm.cdf(x,my_mean,my_std), color=\"black\", linewidth=3)# plot the normal distribution curve\n",
    "\n",
    "# Does the curve fit the data? Try changing cumulative=False into cumulative=True in line \"plt.hist...\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "136.6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFaBJREFUeJzt3XuwXWddxvHvz4RSipRQiwHSYGsJtlWKVC0BLxwtSqhC\nHf6wVkUFRzuj8T5SizO21dERxetUsGLpOMxIZ8TqBKelVGWPN0gptKXQpCRAIBeaCrS0YLEJ+fnH\nWjldZ/ecs3f25bx7r/X9zJzpu/Zee51npTnPefOufYnMRJLUPl9TOoAkaToseElqKQteklrKgpek\nlrLgJamlLHhJaqmBBR8Rb4+IwxFxzyr7/EVE7ImIuyPiRZONKEkaxTAz+BuAbSvdGREXA8/LzC3A\nzwFvnVA2SdIYBhZ8Zv4H8OAqu7wa+Nt6353AhojYOJl4kqRRTWINfhOwv7F9ADhjAseVJI1hUhdZ\no2/b9z+QpMLWT+AYB4HNje0z6tuWiAhLX5JGkJn9k+ihTGIGvwP4SYCI2Ao8lJmHl9sxM+f266qr\nrjp+FvXXA8A5fWd4CvAa4A/rP5bDjf2z7/Gjfo16jKsmlKHEOVw15uNLn0N//nk7h+Xyz9M5rJR/\nds5hte4Zx8AZfES8E3gZcHpE7Kf603pSHeq6zLw5Ii6OiL3Al4HXjZVoiiJG+iW4jKNURb77+JGB\n3wDeCDx9Qt9DksYzsOAz87Ih9tk+mThrYdTfiFcD19TjNwP/WY8DeCdw6XixJGnCfCXr0Bbq/x4C\nfqdx++8wH+W+UDrAGBZKBxjTQukAY1ooHWBMC6UDFBNr9YEfEZGlP1ykWqIZJ0MA24Fr6+3zgQ8x\n/LXqcb//JI5R+vGzkMFzmI0MnsPxY6zWjREx8kXWuSn4ya2fj/s/8ynAo/X2LazyIt8VHu9f6PIZ\nPIfZyOA5HD/GtAp+Ek+TXEOT+MswruPl/q3AKyZwPEmaDtfgR/bzTOYXhiRNhwU/tN2N8cnMx4VV\nSV1mwQ/tHxvjVwKnlgoiSUOx4If2z43xa4qlkKRhzdmzaEpdcX8QOB04Vh/jgXp7rb7/JI9R+vGz\nkMFzmI0MnsPxY0zrWTTO4IfSoyp3gG9ntHKXpLVlwQ/lfY3xRcVSSNKJsOCH8p+N8cuKpZCkE2HB\nD/Ql4O7G9ktKBZGkE2LBD3QHj6+/g28HLGleWPAD3V46gCSNxIIf6EOlA0jSSNb0zcYeeuihtfx2\nE/Lh0gEkaSRrWvAbN5450uMyj042yNC+BOytx+uArxbKIUknbk0L/rHHRp3B7wS2TjLKkO5pjM8B\nPlYggySNxjX4VTUL/gXFUkjSKCz4VTVn7N9SLIUkjcKCX5UFL2l+WfCr2tUYn1cshSSNwoJf0cPA\noXp8EnBWwSySdOIs+BXd1xhvYe4+n1xS51nwK/p4Y/z8YikkaVQW/IoseEnzzYJf0Z7GeEuxFJI0\nKgt+RXsb4+cVSyFJo7LgV/SJxvjsYikkaVQW/LIeAr5Qj58MPKdgFkkajQW/rE81xmfhH5OkeWRz\nLatZ8N9YLIUkjcOCX1b/DF6S5o8Fv6x9jfGZhTJI0ngs+GXta4zPLJRBksYzsOAjYltE7I6IPRFx\nxTL3nx4R74mIuyLioxHx01NJuqY+3Rh/Q7EUkjSOVQs+ItYB1wLbqN4v97KIOLdvt+3AnZn5rcAC\n8McRMefvzGXBS5p/g2bwFwJ7M3NfZh4BbgQu6dvns8Cp9fhU4PNZ7lOyJ+CLVG8VDHAy8MyCWSRp\ndINm2puA/Y3tA8CL+/Z5G/BvEXEIeBrwI5OLV0LzdJ8LRKkgkjSWQQWfQxzjjcBdmbkQEWcDt0XE\nCzPzkSfuenVjvFB/zZrPNMabi6WQ1E29Xo9erzeRYw0q+IMsbbnNVLP4ppcCvweQmZ+IiE8B3wTc\n8cTDXT1izLXUnMFb8JLW1sLCAgsLC4vb11xzzcjHGrQGfwewJSLOjIiTgEuBHX377AZeDhARG6nK\n/ZMjJyqu+fvrucVSSNK4Vp3BZ+bRiNgO3AqsA67PzF0RcXl9/3XA7wM3RMTdVL8w3pCZX1jxoDOv\nOYM/o1gKSRrXwKczZuYtwC19t13XGH8OeNXko5XSnMFb8JLml69kfYJmwW8qlkKSxmXBP8HBxtgZ\nvKT5ZcEv8TDwpXp8MvCMglkkaTwW/BLN2fsmfJGTpHlmwS9xqDH2Y/okzTcLfolmwXuBVdJ8s+CX\ncAYvqT0s+CUseEntYcEvYcFLag8LfonPNsbPLpZCkibBgl/CgpfUHhb8osSCl9QmFvyiR4BH6/Ep\nVB9OJUnzy4Jf1Jy9PwtfxSpp3lnwi+5vjF2ekTT/LPhFzYJ/VrEUkjQpFvwiC15Su1jwiyx4Se1i\nwS9qFvzGYikkaVIs+EWHG2MLXtL8s+AXWfCS2sWCX+QSjaR2seCB6m0KHmhsW/CS5p8FD8CDwNF6\n/DTgKQWzSNJkWPCA6++S2siCByx4SW1kwQNLC/7ri6WQpEmy4AEvsEpqIwseWFrwzuAltYMFDziD\nl9RGFjzgDF5SG1nwgAUvqY0seMCCl9RGFjywtOCfWSyFJE2SBc//AQ/X4/XAMwpmkaTJGVjwEbEt\nInZHxJ6IuGKFfRYi4s6I+GhE9Caecqr+pzE+HX/nSWqL9avdGRHrgGuBlwMHgQ9GxI7M3NXYZwPw\nl8ArMvNARJw+zcCT5/q7pHYaNF29ENibmfsy8whwI3BJ3z4/BvxDZh4AyMzPTT7mNDVn8K6/S2qP\nQQW/Cdjf2D5Q39a0BTgtIt4XEXdExGsnGXD6vMAqqZ1WXaKh+iSMQZ4EXABcBJwCvD8iPpCZe8YN\ntzaaM3iXaCS1x6CCPwhsbmxvpprFN+0HPpeZjwKPRsS/Ay8Elin4qxvjhfqrNJdoJM2OXq9Hr9eb\nyLEic+VJekSsB+6jmp0fAm4HLuu7yHoO1YXYVwBPBnYCl2bmvX3HyuH+QbCcncBWRn/8YopljvEz\nwNvr8V8Bl5/g48f9/mt9jNKPn4UMnsNsZPAcjh9jQA+TmTHKkVedwWfm0YjYDtwKrAOuz8xdEXF5\nff91mbk7It4DfAQ4Brytv9xnmzN4Se206gx+ot9oZmfwLwE+UI//HfjuE3z8uN9/rY9R+vGzkMFz\nmI0MnsPxY0xrBu+rerzIKqmlLHiXaCS1VMcLvvk+NOuADQWzSNJkdbzgmy+69X1oJLVLxxvN5RlJ\n7WXBL7LgJbVLxwu+uURjwUtql44XfP97wUtSe1jwi5zBS2oXC36RBS+pXTpe8P1Pk5Sk9uh4wTuD\nl9ReFvwiC15Su1jwi1yikdQuHS74rwKfb2xb8JLapcMF/yCPv4/z06k+WlaS2qPDBe/6u6R263DB\n+zYFktqtwwXvBVZJ7WbBA87gJbVRhwveJRpJ7dbhgneJRlK7dbjgncFLarcOF7xr8JLazYIHXKKR\n1EYdLniXaCS1W0cLPnGJRlLbdbTgvwx8pR6fDDy1YBZJmo6OFnz/7D1KBZGkqbHgvcAqqaUseNff\nJbWUBe8MXlJLdbTgfYqkpPbraMG7RCOp/Sx4C15SS1nwFryklhpY8BGxLSJ2R8SeiLhilf2+IyKO\nRsRrJhtxGix4Se23asFHxDrgWmAbcB5wWUScu8J+bwLew1y8asiCl9R+g2bwFwJ7M3NfZh4BbgQu\nWWa/XwTexdLmnGHNmF9fLIUkTdOggt8E7G9sH6hvWxQRm6hK/631TTmxdFPxFeCRerwe2FAwiyRN\nz6CCH6as/wz4zcxMquWZGV+i6X+R04zHlaQRrR9w/0Fgc2N7M9UsvunbgBsjAqrGfGVEHMnMHU88\n3NWN8UL9tdZcnpE0u3q9Hr1ebyLHimrivcKdEeuB+4CLgEPA7cBlmblrhf1vAN6dmTctc1+Ovnqz\nE9jK+Ks/QXUdeFu9fRHwLyf4+HEyjPv4WcjgOcxGBs9hNjJM5hwG9DCZOdJSw6oz+Mw8GhHbgVuB\ndcD1mbkrIi6v779ulG9a1gONsTN4Se01aImGzLwFuKXvtmWLPTNfN6FcU+RTJCV1QwdfyeoMXlI3\ndLDgncFL6oYOFnxzBm/BS2qvjhf8xmIpJGnaOl7wrsFLai8LXpJaqoMF/7/1f58MPK1kEEmaqg4W\n/HHPxPehkdRmHS54L7BKajcLXpJayoKXpJbqcMH7DBpJ7dbhgncGL6ndLHhJaikLXpJayoKXpJbq\ncME/q3QASZqqDhX8VxrjdcBppYJI0proUMEfbow30qlTl9RJHWq5+xtjl2cktZ8FL0ktZcFLUkt1\ntOCfXSyFJK2VDhX8ZxtjZ/CS2q+jBe8MXlL7WfCS1FIdLXiXaCS1X0cK/hhLL7I+p1QQSVozHSn4\nzwFH6/FpwMkFs0jS2uhIwR9qjF1/l9QNHSn4g42xyzOSuqGDBb+pWApJWksdLHhn8JK6oSMF31yD\ndwYvqRs6UvAu0Ujqno4U/IHG+IxiKSRpLQ1V8BGxLSJ2R8SeiLhimft/PCLujoiPRMR/RcT5k486\nDgteUvcMLPiIWAdcC2wDzgMui4hz+3b7JPA9mXk+8LvAX0866Oi+DDzY2H5mqSCStKaGmcFfCOzN\nzH2ZeQS4EbikuUNmvj8zv1hv7mSmpskH+rY7siolqfOGabtNwP7G9gFWv1L5M8DN44SarP2Dd5Gk\nFlo/xD457MEi4nuB1wPfufweVzfGC/XXtFnwkuZHr9ej1+tN5FiRuXp/R8RW4OrM3FZvXwkcy8w3\n9e13PnATsC0z9y5znDyB3xV9dgJbGe3x17D0F8uoGQCi8ONnIYPnMBsZPIfZyDCZc1ithyOCzIxR\njjzMEs0dwJaIODMiTgIuBXb0BXguVbn/xHLlXtZnSgeQpCIGLtFk5tGI2A7cCqwDrs/MXRFxeX3/\ndcBvA88A3hoRAEcy88LpxT4Rny4dQJKKGLhEM7FvVGyJZgvQ/EfF/P9zznMo/fhZyOA5zEaG+V+i\nmWPHcIlGUle1vOA/CzxWj7+uZBBJWnMtL/hPNcZnFUshSSV0qODPLBVCkopoecF/sjF2Bi+pW1pe\n8J9ojM8ulkKSSrDgJamlLHhJaqkWF/zDwOF6fBLw3IJZJGnttbjg9zTGZ1O9y4IkdUdHCn5LsRSS\nVEqLC353Y3xOsRSSVIoFL0kt1eKCv7cxtuAldU9LC/4ocF9j+9xSQSSpmJYW/F4efxfJTcCGglkk\nqYyWFvxHGuNvKZZCkkrqQMG/sFgKSSqppQV/V2N8frEUklRSSwv+w43xBcVSSFJJLSz4Q1Qf1Qfw\nVOD5BbNIUjktLPjbG+Nvw/egkdRVLSz4DzTG31EshSSV1sKC/6/G+KXFUkhSaS0r+EdZukTznaWC\nSFJxLSv49/P4K1jPATYWzCJJZbWs4N/bGH9fsRSSNAtaVvC3NMY/UCyFJM2CFhX8Ph5/i4KTcAYv\nqetaVPB/3xhfBDytVBBJmgktKfgE3tHY/pFSQSRpZrSk4G8H7qnHTwFeUzCLJM2GlhT8HzfGlwKn\nlgoiSTOjBQV/N/CuxvavlAoiSTNlzgv+q8DPU63BA7wKP+BDkioDCz4itkXE7ojYExFXrLDPX9T3\n3x0RL5p8zJW8Efjverwe+IO1+9aSNONWLfiIWAdcC2wDzgMui4hz+/a5GHheZm4Bfg5465SyNjwG\n/Brwh43bfruOOC29KR57LfRKBxhDr3SAMfVKBxhTr3SAMfVKByhm0Az+QmBvZu7LzCPAjcAlffu8\nGvhbgMzcCWyIiCm9CcxHgTcD5wJ/2rj9B6lm89PUm/Lxp61XOsAYeqUDjKlXOsCYeqUDjKlXOkAx\n6wfcvwnY39g+ALx4iH3OAA6PnQ6o1tdfW49fsMz9r6b6veMHe0hS06CCzwH3HxfDPO7UU1815OGW\neuSRg+QTjvgM4BrgF5j7a8WSNAWRT2zOx++M2ApcnZnb6u0rgWOZ+abGPn8F9DLzxnp7N/CyzDzc\nd6xhf1lIkhoys38SPZRBM/g7gC0RcSbVp1lfClzWt88OYDtwY/0L4aH+ch8noCRpNKsWfGYejYjt\nwK1Ui9zXZ+auiLi8vv+6zLw5Ii6OiL3Al4HXTT21JGmgVZdoJEnza+pXJ4d5odQsiYjNEfG+iPhY\nRHw0In6pvv20iLgtIj4eEe+NiA2ls64mItZFxJ0R8e56e27yR8SGiHhXROyKiHsj4sVzlv/K+u/P\nPRHxdxHx5FnNHxFvj4jDEXFP47YVs9bntqf+mS7+qTor5P+j+u/O3RFxU0Q8vXHfzOdv3PfrEXEs\nIk5r3HZC+ada8MO8UGoGHQF+NTO/GdgK/EKd+TeB2zLz+cC/1tuz7JeBe3n8GU3zlP/PgZsz81zg\nfGA3c5K/vl71s8AFmfkCqqXNH2V2899A9fPZtGzWiDiP6jrcefVj3hIRpZ/Ctlz+9wLfnJkvBD4O\nXAlzlZ+I2Ax8P/Dpxm0nnH/aJzfMC6VmSmben5l31eMvAbuonuu/+IKu+r8/XCbhYBFxBnAx8Dc8\n/hTWuchfz7a+OzPfDtV1oMz8InOSH3iYapJwSkSsB06heoLCTObPzP8AHuy7eaWslwDvzMwjmbkP\n2Ev1M17Mcvkz87bMPFZv7qR6XQ7MSf7anwBv6LvthPNPu+CXexHUpil/z4mpZ2MvovpLsrHx7KDD\nwJRerTsRfwr8BnCscdu85D8L+J+IuCEiPhwRb4uIpzIn+TPzC1TvX/0ZqmJ/KDNvY07y11bK+hyq\nn+Hj5uHn+fXAzfV4LvJHxCXAgcz8SN9dJ5x/2gU/t1dwI+JrgX8AfjkzH2nel9WV6Zk8t4j4IeCB\nzLyTJ74ADZjt/FTP7LoAeEtmXkD1zKwlyxmznD8izqZ6z+ozqX4gvzYifqK5zyzn7zdE1pk9j4j4\nLeCxzPy7VXabqfwRcQrV+65c1bx5lYesmn/aBX8Q2NzY3szS30AzKSKeRFXu78jMf6pvPhwRz6rv\nfzbwQKl8A7wUeHVEfAp4J/B9EfEO5if/AarZywfr7XdRFf79c5L/24H/zszPZ+ZR4CbgJcxPflj5\n70r/z/MZ9W0zJyJ+mmqZ8scbN89D/rOpJgd31z/DZwAfiur9vU44/7QLfvGFUhFxEtUFgh1T/p5j\niYgArgfuzcw/a9y1A/ipevxTwD/1P3YWZOYbM3NzZp5FdXHv3zLztcxP/vuB/RHx/PqmlwMfA97N\nHOSnuiC8NSKeUv9dejnVxe55yQ8r/13ZAfxoRJwUEWcBW6g+L3OmRMQ2qiXKSzLzK427Zj5/Zt6T\nmRsz86z6Z/gA1QX7w4ySPzOn+gW8EriP6oLAldP+fhPI+11Ua9d3AXfWX9uA04B/oboq/15gQ+ms\nQ5zLy4Ad9Xhu8lN9assHqT6u6ybg6XOW/w1Uv5TuobpI+aRZzU/1r7xDVO/BvZ/qhYorZqVaPthL\n9YvsFTOY//XAHqpnnxz/+X3LHOT/v+N//n33fxI4bdT8vtBJklqq9HNAJUlTYsFLUktZ8JLUUha8\nJLWUBS9JLWXBS1JLWfCS1FIWvCS11P8DaqML5VItc2oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff42135ded0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Next, let's fit the normal distribution to our daily precipitation data - first the probability density function.\n",
    "P_daily = precip_daily.as_matrix()\n",
    "#plt.hist(P_daily[:,3],   bins = 20, normed = True)                      # plot the daily precipitaion data\n",
    "print(P_daily[:,3].max())                                               # What's the largest number we plot?\n",
    "\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "#plt.plot(x,stats.norm.pdf(x,my_mean,my_std), color=\"black\", linewidth=3)# plot the normal distribution curve\n",
    "#plt.show()\n",
    "# Does the curve fit the data?\n",
    "\n",
    "# Second, we fit the cumulative distribution function to the data:\n",
    "\n",
    "plt.hist(P_daily[:,3],   bins = 20, normed = True, cumulative=True)     # Add cumulative=True               \n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "#plt.plot(x,stats.norm.cdf(x,my_mean,my_std), color=\"black\", linewidth=3)# changes stats.norm.pdf into stats.norm.cdf\n",
    "#plt.show()\n",
    "# Does the curve fit the data? Can we conclude that the daily precipitation is normally distributed?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard normal distribution\n",
    "\n",
    "Because $f(x)$ is not analytically integrable, we use a table of probabilities to find the integral up to a value x. This table is developed for the case when the normal distribution has $\\mu=0$ and $\\sigma=1$, that is <br/ >\n",
    "$f(z) = \\frac{1}{\\sqrt{2 \\pi}}\\exp(\\frac{-z^2}{2})$ <br/ >\n",
    "\n",
    "This is called the **standard normal distribution**, and we replace x with _z_ whenever we use it.\n",
    "\n",
    "If you're data is normally distributed with different $\\mu$ and $\\sigma$, transform it into the standard normal by subtracting $\\mu$ from all data, and dividing by $\\sigma$. Then look up in the table: https://en.wikipedia.org/wiki/Standard_normal_table#Cumulative\n",
    "Last, transform you data back to what they were.\n",
    "\n",
    "#### Example 4.3 in the compendium\n",
    "\n",
    "What is the probability for any given year the annual maximum discharge is less than 8 m$^3$/s <br/ > for annual maximum  discharge having mean = 10m$^3$/s and standard deviation = 3 m$^3$/s?\n",
    "\n",
    "First, standardize: <br/ >\n",
    "\n",
    "$$\\frac{8-\\mu}{\\sigma}$$\n",
    "\n",
    "insert values\n",
    "\n",
    "$$\\frac{8-10}{3} = -2/3 = -0.67 $$\n",
    "\n",
    "The Z value corresponding to 8 m$^3$/s is -0.67.\n",
    "\n",
    "When we look up in the table, we see that Z = -0.67 corresponds to a probability of 0.251.\n",
    "Thus, the probability that the annual maximum discharge is less than 8 m$^3$/s is 25 %.\n",
    "\n",
    "Let's expand this example and make a table of different discharge values and their corresponding Z values and probabilities:\n",
    "\n",
    "What is the probability that the discharge is less than 1, 4, 7, 10, 13, 16 and 19 m$^3$/s?\n",
    "\n",
    "|x   | calculation|   = Z  | p(z) |\n",
    "|------|----------|--------|------|\n",
    "|x<1   | (1-1)/3  | = -3   |  0.1 |\n",
    "|x<4   | (1-4)/3  | = -2   |  2.3 |\n",
    "|x<7   | (1-7)/3  | = .... |  ... |\n",
    "|x<10  | (1-10)/3 | = .... |  ... |\n",
    "|x<13  | (1-13)/3 | = .... |  ... |\n",
    "|x<16  | (1-16)/3 | = .... |  ... |\n",
    "|x<19  | (1-19)/3 | = .... |  ... |\n",
    "\n",
    "What is the probability that the discharge is between 7 and 13 m$^3/s$, 4 and 16 m$^3/s$, 1 and 19 m$^3/s$?\n",
    "\n",
    "$P(7 < x < 13) = $ ?<br/ >\n",
    "$P(4 < x < 16) = $ ?<br/ >\n",
    "$P(1 < x < 19) = $ ?<br/ >\n",
    "\n",
    "Learning goal 4: be able to calculate probabilities for normally distributed variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The 68-95-99 rule, or the empirical rule\n",
    "For normal distributions, <br/ >\n",
    "68 % of the observations fall within 1 standard deviation around the mean <br/ >\n",
    "95 % of the observations fall within 2 standard deviations around the mean <br/ >\n",
    "99 % of the observations fall within 3 standard deviations around the mean (almost all observations fall within 3 $\\sigma$).\n",
    "\n",
    "This property can be used to justify the use of the normal distribution even if the hydrological variable is bounded by x = 0. See example 4.1 in the compendium. $\\mu$ = 600, $\\sigma$ = 100.\n",
    "\n",
    "Z = $\\frac{0 - 600}{100} = -6$\n",
    "\n",
    "What does this result tell us?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central limit theorem \n",
    "Are the statistical properties of daily and monthly observations the same? What about annual observations? \n",
    "\n",
    "Monthly averages (normalized sum) of daily data - is found by taking the average of approx. 30 daily data. There are 365 daily data in an annual average. What happens to the variability when we take the average of many observations?\n",
    "\n",
    "$\\bar{X} = \\frac{1}{n}(X_1 + X_2 + ... + X_n)$  is approximately ~N($\\mu$, $\\frac{\\sigma}{\\sqrt(n)}$) <br/ >\n",
    "\n",
    "The central limit theorem says that the average of a variable is normally distributed when $n$ gets large. This has the practical implication that if you can divide a variable into independent subsets that are equally big (think: measuring the precipitation once a day), the total (annual rainfall) is approximately normally distributed. This is why the normal distribution shows up so many places.\n",
    "\n",
    "\"[T]he central limit theorem allows us to use normal probability calculations to answer questions about sample means (read: monthly or annual averages) even when the population distribution (read: daily data) is not normal\" (Moore and McCabe, Introduction to the practice of statistics)\n",
    "\n",
    "In the compendium, it says that the sum (instead of average) of a variable is normally distributed when $n$ gets large: <br/ >\n",
    "$\\bar{X} = (X_1 + X_2 + ... + X_n)$  is approximately ~N(?, ?)\n",
    "\n",
    "######What is the expected mean and average of the sum of variables? <br/ >\n",
    "Hint: use the \"Reproductive properties\" on page 50 in the compendium.\n",
    "\n",
    "Learning goal 5: understand what the central limit theorem means."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normal approximations to other distributions\n",
    "\n",
    "As a result of the central limit theorem, many distributions approach the normal distribution under certain conditions. The hypergeometric, binomial and poisson distributions approach normal when $n$ becomes large. Here, we may use a _continuous_ distribution to find approximate probabilities for _discrete_ distributions. This is very useful, because it is relatively simple to find probabilities for the standard normal distribution.\n",
    "\n",
    "Learning goal 6: understand why the normal distribution is used to approximate probabilities for discrete distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial distribution\n",
    "A sum $Z = X_1+X_2+...+X_n$ of binomial random variables, approaches the normal distribution when $n$ is large.\n",
    "\n",
    "In the binomial distribution, we have two parameters: n and p<br/ >\n",
    "$\\mu = np$ <br/ >\n",
    "$\\sigma = np(1-p)$ \n",
    "\n",
    "In example 4.6 we're asked to find the probability that X is between 5 and 8. We'll approximate the binomial distribution with the normal distribution.\n",
    "\n",
    "The binomial random variables approach the standard normal distribution when we standardize:\n",
    "\n",
    "Z = $\\frac{X-\\mu}{\\sigma} = \\frac{(X - np)}{\\sqrt{(np(1-p))}}$\n",
    "\n",
    "Note that the approximation becomes even better when we adjust the X values according to table 4.2 in the compendium.\n",
    "It means that we include whole rectangles of the histogram. See figure 2 here:\n",
    "http://revisionmaths.com/advanced-level-maths-revision/statistics/normal-approximations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Example 4.6\n",
    "my_n = 25\n",
    "my_p = 0.3\n",
    "my_mean = my_n*my_p\n",
    "my_std = np.sqrt(my_n*my_p*(1-my_p))\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "#plt.plot(x,stats.binom.cdf(x,my_n,my_p, loc=my_mean), color=\"red\", linewidth=3) # plot the binomial distribution curve\n",
    "#plt.plot(x,stats.norm.cdf(x,my_mean,my_std), color=\"black\", linewidth=1)# plot the normal distribution curve\n",
    "#plt.show()\n",
    "\n",
    "\n",
    "#print('The true value of P(5<X<8) = ', stats.binom.cdf(8,my_n,my_p)-stats.binom.cdf(5,my_n,my_p))\n",
    "#print('Approximate value of P(5<X<8) with standard normal = ', 0.587-0.138)\n",
    "#print('Approximate value of P(5.5<X<8.5) with standard normal = ', 0.670-0.192)\n",
    "\n",
    "## The approximation is most accurate when p is close to 0.5. Try changing my_p above into 0.5 and 0.99."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson distribution\n",
    "In the Poisson distribution, we have one parameter: $\\lambda$<br/ >\n",
    "$\\mu = \\lambda$ <br/ >\n",
    "$\\sigma =\\sqrt{\\lambda}$ \n",
    "\n",
    "To transform to the standard normal distribution, we standardize:\n",
    "\n",
    "Z = $\\frac{X-\\mu}{\\sigma} = \\frac{(X - \\lambda)}{\\sqrt{\\lambda}}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_lambda=2\n",
    "my_mean = my_lambda\n",
    "my_std = np.sqrt(my_lambda)\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "#plt.plot(x,stats.poisson.cdf(x,my_lambda), color=\"green\", linewidth=3) # plot the binomial distribution curve\n",
    "#plt.plot(x,stats.norm.cdf(x,my_mean,my_std), color=\"black\", linewidth=1)# plot the normal distribution curve\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lognormal distribution\n",
    "The next distribution we'll go through is not a discrete distribution like the binomial and poisson.\n",
    "The lognormal distribution is used when the observations are not normally distributed, but when the _logarithm_ of the observations are normally distributed.\n",
    "\n",
    "The lognormal distribution is known for its heavy left tail. Because the logarithm of negative numbers is between and 0 and 1, the lognormal distribution is bounded by zero.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "## First, we generate data that are not normal, by applying a mathematical function to normal data\n",
    "my_mean=100\n",
    "generated_norm = np.random.normal(loc=my_mean, scale=my_std, size=500)\n",
    "data = np.exp(generated_norm)\n",
    "\n",
    "## Take the logarithm of the data:\n",
    "log_data = np.log(data)\n",
    "\n",
    "## Plot a new histogram of the logarithm of the dat\n",
    "## plt.hist(data,color=\"red\", bins = 12, normed= True, alpha=0.7)      # plot the non-normal data\n",
    "#plt.hist(log_data,   bins = 12, normed = True)                   # plot the logarithm of those data. They are normal.\n",
    "\n",
    "## Next, we fit the cumulative distribution function to the data:\n",
    "x = np.linspace(my_mean-3*my_std,my_mean+3*my_std,100)                  # define the x values that we want in our plot\n",
    "#plt.plot(x,stats.norm.pdf(x,my_mean,my_std), color=\"black\", linewidth=3)# changes stats.norm.pdf into stats.norm.cdf\n",
    "#plt.show()\n",
    "\n",
    "# If log_data is normally distributed with my_mean and my_std, we say that data is lognormally distributed with my_mean and my_std."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chi square distribution ($\\chi^2$)\n",
    "Up until now, we've looked at the distribution of sums or averages ($\\bar{X}$) of random variables. This next distribution is used to find the distribution of the sum of squares ($X^2$).\n",
    "It is among other things used in calculating confidence interval of variance $\\Big( S^2 = \\frac{1}{n-1}\\sum_{i=1}^n(X-\\bar{X})^2 \\Big)$.\n",
    "\n",
    "If $X_1, X_2, ..., X_n$ are standard normal random variables, the sum\n",
    "\n",
    "$Y = X_1^2 + X_2^2 + ... + X_ n^2$\n",
    "\n",
    "is $\\chi^2$ distributed with n degrees of freedom.\n",
    "\n",
    "The Chi square table is given here: https://en.wikipedia.org/wiki/Chi-squared_distribution#Table_of_.CF.872_value_vs_p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Student's t distribution (the t-distribution)\n",
    "\n",
    "William Gosset published this distribution under the pseudonym \"Student\" in 1908. The statistical methods at the time assumed large sample sizes, which did not approptiately describe the result he obtained from small sample sizes. (Anecdote taken from Løvås, Statistikk for universiteter og høyskoler, p. 165).\n",
    "\n",
    "The sample mean is $\\bar{X} = (X_1 + X_2 + ... X_n)/n$ (approximately normally distributed).\n",
    "The sample variance is  $S^2 = \\frac{1}{n-1}\\sum_{i=1}^n(X-\\bar{X})^2$\n",
    "\n",
    "Then, $Z = \\frac{(\\bar{X}-\\mu)}{\\sigma/\\sqrt{n}}$ follows the standard normal distribuion.\n",
    "\n",
    "In practical situations, we do not know the true standard deviation of the population, $\\sigma$, only $S$.\n",
    "\n",
    "Gosset found the proabaility distribution of   $T = \\frac{(\\bar{X}-\\mu)}{S/\\sqrt{n}}$, which approaches $Z$ for large sample sizes. \n",
    "The T distribution has heavier tails than the standard normal distribution.\n",
    "\n",
    "For small degrees of freedom ($n$), $P(T \\leq 2)$ is much larger tgab $P(Z \\leq 2)$, which means that it is more probable to observe large values of the $T$ variable than the $Z$ variable. The spread of $T$ is larger than the spread of $Z$.\n",
    "\n",
    "The T table is given here: https://en.wikipedia.org/wiki/Student%27s_t-distribution#Table_of_selected_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fischer F distribution\n",
    "Is used to compare the variance of two populations, simply by dividing one by the other.\n",
    "\n",
    "$X = \\frac{U/m}{V/n}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extreme value distributions\n",
    "\n",
    "A standard hydrology task is to calculate the 100-year flood, or 1000-year flood, or the probability of low-flow. In that case, we are interested in the extreme values instead of daily values. We construct a new time series by selecting only the extreme values. This new time series is _very_ skewed (has a long tail), and the appropriate curves to fit tho this data are extreme value distributions.\n",
    "\n",
    "Three types of extreme value distributions exist, each having a heavy tail:\n",
    "\n",
    "|       Gumbel                 |         Fréchet                 | Weibull                          |\n",
    "|------------------------------|---------------------------------|----------------------------------|\n",
    "|              (Type I)        |                 (Type II)       |         (Type III)               |\n",
    "| Heavy tail to left and right |Heavy tail to the right (floods) | Heavy tail to the left (low-flow)|\n",
    "| Annual maximum river flow    |                                 | Low-flow                         |\n",
    "\n",
    "Gumbel:\n",
    "$f(x)= \\frac{1}{\\sigma}e^{-z-e^{-z}}$\n",
    "\n",
    "where $z = \\frac{x-\\mu}{\\sigma}$\n",
    "\n",
    "Fréchet:\n",
    "$f(x)= \\frac{\\alpha}{\\beta}\\big(\\frac{\\beta}{x}\\big)^{\\alpha+1}e^{(\\beta/x)^\\alpha}$,  <br/ >\n",
    "where $\\alpha, \\beta > 0$\n",
    "\n",
    "Weibull:\n",
    "$f(x)= \\frac{\\alpha}{\\beta}\\big(\\frac{x}{\\beta}\\big)^{\\alpha-1}e^{(x/\\beta)^\\alpha}$, <br/ >\n",
    "where $\\alpha, \\beta > 0$\n",
    "\n",
    "(note! If these formulas differ from the compendium, I've made typos)\n",
    "\n",
    "\n",
    "\n",
    "Learning goal 7: know that lognormal, Chi square, Student's t, F distribution and  three extreme value distributions exist: Gumbel, Fréchet and Weibull."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After this lecture, you should be able to ...\n",
    "\n",
    "Learning goal 1: understand what the central, spread and skewness measures are and how they are defined mathematically. (Ex1, Q3) <br/>\n",
    "Learning goal 2: understand what happens if the mean, standard deviation or skewness is changed. <br/>\n",
    "Learning goal 3: understand the difference between the pdf and the cdf. <br/>\n",
    "Learning goal 4: be able to calculate probabilities for normally distributed variables. (Ex1, Q2)<br/>\n",
    "Learning goal 5: understand what the central limit theorem means. <br/ >\n",
    "Learning goal 6: understand why the normal distribution is used to approximate probabilities for discrete distributions. <br/ >\n",
    "Learning goal 7: know that the lognormal, Chi square, Student's t, F distribution and three extreme value distributions exist: Gumbel, Fréchet and Weibull."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Coda\n",
    "__Course__: GEO4310 Statistical and Stochastic Hydrology, University of Oslo\n",
    "\n",
    "__Contributions__\n",
    "\n",
    "* Irene Brox Nilsen, 4. September 2015\n",
    "* John F. Burkhart, 6. September 2017\n",
    "\n",
    "__License__\n",
    "\n",
    "[Creative Commons Attribution-ShareAlike 4.0 International Public License](https://creativecommons.org/licenses/by-sa/4.0/legalcode)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:shyft]",
   "language": "python",
   "name": "conda-env-shyft-py"
  },
  "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}