SU-2016-LAB5-0036478740.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sveučilište u Zagrebu<br>\n",
    "Fakultet elektrotehnike i računarstva\n",
    "\n",
    "## Strojno učenje 2016./2017.\n",
    "\n",
    "<a href=\"http://www.fer.unizg.hr/predmet/su\">http://www.fer.unizg.hr/predmet/su</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laboratorijska vježba 5: Bayesove mreže i naivni Bayesov klasifikator\n",
    "\n",
    "(c) 2015-2016 Jan Šnajder, Domagoj Alagić\n",
    "\n",
    "<i>Verzija: 0.1</i> <br/>\n",
    "<i>Zadnji put ažurirano: 27. prosinca 2016.</i> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Objavljeno: **27. prosinca 2016.**<br>\n",
    "Rok za predaju: U terminu vježbe u tjednu od **2. siječnja 2017.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Upute\n",
    "\n",
    "Peta laboratorijska vježba sastoji se od tri zadatka. Kako bi kvalitetnije, ali i na manje zamoran način usvojili gradivo ovog kolegija, potrudili smo se uključiti tri vrste zadataka: **1)** implementacija manjih algoritama, modela ili postupaka; **2)** eksperimenti s raznim modelima te njihovim hiperparametrima, te **3)** primjena modela na (stvarnim) podatcima. Ovim zadatcima pokrivamo dvije paradigme učenja: učenje izgradnjom (engl. *learning by building*) i učenje eksperimentiranjem (engl. *learning by experimenting*).\n",
    "\n",
    "U nastavku slijedite upute navedene u ćelijama s tekstom. Rješavanje vježbe svodi se na **dopunjavanje ove bilježnice**: umetanja ćelije ili više njih **ispod** teksta zadatka, pisanja odgovarajućeg kôda te evaluiranja ćelija. \n",
    "\n",
    "Osigurajte da u potpunosti **razumijete** kôd koji ste napisali. Kod predaje vježbe, morate biti u stanju na zahtjev asistenta (ili demonstratora) preinačiti i ponovno evaluirati Vaš kôd. Nadalje, morate razumjeti teorijske osnove onoga što radite, u okvirima onoga što smo obradili na predavanju. Ispod nekih zadataka možete naći i pitanja koja služe kao smjernice za bolje razumijevanje gradiva (**nemojte pisati** odgovore na pitanja u bilježnicu). Stoga se nemojte ograničiti samo na to da riješite zadatak, nego slobodno eksperimentirajte. To upravo i jest svrha ovih vježbi.\n",
    "\n",
    "Vježbe trebate raditi **samostalno**. Možete se konzultirati s drugima o načelnom načinu rješavanja, ali u konačnici morate sami odraditi vježbu. U protivnome vježba nema smisla."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%pylab inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Bayesove mreže"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ovaj zadatak bavit će se Bayesovim mrežama, jednim od poznatijih probabilističkih grafičkih modela (*probabilistic graphical models*; PGM). Za lakše eksperimentiranje koristit ćemo programski paket [pgmpy](https://github.com/pgmpy/pgmpy). Molimo Vas da provjerite imate li ovaj paket te da ga instalirate ako ga nemate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (a)\n",
    "\n",
    "Prvo ćemo pogledati udžbenički primjer s prskalicom. U ovom primjeru razmatramo Bayesovu mrežu koja modelira zavisnosti između oblačnosti (sl. varijabla $C$), kiše ($R$), prskalice ($S$) i mokre trave ($W$). U ovom primjeru također pretpostavljamo da već imamo parametre vjerojatnosnih distribucija svih čvorova. Sve varijable modelirane su Bernoullijevom razdiobom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![This](bayes-net-sprinkler.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Koristeći paket [pgmpy](https://github.com/pgmpy/pgmpy), konstruirajte Bayesovu mrežu iz gornjeg primjera. Zatim provedite **egzaktno** zaključivanje za sljedeće upite: $P(w=1)$, $P(s=1|w=1)$, $P(r=1|w=1)$, $P(c=1|s=1, r=1)$ i $P(c=1)$. Provedite zaključivanje na papiru i uvjerite se da ste ispravno konstruirali mrežu.\n",
    "\n",
    "**NB:** Bacite oko na službenu demonstraciju danu na https://github.com/pgmpy/pgmpy/blob/dev/examples/Monte%20Hall%20Problem.ipynb."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pgmpy.models import BayesianModel, NaiveBayes\n",
    "from pgmpy.factors.discrete import TabularCPD\n",
    "from pgmpy.inference import VariableElimination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set(['S', 'R', 'C', 'W'])\n",
      "True\n",
      "P(w=1):\n",
      "+-----+----------+\n",
      "| W   |   phi(W) |\n",
      "|-----+----------|\n",
      "| W_0 |   0.6471 |\n",
      "| W_1 |   0.3529 |\n",
      "+-----+----------+\n",
      "P(s=1 | w=1):\n",
      "+-----+----------+\n",
      "| S   |   phi(S) |\n",
      "|-----+----------|\n",
      "| S_0 |   0.4298 |\n",
      "| S_1 |   0.5702 |\n",
      "+-----+----------+\n",
      "P(r=1 | w=1):\n",
      "+-----+----------+\n",
      "| R   |   phi(R) |\n",
      "|-----+----------|\n",
      "| R_0 |   0.7079 |\n",
      "| R_1 |   0.2921 |\n",
      "+-----+----------+\n",
      "P(c=1 | s=1, r=1):\n",
      "+-----+----------+\n",
      "| C   |   phi(C) |\n",
      "|-----+----------|\n",
      "| C_0 |   0.4444 |\n",
      "| C_1 |   0.5556 |\n",
      "+-----+----------+\n",
      "P(c=1):\n",
      "+-----+----------+\n",
      "| C   |   phi(C) |\n",
      "|-----+----------|\n",
      "| C_0 |   0.5000 |\n",
      "| C_1 |   0.5000 |\n",
      "+-----+----------+\n"
     ]
    }
   ],
   "source": [
    "model = BayesianModel([('S', 'W'), ('R', 'W'), ('C', 'S'), ('C', 'R')])\n",
    "\n",
    "cpd_c = TabularCPD('C', 2, [[0.5, 0.5]])\n",
    "cpd_s = TabularCPD('S', 2, [[0.1, 0.5], [0.9, 0.5]],\n",
    "                  evidence=['C'], evidence_card=[2])\n",
    "cpd_r = TabularCPD('R', 2, [[0.8, 0.2], [0.2, 0.8]],\n",
    "                  evidence=['C'], evidence_card=[2])\n",
    "cpd_w = TabularCPD('W', 2, [[0.99, 0.90, 0.90, 0], \n",
    "                            [0.01, 0.10, 0.10, 1]],\n",
    "                  evidence=['S', 'R'], evidence_card=[2, 2])\n",
    "model.add_cpds(cpd_c, cpd_s, cpd_r, cpd_w)\n",
    "\n",
    "print model.active_trail_nodes('W')\n",
    "print model.check_model()\n",
    "\n",
    "infer = VariableElimination(model)\n",
    "print \"P(w=1):\\n\" + str(infer.query(['W'])['W'])\n",
    "print \"P(s=1 | w=1):\\n\" + str(infer.query(['S'], evidence={'W': 0})['S'])\n",
    "print \"P(r=1 | w=1):\\n\" + str(infer.query(['R'], evidence={'W': 0})['R'])\n",
    "print \"P(c=1 | s=1, r=1):\\n\" + str(infer.query(['C'], evidence={'S': 0, 'R': 0})['C'])\n",
    "print \"P(c=1):\\n\" + str(infer.query(['C'])['C'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** Koju zajedničku vjerojatnosnu razdiobu ova mreža modelira? Kako tu informaciju očitati iz mreže? <br>\n",
    "**Q:** U zadatku koristimo egzaktno zaključivanje. Kako ono radi? <br>\n",
    "**Q:** Ako je trava mokra ($w=1$), što je prema mreži vjerojatnije -- da je pala kiša ($r=1$) ili da je netko upalio prskalicu ($s=1$)? <br>\n",
    "**Q:** Zašto je vjerojatnost $P(c=1)$ drugačija od $P(c=1|s=1,r=1)$ ako znamo da čvorovi $S$ i $R$ nisu roditelji čvora $C$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b)\n",
    "\n",
    "**Efekt objašnjavanja** (engl. *explaining away*) zanimljiv je fenomen u kojem se događa da se dvije varijable \"natječu\" za objašnjavanje treće. S obzirom na to da i prva i druga varijabla mogu samostalno objasniti treću, ono što se događa jest to da, ako znamo da se ostvarila treća varijabla, naše vjerovanje o prvoj varijabli mijenja ovisno o tome je li se ili nije dogodila druga varijabla (i obrnuto).\n",
    "\n",
    "Ovaj fenomen može se primijetiti na gornjoj mreži. U tom se slučaju varijable prskalice ($S$) i kiše ($R$) \"natječu\" za objašnjavanje mokre trave ($W$). Vaš zadatak je pokazati da se fenomen zaista događa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(w=1 | r=1):\n",
      "+-----+----------+\n",
      "| W   |   phi(W) |\n",
      "|-----+----------|\n",
      "| W_0 |   0.9162 |\n",
      "| W_1 |   0.0838 |\n",
      "+-----+----------+\n",
      "P(w=1 | s=1):\n",
      "+-----+----------+\n",
      "| W   |   phi(W) |\n",
      "|-----+----------|\n",
      "| W_0 |   0.9270 |\n",
      "| W_1 |   0.0730 |\n",
      "+-----+----------+\n",
      "------------------------------------------------------------------------\n",
      "P(s | w=0, r=1):\n",
      "+-----+----------+\n",
      "| S   |   phi(S) |\n",
      "|-----+----------|\n",
      "| S_0 |   0.0215 |\n",
      "| S_1 |   0.9785 |\n",
      "+-----+----------+\n",
      "P(s | w=1, r=1):\n",
      "+-----+----------+\n",
      "| S   |   phi(S) |\n",
      "|-----+----------|\n",
      "| S_0 |   0.1945 |\n",
      "| S_1 |   0.8055 |\n",
      "+-----+----------+\n",
      "------------------------------------------------------------------------\n",
      "P(r | w=0, s=1):\n",
      "+-----+----------+\n",
      "| R   |   phi(R) |\n",
      "|-----+----------|\n",
      "| R_0 |   0.0411 |\n",
      "| R_1 |   0.9589 |\n",
      "+-----+----------+\n",
      "P(r | w=1, s=1):\n",
      "+-----+----------+\n",
      "| R   |   phi(R) |\n",
      "|-----+----------|\n",
      "| R_0 |   0.3204 |\n",
      "| R_1 |   0.6796 |\n",
      "+-----+----------+\n"
     ]
    }
   ],
   "source": [
    "print \"P(w=1 | r=1):\\n\" + str(infer.query(['W'], evidence={'R': 0})['W'])\n",
    "print \"P(w=1 | s=1):\\n\" + str(infer.query(['W'], evidence={'S': 0})['W'])\n",
    "print \"------------------------------------------------------------------------\"\n",
    "print \"P(s | w=0, r=1):\\n\" + str(infer.query(['S'], evidence={'R': 0, 'W':1})['S'])\n",
    "print \"P(s | w=1, r=1):\\n\" + str(infer.query(['S'], evidence={'R': 0, 'W':0})['S'])\n",
    "print \"------------------------------------------------------------------------\"\n",
    "print \"P(r | w=0, s=1):\\n\" + str(infer.query(['R'], evidence={'S': 0, 'W':1})['R'])\n",
    "print \"P(r | w=1, s=1):\\n\" + str(infer.query(['R'], evidence={'S': 0, 'W':0})['R'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** Kako biste svojim riječima opisali ovaj fenomen, koristeći se ovim primjerom?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (c) \n",
    "\n",
    "Koristeći [`BayesianModel.is_active_trail`](http://pgmpy.org/models.html#pgmpy.models.BayesianModel.BayesianModel.is_active_trail) provjerite jesu li varijable oblačnosti ($C$) i mokre trave ($W$) uvjetno nezavisne."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# Vaš kôd ovdje...\n",
    "print model.is_active_trail('C', 'W')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Što mora vrijediti kako te dvije varijable ne bi bile uvjetno nezavisne? Provjerite korištenjem iste funkcije."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "# Vaš kôd ovdje...\n",
    "model = BayesianModel([('S', 'W'), ('R', 'C')])\n",
    "\n",
    "cpd_s = TabularCPD('S', 2, [[0.5, 0.5]])\n",
    "cpd_r = TabularCPD('R', 2, [[0.5, 0.5]])\n",
    "cpd_c = TabularCPD('C', 2, [[0.8, 0.2], [0.2, 0.8]],\n",
    "                  evidence=['R'], evidence_card=[2])\n",
    "cpd_w = TabularCPD('W', 2, [[0.8, 0.2], [0.2, 0.8]],\n",
    "                  evidence=['S'], evidence_card=[2])\n",
    "model.add_cpds(cpd_c, cpd_s, cpd_r, cpd_w)\n",
    "print model.is_active_trail('C', 'W')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** Kako možemo na temelju grafa saznati koje dvije varijable su, uz neka opažanja, uvjetno nezavisne? <br>\n",
    "**Q:** Zašto bismo uopće htjeli znati koje su varijable u mreži uvjetno nezavisne?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Naivan Bayesov klasifikator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iako je naivni Bayes specijalna podvrsta Bayesovih mreža i mogao bi se modelirati paketom [pgmpy](https://github.com/pgmpy/pgmpy), ipak ćemo se vratiti na paket `sklearn` (radi jednostavnosti). Za demonstraciju rada Bayesovog klasifikatora s diskretnim značajkama treba nam skup primjera s diskretnim značajkama i oznakama klase. Jedan takav ogledni skup podataka jest *digits*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbf8d1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAFdCAYAAABGoXXzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEX9JREFUeJzt3X+s3XV9x/Hnq1jlhyI34nBOmRaUaTAgt6DOAZtUMERR\nkllF4kIIM6AzxG2JEDFkizF1RiSgGLM5QVCX9g8jOhCobHMtsoYidQzQ6AVhDqotHRCqA+xnf5zT\neL32Ws/t993vPYfnI7l/3EP5nlf643m//Z7T701rDUlS95b0PUCSJpWBlaQiBlaSihhYSSpiYCWp\niIGVpCIGVpKKGFhJKmJgJamIgZWkImMT2CTvS3Jvkp8luTXJsX1vAkhyfJJrk/w4yY4kp/W9CSDJ\nhUk2JHk0yeYkX0ny8r53ASQ5N8mmJI8MP25J8qa+d82V5ILhr+kli2DLxcMtsz/u6nvXTklemOTq\nJFuSbB/++h6zCHbdu4uftx1JLt8bzz8WgU3yDuATwMXAq4FNwA1JDu512MABwB3Ae4HFdGOH44HL\ngdcAK4ClwI1J9ut11cADwAeBY4Bp4Gbgq0le0euqWYZfwN/D4PfaYnEncAjwguHHH/U7ZyDJQcB6\n4P+AU4BXAH8FbOtz19Byfvnz9QLgjQz+nK7eG0+ecbjZS5Jbgf9orZ0//DwM/pBe1lr7u17HzZJk\nB/C21tq1fW+Za/jF6CfACa21dX3vmSvJVuCvW2ufXwRbng1sBM4DPgx8p7X2lz1vuhh4a2ut97PC\nuZKsAl7XWjux7y27k+RS4NTW2l7529yiP4NNspTBWc43dz7WBl8V1gKv62vXGDqIwVfuh/seMluS\nJUneCewPfLvvPUOfBr7WWru57yFzvGx4KeqHSa5J8uK+Bw29Bbgtyerh5ajbk5zT96i5hi05E/jc\n3nrORR9Y4GBgH2DznMc3Mzjl124Mz/gvBda11hbFdbskRyZ5jMFfK68ATm+t3dPzLIaxPxq4sO8t\nc9wKnMXgr+DnAi8FvpXkgD5HDS1jcLb/PeBk4DPAZUne3euqX3c68Fzgqr31hM/YW0+kXl0BvBJ4\nfd9DZrkHOIrBb/g/Bb6Q5IQ+I5vkRQy+EK1orT3Z145daa3dMOvTO5NsAH4ErAT6vqyyBNjQWvvw\n8PNNSY5k8IXg6v5m/Zqzgetbaw/trScchzPYLcAvGFzcn+0QYK/9RI2rJJ8CTgX+uLX2YN97dmqt\nPdVam2mtfae19iEGLyad3/OsaeD5wO1JnkzyJHAicH6SJ4Z/E1gUWmuPAN8HDu97C/AgcPecx+4G\nDu1hyy4lOZTBi71/vzefd9EHdngmsRE4aedjw9/oJwG39LVrHAzj+lbgT1pr9/e9ZzeWAM/qecNa\n4FUMLhEcNfy4DbgGOKotoleEhy/EHc4gbn1bDxwx57EjGJxhLxZnM7iseN3efNJxuURwCXBlko3A\nBuADDF4UubLPUQDDa2CHAzvPbpYlOQp4uLX2QI+7rgDOAE4DHk+y828Aj7TWft7XLoAkHwWuB+4H\nnsPghYcTGVy/601r7XHgV65RJ3kc2Npam3uGtlcl+TjwNQbR+j3gb4AngS/3uWvok8D6JBcyePvT\na4BzgD/vddXQ8ITsLODK1tqOvfrkrbWx+GDwPtP7gJ8xeLV5ed+bhrtOBHYwuIwx++Mfe961q02/\nAP5sEfyc/QMwM/y1fAi4EXhD37vm2XozcMki2PFl4L+HP2f3A18CXtr3rln7TgW+C2wH/gs4u+9N\ns7a9cfh7//C9/dxj8T5YSRpHi/4arCSNKwMrSUUMrCQVMbCSVMTASlIRAytJRUr/oUGS5zG4OcV9\nQK9vbpekjuwLvAS4obW29Tf9wOp/yXUK8MXi55CkPpzJ4B98zKs6sPcVH39irVixorNj3XHHHRx9\n9NGdHe/9739/Z8f6yEc+wkUXXdTJsTZs2NDJcXZavXo1K1eu7ORYl112WSfH2Wn79u3sv//+nRzr\nscce6+Q4T0P37e4HVAfWywILNDU11dmxli5d2unxjjzyyM6OdeCBB3Z2vIce6vbmavvttx+HHtrN\nDaGe8Yxu/6gtWbKk82NqZLvtmy9ySVIRAytJRQysJBUxsE8DXV1HrPDmN7+57wnzOvbYY/ueMK9n\nPvOZfU/Qb8HAPg0s5sCedtppfU+Y13HHHdf3hHkZ2PFgYCWpiIGVpCIGVpKKGFhJKmJgJamIgZWk\nIgsKbJL3Jbk3yc+S3Jpk8b5hUJJ6MnJgk7wD+ARwMfBqYBNwQ5KDO94mSWNtIWewHwA+21r7Qmvt\nHuBcYDtwdqfLJGnMjRTYJEuBaeCbOx9rrTVgLfC6bqdJ0ngb9Qz2YGAfYPOcxzcDL+hkkSRNCN9F\nIElFRg3sFuAXwCFzHj8E6PZ28pI05kYKbGvtSWAjcNLOx5Jk+Pkt3U6TpPG2kG/qcwlwZZKNwAYG\n7yrYH7iyw12SNPZGDmxrbfXwPa9/y+DSwB3AKa21n3Y9TpLG2YK+LWVr7Qrgio63SNJE8V0EklTE\nwEpSEQMrSUUMrCQVMbCSVMTASlIRAytJRQysJBUxsJJUxMBKUhEDK0lFDKwkFTGwklRkQXfTUr1V\nq1b1PWFey5Yt63vCLk1NTfU9YV4PP/xw3xPmtXLlyr4nzGvNmjV9T9gjnsFKUhEDK0lFDKwkFTGw\nklTEwEpSEQMrSUUMrCQVMbCSVMTASlIRAytJRQysJBUxsJJUxMBKUhEDK0lFRg5skuOTXJvkx0l2\nJDmtYpgkjbuFnMEeANwBvBdo3c6RpMkx8g23W2vfAL4BkCSdL5KkCeE1WEkqYmAlqYiBlaQiBlaS\nihhYSSoy8rsIkhwAHA7sfAfBsiRHAQ+31h7ocpwkjbORAwssB/6FwXtgG/CJ4eNXAWd3tEuSxt5C\n3gf7b3hpQZJ2y1BKUhEDK0lFDKwkFTGwklTEwEpSEQMrSUUMrCQVMbCSVMTASlIRAytJRQysJBUx\nsJJUxMBKUpGF3K5wYkxPT/c9YV7Lli3re8K8DjvssL4n7NLMzEzfE+Z100039T1hXov5z8GaNWv6\nnrBHPIOVpCIGVpKKGFhJKmJgJamIgZWkIgZWkooYWEkqYmAlqYiBlaQiBlaSihhYSSpiYCWpiIGV\npCIjBTbJhUk2JHk0yeYkX0ny8qpxkjTORj2DPR64HHgNsAJYCtyYZL+uh0nSuBvpfrCttVNnf57k\nLOAnwDSwrrtZkjT+9vQa7EFAAx7uYIskTZQFBzZJgEuBda21u7qbJEmTYU++ZcwVwCuB13e0RZIm\nyoICm+RTwKnA8a21B7udJEmTYeTADuP6VuDE1tr93U+SpMkwUmCTXAGcAZwGPJ7kkOF/eqS19vOu\nx0nSOBv1Ra5zgQOBfwX+Z9bHym5nSdL4G/V9sP7TWkn6LRlMSSpiYCWpiIGVpCIGVpKKGFhJKmJg\nJamIgZWkIgZWkooYWEkqYmAlqYiBlaQiBlaSihhYSSpiYCWpyJ58T66xNzU11feEeW3cuLHvCfOa\nmZnpe8LYWcy/nqrjGawkFTGwklTEwEpSEQMrSUUMrCQVMbCSVMTASlIRAytJRQysJBUxsJJUxMBK\nUhEDK0lFDKwkFRkpsEnOTbIpySPDj1uSvKlqnCSNs1HPYB8APggcA0wDNwNfTfKKrodJ0rgb6X6w\nrbV/nvPQRUnOA14L3N3ZKkmaAAu+4XaSJcBKYH/g250tkqQJMXJgkxzJIKj7Ao8Bp7fW7ul6mCSN\nu4W8i+Ae4CjgOOAzwBeS/EGnqyRpAox8BttaewrY+U2ZvpPkOOB84Lwuh0nSuOvifbBLgGd1cBxJ\nmigjncEm+ShwPXA/8BzgTOBE4OTup0nSeBv1EsHvAFcBvws8AnwXOLm1dnPXwyRp3I36PthzqoZI\n0qTxXgSSVMTASlIRAytJRQysJBUxsJJUxMBKUhEDK0lFDKwkFTGwklTEwEpSEQMrSUUMrCQVMbCS\nVMTASlKRBX9X2UkwNTXV94R5rV27tu8J6tBi/r22bdu2vidMLM9gJamIgZWkIgZWkooYWEkqYmAl\nqYiBlaQiBlaSihhYSSpiYCWpiIGVpCIGVpKKGFhJKmJgJanIHgU2yQVJdiS5pKtBkjQpFhzYJMcC\n7wE2dTdHkibHggKb5NnANcA5wP92ukiSJsRCz2A/DXyttXZzl2MkaZKM/B0NkrwTOBpY3v0cSZoc\nIwU2yYuAS4EVrbUnayZJ0mQY9Qx2Gng+cHuSDB/bBzghyV8Az2qttS4HStK4GjWwa4FXzXnsSuBu\nYJVxlaRfGimwrbXHgbtmP5bkcWBra+3uLodJ0rjr4l9yedYqSbsw8rsI5mqtvaGLIZI0abwXgSQV\nMbCSVMTASlIRAytJRQysJBUxsJJUxMBKUhEDK0lFDKwkFTGwklTEwEpSEQMrSUUMrCQV2eO7aY2z\nbdu29T1hXtPT031PGDtTU1N9T5jXYv71XLNmTd8TJpZnsJJUxMBKUhEDK0lFDKwkFTGwklTEwEpS\nEQMrSUUMrCQVMbCSVMTASlIRAytJRQysJBUxsJJUxMBKUpGRApvk4iQ75nzcVTVOksbZQu4Heydw\nEpDh5091N0eSJsdCAvtUa+2nnS+RpAmzkGuwL0vy4yQ/THJNkhd3vkqSJsCogb0VOAs4BTgXeCnw\nrSQHdLxLksbeSJcIWms3zPr0ziQbgB8BK4HPdzlMksbdHr1Nq7X2CPB94PBu5kjS5NijwCZ5NoO4\nPtjNHEmaHKO+D/bjSU5I8vtJ/hD4CvAk8OWSdZI0xkZ9m9aLgC8BzwN+CqwDXtta29r1MEkad6O+\nyHVG1RBJmjTei0CSihhYSSpiYCWpiIGVpCIGVpKKGFhJKmJgJamIgZWkIgZWkooYWEkqYmAlqYiB\nlaQiBlaSiizku8pOjJmZmb4nzGt6errvCfN6+9vf3veEXVqsuxa7j33sY31PmFiewUpSEQMrSUUM\nrCQVMbCSVMTASlIRAytJRQysJBUxsJJUxMBKUhEDK0lFDKwkFTGwklTEwEpSkZEDm+SFSa5OsiXJ\n9iSbkhxTMU6SxtlItytMchCwHvgmcAqwBXgZsK37aZI03ka9H+wFwP2ttXNmPfajDvdI0sQY9RLB\nW4DbkqxOsjnJ7UnO2e3/JUlPQ6MGdhlwHvA94GTgM8BlSd7d9TBJGnejXiJYAmxorX14+PmmJEcC\n5wJXd7pMksbcqGewDwJ3z3nsbuDQbuZI0uQYNbDrgSPmPHYEvtAlSb9m1MB+EnhtkguTHJbkXcA5\nwKe6nyZJ422kwLbWbgNOB84A/hP4EHB+a+2fCrZJ0lgb9UUuWmvXAdcVbJGkieK9CCSpiIGVpCIG\nVpKKGFhJKmJgJamIgZWkIgZWkooYWEkqYmAlqYiBlaQiBlaSihhYSSpiYCWpiIGVpCIj365wkszM\nzPQ9YV4XXHBB3xPmtWrVqr4n7NLGjRv7njCv5cuX9z1BPfAMVpKKGFhJKmJgJamIgZWkIgZWkooY\nWEkqYmAlqYiBlaQiBlaSihhYSSpiYCWpiIGVpCIGVpKKjBTYJPcm2bGLj8urBkrSuBr1doXLgX1m\nff4q4EZgdWeLJGlCjBTY1trW2Z8neQvww9bav3e6SpImwIKvwSZZCpwJfK67OZI0OfbkRa7TgecC\nV3W0RZImyp4E9mzg+tbaQ12NkaRJsqDvyZXkUGAF8LZu50jS5FjoGezZwGbgug63SNJEGTmwSQKc\nBVzZWtvR+SJJmhALOYNdAbwY+HzHWyRpoox8Dba1dhO/+o8NJEm74L0IJKmIgZWkIgZWkooYWEkq\nYmAlqYiBlaQiBlaSihhYSSpiYJ8GfvCDH/Q9YV7XXntt3xPmtW7dur4naMwZ2KeBxRzYr3/9631P\nmNf69ev7nqAxZ2AlqYiBlaQiBlaSiizoOxqMYN/i40+sLVu2dHasJ554otPj3XnnnZ0d69FHH+3s\neDMzM50cZ6ft27d3fkxNlN32La21smdP8i7gi2VPIEn9ObO19qXf9AOqA/s84BTgPuDnZU8kSXvP\nvsBLgBtaa1t/0w8sDawkPZ35IpckFTGwklTEwEpSEQMrSUUMrCQVMbCSVMTASlKR/wfmW6zHvgHd\n0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbce68d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "digits = load_digits()\n",
    "\n",
    "plt.gray() \n",
    "plt.matshow(digits.images[0])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### (a)\n",
    "\n",
    "Ako na ovom skupu želimo trenirati Bernoullijev naivni Bayesov klasifikator, značajke trebamo binarizirati. Proučite funkciju [`binarize`](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.binarize.html) i binarizirajte skup podataka.\n",
    "\n",
    "Na binariziranom skupu podataka možemo trenirati naš model. Naravno, da bismo mogli ispitati koliko dobro model radi, trebamo izdvojiti dio skupa za ispitivanje. Proučite funkciju [`train_test_split`](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html#sklearn.model_selection.train_test_split) te načinite podjelu skupa na 70% primjera za učenje i 30% primjera za ispitivanje."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "from sklearn import cross_validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "db = []\n",
    "for im in digits.images:\n",
    "    db.append(binarize(im))\n",
    "db = np.array(db)\n",
    "\n",
    "db = db.reshape(len(db),-1)\n",
    "\n",
    "X_train, X_test, y_train, y_test = cross_validation.train_test_split(db, digits.target, test_size=0.30, random_state=69)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b)\n",
    "\n",
    "Pročitajte opis modula [`naive_bayes`](http://scikit-learn.org/stable/modules/naive_bayes.html). Za klasifikaciju primjera s diskretnim značajkama treba nam tzv. *Bernoullijev naivni Bayes*, implementiran funkcijom [`BernoulliNB`](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.BernoulliNB.html#sklearn.naive_bayes.BernoulliNB]). Proučite tu funkciju. Trenirajte model na skupu *digits* (na skupu za učenje). Izračunajte točnost modela na skupu za učenje i na skupu za ispitivanje koristeći funkciju [`metrics.accuracy_score`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "86.9530628481\n",
      "84.0740740741\n"
     ]
    }
   ],
   "source": [
    "import sklearn.naive_bayes as nbay\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "nb = nbay.BernoulliNB()\n",
    "nb.fit(X_train, y_train)\n",
    "\n",
    "print accuracy_score(y_train, nb.predict(X_train))*100\n",
    "print accuracy_score(y_test, nb.predict(X_test))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** Je li odnos između tih dviju pogrešaka očekivan? Zašto?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (c)\n",
    "\n",
    "Trenirajte model bez zaglađivanja (`alpha=0`) te izračunajte točnost takvog modela na oba skupa podataka."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.78520286396\n",
      "10.1851851852\n"
     ]
    }
   ],
   "source": [
    "nb = nbay.BernoulliNB(alpha=0)\n",
    "nb.fit(X_train, y_train)\n",
    "\n",
    "print accuracy_score(y_train, nb.predict(X_train))*100\n",
    "print accuracy_score(y_test, nb.predict(X_test))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** U čemu je problem? Opišite svojim riječima."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Bayesov klasifikator za kontinuirane značajke"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (a)\n",
    "\n",
    "Učitajte skup podataka *Iris* i podijelite ga skup za učenje i skup za provjeru u omjeru 70:30. Koristite funkciju [`datasets.load_iris`](http://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from sklearn import cross_validation\n",
    "\n",
    "iris = load_iris()\n",
    "\n",
    "X_train, X_test, y_train, y_test = cross_validation.train_test_split(iris.data, iris.target, test_size=0.30, random_state=55)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b)\n",
    "\n",
    "Proučite razred [`GaussianNB`](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html) koji implementira naivan Bayesov klasifikator za kontinuirane značajke. Trenirajte model na skupu za učenje *Iris* i vrednujte njegovu točnost na skupu za ispitivanje."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "94.2857142857\n",
      "95.5555555556\n"
     ]
    }
   ],
   "source": [
    "nb = nbay.GaussianNB()\n",
    "nb.fit(X_train, y_train)\n",
    "\n",
    "print accuracy_score(y_train, nb.predict(X_train))*100\n",
    "print accuracy_score(y_test, nb.predict(X_test))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (c)\n",
    "\n",
    "Naivan Bayesov klasifikator u potpunosti zanemaruje korelacije između varijabli. Postavlja se pitanje bi li klasifikator na skupu *Iris* dosezao veću točnost kada bi modelirao (linearne) zavisnosti između varijabli.\n",
    "\n",
    "Implementirajte Bayesov klasifikator s punom kovarijacijskom matricom za svaku klasu (dakle: nedijeljenom i nedijagonalnom kovarijacijskom matricom). Implementacija se sastoji od dvije funkcije: `train_GB` i `predict_GB`. Funkcija `train_GB(X,y)` neka trenira model Bayesovog klasifikatora sa Gaussovim gustoćama vjerojatnosti i vraća trenirani model. Za procjenu vektora srednjih vrijednosti $\\mu_j$ i kovarijacijskih matrica $\\Sigma_j$ možete koristiti [ove](http://docs.scipy.org/doc/numpy/reference/routines.statistics.html) funkcije. Funkcija `predict_GB(model,X)` neka daje predikcije modela na svim primjerima iz `X`. Alternativno, model možete implementirati kao klasu, a funkcije `train_GB` i `predict_GB` kao njezine metode.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import multivariate_normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class BayesClassifier(object):\n",
    "    \n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def train_GB(self, X, y):\n",
    "        X0 = []\n",
    "        y0 = []\n",
    "        X1 = []\n",
    "        y1 = []\n",
    "        X2 = []\n",
    "        y2 = []\n",
    "        for i in range(0, len(X)):\n",
    "            if(y[i] == 0):\n",
    "                X0.append(X[i])\n",
    "                y0.append(y[i])\n",
    "            if(y[i] == 1):\n",
    "                X1.append(X[i])\n",
    "                y1.append(y[i])\n",
    "            if(y[i] == 2):\n",
    "                X2.append(X[i])\n",
    "                y2.append(y[i])\n",
    "        Xj = [X0, X1, X2]\n",
    "        yj = [y0, y1, y2]\n",
    "        self.py = [len(y0)*1.0/len(y), len(y1)*1.0/len(y), len(y2)*1.0/len(y)]\n",
    "        cov = []\n",
    "        u = []\n",
    "        for j in range(0, 3):\n",
    "            cov.append(np.asmatrix(np.cov(Xj[j], rowvar = False, bias = True)))\n",
    "            u.append(np.median(Xj[j], axis = 0))\n",
    "        self.cov = cov\n",
    "        self.u = u\n",
    "        return\n",
    "\n",
    "    def predict_GB(self, X):\n",
    "        ypred = []\n",
    "        for i in range(0, len(X)):\n",
    "            p = []\n",
    "            x = X[i]\n",
    "            for j in range(0, 3):\n",
    "                u = self.u[j]\n",
    "                cov = self.cov[j]\n",
    "                f1 = 1.0/((2*np.pi)**1.5)/((numpy.linalg.det(cov))**0.5)\n",
    "                f2 = np.exp(-0.5*matmul(matmul(np.transpose(x-u), numpy.linalg.inv(cov)), (x-u)))\n",
    "                pj = f1*f2\n",
    "                p.append(numpy.asscalar(np.log(pj)+np.log(self.py[j])))\n",
    "            ypred.append(np.argmax(p, axis = 0))\n",
    "        return ypred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (d)\n",
    "\n",
    "Izračunajte točnost modela na ispitnome dijelu skupa *Iris* i usporedite je s točnošću naivnog Bayesovog klasifikatora iz podzadatka (b)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99.0476190476\n",
      "95.5555555556\n"
     ]
    }
   ],
   "source": [
    "model = BayesClassifier()\n",
    "model.train_GB(X_train, y_train)\n",
    "print accuracy_score(y_train, model.predict_GB(X_train))*100\n",
    "print accuracy_score(y_test, model.predict_GB(X_test))*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}