Logistic_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy.optimize\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import fetch_mldata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normalize_features(train, test):\n",
    "    \"\"\"Normalizes train set features to a standard normal distribution\n",
    "    (zero mean and unit variance). The same procedure is then applied\n",
    "    to the test set features.\n",
    "    \"\"\"\n",
    "    train_mean = train.mean(axis=0)\n",
    "    # +0.1 to avoid division by zero in this specific case\n",
    "    train_std = train.std(axis=0) + 0.1\n",
    "    \n",
    "    train = (train - train_mean) / train_std\n",
    "    test = (test - train_mean) / train_std\n",
    "    return train, test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First get and preprocess the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get data: contains 70k samples of which the last 10k are meant for testing\n",
    "mnist = fetch_mldata('MNIST original', data_home='./data')\n",
    "\n",
    "# prepare for concat\n",
    "y_all = mnist.target[:, np.newaxis]\n",
    "\n",
    "# intercept term to be added\n",
    "intercept = np.ones_like(y_all)\n",
    "\n",
    "# normalize the data (zero mean and unit variance)\n",
    "train_normalized, test_normalized = normalize_features(\n",
    "    mnist.data[:60000, :],\n",
    "    mnist.data[60000:, :],\n",
    ")\n",
    "\n",
    "# concat intercept, X, and y so that shuffling is easier in a next step\n",
    "train_all = np.hstack((\n",
    "    intercept[:60000],\n",
    "    train_normalized,\n",
    "    y_all[:60000],\n",
    "))\n",
    "test_all = np.hstack((\n",
    "    intercept[60000:],\n",
    "    test_normalized,\n",
    "    y_all[60000:],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I don't think this randomization step is really needed in our case, but let's stick with the ufldl tutorial here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.shuffle(train_all)\n",
    "np.random.shuffle(test_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now prepare the final train and test datasets. Let's only pick the data for the digits `0` and `1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# train data\n",
    "train_X = train_all[np.logical_or(train_all[:, -1] == 0, train_all[:, -1] == 1), :-1]\n",
    "train_y = train_all[np.logical_or(train_all[:, -1] == 0, train_all[:, -1] == 1), -1]\n",
    "\n",
    "# test data\n",
    "test_X = test_all[np.logical_or(test_all[:, -1] == 0, test_all[:, -1] == 1), :-1]    \n",
    "test_y = test_all[np.logical_or(test_all[:, -1] == 0, test_all[:, -1] == 1), -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cost_function(theta, X, y):\n",
    "    h = sigmoid(X.dot(theta))\n",
    "    return -sum(y * np.log(h) + (1 - y) * np.log(1 - h))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gradient(theta, X, y):\n",
    "    errors = sigmoid(X.dot(theta)) - y\n",
    "    return errors.dot(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization took 1.385227918624878 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fhartl/anaconda/envs/py34/lib/python3.4/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: divide by zero encountered in log\n",
      "  app.launch_new_instance()\n",
      "/Users/fhartl/anaconda/envs/py34/lib/python3.4/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: invalid value encountered in multiply\n",
      "  app.launch_new_instance()\n",
      "/Users/fhartl/anaconda/envs/py34/lib/python3.4/site-packages/IPython/kernel/__main__.py:2: RuntimeWarning: overflow encountered in exp\n",
      "  from IPython.kernel.zmq import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "J_history = []\n",
    "\n",
    "t0 = time.time()\n",
    "res = scipy.optimize.minimize(\n",
    "    fun=cost_function,\n",
    "    x0=np.random.rand(train_X.shape[1]) * 0.001,\n",
    "    args=(train_X, train_y),\n",
    "    method='L-BFGS-B',\n",
    "    jac=gradient,\n",
    "    options={'maxiter': 100, 'disp': True},\n",
    "    callback=lambda x: J_history.append(cost_function(x, train_X, train_y)),\n",
    ")\n",
    "t1 = time.time()\n",
    "\n",
    "print('Optimization took {s} seconds'.format(s=t1 - t0))\n",
    "optimal_theta = res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x109f83e80>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEPCAYAAACp/QjLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UXHWd5/H3B0JMIEHIAHkiYJYEIQZNRCMcfGgRkEQE\n",
       "3NkFdHTQlawOGhhcXQkzC60zI7gOrMCoqwHdDCM4UQ8QlaeQoQUUeUyAEKIJEJJA6CAGSCSBAN/9\n",
       "496iq7uqu6s6devWw+d1Tp2693fv/dU3fTr17ft7uooIzMzMiu2SdwBmZtZ4nBzMzKyEk4OZmZVw\n",
       "cjAzsxJODmZmVsLJwczMSmSWHCSNkHS3pOWSVkq6MC3vlLRB0rL0NbvomvmSVktaJem4rGIzM7OB\n",
       "Kct5DpJ2j4iXJA0D7gS+DHwI2BIRl/Q5dxpwNfBuYCJwK3BwRLyeWYBmZlZWps1KEfFSujkc2BXY\n",
       "nO6rzOknAddExI6IWAusAWZlGZ+ZmZWXaXKQtIuk5UA3cFtEPJIemifpQUlXStorLZsAbCi6fAPJ\n",
       "HYSZmdVZ1ncOr0fEDGB/4P2SOoDvAZOBGcBG4OKBqsgyPjMzK29YPT4kIl6Q9CvgXRHRVSiXdAXw\n",
       "i3T3KWBS0WX7p2W9SHLCMDMbgogo16RfVmbJQdI+wKsR8bykkcCxwNckjYuIZ9LTPgY8nG4vBq6W\n",
       "dAlJc9JU4J5ydVfzD7SBSeqMiM6842gF/lnWln+etVXtH9ZZ3jmMBxZK2oWk+eqqiFgq6V8lzSBp\n",
       "MnoC+BxARKyUtAhYCbwKnBleMtbMLBeZJYeIeBh4Z5nyvx7gmm8A38gqJjMzq4xnSFtX3gG0kK68\n",
       "A2gxXXkH0M4ynQSXBUnhPgczs+pU+93pOwczMyvh5GBmZiWcHMzMrISTg5mZlXByMDOzEk4OZmZW\n",
       "wsnBzMxKODmYmVkJJwczMyvRlMlBmn2TNH1O3nGYmbWqpkwOcOOH4chLnSDMzLLRpMkBYMEUOGBe\n",
       "3lGYmbWiJk4OAKNG5h2BmVkravLksHVb3hGYmbWiJk4OZzwG6y7POwozs1bUpMnhE/fC786KWHFD\n",
       "3pGYmbWiLJ8hnaGrb4/AicHMLCNNeufA7LwDMDNrZc2aHPaVODDvIMzMWlWzJoeb8d2DmVlmmjU5\n",
       "3AB4drSZWUYySw6SRki6W9JySSslXZiWj5G0RNIfJN0iaa+ia+ZLWi1plaTjBqj+FqBDYkRW8ZuZ\n",
       "tbPMkkNEbAc+GBEzgLcDH5T0XuBcYElEHAwsTfeRNA04FZgGHA98V1LZ+CJ4DlgBvC+r+M3M2lmm\n",
       "zUoR8VK6ORzYFdgMnAgsTMsXAien2ycB10TEjohYC6wBZg1QvZuWzMwykmlykLSLpOVAN3BbRDwC\n",
       "jI2I7vSUbmBsuj0B2FB0+QZg4gDV34g7pc3MMpHpJLiIeB2YIenNwM2SPtjneEiKgaooVyipEyQ4\n",
       "d39p2ccjbrymdlGbmTU/SR1Ax1Cvr8sM6Yh4QdKvgMOBbknjIuIZSeOBTelpTwGTii7bPy0rV18n\n",
       "QDrXYUxmgZuZNamI6AK6CvuSLqjm+ixHK+1TGIkkaSRwLLAMWAycnp52OnBdur0YOE3ScEmTganA\n",
       "PYN8jPsdzMwykOWdw3hgYTriaBfgqohYKmkZsEjSZ4G1wCkAEbFS0iJgJfAqcGZEDNTkBLAEuEJi\n",
       "ZARevtvMrEY0+PdvY5EUEaGefe4A/imCm3IMy8ysofX97hxMs86QLuZRS2ZmNdYKycH9DmZmNdYK\n",
       "yeFBYA+JqXkHYmbWKpo+OUQQuGnJzKymmj45pNy0ZGZWQ00/Wikp480ky22MjeCl8leambWvdhyt\n",
       "RAQvAPcDHxzsXDMzG1xLJIeU+x3MzGqklZLDDcAciYpvm8zMrLxWSg4rSJ4bcXDegZiZNbuWSQ7p\n",
       "kFaPWjIzq4GWSQ4p9zuYmdVASwxl7TnGniTPgBgfwdb6RmZm1rjacihrQQQvAvcCR+cdi5lZM2up\n",
       "5JC6ATctmZntlFZNDh7Sama2E1oxOTyavh+aaxRmZk2s5ZKDh7Same28lksOKQ9pNTPbCS01lLXn\n",
       "HEYBG4EJEWypT2RmZo2rrYeyFqRzHH4HfCjvWMzMmlFLJoeU+x3MzIaolZPDjcBsD2k1M6teZslB\n",
       "0iRJt0l6RNIKSWel5Z2SNkhalr5mF10zX9JqSaskHbeTIfwe2AFM38l6zMzazrAM694BnBMRyyWN\n",
       "Au6XtAQI4JKIuKT4ZEnTgFOBacBE4FZJB0fE60P58AhC4kaSpqWHd+YfYmbWbjK7c4iIZyJiebq9\n",
       "lWRy2sT0cLmmnpOAayJiR0SsBdYAs3YyDC+lYWY2BHXpc5D0FmAmyQgigHmSHpR0paS90rIJwIai\n",
       "yzbQk0yG6jbgcIk372Q9ZmZtJctmJQDSJqWfAWdHxFZJ3wO+nh7+B+Bi4LP9XF52EoakzqLdrojo\n",
       "Kntx8JLEb4BjgJ9XH72ZWXOS1AF0DPX6TJODpN1IvpT/LSKuA4iITUXHrwB+ke4+BUwqunz/tKxE\n",
       "RHRWEUZhSKuTg5m1jfSP5q7CvqQLqrk+y9FKAq4EVkbEt4vKxxed9jF6OosXA6dJGi5pMjAVuKcG\n",
       "oXhIq5lZlbK8czgK+CTwkKRladl5wMclzSBpMnoC+BxARKyUtAhYCbwKnBk1WNsjgtUSfwbeASzf\n",
       "2frMzNpBS66tVHoNlwEbI7gwo7DMzBqa11Yqz0tpmJlVoV3uHEYC3cCBEWzOJjIzs8blO4cyItgG\n",
       "3AEcm3csZmbNoC2SQ6qwlIaZmQ2iLZqVkuv4T8BvSR4ANKT1mszMmpWblfoRwePACyTLeJiZ2QAy\n",
       "Xz6jsVy9Erp+Ij3/FGzZDusvi1hxQ95RmZk1mrZJDtL0OdAxC34wAZiSlM49SJqOE4SZWW9t06wE\n",
       "k86Cf5nQu2zBFDhgXj7xmJk1rjZKDqNHlC8fNbK+cZiZNb42Sg5btpcv37qtvnGYmTW+NkoO6y+D\n",
       "uWt6l53xGKy7PJ94zMwaV9vMc0iunT4Hpn4Fph8F9y+FdZe7M9rM2kG1351tlRyS6xHwNHBEBE/W\n",
       "LjIzs8blSXCDiCBIHiL0nrxjMTNrVG2XHFJ3A7PyDsLMrFG1a3LwnYOZ2QDars8hqYM3A08Be0Xw\n",
       "am0iMzNrXO5zqEAELwDrgbflHYuZWSNqy+SQctOSmVk/2jk5uFPazKwf7Zwc7sHJwcysrLbskE7q\n",
       "YTiwGRgbwdadj8zMrHE1TIe0pEmSbpP0iKQVks5Ky8dIWiLpD5JukbRX0TXzJa2WtErScVnFBhDB\n",
       "K8BDwOFZfo6ZWTPKsllpB3BORLwNOAL4gqRDgXOBJRFxMLA03UfSNOBUYBpwPPBdSVk3e7lT2sys\n",
       "jMy+fCPimYhYnm5vBR4FJgInAgvT0xYCJ6fbJwHXRMSOiFgLrCH7PgF3SpuZlVFRcpC0h6RDJL1V\n",
       "0h7VfoiktwAzSb6Mx0ZEd3qoGxibbk8ANhRdtoEkmWTJndJmZmX0+wxpSaOBucBpwD4kX+QCxkp6\n",
       "DvgxsCC9K+iXpFHAz4GzI2KL1NMfEhEhaaAe8bLHJHUW7XZFRNdAMQzgMWAPifERbBxiHWZmDUdS\n",
       "B9Ax1Ov7TQ7AdcBPgI8W/aVf+NBxJM1D1wMfGiC43UgSw1URcV1a3C1pXEQ8I2k8sCktfwqYVHT5\n",
       "/mlZiYjoHCDuikUQ0ht3D9fXok4zs0aQ/tHcVdiXdEE112c2lFXJLcJC4LmIOKeo/H+nZd+UdC6w\n",
       "V0Scm3ZIX03yRT0RuBWYEn0CrNVQ1p76+BowLIK/q1WdZmaNptrvzoHuHIorHQNMBd5UKIuI2we5\n",
       "7Cjgk8BDkpalZfOBi4BFkj4LrAVOSetbKWkRsBJ4FTizb2LIyN3Al+rwOWZmTWPQOwdJc4GzSJp5\n",
       "lpMMS70rIo7OPryy8dT6zmFfYDUwJoLXa1WvmVkjyWIS3NkkTT1PRsQHSUYdvTDE+BpOBM8CfwLe\n",
       "mncsZmaNopLksD0itgFIGhERq2i9L1LPdzAzK1JJctggaW+S0UtLJC0m6StoJZ7vYGZWpKrRSum4\n",
       "2T2BmyLilayCGiSGmvY5JHVyFHBpBO+qZb1mZo2i5n0Okq4qbEdEV0QsBq4cYnyN6gFgmsSIvAMx\n",
       "M2sElTQrTS/ekTSMFlvJNIJtwCpgRt6xmJk1gn6Tg6TzJG0BDpO0pfAimdG8uG4R1s/deIVWMzOg\n",
       "snkOF0XEuXWKZ1BZ9Dkk9fIZ4JgI/qrWdZuZ5S2LeQ7nSfqUpPPTDzhAUiuO7PGzHczMUpUkh+8C\n",
       "RwKfSPe3pmWtZhWwn8Rf5B2ImVneKkkO74mIM4FtABHxJ2C3TKPKQQSvAfcB7847FjOzvFWSHF6R\n",
       "tGthR9K+0LJrELlT2syMypLD5cC1wH6SvgH8Brgw06jy45nSZmZUOENa0qH0PNRnaUQ8mmlUA8eS\n",
       "yWilpG4mkqw8u19E+afQmZk1o2q/OytNDrsC40ie/xAAEbFuqEHujCyTQ1I/G4D3R/B4Vp9hZlZv\n",
       "NX/Yj6R5wAUkk99eKzp0WPXhNYVC05KTg5m1rUqeBPe3wFsj4rmsg2kQhU7pn+QdiJlZXirpkF4H\n",
       "vJh1IA3EndJm1vb67XOQ9D/SzWnAIcAvgcIy3RERl2QfXtm4su5z2BN4Gtg7gh1ZfY6ZWT3Vss9h\n",
       "NEnn8zpgPTA8fbW0CF6UWEvSp/JAzuGYmeWi3+QQEZ0Akk6JiEXFxySdknFceSs0LTk5mFlbqqTP\n",
       "YX6FZa3E/Q5m1tb6vXOQNBuYA0yUdBlQaKsaDS3fFn83MC/vIMzM8jLQncPTwP3A9vT9fpKF6RYD\n",
       "H66kckk/lNQt6eGisk5JGyQtS1+zi47Nl7Ra0ipJxw3lH1QjK4AD085pM7O2U8nDfoZHxCsDntT/\n",
       "te8jWeL7XyPisLTsAmBL39FOkqYBV5OsijoRuBU4OCJe73NepqOVej6HO4HzI/iPrD/LzCxrNXvY\n",
       "j6RfSfqvlGl6krSHpFMl3TBQ5RFxB7C5XPVlyk4CromIHRGxFlhDvu3+7ncws7Y1ULPSZ0iGc94n\n",
       "6WFJt0hakjYR3QccCpw+xM+dJ+lBSVdK2istmwBsKDpnA8kdRF6cHMysbQ00lHUTcD5wvqSxwIHp\n",
       "oScjonsnPvN7wNfT7X8ALgY+218Y5QoldRbtdkVE107E05+7SWIzM2s6kjqAjqFeP9Bopa30/+X8\n",
       "Mkmzz99HxK3VfGCadAr1XAH8It19CphUdOr+aVm5Ojqr+cwhWgsMl5gYUT4OM7NGlf7R3FXYT/t7\n",
       "K9Zvs1JEjIqI0eVeJMt3fw64tNqAJY0v2v0YUBjJtBg4TdJwSZOBqSRNO7lIn+dwN25aMrM2VMmq\n",
       "rCUi4lXgQUmXD3SepGuADwD7SFpPsvR3h6QZJHclT5AkGSJipaRFwErgVeDMqORhE9m6h2SF1mtz\n",
       "jsPMrK4qethPI6nXUNbkszge+J8RHF2PzzMzy0omT4JrJHVODmNI+h72juj1oCMzs6ZSs3kOBhH8\n",
       "CegmWbLczKxtODkMzvMdzKztODkMrvDYUDOztuHkMDjfOZhZ23GH9KCfxwjgOWCfCLbV63PNzGrJ\n",
       "HdI1FsF2krkX78w7FjOzenFyqIyblsysrTg5VMbLaJhZW3FyqExhGQ0zs7bg5FCZPwBjJPbNOxAz\n",
       "s3pwcqhABK8D9+KmJTNrE04OlXOntJm1DSeHyrlT2szahpND5e4FZknUbQKemVlenBwqFMFG4M/A\n",
       "QXnHYmaWNSeH6rhpyczagpNDdTzfwczagpNDdXznYGZtwauyVvXZjCJ5MtzeEbySRwxmZkNR7Xfn\n",
       "sCyDaT3T3w+zX4Nn75K6n4X1l0WsuCHvqMzMas3JoULS9Dlw5KXwrdG8sXz33IOk6ThBmFmrcZ9D\n",
       "xSadBQum9C5bMAUOmJdPPGZm2ck0OUj6oaRuSQ8XlY2RtETSHyTdImmvomPzJa2WtErScVnGVr3R\n",
       "I8qXjxpZ3zjMzLKX9Z3Dj4Dj+5SdCyyJiIOBpek+kqYBpwLT0mu+K6mB7my2bC9fvtWPDjWzlpPp\n",
       "l29E3AFs7lN8IrAw3V4InJxunwRcExE7ImItsIaGGja6/jKYu6Z32ZlPw7rL84nHzCw7eXRIj42I\n",
       "7nS7Gxibbk8Afld03gZgYj0DG0jEihuk6cCceUlT0vA94NOT4Lu35x2bmVmt5TpaKSJC0kATLcoe\n",
       "k9RZtNsVEV21jKvfYJJRSW+MTJL4f8A3gLPq8flmZpWS1AF0DPX6PJJDt6RxEfGMpPHAprT8KWBS\n",
       "0Xn7p2UlIqIz2xAr9iXgYYlFEdyZdzBmZgXpH81dhX1JF1RzfR4dvouB09Pt04HrispPkzRc0mRg\n",
       "KslaRg0rgj8BXwCulPCoJTNrGZkunyHpGuADwD4k/QvnA9cDi4ADgLXAKRHxfHr+ecB/A14Fzo6I\n",
       "m8vUmdvyGf2R+HdgbQRfzTsWM7Nyqv3u9NpKNSCxH/AQ8NEI7s07HjOzvqr97mygeQTNK4JNwDnA\n",
       "DyWG5x2PmdnOcnKonZ8AjwPn5R2ImdnOcrNSDUlMAJYDx0TwUN7xmJkVuFkpRxE8TbIcyI8kr3hr\n",
       "Zs3LyaH2fgQ8B3w570DMzIbKzUoZkDgQuA94XwSr8o7HzMzNSg0ggieBTpLRS7vmHI6ZWdWcHLLz\n",
       "PZLJfF/MOxAzs2q5WSlDElOBu4BZETyedzxm1r7crNRAIlgNXAhcIdEUCc3MDJwc6uHbwB7A3LwD\n",
       "MTOrlJuV6kDibSRL574zgvU5h2NmbcjNSg0ogkeAS4Hvu3nJzJqBk0P9fBMYD3wy70DMzAbjZqU6\n",
       "kpgJS/8DLlsObxJs2Q7rL0sfP2pmlplqvzu9/k9dTR8PxwZc39FTNvcgaTpOEGbWSNysVFeTzoL/\n",
       "s3fvsgVT4IB5+cRjZlaek0NdjR5RvnyUnz9tZg3FyaGutmwvX/7nfsrNzPLh5FBX6y+DuWt6l335\n",
       "Jfj8WIm9y19jZlZ/Hq1UZ9L0OUkfw6iRsHUbbPwOLDsamAOcEMGaweowM6tWtd+dTg4NQuJzwNeA\n",
       "UyK4Pe94zKy1eIZ0k4rg+yQT5H4q8emcwzGzNpfbnYOktcCLwGvAjoiYJWkM8O/AgcBa4JSIeL7P\n",
       "dS1551AgcQjwS+CnwN9F8HrOIZlZC2imO4cAOiJiZkTMSsvOBZZExMHA0nS/raSPFT0COIrkLmKP\n",
       "nEMyszaUd7NS3yx2IrAw3V4InFzfcBpDBH8EjgW2AL+WmJBzSGbWZvK+c7hV0n2SCs86GBsR3el2\n",
       "NzA2n9DyF8HLwGeAnwF3S7wz55DMrI3kubbSURGxUdK+wBJJq4oPRkRIKtshIqmzaLcrIrqyCzM/\n",
       "EQRwkcRq4GaJuRFcl3dcZtb4JHUAHUO+vhGGskq6ANhK8rS0joh4RtJ44LaIOKTPuS3dId0fiXcB\n",
       "1wGXwvRHknWaRo/wyq5mVommWJVV0u7ArhGxRdIewHEkY/wXA6eTPPvgdPBfyQUR3CdxBNz0a/jw\n",
       "PnDxnj1HvbKrmdVWLncOkiYD16a7w4AfR8SF6VDWRcABtOlQ1sFIJyyBXx5TemTOTRE3zK5/RGbW\n",
       "DJriziEingBmlCn/E1Dmi8967L5b+XKv7GpmtZP3UFarWn8ru05+t8Q5XsDPzGrByaHplFvZ9YzH\n",
       "IM4HDgcel7gieSSpmdnQNMRopWq0e58DlFvZdd3lhc5oif2AzwKfB54CvgP8LJ03YWZtyquyGgAS\n",
       "w4CPAF8A3g5cCXw/gnXJ8elzPBzWrH04OVgJibcCfwN8CrgDLrgXnv508vzqgrlr4K6znSDMWpOT\n",
       "g/UrXcTvr+Ar/wzfGl16hofDmrWqZlqV1eosgj9H8AN48oHyZ4wdK7FrfaMys0bk5NCW+hsOu+8U\n",
       "oFvixxKfkBhT17DMrGE4ObSl/obD3nAayeTEXwOnAmsl7pA4V+IwqWeJdWn6HGn2TdIpXcn79Dn1\n",
       "/BeYWbbc59CmBhoO23MOI0hWdfwIcALJHxO/gn/shnWfgh8c1HO2O7TNGpk7pC0T6V3DocBH4Evz\n",
       "4ZIyM7HdoW3WqJpibSVrPumzJVYCK6UNHwE+UHrWu46R+E3hvPT1KLA+vb4Xz7Uwa1xODjYE/XVo\n",
       "P3g78HVgWvo6IX0fJfEovZLG6fvBkfP7zLWoeulxJxizbDg52BCsvwzmHtT7i/2Mx+CxiyP4NUmH\n",
       "9hvSxQAPpSdpHA3jPgDfHNG73gVT4NMXSrwMPA1sBF4od9eR1Dt9Dhx5qROMWe05OVjVIlbcIE0H\n",
       "5gzYod1zPpuB36YvAKQnuijbNPUX44C/Ayakr2EST0Ov18bk/bCzeycGSPbnzAMq+nKvRYJxcrFW\n",
       "5ORgQ5J++e3EF2B/TVOPPhDBG53aEqOA8fQki8JrJhx4cPk6Dj9a4nFgW/ra3v/7ez8K/3dy7+sX\n",
       "TIG//KrEncCW/u5ckvhqc/di1micHCwn/TVNrbu8+KwItgKr01cv0oM3AR8urXvFb4EzgBHAyPQ1\n",
       "os97ut3ff4FDZ5Hcpewm8Sz095r5Nzt792LWiJwcLBfVNk2V11+CWf2tCB6rpAbpseOByaVHHuiK\n",
       "YLbESGDffl6TYdyk8jXvvVfl/w6zxuPkYLnZ2aapbBNMcgcTwTZgXfoqIa3Yn7J3LwceLvEAcD2w\n",
       "GFg+UPOUWaPxJDhre5XMFh/42r59Dmc8BsvOgftfAE5KX8NJksRioCuCV2r/LzHrn2dIm9XZYMml\n",
       "aHb5ienrUOAWkruKGyPY7BFPljUnB7MGJzGOZL2qk4AO+OUTsGQCXLpPz1nVr1XlBGMDafrkIOl4\n",
       "4NvArsAVEfHNPsedHKxlSOwOn7oTrppZevRvu+Hb1wKbKD9a6o8RvJrUU3ZIbVUJxsmltTX12kqS\n",
       "dgX+BTgGeAq4V9LiiHg038hal6SOiOjKO45WMJSfZQQvSS+/WP7otueAFcB+JDPL+46WGiPxIvAs\n",
       "fHg/uLjPCKkFU+ATnRIbSBLMG8mkNPbGmW3eU8eL42DPZ5yk8tFQyQGYBayJiLUAkn5Ccuvt5JCd\n",
       "DqAr5xhaRQdD+ln2NyFw/boIvtPfVelT+8YA+8LmH5M8i6OPSVOBq0kSzBiJF0juOjb1fn/Pfyk/\n",
       "X+PkL0vcBbwCvNxfckniqdVs80IdnUDnO/JIUrVNcvnWMVSNlhwmAuuL9jcA78kpFrM6qWxCYF8R\n",
       "vEbaxCRt7C5/1sO/K8w4751M2K/3+55vLn/9Ye8FniAZbfUmJY0SrwAvp6+i7eP3h3/es/f1C6bA\n",
       "WQskfkHpDPUys9bf3lk+Sf3lV9PFG18DXk/f+77S8ncfC0deMtQkVfskl18dO6PRkkNjdYCY1UE9\n",
       "5mskn9OTTEhWx32DtOoIkmVJ+rh/aZ/lTIaRJor0Nbzn/Y8LgcNL63hlC/AgvWep70vprPURMGlK\n",
       "6fUAh8wCbiPpi9wlfe/7SsuP2w3+qU/b+oIp8L9+JfF64cdRdDB6v580rPz1f/+LdFHIKHNdn+2P\n",
       "7gEXDiut47zr0qbAcnX0ieuEveCiN5XWUZ/Z9w3VIS3pCKAzIo5P9+cDrxd3SktqnIDNzJpI045W\n",
       "kjQM+D3wIZJ1be4BPu4OaTOz+mqoZqWIeFXSF4GbSW4Tr3RiMDOrv4a6czAzs8awS94BVErS8ZJW\n",
       "SVot6at5x9PsJK2V9JCkZZLuyTueZiPph5K6JT1cVDZG0hJJf5B0iySvzFqhfn6enZI2pL+jy9IJ\n",
       "sjYISZMk3SbpEUkrJJ2Vllf1+9kUyaFoctzxJJOBPi7p0HyjanoBdETEzIiYlXcwTehHJL+Pxc4F\n",
       "lkTEwcDSdN8qU+7nGcAl6e/ozIi4KYe4mtEO4JyIeBtwBPCF9Puyqt/PpkgOFE2Oi4gdQGFynO0c\n",
       "L0MyRBFxB7C5T/GJwMJ0eyFwcl2DamL9/DzBv6NVi4hnImJ5ur2VZBLxRKr8/WyW5FBuctzEnGJp\n",
       "FQHcKuk+SXPzDqZFjI2IwmS0bmBsnsG0iHmSHpR0pZvpqifpLcBM4G6q/P1sluTgXvPaOyoiZgKz\n",
       "SW4735d3QK0kkpEe/r3dOd8jeUrfDGAjcHG+4TQXSaOAnwNnR8SW4mOV/H42S3J4Cih+HOMkkrsH\n",
       "G6KI2Ji+PwtcS9J0ZzunW9I4AEnjSdYtsiGKiE2RAq7Av6MVk7QbSWK4KiKuS4ur+v1sluRwHzBV\n",
       "0lskDQdOJXmilg2BpN0ljU639wCOAx4e+CqrwGLg9HT7dOC6Ac61QaRfYAUfw7+jFZEk4EpgZUR8\n",
       "u+hQVb+fTTPPQdJsep7zcGVEXJhzSE1L0mSSuwVIJkL+2D/P6ki6BvgAsA9J++35JE92WwQcAKwF\n",
       "TomI5/OKsZmU+XleQLLK7QyS5o8ngM8VtZlbPyS9F7gdeIiepqP5JCtOVPz72TTJwczM6qdZmpXM\n",
       "zKyOnBzMzKyEk4OZmZVwcjAzsxJODmZmVsLJwczMSjg5WFuRtDV9P1DSx2tc93l99n9Ty/rN6snJ\n",
       "wdpNYWLPZOAT1VyYPsZ2IPN7fVDEUdXUb9ZInBysXV0EvC99iMzZknaR9C1J96SrgP53AEkdku6Q\n",
       "dD2wIi27Ll3NdkVhRVtJFwEj0/quSssKdylK6344fcDSKUV1d0n6qaRHJf1bIThJF6UPa3lQ0rfq\n",
       "+pMxo8GeIW1WR18FvhwRHwVIk8HzETFL0puAOyXdkp47E3hbRDyZ7n8mIjZLGgncI+lnEXGupC+k\n",
       "K90WFO5S/jPwDuDtwL7AvZJuT4/NIHmA1UbgN5KOAlYBJ0fEIWlse2bw7zcbkO8crF31fYjMccBf\n",
       "S1oG/A4YA0xJj91TlBgAzpa0HLiLZIXgqYN81nuBq9MFRjcBvwbeTZI87omIp9OVR5cDBwLPA9vT\n",
       "Zxh8DNg25H+l2RA5OZj1+GLRIykPiohb0/I/F06Q1AF8CDgiImYAy4ARg9QblCajwl3Fy0VlrwG7\n",
       "RcRrJMtT/ww4AfDjMa3unBysXW0BRhft3wycWeh0lnSwpN3LXLcnsDkitks6hOQZvQU7+um0vgM4\n",
       "Ne3X2Bd4P8kKmWUfgZkuo75XRNwIfImkScqsrtznYO2m8Bf7g8BrafPQj4DLgLcAD6Tr4W8ieYZA\n",
       "3ydm3QR6/VU+AAAAcklEQVR8XtJK4PckTUsFPwAeknR/RHyqcF1EXCvpyPQzA/hKRGxKH/red1nk\n",
       "IEla10saQZJAzqnJv9ysCl6y28zMSrhZyczMSjg5mJlZCScHMzMr4eRgZmYlnBzMzKyEk4OZmZVw\n",
       "cjAzsxJODmZmVuL/A/+0Tfs1CW9cAAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f34b898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(J_history, marker='o')\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('J(theta)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def accuracy(theta, X, y):\n",
    "    correct = np.sum(np.equal(y, (sigmoid(X.dot(theta))) > 0.5))\n",
    "    return correct / y.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training accuracy: 1.0\n",
      "Test accuracy: 0.9995271867612293\n"
     ]
    }
   ],
   "source": [
    "print('Training accuracy: {acc}'.format(acc=accuracy(res.x, train_X, train_y)))\n",
    "print('Test accuracy: {acc}'.format(acc=accuracy(res.x, test_X, test_y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looking good, right? Well, look closer...\n",
    "I actually had to use the **`L-BFGS-B`** optimization method for it to work.<br>\n",
    "Had I used the expected **`BFGS`** method, **nan** and **inf** values due to `log(0)` would have made trouble.<br>\n",
    "Why? I can think of two reasons:\n",
    "1. Even if being multiplied with `0`, the `log(0)` expression is still evaluated by numpy. And unfortunately: `0 * np.nan = np.nan`.\n",
    "2. Floating point arithmetic limits which don't exist in Mathematics.\n",
    "\n",
    "One way to counteract those issues is to substitute the troubling values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def safe_log(x, nan_substitute=-1e+4):\n",
    "    l = np.log(x)\n",
    "    l[np.logical_or(np.isnan(l), np.isinf(l))] = nan_substitute\n",
    "    return l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cost_function_safe(theta, X, y):\n",
    "    h = sigmoid(X.dot(theta))\n",
    "    return -sum(y * safe_log(h) + (1 - y) * safe_log(1 - h))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: 0.001535\n",
      "         Iterations: 100\n",
      "         Function evaluations: 114\n",
      "         Gradient evaluations: 114\n",
      "Optimization took 5.0561230182647705 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fhartl/anaconda/envs/py34/lib/python3.4/site-packages/IPython/kernel/__main__.py:2: RuntimeWarning: overflow encountered in exp\n",
      "  from IPython.kernel.zmq import kernelapp as app\n",
      "/Users/fhartl/anaconda/envs/py34/lib/python3.4/site-packages/IPython/kernel/__main__.py:2: RuntimeWarning: divide by zero encountered in log\n",
      "  from IPython.kernel.zmq import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "J_history = []\n",
    "\n",
    "t0 = time.time()\n",
    "res = scipy.optimize.minimize(\n",
    "    fun=cost_function_safe,\n",
    "    x0=np.random.rand(train_X.shape[1]) * 0.001,\n",
    "    args=(train_X, train_y),\n",
    "    method='BFGS',\n",
    "    jac=gradient,\n",
    "    options={'maxiter': 100, 'disp': True},\n",
    "    callback=lambda x: J_history.append(cost_function_safe(x, train_X, train_y)),\n",
    ")\n",
    "t1 = time.time()\n",
    "\n",
    "print('Optimization took {s} seconds'.format(s=t1 - t0))\n",
    "optimal_theta = res.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "Notice that the above optimization procedure doesn't converge due to the substitutions (which doesn't allow the gradients to further improve (= get smaller) at some point). Therefore, it used all 100 allowed iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11e1409b0>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAH7lJREFUeJzt3Xu4XXVh5vHvSzACBo2p9iQkoUSSCCG0opVLBYkWQ5LB\n",
       "gNYJ2KedFCHzzIMKOq1jwnSatM8UwV68MMPMSMFJmZJOaitGhZDIeIq1YryAXGIkoaYQNCdeUIkC\n",
       "JvDOH2tts3PYZ5+Tk7PXPnvv9/M85znrttf67ZUn6z2/y1pLtomIiBjKEe0uQEREjG8JioiIaCpB\n",
       "ERERTSUoIiKiqQRFREQ0laCIiIimWhoUkq6S9ICkByVdVS6bImmzpIclbZI0uW77VZK2S9omaWEr\n",
       "yxYRESPTsqCQNB+4HHgt8GvABZJOBFYCm23PBe4q55E0D7gYmAcsAm6QlBpPRESbtfJCfBLwZdtP\n",
       "234W+Efgt4ClwNpym7XAReX0hcA62/ts7wR2AKe3sHwRETECrQyKB4FzyqamY4AlwAygz/ZAuc0A\n",
       "0FdOHwfsqvv8LmB6C8sXEREjcGSrdmx7m6TrgE3AT4H7gGcHbWNJzZ4hkueLRES0WcuCAsD2zcDN\n",
       "AJL+lKKWMCBpqu3dkqYBe8rNHwdm1n18RrnsIMMES0REDMG2RvO5lgaFpF+2vUfS8cBbgTOBWcBy\n",
       "4Lry923l5huAWyX9JUWT0xxgS6P9jvbLdhtJa2yvaXc5xoOciwNyLg7IuTjgcP7IbmlQAJ+Q9EvA\n",
       "PuAK2z+WdC2wXtJlwE5gGYDtrZLWA1uB/eX2qT1ERLRZq5ueXt9g2Q+B84bY/hrgmlaWKSIiDk3u\n",
       "U+hs/e0uwDjS3+4CjCP97S7AONLf7gJ0A3Va644kp48iIuLQHM61MzWKiIhoKkERERFNJSgiIqKp\n",
       "BEVERDSVoIiIiKYSFBER0VSCIiIimkpQREREUwmKiIhoKkERERFNdWRQSIs3SvOXtLscERG9oCOD\n",
       "Au44H876SMIiIqL1OjQoAG6cDce/u92liIjodh0cFACTjm53CSIiul1Lg0LSKkkPSXpA0q2SXihp\n",
       "iqTNkh6WtEnS5EHbb5e0TdLC4Y+w96lWlj8iIloYFJJOAFYAr7Z9KjABuARYCWy2PRe4q5xH0jzg\n",
       "YmAesAi4QVKT8l3+CDx6favKHxERhVbWKH5C8a7sYyQdCRwDfAdYCqwtt1kLXFROXwiss73P9k5g\n",
       "B3B6410v2Qj3XGk/eHvLSh8REUAL35lt+4eS/gJ4FHgKuNP2Zkl9tgfKzQaAvnL6OOCeul3sAqY3\n",
       "3vfti1tU7IiIGKRlQSHpROA9wAnAj4G/k/Q79dvYtqRm72JtuE7SmrrZftv9h1XYiIguI2kBsGAs\n",
       "9tWyoAB+Hfhn2z8AkPQPwFnAbklTbe+WNA3YU27/ODCz7vMzymXPY3tNy0odEdEFyj+g+2vzklaP\n",
       "dl+t7KPYBpwp6WhJAs4DtgKfBpaX2ywHbiunNwCXSJooaRYwB9jSwvJFRMQItLKP4huS/hr4KvAc\n",
       "8HXgY8CxwHpJlwE7gWXl9lslracIk/3AFbabNUtFREQF1GnXYkm2rXaXIyKikxzOtbPD78yOiIhW\n",
       "S1BERERTCYqIiGgqQREREU0lKCIioqmODAqppTcKRkREnY4MCmBiuwsQEdErOjUoXtjuAkRE9IpO\n",
       "DYrUKCIiKpKgiIiIpjo1KNL0FBFRkU4NitQoIiIq0qlBkRpFRERFOjUoUqOIiKhIpwZFahQRERXp\n",
       "1KBIjSIioiItDQpJr5R0b93PjyVdKWmKpM2SHpa0SdLkus+skrRd0jZJC4fYdYIiIqIiLQ0K29+y\n",
       "fZrt04DXAD8DPgmsBDbbngvcVc4jaR5wMTAPWATcIKlRGdP0FBFRkSqbns4Ddth+DFgKrC2XrwUu\n",
       "KqcvBNbZ3md7J7ADOL3BvlKjiIioSJVBcQmwrpzusz1QTg8AfeX0ccCuus/sAqY32FdqFBERFank\n",
       "cd2SJgJvBt4/eJ1tS3KTjzdYt+ht0p2vKGf6bfePQTEjIrqGpAXAgrHYV1XvdVgMfM3298r5AUlT\n",
       "be+WNA3YUy5/HJhZ97kZ5bJBNt5p899bWN6IiI5W/gHdX5uXtHq0+6qq6entHGh2AtgALC+nlwO3\n",
       "1S2/RNJESbOAOcCWBvtLH0VEREVaXqOQ9CKKjuwVdYuvBdZLugzYCSwDsL1V0npgK7AfuMJ2o2ap\n",
       "BEVEREXU+Do8fhX9GV5t8yftLktERKeQZNsazWdzZ3ZERDTVqUGR4bERERXp1KBIjSIioiIJioiI\n",
       "aKpTgyJNTxERFenUoEiNIiKiIp0aFKlRRERUpFODIjWKiIiKdGpQpEYREVGRTg2K1CgiIiqSoIiI\n",
       "iKY6NSjS9BQRUZFODYrUKCIiKtKpQZEaRURERTo1KFKjiIioSKcGRWoUEREVaXlQSJos6ROSvilp\n",
       "q6QzJE2RtFnSw5I2SZpct/0qSdslbZO0cIjdpkYREVGRKmoUHwFut30y8KvANmAlsNn2XOCuch5J\n",
       "84CLgXnAIuAGSY3KmKCIiKhIS4NC0kuAc2zfDGB7v+0fA0uBteVma4GLyukLgXW299neCewATm+w\n",
       "6zQ9RURUpNU1ilnA9yR9XNLXJd0o6UVAn+2BcpsBoK+cPg7YVff5XcD0BvudKDGqd79GRMShObKC\n",
       "/b8aeJftr0j6MGUzU41tS3KTfTRYtxq45o+l/c8B/bb7x6zEERFdQNICYMGY7Mtudo0+zJ1LU4Ev\n",
       "2Z5Vzp8NrAJeAbzB9m5J04DP2z5J0koA29eW228EVtv+ct0+DX4KeLnNT1tW+IiILiLJtkfVEtPS\n",
       "pifbu4HHJM0tF50HPAR8GlheLlsO3FZObwAukTRR0ixgDrClwa6fIR3aERGVaHXTE8C7gb+RNBF4\n",
       "BLgUmACsl3QZsBNYBmB7q6T1wFZgP3CFG1d5fk46tCMiKtHSpqdWKJueHgPOtnm03eWJiOgE47bp\n",
       "qYVSo4iIqEinBkX6KCIiKtKpQZEaRURERTo5KFKjiIioQKcGRZqeIiIq0qlBkaaniIiKdGpQpEYR\n",
       "EVGRTg2K1CgiIirSyUGRGkVERAU6NSieITWKiIhKdGpQpEYREVGRTg2KdGZHRFSkU4MindkRERXp\n",
       "1KBIjSIioiKdGhSpUUREVGREQSHpRZJOkvRKSS86lANI2inpfkn3StpSLpsiabOkhyVtkjS5bvtV\n",
       "krZL2iZp4RC7TWd2RERFhnzDnaRjgRXAJcDLgAFAQJ+kHwB/A9xoe+8wxzCwwPYP65atBDbb/qCk\n",
       "95fzKyXNAy4G5gHTgc9Jmmv7uUH7fAZ4yUi/ZEREjF6zGsVtwJPAm22/wvZZts+0PQu4APgp8KkR\n",
       "HmfwW5WWAmvL6bXAReX0hcA62/ts7wR2AKc32F+aniIiKjJkjcL2bzZZtxv4WPkzHFPUDJ4F/pft\n",
       "G4E+2wPl+gGgr5w+Drin7rO7KGoWg6UzOyKiIkMGRT1JU4A51P0Vb/vuER7jdba/K+nlwGZJ2+pX\n",
       "2nbxHuwhNVqXGkVEREWGDQpJK4ArgRnAfcCZwJeAN47kALa/W/7+nqRPUjQlDUiaanu3pGnAnnLz\n",
       "x4GZdR+fUS4b5LwL4ITZ0k1rgH7b/SMpS0REr5C0AFgwJvuym/0xD5IeBF4LfMn2qySdBHzA9luG\n",
       "3bl0DDDB9pPlaKlNwB8D5wE/sH2dpJXAZNu1zuxbKcJkOvA5YLbrClnUPrwMuNjmbaP50hERvUaS\n",
       "bQ/uLx6RkTQ9PW37KUlIOsr2NkmvHOH++4BPSqod629sb5L0VWC9pMuAncAyANtbJa0HtgL7gSvc\n",
       "OMkyPDYioiIjCYpdkl5KMQpqs6QnKC7uw7L9beBVDZb/kKJW0egz1wDXDLPrdGZHRFRk2KCwXRu6\n",
       "ukZSP/BiYGMrCzUC6cyOiKjIsHdmS7qlNm273/YG4KaWlmp4qVFERFRkJI/wmF8/I+lI4DWtKc6I\n",
       "pUYREVGRIYNC0tWSngROlfRk7YdiKOuGykrYWDqzIyIqMpLhsdfaXllReYZVDo89GbjN5qR2lyci\n",
       "ohMczvDYkTQ9XS3pdyX9UXmw4yU1ev5SldL0FBFRkZEExQ3AWcBvl/N7y2XtlM7siIiKjOQ+ijNs\n",
       "nybpXijugZD0ghaXazipUUREVGQkNYqfS5pQmykf7jf4/RBVS2d2RERFRhIU1wOfBH5Z0jXAF4EP\n",
       "tLRUw3uG1CgiIiox7KgnAEknA7X3U9xl+5stLVXzshg8AXgWOMJu+BjyiIioczijnkYaFBOAqRR9\n",
       "Ggaw/ehoDni4al9WYh9wjM2+dpQjIqKTtPTpsZLeDaymuNHu2bpVp47mgGOo1qGdoIiIaKGRjHp6\n",
       "D/BK2z9odWEOUYbIRkRUYCSd2Y8CP2l1QUYhQ2QjIiowZI1C0u+Xk/8C9Ev6DMXFGYpXXf9lqws3\n",
       "jAyRjYioQLMaxbHAJIoaxWaKi/Kk8ufYkR5A0gRJ90r6dDk/RdJmSQ9L2iRpct22qyRtl7RN0sJh\n",
       "dp0hshERFRiyRmF7DYCkZbbX16+TtOwQjnEVxatNa+GyEths+4OS3l/O196XfTEwj/J92ZLm2h7q\n",
       "5r7UKCIiKjCSPopVI1z2PJJmAEuAvwJqw7KWAmvL6bVA7Q16FwLrbO+zvRPYATR7+GA6syMiKtCs\n",
       "j2IxxUV+uqSPcuBCfywjH5L6IeB9FK9PremzPVBODwB95fRxwD112+2iqFkMJZ3ZEREVaFaj+A7w\n",
       "NeDp8vfXgK9SvLTo/OF2LOkCYI/tezkQMgdxcbdfszv+mq1LjSIiogLN+ii+AXxD0q22fz7Udk38\n",
       "BrBU0hLgKODF5fu3ByRNtb1b0jSKG/kAHgdm1n1+RrnseSStgXfOgK9cLm2R7f5RlC8iomtJWgAs\n",
       "GJN9DfUID0mfBf438FnbPxu07kXABcBy20uGPYh0LvAHtt8s6YPAD2xfJ2klMNl2rTP7Vop+ienA\n",
       "54DZHlTAukd4fBb4HzafOcTvHBHRc1r1CI9LgXcBfyzpWeC7FE1ItWc+/V9g+SEcq3bBvxZYL+ky\n",
       "YCewDMD2VknrKUZI7QeuGBwSg6TpKSKiAiN9KGAf8Cvl7L/WdUZXrq5G8bfAp2zWtassERGdoiU1\n",
       "Ckl7GaIzWdIzFMNX/9D250Zz4DGQGkVERAWadWZPGmqdpCOBUyj6FE5pQblGIsNjIyIqMJIb7p7H\n",
       "9v5yVNT1Y1yeQ5E7syMiKjCqoKix/T/HqiCjkGc9RURU4LCCos1So4iIqEAnB0U6syMiKtDJQZHO\n",
       "7IiICnRyUKRGERFRgU4OitQoIiIq0OlBkRpFRESLdXJQZHhsREQFOjkoUqOIiKhAJwdFOrMjIirQ\n",
       "yUGRzuyIiAp0elCkRhER0WKdHBTpzI6IqEDLgkLSUZK+LOk+SVslfaBcPkXSZkkPS9okaXLdZ1ZJ\n",
       "2i5pm6SFwxwiNYqIiAq0LChsPw28wfargF8F3iDpbGAlsNn2XOCucp7yndkXA/OARcANkpqVL53Z\n",
       "EREVaGnTk+2flZMTgQnAE8BSYG25fC1wUTl9IbDO9j7bOyneoHd6k92nMzsiogItDQpJR0i6DxgA\n",
       "Pm/7IaCv7p3bA0BfOX0csKvu47uA6U12nxpFREQFhnwV6liw/RzwKkkvAe6U9IZB6y2p4Xu5a5s0\n",
       "WihpDcycDL/dJ123wHb/mBU6IqILSFoALBiLfbU0KGps/1jSZ4HXAAOSptreLWkasKfc7HFgZt3H\n",
       "ZpTLGu1vjcR0YJl9bX8Lix4R0ZHKP6D7a/OSVo92X60c9fSy2ogmSUcDbwLuBTYAy8vNlgO3ldMb\n",
       "gEskTZQ0C5gDbGm87/lLYOkt8F9eLi3eWMxHREQrtLJGMQ1YW45cOgK4xfZdku4F1ku6DNgJLAOw\n",
       "vVXSemArsB+4wvYQzVJnfQRunF3OnA8rTpTmYz94ewu/T0RET9KQ1+JxqujTaFTmJRvt2xdXXqCI\n",
       "iA4gybY1ms928p3Zg0w6ut0liIjoRl0UFHufancJIiK6UYcGxYodB89f/gg8en17yhIR0d06tI/i\n",
       "lH8Dx78b5syDp5+BL74nHdkREUM7nD6KjgyK2peVeAfwepvfa2+pIiLGt17uzH4EmD3sVhERMWqd\n",
       "HhQ7gBPbXYiIiG7W6UHxXeAlEpPaXZCIiG7V0UFh8xzwL8Ar2l2WiIhu1dFBUUo/RUREC3VDUKSf\n",
       "IiKihbohKFKjiIhooW4IitQoIiJaqBuCIjWKiIgW6ug7s4t5jgT2Ai+xeaZ9JYuIGL96+c5sbPYD\n",
       "u4AT2lyUiIiu1NKgkDRT0uclPSTpQUlXlsunSNos6WFJm2qvTC3XrZK0XdI2SQtHeKj0U0REtEir\n",
       "axT7gPfaPgU4E3inpJOBlcBm23OBu8p5JM0DLgbmAYuAG8pXqQ4n/RQRES3S0qCwvdv2feX0XuCb\n",
       "wHRgKbC23GwtcFE5fSGwzvY+2zspagqnj+BQj5AaRURES1TWRyHpBOA04MtAn+2BctUA0FdOH0fR\n",
       "31CziyJYhrOD1CgiIlriyCoOImkS8PfAVbaflA50vNt28TKiIT1vnaQ1dbP94NQoIiLqSFoALBiL\n",
       "fbU8KCS9gCIkbrF9W7l4QNJU27slTQP2lMsfB2bWfXxGuewgttccfIyzL4Jz50jb++HJp+Gxj+aN\n",
       "dxHRy2z3A/21eUmrR7uvVo96EnATsNX2h+tWbQCWl9PLgdvqll8iaaKkWcAcYEvzY8xfAif/Gfzp\n",
       "EbD+XLjjfDjrI8XyiIg4XC294U7S2cDdwP0caEJaRXHxXw8cD+wEltn+UfmZq4F3APspmqruHLTP\n",
       "QTfcLd5YhMNgSzbaty8e468UEdGRDueGu5Y2Pdn+J4autZw3xGeuAa4Z+VGOParx8klHj3wfEREx\n",
       "lI6/M7vok2hk71PVliMiojt1QVA89lFYsePgZZc/Ao9e357yRER0l45/KGCxbP4SOP7d8Jrz4IEv\n",
       "wI4/z6iniIgDDqePoiuC4sA6bgc+Zv9iFFVERNDjT48d5MsUz5SKiIgx0m1BcQ8JioiIMdVtTU8v\n",
       "BR4FXlq+pyIiIkjT0y/YPEHxIMFT2l2WiIhu0VVBUUo/RUTEGOrGoEg/RUTEGOrWoDij3YWIiOgW\n",
       "lbyPolqvmQnnz5W+/U/wo7155HhExOHpqqAo7tA+6y/hmgnA64qlK06U5pOwiIgYnS5repp5Jdw4\n",
       "6JWoN84uHu8RERGj0WVBkUeOR0SMtS4LijxyPCJirLX6Vag3SxqQ9EDdsimSNkt6WNImSZPr1q2S\n",
       "tF3SNkkLD/2IjR45viKPHI+IOAytfhXqOcBe4K9tn1ou+yDwfdsflPR+4KW2V0qaB9wKvBaYDnwO\n",
       "mGv7uUH7bHob+oFHjk86GmbMh1ess9+VPoqI6Gnj+jHjkk4APl0XFNuAc20PSJoK9Ns+SdIq4Dnb\n",
       "15XbbQTW2L5n0P5G/GUlXgN3bYIPfR2OeUHRNJXhshHRe8btO7OH0Gd7oJweAPrK6eMobpar2UVR\n",
       "szgM8/tg0UT4TN37uTNcNiLiULT1PgrbltSsStNwnaQ1dbP9tvsbf3zmlfDnkw5eduNsWPJuIEER\n",
       "EV1L0gJgwVjsqx1BMSBpqu3dkqYBe8rljwMz67abUS57HttrRnaoDJeNiN5U/gHdX5uXtHq0+2rH\n",
       "8NgNwPJyejn84rWlG4BLJE2UNAuYA2w5vENluGxExOFq9fDYdcA/A6+U9JikS4FrgTdJehh4YzmP\n",
       "7a3AemArcAdwhQ+7p73RcNnLM1w2IuIQdNUb7hpvXxsuO20avGw23HGxff9nW1nGiIjxZlwPjx1r\n",
       "o/2yEoI7vgW37oVnfpKhshHRSzpteGybzF8MbzgWbplzYFmGykZEDKfLnvXUzMwr4fqpB+bvBvpm\n",
       "wytvkRZvLJqoIiJisB6qUdQPlb0buBP4U4ApwPmpXURENNZDNYr6obKbKEOiTt5bERHRSA8FRf1Q\n",
       "2UYVqbuBI86QlvWnKSoi4oCeGfVUfLY2VPaIM+AzLz2w5qCmqHL+L34GEx6Bp76T0VER0ekyPPaQ\n",
       "9zF/CZz1kQOvTf1D4L+WaxMaEdF9Mjz2ENkP3i7Np3g44KSj4dlTKTq1Obj/ohYanzoGOLX4Sad3\n",
       "RPSWngwKoHahvx1AWrwROL9YU39KGnV6/+5s2LVWWvZQbtqLiF7Qs0FxsMc+CitOLJqi9tctH3x6\n",
       "ajWMO14GnFs2S50jvTXNUhHRtRIUDG6K2j8d3nEi3HzMwaEBaZaKiF7Uk53Zwx+jNjpq/3SYUYYG\n",
       "wJryBw7uAIciODYB3/wh/Owr8NiXYOZZxY1+aaKKiPbKqKeWHq8WGpOOhidPKZudODg0Go2UumUf\n",
       "3PiCA/MZORUR7ZOgqOzY9cNq62sUg2sXhzTcNjWPiGi5rhoeK2kR8GFgAvBXtq9rc5F+Yei+jMGn\n",
       "caiRU4P7Ne4+FW5546CaR9k5/v2fw1GCyU+WAZJAiYi2GFdBIWkC8N+A8yjel/0VSRtsf7O9JTvg\n",
       "4GG185cUoXHEGUDdnd5DjZwaPNx2EweHRC1E7j61QVNWg0D5jd1w5BNDB8ruSQfW1U83264zQ0jS\n",
       "gvIdwT0v5+KAnIuxMa6CAjgd2GF7J4CkvwUuBMZNUNSrhUYRGCvq7vReCKwo+yiaDbcdKkRGGiiX\n",
       "vgKOo3Gg1Dd5NexDabBdzVtPlRbufn7AjDRs2rHd7F+SXv/99pQv52L8nrOci9offyO+sDVie9z8\n",
       "AG8Dbqyb/x3g+kHbuN3lbFz2U5bA4jvg3/YXv09ZXfx+0/1w6U/Bhv/s4nftp35+9RDTg+frP3NO\n",
       "k/0NNT3cun80XD3MdG3+8p+Pn+1+r03HzbkY3+es189F7efy7Ydz7RzVh1oYFL/VqUExshCpD43B\n",
       "/6AjvbCvrps+d9B2q0cwPdy6ww2bdm031qF5KOGaczF+z1mvn4v6H0Z97RxXo54knQmssb2onF8F\n",
       "POe6Dm1J46fAEREdxN0wPFbSkcC3gN8EvgNsAd7ucdSZHRHRa8ZVZ7bt/ZLeRdG7OgG4KSEREdFe\n",
       "46pGERER409HvQpV0iJJ2yRtl/T+dpenKpJmSvq8pIckPSjpynL5FEmbJT0saZOkye0ua1UkTZB0\n",
       "r6RPl/M9eS4kTZb0CUnflLRV0hk9fC5Wlf9HHpB0q6QX9sq5kHSzpAFJD9QtG/K7l+dqe3k9XTjc\n",
       "/jsmKOpuxlsEzAPeLunk9paqMvuA99o+BTgTeGf53VcCm23PBe4q53vFVcBWoFYl7tVz8RHgdtsn\n",
       "A78KbKMHz4WkE4AVwKttn0rRdH0JvXMuPk5xbazX8LtLmgdcTHEdXQTcIKlpFnRMUFB3M57tfUDt\n",
       "ZryuZ3u37fvK6b0UNyBOB5YCa8vN1gIXtaeE1ZI0A1gC/BVQG8XRc+dC0kuAc2zfDEUfn+0f04Pn\n",
       "AvgJxR9Ux5SDYo6hGBDTE+fC9heAJwYtHuq7Xwiss73Pxc3NOyiur0PqpKCYDjxWN7+rXNZTyr+c\n",
       "TgO+DPTZHihXDQB9bSpW1T4EvA94rm5ZL56LWcD3JH1c0tcl3SjpRfTgubD9Q+AvgEcpAuJHtjfT\n",
       "g+eizlDf/TiK62fNsNfSTgqKnu91lzQJ+HvgKttP1q9zMSqh68+RpAuAPbbv5UBt4iC9ci4oRi2+\n",
       "GrjB9quBnzKoaaVXzoWkE4H3ACdQXAgnSfqd+m165Vw0MoLv3vS8dFJQPA7MrJufycGp2NUkvYAi\n",
       "JG6xfVu5eEDS1HL9NGBPu8pXod8Alkr6NrAOeKOkW+jNc7EL2GX7K+X8JyiCY3cPnotfB/7Z9g9s\n",
       "7wf+ATiL3jwXNUP9nxh8LZ1RLhtSJwXFV4E5kk6QNJGiM2ZDm8tUCUkCbgK22v5w3aoNwPJyejlw\n",
       "2+DPdhvbV9ueaXsWRWfl/7P9u/TmudgNPCZpbrnoPOAh4NP02Lmg6MQ/U9LR5f+X8ygGO/TiuagZ\n",
       "6v/EBuASSRMlzQLmUNzcPKSOuo9C0mIOvKviJtsfaHORKiHpbIrHvN7PgSriKop/3PXA8cBOYJnt\n",
       "H7WjjO0g6Vzg920vlTSFHjwXkn6NolN/IvAIcCnF/49ePBf/ieKC+BzwdeBy4Fh64FxIWgecC7yM\n",
       "oj/ij4BPMcR3l3Q18A6Kx1tfZfvOpvvvpKCIiIjqdVLTU0REtEGCIiIimkpQREREUwmKiIhoKkER\n",
       "ERFNJSgiIqKpBEX0HEl7y9+/IuntY7zvqwfNf3Es9x/RDgmK6EW1m4dmAb99KB8sn0zazKqDDmS/\n",
       "7lD2HzEeJSiil10LnFO+AOkqSUdI+jNJWyR9Q9K/B5C0QNIXJH0KeLBcdpukr5YvklpRLrsWOLrc\n",
       "3y3lslrtReW+H5B0v6Rldfvul/R35cuH/k+tcJKuLV/E8w1Jf1bpmYmoM67emR1RsfcDf2D7zQBl\n",
       "MPzI9umSXgj8k6RN5banAafY/tdy/lLbT0g6Gtgi6RO2V0p6p+3T6o5Rq728Ffg1ipcLvRz4iqS7\n",
       "y3WvoniJzHeBL0p6HcWziy6yfVJZthe34PtHjEhqFNHLBj+mfCHw7yTdC9wDTAFml+u21IUEwFWS\n",
       "7gO+RPEkzjnDHOts4FYX9gD/CLyWIki22P5O+Sjo+4BfAX4EPC3pJklvAZ4a9beMOEwJioiDvcv2\n",
       "aeXPibY/Vy7/aW0DSQuA3wTOtP0q4F7gqGH2a54fTLXaxjN1y54FXmD7WYq3jn0CuADYOJovEzEW\n",
       "EhTRy56keLpozZ3AFbUOa0lzJR3T4HMvBp6w/bSkkyjeY16zb4gO7y8AF5f9IC8HXk/x9N+GL18q\n",
       "31Q32fYdwH+kaLaKaIv0UUQvqv0l/w3g2bIJ6ePARynekPb18p0Ge4C3lNvXP2Z5I/AfJG0FvkXR\n",
       "/FTzMeB+SV8r35NhANuflHRWeUwD77O9R9LJPP/tYqYIsE9JOooiTN47Jt88YhTymPGIiGgqTU8R\n",
       "EdFUgiIiIppKUERERFMJioiIaCpBERERTSUoIiKiqQRFREQ0laCIiIim/j/aN/SGIk9sfgAAAABJ\n",
       "RU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f367400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(J_history, marker='o')\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('J(theta)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training accuracy: 1.0\n",
      "Test accuracy: 0.9985815602836879\n"
     ]
    }
   ],
   "source": [
    "print('Training accuracy: {acc}'.format(acc=accuracy(res.x, train_X, train_y)))\n",
    "print('Test accuracy: {acc}'.format(acc=accuracy(res.x, test_X, test_y)))"
   ]
  }
 ],
 "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}