![](\"https://img.icons8.com/dusk/64/000000/mind-map.png\")
CS 236756 - Technion - Intro to Machine Learning\n",
"---\n",
"#### Tal Daniel\n",
"\n",
"## Tutorial 10 - Expectation Maximization\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### ![](\"https://img.icons8.com/bubbles/50/000000/checklist.png\")
Agenda\n",
"---\n",
"* [Demonstrations](#-Demonstrations)\n",
" * [K-Means](#--K-Means)\n",
" * [Gaussian Mixture Models (GMM)](#-Gaussian-Mixture-Models-(GMMs))\n",
"* [Expectation-Maximization](#-Expectation-Maximization-(EM))\n",
" * [Formalization](#-Formalization)\n",
"* [Gaussian Mixture Models (GMM)](#-Gaussian-Mixture-Models-(GMMs)-as-EM)\n",
"* [Bernoulli Mixture Model (BMM)](#-Bernoulli-Mixture-Models-(BMMs)-as-EM)\n",
"* [K-Means as \"Hard GMM\"](#-Relation-to-K-Means:-K-Means-as-\"Hard-GMM\")\n",
"* [Recommended Videos](#-Recommended-Videos)\n",
"* [Credits](#-Credits)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"# imports for the tutorial\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib notebook"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## ![](\"https://img.icons8.com/clouds/100/000000/classroom.png\")
Demonstrations\n",
"---\n",
"* Definitions:\n",
" * **Hard Clustering** - clusters do not overlap, an element either belongs to a cluster or does not.\n",
" * **Soft Clustering** - clusters may overlap. Strength of association between clusters and instances (confidence level)."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### ![](\"https://img.icons8.com/dusk/64/000000/k.png\")
K-Means\n",
"---\n",
"* K-means clustering is one of the simplest and most popular unsupervised machine learning algorithms.\n",
" * Typically, unsupervised algorithms make inferences from datasets using only input vectors without referring to known, or labeled, outcomes.\n",
"\n",
"\n",
"* The objective of K-means: group similar data points together and discover underlying patterns. To achieve this objective, K-means looks for a fixed number (k) of clusters in a dataset.\n",
" * This is also called **hard clustering** as the association of a data point to a cluster is definitive, that is, it belongs to one cluster only and cannot be included in other cluster.\n",
" \n",
"* A cluster refers to a collection of data points aggregated together because of certain similarities."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* Illustration: ![](\"./assets/tut_10_kmeans.jpg\")
\n",
" * *E* and *M* stand for the **E-step** and **M-step**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* Animation: ![](\"./assets/tut_10_kmeans_anim.gif\")
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* **Algorithm**: Starts with a first group of randomly selected centroids, which are used as the starting points for every cluster, and then performs iterative (repetitive) calculations to optimize the positions of the centroids. It halts creating and optimizing clusters when either:\n",
" * The centroids have stabilized — there is no change in their values because the clustering has been successful.\n",
" * The defined number of iterations has been achieved.\n",
"\n",
"Explanations and example from Understanding K-Means Clustering in Machine Learning - towardsdatascience.com"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"# k-means example\n",
"from sklearn.cluster import KMeans\n",
"def generate_random_data():\n",
" # generate random data\n",
" X = -2 * np.random.rand(100,2)\n",
" X[50:100, :] = 1 + 2 * np.random.rand(50,2)\n",
" return X\n",
"\n",
"def plot_data(X):\n",
" fig = plt.figure(figsize=(8,8))\n",
" ax = fig.add_subplot(1,1,1)\n",
" ax.scatter(X[ : , 0], X[ :, 1], s = 50, c = 'b')\n",
" ax.grid()\n",
" ax.set_xlabel(\"x\")\n",
" ax.set_ylabel(\"y\")\n",
" ax.set_title(\"some random data\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHwCAYAAABZrD3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5TddX3n8dd7EpIhM9PFFEx0SEUlx13qMbpJrd3GbQLsCqldLYhTy7q0sifbLrSm0G1isWrXVmna0zVnpbYYW91jcUKjrBxBEI9kZXYXJfEQFkQJ9cghBooK1DsTk5DMZ//43m/ncvneO/fe7+d7P9/P9/t8nDNnZu69c+9nvhN43c+v98eccwIAAHEZCd0AAADQPwIcAIAIEeAAAESIAAcAIEIEOAAAESLAAQCIEAEOYCBm9l0zuzDA6/6amc0M+3WBsiHAAVSWmX3AzD4duh1AEQhwoAbMbGnoNgDwiwAHPDCz7Wb2PTNrmNm3zeyC5u3LzewjZnak+fERM1vevG+TmR02s98zs6fM7Akze6uZbTGzR8zsaTP7/ZbXGDGzHWb292b2QzO72cxWdmhP+tzbzexJSX9jZi8ysy+Y2ffN7Jnm12e3/Mw+M/ugmf3v5u/xJTM7s+X+d5rZY83Xvq7t9bz9nhm/y0+a2a1m9iMz+7qkV7bdv8vMHm/ef8DM3ti8/SJJvy9pysxmzexg8/ZfN7OHm7/jd8zsP/X2VwbKhQAHcjKzV0m6WtLPOOcmJL1J0nebd18n6Q2SXitpnaTXS3pvy4+vljQqaVLS+yR9XNK/l7Re0hslvc/MXtF87G9LequkX5D0UknPSLqhS9NWS1op6WWStir57/1vmt//lKQfS/po28/8qqRfl/RiScsk/W7zdzxP0sckvbP52j8p6eyWn/P5e7a7QdIxSS+R9K7mR6v7mq+7UtJNkv7OzEadc3dI+pCkPc65cefcuubjn5L0Zkk/0fxd/5uZ/csOrw2Ul3OODz74yPEh6VwloXChpNPa7vt7SVtavn+TpO82v96kJESXNL+fkOQk/WzL4w9Iemvz64clXdBy30skPSdpaUabNkk6IWm0S7tfK+mZlu/3SXpvy/f/WdIdza/fJ2m65b6x5vNf6Pv3bGvjkubv+M9bbvuQpJkuv9czktY1v/6ApE8v8vf7n5LeHfrfER989PtBDxzIyTn3qKRtSsLiKTObNrOXNu9+qaTHWh7+WPO21A+dc6eaX/+4+fkfWu7/saTx5tcvk3SLmT1rZs8qCfRTklZ1aNr3nXPH0m/MbIWZ/VVzGPxHkr4q6QwzW9LyM0+2fH205bVfKunxlt95TtIPWx7r8/dsdZakpa2v3fY6MrNrm0Pi/9i8Lv9M0pnqwMwuNrN7m0P3z0ra0u3xQFkR4IAHzrmbnHMblYSsk/QnzbuONG9L/VTztkE8Luli59wZLR+jzrnvdWpW2/fXSnqVkp7vT0j6183brYfXfkLSmvQbM1uhZBg95fP3bPV9SSdbX7v53Gk73ihpu6S3S3qRc+4MSf+ohd/pedegOS//WUl/JmlV8/G3q7drAJQKAQ7kZGavMrPzm+FwTElvMu1tfkbSe83srOaCsPdJGnRb019K+mMze1nzdc8ys7f08fMTzbY921z89v4+fnavpDeb2UYzWybpv+r5///w+Xv+k2av/XOSPtAcQThP0hUtD5lQEvDfl7TUzN6nZG479Q+SzjGztK3LJC1vPv6kmV0s6d/mbScQAgEO5Ldc0vWSfqBkCPrFSlY/S9IfSdov6QFJ/0/SN5q3DWKXpFslfcnMGpLulfSzffz8RySd3mznvZLu6PUHnXMPSbpKySKxJ5TMMx9ueYjP37Pd1UqG15+U9EklC/FSd0r6oqRHlAytH9Pzh9v/rvn5h2b2DedcQ8liwJubv8OvKrmmQHTMufZRNgAAUHb0wAEAiBABDgBAhAhwAAAiRIADABAhAhwAgAhFdULRmWee6c4555xczzE3N6exsTE/Dao5rqVfXE9/uJb+cC396vd6Hjhw4AfOubOy7osqwM855xzt378/13Ps27dPmzZt8tOgmuNa+sX19Idr6Q/X0q9+r6eZPdbpPobQAQCIEAEOAECECHAAACJEgAMAECECHACACBHgAABEiAAHACBCBDgAABEKFuBmNmpmXzezg2b2kJn9Yai2AAAQm5CV2I5LOt85N2tmp0maMbMvOufuDdgmAACiECzAnXNO0mzz29OaHy5UewAAiEnQOXAzW2Jm90t6StJdzrmvhWwPAACxsKQjHLgRZmdIukXSbznnHmy7b6ukrZK0atWq9dPT07lea3Z2VuPj47meAwmupV9cT3+4lv5wLf3q93pu3rz5gHNuQ9Z9pQhwSTKz90uac879WafHbNiwwXEaWXlwLf3ievrDtfSHa+nXAKeRdQzwYHPgZnaWpOecc8+a2emSLpT0J6HaAwBFaTSkPXukQ4ektWulqSlpYiJ0qxC7kKvQXyLpU2a2RMlc/M3OuS8EbA8AeDczI23ZIs3PS3Nz0tiYdM010u23Sxs3hm4dYhZyFfoDkl4X6vUBoGiNRhLejcbCbXNzyectW6QjRySmlzEoKrEBQEH27El63lnm55P7gUGFHEIHgEo7dGihx91ubk569NHhtqdMWBeQHwEOAAVZuzaZ884K8bEx6dxzh9+mMmBdgB8MoQNAQaampJEO/5cdGUnur5vWdQHpG5u5uYXbZ2e7/zwWEOAAUJCJiaRXOTGR9DKl5HN6e6gFbI2GtHu3tH178rl1kV3RWBfgD0PoAFCgjRuT1eZ79iRz3ueem/S8Q4V36OFr1gX4Q4ADQMHGx6UrrwzdinJsa2NdgD8MoQNATZRh+Jp1Af4Q4ABQE2UYvi7ruoAYMYQOADVRluHrQdYFsG/8hQhwAKiJqalkwVqWYQ9f97MuIPTCu7JiCB0AaiLG4Wv2jXdGDxwAPCvzcG/ZtrUtppeFd7325Mv8dxkEAQ4AHpVluLdbWJVlW1svfC28K8vfxScCHAC66KfXVoZ91lK1wsrHwruy/F18Yw4cADqYmZEmJ6Vt26SdO5PPk5PJ7VnKsM+6anPGPvaNl+HvUgQCHAAyDBKEZdhnHUtY9VqP3cfCuzL8XYrAEDoAZBhk8VQZ9lmHCKt+F4f1O8Sfd+FdGf4uRSDAASDDIEFYhn3W3cLqtNOSIGw0/K2+7jeMB52PzrPwrgx/lyIwhA4AGdIgzNKp11aGfdbd5oyfe07au7f7PH4/BplmCDHEX4a/SxHogQNAhkF7baH3WaehtGWLdOqUdPTo8+9Pv/ex+nqQaYZQ89Gh/y5FIMABIENrELYOD4+MLN5rC73POg2rq66SbrpJOnnyhY/ptwhKlkHCeM2a7s959tmDt2cxof8uvhHgANBBzL228XFp9ers8Jb89HarujgsFgQ4AHQRc6+t6IAdZJrh8ce7P+fhw/naVCcsYgOAivJRBKWbQRaHDbI4ENnogQNAReWZx+9Vv9MMVd3SFQIBDgAV5mMef7FCLf1MMwzjTUVdEOAAUHF55vFnZ5N94z4PRol5cWCZEOAAgEyNRtLrLuIUr5gXB5YFi9gAAJm6VUUr08EodUUPHACQ6dAh6cUvzr4v9Cle/R6gUkUEOAAg09q1nY/5DLnlq98DVKqKIXQAQKZuW7rybPnq9SzwTj/b7wEqVUUPHEBlMcyaz8REct0mJvxt+crbex7kAJWqIsABVBLDrH6Mj/vb8jXoWeCtQp1mVkYEOIDKmZ/PHxRY4GvLl4/ec6gDVMo4msMcOIDKefrpxYMCw+ej91x0ffdW6Vz95Zcnq/Hf/W5p505p27akuM3MjL/XGgQ9cACVc/w4w6x5pT3O009PQsxHj9NH77lTKVYz6Td/U/rgB/30kNMpmFOnpKNHn39fWUZz6IEDqJzly+t74lWeFd6pmZmkh7ltm/Tkk/56nL56z2kp1l27pB07pKuvTm6/4QY/PeTWufr28G4VejSHHjiAylm5cnjDrGXiY+Gej4Vmnfg8yCSdl280krBu3T6Wt73d5upbhR7NoQcOoHLSQOjnnOoi+OgN9/NaPvZH97LQLI/23vOuXcn3g+4MKKK93ebqW4UezaEHDqCSQp94NextbL72Rw9jm5bPg0yKaG+3ufpWoUdzCHAAlRXqxKsih6E78RVkobZpDaqI9k5NJW+2OinL+eUMoQOAZ0UPQ2dJgyxLP0E2zG1aPhTR3nSqpXUKZsWKZHHk5ZfnH/b3hR44AHg2zGph6XavBx9Mtjxl6SfI2heaSeXpcWbxuTCuVegpmF4Q4ADg2bCGodvn2UdHk9tHR6VjxwYPstbwGh1NepxlC69WRYVtqCmYXhHgAOBZtzlUX8PQWfPsx44ln51LXv+88wYPsjS89u2TNm3K396ilT1si8AcOAB4ljWH6nsbW7d59qVLk/C+8sry9pqRHz1wAChA0XOonMoFAhwAClLksG5s273gH0PoADBkPiq0xbbdC/7RAweAIfJVoa2o7VOIBwEOAEPiu0JbDHuVURwCHACGxFe98lZ13D6FBHPgADAkrByHT/TAASBDWqL00KFkxffUVDLvnAcrx+ETPXAAaDMzI01OStu2STt3Jp8nJ5Pb82DlOHwiwAGgRetCs7SnPDe3cPvs7ODPPYwKbagPhtABoEURC81asXIcvhDgANBiGAvNWDkOHwhwAGjBQrNyKmJRYeyYAweAFiw0K5+iFhXGjh44ALSgRGm59Fu9rk49dQIcANqw0Kw8+llU6KvOfCwIcADIwEKzcuh1UaHvOvMxYA4cAFBa6aLCLK2LCnvpqVcNAQ4AKK1eFxXWsc48AQ4AKK1eq9f12lOvEubAAQCl1suiwqmpZMFalqpu/yPAAQClt9iiwjpu/yPAAQCVULftf8EC3MzWSPofklZLmpd0o3NuV6j2AADiV6ftfyF74CclXeuc+4aZTUg6YGZ3Oee+GbBNAABEIViAO+eekPRE8+uGmT0saVISAQ6g9oZZErRO5UerpBRz4GZ2jqTXSfpa2JYAQHjDLAlat/KjncT4Jsacc2EbYDYu6X9J+mPn3Ocy7t8qaaskrVq1av309HSu15udndV4VVc0DBnX0i+upz8xX8v5eengweyqYiMj0rp1nQubFPFaR4/2fy3n56Wnn5aOH5eWL5dWrvTX5iLMzibBLSVtT9u6dq3/BXD9/tvcvHnzAefchsw7nXPBPiSdJulOSdf08vj169e7vO6+++7cz4EE19Ivrqc/MV/Lj3/cubEx56QXfoyNObd793Bfq99rec89zk1MLDzv2Fjy/T33+Gu3Tz/6UdK+rGswMeFco+H39fq9npL2uw6ZGOw9kZmZpE9Ietg59+eh2gEAZbJYSdDPflbavl3avfv5B3cU8Vr9lh9tPVAkfd65uYXbZ2fztbcIMddQDzmo8fOS3inpfDO7v/mxJWB7ACC4biVBJekrX5F27pS2bZMmJ5M57CJea5DyozGGYcw11IMFuHNuxjlnzrnXOOde2/y4PVR7AKBdo5H0dH31eHvR7fAOKZlXlvz0bHs9KETq7VrEGIYx11Av8bICAAhnZibp4W7b5q/H24uswzuWLev8+Dw9214PCun1WsQYhv28iSmbUmwjA4AyaZ3LTaU9yy1bknKdRS5yby8Jev/90h13ZD82b892sfKj8/O9X4sYDxSJuYY6AQ4AbXqZyy26XGdrSdDdu6V77skenvbRs+1WfvTpp3u/FoOEYRn2X8daQ50AB4A2ZZvLDdmzPX68v2vRTxiWqYhMjDXUCXAAaJPO5RbV4+1XyGHe5cs7X4vTTkvCutF4fq+5lzAMPU1RBSxiA4A2ZVzYlPZsd+2SduxIPh85UnxPtVsVteeek/buHWxxX4xbzsqGHjgAtMnq8a5YkXz9S78kTU+HmasNMcyb9vK3bJFOnZKOHn3+/en3/faayzZNESN64ACQobXHe/nlSXiNjEg33TS8LWVlkV6Lt71NWtqh29dvrznGLWdlQ4ADQAfj49Lb3y7demuymCvtbZa9PGgRxsel1aulkyez7++311zGaYrYEOAA0AVztQt89pp7LSKDzpgDB4AumKtd4Hs7W6z7r8uCAAeALsq2pSykIrazxbj/uiwIcADoIsbyoEWi11weBDgAdBFzrexBtJc2fcUrXvgYes3lQIADwCLq0uvMKm36wQ8mp6ENu7RpLELWcifAAaAHVe91diptmp5GRmnTFwpdy51tZAAAtsv1qfUNT7rAcdj1AQhwAADb5fpUhjc8BDgAgNKmfSrDGx4CHABAadM+leENDwEOAOhY2rSq2+XyKsMbHlahAwAkZW+Xe/nL2UKWpQz1AQhwAMA/ad8ut29fsKaUXuj6AAQ4AAAdLFaoJWR9AAIcAIAMoQu1LIZFbAAAtClDoZbFEOAAALQpQ6GWxRDgAAC0KUOhlsUwBw4AgYU80QrZ0kItWSFelsp09MABIKCZGWlyUtq2Tdq5M/k8OZncjnDKUKhlMQQ4AAQSw0KpuupUmS69vQyV6RhCB4BAelkoVeUzyMsudKGWxRDgABBIDAul6iiWNQkEOAAMQVYoxLBQqm7KXrylFQEOAAXrFAp795Z/oVSdtK5JSKVvrrZsSYbTyzJ8LrGIDQAK1W2h2tveloR4mRdK1UkMxVta0QMHgAItFgqPP17uhVJ1EtuaBAIcAArUSyiEPNEKC2Jbk8AQOgAUKA2FLGUMhTqLoXhLKwIcAAoUWyjUWQzFW1oxhA4ABUr/59++Cn1kpJyhUHdlL97SigAHgILFFAqQnEs+5ucXvi4jAhwAhoCFanGgkAsAoHRiKREaSmyFXAhwAKiBmHqWocR2uAyr0AGg4ji2tDexFXIhwAGg4mIrERpKbHv2CXAAqLjYepahxLZnnwAHgIqLrWcZCoVcAAClMjWVLFjLUsaeZUgx7dknwAGg4qgG159Y9uwT4ABQAzH1LNEbAhxA5czPS7t3U7CkXSw9S/SGAAdQKTMz0sGD0h/8AQVLUG2sQgdQGWlhknSeV6JgCaqLAAdQGRQseb5GI5lK2L49+dxa4xvxYwgdQGVQsGQBtc+rjx44gMqgYEmC2uf1QIADqIzYSmEWhamEeiDAAVRGWrBkZCSOUphFYSqhHpgDB1ApGzdKJ05Iu3bVt2BJOpWQFeJ1mkqoOgIcQOWMjNS7YAm1z+uBIXQAqJjYTtXCYOiBA0AFUfu8+ghwAKgoap9XG0PoAABEiAAHACBCDKEDANCi0UjWDvR6HG2/j/eFAAcAoKnfGvIha84zhA4AgPqvIR+65jwBDgCotfTY1csuS6r4ZcmqIR+65jxD6ACA2mofAu8kq4Z86JrzQXvgZvbXZvaUmT0Ysh0AgPrJGgLvJKuGfOjja0MPoX9S0kWB2wAAqKFuQ+DtsmrIhz6+NmiAO+e+KunpkG0AANRTtyHwVLca8qFrzjMHDgCopW7Hri5bJl1wgXTppd1ryIesOW/OueJfpVsDzM6R9AXn3Ks73L9V0lZJWrVq1frp6elcrzc7O6txqvl7wbX0i+vpD9fSnypfy/l56eDB7GH0kRFp3brOQ+SD6vd6bt68+YBzbkPWfaXvgTvnbpR0oyRt2LDBbdq0Kdfz7du3T3mfAwmupV9cT3+4lv5U/VouW/bCQiwjI8UVYvF5PUsf4AAAFCXmY1eDBriZfUbSJklnmtlhSe93zn0iZJsAAPUS67GrQQPcOfeOkK8PAECsQu8DBwAAAyDAAQCIEAEOAECECHAAACJEgAMAECECHACACBHgAABEiAAHACBCBDgAABEiwAEAiBABDgBAhAhwAAAiRIADABAhAhwAgAgR4AAARIgABwAgQgQ4AAARIsABAIgQAQ4AQIQIcAAAIkSAAwAQIQIcAIAIEeAAAESIAAcAIEIEOAAAESLAAQCIEAEOAECECHAAACJEgAMAECECHACACBHgAABEiAAHACBCBDgAABEiwAEAiBABDgBAhAhwAAAitGiAm9nVZvaiYTQGAAD0ppce+GpJ95nZzWZ2kZlZ0Y0CAADdLRrgzrn3Slor6ROSfk3SITP7kJm9suC2AQCADnqaA3fOOUlPNj9OSnqRpL1mtrPAtgEAgA6WLvYAM/ttSVdI+oGk3ZL+i3PuOTMbkXRI0u8V20QAANBu0QCXdKakS5xzj7Xe6JybN7M3F9MsAADQzaIB7px7X5f7HvbbHAAA0Av2gQMAECECHACACBHgAABEiAAHACBCBDgAABEiwAEAiBABDgBAhAhwAAAiRIADABAhAhwAgAgR4AAARIgABwAgQgQ4AAARIsABAIgQAQ4AQIQIcAAAIkSAAwAQIQIcAIAIEeAAAESIAAcAIEIEOAAAESLAAQCIEAEOAECECHAAACJEgAMAECECHACACBHgAABEiAAHACBCS0O+uJldJGmXpCWSdjvnrg/ZnqppNKQ9e6RDh6S1a6WpKWliInSrAAA+BAtwM1si6QZJ/0bSYUn3mdmtzrlvDuP1qx5uMzPSli3S/Lw0NyeNjUnXXCPdfru0cWPo1gEA8grZA3+9pEedc9+RJDOblvQWSYUHeNXDrdFIfr9GY+G2ubnk85Yt0pEj0vh4mLYBAPwIOQc+Kenxlu8PN28r1Pz8QriloTY3txB6s7NFt6B4e/Ykv2eW+fnkfgBA3Mw5F+aFzS6T9Cbn3H9sfv9OSa93zv1W2+O2StoqSatWrVo/PT2d63WfeWZW3/3ueGbAjYxIa9ZIZ56Z6yWC+973pCef7Hz/6tXSpIe3SrOzsxqnK+8N19MfrqU/XEu/+r2emzdvPuCc25B5p3MuyIekn5N0Z8v375H0nm4/s379epfXpz99t5Ncx48dO3K/RHAf/7hzY2PZv9/YmHO7d/t5nbvvvtvPE8E5x/X0iWvpD9fSr36vp6T9rkMmhhxCv0/SWjN7uZktk/Qrkm4t+kWXL0/mvLOMjUnnnlt0C4o3NZWMJmQZGUnuBwDELViAO+dOSrpa0p2SHpZ0s3PuoaJfd+XK6ofbxESyIG9iYuHNytjYwu2MhgFA/ILuA3fO3S7p9mG+5shIEmLtq9DT232GW8itauvWSR/6kHTbbZKZ9Iu/KF1xBeENAFURNMBD2bgx2Uq1Z4/06KPJsPnUlN9wC7lVLeu1Z2aSUK/CNjkAQE0DXErC+sori3nukPuw2QMOAPVALfQcGg1p925p+/bkcxqaIfdhswccAOqhtj3wvLoNkR86tNDrbTc3lwzbF2WYrz0/n7xxqWo5WgAoM3rgA2gdps6q5rZmTbitamvXDue1Z2akgwelbduknTuTz5OTye0AgOIR4ANYbJjarL+tap2G4gcxjD3g6RuVdPRBql45WgAoOwJ8AIsNUx8+3Ps+7JmZpOfqqyc7jD3gzLMDQHjMgQ8gHabOCvF0mLqXrWpFrRgveptcyDl+AECCAB/A1FSyYC1L6zD1YlvVeunJDrrVrchtcsOaZwcAdMYQ+gB8DVPH2pMdZq11n+sDAKBK6IEPyMcwdS9D8WWUvlG5776F9hdRjjZkNTsAKDsCPIe8w9S9DsWX0caN0okT0q5dxcyzU1EOALojwANKe7LDOFilCCMjxc2zd1sfcPy4dNll0qWX+i0eE/LwGQDoFwEe2DAOVilaEcHXbX3AiRPSHXdI99zjb0id4XoAsSHAS6DIFeNFKyr4uq0PSPkaUme4HkCMWIWOgc3Pdy8pm6ciW7eV7u1OnZKuumrwlerdhutPnEiG61kBD6BsCHAM7Omni6vIlrVVr5OjR6Wbbhq8kl234frjx5Phemq9AygbAhwDO3682H3s6fqAXbukiy+Wli/v/NiTJxdet98RgG6FaVLp815wgfTRj9IbBxAeAR6pMhQ4Wb68+Ips6fqAPXukZct6/7l+RgD6Ga4/cUK69lp64wDCI8Aj5PsAlEGtXDm8imxZQ+qnndb58f2MAPQzXC8lIc7JawBCI8Ajs9hZ5MMMlHS/epEnn7VqHVLfsUN6xzukFSuyH9vvCEDrc190UW+9fU5eAxAS28giU+QBKIMY9j721i13jYZ0yy3ZjxtkBCB97re/PRnROHGi++PLXK8eQPUR4JEp4wEoofaxF1XJrvV5T5xIFutlKXO9egDVR4BHJtYDUIpS1AhA+ryf+pT0O78jPffcCx9T9nr1AKqNAI9MzAegFKWoEYDx8aRAzLp18darB1BdBHhkYj8AJUa99vI71YTnkBQARSDAI1SFA1Bis1gvv1NN+OuvT1bMc0gKAN8I8EjFfABK1XQ7DOWqq57/WA5JAeAL+8CBnLpt7euEPeQA8iLAgZy6be3rhD3kAPIiwIGcejkMpV0dt/wB8IsAB3Lq5zCUVF23/AHwh0VswCIW2waWtbVvdFQ6diypqd5aknV0NDmEhS1/APIiwIEuOm0Pa98G1rq176GHpI99LLk9q576I49Iq1cPp/0AqoshdKCDfk9+S7f2nXeetGRJ9nMuWSLddlux7QZQDwQ40EEvJ79lKeOBMwCqhwAHOhg0iLutSmf1OQBfCHCgg0GDuNuqdFafA/CFAAc6GDSI01XpExMLbwDGxhZuZ/U5AB9YhQ50kOfkNw6cAVA0AhzoIk8Qc+AMgCIR4MAiCGIAZUSAl9Bilb8AACDAS6bXyl8AgHpjFXqJ9Fv5CwBQXwR4iQxa+QvxaDSk3bul7duTz41G6BYBiBVD6CVCCc5qY3oEgE8EeEEGWYiWVv7KCvEqlOCs8+K81umRVPp33rIl2aoGAP0gwAswaE9raip5XJbYS3DWvffZy/TIK1853DYBiBtz4J7lWYhW1RKcLM5jegSAfwS4Z3kXoqWVv3btknbsSD4fORJ3L5XFeZxQBsA/htA989HTqlrlL3qfvU2P7N8/3DYBiBs9cM9i6mkNa0tTTNekKFWdHgEQDj1wz2JZiDbMRWWxXJOicUIZAJ8IcM/yHEE5LL1safLZzhiuybBUbXoEQDgEeAHK3tPqZVGZ75Ap+zUBgNgQ4AUpc08r1KKyMl8TAIgNi9hqiEVlABA/AryGpqaS+ecsdVpUBgAxI8BriC1NABA/5sBrikVlABA3AnxIyngSF4vKACBeBPgQ1P0kLl9k5OgAAAygSURBVACAf8yBF4yTuAAARSDAC8ZJXC80rBrsAFBlDKEXjJO4no/pBADwgx54wSiasoDpBADwhwAvGEVTFjCdAAD+EOAFo2jKAqYTAMAf5sCHgKIpiXQ6ISvE6zadAAB5EeBDQtGU5E3LNddk31e36QQAyCvIELqZXWZmD5nZvJltCNEGDB/TCQDgT6ge+IOSLpH0V4FeH4EwnQAAfgQJcOfcw5JkZiFeHoExnQAA+TEHjqDKeMgLAMTAnHPFPLHZlyWtzrjrOufc55uP2Sfpd51z+7s8z1ZJWyVp1apV66enp3O1a3Z2VuOM13qR91rOzibBLSX7wNP98mvX1nNInX+b/nAt/eFa+tXv9dy8efMB51zmWrHCeuDOuQs9Pc+Nkm6UpA0bNrhNmzbler59+/Yp73MgkedaNhrS5GR2HfSJiWSevG7/z+Dfpj9cS3+4ln75vJ4UckEQVGUDgHyCzIGb2S9L+u+SzpJ0m5nd75x7U4i2IIy6VWVjrh+Ab6FWod8i6ZYQr41yqFNVNk5gA1AEhtARRF0OeeEENgBFIcARRF2qsjHXD6Ao7ANHMHWoyla3uX4Aw0OAI6g6V2Wr2lw/gOFiCB0oSKMh/cVfdL7frDpz/QCGjwAHCrJnj9St0OFVV1VrugDAcBHgQEG6zX9LSQ8cAAbFHDgKV9ciJnXa6w5g+OiBo1AzM0nN823bpJ07k8+Tk8ntVVeXve4AwiDAUZi6FzGpy153AGEwhI7C9FLEpOpbyOqw1x1AGAR4jQx7LpoiJok673UHUBwCvCZCHKjBIi4AKA5z4DUQai6aRVwAUBwCvAZCHajBIi4AKA5D6DWw2Fz0N78p7d5dzNw4i7gAoBgEeA10m4seHZVuuEFaurS4uXEWcQGAfwyh10C3uehjx6Tjx+u5TxsAYkaAV0yjkQyHb9+efG40Os9Fj44mH1mKnBsHAOTHEHqFLLZVrH0u+sEHpY98JPu5uu3TTveTn3568iahLrXNAaBM6IFXRC9bxdK56A9/OPn80z+90CNv12mfdmtt8yefrFdtcwAoEwK8IgbZKtbvPu08+8mzhvYBAINjCL0iBilbms6Ntw+7j4xk79MetLZ5iCpww1DXY1IBlAMBXhGDli3tZ5/2IG8SWnvtrY+VktuPHIlzT3hV35QAiAcBHrHWHuCaNZJZ9uMWK1va6z7tQd4kVPFEsqq+KQEQFwI8Ulk9wPl5acWKJMgXGw4fxNRU0svM0ulNQhVPJKvimxIA8SHAI9StBzg+Ll1/vXT4sP+ype1z5tLibxKqeCJZFd+UAIgPAR6hbj1A55LiLB/+cDGv3TpnPjoq7drV/U3CIL32sqvimxIA8WEbWYRC9wDTOfPJyeRztx5+FU8k45hUAGVADzxCsfUAq3YiWb/b7wCgCAR4hGIcll5spXtse6qr9qYEQHwI8AhVrQcY655qjkkFEBIBHqmq9ADZUw0AgyHAI1aFHiB7qgFgMLUO8NjmXaso9Ip6AIhVbQM81nnXqoltRT0AlEUt94HnORYTfrGnGgAGU8sAH+TsbBSjioVeAGAYajmEzrxruVRlRT0ADFMtA5x5194Mc5FfFVbUA8Aw1TLAY6xkNmws8vOPXQ8AfKplgFetkplvFFfxjzdEAHyrZYBLzLt2Q3EVv3hDBKAItQ1wiXnXVPvQ7oMPssjPJ94QAShCrQMc2UO7p05Jo6PSsWMvfDyL/PrHrgcARSDAa6zb0G4nLPLrH7seABShloVckOg2tDs6Ki1fTnEVH6g2B6AI9MBrrNvQ7rFjySrp885jkV9e7HoAUAQCvMYWG9o97zwWV/nCrgcAvhHgNUZBm+Fi1wMAn5gDrzEOEgGAeNEDrzmGdgEgTgQ4GNoFgAgxhA4AQIQIcAAAIsQQeoQ4lhIAQIB7MMxA5VhKAIBEgOc2zEDlWMriMboBIBYEeA7DDlSOpSwWoxsAYsIithx6CVRfGg1p716OpSxK65ux9BrPzS3cPjsbtn0A0I4Az2FY5zzPzEiTk9Ldd3d+DMdS5lP0m7FGQ9q9W9q+PfncOmoDAINgCD2HYZzznDVMn4Xa5fkU+WaMoXkARaAHnsMwznnu1jOUkjO7qV2eX/pmLEueN2MMzQMoCgGewzAOA+nWM5Sk889PFsvRk8unqDdjw1wnAaBeGELPqejDQBYbpr/0UnrePqRvutqHukdG8r0ZG9Y6CQD1Q4B7UORhIJzZPTxFvBkbxjoJAPVEgJdcUT3DYYulQIrvN2O8AQNQFAI8ArGf2V3nVdhVeQMGoHwI8EjEemY35V/jfwMGoJwIcBSq2yrsU6fqU/411jdgAMqLbWQoVLdV2EePdq8u5wMV0ABUFQGOQq1dK61Y0fn+vXuLK2aSlqDdtk3auTP5PDmZ3A4AsQsS4Gb2p2b2LTN7wMxuMbMzQrQDxZua6l5JbsmSYoqZUAENQNWF6oHfJenVzrnXSHpE0nsCtQMFm5iQLrmk8/1HjxZTzIQKaACqLkiAO+e+5Jw72fz2Xklnh2gHhmPz5s7D6EUVM6ECGoCqK8Mc+LskfTF0I1CcqalkqDxLUcVMijqcBADKwpxzxTyx2Zclrc646zrn3Oebj7lO0gZJl7gODTGzrZK2StKqVavWT09P52rX7OysxtmA60U/13J2NukVS8kQdnpwyNq1xeyHnp+XDh7MHkYfGZHWret8eEko/Nv0h2vpD9fSr36v5+bNmw845zZk3VdYgC/GzK6Q9BuSLnDOHe3lZzZs2OD279+f63X37dunTZs25XoOJPq9lrOzwy1mklUBLq2AVsYKcPzb9Idr6Q/X0q9+r6eZdQzwIIVczOwiSdsl/UKv4Y34DbuYCRXQAFRZqEpsH5W0XNJdZiZJ9zrnfiNQW1BhVEADUFVBAtw5xxIiAAByKNkyHgAA0AsCHACACBHgAABEiAAHACBCBDgAABEiwAEAiBABDgBAhAhwAAAiRIADABAhAhwAgAgR4AAARIgABwAgQgQ4AAARIsABAIiQOedCt6FnZvZ9SY/lfJozJf3AQ3PAtfSN6+kP19IfrqVf/V7Plznnzsq6I6oA98HM9jvnNoRuRxVwLf3ievrDtfSHa+mXz+vJEDoAABEiwAEAiFAdA/zG0A2oEK6lX1xPf7iW/nAt/fJ2PWs3Bw4AQBXUsQcOAED0ahngZvanZvYtM3vAzG4xszNCtylWZnaZmT1kZvNmxkrVAZjZRWb2bTN71Mx2hG5PzMzsr83sKTN7MHRbYmdma8zsbjN7uPnf+LtDtylWZjZqZl83s4PNa/mHPp63lgEu6S5Jr3bOvUbSI5LeE7g9MXtQ0iWSvhq6ITEysyWSbpB0saTzJL3DzM4L26qofVLSRaEbUREnJV3rnPsXkt4g6Sr+bQ7suKTznXPrJL1W0kVm9oa8T1rLAHfOfck5d7L57b2Szg7Znpg55x52zn07dDsi9npJjzrnvuOcOyFpWtJbArcpWs65r0p6OnQ7qsA594Rz7hvNrxuSHpY0GbZVcXKJ2ea3pzU/ci9Aq2WAt3mXpC+GbgRqa1LS4y3fHxb/k0TJmNk5kl4n6WthWxIvM1tiZvdLekrSXc653Ndyaf5mlZOZfVnS6oy7rnPOfb75mOuUDBP97TDbFpteriUGZhm3sTUEpWFm45I+K2mbc+5HodsTK+fcKUmvba65usXMXu2cy7VWo7IB7py7sNv9ZnaFpDdLusCxl66rxa4lcjksaU3L92dLOhKoLcDzmNlpSsL7b51znwvdnipwzj1rZvuUrNXIFeC1HEI3s4skbZf075xzR0O3B7V2n6S1ZvZyM1sm6Vck3Rq4TYDMzCR9QtLDzrk/D92emJnZWeluJzM7XdKFkr6V93lrGeCSPippQtJdZna/mf1l6AbFysx+2cwOS/o5SbeZ2Z2h2xST5mLKqyXdqWSR0M3OuYfCtipeZvYZSf9X0qvM7LCZXRm6TRH7eUnvlHR+8/+T95vZltCNitRLJN1tZg8oedN+l3PuC3mflEpsAABEqK49cAAAokaAAwAQIQIcAIAIEeAAAESIAAcAIEIEOAAAESLAAQCIEAEOoCMz+xkze6B5nvFY8yzjV4duFwAKuQBYhJn9kaRRSadLOuyc+3DgJgEQAQ5gEc0a7fdJOibpXzVPVQIQGEPoABazUtK4kvMDRgO3BUATPXAAXZnZrZKmJb1c0kucc1cHbhIAVfg8cAD5mdl/kHTSOXeTmS2R9H/M7Hzn3FdCtw2oO3rgAABEiDlwAAAiRIADABAhAhwAgAgR4AAARIgABwAgQgQ4AAARIsABAIgQAQ4AQIT+P5nlzKSKJHGLAAAAAElFTkSuQmCC\n",
"text/plain": [
"![](\"https://img.icons8.com/color/96/000000/normal-distribution-histogram.png\")
Gaussian Mixture Models (GMMs)\n",
"---\n",
"* A Gaussian mixture model is a probabilistic model that assumes all the data points are generated from a mixture of a finite number of Gaussian distributions with unknown parameters. \n",
" * That is, we consider a family of models, $p_{\\theta}(x)$, which is a weighted sum of $k$ Gaussian functions, where each Gaussian can have different mean and co-variance matrix.\n",
"* One can think of mixture models as generalizing k-means clustering to incorporate information about the co-variance structure of the data as well as the centers of the latent Gaussians."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* The parameters of the model are the mean, co-variance and weight of each Gaussian, which are unknown and will be learned by the EM algorithm.\n",
" * The parameters: $\\theta = \\{\\alpha_l, \\mu_l, \\Sigma_l \\}_{l=1}^k$ ($\\alpha$ is the weight vector).\n",
" * In other words- \"with probability $\\alpha_l$ draw the sample from $l^{th}$ Gaussian\".\n",
"* **Goal**: *Soft clustering* the data under the assumption that it is generated by a mixture of Gaussians.\n",
" * The optimization method is called Expectation Maximization (EM) and will be used to achieve this goal."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* Illustration: ![](\"./assets/tut_10_gmm.jpg\")
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* Animation: ![](\"./assets/tut_10_gmm_anim.gif\")
\n",
"* Example: the following example is taken from The Python Data Science Handbook"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"# Generate some data\n",
"from sklearn.datasets.samples_generator import make_blobs\n",
"from sklearn.mixture import GaussianMixture\n",
"def generate_blobs():\n",
" X, y_true = make_blobs(n_samples=400, centers=4,\n",
" cluster_std=0.60, random_state=0)\n",
" X = X[:, ::-1] # flip axes for better plotting\n",
" return X\n",
"\n",
"def plot_blobs(X):\n",
" fig = plt.figure(figsize=(8,8))\n",
" ax = fig.add_subplot(1,1,1)\n",
" ax.scatter(X[:, 0], X[:, 1], s=40, cmap='viridis')\n",
" ax.grid()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAHSCAYAAAAwk8gOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9e5Ad5X3n/X26ey6a0UpgZHRBEsIvIUIStpGGywaToDhKRCyQEVcTp1J+y0UuhLfe4vWycUDMK2axy2yKylvsu6n1etfvbrAsCQkxoI2wJzaywYmQZkZgRkLBiZBmRjcQCClzRnPmdPfz/nGmR+fS19Pd53T3+X6qSKw5ffp2nuf5Pb+7kFKCEEIIIY1BafQNEEIIIc0MBTEhhBDSQCiICSGEkAZCQUwIIYQ0EApiQgghpIFQEBNCCCENRGvERefMmSOXLFnSiEtHTi6XQ2dnZ6NvI9XwHYaH7zA8fIfh4Tt0ZmBg4IyU8tN2nzVEEC9ZsgT9/f2NuHTk7NmzB7fddlujbyPV8B2Gh+8wPHyH4eE7dEYIcczpM5qmCSGEkAZCQUwIIYQ0EApiQgghpIFQEBNCCCENhIKYEEIIaSAUxIQQQkgDoSAmhBBCGggFMSGEENJAKIgJIYSQBkJBTAghhDQQCmJCCCGkgVAQE0IIIQ2EgpgQQghpIBTEhBBCSAOhIK4TubyOIx+OIZfXG30rhBBCEkRD+hE3E7phomfXIWzZPwJNEdBNiQduWISN65ZBU7kPIoSQZoeSIGZ6dh3Ctv4R5HUTuUkDed3Etv4R9Ow6FPm1qHUTQkj6oEYcI7m8ji37i0K4lAsFE1v7R/DY2qXobAv/E1DrJoSQ9MJVOkZOn5+Apgjbz1QhcPr8RCTXqafWTQghJFooiGNk7qx26Ka0/cyQEnNntYe+hqV1XyjYa900UxNCSLKhII6RzjYND9ywCDNayl/zjBYF93ctisQsXS+tmxBCSDzQRxwzG9ctAwBs7R+BKgQMKXFf16Lpv4elHlo3IYSQ+KAgjhlNVbBp/Qo8tnYpTp+fwNxZ7ZFowhaW1r2tv9w8PaNFwX0Rad2EEELig6t0nehs0/CZT8+M5dxxa92EEELig4I4A8StdRNCCIkPrtYZIk6tmxBCSDwwapoQQghpIBTETQhLYRJCSHKgabqJYClMQghJHhTETURpKcz81N+29Y8AADatX9G4GyOEkCaGalCTwFKYhBCSTCiImwSWwiSEkGRCQdwksBQmIYQkEwriJqEeDSgIIYQEh6tvE8FSmIQQkjwoiJsIlsIkhJDkwVW4CWEpTEIISQ70ERNCCCENhII447CcJSGEJBuapjMKy1kSQkg6oCDOKCxnSQgh6SC0aiSEaBdC7BNCvC2EOCiE2BTFjZHaYTlLQghJD1HYKPMAfltK+TkAnwewVghxcwTnJTXCcpaEEJIeQpumpZQSwNjUP1um/rOvpUjqAstZEkJIeogkakcIoQoh3gLwAYA+KeWbUZy3nmQpupjlLAkhJD2IokIb0cmEuATATgCPSCmHKj57CMBDADB37txVW7Zsiey6YTlxbgJnc5PT/760sxULZvvTGsfGxjBzZjKLY4R5LgtTFiOwNVWBg7U7NEl+h2mB7zA8fIfh4Tt0ZvXq1QNSyi67zyIVxAAghOgGkJNS/pXTMV1dXbK/vz/S69ZKd+8QtvWXBzbNaFFwX9ciX9HFe/bswW233RbjHYYjl9drKmdZz/SnpL/DNMB3GB6+w/DwHTojhHAUxFFETX96ShOGEGIGgN8BcDjseetBM0QXW+Usg5qjS9OfcpMG8rqJbf0j6Nl1KKY7JYSQ5iQK1WY+gNeEEL8EsB9FH/GuCM4bO4wutqcZNiiEEJIUooia/iWA6yO4l7qTtOjiWs3IUWNtUPI2n1kbFDaNIISQaGjq8FkrutjJR1wvYZi0cpRJ26AQQkiWafqiwxvXLcN9XYvQ3qKgs1VF+5QQ3rhuWd3uIWn+WKY/EUJI/Wj6FVVTFWxavwKPrV3aELOw5Y/N6/b+2MfWLm2I4LM2Ilv7R6AKAUPKum9QCCGkGWh6QWxhRRfXm6T6Yxu9QSGEkGaBK2uDSbo/tlEbFEIIaRaa3kfcaOiPJYSQ5oarfAKgP5YQQpoXCuIEkEZ/bFJyngkhJO1wBU0QafDHJi3nmRBC0g4FMQlEac6zFem9rX8EAHw1ySCEEFIOVZgUkJReyaxBTQgh0UONOMEkzQyc1JxnQghJMxTECSZpZuCk5zwTQkgaoWk6oSTRDMycZ0IIiR4K4oSS1F7JSWiSQQghWYIqTEKJywwcNv83jTnPhBCSZLiCJpSoeyXrhoknew9i++AoNEXACBn4lYacZ0IISQMUxAkmqtKXumFi3XOv4/CpMQDA5NTftzL/lxBCGg4FcYJxMwMHMTE/2Ts0LYRLmQjZ85hlLgkhJDxcPVNAqRnYT25xqYAEgBcGRh3PrdSQ/5u0/GZCCEkzFMQJp1LrdMst3rhuWZWAvH3FPLQoAgXDIfDLDB74lbT8ZkIISTMUxDESxnRrp3XevXIhtg+MYLJCqFq5xQXDxM4Dx8sE5KtDJ1EwnK9z98qFge7Nym/O6/b5zbWauQkhpFnhihkDlhD94b5hqFMRyl+5cXEg062d1rljcBSmtNdsFSHwwsBoleY7oUuoAmhVgckKgXxJRwueXHdtoGdjmUtCCIkWOvRiYNMrB/H83mOYNCQuFExMGhLP7z2GTa8c9PV9p6paed1EhSI6jWFKtDgUAGnTFFxxaUfV3ycKBr69+7Cve7JgmUtCCIkWCuKIyeV1/ODNYVS6ZA0JbH5z2FdpSreqWpoi0KaVf9amCdz5uQVV17x4bYkTn1RX4ppwKJfp1u3JrczlhuuvwOnzE+zCRAghAaBpOmKOnsnBQWGEIYufL79itus53LROTRHYsPIK7BgchWFK6CZgmMDLb5/AjFa1ynfbrgmsXTEffYdOe5qT/UZDV+Y366aJKy/rwI7B4+h96wSjqAkhJABcJWvE0hqdhG4YLK2zvaX651kypwNPrV+Bu1cugqoUP9dNibxu4pPxgs3xndh053Jf5uRSv3Ru0kBeN7GtfwQ9uw6VfcfKbx54Yg1eeeQLuHvlIhz7aNzze4QQQqqhIA6Ibpjo7h3Cyp4+3PHcG3j35Hl09w5BN4qa6JI5nVCFQ7MGRWDJnE5f19m4bhmWXFbt1z16Jocne4ewY3C0Svu14+hH49BUxbNrUi3dnjrbNMyd1Y4dg6OJ6hJFCCFpgoI4IJVaoyllmfbX2abhwZsWQ63w8aqKwIM3Lvad2pPXTbx/Zrzq7xO6xPbB0arzO2GZnq2uSW2awIwWBW2aKCuXWWu3p6R2iSKNwS2+gBBiD33EAfCbQ9t9xzIoAtiyfwSKEDDlRZ+pX9zShAxDYlJ3SQ4uPbYqkllM/VdOrdHQjKImAKutERIGCuIA+M2hjaJV4Mw2DZOGvenZKTq6kjZNYMP1VxQ3B71D05q8RWk1rFq7PUXdJYqkE1ZbI6R2uFUNQFDtz6oRHUQYWT7oW595DQ61Oxy5+tOdaFUBbepXNUxgx+Bx/OWLv/Tl/7XM1+0tCjpbVbRPCVMvTb7W75FsUEt8ASHkIlRXAlAP7a9UsyhFUwBFETBNOG4Gjn08jqsum4ljH49DN03opoRuyqlUJ/vrRaHJR2EByDpZ7lTFamuEhCNbK0IdqMyhVYQRmfbn5IMGAFVR8KP/81b83l+/7iiIC4bEex9UtzvM686qtZsmH5Rav5dlmsF3yjgBQsKRjZWgjlTm0F47fxY2rV8RyaLqVVHLlLBNQ/KDphRLXZZSmr5kwajXaPGbm51m3KqtVY4vQkg1mZghjTD7WdrfsL8sIl/40SwszdtJc3ZCVQTuXrkQLx4YhSoEDCnLNPlm0NzqTTN1qqq0FFWOL0KIM6leBbImPPz6oDetX4GHV1+NW77zU8c+w6VY39+0fgUe/9K1tpsWRr1GTzP5ToPGCWTZZ05IUFI9A7IoPPxqFpfPaseDNy6uEtrtWrF619GPxm2/b+fHbSbNrZ40o+/UK04ga5tnQqIgtatrVoVHEM3CTmjfPyV087rpW+NoJs2tnjDHuposbp4JCUtqV4K0Cg/LJDezTcNYXncUlH4ikC2h/fDqq/HuyfO4dv4sXD6lZWmq4vv5m1Fzqxf0nV4kq5tnQsKS2lGfNuFhmeR+uG8Ypiy2L9QUAVURNZvmojLztWkKrprTgcOnylOf2jXBqNeQMMf6ImndPBMSN6l1yqQtZcIyyU0aRSEMXGxfWGs6S1SpMT27DuHoR9UNJpbM6WxKzS0OaqmyljXStnkmpF6kVhAD6Smt6FQC0KKWUoBRlRW0zjNhc29Hp3oMExIFads8E1IvUj3y02L2czPJWQQ1zUVl5qO5kNQT+swJqSZ5UqsGkl5a0c0kZ+HXNFca7OXXzOeWs0lzIaknadk8E1JPOAPqgFMai4WfdBa7wKyr5nTg6JkcJkpqSZeey08wF1NsSCMo3TyzuAdpdjjq68TF0pTDMMyLUdOaInyZ5uzyL49+NO5avMNvzubGdctQMIpdmlRFwKS5MFYoeIqwuAchRZp3FagzlSY5rzziUpzyLycKJo5+NI6f/7vVVefym7NpLYY7BkehCqBgmLh31UIuhjGQdcETdIPB4h6EFKEgDkhYbabUJHe5z+94BVSN5fUqH7nfIKwne4ewY3C0rFXizgPH0TK1cSDRkVXBU8sGg8U9CLkIR7pPnBab1bPjv3YtAVVe37mssxV/+eI72LxvpOpzLobRk2XBU8sGg9H6hFwk/fawOuFUPOPEuYnYr11L/qX1nXab72y4fiG6Xz6IHYOjjte0FkMSDW69ptP8rmvNZ2e0PiEXoSD2gdticzY3GagQR62UFi/paFXRpim46/qFjgFVumHClBKTJRqYqggs/lQHXugfxktvnXAt1sHFMFqyKnhq3WCwuAchF6Eg9oHbYmN9HjeaqmDjumXYcP0VKBgmVAHsGBxFz65DODc+iSMfjpVtCHp2HcL2gVGUrv3SlPjnD3NwKPA1TRtrTEdOVgVPmA1GWirjERI3oWe/EGIRgP8JYB4AE8B3pZT/T9jzJgmvghz10mZ6dh3CzgPHUTAkCkbxfp7fewzP7x1Ge4sy7bd+dM01tv5IEwA8CosAwN0rnTVtUjtZrCoVJg+dxT0IKRLFqNcB/F9SykEhxL8BMCCE6JNSBu9ikFDcFptLO9XIFg+3iGynYJ+iPJbITRoAikEy5y4UPEtq2tGmKbh75UJ8a8N1nvcYJP2KFMmq4Am7wUh6ZTxC4ib0KiClPAng5NT//lchxLsArgCQGUEMOC82C2afCX1uP+kfp89PQHW2jk9zoWBi99ApSOmt+ZbSpinT13S7x8o2jooQuGfVQjy1fnmsubBZKoKRNcGT1Q0GIfUi0tkihFgC4HoAb0Z53iTgtNjs2bMn9Ln9pH/MndWOCZ+dkDRFYM2yefjRwVNlGrwqAAgBo8Q83a4JrF0xH0/fdZ3r4lnaxtGiaK6X2LxvGIPDH2PXI7dGLoyzXgQjS2Rtg0FIvRBBNSfHEwkxE8DPADwtpXzR5vOHADwEAHPnzl21ZcuWSK7baMbGxjBzZu2LjymBd0+eh2nzOyhC4Nr5s6CIotA7fPI8/PxaihBYOu/f4PS/5nE2Nzn990s7WgAInB0v+VtnKxbMdvdx66bE4VP/6qllf6qzFVdcMsPHHZbj9g5PnJvA2dxk2ftRhPB1381E2HFI+A6jgO/QmdWrVw9IKbvsPotEEAshWgDsAvAjKeWzXsd3dXXJ/v7+0NdNAnv27MFtt91W8/ePfDiGO557Y9rHW0pnq4qXHr4Fz+89hs37hqcDtNywgmQsTdrOpOvXzDttjt4/UpYG5USrJnBg4+8GNks6vcNcXsfKnj7bNKv2FgUDT6yhCXSKsOOQ8B1GAd+hM0IIR0Ec2rYnhBAA/huAd/0I4WYjl9erUotK8Ur/+P4v3se2/hFPIdzhkP5hmQtLBZbd3+yYNkf7NYlHXJgiq0UwSDlec4SQrBOFOnELgD8E8I4Q4q2pv/2llPLvIjh3avHr2+xs03D3yiuq6j3PaCkW7Cj+3VkQFo+7Al+/9TORBsk4RWm7YSLaVK6sFsEgRej/J6RI6NEupXxDSimklJ+VUn5+6r9MCuEgO3enkpg9uy4Gk+uGie7eIewYPA5jSt5pSrGgxn1di/C1W5a4FhJp04oa8FPrV1RpuGG1DK8iJpW0x1AEJKtFMOpJkrVNP3OklCQ/CyFh4ErmA7edux1+C/yXLkQWqlLM5d20fgVyed1RI2zVBF5/bDUur9AKo9Iy3LTRNk3B76+Yi/81dAqaEDAB3B9TYYosFsGoB0nXNoM0wUj6sxASFgpiH9ilF23ZP4yCYeJ3P1V9vJ/OMnNntdsuRHndxIsHRvH4l671rFpUKYSd77XoY3Yr1FGJ17U3rV+B/1CH3F7mqNZG0lsuBum+lPRnISQs3E564NTwIa9LbN43guOfXIBulH/mx7fpNxApSD1e53s1sXnfMP7yxV9CN0zfJj6va/sN+oqCel4r7dTaEame+PX/p+FZCAkLVzUP3HbuAPDJeAE9uw6V7cz91N/1uxAF0Qi9qm9tHxzF4PBZvH9m3JeJj9poOklDr1+/NarT8CyEhIUasQdeDR9MKW135n60ySCBSKUaoZNG61V9a1KXOHxqzHdwjN21SfJJS7S5H2tPWp6FkDBwZfXAEphb9g+XpReVYrcz96NNBg1E8hO0IiAAX/W3ilwomNiyvxgcA4CabwYI0xGpnviZI2l5FkLCwFHsg43rlqFgmNi8b8T2c7eduVv93aCmX6+gldPnJ9DeothW6XIjr5u447nXMXr2AlpUhVGpGSBN0eZeNarT9CyE1AIFsQ80VcG3NnwWgKgqsKGI8Pmzforl+0n38DKju3HkzDgAYNK42E4RYFRqWsmSfz9Lz0KIHVR3AvDU+uV44IZyn9alna112Zn7ibJ28jsXjwGWzpuJNs1fkQ5GpWaDLPn3s/QshJTCER0Au535/n98oy7mW79BK5VmPN2UuH3FPGy6czk62zQ8uu0tvPz2SV/XZFQqiZIs9ZQmJEo4G2ogaN/VKBYgp6AVq7QkUOzkNHdWu6sZb2Zbi+9rMio1GwQdf1ELTFbGIsQdCuIYiWIBshbFmW0avnrzlcjrJl7oH4HVjKlgAHuPfITrn/pxVaBV5WYhl9exY3DU13UZlZp+go6/uAQmK2MR4g5X2RgJswDphoknew9i+8AoTCmhmxKqAEwJKIoApvpIG1Lin06PAfAOtHIrjqAIQFWAVlVlVGqdqNQ8o9ZEg46/OARmkJrShDQrnAEx8cH5CWzeN1zVR9jPAqQbJtY99zoOnxor+7t1KsMjMtrpGm5+5lZNwc//3WqM5fXIfXilAgZgrrKd5nnVnA4c+TAXWfpYUAEY5PggGwZWxiLEm+ZcCWPAWpwu62zFs33v4YdTjRbsKF2A7Ba1J3uHqoRwUOwWOT9NJC738Yx+hWipwFEFMKGbEBBob1E8O1hlGTvN0/q9o0ofCyoA/Ry/+FMdgU3XrIxFiDcUxBHQ3Ts0vThdKBgQQrhqrYaUuKyztex71qL26Jpr8MKAPz+uG06LXC3FEWr1Hdq1eQTkdMERS9isnl3DA6YUJ82zkrCm26AC0M/xtZiuWRmLEG84C0Jy4twEtk6Vv5zWJqSzELYWoGf73rNd1M5dKKBFEY7atB/cFrlaiiPUsgD7ETiWsPmt357h78EygFcTkVLCmG6DCkCv4wHU7OtlZSxC3KEgDsG58Ul8PJbHhO7vNbZqAvd1FbXeG7/1E9tFbffQSQD+im50tKowTYklczrw/pkcNEXxvcj5TcGqNdjGr8BRhahqI5llglQ/C2u6DSoA3Y4f/ni8Zl8vK2MR4g5nQwi6Xz6I/83nsW2agtcfW43LZ7XjyIdjjouapihYs2wuXh06iQmnJhOKwL2rFuKh3/xMbBG3FrUG2/gVOIaUTZVL6qR5VhKF6TaoAHQ7Pgpfb9D8e0KaheZZASPm3PgkXn77hK9jZ7QoeOCGYjAU4L2obbpzOe6/YTHaWxTM0EQxtUgIzGhR0aYp+OpNi/EfvryirNxfXOX/al2A3cptWlgtHx0qd2aWqvZ/mjJdftSpHWAYrP7Xp89PlJUsdWqnaTeWgrbtJIT4h7OnRrpfPgg3hU8VQHuLfU6ulz9udkdrlWYCNCbtJ0ywTampU0GxyxOEQLtWbkJ/4/Wfx/0YicJJ84zDqmEXaFf0+Ups6x+1Db5zug/6egmJBwriGsjldeweOuX4uSKAN/79alwomI6Lqp9FrdKU1yizXq0LsJ3AAZhHbFH5+8ZhurULtNu8bxiQEoZEWfCdKSUUIRyj4+nrJSQeOItqwCsQ6c7PLcD8Szpcz5GmRS3svSZlQ5F2ghZGcQq0s0utu1AwsfnNYbRoimd0PH29hERLMlf+hOPmN1UVgU13Lvd9riDRy0fP5AAAS+Z0NkRocwFuDH4Ko9jlYgdJlQKKlduMAktRElJvOLNqoNRvWkr7VODK7I7WyK6lGyY2vXIIm98chjGVn6wI4A9uWozuO5Y3VcRxs+KnMMqvd1VHvAVJlXKDpSgJiReu4jXyzduX4srLys3PSy7rwDdvXxrpdXp2HcLmfReFMFBs/LD5zWH07DoU6bVI8rDMy26pThcKJs7mJm2jn+0inVVFQK2Q2+2agCrsw9dZipKQeKEgrpFv7z6MYx+Nl/3t2Efj+Pbuw5Fdw1qE7Xx6hgS27B+uWnwbgVMaDAmPZV72e2wlValSLQoevHExHrypmB7X0aKgVRPYsHIhHrxpcU3pSfz9CQkHTdM1UK/WbqfPT0Bx0FIAQGmwyZAN3+MniHnZTmt1CrTTDROGCbwwMAJVEdg+MIr7b1iEe1YtKv7NR3Q8f39CooGCuAbq1dpt7qx2mC51q80GmwzZ8D1+/FTimtGi4NJO1XXzVxlot+mVg9i6fxiGxHRd881vDuPBmxZj4Ik1vqLj+fsTEg3cttZAHK3d7Mx71iKs2pgmVQE8cMPihkWyOvkuLasAzZTRUWpe7mhRoIqin7e0CteC2f7HXC6v4wdvDqOyr4gxFXsAXEwxczI58/cnJDqoEdeAU9R0LfWBvcx7G9ctKwZn7Rue9hUrAri3DhWN3Co9seG7O1FWyfJTGGXPnj2+z3f0TM6xKpwhgX/5YAw7BkddTc78/QmJDgriGrGEoJI7gs5W+1KWfvAy72mqgp4vr8A3fvcaPL7zHbx68BRaVQU7DxxHm6bE4o/z4/tjw3d7KnN+C6bEvasW4qmp3zIM9SqM8t2f/wt+cvgDV5Oz2++vmxIzmXNMiG9omq4RS0u5dv4svPLIFzDwxJppwekXL/PeB+cnpk2Dz/a9h58c/gC6CYwXTOR1E9v6R2JJYSrdHOQmDdtrsQmAPT27Dk31pzYxXjBRMCQ27xvBuudeT0y7xyVzOh1TlVRF4MeHTnuanN2aehimiVufeQ3dvUOJeWZCkgwFcUgUgZq7HrmlpuiGiS888xrueO4NrOzpw9/uPVYXf1wQ359dakwzNwGw3p1d+8rDp8bwZO/BhtyTXezBgzctroo9UBWBddfNQ4vDZtIyOVuU/v6l41g3EetGkZCs0ZxqS0JwN+8BME1Mepwjan9cEN9fmupl14PT5yeqCmWUsmNwFI9/6dq6vCMv90L3HcugCGDL/hEoQsCUxc8fXXMNXj34E9tzVrocrN//4dVX45bv/LTqeJbHJMQfnB0NxG+TeDei9sfW4vtlDeoic2e1o+CS86sq9Qtisos92FoRe+C0iQra9nIsr6NVVVAwjKrPGLhFiDc0TTeYSvNuqyZ8V1KKwx8bt+83y1WYOts03LtqoePn9cr7dnIvTBRMPP/mMM6NX7SzWJuoyr7DQVwODNwjJBzUiCOklpSVSs1kZpuGW595zXZhUxUBTQE0RYm1KXscDeDdTKVZ4qn1KzA4fBaHT42V/b29jkFsbu4Fw5Tofvkg/vqB6x2/H9Tl4GTZqSWdj5BmhDMkAqIo9Vdq3nVb1Orhj43D9+uWpmXXwi+taKqCXY/ciideGsIL/SPTRTMmdROmlNANM/byj15lMXcPncTT+etc+xhbv7tfk3IcmzdCmgUK4giIutSf26KmqUrd/G1R+X69anP/1m/PCH2NJKGpCto0BS2qgDEVQW1KYPvAKBQhHMeEX4uK13GdbRpuXzEPL711wv7+FMXWbxtmQ8nAPUJqhzMlJKZE5A0gShe1o2dyAIq5n2ktpO8ViZ21XNOLG49yrdRpTJQKwNII5koBGMS8v+nO5Xj57RO2FbSc/LZRbCgZuEdIcCiIQ6IbZiyl/nTDxDOvHs5EZxuvYJ60PY8XQcs/bnrlUFkJUwB4/s1hmBLo+fJFARjEvD+7oxV/ePOV2Lp/uCyv2clvG7SjWJQlPAlpdrK1AjYATVViiRj1U90qLXhFYvsMEk8NQaKIc3kdm98cruo5bZgSm/dd7DftVWjF7nIb1y3D/Tcs9hX97FZcprSQh26Y6O4dwsqevuliM6ygRUg4KIhDoghEnu6Txc42zVSFK0gK2NEzORgOrS4NU+LomRxyeR39Rz92LBbiZN63XBwDT6zxLMPqd/MQ1QYxy2lshASFNqUIiDpiNIudbZotmCeqMfE3e/4Zfe9+AFURGHco+uJl3vfjt/WTghTUfG1HFBkGhGSN7K6EdSRqIRO2QEKS/XfNEszjd0wsmdMJRcCxLeHfT3VBcsISlIo4E/qevTYPUWwQo84wICQLJGuVTjlRCZlaCySE1TaSLMDTiteY6GzT8Ac3LS76iUuEsSoAKYvVsOzoaFFgAtOC8o3Xfx76Xr02D343iE7jKAqNmpAswlGfUGoxbdaqbdBc2Fi671gORQhs2T88nb50+4r56Dt0GrnJ6vrNHa0q/uYPVqJryadiEVxOmwevDWKbpqC7d8hxHGXR5UJIFFAQJ5Sg5u4w2gbNhY3F7rcGgN1Dp2yPN6V0FMJRWjXszrVx3TIUDIkdg6NQleKmwdogeo0j1qQmxB4K4hO65ZIAACAASURBVITj19xdq7ZBc2FyqPytg7on3LTRIFgWkh/uG4aqCBimxFduXIxv3r4U3959uCiEBVAwTNy7aiE2rluGvG76Gkd2z9SuCaxdPi/QPRKSJbjCZoSguauWpkNzYXIJ4p44cW4iMqvGplcOXvRZTzmun997DHuPnMHwxxfKhO3OA8fRoir4o99Y4msclT6TAiCvmygYQN+h09g9dIouEdKUUBBnBD8BXna+4LtXXkFzYULx657I5XWczU3iQkEt+3stVo1cXscPpqp6lWJI4J9O56qOt67x8OqrfY2j0md6fOc7ePXgKUwUzGlfOF0ipBnhtjNDeBXNsCvGsPPAcVw1pyO2/sMkPHY9g0uxql7ZoQqBo2dyvotnHD2Tc0ylcrvGWF4PXNhm99CpqqjwNBetIaRWIlllhRD/HcA6AB9IKbmVbRBuGpSbL/jIh2O4Z9VivHhglC3sUoib1WKiYGDD3/yDL7+xbpj4mz3/Evj6lsYbxJROlwghF4lK3fn/APwnAP8zovORENgFeJ345AIM0z4nddIAxid1vPnNL+Kj3CTziFNGZ5uGSztbMaNFlrklVAFACN9+455dh/D379pHagOAANCmCdcmEn4j/WuNaeC4JFkkEtO0lPLnAD6O4lwkHr7/i6NwKdCEVw+ewrN977maQEljcavPvGB2e5lbok0TAERVMwkn069lMZnQ7YWjKgQevGmxryYSXqZ06xgvUzYbTJBmgStuwolCG8jldewYHHU9ZoLpSg3D6zf2W3ClVBu9MGngvv/yj7YFQexMv26mYqCYu6yIYhxCVKVcvUzZzG8nzYKQDp1fAp9IiCUAdjn5iIUQDwF4CADmzp27asuWLZFct9GMjY1h5sx4fFknzk3gbG5y+t+XdrZiwezgUcyTuolffTAG0+O3VoTAr10+E61afWP44nyHScfPb2wdU/r7KUKUHVv5Dk0JvHvyvO1vrgiBa+fPKms/6Xa80zWjwpRTfb1VZfqeTAkcOnkeduuT3f1HQTOPw6jgO3Rm9erVA1LKLrvP6qb6SCm/C+C7ANDV1SVvu+22el06Vvbs2YM4nqW7d2gqFeliSsqMFon7uuYE1gZyeR2P9PS5Ng8AgPYWBQN3/2bdNeK43mHS8fMb5/I6Vvb0Ia+rVd9v00wMbvwCOts023f4s+nzV6ezPby6egzZHV9Je4vEwBNfiHWM6IaJb7zwNl76ZXW6FAB0tqp4ZU1X5MFczToOo4TvsDaYvpRAouxHbJk97155RZU/rhSmK9UXv7+xZTK2I69LPL7zHUefadAe0NbxbS4WEcusHSc9uw7h1YPOQWPMbydZI6r0pR8CuA3AHCHEKIBuKeV/i+LczUgUqR12fsWr5nTg/TM5aIqCCd0EpESbVt7Fh9QHv7+xW3QxUAyymz2jBatnV38WtF65dfz4pI4XBo7bHhO3EHRKs7No1wQ3jCRzRDKapZRfieI8pEgUxfHtAl2OfTSOu1cuwtdvvWr6HEwLaQx+f2Mrunjr/mHbiGYryO63fnuG47WCtOfM5XW8/PZJx883XL8w1rHiFTS2dsV8bhhJ5qBpOoH4Se1ww83s+eKB0WnB6yfNhMRDkN9447plWLtivuO5VCEiS+lxM4VrCvC1W5ZEch0n3DYobZqCp++6jnWoSebgiE4oQf17pbgtpvXw8UWBW85sVvD7G2tqUQA5+W4NKSMTTm6CUFMVLLjEWfOOArcNygM30CRNsglHdUIJ6t8rJc19X/3mzGaBIL+xV1MPRZyJ5J78NA+pJOrKV0FKZRKSBSiIE04Q/17pd4IupkmhGYs4+P2N3QTUG6//PLL78SsI49o0hdmEEpJGOLozShq1CrfGFKz6Fb2ActJk/V4n7k1TLZtQQtJI865qGSeNWkUzdOSJwowbVkDZabK3r5iHTXcux+yOVl/X4aaJkOjgTMk4adIq0uzb9qJRvm87wW+nyb701gm8/PYJ/OHNV/q6p2bYNBFSLyiISWJIs2/biyjMuEG06XPjk+h++SB2D50qE/yPrrnGsWCGKYGt+4d93VMtmya2MyTEHs4GkijS6Nv2IqwZN4g2bR37t3uPwZKTpYL/3IWCa8GMCV36uqcgm6YorAEU4iTLcETXES4m3qTRt+1FWDNuEG26Z9chbN0/DDtl9ULBxO6hU7YdjYLeE+B/0xTGGtBM6WykeUn3CpcSuJgEJ02+bS/C+L6DaNMXj3UWtJoiMHdWO46cGXc8xq8/3s+mKaw1oBnT2UjzQSlQB0oXk9ykgbxuYlv/CHp2HWr0rZE64FXOEoBjFbEgVdKOnslBEe5Neid0E8fPOgvhWrpwuZVKDVPlza1U65b9wbqQEZJkKIhjJsqWhiS92JWzvGfVQpgSWNnThzueewMre/rQ3TtUVjfajzatGyZOnJvAXf/5F7hQMBzvoV1TIKVE3uGQNs1/GVW/hLEGuLeANF1bQBKSJiiIYyYLdZ9JeCwz7sATa/DKI1/AwBNroAiB7QPulhI/zSF6dh3C2dwkJg1nk7SqCKxdMQ/tDvWqZ7QoePFPfwOb1q+I1F0SpoGJZwvIoZO0KpFMQEEcM1nOjSXBKfV9+7WUuDWHsCwupkMAVqsq8OXPL8DgE7+Dp++6Dk6yWgJYMqczkmespNYGJpYQb2+xX6asCG9alUjaYbBWzGQ5N5bUTpBIaregqOGPxx0tLjNaFGz/k9/A8itmT/+tEWMxTCT8xnXLcO5CAS+9dcL2cxYPIVmAGnEdCNPSkDSOOFsx1mIpsQuKcjuPKWWVlhvVWKzl3dTS/9pqAdmq2m82aFUiWYDqWB3IYm5slqlHullUlhLrPEruSNVnBV3imVcPl923NRYfXn013j15HtfOn4XLAwiyRqTidbZp+MqNi2lVIpmFI7iOZCk3NsvUK3c1qipiG9ctw7Zdo1BgoNTjbAJ4fu8xmFLiL26/FqfPT+CyzlY82/dezYK0UXm9Way4RogFBTEhJdSzq1BUlhJNVTBvVjsgcsWoqxIMCfzt3mFs6x+FpghMFEwAEoZEYEHayI5LtCqRLEMfMSElNCLdzMt36scfO6kbtmUtLawUKUPKqshpvzntcb8bP89Zi5+ZkKTD0UxICUlKN6unP9Yr+jiX13Fh0ojl3bAELGl2KIgJKSFsEFWUjT2c/LEFQ+Lrt15Vdo1WTZ32ndaCkyCtFJIFw4SqCBglAjls0BTrSZNmh4KYkApqCQyKWqtz88du3jeMlw6MwpCYvoYigAdvWozN+4bLhKQf3ASpnZBUIaEqAu2aEjpoqpF+Z0KSAkc4AcAWjaXUEhgUtVbnVvADAMantPUt+4dRMEz87qeA7juKAnnL/hEoQsCUElfN6cDRj8angrSKqIoApER7i+oqSJ2EpCGBNhXY9sf/FkvmdIYaL0FbRHKckizCkdzk0D/njFO6WaUwiEOr86qzbJHXJTbvG8Hym4u/VeUGok1Tij2KS7T7+7sW4dE11+Cj3KSrQCt2c7K/rqYomNGqhhaGfn3yHKcky1AQNzn0z/nHSRh89eYrA2l1fnDyVTvxyXgBPbsOYdP6FVUbCCftXlMVW+2y9DkrNxcWfoOzvDRYvz55jlOSZSiImxj654LhFjwVRzRxqa9aEQLjk84tDk0pXX+zUuHspV2WPqcdfoKzgmiwXj55jlOSdTh6m5ig/rlmxk0YvHhgFHevvAI7DxyPtARjpa/6e68fwY7BUeR1e6Hv9zdz0y4fW7vUVRNu04Sv4KwgGqyXT57jlGQdOleamCTlzCYdr2IWX7vlKttmCo+uuSZ04whLm31q/QrcvXKh43F+fjNrQ+HUfvHomZxLNycVL/7pLZ49i72u4fQunIp1cJySrEONuIlhi0b/eAmDBZfMKNPqrJrON37rJ5EFF2mqgm9t+CwAMaUZX/zNFCFwv4/fzEu7BOD4nBIXuzm5+X6j1mA5TknW4QhuclhM3x9+hYGl1XX3DlWZZrf2j+DchQKevuu6UMLjqfXL0aKKst/s0k4Vf+zjN/PaUCyZ0+n6nG2agu7eIVffbxwaLMcpyTIUxE0Oi+n7x68wcPInTxRMvPTWCfzdOyfxlRsX16wd2/1m+//xDV/n8rOhcHtOP75f6xpb9w9jQo+mAhfHKckyHMkEAFs0+sGvMPAqxjFpyEhSb2r9zSxBu2X/8HThj9INhdNz+o1e1g0TpgQKJUHeqgDuWbUwtAbLcUqyCIO1CAmIVwcgP8U4/HY8ihdR8l81lc/pt/tSz65D2D4wUlb3ukVToAjB4huE2MBZQUjEWKbZGS3u0yvqtop+2ggC5alFFwoG8rqJbf0j6Nl1yPV7Xr7fmW0aDh4/ZxsxPZGIjQchyYSmaUJi4KL5N3x1Ki90w8SJcxNY2dPnGaEdpjiGk3+5vUXBkss6cOszr0ERcHzeKHJ+WWuaZBGOZEJioNTP+vjOd/Dq0MnIApcq6dl1CJfnJpHXVc/iGWFTi+wCuZZc1oGjZ3KOhUYsmqlnMTcMJAgcIYTESGebhr+693OYPaMlltQbS8P982vLhaCThhs2tagykGtmm4Zbn3nNUwg3S8/itG0YSDKgICYkZuJMvfETQFWq4UZVHMMK5PrlyCfQDeemFDNaVEg0T8/itGwYSLJIxuglpAmII/WmFg3XLU/Yr0nV0vz+du8xuAWIf3Hpp/H0XddhdkdrsAcrIS21ptO0YSDJgqOCkBRjabhK7kjZ3900XDsN3epb7Nek2rPrELbuH3YVwgDwk8Mf4LK+90Jpg2mpNZ2WDQNJHnRaEJJiisUzJGRJzq6qCNyzytsUXJonXGpSzU26pzRZmt+Eh18YiCZf2ikdbEaL4qu+dr1Iy4aBJA8KYkJSTLF4xihKl/8WBVAEfAcHBe2W5OaXtiOKfOmN65bZdrdKUq3ptGwYSPLgyCAkpTjWtNZlIJ9kUJOqn8phpUShDaal1jSbU5BaSN5IJoT4IiqfZFCTqlPktR1RtypMeq3ptGwYSLKgaZqQlBKVT7IWk2qlqbhNE1j/uQV44IaFiTYf1wuveuSElMJRQkgJaaqIVKqZllKLFupkUn10zTU48uFY1ftw0/w2rluemndISBLgLCEE6a2IZAlQJXcEna1qzT7JSsF6WWcrnu17Dzd+6yeu78POVJx08zEhSYOCmBCktyKSJUB/+toZvLKmK7QWagnR7t6hVL4PQtJIcrf6hNSJoOk7SUQRiMwnmYX3QUiaoCAmTY/fhvfNAt9HY/HbV5pkB5qmSdOT5YpItQSfZfl9JJm0ximQ8ETy6woh1goh/kkI8c9CiL+I4pyE1IssVkTSDRPdvUNY2dOHO557Ayt7+tDdO+TaKcnC6X20ayK17yMNBCkzSrJFaEEshFAB/L8AbgewDMBXhBDNlzhIUk0aSigGIeyivnHdMtyzaiHUEgt1wQBMCV/CnASDfvnmJoqt7Y0A/llKeQQAhBBbAKwHwG0cSQ1ZqogURTs+TVWgCIEWTYExJRwMKbF9YASKaGzkdJpyvf3Czk3NjSjt2lLTCYS4B8BaKeXXp/79hwBuklL+ecVxDwF4CADmzp27asuWLaGumxTGxsYwcyYnSBj4DsNT+g4ndRO/+mAMps3cVoTAr10+E62auzHMlMC7J887nuPa+bMQoO9DZJw4N4Gzucnpf1/a2YoFs6PxWTdyHCb1fQeFc9mZ1atXD0gpu+w+i2I7aTc8qkaTlPK7AL4LAF1dXfK2226L4NKNZ8+ePcjKszQKvsPwlL7DXF7HIz19VRoxALS3KBi4+zc9NckjH47hkdfeQG7SqPqss1XFK2u66q6hWbnNFwrq9N9mtEjc1zUnEg290ePwZ9PPd/F3s6qkPbw6HbnbjX6HaSWKYK1RAItK/r0QwIkIzksIqYEogs+SFjndDD7UrMUpEP9EoRHvB/BrQoirABwH8ACAByM4LyGkRsK243PqsBR1NyW/NIMPNUtxCiQYoX9lKaUuhPhzAD8CoAL471LKg6HvjBBSM1Es6knqrZs0DT1OWKu7+YhkuyWl/DsAfxfFuQgh0RFmUU+ShpY0Db0ZyGJ0elLh2yWEuJIUDS1JGnqWYYWv+kNBTAhJBUnS0LNMWjuRpRlubwghqcLS0CmEo6cZotOTCAUxIYQQAO6dtxR23ooNCmJCCCEA3KPTxycNfO/1I6w1HgMUxIQQQgA4F4Ox2HngOLtBxQAFMSGEkGk2rluGu65faPsZfcXxQEFMCCFkGk1V8PVbr0KHg1as0lccORTEhBBCypg7qx2GQ2O+rFUySwIUxIQQQsqIonEI8Q/fJiElsKwfIUVYyax+cKUhBCzrR0glrGRWP/hWCQHL+hHiRNha47QyecO3Qpoeq6xfXrcv6/fY2qVcQAgJCK1M/uHbIE2PW1k/pmoQUhulVqbcpIG8bmJb/wgLgthAQUyanmZqOk9IPfjg/AR+6NE8IpfXceTDMRYHAU3ThLDpPCER0t07hM37hlFwSERWhcDjO9/B7qFTNFlPwRWGEDBVg5AoOHFuAtv6RxyFMABM6CZeHTqJvC4ZGDkFBTEhYKoGIW74iXzO5XWczU3iQkF1PE97i4JJ3cREhTW62QMjm++JCXEhbKoGIVkiSOSzV1Bjm6Zg7fJ5+PHBUxgvVLdStAIjm3H+URATUmeYV0n80uixEiS/3i2osVUTeP2x1ehs07B76JTtMc0cGMlVgJAYsFtAmVdJ/JKEseI3v750rF/a2YoZLdI26PHyKSHLwMhqmvOpCYkJtwWU1buIX5IwVqz8+rzNZ6oQOPHJBTy/91jZWO/uErhn1UK8MDDqGPTIwMhqKIgJiRCnBbRgSOwYHGX1LuJJUiq9eeXXf/8XR7HzwGjZWD87bkDpFBh4Yo2jSZ2BkdXQHkZIRFgLqF0Rgx2Do1Dti3exehcpIymV3txaIW64/grsGBytGuumlNg6pbl/5tMzXQWsFRjZ7EIYoCAmJDJcF1BFoJCR6l2siBQvSar0tnHdMtzXtQjtLQo6W1W0T/lyv3bLVYnYLGQFbkUIiQi3BdSUEveuWoidB46nNkglCQFEzUCSKr05mZFzeT0xm4UswNlDSES4mfLu71qEp9avsNUu0hKkwiL+9cNJE23UWKk0IzuNdUUI3J+SjWWS4NsiJELcIkKjDlIpTRuJm6QEEDULaQhoshvrl3aq+OOUbCyTRLJ+WUJSjCUYH1u71HUBDVu9y85E3H2DgG6YsZmIvVJZmrUiUtwEHStxFQCxO6/dZmH/P75BN0UNUBATEpJ6+07tUqTO5gz07DoUW45pkgKISDW1jEE/QtvPeVkWNjwUxISEpJ7FF5xMxFbaSFwm4iQFEJFqgozBIEI7CYVFmgHaEAgJgVvusNUAPUoamWOatAAiUiToGPQbdFfvsd3McBtLSAjq7TttpIk4DQFEzUiQMRgk6C7M2G50s4q0wTdESAjqLRidTMT1TBuhTzBZBBmDQYSr23kndBOXdbbaftbdO8Rc84DwzRASAq/c4TgEo52J+NLOVpqIm5QgYzCI0LbOq9q5QqTEs33vVf35xLkJ5prXAAUxSTVJKLdYb9+pZSIeeGINXnnkCxh4Yg0WzG6nxtHE+B2DQTeOj665BlJWC25DospPnMvrOJubpE+5BmiaJqkkSeUWG+U7pYmYWAQZg0HaEH6Um8SMFhW5SaPqs0pTtlugoAKg/+jH6FryqdBzI4v+52w8BWk6kphWQcFIGo2fMRhEaAcxZbvFQ4wXTPzpDwZhhNgwJ2nzHTXpvnvSlNQrrSIJZm/SXNRzzPlpQxjElN3ZpuHSztaqYy3GQ/qMs1zrnBoxSR1xpwxleedNkkmSx1wQU/aC2e24r2sOtvaPQEFRE66kMlXKj6k567XO03vnpGnxYy4L40dKotmbZJswYy5un2mlKXtmm4axvI68bl/b3Dq2/+jH+NMfDGLcwb984pMLeH7vMV+bj6zXOqcgJqnDrdziPasW4plXD9esWWR9502SR61jLmot2kugt2kK/sc/HPV1vc42DV1LPgXDZcP8/V+8j50HjvvafGS91jntbCSVOKVrACKUH6mWEpL0JZMw1Fq2NKjP1Gmc6oaJ7t4hrOzpwx3PvYGVPX3o7h2CbpRvDIJez82/vOH6hdgxeNx3nEcj8vXrSbrvnjQtdpGfALCypy+QZlGpBQTZecfl18tiegZxphZtL4gW7TVO/ZjFa9XanfzLX735SvS+dTyQqTmIrzptcJYTR9IgEErTNY58OObbj+S2OPntMhS1LznJATskPoJ0trLm5IVJw/dYdxunj61d6kvA1uqjdUqVyuX1wJuPLNc6z8ZTkEhJq0AIolm4LU5+dt5x+JIZJNa8eI05uzlZMKojkoHyse41Tu/rWuRLwIb10VbmN1ubj639I5ioME8vuawDbZrzOpPFfH0KYlJFWgWCX83CjxD12nlHHcXJILHmxkvbs5uTqiKgQsIokY+VY91rnALwJWDj6Ee9cd0yvPn+Rzh8aqzs70fP5NCz61Ci15qoSa56QxpC2nuQ+qm56zc4xq3gQdRRnI3sM0ySg92Yc5qThilhSqBNFY5j3WucLpnT6TsIKkhNdT8BjHndxPtnxqv+PqHLVKw1UcItNikj7fl6fvxIUQjRqDWErKdnkNpxm5MSwJplc/Ho7/667Vj3M079BkHldRN/9BtL8PDqqzGW1x19tH7bIKZ9rYkSCmIyTS6v48KkkQmB4OZHmvZP7R/GhH7xWYMK0coFTDdN/N7yeXh0zTU13a/dgtneomDt8nmBz0eyg9smDQD63j2N79zzuZqbPFib14dXX413T57HtfNn4XKf2QGVAZ2lbRCbPTc4CBTEpGqiFQwTqiLKkvHD+IKShm6YMCVQKCn4owrgnlULA6VCWAvYo2uuQffLB7F76BT6Dp3G7qFTNQW3VS6YE7qJSd3Ejw+eqvmcJP10tmm4fcU8vPTWCdvPNUVx1R69rES1pDdt2XcMr//qQ4yevYAWVYFuSty98gpcMTGJCwW17PpOcQ5x+J2BdGR7VMIZTaoS9U0JQEqoirPvKc307DqE7QMjMEr6rLZoCgwTGP54PLBv6tm+9/Cjg6dCF6Mv7TO8ZtlcaKLoAxwvmJkqcE+Cs+nO5XAIIQjkTrGLeXAr1OHkn84bwJEz45g05PR3dgyOwqZ1MQDnOIcoe3n7LUySREJtF4QQ9wL4vwFcC+BGKWV/FDdF6odTtK4hgVZFYtsf/1ssmdOZmp2lF07PO1EwsXnfMF46MApDwlH7rNxtxxHtrBsmXn77BCqtdoygbl5md7TiD2++MrQ7pZIw6U2V5HUJicblBqc12wMIrxEPAdgA4OcR3AtpAG7RupMG8F9fP+Ka05c23J4XuKh9btk/jCd7h6b/7rTbPvHJhcijnbtfPlglhMOek6SfjeuW4f4bFkeiPVp4ResDzulNdgigar3wU4bST0tGN9Ke7RFqWy2lfBcAhHBe2Eiy8QoEefXgKcye0ZL4HaVfvJ7XIq9LbN43gvMTOp7+8go82/ee7W67YMhIA05yeR27h045fq6bZlMFsZCLxFFZym96U6Uf1wkhBO5euRAvHhgNXIYyjG837RHYQjoZ9YOcRIg9AL7hZpoWQjwE4CEAmDt37qotW7aEvm4SGBsbw8yZyf2B/XDi3ATO5iZhOowFRQhcO3+Wo48qLPV+h17PW4mY+r92ZjdFCFzS0YJPxgtl51OEwKWdrVgw215omrKoZWuqUvZeJ3UTv/pgzPHeLuloxaJLZ1T9PQvjsNE06zu0mw+V49fPnFGEwPxO4FOzZ5WNb8B+rNvdg4Xb3LHDlMC7J8/b3l/c65dfVq9ePSCl7LL7zHPbIYT4ewB2+ROPSyl7/d6ElPK7AL4LAF1dXfK2227z+9VEs2fPHqT9WXTDxDdeeNsxKrOzVcUra7pi21HW+x1aUaJb+0egCGHbL9Uvna0qXnr4Fjy/91hZesj9XYvwxzY+Zq8I1VxexyM2jSuAYiWlwSe+iNkdrVWfZWEcNpoo3mEaI3ZL54Pb+M3ldTy+8x28evBUVVlKVQBfvflKLJ99Zvod+i2V2907NKVxX4y2ntEicV/XnECWuJ9Nn6c6Avvh1cm26HmOFCnl79TjRkjj0FQFT991HXYPnbIVAFnL6as08X3v9SPYMTiKvB7cOmRIiQWXzPBtMvQKKHHLJ76/a5GtECaNJ6312QH/Ju/ONg1/de/nMHtGS0nuvMTtK+Zh053LMbujFXv27Jk+Ps6uTnakuTtTOrZsJHbiyulLMlaAyFNTi8LmfSOBvl/5bryK0ftddOwWlPtTsqA0K2mO2LXw00zBS2ibstgFTRHA5jePodKtHFVXp1ruzY6kWDDCpi/dBeA5AJ8G8L+EEG9JKX8vkjsjdSfNO8owaKqCb234LAAxpRnbB6WoonisptT2bvwuOllu95ZFmrFhR6XQtiwCn86dx0M//hncjEt+uzrppsTMGt6bnw1F0iwYYaOmdwLYGdG9kAbT7ALgqfXL0aIKPP/mcFlVMeCi9hvm3QQt6ZfFdm9ZJO0RuxZhtEPLIvBnS6WrEAb8dXUCAMM0ceszrwUWkH6eI2kWjOZZZYlvmlUAVJesPAlNUcq0X01Van43zWj+bwbSXjM5rHboZBGwo10Ttl2dgKIlTi9JB9TNYrqeXwHp9zmSaMHgzCdNh9eOeXZHK/76gevxdP66yC0DzWr+zzJp32CF1Q7dLAKVrF0xv2qslzaduOU7P636jl8B6fc5kmjBSPYIISRCgu7847AMNLv5P6ukdYPlqh3u96cd+i2S06oKPH3XdY5a9lheR6uqoGBUpxN6CcggWm4SLRhcAUjTEIdfqFa/WrOa/7NKozdYtY7D0+cnquIhLHRT+tIOSy0CTqgC+MqNi13vLYyADKLlJtGCQUFMmoKo/UJRRF0mJXWCREdUGyy/YyPsOJzZpkUS7nNIIQAAIABJREFUtWxp/kruCDSBsoAtVQAP3rTY0zoQRkAGFeJJs2Bw9pOmIGq/UBjtOmmpEyQ5BB0bYa08Y3kdmlIMjKpEU4qfX+7jvi2LwE9fO4Mfr+nCzDYNH/5r8Y7surc5bTRqFZBBhXijLRiVUBCT2IlD8wt6zqj8Qrm8jqNncvjhvmFMGuXnizqohDQfQcZGFFaeubPaoQgB2NRRVxUR2F+qCExvaC+3+W7lRqNgmPj96+ZPV+YqDdx69+R5XDt/lu157KhFiCfFRURBTGIjDs2v1nOG9QuVXlcRokoIW0QZVEKai6BjIworT2ebhq/cuLhu/lK7jcZLb53AK2+fxFdvXoxv3r4U3959uKY1I0ott95uI854EhtxaH5hzhnGL1R6XTeiDCohzYWfsTF3Vvv0/4/KymON/y37hyGEgGlK3HX9Qtt5EUZAueUbG1Jia/8I3nz/Ixz7aDzUmhFGy22U24iCmMRCHJpf2HP63TFXLjYfnJ/A5n3DKDhowRZxBJWQ5sG93KOJ773+PnYMjpYJiPu6FmH7QLk226Yp2HD9wkDzq9i2UMKYaiO4Zf8wVAXovmM5NFUpE1Dq1PXvWbkQT61f7ltAeeUbTxRMHD41VvX3yvkdp7baKLcRBTGJhTg0v6jO6bRjtvNffebTnTjy4birEG7XBCAE7utahEfXXIMjH465drCpR+oEI7LTh9vYuPKyDuw8MFolIO5ZtRD3dS3C1v0j0M1iVSrDNLF9YAQtqvClyfXsOoTN+4ZROsRNCWx+cxiKENi0fsV0m8TSTfDmfcMYHP4Yux65teoaduPPb76xHaoQOPHJBTy/91hs2moj3UacoSQW4tD84tYm7XbDdjv0SkwpsOH6K2BKiRu/9RPPRSLO1AlGZKcbu7Fx1/VXYMfgcVsB8cLAKAaeWIOCYWLH4OhUWUgAkL40OUv42OUSG7KoGT+8+mpHk/LhU2N4sndoqmlKke7eIdvxZ200tu4fxkTAlqOGlPj+L47abka8ntEvjXQbURCTWIhD84tTmwxSL7eSSaO4IEJKGBKei0QcQSXWOsqI7HRjNzZOn59A71snHAXE0TO5KUEdPIr/9PmJqahpexQh8O7J81AV52O2D47i8S8VBe2JcxOu48/aaDy/9xgqjUyW5n/so/Gq+V3cjFR3RotSW22k24hbZBIbG9ctw31di9DeoqCzVUX7lMAMo/nFcU7g4m64VgxTVi0s1iKRy+u237FM5LUsILphort3CCt7+nDHc2/g3ZPn8ZcvvoMf7huu6mLjdR8keZSODS8BAcBx7FqanBNzZ7XDlM7aqSklrp0/y9WkrE1dI5fXcTY36Tr+rI3G4MY1+PLnF6BNK5/HL/3ZLbbz+2u3XFXzM/rF2ujPaCkXizNalKpGFVFDjZjERhxJ85XnnNmmYSyvI6+boUyvYfxXbsRl0qrUfE0psWNw1HFRZUR2evGyBC2Z01mzJmed2671pyqAB25YjMtnteOelQuxed+w7TlMYFpzd6Jy/F1srFLtS7ZbM3J5vS7aaqMqblEQk9gJk07gFHTUpin4H/9wNDJfqF//VXuLgkndhF+ZHYdJy8mM7mZWZ0R2unETEJqqhHLZbFy3rBicta9cGJuyuMHTDRNPrV+OweGPq2ImStsauo0vp/HntDZU/r1eQY6NqrhFQUzqQtAoXq+go6h9obpRFK6FksYvCoBfmzsTxz4ah6oUF7/7uxbBlBLbB0bLFgRVEdM+You4iiK4BZVoioCqoMxfmJZ2fMQZLwERRpPTVAU9X14xHW09OTUHJIDtA6PTkdO7HrkVT/YOYfvgKDQhYAK4v+QanW0aLu1sxYwWGYuwrKe2Wu+KW5yZJFZqjeJ1E7SPrV0aeZpBz65D2D4wMu1zA4DWFgWrrrwUz973eQAXa+bqhglFiLIF4d5ViwBIvDAwGvsi4WZG1xSBDSuvwIsHjieimD2JFicBEVaTy+V17Bg8Pi2ELSrn1Lc2fBaPf2mZ4zUWzG7HfV1zYhGWSasPHSXZeAqSWGrRXL3y+e7rWhRpmoHT9SYKJjbvG8FLb52AUbGBcFoQ/uL2a2NfJLzMdJvWr3BdLEl2sRPUfqxRQdsIus2vuIVlUupDRwlnKImNWhPkvRYFAJEGbnhV/BmfUhMqNxB2C0K9FolKM50ijDLNI4uLVVappfiK23eszy7rbMWzfe/5skZFnbrD8RcMCmISG7UmyHstCkvmdEYauOE3YjpJjRkqtfKjQ/14eDXzhNNELW4bt+8AKPtsomACcM9tLxXo9QiGIvbw7ZLYqHWX7SdCMsrADafr2RHU9B13qUlL8xiuPQWaNIha3DZu37H+d+lnlVibyUfXXFOlLd/XtQj3rFpYlzgHUg4FMYmNMCkHXoI26sCN0uspAMYdBLLTBqJS4LLUJHGjFreN23e27B8GIHxVhhMQeHznEH5y+HSZ0N4+UIy/GHhiTc2mcmvfzVrnweAbIrFSq+bqV9BG5YuqvN73Xn8fOw+Mem4gnASuleLEUpPEjlrcNm7fKZap9GcWuVAwsOudkzZ/v7gJ8DunKsf/ny6dxDN//TO8f2acG9AAUBCTWAmrudY76MO63lPrl6NFFRXF9xfiqzdfiVxen34GO1Ph1v4R26IfSfIxk8ZSi9vG7TtuZSqD4Mf1UqrtPvPq4bLxL6WcLvrBDah/uBqQupC2KMrSDcSJTy7g+78o9oLtfev49C7/0TXXOKY9OcFSkwSozW3j9J32FgVLLuvAe6e9O4V54Ra7YWf9KRjeVeZq3YA2k3k7209HSEg62zQ8v/cYXhwcRV6XZbv8cxcKrmlPdrDUJAGKQuarN1+JgiHx4gH/wVF2rp4ll3Xg6Jmc77KrTnjFbthZf/wSZAPajPEVFMSk7qRpp3tufBJ/u/eYrZl599ApSKcmC4pAi4KyutVMBSF2QubulVfga7dchQWXzPAcG3ZNT2595rWqFohB6GhVYXpsAsK0CQWCbUCbsZUnVwRSN5x2uqtnN/rOnHl85zuOmoamCKxZNg8/Oniqyrx4z6pFUIS/ILU0bUxIOOyEzM4Dx9EyJWD9Yrl6jnw4FtgqY2H1+f36rZ/xHHteRW+8ruN3A5rL6/jhvmFMVvQUtSLDnczbaZ9D6btjklqcdrq/3pW8JFjdMPFk70HseueU8zGmiU13LsfsGS2OXXHcgtSa0QTXzNRaac6NWtp3VmrAfsaa23VURUBTAE1RoAgDS+fNxNEzFxulBMlFLqZAOQSxmbLKvJ2VOURBTOqC2yJ0NmeURSLXev6wO+LKaNAdg6Oux9++Yj5md7S6RoW7Bak1owkuTUStZdVaac4NpwAup25gfjVgv9extN2yCm8P/lbN725mmwYn67duFj8vJStziIKY1AUv05afRahycufy+lRE81HsGByteUdcuasuGCYME2WdmCpRFYFNdy4vu6coGk0wxanxxKVlRV3P2cIugMutG5jXM9jNs9PnJ/DommuqrlN6ztIKb7VmSYzl9el3XommCIzldVxecp9ZmUPpuEuSerxMaG6LkN3CeNWcDhz5MFdsXD41D2vdET/ZexA7Bkd9R4O2agL3dy32XVDfjqDaUdp9YGkiLi0rrub2UXUDc5tnLaoyPcbf/OYX8VFusuax6DaW585qh+oiiEvXiTgsDI2CM5rUBbdF6NJO1XVC2y2MVtEAO7wCOyyKfuAhbN434niMHfesXAhAYlt/7ZWz/GpHWfGBpYW4taw4m9uH7QbmNs8mDfsOZEHwM5aDbFbisjA0As5kUjc2rluG+7oWob1FQWerivapybVgtvOEsRZGr2YMleR1icd3vgPdcP5ez65D2D7gXwi3aQoevHExHv/SMmzrH626J2uxzuV1z3NZC86MlvIpOKNFwf0lC07p4pibNJDXTWzrH0HPrkO+75v4x9Ky7LC0rDBY2uvAE2vwyiNfwMATa7Bp/YqGb6r8zrPiJtffGK88/zdeeBtb9w97jmWndaJys+J3DqWB9NwpST1OJrQ9e/Y4fidM2sSrB09h9owW2927U35wKaoo3rM2Ff15/9RiMPzxeCQmMS/tKEs+sLRQLy0rSZXmcnkd/Uc/huozeSGvm3h85zv4q3s/57mBKNWC7XKQ7caytU48vPpqvHvyPK6dPwuXO7z3OC0M9YSzmNSdIItQLekZFhNTu/eHV19dNZG7Xz7oKoTbNAUP3LDI1u8W1WLtVYc7Sz6wtBCXHzeJlApJVRGOHcfseHXopOMmt5RSi44TlWM5iDsm6i5sjYKm6SYjl9dx5MOxwKalRtHZpuHulQvRptU2VPO6iVu+81N09w7h3Pgkjnw4hg/OT2D3kHN+MADcvXIhNq5bNr1pqExJitIkZncNIFs+sDTh1zSaZPzM81IhOT5pBDr/hC493TB+zd2VY7kWd4zTHEoL6bxrEpg0Bv1Y97x9YGQ6yV9TBDRFYMmcDrx/JgdNUaCbJnSjPGeylIIh8fzeY3h+7zDaWxQUTAnTdF4c1l03D9/acJ3rvdXDJNZM2lmSSLOW5Xeee5Ws7GhRYEpg/iXteP/MuO0xXlYZP26lyrHcrO6Y7D0RsSWNie/WPZeWu1MVYMPKK/CtDZ+1bcfmtPsunkIi57HzVxWBp++qFsKVKRf1Wqyz4gNLI0ny4/rF7zx3E5IdrSr+5g9WomvJpwAA1z/146qSk4C3VcbLrdSmiaqx7OWOOXomhxmtaqo2R37IzpMQR+LcZUaZ31p6LgC295zXJV48cByPf2lZ2UJpTeYwhenbp0zLsztap//mpWHEvVinWTsj8eA054LMc6/exl1LPjV97FduXOxqlcnldRw9k5v67sXzuLVtXLt8Hp6+67pA7piJgoG7/vMvynKak2zRCwJndBMQR9BP1Kbu7t6hsnPdvmJeoHsujbS85Ts/RcHJTl2CIsoXjiWXdeCbty8tOyYploS4BD4LhaQHrzkXZJ4HcXs4WWW+eftSbHzpHfzgzeHpefSN6wy89tIQuu8o3pPdd+93qPClGyaeefUwCg4ph4YEDENGktOcNDjzmoA4gn6iFFAnzk1UnevVoZMoOFiRDSkxs03DkQ/HqgTI5bPa8aDNDt6Oyldy7KNxfHv34en7z7K/Ko0xA82O15wLOs/9uj2crDLdvUPYXCKEAUBCYvO+YSgC0/nRfi061vP5TZLIwjy04IxrAqKO8nWKhgxS0KL0XGdzk1XnKvbxlWivuOf2FgVLLuvArc+8hjueewMre/rQ3TtUVrijMupVFVNF8D2ovP+4izs0EhYKSRd+5lzQeR60uEhpZLLVrtDO8GSYsqroh1dUc62Fe9I+Dy0oiJuEKFMyohRQbse2aUVfUuk9L7msA0fP5FwFSOkC89LDt+DeroWQJQ0cVAFoDnK59P7dNAzdNFObPhTlRorUB79zrpZ5Xkvqz+nzE66bW6WGdcDp+dzIShpfuvV54psog36iNHW7HWsCePqu6/D0Xdfh9PkJzGzTcOszryGvVzcNtzNRdbZpeH7vMfS+daLM3KWpwtGHXHr/Tn40ANAN4JlXD6fSlMtCIenD75yrV3Df3FntMFxsyGYN60DQwj3tKSxl6US6VhASmigS36M0dXe2abi0s9X1XNY9Wy3S7LDTxJ00v7wuIYRAe4VabHf/loZRWf7PkDK1plwWCkkfQedc3AUuOts0fOXGxbZlMVVF4IEbgq8Dds+nKsLxGvdnKI2PgpjURJSm7gWz232dK6gAcTN3tWsK1q6Y73lNTVXw2Nqltma4tJpys1Qsv5lIWsWvjeuW4cGbyoWxgMCDNy6u6Z7snu/BGxfjwZsWT/+tTRP48ucXYPCJ30lEs4yo4IwjNRG1CczPuYJWmrqssxUXHEKvDSnLzN5O19SNYoF7u4IGQDym3HqkFLFQSPpIWk65piro+fJ1+Ivbr53OIz793gH8+era0omi6qucRrL3RKSuRJnf6udcQQTIs33vQQgByHIhqgqUaX5u1+zZdQivHnSuSx2lKbeeKUVJW9SJf0rnSRLywDvbNCy/YjYA4MNfRXO+MH2V0whnHkkVfgWI5R+2DSgRAo+uucbzWl71eNs1EakptxHFQ7K+wGUV5oFnC/5iJJV4BaN4+Yc/yk16XsMrpWLtivmRmXKZUkSCEGceeNo6tGWBUFt5IcR/BHAHgEkA/wLga1LKT6K4MULCUEtkcKWZb+6sdkw4FBiwmkNEpX0wpYj4Ja6Kb9SyG0fYt9sHYIWU8rMA3gPwzfC3REg1QXfpQSKDdcNEd+8QVvb0VVXrknDIbZTBch69YEoR8UtcFd9Yba1xhNKIpZQ/LvnnXgD3hLsdQsoJs0v3G9jl5Js9d6GAdk3BuI1W3KYpkWqp7D1M/BLHpi3LddXTQJRv9n8HsDXC8xESKoDJT2CX2wK0e+hUWWnMUky4VwWrBaYUET/EsWmja6SxCKeFZvoAIf4ewDybjx6XUvZOHfM4gC4AG6TDCYUQDwF4CADmzp27asuWLWHuOzGMjY1h5kwO0DA4vUNTAu+ePA/TZkgpQuDa+bMQpDytKYsatqYq09+b1E386oMxx2vMmtGC8xcKZZ8rQuDSzlYsmB2PudjuPr3gOAxP2t7hiXMTOFsSdBhmTPqZa4D3uEzbO6wnq1evHpBSdtl95imIvRBC/BGAPwHwRSnluJ/vdHV1yf7+/lDXTQp79uzBbbfd1ujbSDVO7/DIh2O447k3kJusLsrR2arilUe+4GuX7mbezusmVvb02aYotbcoePObX8Szfe/56qfaSDgOw5PGdxhlHnF375Ctln3PqkVQBHy5h9L4DuuFEMJREIeNml4L4N8D+C2/QpgQvwTxhbktSE7mbVNKKELYNiK3zHyzO1qnzdtW9aAlczoTJYRJ8xJlHriTa8SUEtv6R+ua395shPUR/ycAbQD6hBAAsFdK+Seh74oQ+POFeQVzufmAN785jBZNqWpEriqizDerGyaeefXw9DUKhonfv24+Nt25HLM7WmN/D4TYEXVVLbuYCgC2FiMGcUVL2Kjpq6O6EULs8Apg8grmcgtCMSRg2EREqwJ4bO3Saa3X7hovvXUCr7x9El+9eXHizNQk28Sd71uqZR/5cIxBXHWAWxmSaNwin/2kXNTS53TSkHh85zv4q3s/h7xuOpa5NKTEVproSEiCarb1LIXK/Pb6wG08SQV2JS39FDZwKuzhxatDJ9Gz65BnmcsJHyUoWTKQ2OFWSMaJepdCZcvM+sC3SFKL3926ZcZ+fu8xOHQzrGJCL2q7D6++2lOjdjLRsWQgcaMWzbYR+b7Mb48frgYktfjdrWuqUubz9YsqBMbyOh64YRHaNWet2MlEx5KBxIlaNdtGmIot99DAE2vwyiNfwMATa7Bp/QpuJiOEb5Kkmo3rluG+rkVob1HQ2aqifSqiunK37mVitsNa2DauW4b7b1gM1ebrTiY6dlMibtRaLzqIqThql4hXxzNSO3yjJNX47U/spkmoikCLUjRHW1SWC9y0fgUeXXMNul8+iN1Dp6Ap7ia6E59csO+FDG8TYhKavZN4CaPZepmK43CJcEzGC98oyQRehQ3ccpKtykFePrDZHa346weux9M+FqXv/+L9wAstfcrNQ5h60V6bzyijqjkm6wMFMWka3DQJy4/sZ9fvJfRzeR07Bo87fr7h+oW2569nWgppPGGDoOzGYdRdlDgm6wMFMWkavDSJqMoFukW2agrwtVuWVP2dbeiaD79ulSBEGVXNMVk/aFsgmcUpWCXuoBM3/5+mKlhwyYyqv8fV7J0knyDj0SsAK8qoao7J+sHtDMkcjfZr1eL/YwUj4obfMR1lr2KOyfpBjZjEQiOrSSUhf9dvWpUFKxgRN4KM6aBjzwmOyfrBN0kipZHaaC6v4+iZXCL8WrX4/1jBiNgR1Fcbpe85yJjM5XVM6iZyeZ1COiB8WyRSGhFlWSr8FQHbBg1AY7rFBAkAiyN4h6SfWgOwogg+9DMmS+ff/7GsgEd6+pjiFBC+JRIZjaomVSr8K69dSlr8WqxgREpJgq/WbUyWzj9TSpZyrQEKYhIZjYiydBL+1dcH7l1ln78bN+y+RMKQZF8tS7lGA7fcJDIasXN3M9uVIQSAYLWmw9Lo6G2SHZIaP9CIblBZhIKYREaUqRN+cRP+pRimxAsDI/iL2+tXhIBViUhUBI0fqFdt6CSYzbMAt+UkUqJKnfCLk9nOjnoWIaDJjsSBV/yAbpjo7h3Cyp4+3PHcG1jZ04fu3iHohrvrJsz9JNVsnib4lkikNCLyt9RspwAYd/AXR71Dd9M6aLIjjaARVpiy+ScE2jSB31s+D4+uuSaW62URCmISC1HVbfZDpfD/3uvvY+eB0djM4358vzTZkXrTqNrQ1vx7dM01ePlHPwEA9B06jd1DpxgT4RO+HZIZLOH/1PrlsZrH/VQ5osmO1JtG14Z+tu89nL9QaGhFu7TC1YA0DMu06yPWKhBxmseDaB1JjXQl2aSRVhhrXvz5teXXZ6cmf/DNkLpTadr9s6WT+FnvUOQmrDjM40F8v6yURepJI7IWLPxo44yJcIarAqk7lQElppQ1B5TUK03Dohato57+ctLcNMoKw5iIcFAQk7oSVUBJo4plNFLrIMSLRllhrHmh5I6U/Z3zwh8M1iJ1JaqAkka2Oqx3rjQhQWlEvfKN65bh0s5Wzosa4DaF1JUoTFiNStOwoO+XkGo0VcGC2e0YeOILnBcBoUZM6koUaT2NTtOwYJckQqrhvAgO3xSpO5UBJYowApmwGBhCCMkSFMSk7lSado8O9ePh1f6jpRkwRQjJElyxSMOwTFjDNXQnZLEMQkhWoCAm/3979xoi51kFcPx/stvcNqgVJZKLvUBpDcU2dZFqURIrUrE0inhDBQUpgpcqSrFKLVoEwSL6QQSJ+sXiqrFaqeINDfipZJtUbE2FkupmTaURtTUp3XS7xw8zU8YwMzsz7yTPvjP/36fMTJg9HHbnzPO85z1PLdkwJWlc+MmlWnNYhqS6s2taY+f00jLHTp7yzF9JteCKWGOj1LQtSarCTyeNjU7Ttn5waIFP/+iPro4lrVkWYo2F1rSt9tuZAJ5ZTn764AmuufM33HHvQyw/t9LlHSSpDAuxxkKvaVuAh5RLWrMsxBoLvaZttbRmUQ+yTW3jl6RzzWYtjYVu07bO1u8h5TZ+STpfLMQaG62pWnOHFlharjaLur3xa6n53A/njwPwhX39j+OUpNX41V5jozVt6/Dtb+KtV29j45AnPHVr/Bpma1uSVmMh1tiZ2TDNXe+4infN7hzqkPK1csyipMng1rTGUpVZ1B6zKOl8ckWssTbMIeWtxq9Nq2xt21EtaRRcEUsd9Dpm0Y5qSaNkIZY66LS1DbDwr6fZ/4fH+MmRRTuqJY2EhVjqYWbDNC9/8ebnV8BTAU93uE+51VF96w1XeC6ypIG4jyZ10H79t/2e4k5FuMWOaknD8Ku71KbT9d9nn1thlemZQO+O6tNLywN3b0uaDH4iSG06TdTqx6bmfcpnF1kbuyStxkIsNbUmai0t939U4ub1U6y0dVSfzVGZklZTqRBHxJ3APmAFeAL4QGaeGEVg0vnWmqjVz0p40wXreNvu7XzodZd23W7uVtht7JLUrure2Fcy85WZeTVwH/D5EcQkFdFrotbUumDDdPzfuMwv7ruy57AQR2VK6kelr+OZ+VTbwxmgj5YWaW3qdpRi6/rvoOMyHZUpqR+Vu0Ui4ksRcRx4L66IVXO337iLd3Y5LGLQcZn9jsqUNNkis/ciNiJ+C7ysw0ufy8x72/7fbcDGzLyjy/vcDNwMsHXr1lfNzc0NHfRacurUKbZs6X3IvHpbizlcyUbH8/TUOrrsLvftxJPP8O/TZ55/fOHMera9cLSr4bWYw7oxh9WZw+727t37QGbOdnpt1ULcr4i4CPh5Zq7aCjo7O5vz8/Mj+bmlHTx4kD179pQOo9YmIYfn+j7iScjhuWYOqzOH3UVE10JcaWs6Ii5re3gT8EiV95PG1TCnQEmaDFU/Fb4cEZfTuH3pb8CHq4ckSdLkqNo1/fZRBSJJ0iRyxp4kSQVZiCVJKshCLElSQRZiSZIKshBLklSQhViSpIIsxJIkFWQhliSpIAuxJEkFWYglSSrIQixJUkEWYkmSChrZecQD/dCIkzROaxoHLwH+WTqImjOH1ZnD6sxhdeawu4sy86WdXihSiMdJRMx3O+xZ/TGH1ZnD6sxhdeZwOG5NS5JUkIVYkqSCLMTVfat0AGPAHFZnDqszh9WZwyF4jViSpIJcEUuSVJCFuIKIuCEi/hIRj0bEZ0rHUzcRsTMifh8RRyPi4Yi4pXRMdRQRUxFxJCLuKx1LHUXEiyLiQEQ80vxdfE3pmOomIj7Z/Bt+KCK+HxEbS8dUJxbiIUXEFPAN4M3ALuA9EbGrbFS1swx8KjNfAVwLfMQcDuUW4GjpIGrs68AvM/MK4CrM5UAiYjvwcWA2M68EpoB3l42qXizEw3s18GhmHsvMM8AcsK9wTLWSmY9n5uHmv/9L4wNwe9mo6iUidgBvAfaXjqWOIuIFwOuBbwNk5pnM/E/ZqGppGtgUEdPAZuBE4XhqxUI8vO3A8bbHi1hEhhYRFwO7gfvLRlI7XwNuBVZKB1JTlwInge82t/f3R8RM6aDqJDP/DtwFLACPA09m5q/LRlUvFuLhRYfnbEEfQkRsAX4MfCIznyodT11ExI3AE5n5QOlYamwauAb4ZmbuBk4D9nsMICIupLEbeAmwDZiJiPeVjapeLMTDWwR2tj3egdsxA4uIC2gU4bsz857S8dTMdcBNEfFXGpdG3hAR3ysbUu0sAouZ2dqJOUCjMKt/bwQey8yTmfkscA/w2sIx1YqFeHiHgMsi4pKIWE+jOeFnhWOqlYgIGtfmjmbmV0vHUzeZeVtm7sjMi2n8/v0uM12JDCAz/wHRfGOqAAAAmklEQVQcj4jLm09dD/y5YEh1tABcGxGbm3/T12PD20CmSwdQV5m5HBEfBX5Fo0vwO5n5cOGw6uY64P3AnyLiweZzn83MXxSMSZPnY8DdzS/Ux4APFo6nVjLz/og4ABymcSfEEZywNRAna0mSVJBb05IkFWQhliSpIAuxJEkFWYglSSrIQixJUkEWYkmSCrIQS5JUkIVYkqSC/gfRw7LYGUAfqAAAAABJRU5ErkJggg==\n",
"text/plain": [
"![](\"https://img.icons8.com/dusk/64/000000/zoom-to-extents.png\")
Expectation Maximization (EM)\n",
"---\n",
"* Probabilistic method for **soft clustering**.\n",
" * The \"soft\" version of K-means\n",
"* Assumes a probabilistic model of clusters that allows computing $Pr(c_j|x)$ for each cluster $c_j$ for a given example $x$.\n",
" * If we had known for each data instance from what distribution it came from, we could have used a parametric estimation.\n",
"* We introduce *unobservable (latent)* variables which indicate source distribution.\n",
"* We run an iterative process:\n",
" * Estimate latent variables from the data and the *current* estimation of distribution parameters.\n",
" * Use current values of latent variables to refine parameter estimation."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### ![](\"https://img.icons8.com/cotton/50/000000/tie.png\")
Formalization\n",
"---\n",
"* Log likelihood for a mixture model (under the *i.i.d* assumption): $$ \\mathcal{L}(X|\\Theta) = \\log \\prod_{i=1}^n Pr(x_i|\\Theta) = \\sum_{i=1}^n\\log \\sum_{j=1}^k Pr(x_i|C_j;\\Theta)Pr(C_j;\\Theta)$$\n",
"* Assume *latent* variables $z$, which when known make the optimization simpler:\n",
" * **Complete** likelihood, $ \\mathcal{L}_c (X,Z|\\Theta) $, in terms of $x$ and $z$\n",
" * **Incomplete** likelihood, $\\mathcal{L}(X|\\Theta)$, in terms of $x$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* The **incomplete** likelihood: $$\\mathcal{L}(X|\\Theta) = \\log \\left(\\prod_{i=1}^n Pr(x_i \\mid \\theta)\\right) = \\sum_{i=1}^n \\log \\left(\\sum_{j=1}^K Pr(x_i \\mid z_i=j, \\theta) Pr(z_i=j \\mid \\theta) \\right) $$\n",
" * Calculating the log-likelihood of the observed data is **hard**!\n",
" * This is because of the sum inside of the logarithm - the derivative w.r.t. to the parameters and equivalenting to 0, creating an equation with **no closed-form solution**.\n",
"* The **complete** likelihood: $$ \\mathcal{L}_c (X,Z|\\Theta) = \\log\\left(\\prod_{i=1}^n Pr(x_i, z_i \\mid \\theta) \\right) = \\sum_{i=1}^n \\log(Pr(x_i \\mid z_i=j, \\theta) Pr(z_i=j \\mid \\theta)) $$\n",
" * Much easier to calculate! But...\n",
"* However, $z$ is *latent*, so we **can't compute** $ \\mathcal{L}_c (X,Z|\\Theta) $\n",
" * But we can compute its **conditional expected value**, given $X$ and old $\\theta^t$: $$ Q(\\Theta;\\Theta^t) = \\mathbb{E}_Z [\\mathcal{L}_c(X,Z|\\Theta)|X,\\Theta^t)] = \\sum_Z Pr(Z|X,\\Theta^t) \\log Pr(X,Z;\\Theta) $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* From a computation viewpoint:\n",
" * **The E-Step**: computes the **posterior probability** $Pr(Z|X,\\Theta^t)$ using the *current* estimates \n",
" * $Pr(z_i=j|x_i, \\Theta)$ - probability point $i$ belongs to model $j$\n",
" * **The M-Step**: updates the **parameter estimates** to get $\\Theta^{t+1}$ by maximizing $Q(\\Theta; \\Theta^t)$\n",
" \n",
"* The EM Algorithm requires an **initial guess** $\\Theta^0$ for the parameters.\n",
"* Each iteration of E-step and M-step is **guaranteed to increase the log-likelihood** of the observed data, $\\log Pr(X|\\Theta)$ until a *local* maximum is reached."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### ![](\"https://img.icons8.com/office/80/000000/baby-footprints-path.png\")
The Steps\n",
"---\n",
"* Iterate the two steps:\n",
" * **E-step**: Estimate $Z$ given $X$ and current $\\Theta$\n",
" * $Q(\\Theta|\\Theta^t) = \\mathbb{E}[\\mathcal{L}(X,Z|\\Theta)|X, \\Theta^t]$\n",
" * **M-step**: Find new $\\Theta$ given $Z,X$ and old $\\Theta$\n",
" * $\\Theta^{t+1} = argmax_{\\Theta} Q(\\Theta;\\Theta^t)$\n",
"* An increase in $Q$ increases the incomplete likelihood $$ \\mathcal{L}(X|\\Theta^t) \\geq Q(\\Theta|\\Theta^t) $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"#### Proof: \n",
"---\n",
"* Denote the prior of $z$ for the $i^{th}$ sample: $q_i (z)$\n",
"$$ \\mathcal{L}(X|\\theta) = \\sum_{i=1}^n \\log \\big[ \\sum_{z_l} p(x_i, z_l|\\theta) \\big] = \\sum_{i=1}^n\\log\\big[ \\sum_{z_l} q_i(z_l) \\frac{p(x_i, z_l|\\theta)}{q_i(z_l)} \\big] = \\sum_{i=1}^n\\log\\big[ \\mathbb{E}_{z \\sim q_i}( \\frac{p(x_i, z|\\theta)}{q_i(z)}) \\big]$$\n",
"* Since $\\log$ is **concave** we use **Jensen's inequality** to derive the following lower bound: $$ \\mathcal{L}(X|\\theta) = \\sum_{i=1}^n\\log\\big[ \\mathbb{E}_{z \\sim q_i}( \\frac{p(x_i, z|\\theta)}{q_i(z)}) \\big] \\geq \\sum_{i=1}^n \\mathbb{E}_{z \\sim q_i}\\log\\big[\\frac{p(x_i, z|\\theta)}{q_i(z)}\\big] = \\mathcal{F}(\\theta, \\{q_i(\\cdot)\\}_{i=1}^n) $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"* This is **true for every choice of $q$**, and in particular for $q_i(z) = p(z|x_i, \\theta)$:\n",
" $$ \\mathcal{L}(X|\\theta) \\geq \\sum_{i=1}^n \\mathbb{E}_{z \\sim p(z|x_i, \\theta)}\\log\\big[\\frac{p(x_i, z|\\theta)}{p(z|x_i, \\theta)}\\big] $$\n",
" * Notice that this equals: $ \\mathcal{F}(\\theta, \\{q_i(\\cdot)\\}_{i=1}^n) = \\sum_{i=1}^n \\mathbb{E}_{z \\sim p(z|x_i, \\theta)}\\log(p(x_i|\\theta)) = \\sum_{i=1}^n \\log p(x_i|\\theta) = \\mathcal{L}(X|\\theta)$ (the first transition is due to the expectation on $z$, not on $x$), which means the **lower bound becomes equality**. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
" * Let's continue:\n",
" $$ \\mathcal{L}(X|\\theta) \\geq \\sum_{i=1}^n \\mathbb{E}_{z \\sim q_i}\\log\\big[\\frac{p(x_i, z|\\theta)}{q_i(z)}\\big] =\\sum_{i=1}^n \\mathbb{E}_{z \\sim q_i}\\big[ \\log p(x_i, z | \\theta) - \\log q_i(z) \\big] $$ $$ = \\sum_{i=1}^n\\big[ \\mathbb{E}_{z \\sim q_i}\\log p(x_i, z | \\theta) - \\mathbb{E}_{z \\sim q_i}\\log q_i(z) \\big] = \\sum_{i=1}^n\\big[ \\mathbb{E}_{z \\sim q_i}\\log p(x_i, z | \\theta) +\\mathcal{H}(q_i(z)) \\big] \\geq Q(\\Theta;\\Theta^t)$$\n",
" \n",
" * $\\mathcal{H}(q_i(z))$ - the **entropy** of $q_i(z)$, which **does not depend on $\\theta$**. Thus, taking the $argmax_{\\theta}$: $$ argmax_{\\theta} \\sum_{i=1}^n\\big[ \\mathbb{E}_{z \\sim q_i}\\log p(x_i, z | \\theta) +\\mathcal{H}(q_i(z)) \\big] = argmax_{\\theta} Q(\\Theta;\\Theta^t) $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"![](\"./assets/tut_10_em_step.jpg\")
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### ![](\"https://img.icons8.com/dusk/64/000000/e.png\")
E Step - EM Receipe\n",
"---\n",
"Estimate $Z$ given $X$ and current $\\Theta$ (the **posterior**):\n",
"\n",
"$$Pr(z_i=j|x_i, \\Theta)= r_{ij} = \\frac{Pr(x_i,z_i=j|\\Theta)}{Pr(x_i|\\Theta)} = \\frac{Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'}Pr(x_i, z_i=j'|\\Theta)} = \\frac{Pr(x_i|z_i=j, \\Theta) \\cdot Pr(z_i=j|\\Theta)}{\\sum_{j'}Pr(x_i|z_i=j', \\Theta) \\cdot Pr(z_i=j'|\\Theta)} $$\n",
"\n",
"* Note, here: $\\Theta = \\Theta^t$ (we use the current estimate of $\\Theta$, and treat it as constant and not as a parameter).\n",
"* Substitute the probabilities with the desired distribution.\n",
"* That means: the probability for sample $x_i$ to be associated with cluster/source $j$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### ![](\"https://img.icons8.com/dusk/64/000000/m.png\")
M Step (Derive Q) - EM Receipe\n",
"---\n",
"$$Q(\\Theta|\\Theta^t) = \\mathbb{E}[\\mathcal{L}(X,Z|\\Theta)|X, \\Theta^t] = \\sum_Z Pr(Z|X,\\Theta^t) \\log Pr(X,Z;\\Theta) = \\sum_i \\sum_{j=1}^k Pr(z_i=j|x_i, \\Theta^t) \\log Pr(x_i, z_i=j, \\Theta) =$$ $$ \\sum_i \\sum_{j=1}^k r_{ij}[\\log Pr(z_i=j|\\Theta) + \\log Pr(x_i|z_i = j, \\Theta)]$$\n",
"\n",
"* $r_{ij} = Pr(z_i=j|x_i, \\Theta)$ (from E step)\n",
"* Substitute $Pr(x_i|z_i=j, \\Theta)$ with the desired probability\n",
"* **Find MLE** (differentiate and compare to 0)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Gaussian Mixture Models (GMMs) as EM\n",
"---\n",
"* One Gaussian: $$ \\mathcal{N}(x| \\mu, \\Sigma) = Pr(x| \\mu_j, \\Sigma_j) = \\frac{1}{(2\\pi)^{\\frac{d}{2}} |\\Sigma_j|^{\\frac{1}{2}}} e^{-\\frac{1}{2} (x-\\mu_j)^T \\Sigma_j^{-1} (x-\\mu_j)} $$\n",
"* Gaussian **Mixture**: $$ Pr(x) = \\sum_{j=1}^k \\alpha_j \\mathcal{N} (x|\\mu_j, \\Sigma_j) $$\n",
" * $\\sum_{j=1}^k \\alpha_j = 1$\n",
" * The parameters of the model are: $\\alpha_j, \\mu_j, \\Sigma_j, \\forall j \\in \\{1,...,k\\}$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* The **log-likelihood** of a GMM: $$ \\mathcal{L}(X|\\Theta) = \\log \\prod_i Pr(x_i|\\Theta) = \\sum_i \\log \\sum_{j=1}^k \\alpha_j \\mathcal{N} (x|\\mu_j, \\Sigma_j) $$\n",
"* **No closed form solution and not convex!**\n",
"* We introduce a **latent random variable** $z$\n",
" * $z \\in \\{0, 1\\}^k$ - a one-hot random variable indicating the source Gaussian the sample belongs to\n",
" * $Pr(z_k) = \\alpha_k$ - the probability of that source\n",
" * Reminder: $\\sum_{j=1}^k \\alpha_k = 1$\n",
"* The marginal probability: $$ Pr(x) = \\sum_z p(z)p(x|z) = \\sum_{j=1}^k \\alpha_j \\mathcal{N} (x|\\mu_j, \\Sigma_j) $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"#### GMM - E-Step \n",
"* The E-step computes the **posterior** probability of the missing data $$ Pr(z_i = j| x_i, \\Theta) = \\frac{Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'|\n",
"\\Theta)} = \\frac{\\alpha_j Pr(x_i|\\mu_j, \\Sigma_j)}{\\sum_{j'} \\alpha_{j'}Pr(x_i|\\mu_{j'}, \\Sigma_{j'})} = \\frac{\\alpha_j |\\Sigma_j|^{-\\frac{1}{2}} e^{-\\frac{1}{2}(x_i-\\mu_j)^T\\Sigma_j^{-1}(x_i-\\mu_j)}}{\\sum_{j'} \\alpha_{j'} |\\Sigma_j'|^{-\\frac{1}{2}} e^{-\\frac{1}{2}(x_i-\\mu_{j'})^T\\Sigma_{j'}^{-1}(x_i-\\mu_{j'})}} $$\n",
"* Denote: $r_{ij} = Pr(z_i=j|x_i, \\Theta)$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"#### GMM - Calculate $ Q(\\Theta;\\Theta^t)$ \n",
"$$Q(\\Theta|\\Theta^t) = \\mathbb{E}[\\mathcal{L}(X,Z|\\Theta)|X, \\Theta^t] = \\sum_Z Pr(Z|X,\\Theta^t) \\log Pr(X,Z;\\Theta) = \\sum_i \\sum_{j=1}^k Pr(z_i=j|x_i, \\Theta^t) \\log Pr(x_i, z_i=j, \\Theta) =$$ $$ \\sum_i \\sum_{j=1}^k r_{ij}[\\log Pr(z_i=j|\\Theta) + \\log Pr(x_i|z_i = j, \\Theta)] = $$ $$ \\sum_i \\sum_{j=1}^k r_{ij}[ \\log \\alpha_j + \\log Pr(x_i|\\mu_j, \\Sigma_j)] = $$ $$ \\sum_i \\sum_{j=1}^k r_{ij} \\log \\alpha_j -\\frac{1}{2} \\sum_{j=1}^k \\log |\\Sigma_j| \\sum_i r_{i,j} -\\frac{1}{2} \\sum_i \\sum_{j=1}^k r_{ij}(x_i-\\mu_j)^T \\Sigma_j^{-1} (x_i-\\mu_j) + Const $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"#### GMM - M-Step \n",
"* To maximize $Q(\\Theta;\\Theta^t)$ with respect to $\\mu_j$, we set the gradient to zero.\n",
"* Reminder: $ \\frac{\\partial}{ \\partial s}(x-As)^TW(x-As) = -2A^TW(x-As)$\n",
"* Derive: $$ \\frac{\\partial}{\\partial \\mu_j} Q(\\Theta; \\Theta^t) = \\sum_{i=1}^n r_{ij} \\Sigma_j^{-1} (x_i-\\mu_j) = 0 \\rightarrow \\hat{\\mu}_j = \\frac{\\sum_{i=1}^n r_{ij}x_i}{\\sum_{i=1}^n r_{ij}}$$\n",
"* Similarly: $$ \\hat{\\Sigma}_j = \\frac{\\sum_{i=1}^n r_{ij}(x_i-\\mu_j)(x_i-\\mu_j)^T}{\\sum_{i=1}^n r_{ij}} $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* To maximize $Q(\\Theta;\\Theta^t)$ with respect to $\\alpha_j$:\n",
" * Use **Lagrange** multiplier:\n",
" * $\\max Q(\\Theta;\\Theta^t) $ s.t. $\\sum_j \\alpha_j = 1 \\iff$\n",
" * $\\mathcal{L} = Q(\\Theta;\\Theta^t) +\\lambda(1 - \\sum_j \\alpha_j) $\n",
" * $\\frac{\\partial \\mathcal{L}}{\\partial \\alpha_j} = \\sum_i \\frac{r_{ij}}{\\alpha_j} - \\lambda = 0$\n",
" * Find an expression for $\\lambda$ by **summing all partial derivatives** of $\\alpha_j$: $$ \\sum_i r_{ij}^{(t)} = \\lambda \\alpha_j \\rightarrow \\sum_j \\sum_i r_{ij}^{(t)} = \\lambda \\sum_j \\alpha_j \\rightarrow \\lambda = n $$\n",
" * $\\sum_j \\sum_i r_{ij}^{(t)} = \\sum_j \\sum_i Pr(z_i=j|x_i, \\Theta^t)= \\sum_i \\frac{\\sum_j Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'| \\Theta)} = \\sum_{i=1}^n 1 = n$\n",
" * Substituting $\\lambda$ back in the *Lagrangian* derivative: $$ \\frac{\\partial \\mathcal{L}}{\\partial \\alpha_j} = \\sum_i \\frac{r_{ij}}{\\alpha_j} - \\lambda = 0 \\rightarrow \\hat{\\alpha}_j = \\frac{\\sum_{i=1}^n r_{ij}}{n} $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* To sum up: $$ \\hat{\\mu}_j = \\frac{\\sum_{i=1}^n r_{ij}x_i}{\\sum_{i=1}^n r_{ij}} $$ $$ \\hat{\\Sigma}_j = \\frac{\\sum_{i=1}^n r_{ij}(x_i-\\mu_j)(x_i-\\mu_j)^T}{\\sum_{i=1}^n r_{ij}} $$
$$ \\hat{\\alpha}_j = \\frac{\\sum_{i=1}^n r_{ij}}{n} $$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "###
![](\"https://img.icons8.com/cotton/64/000000/coin-in-hand.png\")
Bernoulli Mixture Models (BMMs) as EM\n",
"---\n",
"* We have $k$ coins such that:\n",
" * The probability of observing *heads* with the $j^{th}$ coin is $p_j$.\n",
" * We do not observe which coin was used.\n",
" * We only observe $x_i \\in \\{0, 1\\} $, which records whether we see a *heads* or *tails*.\n",
"* Let $z_i \\in \\{1,...,k \\}$ be the **missing information** of which coin was used on each flip (in other words, the *source* like in the GMM case).\n",
" * The probability of using the $j^{th}$ coin is $Pr(z_i=j) = \\alpha_j$ (which is a *parameter*)\n",
"* The complete data is given by $(X,Z)$\n",
" * Using the **law of total probability**, the (marginal) probability of the observed data $X$: $$ Pr(X) = \\sum_j Pr(X|Z=j)Pr(Z=j) $$\n",
" * Thus, the *likelihood* of the full data set (incomplete likelihood) is: $$ \\mathcal{L}(X|\\Theta)= \\prod_i \\sum_j Pr(x_i|z_i=j)Pr(z_i=j) = \\prod_i \\sum_j \\alpha_j p_j^{x_i} (1-p_j)^{1-x_i} $$\n",
" * $\\Theta = (\\alpha, p)$\n",
"\n",
"#### BMM - E-Step \n",
"* The E-step computes the **posterior** probability of the missing data $$ Pr(z_i = j| x_i, \\Theta) = \\frac{Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'|\n",
"\\Theta)} = \\frac{\\alpha_j Pr(x_i|p_j)}{\\sum_{j'} \\alpha_{j'}Pr(x_i|p_{j'})} = \\frac{\\alpha_j p_j^{x_i}(1-p_j)^{1-x_i}}{\\sum_{j'} \\alpha_{j'}p_{j'}^{x_i} (1-p_{j'})^{1-x_i}} $$\n",
"* Denote: $r_{ij} = Pr(z_i=j|x_i, \\Theta)$\n",
"\n",
"#### BMM - Calculate $ Q(\\Theta;\\Theta^t)$ \n",
"$$Q(\\Theta|\\Theta^t) = \\mathbb{E}[\\mathcal{L}(X,Z|\\Theta)|X, \\Theta^t] = \\sum_Z Pr(Z|X,\\Theta^t) \\log Pr(X,Z;\\Theta) = \\sum_i \\sum_{j=1}^k Pr(z_i=j|x_i, \\Theta^t) \\log Pr(x_i, z_i=j, \\Theta) =$$ $$ \\sum_i \\sum_{j=1}^k r_{ij}[\\log Pr(z_i=j|\\Theta) + \\log Pr(x_i|z_i = j, \\Theta)] = $$ $$ \\sum_i \\sum_{j=1}^k r_{ij}[ \\log \\alpha_j + \\log Pr(x_i|p_j)] = $$ $$ \\sum_i \\sum_{j=1}^k r_{ij} \\log \\alpha_j +\\sum_i \\sum_{j=1}^k r_{ij} \\log \\big( p_j^{x_i} (1-p_j)^{1-x_i} \\big) = $$ $$ \\sum_i \\sum_{j=1}^k r_{ij} \\log \\alpha_j + \\sum_i \\sum_{j=1}^k r_{ij}x_i\\log(p_j) +r_{ij}(1-x_i)\\log(1-p_j)$$\n",
"\n",
"#### BMM - M-Step \n",
"* To maximize $Q(\\Theta;\\Theta^t)$ with respect to $p_j$, we set the gradient to zero.\n",
"* Derive: $$ \\frac{\\partial}{\\partial p_j} Q(\\Theta; \\Theta^t) = \\sum_{i=1}^n r_{ij} \\big( \\frac{x_i}{p_j} - \\frac{1-x_i}{1-p_j} \\big) = 0 \\rightarrow \\hat{p}_j = \\frac{\\sum_{i=1}^n r_{ij}x_i}{\\sum_{i=1}^n r_{ij}}$$\n",
"* To maximize $Q(\\Theta;\\Theta^t)$ with respect to $\\alpha_j$:\n",
" * Use **Lagrange** multiplier:\n",
" * $\\max Q(\\Theta;\\Theta^t) $ s.t. $\\sum_j \\alpha_j = 1 \\iff$\n",
" * $\\mathcal{L} = Q(\\Theta;\\Theta^t) +\\lambda(1 - \\sum_j \\alpha_j) $\n",
" * $\\frac{\\partial \\mathcal{L}}{\\partial \\alpha_j} = \\sum_i \\frac{r_{ij}}{\\alpha_j} - \\lambda = 0$\n",
" * Find an expression for $\\lambda$ by **summing all partial derivatives** of $\\alpha_j$: $$ \\sum_i r_{ij}^{(t)} = \\lambda \\alpha_j \\rightarrow \\sum_j \\sum_i r_{ij}^{(t)} = \\lambda \\sum_j \\alpha_j \\rightarrow \\lambda = n $$\n",
" * $\\sum_j \\sum_i r_{ij}^{(t)} = \\sum_j \\sum_i Pr(z_i=j|x_i, \\Theta^t)= \\sum_i \\frac{\\sum_j Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'| \\Theta)} = \\sum_{i=1}^n 1 = n$\n",
" * Substituting $\\lambda$ back in the *Lagrangian* derivative: $$ \\frac{\\partial \\mathcal{L}}{\\partial \\alpha_j} = \\sum_i \\frac{r_{ij}}{\\alpha_j} - \\lambda = 0 \\rightarrow \\hat{\\alpha}_j = \\frac{\\sum_{i=1}^n r_{ij}}{n} $$\n",
"* To sum up: $$ \\hat{p}_j = \\frac{\\sum_{i=1}^n r_{ij}x_i}{\\sum_{i=1}^n r_{ij}} $$
$$ \\hat{\\alpha}_j = \\frac{\\sum_{i=1}^n r_{ij}}{n} $$\n", "\n", "* If all $r_{ij} = \\{0,1\\}$, that is, deterministic, then:\n", " * The component labels, $z_i$, are **known**\n", " * The above update equations reduce to the standard formulas for *binomial distribution*." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##
![](\"https://img.icons8.com/office/80/000000/crystal.png\")
Relation to K-Means: K-Means as \"Hard GMM\"\n",
"---\n",
"* Recall the E-Step for GMM: $$ Pr(z_i = j| x_i, \\Theta) = \\frac{Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'|\n",
"\\Theta)} = \\frac{\\alpha_j Pr(x_i|\\mu_j, \\Sigma_j)}{\\sum_{j'} \\alpha_{j'}Pr(x_i|\\mu_{j'}, \\Sigma_{j'})} = \\frac{\\alpha_j e^{-\\frac{1}{2}(x_i-\\mu_j)^T\\Sigma_j^{-1}(x_i-\\mu_j)}}{\\sum_{j'} \\alpha_{j'} e^{-\\frac{1}{2}(x_i-\\mu_{j'})^T\\Sigma_{j'}^{-1}(x_i-\\mu_{j'})}} $$\n",
"* Let's assume all the Gaussians have the same $\\Sigma = \\epsilon I$ : $$ Pr(z_i = j| x_i, \\Theta) = \\frac{Pr(x_i,z_i=j|\\Theta)}{\\sum_{j'} Pr(x_i,z_i=j'|\n",
"\\Theta)} = \\frac{\\alpha_j e^{-\\frac{1}{2 \\epsilon}||(x_i-\\mu_j)||_2^2)}}{\\sum_{j'} \\alpha_{j'} e^{-\\frac{1}{2 \\epsilon}||(x_i-\\mu_{j'})||_2^2}} $$\n",
"* At the limit $\\epsilon \\to 0$: $Pr(z_i=j|x_i, \\Theta) = 1$ for $j = argmin\\{ x_i - \\mu_j\\}$ and $Pr(z_i=j|x_i, \\Theta) = 0$ for all others."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"* Thus: $$ \\begin{cases}\n",
" r_{ij} = 1 & \\quad \\text{if } j = argmin_j||x_i - \\mu_j||_2^2 \\\\\n",
" r_{ij} = 0 & \\quad \\text{else }\n",
" \\end{cases} $$\n",
"* The GMM equations are now identical to the K-Means' eqautions: $$ \\hat{\\mu}_j = \\frac{\\sum_{i=1}^n r_{ij}x_i}{\\sum_{i=1}^n r_{ij}} $$ $$ \\hat{\\alpha}_j = \\frac{\\sum_{i=1}^n r_{ij}}{n} $$\n", " * The $\\alpha$s are not really required." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "###
![](\"https://img.icons8.com/bubbles/50/000000/video-playlist.png\")
Recommended Videos\n",
"---\n",
"#### ![](\"https://img.icons8.com/cute-clipart/64/000000/warning-shield.png\")
Warning!\n",
"* These videos do not replace the lectures and tutorials.\n",
"* Please use these to get a better understanding of the material, and not as an alternative to the written material.\n",
"\n",
"#### Video By Subject\n",
"\n",
"* EM Motivation:\n",
" * Part 1 (7:53 min)\n",
" * Part 2 (10:39 min)\n",
" * Part 3 (3:05 min)\n",
" * Part 4 (3:29 min)\n",
" * Part 5 (10:53 min)\n",
" \n",
"* EM Lecture - K-Means + GMMs - **VERY RECOMMENDED** - Stanford - Machine Learning (12) - Andrew Ng\n",
" * Andrew Ng is a *l-e-g-e-n-d-a-r-y* ML professor."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"## ![](\"https://img.icons8.com/dusk/64/000000/prize.png\")
Credits\n",
"---\n",
"* Based on **Pattern Recognition and Machine Learning** by *Christopher Bishop* (chapter 9) and slides by *Shai Fine*\n",
"* Icons from Icon8.com - https://icons8.com\n",
"* Datasets from Kaggle - https://www.kaggle.com/"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}