ee046202_tutorial_01_classic_statistics_point_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <img src=\"https://img.icons8.com/dusk/64/000000/artificial-intelligence.png\" style=\"height:50px;display:inline\"> EE 046202 - Technion - Unsupervised Learning & Data Analysis\n",
    "---\n",
    "\n",
    "#### <a href=\"https://taldatech.github.io\">Tal Daniel</a>\n",
    "\n",
    "## Tutorial 01 - Classical Methods in Statistical Inference - Point Estimation\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/checklist.png\" style=\"height:50px;display:inline\"> Agenda\n",
    "---\n",
    "* [Statistical Inference - Classical Methods](#-Statistical-Inference---Classical-Methods)\n",
    "* [Point Estimation](#-Point-Estimation)\n",
    "* [Evaluating Estimators](#-Evaluating-Estimators)\n",
    "* [Point Estimation of Mean and Variance](#-Point-Estimation-of-Mean-and-Variance)\n",
    "* [Non-Parametric Point Estimation Using The Tail Sum Formula](#Non-Parametric-Point-Estimation-Using-The-Tail-Sum-Formula)\n",
    "* [Maximum Likelihood Estimator (MLE)](#-Maximum-Likelihood-Estimation-(MLE))\n",
    "    * [KL Divergence and Asymptotic Consequences](#-Asymptotic-Properties-of-MLEs)\n",
    "* [Recommended Videos](#-Recommended-Videos)\n",
    "* [Credits](#-Credits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# imports for the tutorial\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <img src=\"https://img.icons8.com/doodle/48/000000/statistics.png\" style=\"height:50px;display:inline\"> Statistical Inference - Classical Methods\n",
    "---\n",
    "* **Statistical Inference** - a collection of methods and tools to draw conclusions from data that are usually affected by randomness.\n",
    "    * General Setup - there is an unknown quantity that we wish to estimate by observing given data. There are 2 approaches:\n",
    "        * Classical Inference (**Frequentist**) - the unknown quantity $\\theta$ is **deterministic**. We estimate non-random quantities.\n",
    "        * Bayesian Inference - we assume the unknown quantity $\\theta$ is a **random variable** and we make assumptions on the type of distribution. After observing the data, we can update the distribution using Bayes' rule.  We estimate random variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Examples:\n",
    "    * Predicting the results of an election - we cannot possibly poll the entire population, thus, we pick a random sample from the population to get \"where the wind blows\". Here, the randomness comes from the sampling process. Another source of randomness may be the time frame in which the poll was conducted (one month or one week before the elections).\n",
    "        * In the *classical* approach: $\\theta$ is the percentage of people that vote for candidate A. After polling $n$ randomly chosen voters, where $n_A$ voters said they would vote for candidate A, we may estimate $\\theta$ by: $$ \\hat{\\theta} = \\frac{n_A}{n} $$\n",
    "            * Notice that $\\hat{\\theta}$ is a **random variable** as it depends on the *random* sample. \n",
    "    * Receiver-Transmitter - the receiver may get a corrupted version of messages due to random noise.\n",
    "        * In the *Bayesian* approach, $\\theta \\sim Bernoulli(p)$ is the transmitted bit ({0,1}), where $p$ is the probability to transmit 1. The receiver has to recover $\\theta$ based on the knowledge of the distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <img src=\"https://img.icons8.com/external-wanicon-lineal-color-wanicon/64/000000/external-weight-supermarket-wanicon-lineal-color-wanicon.png\" style=\"height:50px;display:inline\"> Parametric vs. Non-Parametric Estimation\n",
    "---\n",
    "* Parametric models: the model structure (i.e., distribution) is specified a priori. For example, we assume the data is Normal distributed and we try to find its parameters.\n",
    "    * Typical models: Linear/Logistic Regression.\n",
    "* Non-parametric models: the model structure is not specified a priori but is instead **determined from data**. These models do not lack parameters but that the number the parameters is flexible and not fixed in advance.\n",
    "    * Typical models: histogram, kernel denisity estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./assets/param_vs_nonparam.png\" style=\"height:400px\">\n",
    "\n",
    "* Note that non-linear SVM (listed in the table) is a non-parametric method, whereas linear SVM (not listed in the table) is a parametric method as it fits linear classification model (linear classifier).\n",
    "* <a href=\"https://stats.stackexchange.com/questions/268638/what-exactly-is-the-difference-between-a-parametric-and-non-parametric-model\">Source 1</a>, <a href=\"https://probml.github.io/pml-book/book0.html\">Source 2</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/cotton/64/000000/pickup-point.png\" style=\"height:50px;display:inline\"> Point Estimation\n",
    "---\n",
    "* **Assumption** - $\\theta$ is a **fixed**, non-random, quantity.\n",
    "    * Example: $\\theta$ can be the expected value of a random variable $\\theta = \\mathbb{E}[x]$\n",
    "* The data - a random sample $\\{X_i\\}_{i=1}^n$ such that $X_i$'s have the **same distribution** as $X$.\n",
    "* The point estimator - a function of the random sample: $$ \\hat{\\theta} = h(X_1, X_2, ..., X_n) $$\n",
    "    * Note: the estimator may depend *stochastically* on the data, but we will assume the dependence is deterministic.\n",
    "    * Example: if $\\theta = \\mathbb{E}[x]$ then we may choose $\\hat{\\theta}$ to be: $$ \\hat{\\theta} = \\overline{X} = \\frac{X_1 + X_2 + ... + X_n}{n} $$\n",
    "        * Note: this not always the 'best' estimator for the mean\n",
    "* There are many possible estimators for $\\theta$, so how can we make sure we have chosen a good estimator? We need to define ways to evaluate our estimators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/cotton/64/000000/rating.png\" style=\"height:30px;display:inline\"> Evaluating Estimators\n",
    "---\n",
    "* **<a style=\"color:red\"> Bias </a>** - the bias of an estimator $B(\\hat{\\theta})$ is a measure of how far is the estimator $\\hat{\\theta}$ from the real $\\theta$ **on average**. A *scalar* (not R.V). Formal definition: $$ B(\\hat{\\theta}) = \\mathbb{E}[\\hat{\\theta}] - \\theta $$\n",
    "    * Note that the bias cannot be actually computed (why is that?).\n",
    "    * Recall that $\\hat{\\theta}$ is a **random variable**.\n",
    "    * **Unbiased** estimator - we would like the bias to be close to 0, which indicates that on average $\\hat{\\theta}$ is close to $\\theta$. $\\hat{\\theta}$ is an *unbiased* estimator of $\\theta$ if: $B(\\hat{\\theta}) = 0 \\rightarrow \\mathbb{E}[\\hat{\\theta}] = \\theta$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - Bias of an Estimator\n",
    "---\n",
    "Let $X_1, X_2, ..., X_n$ be a random sample. Assume (always) that the samples are independent and identically distributed (iid).\n",
    "1. Show that the sample mean: $\\hat{\\theta}$ to be: $ \\hat{\\theta} = \\overline{X} = \\frac{X_1 + X_2 + ... + X_n}{n}$ is an *unbiased* estimator of $\\theta=\\mathbb{E}[X_i]$\n",
    "2. If we choose $ \\hat{\\theta}_1 = X_1 $, is the estimator unbiased? Is it a good estimator?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "1. $B(\\hat{\\theta}) = \\mathbb{E}[\\hat{\\theta}] - \\theta = \\mathbb{E}[\\overline{X}] - \\theta =\\mathbb{E}[X_i] - \\theta = 0 $\n",
    "2. If we choose $ \\hat{\\theta}_1 = X_1 $ then $ B(\\hat{\\theta}_1) = \\mathbb{E}[\\hat{\\theta}_1] - \\theta = \\mathbb{E}[X_1] - \\theta = 0$. This is an *unbiased* estimation! But is it good, or is it as good as $\\hat{\\theta} = \\overline{X}$? Not necessarily, as there can be many samples with different values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* **<a style=\"color:red\"> MSE </a>** (Mean Squared Error) - the MSE of an estimator: $$ MSE(\\hat{\\theta}) = \\mathbb{E}\\big[(\\hat{\\theta} - \\theta)^2 \\big] $$\n",
    "    * Note that the expression $\\hat{\\theta} - \\theta$ is the *error* we make by estimating $\\theta$ with $\\hat{\\theta}$.\n",
    "    * The MSE is a measure of the expected (squared) error. **Smaller** MSE is generally an indication of a better estimator.\n",
    "    * We define the **variance** of the estimator as follows: $Var(\\hat{\\theta}) = \\mathbb{E}\\big[(\\mathbb{E}[\\hat{\\theta}] - \\hat{\\theta})^2 \\big]$\n",
    "    * It holds: (HW) $$ MSE(\\hat{\\theta}) = Var(\\hat{\\theta}) + Bias^2(\\hat{\\theta}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - MSE of an Estimator\n",
    "---\n",
    "Let $X_1, X_2, ..., X_n$ be a random sample with mean $\\mathbb{E}[X_i] = \\theta$ and variance $Var(X_i) = \\sigma^2$\n",
    " and consider the following estimators:\n",
    " 1. $\\hat{\\theta}_1 = X_1$\n",
    " 2. $\\hat{\\theta}_2 = \\overline{X} = \\frac{X_1 + X_2 + ... + X_n}{n}$\n",
    " \n",
    "Find $MSE(\\hat{\\theta}_1), MSE(\\hat{\\theta}_2)$ and show that for $n>1$ we have $MSE(\\hat{\\theta}_1) > MSE(\\hat{\\theta}_2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "$$ MSE(\\hat{\\theta}_1) = \\mathbb{E}\\big[(\\hat{\\theta}_1 - \\theta)^2 \\big] = \\mathbb{E}\\big[(X_1 - \\mathbb{E}[X_1])^2 \\big] = Var(X_1) = \\sigma^2 $$ \n",
    "\n",
    "* We use the following:\n",
    "    * $\\mathbb{E}[X^2] = Var(X) + (\\mathbb{E}[X])^2 $\n",
    "    * For a constant $b$: $Var(X + b) = Var(X)$\n",
    "    * For a constant $b$: $Var(bX) = b^2Var(X)$\n",
    "    * $\\mathbb{E}[\\overline{X} - \\theta] = \\frac{n\\mathbb{E}[X_i]}{n} - \\mathbb{E}[X_i] = 0$\n",
    "    * $Var[\\overline{X}] = Var[\\frac{1}{n}\\sum_{i=1}^nX_i] = \\frac{1}{n^2}n\\sigma^2 = \\frac{\\sigma^2}{n}$\n",
    "        * Recall that the samples are i.i.d., and the variance of the sum of i.i.d. samples is the sum of the variances.\n",
    "\n",
    "$$ MSE(\\hat{\\theta}_2) = \\mathbb{E}\\big[(\\hat{\\theta}_2 - \\theta)^2 \\big] = \\mathbb{E}\\big[(\\overline{X} - \\theta)^2 \\big] = Var(\\overline{X} - \\theta) + (\\mathbb{E}[\\overline{X} - \\theta])^2 = Var(\\overline{X}) + 0 = \\frac{n\\sigma^2}{n^2} = \\frac{\\sigma^2}{n}$$ \n",
    "        \n",
    "\n",
    "Thus, for $n>1$: $MSE(\\hat{\\theta}_1) > MSE(\\hat{\\theta}_2) $, which means that $\\overline{X}$ is a better estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* **<a style=\"color:red\"> Consistency </a>** - an estimator is consistent if as the sample size $n$ grows, then $\\hat{\\theta}$ converges to the real value of $\\theta$.\n",
    "    * Formally: Let $\\hat{\\theta}_1, \\hat{\\theta}_2, ..., \\hat{\\theta}_n$ be a sequence of point estimators of $\\theta$. We say that $\\hat{\\theta}_n$ is a **consistent** estimator of $\\theta$, if: $$ \\lim_{n \\to \\infty} P(|\\hat{\\theta}_n - \\theta| \\geq \\epsilon) = 0, \\forall \\epsilon >0 $$\n",
    "    * Note: other convergence types yield different definitions of consistency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - Consistency of an Estimator\n",
    "---\n",
    "Let $X_1, X_2, ..., X_n$ be a random sample with mean $\\mathbb{E}[X_i] = \\theta$ and variance $Var(X_i) = \\sigma^2$.\n",
    "\n",
    "Show that $\\hat{\\theta}_n = \\overline{X}$ is a *consistent* estimator of $\\theta$.\n",
    "\n",
    "Reminder:\n",
    "* **Chebyshev's inequality**: let $\\mu = \\mathbb{E}[X], \\sigma^2= Var[X]$. Then: $ P(|X- \\mu| > t) \\leq \\frac{\\sigma^2}{t^2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "$$ P(|\\hat{\\theta}_n - \\theta| \\geq \\epsilon) = P(|\\overline{X} - \\theta| \\geq \\epsilon) $$\n",
    "\n",
    "* Using **Chebyshev's inequality** and recall that $\\overline{X}$ is a random variable with **mean** $\\theta$:\n",
    "\n",
    "$$ \\rightarrow P(|\\overline{X} - \\theta| \\geq \\epsilon) \\leq \\frac{Var[{\\overline{X}]}}{\\epsilon^2} = \\frac{\\sigma^2}{n} \\cdot \\frac{1}{\\epsilon^2} $$\n",
    "$$ \\rightarrow \\lim_{n \\to \\infty} \\frac{\\sigma^2}{n} \\cdot \\frac{1}{\\epsilon^2} = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/cute-clipart/64/000000/75-percents.png\" style=\"height:50px;display:inline\"> Point Estimation of Mean and Variance \n",
    "---\n",
    "* The sample mean, $\\overline{X}$, is often a reasonable point estimator for the mean. But what about the variance? \n",
    "* Before, we assumed that the variance was known, but when we want to estimate it, we need to take a similar approach.\n",
    "\n",
    "By definition, the variance of a distribution $\\sigma^2$ is: $$ \\sigma^2 = \\mathbb{E}[(X - \\mu)^2] .$$\n",
    "* If we define the following *random variable*: $Y = (X -\\mu)^2$ then the number $\\sigma^2$ is the mean of that variable, that is $\\sigma^2 = \\mathbb{E}[Y]$.\n",
    "* But wait, if $\\sigma^2$ is the mean of $Y$, we have already derived a point estimator for the mean! $$ \\hat{\\sigma}^2 = \\hat{Y} = \\frac{1}{n}\\sum_{i=1}^nY_i = \\frac{1}{n}\\sum_{i=1}^n (X_i -\\mu)^2 $$\n",
    "    * By the **linearity** of the expectation, this is an **unbiased** estimator of the variance.\n",
    "    * From the **weak law of large numbers** this is also a **consistent** estimator of the variance.\n",
    "* The problem? What if we do not know the value of $\\mu$? It is often reasonable to replace $\\mu$ with our point estimate of $\\mu$, which transforms the estimate for $\\sigma^2$ to: $$ \\overline{S}^2 = \\frac{1}{n} \\sum_{i=1}^n (X_i - \\overline{X})^2 =... =\\frac{1}{n}\\big( \\sum_{i=1}^nX_i^2 -n\\overline{X}^2\\big)  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\overline{S}^2 = \\frac{1}{n} \\sum_{i=1}^n (X_i - \\overline{X})^2 = \\frac{1}{n}\\sum_{i=1}^n(X_i^2 -2X_i \\overline{X} +\\overline{X}^2)$$\n",
    "$$= \\frac{1}{n}(\\sum_{i=1}^n X_i^2 -2\\overline{X}\\sum_{i=1}^n X_i +\\sum_{i=1}^n \\overline{X}^2)$$\n",
    "$$= \\frac{1}{n}(\\sum_{i=1}^n X_i^2 -2n\\overline{X}^2 +n\\overline{X}^2)$$\n",
    "$$=\\frac{1}{n}\\big( \\sum_{i=1}^nX_i^2 -n\\overline{X}^2\\big)  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - Bias of the Variance Estimator\n",
    "---\n",
    "Let $X_1, X_2, ..., X_n$ be a random sample with mean $\\mathbb{E}[X_i] = \\mu$ and variance $Var(X_i) = \\sigma^2$.\n",
    "\n",
    "Suppose that we use: $$ \\overline{S}^2 = \\frac{1}{n}\\big( \\sum_{i=1}^nX_i^2 -n\\overline{X}^2\\big) $$ as our estimate for $\\sigma^2$. Find the bias of the estimator: $$ B[\\overline{S}^2] = \\mathbb{E}[\\overline{S}^2] -\\sigma^2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "* $\\mathbb{E}[\\overline{X}^2] = (\\mathbb{E}[\\overline{X}])^2 + Var(\\overline{X}) = \\mu^2 +\\frac{\\sigma^2}{n}$\n",
    "* $\\mathbb{E}[X^2] = (\\mathbb{E}[X])^2 + Var(X) = \\mu^2 +\\sigma^2$\n",
    "$$ \\mathbb{E}[\\overline{S}^2] = \\frac{1}{n}\\big(\\sum_{i=1}^n\\mathbb{E}[X_i^2] -n \\mathbb{E}[\\overline{X}^2] \\big) $$ $$ =\\frac{1}{n} \\big( n(\\mu^2 +\\sigma^2) -n(\\mu^2 +\\frac{\\sigma^2}{n})\\big) = \\frac{n-1}{n} \\sigma^2 $$ $$ \\rightarrow B[\\overline{S}^2] = \\mathbb{E}[\\overline{S}^2] -\\sigma^2 = -\\frac{\\sigma^2}{n} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Thus, $ \\overline{S}^2$ is a **biased estimator** of the variance!\n",
    "* If $n$ is very *large*, then the bias is very **small**.\n",
    "* How can we obtain an **unbiased estimator** of the variance?\n",
    "    * By simply multiplying $\\overline{S}^2$ by $\\frac{n}{n-1}$.\n",
    "\n",
    "In conclusion, we define the *unbiased* estimator of the variance, called the **sample variance** to be: $$ S^2 = \\frac{1}{n-1} \\sum_{i=1}^n (X_i - \\overline{X})^2 =... =\\frac{1}{n-1}\\big( \\sum_{i=1}^nX_i^2 -n\\overline{X}^2\\big) $$\n",
    "* The **sample standard deviation** is $S = \\sqrt{S^2}$, which is an **unbiased estimator** of the standard deviation.\n",
    "* Note: a good estimator should be **asymptotically unbiased**. In many cases, for a finite sample, a biased estimator is *better* than an unbiased one! (recall overfitting from ML?). There is nothing special in the lack of bias, *except* asymptotically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/ultraviolet/80/000000/tail-of-whale.png\" style=\"height:50px;display:inline\">Non-Parametric Point Estimation Using The Tail Sum Formula\n",
    "---\n",
    "* Every random variable with possible values $\\{0,1,...,n\\}$ is a **counting variable** representing number of events that occur in some list of $n$ events $A_1,..., A_n$.\n",
    "* To see this, let $A_j$ be the event that $(X \\geq j)$. If $X=x$ for $0 \\leq x \\leq n$, then $A_j$ occurs for $1 \\leq j \\leq x$ and $A_j$ **does not** occur for $x<j\\leq n$. So if $X =x$, the number of events $A_j$ that occur is exactly $x$.\n",
    "* The resulting formula for $\\mathbb{E}[X]$ - **The Tail Sum Formula for Expectation**: $$ \\text{For } X \\text{ with possible values } \\{0,1...,n\\}, $$ $$ \\mathbb{E}[X] = \\sum_{j=1}^n P(X\\geq j) $$\n",
    "* Proof outline:\n",
    "    * Define $p_j = P(X=j)$.\n",
    "    * The expectation $\\mathbb{E}[X] = 1p_1 + 2p_2 + 3p_3 +... + np_n$ is the following sum: $$ p_1 $$ $$ +p_2 +p_2 $$ $$ +p_3 +p_3 + p_3 $$ $$ +...+... +... $$ $$ +p_n +p_n +... +p_n $$\n",
    "    * By the addition rule of probabilities, and the assumption that the only possible values of $X$ are $\\{0,1,...,n\\}$, the sum of the first column of $p$ is $P(X\\geq1)$, the sum of the second solumn is $P(X\\geq 2)$ and so on. The sum of the $j^{th}$ column is $P(X\\geq j), 1 \\leq j \\leq n$. The whole sum is the sum of the column sums: $$ \\sum_{j=1}^n P(X\\geq j) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/dusk/64/000000/maximum-value.png\" style=\"height:50px;display:inline\"> Maximum Likelihood Estimation (MLE)\n",
    "---\n",
    "* The MLE is an estimator that picks the best parameters by maximizing the **likelihood** of the distribution. Recall that when we developed Bayes rule, the likelihood is $p(D|\\theta)$ (which is a function of $\\theta$).\n",
    "    * MLE can be applied to both parametric and non-parameteric models. A parametric example would be estimating the parameters of a distribution (e.g. $\\mu, \\sigma$ of a Gaussian). A non-parametric example would be estimating the density  (like KDE-Kernel Density Estimation), where you don't impose any specific distribution on the data. \n",
    "* Definition: $$\\hat{\\theta}_{MLE} = \\underset{\\theta \\in \\mathcal{R}^{p}}{\\mathrm{argmax}} p(D|\\theta) =  \\underset{\\theta \\in \\mathcal{R}^{p}}{\\mathrm{argmax}} \\log p(D|\\theta)$$\n",
    "    * The last equality is true since the log function is **monotonically increasing**. Therefore if a function $f(x) \\geq 0$, achieves a maximum at $x_1$, then $\\log⁡(f(x))$ also achieves a maximum at $x_1$.\n",
    " \n",
    "* We assume the variables are **I.I.D (independent identically distributed)**. Note that these are the samples.\n",
    "* $L(\\theta) = p(D|\\theta) = p(x_1, x_2, ..., x_n|\\theta) = \\prod_{k=1}^n p(x_k|\\theta)$\n",
    "* $l(\\theta) = \\log \\big(L(\\theta)\\big) = \\sum_{k=1}^{n} \\log p(x_k|\\theta)$\n",
    "    * $\\log(x \\cdot y \\cdot z) = \\log x + \\log y + \\log z$\n",
    "* $\\rightarrow \\hat{\\theta}_{MLE} =  \\underset{\\theta \\in \\mathcal{R}^{p}}{\\mathrm{argmax}} \\{ l(\\theta)\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:50px;display:inline\"> Exercise  - MLE for Univariate Gaussian\n",
    "---\n",
    "Given $\\{x_i\\}_{i=1}^n$ i.i.d samples of $X \\sim N(\\mu, \\sigma^2)$, what is the MLE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "The first thing to ask yourself is, **what are the parameters** in this problem? In our case, the parameters are $\\theta = [\\mu, \\sigma^2]$, it is just a matter of notation.\n",
    "\n",
    "* $p(x_i) = \\frac{1}{\\sqrt{2\\pi \\sigma^2}} e^{- \\frac{1}{2} \\frac{(x_i - \\mu)^2}{\\sigma^2}}$\n",
    "* $L(\\theta) = L(\\mu, \\sigma^2) = p(x_1, x_2, ..., x_n |\\mu, \\sigma^2) = \\prod_{i=1}^n p(x_i|\\theta) = \\frac{1}{(2\\pi \\sigma^2)^{\\frac{n}{2}}} e^{\\frac{-1}{2 \\sigma^2} \\sum_{i=1}^n (x_i - \\mu)^2}$\n",
    "* $l(\\theta) = \\log L(\\theta) = -n (\\log \\pi + \\frac{1}{2} \\log \\sigma^2) - \\frac{1}{2 \\sigma^2} \\sum_{i=1}^n (x_i - \\mu)^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Find the optimal $\\theta$\n",
    "As usual, find the point where the deriviative w.r.t $\\theta$ is 0\n",
    "* $\\frac{\\partial l}{\\partial \\mu} = \\frac{1}{\\sigma^2} \\sum_{i=1}^n (x_i - \\mu) = 0 \\rightarrow \\hat{\\mu}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n x_i$\n",
    "* $\\frac{\\partial l}{\\partial \\sigma^2} = - \\frac{n}{2 \\sigma^2} + \\frac{1}{2\\sigma^4}\\sum_{i=1}^n (x_i - \\mu)^2 = 0 $\n",
    "    * Plug in $\\mu = \\hat{\\mu}_{MLE} \\rightarrow \\hat{\\sigma^2}_{MLE} = \\frac{1}{n}\\sum_{i=1}^n (x_i - \\hat{\\mu}_{MLE})^2 $\n",
    "* Summary: $$ \\hat{\\mu}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n x_i$$ $$\\hat{\\sigma^2}_{MLE} = \\frac{1}{n}\\sum_{i=1}^n (x_i - \\hat{\\mu}_{MLE})^2 $$\n",
    "\n",
    "* Do these look familiar? These are the **empirical** mean and variance!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def plot_normal_mle():\n",
    "    mu_real = 5\n",
    "    var_real = 36\n",
    "    num_samples= 1000\n",
    "    samples = np.random.normal(mu_real, np.sqrt(var_real), size=(num_samples))\n",
    "\n",
    "    mu_mle = np.sum(samples) / num_samples\n",
    "    var_mle = np.sum(np.square(samples - mu_mle)) / num_samples\n",
    "\n",
    "    print(\"mu mle: \", mu_mle)\n",
    "    print(\"var mle: \", var_mle)\n",
    "\n",
    "    x = np.linspace(-30, 30, 10000)\n",
    "    f_x_mle = (1 / np.sqrt(2 * np.pi * var_mle)) * np.exp(-0.5 * (np.square(x - mu_mle)) / var_mle)\n",
    "\n",
    "    # set bins for histogram\n",
    "    n_bins = 100\n",
    "    bins_edges = np.linspace(samples.min(), samples.max() + 1e-9, n_bins + 1)\n",
    "\n",
    "    fig = plt.figure(figsize=(8, 8))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.grid()\n",
    "    ax.set_ylabel('Histogram (PDF)')\n",
    "    ax.set_xlabel('x')\n",
    "\n",
    "    # plot histogram\n",
    "    ax.hist(samples, bins=bins_edges, density=True)\n",
    "    # plot estimation\n",
    "    ax.plot(x, f_x_mle, linewidth=3, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu mle:  4.6470617816063315\n",
      "var mle:  35.021587256422\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAHgCAYAAABJrX+JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwcdZ3/8ddnZjK5JyEHue+DEHKThEsgIEjCYRCQhRUEFgVWQFZBjPtT8VzZdXE9VgUWEVEU8EACQTkkw5WE3Ce5D8h9Qe5jMjPf3x/V0109mcz0JN397a5+Px+PflDfquruT4pM3lPfqvp+zTmHiIiIREuR7wJEREQk/RTwIiIiEaSAFxERiSAFvIiISAQp4EVERCJIAS8iIhJBJb4LSKcOHTq43r17p+3z9u/fT8uWLdP2eflOxyNBxyKZjkeCjkUyHY9k6T4ec+bM2eGc61jXtkgFfO/evZk9e3baPq+8vJxx48al7fPynY5Hgo5FMh2PBB2LZDoeydJ9PMzs/WNtUxe9iIhIBCngRUREIkgBLyIiEkEKeBERkQhSwIuIiESQAl5ERCSCFPAiIiIRpIAXERGJIAW8iIhIBCngRUREIkgBLyIiEkEKeBERkQhSwIuIiESQAl5ERCSCFPAiIiIRpIAXERGJIAW8iIhIBCngRUREIkgBLyIiEkElvgsQEcmk3pOmxJfXPXiZx0pEsktn8CIiIhGkgBcREYkgBbyIiEgEKeBFREQiSAEvIiISQQp4ERGRCFLAi4iIRJACXkREJIIU8CIiIhGkgBcREYkgBbyIiEgEKeBFREQiSAEvIiISQQp4ERGRCFLAi4iIRJACXkREJIIU8CIiIhGkgBcREYkgBbyIiEgEKeBFREQiSAEvIiISQQp4ERGRCFLAi4iIRJACXkREJIIyGvBmNt7MlpvZKjObVMd2M7OfxrYvNLNRoW1fMrMlZrbYzP5gZs0yWauIiEiUZCzgzawY+DkwARgMXG9mg2vtNgEYEHvdBvwy9t5uwBeB0c65IUAxcF2mahWR/NZ70hR6T5riuwyRnJLJM/ixwCrn3BrnXAXwNDCx1j4TgSddYAbQ1sy6xLaVAM3NrARoAWzKYK0iIiKRYs65zHyw2TXAeOfc52LtG4EznHN3hfZ5EXjQOfd2rP0P4KvOudlmdg/wfeAg8Ipz7jPH+J7bCM7+6dSp0+lPP/102v4M+/bto1WrVmn7vHyn45GgY5HM9/FYtHE3AEO7tTnmtmNtTzffxyLX6HgkS/fxuOCCC+Y450bXta0kbd9yNKtjXe3fJurcx8xOIji77wPsAv5oZjc453531M7OPQo8CjB69Gg3bty4Eyo6rLy8nHR+Xr7T8UjQsUjm+3jcHOueX/eZo2u4OdR1X9f2dPN9LHKNjkeybB6PTHbRbwB6hNrdObqb/Vj7XASsdc5td84dAf4CnJ3BWkVERCIlkwE/CxhgZn3MrJTgJrnJtfaZDHw2djf9mcBu59xm4APgTDNrYWYGfBxYmsFaRUREIiVjXfTOuUozuwt4meAu+Medc0vM7I7Y9oeBl4BLgVXAAeCW2LZ3zexPwFygEphHrBteREREGpbJa/A4514iCPHwuodDyw648xjvfQB4IJP1iYiIRJVGshMREYkgBbyIiEgEKeBFREQiSAEvIiISQRm9yU5EJBXhceTXPXhZo7eLyNF0Bi8iIhJBCngREZEIUsCLiIhEkAJeREQkghTwIiIiEaSAFxERiSAFvIiISAQp4EVERCJIAS8iIhJBCngREZEIUsCLiIhEkAJeREQkgjTZjIhEhialEUnQGbyIiEgEKeBFREQiSAEvIiISQQp4ERGRCFLAi4iIRJACXkREJIIU8CIiIhGkgBcREYkgBbyIiEgEKeBFREQiSAEvIiISQQp4ERGRCNJkMyKSl8ITy4jI0XQGLyIiEkEKeBERkQhSwIuIiESQAl5ERCSCFPAiIiIRpIAXERGJIAW8iIhIBCngRUREIkgBLyIiEkEKeBERkQhSwIuIiESQAl5ERCSCFPAiIiIRpIAXERGJIAW8iIhIBCngRUREIkgBLyIiEkEKeBERkQhSwIuIiESQAl5ERCSCFPAiIiIRpIAXERGJoBLfBYiINNqRI5y7di6nbltLs8oKNpadzJt9RrK9VTvflYnkDAW8iOSNouoqePhheOABfrttW9K2aoznB5/Pf55/M1vKOniqUCR3KOBFJC+UHdrHTyf/ENbOqXN7EY5PvVfO+Wvncs8V9wGXZbdAkRyjgBeRtOg9aUp8ed2D6Q3XskP7ePoPX2PwtrXxdVtbteOVAWeyp2lLRm9cyhnrFwPQ7uAefvWn78BzQ+BTn0prHSL5RAEvIrnt0CF+/cdvJYU7X/0q5x0Zy+EmTeOrznp/IT968SG67NtJaXUl/NM/wdSpHgoWyQ26i15EctuXv8zpm5Yl2o8/Dg8+mBTuANN7DeNTNz7EmpO6BiuOHIGrrqLT3h1ZLFYkdyjgRSR3/elP8MtfxpvfveBWuOWWY+6+pawDN/7T99jZvCxYsW0b//W3n4Jzma5UJOco4EUkN+3aBXfdFW++eMrH+NWYKxt828Y2J3PnlZPADIDz187lugUvZ6xMkVylgBeR3PS1r8HWrQBsadWOr024Ox7aDZnRcxh86Uvx9qQ3nqDNwb0ZKVMkVyngRST3LFgAjzwSb37rotvZ27Rl4z7je9+Dvn0BaHtoH//2zu/TWaFIzlPAi0ju+cY3EtfNJ0zg7wPPbvxnNG8OP/xhvHnDvJdgxYo0FSiS+xTwIpJbpk+HF14Ils3gwQdT7po/yqc+xYweQwBoUl0F3/9+mooUyX0KeBHJLd/8ZmL5+uth2LDj/ywz/uv8mxLtp56CNWuO//NE8ogCXkRyxuCta+C114JGURF861sn/Jlzu53KO71ivyRUVcEPfnDCnymSDxTwIpIzPjfruUTjmmtgwIC0fO7Pzr4u0fjNb+J354tEmQJeRHJC5z07uGLpm4kV992Xts+e0WMoc7ueEjSOHIH/+7+0fbZIrtJY9CKSE26cNyW4EQ7g3HNhzJgT+rzw5DeY8cTpVzBq03IAtvzn//Cx3UOpLC6pc2Kcmveme9IckWzSGbyIeFdSVcm1i15NrAgNUpMufzvlHOjcGYDO+z7kkhXT0/4dIrlEAS8i3l24ehYd9+8KGl26wBVXpP07jhQ3gdtvj7dvmP9S2r9DJJco4EXEu39a+EqiccstUJKhq4e33RbcnQ+c9cEiuu/akpnvEckBCngR8WvDBsatmZNo/8u/ZO67unaF8ePjzasXv5657xLxTAEvIn498QTFrhqAt3sNh379Mvt9N98cX7x68T+gujqz3yfiiQJeRPxxDn6fmATm2WEXZ/47r7iCXc1aAdBz91Z4663Mf6eIBwp4EfFn0SJYuhSAA02a8mr/MzP/nc2aMfnU8xPtJ5/M/HeKeKCAFxF//vCH+OJr/c/gYGmzrHztc6ddEGo8Fwx+IxIxCngR8cM5ePrpeDPprDrD5nU9hQ1lHYPGRx/BP/6Rte8WyRYFvIj48e67sG4dALubtuTNPqOy991mvHTKxxLtZ5/N3neLZIkCXkT8CJ29/33g2VSUNMnq108ZFAr4556Dioqsfr9IpingRST7nIO//CXefOHU87JewoIuA9lQdnLQ2LUrMU2tSERoshkRyYqkCVzmz4f164MNbdsyo+fQrNYAgBkvDvoYd8wMftH441f/h6+86Y65vyaekXyjM3gRyb7nn08sX3oplcV+zjX+fso58eUL1syiqGY2O5EIUMCLSPaFA37iRG9lLOgyADp1AqDDgd2M2LzCWy0i6aaAF5HsWrcu6KIHKC1NGhs+25wVweWXx9sXrXrXWy0i6aaAF5Hsmjw5sXzhhVBW5q8WgE9+Mr748VUzPRYikl4ZDXgzG29my81slZlNqmO7mdlPY9sXmtmo0La2ZvYnM1tmZkvN7KxM1ioiWZIj3fNxF10EzYIR9E7Z8QE9NIWsRETGAt7MioGfAxOAwcD1Zja41m4TgAGx123AL0PbfgL83Tk3CBgOLM1UrSKSHa0P74c330ysCJ09e9OiRRDyMeqml6jI5Bn8WGCVc26Nc64CeBqo/ev6ROBJF5gBtDWzLmZWBpwH/ArAOVfhnNuVwVpFJAvOfn8BVFYGjVGjgvnZc8EVV8QXFfASFZkM+G7A+lB7Q2xdKvv0BbYDvzazeWb2mJm1zGCtIpIF56+Zm2h4vLnuKKEb7casf48WFQc9FiOSHuaca3iv4/lgs08DlzjnPhdr3wiMdc7dHdpnCvAD59zbsfY/gPsBA2YA5zjn3jWznwB7nHPfqON7biPo3qdTp06nPx0a/vJE7du3j1atWqXt8/KdjkeCjkWyffv2sXZ34hnyod3aHLXPog27uOlLt1G2YzsA8378Y3YPHx5s27i7/vfWsT287kTUfN6Qm26hwwfrAHjh3q+zbuToOvdriP5uJNPxSJbu43HBBRfMcc6NrmtbJkeX2AD0CLW7A5tS3McBG5xzNX1lfwKOukkPwDn3KPAowOjRo924ceNOuPAa5eXlpPPz8p2OR4KORbLy8nIeent/vL3uM+OO2ue7n3+Yu2PhTuvWjPzCF6BJMP78zeER4+p4b13bw+tORM3nPdJ5FLfHAn5H+QIeKjmzzv0aor8byXQ8kmXzeGSyi34WMMDM+phZKXAdMLnWPpOBz8bupj8T2O2c2+yc2wKsN7NTYvt9HHgvg7WKSIYldc9fdFE83HNFeDa789bO81iJSHpk7AzeOVdpZncBLwPFwOPOuSVmdkds+8PAS8ClwCrgAHBL6CPuBp6K/XKwptY2Eckz560NBfwll/gr5Bhmdx/MwZKmNK88TL8PN9Bt9zY2tjnZd1kixy2jA0A7514iCPHwuodDyw648xjvnQ/UeV1BRPxImjCmMQ4e5Mz1ixLtWMD3bqCbvaHt6XS4pJR3ewxh3No5QPALyR9G5NCNgCKNpJHsRCTz3nqLZpWx+dYHDYLevb2Wcyzhbvpzwz0OInlIAS8imffKK4nlHOyer/Fmn5Hx5Y+9v4BizS4neUwBLyKZ9/rrieWLL/ZXRwNWte/BptYdACg7vJ/hml1O8pgCXkQy68MP47PHVVoRnHuu54LqYcbbvUfEm2e9v9BjMSInRgEvIpn1xhsQG1BrUecB/mePa8D0nsPiy2d9oICX/KWAF5HMCnXPT+811GMhqQkH/OiNSymtPOKxGpHjp4AXkcyaOjW+OK3ncI+FpGZLWQfWntQFgGaVFYzYvNxzRSLHRwEvIpmzdSssWQJARVEJs7uf6rmg1CR10+s6vOQpBbyIZE55eXxxfteBHGrSzF8tjTAjFPBJA/SI5BEFvIhkTqh7fnoedM/XmN4zca/AqI3LaFozSI9IHlHAi0jm5NkNdjW2t2rH6nbdAWhadYRRG5d5rkik8RTwIpIZGzfCypXBcrNmzOs6yG89jRQ+iz9Tj8tJHsroZDMiUsBC3fPvdBrI4ZLStHxstiagmd5zGDfM/xuQ/Dx8+PsbPemOSBbpDF5EMuOtt+KLM3rkT/d8jRmhM/gRm1bAgQMeqxFpPAW8iGRGKOBn9TjNYyHHZ2fLtqxs3wOA0upKmDXLc0UijaOAF5H027EDli4Fguff53cZ6Lmg4zOr++BE4+23/RUichwU8CKSftOmxRcXd+6XN8+/1zZbAS95TAEvIukXCsNZ3fOve77GzHDt06ZBleaHl/yhgBeR9Atdf086C84zG9p0YkurdkFjzx5YpFHtJH8o4EUkrZodOQRz5sTbs7vlx/jzdTJjdvgsXt30kkcU8CKSVsM3r4QjwRSrK9v34KMWbTxXdGJ0o53kKwW8iKTVmA1L4suz8rh7vkbSJYa33gLn/BUj0ggKeBFJqzEb3osv5/MNdjWWduzN3tLmQWPTJrrv3uq3IJEUKeBFJG2KqqsYtXFpvB2FM/jqomLmhu4jCP8CI5LLFPAikjanbl9H64qDQaNrVza06eS3oDQJ/6KigJd8oclmRCRtRofC74WyfmDmsZrUpDJ5Tfg6/GgFvOQJncGLSNokX3/P/+75GvO7DISS4Hxo4M4PaHtwj+eKRBrWYMCbWZGZjTSzy8zsQjOLRp+biKTdyI3L4stzIhTwh5o0g1Gj4u0Rm5Z7rEYkNccMeDPrZ2aPAquAB4HrgS8Ar5rZDDO7xczUAyAiAHTau4Nue7cDcKBJU5Z17O23oHQ788z44qjQLzIiuaq+gP4e8Dugn3PuEufcDc65a5xzw4BPAm2AG7NRpIjkvhGbVsSXF3YeQFVRscdqMuCss+KLI3UGL3ngmDfZOeeur2fbNuDHGalIRPLSqE2Js9p5XQd5rCRDQmfwIzYvp6haE89Ibquvi/4/QssXZ6ccEclXI5MC/hSPlWRIr17QuTMArSsO0n/nes8FidSvvi768aHl/8x0ISKSv6yykmFbVsXbkTyDN9N1eMkruklORE5Yy9WraVZZAcD6Np3Y3uokzxVliK7DSx6pb6Cbk83sy4CFluOccz/KaGUikjfavJd4/j2S3fM1QmfwCnjJdfUF/P8BretYFhFJUlYoAT96NJVWRImrZuDOD2DXLmjb1ndVInWq7y76b2ezEBHJX+GAnxvF6+81WrRg6cl9GLp1ddCeORM+8Qm/NYkcQ73X4M3sAjP7s5ktib3+ZGbjslSbiOSD7dtpvmkTAIeLm/Bep76eC8qsud1Cv8BMn+6vEJEGHPMM3swuA/4X+E7sZcAo4HEzu8s591J2ShSRnPbuu/HFxZ36caS4Sb27pzK5Sy6qqfvKroO4aW7szzBjhseKROpX3zX4rwBXOucWhNbNN7PZwM8ABbyIJJ3FRvr6e0zSJYgZM6C6Gor0QJLknvr+VnauFe4AOOcWAppwRkQCobPYSD7/XssHbTuzs3lZ0Ni1C1asqP8NIp7UF/D7j3ObiBSKqqrgRrOYpOvTUWWm6/CSF+rrou9nZpPrWG9AtO+iEZHUvPce7NsHwJZW7djcuoPngrJjXtdBXLwq9ovN9Olwyy1+CxKpQ30BP7Gebf+d7kJEJA8lXX8fFAznWgCS7jUI9WCI5JL6noN/w8xGAv2AJc65pdkrS0TyQtL19+jfYFdjUecBVGMU4WDxYjhwAFq08F2WSJL6ZpP7JvAMcDUwxcw+n7WqRCQ/hB6RK6SA39e0Bavbdw8aVVUwb57fgkTqUN9Ndv8EjIjNCz8GuC07JYlIXti7F5YGHXuuqIjFnfp7Lii7FnQZmGiom15yUH0Bf8g5dwDAObezgX1FpNDMnQvOAbC/d28OljbzXFB2LegyINFQwEsOSvUueqvVxjn3yYxWJiK5bdas+OLeUwqne75G0hl86FiI5IrG3EWvO+dFJCEc8IMK4Pn3WpZ17AOlpVBRAatXw86d0L6977JE4uq9iz6bhYhIngkF/J5TToG1HmvxoKKkCQwfnjgOs2bB+PF+ixIJqW+ymReAR4G/O+eO1NrWF7gZWOecezyjFYpI7tm+HdbGEr1pU/b37QtrK9Ly0Xk1Gc3YsUcF/KKNu7l50hTWPXjZUbuH/2x1bRdJp/punPs8cC6wzMxmmdlLZva6ma0BHgHmKNxFCtTs2YnlESNwTeqfQS6yxo5NLOtGO8kx9XXRbwHuB+43s95AF+AgsKLm7noRKVDhm8rGjPFXh2/hP/vMmfGnCkRyQX032cU559YB6zJaiYjkDwV84JRToHXrYEyAbdtg/XrfFYnE6dl2EWkc5xTwNYqKjj6LF8kRCngRaZSue7fD1q1Bo3Xr4Cy2kIUDXs/DSw5RwItIowzbvDLRGD06OIstZLrRTnJUgz+ZZna5mc0zsw/NbI+Z7TWzPdkoTkRyz/BwwBdy93yNcMDPno1VV/mrRSQklV+9fwzcBLR3zpU551o758oyXJeI5KhhW1YkGgp46NYNOncOlvfto+3mTX7rEYlJJeDXA4ud0/MfIoXOXDVDN69KrFDAg1nSWXyn1Svq2Vkke1J5TO5+4CUzewM4XLPSOfejjFUlIjmpz4ebKKuIDYPRsSP07Om3oFwxdixMDubi6rR6JYy+xHNBIqmdwX8fOAA0A1qHXiJSYIZvDp2djh0bnL1KUk9GpzWr6tlRJHtSOYNv55z7RMYrEZGcN2yLbrCr0+jR8cUOH6yjaWV6xuUXORGpnMG/ZmYKeBFJPoNXwCe0awf9+wNQXFXJqdsKbGo9yUmpBPydwN/N7KAekxMpYEeOMDgcXAr4ZKEb7YZt1o124l+DAR97LK7IOddcj8mJFLDFi2lW0/Xcq1dwk50khH7hGR6+lCHiSUqTzZjZScAAghvtAHDOvZmpokQkB2n8+fqFjknSaH8injQY8Gb2OeAeoDswHzgTmA5cmNnSRCSnhIdhDY/eJoGRI6G4GKqq6LdzA+zZA2Xq7BR/UrkGfw8wBnjfOXcBMBLYntGqRCT36Ay+fi1awGmnAVCEgzlzPBckhS6VgD/knDsEYGZNnXPLgAKfPkqkwBw4AEuWAFCNwemney4oR2lmOckhqQT8BjNrC/wVeNXMngc02LJIIZk3D6qCSVRWt+8eTBMrR1PASw5p8Bq8c+5TscVvmdlUoA3w94xWJSK5JRRWC7sMYIDHUnKaAl5ySL1n8GZWZGaLa9rOuTecc5OdcxqmSaSQhMJqQWfF+zENHUplkybB8vvvw3bdriT+1BvwzrlqYIGZaUYJkUIWDvguAz0WkuOaNGFHzz6Jts7ixaNUrsF3AZaY2T/MbHLNK9OFiUiO2LULVgbPdR8pKmbZyX0aeENh29q3f6KhgBePUhno5tsZr0JEclfoca9lHXtzuKTUYzG5b1vf0CUMBbx4lMpNdm9koxARSb/ek6YAsO7By47/Q2bPji8u0vX3Bh11Bu+cptUVLxrsoq+ZXKbWa72ZPWdmfbNRpIh4lHT9XQHfkI+6dGNvafOgsW0brF/vtyApWKlcg/8R8BWgG8FwtfcB/wc8DTyeudJEJCfoDL5xiopY3FnX4cW/VAJ+vHPuEefcXufcHufco8ClzrlngJMyXJ+I+LR9e/C4F0DTpqzooAdqUpH0KKECXjxJJeCrzeza2DPxRWZ2bWiby1RhIpIDwuOpjxhBZXFKE1AWvIXhRwkV8OJJKgH/GeBGYFvsdSNwg5k1B+7KYG0i4psmmDkuC8P3KsyeDdXV/oqRgpXKXfRrgCuOsfnt9JYjIjkldP2d0aNhqb9S8smGspOhY8fgEseePcE4Aqdoji7JrlTuou8eu2N+m5ltNbM/m1n3bBQnIp7VDnhJjZnGpRfvUumi/zUwGehKcCf9C7F1IhJlmzYFL4CWLWHQIL/15BsFvHiWSsB3dM792jlXGXs9AXTMcF0i4lv47H3UKCgu9ldLPlLAi2epBPwOM7vBzIpjrxuAnal8uJmNN7PlZrbKzCbVsd3M7Kex7QvNbFSt7cVmNs/MXkztjyMiaaPu+RMTDvh58+DIEX+1SEFKJeD/BbgW2AJsBq4BbmnoTWZWDPwcmAAMBq43s8G1dpsADIi9bgN+WWv7Pei2HhE/FPAn5uSToWds3IBDh2DJEr/1SMFJJeB7OOc+6Zzr6Jw72Tl3JdAjhfeNBVY559bE5o9/GphYa5+JwJMuMANoa2ZdILi5D7gMeCzlP42IpIdzCvh0UDe9eJRKwP8sxXW1dQPCgzBviK1LdZ8fA/cDeoBU5AT1njQl/krJBx8Ej3gBe5q2pM+vlh/zvb0nTWHRxt3pKjVa6gn4Rv3/EDkOx3wO3szOAs4GOprZl0ObyoBU7rapa/qk2iPf1bmPmV0ObHPOzTGzcfV+idltBN37dOrUifLy8hRKS82+ffvS+nn5TscjIV+Oxb1DK49al0rdHd54gyGx5d39+vLlYdWEf9cOf8a9Qyvp1Lzu74q6uo5lzbEoLy+nbZMmjIit3zt1Kvd+6+hZ/fLh79GJyJeflWzJ5vGob6CbUqBVbJ/WofV7CK7DN2QDyV353YFNKe5zDfBJM7sUaAaUmdnvnHM31P6S2Nj4jwKMHj3ajRs3LoXSUlNeXk46Py/f6Xgk5MuxuLmOM8R1nxnX8Btffjm+OKX1QB5alPxPRfgzbp40hXuHVh61TyGo61j+7KnneWhRSbBt5Ei4914AWq9dy//OreJwk6YNfkaU5MvPSrZk83gc8ycyNg/8G2b2hHPufQAzKwJaOef2pPDZs4ABZtYH2AhcB/xzrX0mA3eZ2dPAGcBu59xm4GuxF7Ez+PvqCncRyZDQ9feFmkHu+LVpE4xgt3w5VFVx2rY1zO12qu+qpECkcg3+B2ZWZmYtgfeA5Wb2lYbe5JyrJBir/mWCO+Gfdc4tMbM7zOyO2G4vAWuAVQRT0H7heP4QIpJGtW6wW6g54E9M6Dr8sM0rPRYihSaVPrXBzrk9ZvYZgkD+KjAH+GFDb3TOvRR7T3jdw6FlB9zZwGeUA+Up1Cki6bB6NezaFSy3bx+Mqy7Hb8wY+N3vABi+eYXnYqSQpHIG38TMmgBXAs87546gaWJFoqv243FW172wkrLwGfwWncFL9qQS8I8A64CWwJtm1ovgRjsRiSI9/55eI0ZASdBZ2u/DjZQd2ue5ICkUDQa8c+6nzrluzrlLYwPSvA9ckIXaRMQHzQGfXs2bw5Ah8eaQLas8FiOFpL7n4G9wzv2u1jPwYT/KUE0i4ktVFcydm2iPHg3T5/urJyrGjIH5wXEcvmUl03qPaOANIieuvjP4lrH/tj7GS0SiZsUK2BfrQu7cGbp29VtPVOhOevGgvufgH4n999vZK0dEvKrdPa8b7NJDAS8e1NdF/9P63uic+2L6yxERr3SDXWacdho0awaHDtFt73Y67P+IHS1P8l2VRFx9XfRzQq9P1mrPyXxpIpIt8YlPQgF/86JqTYbSgFQnjOn9jVeY0653vH28Z/GaoEYao74u+t/ULJvZv4XbIhI9xdVVMG9evL1IQ9Sm1cIuAzh90zIAhm9eyev9x3quSKIulefgQQPbiETewB3vw6FDQaNnT3a2bOu3oIhZ0GVgfHnYFo1oJ5mXasCLSMQN3Rx6PlvX39MuPGnPsM0rgzH/RTLomNdxBUwAACAASURBVAFvZnvNbI+Z7QGG1SzXrM9ijSKSBUnDqCrg025tu67sKW0BQPuDe+i+Z5vniiTqjhnwzrnWzrmy2KsktNzaOVeWzSJFJPOGbtEZfCY5K2JRl/7xth6Xk0xTF72IUFp5hFO3rU2sUMBnxMLOoevwmllOMkwBLyKcsn0dpdWVQaNfPzhJz2hnwoIuievwwzWznGSYAl5EdP09SxaGAn7IllWYq/ZYjUSdAl5EdP09Sza17sj22OOHrSsO0nfnRs8VSZQp4EWE4eHrwZoiNnPMWNA53E2v6/CSOQp4kUJ34AADdnwQLJvByJF+64m4heEBb3QnvWSQAl6k0M2dS0nNteBTT4UyPQWbSeEBb4Yr4CWDjjkWvYgUiJkzE8tjNT56uhxrUpjwnfSDt62BigooLc1WWVJAdAYvUugU8Fn1UYs2rG/TCYCmVUdg8WLPFUlUKeBFCp0CPuvCN9oxa5a/QiTSFPAihWz7dlgbjGB3uLgJDB3quaDCEH4eXgEvmaKAFylkoXBZ0qmvrgVniQJeskEBL1LIQt3z4fnKJbMWdepPNRY0Fi+G/fv9FiSRpIAXKWShgJ+vgM+a/U1bsLp996BRXQ3z5vktSCJJAS9SqJzTGbxH6qaXTFPAixSqtWth504AdjVrxbqTunouqLDoTnrJNAW8SKGqffZu5rGYwhMeslYBL5mggBcpVLr+7tXSk/twpKg4aKxaBR995LcgiRwFvEih0vV3rw6XlLKsY+/EitmzvdUi0aSx6EUK0ZEjMHduvJl0w5c0Wnjc+XsbMVbQwi4DGLp1ddCYNQsuvjjNlUkh0xm8SCFavBgOHgyWe/ViR8uT/NZToBZ01nV4yRwFvEgh0vjzOUGPykkmKeBFCpECPies7NATmjcPGhs3wubNfguSSFHAixQiBXxOqCoqhlGjEit0Fi9ppIAXKTR798KSJcFyUVFywEj2jRmTWFbASxop4EUKzdy5wTC1AKedBq1a+a2n0CngJUMU8CKFJtw9f8YZ/uqQQO2Ar/nlS+QEKeBFCo2uv+eW/v2hbdtg+cMPgzkCRNJAAS9SaKZPTywr4P0zg9GjE21100uaKOBFCkiXPduDx7EguPY+ZIjfgiSg6/CSAQp4kQIyauOyRGPsWCgu9leMJCjgJQMU8CIFZOSmUMCfdZa/QiRZ+FLJnDlQVeWvFokMTTYjEjHhiU9qG9VAwNf3Xsmgbt2gS5dgJLv9+2HpUl0+kROmM3iRAlFaeYTTamYuAz0il2vCZ/EzZvirQyJDAS9SIIZsXUXTqkoA1pzUFTp08FyRJDn77MTytGn+6pDIUMCLFIiRoRvs5nUb5LESqZMCXtJMAS9SIEZuWh5fnttVAZ9zTj8dSmK3RS1fDjt3+q1H8p4CXqRAhG+wm6eAzz3NmydP/KPr8HKCFPAiBaDznh103bsDgP1NmrG8Yy/PFUmd1E0vaaSAFykA4bP3BV0GBvOQS+5RwEsaKeBFCsCojUvjy3N1g13uCo9NMHMmVFb6q0XyngJepADoBrs80b079OgRLB84AAsX+q1H8poCXiTiSiuPMGTrqnh7ftdTPFYjDVI3vaSJAl4k4k7bujo+wM3ak7rwYYs2niuSeoW76cNT+4o0kgJeJOLCN9ipez4P6Axe0kSTzYhEXPIIdqd6rEQa0nvSFEqqKlnVvDkcPAjr1jHmrt+yvVU736VJHtIZvEiUOceYje/FmzqDz32VxSVJ88OPCv2CJtIYCniRCOuxeyud9n0IwJ7SFizTADf5IXQd/vTQI44ijaGAF4mwMRuWxJfndD+Vag1wkx9C1+EV8HK8FPAiETZ6Q6J7flb30zxWIo0SOoMfsnUVpZVHPBYj+UoBLxJhY0IBP7v7YI+VSKN07Aj9+wPQtKqSIVtWNfAGkaMp4EUi6qQDuxmwcz0AFUUlLOg8wHNF0ijnnBNfHBu61CKSKgW8SESNDl27XdS5P4ebNPVYjTTaeefFF8co4OU4KOBFIir5+ru65/POuefGF8dseI+i6iqPxUg+UsCLRFT4rG9WD91gl3f694fOnQEoO7yfU3a877kgyTcKeJEIanbkEEO2rI6352gEu/xjlnQWP3a9uumlcRTwIhE0fPNKSquDCWZWtO/JruZlniuS4xLuplfASyMp4EUiKNw9r8fj8ljoRrszNiwG5zwWI/lGAS+Sh3pPmkLvSVOOuX2MbrCLhiFD2NO0JQAd9++i90ebPBck+UQBLxIxRdVVjAo9IqeAz2PFxUn//3QdXhpDAS8SMaduX0frioMAbG7Vng1tOnmuSE5EeIhhDXgjjaGAF4mYMz9YFF+e3X1wcDe25K2ZoUccx65f7LESyTcKeJGICQf89F7DPFYi6bCoc38OlgSjEPbcvZXOe3Z4rkjyhQJeJEqqqpLO8qb3VMDnuyPFTZjX9ZR4W930kioFvEiULFxIm8P7Adjaqh1rT+rquSBJh6Tr8OqmlxQp4EWipLw8vjijx1Bdf4+Id5Ouw+sMXlKjgBeJkqlT44szeg71WIik07yug6goKgFg4M4PYOtWzxVJPlDAi0RFVRW8+Wa8qYCPjoOlzZKuw4d7akSORQEvEhULFsDu3YCuv0dR0hMRr7/urxDJGwp4kajQ9fdIS3oiQgEvKVDAi0RFKOCnq3s+cuZ1HcShktKgsWoVfPCB34Ik55X4LkBE0kDX3/NWXZMG1bWuoqQJs7oN5tz35wcrpk6Fm27KdHmSx3QGLxIFoevvW1q1Y52uv0eSrsNLYyjgRaKg9uNxuv4eSdN6DU80Xn9d88NLvRTwIlHw6qvxRQ1PG12LOvdnb2nzoLFhQ3AtXuQYMhrwZjbezJab2Sozm1THdjOzn8a2LzSzUbH1PcxsqpktNbMlZnZPJusUyWuHDyddf3+790iPxUgmVRUV826PIYkV6qaXemQs4M2sGPg5MAEYDFxvZoNr7TYBGBB73Qb8Mra+ErjXOXcqcCZwZx3vFRGAadPgYDD/O/37s7HNyX7rkYyaHu6mD12aEaktk2fwY4FVzrk1zrkK4GlgYq19JgJPusAMoK2ZdXHObXbOzQVwzu0FlgLdMlirSP4Kdc9z8cX+6pCsmFb7Rjtdh5djyORjct2A9aH2BuCMFPbpBmyuWWFmvYGRwLt1fYmZ3UZw9k+nTp0oT+MQjvv27Uvr5+U7HY8E38fi3qGVAJSXlzPqL3+hLLZ+cZcu3Nu3MqXPCNdf83mpqv3eTs0b/xlRlfFjUd2dI38uo8mePbB9O7Mef5z9/fpl7vtOkO+flVyTzeORyYCv6zbe2r9q1ruPmbUC/gz8m3NuT11f4px7FHgUYPTo0W7cuHHHVWxdysvLSefn5TsdjwTfx+Lm2HPS68YPhRUrgpVFRQy5+24uf/CdlD5j3WfGHfV5qar93nuHVvLQIg2rAWTlWNx9ySXwxz8CMOajjyCHfy59/6zkmmwej0x20W8AeoTa3YFNqe5jZk0Iwv0p59xfMlinSP6aOjXRRTtmDLRt67ceyY5LLkksv/yyvzokp2Uy4GcBA8ysj5mVAtcBk2vtMxn4bOxu+jOB3c65zWZmwK+Apc65H2WwRpH8puvvhSkc8G+9BQcO+KtFclbGAt45VwncBbxMcJPcs865JWZ2h5ndEdvtJWANsAr4P+ALsfXnADcCF5rZ/Njr0kzVKpK3FPCFqXt3GBx7sOjwYXjjDb/1SE7K6IUi59xLBCEeXvdwaNkBd9bxvrep+/q8iMT0/GgzrF0bNFq2hDPP9FuQZNcll8B77wXLL78MEyb4rUdyju6KEclTH6uZdATg/POhtLRR769rQhPJD70nTeG8DW15smaFrsNLHTRUrUieOn/NnERD3fMF590eQzhc3CRoLFum6WPlKAp4kTxUWnmEc95fkFih7tmCc7hJ0+Rha3UWL7Uo4EXy0OgNS2hVERuetm9fGDjQb0HixZt9QvMOvPKKv0IkJyngRfLQBWtmJxqXXqrpYQvUm31GJRqvvQaVGk1QEhTwInnogtW1Al4K0ooOvaBbbJqOXbtg5ky/BUlOUcCL5Js1a+j/4YZguVmznB6mVDLMLHnQmyl6MkISFPAi+eZvf0ssX3ghNG/urxbx7/LLE8svvOCvDsk5CniRfPNSaOwodc/LxRcnxkBYtAjef99vPZIzFPAi+eTgwWAO8BoKeGnVKujJqfHii/5qkZyigBfJJ+XlcOgQACvb94A+ffzWI7lB3fRSBwW8SD4J/eNd3vd0j4VITgkH/NSpsG+fv1okZyjgRfJFdTU8/3y8+Vr/MzwWIzmlVy8YNixYrqhInmVQCpYmmxHJF7Nnw6ZNAHzYvIzZ3Qcnbc7m5DGaqCYHXXEFLFwIwLMP/IL7302efGjdg5el9etq/g6k+3MlfXQGL5IvQmfv/+g3lqqiYo/FSM4JddNfsHo25qo9FiO5QAEvki/++tf44isDNfe71DJ2LJx8MgAdD+xixKYVngsS3xTwIvlgxQp4771guXlz3uo9wm89knuKiuCyRHf5JSumeSxGcoECXiQfhLrnueQSDjVp5q8WyV1XXx1fnLBiGjjnsRjxTQEvkg9C3fNceaW/OiS3XXQRlJUB0GvXFk7btsZzQeKTAl4k123ZAtOnB8tFRcnPPIuENW0a3E0fM365uukLmQJeJNc991yiq/Xcc6F9e7/1SG4LddNfuvwdddMXMAW8SK575pnE8jXX+KtD8sP48eyP3aPR78MNDNjxgeeCxBcFvEgu27wZ3nwzWC4qUsBLw5o3Z2rf0fHmBN1NX7AU8CK57E9/SnSxnn8+dO7stx7JC3875Zz48oTl73isRHxSwIvksnD3/LXX+qtD8srUfqM5VBIMVXvq9nX03bnBc0XigwJeJFdt2ADvxM6+iouTbp4Sqc+B0uZJsw1OfO8Nj9WIL5psRiTDjjUxS4OTdPzxj/HFN3sM47yOHdNZlkTc86eez/gVweOVE98rDy71mCXtowljok1n8CK5KtQ9P2XQxzwWIvno9X5j2FPaAoDeuzbDzJmeK5JsU8CL5KJVq+DddwE4UlTMywPP8lyQ5JvDTZry8sCzEyueespfMeKFAl4kF/32t/HF8r6ns6t5mcdiJF/99bRxicYzz0BlpbdaJPsU8CK5proannwy3vzzkI97LEby2fSeQ9naql3Q2LYN/vEPvwVJVingRXLN22/DunXBctu2vN5vrNdyJH9VFxXzwqBzEyvUTV9QFPAiuSZ09s5111FR0sRfLZL3/nraBYnGX/4Ce/f6K0aySgEvkksOHIBnn020b7rJXy0SCYs79WNF+55BY//+5L9fEmkKeJFc8vzziTOsAQPgjDP81iP5z4xnhl2caD/2mL9aJKsU8CK55PHHE8uf/exRA5OIHI/nhlwITWKXembMgCVL/BYkWaGAF8kVq1bBa68Fy0VFQcCLpMGHLdrAxImJFb/6lb9iJGsU8CK54tFHE8uXXgo9e/qrRaLn1lsTy7/9LRw+7K8WyQoFvEguOHwYfv3rRPv22/3VItF08cXQo0ewvGMHTJ7stx7JOE02I5IGdU3acaxJZur03HPBP7oAPXrQ981qqt9u+P2N+g6JtAb/LhQXwy23wHe+E7QffhjGfDm93yE5RWfwIrngkUcSy5/7HNVFxf5qkei69dbg/g6A119nwPb3/dYjGaWAF/Ft6VIoLw+Wi4uTr5WKpFPPnnDllfHmTXNf9FiMZJoCXsS3n/wksXzFFdCtm79aJPruvju+eNWS1yk7tM9jMZJJCngRn3buTB6a9p57/NUiheH882HIEABaHDnMpxe95rkgyRQFvIhPjz4KBw8GyyNGBP/4imSSWdJZ/I1zpwQzGErkKOBFfKmogP/930T7S1/SyHWSHZ/5DJx0EgC9d22GF17wXJBkggJexJc//hE2bQqWO3eG667zW48UjpYt4fOfT7QffBCc81ePZIQCXsQH5+C//zvRvvNOKC31V48Unnvu4XBxbCiUGTPgrbf81iNpp4AX8eDC1bNg/vyg0by5Rq6T7OvalT8P+Xii/eCD/mqRjFDAi2Sbc9w97ZlE+447oGNHf/VIwXp07FVUWSwG/vY3WLDAb0GSVgp4kSw75/0FjNy8PGiUlsJ99/ktSArWunbd+Nsp5yRW/OAH/oqRtFPAi2TZF6c9nWjceit07eqvGCl4vzzj6kTj2Wdh8WJ/xUhaabIZkSw66/2FnLE+9g9oSQl89atJ2+ubzEMTfUhdGvv3ovb+Szr3h8svhxdfBOd4+arPc/tVX0/aJzyJUiqfner+klk6gxfJFuf46htPJNqf/Sz06uWtHJG47343vnjJyhkM27zCYzGSLgp4kSwZv2IaI2L/cB4ubgIPPOC5IpGYESPg2mvjzfve/K3HYiRdFPAiWVBcXcVXQv9oPjnqsmBmL5Fc8e1vx6eSPW/dPM56f6HnguREKeBFsuDTC1+l34cbANhT2oKfn3VtA+8QybJBg4LLRjFff/0xiqqrPBYkJ0oBL5JhZYf2cd9bibP3R864ml3NyzxWJHIM3/kOB0uaAnDatjVcu/BVzwXJiVDAi2TYv739ezoc2A3AptYdeHz0RM8ViRxDjx788sxr4s373votrQ/v91iQnAgFvEgGDdy+js/OfTHe/t6Fn+NgaTOPFYnU79Gxn2Jj62BkxQ4HdvPFd/7guSI5Xgp4kUxxjm+/9gglLphre1rPYbwUHjVMJAcdatKMH1xwS7x9y+zJGsI2TyngRTLk2oWvctYHiwCotCIeuOh2zfcueeHFQecyo8cQgOAX1M99DiorPVcljaWAF8mEjRv5+uuPxZuPj57Iyo4a1EbyhBlfG393MF4DwOzZ8JOf+K1JGk0BL5JuzsHtt1NWcQCAtSd14UfnfsZzUSKNs7ZdN35yzvWJFd/4Br0+2uSvIGk0BbxIuv3mNzAlMd73/RPu4VAT3Vgn+efRsVfx3sl9gsbBg/z4hYcoqVJXfb7QZDMSSeEJNU5k4otGT56xfDncdVe8+cSoy5kVu5bZ0HeI+HKsv4OVxSV8ZcI9TPn9V+DIEUZuXs4X3/kDPzrvxpQ/e9HG3dw8aYomoPFAZ/AiadK0sgKuuw72B88Nr27Xnf86/ybPVYmcmCWd+8P3vhdv3zX9Wc6I3TwquU0BL5Im/z71VzB/ftBo2pS7P3k/B0qb+y1KJB3uuw8uvBCAIhw/feGHdNz3oeeipCEKeJE0+PTCV7hpbqib86GHeK9TX38FiaRTURE8+SQfxoZY7rTvQx5+7j8orTziuTCpjwJe5ES9/Tbff/kXifbVV8MXvuCvHpFM6NaNuz95P1UWxMbpm5bx7Vd/GTw1IjlJAS9yItauhauuorQ6dmfxsGHwxBMa0EYi6Z3eI/iPcYlR7q5f+Ar89397rEjqo4AXOV6bN8NFF8H27QDsaNEGJk+GVq08FyaSOb8acyXPDR6XWHH//cGjoZJzFPAix6Fk71645BJYswaAw8VNuONT/w69NFqdRJwZkyZ8kXfDj3/eeiu8+OKx3yNeKOBFGuvDDxl2//2wKPaoUHExd06cxOzup/mtSyRLDpeU8vmrvs7Sjr2DFVVVcPXVXLxyhte6JJkCXqQxtm6FceMoW7Ysse7Xv+a1AWf4q0nEgz3NWvHZa78DfWIj3VVU8Iu//oBLl73ttzCJU8CLpGr1ajjvvMSZuxk8/DDcmPqoXiJRsr1VOygvh/79AWhSXcXPJv8X/PznfgsTQAEvkpo33oCxY2HFCgBc7Llgbr/dc2EinvXsGfx8DBoEQLGrDoZrvueeoOtevFHAi9THOW6c+2Jwt/yHsZG7mjZlyQMPwA03+K1NJFd07Qrl5czvMiCx7qc/hQkTaL57l7+6Cpwmm5GckK7JYdL6vTt38uhz3+cT4RuHOneGv/6VHQcPNvh5jflukUw7kb9vKb23Uyeuu/4HPDTlf7hs+TvBuldf5fqZs3l9/L30npTYVRPPZIfO4EXqMnkyDB+eHO4jR8LMmXCGbqgTqcuhJs24a+JX4RvfiA/21HL3R/zumW/wwGuP0PLwAc8VFhYFvEjYhg1w1VUwcSJs3JhYf/fdMG0a9OjhrzaRPOCsCL7zHfj73+Hkk4Fggppb5rzAK7+6k4tWvqvhbbNEAS8ClB3aB1/7GgwcCM89F1+/o0Ubbr36G8H1xGbNPFYokmc+8QmYP5/3h42Kr+q2dzuP/eW7wcx0M2d6LK4w6Bq8FLSTDuzmhnkvcevs5+HQvuSNt97KRS0/zq7YDFoi0khdujD5K99g5bPv8MBrj9L+4J5gfXl5cKlr4sRgqNuzz/ZaZlTpDF4K06JF8K//yrRf/gv3vv0UbcPhPnJk8NjPY48p3EVOlBmTB4/j459/mCdHXsaRouLEtuefh3POCV7PPAOHDvmrM4IU8FI4Nm+GH/84CPBhw+Dhh2leeTixvV8/+MMfYPbsYEAbEUmbXc3L+OYn/pWLb/0FXHNN8sZp0+C666BLF7jzTnjnHT1DnwYKeIms4uoqhm5eCd/+NowZEzyr+6Uvwfz5Sfst7tSPey6/F5YuDf6RKdKPhUimrGvXDf74R1i8GG65BZo0SWzctQt+8Qv42MeCn9dbbw3uidm501/BeUzX4CUanIMtW2DhQpgxg98+/VdGbl5Oq4q6n1enWTO48kqus2HM6DEUzPhJ+B8aEcms006Dxx+H730PHnsMnngC1q5NbN+2Ldj++ONBe8iQoGftzDNhxIhg5Dz9zNYrowFvZuOBnwDFwGPOuQdrbbfY9kuBA8DNzrm5qbxXCtSuXbBuXeK1enVwJrBoUdJv+efW9d6SEjj3XPjnf4ZPfxratGGGBpsR8atrV/jmN+HrX4c334SnnoIXXggmdgpbvDh4/eIXQbu0NPglYfDg4PJazatv3+DxvOLio7+rwGQs4M2sGPg5cDGwAZhlZpOdc++FdpsADIi9zgB+CZyR4nslnzgHlZVw5AgcPgz798O+ffHXhatm0rLiIC2OHIIfLYft249+bd0Ke/c26ms3te5A1ysnwOWXB4/ttG2boT+giJyQoiIYNy54VVfDu+8GA05NnRrcF1P7mnxFBcybF7xqKy6Gjh2DkSdrXu3bQ5s2R7/KyqBFi6BXL/xq2jTvL9dl8gx+LLDKObcGwMyeBiYC4ZCeCDzpnHPADDNra2ZdgN4pvDdzHnoIpk5l6I4dwV8KSB6YoWa5rnUNbT+e92TiM4/jPafv2wctWx69vbo6CO6aV0VFcvvIkSDc6/F4uPH3enc9tlatgm68kSP54oaWzO4+mE1lJ2tYTJF8U1QEZ50VvCA4EZgxA956K7iHZv58+OCDY7+/qiq4ZLdly4nVUVoahH1pafBLw7FeJSV1rzdLvIqKwIzhu3YFTwzEZuDLpEwGfDdgfai9geAsvaF9uqX43syZPx+mTKF91r4wP7T2XQBA8+bQu3fy69RTYejQYFar2G/ck9X1LhIdrVoFEz5ddFFi3UcfBffcrFwZXKqrea1bl5gY6kRVVASvNDoJgl9YsiCTAW91rKs9PuGx9knlvcEHmN0G3BZr7jOz5SlX2LAOwI40fl6+8388Dh4M7nZfujTlt9h/nvh+dWzzfyxyyBd1POIK7VjU9XNTa91RxyPVn8mI6sDIken8+9HrWBsyGfAbgPDA3d2BTSnuU5rCewFwzj0KPHqixdbFzGY750Zn4rPzkY5Hgo5FMh2PBB2LZDoeybJ5PDJ5B8EsYICZ9TGzUuA6YHKtfSYDn7XAmcBu59zmFN8rIiIix5CxM3jnXKWZ3QW8TPCo2+POuSVmdkds+8PASwSPyK0ieEzulvrem6laRUREoiajz8E7514iCPHwuodDyw64M9X3epCRrv88puORoGORTMcjQccimY5HsqwdD3Oal1dERCRy8vspfhEREamTAr4WM/uumS00s/lm9oqZdQ1t+5qZrTKz5WZ2ic86s8XMfmhmy2LH5DkzaxvaVojH49NmtsTMqs1sdK1thXg8xsf+vKvMbJLverLNzB43s21mtji0rp2ZvWpmK2P/PclnjdliZj3MbKqZLY39jNwTW1+ox6OZmc00swWx4/Ht2PqsHQ8F/NF+6Jwb5pwbAbwIfBPAzAYT3M1/GjAe+EVsSN2oexUY4pwbBqwAvgYFfTwWA1cBb4ZXFuLxCA0pPQEYDFwfOw6F5AmC/99hk4B/OOcGAP+ItQtBJXCvc+5U4Ezgztjfh0I9HoeBC51zw4ERwPjY02JZOx4K+Fqcc3tCzZYkBtiZCDztnDvsnFtLcOf/2GzXl23OuVecczXjzM4gGJMACvd4LHXO1TWYUiEej/hw1M65CqBmSOmC4Zx7E6g9bNpE4Dex5d8AV2a1KE+cc5trJgtzzu0FlhKMSlqox8M552qGrGsSezmyeDwU8HUws++b2XrgM8TO4Dn2sLqF5F+Av8WWdTySFeLxKMQ/cyo6xcbzIPbfkz3Xk3Vm1hsYCbxLAR8PMys2s/nANuBV51xWj0dBBryZvWZmi+t4TQRwzv0/51wP4Cngrpq31fFRkXgEoaHjEdvn/xF0wT1Vs6qOjyqY41HX2+pYF4njUY9C/DNLA8ysFfBn4N9q9YgWHOdcVexyb3dgrJkNyeb3Z/Q5+FzlnLuo4b0A+D0wBXiA1IbezUsNHQ8zuwm4HPi4SzxXWbDH4xgiezzqUYh/5lRsNbMuzrnNsdkxt/kuKFvMrAlBuD/lnPtLbHXBHo8azrldZlZOcL9G1o5HQZ7B18fMBoSanwSWxZYnA9eZWVMz60Mwh/3MbNeXbWY2Hvgq8Enn3IHQpoI8HvUoxOOhIaXrNhm4KbZ8E/C8x1qyxswM+BWw1Dn3o9CmQj0eHWueOjKz5sBFBHmSteOhgW5qMbM/A6cA1cD7wB3OuY2xbf+P4Dp0JUH309+O+UERYWargKbAztiqGc65O2LbCvF4fAr4GdAR5Qx8ywAAAVVJREFU2AXMd85dEttWiMfjUuDHJIaU/r7nkrLKzP4AjCOYMW0rQW/fX4FngZ7AB8CnnXNpmr80d5nZx4C3gEUE/34C/DvBdfhCPB7DCG6iKyY4mX7WOfcdM2tPlo6HAl5ERCSC1EUvIiISQQp4ERGRCFLAi4iIRJACXkREJIIU8CIiIhGkgBcREYkgBbyIiEgEKeBF5LiY2RgzWxib97plbM7rrI61LSLHpoFuROS4mdn3gGZAc2CDc+4HnksSkRgFvIgct9gY9LOAQ8DZzrkqzyWJSIy66EXkRLQDWgGtCc7kRSRH6AxeRI6bmU0Gngb6AF2cc3d5LklEYgpyPngROXFm9lmg0jn3ezMrBqaZ2YXOudd91yYiOoMXERGJJF2DFxERiSAFvIiISAQp4EVERCJIAS8iIhJBCngREZEIUsCLiIhEkAJeREQkghTwIiIiEfT/AfqkzNkpaQ/pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# let's see how the MLE performs\n",
    "plot_normal_mle()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - MLE for m-Dimensional Gaussian\n",
    "---\n",
    "Given $\\{x_i\\}_{i=1}^n$ i.i.d samples of $X \\sim N(\\mu, \\Sigma)$, what is the MLE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "The final results are pretty much the same, but with vectors and matrices, though the math is a little more complicated.\n",
    "$$ \\hat{\\overline{\\mu}}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n \\overline{x_i} $$\n",
    "$$ \\hat{\\Sigma}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n (\\overline{x_i} - \\hat{\\overline{\\mu}}_{MLE}) (\\overline{x_i} - \\hat{\\overline{\\mu}}_{MLE})^{T}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/scissors.png\" style=\"height:30px;display:inline\"> Vector & Matrix Deriviatives\n",
    "---\n",
    "* $\\nabla_x Ax = A^{T}$\n",
    "* $\\nabla_x x^{T} A x = (A + A^{T}) x$ \n",
    "* $\\frac{\\partial}{\\partial A} \\ln |A| = A^{-T}$\n",
    "* $\\frac{\\partial}{\\partial A} Tr[AB] = B^{T}$\n",
    "\n",
    "Using the above, we will use the following:\n",
    "1. $\\nabla_{\\mu} {\\mu}^{T} \\Sigma^{-1} x_i = \\Sigma^{-1} x_i$ \n",
    "2. $\\nabla_{\\mu} {\\mu}^{T} \\Sigma^{-1} \\mu = (\\Sigma^{-1} + {\\Sigma}^{-T}) \\mu$\n",
    "3. $\\frac{\\partial}{\\partial \\Sigma^{-1}} \\ln |\\Sigma^{-1}| = \\Sigma^{T} = \\Sigma$\n",
    "4. $\\frac{\\partial}{\\partial \\Sigma^{-1}} Tr[\\Sigma^{-1} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu}) (\\overline{x_i} - \\overline{\\mu})^{T}] = \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu}) (\\overline{x_i} - \\overline{\\mu})^{T}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Solve for the d-dimensional case\n",
    "* $p(x|\\mu, \\Sigma) = \\frac{1}{(2\\pi)^{\\frac{nd}{2}} |\\Sigma|^{\\frac{n}{2}}} e^{- \\frac{1}{2}\\sum_{i=1}^n (x_i - \\mu)^{T} \\Sigma^{-1} (x_i - \\mu)}$\n",
    "* $\\ln p(x|\\mu, \\Sigma) \\propto -\\frac{n}{2} \\ln |\\Sigma^{-1}| -\\frac{1}{2} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu})^{T} \\Sigma^{-1} (\\overline{x_i} - \\overline{\\mu}) $\n",
    "* $\\nabla_{\\mu} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu})^{T} \\Sigma^{-1} (\\overline{x_i} - \\overline{\\mu}) = \\sum_{i=1}^{n} (-2\\Sigma^{-1} \\overline{x_i} + (\\Sigma^{-1} + {\\Sigma}^{-T}) \\mu) = 0 $ $$ \\rightarrow \\hat{\\overline{\\mu}}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n \\overline{x_i} $$\n",
    "\n",
    "* **The Trace Trick** - $\\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu})^{T} \\Sigma^{-1} (\\overline{x_i} - \\overline{\\mu}) = \\sum_{i=1}^n \\textit{Trace}\\big((\\overline{x_i} - \\overline{\\mu})^{T} \\Sigma^{-1} (\\overline{x_i} - \\overline{\\mu})\\big) = \\textit{Trace}\\big(\\Sigma^{-1} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu}) (\\overline{x_i} - \\overline{\\mu})^{T} \\big)$\n",
    "* $\\frac{\\partial}{\\partial \\Sigma^{-1}}\\big( \\frac{n}{2} \\ln |\\Sigma^{-1}| -\\frac{1}{2} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu})^{T} \\Sigma^{-1} (\\overline{x_i} - \\overline{\\mu}) \\big) = \\frac{\\partial}{\\partial \\Sigma^{-1}}\\big( \\frac{n}{2} \\ln |\\Sigma^{-1}| -\\frac{1}{2} Tr[\\Sigma^{-1} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu}) (\\overline{x_i} - \\overline{\\mu})^{T}] \\big) = $ $$ \\frac{n}{2} \\Sigma - \\frac{1}{2} \\sum_{i=1}^n (\\overline{x_i} - \\overline{\\mu}) (\\overline{x_i} - \\overline{\\mu})^{T} = 0 $$ $$ \\rightarrow \\hat{\\Sigma}_{MLE} = \\frac{1}{n} \\sum_{i=1}^n (\\overline{x_i} - \\hat{\\overline{\\mu}}_{MLE}) (\\overline{x_i} - \\hat{\\overline{\\mu}}_{MLE})^{T} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/dusk/64/000000/infinity.png\" style=\"height:50px;display:inline\"> Asymptotic Properties of MLEs\n",
    "---\n",
    "Let $X_1, X_2, ..., X_n$ be a random sample from a distribution with a parameter $\\theta$. Let $\\hat{\\theta}_{ML}$ denote the MLE of $\\theta$. Then (under some conditions):\n",
    "\n",
    "1.  $\\hat{\\theta}_{ML}$ is **asymptotaically consistent**: $$ \\lim_{n \\to \\infty} P(|\\hat{\\theta}_{ML} - \\theta| \\geq \\epsilon) = 0, \\forall \\epsilon >0 $$\n",
    "2. $\\hat{\\theta}_{ML}$ is **asymptotaically unbiased**: $$ \\lim_{n \\to \\infty} \\mathbb{E}[\\hat{\\theta}_{ML}] = \\theta $$\n",
    "3. As $n$ becomes large, $\\hat{\\theta}_{ML}$ is approximately a **normal random variable**, that is, the random variable: $$ \\frac{\\hat{\\theta}_{ML} - \\theta}{\\sqrt{Var(\\hat{\\theta}_{ML})}} $$ **convereges in distribution** to $\\mathcal{N}(0,1) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/task.png\" style=\"height:30px;display:inline\"> Exercise - Asymptotic Consequences of MLE When Choosing The Wrong Model\n",
    "---\n",
    "Definitions:\n",
    "* **Kullback-Leibler (KL) Divergence** - defined to be $$ KL(p(x)||q(x)) = \\underset{x \\sim p(x)}{\\mathbb{E}}[\\log \\frac{p(x)}{q(x)}],$$ it is a way to measure \"distance between distribution\" (how far is $p(x)$ from $q(x)$).\n",
    "* **Entropy** - defined to be $$ H(p(x)) = -\\underset{x \\sim p(x)}{\\mathbb{E}}[\\log(p(x))] $$\n",
    "* **The Weak Law of Large Numbers** - states that if you have a sample of independent and identically distributed random variables, as the sample size grows larger, the *sample mean* will tend toward the population mean $$ \\lim_{n \\to \\infty} P(|\\overline{x}_n - \\mu | \\geq \\epsilon) = 0 $$\n",
    "    * In words: as the sample size $n$ grows to infinity, the probability that the sample mean $\\overline{x}_n$ differs from the population mean $\\mu$ by some small amount $ \\epsilon$ is equal to 0.\n",
    "\n",
    "\n",
    "Let the data $\\mathcal{D}$ be drawn i.i.d. from a distribution $p(x)$ which is not necessarily contained in the parametric model (i.e. **there is no $\\theta$ for which $f(x;\\theta) = p(x), \\forall x$**). What is the *upper bound* for the MLE as $N \\to \\infty$? Use the definitions above to guide you and use the *log* of the estimator ($\\log (\\prod_{i=0}^{N-1}f(x_i;\\theta))$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### <img src=\"https://img.icons8.com/dusk/64/000000/idea.png\" style=\"height:30px;display:inline\"> Solution\n",
    "---\n",
    "We will use the above definitions to derive another form of the estimator. For every $\\theta$ it holds that:\n",
    "$$ \\log(\\prod_{i=0}^{N-1}f(x_i;\\theta)) = N\\big(\\frac{1}{N}\\sum_{i=0}^{N-1}\\log(f(x_i;\\theta))\\big) $$\n",
    "Using the **weak law of large numbers** as $N \\to \\infty$: $$ \\underset{N \\to \\infty}{\\rightarrow} N \\cdot \\mathbb{E}_{x \\sim p(x)}[\\log(f(x;\\theta))] =  $$ $$ = N \\cdot \\mathbb{E}[\\log(\\frac{f(x;\\theta)p(x)}{p(x)})] $$ $$ =N \\cdot \\mathbb{E}[\\log(\\frac{f(x;\\theta)}{p(x)})] +  N \\cdot \\mathbb{E}[\\log(p(x))] =  $$ $$ -N \\cdot KL(p(x)||f(x; \\theta)) -N \\cdot H(p(x)) $$\n",
    "Thus, the MLE takes the form:  $$ \\lim_{N \\to \\infty} \\theta_{MLE} = \\lim_{N \\to \\infty} argmax_{\\theta} \\log(L(\\theta)) $$ $$ = argmax_{\\theta} \\big( -KL(p(x)||f(x; \\theta)) - H(p(x)) \\big) = argmin_{\\theta} KL(p(x)||f(x; \\theta))  $$ \n",
    "\n",
    "* What is the meaning of this?\n",
    "    * In the limit of **large sample size**, the MLE criterion is equivalent to find the model which is *\"closest\"* (in the KL divergence sense) from the true distribution. That is, when the model cannot be really parameterized with $\\theta$ we are not guaranteed to converge to the true distribution, even for infinite number of samples!\n",
    "    * KL divergence satisfies: $KL(p(x)||q(x)) \\geq 0$, where equality holds $\\iff p(x) = q(x)$. That means that for evey $\\theta$ the likelihood $l(\\theta)$ satisfies $ l(\\theta) \\leq -H(p(x)) $. In the case of model mis-specification and in the limit of large sample size, the log-likelihood estimation cannot reach its maximal value (or the negative log-likelihood cannot reach it minimal value), the **entropy** is the lower bound (to the negative log-likelihood) or, the upper-bound of the log-likelihood is the negative entropy.\n",
    "        * Note that we assumed that **we can't** model $p(x)$, and thus it is independent on $\\theta$. If we could model it with $\\theta$, then $p(x)$ would have been dependent on $\\theta$ as well, and the entropy would have been dependent on $\\theta$ (and we couldn't solve this, as we don't know $p(x; \\theta)$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def plot_misspecified_mle():\n",
    "    mu_real_1 = -5\n",
    "    mu_real_2 = 5\n",
    "    var_real = 4.0\n",
    "    num_samples = 1000\n",
    "    samples = [np.random.normal(mu_real_1, np.sqrt(var_real), size=(num_samples // 2)),\n",
    "               np.random.normal(mu_real_2, np.sqrt(var_real), size=(num_samples // 2))]\n",
    "    samples = np.concatenate(samples, axis=-1)\n",
    "\n",
    "    mu_mle = np.sum(samples) / num_samples\n",
    "    var_mle = np.sum(np.square(samples - mu_mle)) / num_samples\n",
    "\n",
    "    print(\"mu mle: \", mu_mle)\n",
    "    print(\"var mle: \", var_mle)\n",
    "\n",
    "    x = np.linspace(-15, 15, 10000)\n",
    "    f_x_mle = (1 / np.sqrt(2 * np.pi * var_mle)) * np.exp(-0.5 * (np.square(x - mu_mle)) / var_mle)\n",
    "\n",
    "    # set bins for histogram\n",
    "    n_bins = 100\n",
    "    bins_edges = np.linspace(samples.min(), samples.max() + 1e-9, n_bins + 1)\n",
    "\n",
    "    fig = plt.figure(figsize=(8, 8))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.grid()\n",
    "    ax.set_ylabel('Histogram (PDF)')\n",
    "    ax.set_xlabel('x')\n",
    "\n",
    "    # plot histogram\n",
    "    ax.hist(samples, bins=bins_edges, density=True)\n",
    "    # plot estimation\n",
    "    ax.plot(x, f_x_mle, linewidth=3, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu mle:  0.04869170509203468\n",
      "var mle:  30.252569170292713\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAHgCAYAAABJrX+JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZxN9R/H8ddnxtgaa0lFISkp/SqSkhpaLC3aNyqV/LTQXuSnfd83pdKmTXspIimkEEqhUhIhohTGOsv398e5c++daWbcGffeM/fM+/l43If7/Z7lfpx7Zj5zzvku5pxDREREgiXN7wBEREQk/pTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCqIrfAcTTDjvs4Jo2bRq3/a1fv57tttsubvtLdToeEToWhel4ROhYFKbjUVi8j8esWbP+dM41KG5ZoBJ806ZNmTlzZtz2N3HiRLKysuK2v1Sn4xGhY1GYjkeEjkVhOh6Fxft4mNnikpbpFr2IiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEUEITvJl1NbP5ZrbAzAYWs7ylmU01s81mdk0xy9PN7Bsz+zCRcYqIiARNwhK8maUDQ4FuQCvgLDNrVWS11cAA4P4SdnM58EOiYhQREQmqRF7BtwMWOOcWOue2ACOBHtErOOdWOudmADlFNzazxsCxwPAExigiIhJIiRyLvhGwJKq8FDi4DNs/DFwH1CptJTPrC/QFaNiwIRMnTixblKXIzs6O6/5SnY5HhI5FYToeEToWhel4FJbM45HIBG/F1LmYNjQ7DljpnJtlZlmlreucexp4GqBt27YunoP4a5KEwnQ8InQsCtPxiNCxKEzHo7BkHo9E3qJfCuwaVW4M/B7jth2AE8xsEd6t/c5m9nJ8wxMREQmuRCb4GUALM2tmZlWBM4FRsWzonBvknGvsnGsa2u5T51yvxIUqIiISLAm7Re+cyzWzy4BxQDrwnHNunpn1Cy0fZmY7ATOB2kC+mV0BtHLOrU1UXCIiIpVBIp/B45wbA4wpUjcs6v0KvFv3pe1jIjAxAeGJiIgElkayExERCSAleBERkQBK6C16EZGgajpwdPj9oruP9TESkeLpCl5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAIooQnezLqa2XwzW2BmA4tZ3tLMpprZZjO7Jqp+VzP7zMx+MLN5ZnZ5IuMUEREJmiqJ2rGZpQNDgaOBpcAMMxvlnPs+arXVwADgxCKb5wJXO+e+NrNawCwzG19kWxERESlBIq/g2wELnHMLnXNbgJFAj+gVnHMrnXMzgJwi9cudc1+H3q8DfgAaJTBWERGRQElkgm8ELIkqL6UcSdrMmgIHANPjEpWIiEglYM65xOzY7DSgi3OuT6h8DtDOOde/mHVvBrKdc/cXqc8EJgF3OOfeKeFz+gJ9ARo2bNhm5MiRcfs/ZGdnk5mZGbf9pTodjwgdi8Iq4/GYs2xN+H3rRnXC7yvjsSiNjkdh8T4enTp1muWca1vcsoQ9g8e7Yt81qtwY+D3Wjc0sA3gbeKWk5A7gnHsaeBqgbdu2Lisrq1zBFmfixInEc3+pTscjQseisMp4PHoPHB1+v6hnVvh9ZTwWpdHxKCyZxyORt+hnAC3MrJmZVQXOBEbFsqGZGfAs8INz7sEExigiIhJICbuCd87lmtllwDggHXjOOTfPzPqFlg8zs52AmUBtIN/MrgBaAfsB5wBzzGx2aJc3OOfGJCpeERGRIEnkLXpCCXlMkbphUe9X4N26L2oKYImMTUREJMg0kp2IiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEUBW/AxCpaJoOHB1+v+juY8u9bdHtt2W/IiJlpSt4ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEAUoIXEREJICV4ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEAUoIXEREJIE02IyJSirJOEjRn2Rp6l7CNJhySZNIVvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIACU0wZtZVzObb2YLzGxgMctbmtlUM9tsZteUZVsREREpWcISvJmlA0OBbkAr4Cwza1VktdXAAOD+cmwrIiIiJUjkFXw7YIFzbqFzbgswEugRvYJzbqVzbgaQU9ZtRUREpGSJTPCNgCVR5aWhukRvKyIiUumZcy4xOzY7DejinOsTKp8DtHPO9S9m3ZuBbOfc/eXYti/QF6Bhw4ZtRo4cGbf/Q3Z2NpmZmXHbX6qrLMdjzrI14fetG9Updp2SjkX0tkW3j2W/qSrI50ZJ31tJ9StXr+GPjZRpmyAL8rlRHvE+Hp06dZrlnGtb3LJEDlW7FNg1qtwY+D3e2zrnngaeBmjbtq3Lysoqc6AlmThxIvHcX6qrLMej0DCjPbOKXaekYxG9bdHtY9lvqgryuVHS91ZS/WOvvM8Dc6qUaZsgC/K5UR7JPB6JvEU/A2hhZs3MrCpwJjAqCduKiIhUegm7gnfO5ZrZZcA4IB14zjk3z8z6hZYPM7OdgJlAbSDfzK4AWjnn1ha3baJiFRFJlKZF7uqIJEtCZ5Nzzo0BxhSpGxb1fgXe7feYthUREZHYaCQ7ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEAUoIXEREJICV4ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEASuhY9CJ+i57oY9Hdx8Ztv3OWrQlP/RnP/YqIxIuu4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCSAleREQkgJTgRUREAkgJXkREJIA02YxIjKInrrm6tY+BSFxEf5+gSYMkeHQFLyIiEkBK8CIiIgGkBC8iIhJASvAiIiIBpAQvIiISQErwIiIiAaQELyIiEkBK8CIiIgGkBC8iIhJASvAiIiIBpAQvIiISQErwIiIiAaTJZkQqkOgJUDT5SeWh710SQVfwIiIiAbTVK3gzSwP+A+wCbATmOef+SHRgIiIiUn4lJngzaw5cDxwF/AysAqoDe5rZBuAp4EXnXH4yAhUREZHYlXYFfzvwJPBf55yLXmBmOwJnA+cALyYuPBERESmPEhO8c+6sUpatBB5OSEQiIiKyzUpsZGdmd0a9Pzo54YiIiEg8lNaKvmvU+3sSHYiIiIjEj7rJiYiIBFBpjex2NLOrAIt6H+acezChkYmIiEi5lZbgnwFqFfNeREREKrjSWtHfksxAREREJH5KfQZvZp3M7G0zmxd6vWVmWUmKTURERMqptJHsjgUeB24NvQw4EHjOzC5zzo1JTogiqSt6EhEJrujv+erW8duXJp6RbVHaFfy1wInOueedc98652Y7554DTsQbwnarzKyrmc03swVmNrCY5WZmj4aWf2dmB0YtuzJ012Cumb1mZtXL+p8TERGprEpL8Ds5574tWumc+w5ouLUdm1k6MBToBrQCzjKzVkVW6wa0CL364g2Ni5k1AgYAbZ1z+wLpwJlb/d+IiIgIUHqCX1/OZQXaAQuccwudc1uAkUCPIuv0AEY4zzSgrpntHFpWBahhZlWAmsDvMXymiIiIUHo3ueZmNqqYegN2j2HfjYAlUeWlwMExrNPIOTfTzO4HfsObovZj59zHMXymiIiIAFZkorjIArMjStvQOTep1B2bnQZ0cc71CZXPAdo55/pHrTMauMs5NyVUngBcBywE3gbOAP4B3gTecs69XMzn9MW7vU/Dhg3bjBw5srSwyiQ7O5vMzMy47S/VpeLxmLNsTfh960Z1tmmb6PqGNeCPjaWvU5qSYilPvBVBqp8bUPbvpKTvOvrcKE0s+0qlc6AkqXhuJFK8j0enTp1mOefaFrestH7wk8zsAKA5MM8590MZP3cpsGtUuTH/vs1e0jpHAb8651YBmNk7wKHAvxK8c+5p4GmAtm3buqysrDKGWbKJEycSz/2lulQ8Hr2jWyT3zNqmbXoXaimdywNzqpS6TmlKiqU88VYEqX5uQNm/k5K+6+hzozSx7CuVzoGSpOK5kUjJPB6lzSZ3I/A6cAow2swuKuO+ZwAtzKyZmVXFayRX9Jb/KODcUGv69sAa59xyvFvz7c2sppkZcCRQ1j8wREREKq3S/sw8A9jfObfBzLYHxuINWRsT51yumV0GjMNrBf+cc26emfULLR8GjAG6AwuADcD5oWXTzewt4GsgF/iG0FW6iIiIbF1pCX6Tc24DgHPuLzMr88xzocFwxhSpGxb13gGXlrDtTcBNZf1MERERib0VvRUp45w7IaGRiYiISLmVluCL9lm/P5GBiIiISPyU2oo+mYGIiIhI/JQ22cwHeA3bxjrncoos2x3oDSwKjU8vknRBmZRDE9JUTrF870E5x8Ufpd2ivwi4CnjYzFYDq4DqQFPgF+Bx59z7CY9QREREyqy0W/Qr8EaVu87MmgI74w0b+1NB63oRERGpmLY+3BLgnFsELEpoJCIiIhI3Ze7bLiIiIhWfEryIiEgAKcGLiIgE0FYTvJkdZ2bfmNlqM1trZuvMbG0yghMREZHyiaWR3cPAycAcV9Lk8SIiIlKhxHKLfgkwV8ldREQkdcRyBX8dMMbMJgGbCyqdcw8mLCoRERHZJrEk+DuAbLxR7KomNhwRERGJh1gSfH3n3DEJj0REUo9zkJOD5eRATg5kZPgdkYiExJLgPzGzY5xzHyc8GpEkKTrRhybyiMEvv8CkSTB7Nnz7LSxZAitWwMaNHFGwTmYm7LQTNGkC++8PBx4IWVmwyy4+Bh4/fk4MpHNWyiqWBH8p3nj0m4EcwADnnKud0MhExH+zZ8OIEfDBB7BgwdbXz8721luwACZMiNTvtx+ccAKcey60aJG4eEUkbKsJ3jlXKxmBiEgFsXkzvPQSPP64d6W+Nenp5JuRlp8P+fnFr/Pdd97r9tuhY0cYMABOOgnS0+Mbu4iExTTZjJnVA1rgNbQDwDk3OVFBiYgPNm+GJ5+E++6D33//9/KaNaFTJ2jf3rv9vueesPPOkJnJ5EmTyDriCPjnH1i+HH780bv6/+IL+Pxz7/l8gc8/91577QWDBkGvXkr0Igmw1QRvZn2Ay4HGwGygPTAV6JzY0EQkKZyDd9+Fa6+FhQsLL6teHU45Bc45x3uWXq1ayfsxg3r1vFerVnDyyV59djZ8/DG8+CKMHg15eV79/PnQuzc8/DA89JC3fxGJm1gGurkcOAhY7JzrBBwArEpoVCKSHIsWwdFHe0k8OrnvtBPce693Nf7yy9ClS+nJvTSZmV6yf/99WLoUbrgB6tSJLJ8927szcNZZ8Oef2/TfEZGIWBL8JufcJgAzq+ac+xHYK7FhiUhCOQfDhkHr1oUbw9WvD488Ar/+6l3R160b38/daSe44w5YvBhuuglq1IgsGzkS9tkH3nknvp8pUknFkuCXmlld4D1gvJm9DxTzgE5EUsLq1XDccXDxxd7tc4C0NK/h288/e/9Wr176PrZVnTpw883w009w9tmR+pUrvbsJF10EGzcmNgaRgNtqgnfOneSc+8c5dzMwBHgWODHRgYlIAsyc6fVNHzMmUteyJUyd6l2516+f3HgaN4ZXXoEPPyzcV374cDj0UK/vvYiUS6kJ3szSzGxuQdk5N8k5N8o5tyXxoYlIXL3wAnTo4N0eL3DVVfDNN9CunW9hAXDssTB3rvccvsDs2dCmDYwf719cIims1ATvnMsHvjWz3ZIUj4jEm3Nw441w/vmwJfS3eZ068N578MADib8dH6t69byr+SefhKqhaS/WrIHu3eHZZ/2NTSQFxdIPfmdgnpl9BawvqHTOnZCwqEQkPrZsgT59vIFrCrRu7TVk22MP/+IqiRn06wcHHQQnnui1us/N9f4Pv/ziNdAz8ztKkZQQS4K/JeFRiEj8bdzodU8bOzZS16ULvPEG1K7gI023aQPTp3uNAb/5xqu76y74+28YOtRrFCgipYplqNpJyQhEpDjlmWAjnhOCVJTJRco8scj69dCjR+EucH36wBNPpM6Mb7vsApMne8/lP/zQqxs2DDZs8G7ZV4lpIM5y8fN7F4mXrf4ZbGbrzGxtkdcSM3vXzHZPRpAiUgbr1nnPraOT+5Ah8PTTqZPcC2RmeqPs9ewZqRsxwutaFz38rYj8Syx/Aj+I1+/9VbyZ5M4EdgLmA88BWYkKTkTKaMMGr0X6559H6u64wxs9LlVVqeINc1ujhtd9DuDNN72GeCNG6Ha9SAli+cno6px7yjm3zjm31jn3NNDdOfc6UC/B8YlIrHJy4LTTCif3++5L7eReID3duwMxYECk7pVX4NJLvV4CIvIvsST4fDM7PdQnPs3MTo9app8skYogPx/OO6/wADb33w/XXONfTPFm5k1M89//RuqGDYOBA5XkRYoRS4LvCZwDrAy9zgF6mVkN4LIExiYisXAO+veH116L1A0eDFdf7V9MiWLmtaKPHt723nvhwQf9i0mkgoqlFf1C4PgSFk+JbzgiUmb33ee1ji9w8cVw223+xZNo6eneqHzZ2TBqlFd37bXQrFlkiloRiakVfeNQi/mVZvaHmb1tZo2TEZyIbMW773q3qAuceSY8/njwB4PJyIDXX/eG3gXvLkavXvDVV/7GJVKBxHKL/nlgFLAL0Aj4IFQnIn76+msvqRU8f+7Y0buyrSytyqtX94bbLRiRb+NGOP54b6pbEYkpwTdwzj3vnMsNvV4AGiQ4LhEpzbJlXjLbsMErN2/uDT9brZq/cSXbDjt4DQu3394rr1zpHZeCaXBFKrFYEvyfZtbLzNJDr17AX4kOTERKsGmTN07777975bp1YfRoL9lVRi1aeFfyBX/czJsHF1yglvVS6cWS4C8ATgdWAMuBU4HzExmUiJRiwABvXnfwBoF5+23Yay9/Y/LbYYd5XeYKvPmm1/hQpBKLJcHv6pw7wTnXwDm3o3PuRGDXRAcmIsV49ll45plI+aGHoHNn/+KpSHr3hsuieu4OGgQff+xbOCJ+i2Wo2seAA2OoE6mUyjMxSbkmM5k1yxu5rUDPnoXL4vWH//ZbbzS//HyvV8GsWV4XukqipHOrzBMWScorMcGb2SHAoUADM7sqalFtID3RgYlIlL/+glNOgc2bvXLr1vDUU8HvDldWGRne7fk2bbyGiH//7SX5zz/3xq4XqURKu0VfFcjE+yOgVtRrLd5zeBFJBufg/PNh8WKvXKeO12J+u+38jauiatgQ3norMp3sV195I/uJVDIlXsGH5oGfZGYvOOcWA5hZGpDpnFubrABFKrtzv/4QPvkgUjFiRKTvtxSvfXu4++7IWPz33w9ZWd5MeyKVRCyN7O4ys9pmth3wPTDfzK5NcFwiAuy9ciGDP3suUnHFFXDCCf4FlEquvBK6d4+UzzvPu20vUknEkuBbha7YTwTGALvhTTgjIglUPWcTj71/L9XycryK/ff3rkolNmlp3jzyu+zilf/6y5ukJi/P37hEkiSWBJ9hZhl4Cf5951wOmiZWJOFunDCcPVYv9Qo1a8LIkZVvpLpttcMO3ix7BcP3Tp6s/vFSacSS4J8CFgHbAZPNrAleQzsRSZCu87/g7G/HRioee0yD2ZTX4YfDjTdGyjfeCLNn+xePSJJsNcE75x51zjVyznV3nsVApyTEJlIp7bD+b+4cNzRc/qBlR68VvZTf4MFw8MHe+5wc6NWLarlb/I1JJMFKTPChMecxs6uKvoABsezczLqa2XwzW2BmA4tZbmb2aGj5d2Z2YNSyumb2lpn9aGY/hPrliwSbc9w5bij1N3o3yX6vtQODu1yq/u7bqkoVeOkl71EHwLx5XD35JX9jEkmw0q7gCzrZ1irhVSozSweGAt2AVsBZZtaqyGrdgBahV1/gyahljwBjnXMtgf8AP2ztM0VS3UnzPuOYn6eFy9d2v4K11TN9jChAWrTwusuF9JnxHu1/+87HgEQSq7R+8E+F/r2lnPtuByxwzi0EMLORQA+8rnYFegAjnHMOmBa6at8ZWA8cDvQOxbAF0P00CbSd1v7JLZ88FS6/dEB3vmi6v48RBVC/fjBqFIwdSxqO+0c/RLcLHmddNQ0aJMFT2lC1j5a2oXNua7fpGwFLospLgYNjWKcRkAusAp43s/8As4DLnXPrt/KZIqnJOe796BFqb/ZO8cV1d+KuLD13jzszeO452HdfWL2axmtXMeiz57mh62Vb31YkxZgrYc5kMzsvqngLcFP0cufci6Xu2Ow0oItzrk+ofA7QzjnXP2qd0cBdzrkpofIE4DrAgGlAB+fcdDN7BFjrnBtSzOf0xbu9T8OGDduMHDmy9P9xGWRnZ5OZqdujBfw4HnOWrSlUbt2oTrHLSqpPlIY14I+N8dvfPp+Oo/Nz3hMqZ8Y7g2/n95b7hJdH//8qoor8s1LcedLgs8/Y59Zbw/Xv3HAby1q1jsvnxfvcKElJ50Qs538yz6eKfG74Id7Ho1OnTrOcc22LW1Zigi+0ktk3zrkDyvKhoUZxNzvnuoTKgwCcc3dFrfMUMNE591qoPB/IwutnP8051zRU3xEY6JwrdZzJtm3bupkF82THwcSJE8nKyorb/lKdH8ej6MxY0TNiRS8rqT5Rrm6dywNzYpmMcet2WbuSj5+9lMwtXlZ45qATuaNzn0LrVPSZwCryz0qx54lzjG3Zga4/TQW8OyZdz3+cjVWrb/PnxfPcKE1J50Qs538yz6eKfG74Id7Hw8xKTPCx9IOH8g1sMwNoYWbNzKwqcCYwqsg6o4BzQ63p2wNrnHPLnXMrgCVmVtDx90gKP7sXCQbnuH3cE+HkvqB+Y+7vqIEiE86MIUdfzJrQs/cm/6zg6s/Vql6CJdYEX2bOuVzgMmAcXgv4N5xz88ysn5n1C602BlgILACeAS6J2kV/4BUz+w7YH7gzUbGK+OX4HybTeaF31ykf47pul7M5Q6PVJcOqzPrcduRF4fIFM0dx4DJ11pHgKK2R3ToiV+41zaxg9DoDnHOu9tZ27pwbg5fEo+uGRb13wKUlbDsbKPa2g0gQ1N24lpsmPB0uv3Rgd75uvLePEVU+b+17JMf/MJkjfv2aNBz3fPQox/Z+lC1VMvwOTWSblXgF75yr5ZyrHXpViXpfK5bkLiKl+9+nz7LDBq9B1PLM7bnv8PO2soXEnRk3dLmM7Ko1AGjx1xL6fxm/hroifkrYLXoRKVmHRbM5de6EcHnIMZeQXa2mjxFVXsvq7Mg9R0T+uOo3/S1arFrsY0Qi8aEEL5Jk1XM2cdfYx8LlD1t25JMWRYeIkGR6+YDuzGjkDbSZkZ/HHR8PxVy+z1GJbBsleJEku3LKq+y25g8A1lTbjluO7OtzROIsjRu6XEpOWjoA7ZZ+z6lzPvE5KpFtowQvkkR7r1zIhTPeC5dv73whqzLr+RiRFPi5QROeaXdSuHzDZ89Tf0PiB00SSRQleJEkMZfPbR8/SZXQrd+pu7XmzdZH+xyVRHv00DNZUqchAPU2rWPQZ8/7HJFI+SnBiyTJqXMm0DbUz3pLWhX+d/Qlmga2gtmUUZ0hR/cLl0+b+wkH/zbHx4hEyk8JXiQJ6mxcx8CJkavBZw86kV922NXHiKQkE5sfxOi9OoTLd4wbSkZejo8RiZRP4gdMFhGum/wi22/0xopaVqsBjx56Zpm2L2188Yo+Tn2ylDQ3QXnceuRFHP7r19TaspE9Vi+l7/R3GHroGdsaokhS6QpeJMH2W/4TZ80eFy7fetRFcZnURBLnj1o78EDUnAD9p75O41DPB5FUoQQvkkBp+Xnc/vETpIVGff5s9zaMa3GIz1FJLEYceCxzGjYHoHruFv736XCfIxIpGyV4kQQ6e/ZY9luxAIDN6RncdFQ/NaxLEflp6dx49MXhctefptLx1699jEikbJTgRRJk+/X/cO3kEeHyE+1P47d6O/sYkZTVN41a8ua+R4XLN3/ytBrcScpQghdJkIETX6DO5vUALKq7M8Pan+pzRFIe92Sdx9qq3jwBzVcv5fyZo3yOSCQ2SvAiCdB6+c+cNjcy1OnNR/2XzVWq+hiRlNef29Xj4cN6hssDvhzJjuv+8jEikdgowYvEm3PcOOGZcHH8Hu2Y2LytjwHJthpx4LHM32E3ADK3bGTQRI1wJxWfErxInB374xQOWvY9ADlp6dzZ6UKfI5JtlZtehZuPioxwd9L3E2m3ZK6PEYlsnRK8SBxVy91S6OruxQOP49f6jXyMSOJlapP9+LBlx3D5lvHDSM/P8zEikdIpwYvE0YUz3qPx2pUArK5Rm0c7nOVzRBJPd3S6gA0Z1QDYe9Uizp79kc8RiZRMCV4kThpkr+aSaW+Gyw8e1pO11TN9jEjibXntBjx+SGTI2msmv0Td0BDEIhWNErxInFwz+SUyt2wE4Kftd+O1/bv6HJEkwvCDTmJRXW88gzqb13PFlFd9iaPpwNHhl0hxlOOntzcAACAASURBVOBF4mCfP37htDmRbnG3d76QvLR0HyOSRNlSJYM7OkcaTvb6ZgwtVi32MSKR4inBi2wr5xgy4ZnwePOf7t6Wybu38TkoSaTxexzMF032A6CKy2fIp8PBOZ+jEilMCV5kG3X5aSrtQ12mci2t0NWdBJQZt3W+iDzzfoUevugbshbO9DkokcKU4EW2QdXcHG6Y+Fy4/NKBx/LL9rv6GJEky487NmPkf44Jl4d8+ixV8nJ9jEikMCV4kW3Qe9YomvyzAoB/qmfycIezfY5IkumBjucUGqf+nG/U4E0qDiV4kXLaYf3f9P9yZLj8cIezWVOjlo8RSbKtrlmHRzucGS5fMeVVqq9TtzmpGJTgRcrpqs9foVaoW9wv9Rvz8gHdfY5I/PBim+P5tV6k29zB74zcyhYiyaEEL1IO2/+2iDO++zhcvq3zheSmV/ExIvFLTnoGd3TqEy7vO2Gsus1JhaAEL1JWztHxledId/kATG56ABN312xxldkne7RjSpP/AJCWr25zUjEowYuU0VELvmLXed8BkGdp3Na5D5j5HJX4yozbjizcba6Tus2Jz5TgRcogIy+HwZ8ND5df2b8bPzdo4mNEUlHMb9CU1/7TJVz+36fDycjL8TEiqeyU4EXK4NyvR9Ps7+UArK22HQ8dpm5xEvFgx15srlnQbW4Z536tbnPiH7UKEolRvQ1ruPyL18LlRw49k79r1vExon+Lnnhk0d3H+hhJcpQ00UpZ6+Nldc06fHXiGXR89XkABnzxGm/v25l/atRO6OeKFEdX8CIxunLKq9TevB6AfxruzIg2x/kckVRE3x3TnYX1dgH8nW1ORAleJAYtVi2m5+yPwuUpZ59PTnqGjxFJRZVfJYO7Ol0QLvf6Zgx7/PmbjxFJZaUEL7I1zjHk0+HhbnFTmvyHXw88yOegpCIrOtvc4M+e9TkiqYyU4EW2ImvhTA5f9A3gdYu7Xd3iZGvMuL1zH/LxzpNOC2dxxMJZPgcllY0SvEgpquTlMuTTyNXX6/sdw487NvMxIkkVP+y4e6HZ5v736XDS8/N8jEgqGyV4kVL0+mYMzVcvBWBd1Ro80LGXzxFJKnmwYy/WVa0BQIu/lnB2VDsOkURTghcpQZ2N67jii0gL6McOPZO/tqvrY0SSav7crh5DDzkjXL5yyqvU3pTtY0RSmSjBi5Tgii9epW7ol/HiujvxQpsTfI5IUtHzbU9gSZ2GANTfuJYBUWMpiCSSErxIMZr/uYRzokYhuzPrArZUUbc4KbvNVapyZ9b54fJ5X39Is9XLfIxIKgsleJFiDP7sWaqEusVN23Vfxu15iM8RSSr7aK8OTG+8DwAZ+Xnc8NlzPkcklYESvEgRhy+cRefQTGD5eLOEqVucbJPQbHMF3eaOXjCdQxfN9jkoCToleJEo6fl5/C+qW9wb+x3NvIbNfYxIgmLuTnvwzr6dw+Uhnw4nTd3mJIE02YxIlLNmj2XPv7xhRbOr1uCBjuf4HJEUlegJYxLp3sPPpfv8KdTM2czeqxZx+nfjGbl/V7/DkoDSFbxISO1N2Vw15ZVw+Yn2p7Eqs56PEUnQrKy1PU8efGq4fPXnL5O5eYOPEUmQKcGLhPT/ciT1N64FYGntHXn2oBN9jkiC6Jl2J7GsVgMAGmz4h0unvuFzRBJUSvAiQLPVy+g964Nw+c5OF7C5SlUfI5Kg2pRRnXuyzguXL5j5Hrv+s8LHiCSolOBFgBs+e46MUIOnrxq3YsxeHXyOSIJs1N5H8PUuewFQLS+XgROf9zkiCSIleKn0Dl00m6MXTA+Xb+usbnGSYGbeeRZy7PwvaLdkro8BSRAlNMGbWVczm29mC8xsYDHLzcweDS3/zswOLLI83cy+MbMPExmnVF5p+XkM+XR4uPzWvkcyZ+cWPkYklcU3jVryXqsjwuUhE57BQoMricRDwhK8maUDQ4FuQCvgLDNrVWS1bkCL0Ksv8GSR5ZcDPyQqRpEzvhvP3qsWAbAhoxr3Hn6uvwFJpXLPEb3ZWKUaAK3/+IVT5n7qc0QSJIm8gm8HLHDOLXTObQFGAj2KrNMDGOE804C6ZrYzgJk1Bo4FhiOSALU2r+fqz18Kl588+FRW1trex4ikslleuwFPtzspXL5u0ovU3LLRx4gkSMw5l5gdm50KdHXO9QmVzwEOds5dFrXOh8DdzrkpofIE4Hrn3Ewzewu4C6gFXOOcO66Ez+mLd/VPw4YN24wcOTJu/4fs7GwyMzPjtr9U58fxmLNsTaFy60Z1il1WUn1pOrz6AgeOeQ+AddvvwMv3DiW3WrWYtm1YA/6oIL+HYzkmiZbMcyPW79cvZT03MjZtpNe1l5L592oAvupxGtNP61mmzyzr+R/UcyMVxPt4dOrUaZZzrm1xyxI5kl1xrZSK/jVR7Dpmdhyw0jk3y8yySvsQ59zTwNMAbdu2dVlZpa5eJhMnTiSe+0t1fhyP3kVGLVvUM6vYZSXVl6TJ37/z37GRph2DDz2fUT9tF3NcV7fO5YE5FWMgyFiOSaIl89yI5fv1U9nPjVr8fMh5PDDmIQD2G/0+V+zSjd9r7xjzHsp6/gf13EgFyTweibxFvxTYNarcGPg9xnU6ACeY2SK8W/udzezlxIUqlc3gz56jan4uADMb7c2ovQ/3OSKpzN7ZtxPf7bQHANVzt3D9xBd9jkiCIJEJfgbQwsyamVlV4ExgVJF1RgHnhlrTtwfWOOeWO+cGOecaO+eahrb71DnXK4GxSiVy6KLZHPPztHD5Vs0WJz5zlsZtnfuEyz1+mMSBy9S+WLZNwu4xOudyzewyYByQDjznnJtnZv1Cy4cBY4DuwAJgA3B+ouKR1FHaZCIlLYt1ApL0/DxunPBMuPzWvkfy3c57li3ACiaVJ1+RiBm77suHex3GcfOnADBkwnBOPuc+nGm4EimfhD5EdM6NwUvi0XXDot474NKt7GMiMDEB4UkldNa342j552IA1mdUV7c4qVDuzurN0QumUS0vlwOWz+eE7yfx/j6d/A5LUpT+NJRKo/ambK76PNKU44n2p6lbnFQoS+vuVGiSo+snvUj1nE0+RiSpTAleKo0BX7xWaLa44ZotTiqgJ9qfzqqadQHYZd2f9P3qXZ8jklSlBC+Vwu5/LeW8ryPd4u7sdAGbM2Lr8y6STNnVanL/4eeEy/2mv0XDdX/6GJGkKiV4qRQGf/ZseLa46Y330WxxUqG92foovt+xGQA1czZz3eQRPkckqUgJXgLv8IWzOPKXGQDkY+oWJxVeflp6odnmTpn7Ka2X/+xjRJKKlOAl0NKLzBb3ZuujmBcaUESkIpvaZD8+btE+XL5xwjOQoKHFJZiU4CXQen4zhhZ/LQEgu2oN7le3OEkhd2adz5Y0rzfzQcu+p/v8L3yOSFKJErwEVt2Na7lyyqvh8tBDTmdVZj0fIxIpm0X1G/Fim8g8W4MmPk+13C0+RiSpRAleAuuayS9Rb9M6AH6r05Dn2hadrVik4nvs0DNZXaM2ALuu+YMLZr7vc0SSKpTgJZD2WbGAs2ePDZdvO/IiNlep6mNEIuWztnomDx4WmT720qlv0CD7bx8jklShBC+BYy6fW8cPIy00O/HEZm0Yv8fBPkclUn6v7d+Vn7bfDYDMLRu5ftIL/gYkKaFiTGgtgRY9Gcqiu4/d6jrb6uS5n9Hm9x8B2JJWhVuO6qtucSlMk+lAXlo6tx3Zh5feuBGAU+dO4NX/dOXrxntv876LHt+SfkYl9egKXgKl1ub1DJz0fLg8vN2J/Fq/kY8RicTH580OZOyeh4TLt41/krTQ4E0ixVGCl0C5fMqrNFj/DwDLM7fn8UPO8Dkikfi5rfNFbKziDbG8z8qF9Jz9kc8RSUWmBC+B0WLVYnrP+iBcvrPTBWyoWsPHiETia1mdHXn8kNPD5Wsnv8T2oT9oRYpSgpdgcI5bPnmKKi4fgGm77ssHex/uc1Ai8fdMu5P5td7OANTevF4N7qRESvASCN3nf8Ghv30HQK6lcdPR/dSwTgJpS5UMbjnyv+Hy6XM+galTfYxIKioleEl5NbZs4n9R482POPA45jdo6l9AIgk2sXlbxkWNU8+ll0KeGtxJYUrwkvL6Tx3JLqH5slfVrMvDh53tc0QiiXfbkRexqWDwpm++gaee8jcgqXCU4CWltVi1mIu+ejdcvveI81hbPdPHiESSY2mdhgxtf1qkYvBgWLXKv4CkwlGCl5RlLp87Ph5KRqgv8FeNW/FW6yN9jkokeZ4++BQW1fUa3PHPPzBwoL8BSYWiBC8p69Q5E2i39HsActLS+d8xl+BMp7RUHpurVPVGaizw3HO0Cf1MiOi3oaSkehvWMGhi1Ih1B53ET2pYJ5XQZ80Pgh6RmRLvGvs4GXk5PkYkFYUSvKSkgRNfoP7GtQAsqdOQRzqc6XNEIj569FHYbjsA9vzrt0LtUqTy0mQzklTxmDjkoCVzOWPO+HD5xqP7sSmj+jbvVyRVNX1iDhcefBZDQt1FB3w5ktEtD2NxvV18jkz8pCt4SSkZeTncMe6JcPmjPQ/1blGKVHIvtDmeOQ2bA1A9dwu3j3sCnPM5KvGTEryklD4z3mPPv34DILtqDW45su9WthCpHPLS0hnUtT95oYamHRfPpsf3E/0NSnylBC8pY9d/VjDgi5Hh8oOH9WJF7R18jEikYpm70x68eOBx4fKQT4dTZ+M6HyMSPynBS2pwjrvGPkaN3M0AzNtxd15sc9xWNhKpfB7o2Ivfa3l/+O5QpLeJVC5K8JISTpsznsMWfwtAnqUxqOtl5KWl+xyVSMWzvlpNb7KlkDO/+5h2S+b6GJH4RQleKrwG2asZ8umz4fLwg07ku5339DEikYptfIv2hSajuWvs41TL3eJjROIHJXip8G4dP4zam9cDsKjuzjykyWREtuqmo/qRXbUGAM1XL+XyL171OSJJNiV4qdC6zv+Cbj99GS4P6tpffd5FYrCi9g7cnXV+uNx3+ju0Xv6zjxFJsinBS4VVZ+M6bhv/ZLj82n7HMLXJfj5GJJJaXtm/K1N3aw1AFZfPfWMe1jC2lYgSvFRY//v0WRqs/weAPzLrc1enC3yOSCS1OEvj+q4D2JBRDYCWfy7m0qlv+ByVJIsSvFRIHX/9mtPmfhIuDzn6Ys3zLlIOv9XbmfsOPzdcvnTqG/Dttz5GJMmiBC8VTu1N2dzz0aPh8od7HcbHex7iY0Qiqe3FA49jZqO9AcjIz4Pzz4cc3aoPOk02IxXOjROeYZd1fwKwukZtbj76vz5HlJqiJ/ZZdPexPkYifstPS+e6bpfz0fP9qZaXA998w71H9eGJQ06PeR86n1KPruClQjn652mcOndCuDz4mEv4c7t6PkYkEgwLt2/Mgx17hsuXf/Eqe65a5F9AknBK8FJh1NuwhjvHPh4uv7/3EXzU8jAfIxIJluEHncTsnVsAUC0vl4c/fECt6gNMCV4qBue4/eMnaLAh0mr+xqjhNkVk2+WlpXNN9yvZnJ4BQKuVv3LFFA2AE1RK8FIhHP/DZI6d/0W4fH3X/qypUcvHiESCacEOu3F3Vu9wud/0t2m7dJ5/AUnCKMGL7xpkr/7XgDYTmx/kY0QiwfZCm+OZ0uQ/AKS7fB788EG227zB56gk3pTgxVfm8nlg9EPU3ZQNwNLaO3JH5z4+RyUSbM7SuKb7laypth0Au635gyGfDvc5Kok3JXjx1QUz3ufwRd8AkI9xzbFXkF2tps9RiQTfito7MOSYi8PlM7/7mKN/nuZjRBJvSvDim33++IXrJ70YLg9rfwrTdtNY8yLJMmrvI/igZcdw+a6xj8GKFT5GJPGkBC++qLFlE4+Ouo+q+bkAfLtTCx46rOdWthKRuDLjf8dcworM+gDssGENnHsu5Of7HJjEgxK8+GLIp8/QfPVSANZnVOfy468hJ9R1R0SSZ02NWlx97FXkY17F+PFw333+BiVxoQQvSddl/pec/e24cPmmo/uxqH4jHyMSqdy+aLo/w9qfEqkYPBimTvUvIIkLJXhJql3WruSesZGJZD5o2ZG39j3Sx4hEBODBw3oxa5eWXiEvD846C/7+29+gZJtoshlJmoy8HIa+d09Ul7gGDO5yKZj5HFnwxTpRiJ8TimgyE3/lplfh8hOuZfTzA6izeT0sXsyY9sdxyYmD9DOaohJ6BW9mXc1svpktMLOBxSw3M3s0tPw7MzswVL+rmX1mZj+Y2TwzuzyRcUpyDPrseQ5YPh+AXEvj8uOv1RzvIhXI0joNub7bgHC5+09f0nP2Rz5GJNsiYQnezNKBoUA3oBVwlpm1KrJaN6BF6NUXKBjOLBe42jm3N9AeuLSYbSWFdJ3/BRfMGhUu353Vm1mN9ZWKVDRj9+rASwd0D5dvnPAM+65Y4GNEUl6JvIJvByxwzi10zm0BRgI9iqzTAxjhPNOAuma2s3NuuXPuawDn3DrgB0CtsFJUk79/594xj4TLH7doz/CDTvIxIhEpze2dLuSHBk0BqJaXw5Pv3UXdjWv9DUrKLJEJvhGwJKq8lH8n6a2uY2ZNgQOA6XGPUBJv40aefO8uam/xxrn+rU5Drul+hZ7piVRgmzOqccmJg1hb1RtVctc1f/DIB/eTlp/nc2RSFuacS8yOzU4Dujjn+oTK5wDtnHP9o9YZDdzlnJsSKk8ArnPOzQqVM4FJwB3OuXdK+Jy+eLf3adiwYZuRI0fG7f+QnZ1NZqaeERco8/Fwjr3uv5+dx4wBIK9KFd686W5WNdsjQREmT8Ma8MdGv6Mon9aN6pS4bM6yNTGtV1Q8flZK+uzo+lSQyudGUc1mTee4h+4Kl7866Qymn3JW0s+NIIn38ejUqdMs51zb4pYlshX9UmDXqHJj4PdY1zGzDOBt4JWSkjuAc+5p4GmAtm3buqysrG0OvMDEiROJ5/5SXZmPxxNPQCi5A9zYuS+vZLeEOfGPLdmubp3LA3NSsxPKop5ZJS7rHd2SvZT1iorHz0pJnx1dnwpS+dz4l6odWNz+NC6d9iYA7d59nSfTW9L/zZtj3oV+jxaWzOORyFv0M4AWZtbMzKoCZwKjiqwzCjg31Jq+PbDGObfczAx4FvjBOfdgAmOURJk0CS6PdH54Z59OvLJ/Nx8DEpHyeKBjLz5vsn+4/PAH98Mvv/gYkcQqYQneOZcLXAaMw2sk94Zzbp6Z9TOzfqHVxgALgQXAM8AlofoOwDlAZzObHXp1R1LD4sVw6qmQGxlnflCXy/TcXSQF5aelc/kJ17KsVgMAr4/8iSfCunU+RyZbk9D7SM65MXhJPLpuWNR7B1xazHZTAGWDVLRhA5x0Evz5p1fecUf6nXQDmzOq+RuXiJTb6pp1uOTEgbzx6vVUy8uFuXPh7LPhvfcgPd3v8KQEGqpW4sc56NMHvvHmdycjA95+m+W1G/gbl4hss2932YtBXftHKj78EAYN8i8g2SoleImf226D116LlB97DA47zL94RCSu3tn3SIYdHDUpzX33wQsv+BaPlE4JXuLj5Zfhppsi5f/+13uJSKDce/i5cPzxkYq+fWHKFP8CkhIFpC+H+GryZLjggkj5yCO9q3ep8JqW0gVNk79IcfLT0uGVV6BDB5gzB3JyvHY306ZB8+Z+hydRdAUv22b+fK9FbU6OV27VCt56y3v+LiLBVKsWfPABNAi1r/nzT+jSBVau9DcuKUQJXspv1Sro3j0yZ3TDht7ANnXr+huXiCRekybw/vtQvbpX/uUXOO44WL/e37gkTAleymfdOi+5L1zolWvU8P6ib9LE37hEJHkOOcRrWJsWSiUzZsDpp4fHwBB/KcFL2W3e7N2WnznTK5vBq6/CQQf5G5eIJN+JJ8LQoZHymDHQr5/XbVZ8pQQvZZOb6w1w8emnkbonn/R+yEWkcurXDwYPjpSffRZuuMG/eARQK3opC+e8H+R3oub+ueMOdYcTEW8cjKVL4cUXvfLdd3uN8Q491N+4KjFdwUtsnGP3p57y/jIvcOWVGslKRDxm8MwzXkO7AoMH0/itt/yLqZJTgpetcw4GD2a311+P1J13Htx/vyaQEZGIjAx4801vLIyQPYYO9RK/JJ0SvJTOORgyBO66K1LXowcMHx5pOSsiUqB6da/7XIcOkbr//tcb7VKSSr+hpXQ33+w9Zy9w3HHwxhtQRc03RKQE220Ho0dDmzZe2Tk499zI83lJCiV4KZ5zcMstcOut4aq/Dj7YG6WualUfAxORlFCnDowbR3azZl7ZOTj/fN2uTyIlePk35+D6672r9wLdujHv1luhmuZ1F5EYbb893z74IPznP17ZOW9ymscf9zeuSkL3WaWwvDy4+OLCf2V36QLvvEP+tGklblbapCVSsZTnu0r2xDM6nyq2spwPOXXreuNmHHMMzJrlVfbv7w2YdfXViQyz0tMVvERs2QI9exZO7j16wHvvRcabFhEpq/r1YcIEb2jbAtdc43Wz1Yh3CaMEL571670pH6O7wp17rvfMXcldRLZV6Jk8hx8eqbv7bujdOzIbpcSVErzAihVwxBHeGNIF+veH559Xa3kRiZ9ateCjj+D44yN1I0bACSdAdrZ/cQWUEnxlN28eHHxw5NkYeP3eH3lE/dxFJP5q1vSGu+7TJ1I3dix07gx//OFfXAGk3+CV2SefeONE//abV05Ph2HDvK5xGqFORBKlShV4+mnvYqLAjBnejJSzZ/sXV8AowVdGznmJvFs3WLvWq8vM9OZz18QxIpIMZt7FxJNPRu4WLlnijYAXPaGVlJsSfGWzeTNcdJHXFS4316vbZRf4/HMv4YuIJFO/fvDhh1C7tlfesAFOOcUbQVMt7LeJEnxlsmyZ15gueka4/feH6dO9f0VE/NCtG0ybBs2bR+r+9z849VRYs8a/uFKcEnxlMWmSNy709OmRup494YsvoHFj/+ISEQHYe2/v91PnzpG6d97xfm/puXy5KMEHXW4u3HRT4Raq6enw8MPw0ktei1YRkYpg++29FvUDBkTqfvkF2rf3BuDSLfsyUYIPsiVLvMR+662Qn+/V7bADjB8Pl1+ulvIiUvFkZHjddF9/3es3D17bob59vbuOf//tb3wpRAk+qN55x5vg4fPPI3VZWd6trk6dfAtLRCQmp58OM2dC69aRutdeg/3284a9la3SMGVB89dfcNllMHJkpC4tzZv6ddAg7/Z8GRVMLHF161x6J3nSEal4SppopLQJYkpaT+eQlGrPPb3GdwMGRBoHL10KRx0FV14Jd96pobRLoSv4IHn3XWjVqnBy3203mDzZa5FajuQuIuKrmjVh+HB4+23vGX2Bhx6CAw6AKVP8i62CU4IPguXL4ayz4OSTYeXKSP1553m35Dt08C82EZF4OPlkmDsXuneP1P34I3Ts6PWl/+cf/2KroJTgU1luLjz6KLRsWfiqfeedvYEjXngB6tXzLTwRkbjaaSfvd9uwYd7omwWeesrrZvfmm2ppH0UJPlVNneqN23z55ZHhZsGb4nXePDhWzzZFJIDMvCG1v/++8Kx0K1Z4DfOOPhrmzPEvvgpECT7V/Pqrdzv+0EMLD/6w117e5DEvvqirdhEJvl13hfff967ad9opUj9hgjcy58UXw6pV/sVXASjBp4q//oKrrvISefTt+Bo1vJak334LRx7pX3wiIslm5g1n+8MPcOmlkYbE+fnebfwWLeCee2D9en/j9IkSfEW3Zo036ULz5l6r0ZycyLJTTvFuUw0aBNWq+RejiIif6taFxx/3LnSOPjpSv2YNDBwIu+/uDZ6zaZN/MfpACb6i+vtvuPlmaNrU6+IWPeHCoYd6Y8i/9Za3XEREYJ99YNw4b+rrFi0i9StXwhVXwB57eNPTbtzoX4xJpARf0SxfDoMHe4n7llsKd/1o0cIboW7KFC/Ji4hIYWZw3HFeY+Onnio8mdayZXDJJd7v1zvugNWrfQszGZTgK4pZs+Ccc6BJE++ZenTL+BYtvC5v8+bBSSdpDHkRka3JyPDGr//5Z+/2fMOGkWUrV3p3RnfbzRsR79df/YszgZTg/bRpk9dg7vDDoW1bePnlws/YW7b06r7/3hu0JiPDv1hFRFJR9ereULe//AIPPlj4in79em9mzebNvQF0Ro2CvDz/Yo0zJXg/fPed13+9USOvy1v0hDDgjTz35pveqE09e0IVTRkgIrJNttvOu1pfuBBGjIB9940scw4++gh69IBmzbwZOBct8i3UeFHmSJZly7yk/cor3gxJRVWpAmec4SX+gw6Kebd+TtxR2uQiItF0rgTX1iYfKpikKpbfTyWdJ3H93ZaR4T0O7dXLS+qPPeY1zCsYAW/JErjpJu912GHeeqedBvXrxy+GJFGCT6QVK7yW7q+/XvKECE2awAUXwIUXelf0IiKSeGbebfnu3b2r+mee8Wasix4cZ8oU79W/P3Tr5o2Hf9xxhSe9qcCU4OPJOW90udGjvdf06cWPi1y1qneiXHghdO7sTecqIiL+2H13uOsur2vyu+/CSy95V/UFz+Nzcrzn86NGeYPpdOwIJ54IJ5zg3dKvoJTgt9Uff8CkSTB+vJfUly8vfr20NOjUybsNf8opKXm7R0Qk0KpVgzPP9F4rV3p3X19+Gb76KrJOXh5MnOi9rrjC6+V01FHeADudOnmD7lQQSvBltWKFN796wRf8ww8lr5uW5v2ld8YZ3hV7dDcNERGpuHbc0bs137+/19Xu3Xfhvfdg2rTCd2Z//tl7Pfmk9zv/oIO8nlGHHuq9dtzRt/+CEnwp0jds8JL4jBneX3AzZsDixaVvVL8+dO3qzebWpUvKPKsREZEStGgB113nvZYv90bKe+89Lz9Ej4qXn+89mp0+PVK3xx6RZN+2LRbdFTrBlOCLc+218NFHHPb991ufWzgjAw4+GLKyvMR+8MHq1iYiElQ77+wNoNO3rzeWyZdfeo9oP/nEG7CsaM5YsMB7jRgBQMeMgyYSUwAABjpJREFUDG/8k5NPTnioykTF+fFHmDePYseLq17dG5SmUycvqbdvDzVrJjlAERHxXfXqXkPpzp29Rnp//eU9wv3yS+81cyZs2VJok7ScHK9RXxIowRenXTv48ENcWhq2777eM5V27bx/991XI8qJiMi/bb+9N5z4SSd55c2b4euvvcnBvvoKZs0ib+lS0vfZJynhJDTBm1lX4BEgHRjunLu7yHILLe8ObAB6O+e+jmXbhOrVCzp14vN16zi8W7ekfayIiARItWpwyCHeK+TL0aO92/RJkLAO2GaWDgwFugGtgLPMrFWR1boBLUKvvsCTZdg2cZo1g8MOI79GjaR9pIiIBF/edtsl7bMSOcJKO2CBc26hc24LMBLoUWSdHsAI55kG1DWznWPcVkREREqQyATfCFgSVV4aqotlnVi2FRERkRIk8hl8cY3Qi/Y5K2mdWLb1dmDWF+/2PkC2mc2POcKt2wH4M477Syi7J7H7H5BixyORdCw8UedcuY5Hos9ZP1Tmc6O477PgeGzLdx2w8yTe50eTkhYkMsEvBXaNKjcGfo9xnaoxbAuAc+5p4OltDbY4ZjbTOdc2EftORToeEToWhel4ROhYFKbjUVgyj0cib9HPAFqYWTMzqwqcCYwqss4o4FzztAfWOOeWx7itiIiIlCBhV/DOuVwzuwwYh9fV7Tnn3Dwz6xdaPgwYg9dFbgFeN7nzS9s2UbGKiIgETUL7wTvnxuAl8ei6YVHvHXBprNv6ICG3/lOYjkeEjkVhOh4ROhaF6XgUlrTjYW5rY62LiIhIyknkM3gRERHxiRJ8EWZ2mpnNM7N8M2sbVd/UzDaa2ezQa1hp+wmKko5HaNkgM1tgZvPNrItfMfrFzG42s2VR50R3v2NKNjPrGvr+F5jZQL/j8ZuZLTKzOaHzYabf8SSbmT1nZivNbG5UXX0zG29mP4f+rednjMlSwrFI6u8MJfh/mwucDEwuZtkvzrn9Q69+SY7LL8Uej9DQwf9v7+5B7CrCOIw//8qIxkKMIY2wijZaREELxRBF1C6mEFKZwkIhqazsBBEUQSwEU9gkzRrSBBcFESNoYaEERfyo/ECCISmsXdjwWpyzsGY/3L26Z3bPeX6w3NkLF16GYV5mzpx3jgH3Ak8D7/Ylhqfm7RVjovWZkUE1Lym9cz3Wj4cpvhp2mm4+WOll4EJV3Q1c6P+fgtOs7gsYcM4wwV+nqn6qqv+zWM6utkF/HAHOVtViVf1K9ybEQ8NGp8YsKa1/qKovgD+v+/oIcKZvnwGeGTSoRtbpi0GZ4LdmLsk3ST5P8mjrYBqznHDnZJLv+u24SWw9ruAYWK2AT5Jc7KtsCvb39U3oP29vHE9rg80Zk0zwST5N8v0afxutPi4Dd1TV/cBLwHySW4aJeHvN2B+bLie8m/1L35wC7gIO0o2Pt5oGO7xJjIEteqSqHqB7bHEiyaHWAWlHGXTO2Nb34Heqqnpiht8sAot9+2KSn4F7gF1/kGaW/mBzpYh3vc32TZL3gA+3OZydZhJjYCuq6o/+82qS83SPMdY6zzMlV5IcqKrL/W2hV1sH1EpVXVluDzFnTHIFP4sk+5YPkSW5k+4O+1/aRtXUAnAsyQ1J5uj646vGMQ2qn6yWHaU7kDgllpReIclNSfYut4Enmd6YWMsCcLxvHwc+aBhLU0PPGZNcwW8kyVHgHWAf8FGSb6vqKeAQ8GqSJeAa8GJVNT1AMYT1+qMvO3wO+BFYAk5U1bWWsTbwZpKDdNvSvwEvtA1nWJaUXmU/cD4JdHPrfFV93DakYSV5HzgM3JbkEvAK8AZwLsnzwO/As+0iHM46fXF4yDnDSnaSJI2QW/SSJI2QCV6SpBEywUuSNEImeEmSRsgEL0nSCJngJUkaIRO8JEkjZIKXNJMkD/aXZuzpq7j9kOS+1nFJ6ljoRtLMkrwG7AFuBC5V1euNQ5LUM8FLmllfg/5r4C/g4QmWK5Z2LLfoJf0XtwI3A3vpVvKSdghX8JJmlmQBOAvMAQeq6mTjkCT1vE1O0kySPAcsVdV8f5Xyl0ker6rPWscmyRW8JEmj5DN4SZJGyAQvSdIImeAlSRohE7wkSSNkgpckaYRM8JIkjZAJXpKkETLBS5I0Qn8Djtb5mAg3tkkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_misspecified_mle()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <img src=\"https://img.icons8.com/bubbles/50/000000/video-playlist.png\" style=\"height:50px;display:inline\"> Recommended Videos\n",
    "----\n",
    "#### <img src=\"https://img.icons8.com/cute-clipart/64/000000/warning-shield.png\" style=\"height:30px;display:inline\"> Warning!\n",
    "* These videos do not replace the lectures and tutorials.\n",
    "* Please use these to get a better understanding of the material, and not as an alternative to the written material.\n",
    "\n",
    "#### Video By Subject\n",
    "\n",
    "* Point Estimation - <a href=\"https://www.youtube.com/watch?v=Ieicfj6LYyQ&t=181s\"> MathNStats - Point Estimates </a>\n",
    "* Maximum Likelihood Estimation (MLE)\n",
    "    * Simple Version (6 min) - <a href=\"https://www.youtube.com/watch?v=XepXtl9YKwc\">StatQuest</a>\n",
    "    * Complete Lecture (50 min) - <a href=\"https://www.youtube.com/watch?v=RIawrYLVdIw&t=2263s\">Cornell CS4780</a>\n",
    "* Evaluating Estimators - Bias & MSE - <a href=\"https://www.youtube.com/watch?v=XqWfeND04vs\"> Actuarial Education </a>\n",
    "* Entropy, Cross-Entropy, KL-Divergence - <a href=\"https://www.youtube.com/watch?v=ErfnhcEV1O8\">by Aurélien Géron </a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## <img src=\"https://img.icons8.com/dusk/64/000000/prize.png\" style=\"height:50px;display:inline\"> Credits\n",
    "---\n",
    "* Examples, exercises and definitions from <a href=\"https://probabilitycourse.com/\">Introduction to Probability, Statistics and Random Processes</a> - https://probabilitycourse.com\n",
    "* Icons from <a href=\"https://icons8.com/\">Icon8.com</a> - https://icons8.com\n",
    "* Datasets from <a href=\"https://www.kaggle.com/\">Kaggle</a> - https://www.kaggle.com/"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}