{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<h2><center> <span style = \"font-family: Babas; font-size: 2em;\"> Implementing Logistic Regression from Scratch </span> </center></h2>\n<h4><center> <span style = \"font-family: Babas; font-size: 2em;\"> Sugata Ghosh </span> </center></h4>","metadata":{}},{"cell_type":"markdown","source":"### Contents\n\n- [Introduction](#Introduction)\n- [Logistic Function](#Logistic-Function)\n- [Log Loss](#Log-Loss)\n- [Cost Function](#Cost-Function)\n- [Gradient Descent](#Gradient-Descent)\n- [Preprocessing](#Preprocessing)\n- [Model Fitting](#Model-Fitting)\n- [Prediction and Evaluation](#Prediction-and-Evaluation)\n- [Regularization](#Regularization)\n- [Acknowledgements](#Acknowledgements)\n- [References](#References)","metadata":{}},{"cell_type":"code","source":"# Importing libraries\nimport time, psutil, os, math\nfrom tqdm.contrib import itertools\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set_theme()\nfrom sklearn.model_selection import train_test_split","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:08.533130Z","iopub.execute_input":"2022-07-25T09:39:08.533796Z","iopub.status.idle":"2022-07-25T09:39:09.931480Z","shell.execute_reply.started":"2022-07-25T09:39:08.533706Z","shell.execute_reply":"2022-07-25T09:39:09.930166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Runtime and memory usage\nstart = time.time()\nprocess = psutil.Process(os.getpid())","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:09.933580Z","iopub.execute_input":"2022-07-25T09:39:09.934025Z","iopub.status.idle":"2022-07-25T09:39:09.939648Z","shell.execute_reply.started":"2022-07-25T09:39:09.933967Z","shell.execute_reply":"2022-07-25T09:39:09.938876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Introduction","metadata":{}},{"cell_type":"markdown","source":"**Classification.** In [statistics](https://en.wikipedia.org/wiki/Statistics) and [machine learning](https://en.wikipedia.org/wiki/Machine_learning), [classification](https://en.wikipedia.org/wiki/Statistical_classification) refers to a type of [supervised learning](https://en.wikipedia.org/wiki/Supervised_learning). For this task, training data with known class labels are given and is used to develop a [classification rule](https://en.wikipedia.org/wiki/Classification_rule) for assigning new unlabeled data to one of the classes. A special case of the task is [binary classification](https://en.wikipedia.org/wiki/Binary_classification), which involves only two classes. Some examples:\n\n- Classifying an email as `spam` or `non-spam`\n- Classifying a tumor as `benign` or `malignant`\n\nThe algorithms that sort unlabeled data into labeled classes are called *classifiers*. Loosely speaking, the [sorting hat](https://en.wikipedia.org/wiki/Magical_objects_in_Harry_Potter#Sorting_Hat) from [Hogwarts](https://en.wikipedia.org/wiki/Hogwarts) can be thought of as a classifier that sorts incoming students into four distinct houses. In real life, some common classifiers are [logistic regression](https://en.wikipedia.org/wiki/Logistic_regression), [k-nearest neighbors](https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm), [decision tree](https://en.wikipedia.org/wiki/Decision_tree_learning), [random forest](https://en.wikipedia.org/wiki/Random_forest), [support vector machine](https://en.wikipedia.org/wiki/Support-vector_machine), [naive Bayes](https://en.wikipedia.org/wiki/Naive_Bayes_classifier), [linear discriminant analysis](https://en.wikipedia.org/wiki/Linear_discriminant_analysis), [stochastic gradient descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent), [XGBoost](https://en.wikipedia.org/wiki/XGBoost), [AdaBoost](https://en.wikipedia.org/wiki/AdaBoost) and [neural networks](https://en.wikipedia.org/wiki/Artificial_neural_network).\n\n**The purpose of the notebook.** Many advanced libraries, such as [scikit-learn](https://en.wikipedia.org/wiki/Scikit-learn), make it possible for us to train various models on labeled [training data](https://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets#Training_data_set), and predict on unlabeled [test data](https://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets#Test_data_set), with a few lines of codes. While it is very convenient for day-to-day practice, it does not give insight into the details of what really happens underneath, when we run those codes. In the present notebook, we implement a logistic regression model manually from scratch, without using any advanced library, to understand how it works in the context of binary classification. The basic idea is to segment the computations into pieces, and write functions to compute each piece in a sequential manner, so that we can build a function on the basis of the previously defined functions. Wherever applicable, we have complemented a function which is constructed using for loops, with a much faster vectorized implementation of the same.\n\n**A problem from particle physics.** We have chosen the particular problem posed in [this competition](https://www.kaggle.com/competitions/higgs-boson). In [particle physics](https://en.wikipedia.org/wiki/Particle_physics), an event refers to the results just after a [fundamental interaction](https://en.wikipedia.org/wiki/Fundamental_interaction) takes place between [subatomic particles](https://en.wikipedia.org/wiki/Subatomic_particle), occurring in a very short time span, at a well-localized region of space. The problem is to classify an event produced in a particle accelerator as *background* or *signal*, based on relevant feature variables. A background event is explained by the existing theories and previous observations. A signal event, however, indicates a process that cannot be described by previous observations and leads to the potential discovery of a new particle. More on this problem is detailed in the introduction section of [this notebook](https://www.kaggle.com/code/sugataghosh/higgs-boson-event-detection-part-1-eda).\n\n**Data.** [The dataset](https://www.kaggle.com/competitions/higgs-boson/data), provided with the competition, has been built from official ATLAS full-detector simulation. The simulator has two parts. In the first, random proton-proton collisions are simulated based on the knowledge that we have accumulated on particle physics. It reproduces the random microscopic explosions resulting from the proton-proton collisions. In the second part, the resulting particles are tracked through a virtual model of the detector. The process yields simulated events with properties that mimic the statistical properties of the real events with additional information on what has happened during the collision, before particles are measured in the detector. Information on $250000$ events are included int he dataset. For each event, it has information on $31$ features ($2$ integer-type features and $29$ float-type features). Additionally, the dataset contains the object-type target variable `labels` and float-type variable `weights`. The target variable can take two possible values: $b$ (indicating a background event) and $s$ (indicating a signal event).","metadata":{}},{"cell_type":"code","source":"# Loading the data\ndata = pd.read_csv('../input/higgs-boson/training.zip')\nprint(pd.Series({\"Memory usage\": \"{:.4f} MB\".format(data.memory_usage().sum()/(1024*1024)),\n                 \"Dataset shape\": \"{}\".format(data.shape)}).to_string())\ndata.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:09.940594Z","iopub.execute_input":"2022-07-25T09:39:09.941315Z","iopub.status.idle":"2022-07-25T09:39:11.763382Z","shell.execute_reply.started":"2022-07-25T09:39:09.941283Z","shell.execute_reply":"2022-07-25T09:39:11.762133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Synopsis of the data**\n\n- Number of observations: $250000$\n- Number of columns: $33$\n- Number of integer columns: $2$\n- Number of float columns: $30$\n- Number of object columns: $1$\n- Number of duplicate observations: $0$\n- Constant columns: None\n- Number of columns with missing values: $0$\n- Memory Usage: $62.94$ MB","metadata":{}},{"cell_type":"markdown","source":"# Logistic Function","metadata":{}},{"cell_type":"markdown","source":"A function $g: \\mathbb{R} \\to \\mathbb{R}$ is said to be a [sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function) if it has the following properties:\n- It is [bounded](https://en.wikipedia.org/wiki/Bounded_function)\n- It is [differentiable](https://en.wikipedia.org/wiki/Differentiable_function)\n- It has nonnegative derivative at each point\n- It has exactly one [inflection point](https://en.wikipedia.org/wiki/Inflection_point)\n\nAn example of a sigmoid function is the standard [logistic function](https://en.wikipedia.org/wiki/Logistic_function) (sometimes simply referred to as the *sigmoid*), which is given by\n\n$$ g(x) = \\frac{1}{1+e^{-x}}, $$\n\nfor $x \\in \\mathbb{R}$. The next two code blocks construct and plot this function.","metadata":{}},{"cell_type":"code","source":"# Logistic function\ndef logistic(x):\n    \"\"\"\n    Computes the logistic function applied to an input scalar/array\n    Args:\n        x (scalar/ndarray): scalar or numpy array of any size\n    Returns:\n        y (scalar/ndarray): logistic function applied to x, has the same shape as x\n    \"\"\"\n    y = 1 / (1 + np.exp(-x))\n    return y\n\nx, x_arr = 0, np.array([-5, -1, 1, 5])\nprint(f\"logistic({x}) = {logistic(x)}\")\nprint(f\"logistic({x_arr}) = {logistic(x_arr)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:11.765780Z","iopub.execute_input":"2022-07-25T09:39:11.766155Z","iopub.status.idle":"2022-07-25T09:39:11.773708Z","shell.execute_reply.started":"2022-07-25T09:39:11.766124Z","shell.execute_reply":"2022-07-25T09:39:11.772548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the logistic function\nplt.figure(figsize = (7.5, 6))\nx = np.linspace(-11, 11, 100)\nplt.plot(x, logistic(x), color = 'red')\nplt.xlabel(\"x\", fontsize = 14)\nplt.ylabel(\"g(x)\", fontsize = 14)\nplt.title(\"Standard logistic function\", fontsize = 14)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:11.775337Z","iopub.execute_input":"2022-07-25T09:39:11.775993Z","iopub.status.idle":"2022-07-25T09:39:12.181110Z","shell.execute_reply.started":"2022-07-25T09:39:11.775950Z","shell.execute_reply":"2022-07-25T09:39:12.179453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Log Loss","metadata":{}},{"cell_type":"markdown","source":"The `loss function`, which corresponds to the *true value* and *predicted value* of a single observation. The `cost function` can be thought of as *expected loss* or *average loss* over a group of observations. Contrary to linear regression, which employs `squared loss`, logistic regression makes use of the `log loss` function, given by\n\n$$ L(y, y') = -y \\log\\left(y'\\right) - \\left(1 - y\\right) \\log\\left(1 - y'\\right), $$\n\nwhere $y$ is the true value of a binary target (taking values $0$ or $1$) and $y'$ is the prediction, which can be thought of as the predicted probability of $y$ being $1$. Observe that the loss is $0$, when the true value and predicted value agree with each other, i.e. $L(0, 0) = L(1, 1) = 0$. On the other hand, the loss explodes towards infinity if the predicted value approaches $1$ when the true value is $0$, or it approaches $0$ when the true value is $1$. Mathematically, $\\lim_{t \\to 1-} L(0, t) = \\lim_{t \\to 0+} L(1, t) = \\infty$. In the next couple of code blocks, we construct the function to compute log loss and plot it for $y = 0$ and $y = 1$. Since the true values (labels) are always $0$ or $1$, we do not need to pay heed to the behaviour of the function $L$ for other values of $y$.","metadata":{}},{"cell_type":"code","source":"# Log loss\ndef log_loss(y, y_dash):\n    \"\"\"\n    Computes log loss for inputs true value (0 or 1) and predicted value (between 0 and 1)\n    Args:\n      y      (scalar): true value (0 or 1)\n      y_dash (scalar): predicted value (probability of y being 1)\n    Returns:\n      loss (float): nonnegative loss corresponding to y and y_dash\n    \"\"\"\n    loss = - (y * np.log(y_dash)) - ((1 - y) * np.log(1 - y_dash))\n    return loss\n\ny, y_dash = 0, 0.6\nprint(f\"log_loss({y}, {y_dash}) = {log_loss(y, y_dash)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.182662Z","iopub.execute_input":"2022-07-25T09:39:12.183021Z","iopub.status.idle":"2022-07-25T09:39:12.191362Z","shell.execute_reply.started":"2022-07-25T09:39:12.182974Z","shell.execute_reply":"2022-07-25T09:39:12.189885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Log loss for y = 0 and y = 1\nfig, ax = plt.subplots(1, 2, figsize = (15, 6), sharex = True, sharey = True)\ny_dash = np.linspace(0.0001, 0.9999, 100)\nax[0].plot(y_dash, log_loss(0, y_dash), color = 'red')\nax[0].set_title(\"y = 0\", fontsize = 14)\nax[0].set_xlabel(\"y_dash\", fontsize = 14)\nax[0].set_ylabel(\"log loss\", fontsize = 14)\nax[1].plot(y_dash, log_loss(1, y_dash), color = 'red')\nax[1].set_title(\"y = 1\", fontsize = 14)\nax[1].set_xlabel(\"y_dash\", fontsize = 14)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.195419Z","iopub.execute_input":"2022-07-25T09:39:12.196430Z","iopub.status.idle":"2022-07-25T09:39:12.603381Z","shell.execute_reply.started":"2022-07-25T09:39:12.196393Z","shell.execute_reply":"2022-07-25T09:39:12.602542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The plots sync with the intuition that loss should be minimum when the predicted value (probability) matches the true value $(0$ or $1)$, and should increase as the two values drift apart.","metadata":{}},{"cell_type":"markdown","source":"# Cost Function","metadata":{}},{"cell_type":"markdown","source":"Let $\\mathbf{y} = (y_1, y_2, \\cdots, y_n)$ be the true values $(0$ or $1)$ and $\\mathbf{y'} = (y_1', y_2', \\cdots, y_n')$ be the corresponding predictions (probabilities). Then, the *cost function* is given by the average loss:\n\n$$ C(\\mathbf{y}, \\mathbf{y'}) = \\frac{1}{m}\\sum_{i = 1}^m L(y_i, y_i'). $$\n\nWe construct the function to compute cost in the following two code blocks (the first one using for loop, the second one using vectorization). An important structural distinction from the log loss function is that here the arguments `y` and `y_dash` are vectors, not scalars.","metadata":{}},{"cell_type":"code","source":"# Cost function - using for loop\ndef cost_func(y, y_dash):\n    \"\"\"\n    Computes log loss for inputs true value (0 or 1) and predicted value (between 0 and 1)\n    Args:\n      y      (array_like, shape (m,)): array of true values (0 or 1)\n      y_dash (array_like, shape (m,)): array of predicted values (probability of y being 1)\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash\n    \"\"\"\n    assert len(y) == len(y_dash), \"Length of true values and length of predicted values do not match\"\n    m = len(y)\n    cost = 0\n    for i in range(m):\n        cost += log_loss(y[i], y_dash[i])\n    cost = cost / m\n    return cost\n\ny, y_dash = np.array([0, 1, 0]), np.array([0.4, 0.6, 0.25])\nprint(f\"cost_func({y}, {y_dash}) = {cost_func(y, y_dash)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.604356Z","iopub.execute_input":"2022-07-25T09:39:12.607339Z","iopub.status.idle":"2022-07-25T09:39:12.617060Z","shell.execute_reply.started":"2022-07-25T09:39:12.607290Z","shell.execute_reply":"2022-07-25T09:39:12.615775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cost function - using vectorization\ndef cost_func_vec(y, y_dash):\n    \"\"\"\n    Computes log loss for inputs true value (0 or 1) and predicted value (between 0 and 1)\n    Args:\n      y      (array_like, shape (m,)): array of true values (0 or 1)\n      y_dash (array_like, shape (m,)): array of predicted values (probability of y being 1)\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash\n    \"\"\"\n    assert len(y) == len(y_dash), \"Length of true values and length of predicted values do not match\"\n    m = len(y)\n    loss_vec = np.array([log_loss(y[i], y_dash[i]) for i in range(m)])\n    cost = np.dot(loss_vec, np.ones(m)) / m\n    return cost\n\ny, y_dash = np.array([0, 1, 0]), np.array([0.4, 0.6, 0.25])\nprint(f\"cost_func_vec({y}, {y_dash}) = {cost_func(y, y_dash)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.618492Z","iopub.execute_input":"2022-07-25T09:39:12.619303Z","iopub.status.idle":"2022-07-25T09:39:12.633030Z","shell.execute_reply.started":"2022-07-25T09:39:12.619261Z","shell.execute_reply":"2022-07-25T09:39:12.631867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us assume that we want to predict $y$ based on $n$ features. In this setup, a logistic regression model is characterized by $n+1$ parameters:\n\n- weight parameters $\\mathbf{w} = (w_1, w_2, \\cdots, w_n)$\n- bias parameter $b$\n\nNote that, the [dot product](https://en.wikipedia.org/wiki/Dot_product) of two vectors $\\mathbf{a} = (a_1, a_2, \\cdots, a_n)$ and $\\mathbf{b} = (b_1, b_2, \\cdots, b_n)$ is given by $\\mathbf{a} \\cdot \\mathbf{b} = \\sum_{i=1}^n a_ib_i$. It is a scalar value and evidently, $\\mathbf{a} \\cdot \\mathbf{b} = \\mathbf{b} \\cdot \\mathbf{a}$. Given the realized values of $n$ features $\\mathbf{x} = (x_1, x_2, \\cdots, x_n)$, the model feeds $\\mathbf{x} \\cdot \\mathbf{w} + b$ to the logistic function $g$, and projects the output as the predicted probability of $y = 1$. Concretely, we have\n\n$$ y' = g\\left(\\mathbf{x} \\cdot \\mathbf{w} + b\\right) = \\frac{1}{1 + e^{-\\left(\\mathbf{x} \\cdot \\mathbf{w} + b\\right)}}. \\tag{1} $$\n\nLet us consider the situation of $m$ observations, with the $i$th observation having feature values $\\mathbf{x_i} = (x_{i,1}, x_{i,2}, \\cdots, x_{i,n})$, true target value $y_i$ and predicted probabilities $y_i' = g\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)$. Stacking them up, we obtain the feature matrix $\\mathbf{X}$, target vector $\\mathbf{y}$ and the vector of predicted probability $\\mathbf{y'}$, as follows:\n\n$$ \\mathbf{X} = \\begin{pmatrix}\n\\mathbf{x_1} \\newline\n\\mathbf{x_2} \\newline\n\\vdots \\newline\n\\mathbf{x_m}\n\\end{pmatrix} = \\begin{pmatrix}\nx_{1,1} & x_{1,2} & \\cdots & x_{1,n} \\newline\nx_{2,1} & x_{2,2} & \\cdots & x_{2,n} \\newline\n\\vdots & \\vdots & \\ddots & \\vdots \\newline\nx_{m,1} & x_{m,2} & \\cdots & x_{m,n}\n\\end{pmatrix},\\;\\;\\;\\;\n\\mathbf{y} = \\begin{pmatrix}\ny_1 \\newline\ny_2 \\newline\n\\vdots \\newline\ny_m\n\\end{pmatrix},\\;\\;\\;\\;\n\\mathbf{y'} = \\begin{pmatrix}\ng\\left(\\mathbf{x_1} \\cdot \\mathbf{w} + b\\right) \\newline\ng\\left(\\mathbf{x_2} \\cdot \\mathbf{w} + b\\right) \\newline\n\\vdots \\newline\ng\\left(\\mathbf{x_n} \\cdot \\mathbf{w} + b\\right)\n\\end{pmatrix}. \\tag{2} $$\n\nNow, we are in a position to rewrite the cost function in terms of model parameters:\n\n$$ J\\left(\\mathbf{w}, b\\right) := C\\left(\\mathbf{y}, \\mathbf{y'} \\,\\vert\\, \\mathbf{X}, \\mathbf{w}, b \\right) = \\frac{1}{m}\\sum_{i = 1}^m L\\left(y_i, \\frac{1}{1 + e^{-\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)}}\\right) = \\frac{1}{m}\\sum_{i = 1}^m \\left[ -y_i \\log\\left(\\frac{1}{1 + e^{-\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)}}\\right) - \\left(1 - y_i\\right) \\log\\left(1 - \\frac{1}{1 + e^{-\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)}}\\right) \\right]. \\tag{3} $$\n\nNote that $J\\left(\\mathbf{0}, 0\\right) = \\log{2}$ for every input data $\\mathbf{X}$ (features) and $\\mathbf{y}$ (target). We construct a function to compute the cost, given data and model parameters, in the general setup of $n$ features (first using for loop, then using vectorization). Then, we visualize the cost function for a simplified setup, consisting of a single feature $x$. In this setup, the model involves two parameters only, the weight parameter $w$ and the bias parameter $b$.","metadata":{}},{"cell_type":"code","source":"# Function to compute cost function in terms of model parameters - using for loops\ndef cost_logreg(X, y, w, b):\n    \"\"\"\n    Computes the cost function, given data and model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values (0 or 1) of target\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash \n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    z = []\n    for i in range(m):\n        s = 0\n        for j in range(n):\n            s += X[i, j] * w[j]\n        z.append(s + b)\n    z = np.array(z)\n    y_dash = logistic(z)\n    cost = cost_func(y, y_dash)\n    return cost\n\nX, y, w, b = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1\nprint(f\"cost_logreg(X = {X}, y = {y}, w = {w}, b = {b}) = {cost_logreg(X, y, w, b)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.637865Z","iopub.execute_input":"2022-07-25T09:39:12.639060Z","iopub.status.idle":"2022-07-25T09:39:12.656507Z","shell.execute_reply.started":"2022-07-25T09:39:12.638978Z","shell.execute_reply":"2022-07-25T09:39:12.655399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The code for logistic function $g$ is constructed in such a way that, if applied on an array, it acts separately on each component and returns an array. Using this, it follows from $(2)$ that\n\n$$ \\mathbf{y'} = g\\left(\\mathbf{X} \\mathbf{w} + b \\mathbf{1} \\right), $$\n\nwhere $\\mathbf{1}$ has the same dimension as that of $\\mathbf{X} \\mathbf{w}$. We use [matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication) in computing $\\mathbf{X} \\mathbf{w}$, [scalar multiplication of a vector](https://en.wikipedia.org/wiki/Scalar_multiplication) in computing $b \\mathbf{1}$ and add them up using [vector addition](https://en.wikipedia.org/wiki/Euclidean_vector#Addition_and_subtraction). This representation leads to a much faster vectorized implementation of computing the cost function in terms of model parameters.","metadata":{}},{"cell_type":"code","source":"# Function to compute cost function in terms of model parameters - using vectorization\ndef cost_logreg_vec(X, y, w, b):\n    \"\"\"\n    Computes the cost function, given data and model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values of target (0 or 1)\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash \n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    z = np.matmul(X, w) + (b * np.ones(m))\n    y_dash = logistic(z)\n    cost = cost_func_vec(y, y_dash)\n    return cost\n\nX, y, w, b = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1\nprint(f\"cost_logreg_vec(X = {X}, y = {y}, w = {w}, b = {b}) = {cost_logreg(X, y, w, b)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.657940Z","iopub.execute_input":"2022-07-25T09:39:12.658880Z","iopub.status.idle":"2022-07-25T09:39:12.674052Z","shell.execute_reply.started":"2022-07-25T09:39:12.658834Z","shell.execute_reply":"2022-07-25T09:39:12.672869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the cost function against model parameters\nfrom mpl_toolkits.mplot3d import Axes3D\nw, b = np.meshgrid(np.linspace(-0.1, 0.1, 21), np.linspace(-1, 1, 21))\nX = np.array([1.56, 0.76 , 0.08, 9.71, 4.65, 4.35, 7.34, 0.91, 9.82, 9.05]).reshape((10, 1))\ny = np.array([0, 1, 0, 0, 0, 0, 1, 1, 1, 0])\ncost = np.array([[cost_logreg_vec(X, y, np.array([w0]), b0) for b0 in b[:, 0]] for w0 in w[0]])\nfig = plt.figure(figsize = (10, 10))\nax = plt.axes(projection = '3d')\nax.plot_surface(w, b, cost)\nax.set(xlabel = \"w\", ylabel = \"b\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:12.675223Z","iopub.execute_input":"2022-07-25T09:39:12.676392Z","iopub.status.idle":"2022-07-25T09:39:13.088563Z","shell.execute_reply.started":"2022-07-25T09:39:12.676349Z","shell.execute_reply":"2022-07-25T09:39:13.087543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The prediction $y'$ in $(1)$ can be converted to a decision by using a threshold value. For instance, suppose we take the threshold to be $0.5$. Then, we may classify the observation to class $1$ if $y' \\geq 0.5$, and to class $0$ otherwise. The problem, however, is that we do not know the model parameters $\\mathbf{w}$ and $b$, and hence cannot compute $y'$ directly. First, we have to fit the model by finding the best fitting parameters. Observe that, given the training data, $C\\left(\\mathbf{y}, \\mathbf{y'}\\right)$ depends on $\\mathbf{w}$ and $b$ only. A reasonable strategy, therefore, is:\n\n$$ \\text{To minimise } J\\left(\\mathbf{w}, b\\right), \\text{ with respect to } \\mathbf{w} \\text{ and } b. $$\n\nWe shall employ the [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) [algorithm](https://en.wikipedia.org/wiki/Algorithm) to solve this [optimization](https://en.wikipedia.org/wiki/Mathematical_optimization) problem.","metadata":{}},{"cell_type":"markdown","source":"# Gradient Descent","metadata":{}},{"cell_type":"markdown","source":"**What is it?** The gradient descent algorithm is a [first order](https://en.wikipedia.org/wiki/Category:First_order_methods) [iterative](https://en.wikipedia.org/wiki/Iterative_method) optimization algorithm for finding a [local minimum](https://en.wikipedia.org/wiki/Maxima_and_minima) of a differentiable function. The idea is to take repeated steps in the opposite direction of the [gradient](https://en.wikipedia.org/wiki/Gradient) (or approximate gradient) of the function at the current point, because this is the direction of steepest descent. Conversely, stepping in the direction of the gradient will lead to a local maximum of that function. This procedure is then known as *gradient ascent*.\n\n**History.** Gradient descent is generally attributed to [Augustin-Louis Cauchy](https://en.wikipedia.org/wiki/Augustin-Louis_Cauchy), who first suggested it in 1847. [Jacques Hadamard](https://en.wikipedia.org/wiki/Jacques_Hadamard) independently proposed a similar method in 1907. Its convergence properties for non-linear optimization problems were first studied by [Haskell Curry](https://en.wikipedia.org/wiki/Haskell_Curry) in 1944, with the method becoming increasingly well-studied and used in the following decades.\n\n**The algorithm.** In the context of minimising the cost function $J$, with respect to the model parameters $\\mathbf{w}$ and $b$, the gradient descent algorithm is given by:\n\n$$ \\begin{align*}\n& \\text{repeat until convergence:}\\; \\{ \\newline\n& w_j := w_j -  \\alpha \\frac{\\partial J(\\mathbf{w},b)}{\\partial w_j},\\; \\text{ for } j = 1, 2, \\ldots, n; \\newline\n& b := b -  \\alpha \\frac{\\partial J(\\mathbf{w}, b)}{\\partial b}. \\tag{4} \\newline\n& \\}\n\\end{align*} $$\n\nwhere $\\alpha$ is the [learning rate](https://en.wikipedia.org/wiki/Learning_rate), and the parameters $\\mathbf{w} = (w_1, w_2, \\cdots, w_n)$ and $b$ are updated simultaniously in each iteration.\n\n**Computing gradient.** Before we can implement the gradient descent algorithm, we need to compute the gradients first! From $(3)$, we can compute the partial derivatives of $J$ with respect to $w_j$ and $b$ as follows:\n\n$$ \\begin{align*}\n& \\frac{\\partial J(\\mathbf{w},b)}{\\partial w_j} = \\frac{1}{m} \\sum\\limits_{i = 1}^m \\left(\\frac{1}{1 + e^{-\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)}} - y_i\\right) x_{i,j},\\;\\; \\text{ for } j = 1, 2, \\ldots, n; \\newline\n& \\frac{\\partial J(\\mathbf{w}, b)}{\\partial b} = \\frac{1}{m} \\sum\\limits_{i = 1}^m \\left(\\frac{1}{1 + e^{-\\left(\\mathbf{x_i} \\cdot \\mathbf{w} + b\\right)}} - y_i\\right).\n\\end{align*} $$\n\nThe next code blocks construct a function to compute these gradients, first using for loops, then using vectorization.","metadata":{}},{"cell_type":"code","source":"# Function to compute gradients of the cost function with respect to model parameters - using for loops\ndef grad_logreg(X, y, w, b):\n    \"\"\"\n    Computes gradients of the cost function with respect to model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values of target (0 or 1)\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n    Returns:\n      grad_w (array_like, shape (n,)): gradients of the cost function with respect to the weight parameters\n      grad_b (float)                 : gradient of the cost function with respect to the bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    grad_w, grad_b = np.zeros(n), 0\n    for i in range(m):\n        s = 0\n        for j in range(n):\n            s += X[i, j] * w[j]\n        y_dash_i = logistic(s + b)\n        for j in range(n):\n            grad_w[j] += (y_dash_i  - y[i]) * X[i,j]\n        grad_b += y_dash_i  - y[i]\n    grad_w, grad_b = grad_w / m, grad_b / m\n    return grad_w, grad_b\n\nX, y, w, b = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1\nprint(f\"grad_logreg(X = {X}, y = {y}, w = {w}, b = {b}) = {grad_logreg(X, y, w, b)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:13.089815Z","iopub.execute_input":"2022-07-25T09:39:13.090235Z","iopub.status.idle":"2022-07-25T09:39:13.104792Z","shell.execute_reply.started":"2022-07-25T09:39:13.090189Z","shell.execute_reply":"2022-07-25T09:39:13.103571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to compute gradients of the cost function with respect to model parameters - using vectorization\ndef grad_logreg_vec(X, y, w, b): \n    \"\"\"\n    Computes gradients of the cost function with respect to model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values of target (0 or 1)\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n    Returns:\n      grad_w (array_like, shape (n,)): gradients of the cost function with respect to the weight parameters\n      grad_b (float)                 : gradient of the cost function with respect to the bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    y_dash = logistic(np.matmul(X, w) + b * np.ones(m))\n    grad_w = np.matmul(y_dash - y, X) / m\n    grad_b = np.dot(y_dash - y, np.ones(m)) / m\n    return grad_w, grad_b\n\nX, y, w, b = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1\nprint(f\"grad_logreg_vec(X = {X}, y = {y}, w = {w}, b = {b}) = {grad_logreg_vec(X, y, w, b)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:13.106384Z","iopub.execute_input":"2022-07-25T09:39:13.107319Z","iopub.status.idle":"2022-07-25T09:39:13.126350Z","shell.execute_reply.started":"2022-07-25T09:39:13.107272Z","shell.execute_reply":"2022-07-25T09:39:13.124607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Armed with the required functions, we can now implement the gradient descent algorithm, given in $(4)$. ","metadata":{}},{"cell_type":"code","source":"# Gradient descent algorithm for logistic regression\ndef grad_desc(X, y, w, b, alpha, n_iter, show_cost = True, show_params = False): \n    \"\"\"\n    Implements batch gradient descent algorithm to learn and update model parameters\n    with prespecified number of interations and learning rate\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): true values of target (0 or 1)\n      w (array_like, shape (n,)): initial value of weight parameters\n      b (scalar)                : initial value of bias parameter\n      cost_func                 : function to compute cost\n      grad_func                 : function to compute gradients of cost with respect to model parameters\n      alpha (float)             : learning rate\n      n_iter (int)              : number of iterations\n    Returns:\n      w (array_like, shape (n,)): updated values of weight parameters\n      b (scalar)                : updated value of bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    cost_history, params_history = [], []\n    for i, j in itertools.product(range(n_iter), range(1)):\n        grad_w, grad_b = grad_logreg_vec(X, y, w, b)   \n        w += - alpha * grad_w\n        b += - alpha * grad_b\n        cost =  cost_logreg_vec(X, y, w, b)\n        cost_history.append(cost)\n        params_history.append([w, b])\n        if show_cost == True and show_params == False and (i % math.ceil(n_iter / 10) == 0 or i == n_iter - 1):\n            print(f\"Iteration {i:6}:    Cost  {float(cost_history[i]):3.4f}\")\n        if show_cost == True and show_params == True and (i % math.ceil(n_iter / 10) == 0 or i == n_iter - 1):\n            print(f\"Iteration {i:6}:    Cost  {float(cost_history[i]):3.4f},    Params  {params_history[i]}\")\n    return w, b, cost_history, params_history\n\nX, y, w, b, alpha, n_iter = np.array([[0.1, 0.2], [-0.1, 0.1]]), np.array([1, 0]), np.array([0., 0.]), 0., 0.1, 100000\nw_out, b_out, cost_history, params_history = grad_desc(X, y, w, b, alpha, n_iter, show_cost = True, show_params = True)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:13.127716Z","iopub.execute_input":"2022-07-25T09:39:13.128326Z","iopub.status.idle":"2022-07-25T09:39:20.676562Z","shell.execute_reply.started":"2022-07-25T09:39:13.128295Z","shell.execute_reply":"2022-07-25T09:39:20.675403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting cost over iteration\nplt.figure(figsize = (9, 6))\nplt.plot(cost_history)\nplt.xlabel(\"Iteration\", fontsize = 14)\nplt.ylabel(\"Cost\", fontsize = 14)\nplt.title(\"Cost vs Iteration\", fontsize = 14)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:20.678209Z","iopub.execute_input":"2022-07-25T09:39:20.678506Z","iopub.status.idle":"2022-07-25T09:39:20.978806Z","shell.execute_reply.started":"2022-07-25T09:39:20.678480Z","shell.execute_reply":"2022-07-25T09:39:20.977728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing","metadata":{}},{"cell_type":"markdown","source":"The columns `EventId` and `Weight` will not be used in this notebook. Thus we drop these columns. Even though there are no `NaN` values in the dataset, we observe that several columns unnaturally contain the value $-999$ for many observations. We suspect that these are missing/corrupted values that have been replaced with $-999$ and we convert them to `np.nan`. We shall impute these values after *train-test split*. Furthermore, we encode the `Label` column as follows: $b \\mapsto 0$ and $s \\mapsto 1$.","metadata":{}},{"cell_type":"code","source":"# Dropping unnecessary columns\ndata.drop(['EventId', 'Weight'], axis = 1, inplace = True)\n\n# Replacing -999 with nan\ndata.replace(to_replace = -999, value = np.nan, inplace = True)\n\n# Encoding the 'Label' column\nlabel_dict = {'b': 0, 's': 1}\ndata.replace({'Label': label_dict}, inplace = True)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:20.980099Z","iopub.execute_input":"2022-07-25T09:39:20.980407Z","iopub.status.idle":"2022-07-25T09:39:21.127482Z","shell.execute_reply.started":"2022-07-25T09:39:20.980380Z","shell.execute_reply":"2022-07-25T09:39:21.126590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we split `data` into two parts:\n\n- `data_train` : The portion of data that we use to train the models\n- `data_test` : The portion of data that we use to test or evaluate the models","metadata":{}},{"cell_type":"code","source":"# Train-test split\ndata_train, data_test = train_test_split(data, test_size = 0.2, random_state = 40)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.128774Z","iopub.execute_input":"2022-07-25T09:39:21.129289Z","iopub.status.idle":"2022-07-25T09:39:21.214120Z","shell.execute_reply.started":"2022-07-25T09:39:21.129256Z","shell.execute_reply":"2022-07-25T09:39:21.213000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we identify $11$ columns with missing values. Among them, $10$ columns have more than $30\\%$ data missing, and hence we shall discard them. The column `DER_mass_MMC` has about $15.25\\%$ data missing. We shall impute the missing values in this column by the median of the rest of the values in the column.","metadata":{}},{"cell_type":"code","source":"# Columns with missing values with respective proportions\n(data.isna().sum()[data.isna().sum() > 0] / len(data)).sort_values(ascending = False)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.215511Z","iopub.execute_input":"2022-07-25T09:39:21.215821Z","iopub.status.idle":"2022-07-25T09:39:21.256629Z","shell.execute_reply.started":"2022-07-25T09:39:21.215793Z","shell.execute_reply":"2022-07-25T09:39:21.255817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Discarding columns with more than 30% missing data\ncols_missing_drop = [\n    'DER_deltaeta_jet_jet',\n    'DER_mass_jet_jet',\n    'DER_prodeta_jet_jet',\n    'DER_lep_eta_centrality',\n    'PRI_jet_subleading_pt',\n    'PRI_jet_subleading_eta',\n    'PRI_jet_subleading_phi',\n    'PRI_jet_leading_pt',\n    'PRI_jet_leading_eta',\n    'PRI_jet_leading_phi'\n]\ndata_train.drop(cols_missing_drop, axis = 1, inplace = True)\ndata_test.drop(cols_missing_drop, axis = 1, inplace = True)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.257910Z","iopub.execute_input":"2022-07-25T09:39:21.258426Z","iopub.status.idle":"2022-07-25T09:39:21.282159Z","shell.execute_reply.started":"2022-07-25T09:39:21.258391Z","shell.execute_reply":"2022-07-25T09:39:21.281170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Median imputation\ndata_train['DER_mass_MMC'].fillna(data_train['DER_mass_MMC'].median(), inplace = True)\ndata_test['DER_mass_MMC'].fillna(data_test['DER_mass_MMC'].median(), inplace = True)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.283185Z","iopub.execute_input":"2022-07-25T09:39:21.284147Z","iopub.status.idle":"2022-07-25T09:39:21.300794Z","shell.execute_reply.started":"2022-07-25T09:39:21.284106Z","shell.execute_reply":"2022-07-25T09:39:21.299793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we shall split the data into features and target.","metadata":{}},{"cell_type":"code","source":"# Features-target split\nX_train, y_train = data_train.drop('Label', axis = 1), data_train['Label']\nX_test, y_test = data_test.drop('Label', axis = 1), data_test['Label']","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.302192Z","iopub.execute_input":"2022-07-25T09:39:21.302759Z","iopub.status.idle":"2022-07-25T09:39:21.322620Z","shell.execute_reply.started":"2022-07-25T09:39:21.302728Z","shell.execute_reply":"2022-07-25T09:39:21.321327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We normalize the data, so that each column has values in the similar scale.","metadata":{}},{"cell_type":"code","source":"# Min-Max normalization\nfor col in X_train.columns:\n    if (X_train[col].dtypes == 'int64' or X_train[col].dtypes == 'float64') and X_train[col].nunique() > 1:\n        X_train[col] = (X_train[col] - X_train[col].min()) / (X_train[col].max() - X_train[col].min())\nfor col in X_test.columns:\n    if (X_test[col].dtypes == 'int64' or X_test[col].dtypes == 'float64') and X_test[col].nunique() > 1:\n        X_test[col] = (X_test[col] - X_test[col].min()) / (X_test[col].max() - X_test[col].min())","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.323998Z","iopub.execute_input":"2022-07-25T09:39:21.324553Z","iopub.status.idle":"2022-07-25T09:39:21.546186Z","shell.execute_reply.started":"2022-07-25T09:39:21.324509Z","shell.execute_reply":"2022-07-25T09:39:21.545063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Fitting","metadata":{}},{"cell_type":"markdown","source":"We fix the initial values of the parameters, based on running the algorithm several times and noting down the final parameter values. It gives us a better *starting point* and helps to achieve a better performance in a limited number of iterations.","metadata":{}},{"cell_type":"code","source":"# Initial values of the model parameters\nw_init = np.array([-5, -15, -10, 9, 4, -6, 3, -10, 1, 14, 0, 0, 15, 0, 0, 7, 0, -3, 1, -8]).astype(float)\nb_init = -1.","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.547868Z","iopub.execute_input":"2022-07-25T09:39:21.548483Z","iopub.status.idle":"2022-07-25T09:39:21.554707Z","shell.execute_reply.started":"2022-07-25T09:39:21.548439Z","shell.execute_reply":"2022-07-25T09:39:21.553397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Learning model parameters using gradient descent algorithm\nw_out, b_out, cost_history, params_history = grad_desc(X_train.to_numpy(),\n                                                       y_train.to_numpy(),\n                                                       w = w_init, # np.zeros(X_train.shape[1]),\n                                                       b = b_init, # 0,\n                                                       alpha = 0.1,\n                                                       n_iter = 2000)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T09:39:21.556388Z","iopub.execute_input":"2022-07-25T09:39:21.556983Z","iopub.status.idle":"2022-07-25T10:34:19.603570Z","shell.execute_reply.started":"2022-07-25T09:39:21.556949Z","shell.execute_reply":"2022-07-25T10:34:19.602126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting cost over iteration\nplt.figure(figsize = (9, 6))\nplt.plot(cost_history)\nplt.xlabel(\"Iteration\", fontsize = 14)\nplt.ylabel(\"Cost\", fontsize = 14)\nplt.title(\"Cost vs Iteration\", fontsize = 14)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:19.605410Z","iopub.execute_input":"2022-07-25T10:34:19.606061Z","iopub.status.idle":"2022-07-25T10:34:19.932722Z","shell.execute_reply.started":"2022-07-25T10:34:19.605986Z","shell.execute_reply":"2022-07-25T10:34:19.931574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Final parameter values\nparams_history[-1]","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:19.934120Z","iopub.execute_input":"2022-07-25T10:34:19.934464Z","iopub.status.idle":"2022-07-25T10:34:19.942065Z","shell.execute_reply.started":"2022-07-25T10:34:19.934432Z","shell.execute_reply":"2022-07-25T10:34:19.940933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction and Evaluation","metadata":{}},{"cell_type":"markdown","source":"First, we construct some functions to compute and display the *confusion matrix*, and to compute *accuracy*, given the true labels and the predicted labels of the target.","metadata":{}},{"cell_type":"code","source":"# Function to compute confusion matrix\ndef conf_mat(y_test, y_pred):\n    \"\"\"\n    Computes confusion matrix\n    Args:\n      y_test (array_like): true binary (0 or 1) labels\n      y_pred (array_like): predicted binary (0 or 1) labels\n    Returns:\n      confusion_mat (array): A 2D array representing a 2x2 confusion matrix\n    \"\"\"\n    y_test, y_pred = list(y_test), list(y_pred)\n    count, labels, confusion_mat = len(y_test), [0, 1], np.zeros(shape = (2, 2), dtype = int)\n    for i in range(2):\n        for j in range(2):\n            confusion_mat[i][j] = len([k for k in range(count) if y_test[k] == labels[i] and y_pred[k] == labels[j]])\n    return confusion_mat","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:19.948313Z","iopub.execute_input":"2022-07-25T10:34:19.948721Z","iopub.status.idle":"2022-07-25T10:34:19.958504Z","shell.execute_reply.started":"2022-07-25T10:34:19.948687Z","shell.execute_reply":"2022-07-25T10:34:19.957451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to print confusion matrix\ndef conf_mat_heatmap(y_test, y_pred):\n    \"\"\"\n    Prints confusion matrix\n    Args:\n      y_test (array_like): true binary (0 or 1) labels\n      y_pred (array_like): predicted binary (0 or 1) labels\n    Returns:\n      Nothing, prints a heatmap representing a 2x2 confusion matrix\n    \"\"\"\n    confusion_mat = conf_mat(y_test, y_pred)\n    labels, confusion_mat_df = [0, 1], pd.DataFrame(confusion_mat, range(2), range(2))\n    plt.figure(figsize = (6, 4.75))\n    sns.heatmap(confusion_mat_df, annot = True, annot_kws = {\"size\": 16}, fmt = 'd')\n    plt.xticks([0.5, 1.5], labels, rotation = 'horizontal')\n    plt.yticks([0.5, 1.5], labels, rotation = 'horizontal')\n    plt.xlabel(\"Predicted label\", fontsize = 14)\n    plt.ylabel(\"True label\", fontsize = 14)\n    plt.title(\"Confusion Matrix\", fontsize = 14)\n    plt.grid(False)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:19.960336Z","iopub.execute_input":"2022-07-25T10:34:19.960739Z","iopub.status.idle":"2022-07-25T10:34:19.972300Z","shell.execute_reply.started":"2022-07-25T10:34:19.960709Z","shell.execute_reply":"2022-07-25T10:34:19.971363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to compute accuracy\ndef accuracy(y_test, y_pred):\n    \"\"\"\n    Computes accuracy, given true and predicted binary (0 or 1) labels\n    Args:\n      y_test (array_like): true binary (0 or 1) labels\n      y_pred (array_like): predicted binary (0 or 1) labels\n    Returns:\n      acc (float): accuracy obtained from y_test and y_pred\n    \"\"\"\n    confusion_mat = conf_mat(y_test, y_pred)\n    num = confusion_mat[0, 0] + confusion_mat[1, 1] # Number of correct predictions\n    denom = num + confusion_mat[0, 1] + confusion_mat[1, 0] # Number of total predictions\n    acc = num / denom\n    return acc","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:19.973620Z","iopub.execute_input":"2022-07-25T10:34:19.974547Z","iopub.status.idle":"2022-07-25T10:34:20.009625Z","shell.execute_reply.started":"2022-07-25T10:34:19.974509Z","shell.execute_reply":"2022-07-25T10:34:20.008788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction and evaluation on the training set and the test set\ny_train_prob = logistic(np.matmul(X_train.to_numpy(), w_out) + (b_out * np.ones(X_train.shape[0])))\ny_test_prob = logistic(np.matmul(X_test.to_numpy(), w_out) + (b_out * np.ones(X_test.shape[0])))\ny_train_pred, y_test_pred = (y_train_prob > 0.5).astype(int), (y_test_prob > 0.5).astype(int)\nprint(pd.Series({\"Training accuracy\": accuracy(y_train, y_train_pred),\n                 \"Test accuracy\": accuracy(y_test, y_test_pred)}).to_string())","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.010909Z","iopub.execute_input":"2022-07-25T10:34:20.011908Z","iopub.status.idle":"2022-07-25T10:34:20.331414Z","shell.execute_reply.started":"2022-07-25T10:34:20.011872Z","shell.execute_reply":"2022-07-25T10:34:20.330319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Confusion matrix for predictions on the test set\nconf_mat_heatmap(y_test, y_test_pred)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.332641Z","iopub.execute_input":"2022-07-25T10:34:20.332945Z","iopub.status.idle":"2022-07-25T10:34:20.603868Z","shell.execute_reply.started":"2022-07-25T10:34:20.332916Z","shell.execute_reply":"2022-07-25T10:34:20.602031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Regularization","metadata":{}},{"cell_type":"markdown","source":"*Regularization* provides a tool to deal with overfitting. Essentially it puts a restriction the parameter values by adding a term in the cost function. This way, it prevents the model from fitting *too well* on the training set and failing to generalize on the test set. The following is one way to do this:\n\n$$ \\text{new } J(\\mathbf{w}, b) := \\text{old } J(\\mathbf{w}, b) + \\frac{\\lambda}{2m} \\sum_{j=1}^n w_j^2, \\tag{5} $$\n\nwhere $\\lambda \\geq 0$ is called the regularization parameter. Evidently, if $\\lambda = 0$, we get back the previous setup. Thus, the new cost function, parameterized by $\\lambda$, gives a generalization of the setup described in the previous sections. The setup in $(5)$ is referred to as the $L_2$ regularization, as the regularization term essentially is $(\\lambda/2m)\\,\\left\\Vert \\mathbf{w} \\right\\Vert_2^2$, i.e. the [$l^2$ norm](https://mathworld.wolfram.com/L2-Norm.html) of $\\mathbf{w}$, multiplied by the term $(\\lambda/2m)$. In this section, we shall employ the $L_2$ regularization on the cost function, given in $(3)$, and check if it improves performance of the logistic regression model, derived thence. We rewrite the functions to compute cost, gradient and to implement gradient descent algorithm, incorporating regularization with the extra parameter $l$.","metadata":{}},{"cell_type":"code","source":"# Function to compute regularized cost function in terms of model parameters - using for loops\ndef cost_logreg_reg(X, y, w, b, l):\n    \"\"\"\n    Computes the cost function, given data and model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values (0 or 1) of target\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n      l (float)                 : regularization parameter\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash \n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    cost = cost_logreg(X, y, w, b)\n    for j in range(n):\n        cost += (l / (2 * m)) * (w[j]**2)\n    return cost\n\nX, y, w, b, l = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1, 1\nprint(f\"cost_logreg(X = {X}, y = {y}, w = {w}, b = {b}, l = {l}) = {cost_logreg_reg(X, y, w, b, l)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.605139Z","iopub.execute_input":"2022-07-25T10:34:20.605439Z","iopub.status.idle":"2022-07-25T10:34:20.614714Z","shell.execute_reply.started":"2022-07-25T10:34:20.605412Z","shell.execute_reply":"2022-07-25T10:34:20.613458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to compute regularized cost function in terms of model parameters - using vectorization\ndef cost_logreg_vec_reg(X, y, w, b, l):\n    \"\"\"\n    Computes the cost function, given data and model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values (0 or 1) of target\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n      l (float)                 : regularization parameter\n    Returns:\n      cost (float): nonnegative cost corresponding to y and y_dash \n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    cost = cost_logreg_vec(X, y, w, b)\n    cost += (l / (2 * m)) * np.dot(w, w)\n    return cost\n\nX, y, w, b, l = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1, 1\nprint(f\"cost_logreg_vec_reg(X = {X}, y = {y}, w = {w}, b = {b}, l = {l}) = {cost_logreg_vec_reg(X, y, w, b, l)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.616135Z","iopub.execute_input":"2022-07-25T10:34:20.616539Z","iopub.status.idle":"2022-07-25T10:34:20.632556Z","shell.execute_reply.started":"2022-07-25T10:34:20.616507Z","shell.execute_reply":"2022-07-25T10:34:20.631467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From $(5)$, we can compute the modified partial derivatives of $J$ with respect to $w_j$ and $b$ as follows:\n\n$$ \\begin{align*}\n& \\text{new } \\frac{\\partial J(\\mathbf{w},b)}{\\partial w_j} = \\text{old } \\frac{\\partial J(\\mathbf{w},b)}{\\partial w_j} + \\frac{\\lambda}{m} w_j,\\;\\; \\text{ for } j = 1, 2, \\ldots, n; \\newline\n& \\text{new } \\frac{\\partial J(\\mathbf{w}, b)}{\\partial b} = \\text{old } \\frac{\\partial J(\\mathbf{w}, b)}{\\partial b}.\n\\end{align*} $$\n\nWe construct a function to compute these gradients, first using for loops, then using vectorization.","metadata":{}},{"cell_type":"code","source":"# Function to compute gradients of the regularized cost function with respect to model parameters - using for loops\ndef grad_logreg_reg(X, y, w, b, l):\n    \"\"\"\n    Computes gradients of the cost function with respect to model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values of target (0 or 1)\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n      l (float)                 : regularization parameter\n    Returns:\n      grad_w (array_like, shape (n,)): gradients of the cost function with respect to the weight parameters\n      grad_b (float)                 : gradient of the cost function with respect to the bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    grad_w, grad_b = grad_logreg(X, y, w, b)\n    for j in range(n):\n        grad_w[j] += (l / m) * w[j]\n    return grad_w, grad_b\n\nX, y, w, b, l = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1, 1\nprint(f\"grad_logreg(X = {X}, y = {y}, w = {w}, b = {b}, l = {l}) = {grad_logreg_reg(X, y, w, b, l)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.634034Z","iopub.execute_input":"2022-07-25T10:34:20.634408Z","iopub.status.idle":"2022-07-25T10:34:20.655074Z","shell.execute_reply.started":"2022-07-25T10:34:20.634374Z","shell.execute_reply":"2022-07-25T10:34:20.653881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to compute gradients of the regularized cost function with respect to model parameters - using vectorization\ndef grad_logreg_vec_reg(X, y, w, b, l):\n    \"\"\"\n    Computes gradients of the cost function with respect to model parameters\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): array of true values of target (0 or 1)\n      w (array_like, shape (n,)): weight parameters of the model      \n      b (float)                 : bias parameter of the model\n      l (float)                 : regularization parameter\n    Returns:\n      grad_w (array_like, shape (n,)): gradients of the cost function with respect to the weight parameters\n      grad_b (float)                 : gradient of the cost function with respect to the bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    grad_w, grad_b = grad_logreg_vec(X, y, w, b)\n    grad_w += (l / m) * w\n    return grad_w, grad_b\n\nX, y, w, b, l = np.array([[10, 20], [-10, 10]]), np.array([1, 0]), np.array([0.5, 1.5]), 1, 1\nprint(f\"grad_logreg_vec_reg(X = {X}, y = {y}, w = {w}, b = {b}, l = {l}) = {grad_logreg_vec_reg(X, y, w, b, l)}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.656582Z","iopub.execute_input":"2022-07-25T10:34:20.657038Z","iopub.status.idle":"2022-07-25T10:34:20.676705Z","shell.execute_reply.started":"2022-07-25T10:34:20.656941Z","shell.execute_reply":"2022-07-25T10:34:20.675569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the regularised cost function and its gradients with respect to the model parameters, we construct a function to implement the gradient descent algorithm for logistic regression, incorporating regularization. Specifically, the algorithm learns the model parameters, with the goal of minimizing the regularised cost function in $(5)$, instead of the usual cost function in $(3)$.","metadata":{}},{"cell_type":"code","source":"# Gradient descent algorithm for logistic regression with regularization\ndef grad_desc_reg(X, y, w, b, l, alpha, n_iter, show_cost = True, show_params = False): \n    \"\"\"\n    Implements batch gradient descent algorithm to learn and update model parameters\n    with prespecified number of interations and learning rate\n    Args:\n      X (ndarray, shape (m,n))  : data on features, m observations with n features\n      y (array_like, shape (m,)): true values of target (0 or 1)\n      w (array_like, shape (n,)): initial value of weight parameters\n      b (scalar)                : initial value of bias parameter\n      l (float)                 : regularization parameter\n      alpha (float)             : learning rate\n      n_iter (int)              : number of iterations\n    Returns:\n      w (array_like, shape (n,)): updated values of weight parameters\n      b (scalar)                : updated value of bias parameter\n    \"\"\"\n    m, n = X.shape\n    assert len(y) == m, \"Number of feature observations and number of target observations do not match\"\n    assert len(w) == n, \"Number of features and number of weight parameters do not match\"\n    cost_history, params_history = [], []\n    for i, j in itertools.product(range(n_iter), range(1)):\n        grad_w, grad_b = grad_logreg_vec_reg(X, y, w, b, l)   \n        w += - alpha * grad_w\n        b += - alpha * grad_b\n        cost =  cost_logreg_vec_reg(X, y, w, b, l)\n        cost_history.append(cost)\n        params_history.append([w, b])\n        if show_cost == True and show_params == False and (i % math.ceil(n_iter / 10) == 0 or i == n_iter - 1):\n            print(f\"Iteration {i:6}:    Cost  {float(cost_history[i]):3.4f}\")\n        if show_cost == True and show_params == True and (i % math.ceil(n_iter / 10) == 0 or i == n_iter - 1):\n            print(f\"Iteration {i:6}:    Cost  {float(cost_history[i]):3.4f},    Params  {params_history[i]}\")\n    return w, b, cost_history, params_history","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.678420Z","iopub.execute_input":"2022-07-25T10:34:20.678698Z","iopub.status.idle":"2022-07-25T10:34:20.697950Z","shell.execute_reply.started":"2022-07-25T10:34:20.678675Z","shell.execute_reply":"2022-07-25T10:34:20.696890Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use the same initial values for the model parameters as in the unregularized implementation.","metadata":{}},{"cell_type":"code","source":"# Initial values of the model parameters\nw_init = np.array([-5, -15, -10, 9, 4, -6, 3, -10, 1, 14, 0, 0, 15, 0, 0, 7, 0, -3, 1, -8]).astype(float)\nb_init = -1.","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.699377Z","iopub.execute_input":"2022-07-25T10:34:20.700358Z","iopub.status.idle":"2022-07-25T10:34:20.718222Z","shell.execute_reply.started":"2022-07-25T10:34:20.700312Z","shell.execute_reply":"2022-07-25T10:34:20.717291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Learning model parameters using gradient descent algorithm\nw_out_reg, b_out_reg, cost_history_reg, params_history_reg = grad_desc_reg(X_train.to_numpy(),\n                                                                           y_train.to_numpy(),\n                                                                           w = w_init, # np.zeros(X_train.shape[1]),\n                                                                           b = b_init, # 0,\n                                                                           l = 1.,\n                                                                           alpha = 0.1,\n                                                                           n_iter = 2000)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T10:34:20.719255Z","iopub.execute_input":"2022-07-25T10:34:20.719955Z","iopub.status.idle":"2022-07-25T11:28:56.056927Z","shell.execute_reply.started":"2022-07-25T10:34:20.719921Z","shell.execute_reply":"2022-07-25T11:28:56.055562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting cost over iteration\nplt.figure(figsize = (9, 6))\nplt.plot(cost_history_reg)\nplt.xlabel(\"Iteration\", fontsize = 14)\nplt.ylabel(\"Cost\", fontsize = 14)\nplt.title(\"Cost vs Iteration\", fontsize = 14)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-25T11:28:56.058653Z","iopub.execute_input":"2022-07-25T11:28:56.059256Z","iopub.status.idle":"2022-07-25T11:28:56.378617Z","shell.execute_reply.started":"2022-07-25T11:28:56.059209Z","shell.execute_reply":"2022-07-25T11:28:56.377314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Final parameter values\nparams_history_reg[-1]","metadata":{"execution":{"iopub.status.busy":"2022-07-25T11:28:56.379780Z","iopub.execute_input":"2022-07-25T11:28:56.380101Z","iopub.status.idle":"2022-07-25T11:28:56.386990Z","shell.execute_reply.started":"2022-07-25T11:28:56.380070Z","shell.execute_reply":"2022-07-25T11:28:56.385922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction and evaluation on the training set and the test set\ny_train_prob_reg = logistic(np.matmul(X_train.to_numpy(), w_out_reg) + (b_out_reg * np.ones(X_train.shape[0])))\ny_test_prob_reg = logistic(np.matmul(X_test.to_numpy(), w_out_reg) + (b_out_reg * np.ones(X_test.shape[0])))\ny_train_pred_reg, y_test_pred_reg = (y_train_prob_reg > 0.5).astype(int), (y_test_prob_reg > 0.5).astype(int)\nprint(pd.Series({\"Training accuracy\": accuracy(y_train, y_train_pred_reg),\n                 \"Test accuracy\": accuracy(y_test, y_test_pred_reg)}).to_string())","metadata":{"execution":{"iopub.status.busy":"2022-07-25T11:28:56.388663Z","iopub.execute_input":"2022-07-25T11:28:56.389301Z","iopub.status.idle":"2022-07-25T11:28:56.698331Z","shell.execute_reply.started":"2022-07-25T11:28:56.389268Z","shell.execute_reply":"2022-07-25T11:28:56.697071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Confusion matrix for predictions on the test set\nconf_mat_heatmap(y_test, y_test_pred_reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T11:28:56.699605Z","iopub.execute_input":"2022-07-25T11:28:56.700558Z","iopub.status.idle":"2022-07-25T11:28:57.102724Z","shell.execute_reply.started":"2022-07-25T11:28:56.700521Z","shell.execute_reply":"2022-07-25T11:28:57.101666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Thus, in this scenario, regularization does not help much. This is expected since the unregularized implementation does not indicate overfitting.","metadata":{}},{"cell_type":"markdown","source":"# Acknowledgements\n\n- [Higgs Boson Machine Learning Challenge](https://www.kaggle.com/competitions/higgs-boson)\n- [The HiggsML challenge documentation](https://higgsml.ijclab.in2p3.fr/documentation/)","metadata":{}},{"cell_type":"markdown","source":"# References\n\n- [AdaBoost](https://en.wikipedia.org/wiki/AdaBoost)\n- [Algorithm](https://en.wikipedia.org/wiki/Algorithm)\n- [Augustin-Louis Cauchy](https://en.wikipedia.org/wiki/Augustin-Louis_Cauchy)\n- [Binary classification](https://en.wikipedia.org/wiki/Binary_classification)\n- [Bounded function](https://en.wikipedia.org/wiki/Bounded_function)\n- [Classification rule](https://en.wikipedia.org/wiki/Classification_rule)\n- [Decision tree](https://en.wikipedia.org/wiki/Decision_tree_learning)\n- [Differentiable function](https://en.wikipedia.org/wiki/Differentiable_function)\n- [Dot product](https://en.wikipedia.org/wiki/Dot_product)\n- [Evaluation metric](https://en.wikipedia.org/wiki/Evaluation_of_binary_classifiers)\n- [First order methods](https://en.wikipedia.org/wiki/Category:First_order_methods)\n- [Fundamental interaction](https://en.wikipedia.org/wiki/Fundamental_interaction)\n- [Gradient](https://en.wikipedia.org/wiki/Gradient)\n- [Gradient descent algorithm](https://en.wikipedia.org/wiki/Gradient_descent)\n- [Haskell Curry](https://en.wikipedia.org/wiki/Haskell_Curry)\n- [Inflection point](https://en.wikipedia.org/wiki/Inflection_point)\n- [Iterative method](https://en.wikipedia.org/wiki/Iterative_method)\n- [Jacques Hadamard](https://en.wikipedia.org/wiki/Jacques_Hadamard)\n- [$k$-nearest neighbors](https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm)\n- [$l^2$ norm](https://mathworld.wolfram.com/L2-Norm.html)\n- [Learning rate](https://en.wikipedia.org/wiki/Learning_rate)\n- [Linear discriminant analysis](https://en.wikipedia.org/wiki/Linear_discriminant_analysis)\n- [Logistic regression](https://en.wikipedia.org/wiki/Logistic_regression)\n- [Machine learning](https://en.wikipedia.org/wiki/Machine_learning)\n- [Mathematical optimization](https://en.wikipedia.org/wiki/Mathematical_optimization)\n- [Matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication)\n- [Maxima and minima](https://en.wikipedia.org/wiki/Maxima_and_minima)\n- [Naive Bayes](https://en.wikipedia.org/wiki/Naive_Bayes_classifier)\n- [Neural network](https://en.wikipedia.org/wiki/Artificial_neural_network)\n- [Particle physics](https://en.wikipedia.org/wiki/Particle_physics)\n- [Random forest](https://en.wikipedia.org/wiki/Random_forest)\n- [Scalar multiplication](https://en.wikipedia.org/wiki/Scalar_multiplication)\n- [Scikit-learn](https://en.wikipedia.org/wiki/Scikit-learn)\n- [Sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function)\n- [Sorting hat](https://en.wikipedia.org/wiki/Magical_objects_in_Harry_Potter#Sorting_Hat)\n- [Statistical classification](https://en.wikipedia.org/wiki/Statistical_classification)\n- [Statistics](https://en.wikipedia.org/wiki/Statistics)\n- [Stochastic gradient descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent)\n- [Subatomic particle](https://en.wikipedia.org/wiki/Subatomic_particle)\n- [Supervised learning](https://en.wikipedia.org/wiki/Supervised_learning)\n- [Support vector machine](https://en.wikipedia.org/wiki/Support-vector_machine)\n- [Test data](https://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets#Test_data_set)\n- [Training data](https://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets#Training_data_set)\n- [Vector addition](https://en.wikipedia.org/wiki/Euclidean_vector#Addition_and_subtraction)\n- [XGBoost](https://en.wikipedia.org/wiki/XGBoost)\n","metadata":{}},{"cell_type":"code","source":"# Runtime and memory usage\nstop = time.time()\nprint(pd.Series({\"Process runtime\": \"{:.2f} seconds\".format(float(stop - start)),\n                 \"Process memory usage\": \"{:.2f} MB\".format(float(process.memory_info()[0]/(1024*1024)))}).to_string())","metadata":{"execution":{"iopub.status.busy":"2022-07-25T11:28:57.104202Z","iopub.execute_input":"2022-07-25T11:28:57.104510Z","iopub.status.idle":"2022-07-25T11:28:57.112725Z","shell.execute_reply.started":"2022-07-25T11:28:57.104483Z","shell.execute_reply":"2022-07-25T11:28:57.111614Z"},"trusted":true},"execution_count":null,"outputs":[]}]}