20191022 27일차 a (K-Means with Nunpy).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K-Means를 numpy로 구현하기\n",
    "- Numpy로 속도를 더 빠르게 할 수 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "K = 3 # cluster 개수\n",
    "N = 100 # data 개수\n",
    "V = 2 # 색인어(Index/Unique/Distinct terms)\n",
    "\n",
    "np.random.seed(0)\n",
    "X = np.random.randint(99, size=(N, V)) +1 # dataset\n",
    "C = np.random.randint(99, size=(K, V)) +1 # initial centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVEElEQVR4nO3dfYxldXnA8e+z4NtKyfKyGoRdFlOiNSSKThC1MQSctCAR/pBmzAZRaOYfW/GlEZr9w/QPIqamSqOlmcBWbCauiqQQY2wmFGL6h7Sz1CovWiiyOwsIK2+2bgISnv5xz9TZ5c7OfTnn3vPy/SSTO/fcO3N/52Weec5znnNOZCaSpHbZNO0BSJLKZ3CXpBYyuEtSCxncJamFDO6S1EIGd0lqoWM3ekNE7AYuBp7KzLOKaScC3wJ2AI8Cf5KZz0ZEADcAFwGHgI9l5r0bfcbJJ5+cO3bsGHEWJKmb9u7d+6vM3NrvtQ2DO/B14KvAN9ZMuxa4MzOvj4hri+fXABcCZxZf7wZuLB6PaseOHSwvLw8wFEnSqojYt95rG5ZlMvOHwDNHTL4EuKX4/hbg0jXTv5E9PwK2RMQpww9ZkjSOUWvub8zMJwCKxzcU008FVta870AxTZI0QWUfUI0+0/pe3yAi5iNiOSKWDx48WPIwJKnbRg3uT66WW4rHp4rpB4Bta953GvB4v1+QmQuZOZOZM1u39j0eIEka0ajB/Q7giuL7K4Db10z/aPScCzy/Wr6RJE3OIK2Q3wTOA06OiAPA54HrgW9HxFXAfuCy4u3fp9cG+TC9VsiPVzBmSdIGNgzumfmRdV66oM97E/jEuIOSJI3HM1QnaGVlhd27d7OysrLxmyVpDAb3CVpaWmLPnj0sLS1NeyiSWm6QM1RVktnZ2cMeh7GyssLS0hKzs7Ns27Zt4x+Q1Glm7hO0bds2rrzyypGCs1m/pGF0InNvQ9Y7TtZfJ21YF1ITdCJzb0PWO07WXydtWBdSE3Qic29L1tsGrgtpMqLXmj5dMzMz6SV/JWk4EbE3M2f6vdaJsowkdY3BXZJayOAuSS1kcJfUOF7KY2MGd0mNY0vtxjrRCimpXWyp3ZjBXVLjrJ7Up/VZlpGkFjK4S1ILGdzVSXZbDK6Oy6qOY6obg7s6yW6LwdVxWdVxTHXjAVVNzSCX/63qEsF2WwyujsuqjmOqGy8cpqnZvXs3e/bsYW5ubt3Oh0HeI3XV0S4cZuY+orUZJeANKEYwSPZlhiaNxuA+otWa36rV780uBzdIr7L9zNJoDO4j6pdR9ssuva2cms5tuJkM7iM6MqNcL7tcm+GbgaqJ3IabyeBeMWvGk2F2WR234WYyuFfMmvFkmF1Wx224mQzuaoUuZZfupWgQnqGqVljNLrsQ7Dw7U4Mwc5capqy9FPcA2s3MXWqYsvZS3ANoNzN3qaO6dJyii8bK3CPi0xFxf0TcFxHfjIjXRsQZEXFPRDwUEd+KiFeXNVhJ5enScYouGjm4R8SpwCeBmcw8CzgGmAO+CHw5M88EngWuKmOg+p2hr2W9uAg7dsCmTb3HxcUqhyepBsatuR8LvC4ijgU2A08A5wO3Fq/fAlw65mfoCEPVShcXYX4e9u2DzN7j/LwBXmq5kYN7Zj4GfAnYTy+oPw/sBZ7LzJeKtx0ATh13kDrc7Owsc3NzG9ZKV1ZW+N+rr4ZDhw5/4dAh2LWrwhGq6bzTUfONU5Y5AbgEOAN4E/B64MI+b+17wfiImI+I5YhYPnjw4KjD6KRBa6VLS0tsfvrp/i/u31/ByNQWdtI03zjdMh8AfpGZBwEi4jbgvcCWiDi2yN5PAx7v98OZuQAsQO9mHWOMQ+uYnZ3l0EkncVy/AL99++QHpMawk6b5xqm57wfOjYjNERHABcADwF3Ah4v3XAHcPt4QNapt27Zx3A03wObNh7+weTNcd910BqVGsJOm+capud9D78DpvcBPi9+1AFwDfCYiHgZOAm4uYZwa1c6dsLAAp58OEXD66Tz9hS+w+4UXrKdKLTZWt0xmfj4z35qZZ2Xm5Zn5QmY+kpnnZObvZ+ZlmflCWYPViHbuhEcfhZdfhkcf5fbjjrOeKrWcZ6h2kPVUqf0M7h3k9bml9vPCYZoo+6cFbgeTYHDXRNk/LXA7mATLMpoo6/0Ct4NJiMzpnz80MzOTy8vL0x6GJDVKROzNzJl+r1mWkaQWMrhLUgsZ3CWphQzuahRb6KTBGNzVKLbQSYOxFVKNYgudNBiDuxrFSydIg7EsI0ktZHCXpBYyuEtSyerQ1WVwl6SS1aGrywOqhZWVFZaWlpidnfW+kVJJuvp3VYeuLjP3Qh3+00pt09W/qzrcYLxzmft6mUQd/tNKZalLxuzf1fR0LnNfL5Oow39aqSx1yZj9u5qezmXuZhLqArdzebMOSWoob9ahiahDb2+TufxUJoO7SlOXOm9TufxUps7V3FWdQeq8deniWKsuY7JOrjKZuas0g3RG1DE7rcuY7CxRmczcNVF1zE7rOCZpXHbLSFJD2S0jSTVUZYeUwV2SpqTK4z3W3KWOqkuXUJdVebzHzF3qqLp0CXVZlR1SYwX3iNgSEbdGxM8i4sGIeE9EnBgRSxHxUPF4QlmDldqgLmeizs7OMjc3V+suobosqyYaN3O/AfhBZr4VeDvwIHAtcGdmngncWTyXVKhLxtyEvvq6LKsmGrnmHhHHA+8HPgaQmS8CL0bEJcB5xdtuAe4GrhlnkFKb2Fc/OJfV6Ebuc4+IdwALwAP0sva9wNXAY5m5Zc37ns3MV5RmImIemAfYvn37u/bt2zfSOCSpq6rqcz8WeCdwY2aeDfyGIUowmbmQmTOZObN169YxhtE91iF71i4Hl4l0uHFaIQ8ABzLznuL5rfSC+5MRcUpmPhERpwBPjTtIHW61Dglw5ZVXTnk007N2OQAuE2mNkYN7Zv4yIlYi4i2Z+XPgAnolmgeAK4Dri8fbSxmp/p91yJ5+y6Hry0RaNda1ZYq6+03Aq4FHgI/TK/V8G9gO7Acuy8xnjvZ7mnJtGU/6kFQnR6u5j3WGamb+GOj3iy8Y5/fWleUQSU3h5QeGYDlEXefea3N4+YEhNOGkD6lKnlTUHAZ3NcfiIuzYAZs29R4XF6c9os5pwiUL1GNZRs2wuAjz83DoUO/5vn295wA7d05vXB2zuveq+jNzVzPs2vW7wL7q0KHedEmvYHBXM+zfP9x0qeMM7mqG7duHmy51nMFdzXDddbB58+HTNm/uTZf0CgZ3NcPOnbCwAKefDhG9x4UFD6YOaZoXWPPibpNlcFdz7NwJjz4KL7/cezSwD22afer2yE9Wa1shPZNOeqVpnmXtGd6T1drg7nVgpFeaZp+6PfKT1drgbpYgtYt748Npbc3d68BI7WLNfjitzdylUZgd1pd748NpbeYujcLssL7cGx+Ombsq16Rs2OxQbWHmrso1KRs2O1RbmLmrctPMhpu01yCVycxdlZtmNtykvQapTGbuJTA7rC9r6OoqM/cSmB3WlzV0dZWZewnMDiXVjcG9BF4zQ1LdWJaRpBYyuEtSCxncR+RdZSTVmcF9RHbISKozD6iOyA6ZZvPcBNVBlduhmfuI7J9uNve8VAdVbodm7uqkLu15jZIdumczGVVuh2bu6qQu7XmNkh26ZzMZVW6HY2fuEXEMsAw8lpkXR8QZwB7gROBe4PLMfHHcz5E0mlGywy7t2bRVZOZ4vyDiM8AMcHwR3L8N3JaZeyLi74H/zMwbj/Y7ZmZmcnl5eaxxSFLXRMTezJzp99pYZZmIOA34IHBT8TyA84Fbi7fcAlw6zmdIkoY3bs39K8DngJeL5ycBz2XmS8XzA8CpY35Ga3jik6RJGTm4R8TFwFOZuXft5D5v7Vv3iYj5iFiOiOWDBw+OOoxG8SCVpEkZ54Dq+4APRcRFwGuB4+ll8lsi4tgiez8NeLzfD2fmArAAvZr7GONoDA9SSZqUsQ+oAkTEecBfFAdUvwN8d80B1Z9k5t8d7ec9oCpJw6vsgOo6rgE+ExEP06vB31zBZ0jqGI9ZDaeUM1Qz827g7uL7R4Bzyvi9krRq9ZgV4M1xBtDoM1T9T65Jaur21tRxH2l2dpa5uTmPWQ2o0cHd7hNNUlO3t6aO+0hdumREGRp94TC7TzRJVWxvk7hAl38n3VRKt8y47JZRV+3evZs9e/YwNzdnHVlDO1q3TKMzd6npzKpVFYO7NEWrdWSpbI0+oCpJ6s/gLkktZHBXJdrSW63pcjsancFdlWhLb7Wmy+1odB5QVWnW9mzbBaIyuB2NzuCu0hx57Q+7QDQuu4lGZ3BXacyyNKhJnJnbddbcVRqv/aFBWUuvnpl7R5gpDc9lVh338qpn5t4RZkrDc5lVx7286pm5d4SZ0vDKXGbuBWjSDO4dYdfB8MpcZt5FSJNmcJcmwD0nTZrBXZoA95w0aR5QlTQSr/tSbwZ3SSOxm6jeLMu0gJ0YmgaPI9SbmXsLmEFpGuxVrzcz9xYwg5J0JIN7C9iJIelIlmUkqYUM7lJNtKW1sC3z0XQGd6km2nJgvC3z0XTW3KWaaMuB8bbMR9NFZk57DMzMzOTy8vK0hyFJjRIRezNzpt9rlmUkTVyZdXlr/P2NHNwjYltE3BURD0bE/RFxdTH9xIhYioiHiscTyhuupDYosy5vjb+/cWruLwGfzcx7I+L3gL0RsQR8DLgzM6+PiGuBa4Frxh+qVB4v2TBdZdblrfH3V1rNPSJuB75afJ2XmU9ExCnA3Zn5lqP9rDV3Tdru3bvZs2cPc3NzngCmxqq85h4RO4CzgXuAN2bmEwDF4xvK+Iyu61JdcaB5XVyEHTtg06be4+LiUJ8xOzvL3Nyc2Z5aa+zgHhHHAd8FPpWZvx7i5+YjYjkilg8ePDjuMFqvS3XFDed1cRHm52HfPsjsPc7PDxXgveiV2m6sPveIeBW9wL6YmbcVk5+MiFPWlGWe6vezmbkALECvLDPOOLqgS3XFDed11y44dOjwaYcO9abv3Fnx6KRmGLnmHhEB3AI8k5mfWjP9r4Gn1xxQPTEzP3e032XNXUPZtKmXsR8pAl5+efLjkaakqpr7+4DLgfMj4sfF10XA9cBsRDwEzBbPpYFtWHPfvn246VIHjVyWycx/BWKdly8Y9fdKqzV3oH8ny3XX9Wrsa0szmzf3pksCvLaMamjDmvtqXX3XLti/v5exX3ddY+rt9thrEgzuqp2Bbj6yc2djgvmRNtwzkUpgcJcmrEudT4Nyb6Z8XjhMmjB77F+pS+dxTEqjM3f/25fL5alpcW+mfI3O3P1vXy6Xp6bFvZnyNTpz9799uUZZnmb70uiq/PtpdObuf/tyjbI8zfal0VX599PozF3T597T4dyTaZ5prrMq/34anblr+tx7Opx7Ms0zzXVW5d+Pmbtap62ZmKrR1nVmcFfrTPMM0IHOrlWttHWdGdzVOm3NxKRhGNzVOm3NxKRheEBV6qhx78vbpfv6NpHBXeqocbtE7AyqN8syUkeNe2zCYxv1NvI9VMvkPVQlaXhV3UNVklRTBndJaiGDuyS1kMFdklrI4C5JLWRwl0rkiT2qC4O7VCJP7FFdeBKTVCJP7BmfNzwph5m7VCJvXjI+937KYeauzjNTrJcm7f3Uedsxc1fnmSnWS5P2fuq87Zi5d0SdM4xpq2um6Dqrv7puO2Dm3hl1zjCmra6Zouus/uq67YCZe2dUlWGYXVanzlmh6q+SzD0i/jgifh4RD0fEtVV8hoZTVYZhdlmdOmeFqr/Sg3tEHAN8DbgQeBvwkYh4W9mfo3qYnZ1lbm7uqNllV8/a7NJ8d2le11O3ZVBF5n4O8HBmPpKZLwJ7gEsq+BzVwCDZZVez+y7Nd5fmdT11WwZV1NxPBdb+6zoAvLuCz1FDdLV23KX57tK8rqduy6D02+xFxGXAH2XmnxbPLwfOycw/P+J988A8wPbt29+1b9++UschSW036dvsHQDW7qOfBjx+5JsycyEzZzJzZuvWrRUMQ5K6q4rg/u/AmRFxRkS8GpgD7qjgcyRJ6yi95p6ZL0XEnwH/DBwD7M7M+8v+HEnS+io5iSkzvw98v4rfLUnamJcfkKQWMrhLUgsZ3CWphUrvcx9pEBEHgWEa3U8GflXRcOrM+e6ers678z2Y0zOzby95LYL7sCJieb3G/TZzvrunq/PufI/PsowktZDBXZJaqKnBfWHaA5gS57t7ujrvzveYGllzlyQdXVMzd0nSUTQuuHflFn4RsS0i7oqIByPi/oi4uph+YkQsRcRDxeMJ0x5rFSLimIj4j4j4XvH8jIi4p5jvbxUXpWuViNgSEbdGxM+K9f6eLqzviPh0sY3fFxHfjIjXtnF9R8TuiHgqIu5bM63v+o2evy3i3E8i4p3Dfl6jgnvHbuH3EvDZzPwD4FzgE8W8XgvcmZlnAncWz9voauDBNc+/CHy5mO9ngaumMqpq3QD8IDPfCryd3vy3en1HxKnAJ4GZzDyL3sUG52jn+v468MdHTFtv/V4InFl8zQM3DvthjQrudOgWfpn5RGbeW3z/P/T+0E+lN7+3FG+7Bbh0OiOsTkScBnwQuKl4HsD5wK3FW1o33xFxPPB+4GaAzHwxM5+jA+ub3gUMXxcRxwKbgSdo4frOzB8Czxwxeb31ewnwjez5EbAlIk4Z5vOaFtz73cLv1CmNZWIiYgdwNnAP8MbMfAJ6/wCAN0xvZJX5CvA54OXi+UnAc5n5UvG8jev9zcBB4B+KctRNEfF6Wr6+M/Mx4EvAfnpB/XlgL+1f36vWW79jx7qmBffoM63V7T4RcRzwXeBTmfnraY+nahFxMfBUZu5dO7nPW9u23o8F3gncmJlnA7+hZSWYfooa8yXAGcCbgNfTK0kcqW3reyNjb/NNC+4D3cKvLSLiVfQC+2Jm3lZMfnJ196x4fGpa46vI+4APRcSj9Mpu59PL5LcUu+3QzvV+ADiQmfcUz2+lF+zbvr4/APwiMw9m5m+B24D30v71vWq99Tt2rGtacO/MLfyKOvPNwIOZ+TdrXroDuKL4/grg9kmPrUqZ+ZeZeVpm7qC3fv8lM3cCdwEfLt7Wxvn+JbASEW8pJl0APEDL1ze9csy5EbG52OZX57vV63uN9dbvHcBHi66Zc4HnV8s3A8vMRn0BFwH/Bfw3sGva46lwPv+Q3m7YT4AfF18X0as/3wk8VDyeOO2xVrgMzgO+V3z/ZuDfgIeB7wCvmfb4KpjfdwDLxTr/J+CELqxv4K+AnwH3Af8IvKaN6xv4Jr3jCr+ll5lftd76pVeW+VoR535Kr5toqM/zDFVJaqGmlWUkSQMwuEtSCxncJamFDO6S1EIGd0lqIYO7JLWQwV2SWsjgLkkt9H90vUSsjdoI8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], c=\"k\", s=2, alpha=0.5)\n",
    "plt.scatter(C[:,0], C[:,1], c=\"r\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K-means with Euclide distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2309.197988874544\n",
      "2190.869848936423\n",
      "2177.278914010543\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n",
      "2175.1491677408344\n"
     ]
    }
   ],
   "source": [
    "for _ in range(10):\n",
    "    # Expectation\n",
    "    R_nk = np.zeros((N, K))\n",
    "    distance = np.zeros((N, K))\n",
    "    # Ck - X 거리\n",
    "    for k in range(K):\n",
    "        distance[:,k] = np.linalg.norm(X-C[k], axis=1)\n",
    "    # 거리 가장 작은 index, R_nk[idnxex] = 1 \n",
    "    for k in range(K):\n",
    "        R_nk[np.where(\n",
    "            np.argmin(distance, axis=1) == k), k] = 1 \n",
    "\n",
    "    # Maximaization\n",
    "    SSE = 0.0\n",
    "    for k in range(K):\n",
    "        C[k] = X[np.where(R_nk[:,k] == 1)].mean(axis=0)\n",
    "        SSE += np.linalg.norm(X[np.where(R_nk[:,k] ==1)] - C[k], \n",
    "                              axis=1).sum()\n",
    "    print(SSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYXklEQVR4nO3df4xlZX3H8fd32cW6s5EFmaXIsjNoN/4IYpUpwmoIgSbKD4U/itHsKLWY3QSroDZK6x/WP0gwNRUaGWUC6NpMRAVSqBobs9XYZuK2s9iAitYtMrsrKzsWWLu7TRey3/5xzsVh9s6de+/59Zzn+bySzZ179s7c5/y43/s93+c5zzF3R0RE4rKq6QaIiEj5FNxFRCKk4C4iEiEFdxGRCCm4i4hESMFdRCRCq1d6gZndA1wFHHT3c/NlpwFfA8aBJ4B3ufszZmbA7cAVwFHgT9394ZXe4/TTT/fx8fEhV0FEJE27d+/+jbuPdvu/FYM78GXg88BXFi27Gdjp7rea2c35808AlwOb839vBr6QP/Y0Pj7O3NxcH00REZEOM5tf7v9WLMu4+w+Ap5csvhrYkf+8A7hm0fKveOaHwHozO3PwJouISBHD1tzPcPcDAPnjhnz5WcC+Ra/bny8TEZEald2hal2WdZ3fwMy2mdmcmc0tLCyU3AwRkbQNG9yf6pRb8seD+fL9wNmLXrcReLLbH3D3aXefcPeJ0dGu/QEiIjKkYYP7Q8B1+c/XAQ8uWv4+y1wIHOqUb0REpD79DIX8KnAJcLqZ7Qc+BdwKfN3Mrgf2AtfmL/822TDIPWRDId9fQZtFRGQFKwZ3d3/PMv91WZfXOvDBoo2KljscOQIjI2DduidERMqhK1Tr4g533gk33JA9ah59EamQgntdjhyB2VkYG8sejxwZ/G+4w+HD+mIQkRUpuNdlZAS2bIH5+exxZGSw31fmLyIDSCO4h5DxmsH27TA1lT0OWnMvI/MPQQj7QiQB8Qf3kDJeM1i3brjO1KKZfwhC2hcikYs/uMeS8RbN/EMQy74QaYH4g3sMGW9Hkcw/BDHtC5HAmQdwajwxMeGVTvmr8eXh0L4QKY2Z7Xb3iW7/18987u3XyXiledoXIrWIvywjIpIgBXcRaRWNpu2PgruItIZG0/ZPwV1EWkOjafun4C4iraHRtP1LY7SMiEShcy3f5KRG065EwV1EWkWjafujsoykSUMu+hLiZgqxTSFScJf0aMhFX0LcTCG2KVQK7tKcflKwKtI0DbnoS4ibKcQ2hUrBXZrRTwpWVZqmIRd9CXEzhdimUKUxcVhVFk+CBZoQaxCHD2dBe2ws+6ROTZ3YS9bPa4alCcz6EuJmCrFNTdHEYVXoZJWzs3DRRdlRNjubpRNtnW+9Tp0UrLPNuqVg/bxmWBpy0ZcQN1OIbQqRMvdhLc4q9+zJjrhXverEDFNpxvL62Tbafo3S5g9br8xdNfdhLS7+XXxx9m9pIVBd+731c/ORtt+gpMV0+LabyjLDWnqpHJx42dzSrv3JSZ1PVkUpZul0+LabMvciFmeV3TJMde3XQylmJXT4tptq7lVTRlm9KkfVBKbuw0mHb9hUc2+SasbVSyTFbOIERYdve6nmLu2XyFSBZdfAlZXHTZm7xCGBFLPMExR1U8RPmbtIS5R5gqKRMPFT5i7SImWdoCTSTZG0Qpm7mX0E+ADgwKPA+4EzgXuB04CHgfe6+7GC7ZSlVDCVAhLppkja0Jm7mZ0FfBiYcPdzgZOAdwOfAT7n7puBZ4Dry2ioLKKCqZQggW6KpBUty6wGXmpmq4G1wAHgUuC+/P93ANcUfA9Zqt9JrXXLGhmCDps4DB3c3f1XwGeBvWRB/RCwG3jW3Z/PX7YfOKtoI2WJfgqmyu5lCDps4lGkLHMqcDVwDvAKYAS4vMtLux4eZrbNzObMbG5hYWHYZqSpUzCdmlp+emHdskaGoMMmHkXKMn8M/NLdF9z9OeABYAuwPi/TAGwEnuz2y+4+7e4T7j4xOjpaoBmJWqlgquEQMgQdNvEYem4ZM3szcA/wR8D/Al8G5oCLgfvd/V4z+yLwiLtP9fpbUc8t06TlRtRopI30oMOjPSqZW8bdd5F1nD5MNgxyFTANfAL4qJntAV4O3D3se0hB3bJ7FVVlBRpFE4dC49zd/VPAp5Ysfhy4oMjflQrp0kSRJOgK1dSoqCqSBM0tk5qmL01UQVfQYVAHZe4paqqoqnq/oMOgLgruUh8NohZ0GNRFwV3qo3q/oMOgLrqHqtRLxVZBh0FZeo1zV4eq1KtT75ek6TConsoy0h6arlCkbwru0g4aYiEyEAV3aQcNsRAZiIK7tIOGWIgMRME9JjMzMD4Oq1ZljzMzTbeoPP3MYS8iL9BomVjMzMC2bXD0aPZ8fj57DrB1a3PtKpOGWIj0TZl7LD75yd8F9o6jR7PlIpIcBfdY7N072HIRqUQoI3YV3GOxadNgy0WkdCGN2FVw7wjl63ZYt9wCa9e+eNnatdlykQa0/SM1jJBG7Cq4Q1hft8PauhWmp7Ojyix7nJ6OpzNVWiWGj9QwQhqxm97EYd1mLDp8ODsKx8ayvTI1pVEZ0mpNT8yV8keqzm1fyQ2yW2m5dCKkr1uRgkLImlP+SIVyg/G0Mvde6UTTqY5ISULJmvWRqp4y945e6UQoX7ciBYWSNesj1ay0MndQOlElbdtCytx82hVpUOa+mNKJaoRQ6G2xsjefDnNJL7hLNfod4Bvg4OcQmhTS+GiJg4K7lKOfQm+A2X0oTQqlTr6cmCccjVV6NXepzkqF3lCGcQTapFDr5EsnHIXs4mddI9c81dylHisVegNMT0NqUqh1ck042k7K3KVeAaanATYpKKtWdS9XmcHx4/W3JxZlHHfK3CUcAaanATYpKJpwtHx19PUouIskqt9RQppwtHx1jI5ScBdJ0CCZoyYcLV8dfT2Fau5mth64CzgXcODPgJ8DXwPGgSeAd7n7M73+jmrukpIQavwhjRJaTgjbqUqh19xvB77j7q8B3gA8BtwM7HT3zcDO/LmIoHH1/QplO1Wp6r6eoYO7mb0MuBi4G8Ddj7n7s8DVwI78ZTuAa4o2UiQWoVyJagbbt2cZ+/bt4WXGoWynNiuSub8SWAC+ZGY/MrO7zGwEOMPdDwDkjxtKaKdIFELKmEMeJRTSdmqroWvuZjYB/BB4i7vvMrPbgd8CH3L39Yte94y7n9rl97cB2wA2bdp0/vz8/FDtSFbsBcmVLF3/Fm2PFjW1UdpOK6uq5r4f2O/uu/Ln9wFvAp4yszPzNz4TONjtl9192t0n3H1idHS0QDMSlEJBspel63/8eKu2R8gZc0i0nYoZOri7+6+BfWb26nzRZcBPgYeA6/Jl1wEPFmphaDSFYPOWrv/Bg2lvD5EuVhf8/Q8BM2Z2MvA48H6yL4yvm9n1wF7g2oLvEY5Oxjg7mxUCm+qJ6hQkO+1IrSC5dP03bEh7e4h0obllBhHS4ODUC5Itrrm3lTZxeDS3TFlC6sJPvSC5dP1T3x4VS72bp42KlmXS0hkcPDmp9EWSsrSbY3IyvCta5cWUuQ9KGaIkKKSTVumPMncRWZFOWttHwV1E+tI5aZV2UFlGRCRCCu4iIhFScBcRiZCCu0hCmpw9I4SZO1Ki4C6SiCYvRNJFUPWLN7grTRB5kSbnm0t9rrsmxBnclSaInKDJC5F0EVT94pw4LKQJvkQC0uTkX0XeW5OWdZfexGFKE0S6anL2jGHfWyfiw4nzClVdKy1DUoYYHk1aNpw4M3fQBF8yMGWIYdKJ+HDizNwlLC1Jh5Uhhkkn4sOJN3OXMLQoHVaGGC6diA9OmbtUq+F0eJCTBmWIEhNl7lKtBtPhYU4alCFKLJS5l6UldeXaNZgOq4YuKVPmXoYW1ZUb0VA6rBq6pEyZexmUIgZJNXRJmTL3MihFDJZq6JIqZe5lUIooIoFRcC+L7h4sIgFRWaYIzRkvIoFScB+WRsiISMAU3IelW8u0mrtz+NhhQrifgaSp6hN/BfdhaYRMa7k7d+6+kxu+dQN37r5TAV5qV8eJv4L7sDojZKamskeNkGmNI88dYXbfLGOnjDG7b5Yjz8V91jVohqiupOrVceKv4F6EBlG30siaEbacvYX5Q/NsOXsLI2viPesaNENUV1I96jjxL3wPVTM7CZgDfuXuV5nZOcC9wGnAw8B73f1Yr79R+j1URVbg7hx57ggja0awiL+cB72dsG4/XJ8ypqOq+h6qNwKPLXr+GeBz7r4ZeAa4voT3ECmVmbHu5HVRB3YYPENUV1J9qj7xL5S5m9lGYAdwC/BR4B3AAvD77v68mV0E/LW7v63X30kmc9fMkdKAQQ87HabtUWXmfhvwceB4/vzlwLPu/nz+fD9wVsH3iIOKmdKQQTNEdSXFYejgbmZXAQfdfffixV1e2jWKmdk2M5szs7mFhYVhm9EeGhcvIjUqkrm/BXinmT1B1oF6KVkmv97MOnPWbASe7PbL7j7t7hPuPjE6OlqgGS2hYqaI1KjwaBkAM7sE+It8tMw3gPvd/V4z+yLwiLtP9fp91dxFpBd9dLqrerTMUp8APmpme8hq8HdX8B7tpGKmyMDUXTWcUoK7u3/f3a/Kf37c3S9w9z9w92vd/f/KeI8eb67L6aRybZ6Lps1tB3VXDavdV6jqK11q0Oa5aNrc9g51Vw2n3cFdX+lSgyrnoqk6q45hHh1N4zScdgd3faVLDaqai6aOrDqWeXTUXTW4UkbLFFVotIy60aUGVcxFc/jYYW741g2MnTLG/KF5pq6cYt3J5U/kkso8Oimqe7RMvfSVLjWoYi6aurLqVObRkRdrf+Yu0mLKqqWIXpn76m4LRaQenaxapGztL8uIiMgJFNylEm2/cEaap2OoGAV3KV0MF85Is3QMFafgLqXpZFqHjx1u/YUz0qwYLr5qmjpUpRSdTGt23ywXbbyILRu3MLt/ttUXzkhzOsNEZ/fpGBqWhkJKKZZekHPHFXdgZhriJ131MwRUw0RXFvdFTBKEpRfkrDt5nS6cka76rafr4qtiVJZJRNVZkJmx/fztTJ43GU2mpcyxGkvr6ZPnTWqsfwWUuSegrpEHMWVaGq1RnVgmMwudgnsCNPJgcGVvM43Z/p3OWd7UlVNsP397FMlAiBTcE6BMaXBlbjOdBZwoprO8UGm0TCJUPx5cWdusrql9JT0aLSPKlIZQ1jbTmZM0QZm7SA1iO3OKbX3aSpm7SMNiOnNSH0I7KLhHQCMxpE4afdUOCu4tpyxK6qY+hHbQFaotp6v9pG4xXo0cI2XuLacsSpoQUx9CrDRaJgIauRCHGPZjDOvQJhotEzllUe0XQ99JDOsQEwV3kQDEMAIlhnWIiYK7SABi6DuJYR1iopq7SCBiqFf3uw5lrWsM26wI1dxFWiCGvpN+1qGs2rxq/L0NHdzN7Gwz+56ZPWZmPzGzG/Plp5nZd83sF/njqeU1V6Qcuqq3OWXV5lXj761I5v488DF3fy1wIfBBM3sdcDOw0903Azvz5yLBUMbXrLJq86rx91Zazd3MHgQ+n/+7xN0PmNmZwPfd/dW9flc19/6kUl+sej01v3rzVHMvR+U1dzMbB94I7ALOcPcDAPnjhjLeI3WpZJt1rKcyvuaV1b8QQz9FVQoHdzNbB9wP3OTuvx3g97aZ2ZyZzS0sLBRtRvRSqS/WsZ66h6ekoFBwN7M1ZIF9xt0fyBc/lZdjyB8Pdvtdd5929wl3nxgdHS3SjCSkkm3WtZ5VZ3wzj84wfts4qz69ivHbxpl5dKaS9xFZztA1d8s+FTuAp939pkXL/wb4b3e/1cxuBk5z94/3+luqufcnlfri8ePHOXj0IBvWbmDVqvaN1p15dIZt/7iNo88dfWHZ2jVrmX7HNFtfv7XBlklsetXciwT3twL/AjwKHM8X/xVZ3f3rwCZgL3Ctuz/d628puEtHp+Y+u2+WLWdvaWXZZPy2ceYPzZ+wfOyUMZ646Yn6GyTR6hXch57P3d3/FVjuU3fZsH9X0hbD/PR7D+3tuTyVMzBpVvvOeSVqMfQtbDpl07LLUxn1JM1TcJegxDCS5ZbLbmHtmrUvWrZ2zVpuueyWZEY9DUJXC1dDwV2C0/axy1tfv5Xpd0wzdsoYhjF2ytgLnakxnJmUSWcy1Wn9rJCqX5ZH27Ie2s6/o6uFi4l2Vkh965dH27I+bT8zKZPOZKrT6uCu+mV5tC2lCTH0sYSq1cFd3/rlGXZbqjNMikr1TKbqz45q7vKCQbdlDBcciTShrM9OtDV3SPdbvwqDbkuVck6kM5l2aWp/1fHZaX1wl+aoLPZi6pRulyb3Vx2fndaXZaRZIZbFmmqThvW1S9P7q4zjNOqyjDQrtLJY7NmYlKfp/VX1Z0eZu0QlhmxM6tP2/aXMXZIRezYm5Yp5fylzl+i0PRurS5HtpG0cBmXukpSYs7GyFOmb0KigdlBwF0lQkXHWur6hHRTcRRJUpG+i6X4N6Y9q7iKJUs29/Sq5h6qItFunb6Lu35V6qCwjIhIhBXcRkQgpuIuUSLNCSigU3EVKovHfEhIFd5GSaPx3cTrzKY+Cu0hJNP67GJ35lEtDISV5ZY3Z7tzsefK8SY3/HsLSM5/J8yaDHm4Z+lh/Ze6StLKzRc1rM7w2nfm04SxDwT0RqmV2F2qdPMX91TnzmbpyKvibrYd63Cym4J6ANmQZTQkxW0x5f7XlzCfE42YpzS2TgCrvThR63bEfoa1D03eTkv6EcNxoPvfEVZVlxJJhhpYttiErlPCOm6UqydzN7O3A7cBJwF3ufmuv1ytzr14VWUa/GWYIGU4TUpl1sU1trUpT26DWzN3MTgLuAC4HXge8x8xeV/b7yGCqyDL6yTBjye4HVXS9Q88KO1Ldv4uFug2qKMtcAOxx98fd/RhwL3B1Be8jDetndEMbRhVUIZX1TmU9ewl1G1QR3M8C9i16vj9fJhFaKcNMtX6cynqnsp69hLoNSq+5m9m1wNvc/QP58/cCF7j7h5a8bhuwDWDTpk3nz8/Pl9oOCUeqNdlU1juV9ewliZo7WaZ+9qLnG4Enl77I3afdfcLdJ0ZHRytohoSiLfXjsqWy3qmsZy8hboMqgvu/A5vN7BwzOxl4N/BQBe8jIiLLKH3iMHd/3sz+HPgnsqGQ97j7T8p+HxERWV4ls0K6+7eBb1fxt0VEZGW6QlVEJEIK7iIiEVJwFxGJkIK7iEiEgpjy18wWgEGvYjod+E0FzQldiuud4jqD1js1w6z3mLt3vVAoiOA+DDObW+7KrJiluN4prjNovZtuR93KXm+VZUREIqTgLiISoTYH9+mmG9CQFNc7xXUGrXdqSl3v1tbcRURkeW3O3EVEZBmtC+5m9nYz+7mZ7TGzm5tuT1XM7Gwz+56ZPWZmPzGzG/Plp5nZd83sF/njqU23tWxmdpKZ/cjMvpk/P8fMduXr/LV8ttHomNl6M7vPzH6W7/eLYt/fZvaR/Pj+sZl91cx+L8b9bWb3mNlBM/vxomVd961l/i6PcY+Y2ZuGec9WBffE7s/6PPAxd38tcCHwwXxdbwZ2uvtmYGf+PDY3Ao8tev4Z4HP5Oj8DXN9Iq6p3O/Add38N8AaybRDt/jazs4APAxPufi7ZLLLvJs79/WXg7UuWLbdvLwc25/+2AV8Y5g1bFdxJ6P6s7n7A3R/Of/4fsg/6WWTruyN/2Q7gmmZaWA0z2whcCdyVPzfgUuC+/CXRrTOAmb0MuBi4G8Ddj7n7s0S+v8lmpn2pma0G1gIHiHB/u/sPgKeXLF5u314NfMUzPwTWm9mZg75n24J7kvdnNbNx4I3ALuAMdz8A2RcAsKG5llXiNuDjwPH8+cuBZ939+fx5rPv8lcAC8KW8JHWXmY0Q8f52918BnwX2kgX1Q8Bu0tjfsPy+LSXOtS24d7uHVdTDfcxsHXA/cJO7/7bp9lTJzK4CDrr77sWLu7w0xn2+GngT8AV3fyNwhIhKMN3kNeargXOAVwAjZCWJpWLc372Ucsy3Lbj3dX/WWJjZGrLAPuPuD+SLn+qcouWPB5tqXwXeArzTzJ4gK7ldSpbJr89P2yHefb4f2O/uu/Ln95EF+5j39x8Dv3T3BXd/DngA2EIa+xuW37elxLm2Bfdk7s+a15rvBh5z979d9F8PAdflP18HPFh326ri7n/p7hvdfZxs3/6zu28Fvgf8Sf6yqNa5w91/Dewzs1fniy4DfkrE+5usHHOhma3Nj/fOOke/v3PL7duHgPflo2YuBA51yjcDcfdW/QOuAP4T+C/gk023p8L1fCvZqdgjwH/k/64gq0HvBH6RP57WdFsrWv9LgG/mP78S+DdgD/AN4CVNt6+idf5DYC7f5/8AnBr7/gY+DfwM+DHw98BLYtzfwFfJ+hWeI8vMr19u35KVZe7IY9yjZKOJBn5PXaEqIhKhtpVlRESkDwruIiIRUnAXEYmQgruISIQU3EVEIqTgLiISIQV3EZEIKbiLiETo/wGfChBS5qZT/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "color = [\"r\", \"g\", \"b\"]\n",
    "for k in range(K):\n",
    "    plt.scatter(X[np.where(R_nk[:, k]==1)][:,0],\n",
    "               X[np.where(R_nk[:, k]==1)][:,1],\n",
    "               c = color[k], s=5, alpha=0.5)\n",
    "    plt.scatter(C[k,0], C[k,1], c=color[k])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K-Means with cosine similarity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "K = 3\n",
    "N = 100\n",
    "V = 2 # 색인어(Index/Unique/Distinct terms)\n",
    "\n",
    "np.random.seed(0)\n",
    "X = np.random.randint(99, size=(N, V)) +1\n",
    "C = np.random.randint(99, size=(K, V)) +1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "98.58979434236662\n",
      "98.77086349429828\n",
      "98.79764883655051\n",
      "98.84424112091719\n",
      "98.88660370350914\n",
      "98.89389377247672\n",
      "98.89389377247672\n",
      "98.89389377247672\n",
      "98.89389377247672\n",
      "98.89389377247672\n"
     ]
    }
   ],
   "source": [
    "for _ in range(10):\n",
    "    # Expectation\n",
    "    R_nk = np.zeros((N, K))\n",
    "    angle = np.zeros((N, K))\n",
    "    # Ck - X 각도\n",
    "    innerproduct = X.dot(C.T)\n",
    "    Cnorm = np.linalg.norm(C, axis=1)\n",
    "    for k in range(K):\n",
    "        Xnorm = np.linalg.norm(X, axis=1)\n",
    "        angle[:,k] = innerproduct[:,k]/(Xnorm*Cnorm[k])\n",
    "    # 각도 유사도 가장 큰 index, R_nk[idnxex] = 1 \n",
    "    for k in range(K):\n",
    "        R_nk[np.where(\n",
    "            np.argmax(angle, axis=1) == k), k] = 1 \n",
    "\n",
    "    # Maximaization \n",
    "    SSE = 0.0\n",
    "    for k in range(K):\n",
    "        C[k] = X[np.where(R_nk[:,k] == 1)].mean(axis=0)\n",
    "        subX = X[np.where(R_nk[:,k] == 1)]\n",
    "        innerproduct = subX.dot(C[k].T)\n",
    "        Xnorm = np.linalg.norm(subX, axis=1)\n",
    "        Cnorm = np.linalg.norm(C[k])\n",
    "        SSE += (innerproduct/(Xnorm*Cnorm)).sum()\n",
    "    print(SSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAY0ElEQVR4nO3dfYxldX3H8fd3d8G6s4FlZaC4sLNgNz7EYJUpwtIQAiYVoUKaUjUMUsTMhrWKD1Vp+UP7BwmmpkojszLhwbVuBAVSiBobs9WQZpR2Fg0+YMsWnd2FlR0rD53dpuyy3/5xzq3DcGfm3nuefuf3+7ySyZ175s7c83S/8z3f8z2/Y+6OiIjEZUXTMyAiIuVTcBcRiZCCu4hIhBTcRUQipOAuIhIhBXcRkQitWu4FZnYncClwwN3fmE9bB9wDbAR+CfyZuz9jZgbcArwDOAT8ubs/stx7nHjiib5x48YBF0FEJE27du36tbsPd/vZssEd+BLwBeDL86bdAOx095vN7Ib8+SeBi4FN+ddbgW3545I2btzI9PR0D7MiIiIdZjaz2M+WLcu4+0PAbxZMvgzYnn+/Hbh83vQve+YHwFozO6X/WRYRkSIGrbmf7O77AfLHk/Lp64G98163L58mIiI1KvuEqnWZ1nV8AzMbN7NpM5uenZ0teTZERNI2aHB/ulNuyR8P5NP3AafNe92pwFPd/oC7T7r7qLuPDg93PR8gIiIDGjS4PwhcnX9/NfDAvOnvtcw5wHOd8o2IiNSnl1bIrwIXACea2T7gU8DNwNfM7FpgD3BF/vJvkbVB7iZrhbymgnkWEZFlLBvc3f09i/zooi6vdeADRWcqWu5w8CAMDYF1Oz0hIlIOXaFaF3e47TbYujV71Dj6IlIhBfe6HDwIU1MwMpI9HjzY3++7w9yc/imISE8U3OsyNASbN8PMTPY4NNT77yrrF5E+pRHcQ8h6zWDLFpiYyB77qbkXzfpDEsK2EElA/ME9pKzXDNas6f9kapGsPyQhbQuRyMUf3GPIeotk/SGJYVuItET8wT2WrHfQrD8ksWwLkRYwD+DQeHR01Csd8lf95eHQthApjZntcvfRbj/rZTz39utkvdI8bQuRWsRflhGRqKjhqjcK7iLSGmq46p2Cu4i0hhqueqfgLiKtoYar3qVxQlVEotC55GNsTA1Xy1FwF5FWUcNVb1SWERGJkIK7pEn9dD0JcTWFOE8hUnCX9KifrichrqYQ5ylUCu7SnF5SsCrSNPXT9STE1RTiPIVKwV2a0UsKVlWapn66noS4mkKcp1ClMXBYFeYPgAUaDKtfc3NZ0B4ZyT6pExMvb4Ho5TWD0gBmPQlxNYU4T03RwGFl62SUU1Nw7rnZHjY1laUSbR5vvU6dFKyz3rqlYL28ZlDqp+tJiKspxHkKkTL3QczPKHfvzva217yme3apNGNxvawbrb9GafWHbanMXTX3Qcwv/J1/fvbVrQioU/tL6+UGJDHcpKSltPu2m8oyg1h4DTR0vx564an9sTEdT1ZB6WUltPu2mzL3Qc3PKBfLLnVqv3pKLyuj3bfdVHOvmrLKalXZUROguncn7b5hU829SaoZVyuh9LKJgxTtvu2lmru0W0JjwJZdA1dWHjdl7tJ+iaSXZR6k6FRF/JS5i7REmQcp6oSJnzJ3kRYp6yAloVMVySqUuZvZR4D3Aw78GLgGOAW4G1gHPAJc5e4vFJxPmU/FUikooVMVyRo4czez9cCHgFF3fyOwEng38Bngc+6+CXgGuLaMGZWciqVSkkROVSSraFlmFfBKM1sFrAb2AxcC9+Y/3w5cXvA9ZL5eB7TW7WpkQNp14jBwcHf3J4HPAnvIgvpzwC7gWXc/kr9sH7C+6EzKPL0US5Xdy4C068SjSFnmBOAy4HTg1cAQcHGXl3bdPcxs3MymzWx6dnZ20NlIT6dYOjGx+PDCul2NDEi7TjyKlGXeBvzC3Wfd/TBwP7AZWJuXaQBOBZ7q9svuPunuo+4+Ojw8XGA2ErRcsVStEDIg7TrxGHhsGTN7K3An8AfA/wBfAqaB84H73P1uM/si8Ki7Tyz1t6IeW6YpS3XUqNtGlqDdoz0qGVvG3R8mO3H6CFkb5ApgEvgk8FEz2w28Crhj0PeQAhbL7lVUlWWoiyYOhfrc3f1TwKcWTH4COLvI35UK6dJEkSToCtXUqKgqkgSNLZOapi9NVEFX0G5QB2XuKWqqqKp6v6DdoC4K7lIfNVEL2g3qouAu9VG9X9BuUBfdQ1XqpWKroN2gLEv1ueuEqtSrU++XpGk3qJ7KMtIeGq5QpGcK7tIOarEQ6YuCu7SDWixE+qLgLu2gFguRvii4t8GOHbBxI6xYkT3u2NH0HNWvl3HsReT/qVsmdDt2wPg4HDqUPZ+ZyZ4DXHllc/PVBLVYiPRMmXvobrzxt4G949ChbLqIyCIU3EO3Z09/00WkUaF07Cq4h27Dhv6mi0hjQurYVXDvCOXf7UI33QSrV7902urV2XSRgIX6kapSSB27Cu4Q1r/bha68EiYns73FLHucnEzvZKq0SsgfqSqF1LGb3sBh3UYsmpvL9sKRkWyrTEyoK0NaK4RBuVL+SNW5/iu5QXYrLZZOhPTvVqSAUDLmlD9SodxgPK3Mfal0IoR0R6SgkDJmfaSqp8y9Y6l0IpR/tyIFhJQx6yPVrLQyd1A6URWt18LKWoXaFOlQ5j6f0onyhVLobbEyV6F2cYEUg7uUr5/m3gCbn0OYpZD6oyUOCu5SXK+F3gAz/FBmKaRaucQhvZq7VKOXQm9IrRwBzpJq5dIv1dyler0UegNMT0OaJdXKpUzK3KVeAaanAc6SJKCM/U6Zu4QjwPQ0wFmSyNVxrkfBXSRRIXQJpaqO7igFd5EEhdIllKo6zvUUCu5mttbM7jWzn5vZY2Z2rpmtM7PvmNnj+eMJZc2sSAxCyJjb0FcfwnqqSh33ey+aud8CfNvdXwe8CXgMuAHY6e6bgJ35cxEhnIw5pC6hbkJZT1Wq+lzPwMHdzI4DzgfuAHD3F9z9WeAyYHv+su3A5UVnUiQWoWTMdWSORYSyntqsSOZ+BjAL3GVmPzSz281sCDjZ3fcD5I8nlTCfIlEIKWMOuUsopPXUVgP3uZvZKPAD4Dx3f9jMbgGeBz7o7mvnve4Zd39Z3d3MxoFxgA0bNpw1MzMz0HwkSY3ZL1kF0K7Voc3XG62n5S3V576qwN/dB+xz94fz5/eS1defNrNT3H2/mZ0CHOj2y+4+CUxCdhFTgflIS6cYOTWVpTQhHlNXbOEqcIfvf789q6OTMcvStJ6KGbgs4+6/Avaa2WvzSRcBPwMeBK7Op10NPFBoDkMSwul7FSNfsgoeeij7Snh1iHRVJHMH+CCww8yOBZ4AriH7h/E1M7sW2ANcUfA9whBKxtwpRnbmI8Fi5PxVcP75L83cE1wdIl1pbJleafjAoLS55t5W2u3Co7FlyhDS6fuQ2xxqMn8VxLY6tl5zJ6uOn8HsKKuOn2HrNXc2PUtJ9J3HRsG9V6E3BksUtl5zJ9u+8i5efH4EWMGLz4+w7SvvajzA61RP+yi49yO2FFGCM3n/RXBkwVHhkaFseoNCOnCV3qjmLhIQs6N0z7mO4t5sLqaae3hUcxdpiZXH7e1rep104NouCu4iARn/k52wakFBe9XBbLpIHxTcRQIycdf7uG7sHlYeNwMcZeVxM1w3dg8Td72v6VmTllHNXSQhTdXNVa+vhmruItJYr7p65JsRZ3APYQwYkcA01auuHvlmxBfclSaIdNVUr7p65JsRX809pDFgRAKjmntc0qq5K00QWVRTvepF31eV1v4VHfI3PJ0xYMbGlCaIRCCU0bbbJr7MHXQpnQxMGWJ4dEJ2MHEGd5EB6Fx8mFRpHUx8ZRkJTltOpi3MEMfGdC4+BKq0DkaZu1SqTdmwMsRwqdLaP2XuUqkms+F+jxiUIUpMlLlLpZrKhgc9YlCGKLFQ5l6GthSVG9BUNqz6uaROmXtRbSoqN6SJbFj1c0mdMveilCIGSfVzSZ0y96KUIgZL9XNJmTL3opQiikiAFNzL0EkRRUQCobLMoDQIiYgETMF9EOqQEZHAKbgPQsPUtZ4OvKRp7s7cC3NUdcMkBfdBqEOm1XTgJU1zd27bdRtbv7mV23bdVkmAV3AfRKdDZmJCdw5ooZQOvAY5QtFRTfUOHj7I1N4pRo4fYWrvFAcPl78TKrgPSk3UrZXKgdcgRyg6qqnH0DFDbD5tMzPPzbD5tM0MHVP+Tlj4BtlmthKYBp5090vN7HTgbmAd8Ahwlbu/sNTfKPUG2SI9SGE4oEHuFa/7y9fH3Tl4+CBDxwxhA+6EVd8g+3rgsXnPPwN8zt03Ac8A15bwHiKlSuHAa5AjlFSOakJgZqw5ds3AgX3Zv18kczezU4HtwE3AR4E/BmaB33X3I2Z2LvBpd/+jpf5OMpl7CumiBGWQXU67aXtUmbl/HvgEcDR//irgWXc/kj/fB6wv+B5xUDFTGjDIEUoKRzUpGDi4m9mlwAF33zV/cpeXdo1iZjZuZtNmNj07OzvobLRHSi0aItK4Ipn7ecA7zeyXZCdQLyTL5NeaWWfMmlOBp7r9srtPuvuou48ODw8XmI2WUDFTRGpUuFsGwMwuAP4y75b5OnCfu99tZl8EHnX3iaV+XzV3EVlKGZ0lMaq6W2ahTwIfNbPdZDX4Oyp4j3ZSMVOkb3VczRmjUoK7u3/P3S/Nv3/C3c92999z9yvc/X/LeI9F3liX0kkt2rqrtXW+56vjas4YtfcKVXWfSE3auqu1db4XquNqzhi1N7ir+0RqUuWuVmVmHctHxMzYctYWJi6ZYMtZW1Rz71F7g7u6T6QmVe1qVWfWMX1Eqr6aM0aldMsUNXC3jLpPpCZV7Gp1jOOij0jc6u6WqY+6T6QmVexqdWTW+oikSzfIFmlI57YAY2PKrKV8Cu4iDepk1iJla3dZRoIUQ2+1NK/qe4zGTsFdShVLb7U0S1elFqfgLqXoZOtzc3H0VkuzdFVqcaq5S2GdbH1qCs49N+v8mJpqf2+1NKdzVerU3ildlTqgdve5SxAW9mvfemt2olAdILKYXkZ51EiQy4u3z12CsLBfe80a9VbL4nqtp+uq1GIU3KWwTr/2xET2qM+iLEX19HoouCei6vbEGK+EVEtnNTTKYz1Uc0/A/BOemzcru+6F1lm1VE8vh2ruiYtl6Nc6lbnOdDHOy6meXj0F9wTENPRrXcpaZ7oYR5qiskwiNPRr/8pYZ3MvzLH1m1sZOX6EmedmmLhkgjXHajAZKYfKMhLlCc+qlbHOdPJQmqLMXaRiMZ48jHGZ2kiZu0iDYjt5qPMI7aDg3nLqxZa66SKkdlBwbzENrytN0HmEdtCokC22sBd7bEx39ZHqmRlbztrC2JljqrkHTJl7i6l/XZoS23mEGClzbzHdYDkuMXSgxLAMsVDm3nLqX49DDB0oMSxDTBTcRQIQQwdKDMsQEwV3kQDE0IESwzLERFeoigQihnp1r8tQ1rLGsM6K0BWqIi0QQwdKL8tQVm1eNf6lDRzczew0M/uumT1mZj81s+vz6evM7Dtm9nj+eEJ5sytSDo2x3pyyavOq8S+tSOZ+BPiYu78eOAf4gJm9AbgB2Onum4Cd+XORYCjja1ZZtXnV+JdWWs3dzB4AvpB/XeDu+83sFOB77v7apX5XNfflpTQee9V1VI2x3jzV3MtRec3dzDYCbwYeBk529/0A+eNJZbxHylIaQ6aOrFoZX/PKOr8Qw3mKqhQO7ma2BrgP+LC7P9/H742b2bSZTc/OzhadjaildA/UOuqonbFRJi6ZYMtZWxQYJEqFgruZHUMW2He4+/355Kfzcgz544Fuv+vuk+4+6u6jw8PDRWYjeimNIVNXVq2MT2I3cM3dsk/FduA37v7hedP/Fvgvd7/ZzG4A1rn7J5b6W6q5Ly+lmvvRo0c5cOgAJ60+iRUr1K0rspilau5FBg47D7gK+LGZ/Sif9tfAzcDXzOxaYA9wRYH3kFxnDJnYuTuTj0wytXeKzadtVtlEZEADB3d3/xdgsU/dRYP+XUnbwpr72Jlj0XWypN7hIfXQMa8EJfZOFvXYS100nrsEJfa7/KRwZNIvHclUQ5m7BCfmTpayjkx2bNvKxo+vYsWnjY0fX8WObVtLntN66EimOu0O7u4wNxf3VT010ngr1Sujx37Htq2MP7mNmTUv4gYza15k/MltrQzwGh+mOu0N7ildtlkDZVD1KXpkcuMTkxw65qXTDh2TTW+b2M+xNKm9wT2lyzZrMEgGpUy/GXuGXuxreshSvlq46s9Pe4N7Spdt1qDfDEqZfnM2HFzZ1/TQxXyOZTF1fH7aG9zNYMsWmJjIHhPaMarQbwalWunL1XUkc9MZ46w+/NJpqw9n06U/TR191vH5aW9wh99etqnAXop+MijVSl+qziOZK6+bYHL9dYzMrcQcRuZWMrn+Oq68bqKy94xRk0efdXx+dA9VGViI/clNzZPGiG+fprdZGfuq7qEqlQitVhp7JiblanqbVf35UeYu0YghE5N6tX2bKXOXJMSeiUn5Yt5mytwlKm3PxET6ocxdkhFzJla2Im2AuoAtfAruIgkqcvJZF7C1g4K7SIKKXESjC9jaQcFdJEFFTj43feJaeqMTqiKJKnLyWSeuw1DVDbJFpMU6J5/r/l2ph8oyIiIRUnAXEYmQgrtISdT7LSFRcBcpgXq/JTQK7iIlUO93OXT0Ux4Fd5ESqPe7OB39lEutkJK8km6awJaztjB25ph6vwe08Ohn7MyxoNstQ+/1V+YuSSszW9SgZcW06einDUcZCu6JUC2zu1Br5Slur35v0t6kUPeb+RTcE9CGLKMpIWaLKW+vthz9hLjfLKSxZRJQ1e3nQq859iq05Wj6doHSmxD2G92sI3FVZBkxZZehZYttyAolvP1moUoydzN7O3ALsBK43d1vXur1ytyrV3aW0U92GUKGU7eiy9ymddamea1KU+ug1szdzFYCtwIXA28A3mNmbyj7faQ/ZWcZvWaXMWX4vSpjmUPPCjtS3L4LhboOqijLnA3sdvcn3P0F4G7gsgreRxrUa2dDG7oKypbSMqe0rIsJdR1UEdzXA3vnPd+XT5PI9JJdplg/TmmZU1rWxYS6DkqvuZvZFcAfufv78+dXAWe7+wcXvG4cGAfYsGHDWTMzM6XOh4QjxZpsSsuc0rIuJomaO1mmftq856cCTy18kbtPuvuou48ODw9XMBsSirbUj8uU0jKntKyLCXEdVBHc/w3YZGanm9mxwLuBByt4HxERWUTpA4e5+xEz+wvgn8haIe9095+W/T4iIrK4SkaFdPdvAd+q4m+LiMjydIWqiEiEFNxFRCKk4C4iEiEFdxGRCAUx5K+ZzQL9XMV0IvDrimYnZFrutGi50zLIco+4e9cLhYII7v0ys+nFrsqKmZY7LVrutJS93CrLiIhESMFdRCRCbQ3uk03PQEO03GnRcqel1OVuZc1dRESW1tbMXUREltC64G5mbzezfzez3WZ2Q9PzUxUzO83Mvmtmj5nZT83s+nz6OjP7jpk9nj+e0PS8ls3MVprZD83sG/nz083s4XyZ78lHG42Kma01s3vN7Of5Nj83kW39kXz//omZfdXMfifG7W1md5rZATP7ybxpXbevZf4+j3GPmtlbBnnPVgX3xO7PegT4mLu/HjgH+EC+rDcAO919E7Azfx6b64HH5j3/DPC5fJmfAa5tZK6qdQvwbXd/HfAmsuWPelub2XrgQ8Cou7+RbBTZdxPn9v4S8PYF0xbbvhcDm/KvcWDbIG/YquBOQvdndff97v5I/v1/k33Y15Mt7/b8ZduBy5uZw2qY2anAJcDt+XMDLgTuzV8S4zIfB5wP3AHg7i+4+7NEvq1zq4BXmtkqYDWwnwi3t7s/BPxmweTFtu9lwJc98wNgrZmd0u97ti24J3l/VjPbCLwZeBg42d33Q/YPADipuTmrxOeBTwBH8+evAp519yP58xi3+RnALHBXXo663cyGiHxbu/uTwGeBPWRB/TlgF/Fv747Ftm8pca5twb3bPayibvcxszXAfcCH3f35puenSmZ2KXDA3XfNn9zlpbFt81XAW4Bt7v5m4CCRlWC6yWvMlwGnA68GhshKEgvFtr2XU8o+37bg3tP9WWNhZseQBfYd7n5/PvnpziFa/nigqfmrwHnAO83sl2QltwvJMvm1+WE7xLnN9wH73P3h/Pm9ZME+5m0N8DbgF+4+6+6HgfuBzcS/vTsW276lxLm2Bfdk7s+a15rvAB5z97+b96MHgavz768GHqh73qri7n/l7qe6+0aybfvP7n4l8F3gT/OXRbXMAO7+K2Cvmb02n3QR8DMi3ta5PcA5ZrY63987yx319p5nse37IPDevGvmHOC5TvmmL+7eqi/gHcB/AP8J3Nj0/FS4nH9Idij2KPCj/OsdZDXoncDj+eO6pue1ouW/APhG/v0ZwL8Cu4GvA69oev4qWN7fB6bz7f2PwAkpbGvgb4CfAz8B/gF4RYzbG/gq2XmFw2SZ+bWLbV+yssyteYz7MVk3Ud/vqStURUQi1LayjIiI9EDBXUQkQgruIiIRUnAXEYmQgruISIQU3EVEIqTgLiISIQV3EZEI/R/guYv3gP78rwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "color = [\"r\", \"g\", \"b\"]\n",
    "for k in range(K):\n",
    "    plt.scatter(X[np.where(R_nk[:,k] == 1)][:,0],\n",
    "                X[np.where(R_nk[:,k] == 1)][:,1],\n",
    "                c=color[k], s=5, alpha=0.5)\n",
    "    plt.scatter(C[k:,0], C[k:,1], c=color[k])\n",
    "plt.show()"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}