13_Logistic Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the purposes of classification, we have so far explored a class of model called the *generative model*, so called because it is possible to generate synthetic based on this model.  Markov models (and HMMs), graphical models, and naive Bayes all fell into this category.  This process was based on the notion that we could specify the joint probability distribution of the class and the feature values:\n",
    "$$P(\\mathbf{x},y),$$\n",
    "where $\\mathbf{x}$ is a vector of features, and $y$ is a class label.\n",
    "\n",
    "It was then easy to use Bayes' rule to infer the class:\n",
    "$$\n",
    "P(y|\\mathbf{x}) \\propto P(\\mathbf{x}|y)P(y).\n",
    "$$\n",
    "However, this required the specification of a statistical model for each class (what is the distribution of $\\mathbf{x}$ given y).  Consider the lobster example for naive Bayes, in which we assumed that the features $\\mathbf{x}$ were normally distributed (independent from one another), with some mean and covariance that depended on which class they were a member of. \n",
    "$$\n",
    "P(\\mathbf{x}|y) = \\begin{cases} \\mathcal{N}(\\mathbf{x}|\\boldsymbol{\\mu}_0,\\Sigma_0) & \\mathrm{if}\\;y=0, \\\\\n",
    "                                \\mathcal{N}(\\mathbf{x}|\\boldsymbol{\\mu}_1,\\Sigma_1) & \\mathrm{if}\\;y=1 \\end{cases}\n",
    "$$\n",
    "where $\\boldsymbol{\\mu}_1$ and $\\Sigma_1$ are the mean and diagonal covariance matrix of the *survived* class (and conversely for $\\boldsymbol{\\mu}_0$, etc.).  Similarly we had a prior\n",
    "$$\n",
    "P(y=1) = \\pi_1 = 1-P(y=0) = 1-\\pi_0\n",
    "$$\n",
    "which is the probability of a lobster's survival without knowing its size.\n",
    "\n",
    "This was straightforward to code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'P(survival|lobster size)')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGL1JREFUeJzt3Xu0JWV95vHvYwvCIDZm6Kwg0DYqiYMBdNIy4yVeEA2KQjQygjKi6DBEiYrR2GoWNhhXEBPRZEGEKF6SIEFRVwNtwGvi3W4NgoBoDyI0qLQXGhBEW37zx66WvY+nz64+3XXqXL6ftfY6VbXfqv3sBX1+p96qet9UFZIkbXafvgNIkmYXC4MkaYSFQZI0wsIgSRphYZAkjbAwSJJGWBgkSSMsDJKkERYGSdKI+/YdYDp23333WrZsWd8xJGnO+NrXvvajqlrSpu2cLAzLli1j7dq1fceQpDkjyffatrUrSZI0wsIgSRphYZAkjbAwSJJGWBgkSSMsDJKkERYGSdIIC4MkaYSFQZI0Yk4++Sx14oz9YeMNfaeAxUvhpCv7TqEFzMIgbbbxBli5se8UsHJx3wm0wNmVJEkaYWGQJI2wMEiSRlgYJEkjLAySpBEWBknSCAuDJGmEzzEsYI877dPcdOtdfcdgz9125gsrDu47hqSGhWEBu+nWu7j+tMP6jsGyFZf0HUHSELuSJEkjLAySpBEWBknSCAuDJGlEp4UhyaFJrk2yLsmKSd5/UZINSS5vXi/tMo8kabzO7kpKsgg4E3gqsB5Yk2RVVV09oem/VtWJXeWQJG2dLs8YDgLWVdV1VfUL4HzgiA4/T5K0HXRZGPYEbhxaX99sm+hPklyR5MNJ9u4wjySphS4LQybZVhPWLwKWVdUBwCeB92/xYMnxSdYmWbthw4btGFOSNKzLwrAeGD4D2Au4ebhBVf24qu5uVv8R+IMtHayqzqmq5VW1fMmSJds9rCRpoMvCsAbYN8k+SXYEjgJWDTdIssfQ6uHANR3mkSS10NldSVW1KcmJwKXAIuDcqroqyanA2qpaBbwiyeHAJuAnwIu6yiNJaqfTQfSqajWwesK2k4eWXw+8vssMkqSt45PPkqQRFgZJ0ggLgyRphIVBkjTCwiBJGmFhkCSNsDBIkkZYGCRJIywMkqQRFgZJ0ggLgyRphIVBkjTCwiBJGmFhkCSNsDBIkkZYGCRJIywMkqQRFgZJ0ggLgyRpxNg5n5PsBDwT+EPgQcBdwDeBS6rqqm7jSZJm2pSFIclK4FnAZ4GvALcAOwG/C5zWFI0/r6oruo0pSZop484Y1lTVyi289/Ykvw0s3b6RJEl9mrIwVNUlw+tJdqmqnw29fwuDswhJ0jzR6uJzkscmuRq4plk/MMlZnSaTJPWi7V1JZwB/BPwYoKq+ATyhq1CSpP60vl21qm6csOlX2zmLJGkWGHu7auPGJI8FKsmOwCtoupUkSfNL2zOGE4CXA3sC64FHNuuSpHmm7RnDzlX1guENSX6ngzxSr5atuGR8o45dv1PfCbTQtS0M303yIeC4qrqr2bYa+O9T7ZTkUOCdwCLg3VV12hbaPRf4EPDoqlrbMpO03V1/2mF9R4CVfQfQQte2K+lK4HPA55M8tNmWqXZIsgg4E3g6sB9wdJL9Jmm3K4NrFl9pG1qS1J22haGq6iwGv8AvSvIsoMbscxCwrqquq6pfAOcDR0zS7s3A6cDPW2aRJHWobWEIQFV9AXgK8Frg4WP22RMYvsV1fbPt3oMmjwL2rqqLxwZIjk+yNsnaDRs2tIwtSdpabQvDMzYvVNX3gYOBQ8fsM1lX06/PMpLch8GDc3/eJkBVnVNVy6tq+ZIlS9rsIkmahnGjqx5TVf/M4PrAZE3+Y4rd1wN7D63vBdw8tL4r8PvAZ5tj/w6wKsnhXoCWpP6Muytpl+bnrtM49hpg3yT7ADcBRwHP3/xmVW0Edt+8nuSzwGssCpLUr3Gjq57d/Dxlaw9cVZuSnAhcyuB21XOr6qokpwJrq2rVdAJLkrrV6jmGJKcDf8Vg9rZ/Aw4EXtV0M21RVa1m8LzD8LaTt9D2SW2ySJK61fbi89Oq6jYGU3yuZzCD22s7SyVJ6k3bwrBD8/MZwAer6icd5ZEk9aztkBgXJfkWg66klyVZgg+kSdK81OqMoapWAI8BllfVL4E7mfwpZknSHNf2jIGq+unQ8s+An03RXJI0R7WewU2StDCMLQwZ2HtcO0nS/DC2MFRVAR+bgSySpFmgbVfSl5M8utMkkqRZoe3F5ycDJyS5nsFF5zA4mTigq2CSenbG/rDxhr5TwOKlcNKVfadYUNoWhqd3mkLS7LPxBli5se8UsHJx3wkWnLbPMXyPwRDaBzfLd7bdV5I0t7T65Z7kTcDrgNc3m3YAphxAT5I0N7X9q//ZwOE0D7VV1c1Mb44GSdIs17Yw/KK5bbUAkuwypr0kaY5qWxguSHI2sFuS/wN8Enh3d7EkSX1pdVdSVf1NkqcCtwG/B5xcVZ/oNJkkqRdtZ3B7a1W9DvjEJNskSfNI266kp06yzWcbJGkemvKMIcmfAi8DHpLkiqG3dgW+0GUwSVI/xnUlnQd8HPhrYMXQ9tud3lOS5qcpu5KqamNVXQ/8JfCD5qnnfYBjkuw2A/kkSTOs7TWGC4FfJXkY8B4GxeG8zlJJknrTtjDcU1WbgOcA76iqk4A9uoslSepL28LwyyRHAy8ELm627dBNJElSn9oWhhcDjwHeUlXfTbIPDqInSfNS22G3rwZeA1yVZH/gpqo6rdNkkqRetH3y+TDgXcD/YzB72z5J/m9VfbzLcJKkmde2K+lvgSdX1ZOq6okMpvo8Y9xOSQ5Ncm2SdUlWTPL+CUmuTHJ5ks8n2W/r4kuStre2heGWqlo3tH4dcMtUOyRZBJzJYOiM/YCjJ/nFf15V7V9VjwROB97eMo8kqSPjhsR4TrN4VZLVwAUM5mQ4Elgz5tgHAeuq6rrmWOcDRwBXb25QVbcNtd+lObYkqUfjrjE8a2j5h8ATm+UNwAPH7LsncOPQ+nrgf0xslOTlwKuBHYGDxxxTktSxKQtDVb14G46dyQ45yWecCZyZ5PkMht44dtKDJccDxwMsXbp0G2JJkqbS6hpDkr2SfDTJLUl+mOTCJHuN2W09sPfQ+l7AzVO0Px/44y29WVXnVNXyqlq+ZMmSNrElSdPQ9uLze4FVwIMYdBFd1Gybyhpg3yT7JNkROKo5xq8l2Xdo9TDgOy3zSJI60uo5BmBJVQ0XgvcledVUO1TVpiQnApcCi4Bzq+qqJKcCa6tqFXBikkOAXwI/ZQvdSJKkmdO2MPwoyTHAB5v1o4Efj9upqlYDqydsO3lo+ZUtP1+SNEPadiUdB/wv4AfA94HnNtskSfNMqzOGqroBOLzjLJKkWWDcA25/zxQPnVXVK7Z7IklSr8adMaydkRSSpFlj3ANu75+pIJKk2WFcV9ILWx7n8qq6YjvkkST1bFxX0j4tj3P9NuaQJM0S47qSTpmpIJKk2aHtsNuTqqqPbN84kqS+bc2w2xMVYGHQNttzt51ZtuKSvmNw/U59J5Bmhy6H3ZZa+cKKWTINx8q+A0izQ9uxkkhyGPAI4Nd/V1XVqV2EkiT1p+18DO8Cngf8GYMJeI4EHtxhLklST9oOovfYqnoh8NPmTqXHMDoJjyRpnmhbGO5qft6Z5EEM5k9o+4yDJGkOaXuN4eIkuwFvA77O4I6kf+wslbSAra/d2Wvl4r5jwGLnVl+o2g67/eZm8cIkFwM7VdXG7mJJC9fj7/47rj/tsL5jaAFre/H5G0nekOShVXW3RUGS5q+21xgOBzYBFyRZk+Q1STzPlKR5qFVhqKrvVdXpVfUHwPOBA4DvdppMktSLrXnAbRmDeZ+fB/wK+ItuIkmS+tSqMCT5CrAD8CHgyKq6rtNUkqTetD1jOLaqvtVpEknSrDBu2O1jquqfgWckecbE96vq7Z0lkyT1YtwZwy7Nz127DiJJmh3GDbt9drN4VlVtmIE8kqSetX2O4YtJLkvykiQP7DSRJKlXbZ9j2Bf4SwbzMXwtycVJjuk0mSSpF23PGKiqr1bVq4GDgJ8A7+8slSSpN23HSnpAkmOTfBz4IvB9BgVi3H6HJrk2ybokKyZ5/9VJrk5yRZJPJXHyH0nqWdvnGL4BfAw4taq+1GaHJIuAM4GnAuuBNUlWVdXVQ83+E1heVXcm+VPgdAZPVkuSejL2jKH5Bf/RqjqpbVFoHASsq6rrquoXwPnAEcMNquozVXVns/plYK+tOL4kqQNjC0NV/Qo4cBrH3hO4cWh9fbNtS14CfHwanyNJ2o7adiVdnmQVg7GSfrZ5Y1V9ZIp9Msm2mrTh4A6n5cATt3iw5HjgeIClSx3xW5K60rYw/BbwY+DgoW0FTFUY1gN7D63vBdw8sVGSQ4A3Ak+sqru3dLCqOgc4B2D58uWTFhhJ0rZrO7Xni6dx7DXAvkn2AW4CjmIwl8OvJXkUcDZwaFXdMo3PkCRtZ22H3X4vk3QDVdVxW9qnqjYlORG4FFgEnFtVVyU5FVhbVauAtwH3Bz6UBOCGqjp867+GJGl7aduVdPHQ8k7As5mkW2iiqloNrJ6w7eSh5UNafr4kaYa07Uq6cHg9yQeBT3aSSJLUq9ZDYkywL+CtQZI0D7W9xnA7o9cYfgC8rpNEkqRete1KcqIeSVog2g6i97gkuzTLxyR5uwPeSdL81PYawz8AdyY5EPgL4HvABzpLJUnqTdvCsKmqisEgeO+sqnfiPNCSNC+1fY7h9iSvB44BntCMuLpDd7EkSX1pe8bwPOBu4CVV9QMGo6S+rbNUkqTeTHnGkCQ18APg7Zu3V9UNNNcYNrfpNqYkaaaMO2P4TJI/SzLyMFuSHZMcnOT9wLHdxZMkzbRx1xgOBY4DPpjkIcBPGYyVtAi4DDijqi7vNqIkaSZNWRiq6ufAWcBZSXYAdgfuqqpbZyKcJGnmjbvGsBNwAvAw4AoGQ2dvmolgkqR+jLvG8H4GU25eCTwD+NvOE0mSejXuGsN+VbU/QJL3AF/tPpIkqU/jzhh+uXnBLiRJWhjGnTEcmOS2ZjnAzs16gKqqB3SaTpI048bdlbRopoJIkmaH6c7gJkmap9oOoid154z9YeMNfaeAxc5WK4GFQbPBxhtg5ca+U0hq2JUkSRphYZAkjbAwSJJGWBgkSSMsDJKkEd6VJGlSjzvt09x06119x+BLOy1hj5WL+w2xeCmcdGW/GWaQhUHSpG669S6uP+2wvmOwbAX95+i7MM2wTruSkhya5Nok65KsmOT9JyT5epJNSZ7bZRZJUjudFYYki4AzgacD+wFHJ9lvQrMbgBcB53WVQ5K0dbrsSjoIWFdV1wEkOR84Arh6c4Oqur55754Oc0iStkKXXUl7AjcOra9vtk1LkuOTrE2ydsOGDdscTpI0uS4LQybZVtM9WFWdU1XLq2r5kiVLtiGWJGkqXRaG9cDeQ+t7ATd3+HmSpO2gy2sMa4B9k+wD3AQcBTy/w8+bU2bDPeJ77rZzr5+vye25284sW3FJ3zH8/2MB66wwVNWmJCcClwKLgHOr6qokpwJrq2pVkkcDHwUeCDwrySlV9YiuMs0ms+Uecc0+X1hxcN8RtMB1+oBbVa0GVk/YdvLQ8hoGXUySpFnCsZIkSSMsDJKkERYGSdIIB9FbyM7YfzDfct8WL+07gaQhFoaFbOMNsHJj3ykkzTJ2JUmSRlgYJEkjLAySpBEWBknSCAuDJGmEhUGSNMLCIEkaYWGQJI2wMEiSRlgYJEkjLAySpBEWBknSCAuDJGmEhUGSNMLCIEkaYWGQJI2wMEiSRlgYJEkjLAySpBEWBknSiPv2HUCSZr3FS2Hl4r5TDHKcdGXnH2NhkKRxZuCXcSszVJzsSpIkjei0MCQ5NMm1SdYlWTHJ+/dL8q/N+19JsqzLPJKk8TorDEkWAWcCTwf2A45Ost+EZi8BflpVDwPOAN7aVR5JUjtdnjEcBKyrquuq6hfA+cARE9ocAby/Wf4w8JQk6TCTJGmMLgvDnsCNQ+vrm22TtqmqTcBG4L92mEmSNEaXdyVN9pd/TaPNoGFyPHB8s3pHkmunmWt34EfT3He7yrZ3nG37dzllVpygzZr/JtuB36UD2/hvZdZ8j+1gd07JdL/Lg9s27LIwrAf2HlrfC7h5C23WJ7kvsBj4yWQHq6pzgHO2NVSStVW1fFuPMxvMl+8yX74H+F1mo/nyPWDmvkuXXUlrgH2T7JNkR+AoYNWENquAY5vl5wKfrqpJzxgkSTOjszOGqtqU5ETgUmARcG5VXZXkVGBtVa0C3gP8U5J1DM4UjuoqjySpnU6ffK6q1cDqCdtOHlr+OXBklxkmsc3dUbPIfPku8+V7gN9lNpov3wNm6LvEnhtJ0jCHxJAkjViQhSHJm5NckeTyJJcleVDfmaYjyduSfKv5Lh9NslvfmaYryZFJrkpyT5I5dwfJuOFf5pIk5ya5Jck3+86yLZLsneQzSa5p/t96Zd+ZpivJTkm+muQbzXc5pdPPW4hdSUkeUFW3NcuvAParqhN6jrXVkjyNwZ1cm5LBnd5V9bqeY01Lkv8G3AOcDbymqtb2HKm1ZviXbwNPZXAL9hrg6Kq6utdg05TkCcAdwAeq6vf7zjNdSfYA9qiqryfZFfga8Mdz8b9LMyLELlV1R5IdgM8Dr6yqL3fxeQvyjGFzUWjswhYeqpvtquqy5olxgC8zeFZkTqqqa6pqug8t9q3N8C9zRlX9B1t4nmguqarvV9XXm+XbgWv4zdEX5oQauKNZ3aF5dfZ7a0EWBoAkb0lyI/AC4ORx7eeA44CP9x1igWoz/It61Izc/CjgK/0mmb4ki5JcDtwCfKKqOvsu87YwJPlkkm9O8joCoKreWFV7A/8CnNhv2i0b9z2aNm8ENjH4LrNWm+8yR7Ue2kUzL8n9gQuBV03oLZhTqupXVfVIBj0DByXprJtv3s7gVlWHtGx6HnAJ8KYO40zbuO+R5FjgmcBTZvtT41vx32SuaTP8i3rQ9MdfCPxLVX2k7zzbQ1XdmuSzwKFAJzcIzNszhqkk2Xdo9XDgW31l2RZJDgVeBxxeVXf2nWcBazP8i2ZYc8H2PcA1VfX2vvNsiyRLNt91mGRn4BA6/L21UO9KuhD4PQZ3wXwPOKGqbuo31dZrhhK5H/DjZtOX5+LdVQBJng38PbAEuBW4vKr+qN9U7SV5BvAO7h3+5S09R5q2JB8EnsRgVNIfAm+qqvf0Gmoakjwe+BxwJYN/6wBvaEZkmFOSHMBg7ppFDP6gv6CqTu3s8xZiYZAkbdmC7EqSJG2ZhUGSNMLCIEkaYWGQJI2wMEiSRlgYJEkjLAyaE5LcMb7Vr9u+L8lzW7bdLcnLpp9symOfkOSF2+E4j0ry7u1wnP2TvG9bj6P5z8KghW43YKsKQwbG/tupqndV1Qemnexeb2Dw8N82qaorgb2SLN32SJrPLAyas5I8OMmnmomKPjXhF94hST6X5NtJntm0f0Qz2cnlzT77AqcBD222va1p99oka5o2pzTbljUTvpwFfJ3RsZFIclqSq5t9/qbZtjLJa5I8qDn+5tevmuxLklzYfNaaJI+b5DvuChxQVd9o1u+f5L1Jrmw+60+a7XckeWuSrzWDFR6U5LNJrkty+NAhL2IwZIe0ZVXly9esfwF3TLLtIuDYZvk44GPN8vuAf2Pwh8++DAa524nBX90vaNrsCOwMLAO+OXTMpzGYcD3N/hcDT2ja3QP8z0ly/BZwLfeOJLBb83Mlg0mHhtu+nMFwBjAYwPHxzfJSBmP6TDz2k4ELh9bfCrxjaP2Bzc8Cnt4sfxS4jMGY/QcyGF5kc/vHARf1/d/T1+x+zdvRVbUgPAZ4TrP8T8DpQ+9dUFX3AN9Jch3wcOBLwBuT7AV8pKq+MxhnbcTTmtd/Nuv3Z1BcbgC+V5PPmHUb8HPg3UkuYVBMfkNzRvBS4A+bTYcA+w1leECSXWswqcxmewAbhtYPYegv/qr6abP4CwbFEAZjA91dVb9MciWDorbZLcCcnMpWM8fCoPmktrAMg0mwzkvyFeAw4NIkLwWum9AuwF9X1dkjGwcTvfxs0g8dTK16EPAUBr+0TwQOnrD/HgxG+jy87p2J6z7AY6rqrim+010MznaG8002wNkvq2rz9nuAu5ts9yQZ/ne+U3NMaYu8xqC57Ivc+9fzCxjMg7vZkUnuk+ShwEOAa5M8BLiuqv6OwbDYBwC3A7sO7XcpcFwzuQtJ9kzy21OFaNoursGona8CHjnh/R2AC4DXVdW3h966jKFJopKM7Ne4BnjYFPs8cKpsk/hdOhrDX/OHhUFzxX9Jsn7o9WrgFcCLk1wB/G/glUPtrwX+ncF0pydU1c+B5wHfzGB6xIczmOz+x8AXmpnk3lZVlzHo+/9S0w3zYUYLx2R2BS5ucvw7cNKE9x8LPBo4ZegC9IOa/Mubi8hXA78xZHpVfQtY3FyEBvgr4IFN3m8wuAaxNZ7MYGIqaYscdlua5ZKcBNxeVdv0LEOS+zEoXI+vqk3bJZzmJc8YpNnvH2iuGWyjpcAKi4LG8YxBkjTCMwZJ0ggLgyRphIVBkjTCwiBJGmFhkCSN+P/3p96ynIBqGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from __future__ import division    # Eliminates annoying default integer division behavior\n",
    "import numpy as np \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_csv('datasets/lobster_survive.dat',header=0,sep=r\"\\s{2,}\",engine='python')\n",
    "X = data['Len'].values\n",
    "\n",
    "# z-normalization\n",
    "X = X - X.mean()\n",
    "X = X/X.std()\n",
    "\n",
    "Y = data['Survive'].values\n",
    "\n",
    "m = len(Y)\n",
    "N = 2\n",
    "n = 1\n",
    "\n",
    "plt.hist(X[y==0],np.linspace(-3,3,10),density=True,histtype='step')\n",
    "plt.hist(X[y==1],np.linspace(-3,3,10),density=True,histtype='step')\n",
    "plt.xlabel('Lobster size (cm)')\n",
    "plt.ylabel('P(survival|lobster size)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we fit our model by finding maximum likelihood estimates of the joint distribution\n",
    "$$\n",
    "P(X,Y|\\Theta) = \\prod_{i=1}^m P(\\mathbf{x}_i|y_i) P(y_i),\n",
    "$$\n",
    "where\n",
    "$$\n",
    "X = \\begin{bmatrix} \\mathbf{x}_1 \\\\\n",
    "                    \\mathbf{x}_2 \\\\\n",
    "                    \\vdots \\\\\n",
    "                    \\mathbf{x}_m \\end{bmatrix},\n",
    "$$\n",
    "and \n",
    "$$\n",
    "Y = \\begin{bmatrix} y_1 \\\\\n",
    "                    y_2 \\\\\n",
    "                    \\vdots \\\\\n",
    "                    y_m \\end{bmatrix}\n",
    "$$\n",
    "are the training data and training labels, and $\\Theta$ is the vector of parameters ($\\Theta = [\\mu_0,\\mu_1,\\Sigma_0,\\Sigma_1,\\pi_0,\\pi_1]$).\n",
    "\n",
    "The maximum likelihood estimators for which for the means and covariances of each Gaussian turned out to be the sample mean and sample variance of the training data for class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'P(survival|lobster size)')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4U3X7x/H33TTdewEyZO+l4gJRcCJTRRBoEdwDxa3o86iIP/dAUUAQB8oSxMFyghMVAWXIEgTZ0NLd0t3790eCT4tAQ5v0NOn3dV25SNKTcz65aHrnnPM991dUFcMwDMM4ws/qAIZhGEbNYgqDYRiGUY4pDIZhGEY5pjAYhmEY5ZjCYBiGYZRjCoNhGIZRjikMhmEYRjmmMBiGYRjlmMJgGIZhlONvdYDKiIuL08aNG1sdwzAMw2usXr36kKrGu7KsVxaGxo0bs2rVKqtjGIZheA0R2enqsuZQkmEYhlGOKQyGYRhGOaYwGIZhGOWYwmAYhmGUYwqDYRiGUY4pDIZhGEY5pjAYhmEY5ZjCYBiGYZRjCoNhGIZRjlde+WwYHjG+A2TusjoFRDaCe9ZbncKoxUxhMIwjMnfB2EyrU8DYSKsTGLWcOZRkGIZhlGMKg2EYhlGOKQyGYRhGOaYwGIZhGOWYwmAYhmGUYwqDYRiGUY4pDIZhGEY55jqGWqzbs8vYm5FndQzqRwWzfMyFVscwDMPJFIZabG9GHn8/28fqGDQes9jqCIZhlGEOJRmGYRjlmMJgGIZhlGMKg2EYhlGOKQyGYRhGOR4tDCLSS0S2iMg2ERlzjJ+PFJEUEVnjvN3oyTyGYRhGxTw2KklEbMBE4BJgD7BSRBao6sajFv1AVe/wVA7DMAzj5Hhyj+EsYJuqblfVQmAOMMCD2zMMwzDcwJOFoT6wu8zjPc7njjZQRNaJyIci0tCDeQzDMAwXeLIwyDGe06MeLwQaq2pH4Gtg+nFXJnKziKwSkVUpKSlujGkYhmGU5cnCsAcouwfQANhXdgFVTVXVAufDN4EzjrcyVZ2qql1UtUt8fLzbwxqGYRgOniwMK4EWItJERAKAIcCCsguISL0yD/sDmzyYxzAMw3CBx0YlqWqxiNwBfAHYgLdVdYOIjANWqeoCYLSI9AeKgTRgpKfyGIZhGK7xaBM9VV0CLDnqucfK3H8YeNiTGQzDMIyTY658NgzDMMoxhcEwDMMoxxQGwzAMoxxTGAzDMIxyTGEwDMMwyjGFwTAMwyjHFAbDMAyjHFMYDMMwjHJMYTAMwzDKMYXB8Bnp+enkFOZYHcMwvJ5HW2IYRnU5lHeIgQsGkl+czxXNryCpTRINI8z0HoZRGWaPwfB6qsp/l/+X3KJcLmhwAXP/nEufj/swetloVh5YierR04D4Ni0poeCvvyg+dMjqKIaXMnsMhtebvXk2y/cu55GzH2Fo66EkH07mgy0fMHfLXL7Z/Q2tY1ozvO1wejXuRYAtwOq4bqWqFB88SN66deSvX0/euvXk//EHpbm52Bs0oOmCT/ELCbE6puFlTGEwvNpfGX/x8uqX6V6/O0NaDQEgISSBO0+7k5s63MSi7YuYsXEG//nxP4xfPZ5rWl3D4FaDiQmKsTh55ZRkZZG3fr2jCKz/g/x16yg+MqOh3U5Q69ZEDhiAf506pIwfT8rEidR54AFrQxtexxQGw2sVlhTy0PcPEWoPZVy3cYiUn002yD+Iq1tezcAWA/l538+8t+k9Jq6ZyJvr3qRvs74ktUmiRXQLi9JXrLSggILNm8lbt5689evIX7eewr///ufnAU2aENr1XII6dCS4YwcCW7fGL+B/e0RFe/aQ9u50Ivv0IahtWwvegeGtTGEwvNaE3yawJX0Lr1/4OnHBccddTkToWr8rXet3ZXvGdmZsmsHCvxby0daPOLfeuSS1TeK8+udZesJNS0sp3LHDcShoVSR5Vw8if8sWKCoCwD8+nqBOHYm84gqCO3YgqH17bBERJ1xnwv33kf3NN+x/9DEaz/0Asdmq460YPsAUBsMr/bL/F6ZvnM7gloO5oOEFLr+uaVRTHjv3MUafNpoPt37I7E2zGbV0FI0jGpMUHka/osOE2KvnmLyWlJD69tvkLv+J/PXrKc3NBcDPP5igLmHEjhxJUMcOBHfsiL1OnZNevy0ykrr/eYS999xL+owZxIwY4e63YPgoUxgMr5NZkMl/fvwPjSMac/+Z91dqHVFBUdzY4UZGtB3BFzu/4P2N7/N/cTFM+PASBrUcxJDWQ6gbWtfNycvLXLiQlJdeJrBNGyIH9P/nkFDAe6cjT7zjlm2E9+pF2CefkvzqBMIvvhh7/fpuWa/h28xwVcOrqCpP/PwEaflpPHf+cwT7B1dpfXabnb5N+zKnzxym7zvIWXXP4p0N73D5/Mt58PsH+ePQH25KXl5pQQEpEyYQ1L49TeZ/SN3HHiPqyisIbNaMo06VVImIUPexR0GVA+OerHVDd43KqXCPQUSCgL5Ad+AUIA/4A1isqhs8G88wyvv0r0/5audX3H363bSNdd8JVRHh9IICTu85nj3Ze5i1eRYfbf2Iz3Z8Ruf4ztzc8Wa6N+jutu2lz5pN8b79nPL004ifZ7+f2evXJ/6u0SQ/+xzZn39OxOWXe3R7hvc74W+kiIwFlgPnAiuAKcBcoBh4VkS+EpGOng5pGAC7s3bzzIpn6FKnCyPbjfTYdhqEN+DBMx/k66u/5qEzHyIlL4VRS0exLX2bW9ZfkpVF6htvEHreeYSec45b1lmRmKQkgtq148BTT1OSmVkt2zS8V0VfVVaq6hmqep+qzlLVr1V1kaq+rKr9gETAt64YMmqk4tJixvw4BpvYePq8p7H5eX6ETVhAGEltk5jTZw7B/sG8se4Nt6w3ddpblGRmknDvPW5ZnyvE3596T46jJD2d5JderrbtGt7phIVBVReXfSwioUf9PFlVV3kimGGUNXXdVNalrOPRcx+lXli9at12VFAUiW0S+fLvL9mavrVK6yo6mEzae+8R0bdvtV9bENS2LTEjRpAxdy6HV5mPrXF8Lh3cFJGuIrIR2OR83ElEJnk0mWE4rUlew5R1U+jXtB+XN7Hm+Pi1ba8lxB7CG2urttdwaOJEtKSE+LtGuynZyYm/YxT2+vXZ/9jjlBYWWpLBqPlcPes1HrgMSAVQ1bXA+Z4KZRhH5Bbl8vAPD1MvtB6PnP2IZTmigqIY1noYX+78kj/T/6zUOgq2bydj/nyihwwhoKE1nV/9QkKoO/ZxCrdvJ3Xqm5ZkMGo+l4dDqOruo54qcXMWw/iXZ1Y8w77cfTzT/RnCAsIszTKi3QjC7GGV3mtIGf8KfkFBxN12q5uTnZyw7t2J6NuX1ClTKPjrL0uzGDWTq4Vht4h0BVREAkTkfpyHlQzDU774+ws+/etTbuxwI6clnGZ1HCIDI0lsk8hXO79iS9qWk3pt3po1ZH/1FTHXX4d/jPUN/Oo8PAYJCWH/44+jpaVWxzFqGFcLw63AKKA+sAfo7HxsGB5xIPcA434eR4e4Dtzaydpv2GUNbzucMHsYU9ZNcfk1qkryiy9hi40lduRIz4U7Cf6xsdR58AHyVq0m48MPrY5j1DCutsQIVtXEsk+IiGf7BRi1VqmW8t8f/0tRaRHPdH8Gu5+92rbdeMziCpcJiDubr4q+ounjb1JacEqFy595YBPjVq3i9Y5XsvjJbytc/u8gV5JWXeRVV5H56QKSX3iRsB49sCckVM+GjRrP1cKwQ0TmAderap7zuSXA6Sd6kYj0Al4FbMA0VX32OMtdDcwDzjTDX433NrzHigMreKLrE5wacWq1bvvvZ/tUuExWYXd6ffgrF/XYwCs9bzrhslpSwo4rp1J6aiNen/kEE+0uFLmxLoatIhGh7hNj2THgCg4+8wwNxo+vng0bNZ6rh5LWAz8AP4pIM+dzJ+zoIiI2YCJwOdAWGCoi/xq4LSLhwGgcV1YbtdzmtM28+vurXNToIq5sfqXVcY4pIiCCpLZJLN21lM1pm0+4bObChRT8+ScJd92FuFIUqllgkybE3XYr2Z99Tva331odx6ghXC0MqqqTcPwBXygi/YCKunGdBWxT1e2qWgjMAQYcY7kngeeBfBezGD4qvzifh75/iOjAaMaeO/ZfE+/UJEltkwi3hzN5zeTjLvNPo7x27Qjv1asa052c2BtuIKB5Mw6MG/dP62+jdnO1MAiAqi4HLgIeAFpX8Jr6QNkhrnucz/1vpSKnAQ1VdVGFAURuFpFVIrIq5chUhoZPeXn1y2zP3M7/nfd/RAVFWR3nhCICIhjedjjLdi9jU+qxB+ilz3Y0yku4/z6PN8qrCgkIoN64Jynet5+UCa9ZHceoAVz9be195I6q7gcuBCr6CnSsr3v/7GWIiB+OC+fucyWAqk5V1S6q2iU+Pt6Vlxhe5Ps93zN782yS2iTR9ZSuVsdxSVLbJMIDwpm89t97DSXZ2aROfoPQbt0IPfdcC9KdnJDTTyNq6BDS3n+fvPWeaTVueI+KuqsmOe8OFZF7j9xwHFLqUsG69wBlL+9sAOwr8zgcaA98KyJ/A+cAC0SkovUaPkZsOTy6/FFaRLfg7jPutjqOy8IDwrm27bV8s/sbNqZuLPezfxrl3XevRelOXsK99+IfG8v+xx5Di4utjmNYqKI9hiNN88KPczuRlUALEWkiIgHAEGDBkR+qaqaqxqlqY1VtDPwC9DejkmoXVSWo3ofkFObwXPfnCLQFWh3ppCS2SXTsNZQ511CUnEza9OlE9OlT7Y3yqsIWHk6dR/9LwaZNpE1/z+o4hoVOOFxVVac4/33iZFesqsUicgfwBY7hqm+r6gYRGQesUtUFJ16DURvM+3Me/uGbueeMh2gR3cLqOCctPCCcEW1H8Pqa19mQuoF2se04NHGSo1He3XdZHe+khV9yCWEXXkjKa68RftmlBDRoYHUkwwKudld9XkQiRMQuIktF5FCZw0zHpapLVLWlqjZT1aeczz12rKKgqj3M3kLtsj1zOy+sfIHinBYMazPM6jiVltgmkYiACCavmUzB9h1kfPgh0ddcY1mjvKoQEeo++l/Ez48DY58wU4HWUq6efL5UVbNwTPG5B2iJY2SSYVRKUUkRY74fQ5B/EPn7B+EnNXfUTkXCAsIY0W4E3+35jr+eH4dfYKDljfKqwl6vHvH33EPujz+StajiK8EN3+Pqp/HIlTm9gdmqmuahPEYt8fqa19mUtomxXceixRFWx6myYa2H0TklBPn2F2JuuB7/2FirI1VJ9LChBHXsyMFnnqEkI8PqOEY1c7UwLBSRzThGIi0VkXjMBWlGJa08sJJ3/niHgS0GclGji6yO4xah9lBG/RRBRggc6HuW1XGqTGw2x1SgWVkcfOEFq+MY1cylwqCqY4BzgS6qWgQc5thXMRvGCWUXZvPIj4/QKKIRD575oNVx3Cb3hx+I3LiHJT1Cmbz1HavjuEVQq1bEXncdmfM/IvcX07GmNjmZiXrSVbXEeT9XVQ94Lpbhq+b/OZ8DuQf4v27/R4g9xOo4bqElJSS/+BL2Ro1onHQTP+z9gXUp66yO5RZxo27H3qgRBx5/nNKCAqvjGNXEe8/4GV6npLSE2Ztn06VOFzondLY6jttkLVrkaJR3910M6ZBEVGAUk9b6xpTofkFB1Bv7OIU7d3LojarNd214jwoLgzh437g7o8b5ds+37MvdR2KbxIoX9hKlhYWkvDqBoLZtCe/Vi1B7KCPbjWT53uWsTVlrdTy3CO3alcgBA0h9cxr5f1ZuvmvDu1RYGNQxkPmTashi+LjZm2ZTN7QuPRr2sDqK22TMnk3Rvn3lGuUNbT2U6MDoE3Ze9TYJYx7CFhbGgcfMVKC1gauHkn4RkTM9msTwaVvTt7LiwAqGtBqCv5+r80PVbCXZ2Rya/AahXbsS2vV/jf9C7CGMbD+S5fuWsyZ5jYUJ3cc/OpqEMQ+Rt2YNGR98YHUcw8Nc/YT2BG51NrvLxdE5VVW1o6eCGb5l1uZZBNoCGdhioNVR3Cb1rbcoycgg/hiN8oa0GsL0DdOZvHYyUy5xfX7oGmV8B8jc9c/DSIXMOrEkPz2W8HW34x9YTXsOkY3gnvXVsy0DcL0wXO7RFIZPyyzIZNFfi+jTtE+Nn2fBVUXJyaS962iUF9yu3b9+HmIP4bp21/HS6pdYk7zGO0+2Z+6CsZn/PBSg7tatbO/Xn4zGzxF3y83Vk2NsZPVsx/iHq9cx7MTRQvtC5/3Drr7WMD7e+jH5JfkMa+29/ZCOdmjSJLS4mPi7Rh93mcGtBhMTFMOkNb4xQgkgsEULQs49h/Q5c0xrbh/mahO9x4GHgIedT9mBGZ4KZfiOktIS5myZQ5c6XWgV08rqOG5RsGMHGfOcjfIaNTruciH2EK5vfz0/7/+Z35N/r8aEnhWTlETx/v1kL11mdRTDQ1z91n8l0B/H+QVUdR8Vz8dgGHy35zv25uz16u6pR0t55VVHo7zbb6tw2UEtBxETFMPENROrIVn1COvRA3v9+qTPMN8NfZWrhaHQOWxVAUQktILlDQOAWZtmUTe0Lj0b9rQ6ilvkrVtH9hdfEHO9a43yjuw1rNi/gtUHV1dDQs8Tm43oYUM5vHIl+Vu2WB3H8ABXC8NcEZkCRInITcDXwDTPxTJ8ga8NUVVVkl98CVtsLDEjR7r8usGtBhMbFOtT1zVEDRyIBAWRPmOm1VEMD3D15POLwIfAfKAV8JiqTvBkMMP7zd4826eGqOb++COHf/2VuNtuwxbm+k5zsH+wY6/hwApWHfCNuahsUVFE9utL5sKFpi23D3L15PNzqvqVqj6gqver6lci8pynwxneK7Mgk4V/LfSZIapaWupolNewIdGDB5306we3GkxccByT1/rOXkN0UhKan0/G/I+sjmK4mauHki45xnPm2gbjuHxtiGrWokUUbNlC/N13IQEBJ/36IP8gbmh/A78e+JWVB1Z6IGH1C2rVipAuXUifNQstKbE6juFGJywMInKbiKwHWonIujK3HYBv9BU23O7IENUz6pzhE0NUSwsLSXnlVYLatiXi8sp/H7q65dXEB8f73F5D0d695Hz3ndVRDDeqaI9hFtAPWOD898jtDFVN8nA2w0sdGaLqK11UM+bNo2jfPuLvu/efRnmVEeQfxA0dbmDlgZU+s9cQfvFF+Neta05C+5gT/paraqaq/g38FzjgvOq5CZAkIt5/4NjwCF8aoqqlpaS99x7BnToR1q1bldd3ZK9h4pqJOEaAezfx9yd6yBByf/qJgu3brY5juImrX3/mAyUi0hx4C0dxmOWxVIbX2pa+jRUHVnBNq2t8Yohq7g8/ULRzF9HDh7tlfYG2QG7ocAOrD672mb2GqMGDELvd7DX4EFcLQ6mqFgNXAa+o6j1APc/FMryVr3VRTZsxE1t8HBGXHmv8ReVc3fJqEoITfGavwT8mhojevcn85BNKcnKsjmO4gauFoUhEhgLXAoucz9k9E8nwVpkFmSzavojeTXoTHRRtdZwqK9ixg9wffiD6miGVGol0PEf2Gn5L/o1fD/zqtvVaKTopidLDh8n86GOroxhu4GphuA44F3hKVXeISBNMEz3jKJ9s+4S84jyf6YuUPms22O1EXzPY7ese2HIgCSEJTFozySf2GoI7tCe4UyfSZ840M7z5AFevfN4I3A9sEJEOwF5VfdajyQyvUlJawuzNszmjzhm0jmltdZwqK8nJJfOjj4jo1Qv/+Hi3rz/QFshNHW7it+TfWHFghdvXb4XopCQKd+4kd/lyq6MYVeTqlc99gL+ACcDrwDYRMRe4Gf/4p4uqj1zQlvnpJ5Tm5hKT5Lkht1e1uIrYoFje3/i+x7ZRnSIuuxRbfBxppuuq13P1UNJLQE9V7aGqF+CY6nN8RS8SkV4iskVEtonImGP8/FYRWS8ia0TkRxFpe3LxjZpi1uZZ1Ampw4WNLrQ6SpVpaSnpM2YS1KEDwZ06eWw7AbYABrcazPd7vmdn1k6Pbae6SEAA0YOvIff7Hyjc6f3vpzZztTAkq+q2Mo+3A8kneoGI2ICJOFpntAWGHuMP/yxV7aCqnYHngZddzGPUINvSt7Fi/wqGtPaNLqq5P/1M4Y4dHt1bOGJQy0H4+/kzZ/Mcj2+rOkRdMxhsNtJnmdHs3qyilhhXichVOM4tLBGRkSIyAlgIVDQI+yxgm6puV9VCYA4woOwCqppV5mEozvkeDO8ya/MsAvwCfGaIavqMGdhiYwmvQvsLV8WHxHNZ48v4eNvH5Bblenx7nmZPSCDissvImP8Rpbne/35qq4r2GI60wAgCDgIXAD2AFKCi8Yj1gd1lHu9xPleOiIwSkb9w7DEcfwJdo0Y6MkS1T9M+PjFEtXD3bnK++46owYPwc+MQ1RNJbJ1IblEun277tFq252nRSYmU5uSQuWCB1VGMSjrhfr+qXleFdcuxVnmMbUwEJorIMBytN0Ycc2UiNwM3AzQ6wTy7RvXyuSGqM2eBzUb0kCHVts0O8R3oGN+R2ZtnM6T1EJeP79ZUwZ07E9SuHWkzZxI1ZAgix/pTYNRkro5KaiAiH4tIsogcFJH5ItKggpftARqWedwA2HeC5ecAVxzvh6o6VVW7qGqXeA8MHzRO3pEhqqcnnO4TQ1RLi4SM+fOJuPQS7HXqVOu2E1sn8nfW3/y076dq3a4niIhj6Oq2vzj8yy9WxzEqwdUvJ+/g6LB6Co7DQQudz53ISqCFiDQRkQBgiHMd/xCRFmUe9gG2upjHqAG+3/O9T3VRzdwZTGl2NtFJ1d84+JJTLyE+OJ4Zm3xjqGdE78uxRUeTZvoneSVXC0O8qr6jqsXO27vACb+2O3sr3QF8AWwC5qrqBhEZJyL9nYvdISIbRGQNcC/HOYxk1EwzN8/0nSGqqqRvDSWwbRuCTzut2rdvt9kZ3Gowy/cuZ4fd+0d2+QUGEjV4MDnffEPhnr1WxzFOkquF4ZCIJImIzXlLAlIrepGqLlHVlqraTFWfcj73mKoucN6/S1XbqWpnVe2pqhsq/1aM6uRrQ1QPr1hBQaadmMQky46JX93yaux+dmaHh1uyfXeLHnINiJAxZ7bVUYyT5GphuB4YDBwA9gNXO58zaqnZm2f71BDVtBkzsAWUENGnt2UZ4oLjuLzJ5XwaHkp2YbZlOdzFXq8e4RddRMa8DynNz7c6jnESXO2VtEtV+6tqvKomqOoVzkl7jFoosyCThdsX+s4Q1T17yVn2DVHNDuMXFGRplmFthnHYz8+nhq6WZGaStWhRxQsbNcYJjwGIyGuc4KIzVTXXHdRCPjdEdfYsECG6ufUXZLWLbcdp+fnM2jyLYW2G4SfePXg15MwzCWzZkrQZM4kcONAMXfUSFf3WrQJWn+Bm1DI+N0Q1L4+MD+cTftFF2ENrRrvoYVk57M7ezY97f7Q6SpU5hq4mUrB5M3mrzZ8Mb1HRnM/TT3SrrpBGzeFzQ1QXLaI0M5PoauiL5KqLcg+TEJLAzE2+MdQzsl8//CIjzdBVL1LRoaRrXVzPGlVd54Y8Rg3nU11UVUmfMZPAli0JOfNMWGJ1Igc7MKTVECb8PoHtGdtpGtXU6khV4hccTNTAgaRNn07RgQPY69a1OpJRgYoOJTVx8RblwYxGDfFXxl/8sv8XnxmimrdqFQVbthCdlFjjjn0PbDmQAL8AZm32jS6l0cOGQmkp6XN8o4usr6uoV9IT1RXEqPlmbfKtLqppM2biFxlJZL9+Vkf5l5igGHo37c2CvxYw+vTRRAREWB2pSgIaNCCsZ08y5s4j7rbb8AsMtDqScQIutd0+3q26QhrWyyrMYuH2hfRu2tsnhqgW7d9P9tdfE3X1QPyCg62Oc0yJbRLJK87j460fWx3FLWKSEilJSyPrs8+sjmJUoKLjASf6KqXAR27MYtRgH2/92DFE1QNTd9aPCqbxmMVuX++JjNi4hEElpfTbVZdk57b/tvYShn9pHdOaM+qcwezNs0lqk4TNz2Z1pCoJOfdcApo1I33GTCIHDKhxh++M//Fk223DR5Qdotomto3b1798TPWeyC4tKGBbj6cIvuhCfn1t+P9+MLZaY7gksU0i9357L9/v+Z6ejXpaHadKRIToxGEcHPck+WvXEty5s9WRjONw+eoZEekjIg+KyGNHbp4MZtQcP+z9gb05e33mgrasxUsoSU+vlqk7q6pnw57UDa3LzM2+MdQzasAA/MLCzNDVGs7V+RjeAK4B7sQxAc8g4FQP5jJqkJmbfKuLatqM9wlo3oyQc86xOk6F/P38GdJqCCv2r2Bruvd3pfcLDSXyqivJ+uILilNSrI5jHIerewxdVfVaIN05Uulcyk/CY/ioI0NUr2l1DXY/u9Vxqizv998p2LiJmCTruqierIEtBhJoC/SdoatDh0JREekfzLU6inEcrhaGPOe/h0XkFKAIx/ULho/7p4tqS98Yopo+YwZ+4eE1cojq8UQFRdG3aV8W/bWIzIJMq+NUWWCTJoR27076B3PQwkKr4xjH4OpVSotEJAp4AfgNx4ikNz2WyqgRsgqzWPDXAno37U1MUIzVcaqs6GAyWV9+RUxiIn6hoVbHOa49GkeDsZHlnhtmtzO/QT0+mtSO6zKrqSV3pOfmVo9JSmT3LbeS9dVXRPbp47HtGJXjUmFQ1Sedd+eLyCIgSFW9/6uLcUKeHKJqhYwP5kBJCdGJNfv9nFcwgb+fLf/HsiVw1hc3MDt0N8PvWuL1V56Hdu+O/dRGjqGrpjDUOK6efF4rIo+ISDNVLTBFwfeVlJYwZ/Mcjw1RrW6lhYWkfzCXsPPPJ6CR574Je9KwNsPYn7uf73Z/Z3WUKhM/P2KGDSPv99/J22AmbqxpXD3H0B8oBuaKyEoRuV9EvPPTZbjkh70/sCdnj88MUc3+/HNKUlOJTkqyOkql9WjQg1NCT/GZoauRV12FhISQboau1jiuzuC2U1WfV9UzgGFAR2CHR5MZlnp3w7s+M0QVIO39GQQ0aUJot65WR6k0m5+Noa2HsvLASrakbbE6TpXZwsOJHNCfrMWLKU5LszqOUcbJXODWWEQeBOYArYEHPZbKsNTKAytZfXA117W/zjeGqK5dS/769UQnJiJ+3j0j2pUtriTbcVh5AAAgAElEQVTYP9hnhq7GJCaihYVkzPvQ6ihGGa6eY1iBoy+SDRikqmep6kseTWZYZsq6KcQGxfpWF9XQUCKvuMLqKFUWGRhJ36Z9Wbx9Men56VbHqbLA5s0JOfcc0ufMQYuLrY5jOLn69WmEqp6uqs+o6naPJjIstSZ5DSv2r+C69tcR5F/DuspVQnFKClmff07klVdiC6u5Q1RPxrDWwygoKWD+1vlWR3GLmKQkivfvJ3vpMqujGE4Vtd0+cqaut4jce/StGvIZ1eyNdW8QHRjNoJaDrI7iFulz50JRUY0fonoymkc355x65zBn8xyKS73/W3ZYjx7Y69cnfcYMq6MYThXtMRz5ihV+nJvhQ9anrGf53uWMaDeCEHuI1XGqTAsLyZjzAaHduxPYxLcu1E9sk8jBwwdZtsv7v2WLzUb0sKEcXrmS/C3ef1LdF5ywMKjqFOfdSar6xNG3ashnVKMp66YQGRjJkNZDrI7iFllffkVxSopXdFE9Wd3rd6dBWANmbvKNoZ5RAwciQUGkvf221VEMXD/H8JOIfCkiN4iI90/fZfzLptRNfLfnO4a3GU6o3TeOxafPmIH91EaEdu9udRS3OzJ09bfk39iUusnqOFVmi4oieuhQMhcuomC7OY1pNVevY2gB/BdoB6wWkUVlzj8YPmDKuimE28N95oK2vD82kLdmDTHDhnn9ENXjuaLFFQT7B/vMXkPsTTciQUEcev11q6PUei5/YlT1V1W9FzgLSAOmeyyVUa3+TP+TpbuWktg2kfAA3zh1lD5jBhISQuRVvjs1eURABP2b9WfJjiWk5qVaHafK/GNiiBk+nKwln5lzDRZz9TqGCBEZISKfAT8B+3EUiIpe10tEtojINhEZc4yf3ysiG0VknYgsFREz+Y8Fpq6bSoh/CEltfGMnsDg1lazFi4kc0B9buG8UuuMZ1mYYRaVFPjN0Nfa6kfiFhZHy2mtWR6nVXN1jWAt0BsapaktVfUhVV5/oBSJiAyYClwNtgaEi0vaoxX4HuqhqR+BD4PmTSm9U2faM7Xz595cMazOMyMDIil/gBTLmzUOLiohJ9L2TzkdrGtmUbqd044PNH1BUWmR1nCqzRUURc91Icr5eSt4fprmeVSosDM4/8B+r6j2q+vNJrPssYJuqblfVQhytNAaUXUBVv1HVw86HvwANTmL9hhtMXT+VIP8ghrcdbnUUt9CiItJnzyG067kENm9udZxqMazNMJLzklm6c6nVUdwiZsQIbJGRpLw2weootVaFhUFVS4BOlVh3fWB3mcd7nM8dzw3AZ5XYjlFJYj/EZzs+45pW1/jERDwA2UuXUnzwoFd3UT1Z59U/j0bhjXzmJLQtLIyYG24g97vvOfz771bHqZVcPZS0RkQWiMhwEbnqyK2C1xxrQl095oKOEU5dcMwQd+yVidwsIqtEZFWKmUTcLQLjvsHuZ2dEuxFWR3GbtBkzsNevT9gFF1gdpdr4iR/D2gxjTcoaNhzyjcMvMUmJ2GJjSZlg9hqs4GphiAFSgQuBfs5b3wpeswdoWOZxA2Df0QuJyMXAf4D+qlpwvJWp6lRV7aKqXeLj412MbRzPnuw9+Ef+zqCWg4gLjrM6jlvkb95M3qrVRA8bhthsVsepVgOaDSDEP8Rnuq76hYQQe9ONHP75F3KTA6yOU+u4eh3Ddce4XV/By1YCLUSkiYgEAEOABWUXEJHTgCk4ikJyZd6AUTnT1k8D9eO69tdZHcVt0mbMQIKCiBrou0NUjycsIIwrml/BZzs+41DeIavjuEX0kCH4JySQsj4c1WMebDA8xNXhqu+IyNtH3070GlUtBu4AvgA2AXNVdYOIjBOR/s7FXgDCgHkiskZEFhxndYYb7c/Zz6d/fUpRRhcSQhKsjuMWhTt3kvnpAiKvGIAtKsrqOJYY2nooRaVFzPtzntVR3MIvKIjYW24mLyWQ3OU/WR2nVnH1UNIiYLHzthSIAHIqepGqLnEOb22mqk85n3tMVRc471+sqnVUtbPz1v/EazTc4a0/3gKgMLWHtUHcKPnFlxC7nfhRo6yOYpnGkY3pXr87c7fMpajE+4euAkQNGoR/SDEpEyaYvYZq5OqhpPllbjOBwUB7z0YzPOFg7kE+2voRA5oNQIt945t17q+/kv3VV8TdfBP+tfz8U2KbRA7lHeLLnV9aHcUt/AICiGuXQ/66deR8+63VcWqNyjaRaQE0cmcQo3q8u+FdSrWUGzvcaHUUt9DSUpKffQ7/evWIGTnS6jiWO/eUc2kc0ZhZm3zjJDRAVJPD2Bs2JGXCa2hpqdVxagVXzzFki0jWkRuwEHjIs9EMdzuUd4h5f86jb9O+NAj3jWsJMxcsIH/jRhLuvQe/IO+fca6qjgxdXXdoHetS1lkdxy3ED+LvGEXBpk1kf/W11XFqBVcPJYWrakSZW0tV9Y3mLLXI9A3TKSot4qaON1kdxS1KDx8mZfwrBHXoQESfPlbHqTH6N+tPmD2M9ze+b3UUt4no25eApk059PpraEmJ1XF8nqt7DN1EJNR5P0lEXjYN77xLWn4aH2z5gMubXM6pEb7xX5f6zjsUHzxInTEP+Wxr7coItYcytPVQPv/7c9Ykr7E6jluIzebYa9i6jawlpkGCp7n6aZoMHBaRTsCDwE7gPY+lMtzu/Y3vk1+cz80dbrY6ilsUHUwmddpbhF92GSFnnGF1nBrnxg43UiekDk+teIqSUt/4hh3eqxeBLVty6PXX0WLvn+u6JnO1MBSrY6zYAOBVVX0VM+ez18gsyGTWpllc1vgymkY1tTqOW6S8+ioUF5Nw371WR6mRQuwhPHDmA2xO2+wz1zWInx/xo+90XLOyYKHVcXyaq4UhW0QeBpKAxc6Oq3bPxTLcacamGRwuPuwz5xbyN24k8+OPiR4+nIBGZnDc8Vx66qWcXfdsJvw+gbT8NKvjuEXYRRcR1K4dhyZNQgsLrY7js1wtDNcABcANqnoAR5fU4za8M2qOrMIsZm6cycWNLqZldEur41SZqnLwueexRUYSd+stVsep0USER85+hLyiPF797VWr47iFiBA/+k6K9uwh46OPrY7js05YGEREAFT1gKq+rKo/OB/vUtX3yi5j1EyzNs0iuyibmzv6xrmFnGXLOLxiBXF33oEtIsLqODVe06imJLVN4qOtH7E+Zb3Vcdwi9PzzCe7cmUNvvEFpwXH7bhpVUNEewzcicqeIlNtfF5EAEblQRKYDvtOz2cfkFuXy/sb36dGgB21i21gdp8q0sJDk518goGlTogcPtjqO17i1063EB8f7zIloESH+rtEUHzhAxlzfOH9S01RUGHoBJcBsEdnvnJ95O7AVGAqMV9V3PZzRqKTZm2eTVZjFLZ1845BL+pw5FO7cSZ2HHkTs5hSXq0LtodzX5T42pG7g422+cfgl5JxzCDnzTA5NmUJpXp7VcXzOCQuDquar6iRV7YajBcZFwOmqeqqq3qSqvjFI2gcdLjrMexveo1v9brSP8/62ViUZGaRMnERo166Enn++1XG8Tu8mvTmjzhm8+turZBZkWh2nyo7sNZQcOkT6rNlWx/E5FZ1jCBKRu0XkdeA6IEVVM6onmlEV8/6cR3pBOrd2vNXqKG5xaPJkSrOzSXjoIcxprZN35ER0dmE2E37zjVnRQrp0IbRbN1KnTaMkJ9fqOD6lokNJ03FMubke6A285PFERpXlF+fzzh/vcHa9s+mc0NnqOFVWsGMHaTNnETVwIEGtvH9klVVaRrdkaOuhzPtzHhtTN1odxy3i7xpNSXo66TN8p/1HTVBRYWirqkmqOgW4GuheDZmMKpq/dT6p+anc0tE3zi0kv/QSfgEBxI++0+ooXu/2zrcTHRTNUyueolS9v1NpcMeOhPXsSerb71CSlWV1HJ9RUWH4Z7YP54xsRg1XUFLA2+vf5ow6Z3Bm3TOtjlNluSt+JefrpcTefHOtn2vBHcIDwrn3jHtZl7KOT7d9anUct4gffSelWVmkvTvd6ig+o6LC0KlMu+1soOOR+87220YN88nWT0jOS+bWTt5/bkFLSjj43LPOuRbMqGh36desH53jO/PKb6+QVej9H+OgNm0Iv/RS0qZPpzg93eo4PqGiUUm2Mq22w1XVv8x9c3VRDVNUUsRbf7xFp/hOnF33bKvjVFnmpwso2LiJhHvvNXMtuJGf+PHI2Y+QUZDBxN8nWh3HLeLvvIPSw4dJe/uEU9EbLjK9in3Igr8WsD93P7d2utXrR+445loYT1DHjkT06W11HJ/TJrYNg1oOYs6WOWxJ22J1nCoLbNGCiD59SJsxk+JDh6yO4/X8rQ5guEdRaRFvrn+TdrHt6HZKN6vjnJzxHSBzV7mnUteHU5wSTv2Om5Bx0dWTI7J2NeS787Q7+fLvL3l6xdO82+tdr/8yETfqdrKWLCH1zWnUeXiM1XG8mikMPmLJ9iXszdnLmLPGeN8HPHMXjP3fRVdFBw+Selkvwnv1IOSV8RYG822RgZHcdfpdjP15LIu2L6Jfs35WR6qSwCZNiBwwgPTZs4m5/jrsdepYHclrmUNJPqCktIQ3179J65jWXNDgAqvjVFnK+FegpISE+++zOorPu7LFlXSI68DLq18mpzDH6jhVFjfqdrS0lNQpU6yO4tVMYfABn//9OTuzdnJLx1u8b2/hKHl/bCDzk0+IGXEtAQ0aWB3H5x05EZ2al8rktZOtjlNlAQ0aEDVwIOnzPqRo716r43gtUxi8XKmWMnXdVJpHNefCRhdaHadKVJXk557DFh1N7C2+cXGeN2gf156rWlzFzE0z2Za+zeo4VRZ36y0IkDLZ+wudVUxh8HJf7fyK7ZnbuaXjLfiJd/93Zn/9NYdXriR+9J3Yws3MsdXprtPvIiwgjGd+fQbHLL7ey16vHlFDhpD58ScU7txpdRyv5N1/SWq5Ui1lyropNIlswiWnXmJ1nCrRwkKSX3iRgGbNiBo0yOo4tU50UDSjTxvNrwd+5Yu/v7A6TpXF3XwTYrdzaNIkq6N4JTMqyYtlFmQSHRjNFc2vwOZnszpOlaTNmkXRrl00nDoF8Te/llYY2GIgH/75IS+seoHzG5zPJS/9wt4M6+c6+DkonnpjI0/qNf5AdOMI0hZ8SmzhNAIjq9jRJ7IR3OMbM+C5wnwCvVh0UDRvXfaW1+/6FxcIhyZNJrRbN0K7mz6NVrH52Xjk7EcY/tlwpqybwt6Mlvz9bB+rY9F4DJXKEZuWRvrFl5DC5TQYW8VhzydZmLydRw8liUgvEdkiIttE5F9XnIjI+SLym4gUi8jVnsziy7x9JNKhDeGU5uSQ8OCDXv9evF3nhM5c0fwK3tv4HhKQYnWcKvGPiSHm2uFkf/Y5+Vu8/+ru6uSxwiAiNmAicDnQFhgqIm2PWmwXMBKY5akcRs1WsH0H6VtDibr6ajPXQg1x9+l3E2wLJqjOAq/fG4297jr8wsNJeXm817+X6uTJPYazgG2qul1VC4E5wICyC6jq36q6DvD+xvBGpSS/+CJ+NjVzLdQgscGxjDptFP5hW1m6a6nVcarEFhlJ3O23k/Pdd44LJw2XeLIw1Ad2l3m8x/lcpYjIzSKySkRWpaR49y6u4ZD7yy/kLFtGbNsc/OPirI5jlHFNq2soya/L8yufJ6/Y+hPQVREzcgRR11xD6tSppL33ntVxvIInC8OxDhZXel9OVaeqahdV7RJvJmzxeo65Fp7H/5R6xLT0/lYMvsbfz5+CAwPYn7ufaeunWR2nSkSEuo89Svgll3Dw6WfIXLzY6kg1nicLwx6gYZnHDYB9Htye4UUyP/mEgk2bSLjvPvzM2LgaqSSvCX2b9uWdP95hV9auil9Qg4nNxikvvkBIly7sG/MwOcuXWx2pRvPkR3Il0EJEmgB7gSHAMA9uz6t0e3aZ5WPE60cFW7Ld0txckl95heBOnYjo3dvxm2L8o35UMI3HWP+ttn5UMPeecS/f7P6GZ399lokXTfTqUWN+gYE0mDSRncOvZe+do2k0fTrBHdpbHatG8lhhUNViEbkD+AKwAW+r6gYRGQesUtUFInIm8DEQDfQTkSdUtZ2nMtUkezPyasQYcSukvvUWJSmHSJgwwav/0HjK8jE1q+fVbZ1u48VVL/Ldnu/o0bCH1XGqxBYRQcOpU9k5bBi7b7mFxrNmEtC4sdWxahyPXsegqktUtaWqNlPVp5zPPaaqC5z3V6pqA1UNVdXY2lIUarOCbdtIffsdInpfTshpp1kdx3DBsDbDaBbZjGd/fZb84nyr41SZvU4CDae9CarsuvEmipKTrY5U45heSUa1yVv/BzuThuMXFkbCfWauBW9h97PzyNmPsDdnL+/88Y7VcdwisEkTGk6dQnFaGrtvvoWS7GyrI9UopjAY1SJ3xa/sGjkSv5AQGs+cgb1+pUcuGxY4q95Z9Grci7f+eIs92XusjuMWwR060GDCBAq2bWPPqDsoLSiwOlKNYQqD4XHZy75h90034V+3LqfOnkXAqadaHcmohPu63Ief+PH8yuetjuI2Yed145RnnuHwr7+y74EH0ZISqyPVCGagYG02voNjvmUPyvw7mH0rogiKKqJhx+X4Tz5G24vIRh7NYLhH3dC63NrpVsavHs/TK57mwTMfxN8HxhpH9utLSVoqB595lgNPPkndxx+v9YMivP9/1ai8zF0wNtNjq0+bNYuDH/wfIV260GDyJGxhYR7bllE9RrQdQXp+Ou9ueJdd2bt48fwXCQvw/v/XmBEjKD50iNQ3p+EfF0/8HaOsjmQpcyjJcDtV5dCUqRwc9yRhF1xAwzenmqLgI2x+Nu7rch9jzx3Lin0rGP7ZcPbm+MbcyvH33kvkFVdw6PXXSZ8zx+o4ljKFwXArVSX5xRdJGT+eiL59afDaBPyCgqyOZbjZwJYDmXzJZA4ePsiwxcNYk7zG6khVJiLUe3IcYRdcwIEnxpH1xZdWR7KMKQyG22hJCQcee5y0t94methQTnn+OcRutzqW4SHn1DuHGb1nEGoP5YYvbuCzHZ9ZHanKxG6n/ivjCe7YkX3330/uil+tjmQJUxgMt9DCQvbefz8Z8+YRe8st1Hn0UcTP/Hr5uqaRTZnZeybt49rz4PcPMnnNZK+f98AvOJgGb0zG3qgRe0aNIn/zZqsjVTvzyTWqrDQvj92j7iD7s89JeOABEu65u9aP6qhNooOiefPSN+nfrD+T1k5izA9jKCjx7msC/KOjaTTtTfzCwth1000U5nj3nOonyxQGo0pKsrLYdcON5P74I3WfHEfsDddbHcmwQIAtgP/r9n/cdfpdLNmxhBu+uIHUvFSrY1WJvV49Gk17Ey0sYte3sRSnevf7ORmmMBiVVpyays4RI8lbv576418metAgqyMZFhIRbuxwIy9d8BJb0raQuCSRbenbrI5VJYHNm9PwjckU5/mx+5ZbKcnJtTpStTCFwaiUon372JmYROGOHTScNJGIXr2sjmTUEJc2vpR3er1DQUkBwz8bzvK93j33Qchpp1G/azr5mzaxd/RotLDQ6kgeZwqDcdIKtu/g78QkilNTafTWNMK6d7c6klHDtI9rz+w+s6kfVp9RS0cxZ7N3XxcQXr+Aek8+Se5PP7Hv4UfQUt+ept4UBuOk5G/cyM6kJLSggFPfm07IGWdYHcmooeqG1mX65dM5r/55PLXiKZ5Z8QzFpcVWx6q0qKuuJP6+e8lavJiDzz7r9aOvTsQUBsNlh1etYue1I5DAQE6dOYOgNm2sjmTUcKH2UF7t+SrD2w5n1uZZ3LnsTnIKvXeO79gbbyRmxLWkv/c+qW9691zYJ2IKg+GSnO+/Z9eNN+EfH0/jWTMJbNLE6kiGl7D52XjwzAd59JxH+Xnfzwz/bDj7crxz+ncRIeGhh4jo04eUl18mY/5HVkfyCFMYjAplLVnC7ttHEdC0CafOeB97vXpWRzK80OBWg5l08SQO5h5k6OKhrE1Za3WkShE/P0555mlCu3Zl/2OPkfOjd59cPxZTGIwTSp87l7333U9wp06cOn06/rGxVkcyvFjXU7oyo/cMQvxDuP7z6/l8x+dWR6oUCQig/oQJRA0cSFC7tlbHcTtTGIzjSp02jQOPPU5o9/NoNO1NbOHhVkcyfEDTqKbM6jOL9nHteeD7B5i81jvbaNjCQqk37gn8o6OtjuJ2pjAY5WhREdnffMOeO0eT/OJLRPS+nIavv45fcLDV0QwfcqSNRt+mfZm0ZhIP//iw17fR8CVmoh4DVSVvzRqyFi4ka8lnlGRkYIuKIvbWW4i/807EVrv6xBjVI8AWwNPnPU3jiMa8vuZ19mbv5cluT9I4srHV0Wo9UxhqsYIsG1kTJpC5cBFFu3cjgYGEX3QhEf36EXbeeaZltuFxIsItnW7h1MhT+e+P/6XfJ/1oHtWcng17clGji2gb63vH772BKQy1TPGhQ2QtWULmgoXk/1EH5A1Czz2HuNtvJ/ySi81Ma4YlejXuRef4zizdtZRlu5bx9h9v8+b6N6kTUofAOk35eV8MXep2we5nvqxUB1MYaoHS3Fyyly4lc8FCcn/+GUpKCGzbhoTOmUS8ugZ7nQSrIxoGdUPrktgmkcQ2iWTkZ/D93u9ZtmsZB3K+5+avfibcHs75Dc+nZ8OenFf/PELtoVZH9lmmMPgoLS4m96efyFywkOylS9G8POynnELsjTcS2a8vgc2bw9hIMEXBqIGigqLo36w//Zv1p/HDH/PObVEs272Mb3d/y+Lti7H72Tmn3jlc2OhCejTsQVxwnGcDRTZyfF6sFtkI7lnv8c2YwuBDVJX89evJXLCQrM8+oyQ1Fb/ISCL79yeyX1+CTz/dzKpmeB8NoGejnvRs1JPi0mLWJK9h2e5lLNu1jB9+/oFxP4+jU3wnLmx0IRc2upBTI051f4Zq+GPskmoqTqYw+IDCnTvJXLiIzIULKNq5CwkIIKxnTyL79SX0/PPxCwiwOqJhuIW/nz9d6nahS90uPNDlAf5M/5Nlu5fxza5veHn1y7y8+mWaRTb7p0i0jW2Ln5gvQyfLo4VBRHoBrwI2YJqqPnvUzwOB94AzgFTgGlX925OZfElJRga7brmF/LXrQISQs84i7uabCb/kEmwREVbHMwyPEhFaxbSiVUwrbut0G/ty9vHN7m/KnbxOCEmgZ8OenFX3LOJD4okJiiEmKIYwe5iZfvYEPFYYRMQGTAQuAfYAK0VkgapuLLPYDUC6qjYXkSHAc8A1nsrka/wiI7En1CHi/vuI6NPH9DAyarVTwk7518nrb3Z9w4K/FvDBlg/KLWv3s/9TJGKCY4gNiv3f46Oejw6KJtAWaNG7soYn9xjOArap6nYAEZkDDADKFoYBwFjn/Q+B10VE1Buvj7eAiNDgtQlWxzCMGqfsyev84nx2ZO4gLT/tn1tqfippeWmkF6STlpfGjowdpOanHvfq6zB72L+KxpH7PRv25JSwU6r5HXqWJwtDfWB3mcd7gLOPt4yqFotIJhALHPJgLsMwapEg/yDaxFY8d4iqklec5yga+Wmk5aX9u5jkp7E7ZzdrU9aSXpBOqZbSJLKJKQwn4VgH8I7eE3BlGceCIjcDNzsf5ojIlkrmiqOGFB55rsqrqPp7eaJGHGetMf8nbmDeiwdU8bPi0ffRjW6eWvWxxPGEVPa9uDxcy5OFYQ/QsMzjBsDRs3McWWaPiPgDkUDasVamqlOBqVUNJSKrVLVLVddTE/jKe/GV9wHmvdREvvI+oPreiyfHca0EWohIExEJAIYAC45aZgEwwnn/amCZOb9gGIZhLY/tMTjPGdwBfIFjuOrbqrpBRMYBq1R1AfAW8L6IbMOxpzDEU3kMwzAM13j0OgZVXQIsOeq5x8rczwcGeTLDMVT5cFQN4ivvxVfeB5j3UhP5yvuAanovYo7cGIZhGGWZa8UNwzCMcmplYRCRJ0VknYisEZEvRcQrByGLyAsistn5Xj4WkSirM1WWiAwSkQ0iUioiXjeCRER6icgWEdkmImOszlMVIvK2iCSLyB9WZ6kKEWkoIt+IyCbn79ZdVmeqLBEJEpFfRWSt87084dHt1cZDSSISoapZzvujgbaqeqvFsU6aiFyKYyRXsYhjpLeqPmRxrEoRkTZAKTAFuF9VV1kcyWXO9i9/Uqb9CzD0qPYvXkNEzgdygPdUtb3VeSpLROoB9VT1NxEJB1YDV3jj/4s4GjuFqmqOiNiBH4G7VPUXT2yvVu4xHCkKTqEc56K6mk5Vv1TVYufDX3BcK+KVVHWTqlb2okWr/dP+RVULgSPtX7ySqn7Pca4n8iaqul9Vf3PezwY24ei24HXUIcf50O68eezvVq0sDAAi8pSI7AYSgccqWt4LXA98ZnWIWupY7V+88g+QrxKRxsBpwAprk1SeiNhEZA2QDHylqh57Lz5bGETkaxH54xi3AQCq+h9VbQjMBO6wNu3xVfQ+nMv8ByjG8V5qLFfei5dyubWLUf1EJAyYD9x91NECr6KqJaraGceRgbNExGOH+Xx2oh5VvdjFRWcBi4HHPRin0ip6HyIyAugLXFTTrxo/if8Tb+NK+xfDAs7j8fOBmar6kdV53EFVM0TkW6AX4JEBAj67x3AiItKizMP+wGarslSFcyKkh4D+qnrY6jy1mCvtX4xq5jxh+xawSVVftjpPVYhI/JFRhyISDFyMB/9u1dZRSfOBVjhGwewEblXVvdamOnnOViKBOGa/A/jFG0dXAYjIlcBrQDyQAaxR1cusTeU6EekNvML/2r88ZXGkShOR2UAPHF1JDwKPq+pbloaqBBE5D/gBWI/jsw7wiLMjg1cRkY7AdBy/X37AXFUd57Ht1cbCYBiGYRxfrTyUZBiGYRyfKQyGYRhGOaYwGIZhGOWYwmAYhmGUYwqDYRiGUY4pDIZhGEY5pjAYXkFEcipe6p9l3xWRq11cNkpEbq98shOu+1YRudYN6zlNRKa5YT8wCpcAAAO6SURBVD0dROTdqq7H8H2mMBi1XRRwUoVBHCr87KjqG6r6XqWT/c8jOC7+qxJVXQ80EJFGVY9k+DJTGAyvJSKnishS50RFS4/6g3exiPwgIn+KSF/n8u2ck52scb6mBfAs0Mz53AvO5R4QkZXOZZ5wPtfYOeHLJOA3yvdGQkSeFZGNzte86HxurIjcLyKnONd/5FbizB4vIvOd21opIt2O8R7DgY6qutb5OExE3hGR9c5tDXQ+nyMiz4nIamezwrNE5FsR2S4i/cusciGOlh2GcXyqam7mVuNvQM4xnlsIjHDevx74xHn/XeBzHF98WuBocheE41t3onOZACAYaAz8UWadl+KYcF2cr18EnO9crhQ45xg5YoAt/K+TQJTz37E4Jh0qu+woHO0MwNHA8Tzn/UY4evocve6ewPwyj58DXinzONr5rwKXO+9/DHyJo2d/JxztRY4s3w1YaPX/p7nV7JvPdlc1aoVzgauc998Hni/zs7mqWgpsFZHtQGvgZ+A/ItIA+EhVtzr6rJVzqfP2u/NxGI7isgvYqceeMSsLyAemichiHMXkX5x7BDcC3Z1PXQy0LZMhQkTC1TGpzBH1gJQyjy+mzDd+VU133i3EUQzB0RuoQFWLRGQ9jqJ2RDLglVPZGtXHFAbDl+hx7oNjEqxZIrIC6AN8ISI3AtuPWk6AZ1R1SrknHRO95B5zo46pVc8CLsLxR/sO4MKjXl8PR6fP/vq/mbj8gHNVNe8E7ykPx95O2XzHanBWpKpHni8FCpzZSkWk7Oc8yLlOwzguc47B8GY/8b9vz4k45sE9YpCI+IlIM6ApsEVEmgLbVXUCjrbYHYFsILzM674ArndO7oKI1BeRhBOFcC4bqY6unXcDnY/6uR2YCzykqn+W+dGXlJkkSkTKvc5pE9D8BK+JPlG2Y2iJh3r4G//f3v2qRBQFcRz/TloQRLYazPsAaxGLT2Ey2EzKbrIKNh/AlzCabBvE7LKIdqvJskEYwznhHv8gsoZ1+X7gpnvDtB9zhntmdRgM+i/WIuK584yBY+AwIqbAAXDS+f4JmFDWnR5l5hzYB2ZR1iMOKMvuX4DbuknuIjNvKGf/d/UY5oo2OL6yDlzXOibA6MP7HWAbOOsMoDdr/cM6RH4APl2ZnpmPwEYdQgOcA/1a7z1lBvEbe5TFVNK3vHZbWnIRMQJeM3OhfxkiokcJrt3MfPuT4rSS7Bik5XdJnRksaAs4NRT0EzsGSVLDjkGS1DAYJEkNg0GS1DAYJEkNg0GS1HgHlD09yCfC4UUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu_0 = X[Y==0].mean()\n",
    "mu_1 = X[Y==1].mean()\n",
    "\n",
    "sigma2_0 = np.cov(X[Y==0])\n",
    "sigma2_1 = np.cov(X[Y==1])\n",
    "\n",
    "from scipy.stats import norm\n",
    "P_X_given_y1 = norm.pdf(X,mu_0,np.sqrt(sigma2_0))\n",
    "P_X_given_y2 = norm.pdf(X,mu_1,np.sqrt(sigma2_1))\n",
    "\n",
    "plt.hist(X[Y==0],np.linspace(-3,3,10),density=True,histtype='step')\n",
    "plt.hist(X[Y==1],np.linspace(-3,3,10),density=True,histtype='step')\n",
    "plt.plot(X,P_X_given_y1)\n",
    "plt.plot(X,P_X_given_y2)\n",
    "\n",
    "plt.xlabel('Lobster size (cm)')\n",
    "plt.ylabel('P(survival|lobster size)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior probabilities of each class (i.e. $P(y=1)=\\pi_1$) can be found by just counting the fraction of the lobsters that survived (or the fraction that died)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_0 = sum(Y==0)/m\n",
    "pi_1 = sum(Y==1)/m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classification was as easy as asking which class had the higher probability, given the data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What's the probability of survival for a lobster that is half a standard deviation below the mean in size.\n",
    "X_new = -0.5\n",
    "\n",
    "posterior_0 = norm.pdf(X_new,mu_0,np.sqrt(sigma2_0))*pi_0\n",
    "posterior_1 = norm.pdf(X_new,mu_1,np.sqrt(sigma2_1))*pi_1\n",
    "psum = posterior_0 + posterior_1\n",
    "\n",
    "posterior_0/=psum\n",
    "posterior_1/=psum\n",
    "\n",
    "print(\"Probability of lobster death is: \",posterior_0)\n",
    "print(\"Probability of lobster survival is: \",posterior_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is all well and good, but sometimes it's not so easy to make an assumption about how the data are distributed!  Why should lobster size be normally distributed?  How do we even know?  This leads us to an alternative way to perform classification called *discriminative classification*, in which rather than model the features as an explicit statistical model, we'll come up with a function that takes as input the features, and outputs a probability for each class.  \n",
    "\n",
    "What function should we use for such an endeavor?  Let's begin by looking at the two-class naive Bayes model again, where we are interested in the survival probability.  By definition, \n",
    "$$\n",
    "P(y=1|\\mathbf{x}) = \\frac{P(\\mathbf{x}|y=1)P(y=1)}{P(\\mathbf{x}|y=1)P(y=1) + P(\\mathbf{x}|y=0)P(y=0)}.\n",
    "$$\n",
    "We can factorize this as follows:\n",
    "$$\n",
    "P(y=1|\\mathbf{x}) = \\frac{1}{1 + \\frac{P(\\mathbf{x}|y=0)P(y=0)}{P(\\mathbf{x}|y=1)P(y=1)}}.\n",
    "$$\n",
    "Then making the substitution\n",
    "$$\n",
    "a(\\mathbf{x}) = \\ln \\frac{P(\\mathbf{x}|y=1)P(y=1)}{P(\\mathbf{x}|y=0)P(y=0)},\n",
    "$$\n",
    "we have that \n",
    "$$\n",
    "P(y=1|\\mathbf{x}) = \\frac{1}{1 + \\mathrm{e}^{-a}}.\n",
    "$$\n",
    "The quantity $a$ is the (logarithm of the) ratio of probabilities, also known as the odds.  Thus $a$ is called the log-odds.  This is a useful quantity because it squishes the odds down to a symmetric function on the complete real line, rather than just a positive value.  \n",
    "\n",
    "This function on the right hand side ends up being common enough that it has a name: the sigmoid (Greek for s-like)\n",
    "$$\n",
    "\\sigma(a) = \\frac{1}{1 + \\mathrm{e}^{-a}}.\n",
    "$$\n",
    "The sigmoid function takes a number on the real line, and squashes it down to a value between zero and one, a proper probability.  Here's what it looks like\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl81NW9//HXJ3tIQgIkQAhh3wWVEllqF21dcLerYLHaqtxabbW3am2xXn/a+9OrbW2t9ipa94WruNGKVVux5aeCBGSHQAgQkgDZQ/ZkknP/SPSXYjADTPKd5f18PHiQmTkm75Hk/Tg5851zzDmHiIiElyivA4iISOCp3EVEwpDKXUQkDKncRUTCkMpdRCQMqdxFRMKQyl1EJAyp3EVEwpDKXUQkDMV49YXT09PdqFGjvPryIiIhae3ateXOuYyexnlW7qNGjSI3N9erLy8iEpLMbK8/47QsIyIShlTuIiJhSOUuIhKGeix3M3vMzErNbPMRHjczu9/M8s1so5l9LvAxRUTkaPgzc38CmPsZj58DjO/8sxD47+OPJSIix6PHcnfO/ROo/IwhFwFPuQ6rgDQzywxUQBEROXqBWHPPAvZ1uV3UeZ+IiHgkENe5Wzf3dXt2n5ktpGPphhEjRgTgS4uIBAfnHI2tbdQ1+TjU5KOu2Udd59/1zT4aWttoaPZR39LGVycN5qTstF7NE4hyLwKyu9weDpR0N9A5txhYDJCTk6PDW0UkKLW1OyrqmymvbaGyvoWK+mYq61uoamilqr6FqoYWahpbP/lT2+TjUGMrvnb/am1wSnxIlPsy4DozWwLMAmqcc/sD8HlFRAKuocVHcVUjxdWN7K9pYn9NEwdqGimtbebgoWbKapuoqG/BddPTZpCaGMuAfnGf/D1qUBKpibGkJMSQkhBL/8QYkuNjSEmIITk+ln5x0STHx5AUH0O/uGgSY6OJiupuwSOweix3M3seOA1IN7Mi4D+AWADn3EPAcuBcIB9oAL7XW2FFRHrinKOivoXd5fXsLq9nT3k9hZUNFFY2sK+ygaqG1n8ZbwYZyfEM6Z9AVloCJ2enkpEcT0ZKPIOS4xmUFMeg5DgG9IsjrV8c0X1QzIHQY7k75+b38LgDrg1YIhERP1XWt7B9/yG2H6hlx8FadpbWkV9aR03j/y/wmCgja0AiIwb2Y+q0TIYPSCQrreNPZloig1PiiY0Ov/dzerZxmIjI0Siva2Z9YTWbimvYUlLD5uJDHDjU9MnjA/rFMn5ICuefmMnYjGRGZyQxJj2JrLREYsKwvHuicheRoNPe7thZWseHeyrJ3VPJusIq9lU2Ah3LKGMzkpkzdhBTMvszKTOFSUP7k5ES73Hq4KJyFxHPOefYXV7Pe7sqeD+/nA8KKqjuXBsfnBLPjJEDWDBrJCdnpzE1K5WkeFVXT/R/SEQ80dTaxvu7ylmxvYx3d5R+MjPPSkvkzMlDmDVmEDNHDSR7YCJmofEiZjBRuYtIn6ltauXv20p5c8sB/rGjjIaWNhJjozl13CAWfmksXxyXzshB/VTmAaByF5Fe1dTaxttbD/LnDSW8u6OMFl87g1Pi+dr0LM6cMoQ5YwcRHxPtdcywo3IXkYBzzrGusJqla4v4y8YSapt8DE6J59KZI7jgpEymZw/okzfyRDKVu4gEzKGmVl5ZV8xzqwvJO1hLYmw050wdyjdmDGf2mEEh8wagcKByF5Hjtru8nsff282LuUU0trYxLSuVu74+jQtOGkayrmzxhP6vi8gxW7u3iv9+dxd/336Q2KgoLjhpGJd/fiQnDu/dTbGkZyp3ETkqzjk+2FXBH97J54OCCgb0i+VHp49jwZyRDE5J8DqedFK5i4jf1u6t5J6/5rF6dyWDU+K59bzJXDprBP3iVCXBRv8iItKjHQdrufuN7byzvZT05Hhuv2AK82aOICFWlzAGK5W7iBxRZX0L9729g+c+LCQpLpqb507kis+P0kw9BOhfSEQ+pa3d8cyqvfz6rTwaWtpYMGsEN5wxgQFJcV5HEz+p3EXkX2wsqmbRK5vZVFzDF8al8x8XTGH8kBSvY8lRUrmLCACNLW385q08HntvN4OS47l//nQuODFT+7yEKJW7iJC7p5Kblm5kd3k935k1gp+dM4n+CbFex5LjoHIXiWAtvnZ++/YOHv7nLrLSEnnuqll8fly617EkAFTuIhGqoKyO65esZ1NxDfNnZnPreVN0CEYY0b+kSAR65aMiFr2ymbiYKB5aMIO5U4d6HUkCTOUuEkGafW3c8eetPLu6kJmjB/L7eSeTmZrodSzpBSp3kQhRUt3INc+sZUNRDf/2pTHcdPZEYqKjvI4lvUTlLhIB1u6t4t+eXktTaxsPLfgcc6dmeh1JepnKXSTMvbS2iJ+/vInMtASWLJzFuMF6Q1IkULmLhCnnHL9+K48HV+xizphB/PE7n9P2ARFE5S4Shlrb2vnZSxt5eV0x807J5s6LpxKr9fWIonIXCTN1zT6ueWYtK3eW8+9nTuBHXxmnLQQikMpdJIxUN7Rw+WMfsrnkEPd840S+fUq215HEIyp3kTBRVtvMZX9aTUF5PQ8vmMEZU4Z4HUk8pHIXCQMl1Y0seHQ1+2uaePyKUzhV+8NEPL9eYTGzuWaWZ2b5ZnZLN4+PMLMVZvaRmW00s3MDH1VEulNS3ci8xasoq23m6StnqtgF8KPczSwaeBA4B5gCzDezKYcNuxV4wTk3HZgH/DHQQUXk0w7UNDH/kVVU1bfw9FWzyBk10OtIEiT8mbnPBPKdcwXOuRZgCXDRYWMc0L/z41SgJHARRaQ7Bw91FHtFXQtPXjmTk7PTvI4kQcSfNfcsYF+X20XArMPG3A68ZWY/ApKAMwKSTkS6VVnfwnceXU3poSaeunImnxsxwOtIEmT8mbl3d4GsO+z2fOAJ59xw4FzgaTP71Oc2s4VmlmtmuWVlZUefVkSoa/ZxxeMfUljZwKOXn8KMkVqKkU/zp9yLgK4Xyw7n08suVwIvADjnPgASgE+9quOcW+ycy3HO5WRkZBxbYpEI1tTaxsKnctlScog/Xvo55owd5HUkCVL+lPsaYLyZjTazODpeMF122JhC4KsAZjaZjnLX1FwkgNraHTcsWc/7uyr49bdO1HXs8pl6LHfnnA+4DngT2EbHVTFbzOwOM7uwc9hPgavNbAPwPHCFc+7wpRsROQ6/en0rf91ygF+eP4WvTR/udRwJcn69ick5txxYfth9t3X5eCtwamCjicjHHl1ZwOPv7eHKL4zmyi+M9jqOhABtEycS5N7YtJ//XL6Nc6YOZdG5k72OIyFC5S4SxDYWVXPD/6xnenYa911yMlFR2t1R/KNyFwlSB2qauPqpXNKT41n83RwSYqO9jiQhRBuHiQShxpY2Fj6dS22Tj5eu+TzpyfFeR5IQo3IXCTLOOW5+aSObimtYfFkOkzP79/wfiRxGyzIiQebRlbv584YSbjxrImfqWnY5Rip3kSDyfn45d73RcWXMD08b63UcCWEqd5EgUVzdyHXPf8SYjGTu/dZJOvdUjovKXSQINPva+OEza2n1tfPwZTNIjtfLYXJ89B0kEgTuWr6dDUU1PLRgBmMzkr2OI2FAM3cRjy3ftJ8n3t/D908dzdypQ72OI2FC5S7ioT3l9dy8dCMnZ6dxyzmTvI4jYUTlLuKRZl8b1z63jugo44FLpxMXox9HCRytuYt45J6/5rGl5BCPfDeH4QP6eR1HwoymCiIeWLG9lD/9v91cPmek3qgkvULlLtLHSg81ceOLG5g0NIWfawtf6SValhHpQ+3tjp++uIH6Fh9L5s/WTo/SazRzF+lDT7y/h5U7y7n1vCmMH5LidRwJYyp3kT6y82Atd/91O1+ZNJjvzBrhdRwJcyp3kT7Q4mvn+iXrSY6P4e5vTNO+MdLrtOYu0gd+97cdbN1/iIcvm8HglASv40gE0MxdpJet3VvFQ//YxbdzhnP2CdpeQPqGyl2kFzW2tHHTixvITE3kl+dP8TqORBAty4j0ol+/lUdBeT3PXjWLlIRYr+NIBNHMXaSXfLi7ksfe281ls0dy6rh0r+NIhFG5i/SChhYfNy3dQPaAftrtUTyhZRmRXvDrN3ewt6KBJQtnk6RTlcQDmrmLBNjavVU8/v5uFswewewxg7yOIxFK5S4SQE2tbdy8dAPDUhO55RxtCibe0e+LIgF0/993squsnie/P1OHXIunNHMXCZDNxTU8/M8CvjljOF+ekOF1HIlwfpW7mc01szwzyzezW44w5ttmttXMtpjZc4GNKRLcfG3t3PLyRgb0i+PW87QcI97r8fdGM4sGHgTOBIqANWa2zDm3tcuY8cDPgVOdc1VmNri3AosEo8ff28Pm4kM8cOl00vrFeR1HxK+Z+0wg3zlX4JxrAZYAFx025mrgQedcFYBzrjSwMUWCV2FFA795O48zJg/mvGmZXscRAfwr9yxgX5fbRZ33dTUBmGBm75nZKjObG6iAIsHMOceiVzcRExXFnRdP1Va+EjT8eTm/u+9W183nGQ+cBgwHVprZVOdc9b98IrOFwEKAESN0WIGEvlfXF7NyZzl3XHQCmamJXscR+YQ/M/ciILvL7eFASTdjXnPOtTrndgN5dJT9v3DOLXbO5TjncjIydDWBhLaq+hbu/Ms2Ts5OY8GskV7HEfkX/pT7GmC8mY02szhgHrDssDGvAqcDmFk6Hcs0BYEMKhJs7npjGzWNrdz19WlERWk5RoJLj+XunPMB1wFvAtuAF5xzW8zsDjO7sHPYm0CFmW0FVgA3Oecqeiu0iNdWFVTwQm4RV31xNJMz+3sdR+RTzLnDl8/7Rk5OjsvNzfXka4scj2ZfG+f+fiXNvnbe/smXSYyL9jqSRBAzW+ucy+lpnN4fLXKUFv+jgF1l9Tz+vVNU7BK0tP2AyFHYU17PH1bkc96JmZw+Ue/Vk+Clchfxk3OOX762mbjoKG7TeagS5FTuIn56fdN+Vu4s56dnTWBI/wSv44h8JpW7iB8ONbVyx5+3MjWrP5fN1jXtEvz0gqqIH3771g7K6pp55Ls5xERrTiTBT9+lIj3YXFzDUx/sYcGskZyUneZ1HBG/qNxFPkN7u2PRq5sZmBTHjWdP9DqOiN9U7iKf4fk1hWzYV82i8yaTmhjrdRwRv6ncRY6goq6Ze/6ax+wxA7n45MN3uRYJbip3kSO4643t1Df7uPMi7dMuoUflLtKNNXsqWbq2iKu+OIbxQ1K8jiNy1FTuIofxtbXzy1c3Myw1gR9/dZzXcUSOia5zFznME+/vYfuBWh6+bAb94vQjIqFJM3eRLg7UNHHf2zv4yqTBnDVliNdxRI6Zyl2kiztf34qv3XH7BSfoRVQJaSp3kU4rd5bx+sb9XHv6OEYM6ud1HJHjonIXoeN0pdte28Lo9CQWfmmM13FEjpteLRKh43Sl3eX1PPX9mSTE6nQlCX2auUvEK6xo4IEV+Zw3LZMvTcjwOo5IQKjcJaI55/iPZZuJiTJuPX+y13FEAkblLhHtzS0HWZFXxk/OnEBmaqLXcUQCRuUuEau+2ccdf97CpKEpXPH5UV7HEQkolbtErPvf2UlJTRO/uniqTleSsKPvaIlIOw7W8qeVu7kkJ5ucUQO9jiMScCp3iTjt7Y5Fr2wiOSGGn50zyes4Ir1C5S4RZ+naItbsqeIX50xmYFKc13FEeoXKXSJKZX0L//eNbZwyagDfnDHc6zgivUblLhHlruXbqGvy8Z9fm0ZUlDYGk/ClcpeIsbqgghfXFnH1l8YwQacrSZhTuUtEaPa18YtXNpGVlsiPvzLe6zgivc6vcjezuWaWZ2b5ZnbLZ4z7ppk5M8sJXESR47f4HwXsKqvnV1+bSmKcNgaT8NdjuZtZNPAgcA4wBZhvZlO6GZcC/BhYHeiQIsejoKyOP6zI5/wTMzl94mCv44j0CX9m7jOBfOdcgXOuBVgCXNTNuDuBe4CmAOYTOS7OOW59dTPxMVHcdv6n5iQiYcufcs8C9nW5XdR53yfMbDqQ7Zz7SwCziRy3l9cV8/6uCn42dxKD+yd4HUekz/hT7t1dL+Y+edAsCrgP+GmPn8hsoZnlmlluWVmZ/ylFjkF5XTN3vr6VGSMHcOnMEV7HEelT/pR7EZDd5fZwoKTL7RRgKvCume0BZgPLuntR1Tm32DmX45zLycjQoQjSu+78y1bqm33c/XVd0y6Rx59yXwOMN7PRZhYHzAOWffygc67GOZfunBvlnBsFrAIudM7l9kpiET+syCvltfUlXHv6OMbrmnaJQD2Wu3POB1wHvAlsA15wzm0xszvM7MLeDihytOqbfdz6ymbGDU7mmtPGeh1HxBN+HZDtnFsOLD/svtuOMPa0448lcuzufTOPkppGlv5gDvExuqZdIpPeoSphJXdPJU9+sIfvzh7JjJHap10il8pdwkZTaxs3L93IsNREbp6rfdolsvm1LCMSCn73t50UlNfzzJWzSIrXt7ZENs3cJSxsLKrmkZUFXJKTzRfGp3sdR8RzKncJec2+Nm58cQPpyXH84rzJXscRCQr63VVC3u/+tpMdB+t4/HunkJoY63UckaCgmbuEtHWFVTz8j11ckpOtHR9FulC5S8hqau1YjhnaP4Fbz9dyjEhXWpaRkHXPX/MoKOu4OiYlQcsxIl1p5i4h6b38ch57bzeXzR6pq2NEuqFyl5BT09DKjS9uYEx6Er84V8sxIt3RsoyEnNuWbaastpmXrvm8zkMVOQLN3CWkLNtQwmvrS/jxV8dzUnaa13FEgpbKXULGvsoGFr28iekj0vihtvIV+UwqdwkJvrZ2rl/yEQD3z5tOTLS+dUU+i9bcJST8/u87WVdYzf3zp5M9sJ/XcUSCnqY/EvRWFVTwwIp8vjVjOBeeNMzrOCIhQeUuQa28rpnrl3zE6EFJ3H7hCV7HEQkZWpaRoNXW7vjJ/6ynuqGVx6+YqT3aRY6CflokaD24Ip+VO8u5++vTmDKsv9dxREKKlmUkKL2fX859f9vB16Zncckp2V7HEQk5KncJOvtrGvnR8x8xJj2JX108FTPzOpJIyFG5S1Bp9rXxg2fW0exr5+HLcrTOLnKM9JMjQeX2ZVvYsK+ahxbMYNzgZK/jiIQszdwlaDz/YSHPf7iPa08fy9ypQ72OIxLSVO4SFFYXVHDba5v50oQM/v3MiV7HEQl5Knfx3L7KBq55dh3ZA/vxh/nTiY7SC6gix0vlLp6qbWrlyifX0Nbu+NPlp5CaqOPyRAJBL6iKZzp2elzPrrJ6nvr+TEanJ3kdSSRsaOYunnDOcduyLbyzvZT/c+EJnDpO56CKBJLKXTzxx3d38dzqQq45bSwLZo/0Oo5I2FG5S5979aNi7n0zjwtPGsZNZ+nKGJHe4Fe5m9lcM8szs3wzu6Wbx//dzLaa2UYz+7uZaSom3VqRV8qNL25g9piB3PutE4nSlTEivaLHcjezaOBB4BxgCjDfzKYcNuwjIMc5dyKwFLgn0EEl9OXuqeSaZ9YycWgKi7+bQ3xMtNeRRMKWPzP3mUC+c67AOdcCLAEu6jrAObfCOdfQeXMVMDywMSXUbS05xPeeWMOw1ESe/P5M+ifokkeR3uRPuWcB+7rcLuq870iuBN7o7gEzW2hmuWaWW1ZW5n9KCWk7D9by3cdWkxwfw1NXziQ9Od7rSCJhz59y725R1HU70GwBkAPc293jzrnFzrkc51xORkaG/yklZOWX1jH/kdWYGc9cNYvhA3S4tUhf8OdNTEVA19MShgMlhw8yszOARcCXnXPNgYknoaygrI5LH1kFwPNXz2JshnZ5FOkr/szc1wDjzWy0mcUB84BlXQeY2XTgYeBC51xp4GNKqMkvrWX+I6toa3c8d/Usxg1O8TqSSETpsdydcz7gOuBNYBvwgnNui5ndYWYXdg67F0gGXjSz9Wa27AifTiLA5uIavv3wKtra4bmrZzNhiIpdpK/5tbeMc245sPyw+27r8vEZAc4lIWrt3iquePxDUuJjePbq2dovRsQj2jhMAmZFXinXPruOwSnxPHv1bLLSEr2OJBKxtP2ABMQLufu46slcxmQk8cIP5qjYRTymmbscF+ccf3gnn9++vYMvjk/nvxfMIFmHWot4Tj+FcsyaWtu45aWNvLq+hK9/Lov/+saJxEbrl0GRYKByl2NSeqiJq59ey4Z91dx09kR+eNpYzLQJmEiwULnLUVtXWMUPn1lHTWMrDy2YwdypQ72OJCKHUbmL35xzPPXBXn71+laGpiaw9Jo5nDAs1etYItINlbv4pa7Zxy9e3sSyDSWcMXkwv/nWyaT2086OIsFK5S49+qiwiuuXrKeoqoGbzp7INV8eq0M2RIKcyl2OyNfWzkP/2MV9f9vJ0P4JvPBvc8gZNdDrWCLiB5W7dGvnwVpuXLqRDfuqueCkYfzq4qmkJmoZRiRUqNzlX7S2tfPoyt3c9/YOkuKj+cP86Zx/YqYucxQJMSp3+cTavZUsemUz2w/UMveEodx58VQyUnRqkkgoUrkLFXXN3PtmHkvW7CMzNYGHFszg7BOGaLYuEsJU7hGsxdfOk+/v4f53dtLQ0sbVXxzNDWdMIEl7w4iEPP0UR6D2dsfrm/bzm7fy2FPRwGkTM1h07mTG61ANkbChco8gzjne2V7Kr9/awbb9h5g4JIUnvncKp00c7HU0EQkwlXsEaG93vL3tIA+8k8+m4hpGDurH7y45mQtOGka03owkEpZU7mGs2dfGnzfs55F/FpB3sJaRg/px99en8Y0Zw7U1r0iYU7mHobLaZp7/sJCnPthLeV0zE4ek8Pt5J3PetExiVOoiEUHlHiacc3xQUMGzqwt5a8sBWtscp03M4KovjOHUcYN0WaNIhFG5h7h9lQ28tK6Il9cVU1jZQGpiLJfPGcX8WSMYm5HsdTwR8YjKPQSVHmri9U37eX3jfnL3VmEGc8YM4oYzxnPutEwSYqO9jigiHlO5h4g95fW8vfUgb209QO7eKpyDSUNTuOnsiVw8PYustESvI4pIEFG5B6mm1jbW7Knk3bwy3s0rZVdZPQCTM/vz46+M54KTMhk3WG86EpHuqdyDRLOvjc3FNawqqOT9XeXk7qmi2ddOXEwUs0YP5DuzRnLmlCFkD+zndVQRCQEqd4+U1jaxvrCaj/ZVs3ZvFRv2VdPsawc6ZucLZo/k1HGDmDMmncQ4raGLyNFRufcy5xzF1Y1s31/LlpJDbC6pYXNxDftrmgCIiTKmDOso81NGDeSUUQMYlKxtdkXk+KjcA8TX1k5xdSMF5fXsKq0jv7SOnaV17DhQS22zDwAzGJOexMzRA5mWlcr0EWmcMCxVV7eISMCp3P3knKOqoZXiqkaKqxsoqmpkX2UDeysbKKxoYF9VA61t7pPxg5LiGDs4mYumD2NyZn8mDe3PxKEpJGs7XRHpAxHfNL62dqoaWqmsb6G8rpmy2mbK65o5eKiJ0tqOvw/UNLG/pumTNfGPJcfHMGJgPyYMSeGsE4YyJj2JUelJjM1I0tKKiHjKr3I3s7nA74Fo4FHn3N2HPR4PPAXMACqAS5xzewIbtXvOOZpa26lr9lHf7KOu2Udtk4/aplZqm3wcamrlUKOP6sYWahpbqWlopaqhheqGViobOu5z7tOfNz4miiH9ExicEs/UrFTOOmEoQ/snMCwtkeEDEslKSyStX6ze1i8iQanHcjezaOBB4EygCFhjZsucc1u7DLsSqHLOjTOzecB/AZf0RuD/WVPIw/8soKG5jfoWHw0tbbS1d9POh0mOjyE1MZbUxFgGJMUyLC2RAf3iGJgUx6Dkjr/Tk+PJSIknPTme/gkxKm4RCVn+zNxnAvnOuQIAM1sCXAR0LfeLgNs7P14KPGBm5lx3c+LjMzApnimZ/UmKiyExLpqk+GiS4mNIjo8hKS6GlIQYkhNiSInvKPKUhI77tBuiiEQSf8o9C9jX5XYRMOtIY5xzPjOrAQYB5V0HmdlCYCHAiBEjjinwmVOGcOaUIcf034qIRAp/prPdrU0cPiP3ZwzOucXOuRznXE5GRoY/+URE5Bj4U+5FQHaX28OBkiONMbMYIBWoDERAERE5ev6U+xpgvJmNNrM4YB6w7LAxy4DLOz/+JvBOb6y3i4iIf3pcc+9cQ78OeJOOSyEfc85tMbM7gFzn3DLgT8DTZpZPx4x9Xm+GFhGRz+bXde7OueXA8sPuu63Lx03AtwIbTUREjpWuDxQRCUMqdxGRMKRyFxEJQ+bVRS1mVgbs9eSLH590DntzVgSItOccac8X9JxDyUjnXI9vFPKs3EOVmeU653K8ztGXIu05R9rzBT3ncKRlGRGRMKRyFxEJQyr3o7fY6wAeiLTnHGnPF/Scw47W3EVEwpBm7iIiYUjlfhzM7EYzc2aW7nWW3mRm95rZdjPbaGavmFma15l6i5nNNbM8M8s3s1u8ztPbzCzbzFaY2TYz22Jm13udqa+YWbSZfWRmf/E6S29QuR8jM8um4+jBQq+z9IG3ganOuROBHcDPPc7TK7ocKXkOMAWYb2ZTvE3V63zAT51zk4HZwLUR8Jw/dj2wzesQvUXlfuzuA26mm0NJwo1z7i3nnK/z5io69vQPR58cKemcawE+PlIybDnn9jvn1nV+XEtH2WV5m6r3mdlw4DzgUa+z9BaV+zEwswuBYufcBq+zeOD7wBteh+gl3R0pGfZF9zEzGwVMB1Z7m6RP/I6OyVm710F6i19b/kYiM/sbMLSbhxYBvwDO6ttEveuznq9z7rXOMYvo+DX+2b7M1of8Oi4yHJlZMvAScINz7pDXeXqTmZ0PlDrn1prZaV7n6S0q9yNwzp3R3f1mNg0YDWwwM+hYolhnZjOdcwf6MGJAHen5fszMLgfOB74axqds+XOkZNgxs1g6iv1Z59zLXufpA6cCF5rZuUAC0N/MnnHOLfA4V0DpOvfjZGZ7gBznXChuQOQXM5sL/Bb4snOuzOs8vaXz/N8dwFeBYjqOmLzUObfF02C9yDpmKE8Clc65G7zO09c6Z+43OufO9zpLoGnNXfzxAJACvG1m683sIa8D9YbOF40/PlJyG/BCOBd7p1OBy4CvdP7bru+c0UqI08xdRCQMaeYuIhKGVO7NET3wAAAAK0lEQVQiImFI5S4iEoZU7iIiYUjlLiIShlTuIiJhSOUuIhKGVO4iImHofwHDdBgKCn2g0wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = np.linspace(-5,5,101)\n",
    "\n",
    "def sigmoid(a):\n",
    "    return 1./(1+np.exp(-a))\n",
    "\n",
    "plt.plot(a,sigmoid(a))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us some insight into what naive Bayes is doing: we're taking the log-odds of the data for the two classes, and turning it into a class probability by running it through the logistic function.  \n",
    "\n",
    "As it turns out, in this univariate case if we use the naive Bayes model, we can write $a$ as a quadratic function of the features (dropping the bold on $x$, since it's just a scalar):\n",
    "$$\n",
    "a(x) = w_2 x^2 + w_1 x + w_0 ,\n",
    "$$\n",
    "where \n",
    "$$\n",
    "w_2 = \\frac{1}{2}\\left[\\frac{1}{\\sigma^2_0} - \\frac{1}{\\sigma^2_1}\\right]\n",
    "$$\n",
    "$$\n",
    "w_1 = \\frac{\\mu_1}{\\sigma^2_1} - \\frac{\\mu_0}{\\sigma^2_0}\n",
    "$$\n",
    "$$\n",
    "w_0 = -\\frac{1}{2}\\frac{\\mu_1^2}{\\sigma_1^2} + \\frac{1}{2}\\frac{\\mu_0^2}{\\sigma_0^2} + \\ln\\frac{\\pi_1}{\\pi_0}.\n",
    "$$\n",
    "It's important to recognize that by using these definitions, we haven't done anything different from normal naive Bayes.  We can show this by computing the probability using the normal naive Bayes thing of evaluating the probability for each class, then normalizing and comparing it to what we get using these formulas\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHOZJREFUeJzt3Xl4VPW9x/H3NyEJEAhbwhIWg8hiBCwYFmvd6lK0Lm29rVLFK2Bja623drnVq7X3tvY+tVZtXarFWhGKIG1t5VYkYt26ECCRNUAgoEAICQFCAiQh2+/+MQFDMpABZnJm+byeZx45M4fwfcbkzY8zZ86Ycw4REYkucV4PICIiwae4i4hEIcVdRCQKKe4iIlFIcRcRiUKKu4hIFFLcRUSikOIuIhKF2o27mf3OzPaY2foTPG5m9pSZFZnZWjMbH/wxRUTkVHQKYJ/ZwDPAnBM8fg0wvPk2CXiu+b8nlZqa6jIyMgIaUkREfPLz8/c659La26/duDvnPjCzjJPsciMwx/muY5BrZj3NbIBzbvfJvm5GRgZ5eXnt/fEiItKCmW0PZL9gHHMfCOxssV3cfJ+IiHgkGHE3P/f5vRqZmWWbWZ6Z5ZWXlwfhjxYREX+CEfdiYHCL7UFAib8dnXOznHNZzrmstLR2DxmJiMhpCkbcFwG3N581MxmobO94u4iIhFa7L6ia2XzgMiDVzIqBHwEJAM6554HFwLVAEVANTA/VsCIiEphAzpaZ2s7jDvhm0CYSEYlW8+bBgw/Cjh0wZAj89Kdw660h+aMCOc9dRETO1Lx5kJ0N1dW+7e3bfdsQksDr8gMiIiHmnKPhgQc+CftR1dW+lXwIaOUuIhJE9Y1NbC0/xIaSKt9tt+/24c5i/79hx46QzKG4i4icpqraejbtPsiGkspjEd9ceoi6xiYAkjrFMap/d64Z3Z/qful0K9vV9osMGRKS2RR3EZF2OOcoqaw9thrf2BzyHfs/OczSOzmR89JTmH5RBpnpKWQOSGFoajKd4puPftc8evwxd4CuXX0vqoaA4i4isaVXLzhw4JPtnj2houLYZn1jE0V7Dn1ySKX5v5U19QCYwdA+yYwZ1IObJwwmc0AKmekp9O2ehJm/N+w3O/qiaQedLWO+Mxk7XlZWltOFw0SkQ7UOO75rpdR2S+Hh333Aht1VbClrdVhlQMqxgGcOSGFU/+4kJ3m3LjazfOdcVnv7aeUuIjHDHTjQ5mJYBnQ+VMW7hXvITO/BZ4ankjkghfPSU8jo0+KwSoRR3EUkqh0+0sB7heUsKSjlqZPsl/fQVR02U0dQ3EUk6lQcruPtjWXkFJTxwZZy6hqa6JOceML9T3KkPGIp7iISFcqqanmroJQlBaXkbttPY5MjvUdnbp00hCnn9Scrozf2RM82x9wB34uqUUZxF5GItX3fYZasLyWnoJQPd/iifXZaMnddcjZTRvdnzMAex5/BUlHR7tky0UJxF5GI4ZyjsOwgS9aXsmR9KZtKDwIwemAK37t6BFNG9+ecvt1P/kWiMOT+KO4iEtaamhyriw+QU1BKzvpSPt5XjRlMOKs3P7wuk6sz+zG4d1evxww7iruIhJ2GxiZWfLSfJQW+Qy5lVUfoFGd8+pxUsi8ZxlWZ/UjrnuT1mGFNcRcRbyQmQn39sc2mhATeWb2TnIJSlm4s40B1PZ0T4rhsRF8+N7ofnx3Vjx5dEjwcOLIo7iLS8VqFHcDq6/nM6IHc9/AbXHluPz53Xn8uHZFGl8R4j4aMbIq7iHQ4V1/v952iSa6R/IeuIrFTZL4rNJzoGRSRDnP4SAPPvlt0wscNFPYg0cpdREKutr6R3+du57n3trLvcB13ez1QDFDcRSRk6hqaeDVvJ8+8s4WyqiN85pxUvnP1COyJhDbH3AFI0AumwaK4i0jQNTQ28dqqXfzq7S3sOlBD1lm9+OXN47hwWB/fDnV1bV9UTUjw3S9BobiLSNA0NTn+b20Jv3x7Cx/tPczYQT343y+N4ZLhqW0/yEIhDynFXUTOmHOOnIIynly6mcKyg4zs153fTLuAqzP7nfzTiSRkFHcROW3OOd7fXM7jb21m3a5Kzk5N5qmp47huzADi4hR1LynuInJalm3dx+NvFZK3vYJBvbrw2L+N5YvjBkbsJxdFG8VdRE5J/vYKnlhayD+L9tEvJYlHvjCar2QN1vnpYUZxF5GArN9VyRNLN/POpj30SU7kh9dlcuukIXRO0OUBwpHiLiIntbnsIE8u3cyb60vp0SWB/5wykn+/MIPkJOUjnOn/joj49fHew/zy7c28vqaE5MRO3HvFcO68eCgpnfVGo0iguIvEMj+nKRbvP8zTfyvijx8WkxBvZF9yNnddMozeJ/mAaQk/irtIrPITdgek907mz/+1mGmTz+Luy4fRt3vnjp9NzpjiLiLHHM39e9+/jPSeXTydRc6Mzl0SkeMYKOxRQHEXiUHOOZzXQ0hIBRR3M5tiZoVmVmRm9/t5fIiZvWtmq8xsrZldG/xRRSQYqusa+M7CNThQ4KNYu3E3s3jgWeAaIBOYamaZrXZ7CFjonBsH3AL8OtiDisiZK9pzkBuf+Sd/Wb2Lp5du9h93p+RHg0BeUJ0IFDnntgGY2QLgRmBDi30ckNL86x5ASTCHFJEz9/rqXTzw2jq6JMQzd8YkPjM8VSGPYoHEfSCws8V2MTCp1T7/DbxlZt8CkoErgzKdiJyx2vpGfvzXDbyyfAcTMnrx9NTx9O+h0xujXSBx93fdztZ/3U8FZjvnHjezC4G5ZjbaOdd03BcyywayAYYMGXI684rIKdixr5q7X8ln/a4q7rr0bL5/9UhdtTFGBBL3YmBwi+1BtD3sMhOYAuCcW2ZmnYFUYE/LnZxzs4BZAFlZWfr3oEgI5RSU8r0/rMGAF27P4qrMfl6PJB0okL/CVwLDzWyomSXie8F0Uat9dgBXAJjZuUBnoDyYg4pIYOobm/jpGxu4a24+Q1OTeePeixX2GNTuyt0512Bm9wA5QDzwO+dcgZn9GMhzzi0Cvgu8YGb34Ttkc4dzeqVGpKPtrqzhnldWkb+9gmmTz+Kh684lqZMuyRuLArr8gHNuMbC41X0Pt/j1BuCi4I4mIqfig83lfPvV1dTWN/LU1HHccH661yOJh3RtGZEI19jk+NXftvD0O1sY0bc7z946nnP6dvN6LPGY4i4SwfYeOsK3F6zmH0V7uWn8IB75wmi6JOowjCjuIhFrxUf7ueeVD6msqefnN43ly1mDMD+X8ZXYpLiLRJimJsesv2/jsZxCBvfqwuzpE8lMT2n/N0pMUdxFIkhldT3f/cNq3t64h2vH9OfRm8bSXR97J34o7iIRYm3xAe6e9yFlVbX86PpM7vh0hg7DyAkp7iJhzjnH73O385O/biS1WyIL77qQcUN6eT2WhDnFXSSMHTrSwAOvreP/1pRw2cg0nvzKp+ilD6qWACjuImGqsPQg35iXz8d7D/P9z43kG5cOIy5Oh2EkMIq7SDhodey8CbjxocV0S0pg3p2TuXBYH2/mkoiluIt4zc+LogZseORa9lbV0Le7rr0up05xFwlDR3OvsMvp0lX7RcKUjq7LmVDcRTyma2NLKCjuIh5as/MADgVegk9xF/HIqh0V3Pbb5Vz66N/8x12fdyNnQHEX8UD+9v1Me3EFvbsl8mr2hcQ554t5y5vIGVDcRTrYyo/3c/uLK0jrnsSC7Mmk9+zi9UgShXQqpEgHWr5tH9Nnr6R/SmfmZ0+mX4pOdZTQ0MpdpIMs27qPO15ayYAenVmgsEuIaeUu0gH+WbSXmS+vZHCvrrzytcmkdU/yeiSJclq5i4TYB5vLmTF7JRl9kpmfrbBLx9DKXSSE3ivcQ/bcfIaldWPenZPorcv1SgfRyl0kRN7ZVEb2nHyG9+3GKwq7dDDFXSQE3t5Qxl1z8xnZvzvz7pykD9iQDqfDMiJBllNQyj2vfEjmgBTmzJxEjy76AGvpeIq7SBC9uW4335q/ijGDevDyjImkdFbYxRs6LCMSJG+s3c0981dx/uCezFHYxWNauYsEwaI1Jdz36mrGD+nJS9Mn0i1JP1riLa3cRc7QX1bt4tsLVnHBWb2YrbBLmFDcRc7An/KLuW/haiYN7cPs6RNIVtglTOg7UeQ0LczbyQ/+tJaLhqXywu1ZdEmM93okkWMUd5HTsGDFDu5/bR0XD/eFvXOCwi7hRYdlRE7RvOXbuf+1dVw6Ik1hl7CllbvIKZiz7GMefr2Az47qy3O3jSepk8Iu4SmglbuZTTGzQjMrMrP7T7DPV8xsg5kVmNkrwR1TxHsv/fMjHn69gCvP7aewS9hrd+VuZvHAs8BVQDGw0swWOec2tNhnOPAAcJFzrsLM+oZqYBEv/Pbv23jkjY187rx+PD11PImddERTwlsg36ETgSLn3DbnXB2wALix1T5fA551zlUAOOf2BHdMEe/M+mArj7yxkWtG9+eZryrsEhkC+S4dCOxssV3cfF9LI4ARZvZPM8s1syn+vpCZZZtZnpnllZeXn97EIh3oufe28r+LN/H5sQN4auo4EuIVdokMgXynmp/7XKvtTsBw4DJgKvBbM+vZ5jc5N8s5l+Wcy0pLSzvVWUU61DPvbOHRJZu44fx0fnXzpxR2iSiBnC1TDAxusT0IKPGzT65zrh74yMwK8cV+ZVCmFOkI9sk6xgF3A1sXrOIXXz6f+Dh/axyR8BXIUmQlMNzMhppZInALsKjVPn8BLgcws1R8h2m2BXNQkZCy4+NtzbfHbxmnsEtEajfuzrkG4B4gB9gILHTOFZjZj83shubdcoB9ZrYBeBf4vnNuX6iGFukIht7lJ5HLnGt9+LxjZGVluby8PE/+bJE27CSrc49+RkT8MbN851xWe/tpYSJC2zMERCKd4i4x7+9bynEo8BJdFHeJaat2VHDX3HyuffJ9/3HXIRmJUIq7xKwtZQeZPnslqd2SmDNjInHO+WLe8iYSoRR3iUnFFdVMe3EFCfFx/H7mJPqmdPZ6JJGgUtwl5uw7dITbX1zB4boG5syYyJA+Xb0eSSToFHeJKQdr67njpZWUVNbwuzsmcO6AFK9HEgkJxV1iRm19I9lz8tm4u4rnbr2ACRm9vR5JJGT0SUwSExoam7h3/iqWbdvHL2/+FJeP0kcOSHTTyl2innOO//rzOt7aUMaPrs/kC+NaX7FaJPoo7hL1frZkEwvzirn3iuFMv2io1+OIdAjFXaLa8+9v5Tfvb2Pa5LO478rhXo8j0mEUd4lar67cwc/e3MT156fzPzech53s4mAiUUZxl6i0ZH0pD7y2jktGpPH4l88nTtdklxijuEvU+dfWvdw7fxXnD+7J87fpA60lNum7XqLK2uIDfO3lPDJSu/LSHRPomqizfSU2Ke4SNbaWH+KOl1bSKzmROTMm0bNrotcjiXhGcZeoUHKghmm/XU6cwdyZk+jfQxcCk9imuEvE23+4jmkvLudgbQOzp09kaGqy1yOJeE4HJCWiHTrSwPSXVlBcUcOcGRMZPbCH1yOJhAXFXSLWkYZGvj43n/UlVTx/2wVMOruP1yOJhA0dlpGI1NjkuO/V1fyjaC8/v2ksV2X283okkbCiuEvEcc7x0F/Ws3hdKQ99/lxuumCQ1yOJhB3FXSLOL94qZP6KHXzz8mHcefHZXo8jEpYUd4kov/37Np59dytTJw7he1eP9HockbCluEvE+GN+MY+8sZFrx/TnkS+M1oXARE5CcZeIsHRDGT/401o+c04qT978KeJ1ITCRk1LcJezlbtvHN1/5kNEDe/CbaReQ1Cne65FEwp7iLmFt/a5KvvZyHkN6+y4Elpykt2aIBEI/KRJ+mo+lO+A8YA1QdqCa3sm6EJhIoLRyl/DS4kVSa3Eb0LOrVxOJRCTFXcKeXjoVOXWKu4QV5/UAIlFCcZewsaeq1usRRKJGQHE3sylmVmhmRWZ2/0n2+zczc2aWFbwRJRbsqarllhdycWj1LhIM7cbdzOKBZ4FrgExgqpll+tmvO3AvsDzYQ0p0Oxr20spa8rbt83+M3Sn5IqcikJX7RKDIObfNOVcHLABu9LPfT4CfA/q3tQSsZdhfnjGRiUN7+0Le+iYipySQuA8EdrbYLm6+7xgzGwcMds79NYizSZTbU1XLLbNyKWsO+4SM3l6PJBI1Aom7v38lH1tKmVkc8CTw3Xa/kFm2meWZWV55eXngU0rUKTsa9qpaZivsIkEXSNyLgcEttgcBJS22uwOjgffM7GNgMrDI34uqzrlZzrks51xWWlra6U8tEa2sqpapzWHXil0kNAKJ+0pguJkNNbNE4BZg0dEHnXOVzrlU51yGcy4DyAVucM7lhWRiiWgtV+wvz5hIlsIuEhLtxt051wDcA+QAG4GFzrkCM/uxmd0Q6gElepRW+sK+R2EXCbmALhzmnFsMLG5138Mn2PeyMx9Lok1pZS1TX8il/OAR5sycyAVnKewioaR3qErItQz7yzMmKOwiHUCX/JWQ8h2KWcbeQ3W8PGMiF5zVy+uRRGKCVu4SMrsraxR2EY8o7hISuytrmDorV2EX8YgOy0jQ+Vbsuew7VMecmRMZP0RhF+loWrlLUJUc8IV9v8Iu4imt3CVoSg7UMPWFT8I+TmEX8YziLkFxdMVecVhhFwkHOiwjZ0xhFwk/WrnLGdl1wHdWTMXhOubeOYlPDe7p9UgiglbucgaOhb1aYRcJN1q5y2nZdcD3BqUD1fXMnamwi4QbxV1OWcuw/37mJM5X2EXCjuIup6S4opqpL+Qq7CJhTnGXgBVXVHPLrFwqaxR2kXCnuEtAjoa9qqaeeXdOYuwghV0knCnucmLm+2x0BwwEPgDW76xQ2EUigE6FFP+aww5gLW5jB+sNSiKRQHGXgFn7u4hImFDcxS/n9QAickYUdzlOXUMTP3p9vddjiMgZUtzlmN2VNdw8axkvL9uOQ6t3kUimuAsA/yray3VP/YPNpQd59qvjiXPO/zF2p+SLRAKdChnjnHM8//42HsvZxNlp3Xj+tgs4p2+3ow96O5yInDbFPYZV1dbzvYVreGtDGZ8fO4Cf3zSW5CR9S4hEA/0kx6hNpVV8fW4+xRU1/PC6TGZclIGZTnYUiRaKewz6y6pd3P/aWlI6JzA/ezITMnp7PZKIBJniHkPqGpp45I0NzFm2nYlDe/PMV8fRt3tnr8cSkRBQ3GPE7soa7p73Iat2HOBrFw/lP6eMIiFeJ0uJRCvFPQb8q2gv35q/itr6Rn5963iuHTPA65FEJMQU9yjmnOO597fyi5zCtqc5ikhUU9yjVFVtPd9duIalG8q4buwAHtVpjiIxRT/tUajlaY4PX5fJdJ3mKBJzFPco8+dVxTzw2jqd5igS4xT3KFHX0MRP/rqBubnbmTS0N0/rNEeRmBbQuXBmNsXMCs2syMzu9/P4d8xsg5mtNbO/mdlZwR9VTmR3ZQ1f+c0y5uZuJ/uSs5l35ySFXSTGtbtyN7N44FngKqAYWGlmi5xzG1rstgrIcs5Vm9k3gJ8DN4diYDne0dMcjzQ08dyt47lGpzmKCIGt3CcCRc65bc65OmABcGPLHZxz7zrnqps3c4FBwR1TWmtqcvz6vSJue3E5vZMTef2eixR2ETkmkGPuA4GdLbaLgUkn2X8m8OaZDCUnV1njO83x7Y1lXH9+Oj/70hid5igixwmkCP7OofN7oW8zuw3IAi49wePZQDbAkCFDAhxRWtq4u4pv/N53muOPrs/kjk/rNEcRaSuQwzLFwOAW24OAktY7mdmVwIPADc65I/6+kHNulnMuyzmXlZaWdjrzxh6zYzdnxsj0HtTUN7IgezLTLxqqsIuIX4HEfSUw3MyGmlkicAuwqOUOZjYO+A2+sO8J/pgxqlW4rfm27MGryNL56yJyEu3G3TnXANwD5AAbgYXOuQIz+7GZ3dC822NAN+APZrbazBad4MvJKfB37MvQB9+KSPsCehXOObcYWNzqvodb/PrKIM8V05xz/H3LXi72ehARiVg6xSLMrN55gEff3MSybfv4yOthRCRiKe5homjPQX6Rs5klBaX0SU7kv6/PxD3q/1QlEZH2KO4eKzlQwy/f3swf84vpmtiJ+64cwcyLh9ItqRM41+ZFVcB3v4jISSjuHtl/uI5fv1vEnNzt4GD6RUO5+7Jh9OmWdPyOCrmInAbFvYMdPtLAi//4iFkfbKO6roGbxg/i21eNYGDPLl6PJiJRRHHvIHUNTcxfsYOn39nC3kN1fO68fnzv6pEM79fd69FEJAop7iHW2ORYtGYXj7+1meKKGiaf3ZtZt49i/JBeXo8mIlFMcQ8R5xzvbNrDYzmFbCo9yHnpKfz0i2O4ZHiqLhkgIiGnuIfAyo/38+ibm8jbXkFGn648PXUcnx8zgLg4RV1EOobiHkQbd1fxWE4h72zaQ9/uSfz0i6P5StZgEuJ1wQAR6ViKexDs2FfNE0sLeX1NCd2TOvGDKaO449MZdEmM93o0EYlRivsZKD94hKff2cL8FTuIjzO+fukwvn7JMHp0TfB6NBGJcYp7oFq8COqAJuDSH77JkYYmbp4wmP+4Yjj9UvSh1CISHhT3QPi5rnocsO4n17Cj/BBDU5M9GUtE5ET0Sl87GpvcCa+rHg8Ku4iEJa3c/Whscqz8eD+L1+3mzfWlrPB6IBGRU6S4N2tscqz4yBf0JQWllB88QlKnOC4f2dfr0URETllMx72xybH8o32+oK8vY++hI3RO8AX92jED+OyoviQndYLbvZ5UROTUxFzcGxqbWPHRft5Yt5ucglL2Hqqjc0Icnx3lC/rlI5uD3pKuqy4iESYm4t7Q2MTyo0FfX8q+w3V0SYj/JOij0uia2M5ToZCLSASJ2rg3NDaRu80X9LcKfEHvmhh/3Apd7yAVkWgVVXFvaGxi2TbfMfScgjL2Nwf9inP78fkx/bl0hIIuIrEhsuI+cCCUlHyynZ5O/Y6dLNt6NOilVFTXk9wc9GvHDOCykWl0TlDQRSS2RE7cW4cdcCUl7O2Zxu3fnENyYjxXZvqCfukIBV1EYlvkxL1V2MH3LtH+h/Yza9oFXKKgi4gcEzlxPwEDrj6vv9djiIiEFV1bRkQkCkVO3NPTT+1+EZEYFjlx37WrbcjT0333i4jIcSLrmLtCLiISkMhZuYuISMAUdxGRKKS4i4hEIcVdRCQKKe4iIlHInEfXKTezcmD7KfyWVGBviMaJVHpO2tJz0paek7Yi+Tk5yzmX1t5OnsX9VJlZnnMuy+s5womek7b0nLSl56StWHhOdFhGRCQKKe4iIlEokuI+y+sBwpCek7b0nLSl56StqH9OIuaYu4iIBC6SVu4iIhKgiIq7mT1mZpvMbK2Z/dnMeno9k9fM7MtmVmBmTWYW1a/+t8fMpphZoZkVmdn9Xs/jNTP7nZntMbP1Xs8SLsxssJm9a2Ybm39u/sPrmUIlouIOLAVGO+fGApuBBzyeJxysB74EfOD1IF4ys3jgWeAaIBOYamaZ3k7ludnAFK+HCDMNwHedc+cCk4FvRuv3SUTF3Tn3lnOuoXkzFxjk5TzhwDm30TlX6PUcYWAiUOSc2+acqwMWADd6PJOnnHMfAPu9niOcOOd2O+c+bP71QWAjMNDbqUIjouLeygzgTa+HkLAxENjZYruYKP2hleAwswxgHLDc20lCI+w+rMPM3gb8feL1g86515v3eRDfP6/mdeRsXgnkORHMz306FUz8MrNuwJ+AbzvnqryeJxTCLu7OuStP9riZ/TtwHXCFi5HzONt7TgTwrdQHt9geBJR4NIuEMTNLwBf2ec6517yeJ1Qi6rCMmU0BfgDc4Jyr9noeCSsrgeFmNtTMEoFbgEUezyRhxswMeBHY6Jx7wut5Qimi4g48A3QHlprZajN73uuBvGZmXzSzYuBC4A0zy/F6Ji80v9B+D5CD70Wyhc65Am+n8paZzQeWASPNrNjMZno9Uxi4CJgGfLa5IavN7FqvhwoFvUNVRCQKRdrKXUREAqC4i4hEIcVdRCQKKe4iIlFIcRcRiUKKu4hIFFLcRUSikOIuIhKF/h/Lmk565lSmTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w_2 = 0.5*(1./sigma2_0 - 1./sigma2_1)\n",
    "w_1 = mu_1/sigma2_1 - mu_0/sigma2_0\n",
    "w_0 = -0.5*mu_1**2/sigma2_1 + 0.5*mu_0**2/sigma2_0 + np.log(pi_1/pi_0)\n",
    "\n",
    "\n",
    "w_nb = np.array([w_0,w_1,w_2])\n",
    "\n",
    "posterior_0 = norm.pdf(X,mu_0,np.sqrt(sigma2_0))*pi_0\n",
    "posterior_1 = norm.pdf(X,mu_1,np.sqrt(sigma2_1))*pi_1\n",
    "psum = posterior_0 + posterior_1\n",
    "\n",
    "posterior_0/=psum\n",
    "posterior_1/=psum\n",
    "\n",
    "a = w_0 + w_1*X + w_2*X**2\n",
    "posterior_logistic = sigmoid(a)\n",
    "\n",
    "plt.plot(X,posterior_1)\n",
    "plt.plot(X,posterior_logistic,'ro')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They're the same, we haven't changed anything by writing the classifier in this way.  \n",
    "\n",
    "The change to logistic regression from naive Bayes comes when we *don't* use these definitions for the parameters defining the log odds $a(x)$.  Instead, what we want to do is find the values of the parameters $\\mathbf{w}$ that maximize the accuracy of the model or (equivalently) maximize the probability of a correct classification.  How do we define a function for maximizing the probability of a correct classification?  Well it should be a distribution on the boolean-valued class $y$ given the data $x$.  The Bernoulli distribution is usually appropriate for this sort of thing:\n",
    "$$\n",
    "P(y|x) = p(x)^y (1-p(x))^{1-y}.\n",
    "$$\n",
    "What should our mean $p(x)$ be?  It should be our prediction:\n",
    "$$\n",
    "p(x) = \\sigma(w_0 + w_1 x + w_2 x^2).\n",
    "$$\n",
    "We can write this more generally as\n",
    "$$\n",
    "\\theta = \\sigma(\\mathbf{w}^T \\boldsymbol{\\phi}(x)),\n",
    "$$\n",
    "where $\\boldsymbol{\\phi}(x) = [1,x,x^2]$, which is basically the expansion of the features as monomials, just like we saw in linear regression.  If we consider all of our data points, we can make a matrix\n",
    "$$\n",
    "\\Phi(X) = \\left[\\begin{array}{ccc} 1 & x_1 & x_1^2 \\\\\n",
    "                      1 & x_2 & x_2^2 \\\\\n",
    "                      \\vdots & \\vdots & \\vdots \\\\\n",
    "                      1 & x_m & x_m^2 \\end{array}\\right]\n",
    "                     $$\n",
    "                           \n",
    "If we assume that all our data points are independent, we can generate a likelihood model that we can maximize:\n",
    "$$\n",
    "P(Y|X) = \\prod_{i=1}^m P(Y_i|X_i),\n",
    "$$\n",
    "where $Y$ is a vector containing all of our class labels.  Substituion of the definition of the Bernoulli in, we get\n",
    "$$\n",
    "P(Y|X) = \\prod_{i=1}^m \\theta(X_i)^{Y_i} (1-\\theta(X_i))^{1-Y_i},\n",
    "$$\n",
    "where \n",
    "$$\n",
    "\\theta(X_i) = \\mathbf{w}^T \\phi_i.\n",
    "$$\n",
    "\n",
    "To find the best values of $\\mathbf{w}$, we can maximize this function with respect to it.  Per usual, it can be tricky to find the gradient of this function that has a product in it.  Instead, we'll find the minimum of the negative log-likelihood:\n",
    "$$\n",
    "\\mathcal{L} = -\\sum_{i=1}^m Y_i \\ln \\theta_i + (1-Y_i) \\ln (1-\\theta_i).\n",
    "$$\n",
    "This log-likelihood shows up fairly frequently, and has it's own name: cross-entropy.  You can think of it as least squares for binary variables, in the sense that it gets big when we mis-classify our data, and it becomes small when we mostly get stuff right.  To minimize it, we'll take the derivative, and try to find the place where it's zero.  The derivative of this function is pretty easy to find, and it turns out to be\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial w_j} = \\sum_{i=1}^m [\\theta_i - Y_i] \\Phi_{ij}.\n",
    "$$\n",
    "\n",
    "Unfortunately, unlike in the case of linear regression, we can't simply set these derivatives to zero and solve a linear system of equations to get the values for $w_i$.  Instead, we'll have to search for the optimal parameter values.  How should we do this?  First let's take stock of the information that we have: first, we have a formula for $\\theta$,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Phi = np.vander(X,N=3,increasing=True)\n",
    "\n",
    "def theta(w,Phi):\n",
    "    a = np.dot(Phi,w)\n",
    "    return sigmoid(a)\n",
    "\n",
    "# The objective function\n",
    "def L(w,Phi,y):\n",
    "    likelihood = -sum(y*np.log(theta(w,Phi)) + (1-y)*np.log(1-theta(w,Phi)))\n",
    "    return likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have the derivative of this function with respect to the values of the weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The gradient of the objective function\n",
    "def dLdw(w,Phi,y):\n",
    "    return np.dot(theta(w,Phi)-y,Phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is it possible to just follow this gradient downhill until we can't go downhill anymore?  Absolutely!\n",
    "\n",
    "This technique is called gradient descent, and the formula is as follows:\n",
    "$$\n",
    "\\mathbf{w}_{t+1} = \\mathbf{w}_t - \\eta \\nabla_{\\mathbf{w}} \\mathcal{L}.\n",
    "$$\n",
    "Note this parameter $\\eta$; this is called the learning rate, and we'll need to adjust it to get good performance.  If it is too small, we won't adjust the parameter very fast and it may take a very long time to reach the minimum.  If it's too large we may skip over the optimal solution without ever seeing it.  There are many variants of this algorithm (we'll explore some later in the course), but for now let's just implement the simplest possible version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.random.randn(3)*1.0\n",
    "eta = 5e-3\n",
    "\n",
    "cross_entropy = [L(w,Phi,y)]\n",
    "\n",
    "g_initial = np.sqrt(sum(dLdw(w,Phi,y)**2))\n",
    "rtol=1e-3\n",
    "converged=False\n",
    "while converged==False:\n",
    "    g = dLdw(w,Phi,y)\n",
    "    w -= eta*g\n",
    "\n",
    "    gnorm = np.sqrt(sum(g**2))\n",
    "    if gnorm/g_initial<rtol:\n",
    "        converged=True\n",
    "        \n",
    "    cross_entropy.append(L(w,Phi,y))\n",
    "\n",
    "plt.plot(cross_entropy)\n",
    "plt.axhline(L(w_nb,Phi,y),color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, the naive Bayes solution is not optimal for this likelihood function (although it's close).  How about the actual training accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sum((theta(w,Phi)>0.5) == y)/len(y))\n",
    "print(sum((posterior_1>0.5) == y)/len(y))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Turns out the same.  This isn't too surprising, since there isn't much flexibility in only 1D.  However, next time, we'll develop the capacity to use this model in a multi-class setting, and see how it can substantially out-perform naive Bayes."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}