GEO4310-Lecture02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 2: Properties of Random Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Outline for the day\n",
    "\n",
    "* Review Lecture 1\n",
    "* Properties of Random variables\n",
    "    - Probability distribution of random variables\n",
    "    - Main descriptor or random variables\n",
    "* Wrap-up review\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review of Lecture 1\n",
    "\n",
    "### Mathematics of Probability, Basic Axioms and the Additive Rule\n",
    "\n",
    "The *Theory of Probability* is based on fundamental assumptions, or **axioms**, as follows:\n",
    "\n",
    "For every event $E$ in the sample space, $S$, there exists a probabilty:\n",
    "\n",
    "**[1]** $ P(E) \\geq 0$\n",
    "\n",
    "Secondly, the probability of the *certain event*, $S$ is:\n",
    "\n",
    "**[2]** $ P(S) = 1.0 $\n",
    "\n",
    "Finally, for two events $E_1$ and $E_2$ that are *mutually exclusive*:\n",
    "\n",
    "**[3]** $ P(E_1 \\cup E_2) = P(E_1) + P(E_2) $\n",
    "\n",
    "\n",
    "These axioms (assumptions) are fundamental to and constitute the basis for *probability theory*. They are not subject to proof.\n",
    "\n",
    "What about $E$ with it's complement $\\bar{E}$? \n",
    "\n",
    "$ P(E \\cup \\bar{E}) = P(E) + P(\\bar{E}) $\n",
    "\n",
    "Recall that the **union** of $E$ with it's complement, $\\bar{E}$ is the sample space $S$. That is,\n",
    "\n",
    "$ P(E \\cup \\bar{E}) = P(S) = 1.0 $\n",
    "\n",
    "Therfore ($\\therefore$), \n",
    "\n",
    "$P(\\bar{E}) = 1 - P(E) $\n",
    "\n",
    "#### The general case for two events\n",
    "\n",
    "We can think of the above relationship more generally for two events $E_1$ and $E_2$ in the case that they are not mutually exclusive:\n",
    "\n",
    "$P(E_1 \\cup E_2) = P(E_1) + P(E_2) - P(E_1 \\cap E_2) $\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional Probability and the Multiplication Rule\n",
    "\n",
    "If the probability of an event depends on the occurence of another even then we have **conditional probability**. \n",
    "\n",
    "What is the probability of $E_1$ given $E_2$? This means we assume $E_2$ has occured, and we want to find the probability of a sample point in $E_1$ that belongs to $E_2$:\n",
    "\n",
    "$P(E_1 | E_2) = \\frac {P(E_1 E_2)} {P(E_2)} $\n",
    "\n",
    "In essence, what we are doing is *normalizing* to a reconstituted sample space where the event has occured. That is, we want to know the likelihood of realizing a sample point in $E_1$ assuming that it belongs to $E_2$. \n",
    "\n",
    "We can rearrange this equation to get the intersection, assuming that $P(E_2) \\neq 0.$\n",
    "\n",
    "$P(E_1 \\cap E_2) = P(E_2) * P(E_1 | E_2) $\n",
    "\n",
    "However, if the two events are *statistically indpendent*, that is to say $P(E_2 | E_1) = P(E_2)$ then we can calculate the **joint probability of two statistically independent events** as their product:\n",
    "\n",
    "$P(E_1 \\cap E_2) = P(E_1) * P(E_2) $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Xu Example 1.4\n",
    "\n",
    "River **a** and river **b** are two adjacent tributaries of a bigger river. River **a** has a yearly discharge smaller than it's mean annual discharge $u_a$ in 50% of the years, while river **b** has a yearly discharge smaller than its mean $u_b$ in 60% of years.At any year the probability of of river **b** having discharge smaller than $u_b$, given that river **a** has a discharge smaller than $u_a$ is 70%.\n",
    "\n",
    "__Event A__: $Q_a < u_a; P(A) = 0.5 $\n",
    "\n",
    "__Event B__: $Q_b < u_b; P(B) = 0.6 $\n",
    "\n",
    "We also know,\n",
    "\n",
    "$P(B | A) = 0.7 $\n",
    "\n",
    "Calculate the probability:\n",
    "\n",
    "[1] That both rivers have a discharge smaller than their mean value:\n",
    "\n",
    "> Since we know that $P(A) \\geq 0$, we can arrange our conditional probability to determine the intersection:\n",
    "\n",
    "$ P(A \\cap B) = P(A) * P(B | A)$\n",
    "\n",
    "> We also know that $P(B|A) \\neq P(B)$, therefore we know they are not **statistically independent**, so that:\n",
    "\n",
    "$ P(A \\cap B) = P(A) * P(B | A) = 0.5 * 0.7 = 0.35 $\n",
    "\n",
    "[2] That *at least* one river has a discharge smaller than it's mean:\n",
    "\n",
    "$P(A \\cup B) = P(A) + P(B) ... $\n",
    "\n",
    "> But we know these are not mutually exclusive and statistically independent, so.\n",
    "\n",
    "$P(A \\cup B) = P(A) + P(B) - P(A \\cap B) = 0.5 + 0.6 - 0.35 = 0.75 $\n",
    "\n",
    "[3] That river **a** has a discharge smaller than it's mean, given river **b** had a discharge smaller than the mean:\n",
    "\n",
    "> We can use the relationship we already determined, nameley\n",
    "\n",
    "$ P(A | B) = P(A \\cap B) / P(B) = 0.35 / 0.6 = 0.5833$\n",
    "\n",
    "[4] That *at least* one river has a discharge higher than it's mean:\n",
    "\n",
    "$P (\\bar{A} \\cup \\bar{B}) = ... $\n",
    "\n",
    "> Here we'll use the full sample space, and that fact that we already know the probability that *at least* one river has a discharge smaller than it's mean [2]. \n",
    "\n",
    "$P (\\bar{A} \\cup \\bar{B}) = P(S) - P(A \\cap B) = 1.0 - 0.35 = 0.65 $\n",
    "\n",
    "[5] That both rivers have a discharge higher than their mean:\n",
    "\n",
    "$P(\\bar{A} \\cap \\bar{B}) = 1.0 - P(A \\cup B) = 1.0 - 0.75 = 0.25 $\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------\n",
    "\n",
    "###Xu Example 1.5 (from Haan Example 2.2)\n",
    "\n",
    "A study of daily rainfall at Station AK in July based on many years of observation has shown that:\n",
    "\n",
    "Probability of a rainy-rainy day (rainy day following a rainy day) is 0.444, \n",
    "\n",
    "a dry-dry day is 0.724, \n",
    "\n",
    "a dry-rainy day is 0.276, \n",
    "\n",
    "and a rainy-dry day is 0.556\n",
    "\n",
    "If it is observed that a certain July day (d1) is rainy, what is the probability that the next **two** days (d2 and d3) will also be rainy?\n",
    "\n",
    "#####Solution\n",
    "Let A be a rainy d2 and B be a rainy d3 following the initial rainy d1. The probability of A is 0.444 since this is the probability of a rainy-rainy day.\n",
    "\n",
    "The probability of two rainy days following a rainy day:\n",
    "\n",
    "$P(A \\cap B) = P(A) * P(B|A) $\n",
    "\n",
    "> We already know the probability of $P(B|A) = 0.444$ since this is the probability of a rainy day, following a rainy day. So to find the probability of the third consecutive rainy day, it is simply:\n",
    "\n",
    "$P(A \\cap B) = P(A) * P(B) = 0.444 * 0.444 = 0.197$\n",
    "\n",
    "The probability of two rainy days following a dry day would be $P(dry-rainy) * P(rainy-rainy)$:\n",
    "\n",
    "$P(dry-rainy) * P(rainy-rainy) = 0.276 * 0.444 = 0.122 $\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Properties of Random Variables\n",
    "\n",
    "Properties of Random Variables\n",
    "------------------------------\n",
    "\n",
    "A random variable is obtained by assigning a numerical value to each\n",
    "**outcome** of an experiment. Random variables may be **discrete** or\n",
    "**continuous**. In some cases, the outcomes of an experiment/observation\n",
    "are already numeric, and these values will form the randome variable of\n",
    "interest (e.g. streamflow, annual discharge, peak discharge, etc.). In\n",
    "other cases, we may have *categorical* data or other non-numeric data\n",
    "for which we need to assign values to create our set of random\n",
    "variables. An example may be ephemeral flow in a stream (e.g.\n",
    "flow/no-flow), or cloud cover (here we could use octals).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vocabulary\n",
    "\n",
    "__Outcome__ - The result of an experiment\n",
    "\n",
    "__Statistical Inference__ - use information from a sample to draw conclusions (inferences) about\n",
    "the population from which the sample was taken.\n",
    "\n",
    "__Population__ - a complete assemblage of all values representative of a random process\n",
    "\n",
    "__Sample__ - a group of values selected for studying the property of population\n",
    "\n",
    "__Parameter__ - a parameter is a value, usually unknown (and which therefore has to be estimated),\n",
    "used to represent a certain population characteristic.\n",
    "\n",
    "__Statistic__ - A statistic is a quantity that is calculated from a sample of data use to estimate the\n",
    "parameters of the population\n",
    "\n",
    "__Estimator__ - An estimator is any quantity calculated from the sample data which is used to give\n",
    "information about an unknown quantity in the population. For example, the sample mean is an\n",
    "estimator of the population mean.\n",
    "\n",
    "__Estimate__ - an estimate is the particular value of an estimator that is obtained from a particular\n",
    "sample of data and used to indicate the value of a parameter. If the value of the estimator in a\n",
    "particular sample is found to be 5, then 5 is the estimate of the population mean μ.\n",
    "\n",
    "__Estimation__ - Estimation is the process by which sample data are used to indicate the value of an\n",
    "unknown quantity in a population\n",
    "\n",
    "__Moments__ - tendency or parameters that describe a distribution\n",
    "\n",
    "\n",
    "## Moments of a Distribution\n",
    "\n",
    "When we consider a statistical analysis of data, and start to think about predictability of the events the data represents (say certain discharge quantities for example), we need to have some mechanism to describe that probability of the event occurring.\n",
    "\n",
    "This is where distributions come in. Say we have a bag that holds a collection of objects (or data) that creates a sample. Let's say that collection of objects are temperature observations in San Francisco, CA over 100 years. If you are asked \"What is the temperature in San Francisco tomorrow?\", how would you proceed?\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pair of rolled dice - example from Wikipedia\n",
    "\n",
    "A random variable can also be used to describe the process of rolling dice and the possible outcomes. The most obvious representation for the two-dice case is to take the set of pairs of numbers n1 and n2 from {1, 2, 3, 4, 5, 6} representing the numbers on the two dice as the sample space, defining the random variable X to be equal to the total number rolled, the sum of the numbers in each pair. In this case, the random variable of interest X is defined as the function that maps the pair to the sum:\n",
    "\n",
    "$ X((n_1, n_2)) = n_1 + n_2 $\n",
    "\n",
    "and has probability mass function $p_X$ given by:\n",
    "\n",
    "$ p_X(S) =  \\tfrac{\\min(S-1, 13-S)}{36}, \\text{for } S \\in \\{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12\\} $\n",
    "\n",
    "<img src=\"http://upload.wikimedia.org/wikipedia/commons/thumb/1/12/Dice_Distribution_%28bar%29.svg/200px-Dice_Distribution_%28bar%29.svg.png\">\n",
    "\n",
    "### Probability Mass Function\n",
    "\n",
    "The probability mass function, or $p.m.f.$ is for **discrete** random\n",
    "variables. We are interested in knowing the probability properties of\n",
    "each outcome in our sample space. Thus, we can assign a value of $p_i$\n",
    "to each of the $x_i$ in the random variable. Note, that each $p_i$ is\n",
    "thus between $0$ and $1$, and the $\\sum_i{p_i} = 1$\n",
    "\n",
    "Repeating:\n",
    "\n",
    "$p.m.f.$ of $X$ is a set of probability values ($p_i$) for each $x_i$ in\n",
    "a **discrete** random variable, where:\n",
    "\n",
    "$$0 \\le p_i \\le 1$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\sum_i{p_i} = 1$$\n",
    "\n",
    "We are interested in the *probability* that a value drawn from the\n",
    "random variable $X$ *takes* the value $x_i$, which is written:\n",
    "\n",
    "$$P(X = x_i) = p_i$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwEAAAF8CAYAAACXPLcjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X24XWV95//3pwk+RErR4ScWiEIFFawP0BYZO46J4hgZ\nBX92ZjDVWnEu5fpNg1Qzig8zndOrM512Wlu1jP5QEXS0TVu0/WGhoq0e+6DyUAmiBCUqmkBBRNFi\n1Cby/f2xduLO5pxk7/Ow19l7vV/Xda7se617rfW9zzrnZH3Xfa97paqQJEmS1B0/1nYAkiRJksbL\nJECSJEnqGJMASZIkqWNMAiRJkqSOMQmQJEmSOsYkQJIkSeoYkwBJkiSpY0wCJEmSpI4xCZDUCUku\nTfIbYz7mY5NsTfKdJJvGeexxSnJrkmf0fX7mAvcz9DkarJvkc0n+9UKOO8kW8/2W1G2r2w5Aksak\nel/j9Frgr6vqyUuxsyS3Ai+rqo8txf6WUA18Xuj3eZRt96tbVT+9wGNOujZ+riVNAZMASa1L8irg\nCOBrVXXRch5qwRsm64H3A0+pqh19yw8U+6OATy7gWKuras8cq4pFtGEhDhDLsh1ymepKkvo4HEhS\nq5L8BPAfgP8P+LuBda9Lsr03nObzSZ4/sP7WJJuT3JDkniRbkjywt+7kJJ/pbbsFeNAiQ/1b4M6B\nBOBAsX8MWAdc2Ivh+CQnJplN8q3e8JXnDbTltUk+C/xTkh8b2N//AR4JfCjJPyV5TZKXJrm8r84t\nSf6kr7wjyRN7n+c99qC5Yhll+4F9HZXkA0m+nuTLSc7rWzf0OTpY3f5hMUnWJvlg75jfSPIHw8Qz\nYuyPTnJ3kpP76t61d0jSgX52e7H+5ySf7Z3Li5McmeQvk3w7yUeTHD5Q/3W9/Xwzybv3/pyPErMk\n9TMJkNS2pwBbq+qaqvr8wLrtwL+qqsOAXwfel+QRfesL+PfAs4HjgCcCL03yAODPgfcADwX+FPgF\nFjds4meAa4eNvaqeQZM4/Eov/q8CHwI+DPxfwHnA+5Oc0LfZC4HnAIdX1X0D+/sl4GvAc6vqx6vq\nd4BPAE+D5uIPOAQ4rVf+KeAhVfXZJIfMc+zHHKC9+2IBVg0R+/30EpkPAdcDRwHPBH41yb8Z5RwN\nWbeASrIK+AvgKzQ9MUcDWw4WzyixA1TVl4ALaH4mHwxcAlxSVX/T28VcP7tH9sX6gt4+Hws8F/hL\n4HXAw2n+b37lQEi/CPwb4NHAY4D/MmrMktTPJEDSSJI8KcnLkvx2krOSvDzJS5I8vrf895K8IMkv\nJHlfktOTPCfJ/0ly/MC+ngKcD6xO8n8PHquqLquqO3qf/wS4BTh1oNpbq+qOqvoWzQXQk2kuhFdX\n1Vuq6odV9QHufwE/2K6jk/xaL9brkjwgyWOS/EaSDcDrgb8eNvY5nEZzUf5bVbWnqj5Oc7H6i3ub\n22vLbVX1gyH2R1V9heZO/cnAvwauAm5P8ljg6cDeC9L5jr1xvl0PxHKw2Ofzc8ARVfXfe9t9BXgX\nTYIxyjkape6pwE8Cr6mq71XVD6rq74eIZ5TYm29S1btoLvavAY4E3ti37mA/u39QVXdV1e00yeKn\nquqG3vf7z4CT++oWcGHvfHwL+B/Mfe5GaZ+kjvOZAEmjOhL4AvDsqrogyUNo7jyeB2wFXlxVrwZI\n8lvAr1XVl5OcCTyO5qIJgKq6Osn3gDfP0QtAkpcArwKO7S06FPgXA9Xu6Pu8i+YO6E8Ctw3U+yrz\njCHvteHPgOdU1d1J/obmrvqfAOuq6p4krwX2PZB7sNjncBSwY2DZV3vL9xpcP4xP0Aw7Or73+R6a\nBOBf9soHOvbRB9hvf/1hYp/Lo4Cjknyrb9kqmuRklHN01Ah11wJfHexJGSKehdZ9F81wsJdX1e69\nC+f52T2ib7s7+z5/b6D8/V79fv3f/68x9/d+lPZJ6jh7AiSNpKo+QjMs4UO9RScD36iqq4BnAe+D\nZsw08OWq+nKv3tOBT82xyxOBmwYXJnkU8A7gV4CHVdVDgc8x3MOg/8j9L3AfxfzDgc4Grququ3tt\n/C7NcI0bewnAg4BDq+quYWKfx+3A2iT98T+K/S9uDzZcaa71nwDW0wwLmuVHScHT+VESMN+xdw55\nrGFin8sO4CtV9dC+r8Oq6rk0yduw52iU8/k14JG9YUFzrZsvnpHrJjkUeDNNIvDrSR7aW76Qn92D\n/Vw/cuDzXN/7A32/JWk/JgGSFuJ0fnSB+cvA7/Yt/0jv87NohqeQ5Kk0wyGOTHLs3p30xkh/o6rm\nuph7CM1F3jeAH0tyDnCwaSD3Xkh9CtiT5JVJDknyApqhEvNZTV8PRZLTaO7g3tDXlk/3hgUNE/tc\nMV1N01Px2l5M62jGgm85yPb97qQZE95vbxLwoN7Qkr8DNgAPo+mhAfj0Io+90O2vphmu9NokD06y\nKslPJ/lZmlmThj1Ho9S9hiZp+K0ka5I8qPfzt3fdfPHMtZ+D1X0LcE1VvQK4Avh/e8sX8rN7IAH+\nU2/I2sNohh398Rz1DvT9lqT9mARIGkmaGXEeBjwjycuBq6vqg72HEg+pqr13l0+hGTcO8EOau8mP\nrapb+3b3FODvmUNV3QS8ieaC/g6ai6i/m6tu/2bNprWb5k7+S4G7aWbw+cABtvsj4OFJnte7wPwx\n4J3AMUmeQzOMYzfw48PEPkdMVNU/A8+jedj2LuBC4Jeq6otD7GOv/wn8lzQz9Ly6t99bgH+iGVdO\nVX0H+BLw93sTlN73Y8HHXuj2vSE5z6V5TuPLvW3fARw2yjkase59vViPp7mbv6NX/4DxjBI7QJKz\naHrE/p/eJq8GTkmycRE/u/2fB8t/SJNgf4kmof7vo8YsSf1y8JtYy3DQ5m7am2nGKr6rqn57YP3j\naGZaOBl4Y1W9qW/d64EXA/cBNwLnDPsQnaTF6z0Ee1pVXbCIffwM8HLgm8AfV9UNB9lkxZjk2DWZ\nknwF+I+18l4SJ2kOSd4N/Fvg61X1hHnqvJXmxsou4KVVdf1c9ZbT2HsCeuM0L6Tprj4J2JjkxIFq\nd9M8ZPi7A9seS/Of7ym9b+oqnPVAGptegv5qmrvmi7m7+EOa8ejfmMCL6EmOXZK0/C6huc6dU5Iz\ngOOr6gTgFcDbxxVYvzZmBzoV2L53SECal76cBWzbW6H38N1dSf7twLbfoemSX5Pkh8AaDv5gmqQl\nUlU305uXfpH72Uozk9DEmeTYJUnLr6r+tv/5tzmcSfPek70zzR2e5MiquvMA2yy5Np4JOJr9pzrb\nyYGnqdunqr5JM87yazTji++pqr9a8gglSVohquo4hwJJU2Wua+Fjxh1EG0nAgh9C6E05+Ks0s3Yc\nBRya5EVLFJckSZI0DoPTAo/9Id02hgPdRvMyl73WcuC5qvv9LPDJvXN5J/kg8FTg/f2Vkoz/aWdJ\nkiR1UlUN8w6bvQavhY+hheHtbSQB1wEn9MZK3U7zkp75Xl0/+A29GfivSR5M80bF02nmcr6fEU+G\nVpAkM1U103YcWhjP3+Ty3E02z99kmp3tzo3Ldeum89psATefLwc2AVt676W5Z9zPA0ALSUBV7Umy\nieYlQquAi6tqW5Jze+svSvII4FqauY3vS3I+cFJV3ZDkvTSJxH3AZ2jmQJYkSZpo03qR3DVJ/ojm\nre1HJNkB/DfgEGiuc6vqyiRnJNkOfBc4p4042+gJoKr+EvjLgWUX9X2+g/27Sfrr/S/gfy1rgJIk\nSdICVNV8I1z662waRywH4huDtRLNth2AFmW27QC0YLNtB6BFmW07AEmTwyRAK05VzbYdgxbO8ze5\nPHeTzfMnaRQmAZIkSVLHmARIkiRJHdPKg8GSJEnDmp3NZmAGOLTlUKSpYU+AJEla6WaY/gTg3rYD\nULeYBEiSpJWuCwnATNtBqFscDiRJkiaGL9SSloY9AZIkSVLHmARIkiRJHWMSIEmSJHWMSYAkSZLU\nMSYBkiRJUseYBEiSJEkdYxIgSZIkdYxJgCRJktQxJgGSJElSx5gESJIkSR1jEiBJkiR1jEmAJEmS\n1DEmAZIkSVLHrG47AEmaFEmq7RgWoqrSdgySpJXFJECSRjJpeYDX/5Kk+zMJkCRpws3OZjMwAxza\nciiSJoTPBEiSNPlm6EYCcG/bAUjTwiRAkqTJ15UEYKbtIKRp4XAgSZKmyLp1Pggu6eDsCZAkSZI6\nppUkIMmGJDcnuSXJBXOsf1ySTyX5fpLNA+sOT3JZkm1Jbkpy2vgilyRJkibf2IcDJVkFXAicDtwG\nXJvk8qra1lftbuA84Plz7OItwJVV9e+SrAYestwxS5IkSdOkjZ6AU4HtVXVrVe0GtgBn9Veoqruq\n6jpgd//yJD8BPK2q3t2rt6eqvj2muCVJkqSp0EYScDSwo6+8s7dsGMcBdyW5JMlnkrwzyZolj1CS\nJEmaYm0kAYt53eZq4BTgbVV1CvBd4HVLEpUkSZLUEW1MEXobsLavvJamN2AYO4GdVXVtr3wZ8yQB\nSWb6irNVNTtamJIkSdJ0aiMJuA44IcmxwO3A2cDGeeruN9dxVd2RZEeSx1TVF2keLv78XBtW1cxS\nBSxJkiRNk7EnAVW1J8km4CpgFXBxVW1Lcm5v/UVJHgFcCxwG3JfkfOCkqrqXZtag9yd5APAl4Jxx\nt0GSJEmaZKlazBD9lSlJVfnGRElLK0kt7rGmNgT/Hk6/2dns+8H0jcHSeE3qdadvDJYkSZI6xiRA\nkiRJ6hiTAEmSJKljTAIkSZKkjjEJkCRJkjrGJECSJEnqGJMASZIkqWNMAiRJkqSOMQmQJEmSOsYk\nQJIkSeoYkwBJkiSpY0wCJEmSpI4xCZAkSZI6xiRAkiRJ6pjVbQcgSdJym53NZmAGOLTlUCRpRbAn\nQJLUBTN0IwG4t+0AJE0GkwBJUhd0JQGYaTsISZPB4UCSpE5Zt67SdgyS1DZ7AiRJkqSOMQmQJEmS\nOsYkQJIkSeoYkwBJkiRpCSXZkOTmJLckuWCO9Uck+XCSrUk+l+Sl447RJECSJElaIklWARcCG4CT\ngI1JThyotgm4vqqeDKwD3pRkrBP2mARIkiRJS+dUYHtV3VpVu4EtwFkDdf4ROKz3+TDg7qraM8YY\nnSJUkiRJWkJHAzv6yjuBpwzUeSfwsSS3Az8O/IcxxbaPPQGSJEnS0qkh6rwB2FpVRwFPBv53kh9f\n3rD2N7VJQJJ1/Z8tW7ZseSnKMNv7mpSyfw/7y1u3+v2wbNny4sq9r5kklya5lPu7DVjbV15L0xvQ\n76nAnwJU1ZeArwCPnWNfyyZVwyQrkyVJVflGSElLK0kNd4NnJQn+PYTZ2ew7cb4xWNJSGrzuTPOA\n7xeAZwK3A9cAG6tqW1+d3wO+XVW/nuRI4B+AJ1bVN8cVt88ESJIkSUukqvYk2QRcBawCLq6qbUnO\n7a2/CPhN4JIkN9CMzHntOBMAaKknIMkG4M0035h3VdVvD6x/HHAJcDLwxqp608D6VcB1wM6qet4c\n+7cnQNKSsydgctkTIGm5TOp159ifCchwc6feDZwH/O48uzkfuInJ+99YkiRJal0bDwYfdO7Uqrqr\nqq4Ddg9unOQY4AzgXcDEZV2SJElS29pIAuaaO/XoEbb/feA1wH1LGZQkSZLUFW08GLzgITxJngt8\nvaqu75+yaZ66M33F2aqaXehxJUmSpGnSRhIwzNyp83kqcGaSM4AHAYcleW9VvWSwYlXNLDZQSZIk\naRq1kQRcB5yQ5FiauVPPBjbOU3e/Mf9V9QaaN6yR5OnAf54rAZDUnmYGnckyibM6LIdJPHfg+ZOk\nhRh7EjDM3KlJHgFcCxwG3JfkfOCkqrp3cHfjjF3SsCbpV9Prx/1N0rkDz58kLYxvDJa0pCZvLv3h\n59GfvLaB7Wv4ngBJy2VSrzvbmB1IkiRJUotMAiRJkqSOMQmQJEmSOsYkQJIkSeoYkwBJkiSpY0wC\nJEmSpI4xCZAkSZI6xiRAkiRJ6hiTAEmSJKljTAIkSZKkjjEJkCRJkjrGJECSJEnqGJMASZIkqWNW\ntx2AJKl9s7PZfMUVsGZN2g5lZLOzqbZjkKRJY0+AJAlgZs2atkMYi3vbDkCSVgKTAEkSwKFtBzAG\n9wIzbQchSSuBw4EkSftZv36SRteEqpq8MUyS1DJ7AiRJkqSOMQmQJEmSOsYkQJIkSeoYkwBJkiSp\nY0wCJEmSpI4xCZAkSZI6xiRAkiRJ6hiTAEmSJKljTAIkSZKkjjEJkCRJkjrGJECSJEnqmNaSgCQb\nktyc5JYkF8yx/nFJPpXk+0k29y1fm+TjST6f5HNJXjneyCVJkqTJtrqNgyZZBVwInA7cBlyb5PKq\n2tZX7W7gPOD5A5vvBl5VVVuTHAr8Q5KPDmwrSZIkaR5t9QScCmyvqlurajewBTirv0JV3VVV19Fc\n9Pcvv6OqtvY+3wtsA44aT9iSJEnS5GsrCTga2NFX3tlbNpIkxwInA1cvSVSSJElSB7SVBNRid9Ab\nCnQZcH6vR0CSJEnSEFp5JoDmOYC1feW1NL0BQ0lyCPAB4H1V9efz1JnpK85W1ezoYUqSJEnTp60k\n4DrghN5wntuBs4GN89TNfoUkwMXATVX15vkOUFUzSxGoJEmSNG1aSQKqak+STcBVwCrg4qraluTc\n3vqLkjwCuBY4DLgvyfnAScCTgRcDn01yfW+Xr6+qD4+9IZIkSdIEStWih+evOEmqqnLwmpKWWpJa\ngsd+xigM+/di8toGw7Zvdjb7GrZ+/SS1cfjzJ0nLYVKvO31jsCRJktQxJgGSJElSx5gESJIkSR1j\nEiBJkiR1jEmAJEmS1DEmAZIkSVLHmARIkiRJHWMSIEmSJHWMSYAkSZLUMSYBkiRJUseYBEiSJEkd\nYxIgSZIkdczqtgOQuiZJtR3DqKoqbccgLdYk/u6Bv3+Sloc9AVIraoK+pGnS9u+Tv39SFyTZkOTm\nJLckuWCeOuuSXJ/kc0lmxxyiPQGSJEnSUkmyCrgQOB24Dbg2yeVVta2vzuHA/waeXVU7kxwx7jjt\nCZAkSZKWzqnA9qq6tap2A1uAswbq/CLwgaraCVBV3xhzjCYBkiRJ0hI6GtjRV97ZW9bvBOBhST6e\n5LokvzS26HocDiRJkiQtnWEe6DkEOAV4JrAG+FSST1fVLcsaWZ+p7QlIsq7/s2XLK6XcmB34vJLL\no/0+/WgfbcU7annU9rUd76jl0dq3devefbQV76hlz59ly5bHW+59zSS5NMml3N9twNq+8lqa3oB+\nO4CPVNX3qupu4G+AJ82xr2WTqumbfSBJOaWaVqokNVmzfmSkKQqnuX2T1zYYtn2zsz+aPnP9+klq\no+dPUrsGrzuTrAa+QHOX/3bgGmDjwIPBj6N5ePjZwAOBq4Gzq+qmccXtcCBJkiRpiVTVniSbgKuA\nVcDFVbUtybm99RdV1c1JPgx8FrgPeOc4EwCwJ0Aau8m7G2lPwL6aE9c2sCegr+YUnz9J7ZnU686p\nfSZAkiRJ0twcDiRJQ5idzeYrroA1aybuZs9+d/klSQJ7AiRpWDNr1rQdwvLbtevQtkOQJI2BSYAk\nDWfqr4537YJLL51pOwxJ0hg4HEiSRjTdD85uXu6AJEkrgD0BkiRJUse0kgQk2ZDk5iS3JLlgjvWP\nS/KpJN9PsnmUbSVJkiQd2NiTgCSraN6QtgE4CdiY5MSBancD5wG/u4BtJUmSJB1AGz0BpwLbq+rW\nqtoNbAHO6q9QVXdV1XXA7lG3lSRJknRgbSQBRwM7+so7e8uWe1tJkiRJtJMELGZajUmakkOSJEla\nkdqYIvQ2YG1feS3NHf0l3TbJTF9xtqpmhw9RkiRJml5tJAHXASckORa4HTgb2DhP3cG5rYfetqpm\nFh2pJEmSNIXGngRU1Z4km4CrgFXAxVW1Lcm5vfUXJXkEcC1wGHBfkvOBk6rq3rm2HXcbJEmSpEmW\nqukbZp+khn1DpjRuzVtZJ+n3bvg3zsL0tm92NvsaNd1vDJ6ktoHtk9S2Sb3u9I3BkiRJUseYBEiS\nJEkdYxIgSZIkdYxJgCRJktQxQ88OlOQhwIuAn6aZmedBwH3AvcCngT+tqvuWI0hJkiSp65byenyo\nJCDJs4CTgL+oqncMrAvwJODVSf6qqrYO2xBJkiRJB7fU1+MHTQKSPAj4SlV9dK711cwxuhXYmuQJ\nQ7RBkiRJ0pAGr8eT/BTwj1X1PVjY9fhBnwmoqu9X1fa+IH4qyYPnqXvjMAeVJEmSNJzB63FgM/AU\ngCT/KsnP99Ud6np8IQ8Gz3tQSZIkScvuGuC4JMdV1d8BDx91BwtJAhZ9UEmSJEkLthb4Ac0zAB8H\nfmbUHSwkCVj0QSVJkiQt2JeBD1TVecC/B7426g4WkgQs+qCSJEmSFuyPgcf3Ph8HHDnqDhaSBCz6\noJIkSZKGk+SBSY7YW66qH1bVZ3qfr62q3+ir+8hh9nnQJGA5DipJkiRpOFX1A+C0JL843yydSR6a\n5BXAo4bZ50HfE1BVP0jyrCSHAX+2dz7SwYPSDA3ahsODJEmSpCVVVX+R5CeBVyV5OM3bgg8Bfgjs\nAnYC76yqbw+zv6HeGLzUB5UkSZI0mqr6R+A3l2JfQyUBS31QSZIkSe1ZyIPBkiRJkiaYSYAkSZI0\nwZKck2RVkmOH3WbRScBCDipJkiRpyTwAeBJwzLAbLEVPwMgHlSRJkjS8JK9N8stJnjfH6q8BRwFP\nG3Z/QyUBSZ5/gDv9Ix9UkiRJ0kieALy3qj40x7qvAlcCQ8/UOWxPwNOBIwCSnLXYg0qSJEkaySeq\nqpIcnuTsJM9O8jCAqrqpqu6rqrcNu7Nhpwj9EPDGJA8CHpzkMcCNwI1VdVOvztAHlSRJkjSSAFTV\nPUn+BHgHcDvwzYXsbNiXhX0M+BhAks3AdcDjgTOTHEXzsrA/qKovLCQISZIkSQf0P5KcAWwFPgvc\nXFU3AiR5SlVdPcrOUlWLjijJC4G1VfU7i97ZEkhSVZW245DmkqRg8b934xNG+X2a1vbNzmZfo9av\nn772wSSeO7B9kto2ruvOJOcDVwOnAT8HPIVmgp5PAA+pqheMsr+h3xh8EP8M3LxE+5IkSZLUp6re\n0vv46b3LkhwBnApsGnV/S5IEVNUHl2I/Euy9Wzd5vFsnqU2T+LfTv5vS4lTVN4Ark4z8XEArbwxO\nsiHJzUluSXLBPHXe2lt/Q5KT+5a/Psnnk9yY5A+TPHB8kWt8asK+JGklaPtvoX83pTZU1acPXmt/\nY08CkqwCLgQ2ACcBG5OcOFDnDOD4qjoBeAXw9t7yY4GXA6dU1ROAVcALxxa8JEmSNAXa6Ak4Fdhe\nVbdW1W5gCzD47oEzgfcA9J50PjzJkcB3gN3AmiSrgTXAbWOLXJIkSZoCbSQBRwM7+so7e8sOWqeq\nvgm8ieYtxbcD91TVXy1jrJIkSdLUaSMJGHYg4P0eFkryaOBXgWOBo4BDk7xo6UKTJEmSpt9STRE6\nituAtX3ltTR3+g9U55jesnXAJ6vqboAkHwSeCrx/8CBJZvqKs1U1u8i4JUmSpKnQRhJwHXBC7yHf\n24GzgY0DdS6nme90S5LTaIb93JnkC8B/TfJg4PvA6cA1cx2kqmaWJXpJkiRpwo09CaiqPUk2AVfR\nzO5zcVVtS3Jub/1FVXVlkjOSbAe+C5zTW7c1yXtpEon7gM8A7xh3GyRJkqRJlqrpm6t3XK9v1vJo\nXngzaT+XGfqlN5PXvuHbBtPbvtnZH72Iaf366WsfTOK5A9vXV3Pi2jfa3xZppZrU685WXhYmSZIk\nqT0mAZIkSVLHmARIkiRJHWMSIEmSJHVMG1OESppSs7PZfMUVsGbNZD0f1f/QryRJXWBPgKSlNLNm\nTdshLK9duw5tOwRJkhbNJEDSUprqK+Rdu+DSS2faDkOStMIl2ZDk5iS3JLngAPV+LsmeJC8YZ3zg\ncCBJy2Ry5tIfdR72zcsdkCRpgiVZBVwInA7cBlyb5PKq2jZHvd8GPgyMfRytPQGSJEnS0jkV2F5V\nt1bVbmALcNYc9c4DLgPuGmdwe5kESJIkSUvnaGBHX3lnb9k+SY6mSQze3ls09u5zkwBJkiRp6Qxz\nQf9m4HVVVTRDgcY+HGhqnwlIsq6qZvd+BrA8GeXGLLCu7zMTUG5MX/tG+33aupUBbcd/sPJo7Ws/\n3lHLtq+/3H68o5anu31t/39j2fJCyj3rgGOZ223A2r7yWpregH4/A2xJAnAE8Jwku6vq8nn2ueTS\nJCDTJUkN+6CfVp7m4ctJ+7kc9eHSSWrf8G3rn29/eh8MnpR27WX79tW0fSvM8G2TVrLB684kq4Ev\nAM8EbgeuATYOPhjcV/8S4ENV9cFxxLvX1PYESJIkSeNWVXuSbAKuAlYBF1fVtiTn9tZf1GqAPfYE\naMWZvLtZ4N26hj0BK5Ht21fT9q0w9gRoOkzqdacPBkuSJEkdYxIgSZIkdYxJgCRJktQxJgGSJElS\nx5gESJIkSR1jEiBJkiR1jEmAJEmS1DEmAZIkSVLHmARIkiRJHWMSIEmSJHWMSYAkSZLUMSYBkiRJ\nUseYBEiSJEkd00oSkGRDkpuT3JLkgnnqvLW3/oYkJ/ctPzzJZUm2JbkpyWnji1ySJEmafGNPApKs\nAi4ENgAnARuTnDhQ5wzg+Ko6AXgF8Pa+1W8BrqyqE4EnAtvGErgkSZI0JdroCTgV2F5Vt1bVbmAL\ncNZAnTOB9wBU1dXA4UmOTPITwNOq6t29dXuq6ttjjF2SJEmaeG0kAUcDO/rKO3vLDlbnGOA44K4k\nlyT5TJJ3JlmzrNFKkiRJU6aNJKCGrJc5tlsNnAK8rapOAb4LvG4JY5MkSZKm3uoWjnkbsLavvJbm\nTv+B6hzTWxZgZ1Vd21t+GfMkAUlm+oqzVTW78JAlSZKk6dFGEnAdcEKSY4HbgbOBjQN1Lgc2AVt6\ns//cU1V3AiTZkeQxVfVF4HTg83MdpKpmliV6SZIkacKNPQmoqj1JNgFXAauAi6tqW5Jze+svqqor\nk5yRZDvNkJ9z+nZxHvD+JA8AvjSwTpIkSdJBpGrYIfqTI0lV1eAzBZoQSWr4R0dWijDsz9zktW/4\nts3OZl/okqnaAAANCklEQVTD1q+flDZO87kD29dX0/atMMO3TVrJJvW60zcGS5IkSR1jEiBJkiR1\njEmAJEmS1DEmAZIkSVLHtDFFqBapefhr8kziQzOSpJVhEv/v8/89rWQmARNr0v4W+ndQkrRYk/R/\nn//vaWVzOJAkSZLUMSYBkiRJUseYBEiSJEkdYxIgSZIkdYxJgCRJktQxJgGSJElSx5gESJIkSR1j\nEiBJkiR1jC8Lk8Zodjabr7gC1qyZrJfIzM5O3ps6JUnS/OwJkMZrZs2atkNYfrt2Hdp2CJIk6QBM\nAqTxmvqr41274NJLZ9oOQ5IkHYDDgaSWrF8/KSNsQlUNPX4pScHm5QxIkiQtkj0BkiRJUseYBEiS\nJEkdYxIgSZIkdYxJgCRJktQxJgGSJElSx5gESJIkSR1jEiBJkiR1jEmAJEmS1DEmAZIkSVLHmARI\nkiRJHWMSIEmSJHVMK0lAkg1Jbk5yS5IL5qnz1t76G5KcPLBuVZLrk3xoPBFLkiRJ02PsSUCSVcCF\nwAbgJGBjkhMH6pwBHF9VJwCvAN4+sJvzgZuAWv6IJUmSpOnSRk/AqcD2qrq1qnYDW4CzBuqcCbwH\noKquBg5PciRAkmOAM4B3ARlb1JIkSdKUaCMJOBrY0Vfe2Vs2bJ3fB14D3LdcAUqSJEnTbHULxxx2\nCM/gXf4keS7w9aq6Psm6A26czPQVZ6tqdugIJUmSpCnWRhJwG7C2r7yW5k7/geoc01v2C8CZvWcG\nHgQcluS9VfWSwYNU1cxSBi1JkiRNizaGA10HnJDk2CQPAM4GLh+ocznwEoAkpwH3VNUdVfWGqlpb\nVccBLwQ+NlcCIEmSJGl+Y08CqmoPsAm4imaGnz+uqm1Jzk1ybq/OlcCXk2wHLgL+03y7G0fMkiRJ\n0rAONh1+khf1psH/bJK/T/LEscdYNX3X0UmqqqZ25qAkNXn5Txj2nExz+2Zns69h69dPShuHP3cw\niefPn819NW3fCjTN7fNvi6bD4HVnbzr8LwCn0wxnvxbYWFXb+ur8S+Cmqvp2kg3ATFWdNs64fWOw\nJEmStHQOOh1+VX2qqr7dK15N8/zrWJkESJIkSUtnmOnw+/1H4MpljWgObcwOJEmSJE2rocetJVkP\nvAz4+eULZ25T2xPQ/x6BJOumqdyYHfi80ss/0vX2bd3afLUb7yjl0X6ffrSPtuIdtTxq+9qOd9Sy\n7bN9K7c8+v9/7ca73O2zPBnl3tdMkkuTXMr9DTMdPmkeBn4ncGZVfWuO/SwrHwyeQJm4h6Nguh9u\nAx8M7qs9cefPn819NW3fCjTN7fNvi6bD4HVnktU0DwY/E7gduIb7Pxj8SOBjwIur6tNjDhlwOJAk\nSZK0ZKpqT5K90+GvAi7eOx1+b/1FwK8BDwXengRgd1WdOs447QmYQJN3NwSm+24W2BPQV3vizp8/\nm/tq2r4VaJrb598WTYdJve6c2mcCJEmSJM3NJECSJEnqGJMASZIkqWNMAiRJkqSOMQmQJEmSOsYk\nQJIkSeoY3xOgFWV2NpuvuALWrJm4mbb2m/5TkiRpJbMnQCvNzJo1bYew/HbtOrTtECRJUoeZBGil\nmfqr41274NJLZ9oOQ5IkdZjDgbRiTc4bdWH0t3puXu6AJEmS5mVPgCRJktQxU9sT0NxtnSzD3kmW\nJEnTYxKvWcDrlkk3tUkATNrvk79HkiR1l9ctGi+HA0mSJEkdYxIgSZIkdYxJgCRJktQxJgGSJElS\nx5gESJIkSR1jEiBJkiR1jEmAJEmS1DEmAZIkSVLHtJYEJNmQ5OYktyS5YJ46b+2tvyHJyb1la5N8\nPMnnk3wuySvHG7kkSZI02VpJApKsAi4ENgAnARuTnDhQ5wzg+Ko6AXgF8Pbeqt3Aq6rq8cBpwK8M\nbitJkiRpfm31BJwKbK+qW6tqN7AFOGugzpnAewCq6mrg8CRHVtUdVbW1t/xeYBtw1PhClyRJkiZb\nW0nA0cCOvvLO3rKD1Tmmv0KSY4GTgauXPEJJkiRpSrWVBNSQ9TLfdkkOBS4Dzu/1CEiSJEkawuqW\njnsbsLavvJbmTv+B6hzTW0aSQ4APAO+rqj+f+xAzfZ/X9b4kSZIktZUEXAec0BvOcztwNrBxoM7l\nwCZgS5LTgHuq6s4kAS4GbqqqN89/iJklD1qSJEmaBq0kAVW1J8km4CpgFXBxVW1Lcm5v/UVVdWWS\nM5JsB74LnNPb/OeBFwOfTXJ9b9nrq+rDY26GJEmSNJFSNezw/MmRpIZ/7GClCFU1+AzE3DWnuH2z\ns9nXsPXrJ6mN03z+hm8bTHf7Jq9tYPv6atq+Fca/LftqTlzbYNTzN82S1CR+L3xjsCRJktQxJgGS\nJElSx5gESJIkSR1jEiBJkiR1TFtThGqBZmez+YorYM2aiXv+ZL+HfiVJktQeewImz8yaNW2HsPx2\n7Tq07RAkSZKmlknA5Jn6q+Ndu+DSS2faDkOSJGlqORxogk33PPqblzsgSZKkzrInQJIkSeoYkwBJ\nkiSpY0wCJEmSpI4xCZAkSZI6xiRAkiRJ6hiTAEmSJKljTAIkSZKkjjEJkCRJkjrGJECSJEnqGJMA\nSZIkqWNMAiRJkqSOMQmQJEmSOsYkQJIkSeoYkwBJkiSpY0wCJEmSpI4xCZAkSZI6xiRAkiRJ6hiT\nAEmSJKljUlVtx7DkkhRMWrtCVWWomrZvBZrm9g3fNpju9k1e28D29dW0fSuMf1v21Zy4tsGo52+a\nJalJ/F7YEyBJkiR1TCtJQJINSW5OckuSC+ap89be+huSnDzKtpIkSVJbFnOtOy5jTwKSrAIuBDYA\nJwEbk5w4UOcM4PiqOgF4BfD2YbeVJEmS2rKYa91xaqMn4FRge1XdWlW7gS3AWQN1zgTeA1BVVwOH\nJ3nEkNtKkiRJbVnote6R4wyyjSTgaGBHX3lnb9kwdY4aYltJkiSpLQu91j1mmePaTxtJwLCPv0/c\nU9aSJEnqvIVe6451iqjV4zxYz23A2r7yWprs50B1junVOWSIbXsmL4dopggbuvbyBbJMbN9+tZcv\nkGUwWttguts3WW0D2zdQe/kCWSbT3D7/tuxXe/kCWSajn7/OWOi17m3LHNd+2kgCrgNOSHIscDtw\nNrBxoM7lwCZgS5LTgHuq6s4kdw+xrfPWSpIkqS0LvtYdZ5BjTwKqak+STcBVwCrg4qraluTc3vqL\nqurKJGck2Q58FzjnQNuOuw2SJEnSXBZzrTtOU/nGYEmSJEnzm6o3BidZm+TjST6f5HNJXtl2TBpN\nklVJrk/yobZj0WiSHJ7ksiTbktzU697UhEjy+t7fzhuT/GGSB7Ydk+aW5N1J7kxyY9+yhyX5aJIv\nJvlIksPbjFHzm+f8/U7vb+cNST6Y5CfajFFzm+vc9a3bnOS+JA9rI7aFmKokANgNvKqqHg+cBvyK\nLxObOOcDNzHmJ+S1JN4CXFlVJwJPBByqNyF641ZfDpxSVU+g6b5+YZsx6YAuoXkJUb/XAR+tqscA\nf90ra2Wa6/x9BHh8VT0J+CLw+rFHpWHMde5IshZ4FvDVsUe0CFOVBFTVHVW1tff5XpqLkKPajUrD\nSnIMcAbwLiZxmoQO6921elpVvRua8ZBV9e2Ww9LwvkNzE2VNktXAGsY8S4WGV1V/C3xrYPG+Fw/1\n/n3+WIPS0OY6f1X10aq6r1e8mjHPF6/hzPO7B/B7wGvHHM6iTVUS0K93Z+tkml8mTYbfB14D3Hew\nilpxjgPuSnJJks8keWeSNW0HpeFU1TeBNwFfo5nJ4p6q+qt2o9KIjuybWeROYKxvHtWSehlwZdtB\naDhJzgJ2VtVn245lVFOZBCQ5FLgMOL/XI6AVLslzga9X1fXYCzCJVgOnAG+rqlNoZjpwOMKESPJo\n4FeBY2l6Tw9N8qJWg9KCVTPjh0MqJ1CSNwL/XFV/2HYsOrjeza43AP+tf3FL4Yxs6pKAJIcAHwDe\nV1V/3nY8GtpTgTOTfAX4I+AZSd7bckwa3k6aOyHX9sqX0SQFmgw/C3yyqu6uqj3AB2l+JzU57kzy\nCIAkPwl8veV4NKIkL6UZEmsCPjkeTXPz5Ibe9csxwD8keXirUQ1pqpKAJAEuBm6qqje3HY+GV1Vv\nqKq1VXUczQOJH6uql7Qdl4ZTVXcAO5I8prfodODzLYak0dwMnJbkwb2/o6fTPKCvyXE58Mu9z78M\neBNsgiTZQDMc9qyq+n7b8Wg4VXVjVR1ZVcf1rl920kywMBFJ+FQlAcDPAy8G1vemmby+94ulyWNX\n9uQ5D3h/khtoZgf6zZbj0ZCq6gbgvTRvudw7rvUd7UWkA0nyR8Angccm2ZHkHOC3gGcl+SLwjF5Z\nK9Ac5+9lwB8AhwIf7V27vK3VIDWnvnP3mL7fvX4Tde3iy8IkSZKkjpm2ngBJkiRJB2ESIEmSJHWM\nSYAkSZLUMSYBkiRJUseYBEiSJEkdYxIgSZIkdYxJgCRJktQxJgGSJElSx/z/USQUs8H5S0gAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa759e817d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# pmf and cdf for a pair of rolled dice.\n",
    "\n",
    "def px(s):\n",
    "    s1 = s-1\n",
    "    s13 = 13-s\n",
    "    return np.min([s1,s13]) / 36.\n",
    "\n",
    "def Fx(pmf):\n",
    "    \"\"\"\n",
    "    $F_x(x) = \\textstyle \\sum_{x}^x_i x_i`\n",
    "\n",
    "    all that is really required:\n",
    "    return np.cumsum(pmf)\n",
    "    \"\"\"\n",
    "    cum = 0\n",
    "    cumulative = []\n",
    "\n",
    "    for fx in pmf:\n",
    "        cum += fx\n",
    "        cumulative.append(cum)\n",
    "\n",
    "    return cumulative\n",
    "\n",
    "#sample space\n",
    "S = np.arange(2,13)\n",
    "\n",
    "pmf = [px(s) for s in S]\n",
    "\n",
    "\n",
    "## Plot it up..\n",
    "f, ax1 = plt.subplots(figsize=[12,6])\n",
    "\n",
    "## Axis 1, pmf\n",
    "ax1.bar(S, pmf)\n",
    "ax1.set_ylabel('$f_x(x)$')\n",
    "\n",
    "## CDF\n",
    "ax2 = ax1.twinx()\n",
    "ax2.step(S+1,Fx(pmf), 'y', lw=3)\n",
    "ax2.set_ylim([0,1.])\n",
    "ax2.set_ylabel('$F_x(x)$')\n",
    "\n",
    "## Title\n",
    "plt.title('$pmf$ and $cdf$ for two rolled dice example')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative Distribution Function -- discrete case\n",
    "\n",
    "Denoted $F(x) = P(X \\le x)$, we can use the *cumulative distribution\n",
    "function*, $c.d.f$ to evaluate the probability properties of a random\n",
    "variable.\n",
    "\n",
    "In the example above, we plotted the $c.d.f$ along with the $p.m.f$.\n",
    "Note, that all is required for the $c.d.f.$ is to calculate the\n",
    "cumulative sum (e.g. cumsum) of the $p.m.f.$.\n",
    "\n",
    "Knowledge of either the $c.d.f.$ or the $p.m.f.$ of a random variable\n",
    "allows the other to be calculated:\n",
    "\n",
    "$$F(x) = \\sum_{y=y\\le x}{P(X = y)}$$\n",
    "\n",
    "where $P(X=y)= p_i$\n",
    "\n",
    "The $c.d.f.$ is an increasing step function, with steps at each value of\n",
    "the random variable. The heights of each step are simply the\n",
    "probabilities associated with each value. We can obtain the $p.m.f.$\n",
    "from the $c.d.f.$ through the relationship:\n",
    "\n",
    "$$P(X = x) = F(X) - F(X^-)$$\n",
    "\n",
    "Where $F(X^-)$ is the limiting value (e.g. 2 in the dice example). So at\n",
    "that point, $F(X) = F(X^-)$ and the $c.d.f.$ is $0$.\n",
    "\n",
    "### Probability Density Function\n",
    "\n",
    "The analogue to $p.m.f.$ for a *continuous* random variable is the\n",
    "$p.d.f.$, or **probability density function**. Continuous random\n",
    "variables are members of the set of real numbers ${\\rm I\\!R}$ and $x$\n",
    "may take any value within a continuous region.\n",
    "\n",
    "Some examples of continuous random variables are instantaneous\n",
    "streamflow, soil moisture values, temperature and humidity. We can still\n",
    "define limits to the random variable, so that they may be within\n",
    "$0 \\le x \\le \\infty$ for example.\n",
    "\n",
    "Note that many continuous random variables are measured at discrete\n",
    "intervals. So, for instance, you might have daily precipitation\n",
    "measurements. But the *value* of the precipitation is what falls on the\n",
    "continuous interval. But this can be further confused by the fact that\n",
    "some measurement devices can only provide a value to a certain\n",
    "sensitivity. Often times continuous random variables are measured in\n",
    "discrete terms (e.g. think about a compass heading). What truly\n",
    "distinguishes a discrete from continuous random variable is how we\n",
    "define the probablistic properties.\n",
    "\n",
    "For continuous random variables, a $p.d.f.$ is used. This function\n",
    "$f(x)$ is must evaluate to a positive value and integrate to 1 over the\n",
    "whole sample space. So we have:\n",
    "\n",
    "$$f(x) \\ge 1$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\int_{\\S} {f(x)dx} = 1$$\n",
    "\n",
    "To find the probability that the value of interest lies within the\n",
    "region of interest, we integrate teh $p.d.f.$ between the two values.\n",
    "\n",
    "$$P(a \\le X \\le b) = \\int_{a}^b{f(x)dx}$$\n",
    "\n",
    "**NOTE** Perhaps the most important way to distinguish between discrete\n",
    "and continuous random variables is that fact that a **continuous random\n",
    "variable cannot take a specific value of x**. Put another way, the\n",
    "probability that $X$ takes the value $a$ is always $0$!\n",
    "\n",
    "Let's look closer:\n",
    "\n",
    "$$P(X=a) = \\int_{a}^a{f(x)dx} = F(a) - F(a) = 0$$\n",
    "\n",
    "Why? We can work it out with an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEWCAYAAACt5MYgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FeXZBvD7mZmc7IEEMJBEZBFkCYbFKC4VvqKoiC0u\nxQUQbBVUkMUqSyuyChUFLVKxSEAUZdGiuABiCwgBCsYFQkD2yCb7FpCAkPf745ykAZIQkpO8Z/Le\nv+vq1eRkzsxDTG4ennlnRpRSICIid7B0F0BERMXH0CYichGGNhGRizC0iYhchKFNROQiDG0iIhdh\naBMRuQhDm4jIRRjaZAQRyRSRNrrrICothjaZQvn+R+RqDG0iPxOR/xORPSJyZRHb9BORl0SkR3nW\nRu7H0KZSE5E4EfmXiOwXkW0i8ozv9boickhEmuXb7oCI3Or7fKCIbBGR4yKSISId8u0zU0SeE5G1\nIpIlIikiEisi80XkmIh8JSKVL9h+oG8/h0VkiogEX069frQMwD6l1M5Cjl8JQEcAcwGk+vnYVMEx\ntKlURMQC8BmA7wHEAWgDoK+ItFVKbQUwAMB0EQkFMBXAVKXUUt/btwC4RSkVBWCYb7tY39cUgPt8\n+7sGQHsA8wEMBHAFvD+7vS8o5xEAbQHUBVAfwAuXU28pvxX5tQDwTRFfvwHAD0qp1UqpDD8elwzA\n0KbSSgZQVSk1Uil1Vim1HcBkAA8BgFJqMrzhvBpALIC/5r5RKfWRUmqv7+PZADYDuD7fvt9QSh1Q\nSu2Bt3tdqZRao5Q6DeBjAM3ybasATFBK7VZKHQHwEoCHL7fe/EQkSUT+KCIvi8jvReQJEXm0oG+C\niNQXkREicieAQQD+U8h2NwDoA8ARkXsL2oaoKI7uAsj1rgIQJyJH8r1mA1ia7/PJ8I4CnlBK/Zr7\noi8A+wGo5XspAkDVfO/bl+/jUxd8nu3bPr/844gd8HbSJak3VyyAjQDuUEoNEJFweDv0d/Nv5Ht9\nNoDWSqmjItIfwKIC9gel1CoROQXgdXbZVBIMbSqtHQC2K6XqF/RFEYkA8Dq8wT1MROYopY6IyFUA\nJgH4LbwdtBKR7wFIEccq6msAUPOCj/dcbr35KaUWisgweMcpgLezP1jApvcBSPcFdgiACKXUgSJ2\n3RDA+ksdn6ggHI9Qaa0GkCUi/UUkVERsEUkUket8X/87gNVKqe4AvgDwlu/1cHhHGgcBWCLyGIDE\nUtQhAJ4WkXgRiYF3DDOzBPVe6DYAX/s+7grgVSDvJGvuXyJVAazxfXw7gP/6xiQXbgffzP6gyvf0\nkQu3ISpKkaEtIsNEpE++z18SkQtP/pDBlFI58J4kbApgG4AD8HbQUSLye3hPDD7l2/xZAM1F5GGl\n1HoAYwGsBLAX3sC+1EoKdcHHF37+AYCFALbCOx8feTn1Xritb5VHDIDfisgTAFYppeb4vvyhbx8A\nMANAgojcBW+A/wogsoDtAO9JyOUXHOrCbYgKJUU9bsz3T9g5SqkWvrPumwAk+070EAUMEdkO4E9K\nqQJnySXc570AWiqlBvhhXy0APAHgMIBZSqk1l3gLUYGKnGkrpX7yrbNtCqA6gO8Y2GQCEWkA778M\ntohIlFLqeCl3eQ7ALgC/MLCpNIpzInIygMfgPZM+pWzLIQoMSqkfAfzGj/v7AcAP/tofmavI8QgA\niEgQgHXwLouqpy71BiIiKjOX7LSVUr+KyCIARxjYRER6XTK0fScgWwJ44BLbMdCJiEpAKVXsJZ9F\nhraINIL3woI5vvtI+O3AuojIUKXUUN11ELkVf4f863Ib3kutHlkP7813iIgoAPCKSCIiFzExtJfo\nLoDI5ZboLsBkl1zyV+wdiSg3zLSJiALJ5WaniZ02EZFrMbSJiFyEoU1E5CIMbSIiF2FoExG5CEOb\niMhFGNpERC7C0CYichGGNhGRizC0iYhchKFNROQiDG0iIhdhaBMRuQhDm4jIRRjaREQuwtAmInIR\nhjYRkYswtImIXIShTUTkIgxtIiIXYWgTEbkIQ5uIyEUY2kRELsLQJiJyEYY2EZGLMLSJiFyEoU1E\n5CIMbSIiF2FoExG5CEObiMhFGNpERC7C0CYichGGNhGRizC0iYhchKFNROQiDG0iIhdhaBMRuQhD\nm4jIRRjaREQuwtAmInIRhjYRkYswtImIXIShTUTkIgxtIiIXYWgTEbkIQ5uIyEUY2kRELsLQJiJy\nEYY2EZGLMLSJiFyEoU1E5CIMbSIiF2FoExG5CEObiMhFGNpERC7C0CYichGGNhGRizC0iYhchKFN\nROQiDG0iIhdhaBMRuQhDm4jIRYwKbRGJ1F0DEVFpGBXaljX9G9v+Nktk8psiw6N110NEdLkMC+24\nyCuvPBQRGZnYHeh52HE+Wycy5n7ddRERFZdRoa2UJS1bnsDixUF2SsoaNGkijUQ6f2jba06JpLwr\nMqiG7hqJiIpiVGgDkvdRUlIlTJ5cQ1JTd8njj+8IiY5u8AjQf4/jfLFF5OXHRIZKETsiItLCsNAG\nRNR5nwcH2+jePQ5ffRVif/DBWiQnqzqW9XCKSMdsy5r6kcig2ppKJSK6iGGhbcmFoZ1f/foRmDCh\nhqxcuV96997uqVatbgdgwDbH+XKnyCu9RToa9v0iokBjWAgVb+Jh24IuXapj3rxwe86cdbjlljPx\njnPfayIvnbasafNEBjQu40KJiApkWmgX2WkXpGbNUIwdGycrVhyxBg7c4sTF1WwLDFrnOP/ZKzJ2\nkEgPp4yKJSK6iGmhDSnh6UXLEtx//xWYOzfK/vzz9WjT5pfYoKD2I0UGnrbt9xeJ/PU6/9ZKRHQx\nA0P78jrtglSvHoxRo+KwYsUJa9iwTdaVV1a/FXj+G8dZclDk9REijwX7oVgioosYFtqWAKUP7Vwi\nQLt2VfGvf1W2Fy78Ee3aHa8SEtJmEDDilG3PXC4y9Ba/HYyICMaFNko8HrmUmJggDBkSj9TUM/ar\nr26QOnViWgK9lznO8qMiE14R6RVeNkcmIpMYFtpFL/nzl9atYzBzZhVr8eLNuPfeQ5XCwm7uB7yQ\nZdsfpYm8dHuZF0BEFZZhoS3w53jkUiIjHQwcGI+lS3Ps8ePT5ZprwpsB3Rfa9jdZIpMmiAyvVG7F\nEFGFYFxol0enXZCbborBe+9dYX399TY89NDuiIiIpCeBp4/a9tw1ImPu0VIUEbmOaaEtZTXTLq7w\ncAfPPnsllixx7EmTfkBiopUo0mWubX93QmTyWyIvxOitkIgCmVGhrVT5jkcupXnzaEyZEmctW7ZD\nHn10Z3hUVOPHgX6HHGduusjo3+uuj4gCj1GhXZIrIstDSIiDXr0SsGiRx548+QckJlqNRR772La/\nOykyaZJI/yq6aySiwGBaaJfZkj9/adq0MlJS4mTZsh3SteuusKioa/8I9D/oOJ+ki4y8V3d9RKSX\naaEdkJ12QUJCbPTsGY9Fi4LslJS1aNIkqLHI4//ydt9vsfsmMpRpoY1AmmkXl/eBDdV93ffusKio\nZr7u++N1IsM76K6PiMqPcaEd6OORoni77zgsWuTYKSnpaNIkuJFIjzm2nXZS5K23RXpz5QlRBWda\naLtmPHIpSUlRmDw5VlJTd0jXrnvDoqKaPQa8eIjdN1HFZlpoa7u4pqwEB1vo2bNGId33xEnsvokq\nFtNCu8J02gU5v/v+OSwqqsUffd13OrtvoorBsNA244/r7b7jsGiRbU+evBaJiSGN882+/ynyXLTu\nGomoZMxIsf/Rfhl7eWvatBJSUq6QZct+ki5dfg6LjGz+J2DgYdv+OF1kBO95QuQyhoV2xZtpF1dI\niI1nnonD4sW53XdwY5Eec2077YTIP98S+Wtl3TUS0aUZFtoW3LhO29+aNq2EKVNiZenSTOnUaU94\nZGSzx4F+R2z7kx9EXm6nuz4iKpxhoY0KfSLycoWG2ujTJx6LF9v2P/+5Bo0aBTUReewL2/42S+Tt\nf/B+30SBx7DQtlx9cU1ZatGiMt55p7q1dGkmHnlkV0RERFIP4Okjtj33O5ExbXXXR0ReRoW2UuLX\nB/tWRKGhNvr2TcCSJbY9ceL30rChkwR0+9K2Vx8TeWusSP8w3TUSmcyo0AZg7InIkkhOjsG0adWt\nJUu24Q9/+DkqLOz6PsCfT9j2RytFhvNJ80QaGBbalnFL/vwhIiIIzz+fgKVLlf3aa2ukbt2IZKDX\nMsdJPSzyxkiRPh7dNRKZwrDQNnfJn7/85jdVMGNGNfs//9mEe+45FB0S0nog8Hy2bc9YIjK4he76\niCo640KbM23/qFQpCC+8kIDU1DP26NHrpWbNarcAf05znMX7Rcb9RaSHo7tGoorItNDmeKQM3H57\nFXz0UbS9YMEGtG2bVS0oqP0IkQHZlvXuApG/NNJdH1FFYlho80RkWapa1YMRI+KxYkWWNXjwJjsu\n7srbgP4ZjvPlbpGXe4t0NO7njcjfDPslsrjkrxyIAL/7XTXMnRtlz52bgVatfo2z7T+8LjL8tGWl\nfCzyfC3dNRK5lWGhzROR5S0+PgRjxtTAypWH5bnntjnVqtW7B/jLdseZt01k1GMcWBFdHtNCmxGh\niWUJHnwwFvPmhduzZ6/DDTfk1LKsLimWteaUyNvTRf5cXXeNRG5gWmiD4xH96tQJw/jxNWTFin3y\n1FM7g2NiGj8E/PVnx/l0g8hLD+iujyiQGRbaPBEZSBxH8Nhj1bFwYbD97rvpSEqyrxH502zb/vYk\nbxdLVDDDQptXRAaqRo0iMGlSdUlN3SmPPronzHu72L6HbXvO9yIv3aG7PqJAYVhoczwS6IKDLfTq\n5X1Yw1tvrZWGDUOvBbovsO1Vx0TefEVkWKjuGol0Mi20xTLsT+xm111XGdOmXWEtWbINDzywNyo0\ntGU/oMdJ2/4wVeRvLXXXR6SDYRHGTtuNIiIc9O+fgGXLlD127BqpXTuyJdB9peOkHhCZMJiXzJNJ\nDAttnoh0u1atqmLWrGr2woU/4s47D1T1eP5vqMigbMv6YIHIkAa66yMqa0aFtlI8EVlRxMQEY9iw\nmli+PNsaPDjDrlEj9jag7wbH+WqHyKtPigzlf2mqkIwKbS922hWJ95L56vj008r2xx+vw803n0qw\n7fvfFHnwtEjKByKDauiukcifDAttiw/2rcCuvDIM48bFy4oVh6RXr61BMTHXdASe3+O9aGfEfbrr\nI/IHw0IbfLCvAWxb0LVrDSxcGGq/9146rr3Wvkak+0e2/c0JkX+8IdInUneNRCVlWGjzhlGmadgw\nAm+/XV2WLv1JHnlkb3h4eMungEHHbHv2f0Ve4HMuyXVMC23emtVQ3qfMx+Hrr2GPH58hV19dORn4\n8zLH+fqQyNgXuWyQ3MK00OZ4hHDTTZXxwQdVrIULf8Rddx2P8XjaDREZeNqypi0QGcRlgxTQDAtt\nnoik/4mJ8WDo0DgsX/6L9eKLm6y4uJq3AQM2OM6CHSKjunPZIAUiw0Ib4HiELiQC3HNPVcydG2V/\n8kkGbr75bIJldX5L5IFs772+n4nVXSNRLsNCmyciqWgJCSEYN66GrFy5X3r23OGJjm78EDBsr+N8\nkiEy4h7d9RGZFtq8IpKKxbYF3brF4quvgu133lmLJk08DUR6zLXtVVki/3hNpFe47hrJTKaFNjge\nocuVmBiJyZNjraVLM+Whh/ZHhIXd9Azw1yzbnrVc5MXrdddHZjEttHkikkosNNTGs8/GYenSHHvc\nuHVSu3Z0S6DfKsdZsl/k9YEiHQ37fSIdDPsh45I/8o9bb43BrFlVrC+//BFt2x6vFhR010sio89Y\n1rtzRYbV1l0fVVymhTY7bfKrKlU8GDEiHsuXZ1kDB260r7jiyruBPtscZ8FWkXFddNdHFY9hoW2B\nM20qC5YluP/+WHzxRZQ9e/YaXH/9mdqW9Ydptp1+SmTK2yLDo3XXSBWDMaEtAvH9P1GZqlMnAm+8\nES+pqT/LH/+4LSQqquFjwFOHbHtOmsjoW3XXR+5mTGgDDG0qXx6PjSefTMCiRR574sQ1Ur9+aDOg\n+9eOk3pIZDwfk0YlYlho5+iugQyVnByN6dOvsP79701o1+5QjMdz21CRgdmW9d7nIn+pq7s+cg/D\nQltxoE1aVa7swZAhCVi+/JT1wgsb7OrVE+4E+m9xnHnbREZ31l0fBT6TQtviSUgKFCLA738fi88+\ni7I//DAdycmqlmV1fde2v/9FZOJbIs9W0l0jBSaTQpvTbApItWuHYcKE6pKaulu6ddsdGhmZ/Dgw\n4Ihtz17NBzXQhUwKbQvIYatNAcvjsfD00zWweLHYEyask6uvrtzC+6CGJQdFXv2LSGtbd42kn0mh\nzU6bXKNly/89qOGOO7KqBAX9boTI26cta+pckb5X6a6P9DEstNlpk7vExHgwfHgcli8/bg0cuNWu\nVq3u3cCQTMf5bLPIsAd110flz6TQNunPShWM94rLapg3L9yeMSMdzZvbdUV6zLDttJPeJ8w/HqG7\nRiofJgUZH+pLFUK9euGYODFWUlN3SqdOe8PCw296Chh53LZnrhB5IVl3fVS2TAptnoikCiU42EKf\nPjXw9de5t4qtegPw3GrHWbRP5NV+fMZlxWRSaPMHmCqsW2+NxqxZ0db8+Rvw29+evMJx7n1V5KEz\nljVllsjzcbrrI/8xLLTZaVPFVq2aB6NHx2HFiiNWv37bnJiY+vcDg3Y7zqcZIqPa666PSs+k0Dbp\nz0qGsyzBI4/E4ssvQ+1p09YiMdFpKPLEZ7a96rjIm2NE/hyiu0YqGZOCjPceISM1bhyJlJRYWbJk\nGx54YF9kSMjNzwL9TlrWzH+LjGisuz66PCaFNu89QkYLD3fQv388UlPP2qNGrbMSEqq0Bp5Z5zj/\n2SHy+hM8cekOJoU2l/wR+bRtWw0ffxxjz52bjptvPpFgWR3+aVl/OCUyJYVP2QlsRoW2CMcjRPnF\nx4dh3LgEWbHigHTvvjU4KqpBV+DpQ7b98TciL/NmVQHIpNDmeISoEI5j4fHH47FoUbA9ceIPUq9e\ncDPg8WW2veyAyIQBIh1NyoqAZtJ/CJ6IJCqG5ORovP9+rP3VV5tw552HqjrObaNERp+2rPc+EhkZ\nr7s+05kU2uy0iS5DdLQHw4dfiRUrTljPPbfeqVo1oQPQc5dtT5iiuzaTmRTa7LSJSsCyBA8+GIf5\n86Ps229fApF6fKK8RqaFtu4aiFzN41EALC4N1Mik0OZ4hKiURHIA2AxtjUwKbY5HiErJshjaupkU\n2uy0iUpJJAdKcTyik0mhzU6bqJQsS4Gdtl6mhbbuGlxv9uzZaNWqFdatW6e7FNLAO9Nmp62TSaFt\nsdMuvbvvvhvBwcFo3Jg3hzORt9NmaOtkUmiz0/aDtLQ0NGvWDCL8vTWRb/WISbkRcEz65jO0/eC/\n//0vRAQLFizAqFGjsGXLFt0lUTlip62fSaHN8chlmjlzJtq0aYPOnTvjp59+AgCsXr0anTp1wp13\n3olbb70Vb775puYqqTxxpq2fSaHNTvsypKWl4bXXXsPrr7+OkydPYsSIEdi7dy+UUmjSpAkA4ODB\ngzh69KjmSqk8cZ22fiaFtsX7aRff+PHjceONN6J+/fpQSiE2NhYbNmxA06ZN87ZZtWoVbrrpJo1V\nUnljp62fSaHNTruY1q1bhw0bNuD2229HcHAwPvnkE7z00ksIDw9HZGQkAGDHjh3YsmULOnfurLla\nKk+WpXhxjWamhTZTuxi++OILALioi05OToZlWfj8888xY8YMTJw4ESEhfKi3SXjvEf0c3QWUI17G\nXkxLly5FnTp1EB19/qMCRQR9+vQBALRv315HaaQZV4/oZ1qnrbuGgLdjxw7s37//vNk1US7fTNuk\n3Ag4Jn3zueSvGL755hsAQGJiouZKKBCx09bPpNBmp10MaWlpAICGDRtqroQCEVeP6GdYaOcwtS/h\nu+++g8fjQe3atXWXQgHIu06boa2TSaFt0p+1RH766SccPnwYdevWhW3busuhAMTVI/qZFGQcj1zC\n999/DwC45pprNFdCgYozbf1MCm2L45GifffddwCAq6++WnMlFKj45Br9TApt/qBdQnp6OoDACO1z\n585peS8VjZ22foaFNjvtwhw5cgS7du2CiKBu3bpaa1m0aBHmzZtX4vdPnTo17y8g8i+u09bPpG8+\nr4gswtq1awEA0dHRqFy5cpkfb+fOnejbty/Gjx+P0aNHQ/mW0H/77bf4/vvvcc8995R43926dUNK\nSgoyMzOLtf1zzz2HTp068SrPYuCSP/1MCm3+oBUhN7TLYzTy66+/olevXmjTpg0OHTqEuXPn4uTJ\nkzhx4gTGjx+PXr16lWr/juNg4MCBGDJkSLFGJS+//DKaN2+Offv2leq4JuB4RD+TQtsS4XikMLnj\nhHr16pX5sVauXIk9e/agefPm6NixI8aPH4+IiAhMnToVd911F4KDg0t9jOrVq6NOnTr4/PPPL7mt\nbduoX79+qY9pAq7T1s+k0OaSv0KcO3cOGzZsAFA+nfa3336L6OhoxMfHo3Hjxrj++utx6tQpfPLJ\nJ2jXrp3fjvPggw9i2rRpftsfcTwSCEwLbaZ2ATIzM5GdnQ0RKZeOMyMjA40aNTrvtdTUVMTFxSEq\nKspvx6lfvz6OHj2KjRs3+m2fpuN4RD/empWwfv16AN4xQZ06dcrsOEOHDsXhw4exZs0a1KpVC717\n90Z8fDwGDBiAVatW4dprry30vRs2bMC8efNgWRb27NmDwYMHY86cOcjKysL+/fvRo0cPJCQknPce\ny7KQlJSElStXnnfB0NatWzF16lRERUUhODgYHo/notvQlvbYFRU7bf1MCm2ORwqRG9q1a9eG45Td\nj8TQoUOxe/dudOjQAT179kTr1q3zvrZp0ybce++9Bb5v165d+PTTTzFgwIC8/XTr1g3Dhg1DTk4O\nnnjiCTRo0ACdOnW66L1XXXUVNm/enPf5mjVr0KdPH4wbNw7NmzcHAPzyyy94+umnIXJxFpXm2BUR\nO239TBqP8NashcgN7QYNGpT5sXJHFReOYfbs2ZP3KLMLvf/+++jdu3fe56dOnUKlSpXQpEkTVK9e\nHZ06dSp0iWBkZCT27NkDAMjJycHQoUORnJycF9gAEBYWhrZt2+YtO/TXsSui3E5bhKuxdDEptNlp\nF+DcuXPYsmULgPK5HeumTZsQERGBuLi4814/ceIEIiIiCnxPly5dEBoamvd5eno6rr/+egBAbGws\n+vTpU+gsvFKlSsjKygLgXda4a9cuJCUlFbve0hy7IvL+YyQHMCs7AopJ33ieiCxAZmYmzpw5AxEp\nl9DeuHFjgSc7RaTAThfAeQGfmZmJAwcO4LrrrivW8ZRSefvNXYd9OSFbmmNXXDkKAG8DqYlJoW35\nOgTKZ9OmTQC8F6SUx939Nm/eXGBoR0ZG4tixY5d8f1paGoKCgs47abl79+5Ctz927Fje2CU2NhYA\nkJ2dfblll+jYFVdODhja2pgU2gLORy6SG9pXX301goKCyvRYx44dw759+woM7bi4uAJDOzs7G3//\n+9/zRjirVq1CvXr18i7AycnJwfTp0ws95vHjxxEfHw8AeXPodevWXbRdQVdOlvbYFVeOglmLGAKK\nSaHNTrsAW7duBQA0bty4zI+VexKyoKsumzZtiu3bt1/0+vLlyzF9+nRs27YNmZmZ2LlzJzweT97X\np06dWuQ9Q7Zv35439rFtGy+++CJSU1PPW1Fy8OBBfPbZZwDO75xLe+yKi+MRnUz625InIguQ20WW\nR2j/+OOPiIyMLLDTvvHGGzF27NiLXm/RogXat2+PDRs2YOPGjXjnnXfwt7/9DaNGjYLH40Hr1q0L\nrf3s2bNYu3bteas/kpOTMX78eLz99tuoUaMGwsLC4DgO2rVrh3fffRd9+/ZFp06d0KFDh1Idu2Jj\naOtkWmgztfPJysrCgQMHICLl8vT1H3/8EcnJybAKuLNns2bNcOjQIRw4cADVqlXLe71y5coYMmTI\nedsOHTq0WMfLyMhAbGzsRZ19YmIixowZc9H23bp1O+/z0hy7YmNo62TYeISZnV/uaCQqKgq1atUq\nk2O888476NmzJwDvevA2bdoUuJ3H40HHjh0xY8YMvx37gw8+QOfOnf22P8rF0NbJpNBmp32Bbdu2\nAfDOk8vK/Pnz4fF4sHnzZgQFBRUa2gDQtWtXrFixAsePHy/1cTMzM7Fnzx6jLnwpPwxtnUwKbct7\nNRflyg3tZs2aldkxunTpgqpVq2LKlCl45ZVXinzKe0hICAYPHoyRI0cWuma7OE6fPo0xY8Zg5MiR\nBV6aTqXF0NbJsJk25yP55S73K8tOu3379pe1wqJx48a47777MGvWLDz00EMlOubUqVPRs2dPXHXV\nVSV6PxVNJCdHKYa2LoaFNjM7vy1btiA0NLRc7jlyOVq2bImWLVuW+P1PPvmkH6uhi3Gdtk5GjUe4\nTvt/fv75Z2RlZSExMbHIkQXRxTge0cmk0OaJyHxyn1TTokULzZWQ+zC0dTIptLnkL5+MjAwAyLtj\nHVHxMbR1Mim02Wnnk56ejvDwcEOv6KPSYWjrZFpo664hIGRnZyMjIwM33HBDgVcnEhXtHENbI5N+\nYy1fh2C8tLQ0nDlzBq1atdJdCrkSO22dTAptYzvtV199FQ8//DDOnj0LAFiwYAGioqKKvDqRqHDn\neD9tjUwKbWNPRK5evRrZ2dnIycnB3r17sWjRIjz88MN594Umuhwi7LR1MmmBvLEnIpOSkhATE4Pj\nx49j+PDhqFmz5kV3tCMqPl5co5NJnbax45GePXsiIyMDHTp0gMfjwRtvvAHHKfh37uzZs5g4cSI+\n+ugjzJw5E/369cOuXbvKuWIKbOy0dTLpb0vL988641SuXBkTJkwo1rajR49GvXr18MADD+Do0aOY\nNGkSEhISyrhCcheGtk7stCnP5s2b8dVXX+G+++4D4L03SfPmzTVXRYGHS/50Mim0LVNn2sW1atUq\nNG3aNO85iKtXr8Z1112HrKwszZVRYGGnrZNJoc1O+xIqVaqEKlWqAAB++eUXLF68GElJSZg/f77m\nyiiwMLRyFtzuAAAChUlEQVR1MmmmzdC+hDvuuAM//PADvvzyS5w5cwZ33HEHli9fHnC3biXduE5b\nJ5NCm1dEXoLH48HgwYN1l0EBjzNtnTgeIaLL4luFZVLDF1BMCm2eiCTyC3baOpkU2uy0ifyCJyJ1\nMi60c59ATkQlo9TqEDC0tTEptC0gR23fvl13HUSuptR3DG2NTAptjkeI/EIBDG1t/HoGWAT3+HN/\nftYQUDh48CCWLl2quxYiV9q8eTN8M+2kAP99r7BE+WlBhYiwjSUiKgGllBR3W7+FNhERlT2TZtpE\nRK7H0CYichGGNhEVi4jUEpF03XWYjqFNROQiRoW2iHwsImkisk5EntBdD5ELOSIyXUTWi8iHIhKq\nuyDTGLV6RESilVJHfD9oqwG0Ukod1l0XkRuISC0A2wDcrJRaKSIpANYrpcZqLcwwRnXaAPqIyA8A\nVgJIAFBPcz1EbrNTKbXS9/F0ALfoLMZExtwTV0RaA2gDoKVSKltEFgMI1lsVkevk/6e5gPeGKHcm\nddpRAI74ArsBgJa6CyJyoZoikvu78wiAZTqLMZFJob0A3pMo6wGMhndEQkTFpwBsBNDT93tUCcBE\nvSWZx6gTkUREbmdSp01E5HoMbSIiF2FoExG5CEObiEgjEUkWkTUiEiwi4b4rthsVuj1PRBIR6SUi\nIwCEAAiF9wKmlwvdlqFNRKSXiAQBSANwCsCNqohg5niEiEi/qgDCAUTA220Xip02EZFmIvIpgA8A\n1AFQQyn1TGHbGnPvESKiQCQijwI4rZSaKSIWgBUi0loptaTA7dlpExG5B2faREQuwtAmInIRhjYR\nkYswtImIXIShTUTkIgxtIiIXYWgTEbkIQ5uIyEX+H+O9D5JZdFwUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa756276350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Example 2/17, Hayler(2012)\n",
    "#Suppose that the probability density function of the amount of milk deposited in a milk\n",
    "#container is\n",
    "from matplotlib.patches import Polygon\n",
    "\n",
    "def pdf(X, f, a, b):\n",
    "    \"\"\"\n",
    "    define a pdf based on func\n",
    "    \"\"\"\n",
    "    val = []\n",
    "    for x in X:\n",
    "        if x < a:\n",
    "            val.append(0)\n",
    "        elif x > b:\n",
    "            val.append(0)\n",
    "        else:\n",
    "            val.append(f(x))\n",
    "    return np.array(val)\n",
    "\n",
    "def f(x):\n",
    "    return 40.976 - 16*x - 30*np.exp(-x)\n",
    "\n",
    "#for 1.95 ≤ x ≤ 2.20 and f (x) = 0 elsewhere. This is a valid probability density function\n",
    "#since it is positive within the state space [1.95, 2.20] and\n",
    "#2.20\n",
    "\n",
    "#(40.976 − 16x − 30e −x ) d x = [40.976x − 8x 2 + 30e −x ] 2.20\n",
    "# 1.95\n",
    "# 1.95\n",
    "# = 54.751 − 53.751 = 1\n",
    "#Figure 2.24 illustrates the area that corresponds to the probability that the actual amount\n",
    "# of milk is less than the advertised 2.00 liters. The area can be calculated to be\n",
    "#2.00\n",
    "#(40.976 − 16x − 30e −x ) d x = [40.976x − 8x 2 + 30e −x ] 2.00\n",
    "#1.95\n",
    "#1.95\n",
    "#= 54.012 − 53.751 = 0.261\n",
    "#Consequently, about 26% of the milk containers are underweight.\n",
    "\n",
    "#for 1.95 ≤ x ≤ 2.20 and f (x) = 0 elsewhere.\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "a, b = 1.95, 2.20 # integral area\n",
    "X = np.arange(0, 3, 0.001)\n",
    "y = pdf(X, f, a, b)\n",
    "plt.plot(X, y, linewidth=1)\n",
    "# make the shaded region\n",
    "ix = np.arange(a, b, 0.01)\n",
    "iy = f(ix)\n",
    "verts = [(a,0)] + list(zip(ix,iy)) + [(b,0)]\n",
    "poly = Polygon(verts, facecolor='0.8', edgecolor='k')\n",
    "ax.add_patch(poly)\n",
    "\n",
    "plt.text(0.5 * (a + b), 1,\n",
    "     r\"$\\int_a^b f(x)\\mathrm{d}x$\", horizontalalignment='center',\n",
    "     fontsize=20)\n",
    "\n",
    "plt.axis([1.9,2.3, 0, 8])\n",
    "plt.figtext(0.9, 0.05, 'x')\n",
    "plt.figtext(0.1, 0.9, 'y')\n",
    "ax.set_xticks((a,b))\n",
    "ax.set_xticklabels(('a','b'))\n",
    "ax.set_yticks([])\n",
    "plt.title(r'example $p.d.f.$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative Distribution Function -- Continuous case\n",
    "\n",
    "The $CDF$ of a continuous random variable $X$ is defined in the exact\n",
    "same manner as the discrete case:\n",
    "\n",
    "$$F(x) = P(X \\le x)$$\n",
    "\n",
    "Characteristics of the $CDF$ are that it is a function $F(x)$ that\n",
    "provides a **nondecreasing** response with an initial value of $0$ prior\n",
    "to the state space and reaching $1$ at the end. Values outside the state\n",
    "space $S$ are $0$.\n",
    "\n",
    "As is the case with the $c.d.f.$ of the discrete probability mass\n",
    "function, knowledge of the $CDF$ or the $p.d.f.$ allows the other to be\n",
    "calculated. Whereas the $c.d.f.$ is the sum of probabilities in the\n",
    "discrete case, the $CDF$, is the integral:\n",
    "\n",
    "$$F(x) = P(X \\le x) = \\int_{-\\infty}^x{f(y)dy}$$\n",
    "\n",
    "#### Example 2/15 Hayler, 2012\n",
    "\n",
    "Suppose that the battery failure time, measured in hours, has a\n",
    "**probability density function** given by:\n",
    "\n",
    "$$f(x) = {2 \\over (x+1)^3}$$\n",
    "\n",
    "for $x \\ge 0$ and $f(x) = 0$ for $x \\lt 0$. This is a valid probability\n",
    "density function because it is positive and because:\n",
    "\n",
    "$$\\int_{0}^\\infty {2 \\over (x+1)^3} = \\left[ {-1 \\over (x+1)^2} \\right]_{0}^\\infty = [0] - [-1] = 1$$\n",
    "\n",
    "If we are interested in calculating the 'chances' the battery will last\n",
    "longer than 5 hours, we can look at the probability that it fails\n",
    "*within* 5 hours:\n",
    "\n",
    "$$P(0 \\le X \\le 5) = \\int_{0}^5 {f(x)dx} =  \\int_{0}^5 {2 \\over (x+1)^3} = \\left[ {-1 \\over (x+1)^2} \\right]_{0}^5 = \\left[ -1 \\over 36 \\right] - \\left[ -1 \\over 1 \\right] = {35 \\over 36}$$\n",
    "\n",
    "Thus, if the chances are the battery will fail are $35/36$, then the\n",
    "**complement** is that the batter will last *longer* than 5 hours, and\n",
    "is $1/36$\n",
    "\n",
    "What is the probability the batter will last between 1 and 2 hours?\n",
    "\n",
    "For this we can use the **cumulative distribution function**, or $CDF$:\n",
    "\n",
    "$$F(x) = P(X \\le x) = \\int_{0}^x {f(y)dy} = \\int_{0}^x {2 \\over (y+1)^3}dy = \\left[ {-1 \\over (y+1)^2} \\right]_{0}^x = \\left[ -1 \\over (x + 1)^2 \\right] - [-1] = 1 - {1 \\over (x + 1)^2}; for x \\ge 0$$\n",
    "\n",
    "In the code below we illustrate the calculation of the $p.d.f.$ and the\n",
    "$c.d.f.$ for calculating the likelyhood the battery will lose powere\n",
    "within 1 and 2 hours:\n",
    "\n",
    "$$P(1 \\le X \\le 2) = F(2) - F(1) = \\left( 1 - {1 \\over (2 + 1)^2} \\right) - \\left( 1 - {1 \\over (1 + 1)^2} \\right) = {8 \\over 9} - {3 \\over 4} = {5 \\over 36}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fa755bd59d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHmxJREFUeJzt3Xu4XHV59vHvTUI4n5FTEtlaOYioHCxEMRIpKlAJ0FqB\nWhSllr6v+BativoqYFtFWluBC+RVBAIKxAKVg4IWRJCDghyiIEEBCUKQAE04nwJ53j/Wb8MwmeyZ\nvffM/NZv9v25rnXNWrNO996EZ9Z+1lqzFBGYmVm5VsodwMzMxseF3MyscC7kZmaFcyE3MyucC7mZ\nWeFcyM3MCudCbrUiaTVJF0t6VNL3Wsw/WtJ3cmTrFkmnS1os6RcdLHuJpIPS+MGSru59QiuNC7kB\nIGmBpKclPSHpwVRs1kjzrpT0jKTHJT0m6UZJR0ia0rD+0ZKWpvWHh0+NIcr7gI2A9SNi/xbzx3zj\nQ6tCKGmOpH8e6zbHkGEmsDuwWUTMaLd8ROwVEUV/cFnvuZDbsADeGxFrATsAbwG+0DDvYxGxNrAJ\n8I/AAcAlTeufExFrNQxfG0OOzYHfRcSyFczXGLbZM5ImjXKVzYEFEfFsL/KsyBhyWkFcyG05EfEA\n8CPgDQ1vK817JiKuAmYDb5X05w3zOyqykl6fjvKXSLpN0t7p/S8BXwT2T0f0H24VD1hV0tz0F8JN\nkt7UsO3PSrorzfuNpH2H9wmcnDI/kfb9UeCvgc+k9y5My24m6XxJD0n6vaSPN2z/aEnnSfqOpMeA\nz0p6StL6DcvskNZ9RfGUdAhwSkOGoyStK+kHafnFqa00tWGdK9N6zb/DIUnLJK3Uatn018e1kv5D\n0iPAUZKmSPqapHvTX10nS1q1k/9mVm8u5NZIAJKmA3sCtzTMe0VLIyLuA24EZo5qB9LKwMVUHxSv\nAj4OnCVpy4g4CvgKMDcd0Z++goz7AP8JrAecDVzQUDTvAt6e/nr4EvBdSRtHxHzg74Gfp22vFxGn\nAGcBx6b39kmF8eL0s28G/BlwuKR3N2SYDZwbEesA/w5cCby/Yf5BVH+dvNj0Ozu1KcOXqP4fPBV4\ndRqeAU5sXI3O20nNy+4E3E3VqvoKcCzwOuDN6XUqcGSH27YacyG3YaIqiEuAq6mK01farPMAVTEd\n9v50pLskHV1u0mKdGcAaEfHViHghIn4K/AA4sCFHuyP7GyPiv1Kh/A9gVeCtABFxXkQ8mMb/E7gT\n2Llh2600vv+nwIYR8S8p3z3At6laScOui4iL0j6eBc4E/gZeamEcAKyor/2KDBGxOCK+HxHPRsST\nVL/zXdv8/J16ICJOSm2q54CPAp+MiEfTvo5p+rmsUJNzB7DaCGCfiLhiFOtMA65pmP5eRHywzTqb\nAfc1vXcv1dFhp+4fHomIkHQ/sCmApA8CnwCG0iJrAhuMYtubA5ulD7Rhk4Cftdp/ciFwsqQhYGvg\nsYi4sZOdSVod+DrwHl7+UFxTkmL832jX+Ht+FbA6cJP00meJ8MHcQHAhtzFJ7ZcdqI7qXnq7g1Uf\nAKY3FarNgTtGsfvpDTlWovpAeUDS5sC3gN2o2hch6ZaGXK0KY/N7fwDuiYgtV7Dv5VodEfGspHOp\njsq3pjpC79Q/AlsCO0XEQ5K2A25OmUcq5E+l19WBJ9N4819Ajes/QtW22SYi/jiKfFYAfxpbp4b7\n56tL2pXqKPT6iLhk5NWW8wvgaaoTjCtLmgW8F5g7im3sKGk/SZOBw4Fn03bXoCpejwArpZOl2zas\ntwiYlvr0je+9tmH6BuAJSZ9RdU37JEnbSnpLmr+iD6szgQ9T9c9Hc7ngmlQF9rF0wvSoTlaKiIeB\nhcBBKeNHgD8ZYfllVCdaj5P0KgBJU5t6/1YoF3Lr1ImSHgcepGoFnAvs0TC/o5NyEbEU2JvqZOrD\nVCf2DoqI33W4nQAuAPYHFgMfAP4iIl6MiNupTj7+POXclle2fn4C/AZ4UNJD6b1TgW1SX/+/UsF7\nL7Ad8PuU8VvA2iPli4hrgWXATelE8Ej5G9c/DliN6sPnOuDSEX7+5nU/Cnw6rbsNcO0IywIcQXUy\n+BfpipvLqP4asMJppDZc+vP5TKqz3gF8KyJOaLHcCVT/Yz4NHBwRtzQvYzboJF0OnB0Rp+XOYhNL\nux75UuATETFP0ppUJ0ouS5dyASBpL+B1EbGFpJ2prtVte8ea2SCR9KdU5wz2yZ3FJp4RWysR8WBE\nzEvjTwLzqa46aDQbOCMtcz2wrqSNe5DVrJYknUHVpjg8Ip5qt7xZt3V81Uq6tGp74PqmWVN55WVO\n91NdRbBonNnMihARH8qdwSa2jk52prbKecA/pCPz5RZpmvYTnc3M+qTtEXm6VOt84LsRcUGLRRbS\ncF0v1dH4whbbcXE3MxuDiBjxHo0RC7mqW8BOBW6PiONWsNhFwGHAXEkzgEcjomVbpV2YOpN0dEQc\nnTvHWJScHZw/N+fPq5OD4HZH5LtQ3a3263SHHMDnqb7ch4j4ZkRcImkvSXdR3W3W6hvrBsFQ7gDj\nMJQ7wDgN5Q4wTkO5A4zTUO4A4zSUO0CvjVjII+IaOuijR8RhXUtkZmaj4js7Ozcnd4BxmJM7wDjN\nyR1gnObkDjBOc3IHGKc5uQP02oh3dnZ1R9V3JBXbIzczy6GT2ukj8g6lL3cqUsnZwflzc/7662sh\nl/y1uWZm3dbX1grEJhG+49PMrFN1bK1s1Of9mZkNPBfyDpXcZys5Ozh/bs5ffy7kZmaF63eP/PAI\nju/LDs3MBoB75GZmE0C/C3mxD5wouc9WcnZw/tycv/58RG5mVrh+98ivj/DzPM3MOuUeuZnZBOBC\n3qGS+2wlZwfnz83566/fhXySxBp93qeZ2UDrd4/8D8CuESzoy07NzApXxx75QxTcXjEzqyMX8g6V\n3GcrOTs4f27OX38u5GZmhet3j/xfgSURfLUvOzUzK5x75GZmE4ALeYdK7rOVnB2cPzfnrz8XcjOz\nwvW7R749cEYEb+7LTs3MCuceuZnZBNDvQv4IsKHU9/2OW8l9tpKzg/Pn5vz119eCGsHzwBPAev3c\nr5nZIOtrjzwiJHEHsF8E8/uyYzOzgtWxRw7uk5uZdZULeYdK7rOVnB2cPzfnrz8XcjOzwuXokR9d\n7Zej+rJjM7OCuUduZjYBuJB3qOQ+W8nZwflzc/76cyE3Mytcjh751sCFEWzVlx2bmRXMPXIzswkg\nRyF/FFhTYkqGfY9ZyX22krOD8+fm/PXX90IewTLgYXxUbmbWFX3vkVfjzAM+EsHNfdm5mVmh6toj\nB/fJzcy6xoW8QyX32UrODs6fm/PXX65CvojCCrmZWV3l6pEfAWwQwWf6snMzs0K5R25mNgG4kHeo\n5D5bydnB+XNz/vprW8glnSZpkaRbVzB/lqTHJN2Shi90sN/iCrmZWV217ZFLmgk8CZwZEW9sMX8W\n8MmImN1mO4098s2BqyN49ViDm5lNBF3pkUfE1cCSdvsaTTDSnZ3SqNczM7Mm3eiRB/A2Sb+SdImk\nbdquEDwNLAXW6sL++6LkPlvJ2cH5c3P++pvchW3cDEyPiKcl7QlcAGzZakFJc4AF1dSXn4Fb9oNz\nz0jzZgFExJV1nAa2k1SbPJ72tKcHczqNH0xlAR3o6DpySUPAxa165C2WvQfYMSIWN73/ij6PxH8D\nX4/g0k6CmplNRH25jlzSxpLSSUztRPXhsLjNalB90mw+3v2bmU10nVx+eA5wHbCVpPskfUTSoZIO\nTYu8D7hV0jzgOOCADve9ABgafeQ8Su6zlZwdnD8356+/tj3yiDiwzfyTgJPGsO97gb3HsJ6ZmTXI\n8l0r1TS7AP8ewYy+BDAzK1BfeuTjsAD3yM3Mxi1nIf8jsJ7EahkzdKzkPlvJ2cH5c3P++stWyNOz\nO+8H36ZvZjYe2Xrk1XtcDvxbBD/uSwgzs8LUvUcOhV2CaGZWR3Uo5EWc8Cy5z1ZydnD+3Jy//upQ\nyIcyZzAzK1ruHvlM4NgI3taXEGZmhXGP3MxsAshdyB8ANpBYNXOOtkrus5WcHZw/N+evv6yFPIIX\nqa4ln54zh5lZybL2yKv3uQI4JoLL+hLEzKwgJfTIwX1yM7NxcSHvUMl9tpKzg/Pn5vz150JuZla4\nOvTI3wF8JYK39yWImVlBSumR34uPyM3MxqwOhXwhsJHEKrmDjKTkPlvJ2cH5c3P++steyCN4gaqY\n+1pyM7MxyN4jr+bxU+BfIvhJX8KYmRWilB45uE9uZjZmdSnkC6h5IS+5z1ZydnD+3Jy//lzIzcwK\nV5ce+SzgnyOY2ZcwZmaFKKlHvoBCHvlmZlY3dSnkC4FNJKbkDrIiJffZSs4Ozp+b89dfLQp5BEuB\nPwLTcmcxMytNLXrk1XyuAr4UwRV9CWRmVoCSeuTgK1fMzMakToX8Xmp8wrPkPlvJ2cH5c3P++qtT\nIV+Aj8jNzEatTj3y3YCjIti1L4HMzApQYo/8tblDmJmVpk6F/F5gfYm1cwdppeQ+W8nZwflzc/76\nq00hj+BFYD7whtxZzMxKUpseebUMc4BrIzilL6HMzGqutB45wK3AG3OHMDMrSd0K+W3AtrlDtFJy\nn63k7OD8uTl//dWtkN9KTQu5mVld1a1HLuB/gNdHsKgvwczMaqy4HnkEQY3bK2ZmdVSrQp7U8oRn\nyX22krOD8+fm/PVXx0LuI3Izs1GoVY+8Wo63A1+LYEYfYpmZ1VpH5xdrWMjXA/4ArBPBst4nMzOr\nr+JOdgJEsAR4lJp9N3nJfbaSs4Pz5+b89de2kEs6TdIiSbeOsMwJku6U9CtJ23ch123U8ISnmVkd\ntW2tSJoJPAmcGRHLFVdJewGHRcReknYGjo+I5frbnbZWqmX5N2BJBF/pZHkzs0HVldZKRFwNLBlh\nkdnAGWnZ64F1JW08mqAt+A5PM7MOdaNHPhW4r2H6fmDaOLdZu9ZKyX22krOD8+fm/PU3uUvbaT7s\nb9mvkTSH6klAUJ3QnBcRV6Z5swDS9Hy4YkvpL3ePWHJ5i/l9nwa2k5Rt/572tKcnxnQaP5jKAjrQ\n0eWHkoaAi1fQI/9/wJURMTdN3wHsGhGLmpbruEdeLc9vgb+M4LZO1zEzGzT9uvzwIuCDaYczgEeb\ni/gY+Q5PM7MOdHL54TnAdcBWku6T9BFJh0o6FCAiLgF+L+ku4JvA/+5Stlqd8Cy5z1ZydnD+3Jy/\n/tr2yCPiwA6WOaw7cV7hNuCgHmzXzGyg1O4W/ZeXZyvgkgj+pIexzMxqrchb9BvcDWwqsWbuIGZm\ndVbbQh7BC8AdwDa5s0DZfbaSs4Pz5+b89VfbQp74yhUzszZq2yOv1uEzwKYRfKJHsczMaq30HjnA\nLcCOuUOYmdVZ3Qv59cAOEqvkDlJyn63k7OD8uTl//dW6kEfwOPBbfFRuZrZCte6RV+txPPBABMf2\nIJaZWa0NQo8c4Grg7blDmJnVVQmF/BpgFylv1pL7bCVnB+fPzfnrr/aFPIIHgf+hJjcGmZnVTe17\n5NW6nA7cEMHJXY5lZlZrg9Ijh6q94j65mVkLpRTy7Cc8S+6zlZwdnD8356+/Ugr5ncCqEq/OHcTM\nrG6K6JFX63M+cH4EZ3cxlplZrQ1SjxyqPvnM3CHMzOqmpEKetU9ecp+t5Ozg/Lk5f/2VVMjnAZtL\nrJ87iJlZnRTTI6+2wWXA8RH8oEuxzMxqbdB65ODryc3MllNiIc9ywrPkPlvJ2cH5c3P++iutkP8C\n2E5itdxBzMzqoqgeebUdbgA+FcHPuhDLzKzWBrFHDnAV8M7cIczM6qLEQn4xsE+/d1pyn63k7OD8\nuTl//ZVYyK8FpkkM5Q5iZlYHxfXIq21xKvDrCI7vxvbMzOpqUHvkABcA++UOYWZWB6UW8suB7SU2\n7NcOS+6zlZwdnD8356+/Igt5BM8AlwF7585iZpZbkT3yant8ANg/gtnd2qaZWd10UjtLLuTrAfcC\nm0XwZLe2a2ZWJ4N8spMIllDdsv+efuyv5D5bydnB+XNz/vortpAnFwD75g5hZpZTsa2VaptMBW4F\nNo5gaTe3bWZWBwPdWgGIYCHwO2DX3FnMzHIpupAnfbk5qOQ+W8nZwflzc/76G4RC/n1gX2kgfhYz\ns1Erukf+8raZDxwSwXW92L6ZWS4D3yNvcDrwd7lDmJnlMEiFfF+J9Xu1g5L7bCVnB+fPzfnrbyAK\neQQPAz8ADs4cxcys7waiR15tn12ojsy3jmBZr/ZjZtZPE6lHDnAd8CywW+4gZmb9NDCFPIIAvgH8\nr15sv+Q+W8nZwflzc/76a1vIJe0h6Q5Jd0o6osX8WZIek3RLGr7Qm6gdOQt4Z7p138xsQhixRy5p\nEvBbYHdgIfBL4MCImN+wzCzgkxEx4veC97pH/vJ++AbwUARH93pfZma91o0e+U7AXRGxICKWAnOB\nfVrta4wZe+Fk4KMSK+cOYmbWD+0K+VTgvobp+9N7jQJ4m6RfSbpE0jbdDDhaEdwK3EOXHwNXcp+t\n5Ozg/Lk5f/1NbjO/k2sTbwamR8TTkvak+hKrLVstKGkOsCBNPgrMi4gr07xZAF2aPhku+Ly03+Ju\nbR/YTlK38nna0572dMvpNH4wlQV0oF2PfAZwdETskaY/ByyLiGNHWOceYMeIWNz0fl965NW+WIV0\nVB7BTf3Yp5lZL3SjR34jsIWkIUlTgP2Bi5p2srEkpfGdqD4cFi+/qf6J4Dngy2kwMxtoIxbyiHgB\nOAz4MXA78L2ImC/pUEmHpsXeB9wqaR5wHHBALwOPwinAllJ3HjpRcp+t5Ozg/Lk5f/2165ETEZcC\nlza9982G8ZOAk7ofbXwieF7iKOAYiV3SDUNmZgNnYL5rpfU+mQT8CvhcBBf3c99mZt3QjR550SJ4\nEfi/wJf9BCEzG1QTobhdBDwFHDiejZTcZys5Ozh/bs5ffwNfyFNv/HPAP0lMyZ3HzKzbBrpH/sr9\n82PggghOzpXBzGy0OqmdE6mQ7wD8EHhTeqKQmVntTfiTnY0iuBn4LnDiWNYvuc9WcnZw/tycv/4m\nTCFPjgTeLPG+3EHMzLplwrRWXs7BDOD7uMViZgVwa6WFCH4BfIca3o1qZjYWE66QJ0cCb5R4f6cr\nlNxnKzk7OH9uzl9/E7KQR/As1ff9niCxUeY4ZmbjMuF65I0kvgq8iep7y1/MncfMrJl75O19EVgF\n+NfcQczMxmpCF/IIlgJ/BewtcchIy5bcZys5Ozh/bs5ff22/j3zQRbBYYm/gZxJ3RXBV7kxmZqMx\noXvkjSTeRXVZ4tsi+H3uPGZm4B75qERwGfBPwMUS6+TOY2bWKRfyBhF8A7gc+JHEuo3zSu6zlZwd\nnD83568/F/LlHQ7cAPxEYoPcYczM2nGPvAUJAccC7wHeFcFDmSOZ2QTlHvkYpacKHQFcCFwpsWnm\nSGZmK+RCvgIRRARHAmcBV0l//te5M41V6T1C58/L+etvwl9H3k4EX5ZYDIefKLE4gh/lzmRm1sg9\n8g5JzATmAt8AjolgWeZIZjYB+JmdXSYxFTgPeBD4UASPZ45kZgPOJzu7SNKsCBYCs4BFwA0SO+dN\n1ZnSe4TOn5fz158L+ShF8FwEfw8cBVwo8TWJ1XPnMrOJy62VcZB4FXAC8BbgkAh+ljmSmQ0Y98j7\nRGJfqmeAXgwc6RuIzKxb3CPvopH6bBFcAGwLPAPcLvEFiTX6la2d0nuEzp+X89efC3mXRLAkgk8A\nOwFvAH4n8beSr9U3s95ya6VHJHaieoTcNOBrwBkRPJM3lZmVxq2VjCK4AXgncDCwF3CPxBf9jYpm\n1m0u5B0aS58tfV/LNRHMpirqQ8CdEqdIzEjfsthzpfcInT8v568/F/I+iWB+BIcA2wB3UT1W7laJ\nT0hsmDedmZXMPfJM0tH4O4BDgNnAtcC5wIURLMmZzczqw9eRF0JibeC9wF8BfwZcA5wPXBrBAzmz\nmVlePtnZRb3ss0XweARnR7AfMBU4k+rpRLdJzJM4RuIdEiuPZful9widPy/nrz9f41wzETxB9XW5\nc9M16DsDewJfB7aQ+DlwFfAz4JcRPJctrJnVglsrBUmXLs6k6q3vCmwF3Ez1sOjr0+sf0qPqzGwA\nuEc+4CTWobqTdHjYGRBVcf9VGuYBv4vgxVw5zWzsXMi7qPo+8rgyd46RpCthpgPbpeHN1XDFNNjt\nDmA+cHt6nQ/cHcGzufJ2qoTf/UicP68ByN+2drpHPkBSS+UPabho+H3pI3vCgoeprmF/PfA3wNbA\nkMRDVNe13wncDSxIwz3AI27TmNWfj8gnsHQydTqwRRpeQ3X36fDrqlQfCvcD9zW8/hF4IL0+5LaN\nWe+4tWLjkq5vn0ZV7BtfNwM2Ta/rA49QPf5uEfBQw+sjLYbHXPjNOteVQi5pD+A4YBLw7Yg4tsUy\nJ1BdIvc0cHBE3DKWMHVWcp+tl9nTte0bARu3eN0A2DANr6Iq+msDjwOL07AEeLTh9VHgsTQ8Xr1+\nYEs46wrgCeDxCJ7vxc/SKyX/2wHnz23cPXJJk4ATgd2BhcAvJV0UEfMbltkLeF1EbCFpZ+BkYMa4\n09fPdsCVuUOMUc+yR7CU6t/Gwk6Wl5gErENV5NcH1gPWbRjWo2rrrENV9NeBrV8LfBFYC1hbYhlV\nUX+yYXgCeIrqYOKppuHphuGZhtfm4dnh1y7/1VDyvx1w/tprd7JzJ+CuiFgAIGkusA/VFQ/DZgNn\nAETE9ZLWlbRxRCzqQd6c1s0dYBxqkz0VyOGj8Y5IRx4d8cWjq3EErEJV5NdIw5pURX6NFsOaVH8h\nrJ6GNYDVmobV0zZXozovsJrEi1SF/VnguYbXxuH5FuPPNw3PwV/MTH+5LE3D8w3jrYYXWky/0DTe\ncohgWae/11Gozb+fMSo9f1vtCvlUqpNbw+6nula53TLTqPqkZl2VrqIZLrA9kT4sVqYq7qs2vTYO\nqwJT0rBK02vDsPIUqg+TKWm7w68rGia3GJ/cMD6pYbrxvZVV/QE+XNhfbPHaPN5q+kVg2cvjH9tc\n4q2t57003vi6rMV087zm96LF/GUrmN9uvOl1tzdKHLDi+USL8U6mm4d285sHWry3MILHGaV2hbzT\nM6HN/ZtBvGRtKHeAcRjKHWCchvq5s/RhMXxE/cR4tyd9byhi7qfHHayjfbES1f/Xk5peh8cnNc1r\nHp9E9R1MDdM/+gzVV0Q0z1upxetKTcuoaV7ze8PTw/Mnt3hveFDT+PByzcs2vbfy9lR/zTRuo9Xr\nSOOtplsNnSwzXC9bvf9J4JLm/67ttCvkC6muVBg2neqIe6RlprGCfqmkogu8pA/lzjBWJWcH589P\n786dYHz0mtwJOvRDjeGSkHaF/EZgC0lDVNcN7w8c2LTMRcBhwFxJM4BHW/XHS75ixcyszkYs5BHx\ngqTDgB9T/RlzakTMl3Romv/NiLhE0l6S7qK6QuDDPU9tZmYv6dsNQWZm1hs9f7CEpD0k3SHpTklH\n9Hp/3STpNEmLJN2aO8tYSJou6aeSfiPpNkn/J3em0ZC0qqTrJc2TdLukY3JnGi1JkyTdIuni3FnG\nQtICSb9OP8MNufOMRroU+jxJ89O/n2Lub5G0VfqdDw+PjfT/b0+PyNMNRb+l4YYi4MDGG4rqTNJM\nqptNzoyIN+bOM1qSNgE2iYh5ktYEbgL2LeX3DyBp9Yh4WtJkqkfgfSoirsmdq1OSPgnsCKwVEbNz\n5xktSfcAO0ZEx9f914WkM4CrIuK09O9njYh4LHeu0ZK0ElX93Cki7mu1TK+PyF+6oSgillI9+Waf\nHu+zayLiaij3QcgR8WBEzEvjT1LdyLVZ3lSjExFPp9EpVOdpiikokqYBewHfZvlLdEtSXHZJ6wAz\nI+I0qM73lVjEk92Bu1dUxKH3hbzVzUJTe7xPayFdebQ91ZOEiiFpJUnzqG4w+2lE3J470yh8Hfg0\n9ORuy34J4HJJN0r6aO4wo/Aa4GFJp0u6WdIpklbPHWqMDgDOHmmBXhdyn0mtgdRWOQ/4h3RkXoyI\nWBYR21Hdn/COUh6kK+m9wEPpC+SKO6JtsEtEbE/1pXgfS+3GEkwGdgC+ERE7UF1R99m8kUZP0hRg\nb+DckZbrdSHv5IYi6yFJKwPnA9+NiAty5xmr9GfxD4G35M7SobcBs1OP+RxgN0lnZs40ahHxx/T6\nMPB9qnZpCe4H7o+IX6bp86gKe2n2BG5Kv/8V6nUhf+mGovTJsj8NT66x3pIk4FTg9og4Lnee0ZK0\noaR10/hqwLuA5b4iuY4i4vMRMT0iXkP1p/EVEfHB3LlGQ9LqktZK42sA7waKuIIrIh4E7pO0ZXpr\nd+A3GSON1YFUBwIj6umj3lZ0Q1Ev99lNks6helr9BpLuA46MiNMzxxqNXage6/ZrScMF8HMR8aOM\nmUZjU+CMdNZ+JeA7EfGTzJnGqsQ248bA96vjASYDZ0XEf+eNNCofB85KB5F3U9jNiunDc3eg7bkJ\n3xBkZla4nt8QZGZmveVCbmZWOBdyM7PCuZCbmRXOhdzMrHAu5GZmhXMhNzMrnAu5mVnh/j9kuj+H\n0CH2PgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa755d5f890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEZCAYAAACKF66QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X28pXO9//HX2wy5GQymTBhGUaFyF+mEpiiDg5NEIzTp\noBvV+XWHzikqcfTrRtIRcnsqzkkSJVLITe4ZxBAxMsbNDDNuQmbM5/xxXXss2957XWuvtdf3+q79\nfj4e12Ota61rXdd7b+Ozv+tz3SkiMDOzfC2VOoCZmbXHhdzMLHMu5GZmmXMhNzPLnAu5mVnmXMjN\nzDLnQm61JOlISXMlzRngvSmSHkyRq1MkfULSo5KekrRKk2UPk3Ry+XyypMWS/P+uLeF/DLaEpL0l\n3SjpaUlzJF0o6Z3le0dIWlgWnqck3S3pB5ImNnx+Sllknm6YfjWMHGsDnwPeFBFrdO4nHLgQSpou\n6cpObqdJhqWB7wDbRcRKETF/qOUj4uiIOKA76SxHLuQGgKTPAd8DjgReA0wCfgjs0rDYWRGxErAK\n8H5gInBTYzEHHoqIFRum3YYRZ23g8Yh4fDg/S0Xq2IqksS1+ZCKwLDCzUxmqGEZOy4QLuSFpZeBr\nwCcj4ryIeC4iXoyI30TEoY2LApTv3QnsBcwFPj+cbUo6U9JjkmZJ+ncVtgd+B6xRjuhPHWIdh5Xt\nl/sl7d3w+s6SbpH0pKS/STq84WNXlI8Lym8WWwE/At5Rbu+Jch2vkvRtSQ9IekTSCZKWLd+bImm2\npC9Jehg4VdLtkv65IcPSkuZJ2rhf5jfwUgFfIOn35evfL7M+WX4r2rrhM0dI+u9BfgezJG030LIN\n3z72l/QA0Let/SXdKekJSReV34AsYy7kBvAOihHiL1v5UEQsBn4FbDOMbf4AWBFYF3gXsB/w0Yj4\nPbAjMKcc0e8/yOcnAqsBawAfAU4qiyTAM8A+EbEysDPwCUl93wz6sq5ctjWuBQ4Crim3t2r5/n8C\n6wEbl49rAl9t2P7qFN9M1gYOBM4E9ml4fyeKbye3NoaOiL8AGzVk2L58fn25rVWAnwE/l7RM38cG\n+R30vRf95vvbFngTMLX8PRxG8Y1qAnAlcNYQ67cMuJAbFAVxXlmYW/UwsGrD/BqS5jdMe/T/gKQx\nFKP5wyLi7xHxAEXPeN++RSpu+ysRsTAirgB+A+wJEBF/jIg7yue3A2dT/LEYbN0ve02SgAOAz0XE\ngoh4Bjga+FDDYouBw8vtPw/8FNhZ0rjy/X2BAUfRA2WIiJ9GxPyIWBwR3wVeBbxxiMyDGWjZI8pv\nWc8DHweOjoi7y//eRwObSJrUwjasZlzIDeBxYMIwj4RYs/x8nzkRsUrDdM4An5kALA080PDa38p1\nVTU/Ip5rmH+AYnSOpLdLuqxs2yygGHGv1sK6Xw0sT9H/ny9pPvDbMnefuRHxQt9MRMwBrgb2kDQe\nmEpR3CuR9IWy3bGg3N7K/bbXjsYjfNYBvt/wc/X9t2vld28140JuANcA/6D4uj2YV3xlLwv/LhRf\nz1sxD1gITG54bW1gdgvrWEXS8g3z6wAPlc9/BpwHrBUR4yl64H3/1gdqPfR/bR7wHLBhwx+k8eWO\n3sE+A3AGRXvlg8CfIuLhKj+IpG2ALwIfLLezCvAk1UbifwdWaJifOMAyjVn/BhzY74/tCmWLyTLl\nQm5ExJMU/d8fStpN0vLlzrodJR1TLrakqEgaK2kDit7qa4Dvtri9F4H/Bb4paZykdYD/B/ykxehf\nK3NuQ9EL/3n5+jiKEfsLkrYE9ualYjaXoi3y+ob1PAKsVR4W2Nf7Pxk4VtKry595TUnva5Lnl8Bm\nwGcoeuZVrQgsAuZJWkbSV4GVmnymzwzgQ+V/k7cBH2DonvqPgC9L2hCW7HT+YAtZrYZcyA2Asi/7\nOeA/gMcoRm6f5KUdoAHsJelpYAHFTs65wOYR8Ujjqipu8tMUo8n7KEb0PwVOq7ieoOjNzwfmUPSi\nDyp3JFLm/rqkp4CvAP/T8HM+C3wTuLpsL2wJXArcATwi6bFy0UOAe4FrJT0JXAL07UwdMF/Zgz6X\n4pvGuU1+/sbPX1ROfwFmUXwb+Fu/ZQfbofkVij9K84EjeGU752U5I+I84Bjg7PLnuh3YoUlWqzk1\nu7FEefjXzsBjEfGWQZY5juJIg2eB6RFxS6eDmuVA0leA9SNiv9RZbPSoMiI/jWLHzYAk7QSsFxHr\nUxyGdUKHspllRdKqwP7ASamz2OjStJBHxJUUX9sGsyvFTh4i4jpgvKTVOxPPLA+SDqBoh/w2Iq5K\nncdGl070yNfk5Yc3zQbW6sB6zbIRESdHxLiI+GTqLDb6dGpnZ//DpHxHZzOzLunERXQeorjAUp+1\neOl43iUkubibmQ1DRAx5TkEnCvn5wMEUhzNtBSyIiEeHE6bOJB0REUekzjFcOefPOTs4f2o9kL/p\nILhpIZd0FsV1KiaouJj/4RSnVxMRJ0bEhZJ2knQvxXHBH20vdm1NTh2gTZNTB2jD5NQB2jQ5dYA2\nTU4doE2Tu7ERCQFjKOrjMuVj4/PG1/qmsf3m74zg7la33bSQR8S0Cssc3OqGzcxaVRbLsRRX61yW\n4uJifY/9n5fTW14n8a8Nry0zxGP/51WmxiIdFJefeKHfY+PU99qiAd47A0agkNsSp6cO0KbTUwdo\nw+mpA7Tp9NQB2nT6YG9IjKG4wNhA03INj8sNMr8cRfHt/zjUtBh4nuL6QP/o93yA+S//g+JSzX2v\nvVA+PktxlnLjay/0m/q/trDfay8r1hG8WPm32kFNz+zs2IakyLlHbpYziWUorukyruFxqGmFhmkc\nRQFufK2vWC9DURAHmp5reKwyPT/AY9/zJQU6gkWd/v3UWZXa6RF5RZKmRMTlqXMMV875c84O7eWX\nGEtxSduVgfHltFLDays1mVYsH8cAT5fTMwM89k1/p7jWTt/zv8OXXgffuval+SXTsxSFtdZHpOX+\n76cKF3KzEbe0JMZTXBN9NYobcfRNqzQ8H1/ONz4uBzxFcVnbBeVj/2ke8NdyuacoinPfY9/0/HAL\nrvT/p0R864rmS1oqbq2YtUjiVRSX7+2bJlDcjKJvWq18bUL5fBWKEewTFDdyeKLfNL+cnqAo1vPL\nxwXA03Uf8drIcmvFrKJyp91rgNeW08QBptXLZZanaD/MbXjsm+4vHx+nGCnPA+ZHsLCLP46NMh6R\nV5R7ny3n/O1ml1iW4uzjSRRnHq9FcY2gxunVFCPihymucf5I+fzR8vmjDdOCVkbJOf/uwflT84jc\nRgWJlYF1KU78mExx27e+aRLFTsGHKC7u9iDFhd3uAv5Qvj4HeMSjZsuVR+RWexJLUbQ71muYXg+8\njqKAv4qipXE/xU2YG6e/AY9FsLj7yc3aV6V2upBbbUisCGwAvLFhegOwPsXRGfeW01/L6b5ymucd\ngtarXMg7qAf6bLXJL7ECsBHwVmDDhmkCxenJd5WP5fS61SPu+22iuG2r0+9+OJw/LffILanyuhhr\nAZuW08bAWyh2Lt5NcePfP1Pc/PhO4IGBTnGW7p/SpchmWfKI3DpGYk1gS2AL4G0UxXsxcAtwM3Ar\ncBtwz2g7zdpsuNxasREjsRywOfBP5bQlxdXfrgduAG6kKN4Pu39tNnwu5B3UA322do/FHg9sA0wB\n3knRIrkD+BNwDXAdRWuk4/+gRvvvPjXnT8s9chu28giSKcC7y8f1gWuBPwKHADdE8GyqfGb2Eo/I\nDVhyivrmwPuA9wKbUYyyLwUuB26M4IVkAc1GKbdWbEgSK1EU7l2AnShOP78YuAS4wiNus/Sq1M6l\nuhUmd5KmpM7Qjr78EhMlPiFxCcWp6h+j2EG5RQRvjuDzEVxUpyLeK7/7XDl//blHPgpIrAHfe7/E\n1yiO5f4NcALw/gieSZvOzNrl1kqPKndW7g7sS9HvvgA4B7gkgudTZjOz6nzUyihTXlxqO2A6sDNw\nBXAicIGLt1nvco+8ojr32SRWlzgUuAf4FsVhgutHsGsEP4/g+Trnbybn7OD8qeWevwqPyDNVXsdk\na+BgYAfgXGAaxfHdPpPSbBRxjzwz5V3V9wA+T3Fz3h8AZ0awIGkwMxsR7pH3EIlxwIHAZylumHAk\n8OuBrhZoZqOLe+QVpeqzSYyTOITiRgpbAXtEsG0Ev2qliOfcJ8w5Ozh/arnnr8Ij8poqb77wKYoW\nymXAuyO4M20qM6sj98hrpuyBHwAcTnGBqq9HcEfaVGaWinvkmZGYCnyH4ponUyOYkTiSmWXAPfKK\nRrLPJrGRxEXA94FDge06XcRz7hPmnB2cP7Xc81fhQp6QxHISR1NcJvZC4C0RXODjwM2sFe6RJyKx\nHcXp8zcBn43gkcSRzKyG3COvIYlVKfrg2wGfiuCCxJHMLHNurVTUiT5bOQq/DXga2KibRTznPmHO\n2cH5U8s9fxUekXeBxDLAN4B9gOkRXJI4kpn1kKY9cklTgWOBMcCPI+KYfu9PAH4CTKT4w/DtiDh9\ngPWMyh65xBuAnwEPA/tHMDdxJDPLSNu3epM0BjgemApsCEyTtEG/xQ4GbomITSjutv4dSR7pAxLT\ngKuB04BdXcTNbCQ065FvCdwbEbMiYiFwNrBbv2UeBlYqn68EPB4RizobM71W+mwSS0kcBRwFbB/B\nD1MfUphznzDn7OD8qeWev4pmI+c1gQcb5mcDb++3zMnApZLmACsCe3YuXn7KO9P/BFgZ2NKjcDMb\nac0KeZVR5JeBGRExRdLrgUskbRwRT/dfUNLpwKxydkH5ucvL96YA1HW+77Whl99pDfjNvwNXwirH\nwYKNIKf89ZyPiMvrlMf565Wv1/KXz6dTmEUFQ+7slLQVcERETC3nDwMWN+7wlHQh8M2IuLqc/wNw\nSETc2G9dPb2zU+KfKO7S87UITkidx8x6Q9s7O4EbgfUlTZa0DLAXcH6/Ze4Cti83uDrwRuC+4UWu\nr6H6bBLvBc6jOLSwlkU85z5hztnB+VPLPX8VQ7ZWImKRpIOBiykOPzwlImZKOqh8/0SKHXqnSbqV\n4g/DlyLiiRHOXRsS76c41X73CK5KncfMRh9fa6UNEvsBxwA7R3Bz6jxm1nt8rZURJPEp4BDgPRHM\nTJ3HzEYvX2ulosY+m8QBwBeAbXMp4jn3CXPODs6fWu75q/CIvEUSHwC+BrwrotqhQWZmI8k98haU\nVy88C9ghgltS5zGz3uceeQdJbEFxiYI9XMTNrE7cI69A4k1wycXAxyL4Y+o8w5FznzDn7OD8qeWe\nvwoX8iYkVgMuhJtOinjFyVBmZsm5Rz4EiTEUN0W+PYIvpM5jZqNPJ07RH+2OpNiPcGjqIGZmg3Eh\nH4TEHsA0YK8IFuXeZ8s5f87ZwflTyz1/FT5qZQASGwEnAFMjmJc6j5nZUNwj70diPHA9cGQEZ6bO\nY2ajW5Xa6ULej8TZwLwIDk6dxczMOztbJLEnsDHwxVe+l3efLef8OWcH508t9/xVuEdekpgIHEdx\nt/vnUucxM6vKrRVAQhR3+Lk9gv9IncfMrI+vtVLdvsBk4IOJc5iZtWzU98gl1gK+DXwkghcGXy7v\nPlvO+XPODs6fWu75qxjVhbxsqZwCHBfBjNR5zMyGY1T3yCU+Anwa2CqCRanzmJn15+PIhyAxDrgb\n2D2C61LnMTMbiI8jH9phwKVVi3jufbac8+ecHZw/tdzzVzEqj1qRWBc4iOLkHzOzrI3K1orEz4Hb\nIvhG6ixmZkPxceQDkNgW2BLYL3UWM7NOGFU98vKOP8cCh7R6Gn7ufbac8+ecHZw/tdzzVzGqCjkw\nHXgW+J/EOczMOmbU9MjLww3/QnFRrBtT5TAza4UPP3y5g4CrXMTNrNeMikIusSzweeCo4a8j7z5b\nzvlzzg7On1ru+asYFYUc+Chws6+nYma9qOd75BJLA/cA0yK4ptvbNzNrh3vkhb2B+1zEzaxX9XQh\nL48bPwz4ZvvryrvPlnP+nLOD86eWe/4qerqQA7sDC4BLUwcxMxspTXvkkqZSnA05BvhxRBwzwDJT\ngO8BSwPzImLKAMt0tUde3jTiFuArEVzQre2amXVS29dakTQGOB7YHngIuEHS+RExs2GZ8cAPgR0i\nYrakCe1H74idAAG/Th3EzGwkNWutbAncGxGzImIhcDawW79l9gZ+ERGzASJiXudjDsthwFERdOSw\nnNz7bDnnzzk7OH9queevolkhXxN4sGF+dvlao/WBVSVdJulGSft2MuBwSGwMrAP8InUWM7OR1uwy\ntlVGs0sDmwHbAcsD10i6NiLuaTdcGw4ATunkfTgj4vJOrSuFnPPnnB2cP7Xc81fRrJA/BExqmJ9E\nMSpv9CDFDs7ngOckXUFx551XFHJJpwOzytkFwIy+X3Lf15925yGuB6bBdp+QLp3S6fV73vOe9/xI\nzpfPp1OYRQVDHrUiaSzFDYq3A+YA1wPT+u3sfBPFDtEdgFcB1wF7RcSd/dbVlaNWJD4C7BnBzp1d\nr5b8UchRzvlzzg7On1oP5G/vqJWIWCTpYOBiisMPT4mImZIOKt8/MSLuknQRcBuwGDi5fxHvsgOB\nbyXcvplZV/XUtVYk3kzxR2edTvbHzcxSGY3XWjkAONVF3MxGk54p5BLLAfsAp4zM+vM+FjXn/Dln\nB+dPLff8VfRMIQf2AK6PqLaX18ysV/RMj1ziSuC7EfxypLZhZtZto6ZHLrEhsB6+roqZjUI9UciB\njwGnRbBwpDaQe58t5/w5ZwfnTy33/FU0O7Oz9iSWAvYE3pc6i5lZCtn3yCW2ojjkcMNOr9vMLLXR\n0iP/AL7KoZmNYlkX8vIuQHsA54z8tvLus+WcP+fs4Pyp5Z6/iqwLObAp8CLFdV7MzEalrHvkEkcB\nS0VwaCfXa2ZWFz3dIy/bKu6Pm9mol20hBzYClgVu7MbGcu+z5Zw/5+zg/Knlnr+KnAv5HsAvOnVz\nZTOzXGXbI5e4Hfh4BFd3ap1mZnXTsz1yiTcCqwHXpM5iZpZaloWcYifnuREs7tYGc++z5Zw/5+zg\n/Knlnr+KXAt5V04CMjPLQXY9conXAdcCr43gxfaTmZnVV6/2yHcDfuUibmZWyLGQbw9c3O2N5t5n\nyzl/ztnB+VPLPX8VWRVyiaWBrYHLUmcxM6uLrHrkElsDx0WwWYdimZnVWi/2yLcDfp86hJlZneRY\nyP+QYsO599lyzp9zdnD+1HLPX0U2hVxiHLAZcFXqLGZmdZJNj1xiR+DQCN7VwVhmZrXWaz1y98fN\nzAaQUyHfnkT9cci/z5Zz/pyzg/Onlnv+KrIo5BKvBtYFbkidxcysbrLokUvsCewbwS4djmVmVmu9\n1CNPdtihmVnd5VLIk/bHIf8+W875c84Ozp9a7vmrqH0hl5gMjAP+nDiKmVkt1b5HLvExYLsI9h6B\nWGZmtdaRHrmkqZLuknSPpEOGWG4LSYsk7T6csENI3lYxM6uzIQu5pDHA8cBUYENgmqQNBlnuGOAi\noO27AL20XgS8hxqcCJR7ny3n/DlnB+dPLff8VTQbkW8J3BsRsyJiIXA2xR16+vs0xT0053Y430bA\nMxE80OH1mpn1jGaFfE3gwYb52eVrS0hak6K4n1C+1Mmm+9uBP3VwfcMWEZenztCOnPPnnB2cP7Xc\n81fRrJBXKcrHAodGsddUdLC1QvGN4PoOrs/MrOeMbfL+Q8CkhvlJFKPyRpsDZ0sCmADsKGlhRJzf\nf2WSTgdmlbMLgBl9fy37+lgvn//1u2Hn0wZ/v6vz/9Y8b63ns83f2OOsQx7nr1e+XstfPp9eRp5F\nBUMefihpLHA3xZmVcyhGx9MiYuYgy58GXBAR5w7wXkuHH0osBzwOrBrB81U/N1IkTcn5K1rO+XPO\nDs6fWg/kb1o7hxyRR8QiSQdT3LV+DHBKRMyUdFD5/okdS/tKmwAz61DEIf8+W875c84Ozp9a7vmr\nqO0JQRKfATaM4OMjGMvMrNZyv2hWrXZ05n4sas75c84Ozp9a7vmrqHMh3wJff9zMrKlatlYkxlMc\nHTM+gkUjm8zMrL5ybq28DbjZRdzMrLm6FvLatVVy77PlnD/n7OD8qeWev4q6FvJa7eg0M6uzuvbI\nZwPbRnDfCMcyM6u1LHvkEmsAywL3p85iZpaD2hVyyv54REevoti23PtsOefPOTs4f2q556+itoU8\ndQgzs1zUrkcu8TvguAh+3YVYZma1VqV21qqQl7d2exzYIIJHuxLMzKzGctzZuR7Frd1qV8Rz77Pl\nnD/n7OD8qeWev4q6FfIt8PHjZmYtqVtr5Vjg4QiO6UooM7Oay7G14hG5mVmLalPIJcZQ3BXoptRZ\nBpJ7ny3n/DlnB+dPLff8VdSmkAPrAnMjeCp1EDOznNSmRy6xC/DJCHbsSiAzswzk1iPfELgzdQgz\ns9y4kFeUe58t5/w5ZwfnTy33/FXUqZBvAMxMHcLMLDe16JGXp+Y/BawdwfyuBDIzy0BOPfK1KE7N\ndxE3M2tRXQp5rfvjkH+fLef8OWcH508t9/xV1KWQuz9uZjZMdemRnwTMiOC/uhLGzCwTOfXIN6Dm\nrRUzs7pKXsjLI1Y2ouatldz7bDnnzzk7OH9queevInkhB14DBPBY6iBmZjlK3iOXmAIcGcHWXQli\nZpaRXHrktT/00MyszupQyLPY0Zl7ny3n/DlnB+dPLff8VdShkG9IzXd0mpnVWR165HOAt0fwYFeC\nmJllpPY9conxwIrA7JQ5zMxyVqmQS5oq6S5J90g6ZID3PyzpVkm3Sbpa0lsrbn8DYGYE3fla0Ibc\n+2w55885Ozh/arnnr6JpIZc0BjgemErRz54maYN+i90HbBsRbwW+AZxUcfvuj5uZtalpj1zSO4DD\nI2JqOX8oQET85yDLrwLcHhFr9Xv9FX0eiW9T3HD5mOH/CGZmvatTPfI14WU7ImeXrw3mY8CFFdYL\nPobczKxtYyssU7l/LendwP7AOwd5/3RgVjm7AC7aDHaYWb43BSAiLq/p/L8BM2qUZ9Tkb+xx1iGP\n89crX6/lL59PLyPPooIqrZWtgCMaWiuHAYsj4ph+y70VOBeYGhH3DrCel309kFgBmAusGMGLVcKm\nJGlK3y89Rznnzzk7OH9qPZC/aWulSiEfC9wNbAfMAa4HpkXEzIZl1gYuBfaJiGurhJHYHDg1go0r\n/jxmZqNOlULetLUSEYskHQxcDIwBTomImZIOKt8/EfgqsApwgiSAhRGxZZNVZ3FqvplZ3SU7s1Pi\nKOD5CL7elQBt6oGvZ9nmzzk7OH9qPZC/1md2ekRuZtYBKUfkfwY+HMGtXQlgZpah2o7Iy9u7rUPF\nQ2vMzGxwqVorqwKLIngy0fZblvv1GnLOn3N2cP7Ucs9fRapCPhmPxs3MOiJJj1xid2C/CP6lKxs3\nM8tUbXvkeERuZtYxqQr5OsADibY9LLn32XLOn3N2cP7Ucs9fhUfkZmaZS9UjvxWYHsEtXdm4mVmm\n6twjz661YmZWV10v5OUNl5cC5nd72+3Ivc+Wc/6cs4Pzp5Z7/ipSjMjXAR7I4YbLZmY56HqPXGI3\n4F8j2KUrGzYzy1hde+ST8RErZmYdk6y1kmC7bcm9z5Zz/pyzg/Onlnv+KjwiNzPLXIoe+c3AQRHc\n0JUNm5llrK49cl+H3Mysg7payCVWBJYF5nVzu52Qe58t5/w5ZwfnTy33/FV0e0TuY8jNzDqsqz1y\niF2AT0awU1c2amaWuTr2yCeT4aGHZmZ1lqK1MqvL2+yI3PtsOefPOTs4f2q5568ixYh8Vpe3aWbW\n07rdI78B+EwE13Zlo2Zmmatrj3xWl7dpZtbTul3IVwQe6/I2OyL3PlvO+XPODs6fWu75q+h2If9b\nBIu7vE0zs57W7R757yLYoSsbNDPrAXXskfsYcjOzDut2IZ/V5e11TO59tpzz55wdnD+13PNX4RG5\nmVnmut0j3zqCq7uyQTOzHlDHHvmsLm/PzKznNS3kkqZKukvSPZIOGWSZ48r3b5W06RCre3jYSRPL\nvc+Wc/6cs4Pzp5Z7/iqGLOSSxgDHA1OBDYFpkjbot8xOwHoRsT5wIHDCYOvL/BjyTVIHaFPO+XPO\nDs6fWu75m2o2It8SuDciZkXEQuBsYLd+y+wKnAEQEdcB4yWt3vGk6Y1PHaBNOefPOTs4f2q552+q\nWSFfE3iwYX52+VqzZdZqP5qZmVXRrJBXPaSl/x7VXryV2+TUAdo0OXWANkxOHaBNk1MHaNPk1AHa\nNDl1gJE2tsn7DwGTGuYnUYy4h1pmrfK1VygOQcyXpI+kztCOnPPnnB2cP7Xc8zfTrJDfCKwvaTIw\nB9gLmNZvmfOBg4GzJW0FLIiIR/uvqNlxkGZmNjxDFvKIWCTpYOBiYAxwSkTMlHRQ+f6JEXGhpJ0k\n3Qv8HfjoiKc2M7MlunZmp5mZjYwRP7OzyglFdSbpVEmPSro9dZZWSZok6TJJd0j6s6TPpM7UCknL\nSrpO0gxJd0o6OnWm4ZA0RtItki5InaVVkmZJuq3Mf33qPK2QNF7SOZJmlv9+tkqdqSpJbyx/533T\nk0P9/zuiI/LyhKK7ge0pdoDeAEyLiJkjttEOk7QN8AxwZkS8JXWeVkiaCEyMiBmSxgE3Af+S2e9/\n+Yh4VtJY4CrgCxFxVepcrZD0OWBzYMWI2DV1nlZIuh/YPCKeSJ2lVZLOAP4YEaeW/35WiIgnU+dq\nlaSlKOrnlhHx4EDLjPSIvMoJRbUWEVcC81PnGI6IeCQiZpTPnwFmAmukTdWaiHi2fLoMxX6arAqK\npLWAnYAf88rDdHORXW5JKwPbRMSpUOzvy7GIl7YH/jpYEYeRL+RVTiiyLiiPPNoUuC5tktZIWkrS\nDOBR4LKIuDN1phZ9D/giZHt5igB+L+lGSQekDtOCdYG5kk6TdLOkkyUtnzrUMH0I+NlQC4x0Ifee\n1Boo2yrnAJ8tR+bZiIjFEbEJxfkJ2+Z0ASRJ/ww8FhG3kOGotvTOiNgU2BH4VNlqzMFYYDPgvyJi\nM4oj6g5NG6l1kpYBdgF+PtRyI13Iq5xQZCNI0tLAL4CfRMR5qfMMV/m1+DfA21JnacE/AbuWfeaz\ngPdIOjOVb1SmAAACiklEQVRxppZExMPl41zglxTt0hzMBmZHxA3l/DkUhT03OwI3lb//QY10IV9y\nQlH5l2UvihOIrAskCTgFuDMijk2dp1WSJkgaXz5fDngvcEvaVNVFxJcjYlJErEvx9fjSiNgvda6q\nJC0vacXy+QrA+4Asjt6KiEeAByW9oXxpe+COhJGGaxrFIGBIzc7sbMtgJxSN5DY7TdJZwLuA1SQ9\nCHw1Ik5LHKuqdwL7ALdJ6iuAh0XERQkzteK1wBnlXvulgP+OiD8kztSO3FqNqwO/LMYDjAV+GhG/\nSxupJZ8GfloOIv9KZicrln88twea7pvwCUFmZpnr9q3ezMysw1zIzcwy50JuZpY5F3Izs8y5kJuZ\nZc6F3Mwscy7k1lPKk8+yOGnFrFNcyM2aKC+BalZbLuTWi8ZIOqm8mcbF5Q0qNpF0raRbJZ3bcOr/\n5ZI2L59PKK+LgqTpks6X9AfgEkkTJV1RXuT/dklbJ/z5zF7Ghdx60frA8RHxZmAB8AHgDOCLEbEx\nxfVCDi+XDQY/dX5T4AMR8W7gw8BF5ZUA3wrMGMH8Zi3xV0brRfdHxG3l85uA1wPjy5uEQFHUh7ws\naOl3EbGgfH49cGp5NcnzIuLWjiY2a4NH5NaL/tHw/EVgfL/3G68NvoiX/j9Ytt9yfXcn6rtT1DYU\nl2Y+XdK+nYlq1j4XchsNngSeaOhr7wtcXj6fxUvXON9jsBVIWhuYGxE/prht26YjktRsGNxasV7U\nv+cdwHTgR+Xtvhovafpt4H8lHUhx44po+EzjeqYAX5S0EHgayOa64tb7fBlbM7PMubViZpY5F3Iz\ns8y5kJuZZc6F3Mwscy7kZmaZcyE3M8ucC7mZWeZcyM3MMvd/FBgpaPT0XqgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa755c00b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Example 2/15 Hayler, 2012\n",
    "#Suppose that the battery failure time, measured in hours, has a probability density function\n",
    "#given by\n",
    "\n",
    "from __future__ import division\n",
    "def f(x):\n",
    "    return 2 / ((x+1)**3)\n",
    "\n",
    "def F(X):\n",
    "    return 1 - (1 / (X + 1)**2)\n",
    "#for x ≥ 0 and f (x) = 0 for x < 0. This is a valid probability density function because it is\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "a, b = 0, 1000 # integral area\n",
    "X = np.arange(0, 7, 0.1)\n",
    "y = pdf(X, f, a, b)\n",
    "plt.plot(X, y, linewidth=1)\n",
    "plt.grid(True)\n",
    "plt.title('PDF of battery failure')\n",
    "\n",
    "plt.figure()\n",
    "Y = F(X)\n",
    "plt.plot(X, Y)\n",
    "plt.xlabel('hours')\n",
    "plt.grid(True)\n",
    "plt.title('CDF of battery failure')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEWCAYAAACt5MYgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH7tJREFUeJzt3Xt0VdW9L/Dvb+2dnfeLIAkkQAISKY/yjCJ60AMWX/QW\nKKK8Kj0dVitW8N6eg/eeq8aKz2qrUUuPrWI9+Kin16Otij212loUEnkFCAESIEAIgbwT8iBk73n+\n2Du5IYQQSLJX5prfzxg6kp2115qbMfky85trziVKKRARkR4suxtARETdx9AmItIIQ5uISCMMbSIi\njTC0iYg0wtAmItIIQ5uISCMMbSIijTC0yQgiUiQis+xuB1FPMbTJFCrwH5HWGNpEvUxE/lFESkRk\naBfHPCAij4vI3cFsG+mPoU09JiJDROT/ichJETkoIj8OvD5SRCpEZFK748pEZEbg+wdFpFBEakUk\nT0TmtjtnkYj8RER2ikidiLwqIokiskFEakTkzyIS1+H4BwPnqRSR10Qk9GLa24v+DuCEUuroea4f\nC2AhgA8AbOzla5PDMbSpR0TEAvBHANsBDAEwC8AqEZmtlDoAYDWA9SISDmAdgHVKqS8Cby8EcK1S\nKgbAo4HjEgM/UwDmB853BYA5ADYAeBDAIPj77v0dmrMYwGwAIwGkA/i/F9PeHv5RtDcFwNdd/Pwq\nADuUUjlKqbxevC4ZgKFNPZUBYKBSao1SqkUpdQjAbwDcAQBKqd/AH845ABIB/GvrG5VSv1dKlQa+\nfhdAAYAr2537RaVUmVKqBP7R6yalVK5S6jSA/wQwqd2xCsBLSqljSqkqAI8DWHSx7W1PRCaIyD+J\nyNMi8h0RuUtEvtfZH4KIpIvIYyJyE4D/DeAv5znuKgArAbhFZF5nxxB1xW13A0h7wwEMEZGqdq+5\nAHzR7vvfwF8KuEspdab1xUAAPgAgNfBSFICB7d53ot3XjR2+bwoc3177csQR+EfSl9LeVokA9gG4\nUSm1WkQi4R+hv9H+oMDr7wK4XilVLSL/AuCzTs4HpVS2iDQCeJ6jbLoUDG3qqSMADiml0jv7oYhE\nAXge/uB+VETeU0pVichwAK8AmAn/CFqJyHYA0sW1uvoZAAzr8HXJxba3PaXUf4nIo/CXUwD/yL68\nk0PnA9gVCOwwAFFKqbIuTv0NAHsudH2izrA8Qj2VA6BORP5FRMJFxCUi40RkauDnLwDIUUr9EMBH\nAH4VeD0S/pJGOQBLRL4PYFwP2iEA7hWRZBEZAH8Z5p1LaG9HNwD4W+DrOwE8C7RNsrb+IzIQQG7g\n628B2Bwok3Q8DoGafblq9/SRjscQdaXL0BaRR0VkZbvvHxeRjpM/ZDCllA/+ScKJAA4CKIN/BB0j\nIt+Bf2LwR4HD/yeAySKySCm1B8BzADYBKIU/sC90J4Xq8HXH798C8F8ADsBfH19zMe3teGzgLo8B\nAGaKyF0AspVS7wV+/B+BcwDA2wBSRORm+AP8DIDoTo4D/JOQX3a4VMdjiM5LunrcWOBX2PeUUlMC\ns+77AWQEJnqI+g0ROQTgB0qpTmvJl3jOeQCmKaVW98K5pgC4C0AlgN8ppXIv8BaiTnVZ01ZKHQ7c\nZzsRQBKAbQxsMoGIjIb/N4NCEYlRStX28JReAMUAGhjY1BPdmYj8DYDvwz+T/lrfNoeof1BK7QXw\nD714vh0AdvTW+chcXZZHAEBEQgDshv+2qFHqQm8gIqI+c8GRtlLqjIh8BqCKgU1EZK8LhnZgAnIa\ngAUXOI6BTkR0CZRS3b7ls8vQFpEx8C8seC+wj0SvXZh6l4hkKqUy7W4HUbDp3vcvdsB7obtH9sC/\n+Q4REfUDXBFJRKQRhrZz/NXuBhDZ5K92NyCYLnjLX7dPJKJY0yYiujgXm50caRMRaYShTUSkEYY2\nEZFGGNpERBphaBMRaYShTUSkEYY2EZFGGNpERBphaBMRaYShTUSkEYY2EZFGGNpERBphaBMRaYSh\nTUTaEpEku9sQbAxtItKSyFPTgF8dtLsdwXbBB/sSEfU3Ik//I7DkUyDZuIGncR+YiPQm8tIC4Ht/\nmT17i5gYYeZ9YiLSlsi/3Qvc/u7ixV9hzZqhRj4pi+URItKCyKuPA/P/z333/R3Ll6ei9VGJEnhe\nl83NCxqGNhH1e5b12zdFbl2UmbkZt946HAAg0jrQfsS+htmAoU1E/ZbIK27LSvrKsv5halbWTrnq\nqiGdHLVHAHCkTURkJ5GXB1rWxF2hoSGD1q8/IsOHJ3RylA9AvBX4wgiciCSifkfk2SkiM48mJp66\nbMMGrzV8eLTdTeo3GNpE1K+IZP0AWJxz5ZX5ng8+SHBFRXm6OFoBCDMqx1geIaJ+QQQi8tZ6YNHi\nH/xgI370o+HdfOcAo279Y2gTke1EXoi3rPFbLWvS8Oef34pp04Z1850KQK1RoW3UrxVE1P+IPDtL\n5MbjAwf6hn30UYU1bdplF/FuBcBlVI4Z9WGJqP/wl0Pe+CWw/NOZM/d4PvwwzpWQEH4JZ4oyaqTN\n8ggRBZ3ICymWNXGjyLRhjz2Wjdmzh19i8CoADUaFNkfaRBRUIi/+UGR+UXJyfconn5TJ7Nk93RI7\nwqjQ5kibiIJC5N/iLCv5U5GFU5Yt+wr339/dycauKACNDG0iot4k8qsfitz0cnz8AVm7dg9GjOiN\nwAb8oW0ZVTFgaBNRnxF5brDLdcWfROaOW7x4o6xaldpuo6feEsqRNhFRT4hkCpD6LLBoVWJinnrp\npX0ybFhaH1xJAWhmaBMRXSqR52+yrDlvW1ZEzKpVOdYddwzt4ytypE1EdNFEnk+zrJHvAUsnXnvt\n39TjjydLeHhfB7YCcJqhTUTUXSJZsSJJbwLLbhk6NNv31FO7MGpUapCCVAEQTkQSEV2IyC/CgPiX\nge8uj40tUA8/nC0zZiS6gt8SD0fapJeQkP+R5/U+lS7S1GxZDaeBhkalmup8vqZKpZorgMZSoOEI\ncOoAUJsPHMpX6t0mu9tNehL5dTgQ/qLId5eHh5fixz/eZt122xAAdux5rQCcYWiTboYMS0/3uL/z\nnUPuykpfRHW1iq+uFtTURKCmJtpXWxvqO3UqGk1NsZbXm2ApFQORUmVZlWdEquqByjKvt7oIqC0E\n6nYDp7YA+3co9e4Zuz8Z9R8iz8QBiS+K3LwoPPwo7rlnq2vRoqEQ6ewRYMGiAIQwtEkvLldIS1JS\nCBYu7HQ5sIW27QrqANShuVmhqKheDhxo8Bw61Ow5ckTiS0oGjjp5MnVmXV28NDcnupSKhmUVeS3r\nRD1QdtzrrdwHVOUCVZuAAxuV+ve64H1CspPIy2NFkl8Gvj8jJmaf7957t7kWLEgG0FsLZHqqhaFN\nuhGIdP+5ph6PID09CunpUWefBHADpwEcQV1dC3btqnHl5TXGFBYi5vDhpFEnToy+pb4+yeXzDRLL\nOui1rON1SpUf9fmq9wDVOUDZXwD3TqUyjXnIqpOJvLDY5Rq7BliclpLyVcsDD+TIjBlJLiDZ7qa1\nowDD9lBiaFOnoqPdmD49AdOnt70UGLGXoaGhFDt21Lp27myI27cPcYcPDx5TVjZ+QWNjqgsAXK6t\njSKlJV5vWR5QFQjzwhyl3jXm4au6EnlpOJDwpMjE+SLzPFOnblP//M87kJqa1I+zguUR0o6g11cG\ndyEiwoXp0+MxfXp860suQEGpQygqakB2dl34rl2+kYWFQ9JKS8ff0tAwbI1SHrhcO5pESku93vLd\n/pF5w6dAWrZStzHMbeSfWPT+L5friruBpSkJCVtabr89371sWQpCQob280BUYHmENNQ/+qwIkJYW\ngbS0CNxxB4C20fkRFBc3YfPmmrDc3DOphYUDh5WUpN/U0DDsUaWixOXa3SRSetzrrcwDarMB36dA\nfA7DvO/4gxorLSvtn4DvXh4aeth3zTVF1r33bkdqarwbiL/gOfoHBcDdP/4CBAlD2xH6f59NSQnD\nggVhWLAAQFuYH8OxY/XYtKk6LDf3TFpBQfiw48dTbmpoGPpTpQaIy5XfBJw47vNV5gF1XwPyGTB2\ns1JTWuz9NHoSWTsSiFjpcg2dC3x3qMdz2JuRccBavvxrmTRpoKv/TCxeLLO6A0PbIS5mIrI/SU6O\nxIIFka1hHliYcRLFxQXYtKk6bOfOlrSCAs/Q48dH3NTQMDRTqURxufafFik94fWW5wM124GGz4Go\nvyl152kbP0q/I7IuHmi8U2TQPJHRGcCi8NjY3Jarrip3L1q0HePHx7mAvtjEKZi4IpI0pFRwa9rB\nkJISi9tui8VttwFo66cVOH78EDZvrgrdubNlWEGBJ7mkZOSsU6eSV/t8Q8SyClosq7Ta56sqUqp2\nL9CwHfB9CWRsdfroXAQCvDgV8MwXGXyDSNoY4LaI0NAC78iRx+SGGwqtefNqER0d5QaiLnxCbQhY\nHiENmdNnBw+Owbx5MZg3D0DbyLwaNTUn8PXXNe7du5sGFhZi4JEjAydWVg644/TpIS6lBoplHWyx\nrJO1SpUf8/mqC4Gm/UDzLiBkGxC3T6f6uci/JwK1M4HQa0USMyxrcDowIga4TaKjD7akp590XXtt\nodx6aw0GDAh3AYPtbnIf4kQkkZZiY0Nxww2DcMMNbS8F+vYJVFQcRG5urXvPnuYBBw9aA4qLE8ZU\nVESq+voBVktLogXEwLKOei2rvBGoqvL5ao4r1VQKNBcDzUeBlkOAOgBEFSj1vZq+/Bwij7qAmBFA\n2DggZDQQdjkQmupyxaYqNWiQz5ccCcwXt/uwLzb2pC8trco1YcJhmTGjEmPGDIBIqBvo6531+hOW\nR0hLziuP9KaEhEjMnBmJmTPbXgqM0P0LiSoq6rF37ynXgQPNUYcP+6JKS11DKypivDU1ob76+ghp\nbo61vN44USpBROoBnFIiDV6RphaRpjNA02n/f82nAXVGKZwBfC1K+byA1wf4FOB2i7jcgMsjYnkA\nlwdwe4CIMKUiPUpFhigV7QL+FUAtXK4TvrCwKl90dJ0aNOi0Ozn5lIwaVY8JE05gzJg4hIS4LSDa\nsme/D7ITQ9sRLm5FJJ0tISES11wTiWuuOetlF9rCvQlAKXy+4ygvb0RZ2WkpLz/trqo6466tVWG1\ntSq6rg5obBScOeNDS4uC1wt4vRZ8PoFSFlyuMwgJOY2QEEFIiEJYmIWICEFcXAUGDLBw2WUhSEoK\nQ1JSBMLDQwDAAkItIDTIfxq64YpIIjoPyxIMGhSBQYMi7G4KncVl1O+ZRv0L5VxG9VmidhQAr1F/\nARjajiC9/4BrIi2YNxFp1IclItIdQ9sROBFJplKAYfVBhrYjGNVniTrgRCQRkSYUAB9Dm3TDxTVk\npkC/NyrHjPqwzsWaNpEpGNpEpDGlYFiOGfVhnYu1ETIZJyJJO6xpk6k4EUlaYk2bTMX7tImIqB9j\naBORxszbmtWoD+tc3DCKTMaJSCIijXAikrTDiUgyFcsjpCWjBhpERmNoE5HGFGDYKgWGtiNwcQ2Z\nSgGwjOr9DG1HYE2bzOQfrHAikohIE5yIJC0ZNdAg6oDlEdIOa9pkMmVU72doE5HGjNtO27BP61BK\ncSKSTMaJSNKOUX2WqB1ORJKWWNMmYylORBIRaUOBE5GkIRHWtMlMXBFJWjKqzxK18W8k77O7GUHF\n0CYijXEikrTEiUgyGcsjpB3ep02m4kQkEZFGuJ82acmoPkvUgVm/ZTK0HYE1bTIVJyKJiDQi4EQk\naYgTkWQqTkSSlozqs0RtAoMVo/4CMLSJiDTC0HYETkSSqTgRSXrihlFkME5EEhFphBORpB2j+ixR\nO5yIJC2xpk1kCoa2I/A+bTJToN8blWNGfVgicho+uYa0ZFSfJWpHFCciSUOsaZOpOBFJRET9GEPb\nETjSJlNxIpK0JDBtI3gioPXuEU5EEhFpxKwBC0PbEVgeIaMZlWNGfVjnYmgTmYKhTUQa4y1/pCGl\nREyr6xH5cSKStMTyCJlJRMAVkURE2uB92qQljrSJTMHQdgQuriFTcSKSNMWRNpmJE5GkJaP6LFEH\nnIgk7bCmTWYK9Hujej9D2zFY0yYyAUObiDTGW/5ISyyPkKkUTIsxsz6tYzG0yVRcEUnkOO+++y6u\nu+467N692+6mUC/zPwTBrBwz6sM6FzeM6sqtt96K0NBQjB071u6mEPUYQ9sRWB7pypYtWzBp0qTA\n5kLkLFwRSeQ4mzdvhojgk08+wRNPPIHCwkK7m0S9RimuiCQNcaTd6p133sGsWbOwdOlSHD58GACQ\nk5ODJUuW4KabbsKMGTPwy1/+0uZWUu/hRCRpiRtGAf4yyC9+8Qs8//zzqK+vx2OPPYbS0lIopTB+\n/HgAQHl5Oaqrq21uKfWWwEQkQ5tIR1lZWbj66quRnp4OpRQSExORn5+PiRMnth2TnZ2N6dOn29hK\n6l0Kpm2XxtB2BJZHdu/ejfz8fHzrW99CaGgo3n//fTz++OOIjIxEdHQ0AODIkSMoLCzE0qVLbW4t\n9R6WR0hLDO2PPvoIAM4ZRWdkZMCyLHz44Yd4++23sXbtWoSFhdnRROoT5q2IdNvdAKLe8MUXX2DE\niBGIj48/63URwcqVKwEAc+bMsaNp1If8gxWOtEk7Zk9EHjlyBCdPnjyrdk2m4IpI0pLZ5ZGvv/4a\nADBu3DibW0L2MKvzM7QdwezQ3rJlCwDgG9/4hs0tIXuwPEKklW3btsHj8SAtLc3uplCQ+e/TFqNy\nzKgP61zSusjAOIcPH0ZlZSVGjhwJl8tld3OI+hxD2xmM+vWwve3btwMArrjiCptbQjYyqv8ztElr\n27ZtAwBcfvnlNreE7MEVkaQlcycid+3aBaB/hLbX67XlvSYL3KdtdzOCiqHtCGbWtKuqqlBcXAwR\nwciRI21ty2effYaPP/74kt+/bt26tn+A6GJwIpJIGzt37gQAxMfHIy4urs+vd/ToUaxatQpZWVl4\n8sknoZT/H8qtW7di+/bt+Pa3v33J516+fDleffVVFBUVdev4n/zkJ1iyZAlXeZpVzgbA0HYEZdZt\nqm1aQzsYpZEzZ87gvvvuw6xZs1BRUYEPPvgA9fX1OHXqFLKysnDffff16PxutxsPPvggHnnkkW6V\nSp5++mlMnjwZJ06c6NF19cetWUlLZta0W8sJo0aN6vNrbdq0CSUlJZg8eTIWLlyIrKwsREVFYd26\ndbj55psRGhra42skJSVhxIgR+PDDDy94rMvlQnp6eo+vqbvAfdpG9X6GNmnJ6/UiPz8fQHBG2lu3\nbkV8fDySk5MxduxYXHnllWhsbMT777+PW265pdeuc/vtt+O3v/1tr53P+bg1K2nJvInIoqIiNDU1\nQUSCMuLMy8vDmDFjznpt48aNGDJkCGJiYnrtOunp6aiursa+fft67ZzOZt5Im1uzOoJRfRYAsGfP\nHgD+MsGIESP67DqZmZmorKxEbm4uUlNTcf/99yM5ORmrV69GdnY2vvnNb573vfn5+fj4449hWRZK\nSkrw0EMP4b333kNdXR1OnjyJu+++GykpKWe9x7IsTJgwAZs2bTprwdCBAwewbt06xMTEIDQ0FB6P\n55xtaHt6bdIDQ5u01BraaWlpcLv7rhtnZmbi2LFjmDt3LlasWIHrr7++7Wf79+/HvHnzOn1fcXEx\n/vCHP2D16tVt51m+fDkeffRR+Hw+3HXXXRg9ejSWLFlyznuHDx+OgoKCtu9zc3OxcuVK/PznP8fk\nyZMBAA0NDbj33nshnQwye3Jt3fAZkaQp8yYiW0N79OjRfX6t1lJFxzJMSUlJ26PMOnrzzTdx//33\nt33f2NiI2NhYjB8/HklJSViyZMl5bxGMjo5GSUkJAMDn8yEzMxMZGRltgQ0AERERmD17dttth711\nbT2Z1fsZ2o5gVk3b6/WisLAQQHC2Y92/fz+ioqIwZMiQs14/deoUoqKiOn3PsmXLEB4e3vb9rl27\ncOWVVwIAEhMTsXLlyvPWwmNjY1FXVwfAf1tjcXExJkyY0O329uTaeuJEJFG/VlRUhObmZohIUEJ7\n3759nU52ikinI10AZwV8UVERysrKMHXq1G5dTynVdt7W+7AvJmR7cm3d8JY/0pRZnXb//v0A/AtS\ngrG7X0FBQaehHR0djZqamgu+f8uWLQgJCTlr0vLYsWPnPb6mpqat7JKYmAgAaGpquthmX9K19WNU\n1wfA0HYIs2raraF9+eWXIyQkpE+vVVNTgxMnTnQa2kOGDOk0tJuamvDCCy+0lXCys7MxatSotgU4\nPp8P69evP+81a2trkZycDABtdejdu3efc1xnKyd7em398BmRpCWzatoHDhwAAIwdO7bPr9U6CdnZ\nqsuJEyfi0KFD57z+5ZdfYv369Th48CCKiopw9OhReDyetp+vW7euyz1DDh061Fb2cblcePjhh7Fx\n48az7igpLy/HH//4RwBnj5x7em3dBMojdjcjqHjLH2mndRQZjNDeu3cvoqOjOx1pX3311XjuuefO\neX3KlCmYM2cO8vPzsW/fPrz++ut46qmn8MQTT8Dj8eD6668/b9tbWlqwc+fOs+7+yMjIQFZWFn79\n619j8ODBiIiIgNvtxi233II33ngDq1atwpIlSzB37tweXVtP0vY/UzC0SSt1dXUoKyuDiATl6et7\n9+5FRkYGLOvcX0onTZqEiooKlJWV4bLLLmt7PS4uDo888shZx2ZmZnbrenl5eUhMTDxnZD9u3Dg8\n88wz5xy/fPnys77vybX1xPu0SUvm1LRbSyMxMTFITU3tk2u8/vrrWLFiBQD//eCzZs3q9DiPx4OF\nCxfi7bff7rVrv/XWW1i6dGmvnc/pTOn37TG0SSsHDx4E4K8n95UNGzbA4/GgoKAAISEh5w1tALjz\nzjvx1Vdfoba2tsfXLSoqQklJicMWvvQ1TkSSlswZabeG9qRJk/rsGsuWLcPAgQPx2muv4Wc/+1mX\nT3kPCwvDQw89hDVr1pz3nu3uOH36NJ555hmsWbOm06Xp1BWz/rxY03YEgSnPyWu93a8vR9pz5sy5\nqDssxo4di/nz5+N3v/sd7rjjjku65rp167BixQoMHz78kt5vOKNSm6FNWiksLER4eHhQ9hy5GNOm\nTcO0adMu+f333HNPL7bGHNwwijRlRnnk+PHjqKurw7hx47osWRA5GUPbEcwI7dYn1UyZMsXmllA/\nY1SOGfVhSW95eXkA0LZjHZGJKyIZ2o5gxkTkrl27EBkZ6bAVfdQz5q2IZGg7gjj+LrGmpibk5eXh\nqquu6nR1IpmKE5GkKaeH9pYtW9Dc3IzrrrvO7qZQP2LSRmmtGNqO4LzEfvbZZ7Fo0SK0tLQAAD75\n5BPExMR0uTqRTCVG5ZhRH9a5nFfTzsnJQVNTE3w+H0pLS/HZZ59h0aJFbftCE/1/zhu0dIWLaxzC\naeWRCRMmYMCAAaitrcVPf/pTDBs27Jwd7Yj8/Z7PiCSy3YoVK5CXl4e5c+fC4/HgxRdfhNvd+Rij\npaUFa9euxe9//3u88847eOCBB1BcXBzkFpM9zJuI5EjbEZy3uCYuLg4vvfRSt4598sknMWrUKCxY\nsADV1dV45ZVXkJKS0sctpP7Aaf2+OzjSdgTn1bS7q6CgAH/+858xf/58AP69SSZPnmxzqyi4zIpu\nhjZpLTs7GxMnTmx7DmJOTg6mTp2Kuro6m1tGwcPQJs0o5bzySHfFxsYiISEBANDQ0IDPP/8cEyZM\nwIYNG2xuGQWDiRORrGk7gllPY2/vxhtvxI4dO/CnP/0Jzc3NuPHGG/Hll1/2u61bqa9wIpJIKx6P\nBw899JDdzSCb+EfaZv2eyfKII5g7EUlkGoa2Izh/wyiizplXHmFoE5G2TBytMLQdwdy7R8h0Ctww\nijTEmjaZScx7BgJDm4h0Z9Z92gxtR2B5hEzFiUjSFEObTMSJSNIYa9pkIgXThiwMbUcwqs8SteGK\nSNIUa9pkLDHtt0yGtiMwtMlUSvE+bSIiTZg4WGFoOwIX15DRjIpuhrYjsDxCZuKKSCIi7XBFJGmH\nW7OSmUSUmFbZZmg7AmvaRKZgaBOR7jjSJt1wIpLMxBWRpC2zui3RWYzq/QxtItKWCLdmJS1xIpLM\nZOJvmAxtR2BNm4xmVO9naDsEQ5vMZVbnZ2g7glmdlqhVYLBi1F8AhrYjsKZNpuJEJGmK5REykYn9\nnqHtCAb2XKI2ZkU3Q9sRuGEUmYkrIklLSrGmTcbiMyKJiHTBFZGkKS6uITIFQ9shGNpkLrN6P0Ob\niLTFiUjSlLTW9oiM4o9rPiOStCMoLMy1uxFENtmSYHcLgomh7RCFhTvtbgJR0PlH2rsZ2qQbo347\nJGrHrNIIALh782Qi+HZvno+6RyTSKis7iS+++MLuphAFVV1dEoCBkSZljyjVOxNYwpkwIqJLolT3\nf2PotdAmIqK+x5o2EZFGGNpERBphaBORlkQkVUR22d2OYGNoExFphKGtORH5TxHZIiK7ReQuu9tD\nFGRuEVkvIntE5D9EJNzuBvU13j2iORGJV0pVBTprDoDrlFKVdreLqK+JSCqAgwCuUUptEpFXAexR\nSj1na8P6GEfa+lspIjsAbAKQAmCUze0hCqajSqlNga/XA7jWzsYEQ6+uiKTgEpHrAcwCME0p1SQi\nnwMItbdVREHVvlRgxHP3ONLWWwyAqkBgjwYwze4GEQXZMBFp7feLAfzdzsYEA0Nbb5/APxGzB8CT\n8JdIiEyhAOwDsCLwdyAWwFp7m9T3OBFJRKQRjrSJiDTC0CYi0ghDm4hIIwxtIiIbiUiGiOSKSKiI\nRAZWN4857/GciCQispeIPAYgDEA4/AuGnj7vsQxtIiJ7iUgIgC0AGgFcrboIZpZHiIjsNxBAJIAo\n+Efb58WRNhGRzUTkDwDeAjACwGCl1I/Pdyz3HiEispGIfA/AaaXUOyJiAfhKRK5XSv210+M50iYi\n0gdr2kREGmFoExFphKFNRKQRhjYRkUYY2kREGmFoExFphKFNRKQRhjYRkUb+G5yEiHl73Z3TAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa755ebd750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# source: http://matplotlib.org/1.2.1/examples/pylab_examples/integral_demo.html\n",
    "# modified jfb to use as pdf example, but currently no accurate (e.g. f(x) does not satisfy all conditions)\n",
    "# implement the example graphs/integral from pyx\n",
    "from matplotlib.patches import Polygon\n",
    "\n",
    "def pdf(X):\n",
    "    \"\"\"\n",
    "    define a pdf based on func\n",
    "    \"\"\"\n",
    "    val = []\n",
    "    for x in X:\n",
    "        if x <= 4:\n",
    "            val.append(0)\n",
    "        elif x >= 8:\n",
    "            val.append(0)\n",
    "        else:\n",
    "            val.append(func(x))\n",
    "    return np.array(val)\n",
    "\n",
    "def func(x):\n",
    "    return (x-3)*(x-5)*(x-7)+85\n",
    "\n",
    "def Fx(pmf):\n",
    "    return np.cumsum(pmf)\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "a, b = 4, 8 # integral area\n",
    "X = np.arange(0, 10, 0.01)\n",
    "y = pdf(X)\n",
    "plt.plot(X, y, linewidth=1)\n",
    "\n",
    "\n",
    "# make the shaded region\n",
    "ix = np.arange(a, b, 0.01)\n",
    "iy = func(ix)\n",
    "verts = [(a,0)] + list(zip(ix,iy)) + [(b,0)]\n",
    "poly = Polygon(verts, facecolor='0.8', edgecolor='k')\n",
    "ax.add_patch(poly)\n",
    "\n",
    "plt.text(0.5 * (a + b), 30,\n",
    "     r\"$\\int_a^b f(x)\\mathrm{d}x$\", horizontalalignment='center',\n",
    "     fontsize=20)\n",
    "\n",
    "plt.axis([2,9, 0, 110])\n",
    "plt.figtext(0.9, 0.05, 'x')\n",
    "plt.figtext(0.1, 0.9, 'y')\n",
    "ax.set_xticks((a,b))\n",
    "ax.set_xticklabels(('a','b'))\n",
    "ax.set_yticks([])\n",
    "plt.title(r'example $p.d.f.$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Main Descriptors of Random Variables\n",
    "------------------------------------\n",
    "\n",
    "The $p.m.f.$ and $c.d.f.$ provide useful tools for describing the\n",
    "probability properties of a random variable, however, there are some\n",
    "general characteristics that are also of interest. These include:\n",
    "\n",
    "-   the *central* or *Expected Value*\n",
    "-   the *Variance*\n",
    "\n",
    "While these characteristics are useful in describing a random variable,\n",
    "it is important to keep in mind that two random variables may have the\n",
    "same mean, median, or otherwise, but their probability distributions may\n",
    "be quite different.\n",
    "\n",
    "### Expected Value of a Random Variable\n",
    "\n",
    "#### Discrete Case\n",
    "\n",
    "The expected value or expectation of a *discrete* random variable with a\n",
    "**probability mass function**, \\$ P(X = x\\_i ) = p\\_i\\$ is:\n",
    "\n",
    "$$E(X) = \\sum_x {p_i  x_i}$$\n",
    "\n",
    "$E(X)$ provides a summary measure of the average value taken by the\n",
    "random variable and is also known as the **mean** of the random\n",
    "variable.\n",
    "\n",
    "It is possible to think of the *expected value* as a *weighted average*\n",
    "withing the sample space, $S$, where the weights are the probabilities\n",
    "$x_i$.\n",
    "\n",
    "#### Continuous Case\n",
    "\n",
    "For the *continuous* case, the expected value is simply the integral of\n",
    "the $p.d.f.$:\n",
    "\n",
    "$$E(X) = \\int_{-\\infty}^\\infty {x \\cdot f(x)dx}$$\n",
    "\n",
    "Some notes for when we cover the *Method of Moments* - The expected\n",
    "value of $(x- \\mu)^i$ is equal to the $i^{th}$ central moment\n",
    "$E[(x- \\mu)^i ] = μ_i$ ; later on we will show that when $i = 2$, $μ_i$\n",
    "is the variance.\n",
    "\n",
    "-   The expected value of $x^i$ is equal to the $i^{ th}$ moment about\n",
    "    origin $E[x^i ] = μ_i$\n",
    "\n",
    "### Measures of Central Tendency\n",
    "\n",
    "#### Arithmatic Mean\n",
    "\n",
    "A sample estimate of the population mean is the arithmetic average.\n",
    "\n",
    "$$\\bar{X} = {1\\over n} \\sum_i {x_i}$$\n",
    "\n",
    "#### Geometric Mean\n",
    "\n",
    "Used for the lognormal distribution (we'll cover this later). Evaluates\n",
    "the *central tendency* by using the product of values (as opposed to the\n",
    "arithmetic mean which uses their sum). The geometric mean is defined as\n",
    "the $n^{th}$ root of the product of $n$ numbers:\n",
    "\n",
    "$$\\bar X_G = \\left(\\prod_{i=1}^n x_i \\right)^{1/n} = \\sqrt[n]{x_1 x_2 \\cdots x_n};$$\n",
    "\n",
    "Where, $\\prod x_i = x_1 x_2 x_3 ... x_n$;\n",
    "\n",
    "For instance, the geometric mean of two numbers, say 2 and 8, is just\n",
    "the square root of their product; that is $\\sqrt{2\\cdot 8}=4$. As\n",
    "another example, the geometric mean of the three numbers 4, 1, and 1/32\n",
    "is the cube root of their product (1/8), which is 1/2; that is\n",
    "$\\sqrt[3]{4\\cdot 1\\cdot 1/32}=1/2$.\n",
    "\n",
    "It is good when comparing, or using together random values of different\n",
    "magnitude, as it naturally *normalizes* the values so that they may be\n",
    "compared.\n",
    "\n",
    "#### Weighted Mean\n",
    "\n",
    "is the quantity\n",
    "\n",
    "$$\\bar{x} = \\frac{ \\sum_{i=1}^n w_i x_i}{\\sum_{i=1}^n w_i},$$\n",
    "\n",
    "which means:\n",
    "\n",
    "$$\\bar{x} = \\frac{w_1 x_1 + w_2 x_2 + \\cdots + w_n x_n}{w_1 + w_2 + \\cdots + w_n}.$$\n",
    "\n",
    "Therefore data elements with a high weight contribute more to the\n",
    "weighted mean than do elements with a low weight. The weights cannot be\n",
    "negative. Some may be zero, but not all of them (since division by zero\n",
    "is not allowed).\n",
    "\n",
    "The common mean $\\frac {1}{n}\\sum_{i=1}^n {x_i}$ is a special case of\n",
    "the weighted mean where all data have equal weights, $w_i=w$. When the\n",
    "weights are normalized then $w_i'=\\frac{1}{n}.$\n",
    "\n",
    "#### Median\n",
    "\n",
    "The median is that value of $x$ for which $P(X≤x) = P(X≥x) = 0.5$. In\n",
    "the case of a **continuous** distribution the median corresponds to an\n",
    "ordinate which separates density curve into two parts having equal areas\n",
    "of 1⁄2 each.\n",
    "\n",
    "$$\\int_{-\\infty}^{\\mu_{md}} {f(x)dx} = \\int_{\\mu_{md}}^{\\infty} {f(x)dx} = 0.5$$\n",
    "\n",
    "Simply, there are equal chances of selecting data above and below the\n",
    "median.\n",
    "\n",
    "For **discrete** random variable, the median is the value halfway\n",
    "through the *ordered data set*, below and above which there is an equal\n",
    "number of data values.\n",
    "\n",
    "$$\\mu_{md} = x_p$$\n",
    "\n",
    "Where; $p$ is determined from $\\sum_{i=1}^p f(x_i) = 0.5$\n",
    "\n",
    "#### Mode\n",
    "\n",
    "Mode is simply a measure of where the most values are located in the\n",
    "distribution. However, it is not necessarily equal to the mean or\n",
    "median. The shape of the distribution has a greater impact on the median\n",
    "and mean, and their lack of equivalence to the mode can indicate extreme\n",
    "values.\n",
    "\n",
    "for the **continuous** case, the mode is equal to the location where the\n",
    "$p.d.f.$ reaches a maximum value:\n",
    "\n",
    "$${df(x) \\over dx } = 0;\\  and \\  {d^2f(x) \\over dx^2} \\lt 0$$\n",
    "\n",
    "For **discrete** variable mode is the x value associated with\n",
    "$Max_{i=1}^n f( x_i )$\n",
    "\n",
    "*A sample or population may have none, one or more than one mode. (Can\n",
    "you give an example of each case?)*\n",
    "\n",
    "### Measures of Dispersion\n",
    "\n",
    "The term dispersion refers to the spread of the data. We are interested\n",
    "in understanding something about the *variability* of the data. If we\n",
    "are told a distribution has an expected value, $E(X)$, of 0.5 --does\n",
    "that give us any information about the *shape* of the distribution? To\n",
    "have an idea of the probability properties of a distribution, we need to\n",
    "understand more information than just the *central tendencies*. Some\n",
    "common measures of dispersion include: **Range**, **Variance**, and the\n",
    "**Coefficient of Variation**.\n",
    "\n",
    "#### Range\n",
    "\n",
    "The simplest measure, it is the difference between the maximum and\n",
    "minumum values\n",
    "\n",
    "$$Range = X_{max}  - X_{min}$$\n",
    "\n",
    "#### Variance\n",
    "\n",
    "Variance provides a more robust measure of the variability of a\n",
    "population. The variance is a non-negative number given us insight into\n",
    "how widely distributed the values of the random variable are likely to\n",
    "be. If we have very high variance, it decreases the 'chance' that we\n",
    "will draw the central value. In other words, the $x_i$ are more\n",
    "scattered within the the distribution. Think of *Variance* \\~\n",
    "*variability*. We symbolize variance with $Var(X) \\ or \\ V(X)$ or\n",
    "$\\sigma^2$\n",
    "\n",
    "$$V(X) = \\sigma^2 = E[X-E(X)]^2 = E(x-\\mu)^2 = \\int (x-\\mu)^2 f(x)dx$$\n",
    "\n",
    "Note that the variance is also the **second central moment** or , $u_2$.\n",
    "\n",
    "For a **discrete** variable, we can calculate variance:\n",
    "\n",
    "$$\\sigma_x^2 = {\\sum{(X_i - \\mu)^2 } \\over n }$$\n",
    "\n",
    "We can make a sample estimation of *population* variance from\n",
    "observations:\n",
    "\n",
    "$$S^2 = {\\sum_{i=1}^n{(X_i - \\bar x)^2 } \\over n-1 }$$\n",
    "\n",
    "Comparing these equations, not two distinct differences: (1) $\\mu$\n",
    "(*population mean*, which we cannot know) is replaced by $\\bar x$\n",
    "(*sample* mean, an estimate of $\\mu$), and (2) $n$ is replaced by $n-1$\n",
    "in order to have an unbiased estimation of population variance.\n",
    "\n",
    "Take it as an exercise to prove this by deriving:\n",
    "$E(S^2) = ... \\sigma^2$\n",
    "\n",
    "Often it is convenient to work with the **standard deviation** rather\n",
    "than the variance because it has the *same units* as the original\n",
    "variable. We can find the standard deviation simply from the variance:\n",
    "\n",
    "$$\\sigma = \\sqrt{V(X)} = \\sqrt { \\sigma^2}$$\n",
    "\n",
    "or from a *Sample Estimation*:\n",
    "\n",
    "$$S = \\sqrt {{\\sum_{i=1}^n {x_i - \\bar x)^2 }} \\over n-1}$$\n",
    "\n",
    "#### Coefficient of Variation\n",
    "\n",
    "The coefficient of variation measures the spread of a set of data as a\n",
    "proportion of its mean. It is often expressed as a percentage. It is\n",
    "dimensionless, so it can be used for cross comparison, i.e. for\n",
    "comparison of the variability of the same variable at different places,\n",
    "etc. It is the ratio of the sample standard deviation to the sample\n",
    "mean:\n",
    "\n",
    "$$Cv = {S \\over \\bar x }$$\n",
    "\n",
    "There is an equivalent definition for the coefficient of variation of a\n",
    "population, which is based on the expected value and the standard\n",
    "deviation of a random variable:\n",
    "\n",
    "$$Cv = { \\sqrt(V(X)) \\over E(X)}$$\n",
    "\n",
    "### Measures of Symmetry\n",
    "\n",
    "Two further *moments* of the distribution exist, these are **skewness**\n",
    "and **kurtosis**. They provide further information regarding the shape\n",
    "of the distribution. Recall, all four of these moments or measures are\n",
    "simply parameters that describe the shape of the distribution function.\n",
    "We are concerned with what we may expect when we 'pull' a sample from a\n",
    "population, given a distribution function, we have a better idea of what\n",
    "to expect... beyond simply an expected value. So, as we think about the\n",
    "shape (so far, central tendency and width), now we look at two measures\n",
    "that give information on the left/right-ness of the distribution and the\n",
    "peakedness.\n",
    "\n",
    "#### Skewness\n",
    "\n",
    "Skewness is a parameter that describes asymmetry in a random variable’s\n",
    "probability distribution.\n",
    "\n",
    "$$skewness = \\int(x-\\mu)^3f(x)dx = \\mu_3$$\n",
    "\n",
    "; it is also the **3rd moment** of a population.\n",
    "\n",
    "**Coefficient of skewness** (dimensionless, more commonly used). Several\n",
    "types of skewness are defined, the terminology and notation of which are\n",
    "unfortunately rather confusing. \"The\" coefficient of skew of a\n",
    "distribution is defined to be:\n",
    "\n",
    "**population skewness** is:\n",
    "\n",
    "$$\\gamma = {\\mu_3 \\over (\\mu_2)^{3/2 }} = {\\mu_3 \\over \\sigma_3}$$\n",
    "\n",
    "Where $\\mu_3$ is the 3rd moment about mean, and $\\sigma$ is the standard\n",
    "deviation.\n",
    "\n",
    "**Sample unbiased estimate** is:\n",
    "\n",
    "$$Cs = { n \\sum (x_i - \\bar x)^3 \\over (n-1)(n-2)S^3 }$$\n",
    "\n",
    "*NOTE*: if $Cs \\gt 0$, Positive skew (long tail in the right, most\n",
    "hydrological variables, e.g. discharge hydrograph, rainfall, etc). If\n",
    "$Cs \\lt 0$, negative skew (tail left).\n",
    "\n",
    "#### Kurtosis\n",
    "\n",
    "Kurtosis is a parameter that describes how peaked is a random variable’s\n",
    "probability distribution.\n",
    "\n",
    "$$kurtosis = \\int(x-\\mu)^4f(x)dx = \\mu_4$$\n",
    "\n",
    "; it is also the **4th moment** of a population.\n",
    "\n",
    "For a **sample** of $n$ values the *sample excess kurtosis* is:\n",
    "\n",
    "$$k_2 = \\frac{\\mu_4}{\\mu_{2}^2} -3 = \\frac{\\tfrac{1}{n} \\sum_{i=1}^n (x_i - \\overline{x})^4}{\\left(\\tfrac{1}{n} \\sum_{i=1}^n (x_i - \\overline{x})^2\\right)^2} - 3$$\n",
    "\n",
    "where $\\mu_4$ is the fourth sample moment about the mean, $\\mu_2$ is the\n",
    "second sample moment about the mean (that is, the **variance**), $x_i$\n",
    "is the $i^{th}$ value, and $\\overline{x}$ is the **sample mean**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References / Notes\n",
    "\n",
    "Concepts and Examples are based on\n",
    "\n",
    "1. Ang H-S.A., Tang W.H.; Probability Concepts in Engineering Planning and Design, Volume 1: Basic Principles <http://books.google.no/books/about/Probability_Concepts_in_Engineering_Plan.html?id=EIRRAAAAMAAJ&redir_esc=y>\n",
    "2. Haan, C.; Statistical Methods in Hydrology 2nd. Edition. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coda\n",
    "__Course__: GEO4310 Statistical and Stochastic Hydrology, University of Oslo\n",
    "\n",
    "__Contributions__\n",
    "\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)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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
}