03. Logistic regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine learning from scratch: Logistic regression\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook, we will build the algorithm for finding a logistic regression model for classification of data in two different classes. This method yields a linear decision boundary, so the data to analyze must be (approximately) linearly separable. We start by importing some scientific libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will generate (pseudo) random data divided in two classes, centered at the points $\\left[-1, 0\\right]$ and $\\left[1, 4\\right]$. These classes have a multivariate Gaussian distribution around the centers, with specified covariance matrices. We also set the number of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_observations = 5000\n",
    "x0 = np.random.multivariate_normal([-1, 0], [[1, .75],[.75, 1]], num_observations)\n",
    "x1 = np.random.multivariate_normal([1, 4], [[1, .75],[.75, 1]], num_observations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAEPCAYAAAD1QTYeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXmcVOWZ9/29q6t3aFzaJaAMjoiSiQtCKzEhQwXbEUFcMiJPtHUyk4Cd93knPqBE4tsYYZTQ0jyTPPMEzWQZwUzaNlGDRmPALqPGaBrEJQmiEgkCGtGIgGzdXdf7x3VO16nqql6gu04v1/fzuT9VZ6lTd51azq+u1YkIhmEYhmEYR0ok7AkYhmEYhjEwMFFhGIZhGEaPYKLCMAzDMIwewUSFYRiGYRg9gokKwzAMwzB6BBMVhmEYhmH0CCYqDMMwDMPoEUxUGIZhGIbRI5ioCAHn3P9yzv3BOfd759xPnHNFYc/JMAzDMI4UZxU1c4tzbgTwLPBJEdnvnGsAHhOR/8r2mPLychk1alSOZmgYhtH/Wb9+/fsiclzY8xhsRMOewCAlChQ755qBEmBHRzuPGjWKdevW5WRihmEYAwHn3J/DnsNgxNwfOUZEtgPLgK3AO8BHIvKr9P2cc7Odc+ucc+t27tyZ62kahmEYRrcxUZFjnHNHA5cBpwDDgVLn3LXp+4nI90RkgohMOO44s+AZhmEYfR8TFbnnQuAtEdkpIs3Ag8AFIc/JMAzDMI4YExW5Zysw0TlX4pxzwBRgY8hzMgzDMIwjxkRFjhGRF4CfAi8Cr6LvwfdCnZRhGIZh9AAmKkJARG4TkTNE5FMiUiUiB8Oek2EY3aC2FuLx1HXxuK43jEGMiQrDMIzuUlEBM2cmhUU8rssVFeHOyzBCxupUGIZhdJdYDBoaVEhUV8OKFboci4U9M8MIFbNUGIYxsMiVayIWU0GxeLHemqAwDBMVhmEMMHLlmojH1UJRU6O36ULGMAYhJioMwxhYBF0TCxfqbVddE121cvhCpaEBFi1KPp8JC2OQY6LCMIyBx+G6Jrpq5WhqShUqvpBpauq512AY/RDrUtoPmDBhglhDMcPoBr4YOJwgyiN5bCZqa1WUBI8Rj6sAmT//8I9rdIhzbr2ITAh7HoMNs1QYhjGwOFLXRFMTTJ2aauU4kkBPSz81BhEmKgzDGFgcqWsiGoX77oOqKrVULF9+ZCLAf/7p0+G661JjPKxgljHAMFFhGMbAwA+ynD8/KSj8i3Ys1jVXQzwOS5bAsmXw+ONqsbjpJliwINV90d201VgMvvAFWLVKj+kLCrNYGAMMK35lGMbAwHcz+FaJaFQFQkODbu9KHEPQyrFrl7pAqqqgpSX7cwUFgv9c6cTjKlKqqtQKArpsBbOMgYaI2OjjY/z48WIYRhdobBQpLxepqhJxTqSuLnV9Y2PXj1NSoscJPq6xUWTp0tRj1tR0fOz0566qEgG9NXoNYJ30gd/vwTbM/WEYRv/Hd0f4qaSrVsGFF2phqu7WqvCtDosXqzVhwQJdTo+tCKatnn125uPU1qZaP4IWiwcftLoWxsAjbFVjo/NhlgrD6ATfGlBXl2qpqKwUAbUmZGPp0lQrw9Klehx/vX+80tLU/WbPFhk2TI89bJhIWVn7xwX3T1/XXeuJ0S0wS0UowywVhmH0f2IxtSjcdJMGQj7+ONxwA6xdm8ziyGYVCKZ81tYmYzH82hJTp6rlY+5cXa6tVavF/feDiK5buBD274d58+D3v89sGbGCWcYgwIpf9QOs+JVhdIHaWr2gr1oFlZXw3HOawjlqlAqFmhp1VbS0pAZrBoXE1KkaSHnDDcnHzZuXPN6iRTBuHFxyiW6PxeBHP4KI9/9sxAh44w19rkWLQjgJho8VvwqJsE0lNjof5v4wDGnvphBJDZz03RF+IOSMGeqSmD5d3QzV1Rp8mekYvosD1GXiu06c08f5j/eXy8pEiop0/7y81MdVVen22bNzc16MjGDuj1CGuT8Mw+gfZKtMuXlz0h1x8CCcc44GUK5erS6JRx+Fk0+GBx6Aa6/NfOyJE9VCMWmSukxGj4Y1a3T/735X3RQ//CEMH66ulNGjIZHQx7a2Ql6e7n/++fClL4FzUF9vgZjGoMPqVBiG0ffI1C8D4Mor2/flAK1WuXix3p83DwoL1XXR3Ky3GzaoSyIW08dfeSXMmqX7X345HDoE554Lzzyjro41a1QoPPCAigRQEbF9O5x5Jrz4YnJO+fn6PKDrr7gCHnpIl5uarA6FMagwUWEYRt+jo+JSJ5ygAsIXCaDWCF9sRKNqscjLU3Fx6JAKhTvvVHFx5ZVqlfjRj1QonHYaHHOMxkycdJIKipEjYetWDcS86CI9ZksLjB8P69erJUJEBYUfT5Gfr881aVJqMKZhDCLM/WEYRt/Dz4yYOVPrTVxxRdIqsWKFZnTU1SXdC8GaET6trTr8LJBIBH7xC3jlFRUdzc26/fXXVVDk5cG2bbrf22/DjBkqOlpa4MABPc6bb+p2ERUYzc16rMpKfXxxMaxbZ24PY9Bi2R/9AMv+MAYtCxeqUCguhilT1D3x2c/q7W23aYZFWZlmZDz2WPKCn0ioYABdV1Sk8RGnn66pp4cOdfy8J54I772XjJsAfY5Nm3TdpEkqHKJRGDtWLSDnn68ZJPX1WthqwYL2mSZGzrDsj3AwS4VhGH2TeFwDMKuqoKAAnnhCAy+ffDIpCkTUFbF6NXziE+qCABUUY8bo/URCm3ndc4+KAh/n2j+n78p49119XDSaHBs26PPecYdaNpxT98qIESpaXnhB97nnHhUUNTXWLMwYdJioMAyj7+HHUCxapIWsFi5MxiwcOKCWgptugssugw8+UFfFtm0wZIi6JMaMUbdGfr66JfxMjFde0WOMGaOCJJ2gZQLU0vCVryQbiiUS8J//qccoKNDnbWxUoTFtGtxyi9bGWLJE4zysvbkx2Ag7p9VG58PqVBj9is7qSXT3GI2NWvehsFAkGtVbEDnzTElpzOUvH3WU3ubn6+Oqq6WtnoRzIhdcoMtdHc5pTYr8fJFIRNdFoyLjxiWf3y+3nd4szEpxhwZWpyKUYZYKwzB6lmz1JLK5AvxmYMHlaDRZvrq+Xt0eBw/qP//CQrU+vPoqlJdrBc3PfEaXTzxRW5ZHIrrPbbfB97+vcRetrbr+t79Nff5IJz+DImqZaG5WS0Venq576SW1mNx3n9a5gNT25tdd171GZoYxADBREQLOuaOccz91zr3mnNvonPt02HMyjB4jmLnRlQ6h6SIkGlXXRtTLeH/55aRLY80aLWSVSGja5/vv6z5+Oui776pbIpFQN8nChXDWWbB7t+7X2tre7ZHu8shEcB//GOPGacBoQYHO/fLL9XWuXKlBoatWadlvExTGYCJsU8lgHMC9wJe9+wXAUR3tb+4Po19SUyOddgj18d0ENTV6e8YZ2hXU7zZaXS1SUKBltkHkpJN0v7IySXFV5OeLTJyoLpDjj5duuTkyDb8Ed6aRny9SXKylwSdOFDnvvNTXkt7ZtCfcQkaXwdwfoQyzVOQY51wZ8DngBwAickhEdoU7K8PoYeJxrSdRU9Nxh1Cfpib9V794sd6OHg0ff6z/9o85Bp56Sl0Q+/Zpeum2berm8C0QPs3N8PzzsGePpoQeKb7LJBPNzWpRWbhQMz+uvjq1SNfKlfDII0krTHfdQobRD7E6FTnGOXcO8D3gj8DZwHrgayLycbbHWJ0Ko18RvLCmV8NsakqW3/ZLcYMWsnrsMS2VvX49HH98ZlFQVtZeSIRFQUEy3fS225IZIunlxeNxfd3z5yfPRbDMuLlHegWrUxEOJipyjHNuAvA88BkRecE5921gt4jUpO03G5gNMHLkyPF//vOfcz9ZwzgcMvXtiMfhrru0OuaSJVoqe/t2rTnR2qrWh+HD4bXXksLBL4Xdl4hENL7i5JM1xuPJJ6GkJJk+2hX8gl7WHr1XMVERDub+yD3bgG0i8oK3/FPg3PSdROR7IjJBRCYcd9xxOZ2gYRwR8+frP/OgyyMWU0FRU6OFoerr4Ve/0mDKlhbN2njtNd13924VGd0VFHl5PfcaMpGfn7SsvPOOumRqatRSUV/ftWN01y1kGP0MExU5RkTeBd52zp3urZqCukIMY+CQKX5gyRL9h75kiQqA5uZkvMLrr6c+fv/+7j+nX5a7tygrg9//XrNOWlr0NS5aBA8/rKJizpyOHx90Ay1alMyQMWFhDCBMVITD/wv82Dn3CnAOcGfI8zGMnsOvHhlMK73iCnV5zJ2rgZgffqj7JBI95+LorN7EkfLBB3q7das+152Br22mkt/pNDWlxlD4qbd+PQ7DGABYTEU/wAI1jT5LpviJ5cvVvH/ttWpxWLVKgxp/+Ut44AE1+1dWajxCV2pE9EVmzNDOp/PmWcBlH8ViKsLBLBWGYRw+Hbk57rtPBUU0qq6Of/93vQDPmAG/+13X/t33VR59VJuULV6smRwmKAwDMFFhGMaRkK165rhxyX0KCjRIc/VqGD9eq19efTVM6Id/IktL9TaR0PbmVVWaDmtxEYYBmKgwDONIicX033rwX3t9vVooamo0ZmLNGi1itX697nP66VowqrfjIHqSSEQLcpWXq7g4+WTt9bF4cecBl+n9TcC6lxoDkmjYEzAMo58TTJP81re0V8dzz8GsWbBhQ7Io1F//qlkfd96pLcr9FNC8vN7P3OgJEgkoKtJ+I8OGwaZNsGyZBp+OG5cMuPQLXQXx3USZCoIZxgCiH/1NMAyjz5GeJjl+vLo5xoxRa8SjjyZTR487TsVDayt89JHWfbj44r4frBm0phw4oFaKjz6CESNUSC1frkLBFw7RaHsLRHebrBlGP8UsFYZhHD7BNMnaWjj2WLU8PPccbNyY3C+RUEtFsErmgQOaEdLXM9ASCRg1CrZsSbpAjjtO02Kvv177f7z0krpCFizQQNVMFoigm6imxgSFMSAxS4VhGEdOba3+Q3/mGbVAQLIWhY/v4vC3B9f1VYqK4Iwz4K23VEgkEmqpmDBBm4U98IAGoa5aBWPHJgVFJsFg1TSNQYCJCsMwDp+gyX/JErjmGrVAZENE3SH9hQsvVIvLZZfBzp0qLD7+WC0XsZhaJn7zG5g0SQXV1KnZBYVV0zQGASYqDMNQDidDwY8VWLJEL6grVvTuHHPNo4/Cpz6lcSIzZsD998PEifo6L7ssWZNj40ZNL73vPo2xSMeqaRqDBRGx0cfH+PHjxTB6ncZGkfJyvc20nI2lS0WqqkRApKBAb/v7OOoovXVOb2fMSD0f1dV6v64u9RzV1YmUlnZ+zoxeB1gnfeD3e7ANs1QYhqEcbobCM8/oP/Tx4+HQodzMtbfZtUtLiYtoa/Ozz9bzceWVuv2731V3SEuLukB8i8PcuRprYRYIY5BiosIwjCSZClkFSXeRxOM6IhEtbFVentv59haRiLY2nzJFhYV/PmbNSo2FqKhQF4jfEh30nKXXqTCMQYKJCsMwksyZA9/+dmqGQjCuwg/MnDNH19fXawrpxRfrhfj99/tXlcxMRCLJ2hmXXKLL0Sh85zu6zupNGEZWrEtpP8C6lBo5IR6Hyy/XWhIPPaTrgsv+hdPfb/9+Xf7Wt/R23rzkscrKYPfu3M29NygtTaa8XnSRVgfdv19FRDyerDexaFG48zQyYl1Kw8GKXxmGoTQ1wcMP6/2ZM9Xc75w2/wr+E/erRz75pC7ffLP+sw/+w+/vggKguFgtLzNmaJzEsmVajru+XpuJ+dacWMwsFYbh0c/tlIZh9Bjz5ycvkH5cxb/+K9xzT2osRTyuAqSwUJd9IdHXy213RCaXzfvvw/HHJwXF3Lm6/sEHrd6EYWTBRIVhGKlkqvxYUQGXXgpf/apeRD/3ObViDBSyCaL33tM6Fb6gsHoThtEhFlPRD7CYCiNnBCs/NjUl25cvXqzb583T/h67d/evypjdxY8JGToU9uxRy813vxv2rIxuYDEV4WCWCsMwkgT/ifvpktdfr5kOoC6PDz4Y2IKiuFgFxYwZ+npnzFCLTaZKmYZhpGCiwjCMzMRiWuzp+9/Xxlrz5iUDMgcy+/drIa+f/1wF1mc+A3V1sHZt2DMzjD6PZX8YhpHEr0PhWytOP12tEh98oNsHsoWisBAOHlTR9OKLGj8yalSykJUfV2EYRlYG+F8OwzA6xc/s8AtcNTRoc7Bzz4VbbtHiVgPdOgEqKEaO1ADUM85Ql0c02nlTNcMw2hgEvxSGYXTI5s1azCoaVSvFhg3q5tiwQS0T+fkqLAY6I0fCRx/pa33rLXV5bNqk5yRYhtswjKyYqDCMwc6sWfrv/Pbb9cI6b56KCd86MWTIwHZ7gAqJbdu0LsWhQ3o+du1K1qSw4laG0SVMVBjGYCcW0zLcLS0aSxBcX1KiRaAGOvn5ap154w2oqlJRka2pmmEYWTFRYRiGErRGRKPapXPfvtCmk1MOHNDb6mr40pf0fn6+NlezapmG0WVMVBjGYCG9bTno8pw52o3Tr5BZVKRWC7+Z1kAlvSLoiBFw1VUaX5Kfr43S0ludG4bRIZZSGhLOuTxgHbBdRKaHPR9jEBBMFw1WyzzzTI0paG2FiRPbl5wONgoLCQFchvuHTWmpZnisX59ct2MH3HRT+66ss2bpOTE3iGF0ilkqwuNrwMawJ2EMMLJZI2prk30qZs6E3/9eAzIvvFDjKE4+Wf+dP/98ewtFH2gU5rLcP2xGjUoVFEVFKrJefFEraKZ3ZfVrVRiG0SEmKkLAOXcSMA34fthzMQYYvjUi2FE0mBLpdyBdtQoqK7UD51lnwWuvaZ2GPkqPdyj6wx+S988/X60xzc1w2mma8WHuDsM4LExUhMO/A/OB8P8CGgOLoDVi4cLU6piQ2oF0wwa1VKxfD8cd16djKHqtH2phoZ6j/HyYPh2+/GUVWhZHYRiHhcVU5Bjn3HTgPRFZ75yb3MF+s4HZACNHjszR7IwBgW+NWLxYxUNQUARFxlFHqQvkE5+Ad94Jd85hUVQEDzyg5cgfeSS53o87sTgKw+gW1vo8xzjnlgBVQAtQBJQBD4rItdkeY63PjW7hi4ezz4Z165JBh7W1Gjewdq1aKJYsgTFj4LnntB7FYEkfDXLBBfr6rbX5gMNan4eDuT9yjIgsEJGTRGQUMAto7EhQGEa3iMfVjL9gAdx6K4jAFVdo2+7Nm9VyMWqU9vS44AJ4+WU4+ujBKSgiERUUM2boOTEM44gxUWEYA4mmJnV7LFyoMRMPP6wBiDffDCtXairlD34Ao0fD6tXqAvnww5xOMTTbqF92fOxYOOEEzWrJy4MTT7TsDsPoIUxUhIiIPGU1KozDIlvqKGiL7kWLtObCj36koiKR0KqRn/609rbY6GUzb9+e23nTi0GXnZFIqHXm9dfhL3/R+hyJhNanMAyjRzBRYRj9kUypo5deqjEToMLi2ms1dbSlRdcVFcGvf61xFIMR393R2qoxFK+8AsuWwS9+oe4hwzCOGMv+MIz+SDB1tLpa00QXLdLgy3HjdJ/779dbEa1J8bvfwe7d2t57oJNeBfSUU7SdOaio8oMy587V27Vrk/cNwzhsLPujH2DZH0ZWFi5Mpo4OGaKWittvV+vEwYP6r3z0aA3SPO00Nf33E/xy3EdUlts5FVVBCgrgl7+0dNEBjmV/hIO5PwyjvxIsZLVihQqKJUtUROzbp0GIpaWwcycce6wKitLSsGfdjmx/a1zabVce00Ykon1M8vOTjcOGDNHX72fEWHErw+hxTFQYRn8kWMhq0SK99QXFhg2aItncrNvOOw/ef19FxjHHaJ+PPkR3RUMmy4UEH1NWptaIq67SzA4RzXLZuxeuuw6eeAKuvrp94zTDMI4Yi6kwjP5IU1Nq+e1YLFmbwjnYskVjK269VbM+hg9Xi8W77+pFt4/iCwaXtg5SLReZ1rWxe7emjc6bp8uVlbBmDYwfrxad0aPhnnt6Y/qGMegxS4Vh9Efmz28fE9DSAnfcAcXFuvzKKyooxo7Vf+l5eWq9+OCDpEugj5FpVh2tS7FQ+IwYAW+8oa9x5EjN+Kiq0g6kY8dqUKZhGL2CiQrDGChUVKgL5JFH9CLa2qqxBRs36r/zYMCiiMYb9HF80ZAuIoLr2omO7dtVcBUUwNatWpti9WoVWzt2aCEwwzB6BRMVhjFQ8F0iAI8/DpMmaVrliSfqv/SDB1Vc+DQ3hzPPbpAuGNJdIw5IRCIIkAhaX556Sl97NKquj0OH4NFHtQ+KxVIYRq9hosIwBgp+qemZMzW+YuNGTSN9913N/sjP19TSE08Md57dJD2+It3d4RIJnh42Aydoga+CAmhuppUILZLn7eQdJRazktyG0YuYqDCMgUJtLdTXq6BYskRvN2/WC+quXWqZcE5FRh+NqchExrgJgGiUXWMncoAixrGBmpJl7DznQohE2H3aubjmgziX0JTbaNTSSA0jB5ioMIyBQG0tPP20iopNm9QNEo+rC8CPpRg+PFllsp8UvZMs919jLG8ffy4//9NZREsLKLt6Kl9YPA73wvP86KrHeGjrBBJFJeSVlqh14uGH9TXX1+f6JRjGoMJEhWEMBCoq4JlnNDizvl4biT36qLoDCgt1fQiNs7oiXSRwK2nrgriiIo2dwHFq5E/8aMc/8I95D/JM7DY2cyrjWpp45NoG/nlVjMK/O5WfXheIofCFxamn9sTLMgwjC1amux9gZbqNLhGPw+WXazXNlhYVE2efrTEUq1cn98tUujokspXgbleiu6gIDhygOVpEpLWZd+UERrCD+wuquCTyODOlgco7YyxZAlOnwn33aa8wa+cxeLEy3eFgxa8MY6DQ1ATnn6/ZDqDppNFoqqCAnAuKdOHQUS+P9G1t9w8cgGgUcXkcEhjOOzBpEjOfuY95LKNgRox587TO1eOPq6BYskQf2tJisZmGkSvM/WEYfZHa2vZBhfG4rs/Gli0qKPLzdRw6pIWffIqKemWqHRG0OKSvp5Nt7VwhLS1Emz8mmgcNBdfCM8/wVPRCaiO3sPTX57cVzpw6VS0UCxZov7WKip5+VYZhZMNEhWH0RSoqNDXUFxZ+r49sV8h4XOMoioq0yNPkyRpHAWqxOOUU/befA4IFqnzSrQ/py+n7ZrOlRID81gN8oaWelVQxuWUN0UQzx55xPBs2aM2v++7TFh9+HTBrRmoYucPcH4bRF4nFNINj5kyortaeFcFeH+ncdZc2EfvKVzR9dPFiLVe9fbtmfLz1Vq9MM5MrI1OzL399tnbm6dsi3m0isF9QaOQlmrkm8hMkAfsjpVz/8lwaHkuenlWrVGCYoDCM3GKWCsPoq8RiKigWL9bbjq6QN98M77yjNRm+8x29om7f3qvTC4qATJaF9KZf6fczZX20evskAvv7j9k9Zjyg2R8OyEu08NvIJKYlHuEr5zQRi6nB5vHH9eX/7GdJj1F3PUmGYRweJioMo68Sj6uFoqZGbzsq3BSLaQDB/v06fvYzbSDWiwR7b2Rr+pVJWGRa9omQFBb7XWnb+pbiIQx7fT1/HFGJoGW5W4nw6cQz3DZjA29cMZ94HL59eZznr6xl5UrNqJ05U2NVu+NJMgzj8DFRYRh9Ef/K19AAixYlXSH+lbG2FubMSRUamzZpb49DhzSt1K+aWV7eI1PKFlAZXM4UT5Ftvb8tKDASQB6wm6GUyMeAQ4C8/Xv56/hK/m77GvJo5btUcxN34YDJq+cxP7qcv9THaXAzOXWWqgXfg9TSkjx9CxcmT6u5Rgyj57GYCsPoi/jNwfwrn3+F9As5bd6sEYn33w9XXw2nnw733qtNwwoKNEgzkdAOncEMkCMkvTto1hTQNDLFWfiPT5CMocgDtnMiw3nXEyDCIfJpjRRwzCtPkXCOg1LIIwVX8Uw0RlEr/FvrLeR997vM+ugjeChVLcRiyUXfk1RTY4LCMHoNEbHRx8f48ePFMFJobBQZNkykpEQkL09Eq0+IRCIidXUipaWp63tgJNJu09d35XGJtPXp296lXBIgrd76nRwte6LD5J3KKhGQeN4UmVbSKL+etlRmz9ZT8K2CGj1mTU2Hp6u8XHcpL9dlY2ADrJM+8Ps92Ia5PwyjPxKLaQlqSKaOglonvv51jacYMaLHni5bIGZHZbSyld3216dbLxLA8bxPqxeIuYNPcCwfEk/8PUPW/IxHz61hfP7LXH89fOEFrWb1n1+MMzuxgkfPrWHf8hU88NV4W2CmfztnTseeJMMweg4TFYbRn2lp0dtIJGXd/kMO2bo1ZdfOREA66YIgUz0JR2oGSPrxXdoAFQ+tXqxEkIi3Tx7CamZwd+GN/On4iUxPrOa/867nmjcX8dOrGrjo7it4aeIcbjg9Tuzumaz9SgNf2rqI701pYPKKmRT8Jp4SoAnZPUmGYfQsJioMoz8yZ47WpD50SPMni4tTNhcf+AiAhNOvuF/voasNz4OCoqOS2p1ZMDIHaTp+x/kcIr9dLQqAZvL5Njfy/Ofms/a9s2jJK+Lq09ZTUQErV0FBgTBiOIxraeLtZQ189YGYVtF8JEbt+AbeWd3E1Kla/KqhAe65p30MRSxmpbsNozcwUWEY/ZEdO9Tt4RwMGZLqAvFwwG4Z2u6i7dOZ1aIrRawyVb9MTzX1604ki1hFqOB3FNDMCqo5SFHbMd50o9lHMQ9xBc1r4vys8h5+9i+PUfzOn7jo2YU8lD+T/1H4MPFZ98D8+YybG6O6WotdffazsGx9jN9Oms+qVZ2X9jAMo+cxUZFjnHMnO+fizrmNzrk/OOe+FvacjH5GbS0MHw4zZmg444oVWUtwD+OjDgtT+fczWRT82/RCVJncIkErSKZja/0JrS+xgxPJI8Eh8vkn7iWPZgT4OTM4UXbwv4fchiDMop41a+DpvBj/fqCa+QcX8+yZ1Xzt4VhbTIRfyqOqCp59Vo03zz6ry52V9jAMo+cxUZF7WoB5IjIWmAj8P865T4Y8J6M/UVEBDz4IN96oV9F0Skuhupr9US0elV4CK70gVSbLgr9feplsybIfgWMEb1u9nxgHtBBhDZWMZDtrqMQhlLIPiHATdVzBz5nBo+zf28I1xQ+zmVMBeON7cb50cAWDBpALAAAgAElEQVS1hTVctHkFf6mPs2AB1NdrzMSCBdqI9fzzYe1auOEGraq5YIEFZBpGrjFRkWNE5B0RedG7vwfYCPRcmL4xMOiotrQfaXjFFfDrX6fuc9ppEIvR+v0f8VLi7JRN6WIiW/EqAvtFAve78thg/Ylm8tr230UZUVr5PI2spIrP8TSJNLkTiUCcGMvz5rNvv66bTJz/bp3JF/MaqHh8EQtPb2Daypk8fbuem4YGHR9/DGedpS3PR41KLXplAZmGkTtMVISIc24UMA54IcO22c65dc65dTt37sz11IywydaldPPm5Lr9+zVQM8gbb8DatSQSwsSEFr3Kls7p0paD+5BlXTrB9R9TlLIughDxnCfbGI4DDlLEBs6hlSgFHGQlVRwin2XcxL8mlgMwqTVOAzN5Ob+C8yNNzKSB3xRocMTxV8e47EAD13+yiQcf1GDMF16ACRNg1ixteT5/vuquigoVFBaQaRg5JOxCGYN1AEOA9cCVne1rxa8GKZkqNjU2asGriRNFolFpKzI1Y4buAyLOya5hJ7cVlAoWnGoh0mEBqkyFq/ZQ3G7/9KJV/tjB8fIqY1PWvcexcoB8uZE6mUyj/Jbz5EPK5Ebq5GaWymQaZQ8lso8iuZ0aeY9ymUyjRCL6tFVVWuSqrEzrek2fLuKcyMkn6/bKyuSpOu883cdfbmzU5alTw34zjVyDFb8K59oW9gQG4wDygSeAuV3Z30TFIKYmQ7XI6mpdV1Ymkp+vggJEiopERowQATlUUpZRPGQSBeJVsMxUKTObeAg+roVI2z77KZQ9lMpB8tu2C8jt1LQd2hcSwaebTKOsYYoIyIrja9oERSSi4qGxUYuEFhXp/TPP1O0nn5yqt0pLdX1paVJQOKe3xuDCREU4w9wfOcY554AfABtFZHnY8zH6MNm6lI4aBePHw+7dcMYZ8NRTWukpkdDcyspK8vftJkFqnER6imgwCbUr9SuCro5EYDlCgvfQpmVFHKSUj8mnmRavnNVKqqhmBZPR+d/FfJ6ifa7n2bzMImr4wnsr+LyLE4noS/rrX/WlR6Nw3XXwwAPw6qtw4omwbRtcdVWyyNWVV+rtxx/DbbfBTTdpnMXcuV14gYZhHDlhq5rBNoDPor/HrwAveeOSjh5jlopBSNB+n77s3/f/rp94ohwsKJVDpZ5/YNiwNsvDPs9ikG51eJ+j2npsZLJCZLNoJEBaQFpwcgh1vxwiTw6S187SsZ9CuZE6eY/yttt0C4VvpQhuCy5Pn55qrKmr0/vV1UnLhHO6PG2a3q+qSnqGJk0K9200wgOzVIQyzFKRY0TkWRFxInKWiJzjjcfCnpfRx8jWpfSWW+DyyzVf8p134Mwz4d13cdE8Nhz4JDJvHvzt34LTMtjFNNOKo8Xrp+FzDLtoIZJSd8Lfvo7xKcvpQZp5wJuMJkqLVxUzQr5n9whaRwTHS4xjJg3k08JMGqigfSpGBRqM6VsvniLGTBr4h6ObGD4cvv1trTtRV6dNWaur1VgTi8Ejj0BJCfz85/DYY2qVOOccrQUWjWrNiuVmDzSM3BG2qrHR+TBLxSBh6dL27TMbG3W9v336dJHSUtlfVCavV9cl/6p71gM/lqHZsxz4loVslofNjBLfZLCLIdKCHqs1bd9DRKSFSNv6VpCHmCEHPWuFv4+A/JJK2UOpHCBfPmRYRutEV0denr7EkhK1RJSVadBm8DT5loyqqtQYisbGZFCnxVQMPjBLRSjDLBWG0VfIlkZaUZHc/vzzEIsRpZWTV3wD+fhjuOQSdv3dBUQQ8r3qlHme5aAV11YNIuHZHnzLgwCnsMWrKRFlKHtpIdpmZQjicERI0EqUQ+TzR8Yyg9X49oyEZ/VYSRXnsoEaFvFffIn7uTqjdaIj8vOT91tbNT5i/3445hhdPv30pAEnGHby+ONqyfBjKHxLxrJlWhTLMIwcELaqsdH5MEvFICJTGmmm7aNGiYDUF1TJD6saZZcbJi3OS5nwRjMRaY0kLQmJNotCXru4iQPkSwLkIHlyiGjb+nWMkwMUtFkt9lIs91Il+yjyRqHcS5W0guyhVCbT2C5GorsjP7/9unHj9La4OHOYSaZlY3CDWSpCGWapMIy+RCymQQOLF2fuiBWLwdSpsGULb4+axEWJx1m5CjYPn0REEm19OhJAHgkkkeCQZ31wwDucQD6t7GEIDtpuo7SwjvE4HB9wNAco4lGmcz+zuIUlgBawAsd1rOKPfJKDFHAJj/N7PsU86mghj1nUt8VEdNdC4dPc3H7dyy9rjERBQXJdtrATq6BpGCEStqqx0fkwS8Ugwk9pqKpq/zd86dJk0EBVlXxcqlkVH1ImeyiVZi+mYQcnpsRPbGSMtIIc9GIitjKizWLRipNfUtm2771Uyd3MbouF8K0ON1LnrS+TvV6hqhupa5fFcTNLDzt+orNRWakv369BYRgdgVkqQhlOz73Rl5kwYYKsW7cu7GkYvY0XQ9E4cQGffWYJBQsXaB3qBQs4tGgJz05awOfjC2HRIuLj5rJ16hxmuftpbJ7E2NaXGcnbvMEYxvA67xaO5MSDW/mYEkrZx2+4gAv4LbsZyjB2e7EWwgGKiNJClBYADlDAQYq5goeooIlT2Uw9swBoYCZ3soAz2ATAlTyYkrXRm/hxFsXFcM01sGWLZnsYRjacc+tFZELY8xhsmPvDMPoKnj3fzZ3LTGng0KIlMHUqLbfWMFMa+JvhLRp5OHcuTU1w1p2zKCwQTvlb+ATv8kfO4ATe5cWTZnDiwa1sYgzF7KeZPM7nBdZzLkexG4BtnMQBCkngyKeF/RSyiiqKOUQ+2k+kiQqu5EFA0z7vZAHfYAn1zOIG7jkiF0d3mTxZG461tsKPfwwnn5y935phGCEStqnERufD3B+Dj8ZGkbuKNVfyruIamT07afL3M09nzxZZMLFR9lIiv6SyLVjybmbL7xgvAnKQfJlLneyjSBIgzbi2wMxm8sQP4txHoXzIMLmdGvmQYXI3s1MKUQV7cvSWe6OzMWmSppT658KCNI2OwNwfoYzQJ2Cj82GiYhDS2Ch7S/RivrekXF6sa2y7aM6erfUb8lQTyL9FVHzcS5V8SJnso9DL5MiXVpz8B9WylxLxBcRGxqTUoNhHoRwgX/ZQ0i6O4maWyu3o8YP9O3IxotFkmxN/VFWlnKIOE2WMwY2JinBG6BOw0fkwUTFwyVTvatu02fJxfplcVtYoNTUil5U1yqHSMtk2bbaUl+uF1b/ITitRAbDIsyT8ksq2jbdTIzdS51kmorKXEtlPgfhiQi0ZUVnNNNlDSZuVA0RupE72UNpWXjvXloqTTkpdzssTOffcZEluvx5Ypn5rhiEiJipCGqFPwEbnw0TFwCWTGf+H+bNlF8PkxTrN+Hi9uk52uWGybfpsqalRl8RPzl0qUyKp9SBupM7L8siXvRTLhwyTG6mTfZ6QuJcq2UNpW2bIS5wprbg28XA3s1MERHq/jiOtP9HR8IqCCkhbh9IhQ5Lr/IqYvuUi2N7cLBVGJkxUhDNCn4CNzoeJin5IZyW301aXl4tMmaIlqGfPFhUU5eWyo7JKEjh5vbpOZs9Wq8XeknKZVtIoX48slbuZ3eay8FNLVzO9Lf3zAAWyh1KvQJWTVzlDEiC/YaK8R7n8B9XSipPVTBeQNlfHGqa0Hdu/sE+mUe5mdq+mjY4Yoe3Njz5al48+WuSMMzSWYsoUPU91dXqOLKbC6AgTFeGM0Cdgo/NhoqIf0s1IQt+MH6wY+VaVrtxRWSUHh5XLXcU1cnBYuSyb3tjWhdO3HvyW82Q101KsC/9BtRwkr00I3EidHCBf/oPqtsf66x9lakpQ5ocMkw8p63UrRUFB8v7QoXo7ZozennKK3o4dq70/gm6Obmg2Y5BioiKcEfoEbHQ+TFT0U7ponw/uVlSkF9AfVjXKTlcub1WpkFgzXAMp3qqqSSljnZcnnpViWJvLw7dc7KGkywWqMrUf94VFLuMpfDfISSfpfT8YNT8/aakwa4TRFUxUhDNCn4CNzoeJin5MJ5GE6QaMurrkBf6HVY3S2Chya0mdJJyT506rkp2uXGYe1/7i7rst9lJ8WCLgZpa2238yjbKGKSL0bOZHem+PsjJpE0h5eUlh4d/m5WkmSGOjuTmMrmOiIpwR+gRsdD5MVPRTumCpSDfjNzaK3Jq/VC7Ma5TiYpHppY1ycFi5PDm9Tn49balclN8oOymXqpNSYx18t4WfOtoTIiCXNSqqqpKtyoOi44QT2rcuNzeH0RVMVIQzQp+Ajc6HiYp+SEcxFVkCAt6cvbRtlylqIJBvRJfKsumNUlenQZwTJ6a6MIJuCz9Ycy8lKfEQRyIoejOmwo+h8K0X1dVJ909JibTFjVRWmnXC6D4mKsIZ1vujH2C9P/ohtbVQUZHaZTQe11LcFRUwcyY0NFDbFKMyGmfckpnc8rcNHH+17n/rrdDSoqWphw2DDz6AvDxdF/zK3kwtTVQA2ptjJg0AzKL+iHpz+McNPnYycSpo4i7mH8YJycwJJ8Bf/gLHHw/vvQdFRXDddXD66bBoEezfr+fgjjv0tc/vuac2BjjW+yMcTFT0A0xUDEDicfZdOpNXLqjm1DUreKq6gfKrYlxyCRw4oF3PX3wRXnhBW363tHR8uFyJgN6gvBz27IFTTtHl0aPh+ee1jTnA8uXQ2AiPPtq+E7xhZMNERThYQzHDCIHaphivXFDNxDWLaRpfzdV3x7j5ZhUUeXnw/e/Dq69CZaUKimi04+Pdxfx2FomniPV5QQHw/vv6mkePhq1bYe1aWLBABcSGDWrgWbxYjTyGYfRtOvmpMgyjJ/G9InnPxBm9dgXPV9Zw3toVTC+N8cj6GEVFKixaW9Xs//TT+rjOLBX9nZYWePxxfc3f+pZ2fH/pJbjvPpg2DebODXuGhmF0BRMVhpFDKirg25fH+e/WmcwsbuCp52JMdDF+slfjIZ46oNaG/Hw4eFAf41xqHMVAIy8PDmm3dfLzNZbirLNg1SooKDBBYRj9CXN/GEYOicWgblYTX8xr4JgvxPj4Y3gyEeNqGqigCed0v+bm5GMGsqAAtcqAunguvFCDM595RpeLi8Odm2EY3cNEhWHkmFPvmc9ZX4uxalVyXZwYdZH5iOg/98HIV76iAZm+1aKgABYu1ESZeDzcuRmG0TVMVBhGjonH4dvf1vgBgNNOg9LS5L9y/5/7QKawMHk/EoGSEg3SbG5WF0hNjVoqFi3SoE0L0jSM/oGJCsPIIfE4XHopfPKTapHIz4c334QpU2DfvrBnlzv8eBFQ8bB4MTz5pJ6PJ55QMfHww+r62bTJ6lMYRn/BRIVh5JCmJr1gvvoqfPnLap2IRmH16oEfOwEwZEjyfjQKVVVqnaip0XiKa69N1qKIxVRYnHpqOHM1DKP7mKgIAefcxc65Tc65N51zt4Q9HyN3zJ+v2QyPPAIPPACXXZYalDnQ2bs3ef+SS2DlSli2DBIJGD4c7rkndf9YzKwUhtGfsJTSHOOcywP+L1AJbAOanHOrReSP4c7M6Ck6qtANsHkzzJoFRx2laZMDPWU0ExdcoMJq+XIVWZs2td/HP2cmKgyj/2CWitxzHvCmiPxJRA4B9cBlIc/J6CFqa9WsH8xYuPRS/VdeUaGC4r77dPlPf9Ltg01QRCLwqU+phWLtWl03axY8+GDynMXjeg4rKsKbp2EY3cdERe4ZAbwdWN7mrUvBOTfbObfOObdu586dOZuccWRUVGg1yAUL9KJ43XXwi19oRseGDXrxTCS0auZgEBN+3Q2fwkKNI6mvh3Hj4LHHdH0spr0+Zs5MppE2NFivD8Pob5ioyD0uw7p2lxcR+Z6ITBCRCccdd1wOpmX0BP7FcckSGDtW3RvXXqulp2+6CX70IxUTkUhSVASDFwc648Zp1ktrqwqLILGYNlJbvFhvTVAYRv/DREXu2QacHFg+CdgR0lyMXiAWg6lTtSrkpEna02LcOBUXq1bBmDFqrfAJBi8ONERg5Eg44wyYPl27rl5zjY70rI54HFas0EyQFSus4JVh9EcsUDP3NAGnOedOAbYDs4AvhjsloydZvlzjJqqqVFAsWABXXKEX2BNP1HRSgLIy2L073Ln2NqWlWtRq+vRkYObatUm3h48fQ+G7PGIxc4EYRn/ELBU5RkRagP8JPAFsBBpE5A/hzsroKeJx/ae9bJmmSzY0wO23axfOSZPggw+S++7ZE948c8XHH2tBK78p2Ny57QUFaJZHUED4biSrpGkY/QsngyFarJ8zYcIEWbduXdjTMLqAn07a1JRMK50zR1uYb96s/SyOPVb/vQ8WfCuFYeQS59x6EZkQ9jwGG2apMIxeoKIimVY6a5ZaKJqbtQPnYBIUMDiyXAzDUExUGMYRUlur1og5c5I1KjZsgIkT4fOf17F7t8YXFBSEPdue4+ijO97unGZ6/OIX8NWvpm6Lx/W8GYYxsLBATcM4Qioq4I47tPZEfr7Wppg3T+/7HDyoGSC/+1148+xpPvww+7byco0ZiUa1lfl//idcdZW6g4JBmYZhDCxMVBjGERKLqYtj5UoNTFyxQv+h+z09ysvh/ffh+efDnWdP4deZ6IjmZrjzztTy2zNnav2JFSssq8MwBiomKgyjB5g1S9NIffyLbl6eukEeeyy1NkV/xbnMgsIXTgBFRZrpsmRJqng44QQtbFVTY4LCMAYqFlNhGD1AfX1m0TBmDDz66MAQFKBBl+mlt0EFxZgxmulx1llqlVmwIJkSaoWtDGNwYKLCMI6QeFzLbx84kFznX3g3bgxnTr1JMJvjlFOS97du1ToUL7ygFoqWFu0wGoyhWLQo2ePDhIVhDDxMVBjGEdLUBBddpP08fAZiGmVeXuprHD8etmyBGTPU5XHhhUmxEIslW5ZbYSvDGDxY8at+gBW/6vvE4yosRDoPYuyvRCLqxnEObrhBBcWFFya7sra0JAt/+YLCMMLCil+FgwVqGkYP0NQEX/mKxgsMJJxLWl0SCe1XcuhQanrolVcmXR1gQZiGMZgxUWEYPcDmzRpXUVSkF92BFJjpWyhAa0/ccIO6ORYsaJ/hYRjG4MZiKgyjB3j6aa3N0No6cASFTyIBQ4fq/ZNOggce0NbuCxeaoDAMIxUTFYZxGNTWpmYvfO5zUFiobgDQoMaBgnNqoSgs1HLjU6fCqlWa6WGCwjCMICYqDMMjXShA9h4VfsOwOXOSTcOi0WT8QWurFoQaCPgukMJC+OIX4cEHrd6EYRiZMVFhGB7BzqKQrK9QUdF+37vu0mDFlSvh4ovhiiugpCR1H7/C5EAgkYDJk+Huu7XWhNWbMAwjEyYqDMPDr58wc6bGC/gFm4Imft+aceGFeoE99lgNzPzoI9i5M7y59yaRiAagbtgAy5YlXTxWb8IwjHRMVBhGgFhMm14tXqy36TEDvjVj3Dg47zzYvj2ceeaCwkIVE2PGaOfV/fv1dQdrUASLXBmGYZioMIwA2XpU+BYK/9/55ZfDiy+GO9fewK+YGYlAZaU2Qnv/fY0ZMauEYRidYaLCMDw66lERDMwEGD062dp8IHH66ZCfrxaKZ57Rdb6YMKuEYRidYcWvDMPD71Hh/xsPxgxUVGjlyPp6uPdeOHhQ9wkWhhoIbNwIdXXq5qivT5bcttRRwzC6gokKY1BTW6uCIfgvfMMGuPRSeOSR5MXUt2Ds3681GiIRmDJFrRiJBIwYMTDiK6JRFRSxmAkJwzC6j7k/jEFNpjTSJUvU/ZGeBVJfr7dVVRrEuGaN1nCYOFFFhU+kn32r/DbtQ4ZoWmx9fbjzMQyj/9LPfv4Mo2fJlkY6d24yC2TqVN23vl5jDYYM0XiKSCRZ5GrPnuQx+4M7JC8vWfVTRIMyi4rgttu0uJXVnjAM43Cw1uf9AGt93vssXKgCoqZGrRR+0ObJJ8NLL8G0aTB8uAYvbtyo9Sk++ABGjoQdO5K1G/o6kYgGm/7gBzpnv5X5k0/qdr9RWLDrqGH0R6z1eTiYpcIY9KSnkS5fnrRY1NVBcTH84he638aN+pgPPtD6DVu3JgWF70boS5x0UvK+X3fi7bc1w8M5deUMHaoVQUFfswkKwzAOF7NU9APMUtF7BNNIYzFdnj5drRYtLckS3RdfrJUz+yt+RscVV8C+fVBQkAxEjce17sasWXDPPWHP1DB6BrNUhINZKnKIc+4u59xrzrlXnHMPOeeOCntOgx0/jdTPdIjF4NFHk4Li0ku11Xc0kCfVHwIxjz46ed852LRJX9tDD6m4uOaa1Nf88MNw6qnhzNUwjIGDWSpyiHPuIqBRRFqcc0sBROTrnT3OLBXhUFsLv/kNrF6tosI5DczsD4GYPuPHw/r1UFqamiJrGAMds1SEQz/4zzVwEJFfiYgf0vc8cFJH+xvhUlEBTzyhWRItLRqf0FcERTZrSVlZ6vJrr6nrIxq1VFHDMHofExXh8c/A42FPwshOLKaplokEDBsGb72l1oqxYzXocexY3S9YoyIXBPtzpLN7t4ogP2h03z51dzz0kLk3DMPofUxU9DDOubXOud9nGJcF9rkVaAF+3MFxZjvn1jnn1u0cqD21+zjxODz/vLY5/+gjvVCLaIOtf/gHbXVeWdm+kmZJSe/OK5HQMXSoLkciOnwhccwxapmIRHS+y5erQKqoUJeOYRhGb2GioocRkQtF5FMZxs8BnHPXA9OBa6SDgBYR+Z6ITBCRCccdd1yupm8EaGqCq66CtWvh3HM1tbSoSNNJV6/WlFK/vgOkWgd6m5EjYe9euOACFRinn54saPXBB5rdMXSoFvBqbEymyfrZLIZhGL2BiYoc4py7GPg6MENEcnDpMY6EigptHrZsmY78/KSVoKQEnnsumRXiWwvGjOn9eUUiWmvik5/UOVRWqtXkV7+CyZN1fnPnqsvjgQfgC1/Q4l7BLBfDMIzewERFbvkPYCiwxjn3knPu7rAnZGSnqUnTS+fO1YvxwoVanvuYY9QaMXSo1q4YOlQv5KeeCq+/3vPCIhJRi0NREZx/vh7/qKPg1Ve14ueaNXD22doI7eWXk0W8NmzQ9atWJV+DYRhGb2KiIoeIyGgROVlEzvHGDWHPyVBqa9v3u6ioSLZB9xuNXXst/PWvmgmyZ49aLPbs0f4fr7+uFSq7Iyr8YMshQ9pvKyrSWz+V9bHHVCRs3Qq7dsGZZ6rForJSYz/mzdMS24sW6e1NN+l6X2RYPw/DMHodEbHRx8f48ePF6F0aG0XKy/U20/LSpSJ1dbquslIERI45JvX2pJNEyspEnNNlf0SjyfvBbXl5ejtqVOr+wf0uuEBk2DCR/HyR0aP11jmR6mp9rokTdXniRJHSUl1XU6OPKS0VmT078+sxjIEOsE76wO/3YBtmqTAMsncr9V0GmzerBeDKKzUT5IIL1GJx9NF6O2wYbNumrhAJhN9ecEFqT5DgttZWjdPYsiV1Ln6WCWjMxN//ve775ptah2LZMl0+dEgDSZct03k88giMHq0lxidM0GW/7Lb/+nzLi2EYRm9gosIwPGKxZLvz6ur2MQgimmXx+9/rxb6oCE44QW9379bqldu36wXeFxLPPadxGEE+8Ynk/eC2SESzOkRSa1CsXq0xG5WVKmB++Uv48Y/hjjvUJTNunLpGNmzQUVWlsRWZXp81CjMMo1cJ21Rio/Nh7o/c4LsIamoyuwr87cOHi4BIVVVyv+pqdU1MmqTbCgtFzjxT2twZkYjItGnJdSUlkuLu8F0h/qirS7pZQMT/CPjrKitT51RVpW6QurrU9ebuMAYrmPsjlGGWCsMgtVup7+a4/PL2wY1nnw07dsCoUZpVUV2t6x94AL78ZQ2MrKyEgwc1O8MnkYAPP9R1Q4em1rIINisDffyiRfDrXyfXvfgiXHSR1syorFQLSDyetK6sWqVBpHPn6v7m7jAMIxTCVjU2Oh9mqeh9li5N/Vff2KjBjsFAx7IyXedbBQoLdd3s2Wqp8AMoS0ulXeClPwoKJCUQ84QT2u+Tl5cM7iwsVOtDfr5ktFD4waPZrCuGMVjBLBWhjNAnYKPzYaIiHILuEF9Q+BfxurpkhkVxsUhRUdL1cN55qVkefsaGv+y7PsaMkRSBMW2aHiP4uLIyXVdYqNkfJSVJ4VBXl7psLg/DSGKiIpxh7g/DSMOvWREM3Dz1VLj6au1W2tCQrFg5cSLs3w+f+Yxui8fhW99K1pjIy9NgTPGyOSIRdX2ceabWtfAzPSortdBWkM99Dm67TetN3HknvPGG7jNzpj5PS4su+wGl5vIwDCNsTFQYRhoVFXrhXr5ci0ZVVcFLL2nmx/z5qVkhfgXLl1/W2Aj/cQUF+rjW1uS+paUaW5GXp7EV+fkqKMaP11iJ666DW29VQTJlih5z0yZNGW1p0WMEhUP6XPztluFhGEZohG0qsdH5MPdH75AeRyGiy36hK+eSGR6+2yObq8Ff9mMrqqrUdeHHQvguDt/l4Y8ZM3Tb+PGS4vJobDR3hmEcCZj7I5Rhlgpj0OJbJPwMDz8DpKJCLQPXXpvM8Jg7N9W10NSUWhzLtyBs2ZJ83MGD6u6YMUOtE5GIujwiERg+PGm5OO88WL9eS3U3N6vLIxYzd4ZhGP2QsFWNjc6HWSp6j2y1KTqrWdHR8UpLRc49V4M3q6uTlo6JE6WtbHdZmW4D3a+wUNrVvjAM4/DBLBWhDLNUGIOaTFU002tW+OW7O2vI5T9u0SJt+nXHHbBypZbSXrIEzjoL6uq0Cdno0XD33XDaaRp3UViosRmPP67NwMw6YRhGfyTa+S6GMXCJxzUY0+/kGYtld200NXXcPjz4uHHjVGBceSXce29qlsauXSpiJk2CZ56B4mJ4+OGky8MXNIZhGP0Ns1QYg5ZsFomKisPLqghmYwQrXc6bl1zvi5iqKnj2WRUfBQWpz2NxFIZh9FdMVBiDlo4sEkdKugUkHk+KmAUL1HSG8IkAABX+SURBVM2xbBm8/XayK6rvXrG0UMMw+itO41mMvsyECRNk3bp1YU/D6CJBC0gwRuPKK2HWLBUtvjUkHk8u+7UnDMM4cpxz60VkQtjzGGyYqOgHmKjom9TWtneVxONw111w883t15toMIzcYaIiHMz9YRiHSbY6F+mCAsylYRjG4MBEhWF0Ab8fSDpXXqlCwo+LCMZoGIZhDDZMVBhGF8hmlZg1q32dC8MwjMGKiQrD6AJ+Zki6VQLaZ3kYhmEMVkxUGEYXSa++CYdXedMwDGOgYqLCMLpIeu2J+vreq3NhGIbRH7Ey3YbRBdJrT/jltGfNSt3P32YYhjEYMUuFYXSB3qy+aRiGMVCw4lf9ACt+ZRiG0T2s+FU4mKUiBJxzNznnxDlXHvZcDMMwDKOnMFGRY5xzJwOVwNaw52IYhmEYPYmJitzzv4H5gPmdDMMwjAGFiYoc4pybAWwXkZfDnothGIZh9DSWUtrDOOfWAidm2HQr8A3goi4eZzYwG2DkyJE9Nj/DMAzD6C0s+yNHOOfOBJ4E9nmrTgJ2AOeJyLudPHYn8OfenWFWyoH3Q3ruw8Xm3Pv0t/mCzTkX9KX5/o2IHBf2JAYbJipCwjm3BZggIn3lC5gR59y6/paWZXPuffrbfMHmnAv623yNnsdiKgzDMAzD6BEspiIkRGRU2HMwDMMwjJ7ELBVGZ3wv7AkcBjbn3qe/zRdszrmgv83X6GEspsIwDMMwjB7BLBWGYRiGYfQIJiqMFJxz9zvnXvLGFufcS1n22+Kce9XbL9RuZ865bzrntgfmfUmW/S52zm1yzr3pnLsl1/MMzOMu59xrzrlXnHMPOeeOyrJf6Oe4s3PmnCv0PjNvOudecM6Nyv0sU+ZzsnMu7pzb6Jz7g3Puaxn2meyc+yjweVkYxlzT5tThe+2U73jn+RXn3LlhzNOby+mBc/eSc263c+7GtH363Dk2coMFahopiMjV/n3nXB3wUQe7x/pQSuz/FpFl2TY65/KA/4v2XdkGNDnnVovIH3M1wQBrgAUi0uKcWwosAL6eZd/QznEXz9m/AB+KyGjn3CxgKXB1+6PljBZgnoi86JwbCqx3zq3J8D4/IyLTQ5hfR3T0Xk8FTvPG+cAK7zbniMgm4Bxo+4xsBx7KsGtfPMdGL2OWCiMjzjkHzAR+EvZceojzgDdF5E8icgioBy4LYyIi8isRafEWn0cLofVFunLOLgPu9e7/FJjifXZCQUTeEZEXvft7gI3AiLDm04NcBqwU5XngKOfcJ8KeFDAF2CwiYRXnM/oYJiqMbEwC/iIib2TZLsCvnHPrvZLiYfM/PbPwD51zR2fYPgJ4O7C8jb5xsfln4PEs28I+x105Z237eELpI+DYnMyuEzxXzDjghQybP+2ce9k597hz7u9yOrHMdPZe99XP7yyy//Hoa+fYyAHm/hiEdNSfRER+7t3/H3RspfiMiOxwzh0PrHHOvSYiT/f0XH066amyAliM/jAvBurQi3XKITI8ttdSn7pyjp1zt6Lm+h9nOUxOz3EGunLOcnpeu4pzbgjwM+BGEdmdtvlFtITzXi/+5mHUrRAmnb3Xfe48O+cKgBmo+y6dvniOjRxgomIQIiIXdrTdORcFrgTGd3CMHd7te865h1BTea9d8Dqbs49z7j+BRzNs2gacHFj2e6/0Cl04x9cD04EpkiWvO9fnOANdOWf+Ptu8z80w4K+5mV5mnHP5qKD4sYg8mL49KDJE5DHn3Hedc+Vhxgd14b3O6ee3i0wFXhSRv6Rv6Ivn2MgN5v4wMnEh8JqIbMu00TlX6gXB4ZwrRTuv/j6H80ufT9C3fEWWuTQBpznnTvH+Yc0CVudifuk45y5GAzNniMi+LPv0hXPclXO2Grjeu/+PQGM2kZQLvHiOHwAbRWR5ln1O9OM+nHPnob+DH+Rulu3m05X3ejVwnZcFMhH4SETeyfFU08lqzexr59jIHWapMDLRzk/qnBsOfF9ELgFOAB7yfjOiwH+LyC9zPssktc65c1Bz8BZgDqTO2cu0+J/AE0Ae8EMR+UNI8/0PoBA1cwM8LyI39LVznO2cOecWAetEZDV6AV/lnHsTtVDMyuUcM/AZoAp41SXTob8BjAQQkbtR8VPtnGsB9gOzwhRCZHmvnXM3QNucHwMuAd5EOx1/KaS5AuCcK0GzguYE1gXn29fOsZEjrKKmYRiGYRg9grk/DMMwDMPoEUxUGIZhGIbRI5ioMAzDMAyjRzBRYRiGYRhGj2CiwjAMwzCMHsFEhWEYhmEYPYKJCsMwDMMweoTDFhXOuZHOub1e61ujn+Kc+4NzbnKOn1Occ6N74bjfdM7d18PHLHbOPeKc+8g590APH/tu51xNTx4zDHrx/cz5Z9NQeus97Ws45/7LOfdvWbb9k3Pu2Q4e+5Rz7su9N7uMz9nn35dORYVzbotzbr8nIPwxXES2isgQEWnNxUQ7mWOBc+6n3lyluz9EufpweCV2/W6a+5xz73rPHVoVQhH5OxF5Kqzn7wf8I1rx8FgRuaonDywiN4jI4q7s29GPX3dJ+06/6x17SE8cuyfpzmfTe01d6g8zGHDOjfJ+C/t81WTn3Kecc0845953zlk1xj6Ec+5/eb8RHzntAF3Y2WO6aqm41BMQ/ujVRjaH+UV4FrgWeLeHp9OTfAe4EZiHtoceAfx/wMVhTsrokL8BXvfaeg8kLhWRIcA5aHvwTJ0mBx09fRHuDxf1PkAz0AD8S9gTMZI45/4BuAWYAowC/ha4vdMHikiHA+2lcGGG9aPQXgtRb/kUtKveHmAt8H+B+7xtk4Ft2Y4LfBP4KXAfsBv4Mip4bgE2o41oGoBjujDfbcDkzvYL7H8H0AocAPaifRluB/6Ptz0f+Bio9ZaLvX2P9pZnAH8AdgFPAWOzPM8Y73kmdDKfLwEbvfP4J2BOYNs/Ac+m7S/AaO/+JcAfvcduB27y1pejnTt3of0ZngEiGd6H84Dfevu9452LgrTnugF4A/jQe49dlteRh/Zc2OzNZz1wcoY5DwNWAjuBP6Miy5/bN/3PUAefuV97x1/jzdf/zBV5n6cPvNfTBJyQZa5jvfdul/dezvDW3w4cQn/09gL/kuGx30Q/u/d783gROLuzY3vb/gv4t+B3BBWc73nn/0vettneHA5583jEW/91733eA2xCO552+zsN1AK/CCwXAsuArcBfgLuB4sD2m7357UBbzLe9nxme6ylgCfA74CPg5wS+x3Tw/aH9b0SD91nZ4z1mgrdtFZBAe0zsBeZ38/3f4p3LV4CDaP+N4Win053AW8C/duM97+7xzgPWob99fwGWB7ZNBJ7zXsPLBH7bvPO1GPiNN49fAeXetq3e+7LXG58+gu/oNGCDN7+3gW8GjpH1PKO/V3/yjv0WcE0nn8vRgHTh8/ttbx67vTlPSntvMn5OvO3jvPdrj/f+1eN9BzM8zz955/b/oJ/d1wh8x7zz/+Uu/lYNQ/vkvIN+Z/8NyDvC385efV+A/wbuDCxPAd7t9P3p7g9QByftt+gPUQHwWe+FdkdUNAOXo2KiGP1H/zza4rcQuAf4SRfm205UAF8EXungMW0fDm/588Cr3v0LvDf3hcC2l737Y1DBUYmKj/low5+CDM9xA7ClC/OfBpwKOODv0eZB5wY+DB2JinfwvmDA0YHHLUEvDPnemESy70vwfRiP/ohFvfd3I3Bj2nM9ChyFNmjaCVyc5XXcDLwKnO69lrNRF0L6nFeiF5qh3nO+jnfxpvMv6m+B5d7n43Pol8T/zM0BHgFK0C/peKAswzzzvffsG+hn9/PecU7PNIcMj/8m+tn9R+9YN6Ff1PwuHPu/SBUVLcAi73GXeO/90en7esunoz8kwwPn5tTOPl8Z3vOTvPfp24Ht/452xTzGe18eAZZ42y5GL3yfAkrRH57ORMX2wP4/C7xHHX5/aP8bccA7L3noZ/r5bL9TXX3/A499CW0tXoz+Bq0HFnrv29+iP8L/0Nl7fpjH+y1Q5d0fAkz07o9ALwqXeMeo9JaPC5zbzd55LPaWv5Xpu3KE39HJwJneHM7y3v/LOzrP3nu9m+Rn/RPA33XyueyqqLgWtfRGURH+LlDU2efEO/d/Bv6X9779o/c+diQqWgL7X42Ki2PSrxt0/lv1MHoNKwWOR0X2nCzP2yfeF1TEXh1YLvee/9gO358u/gDtRdXOLuDh9JOGXmBagJLA4+6je6Li6bTtG0lVhZ/wPgBZvyTeft2yVKR/OLxl3xpxLGot+YZ33CHov9fvePvVAA2Bx0XQH9B2z4/+A38+w1x3ec/1N1nm9jDwtcCHvCNRsdX7MJWl7bMIvXC3++Eni2j0tt0IPJT2XJ8NLDcAt2R57CbgsizbBP0ByUP/yX0ysG0O8FRnX9TAZ640sP2/A5+5f0b/4Z3VyXs/Cf1RigTW/QRP9afPIcPjv0nqxS2CJ+66cOz/IlVU7A9+vlGLxcT0fb3l0d72C/EuZt34vG9Bv9N7vPP5JHCUt82hF/pTA/t/GnjLu/9DvAuXtzyGzkVFcP9PohaXPDr5/tD+N2Jt2nH2Z/scd/X9Dzz2nwPL5wNb0/ZZAPyos/f8MI/3NPq7Up62z9eBVWnrngCuD5zb/y+w7avA/9/e+cdqWZZx/PMlD57lDmCKsgN4WkU5WOVWZi6K0pGgIAw1M4m1+U/WYk3sx5wztFnaD7X+sJq4Od0cOGrkTJEceNoylLVWG4TmrLOj52gHOMgBwYiu/riuh/O8D8/7vM/LXtdp3p+Njfd579/XdV/PdV/3fb9nc3GuVPS75Rxt8t09wN1V44y/vPYDV5CLcrWQQy2noiTfKBEpqtITfOExRC66Gm2vciqK6Z9j3AF8mhpOBX4m600ao33XANsmslxwh3VR7nNX1P/uqnx1z1QsN7Np8W95yfe9wD4zeyP3bLBm2c3S9+F/Dni/pP24k3EMF9BbipkdxsORC3BF7McF9Il41h9Je3HPN8v3H7wfM0uK3Ys7Rvl6ZuHe36m4MUfSYknbJe2Lfl8aaepwRaQfkNQv6cJ4/kN8BbhF0kuSvl2WWdL7JT0WB3MOAN8rqTt/ZuUN3NEqYzaulFWcyfjqIWOA8vEr0guMmtmhQt6Mh3ADvF7SkKQfSOpqUs5gyK7dNmQc190o5+Uot92y91rj2Y2m42tmL+JO31rgn5LWx59Or8tyM+vBnZlzGZfzdHx188fc3Nscz8n6VOhPK4rpu6K+duYPnKh73RVnFurKv6yNfUBv1v8Yg5totD3NZH4y5V2HO2e7Je2QtCSX76pCvvk02pG687GMOnMUSRdI2iZpRNLreNQ105fScY55eXWkHZb0G0nnttG2qvaskfTXODy4H99ayNupZnrSC7xi8YYMWulvWfp25hm4HLvwccjk+As8YlHGRJHLQTy6kZH9f6yqXZ36nYph4F2S3pl7Njv3/0O4oQJAfg11Oo1Y4fMgsDjnzEwzs24ze6VDba6qG9xxuAjfg9sRny/B9z9/F2mGcIUB/HYH3u+yNm4FZkn6aLNGxMnaX+LbSGeb2TTgccLh4MRxnNHQCbMdZrYMV9ZNeCQBMxszszVm9h5gKXCDpItLmvAzfN9wjplNwQ2fStLVYRDfxqliDx596ss9O4fx8WvoL5Dv7zBwuqTTCnkBMLOjZnarmc3Ft7CWAKtK2jAEzJaUnwv5NtThuK5HObOi3E6UnXGCjprZw2Y2Hx8/A+5su1CzfjwK8qN4tAePmMzLzbup5oc6wcc9P7fPoTXF9EejnnbmT8uuNHyoL/+y/IN4ZCZve3rM7NKyPhVk3nZ5ZvY3M7sGn7d3AhtDrwfxSEU+32lmdke749GEOnMUPAL4KL6vPxXfSlW0vek4m9mTZrYQd4J2A/fVqKsSSZ/EIzifw7cGp+FbEnXs1DAwM/Qso5X+lqUvu6hQZasG8UjFmTk5TjGzeU3qnChy2YlvvWR8GHjNzPZWNaojToWZDeAr+7VxvfNC/OWV8QLuLV4Wq4Wb8dV5FT8HbpfUByBpuqRlzRJLOlVSd3ycLKm7oAxVvIbvc+bpx4Wwy8z+RYS6cOMwEmkeAS6TdHH0aw2uPM8UKzCz53HvdL2khfLfP3gHLvCMyfi4jAD/lrQY+Gzu+z8D8ySdF31dm+v/ZEnXSppqZkfxfbNj8d0SSe+L8ciel10F7onvD4b3en3zIWvJOuC7kubI+ZCkM/IJzK8jP4LLuSdkfQO+dQa+L/0p+W+iTCV3QyGnc7dG3+eT0zlJn5H0wRjjA/iLrKzPz+IG4ZuSuuTXkZfiB7jq8hFJK2I19HVcB7Z3qOyMBh2V9AFJF4UjegR3BE72evc9wEJJ58Wq+z7gbklnRV0z5SfBweX1JUlzYxHxnRrlr8ylvw3YmJN9rflTg+L41JV/Gc8BByR9K5un8muP5+fSNJN52+VJWilpeoz9/shzDJ8HSyVdEnm6JX1a0qwafRjBD68W7VqelnM06MEj0UckfQw/o0a0vXScJZ0t6fJwjt7EV72l4x91d+P2j+hns/dDD77tOQKcIukWGlfTVfwh8q6WdIqkFfgisYqzIn2XpKvwg9ePl6SrslXD+CHaH0uaImmSpPdKWtCkzgkhF/y823Uxd0/H39sPNEl7nE7+oua1+N7rXvxk6wa80ZjZ6/h+3zp8FXIIDxdW8RPcC9siaQyfsBdUpH8eN6wz8bDPYWIVFC/bnS3qulLSqKSfxrNn8LMVWVRiF268s8+Zo7ASPx28B39hLA0npIyv4tdK78JvYbyMn96+Gt9zHQNW48Z2FFeSR3P1vYAb5afwGxjFH2b5IvAP+dbFl6NtAHMiz0F8Yt1r5ff/b4w6x/AXy4Ym/ajDXdGPLbhS34+PZ5Gv4frwEt6fh/F9e8zst9GGv+AH3R4r5P0CrhP78Jfbg7nvZuAn9A/gW2f9jDsrxwlZXQ4sxmV4L7DKzHa30ddf4zIcxWWwIlYKnSg7435grjx8ugl3Pu+Icl/Fjd9NUEvfGwgn+UH8jAP4SvBFYHvo0lP4oTHM7AncCdkaabbWqOIh3Bi9ip9KXx1ltTt/qvg+cHOMz43UlH8Z4fAsxa/b/j3atg4Ps2eUyvwky1sE7JR0ELdFnzezI2Y2CCzD5TqCr2C/QQ27HVvRtwO/jzH5eEmyunP0K8BtYYdviTwZzcZ5Eu4kDuHzc0GUU0Yfbq8znT2M2/MyngSewBeqA7hNrrXVHnq1Aj8rMYrL71ctsj2L2889+HheWbZSr2GrVuFO066oeyOF7fAcE0IuZrYZvxm2DR/rAWosIrIbAB1H0gZgt5nVWckkEv+XSFqLH5xa2Srt2xFJT+MH2Nb9r9vSKZLME4nmdCxSIen8COlMkrQI97A3dar8RCKRSCQSE5tO/trbDDyUdAYe1r/ezP7UwfITiUQikUhMYN6y7Y9EIpFIJBJvL9KfPk8kEolEItERklORSCQSiUSiIySnIpFIJBKJREdITkUikUgkEomOkJyKRCKRSCQSHSE5FYlEIpFIJDrCfwFvnrp1xDJsbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113724588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x0[:,0],x0[:,1],'x',color='blue')\n",
    "plt.plot(x1[:,0],x1[:,1],'x',color='red')\n",
    "plt.figtext(0.5, 0.01, 'Figure 1: two Gaussian clouds of points. Red points represent class 1 and blue class 0', \n",
    "            wrap=True, horizontalalignment='center', fontsize=12)\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will assign to each of these points a corresponding class, labeled by $\\{ 0,1\\}$. Then we stack the vectors points in the same vectors, so we have a single vector of training examples ($X$) and a single vector of labels ($Y$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.zeros(len(x0))\n",
    "y1 = np.ones(len(x1))\n",
    "\n",
    "X =  np.vstack((x0,x1))\n",
    "Y = np.append(y0,y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main functions of the algorithm is the sigmoid function. This is defined as on the real numbers by\n",
    "\n",
    "\\begin{align*}\n",
    "g(x) &= \\dfrac{1}{1+\\exp(-x)} ,\\\\ \n",
    "g'(x) &= g(x)(1-g(x)).\n",
    "\\end{align*}\n",
    "\n",
    "We extend its action to vectors of arbitrary length by acting componentwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(v):\n",
    "    return 1/(1 + np.exp(-v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "Suppose we count with a set of examples $X = \\{x^{(1)},\\dots,x^{(m)} \\}$ labeled with their corresponding classes $Y = \\{y^{(1)},\\dots,y^{(m)}\\}$. Our goal is finding a function $h(\\cdot,\\cdot)$ and a parameter $\\theta$ such that for a new example $x^{(m+1)}$, we have\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(y^{(m+1)}=1 \\ |\\ x^{(m+1)}) = h(\\theta,x).\n",
    "$$\n",
    "\n",
    "We will try to separate the two datasets with a linear manifold (straight line), which can be written as $\\ell(\\theta)=\\{ (x_1,x_2) : \\theta_0 + x_2\\theta_1 + x_2\\theta_2 \\}$ for some parameters $\\theta_0,\\theta_1,\\theta_2$. Then, it is reasonable to assume that the probability that a given point $(x_1,x_2)$ is at either side of the line depends linearly on the components of the point $(x_1,x_2)$. Now, trying to adjust $\\mathbb{P}(y^{(m+1)}=1 \\ |\\ \\theta, \\ x^{(m+1)})$ with a linear function of the form $\\theta_0 + x_1\\theta_1 + x_2\\theta_2 $, as this would predict negative values for a probability for vectors in a certain region, as well as values larger than $1$ for vectors in another region. Thus we scale the interval $[0,1]$ to the whole real line, so we can fit a linear function to this scaling. In general, there are many (arbitrary) ways to do that, but we choose one of the most widespread choices:\n",
    "\n",
    "$$\n",
    "h(\\theta,x)= \\mathbb{P}(y =1 \\ |\\ \\theta, \\ x)= \\dfrac{1}{1+\\exp(-(\\theta_0+x_1\\theta_1 +x_2\\theta_2))} = g(\\theta\\cdot (1,x_1,x_2)).\n",
    "$$\n",
    "\n",
    "This transformation makes the log of the odds linear. For classification problems, we will define a threshold $t_0$ such that we say that a new example $(x_1,x_2)$ belongs to the class $1$ if $P(y=1 \\ | \\ \\theta, \\ x)\\geq t_0$ and belongs to the class $0$ otherwise. In general we will take this threshold to be $t_0=1/2$. We introduce the notation $x_0=1$ for all examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we address the question of how to find such parameters. Again, there are many ways to do this. There are two main different approaches that one could follow:\n",
    "+ Minimize a cost/error function of the form $C(\\theta)$ that takes into account the deviations of the predicted class for training examples with respect to the labels of such examples.\n",
    "+ Construct parametric estimators $\\theta$ according to statistical methods of parameter estimations, such as maximum likelihood or the method of moments.\n",
    "\n",
    "In many cases these two approaches are equivalent, but we will follow the second one as it makes the mathematical foundations more clear. In particular we will follow the method of maximum likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum likelihood\n",
    "\n",
    "In our model, given our sets of training examples $X,Y$, if we fix the parameter vector $\\theta$, then for a new training example $x$, $y$ can take the value $1$ with probability $\\mathbb{P}(y=1 \\ | \\ \\theta, \\ x)=h(\\theta,x)$ and the value $0$ with probability  $1-\\mathbb{P}(y=1 | \\ \\theta, \\ x)$. We define the **likelihood function** by\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\theta ) = \\mathcal{L}(\\theta \\ | \\ X,Y) = \\prod_{k=1}^m \\mathbb{P}(y = y^{(k)} \\ | \\ \\theta, X)\n",
    "$$\n",
    "\n",
    "where the product ranges over all the training examples. With the observation above, we can write this function as\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\theta ) =  \\prod_{k=1}^m  (h(\\theta,x^{(k)}))^{y^{(k)}} (1-h(\\theta,x^{(k)}))^{1-y^{(k)}}\n",
    "$$\n",
    "\n",
    "for $t\\in\\{ 0,1\\}$. The **maximum likelihood estimate** is defined by the following relation\n",
    "\n",
    "$$\n",
    "\\theta^*\\in\\mathbb{R}^3 : \\mathcal{L}(\\theta^*  | \\ x) = \\sup_\\theta \\mathcal{L}(\\theta \\ | \\ x).\n",
    "$$\n",
    "\n",
    "The main idea of this definition is that we look for the parameter such that the probability that we observed the input data is maximum (given that it is distributed in the way we assumed). Note that the existence or uniqueness of such parameter are not obvious. If there are more than one parameter achieving the supremum, we pick any of them. The maximization of $\\mathcal{L}$ is equivalent to the minimization of $-\\log \\mathcal{L}$, which can be written as\n",
    "\n",
    "$$\n",
    "-\\log\\mathcal{L} = \\sum_{k=1}^m -y^{(k)}\\log (h(\\theta,x^{(k)})) -(1-y^{(k)})\\log(1-h(\\theta,x^{(k)})).\n",
    "$$\n",
    "\n",
    "This quantity is more convenient to work with, as dealing with products of small numbers can cause computational problems. Let\n",
    "\n",
    "$$\n",
    "C(\\theta)= \\dfrac{1}{m}(-\\log\\mathcal{L}) = \\dfrac{1}{m}\\sum_{k=1}^m -y^{(k)}\\log (h(\\theta,x^{(k)})) -(1-y^{(k)})\\log(1-h(\\theta,x^{(k)})).\n",
    "$$\n",
    "\n",
    "Many sources call this function a **cost function**, and hence our notation, in which the dependence on $X$ and $Y$ is implicit. In this sense, our approach becomes a cost function minimization. This function can be writen in vectorized to make computations more efficient:\n",
    "\n",
    "$$\n",
    "C(\\theta) = \\dfrac{1}{m}(Y^T\\log(g(X\\theta))-(1-Y^T)\\log(1-g(X\\theta))).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Properties of $\\theta^*$\n",
    "\n",
    "The estimator $\\theta^*$ given by the maximum likelihood method has several good properties that make it useful. We will state some of them without precising on the formal meaning or giving the proofs. Assuming that the dataset is linearly separable, then\n",
    "+ Consistency: as we increase the number of training examples, $\\theta^*$ converges (in a precise way) to a parameter $\\theta_0$ that correctly classifies all the observed data.\n",
    "+ Asymptotic efficiency: the fluctuations of the estimated parameter $\\theta^*$ with respect to the *true* parameter of the model $\\theta_0$ satisfy a central limit theorem, i.e., are asymptotically normally distributed. It is also possible to compute the order of such fluctuations.\n",
    "\n",
    "In practice it is not always the case that the data is linearly separable, but this method still yields a reasonable approximation. It is important to evaluate this before choosing this model to make predictions with our data, as there are other models which might yield better predictions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the regression\n",
    "\n",
    "We start by initializing the parameter $\\theta$ to a small random value, to break the symmetry of the system and avoid the algorithm getting stuck."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = 2*np.random.rand(X.shape[1]+1)-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the cost function, that computes the current total cost as a function of the parameter $\\theta$. It also appears to be a function of the dataset, but this parameter will be fixed throughout the whole notebook. This implementation of the cost function is vectorized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost_function(X, Y, theta):\n",
    "    m = X.shape[0]\n",
    "    M = np.hstack(((np.ones(X.shape[0])).reshape(-1,1),X))\n",
    "    h = sigmoid(M.dot(theta))\n",
    "    J = (1/m)*(-(Y.T).dot(np.log(h)) - ((1-Y).T).dot(np.log(1-h)))\n",
    "    return J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use gradient descent as the optimization algorithm, and hence, we will need the gradient of the cost function with respect to the parameter $\\theta$. We proceed componentwise using the chain rule: for $j=0,1,2$, we have\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\dfrac{\\partial C}{\\partial \\theta_j} &= \\dfrac{1}{m}\\sum_{k=1}^m\\dfrac{\\partial }{\\partial \\theta_j} \\left( -y^{(k)}\\log (h(\\theta,x^{(k)}))\\right) - \\dfrac{\\partial }{\\partial \\theta_j}(1-y^{(k)})\\log(1-h(\\theta,x^{(k)}))\\\\\n",
    "&=\\dfrac{1}{m}\\sum_{k=1}^m -y_jx_j\\dfrac{g'(\\theta\\cdot(1,x_1,x_2))}{h(\\theta,X)} + (1-y^{(k)})x_j\\dfrac{g'(\\theta\\cdot(1,x_1,x_2))}{1-h(\\theta,x^{(k)})} \\\\\n",
    "&= \\dfrac{1}{m}\\sum_{k=1}^m -y^{(k)}x_j \\dfrac{h(\\theta,x^{(k)})}{h(\\theta,x^{(k)})} +(1-y^{(k)})x_j\\dfrac{h(\\theta,x^{(k)})(1-h(\\theta,x^{(k)}))}{1-h(\\theta,x^{(k)})}\\\\\n",
    "&=\\dfrac{1}{m}\\sum_{k=1}^m x_j(h(\\theta,x^{(k)})-y^{(k)}).\n",
    "\\end{align*}\n",
    "$$\n",
    "Again, we write this in vectorized form:\n",
    "\n",
    "$$\n",
    "\\nabla_\\theta C = \\dfrac{1}{m}X^T(g(X\\theta)-Y).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient(X, Y, theta):\n",
    "    m = X.shape[0]\n",
    "    M = np.hstack(((np.ones(X.shape[0])).reshape(-1,1),X))\n",
    "    return (1/m)*np.array((M.T).dot(sigmoid(M.dot(theta))-Y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "\n",
    "In order to find the minimum of the cost function, we use [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent). For this, we start with a random value of $\\theta$, and then follow the direction in which the gradient points:\n",
    "\n",
    "$$\n",
    "\\theta \\mapsto \\theta -\\alpha \\nabla_\\theta C.\n",
    "$$\n",
    "\n",
    "Here $\\alpha$ is the **learning rate**, a parameter which controls the weight assigned to the gradient when substracting it from the parameter. A small learning rate is more likely to yield convergence to the optimum, but at the cost of taking more iterations. On the other hand, a large learning rate can converge faster, but it might miss the optimum. The learning rate is a [hyperparemeter](https://en.wikipedia.org/wiki/Hyperparameter) of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic_regression(X, Y, num_steps, learn_rate):\n",
    "    theta = 2*np.random.rand(X.shape[1]+1)-1\n",
    "    for i in range(0,num_steps):\n",
    "        grad = gradient(X,Y,theta)\n",
    "        theta = theta - learn_rate*np.array(grad)\n",
    "    return theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run our model and store the parameters obtained to later plot the line defined by such parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = logistic_regression(X,Y,2000,0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot now our dataset and the corresponding line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEPCAYAAABWc+9sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXt8VOW1//9+5pKEBBKi3OQmAUWsV4Qohx5appCeIohKj8ipoMd+WzT++m17vJ2i32CFUxE0nPb8+hXt9RyhbcRWLKLUgomVlmrDRXpTVAgC4gU03G9J5vn+sWZn75nMJEMyySST9X699mvP3vuZvZ/ZmXz2mvWsZy1jrUVRFEXJHHzp7oCiKIqSWlTYFUVRMgwVdkVRlAxDhV1RFCXDUGFXFEXJMFTYFUVRMgwVdkVRlAxDhV1RFCXDSImwG2P+zRjzN2PMX40xvzDG5KTivIqiKMqZY9o689QYMwj4PfApa+0JY8xK4AVr7X8nek+fPn3ssGHD2nRdRVGU7sbmzZsPWGv7ttQukKLrBYAexpg6IBfY11zjYcOGsWnTphRdWlEUpXtgjHk3mXZtdsVYa98DHgV2A+8Dh6y1v43TobnGmE3GmE379+9v62UVRVGUBLRZ2I0xhcC1QBEwEMgzxsyObWet/YG1dqy1dmzfvi3+klAURVFaSSoGTycDNdba/dbaOuAZYHwKzqsoiqK0glQI+25gnDEm1xhjgEnAGyk4r6IoitIKUuFjfw34JbAF+EvknD9o63kVRVGU1pGSOHZr7QPW2lHW2outtXOstadScV5FUTohS5ZAVVX0vqoq2a90CnTmqaIoZ0ZxMcyc6Yp7VZVsFxent19KI6mKY1cUpbsQCsHKlSLmpaWwbJlsh0Lp7pkSQS12RckkOspNEgqJqC9cKGsV9U6FCruiZBId5SapqhJLvaxM1rEPEyWtqLArSibhdZPMny/rZN0kyVr7zsNi5UpYsMC9nop7p0GFXVEyjda6SZK19qurox8WzsOkujp1n0FpE23O7tgaxo4dazUJmKK0E44gt2Zgsy3vjceSJfJg8J6jqkoeAvfe2/rzdlOMMZuttWNbaqcWu6JkEm11k1RXw5Qp0dZ+WwZfNTQyLaiwK0om0VY3SSAAK1bAnDlisS9d2jYhdq4/bRrcfHO0z18nNbUf1toOX8aMGWMVRUkhixdbW1kZva+yUvYnS2WltX36WFteLus5c6w1Rrbbeq05c6wFWXuvFXsepVmATTYJjdUJSoqSCTguD8c6DwRg0SLZhuT82l5r/+BBccfMmQP19Ymv5VjeznY8qqpg7Vo514oVsm/tWp3U1J4ko/6pXtRiV5R2wLGCYy3tM7WOKyutzc2V83jf57XKnXOWlTV/7thrx1ruyhlBkha7CruidHW8rpGyMvm3LikRcW5JeGOJdcfErr3nca41aVJi14y3b94HT16eumFagQq7onQXEvnGS0rkX7ysLPF7Y/3lixfLeZz9iYR47lxrCwrk3AUF1ubnN32ft33sPvWxtwoVdkXpTpSXi5g77pPS0ujtZFwljqjHc514Hw5z54rQ5+dLu8pKa3NyXBdLvOulYnBX0cFTRelW1NfD7NmwfDmUlMCTT8Ltt8OwYXD55RJuuHChtPMOoFZXw7x5Mvg5ZYoMbt5+u+zfutU939Kl0Ls33Hkn7NsHp07B1VfDNdfAjBmyPWiQtC8razooGm/QNhTSwdP2Ihn1T/WiFruinCEtWbyOa8SxsKdPF4t62jTXgs/NjX8Ox93i+OYdN44x8j7vLwBne/p0aX/WWbIeM8b9hZCfL/1RUg7qilGUDCKRj3ruXHGfFBSIO6S8XMQXrA0GZT16tNs2nrBPmyaiPGGCrM8/3zaJOc/OtnbQINl/xRUi3iNHyrYxsh43Tto6Pnf1n6ecZIVdXTGK0tmIl18FxOURm8cFXDcLwF13QXa2xLHX1cl661bXPTJzppxn1ixpf911cPo0XHEFbNggbpd168Dvh6efhltvlXbhMLz3HlxyCWzZIsffeguKiqCmRtps2QLXXw+rVsl2dbW6WtJFMuqf6kUtdkVphuYiSJwQQ+9gpjemPBCQ436/WNmOW8XvF8t87lxxyQSDsm/UKGvHj5f3DB4s66FDZe3zyflycuQ8Y8bI/sJC9xq5udZmZcm2E/qotBuoK0ZRujCOWE+aJK4NJ/rE8YfH+ssdwXeE3Xnt+MaDQVmPGyeC7LTx+VyRdraNcX3ozlJaKv1wzu+4YRyXTV6eHHP6qrQLKuyK0tVxxLpHD7G2CwqsnTrVjRkvKLB2yBAR4UBALGfHEvcKd26u64v3WteJlgEDXMF3ltGj5Tw5OfILwLHkR4+W41ddJYI+d64bU6+hjCknWWHX7I6K0hmpqpIQwzlzICsLXnwRTpyAl14SnziI5BoDq1fDOedAMCj7Gxpg5Eh5HQ7DF78ITzwB27e75zem6TV9ETn44AN5XyDgLlu3ynW/8x3YuFHen50tIY45OfDaa9LmiSckfLKsTFPzppNk1D/Vi1rsitIM8ab15+a61rMTiui4QKZPl+3+/aPdJI71HgzKOceNa+pGaWlxImwc63/UKDlfXp7rEiovl18SwWDz+WWUNoO6YhQlDaRihmVsfpX8fBkIDQRkDdZecomNCkl0tnv3dkU9P98VZr9fxN8ZKE12MUZcLsGg654JBFwXjFfINTVvu5OssKsrRlFSyZlWDIotIL1kibg+nMIYFRXigjl1SkIHs7Ml1PAvfxF3y/LlcPHFsj1ggKTb9ftFlh94AH70I8jPF/eMzwd//GP09X0tSIC14oKpqxP3jHPu11+HCRNkpuq4cdLWm5o3tqiG0rEko/4tLUBv4JfAm8AbwD80114tdiWjSTalrbet08bJ+eKk3L3qqmj3yUUXyXEnJNFx0fTp09Qiz8tzQxRTufh8MkkJ5BeEN2+MtZqatx2hI10xwP8AX4m8zgJ6N9dehV3JeOLFmyci9kEwapTrw3am8WdluSI+eLC0y893XSNgbVGRuEz69XNdMm1ZvNE1sUswKNE6BQXiu7/yyujPEpsRUpOApYQOE3YgH6gBTLLvUWFXMpozsdgdJk1yrdxp01wBzc219sIL3e0ePaLFPJHoOukEUmGdJzpWVhb9C6O5iVWatjclJCvsRtq2HmPM5cAPgL8DlwGbgW9Ya48les/YsWPtpk2b2nRdRemUeMvExZaNq652UwU4aQNA/OhPPQXHjkn2xaFDYffupufOz4fDhzv28yQiK8sNhXzgAbd8XmwqBG9JPudeeFMiqP/9jDDGbLbWjm2xYTLq39wCjAXqgasi298DFsZpNxfYBGwaOnRouz/ZFCUtJHI5TJnihi/OnSvhgU60SX6+OyPUSajVGRfHeh8yxP2FES9jZHOciYtKaQId6IoZAOzybE8Anm/uPeqKUTKaeOLuxKKXl4uQO64SY6Ljys8+u3Wi294PhGDQHcgNBOQhVFZ2Zil6W+OiUqLoMGGXa7EBuCDy+tvAI821V2FXMppE/mTHYneSaDl5WdJtiSeznH22DIY60TjjxrmfLRlxVx97SkhW2FMVx/6/gZ8ZY/4MXA48lKLzKkrXw1uVaP58Wc+bJ37oKVOgtlbahcMim6mgpXj0tvLxx7LevVuu9ZDnXzxeeoJYqqujfeqhkDvuoKScNg+etgYdPFUygnh506uqZDB0xQr49Kclt3lJieRRGT9e8rqUlEjOl3A4fX1vC9Onw/r1kvtdB0E7lGQHT3XmqaK0lkSzTGfNgltuEVEvKpL1yJEi6tOnw5/+lJyV21lZs0YSiy1cKBEuKuqdDhV2RWktjjvB63Jxqho9/bRY5jU1Iu4bN8KYMbK+8UYY23LEWqcjL0/W4TA884ykDygvj06JoHQKVNgVpS2EQmK1eq1Xx8e+davkU6mpgV69YPNmaXPBBZLmtr394qnE55M4+z59ROCHDJHcMAsXRv9qiUdsPhyQ7SVL2rfP3RiteaoobaGqSvzMZWXw8MNw4ADccANcc43UB62udicWOYOOPXtKMi2QiT5OfvXOTDgsedcPHICCAsnt/uijcOedMHq0OwjqTEby4ris4k3aUtqHZEJnUr1ouKOSEcSG7Dn5zi+8MDqPud/vhgk6S06OTFJqbsp+Z1hi+5eXJ+vBg924fO+9SFQ5SWPYUwJJhjuqxa4orcUbwrdkCZx9tljib7whVYgcGhrgvfdkwNSJQjt5En7zm84fGRMOw7BhsGuX647p21dCNm+5Be6+W1L4rl0r7qdFi+Jb4l6XVVmZDri2M13IyaconRQnh/qGDW55OidW3aGhQdbOce++zkpODowaJWMEffuKyOflycDvc8/JAPHkyZIT/sILXVGPJ9pel9WyZTrg2s6osCtKa3F8x4GAiNpNN4klnghrpWBFV2HyZPn1ce21sH+/iPuxY2LBh0Jiof/hDzJAvGGDTL5KJOqOT33BAjeSSMW9/UjGX5PqRX3sSqejtfnCvfnH0+0Pb4/lootkPX16dN3U6dOj0yQ4ueMdn3sq7q3SBLTmqaKcAa3NZbJ4sSvqWVnpF+JULE6RDiePjSPqzv0oLY0WdW/1J29xDSXlJCvs6opRFEg82ailQb4NGyR9wJgxXSNsMRkOHpTJVdZCbi5cdpncjxkz5Phjj4lrpr5e3DFOqOOdd4rvXfO/pJ9k1D/Vi1rsSqeluXzhsS6FykqxUJ0Scn36pN/aTsXi80ma3kmT3IpNZWWaobETgFrsinKGVFXJFPk5c6IjN5xZks5g6W23ucm+/H74whckFPDAga41mzQePp8bgnn11bIdCMB//Zfsa82vGqXjSUb9U72oxa50OmJzpseuvVaqUygjGJTj5eXRFq9TZLorL3l5MokqJ0d87EOGuPdBqyClDXTwVFHOAK+bxRvpEm8w0CkL57gtvOtMWRy30vTp0cWq587VGaRpJFlh7+K/GxUlRdx7b3QRiNJSmXhz550yGOh1y1RXQ3a2bDtui84+g7Q54rmPDhyAfv1kMNTJCQOS1VHj0Ts9KuyKEkvsLMlAQJJ63XGHCNlnPtO186nHkuih9NFHcPHFrqhrFaQug1ZQUhQvzizJGTOkYEZFhSw33SQif/bZkqmxK80gPVOcbJS9esGRI/Lr5bHH0t0rBa2gpCitw7FKZ80Sgb/gArHO9+wR98vHH2e2qPfoIaI+fbp83unT5YG2dGm6e6acAZrdUVEcYmuYzpgB3/qWWLBr1kgCL284YCZy4oRMtvr1r93xhM9+VmqcOi4ZpdOjwq4oDrEFIS64QKzzjz+W45lsqWdnw6lT8uDaskXGE4YNc4tmqKh3KdQVo3RvnLJtTpm2lSslS+EVV8jgqc+XWQOliTh1CoYOlc86apQ7aKwl7LokKuxK92bHDrjuOhGxmTOlTmk4LOvjx2HaNInsznSGDoVDh2QmbU2NzMDdvl3uSXFxununnCEq7Er3ZtYssVIffFDE7a67xOXixHZv3pze/nUEfj/s3Stx66dPy/04eNCNWdeUAV0OFXalexMKwapVkqlwy5bo/bm5UtIu0wkG5VfK229LnhxjpIRdaamKehdFhV1RIHpgNBCAl18WV0x3wKn6VFoKt94qr4NB+N73dFZpF0WFXekeOIOkXqqqJFPj1Ve7A6Q5OWK9O/VI/f6O7WdHETsgPGgQ3HCDjDcEg/Dww24sv4p7lyNlwm6M8Rtjthpj1qTqnIqSMpxQRifKY+lSGRjdtk3E21oYN65pSGMnGDi1CV63mrw8ifrxsm8f3H23CP6qVRLe+MQTaUsZ8NHfPuKxix/jo7991OHXzgRSabF/A3gjhedTlDMjkVW+ZEl0haS//lUGSSdPFr/6kCFipb76qmupO3SCyUgmwetWM2xY9KBwTo64n7ZskZmmXr96KOTGsncQp4+d5udX/5z9f9/Pz6f+nNPHMqQyVQeSEmE3xgwGpgI/SsX5FKVVeK1ycPO+OOF63qyNJSWSufDSS+HNNyWOu5OS8t8Mf/ub+/qqqyQCqK4Ozj9fImHS7HpZ/eXVHPvoGFg49uExVv+v1WntT1ckVRb7d4F7gfSbN0r3paW6pd6sjVu3isW+eTP07dvUUu9EtNv0qOxsuUfBoLilvvIVedil0a++9Sdbeev5t6g/WQ9A/cl63nruLbb+ZGta+tNVaXNKAWPMNOAja+1mY8zEZtrNBeYCDB06tK2XVZT4OFb5woUi4F5R9wp9797ijjnnHHj//fT2OV3k5MDTT0vqhOeec/c7fvU0hDq+NO8l6o5Fj3PUHa/jpXkvMfrLozu8P12VVFjsnwamG2N2ARXA54wxK2IbWWt/YK0da60d27dv3xRcVlHi4FjlkyZJnU7H8qyuhnnz4JFHZOB00SIYP15EPTc3vX1OFxddJOMKY8ZE70+DX91h0qJJBPOCUfuCuUEmPzw5Lf3pqrRZ2K2186y1g621w4BZQKW1dnabe6YoZ0pVlbgU5s2D+++XiJbrrxch37FDLPhhw8RNc955MlhYWNh94tW9+HywcaMMlg4blu7eNDL6y6MZOXUkgRxxJgRyAoy8ZiSX33p5mnvWtdA4diVzqK4WF8z8+eJDf/ZZGRS85x548kkJ8/vxj6UqUHW1TJ+vre3QLqYteNJJkXDhhdC/v0T7+P0wYEDarPNETP/JdPL65YGBvP55TP/x9HR3qcuRUmG31r5srZ2WynMqShMShTWCxF8vWCAx2T/9qQh7OCyzK//hH0TMX3tNBkvTEMqYtjyR4bC4nt56Cz78EC65RPbt25euHiUkKy+LL73wJfp+qi9fev5LZOVlpbtLXQ612JWuR7ywxmuukVhsEHGfPVvCGusluoKcHPjd72DkyPT0Od04rpeGBhlc/vOfpUj18893yupI/S7qxx1/vYN+F/VLd1e6JFpoQ+l6eMMaS0tlsHTBAhkQHR2JnHjqKVlbKzHrf/qTlHw7dCh9/e4oYqs8FRVJKl6QB5tTv9QpnqHVkTIOLWatdF3mz3fDGnv2FIv9wQfFSj91SqzT886TgdPzzxc3RBfD0gb3jTFNUyJkZcFvfqNZG7soWsxayWy8k42caj+LFomQHz8uA4N5ebB/P5x9toh6Xl66e92ElsyqeKLeoinm80nem2DQTfbVs6d8fidSSBN7ZTQq7ErXwzvZaMECWS9aJNkJt2yRSI+6Ojk2fTocOCACd9ZZkhemE3Gmwh3Pgrfe9+Tni1V+ww1yH6yVyVhHj8LNN8OLL8KNN6YlsZfScaiPXel6VFdHpwoIhSR2vaxMXC5vvy3r+fMlCsbnk+WDD0T4OimOaJuYfXj2mQT7Gjl8WEIa77pLtktKYN06mYS0bJn8onniifbovtKJUItd6Xrce29TH3F9vfjbP/pIcqC8/TYcOyaDiD16iKumrg4+/rjTFqeO16vm9kVZ6g6DBslnN0ZK/W3cKFWRtmwRwV+/PpVdVjopKuxKZlBcLO6YVatg7Vq3QEZDg+RC8Q4iWiv+505OrHBbzxLXWgcp5RcKiTtm926JXV+9Wh5u+/bJZC0l41FhVzKDWPdMVmRSi88n1uqpU+KGcIgtqNEJiRXvWDeNAcI+HxYIe3+FvPyy/FIJBMQNc/o0rFkjDz31rXcLVNiVzMBxz1RVSXm3rCypEhQOi1vC75ewxwED0t3TMyLW3x7rejHhMK8UTMdYZBJWVhbU1dGAj3ob+dXiiH4ak3spHYsKu5IZOGkGKipEyObPh9dfF4s9EBCXjDEygNqF6pjG9aMDBAIcvHAcJ8lhNFspy32U/ZdPBp+Pw+dfgak7hTFhGVAOBDTEsZuhwq50fZYsgVdeEUsdxOVQVSXWurMMHOjOxuxERTVaCm2M9/pNLmRPvyv49c5LCeRlkX/jFL64cDTmtVf56Q0vsGr3WMI5ufjzcsVKf/ZZGVeoqGinT6F0NlTYla5PcTFs2CCCXVEhyb/WrBHXRHa27E9Dsqtk5nTHi3Dxvs9JDGAcXzqGEb6d/HTfP/HP/mfYEHqAHYxgdH01z81eyZeXh8i+aAS/vNnjU3fEfcSIFH0ypbOjKQWUzMDxrR8/LqGP2dlw2WXiU1/tqZkZb5p9mkiULsDZ33g8GIS6OuoCOfga6vjA9mcQ+3gqaw5X+9Yy066k5KEQixbBlCmwYoXk99L0L5lHsikFdIKSkhlUV0th5nXrZNvxra+OKYTcwaIeK97N5X6JPdb4uq4OAgGs8XPawkDehwkTmLlhBXfxKFnTQ9x1l8xFWrtWRH3RInlrfb2Ol3ZH1BWjdD4S5VtfsiTxe3btElEPBmU5fVom5zjk5LRLV5vDa3nH7qeFY03cMvX1BOqOEfDDyqzZsGEDLwcms8T3LRb/7qrGCaZTpoilPm+ejB8XF6f6UyldARV2pfMRL9/6zJmJVaqqSvzqOTkyEWfiRHeA1OeTtLUnT3ZI172TiBxirfDY7URtwfWxg/yzBhtO8sX6Cp5kDhPr1xEI13H2qH5s3SoTTFeskJQwixZJfWpN4tg9UVeM0vmIl2/dO/kolkcekbqdX/0qHDwoqQUGDZJZmOGwm4s8xcRzq8RL0OXsb+I7b+aYDJS6lpdX/P3hOm7y/QIbhhO+PG7ZdicrX3Bvz/LlIvLdUdSttXxy4hN21u6k5mANNbU1XH/h9Yw8u3sVWFFhVzonoZCIupNvvTmVuuceidMuKxM3zJw5om7tSKybJZ6gm2baxHPTNAD+mPbOew6PHEP+W5sJY/Bj8Yfr+YNvAv8n/CBfv7yaUChEVZX42OfMgV/9Cm69VYYeioujb19Vlezvqr7343XH2XVwFzW1NdQcrGkU8Z21O6mpreHI6SNR7YcWDO12wq5RMUrnxHG/JGOxg5R3u/tuEfZAwC20kUaSGSiN/e9zrPQTJo8e9hgGqO/Rk8CJo/xtUAkXvFeJnwbC+DCEeWV6OX/69J0UF8P3rquifFY1I564t/H2zZsnbhnn9nkzHndWi74h3MB7R95rFOpG6zsi3h8c/SCqfY9AD4YXDqeosIjhvSPrwuEU9S6iqLCInlk90/RJUo9GxShdl1j1CYWit5csETPUa45u3y5l8bZskYFTp/5pnz6Sj70VeIW5peiWJlP9W9jvPYfTJoxY7IfoRYE9gsVgsfhPHOWTMSVctHkdFniMUt7hPJZyFxNX38XEz0JFxWhWmplkzVoJuN4sJ4XOmTwj2xtrLbUna5sIt7N+9+C71IXdXD4+42NI/hCGFw7n6vOujhLu4YXD6ZfXD9NJM3amCxV2pfMRL9+6o1KhkOR8+c534IEHYNo0uOUWGTU8fVraOzU/x4+Pjow5QxK5V5xjCcMTkzyP149uEVF/jwEM5IOI2FtOE6TBl8VZf36ZsDGcstk8l3UDGwIhchrgPxq+hf+xx5h16BCsilZs55kIyXu1UsXJ+pON7pJY4d5Zu5PDpw5Hte+T24ei3kVccc4V/POF/xwl3kMLhhL0d/5snJ0JdcUoXY+qKvGpWwv9+0v+cYfSUnjySYmCSaErJtHAZ0s1SRP54uOFPH5EH/pxAIuI/QEKyQmEORqazoB1y3nZP4lHs+/n3lA1Pxt0L089BfNOzOffT0cUe8GCuH04U69WMoRtmH1H9jVa3bHCve9I9EzfnECO6x6JWNpe8e6V3attHeomJOuKUWFXuiZVVWKtHz8evT8QgNxcKQe3e3erTp2o/Fxs1IpDa+qSJoqAaYgMjr7POZzD+6zxTScUXs/LV9zFZ/++jN/cupI7ng6xciVsWVrFrWtnsvGyUj63fRnbF6zk4OgQFRWSPaC4WDIsPPNM63zstSdqowYlveK96+AuTjec9nwew5CCIa5ox4h3/7z+6i5JAepjVzIfr1HiuF/q6zlx2pCze3ezJeaaI56ANyfmzcWsE9MuDNiIeHvxRdr7saxmOpuzP80tBauY9tFqfugv5Z53FvDdG0LMePx6Xp96I4OYxac3zGRG9kp+vyPED28OEbp7JvdlreSPOSHmzxcBnzEjsVdr/IRT4i6JhAXGukwOnjwY1cezepzF8MLhXNb/Mq4fdX2UeA8tGEqWPyuJu6t0BGqxK12P226TCUl+v5v7JWYCkgWs8eGz4ah48GSIFermRD5eZEtzhDG8xlWMYTNBZIDQ8bEb4DRBvsCLBEpCfHHdbXzZ/yTHz7+ULw56jYaXqngh+zp63DKr0SSvIsR110nqgE+frmJ0fTXvz7mXtWuh4qkwo8a+H2V17zzoDljuO7IP6+l9TiCHYb2HRVvcHvHOz+689WK7C+qKUTKXa66R7I0A06fDb38bd2bpQQrI51CT2Z7Qsm/c2yaR0ENTUY/nh/fGpzfgx2IJEOYxSvkyPyUH6fsOcx597EcYDNeximBJiC8XVTHjqZl892Qpc8PL+NceK/nGs6FGC/zQyUPc/dBOfvTLGijcyTmfquH9Ezs5+7wajgZ2carhlKdvhsH5g6PDAT3CPaDnAHxGJ6N3ZjpM2I0xQ4AngQHIL80fWGu/19x7VNiVVrNkiUTFfPBB0wRfMXh91825ZZr7D/BONIptm2iWqXPMu+3Ene9lEIN5jzqC1BMkyCkCNLCa6UxmPY/0XMi/HX2Qp5jF7TzB3NLTFKz6NyblPcbqS6/m46kXserlGs69fCcH6muoPVkb3YkThZxlhnN0bxEzJg7ns5e54j20YCjZgexmPq3S2elIYT8HOMdau8UY0wvYDFxnrf17oveosCutxjv6t2iRm83RIS8Pbr6ZEz98kpz6Ywkt7JZi0p0l1n5tyfL3btfjw08YA5wiwMuE+CfW8SIlhKgii3pOmgBf73kfP+z9eS4sXMNZhZvZ1ceH6fUWewvrIP89pO6dELRBevuHU0gRF/Yfzm8qiuBgEVeeN5xtvyui/mhvgkFJAOadmKRkBh02eGqtfR94P/L6iDHmDWAQkFDYFaVZnAlIiebBr1wp4Y4nTkS/7/zz4YILaPjRT3k9fAXjcGPYY6f4J5pg5I1TTxSz3tzkJKdtPX58kSPvZvfiQOERDvd+iX8pHENB4Uv8ZyG82xtqCus5FVgALOCNyDmyDvehX20v+td8iptrP+H5w1/ntvlT2b50Hw9V38EtWcv4xrNYshgsAAAgAElEQVQhqqvhrQ8k2vPCYngwUiCpokJ87t7Qf6V7kdKoGGPMMGA08FqcY3OBuQBDhw5N5WWVTMPJ7hgbozdjhpvx8cQJd0KSw9tvw549hMOWcWERda/oJhNL3ty+WJz9p/3wZkEWHxSepqY37CyEmsIGdkVef5Lr5C4JA5spOAlFtZB3YAhz336fC2rrWVf7FX5dezcTDu3iV/WzmR38MaMbqnkt/C1qeoS4pADqr4PrN/Tlf19ezcyZIUaPhjfegHHjYNas6MiXrp4PRmkbKRs8Ncb0BH4HfMda+0xzbdUVo7RIvFk1ILHrl14KmzaJWQoygLpxo6QOMIZD+YPJP7Sn8VSOkDfgwxdxjUDLIZAWOEkWOZzmg54i0jsLoaY37PK83psPYY/PJtgAww6KeJ9Tm83Ig3VsrJ3NB7Xjeaj2B1x58h0e5AGC1FNNMc8xDT9hHuEeSlnGTFbyii9EOCwJvVavlsCfBx6Q2/L88zB4MOzZI8U15s2TWzV8ONx4o2RWcJ6LW7fC+vXwwgup+bMo6aVD49iNMUHgV8DPWhJ1RUmKRNkdb7lFhD4/X6z2KVNE+XJyGlP15tYdAppa5r5IdvN4qQEOZ8Ou3lATEesdhY54n2ZXbzgRM6P9nCMwvBYmvAsjakXIh9QGGVwb4NwjdfSw9ZGB21MsoIznkVmh6znEIop5Gdc/cg1ruJ/vMJ+FPN6vjFcOiKj7fPDyy1K69Jpr4P77RaDffRf+8hcYMkSEG0TEr7kG7rpLhhmee06O3X23VFRSuhdtFnYj08l+DLxhrV3a9i4pCmKaLlsmor5smZv4ZNgwGDMGNm+GSy4R5QsEZHLS8uWwaBHBdetoIHrgs94H7xaIcO8shJ0Rq9vZ/jg3+vL5J6HoIIw6AFPeFuu76KCI+ZCDkFfvPiD204e+HMBQh43Eptfjw4flSWZTyjKqCPEyIR4hvm/kMraxgDJKP1rGr/whKiMW+yefyK0IBOCmm+Dpp0XUBwyAvXvh9ttd63zGDLkFx46Jdf/732vt0+5KKqJi/hHYAPwFt+DLfdbahD/+1BWjNEvsvHfvNsjrc85pVLjTnxyBoI9Pvn0nNT9+lHeCx9hVCO/0NuwutNQUwp4Yd0nAcZdEXCZFtSLajngXnogebPVO/QdDGD9B6qmLzCEN0hDl0jlJNvN4iPtYxEPM4z4WMZOVUZY6wESqWMnMxmPe7Z7TxI/u/Gjp3Vss8tJSuOEGsdCPHxdx371brPnZs+EXvxAv1YQJ8Mor7fZXUtKATlBSui4JomKO3n8PNR++yc65M6lZt5Kd5+RQ0/AxO872UZMPJ7PCUacZcEQEe9hBOK/WfT28FvofNmRb2yRaZhNjGMvmZicmVUfa1CH+GWcGqdfFc5IcrkZsm2KqqaaYYqqbWOz3sITqGNfMRKr4p8Jqam64l4oKuPZaKZxx8cXyY2XYMBkUraoScS8ogPffd10ud98tk3IbGtRizzRU2JWuxZIl1I0ZzZ4rRrh5S7a9zM59f6VmQA473/sbBzgW9ZZep1xL+9yDMKzWx/m1YYYc9DHiYJjcOpqkE/AK9E6GMZxdAByiJz05hh9LmGixB3eCkY0MwP6B8VzJn8iiPhLe6CNImBcp4dNsJMhpTpDL9axqYqUnSzAoGRMCAXnGbdgg26tWuc+8+fPFop8zBy6/3PWpjx4ttUeef17FPZPQJGBKp8Nay/7j++MXWAi/wZ5XPqDh9277QAOcmzeQouxPMWNAiKJnKinqP4qhG/7OkI/CDDx+Gt/UqRzcVUvB3zbiegJl3YAhEMmGEsbg82RGsUARuyKiHKAXR6kjiI86LAbjsddrGEYRuxpnje5iKOPZSH3EYhfRhyeZwxTWUsYCRrEdEGv9TITdsbQB6urc9Vlnyf6LLnJFPXYYYvv2aBEPhUTc169XYe9uqMWupJRjp49FJ52KyRh4vC46zW7/vP6NuUqKDhqG//evKZp4PcNXrGHQj1YS+Nxkt7Hja+/ZE3bt4qmsORy/8VZmrLienhzBb11XTD0+fD4fvnB9lEulHj8BovO0nyZIFnXU4cdgCCBhlFsYzcX8jSxOs52RjOIttnEJI3Hyv1ueZiazWc5x8riG5wCifOapYPRoiXDp0UMs8Hjpd7tCyTul7ajFrrQL9eF69hzakzBj4P7j+6Pa98zqSVHvIkacNYKS4SVRCaiG9R5GXlZe9AV2zXdHC72iDqJYU6bA8uXsGTaBz+9dy4zlt3LZoAmMfm9No9tFSsyFCYfhNAEC1OMD9tGfgXzIYXqSz9HGdYB6NjGGy9nGxxRSwBHWM5kNTOBuinmBqxnBTl5hAp9hA9s5n/58yPU8SzHV3EU5D/Ags6jgdp5gJivP2FJvjm3bxB2T5cmK21KRKaV7oxa7EoW1lgPHDzSp/O683n1oNw3WtXj9xs+5vc9NWBmnT26f5AssOKOBM2bA2rXR5mh1tajb3XfD7Nkcf2Yt9x+bxwM8SIAGcjhBgDDvM4ABSLFjA7zJSEbyFvX4CBDmvYg7pR4/fsKsYzKfZ12jK+UEPbiRpxp94xOpYhXXs4EJjONV1jKFOSznMUr5Go81dn0iVXEHR1NFSQl84QviU3/uORXv7ooOnioJOV53PLr6e4zVfawuepCyX16/hGXNBucPJuBLwQ+/iC+hctw8/nHDIrLmz5MsVvPmcXrBIn4/YR6fq5oPCxZQNfpOdk+5jVnmKSrrJnBhwzaGsoe3IyL+QfZQBpzazTFyyeM4f2A84/kjh+lFAYcbqxSdJIcA9Y2ul5NkcYoeXM+qqEiWOgKNYYsS4hhgIWVcw5qUWeXNEYxMjurRQ2LZd+3SmaTdFXXFdGPqw/XsPbw3bi3KmtoaPjz2YVT73GBuo3B/ruhzUeI9rPcwemb1bP9OR3wLhhAzXxnNygUzyZo+hfr7y5iZtYbygdWNpmr1Eih5aBbZD1ZQNATOefsD/s4oBrOPLYOnc8Xe1WxnJOfzNnX4uYrX2MwVFLMZC+xlMP35iDCGIPUcJ5tfMpObWU4Yv3SH4kZfeTHVTWLRX2d0St0tzTFxosSjNzTAz34meWGqqhLnSFMUtdi7INZaPj7xccJCwrsP7aY+XN/Y3m/8DCkYkrAyTt/cvp2qHmVVFWyeOp+7Tyzk0R5lvD1nQWOSKyfEvaICzv5zFfe/Oo3f82lKWMdx8vgZN3EFmylmM6cJMo+H+Q/uJ4eTjZb6aYL4CROgIZIPJptT5PBffJ2v8188xY3czhONk4WWUdqYw6UjhDweEyaIr33WLFl04LR7ohZ7F+d43XGpRxknsmRn7U6Onj4a1b5vbl+KCou4ctCV3HjRjVHiPTh/MEF/MMGVOh8hqrjSLGMBZdxllvHWBSE+P1MKOO/YAffdJ+0aGkLk+e7i/vBCnmQO0/k1N/M/5HCK0wQJUM9w3mkMR/Rj2c5ILuCtxhj1k2TjI0yAOqoIUUWIlczkTS4gSD3LKGU+C1lAWYeKeiAAX/2qhDGCxLDPmQNPPCHbK1c2zZGmoq44qLCniYZwA+8deS9hWOAHRz+Iat8j0KPRwp547kQJD/RY3R3iLkkx8SaYvjftNgp/W8FNPZ7l0rIQN30vxNPzr+P1ibO4fOYTTJnixnlPza1i7vFlLKSM21nGa1zFPyGFNx7mW9TSm6XcxWkCHCMXP/WM4q3G8Mc6AqxnMiGqANMY1fIQ81jIfMpYQCnygPHme2lvBg+WPDCOqPv9cNllsGKFRHo6M0/j5UhTFFBhbzestXxy4pOEMd27D+2mLlzX2N5nfAzJF3fJ1edd7cZ2R8S7X16/TuUuSQXx0q7v+i3MqDM88DCMrl/C2zcFOP64ASNCtmFhFb+4opofvV7MT4+78eKf0Jty7uI0QeoI8HX+iweZz0my6MFpnuJfmMEzWCSH+zYu4RL+SoiXKWMho9jODJ7hQ/pTyjLKWBDlU3cs+fZwxzj1uEEyOu7dKwJ+NPKjbMkSmWB0xx0i9uXliXOkKQqoj71NnKg7Ie6SODHdNQdrOHzqcFT7Prl9JLKksIjhvaPDAocWDO1S7pJmaakCEtG7Z84Ui3TTJsknfvsFVYxeNJP3R09hwLoVvFP6KI823MmHFVX8rH4mN7KSi09WUxTeQQWzAFjFdQRooIoQ+xjIjVTQg5PUEeQZZjCbFfydC7iIN/kj4zifd1jJDZTyOM8zlek8x4PMZz4LWc8kdjCCCmY1ivhEqphFBTsY0W4hjYMGwccfS/RLbS0UFkL//rBvn9zObdsk9/r27fDMM+pj744k62PHWtvhy5gxY2wmMPPpmZZv07jk/EeOvfD7F9qpP5tqv/b81+zSjUvtqjdW2W0fbLOHTx5Od3c7jspKa/v0kXW87RjKyqwFa3v0cJvUzJGd+0rm2FMFfewjPcrsqYI+9tFplTYQkPYTqbQf0cf+kSvtaqbab1JuP6KPnUil/T6l9hR++zhz7T0stt+k3J4kaL9PaeN7nf1rmNJ4rgcps7UU2Fry7UQqo67jbKdqycpyX/fqJeuRI2VdVCTrCy+0NjdXXpeVyb1ZvLjprayslP1KZgNssklorAp7G3hp50t2xbYV9g+7/2DfP/K+DYfD6e5S58ER87KyZkXd2ywnR0TsJ3Mq7X7Tx9bMETFfN3COtWBr5pTZYNAVQ79fRLeWAnuUHraWAjuRSjuRSnuEXPtNyqOE1BHzWIGNFW7nnLXk2wcpaxdRj7cYI+vBg+W13y/bwaC1kyY1exuVboIKu5J+HFPcMTVjiDXky8tdkf3JnEpbWWnt/bnlNmyM3Xj+HLvf9LEz+zYV2AeR6xylR6uE+B4WN2k/kUq7jknWgn2QspSJt/fBBNbm57sPKb/fFXdn7fdbGwjIPWrhh4/SDVBhV9JLEhZ7rEuhstLa+4OL7WR/pe3Rw9ppeZX2VEEf+9K0cvu7qYvt54OVdj997JzBlVEC7LhQjpKbMiH2nre9LfY5c+SzT5sWLfz9+4vAl5dH3yN1uXRfVNiV9NGcjz2Bg/iduYsbm0wSQ9neF1hsH51WacvLrS0osHbcuGh3iteFIu6TfHuU3Cj/eFtEvT197I5P3bHiS0tdV1Rurm0cRygpUStdcUlW2DUqRkk9zUXFeGIcl1SHKAlIBMy3hq+k343S/v77pbSbzyfVgQ4elJKm4bBIoYNTfQjcVLkAs6hgBs+0OjQxUVWjVCf56t8fPvwQ+vWDjz6Setw33wwXXAALFkitbp8PvvMduR+aLkDRJGBK56WqiuPXzOTP40sZsW4ZL5eupM8NIa6+Gk6elHj1LVvgtddkBmZ9ffOn6yghbg/69IEjR6CoSLbPOw9efdUt77p0KVRWwpo1GsqoaLij0olZvNjaP5bIgOfzY8qsMdaOGWMbBwuDQXFHlJTIPsctkalLbq7413NzxR3j+NTLy63Ny5O1+tUVa5N3xejMU6XDcDw0/g1VnLd+Ga+WlHHl+mVMywvx3OYQOTlisTc0iAvilVfkfS1Z7F2d+npJP+/zwcMPS7bi11+XFAJTp2pZO+XMUWFXOoziYvjedVX8vGEmM3us5OWNIcaZEL84Gpmqf1J8DcEgnDol7/FOt89E/H44LVkOCAbFt37ppbB8uVRMUlFXWoOv5SaKkhpCISifVc2X/Cs564shjh2Dl8IhbozkPHdS4dS5KXQyWtTBTWgWCMDkyTJgumGDbPfokd6+KV0XFXalQxnxxL1c+o0Qy5e7+6oIUe67F2vFgu2OfPWrMkjqWO9ZWVIGb+ZMCShSlDNBhV3pUKqq4HvfE38ywPnnQ16ea506Fmwmk53tvvb5IDcXdu+WXyrBoGRsDATELTNvnkSJKsqZoMKudBhOrepPfUos82AQ3nkHJk2C48fT3buOwxk/ABHwhQvhpZfkfrz4ogj6s8+KG2r7do1fV84cFXalw6iuFtH6y1/gK18RKz0QgNWrM9+XDpJj3SEQkIpIdXVioU+eDLNnu7HqoZCI+4gR6emr0rVJibAbY75gjNlujHnHGPOtVJxTyTzuvVeiPJ57Dp5+Gq69NnqgNNM56qlmePXV8OST8OijMqN24EC37J1DKKTWutI62hzuaIzxA/8XKAH2AtXGmNXW2r+39dxK16O5bAIgNUtnzYLevSWkL9PDGeMxfrw83JYulQfd9u1N2ySoS6IoSZEKi/1K4B1r7U5r7WmgArg2BedVuhhLloiLwRvJcc01Yp0WF4uor1gh2zt3yvHuJuo+H1x8sVjq69fLvlmzpCKSc8+cikjFxenrp9K1SYWwDwL2eLb3RvZFYYyZa4zZZIzZtH///hRcVulsFBfLrMl580SYbr4Znn9eIl22bhUBC4dldml3EPTYErXZ2TKuUFEBo0fDCy/I/lBIcsPMnOmGOGqZO6UtpELY41VYbvJva639gbV2rLV2bN++fVNwWaWz4QjUokVw4YXiapk9W6bJ3303/PSnIug+nyvsvXunt88dyejREg3U0CDi7iUUkuRnCxfKWkVdaQupEPa9wBDP9mBgXwrOq3RBQiGYMkVmT06YIDlQRo8WgV++HEaOFKvd4eDB9PW1vbEWhg6FUaNg2jTJVnnTTbLERrtUVcGyZRIhs2yZTkpS2kgymcKaW5AB2J1AEZAFbAMuau49mt0xcykvl6o/c+ZIgQinSEZ+vrUDBrgZDZ2ScJm85OXJeto0995MmdL0np1h7W+lG0OS2R3bbLFba+uBrwEvAm8AK621f2vreZWuR1WVWJyPPiqhfCtXwoMPSvbCCRPg44/dtkeOpK+fHcWxYzLpyEnkdeedrl/dS3V1tE/dcWnpjFOltWihDSVlOKGOTqGkUAhuu03S7+7YIf50Y6RaUHdh2jQJbVSUVJBsoQ2deaqkHKf6XVWVRMJ8/LFMRKqthe4WENUdon+UzocKu9ImliwRq/y229wY9q1bYdw4+NznZDl8WBJ9BQKZI3SFhc0fN0YiYJ5/Hu64I/pYVZXcN0VpL7TQhtImioul2PLJk+JPvvlmuOsuee1w6pRExvzpT+nrZ6qprU18zKljGghIGt4f/hBuuEFcU87kI6emqaK0ByrsSpsIhcTd8uSTMli4bJlYqk4OmD594MABKdCcCThx6M1RVwcPPRSdKmDmTIlPX7ZMJx8p7Y8Ku9JmZs2SVAEOjvD5/eKSeeGF6Nj1roox8UXdeXgB5ORIBNCiRdEC3r+/TD4qK1NRV9of9bErbaaiIr5wjxwJa9ZkhqiDjA/EpgkAEfWRIyUC5tJL5deJt0CGTj5SOhoVdqVNVFVJqoCTJ919jvi98UZ6+tSeeAd/i4rc17t3S5z6a6+JpV5fL5kZvT71BQvcnDAq7kp7osKutInqavj8591Sd5A5kS9e/P7ozzhmDOzaBdOni/tl8mRXsL151HXykZIOdIKS0maqqkTcrc3cmqU+n7iUjIHbbxdRnzzZzWZZX+9OztIc6kp7kewEJR08VdpMdTV89aviP84kvEVAwmHIz5fwRW/o4owZrtsFdGBU6RyoK0ZpMzt2wI9+JC4JXwZ9o6wVF4zDkSNw663iclm6VNazZqmFrnQ+1GJX2swrr7hx65kSAePQ0CDFMU6cgMGDpVbrlClSEOO559RCVzonGWRfKR3FkiXRUR2f+YxUB6qvl22vldvVMUZEPTtbUiNMmSJ55e+8U0Vd6byosCtAU7GGxDlNnCRft93mJvry5oFpaJBJO5mAU/EpOxu+9CWpTarx6EpnR4VdAaIzMkLzBZUfeUQGEJ98Er7wBbj+esjNjW7jzMTMBMJhmDgRHn9cYtE1Hl3p7KiwK0ByBZUdq37yZBG5s8+WKJFDhzI7HW8wKBkrH33UdTdpPLrSmdHBU6URb0HleDlNHKt+5Uq48kqZZZmpZGeLe+nii8Uds3OnZKj03pNQSP3sSudELXalkUQ5TRxL3bFSr7sOtmxJb1/bAydU0+eDkhKJetmxAx5+WK1zpWuhwq4Azec08Q6WApx3nhvemElccIG4XXJyYMMG2ecIujdNgKJ0dtQVowBuThPHKvX6kIuLZYZlRQX8z/9I4Qxwp9lnCm+8AeXl4nKpqHDTA6i7RelqqLB3Y5zi015rdOtWuOaa6Mk3jiV/4oTEcPt8MGmSWPPhMAwaBO+95563qwp+IOD60VXMla6MumK6MfFCHBctEldMbHRMRYWs58yRgcV162RQcdw4EXaHYLBribqTYrhnTwnZrKhIb38UJRWosHdjEoU43nmnGx0zZYq0ragQ33PPnuJf9/nciUhHjrjn7Aq+d7/fnR1rrQyU5uTAAw/IBCSNTVe6Oirs3RxviGNpqZu1cNkysc6XL4f77pPZpQMHyv4BA8QqHzoUfvObrlNQw+eTz+j3u7NkjZG0uytXuil4NfpF6eqosHdzYkMcnayFK1fKzNLSUin19vzzroDv3SsJsXbvdifsxCsZl24GD3ZfZ2eLVb5nj7iLjJEHV69eMnMWoisfKUpXRgdPuzHeEEdnwHDaNLHeHav1scfgnXfEp+7NT753b/S5OmPVJKePTqTL9dfDiy9CVpZ8HufXyXXXiavpiSd00FTJDNpksRtjHjHGvGmM+bMxZpUxpneqOqa0P/HKtq1Z41YDuuYauOMOiZQpKuqc4h2PwkL3tTGwfbt8tlWrROBvuin6Mz/7LIwYkZ6+Kkp70KbSeMaYzwOV1tp6Y8xiAGvtv7f0Pi2N1/lZsgT+8AdYvVrqe27eLAOlXSm5l9PvvDzNna5kBsmWxmuTxW6t/a21NuJl5VVgcHPtla5DcTGsXw8jR4o4DhjQeUQ9Kyu5dm++KW6YQEDDGJXuRSoHT78MrE3h+ZQ0EgqJr/2tt8Ti/eADcWuMGyfx3n37SjtvDHtH4PdDQUFybY8fF9fLqlXqalG6Fy0KuzFmvTHmr3GWaz1t7gfqgZ81c565xphNxphN+zM5x2uG4ExWKimBY8fcgdMdOyRtrxP/7Z1xCk3zsqeahgZJEdyvn2z3jhnVGTpUwhp9Punj0qXykCoujl80RFEykRaF3Vo72Vp7cZzl1wDGmFuAacBNthmHvbX2B9basdbasX0dc0/ptFRXSzGN9evhiiuk7mdODnz8sfjdR46El15y2zvhjsePt3/f+vWTfowfL7ngx48XIff7JRKmRw8JYywthcpKN4QzXtEQRclE2hoV8wXg34Hp1toO+JdWOoriYkn49eijsjipAsJhsco3bhTfNbgW8siR7defYFDWxojF/qlPSR8mTxZ30fr1UuUoHJaZs6tWSeHpL35RZtXGFg1RlEymrT727wO9gHXGmNeNMY+noE9KJ6C6WkIfnaLN8+dLuoCzzhKrvFcvqZ7Uq5eI6YgRIrCpFndntmhWFkyfLpE5vXvDX/4CQ4ZIPPpll0lI5rZt7kSrrVtlvxaeVrojbY2KOc9aO8Rae3lkuT1VHVM6jniFrIuL3UlKjr999mz45BOZ0XnkiFjuR46I2L71lszkPBNhdwpb9OzZ9FhOjqyNEb/6c89JZM6xY3DwIFxyicwiLSmRmbF33SXpABYskPXdd8t+LTytdEustR2+jBkzxiqdh8pKa/v0kXW87cWLrS0vl30lJdaCtWedFb0ePNja/HxrjZFtZwkE3NfeY36/rIcNi27vbTd+vLUFBdYGg9aed56sjbG2tFSuNW6cbI8bZ21enuwrK5P35OVZO3du/M+jKF0VYJNNQmM1V4zSYiHrHTvEEp4xwx2s/OQTmeH5yScSfrh3r7hlvMPn48dH55DxHmtoEL/5rl3RffGmLdi4ET77WWn7zjuQny/+/oYGcQPdcINsFxaKRX/eeRKiOXasbD/xRPTn0+ReSndBhV0B4md59GKtlI77619FcHNyoH9/WR8+LLM833tPRNYR840bm6bxPecc97X3mM8noYrWRtceXb1afPglJfIQef11eOopGTRdtEji1F94QXzqW7eKO2jbtvifT5N7Kd2GZMz6VC/qiul8OO6KsrL4bgvn+MCB4iaZM8dtV1oqbpIJE+RYdra1l1ziulZ8PmunTnX35eZGu14ct4yzlJe7Lh/HnVNZKdcEa7OyZNvp05w54pIpL4/uq7pelEyDJF0xKuxKEyGcO1f81V5hrKy0dtKkaL94WZn7XkfcvYLsXcaPl3WvXk198F5hLykRH3lWlutvz8uTh0FuruzPyXH7VlbmPmhiP9PixR1z/xSlo1BhV5Jm8eKmIl5QED34mJ8v+xzrODtb9s2dK6LuDGrm5cUXdsfS9g6O9u/ftI3f7w64ZmeLFe618L0PE2dAN9GvDEXJNFTYlTbhdc04ou4IaXm5G3nSo4dY0I4b5Moro6NfnEgWZzs7W9YjR0aL/NSpcg7v+5xfDVOnyvakSa6AO4KfKJJHUTKRZIVdB0+VKJyYdu9g6ogRcOONkqfdqYm6apUkBDtxAj79aTlWVQUPP+zGoPv9MkBqI1EuAwbAqVOS2/2tt2Sg1VoZeF2zJrofn/mM5Em//nrYsEEKZKxf70bvbN8u7/HmVdfIF0URVNiVKIqLRTiXLnXrnr7+ukTE3HtvdLSMM9Nz2zZJL+C8LytL3tfQ4LZ1MkSOHAk1NSLutbWyvWUL3Hwz3H+/PBQmTXIjW268Ueqtxgr4iBFNI3c08kVRIiRj1qd6UVdM+on1q1vrDjiWl4uLxIl8cVwwidwe3gFU5335+eI+8bpbRo6U144bxtkeM6ap+0VdK4rSFNQVozSHY5k7U+2d+qfFxeJWmT1b8qyUlorrxevmiFdSb+VKmWzkvO/UKYlDnz5dXDI+n7hfhgyRCU8lJVIMe9QoKeTRs6e4bR54wK2/qq4VRWklyah/qhe12DsHiWLXW4ppb+58eXnWXnGFDKiWlroW/7RpbrRLXp7scyJoHMveGxuvKEpT0KgYJRmcOPCyMtluKW9MImJDEB3hdsS9slJCIx3hN8Z9AOTkuA+R8gvF4nMAABExSURBVHKNP1eURCQr7IF0/2JQ0kdVlQyQOhkQQ6HEbpbq6uZT33rfN3q0uHVmzJCc7k70SigkaQgWLoQJEyTapUcPeP5597iTp0ZRlDaQjPqnelGLPf201jJPlthfAt5rOJOcRo+WePjYyVFqsStKfNDBU6U5mrPM20rsL4GqKndwdt48WLtWsjLu2eNmk3QGcTVkUVHajpGHQMcyduxYu2nTpg6/rtL+OALuPDSc7RkzJB69uloib5xjznZ1tQq6orSEMWaztXZsi+1U2JXWsGSJK9AOVVXwyCNwzz1N96twK0rbSVbY1RWjtIpEcfCxog7qXlGUjkaFXWmReDVRQdwriaouKYqSPlTYlRZJZJ3PmtV81SVFUdKDCrvSIolqokLT6BdFUdKPCruSFLE1UcEV+AULXOFXcVeU9KPCriRFbGx6RUX7xcEritI2NKWA0iKxsenO1P9Zs6LbOccURUkvarErLdKes1QVRUk9OkFJURSli9ChE5SMMXcbY6wxpk8qzqcoiqK0njYLuzFmCFAC7G57dxRFUZS2kgqL/T+Be4GO9+koiqIoTWiTsBtjpgPvWWu3pag/iqIoShtpMdzRGLMeGBDn0P3AfcDnk7mQMWYuMBdg6NChZ9BFRVEU5UxodVSMMeYS4CXgeGTXYGAfcKW19oMW3rsfeLdVF+5Y+gAH0t2JM0T73P50tf6C9rkj6Ij+nmut7dtSo5SFOxpjdgFjrbVd6Q/RLMaYTcmEFnUmtM/tT1frL2ifO4LO1F+doKQoipJhpCylgLV2WKrOpSiKorQetdib5wfp7kAr0D63P12tv6B97gg6TX/TklJAURRFaT/UYlcURckwVNg9GGOeMsa8Hll2GWNeT9BulzHmL5F2ac1mZoz5tjHmPU+/r07Q7gvGmO3GmHeMMd/q6H7G9OURY8ybxpg/G2NWGWN6J2iX1vvc0j0zxmRHvjPvGGNeM8YM6+g+xvRniDGmyhjzhjHmb8aYb8RpM9EYc8jzfZmfjr7G9KnZv7MR/ityn/9sjLkiHf2M9OUCz7173Rhz2BjzzZg26b/H1lpd4ixAOTA/wbFdQJ909zHSl28Dd7fQxg/sAIYDWcA24FNp7PPngUDk9WJgcWe7z8ncM+AO4PHI61nAU2n+LpwDXBF53Qt4K06fJwJr0tnPM/07A1cDawEDjANeS3efPd+RD5DY8k51j9Vij4MxxgAzgV+kuy8p4krgHWvtTmvtaaACuDZdnbHW/tZaWx/ZfBWZ3NbZSOaeXQv8T+T1L4FJke9OWrDWvm+t3RJ5fQR4AxiUrv6kkGuBJ63wKtDbGHNOujsFTAJ2WGs73WRLFfb4TAA+tNa+neC4BX5rjNkcSZWQbr4W+Yn6E2NMYZzjg4A9nu29dJ5/+C8j1lg80nmfk7lnjW0iD6pDwNkd0rsWiLiFRgOvxTn8D8aYbcaYtcaYizq0Y/Fp6e/cWb+/s0hs/KX1Hne70njN5b6x1v468vpfaN5a/7S1dp8xph+wzhjzprX2lVT31aGFfD3LgIXIP8dCxIX05dhTxHlvu4ZDJXOfjTH3A/XAzxKcpkPvcwzJ3LMOv6/JYIzpCfwK+Ka19nDM4S2I6+BoZDzmWeD8ju5jDC39nTvdfTbGZAHTgXlxDqf9Hnc7YbfWTm7uuDEmAMwAxjRzjn2R9UfGmFXIz/Z2E5yW+uxgjPkhsCbOob3AEM+2k9en3UjiPt8CTAMm2YhjMs45OvQ+x5DMPXPa7I18bwqATzqme/ExxgQRUf+ZtfaZ2ONeobfWvmCMecwY08emMRVIEn/nDv/+JsEUYIu19sPYA53hHqsrpimTgTettXvjHTTG5BljejmvkYHAv3Zg/2L74/U1Xp+gL9XA+caYooilMQtY3RH9i4cx5gvAvwPTrbXHE7RJ931O5p6tBm6JvP5noDLRQ6ojiPj3fwy8Ya1dmqDNAGccwBhzJaIBH3dcL5v0J5m/82rg5kh0zDjgkLX2/Q7uaiwJf9V3hnvc7Sz2JGjiNzPGDAR+ZK29GugPrIr83QLAz621v+nwXrosMcZcjvw03QXcBtF9ttbWG2O+BryIjOT/xFr7t3R1GPg+kI387AZ41Vp7e2e6z4numTFmAbDJWrsaEdHlxph3EEt9Vkf1LwGfBuYAfzFuqO59wFAAa+3jyAOo1BhTD5wAZqXzYUSCv7Mx5nZo7PMLSGTMO0g22VvT1FcAjDG5SNW42zz7vP1N+z3WmaeKoigZhrpiFEVRMgwVdkVRlAxDhV1RFCXDUGFXFEXJMFTYFUVRMgwVdkVRlAxDhV1RFCXD6BbCbowZaow5aozxp7svqcZIPvYV6e5HJmGMedkY85UEx/R+RzCSR31y5PV9xpgfdcA1Jxpj4s4Kjxy3xpjzIq8fN8aUtXefOiMZNfPUGLMLmcnW4Nk90lq7G+iZlk7FEJkSvRDJRdMAvAx8vSOmSBtj/hvYa639P+19rXQRyWpYAwQ9qYGVdsZa+1Ay7TryO2itvb29r9FZyUSL/RprbU/P0q7JgiLJn86EQqTo7TDgXOAI8NMUdytjycRfXW2lFd/BdjmH0nnIRGFvgjFmWOQnWiCyXWSMecUYc8QYs94Y83+dn9fxfurF/OT8tjHml8aYFcaYw8C/GmN8xphvGWN2GGM+NsasNMacFa8v1tq11tqnrbWHIwmwvo/k+Ej2sxQZY34X6fs6oE/M8aeNMR8YKc31ionkgjaS5/om4N6IW+q5yH6n30eMMX83xlzfzLWvNMZsMlIO7ENjzFLPsXHGmI3GmING8lBP9Bx72RizyBjzp0i/fu29P4n6HDn238aYZcaYF4wxx4CQMWaqMWZrpB97jDHf9nTTyQp4MPI5/yFyni8bKRlXa4x50RhzrucaJUZK9R0yxnyf+GliveQYKYl3xBizxRhzWeQ89xhjfhVzz/5/Y8x3E9zPCyP35qCRUnbTY+7ZVzzb/2qM+b1n2xpj/j9jzNtAk7oBnu/8XGPMPmPM+8aYuzzHz/h7bIyZY4x5N3Ls/pjrRbmojDH/6Pk+7In0P9F3cKAx5lfGmP3GmBpjzNc95+kR+Q7UGmP+DhQ383eJvQf/bYz5j8jricaYvcaYu4wxH0Xux62ettnGmEeNMbsj3+3HjTE9kr1WpyOd5ZtSvSBJsCbH2T8MSZLllGP7I/AoUvLsH4HDwArrlrXam+i8SCm6OuA65MHYA/gmbiWgbOAJ4BdJ9vmbSBIsZ/tLwJ+baf9HYGnkOp9BLP4VnuNfRsqiZQPfBV73HPtv4D9izncDMDDyWW4EjgHnNHPtOZHXPYFxkdeDkOx1V0fOUxLZ7hs5/jLwHnAxkIeklT2TPh9CHn4+ICfyN7oksn0p8CFwXby/dWTfdUgCqQsR9+P/ATZGjvWJ/P3/GQgC/4bkiP9Kgnvg/P2d9ncTcf0gpemOAb0jbQPAR8CYOOcJRvp0H/I9/Fzkb3mB5559xdP+X4Hfe7YtsA44C+jRzHf+F5F7fgmwn1Z+j4FPAUeR71w28h2sjzmf8z80NPJZ/iXyOc8GLo/3HYxcezMwP3IfhgM7gX+KHH8Y2BD5nEOQzI974/1tPPflvNhrId+ZemBBpE9XIwnFCiPHv4tkkTwL+S4+ByxKt6a1dkl7B1L6YUSAjwIHI8uzMV/yQORLVw/ket63gjMT9ldijr+B5BV3ts+J/NMEWujvpUhWwAlJfj6n73mefT/HI5Ix7XtHPndBZDvqnyrBe14Hrk1w7BXgQWLqUyIpeJfH7HsRuCXy+mXgYc+xTwGnAX+SfX6yhT5/F/jP2L+15/ha4H95tn2Rf+pzgZuJfrAaJP93c8L+asy53nf+hpFrfTXyehrw9wTnmYDUy/R59v0C+LbnnrUk7J9r5p4492GUZ98S4Met+R4jwlvhOZYX+RvGE/Z5wKoE/Yr6DgJXAbtj2swDfhp5vRP4gufYXFov7CdivhcfITVUDfJAHuE59g9ATTL/l51xyURXzHXW2t6R5bo4xwcCn9joPOB74rRrjtj25yKpRw8aYw4i/yANyEBuXIyM3K8FvmGt3ZDkdQcCtdbaY559jfUWjTF+Y8zDkZ/Sh5EHEsS4a2L6cbORSupO3y9upv3/AkYCbxpjqo0x0yL7zwVucM4ROc8/IsLg4L1n7yJWU58k+xx1v40xVxljqiI/3Q8Btzf3GSP9+56nb58g/8yDkHvaeH4r/9UtfR+87cPIg2BgZNf/ALMjr2cDyxOcYyCwJ/J+h3c5s5JvyXxvY+/7wATHoPnvcex9OkbiHONDkELgyXAuMDDmu3Mf7v9O1HXxfN9bwcc2ekD9OPLLsy+QC2z29OE3kf1dkkwU9pZ4HzjLSE5lB291lmPIHxloHKyL/QPH5jreA0zxPFB6W2tzrLXvxetAxL+7HlhorU30j5+o74VGChI4DPW8/hJS+HcyUs1nmHPJeP2O9OOHwNeAs621vZGfunF9zNbat621/wL0AxYDv4z0ZQ9isXs/f5619mHP2733eChiCR5Ios9N+o38SlkNDLHWFgCPJ/qMEfYAt8X07/+1dzYhWpVRHP/9+5CaGIYIQbAkcBE49gERtWgZOJWzGbCFShC1ySiiKCmqScldIC4KoVqkBBpDTYs2Q0UtLJMgFyXSxsX0Nah9GozFdFz8n7fu3N73jprjq+97fsv3ed57z33uec49zznPvefKiPgUj+k/sklSTdZ2VPtfwvyKPpPATZJWY4+9U9m/74Hryv9brMAhK6jpIe3LDJ7ON7fr417dTHAmelwfpwE613edBlZ2aGt3ziO1cw6Gv8lP/bzM1/dzxTHszQ9XZBiKiAtiJ93Z0HeGPVxR/AvgRUlL5OTaaKXLNzg5dq9cZuw5HFNsYiewrZWQk7RUUr2iPaVtOfAR8Er4o/xnI/uWIvudNdkHgZPYkxoA6lvQZnAMs8VVeKIdLbI9gD32tkjaKGlp8TJ/KT/P4VDWqKQ1xQO/oiSrrq38faOkVcUgbAUmImLuNGRuxyBedc3KFWrWV9qOAn/XrnMn8Iz+TSQPSVpX2t4HhiWNycn1x2hvRKvcWun/eJF/P0BEzAIT+OFzILzVth2fY+P9tKTL5WTzKLCntB8ExiQNlNXdgwvI1InnyzGGcYGKvQ19m/R4AlhbkqJL8D3sZD/eAu6SdJ+kyyRdIxeDgf/q4AHgN0mbS6L0UkmrJbWSpG/je3d10adHz/D6F6To82vAdrnuKpKWS1pzrs91vug7w17YgGNox4GXsLKfBIiIX4FNwOvYe/oDL7Wb2IE9yClJv+NJfnuHvg9hxR4vOwNOSDrRapS0QVJTdaP15dg/AePArkrbLrxU/Q44VOSo8gawqiw3JyPiEC5+/RmecDcC+xrOPQJ8XeTdgSvDzEbENPa6n8WGdRp4ivn6tRvHPH/ECdDWzoeFZG7HJmBrGesX8OQHoITYtgH7ynXeERHv4hXGnhLu+QrXrCRch3IdTtIdx0WHm8YA4D2caP4ZVywai4i/Ku1v4rHsuBqLiD9xMeS7scf4KnB/RBwuXbbjGPZMOV4nz38hPsFJ2g+BlyNiqqFvRz0OV9x6BD+wfsDX3nZelIfZPcCTWE8PAjeX5roOzuEH2i04CX0Mz72h0n8L1o8jwBQNY/o/2YzHaX/RkQ+AG2DeC46LsVpYFLKCEiBpL65zOt5tWXoRSR/jxNqiv5l4IVAMwGFgWVQKG59nGa4nX9TqW/rSY5d0m6SV8r7dEextTnZbruTip8TMn8A7SLpi1JOkX982Wwa8g5M/3wIPR8SX3RUpudgpieQZHDoY6bI4SR+ToZgkSZIeoy9DMUmSJL1MGvYkSZIeIw17kiRJj5GGPUmSpMdIw54kSdJjpGFPkiTpMU4ByGUTEJKJ9coAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111522dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x0[:,0],x0[:,1],'x',color='blue')\n",
    "plt.plot(x1[:,0],x1[:,1],'x',color='red')\n",
    "plt.plot(5,4,'d',color='purple')\n",
    "plt.figtext(0.5, 0.01, 'Figure 2: data separated by our predicted line.', \n",
    "            wrap=True, horizontalalignment='center', fontsize=12)\n",
    "plt.axis('equal')\n",
    "x_range = np.array(range(-7,7))\n",
    "plt.plot(x_range,-(theta[0])/theta[2] - (theta[1]/theta[2])*x_range , 'green' )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This algorithm not only provides a binary classification, but also, for a new point of data, it provides a way to assign a probability how likely is that a point belongs to the class $1$ (red cloud)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9774117781336015\n"
     ]
    }
   ],
   "source": [
    "def predict_prob(z,theta):\n",
    "    return sigmoid(theta[0]+z[0]*theta[1]+z[1]*theta[2])\n",
    "\n",
    "P = predict_prob([5,4],theta)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the purple has probability $P$ to belong to the red class according to this model. Finally, we remark that the separating line is not unique in general."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}