1 - Unsupervised Learning/Unsupervised Learning.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/80x15.png\" /></a>&nbsp;| [Emmanuel Rachelson](https://personnel.isae-supaero.fr/emmanuel-rachelson?lang=en) | <a href=\"https://supaerodatascience.github.io/machine-learning/\">https://supaerodatascience.github.io/machine-learning/</a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"font-size:22pt; line-height:25pt; font-weight:bold; text-align:center;\">Unsupervised Learning</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective of this notebook is to discover the beginning of the **Data Analysis Workflow** :\n",
    "1. Collect\n",
    "2. Analyze\n",
    "3. Deploy\n",
    "4. Decide\n",
    "\n",
    "\n",
    "At the same time we will discover **three Unsupervized Learning** tasks are illustrated here:\n",
    "1. [Dimensionality reduction](#dim)\n",
    "2. [Clustering](#clust)\n",
    "3. [Density estimation](#density)\n",
    "\n",
    " The goal is not to explain the algorithms in depth, but rather to understand the tasks they solve and see them in action. This is the occasion to familiarize ourselves with the scikit-learn API. \n",
    "\n",
    "As we have seen in the introduction to Machine Learning, Unsupervised Learning is concerned with learning *structure* from raw inputs, without the help of a *supervisor* who would help map examples to labels. There are three main categories of things one can learn without supervision:\n",
    "- A statistically relevant description of the data. That's the goal of *[Dimensionality reduction](#dim)*. It is intensively used for *data visualization* and exploration. It is strongly linked to previous classes in Statistics and the Introduction to Statistical Learning. It will be covered again in future classes using non-linear methods. \n",
    "- Categories of similar data. This supposes some implicit metric allowing to define distances between data samples. This is the goal of *[Clustering](#clust)*, which has already been covered in the Statistics class.\n",
    "- Generating distributions underlying the observed data. That's the goal of *[Density estimation](#density)*, which will be proposed as a \"going further\" exercice in the SVM notebook (next class)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"dim\"></a> 1. Dimensionality reduction"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img id=\"fig1\" src=\"https://imgs.xkcd.com/comics/real_estate.png\">\n",
    "\n",
    "A Real estate company comes to you : they have been **collecting data** for years, and know they think they can make some business with them.\n",
    "\n",
    "You start the **analysis** by loading and plotting them.\n",
    "\n",
    "**Disclamer** :\n",
    "The Boston housing prices dataset has an ethical problem. As an introduction to Usupervised learning, this course decides to keep using this dataset to educate about biais and ethical issues in data science.\n",
    "\n",
    "\n",
    "If you prefer to work on an other dataset, you can load an alternative dataset with :\n",
    "\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "boston = fetch_california_housing()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "data_url = \"http://lib.stat.cmu.edu/datasets/boston\"\n",
    "raw_df = pd.read_csv(data_url, sep=\"\\s+\", skiprows=22, header=None)\n",
    "X = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])\n",
    "Y = raw_df.values[1::2, 2]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions** :\n",
    "- What does the phrase \"Median Value (attribute 14) is usually the target\" means ?\n",
    "- Using one of the basic plot from the **bokeh** library [https://docs.bokeh.org/en/latest/docs/user_guide/plotting.html#userguide-plotting](here), create one plot showing two interesting features.  \n",
    "- What other feature would you have proposed to collect ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your plot here"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 13 variables in this data set. But only 506 data points. That's very little for a hypercube in 13-dimensional space. Probably the variables describing house prices by neighborhoods are correlated and the data does not really live in $\\mathbb{R}^{13}$, but rather on a manifold. Maybe the data could be more easily described with a few well-chosen variables. If only we had heard of Principal Component Analysis in previous classes!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions** : \n",
    "- In the following cells, what is the parameter \"explained_variance_ratio_\" ?\n",
    "- What is the difference between feature selection and feature engineering ?\n",
    "- After plotting the data using PCA, what kind of applications could you propose to the Real Estate company ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "print(X.shape)\n",
    "\n",
    "boston_pca = PCA()\n",
    "boston_pca.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "boston_pca.explained_variance_ratio_"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, apparently, 80% of the data's variance is already explained by the first principal component. And another 16% can be explained by the second principal component. Let's plot the histogram of explained variance ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.bar(range(X.shape[1]), boston_pca.explained_variance_ratio_, color=\"r\", align=\"center\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's confirm the amount of variance captured by these first two components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.sum(boston_pca.explained_variance_ratio_[:2])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, let's retrain the same PCA so that it provides a projection operator on the two first components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_pca = PCA(n_components=2)\n",
    "boston_pca.fit(X)\n",
    "X_proj = boston_pca.transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_proj.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at these two principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_pca.components_"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now plot the data projected onto these two principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_proj[:,0],X_proj[:,1]);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Going further with PCA in scikit-learn](https://scikit-learn.org/stable/modules/decomposition.html#principal-component-analysis-pca)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"clust\"></a> 2. Clustering"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have filtered out irrelevant dimensions in our data and started exploring it, it seems two types of houses appear. Let's try to group these houses by similarity in the projection space of the two first principal components. This is a task of clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "boston_kmeans = KMeans(n_clusters=2)\n",
    "boston_kmeans.fit(X_proj)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now plot the data by cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = boston_kmeans.predict(X_proj)\n",
    "plt.scatter(X_proj[:,0], X_proj[:,1], c=y_pred);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we could have done this clustering directly in the original $\\mathbb{R}^{13}$ description space. Let's run the k-means algorithm again, predict what samples correspond to what class and then display them in the plane spanned by two first principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_kmeans2 = KMeans(n_clusters=2)\n",
    "boston_kmeans2.fit(X)\n",
    "y_pred = boston_kmeans2.predict(X)\n",
    "plt.scatter(X_proj[:,0], X_proj[:,1], c=y_pred);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions** : \n",
    "- What is the difference between fit and predict ?\n",
    "- How many clusters is too many ?\n",
    "- Using sklearn [examples](http://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html), what would be an interesting alternative to kmeans (and test it) ?\n",
    "- How do you explain to the Real Estate company what are the characteristic of each cluster, with a bokeh plot ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your plot here, to explain your results to the Real Estate company"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See this [example](http://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html) for a great comparison of clustering methods.\n",
    "\n",
    "[Going further with k-means in scikit-learn](https://scikit-learn.org/stable/modules/clustering.html#k-means)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"density\"></a> 3. Density estimation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we would like to estimate the likelihood of seeing new houses similar to those of the first cluster. More specifically, we would like to visualize the probability density underlying this cluster. This is, yet again, another problem of Unsupervised Learning, called Density Estimation.\n",
    "\n",
    "This one has not been seen in previous classes and, for now, we will just admit that \"One Class SVMs\" are a magic class of algorithms that do precisely that: estimate probability densities from data. There are, of course, other density estimation algorithms. Scikit-learn provides [a few of these methods](https://scikit-learn.org/stable/modules/density.html), that are more generally related to the idea of estimating a probability distribution's cumulative distribution function (for example via a histogram)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's fit a One Class SVM model to match our first cluster's data distribution in the first two principal components space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = X_proj[y_pred==0,:]\n",
    "plt.scatter(X1[:,0], X1[:,1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import svm\n",
    "\n",
    "boston_ocsvm = svm.OneClassSVM(gamma = 1e-3)\n",
    "boston_ocsvm.fit(X1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the level set of this density function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = np.min(X1[:,0])\n",
    "xmax = np.max(X1[:,0])\n",
    "ymin = np.min(X1[:,1])\n",
    "ymax = np.max(X1[:,1])\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xmin, xmax, 500), np.linspace(ymin, ymax, 500))\n",
    "\n",
    "Z = boston_ocsvm.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
    "Z = Z.reshape(xx.shape)\n",
    "\n",
    "plt.contourf(xx, yy, Z, cmap=plt.cm.PuBu);\n",
    "#plt.scatter(X1[:,0], X1[:,1]);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also be curious and, since One Class SVMs are a rather versatile tool, we could try to fit the full data set (and not just the first cluster) and then display it on the plane spanned by the two first principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_ocsvm = svm.OneClassSVM(gamma = 1e-3)\n",
    "\n",
    "boston_ocsvm.fit(X_proj)\n",
    "\n",
    "xmin = np.min(X_proj[:,0])\n",
    "xmax = np.max(X_proj[:,0])\n",
    "ymin = np.min(X_proj[:,1])\n",
    "ymax = np.max(X_proj[:,1])\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xmin, xmax, 500), np.linspace(ymin, ymax, 500))\n",
    "\n",
    "Z = boston_ocsvm.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
    "Z = Z.reshape(xx.shape)\n",
    "\n",
    "plt.contourf(xx, yy, Z, cmap=plt.cm.PuBu);\n",
    "plt.scatter(X_proj[:,0], X_proj[:,1]);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Questions**:\n",
    "- Using the previous cells, what threshold would you choose to define outliers ?\n",
    "- Is there on best order to apply the preprocessing steps ? (PCA, clustering, outliers detection ) ?\n",
    "- Produce two files (on a relevant part of the data) : one called \"training_dataset.csv\" and the other one called \"testing_dataset.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your export code here"
   ]
  }
 ],
 "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.11.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}