Lab 03 - Logistic regression.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Mini-lab 03 - Logistic regression.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lyeWNnlxTI22"
      },
      "source": [
        "# Lab 3: Logistic regression\n",
        "\n",
        "Today, you will be experimenting with logistic regression. As described in the lectures, logistic regression is simply a regression model that is passed through a sigmoid function. This forces the output to be in the range [0, 1], which is why our binary labels must be encoded as {0, 1} and also why we can pretend that this output is a probability.\n",
        "\n",
        "Below, you will find code for making data, classifiers, and optimizing the parameters. You will also be provided with some plotting code."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZQNI2sNfmA79"
      },
      "source": [
        "import numpy as np\n",
        "from sklearn.datasets import make_classification, make_moons, make_blobs  # Toy data\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "from mpl_toolkits.mplot3d.axes3d import Axes3D\n",
        "# from sklearn.metrics import confusion_matrix"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eOM_T-7SWt18"
      },
      "source": [
        "## Data\n",
        "\n",
        "Let's make some data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1g-BxSZ_WuCz"
      },
      "source": [
        "X, y = make_blobs(n_samples=200, centers=2, cluster_std=.7, n_features=2, random_state=0)\n",
        "X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1) # Add pseudo input\n",
        "X[:5], y[:5]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2EVcr8ADXp4d"
      },
      "source": [
        "X, y = make_moons(n_samples=200, noise=.1, random_state=0)\n",
        "X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1) # Add pseudo input\n",
        "X[:5], y[:5]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uFHRDNQeXpWE"
      },
      "source": [
        "We want to easily be able to plot the data. The plotting function can take a classifier as a parameter. Python allows us to use any object as a parameter. However, we don't have one yet, but must supply X and y."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gLeJv2oOX1NK"
      },
      "source": [
        "def plot_classification(X, y, classifier=None):\n",
        "  plt.figure(figsize=(5, 5), dpi=100)\n",
        "  markers = ['o', 's', 'd']\n",
        "  for j, label in enumerate(np.unique(y)):\n",
        "    if X.shape[1] == 2:\n",
        "      # No pseudo input\n",
        "      plt.scatter(X[y==label, 0], X[y==label, 1], c=y[y==label], s=50, marker=markers[j], edgecolor='k', cmap='coolwarm', vmin=y.min(), vmax=y.max())\n",
        "    else:\n",
        "      # X includes pseudo inputs\n",
        "      plt.scatter(X[y==label, 1], X[y==label, 2], c=y[y==label], s=50, marker=markers[j], edgecolor='k', cmap='coolwarm', vmin=y.min(), vmax=y.max())\n",
        "  if classifier is not None:\n",
        "    # Create a plotting mesh\n",
        "    a = plt.axis()          # The mesh will adapt plotted data\n",
        "    XX, YY = np.meshgrid(np.linspace(a[0], a[1], 20),\n",
        "                        np.linspace(a[2], a[3], 20))\n",
        "    mesh_points = [np.vstack(XX.ravel()), np.vstack(YY.ravel())]  # Reshape the x and y coordinates to columns\n",
        "    mesh_points.insert(0, np.ones(mesh_points[0].shape))  # Add pseudo input\n",
        "    mesh_points = np.concatenate(mesh_points, axis=1)     # Convert to numpy array\n",
        "    # Run classification\n",
        "    Z = np.apply_along_axis(classifier, 1, mesh_points) \n",
        "    # Plot, contourf will interpolate the colour between our mesh points\n",
        "    plt.contourf(XX, YY, Z.reshape(XX.shape), alpha=.5, cmap='coolwarm', zorder=-1)\n",
        "  plt.tight_layout()\n",
        "  plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plot_classification(X, y)"
      ],
      "metadata": {
        "id": "t5TLhuZzxhcd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s20fb3MZYyGx"
      },
      "source": [
        "## Modelling\n",
        "\n",
        "We can now define the constituent functions and make a model. If you can't remember what a pseudo input is, this might not make sense."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H9konM9xCwUL"
      },
      "source": [
        "def regression(X, theta):\n",
        "  \"\"\"Regression function\n",
        "  Assumes a pseudo input and preforms a dot product between each row of the \n",
        "  data matrix and the parameter vector.\"\"\"\n",
        "  if X.ndim==1:\n",
        "    X = np.vstack(X).T\n",
        "  return X.dot(np.vstack(theta)).ravel()\n",
        "\n",
        "def sigmoid(z):\n",
        "  \"\"\"Sigmoid function\n",
        "  Runs element wise on the input data.\"\"\"\n",
        "  return 1/(1+np.exp(-z))\n",
        "\n",
        "def model(x, theta):\n",
        "  \"\"\"Logistic regression\"\"\"\n",
        "  return sigmoid(regression(x, theta)).ravel()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qYYZR15XnMGu"
      },
      "source": [
        "We can plot the sigmoid to verify that we got it right.\n",
        "\n",
        "Note that the sigmoid function never becomes flat at the extremes (it does come close, though). Think about why?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2WoovN4Qh23n"
      },
      "source": [
        "plt.figure()\n",
        "x = np.linspace(-10, 10, 50)\n",
        "plt.plot(x, sigmoid(x))\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "c0zPRxu-asKt"
      },
      "source": [
        "We can select our initial $\\theta$ in any way we like. Here are two alterantives. We now have everything we need to assemble the model and classify."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jGL0Wt_AaCPg"
      },
      "source": [
        "theta = np.zeros(3)\n",
        "# theta = np.random.normal(size=3)\n",
        "\n",
        "regression(X, theta)[:5], sigmoid(regression(X, theta))[:5], model(X, theta)[:5]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PiflukQ8hkU1"
      },
      "source": [
        "## Choosing parameters\n",
        "\n",
        "To do optimization, we need some loss functions. We only require one, to be strict, but you should change between them to see what happens."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xyW38_AOhjdq"
      },
      "source": [
        "def squared_loss(y, y_hat):\n",
        "  \"\"\"Squared loss\n",
        "  This is normal used for regression, but also works here.\"\"\"\n",
        "  return np.sum((y - y_hat)**2)\n",
        "\n",
        "def logistic_loss(y, y_hat):\n",
        "  \"\"\"Logistic loss\n",
        "  This is the standard loss function for logistic regression.\"\"\"\n",
        "  return np.sum(-y*np.log(y_hat)-(1-y)*np.log(1-y_hat))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xRJ5NkuXEczp"
      },
      "source": [
        "We can try the same optimization algorithms as last time."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qEj4iW6HF2XD"
      },
      "source": [
        "def coordinate_descent(function_to_minimize, initial_theta, n_iterations=50, step=.25):\n",
        "  \"\"\"Coordinate descent (look at wikipedia for definition)\"\"\"\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    theta_changed_during_iter = False\n",
        "    for dim in range(len(initial_theta)):\n",
        "      for direction in [-1, 1]:\n",
        "        new_theta = theta.copy()\n",
        "        new_theta[dim] += direction*step\n",
        "        while function_to_minimize(new_theta) < function_to_minimize(theta):\n",
        "          theta = new_theta\n",
        "          theta_changed_during_iter = True\n",
        "          history.append(theta)\n",
        "          function_values.append(function_to_minimize(new_theta))\n",
        "          new_theta = theta.copy()\n",
        "          new_theta[dim] += direction*step\n",
        "    if not theta_changed_during_iter:\n",
        "      step /= 2\n",
        "      # break\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def finite_difference_gradient_descent(function_to_minimize, initial_theta, n_iterations=20, gamma=.01):\n",
        "  \"\"\"Gradient descent optimization, where the gradient is approximated using\n",
        "  finite difference. Look at scipy.optimize.approx_fprime to learn more.\"\"\"\n",
        "  from scipy.optimize import approx_fprime\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    gradient = approx_fprime(theta, function_to_minimize, epsilon=.001)\n",
        "    theta = theta - gamma*gradient\n",
        "    history.append(theta)\n",
        "    function_values.append(function_to_minimize(theta))\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def random_changes(function_to_minimize, initial_theta, n_iterations=200, noise_sd=.1):\n",
        "  \"\"\"Adds random noise to the parameter vector in the hope of \n",
        "  improving the loss.\"\"\"\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    new_theta = theta + np.random.normal(scale=noise_sd, size=theta.shape)\n",
        "    if function_to_minimize(new_theta) < function_to_minimize(theta):\n",
        "      theta = new_theta\n",
        "    history.append(theta.copy())\n",
        "    function_values.append(function_to_minimize(theta))\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def scipy_minimize(function_to_minimize, initial_theta):\n",
        "  \"\"\"Optimization using external library. Look at the manual of \n",
        "  scipy.optimize.minimize for specifics.\"\"\"\n",
        "  from scipy.optimize import minimize\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  def store_theta_during_optimization(current_theta):\n",
        "    history.append(current_theta)\n",
        "    function_values.append(function_to_minimize(current_theta))\n",
        "  res = minimize(lambda t: function_to_minimize(t), initial_theta, \n",
        "                 callback=store_theta_during_optimization)\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Wyyj7Em_gPeg"
      },
      "source": [
        "# theta, loss, history = coordinate_descent(function_to_minimize=lambda theta: logistic_loss(y, model(X, theta)), \n",
        "#                                           initial_theta=np.random.normal(size=3),\n",
        "#                                           n_iterations=200, step=.25)\n",
        "\n",
        "# theta, loss, history = finite_difference_gradient_descent(function_to_minimize=lambda theta: logistic_loss(y, model(X, theta)), \n",
        "#                                                           initial_theta=np.random.normal(size=3), \n",
        "#                                                           n_iterations=60, gamma=.01)\n",
        "\n",
        "# theta, loss, history = random_changes(function_to_minimize=lambda theta: logistic_loss(y, model(X, theta)),\n",
        "#                                       initial_theta=np.random.normal(size=3),\n",
        "#                                       n_iterations=200, noise_sd=.1)\n",
        "\n",
        "\n",
        "# DON'T WRITE CODE LIKE THIS !!1!\n",
        "# def f(theta):\n",
        "#   y_hat = model(X, theta)\n",
        "#   value = logistic_loss(y, y_hat)\n",
        "#   return value\n",
        "\n",
        "theta, loss, history = scipy_minimize(function_to_minimize=lambda theta: logistic_loss(y, model(X, theta)),\n",
        "                                      initial_theta=np.random.normal(size=3))\n",
        "\n",
        "plt.figure()\n",
        "plt.plot(loss)\n",
        "plt.ylabel(\"Loss\")\n",
        "plt.xlabel(\"Iteration\")\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "msnmoZmfmBm7"
      },
      "source": [
        "plot_classification(X, y, lambda x: model(x, theta))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WawV80hSF67y"
      },
      "source": [
        "import ipywidgets\n",
        "\n",
        "def plot_wrapper(iteration):\n",
        "  plot_classification(X, y, lambda x: model(x, history[iteration]))\n",
        "interact_plot = ipywidgets.interact(plot_wrapper,\n",
        "                                    iteration = ipywidgets.Play(min=0, max=len(history)-1, step=max(1, len(history)//100), value=0));\n",
        "output = interact_plot.widget.children[-1] # This should prevent flickering\n",
        "output.layout.height = '500px'\n",
        "# Note that the play button doesn't always show up in colab. The clickable area is still there though (leftmost gray area)."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2RFg0Vg2oRSC"
      },
      "source": [
        "Making predictions is easy from what we have. By pushing the model output through np.round, we get labels in {0, 1}. This makes is possible to calculate the accuracy."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m0qbkyu6_unK"
      },
      "source": [
        "prediction = np.round(model(X, theta))\n",
        "correct = np.equal(prediction, y)\n",
        "print(\"Accuracy: %.1f%%\" % (100*np.sum(correct)/len(correct)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QSpRUGD9oKdT"
      },
      "source": [
        "## More complex data\n",
        "\n",
        "Logistic regression can be extended to handle multi-class classification. This will be discussed in a future lecture. Can you think of an extension to logistic regression that would give a multi class output?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K29mPutUoL-5"
      },
      "source": [
        "X, y = make_classification(n_samples=200, n_features=2, n_informative=2, \n",
        "                           n_redundant=0, n_repeated=0, n_classes=3, class_sep=1.8,\n",
        "                           n_clusters_per_class=1, random_state=0)\n",
        "X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1) # Add pseudo input\n",
        "\n",
        "plot_classification(X, y)\n",
        "\n",
        "X[:5], y[:5]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jMibt3jxoO98"
      },
      "source": [
        "We now have 3 classes. Can this be done using binary classification?\n",
        "\n",
        "There are several ways of dealing with this. We can try out to simply add two logistic regressions together and, the data willing, we might get a good classification. To do this we have to create a longer $\\theta$ and slice it."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x1-gQwaqHemu"
      },
      "source": [
        "theta = np.random.normal(size=X.shape[1]*2)\n",
        "theta, theta[:len(theta)//2], theta[len(theta)//2:]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i1uOxEnlrVjU"
      },
      "source": [
        "Investigate this model. Are there other ways of setting this up? Is this easilly optimized? Does logistic loss work here?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "alZV9nxqEGct"
      },
      "source": [
        "multilabel_model = lambda X, theta: sigmoid(regression(X, theta[:len(theta)//2])) + sigmoid(regression(X, theta[len(theta)//2:]))\n",
        "multilabel_model(X, theta)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mI3UAqe1LL-c"
      },
      "source": [
        "theta, loss, history = random_changes(function_to_minimize=lambda theta: squared_loss(y, multilabel_model(X, theta)),\n",
        "                                      initial_theta=np.random.normal(size=X.shape[1]*2),\n",
        "                                      n_iterations=200, noise_sd=.1)\n",
        "\n",
        "plt.figure()\n",
        "plt.plot(loss)\n",
        "plt.ylabel(\"Loss\")\n",
        "plt.xlabel(\"Iteration\")\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jGgR0FlALWZB"
      },
      "source": [
        "prediction = np.round(multilabel_model(X, theta))\n",
        "correct = np.equal(prediction, y)\n",
        "print(\"Accuracy: %.1f%%\" % (100*np.sum(correct)/len(correct)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A_9m_5-HM0UQ"
      },
      "source": [
        "plot_classification(X, y, lambda x: multilabel_model(x, theta)) # This will not work with the iris data"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "How does the decision boundaries look? Can you see how they are composed out of two linear boundaries?"
      ],
      "metadata": {
        "id": "q1ka9NtW1tI-"
      }
    }
  ]
}