lectures/Week-2-Thu.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Maximum Likelihood and Applications in Astronomy\n",
    "\n",
    "A.J Connolly 2017 (based on lectures from M Juric with material from G. Richards 2016,Ivezic [Sections 4.0, 4.1, 4.2.1-6, 4.3], Bevington, and Leighly)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Statistical Inference\n",
    "\n",
    "Statistical *inference* is about drawing conclusions from data, specifically determining the properties of a population by data sampling.\n",
    "\n",
    "Three examples of inference are:\n",
    "1. What is the best estimate for a model parameter\n",
    "2. How confident we are about our result\n",
    "3. Are the data consistent with a particular model/hypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Frequentist vs. Bayesian Inference\n",
    "\n",
    "There are two major statistical paradigms which address the statistical inference questions: the classical, or **frequentist** paradigm, and the **Bayesian** paradigm.\n",
    "\n",
    "While most of statistics and machine learning is based on the classical paradigm, Bayesian techniques are being embraced by the statistical and scientific communities at an ever-increasing pace.\n",
    "\n",
    "This week we discuss the three main types of statistical inference from the classical (frequentist) point of view. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Some Terminology\n",
    "\n",
    "* We typically study the properties of some ***population*** by measuring ***samples*** from that population. The population doesn't have to refer to different objects. E.g., we may be (re)measuring the position of an object at rest; the population is the distribution of (an infinite number of) measurements smeared by the uncertainty, and the sample are the measurement we've actually taken.\n",
    "* A ***statistic*** is any function of the sample. For example, the sample mean is a statistic. But also, \"the value of the first measurement\" is also a statistic.\n",
    "* To conclude something about the population from the sample, we develop ***estimators***. An estimator is a statistic, a rule for calculating an estimate of a given quantity based on observed data: thus the rule (the estimator), the quantity of interest (the estimand) and its result (the estimate) can be distinguished. Sometimes *estimator* and the *estimate* are used interchangeably. Much of frequentist statistics concerns itself with development of different estimation rules -- different estimators -- and proving their properties.\n",
    "* Many of these rules have been developed for efficiency and speed of the estimation and make assumptions (usually reasonable) about the properties of the population, model, and any measurement uncertainties (e.g. Gaussianity).\n",
    "* There are ***point*** and ***interval estimators***. The point estimators yield single-valued results (example: the position of an object), while with an interval estimator, the result would be a range of plausible values (example: confidence interval for the position of an object).\n",
    "* Measurements have **uncertainties** (not errors) and we need to account for these (sometimes they are unknown). Data are not variables but fixed values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Maximum Likelihood Estimation (MLE)\n",
    "\n",
    "Let's talk about maximum likelihood estimation ($\\S 4.2$ in the textbook), which is relevant to both Bayesian and Frequentist approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Maximum Likelihood Approach\n",
    "\n",
    "Maximum likelihood estimation consists of the following conceptual steps:\n",
    "\n",
    "1. **Hypothesis**: Formulate a model, a *hypothesis*, about how the data are generated. For example, this could be a statement that the data are a measurement of some quantity that come with Gaussian random uncertainties (i.e., each measurement is equal to the true value, plus a deviation randomly drawn from the normal distribution). Models are typically described using a set of model parameters $\\boldsymbol{\\theta}$, and written as $\\boldsymbol{M}(\\boldsymbol{\\theta})$.\n",
    "2. **Maximum Likelihood Estimation**: Search for the \"best\" model parameters $\\boldsymbol{\\theta}$ which maximize the ***likelihood*** $L(\\boldsymbol{\\theta}) \\equiv p(D|M)$. This search yields the MLE *point estimates*, $\\boldsymbol{\\theta^0}$.\n",
    "3. **Quantifying Estimate Uncertainty**: Determine the confidence region for model parameters, $\\boldsymbol{\\theta^0}$. Such a confidence estimate can be obtained analytically (possibly with some approximations), but can also be done numerically for arbitrary models using general frequentist techniques, such as bootstrap, jackknife, and cross-validation.\n",
    "4. **Hypothesis Testing**: Perform hypothesis tests as needed to make other conclusions about models and point estimates. Possibly GOTO #1.\n",
    "\n",
    "While these steps represent a blueprint for the frequentist approach in general, the likelihood is just one of many possible so-called **objective functions** (also called **fitness functions**, or **cost functions**); other possibilities are explored briefly in §4.2.8 of the textbook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example: Measuring the Position of a Quasar\n",
    "\n",
    "Let's assume we wish to estimate the position $x$ of a quasar from a series of individual astrometric measurements.\n",
    "\n",
    "1. We adopt a model where the observed quasar does not move, and has individual measurement uncertainties \n",
    "2. We derive the expression for the likelihood of there being a quasar at position $x_0$ that gives rise to our individual measurements. We find the value of $\\hat x_0$ for which our observations are maximally likely.\n",
    "3. We determine the uncertainties (confidence intervals) on our measurement.\n",
    "4. We test whether what we've observed is consistent with our adopted model. For example, is it possible that the quasar was really a misidentified star with measurable proper motion?\n",
    "\n",
    "Note: in the text to come, I will use $\\mu$ instead of $x_0$ to denote the true position of the quasar. This is to avoid potential confusion with the first (or zeroth) measurement of $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Likelihood Function\n",
    "\n",
    "**Write down the likelihood function for this case**\n",
    "\n",
    "**Why did you choose the form you did? What assumption did you make?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The Likelihood Function\n",
    "\n",
    "If we know the distribution from which our data were drawn (or make a hypothesis about it), then we can compute the **probability** of our data being generated.\n",
    "\n",
    "For example, for the Gaussian distribution probablity of getting a specific value of $x$ is given by:\n",
    "\n",
    "$$p(x|\\mu,\\sigma) = \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x-\\mu)^2}{2\\sigma^2}\\right).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x124112e80>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAHjCAYAAACabpOIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xt0lNWh/vFnhzsoiWAFFIVAFAwoEG4FATEQElAQwtUL\nrbUgqEd7VFS8nC60Zylaq0t7qqjYqmgFJIgiEkKCAUFBLkEUAQkBRErBEggqKJDs3x+Z8IuYhEwy\nmT3vzPez1ixmJu/MPJE2PHn35TXWWgEAACD0RbkOAAAAgMqhuAEAAHgExQ0AAMAjKG4AAAAeQXED\nAADwCIobAACAR1DcAHiWMSbWGBNdk68zxkQbY7r4nw4AAo/iBiAojDEDjDG5xpgiY8wcY8wMY0yG\nMebeKr7fvZJ2SOoayNcZY+YaYx4vySxpvaTp5R0TDMaYib5brjFmQrA+F0DoMWzACyBYjDHTJd1r\nra3lexwtaaektdba5Cq8X5GkgdbaZYF6nTEmUVK+tXZjqcxdSucr45guks7xN0cls3aRNNdae7Ex\nprOkJjXxOQC8obbrAAAil7W2wBgzV9JE11lKVKYUlT7GGBMrKUvSqBqK1E3SYd/nbqyhzwDgEQyV\nAnCtqaTDxphU39DpRN+f2yXJGNPGGDPdGHOvb4jy9KHVtr7j830lUL7X3WKMebzkdWV8bnmvG2mM\nqbBInnZMgqQYSZOMMRN8XysZDm7sO/5eY8ycCt4vtqzv0RgzUtJoSW2MMVN8Z/rKen3JMHSq73FG\nMIdyAQQPZ9wAOGOMSZCUKukWa+18Y8xMSVbSLZISfEOp6yS1stZ+53tNvjHGWmuf8r3NQBWf7Wor\nab2vsEyX9IKkGGvtd8aYdcaYVGvt/FIfX9br5kp6QMVz4F4uJ3OX0sdYa9OMMZI0w1r7oe+Yt32f\nfaTUS+8v5/2iVTyPrszv0RjTRsXDsE+V9XpJstZm+T5znDGmibV2UHnHAvA2zrgBCDpjzAu+s1z3\nq3iu2Su+L+VLettau8tXssZIOlhSaHxekjSp1OMXrbVHrLU5vq+NstYWSOrqK22xKj4j1ua0GGW9\nLkfFRbFcFRxjSt2fLinJGNPa97ittXZXOW9Zme+xMjJVXEQz/XwdAA/hjBuAoLPW3lrJQ9uW8dwO\n/bKElViv4qFFSYrynUXLlJSn4iHZ8pR+XbVZa3OMMUslTTXGvKiKy6C/32N5dqr4bGW+n68D4CGe\nK27GmFhr7U7XOQAERXkFZkMFr1nnO8uWaa1tIknGmDKHKU9/XRXyVeRJSUsl2TMU1ap8j2UZqeKC\n2k0Sq06BMBX0oVLfZpYzfBNtH/ftk1TR8W/7JvoWGWMKJfm92SaAkGHOfMj/Z619WcULF0rvXZYk\n6bGSQ057yWgVD78mSIo2xjQ2xpQ1TFrW6+4r9TimEvFKH5On4gUEp34+WWuzfM83qehNKvE9nltR\nHmNMjG8Rw4sqPrs4qLxFDAC8z8UZt5clZVhrZ0qSbyXUqLKWuft+CP5HUhcV/8A/XME8EQAhzPdL\n2kjf/RckPVH6/8++8hGr4tWZ60r9TOgq6UVjTFtJByXNtta+4/va/ZLuN8Yk+R5Pt9ZuNMbsVPEZ\nq/UqLjQzJD1gjFnq28qjrNd95ss4UNI5vvJzSMXzxs4xxiRaa5edfozv/eZJekLFBav0IoJ5kjIq\n8Z+nzO/RtxBigIpL4YSSn5un/TedK2mAtfaIMeYl3+PKfCYADwrqBry+FWRrVbxC6ojvuRm+x2PL\nOH6GivcvepHhUQBeY4x5wY/5fABwRsEeKu2q4rNmpZfIr1fxb69laaPijTl3+Jbzx9Z0QACoDt+e\nbPf5rrjwtus8AMKLi+J2+oqnfJUzf8NaO8ha21TFxS5anP4HEPoSVDwU25hLUwEINE/s4+b74ddN\nxfM8mHQLIGRZa9OstU2ttbe5zgIg/AR7cUJZeyU1ke86fBXxXdMwS+WcnTPGBG+yHgAAQDVZa/1a\naS8Fv7it0y+LV1dVfqfvJqpgb6NgLrRAYE2bNk3Tpk1zHQNVwN+dt/H35238/XmX71J5fgvqUKnv\nUjGZJRdC9hkg3/5JvgslT/Tdj/VtDyDf4y6S1rIdCAAAiFQu5riNVvE1/Kb49nK631q72/e1JP3/\nFaYxkp4wxmw3xkyRFMuyegAAEMmCvgGvbyuQMguYtXZqqfs5kuKClQtu9e/f33UEVBF/d97G35+3\n8fcXeYK6AW9NMsbYcPleAABAeDPGVGlxgie2AwEAAICba5UCABAyWrdurd27d5/5QMBPrVq10q5d\nuwL6ngyVAgAimm/IynUMhKGK/rfFUCkAAECYo7gBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAA\ncGrnzp2uI/iloKDA2WdT3AAAQFCNHj1aUVFRioqKUq1atZwWIX9kZWUpLi5OL7/8srMMbMALAACC\npqCgQOeee65ycnJkrVVMTIxat27tOtYZHTlyRN26dVN+fr7THJxxAwAgwsXFxempp5469fjJJ59U\nXFxcjXzW/fffr+joaDVu3FidO3euVmnLzMzUoEGD1KRJk1Nn8ErO4gVa48aNFR0drSZNmgT8vf1B\ncQMAoJKMMSF3C4S8vDy1adPm1ON169apbdu2AXnvsj7r5ZdfVtu2bdWtW7cqz2+bNGmSkpOTFRcX\np5kzZ6pt27Z66aWXNG/ePL399tsBTh06GCoFAKCSwvHSWDt37pQxRgkJCaeey8vLU1JSUrmvmTx5\nsg4dOiTp5/9NSi7xZIxRUlKSJkyY8IvXZmRkSJKWLVumSZMmadCgQdq+fbtfmefNm6eXX35ZWVlZ\nuuqqqyRJBw8e1PTp05Wbm1sjuUMF1yoFAES0SL9WaVpamsaMGaPCwsJTz0VFRWnevHlKTU2t0c8u\nKChQkyZNtHTpUiUmJlb6dd26dVNBQcHPCl9BQYHOOeccvfTSSzVavOLi4jR58mRNmTLljMdyrVIA\nABBQpw+Tlgxdln6upkRHR2vAgAE6fPiwX687fPhwufl27NgRiGghi6FSAAAi2Nq1axUTE3Pq8fr1\n62WMUefOnct9Tekhx9P5O+SYn5//s2HaymjTps0vVncuXbpUxhiNGzcuKLldobgBABDBNmzYoEOH\nDiknJ0d5eXmaPn36qa+lpaVp5MiRv3jNjBkzqvRZO3fu1IYNG069Z05Ojrp3735qZWlWVpby8vI0\nceLECt/n/vvv16BBg7Rs2TIlJibq8OHDmjp1qu6//3516tSp3NdVNXcoYY4bACCiRfoct5L5YpL0\nxBNPqE2bNhowYIAGDhyoOXPmBPSzcnJyNHr0aBljNGnSJLVp0+Zn8+imTp2qtLS0Si1WmDlzpl58\n8cVT23NMnjxZI0aMCGje0goKCjR37lxNnjxZCQkJeuKJJ844L68m5rhR3AAAES2Si1taWpoyMzP1\nwgsvuI5yyvz582t8UUSw1ERxY6gUAIAINWfOnIDtBYfgYFUpAAARyhijsWPHuo5xys6dO/1eqBBp\nGCoFAES0SB4qRc1iHzcAAIAIRnEDAADwCIobAACAR1DcAAAAPILtQAAAEa1Vq1ZsiYEa0apVq4C/\nJ6tKAQAAgoxVpQAAAGGO4gYAAOARFDcAAACPoLgBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAA\ngEdQ3AAAADyC4gYAAOARFDcAAACPoLgBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAAgEdQ3AAA\nADyC4gYAAOARFDcAAACPoLgBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAAgEdQ3AAAADyC4gYA\nAOARFDcAAACPoLgBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAAgEdQ3AAAADyC4gYAAOARFDcA\nAACPoLgBAAB4BMUNAADAIyhuAAAAHkFxAwAA8AiKGwAAgEdQ3AAAADyC4gYAAOARFDcAAACPqB3s\nDzTGREt6QlKupKaSMq21WZV43XRJ1lr7QA1HBAAACElBL26SXpaUYa2dKUnGmFxjzChr7cbyXmCM\nGSBplKS3g5QRAAAg5AR1qNQYkyBppKS5pZ7OlFTuWTRjTKykNpLW12w6AACA0BbsOW5dJR221h4p\n9dx6SQMreM0t1tqXJZkaTQYAABDiXBS3/NOey5cUU9bBxpjHJT1W06EAAAC8wMUct0oxxkyUNNda\n+11lXzNt2rRT9/v376/+/fsHPhiAsHfy5Elt2bJFu3fv1okTJ9SkSRN17txZ0dHRrqMB8Kjs7Gxl\nZ2dX+32Mtbb6aSr7YcVlbLq1tmlFz/mez5UUe/pbqHhlaa0y3tsG83sBEH42bdqkp59+WgsXLtSv\nfvUrtWnTRnXr1tW3336rTZs2KT4+Xrfccouuv/56NWjQwHVcAB5mjJG11u9pYMEeKl2nXw6LdlXx\nAoWfsdbGWWtrldwkpUl6oqzSBgDVsX//fo0bN07JycmKj4/Xxo0btXXrVn3wwQdasGCBVq1apYMH\nD2ratGmaP3++Lr30Ur3zzjvil0UAwRbU4matzZGUaYxJLfX0AEn3ScXbfvjOwAFAUHzwwQe6/PLL\n1bp1a+Xl5em+++7ThRde+Ivj6tatq8GDB2vRokX6+9//rgceeEA333yzjh496iA1gEjl4soJoyUl\nGWOmGGNekHS/tXa372tJKn+FKb/aAgiop59+WhMmTNA777yj6dOnV3r4MzExUevWrdNPP/2kPn36\n6MCBAzWcFACKBXWOW01ijhsAfzz66KOaPXu20tPTddFFF1XpPay1mjZtmubMmaPMzEy1bNkywCkB\nhKuqznEL2VWlAFBTpk+frtmzZ+vDDz9Us2bNqvw+xhg98sgjatiwoZKSkrRy5Uo1bdr0zC8EgCri\njBuAiPLWW29p6tSpWr16tVq0aBGw97333nv18ccfKzMzkxWnAM6oqmfcKG4AIsbatWs1ZMgQZWVl\n6fLLLw/oexcVFem6665To0aN9Morr8gYLvYCoHxe2Q4EAJwoKCjQuHHjNGPGjICXNkmKiorSK6+8\nojVr1mjmzJkBf38AkDjjBiACWGt1/fXXKyYmRi+88EKNftbWrVvVp08frVq1Su3atavRzwLgXSxO\nAIByzJo1S59//rnWrl1b45/Vvn17PfLII/rtb3+rlStXqnZtfswCCBzOuAEIa/v379dll12mjIwM\nde7cOSifWVRUpKSkJA0cOFAPPPBAUD4TgLewOIHiBqAM119/vS688EI98cQTQf3c3bt3q1u3bvro\no4/Uvn37oH42gNBHcaO4AThNenq6brvtNn3xxRdq2LBh0D//6aef1pIlS5Sens4qUwA/w6pSACjl\n+PHj+q//+i89//zzTkqbJN1xxx3as2eP3n33XSefDyD8UNwAhKW//e1vateunVJSUpxlqFOnjv76\n17/q7rvv1rFjx5zlABA+GCoFEHby8/PVvn17ZWdnKz4+3nUcjRw5Ut27d9fUqVNdRwEQIpjjRnED\n4HP33Xfr6NGjmjFjhusokqRt27apT58++uqrr3TOOee4jgMgBFDcKG4AVLyaMyEhQV9++WW1LiAf\naBMmTNB5552nxx57zHUUACGA4kZxAyBp8uTJatKkScgVpD179qhz587avHmzmjdv7joOAMcobhQ3\nIOKVnG3btm2bzj33XNdxfuGuu+5SYWGhnnvuOddRADhGcaO4ARFv0qRJatq0acidbSuxb98+dejQ\nQVu3btV5553nOg4AhyhuFDcgooX62bYSt912m6Kjo/X444+7jgLAIYobxQ2IaHfddZdq1aqlp556\nynWUCu3atUtdu3bVjh07FBMT4zoOAEcobhQ3IGIdPnxYbdq00WeffaYLL7zQdZwz+u1vf6tLLrlE\nDz30kOsoAByhuFHcgIj15JNPatOmTXrjjTdcR6mULVu2qH///tq9e7fq16/vOg4AB7hWKYCIdPz4\ncT333HO65557XEeptEsvvVRdunTR7NmzXUcB4DEUNwCeNmfOHLVv315dunRxHcUvf/jDH/Tss8+K\nkQIA/qC4AfAsa62eeuopTZkyxXUUvyUnJ+vo0aP66KOPXEcB4CEUNwCe9dFHH+n48eNKTk52HcVv\nUVFRp866AUBlsTgBgGddd9116tWrl+68807XUark+++/V+vWrbVu3Tq1bt3adRwAQcSqUoobEFH2\n79+v9u3ba+fOnZ7eD23KlCmKiorSk08+6ToKgCCiuFHcgIjy2GOPKS8vTzNnznQdpVpyc3PVu3dv\n7dmzR/Xq1XMdB0CQsB0IgIhRWFioF198UbfddpvrKNUWFxenyy67TO+++67rKAA8gOIGwHMWL16s\n5s2bKyEhwXWUgLjlllv00ksvuY4BwAMobgA854UXXtCtt97qOkbADB8+XJs2bVJubq7rKABCHHPc\nAHjK3r171bFjR+3du1cNGzZ0HSdgpkyZotq1a2v69OmuowAIAua4AYgIb7zxhkaNGhVWpU2SJk6c\nqFdffVXHjx93HQVACKO4AfAMa61effVV3XTTTa6jBFy7du3Url07LVy40HUUACGM4gbAMz799FMV\nFhaqd+/erqPUiJtvvlmvv/666xgAQhjFDYBnlJxtM8bvaSGekJqaquXLl+vbb791HQVAiKK4AfCE\nH3/8UXPnztX48eNdR6kxZ599toYOHaq33nrLdRQAIYriBsATFixYoK5du+rCCy90HaVG/eY3v2G4\nFEC5KG4APCFcFyWcLjExUfv27dPmzZtdRwEQgihuAELe3r179emnn2r48OGuo9S4WrVq6cYbb9Ss\nWbNcRwEQgihuAELenDlzNHz48LDbu608v/nNb/TGG2+osLDQdRQAIYbiBiDkvfXWW7ruuutcxwia\nDh06qFmzZsrOznYdBUCIobgBCGm5ubn6+uuvddVVV7mOElTjxo3TnDlzXMcAEGIobgBC2uzZszV6\n9GjVrl3bdZSgGjNmjObPn68TJ064jgIghFDcAIQsa23EDZOWaNWqldq1a6fMzEzXUQCEEIobgJD1\n+eef6/vvv1evXr1cR3Fi7Nixmj17tusYAEKIsda6zhAQxhgbLt8LgGIPPPCACgsL9eSTT7qO4sS+\nffsUHx+vffv2qX79+q7jAAggY4ystX5fv48zbgBCkrVWs2fPjshh0hItWrRQly5dtHjxYtdRAIQI\nihuAkLRmzRrVq1dPnTt3dh3FKVaXAiiN4gYgJM2ZM0fjxo2TMX6PJISV1NRULV68WD/88IPrKABC\nAMUNQMix1iotLU2jRo1yHcW5c889V71799Z7773nOgqAEEBxAxBy1q1bpwYNGqhDhw6uo4SEMWPG\nKC0tzXUMACGA4gYg5KSlpWnkyJERP0xaYtiwYVq6dKmOHj3qOgoAxyhuAEJKyTDpyJEjXUcJGU2b\nNlX37t2VkZHhOgoAxyhuAELK559/rpMnTyohIcF1lJCSmpqq+fPnu44BwDGKG4CQkpaWptTUVIZJ\nTzN8+HC9//77On78uOsoAByiuAEIKfPnz1dqaqrrGCHn/PPPV/v27ZWdne06CgCHKG4AQsZXX32l\ngwcPRuy1Sc+E4VIAFDcAISMtLU0jRoxQVBQ/msoyYsQILViwQIWFha6jAHCEn44AQgarSSvWtm1b\nNW/eXJ988onrKAAcobgBCAm7du3S7t271a9fP9dRQhrDpUBko7gBCAkLFizQsGHDVLt2bddRQlpJ\ncbPWuo4CwAGKG4CQ8N5772n48OGuY4S8Dh06qG7dusrJyXEdBYADFDcAzh06dEjr1q3TgAEDXEcJ\necYYjRgxQu+8847rKAAcoLgBcG7x4sXq37+/GjZs6DqKJwwdOlQLFy50HQOAAxQ3AM699957GjZs\nmOsYntGrVy/t3btXX3/9tesoAILM7+JmjGltjEn03TobYxrXRDAAkeH48eNasmSJrrnmGtdRPKNW\nrVoaMmQIZ92ACFSp5VvGmJGSxkqKlZQn6ZDvS00kxRhjrKQNkh631h6piaAAwtOKFSvUrl07NW/e\n3HUUTxk6dKhmzpyp22+/3XUUAEFkKlpSboyJlvSApFxr7cwK38iYLioud59aa4O+yZAxxrI8HvCe\nP/zhD2rWrJkefPBB11E85ciRI2rZsqX27t2rs88+23UcAH4yxshaa/x9XblDpcaYWEkTrbVTz1Ta\nJMlam2OtnSpppzFmgr9BAEQea63ee+89DR061HUUz2ncuLF69eqljIwM11EABFG5xc1au9Na+5S/\nb+grcGcsegDwxRdfSJI6duzoOIk3sboUiDysKgXgTMlqUmP8Hi2AiovbokWLuOg8EEECUtyMMYl+\nHBttjJlhjJlijHncGFPujpvGmFhjTIYxJt8Yk+ubRwcgTLANSPW0atVKLVq00OrVq11HARAklVmc\nMLES7zPWWtu9Uh9ozFxJGSXDqcaYXEmjrLUbyzj2BUn3STKS5kmKtdZeXM77sjgB8JB9+/YpPj5e\nBw4cUJ06dVzH8ayHH35YJ0+e1PTp011HAeCHqi5OqLC4+d44V8XbfuRXcFistbbWGT/MmARJayWd\nU7JtiDFmhu/x2NOObVx6axHfliRTyyuIFDfAW2bOnKnMzEzNnj3bdRRPW7NmjW6++WZt3rzZdRQA\nfqhqcavMPm6TJK2z1hZU8OGVOSsnSV0lHT5tr7f1kn7xq2IZ+8F1V+XO/gHwgPfee0/jxo1zHcPz\nunfvroMHD2rHjh1q27at6zgAatgZ57hZa7MknelU1txKfl5X/fLMXb6kmPJeYIwZ4Dvrd6+kbpX8\nHAAh7Mcff1R2drZSUlJcR/G8qKgoXX311awuBSJEpRYnnOlqCBWdjasuX3EcqOI5bi8aYzrX1GcB\nCI6PPvpIl112mZo0aeI6SlgYOnSo3nvvPdcxAARBpS55FUDrJY0+7bkmkg5X9CJr7S5JY31n3rpJ\n+sVCBkmaNm3aqfv9+/dX//79q54UQI1JT0/nbFsAJSUlafz48Tpy5IgaN+by0UAoys7OVnZ2drXf\n54yLEwLJt53HutILGcpbnFDO6zMkTbfWLivjayxOADyiQ4cO+sc//qEePXq4jhI2kpOTNWnSJKWm\nprqOAqASAn7Jq5pgrc2RlGmMKf2TZYCKt/womc820Xc/uvQeb76tSQ6VVdoAeMeePXu0f/9+de3a\n1XWUsDJ48GB98MEHrmMAqGEurpwwWlKSbwPeFyTdb63d7ftakorns0nFQ6gvGmO2+659OqAyZ+UA\nhLYlS5Zo0KBBqlXrjDsIwQ9DhgzR4sWLxcgDEN78nuNmjJlQavPcU/cry7fQ4dZyvja11P2dkuL8\nzQcgtKWnp3NR+Rpw8cUXq0GDBtq0aZM6derkOg6AGlKVM25J5dwHgAqdPHlSWVlZGjRokOsoYccY\noyFDhjBcCoS5qhQ3U859AKjQmjVr1Lp1a7Vo0cJ1lLA0ePBgLV682HUMADWoKsXNlnMfACrENiA1\nq3///srJydGhQ4dcRwFQQ6p7xg0AKi09PV3JycmuY4StBg0aqF+/flq6dKnrKABqiItVpQAi0Lff\nfquvvvpKvXv3dh0lrDFcCoQ3ihuAoFi6dKmuuuoq1a1b13WUsFayLUhRUZHrKABqQHXnuAFApTC/\nLTjatGmjmJgY5eTkuI4CoAZwxg1AjSsqKtKSJUuY3xYkJWfdAISfqhS3Q+XcB4Aybdy4UTExMYqN\njXUdJSJw+SsgfPld3Ky1k8u6DwDlWbJkCcOkQdSvXz998cUXOnjwoOsoAAKMoVIANY75bcFVr149\nXXXVVcrIyHAdBUCAUdwA1KiCggJt2LBBV155pesoEYXhUiA8UdwA1Khly5apV69eatiwoesoEWXw\n4MFKT09nWxAgzFDcANQo5re50apVK5177rnasGGD6ygAAqjSxc0YU2iMSa3JMADCi7WW+W0OpaSk\naMmSJa5jAAggf864nbpGqTFmgjHmcWNMpxrIBCBMbNu2TYWFhbr00ktdR4lIycnJFDcgzPg7VNrE\nGJMr6SVJ90vaYIx5LPCxAISDkrNtxpgzH4yA69evn3JyclRQUOA6CoAA8be43S/pCUnn+G5jJY0x\nxowIdDAA3sf8NrcaNmyoXr16admyZa6jAAgQf4rbYUn3W2tfttYW+G7zJHWTxEa8AH7m2LFjWrly\npQYMGOA6SkRjuBQIL/4Ut3WSBp7+pLX2cODiAAgXK1asUKdOnRQTE+M6SkQrKW7WWtdRAASAP8Xt\nSUmTjDFXlX7SGJMoqUlAUwHwPFaThoYOHTroxIkT2r59u+soAAKg0sXNWpspaaakLGPMdmPMEmPM\nWklLJb1YUwEBeBPz20KDMYbhUiCM+LU4wVo7SdIYSQWS2ko6JGmQtXZmDWQD4FG7d+/Wf/7zHyUk\nJLiOAjHPDQgnJlzmPRhjbLh8L4DXvfTSS1q+fLnefPNN11EgKT8/X61bt9a3336revXquY4DQMVn\nw621fu+VxCWvAAQc89tCS5MmTRQfH6+VK1e6jgKgmihuAALqxIkTWrZsmQYNGuQ6CkphuBQIDxQ3\nAAG1evVqtW3bVs2aNXMdBaVQ3IDwUO3iZoyZ6Ftlmuq7ta5+LABelZ6eruTkZNcxcJoePXpoz549\n2rdvn+soAKohEGfc8iSlWWvnW2vnq9TF6AFEHua3habatWtrwIABysjIcB0FQDVUu7hZa7OstVNL\nPd5Z3fcE4E0HDhzQjh071KtXL9dRUIbk5GSlp6e7jgGgGiosbsaYWGNM5zKeT2RYFMDpMjIylJiY\nqDp16riOgjIkJydr6dKlKiwsdB0FQBVVWNx8Z8/WlxQ0X5HbruKrJUyWNM8Y80JNhwTgDcxvC20X\nXnihzjvvPG3YsMF1FABVVJmhUiOp5CrRb/v+bGKtHWSt7SbpJWNMao2kA+AZRUVFysjIoLiFOFaX\nAt5W2Tlu1hgTKylB0iRrbcGpL1ibI6lNTYQD4B05OTlq2rSpWrdu7ToKKkBxA7ytMsUtTcXDopMk\nHbLWLivjmEkBTQXAc1hN6g1XXnmlNm7cqIKCgjMfDCDknLG4WWtHS5qq4uHS9ad/3RgzUmwBAkQ8\n5rd5Q4NmayAQAAAgAElEQVQGDdS7d29lZWW5jgKgCqp9kXnfEGqMb8jUGS4yD7hTUFCgli1b6sCB\nA2rQoIHrODiDp59+Wtu2bdOLL77oOgoQsZxdZN5au9N1aQPgVlZWlq644gpKm0eUzHPjl13Ae7hW\nKYBqY36bt8THx+vkyZP66quvXEcB4CeKG4BqsdYyv81jjDGsLgU8iuIGoFq2bNkiY4zat2/vOgr8\nwOWvAG/yu7gZYyaUdR9AZFqyZIlSUlJkDIvLvWTgwIFauXKlfvzxR9dRAPihKmfcksq5DyACMUzq\nTU2aNFGHDh20cuVK11EA+KEqxc2Ucx9AhDl69Kg+/vhjDRgwwHUUVAHz3ADvqUpxs+XcBxBhli9f\nri5duig6Otp1FFQBxQ3wnuqecQMQwUrmt8Gbunfvrm+++UZ79+51HQVAJbGqFECVMb/N22rXrq2B\nAwcqIyPDdRQAlURxA1Alu3btUn5+vrp06eI6CqqB4VLAW6o7xw1AhFqyZImSk5MVFcXvf16WnJys\npUuXqrCw0HUUAJXAT1wAVcJlrsJDy5Yt1aJFC61bt851FACVUJXidqic+wAixPHjx7Vs2TIlJbGV\nYzhguBTwDr+Lm7V2cln3AUSOTz75RBdffLHOO+8811EQABQ3wDsYKgXgN7YBCS/9+vXTpk2bdOgQ\ngyhAqKO4AfAb89vCS/369dWnTx9lZWW5jgLgDChuAPzy73//W3l5eerZs6frKAgghksBb6C4AfBL\nRkaGBgwYoDp16riOggAqKW7WsuMTEMqqXdyMMSMDEQSANzC/LTy1b99exhht2bLFdRQAFah0cTPG\n5BpjlhhjHjfGpBpjOkuStTbNGDOn5iICCBWFhYXKyMjgMldhyBjDcCngAf6ccXtJUoGkJEnzJG0w\nxhQaYwolDayJcABCy4YNG/SrX/1KF110kesoqAEUNyD0Vbq4WWuftNaOsdZ2s9ZGSWoraaykHEkD\naioggNDBatLwNmDAAK1atUrHjh1zHQVAOao8x81au9NaO0/S45LaBC4SgFDF/LbwFhMTo06dOmnF\nihWuowAohz9z3Kb45ra1Pu1LeZKeCGQoAKHn0KFD2rRpk/r16+c6CmoQw6VAaKvtx7GTVXxmzRpj\npOLCdtj3XH7gowEIJVlZWbriiitUv35911FQg5KTk/W73/3OdQwA5fBnjlucb25bnKQxKl6ssFPF\nF5ovMMa8YIyZY4xJrZmoAFxifltk6Nq1q/bv3689e/a4jgKgDFW5yPxOa22atfbPvsUKcdbabpIy\nJRkVL1gAEEastcxvixC1atVSUlISw6VAiArYlRN8ZW6MtZbiBoSZL7/8UrVr19Yll1ziOgqCgHlu\nQOjiklcAzig9PV3JycnyzW9FmBs0aJCysrJ08uRJ11EAnKbCxQm+y1mdU4X3zbfWzq9aJAChJj09\nXbfffrvrGAiS888/Xy1bttSnn36q3r17u44DoJQzrSrNC0oKACHr+++/1+rVqzV/Pr+LRZKUlBQt\nWbKE4gaEmAqHSq21OaVvksac/lx5tyDlB1DDsrOz1a1bN5199tmuoyCImOcGhCZ/57h1Pf0JY0yi\nMSYxQHkAhJjFixdr8ODBrmMgyPr06aMvv/xSBw8edB0FQCl+FTdr7SBjzARjTGtjTGNfYRsk6cWa\niQfAJWstxS1C1atXT3379lVmZqbrKABK8au4+TbXHaTiuW+HVFzYolW8IS+AMLN9+3YdP35cHTt2\ndB0FDpTMcwMQOvwdKn1S0g5JSZIuVvHVE+73Z06bMSbaGDPDd+3Tx40xAyo4NtYYs84YU+T7s4uf\neQFUw+LFi5WSksI2IBGqZJ6btdZ1FAA+/ha3SdbaB6y1WdbaPGvtnyUlGWM6+/EeL0taZ619ylr7\ngKQXy3q9MSZaxRev/72kBN/T64wxjf3MDKCKuMxVZLv44otVt25dbd682XUUAD6mur9JGWNiJO2w\n1jatxLEJktZKOsdae8T33Azf47GnHTtS0npr7S7f42gVD88OtNYuK+O9Lb8VAoFz7NgxnXfeedqz\nZ49iYmJcx4EjkydPVlxcnKZMmeI6ChBWjDGy1vo9nFHtKydYaw+r8nPcuko6XFLafNZLGljG+6aV\nlDbf4wJJVtK6qqcFUFnLly9X586dKW0RjnluQGgpt7j55pelVuZNrLVZpV7X2BgzoZxDu0rKP+25\nfEln/JfBGDNQUuZppQ9ADWE1KSQpMTFRq1ev1tGjR11HAaAKipu1dqekAt8CgkrNYTPGTJT0gLV2\nZqACljJR0ugaeF8AZUhPT6e4QY0bN1aXLl20fPly11EA6AyXvPKdScsyxkw3xjyh4m1Adkg6rOIz\nZU0ktfXdrKQXrbUvV/CW6/XL8tXE937lMsbcK+nxM51tmzZt2qn7/fv3V//+/Ss6HEA58vLyVFBQ\noE6dOrmOghCQnJxMkQeqKTs7W9nZ2dV+H78WJxhjYlU8H61kaPOwpLzSQ6VneH0XFa8orVXquTIX\nJ5T6+kgVL37YeIb3ZnECECB/+9vftHbtWr366quuoyAErF+/XjfccIO2bt3qOgoQNqq6OOFMF5n/\nGd/waUVn1M70+hxjTKYxJtVaW3LF6gHyLU7w7enWpuSsnW9eWxNJG3yl8RwVXy91alUzADiz9PR0\n3XDDDa5jIER06dJFhw4dUl5entq0aeM6DhDRqr2qtApGq3jvtynGmBdUvIHvbt/XkvTzEpchaYak\nXN9traRvgx8ZiBw//vijli9frqSkJNdRECKioqI0ePBgffDBB66jABHP30teTTDGHPTd5hhjWvn7\ngdbaI9baW30b8N5a6sybrLVTS4ZMfZv8Rllra512+4u/nwmg8lauXKkOHTqoadMzbs2ICHL11VdT\n3IAQUJVLXk2XNFWSUfEQ5u8DngqAM2wDgrIkJSXpo48+YlsQwDF/i9tca+2frbUvW2vHqHg1aXdj\nTGINZAPgAKsHUZaYmBglJCToww8/dB0FiGh+F7fSm/Jaaw9bayeruMAB8Livv/5aBw4cUNeuXV1H\nQQhiuBRwz9/iNk/Sy8aYtb7FBSVz3HYEOBcAB9LT05WcnKyoKBfrlhDqhgwZog8++EBsvQS44+9P\n58cl3SIpS9I4STuNMQclPWGMGWGMaV3Zy2QBCD2LFy9WSkqK6xgIUR06dFBRUZG2bNniOgoQsfza\ngLfMNyjea63kliDJlt5gN1jYgBeonp9++knnnXeeduzYoXPPPdd1HISoW2+9VW3bttWUKVNcRwE8\nraob8FZ7PMRam+nbxqObtTZK1digF4A7y5cvV8eOHSltqFDJcCkANwI+kcW3WAGAx7z//vu6+uqr\nXcdAiEtMTNTatWtVUFDgOgoQkZiBDEDWWr3//vu65pprXEdBiGvUqJGuuOIKZWZmuo4CRCSKGwBt\n3bpVJ06c0GWXXeY6CjyA4VLAHYobAC1atEjXXHONjPF7niwiUElxKyoqch0FiDgUNwAMk8IvcXFx\naty4sTZu3Og6ChBxKG5AhDt8+LA2bNigq666ynUUeAhXUQDcoLgBEW7JkiXq16+fGjZs6DoKPGTI\nkCFatGiR6xhAxKG4ARGObUBQFX379tWXX36p//znP66jABGF4gZEsMLCQqWnp1Pc4Ld69eopMTFR\n6enprqMAEYXiBkSwNWvWqEWLFrroootcR4EHXX311Xr//fddxwAiCsUNiGAl24AAVXHNNddoyZIl\nOn78uOsoQMSguAERjG1AUB3NmzdXu3bttGLFCtdRgIhBcQMi1J49e7R371717NnTdRR42LBhw/Te\ne++5jgFEDIobEKEWLVqkwYMHq1atWq6jwMNKipu11nUUICJQ3IAIxTYgCIQOHTooKipKn3/+ueso\nQESguAER6Pvvv9eKFSuUkpLiOgo8zhjDcCkQRBQ3IAJlZGSoZ8+eiomJcR0FYYDiBgQPxQ2IQO++\n+66uvfZa1zEQJvr27avc3Fz961//ch0FCHsUNyDCnDx5UosWLdKwYcNcR0GYqFOnjlJSUtiMFwgC\nihsQYVatWqVWrVpxtQQEFMOlQHBQ3IAIwzApakJKSopWrFihH374wXUUIKxR3IAIYq3VggULKG4I\nuJiYGPXo0UNLly51HQUIaxQ3IIJ88cUXKioq0uWXX+46CsIQw6VAzaO4ARGkZJjUGOM6CsLQ0KFD\n9f7776uwsNB1FCBsUdyACML8NtSk2NhYNW/eXGvWrHEdBQhbFDcgQuzdu1d5eXnq27ev6ygIY8OG\nDdO7777rOgYQtihuQIR47733NGTIENWpU8d1FISxESNGaP78+Vx0HqghFDcgQjBMimBISEjQiRMn\n9MUXX7iOAoQlihsQAY4cOaKPP/5YycnJrqMgzBljlJqaqvnz57uOAoQlihsQARYtWqS+ffvq7LPP\ndh0FESA1NVVpaWmuYwBhieIGRIC0tDSNHDnSdQxEiN69e+vbb7/V9u3bXUcBwg7FDQhzP/zwg5Yu\nXcr8NgRNVFSUhg8fznApUAMobkCYS09PV48ePdS0aVPXURBBmOcG1AyKGxDmGCaFC/3791dubq72\n7NnjOgoQVihuQBj76aeftHjxYo0YMcJ1FESYOnXqaOjQoXrnnXdcRwHCCsUNCGNLly7V5ZdfrmbN\nmrmOggg0cuRIhkuBAKO4AWFs3rx5DJPCmaSkJG3cuFEHDhxwHQUIGxQ3IEydOHFCCxcuVGpqquso\niFD169dXcnIy1y4FAojiBoSpDz/8UJdccolatmzpOgoiGMOlQGBR3IAwxWpShILBgwfr448/1qFD\nh1xHAcICxQ0IQ4WFhVqwYAHFDc6dffbZGjhwIKtLgQChuAFhaPny5WrZsqViY2NdRwE0duxYzZkz\nx3UMICxQ3IAwNHv2bF133XWuYwCSpKuvvlpr1qzRt99+6zoK4HkUNyDMHD9+XPPnz9eYMWNcRwEk\nSY0aNdLgwYOVlpbmOgrgeRQ3IMxkZmaqffv2uuiii1xHAU5huBQIDIobEGZmz56tcePGuY4B/ExK\nSoo2btyoffv2uY4CeBrFDQgjx44d08KFCzVq1CjXUYCfqV+/voYNG6Z58+a5jgJ4GsUNCCOLFy9W\nQkKCmjdv7joK8Atjx47V7NmzXccAPI3iBoQRhkkRygYOHKht27Zpz549rqMAnkVxA8LEd999pyVL\nlnBtUoSsunXravjw4Zo7d67rKIBnUdyAMLFw4UL16dNHTZs2dR0FKNe4ceMYLgWqgeIGhAmGSeEF\n/fv31zfffKNt27a5jgJ4EsUNCAP5+flavny5rr32WtdRgArVrl1b1113nd58803XUQBPorgBYWDO\nnDlKSUlR48aNXUcBzujGG2/UG2+8IWut6yiA51DcgDDw+uuv6ze/+Y3rGECldOnSRQ0aNNCqVatc\nRwE8h+IGeNz27duVl5enQYMGuY4CVIoxRuPHj9esWbNcRwE8h+IGeNysWbN03XXXqU6dOq6jAJV2\n/fXXa968efrpp59cRwE8heIGeFhRUZFmzZrFMCk856KLLtLll1+uRYsWuY4CeArFDfCwVatWqVGj\nRurSpYvrKIDfGC4F/EdxAzzs9ddf1/jx42WMcR0F8NvIkSO1bNky5efnu44CeAbFDfCoY8eOKS0t\nTTfccIPrKECVREdHKyUlhUtgAX6guAEetXDhQiUkJKhly5auowBVNn78eL3++uuuYwCeQXEDPIpF\nCQgHycnJ2rlzp7Zs2eI6CuAJJlx2rjbG2HD5XoAz2bdvn+Lj47Vnzx6dddZZruMA1TJ16lQVFhbq\nz3/+s+soQNAYY2St9XuCMmfcAA967bXXNHLkSEobwsLNN9+s119/XcePH3cdBQh5niluxpho1xmA\nUGCt1SuvvKIJEya4jgIExCWXXKJ27drp/fffdx0FCHlBL27GmGhjzAxjzBRjzOPGmAFnOH6AMSZX\n0sQgRQRC2ooVK1SvXj317NnTdRQgYH7/+9/rlVdecR0DCHkuzri9LGmdtfYpa+0Dkl40xnQu60Bj\nTGNJ6yQ1CWZAIJTNnDlTEyZMYO82hJVRo0bpk08+0d69e11HAUJaUIubMSZB0khJpTftyZT0QFnH\nW2uPWGsLJLE7IyDp8OHDWrhwoW688UbXUYCAatSokcaMGaNXX33VdRQgpAX7jFtXSYettUdKPbde\n0sAg5wA86Z///KeSk5N17rnnuo4CBFzJcGlRUZHrKEDIclHcTj97li8pJsg5AE8qGSYFwlG3bt10\n1llnKTs723UUIGTVdh0gkKZNm3bqfv/+/dW/f39nWYBAW79+vfLz8zVgQIXreQDPMsZowoQJmjlz\nphITE13HAQIqOzs7IL+UBHUDXmPMREnTrbVNK3qujNflSpphrX2qgmPYgBdhbeLEiYqNjdWDDz7o\nOgpQYw4dOqQ2bdpo27ZtOu+881zHAWqMVzbgXadfDot2VfECBQDlOHTokObNm6ff//73rqMANeqc\nc85RamoqW4MA5QhqcbPW5kjKNMaklnp6gKT7pFN7trFfG3Ca1157TYMHD1azZs1cRwFq3O23364Z\nM2aosLDQdRQg5LjYx220pCTfBrwvSLrfWrvb97UklVph6tusd6KkWEljjTFMekDEKSoq0vPPP6/b\nb7/ddRQgKBISEtSiRQstWrTIdRQg5HCReSDELV26VFOmTNHGjRvZdBcRY9asWXrzzTeVnp7uOgpQ\nI7wyxw2An55//nnddtttlDZElNGjR2vDhg3avn276yhASOGMGxDC9uzZo86dO2v37t0666yzXMcB\ngmrq1Kk6ceKE/vKXv7iOAgRcVc+4UdyAEPbwww/ryJEjeu6551xHAYJu586d6t69u77++ms1bNjQ\ndRwgoChuFDeEmWPHjqlVq1b66KOP1K5dO9dxACeGDx+ulJQUTZ482XUUIKCY4waEmVmzZunXv/41\npQ0R7e6779YzzzzD9UsBH4obEIKKior0zDPP6O6773YdBXCqb9++Ovvss/XBBx+4jgKEBIobEILS\n09PVoEEDXXnlla6jAE4ZY3TXXXfpmWeecR0FCAkUNyAEPf3007r77rvZAgRQ8dYg27Zt08aNG11H\nAZyjuAEh5rPPPtOWLVs0ZswY11GAkFC3bl3dcccdnHUDxKpSIOTcdNNNat++vaZOneo6ChAyDh06\npLZt2+qLL77Q+eef7zoOUG1sB0JxQxj45ptvdPnllys3N1dNmjRxHQcIKXfccYcaNWqk6dOnu44C\nVBvFjeKGMHDXXXcpKiqKneKBMuzevVsJCQnKzc3VOeec4zoOUC0UN4obPO7bb79Vu3btGAoCKvC7\n3/1Obdq00f/8z/+4jgJUC8WN4gaPe+ihh5Sfn68XXnjBdRQgZG3btk19+/ZVXl4e1++Fp1HcKG7w\nsIKCArVt21Zr165VbGys6zhASBszZox69uype+65x3UUoMoobhQ3eNhjjz2mrVu36vXXX3cdBQh5\nGzdu1JAhQ5SXl6f69eu7jgNUCcWN4gaPOnr0qGJjY/Xhhx8qPj7edRzAE66++moNHTqUi8/Dsyhu\nFDd41NNPP61Vq1YpLS3NdRTAMz7++GNdf/312rZtm+rVq+c6DuA3ihvFDR703Xff6eKLL1ZmZqY6\nduzoOg7gKYMHD9Y111yj22+/3XUUwG8UN4obPOh///d/tWXLFr355puuowCes379eg0dOlS5ublq\n2LCh6ziAXyhuFDd4zKFDh3TxxRfrk08+0cUXX+w6DuBJI0eOVK9evTRlyhTXUQC/UNwobvCYhx56\nSPv379fMmTNdRwE8a/PmzUpMTNT27dvVuHFj13GASqO4UdzgIQcOHNCll16qnJwcXXTRRa7jAJ52\n44036pJLLtEf//hH11GASqO4UdzgIXfddZdOnjypv/71r66jAJ6Xm5urX//619q2bZuaNm3qOg5Q\nKRQ3ihs8ouQfmc2bN6tZs2au4wBh4bbbblOdOnX07LPPuo4CVArFjeIGjxg5cqS6deumBx54wHUU\nIGwcOHBA8fHxWrVqldq1a+c6DnBGFDeKGzxgxYoVGj9+vLZu3aoGDRq4jgOElSeffFKrVq3Su+++\n6zoKcEZVLW5RNREGwC8VFRXp7rvv1vTp0yltQA2488479fnnn2vZsmWuowA1huIGBMmbb76p2rVr\na9y4ca6jAGGpfv36euKJJ3TPPfeosLDQdRygRlDcgCD4/vvv9eCDD+rpp5+WMX6fGQdQSaNGjVLD\nhg316quvuo4C1AjmuAFBcN999+lf//qX3njjDddRgLC3YcMGDRkyRJs3b2Z7EIQsFidQ3BCiPv/8\ncyUmJuqLL75g+w8gSO68804dO3ZML7/8susoQJkobhQ3hKCioiL169dPN954oyZPnuw6DhAxCgoK\nFB8fr7ffflu9e/d2HQf4BVaVAiHo1Vdf1YkTJ3TLLbe4jgJElOjoaP3lL3/R5MmTdeLECddxgIDh\njBtQQw4ePKgOHTpo8eLF6tKli+s4QMSx1io5OVnJycm65557XMcBfoahUoobQsz48eN17rnn6pln\nnnEdBYhY27dvV69evbR69WrFxcW5jgOcUtXiVrsmwgCR7t1339Unn3yizz77zHUUIKJdfPHFeuih\nh/S73/1Oy5cvV1QUM4TgbfwvGAiwgwcP6tZbb9U//vEPNWrUyHUcIOLdeeedksQF6BEWGCoFAuz6\n669Xs2bNGCIFQkhubq5+/etfcxF6hAyGSoEQMHfuXK1bt04bN250HQVAKXFxcZo2bZpuuukmrVy5\nUrVq1XIdCagSzrgBAbJr1y716NFDH3zwgbp16+Y6DoDTFBUVadCgQerTp4+mTZvmOg4iHKtKKW5w\n6OTJk+rXr59GjBihe++913UcAOXYt2+funbtqjfeeEOJiYmu4yCCsQEv4NCjjz6qs846i72igBDX\nokULvfbaaxo/frwOHDjgOg7gN864AdWUkZGhm266SevXr1eLFi1cxwFQCQ899JDWrVunxYsXs0UI\nnOCMG+DArl27NH78eL311luUNsBDHnnkEf3www967LHHXEcB/MKqUqCKjh07ptTUVE2dOlVXXnml\n6zgA/FC7dm3NmTNHPXv21OWXX65hw4a5jgRUCkOlQBVYa3XTTTfpp59+0ltvvSVj/D7bDSAErFmz\nRtdcc42ys7PVoUMH13EQQRgqBYLo8ccf1+eff66ZM2dS2gAP69mzp55++mkNGzZMBw8edB0HOCPO\nuAF+mjNnju69916tXr1a559/vus4AALgvvvu09q1a5Wenq569eq5joMIwD5uFDcEwccff6xrr71W\nmZmZ6tSpk+s4AAKksLBQY8aMUe3atfXPf/6TKyugxjFUCtSwTZs2acSIEZo1axalDQgztWrV0ptv\nvqkDBw7oD3/4gzgRgFBFcQMqITc3V4MHD9Zzzz2nlJQU13EA1ID69etrwYIFWrlypf70pz+5jgOU\nie1AgDP45ptvlJSUpGnTpmns2LGu4wCoQdHR0UpPT1efPn109tln66677nIdCfgZihtQgV27dikx\nMVF33nmnJk6c6DoOgCBo3ry5li1bpsTERJ04cUL33Xef60jAKRQ3oBzbt2/XwIEDdd999+n22293\nHQdAEF100UXKzs5WYmKiTp48qQcffNB1JEASxQ0o05dffqmkpCQ9+uij+v3vf+86DgAHWrZsqezs\nbA0YMEDHjh3To48+yr6NcI7tQIDTrFixQmPGjNFf/vIX3XDDDa7jAHBs//79uuaaa9SxY0e99NJL\nqlOnjutICANsBwIEwKxZszR69Gi9+eablDYAkqRmzZopOztb//nPf3T11VfryJEjriMhglHcAElF\nRUWaNm2a/vjHP+rDDz/UgAEDXEcCEEIaNWqkd955R3Fxcerbt6/y8vJcR0KEorgh4uXn5+vaa69V\nRkaGVq9erfj4eNeRAISg2rVr629/+5smTJigXr16adGiRa4jIQJR3BDR1q1bp65duyouLk7Z2dlq\n1qyZ60gAQpgxRnfccYfeeecdTZo0SX/84x9VWFjoOhYiCIsTEJGKior07LPP6rHHHtOMGTM0cuRI\n15EAeMz+/fs1btw4FRYW6rXXXlNsbKzrSPAQFicAlbRjxw71799faWlp+uSTTyhtAKqkWbNmyszM\n1LXXXqsePXrolVde4RqnqHEUN0SMwsJC/fWvf1XPnj01fPhwLV++XHFxca5jAfCwWrVq6Z577tGH\nH36o//u//9PVV1/NwgXUKIobIsLHH3+sbt26ad68efroo4909913q1atWq5jAQgTHTt21Jo1a3Tl\nlVeqR48e+tOf/qQff/zRdSyEIYobwtqePXt00003acyYMbrvvvuUnZ2tSy+91HUsAGGobt26uv/+\n+7Vhwwbl5OTosssu09tvv83wKQKK4oawtH//fv33f/+3OnfurPPPP19btmzRddddx+VqANS4iy66\nSPPnz9fzzz+v6dOnq0ePHsrKynIdC2GC4oawsmfPHt17772Kj4+XtVabN2/WY489prPPPtt1NAAR\nJikpSWvXrtW9996ryZMnq2/fvnr//fdVVFTkOho8jOKGsJCTk6Mbb7xRnTp10smTJ5WTk6Nnn31W\nzZs3dx0NQASLiorSmDFjtGXLFt1+++16+OGHdfnll+sf//iHjh496joePIh93OBZ3333nebOnau/\n//3v+vrrr3XnnXdq4sSJiomJcR0NAMpkrVVGRoaee+45rV69WjfccINuueUWdezY0XU0BFlV93Gj\nuMFTjh8/ruzsbL311ltasGCBrrrqKt18881KSUlR7dq1XccDgErbvXu3XnnlFb3yyiu64IILNG7c\nOI0ePVoXXnih62gIAoobxS1sff/998rKylJaWpref/99tW/fXqNGjdKNN96o8847z3U8AKiWkydP\nKisrS3PnztWCBQvUvn17paamKiUlRfHx8SyqClOeKW7GmGhJT0jKldRUUqa1tszlNn4eS3ELEydP\nngA108IAAAfOSURBVNTatWuVmZmppUuXasOGDerZs6dGjBihESNG6IILLnAdEQBqxPHjx7V06VIt\nXLhQ6enpKioqUnJyspKTk9W3b1+upxxGvFTc5krKsNbO9D3OlTTKWruxmsdS3Dzqm2++0Zo1a/Tp\np59qzZo1Wr9+vdq0aaOkpCQNHDhQffv2VaNGjVzHBICgstZq27ZtSk9P19KlS/Xxxx+radOmuuKK\nK9S7d2/16NFD8fHxqlevnuuoqAJPFDdjTIKktZLOsdYe8T03w/d4bFWP9X2N4hbCrLU6dOiQtm7d\nqi+++EKbN28+dTt58qTi4uI0ZMgQ9ezZU927d9c555zjOjIqKTs7W/3793cdA1XE3593FBUVacuW\nLVq1apVWrVql9evX66uvvtIll1yiTp06qVOnToqPj1dcXJxiY2MpdCGuqsUt2LO5u0o6XFLEfNZL\nml7NY+HY0aNHtX///lO33bt3a+fOndq5c6d27dqlnTt3SpLatWunjh07qkOHDho8eLA6dOigli1b\n6pFHHtH//M//OP4uUBX8w+9t/P15R1RUlDp06KAOHTrolltukSQ9/PDDSk1N1aZNm/TZZ5/pww8/\nVG5urr7++mudf/75iouLU9u2bdWyZUtdcMEFP7tFR0czf86DXBS3/NOey5dU1v4N/hyLajhx4oR+\n/PFHHTt27NTt+++/V0FBQbm3w4cPa//+/fr3v/+t/fv368SJE2rWrNmpW6tWrRQbG6s+ffooNjZW\nsbGxiomJ4YcEAARQ7dq1lZCQoISEhJ89f+LECX399dfKzc1Vbm6u9u7dq+zsbO3du/fU7eTJk/rV\nr36lpk2bqkmTJmratOnP7kdHR+uss84q89aoUSM1atSI1fwOhNV/8WeeeUaSTl0Xzp8/q/KaYL9H\nYWGhTp48eep2+uPTb6d//cSJEzp27NgvSpq1Vg0aNDh1q1+/vs466yxFR0crOjpaMTExp+5fcMEF\nio+PV0xMzM+KWuPGjSllABAi6tSpo7Zt26pt27ZKTk4u85jvvvtOBw8ePHXLz88/dX/Xrl06cuSI\nvv/++5/dfvjhh5/dr1WrlurVq6e6deuW+efpz9WqVevULSoqqtzHFd0v+bem9J9lPefPMVV5fUVq\n8t/DYM9xmyhpurW2aUXP+Xus72tMcAP+X3v3lxTFFcVx/Pcr8xwpKhtgzAIisoMhWUAwugBLV2Cw\n3EDQFQT/rABCFgDEvMfRFO+M6AKCQd89eeg71ARwMj0Y7xzm+6miKHqmpg4F1f3re2/fAwBII8Ma\nt55OT3Vek7R7zvdO9MsDAABk8ll7lUbEn5J2bX8/dLgraVWSbHfLqNp/vhcAAGDW1NjH7Us1m+r2\nJV2RtBMRv5bXHkhaGGz3Meq9AAAAs+bCtLwCALRXZjlCUu+szc0BfBq2L0fE0Xk/57NOlX5qti/b\nXrd91/aa7W7tmjA+2wu2e7Y/lO9Xa9eEydh+YHutdh0Yn+3F0o2mHxFPCW15lHPn4Nq3PlhihOlU\nloHtS7p94vhEGSb7diBPdKIllu0zW2Jhugz1ob0lyZKeSurZPu6UgRzKyea6pM3atWA8pTNNT1I3\nIn6vXQ9a25G0EhF7kmT70PZzrn3Tpyz56kmaP+PliTJM2hG3cuJZkbQxdHhX0v06FaGlZUmrEbFX\n/km7agLcUt2y0IbtBUkdNV1NkMeGpE1CW1odScPbYv0l6e9KtWCEiHhXpkf/1VDgPBkmbXDTx1ti\nLVeqBy1ExFZEvB76+UhlnU21ojCJOxHxRE3oRgJlSUJH0pztjXKXv167LrTyi6RN21fLNOnu8PkU\nKUycYTJPldIS6wKxvazm5MM0aRJlTdtPtetAa0tqbpLuRMSbsmzhwPbbiGDGIoGIuGF7W82FfnOw\nEwNSmTjDZB5xw8VyW9IPtYvAeMpd/kZEvK9dC1qbU3On/0Y6Hu1+rGadIhIoYbuvZp3w9bKVFmZE\n5hG3Fzp9oZ8X8/zp2P5R0hqjbanck7Rwoh+fba9GxKVKNWE8r3T6rr6vZvoUOfwm6VZE7NnuqZk2\n3R8sckcKE2eYzCNurVpiYTrZXlGzsTJPQyUSEV9HxKXBl6QtSQ8JbSm8lI6fdjt1HNOtjLYtDp4o\njYgtNWverlUtDG1NnGHSBjdaYuVX1rXNSzoq+xItMuQP/L8i4kDNhX54Pdu3kh7VqQhtlKntt7a/\nGTo8J57sTuU8GSbzVKnUDDM+tN1R0xLr3mDdBqZb2ftrW80i6WEE75xowZJIRNy0/XN5wORQ0h9M\ns6XSlXTf9nNJX0l6wd9vOpUR0huSFiTdtP0yIp6VlyfKMLS8AgAASCLtVCkAAMCsIbgBAAAkQXAD\nAABIguAGAACQBMENAAAgCYIbAABAEgQ3AACAJAhuAAAASRDcAAAAkiC4AQAAJEFwA4CPsN21vT9o\nBG17u/T3BIAq6FUKACOUoHZF0jaNvAHURnADgBFsdyXtSOpExOvK5QCYcUyVAsBoB5JC0mHtQgCA\n4AYAo61IeiVpqXYhAPBF7QIAYBrZnpPUlfRIUkfSd7YVEc/qVgZgljHiBgAnlHVtfUn9iHgn6bGa\nkTcAqIqHEwAAAJJgxA0AACAJghsAAEASBDcAAIAkCG4AAABJENwAAACSILgBAAAkQXADAABIguAG\nAACQBMENAAAgiX8AMZRUQ45n0MsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1060650b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Execute this cell\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=15, usetex=True)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# plot the distributions\n",
    "fig, ax = plt.subplots(figsize=(10, 7.5))\n",
    "dist = norm(5, 1)\n",
    "x = np.linspace(0, 10, 1000)\n",
    "plt.plot(x, dist.pdf(x), c='black',label=r'$\\mu=5,\\ \\sigma=1$')\n",
    "\n",
    "plt.xlim(0, 10)\n",
    "plt.ylim(0, 0.5)\n",
    "\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel(r'$p(x|\\mu=5,\\sigma=1)$')\n",
    "plt.title('Probability of $x$')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "norm?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The Likelihood Function\n",
    "\n",
    "If we want to know the total probability of our *entire* data set (as opposed to one measurement) then we must compute the *product* of all the individual probabilities:\n",
    "$$L \\equiv p(\\{x_i\\}|M(\\theta)) = \\prod_{i=1}^n p(x_i|M(\\theta)),$$\n",
    "where $M$ refers to the *model* and $\\theta$ refers collectively to the $k$ parameters of the model, which can generally be multi-dimensional. In words, this is ***the probability of the data given the model parameters***. Note we're assuming the individual measurements are independent of each other.\n",
    "\n",
    "If we consider $L$ as a function of the model parameters, we refer to it as the ***likelihood of the model parameters, given the observed data***. Note that while the components of $L$ may be normalized pdfs, their product is not.  Also the product can be very small, so we often take the log of $L$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can write this out as\n",
    "$$L = \\prod_{i=1}^n \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}{2\\sigma^2}\\right),$$\n",
    "and simplify to\n",
    "$$L = \\left( \\prod_{i=1}^n \\frac{1}{\\sigma\\sqrt{2\\pi}} \\right) \\exp\\left( -\\frac{1}{2} \\sum \\left[\\frac{-(x_i-\\mu)}{\\sigma} \\right]^2 \\right),$$\n",
    "\n",
    "where we have written the product of the exponentials as the exponential of the sum of the arguments, which will make things easier to deal with later.\n",
    "\n",
    "That is, we have done this: $$\\prod_{i=1}^n A_i \\exp(-B_i) = (A_iA_{i+1}\\ldots A_n) \\exp[-(B_i+B_{i+1}+\\ldots+B_n)]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**NOTE**\n",
    "\n",
    "The Likelihood function is the joint probability density of all measurements. For each measurment we have a probablity density. These probability densities have **units**. They are the reciprical of the units of the measurement (as the integrated probability does not have units but the interval over which we are integrating does). For our case our positions might be in arcsecs so the probability density is in units of 1/arcsec. This means that the probability density is the rate of change in the probability per arcsec. \n",
    "\n",
    "For the Likelihood function we have the probability density of *n* measurements so the units for the likelihood is (1/arcsec)$^n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The argument of the exponential is just\n",
    "\n",
    "$$\\exp \\left(-\\frac{\\chi^2}{2}\\right).$$\n",
    "\n",
    "That is, for our Gaussian distribution\n",
    "$$\\chi^2 = \\sum_{i=1}^n \\left ( \\frac{x_i-\\mu}{\\sigma}\\right)^2.$$\n",
    "\n",
    "So, maximizing the likelihood is the same as minimizing $\\chi^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Core Idea Behind Maximum Likelihood Estimators\n",
    "\n",
    "Let's say that we know that some data were drawn from a Gaussian distribution, but we don't know the $\\theta = (\\mu,\\sigma)$ values of that distribution (i.e., the parameters).\n",
    "\n",
    "Then Maximum Likelihood Estimation method tells us to think of the likelihood as a ***function of the unknown model parameters***, and ***find those that maximize the value of $L$***. Those will be our Maximum Likelihood Estimators for for the true values of the model.\n",
    "\n",
    "Simple as that!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MLE applied to a Homoscedastic Gaussian\n",
    "\n",
    "Let's take a look at our astrometry example, using a model where all the measurements have the same uncertainty, drawn from a normal distribution, $N(0, \\sigma)$.\n",
    "\n",
    "All uncertainties being the same is known as having **homoscedastic** uncertainties which just means \"uniform uncertainties\".  Later we will consider the case where the measurements can have different uncertainties ($\\sigma_i$) which is called **heteroscedastic**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We have an experiment with the set of measured positions $D=\\{x_i\\}$ in 1D with Gaussian uncertainties, and therefore:\n",
    "\n",
    "$$L \\equiv p(\\{x_i\\}|\\mu,\\sigma) = \\prod_{i=1}^N \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(\\frac{-(x_i-\\mu)^2}{2\\sigma^2}\\right).$$\n",
    "\n",
    "Note that that is $p(\\{x_i\\})$ not $p(x_i)$, that is the probability of the full data set, not just one measurement. If $\\sigma$ is both constant and *known*, then this is a one parameter model with $k=1$ and $\\theta_1=\\mu$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Note**\n",
    "\n",
    "For multidimensional distributions, as we will see later in the course, the uncertainties come in the form of a covariance matrix where there uncertainties can be correlated. This still assumes Gaussianity\n",
    "\n",
    "$$ p(D\\; | \\;  \\pmb \\theta\\;) =  \\prod_{k=1}^{n} \\; \\frac{1}{(2\\pi)^{d/2} \\; |\\Sigma|^{1/2}} exp \\bigg[ -\\frac{1}{2}(\\pmb x - \\pmb \\mu)^t \\Sigma^{-1}(\\pmb x - \\pmb \\mu) \\bigg]$$\n",
    "\n",
    "where $\\Sigma$ is the covariance matrix and $d$ is the dimensionality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For practical reasons, it's better to work with the natural logarithm of the likelihood.\n",
    "\n",
    "We define the *log-likelihood function* as ${\\rm lnL} = \\ln[L(\\theta)]$.  The maximum of this function happens at the same place as the maximum of $L$. Given all that, we have:\n",
    "\n",
    "$${\\rm lnL} = {\\rm constant} - \\sum_{i=1}^N \\frac{(x_i - \\mu)^2}{2\\sigma^2}.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We then determine the maximum in the same way that we always do.  It is the parameter set for which the derivative of ${\\rm lnL}$ is zero:\n",
    "\n",
    "$$\\frac{d\\;{\\rm lnL}(\\mu)}{d\\mu}\\Biggr\\rvert_{\\hat \\mu} \\equiv 0.$$\n",
    "\n",
    "That gives $$ \\sum_{i=1}^N \\frac{(x_i - \\hat \\mu)}{\\sigma^2} = 0.$$\n",
    "\n",
    "(note: we should also check that the $2^{\\rm nd}$ derivative is negative, to ensure this is the *maximum* of $L$)\n",
    "\n",
    "(also note: any constants in $\\ln L$ disappear when differentiated, so constant terms can typically be ignored.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Since $\\sigma = {\\rm constant}$, that says \n",
    "$$\\sum_{i=1}^N x_i = \\sum_{i=1}^N \\hat \\mu = N \\hat \\mu.$$\n",
    "\n",
    "Thus we find that\n",
    "$$\\hat \\mu = \\frac{1}{N}\\sum_{i=1}^N x_i,$$\n",
    "which is just the arithmetic mean of all the measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The Sample Mean is an ML Estimator\n",
    "\n",
    "The mean of observations drawn from a $N(\\mu, \\sigma=const)$ distribution is a maximum-likelihood estimator of the distribution's $\\mu$ parameter.\n",
    "\n",
    "We'd intuitively guess that (and we often do), but this derivation clarifies our choice: as an estimator of the real value of $\\mu$, we adopt the value $\\hat \\mu$ with which it's maximally likely for the measured data set to occur.\n",
    "\n",
    "It also exposes the ***assumptions*** behind this conclusion; namely homoscedasticity and gaussianity of uncertainties. For example, if our uncertainties were Cauchy-distributed the mean (or any higher moment) is not defined"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Properties of ML Estimators\n",
    "\n",
    "Assuming the data truly are drawn from the model, ML estimators have the following useful properties:\n",
    "\n",
    "* **They are consistent estimators**; that is, they can be proven to converge to the true parameter value as the number of data points increases.\n",
    "* **They are asymptotically normal estimators**. The distribution of the parameter estimate, as the number of data points increases to infinity, approaches a normal distribution, centered at the MLE, with a certain spread. This spread can often be easily calculated and used as a confidence band around the estimate, as discussed below (see eq. 4.7).\n",
    "* **They asymptotically achieve the theoretical minimum possible variance, called the Cramér–Rao bound**. In other words, they achieve the best possible uncertainty given the data at hand; that is, no other estimator can do better in terms of efficiently using each data point to reduce the total error of the estimate (see eq. 3.33 in the textbook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quantifying Estimate Uncertainty\n",
    "\n",
    "Let's simulate 10 measurements ${x_i}$ drawn from a Gaussian with $\\mu = 5$ and $\\sigma = 1$, and compute the ML estimator for $\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11cc607f0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAowAAAHjCAYAAACpeU7EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VdWh/vF3ISCoEIKCdWQI4oCVIaEik5EhYVCUKcHr\neBUEr61X5Seot/dWqq041KkOICAVrkqYBYGQBAigoEJAHEAUiIrgnAQQWxmyfn/khBtDcnJOzrDO\n8P08D0/JOXvv8wba9GWtvdcy1loBAAAA1anjOgAAAAAiG4URAAAAXlEYAQAA4BWFEQAAAF5RGAEA\nAOAVhREAAABeURgBoBJjTCtjTKtanptgjOkYimv7kSHknwEgvlAYAcQsY8xQY0yRMabUGPOiMaZl\nNcfNNsY8Un6OpJ2Shtbi83pLKpA0sYZrV1soA2WMGRXoZxhjRnl+7TDGjAxeOgDRqq7rAAAQKtba\necaYvpJGWWtv93LoJElFFc7ZVMvPW2GMmatfl7XK1y6pzbX9yDDFGPNobc/3jI6Os9aeZ4zZIKlp\n8NIBiFaMMAKIdTUWNGvtSmvt+6H48FBeO0RS5Pkzs9a+b61dWd2BnhHI/1fh63HGmB1hyAggzBhh\nBBD3PFPFTa21U6p47zZJj6hslHK+57WOkjIltZbURNJoa22hH9duaoyZLamPpDxrbYbn2FaSRkv6\nUVJnSRustY/7+N54STsknSYpoYbvt8prebIOl9TaUwQ3eSuMnu9/V4WvU1Q2HQ4gxlAYAcQ1T/m7\nX2VFZ0ql9xJUNr3c0lp7oMJrj1pr0zxf56hs2jndj2v3tdZmeIrbTmNML5Xd+1ggqUWFzyoyxljP\nuVW+Z619QlKupI7W2gOea47z8v0meLuWMaa1pETPdb39ubWSZCVVnL5v7ckCIMZQGAHENWvtZmPM\nRkmJld76naTflY/+VZAhyXpG4IzKymBrP6+d5Xm/0HNPYxPPdX8sL3EeL0kaI2lfNe+NNsYUSiou\nf6/CNatT3eeMluS1JFbSyfN5n1d67a9+XANAlKAwAkDVmkjqbYzpUOkexGRJu2oagfNTU0lJVby+\nU1IrL++1Vtk0cJEfn+XtWv741XR0hWV8dlV9OIBoxkMvAOKWZxmc6uRIekzSXGNM4wqvl6js3sOK\n10modExtVFfaNqns3sTq3iu/DzEYn+OPzvr1A0XJkmyUPeADwEcURgCxrvJ0sCTJGNNJ0rAKLzWp\nfIy19n6VlaK5FV7OkpRkjJnkWSC7t6SJ1tr9XjIcd+0qPmuKpJJK6x72lfRXa+3U6t6TNE9SQvl7\nnnsQm6j6afJqP8fz+9N8yauy6efWxpiOnodl7it/w/M1gBjClDSAmOUpLr09v1+usunSU1VWpjpK\nGu4pfH0kJXoePkn0vPejZ03FyZImGWOyJI333Jd4m8qeSh6lsjI5yvMZHVVWQsuvZaq4doLK7j3M\nU1lRS/DkmK2yUbrJxpgklY0czrLWLvB8O9W+Z4wZLulRY8xolRXaAkmnGmMaV1Nkq7yWJ39vlRXB\nkZ6iWp0SlT30MtvzZzFS0gpjTJa1NtPLeQCikLHWus4AAIginiLep4bF0AHEEEYYAQD+ylTZ6CKA\nOME9jAAAf1l5lgYCEB+YkgYAAIBXTEkHyLMLAwAAQFSw1hp/z6EwBgGjtNHrwQcf1IMPPug6Bmoh\n1v/ujJFi+UdLrP/9xTr+/qKXMX53RUncw+iTCjsYAAAAxJ2wjzB6Nr5/VGU7F5wqKc9au8LfY2u6\nToDnzpFUvvCsVdmaZQAAAHHJxZT0FEk55QvCGmN2GGOGVbOdlLdja7pOrc71lMkfVLZwr5FUYq39\nPDR/FHAtNTXVdQTUEn930Y2/v+jG31/8CetT0p6tuDZISizffcAYM8nzdaavx0qaqLKdDJpUdZ0A\nz50sqVjSZGttoQ/fk+UeRgDBFuv3MAJwwxhTq4dewn0PY7LKRuwqblVVoLKts/w5NkVSsZfrpARw\nbiuVbfO10xizkfsXAQBAvAv3lHSypKJKrxWp6o3uqzs20YfrdKrtudbaNEny7Ps6WVKOpPOq+4YA\nAPGjZcuW+uKLL1zHALxq0aKFPv/886BeM1qX1Qlkosanc621K40xKZKKjDG9rLUrqzu24tICqamp\n3NsBADHqiy++YCk1RLyKS+fk5+crPz8/4GuGuzAWSBpe6bWmkkr8OLbY816Gl+sEcu4x1tp9xpgV\nqnoE9BjWogIAAJGo8kDWhAkTanWdcN/DuFHHl69kSXl+HltQw3UCObeyppI2VfMeAABAzAtrYbTW\nbpaUZ4wZUuHl3pLGSZIxprcxZlRNx3p5b3wg5xpjWhljytdflDGmo6QNLKsDAADiWViX1ZEkY0xj\nlS2avVNSkqRca+18z3sTJbUqX2KnhmOrfa+253oK4hyV3ec4WdKuites5vthWR0AQceyOpHJsySJ\n6xiAV97+e1rbZXXCXhhjDYURQChQGCMThfH/7Nu3T4888ojuv/9+JSQkuI6DCkJRGNlLGgAA+G34\n8OEqLCzUbbfdFtbPLSyscU+NWh0bDfbt2+fssymMAADAL1OmTNF9992nrKwsZWZmaurUqSH7rOHD\nh6tOnTqqU6eOTjjhBK+lydux+/bt05gxY/TEE0/o/vvv14oVK0KWOdhWrFihNm3aaMqUKc4yROs6\njAAAwJFRo0Yd+/2QIUO8HBmYffv26bTTTtPmzZtlrVWTJk3UsmXLWh07atQopaWlaeTIkZKkNm3a\naO7cuerQoUPI8gfD/v37lZKSoqKiynuOhBcjjAAAwGdt2rTRE088cezrxx57TG3atAnJZ40fP14J\nCQlq3LixOnToUG1ZrOnYTZs2ad68ecrI+L9lmPv06aNHHnmk1tny8vKUlpampk2bHhvVLB/ZDKbG\njRsrISFBTZs2Dep1/UVhBAAAPtu1a5dat2597OuNGzcqKSkpZJ81ZcoUJSUlKSUlxes9id6OLSgo\nUJMmTdS4ceNjryUnJysvr7olmL0bPXq00tPT1aZNG02dOlVJSUl66aWXNHfuXM2ZM6dW14x0TEkD\nABACFbdnixSBPuFdWFgoY4w6dep07LVdu3apb9++1Z4zZswYFRcXH/f55U/yGmPUt2/fY1PFFeXk\n5EiSVq5cqdGjRystLU2fffZZlZ/j7diCgoLjRuiaNm2qkpKqNprzbu7cuZoyZYpWrFihK664QpL0\n448/auLEidqxY0eV5wTyZxApKIwAAIRALC6/s2lT2cZnlad7H3jggWrPmTRpUsCf26tXL23cuFFN\nmzbVypUr1atXr6AcWxsTJ05UUlLSsbIoSRkZGRo9erSmTp1aZekLxp+Ba0xJAwAAn1Seji6f9q34\nWqgkJCSod+/ePo0KVj42OTn5uIdGioqK1KRJ5Z2Ca1ZSUlLt97tz506/rxctGGEEAAA+2bBhw69K\nVkFBgYwxXp80rjgdW5m/07FFRUW/mg739diUlJTjimZBQYH69Onj07Uqat269XHlMzc3V8YYjRgx\nospzgvln4AqFEQAA+GTTpk0qLi7W5s2btWvXLk2cOPHYe/PmzdPQoUOPO6e207GFhYXatGnTsWtu\n3rxZnTt3PjYdvmLFCu3atUujRo2q8diOHTuqT58+mj9//rFlgFasWHHsoZeK16rJ+PHjlZaWdmy6\nu6SkRPfdd5/Gjx+v9u3bV3lOLExJszVggNgaEEAosDVgZIr3rQFTUlKOLYb96KOPqnXr1urdu7f6\n9OmjrKysoH7W5s2bNXz4cBljNHr0aLVu3fpXaz7ed999KiwsVFZWVo3HSmXrGY4fP15JSUnauXOn\n+vbte+yY++67T/Pmzav2gZrKpk6dqsmTJx97kGbMmDEaPHhwkL7zX9u3b59mz56tMWPGqFOnTnr0\n0UdrvC+TvaQjEIURQChQGCNTPBfGefPmKS8vTy+++KLrKCFRcfQx2oWiMDIlDQAAapSVlRWRSwUh\nPHhKGgAA1MgYo8zMTNcxQqKwsNDnh2niFVPSAWJKGkAoMCUdmeJ5ShrRIxRT0owwAgAAwCsKIwAA\nALyiMAIAAMArCiMAAAC8YlkdAAB81KJFC5aWQcRr0aJF0K/JU9IB4ilpAKHAU9IAQoGnpAEAABAS\nFEYAAAB4RWEEAACAVxRGAAAAeEVhBAAAgFcURgAAAHhFYQQAAIBXFEYAAAB4RWEEAACAVxRGAAAA\neEVhBAAAgFcURgAAAHhFYQQAAIBXFEYAAAB4RWEEAACAVxRGAAAAeEVhBAAAgFcURgAAAHhFYQQA\nAIBXFEYAAAB4RWEEAACAVxRGAAAAeEVhBAAAgFcURgAAAHhFYQQAAIBXFEYAAAB4RWEEAACAVxRG\nAAAAeEVhBAAAgFcURgAAAHhFYQQAAIBXFEYAAAB4RWEEAACAVxRGAAAAeEVhBAAAgFcURgAAAHhF\nYQQAAIBXFEYAAAB4RWEEAACAVxRGAAAAeEVhBAAAgFcURgAAAHhFYQQAAIBXFEYAAAB4RWEEAACA\nVxRGAAAAeEVhBAAAgFcURgAAAHhV13UAAIAbJSUl2rZtm77++mvVq1dPLVq00EUXXaS6dfm/BgC/\nxk8FAIgjhw8f1quvvqp//OMf2rRpky688EKdeeaZOnz4sHbs2KG9e/dq8ODBGjNmjC677DLXcQFE\niLAXRmNMgqRHJe2QdKqkPGvtCn+Prek6gZxb4RoTJVlr7f1B+NYBwKns7Gz94Q9/0DnnnKO7775b\n6enpatCgwa+O+e677/Tqq69qxIgRSk5O1pNPPqmWLVu6CQwgYhhrbXg/0JjZknKstVM9X++QNMxa\n+74/x9Z0nUDO9bzWW9JkSXO8FUZjjA33nyGA2GeMFKwfLb/88ovuvvtuLV26VJMnT1Z6enqN5/zz\nn//U3/72Nz377LN64YUXNGzYsOCEAeCUMUbWWuP3eeEsO8aYTpI2SEq01u73vDbJ83Wmr8dKmiip\nQFKTqq4TyLmer1tJ6uP5tYvCCCDcglUYi4qKNHjwYJ166qmaPn26EhIS/Dp/w4YNysjI0C233KI/\n/vGPMsbv/58BEEFqWxjD/ZR0sqSS8qLmUaCyYubPsSmSir1cJyWAcyXpNmvtFEn8ZAQQtX744Qdd\nfvnlSk5O1ty5c/0ui5LUuXNnrVu3TvPmzdPdd98t/oEMxCcXhbGo0mtFkpr4cWyiD9fpVNtzPfct\n/tXbNwEAka6kpETp6ekaOHCg/va3v6lOndr/uD/jjDOUn5+vtWvX6k9/+lMQUwKIFtH6lHQg/8St\n9lxjzChJWdbaA/5c8MEHHzz2+9TUVKWmptY2GwAE7MiRI8rIyFCXLl30yCOPBGUauUmTJlq2bJl6\n9OihM844Q7fffnsQkgIItfz8fOXn5wd8nXAXxgJJwyu91lRSiR/HFnvey/Byndqcu0/SOEmtK/1w\nNcaYcdbaE6rIKOnXhREAXBs7dqzq1KmjZ555Jqj3HDZv3lxLly5V165ddfHFF6tHjx5BuzaA0Kg8\nkDVhwoRaXSfcU9Ibdfz0c7KkPD+PLajhOrU5N9dae5619oTyX5LmSXrUW1kEgEgyZ84cLVmyRLNm\nzQrJAtxJSUmaMWOGMjMz9c033wT9+gAiU1gLo7V2s6Q8Y8yQCi/3VtnInowxvT3Twl6P9fLe+ADO\nHReUbxIAHNm9e7d+//vf6/XXX1eTJlXdGh4c6enpuvXWW3XrrbfyEAwQJ1ysw9hYZYtm75SUpLKR\nvfme9yZKalVheRtvx1b7XqDnVrhGllhWB4AD/i6rc/ToUfXp00dpaWm6//7Q7zVw+PBhdenSRaNH\nj9Ztt90W8s8DEBxRsQ5jLKIwAggFfwvj3//+d82ZM0erVq3SCSeE5y6arVu3qmfPnnr33XeVlJQU\nls8EEBgKoyMURgCh4E9h3Lt3r9q3b6+1a9fqggsuCG2wSh5//HGtWLFCy5YtY1FvIApEy8LdAIAg\nu+uuuzRmzJiwl8Xyz969e7cWLFgQ9s8GED6MMAaIEUYAoeDrCGN2drbuuOMOffTRR2rYsGHog1Uh\nPz9fN910k7Zu3aqTTz7ZSQYAvmGEEQDizJEjRzR27Fg9/fTTzsqiVLbOW/fu3fXXv7JJFhCrGGEM\nECOMAELBlxHGKVOm6LXXXtPKlSud3z+4e/dudejQQR9++KHOPPNMp1kAVI+HXhyhMAIIhZoK48GD\nB9W2bVstXLhQnTt3Dl8wL8aNG6d9+/Zp8uTJrqMAqAaF0REKI4BQqKkwPvTQQ9q6datef/318IWq\nQVFRkc4//3y99dZbOv/8813HAVAFCqMjFEYAoeCtMJaUlKhNmzYRuf7hxIkTVVBQoDlz5riOAqAK\nFEZHKIwAQsFbYfzzn/+swsJCTZ8+PbyhfPDzzz+rdevWWrFihdq1a+c6DoBKKIyOUBgBhEJ1hXH/\n/v1KSkrSunXrdN5554U/mA8effRRffDBB3r11VddRwFQCYXREQojgFCorjD+5S9/0SeffKKZM2eG\nP5SPoqHUAvGKwugIhRFAKFRVGA8ePKhWrVppzZo1TnZ18ceECRP05Zdfatq0aa6jAKiAwugIhRFA\nKFRVGF944QXl5uZGxTZ8xcXFatOmjbZs2aKzzz7bdRwAHhRGRyiMAEKhcmEsLS3V+eefr5dfflk9\nevRwF8wPd911lxo0aKCJEye6jgLAg60BASCGLV68WImJierevbvrKD77wx/+oGnTpungwYOuowAI\nEIURAKLAk08+qXvuucf5FoD+SEpKUvfu3TVjxgzXUQAEiCnpADElDSAUKk5Jb9y4UUOGDNHOnTtV\nr149t8H8tHr1at12223atm2b6tRhjAJwjSlpAIhRTz31lO68886oK4uS1LNnT5188snKzs52HQVA\nABhhDBAjjABCoXyE8bvvvlPbtm1VWFioxMRE17FqZcaMGZo5c6Zyc3NdRwHiHiOMABCDpk+friFD\nhkRtWZSkzMxMbdmyRTt27HAdBUAtURgBIEKVlpZq8uTJGjNmjOsoATnxxBN100036aWXXnIdBUAt\nURgBIELl5uaqSZMm6ty5s+soAbvtttv0j3/8Q7/88ovrKABqgcIIABFq0qRJuv3226NqKZ3qnHfe\nefrtb3+rhQsXuo4CoBZ46CVAPPQCIBSMkRITm+rLL7/UKaec4jpOUGRlZWny5MlauXKl6yhA3OKh\nFwCIMSNGjIiZsihJgwcP1scff6xPP/3UdRQAfqIwAkCEKS0tlSTdeuutjpMEV/369XXzzTdrypQp\nrqMA8BNT0gFiShpAsK1Zs0aXX95TpaU2Ju5frGj79u26/PLL9dVXX6lu3bqu4wBxhylpAIgRr7zy\niiTFXFmUpPPPP1+tWrXS8uXLXUcB4AcKIwBEkIMHD2r+/PmuY4TUTTfddKwUA4gOFEYAiCALFizQ\nZZdd5jpGSGVmZionJ0fFxcWuowDwEYURACLIK6+8optvvtl1jJBKTExUWlqaZs+e7ToKAB9RGAEg\nQuzevVubNm3SoEGDXEcJOaalgehCYQSACDFz5kwNHz5cDRo0cB0l5NLT01VYWMiajECUoDACQASw\n1mrGjBm66aabXEcJi7p16+raa6/V//7v/7qOAsAHFEYAiACbN2/W4cOH1aVLF9dRwubaa6/VrFmz\nxFq2QOSjMAJABJg1a5ZGjBgRk2svViclJUWlpaXatGmT6ygAakBhBADHSktLlZWVpREjRriOElbG\nGI0YMUKzZs1yHQVADSiMAODYO++8o1NOOUUXX3yx6yhhN2LECGVlZR3bPxtAZKIwAoBj8TgdXe7i\niy9WQkKC1q1b5zoKAC8ojADg0NGjRzVnzhxlZma6juIM09JA5KMwAoBDq1ev1plnnqm2bdu6juLM\niBEjNGfOHB05csR1FADVoDACgEOvv/66rr32WtcxnEpKSlKLFi2Un5/vOgqAalAYAcCRQ4cOaf78\n+crIyHAdxbkRI0bo9ddfdx0DQDUojADgSF5eni644AKde+65rqM4l5GRoTfeeEOHDx92HQVAFSiM\nAODI/PnzNWzYMNcxIsLZZ5+tNm3aaPXq1a6jAKgChREAHDhy5IjeeOMNDR482HWUiDFkyBDNmzfP\ndQwAVaAwAoADa9euVYsWLdSyZUvXUSLGkCFDtHDhQh09etR1FACVUBgBwIH58+dryJAhrmNElDZt\n2qh58+Zav3696ygAKqEwAkCYlZaWasGCBRTGKgwZMkTz5893HQNAJRRGAAiz9957TwkJCbrgggtc\nR4k4Q4cO1fz582WtdR0FQAUURgAIM6ajq9euXTvVr19fmzZtch0FQAUURgAII2ut5s2bp6FDh7qO\nEpGMMRo6dChPSwMRhsIIAGH0wQcfyFqr9u3bu44SscqX12FaGogcFEYACKPy6WhjjOsoESslJUX/\n/Oc/tXXrVtdRAHhQGAEgjJiOrpkxhkW8gQhDYQSAMPnss89UVFSkSy+91HWUiDdkyBAtWLDAdQwA\nHhRGAAiTRYsWadCgQapThx+9Nenatat2796t3bt3u44CQBRGAAibRYsW6aqrrnIdIyrUrVtXAwYM\n0OLFi11HASAKIwCExY8//qj3339fvXr1ch0lagwaNIjCCEQICiMAhMGyZct0xRVXqGHDhq6jRI30\n9HS9/fbbOnDggOsoQNyjMAJAGCxevFiDBg1yHSOqNGrUSJdddplycnJcRwHiHoURAELs0KFDysnJ\n0cCBA11HiTqDBg3SokWLXMcA4h6FEQBCbM2aNbrgggt0+umnu44Sda666iotWbJER44ccR0FiGsU\nRgAIMZ6Orr1zzz1X55xzjtavX+86ChDXKIwAEELWWi1evJjCGACmpQH3/C6MxpiWxphenl8djDGN\nQxEMAGLBRx99JEm6+OKLHSeJXhRGwL26vhxkjBkqKVNSK0m7JBV73moqqYkxxkraJOkRa+3+UAQF\ngGhUPrpojHEdJWp16tRJP/30k7Zv367zzz/fdRwgLnktjMaYBEn3S9phrc2o4diOkh4wxrxnrZ0f\nxIwAELUWLVqkhx9+2HWMqGaMOTbKeO+997qOA8SlaqekjTGtJI2y1t5nrZ1a04WstZuttfdJKjTG\njAxmSACIRt988422b9+unj17uo4S9a688kotWbLEdQwgbhlrresMUc0YY/kzBFCVadOmKScnR1lZ\nWX6fa4zEj5b/8/PPP+s3v/mNdu/erYSEBNdxgKhljJG11u97ZHy6hzGYPNPcj0raIelUSXnW2hX+\nHlvTdWp7rmdkdbKkFElFkoZbazcH888AQHx48803NWTIENcxYsJJJ52k7t27KycnR8OHD3cdB4g7\nQRlhNMb0stau9PHY2ZJyyqe5jTE7JA2z1r7vz7E1Xae25xpjXpQ0TpKRNFdSK2vteV6+H0YYARzn\n8OHDatasmT799FM1b97c7/MZYTzec889p4KCAk2fPt11FCBq1XaE0Wth9IzEjfLhOpnW2s41fpgx\nnSRtkJRY/jS1MWaS5+tMX4+VNFFSgaQmVV2nluc2lTSy4lPenqfD7/P2vVEYAVRlzZo1uueee7Rx\n48ZanU9hPN6uXbvUtWtX7d27V3XqsIwwUBshmZK21u4zxoxRWZkq8nJoKx8/L1lSSaWldwpUVuL8\nOTZPUrGX66TU5twqlgTqLN8KMwD8SnZ2ttLT013HiCmtW7dWkyZNtHnzZiUnJ7uOA8QVX/6JNlpl\n07JtqvslaYyPn5es44tnkaQmfhyb6MN1OgVwrowxvT3T1PeqrHwCgF+ys7PVr18/1zFizsCBA3la\nGnCgxsLoeRikpomR2cGJ47NAJmpqPNfzPfdR2T2Mk40xHQL4PABx5ttvv9WuXbvUpUsX11FizoAB\nA7R06VLXMYC449NT0jXt3mKt3efj5xVIqvx4W1NJJX4cW+x5r/JC4hWvE8i5kiRr7eeSMj0jjSmS\njnsop9yDDz547PepqalKTU2t7lAAcSAnJ0e9evVSvXr1XEeJOd27d9e2bdv0/fffq1mzZq7jABEv\nPz9f+fn5AV8n3MvqbNTx08/JKruv0J9jC2q4TiDnVrbL86taFQsjADAdHTonnniievXqpeXLl+v6\n6693HQeIeJUHsiZMmFCr64T1MTPPeoZ5xpiKC5P1VtkyNuX3Do6q6Vgv740P4NxxxpgEY0zv8hc9\nT4kX+7pkEACUlpYqNzeXB15CaODAgUxLA2EW9p1ejDGNVbZo9k5JSZJyy/eeNsZMVNkDNpk+HFvt\ne7U917Nod67K7nN8VFJRTftis6wOgIo2btyoG2+8UVu3bg3oOiyrU709e/bokksu0bfffqu6dcO+\n/wQQ1UKyDiNqRmEEUNHDDz+sH3/8UU899VRA16EwetehQwc9//zz6tatm+soQFSpbWH0e0raGDOy\nqt8DAKTly5dz/2IY8LQ0EF61uYexbzW/B4C4VlJSovfff189e/Z0HSXmcR8jEF61KYymmt8DQFxb\nsWKFunXrpoYNG7qOEvMuvfRSffnll9qzZ4/rKEBcqE1htNX8HgDiGtPR4VO3bl2lpaVp2bJlrqMA\ncSHQEUYAgCRrLesvhhn3MQLhE9Z1GAEgVm3btk3GGJ1//vmuo8SNtLQ0rVy5UocPH3YdBYh5FEYA\nCILy0UVjmIQJl9NPP12tW7fWO++84zoKEPMCvYcRACDuX3SlX79+Wr58uesYQMxjhBEAAvTzzz9r\n3bp16tWrl+socSc9PV3Z2dmuYwAxrzaFsbia3wNAXFq9erU6duyohIQE11HizmWXXaYdO3bou+++\ncx0FiGl+F0Zr7Ziqfg8A8YrpaHfq16+v1NRU5ebmuo4CxDSmpAEgQNnZ2UpPT3cdI25xHyMQesay\nu31AjDGWP0Mgfn3++ef63e9+p2+++UZ16gTv3+DGSPxo8U1hYaG6dOmir7/+Oqh/B0AsMsbIWuv3\ncg78LwsAArB8+XKlp6dTVBxq1aqVmjRpoi1btriOAsQsfsIBQACYjo4MPC0NhJbPhdEYc9QYMySU\nYQAgmhw+fFirVq1SWlqa6yhxj/sYgdDyZ4Tx2Hy3MWakMeYRY0z7EGQCgKiwfv16tWnTRs2bN3cd\nJe5dfvmdvJ2EAAAgAElEQVTlKigo0P79+11HAWKSv1PSTY0xOyS9JGm8pE3GmL8GPxYARL7y7QDh\n3sknn6wuXbpo1apVrqMAMcnfwjhe0qOSEj2/MiVlGGMGBzsYAEQ67l+MLNzHCISOz8vqGGOKJI20\n1s6v9HoTSVnW2rj8qcmyOkB8+vbbb3X++efr+++/V7169YJ+fZbV8d9HH32kq666Srt27ZIxfq8a\nAsSFcCyrs1FSn8ovWmtL/P1QAIh2OTk56t27d0jKImqnXbt2Onz4sD777DPXUYCY409hfEzSaGPM\nFRVfNMb0ktQ0qKkAIMIxHR15jDFKT0/naWkgBHwujNbaPElTJa0wxnxmjFlujNkgKVfS5FAFBIBI\nU1paqpycHApjBOI+RiA0/HroxVo7WlKGpH2SkiQVS0qz1k4NQTYAiEibNm1Ss2bN1KJFC9dRUEmf\nPn20du1a/etf/3IdBYgpfu/0Yq2da61Nsda2sdamWWtXhCIYAEQqpqMjV9OmTdWuXTu99dZbrqMA\nMYWtAQHAT6y/GNnY9QUIPgojAPihpKREW7ZsUc+ePV1HQTW4jxEIPgojAPhhxYoV6tatmxo2bOg6\nCqrRuXNn7d27V3v27HEdBYgZARdGY8woz1PTQzy/WgYeCwAiE9PRke+EE05Qnz59mJYGgigYI4y7\nJM2z1s737ALD8voAYpK1VsuXL6cwRgHWYwSCy+etAVE1tgYE4sfWrVvVv39/ff755yHfeo6tAQOz\nZ88e/fa3v9X333+vE044wXUcIGKEY2tAAIhr5dPR7FMc+c466yydddZZ2rBhg+soQEwIxj2MQ4MR\nBAAiHdPR0YVpaSB4fC6Mxpgdnu0AH/E83NJBkqy184wxWaGLCADu/fzzz1q3bp169erlOgp8RGEE\ngsefEcaXVLYlYF9JcyVtMsYcNcYcldQnFOEAIFKsXr1aHTt2VEJCguso8FGPHj304Ycfqri42HUU\nIOr5XBittY9ZazM82wLWUdle0pmSNkvqHaqAABAJmI6OPg0aNFCPHj2Ul5fnOgoQ9Wp9D6O1ttBa\nO1fSI5JaBy8SAEQe1l+MTkxLA8Hhzz2M/6+ahbl3SXo0mKEAIJIUFhaquLhYHTp0cB0FfurXr5+y\ns7PF8mdAYOr6cewYlY0kWs+SErsklXheKwp+NACIDMuXL1daWprq1GElsmjTtm1b1a1bV1u3blW7\ndu1cxwGilj/3MLbx3LvYRlKGyh6CKZRULGmfMeZFY0yWMWZIaKICgBvcvxi9jDFMSwNBELSdXjzr\nMWZKstbazKBcNAqw0wsQ2w4dOqRmzZrps88+U/PmzcP2uez0Ejzz58/X5MmTKY2Aar/TC1sDBojC\nCMS21atXa+zYsdq4cWNYP5fCGDz79u3T2WefrW+//VYnnXSS6ziAU2wNCAAhwHR09EtISFCHDh20\nZs0a11GAqOX1oRfPNHNiLa5bZK2dX7tIABA5srOz9cwzz7iOgQCV38dI+Qdqp6anpHeFJQUARKBv\nv/1Wu3btUpcuXVxHQYDS09N10003uY4BRC2vhdFauzlcQQAg0uTk5Kh3796qV6+e6ygIUHJysr7/\n/nt9+eWXOvfcc13HAaJOre9hNMY0DmYQAIg02dnZSk9Pdx0DQVCnTh317duXJ6WBWgrkoZcpQUsB\nABGmtLRUOTk5FMYYwnqMQO3VuKyOMabyXlhNJY2WtNHz+1xV2OnFWvt+kDNGNJbVAWLTxo0bdeON\nN2rr1q1OPp9ldYLv66+/1kUXXaTvv/9edev6s9EZEDtqu6xOtf+LMcYkSCpQ2U4ulS88yVo71Rgz\nStJjlc5rJamVtXa/v2EAIFIwHR17zjjjDLVo0ULvvvuuunXr5joOEFW8/RNruKTh3h58sdZOUaWp\naWNME0m3SXoiKAkBwIHs7Gz993//t+sYCLLyaWkKI+Afb/cwFtfmKWlrbYnK9pgGgKhUUlKiLVu2\nqGfPnq6jIMi4jxGonWoLo7V2Xm0vGsi5AODaihUr1L17dzVs2NB1FARZt27dtG3bNv3www+uowBR\nxa+npI0xj4QqCABECu5fjF0nnniiLr/8cuXl5bmOAkQVf5fVSa78gjGmlzGmV5DyAIBT1lq2kItx\nTEsD/vOrMFpr04wxI40xLY0xjT1FMU3S5NDEA4Dw2rZtm+rUqaPzzz/fdRSESHlhZEk0wHf+TkkP\nUVlB3KWy5XYmS0qQlBH8aAAQfuXT0cb4vUwZokSbNm3UsGFDffjhh66jAFHD3ynpxyTtlNRX0nmS\nXpI0nj2nAcSK7OxspqNjnDGGaWnATzXu9PKrg43pba1dUem1oZJ2xtsOL+XY6QWIHT///LNOP/10\nffXVV0pISHCahZ1eQuuNN97Q3//+dx5+Qdyp7U4v/t7DuKKKl1d4fgFAVFu9erU6duzovCwi9K64\n4gq9++67OnjwoOsoQFSotjAaYxr7cgHPQt2/uofR13MBIJIwHR0/GjdurOTkZOXn57uOAkQFbyOM\nt/l6kSpGHn0+FwAiBesvxhfuYwR8520v6c3GmBclzZFUVOH1Emvt5+VfGGM6VDovU9KGoCUEgDDY\ntWuXSkpK1LFjR9dRECbp6em69tprXccAokK1hdFau8IY01plT0ZX1MQYUyypt6RNKltep+LNk7nW\n2vlBTwoAIVQ+HV2njr+LRyBadejQQcXFxSosLFSrVq1cxwEimrcRRllrp0iaUvl1Y0wnlT3o0tda\nWxiibAAQNtnZ2Yw2xZk6deocm5YeM2aM6zhARKvVP6WttZskvURZBBALfvnlF+Xn5ystLc11FIQZ\n9zECvqn13Itn9BEAot5bb72liy66SKeeeqrrKAiztLQ0rVq1SocPH3YdBYho3KwDIO4tW7ZM/fv3\ndx0DDjRv3lxJSUlav3696yhARPO2DmMrz97RfjHGNDbGjAwsFgCED+svxjempYGaVVsYPfcn7jPG\nPFLF0jlVMsaMknS/tXZqsAICQCjt3r1b33zzjVJSUlxHgSMURqBmPu0lbYyZKKmjpF2SdkoqUdna\njE0lJXl+WUmTrbUrQ5Y2ArGXNBDdpkyZolWrVum1115zHeVX2Es6fA4dOqRmzZrps88+U/PmzV3H\nAUKqtntJe11Wp5y19j7Ph7SS1EdSoudXiaQ8a+39fgRNkPSopB2STvWcX+Ve1N6Orek6tT3X8z3O\nkdRJZetMjrLWbvb1+wMQXbKzs3XNNde4jgGH6tevryuuuEK5ubm67rrrXMcBIpJPI4zVnmxMR3/L\nlDFmtqSc8mlrY8wOScOste/7c2xN16nNuZIKVbbu5F9Uthj5VJWNrCZaa/dX8/0wwghEqcOHD6tZ\ns2bavn27Tj/9dNdxfoURxvB68cUXtX79es2YMcN1FCCkajvCGGhhbCVpsqQia+0IH47vpLJtA48V\nMGPMJM/Xmb4eK2mipAJJTaq6Ti3PbSpplqRN5VsfekYiiyX1qW6qncIIRK81a9bonnvu0caNG11H\nOQ6FMbx27dqlrl27au/evez2g5hW28IY0P8qPA/GZEjydbXbZJXtRV1xtK5AZdPc/hybIqnYy3VS\nanFub2vt/Ir7ZFtr96ns3szI+38TAAFbtmwZT0dDktS6dWs1atRIW7ZscR0FiEg1Fsaalsix1pZI\nyvPx85JV9rBMRUWSmvhxbKIP1+kUwLnHGGP6qOz+xiqnowFEt+zsbNZfxDH9+vXjaWmgGr6MME42\nxvy1hqV1fgxWIB8FMlHjz7mjJA0P4LMARKivv/5aX3zxhS699FLXURAhWF4HqJ4vT0kbSfdJGm+M\nKR9NzJE0p8LIW1MfP69Axxewpip72trXY4s972V4uU4g50qSjDH3SnrEl9HFBx988NjvU1NTlZqa\nWtMpABxbvny5+vTpo7p1fVosAnEgNTVVI0aM0IEDB9SoUSPXcYCgyM/PV35+fsDXqfGhF2PMKGvt\nFM/0bF9JvVU25WtVti5jnqQMa22Nm7AaYzpK2mitPaHCa9U99FLtsSp7cKXa6wRyrufroZJ2VvXk\ndhXfEw+9AFEoMzNT6enpuuWWW1xHqRIPvbjRu3dv3Xnnnbr66qtdRwFCImQPvVhrp3j+M89aO95a\nm6Ky4pUpaYXnP6u6B7Gqa22WlFdpy8HeksZJkjGmt2e3GK/HenlvfADnlmfoo7IRx32e7RE7eRYu\nBxAjjhw5otzcXB54wXEGDBigpUuXuo4BRJyAltU5dhFjJpYv7u3DsY1Vtmj2TpXtEJNrrZ1ffh1J\nrSqM9Hk7ttr3anuuMaa3pFwdf5/jOGvt36r5fhhhBKLMunXrdPvtt0f0E7GMMLqxbds2paWl6csv\nv5Qxfg/CABHPyTqMoDAC0eh//ud/dOjQIU2cGLmTBxRGN6y1SkpK0sKFC3XJJZe4jgMEnZN1GAEg\nGrH+IqpjjNHAgQO1ZMkS11GAiEJhBBBXvv32W+3YsUNdu3Z1HQURisIIHI/CCCCuLF26VH369FH9\n+vVdR0GESk1N1QcffKAffwz3EsNA5KIwAogrS5Ys0ZVXXuk6BiJYgwYNlJqaqpycHNdRgIhBYQQQ\nNw4dOqS8vDy2A0SNBgwYwLQ0UAGFEUDcWLt2rS644AI1b97cdRREuAEDBig7O1tHjx51HQWICBRG\nAHFjyZIlGjhwoOsYiALnnnuuzjzzTL377ruuowARgcIIIG68+eabFEb4bODAgez6AnhQGAHEhU8/\n/VQHDx5Ux44dXUdBlGB5HeD/UBgBxIUlS5ZowIABbPcGn3Xp0kVffvml9uzZ4zoK4ByFEUBcYDkd\n+Ktu3bpKS0tjWhoQhRFAHNi/f7/ee+899e7d23UURBmmpYEyFEYAMS83N1ddu3bVKaec4joKoky/\nfv20atUq/fLLL66jAE5RGAHEPJbTQW2ddtppateundasWeM6CuAUhRFATCstLdXSpUspjKg1pqUB\nCiOAGFdQUKCmTZuqdevWrqMgSpUXRmut6yiAMxRGADGNxboRqPbt2+uXX37RJ5984joK4AyFEUBM\nYzkdBMoYo0GDBumNN95wHQVwhsIIIGbt2bNHhYWF6tq1q+soiHKDBg3SokWLXMcAnKEwAohZixYt\nUv/+/VWvXj3XURDlUlNTtW3bNn377beuowBOUBgBxKw33nhDV199tesYiAH169dXWlqaFi9e7DoK\n4ASFEUBM2r9/v9atW6d+/fq5joIYcfXVVzMtjbhFYQQQk7Kzs9WtWzc1atTIdRTEiP79+ys/P18H\nDx50HQUIOwojgJjEdDSCLTExUZ07d1ZeXp7rKEDYURgBxJzDhw9r2bJlGjRokOsoiDFXX301y+sg\nLlEYAcScNWvWqE2bNjrzzDNdR0GMGTRokN58800dPXrUdRQgrCiMAGIO09EIlZYtW+qMM87QO++8\n4zoKEFYURgAxxVpLYURIMS2NeERhBBBTtmzZorp166pdu3auoyBGsesL4hGFEUBMWbRoka6++moZ\nY1xHQYxKTk7WgQMHtH37dtdRgLChMAKIKUxHI9SMMRo0aBDT0ogrFEYAMWP37t364osv1K1bN9dR\nEOO4jxHxhsIIIGa88cYbGjhwoOrWres6CmLcFVdcoa1bt+qbb75xHQUICwojgJgxb948DRkyxHUM\nxIETTzxRAwYM0IIFC1xHAcKCwgggJnz//ffavHmz0tLSXEdBnBg2bJjmzZvnOgYQFhRGADFh4cKF\nSk9PV8OGDV1HQZxIT0/Xhg0b9MMPP7iOAoQchRFATJg7d66GDRvmOgbiyEknnaS0tDQefkFcoDAC\niHrFxcV655131L9/f9dREGeGDh3KtDTiAoURQNRbtGiRevXqpVNOOcV1FMSZgQMH6q233lJJSYnr\nKEBIURgBRL158+Zp6NChrmMgDjVq1EhXXHGFFi9e7DoKEFIURgBR7cCBA8rPz9eVV17pOgriFNPS\niAcURgBRbcmSJerevbuaNGniOgri1FVXXaWVK1fqp59+ch0FCBkKI4CoxnQ0XEtMTFTXrl21dOlS\n11GAkKEwAohaP//8s3JycnT11Ve7joI4N2zYMM2dO9d1DCBkKIwAolZ2drZSUlJ02mmnuY6COHf1\n1Vdr+fLl+uc//+k6ChASFEYAUSsrK0uZmZmuYwBq1qyZkpOTlZ2d7ToKEBIURgBR6aefflJ2draG\nDBniOgogScrIyFBWVpbrGEBIUBgBRKXFixerW7duTEcjYgwbNkzLli3TwYMHXUcBgo7CCCAqMR2N\nSHPaaaepa9euLOKNmERhBBB1SkpKtGrVKl1zzTWuowC/MmLECM2aNct1DCDoKIwAos7ChQt1xRVX\nKCEhwXUU4FeuueYarVq1ir2lEXMojACiTlZWlkaMGOE6BnCchIQE9e7dWwsWLHAdBQgqCiOAqPLD\nDz9o3bp1uuqqq1xHAarEtDRiEYURQFSZN2+e+vXrp5NPPtl1FKBKV155pd5991199913rqMAQUNh\nBBBVmI5GpDvppJM0cOBAtgpETKEwAogaX3/9tTZv3qz+/fu7jgJ4xbQ0Yg2FEUDUeP3113X11Ver\nQYMGrqMAXqWlpenjjz/W7t27XUcBgoLCCCBqzJw5UzfeeKPrGECNTjzxRA0ZMoRRRsQMCiOAqPDh\nhx/qhx9+UGpqqusogE9uvPFGvfLKK7LWuo4CBIzCCCAqzJw5U9ddd53q1OHHFqJDt27ddPDgQb3/\n/vuuowAB4ycvgIh39OhRvfrqq7rhhhtcRwF8VqdOHd1www2aMWOG6yhAwCiMACLeqlWr9Jvf/Ebt\n2rVzHQXwyw033KDXXntNhw8fdh0FCAiFEUDEmzFjBqOLiErnnXeekpKSlJOT4zoKEBDDzbiBMcZY\n/gyB0Dl48KDOOussbd++XaeffrrrOGFjjMSPltgwadIkrVq1SllZWa6jADLGyFpr/D2PEUYAEW3B\nggXq1q1bXJVFxJaMjAxlZ2erpKTEdRSg1iiMACIa09GIdk2bNlXfvn01Z84c11GAWqMwAohYe/bs\n0YYNGzRo0CDXUYCA3HjjjTwtjahGYQQQsV555RVlZGTopJNOch0FCEi/fv20fft27dy503UUoFYo\njAAiUmlpqaZNm6Zbb73VdRQgYPXr19d1112nl19+2XUUoFYojAAi0urVq3XSSSepc+fOrqMAQTFy\n5EhNnz5dR44ccR0F8BuF0QtjTILrDEC8mjp1qkaOHClj/F79AYhI7dq1U8uWLbVkyRLXUQC/hX0d\nRk8Je1TSDkmnSsqz1q7w99iarhPgub0lTZY0yVr7RA3fD+swAkFWXFysVq1aaefOnTr11FNdx3GC\ndRhj0/Tp0zV//nwtXrzYdRTEqdquw+iiMM6WlGOtner5eoekYdba43Zn93ZsTdep7bnGmMaSjKRC\nSX+lMALh9/zzz2vt2rWaNWuW6yjOUBhj08GDB3XOOefogw8+0Nlnn+06DuJQVCzcbYzpJGmopNkV\nXs6TdL8/xxpjOkoaVt11AjnXWrvfWrtPUlEtvkUAQTB16lQedkFMOvnkk5WRkaHp06e7jgL4Jdz3\nMCZLKrHW7q/wWoGkPn4emyKp2Mt1UgI4F4BDmzZtUnFxsXr37u06ChASo0aN0rRp01RaWuo6CuAz\nF4Wx8shdkaQmfhyb6MN1OgVwLgCHXnrpJd1yyy2qU4dn8hCbOnXqpMTEROXl5bmOAvgsWn8iB3Jn\nD3cFARFq3759ysrK0siRI11HAULGGKNRo0ZpypQprqMAPqsb5s8rkDS80mtNJVW1I3t1xxZ73svw\ncp1AzvXbgw8+eOz3qampSk1Nre2lgLg2c+ZM9e3bV2eeeabrKEBIXXfddXrggQf07bff6vTTT3cd\nBzEsPz9f+fn5AV8nrE9Jex442WitPaHCa5MkJVprM309VtJEb9cJ5NwKr+0Qy+oAYWOtVbt27fTC\nCy/wjy7xlHQ8GDlypFq1aqX/+q//ch0FcSQqnpK21m6WlGeMGVLh5d6Sxkll6x8aY0bVdKyX98YH\ncO64oHyTAGpl9erVkqTLL7/ccRIgPO644w5NmjSJnV8QFVysw9hYZYtm75SUJCnXWjvf895ESa3K\nR/pqOLba9wI517Ood4akSZI2SRpvrV3p5fthhBEIgoyMDPXs2VO///3vXUeJCIwwxofu3bvr7rvv\n1tChQ11HQZyImoW7Yw2FEQjc3r171a5dO33++edKSGBHTonCGC+ysrI0adIkrVq1ynUUxImomJIG\ngKpMnTpVmZmZlEXEncGDB2v79u366KOPXEcBvGKEMUCMMAKBOXTokFq3bq0lS5aoffv2ruNEDEYY\n48eECRP0zTff6MUXX3QdBXGAEUYAUWnOnDlq27YtZRFx67bbblNWVpb27dvnOgpQLQojAGestXry\nySd1zz33uI4COHPGGWcoPT2d/aUR0SiMAJxZu3atfvrpJw0YMMB1FMCpu+66S08//TRL7CBiURgB\nOPPkk0/qrrvuYt9oxL1LL71U5557rubNm+c6ClAlHnoJEA+9ALWzY8cOXXbZZfr888918sknu44T\ncXjoJf4sWrRIDz30kN577z0Z4/czCYBPeOgFQFR59tlnNWrUKMoi4HHllVdq//79Wrt2resowHEY\nYQwQI4yA/4qLi5WUlKSPPvpIZ555pus4EYkRxvg0efJkLVmyRIsWLXIdBTGKEUYAUeOFF17QoEGD\nKItAJTfeeKPeffddffLJJ66jAL/CCGOAGGEE/HPw4EG1bt1a+fn5uvDCC13HiViMMMavCRMm6Kuv\nvtKUKVNcR0EMYi9pRyiMgH+eeeYZrV27VnPnznUdJaJRGOPXDz/8oLZt2+rDDz/UWWed5ToOYgyF\n0REKI+C7Q4cOKSkpSQsXLlRycrLrOBGNwhjfxo4dqyNHjuiZZ55xHQUxhsLoCIUR8N20adM0Z84c\nZWdnu44S8SiM8e3rr79Wu3bttHXrVv3mN79xHQcxhMLoCIUR8M3Ro0d1wQUXaNq0aerZs6frOBGP\nwog//OEPatCggR5//HHXURBDKIyOUBgB38yaNUvPPfec1q5dy6LEPqAwYvfu3Wrfvr0+/fRTnXba\naa7jIEawrA6AiHX06FFNmDBB//3f/01ZBHx0zjnnaPjw4Xr66addRwEYYQwUI4xAzWbOnKmXXnpJ\na9asoTD6iBFGSFJhYaFSUlK0Y8cOJSYmuo6DGMCUtCMURsC7w4cP64ILLtDLL7+syy+/3HWcqEFh\nRLmRI0eqefPm+utf/+o6CmIAhdERCiPg3ZQpUzR79mzl5ua6jhJVKIwot3v3bnXo0EEff/wxT0wj\nYBRGRyiMQPX+9a9/qW3btpo9e7a6dOniOk5UoTCionvuuUeHDh3Sc8895zoKohyF0REKI1C9Z599\nVrm5uVq8eLHrKFGHwoiKvv/+e1144YV677331Lp1a9dxEMUojI5QGIGqHThwQG3bttXSpUvVsWNH\n13GiDoURlU2YMEE7duzQzJkzXUdBFKMwOkJhBKr2xz/+UV9++aVmzJjhOkpUojCisgMHDui8885T\nbm6ufvvb37qOgyhFYXSEwggcr/wm/ffff1/nnHOO6zhRicKIqjzzzDNavny5li5d6joKohQLdwOI\nGA888IBuv/12yiIQZLfffrt27typZcuWuY6COMMIY4AYYQR+bePGjRo0aJC2b9+uRo0auY4TtRhh\nRHWWLFmisWPH6sMPP1S9evVcx0GUYYQRgHPWWo0dO1YTJkygLAIhMmDAALVs2VLPP/+86yiII4ww\nBogRRuD/zJ07V3/+85+1efNmnXDCCa7jRDVGGOHN1q1bdfnll2vbtm067bTTXMdBFOGhF0cojECZ\nAwcO6KKLLtJrr72mHj16uI4T9SiMqMmdd96pI0eO6IUXXnAdBVGEwugIhREoM3bsWBUVFWn69Omu\no8QECiNqUlRUpIsuukiLFy9W586dXcdBlKAwOkJhBKQPPvhAffr00ccff6xmzZq5jhMTKIzwxcyZ\nM/Xkk09qw4YNqlu3rus4iAI89ALAidLSUt1+++16+OGHKYtAmF1//fU69dRT9eyzz7qOghjHCGOA\nGGFEvJs2bZqmTJmidevWqU4d/g0aLIwwwlefffaZLrvsMhUUFKhFixau4yDCMSXtCIUR8WzPnj3q\n2LGjcnNz1b59e9dxYgqFEf74y1/+ovXr12vx4sUyxu8ugDjClDSAsLLW6rbbbtMdd9xBWQQcu/fe\ne/X555/r9ddfdx0FMYoRxgAxwoh49corr+ipp57Se++9p/r167uOE3MYYYS/CgoK1L9/f23evFln\nnXWW6ziIUExJO0JhRDwqn4rOyclRhw4dXMeJSRRG1MZDDz2kt956S9nZ2UxNo0pMSQMIi9LSUt18\n88264447KItAhLn//vtVUlKiF1980XUUxBhGGAPECCPizRNPPKGFCxcqPz+fdd9CiBFG1Nb27dvV\nrVs3rV+/Xuedd57rOIgwTEk7QmFEPCm/R2rDhg0s3xFiFEYE4vnnn9fLL7+st99+Ww0aNHAdBxGE\nKWkAIfXTTz/p2muv1bPPPktZBCLcf/zHf6hVq1YaO3as6yiIEYwwBogRRsQDa61GjBihRo0aaerU\nqa7jxAVGGBGoffv2KSUlRQ899JBGjBjhOg4iRG1HGLkBCUCNnnrqKe3cuVNvvfWW6ygAfJSQkKA5\nc+aob9++6tChgy644ALXkRDFGGEMECOMiHWrV69WRkaG3n33XbVs2dJ1nLjBCCOCZcqUKXrmmWe0\nfv16NWrUyHUcOMZDL45QGBHLvvjiC1122WWaPn260tPTXceJKxRGBIu1VmPGjNGePXv0xhtv6IQT\nTnAdCQ7x0AuAoNq/f7+uvPJK3XvvvZRFIIoZY/Tcc8/pn//8p8aNG+c6DqIUI4wBYoQRsejIkSO6\n6qqr1LJlS73wwgvsGOEAI4wItqKiInXp0kXjxo3TyJEjXceBIzz0AiAorLX6z//8T5WWlurZZ5+l\nLEsSUGoAAA3BSURBVAIxomnTpnrzzTfVo0cPtWjRQn379nUdCVGEKWkAv/Lggw9q/fr1mj17turV\nq+c6DoAgatu2rebNm6frrrtO69evdx0HUYTCCOCYp59+WrNmzVJ2drYSEhJcxwEQAt27d9eMGTN0\nzTXX6IMPPnAdB1GCwghAkvTyyy/rqaeeUm5urpo3b+46DoAQ6tevn/7+97+rf//++uyzz1zHQRTg\nHkYAmjx5sh5++GGtWLFC5557rus4AMIgIyND+/fvV69evZSbm8vC3vCKwgjEuWeffVZPPvmk8vPz\nlZSU5DoOgDAaOXKk6tatq169eik7O1uXXHKJ60iIUBRGII49/vjjmjx5slavXq0WLVq4jgPAgZtv\nvlkNGzZUWlqa3nzzTaWkpLiOhAhEYQTi0NGjR3XPPfcoLy9Pq1ev1llnneU6EgCHMjMz1bBhQw0Y\nMEAzZ85ksX4ch4degDhz8OBBDR06VB999JHefvttyiIASdKgQYM0f/583XTTTXrxxRddx0GEoTAC\ncWTv3r1KTU1VYmKili1bpiZNmriOBCCCdO/eXW+//baeeeYZ3XPPPTp69KjrSIgQFEYgTqxatUop\nKSkaPHiwXn75ZdWvX991JAARKCkpSevXr9eWLVvUr18/fffdd64jIQJQGIEYV1paqkceeUT/9m//\nphkzZuiBBx5guz8AXiUmJmr58uX63e9+p+TkZK1du9Z1JDhmLLvbB8QYY/kzRKTau3evbr31Vu3b\nt0+zZ8/W2Wef7ToSfGSMxI8WRIJly5bp3//933XnnXdq3LhxqluX52WjmTFG1lq/Rw0YYQRikLVW\nr732mjp27KguXbpo9erVlEUAtdK/f39t2LBBK1euVLdu3bRt2zbXkeAAI4wBYoQRkWbPnj268847\n9cknn2jGjBlKTk52HQm1wAgjIo21VpMnT9Yf//hHjRs3Tnfffbfq1avnOhb8xAgjEOcOHTqkxx9/\nXO3bt9eFF16ogoICyiKAoDHGaMyYMdqwYYPy8vLUsWNHrVy50nUshAkjjAFihBGuWWu1dOlS3Xvv\nvWrZsqWeeeYZnXfeea5jIUCMMCKSWWu1cOFC3X333ercubMee+wxtWrVynUs+IARRiAOrVmzRj16\n9NC4ceM0ceJELVmyhLIIIOSMMRo8eLC2bt2qiy++WCkpKbr99tv11VdfuY6GEKEwAlHGWqucnBz1\n7dtXN998s0aPHq0PPvhAgwYNYrkcAGF10kkn6U9/+pO2b9+uRo0a6ZJLLtGdd96pwsJC19EQZBRG\nIEocOnRIM2bMUPv27TV27Fhdf/31+uSTT3TDDTfohBNOcB0PQBw77bTT9Nhjj+njjz9WgwYN1Llz\nZw0fPlzr1693HQ1Bwj2MAeIeRoTatm3bNG3aNM2cOVMXX3yx7r33XqWnpzOaGOO4hxHR7MCBA5o+\nfbqefvppnXrqqbrllls0YsQIJSYmuo4W92p7DyOFMUAURoTCN998owULFmjmzJkqLCzUTTfdpFtu\nuUVt27Z1HQ1hQmFELDh69KhycnI0ffp05eTkqH///rr++uvVp08fnXjiia7jxSUKoyMURgTLzp07\n9eabb2revHn68MMPNWDAAGVmZmrAgAHsrBCHKIyINUVFRXr99deVlZX1/9u7v58orzyO4++vBeJg\nOgzWgNOoRbpJ0wwxC7rFCy+aAKm9Fat/QGuTvfGq227/ga2969262ya9rq4mvTJR6npjUhbsthYu\n/AVqq4KmzICxqIH97sU8MxkUZ3gG5ZkHPq/khHl+nZzJhOHDOed5DhcvXuTdd99l3759vPPOOyST\nyaibt2bEJjCaWRPwOXAVeAUYcPfvwp5bqZ5qrw3TvuB8BUapSjab5ezZs5w5c4bTp08zOzvL3r17\n6e/vp6+vT/99r3EKjLKaTU5O8u2333Ly5EnOnz/Pjh076O3tpbe3l+7ubhoaGqJu4qoVp8B4DDjt\n7l8F21eB/e7+Y5hzK9VT7bVh2hccV2CUiubm5vj5558ZHBzk+++/Z3BwkFu3brFnzx76+vro6+sj\nk8loXqIUKTDKWjE7O8v58+cZGBhgYGCAS5cu0dnZSXd3d7Fs2bJF34/PSSwCo5l1AUNAs7vPBPuO\nBtsHl3oucAS4AKQWq6faayvV+4z3pMAYY+fOnePtt99+bvXNz89z48YNRkdHGRkZYXR0lNHRUS5f\nvsxrr73G7t272b17N93d3WQyGQ01L8Pz/uxqzWoPjKv981vtXuTnl8vlGBoaYnBwsFjcnY6ODjKZ\nDB0dHXR0dPDGG2+wadMmBcmQqg2MK/3XaieQK4SxwAXyQS3MuQNAtkw9u6q8tlK9ssqE+dJzd6am\nprh79y6Tk5P8+uuvXL9+fUH55ZdfaG1tJZPJkMlk6Onp4fDhw7z55pu8/PLLL/bNrDEKHPGmzy/e\nXuTnl0qliiMvkP/unZiYYGRkhJGREYaHh/n666+5cuUKjx8/pr29vVja2tp49dVXSafTbN68mXQ6\nTWNj4wtp51oTRWCcemLfFJAKcW7zEurpquLapdQrMeXuPHr0iNnZWR4+fMjDhw958OABN2/e5NSp\nU0xPTy8oMzMz5HI57t27x+TkJHfv3uXevXs0NjbS2tpKS0sLW7dupa2tjbfeeosDBw7Q1tbGtm3b\nWL9+fdRvV0Rk1TAz0uk06XS6GCILcrkc4+PjjI+PMzY2xqVLlzh37hwTExPcuXOHO3fu0NDQQDqd\npqWlhebmZpqbm0mlUk+9TiaTNDY2LigbNmwgkUhoNIiVD4zPy3IGaspdW1W9X3zxBaXD0npd/eu5\nuTnm5uaYn59f9HW5Y/Pz8zx+/HhBKCy8fvToEfX19SQSCdavX08ikSCRSDAzM8Pt27dpamoimUzS\n1NREU1MTra2tJJNJWlpaaGlpKYZE3YgiIlI7UqkUnZ2ddHZ2Lnrc3ZmenmZiYoLJyUlyuRzZbLb4\nc2xsjGw2Szab5f79+/z+++/F8uDBg+Lruro6GhsbSSQSNDQ00NDQQH19PfX19cXXi+2rr6+nrq6O\ndevWPZdSGH5f7s+quPuKFeAQ8FulfZXOBT4oV0+111aq9xnvyVVUVFRUVFRU4lKqyXAr3cM4zNPD\nuzvJzx0Mc+6FCvVUe22lep9SzcRRERERkThZ0bWk3f2/wICZ7SvZ3QN8DGBmPWZ2qNK5ZY59soxr\nyx6r/l2LiIiIxFsUz2FMkn8w9jXgdeCMu58Mjh0BthceYVPh3GceW861leoVERERWWu0NKCIiKy4\nYDTJgWF/xsIIIrJ8Ztbk7tPLrWdFh6RXCzNrMrOjZvaRmX1mZj1Rt0mWzsy2m9mwmf0v+Ln47XVS\n88zsiJl9FnU7ZOnMrCtYQeuau3+lsBgfwXdn4W/f0cIUMqlNwTS/q+Rv3i3dX1WGietjdaL2JU8s\nH2hmz1w+UGpHyVrh7wMGfAUMm1lxVSCJh+BLbj9wPOq2yNIEq3ANAz3u/u+o2yOhnQH63f0nADOb\nMrMh/e2rPcH0umFg4yKHq8ow6mEMKfjC6weOleweAD6NpkUSUi/5G5x+Cn45esgHx13RNkvCMLPt\nQDv5JxtIfBwDjissxlY78ErJ9m9ALqK2SBnuPhMMQy9YjGQ5GUaBMbxnLVnYG1F7JAR3P+Hu10u2\npwnmUUXWKKnGh+7+JfmwLzEQTP1oB1Jmdizo1TgadbsklH8Bx82sMxiOHij9PpVYqDrDaEg6PC0f\nuIqYWS/5Lz0NR8dEMGfxb1G3Q0LbRf6fsw/d/UYwPWTczLLurhGaGHD3A2Z2mnzAOF54oonEStUZ\nRj2MstYdAt6LuhGyNEGvxjF3vx91WyS0FPmejRtQ7N3/J/l5qBIDQci/Rn4e+P7gUXiyRqiHMbwL\nPB0wNqJ5HLFjZn8BPlPvYqx8Amx/Yj1UM7OP3f2liNokSzPG070Y18gPU0s8fAe87+4/mdkw+eHp\nq4WbJyQWqs4w6mEML8zyhlKjzKyf/EPZdXdfjLj7H9z9pUIBTgCfKyzGwg9QvHvzqf1S24Lexa7C\nHdLufoL8nMadkTZMwqo6wygwhqTlA+MvmLe4EZgOnivWpaEVkRfL3cfJB4zS+Yp9wD+iaZGEEUwh\nyJrZH0t2p9CTCmJlORlGQ9LVeQ/43MzayS8f+ElhXo7UtuDZfafJT74vpcAfT1qqKkbc/aCZ/T24\ncWkK+I+GM2OlB/jUzIaATcAFfX61KegRPgBsBw6a2Q/ufjY4XFWG0dKAIiIiIlKWhqRFREREpCwF\nRhEREREpS4FRRERERMpSYBQRERGRshQYRURERKQsBUYRERERKUuBUURERETKUmAUERERkbIUGEVE\nRESkLAVGERERESlLgVFEpIaY2VUz+6hk+2Mzuxplm0REFBhFRGpLOzBWsr0LuBZRW0REAAVGEZGa\nYWbbAQd+KNnd/sS2iMiKU2AUEakdXQDufv2JfUORtEZEJKDAKCJSOxYMRwc9jrBwiFpEZMUpMIqI\n1I4/AbmS7Z2Au/uPEbVHRARQYBQRqSVdQLuZdZpZP/DXwoFgW0QkEgqMIiK1IwdMAcfI3/zyAZAz\ns2/c/USkLRORNc3cPeo2iIiseUEPYq+7/znqtoiIPKku6gaIiAgAB8n3KoqI1BwNSYuI1AYHvom6\nESIii9GQtIiIiIiUpR5GERERESlLgVFEREREylJgFBEREZGyFBhFREREpCwFRhEREREp6/+Y304v\nZMLfPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c99b518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Execute this cell\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=15, usetex=True)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "np.random.seed(seed=42)\n",
    "Nsamples=1\n",
    "measurements = np.random.normal(5, 1, Nsamples)\n",
    "mu = np.average(measurements)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# plot the distributions\n",
    "fig, ax = plt.subplots(figsize=(10, 7.5))\n",
    "dist = norm(mu, 1)\n",
    "C = (1./(2.*np.pi))**(10./2.)\n",
    "x = np.linspace(0, 10, 1000)\n",
    "plt.plot(x, C* dist.pdf(x), c='black',label=r'$\\hat\\mu=%.2f,\\ \\sigma=1$' % mu)\n",
    "#plt.plot([159.9,160.1],[0,0.8], ls=ls, color='k', label=r'%s $\\mu=%i' % (lab, mu))\n",
    "plt.axvline(x=5)\n",
    "\n",
    "plt.xlim(0, 10)\n",
    "plt.ylim(0, 0.5*C)\n",
    "\n",
    "plt.xlabel('$\\mu$')\n",
    "plt.ylabel(r'$L(\\{x\\}|\\mu,\\sigma=1)$')\n",
    "plt.title('Likelihood of $\\mu$')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quantifying Estimate Uncertainty\n",
    "\n",
    "As you've seen on in the previous plot, our ML estimate of $\\mu$ is not perfect. The uncertaintly of the estimate is captured by the likelihood function.\n",
    "\n",
    "We *define* the uncertainty on our MLEs as second (partial) derivatives of log-likelihood:\n",
    "\n",
    "$$\\sigma_{jk} = \\left( - \\frac{d^2}{d\\theta_j} \\frac{\\ln L}{d\\theta_k} \\Biggr\\rvert_{\\theta=\\hat \\theta}\\right)^{-1/2}.$$\n",
    "\n",
    "Taken together, these entries (more accurately, their squares) are know as the *covariance matrix*.\n",
    "\n",
    "This is the same as the **Fisher-information matrix** used to estimate the parameter uncertainities in cosmology\n",
    "\n",
    "The marginal error bars for each parameter, $\\theta_i$ are given by the diagonal elements, $\\sigma_{ii}$. These are the \"error bars\" that are typically quoted with each measurement. Off diagonal elements, $\\sigma_{ij}$, arise from any correlation between the parameters in the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In our example, the uncertainly on the mean is \n",
    "$$\\sigma_{\\mu} = \\left( - \\frac{d^2\\ln L(\\mu)}{d\\mu^2}\\Biggr\\rvert_{\\hat \\mu}\\right)^{-1/2}$$\n",
    "\n",
    "We find\n",
    "$$\\frac{d^2\\ln L(\\mu)}{d\\mu^2}\\Biggr\\rvert_{\\hat \\mu} = - \\sum_{i=1}^N\\frac{1}{\\sigma^2} = -\\frac{N}{\\sigma^2},$$\n",
    "since, again, $\\sigma = {\\rm constant}$.  \n",
    "\n",
    "Then $$\\sigma_{\\mu} = \\frac{\\sigma}{\\sqrt{N}}.$$\n",
    "\n",
    "So, our estimator of $\\mu$ is $\\overline{x}\\pm\\frac{\\sigma}{\\sqrt{N}}$, which is a result that you should be familiar with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is $\\pm \\sigma$? Errors as Gaussian Approximations to the Likelihood Function\n",
    "\n",
    "The result for $\\sigma_{\\mu}$ has been derived by expanding $\\ln L$ in a Taylor series and retaining terms up to second order (essentially, $\\ln L$ is approximated by a parabola, or an ellipsoidal surface in multidimensional cases, around its maximum). If this expansion is exact (as is the case for a Gaussian error distribution), then we've completely captured the error information.\n",
    "\n",
    "In general, this is not the case and the likelihood surface can significantly deviate from a smooth elliptical surface. Furthermore, it often happens in practice that the likelihood surface is multimodal. It is always a good idea to visualize the likelihood surface when in doubt (see examples in §5.6 in the textbook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is $\\pm \\sigma$? Errors as Gaussian Approximations to the Likelihood Function\n",
    "\n",
    "The $(\\hat \\mu - \\sigma_\\mu, \\hat \\mu + \\sigma_\\mu)$ range gives us a **confidence interval**.\n",
    "\n",
    "In frequentist interptetation, if we repeated the same measurement a hundred times, we'd find for 68 experiments the true value was within their computed confidence intervals ($1 \\sigma$ errors)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### MLE applied to a Heteroscedastic Gaussian\n",
    "\n",
    "Now let's look a case where the uncertainties are heteroscedastic.  For example if we are measuring the length of a rod and have $N$ measurements, $\\{x_i\\}$, where the uncertainty for each measurement, $\\sigma_i$ is known.  Since $\\sigma$ is not a constant, then following the above, we have\n",
    "\n",
    "$$\\ln L = {\\rm constant} - \\sum_{i=1}^N \\frac{(x_i - \\mu)^2}{2\\sigma_i^2}.$$\n",
    "\n",
    "Taking the derivative:\n",
    "$$\\frac{d\\;{\\rm lnL}(\\mu)}{d\\mu}\\Biggr\\rvert_{\\hat \\mu} = \\sum_{i=1}^N \\frac{(x_i - \\hat \\mu)}{\\sigma_i^2} = 0,$$\n",
    "then simplifying:\n",
    "\n",
    "$$\\sum_{i=1}^N \\frac{x_i}{\\sigma_i^2} = \\sum_{i=1}^N \\frac{\\hat \\mu}{\\sigma_i^2},$$\n",
    "\n",
    "yields a MLE solution of \n",
    "$$\\hat \\mu = \\frac{\\sum_i^N (x_i/\\sigma_i^2)}{\\sum_i^N (1/\\sigma_i^2)},$$\n",
    "\n",
    "with uncertainty\n",
    "$$\\sigma_{\\mu} = \\left( \\sum_{i=1}^N \\frac{1}{\\sigma_i^2}\\right)^{-1/2}.$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Goodness of Fit\n",
    "\n",
    "The MLE approach tells us what the \"best\" model parameters are, but not how good the fit actually is.  If the model is wrong, \"best\" might not be particularly revealing!  For example, if you have $N$ points drawn from a linear distribution, you can always fit the data perfectly with an $N-1$ order polynomial.  But that won't necessarily perfectly predict future measurements (see cross validation for how to estimate the number of parameters).\n",
    "\n",
    "We can describe the **goodness of fit** as whether or not it is likely to have obtained $\\ln L_0$ by randomly drawing from the data.  That means that we need to know the *distribution* of $\\ln L$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Goodness of Fit\n",
    "\n",
    "For the Gaussian case we have just described, we can write\n",
    "\n",
    "$$z_i = (x_i-\\mu)/\\sigma,$$ then\n",
    "$$\\ln L = {\\rm constant} - \\frac{1}{2}\\sum_{i=1}^N z^2 = {\\rm constant} - \\frac{1}{2}\\chi^2.$$\n",
    "\n",
    "So **$\\ln L$ is distributed as $\\chi^2$ (with $N-k$ degrees of freedom).**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Goodness of Fit\n",
    "\n",
    "What does that tell us?\n",
    "\n",
    "The expectation value for the $\\chi^2$ distribution is $N − k$ and its standard deviation is $\\sqrt{2(N − k)}$. We typically have $N \\gg k$ (where $N$ is the number of data points, and $k$ is the number of parameters in the model). When that holds, it becomes useful to define **$\\chi^2$ per degree of freedom, $\\chi^2_{dof}$**, as:\n",
    "\n",
    "$$\\chi^2_{dof} = \\frac{1}{N-k}\\sum_{i=1}^N z^2_i.$$\n",
    "\n",
    "Therefore, for a good fit we would expect that $\\chi^2_{dof}\\approx 1$ (the expectation value).  If $\\chi^2_{dof}$ is significantly larger than 1, then it is likely that we are not using the correct model.\n",
    "\n",
    "We can also get overly high or low values of $\\chi^2_{dof}$ if our uncertainties are under- or over-estimated. This is particularly interesting when we have correlated uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Goodness of Fit: An Example\n",
    "\n",
    "![Ivezic, Figure 4.1](http://www.astroml.org/_images/fig_chi2_eval_1.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Working with non-Gaussian Likelihoods \n",
    "\n",
    "As an example of MLE with non-Gaussian probability density we can use the same formalism above for a Poisson distribution. In this case we write the probability disrtibution as\n",
    "\n",
    "$$p(x_k|\\mu) = \\frac{e^{-\\theta}\\theta^{x_k}}{x_k!}$$\n",
    "\n",
    "with $\\mu$ the average number of events and $x_k$ is the number of events observed\n",
    "\n",
    "This distribution is particularly useful for characterizing the number of soldiers in the Prussian army killed accidentally by horse kicks\n",
    "\n",
    "We can then write the likelihood  as\n",
    "\n",
    "$$L \\equiv p(\\{x_k\\}|\\mu) = \\prod_{k=1}^{n} \\frac{e^{-\\mu}\\mu^{x_k}}{x_k!}$$\n",
    "\n",
    "and the lnL as\n",
    "\n",
    "$$lnL = \\sum_{k=1}^{n} ln \\bigg( \\frac{e^{-\\mu}\\mu^{x_k}}{x_k!} \\bigg)$$\n",
    "\n",
    "$$= \\sum_{k=1}^{n} -\\mu + x_k \\; ln(\\mu) - ln({x_k!})$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEgCAYAAAB8XO79AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wVOed7/n3V2BiULAEDJ5J/ANakpP41x0jyVuepDLR\nWMiTmt3c3ZhfNdezVTsVYdl3995KbsDYrmyNXZu9BuJMcu+tuRgj193ZGu9dwMpMTbZmYiw8qnXF\n8Q4SOBNj4iCpBcROMhj9ADdgDHr2j/O0fGi61edIre5W9+dV1eXuc87znKebdn/1/DbnHCIiIoVQ\nU+oCiIhI5VBQERGRglFQERGRglFQERGRglFQERGRglFQEalQZpYws7pSl0Oqi4KKlAUzazezQTOb\nNLO9ZvacmR0ws60R0+8zs2fmupwRyhF+Hy+b2S7/Xvaa2bos1xes3OG8/Oc2BLQUMl+RfEzzVKRc\nmNl2YKtzboF/XQckgUPOuT/Mk/Z+YNQ59+bcl3R66fcB1DvnzvljCeAVYNw51xq6NlK5zWwNsMw5\n9+o011yVl5lNAmunSxPlXuX02Ur5U01FypZzbgLYB6yNcO2rZfSjd8b/19IHnHNJoANoDv/VH6Xc\nPiAdzHfTQnwG2e5VZp+tlDkFFSl3K4BxADNrMLPtZrbVN8lMNY2Z2Toz2xx6vd3MOn3T064IxxPZ\n8vbNWQfMbLNPM2pme2fyRnxg6QUei1nuZqAe6PLn1oXKdMDMjmfLy2v014ya2T5/3WbfPPegf73V\nv74n8145ypj1syr05yXzlHNODz3K4gFsB66EXjcDV4CvAXXAKLA0dH4U2AKsAfqBvf54AugPXbcl\nz/H6XHmHnr8M3ODvNQncM8372OrLfUOu9wisjlpu/3wS+IOM8r3s83kwM69Qmr0Z5X4mdO7BjGvv\nybxXljJO+1nN5PPSo7IeqqlI2fGd2/uAbQR9Ai8Am4AzzvdReM8DXc65IwQ/fMBUjaDZ923gnHt2\nuuPAhlx5++ejwCvOubP+XgANM3x7Q6FyRip3iIWejwL7nXMjzrkfZOYVsjtU7ueB9RHLadnKCGxk\n+s8qXbZCfV4yzyioSNlxzj3qnNvonNvknPsHfzjbj9JQjuMQ1HIeNrMzvuN5uuONMfMez/MWptNK\n0Fk/kuN8rnIXwgCwfJZ5xP13gNl9XjLPKKjIfJHrh+tw5gEzq3NBx/Jqgo7+gRzH03+BR857Nsys\nnqBWtDvH+azlLrBstZk4ivJZyfyloCLlxHKdcM7tAcbTncdeB/DvQ6/r/X8bzGyzb355FBgzsxuy\nHB83sxsi5h1Wn+N4zvfhA0ov8I/OuSdz5Jer3ADDBJ3u+SYzhsuWOV9gAx8PEhgmeI+Y2cP+2obQ\nucx71cO0/w7TzWPJ93lJBVlYzJv5L+kOYJBgVE+vcy7vUEnfxuycc0/MJh8pX2bWDqzzz3cBO7I0\nEbUAu82skWDY7v/tnPtrn3YtsMzPqRgDtpuZA5YB25xzZ82MbMcj5J0AOszsJYIfUAdsMrPeUPrw\n+3jYv9xvZsMEQSYB/O/Oub/OuDZvuf3lLxF08teZWdLn12Vm/c65NzPzcsEck23ANjPr8Hlsd879\n1D/fBjxvZmsJ+kMOA/cCP8i415Es+Wb7rH4Qek+RPy+pPEWd/Og7Xw8457r960FgvZtmDLz/ku4m\n6JR8Yqb5iIjI3Cta85eZNRP8JbovdLgXeGKaNAmCKvlA6FjsfEREpDiK2afSQjDqJVz9HWD62dIP\n+zbccBt1evRMnHxERKQIih1URjOOjZKjE8+CpSyydZQ2x8lHRESKpyxHf/klIfZlTLASEZEyV8zR\nXwMEQxrDlpN9YtQ2IOFH66SZmT1GMFIlaj74kTQiIhKTcy7nMP9cillT6efaJqoWgk72qzjnmpxz\nC9IPoIdgiOkCguAUKZ9Qfno4x5/92Z+VvAzl8tBnoc9Cn8X0j5kqWlBxwRpAvemVUb12/GQsv7pp\n5gqrsfMREZHSKerkR4Jmqx1m1kCw3tI259wJf66DYNLUnizpss0MzpWPiIiUSFGDiguGAT+a49zj\n06TbFDUfya2tra3URSgb+iw+ps/iY/osZq/itxM2M1fp71FEpNDMDFfmHfUiIlLhit2nIiIVYvXq\n1Zw4oa7M+WrVqlWMjIwUPF81f4nIjPjmkVIXQ2Yo37+fmr9ERKTkFFRERKRgqqJPZXQ0c/3J6V13\n3XUsXbp0jkojIlK5qiKo/OM/TsS6vqZmlI6OZjLWHhMRkTyqIqisXJmIdf3774/NUUlEKt8bb/yM\n8fFLpS5GLPX1i7jvvrtLXYyKUBVBRUSKZ3z8EitXtpS6GLGcPj2Q/6IyMjExQV1dXamLkZU66kVE\n5omDBw/S1NTEnj3ZlkgsD1VRU0mloLY2+vXnz9ewZ88Qv/Vb51iyZDJyOlWhRWSunD17ltbW1tgD\nj4qtKoJKZyd0d0cLLKkUfPObn+HkyetpaKiJnA7mXxVaRGanqamJRx55hC1btgCwc+dOnn/+eQYH\nBwt+rxtuuAGA5cuXFzzvQqqK5q9kEoaHo107NAQnT17PlSs1sdKJSPUZHh6moaFh6nV/fz+NjY0l\nLFHpVUVNJZGA0L/7tBob4dZbL3Lq1PUkEjWR04lIdUkmk5gZzc3NU8eGh4fp6OjImeaRRx5hbCwY\nXRpeIiW9ZIqZ0dHRQWdn59wVfI5VRVCJ04RVWwvf/e4vGBn5BC0tt8fqixGR6nH48GEgWFgzfOzJ\nJ5/Mmea5556b62KVXMU0f5lZzskocQPDkiWTfO5zKQUUEckps+krmUwCXHWsGhW1pmJmdcAOYBBY\nAfQ65w7muDYB7AZagVFgg9+fPn1+P7DOv3TA/BoYLyLz2qFDh6ivr596PTAwgJlxzz335EwTbv7K\npOavmdkDHHDOdQOY2aCZrXfOvZnl2scIgoYBLwH7gNt8ujrgfWCNPz/unBuZ++KLiAQOHz7M2NgY\nR44cYXh4mO3bt0+d6+npYd26ddekqYbmr6IFFTNrJggS4RDcCzwBbMq49gbn3KOh17uB8B72O4Bx\n4KxzLjlnhRaR2OrrF8274fX19YtmkKYeM2Pjxo3s2LGD7u5u2tvb2bRpE3v37i14GScmJti3bx/J\nZJK9e/fS3NzM/fffX/D7zFYxayotBDWKs6FjA8D2zAszrgG4F9gcet3g83vMzA4TNI0puIiUgWqY\nANzT08O9997Lrl27rjp+5syZObtnXV0dmzdvZvPmzfkvLqFidtS3EPSNhI0C9VmuBcDM2s1sENhK\n0LcCgHPuAefcCmAtUAccKHxxRUSy27t3b9nPbC+Vsh795Tvx1xL0qew2s3syzr9KEGwazKz86oEi\nUpHMjE2bNuW/sAoVs/lrANiQcWw5Qd9ITr4DfpOvsbQCb2acnzCzg0xT49m9+6mp5y0tbbS2tsUo\ntojI1eaiz6TU+vr66Ovrm3U+Nt3G94VkZmuAfufcgtCx54Blzrm8Id/MDgDbfe0k81w/sD7bCDAz\nc/398d7j++8fZnJykhtvbM1/ccjp0wN8+csa2SzVIT0LXOanfP9+/nzsnQqL1vzl55j0mtmDocPt\nBEOH0/0nm/3zOjNrT1/khxCPOedeNbOEma0LnVsDHNKQYhGR0iv2PJUNwA4zawAagW3OuRP+XAeQ\nIJjLspygD8URDB8eDdVm6n0e2wkmRw6Hhx+LiEjpFDWo+KHCWQOAc+7x0PMk0JTjuiO5zomISGmV\n9egvERGZXxRURESkYKpi6XsRKZ43Bt5gPDXtTIGyU19bz30t95W6GBVBQUVECmo8Nc7KppWlLkYs\npwdPl7oIeSWTSbq6uujv72f58uXs37+fNWvWlLpY11Dzl4jIPLBz5056enoYGRmhoaGBjRs3lrpI\nWSmoiIjMUFNTE88+++zU6507d9LUVPjBqWfPnmXXrl0sXbqUG264ga6urqv2ciknCioiIjOUuftj\nf38/jY2NBb/PDTfccNXrQ4cOsWfPnoLfpxDUpyIiMgPJZBIzo7m5eerY8PAwHR0dOdOEd34ML5GS\nXjIl386PBw8epKuri2QySVNT07S7TJaKgoqIyAwcPnwYgNWrV1917Mknn8yZZrY7P7a3t9Pb28u2\nbdvo6uqitbW17AKLmr9ERGYgs+krmQz2CQwfmwurV69m7969JBIJ+vv75/ReM6GaiojIDBw6dOiq\nzvKBgQHMbNqaQ7j5K1OU5q+whoaGOQ9gM6GgIiIyA4cPH2ZsbIwjR44wPDzM9u0f74ze09PDunXr\nrkkz0+aviYkJ+vv7aW9vn3q9bNmyqt+jXkSqQH1t/byYTBhWXxt/eG59fT1mxsaNG9mxYwfd3d20\nt7ezadOmgm/iNTo6SldXF2bGtm3bWL58edluFKagIiIFVQ3LnfT09HDvvfeya9euq46fOXNmTu6X\nSCQYHByck7wLTUFFRCSmvXv3YhZ7U8SqoNFfIiIxmRmbNuXdBb0qqaYiIhJTufZnlIOiBhW/1/wO\nYBBYAfQ65w7muDZBsF1wKzAKbPC7PsbKR0REiqfYNZU9wAHnXDeAmQ2a2Xrn3JtZrn0MWAcY8BKw\nD7htBvmIiEiRFK1PxcyaCYLEvtDhXuCJLNfe4Jx71Dl3zu9rvxsYj5uPiIgUVzE76luAcR8k0gaA\ntZkXZlwDcC+wOW4+IiJSXMUOKqMZx0aBnLOOzKzdzAaBrQR9KzPKR0REiqOshxT7zve1BH0qu82s\nvJbjFBGRqxSzo34A2JBxbDm+ryQX59wIsMnMjhPUVmLns3v3U1PPW1raaG1ti1hkEYnrZ2+8waXx\naf+3LjuL6uu5+77KXwlgOn19ffT19c06n2IGlX6ubaJqIehkjyIJDANjcfPp6noq4i1EZLYujY/T\nsnJlqYsRy8Dp+bVW2cTEBHV1dQXNs62tjba2tqnXTz/99IzyKVrzl59j0mtmD4YOtxMMHU73n2z2\nz+vMrD19kZ+XMuacezVfPiIilergwYM0NTVds5VwMpmktbWVmpoaWltbOXLkSIlKWPx5KhuAHWbW\nADQC25xzJ/y5DiBBMAdlOUEfiiOY5DjqnNsUMR8RkYpz9uxZWltbGR29epzSxMQE27Zt44UXXsA5\nR2dnJ62trYyNjV2zt30xFDWo+GHAj+Y493joeRJomkk+IiLF0tTUxCOPPMKWLVsA2LlzJ88///yc\nrCicDhDLly+/6nhvby87d+6c2tb44MGDLFu2jP7+/pLst6K1v0REZihzS+H+/n4aGxuLWobMzcDq\n6uowM1pbW3OkmFsKKiIiM5BMJjEzmpubp44NDw/T0dGRM014O2Hn3NRxM4u9nXAuvb29rF27tiRN\nX6CgIiIyI4cPHwaYanZKH3vyySdzppnpdsJx7Nmzh/3798/5fXIp68mPIiLlKrPpK5lMAlx1rNi+\n853v8MQTT5SslgKzrKn4hR8z1+kSEal4hw4dor7+4ylzAwMDmBn33JN74Y9w81em2TZ/9fT00NHR\nMe39iyFyUDGzXcCQc+5ZM1tDMLN9zMy2pZegFxGpFocPH2ZsbIwjR44wPDzM9u3bp8719PRc04EO\nc9f81dvby+joKM3NzSSTScbGxti3b99VZSqWODUVc84965/3Arudc4+a2dY5KJeIzFOL6uvn3Qz1\nRfXx16Otr6/HzNi4cSM7duygu7ub9vZ2Nm3aNCc7Q05MTLBv3z6SySR79+6lubmZ+++/n4MHD/LA\nAw9gZlddv3PnzoKXIQoLj0CY9kKzTudct5ltBx4GVjvnzqaPz2kpZ8HMXH9/tPeY9v77h5mcnOTG\nG+MNyTt9eoAvf7klVhqR+So9Yqka9fT00Nvby65du0pdlBnL9+/nz1vOC3KIU1NpMrN+oBlo9gFl\nM7AdKNugIiJSaHv37r2mZiCByEHFOfe4mSX8bPf0elz9aHMsEakyZsamTZvyX1iF4nTUX9XM5Zyb\nMLNhgh0ZS7d6mYhIkc1Fn0mlyBtUzCy9eEyHmQ0B4TpfA8He8M9ek1BERKpOlJrKGHCQYA+TzM2x\nhglWFRYREckfVJxzR8yshaBzvifzvJmVbuqmiIiUlUh9Kr5zPun3iM9cg2CTf4iISJWL01G/l4+b\nv9IbUNcD1TlQXURErhFnnkoH0OK3852S3gK41MLDnUVk7q1atUpzNeaxVatWzUm+cYLKPmAoy/FX\nombg57bsAAaBFUCvc+5gjmsTwH6CyZaHgc3hgGZm+4H04joO0FR2kSIaGRkpdRGkDMUJKgPAfjPL\nDCIdwB9GzGMPcCA938XMBs1svXPuzfBFoeDzNYIhzN1Av5kt8zP564D3gTX+/LhzbiTGexERkTkQ\nJ6hsAFqBzL0yE1ESm1kzQc0ivKZzL8E8l8yO/rXAY+lAYWbtBEObW4FXCQLOOHBWTV4iIuUjziZd\nO5xzy51zTeEH8EDE9C0ENYrw/isDZFnmxTnXE655OOcmCJq4DvlDDQQz+YfMrN83lYmISInFqamc\n8UOKwxqAewkmR+bTAoxmHBslGEE2LTNbS9D/cg7AOfeAP34/sBs4ANwWoQwiIjKH4gSVw1mOOX/8\nicIUJ6fNXDubH+fcq2bWCoya2f3OuVfnuBwiIjKNOEFlm3PuO+EDZvYwwUiuKAa4NjAs5+M5L1n5\nTcCeybVtsV/YMr2MTFa7dz819bylpY3W1rZoJRYRqRJ9fX309fXNOp/Im3TlzMDsuHMub9OT34K4\n3zm3IHTsOWCZcy7rjHwzW0ewhfGb2c6HrusH1mcbAaZNukRE4pvzTbrM7MEshzsIaht5+TXEes3s\nQefcD/zhdnxHvR/h1eCc2+Nfr/V5H/Yd8cuAjQR9KFPrkPlgdUhDikVESi9O89dLWY4NE2wtHNUG\nYIeZNRAMTd7mnDvhz3UQDE/e4wPMAa5dAuYxgmauHX5b493AsHPu0RhlmDPHhv4J/t94e3OfOnWK\nW265Jfa96mvrua/lvtjpRETmUpyg0pWuRcyU7xfJGgCcc4+Hnh9k+uHOTbMpx1xJffgBK5tWxkrT\n/7N+mv+gOfa9Tg/GC14iIsUQeZ6Kc26PmW0xs+NmdsXMDpnZV+eycCIiMr/E6VPZSjB0+Hk+XgPs\nj83MQn0kIiJSxeI0f93rnMvslN9jZlsKWSAREZm/4izTcijH8RWFKIiIiMx/cYLKCjNbGj5gZp0E\nS9OLiIjEav7aDZwwszP+dQPBbHjN9hMRESBGUHHOJf0kxI0Ec0yGZjvEWEREKsu0QcXMVvun9QST\nDCcINtqaOq+Z7IHU5XOcvJIkdel3qV1UGy3NpRQnPzpJ6lIqchoRkXKWr09lmGD48DV7nniPmNkN\nhS3S/JO6fI7O17/IC+f/A50/7CR1KZU/zaUUnT/s5IWJFyKnEREpd1E66h9wzj2bY5XgZ5j7Ze/L\n3tC5t0ieO8okkyTHkwyPDedPMzZEcjwZK42ISLnLF1QO+yVTsvLNYXk32ap0jUvvIrH0TmpYQKI+\nQcOyhvxpljWSqE9QQ03kNCIi5S5fUBnKcx4irlJcyWoXLqX786/RueTf0v2V7kj9I7WLaun+Sjed\ndZ2R04iIlLt8QSVKLUR/YhMEllsWJGIFh9pFtdxy3S0KKCJSMfIFlaSZfTPXSTPbRdCZLyIiknee\nyuME+793EEx+TBI0dzUAXQSz6RvntIQiIjJvTBtUnHPjZtYEvAL08PGmWUYwm36j5qmIiEha3hn1\nzrlhoNFv77uGYAHJQ+ntfOMwszpgBzDo8+nNNbrMz97fT1AbOgxsds4diZuPiIgUT5xlWnqB3lne\nbw9wwDnXDWBmg2a23jn3ZviiUND4GkGtqBvoN7Nlfr5MpHxERKS44qxSPCtm1gysA/aFDveSffLk\nWuAx59xPfaBoJwgurWa2BlgfMR8RESmiogUVgtWMxzNm5g+QZQkY51xPuK/GT7J0QD9wLzAWJR8R\nESmuWQeVGGt/tQCjGcdGiTAXxvfn9PpA0jzTfEREZG7l7FPxHeVr8qRvIOgon+ump83Ahjm+h4iI\nzFK+jvqXCIYOA9Tx8VDitHqCpqcoQWWAawPD8oz8rmFmW4FnQs1dsfPZvfupqectLW20trZFKK6I\nSPXo6+ujr69v1vnkDCp+U66u9EZcZrYd2O2cS6avMbMG4MGI9+rn2iaqFqYZUWZm64BXMkZ1xc6n\nq+upiEUUEalObW1ttLW1Tb1++umnZ5TPtH0qGTs7vh8OKP78MLApyo38HJNeMwsHoXbgMQAzazez\nzekTvh9lOTBhZgkzazaz7fnyERGR0omzR32TmS11zp1LHzCzTuJ1kG8AdvgaTiOwzTl3wp/rABLA\nHjNrBw7w8Qz+tHTg2Ahsz5GPiIiUSJygsgM4YmZjBKOtWgkCSuQOdN8v8miOc4+Hnh9kmlqUH2Kc\nNR8RESmdODPqkwS1lfUEAaUXeCmzSUxERKpXnJoKAM65l4CXzKzO1xhERESAGJMfzWy1mb1sZi+H\nju2KMfmxoqVScHLkd0h9UMxFCkREykucX8DnCeapDMNUv8ZLBIs7VrVUCjo74YW/eJDOr35WgUVE\nqlacX79x59wDXL3T4zDB4o5VbWgIkkmYnFxA8vj1DL+zuNRFEhEpiThBJR1MwsN8t6HthGlshEQC\namqukLjtIg2fvVDqIomIlEScjvq9ZnYAcGbWSLAqcCPw8JyUbB6prYXubnixp4eH/ufbqP3kZKmL\nJCJSEpFrKn4m+waCfhQj2F64Jb1RVrWrrYVbVv1GAUVEqlqsIcW+c34Poc55M1utfepFRARiBBU/\ndHgjQZNX2FqCjbNERKTKxampvEqwf0p/xvHMICMiIlUqbvPX8sxjfs94KaLUpRTHJo7xhQ+/wNJP\nLC11cUREpsQZUpxrv5LMlYRlDqUupej8YSdbj2zli//li5z78Fz+RCIiRRKnpvKMX6LlcOhYPUGf\nym0FLZXkNDQ2RHI8yRV3hbdPv83R00e57+b7Sl0sEREgXlDZR9AhbxnHVxSuOJJP47JGEvUJkmNJ\n7lh5B3euvLPURRIRmRInqDSqT6X0ahfV0v2VbgZ+NkDnA53qUxGRshInqOzOcXysEAWZTpRl9s0s\nUeq9XU6eeo/Xf3IsZppfx75P7aJabq+7XQFFRMpOrNFfWfpUIMY8FTOrI9hBcpCg2azX7/KY6/p2\ngmD2HPBsxrn9wDr/0gEtUcowly5edCyrvz1mmh/PUWlERIovTlDpAJZz7Y93Q4w89gAH0ku7mNmg\nma13zr2ZeaGfbNnv75l5rg54H1hD0Mczrln9IiKlFyeobHfOvZp50Ncm8jKzZoKaRWfocC/wBLAp\n83q/nz1mNpolux3AOHC21E1eIiLysTjzVHab2YOZB6drvsrQQlCjOBs6NkDQfBZXA7AZGDKzfjNL\nzCAPEREpsDhBJcm1/SmY2T0R07cAmbWOUYK5LrE45x5wzq0gCEh1wIG4eYiISOHFaf7aT1Bb2Z9x\nvIsSLSjpnHvVzFqBUTO7P1vz3Eycv5Ii+cFxapd/ltqF5TnC6tgv4o0yA6ivree+Fk2UFJG5Eyeo\nbABauXYByahNTwM+j7DlBH0jM+acmzCzg0xT49m9+6mp5y0tbbS2tuXML3X5HN88+jVOnk/S8O5/\npPvzr5VlYEl9mGJl08pYaU4Pnp6j0ojIfNfX10dfX9+s84k1T8XvUX8VM1uX7eIs+rn2h7+F3GuK\nxbGcLE1zaV1dT0XOaOjcW5w8n+QKl0mee5vhc0e5e5n+uheRytbW1kZbW9vU66effnpG+cTZ+bHH\nzLaY2XEzu2Jmh8zsq865nojpjwC9GZ397cBjEIwiM7PN+fIxs0Q4kPkZ/YcKNaS4celd3LokwUIW\nklh6Bw1LtQyKiEhUcTbp2kow/Pd5YMgf/mMzM+fcDyJmswHYYWYNBM1o25xzJ/y5DoKmtD3+fnUE\nm4IlgE1mdtj3mdT7PLYTTIwcds49GvV95FO7cCnfvfMFRlKDtNy6viybvkREylWc5q97s6z9tcfM\ntkTNwA8nzhoAnHOPZ7y+Zutif/wI0BT1njOxZEEtn/vk3QooIiIxxQkqh3Icr8hVio8N/YzB996P\nleb02LtzVBoRkfkhTlBZYWZLnXNTu0KZWSfQXPhilV7qw3Mkbv8XsdJcunJxjkojIjI/xF2l+ISZ\nnfGvGwiGA5d8IUcRESkPkYOKcy7pl0PZRBBQhpxze/IkExGRKjJtUDGze8IrCPvO8+fnvFQiIjIv\n5aup7DCzbeSZ9a5l50VEBPIHlXsJ1vbK3HXR4SctAjsJ5q+IiEiVyxdUtmX2m/h+lf0Em2NtiDqj\nXkREKt+0QSVLQLmfIKAAtPiJiCIiIkCMtb/8zPlegn1VGhRQREQkU94hxX6v+P0Ea3M975x7ZM5L\nJSIi89K0NRUzW83HW/4+nC2gxNj5UUREKly+msowwUivncCw71MJawQepkQ7P4qISHnJF1TGgW3+\neeaOj+ljywpaIpFydO4cvPUW3HUXLI2xevVM04nMU/mCygbn3MHpLjCzVwpYHpHyc+4cfPGLcPQo\n3HknvPZatAAx03Qi89i0fSr5AkrUa0TmtbfeCgLD5cvw9tvB87lMJzKPxVmluGTMrM6vOyZSfHfd\nFdQ03n4b7rgjeB4j3eTRo1y8+Wbeee89Jn/0o0hJF9XXc/d9982i0CKlUdSg4rcI3gEMEmzu1Ttd\nTcfM2gmW3H8OeHam+YjMytKlQdNVuhkrahOWT/eL7m4+19LCmtrayLccOH16hoUVKa1i11T2AAec\nc90AZjZoZuvDKyGn+fkx/UDmFsax8hEpiKVLYSY1h6VLSd1+O8QIKKRSjP/933Pk/HkmlyyJdTvV\ncKTUihZUzKwZWAd0hg73EixGuSnzer+fPWY2Opt8RK5RziOyUino7KRtaIgFf//30N0dKyCphiOl\nVsyaSgswng4W3gCwvUT5SDWaxYisn73xBpfGp90FIqvksWO0rFwZ7eKhIUgmWTA5CckkDA/D3XfH\nvqdIqRQ7qIxmHBsF6kuUT1VJXUpxbOIYX/jwCyz9RJn9dV5M2UZkRWwuujQ+Hj04hBzv749+cWMj\nJBJcGRrssmrAAAAUVklEQVRiQSIBDQ2x7ydSSvNi9JfMTupSis4fdjI8NsxfvvuXvPanr1VvYJnF\niKxYNY6Zqq2F7m5effFFOh56KF5fjEgZKGZQGQA2ZBxbTp5dJQuRz+7dT009b2lpo7W1LeYt57eh\nsSGS40muuCu8ffptjp4+yn03V2ln7ixGZMWqccxGbS1nbrlFAUWKqq+vj76+vlnnU8yg0s+1TVQt\nBJ3sc5pPV9dTMW9RWRqXNZKoT5AcS3LHyju4c2XEeRaVaiYjskQqXFtbG21tbVOvn3766RnlU7Sg\n4pw7Yma9Zvagc+4H/nA7wQrI6TkpDZkbg8XNR65Vu6iW7q900/N3PXzpM1/ix//fjyOnra+t576W\nKq3ViEhsxe5T2QDsMLMGgsUotznnTvhzHUCCYA5KeoLjRn9sk5kdds69GiEfyaJ2US0r3UpWfW5V\nrHSnBzVEVUSiK2pQ8cOAH81x7vGM1xMEAeaamst0+YhUrVSK2mPH4AtfKL/5N1I1Im8nLFKWzp2D\nn/wk+G8185MmP7N1azAPp9o/DykZBRWZv9ITGX//9/VD6idN1ly5ohWRpaQ0T0Xmr9BExsmjR/lF\nd3cwqiuCosw5KSY/aXIymaQmzkrKIgWmoCLzV2giY00iwedaWiIPEy7anJNiSU+a7OlhxZe+xOSP\no4/w0yKUUkgKKjJ/ZUxkrPp5J7W1vL9yJWtXxRvhp0UopZDUpyLzmyYyipQVBRURESkYBRURESkY\nBRURESkYBRWRapaehV/Nc3ykoDT6S6Ra+Vn4tw0Ncf4//SfeefZZJpcsiZRUw5AlFwUVkWoV2rp4\nyalTrDl7FiIOR9YwZMlFzV8i1Sq9dXFNDWjrYikQBRUpD1oYsvjSs/A7O6G7W3N9pCDU/FVip0+f\n4fWfHIud7uSpX89BaUokvTDk0aPBmlWvvaal24tFWxdLgSmolNCF8zWc/tVtXL/oThYvmYyV9uLF\n6Gs7lb1ZLAwJFbg4pMg8VjFBxcwSzrlkqcsR1YXzNXz765/lVPIpvn3iI771/XdiB5aKMYuFIaEC\nF4cUmceKGlT8FsE7gEFgBdDrnDs4k2vNbD+wzr90QMscFr3gfjmymHdPXI9zNbx7cgHvjiym6Y5U\nqYtVGloYUqRiFLumsgc44JzrBjCzQTNb75x7M861PuC8D6wBDBh3zo0U5y0Uxs2rL3DTqoucSl7H\nTbd+xE2rL5S6SKWlhSHnlaFj8fsBQfNbqkHRgoqZNRPULDpDh3uBJ4BNMa/dCYwBZ+dTk1fY4iWT\nfOv777D7e/vp+saGsm36OvaLmf141NfWc1+LfjwqUirFsmPHaJlBrVLzWypfMWsqLQQ1irOhYwPA\n9hlcmwDWA4+Z2WFgw3wMLouXTFK34njZBhSA1IcpVjbF7wQ/Pagfj4rkZ+HfPzQEfX0aiizXKOY8\nlRZgNOPYKFAf91rn3APOuRXAWqAOODDdjVOXNfdhJlKXUpz86CSpS1Xa1yPXCs3CJ5mE4eFSl0jK\nzLye/OicexVoBRrM7P5c13W+/kUFlphSl1J0/rCTFyZeoPOHnQosEtAsfMmjmM1fA8CGjGPLgfHZ\nXOucmzCzg2Sv8QAw+Lc/48+PfIMbr7+ZlpY2WlvbYhW8Gg2NDZEcTzLJJMnxJMNjw9z923eXulhS\naulZ+C++SMdDD6npq4L09fXR19c363yKGVT6ufaHv4WgA34210IQcA7nunHTv7ybf/f571G7ULO0\no2pc1kiiPsHQ6BCJ+gQNy+L9RXry6DEGzse7pyYxzhOahV+R2traaGtrm3r99NNPzyifogUV59wR\nM+s1swedcz/wh9sJ+kUws3agwTm3Z5prO8wsATQ753p8ujXAoemGFHd//jUFlJhqF9XS/ZVuXvyb\nF3noKw9Ruyj6D0jNByluPHqMljs1iVGk2hR7nsoGYIeZNQCNwDbn3Al/roNgVNeeaa4d8UFkh5lt\nB3YDw865R6e7qQLKzNQuquWW626JHVA++9VO/sXPh+BAn0YHiVSZogYVP0Q4awBwzj0e5Vrn3BGg\nKc59z53L1m2T28WLVT4RcRYW/3yI63+RpCY8Ouhu9cWIVIuKWftrOj97941Y16dOj+Bwc1Saynbh\nc41c/EyCT/x8iAUaHSRh6a2Lv/AFrUJdwaoiqCz7nXidv2dPD85RSSrf5Cdreeevuxn5jy/y1Y0a\nHSSenzT5meFh+Mu/1PYGFWxez1OR8jT5yVr+eZVGB0mInzRZc+UKvP12sNWBVCQFFRGZe37S5OTC\nhXDHHcFmbFKRqqL5S0RKLD1psqeHFV/6EpM/jr7J3MlTp7j1llti31IrIpeGgoqIFEdtLe+vXMna\nVatiJTve309Lc3Ps22lF5NJQ85eIiBSMgorkVPNBipuGT1LzgRaTlBJJpVhx8mQwekzmBTV/SVbh\nmfEf/qiPd/66m8lPajSXFJH2bpmXFFQkq/DM+E+8M8z7//VvONOUiJT2wuRF+t89xL+6/PssXagf\nAZmhbHu3aHWGsqegIlmlZ8Yv+vkQl2+9iVvvvp2blyzOmy515QJ/cvRZjtef4tDrx3jt890KLDIz\n6b1bhrQ6w3yiPhXJKj0z/i/++wc48ey3mIwQUACOn/8lQxfeZdIcb59LcvScdgaUGUoPQ+7sVNPX\nPKKayjx1+vQZXv/JsVhpTp76dazrJz9Zy8nfWRk5oADctuRmGhffxPHUL7mjPsGdS/XXpczCTPdu\n0TpjJaOgMg9dOF/D6V/dxvWL7mTxksnI6S5ejD7hbKZqFyzmr+78Frv/nx/w53/0hJq+pPi0zlhJ\nKajMMxfO1/Dtr3+WU8mn+PaJj/jW99+JFFgWnk/R9KtfU/NBas5HcdUuWMwtl1YqoEhphNYZmzx6\nlF90d5O6/fZISTULf/YUVOaZX44s5t0T1+NcDe+eXMC7I4tpumP6MfwLz6d44OudrBsZ4qO3jmp4\nsFS2jA7+z7VE34FUs/BnT0Flnrl59QVuWnWRU8nruOnWj7hpdf4NxepHhqg7EQzNXHA8yeJ3hkm1\nzO3QzNNnTvP64ddjpTn57sk5Ko1UlXQH/4sv0vFQjO0X1A9TEEUNKmZWB+wABoEVQK9z7mDca+Pk\nU2kWL5nkW99/h93f20/XNzZEavoaX93IxKoES0cG+eBTN/LWb05z+R+i/eCf/c3M/nK7dPkSy25Z\nFvn61JULHD98inOXU2o2k9mL28GvfpiCKXZNZQ9wwDnXDWBmg2a23jn3Zsxr4+RTcZZyjt/jDZby\nZS6T/3+ay0tqOfD9bt76zuP8j1//U+6MMZrrHz68NJuiRhLMbfk2x3/7FIdeH9LcFim+bPu9ROlb\nOXcO3noL7rpLQcgrWlAxs2ZgHdAZOtwLPAFsinqtmW0H1kfJpxJN9Y8kBxk79sf8X11f59Inrs+b\n7sMPr+OH7y3jf2AptVwuQkmjC89tOXp2mP/a/zfctSja7H2YWbPZ28lk7DSVSp8FU/0wBwcH+b2b\nb+ad995j8kc/mjZJzfnzfHbLFq4/dYqaO++MV7up4GBUzJpKCzDunDsbOjYAbI95bS8wFjGfsrbw\nfIo7zpxh4fkUl5dE+8t8qn/EOVac/g3/zfkU47/9qWnTnL+4kMf+4gGS777Ln2x5gL969gC1S8on\nsHw8t+UUjbU30bz6dmoXRKtNzbTZTD+kH9NnwVQ/zP/xjW/Q/r3vsSZKs9k//ROcOgWXL8ev3Xzx\ni8H1cYLRPAlExQ4qoxnHRoH6GNcui5kPQKwfbYBFFy+y4t1fsvAzayKnixsgwjWOcyc6OfD97kjp\npvpHkoOc/52bOPepm/KmOfFePSd/VQcsYOhUHYMn6vnd29/Pmy51fiG/GbuD1PmFsYJQ3HS1Cxaz\nO/EU/+t/+Tv+tz/9I2oXXBftPlcu8K/e+jZDN/6SN14b5vUv7okUWN6buMDwBx/x3sQFPl0XLXi9\nN3GBg7++wO/HSDPTdMW+lz4Ln+ZyDUcXLOG9yzV8OkqCxkY+unUVjCS58OmbGIpQuwGoPXaM2956\niwVXrnDlraMcjzDkueb8eRr/3TdZcuok52+5laE//y6TS5ZEel+j/zzOyE+P8d/9my4+tTrSO5uV\n+Tj6y8VN8KVHHmLPv32MS9fnbyZadPEim//Dn3PjP/+Gib/7u0g/9jMJEOEaR93JJPUjw7x/R/4R\nWen+kde/9W/4b/+n/4Ur1+f/H2bVp8e59VMTJN+9QuMtEzStGs+bJnV+IX+y5QGOj6zj51vORa7d\nzCRd6vxCurb9EcdH/oSun0e/1z+N/orjqfdgwSRvnU3yj/+cpP3Td02b5r2JCzT+7b/m4uK3SPxN\nF/uaNrNy8YJp05y+cIWNg3u4tPwX/NXf/pChf/mfI/1QTd1r+bHI6WaSZtb3WvwWjX/7r8u3fOX6\nWVyu4Xf/8BM0XjSGrl/CT+tvjFS+X12u4eiy6/js6BXeWbaQlbfdyedWrpg2zenXD3H9qREWTsIn\nTo1w83unWfl7rfnv9aszXPfNP+VLoxf5+f+5Ew4NzX1gcc4V5QFsBs7kO5bvWoK+lEj5+HNODz30\n0EOP+I+Z/NYXs6bSz7VNVC0EfSRxrh2IkQ/OOYtdUhERmZGirVLsnDsC9JrZg6HD7cBjAGbWbmab\n812bLx8RESkd801ExbmZ2Q0EkxaHgEbgFefcD/y57UDCObcpwrU5z4mISOkUNahIeTCzhHMuWepy\nFJOZ1TnnJkpdjnIQ5bOoxu+I5Bfle1GRm3SZWZ2ZPWdmW8zsGTNrL3WZSsnM9pvZpH9cAepKXaZi\n8c2qgwSDOcLHq+47kuuz8Oeq5jtiZgkz6/fvtd/M1oTOVdX3YrrPwp+P/b2Yj0OKo6jqZVzC/Dpp\n7wNrACOYVDpS0kIViW8m7QeWZzldVd+R6T6LavqOhNYN/BrBe+0G+s1smZ9QXTXfi3yfxUy/FxVX\nUwkt8bIvdDi9jEs12gFMAGedc29W6o9FNs65s76Z56rJstX4Hcn1WXjV9B1ZSzDg56c+ULQT/GC2\n+r/S11M934ucn4U/P6PvRcUFFXIv8bK2ROUptQaC5o4hX72NvqhW5dJ35GpV8x1xzvWEfxx9oHUE\ntbh7yb4EVEV+L6b5LA75QzP6XlRqUIm1jEslc8494JxbQfA/Rh1woMRFKgf6joRU83fEzNYSbJ1x\nFmimir8Xoc/iHMz8e1GJQUWycM69SlCtbTCz+0tdHik/Vfod2QxsKHUhykTWzyLu96ISg8oA13ZG\nLgfyL3hV4Xz19iBV8pfXNPQdyaGaviNmthV4JtTcVbXfiyyfxVXifC8qMajEWQ6mGi0HDpe6ECWm\n78j0Kv47YmbrCCZNh0d1VeX3IsdnkU2k70XFBRUt4/IxPwZ9Xej1GuBQhY/uyUvfkY9V43fE9x0s\nByb8+282s+3V+L3I9VnM5ntRkTPqtYxLwH8R9hOM6NgNDFfT5+DH2W8EniP4C2ubbx+uuu9Irs+i\n2r4jfjLjKwTvN+wx59x3/ee0nSr4Xkz3WQCvMsPvRUUGFRERKY2Ka/4SEZHSUVAREZGCUVAREZGC\nUVAREZGCUVAREZGCUVAREZGCUVARmYafBFaxG1aJFJqCipQ9v2PhqN99bnXGuc3+3Mtmdk+B77uV\nYBJcSyHzjXH/zf4xaGadGefWhT6TvWa2y+9YeMDMrvjJnSJFp8mPMi+Y2WaC2eDDQEt44Tsz2wXs\nnovd+cxsElibnolfLH6m+z7n3G0+WC7PLIOZPQdsds4tyDjeSbCE+cgclm1ZsT8TmR9UU5H5ZCfB\nxkH7M45PUHkrybbi35PfdS/bD3jW9+yc657DgJIgWK1WJCsFFZlPXgG6gA4z25LtAt9cNJleFNDM\ntvrX9/hmtAP+2FbfrPSymdX5RfT6/SOz6ajRpxs1s32he63x6fb586v98XX+9Wb/3+M5yprw6bf6\nPLam0xPsa9FgZlvi7G3i3+P9oeeTZvZy6Pz29Otpyt8eKv9z/n3v9Vk0E6zk25VukvN5dPprd0Ut\nq1Qo55weepT9g2ADofv98+eAydDr7cDq0LWTwIMZr+/xz/uBl0PnRrO87sxIuxe4AVjjXz+D3wkv\ndN2BbPkCq8NlCZ2v89cszUizxT/fSrAq7HSfyfZQ+fb6MlxJfy7+ml3hfHy+q2OUP/y+7wl9Jn/g\nnyeA/lC6LaX+ruhR2sfC6OFHpDw45x4xs1Zgv/9vHKNcvSdEf5bXjRlpdrugD+eImT0PrCfo23G+\nxmQEHfoNGffZ74JmqJEs5dgInHF+61bveYKa2LMx3o9zzm1Kv8js0CdoMhwys9W+LA3OuRHfR5Wv\n/K+E3jf+fLrfyvzNk6Gl4x93zsUpu1QgBRWZr9oJftj3EzSLFcsAQUBoJlgOfKY/opmBC679YZ+J\nfYR2L/Q/+q8Aj/gmrPRn1UK88k/XZ9UMvOoD1VoX7EsiVUp9KjIvuWB707UEP2jF3kTpEMGPbEf4\noO+biTqUN1cAmdWOi865s74m0h46vBN4GNjoPt4TY5zg85sSs/xTaVww6m41QUAbmHHhpSIoqMh8\n0UjGX/f+L+INWa4dxv/gm9nDBBsNTVcDmG7f7cwx9xuAbQQ/oA2+czrhf8S3uxx7fF+TqXN7gPGM\n5qoO4N/757+Vp1wAy7IdNLNmgia69L2yjdbaSzAAIWr5w2UZJnjvdf6/m30wexQY0xyZKlfqTh09\n9Mj3IGjqGiSoIdyT5fwzXN1Rvw44AxwH7vfpnvH5nPGvVxN0QI/66+4JvT6Uzg/YQtBh/Yx/hDvB\nO33aK/jO/ND9r+Qqbyj96lDeW/Ad+r4c/T6Pzhxp1/nP5IrPYxfBAIZ0Z/3XMq7vDH9GecrfHsp3\nNcEgialrCAYInPFlXuOfdxIMAshaXj2q56HJjyIiUjBq/hIRkYJRUBERkYJRUBERkYJRUBERkYJR\nUBERkYJRUBERkYJRUBERkYJRUBERkYJRUBERkYL5/wGtDmYHaetI2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x124418860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "# generate samples for different values of mu\n",
    "kpts=np.arange(0,25)\n",
    "for mu, c in zip([1,3,12], sns.color_palette()[:4]):\n",
    "    randomVariates = poisson.rvs(mu, size=1000)\n",
    "    plt.hist(randomVariates, normed=True, color=c, alpha=0.2, bins=range(0,26), label='$\\mu=' + np.str(mu) + '$')\n",
    "    plt.plot(kpts, poisson.pmf(kpts, mu), '.', color=c)\n",
    "\n",
    "plt.legend()\n",
    "plt.title(\"Poisson Distribution\")\n",
    "plt.xlabel(\"Number of Events\")\n",
    "plt.ylabel(\"Normed Counts\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximizing the lnL \n",
    "\n",
    "For the Poisson distribution we can solve the for the  maximum liklehood analytically\n",
    "\n",
    "$$\\frac{\\partial \\; L(\\mu)}{\\partial \\; \\mu} = \\frac{\\partial \\; }{\\partial \\; \\mu} \\bigg( \\sum_{k=1}^{n} -\\mu + x_k \\; ln(\\mu)\\bigg)$$\n",
    "\n",
    "$$0 = \\sum_{k=1}^{n} \\bigg( -1 + \\frac{x_k}{\\mu} \\bigg)$$\n",
    "$$\\mu = \\frac{\\sum_{k=1}^{n} x_k }{n}$$\n",
    "\n",
    "For many likelihoods we cannot solve analytically and end up with numerical solutions (see the upcoming lectures on MCMC and robust statistics that account for outliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "Today we heard about:\n",
    "\n",
    "* Classical vs. Bayesian Inference (more next week)\n",
    "* Likelihood function\n",
    "* Maximum Likelihood Estimation\n",
    "  * Homoscedastic Samples\n",
    "  * Heteroscedastic Samples\n",
    "* Assessing the goodness of Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Truncated/Censored Data and Other Cost Functions\n",
    "\n",
    "Note that dealing with missing data points (\"censored data\") adds complications that we don't have time to get into here, but see Textbook 4.2.7  Also maximum likelihood is just one possible \"cost function\", see Textbook, 4.2.8. "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected",
   "theme": "sky"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}