{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Chapter 3 – Classification**\n", "\n", "_This notebook contains all the sample code and solutions to the exercises in chapter 3._" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's import a few common modules, ensure MatplotLib plots figures inline and prepare a function to save the figures. We also check that Python 3.5 or later is installed (although Python 2.x may work, it is deprecated so we strongly recommend you use Python 3 instead), as well as Scikit-Learn ≥0.20." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Python ≥3.5 is required\n", "import sys\n", "assert sys.version_info >= (3, 5)\n", "\n", "# Scikit-Learn ≥0.20 is required\n", "import sklearn\n", "assert sklearn.__version__ >= \"0.20\"\n", "\n", "# Common imports\n", "import numpy as np\n", "import os\n", "\n", "# to make this notebook's output stable across runs\n", "np.random.seed(42)\n", "\n", "# To plot pretty figures\n", "%matplotlib inline\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "mpl.rc('axes', labelsize=14)\n", "mpl.rc('xtick', labelsize=12)\n", "mpl.rc('ytick', labelsize=12)\n", "\n", "# Where to save the figures\n", "PROJECT_ROOT_DIR = \".\"\n", "CHAPTER_ID = \"classification\"\n", "IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, \"images\", CHAPTER_ID)\n", "os.makedirs(IMAGES_PATH, exist_ok=True)\n", "\n", "def save_fig(fig_id, tight_layout=True, fig_extension=\"png\", resolution=300):\n", " path = os.path.join(IMAGES_PATH, fig_id + \".\" + fig_extension)\n", " print(\"Saving figure\", fig_id)\n", " if tight_layout:\n", " plt.tight_layout()\n", " plt.savefig(path, format=fig_extension, dpi=resolution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MNIST" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from sklearn.datasets import fetch_openml\n", "mnist = fetch_openml('mnist_784', version=1)\n", "mnist.keys()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "X, y = mnist[\"data\"], mnist[\"target\"]\n", "X.shape" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "y.shape" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "28 * 28" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "\n", "some_digit = X[0]\n", "some_digit_image = some_digit.reshape(28, 28)\n", "plt.imshow(some_digit_image, cmap = mpl.cm.binary, interpolation=\"nearest\")\n", "plt.axis(\"off\")\n", "\n", "save_fig(\"some_digit_plot\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "y[0]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "y = y.astype(np.uint8)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def plot_digit(data):\n", " image = data.reshape(28, 28)\n", " plt.imshow(image, cmap = mpl.cm.binary,\n", " interpolation=\"nearest\")\n", " plt.axis(\"off\")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# EXTRA\n", "def plot_digits(instances, images_per_row=10, **options):\n", " size = 28\n", " images_per_row = min(len(instances), images_per_row)\n", " images = [instance.reshape(size,size) for instance in instances]\n", " n_rows = (len(instances) - 1) // images_per_row + 1\n", " row_images = []\n", " n_empty = n_rows * images_per_row - len(instances)\n", " images.append(np.zeros((size, size * n_empty)))\n", " for row in range(n_rows):\n", " rimages = images[row * images_per_row : (row + 1) * images_per_row]\n", " row_images.append(np.concatenate(rimages, axis=1))\n", " image = np.concatenate(row_images, axis=0)\n", " plt.imshow(image, cmap = mpl.cm.binary, **options)\n", " plt.axis(\"off\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(9,9))\n", "example_images = X[:100]\n", "plot_digits(example_images, images_per_row=10)\n", "save_fig(\"more_digits_plot\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "y[0]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "X_train, X_test, y_train, y_test = X[:60000], X[60000:], y[:60000], y[60000:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Binary classifier" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "y_train_5 = (y_train == 5)\n", "y_test_5 = (y_test == 5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: some hyperparameters will have a different defaut value in future versions of Scikit-Learn, such as `max_iter` and `tol`. To be future-proof, we explicitly set these hyperparameters to their future default values. For simplicity, this is not shown in the book." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import SGDClassifier\n", "\n", "sgd_clf = SGDClassifier(max_iter=1000, tol=1e-3, random_state=42)\n", "sgd_clf.fit(X_train, y_train_5)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "sgd_clf.predict([some_digit])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_val_score\n", "cross_val_score(sgd_clf, X_train, y_train_5, cv=3, scoring=\"accuracy\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import StratifiedKFold\n", "from sklearn.base import clone\n", "\n", "skfolds = StratifiedKFold(n_splits=3, random_state=42)\n", "\n", "for train_index, test_index in skfolds.split(X_train, y_train_5):\n", " clone_clf = clone(sgd_clf)\n", " X_train_folds = X_train[train_index]\n", " y_train_folds = y_train_5[train_index]\n", " X_test_fold = X_train[test_index]\n", " y_test_fold = y_train_5[test_index]\n", "\n", " clone_clf.fit(X_train_folds, y_train_folds)\n", " y_pred = clone_clf.predict(X_test_fold)\n", " n_correct = sum(y_pred == y_test_fold)\n", " print(n_correct / len(y_pred))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from sklearn.base import BaseEstimator\n", "class Never5Classifier(BaseEstimator):\n", " def fit(self, X, y=None):\n", " pass\n", " def predict(self, X):\n", " return np.zeros((len(X), 1), dtype=bool)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "never_5_clf = Never5Classifier()\n", "cross_val_score(never_5_clf, X_train, y_train_5, cv=3, scoring=\"accuracy\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_val_predict\n", "\n", "y_train_pred = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import confusion_matrix\n", "\n", "confusion_matrix(y_train_5, y_train_pred)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "y_train_perfect_predictions = y_train_5 # pretend we reached perfection\n", "confusion_matrix(y_train_5, y_train_perfect_predictions)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import precision_score, recall_score\n", "\n", "precision_score(y_train_5, y_train_pred)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "4096 / (4096 + 1522)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "recall_score(y_train_5, y_train_pred)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "4096 / (4096 + 1325)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import f1_score\n", "\n", "f1_score(y_train_5, y_train_pred)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "4096 / (4096 + (1522 + 1325) / 2)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "y_scores = sgd_clf.decision_function([some_digit])\n", "y_scores" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "threshold = 0\n", "y_some_digit_pred = (y_scores > threshold)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "y_some_digit_pred" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "threshold = 8000\n", "y_some_digit_pred = (y_scores > threshold)\n", "y_some_digit_pred" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "y_scores = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3,\n", " method=\"decision_function\")" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import precision_recall_curve\n", "\n", "precisions, recalls, thresholds = precision_recall_curve(y_train_5, y_scores)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "def plot_precision_recall_vs_threshold(precisions, recalls, thresholds):\n", " plt.plot(thresholds, precisions[:-1], \"b--\", label=\"Precision\", linewidth=2)\n", " plt.plot(thresholds, recalls[:-1], \"g-\", label=\"Recall\", linewidth=2)\n", " plt.legend(loc=\"center right\", fontsize=16) # Not shown in the book\n", " plt.xlabel(\"Threshold\", fontsize=16) # Not shown\n", " plt.grid(True) # Not shown\n", " plt.axis([-50000, 50000, 0, 1]) # Not shown\n", "\n", "plt.figure(figsize=(8, 4)) # Not shown\n", "plot_precision_recall_vs_threshold(precisions, recalls, thresholds)\n", "plt.plot([7813, 7813], [0., 0.9], \"r:\") # Not shown\n", "plt.plot([-50000, 7813], [0.9, 0.9], \"r:\") # Not shown\n", "plt.plot([-50000, 7813], [0.4368, 0.4368], \"r:\")# Not shown\n", "plt.plot([7813], [0.9], \"ro\") # Not shown\n", "plt.plot([7813], [0.4368], \"ro\") # Not shown\n", "save_fig(\"precision_recall_vs_threshold_plot\") # Not shown\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "(y_train_pred == (y_scores > 0)).all()" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "def plot_precision_vs_recall(precisions, recalls):\n", " plt.plot(recalls, precisions, \"b-\", linewidth=2)\n", " plt.xlabel(\"Recall\", fontsize=16)\n", " plt.ylabel(\"Precision\", fontsize=16)\n", " plt.axis([0, 1, 0, 1])\n", " plt.grid(True)\n", "\n", "plt.figure(figsize=(8, 6))\n", "plot_precision_vs_recall(precisions, recalls)\n", "plt.plot([0.4368, 0.4368], [0., 0.9], \"r:\")\n", "plt.plot([0.0, 0.4368], [0.9, 0.9], \"r:\")\n", "plt.plot([0.4368], [0.9], \"ro\")\n", "save_fig(\"precision_vs_recall_plot\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "threshold_90_precision = thresholds[np.argmax(precisions >= 0.90)] # == 7813" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "threshold_90_precision" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "y_train_pred_90 = (y_scores >= threshold_90_precision)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "precision_score(y_train_5, y_train_pred_90)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "recall_score(y_train_5, y_train_pred_90)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# ROC curves" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import roc_curve\n", "\n", "fpr, tpr, thresholds = roc_curve(y_train_5, y_scores)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "def plot_roc_curve(fpr, tpr, label=None):\n", " plt.plot(fpr, tpr, linewidth=2, label=label)\n", " plt.plot([0, 1], [0, 1], 'k--') # dashed diagonal\n", " plt.axis([0, 1, 0, 1]) # Not shown in the book\n", " plt.xlabel('False Positive Rate (Fall-Out)', fontsize=16) # Not shown\n", " plt.ylabel('True Positive Rate (Recall)', fontsize=16) # Not shown\n", " plt.grid(True) # Not shown\n", "\n", "plt.figure(figsize=(8, 6)) # Not shown\n", "plot_roc_curve(fpr, tpr)\n", "plt.plot([4.837e-3, 4.837e-3], [0., 0.4368], \"r:\") # Not shown\n", "plt.plot([0.0, 4.837e-3], [0.4368, 0.4368], \"r:\") # Not shown\n", "plt.plot([4.837e-3], [0.4368], \"ro\") # Not shown\n", "save_fig(\"roc_curve_plot\") # Not shown\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import roc_auc_score\n", "\n", "roc_auc_score(y_train_5, y_scores)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: we set `n_estimators=100` to be future-proof since this will be the default value in Scikit-Learn 0.22." ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "from sklearn.ensemble import RandomForestClassifier\n", "forest_clf = RandomForestClassifier(n_estimators=100, random_state=42)\n", "y_probas_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3,\n", " method=\"predict_proba\")" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [], "source": [ "y_scores_forest = y_probas_forest[:, 1] # score = proba of positive class\n", "fpr_forest, tpr_forest, thresholds_forest = roc_curve(y_train_5,y_scores_forest)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8, 6))\n", "plt.plot(fpr, tpr, \"b:\", linewidth=2, label=\"SGD\")\n", "plot_roc_curve(fpr_forest, tpr_forest, \"Random Forest\")\n", "plt.plot([4.837e-3, 4.837e-3], [0., 0.4368], \"r:\")\n", "plt.plot([0.0, 4.837e-3], [0.4368, 0.4368], \"r:\")\n", "plt.plot([4.837e-3], [0.4368], \"ro\")\n", "plt.plot([4.837e-3, 4.837e-3], [0., 0.9487], \"r:\")\n", "plt.plot([4.837e-3], [0.9487], \"ro\")\n", "plt.grid(True)\n", "plt.legend(loc=\"lower right\", fontsize=16)\n", "save_fig(\"roc_curve_comparison_plot\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "roc_auc_score(y_train_5, y_scores_forest)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "y_train_pred_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3)\n", "precision_score(y_train_5, y_train_pred_forest)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "recall_score(y_train_5, y_train_pred_forest)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multiclass classification" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "sgd_clf.fit(X_train, y_train) # y_train, not y_train_5\n", "sgd_clf.predict([some_digit])" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "some_digit_scores = sgd_clf.decision_function([some_digit])\n", "some_digit_scores" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "np.argmax(some_digit_scores)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "sgd_clf.classes_" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "sgd_clf.classes_[5]" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "from sklearn.multiclass import OneVsOneClassifier\n", "ovo_clf = OneVsOneClassifier(SGDClassifier(max_iter=1000, tol=1e-3, random_state=42))\n", "ovo_clf.fit(X_train, y_train)\n", "ovo_clf.predict([some_digit])" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "len(ovo_clf.estimators_)" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "forest_clf.fit(X_train, y_train)\n", "forest_clf.predict([some_digit])" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "forest_clf.predict_proba([some_digit])" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [], "source": [ "cross_val_score(sgd_clf, X_train, y_train, cv=3, scoring=\"accuracy\")" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "from sklearn.preprocessing import StandardScaler\n", "scaler = StandardScaler()\n", "X_train_scaled = scaler.fit_transform(X_train.astype(np.float64))\n", "cross_val_score(sgd_clf, X_train_scaled, y_train, cv=3, scoring=\"accuracy\")" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "y_train_pred = cross_val_predict(sgd_clf, X_train_scaled, y_train, cv=3)\n", "conf_mx = confusion_matrix(y_train, y_train_pred)\n", "conf_mx" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "def plot_confusion_matrix(matrix):\n", " \"\"\"If you prefer color and a colorbar\"\"\"\n", " fig = plt.figure(figsize=(8,8))\n", " ax = fig.add_subplot(111)\n", " cax = ax.matshow(matrix)\n", " fig.colorbar(cax)" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "plt.matshow(conf_mx, cmap=plt.cm.gray)\n", "save_fig(\"confusion_matrix_plot\", tight_layout=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "row_sums = conf_mx.sum(axis=1, keepdims=True)\n", "norm_conf_mx = conf_mx / row_sums" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [], "source": [ "np.fill_diagonal(norm_conf_mx, 0)\n", "plt.matshow(norm_conf_mx, cmap=plt.cm.gray)\n", "save_fig(\"confusion_matrix_errors_plot\", tight_layout=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [], "source": [ "cl_a, cl_b = 3, 5\n", "X_aa = X_train[(y_train == cl_a) & (y_train_pred == cl_a)]\n", "X_ab = X_train[(y_train == cl_a) & (y_train_pred == cl_b)]\n", "X_ba = X_train[(y_train == cl_b) & (y_train_pred == cl_a)]\n", "X_bb = X_train[(y_train == cl_b) & (y_train_pred == cl_b)]\n", "\n", "plt.figure(figsize=(8,8))\n", "plt.subplot(221); plot_digits(X_aa[:25], images_per_row=5)\n", "plt.subplot(222); plot_digits(X_ab[:25], images_per_row=5)\n", "plt.subplot(223); plot_digits(X_ba[:25], images_per_row=5)\n", "plt.subplot(224); plot_digits(X_bb[:25], images_per_row=5)\n", "save_fig(\"error_analysis_digits_plot\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multilabel classification" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [], "source": [ "from sklearn.neighbors import KNeighborsClassifier\n", "\n", "y_train_large = (y_train >= 7)\n", "y_train_odd = (y_train % 2 == 1)\n", "y_multilabel = np.c_[y_train_large, y_train_odd]\n", "\n", "knn_clf = KNeighborsClassifier()\n", "knn_clf.fit(X_train, y_multilabel)" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [], "source": [ "knn_clf.predict([some_digit])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Warning**: the following cell may take a very long time (possibly hours depending on your hardware)." ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [], "source": [ "y_train_knn_pred = cross_val_predict(knn_clf, X_train, y_multilabel, cv=3)\n", "f1_score(y_multilabel, y_train_knn_pred, average=\"macro\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multioutput classification" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "noise = np.random.randint(0, 100, (len(X_train), 784))\n", "X_train_mod = X_train + noise\n", "noise = np.random.randint(0, 100, (len(X_test), 784))\n", "X_test_mod = X_test + noise\n", "y_train_mod = X_train\n", "y_test_mod = X_test" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [], "source": [ "some_index = 0\n", "plt.subplot(121); plot_digit(X_test_mod[some_index])\n", "plt.subplot(122); plot_digit(y_test_mod[some_index])\n", "save_fig(\"noisy_digit_example_plot\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "knn_clf.fit(X_train_mod, y_train_mod)\n", "clean_digit = knn_clf.predict([X_test_mod[some_index]])\n", "plot_digit(clean_digit)\n", "save_fig(\"cleaned_digit_example_plot\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Extra material" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dummy (ie. random) classifier" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": [ "from sklearn.dummy import DummyClassifier\n", "dmy_clf = DummyClassifier()\n", "y_probas_dmy = cross_val_predict(dmy_clf, X_train, y_train_5, cv=3, method=\"predict_proba\")\n", "y_scores_dmy = y_probas_dmy[:, 1]" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "scrolled": true }, "outputs": [], "source": [ "fprr, tprr, thresholdsr = roc_curve(y_train_5, y_scores_dmy)\n", "plot_roc_curve(fprr, tprr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## KNN classifier" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [], "source": [ "from sklearn.neighbors import KNeighborsClassifier\n", "knn_clf = KNeighborsClassifier(weights='distance', n_neighbors=4)\n", "knn_clf.fit(X_train, y_train)" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [], "source": [ "y_knn_pred = knn_clf.predict(X_test)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import accuracy_score\n", "accuracy_score(y_test, y_knn_pred)" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [], "source": [ "from scipy.ndimage.interpolation import shift\n", "def shift_digit(digit_array, dx, dy, new=0):\n", " return shift(digit_array.reshape(28, 28), [dy, dx], cval=new).reshape(784)\n", "\n", "plot_digit(shift_digit(some_digit, 5, 1, new=100))" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "X_train_expanded = [X_train]\n", "y_train_expanded = [y_train]\n", "for dx, dy in ((1, 0), (-1, 0), (0, 1), (0, -1)):\n", " shifted_images = np.apply_along_axis(shift_digit, axis=1, arr=X_train, dx=dx, dy=dy)\n", " X_train_expanded.append(shifted_images)\n", " y_train_expanded.append(y_train)\n", "\n", "X_train_expanded = np.concatenate(X_train_expanded)\n", "y_train_expanded = np.concatenate(y_train_expanded)\n", "X_train_expanded.shape, y_train_expanded.shape" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [], "source": [ "knn_clf.fit(X_train_expanded, y_train_expanded)" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [], "source": [ "y_knn_expanded_pred = knn_clf.predict(X_test)" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [], "source": [ "accuracy_score(y_test, y_knn_expanded_pred)" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [], "source": [ "ambiguous_digit = X_test[2589]\n", "knn_clf.predict_proba([ambiguous_digit])" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [], "source": [ "plot_digit(ambiguous_digit)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Exercise solutions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. An MNIST Classifier With Over 97% Accuracy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Warning**: the next cell may take hours to run, depending on your hardware." ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import GridSearchCV\n", "\n", "param_grid = [{'weights': [\"uniform\", \"distance\"], 'n_neighbors': [3, 4, 5]}]\n", "\n", "knn_clf = KNeighborsClassifier()\n", "grid_search = GridSearchCV(knn_clf, param_grid, cv=5, verbose=3)\n", "grid_search.fit(X_train, y_train)" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [], "source": [ "grid_search.best_params_" ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [], "source": [ "grid_search.best_score_" ] }, { "cell_type": "code", "execution_count": 92, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import accuracy_score\n", "\n", "y_pred = grid_search.predict(X_test)\n", "accuracy_score(y_test, y_pred)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Data Augmentation" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [], "source": [ "from scipy.ndimage.interpolation import shift" ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [], "source": [ "def shift_image(image, dx, dy):\n", " image = image.reshape((28, 28))\n", " shifted_image = shift(image, [dy, dx], cval=0, mode=\"constant\")\n", " return shifted_image.reshape([-1])" ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [], "source": [ "image = X_train[1000]\n", "shifted_image_down = shift_image(image, 0, 5)\n", "shifted_image_left = shift_image(image, -5, 0)\n", "\n", "plt.figure(figsize=(12,3))\n", "plt.subplot(131)\n", "plt.title(\"Original\", fontsize=14)\n", "plt.imshow(image.reshape(28, 28), interpolation=\"nearest\", cmap=\"Greys\")\n", "plt.subplot(132)\n", "plt.title(\"Shifted down\", fontsize=14)\n", "plt.imshow(shifted_image_down.reshape(28, 28), interpolation=\"nearest\", cmap=\"Greys\")\n", "plt.subplot(133)\n", "plt.title(\"Shifted left\", fontsize=14)\n", "plt.imshow(shifted_image_left.reshape(28, 28), interpolation=\"nearest\", cmap=\"Greys\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [], "source": [ "X_train_augmented = [image for image in X_train]\n", "y_train_augmented = [label for label in y_train]\n", "\n", "for dx, dy in ((1, 0), (-1, 0), (0, 1), (0, -1)):\n", " for image, label in zip(X_train, y_train):\n", " X_train_augmented.append(shift_image(image, dx, dy))\n", " y_train_augmented.append(label)\n", "\n", "X_train_augmented = np.array(X_train_augmented)\n", "y_train_augmented = np.array(y_train_augmented)" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [], "source": [ "shuffle_idx = np.random.permutation(len(X_train_augmented))\n", "X_train_augmented = X_train_augmented[shuffle_idx]\n", "y_train_augmented = y_train_augmented[shuffle_idx]" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [], "source": [ "knn_clf = KNeighborsClassifier(**grid_search.best_params_)" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [], "source": [ "knn_clf.fit(X_train_augmented, y_train_augmented)" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [], "source": [ "y_pred = knn_clf.predict(X_test)\n", "accuracy_score(y_test, y_pred)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By simply augmenting the data, we got a 0.5% accuracy boost. :)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Tackle the Titanic dataset" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The goal is to predict whether or not a passenger survived based on attributes such as their age, sex, passenger class, where they embarked and so on." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, login to [Kaggle](https://www.kaggle.com/) and go to the [Titanic challenge](https://www.kaggle.com/c/titanic) to download `train.csv` and `test.csv`. Save them to the `datasets/titanic` directory." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's load the data:" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "TITANIC_PATH = os.path.join(\"datasets\", \"titanic\")" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "def load_titanic_data(filename, titanic_path=TITANIC_PATH):\n", " csv_path = os.path.join(titanic_path, filename)\n", " return pd.read_csv(csv_path)" ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [], "source": [ "train_data = load_titanic_data(\"train.csv\")\n", "test_data = load_titanic_data(\"test.csv\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data is already split into a training set and a test set. However, the test data does *not* contain the labels: your goal is to train the best model you can using the training data, then make your predictions on the test data and upload them to Kaggle to see your final score." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a peek at the top few rows of the training set:" ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [], "source": [ "train_data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The attributes have the following meaning:\n", "* **Survived**: that's the target, 0 means the passenger did not survive, while 1 means he/she survived.\n", "* **Pclass**: passenger class.\n", "* **Name**, **Sex**, **Age**: self-explanatory\n", "* **SibSp**: how many siblings & spouses of the passenger aboard the Titanic.\n", "* **Parch**: how many children & parents of the passenger aboard the Titanic.\n", "* **Ticket**: ticket id\n", "* **Fare**: price paid (in pounds)\n", "* **Cabin**: passenger's cabin number\n", "* **Embarked**: where the passenger embarked the Titanic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's get more info to see how much data is missing:" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [], "source": [ "train_data.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, the **Age**, **Cabin** and **Embarked** attributes are sometimes null (less than 891 non-null), especially the **Cabin** (77% are null). We will ignore the **Cabin** for now and focus on the rest. The **Age** attribute has about 19% null values, so we will need to decide what to do with them. Replacing null values with the median age seems reasonable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **Name** and **Ticket** attributes may have some value, but they will be a bit tricky to convert into useful numbers that a model can consume. So for now, we will ignore them." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at the numerical attributes:" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [], "source": [ "train_data.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Yikes, only 38% **Survived**. :( That's close enough to 40%, so accuracy will be a reasonable metric to evaluate our model.\n", "* The mean **Fare** was £32.20, which does not seem so expensive (but it was probably a lot of money back then).\n", "* The mean **Age** was less than 30 years old." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check that the target is indeed 0 or 1:" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [], "source": [ "train_data[\"Survived\"].value_counts()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's take a quick look at all the categorical attributes:" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [], "source": [ "train_data[\"Pclass\"].value_counts()" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [], "source": [ "train_data[\"Sex\"].value_counts()" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [], "source": [ "train_data[\"Embarked\"].value_counts()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Embarked attribute tells us where the passenger embarked: C=Cherbourg, Q=Queenstown, S=Southampton." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: the code below uses a mix of `Pipeline`, `FeatureUnion` and a custom `DataFrameSelector` to preprocess some columns differently. Since Scikit-Learn 0.20, it is preferable to use a `ColumnTransformer`, like in the previous chapter." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's build our preprocessing pipelines. We will reuse the `DataframeSelector` we built in the previous chapter to select specific attributes from the `DataFrame`:" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [], "source": [ "from sklearn.base import BaseEstimator, TransformerMixin\n", "\n", "class DataFrameSelector(BaseEstimator, TransformerMixin):\n", " def __init__(self, attribute_names):\n", " self.attribute_names = attribute_names\n", " def fit(self, X, y=None):\n", " return self\n", " def transform(self, X):\n", " return X[self.attribute_names]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's build the pipeline for the numerical attributes:" ] }, { "cell_type": "code", "execution_count": 113, "metadata": {}, "outputs": [], "source": [ "from sklearn.pipeline import Pipeline\n", "from sklearn.impute import SimpleImputer\n", "\n", "num_pipeline = Pipeline([\n", " (\"select_numeric\", DataFrameSelector([\"Age\", \"SibSp\", \"Parch\", \"Fare\"])),\n", " (\"imputer\", SimpleImputer(strategy=\"median\")),\n", " ])" ] }, { "cell_type": "code", "execution_count": 114, "metadata": {}, "outputs": [], "source": [ "num_pipeline.fit_transform(train_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also need an imputer for the string categorical columns (the regular `SimpleImputer` does not work on those):" ] }, { "cell_type": "code", "execution_count": 115, "metadata": {}, "outputs": [], "source": [ "# Inspired from stackoverflow.com/questions/25239958\n", "class MostFrequentImputer(BaseEstimator, TransformerMixin):\n", " def fit(self, X, y=None):\n", " self.most_frequent_ = pd.Series([X[c].value_counts().index[0] for c in X],\n", " index=X.columns)\n", " return self\n", " def transform(self, X, y=None):\n", " return X.fillna(self.most_frequent_)" ] }, { "cell_type": "code", "execution_count": 116, "metadata": {}, "outputs": [], "source": [ "from sklearn.preprocessing import OneHotEncoder" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can build the pipeline for the categorical attributes:" ] }, { "cell_type": "code", "execution_count": 117, "metadata": {}, "outputs": [], "source": [ "cat_pipeline = Pipeline([\n", " (\"select_cat\", DataFrameSelector([\"Pclass\", \"Sex\", \"Embarked\"])),\n", " (\"imputer\", MostFrequentImputer()),\n", " (\"cat_encoder\", OneHotEncoder(sparse=False)),\n", " ])" ] }, { "cell_type": "code", "execution_count": 118, "metadata": {}, "outputs": [], "source": [ "cat_pipeline.fit_transform(train_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let's join the numerical and categorical pipelines:" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [], "source": [ "from sklearn.pipeline import FeatureUnion\n", "preprocess_pipeline = FeatureUnion(transformer_list=[\n", " (\"num_pipeline\", num_pipeline),\n", " (\"cat_pipeline\", cat_pipeline),\n", " ])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cool! Now we have a nice preprocessing pipeline that takes the raw data and outputs numerical input features that we can feed to any Machine Learning model we want." ] }, { "cell_type": "code", "execution_count": 120, "metadata": {}, "outputs": [], "source": [ "X_train = preprocess_pipeline.fit_transform(train_data)\n", "X_train" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's not forget to get the labels:" ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [], "source": [ "y_train = train_data[\"Survived\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are now ready to train a classifier. Let's start with an `SVC`:" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [], "source": [ "from sklearn.svm import SVC\n", "\n", "svm_clf = SVC(gamma=\"auto\")\n", "svm_clf.fit(X_train, y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great, our model is trained, let's use it to make predictions on the test set:" ] }, { "cell_type": "code", "execution_count": 123, "metadata": {}, "outputs": [], "source": [ "X_test = preprocess_pipeline.transform(test_data)\n", "y_pred = svm_clf.predict(X_test)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we could just build a CSV file with these predictions (respecting the format excepted by Kaggle), then upload it and hope for the best. But wait! We can do better than hope. Why don't we use cross-validation to have an idea of how good our model is?" ] }, { "cell_type": "code", "execution_count": 124, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_val_score\n", "\n", "svm_scores = cross_val_score(svm_clf, X_train, y_train, cv=10)\n", "svm_scores.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, over 73% accuracy, clearly better than random chance, but it's not a great score. Looking at the [leaderboard](https://www.kaggle.com/c/titanic/leaderboard) for the Titanic competition on Kaggle, you can see that you need to reach above 80% accuracy to be within the top 10% Kagglers. Some reached 100%, but since you can easily find the [list of victims](https://www.encyclopedia-titanica.org/titanic-victims/) of the Titanic, it seems likely that there was little Machine Learning involved in their performance! ;-) So let's try to build a model that reaches 80% accuracy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try a `RandomForestClassifier`:" ] }, { "cell_type": "code", "execution_count": 125, "metadata": {}, "outputs": [], "source": [ "from sklearn.ensemble import RandomForestClassifier\n", "\n", "forest_clf = RandomForestClassifier(n_estimators=100, random_state=42)\n", "forest_scores = cross_val_score(forest_clf, X_train, y_train, cv=10)\n", "forest_scores.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's much better!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of just looking at the mean accuracy across the 10 cross-validation folds, let's plot all 10 scores for each model, along with a box plot highlighting the lower and upper quartiles, and \"whiskers\" showing the extent of the scores (thanks to Nevin Yilmaz for suggesting this visualization). Note that the `boxplot()` function detects outliers (called \"fliers\") and does not include them within the whiskers. Specifically, if the lower quartile is $Q_1$ and the upper quartile is $Q_3$, then the interquartile range $IQR = Q_3 - Q_1$ (this is the box's height), and any score lower than $Q_1 - 1.5 \\times IQR$ is a flier, and so is any score greater than $Q3 + 1.5 \\times IQR$." ] }, { "cell_type": "code", "execution_count": 126, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8, 4))\n", "plt.plot([1]*10, svm_scores, \".\")\n", "plt.plot([2]*10, forest_scores, \".\")\n", "plt.boxplot([svm_scores, forest_scores], labels=(\"SVM\",\"Random Forest\"))\n", "plt.ylabel(\"Accuracy\", fontsize=14)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To improve this result further, you could:\n", "* Compare many more models and tune hyperparameters using cross validation and grid search,\n", "* Do more feature engineering, for example:\n", " * replace **SibSp** and **Parch** with their sum,\n", " * try to identify parts of names that correlate well with the **Survived** attribute (e.g. if the name contains \"Countess\", then survival seems more likely),\n", "* try to convert numerical attributes to categorical attributes: for example, different age groups had very different survival rates (see below), so it may help to create an age bucket category and use it instead of the age. Similarly, it may be useful to have a special category for people traveling alone since only 30% of them survived (see below)." ] }, { "cell_type": "code", "execution_count": 127, "metadata": {}, "outputs": [], "source": [ "train_data[\"AgeBucket\"] = train_data[\"Age\"] // 15 * 15\n", "train_data[[\"AgeBucket\", \"Survived\"]].groupby(['AgeBucket']).mean()" ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [], "source": [ "train_data[\"RelativesOnboard\"] = train_data[\"SibSp\"] + train_data[\"Parch\"]\n", "train_data[[\"RelativesOnboard\", \"Survived\"]].groupby(['RelativesOnboard']).mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Spam classifier" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's fetch the data:" ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [], "source": [ "import os\n", "import tarfile\n", "from six.moves import urllib\n", "\n", "DOWNLOAD_ROOT = \"http://spamassassin.apache.org/old/publiccorpus/\"\n", "HAM_URL = DOWNLOAD_ROOT + \"20030228_easy_ham.tar.bz2\"\n", "SPAM_URL = DOWNLOAD_ROOT + \"20030228_spam.tar.bz2\"\n", "SPAM_PATH = os.path.join(\"datasets\", \"spam\")\n", "\n", "def fetch_spam_data(spam_url=SPAM_URL, spam_path=SPAM_PATH):\n", " if not os.path.isdir(spam_path):\n", " os.makedirs(spam_path)\n", " for filename, url in ((\"ham.tar.bz2\", HAM_URL), (\"spam.tar.bz2\", SPAM_URL)):\n", " path = os.path.join(spam_path, filename)\n", " if not os.path.isfile(path):\n", " urllib.request.urlretrieve(url, path)\n", " tar_bz2_file = tarfile.open(path)\n", " tar_bz2_file.extractall(path=SPAM_PATH)\n", " tar_bz2_file.close()" ] }, { "cell_type": "code", "execution_count": 130, "metadata": {}, "outputs": [], "source": [ "fetch_spam_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's load all the emails:" ] }, { "cell_type": "code", "execution_count": 131, "metadata": {}, "outputs": [], "source": [ "HAM_DIR = os.path.join(SPAM_PATH, \"easy_ham\")\n", "SPAM_DIR = os.path.join(SPAM_PATH, \"spam\")\n", "ham_filenames = [name for name in sorted(os.listdir(HAM_DIR)) if len(name) > 20]\n", "spam_filenames = [name for name in sorted(os.listdir(SPAM_DIR)) if len(name) > 20]" ] }, { "cell_type": "code", "execution_count": 132, "metadata": {}, "outputs": [], "source": [ "len(ham_filenames)" ] }, { "cell_type": "code", "execution_count": 133, "metadata": {}, "outputs": [], "source": [ "len(spam_filenames)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use Python's `email` module to parse these emails (this handles headers, encoding, and so on):" ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [], "source": [ "import email\n", "import email.policy\n", "\n", "def load_email(is_spam, filename, spam_path=SPAM_PATH):\n", " directory = \"spam\" if is_spam else \"easy_ham\"\n", " with open(os.path.join(spam_path, directory, filename), \"rb\") as f:\n", " return email.parser.BytesParser(policy=email.policy.default).parse(f)" ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [], "source": [ "ham_emails = [load_email(is_spam=False, filename=name) for name in ham_filenames]\n", "spam_emails = [load_email(is_spam=True, filename=name) for name in spam_filenames]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at one example of ham and one example of spam, to get a feel of what the data looks like:" ] }, { "cell_type": "code", "execution_count": 136, "metadata": {}, "outputs": [], "source": [ "print(ham_emails[1].get_content().strip())" ] }, { "cell_type": "code", "execution_count": 137, "metadata": {}, "outputs": [], "source": [ "print(spam_emails[6].get_content().strip())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some emails are actually multipart, with images and attachments (which can have their own attachments). Let's look at the various types of structures we have:" ] }, { "cell_type": "code", "execution_count": 138, "metadata": {}, "outputs": [], "source": [ "def get_email_structure(email):\n", " if isinstance(email, str):\n", " return email\n", " payload = email.get_payload()\n", " if isinstance(payload, list):\n", " return \"multipart({})\".format(\", \".join([\n", " get_email_structure(sub_email)\n", " for sub_email in payload\n", " ]))\n", " else:\n", " return email.get_content_type()" ] }, { "cell_type": "code", "execution_count": 139, "metadata": {}, "outputs": [], "source": [ "from collections import Counter\n", "\n", "def structures_counter(emails):\n", " structures = Counter()\n", " for email in emails:\n", " structure = get_email_structure(email)\n", " structures[structure] += 1\n", " return structures" ] }, { "cell_type": "code", "execution_count": 140, "metadata": {}, "outputs": [], "source": [ "structures_counter(ham_emails).most_common()" ] }, { "cell_type": "code", "execution_count": 141, "metadata": {}, "outputs": [], "source": [ "structures_counter(spam_emails).most_common()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It seems that the ham emails are more often plain text, while spam has quite a lot of HTML. Moreover, quite a few ham emails are signed using PGP, while no spam is. In short, it seems that the email structure is useful information to have." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's take a look at the email headers:" ] }, { "cell_type": "code", "execution_count": 142, "metadata": {}, "outputs": [], "source": [ "for header, value in spam_emails[0].items():\n", " print(header,\":\",value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There's probably a lot of useful information in there, such as the sender's email address (12a1mailbot1@web.de looks fishy), but we will just focus on the `Subject` header:" ] }, { "cell_type": "code", "execution_count": 143, "metadata": {}, "outputs": [], "source": [ "spam_emails[0][\"Subject\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, before we learn too much about the data, let's not forget to split it into a training set and a test set:" ] }, { "cell_type": "code", "execution_count": 144, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from sklearn.model_selection import train_test_split\n", "\n", "X = np.array(ham_emails + spam_emails)\n", "y = np.array([0] * len(ham_emails) + [1] * len(spam_emails))\n", "\n", "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, let's start writing the preprocessing functions. First, we will need a function to convert HTML to plain text. Arguably the best way to do this would be to use the great [BeautifulSoup](https://www.crummy.com/software/BeautifulSoup/) library, but I would like to avoid adding another dependency to this project, so let's hack a quick & dirty solution using regular expressions (at the risk of [un̨ho͞ly radiańcé destro҉ying all enli̍̈́̂̈́ghtenment](https://stackoverflow.com/a/1732454/38626)). The following function first drops the `` section, then converts all `` tags to the word HYPERLINK, then it gets rid of all HTML tags, leaving only the plain text. For readability, it also replaces multiple newlines with single newlines, and finally it unescapes html entities (such as `>` or ` `):" ] }, { "cell_type": "code", "execution_count": 145, "metadata": {}, "outputs": [], "source": [ "import re\n", "from html import unescape\n", "\n", "def html_to_plain_text(html):\n", " text = re.sub('.*?', '', html, flags=re.M | re.S | re.I)\n", " text = re.sub('', ' HYPERLINK ', text, flags=re.M | re.S | re.I)\n", " text = re.sub('<.*?>', '', text, flags=re.M | re.S)\n", " text = re.sub(r'(\\s*\\n)+', '\\n', text, flags=re.M | re.S)\n", " return unescape(text)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see if it works. This is HTML spam:" ] }, { "cell_type": "code", "execution_count": 146, "metadata": {}, "outputs": [], "source": [ "html_spam_emails = [email for email in X_train[y_train==1]\n", " if get_email_structure(email) == \"text/html\"]\n", "sample_html_spam = html_spam_emails[7]\n", "print(sample_html_spam.get_content().strip()[:1000], \"...\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And this is the resulting plain text:" ] }, { "cell_type": "code", "execution_count": 147, "metadata": {}, "outputs": [], "source": [ "print(html_to_plain_text(sample_html_spam.get_content())[:1000], \"...\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great! Now let's write a function that takes an email as input and returns its content as plain text, whatever its format is:" ] }, { "cell_type": "code", "execution_count": 148, "metadata": {}, "outputs": [], "source": [ "def email_to_text(email):\n", " html = None\n", " for part in email.walk():\n", " ctype = part.get_content_type()\n", " if not ctype in (\"text/plain\", \"text/html\"):\n", " continue\n", " try:\n", " content = part.get_content()\n", " except: # in case of encoding issues\n", " content = str(part.get_payload())\n", " if ctype == \"text/plain\":\n", " return content\n", " else:\n", " html = content\n", " if html:\n", " return html_to_plain_text(html)" ] }, { "cell_type": "code", "execution_count": 149, "metadata": {}, "outputs": [], "source": [ "print(email_to_text(sample_html_spam)[:100], \"...\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's throw in some stemming! For this to work, you need to install the Natural Language Toolkit ([NLTK](http://www.nltk.org/)). It's as simple as running the following command (don't forget to activate your virtualenv first; if you don't have one, you will likely need administrator rights, or use the `--user` option):\n", "\n", "`$ pip3 install nltk`" ] }, { "cell_type": "code", "execution_count": 150, "metadata": {}, "outputs": [], "source": [ "try:\n", " import nltk\n", "\n", " stemmer = nltk.PorterStemmer()\n", " for word in (\"Computations\", \"Computation\", \"Computing\", \"Computed\", \"Compute\", \"Compulsive\"):\n", " print(word, \"=>\", stemmer.stem(word))\n", "except ImportError:\n", " print(\"Error: stemming requires the NLTK module.\")\n", " stemmer = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also need a way to replace URLs with the word \"URL\". For this, we could use hard core [regular expressions](https://mathiasbynens.be/demo/url-regex) but we will just use the [urlextract](https://github.com/lipoja/URLExtract) library. You can install it with the following command (don't forget to activate your virtualenv first; if you don't have one, you will likely need administrator rights, or use the `--user` option):\n", "\n", "`$ pip3 install urlextract`" ] }, { "cell_type": "code", "execution_count": 151, "metadata": {}, "outputs": [], "source": [ "try:\n", " import urlextract # may require an Internet connection to download root domain names\n", " \n", " url_extractor = urlextract.URLExtract()\n", " print(url_extractor.find_urls(\"Will it detect github.com and https://youtu.be/7Pq-S557XQU?t=3m32s\"))\n", "except ImportError:\n", " print(\"Error: replacing URLs requires the urlextract module.\")\n", " url_extractor = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are ready to put all this together into a transformer that we will use to convert emails to word counters. Note that we split sentences into words using Python's `split()` method, which uses whitespaces for word boundaries. This works for many written languages, but not all. For example, Chinese and Japanese scripts generally don't use spaces between words, and Vietnamese often uses spaces even between syllables. It's okay in this exercise, because the dataset is (mostly) in English." ] }, { "cell_type": "code", "execution_count": 152, "metadata": {}, "outputs": [], "source": [ "from sklearn.base import BaseEstimator, TransformerMixin\n", "\n", "class EmailToWordCounterTransformer(BaseEstimator, TransformerMixin):\n", " def __init__(self, strip_headers=True, lower_case=True, remove_punctuation=True,\n", " replace_urls=True, replace_numbers=True, stemming=True):\n", " self.strip_headers = strip_headers\n", " self.lower_case = lower_case\n", " self.remove_punctuation = remove_punctuation\n", " self.replace_urls = replace_urls\n", " self.replace_numbers = replace_numbers\n", " self.stemming = stemming\n", " def fit(self, X, y=None):\n", " return self\n", " def transform(self, X, y=None):\n", " X_transformed = []\n", " for email in X:\n", " text = email_to_text(email) or \"\"\n", " if self.lower_case:\n", " text = text.lower()\n", " if self.replace_urls and url_extractor is not None:\n", " urls = list(set(url_extractor.find_urls(text)))\n", " urls.sort(key=lambda url: len(url), reverse=True)\n", " for url in urls:\n", " text = text.replace(url, \" URL \")\n", " if self.replace_numbers:\n", " text = re.sub(r'\\d+(?:\\.\\d*(?:[eE]\\d+))?', 'NUMBER', text)\n", " if self.remove_punctuation:\n", " text = re.sub(r'\\W+', ' ', text, flags=re.M)\n", " word_counts = Counter(text.split())\n", " if self.stemming and stemmer is not None:\n", " stemmed_word_counts = Counter()\n", " for word, count in word_counts.items():\n", " stemmed_word = stemmer.stem(word)\n", " stemmed_word_counts[stemmed_word] += count\n", " word_counts = stemmed_word_counts\n", " X_transformed.append(word_counts)\n", " return np.array(X_transformed)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try this transformer on a few emails:" ] }, { "cell_type": "code", "execution_count": 153, "metadata": {}, "outputs": [], "source": [ "X_few = X_train[:3]\n", "X_few_wordcounts = EmailToWordCounterTransformer().fit_transform(X_few)\n", "X_few_wordcounts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks about right!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have the word counts, and we need to convert them to vectors. For this, we will build another transformer whose `fit()` method will build the vocabulary (an ordered list of the most common words) and whose `transform()` method will use the vocabulary to convert word counts to vectors. The output is a sparse matrix." ] }, { "cell_type": "code", "execution_count": 154, "metadata": {}, "outputs": [], "source": [ "from scipy.sparse import csr_matrix\n", "\n", "class WordCounterToVectorTransformer(BaseEstimator, TransformerMixin):\n", " def __init__(self, vocabulary_size=1000):\n", " self.vocabulary_size = vocabulary_size\n", " def fit(self, X, y=None):\n", " total_count = Counter()\n", " for word_count in X:\n", " for word, count in word_count.items():\n", " total_count[word] += min(count, 10)\n", " most_common = total_count.most_common()[:self.vocabulary_size]\n", " self.most_common_ = most_common\n", " self.vocabulary_ = {word: index + 1 for index, (word, count) in enumerate(most_common)}\n", " return self\n", " def transform(self, X, y=None):\n", " rows = []\n", " cols = []\n", " data = []\n", " for row, word_count in enumerate(X):\n", " for word, count in word_count.items():\n", " rows.append(row)\n", " cols.append(self.vocabulary_.get(word, 0))\n", " data.append(count)\n", " return csr_matrix((data, (rows, cols)), shape=(len(X), self.vocabulary_size + 1))" ] }, { "cell_type": "code", "execution_count": 155, "metadata": {}, "outputs": [], "source": [ "vocab_transformer = WordCounterToVectorTransformer(vocabulary_size=10)\n", "X_few_vectors = vocab_transformer.fit_transform(X_few_wordcounts)\n", "X_few_vectors" ] }, { "cell_type": "code", "execution_count": 156, "metadata": {}, "outputs": [], "source": [ "X_few_vectors.toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What does this matrix mean? Well, the 64 in the third row, first column, means that the third email contains 64 words that are not part of the vocabulary. The 1 next to it means that the first word in the vocabulary is present once in this email. The 2 next to it means that the second word is present twice, and so on. You can look at the vocabulary to know which words we are talking about. The first word is \"of\", the second word is \"and\", etc." ] }, { "cell_type": "code", "execution_count": 157, "metadata": {}, "outputs": [], "source": [ "vocab_transformer.vocabulary_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are now ready to train our first spam classifier! Let's transform the whole dataset:" ] }, { "cell_type": "code", "execution_count": 158, "metadata": {}, "outputs": [], "source": [ "from sklearn.pipeline import Pipeline\n", "\n", "preprocess_pipeline = Pipeline([\n", " (\"email_to_wordcount\", EmailToWordCounterTransformer()),\n", " (\"wordcount_to_vector\", WordCounterToVectorTransformer()),\n", "])\n", "\n", "X_train_transformed = preprocess_pipeline.fit_transform(X_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: to be future-proof, we set `solver=\"lbfgs\"` since this will be the default value in Scikit-Learn 0.22." ] }, { "cell_type": "code", "execution_count": 159, "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LogisticRegression\n", "from sklearn.model_selection import cross_val_score\n", "\n", "log_clf = LogisticRegression(solver=\"lbfgs\", random_state=42)\n", "score = cross_val_score(log_clf, X_train_transformed, y_train, cv=3, verbose=3)\n", "score.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Over 98.7%, not bad for a first try! :) However, remember that we are using the \"easy\" dataset. You can try with the harder datasets, the results won't be so amazing. You would have to try multiple models, select the best ones and fine-tune them using cross-validation, and so on.\n", "\n", "But you get the picture, so let's stop now, and just print out the precision/recall we get on the test set:" ] }, { "cell_type": "code", "execution_count": 160, "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import precision_score, recall_score\n", "\n", "X_test_transformed = preprocess_pipeline.transform(X_test)\n", "\n", "log_clf = LogisticRegression(solver=\"lbfgs\", random_state=42)\n", "log_clf.fit(X_train_transformed, y_train)\n", "\n", "y_pred = log_clf.predict(X_test_transformed)\n", "\n", "print(\"Precision: {:.2f}%\".format(100 * precision_score(y_test, y_pred)))\n", "print(\"Recall: {:.2f}%\".format(100 * recall_score(y_test, y_pred)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.8" }, "nav_menu": {}, "toc": { "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 6, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 1 }