lectures/Lecture_7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Classical/Frequentist Statistical Inference: III\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "Material in this lecture and notebook is based upon the \"Maximum Likelihood and Applications in Astronomy\" lectures of A. Connolly's & Ž. Ivezić's \"Astrostatistics & Machine Learning\" class at the University of Washington (ASTR 598, https://github.com/dirac-institute/uw-astr598-w18). Also the \"Inference\" and \"Histograms\" lectures of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), J. Bovy's mini-course on \"Statistics & Inference in Astrophysics\" at the University of Toronto (http://astro.utoronto.ca/~bovy/teaching.html). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 4.\n",
    "- [David Hogg's \"Fitting A Model To Data\"](https://arxiv.org/abs/1008.4686)\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "* [Hypothesis Testing](#one)\n",
    "* [Comparing Distributions](#two)\n",
    "* [Nonparametric Modeling & Histograms](#three)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis Testing <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "A common problem in statistics is **hypothesis testing**, i.e., testing whether an observation or observations agree with a certain hypothesis. \n",
    "\n",
    "The simplest example is whether a measured value $x_i$ or an entire dataset $\\{x_i\\}$ are consistent with having been drawn from a Gaussian distribution, $\\mathcal{N}(\\mu,\\sigma)$. In this example, $\\mathcal{N}(\\mu,\\sigma)$ is the ***null hypothesis***, and we are trying to reject it in order to claim that we have measured some effect. \n",
    "\n",
    "A great example in astronomy is ***source detection*** in images that suffer from a substantial noise background (e.g., atmospheric fluctuations in optical images).\n",
    "- Since the background fluctuates, we need the source flux in a given image resolution element (e.g., pixel) to be substantially larger than the background variation in order to give a robust detection.\n",
    "- Here the null hypothesis is that the flux in each pixel is caused by background fluctuations, and if we can reject it then we can achieve source detection. \n",
    "- But note that *rejecting the null hypothesis is not the same as favoring the source hypothesis*! Our understanding of the background could be poor or incomplete, rendering our model of the null hypothesis similarly lacking.\n",
    "\n",
    "In the image below (from the [photutils](https://photutils.readthedocs.io/en/stable/background.html) package) we have sources embedded in Gaussian image noise. \n",
    "- Ideally we would like to be able to find the real sources in here and assign a number that says how much they reject the null background hypothesis. \n",
    "- We would establish a threshold in these numbers above which we could accept sources, knowing that we might be falsely rejecting some real ones, but hopefully not accepting too many spurious noise sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![https://photutils.readthedocs.io/en/stable//background-1.svg](https://photutils.readthedocs.io/en/stable//background-1.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine we have data $\\{x_i\\}$. This could be signal strength in time-domain data, or pixel brightness in image data. Under an assumption about the null background distribution, $h_0(x)$, we can compute the probability that a data point would have a value as large or greater than that observed. This is called the ***$p$-value*** of the data:\n",
    "\n",
    "$$ p_i = \\int_{x_i}^\\infty h_0(x)dx \\equiv 1 - \\int_{-\\infty}^{x_i} h_0(x)dx = 1-H_0(x_i),$$\n",
    "\n",
    "where $H_0(x)$ is the CDF (cumulative density function) of the null distribution.\n",
    "\n",
    "In the image below, we show the distribution of a statistic $S$ under a null hypothesis $\\mathcal{H}_\\mathrm{null}$. The $p$-value of the observed statistic $S_\\mathrm{obs}$ is shown graphically.\n",
    "\n",
    "<img src=\"figures/freq_pvalue.png\" alt=\"\" style=\"width: 500px;\"/>\n",
    "\n",
    "Typically, a threshold $p$-value is adopted, called the **significance, $\\alpha$**. In a classification scheme, data would then reject the null hypothesis if $p_i\\leq\\alpha$. If $\\alpha=0.05$, such data would reject the null hypothesis at the $0.05$ significance level. \n",
    "\n",
    "**NOTE**\n",
    "- If we can't reject the null hypothesis then it doesn't mean the data belongs to the background. It could be that our sample is simply not large enough to detect an effect. \n",
    "- Rejecting the null hypothesis does not necessarily mean that we've proven a signal hypothesis. Our model of the null hypothesis could be bad, or there could be many other things we have not considered.\n",
    "- There is ***alot*** of mis-use and abuse of p-values in the scientific literature. Take a quick scan of [this](https://fivethirtyeight.com/features/statisticians-found-one-thing-they-can-agree-on-its-time-to-stop-misusing-p-values/) and [this](https://www.wikiwand.com/en/Misuse_of_p-values) during class, and read thoroughly in your own time. \n",
    "\n",
    "**EXAMPLE**\n",
    "\n",
    "Let's flip a coin 10 times. \n",
    "- Defining HEADS as a success event, use `scipy.stats` to define a binomial distribution over $10$ coin flips, where the intrinsic success probability is $0.5$. \n",
    "- Plot the probability mass function (the analog of the probability density function for discrete distributions).\n",
    "- <font color='red'>If we get $8$ HEADS, what is the $p$-value with which we reject the hypothesis that the coin is fair?</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import uniform\n",
    "from scipy import optimize\n",
    "from astroML import stats as astroMLstats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code to plot the binomial distribution in here.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code to compute the p-value of 8 tails in here.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We look at our coin more shrewdly. In 20 flips, we get a result with $p$-value of $1$ in $10^3$ (i.e., $p = 0.001$). <font color='red'>How many HEADS did we get?</font> *Hint: check the properties of the `scipy.stats.binom` object for useful functions for this.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When making these kinds of classification decisions, we make two types of errors. These are quite literally described as **Type I** and **Type II** errors by *Neyman & Pearson*.\n",
    "\n",
    "**TYPE I ERRORS**\n",
    "- The null hypothesis is true, but incorrectly rejected.\n",
    "- These are ***false positives*** (or false alarms).\n",
    "- False positive probability is dictated by the significance level $\\alpha$. \n",
    "\n",
    "**TYPE II ERRORS**\n",
    "- The null hypothesis is false, but not rejected.\n",
    "- These are ***false negatives*** (or false dismissals).\n",
    "- False negatives probability is dictated by a variable called $\\beta$, related to $(1-\\beta)$, called the ***detection probability***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at Fig 4.5 of Ivezic, v2. We have two Gaussian distributions; one describing the background noise and one describing the source. These distributions overlap, as they often do in real life.\n",
    "\n",
    "<font color='red'>Execute the following cell. We'll discuss it more below.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ./code/fig_classification_example.py\n",
    "\"\"\"\n",
    "Example of classification\n",
    "-------------------------\n",
    "Figure 4.5.\n",
    "\n",
    "An example of a simple classification problem between two Gaussian\n",
    "distributions. Given a value of x, we need to assign that measurement to one\n",
    "of the two distributions (background vs. source). The cut at xc = 120 leads\n",
    "to very few Type II errors (i.e., false negatives: points from the distribution\n",
    "hS with x < xc being classified as background), but this comes at the cost of\n",
    "a significant number of Type I errors (i.e., false positives: points from the\n",
    "distribution :math:`h_B` with x > xc being classified as sources).\n",
    "\"\"\"\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Generate and draw the curves\n",
    "x = np.linspace(50, 200, 1000)\n",
    "p1 = 0.9 * norm(100, 10).pdf(x)\n",
    "p2 = 0.1 * norm(150, 12).pdf(x)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 7))\n",
    "ax.fill(x, p1, ec='k', \n",
    "        fc='#AAAAAA', alpha=0.5)\n",
    "ax.fill(x, p2, '-k', \n",
    "        fc='#AAAAAA', alpha=0.5)\n",
    "\n",
    "ax.plot([120, 120], \n",
    "        [0.0, 0.04], '--k')\n",
    "\n",
    "ax.text(100, 0.036, r'$h_B(x)$', ha='center', \n",
    "        va='bottom', fontsize=15)\n",
    "ax.text(150, 0.0035, r'$h_S(x)$', ha='center', \n",
    "        va='bottom', fontsize=15)\n",
    "ax.text(122, 0.039, r'$x_c=120$', ha='left', \n",
    "        va='top', fontsize=15)\n",
    "ax.text(125, 0.01, \n",
    "        r'$(x > x_c\\ {\\rm classified\\ as\\ sources})$', \n",
    "        fontsize=15)\n",
    "\n",
    "ax.set_xlim(50, 200)\n",
    "ax.set_ylim(0, 0.04)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$p(x)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we have data $\\{x_i\\}$ drawn from an underlying distribution\n",
    "\n",
    "$$ h(x) = (1-a)h_B(x) + ah_S(x)$$\n",
    "\n",
    "where $h_B(x)=\\mathcal{N}(\\mu=100,\\sigma=10)$ is the background distribution (normalized to unity), $h_S(x)=\\mathcal{N}(\\mu=150,\\sigma=12)$ is the source distribution (normalized to unity), and $a$ is the relative normalization factor that accounts for background noise being $(1-a)/a$ more probable than sources.\n",
    "\n",
    "- ***If we set a classificiation threshold at $x_c=120$, calculate the Type I error (false positive) probability, $\\alpha$.***\n",
    "- ***For the same $x_c$, calculate the Type II error (false negative) probability, $\\beta$.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a sample of size $N$ (containing background noise and sources), the **expected number of spurious sources (Type I / false positives)** is \n",
    "\n",
    "$$ n_\\mathrm{spurious} = N(1-a)\\alpha = N(1-a)\\int_{x_c}^\\infty h_B(x)dx$$ \n",
    "\n",
    "and the **expected number of missed sources (Type II / false negatives)** is\n",
    "\n",
    "$$ n_\\mathrm{missed} = Na\\beta = Na\\int_0^{x_c}h_S(x)dx.$$\n",
    "\n",
    "The **total number of classified sources** (that is number of instances where we reject the null hypothesis) is\n",
    "\n",
    "$$ n_\\mathrm{source} = Na - n_\\mathrm{missed} + n_\\mathrm{spurious} = N[(1-\\beta)a + (1-a)\\alpha].$$\n",
    "\n",
    "The **sample completeness** (or **detection probability**) is defined as\n",
    "\n",
    "$$ \\eta = \\frac{Na - n_\\mathrm{missed}}{Na} = 1-\\int_0^{x_c}h_S(x)dx = 1-\\beta$$\n",
    "\n",
    "Finally, the **sample contamination** is\n",
    "\n",
    "$$ \\epsilon = \\frac{n_\\mathrm{spurious}}{n_\\mathrm{source}}$$\n",
    "\n",
    "where $(1-\\epsilon)$ is sometimes called the **classification efficiency**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing Distributions <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "Hypothesis testing is strongly linked to the **comparison of distributions**. Common questions are whether two samples are drawn from the same distribution, or whether a sample is drawn from some known distribution (like a Gaussian).\n",
    "\n",
    "Sometimes we may be only interested in the location parameter of some underlying distribution. For example, do two gravitational-wave detectors with different measurement uncertainties and local effects measure consistent masses of a merging system of black holes?\n",
    "\n",
    "There is enormous variety in the nonparameteric and parameteric statistics to compare distributions, and we'll only look at a few. See the textbook and wider literature for further examples.\n",
    "\n",
    "### Kolmogorov-Smirnov (KS) Test\n",
    "\n",
    "By far the most popular nonparametric statistic to compare distribution is the KS test. We compute the empirical cdf $F(x)$ for two samples $\\{x1_i\\} = (1,\\ldots,N_1)$ and $\\{x2_i\\} = (1,\\ldots,N_2)$. *Recall that the empirical cdf can be calculated by sorting the samples and dividing by the sample size.*\n",
    "\n",
    "The KS test is based on a statistic that measures the maximum difference in the cdfs of our two samples\n",
    "\n",
    "$$ D = \\mathrm{max}|F_1(x1) - F_2(x2)|$$\n",
    "\n",
    "where $0\\leq D\\leq 1$. \n",
    "\n",
    "> The core idea is to know how often the value of $D$ computed from our data would arise by chance if both samples were drawn from the same distribution (which is our *null hypothesis* here). \n",
    "\n",
    "Amazingly, ***this does not depend on the underlying distribution we care about!*** Kolmogorov published the following result in 1933 and Smirnov published numerical tables in 1948. The probability of obtaining a value of $D$ larger than that observed in our data is\n",
    "\n",
    "$$ Q_\\mathrm{KS}(\\lambda) = 2\\sum_{k=1}^\\infty (-1)^{k-1}\\exp(-2k^2\\lambda^2)$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\\lambda = \\left(0.12 + \\sqrt{n_e} + \\frac{0.11}{\\sqrt{n_e}} \\right)D $$\n",
    "\n",
    "and the ***effective number of data points***, $n_e$, is\n",
    "\n",
    "$$ n_e = \\left( \\frac{1}{N_1} + \\frac{1}{N_2}\\right)^{-1} = \\frac{N_1 N_2}{N_1 + N_2}$$\n",
    "\n",
    "If the probability that $D$ were drawn by chance is very small (low $p$-value) then we can reject the null hypothesis that the two samples were drawn from the same distribution."
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAACECAYAAAA9dKxTAAAgAElEQVR4Ae2dCdAetXn4mbZJk7STphe0aRqadNomlEyn9EjbtEynpTSdpNA/EAMGm9OYyzH3ZcyNQ8DhMOYwGDvcYDBgbozB2MY2GDDG+OAyNjFgiIFwXzn0n9+aZy3L2l1pX+1+7/d+j2be2eOVtNqfnpX0SI+kTYw6JaAElIASUAJKQAkoASWgBJRAAYFNCu7rbSWgBJSAElACSkAJKAEloASUgFGFQYVACSgBJaAElIASUAJKQAkogUICqjAUotE/lIASUAJKQAkoASWgBJSAElCFQWVACSgBJaAElIASUAJKQAkogUICqjAUotE/lIASUAJKQAkoASWgBJSAElCFQWVACSgBJaAElIASUAJKQAkogUICqjAUotE/lIASUAJKQAkoASWgBJSAElCFQWVACSgBJaAElIASUAJKQAkogUICqjAUotE/lIASUAJKQAkoASWgBJSAEohWGKZMmWL0pwxUBlQGVAZUBlQGVAZUBlQGVAb6hwx0qvKowqAKkCqAKgMqAyoDKgMqAyoDKgMqAz0sA32mMHT64IEWXjTwgfbeqd5X+cWRVF5xvFzfys8lUu9aOdbjFhtKOccS8/tXjn4uTd9V7s0STsW39ghDs6/Xe7GnyrDeIxP2RsovjJP4Ul5Cot5R+dXj5oZSji6RZq6VcxquyjENx9hYlHsssTj/qfiqwhDHvbbvVBlWOwH9PKDyi8tA5RXHy/Wt/Fwi9a6VYz1usaGUcywxv3/l6OfS9F3l3izhVHxVYWg2n/LYU2VYHuEAO1F+cRmuvOJ4ub6Vn0uk3rVyrMctNpRyjiXm968c/VyavqvcmyWciq8qDM3mUx57qgzLIxxgJ8ovLsOVVxwv17fyc4nUu1aO9bjFhlLOscT8/pWjn0vTd5V7s4RT8VWFodl8ymNPlWF5hAPsRPnFZbjyiuPl+lZ+LpF618qxHrfYUMo5lpjfv3L0c2n6rnJvlnAqvqowNJtPeeypMiyPcICdKL+4DFdecbxc38rPJVLvWjnW4xYbSjnHEvP7V45+Lk3fVe7NEk7FVxWGZvMpjz1VhuURDrAT5ReX4corjpfrW/m5ROpdK8d63GJDKedYYn7/ytHPpem7yr1Zwqn4qsLQbD7lsafKsDzCAXai/OIyXHnF8XJ9Kz+XSL1r5ViPW2wo5RxLzO9fOfq5NH1XuTdLOBVfVRiazac89lQZlkc4wE6UX1yGK684Xq5v5ecSqXetHOtxiw2lnGOJ+f0rRz+Xpu8q92YJp+KrCkOz+ZTHnirD8ggH2Emv83tq+VPm2COPSZarvc4rGaiCiAY6v9NOOtU8+sijBXTCbw90juGk/D4//OBDs9eQPc0HH3zg9/DJXeVciif4T+UYjKrSY0ydptwrcXbkIRVfVRg6yobwwLEZ9otf/MLwwc1+YJaZMX2GeWLRE+btt94Kf2CAz3fffdc8/dTT5v4Z95tZM2eZZUuXmXffeScg5MZeXli1ysyfO99Mv+se8+iCR8wra17Z2FMHd7qR369+9Ssz876Z5pmnn6n9ZsQBq5NGn2hGH3t87XjcgL3Cq43vwGXHdTfyI11N8/j4o4/Ng7PnmB2328E8ufhJH5qoe73O8fXXXzeLn1hspt893cybO888+8yz5uOPP45iVOT57bffNpdP+nGWF5QTZa7XOacoa8v4yX+9zlHec+3atVmHAPU1cvv8iuflr46Pdeo05R6HffXq1eb2W28PDhTLtyhiVRiKyCS+H5NhM6bfa4bvvZ/ZafsdzYgDDjYjD/x+dr7LTjubiy+4qLK3qSrpL730kjnjtB9kFdGg//c9c9B+B2bnNBLWPeNi8+orr1ZFk/2/fOkyc/ThR2fh99trmDnikMOzOEj7aSefZl568aWgeKo8dRM/0vqzN35mTjr+xOy9p1w3xZv8JU8uMUN22b30973/2ylnf/rJp3njqXOzP/Jy37Pp78B9nn3dbfxIW6c8Tjnx5FJZHDxo11wWKQto/HbqepEjTBY9vsgcOuLQjBff+L577JOz22fo3uaG628w77//vhcfSl9VubDLjjvn8Q3+3q7eeOybvcqZdwwpa20WnZz3Mke4oCiMPeOsrD0xdNch5shDjzRDBw/NZA15fuzRx0rxNVWnKfdy7nam3D/jPkOZQBsr1MXwLYtTFYYyOgn/C82wu+64K/t4URTotRf34uoXM8WBihxzASqdOm7BQw8blAQqOM5lqBsFgUqO+PkN22tf8/prr5c+4tlnnjE0MhDeh+c/nPtl5OKHY87I4hm+z36GXrhOXbfw4z1oLOw9ZK+cVZHCsHjR4tyPcC07jv3h2E4x5eH7I6888caYpr8D+1m+827iR/pS8Bh1zKgoeVz9k9U+NFH3epHjFZMvzzgyIkiHCWUxvaooWGef+aOc8SknnOwtp3/5y1/mfsrKA/kPk6Qq14uceefQsraKT+j/vcqR92fE6vsHjshk7/prrstlE9m94bop2f2ddxiUjZgV8WqqTlPu5dzJD8wTx51zXl52qMJQJKU9cD/kg5j34NxMa6TX+blnn9vorTEfkkpk/Hnnb/R/1Q1GFnYbNDhr4PviJ/y5Y8/Jn8HIQZFj5GDP3ffM/N427daNvNG7Jo3qw75/qHn/PX9v20YBC250Az8q+quvvDrLo/33HZ5zuv7a672pxlSB/Jo65UYze9ZsQ/6ipNGLQ0W49MklmdkZDY2Vz680r732mjeeOjf7Iy95z6a/A3lO2bEb+En6UvE4/phR5qjDjjQP3P+AeXDOg+ahefPNIwseMY8vfDxrJGCSiHndiudWZJ0VP//5zyUJtY+9xnHOrHXmWnTo+Mw3P/roI3PQ8PUjtigXrqOBRrlA5Y/J6dw5c81D8x/KTEQWLVyUmYItX7bc0CGDqUjIKG2vcY4ta13Gda97jaNwoLF5zBHrLAGwLvA5OiKRy9132d2sWrm+s9L221SdptzLuZMflDnkjYwIkVehLoRvSFw6whBCKYGfqgyjUkC7Rwiu/PEVhU88+6yzMz/4u+/e+wr9+f647upr8wLhrTf98yEYsiRu+TEc7DoK8wP3OyDzc+xRx2a9a64frrGPlHioHDtxfc2Pxjw9tHvsNjSbt8CEUHm3IoWB4Vv8pOipjWXXH3nxjm18ByEs+5qfpDElj9HHjTYpR7EkjWXHXuPIqALfNOVekUMhk7KBHl2f439GjVK5XuJcp6xVjuXtgfPPHZfJJGZIb7zxhhfXT174iRET2QP23d9brzdVp/WS/NpwU3CnHYWpOIoeVh9HH35Ulpc6wmCT7rHzqg/i8k+GualIyuYPoGlKZRQ7SXbSpZflYelNLHJosfIMJjC7zm4sY09X5FAsGE4nLsyW6OWo6/qaH5MaUdbefPPN7BUWPvpYzqhIYWAEgXfvRoWhaXmrwwuwTacrVP76Wt4knSl5nNCFCkPK9xNmvmOq5xw+8rDsm8ask5ECn2NUQMpPKnXf6Cr/t6kwpHp/3/va91I8p27ZYaej7nkvfvfvvPNO3hk57pxxpWjOHPPDXHZ9ix40Vacpdz93RiyPO+rYzDpBMk7mjFKGhLoqvqHx6AhDKKkO/ZVlGDawTJRDAA4Ytn/pk6ikpEFPZVQ1z8CObM3LawwTay8af1Fuv2j/L+eYEEmF5xvFOPMHZ+b/r/3pWgnmPY45dUzul+H8uq4b+Nlpx4xDGGEP6nMoZfjpNoWhLXmzmYTw6ot02Wm0z7tB3lLzYJJ+N40wpH4/O//s85TPwWSIkcY7brvDfsQG55gqSdnA0bdiHGV3WwpDyvff4EWdi6aeE1J2OEmpfdmL3z2yKvL4wP0zS9lMu3la7pcFVlzXVJ2m3Mu52/nA8uuSn/b9svMyvmXh3P9aURho5FKYYBO7fPny7MfSc9wbKK4swxY8vCAXgAvPv6ASiazOg9D45g9URlDhYZ891ikvxL/iuQ3nUrC0K71r/MfqSlVOJlPhv8h2sioO/u82fjQceCd+RQoDdsj8320KQ1/IWwivvkhXkex1g7yl5nHy6JO6SmFI/X5FednWc+T5mJFJ2UDnjm8kgjK0LYWhrfdv6jkhZYew7/TYi989KxeKPFbNkxOFAP977rbHRm20puo05b6uc7GIuy3XmENKftr3y87L+JaFc/9rXGFAMTj++OPN9ttvb/74j//YbLLJJtlv8803N8OHDzdLlixx09ST12UZdunFl+QCcOstG08gdoFcfcVVuf8fnZluZR2es2bNmjzuIw89YqPKjgmYIqw/OHWMm7SNru3CnknQdV238WNfDOHA3BCfYx8N/HSbwtAX8hbCqy/S5cs37nWDvKXmwao93TTCkPr9ivKyrefI8xmVlbJh8sTJcnuDY5sKQ1vv39RzQsqODeB2cNFr3z1zFUUWd995t0oyH374YT6PIau7Vm+4UlpTdZpyL+duZ9yoo4/L89TXGWH7lfMyvuIn5NiYwoD9+mOPPWb+/u//PlMQtt12WzNmzBjzwgsvmKlTp5rvfve72f0//MM/NLfddtsGaSXsqFGjzMqVKze4358vyjKMhrd81GwEVuXuuuPO3H/sPIaquH982eQ8bpZQcx0KjaQ1ZKUmJlKJf4bhyds6rtv4yWoRvFuRwiCrWqEwMC+FFVFuuvEmw1ySaTfdkq2UJHMi6jApC9MfeXXTd9AN/FLzOPXEUzKFgU4cGmH0cNP5wEo+TNRlGefUo769yLHsu+M/1ranXGCCKUtZ+hz7LMAf8yVWSGKkmLKX0UquKS9CGwPE38ucQ8paH+M693qNIyueSf3LROYQJ6sfEs6dx9BUnabcTb7qpI+7nW+sdid5GlpGlPG14646b0RhoFF41lln5aMJO+yww0YV0SuvvGK22mqrzM+XvvSlTLmQxE6aNCm7v2xZ8cRc8dtfjmUZdvjI9UOGLHNY5diVWQRmxP4HV3kP/v+FVS8Y2TDoqsuv9IabPHFS/mwquCrHEKiklWPRCg1V8XQbPwpSea8ihYElKvEjS7DSw8PSaO4GWeedfW7hJk9VXIr+74+8uuU7gGk38EvNg2UTMTNgaWWUdzaHpBEhcsyRUUBsxlO5XuRYxubee6ZnPFltpmwTLFY9Ya8b/FHmYt4py1BLfrAELvPOQlwvcw4pa0MYhfjpNY4s4y3yhMVAiLM3cnXnHTZVpyl3s8EGui53O99kpTbytd8rDLzAokWLSpUFefmLLroo9zdixIjs5Qm/yy67mC222MIwPNYrruyDoBKXj5pVCKocSoX4x0Y2hYO1rADCTpBFgnjWGesnPBetDmSnh/0YJK0c625B3238ZLUI3unaApMkmMKLXlx6b4UpvbissW4P4R+8/0Hm5ZdettF1dN4feXXDdyDQu4Ffah6seMY8hvlz5xs2VxTH0smYGrJTO/KMMlGkBEuY0GMvcix6d+YuoIzBsGxSNOEpM+goWLpkqWHERxxx3HHb7fk8MToZQjqReplzSFkr/Do99hrHO29fb43AKmkhDsUCGebnzpFsqk5T7oxMFnO38418lPyRNoX9v++8jK/Pf9G95CMMLAP1jW98I1MEPv/5z2fKQ9HDbcVi0003NYR99dVXze/93u9l8x5SD48XpaON+2UZxpKjIgAhvXt2jwv7AnTqmIx+yoknZ2lgh+aPP1pfgblxs8qSpLVoFMIOg0CLf4727tW2v6rzbuNnTw4rUhiq3on/baUBG/NUrj/y6uvvwGbfDfza5rF69epsBEK+V+yVO3UDhSOrxTFiA7upN0ztFFu2sR6KG/FhIvLee++VxtnLnFOVtaUAP/mz1zjectMtef3L6johzu7BrlJ8ffHVqdOUuzGh3O1Fb/q1wkAD/5RTTslHDZjsXPZCNFR/8zd/M/dPhTVu3Djzmc98pucmQ5d9ECylKpU0vX9Vzt4HgSVQO3HkAUoCz59w4YTKOQbjzxufp3XihEsrH83eC/JuHH27o1ZGUmEi0hf8ln+yZCrv1InCwDczxNr3gv0dUrhuk7cQXn2Rj0Wsu4FfX/CwGxhlO70XcXPvDwSOLG2NaSgTmUPmoLmMiq7tzpmyzTwJ38ucQ8qOIoax93uNI8uoSv176IhDgnBgCidhmE8T6+rUacrdmFDuJ40+Mc+fsva1nW9lfG1/VedJRxhWrVqVN/5ZDWn69Omlz2dE4Stf+UoeZt68eWbLLbc0Bx54YLYEa2ngfvZnWYbZ6+qyaU2Vs+cwUKnUdfCX7eBZ/jTEMaoghck5Y8+uDGLPYaDHtK7rNn6yWkSmMFx1Td3XysKdO/acnGnIsrohD+uPvPrqO/Dx7AZ+fcEDszj5vjlWLcPoY2ff63WO7LGAYkfZtvCxhfard3wu8yHIh6olrHuZc8qytipTeo2jvYfFvnvsU/X62f/2HIZnn3k2KIzrKbZOU+4bzmEo444lgpTRoYvIlPF1867sOqnCMG3atLzxj8JQNQcBTXTrrbfOw1x33XXma1/7mnnppZfK0twv/yvLsDNOX9fDjxCETCS27RJ9m6uEAKLnnyEwhr1jhh1vm3ZbLqwhk6iYSC3C3ckE7W7j9/QnS6bybtd2qDDYK0+xq2MK1x959cV3UMS6G/j1FQ/s5uWbXfT4oiJEQfd7meOLq180w/bcNxshZH361M5e4YZyusxUtJc5pyxrq/Ko1ziufH79ruN8075dx10m9typuouUxNZpyt1sYA5axl3Mx8nPfqswMHHLNkf65je/WfkyrsLwrW99y0yePLnnRhf4IMs+CJbSkwo6ZMSAZTnF/4Oz43dPpuLBDo4VOu6fcb9bXpReL3lySf5sVvupGhJ7eP76VRouvuDi0rjL/uwmfqRTVosgH6658uqypFf+hxmD5Cd20Clcf+TV9ndQxrkb+PUVD3vlJN9O72Xc3P96lSPLnrLKEfMLaNg34X766k/zcoHywbdjtDy3VznzfinLWuFVdOw1jnTa2nOhynquYYLJsMydYUW/ui62TlPu4dzFKoQygTZ0iCvjGxJe/CQbYSDhgwYNykcLWPWoyrkKw9ChQ3tSWYBDWYaxTr80GKmAqoRA1vlmhSTMimIcDXzWd8fedt7ceTFBM7+Et3eCXvl8+V4Zl0/6cf5unUyi7BZ+AoyCV/KsSGHAnOOeu+42rEJT5lgVReI65oijy7wG/9cfebX5HVSB7AZ+KXnw3WKPXLbUpzAZOnhoLo+dmtn0GkcYsRnWgfsdkI0uMO8u1jG6iyJWVXauWrkqzwfKh7Kyvhc5C9eQslb8dnrsRY5s1ij1y+233l6KiJXUxO/UKTdu5LepOk25l3O3M8Ke21TVVpRwZXzFT8gxqcLwne98J1cYLrnkksrnM5ximyTdfffdlWH6q4eqDDvk4JH5h8ruyEXugw8+yHdirGPvzogCBQKbBlU5JkSzoZPrLr340jytLBla5pg4yfMOGn5gmbfK/7qFnyTUtqstUhhk1ZQzx/xQgnmPF4xbP5E8ZDM8byTOzf7Ii1do6ztwcG102S38UvGYMX1G/s2WNVTXrl2b++O7ZdW6TlyvcYTF+eeOy/ZOCFkimlECNsmzHcupwhbTj7IRWspe/PGr2nSrFzkLs5CyVvx2euxFjnQMihwdW2HyyuR6/DLKwMpfrmuqTlPu5dztfLDnMPRrhWHYsGG5wnDzzTfb7+g9p8fEVhjuuecer79euFn1Qdw45cb8ox53znmFr8ykaPmgse10HQLEygis7+067Bf3Gbq3YcfXEIe5E5sJuQ57XSmAqMiKhJYeUsye8HvzjTe50URddwM/O8GwEQY+hQHZlqHgMoUBdgz9SlwMv6dw/Y2XvHOq70Diq3vsFn6peLBnisgYPddFzm6kMvTdqes1jux0S2MqdOlUzDB/dObYDTDKOupVS2LbHQlV5WevcbaBVZW1tt9Oz3uRI2ZJWCPI979mjX8zQEyVZWU237ffZJ2m3NetlOnj7sq03V4oanu5Yar4uv6LrpONMPCAyy67LFcYxowZU/TM7P6DDz5o/vd//9ewVwMTpO1Vleh1oXe7rPelNPIu/LMqw8h40RypkHymAwyFU8nw4U+7eZr3LS8af1FeMDB/wHayxBqVGIWwmCgsXrTYsHQdCgi/xU8sNqyahDJAAeJztn21b8k/8k62MGfyZujkHN+zuNcN/Oy00csoBXDRJleMFuCHZWuLHBOmJZ5LLppQ5C36fn/kxUum+g6igTkBuoVfKh7Y3FOu8GMipM9hOrfXkD0zeaRjoUyx8IX33es1jrKsNKYdc+fMNeyki9kWe+PQEy5lKJutTZ44Oes0IIztxL6bjfIwT/I5RiXIK8oGzBSrJqv2GmebSUhZa/vv5LxXObJ4AWbIyNOoo4/zdvJJXcQqSW+++aYXY1N1mnJftxJaEXfJDFaxk3wkL0PbyFV8Jf6qY1KFAU12xx13zBr/bL728MMbNlhpNGJSc9JJJ5nPfvaz2ejCnXfemSsMF154YZbeq6++OruH0tArLiTDYCNr8aIYoDSIBvnOO++Y0z7ZNG3yxEmFWFi1Qxqg7j4Jl0++PP9P/FQdy9ZuZt8GwlOxYYcvmwuRZnrd+Y9Vf8psbwtfxPmjG/iRJGyWbWWJd0RJo4Ewe9bsPL/wi72n5OfVV15tyENxmH7Qa8gIDPzcHTXFX91jf+Ql75riO5C46h67hR/pT8WD1dCobEYeNDLrFBA2fK8sv8i3ijyz47vPHEH8xxx7jaN8z3AK/bnlMAuEjDtnXYcDnSn2qoAoBiybzRKYxM/S1WWrI0le9Bpn3iumrBUOnR57kaMwmTNrTq6EMnplf+MovjvvMCjrMFjzsn8EgniaqtOU+56mjDtta/LI3g16l512zspt9u567tnnJJu9xxC+3oDOzaQKA3HzYn/7t3+bKwE77bRT1jt82GGHmW222cZsttlmZvPNNzfMcaAhiVJwwAEHZP6//vWvm7POOis7HzlyZMe90s679ullaIYxinDmD87MTXlQHOhhYkWi3QYNNjQ8y7RKe4TBXRIRU6TQSk78lU3CJR1XTL7cyDKMFDisHY8ZE41gJlu9/fbbSbh3Az9exJ4QKozsoyhN8tI0Dhg5IP/4wNloj48ePjTejjjkcMNGfKldf+UlHDr9DiSeusdu4SfpT8UDc8LvHzgiKwfo4aa3UZZR5LulBxEFJZXrJY6Ud7vutEt0GXrl5Vd6cU6/656sgUZZwJLTKGy77LhzFj+9vKFmT0TeS5wFVmxZK+E6OfYiR5sHVgbD91m3Izn1FvWP7LtAXR+y4lcTdZpyL19pzV4YxW5vyHlZJzL5H8rXlhXfeXKFgYegCIwfP94MGTLEbLHFFtluzl/96lezEYVJkyYZbOhsE5UVK1YYdoX+4he/mO3DwApLvTS6UCfD6IFGc6TnmeHvBQ8vCNolmd5ChspDPnyfQNS59/7772eN3rvvvNuwUyzPxwQipYsV+G7jR2ODHjN6D/n4MV9IMfJSxLi/85L3qpuPEr7usVv5peLB8omYIlK+zHtwrmFybhOu1zmmYEae0nuImSnmTXU6WZRzipyIb1jV/R5j6+m6z/FRoe3F8uis1HXTjTdlu5NX9VD74klZpw0E+U3F3ZcXVfdi+RbF14jCIA/jo6Dhz4/Gka0kiB854hfNFb8IYq+5VBnWa1xC30f5hZJa5095xfFyfSs/l0i9a+VYj1tsKOUcS8zvXzn6uTR9V7k3SzgV30YVhmYR9K/YU2VY/3rrdKlVfnEslVccL9e38nOJ1LtWjvW4xYZSzrHE/P6Vo59L03eVe7OEU/FVhaHZfMpjT5VheYQD7ET5xWW48orj5fpWfi6RetfKsR632FDKOZaY379y9HNp+q5yb5ZwKr6qMDSbT3nsqTIsj3CAnSi/uAxXXnG8XN/KzyVS71o51uMWG0o5xxLz+1eOfi5N31XuzRJOxVcVhmbzKY89VYblEQ6wE+UXl+HKK46X61v5uUTqXSvHetxiQynnWGJ+/8rRz6Xpu8q9WcKp+KrC0Gw+5bGnyrA8wgF2ovziMlx5xfFyfSs/l0i9a+VYj1tsKOUcS8zvXzn6uTR9V7k3SzgVX1UYms2nPPZUGZZHOMBOlF9chiuvOF6ub+XnEql3rRzrcYsNpZxjifn9K0c/l6bvKvdmCafiqwpDs/mUx54qw/IIB9iJ8ovLcOUVx8v1rfxcIvWulWM9brGhlHMsMb9/5ejn0vRd5d4s4VR8VWFoNp/y2FNlWB7hADtRfnEZrrzieLm+lZ9LpN61cqzHLTaUco4l5vevHP1cmr6r3JslnIpvbYVBEqDHKfm228pCWagMqAyoDKgMqAyoDKgMqAx0mwx0qpaowjBFhbrbhFrTozKpMqAyoDKgMqAyoDKgMpBOBvpMYej0wQMtvAj9QHvvVO+r/OJIKq84Xq5v5ecSqXetHOtxiw2lnGOJ+f0rRz+Xpu8q92YJp+Jbe4Sh2dfrvdhTZVjvkQl7I+UXxkl8KS8hUe+o/Opxc0MpR5dIM9fKOQ1X5ZiGY2wsyj2WWJz/VHxVYYjjXtt3qgyrnYB+HlD5xWWg8orj5fpWfi6RetfKsR632FDKOZaY379y9HNp+q5yb5ZwKr6qMDSbT3nsqTIsj3CAnSi/uAxXXnG8XN/KzyVS71o51uMWG0o5xxLz+1eOfi5N31XuzRJOxVcVhmbzKY89VYblEQ6wE+UXl+HKK46X61v5uUTqXSvHetxiQynnWGJ+/8rRz6Xpu8q9WcKp+KrC0Gw+5bGnyrA8wgF2ovziMlx5xfFyfSs/l0i9a+VYj1tsKOUcS8zvXzn6uTR9V7k3SzgVX1UYms2nPPZUGZZHOMBOlF9chiuvOF6ub+XnEql3rRzrcYsNpZxjifn9K0c/l6bvKvdmCafi2+8Uhl/96leGX5suxfNSZVib791Nz1J+cbmhvOJ4ub6Vn0uk3rVyrMctNpRyjiXm968c/VyavqvcmyWcim+/Uogr7qYAACAASURBVBiWLl1qtttuO3P88cc3S9eKfcKECWbrrbc2jz/+uHU3/jRVhsU/uTdCKL+4fFRecbxc38rPJVLvWjnW4xYbSjnHEvP7V45+Lk3fVe7NEk7Ft98oDDTYt9hiC7PNNtuYjz76qFm6Vuy//OUvM4Vhk002Mc8995z1T9xpqgyLe2rv+FZ+cXmpvOJ4ub6Vn0uk3rVyrMctNpRyjiXm968c/VyavqvcmyWcim+/UBhWr15tttxyS7PpppuaRx55pFmyntjvvfdeg8Kw+eabm48//tjjo/pWqgyrflJv+lB+cfmqvOJ4ub6Vn0uk3rVyrMctNpRyjiXm968c/VyavqvcmyWcim/XKww00HfdddeswX7hhRc2S7UkdsygUBqGDBlifv7zn5f49P+VKsP8sff+XeUXl8fKK46X61v5uUTqXSvHetxiQynnWGJ+/8rRz6Xpu8q9WcKp+Ha9wjB+/Pisof6v//qvBvOgvnK2adIll1wSPfE6VYb11fv39XOVX1wOKK84Xq5v5ecSqXetHOtxiw2lnGOJ+f0rRz+Xpu8q92YJp+Lb1QoDPflbbbVVpjBcc801zRINiF1Mk/7xH/8x2jQpVYYFJLMnvSi/uGxVXnG8XN/KzyVS71o51uMWG0o5xxLz+1eOfi5N31XuzRJOxbdrFQZ69M8///xMWfj2t7/dp6MLkpWk6d/+7d+yNE2ePDlqlCFVhklaBtpR+cXluPKK4+X67mV+p510qnn0kUfdV27kupc5NgKsZqTKuSY4J5hydIC0dKncmwWdim/XKgzMXfjmN7+ZNc6vuuqqZmlGxH7cccdlafqP//iPqFGGmAxjkvftt94ekapyr2vXrs0aCNPvusfMmzvPPL/i+fIAif997733zBOLnjAzps8wy5cuM+++8070E/qKX2xevP7a6+ahefPNA/c/kHGus6LXO++8Y5YtXWbun3F/ll+vrHklSjkFbl/wQqFe8/Ia88iCRwyy9uiCR8xPXvhJVF5//NHHZtXKVebBOQ9m8rLo8UVm7U/XRsXh8/zmm2+a26bd5vvLe68v+JEQRlWffeaZ7N0fX/i4QZ5SOdg+OHuO2XG7HcyTi5+sjBbuixYuMvfeM92QFhjGur7i6EtnrAzw7T737HMZs5n3zTQrn19pfvGLX/ii7vN7fcmZ7/6p5U+Z2bNmZ2XWC6teiCqvUn3zvSavIUKVor6JreNC0hXrp6/k98MPPzRLnlyS1dfT756e1VmpytzY8qbJtMTwLcu7rlUYrrjiiqxh/nd/93e1JhmXvXQn/82cOTNLFxOgyYRQF5ph98+4zwz+3q5mp+13DI260B+KwtgzzsriGrrrEHPkoUeaoYOHZg2GQ0ccah579DFv2OOPGWWG7LJ71I9GoutuuekWM/LA72fPH7bnvoY00Fjhd8Qhh5sXV7/oBim87gt+oXnBh37h+ReaYXvtm73bQfsdaHbeYVB2Tj7+6MyxBqWpyq1Zs8ZcNP6iPOw+e+xthu+zXxYP7O68/c6qKPL/2+RFg+HuO+82ew3ZM0vrfnsNM7sNGpydk9fI2gP3z8zT5jt59913zY8vm5y/+0HDDzSD/t/3coann3xapmz6wlbdY+PFk0afmMUVugljm/xI/8LHFpqjDz86e//dd97NDN97Xb7Db5+he2dKWNl7nnLiyaXf6+BBu+b5QZzPPvNsYXRLlyw1xIc/fuTFHrutKzcO3O8As/TJJYVh3T/a5ug+X65jZIBv9YbrbzB77r5OnpHlEQccbL73fzuZXXbcOftG63QESFqaOPYV5xnT7zV878gJnI489IiMEzJ8wbjxpZ1qqb75/iKv77/3fuk3WlTnuh1sqeob5DC0jmtCZu0425ZfOqQvu2Si2X2X3TPZpe7ec7c9snPq7B+cOqZWx6a8U0x503RaSFMoX0l/0bErFQYAbrvttlnDfI899ojqrSh60VT3EYQvfvGLWdr23HPPYFOpqgz78IMPzbhzzssElsK3U4Xh7bffNt8/cEQW3/XXXJf3jJH+G66bkt2nUbv4icUboTl85OF5OqTRUHWkV10cvXDjzzs/i+PiCy42Kz7Zv4Jnc37sUcdm/9GImTNrjgQrPbbJLyYv3njjjUwRo3F76y23mldffTV7DxoU9M5KIURDa+XzxSM7q3+yOleoaKzZyhSKxKijj8uYjf3hWEP6qlxbvPhWRx0zKkvbhAsnGN4Dx30awSiGIjvXX3u9N9mvvvJq1iimQYZSBFMcFeWsmbNyLjBmpCrWXXv1tXkakMEQ1xY/0oKyxbufcNzoTGZkFTby/eILLsrSTnlw+aQfF5aFkgfCuuoo+eSyoJeYZxH+yh9fYX72xs9yL7CncUhap94wNb9fdtImx7J0hMrABx98kDV6eX+UtuXLludlJw0+yY/Dvn9oNppW9sw2/+sLzlKPUM/Q+cQ3j6PumfJJHUPnk6+zJNU335/kFQ5V36Xvf0acxaWqb2LqOHl2k8c25Rc5Pe3k07K8oMxlNJ96gXYLdRZ1NflAB877779f67VDy5s20sILVPENfcmuVBiefvrpvBf/jDPOCH2X1vzttttuWfr+4i/+wqDth7iyDMMEgx4stF0ZAUBg6zoKg2OOODoT+jNO+4E3GuyYeQbP5Pm2Q2GgoUsDcMFDD2dDdvQq0pPDx4VZERUpw9DTbroli8eO45orr87unfujc+1o8/OfvvrTrKeF59NjJw3E3IPnpC1+sXkhDXnMNnzu4fkPZyx4VxoZPkdP2/77Ds/8HT7yMK9CQONZennhW+Xa4iUNqLPP/JG3MYtJkoy2wMAd1aKQPnTEIdm703D2OdgSlh+9cCHyIvFQAUgDmPApFIZYGZG0+I40wkkXI3E0SF0Hn6MOOzJ/f1sxt/3SMMMfpnCYc+GPUT9MiegU4Lt95ulnzIrnVpgXVq3yjtryv+TVhedfYEefn1MGSF6EKG9tyWGeQM9JjAxQXvJ+9JbbSrsd7UnHrxut4tvvFtc2Z5RJODHy+fJLL3sxnHriKZkf6gjbpfrm+5u8isIwYv+DzbVXXZN1DmAOU1S3njP27KwjBV7iUtQ3KcsvSVenx7bkl86YMaeOyeSSEXGfmeXTTz2dj26jUNj8Q94ztLxpIy2S3jK+4ifk2IrCQCUNHLS1xYsXm3vuucdMnz7dW2mRaF4Okx9++KvryGh6etHQH3roIbNw4ULzyiuvRAuA+/zrrrsuTx9xhriiDMPWe5eddjZUVNjOHX34UZkwdzLCcP654/LCvKhxRUOOnkIK/QP23X+DhhQKA/Z8IY5081GJI6+l8Xvz1Jvl9kZHGpjS8AixLW+DX2xeMLdA3oHKq8hhUiL+fKMMD81/KP9/2s3TiqIxEydcmvnDTKnKtcGLbwvZ5d2uvqJ4nhGmCfL+KKG2w5Ze/mNExueQKcy8xB8N4hCHSZ6M8EjYThWGWBmpSqeMKtLjVeTuuevu/N3PON3fgTL6uNGG0adO3GWXXJY/R0YFffGJgjfunHG+vze414YcbvBA5yJGBqgnRLmkUVHk5s+dn3Oip7wbXJucmcsh39PcOXMLX/+AYftn/k4cdcIGflJ98/1NXkVheOmllzbg4btgXsfeQ/Yy9qhsivomdfnlS3ude23J741TbsxlF3OsIoe5ksh4Ub3kCxtT3jSdFjt9RXxtPyHnjSsMDLmw6dn//M//ZDs1iyLA8S//8i/N1KlTNzLrOeGEE/IG+bJly0LeYwM/2FTfeeedZvjw4eZf/uVfzKc+9ak8Pp674447mtmzZ28QJuZi0aJFeXys5BTifBlGg+u4o4418x5cX+gyDCaCGhKv64dKT3oJqyr0M8f8MH+WPQly0qWXeXs73WdR8FHBMgohDiWNXmDegcZQkaOHRd6Tj7PKNc2vTl4wIVLe4b57iwsf21zE10N8+eTL83h8JmLCRkZzUPSqGr5N8yJNmB6I0kmva5G7bdqt+fsxj8Z29ggMJg5FTnp+4Y1pTpVDDo898hjDSNrB+x+UP7+Km8TbBj+eNeaU07O0Yf5S5BYvWpynn/kgPofS3qnCICMZfNNlvWqSZrujwJcm7rXF0ff8WBlg5Sj5njEpKHL2dx8zn6MovhT32+TMvDg40aAtkxPMUpEltxMk1Tff3+QVawQ6fUIcc0MYff/Zz9abBNpyV6e+qVPHhaQ1hZ825JeyX5RY5BeTzyLHdy1lAaPoIS6mvGk6LW56fXxdPyHXjSkMAFmxYkXWOKeRvsUWW5gJEyaYZ555xjz//PNZr/+XvvQl82d/9mcbfBQ09mnQi2Lx1ltvhbxH7oeRjBEjRuThBw0alD0X21Q2XNt8883z/xYsWJCH4wTlhlWQVq5cucF99+K1117L49h+++03Unhc/1yHZhiNHBFUXzxV9+647Y48fNVEUwpyeVboR2E/n95i7P3cRtjCRx8zl158iaFHpMjdfONN+bNRHqpcW/zsdFTlBe99x223mysvv9LQI1TkxPwL1jT+XHf2WWfnLIpMmwgjvR6YLVW5tnihDDBMbiucbtrsERS3YUwlxoRxzBdef714RSC7N7FMEZVno/Qyke2tN9/KzP1Ezl1ZFf/usS1+rFBy6cWXeuVC0sTolaQf5dPnUNg6VRiYZC/PKTIz4dkoLfiDcZVri6MvHbEyQO+rvP+5Y8/xRZndQ+nHH41h2768MEALf7TFGTMt3js0/zG3dF2qb77X5FU4UUZhooilgO1S1Td2nFV1nO23yfM25BflS75v5k+W1QU0/qXjEwuBMr/CJaa8aTotkiY5hvIV/0XHRhQG4KIcSKOfzddojItbunSp+e///u/8f8yFxJFR//RP/5T99zu/8ztBjXEJyzOYVyDPHT169EYZzcjDpz/96czPNttss4FZFAoFYZcvXy5Reo8oNb/1W7+V+f2Hf/iH0l4WiSA0w2RCMIJdx9mTTFFsyhx2zfIBYbpR1lvkxkPv8q477WJoMNZxNG7k2Z3aQtvP75RfE3HJ6kFM2vXNeWE1KWFB49tXOHFPVmFiNaEq15a8VaWD/+2GWJmpR1lctlLl9li64Vg6mJEP7O1xzA8Svj62bniuu4nfXXfclae/aP7KyaNP6lhhYCUq4XTV5Vf6sGT26uLHtzKaG6ivONaRAVaOkndjvlCROSfKBP5YEahbXFucMTMVRsyPadJVffO9JK82R+Z5wbju8udV9Y39rJT1pR1v7Hkb8muP0LB4Q5VjBFVkHTO8Mhdb3jSZFl86Q/n6wtr3GlEYmHcgjfavfOUrGzTKefjFF1+c/48/26aPEYK//uu/zv7/+te/7m082S8g54RDAZDnYgblczQYWN1I/F199boJpDSWv/vd72YjIfSAlDniwJyKOL72ta9t9H6+sKEZJpOaENTQxo08j95UEXCWtatyNF7FpIRw9HaGOhq5hPVNGqqKw570jFlUyHu2wc9Ndyd5IXHZw+9FczWY0GwvOTt54uSNlDfWgSePWOIypFezL3jJO7tHsdMn/Qse3nBUz/Xru0Y+ZF4MSw4XNeQIS884ss+ymOKaVBhSyIik03ekXJJ5TeT92wUjrqeccHLHCgMTMKX8oBfZN1dk/Hnr5qOcdcaZXfvddiIDzAURBjQYKFNtx9wvJkQzd4f9MrrFtfW926aBZd9hp1xCvvlekVeXFYoQ84TquJD6xo636fLLflbZeRvySztTvm2+3yqHSbf4Zz+gIlenvGkqLUVpDOVbFF7uJ1cY6OX/kz/5k7xBft99G9t2T5w4Mf+fzdnsBjrnX/7yl7P/2RwtpDGJH3si8tZbb106MjF+/Pj8+TvssEPW4H/ssceyeygaVc/k/3//93/P/P/pn/7pBqMnAtY9hmYYq52IkFalw30GK6BIWCYyhzhZa5xwZWYldlyki0moNFJiHYqdzNOgsLLzviyuNvi5z+8kL4iLJVZFEagaFaC3VvYdIC9YqlBGXmi4oZwxBP/Si9UT5nh2X/By+XFNgwtbXN4JnnUcbESuy2zLkSVWokIu7W+nSYWhUxmp4jF54qTs3RkWLzMTwqSLUTvKX+SGUQkmol8x+fJs5SRWRgoZQcS8TlhzxJyOzfhwMumeeVe+kTLfu7Qth53KAGaU++6xT86AVetgiTyxFC12+8wRK5vs6+PQ9L22OEt9wVEcq+6wpDSLG1xy0YRsbX/kzf4GxW/oMfSb7+/y6vKgMw1lfeqUG92/Kq9j6huJrOnyS55TdWxDfml7iDkdZVuZyTTptedZFi3/Xre8aSItZYxD+ZbFwX9JFQZMdez5AzTGuec6YN1yyy1m2rRpGzUYqfB+93d/N2uM77777kGFDnMONttss1wJuPHG8o9t/vz5ud8tt9wySwOKAvMsqsx45F0GDx6cxfGFL3whqcIw+tjj88oqtsBl8rFU9qHD5fbqM0UfhbyzHNG2eU7ZKgPil/xnTgPzKVgmjqF+nsl6+z7ZkHDuMVTgO+HnPjM2LhrHvCvzGmjA0YtB44pl1kIcw5r0XkoecqTxiyJx3tnnBo0syHP6gpc82z6K+QYNXirDWEfDVCYtkx+UHUWO3m8ae+6oV5MKQ6yMFKVd7rNUMTbyNM6ZH0SDFWXT3bxJ/MuRhj1mhcgPlSJzReg0sGWJxm6IGQmNMLtiReFjjgTHm268qau/2xQywCiCXS7CEFNPZJjODv7vNtfG987olsgT86iYn0AvLLKCGQwLGthyw34yvjkMVexivnni6s/y6rK4+pMlyUNW3+q0vuHZqcsv931Cr9uQX9JCfSwyTEdAkaPtJSO7+C/aNLWT8iZ1WorehfuhfMvi4L+kCgPDLGLqw1HMfaoSYf+PwvDrv/7rWTxHHrnhiiq2Pzmn0XnIIYfkz/3zP//z0kYF4VjeVdL5+7//+9k1k6EZeQhtxJI24iCtIb1toRlm283FKgwItXwMISuYwALFQsKEzkc48wdnZo0H34Y8ki9yZE8IiZ8jPVPYYbNUWdlEYQkvxzb4ybPkGJsX9iorvCvDytjbx5gu0OtBxWszo0ezrHCT9NrHvuBlP59zNlXiPZjrEsPAjoedr4kDpaHMFItVQ2is+EbJmlQYYmXEfjffuaz8wjujKPL+rKJWZIokcTy64BHDPAaW/LQbaWy8RnjZjRdG15WsACTxMUohvckii+zWHrsqUJtymFIGWCQDJV3eXY5M0A8p74VjW8c2ODO6KRwoo5BVRkKZlC91FY1Yu9efPQdiTZdCv3mbbX+UVzv9nDMCiFLKROQQl6K+SV1+haTb56cN+eW5dqfqyINGFn7LmA+LrHNkI0LXdVrepEyLmzb3OpSvG869TqowsL+CNMT/4A/+IKjn3U0QCsPnPve5LJ799682q6HH0Z67cMABG6/a43uGpJPj/fffb771rW9FpZe0EZa0kuYqF5phsikQQiqFcFXc8r89eTa00LF7GFhhqcpR+NOQoQc9xPEO9PiyKgBD+rMfmJVvKkcPapltoB1/G/zs53EemxcoQEwGZz8NJkmxDCo9veQlcZWZlMizGYpnFAa+tjJHHGxXX9VwlHj6gpc8myMKIYoO8wlC89gOzzk7lPPeNPjLRicwiUApsdcst+NqUmGIlRE7Xb5z5Ic8RlYYDZD9SuDI9xlbJsgzmJ9k70nBSEaR4xnMAeE7Z6UW20QHhYNe0BDzJuJvSw5TywCrAbEaFA1eOkiQQ/mxgo1vP5Uinm3cb4Pz08ufyhnAAnMWFCufY78Q4RW6lCjxhH7z9jP7o7za6ZdzRphhxiheiEtR36Quv0LS7fPThvzyXGSFUXuRTawe7NFbOgOQV0Z27e/ezZMU5U2qtPh4uvdC+brh3OtkCgMViD03gP0PQisVO1E0vjfddNOsMb7ddttVVpCMFtjmSCGjGjzDVhg4v+KKK4JHF0gvaSMcaU2qMIxet4soAo1AxTjMfuRDCJ00ZfdosvxllZt6w9TsGaGjEb74KOjs5xbt8GuHDRX4kzrgZz+P8xRxYf8t8xgwUWKeSZGT3VMvvuDiTBbJf+YviNJB3tKICdmuvi94yXvRC83kZHrLyt5X/PuOovyi+NKILnLsjszoAxVf0ffSqMKQUN6K3pFlV+W7/tGZ9TdnE6bEhWmNz9EARO5RDNg1GkclynePUibpYCg+xLUhh6llgIYbyy4ec8TRuewxQmZPEEWBe2HVCyEIWvHTBmd7jhydAWXzqZjDILJCuWfvJ1AEROSz6pu3w/dHebXTb59jwgWzuqsjEVdMfYP/FHWc/Q51z9uQX0kb7VKWkRf5ZH4gHaf8kGsW2KBzk32BxI/dRklZ3nSaFnmnqmMo36p4kikMmPLss88+eUOc86IKvCxRjBh89atfzeJhydKqOJ599tn8mTTg2eOhyrkKw7e//e1KMyY7TtJE2ngeaS3qZbHDhGaYrf2GmkfJc+iRFAGnVzDE2ba6LCtY5nhvRgV4RsyKSr44mSsiPZ40RPhAy1wb/Nznd5IXdlzMa5B8YQjfN9FbNmWj0HLznf0JDjl4ZB4HS5NWfRd9wYt3ZklA8pMVfaomldmM7HNMZuDFe1aZf7BaDwoZduWY4vAtwpdyRDg2qTCkkhH7/d1zKhW7sRoyd8iNg2tGLUQOObrztZApWQXHt3wry9Tu/smmjIQv28ldnt+GHKaUAUZeaDRgjuWa0pAPKGzCkIaFO19G3rvtYxucaYzKu6NMVTlM2MQ/SliZi/nmJZ7+Kq+SfvvIhGWUdDpZqsp2O5zvPKS+kXBtlF/yrLJjG/LrPp/2zl133GnYYJARfcyOKOOkA5gNQkV+7fmdKcsbSVPdtEj4qmMo36p4kikMFKYsS0ojmt/YsfV6wojnb/7mb7I42NhNKv2iF5k7d27+TJ4b0ninMSHp5Pjoo48WRe+9T5pIG2FJK/FVudAMk14GBLXq3d1nMkwuAs4RTbjKSaMd/24F6YaVHWfdzbdcf6HXduVbtQlUG/zcdHeSF25cTDiVvGHpO9thskNlwf9FPfI0TGQjGfzNmD7DjmKj877ghf08ZiyMbmE7X8cxwZf3Yw12yoIyxxKtwjT0iLwzQbNqxam+4Ff2rvRwyTuGbNxXFBc94xKPayrGClT8h0lcEXtMzSQ8MltlZtc0x5QyQPkn3ymr/vgcjTl753bmOXSDa5oz72hPekYZqHL2Pgls1FnkYr55O47+KK92+u1zWZEnZFNKO1zRucgx36pb39hhUtZxdryx523Ib2yafjhmvVkdI2a4lOVNTHp8aYkJH8q3Ks5kCgMF6cknn5w3xG+91V/g2gmiUmJ1IrvBTSOZ+QQ0xn/t136tsOKSeFAQPv/5z+fPteMSP+6RXktbYQgx8bDjIN2kjThIa8gzQzPM3hW4qNK202Kf816YgkiFXjVigO2eNFTpha1y0sC/8PwLqrwG/Y9GL2mlZ73MtcHPfX4neeHGdeKoE/J3dSedjj3jrOy/qqVwxRwMZuxfUeba5sXIAsoCK0uUTU4uS7OsEMKkx5BeNubDsCpQ1U9kHG6237JRsrb5lXHhP3q+5FuBc0iZ44vTXjmJSXvi4C0bPrk7zIofOdqNDFYEK3NNc0wpA/Q2CmN6fIucbcuPAhoiq0VxpbrfNGdJp3Qw8U1VySByJDxZ3tfnYr95iaO/yquk3z5Sz8tozKyZs+y/ap+X1Td2pCnrODve2PO25DcmXVgDIL+YJ0rnbcryptO0xIQP5VsVZzKFgQddf/31eUOc+Qxljgw4/fTTM/+jRq1fn52C4NBDD83jKavUiZ/hI3tkg12kyxyFHM+zFYaqnZ3d+F588cU8/H777RdUYYRmmN0rE6swkE57B+Xbby2vzFlZRQr0qnWf6eGmoYJ/PpoiR/4xF4IGZJXDtEKez0TCMtcWPzsNVXnBClA0uqp2gSROexMY1ioXBy96dOFw2SUT5bb3yHOEF6YhhC1ybfJavnRZtoQsSl+s8i3pZzUp3o1N61K7Jk2SqmSk6l1Q6pGhqm+dZRYl7znapjAx3xyT+SQee7lfe3SSZV3LHOWKxDHmlNPLvAZPeu6UY2kiAnf7lom6zBOqcjZHdm3ta9fW926bsFStfCYdTMiKb0GNTr75XpJXe7M1FswocinqGzfupr8793lF123Jb9Hz3ft0pkp7p6oz0w3LdZ06xxcP9zpNC3GE8i1Kg9xPqjDQw73VVltljekhQ4bktmDyMDmiBAwaNCjzxwiD2IzJ/+zPIA36hx6qnog7a9Ys8+lPfzoLc84550g0GxxRUK666irzz//8z4YVnGTSMs+5++67M79MmD777LPNK6+8skFY94I0SfrYtTrEhWaYXSBXNSJ8z5WVFiik2fa9zMkkW3qL1v50bZnXbGMeaSSUFWqYyoi/qoa0DMPiP1XDo1N+NoSquGTZxZBeRtuMgfkK4pjQK7yqzLLouRe/KAxlri15QwFnJRnmLISYwD391NMbrWPPcszYjdPACHEopDGjGMP2Wm9LXaZk2c9ug59tE17VUy/mgOQ/oyS2C/3m1q5dm8sP8di96JiTiWyVmTDwXBQK8cuchzLXBsey58t/ITKAKR3vFWLyhVIhDDqZpCrp6/TYFmd7Iijld5mzGblKVafffC/JK6vfIUuUoWUuRX3jxl9Vx7n+m7puS35pfDP/hkVXyibiT7/rnvz7rjPqE1LetJUW8iyUb1X+JlUYqIyXLFmSN6aZ+Dx79uxs6PKRRx4xEyZMMCxH+o1vfMOw7wHXvgmgTMb7jd/4jSyekFWPaLSgeNCIZ+flo446yixYsMCwoRtKAiMW//mf/5n9jyLz5JNPGkYJZIM4RhyYx/BXf/VXhs3mqhrq559/fv6Oc+bMqWKc/R+aYbZmWpQOVqdgXX57vXVJBEqbPTFxzZp1u7TK/3JkpaIDhq3b2IlhySrHDrohhZpMXsMvS4+VOXsNaMwBylxb/Ow0VOUFPQ+8JyMEZY6GtMw/QDkTRiWtSQAADKtJREFUe0gJw+gK8WDqUeZY7xx//JgEW+ba4oVtMunBTCPEwdRecYIwyB/mMCEKB2UDq1pUjTxKWiiTpKeIdKZWGKpkhHQUfa/sGSH5WaUwsGSs+MWEzXb2f2XfHKseSRzuN0/lKf9h+lbmJM/x75scbYdtSw7tZ7rnoTLAyB/vRPlJ+VjkiE9MP9nMLmTeXFFcqe63xVmW1YYTClbR98RII6sj4Q8zRdd1+s33iryyb4WUT+zxUeZS1Tf2M0LKL9t/U+dtyS+dTcgkv7KJ+HS24ocV+MQcKfTdQ8ubNtIiaQ7lK/6LjkkVBh4CXFYq+s53vpM3qqU3niPKwpgxY8yqVasKM4JG73/9139l4W1zpaKX4D6Kx7nnnut9Js/ddtttN1BQaIxjNvXZz342D4Oy4FNg3OfKvg9/9Ed/VDlRWMKGZBgTCKXwQFh9hTG9hDSY+L+oJ4zJjBIPDUuf4iG9+6ySZJs3SHrtI71BPI9fVaHG+uX4o2HsTqq040RrlzixFayyh22TH+kMyQvphWACLflS5JjzIe96wbiNTfUmT5yU/1+2MZZt1sR+DWWuLV4onUywIz8pgEkXE2N5D1acQZHgx/+seQ0Hlv8Vx8gW90YfNzrb6I2eQ0zl2Ihp2dJlefjly5Znu23KakFVI2ISP6Zxwh6ZDHVt8OP7FtvlsoY3358saYo9Ld+Y7TBX4t34Fe0PwCR0maPASiw+xeLwkYdnrBgxKxrBoXyURgZpqZrc3gZHm4XvPFQG7I2Uivb0IH7blJIJu93g2uQsG6vxXdmjpTYH8UNPKyNptkv1zfeCvE6/e3pePpWZ+sIvVX0jeRFSx4nfpo9tyS/7+SC3lJXuKnHyjjdOuTH3M3fOXLkdfAwtb9pIiyQ6hK/4LTsmVxjkYTQAGVW47LLLsrkKN9xwg1m8eHFmflSlsVGRjh49OmvIb7/99oWKhTxLjsRL5coGcmeccUZmXnTttdca5jW4Zk+E4d7LL79sLrzwQjNjxgyvH4lbjm+99Zb51Kc+laXtiCOOCE5bWYahIFFZ2Rt10TvDMqk0oOzh3HvvWV/AIPhFG3mxDBgfBX5opNoNLJ6FCQgNCLcwl/e0j3bjnn0Bqhw9kCg1NCgwX2DoTRzKCcuViULDMGtIz3Jb/GLyAr9ipys7nkrPJDLMUrEy3Mz7sgylTwmkEcaybuQVS4TSoLa/EXbUxraf/1nlZsb0ewVn4bENXvSsioyRttAfMi3O3a00NI6yvRmIm7xhp2l72WBGebC7XvLkksLvRtLVBj+ehXIlc1ioqGzFE5tlGhSy4zK761K++Rw24sgYu5cufmJx7oXOAsqR4z7pMaOTwS4Lco/GGBQPmRSNPLujRpQVIqdwRYmrcm1x9KWjjgzICCnlFzsWE4c4vkkabXyDyD2Kvnzv4qevjm1y5p3FxBKZY+RKOKFAyhLRPmUBPqm++V6QV+lMQp5YPrvMwThFfUM8oe2NsvSk/K9N+WUPGeoZWNpWGoz2sDoaecHIob2Uasi7wjW2zmkqLW56y/i6fsuuG1MY3If6GkquH/t62bJl5rd/+7ezuQkoGt3iJk+enCkLn/nMZ8yiRYuCk1WWYZgjlDWUqJjEUdnLCEOVWQoNz+H7rNtpmPiPOOTwvAGFHV/REp7yLDlKjxo9yT7FS/zZR8wtsG2XBiVrlkvDiLTQMKnqUbHja4tfTF5I+ggjvbdUoNiiiskC1wzJr3iufGIkCja7nO6317BMFuhRJr8YEkW5gyOb7Nh25/J837ENXsxHKJPbov/sybbsoFnkr+x+VUONCeRl4Wm0lLk2+Mnz6WnCPA9ZIc3IEnIg3w4KA41X30ihxMGRBrys7EF4ygdZ1YZvl/XGq8xn6HVjFFGWX+WbRX4lPcgipjsoMyGuTY5ueurKAA1g6byhrMVkEHt85o6QP3zfKJ3d5NrmTCcGc+CknKNziB5/kVnKqqLOqJTffH+XVxk1ZsJ9qOu0vqlTx4Wmra6/tuWXRV7olEVe6WRhTgPnlMF8+yzkEevqljdNpMVNexlf12/ZdWsKQ1kifP9ROQ4bNixrnDMnoVvc0KFDszTtu+++lRW4neZUGUaczEuYed/MvFfHfo57Tq8YlRsrsVBQE84esXD9+65REqhEq8wPfGEZPaAhg906P85DRhTcuPqKn5uOsmv40GPMLtjwwjwkVMGSeFGsYUR4RmoY3UGxkx488Vd17A+8qt6hL//vC34ojSwUQN4jQ8hS1d4oPkaM6DFJmjjmPTjXoJDEOhQLRiYY0WRkjJFCTKGqlBb3OX3B0U1D3Ws2A8QkgdV8GFnARK5qZKvuszoN11ec4YHpKXyQE3plq/bm6PRdfeH7q7wyF4PvPcQU2n3vFPWNG2dfXfeF/NI2wgqAUYHbpt2Wfd+x9WwqXk2nJRXfrlUYyIiZM2dmjfNNN93UrFixIlXe1I6HCp3J2syJmDevfOdK9yGpMsyNd6BcK7+4nFZecbxc38rPJVLvWjnW4xYbSjnHEvP7V45+Lk3fVe7NEk7Ft6sVBnqzmMNAA52VjvraycZ0rP6E8hDjUmVYzDN7ya/yi8tN5RXHy/Wt/Fwi9a6VYz1usaGUcywxv3/l6OfS9F3l3izhVHy7WmEAIcuyfuELXzDMGVi4cGGzVEtiZ++Fz33uc2azzTbLlmwt8er9K1WGeSMfADeVX1wmK684Xq5v5ecSqXetHOtxiw2lnGOJ+f0rRz+Xpu8q92YJp+Lb9QoDNt2TJk3KRhkGDx68weoxzSJeHzv2ZTLSMXHixOjRBWJKlWHrUzWwzpRfXH4rrzherm/l5xKpd60c63GLDaWcY4n5/StHP5em7yr3Zgmn4tv1CgMYMU0aO3ZspjSwhwMN+LYcCgv7NWAWxZ4QsaZIks5UGSbxDbSj8ovLceUVx8v1rfxcIvWulWM9brGhlHMsMb9/5ejn0vRd5d4s4VR8+4XCAEqUBhrsNNynT5/eLF0rdpZO/fKXv5ztJF1XWSC6VBlmJW1AnSq/uOxWXnG8XN/KzyVS71o51uMWG0o5xxLz+1eOfi5N31XuzRJOxbffKAzgRGlgUzZ7c6NmMRtz7733ZgpKp6MaqTKs6fft1viVX1zOKK84Xq5v5ecSqXetHOtxiw2lnGOJ+f0rRz+Xpu8q92YJp+LbrxSGZpE2G3uqDGs2ld0bu/KLyxvlFcfL9a38XCL1rpVjPW6xoZRzLDG/f+Xo59L0XeXeLOFUfFVhaDaf8thTZVge4QA7UX5xGa684ni5vpWfS6TetXKsxy02lHKOJeb3rxz9XJq+q9ybJZyKryoMzeZTHnuqDMsjHGAnyi8uw5VXHC/Xt/JzidS7Vo71uMWGUs6xxPz+laOfS9N3lXuzhFPxVYWh2XzKY0+VYXmEA+xE+cVluPKK4+X6Vn4ukXrXyrEet9hQyjmWmN+/cvRzafqucm+WcCq+qjA0m0957KkyLI9wgJ0ov7gMV15xvFzfys8lUu9aOdbjFhtKOccS8/tXjn4uTd9V7s0STsVXFYZm8ymPPVWG5REOsBPlF5fhyiuOl+tb+blE6l0rx3rcYkMp51hifv/K0c+l6bvKvVnCqfjWVhgkAXqcku+xoCyUhcqAyoDKgMqAyoDKgMqAykC3yUCnaokqDFNUqLtNqDU9KpMqAyoDKgMqAyoDKgMqA+lkoHWFodMHangloASUgBJQAkpACSgBJaAE+g+B6BGG/vNqmlIloASUgBJQAkpACSgBJaAEOiWgCkOnBDW8ElACSkAJKAEloASUgBLoYQKqMPRw5uqrKQEloASUgBJQAkpACSiBTgmowtApQQ2vBJSAElACSkAJKAEloAR6mIAqDD2cufpqSkAJKAEloASUgBJQAkqgUwKqMHRKUMMrASWgBJSAElACSkAJKIEeJqAKQw9nrr6aElACSkAJKAEloASUgBLolIAqDJ0S1PBKQAkoASWgBJSAElACSqCHCajC0MOZq6+mBJSAElACSkAJKAEloAQ6JaAKQ6cENbwSUAJKQAkoASWgBJSAEuhhAqow9HDm6qspASWgBJSAElACSkAJKIFOCajC0ClBDa8ElIASUAJKQAkoASWgBHqYgCoMPZy5+mpKQAkoASWgBJSAElACSqBTAqowdEpQwysBJaAElIASUAJKQAkogR4moApDD2euvpoSUAJKQAkoASWgBJSAEuiUgCoMnRLU8EpACSgBJaAElIASUAJKoIcJqMLQw5mrr6YElIASUAJKQAkoASWgBDoloApDpwQ1vBJQAkpACSgBJaAElIAS6GECqjD0cObqqykBJaAElIASUAJKQAkogU4J/H8Tye/qaQb+cgAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For large $n_e$ we have $\\lambda\\approx\\sqrt{n_e}D$. In fact, for $n_e>10$ we can bypass $Q_\\mathrm{KS}$ entirely and compute the value of D that would correspond to a given significance level $\\alpha$,\n",
    "\n",
    "$$ D_\\mathrm{KS} = \\frac{C(\\alpha)}{\\sqrt{n_e}}$$\n",
    "\n",
    "where $C(\\alpha) = \\sqrt{-\\frac{1}{2}\\ln(\\alpha/2)}$. Some instructive significance levels are shown below (from [here](https://www.wikiwand.com/en/Kolmogorov%E2%80%93Smirnov_test)).\n",
    "\n",
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The KS test and variants can be accessed in the `scipy.stats` modules `kstest`, `ks_2samp`, and `ksone.` Let's look at a simple example for comparing a Gaussian sample to a Gaussian null hypothesis distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the kstest function is appropriate for this\n",
    "np.random.seed(0)\n",
    "x = np.random.normal(loc=0, scale=1, size=1000)\n",
    "print(stats.kstest(x, 'norm'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now compare a $t$-distribution with $100$ degrees of freedom to a Gaussian distribution. <font color='red'>Complete the following cell. Discuss with your breakout room colleagues whether the result makes sense.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(987654321)\n",
    "print(stats.kstest(stats.t.rvs(___, size=100), 'norm'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Try again and discuss for a $t$-process with only $3$ degrees of freedom.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many other nonparametric tests, in particular the **Anderson-Darling test** to check whether a sample is consistent with having been drawn from a Gaussian distribution. We don't have time to do a deep dive here, but see the textbook for further details.\n",
    "\n",
    "Read through and execute the following cell to perform the **KS test**, **Anderson-Darling test**, and the **Shapiro-Wilk test** applied to $10,000$ values drawn from a Gaussian distribution, and separately to a mixture of a Gaussian distribution. `scipy.stats` has many of the relevant functions. The print-out statements give statistic values and $p$-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load code/fig_anderson_darling.py\n",
    "\"\"\"\n",
    "Gaussianity Tests\n",
    "-----------------\n",
    "Figure 4.7.\n",
    "\n",
    "The results of the Anderson-Darling test, the Kolmogorov-Smirnov test, and the\n",
    "Shapiro-Wilk test when applied to a sample of 10,000 values drawn from a normal\n",
    "distribution (upper panel) and from a combination of two Gaussian distributions\n",
    "(lower panel).\n",
    "\n",
    "The functions are available in the ``scipy`` package:\n",
    "\n",
    "- The Anderson-Darling test (``scipy.stats.anderson``)\n",
    "- The Kolmogorov-Smirnov test (``scipy.stats.kstest``)\n",
    "- The Shapiro-Wilk test (``scipy.stats.shapiro``)\n",
    "\"\"\"\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "from astroML.stats import mean_sigma, median_sigmaG\n",
    "\n",
    "# create some distributions\n",
    "np.random.seed(1)\n",
    "normal_vals = stats.norm(loc=0, scale=1).rvs(10000)\n",
    "dual_vals = stats.norm(0, 1).rvs(10000)\n",
    "dual_vals[:4000] = stats.norm(loc=3, scale=2).rvs(4000)\n",
    "\n",
    "x = np.linspace(-4, 10, 1000)\n",
    "normal_pdf = stats.norm(0, 1).pdf(x)\n",
    "dual_pdf = 0.6 * stats.norm(0, 1).pdf(x) + 0.4 * stats.norm(3, 2).pdf(x)\n",
    "\n",
    "vals = [normal_vals, dual_vals]\n",
    "pdf = [normal_pdf, dual_pdf]\n",
    "xlims = [(-4, 4), (-4, 10)]\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the statistics and plot the results\n",
    "fig = plt.figure(figsize=(5, 7))\n",
    "fig.subplots_adjust(left=0.13, right=0.95,\n",
    "                    bottom=0.06, top=0.95,\n",
    "                    hspace=0.1)\n",
    "\n",
    "for i in range(2):\n",
    "    ax = fig.add_subplot(2, 1, 1 + i)  # 2 x 1 subplot\n",
    "\n",
    "    # compute some statistics\n",
    "    A2, sig, crit = stats.anderson(vals[i])\n",
    "    D, pD = stats.kstest(vals[i], \"norm\")\n",
    "    W, pW = stats.shapiro(vals[i])\n",
    "\n",
    "    mu, sigma = mean_sigma(vals[i], ddof=1)\n",
    "    median, sigmaG = median_sigmaG(vals[i])\n",
    "\n",
    "    N = len(vals[i])\n",
    "    Z1 = 1.3 * abs(mu - median) / sigma * np.sqrt(N)\n",
    "    Z2 = 1.1 * abs(sigma / sigmaG - 1) * np.sqrt(N)\n",
    "\n",
    "    print(70 * '_')\n",
    "    print(\"  Kolmogorov-Smirnov test: D = %.2g  p = %.2g\" % (D, pD))\n",
    "    print(\"  Anderson-Darling test: A^2 = %.2g\" % A2)\n",
    "    print(\"    significance  | critical value \")\n",
    "    print(\"    --------------|----------------\")\n",
    "    for j in range(len(sig)):\n",
    "        print(\"    {0:.2f}          | {1:.1f}%\".format(sig[j], crit[j]))\n",
    "    print(\"  Shapiro-Wilk test: W = %.2g p = %.2g\" % (W, pW))\n",
    "    print(\"  Z_1 = %.1f\" % Z1)\n",
    "    print(\"  Z_2 = %.1f\" % Z2)\n",
    "\n",
    "    # plot a histogram\n",
    "    ax.hist(vals[i], bins=50, density=True, \n",
    "            histtype='stepfilled', alpha=0.5)\n",
    "    ax.plot(x, pdf[i], '-k')\n",
    "    ax.set_xlim(xlims[i])\n",
    "\n",
    "    # print information on the plot\n",
    "    info = \"Anderson-Darling: $A^2 = %.2f$\\n\" % A2\n",
    "    info += \"Kolmogorov-Smirnov: $D = %.2g$\\n\" % D\n",
    "    info += \"Shapiro-Wilk: $W = %.2g$\\n\" % W\n",
    "    info += \"$Z_1 = %.1f$\\n$Z_2 = %.1f$\" % (Z1, Z2)\n",
    "    ax.text(0.97, 0.97, info,\n",
    "            ha='right', va='top', transform=ax.transAxes)\n",
    "\n",
    "    if i == 0:\n",
    "        ax.set_ylim(0, 0.55)\n",
    "    else:\n",
    "        ax.set_ylim(0, 0.35)\n",
    "        ax.set_xlabel('$x$')\n",
    "\n",
    "    ax.set_ylabel('$p(x)$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonparametric Modeling & Histograms <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "Imagine you have some one-dimensional (\"univariate\") data that you would like to try to understand.  Where by \"understand\" we mean \"know the distribution in the measured space\", i.e., you want to know the probability distribution function. Our constant friend is the histogram, and it's usually the first thing any of us do on new data. Simple, right? Not quite...\n",
    "\n",
    "Let's work through some examples to see what problems we encounter and how we might overcome them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to generate a univariate data array, x\n",
    "# this is the same data used in Ivezic, Figure 6.5\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Generate our data: a mix of several Cauchy distributions\n",
    "#  this is the same data used in the Bayesian Blocks figure\n",
    "np.random.seed(0)\n",
    "N = 1000\n",
    "mu_gamma_f = [(5, 1.0, 0.1),\n",
    "              (7, 0.5, 0.5),\n",
    "              (9, 0.1, 0.1),\n",
    "              (12, 0.5, 0.2),\n",
    "              (14, 1.0, 0.1)]\n",
    "true_pdf = lambda x: sum([f * stats.cauchy(mu, gamma).pdf(x)\n",
    "                          for (mu, gamma, f) in mu_gamma_f])\n",
    "x = np.concatenate([stats.cauchy(mu, gamma).rvs(int(f * N))\n",
    "                    for (mu, gamma, f) in mu_gamma_f])\n",
    "np.random.shuffle(x)\n",
    "x = x[x > -10]\n",
    "x = x[x < 30]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next few cells, make a normalized histogram of this data with 10 bins, 20 bins, and 100 bins. It starts off looking unimodal and Gaussian-ish, but clearly when more finely binned the data breaks up into several modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10 bins\n",
    "plt.hist(x,___,___); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 20 bins\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 100 bins\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How would you characterize this distribution?  Could we reasonably think of it as a normal (Gaussian) distribution that we could characterize by some mean and standard deviation?  Maybe, but even just by looking at this plot we see that it wouldn't be a particularly good description of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find that small changes in parameters to the histogram function *significantly* change the PDF.  That's bad, because the underlying data clearly have **not** changed. One of the problems with histograms is that some bins end up with little (or no) data.  We can fix this by making **variable-width bin sizes** that have the ***same number of objects in each bin***.  How can we do this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Execute this cell\n",
    "a = np.linspace(1,42,num=42)\n",
    "print(a)\n",
    "print(a[::2])\n",
    "print(a[::3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are familiar with slicing, then you know that `[::2]` and `[::3]` say to count by 2 and count by 3.  But that isn't what they really do.  They say to take every other index of the array or every 3rd index of the array.  So, if your array is sorted (like `a` is), then you could use this to instead define the number of values in a bin.  That is for any given value of `M`\n",
    "\n",
    "    bins = np.append(np.sort(x)[::M], np.max(x))\n",
    "    \n",
    "would give bins with `M` objects in each bin.  \n",
    "\n",
    "So if `M=3` and $x = [1,3,5,7,9,11,13,21,29,35]$ then $bins = [1 7 13 35]$.\n",
    "\n",
    "*Note:* you need to add the maximum value to set the right edge of the last bin.  \n",
    "\n",
    "Try it for `M=100` and `M=30` (100 and 30 objects in a bin)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins100 = np.append(np.sort(x)[::____], np.max(x)) #Complete\n",
    "bins30 = np.append(____,____) #Complete\n",
    "print(len(bins100),len(bins30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Note that the underscores here are suppressing the array output\n",
    "#so that we just see the plots\n",
    "_ = plt.hist(____, bins=____, density=True, histtype=\"step\") #Complete\n",
    "_ = plt.hist(____,____,____,____) #Complete\n",
    "plt.xlim(-5,25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, this can look pretty different depending on what the number of objects you choose as the minimum for each bin and compared to the plots above.  And it looks a lot different from the plots above.\n",
    "\n",
    "So, what is the \"right\" way to set the bin size? There is no \"right\" way, but there are some useful rules of thumb."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**\"Scott's rule\"** suggests that the optimal bin width is \n",
    "\n",
    "$$\\Delta_b = \\frac{3.5\\sigma}{N^{1/3}}.$$\n",
    "\n",
    "That's great, but what if we don't know the standard deviation, $\\sigma$ (e.g., if the distribution isn't really Gaussian)?  \n",
    "\n",
    "We can then instead used the **\"Freedman-Diaconis rule\"**: \n",
    "\n",
    "$$\\Delta_b = \\frac{2(q_{75}-q_{25})}{N^{1/3}} = \\frac{2.7\\sigma_G}{N^{1/3}}.$$  \n",
    "\n",
    "Let's try that. Remember that you can compute $\\sigma_G$ using `astroML`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroML import stats as astroMLstats\n",
    "sigmaG2 = astroMLstats.sigmaG(x)\n",
    "print(sigmaG2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Now set the bin size accordingly</font>, using [np.arange](https://numpy.org/doc/stable/reference/generated/numpy.arange.html) and plot.  Make sure that you don't throw away the last object in data set (append that maximum object to the end)! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binsize = ____*_____/(N**(____)) #Complete for Freedman-Diaconis Rule\n",
    "print(binsize)\n",
    "binsG = np.append(np.arange(start=x.____, stop=x.____, step=____) , x.____) #Complete\n",
    "print(len(binsG))\n",
    "print(binsG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(____,____,____,____) #Complete   \n",
    "plt.xlim(-5,25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did you find that tedious? Me too. Try the following shortcut."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.visualization.hist import hist as fancyhist\n",
    "_ = fancyhist(x, bins=\"scott\", histtype=\"step\",density=True)\n",
    "_ = fancyhist(x, bins=\"freedman\", histtype=\"step\",density=True)\n",
    "plt.xlim(-5,25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that even those don't yield quite the same results!  But we can do better!\n",
    "\n",
    "An obvious thing to do is to simply show all of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# execute this cell\n",
    "plt.hist(x,histtype=\"step\")\n",
    "plt.plot(x, 0*x, '|', color='k', markersize=25) #Note markersize is (annoyingly) in *points*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is called a **rug plot** and now we have a better idea of where most of the data and where the gaps really are (as opposed to where the binning makes them *appear* to be).  However, the markers are all piled up, so we have lost all sense of the relative numbers of objects.  Are there ~10 at x=5 or could there be 100?\n",
    "\n",
    "This is where **[Kernel Density Estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation) (KDE)** comes in:\n",
    "- In short the idea here is to represent each data point not as a delta function, but rather as a distribution (e.g., a Gaussian).  \n",
    "- Those individual distributions (\"kernels\") are summed up to produce the PDF.  \n",
    "- One of the advantages of this is that it combines the best of \n",
    "    1. the histogram (tells us the relative height of the distribution) \n",
    "    2. the rug plot (centers the data points at the actual location of the data instead of within some arbitrary bin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just about any distribution can be used as the kernel, but the most common are a **Gaussian kernel** and an **Epanechnikov kernel**.  One downside of the Gaussian kernel is that the tails are technically infinite in extent.  So each point has some finite probability of being *everywhere*.  The Epanechnikov kernel has truncated wings.  \n",
    "\n",
    "One still has the problem of deciding the width of the kernel (e.g., for the Gaussian the *\"mean\"* is fixed at the value of the point, but how wide should you make the Gaussian?). For now, we'll just play with the widths by hand to see what might work best.  N.B. the widths of the kernel distribution are referred to as **\"bandwidth\"**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# execute this cell to load the KDE module\n",
    "from sklearn.neighbors import KernelDensity\n",
    "xgrid = np.linspace(x.min(),x.max(),1000)  # Use this instead of 'x' for plotting\n",
    "\n",
    "def kde_sklearn(data, bandwidth = 1.0, kernel=\"linear\"):\n",
    "    kde_skl = KernelDensity(bandwidth = bandwidth, \n",
    "                            kernel=kernel)\n",
    "    kde_skl.fit(data[:, np.newaxis])\n",
    "    log_pdf = kde_skl.score_samples(xgrid[:, np.newaxis]) # sklearn returns log(density)\n",
    "\n",
    "    return np.exp(log_pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we try the Gaussian and Epanechnikov kernels, let's first start with a tophat using `kernel = \"tophat\"`, which will produce a plot much like the rug plot.\n",
    "\n",
    "Start with `bandwidth=0.01`.  See what happens when you adjust this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDFtophat = kde_sklearn(____,bandwidth=____,kernel=____) #Complete\n",
    "plt.plot(xgrid,____) #Complete"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The defaults give a result that is essentially what you would get if you made a histogram with a really large number of bins.\n",
    "\n",
    "Now let's compare what happens when we adjust the bandwidth (which is just the width of the kernel function).  Try \n",
    "`bandwidth=0.1` and `bandwidth=0.5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDFtophat1 = kde_sklearn(____,____,____) #Complete\n",
    "plt.plot(____,____,label='width=____') #Complete\n",
    "\n",
    "PDFtophat5 = ____(____,____,____) #Complete\n",
    "plt.plot(____,____,____) #Complete\n",
    "    \n",
    "plt.legend(loc=\"upper left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see what we get with the Gaussian `kernel=\"gaussian\"` and Epanechnikov `kernel=\"epanechnikov\"` kernels.  <font color='red'>Play with the bandwidths until you get something that looks reasonable</font> (and roughly matches) for the two kernels.  They need not be the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDFgaussian = kde_sklearn(____,bandwidth=____,kernel=\"____\") #Complete \n",
    "PDFepanechnikov = ____(____,____,____) #Complete\n",
    "plt.plot(xgrid,PDFgaussian,label=\"____\") #Complete\n",
    "plt.plot(____,____,____) #Complete\n",
    "plt.legend(____) #Complete"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is pretty different from the histogram that we started out with, isn't it?\n",
    "\n",
    "**HISTOGRAMS TAKE-AWAY MESSAGE:** \n",
    "\n",
    "Making a histogram is the first-cut we make of data, and it's certainly one of the most sensible things to try to get a feel for the data. But we can't just do it without thinking. We need to explore bin sizes and KDE smoothing bandwidths to tease out the structure in the distributions, and overcome any finite sample effects in bins by potentially having variable bin widths.\n",
    "\n",
    "Finally, the normalized bin height of a histogram can simply be understood as\n",
    "\n",
    "$$ f_k = \\frac{n_k}{\\Delta_b N}$$\n",
    "\n",
    "where $k$ indexes the bin, $n_k$ is the occupancy number of the bin, $\\Delta_b$ is the bin width, and $N$ is the total sample size. If we want to assign **uncertainties** to each bin height (not often done, but its good practice) then we can quote\n",
    "\n",
    "$$ \\sigma_k = \\frac{\\sqrt{n_k}}{\\Delta_b N}$$"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-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.8.5"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected",
   "theme": "sky"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}