notebook_8_implementacion_gmm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementación EM para GMMs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inicialización:\n",
    "Decidir cuantos clusters (c) utilizaremos. Luego, inicializar los parámetros para $c$ normales: $c$ medias ($\\mu_c$), $c$  desvios ($\\sigma_c$), y $c$ pesos ($\\pi_c$).\n",
    "\n",
    "### Paso E:\n",
    "\n",
    "Calcular para cada punto $x^{(i)}$ la probabilidad ($r_{ic}$) de que el punto pertenezca al cluster c:\n",
    "$$r_{ic} = P(c | x^{(i)}) = \\frac{P(c)P(x^{(i)} | c)}{P(x^{(i)})} = \\frac{\\pi_c N(x^{(i)} \\ | \\ \\mu_c,\\sigma_c)}{\\Sigma_{k=1}^K \\pi_k N(x^{(i)} \\ | \\ \\mu_k,\\sigma_k)}$$\n",
    "\n",
    "donde \n",
    "\n",
    "$${\\displaystyle {\\begin{aligned}N(x^{(i)},\\mu_c,\\sigma_c)&{}={\\frac {1}{\\sigma_c {\\sqrt {2\\pi }}}}e^{-{\\frac {(x^{(i)}-\\mu_c )^{2}}{2\\sigma_c ^{2}}}},\\quad \\pi = 3.141592... .\\\\\\end{aligned}}}$$\n",
    "\n",
    "### Paso M:\n",
    "Para cada cluster c, actualizar $\\pi_c$, $\\mu_c$, y $\\sigma_c$ según los datos:\n",
    "\n",
    "$$ N_c = \\sum_{i}{r_{ic}}$$\n",
    "\n",
    "$$ \\mu^{(nuevo)}_c = \\frac{1}{N_c}\\sum_{i}{r_{ic}x^{(i)}}$$\n",
    "$$ \\sigma^{(nuevo)}_c = \\frac{1}{N_c}\\sum_{i}{r_{ic}(x^{(i)} - \\mu^{(nuevo)}_c)^2} $$\n",
    "$$ \\pi^{(nuevo)}_c = \\frac{N_c}{n}$$\n",
    "\n",
    "\n",
    "\n",
    "### Iterar \n",
    "Iterar hasta que la log-likelihood del modelo converja:\n",
    "\n",
    "$$ln \\ p(\\boldsymbol{X} \\ | \\ \\boldsymbol{\\pi},\\boldsymbol{\\mu},\\boldsymbol{\\sigma}) \\ = \\ \\Sigma_{i=1}^N \\ ln(\\Sigma_{k=1}^K \\pi_k N(x^{(i)} \\ | \\ \\mu_k,\\sigma_k))$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.stats\n",
    "from sklearn.cluster import KMeans\n",
    "import math \n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "def gmm_train(X, n_clusters, tolerance = 1e-3):\n",
    "    v_mu, v_sigma, v_pi = inicializacion(n_clusters, X)\n",
    "    old_llk = 0\n",
    "    tolerance = 1e-3 # indica la mínima ganancia en llk para que el algoritmo continúe.\n",
    "    gain = np.inf\n",
    "\n",
    "    while(gain > tolerance):\n",
    "        print_iter(v_mu, v_sigma, v_pi, n_clusters)\n",
    "        r = e(X, n_clusters, v_mu, v_sigma, v_pi)\n",
    "        v_mu, v_sigma, v_pi = m(X, n_clusters, r)\n",
    "        llk = log_likelihood(X, v_mu, v_sigma, v_pi, n_clusters)\n",
    "        print(\"Log likelihood: \", llk)\n",
    "        gain = llk - old_llk\n",
    "   \n",
    "    return v_mu, v_sigma, v_pi\n",
    "\n",
    "\n",
    "def inicializacion(n_clusters, X):\n",
    "    ##################\n",
    "    # COMPLETAR\n",
    "    # Inicializar v_mu, v_sigma, v_pi (vectores que contienen una componente por cluster)\n",
    "    ##################\n",
    "    \n",
    "    return v_mu, v_sigma, v_pi\n",
    "\n",
    "def norm(x, mu, sigma):\n",
    "    return scipy.stats.norm(mu, sigma).pdf(x)\n",
    "\n",
    "def e(X, n_clusters, v_mu, v_sigma, v_pi):\n",
    "    r = np.ones((len(X), n_clusters))\n",
    "    \n",
    "    ##################\n",
    "    # COMPLETAR\n",
    "    # r[i, c] debera contener la probabilidad estimada P(c | x^(i))\n",
    "    ##################\n",
    "    \n",
    "    return r\n",
    "\n",
    "def m(X, n_clusters, r):\n",
    "    N = np.zeros(n_clusters)\n",
    "    v_mu_new = np.zeros(n_clusters)\n",
    "    v_sigma_new = np.zeros(n_clusters)\n",
    "    v_pi_new = np.zeros(n_clusters)\n",
    "    \n",
    "    ##################\n",
    "    # COMPLETAR\n",
    "    # v_mu_new, v_sigma_new, v_pi_new deberán contener los valores actualizados.\n",
    "    ##################\n",
    "    \n",
    "    return v_mu_new, v_sigma_new, v_pi_new\n",
    "\n",
    "\n",
    "def log_likelihood(X, v_mu, v_sigma, v_pi, n_clusters):\n",
    "    suma = 0\n",
    "    for i in range(len(X)):\n",
    "        suma_k = 0\n",
    "        for k in range(n_clusters):\n",
    "            suma_k += v_pi[k] * norm(X[i], v_mu[k], v_sigma[k])\n",
    "        suma += math.log(suma_k)\n",
    "    return suma\n",
    "\n",
    "\n",
    "def print_iter( v_mu, v_sigma, v_pi, n_clusters):\n",
    "    for c in range(n_clusters):\n",
    "        print(\"C{} => mu={} +- {} (pi={})\".format(c+1, \"%.2f\" % v_mu[c], \"%.2f\" % v_sigma[c], \"%.2f\" % v_pi[c]), end=\" | \")\n",
    "    print()\n",
    "\n",
    "def dibujar(X, v_mu, v_sigma, v_pi, n_clusters, title=\"GMMs\"):\n",
    "    plt.figure()\n",
    "    plt.title(title)\n",
    "    colors = sns.color_palette(\"hls\", n_clusters)\n",
    "    plt.plot(X, [1] * len(X), \".\")\n",
    "    for k in range(n_clusters):\n",
    "        col = colors[k]\n",
    "#         plt.plot([v_mu[k], v_mu[k]], [0, 2], label=\"$\\mu_{}$\".format(k), color=col)\n",
    "        plt.plot([v_mu[k] - v_sigma[k], v_mu[k] + v_sigma[k]], [1.1,1.1], color=col, label=\"$\\mu_{}$\".format(k))\n",
    "        plt.plot([v_mu[k] - v_sigma[k], v_mu[k] - v_sigma[k]], [0.9,1.2], color=col)\n",
    "        plt.plot([v_mu[k] + v_sigma[k], v_mu[k] + v_sigma[k]], [0.9,1.2], color=col)\n",
    "    plt.yticks([])\n",
    "    plt.ylim([0, 2])\n",
    "    plt.legend()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DESCOMENTAR lo anterior\n",
      "Salida esperada:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAEICAYAAABoLY4BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADpxJREFUeJzt3X9s3PV9x/HX24njq8CM4ISsi0ku\n2QZrRJyUpJUnMcliGm2gUZmU1dpSKUIpK4Ioq1QqlvJHIEUKnbZVQqqGtjFqte4PVFColGgKILRN\nmtIphtjrCgMT2eZYiM3VgZBhY3Pv/XFf4Obmkvh89/2+8T0fkpXz3ff7/Xy/3/vmma+/tr8xdxcA\nIHstWa8AAKCMIANAEAQZAIIgyAAQBEEGgCAIMgAEQZARhpndb2Y/mO9rwGJBkAEgCIKMpmZmS7Ne\nB+ADBBmZMLN7zex1MztrZv9tZn845/VWM/uRmT1hZsvOM3+3mf27mZ0xs0Ez66l47XYzezFZ9kkz\n+2rFaz1mVkjGf0PSYxXPfd3Mxs3slJnd3sjtB86HICN1ZnadpD2SPuPu7ZI+J2mk4vVPSDokaVrS\nl9z9vTnzr5Z0WNKDkq6SdI+kJ8xsZTLJuKQvSLpC0u2SvmNmN1Qs4jeT+dZK+vOK535D0mpJuyV9\n18yW12mTgUtCkJGF9yW1SdpgZq3uPuLuryavXSHpnyW9Kul2d3//PPN/WdIRdz/i7iV3f1rScUm3\nSJK7H3b3V73sXyQdlfQHFfOXJO1392l3fzd5bkbSAXefcfcjkt6RdF19Nxu4MIKM1Ln7sKSvSbpf\n0riZ/djMfit5uVtSl6SHvPqdr9ZK+pPkcsUZMzsj6UZJn5QkM9tmZsfM7FfJa7dIWlEx/4S7T81Z\nZtHdZys+/19Jly9gM4F5I8jIhLv/0N1vVDmuLunbyUtHJR2U9KyZraoy+2uSvu/uV1Z8XObuD5lZ\nm6QnJP21pFXufqWkI5KscvhGbBOwUAQZqTOz68zspiSeU5LeVfkygiTJ3f9K0g9VjvKK8yziB5K2\nm9nnzGyJmeWSb8x1Slqm8uWQCUmzZrZN0s2N3iagHggystAm6SFJb0p6Q9LVkvZVTuDu31L5G3vP\nmNlVc157TdIXJX1T5fC+Jukbklrc/aykvZIelzQp6c8k/ayRGwPUi3GDegCIgTNkAAiCIANAEAQZ\nAIIgyAAQxLxurLJixQrP5/MNWhUAWJwGBgbedPeVF5tuXkHO5/M6fvx47WsFAE3IzEYvZTouWQBA\nEAQZAIIgyAAQBP9bAoAwZmZmVCgUNDU192Z8Hw+5XE6dnZ1qbW2taX6CDCCMQqGg9vZ25fN5mdnF\nZwjE3VUsFlUoFLRu3bqalsElCwBhTE1NqaOj42MXY0kyM3V0dCzo7J4gAwjl4xjjDyx03QkyAARB\nkAEgCIIMAEEQZACYo6enRy+99JIkqVgs6vrrr09lXH7sDUBIE/19mh4bqesy29bktXLnrotONzw8\nrGuvvVaSNDQ0pI0bN9Z1PaohyABQYXR0VKtXr1ZLS/kCwtDQkLq6unTu3DndddddWrZsmXp6erRz\n5866j02QAYR0KWeyjTA4OKiurq4PPx8YGFBvb6+efPJJ7dixQ9u3b1dvb29Dgsw1ZACocOLEiQ9/\nueOVV17RU089pY0bN6pQKOiaa66RJC1ZsqQhY3OGDAAVBgcHlcvltGnTJnV1dWnDhg3q6+tTPp9X\noVDQ5s2bVSqVGjI2QQaACkNDQ3r++efV3t7+/54/d+6c9uzZo8OHD2v79u0NGZsgA0Di7NmzMrNf\ni7EkXXbZZXrssccaOj7XkAEg0d7erpdffjmz8QkyAARBkAEgCIIMAEEQZAAIgiADQBAEGQCCIMgA\nEARBBoAgCDIAzMEN6gGgQt9Ev0amx+q6zHzbGu1aefHbZmZ1g3rOkAGgQrUb1J88eVK7d+/Wjh07\nGjY2Z8gAQrqUM9lGqHaD+vXr1+vRRx9taJA5QwaACtVuUJ8GggwAFQYHB1UqlbRp0yYdOHDgwxvU\np4FLFgBQodoN6ovFou677z698MILOnjwoPbt21f3sQkyACQudIP6jo4OPfLIIw0dn0sWAJDgBvUA\nAEkEGQDCIMgAEARBBoAgCDIABEGQASAIggwAQRBkAAiCIAPAHNygHgAq9PdNaGxkuq7LXJNv085d\nKy86HTeoB4AAqt2g/tChQ7rjjjvU29uro0ePNmRszpABhHQpZ7KNUO0G9bfeeqtuu+02TU5O6p57\n7tHNN99c97E5QwaAChe7Qf2DDz6ou+++uyFjE2QAqFDtBvXurnvvvVfbtm3TDTfc0JCxuWQBABWq\n3aD+4Ycf1jPPPKO33npLw8PDuvPOO+s+NkEGgMSFblC/d+9e7d27t6Hjc8kCABLcoB4AIIkgA0AY\nBBlAKO6e9SrUbKHrTpABhJHL5VQsFj+WUXZ3FYtF5XK5mpfBT1kACKOzs1OFQkETExNZr0pNcrmc\nOjs7a56fIAMIo7W1VevWrct6NTLDJQsACIIgA0AQBBkAgiDIABAEQQaAIAgyAARBkAEgCIIMAEEQ\nZAAIgiADQBAEGQCCIMgAEARBBoAgCDIABEGQASAIggwAQRBkAAiCIANAEAQZAIIgyAAQBEEGgCAI\nMgAEQZABIAiCDABBEGQACCKTIPf3Tai/byLMcpC9if4+TfT3Zb4MLEzfRL/6JvqzXo26SnOblqYy\nyhxjI9OhloPsTY+NhFgGFmZkeizrVai7NLeJSxYAEARBBoAgCDIABEGQASAIggwAQRBkAAiCIANA\nEAQZAIIgyAAQRCa/qTd+ekZTUyUdfKCwoOWMjU4rl+PflMVgZvy0SlNTKhx8oOZlTI+NqiWXq+Na\nYb5Oz4xrqjSlBwoHs16VuhmdHlOuJZ3jipoBQBCZnCFfvapVkrRvf+eClrPQM2zE0Xr1KklS5779\nNS9jIWfXqI9VrVdLkvZ37st4TeonzbN9zpABIAiCDABBEGQACIIgA0AQBBkAgiDIABAEQQaAIFIJ\n8sDopL773LAGRifDzYP0pfneckyk653pWf3PmXdTeZ/SmqfWbapFw38xZGB0Ujv/8Zjemy1p2dIW\n9X+lu6Z5tqxdXvd5kL4031uOiXQNjE7qxVNvq+TSzp8ca+j7lOY8tWxTrRp+hnzsZFHvzZZUcmlm\ntqRjJ4th5kH60nxvOSbSdexkUSUvP270+5TmPLVsU60aHuTu9R1atrRFS0xqXdqi7vUdYeZB+tJ8\nbzkm0tW9vkMtVn7c6PcpzXlq2aZaNfySxZa1y9X/lW4dO1lU9/oObVm7XEd1bt7z1DIO4knzveWY\nSNeWtcv1qdkr9Pa7M/rWPL60T+uYqHWeWrapVqncXGjL2uXz3pC05kH60nxvOSbSdXnbUl3etlRb\nOmP+fa9lnlq3qRb82BsABEGQASAIggwAQRBkAAiCIANAEAQZAIIgyAAQRCb/yemafFuo5SB7bWvy\nIZaBhcm3rcl6FeouzW0yd7/kibdu3erHjx9v4OoAwOJjZgPuvvVi03HJAgCCIMgAEARBBoAgCDIA\nBEGQASAIggwAQRBkAAiCIANAEAQZAIIgyAAQBEEGgCAIMgAEQZABIAiCDABBEGQACIIgA0AQBBkA\ngiDIABAEQQaAIAgyAARBkAEgCIIMAEEQZAAIgiADQBAEGQCCIMgAEARBBoAgCDIABEGQASAIggwA\nQRBkAAiCIANAEAQZAIIgyAAQBEEGgCAIMgAEQZABIAiCDABBEGQACIIgA0AQBBkAgiDIABAEQQaA\nIAgyAARBkAEgCIIMAEEQZAAIgiADQBAEGQCCIMgAEARBBoAgCDIABEGQASAIggwAQRBkAAiCIANA\nEAQZAIIgyAAQBEEGgCAIMgAEQZABIAiCDABBEGQACIIgA0AQBBkAgiDIABAEQQaAIAgyAARBkAEg\nCIIMAEEQZAAIgiADQBAEGQCCIMgAEARBBoAgCDIABEGQASAIggwAQRBkAAiCIANAEAQZAIIgyAAQ\nBEEGgCAIMgAEQZABIAiCDABBEGQACIIgA0AQBBkAgiDIABAEQQaAIAgyAARBkAEgCIIMAEEQZAAI\ngiADQBDm7pc+sdmEpNEax1oh6c0a511M2A9l7Icy9sNHFvO+WOvuKy820byCvBBmdtzdt6YyWGDs\nhzL2Qxn74SPsCy5ZAEAYBBkAgkgzyH+f4liRsR/K2A9l7IePNP2+SO0aMgDgwrhkAQBBEGQACKIh\nQTazfzKzcTP7RcVzV5nZ02b2SvLn8kaMHUmV/XC/mb1uZieSj1uyXMc0mNk1Zvacmf3SzP7LzP4i\neb6pjokL7IemOibMLGdm/2Fmg8l+eCB5fp2Z/dzMhs3sJ2a2LOt1TVujzpC/J+nzc577S0nPuvvv\nSno2+Xyx+55+fT9I0nfcfXPycSTldcrCrKSvu/sGSd2S7jazDWq+Y6LafpCa65iYlnSTu2+StFnS\n582sW9K3Vd4PvyNpUtLuDNcxEw0Jsrv/q6RfzXn6i5L6ksd9km5rxNiRVNkPTcfdT7n788njs5Je\nlLRaTXZMXGA/NBUveyf5tDX5cEk3Sfpp8vyiPx7OJ81ryKvc/VTy+A1Jq1IcO5o9ZjaUXNJY1F+m\nz2VmeUmflvRzNfExMWc/SE12TJjZEjM7IWlc0tOSXpV0xt1nk0kKasJ/rDL5pp6Xf9auWX/e7u8k\n/bbKX6qdkvQ32a5OeszscklPSPqau79d+VozHRPn2Q9Nd0y4+/vuvllSp6TPSvq9jFcphDSDfNrM\nPilJyZ/jKY4dhrufTg7GkqR/UPlgXPTMrFXlCPW7+5PJ0013TJxvPzTrMSFJ7n5G0nOSfl/SlWa2\nNHmpU9Lrma1YRtIM8s8k7Uoe75L0VIpjh/FBgBJ/LOkX1aZdLMzMJD0q6UV3/9uKl5rqmKi2H5rt\nmDCzlWZ2ZfL4E5L+SOXr6c9J2pFMtuiPh/NpyG/qmdmPJPWofDu905L2Szok6XFJa1S+heeX3H1R\nf8Oryn7oUflLU5c0IumrFddRFyUzu1HSv0n6T0ml5Olvqnz9tGmOiQvshz9VEx0TZtal8jftlqh8\nUvi4ux8ws/WSfizpKkkvSPqyu09nt6bp41enASAIflMPAIIgyAAQBEEGgCAIMgAEQZABIAiCDABB\nEGQACOL/AGjYcPhNzSyUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Datos para ejemplo\n",
    "np.random.seed(134)\n",
    "X = np.array([10, 11, 12, 20, 21, 22, 30, 31, 32, 29, 33])\n",
    "n_clusters = 3\n",
    "\n",
    "# DESCOMENTAR\n",
    "# v_mu, v_sigma, v_pi = gmm_train(X, n_clusters)\n",
    "# dibujar(X, v_mu, v_sigma, v_pi, n_clusters, title=\"implementacion propia\")\n",
    "print(\"DESCOMENTAR lo anterior\")\n",
    "\n",
    "print(\"Salida esperada:\")\n",
    "# SALIDA ESPERADA:\n",
    "sklearn_gmm = GaussianMixture(n_components=3)\n",
    "sklearn_gmm.fit(X.reshape(-1, 1))\n",
    "\n",
    "dibujar(X, sklearn_gmm.means_.ravel(), sklearn_gmm.covariances_.ravel(), sklearn_gmm.weights_.ravel(), n_clusters, title=\"sklearn\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DESCOMENTAR lo anterior\n",
      "Salida esperada:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAEICAYAAABoLY4BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAE8RJREFUeJzt3XtwVOd5x/Hfo+sCFgYECIwQQmlC\ngwET7DR46iZM2hCThLHT0JCETD2u48S1PTQzoeOQdOpLnZJkculkJhO3qe1oHMhl6lsmdjLEiSdu\nG9sN2AhD4wtWhCyCkVhjwLIlJO3bP/YsPVpr0W1X5xF8PzM73j3n7HmfffXuz2dfrV4shCAAQPLK\nki4AAJBFIAOAEwQyADhBIAOAEwQyADhBIAOAEwQy3DCzW8zs+6PdB5wtCGQAcIJAxjnNzCqSrgHI\nIZCRCDO7ycwOmdlJM3vOzP48b3+lmf3AzO41s6ohnr/azH5jZq+aWYuZrYntu9rMfhedu9XMPhPb\nt8bMOqL2X5Z0d2zb58ys08wOm9nVpXz9wFAIZEw4M1si6UZJ7wwh1Eh6v6S22P4pkh6Q1CvpoyGE\nU3nPXyDpIUm3S5olaYuke81sTnRIp6QPSZou6WpJ3zSzVbFTzIuet0jSp2Pbzpe0QNI1kr5tZjOL\n9JKBESGQkYQBSdWSlppZZQihLYTwYrRvuqSfS3pR0tUhhIEhnv9JSQ+HEB4OIWRCCL+QtEvSByQp\nhPBQCOHFkPVrSTsl/Vns+RlJN4cQekMIb0Tb+iTdFkLoCyE8LOk1SUuK+7KBMyOQMeFCCAckfVbS\nLZI6zeyHZnZBtHu1pBWSvhwKr3y1SNJfRdMVr5rZq5IukzRfksxsnZk9YWavRPs+IGl27PldIYSe\nvHOmQwj9scevSzpvHC8TGDUCGYkIIewIIVymbLgGSV+Jdu2UtE3SL82srsDTX5J0TwhhRuw2LYTw\nZTOrlnSvpK9JqgshzJD0sCSLN1+K1wSMF4GMCWdmS8zsvVF49kh6Q9lpBElSCOGrknYoG8qzhzjF\n9yWtN7P3m1m5maWiX8zVS6pSdjqkS1K/ma2TtLbUrwkoBgIZSaiW9GVJRyW9LGmupK3xA0II/6Ts\nL/YeMbNZefteknSFpC8oG7wvSfp7SWUhhJOSNkv6saRjkj4h6SelfDFAsRgL1AOAD1whA4ATBDIA\nOEEgA4ATBDIAODGqhVVmz54dGhsbS1QKAJyddu/efTSEMGe440YVyI2Njdq1a9fYqwKAc5CZHRzJ\ncUxZAIATBDIAOEEgA4AT/GsJANzo6+tTR0eHenryF+ObHFKplOrr61VZWTmm5xPIANzo6OhQTU2N\nGhsbZWbDP8GREILS6bQ6Ojq0ePHiMZ2DKQsAbvT09Ki2tnbShbEkmZlqa2vHdXVPIANwZTKGcc54\nayeQAcAJAhkAnCCQAcAJAhkA8qxZs0bPPvusJCmdTmvZsmUT0i5fewPgUtf2ZvW2txX1nNUNjZqz\n6aphjztw4IDe9ra3SZL27t2r5cuXF7WOQghkAIg5ePCgFixYoLKy7ATC3r17tWLFCnV3d+v6669X\nVVWV1qxZo02bNhW9bQIZgEsjuZIthZaWFq1YseL04927d2vjxo267777tGHDBq1fv14bN24sSSAz\nhwwAMXv27Dn9xx0vvPCCHnzwQS1fvlwdHR1auHChJKm8vLwkbXOFDAAxLS0tSqVSuuiii7RixQot\nXbpUzc3NamxsVEdHh1auXKlMJlOStglkAIjZu3evnnrqKdXU1Aza3t3drRtvvFEPPfSQ1q9fX5K2\nCWQAiJw8eVJm9qYwlqRp06bp7rvvLmn7zCEDQKSmpkbPP/98Yu0TyADgBIEMAE4QyADgBIEMAE4Q\nyADgBIEMAE4QyADgBIEMAE4QyACQhwXqASCmuWu72nrbi3rOxuoGXTVn+GUzk1qgnitkAIgptEB9\na2urrrnmGm3YsKFkbXOFDMClkVzJlkKhBeqbmpp05513ljSQuUIGgJhCC9RPBAIZAGJaWlqUyWR0\n0UUX6bbbbju9QP1EYMoCAGIKLVCfTqf1xS9+UU8//bS2bdumrVu3Fr1tAhkAImdaoL62tlZ33HFH\nSdtnygIAIixQDwCQRCADgBsEMgA4QSADgBMEMgA4QSADgBMEMgA4QSADgBMEMgDkYYF6AIjZ3tyl\n9rbeop6zobFam66aM+xxLFAPAA4UWqD+gQce0LXXXquNGzdq586dJWmbK2QALo3kSrYUCi1Q/8EP\nflBXXnmljh07pi1btmjt2rVFb5srZACIGW6B+ttvv1033HBDSdomkAEgptAC9SEE3XTTTVq3bp1W\nrVpVkraZsgCAmEIL1H/rW9/SI488ouPHj+vAgQO67rrrit42gQwAkTMtUL9582Zt3ry5pO0zZQEA\nERaoBwBIIpABwA0CGYArIYSkSxiz8dZOIANwI5VKKZ1OT8pQDiEonU4rlUqN+Rx8ywKAG/X19ero\n6FBXV1fSpYxJKpVSfX39mJ9PIANwo7KyUosXL066jMQwZQEAThDIAOAEgQwAThDIAOAEgQwAThDI\nAOAEgQwAThDIAOAEgQwAThDIAOAEgQwAThDIAOAEgQwAThDIAOAEgQwAThDIAOAEgQwAThDIAOAE\ngQwAThDIAOAEgQwAThDIAOAEgQwAThDIAOAEgQwATrgJ5O3NXdre3JV0GaMyGWuO69rerK7tze7O\nNR7NXdvV3LU96TJwFpnIMVUxIa2MQHtbb9IljNpkrDmut73N5bnGo623PekScJaZyDHl5goZAM51\nBDIAOEEgA4ATBDIAOEEgA4ATBDIAOEEgA4ATBDIAOEEgA4ATbv5Sr/NIn3p6Mtp2a0fSpYxY+8Fe\npVKT9/9pfZ1HlOnpUce2W8d9rt72gypLpYpQ1fgc6etUT6ZHt3ZsS7oUnCUO9rYrVTYxY3vypgkA\nnGXcXCHPrauUJG29uT7hSkZuMl3ND6Vybp0kqX7rzeM+VzGusouhrnKuJOnm+q0JV4KzxUR+2uIK\nGQCcIJABwAkCGQCcIJABwAkCGQCcIJABwAkCGQCcmJDvIe8+eExPtKa1uqlWFy+a+aZ99z3Voa6j\nVZpTU3362JlTq7T/D8cVJH1kVb0uXjRTO55s113//Xu90TegBeenNGNqlY69fkptR7vV0zeg2vOq\nVVFepqbZ07RmyVzt+8NxHThyUr39GV3aVKsTvf0ySTXVFXq8Na2XXnldr/X2a+n86XrfhfM0Mzrf\n6qZaSdK//vpFPdV+TL39A3r7vOm6ad3bT2/ff/iEGjrqNO/85P867UwK9f2OJ9s169BxlZnpO/c/\no+7efv3y2U6d6s9odk213jJ7mp4/clINs6bqynfUa/8fjqvzZK/m1lTrwgvO16PPder3R7s1a2ql\n/qiuRn/d2y9J+vajBzRzapUefa5TnSd6tPGdDfrEuxpGVNNksePJdv3ot+2qm57SZ97zltOvIf91\nxR9L0h2/flGdJ3p0aVOtaqZUanVTrZ57+aTu+q9WnejpU2V5mRbMmKIZU6s0u6ZaH1mV/U5+7v2Q\n69Pc8184clJ7XnpVDbOm6vgbfaqbntKaJXN1z+Ntaj3arZpUhTasqlfr0W61dr2m13r61dXdqxCk\ncpPqpqf0oRUXqGZKpU6+0afHW9OaOz2laVXleuyFLmUyQX0DQf2ZoHnTU6qbXq3WdLdO9WW0YMYU\n1c+aqkPHXtcr3dn3zFvranTyjT7tP3xC65bN1yfe1XD6/R0kLbvgfN3zeJva0t2aW5PSZW+drb+M\nvcah+q1Q3w4nd3x+Pfn7vY1BCyGM+OBLLrkk7Nq1a1QN7D54TJv+/Qmd6s+oqqJM2z+1elAnf/y7\n2X0rDy+UmbS//pD6BzLKxMqqKjf9zZ8u1h2PtY6q7bEwSZUVZRrIZDSQGbyvzCQz00BU3MrDCyVJ\nH71+xptCx4NCfb/jyXZ94f5n9PX0TknS52rXjrutb7ySPdeW2rWDfnaS9M8fXn66f840Hooh9yX+\nUv1hSK7vcirKTT/69KWSNOh1/eOHLtRtP92vU/0ZVZSXaWAgo4FYv1j03L6Bwu+/inJTmaS+gaCR\nv0v9uO7dTbrrN2061Z8peExFuanMTP0Db+633PiQNKoxkxtjPX2D282Nw9GOwWKMKTPbHUK4ZLjj\nSj5l8URrWqf6swHb15/RE63pQfv6Yj+sEB2T/4buGwj6+f6XS11qtgZla8gPY0nKBJ0O47if7Ttc\n+sLGoFDfl6LeELK3IbpnUHtnGg+TQX7f9Q8EPdGaftPr+tm+w4Me5+dukM4YxrlzT9YwlqSf7395\n0Pt7KP0D4fR7fqh+G6pvhxszuePz5X52nsdgyQN5dVOtqirKVG7ZK8/cx7fcvsqK/y/BomPKbPA5\nKstNl184r9SlSsp2SGVFmcqH6Jkyk8rzi5O0btn80hc2BoX6vhT1mmVvQ3TPoPbONB4mg/y+qyg3\nrW6qfdPrWrds/qDH5Xn9YsqO6zOpKDdVRlfJk9HlF84b9P4eSkW5ne6fofptqL4dbszkjs/v3dzP\nzvMYLPkc8sWLZmr7p1YPOV9z8aKZ+sG1q7NzyL/KziH/w7ULC84hN9ROczWHPCVdrnnnp1xOV0iF\n+z5X77R7fqkyM216V8O455CXtkyXJH3uPUvOOId8pvEwGeRey1BzyPmva8m8mnN+Dvl9F84b9Rxy\nvN8K9e2ZxMfYUHPInsdgyeeQRyq3UM9kXFxoMtUcl1sQqJiLCxXjXONR6jlknHvOqjlkAMDIEMgA\n4ASBDABOEMgA4ASBDABOEMgA4ASBDABOuPlHThsaq5MuYdQmY81x1Q2NLs81Ho3VPv9IB5PXRI4p\nN38YAgBnK/4wBAAmGQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHA\nCQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZ\nAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwg\nkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHA\nCQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZ\nAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwg\nkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHA\nCQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJwgkAHACQIZAJywEMLIDzbrknRwmMNmSzo6\nnqJKxGtdErWNldfavNYlUdtYjbe2RSGEOcMdNKpAHgkz2xVCuKSoJy0Cr3VJ1DZWXmvzWpdEbWM1\nUbUxZQEAThDIAOBEKQL530pwzmLwWpdEbWPltTavdUnUNlYTUlvR55ABAGPDlAUAOEEgA4ATRQ9k\nM7vFzA6Z2Z7o9oFitzGGmi43s+fM7ICZfT7peuLMrM3Mnon6alfCtdxlZp1mti+2bZaZ/cLMXoj+\nO9NJXS7GmZktNLNHzex/zWy/mf1dtN1DvxWqLdG+M7OUmf2PmbVEdd0abV9sZk9G79MfmVnVRNY1\nTG3fM7Pfx/psZUkKCCEU9SbpFklbin3ecdRTLulFSU2SqiS1SFqadF2x+tokzU66jqiWd0taJWlf\nbNtXJX0+uv95SV9xUpeLcSZpvqRV0f0aSc9LWuqk3wrVlmjfSTJJ50X3KyU9KWm1pB9L+li0/Q5J\nf+uotu9J2lDq9s+FKYs/kXQghNAaQjgl6YeSrki4JpdCCI9JeiVv8xWSmqP7zZKunNCiVLAuF0II\nh0MIT0X3T0r6naQF8tFvhWpLVMh6LXpYGd2CpPdK+o9oe1J9Vqi2CVGqQL7RzPZGHzUn/KNangWS\nXoo97pCDQRkTJO00s91m9umkixlCXQjhcHT/ZUl1SRaTx9M4k5k1SnqHsldVrvotrzYp4b4zs3Iz\n2yOpU9IvlP0U+2oIoT86JLH3aX5tIYRcn30p6rNvmll1KdoeUyCb2SNmtm+I2xWSviPpLZJWSjos\n6etFrPdsdFkIYZWkdZJuMLN3J11QISH7Oc7L9yRdjTMzO0/SvZI+G0I4Ed+XdL8NUVvifRdCGAgh\nrJRUr+yn2D+e6BoKya/NzJZJ2qpsje+UNEvSTaVou2IsTwoh/MVIjjOz70r66VjaKKJDkhbGHtdH\n21wIIRyK/ttpZvcrOzgfS7aqQY6Y2fwQwmEzm6/sVUPiQghHcveTHmdmVqls4G0PIdwXbXbRb0PV\n5qnvQgivmtmjki6VNMPMKqKr5MTfp7HaLg8hfC3a3Gtmd0vaUoo2S/Eti/mxhx+WtK/QsRPkt5Le\nGv0Gt0rSxyT9JOGaJElmNs3ManL3Ja1V8v2V7yeSroruXyXpwQRrOc3LODMzk3SnpN+FEL4R25V4\nvxWqLem+M7M5ZjYjuj9F0vuUnd9+VNKG6LCk+myo2p7N9VnUp1eqRH1WitXe7lH2o1BQ9hsEn4nN\npSUi+lrPvyj7jYu7QghfSrKeHDNrknR/9LBC0o4kazOzH0hao+xSg0ck3SzpAWV/+92g7NKrHw0h\nTOgv2ArUtUYOxpmZXSbpPyU9IykTbf6CsnO1Sfdbodo+rgT7zsxWKPtLu3JlLwp/HEK4LXo//FDZ\nKYGnJX0yhNA7UXUNU9uvJM1R9lsYeyRdF/vlX/HaL3YgAwDG5lz42hsATAoEMgA4QSADgBMEMgA4\nQSADgBMEMgA4QSADgBP/BwdMfkI7y7qHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Datos para ejemplo\n",
    "np.random.seed(134)\n",
    "n_clusters = 3\n",
    "X = np.concatenate([np.random.randn(50) * 2, (np.random.randn(30) + 10), (np.random.randn(100) * 3 + 25)])\n",
    "# v_mu, v_sigma, v_pi = gmm_train(X, n_clusters)\n",
    "# dibujar(X, v_mu, v_sigma, v_pi, n_clusters, title=\"implementacion propia\")\n",
    "\n",
    "print(\"DESCOMENTAR lo anterior\")\n",
    "print(\"Salida esperada:\")\n",
    "\n",
    "sklearn_gmm = GaussianMixture(n_components=n_clusters)\n",
    "sklearn_gmm.fit(X.reshape(-1, 1))\n",
    "\n",
    "dibujar(X, sklearn_gmm.means_.ravel(), sklearn_gmm.covariances_.ravel(), sklearn_gmm.weights_.ravel(), n_clusters, title=\"sklearn\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}