{"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":"# Logistic Regression from scratch\nIn this notebook i solved a binary classification problem: \"Titanic - Machine Learning from Disaster\" from Kaggle. To solve the problem i used a Stochastic Gradient Descent implementation of Logistic Regression, and uses the Sigmoid probablity outcome to classify an input by applying a threshold to its probablity of belonging to the positive class.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nfrom sklearn import metrics\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import train_test_split","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-18T20:34:28.538811Z","iopub.execute_input":"2022-07-18T20:34:28.539311Z","iopub.status.idle":"2022-07-18T20:34:28.545797Z","shell.execute_reply.started":"2022-07-18T20:34:28.539274Z","shell.execute_reply":"2022-07-18T20:34:28.544604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Logistic Regression implementation\nUsing stochastic gradient descent","metadata":{}},{"cell_type":"code","source":"class LogisticRegressor:\n    def __init__(self, threshold, lr=0.001):\n        self.lr = lr # Learning rate\n        self.threshold = threshold # Classification threshold\n        self._epoch_loss = np.empty(shape=(1,2)) # Array for stacking the loss by epoch\n    \n    @staticmethod\n    def _log_loss(y, y_pred):\n        return -1 * np.mean(y*np.log(y_pred) + (1-y)*np.log(1-y_pred))\n    \n    @staticmethod\n    def _sigmoid(z):\n        '''Sigmoid function'''\n        return 1 / (1 + np.exp(-z))\n    \n    @staticmethod\n    def _z(X, w, b):\n        '''Linear function to give to sigmoid as param'''\n        return X.dot(w) + b\n    \n    @staticmethod\n    def _sigmoid_dz(z):\n        '''Derivative of sigmoid in respect to z'''\n        return ( (np.exp(-z)) / (1 + np.exp(-z))**2 )\n    \n    @staticmethod\n    def _loss_dy(y, y_pred):\n        '''Derivative of log loss in respect to model predictions i.e sigmoid(z)'''\n        return (-y/y_pred) + ( (1-y)/(1-y_pred) )\n    \n    def _plot_loss(self):\n        plt.plot(\n            self._epoch_loss[1:,0],\n            self._epoch_loss[1:,1]\n        )\n        plt.title('Training Loss')\n        plt.xlabel('Epoch')\n        plt.ylabel('Loss')\n        plt.show()\n    \n    def fit(self, X, y):\n        self.w = np.ones(X.shape[1]) # Init weights\n        self.b = 0 # Init Bias\n        \n        # Stochasitc gradient descent\n        for epoch in range(100): # Epochs\n            \n            # Storing a history of epoch_n - loss pairs for plotting pourposes\n            z_epoch = self._z(X, self.w, self.b)\n            y_pred_epoch = self._sigmoid(z_epoch)\n            epoch_loss = self._log_loss(y, y_pred_epoch)\n            self._epoch_loss = np.vstack(\n                (self._epoch_loss, \n                np.array([epoch, epoch_loss]))\n            )\n                \n            for i in range(X.shape[0]): # Examples\n                z = self._z(X[i], self.w, self.b)\n                y_pred = self._sigmoid(z)\n                dc_dz = self._loss_dy(y[i], y_pred) * self._sigmoid_dz(z)\n                \n                # Updating weights with cost derivative in respect to a given weight\n                self.w = self.w - (self.lr * (dc_dz*X[i]) )\n                self.b = self.b - (self.lr * dc_dz)\n        self._plot_loss()\n                                \n    def predict(self, x):\n        '''Predict probablities'''\n        return self._sigmoid(self._z(x, self.w, self.b))\n    \n    def classify(self, x):\n        '''Classify in function of probabilities'''\n        prob = self.predict(x)\n        # Positive class if prob >= threshold, else negative class\n        return np.array(list(map(lambda a: 1 if a >= self.threshold else 0, prob)))","metadata":{"execution":{"iopub.status.busy":"2022-07-18T20:34:28.661823Z","iopub.execute_input":"2022-07-18T20:34:28.662633Z","iopub.status.idle":"2022-07-18T20:34:28.684219Z","shell.execute_reply.started":"2022-07-18T20:34:28.662589Z","shell.execute_reply":"2022-07-18T20:34:28.682621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data preparation\nTurning raw data into useful features","metadata":{}},{"cell_type":"code","source":"columns_to_exclude = ['Name','Ticket', 'Fare', 'Embarked', 'PassengerId', 'Cabin']\ndf = pd.read_csv('/kaggle/input/titanic/train.csv')\n\n# Dropping columns with no predicting power\ndf.drop(labels=columns_to_exclude, inplace=True, axis=1)\n\n# Filling NaN values with colomn mean for age column.\ndf.fillna(value=df['Age'].mean(), inplace=True)\n\n# Using Z-score normalization for Age feature so as to have a same range as the other features\ndf['Age'] = (df['Age'] - df['Age'].mean()) / df['Age'].std()\n\n# One hot encoding for Sex feature\nsex_one_hot = pd.get_dummies(df['Sex'], prefix='is')\ndf = df.join(sex_one_hot)\ndf.drop(labels='Sex', axis=1, inplace=True)\n\n# One hot encoding for Pclass feature\n# If we let the model learn a proper weight of belonging to every passanger class in dataset\n# this feature will have a more predictive power\npclass_one_hot = pd.get_dummies(df['Pclass'], prefix='is_from_class')\ndf = df.join(pclass_one_hot)\ndf.drop(labels='Pclass', axis=1, inplace=True)\n\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-18T20:34:28.774517Z","iopub.execute_input":"2022-07-18T20:34:28.775528Z","iopub.status.idle":"2022-07-18T20:34:28.810305Z","shell.execute_reply.started":"2022-07-18T20:34:28.775480Z","shell.execute_reply":"2022-07-18T20:34:28.809123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Splitting data into test and validation sets","metadata":{}},{"cell_type":"code","source":"X = df.drop(labels='Survived', axis=1, inplace=False) # Features\ny = df['Survived'].copy() # Label\n\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T20:34:28.871587Z","iopub.execute_input":"2022-07-18T20:34:28.871981Z","iopub.status.idle":"2022-07-18T20:34:28.880879Z","shell.execute_reply.started":"2022-07-18T20:34:28.871949Z","shell.execute_reply":"2022-07-18T20:34:28.880054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training model","metadata":{}},{"cell_type":"code","source":"model = LogisticRegressor(threshold=0.5)\nmodel.fit(X_train.values, y_train.values)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T20:34:28.933420Z","iopub.execute_input":"2022-07-18T20:34:28.934194Z","iopub.status.idle":"2022-07-18T20:34:31.069015Z","shell.execute_reply.started":"2022-07-18T20:34:28.934147Z","shell.execute_reply":"2022-07-18T20:34:31.067625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training model and computing metrics","metadata":{}},{"cell_type":"code","source":"y_pred = model.classify(X_val.values)\ny_probs = model.predict(X_val.values)\n\nz = pd.DataFrame({\n    'True': y_val.values,\n    'Pred': y_pred\n})\n\nprint(\"Accuracy:\", metrics.accuracy_score(y_val.values, y_pred))\nprint(\"Precision:\", metrics.precision_score(y_val.values, y_pred))\nprint(\"Recall:\", metrics.recall_score(y_val.values, y_pred))\n\n# Ploting ROC curve\nfpr, tpr, thresholds = metrics.roc_curve(y_val.values, y_probs)\nauc = metrics.roc_auc_score(y_val.values, y_probs)\nplt.plot(fpr, tpr)\nplt.fill_between(fpr, tpr, alpha=0.1)\nplt.title('ROC curve')\nplt.xlabel('False positive rate')\nplt.ylabel('True positive rate')\nplt.plot([0, 1], [0,1], 'k--', alpha=0.1)\nplt.text(0.3, 0.5, f'AUC = {round(auc, 3)}', fontsize='xx-large')\nplt.show()\n\nz.head(10)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T20:34:31.071632Z","iopub.execute_input":"2022-07-18T20:34:31.072426Z","iopub.status.idle":"2022-07-18T20:34:31.303058Z","shell.execute_reply.started":"2022-07-18T20:34:31.072369Z","shell.execute_reply":"2022-07-18T20:34:31.301796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The final submission of this model to Kaggle competition scored 0.75837 (The higher score possible is 1.0)\n\nhttps://www.kaggle.com/code/sebastianandrade/titanic-challenge?scriptVersionId=100276781","metadata":{}}]}