02_GD_logistic.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "plt.xkcd()\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# learning rate hyperparameter (step size scaling)\n",
    "ALPHA = 1e-4\n",
    "\n",
    "# weight decay hyperparameter (regularization)\n",
    "LAMBD = 1e-3\n",
    "\n",
    "# threshold for convergence test\n",
    "DIFF_THRESHOLD = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1. / (1. + np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def threshold(a):\n",
    "    ''' Threshold given value to return 1 if greater than or equal to 0.5,\n",
    "        or 0 if otherwise.\n",
    "    '''\n",
    "    if a >= 0.5:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "h = np.vectorize(threshold, otypes=[np.dtype('int8')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cost(theta, X, y):\n",
    "    theta = np.matrix(theta)\n",
    "    \n",
    "    m = X.shape[0]\n",
    "    \n",
    "    z = np.dot(X, theta.T).sum(axis=1)\n",
    "    a = sigmoid(z)\n",
    "    y_hat = h(a)\n",
    "    \n",
    "    cost_simple = np.add(np.multiply(y, np.log1p(y_hat)),\n",
    "                         np.multiply((1-y), np.log1p(1-y_hat)))\n",
    "    \n",
    "    reg = (LAMBD/2) * np.square(theta[:, 1:X.shape[1]]).sum()\n",
    "    \n",
    "    cost = (-1./m) * (cost_simple.sum() + reg)\n",
    "    \n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient(theta, X, y):\n",
    "    theta = np.matrix(theta)\n",
    "    \n",
    "    m = X.shape[0]\n",
    "    \n",
    "    z = np.dot(X, theta.T).sum(axis=1)\n",
    "    a = sigmoid(z)\n",
    "    y_hat = h(a)\n",
    "    err = y_hat - y\n",
    "    \n",
    "    derivs = np.multiply(err, X).sum(axis=0).A1 / m\n",
    "    derivs[1:X.shape[1]] += (LAMBD/m) * theta[:, 1:X.shape[1]].sum()\n",
    "    \n",
    "    return derivs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient_descent(theta, X, y):\n",
    "    iters = 1\n",
    "    \n",
    "    m = X.shape[0]\n",
    "    \n",
    "    cost_history = []\n",
    "    \n",
    "    theta_history = [ np.ones(X.shape[1]) ]\n",
    "    \n",
    "    while True:\n",
    "        theta -= ALPHA * gradient(theta, X, y)\n",
    "        \n",
    "        # convergence test\n",
    "        last_theta = theta_history[-1]\n",
    "        diffs = np.sum(np.abs(np.subtract(last_theta, theta)) < DIFF_THRESHOLD)\n",
    "        \n",
    "        if diffs == X.shape[1]:\n",
    "            print('... convergence at {} iterations'.format(iters))\n",
    "            break\n",
    "        else:\n",
    "            theta_history.append(np.copy(theta))\n",
    "            cost_history.append(cost(theta, X, y))\n",
    "            iters += 1\n",
    "\n",
    "    return theta, cost_history, theta_history                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "iris = load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(888)\n",
    "\n",
    "idx_setosa     = np.where(iris.target==0)[0]\n",
    "idx_versicolor = np.where(iris.target==1)[0]\n",
    "idx_virginica  = np.where(iris.target==2)[0]\n",
    "\n",
    "train_setosa     = np.random.choice(idx_setosa, size=40, replace=False)\n",
    "train_versicolor = np.random.choice(idx_versicolor, size=40, replace=False)\n",
    "train_virginica  = np.random.choice(idx_virginica, size=40, replace=False)\n",
    "training_classes = [train_setosa, train_versicolor, train_virginica]\n",
    "\n",
    "idx_train = np.concatenate(training_classes)\n",
    "np.random.shuffle(idx_train)\n",
    "\n",
    "idx_test = np.array([i for i in range(iris.data.shape[0]) if i not in idx_train])\n",
    "np.random.shuffle(idx_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = iris.data[idx_train, :]\n",
    "X = np.insert(X, 0, np.ones(idx_train.size), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Training classifier for Iris-setosa\n",
      "... convergence at 4 iterations\n",
      "\n",
      "Training classifier for Iris-versicolor\n",
      "... convergence at 1625 iterations\n",
      "\n",
      "Training classifier for Iris-virginica\n",
      "... convergence at 139 iterations\n"
     ]
    }
   ],
   "source": [
    "classes = ['Iris-setosa', 'Iris-versicolor', 'Iris-virginica']\n",
    "classifiers = []\n",
    "\n",
    "for idx,cls in enumerate(classes):\n",
    "    print('\\nTraining classifier for {}'.format(cls))\n",
    "    \n",
    "    iters = 1\n",
    "    \n",
    "    y = np.array([e in training_classes[idx] for e in idx_train], dtype=np.int8)\n",
    "    y = y.reshape((idx_train.size,1))\n",
    "    \n",
    "    theta = np.random.normal(scale=1e-3, size=X.shape[1])\n",
    "    \n",
    "    final_theta, cost_history, theta_history = gradient_descent(theta, X, y)\n",
    "    \n",
    "    classifiers.append(final_theta)\n",
    "\n",
    "classifiers = np.vstack(classifiers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_test = iris.data[idx_test, :]\n",
    "X_test = np.insert(X_test, 0, np.ones(idx_test.size), axis=1)\n",
    "\n",
    "Y = iris.target[idx_test]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> test  1: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  2: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  3: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  4: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  5: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  6: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  7: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test  8: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  9: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 10: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 11: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 12: classifier predicted Iris-virginica, actual class is Iris-versicolor\n",
      "> test 13: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 14: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 15: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 16: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 17: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 18: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 19: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 20: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 21: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 22: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 23: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 24: classifier predicted Iris-versicolor, actual class is Iris-virginica\n",
      "> test 25: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 26: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 27: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 28: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 29: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 30: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "\n",
      "... Overall classification accuracy with gradient descent: 0.933\n"
     ]
    }
   ],
   "source": [
    "correct = 0\n",
    "for i,x in enumerate(X_test):\n",
    "    z = np.dot(classifiers, x.T)\n",
    "    g = sigmoid(z)\n",
    "    y_hat = g.argmax()\n",
    "    y = iris.target[idx_test[i]]\n",
    "    \n",
    "    if y == y_hat:\n",
    "        correct += 1\n",
    "\n",
    "    print('> test {0:2d}: classifier predicted {1}, actual class is {2}'.format(i+1,\n",
    "                                                                                classes[y_hat],\n",
    "                                                                                classes[y]))\n",
    "\n",
    "print('\\n... Overall classification accuracy with gradient descent: {:.3f}'.format(correct/X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Training classifier for Iris-setosa\n",
      "success?: True\n",
      "  Local minimum reached (|pg| ~= 0)\n",
      "\n",
      "Training classifier for Iris-versicolor\n",
      "success?: True\n",
      "  Converged (|f_n-f_(n-1)| ~= 0)\n",
      "\n",
      "Training classifier for Iris-virginica\n",
      "success?: True\n",
      "  Converged (|f_n-f_(n-1)| ~= 0)\n",
      "> test  1: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  2: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  3: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  4: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  5: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  6: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  7: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test  8: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  9: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 10: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 11: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 12: classifier predicted Iris-virginica, actual class is Iris-versicolor\n",
      "> test 13: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 14: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 15: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 16: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 17: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 18: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 19: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 20: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 21: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 22: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 23: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 24: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 25: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 26: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 27: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 28: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 29: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 30: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "\n",
      "... Overall classification accuracy with scipy.optimize.minimize: 0.967\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "classes = ['Iris-setosa', 'Iris-versicolor', 'Iris-virginica']\n",
    "classifiers2 = []\n",
    "\n",
    "for idx,cls in enumerate(classes):\n",
    "    print('\\nTraining classifier for {}'.format(cls))\n",
    "    \n",
    "    y = np.array([e in training_classes[idx] for e in idx_train], dtype=np.int8)\n",
    "    y = y.reshape((idx_train.size,1))\n",
    "\n",
    "    res = minimize(cost,\n",
    "                   np.random.normal(scale=1e-3, size=X.shape[1]),\n",
    "                   args=(X,y),\n",
    "                   method='TNC',\n",
    "                   jac=gradient,\n",
    "                   options={'maxiter': 1000, 'disp': True})\n",
    "    \n",
    "    print(\"success?: {}\\n  {}\".format(res.success, res.message))\n",
    "    classifiers2.append(res.x)\n",
    "classifiers2 = np.vstack(classifiers2)\n",
    "\n",
    "correct = 0\n",
    "for i,x in enumerate(X_test):\n",
    "    z = np.dot(classifiers2, x.T)\n",
    "    g = sigmoid(z)\n",
    "    y_hat = g.argmax()\n",
    "    y = iris.target[idx_test[i]]\n",
    "    \n",
    "    if y == y_hat:\n",
    "        correct += 1\n",
    "\n",
    "    print('> test {:2d}: classifier predicted {}, actual class is {}'.format(i+1, classes[y_hat], classes[y]))\n",
    "print('\\n... Overall classification accuracy with scipy.optimize.minimize: {:.3f}'.format(correct/X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> test  1: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  2: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  3: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  4: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test  5: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  6: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  7: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test  8: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test  9: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 10: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 11: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 12: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 13: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 14: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 15: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 16: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 17: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 18: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 19: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 20: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 21: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 22: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 23: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 24: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 25: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 26: classifier predicted Iris-versicolor, actual class is Iris-versicolor\n",
      "> test 27: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 28: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "> test 29: classifier predicted Iris-virginica, actual class is Iris-virginica\n",
      "> test 30: classifier predicted Iris-setosa, actual class is Iris-setosa\n",
      "\n",
      "... Overall classification accuracy with scipy.optimize.minimize: 1.000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s finished\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "model = LogisticRegression(solver='newton-cg',\n",
    "                           multi_class='multinomial',\n",
    "                           verbose=1,\n",
    "                           max_iter=5000)\n",
    "model.fit(X, iris.target[idx_train])\n",
    "\n",
    "Y_hat = model.predict(X_test)\n",
    "\n",
    "correct = 0\n",
    "for i,y_hat in enumerate(Y_hat):\n",
    "    y = Y[i]\n",
    "    if y_hat == y:\n",
    "        correct += 1\n",
    "\n",
    "    print('> test {:2d}: classifier predicted {}, actual class is {}' \\\n",
    "          .format(i+1, classes[y_hat], classes[y]))\n",
    "\n",
    "print('\\n... Overall classification accuracy with scipy.optimize.minimize: {:.3f}' \\\n",
    "      .format(correct/X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}