Softmax_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "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": "markdown",
   "metadata": {},
   "source": [
    "# Data\n",
    "First get and preprocess the data. This time, we will use the complete data set and not just the samples for the numbers **`0`** and **`1`**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get data: contains 70k samples of which the last 10k are meant for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mnist = fetch_mldata('MNIST original', data_home='./data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare for concat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_all = mnist.target[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intercept term to be added"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "intercept = np.ones_like(y_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before the next step, we need to define this util function which normalizes the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "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": [
    "Now, normalize the data (zero mean and unit variance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_normalized, test_normalized = normalize_features(\n",
    "    mnist.data[:60000, :],\n",
    "    mnist.data[60000:, :],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concat **`intercept`**, **`X`**, and **`y`** so that shuffling is easier in a next step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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": [
    "Shuffle the data. As mentioned in the [Logistic_Regression](http://nbviewer.ipython.org/github/HaFl/ufldl-tutorial-python/blob/master/Logistic_Regression.ipynb) notebook already, I don't think it's needed, but let's stick with the ufldl tutorial here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.shuffle(train_all)\n",
    "np.random.shuffle(test_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, get train and test data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_X = train_all[:, :-1]\n",
    "train_y = train_all[:, -1]\n",
    "\n",
    "test_X = test_all[:, :-1]    \n",
    "test_y = test_all[:, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Softmax Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some helpful variables and initial random theta values for all classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m, n = train_X.shape\n",
    "k = np.unique(train_y).size\n",
    "theta = np.random.rand(n, k) * 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This **`indicator_mask`** will come in handy when computing the gradient later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indicator_mask = np.zeros((train_X.shape[0], theta.shape[1]), dtype=np.bool)\n",
    "for i, idx in enumerate(train_y):\n",
    "    indicator_mask[i][idx] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a helper function to keep the code DRY. It computes the probabilities of all classes for all samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def probs(theta, X, y):\n",
    "    if theta.ndim == 1:\n",
    "        theta = theta.reshape((theta.size / k, k))\n",
    "    values = np.exp(X.dot(theta)) \n",
    "    sums = np.sum(values, axis=1)\n",
    "    return (values.T / sums).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cost function of Softmax Regression. We could actually use the **`indicator_mask`** here instead of the loop at the end, but that would be computational overkill."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cost_function(theta, X, y):\n",
    "    log_probs = np.log(probs(theta, X, y))\n",
    "    cost = 0\n",
    "    for i in range(m):\n",
    "         cost -= log_probs[i][y[i]]\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient function of Softmax Regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient(theta, X, y):\n",
    "    gradient_matrix = -X.T.dot(indicator_mask - probs(theta, X, y))\n",
    "    return gradient_matrix.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright, let's run the optimization. 100 iterations are enough here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization took 195.31115102767944 seconds\n"
     ]
    }
   ],
   "source": [
    "J_history = []\n",
    "\n",
    "t0 = time.time()\n",
    "res = scipy.optimize.minimize(\n",
    "    fun=cost_function,\n",
    "    x0=theta,\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.reshape((theta.size / k, k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the evolution of **`J`** (to make sure we did the right thing)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10d5d00f0>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UVfV97/H3RxFEIRpqisiDcnVMRGwktqKJiZOaEqAp\n",
       "apZXSG8JKyE0XaTG5vamFZtbuW3TxHSpUbt0rSqJSKtXE6tgi8SHZGqSGzI1QUUJCaaiMAgkGBGi\n",
       "JhC/94/9O5nNOHM4nNn7nJk5n9das2bv7344v7OXzpff41ZEYGZmVqTDml0AMzMbepxczMyscE4u\n",
       "ZmZWOCcXMzMrnJOLmZkVzsnFzMwKV3pykbRZ0pOS1knqTLGlkram2DpJs3LnL5G0SdJGSTNy8bMk\n",
       "rU/Hrs/FR0i6K8XXSjqx7O9kZmbVNaLmEkB7REyLiLNzsWtTbFpEPAAgaQowF5gCzARukqR0zc3A\n",
       "wohoA9okzUzxhcCuFL8OuLoB38nMzKpoVLOYaoxdCNwZEfsiYjPwDDBd0jhgdER0pvNuBy5K23OA\n",
       "5Wn7HuCCwkptZmZ1aVTN5WFJj0lalItfJukJScskHZtiJwBbc+dsBcb3Eu9KcdLvLQARsR/YLWlM\n",
       "Cd/DzMxq1Ijk8q6ImAbMAj4h6d1kTVyTgTOBF4BrGlAOMzNrkGFlf0BEvJB+/0TSvcDZEfHNynFJ\n",
       "twL3p90uYGLu8glkNZautN0zXrlmErBN0jDgmIh4MV8GSV5AzcysDhHRWxfGQZWaXCQdBRweEXsk\n",
       "HQ3MAP6PpOMjYns67WJgfdpeBdwh6Vqy5q42oDMiQtLLkqYDncB84IbcNQuAtcAlwCO9laXeBzTU\n",
       "SFoaEUubXY6BwM+im59FNz+Lbv35h3nZNZexwL1pwNcw4F8i4kFJt0s6k6w/5lng4wARsUHS3cAG\n",
       "YD+wOLqXbV4M3AaMBFZHxJoUXwaskLQJ2AXMK/k7mZnZQZSaXCLiWbJ+lZ7xD1e55u+Bv+8l/j3g\n",
       "jF7ivwAu7V9JzcysSJ6h33o6ml2AAaSj2QUYQDqaXYABpKPZBRgK1AovC5MU7nMxMzs0/fnb6ZqL\n",
       "mZkVrmWSizRrjTR1drPLYWbWClomucAD74dzr3eCMTMrXwslF4BbToFJlzW7FGZmQ12LJReAUSOb\n",
       "XQIzs6GuBZPL3lebXQIzs6GuxZLLx34Mz9/Y7FKYmQ11pS9cOXDMXgPP3xjx1Opml8TMbKjzJEoz\n",
       "M+uVJ1GamdmA4uRiZmaFc3IxM7PCObmYmVnhnFzMzKxwpScXSZslPSlpnaTOFBsj6SFJP5L0oKRj\n",
       "c+cvkbRJ0kZJM3LxsyStT8euz8VHSLorxddKOrHs72RmZtU1ouYSQHtETIuIs1PsCuChiDiV7J33\n",
       "VwBImgLMBaYAM4GblN6RDNwMLIyINqBN0swUXwjsSvHrgKsb8J3MzKyKRjWL9RwnPQdYnraXAxel\n",
       "7QuBOyNiX0RsBp4BpksaB4yOiM503u25a/L3uge4oPjim5nZoWhUzeVhSY9JWpRiYyNiR9reAYxN\n",
       "2ycAW3PXbgXG9xLvSnHS7y0AEbEf2C1pTOHfwszMataI5V/eFREvSHoL8JCkjfmDERGSSl8mQDri\n",
       "b2D/62m3IyI6yv5MM7PBRFI70F7EvUpPLhHxQvr9E0n3AmcDOyQdHxHbU5PXznR6FzAxd/kEshpL\n",
       "V9ruGa9cMwnYJmkYcExEvPjGkuy7NoKXCvtiZmZDTPpHd0dlX9JV9d6r1GYxSUdJGp22jwZmAOuB\n",
       "VcCCdNoC4L60vQqYJ2m4pMlAG9AZEduBlyVNTx3884GVuWsq97qEbIBAb44u7puZmVk1ZddcxgL3\n",
       "pgFfw4B/iYgHJT0G3C1pIbAZuBQgIjZIuhvYAOwHFkf3ypqLgduAkcDqiFiT4suAFZI2AbuAeX2U\n",
       "xcnFzKxBWmZVZIhpETze7LKYmQ0WXhW5Nq65mJk1iJOLmZkVzsnFzMwK5+RiZmaFc3IxM7PCObmY\n",
       "mVnhnFzMzKxwTi5mZla4VkouRzW7AGZmraKVkotrLmZmDeLkYmZmhXNyMTOzwjm5mJlZ4ZxczMys\n",
       "cE4uZmZWOCcXMzMrXOnJRdLhktZJuj/tL5W0NcXWSZqVO3eJpE2SNkqakYufJWl9OnZ9Lj5C0l0p\n",
       "vlbSiVWK4uRiZtYgjai5XE722uLKKy8DuDYipqWfBwAkTQHmAlOAmcBNSu9HBm4GFkZEG9AmaWaK\n",
       "LwR2pfh1wNVVyuFJlGZmDVJqcpE0AZgN3ApUEoVy23kXAndGxL6I2Aw8A0yXNA4YHRGd6bzbgYvS\n",
       "9hxgedq+B7igSnGOknr9XDMzK1jZNZfrgE8Dr+diAVwm6QlJyyQdm+InAFtz520FxvcS70px0u8t\n",
       "ABGxH9gtaUwfZfkFMLIf38XMzGpUWnKR9AFgZ0Ss48Cays3AZOBM4AXgmrLKcKC/CjhuaerzaW/M\n",
       "Z5qZDR6S2tPfyKWSlvbnXsMKKlNv3gnMkTQbOBJ4k6TbI+LDlRMk3Qrcn3a7gIm56yeQ1Vi60nbP\n",
       "eOWaScA2ScOAYyLixd6L89mfwGdvimBzP7+XmdmQFBEdQEdlX9JV9d6rtJpLRFwZERMjYjIwD/h6\n",
       "RHw49aFUXAysT9urgHmShkuaDLQBnRGxHXhZ0vTUwT8fWJm7ZkHavgR4pEqRXsEjxszMGqLMmkue\n",
       "6B4t9gVJb0/7zwIfB4iIDZLuJhtZth9YHBGVaxYDt5H1mayOiDUpvgxYIWkTsIssifXl5zi5mJk1\n",
       "hLr/fg9dkgLiUWBpBN9odnnMzAYDSRERdY2ybaUZ+j/Hc13MzBqi1ZKLm8XMzBrAycXMzArn5GJm\n",
       "ZoVzcjEzs8K1UnLxPBczswZppeTimouZWYM4uZiZWeGcXMzMrHCtllw8idLMrAFaLbm45mJm1gBO\n",
       "LmZmVjgnFzMzK5yTi5mZFa6VkosnUZqZNUgrJRfXXMzMGqT05CLpcEnrJN2f9sdIekjSjyQ9KOnY\n",
       "3LlLJG2StFHSjFz8LEnr07Hrc/ERku5K8bWSTqxSFCcXM7MGaUTN5XKyVxdXXnl5BfBQRJxK9s77\n",
       "KwAkTQHmAlOAmcBNkipvQLsZWBgRbUCbpJkpvhDYleLXAVdXKccvgMMljijsm5mZWa9KTS6SJgCz\n",
       "gVuBSqKYAyxP28uBi9L2hcCdEbEvIjYDzwDTJY0DRkdEZzrv9tw1+XvdA1zQV1kiCDyR0sysIcqu\n",
       "uVwHfBp4PRcbGxE70vYOYGzaPgHYmjtvKzC+l3hXipN+bwGIiP3AbkljqpTHTWNmZg0wrKwbS/oA\n",
       "sDMi1klq7+2ciAhJ0duxEsqzFD49Ar5yhbT5XyOioxGfa2Y2WKS/1e1F3Ku05AK8E5gjaTZwJPAm\n",
       "SSuAHZKOj4jtqclrZzq/C5iYu34CWY2lK233jFeumQRskzQMOCYiXuytMBGxVOJi+MKXIni8oO9o\n",
       "ZjZkpH90d1T2JV1V771KaxaLiCsjYmJETAbmAV+PiPnAKmBBOm0BcF/aXgXMkzRc0mSgDeiMiO3A\n",
       "y5Kmpw7++cDK3DWVe11CNkCgGjeLmZk1QJk1l54qzV+fB+6WtBDYDFwKEBEbJN1NNrJsP7A4IirX\n",
       "LAZuA0YCqyNiTYovA1ZI2gTsIkti1XgipZlZA6j77/fQJSl177ASuC2Ce5tdJjOzga7yt7Oea1tp\n",
       "hj64WczMrCGcXMzMrHCtmFw8idLMrGStmFxcczEzK5mTi5mZFc7JxczMCtfIeS5NJU2dDb+zCN78\n",
       "G9IPToUtN0Q8tbrZ5TIzG4paJrnAudfDLaekneNh0cnSVJxgzMyK1zKTKLsXCMibvSZi9ayGF8jM\n",
       "bBDwJMq6jRrZ7BKYmQ1FLZ5c9r7a7BKYmQ1FLZRcFj1z4P7HfgzP39icspiZDW0t1Ody+u/DSX8G\n",
       "Z70P/vNr8PyN7sw3M+tbf/pcWia5VB6QxF5gXAR7mlwsM7MBzR36h2YXcFyzC2FmNpTVlFwkHS3p\n",
       "bZLeKqmmGe6SjpT0XUmPS9og6XMpvlTSVknr0s+s3DVLJG2StFHSjFz8LEnr07Hrc/ERku5K8bWS\n",
       "TqyhaD8FfqOW72BmZvXpcxKlpNHAIrK3Ox4H7AAEjJW0C/gX4JaI2Nvb9RHxmqT3RsQr6f3235J0\n",
       "HtmEk2sj4toenzcFmAtMAcYDD0tqS2+jvBlYGBGdklZLmpneRrkQ2BURbZLmAldz8LdR/hTXXMzM\n",
       "SlWt5nIfsAf4g4j4bxFxbkScExGTgQ+QrdO1ssr1RMQraXM4cDjws7TfWxvehcCdEbEvIjYDzwDT\n",
       "JY0DRkdEZzrvduCitD0HWJ627wEuqFaexM1iZmYl6zO5RMQFEXFLROzo5dj2iPiniKj6x1zSYZIe\n",
       "J6v1fCMink6HLpP0hKRlko5NsROArbnLt5LVYHrGu1Kc9HtLKtN+YLekMdXKhJvFzMxKV2ufyxhJ\n",
       "0yW9p/JTy3UR8XpEnAlMAN4jqZ2siWsycCbwAnBNfUWvm2suZmYlO+jClZIWAZ8kSxCPA+cA3wF+\n",
       "t9YPiYjdkv4d+O2I6Mjd+1bg/rTbBUzMXTaBrMbSlbZ7xivXTAK2pX6dYyLixT6+x9Js6/d/B+bu\n",
       "h/m1Ft/MrCWkCkB7EfeqpeZyOXA28FxEvBeYBuw+2EWSjqs0eUkaCfwesE7S8bnTLgbWp+1VwDxJ\n",
       "wyVNBtqAzojYDrycak4iyworc9csSNuXAI/0VZ6IWBoRS+Hfbof5v6jhe5uZtZSI6Kj8rcz+Xtav\n",
       "liX3X4uIVyUh6ciI2CjprTVcNw5YLukwsiS2IiIekXS7pDPJRo09C3wcICI2SLob2ADsBxZH9wzP\n",
       "xcBtwEhgdRopBrAMWCFpE1lz18FGioGbxczMSnfQGfqS7gM+QlaDuYBsxNewiJhdfvGK0WOG/pnA\n",
       "8gje3uRimZkNaA1b/iW1x70JWBMRv6znA5uhR3KZCKyN+PWIMzMz60Wpy79IWlHZTu1xq8iaowar\n",
       "XcBxUq9zbczMrAC1dOhPze+kUVlnlVOc8kXwClmfTk3L2JiZ2aHrM7lIulLSHuAMSXsqP8BOslFa\n",
       "g5k79c3MSlRLh/7nI+KKBpWnFD3bDSW+D/xxBI81sVhmZgNa2UvuXylpvqS/Th82SdLZ9XzYAOIl\n",
       "YMzMSlRLcrkJOBf4w7S/N8UGMzeLmZmVqJZJlNMjYpqkdQAR8aKkI0ouV9lcczEzK1EtNZdfSjq8\n",
       "siPpLcDr5RWpIVxzMTMrUS3J5UbgXuA3Jf098G3gc6WWqnx+YZiZWYkO2iwWEf8s6Xt0v4jrwoj4\n",
       "QbnFKp2bxczMSlRLnwvAj4CX0/khaVJEPF9esUrnZjEzsxLV8j6Xy4CryCZP/ip36IyyCtUArrmY\n",
       "mZWolprLnwFvjYhdZRemgVxzMTMrUS0d+s+TNYkNJT/Fi1eamZWmz5qLpD9Pm/8FdEj6N6CyzH5E\n",
       "xLVlF64sEbwi8TpwFPDzZpfHzGyoqVZzGQ2MIqu5PAQMT/uj0rGqJB0p6buSHpe0QdLnUnyMpIck\n",
       "/UjSg5VXIadjSyRtkrRR0oxc/CxJ69Ox63PxEZLuSvG1kk48hO/upjEzs5L0WXOpvD9Z0qURcXf+\n",
       "mKRLD3bjiHhN0nsj4pW0TP+3JJ0HzAEeiogvSPpL4ArgCklTgLnAFGA88LCktvSq45uBhRHRKWm1\n",
       "pJnpVccLgV0R0SZpLnA1tb3qGLrnujxX4/lmZlajWvpcltQYe4OIeCVtDgcOJ3tF8hxgeYovBy5K\n",
       "2xcCd0bEvojYDDwDTJc0DhgdEZ3pvNtz1+TvdQ/dc3Fq4RFjZmYlqdbnMguYDYyXdAP8uvN7NLCv\n",
       "lptLOgz4PnAycHNEPC1pbETsSKfsAMam7ROAtbnLt5LVYPal7YquFCf93gIQEfsl7ZY0JiJerKF4\n",
       "bhYzMytJtaHI24DvkdUovkeWXALYA3yqlptHxOvAmZKOAb4m6b09joek6i+UKY9rLmZmJanW5/IE\n",
       "8ISkOyLil32dV4uI2C3p38lej7xD0vERsT01ee1Mp3UBE3OXTSCrsXSl7Z7xyjWTgG2pX+eYvmot\n",
       "kpbmdjsgXHMxM8uR1A60F3Gvaq85/ndJ/51eEpCkoyXNlbS6yvXHVUaCSRoJ/B6wjuwVyQvSaQuA\n",
       "+9L2KmCepOGSJgNtQGdEbAdeljRdkoD5wMrcNZV7XQI80ld5ImJp7qcDL15pZnaAiOjI/63sz72q\n",
       "NYt9BPhT4P9I+hXwAlnT2PHpurvo/sPem3HA8tTvchiwIiIeSe+FuVvSQmAzcGn6Uhsk3Q1sAPYD\n",
       "i6P7HcyLgduAkcDqNFIMYBmwQtImsj6UWkeKgZvFzMxKo+6/31VOksYClTkkz+U65AeF3t4DLTED\n",
       "+IsI3tekYpmZDWi9/e2sVbXRYnvJOvB7O/YLsqHCn4mIh+v54AHANRczs5JU69Af1dex1Hl+OnBH\n",
       "+j0I/Y/fgolvk/6rA/a8BltuiHiqzz4kMzOrXa3vczlAROwnG0l2Y8HlaQhp6mx45/+Gzx8JnJ9F\n",
       "F50sTcUJxsys/2rqcxnserYbSrPWwAPvf+OZs9dErJ7VyLKZmQ1U/elzqWX5lyFo9JG9x0eNbGw5\n",
       "zMyGphZNLnte6z2+99XGlsPMbGhq0eSy5QZY9MyBsY/9GJ4flH1IZmYDTUv2uWSxqbPh7C/Asb8B\n",
       "Gx+H5290Z76ZWbf+9Lm0bHLJ4rQDfxfBeY0vlZnZwOYO/fr9EHhrswthZjbUtHpy2Q4MlzxT38ys\n",
       "SC2dXCIIXHsxMytcSyeXxMnFzKxgTi5OLmZmhXNygY04uZiZFcrJJau5vK3ZhTAzG0pKTS6SJkr6\n",
       "hqSnJT0l6ZMpvlTSVknr0s+s3DVLJG2StFHSjFz8LEnr07Hrc/ERku5K8bWSTuTQPANMlupbIdrM\n",
       "zN6o7JrLPuBTEXE6cA7wCUmnkb2E7NqImJZ+HgCQNAWYC0wBZgI3SapM4LkZWBgRbUCbpJkpvhDY\n",
       "leLXAVcfSgEjeJXsFc6T+/NFzcysW6nJJSK2R8TjaXsv8ANgfDrc26zPC4E7I2JfRGwmq1VMlzQO\n",
       "GB0Rnem824GL0vYcYHnavge4oI6iut/FzKxADetzkXQSMA1Ym0KXSXpC0jJJx6bYCcDW3GVbyZJR\n",
       "z3gX3UlqPLAFfv0Ss92Sxhxi8TxizMysQA3pZ5A0CvgqcHlE7JV0M/A36fDfAteQNW+VWYalud2O\n",
       "iOjI7f8QeEeZn29mNtBJagfai7hX6clF0hFkzVX/HBH3AUTEztzxW4H7024XMDF3+QSyGktX2u4Z\n",
       "r1wzCdgmaRhwTES82LMcEbG0SjF/CHyo9m9lZjb0pH90d1T2JV1V773KHi0mYBmwISK+mIuPy512\n",
       "MbA+ba8C5kkaLmky0AZ0RsR24GVJ09M95wMrc9csSNuXAI/UUVT3uZiZFajsmsu7gD8CnpS0LsWu\n",
       "BD4k6UyyUWPPAh8HiIgNku4GNgD7gcXR/U6AxcBtwEhgdUSsSfFlwApJm4BdwLxDL+bUM2HOGOnZ\n",
       "b8FLe2HLDX63i5lZ/Vr6fS7Zsamz4dzr4ZZTuqOLnoHvXO4EY2atzO9z6ZeJnzwwsUC2P+my5pTH\n",
       "zGzwc3Jh9JG9x0eNbGw5zMyGDicX9rzWe3zvq40th5nZ0OHkwpYbsj6WvI/9GJ6/sTnlMTMb/Fq+\n",
       "Qz87PnV21scy9Wx44TlY9xl35ptZq+tPh76TywHnsQT4zQg+1YBimZkNaB4tVpz/AM5vdiHMzAY7\n",
       "11wOOI/hZBMxJ0Swu/ySmZkNXK65FCSCXwKdwHnNLouZ2WDm5PJGj+KmMTOzfnFyeaP/AN7T7EKY\n",
       "mQ1m7nN5w7mMBH4CHB/B3nJLZmY2cPWnz6UhLwsbXKa+F2b8EnZ/W9r2gldINjM7dE4uOd0rJF/7\n",
       "ZuDNwG/BopOlqTjBmJnVzn0uB/AKyWZmRXByOYBXSDYzK0LZrzmeKOkbkp6W9JSkT6b4GEkPSfqR\n",
       "pAclHZu7ZomkTZI2SpqRi58laX06dn0uPkLSXSm+VtKJ9ZfYKySbmRWh7JrLPuBTEXE6cA7wCUmn\n",
       "AVcAD0XEqWTvvL8CQNIUYC4wBZgJ3CSpMlLhZmBhRLQBbZJmpvhCYFeKXwdcXX9xe1sh+fJdXiHZ\n",
       "zOzQNHQosqT7gH9MP+dHxA5JxwMdEfE2SUuA1yPi6nT+GmAp8Bzw9Yg4LcXnAe0R8SfpnKsi4ruS\n",
       "hgEvRMRbenzuIQxFrqyQPGokvH4YvPsd8PVOGHFYVrPx6DEzaw2DYiiypJOAacB3gbERsSMd2gGM\n",
       "TdsnAGtzl20FxpPVgLbm4l0pTvq9BSAi9kvaLWlMRLxYTzlT4lidlXnqbJh0B6x8b/cZHj1mZnYw\n",
       "DUkukkYB9wCXR8Se7pYuiIiQVHr1SdLS3G5HRHQc/KqJn4RrjzkwdsspMPsyUgIyMxsqJLUD7UXc\n",
       "q/TkIukIssSyIiLuS+Edko6PiO2SxgE7U7wLmJi7fAJZjaUrbfeMV66ZBGxLzWLH9FZriYilh156\n",
       "jx4zs9aR/tHdUdmXdFW99yp7tJiAZcCGiPhi7tAqYEHaXgDcl4vPkzRc0mSgDeiMiO3Ay5Kmp3vO\n",
       "B1b2cq9LyAYIFMSjx8zM6lFqh76k88hWGX4SqHzQErJl7e8mq3FsBi6NiJfSNVcCHwX2kzWjfS3F\n",
       "zwJuA0YCqyOiMqx5BLCCrD9nFzAvIjb3KEddnVLdM/YrEysfBa57DbQJXt3mzn0zG8r8muOD6NcD\n",
       "+vXosf3j4cS3wi3Du48uega+c7kTjJkNRU4uB9GfB9R9j1lr4IH3v/HI7DURq2f1595mZgOR30TZ\n",
       "EO7cNzOrlZNLzXp27j8KfAb41RnSrDVZ85mZmYGX3D8EW26ARSdnnfuPAl8DPgswBni/J1eamXVz\n",
       "n8sh3afSuX/YdPi3Nx949FHgcz+F0U97mRgzGwoGxfIvQ0FlaRjp0g7g/O4jlZrMA8d1x12TMbPW\n",
       "5T6XuvTsf3mQ1ESW45eMmVnrcnKpS8+l+XtWACud/Uef685+M2tF7nOp+575pfn3nJ6axOjR2Z/2\n",
       "r3kFDv+xZ/Wb2WDiSZQHUUZyOfD++WViPgP8XTrSM9GAZ/Wb2WDh5HIQZSeX7DMqNZmjz4WvpGX6\n",
       "84mmwqPKzGxwcHI5iEYkl+7Pyi8TszT9VLjJzMwGDy//MqDkO/v39ziWH1VWSTQrj4J/PQOWvB9O\n",
       "/or0wSc9CMDMBjvPcylYxFOrpalkb6vcPx4+ejJ86ajsaP5x95VoOCP78TwZMxu8nFxKUJlsCZW+\n",
       "mNm5UWWkUWV9JZqK+afA1uXSpU/D9lFwpODYPe6nMbPBwMmlZG9MNIvSqLJ8k1lv82QqM/4fPb+X\n",
       "fpp3Sx90P42ZDVhlv+b4S5J2SFqfiy2VtFXSuvQzK3dsiaRNkjZKmpGLnyVpfTp2fS4+QtJdKb5W\n",
       "0ollfp/+ypLAdy6H2Wugcz189JXsSLW+mWr9NA+8H069RZrxPenSDvfVmNlAUXbN5cvAjcDtuVgA\n",
       "10bEtfkTJU0B5gJTgPHAw5LaIhvOdjOwMCI6Ja2WNDMi1gALgV0R0SZpLnA1MK/k79QvvTeZVeub\n",
       "qdZ89ihw2gnw2RO69yu1mp/+0k1pZtYspSaXiPimpJN6OdTb0LYLgTsjYh+wWdIzwHRJzwGjI6Iz\n",
       "nXc7cBGwBpgDXJXi9wD/WGDxS1db30y15rO+ajWPnvHGyZsfPEOasd3JxswaoVl9LpdJ+jDwGPDn\n",
       "EfEScAKwNnfOVrIazL60XdGV4qTfWwAiYr+k3ZLGRMSLZX+BovXdNzMD+CuyRNGz+ayvWs2h1HBe\n",
       "3QZbvgMTz83etunEY2b914zkcjPwN2n7b4FryJq3SiVpaW63IyI6yv7Meh04nHnUSNh+NMwQMPzA\n",
       "5rO+ajW11nA4I6vlrPhduOWI7uNuWjNrRZLagfYi7tXw5BIROyvbkm4F7k+7XcDE3KkTyGosXWm7\n",
       "Z7xyzSRgm6RhwDF91VoiYmkR5W+UfE0m78Dms+1Hw0fHw5fGHZhoaq3hVPbziaXWpjUPjzYbatI/\n",
       "ujsq+5Ku6vPkg2h4cpE0LiJeSLsXA5WRZKuAOyRdS9bc1QZ0RkRIelnSdKATmA/ckLtmAVlz2iXA\n",
       "Iw36Gk3TM+n0Pigg35QG1ftt6mla63UZGzezmdmvlZpcJN1J9mbG4yRtIet8b5d0JtmosWeBjwNE\n",
       "xAZJdwMbyP4aLo7uhc8WA7cBI4HVaaQYwDJghaRNwC4G+EixMvQ9KKDSlHbs3r5rOFBf01phzWxO\n",
       "QmZDlBeubBHdqzbvHw8Tcv02jwIr9mXJIL+Kc88VnZfSvQhnfrvneX29cqC32k7lcyv7lUU8qyUh\n",
       "N8eZNUp//nZ6hn6L6LuGs/dVeH4tzD6n9qa1IprZau3reZTumtChNMflk5ATklmjObm0oL4GC0Ct\n",
       "TWvVhkfX2sxWTxKqtTmuWo0JDjIwodYE5SY9syqcXOwAh554eg6PngEsSs1d1Uaw1ZOEah311p+B\n",
       "CbUkqPx5lf1a+pVqTVauddng5+RiNatteHStzWz1JKFam+PqGZhQ2a8lQfW3Sa/epAZ917qKqIHV\n",
       "c54TnvXOycX6rb5mtnqSUK3NcbXO+ak3QfW3Sa/oWlcRNbB6zqsYaAmv9POcTGvg5GKlqpZ48mpL\n",
       "QrU2x+WTUK0DE6odK7pJr+haV9lJra/zYOAlvLLPg/4n00GRQNN59XNysQGh1iSU13dzXG+DEQ42\n",
       "MAFqT1D9bdIrutZVdlIbTAlvoNceqx0biOfVP4PDycUGrXoSEtRTS+qZoPrbpFfveUUnq/6eBwMv\n",
       "4Q302mO1YwPxvPo5uVjLqTcp9aXYZFVPravspNbXeTDwEt5Arz0WcY9Gnlc/Jxezfio6WVVTfR5S\n",
       "WUntYOcNpIQ30GuPRdyjkefVz8u/mFm/dC8tVEk8R+YSz5E9EtSkcwb3eduPhlNSMq11eaNajw3E\n",
       "80S9fzudXMzMDkH/k+lAT6D58x6Y6eRShZOLmdmh68/fzsOKLoyZmZmTi5mZFa7U5CLpS5J2SFqf\n",
       "i42R9JCkH0l6UNKxuWNLJG2StFHSjFz8LEnr07Hrc/ERku5K8bWSTizz+5iZWW3Krrl8GZjZI3YF\n",
       "8FBEnEpwB7ZuAAAGxUlEQVT2WuIrACRNAeYCU9I1N0mqtPXdDCyMiDagTVLlnguBXSl+HXB1mV9m\n",
       "KJDU3uwyDBR+Ft38LLr5WRSj1OQSEd8EftYjPAdYnraXAxel7QuBOyNiX0RsBp4BpksaB4yOiM50\n",
       "3u25a/L3uge4oPAvMfS0N7sAA0h7swswgLQ3uwADSHuzCzAUNKPPZWxE7EjbO4CxafsEYGvuvK3A\n",
       "+F7iXSlO+r0FICL2A7sljSmp3GZmVqOmduhHNg566I+FNjNrMc1Y/mWHpOMjYntq8tqZ4l3AxNx5\n",
       "E8hqLF1pu2e8cs0kYJukYcAxEfFibx8qyUkskXRVs8swUPhZdPOz6OZn0X/NSC6rgAVkne8LgPty\n",
       "8TskXUvW3NUGdEZESHpZ0nSgE5gP3NDjXmuBS8gGCLyBJ1CamTVWqTP0Jd0JnA8cR9a/8tfASuBu\n",
       "shrHZuDSiHgpnX8l8FGyldMuj4ivpfhZwG3ASGB1RHwyxUcAK4BpwC5gXhoMYGZmTdQSy7+YmVlj\n",
       "DekZ+pJmpgmZmyT9ZbPL00iSJkr6hqSnJT0lqVLb63MS61An6XBJ6yTdn/Zb8llIOlbSVyX9QNIG\n",
       "SdNb+FksSf+PrJd0R5qY3RLPoqhJ7n0ZsslF0uHAP5JNyJwCfEjSac0tVUPtAz4VEacD5wCfSN+/\n",
       "10msLeJyYAPdIxRb9VlcT9a8fBrwW8BGWvBZSDoJWAS8IyLOAA4H5tE6z6K/k9yr5o8hm1yAs4Fn\n",
       "ImJzROwD/i/ZRM2WEBHbI+LxtL0X+AHZQIm+JrEOaZImALOBW+l+MXjLPQtJxwDvjogvQTY/LCJ2\n",
       "04LPAniZ7B9hR6XRpkcB22iRZ1HAJPezq91/KCeXX0+wTCqTMltO+hfaNOC79D2Jdai7Dvg08Hou\n",
       "1orPYjLwE0lflvR9SbdIOpoWfBZp2sI1wPNkSeWliHiIFnwWOYc6yb1PQzm5eKQCIGkU2dI4l0fE\n",
       "nvyxVpnEKukDwM6IWEd3reUArfIsyKYfvAO4KSLeAfycHs0+rfIsJJ0M/BlwEtkfz1GS/ih/Tqs8\n",
       "i97U8N2rPpehnFx6TsqcyIGZd8iTdARZYlkREZX5RDskHZ+O5yexDmXvBOZIeha4E/hdSStozWex\n",
       "FdgaEf+Z9r9Klmy2t+Cz+G3g/0XErrR81L8C59Kaz6Kir/8nepvk3lXtRkM5uTxGtoLySZKGk3VG\n",
       "rWpymRomrSi9DNgQEV/MHapMPIUDJ7EOWRFxZURMjIjJZB22X4+I+bTms9gObJF0agq9D3gauJ8W\n",
       "exZkAxnOkTQy/f/yPrIBH634LCr6+n9iFTBP0nBJk0mT3KvdaEjPc5E0C/gi2SiQZRHxuSYXqWEk\n",
       "nQc8CjxJd/V1Cdl/EL1OYm0Fks4H/jwi5qRFTlvuWUh6O9nAhuHAj4GPkP0/0orP4i/I/oi+Dnwf\n",
       "+BgwmhZ4FkVNcu/z/kM5uZiZWXMM5WYxMzNrEicXMzMrnJOLmZkVzsnFzMwK5+RiZmaFc3IxM7PC\n",
       "ObmYHYSkven3iZI+VPC9r+yx/+0i72/WLE4uZgdXmQw2GfjDQ7kwrbZbzZIDPijiXYdyf7OBysnF\n",
       "rHafB96dXjh2uaTDJP2DpE5JT0j6YwBJ7ZK+KWkl8FSK3SfpsfTitkUp9nlgZLrfihSr1JKU7r1e\n",
       "0pOSLs3du0PSV9LLvv65UjhJn08vvnpC0j809MmY9XCwf1WZWbe/BP5XRPwBQEomL0XE2ZJGAN+S\n",
       "9GA6dxpwekQ8l/Y/EhE/kzQS6JT01Yi4QtInImJa7jMqtaQPAm8ne5nXW4D/lPRoOnYm2UubXgC+\n",
       "LeldZOtkXRQRb0tle1MJ39+sZq65mNWu53L9M4APS1oHrAXGAKekY525xAJwuaTHge+QrS7bdpDP\n",
       "Og+4IzI7gf8Afocs+XRGxLa0JPrjwInAS8BrkpZJuhh4te5vaVYAJxez/vnTiJiWfk6OiIdT/OeV\n",
       "EyS1AxcA50TEmcA64MiD3Dd4YzKr1Gp+kYv9CjgiIn5F9mbArwIfANbU82XMiuLkYla7PWQr5lZ8\n",
       "DVhc6bSXdKqko3q57k3AzyLiNUlvA87JHdvXR6f/N4G5qV/nLcB7yFa07vVlZ+ltksdGxAPA/yRr\n",
       "UjNrGve5mB1cpcbwBPCr1Lz1ZeAGsrcYfj+9D2QncHE6P7/c+BrgTyRtAH5I1jRW8U/Ak5K+l94x\n",
       "EwARca+kc9NnBvDpiNgp6TTe+AbAIEt6KyUdSZaAPlXINzerk5fcNzOzwrlZzMzMCufkYmZmhXNy\n",
       "MTOzwjm5mJlZ4ZxczMyscE4uZmZWOCcXMzMrnJOLmZkV7v8DX2qxSUZpL/8AAAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d5b9f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(J_history, marker='o')\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('J(theta)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even the **`optimal_theta`** still results in a relatively high cost value **`J`**. That's not necessarily a bad sign, but rather the consequence of the model being more complex than e.g., Logistic Regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def accuracy(theta, X, y):\n",
    "    correct = np.sum(np.argmax(probs(theta, X, y), axis=1) == y)\n",
    "    return correct / y.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training accuracy: 0.9435833333333333\n",
      "Test accuracy: 0.9235\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
}