Notebooks/Outlier_Correlation_StatisticalAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objectives:\n",
    "\n",
    "By the end of this lesson, students will be able to:\n",
    "\n",
    "- Describe an outlier and how to remove outliers for univariate data\n",
    "- Visualize the outliers through box-plots\n",
    "- Describe how to use a Gaussian Mixture Model to fit a data distribution\n",
    "- Understand correlation and visualize the correlation amongst features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outlier Detection\n",
    "\n",
    "- **Outliers** are extreme values that can skew our dataset, sometimes giving us an incorrect picture of how things actually are in our dataset. \n",
    "\n",
    "- The hardest part of this is determining which data points are acceptable, and which ones constitute as \"outliers\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activity: find and remove outliers if our dataset is Normal\n",
    "\n",
    "When our sample data is close to normal distribution, the samples that are _three standard deviations away from the mean_ can be considered outliers.\n",
    "\n",
    "**Task**: Write a function that first finds outliers for a normally distributed data, then remove them.\n",
    "\n",
    "**Hint**: Data samples that are below `mean - 3*std` and above `mean + 3*std` are outliers for Normal distribution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def find_remove_outlier(data_sample):\n",
    "    # calculate summary statistics\n",
    "    data_mean, data_std = np.mean(data), np.std(data)\n",
    "    # define cut-off\n",
    "    cut_off = data_std * 3\n",
    "    lower, upper = data_mean - cut_off, data_mean + cut_off\n",
    "    # identify outliers\n",
    "    outliers = [x for x in data_sample if x < lower or x > upper]\n",
    "    # remove outliers\n",
    "    outliers_removed = [x for x in data_sample if x > lower and x < upper]\n",
    "    return outliers, outliers_removed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interquartile range (IQR)\n",
    "\n",
    "We use IQR for finding and removing outliers when the data has any kind of distribution\n",
    "\n",
    "[John Tukey](https://en.wikipedia.org/wiki/John_Tukey) suggested to _calculate the range between the first quartile (25%) and third quartile (75%) in the data_, which is  the IQR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activity:  IQR outlier detection and removal\n",
    "\n",
    "**Task**: write a function to find and remove outliers based on IQR method for this data sample:\n",
    "\n",
    "**Hint**: \n",
    "\n",
    "$Q_1$ is the first quartile (25%) \n",
    "\n",
    "$Q_3$ is the third quartile (75%) \n",
    "\n",
    "<img src=\"Images/iqr.png\"> \n",
    "\n",
    "`x = [norm.rvs(loc=5 , scale=1 , size=100), -5, 11, 14]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n",
      "1.6503050764633622\n",
      "[-5.0, 11.0, 14.0]\n",
      "1.6503050764633622\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "def find_remove_outlier_iqr(data_sample):\n",
    "    # calculate interquartile range\n",
    "    q25, q75 = np.percentile(data_sample, 25), np.percentile(data_sample, 75)\n",
    "    iqr = q75 - q25\n",
    "    print(iqr)\n",
    "    # calculate the outlier cutoff\n",
    "    cut_off = iqr * 1.5\n",
    "    lower, upper = q25 - cut_off, q75 + cut_off\n",
    "    # identify outliers\n",
    "    outliers = [x for x in data_sample if x < lower or x > upper]\n",
    "    # remove outliers\n",
    "    outliers_removed = [x for x in data_sample if x > lower and x < upper]\n",
    "    return outliers\n",
    "\n",
    "y = np.array([-5, 11, 14])\n",
    "x =  np.concatenate((scipy.stats.norm.rvs(loc=5 , scale=1 , size=100), y))\n",
    "print(type(x))\n",
    "print(find_remove_outlier_iqr(x))\n",
    "print(scipy.stats.iqr(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How we can visualy see the outlier?\n",
    "\n",
    "**Answer**: Box plot use the IQR method to display data and outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-2c7f3ccdb2c3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mboxplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.boxplot(x)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation \n",
    "\n",
    "**Correlation** is used to _test relationships between quantitative variables_\n",
    "\n",
    "Some examples of data that have a high correlation:\n",
    "\n",
    "1. Your caloric intake and your weight\n",
    "1. The amount of time your study and your GPA\n",
    "\n",
    "**Question**: what is negative correlation?\n",
    "\n",
    "Correlations are useful because we can find out what relationship variables have, we can make predictions about future behavior. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activity: Obtain the correlation among all features of iris dataset\n",
    "\n",
    "1. Review the iris dataset. What are the features? \n",
    "1. Eliminate two columns `['Id', 'Species']`\n",
    "1. Compute the correlation among all features.\n",
    "    1. **Hint**: Use `df.corr()`\n",
    "1. Plot the correlation by heatmap and corr plot in Seaborn -> `sns.heatmap`, `sns.corrplot`\n",
    "1. Write a function that computes the correlation (Pearson formula)\n",
    "    1. **Hint**: https://en.wikipedia.org/wiki/Pearson_correlation_coefficient\n",
    "1. Compare your answer with `scipy.stats.pearsonr` for any given two features\n",
    "\n",
    "<img src=\"./Images/iris_vis.jpg\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8717541573048714\n",
      "(0.8717541573048712, 1.0384540627941809e-47)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAFBCAYAAAD5QpGoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XeYFFXWx/HvzyGI5JxVBEyggiAqBjAQVcSsq6BrQEysuuKKuiuLAXUVZRVEXHQxB9aAShADBnIQASPIa0CSgCRBYKbP+0fVDD3DDDQwMzXdfT48/UzXrVtVp5qZPn1v3b4lM8M555xLRntFHYBzzjm3uzyJOeecS1qexJxzziUtT2LOOeeSlicx55xzScuTmHPOuaTlScw551zCJD0taYWk+QWsl6R/S1ooaa6kI+PWXSppQfi4tDDi8STmnHNuV/wX6LyD9V2ApuGjF/AEgKRqwF3A0UAb4C5JVfc0GE9izjnnEmZmnwCrd1DlTOBZC0wFqkiqC3QCJpjZajP7DZjAjpNhQjyJOeecK0z1gZ/jlheHZQWV75FSe7oDt+u2rlzkc30Bf2p1U9QhlBjPvXZJ1CGUGLZuVdQhlCjlTu2tPd1Hou85ZWo2vpqgCzDbcDMbvouHyy9e20H5HvEk5pxzqS6WlVC1MGHtatLKazHQMG65AbAkLG+fp3ziHh7LuxOdcy7lWSyxR+EYDfQMRykeA6w1s6XAeKCjpKrhgI6OYdke8ZaYc86lulihJSgkvUTQoqohaTHBiMPSAGY2DBgDdAUWAhuBP4frVku6G5gR7mqAme1ogEhCPIk551yKs6zMwtuX2UU7WW/AdQWsexp4utCCwZOYc86lvsLrKixxPIk551yqS3BgRzLyJOacc6nOW2LOOeeSViEO7ChpPIk551yKK8yBHSWNJzHnnEt13p3onHMuafnADuecc0nLW2LOOeeSlg/scM45l7S8Jeaccy5ZWdbWqEMoMp7EnHMu1XlLzDnnXNLya2LOOeeSlrfEnHPOJS3/nphzzrmk5dNOOeecS1renehSxZ33DeKTSdOpVrUKbz4/LOpwisWf+1/FkSe1YvOmzQy5ZTD/N3/RdnUu6nsJJ559EhUql6fHoRfmlB/S5lAuu+tK9jt4fx694SGmjplcnKEXqklzv+OB58YQi8U4q30rrjijXa71S1eu4c7h/2P9xj+IxWL85fyOnNDiIN6dNIeRYz7Lqffdz8t5+e5rOXi/usV9CoVm0pc/8OCoicFrcVxzLu/YJtf6pavX8fdnx7N+02ZiMaPPmcdzQvNGTPn6R/791mdszcqidEYGN511Am0O2jeis9gFKTywY69EKkm6Q9KXkuZKmiPp6MIKQFJ7Se+Ezy+T9Hhh7TufY1WRdG1+x86nbmlJ90taIGm+pOmSuhRVbMWle9cODBt0T9RhFJuWJ7WibqO63NCuN0/2G8JV91yTb72Z70+n35m3bFe+cslKhvx1MJ+99UlRh1qksmIx7hv5NkP79uSNB/owbso8vv9lRa46T701kU5tmvPqPdfxwHUXcN/ItwE47bgWvHrv9bx67/Xc2/tc6tWoktQJLCsWY+CrHzLkuu68/vdLGTfzW75fuipXnafGTaPjkQfySr9LuP/yrtz3yocAVK1QjsG9z2TUHT25u2cn7hg5LopT2HWxWGKPJLTTlpikY4HTgSPNbLOkGkCZIo+saFQBrgWGJlD3bqAu0Dw879pAu51sU+K1bnEYvyxdHnUYxeaoDm34+H8fAbDg8+8oX6k8VWpVZc2K33LVW/D5d/lu/+vi4I3ekvQPPNv87xfTsHZ1GtSqBkDnYw5j4qyvaVy/1rZKgg1/bAZgw8Y/qFml4nb7GTtlLl2OPbxYYi4q839YRsOaVWhQowoAnVodxMS539O4bvWcOkL8/scWADZs2kzNyuUBOLjhttercd3qbMnMYsvWTMqULtmdWmapO7AjkZZYXWClmW0GMLOVZrZEUitJH0uaJWm8pLoAkiZKelTS5LAF0yYsbxOWfR7+PCjRICV1lDRF0mxJr0mqEJb/IOmfYfk8SQeH5TUlTQjLn5T0Y5h87wcah63Jf4W7ryBplKRvJL2gwD7AVcANcee93MxeDfe/QdID4bm/H57bREmLJHVL9Lxc0atWpzqrlqzMWV61bCXValffwRapacVv66hTrXLOcq1qlVj+27pcda45+xTenfQFHfo8yHUPPcttPU/fbj/jp82j8zHJncRWrNlAnarbEnTtKhVYsWZDrjq9TzuGd2d8Tcc7nuL6oW9y2/knbbef9z9fwMENapb4BAakdEsskST2HtBQ0neShkpqJ6k08Bhwrpm1Ap4G7o3bpryZtSVo9Twdln0DnGhmLYF/APclEmCYfO4ETjWzI4GZwM1xVVaG5U8A2f1BdwEfhuVvANmd1rcB35tZCzPrG5a1BG4EDgUOAI4DmgA/mVnuv/K48wMmhue+HrgH6ACcBQxI5Lxc8ZC0faFZ8QcSsfxOOe9rM3bKXLqd0JIJ/76VIbf05I5ho4jFvbHNXfgze5cpQ9OGtYs63CKV3/9+3tdi3Mxv6XZ0M9679yoev7Y7d44cRyy2bcuFS1Yy+K3PuPOiU4s42kKSlZnYIwnt9COEmW2Q1Ao4ATgJeIXgTbs5MCH8z88AlsZt9lK47SeSKkmqAlQERkpqSvB7VDrBGI8hSDCTwmOVAabErX89/DkLODt8fjxBQsHMxknK3XeU23QzWwwgaQ6wPzB3JzFtAbI7w+cBm81sq6R54fbbkdQL6AUw9OF7uLLnRTs5hNtdnXp25dQLOwCwcO5CqterkbOuep0arF6xOqrQIlO7WiWWrV6bs7xi9Tpq5ekufOPjWTzRtycARzTdl81bM/lt/UaqV64AwPip8+hy7GHFF3QRqV2lAst+W5+zvHzNhpzuwmxvTJ7P0OuDt5MjDqjH5q2ZrPl9E9Uq7sPy39Zz81Nvc3fPTjSsWaVYY99tKTw6MaGBHWaWZWYTzewu4HrgHODLsEXTwswOM7OO8Zvk3QXBNaaPzKw5cAawd4IxCpgQd6xDzeyKuPWbw59ZbEvK+Xz8LtDmuOfZ+1gI7Ctp+4sCga1mOZ9tY9n7MLMYBXwwMLPhZtbazFp7Aita458dQ9+uN9G3603MeG8q7c4JuoKatjyQjet/3+56WDpodkB9flq2isUrVrM1M5NxU+fR7siDc9WpW70y074MRm4u+mUFW7ZmUq1S8OYei8V4b/r8pO9KBGi2Xx1+WvEbv6xcy9bMLMbP+pZ2hx2Qq07dapWY9s1PACxatootmVlUrVCOdRv/4IYn3qRPt+Np2bh+FOHvnhTuTkxkYMdBQMzMFoRFLYCvgY6SjjWzKWH34oFm9mVY5wLgI0nHA2vNbK2kysAv4frLdiHGqcAQSU3MbGF4vaqBmeV/JT7wGXA+8ICkjkDVsHw9QYtwh8xso6QRwL8lXW1mW8JrfqeY2fO7EHuJ0/eu+5nx+VzWrFnHKd0v4dorenDOGZ2iDqvIzP5wFi1Pas1jnwxjy6bNDLnlsZx1/xrzCH273gTAJf0u5fgzT6RMubIMmzqCD16ewGuPvkzjw5vQd3g/yleuQKtTj+L8my7i5g43RHU6u61URgb9ep7ONf8aSSwWo/uJrWjSoDZD/vc+zRrVp/2Rh/DXP3VhwIg3eX7cZCQY0OvsnG62Wd/+QO1qlXIGhiSzUhl7cdv5J3PNkNeJxYwzj21Gk3o1GPrOZA7dtzbtD2/MzWefyIAXJ/DCR7MB8c8enZDEKx9/wU+/rmH42GkMHzsNgGE3nE21ivtEe1I7k8ItMdlOrg+EXYmPEYzsyyRopfQCGgD/BioTJMNHzewpSRMJuvvaAZWAy81sejjKcSTwK/Ah0MPM9pfUHrjFzE6XdBnwOLAmLoRjgAOBB4CyYdmdZjZa0g9AazNbKak18JCZtZdUi6BLsyrwMUFSbRSOMnwROBwYC7ybfezwXB8HZprZfyWVIeg2PRv4A/gd+IeZjZe0wcyyB5f0BzaY2UPhcs66gmxduSj9Lsrk40+tboo6hBLjudcuiTqEEsPWrdp5pTRS7tTeu9KzlK9NY/+d0HtOuS599vhYxW2nSWyXdxgksVvMbGah7njXYigLZJlZZpg8nzCzFlHFk5cnsYAnsW08iW3jSSy3Qkli7z6aWBI77cakS2JJMDZ0t+wLvCppL4JBGFdFHI9zzkUnSUceJqLQk5iZtS/sfe5GDAsIhs4755xL4WtiCY1OdM45l8QKcXSipM6SvpW0UNJt+ax/JJxQYk74/eI1ceuy4taNLoxTS9XuROecc9kKqSUmKQMYQjC5w2JghqTRZvZVzqHMboqrfwO5e8U2Ffb4BG+JOedcqiu8llgbYKGZLTKzLcDLwJk7qH8R4eQXRcWTmHPOpbqsrMQeO1cf+DlueXFYth1J+wGNCL5SlW1vSTMlTZXUfXdPJ553JzrnXKpL/HpXzvR4oeFmNjy+Sj6bFTR8/0JglOWeQn/fcAL5A4APJc0zs+8TCq4AnsSccy7VJZjEwoQ1fAdVFgMN45YbAEsKqHshcF2e/S8Jfy4Kv1PcEtijJObdic45l+oslthj52YATSU1Cmc1uhDYbpRhOF1hVeIma5dUNZyIIvvuJMcBX+Xddld5S8w551JdIU3uG86CdD0wnuDuJU+b2ZeSBhBM2Zed0C4CXrbcU0IdAjwpKUbQgLo/flTj7vIk5pxzqa4Qpxc0szHAmDxl/8iz3D+f7SYDhX4vH09izjmX6jJ92innnHPJKoWnnfIk5pxzKc5iqXvjDE9izjmX6pL0rs2J8CTmnHOpzrsTnXPOJS3vTnTOOZe0fHSic865pFWI3xMraTyJOedcqvOBHc4555KWXxNzhelPrW7aeaU08eKsR6IOoUQoV++EqEMoMSqV3SfqEEqU1et77/lOfHSic4XPE5hzxcMyE7rhZVLyJOacc6nOuxOdc84lLe9OdM45l7S8Jeaccy5p+RB755xzSctbYs4555JWlo9OdM45l6TMuxOdc84lLe9OdM45l7Q8iTnnnEta/j0x55xzSctbYs4555KVZXpLzDnnXLLy0YnOOeeSlncnOuecS1qexJxzziUrM09izjnnkpUP7HDOOZesLIW7E/eKOgDnnHNFLGaJPRIgqbOkbyUtlHRbPusvk/SrpDnh48q4dZdKWhA+Li2MU/OWmHPOpbpC6k2UlAEMAToAi4EZkkab2Vd5qr5iZtfn2bYacBfQGjBgVrjtb3sSk7fEnHMuxVnMEnokoA2w0MwWmdkW4GXgzATD6ARMMLPVYeKaAHTerROK40nMOedSXYLdiZJ6SZoZ9+iVZ0/1gZ/jlheHZXmdI2mupFGSGu7itrvEk1iK+nP/q3js42E8NG4wjZofkG+di/pewhNTRvDcVy/nKj+kzaE88O4gXv7+dY7p2rY4wo3EnfcN4sTTLqT7Jb2jDqXYPDJoAN989RmzZ02gZYvmO6z7xuvPMOfzD3KWHxh4J/PnfczsWRMY9dp/qFy5UlGHW6QGPvh3Zs55n0+nvM3hRxyab53XXh/BJ5NHM3n6GB5+dAB77bXtLfOqq3swbfZ4Jk8fQ/+7by2usHeLZVpiD7PhZtY67jE8z66U3+7zLL8N7G9mhwPvAyN3YdtdtttJTNIdkr4Ms+0cSUfvaTBx+24v6R0FVkqqGpbXlWSSjo+r+6uk6pJ6S+qZz772lzQ/fN5CUte4df0l3VJADHUkvSzpe0lfSRoj6cDCOsei1PKkVtRtVJcb2vXmyX5DuOqea/KtN/P96fQ7c/vTX7lkJUP+OpjP3vqkqEONVPeuHRg26J6owyg2XTqfTNMmjTj40OO55pq/MeTxgQXW7d69Cxs2/J6r7P0PPuGIFidzZKsOLFiwiNv+dn0BW5d8p3ZsR+PG+9G6xanc1OfvPPzIgHzrXX7pXzixbTfatulKjRrV6H5WFwCOP+Foupx2CicccwZt23Tl8cH/Kc7wd10swcfOLQYaxi03AJbEVzCzVWa2OVx8CmiV6La7Y7eSmKRjgdOBI8Nseyq5m4mFwoJv6E0Djg2L2gKfhz+RdBCwMnzRhpnZszvZZQug607qIEnAG8BEM2tsZocCtwO1d+9MitdRHdrw8f8+AmDB599RvlJ5qtSqul29BZ9/x5oV219T/XXxCn765seUvhssQOsWh1G5UsWowyg2Z5zRiedeGAXAtOmzqVylMnXq1NquXvny+3DTX3px38DBuconvP8JWeFt7qdOm039+nWLPugi0vW0U3n5pTcBmDljDpWqVKR27Zrb1Vu/fgMApUqVonSZ0jlfGr78yj8xeNBwtmzZAsDKlauLKfLdU4jXxGYATSU1klQGuBAYHV9BUvwvRjfg6/D5eKCjpKphw6RjWLZHdrclVpcgeWwGMLOVZrZEUitJH0uaJWl89slImijpUUmTJc2X1CYsbxOWfR7+PCifY00iTFrhz0HkTmqTw33ltKrCOL6QNAW4LiwrAwwALghbjheE+zg0jG+RpD5h2UnAVjMblh2Emc0xs0/DVuLHkl6V9J2k+yVdLGm6pHmSGu/ma1poqtWpzqolK3OWVy1bSbXa1SOMyJUE9evVYfHP2z74/rJ4KfXr1dmu3oD+tzLo0SfZuHFTgfv682UXMm78R0USZ3GoW682v/yyNGd5yS/LqFsv/8+oo954mu8WTWXD+t95681xADRu0ohj27ZmwoejeHvsC7Q88rBiiXu3FVJLzMwygesJks/XwKtm9qWkAZK6hdX6hL10XwB9gMvCbVcDdxMkwhnAgLBsj+xuEnsPaBi+iQ+V1E5SaeAx4FwzawU8Ddwbt015M2sLXBuuA/gGONHMWgL/AO7L51iT2ZbE2gBvsq1J2pYgyeX1DNDHzLKTHeFImn8QDP1sYWavhKsOJhg10wa4KzyP5sCsHZz/EcBfgMOAHsCBZtYG+A9wQ34bxF8wXbThhx3ses8FDck8UnjaGZeY/H4v8k5HdMQRzWjcZH/eemtcgfvpd1sfMjMzefHF1ws9xuKSyGuR7dyzLueQpm0pW7YMJ7YL3lJKlcqgcpXKdDj5XO668wGeHjk4321LCosl9khoX2ZjzOzAsJfq3rDsH2Y2Onzez8yamdkRZnaSmX0Tt+3TZtYkfDxTGOe2W98TM7MNkloBJxC0Wl4B7iF4858Q/oJkAEvjNnsp3PYTSZUkVQEqAiMlNSW4wFc6n8NNB1pKKg+UDo+9SFITgiT2cHxlSZWBKmb2cVj0HNBlB6fzbtii3CxpBYl1Gc4ws6Xh8b4nSOoA8whej+2EF0iHA5y335mFnlE69ezKqRd2AGDh3IVUr1cjZ131OjVYvaJkd3e4onFN70u54oqLAZg5cw4NGtbLWVe/QV2WLF2eq/4xR7fiyJaHsfC7qZQqVYpatarzwYTXOKXDeQD06HEep3U9lQ6dzi++kygkV1x1MT0vCzpgPp89N1d3aL36dVi2dEWB227evIWxYz6gy2mnMPGjSSz5ZRnvjA56wmbPmkssZlSvUY1VJbRb0TKjjqDo7PbADjPLMrOJZnYXQfPyHODLsJXTwswOM7OO8Zvk3QVB0/IjM2sOnAHsnc9xNgILgcuB2WHxVIJrW7WAb/NsonyOtSOb455nEST2L9l2MXJn28TilmNE9AXy8c+OoW/Xm+jb9SZmvDeVducEubRpywPZuP73fK99udT3xLCRtD6qI62P6sjo0ePpcfG5ABzd5kjWrV3HsmW537ifHP4s++7fiiYHHkO7k7rz3YJFOQmsU8f29L3lWrqffRmbNv1R7Oeyp0Y89QLtjutGu+O68e4773PhRd0BaH1UC9atXc/y5b/mql++/D4518kyMjLo0LEdC75bBMC777yf0ypr3GR/ypQpXWITGFCYAztKnN0d2HFQ2HrK1oKgf7RmOOgDSaUlNYurc0FYfjyw1szWApWBX8L1l+3gkJOAG4Ep4fIUgu68qZanD8DM1gBr40YwXhy3ej1B629nPgTKSroqu0DSUZLaJbBt5GZ/OIvlPy3nsU+G0fv+63jqzidz1v1rzCM5zy/pdynDpo6gTLmyDJs6gvNuvBCAxoc3YdjUERxz2nH0uu8aBk14rNjPoTj0vet+Lr76Jn74aTGndL+E/729x9eYS7QxYz9g0f/9xLdfT2LYsAe5/obbc9bNnPHeDrYMDH70HipWqMC4sS8zc8Z7DHn8/qIMt0hNGD+RH374mVlffMCjj91D35v756z7eFIwTmGffcrxwivD+HTK23w65W1+/XUVz4x4CYAXnhvFfvs3ZNK0d/nPM49y7dUlfIh9IXYnljTanSn6w67Ex4AqQCZBS6kXwZDJfxMkp1LAo2b2lKSJBImnHVAJuNzMpocJbyTwK0Hi6GFm+0tqD9xiZqeHxzsPeBVoamYLJZUF1gH9zWxgWKc/sMHMHgrjexrYSHAB8lwza65g2pPxBN2WA4FDsrcJ9zEfON3MfpBUD3iUoEX2B/ADQSKtnye2ieHyzLxxF6QouhOT0YuzHtl5pTRRrt4JUYdQYlQqu0/UIZQoq9cvyO/7VbtkxSntEnrPqfXBx3t8rOK2W0lslw8S90Zf5AdLAp7EAp7EtvEkto0nsdwKI4ktPymxJFb7o+RLYj4BsHPOpTpLutyUsGJJYmbWvjiO45xzbnuxTE9izjnnklSyDtpIhCcx55xLcebdic4555KVt8Scc84lLYt5S8w551ySSuWpUz2JOedciotlpu79jz2JOedcivOWmHPOuaTl18Scc84lLR9i75xzLmn5EHvnnHNJKyvmAzucc84lKb8m5pxzLmn56ETnnHNJy1tizjnnklbMRyc655xLVj7E3jnnXNLK8u5E55xzycpbYs4555KWj050heq51y6JOoQSoVy9E6IOocTYtOTTqEMoMWzT+qhDSDk+sMM551zS8u5E55xzSSuVW2KpO6GWc845ALJMCT0SIamzpG8lLZR0Wz7rb5b0laS5kj6QtF/cuixJc8LH6MI4N2+JOedciius7kRJGcAQoAOwGJghabSZfRVX7XOgtZltlHQN8CBwQbhuk5m1KJRgQt4Sc865FBdL8JGANsBCM1tkZluAl4Ez4yuY2UdmtjFcnAo02PMzKJgnMeecS3GGEnpI6iVpZtyjV55d1Qd+jlteHJYV5ApgbNzy3uF+p0rqXhjn5t2JzjmX4mIJfk/MzIYDw3dQJb9+yXz3LukSoDXQLq54XzNbIukA4ENJ88zs+8Siy58nMeecS3FZhdfpthhoGLfcAFiSt5KkU4E7gHZmtjm73MyWhD8XSZoItAT2KIl5d6JzzqW4QrwmNgNoKqmRpDLAhUCuUYaSWgJPAt3MbEVceVVJZcPnNYDjgPgBIbvFW2LOOZfiLN9ewN3Yj1mmpOuB8UAG8LSZfSlpADDTzEYD/wIqAK9JAvjJzLoBhwBPSooRNKDuzzOqcbd4EnPOuRSXYCsrIWY2BhiTp+wfcc9PLWC7ycBhhRgK4EnMOedSXmEmsZLGk5hzzqW4wupOLIk8iTnnXIrLlCcx55xzSSqFbyfmScw551KdXxNzzjmXtGLeneiccy5ZeXeic865pOXdic4555KWj050zjmXtLw70TnnXNKKpW5DzJNYKpo09zseeG4MsViMs9q34ooz2uVav3TlGu4c/j/Wb/yDWCzGX87vyAktDuLdSXMYOeaznHrf/bycl+++loP3q1vcp1CoHhk0gC6dT2bjpk1cccVNfD5nfoF133j9GRo12pcWLU8B4IGBd3La6R3YsmULixb9yBVX3szateuKK/Ric+d9g/hk0nSqVa3Cm88PizqcIvfZ9Dk8MPQZsmIxzu5yCldelPv+jEuW/8o/HnqC1WvWUbliBQb2u4E6NasDsHT5Su4aNIxlv65CwND7+lG/Tq0IziJxqXxNbJduxSIpS9IcSfMlvSZpn53Uvz3B/f4QTs2PpA27EtOuknSZpHr5HTuful3Cu5B+LekbSQ8VZWyFISsW476RbzO0b0/eeKAP46bM4/tfVuSq89RbE+nUpjmv3nMdD1x3AfeNfBuA045rwav3Xs+r917Pvb3PpV6NKkmfwLp0PpmmTRpx8KHHc801f2PI4wMLrNu9exc2bPg9V9n7H3zCES1O5shWHViwYBG3/e36og45Et27dmDYoHuiDqNYZGXFuPexEQy973beGvEIYz+axPc/Ls5V56Enn+OMDify+lMP0bvHuQwe8WLOutsfeJzLzu/G6Kcf4aUhA6lWpXJxn8IuswQfyWhX7ye2ycxamFlzYAvQeyf1E0pixewyoN7OKklqDjwOXGJmhwDNgUVFG9qem//9YhrWrk6DWtUoXaoUnY85jImzvs5dSbDhj+A+dRs2/kHNKhW328/YKXPpcuzhxRFykTrjjE4898IoAKZNn03lKpWpk8+n5vLl9+Gmv/TivoGDc5VPeP8TsrKyAJg6bTb16yd3Ui9I6xaHUbnS9r8HqWjetwvZt14dGtarTenSpejSvi0fTZqRq86iHxdzdMtgwvU2LZrx0eSZAHz/42KysrJo2yr429in3N6U27ts8Z7AbshUYo9ktCc3xfwUaALBbaglTQ9baU9KypB0P1AuLHshrPempFmSvpTUK9EDSaop6X+SZoSP48Ly/pKeljRR0iJJfeK2+XvYepog6SVJt0g6l+B22S+EcZULq98gabakeZIODstuBe41s28guI+OmQ0N9/1fSU9I+ig8brswjq8l/XcPXtM9tuK3ddSptu2TYa1qlVj+W+7ur2vOPoV3J31Bhz4Pct1Dz3Jbz9O328/4afPofEzyJ7H69eqw+OdtN579ZfFS6ters129Af1vZdCjT7Jx46YC9/Xnyy5k3PiPiiROV3xWrFxNnVrVc5Zr16zO8lWrc9U58ID9eP/TaQB88Nl0ft+4iTVr1/PD4iVUrFCeG/s/xHlX38rDTz5HVlbJ76wrxJtilji7lcQklQK6APMkHQJcABxnZi2ALOBiM7uNbS23i8NNLzezVgSJpI+k6vntPx+DgUfM7CjgHOA/cesOBjoBbYC7JJWW1Dqs1xI4OzweZjYKmBnG18LMst+xVpplq5brAAAgAElEQVTZkcATwC1hWXNg1g5iqgqcDNwEvA08AjQDDpPUIsHzKnSWT5+A8gyvHTtlLt1OaMmEf9/KkFt6csewUcRi236F5y78mb3LlKFpw9pFHW6Ry3vuAJbnRTriiGY0brI/b701rsD99LutD5mZmbz44uuFHqMrXnn//wGUZ5b3W67uwcy5X3He1bcyc+5X1KpRjYyMDLKyYsye9zV/7dWDl4YOZPHS5bz13sRiinz3mRJ7JKNdTWLlJM0hSAQ/ASOAU4BWwIxw3SnAAQVs30fSF8BUoCHQNMHjngo8Hu5/NFBJUnbfx7tmttnMVgIrgNrA8cBbZrbJzNYTJJkdyX5nmgXsn2BMb1vw1zAPWG5m88wsBnyZ3z4k9Qqvr80c8cb7CR5i19WuVollq9fmLK9YvY5aeboL3/h4Fp2Obg7AEU33ZfPWTH5bvzFn/fip8+hybKHfu67YXNP7UmbOeI+ZM95jydJlNGi4rfe4foO6LFm6PFf9Y45uxZEtD2Phd1P5+KM3ObDpAXww4bWc9T16nMdpXU+lR8/UvB6WbmrXrM6yFatylpf/uopa1avmqlOrRjUe7X8Lrz35IH0uvwiAihX2oXaNahzcpBEN69WmVEYGJx/Xhq8WlPirDN4Si5PdsmphZjeY2RZAwMi48oPMrH/eDSW1J0hGx5rZEcDnwN67EOexcceoHyYngM1x9bIIRlzu6meK7H1kbw9BMmqVwDaxPDHEyGfUp5kNN7PWZtb6irPyvfFpoWh2QH1+WraKxStWszUzk3FT59HuyINz1albvTLTvgz+8Bb9soItWzOpVql8EHwsxnvT5yd1V+ITw0bS+qiOtD6qI6NHj6fHxecCcHSbI1m3dh3LluUe6PLk8GfZd/9WNDnwGNqd1J3vFizilA7nAdCpY3v63nIt3c++jE2b/ij2c3GFr/lBjfnxl6UsXrqCrVszGTtxMu3bts5V57e163J6J/7z0huc1fmkcNsmrNvwO6vXBF300+bMp/F+DYr3BHZDKiexwhhi/wHwlqRHzGyFpGpARTP7EdgqqbSZbQUqA7+Z2cbwutMxu3CM94DrgX8BSGphZnN2UP8z4ElJAwnO8TTgqXDdeiCRK9j/Al6X9JmZfSdpL+BGMxu0C3EXu1IZGfTreTrX/GsksViM7ie2okmD2gz53/s0a1Sf9kcewl//1IUBI97k+XGTkWBAr7Nzut1mffsDtatVokGtahGfSeEYM/YDOnc+mW+/nsTGTZu48sqbc9bNnPEerY/quMPtBz96D2XLlmXc2JcBmDZtNtddf1uRxhyFvnfdz4zP57JmzTpO6X4J117Rg3PO6BR1WEWiVEYGt99wOb1vu5esWIyzOp9Ek/0b8vh/X6HZgY05qW1rZnzxFYNHvIgQrQ4/hDtuuAKAjIy9+OvVPbiy7wDMjEMPPIBzuxbdh9LCkqwjDxOh/PqHC6wsbTCzCvmUXwD0I2gxbQWuM7Opkh4AugGzgcuBN4H6wLdATaC/mU2U9APQ2sxWSooBS+J2Pwh4FhgCHEKQlD4xs96S+gMbzOyhMI75wOlm9kO47iLgR+BXYKKZPSXpHOA+YBNwLPB13LFbAw+ZWftwf6cD/wT2Ifg9eNfM+oaDN94xs1GS9g+fNw+3yVlX0Ov4x/TXUvl3KmEVjr8x6hBKjE1LPo06hBLDNq3feaU0UqbhEXt8teqRfS9J6D3npp+eT7orY7uUxJKJpApmtiH8LtsnQC8zmx11XOBJLJsnsW08iW3jSSy3wkhiDyeYxP6ahEkslWfsGC7pUILrbiNLSgJzzrnilsqfmlM2iZnZn6KOwTnnSgKfO9E551zSStaRh4nwJOaccynOuxOdc84lrcwUTmOexJxzLsWlbgrzJOaccynPr4k555xLWqk8OnFPbsXinHMuCcSwhB6JkNRZ0reSFkrabg42SWUlvRKunxbOapS9rl9Y/q2kQpnXzJOYc86luKwEHzsjKYNgCsAuwKHAReGkEvGuIJgntwnBLaoeCLc9FLiQ4JZVnYGh4f72iCcx55xLcYXYEmsDLDSzReFdTF4GzsxT50xgZPh8FHCKghnGzwReDm+d9X/AwnB/e8STmHPOpThL8JGA+sDPccuLw7J865hZJrAWqJ7gtrvMk5hzzqW4RO8nFn/z3vDRK8+u8hsikjf/FVQnkW13mY9OdM65FJfooA0zGw4M30GVxUDDuOUG5L51VnydxZJKEdxLcnWC2+4yb4k551yKK8TuxBlAU0mNJJUhGKgxOk+d0cCl4fNzgQ8tuOfXaODCcPRiI6ApMH33zyrgLTHnnEtxWYU0Z4eZZUq6HhgPZABPm9mXkgYAM81sNDACeE7SQoIW2IXhtl9KehX4CsgkuHlyIoMid8iTmHPOpbjCnLHDzMYAY/KU/SPu+R/AeQVsey9wbyGG40nMOedSXaLXxJKRJzHnnEtxqZvCPIk551zK85aYc865pFVYAztKIk9iEbB1q6IOoUSoVHafqEMoMWzT+qhDKDFUrmLUIaQcvxWLc865pGXeEnPOOZesvCXmnHMuacXMW2LOOeeSVOqmME9izjmX8rJSuEPRk5hzzqW41E1hnsSccy7l+ZednXPOJS0fYu+ccy5peXeic865pGU+xN4551yyyvTuROecc8nKr4k555xLWj460TnnXNLya2LOOeeSlo9OdM45l7R82innnHNJy7sTnXPOJS0f2OGccy5p+RB755xzSctviumccy5ppW4K8yTmnHMpL9NHJ7pkMunLH3hw1ERisRhnHdecyzu2ybV+6ep1/P3Z8azftJlYzOhz5vGc0LwRU77+kX+/9Rlbs7IonZHBTWedQJuD9o3oLArPwAf/ToeO7di0aRPX9f4bc7/4ars6r70+gtp1alKqVCmmTJ5J35v7E4sFf/hXXd2DK6++hKzMLN4bP5H+f3+wuE+hUHw2fQ4PDH2GrFiMs7ucwpUXdc+1fsnyX/nHQ0+wes06KleswMB+N1CnZnUAli5fyV2DhrHs11UIGHpfP+rXqRXBWRS9O+8bxCeTplOtahXefH5Y1OEUilQenbjXnu5AUpakOZLmS3pN0j47qX97gvv9QVINSY9IujGufLyk/8QtPyzpZkn1JI0qYF8TJbXOe3xJ+0uav4MYbpH0TXhuX0jqmUjsUcqKxRj46ocMua47r//9UsbN/Jbvl67KVeepcdPoeOSBvNLvEu6/vCv3vfIhAFUrlGNw7zMZdUdP7u7ZiTtGjoviFArVqR3b0bjxfrRucSo39fk7Dz8yIN96l1/6F05s2422bbpSo0Y1up/VBYDjTziaLqedwgnHnEHbNl15fPB/8t2+pMvKinHvYyMYet/tvDXiEcZ+NInvf1ycq85DTz7HGR1O5PWnHqJ3j3MZPOLFnHW3P/A4l53fjdFPP8JLQwZSrUrl4j6FYtO9aweGDbon6jAKVQxL6JGM9jiJAZvMrIWZNQe2AL13Uj+hJBZnMtAWQNJeQA2gWdz6tsAkM1tiZucmsL9Ek2hvoAPQJjy3EwHtSuBRmP/DMhrWrEKDGlUoXSqDTq0OYuLc73PVEeL3P7YAsGHTZmpWLg/AwQ1rUatKBQAa163OlswstmzNLN4TKGRdTzuVl196E4CZM+ZQqUpFateuuV299es3AFCqVClKlymd88n18iv/xOBBw9myJXi9Vq5cXUyRF6553y5k33p1aFivNqVLl6JL+7Z8NGlGrjqLflzM0S0PA6BNi2Z8NHkmAN//uJisrCzatjocgH3K7U25vcsW7wkUo9YtDqNypYpRh1GoLMF/yagwkli8T4EmAJIukTQ9bKU9KSlD0v1AubDshbDem5JmSfpSUq989jmJMIkRJK/5wHpJVSWVBQ4BPo9vVUkqJ+llSXMlvQKUC8u3Oz6QIemp8PjvSSoXlt8OXGtm6wDMbK2ZjQz384Ok+yRNkTRT0pFhC/H7MPlFZsWaDdSpuu0PsHaVCqxYsyFXnd6nHcO7M76m4x1Pcf3QN7nt/JO228/7ny/g4AY1KVM6uXuc69arzS+/LM1ZXvLLMurWq51v3VFvPM13i6ayYf3vvPVm0Apt3KQRx7ZtzYQPR/H22BdoeeRhxRJ3YVuxcjV1alXPWa5dszrLV+VOyAcesB/vfzoNgA8+m87vGzexZu16fli8hIoVynNj/4c47+pbefjJ58jKSt1rLKnIzBJ67ClJ1SRNkLQg/Fk1nzotwvfOL8P36Avi1v1X0v+F79FzJLXY2TELLYlJKgV0AeZJOgS4ADjOzFoAWcDFZnYb21puF4ebXm5mrYDWQB9J1eP3a2ZLgExJ+xIksynANODYcJu5ZrYlTzjXABvN7HDgXqBVuK/8jt8UGGJmzYA1wDmSKgIVzex7CvazmR1LkLj/C5wLHAPk319VTPL7NZRyNyDHzfyWbkc34717r+Lxa7tz58hxxGLbtly4ZCWD3/qMOy86tYijLXp5zx0Kvj5w7lmXc0jTtpQtW4YT2x0LQKlSGVSuUpkOJ5/LXXc+wNMjBxdpvEUlv3NWno6FW67uwcy5X3He1bcyc+5X1KpRjYyMDLKyYsye9zV/7dWDl4YOZPHS5bz13sRiitwVhmLsTrwN+MDMmgIfhMt5bQR6hu+5nYFHJVWJW983fI9uYWZzdnbAwkhi5STNAWYCPwEjgFMIEseMcN0pwAEFbN9H0hfAVKAhQVLJK7s1lp3EpsQtT86n/onA8wBmNheYu4P4/y/uhZoF7E/Qbbiz/9HR4c95wDQzW29mvwJ/5PkPAUBSr7DVNnPEu5/uZNe7r3aVCiz7bX3O8vI1G3K6C7O9MXk+HVsdCMARB9Rj89ZM1vy+Kaj/23pufupt7u7ZiYY1tzuNpHDFVRfz8aTRfDxpNMuWLqd+/bo56+rVr8OypSsK3Hbz5i2MHfMBXU47BQhabu+MHg/A7FlzicWM6jWqFe0JFIHaNauzbMW2a6PLf11Freq5PyTXqlGNR/vfwmtPPkifyy8CoGKFfahdoxoHN2lEw3q1KZWRwcnHteGrBYuKNX63Z7IsltCjEJwJjAyfjwS6561gZt+Z2YLw+RJgBbB9H3+CCvOaWAszuyFsFQkYGVd+kJn1z7uhpPbAqcCxZnYE8Dmwdz7HyL4udhhBd+JUgpZYW4IEl59EP1ZsjnueBZQKuxB/l1RQ4o3fLpZnHzHyGfVpZsPNrLWZtb7itBMSDG3XNduvDj+t+I1fVq5la2YW42d9S7vDcp9G3WqVmPbNTwAsWraKLZlZVK1QjnUb/+CGJ96kT7fjadm4fpHFWNRGPPUC7Y7rRrvjuvHuO+9zYTgKr/VRLVi3dj3Ll/+aq3758vvkXCfLyMigQ8d2LPgueJN+9533c1pljZvsT5kypVmVhNfFmh/UmB9/WcripSvYujWTsRMn075t61x1flu7LmdE5n9eeoOzOp8UbtuEdRt+Z/WadQBMmzOfxvs1KN4TcHsk0Wti8R+2w0d+l3h2pLaZLQUIf+5wCKukNkAZIL7X696wm/GR8JLRDhXVBY8PgLckPWJmKyRVI+ie+xHYKqm0mW0FKgO/mdlGSQcTdMflZxLwV2CRmWUBq8PWTjPgqnzqfwJcDHwkqTlweNy6+OPvyEBgiKQLzGydpErAhWY2PLGXIBqlMvbitvNP5pohrxOLGWce24wm9Wow9J3JHLpvbdof3pibzz6RAS9O4IWPZgPinz06IYlXPv6Cn35dw/Cx0xg+Nrg2MuyGs6lWcYcDTku0CeMn0qFjO2Z98QGbNm3i+mu29W58PGk07Y7rxj77lOOFV4ZRtmwZMjIy+OTjKTwz4iUAXnhuFI8NHcikae+yZctWrr361qhOZY+Uysjg9hsup/dt95IVi3FW55Nosn9DHv/vKzQ7sDEntW3NjC++YvCIFxGi1eGHcMcNVwCQkbEXf726B1f2HYCZceiBB3Bu1+Tvai5I37vuZ8bnc1mzZh2ndL+Ea6/owTlndIo6rD2S6Iwd4fvbDt/jJL0P1Mln1R27EpOkusBzwKVmOc3AfsAygsQ2HPgbO7lEoz29mCdpg5lVyKf8gjCgvYCtwHVmNlXSA0A3YDZwOfAmUB/4lqBJ2d/MJkr6AWhtZislZQC/Af82szvD/f+XoAV3ULi8P/COmTUPB2c8AxwKzCEYbNLHzGbmOf4d2duE+7gFqGBm/RVcTOkLXBHGvxV42MyezxPbZeHz68N95Kwr6DXb9P6w5BwGVMjqn/Vw1CGUGMu+yvfbIWlJ5VJrZOCeKl3jgD0eFd2s9tEJved8uXzaHh1L0rdAezNbGiapidnv0XnqVQImAgPN7LUC9tUeuMXMTt/hMVP5S3AllSexgCexbTyJbeNJLLfCSGKH1GqT0HvO1yum72kS+xewyszul3QbUM3Mbs1TpwwwFnjbzB7Ns65umAAFPAL8EQ7IK1BhD7F3zjlXwhTjwI77gQ6SFhB8z/Z+AEmttW2SivMJBt9dls9Q+hckzSMYMFcD2Om3zpP7S0DOOed2qri+yGxmqwhGo+ctnwlcGT5/nnD0eD71Tt7VY3oSc865FOe3YnHOOZe0knVKqUR4EnPOuRRnhXO9q0TyJOaccykuWWeoT4QnMeecS3GFNPKwRPIk5pxzKS6Vvw/sScw551Kcj050zjmXtHx0onPOuaTl3YnOOeeSlo9OdM45l7SyYj460TnnXJLy7kTnnHNJy7sTnXPOJS1viTnnnEta/j0x55xzScunnXLOOZe0vDvROedc0vIZO5xzziUtb4k555xLWqmcxJTKJ+cKJqmXmQ2POo6SwF+Lbfy12MZfi+SwV9QBuMj0ijqAEsRfi238tdjGX4sk4EnMOedc0vIk5pxzLml5Ektf3te/jb8W2/hrsY2/FknAB3Y455xLWt4Sc845l7Q8iTnnnEtansScc84lLU9izjnnkpZPO5VGJJ0O3A3sR/B/L8DMrFKkgUVAUiPgBmB/4v4OzKxbVDFFSdLhbP9avB5ZQBGTVIncr8XqCMNxO+CjE9OIpIXA2cA8S/P/eElfACOAeUDOzZbM7OPIgoqIpKeBw4Ev2fZamJldHl1U0ZB0NTAA2AQ5U7+bmR0QXVRuRzyJpRFJHwGnmKXwHfISJGmamR0ddRwlgaSvzOzQqOMoCSQtAI41s5VRx+IS492J6eVWYIykj4HN2YVmNii6kCIzWNJdwHvkfi1mRxdSZKZIOtTMvoo6kBLge2Bj1EG4xHkSSy/3AhuAvYEyEccStcOAHsDJxHWhhcvpZiRBIltGkNCzr5UeHm1YkegHTJY0jdwfbvpEF5LbEU9i6aWamXWMOogS4izgADPbEnUgJcDTBAk91/XBNPUk8CH+WiQNT2Lp5X1JHc3svagDKQG+AKoAK6IOpAT4ycxGRx1ECZFpZjdHHYRLnA/sSCOS1gPlCbpJtpLeQ+wnEozIm0HubqO0G2IvaShBQn+b3K9F2g2xl3Qv8CPbvxY+xL6E8iTm0pKkdvmVp+kQ+2fyKU7XIfb/l0+xD7EvwTyJpQFJnYCKZjYqT/mfgF/NbEI0kRU/SU2A2mY2KU/5icAvZvZ9NJE553aHTzuVHv4J5NfC+JDgi53p5FFgfT7lG8N1aUPSg5J651N+k6QHoogpKpIukdQjn/Krwg97roTyllgakDS3oOHSO1qXiiTNN7PmBaybZ2aHFXdMUZH0FdA875ffJe0FzC3odUpFkj4HTjSz9XnKKwEfmVmraCJzO+MtsfSwt6TtRqJKKg2UiyCeKO29g3Xp9lpYfrO3hGWKIJ4oZeRNYABmtg4oHUE8LkGexNLD68BTkspnF4TPh4Xr0skMSVflLZR0BTArgniitFFS07yFYdmmCOKJUun4v49skiriEwOUaN6dmAbCVtg9wJUEw4cFNCSYAPfvZrY1wvCKlaTawBvAFrYlrdYEb1RnmdmyqGIrbpK6AI8R/G7Evxb9gBvNbExUsRU3SbcApwDXmNkPYdn+wBBgopn9K7Lg3A55EksjksoBTcLFhWaWbp+2c0g6Cci+5vOlmX0YZTxRkdQc6EvcawH8y8zmRRdVNMJBLv2ACmHRBuB+M3siuqjczngSSzOS2rL9faOejSygCEnKAGqT+7X4KbqIXEkgqQLBe2N+o1hdCePTTqURSc8BjYE5QFZYbEDaJTFJNwB3AcvJPQFw2ozUzCbpQOAWtv9wk3aTIUsqC3QD9o8fDGVm6fZVlKThSSy9tAYOTfcbYob+AhxkZquiDqQEeI1gkM9/2PbhJl29BawluEa4eSd1XQngSSy9zAfqAEujDqQE+JngzcoFk976dZ9AAzPrHHUQLnGexNKApLcJusoqAl9Jmk6aTnorKXuG8kXAREnvkqY3CJVULXz6tqRrCUZtpvukt5MlHZaOA1uSlSex9PBQ1AGUIBXDnz+FjzJs+x5QunWzziI45+wvNveNW2dA2kx6K2kewTmXAv4saRF+g9Ck4KMT04ikB8zsbzsrSweSzjOz13ZWlg4k7W1mf+ysLJVJ2m9H683sx+KKxe0an7EjvXTIp6xLsUdRMvRLsCwdTE6wLGWZ2Y9horon+3l8WdTxuYJ5d2IakHQNcC1wgKS5casqApPy3yo1hbNUdAXqS/p33KpKQGY0UUVDUh2gPlBOUku2dStWAvaJLLBoNYtfCL9L6JP/lmCexNLDi8BYYCBwW1z5+jS8eL+E4FpQN3LPlbgeuCmSiKLTCbgMaADED2hZD9weRUBRkdSP4JzLSVqXXUwwPdnwyAJzO+XXxNJI3Gi0eOvTae7EbJJKp+N550fSOWb2v6jjKAkkDTSzdO1WTkqexNKIpB8IJv79jeBTZhWC74ytAK4ys5SfxT1uFFq+0nEUWtzXDuKtBWaZ2ZzijicKko7c0Xozm11csbhd492J6WUc8IaZjQeQ1BHoDLwKDAWOjjC24nJ6+PO68Odz4c+LCe7unI5ah4+3w+XTgBlAb0mvmdmDkUVWfB4Of+5N8Fp8QfBB73BgGnB8RHG5nfCWWBqRNNPMWudXJmmOmbWIKrbiJmmSmR23s7J0IGk8cI6ZbQiXKwCjgLMIWmOHRhlfcZL0MnBv9pedw1n+bzGzyyINzBXIh9inl9WS/iZpv/BxK/BbOAJruzv8prjyknI+XYez+293U8Q0sS/BAIZsW4H9wlv1pNv8gQfHz9ZhZvOBtPlwl4y8OzG9/Ilg5vY3CbpKPgvLMoDzI4wrClcAT0uqHC6vAS6PMJ4ovQhMlfRWuHwG8FJ4p+OvogsrEl9L+g/wPMG100uAr6MNye2Idye6tCapEsHfQVpPBiypNXAc4YcbM5sZcUiRkLQ3cA1wYlj0CfBEOs1ekmw8iaURv28USLrEzJ4vYEReWk0AHM9vEOqSlXcnphe/b9S2614Vd1grjeS5QWgW4aS3pNENQiW9ambnF/QVjHT86kWy8JZYGpE0y8zSegodSVXN7Leo4yhJJC0Ejk7nG4RKqmtmSwuaCNgnAC65vCWWXvy+UfCtpF8JJridBEw2s+8ijilqfoNQuEDSJOBzM0urOTSTnbfE0oik/8un2Mwsbe4bBTnXBtvGPWoCU4FJafLF3lwkjQAOAtL5BqEPEfwuHAzMZduHnClp9iEv6XgSc2lNUmOCWe3/AtQ3s3IRh1TsJN2VX7mZ/bO4Y4mapDIEM3a0BY4NH2vS6Qvfyca7E9OIpH2Am4F9zayXpKbAQWb2TsShFZvwS83Zb1ANgUUErbBLgLScHy87WUkqb2a/Rx1PxMoR3IqmcvhYAszb4RYuUt4SSyOSXiG4/UhPM2suqRxBd0nazEggKUaQrAYBb5pZus6XmEPSscAIoIKZ7SvpCOBqM7s24tCKjaThBPcSW08wV+JUYKoPAir5fNqp9NI4vOazFSCcVkg73iTl1APuA44ExkmaLOlxSRdLSqtrg3EeJbi32CoAM/uCbV/2TRf7AmWBZcAvwGKCWVxcCefdiellS9j6Msi5HpRWc+OZ2TLg9fCR3cV6OfBPoBHBFFxpx8x+lnJ9nkmr7xGaWWcFL0Azgu7mvwLNJa0m6K3I97qhi54nsfRyF8HtWBpKeoFgmqHLIo2omIVzJR7LtpGJLYGFBLchmRRhaFH6ObxWaOHAhj6k4XyBFlxbmS9pDcFXDtYS3LqnDcHfjiuB/JpYmpFUHTiGoBtxKlDGzJZEG1XxCb8jNpVgCPVkYHrYrZq2JNUABgOnEvxevAf0Saeh5ZL6EHyoOY6gu30SMCX8Oc/M0u0uD0nDk1iak/STme0bdRyuZJF0o5k9GnUcxUXSIMLvhpnZ0qjjcYnzJJbmJP1sZg2jjqO4SHqbfObGy2Zm3YoxnBIr3T7cSKq2o/Xp1CpNNn5NzKXbp5iHog4gSaTbqNVZBH8L+Z23Aek6crXE8ySWBiQ9Rv7JSkCVYg4nUmb2cdQxJIm0+nBjZo2ijsHtHk9i6WFHNzhM15sfNgUGAocCe2eXp9M8kpLWU/CHm7SbfiubpKpAU3L/XnwSXURuR/yamEtLkj4jGDb9CHAG8GeCvwcfSp3GJF1JMI9mA2AOwUjeKel049hk40ksDfhghu1l31tN0jwzOyws+9TMTog6tuLigxm2F94U8yiCKadaSDoY+KeZXRBxaK4A3p2YHnwww/b+kLQXsEDS9QRTDdWKOKbi5oMZtveHmf0hCUllzewbSQdFHZQrmCexNOCDGfJ1I7APwewUdwMnA5dGGlEx88EM+VosqQrwJjBB0m8EM9m7Esq7E9OID2bYnqRKBDMOrY86lij5YIbtSWpHcDuWsWa2Nep4XP58Fvv08gzwBJAJnAQ8CzwXaUQRkdQ6vP4xF5gn6QtJraKOKwrhYIZPgPEEEyGPB/pHGVNUJOX8PZjZx2Y2Gng6wpDcTngSSy/lzOwDghb4j2bWn6AbLR09DVxrZvub/X979xIiRxVGcfx/mFEUfIFGjSAoBmMUFOcAAAPTSURBVCHxPWRQXEQ0bnQlGAho3CioEERxN0bdiApuBEVkAoGISHwwBN0I2YhMBIkkBiMKrhRdiEoGE4Jvj4tbsVsz0xld9E3lnt9qqnoWHw1VX9e9597yJcBmSpNv0cOUMMNXtm+mbIr8fd2Sqrli+EDSBNDkj5u+yJxYWxJmGDhse/7oge3d3bqpFjUfZpA0AzwGnC7pEIOwy6/A1mqFxXFlTqwhkqYpr9g4hxJmOBt4zvaHVQurQNLzlGDHDkoSbyOwAMwB2N5Xr7rxkrSTsk7uEcqT+QJwiu3bqxZWgaRnbc/UriOWL02sQQkzgKT3RnzsVhe3th5m6EYq7gIutf2UpIuBlbb3VC4tlpAm1hBJaynzPmd2p34E7rW9t15VUZukV23fc7xzLZD0MvAncIvt1V1qc5ft6cqlxRIS7GhLwgwdSRdI2ibp3e54jaT7atdVScIMA9fb3gz8DGB7ATi1bkkxSppYW44JMwCtDilup0TJL+qOv6DMCTVD0kwXZrla0iFJh7vj74C3K5dXy29dEzeApBWUJ7M4QWU4sSEJMwxI+sj2tKSPbV/Xndtv+9ratY1bwgwDku6mXBdTwCvABuBx229VLSyWlIh9W47eoP+9U/uNlKbWUpjhiKRzGfzivoEyR9iiLZI2kTADtl+TtBdYT4nZ32H788plxQh5EosmSZoCXgSuBD4FVgAbbH9StbAKEmYASacBDwKrgAPANtu/160qliNzYg1JmKGslZN0YTd0ehNlgesvwC7gm6rF1ZMwQxk6XEtpYLeRNz/0RppYW7bTeJgBmKXswgBlGHUL8BJlbrDVnRkSZoA1tjfZnqXMg62rXVAsT5pYW86z/SbdDaobLvmjbkljNzH0sseNwFbbc7afoAwltegFYCdwvqSngd3AM3VLGru/F3ZnGLFfEuxoS8IMMCFpsrtRrQfuH/qsyeshYQYArun2TITyHQzvoWjbZ9UrLUZp8qJt2KPAO8Blkj6gCzPULWnsdgDvS/oB+AmYB5C0isYa+iJhhtlWn0JsT9SuIf6fpBMb0G38+7XtbyVNAg8AdwKfAU8ODa81oXsCXUlJ4B3pzl0OnNHYWrk3KMNo85Qww5e2W5sjjZ5LE2uApH3ArbYPSloHvA48RFk3ttp2a09jAUg6YPuq7u9JYI/tqcplRfwnGU5sw6JhBmBO0v6KdUVd/wgzSBr1vxEnpDSxNiTMEItJmCF6LzewNiTMEMdImCFOBpkTa0TCDBFxMkoTi4iI3sqOHRER0VtpYhER0VtpYhER0VtpYhER0VtpYhER0Vt/AWbGsjJ8Uk1SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import seaborn as sns\n",
    "import scipy.stats\n",
    "\n",
    "df = pd.read_csv('Iris.csv')\n",
    "df = df.drop(columns=['Id', 'Species'])\n",
    "sns.heatmap(df.corr(), annot=True)\n",
    "\n",
    "def pearson_corr(x, y):\n",
    "    x_mean = np.mean(x)\n",
    "    y_mean = np.mean(y)\n",
    "    num = [(i - x_mean)*(j - y_mean) for i,j in zip(x,y)]\n",
    "    den_1 = [(i - x_mean)**2 for i in x]\n",
    "    den_2 = [(j - y_mean)**2 for j in y]\n",
    "    correlation_x_y = np.sum(num)/np.sqrt(np.sum(den_1))/np.sqrt(np.sum(den_2))\n",
    "    return correlation_x_y\n",
    "\n",
    "print(pearson_corr(df['SepalLengthCm'], df['PetalLengthCm']))\n",
    "print(scipy.stats.pearsonr(df['SepalLengthCm'], df['PetalLengthCm']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical Analysis\n",
    "\n",
    "We can approxite the histogram of a dataset with a combination of Gaussian (Normal) distribution functions:\n",
    "\n",
    "    Gaussian Mixture Model (GMM)\n",
    "    \n",
    "    Kernel Density Estimation (KDE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the goal for using GMM?\n",
    "\n",
    "- We want to appriximate the density (histogram) of any given data sample with a combination of Normal Distribitions\n",
    "\n",
    "- How many Normal Distributions we need is defined by us (2 or 3 or 4,...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activity : Fit a GMM to a given data sample \n",
    "\n",
    "Task: \n",
    "\n",
    "1. Generate the concatination of the random variables as follows:\n",
    "\n",
    "`x_1 = np.random.normal(-5, 1, 3000)`\n",
    "\n",
    "`x_2 = np.random.normal(2, 3, 7000) `\n",
    "\n",
    "`x = np.concatenate((x_1, x_2))`\n",
    "\n",
    "2. Plot the histogram of `x`\n",
    "3. Obtain the weights, mean and variances of each Gassuian\n",
    "\n",
    "You will need to use the following in your solution: \n",
    "```python\n",
    "from sklearn import mixture\n",
    "\n",
    "gmm = mixture.GaussianMixture(n_components=2)\n",
    "\n",
    "gmm.fit(x.reshape(-1,1))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEHRJREFUeJzt3XGsnfVdx/H3x1aYorIN7uZswdZQ/yjZQrTrTNQ5h7IikW7aapnJOiXpTOxf02gXlbFuf8Ci1sRVsxpQhGBBEmIjdRVH4pKFYS8bwjpWd9dVuOsixSIGF8Y6vv5xno7j2S33ufceetv7e7+Sm/s8v+f3O+d7nl4+5+E5z/k9qSokSW34rsUuQJJ05hj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYsX+wCRl188cW1atWqxS5Dks4pDz/88NNVNTFbv7Mu9FetWsXk5ORilyFJ55Qk/9Gnn6d3JKkhhr4kNcTQl6SG9Ar9JBuSHE4ylWTHDNvfmuSzSU4m2TTUfkWSB5McSvJokl8dZ/GSpLmZNfSTLAN2A1cDa4Hrkqwd6fYE8F7gzpH2rwPvqarLgQ3AnyZ59UKLliTNT5+rd9YDU1V1BCDJXmAj8IVTHarqaLftxeGBVfXvQ8vHkjwFTAD/veDKJUlz1uf0zgrgyaH16a5tTpKsB84DvjzXsZKk8egT+pmhbU73WEzyBuB24Ner6sUZtm9LMplk8vjx43N5aEnSHPQJ/WngkqH1lcCxvk+Q5AeA+4A/qKrPzNSnqvZU1bqqWjcxMesXyiRJ89TnnP5BYE2S1cBXgS3Au/s8eJLzgHuBv6mqv5t3lXpZq3bcN++xR2+6ZoyVSDrbzXqkX1Unge3AAeBx4O6qOpRkZ5JrAZK8Ock0sBn4eJJD3fBfAd4KvDfJI93PFa/IK5EkzarX3DtVtR/YP9J2w9DyQQanfUbH3QHcscAaJUlj4jdyJakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNaRX6CfZkORwkqkkO2bY/tYkn01yMsmmkW1bk3yp+9k6rsIlSXM3a+gnWQbsBq4G1gLXJVk70u0J4L3AnSNjXwt8EHgLsB74YJLXLLxsSdJ89DnSXw9MVdWRqnoB2AtsHO5QVUer6lHgxZGx7wDur6oTVfUMcD+wYQx1S5LmoU/orwCeHFqf7tr6WMhYSdKY9Qn9zNBWPR+/19gk25JMJpk8fvx4z4eWJM1Vn9CfBi4ZWl8JHOv5+L3GVtWeqlpXVesmJiZ6PrQkaa76hP5BYE2S1UnOA7YA+3o+/gHgqiSv6T7AvaprkyQtgllDv6pOAtsZhPXjwN1VdSjJziTXAiR5c5JpYDPw8SSHurEngA8zeOM4COzs2iRJi2B5n05VtR/YP9J2w9DyQQanbmYaeytw6wJqlCSNid/IlaSGGPqS1BBDX5IaYuhLUkMMfUlqSK+rdyQNrNpx37zHHr3pmjFWIs2PR/qS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEO+cJZ0h3nVLZwNDX81ZSPhK5zpP70hSQ3qFfpINSQ4nmUqyY4bt5ye5q9v+UJJVXft3J7ktyWNJHk/ygfGWL0mai1lDP8kyYDdwNbAWuC7J2pFu1wPPVNVlwC7g5q59M3B+Vb0R+HHgfafeECRJZ16fI/31wFRVHamqF4C9wMaRPhuB27rle4ArkwQo4IIky4HvAV4A/mcslUuS5qxP6K8Anhxan+7aZuxTVSeBZ4GLGLwB/C/wNeAJ4I+q6sQCa5YkzVOf0M8MbdWzz3rgW8APAauB307yI9/xBMm2JJNJJo8fP96jJEnSfPQJ/WngkqH1lcCx0/XpTuVcCJwA3g18oqq+WVVPAZ8G1o0+QVXtqap1VbVuYmJi7q9CktRLn9A/CKxJsjrJecAWYN9In33A1m55E/BAVRWDUzpvz8AFwE8AXxxP6ZKkuZo19Ltz9NuBA8DjwN1VdSjJziTXdt1uAS5KMgW8Hzh1Wedu4PuAzzN48/irqnp0zK9BktRTr2/kVtV+YP9I2w1Dy88zuDxzdNxzM7VLkhaH38iVpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDfEmKmcBb+qh2Sz0b8Q7b+kUj/QlqSGGviQ1xNCXpIYY+pLUED/I1TnHD76l+fNIX5IaYuhLUkM8vdO4hZwq8dpv6dzjkb4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQXqGfZEOSw0mmkuyYYfv5Se7qtj+UZNXQtjcleTDJoSSPJXnV+MqXJM3FrKGfZBmwG7gaWAtcl2TtSLfrgWeq6jJgF3BzN3Y5cAfwm1V1OfA24Jtjq16SNCd9jvTXA1NVdaSqXgD2AhtH+mwEbuuW7wGuTBLgKuDRqvo3gKr6r6r61nhKlyTNVZ/QXwE8ObQ+3bXN2KeqTgLPAhcBPwpUkgNJPpvkd2d6giTbkkwmmTx+/PhcX4Mkqac+oZ8Z2qpnn+XATwG/1v1+V5Irv6Nj1Z6qWldV6yYmJnqUJEmajz6hPw1cMrS+Ejh2uj7defwLgRNd+79U1dNV9XVgP/BjCy1akjQ/fUL/ILAmyeok5wFbgH0jffYBW7vlTcADVVXAAeBNSb63ezP4GeAL4yldkjRXs94usapOJtnOIMCXAbdW1aEkO4HJqtoH3ALcnmSKwRH+lm7sM0n+hMEbRwH7q2r+9+eTJC1Ir3vkVtV+BqdmhttuGFp+Hth8mrF3MLhsU5K0yPxGriQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhvb6cJenctmrH/L8If/Sma8ZYiRabR/qS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIV6yqUWxkEsIJc2fR/qS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGtIr9JNsSHI4yVSSHTNsPz/JXd32h5KsGtl+aZLnkvzOeMqWJM3HrKGfZBmwG7gaWAtcl2TtSLfrgWeq6jJgF3DzyPZdwD8uvFxJ0kL0OdJfD0xV1ZGqegHYC2wc6bMRuK1bvge4MkkAkrwTOAIcGk/JkqT56hP6K4Anh9anu7YZ+1TVSeBZ4KIkFwC/B3zo5Z4gybYkk0kmjx8/3rd2SdIc9Qn9zNBWPft8CNhVVc+93BNU1Z6qWldV6yYmJnqUJEmajz7z6U8DlwytrwSOnabPdJLlwIXACeAtwKYkHwVeDbyY5Pmq+tiCK9eic0586dzTJ/QPAmuSrAa+CmwB3j3SZx+wFXgQ2AQ8UFUF/PSpDkluBJ4z8CVp8cwa+lV1Msl24ACwDLi1qg4l2QlMVtU+4Bbg9iRTDI7wt7ySRUuS5qfX7RKraj+wf6TthqHl54HNszzGjfOoT5I0Rn4jV5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSG9ZtmU1K6F3Czn6E3XjLESjYNH+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkN6hX6SDUkOJ5lKsmOG7ecnuavb/lCSVV37zyd5OMlj3e+3j7d8SdJczBr6SZYBu4GrgbXAdUnWjnS7Hnimqi4DdgE3d+1PA79YVW8EtgK3j6twSdLc9TnSXw9MVdWRqnoB2AtsHOmzEbitW74HuDJJqupzVXWsaz8EvCrJ+eMoXJI0d31CfwXw5ND6dNc2Y5+qOgk8C1w00ueXgc9V1TfmV6okaaH6zKefGdpqLn2SXM7glM9VMz5Bsg3YBnDppZf2KEmSNB99jvSngUuG1lcCx07XJ8ly4ELgRLe+ErgXeE9VfXmmJ6iqPVW1rqrWTUxMzO0VSJJ66xP6B4E1SVYnOQ/YAuwb6bOPwQe1AJuAB6qqkrwauA/4QFV9elxFS5LmZ9bQ787RbwcOAI8Dd1fVoSQ7k1zbdbsFuCjJFPB+4NRlnduBy4A/TPJI9/O6sb8KSVIvve6RW1X7gf0jbTcMLT8PbJ5h3EeAjyywRknSmPiNXElqSK8jfUmaj1U77pv32KM3XTPGSnSKR/qS1BBDX5IaYuhLUkM8pz8mCzl3KUlnikf6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ5x7R9JZybn4Xxke6UtSQwx9SWqIoS9JDTH0Jakhhr4kNcSrdyQtOQu9k91SvvrHI31JaoihL0kN6RX6STYkOZxkKsmOGbafn+SubvtDSVYNbftA1344yTvGV7okaa5mDf0ky4DdwNXAWuC6JGtHul0PPFNVlwG7gJu7sWuBLcDlwAbgz7vHkyQtgj4f5K4HpqrqCECSvcBG4AtDfTYCN3bL9wAfS5KufW9VfQP4SpKp7vEeHE/547XQD38kLQ1LeQqIPqG/AnhyaH0aeMvp+lTVySTPAhd17Z8ZGbti3tVK0lnubH/D6BP6maGtevbpM5Yk24Bt3epzSQ73qOt0LgaeXsD4pcR98RL3xUvcFwNn3X7IzQsa/sN9OvUJ/WngkqH1lcCx0/SZTrIcuBA40XMsVbUH2NOn4NkkmayqdeN4rHOd++Il7ouXuC8GWt0Pfa7eOQisSbI6yXkMPpjdN9JnH7C1W94EPFBV1bVv6a7uWQ2sAf51PKVLkuZq1iP97hz9duAAsAy4taoOJdkJTFbVPuAW4Pbug9oTDN4Y6PrdzeBD35PAb1XVt16h1yJJmkUGB+RLR5Jt3emi5rkvXuK+eIn7YqDV/bDkQl+SdHpOwyBJDVkyoZ9kc5JDSV5Msm5kW7NTQSS5MclXkzzS/fzCYtd0Js02hUhLkhxN8lj3dzC52PWcSUluTfJUks8Ptb02yf1JvtT9fs1i1nimLJnQBz4P/BLwqeFGp4IAYFdVXdH97F/sYs6UnlOItOZnu7+D1i5V/GsG//0P2wF8sqrWAJ/s1pe8JRP6VfV4Vc30pa5vTwVRVV8BTk0FoaXv21OIVNULwKkpRNSYqvoUgysLh20EbuuWbwPeeUaLWiRLJvRfxkzTSLQ2FcT2JI92/4vbxP/Cdvy3//8K+KckD3ffgm/d66vqawDd79ctcj1nxDl156wk/wz84Aybfr+q/v50w2ZoW1KXLL3cfgH+Avgwg9f8YeCPgd84c9UtqiX/bz9HP1lVx5K8Drg/yRe7I2A15JwK/ar6uXkM6zUVxLms735J8pfAP7zC5ZxNlvy//VxU1bHu91NJ7mVw+qvl0P/PJG+oqq8leQPw1GIXdCa0cHqn6akguj/mU97F4APvVvSZQqQJSS5I8v2nloGraOtvYSbD08dsBU53tmBJOaeO9F9OkncBfwZMAPcleaSq3uFUEHw0yRUMTmscBd63uOWcOaebQmSRy1osrwfuHdzmguXAnVX1icUt6cxJ8rfA24CLk0wDHwRuAu5Ocj3wBLB58So8c/xGriQ1pIXTO5KkjqEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JD/g/+932J2mV5zgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.09650276]\n",
      " [-4.97020333]]\n",
      "[[[8.57778373]]\n",
      "\n",
      " [[1.07352143]]]\n",
      "[0.68948643 0.31051357]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import mixture\n",
    "\n",
    "# Generate data samples and plot its histogram\n",
    "x_1 = np.random.normal(-5, 1, 3000)\n",
    "x_2 = np.random.normal(2, 3, 7000) \n",
    "x = np.concatenate((x_1, x_2))\n",
    "plt.hist(x, bins=20, density=1)\n",
    "plt.show()\n",
    "\n",
    "# Define a GMM model and obtain its parameters\n",
    "gmm = mixture.GaussianMixture(n_components=2)\n",
    "gmm.fit(x.reshape(-1,1))\n",
    "print(gmm.means_)\n",
    "print(gmm.covariances_)\n",
    "print(gmm.weights_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The GMM has learned the probability density function of our data sample\n",
    "\n",
    "Let's generate the samples from the GMM:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEAVJREFUeJzt3XGsnfVdx/H3x1aYorINumVrwVbp/ijZQuZdMVFxDseKRLrF4spM1hmSzmT9axrtEmWs7g9YVPxjjVkNKLLMgiTERurqNhKXLBv2sk2gYN1drXBXIsUiBhfGOr7+cZ+Gs7PT3efee+i55fd+JTf3eX7P9znne05uPud3n3Oe56SqkCS14Ucm3YAk6cwx9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNWTnpBoZdeOGFtXbt2km3IUlnlQcffPDpqlo1X92yC/21a9cyPT096TYk6ayS5D/71Hl4R5Ia0iv0k2xKcjjJTJKdI7ZfkeSrSU4m2TIwflmSLyc5lOShJO8dZ/OSpIWZN/STrAB2A1cDG4Drk2wYKnsc+ADwmaHxbwPvr6pLgU3Anyd59VKbliQtTp9j+huBmao6ApBkL7AZePRUQVUd7ba9OLhjVf37wPKxJE8Bq4D/WXLnkqQF63N4ZzXwxMD6bDe2IEk2AucA31zovpKk8egT+hkxtqBvXknyBuBO4Ler6sUR27cnmU4yffz48YXctCRpAfqE/ixw0cD6GuBY3ztI8lPAfcAfVtVXRtVU1Z6qmqqqqVWr5v2YqSRpkfqE/kFgfZJ1Sc4BtgL7+tx4V38v8DdV9XeLb1OSNA7zhn5VnQR2AAeAx4C7q+pQkl1JrgVI8rYks8B1wKeSHOp2/03gCuADSb7e/Vz2sjwSSdK8sty+GH1qaqpaOyN37c77lrT/0ZuvGVMnks5WSR6sqqn56jwjV5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkN6hX6STUkOJ5lJsnPE9iuSfDXJySRbhrZtS/KN7mfbuBqXJC3cvKGfZAWwG7ga2ABcn2TDUNnjwAeAzwzt+1rgo8DlwEbgo0les/S2JUmL0WemvxGYqaojVfUCsBfYPFhQVUer6iHgxaF93wV8rqpOVNUzwOeATWPoW5K0CH1CfzXwxMD6bDfWx1L2lSSNWZ/Qz4ix6nn7vfZNsj3JdJLp48eP97xpSdJC9Qn9WeCigfU1wLGet99r36raU1VTVTW1atWqnjctSVqoPqF/EFifZF2Sc4CtwL6et38AuCrJa7o3cK/qxiRJEzBv6FfVSWAHc2H9GHB3VR1KsivJtQBJ3pZkFrgO+FSSQ92+J4A/Zu6F4yCwqxuTJE3Ayj5FVbUf2D80duPA8kHmDt2M2vd24PYl9ChJGhPPyJWkhhj6ktQQQ1+SGmLoS1JDer2RK2nO2p33LXrfozdfM8ZOpMVxpi9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ3xjFw1Zyln1UpnO2f6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIZ4ctaYeMLPmeNzLS2eM31JaoihL0kN6RX6STYlOZxkJsnOEdvPTXJXt/2BJGu78R9NckeSh5M8luQj421fkrQQ84Z+khXAbuBqYANwfZINQ2U3AM9U1SXArcAt3fh1wLlV9Wbg54APnnpBkCSdeX1m+huBmao6UlUvAHuBzUM1m4E7uuV7gCuTBCjgvCQrgR8DXgD+dyydS5IWrE/orwaeGFif7cZG1lTVSeBZ4ALmXgD+D3gSeBz4k6o6scSeJUmL1Cf0M2KsetZsBL4HvBFYB/xukp/5gTtItieZTjJ9/PjxHi1JkhajT+jPAhcNrK8Bjp2upjuUcz5wAngf8Nmq+m5VPQV8CZgavoOq2lNVU1U1tWrVqoU/CklSL31C/yCwPsm6JOcAW4F9QzX7gG3d8hbg/qoq5g7pvCNzzgN+Hvi38bQuSVqoeUO/O0a/AzgAPAbcXVWHkuxKcm1XdhtwQZIZ4MPAqY917gZ+AniEuRePv6qqh8b8GCRJPfW6DENV7Qf2D43dOLD8PHMfzxze77lR45KkyfCMXElqiKEvSQ0x9CWpIV5a+RVgKZcaPnrzNWPsRNJy50xfkhpi6EtSQzy8I50hHobTcuBMX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEL9ERROxlC8UkbR4zvQlqSGGviQ1xNCXpIb0Cv0km5IcTjKTZOeI7ecmuavb/kCStQPb3pLky0kOJXk4yavG174kaSHmDf0kK4DdwNXABuD6JBuGym4AnqmqS4BbgVu6fVcCnwZ+p6ouBd4OfHds3UuSFqTPp3c2AjNVdQQgyV5gM/DoQM1m4KZu+R7gk0kCXAU8VFX/ClBV/z2mvqWmLPXTTkdvvmZMnehs1+fwzmrgiYH12W5sZE1VnQSeBS4A3gRUkgNJvprk90fdQZLtSaaTTB8/fnyhj0GS1FOf0M+IsepZsxL4ReC3ut/vSXLlDxRW7amqqaqaWrVqVY+WJEmL0Sf0Z4GLBtbXAMdOV9Mdxz8fONGN/3NVPV1V3wb2A29datOSpMXpE/oHgfVJ1iU5B9gK7Buq2Qds65a3APdXVQEHgLck+fHuxeCX+f73AiRJZ9C8b+RW1ckkO5gL8BXA7VV1KMkuYLqq9gG3AXcmmWFuhr+12/eZJH/G3AtHAfuryvPvJWlCel17p6r2M3doZnDsxoHl54HrTrPvp5n72KYkacI8I1eSGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpISsn3YDOXmt33jfpFiQtkDN9SWqIoS9JDTH0Jakhhr4kNaRX6CfZlORwkpkkO0dsPzfJXd32B5KsHdp+cZLnkvzeeNqWJC3GvKGfZAWwG7ga2ABcn2TDUNkNwDNVdQlwK3DL0PZbgX9ceruSpKXoM9PfCMxU1ZGqegHYC2weqtkM3NEt3wNcmSQASd4NHAEOjadlSdJi9Qn91cATA+uz3djImqo6CTwLXJDkPOAPgI/9sDtIsj3JdJLp48eP9+1dkrRAfU7Oyoix6lnzMeDWqnqum/iPVFV7gD0AU1NTw7ctaYmWciLd0ZuvGWMnmrQ+oT8LXDSwvgY4dpqa2SQrgfOBE8DlwJYknwBeDbyY5Pmq+uSSO5ckLVif0D8IrE+yDvgWsBV431DNPmAb8GVgC3B/VRXwS6cKktwEPGfgS9LkzBv6VXUyyQ7gALACuL2qDiXZBUxX1T7gNuDOJDPMzfC3vpxNS5IWp9cF16pqP7B/aOzGgeXngevmuY2bFtGfJGmMPCNXkhpi6EtSQwx9SWqIX6LSOL8IRWqLM31JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXEL0aX9EOt3Xnfovc9evM1Y+xE4+BMX5Ia0iv0k2xKcjjJTJKdI7afm+SubvsDSdZ24+9M8mCSh7vf7xhv+5KkhZg39JOsAHYDVwMbgOuTbBgquwF4pqouAW4FbunGnwZ+vareDGwD7hxX45Kkhesz098IzFTVkap6AdgLbB6q2Qzc0S3fA1yZJFX1tao61o0fAl6V5NxxNC5JWrg+b+SuBp4YWJ8FLj9dTVWdTPIscAFzM/1TfgP4WlV9Z/HtvryW8oaVJJ0N+oR+RozVQmqSXMrcIZ+rRt5Bsh3YDnDxxRf3aEmStBh9Du/MAhcNrK8Bjp2uJslK4HzgRLe+BrgXeH9VfXPUHVTVnqqaqqqpVatWLewRSJJ66xP6B4H1SdYlOQfYCuwbqtnH3Bu1AFuA+6uqkrwauA/4SFV9aVxNS5IWZ97Qr6qTwA7gAPAYcHdVHUqyK8m1XdltwAVJZoAPA6c+1rkDuAT4oyRf735eN/ZHIUnqpdcZuVW1H9g/NHbjwPLzwHUj9vs48PEl9ihJGhPPyJWkhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIb4HbmSXjZ+v+7y40xfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ3xgmuSliUv1vbycKYvSQ0x9CWpIa+4wztL+ZdQkl7pes30k2xKcjjJTJKdI7afm+SubvsDSdYObPtIN344ybvG17okaaHmneknWQHsBt4JzAIHk+yrqkcHym4AnqmqS5JsBW4B3ptkA7AVuBR4I/D5JG+qqu+N+4FI0ilL/Y//lfxGcJ+Z/kZgpqqOVNULwF5g81DNZuCObvke4Mok6cb3VtV3quo/gJnu9iRJE9DnmP5q4ImB9Vng8tPVVNXJJM8CF3TjXxnad/Wiu5WkM2BS7w2eif8w+oR+RoxVz5o++5JkO7C9W30uyeEefb3SXQg8Pekmlimfm9F8Xk7vrHhucsuSdv/pPkV9Qn8WuGhgfQ1w7DQ1s0lWAucDJ3ruS1XtAfb0abgVSaaramrSfSxHPjej+bycns/NS/oc0z8IrE+yLsk5zL0xu2+oZh+wrVveAtxfVdWNb+0+3bMOWA/8y3halyQt1Lwz/e4Y/Q7gALACuL2qDiXZBUxX1T7gNuDOJDPMzfC3dvseSnI38ChwEviQn9yRpMnJ3IRcy02S7d1hLw3xuRnN5+X0fG5eYuhLUkO89o4kNcTQX8aS3JTkW0m+3v382qR7mqT5LgfSsiRHkzzc/Z1MT7qfSUpye5KnkjwyMPbaJJ9L8o3u92sm2eMkGfrL361VdVn3s3/SzUzKwOVArgY2ANd3l/nQS36l+ztp/aOJfw1sGhrbCXyhqtYDX+jWm2To62zR53IgElX1ReY+RTho8FIxdwDvPqNNLSOG/vK3I8lD3b+szf5LyujLgXhJj5cU8E9JHuzOcNf3e31VPQnQ/X7dhPuZGEN/wpJ8PskjI342A38B/CxwGfAk8KcTbXayel3So2G/UFVvZe7w14eSXDHphrQ8veK+ROVsU1W/2qcuyV8C//Ayt7Oc9bqkR6uq6lj3+6kk9zJ3OOyLk+1qWfmvJG+oqieTvAF4atINTYoz/WWs++M85T3AI6erbUCfy4E0Kcl5SX7y1DJwFW3/rYwyeKmYbcDfT7CXiXKmv7x9IsllzB3GOAp8cLLtTM7pLgcy4baWi9cD9859hQUrgc9U1Wcn29LkJPlb4O3AhUlmgY8CNwN3J7kBeBy4bnIdTpZn5EpSQzy8I0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWrI/wObYDWdoZRHNAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = gmm.sample(10000)\n",
    "plt.hist(z[0], bins=20, density=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question: Are the samples in z and x the same?\n",
    "\n",
    "- Not their samples, but we know that the _pdf of x and the pdf of z are the same_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel Density Estimation (KDE)\n",
    "\n",
    "**Kernel density estimation (KDE)** is a non-parametric way to estimate the Probability Density Function (PDF) of a data sample. In other words _the goal of KDE is to find PDF for a given data sample._\n",
    "\n",
    "We can use the below formula to approximate the PDF of the dataset:\n",
    "\n",
    "$p(x) = \\frac{1}{Nh}\\sum_{i = 1}^{N} \\ K(\\frac{x - x_i}{h})$\n",
    "\n",
    "where $h$ is a bandwidth and $N$ is the number of data points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activity: Apply KDE on a given data sample\n",
    "\n",
    "**Task**: Apply KDE on the previously generated sample data `x`\n",
    "\n",
    "**Hint**: use \n",
    "\n",
    "`kde = KernelDensity(kernel='gaussian', bandwidth=0.6)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a1cee9668>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4nNWZ9/HvPaPeq1Wt5ibLFVvINpgSCC0BDIRmOiEx2Q2bZMmW8O6GULJ5k90Qv0mWbAKEQGjGIQQMOPQsLbjIFRfJFrKsLktW72XO+4dGRhGyPbZHep6ZuT/X5cujmWdmbo2k35w55zzniDEGpZRSgcFhdQFKKaUmj4a+UkoFEA19pZQKIBr6SikVQDT0lVIqgGjoK6VUANHQV0qpAKKhr5RSAURDXymlAkiQ1QWMlZSUZHJycqwuQymlfMqWLVuajDHJxzvOdqGfk5NDcXGx1WUopZRPEZGDnhyn3TtKKRVANPSVUiqAaOgrpVQA0dBXSqkAoqGvlFIBRENfKaUCiIa+UkoFEA19P/b2ngZK6tutLkMpZSMa+n5qd20bq54q5t//tMvqUpRSNqKh74dcLsP3X9qFy0DxwRYOHu6yuiSllE1o6PuhtcVVbK1s5Z8unIkI/HFrjdUlKaVsQkPfzzR39fPj10soykngm1+YzpnTknhxazUul7G6NKWUDWjo+5kf/3kvnb2D/PDKuYgIX1mcQXVLD5srmq0uTSllAxr6fqS4opm1xdXccVYuM1OiAbhoTiqRIU7+uLXa4uqUUnagoe8nBoZc/PtLu0iPDeNb5804cn1ESBCXzEtj/Sf19PQPWVihUsoONPT9xJN/raCkvoMfXD6HyNC/3SbhK4sy6ewb5M099RZVp5SyCw19P1DX1sPqt/Zxfv4ULixI+dztS3ITyIgL11k8SikNfX/w4Kt7GDKG+y6fg4h87naHQ7hqUQYf7m+kob3XggqVUnahoe/jalp7WP9JPavOnsbUhIijHnflaRm4DPxpm7b2lQpkGvo+rqRueG2dc2YmHfO4vOQoFmXF8cct1Rijc/aVClQa+j6upL4DgBnuKZrH8pXFmew/1MmuGl2ETalA5VHoi8jFIlIqImUi8r1xbj9bRLaKyKCIXD3mtltFZL/7363eKlwN29fQQUZcODFhwcc99tJ56YQEOXTOvlIB7LihLyJO4GHgEqAAWCkiBWMOqwRuA54dc98E4AfAEqAI+IGIxJ962WpEaX0HM1OiPDo2NiKYC2ansG5HLf2DrgmuTCllR5609IuAMmNMuTGmH1gDrBh9gDGmwhizExibJBcBbxljmo0xLcBbwMVeqFsxfELWp42dzEqN8fg+Vy3KoLmrn/f2NU5gZUopu/Ik9DOAqlFfV7uv88Sp3Fcdx4GmLgaGDLNSPWvpA5w9M5mkqBD+uEW7eJQKRJ6E/ucnfoOn0z88uq+IrBKRYhEpbmzUFqinSt2DuLNSPG/pBzsdnJ+fwsflh3UWj1IByJPQrwamjvo6E6j18PE9uq8x5hFjTKExpjA5OdnDh1al9R04HcK0KZEndL+5GTG09QxQ16YnaikVaDwJ/c3ADBHJFZEQ4HpgnYeP/wZwoYjEuwdwL3Rfp7ygtKGD3KRIQoOcJ3S/gvThTwZ7anXqplKB5rihb4wZBO5iOKz3AmuNMbtF5AERuRxARE4XkWrgGuA3IrLbfd9m4EGG3zg2Aw+4r1NeUFrfwSwP5uePNSs1BhHYU6ehr1SgCTr+IWCMWQ+sH3PdvaMub2a462a8+z4OPH4KNapxdPcPUtnczdWLx33ZjykqNIicxEht6SsVgPSMXB+1r6ET4MhmKSeqIC1GW/pKBSANfR+1zz1zJz/1JEM/PYbK5m7aewe8WZZSyuY09H1USX0HYcGOY66seSyz04bfLErqOrxZllLK5jT0fdS+hg5mpkTjdIx3KsTxFaTFArCnts2bZSmlbE5D30eV1HecdH8+QEpMKAmRIezVlr5SAUVD3wcd7uyjqbPvpPvzAUREB3OVCkAa+j6otGG4dX4qLX0YHswtbehgYEhX3FQqUGjo+6BTnbkzoiAthv5BF+WNXd4oSynlAzT0fVBpQwdxEcEkR4ee0uMcWY6hTgdzlQoUGvo+aGT5BZGTm7kzIi8pkpAgh56Z6yVtPQN8d+0Orv31xzz+4QEa2nVBO2U/Hi3DoOzDGMO+hk6uWnTq2xIEOR3kp0brYK4XbKts4R+e20Z9Wy95yZE88OoeHnxtD0U5CVy2IJ1L5qaSGHVqn8yU8gYNfR9T09pDZ98gs06xP3/E7NQY3txTjzHmlD85BCKXy/DYh+X85+ulpMSE8YdvLOO0rHjKDnXy6s5aXtlRy7+/tIsfrNvNBbNT+Nl1C4gI0T87ZR3t3vExn22c4p3QL0iPoaV7gHrtijhhzV393PHkZn60voQvzk5h/bfO4rSs4S2gp0+J4jtfnMnbd5/Dn799Fl9bnssbe+r51nPbGXLp5jXKOtrk8DFHpmt6qaU/Mpi7t66dtNhwrzxmINhYfphvrdlGS9cAD66Yw01Ls8f9pCQizE6LYXZaDGmxYdz3yh4eeGU3910+Rz9ZKUto6PuY0voO0mPDiAkL9srjjUz73FPbznn5KV55TH+3s7qVmx/fRGZcOI/fdjpz0mM9ut9tZ+ZS3dLDYx8eYGpCBF87K2+CK1Xq8zT0fUxpfYfX+vMBosOCyU6M0MFcDx3q6GXV77eQHBXKH76x7IQHZ//Pl2ZT09rDD1/bS3pcOF+alzZBlSo1Pu3T9yEDQy4+bez0WtfOiIK0GJ226YG+wSH+7umttPUM8OgthSc1G8fhEFZft5BFWXF85/ntbDmoG8mpyaWh70MqmroYGDKnfCbuWAVpMVQc7qazb9Crj+tPjDH84OXdbDnYwk+vWXBkLORkhAU7eezW00mPDeNrTxZzoEnPiFaTR0Pfh5TUe2fNnbFGAqxEu3iO6ukNB1mzuYpvfmEaX55/6l0yCZEhPHF7ESLCbb/bRGt3vxeqVOr4NPR9yL6GDpwOYVpylFcf97PlGDT0x/Pxp4e5/5U9nJ8/he9eMMtrj5uTFMmjtxRS29rDPS9+gjE6lVNNPA19H1JS30FOYgRhwU6vPm5qTBhxEcHarz+OquZuvvnsVrITI1h9/UIcJ7lpzdEszo7n7gtm8edd9fxhS7VXH1up8Wjo+5B9DR3kp558X/LRjKytv1db+n+jq2+QVU9tYWDIxaO3FHptmuxYq87OY2leAvev283Bw9q/ryaWhr6P6O4fpLK52+v9+SMK0mIoqe9gUNfWB2BwyMU/PLeN0vp2frHyNPK83KU2mtMh/OzahTgdwrfXbNf9DdSE0tD3EfsbOjEGr87RH60gPYa+QZfOJGF4ps59r+zm3ZJD3L9iLl+YNWXCnzM9Lpz/uHIe26ta+eW7ZRP+fCpwaej7iE8bOwGYPiVyQh5fB3M/8+gH5Ty9oZI7z87j5qXZk/a8ly1I56pFGfz3u/t1/r6aMBr6PqKquQeAzPiICXn8aclRhDh1bf3Xdtbxo/UlfHleGv96cf6kP//9l88hIz6c7zy/nY7egUl/fuX/NPR9RFVLNykxoV6fuTMi2OlgZmpUQLf0txxs5h/XbmdxdjwPXbvA6zN1PBEdFszqaxdS09LDD9btnvTnV/5PQ99HVDV3M3WCWvkjRpZjCMT54geauvjak8Wkx4bx6C2FE/bm6onCnATuOm8GL26t4fVddZbVofyTR6EvIheLSKmIlInI98a5PVREnnffvlFEctzXB4vIkyLyiYjsFZF7vFt+4Khq7mZqwsSG/uy0GA539XOoo29Cn8duWrr6uf13mwB44vYiEiJDLK4IvnXedPJTo3nw1b30DgxZXY7yI8cNfRFxAg8DlwAFwEoRKRhz2B1AizFmOrAa+In7+muAUGPMPGAxcOfIG4LyXP+gi7r23gkP/ZFzAEaWewgELpfh289vp7a1l8duLSQnaWIGyk9UkNPBvZcVUNPaw2MflFtdjvIjnrT0i4AyY0y5MaYfWAOsGHPMCuBJ9+UXgPNleIcIA0SKSBAQDvQDgdtpfJJqW3swBqbGT+wmJyMLuZXWB86P6OG/lPH+vkbuvayAxdkJVpfzN86YlsRFc1L41f9+qpusK6/xJPQzgKpRX1e7rxv3GGPMINAGJDL8BtAF1AGVwE+NMZ+biyYiq0SkWESKGxsbT/ib8HdVLd0AE97Sj48MITUmjJK6wGjp/7WsidVv72PFwnRuXJJldTnj+rcvFTA4ZPjJ6yVWl6L8hCehP94UhrEjfUc7pggYAtKBXOC7IvK57YKMMY8YYwqNMYXJyckelBRYKpsnJ/Rh+OSvvQHQvdPQ3su31mwjLzmKH105z7ZbF2YlRvDV5bm8uLWG7VWtVpej/IAnoV8NTB31dSZQe7Rj3F05sUAzcAPwujFmwBhzCPgIKDzVogNNVXMPwU4hNSZswp8rPy2askMdfr0UwMgSC119Q/zPjYuIDLX3BnJ3nTedpKhQHnhld0DOrFLe5UnobwZmiEiuiIQA1wPrxhyzDrjVfflq4F0z/NtZCZwnwyKBpYB+Tj1BVS3dZMSF45yEeeOzU2MYGDJ+vRzDT9/cx6YDzfzfq+YxY4LWMvKmqNAg/uWiWWytbGXdjrHtLaVOzHFD391HfxfwBrAXWGuM2S0iD4jI5e7DfgskikgZcDcwMq3zYSAK2MXwm8fvjDE7vfw9+L3qSZiuOSI/bTgE/XXFzXf2NvDr9z7lhiVZXHHa2KEp+7p6cSZzM2L48Z9L6O7XHc7UyfNonr4xZr0xZqYxZpox5j/c191rjFnnvtxrjLnGGDPdGFNkjCl3X9/pvn6OMabAGPNfE/et+K/KSQz9vKQoghzil9M2q5q7uXvtDuZmxHDvpWNnHdubwyHce+kc6tp6+c17OoVTnTw9I9fmOvsGaekemPCzcUeEBDmYPiXKL7dO/PeXduFyGX51w2JLz7g9WUW5CXx5fhq/ef9Talt7rC5H+SgNfZurOjJzZ2Ln6I+WnxpNqZ+19N/b18h7+xr59hdnkJU4OW+gE+GeS/IxBla/tc/qUpSP0tC3uSOhP0ktfYD8tBhq23pp6/aPVR6HXIYfvbaXrIQIbl42eUslT4TM+AiuP30qL2+v5VCHnrClTpyGvs2NzNHPmqQ+ffhso5YSPzkzd21xFaUNHdxzST6hQb7XrTPWbWfmMuBy8fSGSqtLUT5IQ9/mqlt6iAoNIi5iYvZnHc9sP1qDp7NvkIfeLOX0nHgunptqdTlekZsUyfn5U3hmw0FdjE2dMA19m6tq7iYzPnxSzxhNiQklLiLYL1r6v/7fT2nq7Offvlxg27NuT8ZXl+dyuKufddt13r46MRr6NlfVMnnTNUeICPmp0T7f0q9t7eHRD8q5fEE6C6fGWV2OVy3LSyQ/NZrHPzqgZ+mqE6Khb2PGGKqaeya1P39EfmoMpfUduFy+Gyj/9UYpBviXi2dZXYrXiQh3LM+lpL6Dv3562OpylA/R0Lexps5+egaGJnxJ5fHkp0bT3T90ZIVPX7OzupU/bavhjuW5E7avsNUuW5BOUlQIv/3wgNWlKB+ioW9jk7Wk8njy04YHc/f64DLLxhh++OpeEiND+Ptzp1ldzoQJC3Zy09Js3i05xKeNnVaXo3yEhr6NVU3ikspjzUyJQgSfPEnrjd0NbKpo5h8vmEl02OTNerLCTUuzCXE6eOKjCqtLUT5CQ9/GrDgxa0RESBA5iZE+N4PH5TI89GYp06dEcf3pU49/Bx+XFBXKioXpvLClmtbufqvLUT5AQ9/Gqpp7SIoKJTzEmhOKZqX43gye9/Y1sv9QJ9/8wjSCnIHx6/3V5bn0DAyxZnPV8Q9WAS8w/ip81PB0zckfxB2RnxZNxeEun1rK95H3y0mNCePS+elWlzJpZqfFcMa0RJ78a4Vfb36jvEND38aqWrot6doZkZ8agzGwv8E3Bgl31bTxcflhbjszh+AAaeWPuGN5LnVtvby+q97qUpTNBdZfhg8ZHHJR29pryRz9EbPTfGsNnsc+KCcyxMnKIntucj6RvjBrCrlJkTz+kU7fVMemoW9TdW29DLmMpd07U+MjiAhx+sS0zbq2Hl7dWcd1p2cRG+7fM3bG43AItyzLZltlK7tq2qwuR9mYhr5NWTlzZ4TDIcxMifaJlv4TH1XgMobbz8yxuhTLXHVaJmHBDp7ZqKtvqqPT0LcpK0/MGm122vAMHjuv79LRO8CzGyu5ZF6a5a+XlWIjgrlsfjovb6+ho9c/9kJQ3qehb1OVzd04HUJabJildeSnxtDaPcChjj5L6ziW5zdX0dE3yKqz8qwuxXI3Lc2mu3+Il3T1TXUUGvo2VdXcQ3pcmOVzzfPdG6rstemeuYNDLn73UQVFOQks8LOVNE/G/MxY5mbE8MyGg7b+dKaso6FvU1ZP1xyRb/MNVf68q56a1h6+dlau1aXYgohw45JsSuo72FrZYnU5yoY09G3KqiWVx4qNCCYtNowSG7b0jTE89kE5uUmRfHF2itXl2MblC9KJDg3iGd1OUY1DQ9+GuvsHaerss82gpF03VNlc0cKO6jbuWJ6Lw+E/u2KdqsjQIK5clMGrn9TR0qXr8ai/paFvQ9UtPQBkWrCO/njy02L4tLGT/kF7neL/6AflxEcE85VFmVaXYjs3LMmif9DFC1uqrS5F2YyGvg1ZuaTyePJToxkYMrZas73ycDdv723gpqXZli1IZ2f5qTEUZsfz7KZKn979THmfhr4NjYS+Hfr0AeZlxAKwrbLV4ko+88zGgzjcg5ZqfDctzeZAUxcfl+t2iuozHoW+iFwsIqUiUiYi3xvn9lARed59+0YRyRl123wR+VhEdovIJyJi7cRzH1DZ3EN4sJPEyBCrSwEgNymSpKhQNh6wR3j0DgzxfHEVFxakkGrxeQx2dvHcVOIjgnl6w0GrS1E2ctzQFxEn8DBwCVAArBSRgjGH3QG0GGOmA6uBn7jvGwQ8DXzDGDMHOBfQUwWPY2RJZRF7DE6KCEvyEthY3myLud+v7qyjtXuAm5dpK/9YwoKdXFM4lTf3NNDQ3mt1OcomPGnpFwFlxphyY0w/sAZYMeaYFcCT7ssvAOfLcGJdCOw0xuwAMMYcNsYMead0/1XVbI85+qMtzU2gvr2XquYeq0vhqY8rmD4limV5iVaXYns3FGUx5DKs1Q1WlJsnoZ8BjP6NqXZfN+4xxphBoA1IBGYCRkTeEJGtIvIv4z2BiKwSkWIRKW5sbDzR78GvGGOobumxzSDuiKLc4YC1uotnR1UrO6rbuHlptm0+CdlZTlIkZ81I4rlNlQzqBisKz0J/vL+ssZ/xj3ZMELAcuNH9/5Uicv7nDjTmEWNMoTGmMDk52YOS/FdL9wCdfYO2C/0ZU6KIjwhm44FmS+t4asNBIkKcXLVobLtDHc2NS7KobevlL6WB3aBSwzwJ/Wpg9A7TmcDY1ZyOHOPux48Fmt3Xv2eMaTLGdAPrgUWnWrQ/+2xJZXvM0R/hcAhFuQmWtvRbuvp5ZUctV56WQXRY4K2Zf7K+ODuFlJhQHdBVgGehvxmYISK5IhICXA+sG3PMOuBW9+WrgXfN8IjfG8B8EYlwvxmcA+zxTun+aWRJ5axEe7X0YbiLp6q5h9pWa/r1/7Clir5Blw7gnqAgp4OVRVm8t6+Rg4e7rC5HWey4oe/uo7+L4QDfC6w1xuwWkQdE5HL3Yb8FEkWkDLgb+J77vi3Azxh+49gObDXGvOb9b8N/jAyU2m0gF2BJbgIAmyzo4nG5DE9vqKQoJ+HIInDKcyuLsnA6RDdYUQR5cpAxZj3DXTOjr7t31OVe4Jqj3PdphqdtKg8cPNxFYmQIkaEe/Wgm1ey0GKLDgth4oJkrTpvcPvX39jdS2dzNP180a1Kf11+kxIRx0ZwU1hZXcfcFMwkL1rOYA5WekWszFYe7yEmKtLqMcTkdwuk51vTrP/XxQZKiQrloTuqkP7e/uGlJNq3dA7y2s87qUpSFNPRtpqKpm5xEe4Y+DHfxlDd2cahj8k72qWru5i+lh7ihaCohQfore7KWTUskLzmSp3RAN6DpX5CN9PQPUd/eS26S/frzRyxxnxC1+cDkbdDxtHudnZVLsibtOf2RiHDz0my2V7Wyq6bN6nKURTT0baTCPbMi28Yt/TnpMUSEOCeti6d3YIi1m6u4YHYKabH2msbqi65alEl4sFOnbwYwDX0bGZlOl2vTPn2AYKeDxdnxbCyfnBk8L22roaV7gFt0mqZXxIYHs2JhOi9vr6WtR5fBCkQa+jZyoGl4jr5dB3JHLMlNoLShY8J3ZXK5DI98UM7cjBiWTdN1drzlpqXZ9AwM8eJW3WAlEGno20hFUxdJUaFE2XC65mgj/fqbKia2tf/23gbKG7tYdfY0XWfHi+ZmxLJwahxPbzhoi1VT1eTS0LeRisNd5NjwTNyx5mfGEhrkmPCTtH7zfjmZ8eF8aa5O0/S2m5dm82mjbrASiDT0bcTOc/RHCw1yclpW3IQO5hZXNLPlYAtfW55LkFN/Tb3ty/PTiNMNVgKS/jXZRHf/IA3tfbYexB1tSW4ie2rbae+dmMHA37xfTlxEMNeePvX4B6sTFhbs5NrCqbyxWzdYCTQa+jZx8PDwIG62D3TvACzJS8BlYEuF9+frlx3q5K09DdyyLIeIEHuPb/gy3WAlMGno20RF0/B0TTufjTvaaVPjCXYKGyagi+exD8oJDXJwq07TnFA5SZGcMS2R54urcLl0QDdQaOjbxAH3HH1f6NMHCA9xsiAzzuuDuYfae3lxaw3XFGaSGBXq1cdWn7eyKIvqlh4+LGuyuhQ1STT0bcJXpmuOtiQvgU+q2+juH/TaYz7x1woGXC6+tjzPa4+pju7COSnERwTz3CZdcjlQaOjbRMXhbluvuTOeotxEBl2GLQe906/f2TfIUxsOcsncVJ/5xOPrQoOcXL04k7f2NNDY0Wd1OWoSaOjbREVTl8/0548ozI4nJMjBy9vH7p55ctZsqqSjd5A7z57mlcdTnrnu9CwGXYYXtugZuoFAQ98GuvoGOdTR53Ot28jQIG5ems2LW6spb+w8pccaGHLx+IcHWJKbwIKpcV6qUHli+pQoinITeH5zpQ7oBgANfRsYma7pay19gL87dxphwU5Wv73/lB5nbXEVtW29fOMcbeVbYWXRVCoOd7NBz9D1exr6NlBxZOaOb/XpAyRFhXL7mTm8sqOWvXXtJ/UYVc3d/Oi1vSzLS+TcWclerlB54pK5acSGB/Ocztn3exr6NnDAx+boj7XqrGlEhwXxs7f2nfB9XS7Dv7ywExHhv66ZrwurWSQs2MmVp2Xwxq56mid49VRlLQ19Gzh4uIvk6FBbbobuidiIYFadlcdbexrYXtV6Qvf9/ccVfFx+mHsvLSAz3vc+6fiTlUVZ9A+5dMllP6ehbwMVTd3k+mgrf8Tty3NJiAzhoTdLPb5PeWMnP369hPPyp3BNYeYEVqc8MSs1mkVZcTy7qVKXXPZjGvo2cOBwl0/2548WFRrE350zjQ/2N7HRg8HAwSEX3/3DDkKDnPz4qnnarWMT1xdlUd7YxeYJWFNJ2YOGvsU6+wZp7Oiz9b64nrp5WTZTokP56Zulx20pPvJBOdsqW3nwirlMiQmbpArV8Vw6P43o0CDW6Bm6fktD32K+sC+up8KCnfzDedPZXNHC+/uPvpZLSX07q9/ax5fnpXHZ/LRJrFAdT0RIECtOS+e1T+po69Y9dP2Rhr7FKpp8d47+eK47PYvM+HAeOkprv3/Qxd3P7yA2PJgHr5ir3To2tLIoi75BFy9u0wFdf+Sb00X8yMgcfV9ZR/94QoIcfPv8GfzzCzv50fq9RIQE0dLdz+Guflq6+qlr6+VAUxeP3lJIQmSI1eWqccxJj2WBew/d287I0TdmP+NRS19ELhaRUhEpE5HvjXN7qIg87759o4jkjLk9S0Q6ReSfvFO2/6ho6mKKD0/XHM+Vp2UwMyWKRz84wM/f2c+6HbXsrW2nf9DF9ClR3HdZARcUpFhdpjoG3UPXfx03aUTECTwMXABUA5tFZJ0xZs+ow+4AWowx00XkeuAnwHWjbl8N/Nl7ZfsPX9kX90QEOR28/M3ldPUPEhcerHvc+qBL56fxw9f28PSGg5wxLcnqcpQXefLXWASUGWPKjTH9wBpgxZhjVgBPui+/AJwv7s+EInIFUA7s9k7J/uWAH8zRH094iJOkqFANfB8VFuzkmsWZvKl76PodT/4iM4DRC3JUu68b9xhjzCDQBiSKSCTwr8D9p16q/+nsG6Sps49sH5+jr/zTjUuyGXQZ1mzS9Xj8iSehP94ozthpGUc75n5gtTHmmOvuisgqESkWkeLGxkYPSvIPI/vi+mNLX/m+nKRIzp6ZzLObDjIw5LK6HOUlnoR+NTB11NeZwNhdM44cIyJBQCzQDCwB/lNEKoDvAP9HRO4a+wTGmEeMMYXGmMLk5MBZZbHCx/bFVYHn5qXZNLT38c7eBqtLUV7iSehvBmaISK6IhADXA+vGHLMOuNV9+WrgXTPsLGNMjjEmB/h/wI+MMf/tpdp93sg6+v4yXVP5n/Pyp5AeG8bTG/QMXX9x3NB399HfBbwB7AXWGmN2i8gDInK5+7DfMtyHXwbcDXxuWqf6vANNXaTEhBIR4j/TNZV/cTqEG5Zk8WFZE5+e4u5oyh48ShtjzHpg/Zjr7h11uRe45jiPcd9J1OfXfHFfXBV4rj19Kj9/Zz/PbKjk3ssKrC5HnSKdT2ehisMa+sr+pkSHcfHcNF7YUkVP/5DV5ahTpKFvkY7eAZo6+3UQV/mEm5Zk0d47yCs7xs7hUL5GQ98iI4O4uTpHX/mAotwEZqZE8fsNFbrBio/T0LfIkX1xtaWvfICIcPPSbHbVtLOjus3qctQp0NC3yMg6+tkJGvrKN1xxWgYRIU6e/GuF1aWoU6Chb5EDTd2kxoQRHuK0uhSlPBIdFszKoizW7ag9cja58j0a+hYpO9RBXrK28pVvufOcPIKdwi/e3W91KeokaehboH/Qxd66DuZlxlpdilInZEp0GDctyealbTWU68laPklD3wJjyEU3AAAQFklEQVSl9R30D7mYnxFndSlKnbA7z5lGSJCD/363zOpS1EnQ0LfAjupWAOZrS1/5oOToUG5ems1L22t0aQYfpKFvgU+q24iPCCYzPtzqUpQ6KXeeM43QIKe29n2Qhr4FdlS3Mj8zTjecVj4rKSqUW5Zl87K29n2Ohv4k6+kfYv+hTu3aUT7v62fnERrk5Bfv6EweX6KhP8n21LUx5DLMz9RBXOXbkqJCueWMbNbtqKXskLb2fYWG/iTb6T6FXVv6yh+sOiuP8GBt7fsSDf1JtrO6jZSYUFJiwqwuRalTlhgVyi3LcnhlZy37GzqsLkd5QEN/ku2sbmWezs9XfmTV2XlEBDv5ubb2fYKG/iTq6B2gvKmLBdq1o/xIQmQIX12ey6s769hc0Wx1Oeo4NPQn0a6adoxBl19Qfufvzp1GRlw4339pF4NDLqvLUcegoT+Jdh45E1e7d5R/iQgJ4vuXFlBS38ETuvSyrWnoT6KdNW1kxoeTEBlidSlKed1Fc1L4wqxkVr+1j/q2XqvLUUehoT+Jdla3skBb+cpPiQj3XT6HAZfhh6/tsbocdRQa+pOkpaufquYe7c9Xfi07MZJvnjudV3fW8eH+JqvLUePQ0J8kO2v0pCwVGO48J4/sxAjufXkXfYNDVpejxtDQnySfuAdx52Zo6Cv/Fhbs5P7L51De1MVjHxywuhw1hob+JNlR3UZeciQxYcFWl6LUhDt31hQumZvKL9/dT1Vzt9XlqFE09CfJJ9VtzNdWvgog37+0AIcI97+ig7p24lHoi8jFIlIqImUi8r1xbg8Vkefdt28UkRz39ReIyBYR+cT9/3neLd83HGrvpb69V+fnq4CSHhfOd744g7f3NvDqzlqry1Fuxw19EXECDwOXAAXAShEpGHPYHUCLMWY6sBr4ifv6JuAyY8w84FbgKW8V7kt0ZU0VqL56Zi4LMmP5/ku7aOzos7ochWct/SKgzBhTbozpB9YAK8YcswJ40n35BeB8ERFjzDZjzMhb/G4gTERCvVG4L9lZ3YpDYE66hr4KLEFOBw9du4Cu/iHuefETjDFWlxTwPAn9DKBq1NfV7uvGPcYYMwi0AYljjvkKsM0YE3Bv9ztr2piZEk14iNPqUpSadNOnRPPPF87i7b0NvLi1xupyAp4noT/eRq5j366PeYyIzGG4y+fOcZ9AZJWIFItIcWNjowcl+Q5jDDur27RrRwW0ry7P5fSceO57ZTd1bT1WlxPQPAn9amDqqK8zgbGjMkeOEZEgIBZodn+dCfwJuMUY8+l4T2CMecQYU2iMKUxOTj6x78Dmqlt6aO7q10FcFdCcDuGn1yxgcMjwr3/Ubh4reRL6m4EZIpIrIiHA9cC6McesY3igFuBq4F1jjBGROOA14B5jzEfeKtqXfKJn4ioFDC/RcM+X8nl/XyPPbao6/h3UhDhu6Lv76O8C3gD2AmuNMbtF5AERudx92G+BRBEpA+4GRqZ13gVMB74vItvd/6Z4/buwsR3VrYQ4HcxKjba6FKUsd9OSbM6Ylsh/vLZHT9qyiNjtY1ZhYaEpLi62ugyvueHRDXT2DbLuruVWl6KULVS3dHPx//uAuRkxPPu1pTgc4w0JqhMlIluMMYXHO07PyJ1AvQNDOoir1BiZ8RF8/9LZbChv5lf/W2Z1OQFHQ38C/aXkEJ19g1w0J9XqUpSylWsLp7JiYToPvbWPd0sarC4noGjoT6CXtteQHB3KGdOSrC5FKVsREX581XwK0mL49nPb+bSx0+qSAoaG/gRp6x7gLyWNXL4gHaf2WSr1OeEhTn5z82KCgxx8/ffFtPcOWF1SQNDQnyDrd9XRP+TiioVjT15WSo3IjI/gVzcuovJwN/+4Zjsul70mlvgjDf0J8tK2GvKSI5mbEWN1KUrZ2tK8RO69rIB3Sg6x+u19Vpfj9zT0J0BNaw8bDzRz5cIMRLRrR6njuXlpNtcWZvLLd8v48yd1Vpfj1zT0J8C67cOrVKzQrh2lPCIiPHjFXE7LiuO7f9jBntp2q0vyWxr6E+ClbTUsyoojKzHC6lKU8hmhQU5+fdNiosOCuOm3G9ld22Z1SX5JQ9/L9ta1U9rQwZWnaStfqROVEhPGmlXLCAtysPKRDWyvarW6JL+joe9lL22vIcghfHl+utWlKOWTcpMief7OZcRFhHDTYxvZXNFsdUl+RUPfi1wuw7rttZwzM5mEyBCry1HKZ01NiGDtncuYEhPKLb/dxEdlTVaX5Dc09L1o44Fm6tp6WaFdO0qdstTYMJ5ftYzsxAhuf2Izfyk5ZHVJfkFD34te2lZDZIiTC2anWF2KUn4hOTqU576+lFkp0ax6qphXd47dv0mdKA19L+kdGGL9rjoumpuqe+Eq5UXxkSE88/UlLMiM465nt3H/K7vpH3RZXZbP0tD3kv8tPURH76Auu6DUBIgJC+aZry/htjNy+N1HFVzz67/qJiwnSUPfS/60bWRFzUSrS1HKL4UGObnv8jn8+qZFlDd18aVffMDru/Ts3ROloe8Fh9p7+UtJI5fNTyfIqS+pUhPp4rlprP/WWeQlRfKNp7fyg5d30Tc4ZHVZPkMT6hQNuQzfXrMdhwNuWppldTlKBYSpCRH84RtncMfyXJ78+CBXPPxXNpYftrosn6Chf4p+/vY+Pi4/zA+vmEdecpTV5SgVMEKCHHz/0gIevaWQ1u5+rntkA3c+VUxFU5fVpdmahv4peH9fI7/8SxnXLM7k6sWZVpejVEC6oCCFd797Lt+9YCYf7G/igtXv8eCre2jr1k1ZxiPG2GvTgsLCQlNcXGx1GcdV39bLl37xAclRobz0zTN1mqZSNnCovZefvbWP54uriA0P5h/Om8H1p08lMjTI6tImnIhsMcYUHvc4Df0TNzjk4oZHN7Krto11dy1n+hTt1lHKTvbUtvOj9Xv5sKyJqNAgVixM54YlWcxJj7W6tAnjaej7/9vfBHjorX1sqmjm59cv1MBXyoYK0mN46o4itla28uzGSl7YUs0zGytZMDWOG4uyuHRBGhEhgRl/2tI/Qe+WNPDVJ4pZWZTF/71qntXlKKU80NY9wIvbqnl2YyX7D3USHRrE+bOncP7sFM6ZlUxMWLDVJZ4y7d6ZANurWrntd5tIiw3nT39/BmHB2o+vlC8xxlB8sIW1m6t4p+QQzV39BDmEJXkJnJ+fwhdnp/js5kca+l5UUt/OQ2/u4609DSRFhfCHb5xBblKk1WUppU7BkMuwrbKFt/ce4p29Dew/1AlAdmIEhdkJFOXGU5iTQF5SpE/sde3V0BeRi4GfA07gMWPMj8fcHgr8HlgMHAauM8ZUuG+7B7gDGAK+ZYx541jPZafQL2/sZPXb+3l1Zy1RoUGsOiuP25fnEhUAMwGUCjQHD3fxzt5DbCg/TPHBFpq7+gFIjAyhMCee07LimZUazayUaNJiw2z3RuC10BcRJ7APuACoBjYDK40xe0Yd8/fAfGPMN0TkeuBKY8x1IlIAPAcUAenA28BMY8xRz5m2OvTbugcoqW/nj1ur+ePWGkKcDm4/M4dVZ+cRF6EboygVCIwxfNrYxeaK5iP/qpp7jtweHRbErJRoZqVGM2NKFJnxEaTHhZMRH05suDXjA96cvVMElBljyt0PvAZYAewZdcwK4D735ReA/5bht8EVwBpjTB9wQETK3I/3saffiDe5XIau/kG6+obo7Bugo3eQqpYeSuraKanvoKSundq2XmD4bL9bl+Xw91+YRlJUqBXlKqUsIiJMnxLF9ClRrCwaXl6ltbuf0voO9jV0UNrQQWl9B+t21NLRO/g3940ODSIjPpy02DDiI0OICw8hLiKYuIhgYsOH/0WEBBEW7CA82ElYsJPQYAdhwU7Cg50ET/D6XZ6EfgZQNerramDJ0Y4xxgyKSBuQ6L5+w5j7TsjawyX17dz17DZcLsOQMQy5DC6XwWVgyBh6+ofo7Bsc975BjuEfcFFuAvlpMeSnRjM/M063PFRKHREXEcKSvESW5H22kq4xhqbOfmpae6hp6aG2tYea1h6qW3qoa+thX0MnbT0DR82esb48L42Hb1w0Ud8C4Fnoj9dxNbZP6GjHeHJfRGQVsAogK+vkFi0LC3IyMyUKhwhOh+AUweEQHAJOhxAeHERUWBBRoU6iQoOPXE6NCWf6lChCgnRFCqXUiRERkqNDSY4OZeHUuKMeNzDkoq1ngNbuAdp6+unuH6J3wEXvwNDwv0EXvf1D5EzCBBFPQr8amDrq60xg7J5lI8dUi0gQEAs0e3hfjDGPAI/AcJ++p8WPlpMUya9uXHwyd1VKqQkV7HSQFBVqi65iT5q3m4EZIpIrIiHA9cC6McesA251X74aeNcMjxCvA64XkVARyQVmAJu8U7pSSqkTddyWvruP/i7gDYanbD5ujNktIg8AxcaYdcBvgafcA7XNDL8x4D5uLcODvoPAN481c0cppdTE0pOzlFLKD3g6ZVNHL5VSKoBo6CulVADR0FdKqQCioa+UUgFEQ18ppQKI7WbviEgjcPAUHiIJaPJSOb5OX4vP6GvxGX0thvnb65BtjEk+3kG2C/1TJSLFnkxbCgT6WnxGX4vP6GsxLFBfB+3eUUqpAKKhr5RSAcQfQ/8RqwuwEX0tPqOvxWf0tRgWkK+D3/XpK6WUOjp/bOkrpZQ6Cr8JfRG5RkR2i4hLRArH3HaPiJSJSKmIXGRVjVYQkftEpEZEtrv/fcnqmiaTiFzs/rmXicj3rK7HSiJSISKfuH8PAmpVQxF5XEQOiciuUdcliMhbIrLf/X+8lTVOFr8JfWAXcBXw/ugr3ZuzXw/MAS4GfuXe7D2QrDbGLHT/W291MZPF/XN+GLgEKABWun8fAtkX3L8HgTZV8QmG//5H+x7wjjFmBvCO+2u/5zehb4zZa4wpHeemI5uzG2MOACObsyv/VwSUGWPKjTH9wBqGfx9UgDHGvM/wXh+jrQCedF9+ErhiUouyiN+E/jGMt7H7hGzObmN3ichO90fcgPgI66Y/+79lgDdFZIt7X+pAl2KMqQNw/z/F4nomhSd75NqGiLwNpI5z078ZY14+2t3Guc6vpiwd63UB/gd4kOHv+UHgIeCrk1edpfz+Z3+CzjTG1IrIFOAtESlxt4BVAPGp0DfGfPEk7ubR5uy+zNPXRUQeBV6d4HLsxO9/9ifCGFPr/v+QiPyJ4e6vQA79BhFJM8bUiUgacMjqgiZDIHTvBPTm7O5f5hFXMjzgHSg2AzNEJFdEQhge0F9ncU2WEJFIEYkeuQxcSGD9LoxnHXCr+/KtwNF6C/yKT7X0j0VErgR+CSQDr4nIdmPMRbo5O/8pIgsZ7taoAO60tpzJY4wZFJG7gDcAJ/C4MWa3xWVZJQX4k4jA8N/9s8aY160tafKIyHPAuUCSiFQDPwB+DKwVkTuASuAa6yqcPHpGrlJKBZBA6N5RSinlpqGvlFIBRENfKaUCiIa+UkoFEA19pZQKIBr6SikVQDT0lVIqgGjoK6VUAPn/HjpudFNCjE4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "kde = KernelDensity(kernel='gaussian', bandwidth=0.6)\n",
    "kde.fit(x.reshape(-1,1))\n",
    "\n",
    "s = np.linspace(np.min(x), np.max(x))\n",
    "log_pdf = kde.score_samples(s.reshape(-1,1))\n",
    "plt.plot(s, np.exp(log_pdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The KDE has learned the probability density function of our data sample\n",
    "\n",
    "Let's generate the samples from the KDE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = kde.sample(10000)\n",
    "plt.hist(m, bins=20, density=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KDE can learn handwitten digits distribution and generate new digits\n",
    "\n",
    "http://scikit-learn.org/stable/auto_examples/neighbors/plot_digits_kde_sampling.html"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}