chapter02_supervised-learning/logistic-regression-gluon.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Binary classification with logistic regression\n",
    "\n",
    "Over the last two tutorials we worked through how to implement a linear regression model,\n",
    "both [*from scratch*](linear-regression-scratch.ipynb)\n",
    "and [using Gluon](linear-regression-gluon.ipynb) to automate most of the repetitive work \n",
    "like allocating and initializing parameters, defining loss functions, and implementing optimizers.\n",
    "\n",
    "Regression is the hammer we reach for when we want to answer *how much?* or *how many?* questions.\n",
    "If you want to predict the number of dollars (the *price*) at which a house will be sold,\n",
    "or the number of wins a baseball team might have, \n",
    "or the number of days that a patient will remain hospitalized before being discharged,\n",
    "then you're probably looking for a regression model.\n",
    "\n",
    "Based on our experience, in industry, we're more often interested in making categorical assignments.\n",
    "*Does this email belong in the spam folder or the inbox*?\n",
    "*How likely is this custromer to sign up for subscription service?*\n",
    "When we're interested in either assigning datapoints to categories\n",
    "or assessing the *probability* that a category applies,\n",
    "we call this task *classification*. \n",
    "\n",
    "The simplest kind of classification problem is *binary classification*,\n",
    "when there are only two categories,\n",
    "so let's start there. \n",
    "Let's call our two categories the positive class $y_i=1$ and the negative class $y_i = 0$.\n",
    "Even with just two categories, and even confining ourselves to linear models, \n",
    "there are many ways we might approach the problem. \n",
    "For example, we might try to draw a line that best separates the points.\n",
    "\n",
    "![](../img/linear-separator.png)\n",
    "\n",
    "A whole family of algorithms called support vector machines pursue this approach.\n",
    "The main idea here is choose a line that maximizes the margin to the closest data points on either side of the decision boundary. \n",
    "In these approaches, only the points closest to the decision boundary (the support vectors) \n",
    "actually influence the choice of the linear separator.\n",
    "\n",
    "With neural networks, we usually approach the problem differently. \n",
    "Instead of just trying to separate the points,\n",
    "we train a probabilistic classifier which estimates,\n",
    "for each data point, the conditional probability that it belongs to the positive class. \n",
    "\n",
    "Recall that in linear regression, we made predictions of the form\n",
    "\n",
    "$$ \\hat{y} = \\boldsymbol{w}^T \\boldsymbol{x} + b. $$\n",
    "\n",
    "We are interested in asking the question *\"what is the probability that example $x$ belongs to the positive class?\"*\n",
    "A regular linear model is a poor choice here because it can output values greater than $1$ or less than $0$.\n",
    "To coerce reasonable answers from our model, \n",
    "we're going to modify it slightly,\n",
    "by running the linear function through a sigmoid activation function $\\sigma$:\n",
    "\n",
    "$$ \\hat{y} =\\sigma(\\boldsymbol{w}^T \\boldsymbol{x} + b). $$\n",
    "\n",
    "The sigmoid function $\\sigma$, sometimes called a squashing function or a *logistic* function - thus the name logistic regression - maps a real-valued input to the range 0 to 1.\n",
    "Specifically, it has the functional form:\n",
    "\n",
    "$$\\sigma(z) = \\frac{1}{1 + e^{-z}}$$\n",
    "\n",
    "Let's get our imports out of the way and visualize the logistic function using `mxnet` and `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import mxnet as mx \n",
    "from mxnet import nd, autograd, gluon\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def logistic(z):\n",
    "    return 1. / (1. + nd.exp(-z))\n",
    "    \n",
    "x = nd.arange(-5, 5, .1)    \n",
    "y = logistic(x)\n",
    "\n",
    "plt.plot(x.asnumpy(),y.asnumpy())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the sigmoid outputs a value between $0$ and $1$,\n",
    "it's more reasonable to think of it as a probability.\n",
    "Note that an input of $0$ gives a value of $.5$. \n",
    "So in the common case, where we want to predict positive whenever the probability is greater than $.5$\n",
    "and negative whenever the probability is less than $.5$,\n",
    "we can just look at the sign of $\\boldsymbol{w}^T \\boldsymbol{x} + b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary cross-entropy loss\n",
    "\n",
    "Now that we've got a model that outputs probabilities,\n",
    "we need to choose a loss function.\n",
    "When we wanted to predict *how much*, we used squared error $(y-\\hat{y})^2$,\n",
    "as our measure our model's performance. \n",
    "\n",
    "Since now we're thinking about outputing probabilities,\n",
    "one natural objective is to say that we should choose the weights \n",
    "that give the actual labels in the training data the highest probability.\n",
    "\n",
    "$$\\max_{\\theta} P_{\\theta}( (y_1, ..., y_n) | \\boldsymbol{x}_1,...,\\boldsymbol{x}_n )$$\n",
    "\n",
    "Because each example is independent of the others, and each label depends only on the features of the corresponding examples, we can rewrite the above as\n",
    "\n",
    "$$\\max_{\\theta} P_{\\theta}(y_1|\\boldsymbol{x}_1)P_{\\theta}(y_2|\\boldsymbol{x}_2) ... P(y_n|\\boldsymbol{x}_n)$$\n",
    "\n",
    "\n",
    "\n",
    "This function is a product over the examples, but in general, because we want to train by stochastic gradient descent, it's a lot easier to work with a loss function that breaks down as a sum over the training examples. \n",
    "\n",
    "$$ \\max_{\\theta} \\log P_{\\theta}(y_1|\\boldsymbol{x}_1) + ... + \\log P(y_n|\\boldsymbol{x}_n)$$\n",
    "\n",
    "Because we typically express our objective as a *loss* we can just flip the sign, giving us the *negative log probability:*\n",
    "\n",
    "$$  \\min_{\\theta} \\left(- \\sum_{i=1}^n \\log P_{\\theta}(y_i|\\boldsymbol{x}_i)\\right)$$\n",
    "\n",
    "If we interpret $\\hat{y_i}$ as the probability that the $i$-th example belongs to the positive class (i.e $y_i=1$),\n",
    "then $1 - \\hat{y_i}$ is the probability that the $i$-th example belongs to the negative class (i.e $y_i=0$). This is equivalent to saying\n",
    "\n",
    "$$ P_{\\theta}(y_i|\\boldsymbol{x}_i) = \\begin{cases}\n",
    "    \\hat{y}_i, & \\text{if } y_i = 1\\\\\n",
    "    1-\\hat{y}_i, & \\text{if } y_i = 0\n",
    "\\end{cases} $$\n",
    "\n",
    "which can be written in a more compact form\n",
    "\n",
    "$$ P_{\\theta}(y_i|\\boldsymbol{x}_i) = \\hat{y}_i^{y_i}(1-\\hat{y}_i)^{1-y_i} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we can express our learning objective as:\n",
    "\n",
    "$$ \\ell (\\boldsymbol{y}, \\boldsymbol{\\hat{y}}) =  - \\sum_{i=1}^n y_i \\log \\hat{y}_i + (1-y_i) \\log (1-\\hat{y}_i).$$\n",
    "\n",
    "\n",
    "If you're learning machine learning for the first time, that might have been too much information too quickly. \n",
    "Let's take a look at this loss function and break down what's going on more slowly. \n",
    "The loss function consists of two terms, $y_i \\log \\hat{y}_i$ and $(1-y_i) \\log (1-\\hat{y}_i)$.\n",
    "Because $y_i$ only takes values $0$ or $1$, for a given data point, one of these terms disappears. \n",
    "When $y_i$ is $1$, this loss says that we should maximize $\\log \\hat{y}_i$, giving higher probability to the *correct* answer. \n",
    "When $y_i$ is $0$, this loss function takes value $\\log (1-\\hat{y}_i)$. That says that we should maximize the value $1-\\hat{y}$ which we already know is the probability assigned to $\\boldsymbol{x}_i$ belonging to the negative class.\n",
    "\n",
    "\n",
    "Note that this loss function is commonly called *log loss* and is also commonly referred to as *binary cross entropy*. It is a special case of negative log likelihood. And it is a special case of cross-entropy, which can apply to the multi-class ($>2$) setting.\n",
    "\n",
    " \n",
    "While for linear regression, we demonstrated a completely different implementation *from scratch* and *with `gluon`*, here we're going to demonstrate how we can mix and match the two.  We'll use `gluon` for our modeling, but we'll write our loss function from scratch."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "As usual, we'll want to work out these concepts using a real dataset. This time around, we'll use the *Adult* dataset taken from the [UCI repository](http://archive.ics.uci.edu/ml/datasets/). \n",
    "The dataset was constructed by Barry Becker from 1994 census data.\n",
    "In its original form, the dataset contained $14$ features, including age, education, occupation, sex, native-country, among others. \n",
    "In this version, [hosted by National Taiwan University](https://www.csie.ntu.edu.tw/%7Ecjlin/libsvmtools/datasets/binary.html),\n",
    "the data have been re-processed to $123$ binary features each representing quantiles among the original features.\n",
    "The label is a binary indicator indicating whether the person corresponding to each row made more ($y_i = 1$) or less ($y_i = 0$) than $50,000 of income in 1994.\n",
    "The dataset we're working with contains 30,956 training examples and 1,605 examples set aside for testing.\n",
    "We can read the datasets into main memory like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_ctx = mx.cpu()\n",
    "# Change this to `mx.gpu(0) if you would like to train on an NVIDIA GPU\n",
    "model_ctx = mx.cpu()\n",
    "\n",
    "with open(\"../data/adult/a1a.train\") as f:\n",
    "    train_raw = f.read()\n",
    "\n",
    "with open(\"../data/adult/a1a.test\") as f:\n",
    "    test_raw = f.read()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data consists of lines like the following:\n",
    "\n",
    "``-1 4:1 6:1 15:1 21:1 35:1 40:1 57:1 63:1 67:1 73:1 74:1 77:1 80:1 83:1 \\n``\n",
    "\n",
    "The first entry in each row is the value of the label. \n",
    "The following tokens are the indices of the non-zero features. \n",
    "The number $1$ here is redundant. But we don't always have control over where our data comes from, so we might as well get used to mucking around with weird file formats.\n",
    "Let's write a simple script to process our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def process_data(raw_data):\n",
    "    train_lines = raw_data.splitlines()\n",
    "    num_examples = len(train_lines)\n",
    "    num_features = 123\n",
    "    X = nd.zeros((num_examples, num_features), ctx=data_ctx)\n",
    "    Y = nd.zeros((num_examples, 1), ctx=data_ctx)\n",
    "    for i, line in enumerate(train_lines):\n",
    "        tokens = line.split()\n",
    "        label = (int(tokens[0]) + 1) / 2  # Change label from {-1,1} to {0,1}\n",
    "        Y[i] = label\n",
    "        for token in tokens[1:]:\n",
    "            index = int(token[:-2]) - 1\n",
    "            X[i, index] = 1 \n",
    "    return X, Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xtrain, Ytrain = process_data(train_raw)\n",
    "Xtest, Ytest = process_data(test_raw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now verify that our data arrays have the right shapes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(Xtrain.shape)\n",
    "print(Ytrain.shape)\n",
    "print(Xtest.shape)\n",
    "print(Ytest.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the fraction of positive examples in our training and test sets. \n",
    "This will give us one nice (necessary but insufficient) sanity check that our training and test data really are drawn from the same distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(nd.sum(Ytrain)/len(Ytrain))\n",
    "print(nd.sum(Ytest)/len(Ytest))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instantiate a dataloader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "batch_size = 64\n",
    "\n",
    "train_data = gluon.data.DataLoader(gluon.data.ArrayDataset(Xtrain, Ytrain),\n",
    "                                      batch_size=batch_size, shuffle=True)\n",
    "\n",
    "test_data = gluon.data.DataLoader(gluon.data.ArrayDataset(Xtest, Ytest),\n",
    "                                      batch_size=batch_size, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "net = gluon.nn.Dense(1)\n",
    "net.collect_params().initialize(mx.init.Normal(sigma=1.), ctx=model_ctx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instantiate an optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.01})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define log loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def log_loss(output, y):\n",
    "    yhat = logistic(output)\n",
    "    return  - nd.nansum(  y * nd.log(yhat) + (1-y) * nd.log(1-yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "epochs = 30\n",
    "loss_sequence = []\n",
    "num_examples = len(Xtrain)\n",
    "\n",
    "for e in range(epochs):\n",
    "    cumulative_loss = 0\n",
    "    for i, (data, label) in enumerate(train_data):\n",
    "        data = data.as_in_context(model_ctx)\n",
    "        label = label.as_in_context(model_ctx)\n",
    "        with autograd.record():\n",
    "            output = net(data)\n",
    "            loss = log_loss(output, label)\n",
    "        loss.backward()\n",
    "        trainer.step(batch_size)\n",
    "        cumulative_loss += nd.sum(loss).asscalar()\n",
    "    print(\"Epoch %s, loss: %s\" % (e, cumulative_loss ))\n",
    "    loss_sequence.append(cumulative_loss)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the learning curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# plot the convergence of the estimated loss function \n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(num=None,figsize=(8, 6))\n",
    "plt.plot(loss_sequence)\n",
    "\n",
    "# Adding some bells and whistles to the plot\n",
    "plt.grid(True, which=\"both\")\n",
    "plt.xlabel('epoch',fontsize=14)\n",
    "plt.ylabel('average loss',fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating accuracy\n",
    "\n",
    "While the negative log likelihood gives us a sense of how well the predicted probabilities agree with the observed labels, it's not the only way to assess the performance of our classifiers. \n",
    "For example, at the end of the day, we'll often want to apply a threshold to the predicted probabilities\n",
    "in order to make hard predictions. \n",
    "For example, if we were building a spam filter, we'll need to either send the email to the spam folder or to the inbox.\n",
    "In these cases, we might not care about negative log likelihood, but instead we want know *how many errors* our classifier makes. Let's code up a simple script that calculates the accuracy of our classifier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_correct = 0.0\n",
    "num_total = len(Xtest)\n",
    "for i, (data, label) in enumerate(test_data):\n",
    "    data = data.as_in_context(model_ctx)\n",
    "    label = label.as_in_context(model_ctx)\n",
    "    output = net(data)\n",
    "    prediction = (nd.sign(output) + 1) / 2\n",
    "    num_correct += nd.sum(prediction == label)\n",
    "print(\"Accuracy: %0.3f (%s/%s)\" % (num_correct.asscalar()/num_total, num_correct.asscalar(), num_total))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This isn't too bad! A naive classifier would predict that nobody had an income greater than $50k (the majority class). This classifier would achieve an accuracy of roughly 75%.\n",
    "By contrast, our classifier gets an accuracy of .84 (results may vary a small amount on each run owing to random initializations and random sampling of the batches). \n",
    "\n",
    "By now you should have some feeling for the two most fundamental tasks in supervised learning: regression and classification. In the following chapters we'll go deeper into these problems, exploring more complex models, loss functions, optimizers, and training schemes. We'll also look at more interesting datasets. And finally, in the following chapters we'll also look more advanced problems where we want, for example, to predict more structured objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next:\n",
    "[Softmax regression from scratch](softmax-regression-scratch.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For whinges or inquiries, [open an issue on  GitHub.](https://github.com/zackchase/mxnet-the-straight-dope)"
   ]
  }
 ],
 "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.2"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}