2015-09-02-fitting-a-mixture-model.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-09-02-fitting-a-mixture-model.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": [
    "# Fitting a Mixture Model with Gibbs Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import stats\n",
    "from collections import namedtuple, Counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we receive some data that looks like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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 0x109c0f290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = pd.Series.from_csv(\"clusters.csv\")\n",
    "_=data.hist(bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that these data exist in three separate clusters. We want to develop a method for finding these _latent_ clusters. One way to start developing a method is to attempt to describe the process that may have generated these data.\n",
    "\n",
    "For simplicity and sanity, let's assume that each data point is generated independently of the other. Moreover, we will assume that within each cluster, the data points are identically distributed. In this case, we will assume each cluster is normally distributed and that each cluster has the same variance, $\\sigma^2$.\n",
    "\n",
    "Given these assumptions, our data could have been generated by the following process. For each data point, randomly select 1 of 3 clusters from the distribution $\\text{Discrete}(\\pi_1, \\pi_2, \\pi_3)$. Each cluster $k$ corresponds to a parameter $\\theta_k$ for that cluster, sample a data point from $\\mathcal{N}(\\theta_k, \\sigma^2)$.\n",
    "\n",
    "Equivalently, we could consider these data to be generated from a probability distribution with this probability density function:\n",
    "\n",
    "$$\n",
    "p(x_i \\,|\\, \\pi, \\theta_1, \\theta_2, \\theta_3, \\sigma)=\n",
    "    \\sum_{k=1}^3 \\pi_k\\cdot\n",
    "        \\frac{1}{\\sigma\\sqrt{2\\pi}}\n",
    "        \\text{exp}\\left\\{\n",
    "            \\frac{-(x_i-\\theta_k)^2}{2\\sigma^2}\n",
    "         \\right\\}\n",
    "$$\n",
    "\n",
    "where $\\pi$ is a 3-dimensional vector giving the _mixing proportions_. In other words, $\\pi_k$ describes the proportion of points that occur in cluster $k$.\n",
    "\n",
    "\n",
    "That is, _the probability distribution describing $x$ is a linear combination of normal distributions_.\n",
    "\n",
    "We want to use this _generative_ model to formulate an algorithm for determining the particular parameters that generated the dataset above. The $\\pi$ vector is unknown to us, as is each cluster mean $\\theta_k$. \n",
    "\n",
    "We would also like to know $z_i\\in\\{1, 2, 3\\}$, the latent cluster for each point. It turns out that introducing $z_i$ into our model will help us solve for the other values.\n",
    "\n",
    "The joint distribution of our observed data (`data`) along with the assignment variables is given by:\n",
    "\n",
    "\\begin{align}\n",
    "p(\\mathbf{x}, \\mathbf{z} \\,|\\, \\pi, \\theta_1, \\theta_2, \\theta_3, \\sigma)&=\n",
    "       p(\\mathbf{z} \\,|\\, \\pi)\n",
    "       p(\\mathbf{x} \\,|\\, \\mathbf{z}, \\theta_1, \\theta_2, \\theta_3, \\sigma)\\\\\n",
    "       &= \\prod_{i=1}^N p(z_i \\,|\\, \\pi)\n",
    "           \\prod_{i=1}^N p(x_i \\,|\\, z_i, \\theta_1, \\theta_2, \\theta_3, \\sigma) \\\\\n",
    "       &= \\prod_{i=1}^N \\pi_{z_i}\n",
    "           \\prod_{i=1}^N \n",
    "           \\frac{1}{\\sigma\\sqrt{2\\pi}}\n",
    "        \\text{exp}\\left\\{\n",
    "            \\frac{-(x_i-\\theta_{z_i})^2}{2\\sigma^2}\n",
    "         \\right\\}\\\\\n",
    "       &= \\prod_{i=1}^N \n",
    "           \\left(\n",
    "           \\pi_{z_i}\n",
    "           \\frac{1}{\\sigma\\sqrt{2\\pi}}\n",
    "        \\text{exp}\\left\\{\n",
    "            \\frac{-(x_i-\\theta_{z_i})^2}{2\\sigma^2}\n",
    "         \\right\\}\n",
    "         \\right)\\\\\n",
    "       &=\n",
    "           \\prod_i^n\n",
    "           \\prod_k^K\n",
    "           \\left(\n",
    "               \\pi_k \n",
    "               \\frac{1}{\\sigma\\sqrt{2\\pi}}\n",
    "        \\text{exp}\\left\\{\n",
    "            \\frac{-(x_i-\\theta_k)^2}{2\\sigma^2}\n",
    "         \\right\\}\n",
    "           \\right)^{\\delta(z_i, k)}\n",
    "\\end{align}\n",
    "\n",
    "### Keeping Everything Straight\n",
    "\n",
    "Before moving on, we need to devise a way to keep all our data and parameters straight. Following ideas suggested by [Keith Bonawitz](http://people.csail.mit.edu/bonawitz/Composable%20Probabilistic%20Inference%20with%20Blaise%20-%20Keith%20Bonawitz%20PhD%20Thesis.pdf), let's define a \"state\" object to store all of this data. \n",
    "\n",
    "It won't yet be clear why we are defining some components of `state`, however we will use each part eventually! As an attempt at clarity, I am using a trailing underscore in the names of members that are fixed. We will update the other parameters as we try to fit the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SuffStat = namedtuple('SuffStat', 'theta N')\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)\n",
    "\n",
    "def initial_state():\n",
    "    num_clusters = 3\n",
    "    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': [None, None, None],\n",
    "        'assignment': [random.choice(cluster_ids) for _ in data],\n",
    "        'pi': [alpha / num_clusters for _ in cluster_ids],\n",
    "        'cluster_means': [-1, 0, 1]\n",
    "    }\n",
    "    update_suffstats(state)\n",
    "    return state\n",
    "\n",
    "state = initial_state()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num_clusters_\n",
      "suffstats\n",
      "data_\n",
      "cluster_means\n",
      "cluster_variance_\n",
      "cluster_ids_\n",
      "assignment\n",
      "pi\n",
      "alpha_\n",
      "hyperparameters_\n"
     ]
    }
   ],
   "source": [
    "for k, v in state.items():\n",
    "    print k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs Sampling\n",
    "\n",
    "The [theory of Gibbs sampling](https://en.wikipedia.org/wiki/Gibbs_sampling) tells us that given some data $\\bf y$ and a probability distribution $p$ parameterized by $\\gamma_1, \\ldots, \\gamma_d$, we can successively draw samples from the distribution by sampling from\n",
    "\n",
    "$$\\gamma_j^{(t)}\\sim p(\\gamma_j \\,|\\, \\gamma_{\\neg j}^{(t-1)})$$\n",
    "    \n",
    "where $\\gamma_{\\neg j}^{(t-1)}$ is all current values of $\\gamma_i$ except for $\\gamma_j$. If we sample long enough, these $\\gamma_j$ values will be random samples from $p$. \n",
    "\n",
    "In deriving a Gibbs sampler, it is often helpful to observe that \n",
    "\n",
    "$$\n",
    "    p(\\gamma_j \\,|\\, \\gamma_{\\neg j})\n",
    "        = \\frac{\n",
    "            p(\\gamma_1,\\ldots,\\gamma_d)\n",
    "        }{\n",
    "            p(\\gamma_{\\neg j})\n",
    "        } \\propto p(\\gamma_1,\\ldots,\\gamma_d).\n",
    "$$\n",
    "\n",
    "The conditional distribution is proportional to the joint distribution. We will get a lot of mileage from this simple observation by dropping constant terms from the joint distribution (relative to the parameters we are conditioned on).\n",
    "\n",
    "The $\\gamma$ values in our model are each of the $\\theta_k$ values, the $z_i$ values, and the $\\pi_k$ values. Thus, we need to derive the conditional distributions for each of these.\n",
    "\n",
    "Many derivation of Gibbs samplers that I have seen rely on a lot of handwaving and casual appeals to conjugacy. I have tried to add more mathematical details here. I would gladly accept feedback on how to more clearly present the derivations! I have also tried to make the derivations more concrete by immediately providing code to do the computations in this specific case. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conditional Distribution of Assignment\n",
    "\n",
    "For berevity, we will use\n",
    "\n",
    "$$\n",
    "p(z_i=k \\,|\\, \\cdot)=\n",
    "p(z_i=k \\,|\\, \n",
    "        z_{\\neg i}, \\pi,\n",
    "        \\theta_1, \\theta_2, \\theta_3, \\sigma, \\bf x\n",
    "        ).\n",
    "        $$\n",
    "        \n",
    "Because cluster assignements are conditionally independent given the cluster weights and paramters,\n",
    "\n",
    "\\begin{align}\n",
    "    p(z_i=k \\,|\\, \\cdot) \n",
    "        &\\propto\n",
    "           \\prod_i^n\n",
    "           \\prod_k^K\n",
    "           \\left(\n",
    "               \\pi_k \n",
    "               \\frac{1}{\\sigma\\sqrt{2\\pi}}\n",
    "        \\text{exp}\\left\\{\n",
    "            \\frac{-(x_i-\\theta_k)^2}{2\\sigma^2}\n",
    "         \\right\\}\n",
    "           \\right)^{\\delta(z_i, k)} \\\\\n",
    "        &\\propto\n",
    "            \\pi_k \\cdot\n",
    "            \\frac{1}{\\sigma\\sqrt{2\\pi}} \n",
    "            \\text{exp}\\left\\{\n",
    "                \\frac{-(x_i-\\theta_k)^2}{2\\sigma^2}\n",
    "             \\right\\}\n",
    "\\end{align}\n",
    "\n",
    "This equation intuitively makes sense: point $i$ is more likely to be in cluster $k$ if $k$ is itself probable ($\\pi_k\\gg 0$) and $x_i$ is close to the mean of the cluster $\\theta_k$.\n",
    "\n",
    "For each data point $i$, we can compute $p(z_i=k \\,|\\, \\cdot)$ for each of cluster $k$. These values are the unnormalized parameters to a discrete distribution from which we can sample assignments.\n",
    "\n",
    "Below, we define functions for doing this sampling. `sample_assignment` will generate a sample from the posterior assignment distribution for the specified data point. `update_assignment` will sample from the posterior assignment for each data point and update the `state` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def log_assignment_score(data_id, cluster_id, state):\n",
    "    \"\"\"log p(z_i=k \\,|\\, \\cdot) \n",
    "    \n",
    "    We compute these scores in log space for numerical stability.\n",
    "    \"\"\"\n",
    "    x = state['data_'][data_id]\n",
    "    theta = state['cluster_means'][cluster_id]\n",
    "    var = state['cluster_variance_']\n",
    "    log_pi = np.log(state['pi'][cluster_id])\n",
    "    return log_pi + stats.norm.logpdf(x, theta, var)\n",
    "\n",
    "\n",
    "def assigment_probs(data_id, state):\n",
    "    \"\"\"p(z_i=cid \\,|\\, \\cdot) for cid in cluster_ids\n",
    "    \"\"\"\n",
    "    scores = [log_assignment_score(data_id, cid, state) for cid in state['cluster_ids_']]\n",
    "    scores = np.exp(np.array(scores))\n",
    "    return scores / scores.sum()\n",
    "\n",
    "\n",
    "def sample_assignment(data_id, state):\n",
    "    \"\"\"Sample cluster assignment for data_id given current state\n",
    "    \n",
    "    cf Step 1 of Algorithm 2.1 in Sudderth 2006\n",
    "    \"\"\"\n",
    "    p = assigment_probs(data_id, state)\n",
    "    return np.random.choice(state['cluster_ids_'], p=p)\n",
    "\n",
    "\n",
    "def update_assignment(state):\n",
    "    \"\"\"Update cluster assignment for each data point given current state \n",
    "    \n",
    "    cf Step 1 of Algorithm 2.1 in Sudderth 2006\n",
    "    \"\"\"\n",
    "    for data_id, x in enumerate(state['data_']):\n",
    "        state['assignment'][data_id] = sample_assignment(data_id, state)\n",
    "    update_suffstats(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conditional Distribution of Mixture Weights\n",
    "\n",
    "We can similarly derive the conditional distributions of mixture weights by an application of Bayes theorem. Instead of updating each component of $\\pi$ separately, we update them together (this is called blocked Gibbs).\n",
    "\n",
    "\\begin{align}\n",
    "p(\\pi \\,|\\, \\cdot)&=\n",
    "p(\\pi \\,|\\, \n",
    "        \\bf{z}, \n",
    "        \\theta_1, \\theta_2, \\theta_3,\n",
    "        \\sigma, \\mathbf{x}, \\alpha\n",
    "        )\\\\\n",
    "&\\propto\n",
    "p(\\pi \\,|\\, \n",
    "        \\mathbf{x}, \n",
    "        \\theta_1, \\theta_2, \\theta_3,\n",
    "        \\sigma, \\alpha\n",
    "        )\n",
    "p(\\bf{z}\\ \\,|\\, \n",
    "        \\mathbf{x}, \n",
    "        \\theta_1, \\theta_2, \\theta_3,\n",
    "        \\sigma, \\pi, \\alpha\n",
    "        )\\\\\n",
    "&=\n",
    "p(\\pi \\,|\\, \n",
    "        \\alpha\n",
    "        )\n",
    "p(\\bf{z}\\ \\,|\\, \n",
    "        \\mathbf{x}, \n",
    "        \\theta_1, \\theta_2, \\theta_3,\n",
    "        \\sigma, \\pi, \\alpha\n",
    "        )\\\\\n",
    "&=\n",
    "\\prod_{i=1}^K \\pi_k^{\\alpha/K - 1}\n",
    "\\prod_{i=1}^K \\pi_k^{\\sum_{i=1}^N \\delta(z_i, k)} \\\\\n",
    "&=\\prod_{k=1}^3 \\pi_k^{\\alpha/K+\\sum_{i=1}^N \\delta(z_i, k)-1}\\\\\n",
    "&\\propto \\text{Dir}\\left(\n",
    "    \\sum_{i=1}^N \\delta(z_i, 1)+\\alpha/K, \n",
    "    \\sum_{i=1}^N \\delta(z_i, 2)+\\alpha/K,\n",
    "    \\sum_{i=1}^N \\delta(z_i, 3)+\\alpha/K\n",
    "    \\right)\n",
    "\\end{align}\n",
    "\n",
    "Here are Python functions to sample from the mixture weights given the current `state` and to update the mixture weights in the `state` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sample_mixture_weights(state):\n",
    "    \"\"\"Sample new mixture weights from current state according to \n",
    "    a Dirichlet distribution \n",
    "    \n",
    "    cf Step 2 of Algorithm 2.1 in Sudderth 2006\n",
    "    \"\"\"\n",
    "    ss = state['suffstats']\n",
    "    alpha = [ss[cid].N + state['alpha_'] / state['num_clusters_'] \n",
    "             for cid in state['cluster_ids_']]\n",
    "    return stats.dirichlet(alpha).rvs(size=1).flatten()\n",
    "\n",
    "def update_mixture_weights(state):\n",
    "    \"\"\"Update state with new mixture weights from current state\n",
    "    sampled according to a Dirichlet distribution \n",
    "    \n",
    "    cf Step 2 of Algorithm 2.1 in Sudderth 2006\n",
    "    \"\"\"\n",
    "    state['pi'] = sample_mixture_weights(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conditional Distribution of Cluster Means\n",
    "\n",
    "Finally, we need to compute the conditional distribution for the cluster means.\n",
    "\n",
    "We assume the unknown cluster means are distributed according to a normal distribution with hyperparameter mean $\\lambda_1$ and variance $\\lambda_2^2$. The final step in this derivation comes from the normal-normal conjugacy. For more information see [section 2.3 of this](http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf) and [section 6.2 this](http://fisher.osu.edu/~schroeder.9/AMIS900/ech6.pdf).)\n",
    "\n",
    "\\begin{align}\n",
    "p(\\theta_k \\,|\\, \\cdot)&=\n",
    "p(\\theta_k \\,|\\, \n",
    "        \\bf{z}, \\pi,\n",
    "        \\theta_{\\neg k},\n",
    "        \\sigma, \\bf x, \\lambda_1, \\lambda_2\n",
    "        ) \\\\\n",
    "&\\propto p(\\left\\{x_i \\,|\\, z_i=k\\right\\} \\,|\\, \\bf{z}, \\pi,\n",
    "        \\theta_1, \\theta_2, \\theta_3,\n",
    "        \\sigma, \\lambda_1, \\lambda_2) \\cdot\\\\\n",
    "    &\\phantom{==}p(\\theta_k \\,|\\, \\bf{z}, \\pi,\n",
    "        \\theta_{\\neg k},\n",
    "        \\sigma, \\lambda_1, \\lambda_2)\\\\\n",
    "&\\propto p(\\left\\{x_i \\,|\\, z_i=k\\right\\} \\,|\\, \\mathbf{z},\n",
    "        \\theta_k, \\sigma)\n",
    "        p(\\theta_k \\,|\\,  \\lambda_1, \\lambda_2)\\\\\n",
    "&= \\mathcal{N}(\\theta_k \\,|\\, \\mu_n, \\sigma_n)\\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "$$ \\sigma_n^2 = \\frac{1}{\n",
    "    \\frac{1}{\\lambda_2^2} + \\frac{N_k}{\\sigma^2}\n",
    "    } $$\n",
    "    \n",
    "and \n",
    "\n",
    "$$\\mu_n = \\sigma_n^2 \n",
    "    \\left(\n",
    "        \\frac{\\lambda_1}{\\lambda_2^2} + \n",
    "        \\frac{n\\bar{x_k}}{\\sigma^2}\n",
    "    \\right)\n",
    "$$\n",
    "\n",
    "Here is the code for sampling those means and for updating our state accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sample_cluster_mean(cluster_id, state):\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",
    "    return stats.norm(posterior_mu, np.sqrt(posterior_var)).rvs()\n",
    "\n",
    "\n",
    "def update_cluster_means(state):\n",
    "    state['cluster_means'] = [sample_cluster_mean(cid, state)\n",
    "                              for cid in state['cluster_ids_']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing each of these three updates in sequence makes a complete _Gibbs step_ for our mixture model. Here is a function to do that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gibbs_step(state):\n",
    "    update_assignment(state)\n",
    "    update_mixture_weights(state)\n",
    "    update_cluster_means(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially, we assigned each data point to a random cluster. We can see this by plotting a histogram of each cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHAFJREFUeJzt3X9w3Hd95/Hne39IsmRLsuLEcmyTGAZTSklCaNKQQBFg\njpBrA0NnQkuhvpK56dxMWw7uaBxuBuzpcYUUhrT01xy/GriUQgMNJuXABirKBUjsxPnh2E7sJHKs\nRD8syZL1a6Vd7fv+2G9cISxpP6v9eqWvXo8ZT1ar73v3FWn98ne/+/1h7o6IiKxsqVoHEBGRpVOZ\ni4gkgMpcRCQBVOYiIgmgMhcRSQCVuYhIApRV5mbWamb3mNlRMztiZr9mZm1mtt/MnjKzfWbWGndY\nERE5v3LXzP8C+I67vxK4AjgG7AL2u/t24AfR1yIiUgO22EFDZtYCHHL3l865/xjwRnfvM7N2oNPd\nfym+qCIiMp9y1sy3AafN7Etm9rCZfc7MmoCN7t4XLdMHbIwtpYiILKicMs8AVwN/4+5XA+PM2aTi\npdV7nRdARKRGMmUs0w10u/uB6Ot7gNuBXjNrd/deM9sE9M8dNDMVvIhIBdzdQpZfdM3c3XuBU2a2\nPbprB/AE8G1gZ3TfTuDeeeaX1Z+PfexjNc+gTMnKpUzKVO0/lShnzRzgj4C7zawOeBr4fSANfN3M\nbgW6gFsqSiAiIktWVpm7+6PANef51o7qxhERkUqsuiNAOzo6ah3hFyhT+ZZjLmUqjzLFa9H9zJf0\n4GYe5+OLiCSRmeHV/gBURESWP5W5iEgCqMxFRBJAZS4ikgDl7mcuIrKgL3zhHxgezlc0e911r+CG\nG66rcqLVRWUuIlVx8OBJ0ul3kU5ng+YGBp7k8ssHY0q1eqjMRaRqWlsvJ5OpD5qZmBgATscTaBXR\nNnMRkQRQmYuIJIDKXEQkAVTmIiIJoDIXEUkAlbmISAKozEVEEkBlLiKSACpzEZEEUJmLiCSAylxE\nJAFU5iIiCaAyFxFJAJW5iEgCqMxFRBJAZS4ikgAqcxGRBFCZi4gkgMpcRCQByroGqJl1AWeBGSDv\n7teaWRvwNeAyoAu4xd2HY8opIstcf38/U1PPk06HXQO0v7+f/n5dA3Spyr2gswMd7j40675dwH53\nv8PMbou+3lXtgCKyMjz5zJOk1j9EOlMXNDfc/wTbuyZjSrV6lFvmADbn65uBN0a37wI6UZmLrFru\nTnP7OjJ1DUFzE+N1uE/ElGr1KHebuQPfN7ODZvafo/s2untfdLsP2Fj1dCIiUpZy18xvcPceM7sY\n2G9mx2Z/093dzLz68UREpBxllbm790T/PW1m/wxcC/SZWbu795rZJqD/fLO7d+8+d7ujo4OOjo6l\nZhYRSZTOzk46OzuX9BjmvvAKtZk1Aml3HzWzJmAfsAfYAQy6+yfNbBfQ6u675sz6Yo8vIsnwqqve\nzvpf/kDwNvOe4wf5zRuG+dQd/zOmZCuPmeHucz+nXFA5a+YbgX82sxeXv9vd95nZQeDrZnYr0a6J\ngXlFRKRKFi1zd38WuOo89w9RWjsXEZEa0xGgIiIJoDIXEUkAlbmISAKozEVEEkBlLiKSACpzEZEE\nUJmLiCSAylxEJAFU5iIiCaAyFxFJAJW5iEgCqMxFRBJAZS4ikgAqcxGRBFCZi4gkgMpcRCQBVOYi\nIglQ1gWdRWR16O3tpaenp6LZydwE2bFx0tl80Fwul2N0bLSi55R/pzIXkXOOHz/OZz97mLVrNwXP\nnh18GXUv9JDJhNVKbmiI0xX+AyL/TmUuIj+nsXE7L3nJW4Ln6hueZ8OaNWTSYVtvz6a0tbca9FMU\nEUkAlbmISAKozEVEEkBlLiKSACpzEZEEUJmLiCSAylxEJAFU5iIiCaCDhkSk5vL5PN3d3RXNtre3\nBx91mkRl/QTMLA0cBLrd/TfNrA34GnAZ0AXc4u7DsaUUkcRKp+Hk6Cl2f3538GxqLMWnPvopmpub\nqx9shSn3n7MPAEeAddHXu4D97n6Hmd0Wfb0rhnwiknB19Rlatq9jS8eW4NlT3z0VQ6KVadFt5ma2\nBbgJ+Dxg0d03A3dFt+8C3hlLOhERKUs5a+afAT4MzH4fs9Hd+6LbfcDGagdb6Q4cOMD4+HhFs21t\nbVxxxRVVTiQiSbZgmZvZbwD97n7IzDrOt4y7u5n5fI+xe/fuc7c7Ojro6DjvwyTOff96H8dnjlO3\npi5obmJkgjdsfoPKXGQV6ezspLOzc0mPsdia+fXAzWZ2E9AANJvZV4A+M2t3914z2wT0z/cAs8t8\ntbn0FZeytm1t0Mxg9yDkYgokIsvS3BXdPXv2BD/GgtvM3f0j7r7V3bcBvw380N3fB+wFdkaL7QTu\nDX5mERGpmtCDhl7cnPIJ4K1m9hTw5uhrERGpkbL3tHf3HwE/im4PATviCiUiImF0OL+ISAKozEVE\nEkBlLiKSACpzEZEEUJmLiCSAzhspIucMDw9z5swUjY3zHgc4v3mPA1/c1HSe/v7w5xw6M0Q+n6/8\niRNEZS4i55w8eZK+R3upey4bPNtQLFb0nBlLkR48zemf/TR49uThcSYnJyt63qRRmYvIz2nJZHhV\nS0vw3LPD9aTMFl9wjvpMhq2Na/nl5nWLLzzHSavsZHZJpG3mIiIJoDIXEUkAlbmISAKozEVEEkBl\nLiKSANqbRUTO6erq4tjIJCN+Inh2wvM0Er5Lo1SHylxEzikUCuQbimQ2hb9pb2MNFeyZKFWiMheR\nn5eCTEZbYFca/cZERBJAZS4ikgAqcxGRBFCZi4gkgMpcRCQBVOYiIgmgMhcRSQCVuYhIAqjMRUQS\nQGUuIpIAKnMRkQRQmYuIJMCCZW5mDWb2gJk9YmZHzOzPovvbzGy/mT1lZvvMrPXCxBURkfNZsMzd\nPQe8yd2vAq4A3mRmrwd2AfvdfTvwg+hrERGpkUU3s7j7RHSzDkgDZ4Cbgbui++8C3hlLOhERKcui\n5zM3sxTwMPAy4G/d/Qkz2+jufdEifcDGGDOKSMLlcnkGBiYWX3COyYlp8vl8DIlWnkXL3N2LwFVm\n1gJ8z8zeNOf7bmY+3/zu3bvP3e7o6KCjo6PisCKSPJlMI889dwWnTs1bI/PqfWGKXC4XQ6oLq7Oz\nk87OziU9RtlXGnL3ETP7F+C1QJ+Ztbt7r5ltAvrnm5td5iIiczU2ttPY2F7R7KkXnq1ymtqYu6K7\nZ8+e4MdYbG+WDS/uqWJma4C3AoeAvcDOaLGdwL3BzywiIlWz2Jr5JuCuaLt5CviKu//AzA4BXzez\nW4Eu4JZ4Y648XSdOMTwwRsPahqC50cFR6prCZkREFixzd38cuPo89w8BO+IKlQTPPjPGQOFqGppa\ngubGhvpovPyFmFKJSFKVvc1cwjW2bGPt+rAdfQr5DKAyF5EwOpxfRCQBVOYiIgmgMhcRSQCVuYhI\nAugDUFlxxsbGmJ6ermg2lUrR2qqTfEryqMxlxbn33u+xb98zZLP1QXPF4gxbtsCf/ukHY0omUjsq\nc1lxCgVoarqRjRtfHTSXy40wMfHFmFKJ1Ja2mYuIJIDWzGMyPj7OBGMU02uC5ibGJxmtG40plYgk\nlco8JiMDA/hYH5mxsHMtT5/p50y+b/EFRURmUZnHaG22jvUNYSfNmk5nIfy0ziKyymmbuYhIAqjM\nRUQSQGUuIpIAKnMRkQRQmYuIJID2ZlmGuntOcvU1Lw+ey6Tr+YvPfI7Xve51MaQSkeVMZb7MtKyt\n5/KtbVxzbdjBRgCdP+mhWCzGkEpEljuV+TKUShmZdPgWMDOLIY2IrAQqc6mJ6elpuru7K5rt7+8n\nn99W5UQiK5vKXGrizJkzfPxvPg7N4bNHf9LNluZXVj+UyAqmMpfaaYCtv741eOypx3piCCOysmnX\nRBGRBFCZi4gkgMpcRCQBVOYiIgmgMhcRSYBF92Yxs63Al4FLKF024X+7+1+aWRvwNeAyoAu4xd2H\nY8x6wR0/fpyvfvX7zMyEz+YmMzStq34mEZHzKWfXxDzwQXd/xMzWAg+Z2X7g94H97n6Hmd0G7Ir+\nJEYul+PEiXVccslbg2ez2Sky6foYUomI/KJFy9zde4He6PaYmR0FNgM3A2+MFrsL6CRhZQ6QyTSw\ndu3GCuZasFTY9T9FRCoVdNCQmV0OvAZ4ANjo7i9eebgPCG88Eak6d2dsbKyi2UKhUOU08RsfH2d0\ndDR4LpPJsGZN+AntlquyyzzaxPIN4APuPjr7pE7u7mZ23ssQ7969+9ztjo4OOjo6Ks0qImW6/fZP\nMzER/qHN0aOTGE0xJIpHcWYDf/7n36Gp6UdBczMzed72tu285z3viilZmM7OTjo7O5f0GGWVuZll\nKRX5V9z93ujuPjNrd/deM9sE9J9vdnaZi8iFMToK27b9t+C53t4H6R/+WQyJ4rGm8a1s3ryT5uaw\nk/z09T1GPn8iplTh5q7o7tmzJ/gxFt010Uqr4F8Ajrj7nbO+tRfYGd3eCdw7d1ZERC6MctbMbwDe\nCzxmZoei+24HPgF83cxuJdo1MZaEEqSrq4vGxsbguaamJrZv3x5DIhG5EMrZm+X/Mf8a/I7qxpEl\naTO+dexb/ODUD4LGpnPTXLnhSj68/cMxBRORuOkUuAnStqWRS197Ka2trUFzI30jFPt1uTmRlUyH\n84uIJIDWzBfQ09ND1/GnGB14OHi2kK/8gKGREefhh8eD505N5dh45UTwmrmIrHwq8wWMjo4y/cLz\nXJwLv7LNRUAmFf7Gp66uhampd9BTwcV0To/dT2Fq5R30ISJLpzJfRF0qwyVN4QdRNA1nyVv4duhM\nZg2ZzJbgOYDUZLaiORFZ+bTNXEQkAbRmvoDu7m6eGu7nbOGp4NmcF2hGZ02cz7Fjx7jvK/eTuS/8\nXU9uoMhw/UMMPhN2/pHCzDiZtQ8GP5/ISqAyX8DMzAwz9UUat4Zvvmgky6zT18gc7k5xqp1t694T\nPHs2N8W1jVtpDzw46mwuz2OTE8HPJ7ISqMzLkEqpleNgGJl0XfBcOl0km8oEf8CcTqVAu9NLQqnM\nZUkGBgbI5XLBc4ODgzGkEVm9VOayJPu+8Q0G77+fNXVha9jPDw2R8fOeNVlEKqAyl6UpFHhbSwvb\nL7ooaOxwNsvddZVdQEFEfpF2TRQRSQCtmYskjLvT29uLWVfwbCWXX5PlQWUukkAnnjtOnz0SPDfm\nY2Tr9YZ9JVKZiyRUa3sFJ1wbBi/os4yVSP8Ei4gkgNbMZUkeevxJDj5ykg1NA0Fzw5MTFIs6MVgc\n3J2zIyP0nDwZPFuYmWFtDJkkfipzWZLBoSnGRn4d7LKgufHpaRrqumNKJYWpKdoqPKd+OvCYAVke\nVOayZA3ZFprrLw6aSdk02Yz2nIhTfUZ/vVcT/bYTZqY4Q6EQdoGKwkwheOZFxWKRmWKRQjHspCcz\ngcuLyMJU5glixSLHDxyg+8lHg+YmR6a4dOLlFT3niVNd9A5s4uh4+KH5hToVuki1qMwTpK2+nmvX\nNdPa0hA0d2pqhJGRyravFr1Iti1F6yVhzyki1aVdE0VEEmBVrJk/+eST5Cv4ZF+naU2emeIMhw8f\nrmj20ksvpa2trcqJRKpjVZT53//9/6Wn5xLS6bD9mru6WjHbFFMqudBSlmYy9yo+85mjwbMTEy/w\noQ9drzKXZWtVlHmhAO3tN9HQEHZ48+joEeobvh9TKrnQsukGLmp5M1u3vi149uTJ+2JIJFI92mYu\nIpIAi66Zm9kXgf8I9Lv7q6P72oCvAZcBXcAt7j4cY04pg3uKp5/Okc2G7TM+eHYKnxqJKZVIfNyL\n/PjfOkmn00FzExPPUiicYOfOd8WU7MIrZzPLl4DPAl+edd8uYL+732Fmt0Vf74ohnwSoq3stL7zw\nK8FzI2NDNGQfiiGRSLzWX5RlbPq54Lnhsyd59pnjMSSqnUXL3N1/bGaXz7n7ZuCN0e27gE5U5jXX\n2Nhe0Vwun9VV62VFuri1saK56ankneSt0m3mG929L7rdB2ysUh4REanAkvdmcXc3s2V9mfXHDh5k\nenQL2UzYyT3Pjo3BKrmC/GhulA9+9IPBc6fHhlbULlFTo6P85LvfDZ47M/ogT+14Bddcc00MqUSW\nrtK/h31m1u7uvVbaEbt/vgV379597nZHRwcdHR0VPmXl8mNjvNKMlsCzyJ3OZpksJv90oNm6FGs2\nNNBwXfgh+Y0PN8DzYR8+1Up9Os3rWiu4+g7w4FSu4pORiSyms7OTzs7OJT1GpWW+F9gJfDL6773z\nLTi7zGupLp0OPiVoNp0mZRZTouXDzEhlUtQ31gfPpjMpiivkZ2RmFZ8WNrUy/hdlhZq7ortnz57g\nxyhn18SvUvqwc4OZnQI+CnwC+LqZ3Uq0a2LwM19A04UCEzPTZFLTQXNTMzMxJZKlmpqZYWw67PcJ\nYECTLr4gCVTO3iy/M8+3dlQ5S2y6B/sYSp2iPjsUPrym+nlkidJwaKwHAq877EBzsZ43b90WSyyR\nWlpJn11VrEiRxkuyrG3UaVqToKWtst9joVAk36N3W5JMOpxfRCQBVOYiIgmgMhcRSQCVuYhIAqjM\nRUQSYFXszSKLG+mf5p5P7qtgLsf6GPKIxG10aIi7Pv3p8MG6Ot7xvvfRWuHRxHFRmQtNDevZ4L8F\np8Nn2wya1qjOZWXJplJsKzpveC789LnfGh5m+t3vjiHV0qjMhUy6jvVrL611DJELJmUp1tXX89L1\n4SsidWOBR6tdINpmLiKSACpzEZEEUJmLiCSAylxEJAH0AaiIrEq5Qp4nBwaC5/pGRpiamooh0dKo\nzEVk1UmnjWdTw9zRc3/w7BPPj/Pus2djSLU0KnMRWXXq6zNcvLWRra8Ouy4wwJGB8RgSLZ22mYuI\nJIDWzGVVyReLHB0MP9S1b3ycBx44xLp1FwfPNjRkuOmm/xA8d+zoUZ45ejR4bqZYDJ6RlU9lLqtG\nKmUU1zvHCb984ED6Ip54Yjv5/IaguZmZPA0NP66ozLufe46Ru+/mZW1tQXPFYpGmFXKRbakelbms\nGqmU0dxcX9Hs8PgG1q//FTZvviporlDIMTT044qeE2BrSwvXbt4cNFMsFqlPH6r4OWVl0jZzEZEE\n0Jq5yDJ17969/PB7D7Gu/nDYoDu5YmXvQFaTU6cynD6dq2CuyPDwcAyJlkZlLrJMjU9MciZ9PWs3\n/FLw7Obwz2lXlcbGdmZm3ot7+Oz01D8xMzNT/VBLpDIXWcbS6ToaG1tqHSNxUqksqVS2wunl+eHy\niinz4eFhBgcHK5r1Sv75FZkjl8sxNBS2J0yhkOPsMjxaUJJnxZT53r17ufOv91BXlw6enS6m0Z5a\nshSpNBw+coDDRw4EzRV9Gp/5HvC/4gkmElkxZQ5w0cYUN1wdfkWcfTM5UqkV9b8qy8yWS5rhkvC5\n6fwUTz9dqH4gkTm0a6KISAIsaXXVzG4E7gTSwOfd/ZNVSXUeZ8+eZeRMjuefD9/+WCymSemfLakV\nhwcffDB47OzZUWBd9fNIIlVc5maWBv4K2AE8Dxwws73uHn4yiTIMDg5y8ulXkx97afBssZiiubm0\nrX14+BitreG7esVJmcq3HHMtlqnBnYG/+7vgx2169lkytjWWTLWgTPFaypr5tcAJd+8CMLN/BN4B\nxFLmAOnMxbS0/OqSHmM5/vKUqXzLMddCmVJmrK3LctNLXhL8uIdOnuKno5W9pVxpP6daWY6ZKrWU\njQ+bgVOzvu6O7hMRkQtsKWvmF3zn7dxkF8dOfGtJjzF29hi5wtIeo9qUqXzLMddCmYozeeoKY7yy\ngs0sRYy+3CSjk09XNVOtJCXTzMzyPG7AKj2gxsyuA3a7+43R17cDxdkfgpqZjtYREamAuwcdHbOU\nMs8ATwJvAV4AHgR+J64PQEVEZH4Vb2Zx94KZ/SHwPUq7Jn5BRS4iUhsVr5mLiMjyUdVDacyszcz2\nm9lTZrbPzFrnWe52M3vCzB43s38ws9hOvhyQqdXM7jGzo2Z2JPpMoKaZomXTZnbIzL4dV55yM5nZ\nVjP71+h3d9jM/jimLDea2TEzO25mt82zzF9G33/UzF4TR47QXGb2u1Gex8zsfjO7otaZZi13jZkV\nzOxdyyGTmXVEr+vDZtZZ60xmtsHMvmtmj0SZ/lPMeb5oZn1m9vgCy4S9xt29an+AO4A/iW7fBnzi\nPMtcDjwD1Edffw3YWc0coZmi790FvD+6nQFaap0p+v6HgLuBvXHlCfjdtQNXRbfXUvrM5JVVzpEG\nTkSvkyzwyNznAG4CvhPd/jXgZ3H+bAJyve7F1w1wY9y5ysk0a7kfAvcBv1XrTEAr8ASwJfp6wzLI\ntBv4sxfzAINAJsZMbwBeAzw+z/eDX+PVPsj9ZkqlSPTfd55nmbNAHmiMPkRtpHQEaVwWzWRmLcAb\n3P2LUPo8wN1HapkpyrWF0i/188R/EuVFM7l7r7s/Et0eo3SAWPiZzxZ27mA0d88DLx6Mdt6s7v4A\n0GpmG6ucIziXu/901uvmAWBLrTNF/gi4Bzgdc55yM70H+Ia7dwO4+8AyyNQDNEe3m4FBd4/tDGnu\n/mPgzAKLBL/Gq13mG929L7rdB/zCk7v7EPBp4DlKe8EMu/v3q5wjKBOwDThtZl8ys4fN7HNm1ljj\nTACfAT4MFGPMEpoJADO7nNKaxQNVzlHOwWjnWybu4gw9SO5W4DuxJiojk5ltplRcfxvdFfeHZOX8\nnF4OtEWb7A6a2fuWQabPAa8ysxeAR4EPxJxpMcGv8eC9WcxsP6W323P9j9lfuLufbz9zM3sZ8F8p\nveUZAf7JzH7X3e8OzVKtTJR+DlcDf+juB8zsTmAX8NFaZTKz3wD63f2QmXVUmqOamWY9zlpKa3of\niNbQq6ncspn7TiXukir78c3sTcD7gRviiwOUl+lOYFf0OzXif4dXTqYspb9vb6H0zvynZvYzdz9e\nw0wfAR5x946oo/ab2ZXuPhpTpnIEvcaDy9zd3zrvM5c26Le7e6+ZbQL6z7PYrwI/cffBaOabwPWU\ntgtXpAqZuoFud3/xygP3UCrzilUh0/XAzWZ2E9AANJvZl93992qYCTPLAt8A/o+731tplgU8D8w+\nu9RWSr+fhZbZQryb6srNRfSh5+eAG919obfRFyrTa4F/LPU4G4C3m1ne3ffWMNMpYMDdJ4FJM/s3\n4EogrjIvJ9P1wMcB3P1pM3sWeAVwMKZMiwl+jVd7M8teYGd0eydwvr/sx4DrzGxNtKawAzhS5RxB\nmdy9FzhlZtuju3ZQ+oCmlpk+4u5b3X0b8NvAD5dS5NXIFP2+vgAccfc7Y8pxEHi5mV1uZnXAu6Ns\nc7P+XpTpOkqb6vqI16K5zOwlwDeB97r7iZjzlJXJ3V/q7tui19E9wH+JscjLygR8C3i9lfbUaqT0\nAV+cHVBOpmOU/t4TbZt+BaUdNWol/DVe5U9o24DvA08B+4DW6P5LgX+ZtdyfUCrLxylt5M9WM0eF\nma4EDlDaXvZN4t2bpaxMs5Z/I/HvzbJoJuD1lLbfPwIciv7cGEOWt1PaU+YEcHt03x8AfzBrmb+K\nvv8ocHWcP5tyc1H6oHpw1s/mwVpnmrPsl4B3LYdMwH+f1QF/XOtMlN61fDt6PT0OvCfmPF+l9Jnh\nNKV3Ku9f6mtcBw2JiCSArr8jIpIAKnMRkQRQmYuIJIDKXEQkAVTmIiIJoDIXEUkAlbmISAKozEVE\nEuD/AwDCNhHQX4aoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1066b8650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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 )\n",
    "    \n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each time we run `gibbs_step`, our `state` is updated with newly sampled assignments. Look what happens to our histogram after 5 steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGGxJREFUeJzt3X2QHPV95/H3R8/Pz4qEHmwBJ8USZQhgHoIfNLGBUpSc\n4JKUgfhBsakkZWIbkjvbEqmyV06hGKd84c4pX8o+ixM5LCIwUQCTRDJ4HBIMEkYCmQVLMshoAa0k\n0AOLHtiH7/0xjW6yWu3O9EzvrHo/r6otenr61/Op3dGHnl9PzygiMDOzfBrS6ABmZpYdl7yZWY65\n5M3Mcswlb2aWYy55M7Mcc8mbmeVYryUvaY2kVknbu63/nKTnJf1M0u1l61dK2inpBUlXZxXazMwq\nM6yP++8Evgnc9c4KSb8BLAPOj4h2SdOT9YuA64BFwGzgh5IWRERXJsnNzKxPvR7JR8RjwMFuqz8D\n/GVEtCfb7E/WXwOsi4j2iNgN7AIurW9cMzOrRpo5+fnAhyQ9Iako6X3J+llAS9l2LZSO6M3MrEH6\nmq453ZjJEXG5pEuA9cA5p9nWn5lgZtZAaUq+BbgfICK2SOqSNA14BZhbtt2cZN1/IMnFb2aWQkSo\n2jFppms2AB8GkLQAGBERB4AHgOsljZB0NqVpnc2nCTrgfr7yla80PIMzOdNgzOVMlf2k1euRvKR1\nwGJgqqQ9wJeBNcCa5G2VbwOfTIq7WdJ6oBnoAG6KWpKZmVnNei35iLjhNHd94jTbrwZW1xrKzMzq\nw1e8JgqFQqMjnMKZKuNMlRuIuZwpW+rvGRVJnsUxM6uSJKKfTryamdkZwiVvZpZjLnkzsxxzyZuZ\n5ViaK17NzCry2muvsfHee6G9vfrBw4dz1e/9HrNmzap/sEHEJW9mmTl+/DjHnnqKqydPrnrspoMH\nOfZbv5VBqsHFJW9mmRo9YgTnpCj50W1tGaQZfDwnb2aWYy55M7Mcc8mbmeWYS97MLMdc8mZmOeaS\nNzPLMZe8mVmOueTNzHLMJW9mlmMueTOzHOu15CWtkdSafGl39/v+q6QuSVPK1q2UtFPSC5KuziKw\nmZlVrq8j+TuBJd1XSpoLXAX8smzdIuA6YFEy5luS/ErBzKyBei3hiHgMONjDXf8d+GK3ddcA6yKi\nPSJ2A7uAS+sR0szM0qn6SFvSNUBLRDzb7a5ZQEvZ7RZgdg3ZzMysRlV91LCkMcCtlKZqTq7uZUik\nCWVmZvVR7efJnwvMA56RBDAH+Kmky4BXgLll285J1p2iqanp5HKhUKBQKFQZw8ws34rFIsViseb9\nKKL3g21J84AHI+K9Pdz3EnBxRLyRnHj9HqV5+NnAD4H/FN0eQFL3VWaWUy+99BL/etttLJ8zp+qx\nd+3Zw/tvvZVzzz03g2RnHklERG8zJz3q6y2U64DHgQWS9kj6VLdNTrZ1RDQD64Fm4J+Am9zmZmaN\n1et0TUTc0Mf953S7vRpYXYdcZmZWB34fu5lZjrnkzcxyzCVvZpZjLnkzsxxzyZuZ5ZhL3swsx1zy\nZmY55pI3M8sxl7yZWY655M3Mcswlb2aWYy55M7Mcc8mbmeWYS97MLMdc8mZmOeaSNzPLMZe8mVmO\nueTNzHLMJW9mlmN9fZH3GkmtkraXrfsrSc9LekbS/ZImlt23UtJOSS9IujrL4GZm1re+juTvBJZ0\nW7cROC8iLgB2ACsBJC0CrgMWJWO+JcmvFMzMGqjXEo6Ix4CD3dZtioiu5OaTwJxk+RpgXUS0R8Ru\nYBdwaX3jmplZNWo90v408HCyPAtoKbuvBZhd4/7NzKwGw9IOlPTnwNsR8b1eNoueVjY1NZ1cLhQK\nFAqFtDHMzHKpWCxSLBZr3k+qkpf0B8BS4CNlq18B5pbdnpOsO0V5yZuZ2am6HwCvWrUq1X6qnq6R\ntAT4AnBNRBwvu+sB4HpJIySdDcwHNqdKZWZmddHrkbykdcBiYJqkPcBXKL2bZgSwSRLATyLipoho\nlrQeaAY6gJsiosfpGjMz6x+9lnxE3NDD6jW9bL8aWF1rKDMzq4/UJ17NbPDYsnkzRw4dqnrcocOH\nwS/oG8olb2Z92rZxI2c1NzNh5Miqxk0H5o8alU0oq4hL3swqcuHMmcyeMKHRMaxK/tgBM7Mcc8mb\nmeWYS97MLMdc8mZmOeaSNzPLMZe8mVmOueTNzHLMJW9mlmO+GMrMBqwXX3yRt956q+pxo0ePZv78\n+RkkOvO45M1sQDpnyBD23XMPb1Y57ujbb3N04ULm33prJrnONC55MxuQPjA73beHvvrmmzzU2Vnn\nNGcuz8mbmeWYS97MLMdc8mZmOeaSNzPLMZe8mVmO9VryktZIapW0vWzdFEmbJO2QtFHSpLL7Vkra\nKekFSVdnGdzMzPrW15H8ncCSbutWAJsiYgHwSHIbSYuA64BFyZhvSfIrBTOzBuq1hCPiMeBgt9XL\ngLXJ8lrg2mT5GmBdRLRHxG5gF3Bp/aKamVm10hxpz4iI1mS5FZiRLM8CWsq2awHSXc1gZmZ1UdMV\nrxERkqK3TXpa2dTUdHK5UChQKBRqiWFmljvFYpFisVjzftKUfKukmRGxV9JZwL5k/SvA3LLt5iTr\nTlFe8mZmdqruB8CrVq1KtZ800zUPAMuT5eXAhrL110saIelsYD6wOVUqMzOri16P5CWtAxYD0yTt\nAb4MfA1YL+lGYDfwUYCIaJa0HmgGOoCbIqK3qRwzM8tYryUfETec5q4rT7P9amB1raHMzKw+/D52\nM7Mc8+fJ96OOjg5ee+211OMnT57MuHHj6pjIzPLOJd+P3nzzTb56x1eJ8dWfqjjRdoJbbriFiy++\nOINkZpZXLvl+1jWyi7mFuX1v2M3Lm1/OII2Z5Z3n5M3Mcswlb2aWYy55M7Mcc8mbmeWYS97MLMdc\n8mZmOeaSNzPLMZe8mVmOueTNzHLMJW9mlmMueTOzHHPJm5nlmEvezCzHXPJmZjnmkjczy7HUJS9p\npaTnJG2X9D1JIyVNkbRJ0g5JGyVNqmdYMzOrTqqSlzQP+EPgooh4LzAUuB5YAWyKiAXAI8ltMzNr\nkLRH8keAdmCMpGHAGOBVYBmwNtlmLXBtzQnNzCy1VCUfEW8A3wBeplTuhyJiEzAjIlqTzVqBGXVJ\naWZmqaT6jldJ5wK3APOAw8C9kj5evk1EhKQev7G6qanp5HKhUKBQKKSJYWaWW8VikWKxWPN+0n6R\n9/uAxyPidQBJ9wO/DuyVNDMi9ko6C9jX0+Dykjczs1N1PwBetWpVqv2knZN/Abhc0mhJAq4EmoEH\ngeXJNsuBDSn3b2ZmdZDqSD4inpF0F/AU0AU8DXwbGA+sl3QjsBv4aJ1ymplZCmmna4iIrwNf77b6\nDUpH9WZmNgD4ilczsxxzyZuZ5ZhL3swsx1zyZmY55pI3M8sxl7yZWY655M3Mcswlb2aWYy55M7Mc\nc8mbmeWYS97MLMdc8mZmOeaSNzPLMZe8mVmOpf6oYetfndHJho0b2PjvG6seO3X8VP7oU3+UQSoz\nG+hc8meImQtn0nasjTbaqhrX0d5BS3NLRqnMbKBzyZ8hRk8YzegJo6se1368nUMcyiCRmZ0JPCdv\nZpZjLnkzsxxLXfKSJkm6T9LzkpolXSZpiqRNknZI2ihpUj3DmplZdWo5kv8fwMMRsRA4H3gBWAFs\niogFwCPJbTMza5BUJS9pIvDBiFgDEBEdEXEYWAasTTZbC1xbl5RmZpZK2iP5s4H9ku6U9LSk70ga\nC8yIiNZkm1ZgRl1SmplZKmnfQjkMuAj4bERskXQH3aZmIiIkRU+Dm5qaTi4XCgUKhULKGGZm+VQs\nFikWizXvJ23JtwAtEbEluX0fsBLYK2lmROyVdBawr6fB5SVvZman6n4AvGrVqlT7STVdExF7gT2S\nFiSrrgSeAx4ElifrlgMbUqUyM7O6qOWK188Bd0saAfwC+BQwFFgv6UZgN/DRmhPaoHLkyBGOHDmS\nevysWbMYMsSXf5i9I3XJR8QzwCU93HVl+jg22G3dupVvf/spRo2aUPXY4cNf5ZvfvNUlb1bGn12T\nwkP/9BBbtm/pe8NuOjo7ONZ+LINE+TJs2IXMnfvhqsf98pe3ZZDG7Mzmkk+h9UAre0bvYcrsKVWP\nnTN0TgaJzMx65pJPadS4UYybMq7RMczMeuXJSzOzHPORvOVGezvcffc/IlV/7HL++edy4YUXZJDK\nrLFc8pYb48cv44knuqoed/DgL5g48VWXvOWSS95yY8aM96Ya19FxHHijvmHMBgjPyZuZ5ZhL3sws\nxzxdYzZIdHZ20tHRkWpsV/T4gbJ2BnDJmw0SO3bs4N5vfIPhnZ3VDz52DM2eXf9QljmXvNkgsiCC\n69/97kbHsH7kOXkzsxxzyZuZ5ZhL3swsx1zyZmY55pI3M8sxl7yZWY75LZSDwFvH3uJv/8/fVj1O\niMWXL+Y973lPVeO6urpYu3Y9J05U/ZC88cYBYFH1A82sRzWVvKShwFNAS0T8Z0lTgL8H3k3yRd4R\ncajmlJba0OFDGftrY9nO9qrHHnjxAO9dkO5Dv378458zdmy673H/lV+Zlmqc2TuOnzjBzp07U42d\nOXMm48ePr3Oixqn1SP5moBl45zeyAtgUEV+X9KXk9ooaH8NqMGToEKa/e3qqsUcPHK3hkcX06Qtr\nGG+WzsihQ5n68stsvv32qse+9vbbLF2xgkWL8vNqMnXJS5oDLAVuA/4sWb0MWJwsrwWKuOTNrB9N\nHTOGj51zTqqx6/fsqXOaxqvlxOtfA18Ayr+lYUZEtCbLrcCMGvZvZmY1SnUkL+m3gX0RsVVSoadt\nIiIk9fjRdU1NTSeXC4UChUKPuzAzG7SKxSLFYrHm/aSdrrkCWCZpKTAKmCDp74BWSTMjYq+ks4B9\nPQ0uL3kzMztV9wPgVatWpdpPqumaiLg1IuZGxNnA9cCjEfEJ4AFgebLZcmBDqlRmZlYX9boY6p1p\nma8BV0naAXw4uW1mZg1S88VQEfFj4MfJ8hvAlbXu08zM6sMfa2BmlmMueTOzHHPJm5nlmEvezCzH\nXPJmZjnmkjczyzGXvJlZjg3aLw1pbW3l9m/ezrH2Y1WPPXriKGMvGJtBKjOz+hq0JR8RHOIQsz4y\nq+qxU5nKkKF+EWRmA9+gLXkADRHDRgzqX4GZ5ZwPR83Mcswlb2aWYy55M7Mcc8mbmeWYzzqaARs2\nbOYHP/hp1ePGjYOvfvVmxo8fn0Eqs9q55G3Qmz37EiIuTjW2peWOOqcxqy+XvA160hAkz1xaPvmZ\nbWaWY6lKXtJcST+S9Jykn0n6fLJ+iqRNknZI2ihpUn3jmplZNdIeybcDfxoR5wGXA38iaSGwAtgU\nEQuAR5LbZmbWIKlKPiL2RsS2ZLkNeB6YDSwD1iabrQWurUdIMzNLp+Y5eUnzgAuBJ4EZEdGa3NUK\nzKh1/2Zmll5NJS9pHPB94OaIeLP8vogIIGrZv5mZ1Sb1WyglDadU8H8XERuS1a2SZkbEXklnAft6\nGtvU1HRyuVAoUCgU0sYwM8ulYrFIsViseT+pSl6SgO8CzRFRfjXIA8By4Pbkvxt6GP4fSt7MzE7V\n/QB41apVqfaT9kj+/cDHgWclbU3WrQS+BqyXdCOwG/hoyv2bmVkdpCr5iPg3Tj+ff2X6OGZmVk/+\nWAM7regKOjs76ezsrGpcV1dXRonMrFoueTutV36+l9WP3sOMGT+qeuzhw/7EDLOBwCVvpxUBQ4de\nybx5/6XRUcz6R1cXm4tFdmzd2ve23YydPJmrli7NIFRtXPJmZolLJkzg8LPPVj2u7e23eXr6dJe8\nmdlAdvbkyanGvX70KE/XOUu9eOLUzCzHXPJmZjnmkjczyzGXvJlZjvnEq9kZ5l8eeoi2AweqHnf4\nyBHG+EK1Qcclb3aG2fn441zU2sr4ESOqHjth4sQMEtlA5pI3OwMtmDqVaWPGNDqGnQFc8nZaEcFL\nLS9x8OimqseOHDaSD13xoQxSDSxvvQV/8RffRar+9NbixQtZuvSqDFKZ/X8ueTutafOmsW/PEJhW\n5cCAg788mEmmgWbevBtJ8wVo+/c3c/Bg9fPqZtVyydtpDR0+jGHDhzF85PCqxpW++XFwGD063RWS\nI0aMA1zylj2/hdLMLMfO+CP5Rx59hF17dlU97sTxE3R0dmSQyMxs4DjjS/65XzzHT9/8KeOmjKtu\n4EiYdN6kbEKZmQ0QZ3zJA0ycMZEps6c0OoaZ2YBT9zl5SUskvSBpp6Qv1Xv/ZmZWubqWvKShwN8A\nS4BFwA2SFtbzMbKye9vuRkc4xUDMdGjvtkZHOMXu3cVGRzjFQMwEUNy9u9ERTuFM2ar3kfylwK6I\n2B0R7cA9wDV1foxMDMRCHYiZXPKV6StTZ2cnx44do62treqfrhreojoQy8uZslXvOfnZwJ6y2y3A\nZXV+DLMz3sGDB/nnu+7i0OPVX03M0aPoXe+qfyjLJdXzwhVJvwssiYg/TG5/HLgsIj5Xtk3U8zE/\nc8tnuP/h+9FQ1bSftv1tjJte5Tt0MtbwTDGEo4few/BRU0+uOnroOcZMOq/vscfFcKq7iCqttrZn\nGTfu/H55rEr1lam9fT9j2c7Y4f174dj+tjamjxtYz/M8ZIoIFp5/Pv+wKcX/tCskiYiouujqXfKX\nA00RsSS5vRLoiojby7YZPJdDmpnV0UAo+WHAz4GPAK8Cm4EbIuL5uj2ImZlVrK5z8hHRIemzwL8A\nQ4HvuuDNzBqnrkfyZmY2sGT+AWWSpkjaJGmHpI2SevwsAUkrJT0nabuk70kaOQAyTZJ0n6TnJTUn\n5xwyU2muZNuhkrZKerDRmSTNlfSj5O/3M0mfzyhLnxfaSfqfyf3PSLowixzVZJL0sSTLs5L+XVLm\nZ4grvSBR0iWSOiT9zkDIJKmQPKd/JqmYdaZKckmaJumfJW1Lcv1BxnnWSGqVtL2Xbap7jkdEpj/A\n14EvJstfAr7WwzbzgBeBkcntvweWNzJTct9a4NPJ8jBgYqN/V2Xb/hlwN/BAozMBM4FfS5bHUTov\ns7DOOYYCu5LnynBgW/fHAJYCDyfLlwFPZPy7qSTTr7/zvKF0kWDDM5Vt9yjwEPC7jc4ETAKeA+Yk\nt6dlmamKXE3AX76TCXgdGJZhpg8CFwLbT3N/1c/x/vio4WWUypLkv9f2sM0RoB0Yk5y8HQO80shM\nkiYCH4yINVA63xARhzPMVFGuJNscSn/s/w3U9t7ROmSKiL0RsS1ZbgOeB2bVOUclF9qdzBoRTwKT\nJM2oc46qMkXET8qeN08CczLMU1GmxOeA+4D9GeepNNPvA9+PiBaAiOiPD9uvJNdrwIRkeQLwekRk\n9vG1EfEY0Ns37lT9HO+Pkp8REa3JcitwSqCIeAP4BvAypXflHIqIHzYyE3A2sF/SnZKelvQdSVl/\nqWYluQD+GvgC0JVxnmoyASBpHqUjkSfrnKOnC+1mV7BNlqVaSaZyNwIPZ5gHKsgkaTalMvtfyaqs\nT8xV8nuaD0xJpv2ekvSJjDNVmus7wHmSXgWeAW7uh1y9qfo5Xpd310jaROkle3d/Xn4jIqKn98lL\nOhe4hdLLpsPAvZI+FhF3NyoTpd/NRcBnI2KLpDuAFcCX02aqRy5Jvw3si4itkgq1ZKlXprL9jKN0\ndHhzckRfT5UWUfdXNlkWWMX7lvQbwKeB92cXB6gs0x3AiuTvKbJ/NVhJpuGU/r19hNIr+Z9IeiIi\ndjY4163AtogoJD21SdIFEfFmhrn6UtVzvC4lHxGn/Tbi5CTCzIjYK+ksYF8Pm70PeDwiXk/G3A9c\nQWnOuVGZWoCWiNiS3L6PUsnXpA65rgCWSVoKjAImSLorIj7ZwExIGg58H/i/EbEhbZZevALMLbs9\nl9LfqLdt5pDttF8lmUhOtn6H0tXgWX/5bSWZLgbuKfU704DflNQeEQ80MNMe4EBEHAOOSfpX4AIg\ny5KvJNcVwG0AEfELSS8Bvwo8lWGu3lT9HO+P6ZoHgOXJ8nKgpwJ4Abhc0ujkyOJKoLmRmSJiL7BH\n0oJk1ZWUTgxlqZJct0bE3Ig4G7geeLSWgq9HpuRv9l2gOSLuyCjHU8B8SfMkjQCuS7J1z/rJJNPl\nlKb9WslOn5kkvQu4H/h4RFT/FWYZZIqIcyLi7OQ5dB/wmQwLvqJMwD8CH1DpXWNjKJ1UzLIDKs31\nAqV/+yRz379K6U0ijVL9czyrs8RlZ4OnAD8EdgAbgUnJ+lnAD8q2+yKlEt1O6cTC8AGQ6QJgC6W5\nuPvJ/t01FeUq234x2b+7ps9MwAconR/YBmxNfpZkkOU3Kb1zZxewMln3x8Afl23zN8n9zwAX9cPz\nu9dMlE6Ov172e9nc6Ezdtr0T+J2BkAn4b2Ud8PmsM1X495sGPJg8n7YDv59xnnWUzku+TenVzadr\nfY77Yigzsxzrj+kaMzNrEJe8mVmOueTNzHLMJW9mlmMueTOzHHPJm5nlmEvezCzHXPJmZjn2/wAa\nW9EfXzD5LQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b2ab710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(5):\n",
    "    gibbs_step(state)\n",
    "plot_clusters(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suddenly, we are seeing clusters that appear very similar to what we would intuitively expect: three Gaussian clusters.\n",
    "\n",
    "Another way to see the progress made by the Gibbs sampler is to plot the change in the model's log-likelihood after each step. The log likehlihood is given by:\n",
    "\n",
    "$$\n",
    "\\log p(\\mathbf{x} \\,|\\, \\pi, \\theta_1, \\theta_2, \\theta_3)\n",
    "\\propto \\sum_x \\log \\left(\n",
    "    \\sum_{k=1}^3 \\pi_k \\exp \n",
    "    \\left\\{ \n",
    "        -(x-\\theta_k)^2 / (2\\sigma^2)\n",
    "    \\right\\}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "We can define this as a function of our `state` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def log_likelihood(state):\n",
    "    \"\"\"Data log-likeliehood\n",
    "    \n",
    "    Equation 2.153 in Sudderth\n",
    "    \"\"\"\n",
    "    \n",
    "    ll = 0 \n",
    "    for x in state['data_']:\n",
    "        pi = state['pi']\n",
    "        mean = state['cluster_means']\n",
    "        sd = np.sqrt(state['cluster_variance_'])\n",
    "        ll += np.log(np.dot(pi, stats.norm(mean, sd).pdf(x)))\n",
    "    return ll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10b6898d0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFyZJREFUeJzt3X+Q3HWd5/HnW0JCdPmVXVeSEELYHW/Nui4QJewqOopQ\nscpKYCuY4C56e6ktl9zJ1ZZ3J4GqY6pur5TzBHFvYcsiHj+MgbgowhqDkWL0XC8GWAQ0IOGKsGSC\nCUYhREAyyfv+6O+E7tAz079merr7+aj6Vn/78/1293s6nXnN9/Pt7+cTmYkkSSPe0O4CJElTi8Eg\nSapgMEiSKhgMkqQKBoMkqYLBIEmqMKWCISKWRMTjEbE9Ij7d7nokqRfFVLmOISKOAn4GfBAYAu4H\nLs7Mx9pamCT1mKl0xHAW8GRm7sjMA8BtwLI21yRJPWcqBcNc4Jmy+zuLNknSJJpKwTA1+rQkqcdN\na3cBZYaAeWX351E6ajgsIgwPSWpAZkat+06lI4YHgL6IODUipgMrgLuO3CkzXVq0XHXVVW2voZsW\n30/fy6m61GvKHDFk5nBE/AfgHuAoYG36jSRJmnRTJhgAMvPbwLfbXcdEev55+NKX4JVXSvcjXn9b\nrW2s20b98z/D5z7X3HOU1zBWfWO1Zb5+OXSoenu1ZeS53vCG196/kfXRbo9sO3K92v3x2p54Au6+\nu1R7I0v5H3bNfB7K35dq67VsG2kba32s7c0u3/8+DA9T4cjP+3j3R9pq/T/WyP+v8f4YH+29rqet\nmkbei3pMqWDodj/4AfzFX8A558CCBbV/UMa6bcasWf3s2dP448f6ZVHrL5LM1/5THvkLu9Zl5HmO\nDJSR9dFuR9arbT+yfby2Q4fgV7/q50tfen1w1LqU/yzlP9NY72+12yPfm2rrtWwbaRtrfazt9S7l\n/+59ff0cc8zrPy+13h9pGy8Mx2ur5RfsePs0EkhjhVMt78VooV2rKXOBWy0iIjup3hHDw/C3fwv/\n8A9w443w4Q+3uyJJvSQiyDpOPnvEMMF27IA//3N405vgoYdg9ux2VyRJY5tK30rqOrfdBmedBX/2\nZ7Bpk6EgqTN4xDABXnwRLrsMfvjDUiCceWa7K5Kk2nnE0GL3318KgqOOggcfNBQkdR6PGFrk0KHS\nVz8//3n4+7+Hiy5qd0WS1BiDoQWGhuBjH4MDB+CBB+CUU9pdkSQ1zq6kJt15JyxaBO9/P9x3n6Eg\nqfN5xNCgl16CT30K7rkHvvEN+JM/aXdFktQaHjE04OGH4Z3vhH37StcmGAqSuonBUKcbb4QPfhCu\nuALWrYPjj293RZLUWg6JUYf9++Hkk2HrVnjrW9tWhiTVpd4hMTxiqMOdd8J73mMoSOpuBkMd1q0r\njY4qSd3MrqQa7d4Nf/AHpWsW3vjGtpQgSQ2xK2mC3HYbLF1qKEjqfgZDjb7yFbuRJPUGg6EGP/tZ\nqQvpAx9odyWSNPEMhhqsWwcrV5ZGTJWkbueQGOPILAXD177W7kokaXJ4xDCOLVtgxgw444x2VyJJ\nk8NgGMdXvlKaszlq/qKXJHU2r2MYw6uvwty5pSEwFiyYtJeVpJbyOoYWuuee0kVthoKkXmIwjMEh\nMCT1IruSRrFvH8ybB089BbNmTcpLStKEsCupRb7+9dJ0nYaCpF7TcDBExEUR8dOIOBgRZx6xbU1E\nbI+IxyPi/LL2RRHxaLHturL2GRFxe9G+JSLmN1pXqzgEhqRe1cwRw6PAhcD3yxsjYiGwAlgILAGu\njzj8Zc8bgFWZ2Qf0RcSSon0VsLdovxa4uom6mrZrF/zLv8CHP9zOKiSpPRoOhsx8PDOfqLJpGbA+\nMw9k5g7gSWBxRMwGjs3MrcV+twAXFOtLgZuL9TuAcxutqxXWr4cLL4RjjmlnFZLUHhNxjmEOsLPs\n/k5gbpX2oaKd4vYZgMwcBl6IiLb17tuNJKmXjTlWUkRsBk6qsumKzLx7Ykpqr5/+FH7xC3jf+9pd\niSS1x5jBkJnnNfCcQ8C8svsnUzpSGCrWj2wfecwpwK6ImAYcn5m/rPbkAwMDh9f7+/vp7+9voMTR\nrVsHF18Mb/D7WpI61ODgIIODgw0/vunrGCLiPuA/ZeaDxf2FwFeBsyh1EX0X+P3MzIj4EXAZsBX4\nFvDFzNwUEauBP8rMSyNiJXBBZq6s8loTeh3DoUOlq5zvvhve8Y4JexlJmlT1XsfQ8LDbEXEh8EXg\nd4BvRcRDmfmhzNwWERuAbcAwsLrst/lq4CZgJrAxMzcV7WuBWyNiO7AXeF0oTIYf/ACOP95QkNTb\nvPK5zCc+AaedBp/+9IS9hCRNunqPGAyGwm9+A3PmwI9/XBoKQ5K6hUNiNGjjxlIXkqEgqdcZDAWv\nXZCkEruSgOefh/nz4emn4YQTWv70ktRWdiU14B//Ec47z1CQJDAYALuRJKlcz3cl/eu/whlnlEZU\nnTGjpU8tSVOCXUl1Wr8eli83FCRpRE8HQybceqvdSJJUrqeD4ZFHYP9+ePe7212JJE0dPR0M69bB\nRz/qSKqSVK5nTz4fPFi6duE734GFC1vylJI0JXnyuUbf+x787u8aCpJ0pJ4NBq9dkKTqerIr6eWX\nYe5c+MlPSiOqSlI3syupBv/0T7BokaEgSdX0ZDDYjSRJo+u5rqS9e+H3fq80FMZxx7WoMEmawuxK\nGsfXvgZLlhgKkjSangsGu5EkaWw91ZX01FOweDEMDcHRR7ewMEmawuxKGsOGDaWRVA0FSRpdTwXD\ntm3wrne1uwpJmtp6KhiGhkoXtkmSRmcwSJIq9FwweLWzJI2tZ4LhxRdheBhOOKHdlUjS1NYzwTDS\njRQ1f2FLknpTzwWDJGlsDQdDRHwuIh6LiIcj4usRcXzZtjURsT0iHo+I88vaF0XEo8W268raZ0TE\n7UX7loiY3/iPVN2uXQaDJNWimSOG7wB/mJl/DDwBrAGIiIXACmAhsAS4PuJwB84NwKrM7AP6ImJJ\n0b4K2Fu0Xwtc3URdVXnEIEm1aTgYMnNzZh4q7v4IOLlYXwasz8wDmbkDeBJYHBGzgWMzc2ux3y3A\nBcX6UuDmYv0O4NxG6xqNwSBJtWnVOYZ/B2ws1ucAO8u27QTmVmkfKtopbp8ByMxh4IWImNWi2kov\nZjBIUk2mjbUxIjYDJ1XZdEVm3l3scyXwamZ+dQLqe52BgYHD6/39/fT399f0OINBUq8YHBxkcHCw\n4cc3NbpqRPxb4K+AczPzlaLtcoDM/GxxfxNwFfA0cF9mvq1ovxh4b2ZeWuwzkJlbImIa8GxmvrnK\n6zU8uurJJ8MPfwinnNLQwyWpY03a6KrFieP/DCwbCYXCXcDKiJgeEQuAPmBrZv4c2BcRi4uT0ZcA\n3yx7zMeL9eXAvY3WVc3Bg7BnD8ye3cpnlaTuNGZX0jj+DpgObC6+dPR/M3N1Zm6LiA3ANmAYWF32\nZ/5q4CZgJrAxMzcV7WuBWyNiO7AXWNlEXa+zezfMmuVw25JUi56YqOf+++Gv/xoefHACipKkKc6J\neqpw8DxJql3PBIPfSJKk2hgMkqQKBoMkqUJPBIMD6ElS7XoiGDxikKTaGQySpApdHwz798OBA07p\nKUm16vpgcEpPSapPzwSDJKk2BoMkqYLBIEmqYDBIkir0RDA4gJ4k1a4ngsEjBkmqncEgSarQ1RP1\nHDwIM2eWLnKbPn0CC5OkKcyJesrs2QMnnmgoSFI9ujoY7EaSpPoZDJKkCgaDJKmCwSBJqmAwSJIq\nGAySpAoGgySpgsEgSarQtcGwfz+8+qpTekpSvbo2GJzSU5Ia03AwRMR/i4iHI+LHEXFvRMwr27Ym\nIrZHxOMRcX5Z+6KIeLTYdl1Z+4yIuL1o3xIR8xv/kUrsRpKkxjRzxPA/MvOPM/N04E7gKoCIWAis\nABYCS4DrIw7/3X4DsCoz+4C+iFhStK8C9hbt1wJXN1EXALt2GQyS1IiGgyEzXyy7+1vAL4r1ZcD6\nzDyQmTuAJ4HFETEbODYztxb73QJcUKwvBW4u1u8Azm20rhEeMUhSY6Y18+CI+O/AJcDLwFlF8xxg\nS9luO4G5wIFifcRQ0U5x+wxAZg5HxAsRMSszf9lobUNDcNppjT5aknrXmMEQEZuBk6psuiIz787M\nK4ErI+Jy4AvAX05AjRUGBgYOr/f399Pf3191v6EhOOecia5GkqaewcFBBgcHG358SybqiYhTgI2Z\n+fYiJMjMzxbbNlE6//A0cF9mvq1ovxh4b2ZeWuwzkJlbImIa8GxmvrnK69Q8Uc/ZZ8M118Cf/mnT\nP54kdbRJm6gnIvrK7i4DHirW7wJWRsT0iFgA9AFbM/PnwL6IWFycjL4E+GbZYz5erC8H7m20rhGe\nY5CkxjRzjuEzEfFvgIPA/wMuBcjMbRGxAdgGDAOry/7MXw3cBMykdISxqWhfC9waEduBvcDKJuri\n4EHYvRtmz27mWSSpN3XlnM/PPgunn14KB0nqdc75jN1IktQMg0GSVKFrg2HOnHZXIUmdqWuDwSMG\nSWqMwSBJqtCVweAAepLUuK4MBo8YJKlxBoMkqULXBcOvf12a0vPEE9tdiSR1pq4LhpGvqjqlpyQ1\npiuDwW4kSWqcwSBJqmAwSJIqGAySpAoGgySpQlcGgwPoSVLjujIYPGKQpMZ11QxuBw/CzJmwfz9M\nnz6JhUnSFNbTM7g99xyccIKhIEnN6KpgsBtJkppnMEiSKhgMkqQKBoMkqYLBIEmqYDBIkioYDJKk\nCgaDJKlC1wTDr38Nv/mNU3pKUrOaDoaI+FREHIqIWWVtayJie0Q8HhHnl7UviohHi23XlbXPiIjb\ni/YtETG/3jqc0lOSWqOpYIiIecB5wNNlbQuBFcBCYAlwfcThX9c3AKsysw/oi4glRfsqYG/Rfi1w\ndb212I0kSa3R7BHDNcB/OaJtGbA+Mw9k5g7gSWBxRMwGjs3MrcV+twAXFOtLgZuL9TuAc+stxGCQ\npNZoOBgiYhmwMzMfOWLTHGBn2f2dwNwq7UNFO8XtMwCZOQy8UN41VYtduwwGSWqFaWNtjIjNwElV\nNl0JrAHOL9+9hXWNamBg4PB6f38//f39QOmI4dRTJ6MCSZraBgcHGRwcbPjxDc3HEBFvB+4FXiqa\nTqZ0BLAY+EuAzPxsse8m4CpK5yHuy8y3Fe0XA+/NzEuLfQYyc0tETAOezcw3V3ndUedjWL4cPvKR\n0iJJes2kzMeQmT/JzLdk5oLMXECpi+jMzNwN3AWsjIjpEbEA6AO2ZubPgX0Rsbg4GX0J8M3iKe8C\nPl6sL6cUOnXxHIMktcaYXUl1OPxnfGZui4gNwDZgGFhd9mf+auAmYCawMTM3Fe1rgVsjYjuwF1hZ\nbwEGgyS1RldM7XnoEBxzjFN6SlI1PTm15549TukpSa3SFcFgN5IktY7BIEmqYDBIkip0TTDMmdPu\nKiSpO3RNMHjEIEmtYTBIkip0RTA4gJ4ktU5XBINHDJLUOh0fDC+9BK+8ArPqGqRbkjSajg8Gp/SU\npNbqimCwG0mSWsdgkCRVMBgkSRUMBklSBYNBklTBYJAkVeiKYHAAPUlqnY6e2nNkSs8XX4QZM9pY\nmCRNYT01teeePXD88YaCJLVSRweDg+dJUut1dDB44lmSWs9gkCRVMBgkSRUMBklSBYNBklTBYJAk\nVTAYJEkVGg6GiBiIiJ0R8VCxfKhs25qI2B4Rj0fE+WXtiyLi0WLbdWXtMyLi9qJ9S0TMH+/1ndJT\nkiZGM0cMCVyTmWcUy7cBImIhsAJYCCwBro84PPHmDcCqzOwD+iJiSdG+CthbtF8LXD3eizulpyRN\njGa7kqr9Wl4GrM/MA5m5A3gSWBwRs4FjM3Nrsd8twAXF+lLg5mL9DuDc8V7YwfMkaWI0GwyfjIiH\nI2JtRJxQtM0BdpbtsxOYW6V9qGinuH0GIDOHgRciYsxOIs8vSNLEGDMYImJzcU7gyGUppW6hBcDp\nwLPA5yeh3sMcJ0mSJsa0sTZm5nm1PElE3AjcXdwdAuaVbT6Z0pHCULF+ZPvIY04BdkXENOD4zPxl\ntdcaGBgAYNMmOPPMfqC/lhIlqWcMDg4yODjY8OMbno8hImZn5rPF+t8A78rMjxYnn78KnEWpi+i7\nwO9nZkbEj4DLgK3At4AvZuamiFgN/FFmXhoRK4ELMnNlldc8PB/DRRfB8uWwYkVD5UtSz6h3PoYx\njxjGcXVEnE7p20lPAZ8AyMxtEbEB2AYMA6vLZtdZDdwEzAQ2Zuamon0tcGtEbAf2Aq8LhSN5jkGS\nJkbHzuA2fz4MDsKCBe2tSZKmunqPGDoyGA4dgpkzYd8+Z2+TpPH0xNSezz0Hxx1nKEjSROjIYPD8\ngiRNHINBklTBYJAkVTAYJEkVOjYYHEBPkiZGxwaDRwySNDE6MhgcQE+SJk5HBoNHDJI0cTouGF5+\nuTSt52//drsrkaTu1HHB4JSekjSxOjIY7EaSpIljMEiSKhgMkqQKBoMkqYLBIEmqYDBIkioYDJKk\nCh03tefRRyf79sExx7S7GknqDF0/tedxxxkKkjSROi4Y7EaSpIllMEiSKhgMkqQKBoMkqYLBIEmq\nYDBIkioYDJKkCk0FQ0R8MiIei4ifRMTVZe1rImJ7RDweEeeXtS+KiEeLbdeVtc+IiNuL9i0RMX+0\n1zQYJGliNRwMEfF+YCnwjsx8O/A/i/aFwApgIbAEuD7i8HxrNwCrMrMP6IuIJUX7KmBv0X4tcDWj\ncErP1hkcHGx3CV3F97N1fC/bq5kjhkuBz2TmAYDMfK5oXwasz8wDmbkDeBJYHBGzgWMzc2ux3y3A\nBcX6UuDmYv0O4NzRXtQpPVvH/3yt5fvZOr6X7dVMMPQB7y26fgYj4p1F+xxgZ9l+O4G5VdqHinaK\n22cAMnMYeCEiZjVRmySpQdPG2hgRm4GTqmy6snjsiZl5dkS8C9gAnNb6EiVJkyozG1qAbwPvK7v/\nJPA7wOXA5WXtm4DFlALmsbL2i4EbyvY5u1ifBjw3ymumi4uLi0v9Sz2/38c8YhjHncAHgO9FxFuB\n6Zn5i4i4C/hqRFxDqYuoD9iamRkR+yJiMbAVuAT4YvFcdwEfB7YAy4F7q71gPcPGSpIa00wwfBn4\nckQ8CrwKfAwgM7dFxAZgGzAMrM7XJn1YDdwEzAQ2Zuamon0tcGtEbAf2AiubqEuS1ISOmqhHkjTx\nOubK54hYUlwwtz0iPt3uejpdROyIiEci4qGI2Dr+IzQiIr4cEbuLo+WRtlkRsTkinoiI70TECe2s\nsZOM8n4ORMTO4vP5UNk1TxpDRMyLiPsi4qfFhceXFe11fT47Ihgi4ijgf1G6YG4hcHFEvK29VXW8\nBPoz84zMPKvdxXSY/03ps1jucmBzZr6V0jmyyye9qs5V7f1M4Jri83lGWbezxnYA+JvM/EPgbODf\nF78r6/p8dkQwAGcBT2bmjuKCutsoXUin5ngyvwGZ+X+AXx3RXH6R5s28dvGmxjHK+wl+PuuWmT/P\nzB8X6/uBxyh9Caiuz2enBMPhC+AKIxfNqXEJfDciHoiIv2p3MV3gLZm5u1jfDbylncV0iU9GxMMR\nsdauufpFxKnAGcCPqPPz2SnB4Bny1nt3Zp4BfIjS4eY57S6oWxTfwvMz25wbgAXA6cCzwOfbW05n\niYjfojS80H/MzBfLt9Xy+eyUYBgC5pXdn0fl8BqqU2Y+W9w+B3yDUnedGrc7Ik4CKMYF29Pmejpa\nZu7JAnAjfj5rFhFHUwqFWzPzzqK5rs9npwTDA5RGYz01IqZTGr31rjbX1LEi4o0RcWyx/ibgfODR\nsR+lcYxcpElxe+cY+2ocxS+vERfi57MmxUjWa4FtmfmFsk11fT475jqGiPgQ8AXgKGBtZn6mzSV1\nrIhYQOkoAUoXOa7z/axdRKwH3kdpCJjdwH8FvklpvLBTgB3ARzLz+XbV2EmqvJ9XAf2UupESeAr4\nRFkfuUYREe8Bvg88wmvdRWsojTZR8+ezY4JBkjQ5OqUrSZI0SQwGSVIFg0GSVMFgkCRVMBgkSRUM\nBklSBYNBklTBYJAkVfj/3sVunhlp/O8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ae490d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "state = initial_state()\n",
    "ll = [log_likelihood(state)]\n",
    "for _ in range(20):\n",
    "    gibbs_step(state)\n",
    "    ll.append(log_likelihood(state))\n",
    "pd.Series(ll).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See that the log likelihood improves with iterations of the Gibbs sampler. This is what we should expect: the Gibbs sampler finds state configurations that make the data we have seem \"likely\". However, the likelihood isn't strictly monotonic: it jitters up and down. Though it behaves similarly, the Gibbs sampler isn't optimizing the likelihood function. In its steady state, it is sampling from the posterior distribution. The `state` after each step of the Gibbs sampler is a sample from the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10b7287d0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGHBJREFUeJzt3X+UVWW9x/HPFwEVdV0lW/JTMaNA0wJMTCGOIYqaBqsI\nrdRcav4M9Ur+TsYy01rXlG6x7kpv/qgMfyQXQkiuMqCVkvwwLkggJsEEiD8RgYCZ7/3jOeMcxzPD\nzDnzzNl79vu11l6cs/c5ez9u93k++3n2s/eYuwsAkD2dKl0AAEBlEAAAkFEEAABkFAEAABlFAABA\nRhEAAJBRJQeAmY0zs2VmVmtmgxstu97MVpnZCjM7qWD+EDNbml92dzkFBwCUp5wWwFJJYyXNL5xp\nZodLGi/pcEmjJf3czCy/eIqk8929v6T+Zja6jO0DAMpQcgC4+wp3X1lk0ZckPeTuO939VUkvSxpq\nZj0l7efuC/Kfe0DSmFK3DwAoT4xrAL0krSt4v05S7yLza/LzAQAV0Lm5hWY2R1KPIotucPcZcYoE\nAGgPzQaAu48qYZ01kvoWvO+jcOZfk39dOL+m2ArMjAcUAUAJ3N12/6mgrbqACjc4XdKZZtbVzA6V\n1F/SAnffIGmzmQ3NXxQ+W9K0plbo7kxtNE2aNKniZegoE/uS/ZnkqbXKGQY61szWSjpW0kwzm5Wv\nuJdLeljSckmzJF3qDSW7VNI9klZJetndZ5e6fQBAeZrtAmqOuz8u6fEmlt0m6bYi8xdKOrLUbQIA\n2g53AmdALperdBE6DPZl22J/VpaV0m8Um5l5EssFAElmZvIKXAQGAKQMAQAAGUUAAEBGEQAAkFEE\nAABkFAEAABlFAABARhEAAJBRBAAAZBQBAAAZRQAAQEYRAACQUQQAAGQUAQAAGUUAAEBGEQAAkFEE\nAABkFAEAABlFAABARhEAAJBRBAAAZBQBAAAZRQAAQEYRAACQUQQAAGQUAQAAGUUAAEBGEQAAkFEE\nAABkFAEAABlFAABARhEAAJBRBAAAZBQBAAAZVXIAmNk4M1tmZrVmNrhgfnczm2tm75rZTxt9Z4iZ\nLTWzVWZ2dzkFBwCUp5wWwFJJYyXNbzR/u6SbJE0s8p0pks539/6S+pvZ6DK2DwAoQ8kB4O4r3H1l\nkflb3f2Pkv5VON/Mekraz90X5Gc9IGlMqdsHAJQn5jUAb/S+t6R1Be9r8vMAABXQubmFZjZHUo8i\ni25w9xlxihRUVVW9/zqXyymXy8XcHACkTnV1taqrq0v+vrk3PlFv5QrM5kq62t0XNZp/rqSj3f3b\n+fc9JT3t7gPz78+SNMLdLy6yTi+3XACQNWYmd7eWfr6tuoCKbfAD89x9vaTNZjbUzEzS2ZKmtdH2\nAQCtVHILwMzGSpos6UBJ70ha7O6n5Je9Kmk/SV0lvS1plLuvMLMhku6TtLekJ9x9QhPrpgUAAK3U\n2hZA2V1AMRAAANB6leoCAgCkDAEAABlFAABARhEAAJBRBAAAZBQBAAAZRQAAQEYRAACQUQQAAGQU\nAQAAGUUAAEBGEQAAkFEEAABkFAEAABlFAABARhEAAJBRBAAAZBQBAAAZRQAAQEYRAACQUQQAAGQU\nAQAAGUUAAEBGEQAAkFEEAABkFAEAABlFAABARhEAAJBRBAAAZBQBAKAi3npLWrNG2rWr0iXJrs6V\nLgCAbFmzRrrzTunBB6Vu3aRNm6SePaV+/aRDDvnwv337Sl27VrjQkbmHf83ad7vm9VtOEDPzJJYL\nQOmWLpV+9CNp5kzpggukK6+UevWSduyQ1q4NwfDqqx/+d/166aMfDYHQOBwOPVQ67DCpUwr7Mtav\nl558UvrDH8K/W7dKvXuHfdK79wen+nm9ejUfhmYmd29xjBAAkbz7rvTPf0qf+ET7p3pHUlsrTZ0a\nzhhHjpRuuUXaa69Klwot5S49+6x0xx3SwoXShAnSJZdI++/f8nXs2iXV1BQPh9WrpS1bpBNOCMfH\nyJHSxz+ezN/cjh3SH/8ozZ4dKv01a6QTT5ROPjlMBxwQ/jtrakLdUf+68P2GDWHfNQ6G+unUUwmA\nRJg8WZo4UTrwQOkLX2iY+vWrdMnSoa5OevRRqaoq/DCuuUb61a+kZcuk++6Tjjmm0iVEc+rqpN//\nXrr9dum116TvfEc699w44b12rfT002F66qlQ+Y8cGX5vI0eGirFSVq9uqPDnzZMGDAiV/ejR4Rju\n3MpO+Lq6sD+LhURNjfTkkwRAIvz4x9LGjdJll4WDsv4A7dat4eA84QSpR49KlzRZ6uqkadOkSZPC\nvvre96STTmo4o3v44XAWed55IRz23LOixUUjO3ZIDz0Uunr22ku69lrpy1+W9tijfbbvLq1aFX5z\nTz0lzZ0buo/qw+CEE6Tu3eNtf8uWsM36Sv+99xoq/BNPDCeEMbVbF5CZjZNUJWmApGPcfWF+/ihJ\nP5TUVdIOSd9x97n5ZUMk3SdpL0lPuPsVTaw79QFw++1hlMMddzTMc5deeqkhEObNCxe/6g/OESPC\n2W4WuUszZoSKv1OnUPGfemrxpvzGjaEbYeXK0Bo4+uh2Ly4a2bJF+sUvpJ/8JHR7XnddOKYr3RVT\nVye9+GJD6+DZZ0MXUX130bBh0r77Nv39Xbuk7dulbduan1avDhX+X/4iffazocI/+WTpqKPadx+0\nZwAMkFQn6b8kXe3ui/LzPyNpg7tvMLMjJP3B3fvkly2QdLm7LzCzJyRNdvfZRdad+gD4wQ9C+t92\nW9Ofqa2VlixpCIQ//Sn8eOoDYdgwaZ99Pvj5HTtaN+3cGS4a7bXXh6e99254veeelbmQ5i7NmiXd\nfHP4sd1yi3TGGbv/0biHM82rrpIuvFD67ndpDVTCpk3ST38qTZkSzq6vvVYaMqTSpWrajh2hkq5v\nISxcKH3yk+F4Klax19aG30mxqf43tPfeoZvppJPCPmguUGJr94vAZjZXBQHQaJlJel1SD0kHSnra\n3Qfml50pKefuFxf5XuoD4HvfC5Xv97/f8u/s2CEtWNAQCAsXhqZzfWVuFirzrl2lLl0aXjc3de4c\nyrF9e9PTtm3Sv/4V1lksKA44QBo0KPywjz46hFS5YeEuzZkTKv733gsV/5gxrV/v+vXSxRdLr7wi\n3X+/NHhweeVC89xDX/Py5dL06dJvfiONGxeud/XvX+nStd7WrWF0UufOxSv5Ll0q34ppjdYGQOz7\nAL4saaG77zSz3pLWFSyrkVTByzNx1da2vt+za9dw1j9sWOgK2b49VPz1FX7MflT3sK1iAbFpk7Ro\nUfjBT5oU3tcHQn0o9O/fssrbPYTbpEnSm2+GfvyvfKX0QOnZM1wz+PWvQ7P7kkukG2/s+OPGY3MP\nF1eXLw8X3pcvb5i6dZMOP1w67rjwPs3Xsbp1k4YOrXQpKqfZADCzOQpn743d4O4zdvPdIyTdLmlU\n6cVLr9ra8rsk6s/A24NZKO+ee0r/9m8fXn7iiQ2v33wzBMILL4TK96abpDfeCGff9YEwZEjoay2s\n2OfNC2f869eHin/8+LYJNTPpG98IXWcXXRRGV9x/v/TpT5e/7o6urk76xz8+XNG/9FLoyjjiiFDZ\nDx0aLrwffnjci6hoX80GgLuXVHmbWR9Jv5N0trv/PT+7RlKfgo/1yc8rqqqq6v3XuVxOuVyulKJU\nTCktgLTo3j0EQmEovPFG6LJauFB67DHp+uvDRfD6UFiyJIzbvvlm6Wtfa/3wt5bo1Su0Uh54QBo1\nSrr88lCOLl3afltptX17GKEyfXro+njppTCuvL6iP+64cE1l4MDsDkhIk+rqalVXV5f8/ba6BjCx\nYBTQ/pLmSZrk7tMaffZ5SRMkLZA0Ux34IvA110gf+Ui4KJZVr7/eEAq9eklf/3r7VcY1NaEi27Ah\ntAaOPLJ9tptEO3aE6y1Tp4aRVoMGhaGZgweHSr9Yiw/p1J6jgMZKmqxwcfcdSYvd/RQzu0nSdZJW\nFXx8lLu/XjAMdG+FYaATmlh36gNg4sTQNzpxYqVLkl3u0i9/GUL4yivDvzFaHi21a1cYtbRoUbjO\nM3y4dPzxcc60d+0K49GnTpUefzyc0Y8fH6639OzZ9ttDMvAoiIS46irp4IPDv6istWvDs2fefDN0\nTx18cPuXYds26cwzw2irq68OjwR49lnp+efD3eHDh4dp2LDw8LNS1NZKzzwTKv3HHgvrPfPMMEqn\n1HUiXZI2CiizOvI1gLTp2zf0e995p/S5z4Uz4vZ8lMTbb4d7G/r0kR55JIxQGpW/urZzp7R4cQiD\n+ruc99mnIQyGDw9n700NRayrk/7851DpP/qodNBB4Uz/ueekj32s/f4bkU60ACK57LLww7388kqX\nBIWmT5fOP1/62c+kr341/vbWrw/DU0eMkO66a/fDXd2lv/0tnMk/80wIhs2bQ1dRfSth0KBwd+vU\nqSE09t03nOmPHx9uakJ20QJICFoAyXTGGeGC6BlnhEdJ3HhjvBt9Vq8Od4eed17Lt2MWHhg2YEC4\niC2FC9r1YfDgg2HkziGHhAp/5kzpU59K181KSA5aAJFceGF4Jsi3vlXpkqCY9etDCAwYIN1zT9s/\nRmLJEum008Kw14suatt1b90a7lKl0kdjrW0BpPDPKKQDLYBk69kz3Ji2fXt47tKmTW237nnzwpn/\n3Xe3feUvhbtXqfzRFgiASAiA5OvWLfSj53LSsceGO2DLNW1aGHXz29+GIZdAkhEAkRAA6dCpk3Tr\nreHZRLlc+NN8pbr33vAsolmzwmMpgKQjACIhANLlnHPC2PlzzgmPNm4N9/B3H269NXT/JPlxyEAh\nRgFFQgCkz/Dh4QatL35RWrEi3Dewu/+HdXXhzx0++WT4bq9e7VNWoC3QAoiEAEinww4LN1YtXx5G\nCW3e3PRnd+6UvvnNcNPV/PlU/kgfAiASAiC99t9feuKJcAfx8cdLa9Z8+DNbt0pjx4bHS8yZw5Mz\nkU4EQCQEQLp16RKuBVxwQXh8xHPPNSx7663wKIfu3cNjJbp1q1w5gXIQAJEQAOlnJl1xRfhj56ef\nHoZ21tRIn/98GDZ63338rQGkGxeBIyEAOo7TTgt/p/n000PXz8SJ4e89cDMW0o4AiIQA6FiOOio8\nuvnFF6WTT650aYC2QQBEQgB0PD16pPsPoAONcQ0gEgIAQNIRAJEQAACSjgCIhAAAkHQEQCQEAICk\nIwAiIQAAJB0BEAkBACDpCIBICAAASUcAREIAAEg6AiASAgBA0hEAkRAAAJKOAIiEAACQdARAJAQA\ngKQjACIhAAAkHQEQCQEAIOkIgEgIAABJRwBEUldHAABINgIgEloAAJKOAIiEAACQdARAJLW1Uif2\nLoAEo4qKhBYAgKQrOQDMbJyZLTOzWjMbXDD/GDNbnJ/+ambjC5YNMbOlZrbKzO4ut/BJRgAASLpy\nWgBLJY2VNL/I/CHuPkjSSZJ+Zmb1VeEUSee7e39J/c1sdBnbTyx3RgEBSL6SA8DdV7j7yiLzt7l7\nXf7t3pLecfdaM+spaT93X5Bf9oCkMaVuP8nq6iSzMAFAUkW5BpDvBlomaZmkf8/P7i1pXcHHavLz\nOhy6fwCkQefmFprZHEk9iiy6wd1nNPW9/Fn+EWY2QNJsM6suq5QpQwAASINmA8DdR5WzcndfYWar\nJX1c4ey/T8HiPgqtgKKqqqref53L5ZTL5copSrsiAAC0h+rqalVXV5f8fXP3sgpgZnMlTXT3hfn3\n/SStc/ddZnaIpGckfcrdN5vZ85ImSFogaaakye4+u8g6vdxyVdI770h9+0qbN1e6JACyxMzk7i2+\n+ljOMNCxZrZW0rGSZprZrPyi4ZKWmNliSY9I+pa711eFl0q6R9IqSS8Xq/w7AloAANKg7BZADGlv\nAWzaJA0cKL3+eqVLAiBL2q0FgKbRAgCQBgRABAQAgDQgACIgAACkAQEQAQEAIA0IgAgIAABpQABE\nQAAASAMCIAICAEAaEAAREAAA0oAAiIAAAJAGBEAEBACANCAAIiAAAKQBARABAQAgDQiACAgAAGlA\nAERAAABIAwIgAgIAQBoQABEQAADSgACIgAAAkAYEQAQEAIA0IAAiIAAApAEBEEFdHQEAIPkIgAho\nAQBIAwIggtpaqRN7FkDCUU1FQAsAQBoQABEQAADSgACIgAAAkAYEQAQEAIA0IAAiIAAApAEBEAEB\nACANCIAICAAAaUAAREAAAEgDAiACAgBAGhAAERAAANKAAIiAAACQBgRABAQAgDQgACIgAACkAQEQ\nAQEAIA1KDgAzG2dmy8ys1swGF1l+sJltMbOrC+YNMbOlZrbKzO4uddtJRwAASINyWgBLJY2VNL+J\n5XdKmtlo3hRJ57t7f0n9zWx0GdtPLAIAQBqUHADuvsLdVxZbZmZjJL0iaXnBvJ6S9nP3BflZD0ga\nU+r2k4wAAJAGbX4NwMz2lXSNpKpGi3pLWlfwviY/r8MhAACkQefmFprZHEk9iiy6wd1nNPG1Kkk/\ncfetZmalFqyqqur917lcTrlcrtRVtTsCAEB7qK6uVnV1dcnfN3cvqwBmNlfS1e6+KP9+vqS++cX7\nS6qT9F1Jv5M0190H5j93lqQR7n5xkXV6ueWqpAkTpMMOk664otIlAZAlZiZ3b/GJd7MtgNZst/6F\nu3++oDCTJL3r7j/Pv99sZkMlLZB0tqTJbbT9RKEFACANyhkGOtbM1ko6VtJMM5vVgq9dKukeSask\nvezus0vdfpIRAADSoOQWgLs/Lunx3XzmlkbvF0o6stRtpgUBACANuBM4AgIAQBoQABEQAADSgACI\ngAAAkAYEQAQEAIA0IAAiIAAApAEBEAEBACANCIAICAAAaUAARFBXRwAASD4CIILaWqkTexZAwlFN\nRUAXEIA0IAAiIAAApAEBEAEBACANCIAICAAAaUAAREAAAEgDAiACAgBAGhAAERAAANKAAIiAAACQ\nBgRABAQAgDQgACIgAACkAQEQAQEAIA0IgAgIAABpQABEQAAASAMCIAICAEAaEAAREAAA0oAAiIAA\nAJAGBEAEBACANCAAIiAAAKQBARABAQAgDQiACAgAAGlAAERAAABIAwIgAgIAQBoQABEQAADSgACI\ngAAAkAYEQBtzl+rqCAAAyUcAtLG6OsksTACQZCUHgJmNM7NlZlZrZoML5vczs21mtjg//bxg2RAz\nW2pmq8zs7nILn0R0/wBIi3JaAEsljZU0v8iyl919UH66tGD+FEnnu3t/Sf3NbHQZ20+kJAZAdXV1\npYvQYbAv2xb7s7JKDgB3X+HuK1v6eTPrKWk/d1+Qn/WApDGlbj+pCICOjX3ZttiflRXrGsCh+e6f\najMblp/XW9K6gs/U5Od1KEkMAAAopnNzC81sjqQeRRbd4O4zmvjaPyX1dfe38tcGppnZEWWWMzUY\nAQQgLczdy1uB2VxJV7v7ouaWS1ov6Wl3H5iff5akEe5+cZHvlFcoAMgod2/xGMRmWwCt8P4GzexA\nSW+5e62ZfUxSf0mvuPvbZrbZzIZKWiDpbEmTi62sNf8BAIDSlDMMdKyZrZV0rKSZZjYrv2iEpBfN\nbLGkRyRd5O5v55ddKukeSasURgrNLr3oAIBylN0FBABIp0TdCWxmo81sRf5GsWsrXZ60M7NXzeyv\n+RFZC3b/DRQys/82s41mtrRgXnczm2NmK83sSTPbv5JlTJMm9meVma0ruHG0w90bFIOZ9TWzufmb\ncf/PzCbk57fq+ExMAJjZHpL+U9JoSYdLOsvMBla2VKnnknL5G/KOqXRhUuiXCsdjoeskzXH3T0h6\nKv8eLVNsf7qkOwtuHKVbuGV2SrrK3Y9Q6Ia/LF9ftur4TEwASDpG4brAq+6+U9JvJX2pwmXqCLig\nXiJ3f0bSW41mnyHp/vzr+9UBb2aMpYn9KXGMtpq7b3D3JfnXWyS9pHBfVauOzyQFQG9Jawver1MH\nvFGsnbmk/zWzF8zswkoXpoM4yN035l9vlHRQJQvTQXzbzF40s3vpUms9M+snaZCk59XK4zNJAcDV\n6LZ3vLsPknSKQhNxeKUL1JF4GEHBcVueKZIOlfQZhXuF/qOyxUkXM9tX0mOSrnD3dwuXteT4TFIA\n1EjqW/C+rz746Ai0kruvz/+7SdLjCt1sKM9GM+shvf98q9cqXJ5Uc/fXPE9hiDjHaAuZWReFyv9B\nd5+Wn92q4zNJAfCCwhNC+5lZV0njJU2vcJlSy8y6mdl++df7SDpJ4QmuKM90SefmX58raVozn8Vu\n5CupemPFMdoiZmaS7pW03N3vKljUquMzUfcBmNkpku6StIeke939hxUuUmqZ2aEKZ/1SuOP71+zP\n1jGzhxRubDxQoT/1Zkn/I+lhSQdLelXSVwtudEQziuzPSZJyCt0/LunvCjeObmxqHQjyD9mcL+mv\naujmuV7hKQstPj4TFQAAgPaTpC4gAEA7IgAAIKMIAADIKAIAADKKAACAjCIAACCjCAAAyCgCAAAy\n6v8B4PRSRLeXXnwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10adc7390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd.Series(ll).plot(ylim=[-150, -100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a future post, I hope to show how we can \"collapse\" the Gibbs sampler and sampling the assignment parameter without sampling the $\\pi$ and $\\theta$ values. This collapsed sampler can also be extended to the model with a Dirichet process prior that allows the number of clusters to be a parameter fit by the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Notation Helper\n",
    "\n",
    "* $N_k$, `state['suffstat'][k].N`: Number of points in cluster $k$.\n",
    "\n",
    "* $\\theta_k$, `state['suffstat'][k].theta`: Mean of cluster $k$.\n",
    "* $\\lambda_1$, `state['hyperparameters_']['mean']`: Mean of prior distribution over cluster means.\n",
    "* $\\lambda_2^2$, `state['hyperparameters_']['variance']` Variance of prior distribution over cluster means.\n",
    "* $\\sigma^2$, `state[cluster_variance_]`: Known, fixed variance of clusters. \n",
    "\n",
    "The superscript $(t)$ on $\\theta_k$, $pi_k$, and $z_i$ indicates the value of that variable at step $t$ of the Gibbs sampler."
   ]
  }
 ],
 "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
}