Plugin Estimator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plugin Estimator\n",
    "This notebook we write code for a method to benchmark our random forest estimate against.\n",
    "We will be using the plugin estimator where we measure the parameters assuming we already know the distributions of X and Y.\n",
    "\n",
    "## Sampling Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAERCAYAAABhKjCtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XlcVNf9+P/XnR1ZZBsWWWRTcANUFKMmGjWbrW3MxyZplqb5uCQ1SdMs5pGYtGkxzdIa86tmqTF+TLSptWnWbxoTYxZN1EQFdxEQlEUFBgQFHIaZuff3B2EiggIycIfhPB8PHg+5c++57zNX5j33nHPPkRRFURAEQRCEC2jUDkAQBEHwTCJBCIIgCO0SCUIQBEFol0gQgiAIQrtEghAEQRDaJRKE4KLWgDZvHEjXk3XyxvdL8EwiQfQhd955J8nJya6fESNGMGnSJH73u99RWFjYat/HH3+cm2++udNlL1++nPfee++S+5xfZllZGcnJyWzdurXrFTlPfn4+CxYscP3+3nvvkZycjM1m61a5HbnwvUxJSWHMmDHceuutrerU1fexqamJP/3pT2zfvt0tcV74fuzevZtFixa5Xl+xYgWTJk26rDLT0tKwWq3t7nPddde55fqqqaSkhPvvv58rrriCzMxMFi5cSGlpqdph9Sk6tQMQumbixIk8+OCDANjtdsrLy1mzZg1z5sxh/fr1pKSkALBw4cIufci+8sor/PGPf7zkPl0tszM+++wzDh8+7Pp96tSpbNiwAYPB4NbztOf891JRFOrr61m3bh333nsv77zzDiNGjOhymZWVlfzzn//k6quvdkuMF74f//nPfygpKXFL2TabjW3btjFjxoxW2/Pz8zl+/LhbzqGWhoYGfv3rXxMYGEhWVhaSJPHqq69y55138vHHH+Pn56d2iH2CSBB9TGBgIOnp6a22zZgxg5tuuomnn36aDRs2ABAbG+v2c/dEmRcKDg4mODi4x88D7b+X48aN46qrrmLDhg1kZWX1ShyX0pPvR3p6Ops3b26TIDZt2sTQoUPJz8/vkfP2hk2bNlFZWcmGDRswm80ApKWlMWXKFD7//HNmz56tcoR9g2hi8gI+Pj7MnTuXvXv3upqaLmwa+eqrr5g9ezapqalMnjyZZ555hsbGRgCSk5MB+OMf/8idd94JwLRp03jppZeYPXs2aWlpfPjhh+02txQVFXHrrbcyatQoZs+ezY4dO1yvtdf8sXXrVpKTkykrK2PFihW8/PLLVFVVkZyczPfff9+mSUWWZd5++21+8pOfkJqayg033MC7777rKq+lqWvLli3cddddpKamcvXVV/Ovf/3rst5Lk8lEXFwcJ0+ebPf1+vp6/vznP3P11VeTmprKbbfdxp49e1yxTJ8+HYD58+fz+OOPtzn+zTffZPTo0TgcDte2G2+8keuuu871uyzLZGZm8s4777R6Px5//HHef/999u3b53oPW7z33nvMmDHDFdPRo0c7rOs111zDV199hdPpbLV906ZNreJpUVhYyLx580hPTyczM5MlS5a4/g9B813YG2+8wcyZMxk5ciQZGRncf//9VFRUuPZJTk7mww8/5P777yc9PZ1Jkybx8ssvXzLO85sCL/y5WLNocHAwd999tys5AJjNZvz8/Dhx4kSH743QTCQIL3HFFVcAsHfv3javFRcX89vf/pbx48ezatUqHnzwQd59911WrFgB4LrrmDdvHk8//bTruNWrV/Pzn/+cpUuXMmHChHbPu3TpUjIzM3n55ZeJiIhg/vz5FBQUdCrmX/ziF8yZM4fAwEA2bNjQbpPOX//6V5577jlmzpzJq6++yuTJk1m8eDFvv/12q/2eeOIJrrjiClauXMmYMWN4+umnO/UheSGHw8GJEyeIiopq85rT6eTuu+9m48aN3HfffSxfvpwBAwbwq1/9ikOHDhEWFub6sHviiSdYuHBhmzImT57MuXPnXM1qZ8+eJS8vj+PHj3P69GkADh8+TG1tLVdeeWWrYxcuXMiUKVMYMmQIGzZsICwsDIDTp0+zatUqHn30UZYuXcrJkyd5+OGHO6zrjBkzOHPmDLt373ZtKykpoaioiGnTprXat6Kigttvv52GhgaWLVvGY489xn//+18ee+wx1z6rVq3i5Zdf5vbbb+f//u//ePjhh/nuu+/4y1/+0qqsZ555htjYWF577TVmzpzJihUr2LJly0Xj3LBhw0V/pk6d2u4xU6ZM4ZFHHmm1bc+ePZw5c4aEhIQO3xuhmWhi8hIhISEAVFdXt3nt4MGDNDU1MXfuXMLCwsjMzMRgMLi+xbY0s0RHR5OUlOQ6buTIkfz617++5HlvuukmHnroIQAmTZrEjBkzWLduXaeaZyIiIoiIiECn07Vp6oHmD75169Zx33338Zvf/AZo/oBtaGhg+fLl3HLLLa3iuPfee4HmpoRPP/2UrVu3tqrPhRRFcb0HsixTXl7O3//+d6qrq5kzZ06b/b/66iv279/PP//5T8aOHQvAlVdeycyZM3nllVd49dVXGTZsGAAJCQntNsklJSURERHBzp07SU1NJTs7m8jISM6ePUtOTg4zZsxg+/btDBkyhIiIiFbHxsbGEhwcTG1tbav3S5Zl/va3vzF06FAALBYLWVlZ1NXV4e/vf9H6R0REMGrUKDZv3kxmZibQ3CeUmZlJQEBAq33feustJEli1apVrvb7QYMG8etf/5rc3FyGDRtGZWUlv/3tb7n99tsBGD9+PEVFRXzxxRetypo8ebIrsUyYMMF1raZMmdJunO393+iqhoYGnn76aWJiYto0qQkXJxJEP5CamorBYODmm2/mpz/9KdOmTWPWrFloNJe+gYyPj++w7Guuucb1b51Ox+TJk8nJyel2zAD79+/Hbrdz/fXXt9o+c+ZM3n33XYqKihgwYADQnBRaDBgwgICAAM6dO3fJ8jdu3MjGjRtbbQsODiYrK4tRo0a12T87Oxuz2exKDgBarZbrrrvOdRfWGZMmTWLXrl3MmzePXbt2MW7cOMrLy1sliAvvHi7F19fXlRygOdEDHSYIaL5+//rXv3jyyScB+Pzzz7npppva7Ldr1y7Gjh2LyWRyJdWW37///nuGDRvGU089BUBVVRWFhYUcPXqUnJwc7HZ7q7LOv1aSJBEeHn7Ja3V+c9yFtFotkiRdso4NDQ3ce++9lJaWsnbt2l4ZAOEtRILwEhaLBcDV7HC+mJgY1qxZw9///nfefPNNVq1aRVRUFE8//fRFv7XBj3cll3LhPkFBQdTV1XUx+vadOXMGgNDQ0HbPWV9f70oQRqOx1T4ajabD5wUmT57M7373O9f+AQEBREdHX/QD5+zZs21igeak0tDQ0Ika/XjeP/zhDzidTnbt2sXNN9/MqVOn2L59O42NjWRnZ7ca+tsRHx+fVr+3xC/LcofHzpgxgxdffJHc3FyCgoI4ePAgr776aqu+BYDa2lr279/fbjNgZWUlAEePHuXJJ59k7969+Pr6MmLECIxGY5vr0NVrdanRZM8991y7Ca3F6dOnmT9/PkVFRbz66qvtJn7h4kSC8BI7d+4EYPTo0e2+npGRwRtvvEFDQwPffvstf//733nooYf47rvvuvWN6uzZs61+r66udo26kSSpTQdoR9/qzzdw4ECg+Rvp+d+Eq6qqWr1+uQICArr0gREQEOA69/mqq6u7FMvEiRNpaGggJyeHw4cPM2bMGCoqKli1ahXbt29Hp9ORkZHR6fK6IyEhgaSkJL744gsCAwMZM2YMoaGhrTrAAfz9/bnuuuuYP39+mzJCQ0ORZZnf/OY3REREsHHjRuLj45Ekib/+9a/dHpb7n//856KvtdwttaeiooJf/epXnD59mtWrVzNmzJhuxdEfiU5qL9DU1MSaNWsYN24cgwcPbvP6Bx98wPTp07Hb7fj6+rr+0BsaGqivrwfosLnpYs5/IKypqYmtW7e6mmAGDBjA2bNnWyWF7OzsVsdf6rypqano9Xo+/fTTVts3btxIUFAQcXFxlxXz5Ro7diwWi6VVHZxOJ5s2bXK1k2u12g7LCQwMZOTIkbzxxhsEBASQmJhIeno6siyzcuVKxo8ff9GkfbnX6VJmzJjBl19+yebNm7n22mvb3WfMmDEUFRUxfPhwRo0axahRowgPD2fZsmWuDvaSkhJ++ctfkpCQgCRJyLLMjh07OnUncykt52vvJygoqN1jbDYb8+bNo7a2lnXr1onkcJnEHUQfU1tb6xqp5HA4KCsr4x//+Afl5eW89NJL7R4zduxYqqqqePjhh7n11luxWq289tprjBkzxvVtPyAggJ07dzJ69GjXw3ad8c9//pPg4GASEhJ46623sFqt3H333UBzU8oLL7zA73//e37xi1+wZ88ePvnkk1bHBwQEcObMGb7++us2dz/BwcHcdtttvPLKK8iyTHp6Olu3buW9995j8eLFnfowdqepU6cycuRIHnzwQR566CFCQ0N5++23OXHiBH/9618BXHc633zzDdHR0RcdMTN58mReffVVVx/OgAEDGD58OHv37uX3v//9RWMICAigpKSEHTt2XPRusauuvfZaVq5ciU6n4/nnn293n7vuuov333+f++67j1tvvRWHw+EaopySkkJgYCCRkZGsXr0aX19fZFlm/fr1HD58GJPJ5JY4u+LNN98kPz+fhx56iMbGxlaj+8LDw4mMjOz1mPoicQfRx2zfvp1bbrmFW265hbvvvpuXX36ZuLg43nnnHYYMGdLuMTExMbz22mtUVFRw//3388QTTzBkyBCWL1/u2uc3v/kNX3/9datpHDrjySef5MMPP2ThwoXU19fz1ltvuYaIDh06lKysLPbs2cOCBQvYvXu364O0xcyZM0lKSuL+++/n22+/bVP+448/zn333cd//vMf7r33XrZv386zzz7Lr371qy7F6Q46nY7Vq1dz1VVX8Ze//IUHH3wQq9XK2rVrSU1NBcDPz4///d//ZcOGDbzwwgsXLWvy5MkArTq8W5qVLtVBffPNN+Pn58eCBQtaPYHeHSNGjGDQoEEMGzaszcipFjExMfzjH/+gqamJ3/3udyxevJjIyEjWrl1LUFAQkiSxfPlyNBoNv/3tb/nDH/6An58fy5Ytw2q1kpeX55ZYO6tl5NRLL73k+ntp+blwiLRwcZJYclQQBEFoj7iDEARBENolEoQgCILQLpEgBEEQhHaJBCEIgiC0SyQIQRAEoV0iQQiCIAjt6rMPytXUNCDLCiEhflRX16sdjtt4W33A++ok6uPZRH3a0mgkgoJ8u3xcn00Qsqwgy4rr397E2+oD3lcnUR/PJurjHqKJSRAEQWiXSBCCIAhCu0SCEARBENolEoQgCILQLpEgBEEQhHaJBCEIgiC0q88OcxX6B/lMBU5LEfLZSrThQ9BGDkXSiP+2gtAbxF+a4JEUWwON2/6B4+iO1i8YfTFNugN90hXqBCYI/YhIEILHcVYVY/3sbyjnajGMnoUuYTwa/xAcJ3Kx7/+Uxi9XIp8+gWHcTUiSaCUVhJ4iEoTgUeQ6C9aNL4JWz4CfP4U27Mc1nfXxY9HFpmHbtpamvR+jyE5ME25RMVpB8G4iQQgeQ7E1YN34EorTwYCfPo42aFCbfSStDuOVd4Okxb5/I7rIZHSD01WIVhC8n7g/FzxG4zdvIp+twOfaB9pNDi0kScJ4xS/RhAzG+vUq5LqqXoxSEPoPkSAEj2A/lo2jaBeGsTeiGzSsw/0lnQGfGQtBdtK4bV0vRCgI/Y9IEILqFFsDtm/XogmJwZB2Q6eP0wwMxzB6Fs6SfThOHO7BCAWhfxIJQlCdbde7KI11mKbM7fIzDoaR1yD5hWD77l8ostxDEQpC/yQShKAqubYce+7X6IdNRRsa1+XjJZ0B4/hfIFeX4CjY5v4ABaEf61SCsNvtZGVlkZmZSWZmJkuXLkW+yLe13Nxc7rjjDtLT05k2bRqrVq1q9Xp5eTkLFixg9OjRzJgxg48++qj7tRD6LNvud0GrxzDm55ddhi4xE01oHLa9H4u7CEFwo04liGXLlrFt2zZef/11XnrpJT744APeeOONNvvV19czd+5chg4dykcffcRTTz3Fa6+9xvvvv+/a57777sNgMPDOO+8wb948Fi9eTE5OjvtqJPQZTsux5o7p1OvRDBh42eVIkoQhbSbKmQocJXvcGKEg9G8dJgibzcb69etZvHgxaWlpTJw4kUceeYS1a9eiKK2XwTt58iRXXHEFTz75JLGxsUybNo2JEyeyc+dOAHbt2kVeXh7PPvssSUlJ3HrrrcyaNYt168QolP7ItutdJJM/htTru12WLn4skn8o9n2fuiEyQRCgEwkiNzcXq9VKRkaGa1tGRgYWi4WysrJW+w4dOpQXX3wRrVaLoijs3r2bXbt2MXHiRAD27NlDcnIyAQEBrcrau3evu+oj9BFOy3GcZQfRp16HZPDpdnmSRoth1HU4KwpwVhx1Q4QCgMMpU1Vr5ey5JhxO0XzX33Q4ZKSiogI/Pz98fX1d28xms+u1mJiYdo8bN24cdXV1XH311cycOdO1f1hYWKv9QkNDqaiouOwKCH1T077/gt4Hw/BpbitTn3wltt3v03TgM3zCk9xWbn+1v7CKf24uoLLGCoDRoOWX04dwZWokkiSpHJ3QGzpMEFarFaPR2GqbwWAAoKmpqd1jZFlmzZo1lJeX88c//pHnnnuOp5566qJlOZ1OHA4HOl3nhziGhPi5/m02+3f6uL7A2+oDrevUVH2SuqLdBE6cTXBUuBvP4o8mfRpnd39K8AAZre/l92t0xNuu0fn1ccoKf/tXDl9llxFl9uXem1JBUdh+4BRvbjxC0ak6HrglHZPBc2fq8ebr05s6vMImk6lNImj53WQytXuMRqNh1KhRjBo1ioaGBp566ikee+wxTCYTtbW1bcrS6/VdSg4A1dX1yLKC2eyPxVLXpWM9mbfVB9rWqXHrO6DVYU+Y4va6OmMnwM6PKf/ucwyp17m17Bbedo0urM/6zQV8lV3GTyfG8bNJcei0zS3RGUNC2fh9Me9tKUJC4a7rU9QK+ZK8/fpcDo1GavWlutPHdbRDREQEdXV1WK1W1zaLxQJAeHjrb39lZWVs2bKl1bYhQ4Zgt9upr68nIiLCdez5ZV3Y7CR4L9l6FnvBdvRDJ3dr5NLFaIOj0ZgTsOdtbTOIQujY57tK+Xx3KddkxHDTVQmu5ADNHzI/uSKO6zNj2bL3JHsLxBxY3q7DBJGSkoKPjw/Z2dmubbt37yYsLIyoqKhW++7bt48HH3yQhoYG17ZDhw4REhJCcHAw6enp5OXlUV9f73o9Ozub0aNHu6MuQh9gz/0anA70I6/psXPoU65CrjmBbDnWY+fwRsXldfzrywLGDDVzy7SL9+HceGUCsWF+rNmYy5mG9puZBe/QYYIwmUzMmTOHrKwscnJy2LFjBy+++CJ33XUXALW1tdTVNd/+TJ06leDgYBYvXkxRURGbN29m2bJlLFy4EGgesZSQkMCjjz5Kfn4+GzZs4OOPP+bOO+/swSoKnkKRHdgPf4k2euQlZ2vtLn1iJmgN2I9s7bFzeBunLPPmp0cIGGDgf2emoNFcvBNar9Mwf9ZwrDYn728t7MUohd7WqQflFi1aRGZmJvPmzePhhx9m9uzZzJ07F4AHHniAP//5zwD4+vqyevVq6uvr+Z//+R+ysrKYO3cud9xxR/PJNBpefvllbDYbc+bM4Y033uC5554jPV3M598fOIp2N68S14N3DwCSwQdd/FjsRTtRnI4ePZe3+GJ3GcXlddx2zVAGmPQd7h9l9uOqtEi2HSjn9NnGXohQUIOk9NGGWtFJ3Xe01KnhgyUotnp8b36ux5cKdZTsxfrp/4fPdb9z+4JC3naNNAYd9zz3BcmxgTw4J7XTQ1irzlh5YuV3TBsTzS9nDOnhKDvP266PR3dSC4I7OKuKkSsLMQyf1ivrSGujRoLRF3vh9z1+rr5u/aY8nLLM7dcM7dLzDaEDfZgwPJwte09w9pzoi/BGIkEIvcKe+zVo9eiHTu6V80laHfq4sTiK96A4xIfXxZSfPsfnO0uYmh6FObDrT7TfMGEwdofM5t2lPRCdoDaRIIQeJ9us2I/uQJc4Hsno2/EBbqJLzAR7I47S/b12zr7m/a1FGHQafjIx7rKOHxTqS1pSKFv3ncIpZtL1OiJBCD2u/tA3YG/EMOzqXj2vdlAKkskfR+HOXj1vX1FcXseuI5X87KpEBvoaLrucSaMiONvQxKFjNW6MTvAEIkEIPe5sziY0ITFowhJ79bySRosuPgNHyV7RzNSO/7f9OD5GHbOndm/eqtTEUHxNOrYfPOWmyARPIRKE0KOcVcdpqjiGPmWqKhO86eLGgKMJ58ncXj+3Jyuz1JOTb2HG2Gj8fDoe1nopep2GzOHh5ORXca7R7qYIBU8gEoTQo+xHvkHSGdAnTVDl/NpBKaA34TguFhI63yc7ijHqtVwzrv3ZmLtq0qhIHE6ZXUcq3VKe4BlEghB6jOJown50B77Jmb3aOX0+SatHFzOquZlJEZ2oABU15/g+t4KrR0d1++6hRVyEP5EhA9h2sNwt5QmeQSQIocc4judA0zn806erGodu8GiUc7XIluOqxuEpPtlRjFaj4brx7rl7gOZlXycMD+do2Rlq621uK1dQl0gQQo+x521F8jdjGjxC1Th0MakgaZoTVj93+mwj2w+Wc2VaJAP9jB0f0AWjhzYvJCZmefUeIkEIPUKur8Z54jD6oZN75cnpS5FMfmgjhuIoFv0Qn+0sRVHghvGxbi87KtQXc6CJPSJBeA2RIIQeYS/YDoB+6ESVI2mmG5zePAV4fbXaoajm7Lkmtuw7wYQR4YRexlPTHZEkidFDzOQWn8ZqE5MkegORIAS3UxQFR/42tJHJaPzNaocDgDYmFQBH6QGVI1HP5t1lNNllbpgwuMfOMXpIKA6nwqFjp3vsHELvEQlCcDvZcgz5TDm6IZ5x9wCgCYxE8g3G2U8TxLlGO19klzFmqJmo0J4bUZYUPRA/Hz05BZaOdxY8nkgQgtvZ87c1T8yXME7tUFwkSWoe7nriMIrc/5o/NmeXYbU5mHWZcy51llajIS0phP1Hq3E4xbDivk4kCMGtFKcDR+H36AaPRjIMUDucVrQxo8BuxVnRv1ZBs9ocfL6rlPSkUAZH+Pf4+dKTQjlnc1B08myPn0voWSJBCG7lLDuAYqtH70HNSy10UcNB0uAsO6h2KL3qi+wyGhodzJoU1yvnGzY4CI0kcVD0Q/R5IkEIbmU/+j2S0Q9tzEi1Q2lDMgxAG57UrzqqzzU6+GxnCamJIcRHBvTKOQeY9CQMCuDQsf47YsxbiAQhuI1it+EozkGXkIGk0akdTru00SORq44jW/tH88enO4tpaHRw45XxvXreEfHBHD9VR71VTN7Xl4kEIbiNo3gPOJrQJaozMV9n6KKbn+p2njyiciQ9r7bexqadpYwfFkZcRO/cPbQYGR+MAhw+LpqZ+rJOJQi73U5WVhaZmZlkZmaydOlS5IusHnX06FHmzp1LRkYGU6ZM4YUXXsBm+3Fulueff57k5ORWP1lZWe6pjaAq+9HvkHyD0EYOVTuUi9KExoHeB+fJw2qH0uM++vYYTlnhpqsSev3ccZH+DDDqRD9EH9epdoBly5axbds2Xn/9dRoaGnjssccICAhgwYIFrfZraGhg/vz5jBs3jg0bNmCxWHjyySdxOp0sXrwYgIKCAhYuXMhtt93mOs7Hx/1PdQq9S2msx1l2AP3Ia1SfWuNSJI0WbWQyjhPevT7EiaoGtu47xdWjowgL6v3RZFqNhmFxQRw6dhpFUVRZC0Tovg7/km02G+vXr2fx4sWkpaUxceJEHnnkEdauXYuiKK323bFjB3V1dTzzzDMkJiYyYcIEHnzwQT766CPXPoWFhYwYMQKz2ez68fPzc3/NhF5lP54NshN9YqbaoXRIFzUc5WyF1067oSgK//gsDx+jllmT41SLY0R8MDV1Nk5Wn1MtBqF7OkwQubm5WK1WMjIyXNsyMjKwWCyUlZW12nfUqFG88sorGAw/rm8rSRL19fUoikJDQwOnTp0iPr53O8yEnuc4thvJ39zchOPhtFHDAHCe8M5mpu0Hy8krreV/piYSMODy15rurhFxwQDkin6IPqvDBFFRUYGfnx++vj8+nm82m12vnS88PJzMzB+/QTqdTtatW0dmZiaSJFFY2PyA0ptvvsnUqVOZOXMmq1evvmh/htA3NDcvHUYXn9EnmhI0QdFIJn8cXrgMaUOjnX9/dZTEQQFclTZI1VjMgT6EDjRxpKRW1TiEy9dhH4TVasVobD1vfMsdQlPTpReCX7JkCUeOHOHf//430Ny8pNFoiIqKYuXKlRw8eJBnn30Wp9PZpj9D6DscxXtAcXrU1BqXIkkS2kHDcJ447FXt44qisO6zPBqsDu68JRmNB9QrJTaIPQUWZEXxiHiErukwQZhMpjaJoOV3k8nU7jFOp5M//elPvPvuu/ztb38jJSUFgBtvvJGpU6cSFBQEQHJyMrW1taxbt67LCSIk5Md+C7O556cP6E19rT6nvshBNzCM8OGpF/2w9bQ6nU0eTVXRTgK1dRhCorp8vKfVB2DzzhJ25lZyxw0pjB3ZtbuHnqrPuJGRfHvgFPVNMonRgT1yjvZ44vXpDrXq02GCiIiIoK6uDqvV6hptZLE0z9QYHh7eZn+73c6jjz7Kl19+yfLly5k+/cflJiVJciWHFklJSVgsli5/k6uurkeWFcxmfyyWuk4f5+n6Wn0UWwPWY/vRj7yGqqr6dvfxxDrJ/s39YJbDezCkdO0ZAU+sT8Xpc/z9vf0kxwQydVRkl+LryfpEBTV/idyx7wQBRm2PnONCnnh9usMd9dFopFZfqjt9XEc7pKSk4OPjQ3Z2tmvb7t27CQsLIyqq7TevP/zhD2zZsoWVK1e2Sg4AL730ErfffnurbYcPHyY+Pt5rbvP7G0fx3ubRSwnj1Q6lS6SB4Ug+AThP5akdSreda3Sw4r0D6LQS82cNR6PxnL+l4AATYUE+HCmuUTsU4TJ0mCBMJhNz5swhKyuLnJwcduzYwYsvvshmJtx3AAAgAElEQVRdd90FQG1tLXV1zdlty5YtvPfeeyxatIghQ4ZgsVhcPwDTpk1jz549vPLKK5SUlPDhhx+yatUq7rnnnh6sotCTHMezkXyD0Jj71sg0SZLQRib3+QThcMq89sEBKk6fY+GNIwkOaL/ZV00psUHkl9XiFINR+pxOPSi3aNEibDYb8+bNw2g0MmfOHObOnQvAAw88QFRUFM8//zwbN24EICsrq83T0fv37yctLY0VK1awYsUKVq5cSVhYGIsWLWLWrFlurpbQGxSHDUfpQfTJV/bJO0BtRDKOol3IdVVo/EPVDqfLFEXhH5vyOHS8hl/fkMKwH4aVepphg4PYuu8kJRX1vTZhoOAeknLh0259hOiDUJ/9eDaNm1bg85PHmqfSvghPrZOzupRz7/4e09ULujQ9uSfUR1EU/rm5gC+yy/jJFYP5nymJl11WT9fnTL2Nh17exi+mJvbocqctPOH6uJNH90EIwsU4juWA0dej5166FE1wFBgG9LlmJkVR2PDlUb7ILuPacTGqzLXUFQP9jESGDCCvVDwP0deIBCFcFkV24CjZiy423WOn9u6IJGnQRgztUwlCVhT++XkBm3aVMn1sNLdMS+oTzXtDYwIpKKtFlvtkg0W/JRKEcFmcp/LB1oAubozaoXSLLjIZ+Uw58jnP/3YrywpvbTzCFzllXD8+lttmDOkTyQEgOSYQq81JaWX7Q6EFzyQShHBZHMV7QKtHF+15K8d1RUvzmLO8QOVILk2WFdZ8kss3+0/xs0lx/OLqxD6THKD5DgIQzUx9jEgQQpcpioKjeC/aqOFIemPHB3gwTchg0OpxVhxVO5SLkmWFNRtz2XawnBsnx3PjlQl9KjlA8/MQoQNN5IsE0aeIBCF0mVxzEqXOgi42Xe1Quk3S6tCa43FWeO4dxL++LGDbgXJ+Pjmen03uW8+bnC85NpD80to2ywQInkskCKHLHCV7ANAN7vsJAkAbMQS5qhjFcenJJ9WweXcpm3eXcU1GDD/vw8kBmpuZ6q12TlY1qB2K0EkiQQhd5ijeiyZ0MBrfoI537gO04UkgO3FajqkdSiv7jlax/osCRg8J5ZZpSWqH023Jsc3/X0QzU98hEoTQJXJjHXJFoVc0L7XQhDd/+HpSM1NtvY3V/80lxuzHglkjPGp+pctlHmgiyN8oOqr7EJEghC5xluwHFHSDR6sdittoTP5oBkZ4zEgmRVH4v//m0mR3cs/PR2A09M4sqD1NkiSGxoh+iL5EJAihSxwl+5AGBKIJ7fkpE3qTNmIIzoqjHvHB9WXOCQ4eO83N05KIDPHt+IA+ZGhMILX1TVhqrWqHInSCSBBCpymyE0fZAXQxF18YqK/ShCeBrQH5zClV4zhTb+M/WwoZGR/M1aO7vpCRpxPPQ/QtIkEIneasOApNVrSxqWqH4nbaiCEAyOXqPg/x/jdFOBwyt10z1OuSMMCgkAH4+egpKD2jdihCJ4gEIXSas3Q/SFp0USPUDsXtNAMjmifuqyxSLYaSijq+2XeKaWOiiQgeoFocPUmSJIZEDxQjmfoIkSCETnOU7EMbORTJ4KN2KG4nSRq0YQk4KwtVi+HfXx1lgEnHzybHqRZDb0iOCaSy1kpNnU3tUIQOiAQhdIpcX418ugxdjPc1L7XQhiUi15ShNPV+B2pBWS2Hj9cwa2IcviZ9r5+/Nw35oR9C3EV4PpEghE5xlB4A8Mr+hxbasERQFJxVx3v93J9+X4KvSceUdO/rmL5QbLgfRoOW/DKRIDydSBBCpzhL9iH5haAJHKR2KD1GG9a88I6zonebmU5VN7CnoIrpY6O95pmHS9FqNAyJGkh+iUgQnk4kCKFDitOB42SuVw5vPZ9k8kMaGIHcy/0Qn35fgl6nYdrY6F49r5qGxgRyoqqBunOeN/+V8CORIIQOOSsKwN6INmaU2qH0OG1YIs7Kwl57YO5MvY0dh8qZnBpJwABDr5zTE7Q8D1FQJoa7ejKRIIQOOUsPgEaLbtAwtUPpcdrwRBTrWZT6ql4537cHTuFwKlyTEdMr5/MU8ZEB6HUa8kQzk0frVIKw2+1kZWWRmZlJZmYmS5cuRZbldvc9evQoc+fOJSMjgylTpvDCCy9gs/04nK28vJwFCxYwevRoZsyYwUcffeSemgg9xlG6H22Edw5vvVBv9kMoisK3+08xNHqg1z73cDF6nYbEQQHkldaoHYpwCZ1KEMuWLWPbtm28/vrrvPTSS3zwwQe88cYbbfZraGhg/vz5hISEsGHDBl544QU2bdrEiy++6Nrnvvvuw2Aw8M477zBv3jwWL15MTk6O+2okuJXcUPPD8Fbvb14C0ARHg9bQKw/MFZSdoaLGypVp3tvxfynJsUGUVtRzrtGudijCRXSYIGw2G+vXr2fx4sWkpaUxceJEHnnkEdauXdumnXbHjh3U1dXxzDPPkJiYyIQJE3jwwQdddwm7du0iLy+PZ599lqSkJG699VZmzZrFunXreqZ2Qrc5SvcDoPXi5x/OJ2l0aEMH47T0fIL4Zt9JTAYtGclhPX4uTzQ0JhAF0Q/hyTpMELm5uVitVjIyMlzbMjIysFgslJWVtdp31KhRvPLKKxgMP3a2SZJEfX09iqKwZ88ekpOTCQgIaFXW3r173VEXoQc4Sw8g+QahCfL+8fktNGEJzSvMyY4eO4fV5mBXXiWZw8P7xdDW9iQOCkCnlcTEfR6swwRRUVGBn58fvr4/TjtsNptdr50vPDyczMxM1+9Op5N169aRmZmJJElUVFQQFtb621JoaGibcgTPoMhOHCcOoYsZ5dXDWy+kNceD0458+kSPnWP3kUqa7DKTUyN77ByezqDXEh8ZIDqqPZiuox2sVitGo7HVtpY7hKamS49hXrJkCUeOHOHf//73JctyOp04HA50ug7DcQkJ8XP922z27/RxfYGn1KexNJf6JitBI8bj182YPKVOnWHXjaL0SxhgPUmAeWS7+3S3PnsLq4kM8SUzNcojkq9a12d0Sjj/+bIAX38TA9w4xUhf+v/WGWrVp8NPZJPJ1CYRtPxuMpnaPcbpdPKnP/2Jd999l7/97W+kpKS49q+tbf1toampCb1e36XkAFBdXY8sK5jN/lgsdV061pN5Un1s+78DScM5/wSs3YjJk+rUGYoyAMnox5miw9hirmjzenfrU2+1s6+giuszY6mqqu9OqG6h5vWJCR2ALCt8t+8EoxJC3FJmX/v/1hF31EejkVp9qe70cR3tEBERQV1dHVbrjxOYWSwWoLlJ6UJ2u52HH36Y999/n+XLlzNjxoxWZbUce35ZFzY7CZ7BUXYAbXgSkqF/DcGUJAlNWALOymM9Un5OvgVZURiXIv7fJ0UNRKuROFIshrt6og4TREpKCj4+PmRnZ7u27d69m7CwMKKi2nZc/uEPf2DLli2sXLmS6dOnt3otPT2dvLw86ut//NaUnZ3N6NHes76xt5DPnUGuKu4XT0+3R2uOR6450SMzu+4+Uok50ERseNe/0Xkbo15L4qAAckWC8EgdJgiTycScOXPIysoiJyeHHTt28OKLL3LXXXcBUFtbS11d8+3Pli1beO+991i0aBFDhgzBYrG4fqB5xFJCQgKPPvoo+fn5bNiwgY8//pg777yzB6soXA5n2UEAr57e+1KaH5hTcFYVu7Xcequd3OIaMlLCPKLvwROkDA6iuKKOc409N2pMuDydavhftGgRNpuNefPmYTQamTNnDnPnzgXggQceICoqiueff56NGzcCkJWVRVZWVqsy9u/fj9Fo5OWXX+b3v/89c+bMITw8nOeee4709HQ3V0voLkfpASSfADQh/WsKiBYaczwAzsoidINS3FbunnwLTlnpt88+tCclNoiPth0nv7SW9CGhaocjnKdTCcJoNLJkyRKWLFnS5rXzH3J7/vnnef755y9ZVnR0NGvWrOlimEJvUmQZR9kBdLHpSFL/nK5L4xOA5G9GdvMDczn5FkICTMRFeNcom+5IjApAp9VwpKRGJAgP0z//+oVLki1FYGtAF5umdiiq0prjcVrc11Ftszs5XNz8ISial36k12lJihL9EJ5IJAihDUfJPpAkdNEj1A5FVVpzPEp9NbL1rFvKO1Jcg90hk5bonuGc3iRlcBCllfXUW8W8TJ5EJAihDUfpfrThQ5CMvh3v7MVa+iFkN91F7C+sxqjXkhwb6JbyvMmwwUEAYrirhxEJQmhFPlf7w/DW/jl66Xza0MGA5JZmJkVR2FdYxfC4IPS6/jn30qXERwZgMmg5fPy02qEI5xEJQmjFWXoAAF2sSBCSwQdNUKRbEkSZpYHTZ22kJYlO2PbotBpSYoM4eEwkCE8iEoTQiqN0P9KAQDTB/XN464U05nhky7FuL0G672jzCnWpov/hokbEB1N1ppHKmnNqhyL8QCQIwUVxOnCUHkQXkypG2fxAa45vXoK0oXvfbPcXVjM4wp9AP2PHO/dTI+KDATh0XPRDeAqRIAQXZ3k+2K1oB4sHF1tozT8sQdqNFebONdopPHnGbZPReavwIB9CAowcEs1MHkMkCMHFUbwHtDp0Uf17eOv5NCExoNF2ayRTbnEtigIj4oLcGJn3kSSJEfHB5BbX4LzImvdC7xIJQgCaR9k4iveiHTQcSS+aQVpIWj2akNhudVQfOn4ao0FLYtRAN0bmnYbHBWO1OTh+ynum6+7LRIIQAJBrTqLUWdANFjPrXqj5ierjKMrlfas9dKyaYbFB6LTiz60jw+OCkUCMZvIQ4n+sAICjZA8AOtH/0IbWHA92K/KZ8i4fW1lzDkttI8NF81Kn+PnoSRgUwP7CarVDERAJQviBo3gvmtDBaHzFB9mFXE9UX8YCQi0jclpG6AgdS00M4fips5xpuPSSxkLPEwlCaH56uqJQNC9dhCZwEOiMl9UPcfjYaYIDjEQE969V+bojNTEUBTgg7iJUJxKEgON4DqCgix+rdigeSdJo0JrjcHZx6m+nLHO4uKa5XV08V9JpseF+BPoZ2F9YpXYo/Z5IEAKOY9lIA8PRBEWrHYrH0pjjkatLUJydX/WspKIeq83BiDjRvNQVkiSRmhjCoeOncTjFcFc1iQTRzymN9ThP5qKPzxDfci9Ba44HpwO5pqzTx7Ssb5AyWPTrdFVaYihWm5OCsjNqh9KviQTRzzmK94Aio4vPUDsUj6Z1LUHa+X6I3OIaokJ9Gehr6KmwvNawuCB0Wsk1h5WgDpEg+jn7sd1IfiFoQuPUDsWjSf5mJKNfp5cgdThlCkprxd3DZTIZdKQMDmJvQVW3J0oULp9IEP2YYmvAWXYInWhe6pAkSWjCOr8EadHJszQ5ZNdCOELXjR1qprLWSmllvdqh9FudShB2u52srCwyMzPJzMxk6dKlyB3MlVJTU8PEiRMpLCxstf35558nOTm51U9WVtbl10C4bPZju0F2oE/MVDuUPkFrTkCuOYHcZO1w39ziGiQQq8d1w+ghZiQJsvMsaofSb+k6s9OyZcvYtm0br7/+Og0NDTz22GMEBASwYMGCdvc/ffo099xzD9XVbccxFxQUsHDhQm677TbXNh8fn8sMX+gOx9Hvmkcv/dC+LlyaNiwBFAXbqSIYEHvJfXOLa4iN8MfXpO+l6LxPgK+BodGBZOdbmH1Vgtrh9Esd3kHYbDbWr1/P4sWLSUtLY+LEiTzyyCOsXbu23bbBbdu2MXv2bOz29hcfLywsZMSIEZjNZtePn59f92sidIlcfxrnySPok64QzUud1JJIbScLLrmfze6k8MQZ0bzkBmOTzZysauBUdYPaofRLHSaI3NxcrFYrGRk/jnLJyMjAYrFQVtZ2yN/XX3/NHXfcwfLly9u81tDQwKlTp4iPF99Y1eYo/A5Q0CddoXYofYbGJwDJ34zt5NFL7nf0xBmcskJKrEgQ3TU2OQwQzUxq6TBBVFRU4Ofnh6+vr2ub2Wx2vXahJ598kvnz56PRtC26pT/izTffZOrUqcycOZPVq1d32J8huJ/96A40YQloBoarHUqfojXHd3gHkVdSg0aSGBItpvfuriB/I4mDAtidV6l2KP1Sh30QVqsVo7H1+gAGQ/O47qamrk2mVVhYiEajISoqipUrV3Lw4EGeffZZnE7nRfszLiYk5MdmKbPZv0vHerqero+t4jh11aWEXDuXgb303nnLNapNGMbpop0E+TjQ+bV/h1B4so4hMYHERvedOwhPvj5TM2JZ/dFBGmWICe9cnJ5cn8uhVn06TBAmk6lNImj53WQydelkN954I1OnTiUoqPkPJzk5mdraWtatW9flBFFdXY8sK5jN/lgs3rO4SG/Up3H7f0GrwxYxulfeO2+6Ro4BUQBYcg+gi2s7uaGtyUl+SQ3XjY/tM3X29OszInYgkgSffFvETZ3orPb0+nSVO+qj0UitvlR3+riOdoiIiKCurg6r9cehfRZLc3tgeHjXmickSXIlhxZJSUlYLBbxMEwvUew27AU70CWMRzKJwQFdpQ0dDJLmohP3/dj/IIa3ukugn5HhccF8d6hcfE70sg4TREpKCj4+PmRnZ7u27d69m7CwMKKiorp0spdeeonbb7+91bbDhw8THx8vRtL0Envhd2C3oh92tdqh9EmSzoghbDDOyvYTxJEf+h+SRP+DW00cEUHVmUaOnhBzM/WmDhOEyWRizpw5ZGVlkZOTw44dO3jxxRe56667AKitraWurnO3P9OmTWPPnj288sorlJSU8OGHH7Jq1Sruueee7tVC6DR77hY0QVFow5PUDqXPMg4agrOyqN0lSI+U1BAf6Y/J0KlHjIROGj00FINew45DbQfGCD2nU09SL1q0iMzMTObNm8fDDz/M7NmzmTt3LgAPPPAAf/7znzt1srS0NFasWMHnn3/OT3/6U1asWMGiRYuYNWvW5ddA6DSn5TiypQj9sKnijq0bTFFDmpcgrT3Vantjk4Pjp+rE/Es9wGTQMWaomV25FdgdYtRjb+nU1xyj0ciSJUtYsmRJm9fWrVvX7jHR0dHk5eW12T59+nSmT5/exTAFd2javxH0PuiHTlI7lD7NGJ0MgLPiKNqgH5tZj5Y19z+I6TV6xqSRkXx3qIKcfAuZw8Xw7N4gJuvrJ+SzFhxFO5vvHgxi+cvu0AdHgtEXubL1PGNHSmrRaiSSokT/Q08YFheEOdDElr0n1A6l3xAJop9oOvAZSBoMo65VO5Q+T5I0aMMScFa07qjOK6khTvQ/9BiNJHFV2iCOlNSKqTd6iUgQ/YDSWI89byu6pCvQ+Ir2cXfQhiUi15xA+WFm18YmB8dO1YnpNXrY5NRBaDUSW/aeVDuUfkEkiH6gad8n4LBjSL1B7VC8RvMoMMW1PkRB2RlkRcy/1NMG+hoYPdTMtgOnsDucaofj9USC8HLyuVqaDm5GlzQBbXDXnlsRLs61BGlF88R9R0pqRP9DL5maPoiGRgffHxbzM/U0kSC8XFPORyA7MWbMVjsUryIZfdEEDsL5Q0d1Xkkt8YMCMBq0Kkfm/YYNDiIq1JdNu0rFk9U9TCQILyafrcSeuwX9sCloAsLUDsfraMIScVYc5VxjU/PzD2J4a6+QJIlrx8VQZqnncHGN2uF4NZEgvJhtx3rQ6jCM+ZnaoXglXcQQsDVQnH9U9D/0sgkjIgjwNbBpZ6naoXg1kSC8lKNkL47iPRjG/BzNAPHNtidoI4YCUFt0EK1GIlH0P/QavU7D9DFRHCiq5kSVGPLaU0SC8EKKo4nGbW+jCYwUzz30IGlgOJLJH01VEQmDAjDqRf9Db5o6OgqDTsOn3xWrHYrXEgnCCzXt/S9KnQXjxDuQtOKhrZ4iSRKKORGz/YRYf1oF/gMMTEmPYsehCiprrR0fIHSZSBBexlldQtOej9ElXYEueoTa4Xi9KkMUZm0dIyLE3YMars+MRaOR+GSHuIvoCSJBeBFFdtD49Wokky+mibd3fIDQbbnWUABiJDENtRqC/I1cmRbJtgOnqD7TqHY4XkckCC/StPcT5OpijJN/JVaL6yU7y4040IKlsOOdhR4xM3MwAJ+Ivgi3EwnCSzirimnK/hBdwnj08Rlqh9MvnG1ooqSqkXrfaJzl+WqH02+FDDRxZWokW/edxCL6ItxKJAgvoDiaaPxqFZKPP6bJv1I7nH7jSEnzQ1r6yGTkqmLXxH1C75s1qXnZ4o++PaZ2KF5FJAgvYNv9PnJNGaar/lc0LfWiI8U1mAxagoeMAkV2zcsk9L4gfyPTx0ax/VA5pRWdWwJZ6JhIEH2co7wA+/5P0adMQRebqnY4/crh4hqGxgSijxgKGi3OU0fUDqlfmzlhMAa9ln98mqt2KF5DJIg+TLHbaPx6FZJ/CMYJt6odTr9SVWulssbKiLhgJL0RrTkBx0nxwaQm/wEGrhsXw/b9pyg6eVbtcLyCSBB9mG3nv1HOVmKaMg/J4KN2OP3KoeOnARgeHwyAdlAKsuW46IdQ2XXjYxnoZ+A/Xx8VM726QacShN1uJysri8zMTDIzM1m6dCmyLF/ymJqaGiZOnEhhYevhf+Xl5SxYsIDRo0czY8YMPvroo8uPvh9zlBdgP/QF+hEz0A1KUTucfufQ8RoC/QwMCmle31sbmfJDP0SBypH1bz5GHbdek8yRkloOFFWrHU6f16kEsWzZMrZt28brr7/OSy+9xAcffMAbb7xx0f1Pnz7NggULqK5ue4Huu+8+DAYD77zzDvPmzWPx4sXk5ORcfg36IcVpx7Z1DZJfCMbxc9QOp99xygq5x083Ny9JEvDDCnMaLc6Toh9CbddNiCMs0Id3vi5ElsVdRHd0mCBsNhvr169n8eLFpKWlMXHiRB555BHWrl3b7i3ctm3bmD17Nna7vc1ru3btIi8vj2effZakpCRuvfVWZs2axbp169xTm36iac/HyLUnMU2+C0lvUjucfqfoRC0NjQ5X8xLwYz+E6KhWnV6n4aYpCZywNLDjULna4fRpHSaI3NxcrFYrGRk/PnyVkZGBxWKhrKyszf5ff/01d9xxB8uXL2/z2p49e0hOTiYgIKBVWXv37r3c+Psd+WwlTXv/iy5xghi1pJK9+RYAhscFt9ou+iE8R0ZKGHER/rz/TZFYu7obOkwQFRUV+Pn54evr69pmNptdr13oySefZP78+Wg0bYuuqKggLKz1ymahoaHtliO0r2URIOOEW9QOpd/am28h2uzHQF9Dq+3aqOGgyGI0kwfQSBK/uDqJ02dtfJF9Qu1w+qwO54K2Wq0YjcZW2wyG5j+MpqamLp3sYmU5nU4cDgc6Xeenpg4J+fGBMLPZv0txeLqL1edc4R7qivcQPO1OAuNiezmq7vGWa9Roc3D42Gl+Ojm+TZ2U4NEc/8yEoTqf0HFTVIrw8njL9WlhNvtjNvvzxZ4TfPJdMbOnDcFvgKHjAz2UWtenw09kk8nUJhG0/G4yda3922QyUVtb26YsvV7fpeQAUF1djywrmM3+WCze8+TkxeqjyE7OfboaaWA4TfFX9ak6e9M12ltQhcMpMySy/TppIpOpK8hBGdt36utN1wda1+fnE+P445FK1v33MHOmJqoc2eVxx/XRaKRWX6o7fVxHO0RERFBXV4fV+mO7qsXS3AYbHh7epZNFRES4jj2/rAubnYS27HnfINeewph5C5JWr3Y4/db+omp8jFqGxLS/jKsueiTK2Urks5W9HJnQnpgwPzJHhLN5dyk1dTa1w+lzOkwQKSkp+Pj4kJ2d7dq2e/duwsLCiIqK6tLJ0tPTycvLo76+3rUtOzub0aNHd6mc/kZx2GjK/gBNeBK6weK9UouiKBworCJ9aBg6bft/OrrokQA4yg72ZmjCJdx4ZQJOWeH/bT+udih9TocJwmQyMWfOHLKyssjJyWHHjh28+OKL3HXXXQDU1tZSV9e525+MjAwSEhJ49NFHyc/PZ8OGDXz88cfceeed3auFl2s6+AXKuVqM43/hGncv9L6TVQ1Un7UxNuXid87SwAgkvxCcZYd6MTLhUsICfbgqfRDf7DtJRc05tcPpUzr1oNyiRYvIzMxk3rx5PPzww8yePZu5c+cC8MADD/DnP/+5cyfTaHj55Zex2WzMmTOHN954g+eee4709PTLr4GXU5qsNO37L9qYVHSRyWqH06/t/+HJ3IxhF28SlSQJXfRIHCcOo8hieKWn+NnEOLRaiQ/FdOBd0qmeYaPRyJIlS1iyZEmb1y72kFt0dDR5eXntbl+zZk0Xw+y/mg5tBlsDxoyb1A6l3ztQWE1MmB8hA30u2WmojR6J/cgWnBVHRVL3EAP9jEwfE82n35fw0yviGBTq2/FBgpisz5Mpdhv2/Z+hjUlFa45TO5x+7Vyjg4KyM6QmhnS4ry56JGh0OIr39EJkQmddnxmLwaAVdxFdIBKEB7Pnfoliq8c45mdqh9LvHTxWjVNWOpUgJIMP2qhhOI7vETOKehD/AQauyYhm15FKyirrOz5AEAnCUymOJpr2fYo2anjzRHCCqnLyLQQM0JM4aGCn9tcNHo1ytgL5zKkejkzoimvHxeJjFHcRnSUShIeyF2xHsZ7BkP5TtUPp9+wOmf2F1aQPMaPRdG4UmS62eeCF47iYZ8yT+PnomTE2hux8C2UWcRfREZEgPJCiyNgPfIYmdDDaQcPUDqffyy0+TWOTkzFDzZ0+RuMXjCZ0ME7RD+FxrhkXg9Gg5WPxXESHRILwQM6S/ci1pzCkXi+ee/AAOfkWTAYtwwYHdek43eDROCuOIlvF8peexM9Hz/Qx0ezKreRUdYPa4Xg0kSA8UNP+jUi+wegSxqkdSr8nywp7CqpITQxBr+van4subgyg4DiW3eG+Qu+6dnwMer2Gj7cXqx2KRxMJwsPYThXiPJWHYdQ1SJquTWAouN/RE2eoO2fvUvNSC01wDJqBETiKdvZAZEJ3BAwwMDU9iu8PV2CpFet3XIxIEB7mzK5PQGdEn9K3pov2VruOVEfzWOQAABzfSURBVKLTahiV0PHw1gtJkoQuMRPnySPI52o7PkDoVdeNj0WS4NOdJWqH4rFEgvAgsvUs9Ye/RT90MpJhgNrh9HtOWWZXbgVpSSH4GC/vbk6XOB5QcBTtcm9wQrcF+RuZODKCb/ef4kxD19a26S9EgvAg9tyvwelAP3K62qEIwJHiWs6eszNheNemtT+fNigKTXA0jkLRzOSJbpgwGIdDZvPuUrVD8UgiQXgIRXZgP/wlPgnpaAMHqR2OAHx3uBwfo7ZTT09fii5hPM6KAuT6ajdFJrhLRPAAxqaE8WVOGecaHWqH43FEgvAQjmPZKOdqGZgxU+1QBMDucJKTb2HMUDN6nbZbZemTJjSXWbDDHaEJbjZzQixWm5Ot+06qHYrHEQnCQ9gPfYHkb8YnSSwI5An2F1ZjtTmZMDyi22VpAsLQRiZjz/tGzM3kgeIiAhg2OIhNu0pwOGW1w/EoIkF4AOfpUpzl+RiGT0OSxCXxBNsPlhPgayBlcPtLi3aVPvkqlLMVOMvz3VKe4F43ZMZSW9/Ed4cq1A7Fo4hPIw9gP/QlaPXok69UOxQBqK23se9oNZNGRqDVuOdPRBefAXoT9iNb3VKe4F4j4oOJNvvx6c4SZHGX5yIShMqUpnPYC7ajS8xEMvmpHY4AbDtwCllRuCrNfYMFJL0RfeIEHEW7UJrEg1meRpIkbpgQy8mqBvYXisEELUSCUJk9fzs4bBiGT1M7FAGQFYWt+06SEhtIeLB7n0XRp1wJzibsR0VntScalxJGcICRT78XD861EAlCRYqiYD/8JRpzPNqwBLXDEYAjxTVYahu50o13Dy005gQ0oYOxH9wsOqs9kE6r4ZqMGPJLayk6KSZYBJEgVOU8dQS59iSGEeLBOE+xdd9JfE06MpK7PvdSRyRJwjDyWuTakzhPHHJ7+UL3XZU2CB+jjk+/F5P4QScThN1uJysri8zMTDIzM1m6dCmy3P5wsLq6Oh555BHGjh3LlVdeyZo1a1q9vm7dOpKTk1v93HPPPd2vSR9kP/QFGH3RJYxXOxQBOH22kew8C5NGRXb72YeL0SWOR/IJoOng5z1SvtA9PkYdU0cPIjvfQmXNObXDUV2nJphZtmwZ27Zt4/XXX6ehoYHHHnuMgIAAFixY0GbfJ598ksrKSt5++22OHz/OE088QVhYGD/5yU8AKCgo4KabbuLhhx92HWM0Gt1Unb5DbqjBcTwH/ajrkHQGtcMRgC+yy5AVhRljo3vsHJJWj37Y1TTlfIh8phzNwO4/ZyG41zUZMXy+q5TPdpZy53XJaoejqg7vIGw2G+vXr2fx4sWkpaUxceJEHnnkEdauXdumHfXEiRNs2rSJZ555hpSUFK6//nrmzp3LW2+95dqnsLCQ5ORkzGaz6ycgIMD9NfNw9tyvQVEwDL9a7VAEoLHJwdd7T5KRHEZooE+Pnks//GrQ6Gja/2mPnke4PIF+Rq4YEcG3B05xtp9P4tdhgsjNzcVqtZKRkeHalpGRgcVioaysrNW+e/fuJTAwkKSkpFb7Hjp0CLvdDjQniPj4eHfF3ycpTjv23K/QxqaiCQhTOxwB+Hb/Kaw2B9f+/+3deVRUV57A8W8VRVHILlKiaBBXVBRxI+6Na1xHu01aYzJmRm3tqJ1Mm8SI3ZzExG7aaGeTZERtY4wTO51xQ2O7NsaIMQJqFCKuuLMIIiUWtVBv/mCoABYW+6uC+znHA+++5f5uFdYP3vK7A9s3eF/KFr64dhuKKeM7LEX3G7w/oeaeiXyqtIhfyi37GzdhdhNEdnY2np6eeHh4WNsCAgKs6ypvq9VW/MALCAjAbDZz79498vLyuH//Pvv27WPUqFGMGTOGNWvWYDQ2ryxtvnoKSV+IuudouUMRKC3rfTD5Jp3b+dCprU+j9KnuMxEkC8az3zRKf0LNtPH3IKJrAP9KvUWxsfkW8bN7DUKv1z92jUCtLj1nXvmD3d62N2+WltT18vIiLi6OzMxMVq5cSWFhIW+//XaNAvf3//mhsoAArxrtK7fbe/6Fa8u2BEY8bbO0hrONpzoceUyHfrhObkEx86f1rnacdR5PgBc5vUZQlH4Uv1EzUHnWT0mPWofjwO9PbdTHeJ5/JpTXPjpGyuV8po7oVA9R1Z5c74/dBKHRaB5LBGXLGo2mRtsOHDiQ77//Hj+/0snfQ0NDAViyZAnLly+3JpPqyMt7iMUiERDgRW6urtr7ya0k5yqGO5dwGzyLe/cenzDd2cZTHY48JnOJha3/vECHQC9CtB7VirO+xiN1H4d0LpG7R75CM2hmnY9XW478/tRGfY2nZQtXurX3Zfu/LjGwa6saz0leX+pjPEqlosIv1dXez94GgYGB6HQ69PqfywPk5uYC0Lp168e2LVtXJicnB1dXV2tSKPtapnPnzpjNZvLz82scvDMynj8Arhpcuw6VOxSB0msP9x4UM3VYRxQKRaP2rfQJRNVlCKb0w2KuCAc1aXAH7usMHD9/V+5QZGE3QYSGhuLu7k5KSoq1LTk5Ga1WS1BQUIVt+/TpQ15eHteuXbO2paSkEBYWhlqtZtu2bYwaNarCMxTp6el4eno+du2iKbI8zMN85RSuoSNQqBv2ThnBPpO5hISkTDoFedOrY0tZYnDrPw1QYDi1XZb+hSfr0cGPkDbefHPiOiVVPPvVlNlNEBqNhunTp7NixQpSU1M5ceIEa9asYfbs2QAUFBSg05X++RMUFERUVBRLly4lPT2d/fv3s3HjRuu2Q4cOJT8/nxUrVpCZmcmRI0dYtWoV8+bNQ1lPVTMdmfH8IcCCOmyM3KEIwIFTN7mvM/BLGf56KKP09EcdNgbzpSRK8kQNIEejUCiYNDiYew+KOZne/EqBV+tBuddffx2DwcDcuXNxc3Nj+vTpzJkzB4DFixcTFBREbGwsALGxscTExDBz5kx8fHx49dVXGT9+PADt2rVjw4YNvPfee0ydOhVvb29mzpzZLJ6klox6TBcSUYUMQOnVSu5wmr37OgN7kq4T0aUV3TvI89dDGXWfiRgvHMXw/TbcJ7wuW7ISbAvv3Ip2AZ7sSbrO0z0CUSqbz/ujkJy0apizXaQ2njuA4cT/0GJqzBML8znLeGrCEce0PiGNUxdyeXdeJNoaPhjXEOMxph3CcPwLNKNfxrWRS6844vtTFw0xnuQLOXyy8zzzJvVgUFjjPv3u0BephbqTSswYf/wnLoFdRdVWB3DxZgEn0rJ5JrJ9jZNDQ3HtPhKl/1MYTmxDMhXLHY5QSd9uAbTXerLru2vNalpSkSAagenScaSifNQRk+UOpdkzmkrYtO8C/t4aJj7dQe5wrBRKJZohLyIV5WNI2SV3OEIlSoWCacM6klOgJ+l8ltzhNBqRIBqYZCnBeGYvylYdcGkXJnc4zd72b6+Snf+I/5wQipu6YSq21pZLYBdcuw3HdG4/JbmZcocjVBLe2Z+QNt4kHL+Gydw8/ooQCaKBma+cRCrMQd13srj4KLNLtwo4eOomURFBsl+Yrorb079G4e5N8dENSCXNt8SDI1IoFPxyeEfyCg0cbiY1mkSCaECSpQRj6m6UfkGogiPkDqdZe6g3Eb87DX8fDc9GyVs24UkUbh5ohs3Gkn8L4+kEucMRKukZ0pJeHf1JSMpE96jp15ATCaIBmTKOYXmQhXrAL23WXBIah0WS2LAnnQdFRn47NQyNulp3d8tGFRyBqvMgjKcTKMm+LHc4QiXPRXWi2Ghm9/FMuUNpcOJTq4FIZiPGlJ0otZ1QBfeVO5xmbd/31/nxSh4zRnUhpI1zzD2iGfoiCs+W6I+sQzKKmc0cSVCAJyPC25J4+jZ38x6vp9aUiATRQIznDyE9KsBt4LPi2oOMUjJy2X70KpE9WhMVEWR/BwehULfAPWo+0sM8io89PjmXIK+pwzri5urClv0ZTfq9EQmiAVgeFWA8nYBL+96o2obKHU6zdfVOIesT0ujY1pv/GB/qdInaJbAL6n5TMV/5HlP6YbnDEcrx9lAzPaoTF24UNOnbXkWCaACG7/8OJSY0g56XO5Rm625eER99fRZvDzWLf9Ubtatj3dJaXeqISbg8FY4h6UvMWZfkDkcoZ3h4Wzq38+HvRy432QvWIkHUM/OdnzBfPoE6fDxKXzEhvRxy7j/ivS9PA/Bfz4Xj7VH9eUYcjUKhxD3qNyi8/Ck+uFaUBXcgSoWC2eO6oTeY2XLgYpM81SQSRD2SzEYM332OwqsV6ohJcofTLJUmhzOYzBZemxFBG38P+zs5OIWbB+5jX0EyG9Hve19ctHYgQQGeTB0WQvKFHL471/TmjBAJoh4ZTn6FpeAummEvoVC52d9BqFc3cx7y5y9SKTaaeW1GBO20NS9O5qhcWgbhPnYxloK76A/GIZWY5A5J+H/jI4MJfcqX/zl4iaz8ppW8RYKoJ+abP2JKO4Rr2BhUoqRGo7tw/T6xW1NRKhW8+UI/ggOb1hzLAKqgHmhG/Aclt9MoPvwpkkU8ae0IlEoFcyf1QOWi4JMd59Abms77IhJEPbA8zKM4cQNKv3a4DXy21scxJO+wuVy53dZ+Rbv/XOVxntSHvWOW396QvMNmjJX71m1dUq0xlF9X/jjl+9FtXVJhuez78v/+dfo2q7el4uNiYEnvXIJaeVQZZ1n7zmNXK/RR1WtUtPvPFcb3pNer7Lg7j1194raHN33yWNvOY1c5u/1vAMRuTbG2lf8KsDe7LW6DZ2HOTKX4SDySpaRafVY+Ttk+tvqwtU35bcvHV/mYlfet6phVxWVv/5qqzhjrQ0tvDfP/rSd37j0ifncaFkvTuB4hEkQdSaZi9Ps/QDKb0Iz+LQpV7S+IGlN32Vyu3G5rP0tWRpXHeVIf9o5Zfntj6i6bMT7Wd1FexeUq+i2/rvxxKvRTlFdhuex7Y+ouHqbs4fNj2WzZn0Go6javqr/C66edNvurPJ7dxzMr9GErLgBLVkaF8VWl/HHtPWE70PTDY227j2fS8d63AFy8+cDaVv5r2ffqsDG4RT6H+eoPFB9cy77jl6v1VG/l45SP1db+5bcpv235+Cofs/K+VR2zqrjs7V9T1RljfQkL8ef5MV04eyWPbUcuNYmL1o5dc8DBSZYSio+sw5J/C/dnfo+Ln/M8iOXs7ph92Vw0jKwSPyY8HcyYjC0oFc7/H7K61OETwEWNIWkrC7zusPHhL+QOSQBG9m1HVv4jDiXfwtVFyfRfdHK652/KEwmilqQSM8VH/hvz9dO4DXkBVftecofULJglJYeLe7Jf3xt3hZHfeh1kwC/+jO5i80kOZdRho1G4e9Hh0DqWeO+lJC8MF//2cofV7M0Y1YUSi8S+kzcosUg8N7IzSidNEuIUUy1IJgP6gx9jvpaM26DnUfccLXdIzUJaZj6rCyfyjT6CcPUN3vRJINS16d1aWBOunSL5WDcOV0UJj3a9g/HC0SZxasOZKRUKXhjTldH92nHg1E3itjvvhWuRIGrIUpDFo50rKLnxI25D/x11r7Fyh9Tk3TS3ZJ1uJGu2ncEoqZjreYTZnsfwUoqpOQEyzQGsfjARF20nDN9uQv/P98UDdTJTKBTMHN2FmaO7cPZyHu9+nsz1LOeb97taCcJkMrFixQoiIyOJjIxk9erVWCy2Z1TS6XQsWbKEfv36MWzYMDZt2lSj9Y5KspRgPH+Ioh1vIT16gPuEJah7jJQ7rCZLkiA9M591upGsLpzENbOWZ6M6scxnF73UzWOylpoolFrgPvF13AbPouTOBYr+/iaGU/+LZNTLHVqzpVAoGNO/Pa/N6MOjYjMrNp9i2+FLTvXXRLWuQfz1r3/l+PHjxMfHU1RUxBtvvIG3tze/+c1vHtt2+fLl5OTksHXrVjIzM1m2bBlarZaJEydWa72jkSQLJdfPYkjZgSXvBi5BPdGM+E+Unv5yh9Yk5RTo+U4fxg+GTuRsO4OHohUT3VMZ5paBNnI9urPNY6rH2lAolKjDxqDq0BfDD//AeDoB4/lDqHtE4dpzlNzhNVuhwX6snBfJ14lXOHDqJsfP3WXsgPaM7NcOD42r3OE9kd0EYTAY+PLLL/nwww8JDw8HYMmSJaxZs4Z58+ZVuEJ/+/ZtDhw4wJ49e+jcuTOhoaFcvnyZzZs3M3HiRLvrHYmlMBfT1VOYM77F8iALhac/mtELUYX0d+q7EhxNsdHMRVMgGaY2XDC15dZ/nwD6EqLKYd6EHnQ7Ho2rQiSFmlB6+uM+cgElvZ7BeHYvxh/3YTy7j4VerTGmK1A91VvuEJudFhpX/v2ZUIaFtyXheCY7jl0jIek6EV1aEdmjNd2D/XB3c7x7huxG9NNPP6HX6+nfv7+1rX///uTm5nLr1i3at//5rokzZ87g6+tL586dK2z76aefYjKZ7K53dW38bCpJEhiKsDzIouT+bSw51yjJysBSUHrxU6nthGbkAlQdB6BQOmdFULlZLBI6i4b8LB25BXqy7z/i9r0irj+YTNb73yJJY1FiIUSVw7NRneiRuoqWLkV4hc1AlySSQ225BHTAffTC0l92Lh3H74fDGL7bjAH4g48X182tuGn259GVp/BX6pBKTChcHPs3WmcX0sab303vzY1sHcd+vMvJ9GxOXcjBRamgQxsvQgK9eaq1F4H+LdD6uuPVQt73w26CyM7OxtPTEw+Pn4ueBQQEWNeVTxDZ2dlotdoK+wcEBGA2m7l3757d9W3atKl24Eqlwub3T2K6lERJzhUkkwGMeiRDEZL+AZSva6PSoArsjCr8GVza90bp2XiT26t8Sl/XsvGofAJQKhXWrzXdr6ptq/N6lR2z8nHLL5e1PbS4kZiUid5gRq+ahBEV5l1pFBvNPOQl9LGHsChmw67z1v18Pd0I8VcSFd4F7U9/5ylVPm4KEx6DZlGU2QJoYR17+TjKL1fVbus1q7xe6+de5etbftnWmG0tl+2j9XO3ud7WtuWVj0fr527za9l2tvZ9Up9KXy2qAdOIT25N7MxOlNy9wIMfUujpks/Tyuvc27eOFcFQ/PV3KNw8UWg8QN2CV9rpMUgqzKgosqg5buj2WFzl+60cp62YbLXXZkz22Iqjqn5q20dddGjjTYc23jw/pguZd3VcuvWAK3cK+enmfX68+vMNBi4uCrw93HBXu9DCTUVosB/Dw9vWuL/ajlEh2bknbufOnaxatYqkpCRrm8VioXv37mzatInBgwdb2+Pi4khMTOQf//iHte3mzZuMHj3aemrpSeuDg4NrNQhBEASh/tm9i0mj0WA0VpwMo2xZo9HUaNuaHEsQBEGQl90EERgYiE6nQ6//+Xa53NxcAFq3bv3YtmXryuTk5ODq6oqfn5/d9YIgCILjsJsgQkNDcXd3JyUlxdqWnJyMVqslKKhi7aE+ffqQl5fHtWvXrG0pKSmEhYWhVqvtrhcEQRAcR7VOMU2fPp0VK1aQmprKiRMnWLNmDbNnzwagoKAAna70CcGgoCCioqJYunQp6enp7N+/n40bN1q3tbdeEARBcBx2L1JD6bMQ7777Lnv37sXNzY3p06fz+9//HoVCwYsvvkhQUBCxsbFAacKIiYnh6NGj+Pj4MGfOnAoJwN56QRAEwTFUK0EIgiAIzY8o1icIgiDYJBKEIAiCYJNIEIIgCIJNTSZBHD16lMmTJ9O7d2+mTZtGcnKy3CHVi9TUVLp3786VK1fkDqXWJEkiPj6e0aNHExERwYwZMzh9+rTcYdVITUreO4OsrCx+97vfERkZyZAhQ4iOjqawsFDusOpFTEwMzz33nNxh1JnZbOa9995j8ODB9O/fnyVLljT6e9QkEkRaWhoLFy5k8uTJJCQkMGDAAF5++WXu378vd2h1YjQaWb58uVN/EAFs2bKFjRs38uabb7Jjxw769u3LnDlzyMrKkju0aitf8v79999n586dbNiwQe6wasVisbBw4UKKiorYvHkzn376KRcuXGDZsmVyh1ZnJ0+e5KuvvpI7jHqxevVqdu/ezQcffMBnn33GxYsXWblyZeMGITUBr7zyirRo0SLrstlsliZOnCgdO3ZMxqjqbvXq1dLzzz8vde3aVbp8+bLc4dTalClTpHXr1lVoGzt2rLRlyxaZIqqZ4uJiKTw8XEpMTLS2bd++XRoyZIhksVhkjKx20tPTpa5du0o5OTnWtuTkZKlbt26STqeTMbK6efTokTRmzBhp5syZ0rPPPit3OHVSWFgo9ezZUzp8+LC1LTExUZo8eXKj/sw1ib8gTpw4wfjx463LLi4u7Nmzh6FDh8oYVd2kpaWxY8cO3nzzTblDqbM//OEPTJky5bF2ZzmlYa/kvbNp06YN69evt1ZlhtLZzyRJ4uHDhzJGVjcffPAB/fr1Y9CgQXKHUmfJycm4uroyfPhwa9uIESPYvXt3o85H4/QJQqfTUVBQgIuLCy+//DKDBg1i1qxZnDt3Tu7Qas1kMhEdHc0bb7zRJGpUDRgwgMDAQOtyYmIimZmZFSoBOzJ7Je+dja+vb4UPHoDPPvuMkJCQCu+TMzlz5gx79+5l6dKlcodSL27cuEFQUBCHDh1iypQpDB8+nJiYmEZP4I43hZENBoOhyvPVZZMMvfPOOyxevJhFixbx9ddfM3v2bPbt2/dYQUFH8KTx+Pj4sHXrVrRaLVOmTHGK31DtjcfX19e6fPnyZZYuXcqkSZPo06dPY4VYJ3q9Hjc3twptZbXDKlcndkbx8fEcPHiQ+Ph4uUOpFaPRSHR0NNHR0RV+1pxZUVERd+7cYf369SxfvhxJknj33XdZtmwZH3/8caPF4RQJ4ty5c8yaNcvmul//+tcATJ061fr9H//4R5KSkti1a5fNebPl9qTxjBs3jpMnT7Jjx45Gjqr2njSeRYsWsXjxYgDOnz/PvHnz6NKlS+NfbKuDplymPi4ujo8++oiYmBiGDRsmdzi1EhcXR3BwMBMmTJA7lHqjUqkoKipi1apVdOrUCSj9JXjGjBnk5uZWOD3YoHE0Si911L9/fzIyMmyuM5vNbN++3foiQun51JCQEG7fvt1YIdbIk8azdu1aDh06ZJ2jW/r/Sii/+tWvWLBgAQsWLGi0OKvrSeMpc+rUKebPn094eDiffPKJU32wli957+5eOkNZVSXvncmf/vQnPv/8c9566y1mzpwpdzi1lpCQQG5uLhEREUDpKdqSkhIiIiLYu3cvbdvWfAY2uWm1WpRKJR07drS2lX1/584dkSCqS6VS0atXL9LS0pg2bRpQehvf1atX6devn8zR1dwLL7zA5MmTrcvZ2dm8+OKLrF27lrCwMBkjq72LFy8yf/58Bg4cyEcffeR0pd3Ll7wvu/GhqpL3zmLt2rV88cUXxMbGMnXqVLnDqZMtW7ZgNpsrLCcnJ/Phhx8+NsWxs4iIiMBisZCenk7Pnj0BuHTpEkqlslF/5pw+QQDMnTuXV199le7du9O/f3+++OIL8vLybN454+h8fX0rnEd1cXEBSu88cdbzq9HR0Wi1WmJiYnjw4IG1vUWLFhUu/Dqq8iXvY2NjMRgMrFmzhpdeeknu0GolIyODuLg45s6dy5AhQypM4tWyZUvrz5yzqPyB6ePjg1qtduopjIODgxk7dizR0dG88847ALz99tuMGzeOVq1aNVocTSJBjBo1ipUrVxIXF8dbb71F9+7d+dvf/ua0vz00Jbdu3bLeURYVFVVh3bx583jttdfkCKvGXn/9dQwGA3PnzrWWvJ8zZ47cYdXKgQMHsFgsxMfHP3Zh+ptvvqlwulaQz1/+8hdiY2OZM2cOkiQxbtw4oqOjGzUGUe5bEARBsMnpn4MQBEEQGoZIEIIgCIJNIkEIgiAINokEIQiCINgkEoQgCIJgk0gQgiAIgk0iQQiCIAg2iQQhCIIg2CQShCAIgmDT/wGJC9ldOhn4CAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "# NOTE: we are storing y's as 1's or 0's but still doing the mean stuff correctly\n",
    "def get_sample(n, mean, var):\n",
    "    x_sample = []\n",
    "    y_sample = []\n",
    "    for i in range(n):\n",
    "        y = np.random.binomial(1, .5)\n",
    "        if (y == 0):\n",
    "            x = np.random.normal(-mean, var)\n",
    "        else:\n",
    "            x = np.random.normal(mean, var)\n",
    "        x_sample.append(x)\n",
    "        y_sample.append(y)\n",
    "    return x_sample, y_sample\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(font_scale=1.3)\n",
    "\n",
    "def split_by_class(sample):\n",
    "    y_one = [pair for pair in sample if pair[1] == 1]\n",
    "    y_neg_one = [pair for pair in sample if pair[1] == 0]\n",
    "    return y_one, y_neg_one\n",
    "\n",
    "X, y = get_sample(100, 2, 1)\n",
    "sample = list(zip(X, y))\n",
    "class_one, class_neg_one = split_by_class(sample)\n",
    "sns.distplot(list(zip(*class_one))[0], hist=False, rug=True)\n",
    "sns.distplot(list(zip(*class_neg_one))[0], hist=False, rug=True)\n",
    "plt.title(\"Distribution Plot with Mean = 2\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We copy below the math needed to solve for H(Y|X) below:\n",
    "    \n",
    "We will do the same math as above to calculate the true conditional entropy when $\\mu = 1$. The math is actually very clean! $H(Y)$ is still the same.\n",
    "\n",
    "$H(X) = -\\int_x (\\frac{1}{2}\\mathcal N(1, 1) + \\frac{1}{2}\\mathcal N(-1, 1))\\log (\\frac{1}{2}\\mathcal N(1, 1) + \\frac{1}{2}\\mathcal N(-1, 1)) dx$  \n",
    "\n",
    "This is messy... BUT notice that $\\mathcal N(1, 1)$ and $\\mathcal N(-1, 1)$ are just mirrors of each other! Besides the different means, the distributions are completely the same. So if we distribute the first summation we should get two identical \"entropies\"! (I did this numerically in wolfram alpha as well just to make sure).  \n",
    "\n",
    "$H(X) = -2(\\int_x \\frac{1}{2}\\mathcal N(1, 1) log (\\frac{1}{2}\\mathcal N(1, 1) + \\frac{1}{2}\\mathcal N(-1, 1))dx = 1.75577$  \n",
    "\n",
    "Now for $H(X|Y)$ . \n",
    "\n",
    "$H(X|Y) = \\frac{1}{2} H(X|Y = 1) + \\frac{1}{2} H(X| Y = -1)$ . \n",
    "\n",
    "Again we use the fact that $H(X|Y = 1) == H(X|Y = -1)$  \n",
    "\n",
    "So, $H(X|Y) = \\int_x \\mathcal N(1, 1)\\log \\mathcal N(1, 1) dx = 1.41894$  \n",
    "\n",
    "So $I(X, Y) = .33683$. $H(Y|X) = .356316$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, we need to approximate the parameters of the distributions of X and Y.\n",
    "For distribution of Y, Y is just Bernoulli (multinomial if k sample). So the MLE estimate of p hat is just n/N, where n is number of 1's. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_p(y):\n",
    "    return y.count(1)/len(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.492"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, y = get_sample(1000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "estimate_p(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For distribution of X (sum of 2 gaussians), we can get the parameters using Gaussian Mixture Model. We can change the covariance_type for the continuous case later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 0.60651298],\n",
       "        [-0.61055456]]), array([[[0.63213827]],\n",
       " \n",
       "        [[0.6256372 ]]]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import mixture\n",
    "X, y = get_sample(100000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "def estimate_x_params(X):\n",
    "    model = mixture.GaussianMixture(n_components = 2, covariance_type = \"full\")\n",
    "    model.fit(X)\n",
    "    return model.means_, model.covariances_\n",
    "    \n",
    "estimate_x_params(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now, for our plug in estimator, we use the following equation:  \n",
    "$\\theta(Y, X) = H(Y) - H(X) - H(X|Y)  = H(Y) - I(X, Y) = H(Y|X)$  \n",
    "So, our estimator is \n",
    "$\\theta(\\hat Y, \\hat X) = H(\\hat Y) - H(\\hat X) + H(\\hat X| \\hat Y)$  \n",
    "X and Y denote the true distributions. $\\hat X, \\hat Y$ denote the distributions with estimated parameters. Sorry for crappy notation right now.\n",
    "\n",
    "To solve this, we require numerial integration solver and then just follow the math from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.418938533204673, 2.8208079448643897e-10)\n",
      "(1.7557693535515035, 1.428584327545427e-08)\n"
     ]
    }
   ],
   "source": [
    "import scipy.integrate as integrate\n",
    "import math\n",
    "\n",
    "def normal_entropy_f(t, mean, var):\n",
    "    return -1/(2*math.pi*var)**(1/2)*np.exp(-1/(2*var)*(t - mean)**2)*np.log(1/(2*math.pi*var)**(1/2)*np.exp(-1/(2*var)*(t - mean)**2))\n",
    "\n",
    "def two_mixture_normals_entropy(t, mean_1, mean_2, var_1, var_2):\n",
    "    return .5*-1/(2*math.pi*var_1)**(1/2)*np.exp(-1/(2*var_1)*(t - mean_1)**2)*np.log(1/(2*math.pi*var_1)**(1/2)*np.exp(-1/(2*var_1)*(t - mean_1)**2)*.5 + 1/(2*math.pi*var_2)**(1/2)*np.exp(-1/(2*var_2)*(t - mean_2)**2)*.5) + \\\n",
    "            .5*-1/(2*math.pi*var_2)**(1/2)*np.exp(-1/(2*var_2)*(t - mean_2)**2)*np.log(1/(2*math.pi*var_1)**(1/2)*np.exp(-1/(2*var_1)*(t - mean_1)**2)*.5 + 1/(2*math.pi*var_2)**(1/2)*np.exp(-1/(2*var_2)*(t - mean_2)**2)*.5)\n",
    "\n",
    "print(integrate.quad(normal_entropy_f, -10, 10, args = (0, 1)))\n",
    "print(integrate.quad(two_mixture_normals_entropy, -10, 10, args = (1, -1, 1, 1)))\n",
    "#NOTE: this doesn't work for mean = 0\n",
    "def plugin_estimate_cat_1D(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    h_x_cond_y = normal_entropy(x_params[1][0])*.5 + normal_entropy(x_params[1][1])*.5\n",
    "    h_x = integrate.quad(two_mixture_normals_entropy, -20, 20, args = (x_params[0][0], x_params[0][1], x_params[1][0].item(), x_params[1][1].item()))[0]\n",
    "\n",
    "    return h_y - h_x + h_x_cond_y\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.418938533204673, 2.8208079448643897e-10)\n",
      "(1.7557693535515035, 1.428584327545427e-08)\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "def normal_entropy_f(t, mean, var):\n",
    "    return -norm.pdf(t, mean, var)*np.log(norm.pdf(t, mean, var))\n",
    "\n",
    "def two_mixture_normals_entropy_f(t, mean_1, mean_2, var_1, var_2):\n",
    "    return -.5*norm.pdf(t, mean_1, var_1)*np.log(.5*norm.pdf(t, mean_1, var_1) + .5*norm.pdf(t, mean_2, var_2)) - .5*norm.pdf(t, mean_1, var_1)*np.log(.5*norm.pdf(t, mean_1, var_1) + .5*norm.pdf(t, mean_2, var_2))\n",
    "\n",
    "print(integrate.quad(normal_entropy_f, -10, 10, args = (0, 1)))\n",
    "print(integrate.quad(two_mixture_normals_entropy, -10, 10, args = (1, -1, 1, 1)))\n",
    "\n",
    "def normal_entropy(var):\n",
    "    return .5*np.log(2*math.pi*math.e*var)\n",
    "\n",
    "#NOTE: this doesn't work for mean = 0\n",
    "def plugin_estimate_cat_1D(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    #h_x_cond_y = integrate.quad(normal_entropy_f, -20, 20, args = (x_params[0][0], x_params[1][0].item()))[0]*.5 + \\\n",
    "    #integrate.quad(normal_entropy_f, -20, 20, args = (x_params[0][1], x_params[1][1].item()))[0]*.5\n",
    "    h_x_cond_y = normal_entropy(x_params[1][0])*.5 + normal_entropy(x_params[1][1])*.5\n",
    "    h_x = integrate.quad(two_mixture_normals_entropy_f, -20, 20, args = (x_params[0][0], x_params[0][1], x_params[1][0].item(), x_params[1][1].item()))[0]\n",
    "    return h_y - h_x + h_x_cond_y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([[-1.00897678],\n",
      "       [ 1.05831844]]), array([[[0.94895086]],\n",
      "\n",
      "       [[0.90439195]]]))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.32401493]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, y = get_sample(10000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "plugin_estimate_cat_1D(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAERCAYAAAB4jRxOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XlclOXeP/DPAMmwiIhCkOJuqCgoLhwFtw4dy1M+5svME09pR0lzbdMUS9JzKk+LPsligZ5cstJy6fc7Zces1FSOuZai4YNmSgkhpCwOwzLz/EEzMczCLPfcy8zn/Xr1esU999xzzQV+7+u+ru91XSq9Xq8HERF5PB+pC0BEROJgwCci8hIM+EREXoIBn4jISzDgExF5CQZ8DyZVApYnJn658zt5Yn2RPDHgS+iRRx5BTEyM8b/Y2FgkJSXhySefxMWLF03OXbJkCaZMmWL3tdeuXYudO3faPKf5NYuLixETE4ODBw86/kWauXDhAh5//HHjzzt37kRMTAy0Wq1L121Ny7rs06cPEhISMHXqVJPv5Gg91tXVYcWKFThy5Igg5WxZH8ePH8eiRYuMr2dmZiIpKcmpa8bHx0Oj0Vg8Z9y4cYL8fqV05coVzJs3D8OHD0diYiLmzJmDq1evSl0sRfGTugDebsSIEVi4cCEAoL6+HiUlJXjnnXcwefJkvP/+++jTpw8AYM6cOQ4FzezsbLz44os2z3H0mvb497//jXPnzhl/HjNmDLZt24Y2bdoI+jmWNK9LvV6P6upqbNmyBbNnz8aHH36I2NhYh6/5yy+/4L333sPYsWMFKWPL+vjoo49w5coVQa6t1Wpx+PBhpKSkmBy/cOECLl++LMhnSKWmpgbTp09HaGgoVq5cCZVKhZycHDzyyCP417/+heDgYKmLqAgM+BILDQ3FwIEDTY6lpKRg0qRJyMjIwLZt2wAAXbp0Efyz3XHNlsLCwhAWFub2zwEs1+XQoUMxatQobNu2DStXrhSlHLa4sz4GDhyIffv2mQX8vXv34s4778SFCxfc8rli2Lt3L3755Rds27YN4eHhAID4+HiMHj0an3/+OR544AGJS6gM7NKRoYCAAMyYMQOnT582du207Ir46quv8MADDyAuLg7Jycn4+9//jtraWgBATEwMAODFF1/EI488AgC46667sGbNGjzwwAOIj4/Hxx9/bLF749KlS5g6dSoGDBiABx54APn5+cbXLHU3HDx4EDExMSguLkZmZiaysrJw/fp1xMTE4OjRo2ZdGDqdDlu3bsWf//xnxMXF4d5778WOHTuM1zN0LR04cADTpk1DXFwcxo4diw8++MCpulSr1ejWrRt+/vlni69XV1fjpZdewtixYxEXF4eHH34Yp06dMpblj3/8IwAgLS0NS5YsMXv/xo0bMWjQIDQ0NBiPTZw4EePGjTP+rNPpkJiYiA8//NCkPpYsWYJdu3bh22+/Ndahwc6dO5GSkmIsU1FRUavf9e6778ZXX32FxsZGk+N79+41KY/BxYsXMXPmTAwcOBCJiYn429/+ZvwbApqektavX4/x48ejf//+GDJkCObNm4fS0lLjOTExMfj4448xb948DBw4EElJScjKyrJZzuZdby3/s9YNGRYWhscee8wY7AEgPDwcwcHB+Omnn1qtG2rCgC9Tw4cPBwCcPn3a7LUff/wRCxYswLBhw5CXl4eFCxdix44dyMzMBADjU8HMmTORkZFhfN+GDRvwX//1X3j99dfxhz/8weLnvv7660hMTERWVhYiIyORlpaG//3f/7WrzA8++CAmT56M0NBQbNu2zWIXymuvvYZXXnkF48ePR05ODpKTk5Geno6tW7eanLd06VIMHz4cb7/9NhISEpCRkWFX0GupoaEBP/30Ezp16mT2WmNjIx577DHs2bMHc+fOxdq1axEYGIhHH30UBQUFiIiIMAavpUuXYs6cOWbXSE5Oxq1bt4zdWJWVlSgsLMTly5dRUVEBADh37hxu3LiBkSNHmrx3zpw5GD16NHr37o1t27YhIiICAFBRUYG8vDw8++yzeP311/Hzzz/j6aefbvW7pqSk4ObNmzh+/Ljx2JUrV3Dp0iXcddddJueWlpYiNTUVNTU1WL16NRYvXoxPPvkEixcvNp6Tl5eHrKwspKam4p///Ceefvpp/Oc//8Grr75qcq2///3v6NKlC9atW4fx48cjMzMTBw4csFrObdu2Wf1vzJgxFt8zevRoPPPMMybHTp06hZs3b6JHjx6t1g01YZeOTHXo0AEAUF5ebvba2bNnUVdXhxkzZiAiIgKJiYlo06aNsZVp6Nbo3LkzevXqZXxf//79MX36dJufO2nSJDz11FMAgKSkJKSkpGDLli12dYdERkYiMjISfn5+Zl0rQFMg27JlC+bOnYsnnngCQFPArKmpwdq1a/HQQw+ZlGP27NkAmh7dP/vsMxw8eNDk+7Sk1+uNdaDT6VBSUoK33noL5eXlmDx5stn5X331Fb777ju89957GDx4MABg5MiRGD9+PLKzs5GTk4O+ffsCAHr06GGxC6xXr16IjIzEN998g7i4OJw4cQJRUVGorKzEyZMnkZKSgiNHjqB3796IjIw0eW+XLl0QFhaGGzdumNSXTqfDm2++iTvvvBMAUFZWhpUrV6Kqqgpt27a1+v0jIyMxYMAA7Nu3D4mJiQCaxlQSExMREhJicu6mTZugUqmQl5dn7P++4447MH36dJw/fx59+/bFL7/8ggULFiA1NRUAMGzYMFy6dAlffPGFybWSk5ONN4o//OEPxt/V6NGjLZbT0t+Go2pqapCRkYHo6GizLiyyjgFfgeLi4tCmTRtMmTIF9913H+666y7cf//98PGx/cDWvXv3Vq999913G//fz88PycnJOHnypMtlBoDvvvsO9fX1uOeee0yOjx8/Hjt27MClS5cQGBgIoCnIGwQGBiIkJAS3bt2yef09e/Zgz549JsfCwsKwcuVKDBgwwOz8EydOIDw83BjsAcDX1xfjxo0zPiXZIykpCceOHcPMmTNx7NgxDB06FCUlJSYBv2Xr3pagoCBjsAeabtwAWg34QNPv74MPPsCyZcsAAJ9//jkmTZpkdt6xY8cwePBgqNVq403S8PPRo0fRt29fPP/88wCA69ev4+LFiygqKsLJkydRX19vcq3mvyuVSoXbb7/d5u+qefdXS76+vlCpVDa/Y01NDWbPno2rV69i8+bNoiQEeAoGfJkqKysDAONjfnPR0dF455138NZbb2Hjxo3Iy8tDp06dkJGRYbVVBfz+1GBLy3Pat2+PqqoqB0tv2c2bNwEAHTt2tPiZ1dXVxoDv7+9vco6Pj0+r+erJycl48sknjeeHhISgc+fOVgNIZWWlWVmApptETU2NHd/o989dvnw5GhsbcezYMUyZMgXXrl3DkSNHUFtbixMnTpikqrYmICDA5GdD+XU6XavvTUlJwRtvvIHz58+jffv2OHv2LHJyckz65gHgxo0b+O677yx2u/3yyy8AgKKiIixbtgynT59GUFAQYmNj4e/vb/Z7cPR3ZStb6pVXXrF4gzKoqKhAWloaLl26hJycHIs3crKOAV+mvvnmGwDAoEGDLL4+ZMgQrF+/HjU1NTh06BDeeustPPXUU/jPf/7jUounsrLS5Ofy8nJjVolKpTIbEGyt1d1cu3btADS1GJu3VK9fv27yurNCQkIcCgAhISHGz26uvLzcobKMGDECNTU1OHnyJM6dO4eEhASUlpYiLy8PR44cgZ+fH4YMGWL39VzRo0cP9OrVC1988QVCQ0ORkJCAjh07mgwIA0Dbtm0xbtw4pKWlmV2jY8eO0Ol0eOKJJxAZGYk9e/age/fuUKlUeO2111xOI/3oo4+svmZ4mrGktLQUjz76KCoqKrBhwwYkJCS4VA5vxEFbGaqrq8M777yDoUOHomvXrmav7969G3/84x9RX1+PoKAg4z/cmpoaVFdXA0Cr3TvWNJ9gVFdXh4MHDxq7PAIDA1FZWWkS5E+cOGHyflufGxcXh9tuuw2fffaZyfE9e/agffv26Natm1NldtbgwYNRVlZm8h0aGxuxd+9eYz+zr69vq9cJDQ1F//79sX79eoSEhKBnz54YOHAgdDod3n77bQwbNszqTdjZ35MtKSkp+PLLL7Fv3z786U9/snhOQkICLl26hH79+mHAgAEYMGAAbr/9dqxevdo44HzlyhX85S9/QY8ePaBSqaDT6ZCfn2/Xk4Yths+z9F/79u0tvker1WLmzJm4ceMGtmzZwmDvJLbwJXbjxg1jJk5DQwOKi4vx7rvvoqSkBGvWrLH4nsGDB+P69et4+umnMXXqVGg0Gqxbtw4JCQnG1nhISAi++eYbDBo0yDh5yx7vvfcewsLC0KNHD2zatAkajQaPPfYYgKaui3/84x944YUX8OCDD+LUqVP49NNPTd4fEhKCmzdvYv/+/WZPJ2FhYXj44YeRnZ0NnU6HgQMH4uDBg9i5cyfS09PtCq5CGjNmDPr374+FCxfiqaeeQseOHbF161b89NNPeO211wDA+CTy9ddfo3PnzlYzQpKTk5GTk2McAwkMDES/fv1w+vRpvPDCC1bLEBISgitXriA/P9/q05yj/vSnP+Htt9+Gn58fVq1aZfGcadOmYdeuXZg7dy6mTp2KhoYGY0ptnz59EBoaiqioKGzYsAFBQUHQ6XR4//33ce7cOajVakHK6YiNGzfiwoULeOqpp1BbW2uSvXb77bcjKipK9DIpEVv4Ejty5AgeeughPPTQQ3jssceQlZWFbt264cMPP0Tv3r0tvic6Ohrr1q1DaWkp5s2bh6VLl6J3795Yu3at8ZwnnngC+/fvN5m2b49ly5bh448/xpw5c1BdXY1NmzYZUxrvvPNOrFy5EqdOncLjjz+O48ePGwOjwfjx49GrVy/MmzcPhw4dMrv+kiVLMHfuXHz00UeYPXs2jhw5gpdffhmPPvqoQ+UUgp+fHzZs2IBRo0bh1VdfxcKFC6HRaLB582bExcUBAIKDg/HXv/4V27Ztwz/+8Q+r10pOTgYAkwFgQzeOrQHbKVOmIDg4GI8//rjJDGVXxMbG4o477kDfvn3NMoMMoqOj8e6776Kurg5PPvkk0tPTERUVhc2bN6N9+/ZQqVRYu3YtfHx8sGDBAixfvhzBwcFYvXo1NBoNCgsLBSmrvQyZQWvWrDH+ezH81zKll6xTcYtDIiLvwBY+EZGXYMAnIvISDPhERF6CAZ+IyEsw4BMReQkGfCIiLyGbiVe//loDna71DNEOHYJRXl4tQomUifVjHevGOtaNbXKsHx8fFdq3D3LoPbIJ+Dqd3q6AbziXrGP9WMe6sY51Y5sn1A+7dIiIvAQDPhGRl2DAJyLyEgz4REReQjaDtkTeJL+gBDsPXER5pRYdQvwxaXRPDI+1vLIlkVAY8IlEll9Qgk17vkddQ9NGIuWVWmza8z0AMOiTW7FLh0hkOw9cNAZ7g7oGHXYeuChRichbMOATiay8UuvQcSKhMOATiaxDiL9Dx4mEwoBPJLJJo3uijZ/pP702fj6YNLqnRCUib8FBWyKRGQZmmaVDYmPAJ5LA8NhIBngSnSBdOiUlJViwYAESExORlJSE9PR0VFZWCnFpIiISiMsBX6fTYe7cuaipqcGmTZuwbt06fP/991i6dKkQ5SMiIoG43KVTWFiIs2fP4tChQwgPDwcALFu2DKmpqaiurkZwcLDLhSQiIte53MKPiopCXl6eMdgDgEqlgl6vR3W1vDYMICLyZi4H/NDQUIwaNcrk2MaNG9G9e3dERnJQiohILlR6vV7QbVxyc3OxZs0a5ObmYuTIkUJeWhw/bAW+XQbcugIEdgHiXwK6p0pdKiIilwmalpmdnY21a9di+fLlDgf78vJqu7YQCw9vi7KyKmeLaJP/te1oe24+VDpN04FbP0J/NA1VVbXQRk1xy2cKzZ31o3SsG+tYN7bJsX58fFTo0MGxMVLBZtq+/PLLyMzMxIsvvojUVGW2iIOKVvwe7H+j0mkQVLRCohIREQlHkBZ+VlYW3n33XaxatQoTJ04U4pKS8Kktdug4EZGSCJKWmZ2djZkzZyIpKQllZWXG18LCwuDr6+vqR4hGp+4M39qrFo8TESmdywF/79690Ol0yM3NRW5urslrn376KXr2VM6CUDW9Mkz78AHofQJQ0ytDwlIREQnD5YA/f/58zJ8/X4iySM4wMBtUtAI+tcXQqTujpleGYgZsiYhs4eJpLWijpjDAE9mB+/IqDwM+ETmM+/IqEzdAISKHcV9eZWLAJyKHcV9eZWLAJyKHcV9eZWLAJyKHcV9eZeKgLRE5jPvyKhMDPhE5hfvyKg+7dIiIvAQDvoz5X9uOsK9j0fHzdgj7Ohb+17ZLXSQiUjB26chUy7X5fWuvou25piUsOBOYxMCZtJ6HLXyZ4tr8JCXDTFpDXr1hJm1+QYnEJSNXMODLFNfmJylxJq1nYsCXKWtr8HNtfhIDZ9J6JgZ8marplQG9T4DJMa7NT2LhTFrPxIAvU9qoKajql4lGdTT0UKFRHY2qfpkcsCVRcCatZ2KWjov8r21324YpXJufpMKZtJ6JAd8FSkmddOdNiTwXZ9J6HnbpuEAJqZOGm5Jv7VWooDfelDiJi8j7sIXvAjFTJ51tpdu6KbGVLx1OaiIpMOC7QKfuDN/aqxaPC8mVriPm88sPtwckqbBLxwVipU660nXEfH754aQmkgoDvgvESp10pZXOfH754aQmkgoDvhOar2IZVLQC2o7joFN3hk9tMYKKVgg+IOpKK535/PLDSU0kFfbhO8hSf3pA8XqofnvdHamZNb0yTD4TcKyVznx+eZk0uqdJHz7ASU0kDrbwHWSxP73FOUKnZrKV7lmGx0Zi2r19jC36DiH+mHZvHw7Yktuxhe8ge7NbhM6CYSvds3BSE0mBLXwH2ZvdYus87mRFRFJgwHeQxayXFufY6l9X+sxXSzcr3sCIlIFdOg4ydKs0n/Wq7TgO/tf/bdcsWKXNfG0+w1fv1x6qxmqo9HUAfhugLpgD6PVQof73YzJcT8gRnAVLnooB3wmW+tNr7Hyvkma+tsxIUjVUmJ1jCP4mx2R8A2sNZ8GSJxO0S0er1eLPf/4zDh48KORlPYqSZr5aehqxlxxvYPbgLFjyZIIFfI1Gg4ULF6KoqEioS3okJc18dSVoy/EGZg/OgiVPJkjAP3v2LCZPnoySEu5o3xo55NTbO8hqT9DWq9pAj9tMj8n0BmYPzoIlTyZIwM/Pz8fo0aPxwQcfCHE5j6eNmoKKkQW4fvdNVIwsED3Y25slVNMrwzyYwwe628J+v1nF5qCq/zqPmRTGrf3IkwkyaJuWlibEZUgEDmcJqVSmeacqP1THvGp2rlIDfEvc2o88mUqv17dMI3dJTEwM8vLyMGrUKCEvK38/bAW+XQbcugIEdgHiXwK6p0pdKnPv+cB85gAAqICHTQcrsbsbcOtH81MDuwITLwtfNiJyK9mkZZaXV0Ona/3eEx7eFmVlVSKUyH4t0xdx60foj6ahqqpW9JZva/UTZmXTlkZ1Z1S0eF/HW1fM1gkCAP2tK7gus9+BPeT4tyMXrBvb5Fg/Pj4qdOgQ7NB7ZBPwlUxJk6kcWXlTrB29yBQnfpG7MOALQEmTqSzNFLY2M9jVZZnJcZz4ZRlvgsJgwBeA0lrC9q686cjNgYRha+KXlAFOyoDLm6BwGPAF4MktYS7LLC45TvySOuDK9SaoRIIH/MLCQqEvKXvOtoSbL0zG1jMBTRO8LAV3eyd+uaMlLnXAleNNUKnYwheIoy1hS1slKn2VSbKPraDsyvaH7mqJSx1wW7sJsn/fflwPXyK2MnuEwDXq5ckQlA0BzBCU8wualiVxZftDdy38JvVyE7ZmP7dWn2SKLXyJuDOzh08P8mVP94iz2x+6qyUu9abrtmY/L8o5zP59BzDgS8RtmT0/bEXbs7OgQqPJYbnOC/A27uwecbX/3xo5LDdh7SYodXeT0jDgS8QdmT3+17YD5xeYBXsDOc4L8DbuCsqAe1vict103Z316YkY8CXijhz3oKIVQOMtq6/LdV6AN3FHUG4+aBkc4Ifb/FSoqW30igFMV+vT2wZ8GfAlJHSOu60WvKfMC1A6obtHWmbmVGsa0MbPB2n39/PowGXgSn1KPb9ACooP+Mxl/521cQE9fBW9Rr2nEbJ7ROoceTlwtj69se4UHfCZjWKqplcGQs4vMOnW0fsEMNh7MDkNWiqte0ROdScWRefhuzuXXWm0UVOAYbkes/sUtU7qHHkDJebDy6XuxKToFr6SVqkUTfdUVARPkLoUXkmKFq7UOfIGSuwekUvdiUnRLXxrWSdCZaNwtqop1od1UrVwXZmZKyQldo/Ipe7EpOgWvjtXqeT4gCnWh21StnBdGQQ2PJVUVGoR5sJTiVLz4eU6v8BdFN3C10ZNQVW/TLf0WXN8wBTrwzYltnCbP5Xo4dpTia31bkg+FN3CB9y3XjvHB0yxPmxTYgtXyKcSOSy/QK1TfMB3F6XtYuVurA/blDgAKPRTibd1jygRA74VnryLlTNYH7YpsYUr96cSpeX1KwEDvhXcz9WU4XsHFy6Gqr4CQFPAp9852sKVOqDJ+anEG5c9EAMDvg3cz9WcqlEDleH/GyqYqeMkOQS05k8lrWXpWLo5Gd7rjhuWEvP6lYABn+xmK1OHAd8xcglohqeS8PC2KCurAmAe3ON6dsDhMyUmN6d//uscVD4qNDTqjceEvGEpMetJCRSdlukJnJnMJNUEKGbqCEeuAc3SBLKvTv1sdnNq1MMY7A2E2E7RwBuXPRADA76EDJOZfGuvQgW9cTKTrQDuzHuE4u6ZzUrlzA1YrgHN0pOHI4S6YTGv3z0Y8CXkzGQmKSdA1fTKMBuo9fZMHWdvwHINaK4GbKFuWN647IEY2IcvIWe6SKTsVmHmkjlnxzXkmsZpLVWzJV8VTPrwAeFvWMzrFx4DvoScmcwk9QQoZi6ZcuUG3FpAc2faZstrT78vFrFdQq2maiYNiMR3F8tFy9Ih92DAl5Azk5k4AUpe3HUDdjRt05Gbg6VrZ334LR69J8bhJw8GeGVhwJeQM10k7FaRF3fdgB1J23T05mDp2tr6RuO12ZXiuRjwJeZMFwm7VeTDXTdgR9I2Hc3pl2tKKLkfAz6JxlM3nHfHDdiRdW4cDeDOrqEj9VIQ5DoGfDfw1MDmCm6g4hhH1rmxlVmz5d/fmw22Wrq2/22+NjNs5LAUBLmOefgCk3JilJzZO3+A2yg2cSQP3Vag/urUz2bbLgIwu/a8B+NtBm5b3UakHIK08Ovr6/HKK6/gk08+AQA8+OCDePrpp+Hj4z33k99b9VeNi4sZcL0Z+9IX+RRgyt7B0+Gxkcj7/+fsuqYhSL82J8nk2s3X0rGE/f6eQZCIvHr1ahw+fBi5ublYs2YNdu/ejfXr1wtxaUUwbdVb5q6JUS1bxPhhq1s+x1X2LMvAbRSd58gMV2eCtFyXgiDHuBzwtVot3n//faSnpyM+Ph4jRozAM888g82bN0Ov17d+AQ9gKVC15Ghetj1dG5a6j/DN47LsBrFnWQYuzuY8S0s1WONMkJbrUhDkGJe7dM6fPw+NRoMhQ4YYjw0ZMgRlZWUoLi5GdHS0qx8hudYGYVsLSI7mZdvbtWHxRtN4S5bdR/akL0o9i1jJLE2YarmsMeB8kJbrUhBKJVXGk8sBv7S0FMHBwQgKCjIeCw8PN76m9IBvT/C1Fqj0AHTqaIezdOxdn8V6i9i8LC2/kxRZRK2lL3IWsWss9fn36hwqWGDhhCxhSJnx5HLA12g08Pc3fURs06YNAKCurs7u63ToEGz3ueHhbe0+12WH/wZYCL4hl/4GxM1oOjDoFeCbx4HGW7+f5BsI1bBc+HZPRYijn2klkPvWFpt+98AuwK0fzc5TQYXw6v8HdE81v8gPW4HzC4xl9a29ipDzC4C2asvniyl8RlM5vl0G3LoCBHaBKv4lhAhYLlH/dkS0/8RVbN5zHtd/1aBj+wA8em9fjBkcjQlj2mLCmN52XcOeurH2Od5AqL+d3YfyLWY87T70g92/K2e5HPDVarVZYDf8rFar7b5OeXk1dLrW+/xbyyYQWsdbVywOxOpvXcF1QzmCJ8C/b615qzl4AuBEWTvc1h4+v+0b21yjujMqml3Pv8cLaHs2DSq0rDc9Gk8tRUXwBLNrhJ1aCt/mNyYAaLxl9XzRBU8AklqUQ6Dft9h/O2Jp2WIs+1WDzO2nUVlVa3eL0Z66EeJzlErIv52yXy2P95X9qnHoM3x8VA41lAEBAn5kZCSqqqqg0WgQENA0KFdWVgYAuP322129vOTs7VcWaral/7XtUNWb/9L1qjZWujYs3yQdHQDlwKi4hOzDFWu7RLlsy6h0zs50FoLLWTp9+vRBQEAATpw4YTx2/PhxREREoFOnTq5eXnJib/oRVLQCKtSbHdf7BpvcUIxjC1au4+juVBwYFY+lbQQ37fke+QUlTl1PrBx55uILQ8qMJ5cDvlqtxuTJk7Fy5UqcPHkS+fn5eOONNzBt2jQhyic5bdQUVPXLRKM6Gnqo0KiORlW/TLcNclpraasafjX52VYqqK0bkifuWqW02blCz1oVK0eeufjCkHI3L0Fm2i5atAharRYzZ86Ev78/Jk+ejBkzZghxaVkQc3VKe7uQbHXB2LohedryykqcnSt0S9mRdXdcIdbneAOpMp5UepnMjpLroK3YWgYwANBDBUBvkuIZ9nWsxRsDAruiLOmMeAWWmLV6aFRHo2JkgckxufztLHjzIKo1DWbHO4T447U5SU5d09UxAXvrxltXzJTL305zkgzakrBMW+BX0ZRk2XQjbN56tZazrop/SfQyS0lpg9D5BSXQ1JoHez9flUstZbFajMzFVzYGfBkydCFZar0aJmAZWq8tu2ZCuqcKlsaoBEqZndu8ZWyJ/20+DKTkdgz4MtZa65U7Xyljdm7L/HVLamobRSwReSvvWb9YgZhC2Tqxs6icYSkrpyVmupAYGPBtkDrdz2IKJQBVQ41byyL193aUNmoKKkYW4PrdN1ExskBWwR5oPfuGmS4kFnbpWCGHdD89fHxYAAAPBUlEQVTD5wQXLoaqvgIqoOm/hgq3lUUO39vT2NqC0EdlmoOvpH58b83YUTK28K2Qy2Yc2qgp0PsGWd1FS2hy+d6exNLMSj9fFXxVgCET2dXZtmITerYwiYMB3wo5pfuJWRY5fW9PYWlmpf9tPmhsMe1ESXvEco9bZWKXjhVySvcTsyxy+t6epGX++l9XfWnxPKWsS8N1dZSJLXwr5LTmjJhlkdP39mRKX5fGmfLnF5RgUc5h/HXVl1iUc5jdPxJgwLdCTul+YpdF7xsAPX7bscsvTHZpjp5A6XvEOlp+9vnLA7t0bJDTxCYxymJpHZ+Wu32RMJS+R6yj5eda+vLAgO8hDPvUorYYYU6ugGnvXrotP9MTVt2UgtLXpXGk/Ozzlwd26diglAlIhpZ502Cr3pg772h5HcnQaf6ZKhc+k7yD0scsPAUDvhVKCmhC5c47spQD8/XJEUofs/AUDPhWKCmgCZU770iGDvP1yREt5yI0n2HMgVvxMOBboaSAJtQia45kA3FhN3LU8NhIY0tfqTOMlY4B3wolBTQhc+ftXYiM+frkDM7QlRYDvhVKCmjNW+YQac6AnOYpkHIwW0daTMu0QmmbfRvy9MPD26JCpB2v5DRPgZTB2sqhzNYRBwO+DQxo5C32n7iKjf8qcPsksEmje5rt/sVsHfEw4CsIJzqRO+QXlGDzZ4XQ1jdts2gYSAWEX59f6TOMlY4BXyG4MQm5y84DF43B3sCdyx4ofYaxknHQViHsmRfgf207sLub7GcGk7xwINV7sIWvEK3NCzA8AUCngQp8AvBkQm8tKNeBVG6hKDwGfIVobWMSRxc+I2GJFZwMywwbBj2F6G+fNLqnSR8+IP1Aqju+J7FLRzFamxegpJnBnkbMtd7dMXFpeGwk5j0Yb7IF47R7+0gaWDlByz3YwleI1uYFWHsCAPQI+zqWGT1uJOZa7+7qbx8zOBqxXUJduoaQOK7gHgz4CmJrXkBNrwzzzUsA9ueLQMzgJNf+dqF5y/cUG7t0PIRhqQMEdoXewutyXenTE4i51ru3LDPsLd9TbAz4HkQbNQWYeBlN7Xpz7M93DzGDU8tlhuXQ3+4O3vI9xSZol052djaKioqwZs0aIS9LDmoto4eEJfbsUW+ZuOQt31NMggX8nTt3IisrC/fcc49QlyQnWerPl+tKn56Cwclz5ReUYPehfJT9qlH8fACXA35tbS1efPFF7NmzB127dhWiTOQipa30SSRXnjYfwOWAX15ejp9++gkfffQRNmzYAK2WaVNywJU+iVwnZsqtGFwO+J06dcKWLVuEKAu5GVfbJHKMp80HaDXga7ValJRYnjHYrl07hIYKM1mjQ4dgu88ND28ryGd6Kov188NW4PwCoPEWgKbc/JDzC4C2aqB7qsgllI5Yfzv7T1zF5j3ncf1XDTq2D8Cj9/bFmMHRony2s+Xhvytz4e0DUParxuJxJdZXqwH/zJkzSE21HBDmzZuH+fPnC1KQ8vJq6HSWMshNhYe3RZlIOzopkbX6CTu1FL6/BXujxltoPLUUFcETRCqdtMT622nZ71v2qwaZ20+jsqpWkm4Ae8rDf1eWTUzubnHDlonJ3SWvLx8flUMNZcCOgD9kyBAUFhY6XSiSB661Ix659fvKrTxKYqif3Yd+YJYOKQdz88Ujt35fuZVHaYbHRmLCmN6St+iFwJm2XqK11TZJOGIutWAPuZWHpMOA7yUMa+00qqOhhwqN6mhU9ctklo4bCLnUQn5BCRblHMZfV32JRTmHnVpymevSkIGgXTqrVq0S8nIkMObmi0OopRaEmvTDjcPJgH34RG4gxFILQg62cukHAtilQyRL+QUlHGwlwTHgE8mMoSvHGg62krMY8EkU/te2I+zrWHT8vB3Cvo6F/7XtUhdJtix15RhwsJVcwT58cjv/a9tNlmvmlou22eqy4SYg5Aq28MntgopWmO+1yy0XrbKVN89gT65gwCe347IOjmHePLkLu3TI7bisg2OYN0/uwoBPbsctFx3HvHlyB3bpkNtxWQcieWALn0TBZR2IpMeATyQT+QUl7Lcnt2KXjgfiJCflMcyuNeTgGxZKc2Z1TCJrGPA9zQ9b0fbcfPjWXoUKeuMkJwZ9ebO1UBqRUBjwPc23yzjJSYG4UBqJgQHf09y6YvEwJznJG3elIjEw4HuawC4WD3OSk7xxdi2JgQHf08S/xL1rFWh4bCSm3dvH2KLvEOLPhdJIcEzL9DTdU1FVVYugohXwqS2GTt0ZNb0ymAOvAJxdS+7GgO+BOMmJiCxhlw4RkZdgwCci8hIM+EREXoIBn4jISzDgExF5CQZ8IiIvwYBPROQlGPCJiLwEAz4RkZdgwPdS3CSFyPu4HPD1ej1yc3ORkpKCQYMGYerUqTh16pQQZSM38b+2nZukEHkhlwP+li1bsGHDBixZsgS7du1CQkICZsyYgZISbs0mV0FFK7hJCpEXcjng79ixAzNmzEBKSgq6deuGxYsXIzw8HPv27ROifOQG1jZD4SYpRJ7N5dUyn3/+eURHR5sdr6ysdPXS5CY6dWf41l61eJyIPJfLLfyhQ4ciMvL3Nbz379+Py5cvY8SIEa5emtykplcGN0kh8kKttvC1Wq3V/vh27dohNDTU+HNRURGee+453HfffRg4cKBwpSRBGdbK5yYpRN5Fpdfr9bZOOH78OFJTUy2+Nm/ePMyfPx8AcPbsWaSlpaFnz55Yv3491Gq18KUlIiKntRrw7XHs2DHMmjUL8fHxyMnJQUBAQOtvaqG8vBo6XetFCQ9vi7KyKmeK6RVYP9axbqxj3dgmx/rx8VGhQ4dgx97j6odeuHABs2bNwrBhw/D22287FeyJiMj9XM7SSU9PR0REBJYvX46bN28ajwcGBiIoKMjVyxMRkUBcCvjFxcU4c+YMAGDs2LEmr6WlpeHZZ5915fJERCQglwJ+586dUVhYKFRZiIjIjbh4GhGRl2DAJyLyEgz4REReggFfBrg2PRGJweW0THKNYW16w3LFhrXpAXCpAyISFFv4EuPa9EQkFrbwJca16cmT5BeUYOeBiyiv1KJDiD8mje6J4bGRrb+RRMEWvsSsrUHPtelJafILSrBpz/cor9QCAMortdi053vkF3D3O7lgwJcY16YnT7HzwEXUNehMjtU16LDzwEWJSkQtMeBLTBs1BVX9MtGojoYeKjSqo1HVL5MDtqQ4hpa9vcdJfOzDlwFt1BQGeFK8DiH+FoN7hxB/CUpDlrCFT0SCmDS6J9r4mYaUNn4+mDS6p0QlopbYwiciQRiycZilI18M+EQkmOGxkQzwMsYuHSIiL8GAT0TkJRjwiYi8BAM+EZGXkM2grY+Pyi3neiPWj3WsG+tYN7bJrX6cKY9Kr9fr3VAWIiKSGXbpEBF5CQZ8IiIvwYBPROQlGPCJiLwEAz4RkZdgwCci8hIM+EREXoIBn4jISzDgExF5CcUH/JMnT6Jv3764eJEbJQOAXq9Hbm4uUlJSMGjQIEydOhWnTp2SuliSqq+vx8qVK5GYmIjExES8/vrr0Ol0rb/RC5SUlGDBggVITExEUlIS0tPTUVlZKXWxZGf58uWYMkX525AqOuDX1dVh2bJl/MfbzJYtW7BhwwYsWbIEu3btQkJCAmbMmIGSkhKpiyaZ1atX4/Dhw8jNzcWaNWuwe/durF+/XupiSU6n02Hu3LmoqanBpk2bsG7dOnz//fdYunSp1EWTlaNHj2L79u1SF0MQig74mZmZCAsLk7oYsrJjxw7MmDEDKSkp6NatGxYvXozw8HDs27dP6qJJQqvV4v3330d6ejri4+MxYsQIPPPMM9i8eTO8fRmpwsJCnD17FqtWrUKfPn0QFxeHZcuW4YsvvkB1dbXUxZMFjUaDF154AQkJCVIXRRCKDfgFBQXYtWsXlixZInVRZOX555/HhAkTzI5762P6+fPnodFoMGTIEOOxIUOGoKysDMXFxRKWTHpRUVHIy8tDeHi48ZhKpYJer2fA/83//M//YPDgwRg+fLjURRGEIgN+fX090tPTsXjxYrRv317q4sjK0KFDERn5+56i+/fvx+XLlzFixAgJSyWd0tJSBAcHIygoyHjMEOBKS0ulKpYshIaGYtSoUSbHNm7ciO7du5v8DXmr06dP45NPPsFzzz0ndVEEI5v18JvTarVW+5zbtWuHrVu3IiIiAhMmTPC6VlprdRMaGmr8uaioCM899xzuu+8+DBw4UKwiyopGo4G/v7/JsTZt2gBoGgOi3+Xm5uLzzz9Hbm6u1EWRXF1dHdLT05Genm7yb0rpZBnwz5w5g9TUVIuvjRs3DkePHsWuXbtELpU82KqbefPmYf78+QCAs2fPIi0tDb1798ZLL70kZhFlRa1WmwV2w89qtVqKIslSdnY21q5di+XLl2PkyJFSF0dy2dnZ6Nq1K8aPHy91UQSluA1QsrKykJOTY2y16fV6aDQaBAQEYPbs2Zg9e7bEJZTesWPHMGvWLMTHxyMnJwcBAQFSF0kyp0+fxkMPPYTTp08b6+Hq1atISUnBl19+iU6dOklcQum9/PLL2Lx5MzIyMvCXv/xF6uLIwl133YWysjL4+TW1ievr69HY2Ai1Wo1PPvkEd9xxh8QldI4sW/i2/Pd//zfuv/9+48+lpaV45JFHkJWVhf79+0tYMnm4cOECZs2ahWHDhmHt2rXG7gtv1adPHwQEBODEiRNITk4GABw/fhwREREM9mhqQL377rtYtWoVJk6cKHVxZGPLli1oaGgw+fn48eN48803ERERIWHJXKO4gB8aGmrSp+br6wugKePAk/ranJWeno6IiAgsX74cN2/eNB4PDAw0Gbj0Fmq1GpMnT8bKlSuxatUqaLVavPHGG5g+fbrURZNcYWEhsrOzMXPmTCQlJaGsrMz4WlhYmPHfljdq2Rho164d2rRpg65du0pUImEoLuCTdcXFxThz5gwAYOzYsSavpaWl4dlnn5WiWJJbtGgRtFotZs6cCX9/f0yePBkzZsyQuliS27t3L3Q6HXJzc80Gaj/99FP07NlTopKRuyiuD5+IiJyjyDx8IiJyHAM+EZGXYMAnIvISDPhERF6CAZ+IyEsw4BMReQkGfCIiL8GAT0TkJRjwiYi8xP8Boq2625toWOUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def get_multivariate_sample(n, d, mean):\n",
    "    x_sample = []\n",
    "    y_sample = []\n",
    "    means = np.zeros(d)\n",
    "    means[0] = mean\n",
    "    for i in range(n):\n",
    "        y = np.random.binomial(1, .5)\n",
    "        if (y == 0):\n",
    "            x = np.random.multivariate_normal(-means, np.identity(d))\n",
    "        else:\n",
    "            x = np.random.multivariate_normal(means, np.identity(d))\n",
    "        x_sample.append(x.tolist())\n",
    "        y_sample.append(y)\n",
    "    return x_sample, y_sample\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def split_by_class(sample):\n",
    "    y_one = [pair for pair in sample if pair[1] == 1]\n",
    "    y_neg_one = [pair for pair in sample if pair[1] == 0]\n",
    "    return y_one, y_neg_one\n",
    "\n",
    "X, y = get_multivariate_sample(100, 2, 2)\n",
    "sample = list(zip(X, y))\n",
    "class_one, class_neg_one = split_by_class(sample)\n",
    "plt.scatter([x[0] for x in list(zip(*class_one))[0]], [x[1] for x in list(zip(*class_one))[0]], c = 'b')\n",
    "plt.scatter([x[0] for x in list(zip(*class_neg_one))[0]], [x[1] for x in list(zip(*class_neg_one))[0]], c = 'orange')\n",
    "plt.title(\"Distribution Plot with Mean = 2\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.22476723,  0.98306512],\n",
       "        [-0.18649446, -0.38599267]]), array([[[ 1.57082944, -0.01301811],\n",
       "         [-0.01301811,  0.77251205]],\n",
       " \n",
       "        [[ 1.56240101,  0.01635323],\n",
       "         [ 0.01635323,  0.47795426]]]))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, y = get_multivariate_sample(100, 2, 1)\n",
    "estimate_x_params(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def multi_normal_entropy(x1, x2, mean, covar):\n",
    "#    return multivariate_normal.pdf([x1,x2], mean, covar)*np.log(multivariate_normal.pdf([x1,x2], mean, covar))\n",
    "\n",
    "#integrate.nquad(multi_normal_entropy, [[-20, 20], [-20, 20]], args = (0, np.identity(2)))\n",
    "\n",
    "def multivariate_normal_entropy(covar):\n",
    "    covar = np.array(covar)\n",
    "    dim = covar.shape[-1]\n",
    "    return multivariate_normal(np.zeros(dim), covar).entropy()\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy_f(t1, t2, t3, t4, mean_1, mean_2, covar_1, covar_2):\n",
    "    t = (t1, t2, t3, t4)\n",
    "    return -.5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2)) - .5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2))\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy(mean_1, mean_2, covar_1, covar_2):\n",
    "    dim = 1\n",
    "    # to handle dimension errors\n",
    "    try:\n",
    "        dim = len(mean_1)\n",
    "    except TypeError:\n",
    "        pass\n",
    "    bounds = [[-5, 5] for _ in range(dim)]\n",
    "    return integrate.nquad(two_mixture_multivariate_normal_entropy_f, bounds, args = [mean_1, mean_2, covar_1, covar_2])\n",
    "\n",
    "def plugin_estimate_cat_multi(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    print(\"integrating\")\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    h_x_cond_y = multivariate_normal_entropy(x_params[1][0])*.5 + multivariate_normal_entropy(x_params[1][1])*.5\n",
    "    h_x = two_mixture_multivariate_normal_entropy(x_params[0][0], x_params[0][1], x_params[1][0], x_params[1][1])\n",
    "    return h_y - h_x + h_x_cond_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integrating\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.29961871, 6.3508492 ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, y = get_multivariate_sample(10000, 4, 1)\n",
    "estimate_x_params(X)\n",
    "plugin_estimate_cat_multi(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integrating\n"
     ]
    }
   ],
   "source": [
    "#def multi_normal_entropy(x1, x2, mean, covar):\n",
    "#    return multivariate_normal.pdf([x1,x2], mean, covar)*np.log(multivariate_normal.pdf([x1,x2], mean, covar))\n",
    "\n",
    "#integrate.nquad(multi_normal_entropy, [[-20, 20], [-20, 20]], args = (0, np.identity(2)))\n",
    "\n",
    "def multivariate_normal_entropy(covar):\n",
    "    covar = np.array(covar)\n",
    "    dim = covar.shape[-1]\n",
    "    return multivariate_normal(np.zeros(dim), covar).entropy()\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy_f(t1, t2, t3, t4, t5, mean_1, mean_2, covar_1, covar_2):\n",
    "    t = (t1, t2, t3, t4, t5)\n",
    "    return -.5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2)) - .5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2))\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy(mean_1, mean_2, covar_1, covar_2):\n",
    "    dim = 1\n",
    "    # to handle dimension errors\n",
    "    try:\n",
    "        dim = len(mean_1)\n",
    "    except TypeError:\n",
    "        pass\n",
    "    bounds = [[-5, 5] for _ in range(dim)]\n",
    "    return integrate.nquad(two_mixture_multivariate_normal_entropy_f, bounds, args = [mean_1, mean_2, covar_1, covar_2])\n",
    "\n",
    "def plugin_estimate_cat_multi(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    print(\"integrating\")\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    h_x_cond_y = multivariate_normal_entropy(x_params[1][0])*.5 + multivariate_normal_entropy(x_params[1][1])*.5\n",
    "    h_x = two_mixture_multivariate_normal_entropy(x_params[0][0], x_params[0][1], x_params[1][0], x_params[1][1])\n",
    "    return h_y - h_x + h_x_cond_y\n",
    "\n",
    "X, y = get_multivariate_sample(10000, 5, 1)\n",
    "estimate_x_params(X)\n",
    "plugin_estimate_cat_multi(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def multi_normal_entropy(x1, x2, mean, covar):\n",
    "#    return multivariate_normal.pdf([x1,x2], mean, covar)*np.log(multivariate_normal.pdf([x1,x2], mean, covar))\n",
    "\n",
    "#integrate.nquad(multi_normal_entropy, [[-20, 20], [-20, 20]], args = (0, np.identity(2)))\n",
    "\n",
    "def multivariate_normal_entropy(covar):\n",
    "    covar = np.array(covar)\n",
    "    dim = covar.shape[-1]\n",
    "    return multivariate_normal(np.zeros(dim), covar).entropy()\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy_f(t1, t2, t3, t4, t5, t6, t7, t8, mean_1, mean_2, covar_1, covar_2):\n",
    "    t = (t1, t2, t3, t4, t5, t6, t7, t8)\n",
    "    return -.5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2)) - .5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2))\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy(mean_1, mean_2, covar_1, covar_2):\n",
    "    dim = 1\n",
    "    # to handle dimension errors\n",
    "    try:\n",
    "        dim = len(mean_1)\n",
    "    except TypeError:\n",
    "        pass\n",
    "    bounds = [[-5, 5] for _ in range(dim)]\n",
    "    return integrate.nquad(two_mixture_multivariate_normal_entropy_f, bounds, args = [mean_1, mean_2, covar_1, covar_2])\n",
    "\n",
    "def plugin_estimate_cat_multi(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    print(\"integrating\")\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    h_x_cond_y = multivariate_normal_entropy(x_params[1][0])*.5 + multivariate_normal_entropy(x_params[1][1])*.5\n",
    "    h_x = two_mixture_multivariate_normal_entropy(x_params[0][0], x_params[0][1], x_params[1][0], x_params[1][1])\n",
    "    return h_y - h_x + h_x_cond_y\n",
    "\n",
    "X, y = get_multivariate_sample(10000, 8, 1)\n",
    "estimate_x_params(X)\n",
    "plugin_estimate_cat_multi(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integrating\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.32172777, 3.48004234])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#def multi_normal_entropy(x1, x2, mean, covar):\n",
    "#    return multivariate_normal.pdf([x1,x2], mean, covar)*np.log(multivariate_normal.pdf([x1,x2], mean, covar))\n",
    "\n",
    "#integrate.nquad(multi_normal_entropy, [[-20, 20], [-20, 20]], args = (0, np.identity(2)))\n",
    "\n",
    "def multivariate_normal_entropy(covar):\n",
    "    covar = np.array(covar)\n",
    "    dim = covar.shape[-1]\n",
    "    return multivariate_normal(np.zeros(dim), covar).entropy()\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy_f(t1, t2, mean_1, mean_2, covar_1, covar_2):\n",
    "    t = (t1, t2)\n",
    "    return -.5*multivariate_normal.pdf(t, mean_1, covar_1)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2)) - .5*multivariate_normal.pdf(t, mean_2, covar_2)*np.log(.5*multivariate_normal.pdf(t, mean_1, covar_1) + .5*multivariate_normal.pdf(t, mean_2, covar_2))\n",
    "\n",
    "def two_mixture_multivariate_normal_entropy(mean_1, mean_2, covar_1, covar_2):\n",
    "    dim = 1\n",
    "    # to handle dimension errors\n",
    "    try:\n",
    "        dim = len(mean_1)\n",
    "    except TypeError:\n",
    "        pass\n",
    "    bounds = [[-10, 10] for _ in range(dim)]\n",
    "    return integrate.nquad(two_mixture_multivariate_normal_entropy_f, bounds, args = [mean_1, mean_2, covar_1, covar_2])\n",
    "\n",
    "def plugin_estimate_cat_multi(X, y):\n",
    "    y_param = estimate_p(y)\n",
    "    x_params = estimate_x_params(X)\n",
    "    print(\"integrating\")\n",
    "    h_y = -y_param*np.log(y_param) - (1 - y_param)*np.log(1 - y_param)\n",
    "    h_x_cond_y = multivariate_normal_entropy(x_params[1][0])*.5 + multivariate_normal_entropy(x_params[1][1])*.5\n",
    "    h_x = two_mixture_multivariate_normal_entropy(x_params[0][0], x_params[0][1], x_params[1][0], x_params[1][1])\n",
    "    return h_y - h_x + h_x_cond_y\n",
    "\n",
    "X, y = get_multivariate_sample(10000, 2, 1)\n",
    "estimate_x_params(X)\n",
    "plugin_estimate_cat_multi(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/sklearn/ensemble/weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.\n",
      "  from numpy.core.umath_tests import inner1d\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.6826752373237107"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.ensemble.forest import _generate_unsampled_indices\n",
    "from sklearn.ensemble import BaggingClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "X, y = get_sample(1000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "\n",
    "def rf_entropy_estimate_bagging(X, y, n_estimators = 200, max_samples = .32, bootstrap = True, depth = 30, min_samples_leaf = 1):\n",
    "    model = BaggingClassifier(DecisionTreeClassifier(max_depth = depth, min_samples_leaf = min_samples_leaf), \n",
    "                              n_estimators = n_estimators, \n",
    "                              max_samples= max_samples, \n",
    "                              bootstrap = bootstrap)\n",
    "    model.fit(X, y)\n",
    "    class_counts = np.zeros((X.shape[0], model.n_classes_))\n",
    "    for tree in model:\n",
    "        # get out of bag indicies\n",
    "        unsampled_indices = _generate_unsampled_indices(tree.random_state, len(X))\n",
    "        \n",
    "        total_unsampled = len(unsampled_indices)\n",
    "        np.random.shuffle(unsampled_indices)\n",
    "        prob_indices, eval_indices = unsampled_indices[:total_unsampled//2], unsampled_indices[total_unsampled//2:]\n",
    "        # get all node counts\n",
    "        node_counts = tree.tree_.n_node_samples\n",
    "        # get probs for eval samples\n",
    "        posterior_class_counts = np.zeros((len(node_counts), model.n_classes_))\n",
    "        for prob_index in prob_indices:\n",
    "            posterior_class_counts[tree.apply(X[prob_index].item()).item(), y[prob_index]] += 1\n",
    "        row_sums = posterior_class_counts.sum(axis=1)\n",
    "        row_sums[row_sums == 0] = 1\n",
    "        class_probs = (posterior_class_counts/row_sums[:, None])\n",
    "        where_0 = np.argwhere(class_probs == 0)\n",
    "        for elem in where_0:\n",
    "            class_probs[elem[0], elem[1]] = 1/(2*row_sums[elem[0], None])\n",
    "        where_1 = np.argwhere(class_probs == 1)\n",
    "        for elem in where_1:\n",
    "            class_probs[elem[0], elem[1]] = 1 - 1/(2*row_sums[elem[0], None])\n",
    "        class_probs.tolist()\n",
    "        partition_counts = np.asarray([node_counts[x] for x in tree.apply(X[eval_indices])])\n",
    "        # get probability for out of bag samples\n",
    "        eval_class_probs = [class_probs[x] for x in tree.apply(X[eval_indices])]\n",
    "        eval_class_probs = np.array(eval_class_probs)\n",
    "        # find total elements for out of bag samples\n",
    "        elems = np.multiply(eval_class_probs, partition_counts[:, np.newaxis])\n",
    "        # store counts for each x (repeat fhis for each tree)\n",
    "        class_counts[eval_indices] += elems\n",
    "    # calculate p(y|X = x) for all x's\n",
    "    probs = class_counts/class_counts.sum(axis = 1, keepdims = True)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    # convert nan to 0\n",
    "    entropies = np.nan_to_num(entropies)\n",
    "    return np.mean(entropies)\n",
    "\n",
    "np.warnings.filterwarnings('ignore')\n",
    "\n",
    "rf_entropy_estimate_bagging(X, y, 500, .32)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEgCAYAAAB1t06HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XlYVFUfwPHvzLCDiArihqYoLqggKribiqLgllmZuaUtkOVWqW1mlqmYy2uamdmiZma5BGq5L+WCa+4sCqKAgKDsyzAz5/2DmBwBHRAR5Xyeh4eZu51z7p2Z373nnnuOQgghkCRJkqRSUD7qDEiSJEmPLxlEJEmSpFKTQUSSJEkqNRlEJEmSpFKTQUSSJEkqNRlEJEmSpFKTQeQ+evbsydq1awtNj4mJoWnTpoSHhxtM/+yzz9i4cSNffPEFrVq1Ijo6utC6W7ZsoXnz5pw5cwaApk2bEhMTA8C7775L586duX37tsE6V65cwc3NjS1btgAQEBBASEjIffPetGnTIv++/PLL+5ZdCMGGDRvIzc0FYNOmTXh5ed13vdI6duwYly5demjbL4mdO3fy/PPP4+bmhoeHB6NGjeLIkSOPLD8hISE0bdqUzMzMUq2v0Wj45ptv8PHxoWXLlnTs2JEpU6YYfD6//PJLhgwZUlZZLpWRI0cyb948APLy8vjpp5/080qav4J95uXlhVarLTR/zZo1NG3aVJ/e40StVjNgwAD27dv3qLMCQrqnHj16iDVr1hSafv36deHi4iLCwsL0086dOyf8/PyEVqsV2dnZomfPnmLs2LEG66WmpoqOHTuK2bNn66e5uLiI69evCyGESElJEZ07dxZTpkzRz8/NzRWDBw8Wb731ln5aaGio6Nevn8jNzb1n3r/++muRmJhY6C8jI+O+ZQ8JCREuLi76ZbOzs0VSUtJ91ystFxcXsXfv3oe2fWPt27dPuLm5iV9//VVcvXpVREREiPnz54sWLVqIEydOPJI8HT161OBYlNS8efNE7969xf79+8X169fFmTNnxGuvvSY6d+4sUlNThRBCZGRkiFu3bpVltkvs9u3bIj09XQghxKZNm4Snp6d+3pIlS8Qzzzxj9LYK9pmrq6s4evRoofkjRowQTZs2FXPnzn3wjJej7OxsERAQUGG+L/JKpAwtXbqUYcOGoVQqsbCw4OOPP+bvv//mjz/+0C+zcOFCLC0tmTx5cpHbqFq1KrNmzWLr1q3s2bMHgMWLF3Pz5k0++eQT/XJNmzbF3t6ebdu23TNP1tbWODg4FPqztra+b3nEXc+hWlhYUKNGjfuu97jbuHEjfn5+DB06lAYNGtC4cWPeeecdPDw8+PXXXx919kpl48aNjB8/nu7du1OvXj1at27N4sWLyczMZPfu3UD+Z6VatWqPNJ92dnbY2NgAhT9/peXp6cmuXbsMpt26dYtz587RokWLMkmjvFy4cIGhQ4cSFxf3qLOiJ4NIGYmJieHgwYP07t1bP61bt274+voyd+5ccnJyuHjxIhs2bOCzzz7D0tKy2G317NmTQYMGMXPmTA4ePMgPP/zA559/XugL3qdPH9asWfNA+U5ISMDf35+2bdvSrl07JkyYQHJyMjExMYwaNQoADw8PNm3aZFCdVVCdd+DAAfr06YObmxsTJ07kxo0bBAQE4Obmhq+vL6dPn9an9ddff/H888/TunVr3NzcGDVqFFFRUfoyA/j7+zN9+nQAoqKiGDduHG5ubvTo0YPAwEDUanWR5Rg+fDhz5swxmPbpp5/yyiuvALB792769+9Pq1at6NmzJ99++22x+0ShUHD27FlSU1MNpn/xxRe8/fbb+vdBQUEMGDCAli1b4uHhQUBAAElJSUB+VUrnzp3Ztm0b3bp1o02bNsyaNYvo6GhGjBhB69atGTp0KJGRkfrl27Vrx6ZNm+jSpQvt2rXjgw8+IDs7u8g83rx5k4kTJ9KmTRu6dOnCBx98QHp6+j3LdPToUTQajX6apaUlW7Zs0X9m76wuGjlyZJHVoJs3bwaMPzY6nQ4vLy+Dk53XX38dX19f/fszZ87g7u5Obm6uvjorJCSE9957j5SUFJo2baqvutXpdAQGBuLl5UXbtm2ZMWMGeXl5xZYbwNvbm7179xpM2717Nx06dCh0MnXgwAEGDRpE69at8fPzY+PGjQbzf/jhB32VoKenJ1OnTiUrKwvIr+4dMmQIK1eupHPnzri7u/P222/r599t+vTpxVY3F1Rv3+3w4cN069aN9evX37PM5epRXwpVdMZWZ61du1b079+/0HKJiYmiXbt2YtmyZeKll14S7733XqFl7qzOKlBQrdWiRQvx8ccfF5m3q1evChcXF3Hz5s0S5f1OAQEBYuzYseLy5cvi0qVLYsiQIeKdd94RGo1G7NixQ7i4uIhr166J7OxssXHjRn31QkH5n3nmGXHu3Dlx+PBh4erqKry8vMTGjRtFRESEePnll8XgwYP1y7u6uopVq1aJa9euidOnT4uBAwcKf39/IYQQycnJwsXFRWzZskWkpaWJnJwc0aNHDzFr1ixx5coVcezYMdG/f3/x0UcfFVmOdevWie7duwudTieEEEKr1YrOnTuLzZs3i6SkJOHq6ipWr14tYmJixB9//CFcXV3F4cOHi9zWsWPHhKurq3BzcxP+/v7ihx9+EJcvXzZY5uTJk8LV1VVs3rxZxMTEiMOHD4vu3buLTz/9VAiRX5XSokULMXr0aBEeHi6CgoKEi4uL6NGjh9i1a5e4dOmSGDBggBg/frx++ebNmwtfX19x8uRJcezYMdGzZ08xbdo0/fw7q7NeeOEFMWHCBBEWFibOnj0rRowYIcaNG1fscV6xYoVwcXERnTp1EtOmTRObNm0q9Lm5s7ro9u3bBtWfEyZMED4+PiI9Pb3Ex2bKlCniww8/FEIIodFohIeHh2jWrJlITk4WQgjx5Zdfitdee00IkV/FNHfuXJGbmyt++OEH4enpKRITE0Vubq5YsmSJcHFxEe+//76IiooSu3fvFq6uruKXX34pMt2CfXbjxg3h6uoqzp8/r583btw4sWnTJn16QggRHh4uWrduLdavXy+io6PFtm3bRPv27cXWrVuFEEIEBQUJDw8PsXfvXhETEyN2794t2rRpI77//nshhBAbN24Urq6uwt/fX0RERIj9+/eL1q1b6+ffLS0trciq5sTERKHRaIo9lgUqSnWWDCL30aNHD+Hq6irc3d0N/tzc3AyCyPTp08XkyZOL3Ma6detEixYtDOqf71RUEBFCiAkTJggXFxexbdu2Irer1WpFy5Yti/0gFZd3d3d3ERsbK4QQYsCAAWLChAkiJydHCJEfmC5cuCCEKPzDVVQQ2bFjhz69oUOHioCAAP377du3C1dXVyGEEFFRUWL16tUG+Vu1apXo1auXwX4oKMtvv/0m+vTpY7D8yZMnRbNmzfR15ne6deuWcHV1FSdPnhRCCHHkyBHRunVrkZGRIS5cuCBcXFzEn3/+qV/++PHjxQZfIYQ4f/68mDx5smjbtq1wcXERLi4uYsyYMfp1zp8/LzZt2mSwzqxZs8SoUaMM9t2dP1wdO3YUn332mf79ypUr9WUsWP748eP6+bt27RKurq4iPT3d4FgcOXJEuLu7G9wPi4+PFy4uLiI8PLzYMv35559i9OjRwtXVVbi4uIjmzZuLWbNm6X+wirvnsH79etGmTRt9IC3psdm8ebPw9vYWQghx5swZ0bNnT+Hr66v/7Lzwwgvip59+EkIIgx/1Oz9vBflr3769yMvL008bPXp0sSdZd+6zsWPHisWLFwsh8n+83d3dRUpKikF6U6dO1Qe7AsuXLxdDhgwRQuR/pnbu3Gkw//XXX9efGG7cuLHQSd348eMN7m+WpYoSREwe9ZXQ4+D1119n4MCBBtMSEhIYOXKk/n1SUhL169cvcv1hw4axbNkyxowZg62trVFp7t69m507d9K+fXtmz55Nx44dC1VnKZVKqlatSnJyconyDlCzZk0gv/po2rRpeHl50aFDB7y9vYtcvjhOTk761xYWFoXeF1Q1PPXUU1haWrJy5UoiIiKIiori0qVL+nzc7fLly1y/fp02bdropwkh0Ol0XL16lZYtWxosX61aNbp06cIff/yBh4cH27dvp0ePHlhbW9O8eXP69OnDhAkTqFu3Lt27d2fgwIHY29sXWy5XV1cWLlyIRqPh7Nmz7Ny5k3Xr1jFt2jRWrVqFq6sr1tbWLFu2jCtXrnDlyhUiIiJo27ZtqfYPgEqlwt3dXf++VatW5OXl6au87tw32dnZRbaUi4yMpEmTJkWWycfHBx8fHzIzMzl27Bhbtmxh7dq1ODo68tprrxW5zqlTp/jss89YsGABzs7O+vRLcmy6du3K9OnTiYuL4+jRo3h6eqJQKDh+/DgdO3bk3LlzLFiwoMj071anTh1MTP772bK1tdW3HrwXb29v1q1bx8SJE9m7dy8eHh5UrVrVYJmIiAjCw8PZunWrfppGo9Gn16FDBy5cuMDixYuJjIzk8uXLREZGMnjwYP3y1tbWBp8rGxubYquzZsyYQXBwcJHztm3bRp06de5bropABhEjVKtWjQYNGhhMU6lUBu+VSiU6na7I9RUKBWZmZlhYWBiVXlJSEh999BEjR47E398fX19fPv30UxYuXFhoWa1WWygv98v7nXx9fenQoQP79u3jr7/+Yvbs2QQHB/Pjjz8aldei9kNRwsLCePHFF+nUqRPt27dn6NChnDlzhp9//rnI5TUaDe7u7oXucwA4OjoWuU7//v0JDAxk6tSp7Nq1i08//RTI3/9ffvkloaGh7Nu3j/3797N+/Xpmz55dqMloZmYmCxcuZPTo0dSvXx8TExM8PDzw8PCgQYMGfPrpp2i1WkJCQvR1++3atWPUqFFs3bqVsLCwUu2fgnl3zi/4PN29DY1GQ506dfj+++8LbaOohg+hoaH8+uuvfPTRR0D+D12PHj3o0aMHkydP5q+//ioyiCQkJDBhwgRefvll+vTpY5B+SY5NjRo1cHV15ejRo4SEhODn54dSqWT16tUcOXKEhg0bUrdu3eJ2i4Gi9p8w4ga8t7c3s2bN4tq1a+zevdugPAW0Wi0jR45k2LBhRW5j06ZNzJw5kyFDhtC1a1cCAgIKNZU3NTU1qhwAEydOZNy4cUXOK+7kqiKSN9bLiL29faFnO0rrww8/xNbWlilTplC9enVmzJjBtm3b2Llzp8FyOp2O1NTUe55R38/ixYuJiYnh2WefZfHixSxdupSjR4+SlJSEQqF40KLobdiwgebNm7N06VJGjx6Np6cnsbGxxf4AODs7Ex0dTa1atWjQoAENGjTg9u3bBAYGFnsjtVevXqSnp7N69Wq0Wi3dunUD8p+x+eyzz2jWrBkBAQH88ssv+Pr6sn379kLbsLCwIDg4mKCgoELzbGxssLW1RaVSsXr1avr168e8efMYPnw47u7uREdHP1CLory8PIPnjs6cOYO5uTkNGzYstG8SExOxtrbW7xsTExPmzJnDrVu3Cm1Xp9Oxdu1ajh07VmSZimqRpVarefPNN2natCmTJk0qlH5Jj023bt04dOgQp0+fpn379rRv357Q0FC2b9+uP053K8vPn4ODA+7u7mzbto1Dhw7h7e1daJmCchWUqUGDBhw5ckT/nNj333/PuHHjmDlzJs899xzNmjV7oGNeo0YNg7Tu/Lvzaquik0GkjLi6uhIaGvrA2/nll184cOAAc+fO1V+5+Pr60rt3b2bOnGnwIxEREYEQAldX12K3l5mZyc2bNwv9paSkAPnVH7NmzeLcuXNER0ezdetW6tatS/Xq1bGysgLymxWW9iG3Ao6OjkRGRnLixAmuX7/ODz/8wK+//mrQosfKyoqIiAhSUlIYOHAgSqWSadOmER4ezqlTp3jvvffIy8ujSpUqRaZhaWlJr169WLp0KX369MHMzAzIbza9ceNGFixYwPXr1zl16hT//PMPrVu3LrQNlUrFW2+9xfLly1m8eDFhYWFERUURHBxMYGAgr776KpB/pnj27FkuXLhAVFQUixYt4uDBg8W2HjPWRx99xIULFzh69CiBgYE899xz+uNQoHPnzjRp0oTJkydz/vx5Ll26xNtvv01sbGyRZ/QtWrSgT58+TJw4kV9//ZVr165x6dIlVq1aRXBwMC+//HKhdT7++GOSk5P1n7mCz016enqpjk23bt3YsWMHtra2ODk5UbduXWrVqsWff/7J008/XeQ6VlZWZGVlcfnyZaOqrO6nd+/erFy5EldX1yKv2MaOHcv+/fv5+uuviY6O5o8//mDevHn6q6uaNWsSEhLC5cuXiYiI4MMPP+Ty5csPfMwfdzKIlJHu3btz9epVbt68WeptXLt2jblz5/Lyyy8b1DcDzJw5E61Wa/CsyIkTJ3Bzc6N69erFbnPhwoV06dKl0N9bb70FwCeffEKDBg145ZVXGDhwIHFxcaxYsQKlUomLiws9evRg7NixbNiwodTlgvwmo506dcLf358hQ4awb98+Zs6cSXJyMvHx8QCMGzeOZcuW8cEHH2BlZcV3331HWloazz33HG+88Qbu7u588cUX90xnwIABZGVl0b9/f/00e3t7vvrqKw4fPkz//v1588036dmzJ/7+/sXmdd68eRw9epThw4czcOBAVq1axYQJE/TVDxMmTKB+/fqMGDGCF198kfDwcKZNm/bAP3h+fn6MGzeOiRMn4uPjo2/ufCelUsny5cuxs7Nj1KhRjBw5EgcHB1auXFls1eaCBQsYOXIkP/zwAwMGDGD48OH89ddffPvtt4U+a5BfdRMbG4u3tzedO3fWf25mz55dqmPTunVrrK2tad++vX6ap6cnNjY2eHh4FLlOx44dadGiBYMHD2b//v332XP317t3bzIzM4usygJo2bIlS5YsYfv27fj5+REYGIi/v7/+mH/wwQcoFAqeffZZXn75ZdRqNa+//joXL1584Lw9zhTiQa6/JQOvvfYaHTp0YOzYsSVar2nTpuzZs4d69eqVaL1hw4YxbNgwgxt70uMpJCSEUaNGcerUKaMeBJWkikJeiZSh8ePHs2HDhiL76SlrFy9eJDk52eCMW5IkqbyVaxAJCwtj2LBhuLu707dvXw4cOFDkclFRUYwZM4Z27drRtWtXFi1aVGZdIDxMbm5udOrUiU2bNj30tP73v/8xa9asx+oGnCRJT55yq85Sq9X07duXkSNHMmLECA4ePMg777yjv5F7p2effZZevXrh7+9PbGwsI0aMYPLkybLaRpIkqYIptyuRkJAQcnJyGDNmDKampvTq1QtPT88iH7aJjIxEp9Pp28krFArMzc3LK6uSJEmSkcotiFy5cgVnZ2eDtt+NGjUiIiKi0LIBAQF89dVXtG7dGm9vbzp16kS/fv3KK6uSJEmSkcotiGRlZRV6YtvCwqLIXkpNTEyYOnUqp0+fZsuWLRw+fLhi9VopSZIkAeXY7YmVlRU5OTkG03Jycgo9SHX+/Hm+/fZbDh06hEKhoHnz5vrnFIrrjqAot29notMJatSwITk5o0zK8LipzGWHyl1+WfbKWXYoffmVSgXVqpW8eXm5BRFnZ+dCYzhERkYWetApLi6uUNcJJiYmJW6FpNMJdDqhf11ZVeayQ+Uuvyx75VWe5S+36iwvLy9UKhXffPMNarWavXv36jtju1PB06uLFy9Go9Fw9epVvv/++0LLSZIkSY9euQURMzMzVq5cyf79++nQoQPz589n0aJFODk5ERQUpL8isbe3Z+XKlRw/flz/9PfQoUMNul2XJEmSKoYnttuT5OQMdDqBg0MVbt4sftjQJ1llLjtU7vLLslfOskPpy69UKqhRw6bE68nHnSXpCSOEIDo6mrS0DOCJPEe8p8TE4sf2qQyKL78CMzMLqlVzKNNu9mUQkaQnTEZGKkLocHSsh0JR+brHMzFRotFU3iBSXPmF0JGSkkRGRipVqtiVWXqV7xMmSU+47OwMbG2rVcoAIhVPoVBSpUo1srPLtvmz/JRJ0hNGp9OiUslKBqkwlcoEna5sexmXQUSSnkBlWectPTkexudCBhFJkp546enppKWlPdQ0cnNzSE5OeqhpVEQyiEiSVGHEx8czevSL9O7dlRUrlpXZdocNG0xcXAwAq1d/x4cfTiuzbRd4441XOXv2nzLfbkUnK04lSaow/vnnJKmpKfzxx74yHXAtNTVV/3rUqJINX218GikPZbsVnQwikiQ9EjduxDFq1At4e/uwb99uXn31DZYty+/uqF+/HqxatZb69RsYrHP27D8sXbqY6OgoateuyxtvTMDTswMAW7duYfXq70lPT6NOnXq88oo/HTt2ZvToFwEYP/5Vpk37kOvXrxEWdonAwMWsWrWCuLhYsrIyOXHiGI6OtZk+/UPWr19LSMgRatWqzYwZn9GkiQu5ubl89dX/CAk5QlLSTapVq86rrwbQp08/pk6dTEJCPJ9++jHXr19j1Kix/P33AVau/JqEhBs89VQjJk16h2bNWpT7fn7YZBCRpCfcnyHX+P1QFLnqsm2VUxRzMxWDOjekr1d9o5bPzs7Gzq4awcG7yMvLw9LSkp9/XsOaNRsKLZuQEM/bb09g6tT36dmzNydPHmfGjOl8++0arK1t+OKLufzww880buzM5s2bCAyczaZN2/jxx5/p0qUdy5atpFmzFqxatcJgu7t372DhwqV8+uk83nvvbd588zXmzFnAxx/P5rPPZvDjj6v47LN5rF+/lvDwUL755kdsbGzYuHEDCxbMpWfP3gQGLmLo0AGMHz+RHj28CQ29yCeffMjcuQtxc2vDnj07efvtt/j5583Y2tqWyb6uKOQ9EUl6wu04fq1cAghArlrLjuPXSrSOt7cPpqamhYaFuNuuXX/SsmUrevfui0qlwtOzA15endi2LQilUoFSqeL33zdx/vw5fH0HsGnTNqNaIzVr1oK2bdtjYmKCm1sbnnqqER07dsbU1BQPj/YkJNwA4JlnhvL55wuwsbEhMTEBCwsLMjMzycrKLLTNrVt/x9u7r367Pj6+1KtXn/3795Ro3zwO5JWIJD3hfNrXL9crEZ/2xl2FFKhevYZRy8XHx/PPP6fo2/dp/TStVku3bj2oWtWOJUuWs3btD0yY8AZmZqY8//xwRo58+b6B5M4rA6VShY2NzR3vlfpu1TMzM1m4cB7nz5+jdu06NGjwFFB0t+vx8Tc4deok+/bt0k/TaDQkJMQbVdbHiQwikvSE6+tV3+jqpUfB2GcX7O3t6dy5G599Nk8/LT4+HisrS9LS0tBqdcyduxDQcfToUd5//x1atmyNh0e7Mkk/MHA2devWY86cBZiYmBAWFsquXX8Wk1cHnn32ecaPn6ifFhsbg51d2XU3UlHI6ixJkh4LvXv35fjxoxw+/Dc6nY6IiDBeeWUkf/11gJSU20yePJ5Tp05gYmJCjRr5VzdVq+b/aJuampKZWbjaqSQyMzMxNzdHqVSSlJSkb4Ks0WgKpdGvX3+2bw/i/PlzCCE4efI4o0a9QFhY6APloSKSVyKSJD0WCq4Cli//kk8++QAbmyoMG/YSfn4DAZg27QMCAz8nOfkmdnbVmTTpXZydGwPg5zeIqVMnMX78pFKnP3HiOwQGfoaPz9PY2dnRv/8gwsNDiY6Owt7eHl/fASxePJ+YmOv4+7/J22+/x/z5s7lx4wbVq9dgypRp970qehzJ8USeYJW57FB5yx8fH029eg0rbU+2shffe5c/Pj6aWrUaFJpe2vFEZHWWJEmSVGrlGkTCwsIYNmwY7u7u9O3blwMHDhRa5sSJE7Rp08bgz9XVFR8fn/LMqiRJkmSEcgsiarWagIAAfHx8OH78OO+++y6TJk0iNjbWYLl27dpx+vRp/d8ff/yBnZ0dH374YXllVZIkSTJSuQWRkJAQcnJyGDNmDKampvTq1QtPT0+Cg4Pvud4HH3zAgAED6Nq1aznlVJIkSTJWubXOunLlCs7OzgZtshs1akRERESx6+zevZuwsDCWLl1aHlmUJEmSSqjcgkhWVhYWFhYG0ywsLMjOzi52neXLl/PKK69gaWlZ4vTubGXg4FClxOs/KSpz2aFylj8xMb+CwcSk8rabqcxlh3uXX6lUlun3otyCiJWVFTk5OQbTcnJyiu0vJzQ0lIiICJ599tlSpSeb+FbeJq4FKmv5dbr85p2VtZmrbOJ77/LrdLoivxcVvomvs7MzUVFRBtMiIyNp3Lhxkcvv2bOHbt26UaVK5TuTlCRJelyUWxDx8vJCpVLxzTffoFar2bt3LyEhIfj5+RW5/JkzZ3B3dy+v7EmS9AhU1iFlnyTlFkTMzMxYuXIl+/fvp0OHDsyfP59Fixbh5OREUFAQbdq0MVg+NjYWBweH8sqeJEmPQGUdUvZJUq59Z7m4uLBu3bpC0wcOHMjAgQMNpm3btq28siVJ0iNSWYeUfZJU7iYMkiQ9MncOKbt69XcEBIwlIGAsvr69CAsLZejQAezbt1u//KpVK5g69b8OFIOCNvPii0Po27cHkyePJybm+qMoRqUne/GVpCdc2uFDpP59sFzSqtqlG7adOhu17J1DymZnZ3Pu3Fm++GIJrVu7YWl571EODxzYy7fffs0XX/yPhg2dWb/+J959dyJr1mzAxMSsLIoiGUleiUiSVCHY2lalQ4dOWFlZ33egqKCgLQwd+gIuLs0wNTVl5MgxZGVlcfr0iXLKrVRAXolI0hPOtlNno68OHiVjh8kFSEi4werVp1i3brV+Wl5eHvHxT97wsxWdDCKSJFUId198KJVKtNr/xoVPS0vVv65Rw4FnnnmOZ599Xj/t2rWr1KxZ66HnUzIkq7MkSXpk7jVsrZNTA/bv34NGoyEy8jL79+/Vz/P17c+GDeuIiopECMGuXX8yatQwEhMTyivr0r/klYgkSY9MwZCy1tbW2NpWNZjn7z+e+fPn4OvbC2dnZ3x9B3DlSn6HrT4+vmRkZPD++++QlJREnTp1mT17PvXrFx6xT3q45PC4T7DKXHaovOWXw+PKvrPk8LiSJEnSY0EGEUmSJKnUZBCRJEmSSk0GEUmSJKnUZBCRJEmSSk0GEUmSJKnUZBCRJEmSSk22XmU8AAAgAElEQVQGEUmSJKnUZBCRJEkqQ1qttlw6goyLi33oaRijXINIWFgYw4YNw93dnb59+3LgwIEil8vIyGD69Ol4enri5eXFjBkzyMvLK8+sSpL0iL355musW7fmgbfTu3dXwsNDyyBHxvn44/fYu3cXAGfOnGbgQJ8yT2Pp0sWsX7+2zLdbGuUWRNRqNQEBAfj4+HD8+HHeffddJk2aRGxs4Wj6/vvvk5aWxt69e9m+fTvnz5/nu+++K6+sSpL0BNm16y9cXJqVW3opKf8N+evm1oagoB1lnkZFGla43IJISEgIOTk5jBkzBlNTU3r16oWnpyfBwcEGyyUmJrJnzx5mzZqFjY0NNWrU4KuvvqJ///7llVVJksrBjRtx9OzZiZ9++hFf314MGNCHNWt+KHLZu69Ktm8PZuTI/7qBX736OwYN8mHQIB+++upLnntuIKdO5Q9Q1aVLO0JDL+pfb9z4C889N5C+fZ9m+vQpZGRkFJlmYmIC7733Dn5+vXjuuUH8/PNaCroaPH/+LOPGjaRv36cZPvxZfvrpRwAWLpzH2bP/sHLlV8yf/zmnTp2gd++uAJw6dYIRI57nm2++ol+/ngwa5MOePTtZuXI5fn69GDTIh507/9Snv2nTr4wePQwfn+707+/Nl18uAuCnn35k584/CArazNtvTwAgKiqSiRPfoG/fHrzwwhB27Nhe4uNRWuXWi++VK1dwdnY2GLGsUaNGREREGCx36dIlHB0dCQ4OZs2aNWg0GgYNGsTEiRPLK6uS9ET56p8vmX98Dpl5Rf9YliVrUxvebf8eb7i/ZdTyarWas2f/4bffgrhx4waTJr2Bo6Mjffr0MzrNP//cxsaNG1i0aBl169bjf/9bwI0bccUuf+jQX6xatZbc3BzeeONVNm/+jZEjxxgso9VqmTZtMm5uHmzatJ3ExASmTp2EnZ0d/fr1JzBwNkOHDmPgwGe4cuUyAQHj6NSpK1OmTCMy8gqdOnVl+PCR+kBW4OrVSLy9+7Bt225++ulHZs36iJdeGs3vv+9gy5bfWLLkC/r06cvZs//w7bdf89VX3/LUUw25dOkC/v5j6dmzNy+9NJqrV6OwtLRkypRpZGVlMXnyeJ5//kUWLFhCZGQEb789kZo1HWnTpq3R+7G0yi2IZGVlYWFhYTDNwsKC7Oxsg2kpKSnEx8dz+fJlgoODuXXrFq+//jrW1tb4+/sbnd6dvVE6OFR5sMw/xipz2aFylj8xMb+CwcQk///XZ74slwACkJmXwddnvmRCu/uf9KlU+fmbPPkdbG2rYGtbhUGDnmH37h34+vqhUChQqRSYmCgNXkN+j7MKRf77P//cxvPPD8PFxQWASZPeZtu2IFQqpX75O18PGzac6tXtAPDy8iI29pp+XoHQ0Atcv36N779fg4mJKQ0bPsXw4SMICtrMgAEDMTc35+DBfdSqVQsPj7bs3n0ApTJ/G3fmtaCMBa8VCgUjR47GzMyE9u09WbFiGSNHjsbCwoxOnTqzePEXaDS5tGjRnDVrfsbRsRa3bt1Crc7B2tqaW7du6veHUpmfRkjIISwtLRk1agwALVq4MmDAIIKDt9C+fftC+12pVJbp96LcgoiVlRU5OTkG03JycrCysjKYZmZmhlar5d1338Xa2hpra2tGjRrFhg0bShREZFfwlbcr9AKVtfw6XX434AXdgfu7vVWuVyL+bm8Z1RW7VqtDqVTi6FhHv7y9vQOHDx9Co9EhhECrFYVeA+h0AiHy3ycmJuDg4KifZ2VlRdWqdmi1Ov20O1/b2trpXyuVKrRadaH8xsbGoVar6du3l36aTiewtbVFo9ExZ84CVq5czuefzyItLY2ePb2ZPPldrK1tDPKq1f53LLRaHebm5qhUpv+ml18rY2lpjUaj49/DhlqtxcTEhO+/X8W+fbuxs6uGi0szdDqBRqPV74/89zri4uKIi4vF27ubwb5t2rRZkcdBp9MV+b0obVfw5RZEnJ2d+fbbbw2mRUZG0qZNG4NpjRo1AvIvcwvcOUSmJEkl84b7W0ZXL5U3nU5HcnIy9vb2AMTHx+Po6FhoOZVKZfA7kJr631C5NWs6kpDwX5PanJycB77xbG/vgK1tVbZu3XVHmink5OSg0Wi4ejWKd955D1NTUy5fjmDmzPdZv/4nxo17/Z7bVdw9BnAx1q9fS2joJX7+eTO2trYIIejXr0exeW3SpCnffps/3ryJiZL4+MRCww0/LOV2Y93LywuVSsU333yDWq1m7969hISE4OfnZ7Bc06ZNadmyJfPmzSMrK4uYmBhWr14tb6xL0hPqm2+WkZubS0REGFu3bqFv38LfdSen+vz99wFyc3OIj49n+/Yg/bz+/QexefNvREVFkpuby7JlSx74xNPVtSV2dnZ8913+79Xt27eZPv1tVq5cjkqlYs6cWaxfvxatVou9vQMKhYKqVfNHZjQzMyMrq+ghf42VmZmJqakJpqam5OTksGLFMjIyMtBoNIDhsMIdO3YhPv4GW7duQaPREBcXy1tvvcaWLRsfKA/GKrcgYmZmxsqVK9m/fz8dOnRg/vz5LFq0CCcnJ4KCggyuSFauXAmAt7c3Q4cOxdvbm9GjR5dXViVJKkdVqtjy/PMDmTZtCi+//Brduxc+4x4x4mUsLCwYONCHadMm4+Pjq5/n7e2Dn99A3nzzVZ57biA2NjaoVCpMTU1LnScTExMCAxcTFnaJIUN8GTFiKE5O9ZkyZRoKhYLZswP5+++D+Pr2ZMSIobRt68ngwUMB6NOnHxs2/MzMmR+UOv1hw0ZgaWnFwIE+vPDCIG7dSsbTswPR0VcB6NmzN0eOHCIgYCy2trYsWrSUHTv+oH//3rz66st07tyN0aPHlTr9kpDD4z7BKnPZofKW/3EZHvfGjTiee24gW7fuxs7OrtTbiYgIw86uGg4ONQFQq3Po2bML69ZtrJRjrlfY4XHj4uKYP38+b7zxBomJiWzatIl//vmnxAlKkiSVpWPHjjJjxnTS0tLIy8tj1apvqFu3Hk5O9R911ioFo4LImTNn8PPzIzQ0lIMHD5Kbm8v58+cZMWIE+/fvf8hZlCRJKt5zz71Iw4bODB/+LP37exMeHkZg4GKjb2JLD8ao6qyXXnqJ7t2789prr9GmTRuCgoJwcnJi6dKl7N69my1btpRHXktEVmdV3uqcApW1/I9LddbDcr/qnCddhazOunjxIn379i00fdCgQURFRZU4UUmSJOnJYFQQqVatGteuXSs0/dy5c9SoUaPMMyVJ0oN5QtvLSA/oYXwujAoiL730EjNmzGDbtm1Afv9Wq1ev5pNPPmHYsGFlnilJkkov/ylszaPOhlQBabUalEpVmW7TqCfWx40bh7W1NQsWLCA7O5sJEyZgb29PQEAAY8aMKdMMSZL0YCwtbUhLu02VKtVRKOS4c1I+IXSkp9/G0rLk9z3uxagb63FxcdSqVQulUklWVhZarZYqVaqg1Wq5dOkSLVu2LNNMlQV5Y73y3lguUFnLL4QgO/s2aWkZQOWr1lIqlfr+wyqj4suvwMzMgmrVHIpsufZQ+87q1asXhw4donr16gYdJsbFxfHSSy9x5syZEicsSdLDoVAoaNCgQaUMoFB5Tx4KlHf5iw0iv/32G+vXrwfyz2zGjRuHSmVYl5aUlETdunUfbg4lSZKkCqvYIOLr66sfbP78+fN06NABa2trg2Wsra3x8Sn78YMlSZKkx0OxQcTKyoo333wTgLp16+Ln54eZmVm5ZUySJEmq+Iy6J/LMM88QGhpKeHi4/oaNEAK1Ws2FCxeYNWvWQ82kJEmSVDEZFURWrFjBokWLsLKyIjs7mypVqpCenn/jpnv37g81g5IkSVLFZVQj8p9//pmpU6dy6tQpHBwc+P333zl48CBubm4VsnmvJEmSVD6MCiJJSUn06dMHgObNm3P69GkcHBx49913CQoKus/akiRJ0pPK6L6zUlLyxyxu2LAhYWFhADg6OpKYmPjwcidJkiRVaEYFkV69ejFjxgwuXbqEl5cXv//+O8eOHeOHH36gTp06RicWFhbGsGHDcHd3p2/fvhw4cKDI5Q4dOkSLFi1o06aN/m/ZsmVGpyNJkiSVD6NurE+bNo05c+YQFhbG4MGD2blzJ6NGjcLa2pqFCxcalZBarSYgIICRI0eyZs0aDh48yKRJk9i6dWuhBxYvXLiAr68vX3zxRclLJEmSJJUbo4KIpaWlQTPeOXPmMG3aNGxsbDAxMWoThISEkJOTw5gxY1AoFPTq1QtPT0+Cg4Px9/c3WPb8+fO4ubmVoBiSJEnSo2BcBAD27dtHREQEarW60LyChxLv5cqVKzg7Oxt0/NWoUSMiIiIKLXvx4kXS09NZu3YtAP369WPSpEnyYUdJkqQKxqggMmPGDDZu3EjDhg0xNzc3mKdQKIwKIllZWVhYWBhMs7CwIDs722CaWq2mdu3a9OnTh8GDB5OYmMjEiRPR6XRMnz7dmOwCGPRG6eBQxej1njSVuexQucsvy155lWf5jQoif/zxB/Pnz8fX17fUCVlZWZGTk2MwLScnx6BXYAAzMzPWrFmjf9+gQQP8/f0JDAwsURCRXcHL3kwrc/ll2Stn2aH05X+oY6zb2NjQpEmTEm/8Ts7OzoXGY4+MjKRx48YG027cuMG8efPQaP4bmS03N1dWZUmSJFVARgWRiRMnMnv2bMLCwsjMzCQ7O9vgzxheXl6oVCq++eYb1Go1e/fuJSQkBD8/P4PlqlatyubNm/n666/RaDRERUWxfPlyhg4dWvLSSZIkSQ+VUdVZ1apV48KFCwwePLjI+ZcuXbrvNszMzFi5ciUzZ87k66+/xtHRkUWLFuHk5ERQUBAff/wxp0+fxsrKipUrVzJ37ly+//57LC0tef755xk7dmzJSiZJkiQ9dEYNj9urVy9atWrFkCFDCt0cB/D09HwomXsQ8p6IrBuuzOWXZa+cZYfyvydi1JXIrVu3ePvtt3FycipxApIkSdKTy6h7Iv369WP37t0POy+SJEnSY8aoKxFbW1sWLVrEli1bcHJyKjTW+v/+97+HkjlJkiSpYjMqiKSmphZqRSVJkiRJRgWROXPmPOx8SJIkSY+hYoPITz/9xNChQzE3N+enn34qdgMKhYLhw4c/lMxJkiRJFVuxQWTVqlX4+vpibm7OqlWrit2ADCKSJEmVV7FBZO/evUW+vpsRj5lIkiRJTyijRzYsGB73TgkJCXTq1KnMMyVJkiQ9Hoq9EtmzZw8nT54EIDY2liVLlhR6Wv3atWsPN3eSJElShVZsEGnWrBk//vijvrrq4sWLmJqa6ucrFAqsrKyYO3fuw8+lJEmSVCEVG0Tq1q3L6tWrAXjvvff44IMPsLEpeb8qkiRJ0pPLqHsic+bM4datW2RkZABw5MgRPvnkE7Zs2fJQMydJkiRVbEYFkS1bttC3b1/OnTtHREQEr7/+OqGhoXz++eesWLHiYedRkiRJqqCMCiIrVqzgo48+omPHjmzatImGDRvy888/s3DhQtavX/+w8yhJkiRVUEYFkZiYGJ5++mkADhw4QI8ePQBo2LAht27demiZkyRJkio2o/rOql27NufPn+fmzZtERkbqA8revXupX7/+w8yfJEmSVIEZdSUybtw4Jk6cyPDhw/Hw8MDd3Z0vv/ySefPmMX78eKMTCwsLY9iwYbi7u9O3b18OHDhwz+V1Oh0jRoxg1qxZRqfxILJzNUTdSCNHrSmX9CRJkh53Rl2JvPDCC7Rq1Yq4uDi6du0KQNeuXenTpw9NmzY1KiG1Wk1AQAAjR45kzZo1HDx4kEmTJrF161bq1q1b5DrLly/n5MmTuLi4GFmcB/PF+n+IupGGpbkJ3d3q0KttPWpULTwcsCRJkpTPqCAC0KJFC8zNzTl48CCdO3emRo0a1KtXz+iEQkJCyMnJYcyYMSgUCnr16oWnpyfBwcH4+/sXWv706dNs27aN3r17G53Gg4q/lQnkX5H8eewaO49fp21TB/q0d8K5btVyy4dUvnRCxxfH53Ik7hCWJpbYmNlQxcwWa1MbqphVwca0yr//bbAxs8HGzBabf+fVsa6Lqcr0/olI0hPKqCCSkZHBlClTOHjwIEqlkh07djB79mxu3LjBN998g6Oj4323ceXKFZydnVEoFPppjRo1IiIiosj0pk+fzvz581m3bl0JivNg3nymFat3hJFwOxsAnRAcD03keGgiznVs6eNZHw8Xe1RKo2oBH5krcalsOxxNZq6G1o2q49XCEfuqlo86WxXWl6cW8cWJ0vW8UNPKkd8Hb8fZrkkZ50qSHg9G/RrOnTuX3NxcDhw4gLm5OQAfffQRVlZWfP7550YllJWVVajvLQsLC7Kzswst+/HHHzNo0CBat25t1LbLSvOnqjP7tQ5MGNqa5g2qGcy7EpfG8i3nmf71Uf4MuUZWTsW7bxIdn87/fj3D7NUn+edyEhHXU9h4IJKpy48w96dTHPgnlsycvEedzQrlSNwh5hz7tNTrJ2YlMGHvG2h12jLMlSQ9Poy6Etm/fz8rVqwwuOKoW7cuM2bMYPTo0UYlZGVlRU5OjsG0nJwcrKysDKZt2bKFuLg4AgMDjdpucWrU+K+LFgeHKiVat3dNW3p3bEhUXCpbDlzh4OlYNFodAMlpOWzYd5ngw1F4ezZgQJdG1La3fqC8Pqhr8Wms2xHGobNxxS4Tfj2F8Osp/LQrAk9XR572cKJd85qYmqjKMafl717HPjEzEf/dY9GJ/GPbyakT0ztPJ12dTlpuGum56UW/VqeTnpvOxZsX0Qotx+ND2HB1NRO8JpRXsYxS0s/9k6Qylx3Kt/xGBZGiriIgv/WUTqczKiFnZ2e+/fZbg2mRkZG0adPGYFpwcDChoaF4eXkB6ANPbGxsiZ6OT07OQKcTODhU4ebNdKPXu5ONqZIR3k0Y0KE+e0/Fsu90LBnZ+Wfy2blagv+KZOtfkXi4OPCidxOq25bvTfiEW1n8fiiKkAsJ3DmqiwLwbOGIV8va7D1xjQtRtygY9kWj1XH47A0On72BtYUJ7ZvVpINrLZrUq2pQ1VjWdDrBmStJ7DsdS2qGmj7tnejYshbKh5jmvY69Vqflha3DuJFxA4AaFjX4qscq6tgU3cijKIHHPtdXg723+z061niap6o2fPCMl4EH+dw/7ipz2aH05VcqFQYn38ZSCCNGlZo4cSIqlYrAwEDat29PUFAQ1tbWTJo0CTs7O5YsWXLfhNRqNb179+all15izJgx/P3330yZMoXg4GCcnJyKXW/69OlYWVkxY8aMEhWsLILI3dR5Wo5eTGDn8evEJWUazLOxNCVgkCvNn6peJmndS1JKNkGHr3L4XDy6uw5fWxcHBnVtSD0HG33ZUzNyCbmUyJEL8UTHF70v7Kta0MG1Fl4tHKlTw6rMAkp2roa/zt5gz8nr3EwxvBJtVMeW4d4uNKpjWyZp3e1ex/6L43MJPP5fVez6/hvpWb9kjTjUWjW9f+3GpVsXAeha72l+G/D7Qw3GxqrMP6SVuexQQYNIUlISb775JuHh4WRnZ1O3bl0SExNp3Lgxy5cvN+rGOkB4eDgzZ84kNDQUR0dHpk6dSo8ePQgKCuLjjz/m9OnThdapSEGkgBCCC1G32Hn8Ouej/ntiX6GAod2d6etV/6H8kNxOz2Xrkasc/CcOrc7wsLV2rsHgrg15qtZ/P8hFlT0uKZOjF+M5cj6B5DTDH/UCNatZ4t7YnjZN7Glcr2qpGhIk3M5iz4kY/jp3g1z1ve8XdGlVm2e7N6KqjXmJ07mX4o79XzEHGBo0EPHv9dvktu/wnlfJPl8FTiecpN+mXvoqsYVPf8mIFsZV8T5MlfmHtDKXHSpoEClw9OhRrly5gkajwdnZmc6dO1eIs66iPMwgcqewa7dZ/vsF0jLV+mltXRwY69ccS3OjW1DfU1qmmu1Ho9l76r97MwWaN6jGM90a0biIJsj3KrtOCC7HpHLkQjzHLyWSlVt0QwEbS1PcnGvg3sSBlg2rY25W/D0UIQQXo2+z+/h1zl5J5u4PlrWFCd3c6qBQKNh5/Boa7X9LWJipGNi5Id7t6mGiKpvWb0WVPyErgZ6/dOZmdiIAnep04beBQZgoS3+sPjn8Ecv++R8AVcxs+WtYSImqxcpankaLlY0lt25loFAoUCoUKBSgVChQKvPHAip4X1G/v6UhhOBc5C2S0nNxrW+HY3Wr+6/0BKrQQeRxUl5BBPKvEJb/fp7LMan6abWqWzF+SCvqPsBN98ycPP4MucbuEzHk5hmezTepV5Vnujai2V2tyO5kbNnzNDrOXkkm5GI85yJvFUqrgIlKSYunqtGmiT3uje31Vw65eVqOXIhnz4kYYu+q5gOoXcOK3u2c6OhaSx+EEm9nsX7PZf65nGSwbK3qVrzo3YRWjWrcN9/3c3f5tTotzwUP4u/YgwDYWzqw7/lDOFrXeqB0sjXZ9PilE5GpVwDo06Ava3x/Kdcf6LRMNf9cTuJU+E0uXr1d6GTjXpQKUCihvqM1T7s70a5pzTI7ASoPOl1+U/xtR6KJuZk/XIWJSkE/rwb4dWyAmemT3XjkbjKIlJHyDCKQf8P6l72X2XMyRj/N3FTFWL/mtG9Ws0Tbys7VsOvEdXYcu072XVcIDWtX4ZmujXBtWP2+P1IlKXtabiqLTy0gS51D9byWaFJqEx9nRkZW0QFFQf49DaeaNhwPTSSziCbPrZ1r0LudEy2eqlZsXs9HJfPz7ghuJGcZTHdvbM8LvRrjWK3kZ5NCCPI0OurWsTMo/9xjn7HwROC/+Vfw68Df6Vbv6RJvv8Dt9Fwys/NwrG7FycSjDNzSVz/vK++VDHV5odTbNsbNlGxOhd/kdPhNImJTKatvsrmpinZNHejcqjYu9e0eauOHB5Gn0XHkQjzbj0aTeLvwowIANe0seamPS5mclFREOp0gOiGd0Gu3CY1O4XpiOq0aOzC6jwtKZcmOmwwidynvIFLgyPl4fvwzFLXmvzNBH08nhj7tfN97C7l5WvadimX70Wh9K7AC9RxseKZrQ9yb2Bt9hluSsvvvGsumiN8MJwoF1XVNaKz0xl7tgTLX/r7bMTdT0aVlbXq1q0ctI6sTNFode0/G8PuhKLJz/wtaJioFvds70b/jUwZnxkII0rPySErNISk1+9//+a+T/32dp9HR0rkGQ7s1or5jFfZd28OwrUP090HeaTedqZ7vG5W/u91Oz+WXvREcu5RfJaZUKHCsbsktrnA2Yz/pymsoLVP5c3gQtWxKdgJxL0IIridmcCr8JqfCk/Rn3UWpYmWGTqdDiPyqSyHy19cJgVanQycECiMeE7OvakGnlrXo3Ko2DnYV44HVXLWWA2fi2HHsGrfTcw3mmZkqqVXDmmt3NSBp19SBF71dqFalbO+7lTedEMQkZhAafZvQaymEXU8pdKIJMGusJ/VqliwgyCByl0cVRACuJ2awbNM5ElP+Oztq6mSH/+CWVLU2K7R8nkbHwTNxbD1yldQMtcG8WtWtGNy1Ie2a1SzxGaGxZY9Jv077ta3RinvfALfSOeKY54mjxpPq2hYo+a+awL6qBd5t69GldR2sLEpXFZKaqWbTgSv8ffaGYZNl0xysaiRRy6wx2VkqklKzUecZV12jUIBnSzuW3hhFgvoakN+KakP/zaiUJavm0Op07DkRw5a/o8i5T2OB/MS1ODlUpa6DNXXtranrYEMde2vMi6peucfXMOF29r+B4yZJqUU3hlCQX8XZxsWBNi4OuDapaXDsr6VFs/rC9/x06UeSc5L/TTN/TQVKzIQtdfO64ZTXiyq6onvmbupkR5fWtWnb1AELs+KPcY5aQ1JKDjdTsklMyeZmSjY3U/KDvFKhwKmmDU6ONtR3rEL9mjZUsSr8nShKVk4ee07GsOtETKGTLCtzE7zb1aNX23o85VSdjbvD+G3/FYN7feZmKgZ3yb/vVt69Tuh0grQsNUqFAhOVElMTBSqV8r7faZ0QxN3M5NK124RG3yb8ekqRV/13atfckVf9mpX4GbCHHkRu3brFuXPnSE5ORqlU4uDgQIsWLahWrfg6+UfpUQYRyP/Arwy+yJkryfppdjZmvPFMK/1NcK1Ox+Fz8QQdiiI5zfCMyr6qBYO6NKSDq2OpP/DGln3WkRksPb0YgObVXWlXqz2RKVeISo0kLjO2yHVMhQ0189piratLquoyM/u9zrMuz5Uqn3eLupHGul3hXIlLK5PtAagVaYSaryWn6ln2vvA3Na1KdoUQfj2FtTvDiLlpeM+nuq05t+46duXFRKWgxVPV8XBxwL2xPbZ3nKA4OFQhITGV/df38N25leyK3qG/CitgqjRlgPMgXm75GgtOzGX/9b0goHuN5xhUfSohFxOK/MEyN1XRrpkDHk0cyMjJyw8Q+mCRTVpWyXpFqFbFnPo1bXByrEIDx/z/DlUt9FfcqZlqdh2/zt5TMYWCt621GT6eTjztXld/tVrwuU/LVLNh32UOn483WKeegw2j+jYtsjFKWdDqdNxIziI6Pp3o+HSuJqRzPSGjyHuNKqUCExMlpiolJqr8AFPwZ2qi4GZKTqGAebeqNmY0r1+NZg2q0ay+HS2a1CQpqfir1OI8lCCi0WjYunUrq1ev5tKlS5iammJra4tOpyM1Nf8mcuvWrRk+fDh+fn4oK1CfUo86iED+WcS2w1fZ8leU/uurUioY1qsJ1hYm/P53lL6frgJ2NmYM6PQUXd3qPHArJWPKnpmXifvq5qTmpgCwxvcXfJ7qp5+flZdFdNpVolIjiUzNDyxR//6Pzbjj/o/KnE2DttK+ltcD5bmAWpvH2HUz0MW1xUIUfvZGp8yhWhVTGjrUxN7OEvuqljhUtaBGVQvsq1qQmqnmtwORnA6/abCefXUFr/t5GN2hZlqmml/3X+bQOcMfoto1rBjRpynNG1QjR63hRnIWMTcziDu8nOAAACAASURBVEvKZHfYcTLSTLAU96/+KylLcxWtne3xcMlvLXf3DXAhBDezb/JH7GaWhiwjOu1qoW3UtanHaNexDG8+Sh9Mw26F8vQvHfVXo9/5rMWngR//XE7m0LkbnItMLrN7LsawNFfhVLMKdjZmnI5IIk9jeOVZw9YC3w716dyqdqEb53d/7sOu3Wb1jrBC9926udVm6NONsbEsfQeaGq2OuKRMohPyA0b0vwFDrTG+YUNJ2Vqb0ay+HU3r22FvryFJRBJ++xLht8O4nn6Nvi59GNv0jRJvt8yDyLFjx/jkk09wcHCgd+/edO3aFScnJ/3ZgRCCiIgIjh8/TnBwMLdv32bWrFn6J80ftfOfB6K+nYKpmQl5j3h8kMwcDXFJmeh0xX8LVSoFNWwtqGpjTgnvhxXLmLLHZ97QtyqyUFnQxrEdxiafp9NwPuks2Zr8L6ep0pRWDu5YqB683jkqNZIbmXGAEnNdVSxMrMjWpqNFjU6RhyD/x87CxJI61nVwsHJEpTAMuumaVEITorHW1UIlDKtMqtqYYW9niUkxO1sAqRm53EzJMThuCiXY21pSzda82P2k0Wk4nXgKjU6LSphjZ+pAVVN71Hla1HnaQk2f70epUGBtYYKNlSmmZpCnVaPW5pKr+/e/Nhd1wTStGl0x1ZJ25tWoZV2bahbVUBSR+//2OZirLGhT0wPlv/tUo82vjknNyL1/VaICTE2UmJkoMTVRYWqi1P8JAblqDTlqLbl5WnLV2hIFJzNTJTVsLahibVbs/i/qcy+AW2k5JKfmGKSnUilwsLPE9o7tCfKrn7Q6UeT/rLxsMtXZCJ0Khc4UhHHfGIUyv9pR8G/tpZHlVqkUmJspUZlo0ClzyBEZZGmyyM7LRiuK/n671/TAyqRkjVLMqtnR8v2pJVoH7hFE3njjDSZNmmT0WB4XL15kyZIlfP311yXOxMNQkYII5N/3iE3KLPTgnVKpoLqtOdWqmJd5K5j7lV0ApxNPkqPJvxpqWLURta3rlCiNHG0OZ2+eQaPLv+S2MrGmlUNrVIrSN6u8kRlHVGqk/n29Kk7Ur9KALE0WcRmx3My6icDwh8xEaUot69rUtq6NqdKUXG0uZ27+82++lFRT1kOVZ2vwA6JUKrCvaoFdFcOAkKPWknArq1DViY2VKTWrWWJqxBXirZxbhP77JDtAixqu2JmXrOo3R5tDYlYC6er0fwNFrv6hRmOplCbUtHKkllUtLE3+397dx9V4/38Af51OpY6ESO6iVVRui1SWjRg1aU1uSjSZzabfGDb23RcLX4Yxdmcmct+GyCr3o4UNLTS560Y37opRoTqdzt3n98dZF0ednHPU6dR5Px+PHs75nOvq+ry7jvO+zvW5q71hXCqX4uI/5yGVK94zXSzt0Nmi+nIPFWIZnpSLUSmRwZhflSye/hjzjdS+EGEAxBI5KiVSVIplXHKRyZQ/lpqZ8tGmpRkszE1e+Ltre99LpHL8U1JR7RaRibERlzxqu9hTl4wngQwVkPJEkPIqIEUFGK+mOvFgBD74PD6MeMaKf8GHEY8PIxhBBinKZY9VXhjUxMzYHH2t3apdVL1InSeRxk4fbmc9TyyRYedvWfgjvRDNTPgYPsAWfh62EJjVz3oUL4r9+M2jCD2oaMdoYWqJS+9ch4Wp5hO3nSs8i7HxARDLFZ0ChnUZjh0jd2s1gO/4zaOYdCiY+7AMdAjChhGbuStiALhffg+bLm/A1qvR3G24KmZ8M4x3CsXVosu4cD8VANCheUckjf8TskoBdp/IRlq28tiUTtbNMfGN7rC1scC+k7k4mXZX6SKxXStzhA7vjj4OmnUTnXYsHL/eiFMcw6IzToWcQwvT2qd4EcvEOJJ3EDuubcWpO8nV2jHUYW5sjr7t+yKkWxjedhwDgYn6V6Sbr2zEf059AkBxQXBu4kW0b95B4zrU5NI/aaiQieDR3lPpfD6PMYaS0krc+qcMD0oq0LmdBZy7tFKrV2KxqAjlxsVowzrVGnda1gP8fDyrWlukNoS8f/CYn4PH/Fw8NsrBY34OxEaPX7xjHbA0bQknK2c4tXZGdysndG/tjLf6+KHskeYXznV+Oys+Ph6BgYG17vzo0SMsXrwYa9eu1fjA9U0fk0iVosciNDc3rrWHS114UezjEgJx8s7vAIAP+36EJd7qTetfk90ZP2NG0tPFxd7v/SGWvabZTMxXH17BqP0jUC5RNAr2t3FHXOBBlVfQZZIy/HxtOzak/4jbpbdq3IbP4+PXtw/Ds4MXV3Y5twg//5ZVrT3KzJSv9O3DmG8E/4FdMdKri1azHT+seIhBv7ijWKSYGie851R8Nbjm/ys3SrKx8/o27Mn8GQ8rHta4DQA0N7FAJ4tO6NC8IzpadEIHi47o2LwTOlp0REeLzujYvCNaNmuFdu0stXrfS+VSDNsziJsPbLzTBPwwTP2JT1X56dIP+OJPRZfqrpZ2CO/5HkJdJqG12cvPNccYw1/3UrDlykYcyImHWC5GcxML+NsHIKjbOLzeeUiNFzSVYhkSzuTh2F+3q00jZN7MGObNjFAqK8L9ypuoRCkkvHJIeGWQ8srRt30PuHd1gcBSCJmREBVSIYT//lRIKlAhFaJCWoEKaQWEknLFv9Kqf58pk5SrdaHQullrOFm5oHtrZzj9myycrJxhI2hfLbnqzWBDV1dXDBkyBIsXL0bLltUbIY8fP47IyEjIZDKcO3dO4wPXN31OIrpSW+wZxdfx+i5F+5URzwh/TbyELpZdX+p4X55bgm8uruaer3j9a7zb63219r1ffg9++4ZyjfW2Lbrg8JgktXpRSeVSHMiJx7q/v8OlB8rzry0cuAQz3GZV20cileO387eR+Gd+jb1metu3wcTh3dBOi8GOz4rLjsWHv03lnv8aeAivdhoEQDHSPTHnV8Rc346zBX9W25cHHoZ2eQMhzhPhZOWCThadXvhNpsrLvO9P3UnG2IS3uOeHx5xAf5sBWv0uAIjN3IX/OzGtWrkZ3wyju43F1N7T0MfaVePfWyYpw76sPdhyZROuFV1RuV1b87YIdAxCULdxcLfxqPahW1YhQdFjEczNjCFoZgwjvhRbrm7EtxdX49Fz33SHd/XFAq/FcGnTQ+P61oQxhkpZpSLBSCoUSUgqhFCiSEjN+M3QrbUTrM2t62V82LPqPInk5OTg008/RVFREZYvXw5vb28AwOPHj7FkyRJu6dpFixahTRv9Gw1KSaT2N9MnyTOx49pWAMAo+0Bs9tvx0seTMznePxaOxJxfASi+BcT4x2Jolzdq3U8oEeLtX9/E3/8mgBamljgY9BucrVw0Oj5jDGcK/sC6tG9xtvAM3ukbhsgBy2u9dVJSWok9v99AyrX7ABRddicM645+3dUf1PmiOr1zOARH8w8DAOwsX8FPw6MRm7ULe7P2VLsdBwAdm3dCqEsYQl3C0LmF6hmua/Oy7/vwwxNxKC8RgOIb4cGg47X+HVU5cfMYwg6HcO0sPPBqvPLubzMA7/Z6H285jkazF3TMyCi+jq1XNmFP5i6USarHaC2wxgPhgxr2VLTzjOk2FkHdxsPJylnpNTmTY2/WbqxIWYo7ZbeVXnNr1w+RA5dyFwD6TG+SCABIJBKsXbsW27ZtQ0hICNzd3bFs2TLIZDIsXLgQI0eO1PiAukJJRPWbqaiiCG7bXSCSKQauJYw+Cq8OA+vkmDUlhENBx6v9h60iZ3JMPfoODuYmAFA/8ahDk3OfV/gE90uEcHO0rnWSSW0UlhVg0C4PlIpVj3vh8/gYYfcmwnpMho/tGxoPhHzey77v8x/nYdAvA7h2rnXDojDOKUSj33HhfirGxAdA+G/vPRerntgTsB8nbv2G6MtRuPzwUrV92pq3xSSXcLzTc4pSAhXLxDiUm4gtVzfV+K1NYCxAULdxCO81FUNdBuHEtdPYlx2LX2/sw73ywhrr17NNbwR1H4cgx7HIKsnEkrNf4GrRZaVt7CxfwQKvRQhweLvRTFapV0mkyunTpzF9+nTIZDK4u7vju+++09tBhlUoiah+M31zYTW+TFkCAOhr7YZjY5Pr9D/I/fJ78N3rww1U7GJphyNjktDWvPq4if+djcT3aU/bCVa+vgZTer1XJ/XQp3Mfc207Zid/VK28q6UdJrlMRojzxJeeCPJZdRH70rOL8F3aGgBA++YdcCb0AixM1PuQySrORMD+ESipLAGguD15MOg3rpGeMYbz9//C5ssbkZjzK5esqhjxjOBrNxITnCch7Z/z2HltOzfz8rMcW3XDlF7vYbzTBLRs1gqAcuwyuQxnC/9EXFYsEnPja/zmV5O25m3xiftnCOsxBaZ89UbU6wu9SyK///47Fi1aBLlcjv79++PYsWOYPHkyZs+eDVNT/f3jUhKp+c0klonRf0cv3BcqBs9pc4WpjssP0xEQ5wuhVDHC26O9F/YFJirdqnj+g/WDvv+H/3kvr7M66NO5Z4xhypFJOJSXCFMjU/jbB2BSj3B4d3pNq9tEL1IXsZeJS+H1cz/8I1Tc6lN33ZWCsrvwjxvOtW+1MWuDA0HH4NCqW43bPxA+QMz1bdh2dbPSAFZV+Dw+RtoHYEqv9+Dd8TW1G5YrZZU4cfM3xGXH4lj+Ye6b+LMExgJ82Pf/8H9uH6vd/qRv9CaJPH78GEuXLkViYiICAwMxf/58WFpa4vjx4/jiiy/QsmVLrFy5En369NH4oLpASaTmN9PerN2IOK5o7G4nsMHFsKv1dqV1JO8QJh+ewN0DH9s9GOuGRYHH4+H0nZMIPjCau1fua/cmtvr9/NK3cZ6lb+deJpfh8kNFBwYrs/ptR6yr2HdlxGBm0nQAilkJ/piQiq6Wdiq3LxEV4639fsgsyQCg6Ca8P/AA3Gz6v/BYUrkUR/MPY/OVjTh9J7na6+2bd8A7PaZgUo/JtXY7Vif2UvETHMo9gH3Ze3Dq32NNdJmMeQM+r9NvhA1Bb5LIoEGDYGRkhMWLF8PHx0fptYcPH2LBggU4ffo03n33XXzyyScaH7i+6SqJXH14BSWVxXi146B6uaJ8Gc/HzhiD794hXHvFfzwWYI675oOLNPHj399j0Zn53PPPPRZilEMgRsa9wd1a6NW2DxJGH1H7Vom69C2J6FJdxS5ncry5byjS/rkIoPZOGEKJEOMSA5F6LwUAYGxkjJiRsfDpMkzj42aXZGHLlY344+4pdGjeEWE9psDX7k2Y8F88pkrT2EtExZAxeY23WxsjvUki8+bNw4IFC2Bpqfor3e7du7Fy5UpcvHhR4wPXN10kkdR7KXhrvx9kTAavDq/i6yHfoVtr9Ub468LzsacUnkPA/hEAFFeVae9cr/f/OIwxfHryY64nGABYm7fj7m+3b94BR8Yk1ctKgJRE6ib21Hsp8I97uv78/sCD8O70mtI2UrkU4YdDcezmEa5s/RubMKb7+DqpgyYM+bwDuk8iKi+dv/rqq1oTCAAEBwcjPj5e7YNlZmYiJCQErq6u8PPzw8mTJ2vcLi0tDcHBwejXrx8GDx6MdevWQR8H1v/49/fchHXnCs/AZ/er+Pr8Sohl4hfs2TCi0n/kHo/tHqyTKy8ej4cVr32N1zoN5sqqEojAWICdI3c36FKy5MUGtPfEmG5Pk8H8Pz6DTP50bA1jDJ8kz1RKIEu9VzRIAiG6pzKJFBQUqPXD5/OVnqsiFosxffp0+Pr6IjU1FXPnzsWsWbNw967yVONCoRAffPABxo0bhwsXLiAmJga//PILjhw5ouI3N4zCsgIcyTuoVCaWi7Hyr2V4I/Y17iu9vrhdeovrRgsA0/poPsuntkz4Joj23Q6HVo5cGQ88rB8erdUgM6J7Cwcu5ib0u1Z0BTuvb+NeW3ZuMX7J2Mk9n+k2B9P66u79RRqWynk3hg4dCh6Pp9E3AB6Ph+vXr9f4WkpKCkQiEcLDw8Hj8TBs2DB4eHggMTERH374dLoMgUCApKQkWFhYgDGGR48eQSaToVWrVhqEVf92Xt/GfQvp3bYvjI343H3jjOLrGBU3AlN6vYf5XpF60csj+nIUNx/V65196mzErbpambVGjH8sxsa/hcLyAiwdtAJvvuKv0zoQ7XW06ISZ/eZgxV9LAQDLU5Yg0GE0dmf+zHUDBoBQ5zDM94psqGqSBqAyiWRkZNTpgXJycuDg4KDUHc/e3h7Z2dnVtrWwUNyXe/XVV1FcXIyAgAC9mWIeUNz/3Xnt6ZXYDLdZCHB4G9GXN+DLlP9BKFXMh7P5ykYcyTuElYPXKK3RoWtlkjKl+n7QZ3qD1MO+pQP+DD2PCqmw3nsnkbo33XUGYq5vx+3SWygWFSP04Dicv/8X97qv3ZtYPeTbRjMoj9QNlUkkPT29TrvvCoVCmJmZKZWZmZmhoqJCxR5AcnIyCgoK8MEHH2DdunWYMWOG2sd7toHI2lrzmWlr82vGr9y6C+2at8Nkz1CY8k0x3+YzTHIPwfSD03H4hmKai4Lyuwg7FIxxPcbhuze/Q3sL3XYftLZugd1/bcMTsWJW0W5W3RDiPqYBe5LV7bl4kbo+941J3cfeAl/7rsb4vYq2jmcTiLetN+JC92o0Y3B9MuTzDug2fpVJJCQkBJMnT8asWbPQrNnLLzIkEAggEikP7hGJRBAIVL/pmjVrhldeeQXh4eH45ZdfNEoi9dk769sz33OPJziF4XFxJQDFlNICWGHr8F3Y/8peLPjjM25G1thrsTiW8xsiB/4PE13eqfFqTSaX4R/hfdwuvY27Zbe5f++U3oZIKsIIOz+E93rvhXMLValaInXNmacjwt/t+QGKHpbXslfTYci9dOor9sHWvhjY0Vtp6hFnKxdsGf4zyh/JUI6G/3sb8nkHdN87S2US+fLLL7Fq1SocP34cS5cufenbSQ4ODti0aZNSWW5uLtzc3JTKcnJyEBERgYSEBC55SSSSF/YU05XcxzmKtaihaBwO6xFebRsej4egbuMwxHYoFp1ZgF0ZMQCAx5WPMCd5BvZm7UZQt3EoKL+LO6WKJHGn7A4Ky+5CIle9nvLpuycRfTkKCwcuwSj7t9S6bXD85lFugSdL05YIdg7VImpCFHg8HpYOWokRsYMhYzJ0trDF7lH70cpMv6dBIvVH5T2Nt99+G0eOHIG3tzfeffddLFy4EGVlmi/+XsXT0xN8Ph9RUVEQi8VISkpCSkoK/P2VG1e7dlVMR/7NN99AIpEgKysLW7Zswbhx47Q+dl3acXUr9/iNriNqnT7dyqwNvhu6HrEB8UqjfM8U/IFPT36MNee/wp7MX3Cm4A/cepJfawKpkv8kD1OPhuGtX/2Qdv/CC7ffkL6eezypx+Q6H9BHDE/vtn2wO2A/5rjPw4GgY+hgodlqmKRpUWsCxmvXrmH58uW4efMmwsPDq93emjhxoloHy8rKwqJFi5CRkQEbGxvMmzcPPj4+SEhIQGRkJNLSFCOpb968iSVLluDSpUuwsrLClClTMGHCBI0Cq4/bWSKpCK7bnblFhmJG7sFwOz+19hVKhFh9fgXWPzO2pCZtzNqgUwtbdLawRecWndG5hS06WdjiduktfHNhVbX1DcZ0G4/5XpE1Thl+T56PPj8p2rWMeEZInZQO2xZd1A230TPk2xoUu2HGDujRiPVnicViREdHY926dWjVqpXSxIs8Hg8nTpzQ+MD1rT6SyLOL69i26IK/Jl7SeK6nyw8uYevVzRBJK/5NEl3QyaIzbFt0QUeLTmhu0lzlviWiYqw5/xWir0Rxc04BisV9Puz7EWb2m620vO3nZ2cjOi0aABS9x3y3a1TXxs6QP0wodsOMHdDDJHL06FGsWLECQqEQs2fPRkhI3c/4Wh/qI4n4xw3nBhH+1/MLzOr/aZ38Xk3lPrqBJWcjuUWDqrQ1t8Z/PBYg1CUMjyofwW27Cypligb/A6N/g0cH/ekmrQuG/GFCsRtm7IAeNazfvn0bS5YswenTpzFy5Ej897//Rdu2TWOCMm1cfXiFSyAmRiYIdXmnwepi38oRW9+MwdmCP/HFn//lloR9WPEAn578GJsu/wQXqx5cAnFr1w8D2ns0WH0JIU2XyiTi7++Pdu3aYePGjXjttddUbWYwtl2N5h772weotfZ3fRvY0RtHx/6OvVm78eW5JdwiUBnF15FR/HTmgGl9ImgAGCGkXqjsnTV58mQcPHiQEggUi/PEZu3mnk/uObUBa6PMiGeE8U4TcCb0Av7jsQACY+U2lfbNOyDA4e0Gqh0hpKlTmUTs7e01GmTIGMPevXvrpFL6Zl92LMoliu7N3Vs74dWOgxq4RtUJTASY4z4PKZP+RliPcG5E+qz+nza65T0JIY2HyiRy7do1BAQEICYmBv/8U31t4ypFRUXYunUr/P3963y+LX3AGMPWK09vZU3u+a5e3xqyEdjg6yHfIXVSOs6/fx5TetbNeuWEEFITlW0i8+fPx5UrV7Bu3TosX74cDg4OcHR0ROvWrcEYQ0lJCTIyMpCfn4/XX38dK1as0Nulcl/G+ft/4WrRZQCAubE5xjtpNl6lodi26GLwvVQIIfVPZRIBgF69emH9+vW4c+cOTp06hcuXL+PWrVswMjJC27ZtERYWhsGDB6Njx6Y7YnXb1c3c49GOY9GymX5NSU8IIQ2p1iRSpXPnzggNNbw5l4pFRYi/Ecc9D++lPw3qhBCiD1QmkRs3bqj9SxwdHV+8USO0K+NnbqyFq7UbXNv1a+AaEUKIflGZREaNGlVtZcOqBuXny1StZtiYyZlcaWxIeC9qoCaEkOepTCLPz4fFGENAQACioqKadBtIldN3TipNoR7oGNTANSKEEP2jMol06tSpxvL27durfK0p2frMt5Bgpwm1ToxICCGGqqHWSNVrhWUFOJJ3kHuuTyPUCSFEn1ASqUHM9e3cmh/eHV9DdyunBq4RIYToJ0oiz5HKpdhxbSv3fHLPdxuuMoQQoudUtonExMRUK5PJZIiPj0fr1srrKau7smFjcCz/CArLCwAA1ubtMNI+oIFrRAgh+ktlEomOjq5W1rZtW8TFxSmV8Xg8tZNIZmYmIiMjkZGRgfbt2+Pzzz/H4MGDq22Xn5+PpUuXIj09HaampvD19cW8efM0mhBSW1uvbuIeT3R5hyYvJISQWqhMIklJSXV6ILFYjOnTpyMsLAw7duzAqVOnMGvWLBw4cKBab6+IiAj4+Pjgxx9/xKNHjxAREYFvv/0W8+bNq9M6PS/vcS6Sbyvi5oGHST0m1+vxCCGksdNZm0hKSgpEIhHCw8NhYmKCYcOGwcPDA4mJyku8FhcXo0OHDoiIiICpqSnatWuHwMBAXLx4sd7ruP3qFu7xG11HoItl13o/JiGENGZqzZ1VF3JycuDg4KA0jbq9vT2ys7OVtrOyslK6lcYYw4kTJ+Ds7FzvddyT+Qv3OJy69RJCyAvpLIkIhUKYmZkplZmZmaGiokLlPnK5HMuWLUN+fj5WrVql0fGeXXDe2rqFWvvIIAUAdG/THcH9g8A34mt0TH2kbuxNlSHHT7EbLl3Gr7MkIhAIIBKJlMpEIhEEAkGN25eVlWHu3LnIz8/Hjh07YG1trdHxiorKIJczjdbUiBkZi9/yjyDEeRKKi4QaHU8fGfp6IoYcP8VumLED2sdvZMRTuvhWez+N99CSg4MD8vLylMpyc3NrnAH4/v37CA4OhkQiwZ49e2Bra6uTOva3GYD/eC6EXctXdHI8Qghp7HSWRDw9PcHn8xEVFQWxWIykpCSkpKTA399faTuxWIz33nsPTk5O2LBhA1q0MOyvpYQQos90lkRMTU2xceNGJCcnw8vLC6tWrcLatWtha2uLhIQEuLm5AQCSk5ORlZWFEydOwN3dHW5ubnBzc8PYsWN1VVVCCCFq4rFnFwdpQrRpE2lqDDl2wLDjp9gNM3agCbeJEEIIaXooiRBCCNEaJRFCCCFaoyRCCCFEa5RECCGEaI2SCCGEEK1REiGEEKI1SiKEEEK0RkmEEEKI1iiJEEII0RolEUIIIVqjJEIIIURrlEQIIYRojZIIIYQQrVESIYQQojVKIoQQQrRGSYQQQojWdJpEMjMzERISAldXV/j5+eHkyZO1bn/79m0MGDAAxcXFOqohIYQQTegsiYjFYkyfPh2+vr5ITU3F3LlzMWvWLNy9e7fG7c+ePYuJEyfiyZMnuqoiIYQQDeksiaSkpEAkEiE8PBwmJiYYNmwYPDw8kJiYWG3b2NhYREZG4qOPPtJV9QghhGjBWFcHysnJgYODA3g8Hldmb2+P7OzsatsOGTIEQUFBKCws1FX1CCGEaEFnSUQoFMLMzEypzMzMDBUVFdW2tba2funjtWlj8czva/HSv6+xMuTYAcOOn2I3XLqMX2dJRCAQQCQSKZWJRCIIBIJ6OV5RURnkcgZr6xZ48KC0Xo6h7ww5dsCw46fYDTN2QPv4jYx4Shffau+n8R5acnBwQF5enlJZbm4uHB0ddVUFQgghdUxnScTT0xN8Ph9RUVEQi8VISkpCSkoK/P39dVUFQgghdUxnScTU1BQbN25EcnIyvLy8sGrVKqxduxa2trZISEiAm5ubrqpCCCGkjvAYY6yhK1EfqE2E7g0bcvwUu2HGDjThNhFCCCFNDyURQgghWqMkQgghRGuURAghhGiNkgghhBCtURIhhBCiNUoihBBCtEZJhBBCiNYoiRBCCNEaJRFCCCFaoyRCCCFEa5RECCGEaI2SCCGEEK1REiGEEKI1SiKEEEK0RkmEEEKI1iiJEEII0ZpOk0hmZiZCQkLg6uoKPz8/nDx5ssbtCgsLMXXqVLi5uWHo0KHYt2+fLqtJCCFETTpLImKxGNOnT4evry9SU1Mxd+5czJo1C3fv3q227axZs+Dg4ICUlBSsXr0aK1aswN9//62rqhJCCFGTzpJISkoKRCIRwsPDYWJigmHDhsHDwwOJiYlK2+Xm5uLy5cuYOXMmTE1N0a9fPwQE2lQwWwAACU9JREFUBNC3EUII0UPGujpQTk4OHBwcwOPxuDJ7e3tkZ2crbZebm4v27dvDwsJCabsDBw5odDwjI16Njw2NIccOGHb8FLvh0iZ+bf9mOksiQqEQZmZmSmVmZmaoqKhQKisvL1druxdp3bo597hNG4tatmzaDDl2wLDjp9gNly7j19ntLIFAAJFIpFQmEokgEAi02o4QQkjD01kScXBwQF5enlJZbm4uHB0dq213//59lJeX17odIYSQhqezJOLp6Qk+n4+oqCiIxWIkJSUhJSUF/v7+StvZ29vDxcUFq1evRmVlJdLS0pCYmIjAwEBdVZUQQoiaeIwxpquDZWVlYdGiRcjIyICNjQ3mzZsHHx8fJCQkIDIyEmlpaQAU40QiIyNx8eJFWFpaIiIiAmPHjtVVNQkhhKhJp0mEEEJI00LTnhBCCNEaJRFCCCFaoyRCCCFEa5RECCGEaK3JJhF1ZwxubNLT0+Hh4cE9F4vFiIyMhKenJzw9PbF8+XLIZDLu9SNHjmDEiBFwdXXFpEmTkJ+fz73WWGZLTk9PR2hoKNzd3TFkyBB8//33YIwZROwA8PvvvyMgIABubm4YPnw4du3aBcAwzn0VoVAIPz8/REdHAwBKS0vx8ccfw93dHd7e3oiKilLaPiYmBoMHD4abmxs+/PBDPHz4kHutMX027NmzBz179oSbmxv3s3//fv0696wJqqysZD4+Pmzz5s1MLBaz48ePM1dXV3bnzp2GrtpLOXjwIOvfvz9zdXXlylavXs1CQ0NZSUkJKywsZIGBgWzDhg2MMcaysrKYq6srS01NZZWVlWzNmjVs5MiRTCaTMcYYGz9+PFu2bBmrrKxkFy5cYO7u7iwtLa1BYlOlrKyMeXl5sR07djCpVMry8vLY0KFD2c6dO5t87IwxdvfuXdarVy92+vRpxhhj169fZ71792aXLl0yiPirfP7558zZ2Zlt2rSJMcbY7Nmz2cyZM5lQKGQ3btxgQ4YMYQcPHmSMMXbq1Ck2cOBAlpmZyYRCIfvss8/Y1KlTGWON77Nh4cKFbM2aNdXK9encN8kkUvUmksvlXNm0adPY+vXrG7BWL2fNmjVs9OjRLDo6WimJeHt7s+TkZO754cOH2YgRIxhjijfazJkzudekUinr378/u3DhAsvJyWEuLi6stLSUe33x4sVswYIFOohGfdnZ2SwiIkKpbPny5WzOnDlNPvYqVfWUyWTszJkzzNXVleXk5BhM/IcOHWITJkxgoaGhbNOmTay8vJz16NGDZWdnc9ts2rSJTZkyhTGmSDArV67kXisqKmJOTk7s3r17je6zISgoiB0+fLhauT6d+yZ5O0vdGYMbk9DQUMTFxaFnz55c2ZMnT/DgwQM4ODhwZfb29rh58ybEYjFycnKUpovh8/no0qULbty4oXK2ZH37Gzk6OmLdunXcc7FYjFOnTsHZ2bnJx17FwsICZWVl6NWrF8LDwxEWFoa2bdsaRPwFBQVYtWoVvvrqKxgZKT6ubt68CblcjldeeYXb7tn6Px+7lZUVWrZsiRs3bjSqzwaJRIKsrCzExcVh0KBBGD58OKKiovD48WO9Ovc6m8VXl9SdMbgxsbGxqVYmFAoBAObm5lyZubk5GGMQiUQ1/h3Mzc0hFArRrFmzRvc3EovFmDNnDkxNTTFy5EisXr3aYGI3NzfH33//jYyMDLz//vtc/Zty/DKZDHPnzsXHH3+Mzp07c+Xl5eUwNTUFn8/nyp6tv6rYKyoqGtVnQ3FxMfr06YOgoCD88MMPuHHjBiIiIiCRSADoz7lvkknEUGYCrnoTPRtr1ZtBIBDA3NwclZWVSvtUVFSgefPmje5v9ODBA8yYMQMAsGXLFu6q1BBiBxRXk3w+H3369MHYsWORnp4OoGnHv379etjY2FSbN08gEEAikUAulyu9D6rqryp2gUDQaGIHFBeOMTEx3PMePXogLCwM8fHxAPTn3DfJ21nqzhjc2LVs2RLW1tbIzc3lynJzc2FnZwdjY2M4OjoqvSaTyXDr1i04ODg0qtmSs7KyMGbMGNjZ2WH79u1o3bq1wcR+9uxZjB8/XqlMIpHA0tKyycd/4MABJCcnw93dHe7u7rhw4QK+/fZb7NixAzweT6nH0bP1fz724uJiPHr0iIu9sXw2ZGZm4ocfflAqq6yshLW1tX6de61bU/RYZWUle/3119mGDRtYZWUlO3HiBOvbty+7detWQ1ftpZ07d06pYX3lypUsODiYPXz4kN27d48FBgayH3/8kTHGWEZGBnN1dWV//PEH10vD19eXSaVSxhhjY8aMYYsWLWIikYhdvHiRubu7s9TU1AaJS5Xi4mLm7e2t1FBapanHzhhjJSUlzNPTk23bto1JpVKWmprK3N3dWUpKikHE/6xJkyZxvbNmzpzJIiIiWGlpKdc7Kz4+njHGWHJyMvPy8mJXrlzhemeFh4czxhrXZ8Pt27dZ79692d69e5lMJmOXLl1iAwcOZIcOHdKrc98kkwhjjGVmZrIJEyYwNzc35ufnx5KSkhq6SnXi+SQiEonY4sWL2auvvso8PDzYsmXLuDcLY4wdPXqU+fn5MVdXVzZx4kSWm5vLvVZQUMDef/991r9/f+bj48NiY2N1Gos6Nm/ezLp378769u3LXF1duZ8ZM2Y0+dirpKens+DgYNavXz82atQodvToUcZY0z/3z3s2iTx69IjNmTOHeXh4MG9vb/bTTz8pbRsTE8OGDh3K3Nzc2LRp09iDBw+41xrTZ8Pp06fZ6NGjmaurK/Px8WE7d+5kjOnXuadZfAkhhGitSbaJEEII0Q1KIoQQQrRGSYQQQojWKIkQQgjRGiURQgghWqMkQgghRGuURAghhGiNkgghhBCtURIhhBCiNUoihOiAk5MT4uLiEBQUhN69e+Ott95CWlpaQ1eLkJdG054QogNOTk7o0KEDFi9eDFtbWyxYsABlZWVISEho6KoR8lLomwghOhIaGorBgwfD3t4eU6dORWZmJsRicUNXi5CXQkmEEB2xs7PjHlctTyqVShuoNoTUDUoihOiIiYlJtTK6m0waO0oihBBCtEZJhBBCiNYoiRBCCNEadfElhBCiNfomQgghRGuURAghhGiNkgghhBCtURIhhBCiNUoihBBCtEZJhBBCiNYoiRBCCNEaJRFCCCFaoyRCCCFEa/8PzcAe7Zb2GV4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm import tqdm_notebook\n",
    "def get_cond_entropy_vs_n(mean):\n",
    "    conditional_entropy = []\n",
    "    plugin_estimate = []\n",
    "    sample_sizes = [i*200 for i in range(1, 26)]\n",
    "    for elem in tqdm_notebook(sample_sizes):\n",
    "        X, y = get_sample(elem, mean, 1)\n",
    "        conditional_entropy.append(rf_entropy_estimate_bagging(np.array(X).reshape(-1, 1), y, 300, .32, depth = 30))\n",
    "        plugin_estimate.append(plugin_estimate_cat_1D(np.array(X).reshape(-1, 1), y).item())\n",
    "    return sample_sizes, conditional_entropy, plugin_estimate\n",
    "\n",
    "#sample_sizes, cond_entropy, plugin_estimate = get_cond_entropy_vs_n(1)\n",
    "#plt.scatter(sample_sizes, cond_entropy, marker = \"None\",)\n",
    "#plt.scatter(sample_sizes, plugin_estimate, marker = \"None\")\n",
    "sns.set(font_scale = 1.2)\n",
    "line_plugin, = plt.plot(sample_sizes, plugin_estimate, label = \"plugin estimate\", linewidth = 3, color = \"green\")\n",
    "line_est, = plt.plot(sample_sizes, cond_entropy, label = \"rf estimate\", linewidth = 3)\n",
    "line_true = plt.axhline(y=.356, color='r', linestyle='-', label = \"true\")\n",
    "plt.xlabel(\"n\")\n",
    "plt.ylabel(\"H(Y|X) estimate\")\n",
    "plt.ylim(ymin = 0)\n",
    "plt.ylim(ymax = .8)\n",
    "plt.legend(handles=[line_est, line_true, line_plugin])\n",
    "plt.title(\"H(Y|X) Estimate vs Sample Size with Mean = 1\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}