smd_pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Feature Extraction, Feature Generation, Dimensional Reduction and PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Introduction\n",
    "\n",
    "### Format of the Lecture Material\n",
    "\n",
    "The more technical lectures will be composed of jupyter notebooks.\n",
    "\n",
    "<https://jupyter.org>\n",
    "\n",
    "Notebooks are a mixture of code, explanatory material, images, results, ...\n",
    "\n",
    "\n",
    "You can execute them yourself on your own machine and “play” with the code\n",
    "\n",
    "<img src=\"./ml/images/jupyter.png\" alt=\"Jupyter Logo\" style=\"width: 200px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### scikit-learn (sklearn)\n",
    "\n",
    "Many examples in the lecture will use the `scikit-learn` library.\n",
    "\n",
    "<http://scikit-learn.org/stable>\n",
    "\n",
    "`scikit-learn` is a library for data mining and “classical” machine learning.\n",
    "\n",
    "The user guide is an excellent resource of examples and mathematical background:\n",
    "<https://scikit-learn.org/stable/user_guide.html>\n",
    "\n",
    "<img src=\"./ml/images/logo.png\" alt=\"Scikit Logo\" style=\"width: 200px;\"/>\n",
    "\n",
    "\n",
    "> Pedregosa, Fabian, et al. \"Scikit-learn: Machine learning in Python.\" the Journal of machine Learning research 12 (2011): 2825-2830."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### pandas\n",
    "\n",
    "To store our data in memory, for preprocessing and simple queries we will mainly use `pandas`.\n",
    "\n",
    "`pandas`' main feature is the `DataFrame` class, a 2d-table object providing many methods for data analysis, visualization, aggregation, ...\n",
    "\n",
    "<https://pandas.pydata.org>\n",
    "\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/e/ed/Pandas_logo.svg\" alt=\"pandas log\" style=\"width: 300px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "These lectures again come with a `conda` environment for all needed packages:\n",
    " \n",
    "Create the environment:\n",
    "```\n",
    "mamba env create -f environment.yml\n",
    "```\n",
    "\n",
    "\n",
    "Update:\n",
    "```\n",
    "mamba env update -f environment.yml\n",
    "```   \n",
    "\n",
    "Activate:\n",
    "```\n",
    "conda activate ml\n",
    "```\n",
    "\n",
    "Start the jupyter lab server:\n",
    "```\n",
    "jupyter lab\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Notation\n",
    "\n",
    "*We try to follow (with some extensions) the notations in  \"Elements of statistical Learning\" by Trevor Hastie et al. https://web.stanford.edu/~hastie/ElemStatLearn/ (Free E-Book)* \n",
    "\n",
    "This means:\n",
    "\n",
    "* Capital letters like $X$ or $Y$ denote a generic random variate\n",
    "* Observations/realizations are small letters, the $i$-th observation of $X$ is $x_i$\n",
    "* Matrices and vectors are capital, bold-face symbols $\\boldsymbol{X}$\n",
    "* Observations/realizations are *rows* of the matrix while different variables are stored in the *columns*\n",
    "\n",
    "Example: If we measure $d=2$ variables – e.g. age and weight – of $N = 100$ people, we get a $N \\times d$-dimensional matrix $\\boldsymbol{X}$.\n",
    "\n",
    "A single observation – one row – is written as $\\boldsymbol{x}_i = (\\mathrm{age}_i, \\mathrm{weight}_i)$.\n",
    "\n",
    "All observations of the *weight* variable are denoted $\\boldsymbol{x}_{\\bullet 1}$, similar to the `numpy` indexing `X[:, 1]` (this is an extension to Hasties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:23.953831Z",
     "start_time": "2018-11-13T16:47:22.205726Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from ml import plots\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:23.953831Z",
     "start_time": "2018-11-13T16:47:22.205726Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "plots.set_plot_style()\n",
    "\n",
    "%matplotlib widget\n",
    "\n",
    "colors = plots.colors\n",
    "cmap = plots.cmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Repetition\n",
    "\n",
    "Last lecture: Linear Fisher Discriminant:\n",
    "\n",
    "> Find the hyper plane optimally separating two populations using the Fisher criterion.\n",
    "\n",
    "Now, we can have a look at a simple example using `sklearn`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:24.481797Z",
     "start_time": "2018-11-13T16:47:23.956854Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "\n",
    "X, y = make_blobs(n_samples=250, centers=2, cluster_std=2.0, random_state=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:24.481797Z",
     "start_time": "2018-11-13T16:47:23.956854Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.set_aspect(1)\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:24.481797Z",
     "start_time": "2018-11-13T16:47:23.956854Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "\n",
    "clf = LinearDiscriminantAnalysis()\n",
    "lda = clf.fit(X, y)\n",
    "\n",
    "projection = lda.transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:24.481797Z",
     "start_time": "2018-11-13T16:47:23.956854Z"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax, ax_proj) = plt.subplots(1, 2)\n",
    "\n",
    "ax.set_aspect(1, 'datalim')\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap=cmap)\n",
    "\n",
    "plots.draw_linear_regression_function(lda, color='gray', ax=ax)\n",
    "\n",
    "\n",
    "ax.set_xlabel('$X_{0}$')\n",
    "ax.set_ylabel('$X_{1}$')\n",
    "\n",
    "for label, color in zip((0, 1), cmap.colors):\n",
    "    ax_proj.hist(projection[y == label], bins=20, range=[-5, 5], color=color, histtype='step')\n",
    "\n",
    "    \n",
    "ax_proj.axvline(0, color='gray')\n",
    "ax_proj.set_xlabel('Projection')\n",
    "    \n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Problems with high-dimensional data\n",
    "\n",
    "> __Curse of dimensionality__\n",
    ">\n",
    "> is a term, intrdocued by Richard Bellman,   \n",
    "> to describe the rapid increase of the volume when adding more dimensions in a mathematical space.\n",
    ">\n",
    "> <https://en.wikipedia.org/wiki/Curse_of_dimensionality> \n",
    "\n",
    "The higher the dimensionality of the space, the more observations are needed to *sufficiently* cover its volume.\n",
    "\n",
    "In the following, we will draw 100 samples from a standard uniform distribution, first in 1D, then in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:25.073068Z",
     "start_time": "2018-11-13T16:47:24.486586Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:25.073068Z",
     "start_time": "2018-11-13T16:47:24.486586Z"
    }
   },
   "outputs": [],
   "source": [
    "sample = rng.uniform(low=0.0, high=1.0, size=100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "hist, edges, plot = ax.hist(sample, range=[0, 1], bins=10, histtype='step', lw=2)\n",
    "\n",
    "ax.scatter(sample, np.full_like(sample, 1.1 * hist.max()), color='C0')\n",
    "density = np.count_nonzero(hist) / hist.size\n",
    "\n",
    "ax.set_ylabel('Number of observations')\n",
    "ax.set_title('Non-empty Bins: {:.2%}'.format(density))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:25.897248Z",
     "start_time": "2018-11-13T16:47:25.078702Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# uniform random numbers in two dimensions\n",
    "sample = rng.uniform(low=0.0, high=1.0, size=(100, 2))\n",
    "\n",
    "# a detailed plotting example. In future, we will focus more on the data processing and machine learning\n",
    "# and many of the plots will be defined as functions in the `ml/plots.py` module, not in the notebook itself\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)\n",
    "\n",
    "\n",
    "# Einzelne Punkte plotten\n",
    "ax1.scatter(sample[:, 0], sample[:, 1])\n",
    "ax1.set_aspect('equal')\n",
    "\n",
    "\n",
    "# das Histogram plotten\n",
    "hist, _, _, plot = ax2.hist2d(\n",
    "    sample[:, 0],\n",
    "    sample[:, 1],\n",
    "    range=[[0, 1], [0, 1]],\n",
    "    bins=10,\n",
    "    cmap='inferno',\n",
    "    vmin=0,\n",
    ")\n",
    "\n",
    "# Anteil besetzter Bins bestimmen\n",
    "density = np.count_nonzero(hist) / hist.size\n",
    "\n",
    "ax2.set_title('Non-empty Bins: {:.0%}'.format(density))\n",
    "ax2.set_aspect('equal')\n",
    "\n",
    "fig.colorbar(plot, ax=ax2)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the 1D example, all bins are non-empty. There are no empty bins.\n",
    "\n",
    "In the 2D case, roughly a third of the bins are empty.\n",
    "\n",
    "\n",
    "Possible mitigation approaches:\n",
    "\n",
    "1. Collect more data and live with the additional cost / time.\n",
    "2. Use larger bins, reduces the detail in each sub-dimension.\n",
    "3. Reduce the number of dimensions, possibly by combining existing dimensions into new ones. \n",
    "\n",
    "Just storing more data is almost always impossible, even today.\n",
    "\n",
    "\n",
    "### Example: IceCube\n",
    "\n",
    "* The IceCube Neutrino Observatory in Antarctica collects about 1 TB of raw data each day.\n",
    "* The satellite uplink is limited to about 100 GB per day.  \n",
    "* The remaing data is shipped once a year on magnetic tapes by container ship.\n",
    "\n",
    "<img src=\"./ml/images/icecube.jpg\" alt=\"Drawing\" style=\"width: 400px;\"/>\n",
    "\n",
    "> Never underestimate the data rate of a container ship full of magnetic tapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# once per year 365 TB is shipped, the ship takes roughly 30 days. \n",
    "# Data rate in MBit/s:\n",
    "365 * 8 * 1024**2 / (30 * 86400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's more than a 1GBit/s connection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "\n",
    "### Example: SKA\n",
    "\n",
    "* The Square Kilometer Array is a planned radio observatory, planned in South Africa and Australia.\n",
    "\n",
    "* It will use tens of thousands of radio antennas.\n",
    "\n",
    "* The expected raw data rate is expected to be multiple  __petabytes per second__.\n",
    "\n",
    "* Storage of the raw data is completely impossible given today's and the near future's technology."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:25.911324Z",
     "start_time": "2018-11-13T16:47:25.900962Z"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "\n",
    "IFrame('https://www.youtube.com/embed/8BBoDw2qVD0?rel=0', width=960, height=540)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Highly-dimensional data also lead to very fundamental, mathematical problems.\n",
    "\n",
    "See this interesting discussians about distance measures:\n",
    "\n",
    "> Why is Euclidean distance not a good metric in high dimensions?  \n",
    ">https://stats.stackexchange.com/questions/99171/why-is-euclidean-distance-not-a-good-metric-in-high-dimensions\n",
    "\n",
    "⇒ We need some form of data reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Data Reduction\n",
    "\n",
    "Two orthogonal, complementarty approaches:\n",
    "\n",
    "1. Feature Extraction\n",
    "    * Hand-crafted new features from higher-dimensional input data:\n",
    "        * This is a task commonly performed by domain experts to extract as much information as possible\n",
    "          from high dimensional data into fewer, more descriptive features \n",
    "    * General algorithmic reduction of dimensions, e.g. to maximize variance (principal component analysis)\n",
    "   \n",
    "1. Feature Selection\n",
    "    * Remove non-descriptive features \n",
    "  \n",
    "\n",
    "\n",
    "The general approach is:\n",
    "\n",
    "* Extract lots of hand-crafted features from the highly-dimensional, often irregular raw data\n",
    "* Perform a feature selection relevant to the given task to further reduce the number\n",
    "\n",
    "\n",
    "This usually happens in a multi-step procedure, creating higher and higher level abstractions from the original raw data.\n",
    "\n",
    "Example from Imaging Air Cherenkov Telescopes (IACTs):\n",
    "\n",
    "* Raw data consists of a charge time-series for each pixel of each telescope camera, possibly with multiple amplification factors (gains).  \n",
    "  For an array of telescopes, not all telescopes participate in each event.    \n",
    "  For each telescope, we have data of the shape $(N_\\mathrm{gains}, N_\\mathrm{pixels}, N_\\mathrm{samples})$  \n",
    "  \n",
    "* Reduce the waveforms by finding the pulses of Cherenkov signal, reducing the pulse to two numbers: \n",
    "    * an amplitude (integral or maximum) converted to the estimated number of Cherenkov photons\n",
    "    * a time (rising edge or mean arrival time)\n",
    "  This gives us two \"images\", one \"brightness\" and one \"time\" image\n",
    "  \n",
    "* Describe these images with a number of features (includes a Principal Component Analysis)\n",
    "\n",
    "* Combine these features per telescope to geometrically estimate the geometry of the shower axis (direction and impact point)\n",
    "\n",
    "* Use this as input for the machine learning for energy estimation and particle type classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Feature Extraction\n",
    "\n",
    "#### Example 1: Simple Data Transformation\n",
    "\n",
    "Given data points $X = (x_1, x_2, ...) , Y = (y_1, y_2, ...) $.\n",
    "\n",
    "Assume it is known or can easily be seen that these coordinates are more naturally expressed using polar coordinates:\n",
    "\n",
    "$$\n",
    " \\begin{align}\n",
    "     r &= x^2 + y^2 \\\\\n",
    "     \\phi &= \\operatorname{arctan2}(y, x)    \n",
    " \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:27.177886Z",
     "start_time": "2018-11-13T16:47:25.917299Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_circles\n",
    "\n",
    "def transform(X):\n",
    "    r = np.sqrt(X[:, 0]**2 + X[:, 1]**2)\n",
    "    phi = np.arctan2(X[:, 1], X[:, 0])\n",
    "    return np.column_stack([r, phi])\n",
    "\n",
    "\n",
    "X_original, y = make_circles(n_samples=500, noise=0.15, factor=0.4, )\n",
    "X_transformed = transform(X_original)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:27.177886Z",
     "start_time": "2018-11-13T16:47:25.917299Z"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "\n",
    "\n",
    "for X, ax, title in zip((X_original, X_transformed), axs, ('Original', 'Transformed')):\n",
    "    \n",
    "    ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap)\n",
    "    ax.set_xlabel('X1')\n",
    "    ax.set_ylabel('X2')\n",
    "    ax.set_title(title)\n",
    "\n",
    "axs[0].set_aspect(1, adjustable='datalim')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In this example, all information is contained in the new $r$ attribute, we can remove $\\phi$ reducing the number of dimensions from 2 to 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Example 2: Data representation\n",
    "\n",
    "Often, data is not in a suitable representation for applying the machine learning methods we learn about here.\n",
    "\n",
    "In general, everything has to be transformed into plain numbers and we need to arrive at the 2-d regular matrix $X$ of shape $(N_\\mathrm{observations}, N_\\mathrm{features})$.\n",
    "\n",
    "In this example, we look at *text* data from internet forums.\n",
    "\n",
    "We somehow have to transform this textual data into numerical features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:31.100707Z",
     "start_time": "2018-11-13T16:47:27.181136Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_20newsgroups\n",
    "\n",
    "atheist_texts = fetch_20newsgroups(remove=('headers', 'footers', 'quotes'), categories=['alt.atheism'])\n",
    "religous_texts = fetch_20newsgroups(remove=('headers', 'footers', 'quotes'), categories=['talk.religion.misc'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:31.100707Z",
     "start_time": "2018-11-13T16:47:27.181136Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:31.100707Z",
     "start_time": "2018-11-13T16:47:27.181136Z"
    }
   },
   "outputs": [],
   "source": [
    "print('** Random sample from the atheist forum: **\\n')\n",
    "print(rng.choice(atheist_texts.data).lstrip())\n",
    "print('\\n'*3)\n",
    "\n",
    "print('** Random sample from the religon forum: **\\n')\n",
    "print(rng.choice(religous_texts.data).lstrip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Hypothesis 1__: Atheists use more words\n",
    "\n",
    "Let's look at the length of the texts and see if we can separate the two populations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:31.680634Z",
     "start_time": "2018-11-13T16:47:31.103985Z"
    }
   },
   "outputs": [],
   "source": [
    "def number_of_words(text: str):\n",
    "    return len(text.split())\n",
    "\n",
    "atheist_lengths = list(map(number_of_words, atheist_texts.data))\n",
    "religous_lengths = list(map(number_of_words, religous_texts.data))\n",
    "\n",
    "\n",
    "bins =  np.arange(0, 2000, 25)\n",
    "\n",
    "hist_options = dict(\n",
    "    bins=bins,\n",
    "    alpha=0.5,\n",
    "    density=True,\n",
    "    histtype='step',\n",
    "    lw=4,\n",
    ")\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(atheist_lengths,  label='Atheist', **hist_options)\n",
    "plt.hist(religous_lengths, label='Religous', **hist_options)\n",
    "plt.xlabel('Text Length in Words')\n",
    "plt.legend()\n",
    "None\n",
    "\n",
    "print(f'Median      : {np.median(atheist_lengths):3.0f}, {np.median(religous_lengths):3.0f}')\n",
    "print(f'75%-Quantile: {np.percentile(atheist_lengths, 75):3.0f}, {np.percentile(religous_lengths, 75):3.0f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No significant difference visible.\n",
    "\n",
    "__Hypothesis 2__: Atheists use longer words\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:32.289809Z",
     "start_time": "2018-11-13T16:47:31.683907Z"
    }
   },
   "outputs": [],
   "source": [
    "def max_word_length(s):\n",
    "    words = s.split()\n",
    "    if not words:\n",
    "        return 0\n",
    "    return max(map(len, words))\n",
    "\n",
    "atheist_lengths = list(map(max_word_length, atheist_texts.data))\n",
    "religous_lengths = list(map(max_word_length, religous_texts.data))\n",
    "\n",
    "\n",
    "bins =  np.arange(0, 100, 1)\n",
    "hist_options = dict(\n",
    "    bins=bins,\n",
    "    alpha=0.5,\n",
    "    density=True,\n",
    "    histtype='step',\n",
    "    lw=4,\n",
    ")\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(atheist_lengths, label='Atheists', **hist_options)\n",
    "plt.hist(religous_lengths, label='Religous', **hist_options)\n",
    "plt.legend()\n",
    "None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where do this enourmous word length come from? It's not German...\n",
    "\n",
    "We are on to an important point: data preprocessing and \"cleaning\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:32.306869Z",
     "start_time": "2018-11-13T16:47:32.294360Z"
    }
   },
   "outputs": [],
   "source": [
    "atheist_texts.data[2].split()[40:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Hypothese 4__: Atheist and religous people use different words."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "string.punctuation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "remove_punctuation = str.maketrans({c: \" \" for c in string.punctuation})\n",
    "\"-----------------------------------------------------------------------------\".translate(remove_punctuation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:33.752517Z",
     "start_time": "2018-11-13T16:47:32.311675Z"
    }
   },
   "outputs": [],
   "source": [
    "import string\n",
    "from collections import Counter\n",
    "\n",
    "\n",
    "def extract_words(texts):\n",
    "    return (\n",
    "        ' '.join(texts.data)           # join all texts to one giant string\n",
    "        .lower()                       # all to lower case\n",
    "        .translate(remove_punctuation) # remove punctuation ,.- etc.\n",
    "        .split()                       # split sring into list of words at whitespace\n",
    "    )\n",
    "\n",
    "\n",
    "def most_common(words, n, min_length=5):\n",
    "    counter = Counter(filter(lambda w: len(w) > min_length, words))\n",
    "    return dict(counter.most_common(n))\n",
    "\n",
    "\n",
    "atheist_words = extract_words(atheist_texts)\n",
    "common_atheist_words = most_common(atheist_words, 25)\n",
    "\n",
    "religous_words = extract_words(religous_texts)\n",
    "common_religous_words = most_common(religous_words, 25)\n",
    "\n",
    "\n",
    "f, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "ax1.set_title('Atheists')\n",
    "ax1.barh(list(common_atheist_words.keys()), list(common_atheist_words.values()))\n",
    "\n",
    "ax2.set_title('Religous')\n",
    "ax2.barh(list(common_religous_words.keys()), list(common_religous_words.values()), color='C1')\n",
    "\n",
    "\n",
    "atheist_set = set(common_atheist_words.keys())       \n",
    "religous_set = set(common_religous_words.keys())\n",
    "s = atheist_set | religous_set\n",
    "print(\"Words atheist use more often:\", *(atheist_set - religous_set))\n",
    "print(\"Words religous people use more often:\", *(religous_set - atheist_set))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This kind of data reduction\n",
    "\n",
    "* requires (the obtaining) of large amounts of expert domain knowledge\n",
    "* gets harder and harder as data volume and dimensionality rise\n",
    "* is a very tedious process\n",
    "* usually obtained the best results... until relatively recently"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Principal Component Analysis (PCA)\n",
    "\n",
    "\n",
    "PCA searches for a basis in the data space that maximizes the variance along the basis vectors.\n",
    "\n",
    "Given $N$ data points in $d$ dimensions that shall be transformed to $k \\leq d$  dimensions.\n",
    "\n",
    "We do this by defining a new basis, combining the original features.\n",
    "\n",
    "High-Level approach:\n",
    "\n",
    "0. Centralize data points with respect to their mean.\n",
    "1. Calculate the covariance matrix $\\mathrm{Cov}(X)$ of the data matrix $\\boldsymbol{X}$\n",
    "2. Compute Eigenvalues und Eigenvectors of the covarinace matrix\n",
    "3. Choose the $k$ largest Eigenvalues and the corresponding Eigenvectors. \n",
    "4. Fill the $d \\times k$ matrix $\\boldsymbol{W}$ with the $k$ Eigenvectors.\n",
    "5. Apply $\\boldsymbol{W}$ to each row $x_i$ of $\\boldsymbol{X}$ ⇒ $x^\\prime = \\boldsymbol{W}^T \\cdot x^T $ \n",
    "    \n",
    "\n",
    "##### 1. Centralization \n",
    "\n",
    "Compute mean vector $\\mu$:\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\mu} = \\begin{pmatrix}\n",
    "    \\bar{\\boldsymbol{x}}_1 \\\\\n",
    "    \\ldots \\\\\n",
    "    \\bar{\\boldsymbol{x}}_d \\\\\n",
    "\\end{pmatrix}\n",
    " = \\frac 1 N\n",
    " \\begin{pmatrix}\n",
    "    \\sum_{i=0}^{N} \\boldsymbol{x}_{1, i} \\\\\n",
    "    \\ldots \\\\\n",
    "    \\sum_{i=0}^{N} \\boldsymbol{x}_{d, i} \\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Alternative formulation:\n",
    "$$\n",
    "\\boldsymbol{\\mu} = \\begin{pmatrix}\n",
    "    \\text{Mean over all observations of attribute 1}\\\\\n",
    "    \\ldots \\\\\n",
    "    \\text{Mean over all observations of attribute 1}\\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "New data point:\n",
    "\n",
    "   $$\n",
    "      x^{\\prime}_i = x_i - \\boldsymbol{\\mu}\n",
    "   $$\n",
    "\n",
    "##### 2. Covariance\n",
    "\n",
    "The covariance matrix of a random variable $X$ in arbitrary dimensions:\n",
    "$$\n",
    "\\operatorname {Cov} (X)=\\operatorname {E} {\\bigl [}(X-\\operatorname {E} (X))\\cdot (X-\\operatorname {E} (X))^T{\\bigr ]}\n",
    "$$\n",
    "\n",
    "Estimation of the covariance matrix is possible through *simple* matrix operations.\n",
    "\n",
    "##### 3. Eigenvalues und Vectors\n",
    "\n",
    "Compute the $d$ different Eigenvalues of $ \\operatorname {Cov} (\\boldsymbol {X} )$.\n",
    "\n",
    "Obtain Eigenvalues $\\lambda_1, \\ldots, \\lambda_d$ with corresponding Eigenvectors $v_1, \\ldots, v_d$\n",
    "\n",
    "\n",
    "##### 4. Sorting and Choosing a Subset\n",
    "\n",
    "Sort indices of eigenvalues and -vectors so that:\n",
    "$$\n",
    "\\lambda_1 > \\lambda_2 > \\lambda_3 \\ldots > \\lambda_d\n",
    "$$\n",
    "\n",
    "Choose the $k$ largest eigenvalues and discard the rest.\n",
    "\n",
    "\n",
    "##### 5. Create the Transformation Matrix\n",
    "\n",
    "Use the $k$ chosen eigenvectors as columns of the matrix $\\boldsymbol{W}$\n",
    "\n",
    "$$\n",
    "\\boldsymbol{W} = \\begin{pmatrix}\n",
    "    v_1, \n",
    "    \\ldots,  \n",
    "    v_k\n",
    "\\end{pmatrix}\n",
    "= \\begin{pmatrix}\n",
    "    v_{1,1}, \n",
    "    \\ldots,  \n",
    "    v_{k, 1} \\\\\n",
    "    \\ldots \\\\\n",
    "        v_{1,d}, \n",
    "    \\ldots,  \n",
    "    v_{k, d}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "##### 6. Transformations\n",
    "\n",
    "Multiply the transformation matrix $\\boldsymbol{W}$ to each observatoin $x_i$ in $\\boldsymbol{X}$ to obtain the new observations reduced to $k$ dimensions:\n",
    "\n",
    "$$\n",
    "\\boldsymbol{X^\\prime} = \\boldsymbol{X} \\boldsymbol{W}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:33.911709Z",
     "start_time": "2018-11-13T16:47:33.755670Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "\n",
    "N = 8 # Number of observations\n",
    "d = 4 # Number of dimensions\n",
    "k = 2 # Number of dimensions after reduction\n",
    "\n",
    "X = rng.normal(size=(N, d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:33.911709Z",
     "start_time": "2018-11-13T16:47:33.755670Z"
    }
   },
   "outputs": [],
   "source": [
    "# centralize\n",
    "X = X - X.mean(axis=0)\n",
    "\n",
    "# compute covariance\n",
    "c = np.cov(X, rowvar=False)\n",
    "\n",
    "# compute eigenvalues and eigenvectors\n",
    "# eigh already sorts the eigenvalues, but ascending (smallest first)\n",
    "l, W = np.linalg.eigh(c)\n",
    "\n",
    "# Invert order\n",
    "l = l[::-1]\n",
    "W = W[:, ::-1]\n",
    "\n",
    "# Select first k\n",
    "l = l[:k]\n",
    "W = W[:, :k]\n",
    "\n",
    "X_prime = X @ W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, X_prime.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:33.911709Z",
     "start_time": "2018-11-13T16:47:33.755670Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=k)\n",
    "X_prime_sklearn = pca.fit_transform(X)\n",
    "\n",
    "\n",
    "with np.printoptions(precision=2):\n",
    "    print('By Hand:\\n', X_prime)\n",
    "    print('sklearn:\\n', X_prime_sklearn)\n",
    "\n",
    "# testen aller einträge auf gleichheit (bis auf vorzeichen)\n",
    "print('\\n All close:', np.allclose(np.abs(X_prime), np.abs(X_prime_sklearn)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Example in 3D\n",
    "\n",
    "An artifical data set using $d = 3$ Dimensionen is reduced to $k=2$ dimensions.   \n",
    "\n",
    "The data set is composed of two normally-distributed populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:35.613041Z",
     "start_time": "2018-11-13T16:47:33.914685Z"
    }
   },
   "outputs": [],
   "source": [
    "X, y = make_blobs(n_samples=500, n_features=3, cluster_std=3, random_state=2)\n",
    "plots.plot_3d_views(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:37.396225Z",
     "start_time": "2018-11-13T16:47:35.616379Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=3)\n",
    "transformed = pca.fit_transform(X)\n",
    "\n",
    "plots.plot_3d_views(transformed, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:38.050103Z",
     "start_time": "2018-11-13T16:47:37.399569Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "transformed = pca.fit_transform(X)\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(transformed[:, 0], transformed[:, 1], c=y, cmap=cmap)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Facial Recognition \n",
    "\n",
    "Given each pixel of an image as one feature, where the value is a gray-scale value between 0 (black) and 1 (white).\n",
    "\n",
    "\n",
    "Each image is stored as a 1d vector, row after row.\n",
    "Multiple images give us the matrix $\\boldsymbol{X}$\n",
    "\n",
    "If our images have a resolution of $64 \\times 64$ pixels, we get a vector of length $4096$.\n",
    "\n",
    "##### Simple Facial Recognition \n",
    "\n",
    "Assume the task is the assigment of photos of all students of TU Dortmund University to their names.\n",
    "\n",
    "\n",
    "We look for a function that transforms an image into a name.\n",
    "\n",
    "\n",
    "Idea:\n",
    "1. Store images of all students in a matrix $\\boldsymbol{X}$ of dimensionality $N_\\text{studens} \\times N_\\text{pixels}$ and a vector of labels of length $N_\\text{students}$ containing the names (or matriculation ID)\n",
    "2. Compute the distance $D$ between a new photo $x_{\\text{neu}}$ to all photos stored in $\\boldsymbol{X}$\n",
    "3. Return $y_i$ corresponding to the $i$ where $D(x_{\\text{neu}}, x_i)$ is minimal.\n",
    "\n",
    "Problems: \n",
    " - Storing so many images in memory is difficult to impossible. \n",
    " - Computing the distance will take a long time.\n",
    " - As discussed above, chosing the right distance measure for this very high dimensional data set might be difficult.\n",
    " \n",
    "##### Eigenfaces \n",
    "\n",
    "*Original Article from 1991 by Turk and Pentland http://www.mitpressjournals.org/doi/10.1162/jocn.1991.3.1.71*\n",
    "\n",
    "The input stays the same, the matrix of all photos $\\boldsymbol{X}$. However, this time, we don't store the full image information.\n",
    "\n",
    "Idea:\n",
    "1. Apply a PCA to $\\boldsymbol{X}$. \n",
    "2. Obtain transformation matrix $\\boldsymbol{W}$ of shape $d \\times k$\n",
    "3. Calculate weights $g_m = \\boldsymbol{v}_m^T \\cdot(x_i - \\boldsymbol{\\mu}) $ for each iamge $x_i$ and each eigenvector $\\boldsymbol{v}_m$ with $m \\in \\{1, \\ldots, k\\}$ and obtain an weight vector $G$ of length $k$.\n",
    "4. Compute distance $D$ between weight vectors of all \"training images\" and the new image $G_{\\text{new}}$\n",
    "5. Return  $y_i$ with $i$ for which $D(G_{\\text{neu}}, G^{i})$ is minimal.\n",
    "\n",
    "In reality, computing a PCA on very large, high-dimensional data sets is not as trivial as it seems above.\n",
    "\n",
    "##### Python Example for Eigenfaces\n",
    "\n",
    "The LFW (Labeled Faces in the Wild) data set is a common data set to benchmark algorithms for facial recognition:   \n",
    "http://vis-www.cs.umass.edu/lfw/\n",
    "\n",
    "It comprises ca. 13,000 images of famous (as of ca. 2003) people.\n",
    "We choose a subset here that requires a minimum number of images for each person."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:39.027898Z",
     "start_time": "2018-11-13T16:47:38.055271Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_lfw_people\n",
    "\n",
    "lfw_people = fetch_lfw_people(min_faces_per_person=20, resize=0.8)\n",
    "\n",
    "h, w = lfw_people.images[0].shape\n",
    "\n",
    "X = lfw_people.data\n",
    "y = lfw_people.target\n",
    "\n",
    "# remove border pixels\n",
    "border = 10\n",
    "width = w - 2 * border\n",
    "height = h - 2 * border\n",
    "X = X.reshape((-1, h, w))[:, border:-border, border:-border].reshape((-1, height * width))\n",
    "\n",
    "\n",
    "names = lfw_people.target_names\n",
    "\n",
    "X.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(lfw_people.target_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:39.027898Z",
     "start_time": "2018-11-13T16:47:38.055271Z"
    }
   },
   "outputs": [],
   "source": [
    "index = rng.choice(len(lfw_people.images))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plot = ax.imshow(X[index].reshape(height, width), cmap='gray')\n",
    "title = ax.set_title(rf'Image {index} in $\\mathbf{{\\mathit{{X}}}}$ ({y[index]} = {names[y[index]]})')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:39.598180Z",
     "start_time": "2018-11-13T16:47:39.031942Z"
    }
   },
   "outputs": [],
   "source": [
    "mean_face = X.mean(axis=0).reshape(height, width)\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(mean_face, cmap='gray')\n",
    "plt.title('Mean Face (Lord Voldemort)')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:40.141860Z",
     "start_time": "2018-11-13T16:47:39.601653Z"
    }
   },
   "outputs": [],
   "source": [
    "n_components = 100\n",
    "pca = PCA(n_components=n_components)\n",
    "pca.fit(X)\n",
    "\n",
    "eigenfaces = pca.components_.reshape((n_components, height, width))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:41.150702Z",
     "start_time": "2018-11-13T16:47:40.147771Z"
    }
   },
   "outputs": [],
   "source": [
    "f, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "ax1.imshow(eigenfaces[0], cmap='gray')\n",
    "ax1.set_title('Eigenface of Eigenvalue 0')\n",
    "\n",
    "ax2.imshow(eigenfaces[n_components - 2], cmap='gray')\n",
    "ax2.set_title(f'Eigenface of Eigenvalue {n_components - 2}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Facial Recognition\n",
    "\n",
    "We choose a single test image and try to find other images that are \"similar\" to it, via the distance of the PCA weight vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 28\n",
    "test_img = X[index]\n",
    "test_name = names[y[index]]\n",
    "\n",
    "X_rest = np.delete(X, index, axis=0)\n",
    "y_rest = np.delete(y, index)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f'Image {index} ({y[index]} = {test_name})')\n",
    "plt.imshow(test_img.reshape(height, width), cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for this large an input, sklearn uses a randomized pca by default\n",
    "# make sure to use the normal one we introduced above (slower)\n",
    "pca = PCA(n_components=n_components, svd_solver='full')\n",
    "\n",
    "X_trafo = pca.fit_transform(X_rest)\n",
    "test_trafo = pca.transform(test_img[np.newaxis, :])\n",
    "\n",
    "distance = np.linalg.norm(X_trafo - test_trafo, ord=2, axis=1)\n",
    "\n",
    "best_matches = np.argsort(distance)\n",
    "worst_matches = best_matches[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(4, 4, figsize=(9, 9))\n",
    "\n",
    "for ax, index in zip(axs.flat, best_matches):\n",
    "    prediction = names[y_rest[index]]\n",
    "    ax.imshow(X_rest[index].reshape(height, width), cmap='gray')\n",
    "    ax.set_title(f'{index} ({prediction})', fontsize=10)\n",
    "    ax.set_axis_off()\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Disadvantages of the PCA\n",
    "\n",
    "Interpretability?\n",
    "\n",
    "Especially in scientific contexts, the new features are hard to understand.\n",
    "\n",
    "\n",
    "Assumption: we observe an energy $E$, time $t$ and coordinates $(x, y)$.\n",
    "\n",
    "* What is a meaning of a principal component, which is computed as:\n",
    "  $$\n",
    "    0.78 \\cdot E - 0.23 \\cdot t + 0.8 \\cdot x - 0.2 \\cdot y ?\n",
    "  $$\n",
    "  \n",
    "\n",
    "* Units? Is it just unitless? Or did we find a real physical \"feature\" with units? \n",
    "\n",
    "    $$\n",
    "        \\frac{0.78}{\\mathrm{GeV}} \\cdot E\n",
    "        - \\frac{0.23}{\\mathrm{s}} \\cdot t\n",
    "        + \\frac{0.8}{\\mathrm{m}} \\cdot x\n",
    "        - \\frac{0.2}{\\mathrm{m}} \\cdot y\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Selection \n",
    "\n",
    "* The second approach to dimensional reduction.\n",
    "* Instead of transforming the data, we \"just\" discard features\n",
    "* Goal: keep as much relevant information as possible\n",
    "* This means it is highly problem specific, which features are \"relevant\"\n",
    "\n",
    "General idea:\n",
    "\n",
    "> Good feature subsets contain features highly correlated with the classification, yet uncorrelated to each other.\n",
    "> \n",
    "> -- Mark Hall\n",
    "\n",
    "### Univariate Feature Selection\n",
    "\n",
    "* Only look at one feature (and the target) at a time\n",
    "* Neglects correletions between features\n",
    "\n",
    "#### Correlation with the target quantity\n",
    "\n",
    "* Assume we want to estimate a quantity $y$ from an $N$-dimensional dataset.\n",
    "* However, parts of this dataset have high noise or have no causal or statistical relationship to $y$ \n",
    "\n",
    "\n",
    "* Simplest approach: look for $k$ features with highest pearson correlation, discard the rest\n",
    "\n",
    "\n",
    "Example: a 4d dataset where only two features are relevant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:42.531810Z",
     "start_time": "2018-11-13T16:47:41.154671Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_regression\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "X, y = make_regression(n_samples=300, n_features=4, n_informative=2, n_targets=1, random_state=0, noise=0.1)\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, figsize=(9, 9), sharex=True, sharey=True)\n",
    "\n",
    "for i, (ax, col) in enumerate(zip(axs.flat, X.T)):\n",
    "    ax.scatter(col, y)\n",
    "    ax.set_ylabel('Target Variable')\n",
    "    ax.set_xlabel('X{}'.format(i))\n",
    "    \n",
    "    r, _ = pearsonr(X[:, i], y)\n",
    "    ax.set_title(rf'Correlation $\\rho = {r:.3f}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Feature Selection\n",
    "\n",
    "Multivariate approaches try to take correlations between features into account.\n",
    "\n",
    "Some criteria that might be used for subsets of features\n",
    "\n",
    "* [Mutual information](https://en.wikipedia.org/wiki/Mutual_information)\n",
    "* [Cross entropy](https://en.wikipedia.org/wiki/Cross_entropy)\n",
    "* [Minimal description length](https://en.wikipedia.org/wiki/Minimum_description_length)\n",
    "* Performance of the estimator\n",
    "\n",
    "\n",
    "In general, it's impossible to test all combinations of attributes, as\n",
    "the number of combinations grows exponentially with the number of attributes.\n",
    "\n",
    "Given $n$ attributes, the binomial theorem gives us the number of combinations:\n",
    "\n",
    "$$\n",
    "N = \\sum_{k = 1}^{n} \\begin{pmatrix}\n",
    "    n\\\\\n",
    "    k\\\\\n",
    "    \\end{pmatrix} = 2^n - 1\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#### Correlated Features\n",
    "\n",
    "In the example below, it is obvious that two features are strongly correlated.\n",
    "\n",
    "One of these features is thus superfluous or redundant.\n",
    "\n",
    "To find correlated attributes, we \"only\" need to compute the correlation coefficient for pairs of features, yielding a quadratic runtime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-13T16:47:44.534634Z",
     "start_time": "2018-11-13T16:47:42.537230Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_classification\n",
    "from scipy.stats import pearsonr\n",
    "from itertools import combinations\n",
    "\n",
    "X, y = make_classification(n_samples=300, n_features=3, n_informative=2, n_redundant=1, n_repeated=0, n_classes=2, n_clusters_per_class=2, random_state=0)\n",
    "plots.plot_3d_views(X, y)\n",
    "\n",
    "for i, j in combinations(range(3), 2):\n",
    "    r, p = pearsonr(X[:, i], X[:, j])\n",
    "    print('Correlation between feature {} and {} : {:.2f}'.format(i + 1, j + 1, r))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of \"greedy\" strategies to reduce the number of combinations to be tested.\n",
    "\n",
    "Most of them only look at pairs of features and are thus called \"bivariate\"\n",
    "\n",
    "The two simplest and most commonly employed are *Forward* and *Backward* Selection\n",
    "\n",
    "\n",
    "###### Forward Selection:\n",
    "\n",
    "Start with a single attribute and add further attributes according to some to be defined criterion\n",
    "until a termination condition is reached (e.g. the wanted number of features).\n",
    "\n",
    "###### Backward Selection:\n",
    "\n",
    "Start with all attributes and consecutively *remove* features until a condition is reached.\n",
    "\n",
    "\n",
    "#### Max-Relevance, Min-Redundancy (mRMR)\n",
    "\n",
    "Original Publication by Peng et al. (2005): [ieeexplore.ieee.org/document/1453511/](ieeexplore.ieee.org/document/1453511/)\n",
    "\n",
    "Choose the subset of features $S_k = \\{f_1, f_2, \\ldots, f_k\\}$ that as a whole has the highest relevance towards the target variable $y$ and at the same time has the lowest correlation between the features in $S_k$.\n",
    "\n",
    "The measure of relevance is often a correlation measure or the Mutual Information (more on that later):\n",
    "\n",
    "For the search for $S_k$, it should hold $\\max _{S_{k}}(D - R)$, with:\n",
    "\n",
    "\\begin{align}\n",
    "D(S, y) =& {\\frac {1}{|S|}}\\sum _{f_{i}\\in S}I(f_{i}, y) \\\\\n",
    "R(S)   =& {\\frac {1}{|S|^{2}}}\\sum _{f_{i},f_{j}\\in S}I(f_{i}, f_{j}) \n",
    "\\end{align}\n",
    "\n",
    "The subsets are created via forward selection. The next attribute is choses as \n",
    "\n",
    "$$\n",
    "\\mathrm {mRMR} =\\max _{S}\\left[{\\frac {1}{|S|}}\\sum _{f_{i}\\in S}I(f_{i}, y)-{\\frac {1}{|S|^{2}}}\\sum _{f_{i},f_{j}\\in S}I(f_{i}, f_{j})\\right].\n",
    "$$\n",
    "\n",
    "The mRMR algorithm belongs to a class of algorithms that try to maximize relevance while minimizing redundancy, it has some interesting properties and is especially interesting for applications in biology and genetics, where it is common to have more features than observations.\n",
    "\n",
    "Think few participants to a medical study but tens of thousands of genes.\n",
    "\n",
    "\n",
    "**Feature Selection is most relevant when the number of attributes is larger than the number of observations $d > N$**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_classification\n",
    "from mrmr import mrmr_classif\n",
    "\n",
    "\n",
    "X, y = make_classification(\n",
    "    n_samples =10000,\n",
    "    n_features=50,\n",
    "    n_informative=10,\n",
    "    n_redundant=20,\n",
    "    random_state=0,\n",
    ")\n",
    "\n",
    "X = pd.DataFrame(X)\n",
    "y = pd.Series(y)\n",
    "\n",
    "# select top 10 features using mRMR\n",
    "selected_features, relevance, redundancy = mrmr_classif(X=X, y=y, K=10, return_scores=True, show_progress=False)\n",
    "selected_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, layout=\"constrained\")\n",
    "\n",
    "idx = np.arange(len(relevance))\n",
    "mask = np.isin(idx, selected_features)\n",
    "\n",
    "ax1.bar(idx[mask], relevance[mask], color=\"C1\", label=\"selected\")\n",
    "ax1.bar(idx[~mask], relevance[~mask], color=\"C0\", label=\"not selected\")\n",
    "ax1.legend()\n",
    "\n",
    "ax1.set(\n",
    "    ylabel='Relevance',\n",
    "    xlabel=\"Feature\",\n",
    ")\n",
    "\n",
    "\n",
    "mat = ax2.matshow(redundancy.T)\n",
    "fig.colorbar(mat, label=\"Redundancy\", ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problems\n",
    "\n",
    "Algorithms like mRMR that rely on an iterative approach like forward selection are called \"greedy\" heuristics.\n",
    "\n",
    "\n",
    "It is not guaranteed, that the global optimum is reached using a greedy heuristic.\n",
    "\n",
    "This of course depends on the algorithm and the target function.\n",
    "\n",
    "\n",
    "However, it is almost always difficult to impossible to find *the* optimum in finite time/resources.\n",
    "\n",
    "Interresanting article:\n",
    "\n",
    "https://en.wikipedia.org/wiki/Multivariate_mutual_information\n",
    "\n",
    "\n",
    "\n",
    "All approaches here assume that completely unusable attributes have already been removed in preprocessing, e.g.:\n",
    "\n",
    "* Labels from the simulation not available on observed data\n",
    "* Attributes that have large mismatches between observed and simulated data"
   ]
  }
 ],
 "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": 4
}