2015-10-14-collapsed-gibbs-sampling-for-mixture-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Note__: This is best viewed on [NBViewer](http://nbviewer.ipython.org/github/tdhopper/notes-on-dirichlet-processes/blob/master/2015-10-14-collapsed-gibbs-sampling-for-mixture-models.ipynb). It is part of a series on [Dirichlet Processes and Nonparametric Bayes](https://github.com/tdhopper/notes-on-dirichlet-processes).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Collapsed Gibbs Sampling for Bayesian Mixture Models (with a Nonparametric Extension)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[In an earlier notebook](http://nbviewer.ipython.org/github/tdhopper/notes-on-dirichlet-processes/blob/master/2015-09-02-fitting-a-mixture-model.ipynb), I showed how we can fit the parameters of a Bayesian mixture model using a Gibbs sampler. The sampler defines a Markov chain that, in steady state, samples from the posterior distribution of the mixture model. To move the chain forward by one step we:\n",
    "\n",
    "* Sample the cluster assignment $z_i$.\n",
    "* Sample the mixture weights $\\pi$\n",
    "* Sample the cluster means $\\mu_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that we can derive a Gibbs sampler that _just_ samples the assignments instead of the mixture weights and cluster means. This is known as a _collapsed_ Gibbs sampler. If we integrate out the cluster means $\\theta_k$ and mixture weights from the margial distribution of cluster assignment \n",
    "$$p(z_i=k \\,|\\, \n",
    "        z_{\\neg i}, \\pi,\n",
    "        \\theta_1, \\theta_2, \\theta_3, \\sigma, \\mathbf{x}, \\alpha\n",
    "        )$$ we are left with \n",
    "$$p(z_i\\,|\\, z_{\\neg i}, \\sigma, \\mathbf{x}, \\alpha).$$\n",
    "        \n",
    "By the conditional independence, we can factorize this marginal distribution \n",
    "$$\n",
    "\\begin{align}\n",
    "p(z_i=k\\,|\\, z_{\\neg i}, \\sigma, \\mathbf{x}, \\alpha)\n",
    "    &\\propto\n",
    "    p(z_i=k\\,|\\, x_i, z_{\\neg i}, \\sigma, \\mathbf{x}_{\\neg i}, \\alpha)\\\\\n",
    "    &=\n",
    "    p(z_i=k\\,|\\, z_{\\neg i}, \\sigma, \\mathbf{x}_{\\neg i}, \\alpha)\n",
    "    p(x_i \\,|\\, z, \\sigma, \\mathbf{x}_{\\neg i}, \\alpha)\\\\\n",
    "    &=\n",
    "    p(z_i=k \\,|\\, z_{\\neg i}, \\alpha) p(x_i \\,|\\, z, \\mathbf{x}_{\\neg i}, \\sigma)\\\\\n",
    "    &=\n",
    "    p(z_i=k \\,|\\, z_{\\neg i}, \\alpha)p(x_i \\,|\\, z_i=k, z_{\\neg i}, x_{\\neg_i}, \\sigma)\\\\\n",
    "    &=\n",
    "    p(z_i=k \\,|\\, z_{\\neg i}, \\alpha)p(x_i \\,|\\, \\left\\{x_j \\,|\\, z_j=k, j\\neq i\\right\\}, \\sigma).\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The two terms have intuitive explanations. $p(z_i = k \\,|\\, z_{\\neg i}, \\alpha)$ is the probability point $x_i$ will be assigned to component $k$ given the current assignments. Because we are using a symmetric Dirichlet prior, this is the predictive likelihood of a Dirichlet-categorical distribution. This is given by:\n",
    "$$p(z_i=k \\,|\\, z_{\\not i}, \\alpha)=\n",
    "    \\frac{N_k^{-i}+\\alpha / K}{N-1+\\alpha}$$\n",
    " where $N_k^{-i}=\\sum_{j\\neq i} \\delta(z_j, k)$ is the number of observation assigned to $k$ (except $x_i$). We also need to define $\\bar{x}_k^{-i}$ to be the mean of all observations assigned to component $k$ (except $x_i$).\n",
    " \n",
    "The second term is the predictive likelihood that point $x_i$ is distributed according to cluster $k$ (given the data currently in cluster $k$). For our example, we are assuming unknown cluster means are distributed according to a normal distribution with hyperparameter mean $\\lambda_1$ and variance $\\lambda_2^2$ and known cluster variance $\\sigma^2$.\n",
    "\n",
    "Thus, \n",
    "$$\n",
    "\\begin{align}\n",
    "p(x_i \\,|\\, \\left\\{x_j \\,|\\, z_j=k, j\\neq i\\right\\}, \\sigma)\n",
    "    &= \\mathcal{N}(x_i \\,|\\, \\mu_k, \\sigma_k^2+\\sigma^2)\n",
    "\\end{align}\n",
    "    $$\n",
    "where\n",
    "$$\\sigma_k^2 = \\left( \\frac{N_k^{-i}}{\\sigma^2} + \\frac{1}{\\lambda_2^2} \\right)^{-1}$$\n",
    "and\n",
    "$$\\mu_k = \\sigma_k^2 \\left( \n",
    "    \\frac{\\lambda_1}{\\lambda_2^2}+\\frac{N_k^{-i}\\cdot \\bar{x}_k^{-i}}{\\sigma^2}\n",
    "\\right).$$\n",
    "This is derived in Kevin Murphey's fantastic article [Conjugate Bayesian analysis of the Gaussian distribution](http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf). \n",
    "\n",
    "\n",
    "At each step of the collapsed sampler, we sample each $z_i$ as follows:\n",
    "\n",
    "* For each cluster $k$, compute $$f_k(x_i) =p(x_i \\,|\\, \\left\\{x_j \\,|\\, x_j=k, j\\neq i\\right\\}, \\lambda).$$ This is the predictive probability that $x_i$ is in cluster $k$ given the data currently assigned to that cluster.\n",
    "* Sample $$z_i\\sim \\frac{1}{Z_i}\\sum_{k=1}^K(N_k^{-i}+\\alpha/K)f_k(x_i)\\delta(z_i,k)$$\n",
    "\n",
    "where the normalizing constant is $Z_i=\\sum_{k=1}^K(N_k^{-i}+\\alpha/K)f_k(x_i)$.\n",
    "\n",
    "Let's write code for this Gibbs sampler!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from collections import namedtuple, Counter\n",
    "from scipy import stats\n",
    "from numpy import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(12345)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, load the same dataset we used [previously](http://nbviewer.ipython.org/github/tdhopper/notes-on-dirichlet-processes/blob/master/2015-09-02-fitting-a-mixture-model.ipynb):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHF9JREFUeJzt3X2QJPV93/H3xxySjRVpIaSOR2dB4XhQSVo9gIgiohUC\nCqscUDlVFthGtyZJuazowXYicSdVCVeljBEpx6pYSaqCfVqwBRZGFIFEUjgsWiZB4iFmEOK4HCS6\nMieFlXiSTZUTHXXf/DG9d6Ppud3enu75/Wbu86ra0nRP98xHv25+N/OZ2V5FBGZmNpt+InUAMzPr\njid5M7MZ5knezGyGeZI3M5thnuTNzGaYJ3kzsxm25iQvaYekFUmPD63/iKQnJX1b0mcG1m+X9JSk\n3ZIu6Sq0mZnVs2md+z8P/D5w8+oKSe8BLgPeFBH7Jf2dcv05wAeAc4CTgXslbYmIA50kNzOzda35\nSj4i7gdeHFr9a8DvRMT+cpsflOsvB26NiP0RsRd4Gjiv3bhmZrYRTTr5M4B/KOmbkgpJby/XnwTs\nG9huH/1X9GZmlsh6dc3h9jk2Is6XdC5wG3D6Ybb1NRPMzBJqMsnvA+4AiIiHJR2QdDzwXeDUge1O\nKdf9GEme+M3MGogIbXSfJnXNncCFAJK2AK+KiOeAu4ArJL1K0mn0a52HDhM0u59rr702eQZncqYj\nMZcz1ftpas1X8pJuBd4N/G1JzwCfBnYAO8qvVf4I+GA5ce+SdBuwC3gF+FCMk2zC9u7dmzpChTPV\n40z15ZjLmbq15iQfEVce5q6rDrP9dcB144YyM7N2+DdeS0tLS6kjVDhTPc5UX465nKlbmnSjImma\nWhwzsyxIIib0wetMKooidYQKZ6rHmerLMZczdcuTvJnZDHNdY2Y2BVzXmJlZhSf5Uo4dnDPV40z1\n5ZjLmbrV5LIGZma1SRtuGH6M693xuJM3s071J/mm/83Lk3zJnbyZmVV4ki/l2ME5Uz3OVF+euYrU\nASryHKdmPMmbmc0wd/Jm1il38u1wJ29mZhWe5Es5dnDOVI8z1ZdnriJ1gIo8x6kZT/JmZjPMnbyZ\ndcqdfDvcyZuZWYUn+VKOHZwz1eNM9eWZq0gdoCLPcWpmzUle0g5JK+Uf7R6+719IOiDpuIF12yU9\nJWm3pEu6CGxmZvWt2clLugB4Gbg5It44sP5U4EbgTOBtEfGCpHOAW4BzgZOBe4EtEXFg6DHdyZsd\nQdzJt6OTTj4i7gdeHHHXvwE+MbTucuDWiNgfEXuBp4HzNhrIzMzas+FOXtLlwL6I+NbQXScB+waW\n99F/RT8VcuzgnKkeZ6ovz1xF6gAVeY5TMxu6nrykY4BPAhcPrl5jF7/PMjNLaKN/NOT1wDzwWPmH\nAE4B/oekdwDfBU4d2PaUcl3F0tIS8/PzAMzNzbGwsMDi4iJw6F9QLy+yuLiYVZ5VRVFkk2f4FVcu\neXJdXl2X6vgcetW+OLR8uPvzPN8msVwUBcvLywAH58sm1v1lKEnzwN2DH7wO3Pcdqh+8nsehD17/\n3vCnrP7g1ezI4g9e29HJB6+SbgUeALZIekbSrwxtcnD0I2IXcBuwC/gK8KFpms2rrzjSc6Z6nKm+\nPHMVqQNU5DlOzaxZ10TElevcf/rQ8nXAdS3kMjOzFvjaNWbWKdc17fC1a8zMrMKTfCnHDs6Z6nGm\n+vLMVaQOUJHnODXjSd7MbIa5kzezTrmTb4c7eTMzq/AkX8qxg3OmepypvjxzFakDVOQ5Ts14kjcz\nm2Hu5M2sU+7k2+FO3szMKjzJl3Ls4JypHmeqL89cReoAFXmOUzOe5M3MZpg7eTPrlDv5driTNzOz\nCk/ypRw7OGeqx5nqyzNXkTpARZ7j1IwneTOzGeZO3sw65U6+He7kzcyswpN8KccOzpnqcab68sxV\npA5Qkec4NbPeH/LeIWlF0uMD6/61pCclPSbpDkmvG7hvu6SnJO2WdEmXwc3MbH1rdvKSLgBeBm6O\niDeW6y4G/iwiDki6HiAitkk6B7gFOBc4GbgX2BIRB4Ye05282RHEnXw7OunkI+J+4MWhdTsHJu4H\ngVPK25cDt0bE/ojYCzwNnLfRQGZm1p5xO/mrgS+Xt08C9g3ct4/+K/qpkGMH50z1OFN9eeYqUgeo\nyHOcmtnUdEdJnwJ+FBG3rLHZyPdZS0tLzM/PAzA3N8fCwgKLi4vAocGd9PKqVM8/Lcu9Xi+rPEVR\n0Ov1ssozKJc8qY/fIavLiwPLvaHlwfv7j3EkHr+iKFheXgY4OF82se735CXNA3evdvLluiXgnwHv\njYj/W67bBhAR15fLXwWujYgHhx7PnbzZEcSdfDsm9j15SZcCHwcuX53gS3cBV0h6laTTgDOAhzb6\n+GZm1p71vkJ5K/AAcKakZyRdDfw+8Bpgp6RHJf17gIjYBdwG7AK+Anxoml6yV99WpudM9ThTfXnm\nKlIHqMhznJpZs5OPiCtHrN6xxvbXAdeNG8rMzNrha9eYWafcybejaSff+Ns1ZnZk6U/WNm187ZpS\njh2cM9XjTPWNnysa/KybasxM7cv1+DXhSd7MbIa5kzezWpp36+7k2+DryZuZWYUn+VKOHZwz1eNM\n9eWZq0gdoCLPcWrGk7yZ2QxzJ29mtbiTT8udvJmZVXiSL+XYwTlTPc5UX565itQBKvIcp2Y8yZuZ\nzTB38mZWS6pOvqlZm2d87Rozm1FN/2ExcF1zUI4dnDPV40z15ZmrSB2gIs9xasaTvJnZDHMnb2a1\npOvkmz3nrM0z/p68mZlVeJIv5djBOVM9zlRfnrmK1AEq8hynZtb7Q947JK1Ienxg3XGSdkraI+ke\nSXMD922X9JSk3ZIu6TK4mZmtb81OXtIFwMvAzRHxxnLdDcBzEXGDpGuAYyNim6RzgFuAc4GTgXuB\nLRFxYOgx3cmbTSF38ml10slHxP3Ai0OrLwNuKm/fBLy/vH05cGtE7I+IvcDTwHkbDWRmZu1p0slv\njoiV8vYKsLm8fRKwb2C7ffRf0U+FHDs4Z6rHmerLM1eROkBFnuPUzFi/8RoRIWmt90Qj71taWmJ+\nfh6Aubk5FhYWWFxcBA4N7qSXV6V6/mlZ7vV6WeUpioJer5dVnkG55Gnr+B2akDe6zBr399bYf3Vd\ns+eb5uNXFAXLy8sAB+fLJtb9nrykeeDugU5+N7AYEc9KOhG4LyLOkrQNICKuL7f7KnBtRDw49Hju\n5M2mkDv5tCb5Pfm7gK3l7a3AnQPrr5D0KkmnAWcADzV4fDMza8l6X6G8FXgAOFPSM5J+BbgeuFjS\nHuDCcpmI2AXcBuwCvgJ8aJpesg+/TcuBM9XjTPXlmatIHaAiz3FqZs1OPiKuPMxdFx1m++uA68YN\nZWZm7fC1a8ysFnfyafnaNWZmVuFJvjSpDk7SWD+p5dhVOlN9eeYqUgeoyHOcmvEkn0TU/LlvaNnM\nbGPcyU9Y814TZrFntOnhTj4td/JmZlbhSb6UZwdXpA5QkeM4OVN9eeYqUgeoyHOcmvEkb2Y2w9zJ\nT5g7eZtW7uTTcidvZmYVnuRLeXZwReoAFTmOkzPVl2euInWAijzHqRlP8mZmM8yd/IS5k7dp5U4+\nLXfyZmZW4Um+lGcHV6QOUJHjODlTfXnmKlIHqMhznJrxJG9mNsPcyU+YO3mbVu7k03Inb2ZmFZ7k\nS3l2cEXqABU5jpMz1ZdnriJ1gIo8x6mZxpO8pO2SnpD0uKRbJL1a0nGSdkraI+keSXNthjUzs41p\n1MlLmge+BpwdEf9P0heBLwNvAJ6LiBskXQMcGxHbhvZ1J+9O3qaQO/m0Jt3J/xWwHzhG0ibgGOB7\nwGXATeU2NwHvb/j4ZmbWgkaTfES8APwu8Jf0J/eXImInsDkiVsrNVoDNraScgDw7uCJ1gIocx8mZ\n6sszV5E6QEWe49TMpiY7SXo98OvAPPBD4E8l/fLgNhERkka+X1paWmJ+fh6Aubk5FhYWWFxcBA4N\n7qSXV3X9fOWzAIsDt6m9nGp8Vpd7vV7S5x+13Ov1ssozKJc8bR2/jZ6v1Ql81P29NfZfXdfs+ab5\n+BVFwfLyMsDB+bKJpp38B4CLI+KflstXAecDFwLviYhnJZ0I3BcRZw3t607enbxNIXfyaU26k98N\nnC/pp9Q/8hcBu4C7ga3lNluBOxs+vpmZtaBpJ/8YcDPwCPCtcvV/BK4HLpa0h/6r+uvbCDkJw2/T\n8lCkDlCR4zg5U3155ipSB6jIc5yaadTJA0TEDcANQ6tfoP+q3szMMuBr10yYO3mbVu7k0/K1a8zM\nrMKTfCnPDq5IHaAix3FypvryzFWkDlCR5zg140nezGyGuZOfMHfyNq3cyaflTt7MzCo8yZfy7OCK\n1AEqchwnZ6ovz1xF6gAVeY5TM57kzcxmmDv5CXMnb9PKnXxa7uTNzKzCk3wpzw6uSB2gIsdxcqb6\n8sxVpA5Qkec4NeNJ3sxshrmTnzB38jat3Mmn5U7ezMwqPMmX8uzgitQBKnIcJ2eqL89cReoAFXmO\nUzONrydvafTfMm/crL11NbN63MlP2LidvPtJS8WdfFru5M3MrMKTfCnPDq5IHaAix3FypvryzFWk\nDlCR5zg103iSlzQn6XZJT0raJekdko6TtFPSHkn3SJprM6yZmW1M405e0k3A1yNih6RNwE8DnwKe\ni4gbJF0DHBsR24b2cyfvftKmkDv5tJp28o0meUmvAx6NiNOH1u8G3h0RK5JOAIqIOGtoG0/yPuFt\nCnmST2vSH7yeBvxA0ucl/YWkGyX9NLA5IlbKbVaAzQ0ff+Ly7OCK1AEqchwnZ6ovz1xF6gAVeY5T\nM02/J78JeCvw4Yh4WNJngR+rZSIiJI38p3RpaYn5+XkA5ubmWFhYYHFxETg0uJNeXtX185XPAiwO\n3GYCyzTKO7zc6/XG2r+L5V6vl1WeQbnkaev4jXv+jb6/t8b+q+uaPd80H7+iKFheXgY4OF820bSu\nOQH4RkScVi6/C9gOnA68JyKelXQicJ/rmh/nusamleuatCZa10TEs8AzkraUqy4CngDuBraW67YC\ndzZ5fDMza8c435P/CPAFSY8BbwJ+G7geuFjSHuDCcnkq5NnBFakDVOQ4Ts5UX565itQBKvIcp2Ya\nX7smIh4Dzh1x10XN45g1vz7Pqll7m242Dl+7ZsLcya/P19zPkzv5tJp28r4KZUPjvto0M5sEX7um\n1KyDiwY/G0rVIFO38uwqi9QBKvIcp1xzFakDVOQ5Ts14kjczm2Hu5BtyP9kdd/J58jmfljt5s5L/\nepbZIa5rSnl2cEXqABXTMU5df1ZSI1GW45RrriJ1gIo8x6kZT/JmZjPMnXxD7ie7498lyJPP+bT8\nN17NzKzCk3wpzw6uSB2gwuNUT57j1H812PSnO0WHj91MrsevCX+7xuyIM051YtPGnXxD7ie7406+\nO6nG1sdzfO7kzcyswpN8Kc8OrkgdoMLjVE+e45SrInWAilk6fp7kzcxmmDv5htzJd8edfHfcyU8v\nd/JmZlbhSb6UZwdXpA5Q4XGqJ89xylWROkDFLB2/sb4nL+ko4BFgX0T8I0nHAV8E/i6wF/iFiHhp\n7JQ2tnF+maXJ217/5SyzPIzVyUv6TeBtwN+KiMsk3QA8FxE3SLoGODYitg3t405+yjrR5pP8pP9/\njrPv7HW4oxxJnfw4cjwXJt7JSzoFeB/wBxwa0cuAm8rbNwHvb/r4ZmbjaXLJ6fwm93GN08n/HvBx\n4MDAus0RsVLeXgE2j/H4E5VnB1ekDlDhcaonz3HKVZE6QMUsHb9GnbyknwO+HxGPSloctU1EhKSR\n/ywuLS0xPz8PwNzcHAsLCywu9h9mdXAnvbxqo9sfOkEXay6vrqu7fVvLrHP/4Zb7/59X///3er3+\nvZ2NT5P9e4w/vuVSovNpUstlKvI6/3ocfv/VdW0+X/39Ux6voihYXl4GODhfNtGok5d0HXAV8Arw\nk8BrgTuAc4HFiHhW0onAfRFx1tC+7uSnrBN1Jz87jqxOfrb+TvBEO/mI+GREnBoRpwFXAF+LiKuA\nu4Ct5WZbgTubPL6ZmbWjre/Jr/6zdz1wsaQ9wIXl8lTIs4MrUgeo8DjVk+c45apIHaBilo7f2NeT\nj4ivA18vb78AXDTuY5qZWTt87ZqG3Mmvs5c7+Sy5k6+3b47ngq9dY2ZmFZ7kS3l2cEXqABUep3ry\nHKdcFakDVMzS8fMkb2Y2w9zJN+ROfp293MlnyZ18vX1zPBfcyZuZWYUn+VKeHVyROkCFx6mePMcp\nV0XqABWzdPw8yZuZzbAjupMf/w9buJ887F7u5LPkTr7evjmeC007+bF/43X6jXPymZnlzXXNQUXq\nACMUqQNU5NlVFqkDVOQ5TrkqUgeomKXj50nezGyGuZOfsq5wWvpJd/J58jlfb98czwV/T97MzCo8\nyR9UpA4wQpE6QEWeXWWROkBFnuOUqyJ1gIpZOn7+do1ZaZyv1Ob49t4M3MkzbV3htPST09jJz1qH\nO4rP+Xr75ng8/T1568z4vzRmZqk06uQlnSrpPklPSPq2pI+W64+TtFPSHkn3SJprN26XitQBRihS\nByjFwM99Q8uH+5mkYsLPt75Z6nS7V6QOUDFLx6/pB6/7gd+IiDcA5wP/XNLZwDZgZ0RsAf6sXDYz\ns0Ra6eQl3Ql8rvx5d0SsSDoBKCLirKFt3cn7OTPcd/Y63FF8ztfbN8fjmex78pLmgbcADwKbI2Kl\nvGsF2Dzu45uZWXNjTfKSXgN8CfhYRPz14H3ly/X8/jk8rCJ1gBGK1AFGKFIHGKFIHaBiljrd7hWp\nA1TM0vFr/O0aSUfTn+D/KCLuLFevSDohIp6VdCLw/VH7Li0tMT8/D8Dc3BwLCwssLi4ChwZ3UsvV\nE2x1efj+cbdfXV5dV3f7tpZZ5/7DLa+uW13udfx8TfbvMf74buT5Dr//8OQw6fN5/fN9NfN6/3/a\nXmaN+3scfv/VdW0+X/39Ux6voihYXl4GODhfNtGok1e/2LsJeD4ifmNg/Q3lus9I2gbMRcS2oX3d\nyfs5M9x39jrcUXzO19s3x+PZtJNvOsm/C/hz4FscGsntwEPAbcDPAHuBX4iIl4b29STv58xw39mb\nFEbxOV9v3xyP50Q/eI2I/xYRPxERCxHxlvLnqxHxQkRcFBFbIuKS4Qk+b0XqACMUqQOMUKQOMEKR\nOkDFLHW63StSB6iYpePnC5SZmc0wX7tmyt5G+jm72nf23t6P4nO+3r45Hk9fT97MzCo8yR9UpA4w\nQpE6wAhF6gAjFKkDVMxSp9u9InWACkmNfnLkSd7MrKLuhfhSXZSvPnfyU9YV+jm72nf2OtxRfM53\n+5xdngfu5M3MrMKT/EFF6gAjFKkDjFCkDjBCkTpAhTv5jShSBxihSB2gNZ7kzcxmmDv5I6QrPDKe\nc5x93cnX2LvhvkfOc+bYyftvvJpNoVy/rmf5cV1zUJE6wAhF6gAjFKkDjFCkDlAxmU5+o1/xy/Xd\nRpE6wAhF6gCt8SRvZjbD3MkfIV3hkfGc4+w77nM20/S/hebn7jSO7fQ8pzt5s5nVdCIy65brmoOK\n1AFGKFIHGKFIHWCEInWAEYrUAaZIkTrACEXqAK2Z+lfyn/70v+KBB3rrb2hmdgSa+k7+/PMv5cEH\n3w68ZYN7fgf4OEdKV3hkPOc4+05Xh+tOPs/ndCffmXcBl25wH7/6N7PZ13onL+lSSbslPSXpmrYf\nvztF6gAjFKkDjFCkDjBCkTrACEXqAFOkSB1ghCJ1gNa0OslLOgr4HP2X1ecAV0o6u83n6E6Or+yd\nqR5nmm45jlWOmZpp+5X8ecDTEbE3IvYDfwJc3vJzdOSl1AFGcKZ6nGm65ThWOWZqpu1O/mTgmYHl\nfcA7Wn4Os5nha9BY19qe5Cf+q6ybNsExx1zLpk3/bkP7HTjwQ15+eXDN3jZjtWRv6gAj7E0dYIS9\nqQOMsLfmdpP/Ldv87E0dYIS9qQO0ptWvUEo6H/itiLi0XN4OHIiIzwxsk8c1DczMpkyTr1C2Pclv\nAv4n8F7ge8BDwJUR8WRrT2JmZrW1WtdExCuSPgz8V+Ao4A89wZuZpTPx33g1M7PJ6fwCZZKOk7RT\n0h5J90iaO8x22yU9IelxSbdIenUGmeYk3S7pSUm7ys8cOlM3V7ntUZIelXR36kySTpV0X3n8vi3p\nox1lWfcX7ST92/L+xyRt9FoXrWeS9Etllm9J+u+S3pQ608B250p6RdLP55BJ0mJ5Tn9bUtF1pjq5\nJB0v6auSemWupY7z7JC0IunxNbbZ2DkeEZ3+ADcAnyhvXwNcP2KbeeB/A68ul78IbE2ZqbzvJuDq\n8vYm4HWpx2pg298EvgDclToTcAKwUN5+Df3PZc5uOcdRwNPluXI0/d9WOXtom/cBXy5vvwP4Zsdj\nUyfT3189b+j/kmDyTAPbfQ34z8A/Tp0JmAOeAE4pl4/vMtMGcv0W8DurmYDngU0dZrqA/oW4Hj/M\n/Rs+xydxqeHL6E+WlP/7/hHb/BWwHzim/PD2GOC7KTNJeh1wQUTsgP7nDRHxww4z1cpVZjuF/sH+\nA7r/Lt26mSLi2YjolbdfBp4ETmo5R51ftDuYNSIeBOYkbW45x4YyRcQ3Bs6bB4FTOsxTK1PpI8Dt\nwA86zlM30y8CX4qIfQAR8Vwmuf4P8Nry9muB5yPila4CRcT9wItrbLLhc3wSk/zmiFgpb68AlUAR\n8QLwu8Bf0v9WzksRcW/KTMBpwA8kfV7SX0i6UdIxHWaqmwvg9+hfQvNAx3k2kgkASfP0X4k82HKO\nUb9od3KNbbqcVOtkGvRPgC93mAdqZJJ0Mv3J7D+Uq7r+YK7OOJ0BHFfWfo9IuqrjTHVz3Qi8QdL3\ngMeAj00g11o2fI638u0aSTvpv2Uf9qnBhYiIUd+Tl/R64Nfpv236IfCnkn4pIr6QKhP9sXkr8OGI\neFjSZ4FtwKebZmojl6SfA74fEY9KWhwnS1uZBh7nNfRfHX6sfEXfproT0fA7my4nsNqPLek9wNXA\nP+guDlAv02eBbeXxFN2/G6yT6Wj6/729l/47+W9I+mZEPJU41yeBXkQslvPUTklvjoi/7jDXejZ0\njrcyyUfExYdN0/8Q4YSIeFbSicD3R2z2duCBiHi+3OcO4J30O+dUmfYB+yLi4XL5dvqT/FhayPVO\n4DJJ7wN+EnitpJsj4oMJMyHpaOBLwB9HxJ1Ns6zhu8CpA8un0j9Ga21zCt3WfnUyUX7YeiNwaUSs\n9VZ8UpneBvxJf37neOBnJe2PiLsSZnoGeC4i/gb4G0l/DrwZ6HKSr5PrncBvA0TE/5L0HeBM4JEO\nc61lw+f4JOqau4Ct5e2twKgJYDdwvqSfKl9ZXATsSpkpIp4FnpG0pVx1Ef0PhrpUJ9cnI+LUiDgN\nuAL42jgTfBuZymP2h8CuiPhsRzkeAc6QNC/pVcAHymzDWT9YZjqffu23QnfWzSTpZ4A7gF+OiKc7\nzFI7U0ScHhGnlefQ7cCvdTjB18oE/CfgXep/a+wY+h8qdjkH1M21m/5/+5Td95n0vySSysbP8a4+\nJR74NPg44F5gD3APMFeuPwn4LwPbfYL+JPo4/Q8Wjs4g05uBh+l3cXfQ/bdrauUa2P7ddP/tmnUz\n0f+rLQfofzvh0fLn0g6y/Cz9b+48DWwv1/0q8KsD23yuvP8x4K0TOL/XzET/w/HnB8blodSZhrb9\nPPDzOWQC/uXAHPDRrjPVPH7HA3eX59PjwC92nOdW+p9L/oj+u5urxz3H/ctQZmYzbBJ1jZmZJeJJ\n3sxshnmSNzObYZ7kzcxmmCd5M7MZ5knezGyGeZI3M5thnuTNzGbY/wcVN926vQHvhAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10af51050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = pd.Series.from_csv(\"clusters.csv\")\n",
    "_=data.hist(bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we want to define a state object and a function for updating the sufficient statistics of the state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SuffStat = namedtuple('SuffStat', 'theta N')\n",
    "\n",
    "def initial_state(num_clusters=3, alpha=1.0):\n",
    "    cluster_ids = range(num_clusters)\n",
    "\n",
    "    state = {\n",
    "        'cluster_ids_': cluster_ids,\n",
    "        'data_': data,\n",
    "        'num_clusters_': num_clusters,\n",
    "        'cluster_variance_': .01,\n",
    "        'alpha_': alpha,\n",
    "        'hyperparameters_': {\n",
    "            \"mean\": 0,\n",
    "            \"variance\": 1,\n",
    "        },\n",
    "        'suffstats': {cid: None for cid in cluster_ids},\n",
    "        'assignment': [random.choice(cluster_ids) for _ in data],\n",
    "        'pi': {cid: alpha / num_clusters for cid in cluster_ids},\n",
    "    }\n",
    "    update_suffstats(state)\n",
    "    return state\n",
    "\n",
    "def update_suffstats(state):\n",
    "    for cluster_id, N in Counter(state['assignment']).iteritems():\n",
    "        points_in_cluster = [x \n",
    "            for x, cid in zip(state['data_'], state['assignment'])\n",
    "            if cid == cluster_id\n",
    "        ]\n",
    "        mean = np.array(points_in_cluster).mean()\n",
    "        \n",
    "        state['suffstats'][cluster_id] = SuffStat(mean, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define functions to compute the two terms of our marginal distribution over cluster assignments (as we derived above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def log_predictive_likelihood(data_id, cluster_id, state):\n",
    "    \"\"\"Predictive likelihood of the data at data_id is generated\n",
    "    by cluster_id given the currenbt state.\n",
    "    \n",
    "    From Section 2.4 of \n",
    "    http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf\n",
    "    \"\"\"\n",
    "    ss = state['suffstats'][cluster_id]\n",
    "    hp_mean = state['hyperparameters_']['mean']\n",
    "    hp_var = state['hyperparameters_']['variance']\n",
    "    param_var = state['cluster_variance_']\n",
    "    x = state['data_'][data_id]\n",
    "    return _log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x)\n",
    "\n",
    "\n",
    "def _log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x):\n",
    "    posterior_sigma2 = 1 / (ss.N * 1. / param_var + 1. / hp_var)\n",
    "    predictive_mu = posterior_sigma2 * (hp_mean * 1. / hp_var + ss.N * ss.theta * 1. / param_var)\n",
    "    predictive_sigma2 = param_var + posterior_sigma2\n",
    "    predictive_sd = np.sqrt(predictive_sigma2)\n",
    "    return stats.norm(predictive_mu, predictive_sd).logpdf(x)\n",
    "\n",
    "\n",
    "def log_cluster_assign_score(cluster_id, state):\n",
    "    \"\"\"Log-likelihood that a new point generated will\n",
    "    be assigned to cluster_id given the current state.\n",
    "    \"\"\"\n",
    "    current_cluster_size = state['suffstats'][cluster_id].N\n",
    "    num_clusters = state['num_clusters_']\n",
    "    alpha = state['alpha_']\n",
    "    return np.log(current_cluster_size + alpha * 1. / num_clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given these two functions, we can compute the posterior probability distribution for assignment of a given datapoint. This is the core of our collapsed Gibbs sampler.\n",
    "\n",
    "To simplify the computation of things like $N_k^{-i}$ (where we remove point $i$ from the summary statistics), we create two simple functions to add and remove a point from the summary statistics for a given cluster. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cluster_assignment_distribution(data_id, state):\n",
    "    \"\"\"Compute the marginal distribution of cluster assignment\n",
    "    for each cluster.\n",
    "    \"\"\"\n",
    "    scores = {}\n",
    "    for cid in state['suffstats'].keys():\n",
    "        scores[cid] = log_predictive_likelihood(data_id, cid, state)\n",
    "        scores[cid] += log_cluster_assign_score(cid, state)\n",
    "    scores = {cid: np.exp(score) for cid, score in scores.iteritems()}\n",
    "    normalization = 1.0/sum(scores.values())\n",
    "    scores = {cid: score*normalization for cid, score in scores.iteritems()}\n",
    "    return scores\n",
    "\n",
    "def add_datapoint_to_suffstats(x, ss):\n",
    "    \"\"\"Add datapoint to sufficient stats for normal component\n",
    "    \"\"\"\n",
    "    return SuffStat((ss.theta*(ss.N)+x)/(ss.N+1), ss.N+1)\n",
    "\n",
    "\n",
    "def remove_datapoint_from_suffstats(x, ss):\n",
    "    \"\"\"Remove datapoint from sufficient stats for normal component\n",
    "    \"\"\"\n",
    "    return SuffStat((ss.theta*(ss.N)-x*1.0)/(ss.N-1), ss.N-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we're ready to create a function that takes a Gibbs step on the state. For each datapoint, it\n",
    "1. Removes the datapoint from its current cluster.\n",
    "2. Computes the posterior probability of the point being assigned to each cluster (given the other current assignments).\n",
    "3. Assigns the datapoint to a cluster sampled from this probability distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gibbs_step(state):\n",
    "    pairs = zip(state['data_'], state['assignment'])\n",
    "    for data_id, (datapoint, cid) in enumerate(pairs):\n",
    "\n",
    "        state['suffstats'][cid] = remove_datapoint_from_suffstats(datapoint, \n",
    "                                                                  state['suffstats'][cid])\n",
    "        scores = cluster_assignment_distribution(data_id, state).items()\n",
    "        labels, scores = zip(*scores)\n",
    "        cid = random.choice(labels, p=scores)\n",
    "        state['assignment'][data_id] = cid\n",
    "        state['suffstats'][cid] = add_datapoint_to_suffstats(state['data_'][data_id], state['suffstats'][cid])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's our old function to plot the assignments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_clusters(state):\n",
    "    gby = pd.DataFrame({\n",
    "            'data': state['data_'], \n",
    "            'assignment': state['assignment']}\n",
    "        ).groupby(by='assignment')['data']\n",
    "    hist_data = [gby.get_group(cid).tolist() \n",
    "                 for cid in gby.groups.keys()]\n",
    "    plt.hist(hist_data, \n",
    "             bins=20,\n",
    "             histtype='stepfilled', alpha=.5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomly assign the datapoints to a cluster to start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHS9JREFUeJzt3X1wXHd97/H3d7W7erItWU4iO3YenEASYAjhIQFCuAhw\nWkNJmsKdXCgtviXTYcLl4cLcNoEyxZ4OLWWG28xte+kMBcY35bGEhtD0kpiQDYQb8tA4IYntOLZj\nxw/Ryras56fd1ff+oROjClna32qPdnX0ec14cnZ1frufSKuPzv7O2XPM3RERkaUtVesAIiKycCpz\nEZEEUJmLiCSAylxEJAFU5iIiCaAyFxFJgHnL3MwuNbOd0/71m9knzKzDzHaY2V4zu9fM2hcjsIiI\n/CYLOc7czFLAUeAq4OPACXf/kpndAqx291vjiSkiInMJnWbZBOxz98PA9cD26P7twA3VDCYiIuUL\nLfP3A9+OljvdPR8t54HOqqUSEZEgZZe5mWWB64B/nvk1n5qr0XkBRERqJB2w7ruAf3f349HtvJmt\ndfduM1sH9MwcYGYqeBGRCri7hawfMs3yAX49xQJwF7AlWt4C3HmGQHX17/Of/3zNMyhTsnIpkzJV\n+18lyipzM2tlaufnD6bd/UXgWjPbC7wjui0iIjVQ1jSLuw8DZ824r5epghcRkRpbdp8A7erqqnWE\n36BM5avHXMpUHmWKV9CHhoIf3MzjfHwRkSQyMzzGHaAiIlKnVOYiIgmgMhcRSQCVuYhIAqjMRUQS\nQGUuIpIAKnMRkQRQmYuIJIDKXEQkAVTmIiIJoDIXEUkAlbmISAKozEVEEkBlLiKSACpzEZEEUJmL\niCSAylxEJAFU5iIiCaAyFxFJAJW5iEgCqMxFRBIgXc5KZtYO/CPwKsCBPwKeA74LXAAcBG509754\nYopIvXvxxRdx94rGtrS00N7eXuVEy4uV8803s+3AA+7+dTNLA63AnwEn3P1LZnYLsNrdb50xziv9\n4YrI0vKpz32KfusnlQp7wz8+Os57r3kvv3f978WUbOkxM9zdQsbMu2VuZm3AW919C4C7F4F+M7se\neFu02nYgB9w664OISOKNFcY4953nks6W9Yb/tKO7jzI5ORlTquWjnD+hG4HjZvYNM3vczL5qZq1A\np7vno3XyQGdsKUVEZE7l/AlNA68DPubuj5rZbczYAnd3N7NZ51O2bt16ermrq4uurq6Kw4qIJFEu\nlyOXyy3oMeadMzeztcBD7r4xun0N8BngIuDt7t5tZuuA+939shljNWcuskzcfMvNrHnnmoqmWTav\n28z7bnhfTMmWnkrmzOedZnH3buCwmV0S3bUJeAb4EbAlum8LcGfIE4uISPWU+yf048A3zSwL7Gfq\n0MQG4HtmdhPRoYmxJBQRkXmVVebu/iRw5Sxf2lTdOCIiUgl9AlREJAFU5iIiCaAyFxFJAJW5iEgC\nqMxFRBJAZS4ikgAqcxGRBFCZi4gkgMpcRCQBVOYiIgmgMhcRSQCVuYhIAqjMRUQSQGUuIpIAKnMR\nkQRQmYuIJIDKXEQkAVTmIiIJoDIXEUkAlbmISAKozEVEEkBlLiKSACpzEZEESJezkpkdBAaAElBw\n96vMrAP4LnABcBC40d37YsopIiJzKHfL3IEud3+tu18V3XcrsMPdLwHui26LiEgNhEyz2Izb1wPb\no+XtwA1VSSQiIsFCtsx/YmaPmdkfR/d1uns+Ws4DnVVPJyIiZSlrzhx4i7u/aGZnAzvMbM/0L7q7\nm5nPNnDr1q2nl7u6uujq6qowqohIMuVyOXK53IIew9xn7eAzDzD7PDAE/DFT8+jdZrYOuN/dL5ux\nroc+vogsTTffcjNr3rmGdLbcbcQpR3cfZfO6zbzvhvfFlGzpMTPcfebU9pzmnWYxsxYzWxkttwK/\nBTwF3AVsiVbbAtwZFldERKqlnD+hncC/mNlL63/T3e81s8eA75nZTUSHJsaWUkTq3qMPPos/04al\ngjYoGe4b5uz/9IS2zBdo3jJ39+eBK2a5vxfYFEcoEVl6RocbadvwR6SzTUHjBk7sZGJiMKZUy0fY\n5JaIyBxSDVlSDdmgMWYNMaVZXvRxfhGRBFCZi4gkgMpcRCQBVOYiIgmgMhcRSQCVuYhIAqjMRUQS\nQGUuIpIAKnMRkQRQmYuIJIDKXEQkAVTmIiIJoDIXEUkAlbmISAKozEVEEkBlLiKSACpzEZEEUJmL\niCSAylxEJAF0DVARqYrBwQFKx3tIpxuDxg0NDNDT0xNTquVDZS4iVTE2PMTqvn7SmbALOjM0xFB/\nXzyhlhGVuYhUTWs2SzYTtmV+PKXZ3mooq8zNrAF4DDji7teZWQfwXeAC4CBwo7vrT6vIEnf06FGO\nHDlS2WCvbhYJU+6W+SeBXcDK6PatwA53/5KZ3RLdvjWGfCKyiA4cOMBtP7qNtnPagsdOtk5iFkMo\nKcu8ZW5mG4B3A18APh3dfT3wtmh5O5BDZS6SCCvXruT8150fPK5x1ZOY2rxmytky/xvgT4BV0+7r\ndPd8tJwHOqsdTESWj4mJCQ4dOlTR2PXr15NOa/ffnN8BM3sP0OPuO82sa7Z13N3N7IyzZVu3bj29\n3NXVRVfXrA8jIstUKtXIC4fX8Rd/8dPgsWZH+PKXP8mqVavmX7mO5XI5crncgh5jvj9nVwPXm9m7\ngSZglZndDuTNbK27d5vZOuCMB4lOL3MRkZkam9bStuptbNjwW8FjDx/+nzEkWnwzN3S3bdsW/Bhz\nHhPk7p919/PcfSPwfuCn7v6HwF3Almi1LcCdwc8sIiJVE3qA50vTKV8ErjWzvcA7otsiIlIjZe81\ncPcHgAei5V5gU1yhREQkjD56JSKSADqeR0RO2717N4/96zGee2Q4eGxhtOnXHysM1JfP88sdO8LH\njTzO4ODgkj+apRpU5iJy2sTEBCP9b6Ct883BY5va9mEWXinN6QwXZjJcGjwS7h8awl3nEQCVuYjM\nkGrIkmkM38RuSLdW9AnQlBnZhgaaM5ngsTp/wK+pzEXktO7ubgb7h+k5dix4rJdKUEkhS1WozEXk\ntKGhIRgZo3V4KHhsa7oBbSfXjspcRP6DhpTRGnqBCak5HZooIpIAKnMRkQRQmYuIJIDmzGNy8OBB\nxsfHKxrb2trKhg0bqpxIRJJMZR6T2++4nX39+0hnw77F46PjXH3x1Xz0wx+NKZmIJJHKPCbFUpGz\nrzibFR0rgsadPHKS0lgpplQiklSaMxcRSQCVuYhIAqjMRUQSQGUuIpIAKnMRkQRQmYuIJIDKXEQk\nAXSceUx2PbmP47vHaGxpDBo30jdM4SydSFREwqjMY5LPF+izd9C8anXQuIHjRzgxuS+mVCKSVCrz\nGGVbzqJ5ZWfQmNHB8AvpiojMOWduZk1m9rCZPWFmu8zsr6L7O8xsh5ntNbN7zax9ceKKiMhs5ixz\ndx8D3u7uVwCXA283s2uAW4Ed7n4JcF90W0REamTeo1ncfSRazAINwCngemB7dP924IZY0omISFnm\nnTM3sxTwOHAx8BV3f8bMOt09H62SB8ImhkVEphktFMgPhV9EemR8nGKxGEOipWfeMnf3SeAKM2sD\n7jGzt8/4upuZn2n81q1bTy93dXXR1dVVcVgRSZ5MJsXhhgGOjA4Ej+0e6WN0dDSGVIsrl8uRy+UW\n9BhlH83i7v1mdjfweiBvZmvdvdvM1gE9Zxo3vcxFRGZqbs7Q3JypaOzhE1UOUyMzN3S3bdsW/Bjz\nHc1y1ktHqphZM3AtsBO4C9gSrbYFuDP4mUVEpGrm2zJfB2yP5s1TwO3ufp+Z7QS+Z2Y3AQeBG+ON\nKSIic5mzzN39KeB1s9zfC2yKK5SIiITRibZERBJAH+ePyYtHjjCcfZ5s/lTQuNH+F2keeD6mVCKS\nVCrzmHipxNnpBtobw86a2G0pvFSKKZWIJJXKPFaGmU5nKyLx05y5iEgCqMxFRBJAZS4ikgAqcxGR\nBFCZi4gkgMpcRCQBdGhiHSoUCjz66KMVjb3ssstYuXJllRNV38jICM8880zF4y+66CLWrFlTxUQi\nS5vKvM40pI2BxiG+kvtK8NiJ/ARf+PgXlkSZDw4O8rff+VvSa8NfgqM9o3z2Q59VmYtMozKvM+l0\nikxnC+dfeX7w2MO5wzEkik+mOcN5V54XPO7wQ0vr/1NkMajM61CpWKK/pz943FD/EOPj4zEkEpF6\npzKvM5bKMJDfSG57eCkPHxvk+H8+ziWXXBJDMhGpZyrzOtPY2EFr8+/QdvaG4LHDx7bHkEhElgId\nmigikgAqcxGRBFCZi4gkgMpcRCQBtAN0Dnv37uWe22+HYjF4bHF0FOr/szsikhAq8zmMj4/TcegQ\nv712bfDYb2ezjJve+IjI4lCZz6MxneaslpbgcemGBiZ0xTgRWSTzbjqa2Xlmdr+ZPWNmT5vZJ6L7\nO8xsh5ntNbN7zaw9/rgiIjKbcuYBCsCn3P1VwJuA/2ZmrwBuBXa4+yXAfdFtERGpgXnL3N273f2J\naHkI2A2sB64HXvrI4XbghrhCiojI3ILmzM3sQuC1wMNAp7vnoy/lgc6qJhOZw549e+jr6wsaMz4+\nzv79+3nlK19Z0XO+7GUvY926dRWNXUzuzoMPPljR2N7eXqC1uoFkUZRd5ma2ArgD+KS7D5r9eu+e\nu7uZ+Wzjtm7denq5q6uLrq6uSrOKANC4tpG7X7gbXggbNzoyyrG7n+Lzl10V/Jz7hoZouuWWJVHm\nAPd97WtcU8G4tmefJZN6Q9XzyNxyuRy5XG5Bj1FWmZtZhqkiv93d74zuzpvZWnfvNrN1QM9sY6eX\nuUg1nLPxHNgYPq7/ZD/5/7uLTRdeGDx27NCh8CesJXc2bQz/Ju3p6SFzoiGGQDKXmRu627ZtC36M\nco5mMeBrwC53v23al+4CtkTLW4A7Z44VEZHFUc6W+VuAPwB+ZWY7o/s+A3wR+J6Z3QQcBG6MJaGI\niMxr3jJ39wc58xb8purGERGRSugToHMoFAqMFgoMTUwEj/VZdwdLNRSLRSYnJ4PHFQqFGNKI1AeV\n+Ryef/55vvHsw9yffzp4bO9kieYYMgns3bWLkwf20xB47pvCeImMCl0SSmU+n7Oc869qCx62f3yU\nQkFHBcRicpJLMhk6W1cEDRsbK1LoKMUUSqS2dFo/EZEE0JZ5HRodGuTwvueCxw2cOE5Pz6yH+8fm\nW9/6AQ8/fCB43MmTJ/n54z/j0mJH8NjixARrGvSuR2Q6lXmdaU6nOa/CohoqliraMbgQp06NMT5+\nLR0dFweNO3lyN8Whu3hDJhv+pJksDSmdX1hkOpV5nTEzGqzCoqpRv2UyzWSzYfPX6XQzZpDVFrZI\nVajM59DT08ORQ87DPhI8dmiogcbGGELNxdPcd9/jHD06Gjx048azue66344hlIgsBpX5HEZHR+k/\n9Wp6V7wpeGw6Del0+BWKFqIhcyU9PV2k0+cGjRse7mF0dA/XXRdTMBGJncp8HqlUlubmpXF233R6\nNatWXUhHx/lB4xoassCeeEKJyKJQmcuCHD78AodeOJtTfWFHuR7v6aZUKrFvX39MyX7TxEQJ613c\nHcS14O7sP/Yi91RwSP0L/YPVDxSzh375EC8cDjwfMnDuunN59atfHUOi2lCZy4IcOPw8u7qbWTEx\nHDRuYmKIyexmnn568aaiisVxmoZ+tmjPV0u/OlDiO8dfHzwuP3Qp2ezSORdFqhXuOXYPLYNhr6Oh\n3iHec+l7VOYi07WsaqG9M+x63hMTLRTGLqO9sSmmVL9pbGyI4jIpc4AL2i+vYFQfx7y76lnikmlJ\ns/4V61nVsSpoXH5/fv6VlhiVecIMDAwEf3BocPAEmczJmBLJYnN3xiYm6BkOe7cEMFwoLLlWOHny\nJGPFsaAxvb299I73xpSoNpbYj03m0giM7t7N8f37g8b1j3YzelH4ycSkfh0f6eeh4cMVjU1nls5Z\nPhrd6X9iJyONYZ9X6HlxkANrws+5VM9U5gmyIpPh5StX0t4UNnVxmH76dc7eRHGctnMWbwqrVtqy\njbxi1SoaG8OqzHrHE3eeapV5wowXi4wEnuZ1rFTUub5lyRodLVIqhRXz+FiRiUz4dQrqmco8QSbT\nzi/7jgSPGxo7Tstk+DiRWisWO3nggUzwuFMDKUYvWtyT0sVNZZ4gq8+u7HIYE70NlCaT9ZZTloc1\na66taNzAyF4mJ5+scpraWjp7OkRE5IyWxZb5X375LznYfTB43IHnDjBp4efbFhFZbMuizE8Nn6Lt\nLW00rQjbu9+7ppfm58Ln40REFtuyKHOAhnQDDemwY1FTDSmsVicJFxEJMG+Zm9nXgd8Betz91dF9\nHcB3gQuAg8CN7t4XY06pU/ljx+g9mWe8GHbFoEl30E5XkaopZwfoN4DNM+67Fdjh7pcA90W3ZRka\nHx6hsTDB6kIh6N+aYpGzsxVcMk5EZjXvlrm7/9zMLpxx9/XA26Ll7UAOFfqylbYUTellM2MnUpcq\n/Q3sdPeXTjuWB+r66g2HDh5ksqVAtiVsS/BEPg9siCdUnRkvjHPHnXcEjxsZC7+knohU34I3p9zd\nzeyMk59bt249vdzV1UVXV9dCnzLYiUMHWds4RmtL2JEpq4+P0p26IKZU9aMhnaLYXOTH3T8OHju6\nepzsKX1cQWQhcrkcuVxuQY9RaZnnzWytu3eb2TrgjJ+LnV7mtXROayvtgYcmpkdSHAg8AmYpamgw\nmtuaWH/Z+uCxrat3MTmoMhdZiJkbutu2bQt+jErL/C5gC/DX0X/vrPBxFsXQ0Aj5fIqRgbCTSQ0O\njseUqP6UCkWOPXsseFxhrEDy/9yJ1L9yDk38NlM7O88ys8PAnwNfBL5nZjcRHZoYZ8iF6u7OMjj0\nWpqyjcFjU6lzYkhUX1KpRkb63sijPzw7eOzQ8SY6Gio7J4yIVE85R7N84Axf2lTlLPHxFK2tV7Ci\nJezSUstFJrOCltaraDvn/OCxY4WjZEa0E1Sk1jTZKSKSACpzEZEEUJmLiCSAylxEJAFU5iIiCaAy\nFxFJAJW5iEgC6FR3IrIs5fNFPve5vw8e19QEH/3oB+joqK9LSqrMRWTZyTaeQzr1QcbHZ16qYX5H\njnybYrEYQ6qFUZmLyLKTSmXIZNppbQ0/XUdvb33WZn2mkkVXmJjgeHd38LiJsbEY0ojEb3R4mF1P\nPBE87uSpfQwODnLOOfV13iaVuZBJpegoFmFgIHhsM5DVVYZkicmkGuiYmGDNkSPBYw/2HWGsDjdi\n9FsoNKRSrND1OGUZSadSrGpspHPFiuCxDQP1edJnlblInfrFL37B008/HTyuVCrFkCZ5RgsFDvf3\nB48bHB1ldHQ0hkQLozIXqVO5B37O328/QOuqsEvsOkDrW2LJlBSZTIp80zD50nDw2Px4PyN1eNpn\nlblIHWtdvZH1l705aIy7U2zYE1OiZGhsTNPYWFn9HTvjRTJra8mU+cDAACdOnKhorPsZrzcty0xp\ncpLcwYPB457q7ubQI48wPBy+JZdOp7nmmmuCx4mEWDJlvnfvXm77zm00tYVdlBmgkC1iFkMoWVJS\nqRR9fjEfuT98nnRoYpTOU19iZeDrb7I0iU+s5MGfPx78nCIhlkyZA2Q7s2y4ckPwuOa7n6QhpdPQ\nLHfZbAsXXPzuisb29T3C6y8fZP36sEsPjk0UyT1wqqLnFAmxZMq8WCxSnChSGCsEj9Usi1RDqeQU\nCmFHihQKJUqlyQqfr4RPmo5OkbIsmTLfu3cvD/1LD089GF7mYwOrYXUMoWQZaWLnzjZ27gwbVSwZ\nJ05Wdh6P5w/s51TPSiZTzwWPTekPwLKzZMocALucts7wt8mDQ/tJa9JcFqC9/XLg8uBxE4Vxenr+\nd8XP22xwflP4fiJZfhY0kWxmm81sj5k9Z2a3VCuUiIiEqXjL3MwagL8DNgFHgUfN7C53312tcHHo\n69tDe/tltY7xHyhT+eox13yZil7i01s/Hfy4u5/fA1wZS6ZaUKZ4LWSa5Spgn7sfBDCz7wC/C6jM\nAylT+eox11yZzKCxM0v2DeHnvmnMZ2nKV/YrutS+T7VSj5kqtZAyXw8cnnb7CPDGhcURSRYzI5Uy\nmlaEz3unGzNYSvt6pDwLKfNFP+Dv+JF/p/f48+EDSyUGSpOkUsbQwB7Gij+sfrgFUKby1WOuuTJN\nThaYTA/xDzf/Q/gDuzHSM8pg//6qZqqVpGQqlcJPFb0YrNKPupvZm4Ct7r45uv0ZYNLd/3raOjrC\nW0SkAu4e9LZsIWWeBp4F3gkcAx4BPlDvO0BFRJKo4mkWdy+a2ceAe4AG4GsqchGR2qh4y1xEROpH\nVc8+ZWYdZrbDzPaa2b1m1n6G9T5jZs+Y2VNm9i0za6xmjgoztZvZ981st5ntivYJ1DRTtG6Dme00\nsx/FlafcTGZ2npndH/3snjazT8SUZd4Po5nZ/4q+/qSZvTaOHKG5zOyDUZ5fmdkvzCz8I6NVzjRt\nvSvNrGhm762HTGbWFb2unzazXK0zmdlZZvZjM3siyvRfY87zdTPLm9lTc6wT9hp396r9A74E/Gm0\nfAvwxVnWuRA4ADRGt78LbKlmjtBM0de2Ax+OltNAW60zRV//NPBN4K648gT87NYCV0TLK5jaZ/KK\nKudoAPZFr5MM8MTM5wDeDfxbtPxG4Jdxfm8Ccr35pdcNsDnuXOVkmrbeT4F/Bd5X60xAO/AMsCG6\nfVYdZNoK/NVLeYCTQDrGTG8FXgs8dYavB7/Gq31e2OuZKkWi/94wyzoDQAFoiXaitjD1CdK4zJvJ\nzNqAt7r712Fqf4C7h5/0uoqZolwbmPqh/iMQ9wHH82Zy9253fyJaHmLqA2LnVjnH6Q+juXsBeOnD\naLNmdfeHgXYzC7u2Wgy53P2haa+bh4Hw8zVXOVPk48D3geMx5yk30+8Dd7j7EQB3r+yqM9XN9CLw\n0vmNVwEn3b2yM6SVwd1/Dsx1buTg13i1y7zT3fPRch74jSd3917gy8ALTB0F0+fuP6lyjqBMwEbg\nuJl9w8weN7OvmllLjTMB/A3wJ0Bl51CNJxMAZnYhU1sWD1c5x2wfRltfxjpxF2c5uaa7Cfi3WBOV\nkcnM1jNVXF+J7op7J1k536eXAx3RlN1jZvaHdZDpq8CrzOwY8CTwyZgzzSf4NR58NIuZ7WDq7fZM\nfzb9hrv7bMeZm9nFwH9n6i1PP/DPZvZBd/9maJZqZWLq+/A64GPu/qiZ3QbcCvx5rTKZ2XuAHnff\naWZdleaoZqZpj7OCqS29T0Zb6NVUbtnMfKcSd0mV/fhm9nbgw0DcV1UuJ9NtwK3Rz9SI/x1eOZky\nTP2+vZOpd+YPmdkv3T38XL/Vy/RZ4Al374o6aoeZvcbdB2PKVI6g13hwmbv7tWd85qkJ/bXu3m1m\n64DZLn36BuD/ufvJaMwPgKuZmheuSBUyHQGOuPuj0e3vM1XmFatCpquB683s3UATsMrM/o+7f6iG\nmTCzDHAH8E/ufmelWeZwFDhv2u3zmPr5zLXOBuKdqis3F9FOz68Cm9097ksMlZPp9cB3pnqcs4B3\nmVnB3e+qYabDwAl3HwVGzexnwGuAuMq8nExXA18AcPf9ZvY8cCnwWEyZ5hP8Gq/2NMtdwJZoeQsw\n2y/7HuBNZtYcbSlsAnZVOUdQJnfvBg6b2SXRXZuY2kFTy0yfdffz3H0j8H7gpwsp8mpkin5eXwN2\nufttMeV4DHi5mV1oZlngv0TZZmb9UJTpTUxN1eWJ17y5zOx84AfAH7j7vpjzlJXJ3S9y943R6+j7\nwM0xFnlZmYAfAtfY1JFaLUzt4IuzA8rJtIep33uiuelLmTpQo1bCX+NV3kPbAfwE2AvcC7RH958L\n3D1tvT9lqiyfYmqSP1PNHBVmeg3wKFPzZT8g3qNZyso0bf23Ef/RLPNmAq5hav7+CWBn9G9zDFne\nxdSRMvuAz0T3fQT4yLR1/i76+pPA6+L83pSbi6kd1SenfW8eqXWmGet+A3hvPWQC/se0DvhErTMx\n9a7lR9Hr6Sng92PO822m9hlOMPVO5cMLfY3rQ0MiIgmgS9aLiCSAylxEJAFU5iIiCaAyFxFJAJW5\niEgCqMxFRBJAZS4ikgAqcxGRBPj/1ebIh1tZkMEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b3607d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "state = initial_state()\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look what happens to the assignments after just one Gibbs step!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG7xJREFUeJzt3X2QXXWd5/H3px/y0HkkJOQZE5EIUVBRwEFdriO6DOMG\ny6lRnBEz4rgzWD6O60iwSjs1JePjDKtTbtWgoWBWI4gMG1xHCbhXcUGIPEMIJIFAOjGd585DJ+mH\n+90/7oFtOp1On9P39O2cfF5VXbn33PO795Pum09O/8655ygiMDOzYmqodwAzM8uPS97MrMBc8mZm\nBeaSNzMrMJe8mVmBueTNzAps0JKXtEJSu6Qn+i3/lKSnJT0p6et9li+TtF7SOknvySu0mZkNTdNx\nHr8R+C5w80sLJL0TWAKcGxHdkmYkyxcDHwQWA3OBuyUtiohKLsnNzOy4Bt2Sj4h7gT39Fl8N/GNE\ndCfr7EiWXw6sjIjuiNgEbAAuqG1cMzNLI8uc/JnAf5L0O0llSW9Jls8B2vqs10Z1i97MzOrkeNM1\nxxpzSkS8VdL5wK3Aq4+xrs+ZYGZWR1lKvg24HSAi1kiqSJoObAHm91lvXrLsFSS5+M3MMogIpR2T\nZbrmDuCPASQtAsZExE5gFXCFpDGSFlKd1nnwGEFH3ddXvvKVumdwJmc6GXM509C+shp0S17SSuBi\n4FRJm4EvAyuAFclhlV3AR5LiXivpVmAt0AN8IoaTzMzMhm3Qko+IDx3joSuPsf51wHXDDWVmZrXh\nT7wmSqVSvSMcxZmGxpmGbjTmcqZ8aaRnVCR5FsfMLCVJxAjteDUzsxOES97MrMBc8mZmBeaSNzMr\nMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVv\nZlZgWS7kbWY2JDt37uTB3/4WKpX0gxsauODtb2f69Om1D3YSccmbWW7279/Pxp/8hAtbWlKPfeDg\nQV57zjku+WE63oW8VwB/CmyPiHP6PfZ54JvA9IjYnSxbBlwF9AKfjoi7ckltZieMyePHc8HcuanH\nrdu8OYc0J5/jzcnfCFzaf6Gk+cC7gRf6LFsMfBBYnIz5niTP+ZuZ1dGgJRwR9wJ7Bnjon4C/77fs\ncmBlRHRHxCZgA3BBLUKamVk2qbe0JV0OtEXE4/0emgO09bnfBqT/Hc3MzGom1Y5XSS3AtVSnal5e\nPMiQyBLKzMxqI+3RNWcAC4DHJAHMAx6SdCGwBZjfZ915ybKjtLa2vny7VCpRKpVSxjAzK7ZyuUy5\nXB728yhi8I1tSQuAO/sfXZM89jzw5ojYnex4/RHVefi5wN3Aa6LfC0jqv8jMCur555/nN1/9Kkvn\nzUs99ubNm3nbtddyxhln5JDsxCOJiBhs5mRAg87JS1oJ3AcskrRZ0kf7rfJyW0fEWuBWYC3wH8An\n3OZmZvU16HRNRHzoOI+/ut/964DrapDLzMxqwMexm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnk\nzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MC\nc8mbmRWYS97MrMBc8mZmBXa8C3mvkNQu6Yk+y74p6WlJj0m6XdKUPo8tk7Re0jpJ78kzuJmZHd/x\ntuRvBC7tt+wu4HUR8QbgWWAZgKTFwAeBxcmY70nybwpmZnU0aAlHxL3Ann7LVkdEJbn7ADAvuX05\nsDIiuiNiE7ABuKC2cc3MLI3hbmlfBfw8uT0HaOvzWBswd5jPb2Zmw9CUdaCkLwFdEfGjQVaLgRa2\ntra+fLtUKlEqlbLGMDMrpHK5TLlcHvbzZCp5SX8FXAa8q8/iLcD8PvfnJcuO0rfkzczsaP03gJcv\nX57peVJP10i6FPgCcHlEHO7z0CrgCkljJC0EzgQezJTKzMxqYtAteUkrgYuB6ZI2A1+hejTNGGC1\nJID7I+ITEbFW0q3AWqAH+EREDDhdY2ZmI2PQko+IDw2weMUg618HXDfcUGZmVhs+jt3MrMBc8mZm\nBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnk\nzcwKLPOVoczs5LHiu99l6+OPpx4XEZze2ZlDIhsql7yZHVdvRwdXTpjAnEmTUo9NrjthdeKSN7Mh\naWpooLmxsd4xLCXPyZuZFZhL3syswDxdY2aj1oYNG+jo6Eg9rqWlhbPOOiuHRCee413IewXwp8D2\niDgnWTYNuAV4FbAJ+EBE7E0eWwZcBfQCn46Iu/KLbmZF9pqGBnbccgubU47r7Opi/9lnc9aXvpRL\nrhPN8bbkbwS+C9zcZ9k1wOqI+IakLyb3r5G0GPggsBiYC9wtaVFEVHLIbWYFd9HcuZnGbd2/n59V\nXDsvGXROPiLuBfb0W7wEuCm5fRPwvuT25cDKiOiOiE3ABuCC2kU1M7O0sux4nRkR7cntdmBmcnsO\n0NZnvTaqW/RmZlYnw9rxGhEhKQZbZaCFra2tL98ulUqUSqXhxDAzK5xyuUy5XB7282Qp+XZJsyJi\nm6TZwPZk+RZgfp/15iXLjtK35M3M7Gj9N4CXL1+e6XmyTNesApYmt5cCd/RZfoWkMZIWAmcCD2ZK\nZWZmNXG8QyhXAhcD0yVtBr4MfA24VdLHSA6hBIiItZJuBdYCPcAnImKwqRwzM8vZoCUfER86xkOX\nHGP964DrhhvKzMxqw6c1MDMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXm\nkjczKzCXvJlZgbnkzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3MyuwLBfytow6Ojr4\n1rdW0NmZfmxjI1x55SWcc845tQ9mZoXlkh9BlUqFtrZg5syrUo/duvVuurq6ckhlZkWWueQlLQM+\nDFSAJ4CPAhOAW4BXkVzkOyL2Dj9mcTQ0NDJu3NTU45qaxuaQxsyKLtOcvKQFwMeB8yLiHKARuAK4\nBlgdEYuAe5L7ZmZWJ1l3vO4DuoEWSU1AC7AVWALclKxzE/C+YSc0M7PMMpV8ROwGvg28SLXc90bE\namBmRLQnq7UDM2uS0szMMsk0Jy/pDOCzwAKgA/iJpA/3XSciQlIMNL61tfXl26VSiVKplCWGmVlh\nlctlyuXysJ8n647XtwD3RcQuAEm3A38EbJM0KyK2SZoNbB9ocN+SNzOzo/XfAF6+fHmm58k6J78O\neKuk8ZIEXAKsBe4ElibrLAXuyPj8ZmZWA5m25CPiMUk3A7+negjlw8C/ApOAWyV9jOQQyhrlNDOz\nDDIfJx8R3wC+0W/xbqpb9WZmNgr43DVmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZgLnkzswJzyZuZ\nFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MCc8mbmRWYS97MrMBc8mZmBeaS\nNzMrMJe8mVmBZS55SVMl3SbpaUlrJV0oaZqk1ZKelXSXpKm1DGtmZukMZ0v+vwM/j4izgXOBdcA1\nwOqIWATck9w3M7M6yVTykqYA74iIFQAR0RMRHcAS4KZktZuA99UkpZmZZZJ1S34hsEPSjZIelnSD\npAnAzIhoT9ZpB2bWJKWZmWXSNIxx5wGfjIg1kq6n39RMRISkGGhwa2vry7dLpRKlUiljDDOzYiqX\ny5TL5WE/T9aSbwPaImJNcv82YBmwTdKsiNgmaTawfaDBfUvezMyO1n8DePny5ZmeJ9N0TURsAzZL\nWpQsugR4CrgTWJosWwrckSmVmZnVRNYteYBPAT+UNAbYCHwUaARulfQxYBPwgWEntJft2bOHtra2\n1OMaGhqYM2dODonMbLTLXPIR8Rhw/gAPXZI9jh2LNJWVK9cBz6caF1Fh2rS9fOtbX8wnmJmNasPZ\nkrcRNH/+O4B3pB7X3d3J3r3frX0gMzsh+LQGZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MC\nc8mbmRWYS97MrMD8YagMuru7qVQqqccdOXIkhzRmZsfmks/gxz++g7vvfoaGhsbUYzs7fUVEMxs5\nLvkMurpgwoT3cdppr693FDOzQXlO3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MC\nG1bJS2qU9IikO5P70yStlvSspLsk+ZM/ZmZ1NNwt+c8Aa4FI7l8DrI6IRcA9yX0zM6uTzJ94lTQP\nuAz4KvB3yeIlwMXJ7ZuAMi56s1Fh/fr1rLrhBujuTj324M6dMGNGDqnyse2ZZ/j25z6XfmBzM//l\nr/+aRYsW1T5UnQzntAb/DHwBmNxn2cyIaE9utwMzh/H8ZlZDPT09zNy+nSWzZ6cfPG8eE5qbax8q\nBzMnTOCz8+dnGnvnH/5AT09PjRPVV6aSl/ReYHtEPCKpNNA6ERGSYqDHWltbX75dKpUolQZ8CjOr\nsabGRiaPHVvvGLlqbGjI/Hdsbkx/0sG8lMtlyuXysJ8n65b8RcASSZcB44DJkv4NaJc0KyK2SZoN\nbB9ocN+SNzOzo/XfAF6+fHmm58m04zUiro2I+RGxELgC+FVEXAmsApYmqy0F7siUyszMaqJWx8m/\nNC3zNeDdkp4F/ji5b2ZmdTLs88lHxK+BXye3dwOXDPc5zcysNvyJVzOzAnPJm5kVmC//d1II9uzZ\nk2nk+PHjGTduXI3z5GPv3r1EDHjU7nGNHTuWlpaWGicyqz+XfOGJnTvH8fnP35R6ZE/PIf72by/m\noosuyiFX7S3/5nI6ejtA6cZ1d3Xz5+/8c5a8d0k+wczqyCVfcM3N4zn77M9mGvvCC3fVOE2+DnYd\nZO5/nktjU7oPtLStbaMSlZxSmdWX5+TNzArMJW9mVmAueTOzAvOcvI0qGzdu5Olnns40tqe3WGcP\nNKsFl7yNKm1tbfzg3h8wbe601GP1aiGlPLTGrOBc8jbqTJ4xmfmvz3Y+cDN7Jc/Jm5kVmLfkreYi\ngscffzzT2A0bNnBo36EaJzI7ebnkreYigu9859+JODf12Oeea6P5rI4cUpmdnE7akt+2bRvXXfd9\nDmXYaDx8uJdJk86ufagC6e1tYMGC96cet3VrFxse/g272h5LPXZsyxEu+a9vTv2JV7MiO2lLHuDg\nwVOZO/fjmcZK3p2Rh5kzL2BrbxOTpp2SemzH9u/lkMjsxHZSlzyIhoaT/FswykgNqKERNWTYGvfh\nk2ZH8eaomVmBZdqMlTQfuBk4jer1Xf81Ir4jaRpwC/AqYBPwgYjYW6OsZrm59/f38uSGJ1OPG988\nnquvutrnordRK+tcRTfwuYh4VNJE4CFJq4GPAqsj4huSvghck3yZjVqnLTyNwzMOs5vdqcfufmA3\nvb29OaQyq41MJR8R24Btye0Dkp4G5gJLgIuT1W4CyrjkbZQbM34MY8aPyTR2b4N/UbXRbdhz8pIW\nAG8CHgBmRkR78lA7MHO4z29mZtkNq+STqZqfAp+JiP19H4vqxTazXXDTzMxqIvPxg5KaqRb8v0XE\nHcnidkmzImKbpNnA9oHGtra2vny7VCpRKpWyxjAzK6RyuUy5XB7282Q9ukbAD4C1EXF9n4dWAUuB\nryd/3jHA8FeUvJmZHa3/BvDy5cszPU/WLfm3AR8GHpf0SLJsGfA14FZJHyM5hDLj85uZWQ1kPbrm\ntxx7Pv+S7HFstNm8eTNr1qxJNaZSqeSUxszS8mf67ZhaWubxy192c4xdK8cUEfT2nueyNxsFXPJ2\nTDNmLAYWpx4XEfzi7l/wxPpV2V54SrZhZnY0l7zl4kjvEaaeNXVEX7PS08CTv3o+03VeT1swkVln\n+mMdJ71KhQfLZZ595JHjr9vPhFNO4d2XXZZDqOFxyVthjBn/HjanP/0Mh/ZvYcz4jcw6s/aZ7MRy\n/uTJdGS4qtmBri4enjHDJW+Wp4nTXpNpXKW3C9hY2zB2Qlp4SvrrGADs6uzk4RpnqRWfatjMrMBc\n8mZmBebpmhEUEcM6rFASDQ3+f9nMhs4lP4IOHz7MPffeQ0/0pB4rxPnnnM+cOXNySGZmReWSH2G9\njb1MPSP9oYUdWzpySGMnol/+7Gcc2Lkz9biOffto8QfUTjouebMTzPr77uO89nYmjUl/oZPJU/xJ\ns5ONS97sBLTo1FOZ7uvK2hC45E8Qlajw+LrHeWrDUyP2mkK87szXMXv27BF7zRPNwZ6D/MP1/4Aa\n0n/K9uI3X8xll46+D89YsbjkM9i5cyednZ2px3V3d2d+zckzJ4/4Cb/2bd/Hjh07UueuXhTs5LDg\nnQsyXf9sxws72LNvT63jmB3FJZ/Bi1tfZOOujTSPbU49VhPTb/EBNDY30khjprFZjZs8jucOPAcH\n0o+NySdH0Y+fND7TuKwXDh+O3kqFrt7ezOPHNjXRkOG8QFZfLvksAlqmtjDhlAn1TpKrlikttEzx\nvG9RPLtrF99aez/RmP4/4OaeBr583sXMnTw5h2SWpxO+5O+559ds2LAt9bg9e3bz9DMbaN/1+9Rj\n93bsBb/X7UR0anD669MfYfP8Q3u57bm1jGtKXxkLJ5zCexaekXqc1cYJX/JPPbWZhx6aw8SJs1KN\n6+zsZPch6KY9/YtOgbEtY9OPs1wdOnQo07gjR45w+PDhGqcplsYZsKZrS+pxnZ09vOqFKTRGtmme\n0yZMYHxz+mnRsY2NzPfhokABSh5gypTTmZbyDIQHDhxg/MQtTJha7CmXk8n2tjbG9vamPp/8gY6t\n7M7w4aKTyYsvNrJjx2tobExXuN3dvTyxv4efbzqY+jUrlXbedNYeJk9Mt0HV3dPL2TGDa9789tSv\nWUQ1L3lJlwLXA43A9yPi67V+DbMBVSqcNnZs6p2DlYaR3aF9IqpUmmhpeRtjx2Y4Fe+p2V6zo+M/\nmHVWN9Ompdu5vX//Ebqfzb6DuWhqWvKSGoF/oXox7y3AGkmrIuLpWr5OHvZue5Sps95Y7xiv4ExD\nU4tMhzo62fTopkxj5y2eR9OYV/5T2vToJha8ccGwMtXankOHeHL9dnonpt/x2tWV34nx9u5dx9Sp\nZ+X2/FmUN22itGBBvWPURK235C8ANkTEJgBJPwYuB1zyGfTNVKlU6nL8ef8zX+b9fcry99y99WEm\nzTgn82s2j5nCjufezUM70+9Nj3iAWWf2pC753t5eDh0+xIED6Y9PrWR8H+w+dIhnn5/AYd6Tafyk\nSeMyjTuewUq+p6dCV1e6rfKurl6O9PRwoKsrU55xTU0u+UHMBTb3ud8GXFjj1zgp7d61i4O7dtLA\nyB2nXCGYOnsOU0ZwB9a+ffvo+MMfUk257N+9m63r19PQ2wsZzucydtw0DlVOoXN/6qEcPLidrq7T\nGTchXQHu2bOHX/zkZvaW70v/op2d6PTT048DGhrGcMop6S/OXg+9vZO5//5Jqcd1dzcwpqORrdte\nSD1WdPL5i4p1plfVcutQ0p8Bl0bEx5P7HwYujIhP9VknavmaV199DbffvhEp/T/uIxyhYVy1TDr3\nPkXL1NfVLFctvCJTAEeO0Nw0cvPH3T29MO6VO71G5Pt0ON3f88C+p5k4+ewcAx1bV2UDY6bsOGr5\ngR0HmDhj4jHHdR/qZsJBmJByR+ZwbepoYMqpo2u7q9Y/v97uCmN7mmjKcE2k7t4tTB73HDsOHGDG\nxGP//PqLCM4+91z+ffXq1K85VJKISH+YUq1L/q1Aa0RcmtxfBlT67nyVdHJ8FNLMrMZGQ8k3Ac8A\n7wK2Ag8CHzoRdryamRVRTefkI6JH0ieBX1I9hPIHLngzs/qp6Za8mZmNLrlfFVrSNEmrJT0r6S5J\nA177TtIySU9JekLSjyTldt6AFJmmSrpN0tOS1ib7HHIz1FzJuo2SHpF0Z70zSZov6f8kP78nJX06\npyyXSlonab2kLx5jne8kjz8m6U155EiTSdJfJlkel/R/JZ1b70x91jtfUo+k94+GTJJKyXv6SUnl\nvDMNJZek6ZJ+IenRJNdf5ZxnhaR2SU8Msk6693hE5PoFfAP4++T2F4GvDbDOAuA5YGxy/xZgaT0z\nJY/dBFyV3G4CptT7e9Vn3b8DfgisqncmYBbwxuT2RKr7Zc6ucY5GYEPyXmkGHu3/GsBlwM+T2xcC\nv8v5ezOUTH/00vsGuHQ0ZOqz3q+AnwF/Vu9MwFTgKWBecn96nplS5GoF/vGlTMAuoCnHTO8A3gQ8\ncYzHU7/Hc9+SB5ZQLUuSP983wDr7gG6gJdl520L1E7N1yyRpCvCOiFgB1f0NEZH31bSH8r1C0jyq\nP+zvQ+4Hzh83U0Rsi4hHk9sHqH74rdYHG7/8QbuI6AZe+qDdgFkj4gFgqqSZNc6RKlNE3N/nffMA\nMC/HPEPKlPgUcBtw9PGf9cn0F8BPI6INICJG4mRCQ8n1B/7/OWcnA7sioievQBFxLzDY1WRSv8dH\nouRnRsRLp3psB44KFBG7gW8DL1I9KmdvRNxdz0zAQmCHpBslPSzpBkl5n1x9KLkA/hn4AjASl4oa\naiYAJC2guiXyQI1zDPRBu7lDWCfPUh1Kpr4+Bvw8xzwwhEyS5lIts/+RLMp7x9xQvk9nAtOSab/f\nS7oy50xDzXUD8DpJW4HHgM+MQK7BpH6P1+ToGkmrqf7K3t+X+t6JiBjoOHlJZwCfpfprUwfwE0l/\nGRE/rFcmqt+b84BPRsQaSdcD1wBfzpqpFrkkvRfYHhGPSCoNJ0utMvV5nolUtw4/k2zR19JQi6j/\nbzZ5FtiQn1vSO4GrgLflFwcYWqbrgWuSn6fI/7fBoWRqpvrv7V1Uf5O/X9LvImJ9nXNdCzwaEaWk\np1ZLekNEZPh8dM2keo/XpOQj4t3HTFPdiTArIrZJmg1sH2C1twD3RcSuZMztwEVU55zrlakNaIuI\nNcn926iW/LDUINdFwBJJlwHjgMmSbo6Ij9QxE5KagZ8C/zMi7siaZRBbgPl97s+n+jMabJ155Dvt\nN5RMJDtbb6D6afC8L+w6lExvBn5c7XemA38iqTsiVtUx02ZgZ0QcAg5J+g3wBiDPkh9KrouArwJE\nxEZJzwOvBdJfbag2Ur/HR2K6ZhWwNLm9FBioANYBb5U0PtmyuARYW89MEbEN2CxpUbLoEqo7hvI0\nlFzXRsT8iFgIXAH8ajgFX4tMyc/sB8DaiLg+pxy/B86UtEDVc1h8MMnWP+tHkkxvpTrtl+GqMLXL\nJOl04HbgwxGxIccsQ84UEa+OiIXJe+g24OocC35ImYD/Bbxd1aPGWqjuVMyzA4aaax3Vf/skc9+v\npXqQSL2kf4/ntZe4z97gacDdwLPAXcDUZPkc4H/3We/vqZboE1R3LDSPgkxvANZQnYu7nfyPrhlS\nrj7rX0z+R9ccNxPwdqr7Bx4FHkm+Ls0hy59QPXJnA7AsWfY3wN/0WedfkscfA84bgff3oJmo7hzf\n1ef78mC9M/Vb90bg/aMhE/Df+nTAp/PONMSf33TgzuT99ATwFznnWUl1v2QX1d9urhrue9wfhjIz\nK7CRmK4xM7M6ccmbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmD/D6/iAGl7dGkt\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b839f90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gibbs_step(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGP9JREFUeJzt3X2QHPV95/H3R08gsXoWJ6EHWzxIDnIBAWxMsH2MzUMp\nSg6oJGUgxlZiKkmZGNvni21JqbJXzplYXPnCXVLO1TkWJ3KALTCnA+zEWuMMxodBknk0CycJI9AC\nWskgCa0kpH343h/T6Car1e5MT7dmt/V5VW2pp6d/PR92h0/1/LpnRhGBmZkV06hmBzAzs/y45M3M\nCswlb2ZWYC55M7MCc8mbmRWYS97MrMAGLXlJqyV1Snq23/qbJT0v6ZeSVlWtXy5pi6QXJF2ZV2gz\nM6vNmCHuvx34W+COd1ZI+ghwFXBuRHRLOjVZvwi4FlgEzAF+LGlhRPTlktzMzIY06JF8RDwC7O63\n+tPAX0dEd7LNrmT91cDdEdEdEduArcBF2cY1M7N6pJmTXwD8W0mPSSpLel+yfjbQUbVdB5UjejMz\na5KhpmuONWZqRFws6f3AWuCMY2zrz0wwM2uiNCXfAdwHEBEbJfVJmgG8Csyr2m5usu5fkeTiNzNL\nISJU75g00zXrgI8CSFoIjIuIXwP3A9dJGifpdCrTOhuOEXTY/Xz1q19tegZncqYTMZcz1faT1qBH\n8pLuBi4FpkvaDnwFWA2sTi6rPAx8MinudklrgXagB7gpGklmZmYNG7TkI+L6Y9z1iWNsfwtwS6Oh\nzMwsG37Ha6JUKjU7wlGcqTbOVLvhmMuZ8qXjPaMiybM4ZmZ1kkQcpxOvZmY2QrjkzcwKzCVvZlZg\nLnkzswJL845XM7OavP7666y/5x7o7q5/8NixXPEHf8Ds2bOzD3YCccmbWW7efvttDm7axJVTp9Y9\ntm33bg7+zu/kkOrE4pI3s1yNHzeOM1KU/PiurhzSnHg8J29mVmAueTOzAnPJm5kVmEvezKzAXPJm\nZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZgLnkzswIbtOQlrZbUmXxpd//7/oOkPknTqtYtl7RF0guS\nrswjsJmZ1W6oI/nbgcX9V0qaB1wBvFy1bhFwLbAoGfMtSX6lYGbWRIOWcEQ8Auwe4K7/DHyp37qr\ngbsjojsitgFbgYuyCGlmZunUfaQt6WqgIyKe6XfXbKCj6nYHMKeBbGZm1qC6PmpY0gRgBZWpmiOr\nBxkSaUKZmVk26v08+TOB+cDTkgDmAr+Q9AHgVWBe1bZzk3VHaW1tPbJcKpUolUp1xjAzK7ZyuUy5\nXG54P4oY/GBb0nzggYg4Z4D7XgIujIg3kxOvd1GZh58D/Bg4K/o9gKT+q8ysoF566SV++vWvs3Tu\n3LrH3rF9Ox9csYIzzzwzh2QjjyQiYrCZkwENdQnl3cCjwEJJ2yX9cb9NjrR1RLQDa4F24J+Am9zm\nZmbNNeh0TURcP8T9Z/S7fQtwSwa5zMwsA76O3cyswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZW\nYC55M7MCc8mbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEve\nzKzAXPJmZgXmkjczK7Chvsh7taROSc9WrftPkp6X9LSk+yRNrrpvuaQtkl6QdGWewc3MbGhDHcnf\nDizut2498N6IOA/YDCwHkLQIuBZYlIz5liS/UjAza6JBSzgiHgF291vXFhF9yc3HgbnJ8tXA3RHR\nHRHbgK3ARdnGNTOzejR6pP0p4IfJ8mygo+q+DmBOg/s3M7MGjEk7UNJfAocj4q5BNouBVra2th5Z\nLpVKlEqltDHMzAqpXC5TLpcb3k+qkpf0R8AS4LKq1a8C86puz03WHaW65M3M7Gj9D4BXrlyZaj91\nT9dIWgx8Ebg6It6uuut+4DpJ4ySdDiwANqRKZWZmmRj0SF7S3cClwAxJ24GvUrmaZhzQJgng5xFx\nU0S0S1oLtAM9wE0RMeB0jZmZHR+DlnxEXD/A6tWDbH8LcEujoczMLBupT7ya2Ylj44YNvLVnT93j\n9uzdC35B31QueTMb0lPr13NaezuTTjqprnGnAgtOPjmfUFYTl7yZ1eT8WbOYM2lSs2NYnfyxA2Zm\nBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzA/GYoMxu2fvWrX7F///66x40f\nP54FCxbkkGjkccmb2bB0xqhR7Pzud9lX57gDhw9z4OyzWbBiRS65RhqXvJkNSx+ak+7bQ1/bt48H\ne3szTjNyeU7ezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKbNCSl7RaUqekZ6vWTZPUJmmzpPWSplTd\nt1zSFkkvSLoyz+BmZja0oY7kbwcW91u3DGiLiIXAQ8ltJC0CrgUWJWO+JcmvFMzMmmjQEo6IR4Dd\n/VZfBaxJltcA1yTLVwN3R0R3RGwDtgIXZRfVzMzqleZIe2ZEdCbLncDMZHk20FG1XQeQ7t0MZmaW\niYbe8RoRISkG22Sgla2trUeWS6USpVKpkRhmZoVTLpcpl8sN7ydNyXdKmhUROySdBuxM1r8KzKva\nbm6y7ijVJW9mZkfrfwC8cuXKVPtJM11zP7A0WV4KrKtaf52kcZJOBxYAG1KlMjOzTAx6JC/pbuBS\nYIak7cBXgG8AayXdCGwDPgYQEe2S1gLtQA9wU0QMNpVjZmY5G7TkI+L6Y9x1+TG2vwW4pdFQZmaW\nDV/HbmZWYP48+eOop6eH119/PfX4qVOn0tLSkmEiMys6l/xxtG/fPr72tf9BxOy6xx469Caf//xH\nufDCC3NIZmZF5ZI/zvr6JjFv3o11j3vllQdySGNmRec5eTOzAnPJm5kVmEvezKzAXPJmZgXmkjcz\nKzCXvJlZgbnkzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCix1\nyUtaLuk5Sc9KukvSSZKmSWqTtFnSeklTsgxrZmb1SVXykuYDfwJcEBHnAKOB64BlQFtELAQeSm6b\nmVmTpD2SfwvoBiZIGgNMAF4DrgLWJNusAa5pOKGZmaWWquQj4k3gm8ArVMp9T0S0ATMjojPZrBOY\nmUlKMzNLJdV3vEo6E/g8MB/YC9wj6YbqbSIiJMVA41tbW48sl0olSqVSmhhmZoVVLpcpl8sN7yft\nF3m/D3g0It4AkHQf8FvADkmzImKHpNOAnQMNri55MzM7Wv8D4JUrV6baT9o5+ReAiyWNlyTgcqAd\neABYmmyzFFiXcv9mZpaBVEfyEfG0pDuATUAf8ATw34GJwFpJNwLbgI9llNPMzFJIO11DRNwK3Npv\n9ZtUjurNzGwY8DtezcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3M\nCswlb2ZWYC55M7MCc8mbmRWYS97MrMBc8mZmBZb6o4bt+OrthXXrfsb69U/VPXb69JP40z+9YegN\nzaxwXPIjxKxZH6Sr6zy6uuob19NziI6O+/IJZWbDnkt+hBg/fhrjx0+re1x39wH27MkhkJmNCJ6T\nNzMrMJe8mVmBpS55SVMk3SvpeUntkj4gaZqkNkmbJa2XNCXLsGZmVp9GjuT/C/DDiDgbOBd4AVgG\ntEXEQuCh5LaZmTVJqpKXNBn4cESsBoiInojYC1wFrEk2WwNck0lKMzNLJe2R/OnALkm3S3pC0rcl\nnQLMjIjOZJtOYGYmKc3MLJW0l1COAS4APhMRGyXdRr+pmYgISTHQ4NbW1iPLpVKJUqmUMoaZWTGV\ny2XK5XLD+0lb8h1AR0RsTG7fCywHdkiaFRE7JJ0G7BxocHXJm5nZ0fofAK9cuTLVflJN10TEDmC7\npIXJqsuB54AHgKXJuqXAulSpzMwsE4284/Vm4E5J44AXgT8GRgNrJd0IbAM+1nBCMzNLLXXJR8TT\nwPsHuOvy9HHMzCxL/uyaFB58cD0bN75Y97ienh4OHswhUIFs2rSJB37yQKqxE8ZN4At//gXGjh2b\ncSqzkcsln0Jn5162bz+XadPOqnvs3Lmjc0hUHPv372dLzxZm/8bsuse+/LOXc0hkNrK55FM6+eTJ\ntLT4bQB5GHvyWFqmtdQ97g29kUMas5HNH1BmZlZgPpK3wuju6+bOe+5EUt1jzz37XM7/zfNzSGXW\nXC55K4yJ50zksQOP1T1u92u7mdwy2SVvheSSt8KYeUa6cyQ9h3syTmI2fHhO3syswFzyZmYF5uka\nsxNEb28vPT3ppqb6YsAPlLURwCVvdoLYvHkz93zzm4zt7a1/8MGDaM6c7ENZ7lzyZieQhRFc9+53\nNzuGHUeekzczKzCXvJlZgbnkzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswBq6Tl7SaGAT0BER/07S\nNOB7wLtJvsg7IvY0nNIa0tfXy5YtW1KNnTFjBlOnTs04kVm+3j50KPVzftasWUycODHjRM3T6Juh\nPge0A+/8RpYBbRFxq6QvJ7eXNfgY1gBpNPv2vZtVqzbUPXb//l3cfPNFXHLJJTkkM8vHSaNHM/2V\nV9iwalXdY18/fJgly5axaNGiHJI1R+qSlzQXWAJ8HfhCsvoq4NJkeQ1QxiXfVGPGnMRZZ3081diX\nX16fcRqz/E2fMIGPn3FGqrFrt2/POE3zNTIn/zfAF4G+qnUzI6IzWe4E/CWoZmZNlOpIXtLvAjsj\n4klJpYG2iYiQNOBH17W2th5ZLpVKlEoD7sJGqIhg3759qcYeOnQo4zRmI1O5XKZcLje8n7TTNZcA\nV0laApwMTJL0j0CnpFkRsUPSacDOgQZXl7wVT0Sw4j+u4IAOpBo/eu7ojBOZjTz9D4BXrlyZaj+p\nSj4iVgArACRdCvxFRHxC0q3AUmBV8u+6VKlsxOvq7mL+1fObHcPshJfVdfLvTMt8A7hC0mbgo8lt\nMzNrkoY/Tz4iHgYeTpbfBC5vdJ9mZpYNv+PVzKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZg\nLnkzswJr+Dr5kaqzs5NVq1Zz8GD9Yw8c6OaUU87OPpSZWcZO2JKPCPbsmczs2TfWPXb6dBg16oT9\n1ZnZCHJCN5U0ijFjTmp2DDOz3HhO3syswFzyZmYF5pI3Myswl7yZWYGd0Cdezd6x7qF1/ODhH9Q9\nrmVcC19b/jUmTpyYQyqzxrnk7YQ35zfmEO8Z8OuIh9Txo46M05hlyyVvJzyNEkLNjmGWC8/Jm5kV\nWKqSlzRP0r9Iek7SLyV9Nlk/TVKbpM2S1kuakm1cMzOrR9oj+W7g30fEe4GLgT+XdDawDGiLiIXA\nQ8ltMzNrklQlHxE7IuKpZLkLeB6YA1wFrEk2WwNck0VIMzNLp+E5eUnzgfOBx4GZEdGZ3NUJzGx0\n/2Zmll5DJS+pBfg+8LmI2Fd9X0QEkO66NDMzy0TqSygljaVS8P8YEeuS1Z2SZkXEDkmnATsHGtva\n2npkuVQqUSqV0sYwMyukcrlMuVxueD+pSl6SgO8A7RFxW9Vd9wNLgVXJv+sGGP6vSt7MzI7W/wB4\n5cqVqfaT9kj+g8ANwDOSnkzWLQe+AayVdCOwDfhYyv2bmVkGUpV8RPyMY8/nX54+jpmZZcnveDUz\nKzCXvJlZgfkDyszM3tHXx4Zymc1PPjn0tv2cMnUqVyxZkkOoxrjk7ZgOHTrEjh1dvPLKK3WN6+vr\nyymRWb7eP2kSe595pu5xXYcP88Spp7rkbWTZtWsn/+3e7zHvmXl1jz10yqEcEpnl6/SpU1ONe+PA\nAZ7IOEtWXPI2qHGzxjGvVH/Jm9nw4BOvZmYF5pI3Myswl7yZWYF5Tt6OKYCdL+3nyX/aWvfYMeP6\nOOeyhdmHMn704IN0/frXdY/b+9ZbTPCVTyccl7wd08SJ76Fj1zX09U2sb2D00dv9A865LJ9cJ7ot\njz7KBZ2dTBw3ru6xkyZPziGRDWcueTum8eP/DeMnzqdlan3FENHHWwN+yLRlZeH06cyYMKHZMWwE\ncMmbNWB/z37+6ra/QqNU99hLL7yUJYuH35tnrFhc8mYNmP+R+am+/2zXy7vY/dburOOYHcUlb9aA\n8RPHpxo3bnz98+lmafgSSjOzAhvxR/IPPfQwW7fuqHvcoUNv09OTQyAzs2FkxJf8c89t5xe/mE1L\ny6y6x06ZcnIOiawiePPVN1ONPLnlZCZM9pUjZlkY8SUPMHnyu5g27axmx7AqPd0L+dld9b9UOvz2\nPi5YspczLjw9h1RmJ57MS17SYuA2YDTwDxGxKuvHsOFNGsX0uekuDdyz4ymg/i9sMLOBZXriVdJo\n4O+AxcAi4HpJZ2f5GHnZtq3c7AhHGY6ZKiU8vGx7aluzIxxlOGYCKG/b1uwIR3GmfGV9dc1FwNaI\n2BYR3cB3gaszfoxcDMdCzSJTRHD48OFUP729vUft73iUfG93L4cPHq7558WNLx5Zjkhx0XoOhir5\n3t5eDr59kK6urrp/+hr4bxyO5eVM+cp6umYOsL3qdgfwgYwfw+rQ09NDW7mNHtU/Px4Eo6Ye36ts\nNWos7Q8H7T99ueYxr2/ZQ9fulxkzbjeLP3Mho8eMzjFhNnbv3s0/33MHe8qP1j/4wAH0rndlH8oK\nSVke+Uj6fWBxRPxJcvsG4AMRcXPVNpHlY37608u4774XkRp7c0lX1zO0tJybUapsZJWpm8OQ0YVE\nB/Y8x4Qp781mZxn5/5k2M3Z8Z7PjANC1q4uWU1uOeX/3wW5O2Q+njB57HFPBrq4uTm05dq5mKEKm\niODsc8/lf7W15ZZJEhFR9+dnZF3yFwOtEbE4ub0c6Ks++SppeLyeNjMbYYZDyY8B/i9wGfAasAG4\nPiKez+xBzMysZpnOyUdEj6TPAD+icgnld1zwZmbNk+mRvJmZDS+5XzohaZqkNkmbJa2XNOUY2y2X\n9JykZyXdJemkYZBpiqR7JT0vqT0555CbWnMl246W9KSkB5qdSdI8Sf+S/P1+KemzOWVZLOkFSVsk\nffkY2/zX5P6nJZ2fR456Mkn6eJLlGUn/R1LuZ/dr+T0l271fUo+k3xsOmSSVkuf0LyWV885USy5J\nMyT9s6Snklx/lHOe1ZI6JT07yDb1PccjItcf4FbgS8nyl4FvDLDNfOBXwEnJ7e8BS5uZKblvDfCp\nZHkMMLnZv6uqbb8A3Anc3+xMwCzgN5PlFirnZc7OOMdoYGvyXBkLPNX/MYAlwA+T5Q8Aj+X8u6kl\n02+987yh8ibBpmeq2u4nwIPA7zc7EzAFeA6Ym9yekWemOnK1An/9TibgDWBMjpk+DJwPPHuM++t+\njh+Pi6CvolKWJP9eM8A2bwHdwITk5O0E4NVmZpI0GfhwRKyGyvmGiNibY6aaciXZ5lL5Y/8DUP9X\nEmWcKSJ2RMRTyXIX8DwwO+MctbzR7kjWiHgcmCJpZsY56soUET+vet48DszNMU9NmRI3A/cCu3LO\nU2umPwS+HxEdABFR/zeV55PrdWBSsjwJeCMicvv82oh4BBjs22Tqfo4fj5KfGRHvXLzcCRwVKCLe\nBL4JvELlqpw9EfHjZmYCTgd2Sbpd0hOSvi0p749GrCUXwN8AXwT6cs5TTyYAJM2nciTyeMY5Bnqj\n3ZwatsmzVGvJVO1G4Ic55oEaMkmaQ6XM/j5ZlfeJuVp+TwuAacm03yZJn8g5U625vg28V9JrwNPA\n545DrsHU/RzP5OoaSW1UXrL395fVNyIiBrpOXtKZwOepvGzaC9wj6eMRcWezMlH53VwAfCYiNkq6\nDVgGfCVtpixySfpdYGdEPCmp1EiWrDJV7aeFytHh55Ij+izVWkT9X9nkWWA171vSR4BPAR/MLw5Q\nW6bbgGXJ31Pk/2qwlkxjqfz/dhmVV/I/l/RYRGxpcq4VwFMRUUp6qk3SeRGxL8dcQ6nrOZ5JyUfE\nFcdMUzmJMCsidkg6Ddg5wGbvAx6NiDeSMfcBl1CZc25Wpg6gIyI2JrfvpVLyDckg1yXAVZKWUHkf\n6yRJd0TEJ5uYCUljge8D/zMi1qXNMohXgXlVt+dR+RsNts1c8p32qyUTycnWb1N5N3jeX+xaS6YL\nge9W+p0ZwG9L6o6I+5uYaTvw64g4CByU9FPgPCDPkq8l1yXA1wEi4kVJLwHvATblmGswdT/Hj8d0\nzf3A0mR5KTBQAbwAXCxpfHJkcTnQ3sxMEbED2C5pYbLqcionhvJUS64VETEvIk4HrgN+0kjBZ5Ep\n+Zt9B2iPiNtyyrEJWCBpviqfYXFtkq1/1k8mmS6mMu2X5+ccDJlJ0ruA+4AbImJrjllqzhQRZ0TE\n6clz6F7g0zkWfE2ZgP8NfEiVq8YmUDmpmGcH1JrrBSr/75PMfb+HykUizVL/czyvs8RVZ4OnAT8G\nNgPrgSnJ+tnAD6q2+xKVEn2WyomFscMg03nARipzcfeR/9U1NeWq2v5S8r+6ZshMwIeonB9458Pg\nn6Ry1Jp1lt+mcuXOVmB5su7PgD+r2ubvkvufBi44Ds/vQTNROTn+RtXvZUOzM/Xb9nbg94ZDJuAv\nqjrgs3lnqvHvNwN4IHk+PQv8Yc557qZyXvIwlVc3n2r0Oe43Q5mZFdjx/RxZMzM7rlzyZmYF5pI3\nMyswl7yZWYG55M3MCswlb2ZWYC55M7MCc8mbmRXY/wOc2FQxfK3SswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b82e0d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gibbs_step(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After just two steps, our assignments look _really_ good. We can run it a few more times and see the assignments again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGQJJREFUeJzt3XuQHOV97vHvo9Wd1RUoCV3wCiz5IBdgwAZi4zC2gVJI\nIqgkZSDGVmIqlTLx7eT4IpEqe+WUseGUT/Ah5VTZsbDIwTgCEx3ATiwZZxwczM1IgBFYEmZBK9AK\nC10Quu1qf/ljGmW8Wu3O9HRrdnufT5WKnp5+e5/aHZ7qebt7RhGBmZkV06hmBzAzs/y45M3MCswl\nb2ZWYC55M7MCc8mbmRWYS97MrMAGLHlJKyR1SXq6z/pPSHpW0i8l3VS1fpmkTZKek3RZXqHNzKw2\nowd5/jbgVuD2N1dIeh+wGDgrIrolnZysXwhcBSwEZgM/lrQgInpzSW5mZoMa8Eg+Ih4EdvZZ/THg\nKxHRnWzzarL+CuDOiOiOiA5gM3B+tnHNzKweaebk5wO/K+lhSWVJ70zWzwI6q7brpHJEb2ZmTTLY\ndM2xxkyLiAslvQtYBZx2jG39mQlmZk2UpuQ7gXsAIuIxSb2STgK2AnOrtpuTrPstklz8ZmYpRITq\nHZNmumY18H4ASQuAsRHxG+Be4GpJYyXNozKt8+gxgg65f1/84hebnsGZnGkk5nKm2v6lNeCRvKQ7\ngYuBEyVtAb4ArABWJJdVHgI+khT3BkmrgA1AD3B9NJLMzMwaNmDJR8Q1x3jqw8fY/kbgxkZDmZlZ\nNnzHa6JUKjU7wlGcqTbOVLuhmMuZ8qXjPaMiybM4ZmZ1kkQcpxOvZmY2TLjkzcwKzCVvZlZgLnkz\nswJLc8ermVlNXnnlFdbcdRd0d9c/eMwYLv2TP2HWrFnZBxtBXPJmlpsDBw6w//HHuWzatLrHrt25\nk/2///s5pBpZXPJmlqsJY8dyWoqSn7B3bw5pRh7PyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG5\n5M3MCswlb2ZWYC55M7MCc8mbmRWYS97MrMAGLHlJKyR1JV/a3fe5/yWpV9L0qnXLJG2S9Jyky/II\nbGZmtRvsSP42YFHflZLmApcCL1atWwhcBSxMxnxDkt8pmJk10YAlHBEPAjv7eer/AJ/rs+4K4M6I\n6I6IDmAzcH4WIc3MLJ26j7QlXQF0RsRTfZ6aBXRWPe4EZjeQzczMGlTXRw1LmgjcQGWq5sjqAYZE\nmlBmZpaNej9P/nSgDXhSEsAc4BeSLgC2AnOrtp2TrDtKe3v7keVSqUSpVKozhplZsZXLZcrlcsP7\nUcTAB9uS2oD7IuLMfp57ATgvIl5LTrx+l8o8/Gzgx8Bbo88PkNR3lZkV1AsvvMB/fPnLLJkzp+6x\nt2/ZwntuuIHTTz89h2TDjyQiYqCZk34NdgnlncBDwAJJWyT9eZ9NjrR1RGwAVgEbgH8Frnebm5k1\n14DTNRFxzSDPn9bn8Y3AjRnkMjOzDPg6djOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnk\nzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MC\nc8mbmRWYS97MrMAG+yLvFZK6JD1dte5/S3pW0pOS7pE0peq5ZZI2SXpO0mV5Bjczs8ENdiR/G7Co\nz7o1wNsj4mxgI7AMQNJC4CpgYTLmG5L8TsHMrIkGLOGIeBDY2Wfd2ojoTR4+AsxJlq8A7oyI7ojo\nADYD52cb18zM6tHokfZHgR8my7OAzqrnOoHZDe7fzMwaMDrtQEl/AxyKiO8OsFn0t7K9vf3IcqlU\nolQqpY1hZlZI5XKZcrnc8H5SlbykPwMuBz5QtXorMLfq8Zxk3VGqS97MzI7W9wB4+fLlqfZT93SN\npEXAZ4ErIuJA1VP3AldLGitpHjAfeDRVKjMzy8SAR/KS7gQuBk6StAX4IpWracYCayUB/Dwiro+I\nDZJWARuAHuD6iOh3usbMzI6PAUs+Iq7pZ/WKAba/Ebix0VBmZpaN1CdezWzkeOzRR9mza1fd43bt\n3g1+Q99ULnkzG9T6NWs4ZcMGJo8bV9e4k4H548fnE8pq4pI3s5qcM3MmsydPbnYMq5M/dsDMrMBc\n8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmG+GMrMh69e//jVvvPFG3eMmTJjA\n/Pnzc0g0/LjkzWxIOm3UKLZ/73u8Xue4fYcOse+MM5h/ww255BpuXPJmNiRdNDvdt4e+/Prr3H/4\ncMZphi/PyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYENWPKSVkjqkvR01brpktZK2ihpjaSpVc8t\nk7RJ0nOSLsszuJmZDW6wI/nbgEV91i0F1kbEAuCB5DGSFgJXAQuTMd+Q5HcKZmZNNGAJR8SDwM4+\nqxcDK5PllcCVyfIVwJ0R0R0RHcBm4PzsopqZWb3SHGnPiIiuZLkLmJEszwI6q7brBNLdzWBmZplo\n6I7XiAhJMdAm/a1sb28/slwqlSiVSo3EMDMrnHK5TLlcbng/aUq+S9LMiNgm6RRge7J+KzC3ars5\nybqjVJe8mZkdre8B8PLly1PtJ810zb3AkmR5CbC6av3VksZKmgfMBx5NlcrMzDIx4JG8pDuBi4GT\nJG0BvgB8FVgl6TqgA/ggQERskLQK2AD0ANdHxEBTOWZmlrMBSz4irjnGU5ccY/sbgRsbDWVmZtnw\ndexmZgXmz5M/jnp6enjllVdSj582bRqtra0ZJjKzonPJH0evv/46X/rSd4iYVffYgwdf49Offj/n\nnXdeDsnMrKhc8sdZb+9k5s69ru5xL710Xw5pzKzoPCdvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF\n5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MCc8mbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmBpS55\nScskPSPpaUnflTRO0nRJayVtlLRG0tQsw5qZWX1SlbykNuAvgHMj4kygBbgaWAqsjYgFwAPJYzMz\na5K0R/J7gG5goqTRwETgZWAxsDLZZiVwZcMJzcwstVQlHxGvAV8DXqJS7rsiYi0wIyK6ks26gBmZ\npDQzs1RSfcerpNOBTwNtwG7gLknXVm8TESEp+hvf3t5+ZLlUKlEqldLEMDMrrHK5TLlcbng/ab/I\n+53AQxGxA0DSPcDvANskzYyIbZJOAbb3N7i65M3M7Gh9D4CXL1+eaj9p5+SfAy6UNEGSgEuADcB9\nwJJkmyXA6pT7NzOzDKQ6ko+IJyXdDjwO9AJPAN8EJgGrJF0HdAAfzCinmZmlkHa6hoi4Gbi5z+rX\nqBzVm5nZEOA7Xs3MCswlb2ZWYC55M7MCc8mbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArM\nJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczK7DUnydvx9+WLVtoaWmpe1xLSwtnnnlmDonM\nbKhzyQ8TY8fO4b77XqTyhVu1O3y4hylTnufrX3fJm41ELvlhYubMc4Bz6h7X3b2PXbtuzT6QmQ0L\nnpM3Myuw1CUvaaqkuyU9K2mDpAskTZe0VtJGSWskTc0yrJmZ1aeRI/mvAz+MiDOAs4DngKXA2ohY\nADyQPDYzsyZJVfKSpgDvjYgVABHRExG7gcXAymSzlcCVmaQ0M7NU0h7JzwNelXSbpCckfUvSCcCM\niOhKtukCZmSS0szMUkl7dc1o4Fzg4xHxmKRb6DM1ExEhKfob3N7efmS5VCpRKpVSxjAzK6ZyuUy5\nXG54P2lLvhPojIjHksd3A8uAbZJmRsQ2SacA2/sbXF3yZmZ2tL4HwMuXL0+1n1TTNRGxDdgiaUGy\n6hLgGeA+YEmybgmwOlUqMzPLRCM3Q30CuEPSWOB54M+BFmCVpOuo3Jr5wYYT2oiyZ88e9uzZk3r8\nrFmzGDXKt3+YvSl1yUfEk8C7+nnqkvRxbKRbt24d37z3m4xvHV/32DFvjOHWr9zqkjer4o81SOH+\n+9fw2GPP1z2up6eH/ftzCFQwo2ePZu65c+se9+J9L+aQxmx4c8mn0NW1my1bzmL69LfWPXbOnPo/\nRdLMLC2XfErjx0+htdW3AZjZ0ObJSzOzAvORvBVGd283d9x1B5LqHnvWGWdxzjvq/yhns6HOJW+F\nMenMSTy87+G6x+18eSdTWqe45K2QXPJWGDNOS3eOpOdQT8ZJzIYOz8mbmRWYS97MrMA8XWM2Qhw+\nfJiennRTU73R7wfK2jDgkjcbITZu3MhdX/saYw4frn/w/v1o9uzsQ1nuXPJmI8iCCK5+y1uaHcOO\nI8/Jm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZgDV0nL6kFeBzojIg/lDQd+Gfg\nLSRf5B0RuxpOacNKRLB58+ZUY3fs2JFxGhuJDhw8yKZNm1KNnTlzJpMmTco4UfM0ejPUp4ANwJu/\nkaXA2oi4WdLnk8dLG/wZNsxEBDf/w83E9HS3wk+YOSHjRDaSjGtp4cSXXuLRm26qe+wrhw5x+dKl\nLFy4MIdkzZG65CXNAS4Hvgz8dbJ6MXBxsrwSKOOSH5G6o5u2i9qaHcNGoBMnTuRDp52WauyqLVsy\nTtN8jRzJ/x3wWWBy1boZEdGVLHcB/hLUIaCnp5tyuZxq7Omnn87cuXOzDWRmx02qkpf0B8D2iFgn\nqdTfNhERkvp9v97e3n5kuVQqUSr1uwvLwKhRYzh48CK+8536x+7Z8zyf+cxYl7xZE5TL5dQHZ9XS\nHsm/G1gs6XJgPDBZ0j8BXZJmRsQ2SacA2/sbXF3ylq+WljG0tZVSjX3xxUPZhjGzmvU9AF6+fHmq\n/aS6hDIiboiIuRExD7ga+ElEfBi4F1iSbLYEWJ0qlZmZZSKr6+TfnJb5KnCppI3A+5PHZmbWJA1/\nnnxE/BT4abL8GnBJo/s0M7Ns+I5XM7MCc8mbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzAqs\n4evkh6uuri5uumkF+/fXP3bfvm5OOOGM7EOZmWVsxJZ8RLBr1xRmzbqu7rEnngijRo3YX52ZDSMj\nuqmkUYwePa7ZMczMcuM5eTOzAnPJm5kVmEvezKzAXPJmZgU2ok+8mr1p9QOr+cFPf1D3uNaxrXxp\n2ZeYNGlSDqnMGueStxFv9v+YTbyt368jHlTnjzozTmOWLZe8jXgaJYSaHcMsF56TNzMrsFQlL2mu\npH+X9IykX0r6ZLJ+uqS1kjZKWiNparZxzcysHmmP5LuB/xkRbwcuBP5K0hnAUmBtRCwAHkgem5lZ\nk6Qq+YjYFhHrk+W9wLPAbGAxsDLZbCVwZRYhzcwsnYbn5CW1AecAjwAzIqIreaoLmNHo/s3MLL2G\nSl5SK/B94FMR8Xr1cxERQLrr0szMLBOpL6GUNIZKwf9TRKxOVndJmhkR2ySdAmzvb2x7e/uR5VKp\nRKlUShvDzKyQyuUy5XK54f2kKnlJAr4NbIiIW6qeuhdYAtyU/Hd1P8N/q+TNzOxofQ+Aly9fnmo/\naY/k3wNcCzwlaV2ybhnwVWCVpOuADuCDKfdvZmYZSFXyEfEzjj2ff0n6OGZmliV/rIEdU0Qvhw8f\n5vDhw3WN6+3tzSmRmdXLJW/HtHXrVm689W5m3FP/lbC7D+7OIZGZ1cslbwNqmdNC2+K2ZscwOz56\ne3m0XGbjunWDb9vHCdOmcenll+cQqjEueTOzxLsmT2b3U0/VPW7voUM8cfLJLnkzs6Fs3rRpqcbt\n2LePJzLOkhV/1LCZWYG55M3MCswlb2ZWYC55M7MC84lXs2HmR/ffz97f/Kbucbv37GGib1QbcVzy\nZsPMpoce4tyuLiaNHVv32MlTpuSQyIYyl7zZMLTgxBM5aeLEZsewYcAlb8fU2wsvrDvArlc21D12\n7AmH+N1r35FDqqHljZ43+Ntb/haNUt1jLz7vYi5fNPRunrFiccnbMZ140jvZvmMcweT6Bkawc+t3\n8og05LS9ry3V95+9+uKr7NyzM+s4ZkdxydsxjR49gdFjJzFmXH0lHzFyTu5NmDQh1bixE+qfTzdL\nw5dQmpkV2LA/kn/ggZ+yefO2uscdPHiAnp4cApmZDSHDvuSfeWYLv/jFLFpbZ9Y9durU8TkkMjMb\nOoZ9yQNMmXIq06e/tdkxzMyGnMxLXtIi4BagBfjHiLgp659htevu7uZnP/8Zh3oP1T22p6eHei+s\nMbOhJdMTr5JagL8HFgELgWsknZHlz8hLR0e52RGOklWm3Qd20zKrpe5/404dx6QTJ/3WvnZtW59J\npix1rO9odoSjDMVMAOWOjmZHOIoz5SvrI/nzgc0R0QEg6XvAFcCzGf+czHV0lGlrKzU7xm+pztTT\n05PqC7K7u7sBGD02mz/1rm3rmTqzlpucgkP763/3ADCqZVRdeTvWd9D2jrZUPysvg2U6fPgw+w/s\nZ+/evXXvuzdSXJifKHd0UGo7dq5mcKZ8ZV3ys4EtVY87gQsy/hkj0sbNG/nVi79i1Kj633x1051D\nooEd3Dedf731xbrH9R7u4ezLxnPaefNySDV07Ny5k3+763Z2lR+qf/C+fejUU7MPZYWkaOCo4Kid\nSX8MLIqIv0geXwtcEBGfqNomsvyZH/vYUu6553mkxm4u2bv3KVpbz8ooVTaqMwXBAfYzenxzz5Xv\n2/UME6e+Pbf9dx/YwcRpz1LPbaR7X91L68mtuWVKY7BM3fu7OeENOKFlzHFMBa/u3cvJrUPrd1WE\nTBHBGWedxb+sXZtbJklERN2fn5F1yV8ItEfEouTxMqC3+uSrpOx+oJnZCDIUSn408CvgA8DLwKPA\nNREx5OfkzcyKKNP3/hHRI+njwI+oXEL5bRe8mVnzZHokb2ZmQ0vuH1AmabqktZI2Slojaeoxtlsm\n6RlJT0v6rqRxQyDTVEl3S3pW0obknENuas2VbNsiaZ2k+5qdSdJcSf+e/P1+KemTOWVZJOk5SZsk\nff4Y2/zf5PknJZ2TR456Mkn6UJLlKUn/KSn3s/u1/J6S7d4lqUfSHw2FTJJKyWv6l5LKeWeqJZek\nkyT9m6T1Sa4/yznPCkldkp4eYJv6XuMRkes/4Gbgc8ny54Gv9rNNG/BrYFzy+J+BJc3MlDy3Evho\nsjwamNLs31XVtn8N3AHc2+xMwEzgHclyK5XzMmdknKMF2Jy8VsYA6/v+DOBy4IfJ8gXAwzn/bmrJ\n9Dtvvm6o3CTY9ExV2/0EuB/442ZnAqYCzwBzkscn5ZmpjlztwFfezATsAEbnmOm9wDnA08d4vu7X\n+PH4qOHFVMqS5L9X9rPNHqAbmJicvJ0IbG1mJklTgPdGxAqonG+IiN05ZqopV5JtDpU/9j8C9X8l\nUcaZImJbRKxPlvdSufltVsY5jtxoFxHdwJs32vWbNSIeAaZKmpFxjroyRcTPq143jwBzcsxTU6bE\nJ4C7gVdzzlNrpj8Fvh8RnQARUf83leeT6xX++8M9JgM7IiK3z6+NiAeBgb5Npu7X+PEo+RkR0ZUs\ndwFHBYqI14CvAS9RuSpnV0T8uJmZgHnAq5Juk/SEpG9JyvtLNWvJBfB3wGeB4/HtHLVmAkBSG5Uj\nkUcyztHfjXaza9gmz1KtJVO164Af5pgHasgkaTaVMvuHZFXeJ+Zq+T3NB6Yn036PS/pwzplqzfUt\n4O2SXgaeBD51HHINpO7XeCZX10haS+Ute19/U/0gIqK/6+QlnQ58msrbpt3AXZI+FBF3NCsTld/N\nucDHI+IxSbcAS4EvpM2URS5JfwBsj4h1kkqNZMkqU9V+WqkcHX4qOaLPUq1F1PedTZ4FVvO+Jb0P\n+CjwnvziALVlugVYmvw9Rf7vBmvJNIbK/28foPJO/ueSHo6ITU3OdQOwPiJKSU+tlXR2RLyeY67B\n1PUaz6TkI+LSY6apnESYGRHbJJ0CbO9ns3cCD0XEjmTMPcC7qcw5NytTJ9AZEY8lj++mUvINySDX\nu4HFki4HxgOTJd0eER9pYiYkjQG+D/y/iFidNssAtgJzqx7PpfI3GmibOeQ77VdLJpKTrd+icjd4\n3l/sWkum84DvVfqdk4Dfk9QdEfc2MdMW4DcRsR/YL+k/gLOBPEu+llzvBr4MEBHPS3oBeBvweI65\nBlL3a/x4TNfcCyxJlpcA/RXAc8CFkiYkRxaXABuamSkitgFbJC1IVl1C5cRQnmrJdUNEzI2IecDV\nwE8aKfgsMiV/s28DGyLilpxyPA7Ml9SmymdYXJVk65v1I0mmC6lM+3WRn0EzSToVuAe4NiI255il\n5kwRcVpEzEteQ3cDH8ux4GvKBPx/4CJVrhqbSOWkYp4dUGuu56j8v08y9/02KheJNEv9r/G8zhJX\nnQ2eDvwY2AisAaYm62cBP6ja7nNUSvRpKicWxgyBTGcDj1GZi7uH/K+uqSlX1fYXk//VNYNmAi6i\ncn5gPbAu+bcohyy/R+XKnc3AsmTdXwJ/WbXN3yfPPwmcexxe3wNmonJyfEfV7+XRZmfqs+1twB8N\nhUzAZ6o64JN5Z6rx73cScF/yenoa+NOc89xJ5bzkISrvbj7a6GvcN0OZmRXY8ZiuMTOzJnHJm5kV\nmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZg/wVSUN8JnvTZ+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b3cb7d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(20): gibbs_step(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonparametric Mixture Models!\n",
    "\n",
    "It turns out, the collapsed Gibbs sampler for mixture models is almost identical in the context of a _nonparametric_ model. This model uses a _Dirichlet process prior_ instead of a _Dirichlet distribution prior_. It doesn't require us to specify how many clusters we are looking for in our data.\n",
    "\n",
    "The cluster assignment score changes slightly. It is proportional to $N_k^{-i}$ for each known cluster. We assign a datapoint to a _new_ cluster with probability proportional to $\\alpha$ (which is now the DP dispersion parameter)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def log_cluster_assign_score_dp(cluster_id, state):\n",
    "    \"\"\"Log-likelihood that a new point generated will\n",
    "    be assigned to cluster_id given the current state.\n",
    "    \"\"\"\n",
    "    if cluster_id == \"new\":\n",
    "        return np.log(state[\"alpha_\"])\n",
    "    else:\n",
    "        return np.log(state['suffstats'][cluster_id].N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predictive likelihood remains the same for known clusters. However, we need to know the likelihood of assigning a datapoint to a new cluster. In this case, we fall back on the hyperparameters to get:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(x_i \\,|\\, z, x_{\\neg_i}, \\sigma)\n",
    "    &= \\mathcal{N}(x_i \\,|\\, \\lambda_1, \\lambda_2^2+\\sigma^2)\n",
    "\\end{align}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def log_predictive_likelihood_dp(data_id, cluster_id, state):\n",
    "    \"\"\"Predictive likelihood of the data at data_id is generated\n",
    "    by cluster_id given the currenbt state.\n",
    "    \n",
    "    From Section 2.4 of \n",
    "    http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf\n",
    "    \"\"\"\n",
    "    if cluster_id == \"new\":\n",
    "        ss = SuffStat(0, 0)\n",
    "    else:\n",
    "        ss = state['suffstats'][cluster_id]\n",
    "        \n",
    "    hp_mean = state['hyperparameters_']['mean']\n",
    "    hp_var = state['hyperparameters_']['variance']\n",
    "    param_var = state['cluster_variance_']\n",
    "    x = state['data_'][data_id]\n",
    "    return _log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this, we can define the marginal distribution over cluster assignment. The only change is that the \"`new`\" state enters in the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cluster_assignment_distribution_dp(data_id, state):\n",
    "    \"\"\"Compute the marginal distribution of cluster assignment\n",
    "    for each cluster.\n",
    "    \"\"\"\n",
    "    scores = {}\n",
    "    cluster_ids = state['suffstats'].keys() + ['new']\n",
    "    for cid in cluster_ids:\n",
    "        scores[cid] = log_predictive_likelihood_dp(data_id, cid, state)\n",
    "        scores[cid] += log_cluster_assign_score_dp(cid, state)\n",
    "    scores = {cid: np.exp(score) for cid, score in scores.iteritems()}\n",
    "    normalization = 1.0/sum(scores.values())\n",
    "    scores = {cid: score*normalization for cid, score in scores.iteritems()}\n",
    "    return scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to be able to create a new cluster when \"`new`\" is drawn, and destroy a cluster if its emptied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_cluster(state):\n",
    "    state[\"num_clusters_\"] += 1\n",
    "    cluster_id = max(state['suffstats'].keys()) + 1\n",
    "    state['suffstats'][cluster_id] = SuffStat(0, 0)\n",
    "    state['cluster_ids_'].append(cluster_id)\n",
    "    return cluster_id\n",
    "\n",
    "def destroy_cluster(state, cluster_id):\n",
    "    state[\"num_clusters_\"] = 1\n",
    "    del state['suffstats'][cluster_id]\n",
    "    state['cluster_ids_'].remove(cluster_id)\n",
    "    \n",
    "def prune_clusters(state):\n",
    "    for cid in state['cluster_ids_']:\n",
    "        if state['suffstats'][cid].N == 0:\n",
    "            destroy_cluster(state, cid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can define the `gibbs_step_dp` function. It's nearly identical to the earlier `gibbs_step` function except\n",
    "* It uses `cluster_assignment_distribution_dp`\n",
    "* It creates a new cluster when the sampled assignment is \"`new`\".\n",
    "* It destroys a cluster any time it is emptied.\n",
    "\n",
    "For clarity, I split out the code for sampling assignment to its own function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def sample_assignment(data_id, state):\n",
    "    \"\"\"Sample new assignment from marginal distribution.\n",
    "    If cluster is \"`new`\", create a new cluster.\n",
    "    \"\"\"\n",
    "    scores = cluster_assignment_distribution_dp(data_id, state).items()\n",
    "    labels, scores = zip(*scores)\n",
    "    cid = random.choice(labels, p=scores)\n",
    "    if cid == \"new\":\n",
    "        return create_cluster(state)\n",
    "    else:\n",
    "        return int(cid)\n",
    "\n",
    "def gibbs_step_dp(state):\n",
    "    \"\"\"Collapsed Gibbs sampler for Dirichlet Process Mixture Model\n",
    "    \"\"\"\n",
    "    pairs = zip(state['data_'], state['assignment'])\n",
    "    for data_id, (datapoint, cid) in enumerate(pairs):\n",
    "        state['suffstats'][cid] = remove_datapoint_from_suffstats(datapoint, state['suffstats'][cid])\n",
    "        prune_clusters(state)\n",
    "        cid = sample_assignment(data_id, state)\n",
    "        state['assignment'][data_id] = cid\n",
    "        state['suffstats'][cid] = add_datapoint_to_suffstats(state['data_'][data_id], state['suffstats'][cid])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we will start by randomly assigning our data to two clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG2VJREFUeJzt3X1wXXd95/H3R8+yZVt24siO48QmQ6DQPLeQElguEDoJ\nBW8mTCmlUC9kGEpLCctsG4fdWWSmnSZpWWi32d0hEMawPKUJCQ7QEidwsy1J80Cc4MQxTpwn+Umy\nZVuy/KSn7/5xT4LGlqV7ru7xlY4+rxmNzj33/M79RLr+5Ojc86CIwMzMZra6WgcwM7Opc5mbmeWA\ny9zMLAdc5mZmOeAyNzPLAZe5mVkOTFrmkm6Q9LSkTZK+LalZ0iJJGyRtlXSvpPZTEdbMzMY3YZlL\nWgF8HLgkIs4H6oEPAmuADRFxHnB/8tjMzGpksi3zfmAImCOpAZgD7ARWAeuSZdYBV2eW0MzMJjVh\nmUfEPuCLwMuUSvxARGwAOiKiO1msG+jINKWZmU1ost0s5wKfAVYAZwJtkj48dpkoXQ/A1wQwM6uh\nhkme/y3gwYjoBZD0feB3gN2SlkTEbklLgZ7xBktyyZuZVSAilGb5yfaZbwEuk9QqScAVwGbgHmB1\nssxq4O4JAk2rr89//vM1z+BM+crlTKWvvr6+Cb/WrFlz0ucOHz48a35O5XxVYsIt84h4UtI3gMeA\nUeBx4CvAPOB2SdcCLwIfqOjVzSw3vvCFW9i3r4nSdt+JHn/8cXp6vnrC/KGhQa655nyuvvr3so6Y\na5PtZiEibgZuPm72Pkpb6WZmABw6BGee+SkaGprHfX7btn6WL//sCfN37HiEkZE9WcfLvVl3Bmih\nUKh1hBM4U/mmYy5nKs+KFYVaRzjBdPw5VUqV7p8pa+VSZLl+M5s+PvnJv+G00z570i3zk9mx4xGu\nvHIP73+/d7O8QhJR5Q9AzcxsBnCZm5nlgMvczCwHXOZmZjngMjczywGXuZlZDrjMzcxywGVuZpYD\nLnMzsxxwmZuZ5YDL3MwsB1zmZmY5MOklcM3MytHb28vw8G7q69NdaKu3t5fe3t6MUs0eLnMzq4rN\nz26mrv0R6uubUo3r27OZc7cdzijV7OEyN7OqiAjmL51HQ1NLqnGHjzRWfKs0+zXvMzczywGXuZlZ\nDkxa5pJeJ2njmK8+SZ+WtEjSBklbJd0rqf1UBDYzsxNNWuYR8auIuDgiLgYuBQ4DdwFrgA0RcR5w\nf/LYzMxqIO1uliuA5yKiC1gFrEvmrwOurmYwMzMrX9oy/yDwnWS6IyK6k+luoKNqqczMLJWyy1xS\nE/A+4J+Ofy5KxxX52CIzsxpJc5z5VcAvImJP8rhb0pKI2C1pKdAz3qDOzs5XpwuFAoVCocKoZmb5\nVCwWKRaLU1qHyj1YX9J3gX+OiHXJ45uB3oi4SdIaoD0i1hw3JnwygNns8MaLrmLhG65LfdLQrmcf\n432XH+Dvbv6rjJLNPJKICKUZU9ZuFklzKX34+f0xs28E3i1pK/DO5LGZmdVAWbtZIuIQcPpx8/ZR\nKngzM6sxnwFqZpYDLnMzsxxwmZuZ5YDL3MwsB1zmZmY54DI3M8sBl7mZWQ64zM3McsBlbmaWAy5z\nM7MccJmbmeWAy9zMLAdc5mZmOeAyNzPLAZe5mVkOuMzNzHLAZW5mlgMuczOzHHCZm5nlQLk3dG6X\ndIekZyRtlvRmSYskbZC0VdK9ktqzDmtmZuMrd8v874EfR8RvABcAW4A1wIaIOA+4P3lsZmY1MGmZ\nS1oAvC0ibgOIiOGI6ANWAeuSxdYBV2eW0szMJlTOlvlKYI+kr0t6XNKtkuYCHRHRnSzTDXRkltLM\nzCbUUOYylwCfiohHJX2Z43apRERIivEGd3Z2vjpdKBQoFAoVhzUzy6NisUixWJzSOhQxbgf/egFp\nCfBQRKxMHr8VuAF4DfCOiNgtaSnws4h4/XFjY7L1m1k+vPGiq1j4hutoaGpJNW7Xs4/xvssP8Hc3\n/1VGyWYeSUSE0oyZdMs8KesuSedFxFbgCuDp5Gs1cFPy/e4KMpvZNNLb28uePXsqGhtUvuF2+PAh\ntmzZUtHYc889l8bGxopfOy/K2c0C8OfAtyQ1AduAjwL1wO2SrgVeBD6QSUIzO2U2b97MLetvYW77\n3NRjhxqHINW2ZEl9Yz3P9G7hb+/82/SD98GXPv8llzlllnlEPAn89jhPXVHdOGZWay1LW1h+yfLU\n4xY8vJ+6uvTnIbbMbWHx6xez/PL0r9n1L12px+SVzwA1M8sBl7mZWQ64zM3McqDcD0AtpTvuuoOu\n7sr2573+Na/nqt+9qsqJzCa3bds2Nj+wm13bBlOPHTws5izIIJSVxWWekU3PbmLXgl20tKU75vbg\n3oM0vOxfi9XGwMAAe14+H9VflHps63xQnd+7teKffIbmL55P26K2VGMiAo5mFMisDPWNc2mdv6zW\nMSwl7zM3M8sBl7mZWQ64zM3McsBlbmaWAy5zM7MccJmbmeWAy9zMLAdc5mZmOeAyNzPLAZe5mVkO\nuMzNzHLAZW5mlgMuczOzHCjrqomSXgT6gRFgKCLeJGkR8D3gHJIbOkfEgYxympnZBMrdMg+gEBEX\nR8SbknlrgA0RcR5wf/LYzMxqIM1uFh33eBWwLpleB1xdlURmZpZami3z+yQ9JunjybyOiOhOpruB\njqqnMzOzspR7p6HLI2KXpMXABklbxj4ZESEpxhvY2dn56nShUKBQKFQY1cwsn4rFIsVicUrrKKvM\nI2JX8n2PpLuANwHdkpZExG5JS4Ge8caOLXMzMzvR8Ru6a9euTb2OSctc0hygPiIOSpoL/C6wFlgP\nrAZuSr7fnfrVzcyAwwcG2LFlR+pxe3btYXBwMINEM085W+YdwF2SXln+WxFxr6THgNslXUtyaGJm\nKc0st5paF7L9mYvZvmXcPbUTOrR9B0eP+g7oUEaZR8QLwEXjzN8HXJFFKDObPVraltDStqSisYd3\nbKxympnLZ4CameWAy9zMLAdc5mZmOeAyNzPLAZe5mVkOuMzNzHLAZW5mlgMuczOzHHCZm5nlgMvc\nzCwHXOZmZjngMjczywGXuZlZDrjMzcxywGVuZpYDLnMzsxxwmZuZ5YDL3MwsB1zmZmY5UFaZS6qX\ntFHSPcnjRZI2SNoq6V5J7dnGNDOziZS7ZX4dsBl45fbZa4ANEXEecH/y2MzMamTSMpd0FvAe4KuA\nktmrgHXJ9Drg6kzSmZlZWcrZMv8S8BfA6Jh5HRHRnUx3Ax3VDmZmZuVrmOhJSe8FeiJio6TCeMtE\nREiK8Z4D6OzsfHW6UChQKIy7GjOzWatYLFIsFqe0jgnLHHgLsErSe4AWYL6kbwLdkpZExG5JS4Ge\nk61gbJmbmdmJjt/QXbt2bep1TLibJSI+FxHLI2Il8EHgpxHxEWA9sDpZbDVwd+pXNjOzqkl7nPkr\nu1NuBN4taSvwzuSxmZnVyGS7WV4VEQ8ADyTT+4ArsgplZmbplF3mls7mTdvo/dVRmuc0pxp36MAh\nRk7T5AuamY3hMs9I965BDujttM5blGpc397t9AxvyyiVmeWVyzxDTXPOoHV+ukPwjwwcBlzmZpaO\ny3waevnFAf7sz25KPa6lBT7zmQ+xfPnyDFKZ2XTmMp9mmucuob3xE7S3vyn12F27vsXo6OjkC5pZ\n7rjMpxmpnvr6FhobWysY6ysam81W/tdvZpYDLnMzsxzwbpYJHDhwgK6urorGjgwPVTmNmdnJucwn\n0NXVxY033kdLy9mpxx7s/01aFjZlkMrM7EQu80k0Ny9n+fL3px733MuLGGry1rmZnRreZ25mlgMu\nczOzHHCZm5nlgMvczCwHXOZmZjngMjczywGXuZlZDrjMzcxyYMIyl9Qi6WFJT0jaLOlvkvmLJG2Q\ntFXSvZLaT01cMzMbz4RlHhFHgXdExEXABcA7JL0VWANsiIjzgPuTx2ZmViOT7maJiMPJZBNQD+wH\nVgHrkvnrgKszSWdmZmWZtMwl1Ul6AugGfhYRTwMdEdGdLNINpLvRpZmZVdWkF9qKiFHgIkkLgJ9I\nesdxz4ekONn4zs7OV6cLhQKFQqHisGZmeVQsFikWi1NaR9lXTYyIPkk/Ai4FuiUtiYjdkpYCPScb\nN7bMzczsRMdv6K5duzb1OiYsc0mnA8MRcUBSK/BuYC2wHlgN3JR8vzv1K9tJHTx4kKc2P5V63L7e\n59m//xLOOeecDFJVV39/Pz+57yeMxEjqsUJc/ubLOeusszJIZjYzTbZlvhRYp9KdguuAb0bE/ZI2\nArdLuhZ4EfhAtjFnj6aWJo60HeGFIy+kHtvX8zKHDh3KIFX1HTlyhPUPrWfOa+ekHtv/Qj9vOO8N\nLnOzMSYs84jYBFwyzvx9wBVZhZrNGpoamHfavIrGHtxeX+U02WpqbmLpeUtTjxvuHc4gjdnM5jsN\nWU3s27ePJ3/+Es91DaQe29/Vzy0v3UlHx7+lHrt0aRt/8icfST3ObLpzmVtNjIyM0LdvEU1n/EHq\nsUNDBzly5FIOHVqcatzg4EEGBu5J/XpmM4HL3GpGqqepdVHqcY1N9cyZs5i2tnSnNxw92sKQb8tq\nOeULbZmZ5YDL3MwsB7ybJWf27t3LSy+9lHpcc3MzS5YsySCRmZ0KLvMciYZRvv2zb7N+4/pU44YH\nh7nwrAv57J9+NqNkZpY1l3mOtC2ex9LLFtK+JN3l5fu6+xjq8SeDZjOZy9xmpN7eXgYHB1ONGRzs\nJ2J3Rommj4jgoYceqmhsd3c3pStd20zjMrcZp35uPVsPbIUD6cYNHTtI69CWbEJNM7fefissTz+u\ne183jS2XVz+QZc5lbjNO26K2isYdHRDDO6scZpoajVFWXroy/cBG2LmtufqBLHMuc7OciQi2bdnO\n/h+lvyLl4f5DEMoglWXNZW6WQ7tehmMt76po7Jz2BVVOY6eCy9wsp+YurGA3i81YLnObktvvvJ1H\nn3409bj+/n6GR30pW7NqcZnblOw5sIf+Jf20L013bHvj/kbm/bKyDzLN7EQuc5uyptYmWue1phoz\nMjRCXcPRjBKZzT6+0JaZWQ64zM3McmDS3SySlgPfAM4AAvhKRPyDpEXA94BzSG7qHBEpz8mb3g4d\nOkR/fx+9vb2px46MpD/G18zSe/DBB1m4cGHqcR0dHVxwwQUZJKqNcvaZDwH/OSKekNQG/ELSBuCj\nwIaIuFnS9cCa5Cs3duzYwcYtT/D8vnS3JwMYGR2hrc4f8Jllabh+Jf/tpvtTjxs8eogrC2fwf275\nHxmkqo1JyzxKVybanUwPSHoGWAasAt6eLLYOKJKzMgeoa65j/tnzax3DzMZx9oWrKhrX8+ImRkaf\nr3Ka2kp1NIukFcDFwMNAR0R0J091A+luyGi5sGP7DrYf287B3oOpxg0eGQTmZBPKbBYqu8yTXSx3\nAtdFxEHp19dviIiQFOON6+zsfHW6UChQKBQqzWrT0HNbe9n07BtoW5R+V1RdQ0sGicxmnmKxSLFY\nnNI6yipzSY2UivybEXF3Mrtb0pKI2C1pKdAz3tixZW75NGfBCtqXvLbWMcxmrOM3dNeuXZt6HZMe\nmqjSJvjXgM0R8eUxT60HVifTq4G7jx9rZmanRjlb5pcDHwZ+KWljMu8G4EbgdknXkhyamElCMzOb\nVDlHs/wbJ9+Cv6K6cczMrBI+A9TMLAdc5mZmOTArrprY19fH6Oho6nFHjhzJIE22Bo8McuRgutxH\nB45y7OixjBKZ2akwK8r8i//ri3Qd6EJ16e5tuK9nH9SfnVGq6hsdbueh218GBlKNO3b4KOee1pdN\nKDM7JWZFmR8ZOkLH2zpoaUt3kkrz8828tGvmXF9l4ZmVfR7d1/Myw8P3VTmNmZ1K3mduZpYDLnMz\nsxxwmZuZ5cCs2Gdu9oqh4WF++OMfVjT2wvMvZPny5VVOZFYdLnObNerqxWjbCHc9f1fqsX3b+7h+\n8fUuc5u2XOY2a9TV19F2+jyW/2b6Qh49nP48BbNTaVaU+VMbt3Hk+XoamtL95w4fG2R05LSMUpmZ\nVc+sKPO9PaOMDL6Hptb0x4y3LWzNIJGZWXXNijIHaGo9nZa29lrHMDPLhA9NNDPLAZe5mVkOzJrd\nLDaxwaFBnnrqqdTjjh07mkEaM0vLZW7U1Tfw/Av1XPPhv049dmSkjvlnNWWQyszSmLTMJd0G/B7Q\nExHnJ/MWAd8DziG5/2dEHMgwp2Vo3mlnMu+0j9c6hplNQTlb5l8H/ifwjTHz1gAbIuJmSdcnj9dk\nkM+sqkZHhtm1dVfqcb07etm/f38Gicyqo5wbOv+rpBXHzV4FvD2ZXgcUcZnbNKe6Bo4OXMwjd6U/\nm/PAzsPsuiT9/wTMTpVK95l3RER3Mt0NdFQpj1lm6htaWLj0P1Q0dmDP4SqnMauuKR+aGBEBRBWy\nmJlZhSrdMu+WtCQidktaCvScbMHOzs5XpwuFAoVCocKXNDPLp2KxSLFYnNI6Ki3z9cBq4Kbk+90n\nW3BsmZuZ2YmO39Bdu3Zt6nVMuptF0neAB4HXSeqS9FHgRuDdkrYC70wem5lZjZRzNMsfnuSpym4F\nb2ZmVeczQM1sVoqAkZGRisbW1dUhqcqJpsZlbmaz0tOb9nPttekvYdHaOsratX/KGWeckUGqyrnM\nzWzWmTN/JR2Nb2PFiotSj+3quiWDRFPnMjezWWloaIi+vr7U4wYGBhgaGsog0dS4zM1s1qlvqGdH\n7w52btyZeuyhvRs5cOC9LFu2LINklXOZm9ms0zq/ldb5ld3f99A0vT6sy9xsmnrppZfYuTP9lmPp\nChs227jMzaapLb/awld++hXmnz4/1biIgDkZhbJpy2VuNo0tWLaAs88/O9WYiGDjfdN0X4BlxmVu\nVoYANj61kfh2+l0Yc1vm8vvX/H71Q5mN4TI3K0PzvBa6WrZx7OixVONGhkZo2NTgMrfMuczNytDY\n0sCiZYvoeE26+7AMDw6z78V9GaUy+7UZU+abNm3iO/d8h5FIfy2FgSMDLMggk1mWNm16mkf/eSdb\nHxxINS4Iho42ZZTKpqsZU+bHjh1je2xn6flLU4+d99g86hvqM0hllp1jxwY53HcZ8xdfknrswjMz\nCGTT2owpc4CGpoaKDvRvaG5AddPrCmc284wOjzJ0LN1p3MODwwwPDVf8mnX1zTQ2pzs00WanGVXm\nZrUSo808/qNDwLZU40ZHRqjr780mlNkYLnOzMixc+mbgzanHDQ8eZf/ev69+ILPjTHrbODMzm/5c\n5mZmOTClMpd0paQtkp6VdH21QpmZWToV7zOXVA/8I6UbO+8AHpW0PiKeqVa4sfbs2cPubbsZbRxN\nPXZkaJi65tL0gd1P0L4k/d1FsuRM5ZuOuSbLFKPBAw88kHq9PT09wOJMMtWCM2VrKh+Avgl4LiJe\nBJD0XeA/ApmV+eafL2PHi69PPzhEU2sjMD1/ec5UvumYa8JMgpG5I9z2i9tSr3fH4d00tby7+plq\nxJmyNZUyXwZ0jXm8nUo+7k+hsWUhC5dcmuVLmFVVXV0d8xfPZ+UlK1OPHRmG3q7mDFJZHk2lzE/5\nFfCHBp5n56bvTWkdB7ufmvI6qs2Zyjcdc02UaXR0iIbGvfzgCz9Ivd4YFb27+jnUs7mqmWolL5li\nOP19Q08FVXpXEkmXAZ0RcWXy+AZgNCJuGrOMb3liZlaBiEh12vpUyrwB+BXwLmAn8Ajwh1l9AGpm\nZidX8W6WiBiW9CngJ0A98DUXuZlZbVS8ZW5mZtNHVc8AlbRI0gZJWyXdK6n9JMvdIOlpSZskfVtS\nZh/Zp8jULukOSc9I2px8JlDTTMmy9ZI2SronqzzlZpK0XNLPkt/dU5I+nVGWSU9Gk/QPyfNPSro4\nixxpc0n6oyTPLyX9XNIFtc40ZrnfljQs6ZrpkElSIXlfPyWpWOtMkk6X9C+Snkgy/aeM89wmqVvS\npgmWSfcej4iqfQE3A3+ZTF8P3DjOMiuA54Hm5PH3gNXVzJE2U/LcOuBjyXQDsKDWmZLnPwt8C1if\nVZ4Uv7slwEXJdBulz0x+o8o56oHnkvdJI/DE8a8BvAf4cTL9ZuDfs/zZpMj1O6+8b4Ars85VTqYx\ny/0U+CHw/lpnAtqBp4GzksenT4NMncDfvJIH6AUaMsz0NuBiYNNJnk/9Hq/2tVlWUSpFku9Xj7NM\nPzAEzEk+RJ1D6QzSrEyaSdIC4G0RcRuUPg+IiCyPPyrn54Sksyj9Ur8KZH1B9kkzRcTuiHgimR6g\ndIJYtW+D8OrJaBExBLxyMtq4WSPiYaBdUrr7uWWQKyIeGvO+eRg4q9aZEn8O3AHsyThPuZk+BNwZ\nEdsBImLvNMi0C3jlwvHzgd6IqPxC9JOIiH8F9k+wSOr3eLXLvCMiupPpbuCEF4+IfcAXgZcpHQVz\nICLuq3KOVJmAlcAeSV+X9LikWyXNqXEmgC8BfwGkv4ZBdpkAkLSC0pbFw1XOMd7JaMvKWCbr4iwn\n11jXAj/ONFEZmSQto1Rc/zuZlfWHZOX8nF4LLEp22T0m6SPTINOtwBsl7QSeBK7LONNkUr/HUx/N\nImkDpT+3j/dfxz6IiBjvOHNJ5wKfofQnTx/wT5L+KCK+lTZLtTJR+jlcAnwqIh6V9GVgDfDfa5VJ\n0nuBnojYKKlQaY5qZhqznjZKW3rXJVvo1VRu2Rz/l0rWJVX2+iW9A/gYcHl2cYDyMn0ZWJP8TkX2\nf+GVk6mR0r+3d1H6y/whSf8eEc/WMNPngCciopB01AZJF0bEwYwylSPVezx1mUfESS8WkezQXxIR\nuyUtBXrGWey3gAcjojcZ833gLZT2C1ekCpm2A9sj4tHk8R2UyrxiVcj0FmCVpPcALcB8Sd+IiD+u\nYSYkNQJ3Av83Iu6uNMsEdgDLxzxeTun3M9EyZ5Htrrpyc5F86HkrcGVETPRn9KnKdCnw3VKPczpw\nlaShiFhfw0xdwN6IOAIckfT/gAuBrMq8nExvAf4aICK2SXoBeB3wWEaZJpP6PV7t3SzrgdXJ9Gpg\nvH/sW4DLJLUmWwpXAOnPV65ipojYDXRJOi+ZdQWlD2hqmelzEbE8IlYCHwR+OpUir0am5Pf1NWBz\nRHw5oxyPAa+VtEJSE/AHSbbjs/5xkukySrvqusnWpLkknQ18H/hwRDyXcZ6yMkXEayJiZfI+ugP4\nZIZFXlYm4AfAW1U6UmsOpQ/4suyAcjJtofTvnmTf9OsoHahRK+nf41X+hHYRcB+wFbgXaE/mnwn8\naMxyf0mpLDdR2snfWM0cFWa6EHiU0v6y75Pt0SxlZRqz/NvJ/miWSTMBb6W0//4JYGPydWUGWa6i\ndKTMc8ANybxPAJ8Ys8w/Js8/CVyS5c+m3FyUPqjuHfOzeaTWmY5b9uvANdMhE/BfxnTAp2udidJf\nLfck76dNwIcyzvMdSp8ZDlL6S+VjU32P+6QhM7Mc8G3jzMxywGVuZpYDLnMzsxxwmZuZ5YDL3Mws\nB1zmZmY54DI3M8sBl7mZWQ78f0vr+1GshcH3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b3cead0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "state = initial_state(num_clusters=2, alpha=0.1)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what happens when we run our Gibbs sampler once.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGphJREFUeJzt3X2QXXWd5/H3J/2QR/JEsumERIMISlyfQCXjw3JVxAhs\nsGYcHhw0jtTM1Lg+7bgqYaq0M1Uy4tTMMLNTriUaBmcEiUCxwUci7nV1JAhKSEyIJEBCmqQ7PCQh\nD510d/q7f9xD9trpdPc5956+nZPPqyrFueee37mf6r794dzzcI8iAjMzK6ZxjQ5gZmb5ccmbmRWY\nS97MrMBc8mZmBeaSNzMrMJe8mVmBDVnyklZK6pK0YcD8T0h6TNJvJd1YNX+5pC2SNku6OK/QZmY2\nMs3DPH8L8D+Bb700Q9I7gaXA6yKiV9LsZP4i4EpgEXAG8BNJ50REfy7JzcxsWENuyUfEz4E9A2b/\nJfC3EdGbLPNsMv9y4PaI6I2IbcBW4C31jWtmZmlk2Sd/NvBfJK2VVJb0pmT+PKCjarkOKlv0ZmbW\nIMPtrjnRmBkRsVjSm4FVwCtOsKy/M8HMrIGylHwHcDdARDwkqV/SLOAZYEHVcvOTeb9HkovfzCyD\niFDaMVl219wDvAtA0jlAa0Q8B6wGrpLUKulMKrt1fnWCoGPu3xe/+MWGZ3AmZzoVcznTyP5lNeSW\nvKTbgQuB0yXtAL4ArARWJqdV9gAfTop7k6RVwCagD/hY1JLMzMxqNmTJR8TVJ3jqQydY/gbghlpD\nmZlZffiK10SpVGp0hOM408g408iNxVzOlC+N9h4VSd6LY2aWkiRilA68mpnZScIlb2ZWYC55M7MC\nc8mbmRVYlitezcxGZNeuXXz3+9+nN8PYFuADl1zCvHnz6h3rlOKSN7PcHD58mIeff54ZixenHrvn\nwQe5tLs7h1SnFpe8meWqdcIEZrzsZanHHXj00RzSnHq8T97MrMBc8mZmBeaSNzMrMJe8mVmBueTN\nzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgU2ZMlLWimpK7lp98DnPiOpX9LMqnnLJW2RtFnS\nxXkENjOzkRtuS/4WYMnAmZIWAO8BtlfNWwRcCSxKxnxVkj8pmJk10JAlHBE/B/YM8tQ/AJ8bMO9y\n4PaI6I2IbcBW4C31CGlmZtmk3tKWdDnQERHrBzw1D+ioetwBnFFDNjMzq1GqrxqWNAm4nsqummOz\nhxgSWUKZmVl9pP0++bOAhcCjkgDmA7+WdAHwDLCgatn5ybzjtLe3H5sulUqUSqWUMczMiq1cLlMu\nl2tejyKG3tiWtBC4NyJeO8hzTwHnR8QLyYHX26jshz8D+AnwyhjwApIGzjKzgnrqqaf40urVzF+6\nNPXYHffey/WXXspZZ52VQ7KTjyQiYqg9J4Ma7hTK24FfAudI2iHpTwcscqytI2ITsArYBPwQ+Jjb\n3MyssYbcXRMRVw/z/CsGPL4BuKEOuczMrA58HruZWYG55M3MCswlb2ZWYC55M7MCc8mbmRWYS97M\nrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCX\nvJlZgbnkzcwKzCVvZlZgw93Ie6WkLkkbqub9naTHJD0q6W5J06qeWy5pi6TNki7OM7iZmQ1vuC35\nW4AlA+bdB7wmIl4PPA4sB5C0CLgSWJSM+aokf1IwM2ugIUs4In4O7Bkwb01E9CcPHwTmJ9OXA7dH\nRG9EbAO2Am+pb1wzM0uj1i3tjwI/SKbnAR1Vz3UAZ9S4fjMzq0Fz1oGS/hroiYjbhlgsBpvZ3t5+\nbLpUKlEqlbLGMDMrpHK5TLlcrnk9mUpe0keAS4B3V81+BlhQ9Xh+Mu841SVvZmbHG7gBvGLFikzr\nSb27RtIS4LPA5RFxuOqp1cBVklolnQmcDfwqUyozM6uLIbfkJd0OXAjMkrQD+CKVs2lagTWSAB6I\niI9FxCZJq4BNQB/wsYgYdHeNmZmNjiFLPiKuHmT2yiGWvwG4odZQZmZWHz6P3cyswFzyZmYF5pI3\nMyswl7yZWYFlvhjKzE4dTz75JN3d3anH7d69O4c0loZL3syGtepHP2JjTw+tEyemHzxnTv0D2Yi5\n5M1sWH0RzL7gAqa2tTU6iqXkffJmZgXmkjczKzCXvJlZgbnkzcwKzAdezWzMevLJJzl48GDqcRMn\nTuTss8/OIdHJxyVvZmPSuLlz+c62bbBtW6pxPd3dnDtuHNe75AGXvJmNUWe86U2Zxu3v6uLoww/X\nOc3Jy/vkzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswIYseUkrJXVJ2lA1b6akNZIel3SfpOlVzy2X\ntEXSZkkX5xnczMyGN9yW/C3AkgHzrgPWRMQ5wP3JYyQtAq4EFiVjvirJnxTMzBpoyBKOiJ8DewbM\nXgrcmkzfCrw/mb4cuD0ieiNiG7AVeEv9opqZWVpZtrTnRERXMt0FvHRHgHlAR9VyHcAZNWQzM7Ma\n1XTFa0SEpBhqkcFmtre3H5sulUqUSqVaYpiZFU65XKZcLte8niwl3yWpLSI6Jc0FXrqJ4zPAgqrl\n5ifzjlNd8mZmdryBG8ArVqzItJ4su2tWA8uS6WXAPVXzr5LUKulM4GzgV5lSmZlZXQy5JS/pduBC\nYJakHcAXgC8DqyRdC2wDrgCIiE2SVgGbgD7gYxEx1K4cMzPL2ZAlHxFXn+Cpi06w/A3ADbWGMjOz\n+vB57GZmBebvkx9FfX197Nq1K/P4GTNmMGXKlDomMrOic8mPov379/Ovf/M3zMtwqOKFI0d416c/\nzfnnn59DMjMrKpf8KJva38+1CxYMv+AA9z79dA5pzKzovE/ezKzAXPJmZgXmkjczKzCXvJlZgbnk\nzcwKzCVvZlZgLnkzswJzyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MC\ny1zykpZL2ihpg6TbJI2XNFPSGkmPS7pP0vR6hjUzs3QylbykhcCfAedFxGuBJuAq4DpgTUScA9yf\nPDYzswbJuiX/ItALTJLUDEwCdgJLgVuTZW4F3l9zQjMzyyxTyUfEC8DfA09TKfe9EbEGmBMRXcli\nXcCcuqQ0M7NMMt3jVdJZwKeBhcA+4LuSrqleJiJC0qB3rG5vbz82XSqVKJVKWWKYmRVWuVymXC7X\nvJ6sN/J+E/DLiHgeQNLdwB8AnZLaIqJT0lxg92CDq0vezMyON3ADeMWKFZnWk3Wf/GZgsaSJkgRc\nBGwC7gWWJcssA+7JuH4zM6uDTFvyEfGopG8BDwP9wG+ArwOnAaskXQtsA66oU04zM8sg6+4aIuIr\nwFcGzH6Byla9mZmNAb7i1cyswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MCc8mbmRWY\nS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmCZv2rYRtnRo/zinntYd999qYeOP/10rvnz\nP88hlJmNdS75k8Tb2tp4/YEDcOBAqnFH+vq4u6Mjp1RmNta55E8SMydOZObEianHHerthb17c0hk\nZicD75M3Myswl7yZWYFlLnlJ0yXdKekxSZskXSBppqQ1kh6XdJ+k6fUMa2Zm6dSyJf9PwA8i4lzg\ndcBm4DpgTUScA9yfPDYzswbJVPKSpgHviIiVABHRFxH7gKXArclitwLvr0tKMzPLJOuW/JnAs5Ju\nkfQbSTdLmgzMiYiuZJkuYE5dUpqZWSZZT6FsBs4DPh4RD0m6iQG7ZiIiJMVgg9vb249Nl0olSqVS\nxhhmZsVULpcpl8s1rydryXcAHRHxUPL4TmA50CmpLSI6Jc0Fdg82uLrkzczseAM3gFesWJFpPZl2\n10REJ7BD0jnJrIuAjcC9wLJk3jLgnkypzMysLmq54vUTwLcltQJPAH8KNAGrJF0LbAOuqDmhmZll\nlrnkI+JR4M2DPHVR9jhmZlZP/u6aDO773vd44qGHhl9wgL6+PujuziGRmdngXPIZ7Ovq4nU7dvDK\nmTNTj22aPz+HRGZmg3PJZzRtwgTmTJnS6BhmZkPyF5SZmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmB\nueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZg/j75\nU0DPwYN852tfSz9Q4g0XXsirX/3q+ocys1FRU8lLagIeBjoi4r9KmgncAbyc5EbeEbG35pSWWWtT\nEx+YPBk2bEg9dt1zz/HCa1+bQyqzfB05fJgtW7ZkGtvW1sZpp51W50SNU+uW/KeATcBLP5HrgDUR\n8RVJn08eX1fja1gNmseN49zZszON3XHoUJ3TmOWvqbWVp8eP58Yf/zj12J7nnuO6K65g0aJFOSRr\njMwlL2k+cAnwJeCvktlLgQuT6VuBMi55S2Ht2rXc8cM7INKPnTJ+Cl/43BdoaWmpf7AC2LJlCzff\neSe9kf6H+9yLLzL7vPNySFV/k2bM4BWXXZZp7I4M/2MY62rZkv9H4LPA1Kp5cyKiK5nuAubUsH47\nBR05coRnJz3Lgv+8IPXYHffvyCFRcfT19bF78mTmlkqpx84HWiZNqnsmy1+mkpd0GbA7Ih6RVBps\nmYgISYNuMrS3tx+bLpVKlDK86ay4mlqaGD95fOpxknJIUyxNzc2ML9D+5iIrl8uUy+Wa15N1S/6t\nwFJJlwATgKmS/g3oktQWEZ2S5gK7BxtcXfJmZna8gRvAK1asyLSeTOfJR8T1EbEgIs4ErgJ+GhEf\nAlYDy5LFlgH3ZEplZmZ1Ua+LoV7aLfNl4D2SHgfelTw2M7MGqfliqIj4GfCzZPoF4KJa12lmZvXh\nrzUwMyswl7yZWYH5u2usUPbu3Utzc/q39fjx45nk88CtgFzyVhgHxx3k8//w+dTjent6+eN3/jFL\nL1uaQyqzxnLJW2EsujTb9410bOqgP/rrnMZsbPA+eTOzAnPJm5kV2Cm7u6arq4uVN94I3d2px/Ye\nOsS5kyfnkMrMrL5O2ZKPCKbt3cu18+alH3z66TSP84cgMxv7TtmSBxgnMT7D6XZmZicLb46amRWY\nS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArM5w9a3UUEu3cPenvfYR04cKDOacxObS55q7uIoP3v\n2umZ2JNpfOu81jonMjt1ZSp5SQuAbwH/icr9Xb8eEf8saSZwB/ByYBtwRUTsrVNWO4l0H+1m4bsX\nNjqG2Skv6z75XuC/R8RrgMXAf5N0LnAdsCYizgHuTx6bmVmDZCr5iOiMiHXJ9AHgMeAMYClwa7LY\nrcD76xHSzMyyqfnsGkkLgTcCDwJzIqIreaoLmFPr+s3MLLuaDrxKmgLcBXwqIvZLOvZcRISkqDGf\n2ajYvGUzq+5elXpc07gmLl1yKRMmTMghlVntMpe8pBYqBf9vEXFPMrtLUltEdEqaCwx6Hl17e/ux\n6VKpRKlUyhrDcrZv3z527tyZakx/fz+HDh3KKdGJ7X9+f+U0gJRaxrewuWczTz77ZOqxh393mIvf\nfbFL3uquXC5TLpdrXk/Ws2sEfBPYFBE3VT21GlgG3Jj8955Bhv9eydvYNbWpifV33cWOu+5KNe5o\nfz9bHlvPoquz3XM1q1/c9jg9h+YipdsL2dd7mEWloyx8Q/p7Czz9xNOpx5iNxMAN4BUrVmRaT9Yt\n+bcB1wDrJT2SzFsOfBlYJelaklMoM67fxoDF8+ezOMO4vv5+vv7bvrrnGfZ1jzRz2uwrGDcu3dv6\nxWc3QazNKZVZY2Uq+Yj4BSc+aHtR9jhmZlZPvuLVxpToDyKyHa/3UX6z47nkbUzZvn476364D9Cw\nyw505NAEpsysfyazk5lL3sac0NuZ0XZ+o2PYKag/gvIDD/DIpk2px8447TQuee97c0hVG5e8mVli\n6mtew/r9+1OP6zl4kNnr17vkzczGshkLFmQad2jPHti+vc5p6sMlb3XXH0Hnzj5++o30H3l7e3og\n0u+PN7PBueQtF0cOT6Wv74Opx2kcTJ01PodEZqcml7zlRLROmN7oEGanPN/j1cyswFzyZmYF5pI3\nMyswl7yZWYH5wKvZKeLpp5/m16vX8rv1z6Ue29QavPma9zFl9uwcklmeXPJmp4je3l4O7pnP5Jdf\nk3rsvp130n/0aA6pLG8uecvN0YylIIlx40Z3T2L/0X56j/SmHtfX20d3dzfNzen/lJqammhtbU09\nrhZqaqZlwoz048Y15ZDGRoNL3nKh/n6e2bIl9biIYGpbGzNmpC+irKQmHn/gKFvWPpF67Iu7tvOZ\n/s/Q0tqSatzRo0e55K2XcOUHrkz9mmZpnPQl/7P776dz69bU4w4fOQJ9o3/3opPJbzp38dBzz6Qe\n1x/BtNYWXpbhvqd7Dx8e9e+FP23Wq4BXZRx9M/PfO5/xk9Jdpdu5tZO+fr//LH8nfcnv2LiReb/+\nNW1TpqQeO2G6r8gcyq6D+/lp0zZmzZqYalxE0DrdH+/zsn37dnp70+9a2rlzJ309PTkksrHspC95\ngJdNm8YrZ/puEXmYMqWF2bMnpxoTEUyY4DLJy8rvrGT7oe00t6T7831i4xMc3HdOTqlsrKp7yUta\nAtwENAHfiIgb6/0aZqeyO+74IXuaxLiWdJ+Weg8chn5/wjrV1LXkJTUB/0LlZt7PAA9JWh0Rj9Xz\ndfJQ3raN0sKFjY7xe8Zipm3b9rJw4djazbW3cx3T297QsNfv2NRBU/Pvl2fn1k7aXtl2wjEvPPMC\nO+ftzPR6Bw400T/vGlqmzUo1rn/8YcaN68z0mt179rDjkUeYnOET84SpU2k799wTPr9t7VoWLl6c\nKVdexmKmrOq9Jf8WYGtEbAOQ9B3gcsAln8FYzOSSH2gx639y/K6pzq07aXvla084av+zu5h7wU4O\nHDiQ+hWDoKl1Cs0TpqYa19/TwpH95dSvB7DvuU7WPfofjJ+V7jWPHjrCvKmvc8k3UL1L/gxgR9Xj\nDuCCOr+G2Zgxbc7gRb6vc8OQ96ntO7yRTRs38slPfi31ax45OJOWDDc6r0X0nca+306laUK6U0V7\nDx0lTn8CPpJPLhueIup3wpqkPwKWRMSfJY+vAS6IiE9ULRP1fM1rr7qKH917L+NU25t+X08P00b5\nwpThNDpTAAcn9NLS+v8vTDq4v5fJpw3/h97bMxd4derX7It++pqaSdNhvQcep2XK2DqgOFym/p5+\nWvqaM90F66iOcrQ1UEu6C8b6e/uJA0/QMv7EW9UnEuqnv+Uoak6Xt//wHpr1FE3qPuEyPd3dtE5M\ndwZX3tJmighefe65/PqXv8wtkyQi0r9h6l3yi4H2iFiSPF4O9FcffJU02qdBm5kVwlgo+Wbgd8C7\ngZ3Ar4CrT4YDr2ZmRVTXffIR0Sfp48CPqZxC+U0XvJlZ49R1S97MzMaW3L/qT9JMSWskPS7pPkmD\nnn8nabmkjZI2SLpNUrovA8kn03RJd0p6TNKm5JhDbkaaK1m2SdIjku5tdCZJCyT9n+T391tJn8wp\nyxJJmyVtkfT5Eyzzz8nzj0p6Yx450mSS9CdJlvWS/kPS6xqdqWq5N0vqk/SHYyGTpFLynv6tpHLe\nmUaSS9IsST+StC7J9ZGc86yU1CVpwxDLpHuPR0Su/4CvAJ9Lpj8PfHmQZRYCTwLjk8d3AMsamSl5\n7lbgo8l0MzCt0T+rqmX/Cvg2sLrRmYA24A3J9BQqx2XOrXOOJmBr8l5pAdYNfA3gEuAHyfQFwNqc\nfzYjyfQHL71vgCVjIVPVcj8Fvgf8UaMzAdOBjcD85PGsPDOlyNUO/O1LmYDngeYcM70DeCOw4QTP\np36Pj8aXdi+lUpYk/33/IMu8CPQCk5KDt5OoXDHbsEySpgHviIiVUDneEBH7csw0olxJtvlUftnf\nINXJhvlkiojOiFiXTB+gcvHbvDrnOHahXUT0Ai9daDdo1oh4EJguaU6dc6TKFBEPVL1vHgTm55hn\nRJkSnwDuBJ7NOc9IM30QuCsiOgAiIv3tq/LJtQt46QqwqcDzEZHb14dGxM+BPUMskvo9PholPyci\nupLpLuC4QBHxAvD3wNNUzsrZGxE/aWQm4EzgWUm3SPqNpJslTcox00hzAfwj8FmgP+c8aTIBIGkh\nlS2RB+ucY7AL7c4YwTJ5lupIMlW7FvhBjnlgBJkknUGlzP5XMivvA3Mj+TmdDcxMdvs9LOlDOWca\naa6bgddI2gk8CnxqFHINJfV7vC5n10haQ+Uj+0B/Xf0gImKw8+QlnQV8msrHpn3AdyX9SUR8u1GZ\nqPxszgM+HhEPSboJuA74QtZM9cgl6TJgd0Q8IqlUS5Z6ZapazxQqW4efSrbo62mkRTTwk02eBTbi\ndUt6J/BR4G35xQFGlukm4Lrk9yny/zQ4kkwtVP7e3k3lk/wDktZGRPo7z9Q31/XAuogoJT21RtLr\nI2J/jrmGk+o9XpeSj4j3nDBN5SBCW0R0SpoL7B5ksTcBv4yI55MxdwNvpbLPuVGZOoCOiHgoeXwn\nlZKvSR1yvRVYKukSYAIwVdK3IuLDDcyEpBbgLuDfI+KerFmG8AywoOrxAiq/o6GWmU++u/1Gkonk\nYOvNVK4GH+qj+GhlOh/4TqXfmQW8T1JvRKxuYKYdwHMR0Q10S/q/wOuBPEt+JLneCnwJICKekPQU\nlTvMPJxjrqGkfo+Pxu6a1cCyZHoZMFgBbAYWS5qYbFlcBGxqZKaI6AR2SHrp2vSLqBwYytNIcl0f\nEQsi4kzgKuCntRR8PTIlv7NvApsi4qaccjwMnC1poaRW4Mok28CsH04yLaay26+L/AybSdLLgLuB\nayIi/S3McsgUEa+IiDOT99CdwF/mWPAjygT8b+Dtqpw1NonKQcU8O2CkuTZT+dsn2ff9KioniTRK\n+vd4XkeJq44GzwR+AjwO3AdMT+bPA75ftdznqJToBioHFlrGQKbXAw9R2Rd3N/mfXTOiXFXLX0j+\nZ9cMmwl4O5XjA+uAR5J/S3LI8j4qZ+5sBZYn8/4C+IuqZf4lef5R4LxReH8PmYnKwfHnq34uv2p0\npgHL3gL84VjIBPyPqg74ZN6ZRvj7mwXcm7yfNgAfzDnP7VSOS/ZQ+XTz0Vrf474YysyswEZjd42Z\nmTWIS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAvt/RJGFSN+QcUkAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10af45190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gibbs_step_dp(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We went from 2 to 4 clusters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After 100 iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGKlJREFUeJzt3X2QHPV95/H3Bz2ApNUzIKEHWwJLKYkyBDAPwXaYxEAR\nOSe4JGUgflBsKsmZM4ZLYlsiVfbKKYjxxRcSp3xX9lmUyGEcgTkFMIkl44xNgkE8CJBZEQmDQAvS\n8iQhLXra1X7vj2l0k9Vqd6ane2e39/Oq2qKnp389n9odfej5dc+MIgIzMyum45odwMzM8uOSNzMr\nMJe8mVmBueTNzArMJW9mVmAueTOzAuu35CWtktQhaVOv9ddJ2izpF5JuqVq/QtJWSc9JujSv0GZm\nVpvRA9x/G/BN4PZ3V0j6DWApcEZEdEk6KVm/GLgSWAzMBn4saWFE9OSS3MzMBtTvkXxEPATs6rX6\ns8BfRkRXss3ryfrLgTsjoisitgHPA+dlG9fMzOqRZk5+AfDrkh6RVJb0gWT9LKC9art2Kkf0ZmbW\nJANN1xxrzNSIuEDSucAa4NRjbOvPTDAza6I0Jd8O3AMQEY9J6pF0IvAKMLdquznJuv9AkovfzCyF\niFC9Y9JM16wFfhNA0kJgbES8AdwLXCVprKT5VKZ1Nhwj6JD7+cpXvtL0DM7kTCMxlzPV9pNWv0fy\nku4ELgKmS9oOfBlYBaxKLqs8BHwqKe42SWuANqAbuDYaSWZmZg3rt+Qj4upj3PXJY2x/M3Bzo6HM\nzCwbfsdrolQqNTvCUZypNs5Uu6GYy5nypcGeUZHkWRwzszpJIgbpxKuZmQ0TLnkzswJzyZuZFZhL\n3syswNK849XMrCY7duxg3V13QVdX/YPHjOGS3/s9Zs2alX2wEcQlb2a5OXDgAPsff5xLp06te+z6\nXbvY/9GP5pBqZHHJm1muxo0dy6kpSn5cZ2cOaUYez8mbmRWYS97MrMBc8mZmBeaSNzMrMJe8mVmB\nueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzA+i15SaskdSRf2t37vj+V1CNpWtW6FZK2SnpO0qV5\nBDYzs9oNdCR/G3BZ75WS5gKXAC9VrVsMXAksTsZ8S5JfKZiZNVG/JRwRDwG7+rjrfwBf7LXucuDO\niOiKiG3A88B5WYQ0M7N06j7SlnQ50B4Rz/S6axbQXnW7HZjdQDYzM2tQXR81LGk8cCOVqZojq/sZ\nEmlCmZlZNur9PPnTgHnA05IA5gBPSDofeAWYW7XtnGTdUVpbW48sl0olSqVSnTHMzIqtXC5TLpcb\n3o8i+j/YljQPuC8i3t/HfS8C50TEW8mJ1+9RmYefDfwYeF/0egBJvVeZWUG9+OKL/Oymm1g2Z07d\nY2/fvp0P3ngjp512Wg7Jhh9JRER/Myd9GugSyjuBh4GFkrZL+nSvTY60dUS0AWuANuCfgGvd5mZm\nzdXvdE1EXD3A/af2un0zcHMGuczMLAO+jt3MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAu\neTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKzCVvZlZgLnkzswJzyZuZFZhL3sys\nwFzyZmYF5pI3Myuwgb7Ie5WkDkmbqtb9d0mbJT0t6R5Jk6vuWyFpq6TnJF2aZ3AzMxvYQEfytwGX\n9Vq3Djg9Is4EtgArACQtBq4EFidjviXJrxTMzJqo3xKOiIeAXb3WrY+InuTmo8CcZPly4M6I6IqI\nbcDzwHnZxjUzs3o0eqT9GeCBZHkW0F51Xzswu8H9m5lZA0anHSjpz4FDEfG9fjaLvla2trYeWS6V\nSpRKpbQxzMwKqVwuUy6XG95PqpKX9AfAEuAjVatfAeZW3Z6TrDtKdcmbmdnReh8Ar1y5MtV+6p6u\nkXQZ8AXg8og4UHXXvcBVksZKmg8sADakSmVmZpno90he0p3ARcCJkrYDX6FyNc1YYL0kgJ9HxLUR\n0SZpDdAGdAPXRkSf0zVmZjY4+i35iLi6j9Wr+tn+ZuDmRkOZmVk2Up94NbOR47ENG9ize3fd43a/\n/Tb4BX1TueTNbEBPrVvHKW1tTDr++LrGnQQsOOGEfEJZTVzyZlaTs2bOZPakSc2OYXXyxw6YmRWY\nS97MrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmAueTOzAvObocxsyHrhhRd455136h43btw4\nFixYkEOi4cclb2ZD0qnHHcdr3/8+e+sct+/QIfYtWsSCG2/MJddw45I3syHpQ7PTfXvoq3v3cv/h\nwxmnGb48J29mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXWb8lLWiWpQ9KmqnXTJK2XtEXSOklTqu5b\nIWmrpOckXZpncDMzG9hAR/K3AZf1WrccWB8RC4EHk9tIWgxcCSxOxnxLkl8pmJk1Ub8lHBEPAbt6\nrV4KrE6WVwNXJMuXA3dGRFdEbAOeB87LLqqZmdUrzZH2jIjoSJY7gBnJ8iygvWq7diDduxnMzCwT\nDb3jNSJCUvS3SV8rW1tbjyyXSiVKpVIjMczMCqdcLlMulxveT5qS75A0MyJ2SjoFeC1Z/wowt2q7\nOcm6o1SXvJmZHa33AfDKlStT7SfNdM29wLJkeRmwtmr9VZLGSpoPLAA2pEplZmaZ6PdIXtKdwEXA\niZK2A18GvgaskXQNsA34GEBEtElaA7QB3cC1EdHfVI6ZmeWs35KPiKuPcdfFx9j+ZuDmRkOZmVk2\nfB27mVmB+fPkB1F3dzc7duxIPX7q1Km0tLRkmMjMis4lP4j27t3LV2/9KjGx/lMVBzsPcsPVN3DO\nOefkkMzMisolP8h6ju9hbmnuwBv28vKGl3NIY2ZF5zl5M7MCc8mbmRWYS97MrMBc8mZmBeaSNzMr\nMJe8mVmBueTNzArMJW9mVmAueTOzAnPJm5kVmEvezKzAXPJmZgXmkjczKzCXvJlZgbnkzcwKLHXJ\nS1oh6VlJmyR9T9LxkqZJWi9pi6R1kqZkGdbMzOqTquQlzQP+EDg7It4PjAKuApYD6yNiIfBgctvM\nzJok7ZH8HqALGC9pNDAeeBVYCqxOtlkNXNFwQjMzSy1VyUfEW8A3gJeplPvuiFgPzIiIjmSzDmBG\nJinNzCyVVN/xKuk04AZgHvA2cJekT1RvExEhqc9vrG5tbT2yXCqVKJVKaWKYmRVWuVymXC43vJ+0\nX+T9AeDhiHgTQNI9wK8BOyXNjIidkk4BXutrcHXJm5nZ0XofAK9cuTLVftLOyT8HXCBpnCQBFwNt\nwH3AsmSbZcDalPs3M7MMpDqSj4inJd0OPA70AE8C3wYmAmskXQNsAz6WUU4zM0sh7XQNEfF14Ou9\nVr9F5ajezMyGAL/j1cyswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYC55M7MCc8mbmRWYS97M\nrMBc8mZmBeaSNzMrMJe8mVmBueTNzArMJW9mVmCpP2rYBtfhOMzadWtZ92/r6h47feJ0/ujTf5RD\nKjMb6lzyw8TMRTPp3N9JJ511jevu6qa9rT2nVGY21Lnkh4lxk8YxbtK4usd1HehiN7tzSGRmw4Hn\n5M3MCswlb2ZWYKlLXtIUSXdL2iypTdL5kqZJWi9pi6R1kqZkGdbMzOrTyJH83wAPRMQi4AzgOWA5\nsD4iFgIPJrfNzKxJUpW8pMnAhyNiFUBEdEfE28BSYHWy2WrgikxSmplZKmmP5OcDr0u6TdKTkr4j\naQIwIyI6km06gBmZpDQzs1TSXkI5Gjgb+FxEPCbpVnpNzURESIq+Bre2th5ZLpVKlEqllDHMzIqp\nXC5TLpcb3k/akm8H2iPiseT23cAKYKekmRGxU9IpwGt9Da4ueTMzO1rvA+CVK1em2k+q6ZqI2Als\nl7QwWXUx8CxwH7AsWbcMWJsqlZmZZaKRd7xeB9whaSzwS+DTwChgjaRrgG3AxxpOaCPKnj172LNn\nT+rxs2bN4rjj/PYPs3elLvmIeBo4t4+7Lk4fx0a6jRs38u1vP84JJ0yqe+yYMa/yzW/e6JI3q+LP\nrknh/n+6n8c2PTbwhr10H+5mf9f+HBIVy+jRZzF37m/WPe6ll27KIY3Z8OaST6HjjQ62j9vOtNnT\n6h47Z9ScHBKZmfXNJZ/SCS0n0DKtpdkxzMz65clLM7MC85G8FUZXF9xxxz8i1X/scsYZp3HWWWfm\nkMqsuVzyVhgTJy7lkUd66h63a9cvmTz5VZe8FZJL3gpjxoz3pxrX3X0AeCvbMGZDhOfkzcwKzCVv\nZlZgnq4xGyEOHz5Md3d3qrE90ecHytow4JI3GyG2bNnCXd/4BmMOH65/8P79aPbs7ENZ7lzyZiPI\nwgiueu97mx3DBpHn5M3MCswlb2ZWYC55M7MCc8mbmRWYS97MrMBc8mZmBeZLKC1zPT09rF69hoMH\n6x/71ltvAIszz2Q2UjVU8pJGAY8D7RHxnyRNA/4BeC/JF3lHxO6GU9qw89Of/jsTJqT7HveTTz4x\n4zQ20hw4eJCtW7emGjtz5kwmTpyYcaLmafRI/nqgDXj3N7IcWB8RX5f0peT28gYfw4YlcdJJi5od\nwkag40eNYvrLL7PhllvqHrvj0CGWLF/O4sXFeTWZuuQlzQGWADcBf5KsXgpclCyvBsq45M1sEE0f\nP56Pn3pqqrFrtm/POE3zNXLi9a+BLwDV39IwIyI6kuUOYEYD+zczswalOpKX9NvAaxGxUVKpr20i\nIiT1+dF1ra2tR5ZLpRKlUp+7MDMbscrlMuVyueH9pJ2uuRBYKmkJcAIwSdLfAx2SZkbETkmnAK/1\nNbi65M3M7Gi9D4BXrlyZaj+ppmsi4saImBsR84GrgJ9ExCeBe4FlyWbLgLWpUpmZWSayuk7+3WmZ\nrwFrJF1DcgllRvu3BvQc7mHz5s2pxp588slMnz4940RmNlgaLvmI+Cnw02T5LeDiRvdp2dFxYu+E\nvfzVPX9V99h9u/dx3X++jgsvvDCHZGY2GPyO14IbPXY07/v196Ua+9ITL2WcxswGmz+7xsyswFzy\nZmYF5pI3Myswl7yZWYG55M3MCswlb2ZWYCP2EsqOjg5u+eYt7O/aX/fYfQf3MeHMCTmkMjPL1ogt\n+YhgN7uZ9ZFZdY+dznSOG+UXQWY29I3YkofKu0FHjx3RvwIzKzgfjpqZFZhL3syswFzyZmYF5pI3\nMyswn3U0A9au3cAPf/hE3eNaWuCrX72eiRMn5pDKrHEueRvxZs8+l4hzUo1tb7814zRm2XLJ24gn\nHYfkmUsrJj+zzcwKLFXJS5or6V8kPSvpF5I+n6yfJmm9pC2S1kmakm1cMzOrR9oj+S7gv0XE6cAF\nwH+VtAhYDqyPiIXAg8ltMzNrklQlHxE7I+KpZLkT2AzMBpYCq5PNVgNXZBHSzMzSaXhOXtI84Czg\nUWBGRHQkd3UAMxrdv5mZpddQyUtqAX4AXB8Re6vvi4gAopH9m5lZY1JfQilpDJWC//uIWJus7pA0\nMyJ2SjoFeK2vsa2trUeWS6USpVIpbQwzs0Iql8uUy+WG95Oq5CUJ+C7QFhHV7wa5F1gG3JL8d20f\nw/9DyZuZ2dF6HwCvXLky1X7SHsl/EPgE8Iykjcm6FcDXgDWSrgG2AR9LuX8zM8tAqpKPiH/l2PP5\nF6ePY2ZmWfI7Xs3MCswlb2ZWYP6AMjumt3bs4kc/+ldefLHPi6SOKSLo7s4plFmeenrYUC6zZePG\ngbftZcLUqVyyZEkOoRrjkrdjemfXAR56ciIdHfPqHtvSMj/7QGY5O3fSJN5+5pm6x3UeOsSTJ53k\nkrfhZ/z4mcyc+avNjmE2KOZPnZpq3Jv79vFkxlmy4jl5M7MCc8mbmRWYS97MrMBc8mZmBeYTr2bD\nzI/uv5/ON96oe9zbe/Ywvqcnh0Q2lLnkrV8HDx5k165ddY+TxJQp/vbHPGx9+GHO7uhg4tixdY+d\nNHlyDolsKHPJ2zGNPmE0L7zxEm88/bP6x3aN5qOXfDSHVAawcPp0Thw/vtkxbBhwydsxTZsznY5T\nJjD55PqO/iKCPVv25JRqaHnnHfiLv/guUv2nty66aBFLllySQyqz/88lb9aAefOuIc0XoL3+ehu7\ndtU/r25WL5d8wUUE3YfSfZBMz+HDGacpnnHj0r1DcuzYFsAlb/lzyRdc98Fu1v2vjXQfnFT32Ahx\n3KhROaQys8Ey7Ev+wZ88yPPbn6973MEDB+k+PDI+KrH74CQmz/gvzY5hZk0w7Ev+2V8+yxN7n6Bl\nWkt9A4+HKaf7Ej8zK7ZhX/IAk2dMZtrsac2OYWY25GT+sQaSLpP0nKStkr6U9f7NzKx2mZa8pFHA\n3wGXAYuBqyUtyvIx8rLtqW3NjnCUoZhp986nmh3hKNu2lZsd4ShDMRNAedu2Zkc4ijPlK+sj+fOA\n5yNiW0R0Ad8HLs/4MXIxFAt1KGZyyddmoEyHDx9m//79dHZ21v3TE/Vfl/+uoVhezpSvrOfkZwPb\nq263A+dn/Bhmw96uXbv459tvZ/fD6+sfvG8fes97sg9lhaRo4KjgqJ1JvwtcFhF/mNz+BHB+RFxX\ntU1k+ZifveGz3PPAPWiUGtpP5+udtJxU5xU6Ocsq06H9MxELMkgE+3Y/y/gppw+84QExhjGZPOZA\nOjufoaXljEF5rFoNlKmr63UmsIkJY7L7t1CL1zs7OallaD3Pi5ApIlh0xhn83/Up/qddI0lERN1F\nl3XJXwC0RsRlye0VQE9E3FK1zeA+q83MCmIolPxo4N+BjwCvAhuAqyNic2YPYmZmNct0Tj4iuiV9\nDvgRMAr4rgvezKx5Mj2SNzOzoSX373iVNE3SeklbJK2T1OdnCUhaIelZSZskfU/S8UMg0xRJd0va\nLKktOeeQm1pzJduOkrRR0n3NziRprqR/Sf5+v5D0+ZyyDPhGO0l/m9z/tKSz8shRTyZJH0+yPCPp\n3yTlfoa41jckSjpXUrek3xkKmSSVkuf0LySV885USy5JJ0r6Z0lPJbn+IOc8qyR1SNrUzzb1Pccj\nItcf4OvAF5PlLwFf62ObecALwPHJ7X8AljUzU3LfauAzyfJoYHKzf1dV2/4JcAdwb7MzATOBX02W\nW6icl1mUcY5RwPPJc2UM8FTvxwCWAA8ky+cDj+T8u6kl06+9+7yh8ibBpmeq2u4nwP3A7zY7EzAF\neBaYk9w+Mc9MdeRqBf7y3UzAm8DoHDN9GDgL2HSM++t+jud+JA8spVKWJP+9oo9t9gBdwPjk5O14\n4JVmZpI0GfhwRKyCyvmGiHg7x0w15UqyzaHyx/7fQGPXjmaQKSJ2RsRTyXInsBmYlXGOWt5odyRr\nRDwKTJE0I+McdWWKiJ9XPW8eBebkmKemTInrgLuB13POU2um3wd+EBHtABExGB+2X0uuHcC7n9M9\nCXgzInL7+NqIeAjo70uV636OD0bJz4iIjmS5AzgqUES8BXwDeJnKVTm7I+LHzcwEzAdel3SbpCcl\nfUdS3l+qWUsugL8GvgD05JynnkwASJpH5Ujk0Yxz9PVGu9k1bJNnqdaSqdo1wAM55oEaMkmaTaXM\n/meyKu8Tc7X8nhYA05Jpv8clfTLnTLXm+g5wuqRXgaeB6wchV3/qfo5ncnWNpPVUXrL39ufVNyIi\n+rpOXtJpwA1UXja9Ddwl6eMRcUezMlH53ZwNfC4iHpN0K7Ac+HLaTFnkkvTbwGsRsVFSqZEsWWWq\n2k8LlaPD65Mj+izVWkS9X9nkWWA171vSbwCfAT6YXxygtky3AsuTv6fI/9VgLZnGUPn39hEqr+R/\nLumRiNja5Fw3Ak9FRCnpqfWSzoyIvTnmGkhdz/FMSj4ijvltxMlJhJkRsVPSKcBrfWz2AeDhiHgz\nGXMPcCGVOedmZWoH2iPiseT23VRKviEZ5LoQWCppCXACMEnS7RHxqSZmQtIY4AfA/4mItWmz9OMV\nYG7V7blU/kb9bTOHfKf9aslEcrL1O1TeDd7fS/HBynQO8P1Kv3Mi8FuSuiLi3iZm2g68ERH7gf2S\nfgacCeRZ8rXkuhC4CSAifinpReBXgMdzzNWfup/jgzFdcy+wLFleBvRVAM8BF0galxxZXAy0NTNT\nROwEtktamKy6mMqJoTzVkuvGiJgbEfOBq4CfNFLwWWRK/mbfBdoi4taccjwOLJA0T9JY4MokW++s\nn0oyXUBl2q+D/AyYSdJ7gHuAT0RE/V9hlkOmiDg1IuYnz6G7gc/mWPA1ZQL+EfiQKleNjadyUjHP\nDqg113NU/u2TzH3/CpWLRJql/ud4XmeJq84GTwN+DGwB1gFTkvWzgB9WbfdFKiW6icqJhTFDINOZ\nwGNU5uLuIf+ra2rKVbX9ReR/dc2AmYAPUTk/8BSwMfm5LIcsv0Xlyp3ngRXJuj8G/rhqm79L7n8a\nOHsQnt/9ZqJycvzNqt/LhmZn6rXtbcDvDIVMwJ9VdcDn885U49/vROC+5Pm0Cfj9nPPcSeW85CEq\nr24+0+hz3G+GMjMrsMGYrjEzsyZxyZuZFZhL3syswFzyZmYF5pI3Myswl7yZWYG55M3MCswlb2ZW\nYP8PeHPoNmTwMCsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b442a50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(99): gibbs_step_dp(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After 100 iterations, our assignment looks correct! We went back to 3 clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can sample the mixture weights, if we need them, using the \"Conditional Distribution of Mixture Weights\" derived [here](http://nbviewer.ipython.org/github/tdhopper/notes-on-dirichlet-processes/blob/master/2015-09-02-fitting-a-mixture-model.ipynb#Conditional-Distribution-of-Mixture-Weights)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.21330625,  0.29838101,  0.48831275])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss = state['suffstats']\n",
    "alpha = [ss[cid].N + state['alpha_'] / state['num_clusters_'] \n",
    "         for cid in state['cluster_ids_']]\n",
    "stats.dirichlet(alpha).rvs(size=1).flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also sample the cluster means using [the method we derived earlier](http://nbviewer.ipython.org/github/tdhopper/notes-on-dirichlet-processes/blob/master/2015-09-02-fitting-a-mixture-model.ipynb#Conditional-Distribution-of-Cluster-Means):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster_id: 1 mean -0.0176257860235\n",
      "cluster_id: 2 mean -0.400581819532\n",
      "cluster_id: 3 mean 0.600302879661\n"
     ]
    }
   ],
   "source": [
    "for cluster_id in state['cluster_ids_']:\n",
    "    cluster_var = state['cluster_variance_']\n",
    "    hp_mean = state['hyperparameters_']['mean']\n",
    "    hp_var = state['hyperparameters_']['variance']\n",
    "    ss = state['suffstats'][cluster_id]\n",
    "\n",
    "    numerator = hp_mean / hp_var + ss.theta * ss.N / cluster_var\n",
    "    denominator = (1.0 / hp_var + ss.N / cluster_var)\n",
    "    posterior_mu = numerator / denominator\n",
    "    posterior_var = 1.0 / denominator\n",
    "\n",
    "    mean = stats.norm(posterior_mu, np.sqrt(posterior_var)).rvs()\n",
    "    print \"cluster_id:\", cluster_id, \"mean\", mean\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much thanks to Erik Sudderth's excellent introduction to nonparametric Bayes in [Chapter 2 of his dissertation](http://cs.brown.edu/~sudderth/papers/sudderthPhD.pdf). Algorithms 2.2 and 2.3 in that piece are the clearest formulation of collapsed Gibbs sampling for mixture models that I have come across."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}