ISDAwPython_3_1_NB_1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.10"
    },
    "colab": {
      "name": "ISDAwPython-3.1-NB-1.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/sigvehaug/ISDAwPython_day3.1/blob/main/ISDAwPython_3_1_NB_1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4EjvQpwJdKUT"
      },
      "source": [
        "ISDAwPython-3.1-NB-1, Statistics with Python, S. Haug, University of Bern. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ome17cCJdKUa"
      },
      "source": [
        "\n",
        "# Statistics with Python\n",
        "\n",
        "In this notebook we refresh some basic statistics and practice how to use Python to get them."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lbgBKkuoduzs"
      },
      "source": [
        "## Warm up and discussion (homework)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "teUxbdU6eXu3"
      },
      "source": [
        "Let's look at over- and underbinning of histograms. So first take a look at the documentation page of the Scipy.stats module for the normal distribution: \n",
        "\n",
        "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm\n",
        "\n",
        "and use the norm.rvs() function to generate 5000 random values from the normal distribution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y59ULkTCeo7F"
      },
      "source": [
        "from scipy.stats import norm\n",
        "import matplotlib.pyplot as plt\n",
        "# Generate 5000 random values\n",
        "# Plot the generated data as histograms with different binnings\n"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3culwpXefOpb"
      },
      "source": [
        "How do you decide on which binning to use?\n",
        "\n",
        "\n",
        "*   Not smaller bins as the accuracy of your data\n",
        "*   Or let Python decide\n",
        "*   Plot with correct uncertainty bars (the best) \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OWL0J7zVgM2g"
      },
      "source": [
        "# 1. Descriptive Statistics\n",
        "\n",
        "We divide the statistics field into descriptive and inferential statistics. With descriptive statistics we describe our data. With inferential statistics we draw conclutions from the data. At the heart of statistics there are random variables and probability density/mass functions."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UqpHakxtdKUh"
      },
      "source": [
        "## 1.0 Random variables and probability density functions (pdf)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pl0cfeQjPAgV"
      },
      "source": [
        "In practice the measurement process, the data taking, is a **random, or stochastic, process**. The outcome varies from measurement to measurement. There are three (at least) reasons:\n",
        "\n",
        "- Measurements are normally on a sample, not the full population. Samples fluctuate.\n",
        "- Sensors have limited resolution, measurements on the same sample vary within the resolution \n",
        "- According to quantum mechanics, i.e. at smallest distances, measurements are by nature stochastic\n",
        "\n",
        "Our Iris dataset is a sample from 50 flowers in each class. So in each class there are 50 varying measurements for each of the four observables, sepal and petal length and width. This is due to the first and maybe the second reason (quantum mechanics can be neclected at scales larger than molecules). In descriptive statistics the observables are therefore called **random** variables. Let us call one x for examplification.\n",
        "\n",
        "If x can take on any value from a continuous range, we write $f(x;\\theta)dx$ as the probability that the measurement’s outcome lies between x and $x + dx$. The function $f (x; \\theta)$ is called the **probability density function (p.d.f.)**, which may depend on one or more parameters $\\theta$ (for example the Iris class).\n",
        "\n",
        "A random variable can be discrete or continuous. If discrete, then we use $f(x;\\theta)$ to denote the probability to find the value x (in python the term probability mass fundtion, pmf, is then used). In the following the term p.d.f. is often taken to cover both the continuous and discrete cases, although technically the term density should only be used in the continuous case.\n",
        "\n",
        "The p.d.f. is always normalized to **unity** (the number 1), i.e. the integral, i.e. the surface under the curve equals one. Both x and $\\theta$ may have multiple components and are then often written as vectors. If $\\theta$ is unknown, we may wish to estimate its value from a given set of measurements of x; this is a central topic of statistics (see next notebook on parameter estimation and regression).\n",
        "\n",
        "The p.d.f. should be chosen to describe the fluctuation of the random variable in a best possible way. In other words, we should always choose an approprate p.d.f to describe our data. The most used p.d.f. is the normal distribution. Examples of other pdfs can be found in the appendix."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LvcFjoRudKUi"
      },
      "source": [
        "## 1.1 The normal pdf\n",
        "\n",
        "The normal (or Gaussian) probability density function is probably the most used one (informally the \"bell curve\"). It \n",
        "derives its importance in large part from the *central limit theorem*: \"In most situations, when independent random variables are added, their properly normalized sum tends toward a normal distribution (informally a \"bell curve\") even if the original variables themselves are not normally distributed.\" https://en.wikipedia.org/wiki/Central_limit_theorem \n",
        "\n",
        "**Example:** If one flips a coin many times the probability of getting a given number of heads in a series of flips will approach a normal curve, with mean equal to half the total number of flips in each series. (In the limit of an infinite number of flips, it will equal a normal curve.)\n",
        "\n",
        "This means that in many or most cases it is sufficient to know the characteristics of the normal p.d.f. Others can be looked up if needed. Also often unspecified statements like the *error*, or better, the *uncertainty* refer to their meaning on the normal p.d.f.\n",
        "\n",
        "As a formula the normal distribution function looks like (in one dimension)\n",
        "\n",
        "$$ f(x;\\mu,\\sigma) = \\frac{1}{\\sigma \\sqrt{2\\pi}} \\exp(\\frac{-(x-\\mu)^2}{2\\sigma^2}) $$\n",
        "\n",
        "It reads, given the distribution parameters mean $\\mu$ and standard deviation $\\sigma$, x follows this function.\n",
        "\n",
        "### Example 1.1 Plot the normal p.d.f\n",
        "\n",
        "Plot the normal distribution with mean 0 and standard deviation 5 for 400 x values between -20 to 20. Repeat this for two other means and standard deviations. How big is the surface under the curves ?\n",
        "\n",
        "(See also scipy.stat.norm https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.stats.norm.html)  \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1wxogyGodKUj"
      },
      "source": [
        "import numpy as np\n"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n7vfyQRiRUct"
      },
      "source": [
        "The normal pdf is only one of many built in pdfs and pmfs in Python (https://docs.scipy.org/doc/scipy/reference/stats.html). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u7KkgBZpdKUm"
      },
      "source": [
        "## 1.3 The PDF of your data - histogram\n",
        "\n",
        "Histograms can be used as approximate pdfs of the data. Let's look at the histograms of the Iris data.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1tVXLQabdKUg"
      },
      "source": [
        "# Load the Iris dataset into a dataframe an study the dataframe content\n",
        "# If you are not used to the Iris plant, look up on Wikipedia.\n",
        "import pandas as pd\n",
        "url = 'https://www.openml.org/data/get_csv/61/dataset_61_iris.arff'\n",
        "\n"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "id": "jtznAZBZdKUg"
      },
      "source": [
        ""
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F0yqLrBFdKUy"
      },
      "source": [
        ""
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eo8sm-kxUK0E"
      },
      "source": [
        "It is hard to see by eye if these data are normally distributed. There are statistical tests checking this better. More on that later."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WLOJ4il1dKUn"
      },
      "source": [
        "## 1.4 p-Value, Significance and Confidence Level\n",
        "\n",
        "The p-value of a (measured) value is the surface (integral) under the pdf above that value, i.e. the probability of getting such a value or higher. Thus, the p-value is a measure for significance. Most implementations of statistical tests return the p-value and the value of the test statistic."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yZNe1XCRQHkP"
      },
      "source": [
        "# For pdfs in scipy.stats, we can calculate p-values with cumulated pdf (cdf) \n",
        "# The p-value of 3 on a normal pdf with mean of 2 and SD of 5\n"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k4-Bj7yMRbeO"
      },
      "source": [
        "# On our data set we can calculate the p-value manually\n",
        "# The p-value of measuring a setosa slength longer than 5.2\n"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TARMMegtVSph"
      },
      "source": [
        "What is significant for you ?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0J94oO0ndKUu"
      },
      "source": [
        "## 1.5 Mean, variance, standard deviation, skew, kurtosis\n",
        "\n",
        "The $n^{th}$ moment of a random variable x with p.d.f. $f(x)$ is\n",
        "$$ \\alpha_n \\equiv E[x^n] = \\int_{-\\infty}^{\\infty} x^nf(x)dx   $$\n",
        "In the discrete case for $n=1$ this integral becomes the sum known as the arithemtic mean:\n",
        "$$ \\mu = \\frac{1}{n} \\sum_{i=1}^n x_i$$\n",
        "\n",
        "<!--\n",
        "and the $n^{th}$ central moment of x (or moment about the mean, $α_1$) is\n",
        "\n",
        "$$ m_n \\equiv E[(x-\\alpha_1)^n] = \\int_{-\\infty}^{\\infty} (x-\\alpha_1)^nf(x)dx   $$\n",
        "-->\n",
        "\n",
        "The most commonly used moments are the **mean $\\mu$ (or expectation value) and variance $\\sigma^2$**:\n",
        "\n",
        "$$\\mu \\equiv \\alpha_1 $$\n",
        "$$\\sigma^2 \\equiv V[x] \\equiv \\int_{-\\infty}^{\\infty} (x-\\mu)^2 f(x)dx = ... = \\alpha_2 - \\mu^2 $$\n",
        "\n",
        "The mean is the location or the “center of mass” of the p.d.f., and the variance is a measure of the square of its width. It is often convenient to use the **standard deviation (SD)** of $x$, $\\sigma$, defined as the square root of the variance. In the discrete case the variance becomes\n",
        "$$ \\sigma^2 = \\frac{1}{n} \\sum_{i=1}^n (x_i - \\mu)^2$$\n",
        "\n",
        "For the normal p.d.f the standard deviation is its width. \n",
        "\n",
        "Based on higher moments other distribution descriptors are formed. **Skewness** and **Kurtosis** you may encounter. The skewness is a number indicating the deviation from a symmetric form. Kurtosis is a number indicating if the tails of the distribution is larger or smaller then the tails of the normal distribution.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hXxW3QldYMch"
      },
      "source": [
        "# Get the descriptive statistics from a pdf in scipy.stats\n"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YQwudiIsYtJJ"
      },
      "source": [
        "# Print the summary statistics estimated from the setosa sample\n"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7EPOlNl-fiWg"
      },
      "source": [
        "Always remember this plot:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BMx1QRebe_wa"
      },
      "source": [
        "![Standard_deviation_diagram.svg.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAooAAAFFCAYAAABrIsxwAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3gc1dXH8e+9M7Or3qslWe69994AF2xsY8CYagzGmN5C6AECIaG3JPQSQgspQELemBpCSegJvTdTDBjjXqXdff+YBXdbWq80Wun3eZ59JO3enTks3pkzd+49F0RERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERESkwTkNtN3pQBrwdQNtX0RERERSVAy4JeggRERERCRxbgLvmVjHdpVbtF2QwL5EREREJCAmgffEGnFfIiIiIhKQRHoUAVYDDwCR7bw+G/gA+HeC2xcRERGRFHQ2sB54Hui2nTYaoygiIiLSQvUCXgXWAWexdc+kEkURERGRFswDLgBqgFfwk8cfKFEUEREREfoBbwIbgPPxE0gliiIiIiICQBj4FVALvI4SRREREZGUZ5O0nfXAGcBwIJSkbYqIiIhIM+MAWfi9jCIiIiIiIiIiItKcNPRqKRlATvz3rxt4XyIiIiKSQk7Cn9iS6LJ/IiIiIhKQRJfwq6u1wJJ6tD8evxeyLv4IfFzviEREGkksFmsFlANvGGM2BB2PiEh9NfSt5/r6E5Bbh3aF+PUaezZsOCIidWaAtvi1ZQdhner27duNyMvJafXm2+++t75mw2JqI+9C9CPgRfxlUNcEGbCIyM40tUSxrmYA1wEVQQciIi1aObA/jjuBWGwU0Wg6GXnrTE6xa7IKvN2H9mefEb04+75n+H7Jd7Dyu2hs1XfrYku/8ojUGBzvFWo33IJ/kbws4P8WEZGtJOvWc4iNPYHL8VdpERFpjsLAQTjeXCK1Q8jMX2+reoZNeWdjituAG07/oeHI6cOYN2MIl7xpWLZkJfi1azOIxYgt+4rYwtcGRz96sQ/rV1+PMX8gGrkSeDWY/ywRka3tSqI4DJgHjAFas7F3MgYsBJ4CbgL+vQv7EBFpKjKAY7Hu6Thulu04NGSq+xmT3yqt3lsyBpNfgcmvwPaaFI4teo/ou0/tF/v6gwNxnL8SiZwOvJf0/wIRkXpKJFE0wDX4E08AosBX+D2J4PcsVgGz44/rgBPRzGcRSU0OcDjWvZhwRrbtNTFs2/QHJ0k3ZIzBtOqC06pLKLZsEdH/PTwptui9vTD2DqKRnwBLk7MjEZH6S+RIdxR+kvgacD7wGLB6izaZwHjggnjbt4AbE45SRCQY3bDuXVjbw/aa5NlOw8A2XLEIk1eOM+bIUOzbj4i+cP+hsTXLZhCpPRJ/DKOISKNLZK3no4FPgBHAg2ydJBJ/7gH8tZ8/jb9HRCSVnIgxr5mKrj3cvc70bJdRDZokbsqUtMeZ/FPXdt89D+v8ATf8CJq8JyIBSCRR7AQ8DKyqQ9uV8badE9iPiEgQMnG8+3G8y51hB7nOyMM80nN2/q5ksw62xx44k39qTUHlWKzzATCz8QMRkZYskURxLVBcj/bFqFaYiKSGMqz7b9JypjqTTnFNdd+g48FkF+HsfrRn++6VjnXuwfFuAeo/gUZEJAGJJIrP4Ncx3KMObScAe8ffIyLSlHXFcf9rCiq6uRNPDJuckqDj2YTBdh6JM+Ekh1D6wVj330CroKMSkeYvkQE35+IniY8AjwILgA/YfNZzJ2BivN26+HtERJqq0Vjn76aiW7oz9CCbtBnNSWbyW+FOOjUceeqW7rFlX71CNDIS+DDouESk+UrkaPg6MA64Fb/HcMIO2r4FHAG8kcB+REQaw74Ye6/tPMravpNtk1+wKi0LZ/xxocizdxXFvnrnBaK1o/CPtSIiSZfoZfPz+OssD8cvuN2JzVdmeR+/4PZz+HUWRUSaopkYc6/tPcnabuOCjqXurIsz8lA38u97cmKfv/5vopFxwCtBhyUizc+u3F+J4o891PhDEUlF+2PMPbbvVGu7jAo6lvozFmfYgW70BTcr+umrTxON7IFWwhKRJEtkMouISKqbhTH32P7TUzNJ/IGx2CH7W9t+cBrWWQB0DzokEWlelCiKSEuzD8beZQfua22nEUHHkgQGO3CGNdV9MnHcp4DqoCMSkeZDiaKItCR7Yux9ztBZju0wJOhYksjgDJllTXnXPKz7HFAedEQi0jwoURSRlmIg1vmz7TXBMW36Bx1L8sXHLJq8slJc71EgPeiQRCT1KVEUkZagK9Z9wnYcHrbdd2/i9W92gRvCGTvPJS2nE457F02+1o+INHVKFEWkuSvCOo+bqp7ptv/U5p84hTNxxh0Vwrp7Ye2FQYcjIqlNiaKINGdhHHeBKaoucoYe4LaUDjaTVYgz9kgPOAM4Muh4RCR1KVEUkebLOtfhpfd0RhwawjpBR9OoTFEbbP+9HYz9DdCcZu6ISCNSoigizdV8MHOcsUeGSMsOOpZA2I7DsB2HW6zzV6As6HhEJPUoURSR5mgIxvzaGXaQa/Irgo4lULb/VMcUtsnD9R4CvKDjEZHUokRRRJqbIqz7kO00wpjWvYKOJXjG4oye4xHK6IPj/CbocEQktShRFJHmxXFvM9lFebbPFB3ffhBKxxl5WIgYhwMHBh2OiKQOHUhFpDk5GphkR84O4bhBx9KkmMLW2P7THay9FegcdDwikhqUKIpIczEYY65zhh7ompySoGNpkmzHYZiq3h6O93e0couI1IESRRFpDnKw7p9su0GY1r2DjqVJcwbt6xDOqMI61wQdi4g0fUoURST1Oe4dJquw2A6Y0bKKJSbCS8MZOSdELHY4sH/Q4YhI06ZEUURS3SHEYlPt8IPDGpdYN6awCttzD4t1bgJ0n15EtkuJooikso5Ye5PtP80x+a2CjiWl2O67G5NXkYbj3hJ0LCLSdClRFJFU5WC9P5jKHo7tODzoWFKPsdgRh4SAScAhQYcjIk2TEkURSVH2FFy3uzNwX602kiCTVYDtM8XFOjcAlUHHIyJNjxJFEUlF3TBc7AzcN0Q4M+hYUprtNAJTUOXqFrSIbIsSRRFJNQ7W+72p7B4z1X2CjiX1GYMddlAI2A3dghaRLShRFJFUcyKu28MZtJ9uOSeJbkGLyPYoURSRVNIRY3/pDJihW85JplvQIrItShRFJFUYXO8OU9bJmjb9go6l+TEGO2T/ELHY7sCsoMMRkaZBiaKIpIr5YAY6g2eqqnYDMTkl2F4TLda5HsgPOh4RCZ4SRRFJBRVYe4XtP90jIzfoWJo122W0MZn5GVjnV0HHIiLBU6IoIk2f495gito4tv2goCNp/qyDHXZwiFh0LjA06HBEJFhKFEWkqZtMNLqnHbhvCEzQsbQIprAK225gFOveDuhWv0gLlmiimA4cBpwODNzk+S7AncCrwJPAiYBKWIhIotKw7g2262hjckuDjqVFsX33cnHctsAxQcciIsFJ5EoxE3gG6Bv/O4afNL4EPA/kbNJ2LDABmAJEE45SRFqqMwhnFNse49WV2NhCGdh+00LRl/74S6LRPwFfBR2SiDS+RHoU5+Inia8CZwDPApcB58W3dzhQjZ8kvoW/4PzByQhWRFqUThhztjNwnzBuKOhYWiTbfiCmsLWDda4KOhYRCUYiieJBwCpgHHAJMD6+nf2Bs4DbgYXAU8BU/J5EJYoiUj+Oe4Mp7xw1lT2CjqQFM9hB+4WJRffFP+aLSAuTSKLYFr8XcXn873XA0/HfH9yi7cfAy0DvhKITkZZqCrHYKNt/b3UlBszklmE7johhvevRmHORFieRRDEb+H6L576L//xyG+2/QIVbRaTu0nHcG2y33azJLgo6FgFs74kOrlcNHBt0LCLSuBJJFL8DCrZ4Lhb/ua0JKzls7H0UEdkJezqhjELbfTdNYGkqvDRsv73CWOcioDjocESk8SSSKL4HdNjiuZPZfq9ha/wxiyIiO1ON4SxnwIw0HJXva0psu4GY/AoH61wedCwi0ngSSRRfwE8UN72qXAss20bbKqAj8K8E9iMiLY11LjPFbaOmqmfQkchWDHbgPmnEogcDWiJHpIVIJFG8ECgHltahbSVwMXBHAvsRkZZlBLHovrb/tHDQgci2mYJKbJv+Uax3A1omR6RFSCRRXAt8DdTWoe1/gHOA1xPYj4i0HAbrXWfbDIia/IqgY5EdsH2nuBh64JdKE5FmTms9i0hTMAtDd9t3shN0ILITadnYnuNdrHslkBV0OCLSsBo6UcwAyuIPEZFtycC619juu7mkZQcdi9SB7TLakJ6dg7/Igog0Yw2dKM4DFsUfIiLbYE8hnJ5ju47RmLdUYR2cAXuHMeY0/AmLItJMNXSiuBZYEn+IiGypBGvOcnpPDuNo0Y9UYiq6Y4rbRjDOL4OORUQaTkMXKrsx/qirx6jbKi55gJb2Ekl11rnQZJc4pu2AoCORBNj+08ORBVfNAIYAzwcdj4gkX1OraPtX/HGNO9MD2KOBYxGRhtWdaHSuHTjDYnTXORWZ/Aps20HR6Kev3kC0pi8bV+kSkWYiWYliCMiN/74c2JDgdq6rY7sZwLgE9yEiTYHjXmHKOkZMSTtVX0hhts+eTvSzV7sBs4B7g45HRJJrVw7Qw/ALaX8KrAO+jT/WxZ+7I95GRGRL44hG9rB9p2pgYqpLy8J2H+dinaup2x0hEUkhiSSKBrgWeA6Yjb9M31fA2/HHV/HnZsfbXIsq+IvIRgbXu9K2G4jJKQk6FkkC23WsIZSRC/aEoGMRkeRKJFE8CjgeeA3YG8jBX6qve/xRGX9uBvBGvO28ZAQrIs3CvsRi3W2vSbrl3Fw4Hk6fyWEMPwNKgw5HRJInkQP10cAnwAjgQWD1NtqsBh4AhuPfhj46wfhEpHkJYd2rbJcxLuk5QcciSWTaDsDkllqsc0HQsYhI8iSSKHYCHgZW1aHtynjbzgnsR0Sanzk4TqntOiboOCTZjMH23StMLDoXFeEWaTYSSRTXAsX1aF8MrElgPyLSvGRh3V/anhNcQulBxyINwJR3xpR2iOC4VwQdi4gkRyKJ4jP44w/rUsdwAv44xmcS2I+INCv2dNKyMm2nEUEHIg3IDpgRIhqZDIwJOhYR2XWJ1FE8Fz9JfAR4FFgAfIBfPxH8eoqdgInxduvi7xGRlqsSw0+dvnuFsE7QsUgDMjkl2DYDYtHP/ncl0Zr+qAi3SEpLJFF8Hb/Y9a34PYYTdtD2LeAI/NnPItJSOc4FJrcVprp30JFII7C9JznRz17tCUzHn9goIikq0ZVZngd64s9qHoPfg7jpyizvA0/h11GM7lKEIpLquhOJHmb7TbUqqdpCZORiO490ou89dyXRmr8BtUGHJCKJ2ZUl/KL4Yw81/lBEts9xLzclHaJaqq9lsd13N9EP/9OKKHOAm4OOR0QSowO3iDSksUQj423/qclaV15SRSgd22OPENb9JZAZdDgikhgliiLSUAyud5VtOwCTo8U6WiLbaSSE0rPBnhR0LCKSGCWKItJQphKN9rA9J+g401I5LrbXxBDWnAnkBx2OiNSfDuAi0hBcHO8q22mkJVP5QUtm2w+CjDwX7BlBxyIi9adEUUQawuFYW2l77K5pzi2dsTh9JocxnAS0CjocEakfJYoikmwZWPdi230PT0v1CYBp3QuTWwbWOS/oWESkfpQoikiS2VMIpWfbziODDkSaDIMdsHeIaHQu0DXoaESk7pQoiki9xGKxAbFY7PtYLDZnGy8XYs2ZtveeIRxVxJGNTEk7TEn7CI570bZej8ViT8ZisecaOy4R2TEliiJSXy7+DNbw1i/ZM8nMd227AY0dk6QA22+KRySyN9BvGy9nAzmNHJKI7IQSRRFJlioMJzh9poQwOrTI1kxBFaaiewTHvTToWESkbnQ0F5HksM75Jq9V1FT1DDoSacJs38ku0chYYEzQsYjIzilRFJFk6Eo0epjtP3Ubt6NFNjI5JZjq/jEc77KgYxGRnVOiKCK7znF/YUraR0xJ+6AjkRTg9J7oEI30BSYHHYuI7JgSRRHZVYOJRKbbfnt5QQciKSIzH9txqMV6V6LzkEiTpi+oiOwax/2VqewRMQWVQUciKcT2HG8wtANmBh2LiGyfEkUR2RW7EY2Msn0nq2ii1E84C9tllIvjXQKoN1qkiVKiKCIJiUajButdadsNxGQXBx2OpCDbbRxYWw5sq3i7iDQBShRFJCFnnXVWb4h2s70m6TgiifHSsN1397DuxatWrdK/I5EmSF9MEam32tpabrr5ln1tpxEO6VpMQxJnO48EL5x96aWXqltapAlSoigi9fb73/+e1WvW5Nnuu5mgY5EU53jYXhNC1157XaulS5fqnCTSxOhLKSL18vnnn3vnn38+fSbOtIQzgw5HmgHbfgjZubn2oosuKgo6FhHZnBJFEamXWbNmTVuzZg29dp+h3kRJDutw+LEnm5tvvrkIKAw6HBHZSImiiNRHzkuvvHLUz372M0LpGUHHIs3IlBkzqaxuY7DOeUHHIiIbKVEUkfo43kvPDh155JFBxyHNjLWW4087C0NsPqDq7SJNhBJFEamrUmOdc3ofcEIoLS0t6FikGdpj8l6kt+4as473y6BjERGfEkURqRtrzw4VV5o2wycGHYk0Y6WT54Wi0doDge5BxyIiShRFpG7aEePo0snzwsbosCENJ6NtT7I6D4gYN3Rp0LGIiBJFEakD63gXpFd3jWR3Gxp0KNIClE4+yotFaiYBw4KORaSlU6IoIjvTOxqpPahs8rxw0IFIy5BW3o6cPmOjxgtdGXQsIi2dEkUR2SHjhi7O6jKgNqNtz6BDkRakdNIRDpHIQECDYkUCtCuJYiWwHzAVyNtBu4nAT3ZhPyISnGGxSM2k0slHeUEHIi1LqKCcvMGTMW7oKtSpIRKYRL98ZwCfAPcDDwGfAydtp+0s4LIE9yMiATKud2Vun3GRtPJ2QYciLVDJHodYiHUAZgYdi0hLlUiiuCfwSyAG/Am4DdgAXAX8LsFtikjTsyfR6IDSSUe4QQciLZObU0jR6P0c44UuAdSrLRKARJK6E+I/J+Pfej4C6AA8CBwK3J7gdkWk6bDGC12ZP2Qv4xWUBR2LtGCFYw8wxrpl+OcaEWlkiSR0A4H/AI9t8txSYB/gRpQsijQHs4B2xeMP1fdYAuWkZVK8+8Eh43gXAVpgXKSRJXISyAI+3cbzUeBoNiaLtyW4fREJVsg4oUuLRs/03KwdzVMTaRyFI2fgpGdmsfGOlog0kkQSua+A0u28FsNPFm8HZgM3Ayax0EQkIHON65YUjtk/6DhEADBuiJKJh4eN45wDFAQdj0hLkkii+DowmO3fAogBc4F7gMPRbDWRVJJpXO/nJeMP9Zy0zKBjEflR3qA98fJKXeCnQcci0pIkkij+BcgEdtTdEMW//fxHIC2BfYhIEKw9yUnPzi4YPj3oSEQ2Y6ylZM+5YWPtyfh1fEWkESSSKD4InAms2Um7CHAg8HPgNwnsR0QaV6Ex9qySSYeHjBsKOhaRreT2Gk24vF0U65wTdCwiLUUiieJy4FfAH+rQthY4Dzgugf2ISGOy9jQ3p8jN6z8+6EhEts0YSicekUYsdgSgKvAijUCzkkUEoMrAKWV7zQ8ZR/W1penK6jqYzPa9o8Z1rww6FpGWoKETxQygLP4QkabKcc4JlbaJ5fQcGXQkIjtVOuWoUCwSmYpf11dEGlBDJ4rzgEXxh4g0TZ2JRueWTz0mhFE1K2n60is7kdtnbNR4oWuCjkWkuWvoe0xrgSX1aF8NOHVoV4rqM4okhXG96zI79o9kduynoSiSMkqnzHdWvP7MQGBv4IGg4xFprho6Ubwx/qirl4DiOrb9vv7hiMgWdo9FI7uVTT1aSaKkFC+3iMJR+zpLnv3zVbGaDQ8DNUHHJNIcNbVR652oW0/hXvgzr0UkcdZ4oevy+k8gXFzVKDtMD7mcNX0wgzuU07E8n+LsDFas3cCni5dz17Nvc8dTb7JmQ+0Ot3HGtMFcOHM4AFMu/QuPvPZpnfZ9wLAujOxSSd+2JfSsKibsOZz9h2e59K8vbrP93HG9OHnP/pTnZ/LGwu849fdP8fLHX2/V7uAR3bhp3ngGnnUXb33xXZ1i2RF9RnVXvPtBZunzfyuL1GyYC1yftA2LyI+S1YsQwu8JLI7/nqhlwNI6PFbhrwAjIok70GDal0yc02i9ibkZYc6YNph2JXm8++X3PPDSB7z6yTd0rSjkmtnj+OfP9ifN2/71a4+qIs6dMYTV6+vfeXT+vsM4crdetCvJ4+tlq3fYdtawLlx/xO48+NIHzLrmYdZuqOWxs/ejPG/z1WqKstO57ODRXPa3l5KWAOkzqjsbzqBkwpywcdyLgZykblxEgF3rURyGP1llDNCajT2BMWAh8BRwE/DvXdiHiDSMDON6VxSNO9Bzs/IabadLVq6l8pgb+Gb55vX6czPCPHb2fvRrW8rMIZ2485m3t3qv51hunz+R1z5bzPuLlnLQiK712vfRtz7OR98u47PFKzhhYj+uOGTMdtvO260XL3/8NWf/4VkA/vfpt3z+26M4ZFT3zXrXLj94DEtXrePiB5+vVyw7os+ofvKH7sV3/7o/vWbpN6fi1+0VkSRKpCfBANcCzwGzgSrgK+Dt+OOr+HOz422uRRNPRJoWa4930rLyCsc07lLsNZHoVgkQwPI16/nj8+8BUJq37TWmz5o+hK4VhRxx4yNEYvW/ofDkWwv5bPGKOrVtU5zL219snIf37Yo1LF6xlrbFGzutdu9ZzYHDuzL/1sdYXxOpdzzbo8+ofozjUrbX/LCxzulAeYPsRKQFSyRRPAo4HngNf7ZZDv66m93jj8r4czOAN+Jt5yUjWBFJiiJj7Lkle84NWS8cdCwAhD2HSX3aEovBP9/6fKvX+7Yp4Yxpg7joged558v6FFJIzDfLV1O2STKW5rnkZ4b5epmfwGWEXH57+O7c8a83efqdLxo8HtBntCM5vUaTVtEBHOfCBt2RSAuUyK3no4FPgBH4YwW3ZTV+uYLHgdfj76nP7GcRaSjWOS9UXOnmD5wYWAjpIZcrDh6DMZCflcbILpUAHHHjgq0mQ4Rch9uPnsgbn3/H5Q+/1Cjx/eXFD7hw5gim9m/Ps+99yVnTB+M6lgdf/gCAn+0zjIywx+n3PN1gMegzqp/y6ceHP/718XOAq4C3GmWnIi1AIoliJ+Bmtp8kbmol8DAwN4H9iEjydSQWO7psynwnyOLaYdfhyN16bfbc/c+/x4sfbT1j9rx9htKpvIAh59xNbSTaKPFdu+BVOpXn8+dTpgGwen0N8295jNc+W0zv6mJOnNSP2b/9ByvWrudXB4xizpgeZKV5/POtzzn61sf4fMnKXY5Bn1H9pFd3I7vb0Miq916+PFa7YVJSNy7SgiWSKK6l7rUOibfdesCNiDQ647qXZbTpEcnqMqguhe0bzLI16/EOuhJrDBUFWUwb0IELZ45gYu+2DPvZPbz3lV8mdVCHck6dMpCLH3ye1xcubrT41tdEOPKmRzn1909RlpfJJ98upyYSxbGGG+eO59HXP+X+59/j1CkDOGnP/pxx79N8sGgp18wex70nTGHEeffucgz6jOqvbMp874N3DhsPjAOeTPoORFqgRMYoPoM//nCPOrSdgD+O8ZkE9iMiyTU8FolMLZt67K6UsEqqaCzG50tW8utH/stJdz5JTnqIM6YNAsCxhtvnT+SdL5fwywdfCCS+FWs38P6ipdTEe+mOn9iPzq0KOP72J378+88vvs/V//cKf//vx5xx79MM7lDO0E6tkhaDPqO6CxVXkj94CsYLX4kmUYokRSI9iufiJ4mPAI8CC4APgOXx13Pxb09PjLdbF3+PiATHGC90bU6v0dG0Vu0D7U3cnn++uRCAHpVFgD9Gr1N5PgBr7jxpm+95+KczANj3qr/y0MsfNmh81cU5nL/vMM69/zkWLllJZtijIj+LVz7+5sc2L8d/71xewH/e/yrpMegz2rmSCYfZ5S8/0i0G+wH3J30HIi1MIoni6/jd+rfi9xhO2EHbt4Aj8Gc/i0hwZhCN9Srd88gmmSQCtC3x6zkuW7Me8MvE3P7Um9tsO7JLBR3K8nn09U/58vtVfLp4+TbbJdNv5uzO218s4TeP/Hez5xtzqKc+o51zs/IoGnegu/iJu6+K1W54CFjfeHsXaX4SLbj9PNATGI5fcLsTfk8i+D2L7+MX3H4OaJyR1SKyPSHjhK4uHLWP4+UWBRrI8M4VLF6xhvcXLd3s+bYluVx72DgA/vrKR4A/Bm7ezY9uczu3zp9Ih7J8rl3w6lbL0x25Wy/SPZdbnnx9p0vd1dUBw7qwW4/WDDrnbqLx+oSr19fw+ZKVDGhX9mO7ge1LAXj3q8TL0+gz2nWFY2aa7597oLB21Yaj8Gv5ikiCdmVllij+2EONPxRp2g43rltSNG5W4GO2JvRqw5nTB/Pul9+zcMkKVq2robIgi/7tynCsYcFrn3DDY//bpX1cOHMEhVlp3P/8e5slQQeN6Mq4HtUAdKsoAGDfwZ3o3Mr//el3Pud3/9q6qkpBVhqXHzKGK//+Mm9sMVnkugWv8ssDRnHqlAF8+PUyfnXgaF74cBHPf7Ao4fj1Ge0664UpmTAnvOihX18Qq625Hb8Ch4gkYFcSRRFp+rKN411cMnFOyEnPDjoW/vjC+2SnhxjdtZKB7cvITguxbM16nnr7c+559m3uevadH3ujkm1wh3IOHdlts+f6timhb5sSACKR6DaToMsOGs2qtRu46IGtl6C7dsGrFGWnc+rkgWSnezz55kKOve2JXYpTn1Fy5A3ek8VP3ptes/zb84lGT23wHYo0U4H3MCRoBnAdUBF0ICJN3K9Cha1O7HD679KMTc7wxN7lOdx9YH+Ove1xbnri9aRss7k6b59hnDNjCO1PuJmFSa4b2Nw8f+FB5OZmxabf8WLSzkur3nmBz247K0os1hvY9mBOEdmhRMrjiEhqqDDWOY43EnAAACAASURBVKlkz7lJSxJFUklW18Fktu8dMW7o8qBjEUlVShRFminjuNekt+1Jbu8xQYciEpjyGSd5sWjtHsCeQccikoqUKIo0T6Nj0ejerWacGA46EJEghUtaUzh8BsYL/Rrwgo5HJNUoURRpfhzjhn5bMHRKLFxaHXQsIoErmXCYtV64Ajgm6FhEUo0SRZHmZ65x3A4lk47QwEQRwKZlUDppbsg47kVAsMVERVKMEkWR5qXQOO5lpXvObRLlcESaivwhUwiVVLlY58KgYxFJJUoURZoTa8/zClt5+UOnBh2JSNNiDK32OTmNWHQe/spiIlIHShRFmo9OxDimfNqxacbqqy2ypYw2PcjuOrjWuN7VQccikip0NhFpJozr3ZzdbUg0q/PAoEMRabLK9zk5BGYksH/QsYikAiWKIs3DJKLR4WV7Ha3yHyI74OUWUzzuQNe43jVARtDxiDR1ShRFUl+accM3FY07wIaKtKqlyM4UjZtl3Mz8fKz9WdCxiDR1ShRFUp2157k5BUXFexyaqmu3izQq44Zotf9pIWKx09DEFpEdUqIoktq6EuMnFfudkmYcN+hYRFJGVqf+ZHcbGjVu6GZAF1ki26FEUSR1GeOGbsvpOSKW2bF/0LGIpJzyvU9wIdYfmBV0LCJNlRJFkdR1kLG2X/neJ2gCi0gCvLwSinc7yDGO+2sgN+h4RJoiJYoiqSnHuN7VxeNne252QdCxiKSswjH7GzcrLwtrzwo6FpGmSImiSCqyzkVeQVlW4ch9NLZKZBdYL0z5jBNDxDgF6BZ0PCJNjRJFkdQzmFj02MoDzgprAovIrsvuPpycXiON8UL3Ak7Q8Yg0JUoURVKLNW7ot3n9d4+kV3UOOhaRZqN82nGOMbYbcGTQsYg0JUoURVKJtScYx+1ZOvkoTWARSSI3p5DSKfNd47hXAK2CjkekqVCiKJI6OhrMJa32PVkTWEQaQMHQvUiv7uYa17sj6FhEmgoliiKpwRgvdGdWl0Emt+9uQcci0jwZQ8XMn4SIxcYB+wUdjkhToERRJDUcZYzt32rfU3TLWaQBhYoqKd79YGsc7wYgP+h4RIKmRFGk6asyjnNl2dRjPDenMOhYRJq9onEHGi+/NAvHuSzoWESCpkRRpIkzbuiGtMpOTv7gyUGHItIiGMel1X6nhohGDweGBB2PSJCUKIo0bTOJRSdUzDojhFFtbZHGktm+N/mDJ0eNF7oLSAs6HpGgKFEUabpKjePeXLrnXCdcXBV0LCItTtm0Yx03q6AS61wadCwiQVGiKNJEWce7JlRcGS4YuU/QoYi0SNYLUzHrp2FisWOBoUHHIxKExkgU5wC3NMJ+RJqT/WPE9q06+GdhY7WimEhQMtv3oXDkPhg39CcgO+h4RBpbYySKo4EjGmE/Is1FmXHcm0smHOaEy9oEHYtIi1e651zr5ZcU4ni6BS0tjm49izQxxg3dkVbZKVw09oCgQxERwLgeVQefGyYamQdMCDoekcbkJvCeM+rZvmcC+xBpqeYAu1UeeLarWc4iTUdaRUeKxuxvljzz5ztitRs6AyuCjkmkMSSSKP4y6VGICEBr4zi/KdtrvhsqLA86FhHZQsmEw8yKt57Nr13y9dXRSM3hQccj0hgSSRTXA4uAa+rY/kBgYAL7EWlJjPHCt6VXdnQLhk0LOhYR2QbjelTOOiP88XXHzwb+BPxf0DGJNLREEsU3gLbA1XVs3wcliiI7c4SB0RUHnKlbziJNWHrrrpTscahZ/MRd98Qitd2Ar4KOSaQhJTKZ5VWgEKhOciwiLVW1cdxry6Yf54YKdMtZpKkr2v1gk966S4Zxw/ejSaHSzCXyD/xp4COgQx3bPw78JoH9iLQEjvFC92e06+XmD9oz6FhEpA6MtVQefK5nXGcQ1p4VdDwiDSmRRPFu/CTxiTq2vws4LoH9iDR/1r3YhtL7Vh50jqdbziKpw8stpurgn3nEYhcA44KOR6ShqMtcJDh7EIv+pOrQ8z03Ky/oWESknrK6DKJwxD4Yx70PKAo6HpGGoERRJBhlxnHvL97tIJPZvnfQsYhIgkqnzLPhsjbZxnHvBnRbQJqdhk4UM4Cy+ENEfMZ44XvTq7tllEw4TCcWkRRmHJeqQ89PM9YZCxwbdDwiydbQieI8/JqLixp4PyKpw9rjjbHDK2edEdK4RJHUFypsRdnex3vGOlcA3YOORySZEqmjWB9rgSX1aH88fi/kzvQAnIQiEglWb4O5omLWT12vQB3tIs1F/qA9WfXeS2bl2/95IFazoT+wMuiYRJKhoRPFG+OPupoK5NehXR5KFCX1ZJpQ+KG8fuNtTq/RQcciIklWsf/p3sfXzG+9YcnXd8dqN0wDYkHHJLKrGjpRrK896thuBnBdQwYikmTGuO59Xm5xWdm0YzSJTKQZsqE0Wh9+cfijK+dNNNHac6PR6M+DjklkVyXrhBUCiuOPUJK2KdJ8WPsL44YntDny0rD1wkFHIyINJFTYitZzLvSiMc4D9gs6HpFdtSuJ4jDgDuBTYB3wbfyxLv7cHfE2Ii3d/sQ4verQ8z2NSxRp/jI79KVs8jxjrHMn/ph6kZSVSKJogGuB54DZQBX+ouhvxx9fxZ+bHW9zLaotJS1XL2Od35XuOddkdeofdCwi0kgKx8w02T1HuMYJPYw/rl4kJSWSKB6FPzv5NWBvIAeoxC8J0D3+ew7+OMI34m3nJSNYkRRTZJzQguxeo9yisbN0sSTSwlTMOsMNlVSUGcf9M5qAKSkqkUTxaOATYATwILB6G21WAw8Aw/FvQx+dYHwiqco1rveXcGllYcX+P9UJQqQFsl6Y6jm/CBsvPBJrNbFFUlIiiWIn4GFgVR3aroy37ZzAfkRSl7UXGccdWnXoBSFNXhFpubyCMioPPMsjFjsTmBZ0PCL1lUiiuBZ/dnNdFQNrEtiPSKqaSYzTqg493w0VVQQdi4gELLvbUMqmHmOMde4HhgQdj0h9JJIoPoM//rAuNQ8n4I9jfCaB/Yikov7GOr8rmzzPZHUeGHQsItJEFI7ch4Lh0z3jeo8AHYOOR6SuEim4fS5+kvgI8CiwAPgAWB5/PRf/9vTEeLt18feINGmxWOw24F1jzKUJbqKTcdwnC4ZNDxWOmanJKyKymbKpx5jaNSsyV/zvqWdjkZoBwOf13UYsFssFbgb+aoy5K+lBimwhkR7F14FxwDv4PYZX4Y9DfCb+eBi4EhiPXy5nLP7sZ2kZcoCDgb/gl0qqAb4D/g//30R9HAjcALyIf8ERA87YQft5wLvACuBZYHtdeocAG9i6vtk0/AlYiSg1TuiJnJ4j08umHq2VV0Rka8ZQMfM0J6Ntj3zjhZ4gsbI5YfxC3j13IZJUPk5LI0v0hPY8/j/SUcDPgLuAv8Ufd8WfGwX0Al7Y9TAlhRwD/B6/R/k94E/A+8Ak/F7o0+qxrZ/jl2NqDyzaSdsD8NcVfwD/ILoWeBIo36JdEf6FzCXAm/WIZUfCxkv7a7i8uqTV/j/1MOpMFJFtM45L1ezzPS+3pNp44fsJZindlniclgTtSs9HFL8H8UL8zH9q/HFI/Lln4m2kZfkEmIs/iWks/oFhGP5Vag1wMX6tzbo4CmgLFALX7KTtfOAl4Ez8A91BQAZ+4fdNXQV8D1xUxxh2xjFu6I9uVm7v6rmXaIaziOyUk55Nm6MuC9lQ2hjreDcGEEJLO07LLtAtMkm2PwC3snV9zceAh/CvnofWcVtP4NfhrIu2wFub/P0tsDj+/A/2wD8wzQPW13G7O+IYx/2Lk5E7vu1x14XdLC2+ICJ14+WX0vaYqz1c9xCs/UUj774lHadlFylRlMa0Iv6zIb78X7P57Ys0ID/+PPhXrTcAtwH/SsoeHe9qG0qb1Gb+ZWEvtygpmxSRliNcWk3V7J97BnM6cHjQ8cQ1r+O07DIlitJY8vAni6wCnm6A7f8Z2D2+jwL8Wycu/lgYgPOBTOo39mb7rL3EOvaoNvOv8MIlrZOySRFpebI69afykJ85xjo34U8MCVLzOk5LUihRlMZg8a8QC/FLJS1rgH1cDfwOf1nJJfi3LeYB/wP6ACcDJ+JfLV8ab7Me+AdQVa89Wfcq43inVM+7zEurUDk0Edk1OT1HUnHAmY6xzp3APgGF0byO05I0ShSloRngt/iF13+Pf6BoCOuBI/DreHbBv51xK+Dg1xx7BH9czinxxy/wD8hdgPvrvBfrXmUd59g2R13uZrRR1QYRSY7cvuNotd+pjrHOfcBejbz75nWclqQKYlq+tBwGuB5/VtzdwJxG2OcKNo6xAf/qtAsbr9JPxC8FcWX873T8A9CwnW7Zuldbxzmmev4VXkZ1t6QFLCICkDdwItGa9c6iB659kFjsQPykqaGl2nH6340Qn2xCPYrSUCxwCxsPPrOBSCPH0Aa/xtc5wEL8sS8VwMubtPnh9y472I7FOr9VkigiDa1g2DRTNmW+xdh78Gf/NqTmdJyWBqIeRWkIDv44lIOAO/GvUIOoqXk9fimG67Z4vj4VsUPGde8xjjeteu6vXCWJItLQCkfvh/FCdtED191JLOriH0+TrTkdp6UBKVGUZHOBe/Cr7t+OX9R1Zwefo/BvLdwErElSHAfiz67rv8n+V+OvrbrpklE//P7ONraRZVzvYRvOGNZm/hVuWnm7JIUmIrJjBcOmYdyQ/epPV9yGMS6RyK1J3HxzOk5LA1OiKMl2Lv7BZx3+Fett22hzH7Bgk79/gT/T7j42PwAdjH8QAfihK28/Nt5++Bf+QW5LBfiV/S/HX5t8U9fgLwt1GvABcBn+kpT/2aJdsXG9x9zsoq5tj77C8wrKtrEbEZGGkz9oEjYUtl/e88sbjbVF0Wj0kiRturkcp6URKFFsGD3xlzEciV9A9D387vWb8BdMr4sD8dfL7oe/ZnYYf9mjX22n/Tz8WWKt8L90J+MvlbSlQ/BnmfWjYdbQLIz/TAMO3U6bN9n8ALQ9Q9h6aad+8QdALds+AF0BrMQf97Klq/HXET0NyMZfVWD+Fm06Gzf0eFqrdqWtD7/Y04orIhKU3D7jcDPznIW3n3ORNbHWixcvvrC4uHhXN9scjtPJ0JLP1XWmySzJNxh4AX9x9cfwBwrn4leb39k6mJtK1YXWj8MfW7Kjx+XbiMuwsTp/Xbc1dzsxzAE64H8OW4rgf4lL8G+jTAG++OHFO++8s9RY5395gyaUtz3uOiWJIhK4zI79aHfi9a4NZR7ZtmOne9eu3dahrV5S+jidJC39XF1nShSTy+JfAaQBU4FZwPH4VxlvxH/feRkWnxZab2S33XZb+MgjjxxUvPsh4VZ7n+QY6wQdkogI4C/31/a467yImz5i3LhxvPPOO5lBx5TCdK6uByWKyTUc6I5/dfLIJs+vBM6L/z6vjtvSQuuNx7OOc8P8+fMzT7/sN6Z4/KEGowl3ItK0hArL6f2TG9xYLMaIkSPn4K9mIvWnc3U9KFFMrjHxn9sa1/Eo/qyuMdt4bVc1h4XWR+N/MRo7Qyv3QuEXMnIK5zz93H+iYydPb+Tdi4jUXSi3kGeffZbRE6enG2tfwh/L1tjGAkcGsN9kGRP/qXN1HWgyS3J1iv98fxuvrcYfv9AaCAEbkrjfP+PPSJsGPINfuDTVFlo/FDgcf5xIXQcR76o+Xlr6I4UV7Qunn3al07tPn+jzC5fWeyOTu5bSvzKPbiVZdCrOIuRYrn7mY2558bOt2oZdy/whbehVnkOb/AzyMzxWb4jwxfK1/O3tb/jLm4tYV1O/ere7dyxm/94VdCvNIs11+HbVel79cjnXPfcxX6/ceEG6X69WzO5fRXFWmPcXr+KSpz7kza9XbLW9vbqVceGELuz3+5f44LvV9f48RBpLS/3uua7LCedcbL43Oe7T91xzRwy6EY2eiz9xpDEchj/b+eZG2l+y6VxdD+pRTK7c+M/l23l9OX6PWe52Xk+UFlqvH2OtPdk6zktdh08qOeCCW52s/MRnER43rC0ze7WiKi+dxat3fEzJDrscObiayrx0Pv5+NY9/sJi3vllJ+8JMzhrXkTv370vYrdvX0hrDBeO7cPXUHvQsz+b1RSt48sPv+GbVeiZ0KqYqL/3Htnt2KeW8PTrz+IeLOeVvb7KuNsptM/tQnBnabJv56R4/HdOBW1/8TEmiNHkt/bs3cK9D2fes39q0jOxTvVD6c/gzaWXndK6uB/UoNg8/LLR+Mn639sdADVsvtH4a/rT8n+JfSV2Hv37m0MYPOTAFbjjtbmKxPcbPO9fpNnLyLm/wgsfe4/Nla/lyxToO6V/F6WM6bLftsrU1jL7+OZas2fyklhV2uW2/PnQrzWZi5xIeemvLiYVbO6RfJfv0LOepj77jjH+8w6r1GzsTQo7F2o138Wf2bsWbX6/g6mc+BuCdb1fx1PxhTOtevlnvy+ljO7J8XQ03Pr91j4xIU6PvHlT3HMycy//k/e2a0/su+vDNtyI1G/bBn0krTU9KnqvVo5hcP1ydbO8qJBf/tur2rmJ21Qr8OlA18b9/WGj9mE3+/mGh9Yfx/xEOoe6zu1LdUOu6b+YWtxp76CX3JSVJBHh+4VK+XLGuTm1ro7GtTlQAq9bX8sj73wJQtEVPw7akeQ5HD23Dd6s3cNrf397sRAWwIRLd7DZaZW46Hy7ZWCP3+zUbWLq2horctB+fG1pdwOSupZz/6HtsiASxkpdI/ei758vML2L/n93sDZp2WK6xzqPW2l/h1/OTbdO5uh7Uo5hcP4x36LiN1zLxryAWktwxD9vTBr++09nUbaH1fzdCTEFxrLVnxuCCnmP3Zuzsn1jX2/kJoTGFHMvItoXEgBcWLttp+xFtCsgKuzzw5udYA1O7lVGVl87ydTU8/9lSPlyy+a2r71av3+xWV9i15IRdvovfrkvzHM7boxMPvLGIl7/Y+f5Fmovm8t0z1jJ8v6NNRec+zt+vPfPU2vXr967ZsG4/tl71RHSurhclisn1VPznRPyxCJvaA78H9ykahxZa97Xz0tLvdxyv16Rjfm7bDxgddDyAf7I4fWxHfxBMmseASv/C9ux/vLPNQe5b6l6aDfi9JH+bM5iSrI2dBzHg/te+5BdPfEA05s8LeuyDxZw4oh3jOhTxyhfLOWpINY41PPHhYgCOHdqGdM/h8qc/TO5/qEgT09y/e216DWXOlQ+4j93yi/YfvvTPV2JwNtHo5ex8LeeW5Kn4T52r60CJYnI9h/8/fDz+P7bH4s9n4c9mAn9poE1pofWGc4R13F9Xde3vTDz6Ai8jtyDoeH4Uciwze20+7nzBe9/yRh1OVACFGX4PxewBVbyxaAXHPfgGn36/hm6l2Zw/vjP7967gy+XruO2lhQD8/pUvaJOfwbXTegKwtibC+Y+9x7vfrqJLSRaH9q/i9P97m1XrI5w6qj0zepaT4Tm88Pkyzn/03c1mcIqkspbw3cvIyWfaKZc77z63gEdvuvDCaDQ6o3bDukPw100WnavrRYlickXxlyt6En9cwV/wi2lOwS+oeR1bdxtrofXk6+aF038fi0V77zbndKfnuL2DjmcrK9fX0uOKf2KNoSQrxG4dizlxRDtGtC3ggLtf4ZPvd3wc+qEe+NqaKMc9+AZL1/pDXV7+Yhkn//VN/nzoQGYPqPrxZLUhEuXcR97lkn9+SFFWiC+WraU2Gvtx9uZzn37Pgve+Zc7A1sweUMUV//qIT5eu4azdOnHFXj046J5XGvTzEGksLem712X4RKq6DQg9ccclAz586Z9vY8xF0draX9I4t1SbMp2r60GJYvI9jz/o9EJgAv4VyPvA0fhrPNZVKi+0HpQ0sGcbyxnVPQYy7vAznZyisqBj2qFoLMbXK9dz96tfsGZDhAsndGHe4GrO/MeOLxx/GED/+qLlP56ofvDBd6v5Yvk6WuelU5od5ptNeiRWbahl1fcbB98f0q+StgUZnPTXN3/8+9H3F/O7Vz4HIM39iCv26k6fVrn876uGGtct0vhayncvM7+IqSdf5nz86jM8cuOFZ25Ys/KQmg3rZuP3qrVkOlfXkRLFhvE6fkHNuijazvPHxR/1NWcHr/2w0PqZCWy3qRvnhdPuCGfklO5x5Flu+/5NYyxiffxQ7LtDUdZO236ydC0AK9dvu0DwyvjJLG0HdeEqctI4bnhbrnn2YxatWEe651CSFebNr1f+2ObNb/zbcW0LMpQoSrPVEr577fqNZO61D4Wfve837V/9x73POKHQgtr1644HPtrljacunavrQOVxJNV1dsNp/zDGPt5n/P5VR1zzUCgVk0SAyni5jJXra3bSEl6Mn9jaF2Zs9VrIsVTnpVMbjfHNqu3fYTp3j858uGQ19/z3y82e1zLX0tK0lO+eF05n7OyfMPuyP5hWHXvtbqx911p7OckvLC3NiBLFpqUnG5cWkh0rttb+1hj7VkWn3uNmX36/GX3wSXjhtJ2/M0D9KnJpU7D1CaYyN52zd/P/1z/54XebvTazVysO6V9Fmuf8+NynS9fwwsKltC/MZHr3zW+vzx/ahqywy78+XrLdJckmdy1laOt8znv0vR9nZ66tifD1yvX0KMv+sV2PshwAPt7JuC2Rpk7fvY2Kqjow89wbvemnXeVmF7c63vFCn2Ht0YCX9J01P63xJ5+0mPxJt56bloeBD4Hdgg6kCcvF2pMca0/PK692x80+zanuOdjZ+dsazpSupQyt9mdU/9DLML5TMW3jJ6WXv1jGA28uAmBE20LmDa7m4+/XsGjFOtbURCjNCtOjLBtrDM98soT7/rd5L8MJI9qRl+7xj3e/2ezk8/PH3+OuA/pz0cSu7NmllE+X+jMv+7TKZfHqDfzqn9ue4Jib5nH6mA7c/vJC3l+8arPX7nr1c04Z1Z45A1uzcOkafjKqPa8tWsFruu0sTZC+e7umfb9RtO09LPS/R/8Yeu6P119du379WZHaDecAd+Hf/pStnYk/XjAdqFu19xSXrEQxxOZrJ7b0GVWSfHnASdZxTs3ML3ZGHXBCepfhEzAm+Iu63q1ymbZFz0K30my6xeutRaKxH09WC977lkzPYWBVHj3LcsgMOaxYX8sLC5fx8Dtf87e3v/mxl2FnPlu6lv1+/xLHDW/HiDYFDG6dz5I1G7j/9a+44T+f8u2qbZfVOG1MB1ZviHD9fz7d6rU7X/mCvHSPwwe29kt0LFzKBY+/V+fPQqQx6bu366zj0m/SAfQcOy308t/vrnzxodtvgdgFNevXnYm/nJzqL7ZwuzIiYhj+YtZj8Ltif9hWDL+6+FP49YYaoor4DPzp6xUNsO0gfUbL7VGcBPTCXwh906N1qbX2BGPtyRm5Rc6I/Y8JdRs5GWOTnyD+fO+e0ecXLrXHP/hG0rfdnPQuz+HuA/tz7G2Pc9MTWvRhR87bZxjnzBhC+xNuZuGSlTt/Qwv2/IUHkZubFZt+x4saJbsDBRkhnj56OE+9+y1/f/2rpG9/7crlvPjgbby64N4I1n4e2bD+AuBuNi43BzAZ6I5/vG5prkc9ijtl8Gv8HB//Owp8xeZrJ1bhTxefjZ/QncjmJ3+RLf0j/vhBb8fzfhKNRGZlF7eKDJl+RLj76ClYR6MlREQaSnp2LqMPOZmBU2c7/33kvjav/uPemyOR2ktrazb8gmj0Zvwagn+PP6QFSOSsexR+kvgafgXzx/AriW8qE7/i+QXxtm9Rv7pE0jKFgL1C6RmnbVi7ZnBFl361A6cc4rbtPczVVFwRkcaTkVvA8JnHMGjaHPeNfz5U/OJDt/9q7YqlP49FIzdFo9EbaNlldVqURBLFo4FPgBHAqu20WQ08ADyOX6eovgUsJUlisVgmEDLGLA06lh3ob60zx1h7qHXctC4j9nT6TZxFUWV7dR+KiATIC6fTb+Is+o6fmfbhy0+lvfKPe0768p1XT3XDGc/UrFt9Jf4kzCY58SUWi4XwO65WGmNqd9Zeti2RE3En4Ga2nyRuaiX+P6K5CexHkuMGYBZNr+xBG2B/1wsfEand0KGqW/8NPXebEe4wcCyuFwo6NhER2YSxlo6DxtFx0Dj3uy8+4tUF9414+6m/DYPYymg0clc0Evkd0NTW+jwQf1WUsfjzJiQBiSSKa4HierQvJnkLaEtqa421+4TTM49Yv3pl9+LWHdd3GzUl3HX4RLIKSsJBByciIjtXVNme8XPPtmMPPtl+8OI/89965u/HLHzjheO8tPQva2vW3xytrf0jsOO1ECVlJJIoPoM/63gP/PGJOzIB2JvNJylIy+HiD1GY4oTC+0Y2rK8uKK1c123UlLTOQyeQX95ayaGISIry0jLoNmoy3UZNdlYvW8K7/15Q+c6z/zjrm4/fPt8Jp30Z3bD+vmg0+iDwH1RmJ2Ulkiiei58kPgI8CiwAPmDzWc+dgInxduvi75GWoR0w1jrOBDCTYrFoRqsOPWs6DBoXbt9/FAWt2jTtpVNERKTeMvMK6b/nQfTf86DwmuVL+PDlf1V+8MITJ3725osnG2tXG2OeqN2w/v+AJ/DnOUiKSCRRfB0YB9yK32M4YQdt3wKOAFSYrvnqCAxxvNB4Y53xtevXlmTkFmxo03uY167PcNOmzzDSMnPUcygi0kJk5BbSa7cZ9Npthlezfi0L33wx55P//Xv6R68+PWnld1+nOaHwV7HamgXRaPRp4Fk0g7pJS3RW6fP46xIPxy+43YnNV2Z5H3/g6HOou7k5aQX0BQa6XnhkNBYZGK2tzU7Pzl3fuvsgp7rnYLeq+0Dyy1trNoqIiOCF02nffzTt+482ux9xZtqyb77g09f+0+qLd1+d88Xbrx60aum3YddLW4qJPVu7Yf3TwH/jj+8DDl3idqX8SBR/vOIzSYpFmo4c/OS/C9DLC6cNiUYivSK1NbluKFxb0qZzbWWXfmnlHXtQ1r4H2YWl6jEUEZGdyiutpM/4/egzfj8DhNcsX8JXH7yR/8XbLNgzoQAAC+dJREFUr0z57I0Xdl/y5SehaKTWcb3wN8RiL9fWbngNv/Pp3fjPplzqrVlSnbqWy8UvUdMZ6Iy1XbxQuFc0EukYqdlQYIyJZeQVrSup7uSUtusWKqnuREnbLuSVVLgYo383IiKyyzJyC+kwYAwdBowxQHo0EuH7Lz/mm0/fK1382fuTl3z+0YTvvvw4sur7b0OxaNQ4Xmil43kfR2pqXo/UbHgLP3l8H3/pYK2T2QAa+oSfgd87BfB1A++rOfg7yfmcXKAUqLz++uuramtrLfALsFVuOFRtMNW1NesrYtGo66Wl1+SXVtUWV3cKFVS0dfLLqykob/3/7d1rjCRHYcDxf3fP7N7dvu727nx33IMzkgEBNkZC2DKOYnQ4RkFCIha2BAgESECiAEYi4ZQPQUgQKV9CeCiEgA4nyDGKEZgoMTbIYAmBwWBewoiHY9/5cfa9dm93z7ev6e58qB52bj27Mzu3O4+b/09q9Ux1z2x1bXVVTXd1Fdv2vDBKygOb1yEukiQ1JU4Sdhy4gh0HrqgGlYBSWllk6sTTTBx/fGTimWOvnDh+7MozTzw6P/HMsXj+/Mxg8dm5pFR+NocnK/MLRw8dOrT15ptv5sEHH7yO0Ih8GjgNXMzg2z8D7qJLBxnfCBvdUHwv8KnitXOwNfZXdcISQmN7BNhBGJdye+2SlMuXDZQG9gC7K2m6vbI4vy3P8wjg8OHD2f4DB+IDV15z29hlezeNju+KR3fuYeuufWx7wUG2jG4r032DcUuS9EdJqcz43oOM7z1YDYqBzQCLc+eZPnOCcxMnN52bPHVw5vSzB2cmT7727HNnKkeOHOHR/3vsYxB9EnIASqXyuaRUmiTnVCVdfCatVE4Ap4rldLHMFMsUMF28nidMOPLFth14F9johuIscGaD/8ZGKQPDNe9LhMZaVbUBt9L7mPCATzW8+vnq9w6Uy+WxJCltSZLSaJalY1Ecj+ZZtjnP86FKpTKSZelglmUXNOIGBzdVhoZHK8Mjo/no2LZ469bxZGzbeGlkdCujo1vZNr6DHZftYcfO3eza/QJe9fJ98a7xofxv7/rFlvVKGEmSukV50xa2772c7Xsvrw2OX31wfODW1xzg5789OXBq4hzHnzrG5MQppqYmh6fDsn96apLpqbPZ5MTpxanJiWx6aoKZmelkbvZ8UqksJhd8YRxX4iSZK5XK5+IkORcRzeTkc1EUna0szJ+vVCpn0zStNiingAVCA3OWMFQghKuZtbfIz1JtwYbpjxeK1wvF+6qO9c3stqt8X2Pp6enV7CMMy7LQaEdCgy1puFdIi2b2W/kLonjpUnQEcRRnUQRRHOdRFBHHSXVNHMdRHCd5HEd5FMUkSZxFUUwcx3kcx3mclPI/vo6TPI5jkiTJomht/7LBgWQoz9Py6dmmplzsa9s3RcMLeZQ9ObPgTEKrGBosJ+PxwuapxXzhxNnZZs7BvrV3fGhwpET58cn58/OLFUeAWMUVu0eH4jjhiemF5xrv3b9KSRzt25wP5fFAZWounWv8if41PFAaKGVzA+XBLbOLlWzNt4rzPCdNK1GWpVGaplGaVqIsTaM0TcmyNMrSLMrzjDRNozzPybIsrNOULM/yPMuiLMvI84zqXb48hzzP4urfyLKl12uQQpQ33o2UpUboqkdavLiTMKThBbqtofgBQr/GRoaBq1m5oTjPysPyzK4QnhWfa/V9t7oeuAr4l05HpAe8BzhKGBBWK4uBfwA+BzzV4bh0u5cRZrL6RKcj0gPeTLjzclenI9IDDgP3EMY11sp2AR8G/p7mLix12gBLF6yS4j2EtlrtZBWDhHJ4ueX71SrT+C7y3cBXmoppCwYIfed2snRg6g63EaZPUmP3EQoUra5E+JX66k5HpAe8CTjR6Uj0iM8DX+50JHrEMeCWTkeiB7yCUFYNN9pRK2vlkmfVdcDthCswc8DJYpkrwm4v9pEkSVIPauVhlgj4NOE2MYRbsMe5cK7n/cA7i+WzwIdo6j65JEmSukUrVxTfR2gk/pLQp2SU8HDJy4tlXxH2F4Q5nj9AGCZHkiRJPaSVhuJfAo8THpK4mwsf3656DvgGYS7oo8VnJEmS1ENaaSi+GPgfaGq4lZli35e08HckSZLUQa00FGcJTzc3ayfguHSSJEk9ppWG4vcJ/Q9vbGLfmwj9GL/fwt+RJElSB7Uy4PZVwI8Igzp+G7gX+AMXPvX8YuANhMbkHHAN4cEWtd8YYQyppzsdkR6wh5BfOzZVUo+IgMsJeaoXBp3vpC2EuyrHOh2RHrCDkLdOdToiPeAAYXpcZ7FZ3QDhAdujrDwJhzbItcAjhCFvVlt+TWgkSpIkqcdczBR+MeGp5hsIVxCrczRPAb8HHgB+gK14SZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSWu0GfgkcD/wJGFaupPAj4G/JkwrpmAUeDvwdeA4sAicBu4B/qyD8eoFbwX+FXiIkMdy4HBHY9SdrgTuJky5Ngv8AngfFzfhQS8z3zTP8qk51nnSGu0mFL6PE+bj/k/gW8BMEf5TwlzdChVUDpwHvgfcCfyQpWko/6ZzUet6jxLS6Awhr1nhP981hLw1D3wV+CxLafWZDsark8w3zbN8ao51nrRGZWBXnfAx4GHCifPOtsaoe90KvAcYWhZ+I7BA+AW/r92R6hGHgIPF69uwwl8uJsx3nwE31YSPAL8ipNd1HYhXp5lvmmf51BzrPGkdfZRw0nx0WfgI8I/A7wi3x/I6y4/aF82ucBfhuN+yLNy0er5mKvx+S7c/IRzTfXW2vbnYdns7I9SFzDets3xqjnXeRSp1OgJqq0HgzwknwHdrwkuEy/XXEvp0fIFwNeRtwNWEX64/AR5pZ2S7wHSxnq8JM61a04/pdkOxvrfOtm8TrjTeUGeblvRjvmmW5VNj1nlSA5sJnca/QPj1eaJY3rFsv3ezdCIlNeGDwM/pz1tkWwmdxmeK11WmVX2Nrgz1Y7p9hXBcb1xh+1OExuJA22LUfcw3rbF8qs86T1qjrTz/UvpXgZcu2+9/i21/Wuc7bim2/dvGRbPrxISnDHNCRVbLtKqvUYXfj+n234Tjun6F7Y8U23e2LUbdx3yzdpZPK7POU195Cji3huWGVb4rBvYDHyT8Ap3iwhPnMcKVjXpPhe0nnDQPtXwkG2890yoi/CLNgf+os73X06rWeqZbowr/Ukq3ZtlQbMx8szb9VD5djEu9zmsr+yh2r/t5/pNuqzmzyraMMLbUZwgnzRHg71i6HL+ZMPzCXJ3PnirWI2uIS7utV1pFwOcJY9zdAbyrzj69nla11jOPNXIppVuzpor12ArbxwgV0tQK29Wf+WYl/VY+XYxLvc5rKxuK3WujHuW/v1i/oibsLGEcqhHCSVVrb7Ge2KD4rIf1SKsY+CKh78odxXemdfbr9bSq1c7hIi6ldGvW74v1FXW2DQF7gCcIHedVXz/mm3r6sXxaL5dinddWcacjoLZ7UbE+WxP202L9+jr7HyrWD29YjDovIdzGeXexfgf1C2EwrVrVj+n2QLF+Q51tNxLK3wfqbNOSfsw3y1k+XRzrPKmO64GX1Al/EWEA4Bz4UE3464qwXwHDNeE7CCPdZ8CrNiSmnVcC/otw/Edo/OOpn9NqNY36mvVjutUOuH1jTfgwYRo/n6w03zRi+dQc6zxpjT5BOAl+QxjD7WuEgUMrRfg9hJHsa/1zse0o8E/Ap4FnirCPtSPSHfJxwjHOAv9OGAB5+bL8ilC/ptVyb2cpjR4iHP/DNWHL+1D1Y7pdy9IUfncS+kw9Rjjmfp3Cz3zTPMun5ljnSWt0FSHT/5LQz2KR0EH3O4S+LSv9Kn0XoeB+rlh+QBgq4FL2OeqPyF+7fKTO5/oxrZZrlHZfqvOZfky3q4BvEs7FWcJ5+X7Cwwn9yHzTPMun5ljnSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkte7/AQPHPw7DDzvaAAAAAElFTkSuQmCC)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dhcvnf4jdKUv"
      },
      "source": [
        "## 1.6 Standard Error of the Mean (SEM)\n",
        "\n",
        "The standard error or standard error of the mean is defined:\n",
        "$$SE = \\frac{\\sigma}{\\sqrt{n}}$$\n",
        "where $\\sigma$ is the standard deviation of the population and $n$ is the sample size. If you report SE and not SD, be sure you know the difference and be clear on what you report. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a7xVg4-DZOJP"
      },
      "source": [
        "# Get the SEM estimates of the setosa sample\n"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z_KPqXASdKUu"
      },
      "source": [
        "## 1.7 Quantiles, median and mode\n",
        "\n",
        "The **quantile $x_{\\alpha}$** is the value of the random variable x at which $\\alpha$% of the area is below x. An important special case is the **median, $x_{med} \\equiv x_{50}$**. At the median half the area lies above and half lies below. \n",
        "For the normal p.d.f. the median equals the mean. The most probable value of a distribution is called **mode**. \n",
        "\n",
        "Special quantiles are the quartiles and percentiles. The first quartile is the $x_{25}$, the second the $x_{50}$ etc. Percentiles are for example $x_{13}$ etc. \n",
        "\n",
        "<!--\n",
        "Any odd moment about the mean is a measure of the **skewness** of the p.d.f. The simplest of these is the dimensionless coefficient of skewness $\\gamma_1 = m_33/\\sigma^3$.\n",
        "\n",
        "The fourth central moment $m_4$ provides a convenient measure of the tails of a distribution. For the Gaussian distribution, one has $m_4 = 3\\sigma^4$. The **kurtosis** is defined as $\\gamma_2 = m_4/\\sigma^4 − 3$, i.e., it is zero for a Gaussian, positive for a leptokurtic distribution with longer tails, and negative for a platykurtic distribution with tails that die off more quickly than those of a Gaussian.\n",
        "\n",
        "The **quantile $x_{\\alpha}$** is the value of the random variable x at which the cumulative distribution is equal to $\\alpha$. That is, the quantile is the inverse of the cumulative distribution function, i.e., $x_{alpha} = F^{−1}(\\alpha)$. An important special case is the **median, $x_{med}$**, defined by $F(x_{med}) = 1/2$, i.e., half the probability lies above and half lies below $x_{med}$. (More rigorously, $x_{med}$ is a median if $P(x \\geq x_{med}) \\geq 1/2$ and $P(x \\leq x_{med}) \\geq 1/2$. If only one value exists, it is called ‘the median.’)\n",
        "\n",
        "Under a monotonic change of variable $x \\rightarrow y(x)$, the quantiles of a distribution (and hence also the median) obey $y_{\\alpha} = y(x_{\\alpha})$. In general the expectation value and **mode** (most probable value) of a distribution do not, however, transform in this way.\n",
        "\n",
        "Let us look at median and quantile$_68$ of the normal pdf: \n",
        "-->"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mEXXg2DbahK6"
      },
      "source": [
        "# Let's calculate some quantiles\n"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eW6ny3DwdKUz"
      },
      "source": [
        "## 1.8 Box and scatter plots\n",
        "\n",
        "Often box plots are used to visually inspect differences between categories. Let's display the setal length of the three Iris categories with bx plots (https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.boxplot.html):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5unmWh2fdWTC"
      },
      "source": [
        "# Plot the boxplots for the sepal lenght for the 3 Iris species\n"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bv2uyf7ydpjq"
      },
      "source": [
        "A visual inspection of correlations between columns is often useful. Scatter plots provide this."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sKzgV4PYdZh9"
      },
      "source": [
        "# Plot all scatter plots for the setosa species\n",
        "# https://pandas.pydata.org/pandas-docs/stable/generated/pandas.plotting.scatter_matrix.html\n",
        "from pandas.plotting import scatter_matrix\n"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ETxa8KgodKUt"
      },
      "source": [
        "## 1.9 Uncertainties\n",
        "\n",
        "All data have uncertainties. These should always be communicated when showing scientific numbers or plots. We distinguish between two types.\n",
        "\n",
        "- Statistical uncertainties \n",
        "    - Fluctuations, can be made smaller by taking more data, i.e. get more statistics (see definition of standard error)\n",
        "- Systematic uncertainties \n",
        "    - Shift of data in one direction due to some \"mistake\" in the measurement, e.g. wrongly calibrated instrument showing all measured values systematically higher as they really are. Or for instance uncertainty due to the choice of methods and tools\n",
        "    \n",
        "The statistics tools can mostly handle the statistical uncertainties. There is no mathematical recipe for dealing with systematical uncertainties. You have to think through your experiment and try to estimate the influence of everything that can go wrong.\n",
        "\n",
        "When uncertainties are stated on numbers or in graphs as error bars or error bands they generally show one standard deviation. If the data are well described by a normal p.d.f, the interpretation of one standard deviation is clear: if the measurement is repeated many times, 32% (or about 1/3 of the measurements) **should** be outside the error bars. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AoeTElgmdKUz",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "outputId": "e552aa47-c0c3-4ed9-8a34-e8c073f7dae4"
      },
      "source": [
        "# If we have time, we study this code which illustrates the above statement\n",
        "import scipy\n",
        "n = scipy.stats.norm.rvs(3.418,0.318,100)\n",
        "# Draw a histogram which is not normalised\n",
        "entries1, edges, patches = plt.hist(n, bins=10, histtype='step')\n",
        "# Draw a histogram which IS normed\n",
        "entries2, edges, patches = plt.hist(n, bins=10, histtype='step',density=True)\n",
        "# Close plt so that the previous histograms are not shown\n",
        "plt.close() \n",
        "# Calculate the poisson standard deviation and scale down to second histogram\n",
        "errors = np.sqrt(entries1) * entries2/entries1\n",
        "# calculate bin centers\n",
        "bin_centers = 0.5 * (edges[:-1] + edges[1:])\n",
        "# draw errobars, use the sqrt error.\n",
        "plt.errorbar(bin_centers, entries2, yerr=errors, fmt='r.')\n",
        "# Draw a normal distribution\n",
        "x = np.linspace(2.4,4.2,100)\n",
        "plt.plot(x,scipy.stats.norm.pdf(x,3.418,0.318))\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZd7G8e8vnRpKQpEWCEiRIhCabcGKvaGiAmJZ7Lvrur6767ruvvpuL7q6KiKwimBlxYquiqKitASp0iIdAoQWEtKT5/1jBjdGMBOY5ExO7s91zcWZmZM5t8fJnZnTHnPOISIi/hPldQAREakZKngREZ9SwYuI+JQKXkTEp1TwIiI+FePVgpOSklxKSopXixcRqZMyMjL2OOeSQ5nXs4JPSUkhPT3dq8WLiNRJZrY51Hm1iUZExKdU8CIiPqWCFxHxKRW8iIhPqeBFRHxKBS8i4lMqeBERn1LBi4j4lApeJJyGDw/cRCKACl5ExKdU8CIiPqWCFxHxKRW8iIhPqeBFRHxKBS8i4lMqeBERn1LBi4j4lApeRMSnVPAiIj6lghcR8SkVvIiIT6ngRUR8SgUvIuJTKngREZ9SwYuI+JQKXkTEp6oseDPrYGYfm9lXZrbKzH58hHnMzB4zs0wzW25mA2omroiIhComhHlKgXudc0vMrAmQYWYfOOe+qjDP+UC34G0I8FTwXxER8UiVn+Cdc1nOuSXB6VxgNdCu0myXAtNcwAKgmZm1DXtaEREJWbW2wZtZCtAfWFjpqXbA1gr3t/HdPwKY2QQzSzez9Ozs7OolFRGRagm54M2sMfBv4CfOuYPHsjDn3CTnXJpzLi05OflYXkJEREIUUsGbWSyBcp/hnHvtCLNsBzpUuN8++JiIiHgklKNoDJgCrHbO/f0os70JjAseTTMUyHHOZYUxp4iIVFMoR9GcCowFVpjZ0uBj9wMdAZxzE4HZwAVAJpAP3Bj+qCIiUh1VFrxzbh5gVczjgDvDFUpERI6fzmQVEfEpFbyIiE+p4EVEfEoFLxJOOTmwZQvMn+91EhEVvEjYzJ8Py5fDxo1w1lkqefGcCl4kXObOhfLywHRxceC+iIdU8CLhMnw4RAV/peLiAvdFPKSCFwmXYcOgb1/o3BnmzAncF/FQKGeyikioEhMDN5W7RAB9ghcR8SkVvIiIT6ngRUR8SgUvIuJTKngREZ9SwYuI+JQKXkTEp1TwIiI+pYIXEfEpFbyIiE+p4EVEfEoFLyLiUyp4kbpk+HBdhlhCpoIXEfEpFbyIiE+p4EVEfEoFLyLiUyp4ERGfUsGLiPiUCl5ExKdU8CIiPqWCFxHxKRW8iIhPqeBFRHxKBS8i4lMqeBERn6qy4M1sqpntNrOVR3l+uJnlmNnS4O3B8McUEZHqiglhnmeBfwLTvmeez5xzF4UlkYiIhEWVBe+c+9TMUmo+iogPzJ3rdQKRb4RrG/wwM1tmZu+a2UlHm8nMJphZupmlZ2dnh2nRIiJyJOEo+CVAJ+dcP+Bx4PWjzeicm+ScS3POpSUnJ4dh0SIicjTHXfDOuYPOubzg9Gwg1sySjjuZiIgcl+MueDNrY2YWnB4cfM29x/u6IiJyfKrcyWpmLwLDgSQz2wb8BogFcM5NBEYBt5tZKVAAjHbOuRpLLCIiIQnlKJprq3j+nwQOoxQRkQiiM1lFRHxKBS8i4lMqeJG6JCcHtmyB+fO9TiJ1gApepK6YPx+WL4eNG+Gss1TyUiUVvEhdMXculJcHpouLdVkEqZIKXqSuGD4cooK/snFxgfsi30MFL1JXDBsGfftC584wZ07gvsj3COVywSISKRITAzeVu4RAn+BFRHxKBS8i4lMqeBERn1LBi4j4lApeRMSnVPAiIj6lghcR8SkdBy8SJqVl5Wzam8/anblsyM4jp6CE3MJSDhWXkhAbTdOEWJo2iKFzUiN6tGlKl+RGxEbrM5bUHBW8yHHYsjefj9bs4uO12SzcuJfCkvJvnmsUF03jhBgaxcVQUFJGbmEpeUWl3zwfFx1F/47NOLNHK87s0YqurRoTHP1SJCxU8CLVVFBcxtvLd/DCoi18ueUAAF2SGjF6UEd6t0ukR5smdG3VmITY6O/8bHFpORv25LF2Zy6rdhzk03XZ/OHdNfzh3TX0aNOE64d05LL+7WiSEFvb/1niQ+bV8KlpaWkuPT3dk2VLPXb4Al3HcCXGnPwSnvlsA9Pmb+JgYSldkhtx7aCOnNOrNSlJjY450vYDBcxZvYuXF29l1Y6DNIyL5uq0DtwxPJVWTRPCll/8wcwynHNpocyrT/AiVcgrKmXKZxuZPG8DuYWlnN+7DTecksKQzi3CskmlXbMGjBuWwtihnVi+LYdp8zfz/ILNvLR4C+OGpXDbD1Jp0SguDP8lUt+o4EWOwjnHuyt38r9vrWLXwSLO7dWae845kZ5tm9bI8syMfh2a8bcOzbj7zK48Nmc9z3y2gZcXb+WX5/fg6rQOOuxNqkUFL3IEW/fl8+s3VjJ3bTa92jblqTEDGdCxea0tPyWpEX+/5mRuG57KA6+v5BevreDVjG38vkES3Qv21FoOqdtU8CIVOOeY9eV2HnxjFc45HryoF+OGdSLGo8MZT2zdhJcnDOXfS7bz+9mrubjPWH6+5VNuLHdERemIG/l++sYnEnSwsIQfv7SUn76yjF5tm/Kfe87gptM6e1buh5kZowa254N7zuCMA5t4OOVMxj+7mN25hZ7mksinghcB1u3K5eLH5/HOiizuPedEXpwwlPbNG3od61taNo7nmXWzeHjjByzcsJcL/jGP9E37vI4lEUwFL/XeeyuzuOyJz8kvLuPlCUO5+6xuREfo5g8Dxu5aypt3nUbj+GiufWYBLyzc4nUsiVAqeKm3nHM8+uE6bpu+hBNbN+Htu08jLaWF17FC0r1NE9648zSGpSZx/6wVPPD6CkrLyqv+QalXVPBSL5WUlfM/M5fz6IfruXJAe16aMJTWlU8qinCJDWP51/hB3HpGF6Yv2MKtz2eQX1xa9Q9KvaGCl3rnUFQstzyXzqsZ2/jxWd3461V9j3hZgbogOsr45QU9efiy3ny8djfXPbOQvXlFXseSCKGCl3rlQHQC1/a6hnmZe/jjFX2455wTfXGBr7FDO/HUmIGszjrIVRPns+NAgdeRJAKo4KXe2JNXxOhe17CmYTJPjxnI6MEdvY4UVued1IYZtwwhO6+Iq5+ez9Z9+V5HEo+p4KVe2HWwkNGTFrApoTlT1r7G2b1aex2pRqSltOCFW4aSV1TKVRPnsyE7z+tI4iEVvPje7mC5Zx0o4Lk1Mzk9Z7PXkY7d3LlVXkmyT/tEXvzhUErKyrn66QVk7lbJ11cqePG1PXlFXDd5IbsPFjLt5sEMyd3mdaRa0bNtU16+dSjguH7yAjbvPeR1JPGACl58a/+hYsZMXsi2/flMHT+IgZ3qxjHu4dK1VROm3zKEotJyrnsmsB6kfqmy4M1sqpntNrOVR3nezOwxM8s0s+VmNiD8MUWqJ7ewhHFTF7FhzyEmjxvEkC4tvY7kiR5tmjL95iEcLCzh+uA3Gak/QvkE/yww8nuePx/oFrxNAJ46/lgix66wpIwJ0zJYnXWQp64fwGndkryO5Kne7RKZdtNgsnOLGDd1ETkFJV5HklpSZcE75z4Fvu+KRpcC01zAAqCZmbUNV0CR6igrd/zkpaXM37CXv17Vj7N6+vNomerq37E5T48dyNfZefzwuXQKS8q8jiS1IBzb4NsBWyvc3xZ8TKRWOed44PUVvLdqJw9e1IvL+uttWNHp3ZJ55JqTWbx5H3e98KWuXVMP1OpOVjObYGbpZpaenZ1dm4uWeuCxOZm8uGgrd45I5abTOh95ppwc2LIF5s+v3XAR4qK+J/DQJSfx4epd/Do4qMl3DB/+38G9pU4LR8FvBzpUuN8++Nh3OOcmOefSnHNpycnJYVi0SMDMjG088uE6rhzQnp+d2/3IM82fD8uXw8aNcNZZ9bbkxw5L4Y7hqby4aAtPffK113GkBoWj4N8ExgWPphkK5DjnssLwuiIh+Wx9Nr/493JO65rEH67oc/Rry8ydC+XBzRLFxVWeMORnPzu3O5eefAJ/fm8tbyw94ucx8YEqx2Q1sxeB4UCSmW0DfgPEAjjnJgKzgQuATCAfuLGmwopUtnZnLrdPX0LXVo15cswA4mK+5zPL8OEQFRUo+bi4er0ZIirK+POovuzMKeRnry6jddMEhtbTQ0n9zI64Da4WpKWlufT0dE+WLf6QnVvEZU98TklZOW/cdSptExtU/UP9+we2w8+YAcOG1XzICJeTX8IVT33O3kPFzLrjVDonNfrvH756/A0nkplZhnMuLZR5dSar1EmFJWX8cFo6+w4VM+WGQaGVO0BiInTsqHIPSmwYy9TxgzDg5mcXcyC/2OtIEkYqeKlznHP87NVlLNt2gEeuOZk+7RO9jlSndWrZiEnj0ti2v4Dbpy+hxFQLfqH/k1LnPDYnk7eXZ/HzkT0Y2buN13F8YVBKC/40qg/zN+zlwZSz8WbDrYRblTtZRSLJuyuyvjkc8tYzungdx1cu79+e9bvyeHIu9MjP5gavA8lx0yd4qTNWbs/hp68sY0DHZvz+it6+GGov0vzs3O6cvW89D6Wcybz1e7yOI8dJBS91QnZuEROmpdOsYSwTxw4kPqZuDpId6aKijEcz36FrwV7umJHBxj26jnxdpoKXiFdcWs7t0zPYl1/MM+PSaNUkwetIvta4vITJa2cRHWX8cFo6uYW6+mRdpYKXiPebN1eRvnk/fxnVj97tdMRMbehQlMMT1w9g455D3PPyUsrLtdu1LlLBS0SbvmAzLy7awh3DU7m43wlex6lXTklN4tcX9uTD1bt59MN1XseRY6CjaCRiLdywl9++uYoR3ZO592gXEJMadcMpKXyVdZDHPsqkZ9umnN9HQz3UJfoELxFpx4EC7pixhI4tGvKPa/sTHaUjZrxgZjx8WW/6d2zGva8uY+3OXK8jSTWo4CXiFJaUcdv0DIpKy5k0Lo2mCbFeR6rX4mOimThmII3iY5jwfDo5+drpWleo4CWiOOf41ayVLN+WwyPXnEzXVo29jiRA66YJTBwzgB0HCrj7pS8p007XOkEFLxHl2S828e8l27jn7BM5p5fGU/XEUUa9GtipBQ9d2ptP12Xzl/+s9SicVIcKXiLGgg17+b93VnNOr9bcfWZXr+PUT1WMenXt4I5cN6QjEz/5mreX7/AopIRKBS8RYfuBAu6csYSUlg35+9X9iNJOVW+EMOrVby7uxYCOzbjv1eWs2XmwVuNJ9ajgJXQ1NBhzYUkZtz2fQXFOLpM++AdNtFPVO4dHvYKjjnp1eKdrk4QYJkzL0DXkI5gKXjzlnOP+WStYsT2HRzLfIbVwv9eR6rdhw6BvX+jcGebMOerAKK2aJvDUmIFk5RRw94va6RqpVPDiqWe/2MRrS7Zzz9kncvaBr72OIxDyqFcDOzXnoUt789n6PdrpGqFU8OKZL77ew/+9s5pztVO1zrp2cEeuD+50fWuZdrpGGl2qQDyxbX8+d73wJZ2TGvG32typqoGkw+43F5/E2p253DdzGanJjel1QlOvI0mQPsFLrSsoLmPCtAxKSsuZNHagdqrWcXExUTw5ZgDNGsR9MxC6RAYVvNQq5xz3zVzG6p0Heeza/nRJ1pmqftCqSQITxw4kO6+IO2csoaSs3OtIggpeatnETzbw9vIs7juvOyN6tPI6joTRyR2a8YfLAwN3/+6d1V7HEbQNXmrRx2t38+f/rOGivm25/QepXseRGnDlwPas2nGQqZ9vpNcJTbk6rcN3Zzp8bL32h9Q4fYKXWpG5O48fvfAlPds05S+j+mnAbB+7/4IenNq1JQ/MWknG5n1ex6nXVPBS43LyS/jhtHTiY6N45oY0GsRpwGw/i4mO4onrBtC2WQK3Pr+EHQcKvI5Ub6ngpUaVlpVz14tL2LY/n4ljBtKuWQOvI0ktaNYwjsnj0igsKeOH09IpKC7zOlK9pIKXGvW72av5bP0efndZH9JSWngdR2pRt9ZNePza/nyVdZCfvbpMA3d7QAUvNWbGws386/NN3HRqZ64edISdbeJ7I3q04pfn9+CdFVk8Ome913HqHR1FIzXi88w9PPhGYMDsX13Y0+s44qEfnt6F9bvyeGzOelKTG3Gp14HqEX2Cl7D7OjuP26dn0DW5MY9pwOx6z8z43eV9GNy5BffNXM6Sxm29jlRvqOAldEcZyq2ifYeKufnZxcRGRzH5hrTqXYYghNeXuikuJoqJYwbSpmkCE068nK3xul5NbVDBS2iqGMoNAgN3TJiWzo6cQiaNG0iHFg3D+vpSt7VoFMfU8YMojormxu5XklNQ4nUk31PBS2iqGMqtvNxx38zlpG/ezyNXn8zATtU8YiaEoeKk7uvaqjFPr3udzQnNuX16BsWlumZNTVLBS2iqGMrtbx+s5a1lO/j5yB5c2PcYtrGGMFSc+MOwg1v544b/8MXXe7l/1gqc0+GTNSWkgjezkWa21swyzewXR3h+vJllm9nS4O2W8EcVT33PUG7TF2zmiY+/5trBHbjtB13C/vpSy+bOrfFvUFfuWcVPzu7GzIxtPPKhDp+sKVUeJmlm0cATwDnANmCxmb3pnPuq0qwvO+fuqoGMEikSEwO3CuX7/qqdPPjGSs7s0YqHL+19fNeYOcLri3/9+Kxu7DhQwGNz1tOmaQLXDenodSTfCeU4+MFApnNuA4CZvQRcClQueKlnMjbv5+4Xv6RPu0T+eV1/YqK1xU9CkJMDOTnYggX87vIhZOcW8cDrK2jVJJ6ze7X2Op2vhPIb2Q7YWuH+tuBjlV1pZsvNbKaZHfG0RTObYGbpZpaenZ19DHElUqzflcvNzy2mbWICU8YPomGczpmTEFQ6Wip20UKeuH4AfdolcteLS0jfpKtPhlO4PnK9BaQ45/oCHwDPHWkm59wk51yacy4tOTk5TIuW2rZtfz5jpywiNjqK524aTFLjeK8jSV1xhKOlGsbFMGX8IE5IbMBNzy5mddZBTyP6SSgFvx2o+Im8ffCxbzjn9jrnioJ3JwMDwxNPIs2emIaMnbKI/OJSpt00mE4tG3kdSeqSoxwtldQ4nmk3D6ZhXAzjpi5i895DnkX0k1AKfjHQzcw6m1kcMBp4s+IMZlbxuLhLAI3X5UM50fHc0HMUWTkFTB0/iJ5tdTaiVNP3HC3VvnlDnr95MCVl5YydsohdBws9DOoPVRa8c64UuAv4D4HifsU5t8rMHjKzS4Kz/cjMVpnZMuBHwPiaCizeyCsq5cYeV7KuQRJPjRmoS//KsUtMhI4dj3i0VLfWTXj2xsHszSviumcWsCev6AgvIKEKaRu8c262c+5E51yqc+53wccedM69GZz+pXPuJOdcP+fcCOfcmpoMLbWroLiMW55bzLLGbXl8/VuM6K7BsqXmnNyhGVPHD2L7gQLGTF7IgfxiryPVWTquTb5XYUkZt07PYOHGffw9czYj9+ukFKl5Q7q0ZPK4QWzYc4ixUxbpujXHSAUvR1VYUsatz2fw6bps/nRFXy7dq10rUntO65bExDEDWLPzIOOmLFTJHwMVvBzR4bE0P12fzZ+v7KsRmcQTZ/ZozVPXD+SrrIOMnbKQnHyVfHWo4OU7CooD5T4vcw9/UrmLx87u1Zqnxw5kTVYu109ZwP5D2iYfKhW8fMvBwhLGTV3IvMw9/GVUP65OU7mL987s0Zqnxw1k3a48Rk9awG4dQhkSFbx84/ChaV9uOcDj1/Zn1MD2XkcS+caI7q341/hBbN2fz1VPz2frvnyvI0U8FbwAsONAAddMWsD6XXk8My6Ni/qe8N2ZavoysrVwmVqp207tmsSMW4ZwIL+EqybOJ3N3rteRIpoKXli7M5crnvyCXTmFPHfTYEb00HHuErn6d2zOSxOGUlruuPKp+SzWBcqOSgVfzy3YsJdRE7/A4XjltmEM7dLS60giVerZtimz7jiFlo3iuH7yQt5bmeV1pIikgveT4cOrNdTdG0u3M27KIlo3TeC1O07VtWWkTunQoiEzbz+Fk05oyu0zljBl3kYN/1eJCr4eKi93/O39tfz4paX079iMmbcNo12zBl7HEqm2Fo3ieOGWoZzbqzUPv/0V989aSUmZBvI+TKM01DMFxWXc++pSZq/YyTVpHXj4st7ExejvvNSiMO9IbxAXzVPXD+Sv76/lyblfs2nPIZ4aM4BmDePCupy6SL/Z9cjmvYe4/MnPeXflTh64sCd/vLKPyl18ISrK+J+RPfj71f3I2Lyfi/85j1U7cryO5Tn9dtcTH6/ZzcWPzyMrp5B/jR/ELad3Ob4BskUi0BUD2vPSrUMpKXVc8eQXvLZkm9eRPKWC97nSsnL+9v5abnpuMe2bN+Ttu09juC73Kz42oGNz3rr7NPp3bMZPX1nGr2atoLCkzOtYnlDB+9j2AwWMnrSAxz/KZNSA9vz79lPo0KKh17FEalxyk3im3zyEW8/owoyFW7jsic/r5UlRKnifendFFhf84zPW7MzlH6NP5i9X9aNBXLTXsURqTUx0FL+8oCf/unEQ2blFXPT4PF5YuKVeHUqpgveTnBwKt27nkYee5fYZS+jUMrBJ5tKT23mdTMQzI7q34t0fn87ATs25f9YKbnx2cb0Z79W8+muWlpbm0tPTPVm2L82fT9lpp+OcoyQ6hjf/9jxX3DmK2Gj9DReBwPkf0+Zv4o/vrSE+JprfXtKLy05uV+cONjCzDOdcWijz6rffB3bnFvLGozNwzhHjyklwZVxz6GuVu0gFUVHG+FM7M/tHp5Oa3Ih7Xl7GuKmL2LLXv1elVAPUYWXljukLNnP23z5hRkIXXFQUDrC4uGpdskCkPumS3JhXbzuF/73kJJZs3s+5j37Ck3MzKSr135E22kRTRy3csJffvvUVq7MOMrRLC353eR9Szz0dcnJgxgwYNszriCIRLyungN+8sYr3v9pFSsuGPHBhL87q2SqiN9tUZxONLlVQx2zIzuOv769l9oqdnJCYwD+v68+FfdoG3pCJiYGbyl0kJG0TGzBpXBqfrMvmobdWccu0dE7vlsTPR/agd7tEr+MdNxV8HbHrYCGPfrieV9K3Eh8TxU/O7satZ6Tq0EeRMPjBicm895MzeO6LTTz+USYXPT6Pi/udwL3nnEhKUiOv4x0zFXyE23GggImffM1Li7finGPs0E7cOaIryU3ivY4m4iux0VHccnoXrkrrwDOfbmDKvI3MXpHFZSe3444RqaQmN/Y6YrWp4CPU+l25TJm3kX8Hr6Vx5YD23Dmiq85EFalhiQ1i+dl53Rl3Siee/mQDMxZu5rUvt3Fhn7ZMOKMLfds38zpiyFTwEaS83DEvcw9TP9/I3LXZxMdEMXpQR24bnqrrtYuE6vARZMd5WeJWTRL49UW9uH14KpM/28j0BZt5e3kWg1NacNN/pnL2/kxi5n583HFrkgo+AuzNK+LVjG28uGgLm/fmk9Q4nnvPOZHrh3aiRSNd01rES0mN4/nF+T24c0QqLy/eyrNfbOK27pfRpiiXqz9Yx+hBHTghQj+AqeA9UlhSxkdrdvPaku3MXbub0nLH4M4t+Ok5JzKydxviY7TzVCSSNEmI5ZbTuzD+lBTmjJrAC6378fhH63n8o/WcmprE5f3bMbJ3GxrFR06tRk6SeqCwpIx56/cwe2UWH3y1i9zCUlo1ieem0zpzdVp7urZq4nVEEalCTHQU5+3P5Lz9mWx9bTavpm9l1tLt3PvqMh54fSUjeiRzfu+2nNmjledlr4KvYbtzC5m7JpuP1uzms/XZHCouI7FBLOf2asNl/U/glNQkoqPCdFJFmIdCE6mTcnICt/nza+6ckOAyOqxdxk/PHcY955xIxub9vL50O++t3MXsFTuJj4liWGpLzuzRihHdW3lygITOZA2znIIS0jft4/PMvXzx9R7W7Axcg7ptYgIjerRi5EltGJbaUteJEakJ8+fDaadBeTk0aABz5oS/5KtYRlm5I33TPt5btZOP1uxmc/BaN12SGnFK15acmprEsNSWxzxmrM5krSVl5Y6vs/NYtvUAS7ceIGPzftbuysU5iI+JIi2lOfed150R3VvRs22TiD79WcQX5s4NFC9AcXHgfrgLvoplREcZQ7q0ZEiXljx4US827DnEx2t283nmHmYt2c70BVu48dQUfnPxSeHNdQQq+BDtP1RMZnYea3fmsmbnQVZn5bI66yD5xYELFDWOj2FAp+Zc2KctaSkt6N+xGQmx2lEqUquGD4eoqEAB19RF96qxDDMjNbkxqcmNueX0LpSUlbN824Fj/vReXSr4oPJyx+7cIrYfyGfb/gK27M1n0958tuw7xIbsQ+w9VPzNvE3iY+jZtilXp3WgT7tE+nVIpEtSY6LCtS1dRI7NsGHQt2/NXnTvOJYRGx3FwE4twp/pKEIqeDMbCfwDiAYmO+f+WOn5eGAaMBDYC1zjnNsU3qjVV1xaTk5BCTkFxew7VMK+Q0XsO1TCnrwi9uQVkZ1bxK6DhezMKWR3bhGl5d/eH9G6aTydWjbi7J6t6dqq8Te39s0baHOLSKSqjYvu1ZEL+1VZ8GYWDTwBnANsAxab2ZvOua8qzHYzsN8519XMRgN/Aq6picDrduXyzvIsCkrKyC8uJb+4jENFpRwqKiOvqJTcwhJyC0vJLSyl4HtGUk9sEEtS4zhaN01gaGpL2jRNoG2zBrRv3oD2zRrQrnkDGsbpC46I1F2hNNhgINM5twHAzF4CLgUqFvylwG+D0zOBf5qZuRo4RCdzdx7/mLOehNgoGsbF0CA2msbxMTSKj6ZJQgwnNEugSXwsTRJiSGwQS7OGsTRtEEuLRnHfuulEIhHxu1AKvh2wtcL9bcCQo83jnCs1sxygJbCn4kxmNgGYANCxY8djCnzeSW3Y8PsLtL1bRI6sNs4HqSPnnNTqwdjOuUnOuTTnXFpycvIxvUZ0lKncRURCEErBbwc6VLjfPvjYEecxsxggkcDOVhER8UgoBb8Y6GZmnc0sDhgNvFlpnjeBG4LTo4CPamL7u4iIhK7KbfDBbep3Af8hcJjkVOfcKjN7CEh3ztO3LlIAAAW0SURBVL0JTAGeN7NMYB+BPwIiIuKhkI4DdM7NBmZXeuzBCtOFwFXhjSYiIsdDV7wSEfEpFbyIiE+p4EVEfEoFLyLiU54N+GFm2cBmTxYekESlM20jWF3JqpzhVVdyQt3J6oecnZxzIZ0p6lnBe83M0kMdFcVrdSWrcoZXXckJdSdrfcupTTQiIj6lghcR8an6XPCTvA5QDXUlq3KGV13JCXUna73KWW+3wYuI+F19/gQvIuJrKngREZ/yXcGbWQcz+9jMvjKzVWb24yPMM9zMcsxsafD2YIXnRprZWjPLNLNfeJzzvgoZV5pZmZm1CD63ycxWBJ9Lr6mcwWUlmNkiM1sWzPq/R5gn3sxeDq63hWaWUuG5XwYfX2tm53mc86fBdb7czOaYWacKz5VVWN+VL4ld2znHm1l2hTy3VHjuBjNbH7zdUPlnaznnIxUyrjOzAxWeq5X1WWF50Wb2pZm9fYTnPH9/ViNr+N6jzjlf3YC2wIDgdBNgHdCr0jzDgbeP8LPRwNdAFyAOWFb5Z2szZ6X5LyZwnf3D9zcBSbW0Tg1oHJyOBRYCQyvNcwcwMTg9Gng5ON0ruB7jgc7B9RvtYc4RQMPg9O2Hcwbv50XQ+hwP/PMIP9sC2BD8t3lwurlXOSvNfzeBy4nX6vqssLyfAi8c5Xfb8/dnNbKG7T3qu0/wzrks59yS4HQusJrAmLGh+GaAcedcMXB4gPFIyHkt8GJNZKmKC8gL3o0N3irvnb8UeC44PRM4y8ws+PhLzrki59xGIJPAevYkp3PuY+dcfvDuAgIjlNWqENfn0ZwHfOCc2+ec2w98AIysgZjHktOz96iZtQcuBCYfZRbP35+hZg3ne9R3BV9R8GtYfwKfPCobFvzq+a6ZnRR87EgDjIf6x+GYVZETM2tI4Jf43xUedsD7ZpZhgcHMazpjtJktBXYTKJjKWb818DpweOD1Wl2nIeSs6Gbg3Qr3E8ws3cwWmNllNZWxGjmvDH5Nn2lmh4fNjMj1GdyM0Bn4qMLDtbY+gUeB/wHKj/J8RLw/g6rKWtFxvUd9W/Bm1phAIf7EOXew0tNLCFzPoR/wOPB6bec7rIqch10MfO6c21fhsdOccwOA84E7zeyMmszpnCtzzp1M4NPEYDPrXZPLO1ah5jSzMUAa8JcKD3dygdPDrwMeNbNUD3O+BaQ45/oS+JT+XOXXqA3V+P8+GpjpnCur8FitrE8zuwjY7ZzLqInXD6fqZA3He9SXBW9msQRKc4Zz7rXKzzvnDh7+6ukCo1XFmlkSoQ0wXms5KxhNpa++zrntwX93A7Oo4a+VFZZ7APiY724WONrA67W6TkPIiZmdDfwKuMQ5V1ThZw6v0w3AXALfqjzJ6ZzbWyHbZGBgcDri1mfQ971Ha3p9ngpcYmabCGxWPdPMpleaJ1Len6FkDd97NJw7DiLhRmDH0DTg0e+Zpw3/PclrMLAl+HMxBHZadea/O1lP8ipncL5EAuPcNqrwWCOgSYXpL4CRNbhOk4FmwekGwGfARZXmuZNv78R6JTh9Et/eibWBmtvJGkrO/gR2pHWr9HhzID44nQSsp+Z2sIeSs22F6cuBBcHpFsDGYN7mwekWXuUMPteDwE5/82J9VsoynCPvuPT8/VmNrGF7j4Y0JmsdcyowFlgR3HYIcD/QEcA5NxEYBdxuZqVAATDaBdbaEQcY9zAnBH6533fOHarws62BWYF9RMQALzjn3quhnBA44uc5M4sm8K3vFefc2xbCwOsuMED7K8BXQClwp/v21/jazvkXoDHwanD9bXHOXQL0BJ42s/Lgz/7ROfeVhzl/ZGaXEFhn+wgcVYNzbp+ZPQwsDr7WQ+7bm+5qOycE/l+/FPwdOqw21+cRReD7M9SsYXuP6lIFIiI+5ctt8CIiooIXEfEtFbyIiE+p4EVEfEoFLyLiUyp4ERGfUsGLiPjU/wNGrmGRjrGquAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G_NqjEomdKUt"
      },
      "source": [
        "## 1.10 Recapitulation\n",
        "\n",
        "- Collected data almost always contain random variables. Why?\n",
        "- Random variables can be modelled with pdfs. The scipy.stats has many built in pdfs.\n",
        "- There are many pdfs (infinite). We know the most important PDF\n",
        "- We know the p-values for different 1, 2 and 3 standard deviations on the normal PDF.\n",
        "- We know how to print descriptive stats of a dataset and plot histograms and scatter plots (correlations)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XfC99sWpe_ZY"
      },
      "source": [
        "## 1.11 Exercise\n",
        "\n",
        "Now use the Iris virginica class.\n",
        "\n",
        "- Print the descriptive statistics\n",
        "- Plot the histograms and scatterplots\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q_XB84FUdKUj"
      },
      "source": [
        "## Annex 1 The Poisson pdf\n",
        "\n",
        "The Poisson distribution is popular for modelling the number of times an event occurs in an interval of time or space. For example:\n",
        "\n",
        "- The number of meteorites greater than 1 meter diameter that strike Earth in a year\n",
        "- The number of patients arriving in an emergency room between 10 and 11 pm\n",
        "\n",
        "The probability mass function is $$f(k;\\lambda) = \\frac{\\lambda^k exp(-\\lambda)}{k!}$$\n",
        "For large k the normal distribution is an excellent approximation of the poisson p.d.f. For k below 20 one should be careful using statements based on the normal distribution.\n",
        "\n",
        "The standard deviation of the Poisson pdf is simply $\\sqrt{\\lambda}$. This is very convenient, in particular for higher $\\lambda$ where the normal pdf is a good approximation and the probabilistic interpretation is easy (see below)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rfa41OH3dKUk"
      },
      "source": [
        "**Plot the Poisson**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dYggRMKydKUl"
      },
      "source": [
        "x = np.arange(-1, 40)\n",
        "plt.plot(x,scipy.stats.poisson.pmf(x,2),drawstyle='steps-mid')\n",
        "plt.plot(x,scipy.stats.poisson.pmf(x,15),drawstyle='steps-mid')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FjE0O77tdKUl"
      },
      "source": [
        "## Annex 2 The Binomial PDF\n",
        "\n",
        "The binomial pdf models the probability of getting k positives when drawing n times, each time asking positive or false. For example,   \n",
        "\n",
        "$$f(k;n,p) = \\binom{n}{k} \\cdot p^kq^{n-k}$$\n",
        "\n",
        "The binomial distribution converges towards the Poisson distribution as the number of trials goes to infinity while the product np remains fixed or at least p tends to zero. Therefore, the Poisson distribution with parameter λ = np can be used as an approximation to B(n, p) of the binomial distribution if n is sufficiently large and p is sufficiently small. According to two rules of thumb, this approximation is good if n ≥ 20 and p ≤ 0.05, or if n ≥ 100 and np ≤ 10.\n",
        "\n",
        "For n>20 and p not too close to 1 or 0, the normal distribution is also here a good approximation. \n",
        "\n",
        "**Plot the binomial pdf for various parameters**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PYzb7ONudKU1"
      },
      "source": [
        "## Annex 3 Probability\n",
        "\n",
        "An abstract definition of probability can be given by considering a set $S$, called the sample space, and possible subsets $A,B,...$ the interpretation of which is left open. The probability $P$ is a real-valued function defined by the following axioms due to Kolmogorov (1933) [9]:\n",
        "\n",
        "- For every subset $A$ in $S$, $P(A) ≥ 0$;\n",
        "- For disjoint subsets (i.e., $A ∩ B = ∅$), $P(A ∪ B) = P(A) + P(B)$; \n",
        "- $P(S)=1$.\n",
        "\n",
        "From this further properties can be derives, e.g.\n",
        "\n",
        "- $P(\\bar{A}) = 1 - P(A)$\n",
        "- $P(A \\cup \\bar{A}) = 1$\n",
        "- $P(\\emptyset) = 0$\n",
        "- if A in B, then $P(A)\\leq P(B)$\n",
        "- $P(A \\cup \\bar{A}) = P(A) + P(B) - P(A\\cap B)$\n",
        "\n",
        "#### Conditional probability\n",
        "In addition, one defines the conditional probability $P(A|B)$ (read as $P$ of $A$ given $B$) as $$P(A|B) = \\frac{P(A ∩ B)}{P(B)}$$\n",
        "\n",
        "As an example, when throwing the dice, consider obtaining more than 3 eyes given only trows with even number of eyes outcomes. We calculate the (conditional) probability:\n",
        "\n",
        "$$P(n>3|n\\; even) = \\frac{P(n>3 \\cap n\\; even)}{P(even)} = \\frac{2/6}{3/6} = \\frac{2}{3}$$\n",
        "\n",
        "#### Independence\n",
        "\n",
        "If A and B are independent, then \n",
        "\n",
        "$$P(A|B) = \\frac{P(A ∩ B)}{P(B)} = \\frac{P(A)P(B)}{P(B)} = P(A)$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I6nnbYcJdKU1"
      },
      "source": [
        "### Some words on probability and bayesian versus frequentist statistics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9Du390uOdKU1"
      },
      "source": [
        "Since data taking is data taking of random variables, we need to define and talk about probability. In mathematichs probability is defined in a rather abstract manner (see Annex below). For our purposes we go directly to the interpreation as either **relative frequency** or **subjective probability**. If A is a possible outcome of an experiment repeated n times, then the probability of A is the realtive frequency\n",
        "\n",
        "$$P(A) = \\lim_{n\\rightarrow \\infty} \\frac{times \\; outcome \\;is\\;A}{n}$$\n",
        "\n",
        "The subjective probability is\n",
        "\n",
        "$$ P(A) = degree\\;of\\;belief\\;that \\;A \\; is\\; true$$\n",
        "\n",
        "Both concepts are consistent with the abstract mathematical definition .\n",
        "\n",
        "#### Bayes' theorem\n",
        "\n",
        "From this definition and using the fact that $A ∩ B$ and $B ∩ A$ (intersection) are the same, one obtains\n",
        "Bayes’ theorem\n",
        "\n",
        "$$P (A|B) = \\frac{P(B|A)P(A)}{P(B)}$$\n",
        "\n",
        "first published by the Reverend Thomas Bayes (1702-1761). Statistics based on the relative frequency interpretation of probability is called frequentist statistics, on bayesian theorem bayesian statistics."
      ]
    }
  ]
}