aula16-LDA-naiveBayes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aula 16 &mdash; Análise de discriminante e Naïve Bayes\n",
    "\n",
    "Renato Vimieiro\n",
    "\n",
    "rv2 {em} cin.ufpe.br\n",
    "\n",
    "maio 2017"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introdução\n",
    "\n",
    "Na aula passada vimos o algoritmo de Regressão Logística que é usado para o problema de classificação. O uso de regressão logística, contudo, é recomendado apenas quando o número de classes é dois. Isso porque o modelo apresenta problemas de estabilidade quando o número de amostras é pequeno ou quando as classes são bem separadas (ver ISLR). Nesses casos, outros classificadores são mais recomendados. Nessa aula veremos dois classificadores que se assemelham ao modelo de regressão logística, mas que estimam a probabilidade condicional das classes dados os atributos de maneira diferente.\n",
    "\n",
    "Como dito, o objetivo central de regressão logística era ajustar um modelo para predizer $P(Y=1~|~X)$. Pelo teorema de Bayes, sabemos que tal probabilidade condicional pode ser reescrita como\n",
    "\n",
    "$$\n",
    "P(Y=k~|~X) = \\frac{P(X~|~Y=k)P(Y=k)}{P(X)}\n",
    "$$.\n",
    "\n",
    "Dessa forma, podemos enxergar o problema como a estimativa da probabilidade de se observar $X$ na classe $k$, a probabilidade de uma amostra ser da classe $k$ e a probabilidade de se observar o valor $X$. Em outras palavras, invertemos nossa análise para a probabilidade de se observar os valores nas diferentes classes. Assim, classificamos uma amostra como da classe $k$ quando seus atributos tiverem maior probabilidade de ocorrer nesta classe. Essa é a intuição por trás da Análise de discriminante e Naive Bayes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Análise de discriminante\n",
    "\n",
    "A análise de discriminante supõe que os atributos sigam uma distribuição normal em cada uma das classes. No caso mais simples, linear, é assumido ainda que a variância dos atributos é a mesma em todas as classes. Portanto, o problema se resume a estimar as funções de densidade $f_k(X)$.\n",
    "\n",
    "Se expressarmos $P(Y=k) = p_k$, podemos reescrever a equação acima como\n",
    "\n",
    "$$\n",
    "P(Y=k~|~X) = \\frac{p_k~f_k(X)}{\\sum p_i~f_i(X)}\n",
    "$$\n",
    "\n",
    "Iniciamos nossa análise pelo caso mais simples onde o modelo contém apenas um atributo. Como a suposição é de que os atributos sigam uma distribuição normal em cada classe com uma mesma variância $\\sigma^2$, temos que\n",
    "\n",
    "$$\n",
    "f_k(X) = \\frac{1}{\\sigma\\sqrt{2\\pi}}\\exp\\left({-\\frac{1}{2\\sigma^2}(X-\\bar{X_k})^2}\\right)\n",
    "$$\n",
    "\n",
    "e consequentemente\n",
    "\n",
    "$$\n",
    "P(Y=k~|~X) = \\frac{p_k~\\frac{1}{\\sigma\\sqrt{2\\pi}}\\exp\\left({-\\frac{1}{2\\sigma^2}(X-\\bar{X_k})^2}\\right)}{\\sum p_i~\\frac{1}{\\sigma\\sqrt{2\\pi}}\\exp\\left({-\\frac{1}{2\\sigma^2}(X-\\bar{X_i})^2}\\right)}\n",
    "$$\n",
    "\n",
    "A classe da amostra $X$ é aquela com a maior probabilidade. Com alguma manipulação algébrica, pode-se verificar que o máximo da equação acima é equivalente a\n",
    "\n",
    "$$\n",
    "argmax_k(P(Y=k~|~X)) = argmax_k\\log(P(Y=k~|~X)) = argmax_k \\delta_k(X) = X\\frac{\\bar{X_k}}{\\sigma^2}-\\frac{\\bar{X_k}^2}{2\\sigma^2} + log(p_k)\n",
    "$$\n",
    "\n",
    "Assim, precisamos estimar $p_k = n_i/n$, a média de $X$ e sua variância.\n",
    "\n",
    "Vamos retomar o exemplo da aula passada e ajustar um modelo para os dados."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "prostate = pd.read_csv(\"https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data\",\n",
    "                       header=0,index_col=0,sep='\\t')\n",
    "prostate['prognostico'] = prostate.gleason.map(lambda x: x >=7 and 'mau' or 'bom')\n",
    "treino = prostate[prostate.train=='T'].drop('train',axis=1)\n",
    "teste = prostate[prostate.train=='F'].drop('train',axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "lpk = np.log(treino.prognostico.value_counts()/treino.shape[0])\n",
    "avgPSA = treino.groupby('prognostico').lpsa.mean()\n",
    "varPSA = treino.lpsa.var()\n",
    "deltaK = [treino.lpsa*(avgPSA[i]/varPSA)-((avgPSA[i]**2)/(2*varPSA))+lpk[i] for i in ['bom','mau']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1106225c0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XlclXXa+PHPxY4KuIALm6CCgqKiuJdmLmmZlm22jS3TMlNPMzXzzNRMTz3TTDVTv5rmmZZps3XM1DYrzVyyrNxQUVlEUVFAEERFEFnP9/fHfZiIVI5wDvd9zvm+Xy9e55z73MsFysV9vsv1FaUUmqZpmnfwMTsATdM0rePopK9pmuZFdNLXNE3zIjrpa5qmeRGd9DVN07yITvqapmleRCd9TdM0L6KTvqZpmhfRSV/TNM2L+JkdQEvh4eEqLi7O7DA0TdPcytatW48qpSJa28+hpC8iM4B/AL7Aa0qpv55lv6uApcAopVS6fdtDwO1AI3CfUmrlua4VFxdHenq6I2FpmqZpdiJy0JH9Wk36IuILvABMAwqBLSKyTCmV3WK/EOBXwKZm25KBecBgIBJYLSKJSqlGR78RTdM0zXkcadMfDeQppfYrpeqARcCcM+z3Z+BvQE2zbXOARUqpWqXUASDPfj5N0zTNBI4k/SigoNnrQvu2/xCREUCMUurz8z3WfvydIpIuIullZWUOBa5pmqadv3Z35IqID/AscEtbz6GUegV4BSAtLU3XetY0rU3q6+spLCykpqam9Z3dVFBQENHR0fj7+7fpeEeSfhEQ0+x1tH1bkxBgCLBORAB6A8tEZLYDx2qapjlNYWEhISEhxMXFYc9HHkUpRXl5OYWFhcTHx7fpHI4072wBEkQkXkQCMDpmlzULokIpFa6UilNKxQEbgdn20TvLgHkiEigi8UACsLlNkWqaprWipqaGHj16eGTCBxARevTo0a5PMq3e6SulGkTkXmAlxpDNBUqpLBF5DEhXSi07x7FZIrIYyAYagHv0yB1N01zJUxN+k/Z+fw616SullgPLW2x75Cz7XtTi9ePA422MT9O8g80GPnqCvOZ6lpuRq2leofoYbH4VDm+D4p1QfRT6T4Ehc2HgTAgMMTtCrY3y8/OZNWsWmZmZZodyRjrpa1pHy1sDH/8STpVCeCLEXQBBobB7OexZAcHdYd5C6DvO7Eg1D6Q/T2paR2lsgBUPwrtzIbgr3LkO7tkEV70Klz0D92fBrSugU3d4ew5kf2J2xFobNTQ0cOONN5KUlMTVV19NdXU1a9asITU1lZSUFG677TZqa2sBo/TMQw89xPDhw0lLS2Pbtm1ccskl9O/fn3/9619Oj03f6WtaR/nyYdj0Eoy+C6b9CfyDf/y+jw/0HQ+3r4L35sHi+XDp0zD6DnPidXN/+jSL7MMnnXrO5MhQHr18cKv75ebm8vrrrzNhwgRuu+02nn32WV5++WXWrFlDYmIiP/vZz3jppZf49a9/DUBsbCwZGRncf//93HLLLXz33XfU1NQwZMgQ7r77bqd+D/pOX9M6QvoCI+GP/SVc+tRPE35znbrDzz4x2vaX/zcc3NBxcWpOERMTw4QJEwC46aabWLNmDfHx8SQmJgIwf/58vvnmm//sP3v2bABSUlIYM2YMISEhREREEBgYyIkTJ5wam77T1zRX2/+1kbwHTIPpf3HsGP9gmPsK/OsC+OhOuPs7o91fc5gjd+Su0nJYZdeuXSkvLz/r/oGBgQD4+Pj853nT64aGBqfGpu/0Nc2VqkphyXzoMQCuXgA+vo4fGxgCV74CFYXwxUOui1FzukOHDrFhg/EJbeHChaSlpZGfn09eXh4A77zzDpMmTTIlNp30Nc2VVj0CtVVw7Tttu1OPHQMXPAAZ70LOp86PT3OJgQMH8sILL5CUlMTx48e5//77eeONN7jmmmtISUnBx8fH6W31jhKlrFXfLC0tTelFVDSPcPB7eGOmkbSnPtr28zTUwWtToLoc7tsOfoGtH+OlcnJySEpKMjsMlzvT9ykiW5VSaa0dq+/0Nc0VGhvg899CWAxM/G37zuUXANMeg5NFsO1t58SneS2d9DXNFTa/DKVZMONJCOjc/vP1uwhixsL6Z6Hec8sGa66nk76mOdvpE7DubzBgKgya5ZxzisDkh6DyMGx/xznn1LySTvqa5mxbXoXaCpjyiJGsnSV+EsSOh/XP6Lt9rc100tc0Z6o7BRteNMbk9xnm3HP/526/WLfta22mk76mOdPWN+H0sfZ33p5N/ESIGWP0GVhs5J3mHnTS1zRnaaiF7/8JfS+A2LGuu87IW6E8Dw5+57praB5LJ31Nc5aMfxtNLxN/49rrJM+BwDDY+pZrr6N5JIeSvojMEJFcEckTkQfP8P7dIrJLRDJE5FsRSbZvjxOR0/btGSLi/DqhmmYFNptxlx85AvpNdu21AjrB0GuN0svVx1x7Le285efnM2jQIG655RYSExO58cYbWb16NRMmTCAhIYHNmzezefNmxo0bR2pqKuPHjyc3NxeAN998k3vvvfc/55o1axbr1q1zanytFlwTEV/gBWAaUAhsEZFlSqnsZrstVEr9y77/bOBZYIb9vX1KqeFOjVrTrObAOji2H+a+6twRO2czcr4xSmjnYhhrznR+y1vxIJTscu45e6fAzL+2ulteXh5LlixhwYIFjBo1ioULF/Ltt9+ybNkynnjiCd5++23Wr1+Pn58fq1ev5g9/+AMffPCBc2M9C0eqbI4G8pRS+wFEZBEwB2OxcwCUUs2LVncGdA+T5l3SFxgrXiXN7pjr9U4xPlVsfRPG3NUxf2g0h8XHx5OSkgLA4MGDmTJlCiJCSkoK+fn5VFRUMH/+fPbu3YuIUF9f32GxOZL0o4CCZq8LgTEtdxKRe4AHgADg4mZvxYvIduAk8LBSan3bw9W0s7PZFLuKKujWKYDobsH4+HRQIjxZbCx1OO6X4B/UMdcEGHkLfHofFG6BmNEdd1134cAduau0LI/cvHRyQ0MD//M//8PkyZP56KOPyM/P56KLLgLAz88Pm832n2Nrapw/H8Np9fSVUi8AL4jIDcDDwHygGIhVSpWLyEjgYxEZ3OKTASJyJ3AnGCvIaNr5OFFdx+L0At7deIhDx6oBCPb3ZVhMGI/NGUJiLxcvMr79HVCNxqiajjTkKqPkcsZCnfTdTEVFBVFRUYDRjt8kLi6OF198EZvNRlFREZs3b3b6tR3pyC0CYpq9jrZvO5tFwBUASqlapVS5/flWYB+Q2PIApdQrSqk0pVRaRESEo7FrGgeOnmLGc+t5YvlueocF8cw1w/jbVSlcPzqWvNIqZj//LQs3HcJl1WQbG4wmln6ToUd/11zjbAK7QOIlRsllW2PHXltrl9/97nc89NBDpKam/miRlAkTJhAfH09ycjL33XcfI0aMcPq1Wy2tLCJ+wB5gCkay3wLcoJTKarZPglJqr/355cCjSqk0EYkAjimlGkWkH7AeSFFKnXXIgS6trDnqwNFTzHtlA/WNigW3jGJ4TNcfvV9aWcNvFu9g/d6jXJcWw1+vSvnJikbttns5LLreqJef3EHt+c1lfWws0jL/M4i/sOOvbzG6tLITSisrpRqAe4GVQA6wWCmVJSKP2UfqANwrIlkikoHRrj/fvn0isNO+fSlw97kSvqY56mD5Dwn/vTvG/iThA/QMCeKtW0fzi4v68356Aa+tP+D8QLa9BV16G+vZmiFhGvgFG8M3Nc0BDrXpK6WWA8tbbHuk2fNfneW4D4COGYekeY1Gm+JXizKobbDx/p3jGNj77G32Pj7C7y4ZyMHyUzy5IodBfUK4MMFJTYhVZbB3FYy/F3z9nXPO8xXQ2Uj8Octg5lPgo+dbauem/4dobuet7/PJKDjBn2YPPmfCbyIiPH31MBJ7hXDvwu0cKq92TiCZHxgduEPnOed8bZU8B6qOQMFGc+OwCKutBuhs7f3+dNLX3ErBsWr+35e5XDQwgtnDIh0+rnOgH6/cnIbNpnj4k0znJIad7xvj5Xslt/9c7ZF4CfgF6SYeICgoiPLyco9N/EopysvLCQpq+9Bgpw3Z1DRXU0rxx48zAXj8yvPvlI3t0Yn7pyXy2GfZrMo+wvTBvdsezNG9cHgbTP9L28/hLIEhxoIt2cvgkie9uoknOjqawsJCysrKzA7FZYKCgoiOjm7z8Trpa25j7e5SvtlTxqOXJxPVNbhN57h5XF8WbTnEY59lMzExgiB/37YFs2MRiA+kXNO2450teQ7s/syYqBX7k7mTXsPf35/4+Hizw7A0770l0NyKUop/rs0julswN43t2+bz+Pv68L+zB1N4/DQvf72/bSex2WDXYmPd2pB2fFpwpsRLwMcPcpe3vq/m1XTS19zCd3nlZBSc4BcX9cfft33/bcf3D+eyoX14cV0eJRVtmOZesBFOHDK/A7e5oDCIHQd5q82ORLM4nfQ1t/B/a/fSOzSIq0e2vS2zuQdnDKLBpnj5m33nf/DOxeDfCQZd5pRYnGbAVDiSCScPmx2JZmE66WuWt2l/OZsPHOOuSf0I9GtjG3wLMd07ccXwKN7bfIiyylrHD2xsMMbEJ84wyiBYScI041Hf7WvnoJO+ZnkvrNtHeJcA5o1ybjG+X07uT22Djde/PY+Zuge/hepyGHyFU2Nxip7JEBJpTBjTtLPQSV+ztANHT/HNnjJuGR9HcIBz7vKb9I/owmUpfXhnQz4nquscOyjrY6NpZ8A0p8biFCKQMBX2r4PGjqvPrrkXnfQ1S3tv8yH8fIRrR8W0vnMb3HvxAE7VNbLgu/zWd27etBPQySXxtFvCdKg9CQXOL8mreQad9DXLqm1oZOnWQqYm9aJniGsWJxnUO5Rpyb14e0M+NfWtlCe2ctNOk/hJxtDNvV+aHYlmUTrpa5a1MusIx07VccMY1y6sc+uEOE5U17Mso5VRL1Zu2mkSFKqHbmrnpJO+ZlkLNx0ktnsnLhgQ7tLrjOvXg4G9Qnjz+/yz12xxh6adJnropnYOOulrlrSvrIqN+48xb3SMy9e6FRF+Nr4v2cUnST94/Mw7uUPTTpMBU43HfV+ZG4dmSTrpa5a0eEsBfj7CNSNd04Hb0pWpUYQG+fHm9/ln3iF7mfWbdpr0TIZO4XDgG7Mj0SxIJ33Nchptik8yDnPRwJ5EhAR2yDU7BfhxbVoMX2SW/LQ0g80Guz+HAVOs37QDRpXN+AuNpO+hJYa1ttNJX7OczQeOUXKyhjnDHa+X7ww/GxeHTSkWbj704zeKtkJVCQy6vEPjaZf4iVB5GMrbUGZC82gOJX0RmSEiuSKSJyIPnuH9u0Vkl4hkiMi3IpLc7L2H7Mflisglzgxe80yfZBTROcCXqUm9OvS6sT06cWFCBEvSC2i0NbtD3v2pMQwycXqHxtMu8ZOMxwNfmxuHZjmtJn0R8QVeAGYCycD1zZO63UKlVIpSajjwFPCs/dhkYB4wGJgBvGg/n6adUW1DI8t3FXPJ4N5On4HriOtHxVBcUcM3e+2LcCgFOZ9B3IUQ3K3D42mz7v0gNEq362s/4cid/mggTym1XylVBywC5jTfQSl1stnLzkDTbdIcYJFSqlYpdQDIs59P087o69wyTtY0MLuDm3aaTEnqRY/OAby/ucDYULYbju2DpFmmxNNmIkYTT/56o09C0+wcSfpRQEGz14X2bT8iIveIyD6MO/37zudYTWvyScZhenQOcPnY/LMJ8PPhqpHRrM45YlTfzPnMeGOgxcooOyJ+ojHMtDTb7Eg0C3FaR65S6gWlVH/g98DD53OsiNwpIukiku7Ja1tq51ZZU8/qnCPMGtoHv3YulNIe16bF0GBTfLit0GjPjx4FoX1Mi6fN4icaj7qJR2vGkd+sIqD5YOlo+7azWQQ0zWBx6Fil1CtKqTSlVFpERIQDIWmeaHXOEWobbKY17TQZ0LMLo+K6sW7TVijeAYPcrGmnSVg0dO+vk772I44k/S1AgojEi0gARsfssuY7iEhCs5eXAXvtz5cB80QkUETigQRAl//TzmjFrhJ6hQaSGmN+h+l1o2IZWLHeeOGuSR+Mu/2D3xllJDQNB5K+UqoBuBdYCeQAi5VSWSLymIjMtu92r4hkiUgG8AAw335sFrAYyAa+AO5RSrVSylDzRqdqG/h6Txkzh/RxedkFR8wc0psZfts4EtgXwgeYHU7bxV9olFou3mF2JJpF+Dmyk1JqObC8xbZHmj3/1TmOfRx4vK0Bat5hXW4ZtQ02ZgzpbXYoAHS2VTJKcnizZhY31DWaMnzUKfpOMB4PfQ/RI82NRbMEPSNXs4TlmcWEdwlgVFx3s0Mx7F2FL418XpfKl9klZkfTdiG9jTH7B783OxLNInTS10xXU9/IV7tLmZbcG18LNO0AkLsc1bknpaEpLN1aaHY07dN3vJH09Xh9DZ30NQv4Zk8Z1XWNXJpijaYdGmph72pk4AyuHBnDd3lHf1qEzZ3EjoeaE8ZEM83r6aSvmW5FZglhwf6M7dfD7FAM+euhrhIGXsbcEdHYFHy0/VyjlC2u73jj8eB35sahWYJO+pqp6hpsrM45wrTkXvibOCHrR3YvN2rn95tEfHhnRvbtxgfbCs++qpbVdYuDkEg4tMHsSDQLsMhvmeatNh0op7KmgUsGW6RpRynIXQH9Lwb/YADmjogir7SKzKKTrRxsUSLQd5zRru+uf7g0p9FJXzPVquwjBPv7cmGCObV2fqI4w6hDP+iHWjuzUiIJ8PVx/yaeymI4fsDsSDST6aSvmUYpxZdZR5iYGE6Qv0XGwe9eDuIDCT8s/RDWyZ+LB/Vk2Y7DNDS66QiYpvH6B3UTj7fTSV8zza6iCkpO1jAt2SJNO2A07cSMhc4/7lS+IjWKo1W1rM87alJg7RQ+0FgPQI/X93o66WumWZV9BB+BKYN6mh2K4fhBOLILBl36k7cmD4ogLNifj921icfHxxi6qUfweD2d9DXTfJl1hFFx3enWOcDsUAx7vjAeB/406Qf6+TJraB9WZpVQVeumxcv6jjPa9CvdeIax1m466WumOFh+itwjlUy3yqgdgNzlEJ4IPfqf8e25I6KoqbfxRaabJs3YccbjoY3mxqGZSid9zRSrso8AMD25Yxc/P6vTJyD/2zPe5TcZEduN2O6d3LeJp/dQ8AvSSd/L6aSvmeLLrCMM6h1CTPdOZodiyFsNtoZzJn0R4YrUKL7bd5QjJ92wLINfAESlQYFO+t5MJ32twx07VUf6wWPWucsHY9RO5wiITjvnblcMj0QpWJZxuIMCc7LYMVC8E+pOmR2JZhKd9LUOt3Z3KTaFdYZqNtbD3lWQeAn4nHu+QL+ILgyL6eq+E7Vix4FqhMJ0syPRTKKTvtbhVmWX0Ds0iCFRoWaHYsj/Fmorztm009yVwyPJLj5JbkmliwNzgehRgEDBJrMj0Uyik77WoWrqG/lmz1GmJvdExDq18/ELhn6THdp91rBIfH3EPe/2g7tCzyTdmevFHEr6IjJDRHJFJE9EHjzD+w+ISLaI7BSRNSLSt9l7jSKSYf9a1vJYzbt8l3eU0/WN1mnaUcoovdB/MgQ41qkc3iWQiQnhfJJRhM3mhgXMYsdCwWaw6eWqvVGrSV9EfIEXgJlAMnC9iCS32G07kKaUGgosBZ5q9t5ppdRw+9dsNK+2KvsIXQL9GNvPIssiluyEk4UON+00uSI1iuKKGjYeKHdRYC4UM9ZYL6A02+xINBM4cqc/GshTSu1XStUBi4A5zXdQSn2llKq2v9wIRDs3TM0T2GyK1TmlTBoYQaCfhQqsIZA447wOm57cm84Bvu45Zj92rPGom3i8kiNJPwooaPa60L7tbG4HVjR7HSQi6SKyUUSuONMBInKnfZ/0srIyB0LS3FFG4QmOVtVabKjm5xAzBrpEnNdhwQG+zBjShxW7Sqipd7Nmkq6xENJHJ30v5dSOXBG5CUgDnm62ua9SKg24AXhORH4yx10p9YpSKk0plRYRcX6/fJr7+DLrCH4+wkWJFimwduIQlJy5wJoj5o6IorK24T+zi92GiPGHTo/g8UqOJP0iIKbZ62j7th8RkanAH4HZSqnapu1KqSL7435gHZDajng1N/Zldglj+/UgrJO/2aEYdi83Hgdedu79zmJsvx70Dg1y3yaeigKocMPYtXZxJOlvARJEJF5EAoB5wI9G4YhIKvAyRsIvbba9m4gE2p+HAxMA3XvkhfJKq9hfdopLBlusaSc8EcIHtOlwXx9hTmokX+8po7yqtvUDrCRmjPGoSzJ4nVaTvlKqAbgXWAnkAIuVUlki8piINI3GeRroAixpMTQzCUgXkR3AV8BflVI66XuhlVlGZcqpVmnPP33cWFDkPEfttDQ3NZoGm+LTHW5WlqF3irH4+yHdxONt/BzZSSm1HFjeYtsjzZ5PPctx3wMp7QlQ8wxfZh9hWHQYfcKCzQ7FsGelUWAt6fJ2nWZg7xCS+4Ty0fYibpkQ76TgOoCvP0SN1O36XkjPyNVcrqSihh0FJ6xVOz/nU2MES+SIdp9q7ogodhRWkFda5YTAOlDMGKMju9bN4tbaRSd9zeVW5Visdn5dNeStgUGXGcsIttPs4UZZhg+3FTohuA4UO9Yovla01exItA6kk77mcl9mlRAf3pkBPbuYHYph/1fQcBoGzXLK6XqGBDExIZyPtrtZWQZdfM0r6aSvuVRFdT0b95czPbmXdQqs5XwGQV0h7gKnnfKqkdEUV9SwYb8blWVoKr6mk75X0Ulfc6nVOUeob1TMTOljdiiGxgbYs8Iou+DrvPkCU5N6ERLkxwfu1sQTMwYKtoDNZnYkWgfRSV9zqRWZxUSGBTEsOszsUAwHvzOGayY5p2mnSZC/L7OGRvJFZgmnahucem6XihljrCVQlmN2JFoH0Ulfc5nKmnq+2XOUGUP6WKdpZ/fnxuLg/S92+qmvGhFFdV0jX2SWOP3cLhNrn6Sl6/B4DZ30NZdZu7uUukYbl6ZYZKimzQa7P4P+UyCgs9NPP7JvN/r26MTSrW7UxNMtHjr3NOrra15BJ33NZVbsKqFnSCAjYruZHYrh8DY4WQTJrlnWQUS4ekQ0G/aXU3CsuvUDrEDEuNs/tMHsSLQOopO+5hKnahv4KreUGUN64+Njkaad7E/Ax/+8a+efj6tGRiOCe93tx4yFEwfhZLHZkWgdQCd9zSXW5ZZR22Bj5hCLjNpRCnKWQb9JxlBFF4nsGswFA8JZurXQfcbsx44zHnXxNa+gk77mEst3FdOjcwCj462yLOIuOJ4PSa5fsfPqkdEUnTjNRncZs99nqLEwvC6+5hV00tecrqq2gTW7j3BpSh98rdK0k7MMxMcoveBilwzuTUiQH0vcpYnH1x+i03S7vpfQSV9zulXZJdTU25gzPNLsUH6QvQz6ToDO4S6/VJC/L5cPi2RFZjEna+pdfj2niB2ri695CZ30Naf7JOMwUV2DrTNqpywXjuZC8pwOu+Q1I6Opqbfx2Q436RyNaSq+lm52JJqL6aSvOVV5VS3r9x7l8mGRFhq1swwQpxVYc8TwmK4k9urC++kFHXbNdomxF1/Tk7Q8nk76mlMtzyyh0aaYPcxKTTufGOUGQjtuJJGIcN2oWHYUnCCn+GSHXbfNgsKg12Cd9L2AQ0lfRGaISK6I5InIg2d4/wERyRaRnSKyRkT6NntvvojstX/Nd2bwmvV8mnGYhJ5dSOoTYnYohqN5cGQXDL6iwy99ZWoUAb4+vL/FTe72Y8dC4RajKJ3msVpN+iLiC7wAzASSgetFJLnFbtuBNKXUUGAp8JT92O7Ao8AYYDTwqIhYpKFXc7aiE6fZnH+MOcMjrVNrJ/sj47ED2/ObdO8cwPTBvfg4o4ia+sYOv/55ix0HdVVwJNPsSDQXcuROfzSQp5Tar5SqAxYBP/oNUkp9pZRqmne+EYi2P78EWKWUOqaUOg6sAlw3HVIz1cfbiwC43EpNO1kfG52UoebENG9ULCeq6/ky+4gp1z8vMfbia7q+vkdzJOlHAc0/nxbat53N7cCKNh6ruSmlFEvSCxgd352+PZxfzKxNju417loHX2laCOP79yC6WzDvbzlkWgwO6xoDodFw8HuzI9FcyKkduSJyE5AGPH2ex90pIukikl5WVubMkLQOsiX/OPnl1VybFmN2KD/I+th4dFGBNUf4+AjXpsXwXV45+UdPmRaHw/qOMyZpKTcpIaGdN0eSfhHQ/Dc52r7tR0RkKvBHYLZSqvZ8jlVKvaKUSlNKpUVERDgau2Yhi9ML6Bzga50yygBZH5natNPkulEx+PoI77nD3X7f8VB1BI7tNzsSzUUcSfpbgAQRiReRAGAesKz5DiKSCryMkfBLm721EpguIt3sHbjT7ds0D1JV28DyXcVcPiySTgF+ZodjKNsDpVmmNu006RUaxJRBPVmaXkhdg8WXJew7wXg8+J25cWgu02rSV0o1APdiJOscYLFSKktEHhORps/NTwNdgCUikiEiy+zHHgP+jPGHYwvwmH2b5kGW7yymuq6Ra6zUtJNtftNOczeMiaX8VB0rsyy+qlZ4InTqodv1PZhDt2VKqeXA8hbbHmn2fOo5jl0ALGhrgJr1LU4voH9EZ0bEuq5k8XnL/NAYgmhy006TiQkRRHcLZuGmQ9Ya3dSSiNHEo+/0PZaekau1S25JJekHj3NtWox1xuYfyTYW+h5yldmR/IePj3D96Fg27C9nf5nFi5r1vQBOHIITbjKpTDsvOulr7fLWhnwC/XysNWon8wOjjLIJE7LO5Zq0aPx8hPc2W7xDt+9441GXWvZIOulrbVZRXc9H24qYMzySbp0DzA7HoJSR9OMnQpeeZkfzIz1Dgpg+uBdLthZae4Zur8EQGKabeDyUTvpamy1OL+B0fSPzx8eZHcoPijPg+AFLNe00d/PYOE5U17Nsx2GzQzk7H1+jDo/uzPVIOulrbdJoU7y9MZ/Rcd0ZHBlmdjg/yPzAWPy8A8son4+x/bqT2KsLb2/IR1l5AlTf8XB0D1TpyZKeRid9rU3W7i6l4Nhpa93l22yQ+REMmAKdLLI2bwsiws3j4sgsOsn2ghNmh3N2TeP1D+m7fU+jk77WJgu+PUDvUKON2jIKN8NboJuNAAAgAElEQVTJQss27TS5MjWKLoF+vLPhoNmhnF3kcPDvBPm6Xd/T6KSvnbetB4+zYX85t18Qj7+vhf4L7VoKfkEwcKbZkZxTl0A/rhoRxec7izlaVdv6AWbw9Tfa9fPXmx2J5mQW+o3V3MULX+XRrZM/N4yJNTuUHzQ2GLV2EmdAoEUWcDmHm8f1pa7RxiIrD9+Mnwil2VBV2vq+mtvQSV87L5lFFazdXcrtF8TTOdAidXYADqyD6qOQco3ZkThkQM8QLkwI552NB61bjyd+ovGo7/Y9ik762nl5cV0eIYF+3DwuzuxQfmzXUmNsecI0syNx2G0XxHPkZC0rMovNDuXMeg8zfqYHvjE7Es2JdNLXHJZXWsmKzBLmj48jLNjf7HB+UH8acj6F5MvBL9DsaBw2KSGCfuGdWfDtAWsO3/T1M4Zu6qTvUXTS1xz2/1buoZO/L7ddEG92KD+2Z6WxtqubNO008fERbp0Qx47CCrYdsujwzfiJRm39ikKzI9GcRCd9zSFb8o/xRVYJd0/qT3erlFxosmsJdOkFcReaHcl5mzsimtAgPxZ8d8DsUM6sqV3/gG7X9xQ66WutUkrxl89z6BUayM8v7Gd2OD92+gTs/dIYm+/ja3Y0561zoB/zRsfyRWYJhcerzQ7np3omG/X1dROPx9BJX2vVZzuL2VFwgt9OH0hwgMUSa86n0FgHQ642O5I2u2V8HAIs+Dbf7FB+ysfH+AR14Bu9bq6H0ElfO6ea+kb+9sVukvqEMndEtNnh/NSuxdC9H0SNMDuSNovsGszlwyJZtOUQFdX1ZofzU/ETjZnOet1cj6CTvnZOz6/No/D4aR6+LAlfH4ssktLk5GGjrTnlWmPFJzd2x4X9qK5r5N1NFizNED/JeDzwtblxaE7hUNIXkRkikisieSLy4Bnenygi20SkQUSubvFeo33d3P+snau5h90lJ/nX1/uYOyKKCQPCzQ7np3YtBRQMvdbsSNotOTKUCxPCefP7fGobLFZrv0d/CI2GfWvNjkRzglaTvoj4Ai8AM4Fk4HoRSW6x2yHgFmDhGU5xWik13P5ljVWqtVY12hQPfrCL0GB/Hr6s5T+3RexaDFEjjaTkAe6c2I+yylo+3l5kdig/JgIDLob9X0OjBZuftPPiyJ3+aCBPKbVfKVUHLAJ+tA6dUipfKbUTsOh8cu18vbMhn4yCEzwyK9l6QzQBSndDyS6jacdDXDAgnOQ+obz8zX4abRbrNO0/BWpPQmG62ZFo7eRI0o8Cmq+QXGjf5qggEUkXkY0icsV5RaeZIq+0kr9+sZtJiRHMGR5pdjhntmsxiC8MmWt2JE4jIvxycn/2l53ii8wSs8P5sX4XGT/vfWvMjkRrp47oyO2rlEoDbgCeE5GffBYXkTvtfxjSy8r0Sj1mqqlv5N6F2+kU4MfTVw9FrNhBarPBziXQf7Ll1sFtr5lD+tAvojPPf5VnrdIMwV0hOg3ydNJ3d44k/SIgptnraPs2hyiliuyP+4F1QOoZ9nlFKZWmlEqLiIhw9NSaCzy5PIfdJZU8c80weoYGmR3OmRVsgopDHtW008TXR/jFpP7kFJ/kq1yLlTTuPwUOb4dT5WZHorWDI0l/C5AgIvEiEgDMAxwahSMi3UQk0P48HJgAZLc1WM21vsgs4a0NB7n9gngmD7LwHfTORcaqToMuMzsSl7giNYqorsE8v9Zid/sDpgIK9n9ldiRaO7Sa9JVSDcC9wEogB1islMoSkcdEZDaAiIwSkULgGuBlEcmyH54EpIvIDuAr4K9KKZ30LSi3pJLfLM5gWHQYv5sx0Oxwzq6h1lgsJelyCOxidjQu4e/rw90X9WfboRNs2Gehu+rI4RDcTTfxuDmHVsFQSi0HlrfY9kiz51swmn1aHvc9kNLOGDUXO36qjjveTqdToB8v35xGoJ/FSi00t2cl1FTA0OvMjsSlrhkZzfNr9/Lc6r2M69/DGn0rPr7Qb7LRmauU20+I81Z6Rq6Xq2+0cc/CbZRU1PDyzSPpHWbRdvwmO983Kmo2zRL1UEH+vtw7eQCb84/xbd5Rs8P5wYCpUHUEjmSaHYnWRjrpezGlFA99uIvv95XzxNwURsR2Mzukc6s+Ztzpp1xjLPDh4a4dFUNkWBDPrtpjnbb9AVOMxz0rzY1DazOd9L3Y31ftYenWQn41JYGrR1qwmFpLWR+Brd4jyi44ItDPl3svTmD7oROs22ORocwhvSFyBOQub31fzZJ00vdS/950kP9bm8d1aTH8emqC2eE4Zuf7EJEEvYeaHUmHuSYtmpjuwfzdSnf7gy6Foq1QabEJZJpDdNL3Qst3FfM/H2cyeWAEj185xBqdhK05tt8Ynz/U/Stqng9/Xx/+6+IEdhZWWGeW7sBLjcfcFebGobWJTvpeZv3eMn61aDsjYrvx4o0j8fN1k/8CO94HxGuadpqbmxpFQs8uPLUyl/pGC5S36pkMXfvqpO+m3OQ3XnOG7YeOc9c7WxnQM4TXbxllvVWwzkYp2PGesZhHmBv0PTiZn68Pv58xiANHT/H+loLWD3A1EeNuf/86qDtldjTaedJJ30vkllRyyxtbiAgJ5K3bRhEW7G92SI47tBFOHIRh15sdiWmmJPVkdFx3nlu9l1O1DWaHAwNnQmOtrrHvhnTS9wKHyqu5+fVNBPn78O7tY+gZYvGx+C3teM8ou5B0udmRmEZE+P3MQRytquW19QfMDgf6joegMN3E44Z00vdwpZU13LxgE3WNNt65fQwx3TuZHdL5qT8NWR9D0myPLbvgqJF9uzFzSG9e/mYfJRU15gbj6w8J02HPF2Cz2Epf2jnppO/BqmobuO3NLZSerOWNW0aR2CvE7JDOX+4KqK2A4d7btNPcQzOTaLAp/vbFbrNDMdr1q8vh0AazI9HOg076HqquwcYv3t1KTnElL940glSrz7Y9mx2LIDQK4i40OxJLiO3RiTsujOej7UVsPXjc3GASpoNfsDFpTnMbOul7IKUUD364k/V7j/LXuSlMHmjhMsnnUnkE8lYbZRd83GSkUQf45UUD6BUayJ8+zcJm5rKKgV0g8RKj+a3RAp3LmkN00vdAz6/N48NtRTwwLZFr0mJaP8Cqdi0G1QjDbzA7EkvpHOjHgzMHsbOwgqVbC80NZshVUH0U8tebG4fmMJ30PcznO4t5ZtUe5qZG8V8XDzA7nLZTCjIWQlQaRFi4vr9Jrhgexai4bjy5Iodjp+rMCyRhGgR0gcwPzItBOy866XuQnYUn+M2SDEb27caTV6W4R3mFsynOgNJsSL3R7EgsSUR4/MoUKmsaeHJ5jnmB+AcbK5jlLIMGE//4aA7TSd9DlFfVctc7W+nROZCXbx5p7YVQHJGxEHwDYfBcsyOxrMReIdwxsR9Lthayab+JK2wNnmssbKOXUXQLOul7gIZGG//13naOnarj5ZtHEt4l0OyQ2qehFnYtgaRZENzV7Ggs7b6LE4juFswfP86ktsGk8fL9LzYmaukmHrfgUNIXkRkikisieSLy4Bnenygi20SkQUSubvHefBHZa/+a76zAtR88s2oP3+8r589XDGFIVJjZ4bRf7go4fVx34DogOMCXP18xhLzSKp5fm2dOEH4BxuS53cuNyXSapbWa9EXEF3gBmAkkA9eLSHKL3Q4BtwALWxzbHXgUGAOMBh4VETcdMG5Nq7OP8NK6fVw/OpZr3XmkTnMZCyEk0liPVWvV5IE9uWpENC+u20dmUYU5QQy9FuoqIedTc66vOcyRO/3RQJ5Sar9Sqg5YBMxpvoNSKl8ptRNoWff1EmCVUuqYUuo4sAqY4YS4NaCkoob/XrqD5D6hPHp5y7/DburkYchbBcPm6bH55+GRWcn06BzAb5fsoK7BhPLLfS8wyi1ve7vjr62dF0eSfhTQvJ5roX2bIxw6VkTuFJF0EUkvK7PIsnAW12hT3P9+BjX1Nv55QypB/h6SIDMWgrJB6k1mR+JWwjr588SVKewuqeT5tXs7PgAfHxhxszFe/9j+jr++5jBLdOQqpV5RSqUppdIiIiLMDsct/OvrfWzYX86f5gymf4SHFCKz2WD7O0bJhR79zY7G7UxN7sXcEVG8sG6fOSUaht8I4gPb3+34a2sOcyTpFwHNG4uj7dsc0Z5jtbPYWXiCZ1ftYdbQPlzjDguaO+rgt3A8H0b8zOxI3Nb/zh5M79Ag7n8/g6qOrrsfGgkDphqf1nRZBstyJOlvARJEJF5EAoB5wDIHz78SmC4i3ewduNPt27Q2qqlv5IHFO4joEsjjV7j5BKyWtr0NgWFeXTe/vUKD/Hlu3nAKj1fzv8uyOj6A1JuhstiomaRZUqtJXynVANyLkaxzgMVKqSwReUxEZgOIyCgRKQSuAV4WkSz7sceAP2P84dgCPGbfprXRM1/mkldaxd+uHkpYJzda/ao1p49D9jJjFIh/sNnRuLVRcd25Z/IAlm4t5NMdhzv24okzoHOE0UynWZKfIzsppZYDy1tse6TZ8y0YTTdnOnYBsKAdMWp2mw8c47VvD3DjmFgmJXpY38fOJcbye7ppxynum5LA+r1HeejDXaREhREX3rljLuwXYCxrueEFqCj0yjWNrc4SHbla607XNfLfS3cQ060Tf7g0yexwnEsp2Pom9BkGfYaaHY1H8Pf14fkbUvERuGfhNmrqO3C27qifAwo2v9Jx19QcppO+m/j76j0cLK/mr1el0DnQoQ9o7qNgE5RmwchbzY7Eo0R368Qz1w4n6/BJnujIomzd+hozdNPfhNqqjruu5hCd9N3AjoITvLZ+P9ePjmV8/3Czw3G+La9DQIixWIrmVNOSe3HHhfG8veEgn2R04MC5cfcay1xm/Lvjrqk5RCd9i6trsPH7D3bSMySIhy4dZHY4znfqKGR/bMzA9fKFz13ldzMGMTquOw9+sIvdJSc75qIxoyB6NGx8US+cbjE66Vvcy1/vY3dJJX+5YgihQR40WqfJ9nehsQ5G3W52JB7L39eH529MJSTIj7ve2UrF6fqOufC4e4x5F7nLW91V6zg66VvY/rIq/vlVHrOG9mFqci+zw3E+mw22vgF9J0BPD+uctpieIUG8dNMIDp84zf3vZ9DYEWvrDpoFXWPh+38anfWaJeikb1FKKf7w0S6C/Hx4xFOKqbW0b61xJ5h2m9mReIWRfbvzyKxk1u4u5Zkvc11/QV8/GH+f0VG/b63rr6c5RCd9i1qytZCN+4/x0KVJ9AwJMjsc19j8ijGRJ2m22ZF4jZvG9uX60bG8uG5fx3TsjvgZhMXC2j/ru32L0EnfgsqranlieQ6j4rpxnafUyG/p6F7Yu9K4y/cLMDsaryEi/Gn2YEbHded3S3eys/CEay/oFwiTfgeHt+u2fYvQSd+CHv88h1O1DTxxZQo+Ph5UW6e5jS+Bb4B9Io/WkQL8fHjxphGEdwnkjrfTKamoce0Fh10P3fvD2seNfhzNVDrpW8z3eUf5cHsRd0/qT0KvELPDcY3qY0YlxqHXQpeeZkfjlcK7BPLa/DSqahr4+dtbqK5zYVVMXz+Y/AdjAl7Wh667juYQnfQtpKa+kT9+nElcj07cM3mA2eG4TvoCaDgNY+8xOxKvltQnlH/MSyXr8EkeeH8HNleO6Bk8F3oOhjV/grpq111Ha5VO+hby4rp9HDh6ir9ckeI5K2G11FAHm1811r/t5aGjktzI1ORe/PHSJL7IKuGplS4c0ePjAzP/BicOwfpnXHcdrVU66VvE3iOVvLQujyuGR3JBggeWWmiS+QFUlRjT9DVLuP2CeG4YE8u/vt7Hos2HXHeh+Ath6HXw3T+gbI/rrqOdk076FmCzGWPyOwf68fAsD777tTXCt3+HnskwYIrZ0Wh2IsJjswczMTGChz/O5Nu9R113sel/gYBOsPw3eginSXTSt4BFWwrYkn+cP1yaRHiXQLPDcZ3sj+FoLkz8b/CkFb88gJ+vDy/ckMqAnl34xbtbXVejp0tPmPIIHPgGdi1xzTW0c9JJ32SllTU8uSKHsf26e9Z6ty3ZbPD10xA+EJLnmB2NdgYhQf68fssoggN8ufWNLa4byjnyVogeBZ//Fo4fdM01tLNyKOmLyAwRyRWRPBF58AzvB4rI+/b3N4lInH17nIicFpEM+9e/nBu++3v0kyxqG2w8caWHrXfbUs4yKMsxJur4eGgntQeI6hrMgltGcfJ0Pbe+ucU1i6v7+MLcVwEFH/wcGjuoAJwGOJD0RcQXeAGYCSQD14tIy4bn24HjSqkBwN+BvzV7b59Sarj9624nxe0RVuwqZkVmCb+akkC/CA8uK2yzwTdPQ48EGHyl2dForRgSFcYLN45gz5FKfvHuVuoaXDChqns8XP4cFG6Gr55w/vm1s3LkTn80kKeU2q+UqgMWAS0/n88B3rI/XwpMEY++bW2/E9V1/M8nWQyODOXOif3MDse1dn8GRzKNtnx9l+8WLhrYkyevTGH93qM8+MFOlCs6XYdcZdTm+fbvkLfa+efXzsiRpB8FFDR7XWjfdsZ9lFINQAXQw/5evIhsF5GvReTCdsbrMf78WQ4nqut46uqh+Pt6cNdKQy2sesRoyx9yldnRaOfh2lExPDAtkQ+3F7luDP+MvxmjuRbfAsU7XXMN7UdcnW2KgVilVCrwALBQREJb7iQid4pIuoikl5WVuTgk863JOcIH2wq5e1J/BkeGmR2Oa21+BY4fgBlPGNPxNbfyXxcP4IYxsby0bh9vfHfA+RcI6AQ3LoGgUPj3NcbkLc2lHEn6RUDzUo/R9m1n3EdE/IAwoFwpVauUKgdQSm0F9gGJLS+glHpFKZWmlEqLiIg4/+/CjZRX1fL7D3aR1CeU+6YkmB2Oa506Cl8/BQnTYcBUs6PR2qBpDP/05F786dNs15RjDouCmz4wSnO8e5VRm0lzGUeS/hYgQUTiRSQAmAcsa7HPMmC+/fnVwFqllBKRCHtHMCLSD0gA9jsndPejlOKPH2Vy8nQ9z147jAA/D27WAfjqcag7BdMfNzsSrR38fH34v+tTGRPfnd8s3sG63FLnX6RnEsx7zxjC+cZMqOjARdy9TKtZx95Gfy+wEsgBFiulskTkMRFpWv3idaCHiORhNOM0DeucCOwUkQyMDt67lVJe+2f8o+1FfJFVwgPTE0nq85NWLs9SvBO2vgmj74CIn3y409xMkL8vr85PI7FXCL94dxtb8l3waxw3wbjjryiCBZcYay5oTicu6ZVvh7S0NJWenm52GE6Xf/QUs/75LYN6h/D+XePw9dQ6+WB03r56MVSVwj2boFN3syPSnKSsspbrXt5AWWUt7905liFRLuiTKt5hNPMoG1zzllGzR2uViGxVSqW1tp+Hty9YQ21DI/e+tw1fH+G5ecM9O+EDrHvSGKI5+5864XuYiJBA3v35GEKD/bn59U3sPVLp/Iv0GQa3rYTg7vD2bGMmt158xWl00u8ATy7fTWbRSZ6+eijR3TqZHY5rHdpkVFFMvRkGzjA7Gs0FIrsG8++fj8HP14cbXttEXmmV8y/Soz/c+ZUxzPerv8C7c6Gi0PnX8UI66bvYil3FvPl9PrdNiGf64N5mh+NatZXw0V0QGg2X6FmWniwuvDMLfz4GpRTXv7qRfWUuSPyBIUa5hlnPQcEmeH40bHjRqNaqtZlO+i6UU3yS3yzZwbCYrvx+5kCzw3GtxnpYPN8YZ33lS8a4a82jJfQK4b07xhqJ/xUXJX4RSLsVfrnR6Ohd+RC8Otmo0qm1iU76LlJeVcvP30onJMiPV24eSaCfB5cfUAo+fwD2rTHqqcRdYHZEWgdJ6BXCwjvGYlOK617eQE6xi0oyd+sLNyyGa940xvG/dbnR2VuyyzXX82A66btAXYONX/57G2VVtbxycxq9QoPMDsm11j8D296GC39r1FLRvEpiL2NEmr+vD9e9vIHth4675kIiRsG+e9Nh2p+hMB3+dQG8f5Mu4XAedNJ3skab4oHFGWw6cIynrhrKsJiuZofkOjYbrHkM1v4ZUq6Fix82OyLNJP0jurD4rnF06xzAja9t4itXTOBq4h8EE+6DX2XApN/D/m/g5Qth4Two2Oy663oInfSdSCnF/y7L4rOdxTw0cxBXpLasS+dB6qphyXzjLn/EfLjiRb0alpeL6d6JJXeNI65HZ37+Vrpr19sFCO4Gk/8Av94JF/0BCjbC69Pgjctg72q9HONZ6KTvJEopnvlyD+9sPMhdk/px16T+ZofkOgc3wOvTIedTo8TC5f8AX3+zo9IsoGdoEIvvHseEAeE8+OEunvpiNzabi5NvcFe46PdwfxZc8qRR4O/fV8HLEyHzQz3apwU9I9cJlFI8sTyHV9cfYN6oGJ6c66GrYJXvM5pzsj+GkD7GUDo9Fl87g/pGG//zcSaLthQweWAEz81LJSy4g24MGupg5/vw3XNQngfd+8MF98PQ68AvoGNiMIGjM3J10m+nRpvi4Y938d7mAuaP68ujlw/Gp6Nn3NpscOIgHN0DZblQdQROHze+lA18/IyvoDBjhmynHtClN4T0NpJ35x4QGAY+zT74KWWUUTieD/nrjeUOi3eAXzBM+JXRphrQuWO/T82tKKV4d9Mh/rQsi+huwbx008iOrTllazQ+ja5/Bkp2GvNHLvi1MXHQ3/MGV+ik3wEqa+p5YPEOVmUf4Z7J/fnt9IEdc4evFJTthrw1kP8tHPoeaip+eN+/s9HeGRRmrFRlawRbvbFPdTnYzrDuqfgaY+vFBxCorza+mkSPgqTZkHINhPZx+beoeY70/GP84t/bqDhdz4MzBnHL+LiOvTFSyliZ65unjUleXXrB+PuM8f8edOOik76L7S+r4o6308kvr+bhy5K4dUK86y9asgsyPzDuXsrzjG3d+xnj4qNGQkSSUdEyuNvZz6EU1J6EyhKoLIbKI1B91Kh9X3vS3vmlwL8TdO0L3eKg9xAIjXT996d5LGMdiZ2szinlwoRw/nrVUKK6BndsEEoZN0nfPA0HvjY+8Y79BYy6w+gXcHM66buIUoqPthfx6CdZ+Pv58PwNqYzvH+66C1YeMdondyyC0iyjmSbuQkiaBYkzjQUoNM0NNDX3PPF5DiLwwLREbhkfh58Zy4UWbIZv/h/sXQmBoTDqdhjzCwjp1fGxOIlO+i5QerKGP3y0i9U5pYzs241/zBvumgJqjQ3G7NZtb0PuClCNRvPKsHkweK6uXKm5tcLj1Tz6SRZrdpcyqHcIv58xiIsGRpgz+KF4B6x/FrI/Ad8AGH4DjP2lW64BoZO+E9XUN7LguwO89NU+6hpt/PclA7l1QrzzSyQfOwDb34WMfxtNL50jYNj1RseTG/4n1LSzUUqxMquEJ1fs5mB5NaPju3P/1ETG9utuTvIv3wff/x9kvAeNtTBgGoy9G/pd/OMBDhamk74TVNc18MHWQl5ct4/iihqmJvXkD5cm0S+ii/Mucvo4ZH1sNOEc2mB0pA6YaiT6xBkePcRM0+oabLy/5RD/WJPH0apahkSFcseF/ZgxpLc59aqqymDrG7D5VThVCl1jIfVnkHqj5fu1dNJvh7zSSpZsLWTR5gIqTtczPKYrD84cxNh+PZxzgVPlkPs5ZC+D/euMkTXhA2HotcadvW6n17xMTX0jH20v4tX1+9lfdoqwYH9mD4vkyhFRDI/u2vHDoBtqjQET296yV/QU6DsBhsw1RrF1iejYeBzg1KQvIjOAfwC+wGtKqb+2eD8QeBsYCZQD1yml8u3vPQTcDjQC9ymlVp7rWmYk/UabYldRBev3lPH5rmJ2l1TiIzBjSG9umxDPyL7d2veRs64airYaIwby1sDh7YAy7iKSZkPK1dBnuC5joHk9m03xbd5Rlm4tZGVWCbUNNsK7BDJlUE8uTAxndFx3enZ0AcNj+2HnEshcasyFQSAyFRKmQ79JEDnCEuP+nZb0RcQX2ANMAwqBLcD1SqnsZvv8EhiqlLpbROYBVyqlrhORZOA9YDQQCawGEpVSZ50X7eqkX99oo+BYNXmlVWQWVbCzqIJtB49zssYYu54a25XZwyK5LKXP+f/nUsoYB390rzHS5kiWkeBLdhlj48UHotJgwBSj6abPMJ3oNe0sKk7X89XuUlblHOGb3DIqa43f0bgenRgSFUZyZChJvUOJC+9MdLdg/F09Ckgp43c6dwXs/RIKtwDK6ACOTDVu3HoPgV6DjVnAHTwM1JlJfxzwv0qpS+yvHwJQSj3ZbJ+V9n02iIgfUAJEAA8237f5fme7XluTfl2DjV1FJ6isaaCqtoHKmgaOV9dx/FQdZZW1FFfUUHKyhsLjp2m01wLxEaMs7PCYrkwYEM64/j0I7xJo/OOePm58xGuoMR7rT0HdKaitgpoTcPoEnD5mjHevKoWTh43Zq3XN1gwNDIXeQyF2DMSMgZjR5x5Dr2naGdU32sg6fJLNB8pJzz9O1uGTFJ04/Z/3/XyEXqFB9AkLoldYEOGdA+jWOYCuwf6EBPkTEuRHl0A/ggJ8Cfb3JbJrcPvLQlQfg0Mbjb64gk1QkmnkiSbB3Yx5LiGRxlDQLr2MdX+DuxoTJwO6QEAnYzKlXwD4BYF/sPFeGzia9P0cOFcUUNDsdSEw5mz7KKUaRKQC6GHfvrHFsS5psD5ZU89VL/30b0mQvw/hXQLpExZESlQYs4b2IT68C/0jOjOodyjBAWfpLHqqH9BK05f4GCNsuvQ0Onn6jjf+kXv0h57JEBat7+Q1zQn8fX0YHtOV4TFduXOisa2iup69pZUcOHqK/PJTHD5RQ0lFDTmHT1J+qo6K0/VnPd9f56Ywb3Rs+4Lq1B0GXWp8gb0cSr7xaeDYAaPw2/GDxs1gwUajFaA1kSOMtYFdyJGk73Iicidwp/1llYjkOvP8uRAOHHXmOX9wAtjrmlO7hgt/Fm5H/yx+4FU/i+v/Btef/W0Tfxbr4K423yj2dWQnR5J+ERDT7HW0fYUT/BEAAAMvSURBVNuZ9im0N++EYXToOnIsSqlXgFccCbgtRCTdkY893kD/LH6gfxY/0D+LH3j6z8KRno8tQIKIxItIADAPWNZin2XAfPvzq4G1yugsWAbME5FAEYkHEgC9tI2maZpJWr3Tt7fR3wusxBiyuUAplSUijwHpSqllwOvAOyKSBxzD+MOAfb/FQDbQANxzrpE7mqZpmmtZbnKWK4jInfYmJK+nfxY/0D+LH+ifxQ88/WfhFUlf0zRNM7hHJSFN0zTNKbwm6YvI0yKyW0R2ishHIuL+qyacJxGZISK5IpInIg+aHY9ZRCRGRL4SkWwRyRKRX5kdk9lExFdEtovIZ2bHYiYR6SoiS+25Isc+OdWjeE3SB1YBQ5RSQzHKSjxkcjwdyl5O4wVgJpAMXG8vk+GNGoDfKKWSgbHAPV78s2jyKyDH7CAs4B/AF0qpQcAwPPBn4jVJXyn1pVKqaXHYjRhzBrzJaCBPKbVfKVUHLALmmByTKZRSxUqpbfbnlRi/2F5b2lREooHLgNfMjsVMIhIGTMQYjYhSqk4pdcLcqJzPa5J+C7cBK8wOooOdqZyG1ya6JiISB6QCm8yNxFTPAb8DbGYHYrJ4oAx4w97U9ZqIeM7K6XYelfRFZLWIZJ7ha06zff6I8fH+3+ZFqlmBiHQBPgB+rZQ6aXY8ZhCRWUCpUmqr2bFYgB8wAnhJKZUKnMJeNNKTWKL2jrMopaae630RuQWYBUxR3jdW1aGSGN5CRPwxEv6/lVIfmh2PiSYAs0XkUiAICBWRd5VSN5kclxkKgUKlVNOnvqV4YNL3qDv9c7EvBPM7YLZSqtrseEzgSDkNryDGijivAzlKqWfNjsdMSqmHlFLRSqk4jP8Ta7004aOU+v/t3aENQkEQRdE7ilIoAkkdCByGBqiDFmgBQwNIEjwGQxeDYC3BEIb8uUeOeuplssnuPoB7RMzHaMnrNYFJmdSm/8EemAGn8QvWOTM3tZF+591zGsWxqiyAFXCNiMuY7TLzWJhJ/2ELHMZidAPWxXm+zhu5ktRIm+MdSZKlL0mtWPqS1IilL0mNWPqS1IilL0mNWPqS1IilL0mNPAGKkseTGQE/VAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fd91cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.kdeplot(treino[treino.prognostico=='bom'].lpsa,label='bom')\n",
    "sns.kdeplot(treino[treino.prognostico=='mau'].lpsa,label='mau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,\n",
       "       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
       "       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
       "       1])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "classes = np.fromiter(map(lambda x: np.argmax(np.fromiter(x,dtype='float')),zip(*deltaK)),dtype='int')\n",
    "classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>12</td>\n",
       "      <td>38</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0             13    4\n",
       "1             12   38"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(classes,treino.prognostico)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Caso multivariado\n",
    "\n",
    "No caso em que $X$ corresponde a um vetor de atributos, usamos a distribuição normal multivariada, cuja função de densidade é dada por\n",
    "\n",
    "$$\n",
    "f_k(X) = \\frac{1}{2\\pi^{p/2}det(\\Sigma)^{1/2}}\\exp\\left({-\\frac{1}{2}(X-\\mu)^T\\Sigma^{-1}(X-\\mu)}\\right)\n",
    "$$\n",
    "\n",
    "Substituindo a função de densidade multivariada na equação de probabilidade condicional, podemos novamente constatar que seu argumento máximo é equivalente ao da função discriminante\n",
    "\n",
    "$$\n",
    "\\delta_k(X) = X^T\\Sigma^{-1}\\mu_k-\\frac{1}{2}\\mu_k^T\\Sigma^{-1}\\mu_k+\\log(p_k)\n",
    "$$\n",
    "\n",
    "Voltando ao exemplo, podemos agora incluir no modelo outros atributos como o volume e o peso da próstata. Felizmente não precisamos executar os cálculos manualmente. Podemos usar o modelo `sklearn.discriminant_analysis.LinearDiscriminant`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearDiscriminantAnalysis(n_components=None, priors=None, shrinkage=None,\n",
       "              solver='svd', store_covariance=False, tol=0.0001)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "attrs = ['lpsa','lcavol','lweight','lcp','lbph']\n",
    "model = LinearDiscriminantAnalysis()\n",
    "model.fit(treino[attrs],treino.prognostico.map(lambda x: x=='mau' and 1 or 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>18</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>7</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0             18    9\n",
       "1              7   33"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(treino[attrs]),treino.prognostico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0              6    4\n",
       "1              4   16"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(teste[attrs]),teste.prognostico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.77402017, 0.22597983],\n",
       "       [0.90268906, 0.09731094],\n",
       "       [0.80471737, 0.19528263],\n",
       "       [0.613118  , 0.386882  ],\n",
       "       [0.10543924, 0.89456076],\n",
       "       [0.61839164, 0.38160836],\n",
       "       [0.40897862, 0.59102138],\n",
       "       [0.71431772, 0.28568228],\n",
       "       [0.62182098, 0.37817902],\n",
       "       [0.65799105, 0.34200895],\n",
       "       [0.61277706, 0.38722294],\n",
       "       [0.6590894 , 0.3409106 ],\n",
       "       [0.33489838, 0.66510162],\n",
       "       [0.35651756, 0.64348244],\n",
       "       [0.49969552, 0.50030448],\n",
       "       [0.25928121, 0.74071879],\n",
       "       [0.19195958, 0.80804042],\n",
       "       [0.10048359, 0.89951641],\n",
       "       [0.37545518, 0.62454482],\n",
       "       [0.08787725, 0.91212275],\n",
       "       [0.03714086, 0.96285914],\n",
       "       [0.03814711, 0.96185289],\n",
       "       [0.4627366 , 0.5372634 ],\n",
       "       [0.23620869, 0.76379131],\n",
       "       [0.30420294, 0.69579706],\n",
       "       [0.0350234 , 0.9649766 ],\n",
       "       [0.13845334, 0.86154666],\n",
       "       [0.02600525, 0.97399475],\n",
       "       [0.0040018 , 0.9959982 ],\n",
       "       [0.00243055, 0.99756945]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict_proba(teste[attrs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Análise discriminante quadrática\n",
    "\n",
    "Se abrirmos mão da suposição de que a matriz de variância e covariância é a mesma em todas as classes e tomarmos uma para cada classe, a função discriminante vista acima se transforma em uma função quadrática. A vantagem desse novo modelo é que a fronteira entre classes passa a ser não-linear (no caso mais simples uma parábola ao invés de uma reta). Por outro lado, o modelo contém muito mais parâmetros uma vez que cada matriz de covariância contribui com a estimativa de $p(p+1)/2$ parâmetros, levando a um total de $Kp(p+1)/2$ parâmetros. Essa flexibilidade pode contribuir para o aumento da variância. A escolha de qual modelo é melhor depende sempre do conjunto de dados, mas via de regra o modelo linear deve ser a escolha quando o número de exemplos no conjunto de treinamento é pequeno.\n",
    "\n",
    "Para exemplificar, vamos treinar o modelo quadrático com os mesmos dados e ver se há uma melhora na acurácia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QuadraticDiscriminantAnalysis(priors=None, reg_param=0.0,\n",
       "               store_covariance=False, store_covariances=None, tol=0.0001)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis\n",
    "model = QuadraticDiscriminantAnalysis()\n",
    "model.fit(treino[attrs],treino.prognostico.map(lambda x: x=='mau' and 1 or 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>21</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0             21    9\n",
       "1              4   33"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(treino[attrs]),treino.prognostico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>9</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0              9    4\n",
       "1              1   16"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(teste[attrs]),teste.prognostico)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Bayes\n",
    "\n",
    "O classificador Naive Bayes é muito similar à análise de discriminante. Nas funções discriminantes acima, nosso maior interesse foi em avaliar o valor máximo de $P(X|Y=k)P(Y=k)$. Isso pois o denominador é independente da classe. O numerador, por outro lado, é equivalente a probabilidade conjunta de $X$ e $Y$. Assim, usando as propriedades de probabilidades condicionais e conjuntas, temos que\n",
    "\n",
    "\\begin{align}\n",
    "P(X|Y=k)P(Y=k) & = P(X,Y=k) = P(X_1,X_2,\\dots,X_p,Y=k) \\\\\n",
    "& = P(X_1~|~X_2,\\dots,X_p,Y=k)P(X_2,\\dots,X_p,Y=k) \\\\\n",
    "                                                     & = P(X_1~|~X_2,\\dots,X_p,Y=k)P(X_2|X_3,\\dots,X_p,Y=k)P(X_1~|~X_2,\\dots,X_p,Y=k)P(X_3,\\dots,X_p,Y=k) \\\\\n",
    "                                                     & = P(X_1~|~X_2,\\dots,X_p,Y=k)\\ldots P(X_p|Y=k)P(Y=k)\n",
    "\\end{align}\n",
    "\n",
    "O classificador Naive Bayes assume que os atributos são independentes entre si. Logo, o modelo pode ser reescrito como\n",
    "\n",
    "$$\n",
    "P(X|Y=k)P(Y=k) = P(X,Y=k) = P(Y=k)\\prod P(X_i~|~Y=k) \n",
    "$$\n",
    "\n",
    "Portanto, temos que a probabilidade condicional de $Y$ dado $X$ é\n",
    "\n",
    "$$\n",
    "P(Y=k~|~X) = P(Y=k)\\prod P(X_i~|~Y=k) \n",
    "$$\n",
    "\n",
    "Dessa forma, basta usarmos funções de densidade para estimarmos $P(X_i~|~Y=k)$ para ajustarmos o modelo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussianNB(priors=None)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.naive_bayes import GaussianNB \n",
    "model = GaussianNB()\n",
    "model.fit(treino[attrs],treino.prognostico.map(lambda x: x=='mau' and 1 or 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0             20   10\n",
       "1              5   32"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(treino[attrs]),treino.prognostico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prognostico</th>\n",
       "      <th>bom</th>\n",
       "      <th>mau</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prognostico  bom  mau\n",
       "row_0                \n",
       "0              8    4\n",
       "1              2   16"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(model.predict(teste[attrs]),teste.prognostico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9.53588905e-01, 4.64110947e-02],\n",
       "       [9.95408406e-01, 4.59159404e-03],\n",
       "       [9.76165352e-01, 2.38346477e-02],\n",
       "       [7.37601945e-01, 2.62398055e-01],\n",
       "       [1.79112565e-03, 9.98208874e-01],\n",
       "       [8.97149447e-01, 1.02850553e-01],\n",
       "       [7.18938180e-01, 2.81061820e-01],\n",
       "       [9.60335980e-01, 3.96640200e-02],\n",
       "       [8.88610701e-01, 1.11389299e-01],\n",
       "       [9.46200543e-01, 5.37994570e-02],\n",
       "       [5.72224579e-01, 4.27775421e-01],\n",
       "       [6.60351546e-01, 3.39648454e-01],\n",
       "       [1.32114260e-02, 9.86788574e-01],\n",
       "       [2.76382583e-01, 7.23617417e-01],\n",
       "       [5.27991561e-01, 4.72008439e-01],\n",
       "       [4.42831638e-01, 5.57168362e-01],\n",
       "       [2.53924453e-01, 7.46075547e-01],\n",
       "       [2.38180110e-04, 9.99761820e-01],\n",
       "       [2.81004421e-01, 7.18995579e-01],\n",
       "       [5.51875580e-02, 9.44812442e-01],\n",
       "       [2.65690524e-05, 9.99973431e-01],\n",
       "       [6.94959536e-06, 9.99993050e-01],\n",
       "       [3.79873582e-01, 6.20126418e-01],\n",
       "       [1.91444839e-01, 8.08555161e-01],\n",
       "       [8.04275552e-02, 9.19572445e-01],\n",
       "       [1.10724257e-03, 9.98892757e-01],\n",
       "       [9.95984194e-03, 9.90040158e-01],\n",
       "       [1.96499576e-05, 9.99980350e-01],\n",
       "       [1.29463078e-07, 9.99999871e-01],\n",
       "       [2.21829688e-09, 9.99999998e-01]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict_proba(teste[attrs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gerando conjuntos de treinamento e teste com sklearn\n",
    "\n",
    "Uma funcionalidade bastante interessante em `sklearn` é a geração de conjuntos de treinamento e teste. Obviamente, você pode separá-los manualmente, mas o pacote auxilia nessa tarefa, permitindo inclusive a amostragem estratificada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(67, 5) 41 26\n",
      "(30, 5) 21 9\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "dados = prostate.drop('train',axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    dados[attrs], dados.prognostico.map(lambda x: x=='mau' and 1 or 0), \n",
    "    test_size=0.3)\n",
    "\n",
    "print(X_train.shape, np.sum(y_train), y_train.shape[0]-np.sum(y_train))\n",
    "print(X_test.shape, np.sum(y_test), y_test.shape[0]-np.sum(y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validação cruzada com sklearn\n",
    "\n",
    "Outra funcionalidade de `sklearn` é a avaliação do modelo usando validação cruzada. O pacote fornece mecanismos para auxiliar na partição automática do conjunto de dados e na avaliação do modelo em cada partição. Em seguida, a função retorna os resultados obtidos para avaliação do usuário."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.5       ,  0.75      ,  0.78947368,  0.89473684,  0.63157895])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "model = GaussianNB()\n",
    "scores = cross_val_score(model, dados[attrs], \n",
    "                         dados.prognostico.map(lambda x: x=='mau' and 1 or 0), \n",
    "                         cv=5)\n",
    "scores "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por padrão, a função `cross_val_score` avalia a acurácia em cada fold. Contudo, o pacote permite que outras funções de avaliação sejam usada. Basta indicar uma das funções disponíveis em http://scikit-learn.org/stable/modules/model_evaluation.html#scoring-parameter como parâmetro. No exemplo a seguir, usamos a sensitividade como função de avaliação."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.53846154,  0.84615385,  0.83333333,  0.91666667,  0.83333333])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "scores = cross_val_score(model, dados[attrs], \n",
    "                         dados.prognostico.map(lambda x: x=='mau' and 1 or 0),\n",
    "                         scoring='recall',\n",
    "                         cv=StratifiedKFold(5,shuffle=True,random_state=18))\n",
    "scores "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}