lecture2-bnp.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:f046c4b36ed613b6d103a67f7ad4c738a10bf0626e51044ecbee8eda3c5cc00c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bayesian Nonparametric Model\n",
      "============================\n",
      "\n",
      "Dirichlet process\n",
      "-----------------\n",
      "Definition:\n",
      "$$\n",
      "G_0 \\sim \\text{DP}(\\alpha, H)\n",
      "$$\n",
      "where $\\alpha$ is a concentration parameter, and $H$ is a base distribution."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Stick breaking process"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Stick breaking process is the constructive definition of Dirichlet process.\n",
      "\n",
      "$$\n",
      "\\text{for k=1 ... }\\infty \\\\\n",
      "V_k \\sim Beta(1, \\alpha)\\\\\n",
      "\\pi_k = V_k\\prod_{j=1}^{k-1}(1-V_j) \\\\\n",
      "\\theta_k \\sim H \\\\\n",
      "G_0 = \\sum_{k=1}^{\\infty} \\pi_k \\delta_{\\theta_k}\n",
      "$$\n",
      "\n",
      "Draw results of stick breaking process with truncation level T and concentration parameter alpha"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#stick breaking process, return stick lengths\n",
      "def sbp(_alpha, _T):\n",
      "    V = np.random.beta(1,_alpha,size=_T)\n",
      "    V[_T-1]=1\n",
      "    oneV = np.ones(_T)\n",
      "    oneV[1:] = 1-V[:_T-1]\n",
      "    return V*np.cumprod(oneV)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# concentration parameter alpha\n",
      "alpha_1 = 10\n",
      "alpha_2 = 1\n",
      "alpha_3 = 0.1\n",
      "# truncation level\n",
      "T = 50\n",
      "\n",
      "# draw \n",
      "p1 = sbp(alpha_1, T)\n",
      "p2 = sbp(alpha_2, T)\n",
      "p3 = sbp(alpha_3, T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#subplots, figsize = col_size x row_size\n",
      "plt.subplots(1,3, figsize=(18,4))\n",
      "\n",
      "plt.subplot(1,3,1)\n",
      "plt.bar(range(T), p1, )\n",
      "plt.title('alpha = %f' % alpha_1)\n",
      "\n",
      "plt.subplot(1,3,2)\n",
      "plt.bar(range(T), p2)\n",
      "plt.title('alpha = %f' % alpha_2)\n",
      "\n",
      "plt.subplot(1,3,3)\n",
      "plt.bar(range(T), p3)\n",
      "plt.title('alpha = %f' % alpha_3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.text.Text at 0x10df4b250>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAABBkAAAEKCAYAAABXBY1gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH/NJREFUeJzt3XuwrXddHvDnS0JQAhhoIGgIohIqFwNIG1PBchCLh1iB\n0ZliRJGLGKsRRm0F7FQObZWiZWQc2hg0UkRK6ghoaKMhKKdGjIHYQIAkQsTYJMQQrgaQNjHf/rHe\nQ1bW2Ze1z37P3muv9fnMrMla672s32/vfZ7s85z3Ut0dAAAAgO26x24PAAAAAFgOSgYAAABgFEoG\nAAAAYBRKBgAAAGAUSgYAAABgFEoGAAAAYBRKBpIkVfW8qrp07HUBOHKyGWB3yWHYOiUDe1JVPaaq\nLq6qW6vqzjWWP6Cq3l5Vn6+q66vqrE3295NVdXNVfa6qzq+q4+bdV1U9taquraovVNUfVdVDZ5a/\nuqo+OTz+48yyh1XVu4dtr6mqp84s//6q+uvhs99eVfffytcJYCdtls1rrP+4qvrzIQOvqKrHzixf\nyGwGWGQbZeca675+yMq/r6of2sq+5DDrUTKwV/2/JBckeeE6y/9zki8leVCS5yQ5t6oetdaKVfWd\nSV6a5NuTfG2Sr0/yynn2VVUnJnlrkn+T5P5Jrkjy36f2fXaSZyY5bXh89/DeIW9J8udJHjDs43eG\nfaaqHp3kV4fPPCnJF5P8l42/LAC7arNs/rLhF9XfS/KbSU5I8sYkv1dV9xyWL2Q2AyyyObJz1vuT\n/FiS/52kt7gvOczauttjRR5JXpbkuiR/m+TDSZ41tex5SS6den1nkp9I8pdJbk3yi0lqet0kv5Tk\n00k+lmT/1LbPT3L18Dl/meRHjuKcHp7kzpn3jk/yf5M8fOq9NyZ51Tr7+G9J/sPU66ckuXmefSX5\nkSR/MrXs3pmUAY8YXv9pkh+e+dpcNjx/RCbBfPzU8v+V5Ozh+S8k+a2pZV8/jOX4tebh4eGxNx+r\nks1rrPO0JDfOvPfXSZ42PF/IbPbw8Fi+xzLl8EbZucl2lyZ57rz7ksMeGz0cybBarkvypO6+XyYt\n5G9V1UkbrP+sJE9I8s2ZNI0vmFr2LUmuTfIPMgnX86eW3ZLku4bPeX6SX66qx6/1AVX1pKr6zAaP\nbz2CeT4iyR3dfd3Uex9I8uh11n/UsPyQq5KcNJyasNm+Hj29bXd/MZOv86Hla+17etuPdfcX5tz3\nxzIJ80esMw9gb1qVbJ716Ewycdp0Bi5qNgPLZ5lyeKPs3Co5zBFRMqyQ7v6d7v6b4flvJ/loJkG4\nnld392e7+4Ykr00yfZ7VX3f3+d3dmRzq+tVV9aBh3xd1918Nz/84yTuTfNs6Y/qT7r7/Bo8/PYKp\n3ieThnjabUnuu8H6n5t6fWjb+86xr7WW/+3M8tl932edZYf2fWj58Wssn943sARWKJtnrZWBm+Vn\nsjvZLHthiS1ZDm+UnVslhzkiSoYVUlXPraorDzWgSR6TScu6nhumnv+fJF8z9fpvDj0ZmslkCIaq\nenpV/VlVfWr4nDM3+ZyxfT7J/Wbe+6pMgm+e9b9q+O9tG+zrb6fW2eiz1tr35+cc5+enxjLPPIA9\naIWyedZa+XlCNs7PQ9vtdDafkMN/WQaWxJLl8EbZOea+5DDrUjKsiKr62iSvT/LjSR7Q3fdP8qEk\ntcFmD515ftMcn3OvTC7y8otJHjR8zkXrfU5VfVtV3bbB44lzTfDuPpLk2Kp6+NR7j81kvmv5cJLH\nzax7S3d/ZoN9fXhq2y9fDb2qjk/yDTPLZ/f9oallX19V95lZvt6+vyHJccOYgCWwYtk868OZXOxr\n2jdl4/xchGwGlsgS5vBG2blVcpgjomRYHcdncsXYTya5R1U9P5OWdiP/qqpOqKpTkrw4U1eE3cBx\nw+OTSe6sqqdncnGvNXX3pd193w0e71lv26r6iuGzUlX3GsI7w/lbb0vy76rq3lX1pCTfneRN6+zq\nN5O8sKoeOZxj9m+TvGHOfb09yWOq6nuG8bwiyfu7+yNT+/6pqvqaqjo5yU8l+a/Dvj+SyRV9X1FV\nX1FV35PJ9+Stw7ZvzuRKu08agvnfJ3nrzPlpwN62Mtm8hoNJ/r6qXjys9+JMLqj2R8PyRc1mYLks\nWw6vm51rqap7Djl5jyTHDblXm+1LDrMRJcOK6O6rk7wmyWWZHMb1mCR/Mr1KZm5bk8mtxf48yZVJ\n/kfuunDNWuv28Dm3ZRK2v53JVXXPGvYzqqp6WCZXqP3Q8Nl/l+SaqVV+LMlXJvlEkt9K8qPdfc2w\n7UOHBvghw5gvzqRVfneS6zO52u8r5tlXd38yyfcm+flhvv8oyfcd2rC7z0vyjiQfzOSCNu/o7tdP\n7fv7hm0+Pezje7v7U8O2Vyf50UzKhluGMfzY1r9awKJatWyuqouq6mXDmG7P5OJpz03ymeG/z+ru\nO4blC5nNwHJZthzeLDunc3hwSSa5fUYmR3R8McN1IuQwR+rQ7VbWX6FqfyYXNDkmya9396tnlj8n\nyc9kcqjPbUn+ZXdfNSy7PpNzZ/4+ye3dffrYE+DoqKo7M7klzcd2eyzA+qrqN5J8V5JPdPc3rbPO\nryR5eia/ODyvu6/cwSEyItkMi0cOrxY5DJvb8EiGqjomyeuS7M/kNiNnVdUjZ1b7WJJ/2t2nZXI4\n93QD1Un2dffjFQwAR8UbMsnoNVXVmZn8MnRqJvesPnenBgawIuQwwJTNTpc4Pcl13X39cFjjBZnc\nC/bLuvuy7j50i5HLkzxkZh8bXTSFxbXxIS7AQujuSzM51Hw9z0jyxmHdy5OcUBvf+5vFJpthwcjh\nlSOHYROblQwn5+63aLlxeG89L8zkKqmHdJJ3VdUVVfWiIxsiu6G7j3EYGCyFtXJ8tgxmj5DNsCfJ\n4SUih2Fzx26yfO6mrqqekuQFSaZvp/LE7r65qh6Y5JKqunZoewHYObNHlPlXGICdJYeBlbFZyXBT\nklOmXp+SSft6N1V1WpJfS7J/+h6s3X3z8N9bq+rtmZx+cenMtkIWWFjdvddP+ZrN8Ydkjft5y2Jg\nUclhgN23lSze7HSJK5KcWlUPq6rjkjw7yYXTK1TVQzO5R+oPdPd1U+/fu6ruOzw/PpP7wH5wnQGv\nxOMVr3jFro/BXM3VfOd/LIkLM7k1YKrqjCSf7e5b1lpxt7/efobN1VzNdfaxJObO4bvugHjXnRB3\n+3vgZ9hczdV8t2rDIxm6+46qOifJxZncwvL87r6mqs4elp+X5OeS3D/JuVWV3HWrygcnedvw3rFJ\n3tzd79zyCAFYV1W9JcmTk5xYVTdkcv/qeyaTjO7ui6rqzKq6LskXkjx/90YLsHzkMMDdbXa6RLr7\n95P8/sx75009/+EkP7zGdh9L8rgRxgjAOrr7rDnWOWcnxgKwiuQwwN1tdroEI9q3b99uD2HHmOvy\nWrX5snxW6WfYXJfTKs2V5bRKP8PmurxWbb5bUUdyjsWoA6jq3R4DwFqqKr33Lzg2F1kMLKJVy+HD\nbzpRR3Q+NMCYtprFjmQAAAAARqFkAAAAAEahZAAAAABGoWQAAAAARqFkAAAAAEahZAAAAABGoWQA\nAAAARqFkAAAAAEZx7G4PAIDFUFV3e93duzQSAAD2KiUDAIPpUqHWXQsAANbjdAkAAABgFAtxJMPs\nIbqJw3QBAABgr1mIkuHuh+gmDtMFAACAvcfpEgAAAMAolAwAAADAKJQMAAAAwCiUDAAAAMAolAwA\nAADAKJQMAAAAwCiUDAAAAMAolAwAAADAKJQMAAAAwCiUDAAAAMAolAwAAADAKJQMAAAAwCiUDAAA\nAMAolAwAAADAKJQMAAAAwCiUDAAAAMAolAwAAADAKJQMAAAAwCiUDAAAAMAolAwAAADAKJQMAAAA\nwCg2LRmqan9VXVtVH62ql66x/DlV9YGquqqq3lNVp827LQDbM0dGn1hVf1BV76+qD1XV83ZhmABL\nTRYD3KW6e/2FVcck+Ysk35HkpiTvS3JWd18ztc4/SXJ1d3+uqvYnOdDdZ8yz7bB9J7NjqGw0LoCd\nUFXp7trtcaxnzow+kORe3f3yqjpxWP+k7r5jZl8zWSyHgd236DmcjJfFficGFtVWs3izIxlOT3Jd\nd1/f3bcnuSDJM6dX6O7Luvtzw8vLkzxk3m0B2JZ5cvbmJPcbnt8vyadmCwYAtkUWA0zZrGQ4OckN\nU69vHN5bzwuTXHSE2wKwNfPk7K8leXRVfTzJB5K8ZIfGBrAqZDHAlGM3WT738VlV9ZQkL0jyxK1u\nmxyYer5v/s0ARnTw4MEcPHhwt4exFfPk7M8meX9376uqb0hySVU9trtvO3zVA+OODmCL9mAOJ6Nm\n8YGp5/vGGh/Almw3ize7JsMZmVxjYf/w+uVJ7uzuV8+sd1qStyXZ393XbXFb558BC2nRzwWeJ2er\n6qIkP9/d7xle/2GSl3b3FTP7ck0GYOEseg4n42Wx34mBRTX2NRmuSHJqVT2sqo5L8uwkF8584EMz\nKRh+4FDBMO+2AGzLPDl7bSYXI0tVnZTkHyb52I6OEmC5yWKAKRueLtHdd1TVOUkuTnJMkvO7+5qq\nOntYfl6Sn0ty/yTnVlWS3N7dp6+37VGcC8BKmTOjfyHJG6rqA5kUyz/T3Z/etUEDLBlZDHB3G54u\nsSMDcGgYsKD2wmG6Y3G6BLCIVjuHE1kMLIKxT5cAAAAAmIuSAQAAABiFkgEAAAAYhZIBAAAAGIWS\nAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIB\nAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEA\nAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAA\nABiFkgEAAAAYhZIBAAAAGIWSAQAAABiFkgEAAAAYhZIBAAAAGIWSAQAAABjFpiVDVe2vqmur6qNV\n9dI1ln9jVV1WVV+qqp+eWXZ9VV1VVVdW1XvHHDgAm2f0sM6+IYc/VFUHd3iIAEtNDgPcXXX3+gur\njknyF0m+I8lNSd6X5KzuvmZqnQcm+dokz0ryme5+zdSyv0ryhO7+9Aaf0cnsGCobjQtgJ1RVurt2\nexzrmTOjT0jyniTf2d03VtWJ3f3JNfY1k8VyGNh9q53DiSwGFsFWs3izIxlOT3Jdd1/f3bcnuSDJ\nM6dX6O5bu/uKJLevN6Z5BwPAlmya0Um+P8lbu/vGJFnrF1sAjpgcBpixWclwcpIbpl7fOLw3r07y\nrqq6oqpetNXBAbCheTL61CQPqKp3D1n8gzs2OoDlJ4cBZhy7yfLtHp/1xO6+eTil4pKqura7L93m\nPgGYmCej75nkm5M8Ncm9k1xWVX/W3R89qiMDWA1yGGDGZiXDTUlOmXp9SiYN7Vy6++bhv7dW1dsz\nOaRsjZLhwNTzffPuHmBUBw8ezMGDB3d7GFsxT0bfkOST3f13Sf6uqv44yWOTrPHL7YGjMkiAecnh\nA1PP9404TID5bTeLN7vw47GZXMzmqUk+nuS9mbmYzdS6B5LcdujCj1V17yTHdPdtVXV8kncmeWV3\nv3NmOxe5ARbSHrjg2KYZXVXfmOR1Sb4zyb2SXJ7k2d199cy+XPgRWDirncOJLAYWwVazeMMjGbr7\njqo6J8nFSY5Jcn53X1NVZw/Lz6uqB2dyJd37Jbmzql6S5FFJHpTkbVV16HPePFswAHDk5sno7r62\nqv4gyVVJ7kzya7O/2AJwZOQwwOE2PJJhRwagtQUW1KL/C9qYHMkALKLVzuFEFgOLYOxbWAIAAADM\nRckAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxC\nyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJ\nAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskA\nAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAAAACMQskAAAAAjELJAAAAAIxCyQAA\nAACMQskAAAAAjELJAAAAAIxi05KhqvZX1bVV9dGqeukay7+xqi6rqi9V1U9vZVsAtmfenK2qf1xV\nd1TV9+zk+ACWnRwGuLsNS4aqOibJ65LsT/KoJGdV1SNnVvtUkp9I8p+OYFsAjtC8OTus9+okf5Ck\ndnSQAEtMDgMcbrMjGU5Pcl13X9/dtye5IMkzp1fo7lu7+4okt291WwC2Zd6c/Ykkv5Pk1p0cHMAK\nkMMAMzYrGU5OcsPU6xuH9+axnW0B2NymOVtVJ2fyC++5w1u9M0MDWAlyGGDGZiXDdkJQgAIcXfPk\n7GuTvKy7O5NDdB2mCzAeOQww49hNlt+U5JSp16dk0tDOYwvbHph6vm/O3QOM6+DBgzl48OBuD2Mr\n5snZJyS5oKqS5MQkT6+q27v7wsN3d+CoDBJgXnL4wNTzfSMOE2B+283impSq6yysOjbJXyR5apKP\nJ3lvkrO6+5o11j2Q5Lbufs1Wtq2qPrwErmw0LoCdUFXp7oX9F6etZPSw/huSvKO737bGspkslsPA\n7lvtHE5kMbAItprFGx7J0N13VNU5SS5OckyS87v7mqo6e1h+XlU9OMn7ktwvyZ1V9ZIkj+ruz6+1\n7ZFNC4BZ82T0rg4QYMnJYYDDbXgkw44MQGsLLKhF/xe0MTmSAVhEq53DiSwGFsGoRzKwseHcusP4\nnwEAAACrSMmwbYc3zgAAALCKNruFJQAAAMBcHMnAYdY6DcQpIAAAAGxGycA67n4BOAAAANiM0yUA\nAACAUSgZAAAAgFEoGQAAAIBRKBkAAACAUSgZAAAAgFEoGQAAAIBRKBkAAACAUSgZAAAAgFEoGQAA\nAIBRKBkAAACAURy72wNYT1Wt+X537/BIAAAAgHksbMkwMVsorF08AAAAALvP6RIAAADAKJQMAAAA\nwCiUDAAAAMAoFvyaDIzBRTQBAADYCUqGleEimgAAABxdTpcAAAAARqFkAAAAAEahZAAAAABGoWQA\nAAAARqFkAAAAAEahZAAAAABGoWQAAAAARqFkAAAAAEahZAAAAABGoWQAAAAARqFkAAAAAEahZAAA\nAABGoWQAAAAARqFkAAAAAEaxaclQVfur6tqq+mhVvXSddX5lWP6Bqnr81PvXV9VVVXVlVb13zIED\nsHlGV9Vzhmy+qqreU1Wn7cY4AZaZLAa4y7EbLayqY5K8Lsl3JLkpyfuq6sLuvmZqnTOTPLy7T62q\nb0lybpIzhsWdZF93f/qojB5ghc2T0Uk+luSfdvfnqmp/ktfnrowGYJtkMcDdbXYkw+lJruvu67v7\n9iQXJHnmzDrPSPLGJOnuy5OcUFUnTS2vsQYLwN1smtHdfVl3f254eXmSh+zwGAGWnSwGmLJZyXBy\nkhumXt84vDfvOp3kXVV1RVW9aDsDBeAw82T0tBcmueiojghg9chigCkbni6RSUkwj/WOVnhSd3+8\nqh6Y5JKqura7Lz18tQNTz/fN+ZEA4zp48GAOHjy428PYinkzOlX1lCQvSPLE9dc6sO0BAWzHHszh\nZNQsPjD1fN82hgRw5LabxdW9fi5W1RlJDnT3/uH1y5Pc2d2vnlrnV5Mc7O4LhtfXJnlyd98ys69X\nJPl8d79m5v0+PJsPdRaHv7/ReHdaVWXRx5hsfZyHr794c4KdUFXp7oU95WuejB7ePy3J25Ls7+7r\n1tlX+3MPLJpFz+FkvCxe73diWQzstq1m8WanS1yR5NSqelhVHZfk2UkunFnnwiTPHT78jCSf7e5b\nqureVXXf4f3jkzwtyQfnHRgAm9o0o6vqoZn8UvsD6xUMAGyLLAaYsuHpEt19R1Wdk+TiJMckOb+7\nr6mqs4fl53X3RVV1ZlVdl+QLSZ4/bP7gJG+b/Kt4jk3y5u5+59GaCMCqmSejk/xckvsnOXfI49u7\n+/TdGjPAspHFAHe34ekSOzIAp0scdU6XgCOzFw7THYvTJYBFtNo5nMhiYBGMfboEAAAAwFyUDAAA\nAMAolAwAAADAKDa88CPsNcPFlA7jfEYAAICjT8nAElrvQqIAAAAcTU6XAAAAAEahZAAAAABGoWQA\nAAAARqFkAAAAAEahZAAAAABGoWQAAAAARrHnbmFZtfbtCLtnb1sIAAAA7KQ9VzJMzBYKaxcPG1mr\nrFBUAAAAwJHboyXDWKZLha0XFQAAAMBdlqZkcBoFAAAA7K6lKRkmtn8aBeNQ+gAAAKyeJSsZWCxK\nHwAAgFXiFpYAAADAKJQMAAAAwCiUDAAAAMAoXJMB1rDWhStdtBIAAGBjSgYWwmLejWL6s120EgAA\nYDNKhqNgMf/CvBe4GwUAAMBepmQ4ag7/C/N65cO6e1BKLBynUQAAAKxPybDj1vvX+vlLCX+p3W1O\nowAAAFiLkmHhLdcpBFs9mgMAAIC9Q8nALliu4gQAAIAJJQN70l4+IsJpMAAAwLJSMrCH7eUjIvby\n2AEAANZ2j90eAAAAALAclAwAAADAKJwuMae9fA0AAAAA2AlKhi2ZPo9+9UqHtYoWFysc70KOLggJ\nAADsdUoGtmi1i5b1jXUhRxeEBAAA9i4lAwBb5sgmAADWsumFH6tqf1VdW1UfraqXrrPOrwzLP1BV\nj9/KtmxdVa35WEarNNf1+Bqwke1k9Pb11ANgde1uFgMslg1Lhqo6JsnrkuxP8qgkZ1XVI2fWOTPJ\nw7v71CQ/kuTcebfdCcv7F7TOIvyCvzNf3+3Pde//HBz+NVhvTkd7rgcPHhxtX2zPdjJ6C5+xh//c\nrG2VfobNdTmt0lz3gp3I4mWzSj/D5rq8Vm2+W7HZkQynJ7muu6/v7tuTXJDkmTPrPCPJG5Okuy9P\nckJVPXjObXfI/H9JXcZfqI++xSg8NrdXxrkV683p6M11GQN1D/+5P9KMPmlrH3Pk+bmIX8tl/Ble\nj7kup1Wa6x6xQ1m8PFbpZ9hcl9eqzXcrNisZTk5yw9TrG4f35lnna+bYdkHtzb+I7pVf8GG3bPxn\nZE/+uT/SjH7I0RvSMpZ5ABtawCwG2D2blQzz/obob7ILwy/4W3E0S5lFK322Op711n3lK1+55VM1\ntvOZR/I1m79M2PN/Ro40o3d04jvzfQXYNXsiiwF2Sm10NfCqOiPJge7eP7x+eZI7u/vVU+v8apKD\n3X3B8PraJE9O8nWbbTu8L2CBhdXdC/u32O1kdHffMrMvWQwspEXO4WS8LJbDwCLbShZvdgvLK5Kc\nWlUPS/LxJM9OctbMOhcmOSfJBUPIfra7b6mqT82x7cL/jwNggR1xRs/uSBYDHLFRslgOA8tiw5Kh\nu++oqnOSXJzkmCTnd/c1VXX2sPy87r6oqs6squuSfCHJ8zfa9mhOBmCVbCejARiHLAa4uw1PlwAA\nAACY12YXfjyqqmp/VV1bVR+tqpfu5ljGVlW/UVW3VNUHp957QFVdUlUfqap3VtUJuznGsVTVKVX1\n7qr6cFV9qKpePLy/dPOtqq+oqsur6v1VdXVVvWp4f+nmekhVHVNVV1bVO4bXSznXqrq+qq4a5vre\n4b2lnOs0Obwc31M5vJxzPWRVcjiRxbJ4b5PFyznXQ1Yli8fI4V0rGarqmCSvS7I/yaOSnFVVj9yt\n8RwFb8hkbtNeluSS7n5Ekj8cXi+D25P8ZHc/OskZSX58+F4u3Xy7+0tJntLdj0tyWpKnVNWTsoRz\nnfKSJFfnrqtgL+tcO8m+7n58d58+vLesc00ih7Nc31M5vIRznbIqOZzIYlm8t8niJZzrlFXJ4m3n\n8G4eyXB6kuu6+/ruvj3JBUmeuYvjGVV3X5rkMzNvPyPJG4fnb0zyrB0d1FHS3X/T3e8fnn8+yTWZ\n3A96Wef7xeHpcZmce/mZLOlcq+ohSc5M8uu569ZbSznXwexFt5Z5rokcXprvqRyWw7swtKNJFsvi\nPUkWy+JdGNrRsq0c3s2S4eQkN0y9vnF4b5mdNHUl4VuSnLSbgzkaanJl5ccnuTxLOt+qukdVvT+T\nOb27uz+cJZ1rkl9O8q+T3Dn13rLOtZO8q6quqKoXDe8t61wPkcPL9z2VwxNLM9esVg4nsjiRxUtB\nFidZorlmtbJ42zm82S0sj6aVvuJkd3ct2f2Qq+o+Sd6a5CXdfVvVXQXYMs23u+9M8riq+qokF1fV\nU2aWL8Vcq+qfJ/lEd19ZVfvWWmdZ5jp4YnffXFUPTHJJTe5h/mVLNtdDlm0+W7KM31M5/OXlSzHX\nFczhRBavnGX8nsriLy9firmuYBZvO4d380iGm5KcMvX6lEya22V2S1U9OEmq6quTfGKXxzOaqrpn\nJmH6pu7+3eHtpZ1vknT355L8zyRPyHLO9VuTPKOq/irJW5J8e1W9Kcs513T3zcN/b03y9kwOX13K\nuU6Rw0v0PZXDSznXlcrhRBYPZPEeJouXcq4rlcVj5PBulgxXJDm1qh5WVccleXaSC3dxPDvhwiQ/\nNDz/oSS/u8G6e0ZN6tnzk1zd3a+dWrR0862qEw9dTbWqvjLJP0tyZZZwrt39s919Snd/XZLvS/JH\n3f2DWcK5VtW9q+q+w/PjkzwtyQezhHOdIYeX5Hsqh+Xwbo1xTLJYFu91slgW79YYxzJWDlf37h3V\nUVVPT/LaTC4Ucn53v2rXBjOyqnpLkicnOTGT81Z+LsnvJfntJA9Ncn2Sf9Hdn92tMY6lJleS/eMk\nV+WuQ/5enuS9WbL5VtU3ZXKxk3sMjzd19y9V1QOyZHOdVlVPTvLT3f2MZZxrVX1dJk1tMjmN7M3d\n/aplnOssObwc31M5LIeXYa6yWBbv1hjHIotl8V6f61g5vKslAwAAALA8dvN0CQAAAGCJKBkAAACA\nUSgZAAAAgFEoGQAAAIBRKBkAAACAUSgZAAAAgFEoGQAAAIBRKBkAAACAUfx/5/BWS3hQ02oAAAAA\nSUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10dbec3d0>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Normalized gamma process"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#normalized gamma process\n",
      "def ngp(_alpha, _T):\n",
      "    V = np.random.gamma(1,_alpha,size=_T)\n",
      "    return (V, V/np.sum(V))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha = 1\n",
      "T1 = 50\n",
      "T2 = 100\n",
      "T3 = 200\n",
      "\n",
      "v1,p1 = ngp(alpha, T1)\n",
      "v2,p2 = ngp(alpha, T2)\n",
      "v3,p3 = ngp(alpha, T3)\n",
      "\n",
      "#subplots, figsize = col_size x row_size\n",
      "plt.subplots(1,3, figsize=(18,4))\n",
      "\n",
      "plt.subplot(1,3,1)\n",
      "plt.bar(range(T1), v1, )\n",
      "plt.title('T = %f' % T1)\n",
      "\n",
      "plt.subplot(1,3,2)\n",
      "plt.bar(range(T2), v2)\n",
      "plt.title('T = %f' % T2)\n",
      "\n",
      "plt.subplot(1,3,3)\n",
      "plt.bar(range(T3), v3)\n",
      "plt.title('T = %f' % T3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<matplotlib.text.Text at 0x10f450610>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAABBYAAAEKCAYAAACmDtbtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XvQJXV95/HPZxiQmzIakonCuGCElFhGUJewSMIhIiqr\nuJuyslBRE7airBsjoWqNl9XwWLsVc9GoJFFJBIKWgU0wKiawXjmslBZemOE2g4oryYAwuHITRrOw\nfPeP7memn/OcS/c53adv71dVz5ynT58+v9/vdH9P9/f8ft2OCAEAAAAAAMxjQ90FAAAAAAAA7UVi\nAQAAAAAAzI3EAgAAAAAAmBuJBQAAAAAAMDcSCwAAAAAAYG4kFgAAAAAAwNxILAAAAAAAgLmRWIBs\nP2z7R+n0uO3dmb/PKvF9/tr2v2TW/ZBtZ54/1vY3bT9i+xu2nztlXU+wfbHtB23fbfu8keenrsv2\neenrHrR9ke39Ms89xfYn03a5Y7QNbL/I9m3pur9k++kjz/+R7f+TTn84b3sB6LYlxt5fs/2VNGZd\nM+b5ueNlmeuqMvbaPsL2Nelrd9h+Ub7WA9BlS4zD77H97fTYd4ft14w8TxxG65FYgCLi4Ih4YkQ8\nUdI/SXr56t8RcVmZbyXpjzLrflJEhCSlQe3Tkj4qaZOkSyV92va+E9a1IunnJD1d0imSfs/2S/Ks\nK13uLZJ+RdK/kvQMSe/KrPsvJP1E0s9I+nVJH7J9TPraQyV9QtJ/lfRkSd+Q9D9WX2j7HEmvlPQL\n6fSKdB4ArLHE2PtDSX8qaV2is4R4Wea6qoy9l0n6pqSnpOu4Il0ngB5bYhx+OF33kyT9hqQP2P43\nEnE4bwOiBSKCiWnPJOl7kn6lonVfIum/TXjuNEl3jsz7J0kvmbD8XZJOzfz9LkmXzVjXaenjv5H0\n3zPPnSLp7vTxQZL+RdIzM89fKund6ePXS7ou89yBknZLOjr9+yuSfivz/NmSvlr358rExNTsqcrY\nm3mP35J0zci8uePlmPU3MvZKOlrJgfJBmeevlXRO3Z87ExNTc6ZlxOHMe31a0nnpY+IwUycmeiyg\nMNtvtX3/hOm+GS//z7Z/mHbN+tXM/GdLumlk2RvT+aPv/2RJT02fX3VTZtlZ6zpmzGs3p+s9WtJj\nEXH7hNc+O/vaiNgt6fYZ615XBwAoasHYO8ki8bLMdVUZe58t6X9HxCMT1g0AuZQRh20fIOlfS7o1\nnUUcRieQWEBhEfGHEfHkCdNTprz0AknPlPTTkt4p6a9tn5g+d7CkB0eWf0jSE8es5+D0/wcnLDtr\nXaPPP5T+/8T0uYe01o9GXjv6/Kx1HywAWNACsXeaReJlmeuqMvYW+X4BgIlKisMflrQtIj6X/k0c\nRieQWMDSRMTWiLg/Ih6PiKslfVzSaq+FhyU9aeQlh2h9AFtdViPLH6Ik+E1a16Ypzx+S/v+jHOX4\n0YTnp637YQFAM42LaXnjZZnrqjL2TvpOGPf9AgCVsf0nSn7Z/7XMbOIwOoHEAgqz/fbMFXNHp3kD\nxK1KLvSS9Qva201sj4i4X9Ldko7NzH6upFumrOs5mXXdOua1u9L1flvSRtvPHHk++9o9V9e1fZCS\ni0hOW/ctAoAFlRB7Y8y8ReJlmeuqMvbeKukZtg8eeX7d9wsATLNIHLb9LkkvUXK9g+yPTsRhdEPd\nF3lgatakai/e+ColXaE2KLm4zEOSfjl9bj9Jd0h6k6QnpP9/T9LGCet6t6Shkmzns5QkGk7Lsy4l\nQf3u9HVPTtfzB5l1X6bk4jYHSjpJ0gOSnpU+d2j6969K2l/SH0v6Sua150jaLulpkg5TEjBfX/fn\nysTE1Oyp4ti7IY1X/0nJxbKeIGnf9Ll9F4mXI++z0LqqjL2SvirpT9LX/qqk+yX9VN2fOxMTU3Om\niuPw25ScuG8e8xxxmKkTU+0FYGrWVHFQ/V9pQHpQ0lZJvzby/LFKbl2zO/3/uZnnfl3SLZm/95N0\nUbqueyT9bt51pc+fl77uwXQ9+2aee7KkTyrptnWHpDNHXvsiSTvSdX9J0tNHnv8jJbd3+6GkP6z7\nM2ViYmr+VHHs/U1Jj49MF2eeXyReXiXprSWtq7LYq+S2atekr91RVVszMTG1d6o4Dj8u6cdKhg2s\nTmXFTuIwUyMmpx/0VLb3STfMOyPiFWOev0DSy9IN5TcjYuvMlQIASmf7bZJereQg5mZJZ0fEv9Rb\nKgDoB9s/L+nyzKxnSHpnRFxQU5EAYCnyXmPhXCVdW9ZlIWyfruR+p0cpub/ph8orHgAgL9tHSHqd\npOdFxHMk7SPpzDrLBAB9EhHfiojjIuI4Sc9X8qPbJ2suFgBUbmZiwfbhkk6X9BFJHrPIGZIulaSI\nuF7SJtubyywkACCXhyQ9KulA2xuVjJG8q94iAUBvnSrpuxGxs+6CAEDV8vRYeJ+kNyvpVjvOYZKy\nAfNOSYcvWC4AQEERcZ+k90r6Z0nfl/RARHyh3lIBQG+dqeRCeADQeVMTC7ZfLune9JoJ43or7Fl0\n5O/ZF24AAJTK9s9J+l1JRyi5KvPBtn+91kIBQA/Z3k/SKyT9Xd1lAYBl2Djj+RMlnZFeR2F/SU+y\n/dGIeG1mmbskbcn8fbjGdL21TbIBQCNFxLTEaZu8QMmtn34oSbb/Xkkc//jqAsRiAE3VoVgsJRc1\n/2ZE/GD0CeIwgKZaJA5P7bEQEW+PiC0RcaSS7lxfGkkqSNKVkl4rSbZPUNL1dteE9fViOv/882sv\nA3WlrtQ139Qxt0k6wfYBtq1kfO/20YXqbnO2Y+pKXanr6NRBZ0m6bNKTdbd3G6Y+bf+0E+3UhGlR\ns3osrIuDkmT7nDQoXhgRV9k+3fbtkh6RdPbCpQIAFBYRN9r+qJLbAz8u6QZJf1lvqQCgX2wfpCSx\n+7q6ywKsSn5vUCknkMA4uRMLEXGtpGvTxxeOPPfGkssFAJhDRPyxpD+uuxwA0FcR8YikQ+suBwAs\nU567QqCgwWBQdxGWhrp2U5/qiu7q03ZMXbupT3UFRrH950M75UM7Vc/L6g5jO+h6A6BpbCu6dcGw\nqYjFAJqoT7GYOIw6MBQCsywah+mxAAAAAAAA5lb04o0AAADokNVfMlfxiyYAoCh6LAAAAPRepBMA\nAMWRWAAAAAAAAHMjsQAAAAAAAOZGYgEAAAAAAMyNizc2wOhFk1Zx8SQAAAAAQNORWGiM0SRCL27l\nDAAAAABoOYZCAAAAAACAuZFYAAAAAAAAcyOxAAAAAAAA5kZiAQAAAAAAzI3EAgAAAAAAmBuJBQAA\nAAAAMDcSCwAAAAAAYG4kFgAAAAAAwNxILAAAAAAAgLnNTCzY3t/29ba32d5u+91jlhnYftD21nR6\nRzXFBQBMYvvnM3F4axqX31R3uQAAANBtG2ctEBE/sX1KROy2vVHSdbZPiojrRha9NiLOqKaYAIBZ\nIuJbko6TJNsbJN0l6ZO1FgoAAACdl2soRETsTh/uJ2kfSfeNWcxlFQoAsLBTJX03InbWXRAAAAB0\nW67Egu0NtrdJ2iXpmojYPrJISDrR9o22r7J9TNkFBQAUcqakv6m7EAAAAOi+vD0WHo+IYyUdLumX\nbQ9GFrlB0paIeK6kP5P0qVJLCQDIzfZ+kl4h6e/qLgsA9I3tTbavsL0jvT7ZCXWXCQCqNvMaC1kR\n8aDtf5T0AknDzPwfZR5fbfuDtp8SEWuGTKysrOx5PBgMNBgM5is1AMxpOBxqOBzWXYyqvUzSNyPi\nB+OeJBYDqFvHY/EHJF0VEa9Kr092UN0FAoCqOSKmL2AfKumxiHjA9gGSPivpXRHxxcwymyXdGxFh\n+3hJfxsRR4ysJ2a9V1/ZVjKaZM1c0V5A9WwrIjp1jRjbl0u6OiIuHfMcsRjAGmuPQ+o5/uhKLLZ9\niKStEfGMKcsQh7F0yX4uzi8w0aJxOE+PhadKujS9wvgGSR+LiC/aPkeSIuJCSa+S9Abbj0narWRs\nLwBgyWwfpOTCja+ruywA0ENHSvqB7UskPVfSNyWdm7kQOgB00sweC6W9EdnZieixANSnK7+S5UUs\nBjCKHgvlsf0CSV+VdGJEfN32+yU9FBG/n1kmzj///D2vYUgaloEeCxg1OiTtXe9610JxmMRCA5BY\nAOrTlYPZvIjFAEaRWCiP7Z+V9NWIODL9+yRJb42Il2eWIQ5j6UgsYJZF43Cuu0IAAAAAmC4i7pG0\n0/bR6axTJd1aY5EAYCkK3RUCAAAAwFS/I+nj6a1/vyvp7JrLAwCVYyhEAzAUAqhPV7rf5kUsBhaz\n2p14VRf2J4ZCLBdxGHVgKARmYSgEAADAUoXW/yAAAEB/kVgAAAAAAABzI7EAAAAAAADmRmIBAAAA\nAADMjcQCAAAAAACYG4kFAAAAAAAwNxILAAAAAABgbiQWAAAAAPSCbdmuuxhA55BYAAAAAAAAcyOx\nAAAAAAAA5kZiAQAAAAAAzI3EAgAAAAAAmBuJBQAAAAAAMDcSCwAAAAAAYG4kFgAAAAAAwNymJhZs\n72/7etvbbG+3/e4Jy11g+zu2b7R9XDVFBQDMYnuT7Sts70jj9gl1lwkAAADdtnHakxHxE9unRMRu\n2xslXWf7pIi4bnUZ26dLemZEHGX7FyV9SBIHsgBQjw9IuioiXpXG7YPqLhAAAAC6beZQiIjYnT7c\nT9I+ku4bWeQMSZemy14vaZPtzWUWEgAwm+1DJP1SRFwsSRHxWEQ8WHOxAAAA0HEzEwu2N9jeJmmX\npGsiYvvIIodJ2pn5+05Jh5dXRABATkdK+oHtS2zfYPuvbB9Yd6EAAADQbVOHQkhSRDwu6dj0l7DP\n2h5ExHBkMY++bNy6VlZW9jweDAYaDAZFygoACxsOhxoOh3UXoyobJT1P0hsj4uu23y/prZJ+P7sQ\nsRhA3ToeiwGgdxwxNgcwfmH7nZJ+HBHvycz7sKRhRFye/n2bpJMjYtfIa6PIe/WJba3PxVi0F1A9\n24qI0eRoK9n+WUlfjYgj079PkvTWiHh5ZhliMbCAtd/Z3fiubkKduhSLZyEO1yvZ3tWJfbeIvtYb\n+S0ah2fdFeJQ25vSxwdIerGkrSOLXSnptekyJ0h6YDSpAACoXkTcI2mn7aPTWadKurXGIgEAAKAH\nZg2FeKqkS21vUJKE+FhEfNH2OZIUERdGxFW2T7d9u6RHJJ1dbZEBAFP8jqSP295P0ndFTAYAAEDF\nCg2FWOiN6PY1EUMhgPr0qfutRCwGFtWEYQNla0KduhSLbd8h6SFJ/0/SoxFx/MjzxOEa9XVIQF/r\njfwWjcMzL94IAAAAILeQNIiI0Vu0A0BnzbzdJAAAAIBCOtH7AgDyIrGwRLbHTgAAAOiMkPQF29+w\n/bq6CwMAy8BQiKVbfy0FAAAAdMYLI+Ju2z8t6fO2b4uIL2cXWFlZ2fN4MBhoMBgst4QAem84HGo4\nHJa2Pi7euESTLtKY4OKNQB26dMGwPIjFwGKacKHDsjWhTl2NxbbPl/RwRLw3M484XKM+XcQwW9c+\n1RvzWTQOMxQCAAAAKIHtA20/MX18kKTTJN1cb6kAoHoMhQAAAADKsVnSJ9NfhzdK+nhEfK7eIgFA\n9UgsAAAAACWIiO9JOrbucgDAsjEUAgAAAAAAzI3EAgAAAAAAmBuJBQAAAAAAMDcSCwAAAAAAYG4k\nFgAAAAAAwNxILAAAAAAAgLmRWAAAAAAAAHMjsQAAAAAAAOZGYgEAAAAAAMxtY90FAFAf22PnR8SS\nSwIAAAA03+rxM8fLa5FYAHpvNCiOTzYAAAAAwDgzh0LY3mL7Gtu32r7F9pvGLDOw/aDtren0jmqK\nCwCYxvYdtm9KY/HX6i4PAAAAui9Pj4VHJZ0XEdtsHyzpm7Y/HxE7Rpa7NiLOKL+IAIACQtIgIu6r\nuyAAAADoh5k9FiLinojYlj5+WNIOSU8bsyj9pwGgGYjHAAAAWJpCd4WwfYSk4yRdP/JUSDrR9o22\nr7J9TDnFAwAUFJK+YPsbtl9Xd2EAAADQfbkv3pgOg7hC0rlpz4WsGyRtiYjdtl8m6VOSjh5dx8rK\nyp7Hg8FAg8FgjiIDwPyGw6GGw2HdxajSCyPibts/Lenztm+LiC9nFyAWA6hbD2IxAPSK89wmw/a+\nkv5B0tUR8f4cy39P0vOzY3xtR99vyZHcmmTSFfjXz+97e6F6k7bJPm17thURnRw6YPt8SQ9HxHsz\n83ofi4FFrI2b3YiXTahTl2PxKOJwvfp0q8BsXftU76p1tS0XjcN57gphSRdJ2j4pqWB7c7qcbB+v\nJGGR68JhtsdOAIBibB9o+4np44MknSbp5npLBQAAgK7LMxTihZJeLekm21vTeW+X9HRJiogLJb1K\n0htsPyZpt6QzixVj0q/4AIACNkv6ZJqc3Sjp4xHxuXqLBAAAgK7LNRSilDea0O2rT12xGQqBpunT\n/jdJn7rfSnTBBRbVhGEDZWtCnfoUi4nD9epqN/ZxGApRja62ZeVDIQAAAAAAACYhsQAAAABgLlwf\nDYBEYgEAAAAAACyAxAIAAABQItv72N5q+zN1lwUAloHEAgAAAFCucyVt1/orJLcCt4AHUBSJBQAA\nAKAktg+XdLqkj6jV91BvZU4EQE1ILAAAAADleZ+kN0t6vO6CAMCybKy7AAAAAEAX2H65pHsjYqvt\nwaTlVlZW9jweDAYaDCYuCnRadrhNBL1klmk4HGo4HJa2Pi/rA7Qd494r2ZhG57uTG9akuib60QZo\nlj7tf5PYVkS0uKtqMZNiMYB81sbNbsTLJtSpK7HY9h9Ieo2kxyTtL+lJkj4REa/NLNP4OLx3m5i9\nPayeGDa9TqvaVt5FZOva1HoX2daaoqltuahF4zCJhSUisYCm6dP+N0lXDmbzasMBLdBkTTgJL1sT\n6tTFWGz7ZEn/JSJeMTK/8XGYxEI3kFioRlPbclGLxmGusQAAAABUo1tnHgAwAddYAAAAAEoWEddK\nurbucgDAMpBYQC+Muw9z17ovAQAAAEAdSCygR7KJhMnDh8YlISQSEQAAAAAwDokFYKxJF9kEAADo\np65etA7A4rh4IwAAAAAAmBuJBQAAAAAAMDeGQgAAACxB9ho+dCUHgG7p+1AheiwAAAAsTT8POLHW\npAtFA0BbzUws2N5i+xrbt9q+xfabJix3ge3v2L7R9nHlFxUAMIvtfWxvtf2ZussCAACAfsjTY+FR\nSedFxLMlnSDpt20/K7uA7dMlPTMijpL0ekkfKr2kAIA8zpW0XfwsCqCBbO+ZAADdMTOxEBH3RMS2\n9PHDknZIetrIYmdIujRd5npJm2xvLrmsAIApbB8u6XRJHxH3SAXQWOQ9AaBrCl1jwfYRko6TdP3I\nU4dJ2pn5+05Jhy9SMABAYe+T9GZJj9ddEAAAAPRH7rtC2D5Y0hWSzk17LqxbZOTvdenolZWVPY8H\ng4EGg0HetweAUgyHQw2Hw7qLUTrbL5d0b0RstT2YtiyxGEATZGMRAKDdnOd2GLb3lfQPkq6OiPeP\nef7DkoYRcXn6922STo6IXZllYtx7JWPsRue7k7fpmFTXRD/aoC7r235y+7JNdrOuk9hWRLR+2IDt\nP5D0GkmPSdpf0pMkfSIiXjuy3NhYDCCftXGzWLzc+9pmxdlF6jT/e619n67E4jxsh9TsW9JN+5wk\nzZzXZG0r7yKydW1qvZsaF6eZ1JZNbeO8Fo3Dee4KYUkXSdo+LqmQulLSa9PlT5D0QDapAACoVkS8\nPSK2RMSRks6U9KXRpAIAAGgGLmKKrskzFOKFkl4t6SbbW9N5b5f0dEmKiAsj4irbp9u+XdIjks6u\npLQAgLzamS4HAAAztf3XcXRPrqEQpbwRQyEYClEjhkKMN32bXK+rbdCX7rcSQyGARTEUoqz3YihE\nk7aBUQyFaH458ryeoRDVYCjEeLkv3gigb/InHAAAAAD0V6HbTQIAAAAAAGSRWAAAAAAANB4XvWwu\nEgsAAAAAgNYi4VA/EgsAAAAAAGBuJBYAAAAAAMDcuCsEAABopWy317be3gsAgC4gsQAgt0lj1zig\nR5+N7hfsD8uW3P8caALb+0u6VtITJO0n6dMR8bZ6SwUA1SOxAKCg0ZMmDuiBvfsF+wPQZxHxE9un\nRMRu2xslXWf7pIi4ru6yAUCVuMYCAAAAUJKI2J0+3E/SPpLuq7E4ALAUJBYAAACAktjeYHubpF2S\nromI7XWXCQCqxlAIAAAAoCQR8bikY20fIumztgcRMRxdbmVlRZI0GAw0GAyWWUQA0HA41HA4LG19\nXtZFpmxPeaP1Y7a7ePGr5AJfk8an96MN6rK+7Se376TPqYufR9FtctL8NreNbUVEbwbG2442f15N\ntHY/avf+0DZ723557b7I511HefNY5jY8qQ26Gottv1PSjyPiPZl5ITX7Qq/TPidJM+c1WVPKu2g5\n8rw+u0xT6j2qaFycVI9l1q8JZajConF4yUMhYmQCAAAAusH2obY3pY8PkPRiSVvrLRUAVI+hEAAA\nAEA5nirpUtsblPyA97GI+GLNZQKAypFYAAAAAEoQETdLel7d5QCAZSOxAAAAOmt1zKvU3nGvAAA0\nHYkFAADQccmFwQAAKILkdH5LvngjAAAAAABtQUIhj5mJBdsX295l++YJzw9sP2h7azq9o/xiAugS\n22MnLMb2/ravt73N9nbb7667TAAAAOi+PEMhLpH0Z5I+OmWZayPijHKKBKAfRrO/JBYWFRE/sX1K\nROy2vVHSdbZPiojr6i4bAAB9tPrDCd3o0XUzeyxExJcl3T9jMc4IAKABImJ3+nA/SftIuq/G4gAA\nAKAHyrjGQkg60faNtq+yfUwJ6wQAzMH2BtvbJO2SdE1EbK+7TAAAAOi2Mu4KcYOkLWnX25dJ+pSk\no8cvupJ5PCjhrQGgmOFwqOFwWHcxKhMRj0s61vYhkj5rexARw+wyKysrex4PBgMNBoNlFhEAJK2N\nRQCAdnOe8T62j5D0mYh4To5lvyfp+RFx38j8mDymev38Lo5DSsZY9bsN6rK+7Se376TPqYufR9Ft\nctL8om3TpDa2rYjo5HAu2++U9OOIeE9mXnRxW67T2u25m7Giqfa2fZ6YXs5ns8jnXXZZyrLMbXhS\nG3Q5Fo9KjombPeZ+2uckaea8Jlt2eSe936LlyPP67DJN/ZyKxsWq2rNI+ZZdhmVZNA4vPBTC9man\nrWj7eCXJCsb0AsCS2T7U9qb08QGSXixpa72lAgAAQNfNHAph+zJJJ0s61PZOSedL2leSIuJCSa+S\n9Abbj0naLenM6ooLAJjiqZIutb1BSeL4YxHxxZrLBAAAgI7LNRSilDdiKARDIWrEUIjxyhoKMUkb\n2rhP3W8lhkJUgaEQ9WEoRDkYCrFcDIWoF0MhmvU5MRSiORaNw2VcvBEAVCThAAAA5tP2kxcA3URi\nAY2w+iU5ii9NAAAAAGg2EgtoEH7xBgAAAIC2WfiuEAAAAAAAoL9ILAAAAAAAgLmRWACAnrE98bom\nAAAAQFEkFgCgd7goKgAAAMpDYgEAAAAAAMyNxAIAAAAAoDCGV2IViQUAAAAAADC3jXUXAJjHuMxo\nBOPGAQAAAGDZ6LGAFovMBABA86x2E6arcD/Y3mL7Gtu32r7F9pvqLhMALAM9FgAAACoVkkgs9MSj\nks6LiG22D5b0Tdufj4gddRdsXqtJMXqGApiGHgsAAGAh2V/l+WUefRYR90TEtvTxw5J2SHpavaUC\n+qXN30Vt/i4lsQAAAErA0DQgy/YRko6TdH29JQHQLu38LmUoBAAAAFCidBjEFZLOTXsurLOysiJJ\nGgwGGgwGSysbUCeG1jTHcDjUcDgsbX1e1odqO9ZnX1a7eKyf38WNLdmR+t0Gk0xqm0ltsH756e1V\nZPmiZWmzottk0fltaGPbioj29Teb095Y3M1tug5rt+dut+to18zVutbVBnvfN0+8KadcReuaff9E\n8/a/ZX5+kz6PLsVi2/tK+gdJV0fE+8c8H9L8J1bLODFbv90m7zfuvdt2orjs8k56v0XLkef12WWq\nqnfR9Y4uXzRG52nPMus6rnxVfaaT3m9ZFo3D9FgAACxN9sS0LQehy9Kettl7AgpgLSc78kWSto9L\nKvRV25IPAIrjGgsAgCXjwHIy2gZouRdKerWkU2xvTaeX1l0otF9bL+jXJbT/dDN7LNi+WNK/lXRv\nRDxnwjIXSHqZpN2SfjMitpZaSgDATLa3SPqopJ9Rcob6lxFxQb2lAoD+iIjrxA93GINeG9XhhL8Z\n8gS+SyRNzLTaPl3SMyPiKEmvl/ShksoGAChm9f7pz5Z0gqTftv2smssEAADmQC8FtMnMxEJEfFnS\n/VMWOUPSpemy10vaZHtzOcVD343eG50AC0zG/dP7hZgIAFg2vnfaoY7PqYyuWodJ2pn5+05Jh5ew\nXiAVIxOAWbh/ev2Wc+JPTAQANFeV34NVf8eSQCmmrLtCjLb6hCOdlczjQUlvDQD5lX3P3iaaff/0\nleTflZWp90+fdGtBTLe23dberg3AXisrK3UXAQBQEue8X+gRkj4z7uKNtj8saRgRl6d/3ybp5IjY\nNbJcjLtvfaIZ97Ov2tr7RO+Zm/7fjzaYpGjbrJ8/vb3Wrz/PPc/XvmcXP4/F2336/Da0cZfunS7l\nvX96vvsjV3Fv+zrvz1yV0XZam0woem/u6pYvw6RtooptpVh58sSbsrdh5Vpn9v0Tzdv+l/n5Tfo8\nuhaLp0ni8PzJ2mVcBHD9dpu837j3zlOeJl24cNllmZSkL9puo8sv+vp5FS3HrOXHbWt51zeuR0EV\ndc2Wb9r6y3jfsr6z5tkvF43DZQyFuFLSa9PCnCDpgdGkQpkYbw8A41V9/3TiLgAAqNbyEj4cz5Qr\nz+0mL5N0sqRDbe+UdL6kfSUpIi6MiKtsn277dkmPSDq7ygIn1v7yDCzDpODThOw7kFq9f/pNtldv\n+/u2iPif5aw++6t8ubL7F/sUAADVaFLvEXTLzMRCRJyVY5k3llMcLAsnyfOaNAwAqF/775/O9Qj6\nokuJpC7VBc0ybtvq6y+snAwDzVfWxRvRYJO/hDhJBgCpPReqLPukot56tyeRNL3dq+vJ01Rt2V+6\noT37Cfq//jvIAAAQ1klEQVSnr4kujEdioTeqGT5CzwcA3bH+BLGZJ1Bln2j078R4PrTTWrRH29EL\nAECZSCygBPR8GKfqpMukK+E2CZlsVGW53c85gQIAAJiGxAJQqaqTLm24kCmJJ1SFLsIA2oseA+3D\nZwZM1uKLfAEAAABtxckpgO6gxwIAAGg97s6APuOX9MXRhsBiGptYYGw2AAAohuExQNdNSgCQGADq\n1djEQoKx2QDQdstMFPOrNaZh+wBQFhIZmEeXv4canlioRttvkdiGuwEAwFrL/CV58fei11yX0ash\nr7Jut8r+tDjbHOsBaLReJhYSbe8N0Ya7AQBAW3GLSSBR1r5AQgdoInpeoCydSSzwKz4AAOiSLneZ\nxXj07ugWPs/m61tipcr6diaxkOBXfABAoisnZRyY9h2/9PcPnzlQhfEn1exvZdlQdwEAAKhOexMK\na3WlHgD6jEQpls02292SdKzHAgAAKFPXD8jKukAhgMWN/qLc9fjTNH0bFtB0bet5SWIBAADM0PWL\nWba3fk048GxCGbpg1kl0f0+y6areR02NK8vfD9uz/TMUAgAANMJql9W+nUAtXu+6D7qjAWVoDtsX\n295l++a6ywK0G3GlTUgsAGiM7MF1X08wUBzbStf09UCyr/XupEskvbTuQnQdMR9oFoZCAA3Xvy/O\n0YPrvtUf82lPV0EAia5+v0XEl20fkXf5LrVDk8boN6ksy9LHOvdFGz7bXD0WbL/U9m22v2P7LWOe\nH9h+0PbWdHpH+UUF+ixGJsyj670h6H4LYJbmxT++0xLF26FZnyNQXBXbL/tEfWb2WLC9j6Q/l3Sq\npLskfd32lRGxY2TRayPijArKCKDnJn1JzJe1zb6mc18+l0j6M0kfrbsgAJqMHj7NsFJ3AQop+4Rt\nmRfna8OvvUBdVlZWSllPnh4Lx0u6PSLuiIhHJV0u6ZVjluvsN1TRcd+MEweqQK+NWSLiy5Lur7sc\nAIA8VuouQANM/j7n2Bkox6xz0ZWVlVKSC3kSC4dJ2pn5+850XlZIOtH2jbavsn3MwiVrnKInNf0+\nCSK5AgDjERMBtB0xDChHl/alPImFPGfFN0jaEhHPVdIN91MLlQod0e/kCgBMtpyYSBIDEtvBstm+\nTNJXJB1te6fts+suE7prWfs2MQSz5LkrxF2StmT+3qKk18IeEfGjzOOrbX/Q9lMi4r61q1rJPB4U\nLCqARfBlsGqYTn22kvy7sqLBYKDBYFBradqoXfsT4+khjdsOmrEdr9RdgNJFxFl1lwEoQzNiBNrC\nsy5iYnujpG9JepGk70v6mqSzshdvtL1Z0r0REbaPl/S3EXHEyHpi8m3kFplvRUS64a+fP6FOY9dd\n9fLryz59fv71T1520vJF3zP/uudb/+JtVrwNyipLWdvNJPnLU9520K75a83TBhHRmW/O9BZnn4mI\n50x4PlZPMHLEf+1tr7WPi2zHe9eT/SzXriO7TJ51T1rnpMeL1HVcvWcvX7yueeYXrd+sdh2tR55l\n8rRH/vYet+stvh3kab8iZS+6fSxSxqL7RdFlsmWY93MtGhvy7q9disXTJHFYKhLDptl7HFzs4oTT\nXjN9H5383pPWOe5CjeO23XHrGtlGJl70cfS9J5VrVhlnH0fkb/NxJ+bjyjfP+42ue9oyo+02a5lx\nn9Es0+PPrDLl/46bVqdxZZq1Hcz+jhwt09r1T6vrOKPb7+w2X/veeffxPNv/uPdeJA7P7LEQEY/Z\nfqOkz0raR9JFEbHD9jnp8xdKepWkN9h+TNJuSWfOW6A2IYsHNMH6REFfpd1vT5b0U7Z3Svr9iLik\n5mIBOeQ7cC0T3+EAAJQnz1AIRcTVkq4emXdh5vFfSPqLcovWFvl+RQVm4SAXi6L7LVBE9hd3oDvm\n6b2AYqo+ZuOYcLay2yj9tb7UdfZNrsQCgGUhUQUAANqNk7QyLL8nF/qprCQNiQUAAHqEX8IwC9sI\nAOTTjR5C5SSxSCwAAFqBk50yMQygSu3fVtk+umqRCxf2SR/ag14lKBuJBQBAi9A1FNWadMX54thW\ngWUoksibN+nX/mRhu82b6Kn7c2tygqqKtiGxAKBzpt3eCeWadPuqstbJ54bl49d6oM/qPhkF2mpD\n3QUAgGpEZkK1qmhnPjcAQB32fv/Y7kSioQl1aEIZUC0SCwAAoHKrB+gcXAJoj/KT3IvGwDJjKPEY\nZWIoBAAAWJL5rjtQ1/AYDrrRBpO2Uy7Oh9m4FgzKQ2IBAFA6TsjyoZ2KqOMAmOstAPNq8oXrqrBI\nffvWVk3Whs+iqWUksQAsiBMDYBJ+CcmHdgLQHGUf1zS9636TjuMWKUuT6oH1lvH51J1w4BoLQClC\nXCwQAAA0Sfa6JlWf2LT7xLa8Y7fFriXDMSTaix4LyKXdXxYAgDzacKtPvo8ma8Pnh2Zj/0Kfsf0v\nhsQCChg9SGHnA7pmmV+qfIFPVl/btOWaAgwfmYy2wSLq2X74Piimr+1Vfb3ri59d+ExJLKBTxu2U\n/GoDFLXMk8vpX+Jd+KKdHyeIAKrRzNjarZjXzDaG1NTPpv3bP4kFrNHMHa2obCKhefXpRhsDy7T4\nly1dxNE0fBegXu0/iUGzEeP6h8QCxmDIQ/VoY2D5OJDui/Yc0LJNAlXIGwPaEyvQdGxLJBYAAEAD\nLX6Q1u+Tdg5ykdWn7aFtdW1bedEcdd9echSJBQBAYRwItU87P7N+JwcWR/sBkzQnJrKfLkMVn3dz\ntqG16ioXiYUGa+rGuogu1qlpuIAllqMtdw/oh/zXkOAAtkp8xwHLtdg+V048tM1xFiBpw6wFbL/U\n9m22v2P7LROWuSB9/kbbx5VfzD6LkakLulinpqF9+ypPzO4723umbkn29+7Wrw2IuyAOl6HM+EUs\nXI483zt8Fmstoz2W2eZTEwu295H055JeKukYSWfZftbIMqdLemZEHCXp9ZI+VFFZ55I9wJq1wY9b\nlh2gu/i80TV5YvbI8j3d7vtw8rd4/fIeJOb9XkV5aNfmKhqH+yjfttv1GF1MHSegi7xnntd2PX7l\nq99yt/Oq23xWj4XjJd0eEXdExKOSLpf0ypFlzpB0qSRFxPWSNtneXHpJF1Lk11t+Te8XPm90Sp6Y\nnbF3uy/rJIWTnfHaeSKYJy7mXQaLWr8N0a4NVTAOl6tdMQb5Fdvf59sOlh1TiGFdMyuxcJiknZm/\n70znzVrm8MWLhqL49R3ovTwxe4KykmvtStItN1a2p13QVM3fhjgGWSQOo626tr13rT7LVkf7NeEz\nm3XxxrzfYKM1af43X2dlm77+DQzAUhF7Cyt2EcomfHG3Ae00vya3XbGuvb1NLrQuDlfd5R3rNb/d\n5r+wZfPrlg/X+ShuVmLhLklbMn9vUZJ5nbbM4em8MSY16mLz935Y+eZXWZa65tMGtIFEG0jztkFn\n5InZ2tsGHjOvD48XWYdyLNOGx1W/j3IsU/fjut9/2mPlWKaONiharkWWb62ccVharP0XVfRzqkpT\n2qDO7bOacuQ7maX9i8evZbz35DKN+1wXT1yUU79ZiYVvSDrK9hGSvi/pP0g6a2SZKyW9UdLltk+Q\n9EBE7BpdUUR07psDABpmZswmFgNApYjDAHppamIhIh6z/UZJn5W0j6SLImKH7XPS5y+MiKtsn277\ndkmPSDq78lIDANaZFLNrLhYA9AZxGEBfOaJ1Q8EAAAAAAEBDzLorRClsv9T2bba/Y/sty3jPZbF9\nse1dtm/OzHuK7c/b/rbtz9neVGcZy2J7i+1rbN9q+xbbb0rnd66+tve3fb3tbba32353Or9zdV1l\nex/bW21/Jv27k3W1fYftm9K6fi2d18m6ZnU8DvcmNq3q0f66yfYVtneksfgXO1zXt6Xb8M22/8b2\nE7pS16LHSmlbfCeNWafVU+pqdDkWL6Kv382zsO/kM6GdVmzfmW5TW22/LPNcX9up8PFSkbaqPLFg\nex9Jfy7ppZKOkXSW7WdV/b5LdImSumW9VdLnI+JoSV9M/+6CRyWdFxHPlnSCpN9OP8vO1TcifiLp\nlIg4VtIvSDrF9knqYF0zzpW0XXuvaN3VuoakQUQcFxHHp/O6WldJvYjDvYlNGX3ZXz8g6aqIeJaS\nWHybOlhXJ+PxXyfpeRHxHCVd6M9Ud+qa+1jJ9jFKrktwTPqaD9peyg9hVetBLF5E776bc2LfyWdc\nO4WkP023qeMi4mqp9+1U6HipaFstoxGPl3R7RNwREY9KulzSK5fwvksREV+WdP/I7DMkXZo+vlTS\nv1tqoSoSEfdExLb08cOSdii5N3NX67s7fbifkoO8+9XRuto+XNLpkj6ivZeG7WRdU6MXzupyXaXu\nx+Fexaa+7K+2D5H0SxFxsZSMXY+IB9XBukp6SMkB34G2N0o6UMmF/zpR14LHSq+UdFlEPBoRd0i6\nXUkM64JOx+IS9O27eSb2nXwmtJM0/nYHfW6nosdLhdpqGYmFwyTtzPx9ZzqvyzZn7oyxS9LmOgtT\nhfTXleMkXa+O1tf2BtvblNTpmoi4VR2tq6T3SXqzpMcz87pa15D0BdvfsP26dF5X67qqN3G4D7FJ\n/dlfj5T0A9uX2L7B9l/ZPkgdrGtE3CfpvZL+WUlC4YGI+Lw6WNeMSXV7mtbenrFL8ao3sXgOffxu\nnlcf9515/Y7tG21flOneTzsp9/FSobZaRmKh11eHjOTqmJ1qA9sHS/qEpHMj4kfZ57pU34h4PB0K\ncbikX7Z9ysjznair7ZdLujcitmrCjWy7UtfUCyPiOEkvU9IF7JeyT3asrqu6Vp+x+hCbera/bpT0\nPEkfjIjnKbnz1Jqu0F2pq+2fk/S7ko5QciB3sO1XZ5fpSl3HyVG3rtS7K/WoQh+/mxfWo31nHh9S\nkqA+VtLdSpK3k/SqnRY8Xpr43DISC3dJ2pL5e4vWZj66aJftn5Uk20+VdG/N5SmN7X2VbIgfi4hP\npbM7W19JSrve/qOk56ubdT1R0hm2vyfpMkm/Yvtj6mZdFRF3p///QNInlXTp6mRdMzofh3sUm/q0\nv94p6c6I+Hr69xVKEg33dLCuL5D0lYj4YUQ8JunvJf0bdbOuqyZts6Px6vB0Xhd0PhbPq6ffzfPq\n475TWETcGyklQwdXu/D3up0KHi8VaqtlJBa+Ieko20fY3k/JBSCuXML71ulKSb+RPv4NSZ+asmxr\n2LakiyRtj4j3Z57qXH1tH7raZcr2AZJeLGmrOljXiHh7RGyJiCOVXCzsSxHxGnWwrrYPtP3E9PFB\nkk6TdLM6WNcRnY7DfYpNfdpfI+IeSTttH53OOlXSrZI+o47VVclFKU+wfUC6PZ+q5OKcXazrqknb\n7JWSzrS9n+0jJR0l6Ws1lK8KnY7F8+rxd/O8+rjvFJaeIK/690q2KanH7TTH8VKxtoqIyicl3Zq+\npeSCD29bxnsua1Lyi9H3Jf1fJePmzpb0FElfkPRtSZ+TtKnucpZU15OUjOndpuQke6uSK4R2rr6S\nniPphrSuN0l6czq/c3UdqffJkq7sal2VdInblk63rMajLtZ1TN27HId7E5tG6t3p/TWt13MlfV3S\njUp+xT+kw3X9PSWJk5uVXDxr367UteixkqS3p7HqNkkvqbv8JbdFZ2PxAm3S2+/mHG3DvjNfO/1H\nSR9Vcgx/o5IT5c20U/HjpSJt5fQFAAAAAAAAhfXlnp0AAAAAAKACJBYAAAAAAMDcSCwAAAAAAIC5\nkVgAAAAAAABzI7EAAAAAAADmRmIBAAAAAADMjcQCAAAAAACYG4kFAAAAAAAwt/8PFd6Mn1kgU9QA\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10dc9f5d0>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Dirichlet Process Mixture#"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Dirichlet Process Mixture example\n",
      "\n",
      "Likelihood model = two-dimensional normal distribution\n",
      "\n",
      "###Test data generate from two-dimensional normal distribution###"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import multivariate_normal \n",
      "\n",
      "n = 60 #n should be multiple of three\n",
      "x = np.zeros([n,2])\n",
      "color = ['r','b','g','m','c','y','k']\n",
      "n_c = n/3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First cluster:\n",
      "\n",
      "$x \\sim \\text{Mulrivariate-Normal}([2,2], 0.5*I)$\n",
      "\n",
      "Second cluster:\n",
      "\n",
      "$x \\sim \\text{Mulrivariate-Normal}([-2,-2], 0.5*I)$\n",
      "\n",
      "Third cluster:\n",
      "\n",
      "$x \\sim \\text{Mulrivariate-Normal}([2,-2], 0.5*I)$\n",
      "\n",
      "\n",
      "We randomly generate 20 data points for each cluster."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# we assume three clusters\n",
      "# mean and variance of original data \n",
      "mu1 = [2,2]; mu2 = [-2,-2]; mu3 = [2,-2]\n",
      "sigma = np.identity(2) * 0.5\n",
      "\n",
      "# randomly generate n/3 data point per cluster\n",
      "x[:n_c,:] = np.random.multivariate_normal (mu1, sigma, size=n_c)\n",
      "x[n_c:n_c*2:] = np.random.multivariate_normal (mu2, sigma, size=n_c)\n",
      "x[n_c*2:,:] = np.random.multivariate_normal(mu3, sigma, size=n_c)\n",
      "\n",
      "y = [0]*n_c + [1]*n_c + [2]*n_c"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# plot randomly generated data to 2-dim space\n",
      "c_mask = [color[y[i]] for i in xrange(n)]\n",
      "plt.scatter(x[:,0],x[:,1], c=c_mask)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "<matplotlib.collections.PathCollection at 0x10dfd6710>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEACAYAAACqOy3+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX+x/H3ycwkM5OQBEMSehFIIIB0RBEcRJReRFdc\nlLUtKiqsa28r/kRddYW1gA0sqCiIUiIiIBJ6CyWBUAQESUApAYVkUifn90dGFiWkzWTuTPJ9PU+e\nJ3fmlk+AfDlz7jnnKq01QgghAk+Q0QGEEEJUjhRwIYQIUFLAhRAiQEkBF0KIACUFXAghApQUcCGE\nCFBeKeBKKZNSaqtSKtEb5xNCCFE2b7XAxwM7ARlULoQQPuJxAVdKNQQGANMA5XEiIYQQ5eKNFvhk\n4GGgyAvnEkIIUU4eFXCl1CDgmNZ6K9L6FkIIn1KerIWilHoBuAUoBKxAOPCl1nr0OftIv7gQQlSC\n1rrUhrFHLXCt9RNa60Za62bASOD7c4v3Ofv5/dczzzxjeAbJKTkDNaPk9P5XeXh7HLi0toUQwkfM\n3jqR1noFsMJb5xNCCFE6mYnp5nA4jI5QLpLTuwIhZyBkBMlpBI9uYpbrAkrpqr6GEEJUN0opdFXe\nxBRCCGEcKeBCCBGgpIALIUSAkgIuhBABSgq4EEIEKCngQggRoKSACyFEgJICLoQQAUoKuBBCBCgp\n4EIIEaCkgAshRICSAi6EEAFKCrgQQgQoKeBCCBGgpIALIUSAkgIuhBABSgq4EEIEKI8KuFLKqpTa\noJTappTaqZR60VvBhBBClM6jhxprrXOVUr211k6llBlYrZS6Qmu92kv5hBBCXIDHXShaa6f722DA\nBJz09JxCCCHK5nEBV0oFKaW2AUeB5VrrnZ7HEkIIURZvtMCLtNYdgIZAL6WUw+NUQgghyuRRH/i5\ntNa/KaUWAl2ApHPfmzBhwtnvHQ4HDofDW5cVQohqISkpiaSkpAodo7TWlb6gUqoOUKi1/lUpZQMW\nA89qrZeds4/25BpCCFETKaXQWqvS9vG0BV4P+EgpFURxd8zH5xZvIYQQVcejFni5LiAtcCGqraKi\nIj7++GO2b9lCy4QE7rjjDsxmr/XM1mjlaYFLARdCVNqYW24h9auvGO508q3dTh2Hg9lff41SpdYd\nUQ5SwIUQVSY9PZ2OcXEczM0lDMgD4u12FqxbxyWXXGJ0vIDniz5wIUQNobXmvXfeYe2yZdRr0oRh\nN9xAhMlEqPv9ECDabCYrK8vImDWKtMCFEOXy8LhxrHj/fe7KzmZDcDBr6tdHASPS07nZ5WJBUBBT\noqPZvn8/oaGhZZ5PlE66UIQQXlFQUEAtu50jhYVcBGjg6rAwRk6axNeff05qairxcXFMmTGD5s2b\nGx23WpAuFCGqiaKiIo4fP05ERARWq9Xn13e5XGitsbu3FRAO1KpVi/nLlrF7925+/PFHTCaTz7PV\nZLIeuBB+bt++fbRt1ow2TZtSJyKCd6ZO9XkGq9XK0H79GGW1shaYFBTERouFPn368J8XXsDRqROv\n/fWvdE1I4NOPP/Z5vppKulCE8HOd4+MZvXcv47VmP9DLbmf+ihV06dLFpzlycnJ46qGHWL1sGfUb\nNeKlKVNQStGjfXu25eRQH9gJXG61cujoUcLDw32ar7qRLhQhAlxhYSHb9u5lo7sR1BzoDyQnJ/u8\ngNtsNl6dMuUPr3333Xe0CQ6mfk4OAAlAlNnMzz//LAXcB6QLRQg/ZjabqVe7Nqvc2znAxqAgmjRp\nYmSss1q3bs32ggK2uLcXA1lBQTRu3NjIWDWGFHAh/Nz7n33GDXY7Q8LDuSQ0lC4DB9KvXz+jYwHQ\noEED3p0xgz42G41DQxkdHs4XiYnYbDajo9UI0gcuRABIT09n8+bNxMTEcNlll/ndVHWn08nRo0ep\nX78+ISEhRsepFmQcuBBCBKjyFHDpQhFCiAAlBVwIIQKUFHAhhAhQUsCFECJASQEXQogAJQVcCCEC\nlMcFXCnVSCm1XCmVppTaoZQa541gQgghSueNFngB8IDWug3QHbhXKdXaC+cVQhhgwYIFtL/4YppG\nRzPu738nLy+vzGNOnz7N8889x/1jxjB79mxk7odveFzAtda/aK23ub/PAnYB9T09rxDC9zZu3Mjf\nR47kPwcOsPTECfZ/+ikPjh1b6jFOp5Mru3Rh5/PP0/y995hw2208P2GCbwKX4ZMZM+h5ySX0at+e\nz2bONDqO92mtvfYFNAV+AsLOeU0LIQLDv556Sj+hlNagNej9oBvWrl3qMbNnz9ZXhYXpIvcxGaBt\nFot2uVw+Sl2yWZ9/rpva7XoR6G9AN7bb9Zdz5hiaqSLctbPUmuu15WSVUmHAHGC8Lm6JnzXhnP+N\nHQ4HDofDW5cVQnhRrYgIUi0WyM8HIB0IK+P5lk6nk2it+X3OdxTgKirC5XIRFHThD/kul4vnnn6a\nLz7+GLvdzpMvvcSwYcO884MAn0ydyktOJ78v+/Wi08mnb7/NdSNGeO0a3pSUlERSUlLFDiqrwpfn\nC7BQvJLkP0p4zwf/VwkhvCEzM1O3qF9f/y04WE9QSte12/UXs2eXekx6erqOrlVLTwedAvqmkBB9\n3bXXlnmtCU88oS+32/UW0ItA17Xb9cqVK731o+gR11yj33N/KtCgp4IeOWiQ185f1ShHC9zjxaxU\n8bJoHwGZWusHSnhfe3oNIYTvZGZmMm3aNE6fOsWAwYPp0aNHmcds27aNB8eM4ecjR+jVpw//mTKF\nsLCwUo9p1aABE48cYSBgA/4NHL/vPl59443z9l26dCm7d+8mISGBPn36lOvnWL16NcOvuYZHc3Io\nAl6x20lctozu3buX63ij+WQ1QqXUFcBKIJXih1UDPK61/tb9vhRwIUqwdOlSPnzzTUxmM3c/+CCX\nX3650ZF85u033+Sf48ZRT2ucwFxglslE2COP8NwLL/xh30fHj2fe9On0KSzkO7OZ6++6ixdefbVc\n19mwYQMfvvUWKiiI28eO9flTjDwhy8kK4ae++eYb7rjhBp5zOskDJthszP/uuxpRxNPS0ujTtSvr\ncnJoBiQCowFbZCQbUlNp1KjR2X0PHDhAt4QEfsjNpTZwEoizWtmyZ0+1f+qPLCcrhJ+a8uKLTHY6\nuRO4F3gmJ4e3y9mqDHRpaWn0MJtp5t4eDLhMJhYlJf2heAOcOHGCRsHB1HZvXwQ0CA7mxIkTPkzs\nv6SAC2EAV2EhwedsBwNFLpdRcUpUVFTEsWPHKCws9Op5W7RowQaXi2Pu7dVAsM1G27Ztz9u3devW\nHDOZmAnkA58AJ81m4uPjvZopUEkBF8IAdzzwAA/Y7XwJzAT+ZbNx2/33V+k1ly9fzqBevbi6a1em\nv/deqbMlN2/eTNPYWFo3bkx0RATz583zWo5OnTpx94MP0tZmo2dEBMNDQ/l0zhxMJtN5+4aFhZG4\nbBkvNG2KPSiIl5o14+vvvye0jKGNNYX0gQthkDlffMH0yZMxmc3c+/jj9O/fv8qutWHDBgZfdRWT\nnE5qAw/a7TzwyivcVcIsy4KCAi6uV49JmZncAGwCBtjtbNm9+7wuDk/s27ePjIwMEhISiImJKXN/\nrbXfPQu0KslNTCEEAPePGUOj997jEfd2EvBofDwbdu8+b9+DBw9yRZs2ZDidZ1+7NiKC8TNnMmDA\ngCrNmZuby4EDB4iJiSEqKqpKr+Xv5CamEAIAk8lE7jmt11z3ayWJjo7mtMvF76X9JLCjoMCrre+S\nbN68mRYNGjC0WzcubtCAyS+/XKXXqw6kgAtRA9x57728YbfzCjAd+Lvdzvinny5x39DQUN58+22u\ntNsZXqsWHex2bhs7lnbt2lX4upmZmaSlpeE8pzVfEq01NwwcyKSTJ/khK4u0vDxeffZZNm/eXOFr\n1iTShSJEDZGSksIbL79MntPJyDvvZODAgaXuv2fPHlJTU2nWrFmlJsC8+d//8uRjj1E3OJgzJhPz\nFi+mW7duJe6bnZ1NVEQEOS7X2TVVbg4N5eo33+TWW2+t8LWrA+kDF0IYIjU1lWu7d2d9Tg5NKJ5p\nOS4qikPHj5d4I1JrTYOoKD46dYq+FHfbdLbb+WTJknJN5a+OpA9cCGGItLQ0eprNNHFvDwd+O3OG\nX3/9tcT9lVLM/OorRoWFcUVEBAk2G6PGjq2xxbu8vLacrBDCP508eZJ7Ro9m9Zo11I+N5fUPPuCy\nyy6r0mu2aNGC9S4Xx4Foike9WK1WIiIiLniMw+Fg54EDpKWlUa9ePeLi4qo0Y3UgXShCVHPX9OhB\ni+RkHsvPZwNwb1gYyWlpVb6WyLNPPMGUyZOJCwlhT2Ehn8+fX+6VBIX0gQtR4+Xk5BBZqxbZLtfZ\nj9sjw8IY9NZb3HzzzVV+/b1793LkyBESEhKIjo6u8utVJ+Up4NKFIkQ1FhwcjCkoiCMuF42BIuAQ\nUKtWrRL3X7duHevXr6dBgwaMGDHigmPFy6tly5a0bNnSo3OIC5ObmEJUYyaTiYkTJ9LbbudZYIjN\nhjkursRp+29PmcINV1/NgcceY9Ltt3P9gAEUFRX5PrQoN+lCEaIGWLRoEWtWraJ+w4bcfvvtWK3W\nP7xfWFhIRGgoqfn5NAcKgM5hYbz61Vf07dvXkMw1nXShCCEA6N+/f6mLZWVnZ0NRERe7ty1AvFKy\n7rafky4UIQQREREktGzJsyYT2cBSIMnlqvLhhsIzUsCFEADMXbKEFR07EmU2c3dsLJ/Nn0/Tpk2N\njiVK4Y2HGr8PDASOaa3PW+1G+sCFEKLifDWV/gOgnxfOI4QQogI8LuBa61XAKS9kEUIIUQHSBy6E\nEAHKJ8MIJ0yYcPZ7h8OBw+HwxWWFEAHg4MGDbNu2jYYNG1Zq3fGqoLVmxowZpGzaRPPWrRkzZgwW\ni6VKr5mUlERSUlKFjvHKRB6lVFMgUW5iCiEqYv68edw5ahTdzWa2FxYyfPRoJr/1ltGxuPf229k0\naxY3OJ0stdkI6d6d+d99R1CQ7zotfLaYlRRwIURFuVwu6oSHs8TppCtwGugYGsqMxYsNXQf8l19+\noVWTJhzKzyec4lmpbUJD+eT77y/4RKGq4JNRKEqpz4C1QJxSKl0pdZun5xRCVH+nT5/GVVhIV/d2\nONA5KIhDhw4ZGQun00kts5nfl/uyANEmU/FsVT/jjVEoN2mt62utQ7TWjbTWH3gjmBCieouMjCQm\nKoqP3Nu7gBWFhbRv397IWDRp0oTYxo15zGxmDzA5KIjDISF07tzZ0FwlkVEoQghDKKWYu3gxz8bG\nEmO1cmlICP956y0SEhIMzWUymfh6+XL29enDwNhYvu3enaVr1hAeHm5orpLIaoRCCEO5XC6OHTtG\n7dq1z1slsSaTJ/IIIUSAkqfSCyFENSYFXAghApQUcCGECFBSwIUQZSooKCAtLY2ffvrJ6CjiHFLA\nhRClOnz4MJ3i4xnWvTtdW7XitpEj5WHHfkIKuBCiVPf+7W8MP3SIH7KyOJCbyw+JiXzwgczX8wdS\nwIUQpdqxYwejXC4UEAoMczrZsWWL0bEEUsCFEGWIi4tjnnsVvjxgkd1OfLvz1q0TBpCJPEKIUh04\ncIBrrriC8KwsThQW0s3h4LP58zGbffI4gRpLZmIKIbzC6XSSmppKaGgobdu2RalS64rwAingQggR\noGQqvRBCVGNSwIUQIkBJARdCiAAlBVwIIQKUFHAhhAhQ3niocT+l1G6l1F6l1KPeCCWEEKJsHg0j\nVEqZgD3A1cBhYBNwk9Z61zn7yDBCIYSoIF8MI+wG7NNaH9RaFwCfA0M9PKcQQohy8LSANwDSz9nO\ncL8mhBCiinm6mEG5+kYmTJhw9nuHw4HD4fDwskIIUb0kJSWRlJRUoWM87QPvDkzQWvdzbz8OFGmt\nXzpnH+kDF0KICvJFH3gy0FIp1VQpFQzcCCzw8JxCCCHKwaMuFK11oVLqPmAxYAKmnzsCRQghRNWR\n1QiFEMIPyWqEQghRjUkBF0KIACUFXAghApQUcCGECFBSwIUQIkBJARdCiAAlBVwIIQKUp2uh1Dib\nN29m2bJl1K5dm1GjRmG3242OJISooWQiTwV89dVX3HzzPRQUjCI4+AcaN/6FzZtXShEXQnhdeSby\nSAGvgNjYizl2bAZwBaCx24cwefJgxowZY3Q0IUQ1IzMxvezMmVNAvHtLkZ8fz8mTJ42MJISowaSA\nV0CfPtcSEvIQcAJYi8XyCX369DE6lhCihpICXgGffvouffvmYbO1IDp6FB999CZdu3Y1OpYQooaS\nPnAhhPBD0gcuhBDVmBRwIYQIUFLAhRAiQEkBF0KIACUFXAghAlSlC7hS6galVJpSyqWU6uTNUMK7\nsrOzmTp1KhMnTmTdunVGxxFCeEmlhxEqpVoBRcA7wINa6y0X2E+GERooOzubTp16kp7emLy8VoSE\nfMS0aZP4619vMjqaqGaWL1/O3PlzqR1Zm7H3jCU2NtboSAHNJ2uhKKWWIwXcb02bNo3x4+fhdCYC\nCthInTo3cPz4T0ZHE9XIzJkzufP+O8npmIP5jJmLMi5i+9btxMTEGB0tYMk4cMGvv/5KQUFzios3\nQHOysn4t9RitNe++O52ePQcxcOCNbNlS4v/NQpz1yNOPkDM0B3pC4YBCfq3/Kx9++KHRsaq9UtcD\nV0otBeqW8NYTWuvE8l5kwoQJZ793OBw4HI7yHio81LdvX/71r74UFAwFWhES8gh9+w4o9ZhJk17n\nmWfeJTt7InCEFSuuZdOmlbRu3donmYVvLFy4kFlfziIyIpKHHniIxo0bV/pcuc5cCPvfdoG9gDNZ\nZ7yQsuZISkoiKSmpYgdprT36ApYDnUp5XwtjJSYm6gYN4nVYWLQePnyUPn36dKn7N2zYWsNGDVqD\n1ko9rh999AkfpRW+MP396dpex67pjzZdYdK1Y2rrjIyMSp/vnvvv0bY4m+YeNDehbRE2nZyc7MXE\nNY+7dpZaf731RJ5S+2mEsQYNGkRGxqAKHKEovj/9OxdKWbycShjpXxP/hXOIExqDCxdnCs7w4Ycf\n8uSTT1bqfK+9+hoWi4Uv535JrfBaTP58Mp07d/ZyavFnngwjHK6USge6AwuVUou8F0sY6eGHx2K3\njwZmodQk7Pb3ue220UbHKrf8/HymTp3Kgw8+ypw5c37/JCjOkZ+XD7b/bbuCXeTm5Vb6fBaLhdde\nfY2MHzPYtW0X/fr180JKUZZKt8C11nOBuV7MIjyktUYpzz8MjRt3LxER4Xz88SwiIsJ45pllxMXF\neSFh1XO5XPTuPYitW4PIyelFaOizbNiwlVdeed7oaH7l1ltuZcrsKTgdTvgNbKk2rn/jeqNjiQqS\n5WSrgV9++YVhw24mOXklERExfPDBVIYMGWJ0LEMsX76cIUPGk5W1FTABJ7BYmnDy5FHCwsLKOrzG\ncLlc/N/z/8fnX35OrbBavPrCq1x55ZVGxxLnkGdiVmPffPMNCxcuJTY2irlzF7NjRw8KC58BtmCz\nDWPz5po5aiQxMZGbb57C6dPful8pIiSkDocO7ZYxySKglKeAe+smpvChKVPe5pFHXsLpHIvFkkZB\nwVqKBwOZgR4oNZjVq1dfsICfOnWKzz//HKfTyaBBg4iPjy9xv0B02WWXERR0FzAdcGCxTKFVq9ZE\nR0cbHU1UU/n5+Rw8eJCoqCiioqJ8e/Gyhql4+oUMI/S6iIi6Gra7h/kVabBrSHFvF+qwsG567ty5\nJR57/PhxXb9+C22z/UVbLPdqu72OXr16tY9/gqqVkpKiO3ToqevUaaKvvXaEPnbsmNGRRClcLpfe\ntm2b3rRpk87LyzM6ToXs3LlTxzaM1WExYTrEHqInPDfBa+emHMMIpQslAFmt4eTl/QjUAcBk6ktQ\n0DaUGoHJtIlOneqwYsU3mEym84596qlnePnlXygoeMf9yud06DCVrVtX+u4HEMItJyeHqwdcTcqu\nFIIsQdQNr8uapDUB84kprm0c+5rtQ3fRcAbsM+wsmrOIXr16eXxumUpfTQ0ffj1W6x3ADuALQkK2\n8fzzj2CxfIXLdYTk5HV8+OHHJR577NhJCgpanfNKKzIzT/oitvCykydP8tprr/Hiiy+yfft2o+NU\nyr9f/jdbMreQfVc2Z+44w8HaB7nvgfuMjlUuWmv2796P7uBuoNaCohZFpKSk+CyDFPAA9MEHU7n5\n5ibUr3897dr9l2++mcN//jOF7OzXyc//mby8Tdx//6Ps2rXrvGMHD74Gu/0NYDtwFJvtSQYPvtbn\nP4PwzIkTJ2jToQ2PffQYTy98mu49u/P9998bHavCUtJSyG2eW1yJFBTEFbBj5w6jY5WLUoq6jerC\nXvcLeWA6ZKJ58+Y+yyAFPABZrVbee+91Dh/eTWrqGtq0acNvv/0GjHTvEY/Z3JPU1NTzjh08eDAv\nvfQwkZH9sNlacf31DZk06QWf5heee3PKm2TWzSR3aC6uvi6c/Zzc/9D95+23atUqGjdvTLAtmK49\nupKenm5A2gvr3L4zth9sUAhoCN4ZTIdLOhgdq9zmzJxDraW1iPgsAvu7dv4y4C/079/fZ9eXPvBq\noLCwkIiIGJzOb4FuwCns9o4sXz6bbt26/WFfl8vF6tWryc7Opnv37lx00UWGZBaeuXfcvUzdNRWu\ncL9wFBp824CMHzPO7nP48GFatWtFVr8saAKmjSYu/uVi9uzY45UJX96Ql5fHoOGDWLNhDUHmIJo1\naMaK71YE1L/LzMxMUlJSiI6Opl27dl47rwwjrCJ5eXlMnPgS69enkJDQnOeee4rw8HDD8pjNZmbO\n/JC//nUgFksnCgrSGDPmlvOKd35+Pr17DyI19WeCgupiMu1k1aoltGnTxqDkorKGDR7Gh599iLOZ\nE8LAttzG0MFD/7DPxo0bCWoYBO5Roq6eLg69eogTJ074zU3CkJAQlixcwr59+ygsLCQuLq7Em+/+\nLCoqiquuusqQa0sLvIK01vTrdx2rVhWRkzOKkJBFtGy5iy1bVmGxGLvgU3p6Otu3b6dhw4Zccskl\n573/5ptv8sgjC8nJ+RowodTbdO48m02bAq/vVMD7H7zP408/Tm5uLiNGjODtN94mODj47PsrVqxg\n0F8HkXV7VnFT7TewvGXh9KnTWK1W44KLcpGZmFUgIyODli07kpt7GAgGNLVqtefbb9/m8ssvNzpe\nqcaPf4jXX48GHnW/spfo6Gs5cuQHpk2bRmrqbjp2bMPtt98ecK0gcb6ioiIGXzeYFakryK+fj+UH\nC88++iwPPfiQ0dFEOUgXShVwuVwoZeKP938tuFwuoyKV22WXdWH69H+Tnf13IBKLZSpdunRh+PBR\nfP/9MZzOQdjtM/j225XMmTPDb/pJazqXy8WE5ybw6axPCQ0N5ZWJr5Rrtb+goCAWfLmAOXPmcOjQ\nIbp16ybrnVQ3Zc308fSLajYTs6ioSF9+eV9ttY7SsFhbLON18+aX6JycHKOjlamoqEiPG/ewtlhC\ntdUapdu1665XrVql7faGGnLdMzmd2maL1fv27TM6rnB77MnHtP1iu+ZONCPR9ki7Xr9+vdGxRBWj\nHDMxZRhhBSmlWLz4K267LYbOnV9i5Mgc1q9fFhB9ikopXnvtZU6cOMKBAztISVmL1WrFbL4ICHHv\nZcVsjiQ7O9vIqOIcM2bOwHmNExoCrcDZ0cnsObONjiX8gHShVEJYWBhTp04yOkalhYeHnx0107Zt\nW8LDc8nOfh6XawQm0+dcdJGJVq1alXEW4StWqxWc/9s25ZoItYcaF0j4DWmB13BWq5XVq5fQs+d6\nYmMH43BsY/XqxX8YzSCM9cIzL2D/2g5rwbTURK19tRjz9zE+zbBy5UoaXdyIYFsw3a7oRkZGRtkH\nVbGioqKyd6rmZBSKEAFg2bJlzJozi/Ba4Yy7b5xHT5CvqIyMDFq3a01Wf/eEoA0mWpxowa6UXYbc\n6F6wYAF/u/NvnD55mg5dO7BgzgIaNGjg8xxVTYYRCiE8NmfOHO547g5OX3e6+AUNwf8J5shPR3y+\n/vWePXvoeGlHckbkQH0wrTbRJqsNKZt8t4CUr1T5aoRKqVeUUruUUilKqa+UUhGenE8I4X+ioqIo\nyiyC30fK/gbapQ15RN3atWsJahkEjQEzuHq52LFtB3l5eT7P4g887QNfArTRWrcHfgAe9zySEIHL\n5XIxa9YsXnnlFVatWmV0HK+48sor6dGxB6GfhmJZasH+qZ0Xn3+RkJCQsg/2spiYGNQx9b//TI5D\niDWkxt6z8VoXilJqODBCa33zn16XLhRRIxQVFdF/aH/W7FxDfr18LHssTHxqIg+Mf8DoaB5zuVx8\n8cUXpKenGzohyOVy0W9IP9bvWo8rxgV74O3X32b0LaMNyVOVfNoHrpRKBD7TWs/80+tSwN1ycnKw\nWq0yw7GaWrZsGcNuHUbWbVlgAn4tXnsk+0y24evkVCcul4v58+fz888/c9lll9GpUyevnLewsBCT\nyeQ3v59emUqvlFoK1C3hrSe01onufZ4E8v9cvH83YcKEs987HA4cDkdZl61W9u3bx4ABN7B/fxp2\nezgffzydYcOGln2gCCgnT54kKCqouHgDRIAKUmRlZVG7dm1Ds1UnJpOJ6667zmvnO3LkCIOuG8S2\nTdsIDQ9l2tvTuPHGG712/vJKSkoiKSmpQsd43AJXSt0K/B3oo7XOLeH9Gt0C11rTvPklHDx4B1qP\nBzZhsw0kJWUtLVu2NDqe8KL09HRaX9Ka7AHZ0Lh4uF1cZhxp29L8plUnztf5ss6khKTg6uWCo2Cb\nbWP9ivUlrujpS74YhdIPeBgYWlLxFvDbb7+RkXEQrf8BKKAbZnNvkpOTjY4mvKxRo0YsnLeQBmsb\nYPmvhU75nViycIkUbz/mcrnYtnFbcfE2AfWBeFizZo3R0crF06n0b1C8pupS9z/SdVrrsR6nqkZq\n1aqFyaQoKNgJJAA5FBVtp149+WPyNpfLxfHjx6lTpw5mszGrRFx55ZVkHDB+lqIoH5PJRFhkGKd/\nOQ0NABeYjpmIjY01Olq5eNQC11q31Fo30Vp3dH/V2Kp07NgxBg78C7GxzenWrQ87d+4Eiv+BvPPO\nVOz2qwiRLklOAAAPdUlEQVQNHU1YWFcGDrxUlvX0svXr1xMd3ZhmzdoRGRnLwoULjY4kAsT0d6Zj\nm23D/o2dsE/C6BrXlaFDA+MeVY2ciZmfn8+DDz7J3LkLiYiI4PXXJ9KnT59Kn09rTfv2l7N792UU\nFNyNUsuoXft59u3bfvbm1fbt20lOTqZBgwb07dtXPlZ7UW5uLnXrNuO3394BhgDrCQ0dzN69qdSr\nV8/oeCIA7Nixg7Vr1xITE8PgwYP94oEmMpX+AsaMGccnn+whJ+cV4Efs9r+zdu13tG/fvlLnO3Lk\nCM2bdyA39yiQDWjCw4cxa9bD5Vp4X3hmz549dOkykKysfWdfi4hwMGfOU1x99dUGJhOi8qr8Jmag\nmjXrC3Jy3gEuAYaRl3criYlfV/p8drsdl8sJ3ALEALFkZ++Qsb8+EhsbS0HBcWCv+5Vj5OfvolGj\nRkbGEqLK1cgCbrPZgWNnt83mo4SG2it9vsjISDp16gIcBI4DJ9G6A99+Kw8L9oXIyEhef30SNtsV\nhIcPw27vxEMP3U98fLzR0WqczMxMvvjiC+bNm4fT6Sz7AOGRGtmF8sEHH3HffU/hdI7DbP6RqKjF\n7NixkTp16lT6nFddNZzly0cB17tfWUS3bpPZsGGJVzKLsu3evZu0tDSaN29Ohw4djI5T4+zbt4/u\nV3QnPzofna+JNcWSvDaZyMhIo6MFJOlCuYDbbvsbc+dOZ+zYIzz+eCypqes9Kt4ALVo0wmL53+JF\nZvNqmjVr6GlUUQGtWrVixIgRUrwNMvYfYznV/hRnrj9D1k1ZpIem8/yLz1foHO+++y51G9floroX\n8Y+H/kFhYWEVpa0eamQLvCqcOHGCLl16cfJkDGAmLOwgyckrqV+/vtHRhPCJVh1asafznuKlXgG2\nwHX26/jysy/LdXxiYiIj7xiJc5gTrGD/xs4Dox5g4rMTqy60H5MW+AWcOnWKO+64j0svvYaxY//J\nmTNnPD5nnTp1SEvbxKefPsgnn4xj9+4tUryF33K5XGXvVEG9e/bGmmyFQiAH7Nvt9OlZ/uG5X8z9\nAmdXZ/GEmihwOpzM/koe3lyaGlfACwoKuOKKa/nkk0I2bnyA99/PpE+fIV55vl5oaCiDBw9myJAh\nZx8aLIQ/+f7774muF40l2EJ8u3j27t1b9kHl9OpLr+Jo4sD8shnTJBOjB4zm7rvvLvfxF0VehOm3\nc8Zfn0L6z8tQ47pQkpOTueqqv3HmzA6K1yZxYbc3Y+vW74iLizM6nhBnFRUVkZycTHZ2Np07d/a4\nUXD48GHi28aTPTgbmoJKVjTc05CDew8SFOS9tlx2djZms7nCD3zIyMigfZf2nGl8BleIC+t2K98m\nfkvPnj29li2QeGU52eqm+A/FBWiKC7gGimRmpPArBQUFXDPoGjZt34QpzERIdghrktZ4tILl5s2b\nMTU0QfPibX2p5sS6E/zyyy9e7e4LDQ2t1HENGzYkbVsaM2bMIDc3l+veuo62bdt6LVd1VOMKePv2\n7WnRog67dt1GXt4QrNZZdOrUhhYtWvgsg9YarbVXWz2ienn33XfZcGgDOXfmgAmC1gdx65hbWbO8\n8qvkxcTE4DrugnyKl6A7Ba48l191U9StW5dHHnnE6BgBo8ZVELPZzMqVi7jrrlh69/6Y+++PY8mS\neT5pgWuteeyxf2G1hhESYuemm24nPz+/yq8rAs/OPTvJaZRz9uEQRS2K2L9/v0fnvPTSSxlyzRBC\nZ4RiW2TD/omdl196Gbu98pPYhLFqXB+4kaZP/4Dx418nO3shEIbNdhP33NOeV199wehows989NFH\njJ0wFudIJ4SA5XsLV0VexbcLvvXovFprFi1axE8//USnTp249NJLvZRYeJssZuVnhgy5icTEGKA/\ncCWQTOvWj7Bz5zqDkwl/o7Xmjrvv4NNPP8VsNdO4fmOSliYFzDrVwnNSwP1IRkYG8fGdcDpbUHzz\nNBu4iT591vHdd/MMTif81bFjx8jOzqZx48Z+scSpUTZt2sTevXtJSEioMTNtZSKPHxk//gny8u4C\n1gKrgU5YLC/y+uvSfSIuLCYmhmbNmlWqeGuteePNN2jZriWtO7Rm5swSnznu956e8DSOAQ7u/s/d\n9LiqB5MmTzI6kt+QFriPdOhwJSkpzwBXuV+ZicPxOcuXLzAylqjG3nn3Hf454Z84+zmhEOyL7Hw2\n/TOGDBlS5rGFhYWYTCbDh9fu37+ftp3akjsmF8KA3yDk3RAO/XiImJgYQ7NVtSptgSulnlNKpSil\ntimllimlZPHlUvTs2Q2r9S2Kx3BlY7dP59prexkdS1Rj02ZMw9nbCU2A5uC83Mn7n7xf6jFHjx6l\n2xXdCA4JplbtWsz4eIZvwl7AkSNHCIkOKS7eABEQHBnM0aNHDc3lLzzpQnlZa91ea90BmAc846VM\n1dLLL/8fvXoVYrFEYbHEMGRIEx566B9GxxLVmM1mg5z/bascRai99Ek2w28czla1Ff2kJvuv2dwz\n/h6Sk5OrOOmFJSQkUJhZCD+6X9gNQXlBXHzxxYZl8ieVnsijtT53Bagw4ITncaovm83G4sVzOXXq\nFCaTSdZKEVXuuaeeY8CwAThPO8EF9m12Hlv5WKnHbFy7EdfDruLx57HgauVi9erVdOnSxTeh/yQq\nKorErxIZfsNwnE4n4RHhfJ34daVne1Y3Hs3EVEo9T/FzxJxAd68kquZ+f8ixEFXtyiuvZPni5bz/\n0fuYzWbGThlLQkJCqcdERkWS+XNmcbdLEZiPm6lbt65Xc2VmZvLEv57gh/0/0LN7T5564imCg4Mv\nuH/v3r05dfwUv/76K5GRkYb3y/uTUm9iKqWWAiX97T2htU48Z7/HgHit9W0lnMNrNzFPnz5NXl4e\nderUkb9EIapAYmIiN95yIypeEXQiiA5NOrB8yXLMZu+supGTk0PbTm1Jj0ynoHEBth02rm51NQu+\nlJv5f+bxYlZa677lvNZM4JsLvTlhwoSz3zscDhwORzlPW6yoqIi77/4HH344naCgYDp27MyiRXP8\nag0HIaqDwYMHs3ndZlavXk2dOnUYPHiw14o3wOrVqzlecJyCawtAQU7LHBZPXkxmZiZRUVFeu04g\nSkpKIikpqULHVHoYoVKqpdZ6r/v7+4FuWutbStjP4xb4tGnTGT/+XZzOJUAYwcF3c911RXz22XSP\nziuE8K0lS5Zww703cPrm08UvFELw5GAOHzzs8WMNq5uqXk72RaVUPOAC9gP3eHCuUq1atQmnczQQ\nAUB+/j2sX39eb40Qws/16NGDSB2J8zsnhY0LsaZa6eXoJcW7kio9jFBrfb3Wup3WuoPWeoTW+pg3\ng50rLq4JVmsSUPzUnKCgJJo2bVzqMUII/xMaGsqmtZsY2WIkl2Zcyr0D7mX+nPlGxwpYATET0+l0\n0qPHNezbl0NQ0EVYLHtYt26ZR4vbCyGEP6tWi1kVFBSwcuVKcnJy6NGjhwzHE0JUa9WqgAshRE0i\nqxEKIUQ1JgVcCCEClBRwIYQIUFLAhRAiQEkBF0KIACUFXAghApQUcCGECFBSwIUQhtJac+jQIX78\n8UeKioqMjhNQpIALIQyTl5dH/yH9ib8knrZd29K9Z3dOnz5tdKyAIQVcCGGYf7/8b1b+uJLc+3PJ\nuS+H1PxU/vnoP42OFTC8t1K7EEJU0IbNG8hplXO2EuUl5JG8xbiHKAcaaYELIQzTplUbQg6EFK8U\nrcGy30Lr+NZGxwoYspiVEMIwWVlZ9LyqJ/uO7EOZFdG2aNatWEdMTIzR0QwnqxEKIfxeYWEhycnJ\nuFwuOnfujNVqNTqSX5ACLoQQAUqWkxVCiGrM4wKulHpQKVWklLrIG4GEEEKUj0cFXCnVCOgL/OSd\nOMZJSkoyOkK5SE7vCoScgZARJKcRPG2BTwIe8UYQowXKX6rk9K5AyBkIGUFyGqHSBVwpNRTI0Fqn\nejGPEEKIcip1JqZSailQt4S3ngQeB645d3cv5hJCCFGGSg0jVEq1BZYBTvdLDYHDQDet9bE/7Stj\nCIUQohJ8Mg5cKXUA6Ky1PunxyYQQQpSLt8aBSytbCCF8rMpnYgohhKgaPp2J6e+TfpRSzymlUpRS\n25RSy9zj3P2OUuoVpdQud9avlFIRRmf6M6XUDUqpNKWUSynVyeg8f6aU6qeU2q2U2quUetToPCVR\nSr2vlDqqlNpudJbSKKUaKaWWu/++dyilxhmdqSRKKatSaoP793unUupFozNdiFLKpJTaqpRKLG0/\nnxXwAJn087LWur3WugMwD3jG6EAXsARoo7VuD/xA8Yggf7MdGA6sNDrInymlTMCbQD8gAbhJKeWP\na5h+QHFGf1cAPKC1bgN0B+71xz9PrXUu0Nv9+30J0FspdYXBsS5kPLCTMrqnfdkC9/tJP1rrM+ds\nhgEnjMpSGq31Uq317w8P3EDxKCC/orXerbX+wegcF9AN2Ke1Pqi1LgA+B4YanOk8WutVwCmjc5RF\na/2L1nqb+/ssYBdQ39hUJdNa/z5yLhgwAX438EIp1RAYAEyjjOHZPinggTTpRyn1vFLqEPA34N9G\n5ymH24FvjA4RYBoA6edsZ7hfEx5SSjUFOlLcsPA7SqkgpdQ24CiwXGu90+hMJZgMPEzxYy5K5bVH\nqgXKpJ9Scj6htU7UWj8JPKmUeoziP8jbfBrQrayc7n2eBPK11jN9Gs6tPBn9lNy5rwJKqTBgDjDe\n3RL3O+5Prh3c940WK6UcWuskg2OdpZQaBBzTWm9VSjnK2t9rBVxr3fcCgdoCzYAUpRQUf9zfrJQ6\nb9KPL1woZwlmYmDLtqycSqlbKf6Y1ccngUpQgT9Lf3MYOPcGdSOKW+GikpRSFuBL4BOt9Tyj85RF\na/2bUmoh0AVIMjjOuS4HhiilBgBWIFwpNUNrPbqknau8C0VrvUNrHau1bqa1bkbxL0onI4p3WZRS\nLc/ZHApsNSpLaZRS/Sj+iDXUfWPG3/nbMgvJQEulVFOlVDBwI7DA4EwBSxW3zKYDO7XW/zU6z4Uo\npeoopSLd39soHlThV7/jWusntNaN3LVyJPD9hYo3GPNAB3/++PqiUmq7u4/MATxocJ4LeYPim6xL\n3UONphod6M+UUsOVUukUj0pYqJRaZHSm32mtC4H7gMUU3+mfpbXeZWyq8ymlPgPWAnFKqXSllCHd\neeXQA7iZ4lEdW91f/jh6ph7wvfv3ewOQqLVeZnCmspRaL2UijxBCBCh5pJoQQgQoKeBCCBGgpIAL\nIUSAkgIuhBABSgq4EEIEKCngQggRoKSACyFEgJICLoQQAer/AWHvM/IZA+dIAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d73fcd0>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Sampling Code ###"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "table = np.zeros(n)    #currently assigned table number for each data\n",
      "\n",
      "alpha = 0.5 #hyperparameter (concentration parameter)\n",
      "\n",
      "#set table of first data point\n",
      "table[0] = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Sampling table number of second data point"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "candidate = list()  # currently allocated topic list\n",
      "prob = list()       # probability of being assigned to each table\n",
      "\n",
      "\n",
      "# sampling table number of second data point\n",
      "for _t in np.unique(table):\n",
      "    if _t != 0:\n",
      "        candidate.append(int(_t))\n",
      "        # this is not correct marginalized probability, see http://en.wikipedia.org/wiki/Conjugate_prior for correct posterior analysis\n",
      "        mean = np.mean(x[table == _t, :], axis=0)\n",
      "        prob.append( np.sum(table == _t)/(np.sum(table != 0)+alpha) \\\n",
      "                    * multivariate_normal.pdf(x[1,:],mean, np.identity(2)) )\n",
      "\n",
      "# get new table number\n",
      "def get_new_table_no():\n",
      "    u_t = set(np.unique(table))\n",
      "    for i in xrange(1,1000):\n",
      "        if not i in u_t:\n",
      "            return i\n",
      "\n",
      "# compute probability of being assigned to new table\n",
      "new_table = get_new_table_no()\n",
      "candidate.append(new_table)\n",
      "prob.append( alpha/(np.sum(table != 0)+alpha) \\\n",
      "            * multivariate_normal.pdf(x[1,:], np.zeros(2), np.identity(2))  )\n",
      "\n",
      "# normalize probability\n",
      "prob = np.array(prob)\n",
      "prob /= np.sum(prob)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Candidate tables and corresponding probabilities"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print candidate\n",
      "print prob"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[1, 2]\n",
        "[ 0.99794869  0.00205131]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# draw a random sample from multinomial distribution parameterized by prob\n",
      "sample = np.random.multinomial(1, prob)\n",
      "# new table number for second data point\n",
      "table_idx = sample.argmax()\n",
      "table[1] = candidate[table_idx]\n",
      "\n",
      "# table assignment of first two data points\n",
      "table"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "array([ 1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "max_iter=100\n",
      "number_of_clusters = list()\n",
      "\n",
      "for s_iter in xrange(max_iter):\n",
      "    for idx in xrange(n):\n",
      "        candidate = list()\n",
      "        prob = list()\n",
      "        \n",
      "        table[idx] = 0\n",
      "\n",
      "        for _t in np.unique(table):\n",
      "            if _t != 0:\n",
      "                candidate.append(int(_t))\n",
      "                # this is not correct marginalized probability, see http://en.wikipedia.org/wiki/Conjugate_prior for correct posterior analysis\n",
      "                mean = np.mean(x[table == _t, :], axis=0) #mean of currently assigned data\n",
      "                prob.append( np.sum(table == _t)/(np.sum(table != 0)+alpha) \\\n",
      "                            * multivariate_normal.pdf(x[idx,:], mean, np.identity(2)) )\n",
      "\n",
      "        # compute probability of being assigned to new table\n",
      "        new_table = get_new_table_no()\n",
      "        candidate.append(new_table)\n",
      "        prob.append( alpha/(np.sum(table != 0)+alpha) \\\n",
      "                    * multivariate_normal.pdf(x[idx,:], np.zeros(2), np.identity(2)) )    \n",
      "\n",
      "        # normalize probability\n",
      "        prob = np.array(prob)\n",
      "        prob /= np.sum(prob)\n",
      "\n",
      "        #sample table index\n",
      "        sample = np.random.multinomial(1, prob)\n",
      "        table_idx = sample.argmax()\n",
      "        table[idx] = candidate[table_idx]\n",
      "        \n",
      "    #for tracking purpose\n",
      "    number_of_clusters.append(len(np.unique(table)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Plot Results###\n",
      "Plot assignment of final samples"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "c_mask = [color[int(table[i])%len(color)] for i in xrange(n)]\n",
      "plt.scatter(x[:,0],x[:,1], c=c_mask)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "<matplotlib.collections.PathCollection at 0x10dffda90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEACAYAAACqOy3+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0VNX6//H3zqRMJiEhhITeJUDoPQjIQETpRfR7vcq1\nYLugAl4sIHoNPzsKVkBFBBFBBaVEQEQgQGgSOqEICJKAEpoGMqmT/fsjkYsS0mYyZyZ5XmtlrZkz\np3xA87Bnn733UVprhBBCeB4vowMIIYQoHSngQgjhoaSACyGEh5ICLoQQHkoKuBBCeCgp4EII4aGc\nUsCVUial1C6lVKwzzieEEKJozmqBjwEOADKoXAghXMThAq6Uqg30Az4GlMOJhBBCFIszWuBvAU8B\nuU44lxBCiGJyqIArpQYAKVrrXUjrWwghXEo5shaKUuoV4F9ADmAGgoCvtdb3XLWP9IsLIUQpaK0L\nbRg71ALXWj+rta6jtW4A3Amsvbp4X7Wf2/+88MILhmeQnJLTUzNKTuf/FIezx4FLa1sIIVzE21kn\n0lqvB9Y763xCCCEKJzMx81mtVqMjFIvkdC5PyOkJGUFyGsGhm5jFuoBSuqyvIYQQ5Y1SCl2WNzGF\nEEIYRwq4EEJ4KCngQgjhoaSACyGEh5ICLoQQHkoKuBBCeCgp4EII4aGkgAshhIeSAi6EEB5KCrgQ\nQngoKeBCCOGhpIALIYSHkgIuhBAeSgq4EEJ4KCngQgjhoaSACyGEh5ICLoQQHsqhAq6UMiultiml\ndiulDiilXnVWMCGEEIVz6KHGWusMpVRPrbVNKeUNxCulummt452UTwghxHU43IWitbblv/QFTMAF\nR88phBCiaA4XcKWUl1JqN3AGWKe1PuB4LCGEEEVxRgs8V2vdBqgN3KSUsjqcSgghRJEc6gO/mtb6\nD6XUcqADEHf1ZzExMVdeW61WrFarsy4rhBDlQlxcHHFxcSU6RmmtS31BpVRVIEdr/btSyh9YBUzS\nWq+5ah/tyDWEEKIiUkqhtVaF7eNoC7wG8KlSyou87pjPri7eQgghyo5DLfBiXUBa4EKUW7m5uXz2\n2Wfs3LmPyMjGPPDAA3h7O61ntkIrTgtcCrgQotT+9a+H+eabvdhsQ7FYvsNqrcq3336FUoXWHVEM\nUsCFEGUmKSmJiIi2ZGScAAKBTCyWJmzZsoxWrVoZnM7zuaIPXAhRQWit+fDDmaxZs5l69Wpwxx1D\nMJmCgYD8Pfzw9g7j8uXLRsasUKQFLoQoltGjn+KTT9aTlvYIvr7bqFlzE6BIShqG3T4cL69lhIVN\n49ixfQQEBBR5PlE46UIRQjhFdnY2FkslcnJOA1UATWDgzUydeidffPEte/fuJSKiCXPnTqNRo0ZG\nxy0XpAtFiHIiNzeXs2fPEhwcjNlsdvn17XY7eQ0xS/4WBQRRqVIl1qxZyqFDh/j5558xmUwuz1aR\nyXrgQri5o0eP0qBBC+rXb05wcFWmT//Q5RnMZjN9+gzGbL4b2IyX11R8fH4kOjqaV155k3btrNx1\n1ztERnbks88+d3m+ikq6UIRwc02atOfIkXvQegxwDIvlJtavX0qHDh1cmiM9PZ0nn3yONWviqVOn\nJtOmvY5Sitatu5KevhuoCRzAbL6RM2dOEhQU5NJ85Y10oQjh4XJycjhyZDda/5i/pRHQl4SEBJcX\ncH9/f6ZNm/KXbT/88AO+vs1JT6+ZvyUSb+9Qfv31VyngLiBdKEK4MW9vb0JCagAb87ek4+X1I/Xq\n1TMy1hXNmjUjO3sfsDN/yyq8vC5Tt25dI2NVGFLAhXBzCxZ8gsVyB0FBgwgIaEX//h3o06eP0bEA\nqFWrFnPnfoS/fzQBAXUJCrqH2NiF+Pv7Gx2tQpA+cCE8QFJSEjt27CA8PJwuXbq43VR1m83GmTNn\nqFmzJn5+fkbHKRdkHLgQQnio4hRw6UIRQggPJQVcCCE8lBRwIYTwUFLAhRDCQ0kBF0IIDyUFXAgh\nPJTDBVwpVUcptU4plaiU2q+UGu2MYEIIIQrnjBZ4NvCE1ro5EAU8qpRq5oTzCiEMsGzZMho2bE1Y\nWH0eemg0mZmZRR6TmprKiy++zMMPP85XX32FzP1wDYcLuNb6N6317vzXl4GD5C1LJoTwMD/++CN3\n3vkQx4+/yblzq/n882OMGjWu0GNsNhsdOvTg5ZcPMHNmI+6/P4aYmJddlLhwc+fOo1Wr7rRufRPz\n5y8wOo7zaa2d9gPUB34BAq/apoUQnuG55/6rlXpWg87/OaZDQmoXesxXX32lAwN7acjNPyZZ+/j4\na7vd7qLUBfviiy+1xVJfw0oNK7TFUlcvWvS1oZlKIr92FlpznbacrFIqEFgEjNF5LfErYmJirry2\nWq1YrVZnXVYI4UTBwZXw8dlLVtafW5IICAgs9BibzYbWYeQ9pQcglNxcO3a7HS+v63/Jt9vtPP/8\ni3z22UIsFguvvz6RIUOGOOOPAcD06fOw2V4H+uTnfJUPPvicYcNuc9o1nCkuLo64uLiSHVRUhS/O\nD+ADrALGFvCZC/6tEkI4w/nz53XNmjdoX997tVIx2mKprr/6amGhxyQlJelKlcI0zNKwR/v5/VPf\neuttRV7r2WdjtMVyo4adGlZqi6W63rBhg7P+KPqWW4ZpmHnVt4npesCAO512/rJGMVrgDi9mpfKW\nRfsUOK+1fqKAz7Wj1xBCuM758+f5+OOPuXgxlYED+9G1a9cij9m9ezcPPzyO06d/JTr6JqZNe5PA\nwMJb7rVqNeX06ZeA/oA/8BqPPXaW996bcs2+q1ev5tChQ0RGRhIdHV2sP0d8fDy33DKU9PRngFws\nljdYsyaWqKioYh1vNJesRqiU6gZsAPYCf55sgtb6u/zPpYALUYDVq1fz/vtz8PY2MW7cv7nxxhuN\njuQy77//AaNH/wetawA2YDEm05c8/XQgr7zy4l/2HTPmGWbNWkJOTjTe3j/wyCO3M2XKK8W6zrZt\n25gxYw5eXopRo0a4/ClGjpDlZIVwUytWrOCOOx7AZnsRyMTfP4YfflhaIYp4YmIiHTtGk56+BWgA\nxAL3ULmyP3v3bqNOnTpX9j1+/DiRkZ3IyPgJCAEuYDZHcPjwznL/1B9ZTlYIN/Xqq9Ow2d4CHgQe\nJT39BaZM+cDoWC6RmJiIt3dX8oo3wEBMJjtxcSv/UrwBzp07h69vHfKKN0AVfH1rce7cORcmdl9S\nwIUwQE6OHfC9aosvdnuuUXEKlJubS0pKCjk5OU497w033IDdvg1Iyd8Sj7+/Ly1atLhm32bNmmEy\npQDzgSxgHt7eF2jSpIlTM3kqKeBCGOCJJx7AYnkC+BqYj7//f3n88fvL9Jrr1q3jppsG0LHjzcyc\nOavQ2ZI7duygWrX61K3bjODgMJYsWeq0HO3atWPcuH/j79+C4ODuBAQMZdGizzGZTNfsGxgYyJo1\nsdSv/wpeXhYaNHidtWu/JSAgwGl5PJn0gQthkIULF/HWW7Pw9jYxYcKj9O3bt8yutW3bNnr1GojN\nNhUIwWIZxxtvPMGoUY9cs292djY1ajTk/PmpwB3AdiyWfhw6tPOaLg5HHD16lOTkZCIjIwkPDy9y\nf6212z0LtCzJTUwhBAAPP/w4M2fWAZ7O3xJHkybPcOjQtmv2PXHiBM2bd8NmS76yLTj4VubPH0O/\nfv3KNGdGRgbHjx8nPDyc0NDQMr2Wu5ObmEIIAEwmE0plXLUlo8AuC4CwsDDs9lTgUP6WC2Rn73dq\n67sgO3bsoFatG+jUaTC1ajVk8uS3yvR65YEUcCEqgEcffRCL5T3gDWAWFstDPP/8mAL3DQgI4IMP\n3sdi6UGlSkOxWNowatT9tGzZssTXPX/+PImJidhstkL301rTv/8dXLgwlcuXfyIzM5FJk6awY8eO\nEl+zIpEuFCEqiD179jB58nvYbJk8+OCd9O/fv9D9Dx8+zN69e2nQoEGpJsC8/fb7jB8/EV/f6phM\nl1i1agmdOnUqcN+0tDSCg0Ox29P5c02VgIDhvP/+zdx3330lvnZ5IH3gQghD7N27l6ioW0lP3wrU\nAxYTGjqas2dPFngjUmtNaGgtLl78FOgNXMBiac/3388r1lT+8kj6wIUQhsibrNOdvOINMJRLl/7g\n999/L3B/pRTffDOfwMC7CQ7uhr9/JKNG3V1hi3dxOW05WSGEe7pw4QL33DOSTZviqVatJrNnv0uX\nLl3K9Jp5k3W2AmeBMCAOs9lMcHDwdY+xWq0cP36AxMREatSoQURERJlmLA+kC0WIcq5r11tISLiB\nrKzxwDYCAx8lMTGhzNcSefbZSbz11jT8/CLIyTnM0qVfFHslQSF94EJUeOnp6VSqVBm7PY0/v3AH\nBt7JjBkDGD58eJlf/8iRI5w+fZrIyEjCwsLK/HrlSXEKuHShCFGO+fr64uVlwm4/DdQFcoGTVKpU\nqcD9t2zZwtatW6lVqxbDhg277ljx4mrcuDGNGzd26Bzi+uQmphDlmMlk4qWXXsJi6QlMwt9/EBER\n3gVO25827QNuvvkOxo8/zogRU+nX73Zyc91rgS3xV9KFIkQFsHLlSjZu3ETt2jUZMWIEZrP5L5/n\n5OQQEBBMVtZeoBGQTWBge775Zgq9e/c2JHNFJ10oQggA+vbtW+hiWWlpaeQ1thvmb/FBqSay7rab\nky4UIQTBwcE0bhyJyTQJSANWY7fHlflwQ+EYKeBCCAC+/34xbduux9s7lGrV/s3SpQuoX7++0bFE\nIZzxUONPyHusdIrW+prVbqQPXAghSs5VU+lnA32ccB4hhBAl4HAB11pvBC46IYsQQogSkD5wIYTw\nUC4ZRhgTE3PltdVqxWq1uuKyQggPcOLECXbv3k3t2rVLte54WdBaM3fuXLZv30OzZo14+OGH8fHx\nKdNrxsXFERcXV6JjnDKRRylVH4iVm5hCiJJYsmQpd9/9IN7eUeTk7OOee4YyY4bxj1IbMeJRvvxy\nOzbbHfj7ryYqyo8ffliKl5frOi1ctpiVFHAhREnZ7XaCgqpis30PdARSCQhoy6pVcw1dB/y3336j\nXr2mZGWdBIKAbAICmrN27bzrPlGoLLhkFIpSagGwGYhQSiUppe539JxCiPIvNTWVnBw7ecUbIAgv\nr/acPHnSyFjYbDa8vSsBfy745YPJFEZaWpqRsQrkjFEo/9Ra19Ra+2mt62itZzsjmBCifKtcuTKh\noeHAp/lbDpKTs57WrVsbGYt69epRt241vL3HA4fx8noLP79TtG/f3tBcBZFRKEIIQyilWLVqMdWq\nTcJsDsfPrzMzZrxJZGSkoblMJhPr1n1LdPRRqlXrT1TUd2zatJqgoCBDcxVEViMUQhjKbreTkpJC\nSEjINaskVmTyRB4hhPBQ8lR6IYQox6SACyGEh5ICLoQQHkoKuBCiSNnZ2SQmJvLLL78YHUVcRQq4\nEKJQp06dokmTdkRFDaFp047ceef98rBjNyEFXAhRqHvvfZSTJ4dy+fJPZGQcJzb2J2bPlvl67kAK\nuBCiUPv378duvxtQQAA22xB27txvdCyBFHAhRBEiIiLw8lqS/y4Ti2UlLVs2MTSTyCMTeYQQhTp+\n/Djdut3C5ctB5OScw2rtxNKlC/D2dsnjBCosmYkphHAKm83G3r17CQgIoEWLFihVaF0RTiAFXAgh\nPJRMpRdCiHJMCrgQQngoKeBCCOGhpIALIYSHkgIuhBAeyhkPNe6jlDqklDqilHrGGaGEEEIUzaFh\nhEopE3AYuBk4BWwH/qm1PnjVPjKMUAghSsgVwwg7AUe11ie01tnAF8BgB88phBCiGBwt4LWApKve\nJ+dvE0IIUcYcXcygWH0jMTExV15brVasVquDlxVCiPIlLi6OuLi4Eh3jaB94FBCjte6T/34CkKu1\nfv2qfaQPXAghSsgVfeAJQGOlVH2llC/wD2CZg+cUQghRDA51oWitc5RSjwGrABMw6+oRKEIIIcqO\nrEYohBBuSFYjFEKIckwKuBBCeCgp4EII4aGkgAshhIeSAi6EEB5KCrgQQngoKeBCCOGhHF0LpcLZ\nsWMHa9asISQkhLvvvhuLxWJ0JCFEBSUTeUrgm2++YfgDw8luno3v777UNdVlx5YdUsSFEE4nE3mc\nbOSYkaQPTSendw62222ctJ9k3rx5RscSQlRQUsBL4NIfl6Bq/hsFWSFZXLhwwdBMQoiKSwp4CUTf\nHI3fGj9IA06Cz34foqOjjY4lhKigpICXwOdzPqd3vd74z/AnbFUYn370KR07djQ6lhCigpKbmEII\n4YbkJqYQQpRjUsCFEMJDSQEXQggPJQVcCCE8lBRwIYTwUKUu4EqpO5RSiUopu1KqnTNDCedKS0tj\n+vTpvPTSS2zZssXoOEIIJyn1MEKlVFMgF/gQGKe13nmd/WQYoYHS0tJoF9WOJJLIDMnEb78fH7//\nMXfddZfR0UQ5s27dOpYtXkZwSDAjR42kWrVqRkfyaGU6jFBrfUhr/VNpjxeusWDBApJ1MunD0smN\nziV9WDpjnhxjdCxRzsyfP587+99J2ntp7H5lNx1bdSQlJcXoWOWe9IGXc7///jvZQdnw57/jVeBy\n6uVCj9Fa89HMj+h+c3f6D+3Pzp0FfrkS4oqYp2OYmD6Ru7iLsTljafV7K+bMmWN0rHKv0PXAlVKr\ngeoFfPSs1jq2uBeJiYm58tpqtWK1Wot7qHBQ7969+e9L/yW7cTZUBb91fvS+tXehx0x9eyovvPkC\nad3T4BKsj17P9s3badasmYtSC1dYvnw5i79cTFDlIMY+OZa6deuW+ly2DBtVqHLlfUh2CJcvFd5Q\nEH8VFxdHXFxcyQ7SWjv0A6wD2hXyuRbGio2N1bUa1NKBIYF66P8N1ampqYXuX7thbc1DaGLyflR3\npZ8Z/4yL0gpX+GTWJ7qGpYYezWj9T9M/dfWQ6jo5ObnU5xszcozu4N9Bz2KWfpmXdRX/KjohIcGJ\niSue/NpZaP111hN5Cu1oF8YaMGAAyQOSi3+AAq6+76zzbqiI8uPV/77Ks7ZnaUELsEP6pXTmzJnD\nxIkTS3W+N955gwk+E3jt69eoVKkSn7/1Oe3bt3dyavF3jgwjHKqUSgKigOVKqZXOiyWM9NSYp7B8\na4H9oLYoLPst3H/f/UbHKrasrCymT5/OuKfGsWjRoj+/CYqrZGZlUolKV94H2APIzMgs9fl8fHx4\n8503+Sn5J3Yc3EGfPn2cEVMUodQtcK31YmCxE7MIB2mtndJSHv34aIKDg/nsy88IrhzMC+tfICIi\nwgkJy57dbqfnrT3Z9dsu0mumE/BFANu2b+ON198wOppbGX7fcKZOm8pDtodIIYWV/itZe/tao2OJ\nEpLlZMuB3377jSH/N4SELQkEVwlm9szZDBo0yOhYhli3bh2D7h3E5fsv532/TAOf93y4cPYCgYGB\nRsdzG3a7nVf+3yss/mIxgZUCeXHKi/To0cPoWOIqxRkHLgXcQ61YsYLl3y2nWlg1Fn+7mP1++8m5\nKQd+Bf+v/dmxZUeFHDUSGxvL8AnDSb0jNW9DLvhN9ePksZOEh4cbG06IEihOAXfWTUzhQtOmT+Pp\nmKextbHhs8mH7F3Z8BxgAuqCaqKIj4+/bgG/ePEiX3zxBTabjQEDBtCkSROX5i9LXbp0wSvFC3YC\n9cFnhw9NmzYlLCzM6GiinMrKyuLEiROEhoYSGhrq0mvLRB4PNPGFidhut0FXyB6YnffP8Nn8D3PB\n66zXdQvWuXPnaNG2BeM+GseERRNo17kdmzZtcln2sla1alXW/7CeNmfaUHVhVXpV6cXqFatlFI0b\ny83NZc+ePSQkJJCVlWV0nBI5ePAgEXUjiG4fTf1a9Xl50ssuvb50oXggc4CZzEczISDvvWmeCa8z\nXqhIhelXE+0atWP96vWYTKZrjn3u+eeYvGoy2f2z8zbsgzZJbdi1bZcL/wRC5ElPT2fAzQP4ac9P\n+Hn5YaluYc2mNR7zjal1RGuij0YzSA/iAhcYbRnN/JXzuemmmxw+tzxSrZwaettQzCvMcAZIBL+z\nfrz8/Mv4HPbB/oedhG0JzPl0ToHHppxLIbtK9v82VIXzF867JLdwrgsXLvDOO+/w6quvsm/fPqPj\nlMrk1yZj32lnTtocZl6aSbMTzRj32DijYxWL1prEY4n01X0BqEIVOuZ2ZM+ePS7LIAXcA82eOZvh\nNw2n5nc1aflzS1YsXcGbb79JWnQaWU9kkTkik8fHPc7BgwevOXZgv4FYdlnyiv9l8N/oz8B+A13/\nhxAOOXfuHO2bt2fF+BXsfn431igra9d63jDAQ3sOEZURhQkTCkXX7K4c2n/I6FjFopSiXvV6bGMb\nAOmks9e0l0aNGrksgxRwD2Q2m5k5fSanfj7F3oS9NG/enD9S/4CW+TtUBe/63uzdu/eaYwcOHMjr\n/32dyosq4/+hP7d3vp2pk6e69g8gHDbt/Wm0PN+SCRkTGGkfyVjbWMY/Pv6a/TZu3EhE3QgsvhZu\n6ngTSUlJBqS9vpbtWxLvH0822Wg063zX0aJNC6NjFdu8RfN4p9I7jAsexwjLCHr/X2/69u3rsuvL\nKJRyoHLlypi8TJAM1AbSwX7KToMGDa7Z126307JFSz6f8zlRUVFUqVLlmn2E+7tw9gI1smtceV+L\nWvz+x+9/2efUqVMM7TeUcZfH0YpWfLPrG/pH92fP4T1uc1N33NPj2LxhM8M3DcfXy5fqDaoz671Z\nRscqti5dunDo+CH27NlDWFgYLVu2LPogJ5ICXgqZmZm89MpLbN2xlcgmkbz4wosEBQUZlsfb25v5\nc+dz17134VPLh+wz2Tx8/8N06tTpL/tlZWXR89ae7D26F69KXpjOm9i4diPNmzc3KLkorQFDBnDP\nnHtoa2tLFaow038m/Qf3/8s+P/74I5FekXShCwDD7cP5+uTXnDt3zm1uEvr5+RH7fSxHjx4lJyeH\niIiIAm++u7PQ0FB69eplyLVlFEoJaa3pM7APG3/ZSHrTdPyO+9E4tzE7t+3Ex8fH0GxJSUns27eP\n2rVr06pVq2s+f//993l6+tOk35EOXqASFO0vtGf7pu0GpBWOmv3JbGImxJCekc5tw27j3Q/exdfX\n98rn69evZ8SAEXxw+QN88CGFFO7zuY8LqRcwm80GJhfFITMxy0BycjKNmzcmY3RG3vcXDZU+qcR3\nX3zHjTfeaHS8Qo35zxje3fcudMvfcB7Cvgnj9InTfPzxx+xN3EvbVm0ZMWKEx7WCxLVyc3O5feDt\nHFl/hMisSOJ94nly0pM88eQTRkcTxSAzMcuA3W5Heam/LqBrytvu7rp06sKsRbNIa5cGZvDZ6UOH\nDh0Y+n9DWbt/LbYGNiwrLXy35jsWLVjkNv2kFZ3dbuelmJdY+PlCAgICmPTGpGKt9ufl5cXCZQtZ\ntGgRJ0+eZFSnUbLeSTkjLfAS0lrTrVc3dv6xk4zmGfj87EPdi3XZv3O/238t1VozdtxYZsyYgcnX\nROPGjZn+9nRuHXIrtpG2vH/Os8F/mj/7Eva5dDiUuL7nJzzP0neXMtI2kvOc523L2yxfu5zOnTsb\nHU2UIZnIUwaUUqyKXcX93e6n/c/tubPJnWzdsNXtizfkZX9n6jucO3OO44ePs2f7HsxmM94B3v/7\nLuYN3hZv0tLSDM0q/ufLuV8y2jaaZjSjG90YZBvEoq8WGR1LuAHpQimFwMBApr873egYpRYUFHRl\n1EyLFi0IMgWRtjENe1M7pgMmqgRUoWnTpganFH8ym82kknrlfaoplToBdQxMJNyFtMArOLPZTPy6\neLqbulMtthpWs5X4tfF/Gc0gjDXxlYm8bnmdr/iKGaYZxFeK56GHH3Jphg0bNtC4TmMsvhZ6dOpB\ncnIJHtFXRnJzc42OYDjpAxfCA6xZs4ZvvvyGgKAAHhv9mENPkC+p5ORk2jRrc2VC0CLTInbdsItd\nB3cZcqN72bJlPHzvw5xLPUdUmyi+XPYltWrVcnmOsibDCIUQDlu0aBHvPfAek1InAaDRDPYdzPHT\nx12+/vXhw4e5se2NTEqfRBOaMM80j8PND7N1z1aX5nCFMr+JqZR6Qyl1UCm1Ryn1jVIq2JHzCSHc\nT2hoKKdyT5FDDgAppJCjcwx5RN3mzZvp5NWJFrTABx/utd/Lzv07ycws/QOZPZmjNzG/B57RWucq\npV4DJgDXrqgjRAVht9uvjLuOioqie/fuRkdyWI8ePWjetTn/if8PTbOasslnE6/8v1fw8/NzeZbw\n8HBOqBPkkIM33vzCL/j7+VfYezZO60JRSg0Fhmmth/9tu3ShiAohNzeXvoP7sunAJrJqZOFz2IeX\nnnuJJ8Z4/sxHu93OwoULSUpKolOnToZNCLLb7QzpM4TjW4/TyN6IzWxm6gdT+dc9/zIkT1lyaR+4\nUioWWKC1nv+37VLA86Wnp2M2m2WGYzm1Zs0ahtw3hMv3X857Punv4DPDh7RLaYavk1Oe2O12li5d\nyq+//kqXLl1o166dU86bk5ODyWRym99Pp0ylV0qtBqoX8NGzWuvY/H0mAll/L95/iomJufLaarVi\ntVqLumy5cvToUfoN7sexn45hCbTw2ezPGDJkiNGxhJNduHABr1CvvOINEAzKS3H58mVCQkIMzVae\nmEwmbrvtNqed7/Tp09wx4A627d5GUEAQMz6ewT/+8Q+nnb+44uLiiIuLK9ExDrfAlVL3AQ8B0Vrr\njAI+r9AtcK01jZo14kSDE+jOGk6B/0J/9mzfQ+PGjY2OJ5woKSmJZq2akdYvDeqCaZuJiPMRJO5O\ndJtWnbhWt/bdaLinIffa7+UYx3jW/1nWbl1b4IqeruSKUSh9gKeAwQUVbwF//PEHySeT0VE6bwGs\n2uDd0JuEhASjowknq1OnDsuXLKfW5lr4vO1Du6x2fL/8eynebsxut7N191busd+DCRMRRHAjN7Jp\n0yajoxWLo6NQ3gN8gdX5/5Nu0VqPcjhVOVKpUiVMykR2SjaEA9mQeyaXGjVqFHmsKBm73c7Zs2ep\nWrUq3t7GrBLRo0cPko8bP0tRFI/JZCIkMISjqUdpSlPs2PnZ9DPVqlUzOlqxONQC11o31lrX01q3\nzf+psMU12DdOAAAPxElEQVQ7JSWF/kP6U61ONTp178SBAweAvP9BPvzgQywLLAR8G0DgnED6W/vL\nsp5OtnXrVmqH1SayQSRhlcNYvny50ZGEh5gxawYT/ScyxTKFMYFjqNuxLoMHDzY6VrFUyJmYWVlZ\njHtmHIuXLSY4OJh333iX6OjoUp9Pa03rDq055H+I7LbZqOOKkO0hHD109MrNq3379pGQkECtWrXo\n3bu3fK12ooyMDOpVr8foP0bTla4c4ADPBzzP/iP75ZuOKJb9+/ezefNmwsPDGThwoFs80EQe6HAd\nj419jHnr5pF+SzqnLp5i0LBBbF6/mdatW5fqfL/++itHjh0he2w2ZIFupck5lsO2bduuLLzfsmVL\nlz/wtKL45Zdf8LP70ZWuAEQSSQPvBiQmJkoBF8XSokULWrRoYXSMEquQBfzLr74kfXg6hADVITM5\nk9jY2FIXcIvFgj3TDt8AB/O2pfnL2F9XqVatGr9n/04yydSmNhe5yImsE9SpI0uuivKtQi4n6+/v\nD1c9r8A73ZuAgIBSn69y5cq0a98OfgeeBp4BXU3z3ervHM4qila5cmWmvDuFsf5jeSHoBUZaRjL6\nydE0adLE6GgVzvnz51m4cCFLlizBZrMZHafcq5B94LPnzOaxJx/D1s6G9x/ehJ4OZf+u/VStWrXU\n5+zVtxfrKq2D5vkbjkCnk53YtmGbc0KLIh06dIjExEQaNWpEmzZtjI5T4Rw9epQeUT1omNUQm7aR\nWS2T+IR4KleubHQ0jySPVLuO+++7n8XzFzOq6Sgm9JnA3h17HSreADc0uAGf5P91mXgne9OgXgNH\no4oSaNq0KcOGDZPibZBxo8Yx6OIgXrz0Im9efpP6SfV57eXXSnSOmR/NpEH1BtSsUpOnxj5FTk5O\nGaUtHypkC7wsnDt3jg5dOnDB6wJ4QaAtkIQtCdSsWdPoaEK4RPum7bnv8H20JO9m/QpWcPq208z7\nel6xjo+NjWXknSN5zvYcgQQyxTKFwU8MJualmLIL7cakBX4dFy9e5IF/P0DnmzozavQoLl265PA5\nq1atSuLuRD5/83PmvT6PQ/sOSfEWbstutzv9nN16dmOxeTFZZHGJS6y0rKRrdNdiH7904VKG2YbR\nlKbUpjYP2h5k6VdLnZ6zPKlwBTw7O5tuPbsxb/c8fqz7I59s/oToPtFOeb5eQEAAAwcOZNCgQVce\nGiyEO1m7di21w2rj6+NLmyZtOHLkiNPO/eqUVwmyBjHIexC3m26n5z09eeTfjxT7+OAqwZwxnbny\n/ld+JbiyPCOmMBWuCyUhIYFeQ3px6cFLeWuT5IJluoVdm3YRERFhdDwhrsjNzSUhIYG0tDTat2/v\ncKPg1KlTtGrSimfTnqUNbViqlrKy9koOnziMl5fz2nJpaWl4e3uX+IEPycnJdGrdifaX2hNoD+R7\n8/cs+W5JuXgoRmnIRJ4C5P+lgCavgOu8H5kZKdxJdnY2g24ZxMHtBwkxhXDe7zxrNq1xaAXLHTt2\n0MzUjPa0B+A2fRvzz83nt99+c2p3X2mH5NauXZudiTuZO3cuGRkZTLxtokdOrnGlClfAW7duzQ21\nbuDg8oNkNsrEfNhMu1btuOGGG1yWQWuN1tqprR5Rvnz00Udc3HaRWemzMGFikdciRt43kh82/VDq\nc4aHh3PSfpIMMjBj5jd+I92e7lbD/KpXr87TTz9tdAyPUeEqiLe3NxvWbOCRHo/QM7Unj/d93GVL\nfmqtGT9xPGaLGT9/P/75r3+SlZVV5tcVnuenAz/RNr0tpvynQ3TK7cSxY8ccOmfnzp3pNagXjwc8\nzlv+bzHWMpbXJ7+OxWJxRmRhgArXB26kWZ/MYkzMGNLuSANf8F/mz8j+I5kyeYrR0YSb+fTTT5k8\najKTbZOxYOEjn4/I6JXB4u8WO3RerTUrV67kl19+oV27dnTu3NlJiYWzufSZmIWEkAKeb9Dtg4hN\njoUbgPrAaWi2uxkHdh4wOJlwN1prRj4wkvmfz8fibaFG3Rp8F/edx6xTLRwnBdyNJCcn06RFE2yV\nbHk3T7OAlhDtE80PK0rfrynKt5SUFNLS0qhbt65bLHFqlO3bt3PkyBEiIyMrzExbmcjjRsY8OYbM\ntpnwIDACqAE+m3149813jY4m3Fh4eDgNGjQoVfHWWjPtvWm0adyGDs06MH9+gc8cd3uTnp/EIOsg\nPvn3J9zS9RbemfqO0ZHchrTAXaRN5zbsaboHGuZv2AvWNCvrvltnaC5Rfn304Ue89p/XGGsbSxZZ\nTLFM4cMFHzJo0KAij83JycFkMhk+vPbYsWN0bNGRWRmzCCGEFFJ40O9Bjp48Snh4uKHZylqZtsCV\nUi8qpfYopXYrpdYopWTx5UJ0v7E75t1myAGywLLfwq29bjU6lijH5n88n4dtD9OKVnSgA8Ntw1nw\nyYJCjzlz5gw9OvXA7GsmtFIon839zEVpC3b69Gnq+NUhhLwnW4UTTphvGGfOnCniyIrBkS6UyVrr\n1lrrNsAS4AUnZSqXJr8ymZvq3ITPFB98pvgwqNMgnvzPk0bHEuWY2d9MKqlX3qeqVPwD/As95q6h\nd1FzV01W6VW8mfYm40aOIyEhoayjXldkZCTJOcnsYAcAm9jEZa/LNGzYsIgjK4ZST+TRWl+9AlQg\ncM7xOOWXv78/q75dxcWLFzGZTLJWiihzE16cwG39buOc7RxZZBFriSVufFyhx8T/GM8S+xJMmGhI\nQ7rbuxMfH0+HDh1cE/pvQkNDWRS7iH8M/Qc2m43goGCWfrvUoQewlCcO3cRUSr2slDoJ3AuUbOHf\nCiokJESKt3CJHj16sHLdSvwe8aPyo5XZ+OPGIp/LGlY5jKMcBcCOnZ+9f6Z69epOzXX+/Hkef+Rx\n+vfsz6T/TipyMlvPnj05c/EMSWeSSD6bTFRUlFPzeLJCb2IqpVYDBf3Xe1ZrHXvVfuOBJlrr+ws4\nh9NuYqamppKZmUnVqlUNv7kiRHkUGxvLff+4jy6qC794/UL1NtVZuW4l3t7OWXUjPT2dDi06EJEU\nQdvstqzyX0XNm2uycNlCp5y/PHF4MSutde9iXms+sOJ6H8bExFx5bbVasVqtxTxtntzcXP792L+Z\nM3sOXiYv2rZry8plK91qDQchyoOBAwcSvyOe+Ph4qlatysCBA51WvAHi4+MxnTUxOns0CkVUehTD\nVg3j/PnzhIaGOu06niguLo64uLgSHVPqYYRKqcZa6yP5rx8HOmmt/1XAfg63wD/++GPGvDwG2502\n8AXflb7c1vQ2Fswt/I66EMK9fP/99zxzxzO8lfoWAFlkcbvv7fx86meHH2tY3pT1crKvKqWaAHbg\nGDDSgXMVauOWjdgibWDOe5/VNout67eW1eWEEGWka9eupFVO4wPbB7TJacNK80qsN1mleJdSqW9i\naq1v11q31Fq30VoP01qnODPY1SIaRWBONkP+Q3O8fvGifr36ZXU5IUQZCQgIIH57PAF3BrC682pu\nfPRGvlr6ldGxPJZHzMS02Wx0tXblaMpRvCxe+Fz0YcuGLQ4tbi+EEO6sXC1mlZ2dzYYNG0hPT6dr\n166EhIQ4IZ0QQrinclXAhRCiIpHVCIUQohyTAi6EEB5KCrgQQngoKeBCCOGhpIALIYSHkgIuhBAe\nSgq4EEJ4KCngQghDaa05efIkP//8M7m5uUbH8ShSwIUQhsnMzGRI3yG0a9KOLi26YI2ykpqaWvSB\nApACLoQw0OTXJnNuwzkWZCxgfvp8Ku2txPj/jDc6lseQAi6EMMyebXvomd4TH3wwYSI6M5o9CXuM\njuUxpIALIQzTuHljtvttJ5dcNJptPtto3ExWGS0uWcxKCGGYy5cv07t7b84ePYuP8sErzIt1W9YR\nHh5udDTDyWqEQgi3l5OTQ0JCAna7nfbt22M2m42O5BakgAshhIeS5WSFEKIcc7iAK6XGKaVylVJV\nnBFICCFE8ThUwJVSdYDewC/OiWOcuLg4oyMUi+R0Lk/I6QkZQXIawdEW+FTgaWcEMZqn/EeVnM7l\nCTk9ISNITiOUuoArpQYDyVrrvU7MI4QQopi8C/tQKbUaqF7ARxOBCcAtV+/uxFxCCCGKUKphhEqp\nFsAawJa/qTZwCuiktU75274yhlAIIUrBJePAlVLHgfZa6wsOn0wIIUSxOGscuLSyhRDCxcp8JqYQ\nQoiy4dKZmO4+6Ucp9aJSao9SardSak3+OHe3o5R6Qyl1MD/rN0qpYKMz/Z1S6g6lVKJSyq6Uamd0\nnr9TSvVRSh1SSh1RSj1jdJ6CKKU+UUqdUUrtMzpLYZRSdZRS6/L/e+9XSo02OlNBlFJmpdS2/N/v\nA0qpV43OdD1KKZNSapdSKraw/VxWwD1k0s9krXVrrXUbYAnwgtGBruN7oLnWujXwE3kjgtzNPmAo\nsMHoIH+nlDIB7wN9gEjgn0qpZsamKtBs8jK6u2zgCa11cyAKeNQd/z611hlAz/zf71ZAT6VUN4Nj\nXc8Y4ABFdE+7sgXu9pN+tNaXrnobCJwzKkthtNartdZ/PjxwG3mjgNyK1vqQ1vono3NcRyfgqNb6\nhNY6G/gCGGxwpmtorTcCF43OURSt9W9a6935ry8DB4GaxqYqmNb6z5FzvoAJcLuBF0qp2kA/4GOK\nGJ7tkgLuSZN+lFIvK6VOAvcCrxmdpxhGACuMDuFhagFJV71Pzt8mHKSUqg+0Ja9h4XaUUl5Kqd3A\nGWCd1vqA0ZkK8BbwFFDkE54LnchTEp4y6aeQnM9qrWO11hOBiUqp8eT9Rd7v0oD5isqZv89EIEtr\nPd+l4fIVJ6Obkjv3ZUApFQgsAsbkt8TdTv431zb5941WKaWsWus4g2NdoZQaAKRorXcppaxF7e+0\nAq617n2dQC2ABsAepRTkfd3foZS6ZtKPK1wvZwHmY2DLtqicSqn7yPuaFe2SQAUowd+luzkFXH2D\nug55rXBRSkopH+BrYJ7WeonReYqitf5DKbUc6ADEGRznajcCg5RS/QAzEKSUmqu1vqegncu8C0Vr\nvV9rXU1r3UBr3YC8X5R2RhTvoiilrn4Y32Bgl1FZCqOU6kPeV6zB+Tdm3J27LbOQADRWStVXSvkC\n/wCWGZzJY6m8ltks4IDW+m2j81yPUqqqUqpy/mt/8gZVuNXvuNb6Wa11nfxaeSew9nrFG4x5oIM7\nf319VSm1L7+PzAqMMzjP9bxH3k3W1flDjaYbHejvlFJDlVJJ5I1KWK6UWml0pj9prXOAx4BV5N3p\n/1JrfdDYVNdSSi0ANgMRSqkkpZQh3XnF0BUYTt6ojl35P+44eqYGsDb/93sbEKu1XmNwpqIUWi9l\nIo8QQngoeaSaEEJ4KCngQgjhoaSACyGEh5ICLoQQHkoKuBBCeCgp4EII4aGkgAshhIeSAi6EEB7q\n/wNvF9zxlJjA5gAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10dc0ba10>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot number of clusters through the iterations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(range(max_iter), number_of_clusters)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "[<matplotlib.lines.Line2D at 0x111a8af90>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+wJXV55p/n3jv3zgwRiEXFZR1Ssko2P0oEljWsxeph\nFYMTQcs/srqRuCRuJsm64lprTKioM0k21KZMzCZu2PFHLJaqDUkxVnYoKY1LOOOSxQHiDAxOhsCS\nlBBlDEEQnCHeH+/+0aen+/b0r3NOd3/fc/r5VN2659zTt/s979v3PU+/T3dfmhmEEELMHwuhAxBC\nCNEOavBCCDGnqMELIcScogYvhBBzihq8EELMKWrwQggxp9Rq8CQXSR4ieXvOaz9J8gGSD5L8c5IX\nNh+mEEKIcVmqudz1AI4CeFHOa48BeK2ZPUvyKgCfAHBZQ/EJIYSYkEoFT3IHgJ0APgWA2dfN7B4z\ne3b09CCAHY1GKIQQYiLqjGg+BuADADZqLPszAO6YKiIhhBCNUNrgSb4ZwDfN7BBy1Htm2SsA/DSA\nDzYXnhBCiElh2b1oSP4GgGsBrAHYCuBMAPvM7Kcyy10I4LMArjKzRwvWpZveCCHEBJhZqcAuolTB\nm9kNZnaemZ0P4O0A/iynuX8/oub+zqLmnlqfvszwkY98JHgMXr6UC+VCuSj/moa6Z9Gc6tEAQHLX\nqGHvBfBhAN8L4CaSALBqZq+eKiohhBBTU7vBm9kBAAdGj/emfv5uAO9uPjQhhBDToCtZAzAYDEKH\n4AblIkG5SFAumqHUZG10Q6R1tS0hhJgXSMLaMFmFEELMLmrwQggxp6jBCyHEnKIGL4QQc4oavBBC\nzClq8EIIMaeowQshxJyiBi+EEHOKGrwQQswpavBCCDGnqMELIcScogYvhBBzihq8EELMKWrwQggx\np6jBCyHEnKIGL4QQc4oavBBCzClq8EIIMaeowQshxJyiBi+EEHOKGrwQQswpavBCCDGnqMELIcSc\nogYvhBBzihq8EELMKWrwQggxp9Rq8CQXSR4ieXvB679L8hGSD5C8uNkQhRBCTEJdBX89gKMALPsC\nyZ0AXmFmFwD4WQA3NReeEEKISals8CR3ANgJ4FMAmLPINQBuBgAzOwjgbJIvaTJIIYQQ41NHwX8M\nwAcAbBS8/lIAj6eePwFgx5RxCSGEmJKlshdJvhnAN83sEMlB2aKZ56eNcor49V8HDh5Mnl93HfC2\nt52+3Mc/Drz1rcCOzEfH2hpw7bXA88/X3WI+O3cCP//z061jWj70IeDw4eT5z/0c8OM/3m0MDz4I\n/MqvAFajgm96E/ALv9BOHL/2a8D73w+ccUY7658HvvGNKP9ra9XLXnwx8Ku/2n5MX/kK8OijwE/8\nRPvb8s6Xvww8+WTUt0JR2uABvAbANaM5+1YAZ5L8H2b2U6ll/hbAeannO0Y/O43du3efejwYDDAY\nDLBvH/CudwEvfzlw++3AgQP5Df6P/gj4wR88vcE/9xywfz9w660V76SE++4D7rgjfIP/4z8G3ve+\n6D3edhtw993dN/jDh6Ocvv/95cvFOWurwX/iE8A73gG84hXtrH8eeOwx4Ngx4Dd/s3y5xx8Hfu/3\numnw990H3HOPGjwQ5eGRR8Zv8MPhEMPhsJEYShu8md0A4AYAIPk6AP8p09wBYD+A9wC4leRlAJ4x\ns+N560s3+JjVVeD1rwde+cpoR3zoofxYVlejr7yfb98OXH112TspZ8sW4N57J//9plhdBd74xujD\n7uGHo0//EDGcf351PpeXNx95tRFHXr1Fwuoq8H3fV12rRx4Bfvu3u4tJdYuYNBex+I3Zs2fPxDFU\nKfgsBgAkdwGAme01sztI7iT5KIDvALhunBWurQFLoyiWlooPN9fW8l9L//6klG23S+rmoqsYymg7\nvqJ6iwQvtUqjuiV4yEXt1mhmBwAcGD3em3ntPZMGsLoaKWgg+l70iVem4OPfn5Sy7XZJ3Vx0FUMZ\nbccnJViNl1qlUd0SPOQi+JWsa2tq8DFq8Jvj8FATz3ipVRrVLWF1NbyCD97gV1c1oonRiGZzHB5q\n4hkvtUqjuiWsrYX/sHPR4KXgI6TgI8ykBOvgoVZZVLcED7kI3uBlsiZIwUdsbCSxiGI81CqLFHyC\nh1wEb/BS8AlS8EkM6e8iHw+1yuJBtXrBQy7U4Cu22xXr69H3xcXouxp8+Jp4p26tlpai/avO1clN\nxKS6RXjIRfAGrxFNRPZ99HlEE683dE28U7dWZCQcusinh7GEFzzkImiD39iIvuqo1nlX8Nn3IQUf\nvibeGWff72p/8qBaveAhF0EbfKxAOLpVmRR88lwKPnxNvDPOvt/V/uRBtXrBQy6CNvi6qrXstDkp\n+HbjKEIKPjxS8L7xkIuZaPCxAakG330cRajBh0cN3jcecuFiRBNTdBhZdsiuEU27cRShEU14NKLx\njYdczISCL1N0UvDtxlGEFHx4pOB94yEX7hR8XkLaVvALC8kZPaGYNQW/MNpz2siZFHw9pOB94yEX\n7hR8XkLaVvBk8ba7YtYUPNBejFLw9fBQqyweVKsXPOTCXYMPMaIp23ZXqMFvjiH9XeTjoVZZPDQ1\nL3jIhbsRTQiTtWzbXTFrIxqgvRg1oqmHh1pl8TCW8IKHXEjBV2y7K6TgN8eQ/i7y8VCrLKur3d33\nxjtS8Gv1mlqVgm+iwXtQ8On3EVLB182nFHxYxqlVVx6TapcgBb9abyxRpeCbGNF4UPDp9xFSwdfN\npxR8WMapVdEZak2j2iWsrkb9LOTRTPAGn6daswnRiCZcHGWowYfFQ62yqHYJcQ5CqvjgI5q0Aolv\naxrfmiC9XPp72TomxcOIRiZrEkP6u8jHQ62yqHYJHnLhSsED+Uqjjwo+PqTu+vDOgyqUCqyHh1pl\nUe0SPOQiuIKv0+C7Mlk9NfiFhegrezTTdRxltJUzD38Ys4AavG96P6LJM4nyDiVXV4GVlfZNVk8j\nGiBMTOMc9rcV39paVG8d5pfjdUSj2kXEueitgh9nRLN9e79GNECYmDyowrJ6iwQPtcqi2iV4yEXw\nEU0dBb+2Bmzb1i+TFQgTkwdVWFZvkeChVllUuwQPuZgZBb9tmxR8qDiKaFPBF9VbJHioVRbVLsIs\nauxbt/ZcwXtp8N5MViBMTF5MVjWJatTg/bK2Fp3yvbzsWMGT3EryIMnDJI+SvDFnmXNIfn60zEMk\n/23djdc1WbsY0chkLY6jiDZN1tCHtrOAtxFNWrX2vXZxbUILx9IGb2YvALjCzC4CcCGAK0henlns\nPQAOjZYZAPgtkrV2O5msCRrRbI4htDk1C3ioVZr19ejU3tBnjnggrk3ovlI5ojGzE6OHywAWATyd\nWeQbAM4cPT4TwN+bWa3P7yJjMe88eJms7RMrsNCqUAq+Hh5qlY0nbmp9r52XXFQ2eJILJA8DOA7g\nLjM7mlnkkwB+hOTXATwA4Pq6Gy9SrXnnwctkbZ9YgS3UdGZksobFm4KPR66hxxIe8JKLys9/M9sA\ncBHJswB8geTAzIapRW4AcNjMBiRfDuCLJF9lZs9l17V79+5TjweDAdbWBjJZR8QXc6XpOqZxcymT\nNSweG7yHsYQH0rkYV8EPh0MMh8NG4qg93DCzZ0l+DsClANJbfw2A/zxa5v+R/GsA/xTA/dl1pBs8\nAHzpSzJZY9bWgDPO2PyzrmMaN5cyWcMy7ojmH/6hm3hCjzs9EI9oJhFBg8EAg8Hg1PM9e/ZMHEfV\nWTTnkDx79HgbgCsBHMosdgzAG0bLvARRc3+szsZlsiZ4GNGMm0uZrGGRgvdLPKIJnYuqz/9zAdxM\ncgHRh8EtZnYnyV0AYGZ7AfwGgM+QfGC0zC+aWdaIzSVWapsCksl6iq5jGjeXbZqsZ54pFViFTFa/\neMlF6e5hZkcAXJLz872px08BuHqSjY9jsm7dmvyvR7J8HZMQ+pNWCn5zHJrBV+NRwXswFj3gJRcz\ncyVr0TxrnkzW0FeyejJZNaKpxmOD14gmYhqTtUmC34umrslalKx5MllDX8kqk3V22NiIvhYX6y2v\nEU23pHPRWwU/jsladLijEU27MZShEU044g/j9LiyDI1ousVLLmZiRFOl4Jsa0YRW8Hkjmq4V/Lgj\nGin4MIxbqy5UtRR8gpdcBFfwdf+jU5mCb2pEE1rB541oulbw445opODDMG6tulCSXlSrB7zkYiYU\nfJl5I5O13RjKkMkaDi/jtDQyWRNksqJYwY87opHJ2l4MZbRpsq6sJEaiOB0v1yyk8TKW8IBMVox3\nHrxM1jAxlNHmiCY+LbbvjaIIL7VK42Us4QEvuZiJEU1XJqu3Bt/nEY0H9eMZrw1edYvQiAb+TFaN\naPyMaHTTqnK8jmhUt4hpbjbWJDOh4IuUgVl0+4J5OYtGI5rNcYSuiWe81CqN6pYQC8/QwtGlgq87\nollbi67kq3uxRxmhVYduNrY5Dg+Ht57xUqs0qluCTFZMb7I2ZbDG25WC96EKvRhUnvFSqzSqW4KX\nXMzEiKZMwTfV4EMXQibr6XGE/tD1jNcGr7pFyGTF9CZrUwYrEL4QMllPjyP02MwzXkc0qluETFZM\nb7JqRNN+DGXIZA2Hl1qlUd0SZLKiGZO1KQUfWnXIZN0ch4fDW894qVUa1S1BJitksqaRgt8chweD\nyjNeapVGdUvwkguNaEaELoRM1tPjCP2h6xmvDV51i5DJivFGNHnmTZMjmtCFkMl6ehyhx2ae8Tqi\nUd0iZLKiWMHnjWhksoaJoYw24ktfnRy6Jp7xUKssUvAJMlkhkzWNTNbNMZDh/zg846FWWWSyJshk\nhUzWGLP6fkSbeFCF6Q/90Ie3nvFQqyxejEUPxPUJnQuXI5q+maxF99Tpo8majiG0+vHMuLVaXIy+\nr6+3Ew+gEU2a3o9ozHyNaEIWouh99NFkTX/oh/7j8Mwk+37bYxqNaBJ6P6JZXwcWFqKvNH0c0RS9\nD41opASLmGTfb3t/0ogmwUsugjX4ohuF5e2EXd1sLKSCz3sfIUzWcUc0UvBhmGTfbzufUvAJXnIR\nrMEX3Sgs2zTi0+YWF9u/2VhIBV80oulawY87opGCD8Mk+37b+fSiWj0wEyYrya0kD5I8TPIoyRsL\nlhuQPETyIZLDOhuuq+DjnSY+bW4eTdai9yGTVY2iCK8jGg9zZw94MVlLNYCZvUDyCjM7QXIJwN0k\nLzezu+NlSJ4N4L8B+DEze4LkOXU2XKbg0ztH2SG7TNZu4ihCJms4ZLL6ZmZMVjM7MXq4DGARwNOZ\nRf4NgH1m9sRo+afqbLjMWEzvHGWH7G0oeLNm1jcOs2qyxg2jyZxpRFMPrwpeI5oIL7mobPAkF0ge\nBnAcwF1mdjSzyAUAXkzyLpL3k7y2zobrjmiqFHxTDT4+o6fN84SLmNURDRl5I02qNY1o6uG1wXtQ\nrR5I58LtiAYAzGwDwEUkzwLwBZIDMxumFtkC4BIArwewHcA9JL9sZo9k17V79+5Tjy+4YIClpcHp\nAS2Np+CbGtEASTGaXGcdZnVEAyQxNvVBm44h9O0jPON1RKObjUVMc7Ox4XCI4XDYSBy1dxEze5bk\n5wBcCiC99ccBPGVmJwGcJPklAK8CUNrgjx2rb7IWKbomRzRAUoytW5tbZx1mVcEDSYzbtjUfg5Rg\nMVLwvpnGZB0MBhgMBqee79mzZ+I4qs6iOWdkooLkNgBXAjiUWex/Abic5CLJ7QB+FEB2jHMa3kzW\nvPV3xTwo+CZjkMlajVcF72Es4QEvJmvVLnIugJtJLiD6MLjFzO4kuQsAzGyvmR0j+XkADwLYAPDJ\nnDn9aXgzWeNthyjGrJqsQPMxZut98mRz654nPNQqS1y7+Or0+PqVPuLFZK06TfIIovl69ud7M88/\nCuCj42zYm8kKhCvGPIxo2oghtPrxjNcGn61dnxu8h6MZl1eyFhmpXZmsXaMRTX4MMuuK8TqiUe0i\nYvG5sABsbERfIXB3L5r4Ez9OSAiTtWuk4PNjkIIvZlYUfF/JXoEf6sPOnYIHNjcNmazh4yhDJmsY\nvCp41S4im4tQH3ZBG3yRAknvHH02WRcXI6Oqq6trPahCXclaDw+1yqLaJaTrEzIX7kY0wOYdsc8m\nK9ltTBrRzA5eG7xqF5H+sNOIJkO6afTZZAW6jcnLiEZGXTVeRzSqXURafErBZ8iOaPpqsgJS8H1W\ngWVIwftGCl4m6ymk4DfHIKOuGq8KXrWLkMkqk/UUVbnwrOBlsobBQ62yqHYRGxvRiRHxKd8a0WRI\n74ga0fhu8BrRhMFrg1ftTp9OaESTITuiKTJuNKLpLo4iZLKGwduIJr5aM61a+1q7rHiVgs8gkzVB\nCr65dc8T3hR8/IFDtr8t70jBo/7cucvz4EMq+Dofdm2SVWB1kckahkn2/TbzmY2nz7XLy0XvFLzX\n8+BDKfiyEU0XMWUVWF1ksoZhkn2/zXxm4+lz7bLiVSOaDBrRJHQV06S51IgmDN5GNNl4+lw7jWig\n8+DTeDBZJ82lRjRh8GayakSToBENxlPwXZ0HLwU//u+1oeA1oqnGo4LXiCYib0TTSwXv0WT11uC7\nimnSBt/GDF4jmmo8NngvqjU0eSOa3in4SUzW9P96rFrHJGhEM/7v6Tz4MHgc0WQVfF9rp/PgMZnJ\nGr+Wd5VrE2hEM/7vyWTtHjMpeM/IZMVkJiuwOVkyWbuJoQyZrN2zsRGdzrow5l+vTNZukMmKyUxW\n4PTxjRR8+zGUIZO1e7z4JWlksibIZMVkJmv8WlrBy2RtP4YyZLJ2j5dapdGIJkEmKyYzWatemxaN\naMb/PZms3TNprWSydoNMVshkTaMRTX4cca27+qfjs4IUvG9kskImaxop+M1xxPVeWIi+NjaaW/88\n4FXBy2SNkMmK8f7hh0zWsDGU0abJ2sb65wGvCl51i5DJimpjschIlcnafQxltGmytrH+ecBLrdKo\nbgkzY7KS3EryIMnDJI+SvLFk2X9Oco3k2+psWCZrgkY0xXH02awrwuuIRnWL8GSylu4mZvYCySvM\n7ATJJQB3k7zczO5OL0dyEcB/AfB5ALXuKF41ovnud6PHRcogvl3BuP+gogyNaMb/vTZNVqDfSrAI\nKXjf5OXi5MkwsVSOaMzsxOjhMoBFAE/nLPYfANwG4O/qbnhak7VpgzW97q6Rgt8ch8y6crwqeNUt\nIlsftyMaACC5QPIwgOMA7jKzo5nXXwrgLQBuGv2o1klt05qsTRus2e12iWbwm+OQWVeOl1qlUd0S\nPJmslTrAzDYAXETyLABfIDkws2Fqkd8B8EtmZiSJkhHN7t27Tz1++ukBtmwZ5C5Xx2Rt2mAF/I5o\nuji804hmdvDa4FW3iGlN1uFwiOFw2EgstQ/0zOxZkp8DcCmA9Nb/GYBbo96OcwC8ieSqme3PriPd\n4Pftm85kbdpgBTSi8TKikVlXjtcRjeoWMa3JOhgMMBgMTj3fs2fPxLGU7iYkzwGwZmbPkNwG4EoA\nm7ZmZv8ktfxnANye19yzTHslaxsjGq8K3vOIRgq+e6TgfZOXC68jmnMB3ExyAdG8/hYzu5PkLgAw\ns72Tbrju3LnLEU1IBV/nwy5UDGXIZO0eL7VKo7oleDJZq06TPALgkpyf5zZ2M7uu7obLDjNDjmhC\nKfiyEU1XCn7SEY1M1m6ZtFZt5lJ1S1hdBZaXk+e6kjWDRjQJGtH0t1EUoRGNb2bmStY2Gec8+Dzz\nRufBdxdDGTJZu0cmq288Xcnq9mZjUvARfVLweVcn91kJFiEF7xtPJqvbEU0ok9Vbg+/ThU55MfS5\nURThoVZZ1OATNKKBX5O160/ajY3oH1oU/QPlPo1o8mLo86F+EdOMaNbX2/kHKhrRJOSNaKTgU/Rp\nRBO/DxZc/9unEY0UfD0mrRUZjb/aaDZS8Am9V/BmkZIouhNkn242VvU++qbg8xp8X5VgEdPs+22p\nSZ0Hn9D7/+gU76BFqjWdkLIRzTwp+CL6puDzRjR9VYJFTLPvt9VsdB58QrY+vRvRVO2g6U//eTdZ\n6+TCc4OXydo9Xhu8F9Uamt6PaKoM0vjTPz5tLm1AzpvJqhFNeQx9NuuK8Dqikcka0fvz4KvUd9w0\nyhSdRjTdxlGETNbukYL3Te/Pg6+r4MtMN5ms3cZRhEzW7vGq4L00tdB4utmYWwVfNIZp02RdXGzv\nPOEiZl3BLyxE5/JvbDQTg0zWarwqeJmsETJZa5qsVQq+6QZPdr9jzrrJSjYXo0Y09fDa4DWiiZDJ\nWnNEU6Xgmx7RAN0fWs76iAZoLkaZrPXwOqKRyRohk9WpyQr4U/DeRzRAczFKwddDCt43Mlmdmqzp\n9XeFFPzmGGSyVuNVwWdVa9d+lhdksjo1WdPr7wop+M0xyGStxquCT9euzfveeEcmq1OTNV6/pwbv\n3WQFZLJ2jYdaZVHtEorOKApxNON2RBPP4GWy9mtEI5O1Gq8jGtUuIis+Fxair/jK/C5xOaKJD+9O\nnpTJ2rcRjVRgNVLwvinKRYgPO5cjGiB6vajBtz2i6VrB1xlXhY6jDJms3eKhVllUu4S8o5lQH3ZB\nbxdcxtJS1OC7HtGEUPBl7yO+0Vrbh3fT5FMma7d4qFUW1S4hT8CGGldJweds19OIBugmJg+H/TrM\nr4eHWmVR7RI85cKlyQqEU/DeTFagm5hkss4OMll9U5SL3jT4Oup7yxbgxAmZrEA3MclknR2k4H0j\nk1UjmlNoRFMeQ1+bRBkeapVFtYswk8nq3mT1NqLpIiYPh/06zK+Hh1plUe0i1tejU7yz/2/arclK\ncivJgyQPkzxK8sacZX6S5AMkHyT55yQvLFtnXdUaYkQjBT8+UvDd4qFWacxUu5ii2oTKRaUOMLMX\nSF5hZidILgG4m+TlZnZ3arHHALzWzJ4leRWATwC4rGid05qsfbrZGCCTtW8qsApvCn5jI1KsCxm5\n2MfaFdUmlMlaazcxsxOjh8sAFgE8nXn9ntTTgwB2lK2vrsmqK1kj2o6pSIHVpUmTdWVl88/6qAKr\n8KbgvanWkJTlwuWIBgBILpA8DOA4gLvM7GjJ4j8D4I6y9clkTfAwollfT+6XMQka0XSLGrxfiqYT\nrk1WM9sws4sQKfPXkhzkLUfyCgA/DeCDZeuTyZrgwWSddtwlk7VbvI1oysYSfatdkfAMlYuxdpPR\njP1zAC4FMEy/NjJWPwngKjP7Vt7v7969GwBw4ACwY8cAwKBwW30yWbdvL1+m7ZimzaUUfLdMq+C/\n851u4ulj7ZrIxXA4xHA4bCSeygZP8hwAa2b2DMltAK4EsCezzPcD+CyAd5rZo0Xrihv8hz8cnUpU\nGthS9YhGJmt3MZTRpMnq5QIRz3hU8J7mziEpqs04DX4wGGAwGJx6vmfPnuKFK6izm5wL4GaSC4hG\nOreY2Z0kdwGAme0F8GEA3wvgJkYngK6a2auLVri2BmzdWr7RLVuA557Tf3QC2o9p2lzqZmPd4nEG\n7+nMkZAU1cbtiMbMjgC4JOfne1OP3w3g3XU3Oo3JurQUmYLzNKIJbbJqRDNbeGzwGtFEeMvFzF3J\nGv8zkBdekMnaZQxlyGTtFo8jGpmsEd5y4fpeNHkma9Vr0yIFPz5S8N1hFh3BTnNRmhR8e3jLhesG\nnzeiqXptWkKYrHVy0baCn7bBy2TthlghZu91Upc28imTNaEJk7VJZm5EU/XatIQwWevkom0FP+2I\nRiZrN3ipVRqZrAneTFYp+Jx1a0QzHhrRdIeXWqXxNpYIibdcuFbwVYaFTNbuYihDJmt3eKlVGm/G\nYki85cK1gk9/r/vatEjBj0+bCj4+LdZs+vXPA15qlcabag2Jt1yoweesWw1+PNps8GQ/Z7lFeKlV\nGm9NLSQzebOxpqk7lkh/r/vatGhEMz5tjmiaXP884KVWaVS3BG83G5OCz1m3FPx4tKngm1z/POCl\nVmlUtwRvuZCCz1m3FPx4SMF3h5dapVHdErzlYmYVfN4/tm0CKfjxkYLvDi+1SqO6JXjLxcw2+DbG\nM/G61eDHQw2+O7zUKo3qliCTFc2MaNoYz8Tr1ohmPDSi6Q4vtUqjuiXIZIUUfBop+Oo4+qgEi/BS\nqzSqW4K3XEjBZ+j6Bkn6j07VcfRRCRYhBe8bb7mQgs8Q4mZjVe+li5uNTZPPJm825kn9eEQK3jfe\ncqEGn7Nubw1eI5r+NYoivNQqjeqW4C0XGtHkrNvbiEYma/8O9YvwUqs0qluCt1xIweesWwp+PKTg\nu8NLrdKobgneciEFn0Em6/g0EV/Zv6LroxIsQgreN95yIQWfQSbr+DQR39pa8dXJfVSCRUjB+8Zb\nLtTgc9bdVSHM6it4zw2+ifjKYuhjoyjCQ62yeGtqIfGWC41octbd1aHU+jqwsBB9hYzJw2F/WQx9\nPNQvYtpaLSxEwmJjo/2Y+lg3b7mQgs9Zd1eftHXVmBR8/5RgEdPWimw+n95Ua0i85aLzBr+xEe1k\nVao1ZIPv6pO26L4VWbowWadt8E0o+LIG3zclWMS0tQKaz2dRTH2sW5GCD5WLzht80d3WsoQc0XSp\n4Ovmom0FP+2IpgkFXzai6ZsSLGLaWgHN57Mopj7WrUjBh8pFkAZfV7UCxZ+GGtF0H0cRGtF0x7S1\nAjSiaRNvuSht8CS3kjxI8jDJoyRvLFjud0k+QvIBkheXrbOuSbS0VHza3LyYrOPkQibrdOufF6at\nFdB8Pr0ZiyHxlovSBm9mLwC4wswuAnAhgCtIXp5ehuROAK8wswsA/CyAm8rWOY5qLVN0s6zgh8Mh\nACl4oF4u+qIE41yU0RcFXycXHpkpBQ8AZnZi9HAZwCKApzOLXAPg5tGyBwGcTfIlReublQZv1s76\ngflr8LE6mSRnavAJavAJavDNUNngSS6QPAzgOIC7zOxoZpGXAng89fwJADuK1jfOWKLskL2tEc3C\nQjQWavI84SLmZUQTn8u/vt5ODH081C9CIxrfeMtF5a5iZhsALiJ5FoAvkByY2TCzWHZSnqvlrr4a\nOHECWF6uDmzrVmBlZfzXmmBlBbjmmupTOSfl4YeBv/gL4NvfrpeLlRXga1+L8tcG994LXHfddOuI\nc7a4ON7au5nsAAAEYUlEQVTvxbn41reKc7GyAnz608CBA9PF6J04F2V89avAe9873XZWVoBdu4AX\nvWi69cTcf39+7ZaXgeefn2y/rZMLjxw5kp+LlRXgqafq5WJ5Gdi3r5l4aGMcV5P8EICTZvbR1M/+\nO4Chmd06en4MwOvM7Hjmd1scegghxPxiZjmnm1RTquBJngNgzcyeIbkNwJUA9mQW2w/gPQBuJXkZ\ngGeyzX2aAIUQQkxG1YjmXAA3k1xANK+/xczuJLkLAMxsr5ndQXInyUcBfAfAlAf7QgghmmCsEY0Q\nQojZoZMrWUleRfLY6GKoD3axTQ+QPI/kXSS/SvIhku8d/fzFJL9I8q9I/inJs0PH2hUkF0keInn7\n6Hkvc0HybJK3kfzL0UWEP9rjXPzy6G/kCMn/SXKlL7kg+Qckj5M8kvpZ4Xsf5eqRUT99Y9X6W2/w\nJBcBfBzAVQB+GMA7SP5Q29t1wiqA/2hmPwLgMgD/fvTefwnAF83sBwDcOXreF64HcBTJmVZ9zcV/\nBXCHmf0QoosIj6GHuSD5MgD/DsAlZvZKRNfavB39ycVnEPXGNLnvneQPA/jXiProVQB+fzQ+L6QL\nBf9qAI+a2d+Y2SqAWwG8pYPtBsfMnjSzw6PHzwP4S0TXDZy6OGz0/a1hIuwWkjsA7ATwKSSn1vYu\nF6NTjv+lmf0BAJjZmpk9ix7mAsC3EQmh7SSXAGwH8HX0JBdm9n8AfCvz46L3/hYAf2hmq2b2NwAe\nRdRfC+miweddCPXSDrbripFSuRjAQQAvSZ1pdBxA4ZW/c8bHAHwAQPoysj7m4nwAf0fyMyS/QvKT\nJM9AD3NhZk8D+C0AX0PU2J8xsy+ih7lIUfTe/zGi/hlT2Uu7aPC9d3FJfg+AfQCuN7Pn0q9Z5HLP\nfY5IvhnAN83sEE6/MA5Af3KB6Oy1SwD8vpldgujss00jiL7kguTLAbwPwMsQNbDvIfnO9DJ9yUUe\nNd57aV66aPB/C+C81PPzsPlTaK4huQVRc7/FzP5k9OPjJP/R6PVzAXwzVHwd8hoA15D8awB/COBf\nkbwF/czFEwCeMLP7Rs9vQ9Twn+xhLi4F8H/N7O/NbA3AZwH8C/QzFzFFfxPZXrpj9LNCumjw9wO4\ngOTLSC4jMgn2d7Dd4JAkgE8DOGpmv5N6aT+Ad40evwvAn2R/d94wsxvM7DwzOx+RifZnZnYt+pmL\nJwE8TvIHRj96A4CvArgdPcsFInP5MpLbRn8vb0BkwvcxFzFFfxP7Abyd5DLJ8wFcAODe0jWZWetf\nAN4E4GFEpsAvd7FND18ALkc0bz4M4NDo6yoALwbwvwH8FYA/BXB26Fg7zsvrAOwfPe5lLgC8CsB9\nAB5ApFrP6nEufhHRB9wRRKbilr7kAtHR7NcBfBeRV3ld2XsHcMOojx4D8GNV69eFTkIIMad0/i/7\nhBBCdIMavBBCzClq8EIIMaeowQshxJyiBi+EEHOKGrwQQswpavBCCDGnqMELIcSc8v8BacWgtjQE\nQBUAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d735cd0>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.mean(number_of_clusters)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "3.2200000000000002"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.unique(table)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "array([ 1.,  2.,  3.])"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for t in np.unique(table):\n",
      "    print 'Table', int(t)\n",
      "    print 'number of data points assigned to this table', np.sum(table==t)\n",
      "    print np.mean(x[table==t,:], axis=0)\n",
      "    print"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Table 1\n",
        "number of data points assigned to this table 20\n",
        "[ 1.99834569  1.75685838]\n",
        "\n",
        "Table 2\n",
        "number of data points assigned to this table 20\n",
        "[-1.99439971 -1.80973143]\n",
        "\n",
        "Table 3\n",
        "number of data points assigned to this table 20\n",
        "[ 1.89358224 -2.15001022]\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}