zokuwaka12_algorithm2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg\n",
    "np.random.seed(13)\n",
    "from collections import Counter\n",
    "from math import pi,sqrt, gamma, isnan\n",
    "from scipy import special\n",
    "from scipy.stats import wishart\n",
    "from scipy.misc import factorial\n",
    "from scipy  import linalg\n",
    "import copy\n",
    "import matplotlib as mpl\n",
    "from hoge import DirichletProcessMixtureGaussianDistribution\n",
    "from sklearn  import mixture\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEWCAYAAAA6maO/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYVGXiB/DvMCACaQ4Ccp81BUVTLFfEUiNdQ01NvJQo\nLal4QcjLes9oVTRvq5FBpYEtu1r9yuhi6mrbk6wbCFYquq7bqD8Gb9yUMWEUdZjfHy7zc2S4OjNn\n3uH7eR6fx86cmfmOAV/Oe97zHplGo9GDiIjIxjlIHYCIiKgpWFhERCQEFhYREQmBhUVEREJgYRER\nkRBYWEREJAQWFhERCcEihZWTk4Po6Gj06NEDCoUCH3/8sdHjc+bMgUKhMPrz3HPPWSIKERHZCUdL\nvGhVVRV69uyJ6OhoxMfHm9zn2Wefxfbt26HX37tu2cnJyRJRiIjITliksIYNG4Zhw4YBuHc0ZUqb\nNm3g4eFhibcnIiI7JNk5rCNHjiAoKAi//e1vMW/ePJSXl0sVhYiIBGCRI6zGDBs2DGPGjIFSqURR\nURGSk5MxZswYZGdnc2iQiIhMkqSwoqKiDH8PCQlBaGgoevXqhQMHDmDUqFFSRCIiIhtnE9Pavb29\n4evri/Pnz0sdhYiIbJRNFFZ5eTmuXLmCTp06SR2FiIhslEUKq6qqCidPnkRBQQFqampw8eJFnDx5\nEhcvXkRVVRWSkpJw9OhRFBUV4fDhw5g8eTK8vLyEHg5UqVRSR2gS5jQv5jQfETICzCklixTWsWPH\nMHjwYERERODWrVtYt24dnnnmGaxbtw5yuRynT5/GlClT0K9fPyQkJCA4OBgHDx6Em5ubJeIQEZEd\nsMiki4EDB6KioqLexz///HNLvC0REdkxmziHRURE1BgWFhERCYGFRUREQmBhERGREFhYREQkBBYW\nEREJgYVFRERCYGEREZEQWFhERCQEFhYREQmBhUVEREJgYRERkRBYWEREJAQWFhERCYGFRUREQmBh\nERGREFhYREQkBBYWEREJwVHqAESiUavV2JS+E+U378L57k2sXhAPpVIpdSwiu8fCImoGtVqN2ORU\n6CITIHd2ha5ai9jkVGQmJbK0iCyMQ4JEzbApfaehrADcK63IBGxK3ylxMiL7x8Iiaobym3cNZVVL\n7uyK8pt3JUpE1HqwsIiawcPFEbpqrdE2XbUWHi4cXSeyNBYW0QPUajUSk9Zi0qJVSExaC7VabXhs\ncVwM5AfSDKWlq9ZCfiANi+NipIpL1Grw10Ki+zw4qeLSA5MqlEolMpMSjWcJcsIFkVWwsIju09Ck\nitTkFQAApVJp+LtKpWJZEVkJhwSJ7sNJFUS2i4VFdB9OqiCyXSwsovtwUgWR7bJIYeXk5CA6Oho9\nevSAQqHAxx9/XGefdevWISQkBD4+Phg1ahTOnDljiShEzVI7qSK0YBe8//EeQgt2cRULIhthkXGO\nqqoq9OzZE9HR0YiPj6/zeEpKCt577z28++676Nq1KzZs2ICoqCj8+OOPcHNzs0Qkoia7f1IFEdkO\nixxhDRs2DK+//jrGjBkDmUxW5/H3338fCxYswKhRo9C9e3e89957qKysxO7duy0Rh4iI7IDVz2EV\nFhaipKQEzz77rGFb27Zt8dRTTyEvL8/acYiISBBWL6zS0lLIZDJ4enoabff09ERpaam14xARkSA4\nV5fs1v33rfJwccTiuBhOniASmNULy8vLC3q9HmVlZfDz8zNsLysrg5eXV4PPValUlo73UGw9X63W\nkPPS5St4LSML8tHzDUssRSe9hbmjBmJ/zo+ouK2Hoo0MU6NGws/XR7Kc1iRCThEyAsxpLkFBQc3a\n3+qF9Zvf/AadOnXC999/jz59+gAAbt26hdzcXKxZs6bB5zb3w1mTSqWy6Xy1WkvOt//yqaGsgHur\nVdx+6kUkf7wD7pNfh9zZFRXVWqz8a9pDTVtvLf+e1iBCRoA5pWSRc1hVVVU4efIkCgoKUFNTg4sX\nL+LkyZO4ePEiACA+Ph4pKSnYs2cPTp8+jTlz5uCRRx7B+PHjLRGHWiFTSyxdz9tvKCuAN18kEo1F\njrCOHTuG0aNHG6a0r1u3DuvWrUN0dDTS0tIwb9483Lp1C0uWLIFGo0Hfvn2RlZXFa7DIbDxcHHGp\nWmtUWjW6O1wnkEhgFimsgQMHoqKiosF9li5diqVLl1ri7YmwOC7G6DYhumotHK6o7i21dF9pcZ1A\nInFwLUGyS6aWWMpYtYjrBBIJjL9akt0ytcRSpp+v8VR3rhNIJAwWFrUqXCeQSFwcEiQiIiGwsIiI\nSAgcEqRWics2EYmHR1jU6qjVasQmp+JE7ykoHhyPE72nIDY5FWq1WupoRNQAFha1OpvSdxquzwK4\n4gWRKFhY1OqYWraJK14Q2T4WFrU6Hi6OhouHa3HFCyLbx8KiVmdxXAxXvCASEH+lpFandtkmrnhB\nJBYWFrVKXPGCSDwcEiQiIiGwsIiISAgsLCIiEgILi4iIhMDCIiIiIbCwiIhICCwsIiISAguLiIiE\nwMIiIiIhsLCIiEgILCwiIhICC4uIiITAwiIiIiGwsIiISAi8vQgRkeDUarXx/d3s9GakPMIiIhKY\nWq1GbHIqTvSeguLB8TjRewpik1Nx6fIVqaOZnSSFtX79eigUCqM/3bt3lyIKEZHQNqXvhC4yAXJn\nVwCA3NkVusgEfPjFPomTmZ9kQ4LBwcHYu3cv9Ho9AEAul0sVhYhIWOU37xrKqpbc2RWa23qJElmO\nZIUll8vh4eEh1dsTEdkFDxdHXKrWGpWWrlqLDm1kEqayDMnOYanVaoSEhCA0NBTTp09HYWGhVFGI\niIS1OC4G8gNp0FVrAdwrK/mBNEyNGilxMvOT5AirX79+ePfddxEUFISysjJs2rQJkZGRyMvLQ4cO\nHaSIREQkJKVSicykRONZgkmJuH37ttTRzE6m0WgkH+jUarUIDQ3FggULMGfOnHr3U6lUVkxFRESW\nFBQU1Kz9beI6LFdXV3Tv3h3nz59vcL/mfjhrUqlUNp2vFnOaF3OajwgZAWlymrrOSqlUNvgcUf49\nm8MmCuvWrVtQqVQYPHiw1FGIiGxK7XVWtVPXL1VrEZucisykxEZLqznv0dxClIIkhZWUlIThw4fD\n39/fcA5Lq9UiOjpaijhERDarvuusVr69DW5ubigqvYaSK5fhE/gb+CseaXbZWKMQzUWSWYKXL1/G\njBkzEBYWhtjYWLRt2xZ///vf4e/vL0UcIiKbVd91Vjn/uYCjfs/i37/q0eblddA8t8CwyoVarW7y\n69dXiJvSd5r1c5iDJEdYGRkZUrwtEZFw6rvOSt/RH9fz9sN73FyTZTPv9y826fXLb96F7kYFyvZl\nQK+vgUzmAPdB41B+865FPs/D4FqCREQ2zNR1VmW7kuER8SL0+hqTR1/NKRvnO1UoO5AJz5HT4TNu\nHjxHTkfZgUw436ky6+cwB5uYdEFERKaZus7qhrIjzrdTQCZzuHeh8ANHXx4uTf/R7iB3hPeEOUZH\nad4T5sMhd7vZP8vDYmEREdk4pVKJ1OQVhv+unSjxaP8oFGdtNQwL1q5y0ZwLh286OJs8Srvp4GzW\nz2AOLCwiIsEolUqsnTYOf1i/FY63a3A5ZRZcFR6Q19xFaFDzZvbVd46sOUdp1sJzWEREglGr1Vix\nIwttXl6HTnPege/8bbjl5gnXCctwfuCcZt0Pq761CG3xJpAsLCIiwZiaiu49YT6uHc5q9v2was+R\nhRbsgvc/3kNowS6bvAYL4JAgEZFw6rs2S6+vMfy9OffDevAcma3iERYRkWA8XBwNQ3i1dNVayGQO\nhr/zflhERGQWarUaiUlrMWnRKiQmrW3W6hSmzjsV706B+6BxvB8WERGZz8Ou3/fgtVkuNdV4zLsN\nqk9/bdf3w2JhERFZWUPr9zX1XFJj553s8f6BHBIkIrKy+iZN2OL6fbaEhUVEZGX1TZqwxYt1bQkL\ni4jIyqx5se7DTO6wNaxzIiIrM7Wg7WILXKwr0s0Zm4KFRUQkAWtcrGuOyR22hEOCRER2yt4md7Cw\niIjslL1N7mBhEZFdsqfJBi0l0krsTSFmzRIRNeD+yQa6GxX4+dCn+HbO63iqWwBWzpsl5ISDlrDW\n5A5rYWERkbDUarXxD+O4GCiVSsNkA92NCpR9+1fDHXnPCj5LriVEWYm9KVhYRCQkU1O2R85PxvbF\n0w2TDcr2ZRjKCjCeJbc4LsZk2ZHtYmERkZBMTdl2n/w6pv9xIUK7BEJXrcXdKg3K9mVAr6+BTOYA\n90Hj0MbDDxcrKjH5jRQ4jJpnKLvJb6Tgo9XzWVo2jIVFREKqPYq6XX4J1w5nGUrpptwFR06egf4/\nC+HQ1hU+ExfeO7Kq1qI4ayvcn5mAwmPH0KbLE/C4UQG5s+u90hs1Dyvf3oYPt7wp9UejerCwiEgI\nD56vcr5ThZuXz+Ja9m7DsJ+uWosLGSvgPmomtL/8DM+R041vIz9uLtRpC+DoGQjFU6NR9u1f4Tns\nZbTx8IPc2RUnCq9I/CmpIZzWTkQ2r/Z81YneU1A8OB4nek/BWc0dFP81uc45qoDpa1F+8K/QVd8y\nedGss28X+E15Ddfz9sN73FxcO5wF4N6U75pb2jrvTbaDhUVENs/U+SqHUfPQXtHRZCm1DQjG7VK1\nyYtmb5cWoXTvB7hZdAa6GxWoLinC5d1boE5bgK4+7lb7TNR8HBIkIptX3xJDuHvr3sWw9z2mq9bC\nQe4EJ/dOhnNW1/P2o0Z3B9rzp+D1/HS07/k0bl4+i4t/WQ0HZ1c4yJ3gM2kxruR+BrVazYkXNopH\nWERk8+pbYsipRoeLf1lttJJDcdZWPNp/BBzdOqBd70Eo+TINniOnw3fCH9B5XipuFBxG5S8/4Vr2\nbigT3sJvEt6C58jpuJa9G3cHTMSm9J1SfERqAkkLKz09HaGhofD29kZERARyc3OljENENqq+JYb8\nOz+Gjs9MhPrdP+Dy7i0o25cB92cmoPzr9xDcRgun3E8RMH1tnYkXZX/7sM65L+9xc3E9b7+wC8O2\nBpIVVlZWFpYvX45Fixbh8OHDCAsLw8SJE3Hp0iWpIhGRjapdYii0YBe8//EeQgt2ITMpEYFe7nBR\nhsD/5SQ4yJ2g19egImcPBgb54Ottm9E5uLvJoUSHNi4mt9fo7gi7MGxrINn/mXfffRcxMTF4+eWX\nAQAbN27Ed999hx07diApKUmqWERko0wtMbQ4LgaxyalAZAK8o141HHmt/kMigHtDiZdMnOOquX3T\n9LmvKyosTt5inQ9EzSbJEdadO3dw/PhxREREGG0fMmQI8vLypIhERAKq78irdtKEqaHEax+tQbB3\nB1z7aE2d7RmrFnHChQ2T5Ajr6tWr0Ol08PLyMtru6emJ7OxsKSIRkaAaWtzV5GrlKUlQKpV1F879\n73ayXTKNRqO39psWFxcjJCQE+/btw4ABAwzbN27ciN27dyM/P9/k81QqlbUiEhGRhQUFBTVrf0mO\nsDp27Ai5XI7S0lKj7WVlZXWOuu7X3A9nTSqVyqbz1WJO82JO8xEhI8CcUpLkHJaTkxP69OmDQ4cO\nGW3//vvvER4eLkUkIiKycZLNEkxISMDs2bPxxBNPIDw8HBkZGSgpKcErr7wiVSQiIrJhkhVWVFQU\nKioqsHnzZpSUlCAkJASfffYZ/P39pYpEREQ2TNIr5KZNm4Zp06ZJGYGIiATBtQSJiEgILCwiIhIC\nF80iElidi1/jYnjxK9ktFhaRoGrvwlt7Y8NL1VrEJqcaLU3UnNdi8ZGtY2ER2bj6ysTUXXh1kQl4\nY0sa2rV/1Gj/xl7fXMVHZEk8h0Vkw2rL5ETvKSgeHI8TvacgNjkVarXa5F14dTcqkKe+Wmf/S5ev\n1Pse9RUfb2RItoZHWERW1tDw24OPVVVV1Vsmpm6dUX7oU3SaklRn/w+/eB8Rzww2mae+28/zRoZk\na1hYRFbU0PAbgDqPlf45CT4Rpstkw6uvGO2vq9ZCdvWiyfLR3K5/jev67hll6kaGPNdFUmJhEVlR\nY8NvDz6m7+hv8kaDHi6OJm+dURbggcsm9tdVapCYtNZk0dTeBPH+4pMfSMPi/5ZoLZ7rIqnxHBaR\nFTU0/GbqMY+IF1G2K9noRoPyA2mGiRS194L65E9/RGryCrR1dkLx7hSj/Yt3p+D0+SKT58FqX6Oh\nmyDW4rkukhqPsIisqLHhtwcfk7dToL+yI9oV7Pr/o6MGjmiqndzgGRmNsn0Z0OtrIJM5wDMyFuWH\n/sdk0dTe+LChmyDW4rkukhoLi8iKGht+M/XY6mYMuXm4OOJSOwW8o141bNNVa+EgdzLaryVF05xz\nXUSWwCFBIitqaPitqUNzDVkcFwP5gTSjIcFrH63Bo/1HGO3XkqIx9dr3D08SWRp/NSKysoaG35oy\nNNfYaz84ESN68XQsfv9/oBs9v8FJFS157YaGJ4nMjYVFZGdMld6bd+5g99+bdh6sua9NZC0sLKJW\nwM/Xh0VDwuM5LCIiEgKPsIjsCFeiIHvGwiKyE40t+0QkOg4JEtkJrkRB9o6FRWQnuBIF2TsWFpGd\n8HBxNFzUW4srUZA9YWER2QmuREH2jr96EdmJhlaiUKlUUscjemgsLCI7wpUoyJ5xSJCIiITAwiIi\nIiGwsIiISAgsLCIiEoIkhfX8889DoVAY/ri7uyMuLk6KKEREJAhJZgnKZDLExMTgj3/8I/R6PQCg\nbdu2UkQhIiJBSDat3cXFBR4eHlK9PRERCUayc1hZWVno0qULBgwYgKSkJFRWVkoVhYiIBCDJEdaL\nL76IgIAAeHt748yZM1i5ciVOnz6Nzz//XIo4REQkAJlGo9Gb44XWrFmDzZs31/9GMhn27NmDp59+\nus5jx44dw5AhQ5CdnY3evXvX+xpcXoaIyH4EBQU1a3+zFVZFRQWuXr3a4D7+/v4mJ1fo9Xp4enoi\nPT0dY8eONUccq1OpVM3+x5cCc5oXc5qPCBkB5pSS2YYEa6eot8SpU6eg0+nQqVMnc8UhIiI7Y/Vz\nWIWFhfj000/x3HPPwd3dHWfOnEFSUhL69OmD8PBwa8chIiJBWL2wnJyckJ2djW3btqGqqgp+fn6I\njIzEkiVLIJPJrB2HiIgEYfXC8vPzw969e639tkREJDiuJUhEREJgYRERkRBYWEREJAQWFhERCYGF\nRUREQmBhERGREFhYREQkBBYWEREJgYVFRERCYGEREZEQWFhERCQEFhYREQmBhUVEREJgYRERkRBY\nWEREJAQWFhERCYGFRUREQmBhERGREFhYREQkBBYWEREJgYVFRERCYGEREZEQWFhERCQEFhYREQmB\nhUVEREJwlDqApajVamxK34nym3fh4eKIxXExUCqVUsciIqIWsssjLLVajdjkVJzoPQXFg+NxovcU\nxCanQq1WSx2NiIhayC4La1P6TugiEyB3dgUAyJ1doYtMwKb0nRInIyKiljJ7YWVmZmL06NFQKpVQ\nKBS4cOFCnX00Gg1mzpyJwMBABAYGYtasWbh+/brZMpTfvGsoq1pyZ1eU37xrtvcgIiLrMnthabVa\nDB06FMuXL4dMJjO5T1xcHE6dOoUvvvgCWVlZKCgowOzZs82WwcPFEbpqrdE2XbUWHi52e8qOiMju\nmf0neHx8PADg+PHjJh//5Zdf8N133+HgwYPo27cvAOCtt97CiBEjcO7cOXTp0uWhMyyOi0Fscqph\nWPDm5bO4/vlmXOwWjMSktYYJGJyYQUQkDqsfcuTn56Ndu3bo16+fYVt4eDjc3NyQl5dnlsJSKpXI\nTErEpvSdKCq9hrLS6/CcuRkaZ1dcrdYiNjkVa6eNw4odWYZSu/Tf7ZlJiSwtIiIbZPVJF6WlpejY\nsWOd7R4eHigtLTXb+yiVSqQmr0CglzvcJ79eZwLGH9Zv5cQMIiKBNKmw1qxZA4VCUe8fd3d3/PDD\nD5bO2iL1TcDQypw5MYOISCBNGhJMSEjApEmTGtzH39+/SW/o5eWFq1ev1tleXl4OLy+vBp+rUqma\n9B73c757E7pqrVE56aq1aHO70uR257s3W/Q+Lc0nBeY0L+Y0HxEyAsxpLkFBQc3av0mFVXskZQ5h\nYWGorKzE0aNHDeex8vLyoNVq0b9//waf29wPBwCrF8QbTcDQVWshP5CGrUkLsWJHWp3tq1t4Dkul\nUrUon7Uxp3kxp/mIkBFgTimZfdJFaWkpSkpKoFKpoNfrcebMGWg0GgQEBKBDhw4IDg7G0KFDMX/+\nfKSkpECv12PBggUYPny4WSZcPOj+CRiG2YD/LaVMP1+T24mIyPaYvbB27NiBDRs2QCaTQSaT4aWX\nXgIApKWlITo6GgCQnp6OJUuWYPz48QCAkSNHYuPGjeaOYlA7AaOp24mIyPaYvbCWLVuGZcuWNbjP\no48+im3btrXo9XntFBFR6yTUWoJc1JaIqPUSqrC4qC0RUeslVGFxUVsiotZLqMLiorZERK2XUIW1\nOC4G8gNphtKqvXZqcVyMxMmIiMjShDo0aeiaKiIism9CFRbAa6eIiForoYYEiYio9WJhERGREIQb\nErR1XImDiMgyeIRlRlyJg4jIclhYZsSVOIiILIeFZUZciYOIyHJYWGbElTiIiCyHhWVGXImDiMhy\n+Ku/GXElDiIiy2FhmRlX4iAisgwOCRIRkRBYWEREJAQWFhERCYGFRUREQmBhERGREFhYREQkBBYW\nEREJgYVFRERCYGEREZEQWFhERCQEFhYREQmBhUVEREIwe2FlZmZi9OjRUCqVUCgUuHDhQp19evXq\nBYVCYfjj7u6O1atXmzsKERHZEbOv1q7VajF06FA8//zzeO2110zuI5PJsGzZMkyfPh16vR4A4Obm\nZu4oRERkR8xeWPHx8QCA48ePN7ifm5sbPDw8zP32RERkpyQ7h5WamorHHnsMgwYNwubNm3Hnzh2p\nohARkQAkuYHj7Nmz0bt3b7i7u+Onn37CypUrUVRUhLfffluKOEREJACZRqPRN7bTmjVrsHnz5vpf\nRCbDnj178PTTTxu2HT9+HEOGDMGJEycQEBDQ4Ot/9dVXmDp1Ks6fP48OHTo0Iz4REbUWTTrCSkhI\nwKRJkxrcx9/fv8UhnnzySej1epw/fx5PPvlki1+HiIjsV5MKq3b6uaUUFBRAJpOhU6dOFnsPIiIS\nm9nPYZWWlqKkpAQqlQp6vR5nzpyBRqNBQEAAOnTogKNHj+Lo0aMYNGgQ2rdvj59//hkrVqzAyJEj\n4efnZ+44RERkJ5p0Dqs51q9fjw0bNkAmkxltT0tLQ3R0NE6cOIFFixZBpVLh9u3bCAgIwPjx4zF3\n7ly0bdvWnFGIiMiOmL2wiIiILMFm1xJsbImnoqIivPrqq+jTpw98fHzQp08frF69Grdu3bKpnACg\n0Wgwc+ZMBAYGIjAwELNmzcL169etmvNBpaWlmDlzJrp16wZfX18MHDgQn332maSZ6vPTTz8hKioK\n/v7+CAgIwPDhw1FRUSF1rHpNmDABCoUCX3/9tdRRjGg0GixZsgRhYWHw8fHB448/joULF9rEv2V6\nejpCQ0Ph7e2NiIgI5ObmSh3JyJYtWzBkyBAEBgaia9eumDRpEv79739LHatBW7ZsgUKhwJIlS6SO\nUkdJSQni4+PRtWtXeHt7Y8CAAcjJyWn0eTZbWLVLPC1fvrzO8CIAqFQq1NTUICUlBUeOHMGmTZvw\nySefYPny5TaVEwDi4uJw6tQpfPHFF8jKykJBQQFmz55t1ZwPmjVrFs6ePYtPPvkEubm5mDRpEmbN\nmmVzPyh+/PFHjBs3DoMHD8Z3332H7OxsJCYmwtFRkksIG/XOO+9ALpfX+7UgpStXrqC4uBjJycnI\nzc3F9u3bkZOTg7i4OElzZWVlYfny5Vi0aBEOHz6MsLAwTJw4EZcuXZI01/1ycnIwY8YMHDx4EHv2\n7IGjoyPGjh0LjUYjdTSTjh49iszMTDz++ONSR6nj+vXriIyMhEwmw+7du5Gfn48NGzbA09Oz0efa\n/JBgc67nysjIwJtvvolz585ZKd3/qy/nL7/8gv79++PgwYPo168fAODIkSMYMWIEfvzxR3Tp0sXq\nWYF7lyFs3LgRkydPNmzr1asXZs2ahcTEREkymRIZGYnBgwdjxYoVUkdp1M8//4zf//73yM7ORteu\nXZGZmYkxY8ZIHatB3377LSZNmgS1Wo1HHnlEkgy/+93v0KtXL7z11luGbX379sXYsWORlJQkSabG\nVFVVITAwEB999BEiIyOljmPk+vXriIiIwDvvvIP169ejR48e2Lhxo9SxDFavXo3c3Fzs37+/2c+1\n2SOslvj1119t7sLj/Px8tGvXzlBWABAeHg43Nzfk5eVJlmvAgAH48ssvUVFRAb1ej7179+LatWuI\niIiQLNODysvLkZ+fDy8vL4wYMQJBQUEYMWIEsrOzpY5Wx40bNzBjxgxs3boVHTt2lDpOk/36669w\ndnaGq6urJO9/584dHD9+vM7X3ZAhQyT9/mjMjRs3UFNTY3M/bwBg/vz5iIqKwsCBA6WOYtK+ffvQ\nt29fTJs2DUFBQRg0aBA++OCDJj3XbgqrqKgIqampkg9vPKi0tNTkDzAPDw+UlpZKkOieHTt2AAAe\ne+wxeHl5Yfbs2UhPT7epIYTCwkIA92aevvzyy8jKysJTTz2F8ePH41//+pe04R6wcOFCDBs2DEOG\nDJE6SpNpNBq8+eabiI2NhYODND8Krl69Cp1OBy8vL6Ptnp6ekn5/NGbZsmUIDQ1FWFiY1FGMZGZm\norCwEK+//rrUUepVWFiIjIwMdO7cGVlZWYiPj8eqVauQnp7e6HOteiKgJUs8NUVpaSkmTpyIoUOH\nGlaLfxiWymlpzcmdnJyMa9eu4euvv4a7uzv27t2LWbNmYf/+/ejZs6dN5HRycgIATJ061TB02atX\nLxw+fBgffvgh/vSnP9lEzgsXLuDUqVM4dOiQRfPUpyVfr1VVVYiOjoafnx9WrVpljZh247XXXkN+\nfj7+9re/2dS5yrNnzyI5ORkHDhyQ7BeQpqipqUHfvn0Nw729evXCuXPnkJ6e3ugBh1ULyxJLPJWU\nlOCFF17erkauAAAD3UlEQVRAz5498f777z9MPANz5vTy8sLVq1frbC8vL6/zW+XDamruwsJCfPDB\nB/jhhx/Qo0cPAEDPnj2Rk5OD7du3W3wR4qbmLCkpAQB069bN6LFu3bqZnI1pbk3J6efnh127duE/\n//kPfH19jR6bOnUqwsLCWjRW3xzN/XqtqqrChAkT4ODggE8++QRt2rSxaL6GdOzYEXK5vM7RVFlZ\nmdm/P8xh+fLl+PLLL/HNN98gMDBQ6jhG8vPzce3aNfTv39+wTafTIScnBx9++CEuX75s+CVQSp06\ndUJwcLDRtuDgYGzbtq3R51q1sMy9xFNxcTHGjBmDHj16ID093Wy/VZgzZ1hYGCorK3H06FHDeay8\nvDxotVqjLyxzaGpurVYLmUxW599LLpejpqbGrJlMaWpOpVIJHx8fqFQqo+1nz561ytBlU3O+8cYb\nmDt3rtG2AQMGYO3atRgxYoSl4hk05+u1srISEydOBAB89tlnkp27quXk5IQ+ffrg0KFDeOGFFwzb\nv//+e4wdO1bCZHUtXboUX331Fb755hvJJks1ZNSoUXXWYp0zZw66du2KhQsX2kRZAffO4T/4Pa1S\nqRqdVAdIdHuRpmhsiafi4mI8//zz8PX1xdq1a1FeXm54roeHh9UOiRvLGRwcjKFDh2L+/PlISUmB\nXq/HggULMHz4cMm+6IODg9G5c2csXLgQycnJcHd3x549e3Do0CF8/PHHkmSqz6uvvor169ejZ8+e\n6N27N7KysvDTTz81OARmbd7e3vD29q6z3dfXF0qlUoJEplVWViIqKgpVVVXYtWsXKisrUVlZCeBe\n6Un1Ay0hIQGzZ8/GE088gfDwcGRkZKCkpASvvPKKJHlMWbRoET799FPs2rUL7du3NxwRurm52czd\n0tu3b4/27dsbbXN1dUWHDh3qjFJIac6cOYiMjMTmzZsxbtw4nDhxAtu3b8fKlSsbfa7NTmtvbImn\njz76qM70a71eD5lM1qQp8NbKCdybZrpkyRLD0NDIkSOxcePGOl9c1vS///u/WLlyJY4cOYKqqip0\n7twZiYmJeOmllyTLVJ+tW7figw8+QEVFBbp374433ngDgwcPljpWg9zd3fHnP//Zpqa1//Of/6yT\np/Z7Rupzsjt27MDbb7+NkpIShISEYN26dQgPD5csz4MUCoXJ81VLly7F0qVLJUjUNKNHj0ZISIhN\nTWsH7l1OsWrVKpw7dw7+/v6YOXMmZsyY0ejzbLawiIiI7me7U0mIiIjuw8IiIiIhsLCIiEgILCwi\nIhICC4uIiITAwiIiIiGwsIiISAgsLCIiEgILi4iIhPB/IQhApH95OXMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c8c3e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_size_list = np.array([20,10,10,5,5])\n",
    "sigma = np.diag([1,1])\n",
    "x_means = [1,-2,4,-10,-5]\n",
    "y_means = [1,-2,4,-10,10]\n",
    "\n",
    "sigmas = [sigma*0.1, sigma*0.4, sigma*0.9, sigma*2, sigma*0.85]\n",
    "\n",
    "x = []\n",
    "y = []\n",
    "for i in range(len(data_size_list)):\n",
    "    np.random.seed(13+i)\n",
    "    data = np.random.multivariate_normal([x_means[i], y_means[i]], sigmas[i], data_size_list[i])\n",
    "    x.extend(list(data[:,0]))\n",
    "    y.extend(list(data[:,1]))\n",
    "\n",
    "plt.plot(x,y,\"o\")\n",
    "data = []\n",
    "for x,y in zip(x,y):\n",
    "    data.append([x,y])\n",
    "X = np.array(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "from math import pi, gamma, inf\n",
    "import copy\n",
    "import numpy as np\n",
    "from numpy.linalg import inv, det\n",
    "from scipy.special import gamma\n",
    "from scipy.misc import factorial\n",
    "\n",
    "np.random.seed(13)\n",
    "\n",
    "class DirichletProcessMixtureGaussianDistribution():\n",
    "    def __init__(self, data):\n",
    "        self.alpha = 1.\n",
    "        self.beta = 1/3\n",
    "        self.v = 15\n",
    "        self.S = np.diag([1, 1]) * 0.1\n",
    "        self.p_max = -inf\n",
    "        self.X = data\n",
    "        self.n = len(data)\n",
    "        self.mu_0 = np.mean(data, axis=0)\n",
    "        self.K = data.shape[1]\n",
    "        self.d = data.shape[1]\n",
    "        self.drop_count = 0\n",
    "        self.stop_count = 20\n",
    "\n",
    "        self.class_list = np.zeros(self.n, dtype=np.int32)  # データの付与されたクラスのリスト，indexがXと対応\n",
    "        self.class_freq = Counter(self.class_list)  # クラスの頻度\n",
    "\n",
    "\n",
    "    def fit(self, ite=20):\n",
    "        for i_ite in range(ite):\n",
    "            previous_class_list = copy.deepcopy(self.class_list)\n",
    "            previous_class_freq = Counter(previous_class_list)  # クラスの頻度\n",
    "            for k in range(self.n):\n",
    "\n",
    "                x_k = self.X[k]\n",
    "                s_k = self.class_list[k]\n",
    "\n",
    "                self.class_freq[s_k] -= 1\n",
    "                self.class_freq += Counter()\n",
    "\n",
    "                pro = []\n",
    "                for s_i, s_i_freq in self.class_freq.items():\n",
    "                    # 12.47\n",
    "                    x_j_index = np.where(self.class_list == s_i)[0]\n",
    "                    x_j_index = x_j_index[x_j_index != k]\n",
    "                    x_ = np.mean(X[x_j_index], axis=0)\n",
    "\n",
    "                    # 12.48\n",
    "                    mu_c = (s_i_freq*x_ + self.beta*self.mu_0) / (s_i_freq + self.beta)\n",
    "\n",
    "                    # 12.45\n",
    "                    _second_term = np.zeros((self.d, self.d))\n",
    "                    for x_index in x_j_index:\n",
    "                        _second_term += np.tensordot(X[x_index] - x_, X[x_index] - x_, axes=0)\n",
    "\n",
    "                    S_q_ = inv(self.S) + _second_term + \\\n",
    "                           (((s_i_freq*self.beta)/(s_i_freq+self.beta)) * \\\n",
    "                           np.tensordot(x_ - self.mu_0, x_ - self.mu_0, axes=0))\n",
    "\n",
    "                    # 12.46\n",
    "                    S_r_ = S_q_ + (s_i_freq+self.beta)/(1+s_i_freq+self.beta)* \\\n",
    "                                  np.tensordot(x_k-mu_c, x_k-mu_c, axes=0)\n",
    "\n",
    "                    # 12.44\n",
    "                    collapsed_term = pow(((s_i_freq+self.beta)/((1+s_i_freq+self.beta)*pi)), self.d/2)\n",
    "                    collapsed_term *= (pow(det(inv(S_r_)), ((self.v+s_i_freq+1)/2))*gamma((self.v+s_i_freq+1)/2))\n",
    "                    collapsed_term /= (pow(det(inv(S_q_)), ((self.v+s_i_freq)  /2))*gamma((self.v+s_i_freq+1-self.d)/2))\n",
    "\n",
    "                    # 12.13\n",
    "                    p = s_i_freq/(self.n-1+self.alpha) * collapsed_term\n",
    "                    pro.append((s_i, p))\n",
    "\n",
    "                # cal new class\n",
    "                S_b_ = inv(self.S) + (self.beta / (1 + self.beta) * np.tensordot(x_k - self.mu_0, x_k - self.mu_0, axes=0))\n",
    "\n",
    "                p = (self.alpha/(self.n-1+self.alpha)) * \\\n",
    "                    pow((self.beta/(1+self.beta)*pi), self.d/2) * \\\n",
    "                    (pow(det(inv(S_b_)), (self.v+1)/2 ) * gamma((self.v+1       ) / 2)) / \\\n",
    "                    (pow(det(self.S),    self.v    /2)   * gamma((self.v+1-self.d) / 2))\n",
    "\n",
    "                pro.append((max(self.class_freq.keys()) + 1, p))\n",
    "                s_k = self._sample(pro)\n",
    "                self.class_freq[s_k] += 1\n",
    "                self.class_list[k] = s_k\n",
    "\n",
    "            # cal posterior\n",
    "            ## 12.24\n",
    "            ### 11.11\n",
    "            V = np.log2(pow(self.alpha, len(self.class_freq))/self._af())\n",
    "            for s_i, freq in self.class_freq.items():\n",
    "                V += np.log2(factorial(freq-1))\n",
    "                V += np.log2(pow(self.beta/((freq+self.beta)*pow(pi, freq)), self.d/2))\n",
    "        \n",
    "                x_j_index = np.where(self.class_list == s_i)[0]\n",
    "                # 12.38 S_q\n",
    "                ## 12.39\n",
    "                x_ = np.mean(X[x_j_index], axis=0)\n",
    "                S_q_ = inv(self.S) + ((freq*self.beta)/(freq+self.beta))* \\\n",
    "                                     np.tensordot(x_-self.mu_0, x_-self.mu_0, axes=0)\n",
    "\n",
    "                for index in x_j_index:\n",
    "                    S_q_ += np.tensordot(X[index] - x_, X[index] - x_, axes=0)\n",
    "\n",
    "\n",
    "                numerator = pow(det(inv(S_q_)), (self.v+freq)/2)\n",
    "                denominator = pow(det(self.S), self.v/2)\n",
    "\n",
    "                for j in range(1, self.d+1):\n",
    "                    numerator *= gamma((self.v + freq + 1 - j) / 2)\n",
    "                    denominator *= gamma((self.v + 1 - j) / 2)\n",
    "\n",
    "                V += np.log2(numerator)\n",
    "                V -= np.log2(denominator)\n",
    "\n",
    "            if V > self.p_max:\n",
    "                self.drop_count = 0\n",
    "                self.p_max = V\n",
    "                print(i_ite, \"update\")\n",
    "            else:\n",
    "                self.drop_count += 1\n",
    "                self.class_list = previous_class_list\n",
    "                self.class_freq = previous_class_freq\n",
    "                if self.drop_count == self.stop_count:\n",
    "                    return self.class_list\n",
    "\n",
    "        return self.class_list\n",
    "    def _sample(self, pro:list):\n",
    "        p = np.array([pair[1] for pair in pro])\n",
    "        rand = np.random.random() * sum(p)\n",
    "        add_p = 0.\n",
    "        for index, p_value in enumerate(p):\n",
    "            add_p += p_value\n",
    "            if rand < add_p:\n",
    "                return pro[index][0]\n",
    "\n",
    "    def _af(self):\n",
    "        result = 1.\n",
    "        for i in range(self.n):\n",
    "            result *= (self.alpha + i)\n",
    "        return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 update\n",
      "1 update\n",
      "2 update\n",
      "3 update\n",
      "6 update\n",
      "7 update\n",
      "13 update\n",
      "[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 7 7 7 7 7 7 7 7 7 7 3 3 3 3 3 3 3\n",
      " 3 3 3 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "model = DirichletProcessMixtureGaussianDistribution(X)\n",
    "s = model.fit(ite=100)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEWCAYAAAAHC8LZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4FPW9P/D37OSeAEnIJlxCEoVEwi3cwTvCAbyAR0pt\nzRGLFqpcPK2cAt6aFgsWJFKpPy7SB6iceqTHSmyPp8dqa5HHCs0FgRBizEJIuJpbs0j2lt2Z+f2x\nEAk7STZhd2Yv79fz8Dyys9n9ZFzyzsx8vp8RzGazAiIioiBj0LsAIiKi3mCAERFRUGKAERFRUGKA\nERFRUGKAERFRUGKAERFRUGKAERFRUPJLgB08eBD5+fkYMWIEkpKSsHfv3g7bly1bhqSkpA5/Zs2a\n5Y9SiIgoREX440UtFgtGjhyJ/Px8LF26VPU599xzD379619DUdzrqCMjI/1RChERhSi/BNjMmTMx\nc+ZMAO6jLTVRUVFISUnxx9sTEVEY0O0a2D/+8Q9kZ2dj4sSJ+NGPfoSmpia9SiEioiDklyOw7syc\nORMPPvggMjMzcebMGaxduxYPPvggDhw4wFOJRETkFV0CbN68ee3/nZubi7y8PIwePRoffvgh5syZ\no0dJREQUZAKijX7AgAEYNGgQampq9C6FiIiCREAEWFNTEy5evIi0tDS9SyEioiDhtzb6mpoaKIoC\nWZZx7tw5HD9+vH3N14YNG/Dggw8iLS0NdXV1WLt2LVJTU3n60AdMJhOys7P1LiMo+GtfCefOQSwu\nBuLiOm5wuaD07w/pttt8/p7+xs+V97ivtOOXI7AjR47grrvuwrRp02C327F+/XrcfffdWL9+PURR\nRGVlJR599FFMmjQJy5cvR05ODj766CPEx8f7oxwi7cgyxM8/9wwvAHC5IE2YoH1NRCHKL0dgd9xx\nB1paWjrdvm/fPn+8LZHuhBMnAFn23OBwQMrJAaKjtS+KKEQFxDUwopDgcECsrlYPKVGEkpurfU1E\nIYwBRuQjYlmZenhZrZDGjQMM/OdG5Ev8F0XkCy0tEC5eBESx4+OyDKVfPyjp6frURRTCGGBEPhBR\nUqLeuGG1Qpo8WfuCiMIAA4zoBgmnTwMWCyAIHTe4XFAyMoC+ffUpjCjEMcCIboTLBfHYMSA2VnWb\nNG6c9jURhQkGGNENEMvLPY+8AMBuhzR8OBAVpX1RRGGCAUbUW1YrhJoa9ZCKioJyyy3a10QURhhg\nRL0klpSonzpk2zyRJvgvjKgXhPp6GJqaPENKlqEkJkIZNEifwojCCAOMqKcUBWJZGRS12Z1smyfS\nDAOMqIeE6mqgrc1zg9MJJSsL6NNH85qIwhEDjKgnnE6IlZVATIznNkmCNHas9jURhSkGGFEPiJ9/\nDkSo3MTBboc0YgQQGal9UURhigFG5K1LlyCcOaMeYFFRUHJytK+JKIwxwIi8JJaWqs47FCwW940q\n1RY0E5HfMMCIvCCcPQvh0iXVtnk5ORnKgAH6FEYUxhhgRN2RZYhHjqhPm7fZ2DZPpBOVk/lEdC3h\nxAlAllHb0IiNB4rRqIgwChJW3zEBmXl5QEKC3iUShSUGGFFX7HaIJhNqL7diwZ+LIX1rFcToOJxz\nWLFg30bsuftuZOpdI1GY4ilEoi6IZWVAVBQ2HvgmvABAjI6DNH81Ct/8nc4VEoUvBhhRZ1paIFy8\nCIgiGhWxPbyuEqPj0GRz6VQcETHAiDoRUVwMXJl3aBQkSA5rh+2Sw4qUWJ6FJ9IL//URqaj7+2fY\ntPcPaDREwShI+LfR2SgvKmw/jSg5rBA/3IpVBU/rXSpR2GKAEV2n7vRpLNzyX5DmP9fesFFeVIgN\nk7Px9vvr0dA/DSkJ0VhV8DQyM9nCQaQXBhjRdV7dtBXS/NUdGza+tQpv/+kVbPmPpZA5sJcoIPAa\nGNG1bDY0fW1XbdhoVETIo0frVBgRXY8BRnQNsaQExghZvWEjuQ8gijpVRkTXY4ARXSE0NMDQ1ITV\nd0+FWFTYHmKSwwqxqBArn1mic4VEdC2/BNjBgweRn5+PESNGICkpCXv37vV4zvr165Gbm4uBAwdi\nzpw5qKqq8kcpRN5RFIhlZVDi4pCVasRb907BuA8KMfC9dRj3/nrsef4pNmwQBRi/BJjFYsHIkSOx\nYcMGxKkMQN28eTO2b9+OwsJC7N+/H0ajEfPmzYPFYvFHOUTdEkwmwOFo/3tWqhHbHp6D38+fha1P\nPobMMWN0rI6I1PglwGbOnImf/OQnePDBByGo3CPpjTfewIoVKzBnzhwMHz4c27dvR2trK959911/\nlEPUNacTYmUlEBPjua2tDdLEidrXRETd0vwaWG1tLerr63HPPfe0PxYTE4PbbrsNxcXFWpdDBPHY\nMc/7fAGAwwEpO1s92IhId5oHWENDAwRBgNFo7PC40WhEQ0OD1uVQuLt8GUJtLRAZ6bnNYIAycqTm\nJRGRd7iQmcJGXV0dCne+hSabCymxEfj2v9yJ4WfOALGxnk+2WiFNmqR+ZEZEAUHzAEtNTYWiKGhs\nbMTgwYPbH29sbERqamqXX2symfxdXkjgfvJ0/sJFvLCrCOLcZyBGx+G8w4rizS8hW7HBHpeIZMWB\npeNHYEhyMqAoUCIi8LXDAXBftuPnynvcV97Jzs6+oa/XPMCysrKQlpaG/fv3Y+yVkTx2ux2HDh3C\nunXruvzaG/1mw4HJZOJ+UvGr/3ynPbwAQLrcAnNUP5z99s/cUzYcVvy4qBBv3ZuGrLhYuGbORFpi\nos5VBw5+rrzHfaUdv7XRHz9+HOXl5ZBlGefOncPx48dx7tw5AMDSpUuxefNmvP/++6isrMSyZcuQ\nkJCA+fPn+6McIjTZXB3GQ/3z0yIM+PYzHvMON35yCEp6OsDwIgp4fjkCO3LkCObOndveQr9+/Xqs\nX78e+fn52Lp1K370ox/Bbrdj9erVMJvNmDBhAoqKihB/5d5LRL6WEhuB8w5re2Apiqw+71ASIU2Y\noEeJRNRDfgmwO+64Ay0tLV0+59lnn8Wzzz7rj7cn8rBq8QIsXLsF0uzlEKPjoMjueYfXhpjksCKl\nXxwQFaVjpUTkLbZYUVjIzMzEnoKnkVf2Gwx89yXc3nYBrt+9fN28w41Y+eNlOldKRN5iGz2FjczM\nTLzxL3cATicgCCiuqMCeDwrRqIgwSg78+D8WIzMrS+8yichLDDAKG8KZM8DXXwNX5nMOSU7GtodH\nAbIMJTrave6LiIIGTyFSeJBliEePtodXB1YrpClTtK+JiG4IA4zCgnDiBCDLnhucTihZWUCfPprX\nREQ3hqcQKfTZ7RCrq9WPvmQZNcnJKCx4uX3E1KrFC3jvL6IgwCMwCnliWRkQHe3xuKGtDacSk7Bw\n/Rs4NuZRfHXXUhwb8ygWrt2Curo6HSolop5ggFFo++c/IXz1FSCKHpvkiAgU/vXT9rVhwJWJHLOX\no3DnW1pXSkQ9xACj0KUoiCgpUT11KFgssIwahSa7pDqRo8nm0qpKIuolBhiFLKG2FrBagevvCi7L\nkFNS4EpJQUpsRPti5qskhxUpsbw8TBToGGAUmlwuiOXl6vf6stna13ytWrwA4odbO07k+HArVi1e\noGW1RNQL/DWTQpLh+HH1DW1tUG6+GbgyOPrqiKlrb3S5quBpdiESBQEGGIUeqxWGmhr1tnkA0pgx\nHf6emZmJLWtf1KIyIvIhnkKkkCOWlgIxMZ4bbDZIo0YBEfy9jSgUMMAopAgNDTA0NgKG6z7aigLE\nxrpPHxJRSGCAUehQFIhlZVBUbowqWCxwTZ7s2ZFIREGLAUYhQzCZAIfDc4MkQR4wAOjfX/uiiMhv\nGGAUGpxOiJWV6te+HA7eKoUoBDHAKCSIR496XvcCAIcD8rBh6sFGREGNAUbB7/Jl99SNyEjPbQYD\n5FGjNC+JiPyPAUZBT+xk3iFsNkh5eaqDfIko+DHAKKgJFy5AaGlRb5tPSIDCiRpEIYsBRsFLliF+\n/nn7WKgOrFZ32zwRhSwGGAUtoaoKcKnc9sTlgjJoEJCYqH1RRKQZBhgFp7Y2iF9+qXqnZTidkCZM\n0L4mItIUA4yCklhWpj7T0OGAlJ2tHmxEFFIYYBR8zGYIFy6oB5jBAGXkSO1roqBnNpvxwo+fhNls\n1rsU8hIDjIJORGmpetu81Qpp/Hj1Bc1EXTCbzVizagEem1iNNasWMMSCBP+lU1ARzpwBvv7acyiv\nLEPp1w9Kero+hVHQuhpeq+bakZEWg1Vz7QyxIKFLgG3YsAFJSUkd/gwfPlyPUiiYSJJ7ZFRnR19s\nm6ceuja8EhPcp6QTEyIYYkFCtyOwnJwcmEwmVFdXo7q6GgcPHtSrFAoShhMn3AuUr+dyQcnIAPr2\n1b4oCmob167GU9Nb28PrqsSECDw1vRUb167WqTLyhm4BJooiUlJSYDQaYTQakZycrFcpFAzsdhhM\nJiAqynOby+W+9kXUQ6sLNmLH3xJgbu24ntDc6sKOvyVgdcFGnSojb+gWYHV1dcjNzUVeXh4WLVqE\n2tpavUqhICCWlqq3xtvtkEaMUB/kS9SNxMRErCl8C4Xvx7SHmLnVhcL3Y7Cm8C0kcjF8QNMlwCZN\nmoRt27Zh3759eP3111FfX4/Zs2fzfDOpa26GUF+vPpQ3KgpKTo72NVHIuDbEztTbGV5BRDCbzSoX\nFbRltVqRl5eHFStWYNmyZZ0+z2QyaVgVBQRFQd9PP4WgKB6dhwabDa3jx8NpNOpUHAWjy5cv4zc7\nXsUTT61Enz59un2c/Cc7O/uGvl5lJaj24uLiMHz4cNTU1HT5vBv9ZsOByWQKqf0knDoFMSkJiI3t\nuEGWocTGIuW223r92qG2r/wpVPaV2WzG1o0rsGx6K3Zs/6nHkdb4HW/36LU2rl2N1QUbO7xGqOyr\nYBAQ68DsdjtMJhPS0tL0LoUCicsFsaLCM7wA972+2DZPPdDZeq+6ujr8+OnHsfKHT3h9GYMLnwOD\nLgFWUFCAzz77DHV1dSgrK8PChQthtVqRn5+vRzkUoMTycvUNTieUrCwgIUHTeih4dbXe65lFs9BS\n9ykW3XoaBSvyuw0jLnwOHLoE2IULF/CDH/wAkydPxsKFCxETE4O//vWvSOcUBbrKYoFQU6PeNq8o\nkMaO1b4mClpdrff6xaKB6JcQgYy0GDw3z9lliF0bXoIAvPSbOggCGGI60SXAdu3ahcrKStTX1+PE\niRPYs2cPcthJRtcQS0s7P3U4YoT6IF+iTnS13mtL0QWs+rchANyB9tw8J37yjHqIXQ1CQQA2vHUW\nT9yfhg1vnYUggAufdRAQ18CIriXU18PQ1KQ+lDc2FsqwYdoXRUGts/Vea9+sw4vfy0C/+G9+IUpM\niMDy2TbVMFpdsBGv/zkGL//nGTz76BBkpMXg2UeH4OX/PIPX/xzDhc8aY4BRYFEUiGVlUOLjPTYJ\nFgukiRM9B/kSeeH69V5Pv34Wkqx4TCczt7qw/vdxnYZRhEHAC49ldLiW9sJjGYgw8HOpNQYYBRTB\nZALa2jw3yDJkoxEK13zRDbgaYr8ty8HmXR8hot8wrH2zrsNR2XNvinjgu/+tupB549rVWDbLpnot\nbdks9aM28h8GGAUOpxPiiRNATIznNrudbfPkE4mJifjFpl8jMzMTz/18B863AGvfrMOZejvW7K6F\nrBgQF6d+NMXZiYGFAUYBQzx6VL05o60N8tCh6k0dRL1kNpvx2rrleGNFJn76RCZ+/T8XIRoE/OQ7\nNnzwe/UmDs5ODCwMMAoMly9DqKvrtLtQHj1a44IolF2/LqxffATW/eAmFDyeie1/uIh/v++yR1t8\nba37qIyzEwMHA4wCglhc3Hnb/Jgx6oN8iXqpq3VhSx8aiN/8X32HtvjaWgG33trHI8R+W5bD8NIR\nA4x0J5w7B+HSJc+2eUUB4uOh3HSTPoVR0DKbzXjhx092urC4q2tZ2/9wEU/cn9bhmlZWloJDhy4j\nK+ublsWr19IYXvphgJG+ZBnikSNAXJznNosFLjZuUA95M6ews2tZr/zXWSx9aCB2/C3B48jq2vCi\nwMAAI10JX3wBSJLnBkmCMnAgkJSkfVEUtHoyp/D6a1n//v/OYsEso2p4UWBigJF+HA6IX36pfqfl\ntjb3omUiL3U1sLe7EPttWQ5e2/kR/lA5guEVRBhgpBvx8GEgMtJzg8MBKTtbfT0YUSe6aszoak7h\ntevCeE0ruDDASB9mM4QLF9Tb5g0GKCNHal8TBTV/LTLuriGE9MMAI11ElJSoN25Yre5bpagN8iXq\ngj8WGfPGlYGNPyVIc0JdHdDa6jmUV1GAvn2hZGToUxgFPV8uMuaNKwMfA4y0JUkQjx1TX7RsscA1\naZL2NVFI8cUi4940hJD2GGCkKUNFBSDLnhtcLijp6QAvoJMP3Ogi4942hJC2GGCkHbsdhlOn1Nvm\nnU5IEyZoXxORCk6dDw4MMNKMWFoKREV5brDbIQ0frr6NgpZ49BBgudzxQctl9+MBjlPngwMDjLTR\n3AxDfb36UN7ISCjDh2tfE/mVlD0KUft2QfzH39xBZrmMqH27IGWPCoog49T5wMcAI/9TFESUlEDp\nrG1+/Hi2zQc51aMtAFLOGIgnyhD91uuI+t12tM1fBADfBFmA49T5wMafGuR3Qk0NYLd7ts3LMpTE\nRCiDBulTGPmMlD0KMTteBhovuh+4erQ1NBfysFGIPPgXCC4nBGsrYna8jLb7vgvE9/nmBQL4iIxT\n5wMXA4z8y+WCWFGhPhbKaoXEafOhIb4P7I/9CLGvvQDhzEnE7HgZ0sAMRBXthuFkBSyv7oXiciJ+\nZT7kJCOi/ue33xyxWS4j6nfbAbtN3++Bgg4DjPxKLC9X3+B0uu/z1aeP+nYKOuL5WtgXrkB8wWK4\nbhqOqA/fhXjyBJToWCg2C8QvyyH1TQKslwG7FdFbX4Jw5iSit74EAJBGcw0g9QwDjPzHYnGfPlTr\nLpRlSHl52tdEPqF6zavhAmLeWAvrf2xAzB/ehBLfB5BlRJQXI27dvwOR0ZCG56FhyGCIZ04CrV8j\nvmAxlPgEtD2ytOMpRSIvMMDIb8TSUvWJG3Y7pBEj1CfRU1C42mHYHmKNFxH1tz9CyroF0Xu34fSS\nFTDUnwNkBWL9ORgcNkgZQ1E3bTou/2k3WpOSYTA3AQCECH4OqHcYYOQXQn09DM3N6t2F0dFQsrO1\nL4p8J74P2uYvQtS+XRAaLyJ2yxrYnlgJJPTF5bQBuOmN19CamQXYWgEA54feDPFUJbI3/hQZxlvQ\n74tySLeMgWXtTuBSC6Leet3dbk/UAwww8j1FgXj4sGrbvGCxuCduXN+RSMEnvg+c930X8SvzYV+0\nCjF/ehttM7+F+Lh+AICE06eByCjYH34KxrShgOgeyxTb8BWkBPdzIv9SBMcjSyCe+gKGkyd0+1Yo\nODHAyOeE6mrA4fDcIMuQU1KgpKVpXxT5nuUyIj/4b1he3YvIT/4X9m8vRuy2nwMuJ1w5o6FER0Ma\nNhKG8zWAywlERMA1ZBgESQIioyDWVkOwWxH1wX9fabfnPeCoZ3QNsJ07dyIvLw8DBgzAtGnTcOhQ\nYK4DoR5wOiF+8YV627zNxrb5UHFlnVfb/EVQjAPRNn8Ron/3BuSBGYgs+QSuO+6F7afbgIS+UKJi\nIJ6tgfPuuUDfRFhXb4LY0gjE90VkyScw/LMRbQt+CGnqdL2/KwoyugVYUVERnn/+eaxcuRKffvop\nJk+ejIcffhjnz5/XqyTyAfHIEfVxUW1tUG6+Wf0mlhR0RFOFe6rG1c7B+D5wTbkHaLPD8upeGOpM\nQNyV7sJ+ybAv+ylifr8DjkeWILL4b7C88DrQ4m7ikJONOn4nFMx0C7Bt27ZhwYIFeOyxx5CdnY2N\nGzciLS0Nu3fv1qskulFff+2+WWVEhOpmacwYjQsif5HG3uo5SaPOBMeyn7UfkUXt2wUAcM7+NiI/\n+V9Y1u5EzBsvo23aHET+/c+QRoxD220zAcC9kFllFBVRV3QJMKfTiaNHj2LatGkdHp8+fTqKi4v1\nKIl8QCwpUT/CstkgjRrVabBR8FM7Imubvwji8dL2U42GfzbCtuIXiPndG4DLibYFP3SfOsy71f0a\nx0t1/A4oGOnyE6W5uRmSJCE1NbXD40ajEQcOHNCjJLpBwrlzEC5d8gwwRQFiY92nDylkSWNv9Xww\nvg8QE9sebFef47zr/m+2A5CmToc0ehJEU4VW5VKICKpfiU0mk94lBAXN95MsI/HAASiiCFx3q3WD\n1Yqvp0yBdPKktjV5iZ8p7/VqX8WnABe+AvDVN4/1H3L1BT2fGyL/P/i58k72Da4H1SXA+vfvD1EU\n0dDQ0OHxxsZGj6Oya93oNxsOTCaT5vtJqKiAmJrqeadlSYLSrx9SArTzUI99Fay4r7zHfaUdXa6B\nRUZGYuzYsfjkk086PL5//35MnTpVj5KotxwOiNXVnuF1ZZs0iQNaicg/dDuFuHz5cixZsgTjxo3D\n1KlTsWvXLtTX1+Pxxx/XqyTqBfHwYfWZhg4H5Oxs9fVgREQ+oFuAzZs3Dy0tLdi0aRPq6+uRm5uL\n3//+90hPT9erJOopsxnChQtAfLznNoMB8khOViAi/9G1ieP73/8+vv/97+tZAt2AiK7a5idMUF/Q\nTETkI5yFSL0i1NUBra2eQ3kVBUhIgJKZqU9hRBQ2GGDUc5IE8dgx9Xt9Wa1wBWjXIRGFFgYY9Zih\nogKQZc8NLheUQYOAxETtiwpBpy9+oXcJRAGNAUY9Y7fDcPKketu8y+W+9kU+cfqrqp5/DUOPwggD\njHpELC3tfM1XTo76NupWd8HjbTD1JvSIghUDjLzX3AyhoUG9u1AUoeTmal9TiOgueBhMRJ4YYOQd\nRXG3zXfSuCGNHw8Y+HFSo3b0xFN9RDeOP3HIK0JNDWCzebbNyzKUfv2gDB6sT2FBQO3oqadHVC2t\njTdeB0OTQgwDjLrnckGsqOj86GvKFO1rCmEOWcah1lY0JqTjf8xm/K/ZjBrzVzhiscBxpfuzN2F0\n5DxvV0KhhQFG3RLLy9U3uFzuBct9+qhvp14529aGM21tSDBmI1IQIAoCFACn29rwucUCoPMjuJsG\nDFd9vNbhwK7GJtQ6HP4qm0hzQXU/MNKBxeI+fag271CSII0bp31NIeja4BkWE4N+oojTDgdssgwZ\nQEJ8f9wSE4Ps6OguQ+imgeqNNFnR0VhkTEEWu0QphDDAqEtiaan6qUO7HdKIEeqT6KnHrg8eY2Qk\njNfsWyVpEEbGxqLW4cCtX3yBrZKrx++RKPKfO4UWnkKkTgkNDTA0Nal3F0ZFQcnJ0b6oIKR2Wq+z\nU33dyYqOxqHcXIYRERhg1BlFgXj4MBSVU4eCxQJp4kTPjkRSpXZar7NTfZ2+xjWB19vTgL0NTaJA\nxQAjVYLJBNjtnhtkGXJKCpS0NO2LCmPXB15vwqinoUkU6Bhg5MnphFhZqX43ZZsNEqfN645hRMQA\nIxXi0aPq173a2qDcfLP6TSyJiDTGAKOOLl+GUFvbaXehNGaMtvUQAE7RIFLDAKMOxJIS9SMsmw3S\nqFFABLvf9MBhvkSeGGDUTrhwAYLZ7Hn6UFGAuDj36UMiogDBACM3WYb4+efqR19WK1yTJrFtnogC\nCgOMAABCVRXgUpnuIEnulvn+/bUvioioCwwwAtraIFZVdX6n5UmTtK+JiKgbDDCCePiwetdhWxvk\n7Gz19WCkKU7RIPLEAAt3ZjOEc+fUuwsFAfLIkdrXRB64cJnIEwMszEWUlanfKsVqhTR2LCCK2hdF\nROQFBlgYE86eBS5d8uwuVBSgb18oGRn6FEZE5AUGWLiSZYhHjnTdNk9EFMAYYGFKOHECkGXPDS4X\nlEGDgMRE7YsiIuoBXQLsgQceQFJSUvuf5ORkLF68WI9SwpPdDtFkUm+bdzrd9/oiIgpwugy2EwQB\nCxYswM9+9jMoigIAiGGrtmbEw4eBqCjPDQ4HpFtuUd9GRBRgdJvMGhsbi5SUFL3ePny1tEC4cAFI\nSPDcFhEBJZft2kQUHHS7BlZUVIShQ4fi1ltvRUFBAVpbW/UqJaxElJR03jY/frz6fcCIiAKQLkdg\n3/nOdzBkyBAMGDAAVVVVWLNmDSorK7Fv3z49ygkbwunTgMUCxMZ23CDLUBIT3c0bRERBQjCbzYov\nXmjdunXYtGlT528kCHj//fdx++23e2w7cuQIpk+fjgMHDmBMFzdMNJlMvig1PEkSEg8cgKIycUO0\n2WC+4w7IakdmRER+kp2dfUNf77MAa2lpQXNzc5fPSU9PV23WUBQFRqMRO3fuxEMPPeSLcsKWyWRS\n/VAYjh6Foa7Os0HD6YQyaFBYdh52tq/IE/eV97ivtOOzU4hXW+J7o6KiApIkIS0tzVfl0LVsNhhO\nnlS/9iXLkPLytK+JiOgGaX4NrLa2Fu+88w5mzZqF5ORkVFVVoaCgAGPHjsXUqVO1LicsiCUlnte9\nAMBuhzRypPokeiKiAKd5gEVGRuLAgQPYsWMHLBYLBg8ejNmzZ2P16tUQeMdfnxOammBobISidvQV\nHQ2FpzqIKEhpHmCDBw/Gn/70J63fNjwpCsTSUtXwEiwWuO66y3OQLxFRkOCinxAmnDwJ2O2eG2QZ\nstEIJTVV+6KIiHyEARaqXC6IlZXqd1O22yFx2jwRBTkGWIgSjx1TPz3Y1gb55pvVb6NCRBREGGCh\nqLXVPXWjk+5CefRojQsiIvI9BlgI6rRt3maDNGYMoDKNg4go2DDAQkxkYyMM//yn51BeRQHi4qBk\nZelSFxGRrzHAQomiIO7ECfU1X1YrXJMmsW2eiEIGAyyECF9+CcHl8twgSVAGDAD699e+KCIiP2GA\nhYq2NohffAGlszsth+GwXiIKbQywECEeOaLenOFwQM7OVl8PRkQUxBhgoeDrryGcOaMeYAYD5FGj\ntK+JiMjPGGAhQCwpUV+YbLVCGjvWsyORiCgE8CdbkBPOnYNw6ZJ623zfvlAyMvQpjIjIzxhgwUyW\n3de+1I6+LBZ32zwRUYhigAUxobISkCTPDS4XlPR0IDFR+6KIiDTCAAtWDgfE6mogOtpzm9MJacIE\n7WsiItJ05gcUAAAPmUlEQVQQAyxIiYcPqw7rFRwOSLm5gNp6MCKiEMIAC0ZmM4QLF1Tb5pWICCi3\n3KJDUURE2mKABaGILtrmrSNHsm2eiMICf9IFGaGuDmht9RzKK8tQEhPhTE3VpzAiIo0xwIKJJLnv\ntNzZvb6mTNG+JiIinTDAgoihogKQZc8NTqf7Pl8JCZrXRESkFwZYsLDbYTh1Sr1tXpYh5eVpXxMR\nkY4YYEFCLC1Vb4232yGNGKHaUk9EFMoYYMGguRmG+npAFD23RUdDyc7WviYiIp2p3H+DtFJXV4fC\nnW+hyeZCSmwEVi1egMzMzI5PUhRElJRAUWmbF6xWuO6807MjkYgoDPAITCd1dXVYuHYLjo15FF/d\ntRTHxjyKhWu3oK6ursPzhJoawGZTbZuXU1KgsG2eiMIUA0wnhTvfgjR7OcRo95GVGB0HafZyFO58\n65snuVwQKyrU2+btdkicNk9EYcznAbZnzx7MnTsXmZmZSEpKwtmzZz2eYzab8eSTTyIjIwMZGRl4\n6qmncOnSJV+XEtCabK728LpKjI5Dk831zd/Ly9W/uK0N8tCh6tM4iIjChM8DzGq1YsaMGXj++ech\ndHJtZvHixaioqMB7772HoqIilJeXY8mSJb4uJaClxEZAclg7PCY5rEiJvXJZ0mJxnz5U6zxUFMij\nR2tQJRFR4PJ5E8fSpUsBAEePHlXdXl1djY8//hgfffQRJly55cdrr72G++67D6dOncLQoUN9XVJA\nWrV4ARau3QJp9nJIl1vQ9Mk7EJrPwXLLEHx28BB+9+ZeNLUBRoOE1XdPQVaq0f2FNhukcePUOxKJ\niMKI5l2IJSUl6NOnDyZdc/1m6tSpiI+PR3FxcdgEWGZmJvYUPI2f/nIriuuakfZoAcToOJx0WPH4\n+peQ/MAPEDtoGM45rFhQVIi37p2CLGMKEB8P5aab9C6fiEh3mjdxNDQ0oH///h6Pp6SkoKGhQety\nfEqore3R8zMzM9Gnbz8Yr4QX4L4ONuB7P8Ol4g/a/y59axU2HigGLBa4Jk/2ddlEREHJqwBbt24d\nkpKSOv2TnJyMzz77zN+1BjShthZ9br21xyHWWTOHosgd/t4oG6AMHAgkJfmiXCKioOfVKcTly5fj\nkUce6fI56enpXr1hamoqmpubPR5vampCajdrmkwmk1fvoZeot99Gm9MJ9KDOaJcNksPaIcQkhxWC\nYOjw93jrJXyZmAjFi9cO9P0USLivvMd95T3uK+9k3+AUIa8C7OqRli9MnjwZra2tKC0tbb8OVlxc\nDKvViind3A7kRr9Zv+tFfT9fsbS9mUOMjoPksOKrPWuQPOdJAO7wEvdtRMHix5A5alS3r2cymQJ/\nPwUI7ivvcV95j/tKOz5v4mhoaEB9fT1MJhMURUFVVRXMZjOGDBmCxMRE5OTkYMaMGXjmmWewefNm\nKIqCFStW4N577w2bBo5rXW3mKNz+GzRfaIQxAvi3GWPxdsl/oVERYRQkrJ41Een/MkPvUomIAorP\nA2z37t145ZVXIAgCBEHAd7/7XQDA1q1bkZ+fDwDYuXMnVq9ejfnz5wMA7r//fmzcuNHXpQSNzMxM\nbJ95JwS7HTC4Tx3eMfwW90arFdLkyVAMHJpCRHQtnwfYc889h+eee67L5/Tr1w87duzw9VsHJG8G\n9goXLkBoaQHi4zt+saJA6dcPypAhGlZMRBQc+Gu9H3k1sFeWIX7+uWd4Ae6jr4kTtSuYiCiIMMD8\nyJuBvcIXXwAul+cXu1xQBg8GEhN7/f49beknIgomDDA/6nZgb1sbxC+/BKKjPb/Y5YJ0ZdRWb/R2\nXRoRUbBggPlRdwN7xcOHgchIzy90OCANH64+yNdLSlYWLh86BCUrq9evQUQUyBhgfrRq8QKIH25t\nDzHJYYX44VasWrwAMJshnDsHRKj00URGQrnllht+f4YXEYUyzYf5hpP2NV7XdiEWPI3MzExE/OUv\nqo0bgsUC1+23t7fTExGROgaYn2VmZmLL2hc7PCacOQN8/bXnDSllGXJysnvmIRERdYm/5mtNliEe\nPap+N2WbDRKnzRMReYUBpjHhxAlAlj03OJ3ua1YJCZrXREQUjHgKUUt2O0STCYiN9dymKJDy8nzy\nNt5M/yAiCnYMMA2JZWXqrfF2O6RRo9Rb6nvo/IWLWPPb99sXUJ93WLFw7RbsudI8QkQUKngKUSst\nLRAuXgRE0XNbTAyUYcN88ja/ee//up3+QUQUChhgGokoKem0bV6aNAkQBJ+8T0ub0vX0DyKiEMEA\n04Bw+jRgsXiGlCxDNhqhpKT47L2SooQup38QEYUKBpi/SRLE8nL1xg273edt80/Mu7/z6R9ERCGE\nv5b7meH4cfUNbW2Qhw5VD7YbMHjQwE6nfxARhRIGmD/ZbDCcOqW+aBmAPHq0X95WbfoHEVGo4SlE\nPxJLS4GYGM8NNpt7zZdaRyIREXmFAeYnQlMTDA0NnkN5FQVISOCkeCKiG8QA8wdFgVhaCkWlbR5W\nK1ycd0hEdMMYYH4g1NQAdrvnBpcLyqBBQGKi9kUREYUYBpivuVwQKyrUr325XJAmTtS+JiKiEMQA\n8zGxvFx9g8MBKSdHfRYiERH1GAPMl6xW99QNtZCKiICSm6t9TUREIYoB5kOdts1brZDGj/fsSCQi\nol7jT1QfERoaYGhq8gwpWYaSmOhu3iAiIp9hgPmCokA8fBiK2sQNmw3SlCna10REFOIYYD4gmEzq\nbfNOp3vBckKC5jUREYU6BtiNcrkgVlaqX/tSFPfIKCIi8jmfB9iePXswd+5cZGZmIikpCWfPnvV4\nzujRo5GUlNT+Jzk5GT//+c99XYomxGPH1JszbDZII0YAkZHaF0VEFAZ8Po3earVixowZeOCBB/DC\nCy+oPkcQBDz33HNYtGgRFEUBAMSrjV0KdK2t7rZ5tdpjY6EMG6Z9TUREYcLnAbZ06VIAwNGjR7t8\nXnx8PFJ8eCdiPYilpar38xIsFrimTfO8AzMREfmMbtfAtmzZgptvvhl33nknNm3aBKfTqVcpvSLU\n18PQ3KzaNi+npUEJ8nAmIgp0utzQcsmSJRgzZgySk5Nx+PBhrFmzBmfOnMGvfvUrPcrpOUWBWFam\nPm3e4YA0aZL2NRERhRnBbDYr3T1p3bp12LRpU+cvIgh4//33cfvtt7c/dvToUUyfPh3Hjh3DkCFD\nunz9P/7xj3jiiSdQU1ODRE5qJyIiL3h1BLZ8+XI88sgjXT4nPT2910WMHz8eiqKgpqYG48eP7/Xr\nEBFR+PAqwK62u/tLeXk5BEFAWlqa396DiIhCi8+vgTU0NKC+vh4mkwmKoqCqqgpmsxlDhgxBYmIi\nSktLUVpaijvvvBN9+/bF559/jhdffBH3338/Bg8e7OtyiIgoRHl1DawnNmzYgFdeeQXCdS3kW7du\nRX5+Po4dO4aVK1fCZDKhra0NQ4YMwfz58/HDH/4QMWrTLIiIiFT4PMCIiIi0ELCzEMNtJNWN8mZ/\nmc1mPPnkk8jIyEBGRgaeeuopXLp0SYdqA8sDDzzg8TlavHix3mUFjJ07dyIvLw8DBgzAtGnTcOjQ\nIb1LCjgbNmzo8BlKSkrC8OHD9S4rIBw8eBD5+fkYMWIEkpKSsHfvXo/nrF+/Hrm5uRg4cCDmzJmD\nqqoqr147YAPs6kiq559/3uN05FVXR1KZTCZUV1fjyy+/xMqVKzWuNDB4s78WL16MiooKvPfeeygq\nKkJ5eTmWLFmicaWBRxAELFiwoMPn6LXXXtO7rIBQVFSE559/HitXrsSnn36KyZMn4+GHH8b58+f1\nLi3g5OTktH+GqqurcfDgQb1LCggWiwUjR47Ehg0bEKdyy6nNmzdj+/btKCwsxP79+2E0GjFv3jxY\nLJZuX1uXhczeCKeRVL7Q3f6qrq7Gxx9/jI8++ggTJkwAALz22mu47777cOrUKQwdOlSzWgNRbGws\nP0cqtm3bhgULFuCxxx4DAGzcuBEff/wxdu/ejYKCAp2rCyyiKPIzpGLmzJmYOXMmAGDZsmUe2994\n4w2sWLECc+bMAQBs374d2dnZePfdd7Fw4cIuXztgj8C8FewjqbRSUlKCPn36YNI1U0KmTp2K+Ph4\nFBcX61hZYCgqKsLQoUNx6623oqCgAK2trXqXpDun04mjR49i2rRpHR6fPn06PzMq6urqkJubi7y8\nPCxatAi1tbV6lxTwamtrUV9fj3vuuaf9sZiYGNx2221efcYC9gjMG0E/kkpDDQ0N6N+/v8fjKSkp\naGho0KGiwPGd73wHQ4YMwYABA1BVVYU1a9agsrIS+/bt07s0XTU3N0OSJKSmpnZ43Gg04sCBAzpV\nFZgmTZqEbdu2ITs7G42NjSgsLMTs2bNRXFzM6UJdaGhogCAIMBqNHR43Go346quvuv16TQOsNyOp\nunLt4eiIESPQt29fPPHEE3jppZdC4kPj6/0VTnqy7773ve+1P56bm4usrCxMnz4d5eXlGDNmjBbl\nUpCbMWNGh79PmjQJeXl5ePvtt1VPm5FvaBpgHEnVM77cX6mpqWhubvZ4vKmpyeM37FBwI/tu7Nix\nEEURNTU1YR1g/fv3hyiKHkfojY2NIfmZ8aW4uDgMHz4cNTU1epcS0FJTU6EoChobGzsMsvD2M6Zp\ngHEkVc/4cn9NnjwZra2tKC0tbb8OVlxcDKvViilTpvjkPQLJjey7iooKSJIUMp+j3oqMjMTYsWPx\nySef4F//9V/bH9+/fz8eeughHSsLfHa7HSaTCXfddZfepQS0rKwspKWlYf/+/Rg7diwA9747dOgQ\n1q1b1+3XB+w1MI6k6pnu9ldOTg5mzJiBZ555Bps3b4aiKFixYgXuvffesO5ArK2txTvvvINZs2Yh\nOTkZVVVVKCgowNixYzF16lS9y9Pd8uXLsWTJEowbNw5Tp07Frl27UF9fj8cff1zv0gJKQUEB7r33\nXqSnp7dfA7NarcjPz9e7NN1ZLBbU1NRAURTIsoxz587h+PHjSEpKQnp6OpYuXYpf/vKXGDZsGIYO\nHYpXX30VCQkJmD9/frevHbCTODiSqme6218AcOnSJaxevRoffPABAOD+++/Hxo0b0bdvX83rDRTn\nz5/Hk08+iaqqKlgsFgwePBizZ8/G6tWrQ+I6qi/s3r0bv/rVr1BfX4/c3FysX7+e4X6dRYsW4dCh\nQ2hubkZKSgomTpyIF198ETk5OXqXpru///3vmDt3rsfPpvz8fGzduhUA8Morr+DNN9+E2WzGhAkT\n8Oqrr3q1EDxgA4yIiKgrQb8OjIiIwhMDjIiIghIDjIiIghIDjIiIghIDjIiIghIDjIiIghIDjIiI\nghIDjIiIghIDjIiIgtL/ByBIpM2862lQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ee44dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label = [\"o\",\"x\",\"D\",\"|\", \"_\"]\n",
    "color_iter = ['r', 'g', 'b', 'c', 'm']\n",
    "# print(s)\n",
    "splot = plt.subplot(1, 1,1)\n",
    "for i, c in enumerate(set(s)):\n",
    "    color = color_iter[i]\n",
    "    indexs = np.where(s == c)\n",
    "    x = X[indexs]\n",
    "    plt.plot(x[:,0],x[:,1], label[i])\n",
    "    covar = np.cov(x, rowvar=0)\n",
    "    v, w = linalg.eigh(covar)\n",
    "    u = w[0] / linalg.norm(w[0])\n",
    "    plt.scatter(x[0], x[1], .8, color=color)\n",
    "\n",
    "    mean = np.mean(x, axis=0)\n",
    "    angle = np.arctan(u[1] / u[0])\n",
    "    angle = 180 * angle / pi  # convert to degrees\n",
    "    ell = mpl.patches.Ellipse(mean, v[0], v[1], 180 + angle, color=color)\n",
    "    ell.set_clip_box(splot.bbox)\n",
    "    ell.set_alpha(0.3)\n",
    "    splot.add_artist(ell)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}