{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Generative Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Preliminaries\n",
    "\n",
    "- Goal \n",
    "  - Introduction to linear generative classification with a Gaussian-categorical generative model\n",
    "  \n",
    "- Materials        \n",
    "  - Mandatory\n",
    "    - These lecture notes\n",
    "  - Optional\n",
    "    - Bishop pp. 196-202 (section 4.2 focusses on binary classification, whereas in these lecture notes we describe generative classification for multiple classes).    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Challenge: an apple or a peach?\n",
    "\n",
    "- **Problem**: You're given numerical values for the skin features roughness and color for 200 pieces of fruit, where for each piece of fruit you also know if it is an apple or a peach. Now you receive the roughness and color values for a new piece of fruit but you don't get its class label (apple or peach). What is the probability that the new piece is an apple?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- **Solution**: To be solved later in this lesson.\n",
    "\n",
    "-  Let's first generate a data set (see next slide)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAG2CAYAAAB4e1KRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABafklEQVR4nO3deXgUVdo28LvBJGxJWMISSFgcZRGNmOASZVcRUBSREREZ3CAuoMjo6yiEhMXBZXCZV8UQEXVwHSCOjoowQhBEFBJ44ZPFDQWlHYJCgogJJOf7o6zel+ru2uv+XVeuTldXV53qajhPznnOOS4hhAARERGRwzQyugBERERERmAQRERERI7EIIiIiIgciUEQERERORKDICIiInIkBkFERETkSAyCiIiIyJEYBBEREZEjMQgiIiIiR2IQRERERI5k2SBo/vz5cLlcmDZtmtFFISIiIguyZBC0efNmLFq0CDk5OUYXhYiIiCzKckHQL7/8gvHjx6O0tBStWrUyujhERERkUacYXYBY3Xnnnbj88stxySWXYN68eRH3ra2tRW1tred5Q0MDfv75Z7Rp0wYul0vrohIREZEKhBA4evQoOnbsiEaN1Gu/sVQQ9Prrr6OyshKbN29WtP/8+fMxe/ZsjUtFREREeti/fz+ysrJUO55lgqD9+/fj7rvvxqpVq9CkSRNF73nggQcwffp0z/Pq6mp07twZ+/fvR1pamlZFJSIiIhXV1NQgOzsbqampqh7XJYQQqh5RI2+99RauvvpqNG7c2LOtvr4eLpcLjRo1Qm1trd9rodTU1CA9PR3V1dUMgoiIiCxCq/rbMi1BF198MXbs2OG37aabbkLPnj1x//33Rw2AiIiIiHxZJghKTU3FmWee6betefPmaNOmTdB2IiIiomgsN0SeiIiISA2WaQkKpby8XPVjCiFw8uRJ1NfXq35somiSkpLYtUtEpBNLB0Fqq6urg9vtxq+//mp0UcihXC4XsrKy0KJFC6OLQkRkewyCftfQ0IC9e/eicePG6NixI5KTkzmhIulKCIGqqip8//33OP3009kiRESkMQZBv6urq0NDQwOys7PRrFkzo4tDDtW2bVt8++23OHHiBIMgIiKNMTE6gJrTcRPFiq2PRET6YY1PREREjsQgiIiIiByJQRApVlxcjD59+hhdDCIiIlUwCNKK2w0UF0uPOtm4cSMaN26MYcOG6XZOIiIiq2IQpBW3G5g9W9cg6IUXXsDUqVOxYcMG7Nu3T7fzEhERWRGDIJs4duwY3nzzTdx+++244oor8OKLL3peKy8vh8vlwrvvvouzzz4bTZo0wfnnn++3IO2LL76Ili1b4q233kL37t3RpEkTXHrppdi/f3/E8y5ZsgS9evVCkyZN0LNnTzz77LOe1+rq6jBlyhRkZmaiSZMm6Nq1K+bPn6/6tRMREcWDQZCa3G6gstL7A/g/17BV6I033kCPHj3Qo0cP3HDDDViyZAmEEH773Hffffjb3/6GzZs3o127drjyyitx4sQJz+u//vorHnroIbz00kv4+OOPUVNTg+uuuy7sOUtLSzFjxgw89NBD2LVrF/7617+isLAQL730EgDg73//O95++228+eab2LNnD5YuXYquXbtqcv1ERESx4mSJaiopkbrAfE2a5P29qEjKE9LA4sWLccMNNwAAhg0bhl9++QUffvghLrnkEp/TF+HSSy8FALz00kvIyspCWVkZrr32WgDAiRMn8PTTT+P888/37NOrVy989tlnOO+884LOOXfuXCxYsACjR48GAHTr1g07d+5ESUkJJk6ciH379uH0009Hv3794HK50KVLF02unYiIKB5sCVJTQQFQUSH9lJZK20pLvdsKCjQ57Z49e/DZZ595Wm1OOeUUjB07Fi+88ILffvn5+Z7fW7dujR49emDXrl2ebaeccgr69u3red6zZ0+0bNnSbx9ZVVUV9u/fj1tuuQUtWrTw/MybNw9ff/01AODGG2/Etm3b0KNHD9x1111YtWqVqtdNRESUCLYEqSkzU/rxlZsr/Who8eLFOHnyJDp16uTZJoRAUlISDh8+HPG9gTMUh5qxONS2hoYGAFKXmNxyJJOXe8jNzcXevXvx/vvv4z//+Q+uvfZaXHLJJVi2bJmyCyMiItIQgyCLO3nyJF5++WUsWLAAQ4cO9XvtmmuuwSuvvIIzzzwTALBp0yZ07twZAHD48GF88cUX6Nmzp9+xtmzZ4un62rNnD44cOeK3j6x9+/bo1KkTvvnmG4wfPz5s+dLS0jB27FiMHTsWY8aMwbBhw/Dzzz+jdevWCV87ERFRIhgEaSUzU8oBCmwZUtm///1vHD58GLfccgvS09P9XhszZgwWL16MJ554AgAwZ84ctGnTBu3bt8eMGTOQkZGBUaNGefZPSkrC1KlT8fe//x1JSUmYMmUKLrjggpD5QIA0eeJdd92FtLQ0DB8+HLW1tdiyZQsOHz6M6dOn44knnkBmZib69OmDRo0a4Z///Cc6dOiAli1bavVxEBERKcacIK1kZkpJ0BoHQYsXL8Yll1wSFAABUkvQtm3bUPn7SLWHH34Yd999N/Ly8uB2u/H2228jOTnZs3+zZs1w//334/rrr0d+fj6aNm2K119/Pey5b731Vjz//PN48cUXcdZZZ2HgwIF48cUX0a1bNwBAixYt8Mgjj6Bv374499xz8e233+K9997jIrVERGQKLhE4jtrGampqkJ6ejurqaqSlpfm99ttvv2Hv3r3o1q0bmjRpYlAJtVFeXo7Bgwfj8OHDYVthXnzxRUybNg1HjhzRtWzkz87fQyKieEWqvxPBP8mJiIjIkRgEERFRkOJiYO7c0K/NnavZlGdEumIQ5ACDBg2CECJiQvKNN97IrjAi8mjcGJg1KzgQmjtX2v77TBhElsbRYUREFKSwUHqcNcv7XA6A5szxvk5kZQyCiIgoJN9AaN48oK6OARDZC7vDiIgorMJCIDlZCoCSkxkAkb0wCCIiorDmzvUGQHV14ZOliayIQRAREYXkmwNUWys9hkqWJrIq5gQREVGQUEnQoZKliayMQRAlzOVyoayszG8dMiKytvr60EnQ8vP6ev3LRKQ2doepyKjJxW688Ua4XC64XC4kJSXh1FNPxb333otjx45pc0Iisr3i4vAtPYWFnCyR7IFBkIqMnFxs2LBhcLvd+OabbzBv3jw8++yzuPfee7U7IRERkcUxCFJRYWFw4qBek4ulpKSgQ4cOyM7OxvXXX4/x48fjrbfeghACjz76KE499VQ0bdoUZ599NpYtW+Z5X319PW655RZ069YNTZs2RY8ePfDUU08FHf+FF15A7969kZKSgszMTEyZMsXv9UOHDuHqq69Gs2bNcPrpp+Ptt9/2e33nzp0YMWIEWrRogfbt22PChAk4dOiQ5/Vly5bhrLPOQtOmTdGmTRtccsklbMkiIiJNMQhSmW8glJJi3OyqTZs2xYkTJzBz5kwsWbIECxcuxOeff4577rkHN9xwA9atWwcAaGhoQFZWFt58803s3LkTs2bNwoMPPog333zTc6yFCxfizjvvxOTJk7Fjxw68/fbbOO200/zON3v2bFx77bXYvn07RowYgfHjx+Pnn38GALjdbgwcOBB9+vTBli1bsHLlSvz3v//Ftdde63l93LhxuPnmm7Fr1y6Ul5dj9OjREELo9GkREZEjCQeprq4WAER1dXXQa8ePHxc7d+4Ux48fV+VcyclCANKj1iZOnCiuuuoqz/NPP/1UtGnTRowZM0Y0adJEbNy40W//W265RYwbNy7s8e644w5xzTXXeJ537NhRzJgxI+z+AMTMmTM9z3/55RfhcrnE+++/L4QQorCwUAwdOtTvPfv37xcAxJ49e0RFRYUAIL799ltF12tnan8PiYjsIFL9nQiODtNAqMnFtG4J+ve//40WLVrg5MmTOHHiBK666irce++9WLZsGS699FK/fevq6nDOOed4nj/33HN4/vnn8d133+H48eOoq6tDnz59AAAHDx7EgQMHcPHFF0c8f05Ojuf35s2bIzU1FQcPHgQAVFRUYO3atWjRokXQ+77++msMHToUF198Mc466yxcdtllGDp0KMaMGYNWrVrF+3EQERFFxSBIZYE5QPJzQNtAaPDgwVi4cCGSkpLQsWNHJCUl4dNPPwUAvPvuu+jUqZPf/ikpKQCAN998E/fccw8WLFiA/Px8pKam4rHHHvO8t2nTporOn5SU5Pfc5XKhoaEBgNTlNnLkSDzyyCNB78vMzETjxo2xevVqbNy4EatWrcL//u//YsaMGfj000/RrVu32D4IIiIihRgEqcjIycWaN28elKdzxhlnICUlBfv27cPAgQNDvm/9+vW48MILcccdd3i2ff31157fU1NT0bVrV3z44YcYPHhwXGXLzc3F8uXL0bVrV5xySuivnMvlwkUXXYSLLroIs2bNQpcuXVBWVobp06fHdU4iIqJoGASpyGyTi6WmpuLee+/FPffcg4aGBvTr1w81NTXYuHEjWrRogYkTJ+K0007Dyy+/jA8++ADdunXDP/7xD2zevNmvBaa4uBi33XYb2rVrh+HDh+Po0aP4+OOPMXXqVEXluPPOO1FaWopx48bhvvvuQ0ZGBr766iu8/vrrKC0txZYtW/Dhhx9i6NChaNeuHT799FNUVVWhV69eWn00REREDILUFGnyMKOml587dy7atWuH+fPn45tvvkHLli2Rm5uLBx98EABw2223Ydu2bRg7dixcLhfGjRuHO+64A++//77nGBMnTsRvv/2GJ554Avfeey8yMjIwZswYxWXo2LEjPv74Y9x///247LLLUFtbiy5dumDYsGFo1KgR0tLS8NFHH+HJJ59ETU0NunTpggULFmD48OGqfx5EREQylxDOGYdcU1OD9PR0VFdXIy0tze+13377DXv37kW3bt3QpEkTg0pITsfvIRFRsEj1dyI4TxARkQMYtawPkZkxCCIicgAjl/UhMivmBBEROUCokap6LetDZFYMgoiIHMI3EJo3T5rMlQEQORm7w4iIHKSw0DubfXIyAyByNgZBREQOEmpZHyKnYhBEROQQvjlAtbXSY6hkaSKnYE4QEZEDGLmsD5FZMQgiIjKR4mJpuHqogGTuXGn5nXjm9DHbsj5mptU9IPNhd5gG6uvrUV5ejtdeew3l5eWot9H/LsXFxejTp4/RxTAVfiakpnjn84k2GSIQvqWnsDBype60iRY5p5JzMAhS2YoVK9C1a1cMHjwY119/PQYPHoyuXbtixYoVmp9748aNaNy4MYYNG6b5uSg2Dz/8MHr37o1mzZqhe/fuePXVV40uEplUYWFwro6S+Xy0rLidFhTEew/IgoSDVFdXCwCiuro66LXjx4+LnTt3iuPHj8d9/OXLlwuXyyUA+P24XC7hcrnE8uXLEyl+VLfccou4++67RfPmzcV3332nyTmKiorE2Wefrcmxw2loaBAnTpzQ9ZyxUPKZjBgxQnzwwQfi66+/FrNnzxaNGzcWX3/9ddB+anwPyR7mzBECECI5WXqcM0f5e+R9A5+rUR4tjm1W8dwD0kak+jsRDIJ+l2jlc/LkSZGVlRUUAPkGQtnZ2eLkyZOJXkZIv/zyi0hNTRW7d+8WY8eOFbNnz/Z7fe3atQKA+Pe//y1ycnJESkqKOO+888T27ds9+yxZskSkp6eLsrIycfrpp4uUlBRxySWXiH379nn2CVXhv/DCC6Jnz54iJSVF9OjRQzzzzDMRy/rbb7+JqVOnirZt24qUlBRx0UUXic8++yyorCtXrhR5eXkiKSlJrFmzJug4e/fuFQDEa6+9JvLz80VKSoo444wzxNq1a/32+/zzz8Xw4cNF8+bNRbt27cQNN9wgqqqqPK+///774qKLLhLp6emidevW4vLLLxdfffWV3zH2798vxo4dK1q1aiWaNWsm8vLyxKZNm/w+k5dffll06dJFpKWlibFjx4qampqQ1//TTz8JAGL9+vVBrzEIIl9y5ZucrPw9WlbcRgYFRUXhzzdnjvS6FuK5B6Q+BkEq0DIIkivuaD+BFbRaFi9eLPr27SuEEOKdd94RXbt2FQ0NDUHl69Wrl1i1apXYvn27uOKKK0TXrl1FXV2dEEIKgpKSkkTfvn3Fxo0bxZYtW8R5550nLrzwQs9xAoOgRYsWiczMTLF8+XLxzTffiOXLl4vWrVuLF198MWxZ77rrLtGxY0fx3nvvic8//1xMnDhRtGrVSvz0009+Zc3JyRGrVq0SX331lTh06FDQceQgKCsrSyxbtkzs3LlT3HrrrSI1NdWz/4EDB0RGRoZ44IEHxK5du0RlZaW49NJLxeDBgz3HWbZsmVi+fLn44osvxNatW8XIkSPFWWedJerr64UQQhw9elSceuqpon///mL9+vXiyy+/FG+88YbYuHGj5zNp0aKFGD16tNixY4f46KOPRIcOHcSDDz4YVOaGhgZx4403ijPPPFPU1tYGvc4giGSJBBxaVtxGBQXhWp60bJFiS5B5MAhSgZZB0KuvvqooCHr11VcTvYyQLrzwQvHkk08KIYQ4ceKEyMjIEKtXr/a8LgcWr7/+umfbTz/9JJo2bSreeOMNIYQUBAHwtHAIIcSuXbsEAPHpp58KIYKDoOzs7KBrmjt3rsjPzw9Zzl9++UUkJSWJV155xbOtrq5OdOzYUTz66KN+ZX3rrbciXrMcBD388MOebSdOnBBZWVnikUceEUIIUVhYKIYOHer3vv379wsAYs+ePSGPe/DgQQFA7NixQwghRElJiUhNTfUEaYGKiopEs2bN/Fp+7rvvPnH++ecH7XvzzTeL7t27i++//z7ksRgEkRCJdT3ZtSXI9/yBn8uQIeH3j7eFyIndf3qKtWWPQZAK7NoStHv3bnHKKaeIH3/80bPtzjvvFOPGjQsqX2CuUJ8+fURxcbEQQgqCTjnllKAuu5YtW3padnyDIDlYaNq0qWjevLnnJyUlRbRr1y5kWf/v//5PABDffvut3/ZRo0aJm266ya+s4QIFmRwErVu3LuhYN954oxBCysVJSkryK1/z5s0FAPHee+8JIYT46quvxLhx40S3bt1Eamqq5/V3331XCCHE7bffLgYMGBC2HEVFReKMM87w2/b444+Lbt26hbz23bt3hz0WgyBKpMXDCTlBgYHYkCHqtxAZ0erkNLF+xloFQZwnSCX9+/dHVlYWfvjhBwghgl53uVzIyspC//79VT/34sWLcfLkSXTq1MmzTQiBpKQkHD58GK1atYr4fpfLFfF5uG0NDQ0AgNLSUpx//vl+rzUOM1xE/mwCjyeECNrWvHnziOWORD5WQ0MDRo4ciUceeSRon8zMTADAyJEjkZ2djdLSUnTs2BENDQ0488wzUVdXBwBo2rRp1PMlJSUFnV/+fGR79+4FAPTo0SP2CyLHiHc+n3gmQ1Q6H46ZJlosLPQu/pqcDHz4obd88uuJjuTinEraC/X9MWIEHoMglTRu3BhPPfUUxowZA5fL5RcIyRXyk08+GTY4iNfJkyfx8ssvY8GCBRg6dKjfa9dccw1eeeUVTJkyxbNt06ZN6Ny5MwDg8OHD+OKLL9CzZ0+/423ZsgXnnXceAGDPnj04cuSI3z6y9u3bo1OnTvjmm28wfvx4ReU97bTTkJycjA0bNuD6668HAJw4cQJbtmzBtGnTYrp232saMGCAp/wVFRWea87NzcXy5cvRtWtXnHJK8Nf9p59+wq5du1BSUuIJUDds2OC3T05ODp5//nn8/PPPaN26dVxlBICBAwdi8+bNcb+fnCHSnDuRKoZ4Km556HvgsX0ro3iPrZVQa5/5VqhygJRIRRrvPaDYqH3f4qJqu5LJaT1EXghpmHzgKLHs7GzNhseXlZWJ5ORkceTIkaDXHnzwQdGnTx8hhLeLqXfv3uI///mP2LFjh7jyyitF586dPQm6cmL0eeedJzZt2iQqKipEfn6+uOCCCzzHDMwJKi0tFU2bNhVPPvmk2LNnj9i+fbt44YUXxIIFC8KW+e677xYdO3YU77//vl9i9M8//+xX1sOHD0e8drk7rHPnzmLFihVi165dYvLkyaJFixae0V8//PCDaNu2rRgzZoz49NNPxddffy0++OADcdNNN4mTJ0+K+vp60aZNG3HDDTeIL7/8Unz44Yfi3HPPFQBEWVmZEEKI2tpa0b17d9G/f3+xYcMG8fXXX4tly5b5JUYHjph74oknRJcuXfy2rVixQvTo0SPiNbE7jPRmlm4uJaKVlSO5rEnJfWNOkAr0CIKEkIbLr127Vrz66qti7dq1mg2LF0KIK664QowYMSLkaxUVFQKAqKio8AQW77zzjujdu7dITk4W5557rti2bZtnf3mI/PLly8Wpp54qkpOTxZAhQ/zyd0JV+K+88oro06ePSE5OFq1atRIDBgwQK1asCFvm48ePi6lTp4qMjIyIQ+SVBkGvvvqqOP/880VycrLo1auX+PDDD/32++KLL8TVV18tWrZsKZo2bSp69uwppk2b5hk9t3r1atGrVy+RkpIicnJyRHl5uV8QJIQQ3377rbjmmmtEWlqaaNasmejbt2/YZHEhQgdBcuJ5JAyCyAhGJzwrES2HRM4NMvM1UDCl3z0GQSrQKwgyIyWBhRwEWYUcBG3dutXooqjG7t9DMi8tWlHUnNsn0rECk6PN3JpFXrG0QjIxmoiIVBGYEB2YZzN4MLB2beLnUZpzpLTMocydC6xZY46kbVLOLMn2DIKIiBzGNzgB/AOSWbOA8nL/hON46TECyExJ26ScWe6bS4gQ47ltqqamBunp6aiurkZaWprfa7/99hv27t2Lbt26oUmTJgaVkJzOyd9DpcO1SR2+w8p9A6DA39UIVORzyS1NXISUYhWp/k6EpVaRX7hwIXJycpCWloa0tDTk5+fj/fffN7pYRKQCp61UbrTCQmDQIOn3efP8gx55FXW1/hovLPQGQMnJDIDIPCzVHZaVlYWHH34Yp512GgDgpZdewlVXXYWtW7eid+/eqpzDQQ1jZEJO/v6ZZfI0J1m7FkhJCR2cqPl5R5rbh8hIlmoJGjlyJEaMGIHu3buje/fueOihh9CiRQts2rQp4WPLM/7++uuvCR+LKF7yLNVqT6ppFXILxKxZUuXMAEhboYITLc4h38faWu/91eJcRLGyVEuQr/r6evzzn//EsWPHkJ+fH3Kf2tpa1NbWep7X1NSEPV7jxo3RsmVLHDx4EADQrFmzkEtFEGmloaEBVVVVaNasWcjZrZ0icFkEBkDaCGxlC1x6Qotz+B6bI7fIDCz3P+2OHTuQn5+P3377DS1atEBZWRnOOOOMkPvOnz8fs2fPVnzsDh06AIAnECLSW6NGjdC5c2dHB+DsOtGeXsGJWUYAEYVjudFhdXV12LdvH44cOYLly5fj+eefx7p160IGQqFagrKzs6Nml9fX1+PEiROalJ8okuTkZDRqZKlealWFa51gl5i6OBKPrEar0WGWC4ICXXLJJfjDH/6AkpKSqPtq9SESUeLCBTwMhIhIq/rbct1hgYQQfq09RGRN7DohIr1ZKgh68MEHMXz4cGRnZ+Po0aN4/fXXUV5ejpUrVxpdNCJKUKTuF7YAEZEWLBUE/fe//8WECRPgdruRnp6OnJwcrFy5EpdeeqnRRSMiIiKLsVQQtHjxYqOLQERkW0yYJqdx7jAUIiLyw6VLvIqLw0/oOHcug0G7YBBEREQA/GfslgMAp47Oc3pA6JQg0FLdYUREpC3fSRPlmbudFgABXMtODgKB8FNW2IHl5wmKBecJIiJSxndhVSfPQiJX+vIM5k4IgGRmmryUkyWqgEEQEVF0Tq74Q3FyQGiW74JW9TdzgoiIyIOrvvsLtZadkxQWeq/djgsaMwgiIiIA4RdWdWogxIDQ/kEgE6OJiAgAly7xFS4gBEInDNtRuJwgwD7XziCIiIh0Y5UJGZ0eEDolCGQQREREAPQZFm2VoddOX8vOKUEggyAiIp2ZtTVEj7lxnD7/jlU4JQhkEEREpDMzt4boMVkiJ2Qks+A8QUREBjDTRHSh6DE3jpPn36HYcJ4gIiIb8R16npJirgBIj2HRdh96TdbAIIiIyCBmnIhOj7lxAs8xaFD4c9hpsU4yHwZBREQGMao1JNwK4XJwMmiQdpMlhur2GzJEegw8h1NWbDcrJ6wkzyCIiMgARs5GLCdmB55rzRrpUQ5KZHIgpMaw6FBDr+Xj+5YhUo6UFSpnK5QxmnDfE1sFp8JBqqurBQBRXV1tdFGIyMHmzBECkB6VbNejDHqeO1qZkpMjl8UMn180ViijEmb5nmhVfzMIIiLSWVFR5Aq+qCix/ZVSGnToSS5LcnLk/cxSOUcSTxm1uteJMMP3hEGQChgEEZEVadmqoDTo8GWWoMwMlXM08V6T2VqQ4vmeqIlBkAoYBBGZhxn/4jUzLVo+4g0itKio470+oytnJWIto9laucwQbDIIUgGDICLzMOtfvGamZmWUaEWrZkUd73fBDJVzNIkGmkZfm1kCMgZBKmAQRGQuZvkP1krUaPlQKwBVq6KOp1XQCt+dRMtodCuXmf5QYRCkAgZBROZjlr94rcDIoCMcIypqM1XO4SRaRjP8uzBTlzWDIBUwCCIyJ6P/4rUCM7Z8GFVRm6lyDieRMprxXhuNQZAKGAQRmY/Rf/FaoUI1Y8sHK2ptmPFem4FW9fcpek/OSEQkC7eSOqDfOlryrLiB5/Qtm9FCzbIMeJ+rMZNzLELN5iw/6n3/7MZs99ruGAQRkSHMUpGGOmekJRuMEGmJBSPKx4paO2a713bnEkIIowuhl5qaGqSnp6O6uhppaWlGF4fI0YqLpVaYUP+xz50rVaSJrK8U6/HlwEdezNQsARARaVd/cwFVIjJEcXH4IKOwMPEFJmNd/LGw0BsAJSfbOAByu6UP1+02uiQeVlps1EplpegYBBGRLckrk/sGQpG6uebO9QZAdXX6rOZuCLcbmD3bVEGQlVYrt1JZSQFV06xNjqPDiJxHyegzR410qqiQLq6iwuiS+LHSPbBSWe1Cq/qbOUFEZHspKd5Wntpa/9fCtQ6ZLTkaSCCPyu32tvxUVgKTJgGlpUBurrQtM1P6MZiV8rKsVFY7YE4QEVEI0XI0Bg+O3M0VaaTTnDnmGukUd1dMSQmQlyf9TJokbZs0ybutpETTcitlpbwsK5WVIlC1Xcnk2B1GZD/RJpezW7dFXF0xBw5I3V8VFUKUlkpvKC31bjtwQJeyR2P0xJmxsFJZ7YAzRquAQRCRPYULDOw68268FXBRkRBzbvshZE6Q0bNjWynPxkpltQsGQSpgEERkX4GBwaBB5l8OIxHxrLfmqawx0y8IMroSt9JSEVYqq51w2QwioggKC4F587w5GmvXRt7XykIN51dyTYWFAI4exazH5gKvH0VhrjkSwK00A7WVykoKqBpSmRxbgoiMp9WCpVrnaJhloVU1umKsks9ils88VlYtt5lpVX9zdBgR6UqLyeZ8WzNqa4MnSVSDGSbJC7feWqzXa5WRTWb4zONh1XI7kqohlcmxJYjIHNRMLNUzR0PXhNgDB6QmA5+RW2q1MFilJUgI6yYhW7XcZsXEaBUwCCIyD7UqYr27HnQLIDSa2dmKlbOVgjZfVi23GXHGaBVwxmgic4k0k7OZ6VLuykppIsOKCu/Mzgmy0uzYgfhdcTbOGE1EtmLVBUs1LbfbLQU/8g/g/zzBRU+Vzo5ttpXS+V0hzajarmRy7A4jMgcju2QS6T7TvNxFRd5ZHkP96DSsSO+5cCLdkyFDrNd9J4Q1ux3NjDlBKmAQRGQ8oyebi/f8upTbRMtb6FmJhzu2HAANGaJsf7Mw+jtuR5wskYhswejJ5uTzzJrlfa4kJ0aXcodazT03V7WcoFj4fk7yJJRa5QyFuydr1gBDhgAffhh6f7NOTGj0d5yUY2I0ETmSHPjI+RqmSwrWIDE6Hnom9pr+npBhmBhNRJQot1vK6nW7zT9hYGYmUFQU3DKkI70Te01/TygksyXSx4JBEBE5h9sNzJ4NuN2qVvCaVAKZmdIbDQqC9JiFO9Q5OZrKeiw9Q7aqGUYmx8RoIof7fQLCObf9oGrSr90SYY24Ho6msjat7x8To4mI4uF2e+fXqazEXMzErOc6Ys5tB1B4+Y+AOxOFhVJri29ibiziTbY2K70Te8OtiQbEf09IX3om0quJidFEZG/FxVIXmPwURWiMehRinrShqMjTXzV3rlTBx5vDIFfmjRoBDQ2hK4FEz2FHxcVSl0moCpOfl7VolUivVf3NIIiI7C2gJQiTJgGlpd4RV6GGpSdArgQAay5PQRQvLUf3aVV/szuMiOxNx7l3AhN77dI9RhRN4Pdbfg6Y+/vO0WFERCoINZoKkLalpDAAIvsKl9Olx4jCRLEliIicQ6O5d6Il9nLeG7IzK8+QzSCIiJxDnntHZeEqAVnjxt55bxgIkd1E+idl9u87gyAiogSFqgSsmiNB5CQMgoiIVMZ5b2zI7QZKSoCCAkOXMiF1MQgiIk04ee4XK+VIxHKfnHxPPUuuXHklgyAb4egwItKEpdcTSlBxcfiWnsJCcwUKsdwnJ99Tsie2BBGRJuy2lISlxNB1E8t9ctw9DZxo0/cRUH2iTTKAqiuRaeyvf/2r6Nu3r2jRooVo27atuOqqq8Tu3bsVv58LqBLpT15IMTmZC2Lq5veFYkVFheK3xHKfHHNPi4qkCwz3U1RkdAkdQ6v621LLZgwbNgzXXXcdzj33XJw8eRIzZszAjh07sHPnTjRv3jzq+7lsBpExtFpPiMKorATy8oCKiphmxo7lPjninuq85AqFx2UzAKxcudLv+ZIlS9CuXTtUVFRgwIABBpWKiCIJXEqCc+VoJMGum1juk2PuqY5LrpAxLJ0YXV1dDQBo3bp1yNdra2tRU1Pj90NE+gm1lITZp9G3rJISqfUnL09qsQCkR3lbSUnYt8Zyn3hPyVZU7VzTUUNDgxg5cqTo169f2H2KiooEgKAf5gQRaU/OGwnMFwm3nRJ04ICUA1RRIURpqfQhl5Z6tx04EPJtsdynkPseOCDmDPow4j0tKgr/2pw5FkmtOXBAKmiYz5G0pVVOkKW6w3xNmTIF27dvx4YNG8Lu88ADD2D69Ome5zU1NcjOztajeETO9fvIpPqjf8acOamWmCvHFuLsuollTqOQ+7rdKCy/GLjtB9TXdwx5Dnlove9xAf9WJdPTaMkVMpiqIZVOpkyZIrKyssQ333wT0/s4OoxIB3GMTDKMXf+61/MeKDxXYCsSWwQpFmwJAiCEwNSpU1FWVoby8nJ069bN6CIRkZXZdRbgzEygqEi7a4ojCdt3jqF586SEalvOLUSWYqkh8nfccQdeffVV/Otf/0KPHj0829PT09G0adOo7+cQeSKNWHUocZxDyR2vuFgKHsMpKgrbdeSIofWkOg6RB7Bw4UIAwKBBg/y2L1myBDfeeKP+BSIiSUlJcKUoj1ACIlaKqohlcUvOApy4ggKp9QwIH/SG4Jih9WQZlgqCLNRoReQscVaKqomlW8vogM1s4lkdPY4k7MDlNeTngHaBkKMXfCVFLBUEEZFBolWUVppUzuiAzWx0yIsKtb5YqHXI1GaLUWmkKQZBRBSdGROIw3RrFZdkonEjgcJZrqCySn/9Z6K42OCALVpQGU/rjEaitqYcOR3FUZKwYxmGrybHLfhKMWMQRETq0npkkixMt1ZjzMQszAV2r0HhWm8ZTPXXf7SgUuugM4a8qOitKalAYXHE00XqctI6EOGoNIpI1QH3Jsd5gohiEOcMxGYo35zbflA2J41R8wRFm1tH63l+Ylwd3Q5z/CQnNwhAeiTr4TxBRKQvsycQR8hDKlwIoKOCv/5jnAU4oUTbaK0vjRoBDQ3hXw83ai3WkXElJcCoUTHlRVm9NUUaleZCMmpRV5fCUWnkpWpIZXJsCSKKgdlbgnyFaTlJTha///WvzmkSWg8tWuvLwIExtc54xNJqFGrfGN6v9uepB8+9ue0HIYCgVkKyBrYEEZG+rDTiK0QekhZz0sSdaOt2A0ePAitXAm3bhm59CWwJMtmoNSvO8TP3f45i1mOpmHPbARTmvQcA0uNtIzBrVkfg6FEUPppqcCnJSAyCiMj6Arq1tJyTJq6uIbcbePxxYPx4/yAyWlAZ7vVYJnxU0g0XJZE90c/TqPl66j+rwBx8iMLn5nk3TpoEqRgzUf/ZxQAGqX9isg5V25VMjt1hRHHSOIG4qCh898ScOeF7gsLtH3eXVQxi6hoK7HJKNDE6lsTmGJOgA6nxeep1T4JYqUuXImJ3GBEZJ8YE4ljFNaldmIRgPeakUdQ1FKkFpqoKmD49fOtLtGkGYpnwMcHJIdX4PA2br8dKXbpkDFVDKpNjSxCRecU8DDvRYeRxtm4pLmeCLTCKJZoYrSP5s5Jb0XRNTjb42ikxbAkiIlvTfRh2HBMSxrT8QwItMHZd86qw0Htvk5N1TqzWaxJPspRGRheAiEhWWOjtYkpOBgpvdUu1vdyt5HZLAYX8A/g/l/fTSKSuoTlzArqGMjO9XS9y4OP7PEJlLHcPzp3rv10Owho3DjiP0srd4EAgVDei5txub8RYXMwgiPyp2q5kcuwOIzK3oO6S3+d28SSxxjuXjkxhoqyaidpCiLi6YuwwS7Mvw66H3WC2oFX9zSCIiEwhbCWJmd4gBRBi6dL4R/sozNNRfTRTgvlHhuTQqMiw0WFCMAiyCeYEEZFt+eXa3OoGKt0ovBxARQNm/WsuMGcLCv+4R9o5IyO2uXZ8KczTUX00U5yj6wzNoVGR2iP2oq9sfxTF47+UNsSyBAk5j6ohlcmxJYjInPy6nwJaa+ZgpiiCz7aRI6W/6pcuTewvfAUtBEa3xBh9frOK2rI06EN9RuaRbtgdpgIGQUQW4Ju3M3Jk9PW24p3wTmE3iVHrZdktJ0htET8fTpJoO+wOIyJtxbIauZZ8uyruvBN45x1g6VLg+PHQ3VfxllXBSCmj1suKaSi+Q0WeUoGTJJIyHCJPRBJ53px4h5m7A4azq6FtW+mxV6+Yh5lHJefphDmGbyBSWys9hhq2roWYhuI7WNCUCg4PDCl2DIKISB2JBlGhGDSvTbiWGL0CoeLi8BV64a1uFKNY8zmRrEDRvEOcJJEiYHcYkZPFshq5EQJHVelUmemx/ljc4pjp2o4Ur2yv8bp3ZG0MgoicrKREqlB9TZrk/b2oKHIFomcQpWNlFuk0pupyMUsel86YM0VqYRBE5GQJrjCecBDlBEoDlUj7hQs2mzaVPv/8fEcFQaZuqSNLcQkhRKxvOn78OH7++Wd06tTJb/vnn3+O3r17q1Y4tdXU1CA9PR3V1dVIS0szujhE5lJZCeTlARUVykfRBFbOao7eMorarSsffAAMGwasXAlcdln4/SJ9/sXFwcGmr8mTpTIT2ZRW9XfMidHLli1D9+7dMWLECOTk5ODTTz/1vDZhwgTVCkZEFpDAIqGmpXaC96FD/o/xKCiQgqOKCuCxx6RtM2dKPwBw+unKFpHVYgQfAQDq6+tRXl6O1157DeXl5ahnc5QlxNwdNm/ePFRWVqJt27bYsmULJk6ciBkzZuD6669HHI1KRGQWHEWjHt8Wst27vY9yV5bcQqY0p8q3RU1u8Zk3z7vfffd5f4/UBcmkak2sWLECd999N77//nvPtqysLDz11FMYPXq0gSWjqGKdXfGMM87we37o0CExYMAAMXv2bHHOOeeoM4WjRjhjNJGGYl0kNM5FRRMS7pxqzzA8fXrkma6nT5f2U7igq5+VK70LycZaVjMvJmrE90EFy5cvFy6XSwDw+3G5XMLlconly5cbXURb0Kr+jrk7rF27dti+fbvneZs2bbB69Wrs2rXLbzsR2ViobpUokw+GPIba8wrFe86SEikfJy/Pm9g9aZJ3m1b5Nr7dXKWl0rbSUu+2goLg9+TkSK09Q4Yo64J0u71dZb6tTUq6z/RkxPchQfX19bj77rtD9oLI26ZNm8auMRNTHAQdPXoUAPCPf/wD7dq183stOTkZr732GtatW6du6YhIf0ryRqxWYcnXVFUV+vV4gpFI7r3X+145b2fmTO+2e++VtsWTUxUu2Ax334wK8Bxg/fr1fl1ggYQQ2L9/P9avX69jqSgWinOC+vfvj5UrVyIrKyvsPhdddJEqhSIijUUbjq1V3ogRkzO63cCaNdI1yQFJuLwbX4msNeV7vF27pMeePdVfu8o3jyvcfYs0DUJVFbBqlfReI3KEzD5ZZxRuhX8EKN2P9Ke4Jahv3744//zzsVtO8vvd1q1bMWLECNULRkQaiqclR41uFSNaJUpKgBtukH6Xk4n1bAnJyPB/DCeexHQlXZCRWpvatgUef9y4Fj2Lt1JlKrxXSvcj/SkOgp5//nncfPPN6NevHzZs2IAvvvgC1157Lfr27YuUlBQty0hEWqmqkirRbduiBzhqVFihup1GjpTm0FHS7RTLEG85aMvP97YA/elP0mP//sArr4Q+p9qj5OQcnpycyPvFmlMF6J/vo/YQe7W7IXXWv39/ZGVlweVyhXzd5XIhOzsb/fv317lkpFismdR//etfRZMmTURSUpK44oorRIUZRxmEwdFh5GihRkDNnCk9jh4dfZSS2iOo5JFKsYxWimV0U7SRV7m56o5EMmJ0U6yjyw4ckEamrVwZ333UcnSZmUeuRSCPDgscIcbRYerSqv5WHAQdOHBATJ06VTRt2lTk5uaKZs2aiaVLl6paGK0xCCJHi1ZhTp6svGJUo8LSOgiKFPTFel61y6aWeALTeIblyxgEhbR8+XKRlZXlFwRlZ2czAFKRVvW34sToU089FT179sQ///lPXH755fjggw9w7bXX4vvvv8f999+vfhMVEalLTpCtqgI2bJDyY/70J+Dll6Xuon79pP2ys6XHRBKDI5GTYauqpK6wd96JnAwbb/Ks73Z5VFj79t5zAt6kZZMn4IYVT0J3rOvF6ZW8bOHJOkePHo2rrroK69evh9vtRmZmJvr374/GjRsbXTSKRmm09NprrwVtq6ioEJmZmeL2229XNTLTCluCiISyFqFof5En0vUTa0tEIi0XMvmaIl1zPNeidhdhIuJpSVHyHjU+f6IEGd4SdN111wVty83NxcaNGzk6jMhKCgqktaZuuEFqAZo3z78loFGj6H/dy0m8sXK7gaNHpUTotm2VtUQkutI9ANx+O7Bli38Lhq9Fi+K7ppKS4IVN5aRxIPISFiooLgYaN/599fSAlpS5c6XV1BM+vRqfP5FJxbx2WKCuXbvi448/VqMsRNYVbeVxtVcmT0RmpjTbcFGRNHIKCO5C6dNHm3O73dKQ7PHjvfPUAFIXXLguHDXm8OnTB/j3v6Xz79rlHTKfaGWuRoCwbRswbRrw5JMxf+6NGwOzZkm/FxZ6g7i5c6Xtc+ZEOYCSLii151AiMpGEgyAAaNWqlRqHIbKuaBMMmm3hSrnVI1zLiF7UWGFdKS0qczWO+fnnwLp10mOMQVBhofToDYT8AyD59bDibdELxUyBPpFCqgRBRGRRWiajypXiqFFAQ4O0LTCx9vBh6dF3IsFIlaka5c3MlFqCli4Nv4yGhfgGQvPmAXV1CgOgeET6/M0W6BMpwCCI7EmPv0qjjZpp1Ch85Q+YY0SSmi0BMvmzz8/3zkq9aJH/Pr55MwCwf7/3s6mqCl+ZqlFe3yBITbEEaNu2SS0/gLRshe8jAPTuHVOrUGGhNwBKTtYoAAK0+b4QGUnVNGuT4+gwB9FjzpFoo2YGDnTmqBr5s1+6VHqUJ+bzHUGlZHRavCO2Yimj/P3Qe6LDaN+NgQNjOtycOdLbkpOlxzlzNCl1MDONjiNbM3x0GBEFiJYUG9gS5IRRNXLiMQDI6wzu3y+NBAO8cxAtXQr06uX9XB57zPt5yK1rixYBAwZI+8mtZom08EVqudu1S2p9ys/XJ6/lySf9W4Jeflmas2noUGlb796KDyXnAA0aJOW7A/45QvI+qowUC2Tw6DiiRDEIIvvQe0XqeJJi7TqqRv7sS0q8XV++i5XKJk+WHnv18v8cvvwSuO++4OPKo7jkyjSRvJNoFTYgJWhrndfidgNvveUfaL38shQAjR8f06ECR4HJv8+Z4w2EfLerjsPnyeIYBJF98K9S44T67H1NnixVmOHmIBo9Wvo5dEhqQZIDqJkzgZ49pcTpRBftDFVhy8eXzym3Xu3aFX/QrGS6BDUCLbcb9R/uwZz78lBYmOrZ7BsIrVkDlJdrnCjN4fNkZap2rpkcc4Jszsj8hGg5JUYsrqmnSOt0LV2q7HNRMpO1WvfVd92yaGupyccPvIfh7mm0fLTA17dulXKAtm6N7xoCzmNYfpCF1/4i8zN8AVU7YBDkIPwP2TiBidFK74EcSMnvk4OdaEtexJNg7ltGpccP/E6F+46F2q5FgB7hOy4HQMnJsR82bkqDRKI4MDGaiKwlIyO2OX3krpXMTKn7bNEiqVvl8sulrqWqKuCZZ6TFTxPNO5GHsw8ZIv0MGBB+GRElx4+Wj/bKK9JM2b7i6apVkPc29/lMz1D5ujopb0jVrrBw3X2Bw+c5bxBZAIMgsicLr0htefJnn5MDXHZZfO8vKPAmWMuBUWWld/V3NWZ69q2we/WSHnv29D++HHS43d5go6xMyhmS84cqK4HNmyPPhTR9OlBR4d1ffk0eJaf0exol723uoA8xqzzTkwMkJ04DKgZCDG7IRhgEkT1xUjfjqDWhoZ5BbKNGwMCBgMslPS8pka4hVNAhJ23LfIOd3Fxg3DhppFtga1KoawkcJRdNhNFYc0s7YNZzHf2SoEMtqxFRopOM6j1CkyhBDIKI7MJOazfJgZRvC4z8OHKk1DVWWalepdrQIK3f9cAD3q64ggJlI8rkYEdemHX6dGn/wNYqOUCQ51HyvSb5mqNdS4TRWPXtQ48Ck5/X1yv4HMK18igNbjhCk6xG1Qwjk2NiNNmaXsngeia8Rhsxptas276fndKE53DPr7469PvVvhYt7ne4YyotO2eQJo0wMZqIzCGwtUBJC1S8rVRyS4zcygKoNxmfb+vGmjXSY1mZ93V5W7TzyC07vjlDo0dLrVVut/e9ak8sqFaXoZJWHqVl57xBZDWqhlQmx5Ygsh0j/vJWOlQ80nviPaeaLR/RWjcCWzpCDQGPd304M03hEGsLldKym+kayfLYEkREwfTKwYi27pYWQp3T93yJ5gP5tm6sWSMlM8+cKT2fN09az+yss4AVK4BRo4ITvjMzgddeA7Zvjz7TtZlztLRa+oIjNMkCGAQRWZleazcpWXcrsAsFiH+kkNstjbJaty74tcD1xOIV6vxXXy09zpvnXY1UTpIOd4xoI8hCldNMAUKsXVhKy84RmmQBDIKIrEyvHIxQwVagwBaoo0fjnyDQ7ZYCoMDV5pUGeGqPlJOHzIc6Vjx5S1YOEKxcdqIADIKIrMLIIfChgq1oAcr27VIQtHQpcPx4fK1UgfPoKA3wIg31Djfbsdy68d//SsPkq6qA/ful1xctkmaVlic2DHxvqJmu4w1Ejb7PZmmhItIBgyAiq4g2U6/eFVi0AEXuCuvVSwooACA7Ozg4kCv9UaOk+XoA7fKOtm+XPsP8/PBLPpSUSIFM4AzQ0brhAme6jpeRMzKzlYcchkEQkV1EqsDUbF2IFGyFS6D+8Ufp9y+/DF5KQ6703e7IS08MHBi9CyxaDtKhQ9Lv8mMoBQVSkBRPsrMegaidJsUkMhiDICIzU2sZAjVbF0KNkpIr/mgJ1Js2AVOmhD7u6NHeBORwXWyRyh7t3JMnA+3aSb/v3u39HEN1b8WT7Cy/N56WlFjuM9fuIlKNSwghjC6EXmpqapCeno7q6mqkpaUZXRyi6IqLgytjX0pHSFVWAnl50iKeWk5cFzgBoTzs/OhR4KmnpDUcRo3ydo+1bRs64KmqAoYNi628gYGEfMzycmkV93CmTwcWLAh9rEjJzmoGILHcZ73uJZGJaFV/W6ol6KOPPsJjjz2GiooKuN1ulJWVYdSoUUYXi0g7iQyBN2IxS99jlpRIj74tKHPnSj+hBLbaxCJcF1F2NvD115Hf+8svUoDh+161k52jiXafGzUKvodcmJQoYZYKgo4dO4azzz4bN910E6655hqji0OkvUSGwBu9mOXo0VLwcPXVQOfOUkvQzJnS82gtQY0axVaxh+siOnRI6oKTR7KVlUlBmVwOwNvqFKp7yTfZOdIw+URFu8+hWoq4MClRwiwVBA0fPhzDhw9XvH9tbS1qa2s9z2tqarQoFpE5aTWRotLE3PbtpUBoxQrgllv8X2vbNrjiz8315uMUFCir1OWy5Of7b5fzlDIypOfySLaNG71lk88V+N5Avq1BWicjy9cTSL6XVVXAM88A77yjzaSYRA7TyOgCaGn+/PlIT0/3/GRnZxtdJKL4xTryKDPT25ogV5a+zxMJguTRXJG89ZYUAAHA4sXS47x5Uj5LXl7oyl7pseV916yR9t+wQdpWWQl88IG0PT/fO8/PmjVSICEPta+uDv1e+cf3/HJrkB7kEXKTJwcna+fmSsHjO+9I29S4l0QOZ6mWoFg98MADmD59uud5TU0NAyGyLqvM4SLnIuXnS91O8+YBY8YAy5YBd98NDB/ubQlyu6Wk5enTvc+V8u3uk/OOQs1kDUgJ2r527fJfJyzwvUVFUuATS06VmkPXOfydSBe2DoJSUlKQkpJidDGIjJfo/DWxJFmHykVatkx6fOopoGVL/5FO8qzSbrfyYCMwyPrTn4CXX5aed+ggtfKMGQMcOyYFN7fcIrVGjRsnLXqakeHdJgdoM2cC/fp5A7RYc6riHboe7bNt1Ch4EsmRI6WuscpKJkUTJULVNel1BECUlZXF9J7q6moBQFRXV2tTKCK7KioSAgj/U1Tk3ffAASEqKqSf0lLp9ZkzpcelS6XXZRUVkY8beGwlZZk8WXpcuVI6X7Tjx3IdpaXebb7X4XstFRXqfrYDB8b2+RDZkFb1t61bgohIJbEkWYdqmejXT2o5kVdmD2zxkGdjlmdojpT0G6oscovQ0qVSK8+iRVI+UqQlLPr3B9avl7ronnpKeq9cvnDXEWppkESnIVAyPN63JUitBHcislZ32C+//IKvvvrK83zv3r3Ytm0bWrdujc6dOxtYMiKTSzRfJdHV6tu29XYdhRru7TuXULRjhypLr15SMnFGhjcZ+vTTgf/9X6lr7MwzpfmJAoOlYcOkHKWWLaUAKNbPRo1pCOL5bLWcs4jIQSwVBG3ZsgWDBw/2PJeTnidOnIgXX3zRoFIRWYBRSy2EykWK1PLhO0NzLNatC1701DcZum1b6bFnT+mxVy//16KtuTZqVOicKq2mISAiXVgqCBo0aBCEc1b5IDKnRJOsI7V8RDt2YIuWvP+oUcHrjj32mLRg6+jR0vZFi6TWH9/jR7sO3+Ax3FphibSQhTpepDLF+9lz0dVg/EwIsG5idDyYGE2OojSx98ABKbk2MNFXDZGShQ8c8CYxK00mVpJ8LO8jJ0XL1xrPNcaS7BxvYrQezFw2pdT+ntrhM3EQJkYTUWyU5qsY1VUWbmJALcQyx1K8yc6JtpBRZEZ9T8nWGAQR2ZVR+SqxBhHRuiNiOZ7bLc2fM3ly8HtCnTuUeJOdzTaZpREL6JodPxMKwCCIyK4i5avIlYHSCQpjES2ImDxZCnyUnldJUFJQIO139Kg0+WKkfaMFKkYnO6uVq2L0ArpqUDtoscNnQqpyCeGcTOOamhqkp6ejuroaaWlpRheHSD+VldJ6XRUVUmUeapi6r0Qqg8CKyzeIKCmJPHePzzI3uPde71Ia4Y4HePfJywNWrvSOBAs1h1Csw+ADPzc9qHVOJZ+b2Vs91P6e2uEzcSit6m+2BBE5QWC+itza4TskXa3WjkgtUMXFwaO4fM9bVSXN3QMA48d7jxV4vM2bgcsv9+8Gk739tnQO+Zi+w+LtUMEpbSlSe+SaEdRulbPDZ0KqYhBE5ARyvorb7Q1EAv/j16MyiFYJ+XZ1RLJokTT0PbCrZMMGqdUnI8O77+7dwcdW+he/XsnOseY9OSVBmEELaYxBEJGTyBXo6adLLSO+Fe2uXdKjml0CSoKIqirggw+AQ4e8AQsAlJVJZcrIAHJyvOWaPDn8shjyzNNTpwZviyf3Q69kZ61zVThyLRg/EwJzgoicRc43iUSv5FC5WydUMnOg6dOl7jHAfzLEzEzg8GHgq6+k9b9kd98NnHYa0KqVdJ777jN37ke0XBUl64eZ6Xq0wMkNHU2r+ptBEJHdhapgoy1YqqSSUatScruB7du9LUFyy41cxowMYNWq6IFSOHLLkZ7JzYkIlRgdLUF4+nQgNZUBAtkWE6OJKD6hulpiWbA0HLlrTc4zUmMZjcpKb9muvtpbppyc4Jag0lIgO1taGmP/fuDRR73H9A2ggMgj0qwgWoKwnFDuhDwhIhUxCCKyOy0WLPW1aJH2LRDhEmTffjt0C4kcSMlzCFkp9yNUropaCeVE5IdBEFEkdshDSGTB0kChRjEB6iVVZ2Z65wpScpxQAR4gzQkkD4kPTG7W657Gex6lydhVVcEjyTj7MVFsVF2JzOS4gCrFzG6LLCZ6PUVF0vvD/UyerM1CrL7CLaQpX1u0Muh1T7U6j3z906dHvhdFReqel8hAXECViBKX6LDgggLvwqehLFoUvSUj0ZaYaMe3cqudEr5zPoXKk9JreQ8iG2AQROZjdBeUnRdZTHTeG/n9BQXxzzat1WR/kQI8ve6pnt8dtScSNPrfHZERVG1XMjl2h1mE0V1Q0bp89OxmCNf1YwbyfYr1Xhlxf/W6p0rPo/Z9VeMzNfrfHVEE7A4jZ5D/GjWS0auI+zLzEgm+szdHY3Trml73VOl51Lqv8r+Xiy4CBg6UJlUkIsUYBJE5yJXkrl3eStWoLiiuV6SM3DWm5N5ovSxENHrdU72/O3IwtXQpsG6dd1bpWN5v165fIgUYBJE5GF1JmonZKqZIuSJKc4zM1LpmFLPdV4D/7sjxGASR8dxuID9f+ms23LIJOTnGlM2IRRbNVjFt3y6VJz9fnVmhZWq0kMSTzKvXPQ08j1r3VQ6mqqqADRukbatWSY9lZdL2tm2VBVUMTsnhuHYYGS/aukhO+2s02mKaercYvPKKNAps6VLvkOxEhFobywzH0ppa9zXavxdZrP9urPRZkuNw7TCyr0iz/sqLZ7rdzvmrVMu8EqUtJ74V9u7d3ke5+yaRQEytlhgzJNHHQq37Kv97kVuC5s0D/vQn4OWXpdbTfv28LUFEFBGDIDJeqMph8mRgyBCponv8cakFgv+pJ853VBIQPiD629+CV22fN8/bVTl9OrBgQXxlSHSuIjMk0Rs5p47v9bVtK92ToUOlIMh30dl4jmulNdaIVMDxlGROnLBNomXF5LsKvJWUlEjdNr4Lv06aJG3Ly9OndSjRz86MAYfvaD8ih2BLEJmLvICmvDikWUbRGCXRVhMg/Kikpk2l36uqgvcHgJUrpZaGsjKptWHmTKmlQS6XEdRKojd6dmQ17qt8nKIioHdv8wVVRBbAIIjMJTMTSE0Fhg3z385hu/GLNirpmWekYAeQPn/fLsjcXO8K8T17Gp8wG+paAG8wVFQEXHZZ9OPEM1mhGYe4+wZTffroe24iG2B3GJlPQYE0QqWiQho9A0iP8raCAmPLZzS327uAphK+n+fIkcGvv/NO5K6kjAz/RyOF+m4AUsuQ1t8NuRsuL88bROrdDUdEqmJLEJkPZ2yOLNZWDN/P8847paBn6VLg+HGpEp85U2pJmTsXOPfc0K0c06dL3Uxm6EYKl0SvZOHWcC05VVXSKMR77w1/HDXm1DH68yMiPwyCiKwk0WHhcrdXr17ebT17So+Fhf77BnZBZmZKlb/Z1jJTGlBE6xYEIo9CVCM4902o1iMJmUEXUUQMgsjczDiKxghqDQuXP89GjYDPP5e2yfMAycnFcsKx2WcOjvW7EaklZ9cu/9FmWlu0SJ/AxMwL8BKZAIMgMje1RtFYnVpLLsif55//7J0HSE4qlh9lciuHHIC53eZJCJbPG8t3I1Q55RFyx49Lj0qvK5YALFQ3HOBNOHfSaEcik+GyGeRMVuomcLul9bsOHYo8LDyW6/jgA2kE3ujR0pDz++4LbhWRl0+w47Im8hIRkah1XdE+v8mT1e0aM9uyK0Qq4LIZRGqyUjeBWsPCfcm5QStWSIEQ4G35CWzlsOMim/J8VEOHSp+FltdVUCB93+RuzECLFqnb4mm2BXiJTIxBEJHZRVpbrVev0JV1qJaucN0ye/dKj/KkiYEVsl6j9fRsncvMDL3shxbXJX+eBQX+uUdaBZJ2DFqJNMIgiJzDqMnuEq3c4xkWHqqlK1yLkjwqbOlSb7eangGJfK78fOu0zsVKz2kfOMUEkWKcLJGcw6jJ7rRYoyue4KSgQAqewlm61PsZhCuzFqP15HMdOqTeMZWc03fCSb1GIWZmRr4HRKQrtgSRc9ihmyBSZa2kpSvRbhm1R+vJw/4B71B9vVrnfFud4r2uWFvMfBcp1eP7xikmiCJiEETG06vrRc9uAq263iJV1koTYsN9BnoOhZfPU1LiTRiWE72tlMQbT4K9ntM+cIoJoogYBJHxrDRSSykjRujE0tIld8v4jljSs8zh8pNkkyerHxSbcQFUIjIUgyCyv1AtTVp3ExjR9RZLS1eobhk9yxzqXPIaZkuXKlsLLFZqBXkMpohsg0EQGUPPiiRUS5PW3QRWGKFj1FD4cOeS1zALN+w/UWoFeZyHh8g2GASRMViRaMuKCbEZGdqWWa0gzw4J9kQEgEEQGUXrikSvliYlSd1GBCSJtnTpWWb5XDk5sc98bQQrtPIRkSJcO4yMJ6/jJK9VpQa91rvSoux6stIaampS67qtfv+JLIJrhxHFgl0WythxZJ4SauWEWbHbkYg8GASR8bSoSLTssuDoIJJxHh4iS2MQRMazWkVi9aRupUGcU7vKiMgxuHYY2Z/aLU0FBVIOSEWF1MUGSI/ytoICdc6jFaVrqGmx5hnFL3C9MyJKGFuCyP7Ubmmy+ugg5ktZk1Pzt4g0xCCIyGkiBXF6rh+mFLvliEgjDIKIEmG30UFmzHdycgsIk/CJNMUgiCgRVkvqDhQYxLGrzFzMGJQS2QiDICInM3L9sEjYAiJhUEqkKQZBRGS+vBu2gEjMEpQS2RSDICIKnXdjZL4TW0CISAcMgojMxEwtMkbmO7EFJJjdkvCJTIBBEJGZ6DkSink31mL1JHwiE2IQRORUVsm7YQsIEWmEQRCR0YxqkbFK3g1bQIhIIwyCiIxmVIsM826IyOEYBBEZzSotMkRENsMgiMhoZmiRYd4NETlQI6MLEKtnn30W3bp1Q5MmTZCXl4f169cbXSQi65PzbhgEEZGDWCoIeuONNzBt2jTMmDEDW7duRf/+/TF8+HDs27fP6KIRqYMtMkREunEJIYTRhVDq/PPPR25uLhYuXOjZ1qtXL4waNQrz588P2r+2tha1tbWe5zU1NcjOzkZ1dTXS0tJ0KTMRERElpqamBunp6arX35ZpCaqrq0NFRQWGDh3qt33o0KHYuHFjyPfMnz8f6enpnp/s7Gw9ikpkfW631D0mD90nIrIhywRBhw4dQn19Pdq3b++3vX379vjxxx9DvueBBx5AdXW152f//v16FJXI+uSZqxkEEZGNWW50mMvl8nsuhAjaJktJSUFKSooexSIiIiKLsUwQlJGRgcaNGwe1+hw8eDCodYiI4sC1xIjIYSzTHZacnIy8vDysXr3ab/vq1atx4YUXGlQqIhspKQHy8qQfecbqSZO820pKjC0fEZHKLNMSBADTp0/HhAkT0LdvX+Tn52PRokXYt28fbrvtNqOLRmR9nLmaiBzGUkHQ2LFj8dNPP2HOnDlwu90488wz8d5776FLly5GF43I+swwczURkY4sFQQBwB133IE77rjD6GIQERGRxVkmJ4goLpzvJj6cuZqIHIBBENkb57uJD9cSIyIHYBBEREREjmS5nCCiqDjfjXrcbmlofEEBPzMish22BJH9cL4b9bA7kYhsjC1BZD+c74aIiBRgEET2w/luEsPuRCJyCAZBROSvpETqAvMldysC0tD54mJdi0REpAUGQWRvnO8mduxOJCKHYBBE9ibPd0PKsTuRiByCo8OIiIjIkRgEEVF47E4kIhtjdxgRhcfuRCKyMbYEEdkRF44lIoqKQRCRHXGmZyKiqBgEERERkSMxJ4jILjjTMxFRTBgEEdkFZ3omIooJgyAiu+BMz0REMWEQRGQXnOmZiCgmTIwmIiIiR2IQRGRHnOmZiCgqdocR2RFneiYiiootQURERORIDIKIiIjIkRgEERERkSMxCCIiIiJHYhBEREREjsQgiIiIiByJQRARERE5EoMgIiIiciQGQURERORIDIKIiIjIkRgEERERkSMxCCIiIiJHYhBEREREjsQgiIiIiByJQRARERE5EoMgIiIiciQGQeRsbjdQXCw9EhGRozAIImdzu4HZsxkEERE5EIMgIiIicqRTjC4Ake7cbm/LT2Wl/yMAZGZKP0REZGsMgsh5SkqkLjBfkyZ5fy8qkvKEiIjI1hgEkfMUFABXXin9XlkpBUClpUBurrSNrUBERI7AIIicJ1R3V26uNwiyC7dbavUqKGBgR0QUAhOjieyKI9+IiCJiEETOlpkp5QCxpYSIyHHYHUbOlplpryRojnwjIlKMQRCRnXDkGxGRYgyCiOyEI9+IiBRjEERkJ04Z+UZEpAImRhMREZEjMQgisiuOfCMiiojdYUR2ZbeRb0REKmNLEBERETkSgyAiIiJyJAZBRERE5EgMgoiIiMiRGAQRERGRIzEIIiIiIkdiEERERESOxCCIiIiIHMkyQdBDDz2ECy+8EM2aNUPLli2NLg4RERFZnGWCoLq6Ovzxj3/E7bffbnRRiIiIyAYss2zG7NmzAQAvvvii4vfU1taitrbW87y6uhoAUFNTo2rZiIiISDtyvS2EUPW4lgmC4jF//nxP8OQrOzvbgNIQERFRIn766Sekp6erdjxbB0EPPPAApk+f7nl+5MgRdOnSBfv27VP1QzS7mpoaZGdnY//+/UhLSzO6OLrhdfO6nYDXzet2gurqanTu3BmtW7dW9biGBkHFxcUhW2p8bd68GX379o3r+CkpKUhJSQnanp6e7qgvjywtLY3X7SC8bmfhdTuLU6+7USN1U5kNDYKmTJmC6667LuI+Xbt21acwRERE5CiGBkEZGRnIyMgwsghERETkUJbJCdq3bx9+/vln7Nu3D/X19di2bRsA4LTTTkOLFi0UHSMlJQVFRUUhu8jsjNfN63YCXjev2wl43epet0uoPd5MIzfeeCNeeumloO1r167FoEGD9C8QERERWZplgiAiIiIiNVlmxmgiIiIiNTEIIiIiIkdiEERERESOxCCIiIiIHMn2QdBDDz2ECy+8EM2aNUPLli0VvUcIgeLiYnTs2BFNmzbFoEGD8Pnnn2tbUJUdPnwYEyZMQHp6OtLT0zFhwgQcOXIk4ntuvPFGuFwuv58LLrhAnwLH6dlnn0W3bt3QpEkT5OXlYf369RH3X7duHfLy8tCkSROceuqpeO6553Qqqbpiue7y8vKg++pyubB7924dS5y4jz76CCNHjkTHjh3hcrnw1ltvRX2PHe53rNdth/s9f/58nHvuuUhNTUW7du0watQo7NmzJ+r7rH6/47luO9zvhQsXIicnxzMLdn5+Pt5///2I71HrXts+CKqrq8Mf//hH3H777Yrf8+ijj+Lxxx/H008/jc2bN6NDhw649NJLcfToUQ1Lqq7rr78e27Ztw8qVK7Fy5Ups27YNEyZMiPq+YcOGwe12e37ee+89HUobnzfeeAPTpk3DjBkzsHXrVvTv3x/Dhw/Hvn37Qu6/d+9ejBgxAv3798fWrVvx4IMP4q677sLy5ct1LnliYr1u2Z49e/zu7emnn65TidVx7NgxnH322Xj66acV7W+X+x3rdcusfL/XrVuHO++8E5s2bcLq1atx8uRJDB06FMeOHQv7Hjvc73iuW2bl+52VlYWHH34YW7ZswZYtWzBkyBBcddVVYRsfVL3XwiGWLFki0tPTo+7X0NAgOnToIB5++GHPtt9++02kp6eL5557TsMSqmfnzp0CgNi0aZNn2yeffCIAiN27d4d938SJE8VVV12lQwnVcd5554nbbrvNb1vPnj3FX/7yl5D7/8///I/o2bOn37aCggJxwQUXaFZGLcR63WvXrhUAxOHDh3UonT4AiLKysoj72OV++1Jy3Xa83wcPHhQAxLp168LuY8f7reS67Xi/hRCiVatW4vnnnw/5mpr32vYtQbHau3cvfvzxRwwdOtSzLSUlBQMHDsTGjRsNLJlyn3zyCdLT03H++ed7tl1wwQVIT0+Peg3l5eVo164dunfvjkmTJuHgwYNaFzcudXV1qKio8LtPADB06NCw1/jJJ58E7X/ZZZdhy5YtOHHihGZlVVM81y0755xzkJmZiYsvvhhr167VspimYIf7nQg73e/q6moAiLiCuB3vt5LrltnlftfX1+P111/HsWPHkJ+fH3IfNe81g6AAP/74IwCgffv2ftvbt2/vec3sfvzxR7Rr1y5oe7t27SJew/Dhw/HKK69gzZo1WLBgATZv3owhQ4agtrZWy+LG5dChQ6ivr4/pPv34448h9z958iQOHTqkWVnVFM91Z2ZmYtGiRVi+fDlWrFiBHj164OKLL8ZHH32kR5ENY4f7HQ+73W8hBKZPn45+/frhzDPPDLuf3e630uu2y/3esWMHWrRogZSUFNx2220oKyvDGWecEXJfNe+1ZdYO81VcXIzZs2dH3Gfz5s3o27dv3OdwuVx+z4UQQdv0pvS6geDyA9GvYezYsZ7fzzzzTPTt2xddunTBu+++i9GjR8dZam3Fep9C7R9qu9nFct09evRAjx49PM/z8/Oxf/9+/O1vf8OAAQM0LafR7HK/Y2G3+z1lyhRs374dGzZsiLqvne630uu2y/3u0aMHtm3bhiNHjmD58uWYOHEi1q1bFzYQUuteWzIImjJlCq677rqI+3Tt2jWuY3fo0AGAFGlmZmZ6th88eDAo8tSb0uvevn07/vvf/wa9VlVVFdM1ZGZmokuXLvjyyy9jLqvWMjIy0Lhx46DWj0j3qUOHDiH3P+WUU9CmTRvNyqqmeK47lAsuuABLly5Vu3imYof7rRar3u+pU6fi7bffxkcffYSsrKyI+9rpfsdy3aFY8X4nJyfjtNNOAwD07dsXmzdvxlNPPYWSkpKgfdW815YMgjIyMpCRkaHJsbt164YOHTpg9erVOOeccwBIeRjr1q3DI488osk5lVJ63fn5+aiursZnn32G8847DwDw6aeforq6GhdeeKHi8/3000/Yv3+/XzBoFsnJycjLy8Pq1atx9dVXe7avXr0aV111Vcj35Ofn45133vHbtmrVKvTt2xdJSUmallct8Vx3KFu3bjXlfVWTHe63Wqx2v4UQmDp1KsrKylBeXo5u3bpFfY8d7nc81x2K1e53KEKIsKkYqt7rmFOpLea7774TW7duFbNnzxYtWrQQW7duFVu3bhVHjx717NOjRw+xYsUKz/OHH35YpKenixUrVogdO3aIcePGiczMTFFTU2PEJcRl2LBhIicnR3zyySfik08+EWeddZa44oor/Pbxve6jR4+KP//5z2Ljxo1i7969Yu3atSI/P1906tTJtNf9+uuvi6SkJLF48WKxc+dOMW3aNNG8eXPx7bffCiGE+Mtf/iImTJjg2f+bb74RzZo1E/fcc4/YuXOnWLx4sUhKShLLli0z6hLiEut1P/HEE6KsrEx88cUX4v/9v/8n/vKXvwgAYvny5UZdQlyOHj3q+fcLQDz++ONi69at4rvvvhNC2Pd+x3rddrjft99+u0hPTxfl5eXC7XZ7fn799VfPPna83/Fctx3u9wMPPCA++ugjsXfvXrF9+3bx4IMPikaNGolVq1YJIbS917YPgiZOnCgABP2sXbvWsw8AsWTJEs/zhoYGUVRUJDp06CBSUlLEgAEDxI4dO/QvfAJ++uknMX78eJGamipSU1PF+PHjg4ZQ+l73r7/+KoYOHSratm0rkpKSROfOncXEiRPFvn379C98DJ555hnRpUsXkZycLHJzc/2Gkk6cOFEMHDjQb//y8nJxzjnniOTkZNG1a1excOFCnUusjliu+5FHHhF/+MMfRJMmTUSrVq1Ev379xLvvvmtAqRMjDwUO/Jk4caIQwr73O9brtsP9DnW9gf9P2/F+x3PddrjfN998s+f/s7Zt24qLL77YEwAJoe29dgnxezYRERERkYNwiDwRERE5EoMgIiIiciQGQURERORIDIKIiIjIkRgEERERkSMxCCIiIiJHYhBEREREjsQgiIiIiByJQRARERE5EoMgIrKc1157DU2aNMEPP/zg2XbrrbciJycH1dXVBpaMiKyEy2YQkeUIIdCnTx/0798fTz/9NGbPno3nn38emzZtQqdOnYwuHhFZxClGF4CIKFYulwsPPfQQxowZg44dO+Kpp57C+vXrPQHQ1VdfjfLyclx88cVYtmyZwaUlIrNiSxARWVZubi4+//xzrFq1CgMHDvRsX7t2LX755Re89NJLDIKIKCzmBBGRJX3wwQfYvXs36uvr0b59e7/XBg8ejNTUVINKRkRWwSCIiCynsrISf/zjH1FSUoLLLrsMhYWFRheJiCyIOUFEZCnffvstLr/8cvzlL3/BhAkTcMYZZ+Dcc89FRUUF8vLyjC4eEVkIW4KIyDJ+/vlnDB8+HFdeeSUefPBBAEBeXh5GjhyJGTNmGFw6IrIatgQRkWW0bt0au3btCtr+r3/9y4DSEJHVcXQYEdnOZZddhsrKShw7dgytW7dGWVkZzj33XKOLRUQmwyCIiIiIHIk5QURERORIDIKIiIjIkRgEERERkSMxCCIiIiJHYhBEREREjsQgiIiIiByJQRARERE5EoMgIiIiciQGQURERORIDIKIiIjIkRgEERERkSP9f5L9yLHE24nDAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using Distributions, PyPlot\n",
    "N = 250; p_apple = 0.7; Σ = [0.2 0.1; 0.1 0.3]\n",
    "p_given_apple = MvNormal([1.0, 1.0], Σ) # p(X|y=apple)\n",
    "p_given_peach = MvNormal([1.7, 2.5], Σ) # p(X|y=peach)\n",
    "X = Matrix{Float64}(undef,2,N); y = Vector{Bool}(undef,N) # true corresponds to apple\n",
    "for n=1:N\n",
    "    y[n] = (rand() < p_apple) # Apple or peach?\n",
    "    X[:,n] = y[n] ? rand(p_given_apple) : rand(p_given_peach) # Sample features\n",
    "end\n",
    "X_apples = X[:,findall(y)]'; X_peaches = X[:,findall(.!y)]' # Sort features on class\n",
    "x_test = [2.3; 1.5] # Features of 'new' data point\n",
    "\n",
    "function plot_fruit_dataset()\n",
    "    # Plot the data set and x_test\n",
    "    plot(X_apples[:,1], X_apples[:,2], \"r+\")   # apples\n",
    "    plot(X_peaches[:,1], X_peaches[:,2], \"bx\") # peaches\n",
    "    plot(x_test[1], x_test[2], \"ko\")           # 'new' unlabelled data point\n",
    "    legend([\"Apples\"; \"Peaches\"; \"Apple or peach?\"], loc=2)\n",
    "    xlabel(L\"x_1\"); ylabel(L\"x_2\"); xlim([-1,3]); ylim([-1,4])\n",
    "end\n",
    "plot_fruit_dataset();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Generative Classification Problem Statement\n",
    "\n",
    "- Given is a data set  $D = \\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$\n",
    "  - inputs $x_n \\in \\mathbb{R}^M$ are called **features**.\n",
    "  - outputs $y_n \\in \\mathcal{C}_k$, with $k=1,\\ldots,K$; The **discrete** targets $\\mathcal{C}_k$ are called **classes**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- We will again use the 1-of-$K$ notation for the discrete classes. Define the binary **class selection variable**\n",
    "$$\n",
    "y_{nk} = \\begin{cases} 1 & \\text{if  } \\, y_n \\in \\mathcal{C}_k\\\\\n",
    "0 & \\text{otherwise} \\end{cases}\n",
    "$$\n",
    "  - (Hence, the notations $y_{nk}=1$ and $y_n \\in \\mathcal{C}_k$ mean the same thing.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  The plan for generative classification: build a model for the joint pdf $p(x,y)= p(x|y)p(y)$ and use Bayes to infer the posterior class probabilities \n",
    "\n",
    "$$\n",
    "p(y|x) = \\frac{p(x|y) p(y)}{\\sum_{y^\\prime} p(x|y^\\prime) p(y^\\prime)} \\propto p(x|y)\\,p(y)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  1 - Model specification \n",
    "\n",
    "##### Likelihood\n",
    "\n",
    "- Assume Gaussian **class-conditional distributions** with **equal covariance matrix** across the classes,\n",
    " $$\n",
    " p(x_n|\\mathcal{C}_{k}) = \\mathcal{N}(x_n|\\mu_k,\\Sigma)\n",
    " $$\n",
    "with notational shorthand: $\\mathcal{C}_{k} \\triangleq (y_n \\in \\mathcal{C}_{k})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "##### Prior\n",
    "\n",
    "- We use a categorical distribution for the class labels $y_{nk}$: \n",
    "$$p(\\mathcal{C}_{k}) = \\pi_k$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Hence, using the one-hot coding formulation for $y_{nk}$, the generative model $p(x_n,y_n)$ can be written as\n",
    "\n",
    "$$\\begin{align*}\n",
    " p(x_n,y_n) &= \\prod_{k=1}^K p(x_n,y_{nk}=1)^{y_{nk}} \\\\\n",
    "   &= \\prod_{k=1}^K \\left( \\pi_k \\cdot\\mathcal{N}(x_n|\\mu_k,\\Sigma)\\right)^{y_{nk}}\n",
    "\\end{align*}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- We will refer to this model as the **Gaussian-Categorical Model** (<a id=\"GCM\">GCM</a>). \n",
    "  - N.B. In the literature, this model (with possibly unequal $\\Sigma_k$ across classes) is often called the Gaussian Discriminant Analysis  model and the special case with equal covariance matrices $\\Sigma_k=\\Sigma$ is also called Linear Discriminant Analysis. We think these names are a bit unfortunate as it may lead to confusion with the [discriminative method for classification](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/Discriminative-Classification.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- As usual, once the model has been specified, the rest (inference for parameters and model prediction) through straight probability theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computing the log-likelihood\n",
    "\n",
    "- The <a id=\"generative-classification-llh\">log-likelihood</a> given the full data set $D=\\{(x_n,y_n), n=1,2,\\ldots,N\\}$ is then\n",
    "$$\\begin{align*}\n",
    "\\log\\, p(D|\\theta) &\\stackrel{\\text{IID}}{=} \\sum_n \\log \\prod_k p(x_n,y_{nk}=1\\,|\\,\\theta)^{y_{nk}}  \\\\\n",
    "  &=  \\sum_{n,k} y_{nk} \\log p(x_n,y_{nk}=1\\,|\\,\\theta) \\\\\n",
    "     &=  \\sum_{n,k} y_{nk}  \\log p(x_n|y_{nk}=1)  +  \\sum_{n,k} y_{nk} \\log p(y_{nk}=1) \\\\\n",
    "   &=  \\sum_{n,k} y_{nk}  \\log\\mathcal{N}(x_n|\\mu_k,\\Sigma)  +  \\sum_{n,k} y_{nk} \\log \\pi_k \\\\\n",
    "   &=  \\sum_{n,k} y_{nk} \\underbrace{ \\log\\mathcal{N}(x_n|\\mu_k,\\Sigma) }_{ \\text{see Gaussian est.} } + \\underbrace{ \\sum_k m_k \\log \\pi_k }_{ \\text{see multinomial est.} } \n",
    "\\end{align*}$$\n",
    "where we used $m_k \\triangleq \\sum_n y_{nk}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2 -  Parameter Inference for Classification\n",
    "\n",
    "- We'll do Maximum Likelihood estimation for $\\theta = \\{ \\pi_k, \\mu_k, \\Sigma \\}$ from data $D$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Recall (from the previous slide) the log-likelihood (LLH)\n",
    "\n",
    "$$\n",
    "\\log\\, p(D|\\theta) =  \\sum_{n,k} y_{nk} \\underbrace{ \\log\\mathcal{N}(x_n|\\mu_k,\\Sigma) }_{ \\text{Gaussian} } + \\underbrace{ \\sum_k m_k \\log \\pi_k }_{ \\text{multinomial} } \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Maximization of the LLH for the GDA model breaks down into\n",
    "  -  **Gaussian density estimation** for parameters $\\mu_k, \\Sigma$, since the first term contains exactly the log-likelihood for MVG density estimation. We've already done this, see the [Gaussian distribution lesson](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#ML-for-Gaussian).\n",
    "  - **Multinomial density estimation** for class priors $\\pi_k$, since the second term holds exactly the log-likelihood for multinomial density estimation, see the [Multinomial distribution lesson](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Multinomial-Distribution.ipynb#ML-for-multinomial). \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " - The ML for multinomial class prior (we've done this before!)\n",
    "$$\\begin{align*}   \n",
    "\\hat \\pi_k = \\frac{m_k}{N} \n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Now group the data into separate classes and do MVG ML estimation for class-conditional parameters (we've done this before as well):\n",
    "$$\\begin{align*}\n",
    " \\hat \\mu_k &= \\frac{ \\sum_n y_{nk} x_n} { \\sum_n y_{nk} } = \\frac{1}{m_k} \\sum_n y_{nk} x_n \\\\\n",
    " \\hat \\Sigma  &= \\frac{1}{N} \\sum_{n,k} y_{nk} (x_n-\\hat \\mu_k)(x_n-\\hat \\mu_k)^T \\\\\n",
    "  &= \\sum_k \\hat \\pi_k \\cdot \\underbrace{ \\left( \\frac{1}{m_k} \\sum_{n} y_{nk} (x_n-\\hat \\mu_k)(x_n-\\hat \\mu_k)^T  \\right) }_{ \\text{class-cond. variance} } \\\\\n",
    "  &= \\sum_k \\hat \\pi_k \\cdot \\hat \\Sigma_k\n",
    "\\end{align*}$$\n",
    "where $\\hat \\pi_k$, $\\hat{\\mu}_k$ and $\\hat{\\Sigma}_k$ are the sample proportion, sample mean and sample variance for the $k$th class, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Note that the binary class selection variable $y_{nk}$ groups data from the same class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  3 - Application: Class prediction for new Data\n",
    "\n",
    "-  Let's apply the trained model to predict the class for given a 'new' input $x_\\bullet$:\n",
    "$$\\begin{align*}\n",
    "p(\\mathcal{C}_k|x_\\bullet,D ) &= \\int p(\\mathcal{C}_k|x_\\bullet,\\theta ) \\underbrace{p(\\theta|D)}_{\\text{ML: }\\delta(\\theta - \\hat{\\theta})} \\mathrm{d}\\theta \\\\\n",
    "&= p(\\mathcal{C}_k|x_\\bullet,\\hat{\\theta} ) \\\\\n",
    "&\\propto p(\\mathcal{C}_k)\\,p(x_\\bullet|\\mathcal{C}_k) \\\\\n",
    "&= \\hat{\\pi}_k \\cdot \\mathcal{N}(x_\\bullet | \\hat{\\mu}_k, \\hat{\\Sigma}) \\\\\n",
    "  &\\propto \\hat{\\pi}_k \\exp \\left\\{ { - {\\frac{1}{2}}(x_\\bullet - \\hat{\\mu}_k )^T \\hat{\\Sigma}^{ - 1} (x_\\bullet - \\hat{\\mu}_k )} \\right\\}\\\\\n",
    "  &=\\exp \\Big\\{ \\underbrace{-\\frac{1}{2}x_\\bullet^T \\hat{\\Sigma}^{ - 1} x_\\bullet}_{\\text{not a function of }k} + \\hat{\\mu}_k^T \\hat{\\Sigma}^{-1} x_\\bullet - {\\frac{1}{2}}\\hat{\\mu}_k^T \\hat{\\Sigma}^{ - 1} \\hat{\\mu}_k  + \\log \\hat{\\pi}_k  \\Big\\}  \\\\\n",
    "  &\\propto  \\frac{1}{Z}\\exp\\{\\beta_k^T x_\\bullet + \\gamma_k\\} \\\\\n",
    "  &\\triangleq \\sigma\\left( \\beta_k^T x_\\bullet + \\gamma_k\\right)\n",
    "\\end{align*}$$\n",
    "where \n",
    "$\\sigma(a_k) \\triangleq \\frac{\\exp(a_k)}{\\sum_{k^\\prime}\\exp(a_{k^\\prime})}$ is called a [**softmax**](https://en.wikipedia.org/wiki/Softmax_function) (a.k.a. **normalized exponential**) function, and\n",
    "$$\\begin{align*}\n",
    "\\beta_k &= \\hat{\\Sigma}^{-1} \\hat{\\mu}_k \\\\\n",
    "\\gamma_k &= - \\frac{1}{2} \\hat{\\mu}_k^T \\hat{\\Sigma}^{-1} \\hat{\\mu}_k  + \\log \\hat{\\pi}_k \\\\\n",
    "Z &= \\sum_{k^\\prime}\\exp\\{\\beta_{k^\\prime}^T x_\\bullet + \\gamma_{k^\\prime}\\}\\,. \\quad \\text{(normalization constant)} \n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- The softmax function is a smooth approximation to the max-function. Note that we did not a priori specify a softmax posterior, but rather it followed from applying Bayes rule to the prior and likelihood assumptions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Note the following properties of the softmax function $\\sigma(a_k)$:\n",
    "  - $\\sigma(a_k)$ is monotonicaly ascending function and hence it preserves the order of $a_k$. That is, if $a_j>a_k$ then $\\sigma(a_j)>\\sigma(a_k)$. \n",
    "  - $\\sigma(a_k)$ is always a proper probability distribution, since $\\sigma(a_k)>0$ and $\\sum_k \\sigma(a_k) = 1$. \n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Discrimination Boundaries\n",
    "\n",
    "-  The class log-posterior $\\log p(\\mathcal{C}_k|x) \\propto \\beta_k^T x + \\gamma_k$ is a linear function of the input features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Thus, the contours of equal probability (**discriminant functions**) are lines (hyperplanes) in the feature space\n",
    "$$\n",
    "\\log \\frac{{p(\\mathcal{C}_k|x,\\theta )}}{{p(\\mathcal{C}_j|x,\\theta )}} = \\beta_{kj}^T x + \\gamma_{kj} = 0\n",
    "$$\n",
    "where we defined $\\beta_{kj} \\triangleq \\beta_k - \\beta_j$ and similarly for $\\gamma_{kj}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  How to classify a new input $x_\\bullet$? The Bayesian answer is a posterior distribution $ p(\\mathcal{C}_k|x_\\bullet)$. If you must choose, then the class with maximum posterior class probability\n",
    "$$\\begin{align*}\n",
    "k^* &= \\arg\\max_k p(\\mathcal{C}_k|x_\\bullet) \\\\\n",
    "  &= \\arg\\max_k \\left( \\beta _k^T x_\\bullet + \\gamma_k \\right)\n",
    "\\end{align*}$$\n",
    "is an appealing decision. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <a id=\"code-generative-classification-example\">Code Example</a>:  Working out the \"apple or peach\" example problem\n",
    "\n",
    "We'll apply the above results to solve the \"apple or peach\" example problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p(apple|x=x∙) = 0.7996875392706204\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAG2CAYAAAB4e1KRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpZUlEQVR4nO3dd3hUVf4/8PekUtIzpBASCIggYoAASlyQJgg2QKwLiC7xh66wKroqIAqCgm3Vr6wYCasirmWNEV0Uy0pbBQ1txRIsBNImZdJpqff3x+VOy52aO+XOvF/PwxMyc3Pn3Bn0vnPO55yjEQRBABEREVGACfJ2A4iIiIi8gSGIiIiIAhJDEBEREQUkhiAiIiIKSAxBREREFJAYgoiIiCggMQQRERFRQGIIIiIiooDEEEREREQBiSGIiIiIApJqQ9DatWuh0Whw7733erspREREpEKqDEEFBQV49dVXkZGR4e2mEBERkUqpLgSdPHkSc+bMwcaNGxEbG+vt5hAREZFKhXi7Ac66++67cdVVV+Hyyy/HmjVrbB7b3NyM5uZmw/cdHR2ora1FfHw8NBqNu5tKREREChAEAU1NTejduzeCgpTrv1FVCHrnnXdw8OBBFBQUOHT82rVrsWrVKje3ioiIiDyhpKQEffr0Uex8qglBJSUluOeee/D555+jW7duDv3M0qVLsWTJEsP3DQ0NSEtLw9tvlyA5OcpdTSUiIiIFnTzZiKuvTkVkZKSi51VNCDpw4ACqqqowcuRIw2Pt7e3YvXs31q9fj+bmZgQHB5v9THh4OMLDwzudq0ePKEREMAQRERGpgSCIX5UuZVFNCJo8eTKOHDli9tjtt9+OwYMH46GHHuoUgIiIiIhsUU0IioyMxNChQ80e69mzJ+Lj4zs9TkRERGSP6qbIExERESlBNT1Bcnbu3OmGswoA2gC0u+HcRPaEAuDQLhGRJ6g6BCmvBYAOwcGnAQBcSog8SRCAjg4NBKEPgAhvN4eIyO8xBBl0QKMpQvfuwYiL643Q0DAATEHkOYIgoL6+Gg0NpRCEgWCPEBGRezEEGbQgKKgDCQmp6Nath7cbQwEqJqYXmpqOo729FQxBRETuxcJoExoNoNHwLSHv0Wg0HIYlIvIQ3vGJiIgoIDEEERERUUBiCCKHPfHESmRlDfd2M4iIiBTBEOQmmgodQp9YCU2FzmOvuW/fN4iKCsbMmdM89ppERERqxRDkJpoKHcLWrvJoCHrzzX/gzjsXY+/e/6KkpNhjr0tERKRGDEF+4tSpU/jgg/eQnX0Xpk27Glu2vG54bvfunYiI0GD79m0YM2YY4uO7YcKES/DDD8YNabdseR0pKTH4+OMPMXz4+YiP74ZrrpmC0tISm6/75puvITPzAsTHd8OIEYPx6qsvG55raWnBkiWLMGBAMuLju2HIkH549tm1il87ERGRKxiCFKSp0CHo8EHDHwBm37uzVygv710MHDgI558/CDffPBdbtrwGQRDMjlm+/K948slnsWtXAXr1SsBNN12L1tZWw/OnT5/GM888gZycN/Dll1+jqakRt912s9XXfO21jVi1ajkee+wJHDjwM1aufBJr1qzAW2+9AQDYsOH/8MknH2Hz5vdw6NBR5OZuQVpaP7dcPxERkbO4WKKCQjblIGztKrPHwhfdYfh7y9LH0Lp8pVtee/PmTbjpprkAgClTpuHkyZPYufM/mDjxcsMxS5c+hkmTpgAAcnLewKBBffDRR/mYPftGAEBrayuee249Ro++xHDMyJEXYP/+7zBq1MWdXvOpp1bjySefw4wZ1wEA+vVLR2HhT/jHP3IwZ858lJQUY8CAgbj00rHQaDRIS+vrlmsnIiJyBUOQgtoWLET7VdcCEHuAwhfdgeb1G9ExPBMAICQlu+V1f/nlKPbv/w5vvfUBACAkJASzZ9+EzZv/YRaCLrkky/D3uLg4DBw4CEeP/mx4LCQkBJmZowzfDxo0GDExMTh69OdOIai6uhqlpSW4++4FWLzYGPTa2toQFRUNAJg79zZce+0UjBgxCJdfPg3Tp1+NyZOnKnvxRERELmIIUpCQlNwp6HQMzzSEIHfZvHkT2tracP75Kca2CAJCQ0NRV1dn82c1FssTW35v7bGOjg4AwPr1GzFq1CVmzwUHi9s9DB+eiR9+KMLnn3+KnTu/xK233ogJEy7HW2+979iFERERuRFDkMq1tbXhn//cjLVrn8OkSea9LHPnzsa7776FIUOGAgC++24fUlPTAAB1dXX47bdfcP75g83OdfDgfkOvzy+/HEV9fb3ZMZLExET07p2CoqJjuOmmOVbbFxUVheuvvwnXX38TZs68HjNnTkNtbS3i4uK6fO1ERERdwRDkJkJSMlqWPua2ITDJp5/+G/X1dbj11gWIjo42e27GjOuxefMmrFv3PABg3brHERcXj4SERDz++HLEx2txzTUzDceHhobigQcW45ln/g+hoaG4//5FuPjiMbL1QACwbNlK/PWvf0FkZBSmTp2O5uZmHDq0H/X1dVi8eAnWr38eiYnJyMgYjqCgIOTn/wuJiUmIiYlx19tBRETkMIYgNxGSkt1WBG1q8+ZNmDjx8k4BCABmzpyNZ599Ev/7nzhT7fHH1+HBB+/B77//iosuGob33vsIYWFhhuN79OiB++57CH/60x9RVlaKrKyx2LDhH1Zf+7bbstGjRw+88MIzWLHiQfTs2RNDhlyEu+++FwDQs2cEnn/+Kfz++68IDg5GZuZofPDBJwgK4qREIiLyPo1gOY/ajzU2NiI6OhpbtzYgJSXK4tmzCAkpQmpqOsLDu3mlfe6ye/dOXHnlRJSW1lnthdmy5XU89NC9KCur92jbyFxz81mUlBShrS0dgH/9OyQiclVTUyMmToxGQ0MDoqIs79+u46/kREREFJAYgoiIqJOcHCA3V/653FzxeSK1YwgKAJddNgEnTwo2C5Lnzr2NQ2FEZBAcDLzySucglJsrPn5uJQwiVWNhNBERdZKdLX595RXj91IAuvNO4/NEasYQREREskyD0KZNQGsrAxD5Fw6HERGRVdnZQGioGIBCQxmAyL8wBBERkVW5ucYA1NpqvViaSI0YgoiISJZpDdDeveJXuWJpIrViTRAREXUiVwQtVyxNpGbsCaIui4jQ4OOPP/R2M4hIQe3t8kXQ2dni4+3t3mkXkZIYghT0xBPAunXyz61bJz7vDgsX3oaICA0iIjSIiQnF0KH9sWzZAzh16pR7XpCI/N7ChdZ7erKzxeeJ1I4hSEHBwcCaNZ2D0Lp14uPuXFxsypRp+P13HX744RgefXQNNm58GcuWPeC+FyQiIlI5hiAFPfww8Mgj5kFICkCPPCI+7y7h4eFITExCnz6puPHGP+Kmm+bg3//+EIIg4Pnnn8bQof2h1XbHmDHDkJ//vuHn2tvb8ec/L8CFF6ZDq+2OESMG4e9/f7HT+Tdv/gdGjboQcXHhGDAgGUuWLDJ7vqZGj5tvnoVevXpg2LCB2LbtI7Pnf/75J1x33ZVITIxAenoisrPnQa/XG57Pz38fF198EbTa7khLi8fVV1/OniwiInIrhiCFmQahuDjPBCA53bp1R2trK1ategRvvvkaXnhhAwoKfsSiRfchO3su9uzZBQDo6OhA7959sHnze9i//yc8/PCjWLVqGfLy3jOca+PGDViy5G7cfvv/w7ffHsF7732EAQPOM3u9tWtX4brrbsS+fd9j6tQrsWDBHNTW1gIAKip0mDZtPDIyhmP37v348MPtqKqqxK233mh4/vbbb8Gtt/4JBw78jE8+2Ylrr70OgiB46N0iIqJApBEC6E7T2NiI6OhobN3agJSUKItnzyIkpAipqekID+/W5deKiwNaWoCwMOBcFnCbhQtvQ0NDPd5550MAwP793+G6667EZZdNxPbt/8a2bV/hkkuyDMfffXc2Tp8+jdde+6fs+e67725UVVXirbfEHqOBA1Mwd+7teOyxNbLHR0Ro8OCDj+DRR1cDAE6dOoWkpEh88MEnmDJlGlavfhT793+LrVs/M/xMWVkpBg1KxaFDR3Hq1EmMHTsSP/10HGlpfZV4S1SrufksSkqK0NaWDqDr/w6JiPxBU1MjJk6MRkNDA6KiLO/fruMUeTdYt84YgFpaxO/d3RP06af/RmJiBNra2tDa2oqrrpqBe+55AB9++D6uvXaK2bEtLS0YNmyE4fvc3Ffwxhu5KC4+gbNnz6ClpQUZGcMBAFVVVdDpyjFhwmSbrz90aIbh7z179kRkZCSqq6sAAIcPH8Du3TuQmBjR6eeKin7H5MlTMWHCZFxyyUWYPPkKTJ48FTNnXo/Y2FhX3w4iIiK7GIIUZlkDJH0PuDcIXXbZRLzwwgaEhoYiObk3QkNDUVDwLQDg/fe3oXfvFLPjw8PDAQB5ee/h4Yfvw5NPPodLLslCREQkXnjhGezfL/5s9+7dHXr90NBQs+81Gg06OjoAiENu06dfg9Wrn+r0c0lJyQgODsbHH3+Bffu+wX/+8zleeeUlrFq1HDt3fot+/dKdeyOIiIgcxBCkILkiaOmru4NQz549O9XpDB48BOHh4SgtLca4ceNlf+6bb/bgkksuxf/7f382PFZU9Lvh75GRkejbtx927vwPxo+f6FLbhg3LxNateejbtx9CQuT/yWk0GmRl/QFZWX/A0qWP4oIL+uLjj/OxePESl16TiIjIHoYgBbW3yxdBS997enGxyMhI/OUvD+Chh+5DR0cHsrLGoqmpEfv2fYOIiAjMmTMf/fufh7ff3owvv/wMffum4+2338TBgwXo29fYA7Ns2Urcc8+d6NUrAVOnTsfJk03Yu/dr3HXXYofasXDh3Xj99Y247bZbcO+9f0V8vBbHjv2G999/B+vXb8TBg/uxc+d/MHnyVPTqlYCCgm+h11dj0KAL3PXWEBERMQQpafly6895enaY5NFHV6NXrwQ8++xaHD9+DNHRMRg+PBMPPLAMAJCdfSeOHDmM+fNvgkajwfXX34I77vgzPv/8U8M55syZj7Nnz2L9+uexfPkDiI/XYubM6x1uQ3Jyb3z55ddYseIhzJx5BZqbm5Ga2hdTpkxDUFAQIiOj8PXXu/H3v7+ApqZGpKX1xZNPPoepU6cr/n4QERFJODvMQNnZYUSu4OwwIqLO3DU7jOsEEREFgJwc67u/5+aKzxMFGoYgIqIAEBws7v5uGYSk3eLdua0Pka9iTRARUQCQNkN95RXj91IAktstnigQMAQREQUI0yC0aRPQ2soARIGNw2FERAEkOxsIDRUDUGgoAxAFNoYgIqIAkptrDECtrdaLpYkCAUMQEVGAMK0B2rtX/CpXLE0UKFgTREQUAOSKoOWKpYkCCUMQEZEPyckRp6vLBZLcXHH7nYULnT9ve7t8EbT0vae39fFl7voMyPdwOMwN2tvbsXv3Trz33tvYvXsn2v3o/y5PPLESWVnDvd0Mn8L3hJTk6no+9hZDBKz39GRn276pB9pCi1xTKXAwBCls69YPMGRIP1x55UT86U9/xJVXTsSQIf2wdesHbn/tffu+QVRUMGbOnOb21yLnPPfcOowadSF69eqB4cPPx3vv/dPbTSIflZ3duVbHkfV83HnjDrRQ4OpnQOrD4TAFbd36AebOvR6W27GVl5dh7tzrsWXL+5gx4zq3vf6bb/4Dd965GG+8kYuSkmKkpqa57bU8SRAEtLe3IyREvf9cv/lmD5566nn0738e3nlnC+6441aMHj0G6en9vd008kGurOfjzsUQA3GhRa6pFBjYE6SQ9vZ2PPjgPZ0CEADDYw89dK/bhsZOnTqFDz54D9nZd2HatKuxZcvrZs/v3r0TEREabN++DWPGDEN8fDdMmHAJfvjhiOGYLVteR0pKDD7++EMMH34+4uO74ZprpqC0tMTma7/55mvIzLwA8fHdMGLEYLz66ss2j29ubsYDD/wF/folID6+G6ZMGYsDBwo6tfXLLz/DuHGjEBcXjq+/3tPpPCdOHEdEhAb/+tc7mDz5UsTHd8OoURdi9+6dZsf9/PNPuO66K5GYGIH09ERkZ8+DXq83PP/FF9sxZcpYpKTEIC0tHtdffzWOHfvd7BxlZaWYP/9mpKbGISGhJ8aNG4WCgm/Njnn77TcxZEg/9O4djfnzb0ZTU5Phuby8bZg8eSrS0/tj4cJFaG9vh05XbvN9osDmyno+pj0YWVnKhhR3ntsR3hiS45pK/o8hSCFff70HZWWlVp8XBAGlpSWyN3Ml5OW9i4EDB+H88wfh5pvnYsuW12QD2fLlf8WTTz6LXbsK0KtXAm666Vq0trYanj99+jSeeeYJ5OS8gS+//BpNTY247babrb7ua69txKpVy/HYY0/gwIGfsXLlk1izZgXeeusNqz/zyCMPYuvWPLz66hv4738Pon//8zBz5hWora3tdNyqVWtx4MDPGDo0w8b5/orFi+/H118fwpgxl+Kmm65FTU0NAKCiQodp08YjI2M4du/ejw8/3I6qqkrceuuNhp8/deoUFi1agl27CvDvf/8HQUFBuOWWWejo6AAAnDx5EtOmjYdOV4533/0Ie/f+D/fe+6DheQAoKvodH3/8If71r3/jX//6N/7731147rl1ndoqCAKWLbsfQ4YMxahRF1u9JiJX1/Nx543bm6HAG0NyXFPJ/6l3fMHHVFToFD3OWZs3b8JNN80FAEyZMg0nT57Ezp3/wcSJl5sdt3TpY5g0aQoAICfnDQwa1AcffZSP2bPFUNDa2ornnluP0aMvMRwzcuQF2L//O9mb9lNPrcaTTz5nGObr1y8dhYU/4R//yMGcOfM7HX/q1Cnk5m5ATs7rmDp1OgBg/fqNGDLkC2zevAn33vtXw7GPPPK4oa22LFy4CDNnzgYAvPDCBnzxxXZs3rwJ9933IDZu3IDhwzOxcuWThuM3bPgHBg1Kxa+//oKBA883/Kzk73/fhPT0BPz880+48MKheO+9f0Kvr8auXQWIi4sDAAwYcJ7Zz3R0dCAn53VERkYCAG65ZR527foPgCfMjrv77mx8++032LbtK4SFhdm9NgpMlkNN0veA/eAhd+NWKqy489z22BqSGz1a+ZlcXfkMyD5fmYHHniCFJCUlK3qcM3755Sj27/8O118v9tiEhIRg9uybsHnzPzode8klWYa/x8XFYeDAQTh69GfDYyEhIcjMHGX4ftCgwYiJiTE7RlJdXY3S0hLcffcCJCZGGP48/fSaTsNJkqKi39Ha2ooxY/5geCw0NBQjR17c6TVM22HLxRcbrykkJAQjRowynOvw4QPYvXuHWfsyMwcb2gIAx479jttv/yOGDu2P5OQoDB2aDgAoLS0GAHz//WFkZIwwBCA5aWn9DAEIED/n6uoqs2N++OF7bN78D7z77kfo3TvFoWujwGNtPR9HFjZ052KIvrDQotyQ3OjRQEGBsj1EXfkMyDG+UmzPniCF/OEP45CS0gfl5WWyw1AajQYpKX3whz+MU/y1N2/ehLa2Npx/vvHGKggCQkNDUVdXh9jYWJs/r9FobH5v7TFpOGj9+o0YNeoSs+eCrfwLlt6bzucTOj3Wo0dPm+22RTpXR0cHpk+/BqtXP9XpGCmQ3nDDNejTJxUvvbQRycm90dHRgYsvHoqWlhYAQPfu3e2+XmhoaKfXNx0uA4Djx4sAAOefP8j5C6KA4ep6Pq4shujob+O+tNBidraxUDk0FNiwoXMvTVeLtrmmkvv5SrE9Q5BCgoOD8fTTL2Lu3Ouh0WjMgpB0Q37qqReshgNXtbW14Z//3Iy1a5/DpElTzZ6bO3c23n33Ldx55yLDY999t88wa6yurg6//fYLzj9/sNn5Dh7cbxj6+uWXo6ivrzc7RpKYmIjevVNQVHQMN900x6H29u9/HsLCwrB373+RmvpHAOIQ3MGD+3H33fc6de2SgoJ9GDv2MkP7Dx8+gIULxWseNiwTW7fmoW/ffrKzy2pqanD06M/4v//LMQTUb775r9kxQ4dm4I03clFbW2uzN8iesWPHY/fuAvsHUkCzNQRg68bgyo1b+m3c8tymNyNXz+0utobklJrJ5epnQM7xhRl4DEEKmjHjOmzZ8j4efPAesyLplJQ+eOqpF9wyPf7TT/+N+vo63HrrAkRHR1u053ps3rzJLAStW/c44uLikZCQiMcfX474eC2uuWam4fnQ0FA88MBiPPPM/yE0NBT3378IF188xmoR77JlK/HXv/4FkZFRmDp1Opqbm3Ho0H7U19dh8eIlnY7v2bMnsrPvwvLlf0VsbBz69EnDCy88jTNnTuPWWxe49B68+urfMWDAQAwadAHWr38e9fV1mDfvTwCAhQvvxuuvb8Rtt92Ce+/9K+LjtTh27De8//47WL9+I2JjYxEXF4/XXnsVSUnJKCkpxmOPPWx2/htuuAXPPvskbr55JlatWoukpGT873+HkJzc22x40Z7du3fgsceW4tChQpeuk8gWV27cjv427iuhwF6djmkPEcOKOnj7c2MIUtiMGdfh6qtn4Ouv96CiQoekpGT84Q/jFO8BkmzevAkTJ17eKQABwMyZs/Hss0/i8OGDhscef3wdHnzwHvz++6+46KJheO+9j8wKdHv06IH77nsIf/rTH1FWVoqsrLHYsKFzbZHkttuy0aNHD7zwwjNYseJB9OzZE0OGXGSzV+fxx9eho6MD2dnzcPJkEzIzR+HDDz+zO2xnzapV6/D880/hf/87hPT0AXjnna3QarUAgOTk3vjyy6+xYsVDmDnzCjQ3NyM1tS+mTJmGoKAgaDQavPHGO3jggb/g4ouHYuDAQXjmmf/D9OkTDOcPCwvD1q2fY+nS+zF79pVoa2vD4MFD8Le//d2pdjY2NuDXX4+6dI1E7uILv407wt6Q3IED3ivaJtd5s9geADSCXAGLn2psbER0dDS2bm1ASkqUxbNnERJShNTUdISHd/NK+9xp9+6duPLKiSgtrUNMTIzsMVu2vI6HHroXZWX1Hm2bq06cOI4LL0zHN98cQkbGcG83RxHNzWdRUlKEtrZ0AP7375B8V1aW8Wa0d68y51RyBpCtc911l1gcbdlD5Ithjoys9ezJfW5NTY2YODEaDQ0NiIqyvH+7jj1BREQBxjJQWP42vnChMosPOlpz5AhrYSk31zwAmb4Wp7T7Ll8ptmcIIiIKMKbhBDAPJNLQkhLDEp6YAeRLRdvkOF/53DgcZuDfw2GkDoE8HOYri6cFCtOiYtMAZPl3JYKK9FpSTxOHqchZ7hoOU9ViiRs2bEBGRgaioqIQFRWFrKwsfPrpp95uFhEpwFcWTwsU2dnAyJHi3zdtMg890sKASv02zj24yFepajisT58+WLduHc47T9yy4I033sCMGTNw6NAhXHjhhV0+vyBAdqFDIs8REKj/BH1l8bRAkpNjXhBt+h4r+X57ewYQkTWqCkHXXHON2fdPPPEENmzYgH379ikQgkIhCEBz82l062Z/hWAid2htbTn3t8Ds9lDLdG1/4Ylwwj24yJepKgSZam9vx7/+9S+cOnUKWVnyC9Y1NzejubnZ8H1jY6ONMwajoyMGer2431N4eA/ZrSKI3EUQOlBbW4329h5Q8X+aXebtxdMChSfCia/MACKyRnX/pz1y5AiysrJw9uxZREREID8/H0OGDJE9du3atVi1apUTZ09CaytQWVkFjQZgBiJPEgSgoyMIQBqAwP3Hx6ET9/NUOPGVGUBE1qhudlhLSwuKi4tRX1+PvLw85ObmYteuXbJBSK4nKDU11crsMFPtAFqVbzyRXWFQ2XwFRTmzeBq5jjPxSG3cNTtMdSHI0uWXX44BAwYgx4GVvWxPkScib7IWeBiEiIgrRlshCIJZbw8RqROHTojI01QVgpYtW4bp06cjNTUVTU1NeOedd7Bz505s377d200joi7ylZ3KiShwqCoEVVZWYt68edDpdIiOjkZGRga2b9+OKVOmeLtpREREpDKqCkGbNm3ydhOIiPwWC6Yp0ATuNBQiIjLDrUuMcnI6vw+S3FzxeVI/hiAiIgJg3DPMNAgF6uy8QA+EgRICVTUcRkRE7sWtS0SBvpedFAIB60tW+AOGICIiMsOtS0SBHAgDJQRyOIyIiMzIbV0SqLKzje9DoAVC0+HRrCz/C0AAQxAREZkw/W1/797ONUKBJtADob+HQIYgIiICYH1j1UANQgyE/h8CWRNEREQAuHWJKWuBEJAvGPZH1jY0Bvzn2hmCiIjIY9SyIGOgB8JACYEMQUREBMAz06LVMvU60PeyC5QQyBBERORhvtob4olp0YEy9VrtAiUEMgQREXmYL/eGeGJtnEBef4d8C2eHERF5mK9vT+GJadH+PvWa1IEhiIjIC3x5ITpPTIv296nXpA4MQUREXuKLvSGeWBvH8jVGjrT+Gv60WSf5HoYgIiIv8VZviLUdwqVwMnKk+xZLlBv2Gz1a/Gr5GoGyY7uvCoSd5BmCiIi8wJurEUuF2ZavVVAgfpVCiUQKQkpMi5abei2d37QNtmqk1HBzVkMb7bH278SfwilnhxEReZi3F6KzNk39wAHrdUlKtcfa1GvTNmVl2Z4x5suz6yRqaKM9gbCcAUMQEZGHObsQnTvWFfLFaerZ2ca22KqRUsPN2dU2+toaUr7470RJDEFERB7m7EJ07upVcDR0WHLXjVquRsqRIOSrN2dX2uiLPUiu/jtRA9YEEZFX+EPNhKe4a10hVwuz3VEr4kqNlC/OrrPkbBt9cQ0pf17OgD1BROQVvvgbry9TuuejKzuEKz0c5WqNlDM9R97iSht9qZfL33eSZwgiIq9QQ12Hr1FqWEKJwmwlb9SubNaphptzV4Omt4egvF3A7wkMQUTkNb70G68aKNXzodQO4UrdqJ2tkVLDzbmrbfSFXq5A2EmeIYiIvMoXfuNVAyV7PpTaIdxbN2o13Jy70kZf6eUKhJ3kGYKIyKu8/Ruvr01JluOLPR/evFGr4ebsaht98bP2ZwxBROQ1vvAbrxoKtH2t54M3avfxtc/a3zEEEZFX+MqNVA0F2r7W88Ebtfv42mft7xiCiMgr3H0jdWaYiwXazuGNmvwFF0skIq9YuND2asBdrcNxdkE/NSy8p4QQvQ7JOSsRotd5uykGalo4U01tJfsYgojILzm78q4/r4prKlSvQ++NqxDqQyFITbuVq6mtZB+Hw4jIbzk6zOULBdqBTA11WRI1tZXsYwgiIr9mbx0iXynQdoSr0/lD9DpDz0+PwoNmXwGgVZuMNm2yW9rsKDXVZamprWRbQA6HFRYWoLq6HB0dHd5uChF1kb0ajYULbQ9z2SrQvvNO35rp5OpQTK+8HAyZOxJD5o5EvzV3AAD6rbnD8FivPN8oZFFTXZaa2krWBWRP0EMPXQ4ACAkJRUKvPkhKTkViYhqSkox/EhNTkZSUhoiIKC+3lohssbfOD2B7mEtNM51cHYqpnr0Q9eOvBSD2APVbcweOP7IRpwdnAhB7gnyBtxfOdIaa2krWBWQI6hXVCzUna9HW1opyXRHKdUVWj43oGS2GoqRUs5AkBaWEhBSEhIR6sPVEZMpWMADUMczlDFeGYtrODXfl5ADhtUl4CcDpwZk4cy4EAd5fHVtNdVlqaivZFpAh6Lc3cxEarUVxaSWOlepQXF6B4ooKlFRWoqSmGmU1Vaisr0bjmUacPNWA334/gt9+PyJ7Lo1GA622N5KT0871JnXuVYqOjoNGo/HwVRIFDrlgMHIkMHq0fy7o5+p+a8HBwPq83kjAI5ht8ri3V8dWU12WmtpK9gVkCAKA7hHhGDQ4DYMGp8k+L3QIqNU3oqhEh6KScpyoqBRDUpUYlHS11ahqqEJreyuqq8tQXV0GYK/sucLDuyMpMQ1JyZ2H26S/h4d3c+PVEvk/y2Bga70Wtd+kXB2Kyc4Ggk414dE3V6Pisyb8abBvzGxS0wrUamor2acRBEHwdiM8pbGxEdHR0WjYuhVRKSldPl/L6RaUV9aiqKQcRWU6FOsqUFJVhZLqKpTVVkFXV426k7UOnSs2NgFJiamGoGTsTRLDUlxcIoKCArKOnfyMuzYslW7mUjBQ+qbuKxutWhuKceZ63f1eKcVX3nNnqbXdvqypqRETJ0ajoaEBUVHK1eoGbE+QEsJ6hKFfehL6pSdhoszzQoeAk41nUFSsQ1FpOU6UV+BERSVKq6tQqq9GeW0VKhqqcLblDOrqqlBXV4WfCw/IvlZISCgSE8SQZNqLZNqb1LNnpHsvmEgB7tiw1BM1Gr6w0apSQzGuDqd5mi+8565Qa7sDEUOQG2mCNIiM6YGMmAHIyBgge0x7azsqK+tQVKLD8TIdTujE+qRSfTVKa6qhq6uCvlGPtrZWlJUfQ1n5MauvFxkRg8SktHP1SZ2DklbbGyEh/MjJu5RebM5TNRqeXiQvRK9Dr7wcVM9eaFjDR6mhGLXMbFLrwoRqbXcg4nCYCpw52YwTpRUoKq3AiXNBqaSqEqU1erE3qb4KTWea7J4nKCgIWm1v2bok6U9UVCyLuMkjlBqS8fTQg6eGkroXHsSQuSPx05YDZrO4ukqJ4TRPU8vwnSW1ttsXuWs4jCHID3S0daC2thFFxTocLz3Xm3SuiLu0Rg9dnbGI255u3XqYFXFbBqWEhD4s4ibFZGUZeyT2ys8r8EmeaLc7QpC1wKOGIMR/K4GNNUFkVVBIELQJMdAmxGD0qAtkj2k53YLSihpxpluZDsUVlWJQqq5CWW01KurFIu6zZ0/j+IlCHD9RaPX1YmMTZBeWlP7ExSWwiJvsUsuQjCV3ttvd21s4Opzma4W9/LdC7sIQFCDCeoShf/9k9O8v/z9QoUNAU8Ppc0sC6HCiXIfiyiqUVFcairgr66twtvWssYj75/2y5woNDTu3Erf1JQFYxB3YvLnYXFdu8O5ud6+8HPTeuMrsMWmbCwAov+Mx6BaudPn8jq6O7enCXlufyV13AQUF6luYkAsqqgNDEAEQi7ijYntiWOx5GJZxnuwxbc1tqKqux7HichwvrzgXlM7NdjtXxF3TVIPW1hb7RdyRsUhMTEWy2ZIAqSziDgDeXmzO1Ru8J9rtK9tbeLqw19pnIgUg00Uv1bAwobf/jZPjeJchh4WEh6B3Hy1699FiLDJkjzl9shknSipwvPTcbLcKsYi7RF8NXZ0eFXWVOHn2JJqa6tDUVIfffvte9jxBQUHopU3ptF2JaY8Si7jVyduLzbl6g/dEu9tkhrsst7fwFFe251DitaTvc3ONAWjDBvnjfXVhQm//GyfHsTCaPKqjrQM1NY0oKik/V8RdKW5ZYuhNEou429rb7J6re7eeYihK7rxmEou4yR5fn7njrtlhzvJkYa+vfybkPZwdpgCGIHVoOd2CEp3eEJRKKitxoqICpTV6w75udafqHDpXXFyi1cUlWcQdeCzX3vHlmTty6wR5mjdCiS9/JiTPE4X0nB1GASOsRxgGDOiNAQN6yz4vdAhobDh1biVuHU6U6VBSVYXiKpMi7oYqNLc2o7a2ErW1lTaLuBMTUjsNu5kGpR49Itx5ueRBoXodem9chfrx1+KVD5MVm7njjptAmza5S0XQXeWNwl7OplInNa+QzRBEqqMJ0iA6NgLDYwdi+LCBsse0NbehsqoOx0rKcbysAid0OpRUntvXrUbc100q4i4t+x2lZb9bfb3IyFizfdyM+7qJj7GIW3025CfhlTzlbvBqvgnI8UZhL2dTqZeaV8jm/7nJL4WEhyAltRdSUnthHIbJHiMVcReVlON4ubjIZGl1FUr01SivrUZFfRVOmRRx//rr/2TPIxVx29rXjUXc3mO59s5qPIL1eb2xaHY57hpbgVZ9MrKzxeEmV2+6ar4JyPF0YS9nU6mfJwvplcSaICIrTIu4i0qMK3FLSwKU11ahqqEK7R327wjdu/UUQ5FMUEpMTGURtxsl56w0W3tnJR5DMNqxAmsAmK+909X6BelmHhQEdHTI3wS4i3hnvrY4I7nOXTVdLIxWAEMQKa35VDNKdXoUlZ7bsqSiAiWVlYYi7or6KtSfqnfoXHFxiTb3dYuN7cUibhdY9gTJrb2jZOGxdBMA1Lk9BZGr3FlIz8JoIh8U3jMcA85LwYDz5EO1VMR97MS52qQy3bl93Yy9SZZF3D/9VCB7LqmI29ZK3Czi7syTa+9YFvb6y/AYkT1qreliCCJyI6mIe0Ts+Rgx/HzZY9qa26CrrMWx4nKcKK/AifJzC0xK+7rVVUPfpHeoiDsqKs7mkgBabTKLuN3E2k1AbTUSRM5Sc00X/29I5GUh4SFITUtAalqC7PNCh4Azp1twvFgccjMUcVdVoURfhfI6vaGIu7GxFo2NtVaLuIODgtGrl+0i7sjIGL8t4m7VJqP8jscU337C3k1A6h3y1RsBUVeoeYVs1gQR+YGOtg5UV9Wf2/y2wlDEXVJVabYSt0NF3N0jkJSUhuRzQcl8SYA0JCSkICws3ANXpR7WCnsdKZQmIvtYE0REVgWFBCGxdxwSe8dhDC6UPUYq4j5WIi4wKRVxl+irUVYrrsRdf6oeZ86cRFHRTygq+snq68XHJVmtTZKKuP21N0mO3MwltdZIEAUShiCiAOFIEXdD/SkUFZejqFSH4vIKQxF3SU0VymvF3qTm1mbU1FagprYCP/74ney5wkLDkZCYimSzoGS64GQqunfv6c7L9So110iQPF/YxoSUxxBERADEIu6YuAiMiHOsiPt4mQ7FukpDEXdpTRUq6qpRc7IGLa3NKC39DaWlv1l9vfDwePTta7uIOzg42F2X61ZqqpFwZo2eQF7Px3TLFYYg/8EQREQOMy3iHo/hnZ4XOgScPiWuxP3k5ja89WUrxg45jLT4oyjRi7PdSvR6tLY3obm5Br/8UoNffjks+1rBQcHoldDH0HskN/QWERHtk8NutoKAr/UAObPlh79tD0LEEEREitEEadAzshuGDOmHLeuAQblJePSVuXj8zjK8lV2B1blJePSVFDw45wfMuvQgisoqUKzTma3EraurQlVDNdo72lFRcQIVFSesvl6PHpGGPdxMC7iloJSY2AehoWEefAd8gzNDN85s+eFv24PYY7nQpulXQPmFNsnzVDU7bO3atfjggw9QWFiI7t2749JLL8VTTz2FQYMGOfTznB1G5HlS8AkL7UBLaxAev7MMK7IrbP5M86lmlOj0KCrR4XhZOU7oKlFSVWHc162uCg2nG+y+tkajMRRxy21X4q9F3N0LD2LI3JH4acsBhxeFdGa1X3euDOxLLLdcsWS65Qq5F7fNADBt2jTcfPPNGD16NNra2rB8+XIcOXIEP/30E3r2tF9kyRBE5B3hWSPQ0hqEsNAONO891OXzSUXc4krc4rIAxRUVYm2SvhrlddWorK9ES1uL3XOFhXVD4rkibmMPUppJD5P6irhdCUGAc/s+uWuPKF/i6S1XyDpOkQewfft2s+9fe+01JCQk4MCBA7jsssu81CoismV1bpIhALW0BmF1bpLdniB7pCLuzLjzkTnCehF3eUUNikp04mw3aUmAcytx6+qqUHuyFi0tZ1FS8itKSn61+nrRUfGGJQHkepR8oYi7q0M3llt+5Oba7gly9Fg18+SWK+QdqgpBlhoaxO7wuLg42eebm5vR3Nxs+L6xsdEj7SIikTQUJg2BSd8D6HIQsickPARpfROR1jfRZhH38WIxJB0vE4NSaXXVuWG3KujqKnGm5QwaGmvQ0FiDo0fle7GCg0PQq1eKYUkAuaDk7iLuXnk5nYZu+q25w/B3W0M3zqxpxPWPyJ+oNgQJgoAlS5Zg7NixGDp0qOwxa9euxapV1sdzich9LAMQYAw+ngpCtkhF3BdemI4LL0yXPaa9tR3V1eJK3EVlFThRrkNJZSVK9dUo0YtLAlQ1VqO9vc2JIm75LUsSElK6VMRdPXsh6sdfK76WlaEbOc6saSR37J0zdYgo+BnPvjLJ7FhT/jC13l1brpB3qTYELVq0CN9//z3++9//Wj1m6dKlWLJkieH7xsZGpKameqJ5RIFLrwfy8tB++kHZImjp+/Z23y9GDg4NRlLveCT1jkcW5H/Zaj7VjOLyarE3qbQcxRWV4iKT+iqU14r7ujWebsDp0004duxHHDv2o+x5TIu4rQWlmBit1d4kV4dunFnTSO7YUL0OzxyYjLOzy9Dc3lv2Nfxhan2bNplF0H5IlSFo8eLF+Oijj7B792706dPH6nHh4eEID+ceR0QepdcDGzdi5ZbxwODBsod4swfIzLnAhtmzAa3WpVOE9wzHwIF9MHCg/P+LhA4BdTVN4ky3Ul2nfd3Ka6tR2VCF1rYW6Gt00Nfo8MMP38qey7yIWy4spaJbtx5Otd+ZNY1sHXvXrAqcGSwfggJtaj2ph6pCkCAIWLx4MfLz87Fz506kp8t3YRMROeRcYMP48S6HIHs0QRrE9YpCXK8ojMyUX86jrbkNZYYlASpwQqczrMRdXquHrq4KNU01DhdxJ/dKRt++gxD33t+RkH6BWVCKj0/qchG3K0XYpkFo0yb/nlpP6qGqEHT33Xfjn//8J7Zu3YrIyEhUVIi/TUZHR6N79+5ebh1RANPrxT8AUFho/hUQA4abQoY/CAkPQd9+SejbL0n2eaFDwKmTZ3G8RIeiEpMlAarEDXB1ddWdirgLAeDE0U7nsizittzXTVqJ2xZXi7Czs40BKDSUAYi8T1XrBFkbC3/ttddw22232f15rhNE5CY5OWKPijV33OHeyldnhrUsA9uaNcAjjxiH7lQa2KQibnFftwqc0ElLAlSitEZvKOLu6LC/cVnPnlGyNUnSY72DgtCjoQaAc+vnBMoii6Q8rhMEcTiMiHzQ7NnikBJgPVi4kzPDWnl5nQPbmjXGv7s7sLmJaRH3pbhI9hipiPtYSTmOl56b6VZyHBW//YTjmhDoGmvReKYRp0412i3i1sYni4GoRwSGAgg7cRTxMVoxLAWHIEYQzH5x9cbUen+YlUbupaoQREReYq+nRa73ZPBgq4XRXuXtwOZFskXchYXA3LnAli0Qzh+EupomsTepVIfjugqUVFahtKoSJef2datsqEZrWwuq9eWo1pfjCIDPAODNZ8U/54SFdTP0IDU2pqGwMBUTJqThoovScOJEGubOTQXQ3a1ByB9mpZF7MQQRkX0eKCB2mpU6pJV5FyE4SMCKOyo6tXV1bhLa25OxcqHO/FyeDmz2QqUCs9ZcYVrEPWqk+H6szEmGNkVA3rkZfaZF3C++34r6mhMY0vIRfu8ehdLGBrOVuIuLf0Fx8S+G8+/cKf6RxERrodWm4dNP01Bf33lZgPj4JAQFBbl8PZyVRvYwBBGRsrRacUjJ3TdvK8NawQAexWrg+HGsyDE+Zbp4o9fZC5XuDp1OFLIHBwtmi1tKRdybvxyOD7+WFsO8wfCjQoeAk41nxCLuUnG2W0ml6Urc1aioF4u46xv0APTQ6w/ihMw6kyEhoYYibuN+buZhKSLCdn0IZ6WRLQxBRCTP1RlfWq1nCi2sDGutGDwYyP8fHs2bAeSab9fRafFGTwU2X+NEXZTcKt9W30+IvUmRMT1wUcwAXHTRANmXtyziPl5eIQ65nSvi1tVWobpJj7a2Vuh0x6HTHbd6KRE9o8VglCRfyJ2QkILs7FBs2iSgtVWD0FAB2dm+v1AneQZDEBHJ8/UCYht1SCuWtgG9yvDoKylYsykZLa1BsjdsZwPbypxkBAcLsos9ikNtms5DbRJ7oVKjAaTJH86ETmdnxuXlARMmOFUXZRqEbL6fDnKkiPvMyWaU6oxF3OKSAFXnVuKuRkVdFRrPNOLkqQb89vsR/Pb7EdnzaDQa9OjeG62taQhCClpb++Gee9Iwa5Zxgcno6Hi37utGvoshiIjkqbyAeEV2heGGHRbaocgq1ZZDQxKHhtrshcrMTODgQevPWwudzgydmR5rWQNlpy7KHe+nLd0jHF+Ju6ik/FwRdyVKq6o6FXGfOl0GoAwd537266/FP5Lw8O5ISkyzumVJQkIfdOvGtej8EUMQEclT04wvmWGt1blJhht2S2sQVucmdfnG7ezQkIFeD5w+Dbz0EhAbKx8qLXuCfCx0uuP97ApHVuJ+7LmeePztWNx+6eeYHPE+qj//CIUXZOKrqhj8WnMS3UOP4UyrHs3NZ3Ci+ChOFHdeXFISG9NLDEVWglJcXGKXirjJOxiCiEj9LIa1LIOJ9D3Q9X3LXBoa0uuBt94Cpk83D5H2QqW1552p13JkGM5OXVRX388uDSN2gebor3gc32HFNyY9aj+LvW2r8Qjahl6MB56JQVGxDsfLxCLu4spz9Un6aujq9NDVV+JsyxnU1Vejrr4aPxcekH2tkJBQJPTqYzUkJSam2i3iJs9jCCIi+9xcQKzkTVKuZ0auB6crPD001Ikz9VpdrO1S4v3s0jBiF6x8MgzQDwawpVPv2goA0IYBMT2QETMAGRnWi7irqupw7FxQOqGrREllBUqqq1B6bthN31SDtrZWlOuKUK4rstoeqYjbWlDq1as3QkJC3fJekDyGICKyz80zvly6SVopCG5v18j2zEjft7d3vQDWoaEhWz0wdXXAnDnWQ6W90OlMvVYXa7uUeD9dHkbsKgWGdINDg5GcokVyihZ/sFHEXVxaiaIyMSgVV1SgtFKsTSqrqUJFfRWazjTZLeIOCgpCfHyy2ZIAlkEpOjqORdwKYggiIq9z6SZppSDYVo+R2XlcXJDQ4aEhR3pgbIUgW6HTmZt7F4OAw++nHUrPMPMl3SPCMWhwGgYNTpN9XugQUKtvNBRxn9BVoLiyUlw7qaYautpqVDVUobW9FdXVZaiuLgOwV/Zc3br1QGJCqqE3yTIoJSamIjy8mxuv1r8wBBGRT/D4TdKFBQmdGhrqQg+Mt2po3M2rw4heXBNKE6RBfEI04hOiDStxW2o924pSnV5cZLJMh+JzQamkugrlddXQ1VWj7mQtzp49bb+IOzYBSYmpJsNu4lIALOLujCGIiHxGp5vkzB+AHJPeGlcXcFSIU0NDXeiBcWp40Jmbu5cXh/TKDDPTHj8f3i01tFso0tOTkZ6ejIkyzwsdApoaTp/rTRJD0omKcwtM6qtRXluFyvoqnG09i7q6KtTVVdks4k5MSO20wKRpj1LPnpHuvWAfwRBERD6j001yYxJW5G009qgsW+baWjoSB0OUtZ6YlQt1WJ2bhJU5nfcfU/Jm7tTwoDP1Wp5azVuGO2fs2eSL+965QBOkQVRsTwyLPQ/DMs6TPaa9tR2VlXXiStzlFWJvUoV5EXfNuSLusvJjKCs/ZvX1IiNixFCU3HkVbn8q4mYIIiKfIH+THAbgEXEmj14vBqDVq4H0dNfW0nFwppTis5lc6IHxpxoaT8zYI7GIu3cfLXr30WIsMmSPOXOyGSdKK1BUcq6Iu7ISpVWVKNZXo/xcEffJsyfRdLIeTSfrbRZxa7W9kZSUaravm2lQUkMRN0MQEXmd2U1y5g9AoR4rxhYCPyfh0V2rgY3/wYrLzxWKxsQ4t9aOKQfrdBSfzeRiD4zXp+IrROkZe3ZrpprOYOX0b8UHvDBs6su6R4Rj8OC+GDy4r+zzHe0CamsaUFSsQ1Fp+bklASpRUiXu61ZeV4Xqhmq0treiqqoUVVWlsFXEnZTYeV83KST5QhE3QxAReZ3ZTTLH2FuzAgDwCNp3BQO7VokHv/eeuOpykfX1WKxyok7HF3pifG2VZlcpNcNMYrenbuRW4K255j/kS/ve+bCgYA20CTHQJsRg9KgLZI9pOd2C0ooaHC8pP7d20rl93aqrUHZuX7e6U3U4e/Y0jp8oxPEThbLnAc4VcScZC7eNvUqpHini1giCtE67/2tsbER0dDQatm5FVEqKt5tDRHJM63ZycoA9e6wfm5kJPPmka7/VFxYCc+cCW7bY7EkKzxphCCLNew85/zouslZDo9YhMaXZfH9m/mBe+yXX4xegPUGeYFnEfaJch+LKKtkibnukIu5eCb1x+PB/0dDQgKgo5VbeZk8QEYlcXDdHcaY3qBtvFEPQ6tVAc7OyNzMH6nS81RPDGhr7bPfUqWjfOz/kSBF3W3MbqqrrDUXcJ8p14rCbC0XcXcEQRESirs6icUeIio0Vv6anGx9T6mZmp07Ha7OZ4JlVr/2Bv9RMBaKQ8BCnirh/LjqO+199VvF2cLUkIlKGFKKkYQgleGldG2s9MY/fWYZHX0nB6twkt77+yoU6qzf0FTN/wEqsVPZ9Vim5nrpOvLw2ErlOKuKePmUMsm+5yi2vwZ4gokDm5cUH7bLsrfHQzcyne2L8ZN2brnK4p86LayOR72MIIgpkXdxh3KMhyoM3M6VnM7mNr9RxeRhrpkgpDEFEgayLO4x3OUQFAkeDiq3jrIXN8HDx/c/ICKgQ5NM9daQqLoWgM2fOoLa2FikW08x//PFHXHjhhYo0jIg8oIs7jHc5RPkipXtXfv3VsaBia5jLXtj86isgK6vrbVUJ1fTUkc9zujD6/fffx/nnn48rr7wSGRkZ+Pbbbw3PzZs3T9HGEZGP02qNoUkKPqbfqzUEKVngXV9v/tUVs2eL6xlt2QLcc4/42IIF4h8ASEsTQ2hhoe126/Xi2kssqlZce3s7du7fj7e3b8fO/fvR3t7u7SaRA5zuCVqzZg0OHjyIXr16Yf/+/Zg/fz6WL1+OP/7xjwigdReJ/A9n0SjHdPjq+HHjV2koS+qBc7SmyrTHLi9P/Lppk/G4F180/t3WECSLqt3ig6++wj3PPovSqirDY30SEvDiAw/gukmTvNgyssfpENTa2opevXoBAEaNGoXdu3fjuuuuw2+//ebzG6URkQ1dLTx2NkR5o6jX2msqXeD95pvAW2+ZP7ZpkzG4zJkD3HefazVVkyYB+fm2F5BUI5UWeX/w1Ve4/sEHYdkFUFZVhesffBDvP/00g5APczoEJSQk4Pvvv0dGhri4UXx8PL744gvMnz8f33//veINJCIfJHfDcjZEeaNXwtpreqvA25WaqoEDxfaMHm0MbrbquHx9GQSJCnup2tvbcc+zz3YKQAAgANAAuPe55zBj/HgEBwd7uHXkCIdDUFNTEyIjI/Hmm28iJMT8x8LCwvD2229j0aJFijeQiDzMkd/I1XbDkq4pQ35lWsULvOfNA6ZPF/++Y4fYA7RgATBxovn5XClMNw2bprU91j43zuBzmz2HDpkNgVkSAJRUVmLPoUOYMGqU5xpGDnM4BI0bNw7bt29Hnz59rB7zhz/8QZFGEZGb2ZuO7a6A441eCb0eKCgQr0kqJLZWd2OqK9tzmJ5P2u2+Xz/l964yHYK09rnZCnh1dcC+feLPeiPMqqWXygqdgwXmjh5HnudwCBo1ahQuueQSfPbZZxhs8h/yoUOHsHz5cnzyySduaSARuYErQUeJG5Y3eiVMX1OqyfFkT0hMjPlXa1wpTLfWK2R5jLWAV1go1i5Nn+6dsKHyXqpkB98zR48jz3M4BOXm5mLVqlUYO3YsPvzwQyQkJOCRRx5BXl4err32Wne2kYjcpa5OnDI9YQIgze60FnCUuGHJ9UqMGyfuFh8ba/9G7EzxrBTaMjLEHqBNm4CrrgK2bQNGjBDP0a9f5/MoPUtOquEZOND2ca4Upnu6J0Xp4mWVrzM1bsQI9ElIQFlVlWxdkAZAn8REjBsxwtNNIwc5VRj92GOPISwsDFOmTEF7ezuuuOIKFBQUIDMz013tIyKlyN0wDx8Ww8Fvv4m1K6YsA44SNyy5m/KePeLN35GhImd6sORC27Zt4tdDh4AzZ4AXXpAPQa70PlgLCO7c7sPZYKrVijPT6uqM6woBjgcnpYdKlR6G9LDg4GC8+MADuP7BB6EBzIKQNFf6hfvvZ1G0D3M4BOl0Oqxduxa5ubkYMmQICgsLcfPNNzMAEamF3A1TGh7asQOYNUu8gVsLOGq7YcmFNqlHSHpMyVoYbxSLOxtMtVqgRw9g8WLzx1U0BOVrrps0Ce8//XTndYISE/HC/fdzeryPczgE9e/fH4MHD8a//vUvXHXVVfjss89w4403orS0FA899JA720hESpBumHV1xh4gaXhowQJg+HDxuMRE8au7Ao7UI1VXJw6F7dljuyfC1SEf08fr6sSv8fHG1wSMRcs+XoBrlSvB1Nng5KkhNxUv1nndpEmYMX489hw6BJ1ej2StFuNGjGAPkAo4HIJee+013HzzzYbvr7jiCuzYsQNXX301Tpw4gZdfftktDSQihUg3q5wcY2+INDxkuvrwrFmOncvVG5azQzhK1CJ99ZX49emnzR9fsUL8OmuWeA5nr0eNs5ucDU6eKl5257ChBwQHB3MavAo5HIJMA5AkMzMT33zzDa688kpFG0VEbjR7trjX1IoVxuEh054Ajcb+zbsrdTOnTwMvvSQWQjvSE6FELdL11wM//2weUEzl57t2TV6e3bQyJxnBwYK4aahFMF2dm4T2do3NzUYdovLiZSJbnCqMltOvXz98/fXXSrSFSL3szZrxpS0BtFpxteE77jAuHmjZEzBokHteW683TsmW1qkBxCE4az0RStQiDRokFkHr9eIQmNQD1NWbuRIB4ehR4LnngPvvd/p9Dw4W8OgrKQCAFdkwBK7VuUl49JUUPH5nme0TONKjp7ZaMCIndDkEAUBsbKwSpyFSL3tFsb62wrLU62GtZ8RTlNhh3VHuuJkrcc5jx4CDB8WvToagFdkVAGAShCrMApD0vFVKDkH5UtAncpAiIYiIVMqdxajSTdHWGkRNTeJX04UEbd1MlWivViv2RH36qbEnSsVMg9CaTcloaQ1yLAC5wtb772tBn8gBDEHknzzxW6m9oliNxv4ChN6+WbijGNV0n66NG8Xv8/PNjzGtmwGAykrje1NXZ/1mqkR7tVrgyivFEKQkZwLa0aNizw8gblth+hUA+vd3qldoRXaFIQCFhXa4JwABqi9eJrLEEET+yRO/ldoris3MFIc5rD3vr+uxSO/96tXi95MmiWEUMNbNWDJ9TJqdlpfn2owtR0hD+NJXJUKzMwHhuec6/9vYts04Wy8zE3j1VYdfenVukiEAtbQGYXVukvuCkCk1zo4jMsEQROQqe0Wxlj1BgTCrRio8BoDjx8WvlZXGsCGtQbR6NZCebnxf7rnH+H5IN9X8fDEMpKcbb6ZdCSu2bthFRWJwy8jwTF3L/feb9wRt2yau2TRmjPhY//4On0qqAZowshGTRovDi6Y1QtIxiswUs6Tyvb+IGILIf3j6t1JXimL9dVaN9N7n5RmHvuQ2K5V6edLTzd+H4mLgxRc7n1eaxSXdTLvSw2fvhg2IBdru7kHU64GdO82D1rZtYgCaPt2pU1nOApP+/vidZYYgZPq44jh9nlSOIYj8B38r9R65996UtCWHtTWIJk0S/9TXiz1IUoBasEDc5DQmxvou6Y6yto1Gv37G15R6r4qKXA/NjiyXoETQ0uvR/l05Hp93GiuyGw0PmwahrwoisfNAlHsLpTl9nlSMIYj8hzd/K7VXFKviLQEcYmufrtWrxXWJpGs3LfiV3peBA23vbQaIQerCC42vYXoOR95XueNMz2/6/YoVxuBmbSjOWthxNuT07y8O+zkxBCa9zsqDc4ElWwCI/8ZNZ4pJ9UFuC0BEfoAhiPyHN38rtVcU6++zauTe+379xK9STY+1n5PeFylIWS5m+OOP4hCb9AdQrodv9WqxQNly9hpgfD1rQ3HOhB17Q7VPPqlYQPbYTDFLlkGf6waRCjAEEZF7xMQ41/slBSmtVuyFyc8XA+zYseKNtK4OeO89cfPTrvbwSTfs0aPFP5mZ1rcRceT89kLOp5+KK2WbciXIOVD3tvrDoe6dKWYt3FgGfa4bRCrAEET+yd+Hn3yZ6RBXVpZrPz97trF3RgpGhYXG3d+VWOnZ9Iadni5+lXqvpPNLoUOvN4aNHTvE3iqpfqiw0NhbZco05MyZA2zZYjxeek6aJefov1M7dW+rR36IRw8YV4uWCqcBKBeEGG7IjzAEkX/y9+EnX6bUgoaeDLEajdgbpNGI30trFNmrUwLMw87gwcAVV4gz3Sx7k+SuxXKWnD026t5W5w/Fo3nDzGqA5LbVsKmrQ1hcN4hUhiGIyF/4Uw2GFKRMe2Ckr+PGiUNjhYXK3VQFQawNuv1241Dc7NmOzSiTwo5UyzRnjni8ZW+VFBCkdZRMr0m6ZnvXYqPurT2ul2wRtPR9e7vG/vtgrZfH0XDDGZqkMgxBRP7CU8MUngxbcjfVPXuMw2JK31RjY+WH4kxNnCgGj8JCMQRZhp2vvpI/t5sDgq2FELs8FOZo27luEKkMQxAROUdulpS9UORqcLI2Y0yJm6pp70ZBgfh1xw7j89Jj9l5H6tkxrRmaOFHsrdLrjT+rdEBQasjQkV4eR9vOdYNIZRiCiNTMF2owHOmBcrWXyp03VXv1PqYrWN9xhxgETEOHVmucVWZpxw7xj2nvjtLXolTdm6O9PAw35IcYgojUzFM1GPb23XIHudc0fb2uBjzT3o2CAjH0LFggfr9pk7if2XnnicNbEyZ0Dh1arbi+z6+/2l/p2peHgdw1hMUZmqQCDEFEauapGgxH9t2y7IECXO+l0uuBZcs677QOdN5PzFXW6n0AMcyMHi3+XSqStnYOez1Kcu30pYDgbA+Vo23nDE1SAYYgIjXzVA2GXNiyZNkDdfq06wsE6vViALLcbd7RgKd08bY0ZV7uXK7ULak5IKi57UQWGIKI1MKbU+Dlwpa9gPLrr2IIWr0aaG52rZfKch0dRwOerane1lY7lno3amrEafJ1dUBlpfh8fr5Y/yMtbGj5s3IrXbsaRL39OftKDxWRBzAEEamFveJiT9/A7AUUaSgsPV0MFACQmNg5HEg3/QkTxPV6APfVHf36q/geZmRY3/IhL898nzKJvWE4y5WuXeXNFZnZy0MBhiGIyF/YuoEp2btgK2xZK6CuqRH/XlLSeSsN6aav19veeiIz0/4QmL0apPp68e/SVzmzZ4shyZViZ08EUX9aFJPIyxiCiHyZUlPglexdkJslJd347RVQHzkC3Hij/HknTTIWIFsbYrPVdnuvPWsWEBcn/v34ceP7KDe85Uqxs/SzrvSkOPM5c+8uIsUwBBH5MjVsQ2B647c27fz0aeDtt4GUFPHmLg2PxcYab/aVleL3gDhsBjhXW2NtptyBA+Iu7qa9TJs2GUPNnDnAfffJn8sdizTKUcPnTOSHVBWCdu/ejWeeeQYHDhyATqdDfn4+Zs6c6e1mEblPV6bAe2MhRdNz5uWJX017UHJzxT9yLHttnGFtiCgxESgttf2zZ84AOTnmP6t0sbM99j5njabzZ8iNSYm6TFUh6NSpUxg2bBhuv/12zLa2bgeRP+nKFHhv9y5MmiSGh4kTgaQksSdowQLjlhKAsSdI7qbvzI3d2hBRfb04BCfNZNuxQwxlUjsAsS2LF8sPL5kWO9uaJt9V9j7nnBz2FBG5gapC0PTp0zF9+nSHj29ubkZzc7Ph+8bGRnc0i8g3uWshRUcLc+PixKCxYwcwY4b5c7GxnW/8gwcb63Fmz3bspi61JSPD/HGpTikmRvxemsn2/ffi9/Hxxtey/FlLpr1B7i5Glq7HkvRZ1tUB770nbiDLjUmJuizI2w1wp7Vr1yI6OtrwJzU11dtNInKdszOPtFpjb4J0szT9vishSJrNZcvOncYNSbduFb9u2gTMnSv+kbvZO3pu6diCAvH4w4fFxwoLgb17xcczMozr/BQUAE88YZxqf/Kk/M9Kf0xfX+oN8gRphtysWZ2LtQcPFsPjnj3iY0p8lkQBTlU9Qc5aunQplixZYvi+sbGRQYjUSy1ruEi1SBkZ4rDTpk3A5MnAf/4D3HILcOmlxp4gvV4sWp4zx/i9o0yH+6S6I7mVrAHzzVABMQy9/LL1n5U2THWmpkrJqeuc/k7kEX4dgsLDwxEeHu7tZhB5X1fXr3GmyFquFuk//xG/vv02EBFhDHOFhcZVpfV6x8OGZci66ipg2zbxe60W+O47MXidPSuGmxkzxN6oK64APvtMHCaTHpMC2oIFwPDhxoDmbE2Vq1PX7b23Gk3nRSTHjROHxgoLWRRN1AUaQZD+61IXjUbj9OywxsZGREdHo2HrVkSlpLivcUT+Rq4w15RpILC8qa9ZYwwrq1eLG5NKN+3CQnFozBbLsGGvLVL9zksviYXR0hR3R9i7DlvrFknXsmWLc7PI7F1PZqb8RrJybSbyU41NTYieOBENDQ2IiopS7Lx+3RNERApxpsharmdi+HDxZi3tzG7Z4yGtxiyt0Gyr6FeuLaYhKyZGDEFffWV7C4sRI4BDh8QhurffNgY0W9chtzVIV5chcGR6vGlPkFIF7kSkrhB08uRJ/Pbbb4bvi4qKcPjwYcTFxSEtLc2LLSPycV2tV+nKVH1AHGKSeivkej5M1xKyd265tqSniz1AMTHGYui0NODBB8UC6AEDxPWJLMPS4sVijVJEhHkPlaOUWIbAlffWnWsWEQUQVYWg/fv3Y6K0tgdgKHqeP38+Xn/9dS+1ikgFvLXVglwtkq2eD9MVmp1x4EDnTU9Ni6GlqfL9+olf09ONz5kGNEumm7vK1VS5axkCIvIIVYWgCRMmQKUlTET+o6tF1rZ6Puyd27JHSzp+woTO+47dcw9QXCwu2giIASkmxvz89q7DNDxa2yusKz1kcuez1SZX33tuutoZ3xOCykIQETnB0XoVZ28GzkzVt9UDJbcwoL1zW57P1vFardgbJF3XHXcAAwea72LvawXF9q6/Kxu0qn3TVaVDiz+8J9RlDEFE/srRehVv3QysLQzoDs4GN1eKnbvaQ0a2MbSQGzAEEfkrb9WrOBsi7P1m78z59Hpx/ZxZszr/jNxry3G12NnXFrP0xga6vo7vCVlgCCLyV7bqVaSbgaMLFDrDXoiYNUsMPo6+riOhZPZs8bjTp8XFF20day+oeLvYWalhH29voKsEpUOLP7wnpCjVLpboCi6WSAHLciE/ZxY/dJatRQbz8myv3TNnjvHv8+YZa5bsLVqo14vX99JL4mwv02OtLdToCFcXQOwKpV7T2cUefZHS/0794T0JUFwskYhcZ1mvIvV2mE5JV6q3w1YP1MKFnWdxmb5uXZ24dg8ATJ9uPJfl+X78ERg71nwYTLJrl/ga0jlNp8X7ww3O0Z4ipWeueYPSvXL+8J6QohiCiAKBVK+i14u/XZuGBIknbgb2bkKmQx225OeLU98th0oOHxZ7faTeIEBchdry3I7+xu+pYmdn654CpUCYoYXcjCGIKJBIN9C0NLFnxPRGW1QkflVySMCREFFXB+zdK+7zJQUWANixQ2xTTIw4tV1ql7Q3mNy2GNLK008/3fkxV2o/PFXs7O5aFc5c64zvCYE1QUSBxZUNS91FGtaRK2a2NGeOODwGmC+GqNUCTU1ASYm4/5fklluA1FQgMlJ8nRdf9O3aD3u1Ko7sH+ZL1+MOXNwwoLEmiIhcIzfU4syGpbbO29X9yKQhujFjjD1BUs+N1MaYGGDfvs7hzXRbDEumgWjWLPGrLw+j2Bv2kSsQNu0pmjMH6NHDvwOCry1BQH6BIYjI38kNtTizYak10tCaXi/enJTYRqOw0Ni2iRONbRo4sHNP0COPAImJYi9QRQWwebPxnKYBCrA9I00N7BUISwXlgVAnRKQghiAif+eODUtN5ee7vwfCWk/Jrl3yU6ilICWtIaSm2g+5WhWlCsqJyAxDEJEt/lCHYOsG6mxxqNzQGqBcUbVWa1wryJHzyAU8QFwTSJoSbzmM4qnP1NXXcXTYp66u80wyrn5M5BSGICJb/H06srN1FnJDa4CxN2nWrK4Pjd13n/1jpOAmd6OfNcv2ooie+kzd9TrS9e/b1/XVsYkCHEMQUSDp6rTg2bONG5/Kyc+3H6yUKqi21UZ/DKwS04JyuTopT23vQeQHGILI93h7CMqfN1ns6gwb6ednz3Z9tWl395DIndNTn6kn/+0ovZCgt/+7I/IChiDyPd4egvKlTRZ98cbkq6v42gp4nvpMHX0dX/xcvf3fHZEXMASRb5FuDt7k7V3ETfnyjcl09WZ7vN275qnP1NHXUepzlf57GTYMyMwUF1UkIocxBJFvkG6SRUXGm6q3hqB8tafD10g9L458Nt7uXfPUZ+rpfztSmFq9Gjh40LiqtDM/769Dv0QOYAgi3+Dtm6Qv8bUbk62hG0drjHypd81bfO1zBfjfHQU8hiDyPr0eyMgQf5u1tm3CwIHeaZs3Nln0tRvTr7+K7cnIUGZVaIkSPSSu1NZ46jO1fB2lPlcpTNXVAYcPi4/t2yd+3bFDfDw21rFQxXBKAY4hiLzP2tozpqv+ZmV5tk0Sb+xX5Gs3pvp686++xJXaGk99ppavo9TnKvffy7Zt4lfT7VAcCVUc+qUAxxBE3mdr1V9p80y9PnB+K3XnjcnRnhPToZvjx41fpeGbrgzdKNUT4wtF9M5Q6nOV/nuReoI2bQKuukoMQgsWAMOHG3uCiMgmhiDyPlur/ur14qq406fzf+pKMO05AawHojff7Lwa8aZNxp6GOXPsr+xsTVd7YnyhiN6bU9xNry82VvxMxowRQ5DpprOunFdNe6wRKYAhiHyTdHOReiMClTtvTL48/d4WX6iZ6up754uBwxtDv0RexhBEvkWrFXsZpM0hfWUWjbcocWOyNispPFz8e11d5+MB4KWXxJ6GHTvE3oYFC8SeBqld3qBUEb23FytUKnBIYap/f98LVUQqwBBEvkWrBXr0ABYvNn+c03ZdZ6/n5L33xLADGHvfpCHIwYONO8T36+f9glmliuhd6cnxxSnupmFq0CDPvjaRH2AIIt/ja7OjfI2zvRim72dODrBnj/nze/YYH7vjDuOxkpgY86/eZKuIPj3dvf82fGEYjogUxRBEvofTdm1zthfD9P288UYx8KxeDTQ3izfxBQvEnpS77gKGDJHv5ZgzRxxm8oVhJGtF9I5s3GqtJ6euTpyFOG+e9fMoEc69/f4RkRmGICI16eq0cGnYKz3d+Fi/fuLXDRvMj7Xs5dBqxZu/rxVTOxoo7PXkALZnISoRzqUAq9cbtxxxJ4YuIpsYgsi3+eIsGm9Qalq49H5qNMCxY+Jj0jpAUnGxVHDs60OQzv7bsNWTU1QErFjhvrZays/3TDBR6wxAIg9hCCLfxmm7IqXqUaT38/nnjesASUXFpqsNA8ZeDimA6fW+UxAsva4z/zbk2inNkGtuFr86el3OBDC5YTjAWHAeSLMdiXwMQxAFJjUNE7hjb7UxY8QQNHGieO4XX7TeK+LPBcGWvT+OXpczAczajDbptWfNUnZozBdnsRH5KIYgCkxqGiZwx95qUm3Qjh3ApEni36WeH8teDn+crSetRzVmjPheuPO6Zs8W/71Jw5iW8vOV7fH059BKpDCGICJf58q0cLmeLmvDMuXl4ldp0UTLG7KnZut5sndOq5Xf9sMd1yW9n7Nnm/eyuStI+mNoJXIThiAKHN4aJujqzd2VaeFyPV3WepSkWWGffCIOq0kLJnoqkEivlZGhnt45Z3ly2QcuMUHksCBvN4DIY/LygLlzxT9Sb8qaNcbH3LUjuem0aKW4Ek5mzxbDkzWffmp8D6y12R2z9aTXqq9X7pyOvGZOjvH6PDULUau1/RkQkUexJ4gChz8ME9i6WTvS09XVYRmlZ+tJ0/4B41R9T/XOmfY6uXpdzvaYSa/jqeJkLjFBZBNDEHmfp4ZePDlM4K6hN1s3a0cLYq29B56cCi+9Tl6esWBYKvRWUxGvKwX2nlz2gUtMENnEEETep6aZWo7yxgwdZ3q6pGEZ0xlLnmyztfokyaxZyodiTh0nIgsMQeT/5Hqa3D1M4I2hN2d6uuSGZTzZZrnXkvYwW73asb3AnKVUyGOYIvIbDEHkHZ68kcj1NLl7mEANM3S8NRXe2mtJe5i5azd4pUIe1+Eh8hsMQeQdvJG4lxoLYmNi3NtmpUKePxTYExEAhiDyFnffSDzV0+RIUbc3AklXe7o82WbptQYOdH7la29QQy8fETmEIYi8w903Ek/1NDlS1O3LM3SshTh/n8Gkxp4yIlIcQxD5Jw5ZOMYfZ+Y5QqngxTBFpGoMQeR97riRuLOnibODSOLLvXxEZBdDEHmf2m4kai/qdjTEeXL/MCIiL2AIIv+ndE+T2ofaHA1xgTpU5qsYSokUxxBE/k/pnia1zw5Se4gLVAylRIpjCCIKNLZCnCf3D3MUe0CIyE0Ygoi6wt9mB/livVMg94CwCJ/IrRiCiLpCbUXdlixDHIfKfIsvhlIiP8IQRBTIvLl/mC3sARExlBK5FUMQEfle3Q17QES+EkqJ/BRDEBHJ1914s96JPSBE5AEMQUS+xJd6ZLxZ78QekM78rQifyAcwBBH5Ek/OhGLdjbqovQifyAcxBBEFKrXU3bAHhIjchCGIyNu81SOjlrob9oAQkZswBBF5m7d6ZFh3Q0QBjiGIyNvU0iNDRORnGIKIvM0XemRYd0NEASjI2w1w1ssvv4z09HR069YNI0eOxJ49e7zdJCL1k+puGIKIKICoKgS9++67uPfee7F8+XIcOnQI48aNw/Tp01FcXOztphEpgz0yREQeoxEEQfB2Ixx1ySWXIDMzExs2bDA8dsEFF2DmzJlYu3Ztp+Obm5vR3Nxs+L6xsRGpqalo2LoVUSkpHmkzERERdU1jUxOiJ05EQ0MDoqKiFDuvanqCWlpacODAAUydOtXs8alTp+Kbb76R/Zm1a9ciOjra8Cc1NdUTTSVSP70eyMkxTt0nIvJDqglBer0e7e3tSExMNHs8MTERFRUVsj+zdOlSNDQ0GP6UlJR4oqlE6ietXM0QRER+THWzwzQajdn3giB0ekwSHh6O8PBwTzSLiIiIVEY1IUir1SI4OLhTr09VVVWn3iEicgH3EiOiAKOaEBQWFoaRI0fiiy++wKxZswyPf/HFF5gxY4YXW0bkJ9SylxgRkUJUE4IAYMmSJZg3bx5GjRqFrKwsvPrqqyguLsadd97p7aYRqR9XriaiAKOqEHTTTTehpqYGjz/+OHQ6HYYOHYpPPvkEffv29XbTiNTPF1auJiLyIFWFIAD485//jD//+c/ebgYRERGpnGqmyBO5hOvduIYrVxNRAGAIIv/G9W5cw73EiCgAMAQRERFRQFJdTRCRXVzvRjl6vTh1fvZsvmdE5HcYgsj/cL0b5UjDiePHMwQRkd9hCCL/w/VuiIjIAQxB5H+43k3XcDiRiAIEQxARmeNwIhEFCIYg8m9c78Z5HE4kogDBEET+TVrvhhzH4UQiChBcJ4iIiIgCEkMQEVnH4UQi8mMcDiMi6zicSER+jD1BRP6IG8cSEdnFEETkj7hxLBGRXQxBREREFJBYE0TkL7jSMxGRUxiCiPwFV3omInIKQxCRv+BKz0RETmEIIvIXXOmZiMgpLIwmIiKigMQQROSPuNIzEZFdHA4j8kdc6ZmIyC72BBEREVFAYggiIiKigMQQRERERAGJIYiIiIgCEkMQERERBSSGICIiIgpIDEFEREQUkBiCiIiIKCAxBBEREVFAYggiIiKigMQQRERERAGJIYiIiIgCEkMQERERBSSGICIiIgpIDEFEREQUkBiCiIiIKCAxBFFg0+uBnBzxKxERBRSGIApsej2wcSNDEBFRAGIIIiIiooAU4u0GEHmcXm/s+SksNP8KAFqt+IeIiPwaQxAFnrw8cQjM1Jo1xr/fcQewcKFn20RERB7HEESBZ/ZsYPx48e+FhWIAeuQRYPBg8TH2AhERBQSGIAo8csNdgwcbQ5C/0OvFXq/ZsxnsiIhksDCayF9x5hsRkU0MQRTYtFqxBog9JUREAYfDYRTYtFr/KoLmzDciIocxBBH5E858IyJyGEMQkT/hzDciIocxBBH5k0CZ+UZEpAAWRhMREVFAYggi8lec+UZEZBOHw4j8lb/NfCMiUhh7goiIiCggMQQRERFRQGIIIiIiooDEEEREREQBiSGIiIiIAhJDEBEREQUkhiAiIiIKSAxBREREFJBUE4KeeOIJXHrppejRowdiYmK83RwiIiJSOdWEoJaWFtxwww246667vN0UIiIi8gOq2TZj1apVAIDXX3/d4Z9pbm5Gc3Oz4fuGhgYAQOPp08DJk4q2j4iIiNyj8dQpAIAgCIqeVzUhyBVr1641hCdTqbfc4oXWEBERUVfU1NQgOjpasfP5dQhaunQplixZYvi+vr4effv2RXFxsaJvoq9rbGxEamoqSkpKEBUV5e3meAyvm9cdCHjdvO5A0NDQgLS0NMTFxSl6Xq+GoJUrV8r21JgqKCjAqFGjXDp/eHg4wsPDOz0eHR0dUP94JFFRUbzuAMLrDiy87sASqNcdFKRsKbNXQ9CiRYtw88032zymX79+nmkMERERBRSvhiCtVgutVuvNJhAREVGAUk1NUHFxMWpra1FcXIz29nYcPnwYAHDeeechIiLCoXOEh4fjsccekx0i82e8bl53IOB187oDAa9b2evWCErPN3OT2267DW+88Uanx3fs2IEJEyZ4vkFERESkaqoJQURERERKUs2K0URERERKYggiIiKigMQQRERERAGJIYiIiIgCkt+HoCeeeAKXXnopevTogZiYGId+RhAErFy5Er1790b37t0xYcIE/Pjjj+5tqMLq6uowb948REdHIzo6GvPmzUN9fb3Nn7ntttug0WjM/owZM8YzDXbRyy+/jPT0dHTr1g0jR47Enj17bB6/a9cujBw5Et26dUP//v3xyiuveKilynLmunfu3Nnpc9VoNCgsLPRgi7tu9+7duOaaa9C7d29oNBp8+OGHdn/GHz5vZ6/bHz7vtWvXYvTo0YiMjERCQgJmzpyJo0eP2v05tX/erly3P3zeGzZsQEZGhmEV7KysLHz66ac2f0apz9rvQ1BLSwtuuOEG3HXXXQ7/zNNPP42//e1vWL9+PQoKCpCUlIQpU6agqanJjS1V1h//+EccPnwY27dvx/bt23H48GHMmzfP7s9NmzYNOp3O8OeTTz7xQGtd8+677+Lee+/F8uXLcejQIYwbNw7Tp09HcXGx7PFFRUW48sorMW7cOBw6dAjLli3DX/7yF+Tl5Xm45V3j7HVLjh49avbZDhw40EMtVsapU6cwbNgwrF+/3qHj/eXzdva6JWr+vHft2oW7774b+/btwxdffIG2tjZMnToVp87tJC7HHz5vV65boubPu0+fPli3bh3279+P/fv3Y9KkSZgxY4bVzgdFP2shQLz22mtCdHS03eM6OjqEpKQkYd26dYbHzp49K0RHRwuvvPKKG1uonJ9++kkAIOzbt8/w2N69ewUAQmFhodWfmz9/vjBjxgwPtFAZF198sXDnnXeaPTZ48GDh4Ycflj3+wQcfFAYPHmz22MKFC4UxY8a4rY3u4Ox179ixQwAg1NXVeaB1ngFAyM/Pt3mMv3zephy5bn/8vKuqqgQAwq5du6we44+ftyPX7Y+ftyAIQmxsrJCbmyv7nJKftd/3BDmrqKgIFRUVmDp1quGx8PBwjB8/Ht98840XW+a4vXv3Ijo6GpdcconhsTFjxiA6OtruNezcuRMJCQk4//zzcccdd6CqqsrdzXVJS0sLDhw4YPY5AcDUqVOtXuPevXs7HX/FFVdg//79aG1tdVtbleTKdUtGjBiB5ORkTJ48GTt27HBnM32CP3zeXeFPn3dDQwMA2NxB3B8/b0euW+Ivn3d7ezveeecdnDp1CllZWbLHKPlZMwRZqKioAAAkJiaaPZ6YmGh4ztdVVFQgISGh0+MJCQk2r2H69Ol466238NVXX+G5555DQUEBJk2ahObmZnc21yV6vR7t7e1OfU4VFRWyx7e1tUGv17utrUpy5bqTk5Px6quvIi8vDx988AEGDRqEyZMnY/fu3Z5ostf4w+ftCn/7vAVBwJIlSzB27FgMHTrU6nH+9nk7et3+8nkfOXIEERERCA8Px5133on8/HwMGTJE9lglP2vV7B1mauXKlVi1apXNYwoKCjBq1CiXX0Oj0Zh9LwhCp8c8zdHrBjq3H7B/DTfddJPh70OHDsWoUaPQt29fbNu2Ddddd52LrXYvZz8nuePlHvd1zlz3oEGDMGjQIMP3WVlZKCkpwbPPPovLLrvMre30Nn/5vJ3hb5/3okWL8P333+O///2v3WP96fN29Lr95fMeNGgQDh8+jPr6euTl5WH+/PnYtWuX1SCk1GetyhC0aNEi3HzzzTaP6devn0vnTkpKAiAmzeTkZMPjVVVVnZKnpzl63d9//z0qKys7PVddXe3UNSQnJ6Nv37749ddfnW6ru2m1WgQHB3fq/bD1OSUlJckeHxISgvj4eLe1VUmuXLecMWPGYMuWLUo3z6f4w+etFLV+3osXL8ZHH32E3bt3o0+fPjaP9afP25nrlqPGzzssLAznnXceAGDUqFEoKCjAiy++iJycnE7HKvlZqzIEabVaaLVat5w7PT0dSUlJ+OKLLzBixAgAYh3Grl278NRTT7nlNR3l6HVnZWWhoaEB3333HS6++GIAwLfffouGhgZceumlDr9eTU0NSkpKzMKgrwgLC8PIkSPxxRdfYNasWYbHv/jiC8yYMUP2Z7KysvDxxx+bPfb5559j1KhRCA0NdWt7leLKdcs5dOiQT36uSvKHz1spavu8BUHA4sWLkZ+fj507dyI9Pd3uz/jD5+3KdctR2+ctRxAEq6UYin7WTpdSq8yJEyeEQ4cOCatWrRIiIiKEQ4cOCYcOHRKampoMxwwaNEj44IMPDN+vW7dOiI6OFj744APhyJEjwi233CIkJycLjY2N3rgEl0ybNk3IyMgQ9u7dK+zdu1e46KKLhKuvvtrsGNPrbmpqEu6//37hm2++EYqKioQdO3YIWVlZQkpKis9e9zvvvCOEhoYKmzZtEn766Sfh3nvvFXr27CkcP35cEARBePjhh4V58+YZjj927JjQo0cP4b777hN++uknYdOmTUJoaKjw/vvve+sSXOLsdT///PNCfn6+8Msvvwg//PCD8PDDDwsAhLy8PG9dgkuampoM//0CEP72t78Jhw4dEk6cOCEIgv9+3s5etz983nfddZcQHR0t7Ny5U9DpdIY/p0+fNhzjj5+3K9ftD5/30qVLhd27dwtFRUXC999/LyxbtkwICgoSPv/8c0EQ3PtZ+30Imj9/vgCg058dO3YYjgEgvPbaa4bvOzo6hMcee0xISkoSwsPDhcsuu0w4cuSI5xvfBTU1NcKcOXOEyMhIITIyUpgzZ06nKZSm13369Glh6tSpQq9evYTQ0FAhLS1NmD9/vlBcXOz5xjvh73//u9C3b18hLCxMyMzMNJtKOn/+fGH8+PFmx+/cuVMYMWKEEBYWJvTr10/YsGGDh1usDGeu+6mnnhIGDBggdOvWTYiNjRXGjh0rbNu2zQut7hppKrDln/nz5wuC4L+ft7PX7Q+ft9z1Wv5/2h8/b1eu2x8+7z/96U+G/5/16tVLmDx5siEACYJ7P2uNIJyrJiIiIiIKIJwiT0RERAGJIYiIiIgCEkMQERERBSSGICIiIgpIDEFEREQUkBiCiIiIKCAxBBEREVFAYggiIiKigMQQRERERAGJIYiIVOftt99Gt27dUFZWZngsOzsbGRkZaGho8GLLiEhNuG0GEamOIAgYPnw4xo0bh/Xr12PVqlXIzc3Fvn37kJKS4u3mEZFKhHi7AUREztJoNHjiiSdw/fXXo3fv3njxxRexZ88eQwCaNWsWdu7cicmTJ+P999/3cmuJyFexJ4iIVCszMxM//vgjPv/8c4wfP97w+I4dO3Dy5Em88cYbDEFEZBVrgohIlT777DMUFhaivb0diYmJZs9NnDgRkZGRXmoZEakFQxARqc7Bgwdxww03ICcnB1dccQVWrFjh7SYRkQqxJoiIVOX48eO46qqr8PDDD2PevHkYMmQIRo8ejQMHDmDkyJHebh4RqQh7gohINWprazF9+nRce+21WLZsGQBg5MiRuOaaa7B8+XIvt46I1IY9QUSkGnFxcfj55587Pb5161YvtIaI1I6zw4jI71xxxRU4ePAgTp06hbi4OOTn52P06NHebhYR+RiGICIiIgpIrAkiIiKigMQQRERERAGJIYiIiIgCEkMQERERBSSGICIiIgpIDEFEREQUkBiCiIiIKCAxBBEREVFAYggiIiKigMQQRERERAGJIYiIiIgC0v8HzBK2E3eSgcsAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make sure you run the data-generating code cell first\n",
    "\n",
    "# Multinomial (in this case binomial) density estimation\n",
    "p_apple_est = sum(y.==true) / length(y)\n",
    "π_hat = [p_apple_est; 1-p_apple_est]\n",
    "\n",
    "# Estimate class-conditional multivariate Gaussian densities\n",
    "d1 = fit_mle(FullNormal, X_apples')  # MLE density estimation d1 = N(μ₁, Σ₁)\n",
    "d2 = fit_mle(FullNormal, X_peaches') # MLE density estimation d2 = N(μ₂, Σ₂)\n",
    "Σ = π_hat[1]*cov(d1) + π_hat[2]*cov(d2) # Combine Σ₁ and Σ₂ into Σ\n",
    "conditionals = [MvNormal(mean(d1), Σ); MvNormal(mean(d2), Σ)] # p(x|C)\n",
    "\n",
    "# Calculate posterior class probability of x∙ (prediction)\n",
    "function predict_class(k, X) # calculate p(Ck|X)\n",
    "    norm = π_hat[1]*pdf(conditionals[1],X) + π_hat[2]*pdf(conditionals[2],X)\n",
    "    return π_hat[k]*pdf(conditionals[k], X) ./ norm\n",
    "end\n",
    "println(\"p(apple|x=x∙) = $(predict_class(1,x_test))\")\n",
    "\n",
    "# Discrimination boundary of the posterior (p(apple|x;D) = p(peach|x;D) = 0.5)\n",
    "β(k) = inv(Σ)*mean(conditionals[k])\n",
    "γ(k) = -0.5 * mean(conditionals[k])' * inv(Σ) * mean(conditionals[k]) + log(π_hat[k])\n",
    "function discriminant_x2(x1)\n",
    "    # Solve discriminant equation for x2\n",
    "    β12 = β(1) .- β(2)\n",
    "    γ12 = (γ(1) .- γ(2))[1,1]\n",
    "    return -1*(β12[1]*x1 .+ γ12) ./ β12[2]\n",
    "end\n",
    "\n",
    "plot_fruit_dataset() # Plot dataset\n",
    "x1 = range(-1,length=10,stop=3)\n",
    "plot(x1, discriminant_x2(x1), \"k-\") # Plot discrimination boundary\n",
    "fill_between(x1, -1, discriminant_x2(x1), color=\"r\", alpha=0.2)\n",
    "fill_between(x1, discriminant_x2(x1), 4, color=\"b\", alpha=0.2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Recap Generative Classification\n",
    "\n",
    "- Gaussian-Categorical Model specification:  \n",
    "\n",
    "$$p(x,\\mathcal{C}_k|\\,\\theta) = \\pi_k \\cdot \\mathcal{N}(x|\\mu_k,\\Sigma)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- If the class-conditional distributions are Gaussian with equal covariance matrices across classes ($\\Sigma_k = \\Sigma$), then\n",
    "    the discriminant functions are hyperplanes in feature space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- ML estimation for $\\{\\pi_k,\\mu_k,\\Sigma\\}$ in the GCM model breaks down to simple density estimation for Gaussian and multinomial/categorical distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Posterior class probability is a softmax function\n",
    "$$ p(\\mathcal{C}_k|x,\\theta ) \\propto \\exp\\{\\beta_k^T x + \\gamma_k\\}$$\n",
    "where $\\beta _k= \\Sigma^{-1} \\mu_k$ and $\\gamma_k=- \\frac{1}{2} \\mu_k^T \\Sigma^{-1} \\mu_k  + \\log \\pi_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<!--\n",
       "This HTML file contains custom styles and some javascript.\n",
       "Include it a Jupyter notebook for improved rendering.\n",
       "-->\n",
       "\n",
       "<!-- Fonts -->\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<!-- Custom style -->\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.container {\n",
       "    min-width: 960px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "    margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running {\n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256);\n",
       "    border-radius: 0px;\n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Alegreya Sans' sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 125%;\n",
       "    font-weight: 400;\n",
       "    width:800px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: 400;\n",
       "    font-size: 45pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    margin-top:16px;\n",
       "    font-size: 22pt;\n",
       "    font-weight: 600;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: rgb(102,102,0);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-size: 14pt;\n",
       "    text-align: center;\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Nixie One', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 16pt;\n",
       "    color: rgb(163,0,0);\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.8em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"PT Mono\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       ".boxed { /* draw a border around a piece of text */\n",
       "  border: 1px solid blue ;\n",
       "}\n",
       "\n",
       "h4#CODE-EXAMPLE,\n",
       "h4#END-OF-CODE-EXAMPLE {\n",
       "    margin: 10px 0;\n",
       "    padding: 10px;\n",
       "    background-color: #d0f9ca !important;\n",
       "    border-top: #849f81 1px solid;\n",
       "    border-bottom: #849f81 1px solid;\n",
       "}\n",
       "\n",
       ".emphasis {\n",
       "    color: red;\n",
       "}\n",
       "\n",
       ".exercise {\n",
       "    color: green;\n",
       "}\n",
       "\n",
       ".proof {\n",
       "    color: blue;\n",
       "}\n",
       "\n",
       "code {\n",
       "  padding: 2px 4px !important;\n",
       "  font-size: 90% !important;\n",
       "  color: #222 !important;\n",
       "  background-color: #efefef !important;\n",
       "  border-radius: 2px !important;\n",
       "}\n",
       "\n",
       "/* This removes the actual style cells from the notebooks, but no in print mode\n",
       "   as they will be removed through some other method */\n",
       "@media not print {\n",
       "  .cell:nth-last-child(-n+2) {\n",
       "    display: none;\n",
       "  }\n",
       "}\n",
       "\n",
       "footer.hidden-print {\n",
       "    display: none !important;\n",
       "}\n",
       "    \n",
       "</style>\n",
       "\n",
       "<!-- MathJax styling -->\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"],\n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "open(\"../../styles/aipstyle.html\") do f\n",
    "    display(\"text/html\", read(f,String))\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Julia 1.8.2",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}