practicals/classification/data414_classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical: Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Herman Kamper, 2023"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from scipy.spatial import distance\n",
    "from sklearn import datasets, linear_model, neighbors, preprocessing\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. $K$-nearest neighbours (KNN) classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we apply $K$-nearest neighbours (KNN) classification to the `Default` dataset in order to predict whether a person will default on his or her credit card payment, based on their annual income and monthly credit card balance [ISL, $\\S$4.1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Pre-process and visualise the data\n",
    "\n",
    "\n",
    "**Questions:**\n",
    "\n",
    "- Load the data from `default.csv`. You can use `np.loadtxt`.\n",
    "- We will only be using the `default`, `balance` and `income` columns, so you can disregard the `student` column when loading the data.\n",
    "- Split the data so that the first 8000 entries are used as training, the next 1000 as validation, and the last 1000 as testing.\n",
    "- Plot the training set with `balance` on the $x$-axis and `income` on the $y$-axis, distinguishing the individuals who default on their credict payments from those who do not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you haven't downloaded the data already with the notebook: Download data\n",
    "import urllib\n",
    "urllib.request.urlretrieve(\n",
    "    \"https://raw.githubusercontent.com/kamperh/data414/main/practicals/classification/default.csv\", \"default.csv\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Implement KNN\n",
    "\n",
    "In KNN classification, a test point is classified by determining the $K$ closest points in the training data, and then assigning it the most common class among these neighbours.\n",
    "\n",
    "**Question:** Using NumPy, write a function `knn_predict` taking arguments `X_new`, `X_train`, `y_train` and `K`. The `X_train` and `y_train` arguments should be the training set inputs and outputs, respectively. The `X_new` argument is an unseen set on which predictions should be performed; it's shape will be `[N_new, D]` where `N_new` is the number of unseen points and `D` is the dimensionality of a single input (2 in this case). `K` is the number of neighbours that should be considered to make the KNN prediction. Use the `Euclidean` distance metric. The function should return a single vector with `N_new` predictions of either 0 (do not default) or 1 (default).  \n",
    "*Hint:* Several functions can be useful here, including `scipy.spatial.distance.cdist`, `np.argsort`, `np.mean` and `np.where`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer: Complete the function\n",
    "\n",
    "def knn_predict(X_new, X_train, y_train, K=1):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions:**\n",
    "- Apply your `knn_predict` function to the validation data with $K = 1$.\n",
    "- Calculate the validation accuracy.\n",
    "- Apply your `knn_predict` function to the first 1000 entries in the training data, again with $K = 1$. Before running the code, what do you guess this training accuracy will be?\n",
    "- Change the value of $K$ and observe the effects on both the validation and training accuracies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$K$ can be seen as a regularisation parameter, with a higher $K$ resulting in a smoother decision boundary while a lower $K$ results in a prediction which is closer to the training data.\n",
    "\n",
    "**Question:** To illustrate this, vary $K$ from 1 to 15 and calculate the validation and training accuracies as  above. Then plot the validation and training accuracies on a single plot with $K$ on the $x$-axis and accuracy on the $y$-axis. What do you observe? What value of $K$ should you use on the held-out test data?  \n",
    "*Note:* Depending on your implementation of KNN, it might take a while to calculate all these accuracies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn also has a `KNeighborsClassifier` class in `sklearn.neighbors`, which you can use from this point onwards. This is (probably) a more efficient implementation than your own, includes a number of extensions, and allows for data with more than just two classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, you will build a logistic regression model to predict whether a student is admitted to a university or not. Suppose that you are the administrator of the university and you have to determine each applicant's chance of admission based on their results on two exams. You have historical data from previous applicants that you can use as training data. For each training example, you have the applicant’s scores on the two exams and the admissions decision. Your task is to build a classifier that estimates an applicant’s probability of admission based the scores from these two exams."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Pre-process and visualise the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions:**\n",
    "\n",
    "- Load the data from `admissions.csv`.\n",
    "- Plot the data, distinguishing the students who were admitted from those who were not.\n",
    "- Normalise the data to have zero mean and unit standard deviation. You can use `(X - X.mean(axis=0))/X.std(axis=0)`. Why would we do this?\n",
    "- Let us assume that the first 80 lines are students that were admitted/not admitted in the past (our training data), that the next 10 lines are data we can use for validation, and that the last 10 lines are data that we want to use for final testing. Split the data accordingly.\n",
    "- Plot the normalised training set, again distinguishing the students who were admitted from those who were not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you haven't downloaded the data already with the notebook: Download data\n",
    "import urllib\n",
    "urllib.request.urlretrieve(\n",
    "    \"https://raw.githubusercontent.com/kamperh/data414/main/practicals/classification/admissions.csv\",\n",
    "    \"admissions.csv\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Implement logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For logistic regression, we use the following prediction function:\n",
    "\n",
    "$$f(\\mathbf{x}; \\mathbf{w}) = \\sigma(\\mathbf{w}^\\top\\mathbf{x}) = \\frac{1}{1+e^{-\\mathbf{w}^\\top\\mathbf{x}}}$$\n",
    "\n",
    "Assuming we have binary labels $y \\in \\{ 0, 1 \\}$, we interpret this function as $f(\\mathbf{x}; \\mathbf{w}) = P(y=1 \\,|\\, \\mathbf{x}; \\mathbf{w})$. We fit it using maximum likelihood estimation, specifically by minising the negative log likelihood:\n",
    "\n",
    "$$J(\\mathbf{w}) = - \\sum_{n = 1}^N \\left[ y^{(n)}\\log f(\\mathbf{x}^{(n)}; \\mathbf{w}) +  (1 - y^{(n)})\\log(1 - f(\\mathbf{x}^{(n)}; \\mathbf{w})) \\right]$$\n",
    "\n",
    "One way to optimise this loss is by using gradient descent:\n",
    "\n",
    "$$\\mathbf{w} \\leftarrow \\mathbf{w} - \\eta \\frac{\\partial J (\\mathbf{w})}{\\partial \\mathbf{w}}$$\n",
    "\n",
    "In this case, the gradient is given by \n",
    "\n",
    "$$\\frac{\\partial J (\\mathbf{w})}{\\partial \\mathbf{w}} = - \\sum_{n=1}^N (y^{(n)} - f(\\mathbf{x}^{(n)}; \\mathbf{w})) \\mathbf{x}^{(n)}$$\n",
    "\n",
    "**Questions:**\n",
    "\n",
    "- Complete the code below to calculate the negative log likelihood for a given dataset.\n",
    "- Complete the code below to perform logistic regression on he admissions data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer: Complete the `get_loss` function\n",
    "\n",
    "# Add bias terms\n",
    "X_train_bias = np.c_[np.ones((X_train.shape[0], 1)), X_train]\n",
    "X_val_bias = np.c_[np.ones((X_val.shape[0], 1)), X_val]\n",
    "X_test_bias = np.c_[np.ones((X_test.shape[0], 1)), X_test]\n",
    "\n",
    "# Function for calculating NLL\n",
    "sigmoid = lambda x: 1./(1. + np.exp(-x))\n",
    "def get_nll_loss(X, y, w):\n",
    "    loss = 0.0\n",
    "    # Answer: Add your own code here\n",
    "    return loss\n",
    "\n",
    "# Calculate NLL for random weights\n",
    "w = np.random.randn(X_train.shape[1])\n",
    "loss = get_nll_loss(X_train_bias, y_train, w)\n",
    "print(\"NLL for random weights: {:.4f}\".format(loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer: Add the gradient calculation below\n",
    "\n",
    "# Training parameters\n",
    "n_iterations = 20\n",
    "learning_rate = 0.1\n",
    "\n",
    "# Gradient descent\n",
    "w = np.random.randn(X_train_bias.shape[1])  # initialise weights\n",
    "for iteration in range(n_iterations):\n",
    "    gradients = 0\n",
    "    \n",
    "    # Answer: Add your own code here\n",
    "    \n",
    "    w = w - learning_rate*gradients\n",
    "    train_loss = get_nll_loss(X_train_bias, y_train, w)\n",
    "    val_loss = get_nll_loss(X_val_bias, y_val, w)\n",
    "    print(\n",
    "        \"Epoch {:2d}: training loss {:.2f}, validation loss {:.2f}\".format(\n",
    "        iteration, train_loss, val_loss)\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Plot the decision boundary and make predictions\n",
    "\n",
    "**Question:** Plot the resulting decision boundary on top of the training data.  \n",
    "*Hint:* The decision boundary is where $w_0 + w_1x_1 + w_2x_2 = 0$. Write this equation as $y = mx + c$ (with $x$ and $y$ now indicating the plot axes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions:**\n",
    "\n",
    "- Calculate the validation accuracy.\n",
    "- You can go back and change the number of training iterations and learning rate. Note how this affects the validation accuracy and try to find values that give consistent performance.\n",
    "- Only once you have decided on the final values for these hyperparameters, can you apply your model to the test set. Calculate the test accuracy for your final model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Optional question:** You can go back to the start of Section 2.1 and repeat the entire exercise directly on the unnormalised data. Although the scales of the $x_1$ and $x_2$ features are similar, you will observe that the scales of $w_0$, $w_1$ and $w_2$ are actually very different now, making the algorithm much more sensitive to the initialisation of $\\mathbf{w}$ since the same learning rate $\\eta$ is used for all three weights. With a bit of effort, you could still get it to work, but you will probably have to initialise $\\mathbf{w}$ very carefully (maybe try an initial $\\mathbf{w} = \\begin{bmatrix} -50 & 0 & 0 \\end{bmatrix}^\\top$ and use a very small learning rate such as $\\eta = 0.00001$ and play around with these values).\n",
    "\n",
    "Scikit-learn also has a `LogisticRegression` class in `sklearn.linear_model`, which you can use from this point onwards. This is (probably) a more efficient implementation than your own, includes a number of extensions, and allows for data with more than just two classes. It also provides the option of specifying different optimisation methods from gradient descent, which might be more robust to some of the initialisation and update issues described above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Non-linear decision boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During quality assurance in a microchip fabrication plant, each microchip goes through various checks to ensure it functions correctly. Suppose you are the product manager of the factory and you have the results for some microchips on two different quality checks. From these two checks, you would like to determine whether a new microchip should be accepted or rejected. I.e., you have a dataset of check results on past microchips, from which you can build a prediction model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. KNN classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below loads the data from `microchip.csv`, applies KNN classification to the data, and plots the resulting decision boundary.\n",
    "\n",
    "**Question:** Change the number of nearest neighbours `K` and observe the change in the decision boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you haven't downloaded the data already with the notebook: Download data\n",
    "import urllib\n",
    "urllib.request.urlretrieve(\n",
    "    \"https://raw.githubusercontent.com/kamperh/data414/main/practicals/classification/microchip.csv\",\n",
    "    \"microchip.csv\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load and plot the data\n",
    "data = np.loadtxt(\"microchip.csv\", delimiter=\",\")\n",
    "X = data[:, :2]\n",
    "y = data[:, 2]\n",
    "plt.scatter(X[y==0, 0], X[y==0, 1], c=\"C0\", marker=\"s\", edgecolor=\"white\")\n",
    "plt.scatter(X[y==1, 0], X[y==1, 1], c=\"C1\", marker=\"o\", edgecolor=\"white\")\n",
    "plt.xlabel(\"Check 1\")\n",
    "plt.ylabel(\"Check 2\")\n",
    "\n",
    "# KNN classification\n",
    "K = 2\n",
    "knn = neighbors.KNeighborsClassifier(n_neighbors=K)\n",
    "knn.fit(X, y)\n",
    "\n",
    "# Create mesh\n",
    "def make_meshgrid(x, y, h=.01):\n",
    "    \"\"\"\n",
    "    Based on: https://stackoverflow.com/questions/51297423\n",
    "    \"\"\"\n",
    "    x_min, x_max = x.min() - h, x.max() + h\n",
    "    y_min, y_max = y.min() - h, y.max() + h\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "    return xx, yy\n",
    "X0, X1 = X[:, 0], X[:, 1]  # modify to match your training data\n",
    "xx, yy = make_meshgrid(X0, X1)\n",
    "X_new = np.c_[xx.ravel(), yy.ravel()]\n",
    "\n",
    "# Plot the decision boundary\n",
    "predictions = knn.predict(X_new).reshape(xx.shape)\n",
    "plt.contour(xx, yy, predictions, levels=[0.5], linewidths=[1.5], colors=[\"black\"]) #, linestyles=[\"dashed\"])\n",
    "# plt.tight_layout()\n",
    "\n",
    "plt.savefig(\"fig/question3.1.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Regularised logistic regression with polynomial features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will now apply logistic regression to this data using the `LogisticRegression` class from `sklearn.linear_model`. \n",
    "\n",
    "The loss for regularised logistic regression includes an additional term:\n",
    "\n",
    "$$J(\\mathbf{w}) = - \\sum_{n = 1}^N \\left[ y^{(n)}\\log f(\\mathbf{x}^{(n)}; \\mathbf{w}) +  (1 - y^{(n)})\\log(1 - f(\\mathbf{x}^{(n)}; \\mathbf{w})) \\right] + \\lambda \\sum_{d = 1}^D w_d^2$$\n",
    "\n",
    "The $\\lambda$ hyperparameter controls the extent of the regularisation. For the `LogisticRegression` class, the `C` argument is proportional to $1/\\lambda$, i.e. less regularisation is applied with a large `C`. Regularisation can be turned off by setting `penalty=\"none\"`. Use the `solver=\"lbfgs\"` optimiser in this section.\n",
    "\n",
    "**Questions:**\n",
    "\n",
    "- Using scikit-learn classes and functions, apply standard logistic regression to the microchip data and plot the decision boundary.\n",
    "- Experiment with no regularisation (`penalty=\"none\"`) and with adding more and less regularisation (small and large values for `C`).\n",
    "- Does regularisation help?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will now use the `PolynomialFeatures` class from `sklearn.preprocessing` to fit a logistic regression model with $6^{\\text{th}}$ order polynomial features, i.e.,\n",
    "\n",
    "$$\\boldsymbol{\\phi}(\\mathbf{x}) = \\begin{bmatrix} 1 & x_1 & x_2 & x_1^2 & x_1x_2 & x_2^2 & x_1^3 & \\ldots & x_1x_2^5 & x_2^6 \\end{bmatrix}^\\top$$\n",
    "\n",
    "**Questions:**\n",
    "\n",
    "- First fit the logistic regression models without any regularisation, and plot the decision boundary.\n",
    "- Add some regularisation, and plot the decision boundary again.\n",
    "- Increase regularisation (small `C`, e.g. `0.01`) and observe the effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Answer: Add code and new cells here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Acknowledgements\n",
    "\n",
    "- Section 1 is based on examples from [An Introduction to Statistical Learning (ISL)](http://faculty.marshall.usc.edu/gareth-james/ISL/).\n",
    "- Sections 2 and 3 are based on examples from the [Coursera Machine Learning](https://www.coursera.org/learn/machine-learning) course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}