{ "cells": [ { "cell_type": "markdown", "id": "556bce9e-f718-401f-a3f2-3f9cf6fad8ca", "metadata": {}, "source": [ "# Logistic Regression From Scratch \n", "--- \n", "\n", "This is a full on implementation of the classification algorithm **Logistic Regression** from scratch using only **Numpy**, including :\n", "- Logistic Function\n", "- Gradient Descent\n", " - Full Batch\n", " - Mini-Batch\n", " - Stocastic\n", "- Cross entropy Loss (Cost function)\n", "- Prediction and accuracy\n", "- Ridge Logistic Regression **L2**\n", "- Lasso Logistic Regression **L1**\n", "- Cross Validation\n", "\n", "**Note**: \n", "- This notebook only covers the basic concepts and functions for **Logistic Regression** \n", "- There will be another jupyter notebook for the **Heart Disease**, **Breast Cancer** and other data set + **sklearn** benchmark)" ] }, { "cell_type": "markdown", "id": "4d196567-b405-453c-b82f-3276a12bf719", "metadata": {}, "source": [ "## Implementation " ] }, { "cell_type": "code", "execution_count": 2, "id": "5406c920-c19d-4d9e-865c-cad50fab47ef", "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "id": "2d0278fa-d060-4eda-97d0-af8eec603aee", "metadata": {}, "source": [ "- We only need **numPy** to implement every logistic regression functionality in this project\n", "- The class implementation based can be found in the `Logistic_Regression_Scratch.py` File which contain better code practices similar to `sklearn` library " ] }, { "cell_type": "markdown", "id": "de08a418-a655-4c79-b7ae-c8f9b6a7f9d1", "metadata": {}, "source": [ "---\n", "### Generating Data" ] }, { "cell_type": "code", "execution_count": 2, "id": "b8de0333-6c6c-481f-8d66-2e0684da28a7", "metadata": {}, "outputs": [], "source": [ "def generate_dummy_data(n=1000,p=3,seed=12):\n", " np.random.seed(seed)\n", " rng = np.random.randomState(seed)\n", " \n", " X= np.random.randn(n,p)\n", " noise = rng.normal(0,0.2,size=X.shape)\n", " X+=noise\n", " X =(X-np.mean(X,axis=0))/np.std(X,axis=0)\n", " \n", " \n", " coeff_true = np.array([0.6]*p)\n", " intercept = 3 \n", "\n", " P_x = np.exp(X @ coeff_true + intercept)/(1+ np.exp(X@coeff_true + intercept))\n", " \n", " y = np.random.binomial(1,P_x)\n", "\n", " return X,coeff_true, intercept , y" ] }, { "cell_type": "markdown", "id": "6060af39-8c97-4579-ba39-a9e31ad1fca8", "metadata": {}, "source": [ "- This function serves as a way to generate dummy data to set our function\n", "- `n` number of **observations** and `p` number of **features** or predictors\n", "- `noise` to simulate real world data and normalized $X$ for stability \n", "- Binary Logistic Regression follow's a **Binomial** distribution so the **true responses** are randomly drawn from it with a probability $P(x)$\n" ] }, { "cell_type": "markdown", "id": "8875e661-f439-42a9-b7b0-ac9e1775d8e5", "metadata": {}, "source": [ "$$ P(X)=\\frac{e^{\\beta_{0}+X\\beta}}{1+e^{\\beta_{0}+X\\beta}}=\\frac{1}{e^{-(\\beta+X\\beta)}+1}$$" ] }, { "cell_type": "markdown", "id": "a976987b-0b4d-4fa5-92d9-8336415a5a49", "metadata": {}, "source": [ "- This is called the **Logistic Function** which gives us results between $0$ and $1$ also called the sigmoid function\n", "- $P(x)$ gives us probability results for each observation " ] }, { "cell_type": "markdown", "id": "6a5b656a-2f60-4347-8a70-2134c260d09b", "metadata": {}, "source": [ "---\n", "### Logistic Function" ] }, { "cell_type": "code", "execution_count": 3, "id": "4858dcd1-6743-4efa-9626-dc40064241e9", "metadata": {}, "outputs": [], "source": [ "def sigmoid_function(X,intercept,beta):\n", " logit = X @ beta + intercept\n", "\n", " P_x = np.where(\n", " logit >= 0,\n", " 1 / (1 + np.exp(-logit)),\n", " np.exp(logit) / (1 + np.exp(logit))\n", " )\n", " return P_x" ] }, { "cell_type": "markdown", "id": "b2d6b095-0bd4-47ef-a0f5-1da0eae35de9", "metadata": {}, "source": [ "- This `sigmoid_function` calculated the probability of an observation as stated above\n", "- The `logit` is just log of the **odds**\n", "$$ odds =p(X)/(1-p(X))$$" ] }, { "cell_type": "markdown", "id": "a2492e62-735d-4553-a51b-2cf16018c617", "metadata": {}, "source": [ "- The`logit` is just a **Linear Regression** equation which allow us to do **inference** and statistical analysis on **Logistic Regression** " ] }, { "cell_type": "markdown", "id": "ff3ed4e6-c818-41e9-925a-a26986664320", "metadata": {}, "source": [ "---\n", "### Cross Entropy Loss (Cost Function)" ] }, { "cell_type": "code", "execution_count": 4, "id": "eecc8381-8d4c-41b4-bd6a-00880e555eff", "metadata": {}, "outputs": [], "source": [ "def cross_entropy_loss(y,X,beta,intercept):\n", " eps = 1e-15\n", " sigmoid_fn = sigmoid_function(X,intercept,beta)\n", " sigmoid_fn = np.clip(sigmoid_fn, eps, 1 - eps)\n", " \n", " cost_fn = -np.mean(y.T*np.log(sigmoid_fn)+(1-y).T*np.log(1-sigmoid_fn))\n", " return cost_fn" ] }, { "cell_type": "markdown", "id": "d04d9c7f-89ab-4228-80ed-99620b1e3111", "metadata": {}, "source": [ "- The cost function for the **Logistic Regression** is called **cross_entropy_loss** given by : $$\\mathcal{l}(\\beta)=-[y^T \\log(p(X))+(1-y)^T\\log(1-p(X))] $$" ] }, { "cell_type": "markdown", "id": "a0a36344-00fa-4579-a29b-b110582b13b1", "metadata": {}, "source": [ "- it's simply the log likelihood of the **maximum likelihood** function\n", "- The **ML** is similar to the binomial distribution **PMF**\n", "- The `-` on the equation is simply for optimizaiton purpose to apply the **Gradient Descent**
\n", "(more information and details on the documentation pdf)" ] }, { "cell_type": "markdown", "id": "710ddf78-3f24-4474-8306-926df926de1a", "metadata": {}, "source": [ "---\n", "### Gradient Descent" ] }, { "cell_type": "markdown", "id": "3d8d4df9-701f-4c43-b837-927f01fdfcd3", "metadata": {}, "source": [ "- Time to fit our logistic regression and estimate the coefficients $beta$, This function will apply all 3 types of known gradient descent " ] }, { "cell_type": "code", "execution_count": 5, "id": "23408ef0-062a-43d5-b032-3b95ca784667", "metadata": {}, "outputs": [], "source": [ "def gradient_descent(lr,n_itr,batch_size,Y,X,n):\n", " p = X.shape[1]\n", " beta_est = np.zeros((p,1))\n", " intercept_est = 0\n", "\n", " for i in range(n_itr):\n", " idx = np.random.choice(n,size = batch_size , replace = False)\n", " X_GD = X[idx]\n", " Y_GD = Y[idx].reshape(-1,1)\n", "\n", " sigmoid_fn = sigmoid_function(X_GD,intercept_est,beta_est)\n", "\n", " gradient_cel = (1/batch_size)*(X_GD.T@(sigmoid_fn-Y_GD))\n", "\n", " gradient_intercept = (1/batch_size)*np.sum(sigmoid_fn-Y_GD)\n", "\n", " beta_est = beta_est - (lr*gradient_cel)\n", " intercept_est = intercept_est - (lr*gradient_intercept)\n", "\n", " return beta_est , intercept_est" ] }, { "cell_type": "markdown", "id": "ecb8639f-a987-43d9-a8e8-7bd715be9e66", "metadata": {}, "source": [ "- The **Logistic Regression** has no closed form solution unlike the **Linear Regression** OLS, so gradient descent is the only way to estimate the coefficients of the model, the gradient of the cost function (cross entropy loss) is : $$ \\nabla J(\\beta)=\\frac{1}{n}X^T(\\sigma(X\\beta)-y)$$\n", "\n", "And for the intercept it's : $$ \\nabla J(\\beta_{0})=\\frac{1}{n}\\sum(\\sigma(X\\beta)-y)$$" ] }, { "cell_type": "markdown", "id": "af9b6474-47d6-4225-9265-1bbeedd4573b", "metadata": {}, "source": [ "- They are simply the pratial derivaiton with respect to $\\beta$ and for the intercept for $\\beta_{0}$" ] }, { "cell_type": "markdown", "id": "eee85ef8-06ff-4c04-9198-8a2a7cd816a7", "metadata": {}, "source": [ "- The `idx` is simply randomly samples take batches from the data `n`\n", "- Both of `X_GD` and `Y_GD` are samples to used to calculate the gradient for the next step in the **Gradient Descent** algorithm " ] }, { "cell_type": "markdown", "id": "c56b74a3-20e8-4362-82ec-1cd41c209e9f", "metadata": {}, "source": [ "---\n", "### Prediction & Accuracy " ] }, { "cell_type": "code", "execution_count": 6, "id": "d5a854ce-4e23-4687-a7c6-4d4c1f52c67f", "metadata": {}, "outputs": [], "source": [ "def predict(X,intercept,beta,threshold=0.5):\n", " predicted_probability = sigmoid_function(X,intercept,beta)\n", "\n", " predicted_class = (predicted_probability>= threshold).astype(int)\n", " return predicted_class, predicted_probability" ] }, { "cell_type": "markdown", "id": "c7c66d1e-5176-4165-a8e4-a0ba95a0213b", "metadata": {}, "source": [ "- This function simply calculate the probability of each observation using the estimated coefficients from the `gradient_descent` function\n", "- Classify based on a `threshold` usually set to $0.5$ to either $0$ or $1$" ] }, { "cell_type": "code", "execution_count": 7, "id": "75f1cc5c-23a6-4044-a40e-41bf5d387d43", "metadata": {}, "outputs": [], "source": [ "def accuracy(Y,Y_pred):\n", " return np.mean(Y_pred==Y)" ] }, { "cell_type": "markdown", "id": "12f52e53-8438-4466-83ec-0f7e6a6e92c6", "metadata": {}, "source": [ "- comparing the true values of the response `Y` and the predicted values of `Y_pred`\n", "- This will come in handy when we compare different regularizations and hyperparameters" ] }, { "cell_type": "markdown", "id": "eba8fa77-8433-43d0-95bc-1f7defa27dcd", "metadata": {}, "source": [ "---\n", "### Ridge Logistic Regression (L2)" ] }, { "cell_type": "markdown", "id": "e67a0904-baef-41a8-adbb-6c49b250116f", "metadata": {}, "source": [ "The same concept as in ridge regression, it **shrinks** the coefficients to prevent overfitting and the model being senstive by introducing some **bias** on the training phase for better **variance** on the predictions, which improve prediction accuracy and reduce colinearity " ] }, { "cell_type": "code", "execution_count": 8, "id": "1f37cd6f-fd9e-4e75-83f3-b6b23a76fd4f", "metadata": {}, "outputs": [], "source": [ "def gradient_descent(penalty,lr,n_itr,batch_size,pen,Y,X,n):\n", " p = X.shape[1]\n", " beta_est = np.zeros((p,1))\n", " intercept_est = 0\n", "\n", " for i in range(n_itr):\n", " idx = np.random.choice(n,size = batch_size , replace = False)\n", " X_GD = X[idx]\n", " Y_GD = Y[idx].reshape(-1,1)\n", "\n", " sigmoid_fn = sigmoid_function(X_GD,intercept_est,beta_est)\n", " if penalty.lower() == \"l2\":\n", " penalty_term =pen* beta_est\n", " gradient_cel = (1/batch_size)*((X_GD.T@(sigmoid_fn-Y_GD))+penalty_term)\n", " else:\n", " gradient_cel = (1/batch_size)*(X_GD.T@(sigmoid_fn-Y_GD))\n", " \n", " gradient_intercept = (1/batch_size)*np.sum(sigmoid_fn-Y_GD)\n", "\n", " beta_est = beta_est - (lr*gradient_cel)\n", " intercept_est = intercept_est - (lr*gradient_intercept)\n", "\n", " return beta_est , intercept_est" ] }, { "cell_type": "markdown", "id": "77979c11-bcff-4cf7-9e49-fde89f9a617c", "metadata": {}, "source": [ "The ridge logistic regression is simply the cross entropy loss function plus the **penalty term** also known as **regularization term** :\n", "$$-l_{ridge}(\\beta,\\lambda)=l(\\beta)-\\frac{\\lambda}{2}\\lVert \\beta \\rVert_{2}^2 $$\n", "- $l(\\beta)$ is the cost function\n", "- $\\lambda$ is the hyperprameter, it's selected using the **Cross Validaiton** method\n", "\n", "The gradient of the **ridge logistic regression cost function** is : $$ \\nabla J(\\beta;\\lambda)=-\\frac{1}{n}X^\\intercal(y-\\sigma(X\\beta))+\\lambda \\beta $$" ] }, { "cell_type": "markdown", "id": "c1cd058f-308e-464a-af0d-df6261074922", "metadata": {}, "source": [ "**Note** : The regulzarization is only applied to the coefficients and not the intercept(bias), cause the ridge regression introduce **bias** to the ceofficients and we dont want our baseline $\\beta_{0}$ to be bais and most of the time it doesnt cause overfitting or get corelated with other coefficients $\\beta$" ] }, { "cell_type": "markdown", "id": "c701feb4-92be-4f36-adae-0fc512c7a47b", "metadata": {}, "source": [ "--- \n", "### Lasso Logistic Regression (L1)" ] }, { "cell_type": "markdown", "id": "0daf0aaf-a4cc-499d-ac7d-4e9b7bad8294", "metadata": {}, "source": [ "The **Lasso Logistic Regression** is yet similar to the **Lasso Regression** which uses the $L_{1}$ Norm as a penalty term which is the **absolute value** of the coefficients $\\beta$, the lasso logistic regression cost function is given : $$ l_{ridge}(\\beta,\\lambda)=l(\\beta)-{\\lambda}\\lVert \\beta \\rVert_{1} $$\n", "\n", "With : $${\\lambda}\\lVert \\beta \\rVert_{1} = \\lambda \\sum_{j=1}^{p} |\\beta_{j}| $$\n", "\n", "- The special thing that the lasso offers is **Sparsity** due to the $L_{1}$ it set some of the coeffcients to **Zero**, which make it useful for **Feature** selection and promoting sparse and interpretble models with less complexity\n" ] }, { "cell_type": "markdown", "id": "40cc239c-21d7-471b-8d82-f5ff30022a6f", "metadata": {}, "source": [ "The gradient of the **Lasso logistic regression** cost function is :$$ \\nabla J(\\beta;\\lambda)=-\\frac{1}{n}X^\\intercal(y-\\sigma(X\\beta))+\\lambda \\text{Sign }(\\beta)$$\n", "\n", "- The **absolute value** doesn't have a solution at zero, so we use **Subgradient** and **KKT** stationarity condition to prove that it set some coefficients to zero which promote sparsity in the model\n", "- The sub-gradient of $\\partial \\lVert \\beta \\rVert_{1} = \\begin{cases}\n", "\\text{sign}(\\hat{\\beta_{j}}) & \\quad , \\hat{\\beta}_{j} \\neq 0 \\\\\n", "\\in[-1,1] & \\quad ,\\hat{\\beta}_{j}=0\n", "\\end{cases}$\n", "\n", "**Intuition** : \n", "- there exist many slops for the point zero and if we draw them all we will get the interval $[-1,1]$ when the coefficient $\\beta =0$" ] }, { "cell_type": "code", "execution_count": 9, "id": "ea8c8ab4-8f5b-42f5-9772-1c108cee9e84", "metadata": {}, "outputs": [], "source": [ "def gradient_descent(penalty,lr,n_itr,batch_size,pen,Y,X):\n", " p = X.shape[1]\n", " beta_est = np.zeros((p,1))\n", " intercept_est = 0\n", " n = X.shape[0]\n", "\n", " for i in range(n_itr):\n", " idx = np.random.choice(n,size = batch_size , replace = False)\n", " X_GD = X[idx]\n", " Y_GD = Y[idx].reshape(-1,1)\n", "\n", " sigmoid_fn = sigmoid_function(X_GD,intercept_est,beta_est)\n", " if penalty.lower() == \"l2\":\n", " penalty_term =pen* beta_est\n", " gradient_cel = (1/batch_size)*((X_GD.T@(sigmoid_fn-Y_GD))+penalty_term)\n", " elif penalty.lower() == \"l1\":\n", " penalty_term = pen * np.sign(beta_est)\n", " gradient_cel = (1/batch_size)*((X_GD.T@(sigmoid_fn- Y_GD))+penalty_term)\n", " else:\n", " gradient_cel = (1/batch_size)*(X_GD.T@(sigmoid_fn-Y_GD))\n", " \n", " gradient_intercept = (1/batch_size)*np.sum(sigmoid_fn-Y_GD)\n", "\n", " beta_est = beta_est - (lr*gradient_cel)\n", " intercept_est = intercept_est - (lr*gradient_intercept)\n", "\n", " return beta_est , intercept_est" ] }, { "cell_type": "markdown", "id": "9bdacee5-30d0-498b-9375-3277e3ee4b77", "metadata": {}, "source": [ "- The `np.sign` returns $1$ if `beta_est` is positive\n", " - $-1$ if `beta_est` is negative\n", " - $0$ if `beta_est` is zero " ] }, { "cell_type": "markdown", "id": "2b4fa553-c48d-480c-9f9f-1184a2d26eca", "metadata": {}, "source": [ "---\n", "### Cross-Validation" ] }, { "cell_type": "markdown", "id": "c3a12315-57e4-4a54-bccf-b08cc8f2fb25", "metadata": {}, "source": [ "**Cross-Validaiton** plays an important role, Since it estimates the **test error rate** which helps us in :\n", "- Model Selection\n", "- Detect Overfitting\n", "- Hyperparameter Tunning ($\\lambda$ for the regularization)\n", "- Assessing Model Stability and Detect unsual behavoirs" ] }, { "cell_type": "code", "execution_count": 23, "id": "19916654-08c8-4cc7-8739-715814e9e159", "metadata": {}, "outputs": [], "source": [ "def K_folds(n,K):\n", " n = np.arange(n)\n", " folds = []\n", " fold_size= len(n)//K\n", " #print(\"this is fold size \" , fold_size)\n", " for i in range(K):\n", " start = i*fold_size\n", " end = (i+1)*fold_size if i< K-1 else len(n)\n", " test_idx = n[start:end]\n", " train_idx = np.concatenate([n[:start],n[end:]])\n", " folds.append((train_idx,test_idx))\n", "\n", " return folds" ] }, { "cell_type": "code", "execution_count": 28, "id": "b670ba5d-e880-4556-a1a2-e06bf08bb411", "metadata": {}, "outputs": [], "source": [ "def strat_Kfolds(n,K,y):\n", " class_idx = {}\n", " for class_label in np.unique(y):\n", " class_idx[class_label] =np.where(y==class_label)[0]\n", "\n", " #print(class_idx[1])\n", " #print(class_idx[0])\n", " folds =[([],[]) for _ in range(K)]\n", "\n", " for class_label,idx in class_idx.items():\n", " class_samples = len(idx)\n", " #print(f\"the class {class_label} has {class_samples}\")\n", " fold_size_class = class_samples //K\n", " #print(f\"the class {class_label} has fold size of {fold_size_class}\")\n", "\n", " class_fold = []\n", " for i in range (K):\n", " start = i* fold_size_class \n", " end = (i+1)*fold_size_class if i < K-1 else class_samples\n", " test_idx_class = idx[start:end]\n", " train_idx_class = np.concatenate([idx[:start],idx[end:]])\n", " folds[i][0].extend(train_idx_class)\n", " folds[i][1].extend(test_idx_class)\n", " #print(folds[1][0])\n", " folds = [(np.array(train),np.array(test)) for train , test in folds]\n", "\n", " return folds\n" ] }, { "cell_type": "markdown", "id": "6d8a4bd5-4056-4cd4-96cf-2aaf80a202a0", "metadata": {}, "source": [ "**K-Fold Cross Validation** is the most widely used approach, Since it's a generalized case of the **LOOCV** by splitting the data into :\n", "- $K$ Folds\n", "- Train on $\\to K-1$\n", "- Validate and test on $\\to 1$\n", "\n", "Here we calculate `start` and `end` to construct our `test_idx` indcies for the test data and `train_idx` for the training indices:\n", "- The indices are stored in `folds` list which will be used on the `cross_validaition` function\n" ] }, { "cell_type": "code", "execution_count": 34, "id": "c3145488-e802-41cf-8478-2a44342779b3", "metadata": {}, "outputs": [], "source": [ "def cross_validation(K,X,Y):\n", " n = X.shape[0]\n", " #folds = strat_Kfolds(n,K,Y)\n", " \n", " folds = K_folds(n,K)\n", " losses = []\n", " for i, (train,test)in enumerate(folds):\n", " X_train = X[train]\n", " Y_train = Y[train].reshape(-1,1)\n", " X_test = X[test]\n", " Y_test = Y[test].reshape(-1,1)\n", " # print(\"this is shape of X and Y train\",X_train.shape,Y_train.shape)\n", " # print(\"this is shape of X and Y test\",X_test.shape,Y_test.shape)\n", " beta_est , intercept_est = gradient_descent(penalty=\"None\",pen=\"1\",lr=0.001,n_itr=200,batch_size=2,Y=Y_train,X=X_train) \n", " loss = cross_entropy_loss(Y_test ,X_test,beta_est,intercept_est)\n", " losses.append(loss)\n", " return losses , np.mean(losses) " ] }, { "cell_type": "markdown", "id": "864ea739-0523-4163-8107-901ad92b03cb", "metadata": {}, "source": [ "The **Cross-Validation** algorithm steps are : \n", "- Split the data into **K-Folds**\n", "- Fit(Train) the model on the training folds\n", "- Validate by calculating the loss\n", "- Repeat for all $\\text{K-Folds}$ and each time on Fold act as **test fold** while the others as training sets $K-1$" ] }, { "cell_type": "code", "execution_count": 41, "id": "1312fa2d-afea-4152-80cd-09c582c09747", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1000, 3)" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X,beta,inter,y = generate_dummy_data()\n", "X.shape" ] }, { "cell_type": "code", "execution_count": 1, "id": "10eeeb6d-7cb2-4a4c-a39d-5123a4cfc0bb", "metadata": { "scrolled": true }, "outputs": [ { "ename": "NameError", "evalue": "name 'cross_validation' is not defined", "output_type": "error", "traceback": [ "\u001b[31m---------------------------------------------------------------------------\u001b[39m", "\u001b[31mNameError\u001b[39m Traceback (most recent call last)", "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m losses,avg_loss = \u001b[43mcross_validation\u001b[49m(\u001b[32m1000\u001b[39m,X,y)\n\u001b[32m 2\u001b[39m avg_loss\n", "\u001b[31mNameError\u001b[39m: name 'cross_validation' is not defined" ] } ], "source": [ "losses,avg_loss = cross_validation(1000,X,y)\n", "avg_loss" ] }, { "cell_type": "code", "execution_count": null, "id": "42c6a48c-e26b-4931-a937-c8202c0338ab", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }