{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e268ad08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       ".output_png {\n",
       "    display: table-cell;\n",
       "    text-align: center;\n",
       "    vertical-align: middle;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "\n",
    "HTML(\n",
    "    \"\"\"\n",
    "<style>\n",
    ".output_png {\n",
    "    display: table-cell;\n",
    "    text-align: center;\n",
    "    vertical-align: middle;\n",
    "}\n",
    "</style>\n",
    "\"\"\"\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "58fa497c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 8)\n",
    "sns.set_theme()\n",
    "sns.set_context(\"talk\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51539fa1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A second partitionning-based algorithm: Gaussian Mixture Models\n",
    "\n",
    "## Outline\n",
    "\n",
    "- Recap of Kmeans and limitations\n",
    "- Derivation of GMM in 1D\n",
    "- Implementation in scikit-learn\n",
    "- Parameter selection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48cd7431",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Kmeans: recap and limitations\n",
    "\n",
    "### Recap\n",
    "On a simple 1D example, present how the Kmeans algorithm works\n",
    "- Training\n",
    "- Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf998b0",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Build a 2D toy example containing 2 natural clusters (coming from two multi-variate Gaussian distributions). The clusters should have different scale (one very dense and one more spread out). Each cluster contains 500 points.\n",
    "\n",
    "_Use the code from the Agglomerative clustering slides as a starting point_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d02121b3",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 64732 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([4.5, 4.5])]\n",
    "covariance = [np.array([[0.1, 0], [0, 0.1]]), np.array([[2, 0], [0, 2]])]\n",
    "n_points = 500\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance[i], size=n_points)\n",
    "    for i in range(2)\n",
    "]\n",
    "data = np.r_[data[0], data[1]]\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1])\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac0dab3e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "- Build a 2D toy example containing 2 natural clusters (coming from two multi-variate gaussian distributions). The clusters should have different scale (one very dense and one more spread out). Each cluster contains 500 points.\n",
    "\n",
    "- __Try to cluster these points using Kmeans? What do you observe?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bd1dfe95",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 113760 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "km = KMeans(n_clusters=2)\n",
    "clusters = km.fit_predict(data)\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f829a11",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "- Build a 2D toy example containing 2 natural clusters (coming from two multi-variate gaussian distributions). The clusters should have different scale (one very dense and one more spread out). Each cluster contains 500 points.\n",
    "\n",
    "- Try to cluster these points using Kmeans? What do you observe?\n",
    "\n",
    "- __Can we solve the problem by normalizing the data?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6aa5c368",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 128352 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "ss = StandardScaler()\n",
    "data_scaled = ss.fit_transform(data)\n",
    "\n",
    "km = KMeans(n_clusters=2)\n",
    "clusters = km.fit_predict(data_scaled)\n",
    "\n",
    "sns.scatterplot(x=data_scaled[:, 0], y=data_scaled[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "520d81fb",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian Mixture Models\n",
    "\n",
    "Derivation on a 1D example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "47e24912",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 16740 characters]",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [1, 4.5]\n",
    "scale = [0.1, 1.5]\n",
    "n_points = 50\n",
    "\n",
    "data1D = [\n",
    "    np.random.normal(loc=means[i], scale=scale[i], size=n_points) for i in range(2)\n",
    "]\n",
    "data1D = np.r_[data1D[0], data1D[1]]\n",
    "labels = [0] * n_points + [1] * n_points\n",
    "\n",
    "plt.figure(figsize=(8, 2))\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be07984f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's suppose we manage to find the true cluster centers (using KMeans for example). Then, the decision boundary is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e4ab9fe2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17684 characters]",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 2))\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "\n",
    "plt.scatter([1, 4.5], [0, 0], color=\"k\", marker=\"x\", s=200)\n",
    "plt.axvline(x=2.75, color=\"k\", linestyle=\"--\")\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4bbef06",
   "metadata": {},
   "source": [
    "Does not allow to properly cluster the points. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa7b582b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now if we also have access to the true variances of the distributions of the centers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ce7d1263",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30424 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8fd29f9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question: In which cluster should we place the black point ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4c49e48a",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30820 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([2], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87da1c72",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can correctly classify a point $x$ (cluster with maximum likelihood)\n",
    "    - Likelihood: $p(x|C_k) = \\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} \\exp \\left(\\frac{||x - \\mu_k||^2}{2\\sigma_k^2}\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d11296f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question: In which cluster should we place the black point ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "87e7c92c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30824 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([0], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03d5be72",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question: In which cluster should we place the black point ?\n",
    "\n",
    "Let's zoom in!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fdacd9f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18980 characters]",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-0.7, 0.7, 100)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([0], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48509178",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question: In which cluster should we place the black point ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "65280fcc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30788 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([1.315], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13867b50",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question: In which cluster should we place the red point ?\n",
    "\n",
    "Let's zoom in, and compute the likelihoods. __Can we really conclude that the point is in the blue cluster?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b50b205b",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21732 characters]",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(1.3, 1.33, 100)\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([1.315], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f08f136c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.02794258414879452 0.027912633932334185\n"
     ]
    }
   ],
   "source": [
    "p_x_blue = stats.norm.pdf(1.315, means[0], scale[0])\n",
    "p_x_orange = stats.norm.pdf(1.315, means[1], scale[1])\n",
    "\n",
    "print(p_x_blue, p_x_orange)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62cb755e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can compute the probability that a point $x$ belongs to a cluster $C_k$\n",
    "    - Likelihood: \n",
    "    $$p(x|C_k) = \\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} \\exp \\left(\\frac{||x - \\mu_k||^2}{2\\sigma_k^2}\\right)$$\n",
    "    - Probability: \n",
    "    $$p(C_k|x) = \\frac{p(x|C_k)}{\\sum_j p(x|C_j)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "383c80a2",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18980 characters]",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P blue   =  2.6043079246504417e-19\n",
      "P orange =  1.0\n",
      "Sum      =  1.0\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(-0.7, 0.7, 100)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([0], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "p_x_blue = stats.norm.pdf(0, means[0], scale[0])\n",
    "p_x_orange = stats.norm.pdf(0, means[1], scale[1])\n",
    "\n",
    "p_blue_x = p_x_blue / (p_x_blue + p_x_orange)\n",
    "p_orange_x = p_x_orange / (p_x_blue + p_x_orange)\n",
    "\n",
    "print(\"P blue   = \", p_blue_x)\n",
    "print(\"P orange = \", p_orange_x)\n",
    "print(\"Sum      = \", p_blue_x + p_orange_x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "90a00923",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21732 characters]",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P blue   =  0.5002681058054131\n",
      "P orange =  0.49973189419458686\n",
      "Sum      =  1.0\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(1.3, 1.33, 100)\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], scale[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], scale[1]))\n",
    "\n",
    "plt.scatter([1.315], [0], color=\"k\", marker=\"x\", s=200)\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "p_x_blue = stats.norm.pdf(1.315, means[0], scale[0])\n",
    "p_x_orange = stats.norm.pdf(1.315, means[1], scale[1])\n",
    "\n",
    "p_blue_x = p_x_blue / (p_x_blue + p_x_orange)\n",
    "p_orange_x = p_x_orange / (p_x_blue + p_x_orange)\n",
    "\n",
    "print(\"P blue   = \", p_blue_x)\n",
    "print(\"P orange = \", p_orange_x)\n",
    "print(\"Sum      = \", p_blue_x + p_orange_x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1f46623",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can compute the probability that a point $x$ belongs to a cluster $C_k$\n",
    "    - Likelihood: \n",
    "    $$p(x|C_k) = \\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} \\exp \\left(\\frac{||x - \\mu_k||^2}{2\\sigma_k^2}\\right)$$\n",
    "    - Probability: \n",
    "    $$p(C_k|x) = \\frac{p(x|C_k)}{\\sum_j p(x|C_j)}$$\n",
    "    \n",
    "#### How do we find these papameters ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5f5da58f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13264 characters]",
      "text/plain": [
       "<Figure size 864x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 2))\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0])\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43506b9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can compute the probability that a point $x$ belongs to a cluster $C_k$\n",
    "    - Likelihood: \n",
    "    $$p(x|C_k) = \\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} \\exp \\left(\\frac{||x - \\mu_k||^2}{2\\sigma_k^2}\\right)$$\n",
    "    - Probability: \n",
    "    $$p(C_k|x) = \\frac{p(x|C_k)}{\\sum_j p(x|C_j)}$$\n",
    "    \n",
    "#### If we knew the correct cluster assignements for the points, would we be able to compute the Gaussian parameters?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5a4da38d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18240 characters]",
      "text/plain": [
       "<Figure size 864x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 2))\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02e988d0",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can compute the probability that a point $x$ belongs to a cluster $C_k$\n",
    "    - Likelihood: \n",
    "    $$p(x|C_k) = \\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} \\exp \\left(\\frac{||x - \\mu_k||^2}{2\\sigma_k^2}\\right)$$\n",
    "    - Probability: \n",
    "    $$p(C_k|x) = \\frac{p(x|C_k)}{\\sum_j p(x|C_j)}$$\n",
    "    \n",
    "- If we know the correct cluster assignements of the points we can compute the mean and variance of the Multi-Variate Gaussians\n",
    "    - Mean:\n",
    "    $$\\mu_k = \\frac{\\sum_i p(C_k|x_i).x_i}{\\sum_i p(C_k|x_i)}$$\n",
    "    - Variance:\n",
    "    $$\\sigma_k^2 = \\frac{\\sum_i p(C_k|x_i).(x_i-\\mu_k)^2}{\\sum_i p(C_k|x_i)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a208100",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Gaussian Mixture Model algorithm\n",
    "\n",
    "- If we know the correct parameters of the clusters (mean and variance): we can compute the probability that a point $x$ belongs to a cluster $C_k$\n",
    "    \n",
    "- If we know the correct cluster assignements of the points we can compute the mean and variance of the Multi-Variate Gaussians\n",
    "\n",
    "__Does it ring a bell?__\n",
    "\n",
    "__How can we compute the means and variances ?__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d17f69b1",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Expectation Maximization Algorithm for 1D GMM\n",
    "\n",
    "- Fix K, the desired number of clusters\n",
    "- Initialize the K means and variances at random\n",
    "- Repeat:\n",
    "    - Compute the soft allocations of data points under current parameters of the Gaussians\n",
    "    - Update the parameters using the new soft allocation values\n",
    "- Stop when the total likelihood is not changing anymore\n",
    "\n",
    "__Exercise:__ Try to implement EM to fit a GMM to the 1D toy example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "800f6b20",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33628 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialization\n",
    "\n",
    "means = [0.5, 0]\n",
    "stds = [np.sqrt(0.5), np.sqrt(0.5)]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], color=\"k\")\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5053cfe5",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36484 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  23.403336833443856\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "f554246b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40380 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0a3549db",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40448 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  40.54875472662691\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "2e683d8d",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34972 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e6914262",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35044 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  75.62518759531767\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7c983975",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29084 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "391e5c32",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29160 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  106.90526584779028\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4d1af9dc",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29952 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "de2c408b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29960 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  132.1856163952867\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "2f7810ef",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30256 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "4e622b8c",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30256 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  148.44003898677062\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "b23cdf1d",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29800 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters update\n",
    "\n",
    "means = [np.sum(p_Ck_x[k] * data1D) / np.sum(p_Ck_x[k]) for k in range(2)]\n",
    "stds = [\n",
    "    np.sqrt(np.sum(p_Ck_x[k] * (data1D - means[k]) ** 2) / np.sum(p_Ck_x[k]))\n",
    "    for k in range(2)\n",
    "]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "77c8bf21",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29800 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total likelihood =  149.84601632054645\n"
     ]
    }
   ],
   "source": [
    "# Allocation update\n",
    "\n",
    "p_x_Ck = [stats.norm.pdf(data1D, means[k], stds[k]) for k in range(2)]\n",
    "\n",
    "p_Ck_x = [p_x_Ck[k] / np.sum(p_x_Ck, axis=0) for k in range(2)]\n",
    "\n",
    "labels = [np.argmax([p_Ck_x[0][i], p_Ck_x[1][i]]) for i in range(len(p_Ck_x[1]))]\n",
    "\n",
    "sns.scatterplot(x=data1D, y=[0] * data1D.shape[0], hue=labels)\n",
    "x = np.linspace(-1, 9, 1000)\n",
    "plt.plot(x, stats.norm.pdf(x, means[0], stds[0]))\n",
    "plt.plot(x, stats.norm.pdf(x, means[1], stds[1]))\n",
    "\n",
    "plt.legend(\"\")\n",
    "plt.show()\n",
    "\n",
    "print(\n",
    "    \"Total likelihood = \",\n",
    "    np.sum([np.max([p_x_Ck[0][i], p_x_Ck[1][i]]) for i in range(len(p_x_Ck[1]))]),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0366af23",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Extension to multi-dimensional data\n",
    "\n",
    "- The mean value becomes a mean vector\n",
    "- The variance becomes a covariance matrix\n",
    "\n",
    "The principle remains the same: EM-algorithm\n",
    "\n",
    "- Use current estimate of the parameters to allocate soft clusters to the data points\n",
    "- Use the soft allocations to compute the new estimates for the parameters\n",
    "- Stop when the total likelihood stops changing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b1f7634",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise: Understanding covariance matrices\n",
    "\n",
    "Play with the parameters of a covariance matrix to understand what they represent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "da23dc88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 45488 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean = np.array([5, 5])\n",
    "cov = np.array([[1, -0.5], [-0.5, 1]])\n",
    "\n",
    "n_points = 500\n",
    "\n",
    "data_test_cov = np.random.multivariate_normal(mean=mean, cov=cov, size=n_points)\n",
    "\n",
    "sns.scatterplot(x=data_test_cov[:, 0], y=data_test_cov[:, 1])\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7572f6f5",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## GMM using sklearn\n",
    "\n",
    "__Exercise:__ Solve the initial example problem using GMM.\n",
    "\n",
    "- from sklearn.mixture import GaussianMixture\n",
    "- GaussianMixture(n_components=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "1bb3f99f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 65248 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([4.5, 4.5])]\n",
    "covariance = [np.array([[0.1, 0], [0, 0.1]]), np.array([[2, 0], [0, 2]])]\n",
    "n_points = 500\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance[i], size=n_points)\n",
    "    for i in range(2)\n",
    "]\n",
    "data = np.r_[data[0], data[1]]\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1])\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "c1d47c17",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 128264 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "gm = GaussianMixture(n_components=2)\n",
    "clusters = gm.fit_predict(data_scaled)\n",
    "\n",
    "sns.scatterplot(x=data_scaled[:, 0], y=data_scaled[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcda2826",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## GMM using sklearn\n",
    "\n",
    "__Exercise:__ What is the probability of the point $[-0.5, -0.5]$ to be in the blue cluster? The orange?\n",
    "\n",
    "- gm.predict_proba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "3fa2d2e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 128712 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "new_point = np.array([-0.45, -0.45])\n",
    "\n",
    "sns.scatterplot(x=data_scaled[:, 0], y=data_scaled[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.scatter(new_point[0], new_point[1], color=\"k\", marker=\"x\", s=200)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "e90c9136",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 128712 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proba new point in blue cluster   =  0.5714547896575101\n",
      "Proba new point in orange cluster =  0.42854521034249016\n"
     ]
    }
   ],
   "source": [
    "sns.scatterplot(x=data_scaled[:, 0], y=data_scaled[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.scatter(new_point[0], new_point[1], color=\"k\", marker=\"x\", s=200)\n",
    "plt.show()\n",
    "\n",
    "probas = gm.predict_proba(new_point.reshape(1, -1))\n",
    "\n",
    "print(\"Proba new point in blue cluster   = \", probas[0][0])\n",
    "print(\"Proba new point in orange cluster = \", probas[0][1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a083be9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "- How many parameters do we need to estimate for a GMM with $K$ components on data in dimension $N$?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d52ccd83",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "- How many parameters do we need to estimate for a GMM with $K$ components on data in dimension $N$?\n",
    "\n",
    "$$n_{params} = K \\frac{N^2 + N}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31465c4f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "To reduce this number of parameters, we can place constraints on the Covariance matrices\n",
    "\n",
    "- No constraint: _covariance_type=\"full\"_\n",
    "$$n_{params} = K \\frac{N^2 + N}{2}$$\n",
    "\n",
    "- Force all components to share same covariance matrix: _covariance_type=\"tied\"_\n",
    "$$n_{params} = \\;\\; ?$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ad72db",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "To reduce this number of parameters, we can place constraints on the Covariance matrices\n",
    "\n",
    "- No constraint: _covariance_type=\"full\"_\n",
    "$$n_{params} = K \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "- Force all components to share same covariance matrix: _covariance_type=\"tied\"_\n",
    "$$n_{params} = \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "- Force diagonal covariance matrices: _covariance_type=\"diag\"_\n",
    "$$n_{params} = \\;\\; ?$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a5e49cd",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "To reduce this number of parameters, we can place constraints on the Covariance matrices\n",
    "\n",
    "- No constraint: _covariance_type=\"full\"_\n",
    "$$n_{params} = K \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "- Force all components to share same covariance matrix: _covariance_type=\"tied\"_\n",
    "$$n_{params} = \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "- Force diagonal covariance matrices: _covariance_type=\"diag\"_\n",
    "$$n_{params} = 2KN$$\n",
    "\n",
    "- Force covariance matrix to be diagonal, all dimensions with equal variance: _covariance_type=\"spherical\"_\n",
    "$$n_{params} = \\;\\; ?$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01254e31",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "To reduce this number of parameters, we can place constraints on the Covariance matrices\n",
    "\n",
    "- No constraint: _covariance_type=\"full\"_\n",
    "\n",
    "$$n_{params} = K \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "\n",
    "- Force all components to share same covariance matrix: _covariance_type=\"tied\"_\n",
    "\n",
    "$$n_{params} = \\frac{N^2 + 3N}{2}$$\n",
    "\n",
    "\n",
    "- Force diagonal covariance matrices: _covariance_type=\"diag\"_\n",
    "\n",
    "$$n_{params} = 2KN$$\n",
    "\n",
    "\n",
    "- Force covariance matrix to be diagonal, all dimensions with equal variance: _covariance_type=\"spherical\"_\n",
    "\n",
    "$$n_{params} = \\;\\; K(N+1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d262990f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 1:__ Can we use _diag_?, _tied_? _spherical_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "42f02614",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 65064 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([4.5, 4.5])]\n",
    "covariance = [np.array([[0.1, 0], [0, 0.1]]), np.array([[2, 0], [0, 2]])]\n",
    "n_points = 500\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance[i], size=n_points)\n",
    "    for i in range(2)\n",
    "]\n",
    "data = np.r_[data[0], data[1]]\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1])\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a50c8fb",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example :__ Can we use _diag_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "14578c29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 120384 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=2, covariance_type=\"diag\")\n",
    "clusters = gm.fit_predict(data)\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0dfbfb3",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _tied_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "f98428b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 120632 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=2, covariance_type=\"tied\")\n",
    "clusters = gm.fit_predict(data)\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6ead7ef",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _spherical_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "babfc60f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 123764 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=2, covariance_type=\"spherical\")\n",
    "clusters = gm.fit_predict(data)\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "867ab7a2",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _diag_?, _tied_? _spherical_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "f18d6086",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 39916 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([1, 4]), np.array([1, 7])]\n",
    "covariance = np.array([[5, 0], [0, 0.05]])\n",
    "n_points = 100\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance, size=n_points)\n",
    "    for i in range(3)\n",
    "]\n",
    "all_data = np.r_[data[0], data[1], data[2]]\n",
    "\n",
    "sns.scatterplot(x=all_data[:, 0], y=all_data[:, 1])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3a2fc2",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _diag_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "3619b61a",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 79972 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=3, covariance_type=\"diag\", n_init=50)\n",
    "clusters = gm.fit_predict(all_data)\n",
    "\n",
    "sns.scatterplot(x=all_data[:, 0], y=all_data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa6d0f23",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _tied_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "69dcda32",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 79812 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=3, covariance_type=\"tied\", n_init=50)\n",
    "clusters = gm.fit_predict(all_data)\n",
    "\n",
    "sns.scatterplot(x=all_data[:, 0], y=all_data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "625a8ba7",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Constraints on the covariance matrices\n",
    "\n",
    "__Example 2:__ Can we use _spherical_?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "2201e847",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 78468 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm = GaussianMixture(n_components=3, covariance_type=\"spherical\", n_init=50)\n",
    "clusters = gm.fit_predict(all_data)\n",
    "\n",
    "sns.scatterplot(x=all_data[:, 0], y=all_data[:, 1], hue=clusters, palette=\"tab10\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72aa50d3",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Selecting covariance type / number of components\n",
    "\n",
    "- Increasing the number of parameters (number of components/no constraints on covariance) generally results in a higher total likelihood.\n",
    "\n",
    "- However, models that are too large are not very informative (e.g., each point its cluster)\n",
    "\n",
    "- Hence, we want to find a good __balance__ between __number of parameters__ and __high total likelihood__.\n",
    "\n",
    "### The BIC criterion\n",
    "\n",
    "$$BIC = -2 \\log(L) + \\log(N)$$\n",
    "\n",
    "$L$ total likelihood, $N$ number of parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb9efb1",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Selecting covariance type / number of components\n",
    "\n",
    "__Example 1:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "69e56569",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 64728 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([4.5, 4.5])]\n",
    "covariance = [np.array([[0.1, 0], [0, 0.1]]), np.array([[2, 0], [0, 2]])]\n",
    "n_points = 500\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance[i], size=n_points)\n",
    "    for i in range(2)\n",
    "]\n",
    "data = np.r_[data[0], data[1]]\n",
    "\n",
    "sns.scatterplot(x=data[:, 0], y=data[:, 1])\n",
    "plt.axis([-1, 11, -1, 11])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "aafff2c4",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "spherical__2 5475.201351795416\n"
     ]
    }
   ],
   "source": [
    "n_components_range = range(1, 5)\n",
    "cv_types = [\"spherical\", \"tied\", \"diag\", \"full\"]\n",
    "\n",
    "bic = []\n",
    "names = []\n",
    "for cv in cv_types:\n",
    "    for n_comp in n_components_range:\n",
    "        names.append(cv + \"__\" + str(n_comp))\n",
    "        gmm = GaussianMixture(n_components=n_comp, covariance_type=cv)\n",
    "        gmm.fit(data)\n",
    "        bic.append(gmm.bic(data))\n",
    "\n",
    "index = np.argmin(bic)\n",
    "print(names[index], bic[index])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bb70f17",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Selecting covariance type / number of components\n",
    "\n",
    "__Example 2:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "c3d414c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40108 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "means = [np.array([1, 1]), np.array([1, 4]), np.array([1, 7])]\n",
    "covariance = np.array([[5, 0], [0, 0.05]])\n",
    "n_points = 100\n",
    "\n",
    "data = [\n",
    "    np.random.multivariate_normal(mean=means[i], cov=covariance, size=n_points)\n",
    "    for i in range(3)\n",
    "]\n",
    "all_data = np.r_[data[0], data[1], data[2]]\n",
    "\n",
    "sns.scatterplot(x=all_data[:, 0], y=all_data[:, 1])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "dade282e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tied__3 1941.9421774610244\n"
     ]
    }
   ],
   "source": [
    "n_components_range = range(1, 5)\n",
    "cv_types = [\"spherical\", \"tied\", \"diag\", \"full\"]\n",
    "\n",
    "bic = []\n",
    "names = []\n",
    "for cv in cv_types:\n",
    "    for n_comp in n_components_range:\n",
    "        names.append(cv + \"__\" + str(n_comp))\n",
    "        gmm = GaussianMixture(n_components=n_comp, covariance_type=cv)\n",
    "        gmm.fit(all_data)\n",
    "        bic.append(gmm.bic(all_data))\n",
    "\n",
    "index = np.argmin(bic)\n",
    "print(names[index], bic[index])\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
