Explaining 'Explaining the Gibbs Sampler'

gibbs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explaining 'Explaining the Gibbs Sampler'\n",
    "\n",
    "Eamonn Bell, Columbia University <epb2125@columbia.edu>\n",
    "\n",
    "---\n",
    "\n",
    "Example taken from Casella, George, and Edward I. George. \"Explaining the Gibbs Sampler.\" The American Statistician 46, no. 3 (1992): 167-74.\n",
    "\n",
    "https://www.jstor.org/stable/2685208\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.special import comb, betaln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a joint of $X, Y$:\n",
    "\n",
    "$$f(x,y) \\propto \\binom{n}{x} y^{x + \\alpha - 1}(1 - y)^{n - x + \\beta -1 }$$\n",
    "\n",
    "We can get the marginal $f(x)$ analytically (the beta-binomial)\n",
    "\n",
    "$$f(x) = \\binom{n}{x} \\frac{\\Gamma \\left( \\alpha + \\beta \\right)}{\\Gamma \\left( \\alpha \\right)\\Gamma \\left( \\beta \\right)} \\frac{\\Gamma \\left( x + \\alpha \\right)\\Gamma \\left( n - x + \\beta \\right)}{\\Gamma \\left( \\alpha + \\beta + n \\right)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# http://www.channelgrubb.com/blog/2015/2/27/beta-binomial-in-python\n",
    "    \n",
    "def BetaBinomAnalytic(alpha,beta,n,x):\n",
    "    # log sum exp\n",
    "    part_1 = comb(n,x)\n",
    "    part_2 = betaln(x+alpha,n-x+beta)\n",
    "    part_3 = betaln(alpha,beta)\n",
    "    \n",
    "    result = (np.log(part_1) + part_2) - part_3\n",
    "    \n",
    "    return np.exp(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ALPHA = 2\n",
    "BETA = 4\n",
    "N = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.arange(0, 16, 1)\n",
    "y = [BetaBinomAnalytic(alpha=ALPHA, beta=BETA, n=N, x=x_i) for x_i in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f39a3a171d0>"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEECAYAAAA2xHO4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wFOXhB/DvXi5Rklvg7pJoQjyjCbbNRRAPSmocYs6Z\nOhXGSdXeVKmtnTjVAiVgX4SJpbVGoZVIhQy+psYxWj20YbRMxzrDy4wndQ7K0XC81DPEKDGS5BQv\nyQG5u+f3B79cOQmbDXvJbcL3M+OY5Z5Nvoly3+yzu89KQggBIiKi8zCkOgAREekbi4KIiBSxKIiI\nSBGLgoiIFLEoiIhIEYuCiIgUJaUofD4fVqxYgZqaGmzduvWc1w8dOoSHHnoId911Fz744INzXg+H\nw3jggQfwl7/8JRlxiIgoiTQXRSwWQ2NjI2pra1FfXw+Px4Njx44ljMnJycHSpUtx4403Dvs5Xn/9\ndZSUlIzq6/r9/gvOPFaYSR09ZgL0mYuZ1GEm9S4kl+aiCAQCyMvLQ05ODoxGI8rLy+H1ehPGZGdn\nw2azQZKkc/Zva2vDiRMnMHv27FF9XT3+R2AmdfSYCdBnLmZSh5nUS0lRBINBWK3W+LbFYkEwGFS1\nrxACL7/8Mu655x7wBnEiIn0ak5PZwx05DOedd97BnDlzYLFYxiIGERElgaR1raf//ve/2LJlC2pr\nawEgfjK7qqrqnLGbN2+Gw+HA/PnzAQAbN27EkSNHIEkSwuEwotEovvvd7+Luu+8+Z1+/359wyORy\nubTEJiK6aLnd7vjHdrsddrtdcbxR6xcsLi5GV1cXuru7YTab4fF4UFNTc97xZ/fS8uXL4x/v3LkT\nbW1tw5YEMPw309nZqTF9csmyjFAolOoYCZhJPT3mYiZ1mEm9/Pz8Uf+irbkoDAYDqqurUVdXByEE\nnE4nCgoK4Ha7UVRUBIfDgY8++gjr169Hf38/9u7diy1btqC+vl7rlyYionGgeeoplXhEMTJmUk+P\nuZhJHWZSLz8/f9T78M5sIiJSxKIgIiJFLAoiIlLEoiAiIkUsCiIiUsSiICIiRSwKIiJSxKIgIiJF\nLAoiIlLEoiAiIkUsCiIiUsSiICIiRSwKIiJSxKIgIiJFLAoiIlLEoiAiIkUsCiIiUsSiICIiRSwK\nIiJSZEzGJ/H5fGhqaoIQApWVlaiqqkp4/dChQ2hqakJHRwdWrFiB+fPnAwDa29vxwgsvIBwOw2Aw\n4Pvf/z5uuOGGZEQiIqIk0VwUsVgMjY2NWLNmDcxmM1avXo158+ZhxowZ8TE5OTlYunQp3n777YR9\nL730UixbtgyXX345vvjiC6xatQrXXXcdMjMztcYiIqIk0VwUgUAAeXl5yMnJAQCUl5fD6/UmFEV2\ndjYAQJKkhH0vv/zy+MdmsxlTp07FV199xaIgItIRzecogsEgrFZrfNtisSAYDI768wQCAUSj0YTy\nICKi1BuTk9lfP3IYyRdffIGGhgYsWbJkLOIQEZEGmqeeLBYLenp64tvBYBBms1n1/uFwGOvWrcNd\nd92F4uLi847z+/3w+/3xbZfLBVmWLyz0GMnIyGAmFfSYCdBnLmZSh5lGx+12xz+22+2w2+2K4zUX\nRXFxMbq6utDd3Q2z2QyPx4OamprzjhdCxD+ORCJ44oknUFFREb8S6nyG+2ZCoZC28EkmyzIzqaDH\nTIA+czGTOsyknizLcLlco9pHc1EYDAZUV1ejrq4OQgg4nU4UFBTA7XajqKgIDocDH330EdavX4/+\n/n7s3bsXW7ZsQX19PXbv3o3Dhw+jv78fO3fuhCRJWLJkCa688kqtsYiIKEkkcfav+BNMZ2dnqiMk\n0ONvEMyknh5zMZM6zKRefn7+qPfhndlERKSIRUFERIpYFEREpIhFQUREilgURESkKCmrxxKpEYtJ\nOHRoCjo6DLDZMlFSEoYkTdiL7oguGiwKGjeHDk3BwoXTMDgoIT1dYNs2wG4fSHUsIhoBp55oRLGY\nBL8/E9u2yfD7MyHE6NbyGtLenobBwTP7Dg5KaG9PS2ZMIhojPKKgESXrSKCwMIr0dBH/PIWF0TFI\nS0TJxqKgEQ13JDDCGmLDKikJY9s2oKMjHTbbIEpKwklOSkRjgUVBI0rWkYAkCdjtAygrkxEK8dwE\n0UTBoqARDR0JtLenobAwqosjgaErqM7OxCuoiMYGi4JGNHQkcCHTTWOFV1ARjR9e9UQTEq+gIho/\nLAqakIbOmwDgFVREY4xTTzQh6fG8CdFkxaKgCUmP502IJitOPRERkSIWBRERKUrK1JPP50NTUxOE\nEKisrERVVVXC64cOHUJTUxM6OjqwYsUKzJ8/P/7azp070dLSAgC4/fbbUVFRkYxIRESUJJqPKGKx\nGBobG1FbW4v6+np4PB4cO3YsYUxOTg6WLl2KG2+8MeHP+/r68Oabb2Lt2rV4/PHH8cYbb2BggNfC\nExHpieaiCAQCyMvLQ05ODoxGI8rLy+H1ehPGZGdnw2azQZISVx3dv38/Zs2ahczMTGRlZWHWrFnw\n+XxaIxH+t+Kr223QtOIrEZHmqadgMAir1RrftlgsCAQCF7xvMBjUGonAO5eJKHnG5PLYrx85nI8Q\n6tfm8fv98Pv98W2XywVZlkedbSxlZGToJlNHhyHhzuWOjnSUlekjm55+TmfTYy5mUoeZRsftdsc/\nttvtsI9wnbnmorBYLOjp6YlvB4NBmM1mVftardaEN//e3l6UlpYOO3a4byYUCl1A4rEjy7JuMtls\nmQkrvtpsg7pZsVVPPyfg7Ee0RmCz6WuBQb39rABmUkuPmYAzuVwu16j20VwUxcXF6OrqQnd3N8xm\nMzweD2pqas47/uyjiNmzZ+O1117DwMAAYrEYWltbsXjxYq2RCHz2w2hwmo5ImeaiMBgMqK6uRl1d\nHYQQcDqdKCgogNvtRlFRERwOBz766COsX78e/f392Lt3L7Zs2YL6+nqYTCbccccdWLVqFSRJwp13\n3omsrKxkfF8XPT77Qb1kPZiJaLKSxGhOFOhMZ2dnqiMk0OOhJjONzO/P/NoRxQndHFHo7WcFMJNa\neswEAPn5+aPeh2s90UWP03REylgUdNHjNB2RMq71REREilgURESkiEVBRESKWBRERKSIRUFERIpY\nFEREpIhFQUREilgURESkiEVBRESKWBRERKSIS3gQJcnQcy3a29NQWKiv51oQacGiIEoSPteCJitO\nPRElyXDPtSCaDFgURElSWBhFevqZqab0dIHCwmiKExElB6eeiJJk6LkWZ5+jIJoMWBRESTL0XAs+\nRpUmG049ERGRoqQcUfh8PjQ1NUEIgcrKSlRVVSW8HolE0NDQgLa2NsiyjJUrVyI7OxvRaBTPPPMM\njh49ilgshgULFpyzLxERpZbmI4pYLIbGxkbU1taivr4eHo8Hx44dSxizfft2mEwmbNy4EQsXLkRz\nczMAYPfu3YhEIli/fj3WrVuHd999Fz09PVojERFREmkuikAggLy8POTk5MBoNKK8vBxerzdhjNfr\nRUVFBQCgrKwMBw4cAABIkoRTp04hFovh1KlTSE9Px5QpU7RGmtBiMQl+fya2bZPh92dCCCnVkYjo\nIqd56ikYDMJqtca3LRYLAoHAeccYDAZkZmair68PZWVl8Hq9+NnPfobTp0/jJz/5CbKysrRGmtB4\n0xYR6c2YXPUkScq/BQtx5lrzQCCAtLQ0PPfcc+jr68OaNWtw7bXXIjc395x9/H4//H5/fNvlckGW\n5eQG1ygjI0Nzpo4OQ8JNWx0d6Sgru/DPmYxMyabHTIA+czGTOsw0Om63O/6x3W6HfYRL9TQXhcVi\nSTivEAwGYTabE8ZYrVb09vbCYrEgFoshHA7DZDLhvffew3XXXQeDwYCpU6fiG9/4Btra2oYtiuG+\nmVAopDV+UsmyrDmTzZaJ9HQRP6Kw2QYRCl34EUUyMiWbHjMB+szFTOowk3qyLMPlco1qH83nKIqL\ni9HV1YXu7m5EIhF4PB7MnTs3YYzD4cCuXbsAnDmBXVpaCgDIzs6On684efIkPvzwQ+Tn52uNNKGd\nuWnrBJ57LoRt207wpi0iSjlJDM0DaeDz+fDiiy9CCAGn04mqqiq43W4UFRXB4XBgcHAQmzZtQnt7\nO2RZRk1NDXJzc3Hy5Els3rw5fpVUZWUlFi1apPrrdnZ2ao2eVHr8DYKZ1NNjLmZSh5nUu5BfxpNS\nFKnCohgZM6mnx1zMpA4zqXchRcE7s4mISBGLgoiIFLEoiIhIEVePJdKhoceqdnQYYLNl8rGqlFIs\nCiId4h36pCeceiLSIT5WlfSERUGkQ3ysKukJp56IdGjosaodHemw2QZ5hz6lFIuCSIeGHqtaViZr\nWuuLKBk49URERIpYFEREpIhFQUREilgURESkiEVBRESKWBRERKSIRUFERIpYFEREpIhFQUREilgU\nRESkKClLePh8PjQ1NUEIgcrKSlRVVSW8HolE0NDQgLa2NsiyjJUrVyI7OxsA8PHHH+P5559HOByG\nwWDA2rVrYTRyZREiIr3Q/I4ci8XQ2NiINWvWwGw2Y/Xq1Zg3bx5mzJgRH7N9+3aYTCZs3LgR77//\nPpqbm7FixQrEYjE0NDTgF7/4BWw2G/r6+pCWxuWUiYj0RPPUUyAQQF5eHnJycmA0GlFeXg6v15sw\nxuv1oqKiAgBQVlaGAwcOAAD279+PK6+8EjabDQBgMpkgSZLWSERElESajyiCwSCsVmt822KxIBAI\nnHeMwWBAZmYm+vr68NlnnwEAHnvsMYRCIdxwww247bbbtEYiIqIkGpOTASMdFQhx5oEs0WgUR44c\nwdq1a5GRkYE//OEPuPrqq1FaWnrOPn6/H36/P77tcrkgy3Jyg2uUkZHBTCroMROgz1xaM0UiAvv2\nSTh61ICrrorh+usF0tK0HbVPxp/TWNBjpiFutzv+sd1uh91uVxyvuSgsFgt6enri28FgEGazOWGM\n1WpFb28vLBYLYrEYwuEwTCYTrFYrvvWtb8FkMgEA5syZg6NHjw5bFMN9M6FQSGv8pJJlmZlU0GMm\nQJ+5tGby+zO/9uztE5qfvT0Zf05jQY+ZgDO5XC7XqPbRfI6iuLgYXV1d6O7uRiQSgcfjwdy5cxPG\nOBwO7Nq1CwCwe/fueBHMnj0bHR0dOH36NKLRKA4ePIiCggKtkYjo//HZ25QMmo8oDAYDqqurUVdX\nByEEnE4nCgoK4Ha7UVRUBIfDAafTiU2bNmH58uWQZRk1NTUAgKysLCxatAirV6+GJEm4/vrrMWfO\nHM3fFBGdMfTs7aEjCj57my6EJIZOGExAnZ2dqY6QQI+Hmsyknh5zac0khISDB6egvT0NhYVRlJSE\nIUna/spPxp/TWNBjJgDIz88f9T68s41oEht69vYI5yqJFHEJDyIiUsSiICIiRSwKIiJSxKIgIiJF\nLAoiIlLEoiAiIkUsCiIiUsSiICIiRSwKIiJSxDuzkyAWk3Do0BR0dBhgs2UmZZkEIiK9YFEkwaFD\nU762lDM0L+VMRKQXnHpKAi7lTESTGYsiCYaWcgbApZyJaNLh1FMSlJSEsW0b0NGRDpttECUl4VRH\nIko6nou7eLEokmBoKeeyMhmhEM9N0OTEc3EXL049EZEqPBd38WJREJEqPBd38eLUExGpwnNxF6+k\nFIXP50NTUxOEEKisrERVVVXC65FIBA0NDWhra4Msy1i5ciWys7Pjr/f09ODBBx+Ey+XCokWLkhGJ\niJKM5+IuXpqnnmKxGBobG1FbW4v6+np4PB4cO3YsYcz27dthMpmwceNGLFy4EM3NzQmvv/TSS5gz\nZ47WKERENAY0F0UgEEBeXh5ycnJgNBpRXl4Or9ebMMbr9aKiogIAUFZWhtbW1oTXLrvsMlxxxRVa\noxAR0RjQXBTBYBBWqzW+bbFYEAwGzzvGYDAgKysLfX19OHXqFN566y384Ac/gBC8HpuISI/G5GS2\nJEmKrw+VgtvtxsKFC3HJJZck/Plw/H4//H5/fNvlckGW5SSkTZ6MjAxmUkGPmQB95mImdZhpdNxu\nd/xju90Ou92uOF5zUVgsFvT09MS3g8EgzGZzwhir1Yre3l5YLBbEYjGEw2GYTCYEAgF88MEHaG5u\nRn9/PwwGAzIyMnDLLbec83WG+2ZCoZDW+EklyzIzqaDHTIA+czGTOsyknizLcLlco9pHc1EUFxej\nq6sL3d3dMJvN8Hg8qKmpSRjjcDiwa9cuzJw5E7t370ZpaSkA4JFHHomP2bJlC6ZMmTJsSRARUepo\nLgqDwYDq6mrU1dVBCAGn04mCggK43W4UFRXB4XDA6XRi06ZNWL58OWRZPqdIiIhIvyQxgc8id3Z2\npjpCAj0eajKTenrMxUzqMJN6+fn5o96HS3gQEZEiFgURESliURARkSIuCkhE427oIUjt7WkoLIzy\nIUg6x6IgonHHhyBNLJx6IqJxx4cgTSwsCiIad3wI0sTCqSciGndDD0E6+xwF6ReLgojG3dBDkEZY\ni450glNPRESkiEVBRESKWBRERKSIRUFERIpYFEREpIhFQUREilgURESkiEVBRESKWBRERKSId2YT\n0YQ1tFx5R4cBNlsmlysfI0kpCp/Ph6amJgghUFlZiaqqqoTXI5EIGhoa0NbWBlmWsXLlSmRnZ+M/\n//kPXn31VUSjURiNRixevBilpaXJiEREFwEuVz4+NE89xWIxNDY2ora2FvX19fB4PDh27FjCmO3b\nt8NkMmHjxo1YuHAhmpubAQBTp07FqlWr8MQTT2DJkiVoaGjQGoeILiJcrnx8aC6KQCCAvLw85OTk\nwGg0ory8HF6vN2GM1+tFRUUFAKCsrAytra0AgMLCQkyfPh0AcMUVV2BwcBCRSERrJCK6SHC58vGh\neeopGAzCarXGty0WCwKBwHnHGAwGZGVloa+vDyaTKT7mX//6F6666ioYjTxtQkTqDC1X3tGRDptt\nkMuVj5ExeVeWJEnxdSESTzZ98sknePXVV/Hwww+fdx+/3w+/3x/fdrlckGVZW9Aky8jIYCYV9JgJ\n0GcuZhpZWRmwYIERp0/HAJhGHD9e9PZzOpvb7Y5/bLfbYR9hvXfNRWGxWNDT0xPfDgaDMJvNCWOs\nVit6e3thsVgQi8UQDofjRxO9vb1Yv349li1bhtzc3PN+neG+mVAopDV+UsmyzEwq6DEToM9czKQO\nM6knyzJcLteo9tF8jqK4uBhdXV3o7u5GJBKBx+PB3LlzE8Y4HA7s2rULALB79+74lU39/f1Yt24d\nFi9ejGuuuUZrFCIiGgOajygMBgOqq6tRV1cHIQScTicKCgrgdrtRVFQEh8MBp9OJTZs2Yfny5ZBl\nGTU1NQCAd955B59//jnefPNNvPHGG5AkCbW1tZg6darmb4yIiJJDEl8/YTCBdHZ2pjpCAj0eajKT\nenrMxUzqMJN6+fn5o96HS3gQEZEiFgURESniTQtERPjfulHt7WkoLIxy3aizsCiIiMB1o5Rw6omI\nCFw3SgmLgogIXDdKCaeeiIjwv3Wjzj5HQWewKIiIAEiSgN0+gBGWPbooceqJiIgUsSiIiEgRi4KI\niBTxHAURURIN3bjX0WGAzZY5KW7cY1EQESXRZLxxj1NPRERJNBlv3GNREBEl0WS8cY9TT0RESTR0\n415HRzpstsFJceMei4KIKImGbtwrK5MRCk3scxNDWBRERDqll6XPWRRERDqllyuoklIUPp8PTU1N\nEEKgsrISVVVVCa9HIhE0NDSgra0Nsixj5cqVyM7OBgC0tLRgx44dSEtLw7333ovZs2cnIxIR0YQ3\n3BVUqViLSvNVT7FYDI2NjaitrUV9fT08Hg+OHTuWMGb79u0wmUzYuHEjFi5ciObmZgDAp59+it27\nd2PDhg1YvXo1XnjhBQgxsW9MISJKFr1cQaW5KAKBAPLy8pCTkwOj0Yjy8nJ4vd6EMV6vFxUVFQCA\nsrIyHDhwAACwZ88e3HDDDUhLS0Nubi7y8vIQCAS0RiIimhTOXEF1As89F8K2bSdSdgWV5qmnYDAI\nq9Ua37ZYLOe82Z89xmAwIDMzE319fQgGg7jmmmsS9g0Gg1ojERFNCnpZ+nxMTmZLkqRq3HDTTOfb\n1+/3w+/3x7ddLhdkWb6wgGMkIyODmVTQYyZAn7mYSR1mGh232x3/2G63wz5CE2kuCovFgp6envh2\nMBiE2WxOGGO1WtHb2wuLxYJYLIaBgQGYTCZYrdaEfXt7e8/Zd8hw30woFNIaP6lkWWYmFfSYCdBn\nLmZSh5nUk2UZLpdrVPtoPkdRXFyMrq4udHd3IxKJwOPxYO7cuQljHA4Hdu3aBQDYvXs3SktLAQBz\n587F+++/j0gkguPHj6OrqwvFxcVaIxERURJpPqIwGAyorq5GXV0dhBBwOp0oKCiA2+1GUVERHA4H\nnE4nNm3ahOXLl0OWZdTU1AAACgoK8J3vfAcrV66E0WjEfffdp3raioiIxockJvD1qJ2dnamOkECP\nh5rMpJ4eczGTOsykXn5+/qj34eqxRESkiEVBRESKWBRERKSIRUFERIpYFEREpIhFQUREilgURESk\niEVBRESKWBRERKSIRUFERIpYFEREpIhFQUREilgURESkiEVBRESKWBRERKSIRUFERIpYFEREpIhF\nQUREijQ9M7uvrw9//vOf0d3djdzcXKxcuRKZmZnnjNu5cydaWloAALfffjsqKipw+vRpPPnkk/j8\n889hMBjgcDhw9913a4lDRERjQNMRxdatW3Httdfiqaeegt1uj5fB2fr6+vDmm29i7dq1ePzxx/HG\nG29gYGAAAHDbbbdhw4YN+NOf/oQjR47A5/NpiUNERGNAU1Hs2bMHFRUVAICbbroJXq/3nDH79+/H\nrFmzkJmZiaysLMyaNQs+nw8ZGRkoKSkBAKSlpeGqq65CMBjUEoeIiMaApqI4ceIEpk+fDgCYPn06\nvvrqq3PGBINBWK3W+LbFYjmnEPr7+7F3716UlpZqiUNERGNgxHMUjz76KE6cOBHfFkJAkiT88Ic/\nVPUFhBCKr8diMWzcuBG33norcnNzVX1OIiIaPyMWxW9/+9vzvjZ9+nR8+eWX8X9PmzbtnDFWqxV+\nvz++3dvbm3Dk8OyzzyIvLw/f+973FHP4/f6Ez+NyuZCfnz9S/HEny3KqI5yDmdTTYy5mUoeZ1HO7\n3fGP7XY77Ha74nhNU08OhwM7d+4EcObKprlz554zZvbs2WhtbcXAwAD6+vrQ2tqK2bNnAwBee+01\nhMNh3HvvvSN+LbvdDpfLFf/n7G9UL5hJHT1mAvSZi5nUYSb13G53wnvpSCUBaLw8tqqqChs2bMCO\nHTuQnZ2NBx98EADQ1taGd999F/fffz9MJhPuuOMOrFq1CpIk4c4770RWVhaCwSBaWlowY8YM/OY3\nv4EkSbjlllvgdDq1RCIioiTTVBQmk2nYqamrr74a999/f3z7pptuwk033ZQwxmKx4PXXX9fy5YmI\naByk/f73v/99qkNcKD2e/GYmdfSYCdBnLmZSh5nUG20uSYx0WRIREV3UuNYTEREpYlEQEZEiTSez\nU6m5uRl79+6F0WjEZZddhiVLlgy7IOF48Pl8aGpqghAClZWVqKqqSkmOs/X29qKhoQFffvklDAYD\nbr75Ztx6662pjgXgzE2Wq1evhsViwUMPPZTqOBgYGMAzzzyDTz75BJIk4ec//zlmzpyZ0kx///vf\nsWPHDkiSBJvNhiVLlsBoHP+/rk8//TT+/e9/Y9q0aVi/fj0A9YuBjmemVL8fDJdpyFtvvYVXXnkF\njY2NMJlMKc/0j3/8A++88w7S0tJw/fXXY/HixSN/MjFB7d+/X0SjUSGEEM3NzeKVV15JSY5oNCqW\nLVsmjh8/LgYHB8WvfvUr8emnn6Yky9m++OILcfToUSGEEOFwWCxfvlwXuYQQ4u233xZPPfWUWLdu\nXaqjCCGEaGhoENu3bxdCCBGJRER/f39K8/T29oqlS5eKwcFBIYQQTz75pNi5c2dKshw6dEgcPXpU\n/PKXv4z/2csvvyy2bt0qhBCipaVFNDc3pzxTqt8PhsskhBA9PT2irq5OLFmyRIRCoZRnOnDggHj0\n0UdFJBIRQghx4sQJVZ9rwk49zZo1CwbDmfgzZ85Eb29vSnIEAgHk5eUhJycHRqMR5eXlwy6OON6m\nT5+OwsJCAMCll16KGTNm6GLRxd7eXuzbtw8333xzqqMAAMLhMA4fPozKykoAZxaoTNWR6dlisRhO\nnjyJaDSKU6dOwWw2pyTHN7/5TWRlZSX8mZrFQMc7U6rfD4bLBAAvvfQS7rnnnnHNMmS4TP/85z9R\nVVWFtLQ0AMDUqVNVfa4JO/V0th07dqC8vDwlX3u4RQ8DgUBKspzP8ePH8fHHH6d8OgX431+coaXm\nU+3zzz+HLMvYvHkzPv74Y1x99dX46U9/ioyMjJRlslgsWLRoEZYsWYJLLrkEs2bNwqxZs1KW5+vU\nLAaaSql8Pzjbnj17YLVaYbPZUh0l7rPPPsPBgwfx17/+FRkZGfjRj36EoqKiEffTdVEoLUg4tFzI\n3/72N6SlpeHGG29MVcxzSJKU6ghxJ0+exJNPPol7770Xl156aUqzDM2XFhYWwu/3j7hg5HiIxWI4\nevQoqqurUVRUhKamJmzduhUulytlmfr7+7Fnzx5s3rwZmZmZqK+vx3vvvaer/8f1Si/vB6dPn0ZL\nSwsefvjh+J/p4f/3aDSKgYEBPPbYYwgEAtiwYQMaGhpG3E/XRaG0ICFwZn2pffv2Yc2aNeOU6FwW\niwU9PT3x7WAwmLJpgq+LRqOor6/HggULMG/evFTHweHDh7Fnzx7s27cPp0+fRjgcRkNDA5YtW5ay\nTBaLBVarNf5bVVlZGbZu3ZqyPADQ2tqK3Nzc+InP+fPn48iRIyl/8xuiZjHQVNDD+8GQrq4uHD9+\nHL/+9a8hhEAwGMSqVavw+OOPp/TnlZ2djW9/+9sAgOLiYkiShFAoNOLihbouCiU+nw9vvfUWHnnk\nEaSnp6csR3FxMbq6utDd3Q2z2QyPx4OampqU5Tnb008/jYKCAt1c7XT33XfHH3d78OBBvP322ykt\nCeDMm54ZZ4+GAAABbElEQVTVakVnZyfy8/PR2tqKgoKClGbKzs7Ghx9+iNOnTyM9PR2tra2qpgfG\nihAi4bfhocVAq6qqzrsY6Hhn0sP7wdmZbDYbnn/++fhrS5cuxR//+Mdxverp65kAYN68eThw4ABK\nSkrQ2dmJaDSqaoXbCXtn9vLlyxGJROLf5MyZM3HfffelJIvP58OLL74IIQScTqcuLo89fPgwfve7\n38Fms0GSJEiShLvuugvXXXddqqMB+F9R6OHy2Pb2djz77LOIRCIpv9R6yJYtW/D+++8jLS0NhYWF\neOCBB1JyeexTTz2FgwcPIhQKYdq0aXC5XJg3bx42bNiAnp6e+GKgw53IHc9MLS0tKX0/GC7T0AUS\nALBs2TKsW7duXItiuEwLFizA5s2b0d7ejvT0dPz4xz+OP2lUyYQtCiIiGh8T9vJYIiIaHywKIiJS\nxKIgIiJFLAoiIlLEoiAiIkUsCiIiUsSiICIiRSwKIiJS9H/2cRkZwzwvJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f39a3876eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We might also want to use `scipy` rvs interface to construct something that will give us samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def BetaBinomialSampler(alpha, beta, n, size):\n",
    "    \n",
    "    # draw size samples from the beta \n",
    "    beta_samples = sp.stats.beta.rvs(a=alpha,b=beta,size=size)\n",
    "    \n",
    "    beta_binomial_samples = []\n",
    "    \n",
    "    for p_from_beta in beta_samples:\n",
    "        sample = sp.stats.binom(n=n, p=p_from_beta).rvs()\n",
    "        beta_binomial_samples.append(sample)\n",
    "        \n",
    "    return np.array(beta_binomial_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEECAYAAAAifS8cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFz5JREFUeJzt3XtsU/fdx/HPsUMShbhxLg1qiKL0QoWwoKyE0ctGuKxC\nQkjNH5ulwrpRbdIzLi2DboOKIUCjHSmXJBqDVqhqq4EmUmmJCtq0/UFcbZ2qkgEbGNKOliJaCgmB\nJA4hhNjn+QM1z0NxsJOcE8e/vl//xDY+X38MyYeTn4+PLdu2bQEAjONJdQAAgDsoeAAwFAUPAIai\n4AHAUBQ8ABiKggcAQ2Ukc6eenh699tprOn/+vCzL0rJly3TfffeptrZWbW1tKi4u1urVq5WTk+N2\nXgBAkpLag3/zzTf1rW99SzU1Ndq2bZsmTpyoxsZGTZ06VXV1dQoEAmpoaEjqAcPh8IgCjxZyOicd\nMkrkdBo5nTWcnAkL/vr162ppadHcuXMlSV6vVzk5OWpublZlZaUkac6cOTpy5IhrIVOBnM5Jh4wS\nOZ1GTmcNJ2fCJZpLly7J5/Np9+7dOnfunB544AEtXbpUnZ2d8vv9kiS/36+urq6hJwYAuCbhHnws\nFtPZs2e1YMECVVdXKysrS42NjaORDQAwAlaic9F0dHTo17/+tXbt2iVJamlpUWNjoy5duqSNGzfK\n7/ero6NDmzdvVk1NzR3bh8Ph2361CAaDDj8FAPhmqK+vH7gcCAQUCATuev+ESzR+v1+FhYW6cOGC\nSkpKdOLECZWWlqq0tFShUEhVVVUKhUKqqKiIu328EBcuXEjmuaSUz+dTJBJJdYyE0iFnOmSUyOk0\ncjqrpKRkyDvISR0m+dxzz+l3v/ud+vv7NWHCBC1fvlyxWEw1NTVqampSUVGR1qxZM6zQAAB3JFyi\ncQN78M5Jh5zpkFEip9PI6aySkpIhb8M7WQHAUBQ8ABiKggcAQyX1Iiuc4b16WbrS5ujMmxNKpNw8\nR2cCMAMFP5qutKlv61pHR3rX76DgAcTFEg0AGIqCBwBDsUST5myPR95PTjs3sOBeRfOLnJsHIGUo\n+HQX6VRf7SbHxmWuq5YoeMAILNEAgKEoeAAwFAUPAIai4AHAUBQ8ABiKggcAQ1HwAGAoCh4ADEXB\nA4CheCcrbmNlZDh76gNxSmMgVSh43C7Spb66zY6O5JTGQGqwRAMAhqLgAcBQFDwAGIqCBwBDUfAA\nYCgKHgAMRcEDgKEoeAAwVFJvdFqxYoVycnJkWZa8Xq9++9vfqru7W7W1tWpra1NxcbFWr16tnJwc\nt/MCAJKUVMFblqWNGzcqNzd34LbGxkZNnTpVTz/9tBobG9XQ0KAlS5a4FhQAMDRJLdHYti3btm+7\nrbm5WZWVlZKkOXPm6MiRI86nAwAMW9J78C+//LIsy9L3vvc9zZ8/X52dnfL7/ZIkv9+vrq4uV4MC\nAIYmqYLfsmXLQIlv2bJFJSUlbucCAIxQUgX/1Z76Pffco5kzZ+rMmTPy+/3q6OgY+JqXF/9sgeFw\nWOFweOB6MBiUz+dzILq7MjMzHc95w+vCyTsth8dZDg+U5PFY39h/czeQ01npklOS6uvrBy4HAgEF\nAoG73j9h49y4cUO2bSs7O1u9vb36z3/+o+9///uaMWOGQqGQqqqqFAqFVFFREXf7eCEikUgyzyWl\nfD6f4zm90X5H50mS7MR3GdI42+GBkmIx+xv7b+4GcjornXIGg8EhbZOw4Ds7O7Vt2zZZlqVoNKrv\nfve7euSRR/Tggw+qpqZGTU1NKioq0po1a4YdHADgvIQFX1xcrG3btt1xe25urjZs2OBKKADAyPGJ\nTnCd7fE4/jGAKrhX0fwiZ2cChqHg4b5Ip/pqNzk6MnNdtUTBA3fFuWgAwFAUPAAYioIHAENR8ABg\nKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQvJMVacnKyHD89Ac3J5RIufFPew2kIwoe6SnSpb66zY6O\n9K7fQcHDKCzRAIChKHgAMBQFDwCGouABwFAUPAAYioIHAENR8ABgKAoeAAxFwQOAoSh4ADAUBQ8A\nhqLgAcBQFDwAGMqIs0l6r16WrrQ5OpNTxwJId0YUvK60qW/rWkdHcupYAOmOJRoAMFTSe/CxWEwv\nvfSSCgoKtHbtWrW2tqqurk7d3d26//779fzzz8vr9bqZdVTZHo/jnxhk9d90dB4A3E3SBf/nP/9Z\nEydO1PXr1yVJ+/fv16JFi/T4449r7969Onz4sJ566inXgo66SKf6ajc5OjJr1UZH5wHA3SS1RNPe\n3q5jx45p/vz5A7edPHlSs2bNkiRVVlbqww8/dCchAGBYkir4t99+W88++6wsy5IkRSIR5ebmyuO5\ntXlhYaGuXr3qXkoAwJAlXKI5evSo8vLyVF5ernA4LEmybVu2bd92v6/K/+vC4fDAdpIUDAbl8/lG\nkvkON7wuHAwU/+mMbOQgf0cjG+rwuDTIKLmT0+OxHP/edENmZiY5HZQuOSWpvr5+4HIgEFAgELjr\n/RM2Y0tLi5qbm3Xs2DH19fXp+vXreuutt9TT06NYLCaPx6P29nbl5+fH3T5eiEgkksxzSZo32u/o\nPEmSnfguQx5puzHU4XFpkFFyJ2csZjv+vekGn89HTgelU85gMDikbRIW/OLFi7V48WJJ0qlTp3Tw\n4EG98MILqqmp0QcffKAnnnhC7733nioqKoaXGgDgimEfB79kyRIdOnRIq1atUnd3t+bNm+dkLgDA\nCA1p8XrKlCmaMmWKJKm4uFivvPKKK6EAACPHO1kBwFAUPAAYioIHAENR8ABgKAoeAAxFwQOAoSh4\nADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKEoeAAwFAUPAIai4AHAUBQ8ABiKggcA\nQ1HwAGAoCh4ADEXBA4ChKHgAMBQFDwCGouABwFAUPAAYioIHAENR8ABgqIxEd7h586Y2btyo/v5+\nRaNRPfbYY/rBD36g1tZW1dXVqbu7W/fff7+ef/55eb3e0cgMuML2eOT95LSzQwvuVTS/yNmZQJIS\nFvy4ceO0ceNGZWVlKRaLacOGDZo+fboOHTqkRYsW6fHHH9fevXt1+PBhPfXUU6ORGXBHpFN9tZsc\nHZm5rlqi4JEiSS3RZGVlSbq1Nx+NRmVZlsLhsGbNmiVJqqys1IcffuheSgDAkCXcg5ekWCymdevW\n6dKlS1qwYIEmTJig8ePHy+O59f9DYWGhrl696mpQAMDQJFXwHo9Hr776qnp6erR9+3Z98cUXd9zH\nsqy424bDYYXD4YHrwWBQPp9vmHHju+FN6mkMTfynM7KRg/wdjWyow+PSIKOUPjm93gzlOPz9npmZ\n6fjPkBvI6bz6+vqBy4FAQIFA4K73H1Iz5uTkaMqUKfr444917do1xWIxeTwetbe3Kz8/P+428UJE\nIpGhPGxC3mi/o/MkSbYLI203hjo8Lg0ySumTMxrtd/z73efzOT7TDeR0ls/nUzAYHNI2Cdfgu7q6\n1NPTI0nq6+vTiRMnVFpaqkAgoA8++ECS9N5776miomIYkQEAbkm4B9/R0aHf//73isVism1bTzzx\nhB599FGVlpaqtrZWBw4cUHl5uebNmzcaeQEASUpY8GVlZaqurr7j9uLiYr3yyiuuhAIwuJtffi7v\npQvODuV4fSO58OokADfFLreqb+taR2dyvL6ZOFUBABiKggcAQ1HwAGAoCh4ADEXBA4ChKHgAMBQF\nDwCGouABwFAUPAAYioIHAENR8ABgKAoeAAxFwQOAoTibJOAiKyND3k9OOzrTjkYdnQdzUfCAmyJd\n6qvb7OjI7J9vcnQezMUSDQAYioIHAENR8ABgKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQKXmjk+Xk\nO/Esy7lZAGCQlBR87NV1js0aV7VEdmaWY/MAwBSpKfhPP3Jslh3plAqLHZsHAKZgDR4ADJVwD769\nvV27du1SR0eHPB6P5s+fr4ULF6q7u1u1tbVqa2tTcXGxVq9erZycnNHIDABIQsKC93q9+vGPf6zy\n8nL19vZq7dq1euSRR9TU1KSpU6fq6aefVmNjoxoaGrRkyZLRyAwASELCJRq/36/y8nJJUnZ2tiZO\nnKj29nY1NzersrJSkjRnzhwdOXLE1aAAgKEZ0hp8a2urzp07p4cfflidnZ3y+/2Sbv0n0NXV5UpA\nAMDwJH0UTW9vr3bu3KmlS5cqOzs76QcIh8MKh8MD14PB4NASJmBZHnm8LhwM5MLh9ZYbx+w7PDId\nMkrkdJrXm6Ecn8/RmbFLF5TZetHRmZ6iYo27r9TRmZmZmfI5/NzdUl9fP3A5EAgoEAjc9f5JNWM0\nGtWOHTs0e/ZszZw5U9KtvfaOjo6Br3l5eXG3TSbESNh2TNFovwuDXRhpuzHU4XFpkFEip9Oi0X5F\nIhFHZ2a2XtT1l190dua6avXmxu+a4fL5fI4/dzf4fL4h7yAntUSzZ88elZaWauHChQO3zZgxQ6FQ\nSJIUCoVUUVExpAcGALgr4R58S0uL/v73v6usrEy/+tWvZFmWnnnmGVVVVammpkZNTU0qKirSmjVr\nRiMvACBJCQt+8uTJOnDgQNw/27Bhg+OBAADO4J2sAGAoCh4ADEXBA4ChKHgAMBQFDwCGouABwFAU\nPAAYioIHAENR8ABgKAoeAAyVkg/dBjC2WBkZ8n5y2tGZdjTq6DwMHQUPQIp0qa9us6Mjs3++ydF5\nGDqWaADAUBQ8ABiKggcAQ1HwAGAoCh4ADEXBA4ChKHgAMBQFDwCGouABwFAUPAAYioIHAENR8ABg\nKAoeAAxFwQOAoSh4ADAUBQ8Ahkr4gR979uzR0aNHlZeXp+3bt0uSuru7VVtbq7a2NhUXF2v16tXK\nyclxPSwAIHkJ9+Dnzp2r9evX33ZbY2Ojpk6dqrq6OgUCATU0NLgWEAAwPAkLfvLkyRo/fvxttzU3\nN6uyslKSNGfOHB05csSddACAYRvWGnxnZ6f8fr8kye/3q6ury9FQAICR40VWADBUwhdZ4/H7/ero\n6Bj4mpeXN+h9w+GwwuHwwPVgMDichxyUZXnk8Q7raSQY7MJIy42hDo9Lg4wSOR0fmSY5vd4M5fh8\njs7MzMyUz+GZbqmvrx+4HAgEFAgE7nr/pJrRtm3Ztj1wfcaMGQqFQqqqqlIoFFJFRcWg2yYTYiRs\nO6ZotN+FwS6MtN0Y6vC4NMgokdPxkWmSMxrtVyQScXSmz+dzfKYbfD7fkHeQExZ8XV2dTp06pUgk\nomXLlikYDKqqqko1NTVqampSUVGR1qxZM+zQAAB3JCz4VatWxb19w4YNjocBADiHF1kBwFAUPAAY\nioIHAENR8ABgKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKFc\n+CgkAHCHlZEh7yenHZ15c0KJlDv4p9KlMwoeQPqIdKmvbrOjI73rdxhb8CzRAIChKHgAMBRLNAC+\n0WyPx/F1fRXcq2h+kbMzh4GCB/DNFulUX+0mR0dmrquWxkDBs0QDAIai4AHAUBQ8ABiKggcAQ1Hw\nAGAoCh4ADEXBA4ChKHgAMNSI3uh0/PhxvfXWW7JtW3PnzlVVVZVTuQAAIzTsPfhYLKY33nhD69ev\n144dO/T+++/riy++cDIbAGAEhl3wZ86c0X333ad7771XGRkZevLJJ3XkyBEnswEARmDYBX/lyhUV\nFhYOXC8oKNCVK1ccCQUAGDlHTzZmWVZS9xu3+H8ce0zP/Q8r1tXh2DwAMIVl27Y9nA0//vhjvfPO\nO1q/fr0kqbGxUZLueKE1HA4rHA4PXA8Gg8PNCgDfaPX19QOXA4GAAoHA3TewhykajdorV660W1tb\n7Zs3b9q/+MUv7PPnzyfc7sCBA8N9yFFFTuekQ0bbJqfTyOms4eQc9hKNx+PRT37yE23ZskW2bWve\nvHkqLS0d7jgAgMNGtAY/ffp01dXVOZUFAOAg76ZNmzaN9oMWFxeP9kMOCzmdkw4ZJXI6jZzOGmrO\nYb/ICgAY2zgXDQAYioIHAEM5+kanu0mHE5O1t7dr165d6ujokMfj0fz587Vw4cJUxxpULBbTSy+9\npIKCAq1duzbVceLq6enRa6+9pvPnz8uyLC1btkyTJk1Kdaw7HDp0SE1NTbIsS2VlZVq+fLkyMkbt\nx2NQe/bs0dGjR5WXl6ft27dLkrq7u1VbW6u2tjYVFxdr9erVysnJGXM59+3bp3/961/KyMjQhAkT\ntHz58jGZ8yvvvvuu9u/frzfeeEO5ubkpSjh4xr/85S/661//Kq/Xq0cffVRLlixJPMzhQzXjinfM\n/Oeffz4aDz0kV69etc+ePWvbtm1fv37dfuGFF8Zkzq8cPHjQrqurs7du3ZrqKIPatWuXffjwYdu2\nbbu/v9++du1aihPdqb293V6xYoV98+ZN27Zte+fOnXYoFEpxqltOnz5tnz171n7xxRcHbvvDH/5g\nNzY22rZt2w0NDfa+fftSFW9AvJz//ve/7Wg0atu2be/bt8/ev39/quINiJfTtm378uXL9pYtW+zl\ny5fbkUgkReluiZfx5MmT9m9+8xu7v7/ftm3b7uzsTGrWqCzRpMuJyfx+v8rLyyVJ2dnZmjhx4pg9\nv057e7uOHTum+fPnpzrKoK5fv66WlhbNnTtXkuT1elO+BzeYWCym3t5eRaNR3bhxQ/n5+amOJEma\nPHmyxo8ff9ttzc3NqqyslCTNmTNnTPwsxcs5bdo0eTy3KmbSpElqb29PRbTbxMspSW+//baeffbZ\nFCS6U7yMf/vb31RVVSWv1ytJuueee5KaNSq/g8Y7MdmZM2dG46GHrbW1VefOnRuTywnS/31D9vT0\npDrKoC5duiSfz6fdu3fr3LlzeuCBB/Tcc88pMzMz1dFuU1BQoEWLFmn58uXKysrStGnTNG3atFTH\nGlRnZ6f8fr+kWzslXV1dKU6UWFNTk5588slUx4irublZhYWFKisrS3WUQX355Zc6deqU/vjHPyoz\nM1M//OEP9eCDDybcLmUvsiZ7YrJU6O3t1c6dO7V06VJlZ2enOs4dvlqfKy8vl23bssfoka6xWExn\nz57VggULVF1draysrIFzFo0l165dU3Nzs3bv3q3XX39dvb29+sc//pHqWMb405/+JK/Xq+985zup\njnKHvr4+NTQ03HaOrLH48xSNRtXT06OXX35ZS5YsUU1NTVLbjUrBFxQU6PLlywPXr1y5MmZ+Bf66\naDSqHTt2aPbs2Zo5c2aq48TV0tKi5uZmrVy5UnV1dQqHw9q1a1eqY92hoKBAhYWFA3sajz32mD79\n9NMUp7rTiRMnVFxcrNzcXHk8Hs2aNUsfffRRqmMNyu/3q6Pj1hlUOzo6lJeXl+JEgwuFQjp27JhW\nrVqV6ihxXbx4Ua2trfrlL3+pFStW6MqVK1q3bp06OztTHe02RUVF+va3vy1Jeuihh2RZliKRSMLt\nRmWJ5qGHHtLFixfV1tam/Px8vf/++2P2H3zPnj0qLS0d00fPLF68WIsXL5YknTp1SgcPHtTKlStT\nnOpOfr9fhYWFunDhgkpKSnTixIkxeb6ioqIi/fe//1VfX5/GjRunEydOJPXr72j5+m9pM2bMUCgU\nUlVVlUKhkCoqKlKY7v98Pefx48f17rvvavPmzRo3blwKk93u/+csKyvT3r17B/5sxYoVqq6uTulR\nNNKdf5czZ87UyZMnNWXKFF24cEHRaFQ+ny/hnFF7J+vx48f15ptvDpyYbCweJtnS0qKNGzeqrKxM\nlmXJsiw988wzmj59eqqjDeqrgh+rh0l+9tlnev3119Xf3z9mDpWL55133tE///lPeb1elZeX62c/\n+9mYOEyyrq5Op06dUiQSUV5enoLBoGbOnKmamhpdvnxZRUVFWrNmTdwXDlOds6GhQf39/QNFNGnS\nJP30pz8dczm/OghAklauXKmtW7emtODjZZw9e7Z2796tzz77TOPGjdOPfvQjTZkyJeEsTlUAAIbi\nnawAYCgKHgAMRcEDgKEoeAAwFAUPAIai4AHAUBQ8ABiKggcAQ/0veqAZaNbeuxwAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f39a3baf240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(BetaBinomialSampler(alpha=ALPHA, beta=BETA, n=N, size=500),bins=np.arange(16))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, it is possible to get the marginal analytically ($p(x) = \\frac{p(x,y)}{p(y \\mid x)}$). We had to do that to compute the first function above.\n",
    "\n",
    "The idea behind Gibbs sampling is to get draws from a distribution that is \"similar enough\" to this marginal, so that when it's hard to compute the distribution analytically we aren't stuck.\n",
    "\n",
    "### Gibbs sequence\n",
    "\n",
    "To do this, we consider a \"Gibbs sequence\" of random variables that are related to the conditional distributions $f(x \\mid y)$ and $f(y \\mid x)$ in a special way.\n",
    "\n",
    "$$Y'_0, X'_0, Y'_1, X'_1, \\ldots, Y'_k, X'_k$$\n",
    "\n",
    "Where some initial value $Y'_0 = y'_0$ is specified, and the following relations hold:\n",
    "\n",
    "$$\n",
    "X'_j \\sim f(x \\mid Y'_j = y'_j)\\\\ \n",
    "Y'_{j + 1} \\sim f(y \\mid X'_j = x'_j)\n",
    "$$\n",
    "\n",
    "The idea is that it is generally easier to specify these conditionals analytically than the marginal.\n",
    "\n",
    "It turns out that as $k \\rightarrow \\infty$, $X'_k = x'_k$ is effectively a sample from $f(x)$. This is what the Gibbs sampler does.\n",
    "\n",
    "### Applied\n",
    "\n",
    "In the present case: $f(x \\mid y)$ is $\\textrm{Binomial}(n, y)$, while $f(y \\mid x)$ is $\\textrm{Beta}(x + \\alpha, n - x + \\beta)$. Lets wrap the `scipy` rvs interface with our parameters to get some samplers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def x_given_y(y):\n",
    "    return sp.stats.binom(n=N, p=y).rvs()\n",
    "    \n",
    "def y_given_x(x):\n",
    "    return sp.stats.beta(x + ALPHA, N - x + BETA).rvs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sampler implements the relationships above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def GibbsSampler(k, initial_state):\n",
    "    y_0 = initial_state\n",
    "    \n",
    "    y_current = y_0\n",
    "    x_current = 0\n",
    "    \n",
    "    for _ in range(k):\n",
    "        # implements relations above\n",
    "        x_current = x_given_y(y_current)\n",
    "        y_next = y_given_x(x_current)\n",
    "    \n",
    "        # update for next iteration\n",
    "        y_current = y_next\n",
    "    \n",
    "    return x_current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we generate 500 samples each from the analytic solution (the Beta-Binomial) and the Gibbs-sampled marginal (with k iterations), respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analytic_samples = BetaBinomialSampler(alpha=ALPHA, beta=BETA, n=N, size=500)\n",
    "\n",
    "gibbs_initial = 0.01\n",
    "gibbs_samples = [GibbsSampler(10, initial_state=gibbs_initial) for _ in range(500)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEECAYAAAAifS8cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGrxJREFUeJzt3X9sG/Xh//GX7TTNJ80RJzGpKFE+oYUJ1Sp0NIXyY21K\nh0B8kYgmzRPN2EBMX422UKX7kaJ+qhStQMOvJFrXgvpFwNZqX4JEIqi2MX1EgjYQWjPKFtxlrKL0\nAxSaX03iJE3T2Pf5o0rWH7Hv4p5j5/p8/JPYd/e+lx3nlcv5fOcxTdMUAMB1vOkOAABIDQoeAFyK\nggcAl6LgAcClKHgAcCkKHgBcKstqhmPHjqmhoUEej0emaer48eP63ve+p5UrV6qhoUHd3d0qLi5W\ndXW1cnNzZyIzAMAGyy34BQsW6Omnn1ZdXZ127NihnJwc3XjjjWppadGSJUvU2NioYDCo5uZmWysM\nh8MXHXomkNM5syGjRE6nkdNZyeSc1i6ajo4OzZ8/X4FAQO3t7Vq1apUkqaKiQgcOHEhZyHQgp3Nm\nQ0aJnE4jp7NSXvDvv/++brvtNknSwMCA/H6/JMnv92twcHDaKwcApI7tgh8fH1d7e7tWrFiRyjwA\nAId47J6Lpr29XW+//ba2bNkiSaqurlZtba38fr/6+/v1+OOPq76+/oLlwuHwOf9ahEIhh6IDwKWl\nqalp8vtgMKhgMJhwftsF39DQoKVLl6qiokKStHfvXuXl5amyslItLS0aHh5WVVWVrZDHjh2zNV86\nGYahSCSS7hiWZkPO2ZBRIqfTyOmsBQsWTHsZW7toxsbG1NHRoZtuumnyvsrKSnV0dGjjxo3q6OhQ\nZWXltFcOAEgd21vwTmIL3jmzIedsyCiR02nkdFbKtuABALMPBQ8ALkXBA4BLUfAA4FIUPAC4lOXZ\nJDE7+E70SH3d8WcovFzRgsDMBQKQdhS8W/R1a2xHTdzJ2ZvrJAoeuKSwiwYAXIqCBwCXouABwKUo\neABwKQoeAFyKggcAl6LgAcClKHgAcCkKHgBcioIHAJei4AHApSh4AHApCh4AXIqCBwCXouABwKUo\neABwKVsX/BgZGdELL7ygzz//XB6PRw8//LCuuOIKNTQ0qLu7W8XFxaqurlZubm6q8wIAbLJV8C+/\n/LK++c1vatOmTYpGozp16pTeeOMNLVmyRPfee69aWlrU3NysqqqqVOfNOD2jMfUMjyecJzAvS4Ec\n/lkCMLMsC/7kyZPq7OzU+vXrJUk+n0+5ublqb2/Xtm3bJEkVFRXatm3bpVnww+Oq+cOnCeepu2uh\nAjnZM5QIAM6wLPjjx4/LMAzt2rVLR48e1cKFC/XAAw9oYGBAfr9fkuT3+zU4OJjysAAA+ywLPhaL\n6ciRI3rooYe0aNEivfLKK2ppabG9gnA4rHA4PHk7FArJMIzk0s6g7OxsWzl9/f3W8/h8KXvMEzlP\n+RL/KH2+LOWm6Xm3+1ymGzmdRU7nNTU1TX4fDAYVDAYTzm9Z8IWFhSoqKtKiRYskSStWrFBLS4v8\nfr/6+/snv+bn50+5/FQhIpGI5QNJN8MwbOWMRqO25knVY57I6Ysmfh8gGh1P2/Nu97lMN3I6i5zO\nMgxDoVBoWstYvvPn9/tVVFSkY8eOSZI6OjpUUlKiZcuWqa2tTZLU1tam8vLy6ScGAKSMraNoHnzw\nQf3yl7/U+Pi45s+fr3Xr1ikWi6m+vl6tra0KBALatGlTqrMCAKbBVsGXlZXpqaeeuuD+rVu3Oh4I\nAOAMDs4GAJei4AHApSh4AHApCh4AXIqCBwCXouABwKUoeABwKQoeAFzK1ged4CzfiR6przv+DIWX\nK1oQmLlAAFyJgk+Hvm6N7aiJOzl7c51EwQO4SOyiAQCXYgv+EtGVN1/dvWNxp3NZQcB9KPhLRPdp\nr2r+O/6lBbmsIOA+bLIBgEtR8ADgUhQ8ALgUBQ8ALkXBA4BLUfAA4FIUPAC4FAUPAC5FwQOAS1Hw\nAOBStk5VsH79euXm5srj8cjn8+mpp57S0NCQGhoa1N3dreLiYlVXVys3NzfVeQEANtkqeI/Ho9ra\nWuXl5U3e19LSoiVLlujee+9VS0uLmpubVVVVlbKgmP0sz4MvcS58wEG2Ct40TZmmec597e3t2rZt\nmySpoqJC27Zto+CRmMV58CXOhQ84yfYW/BNPPCGPx6Nvf/vbWrNmjQYGBuT3+yVJfr9fg4ODKQ0K\nAJgeWwW/ffv2yRLfvn27FixYYHsF4XBY4XB48nYoFJJhGNNPOsOys7Nt5fT191vP4/OdM9YpX+Kn\n3efLUq7N52gip9WYHo9nWhkl6fRXXyjW0xV3GW+gWHOuKHEs45kc9h+70+z+zNONnM6aLTklqamp\nafL7YDCoYDCYcH5bBT+xpX7ZZZdp+fLlOnz4sPx+v/r7+ye/5ufnT7nsVCEikYid1aaVYRi2ckaj\nUVvznD2WLzpuMf+47edoIqfVmOfvYrPKKEm+48csLy04mjf1zz2ZjGdy2H/sTrP7M083cjprNuUM\nhULTWsbyMMlTp05pdHRUkjQ6Oqq///3vKi0t1bJly9TW1iZJamtrU3l5+fQTAwBSxnILfmBgQM88\n84w8Ho+i0ai+9a1v6frrr9eiRYtUX1+v1tZWBQIBbdq0aSbyAgBssiz44uJiPfPMMxfcn5eXp61b\nt6YkFGYHrvMKZDauyYqkcZ1XILOxeQUALsUWPDIKu30A51DwyCjs9gGcw6YQALgUBQ8ALkXBA4BL\nUfAA4FIUPAC4FAUPAC5FwQOAS3EcfAbiwz4AnEDBZyA+7APACWwGAoBLUfAA4FLsooHr9YzG1DMc\n/3KBvKcBt6Lg4Xo9w+Oq+QPvaeDSw2YLALgUW/DAeXwneqS+7sQzFV6uaEFgZgIBSaLggfP1dWts\nR03CWbI310kUPDIcu2gAwKUoeABwKQoeAFzK9j74WCymxx57TIWFhaqpqVFXV5caGxs1NDSkq666\nSo888oh8Pl8qswIApsF2wf/ud7/TlVdeqZMnT0qS9u3bp3vuuUc333yz9uzZo3feeUd33HFHyoLa\nYXn0A0c+ALiE2Cr43t5eHTx4UN/5zne0f/9+SdLHH3+sjRs3SpJWrVql119/Pe0Fb3X0A0c+ALiU\n2NoH/+qrr+r++++Xx+ORJEUiEeXl5cnrPbN4UVGRTpw4kbqUAIBps9yC//DDD5Wfn6+ysjKFw2FJ\nkmmaMk3znPkmyv984XB4cjlJCoVCMgzjYjLHdcqX+OH0XHaFevqjCecpzstWSUGusrOzbeX09fdb\nz+PznTOWVc54z+VU403kdHJMp3PazZhsTitWP6Nkclq9liZeR6lg97WZbuR0XlNT0+T3wWBQwWAw\n4fyWr+TOzk61t7fr4MGDGhsb08mTJ/XKK69oZGREsVhMXq9Xvb29KigomHL5qUJEIhE7j2XafNH4\nJ5SSpOOnpJo//ivhPHV3LVR+VlSGYdjKGY0m/oMxMc/ZY1nlPP+PZ6LxJnI6OabTOe1mTDanFauf\nUTI5rV5LE6+jVLD72kw3cjrLMAyFQqFpLWNZ8GvXrtXatWslSYcOHdJbb72lRx99VPX19frggw90\nyy236N1331V5eXlyqQEAKZH0cfBVVVXav3+/Nm7cqKGhId1+++1O5gIAXKRpnYtm8eLFWrx4sSSp\nuLhYTz75ZEpCAQAuHp9kBQCXytizSfaOmho5HYs7fd4crwpzEh9xAQCXsowt+M6ek3r63f+JO33L\n6v/UjQv+I+U5LD8dm1ea8gwAkIyMLfiMYXVu8Npfz1wWAJgG9sEDgEtR8ADgUhQ8ALgUBQ8ALkXB\nA4BLUfAA4FIUPAC4FAUPAC5FwQOAS1HwAOBSFDwAuBQFDwAuRcEDgEtR8ADgUhQ8ALgUBQ8ALkXB\nA4BLcUUnzGqWl1SUuKwiLlkUPGY3q0sqSlxWEZcsy4I/ffq0amtrNT4+rmg0qhUrVui73/2uurq6\n1NjYqKGhIV111VV65JFH5PP5ZiIzAMAGy4KfM2eOamtrNXfuXMViMW3dulVLly7V/v37dc899+jm\nm2/Wnj179M477+iOO+6YicwAABtsvck6d+5cSWe25qPRqDwej8LhsG666SZJ0qpVq/SXv/wldSkB\nANNmax98LBbT5s2bdfz4cd15552aP3++5s2bJ6/3zN+HoqIinThxIqVBAQDTY6vgvV6vnn76aY2M\njOjZZ5/Vl19+ecE8Ho9nymXD4bDC4fDk7VAoJMMwrNfpGbXMdP44p3yJH068jGfz+XwyDEPZ2dky\nDMPRMZ3KefZ4TuU8P6OTOe1mTCank2OmMqeTJnJmOnI6r6mpafL7YDCoYDCYcP5pHUWTm5urxYsX\n65NPPtHw8LBisZi8Xq96e3tVUFAw5TJThYhEIpbripmxxNNjsQvG8UXHEy5jmqbleqPRqCKRiAzD\nUCQScXRMp3KePZ5TOc/P6GROuxmTyenkmKnM6aSJnJmOnM4yDEOhUGhay1jugx8cHNTIyIgkaWxs\nTB0dHSopKVEwGNQHH3wgSXr33XdVXl6eRGQAQKpYbsH39/frV7/6lWKxmEzT1C233KIbbrhBJSUl\namho0GuvvaaysjLdfvvttlfqGT8t7+FD0ujJ+DMVJv7XAwCQmGXBl5aWqq6u7oL7i4uL9eSTTya1\nUk8sqvH///8U+/Jo/Hn+a19SYwMAzuBcNADgUhQ8ALgU56IBZoCdk6J1Bf5T3ePxfyUD87IUyGGb\nDPZR8MBMsHFStO7aX6vm3f+JO73uroUK5GQ7nQwuxuYAALgUBQ8ALkXBA4BLUfAA4FK8yQpgkuXR\nPoWXK1oQmLlAuCgUPIB/szjaJ3tznUTBzxrsogEAl6LgAcClKHgAcCkKHgBcioIHAJei4AHApSh4\nAHApCh4AXIqCBwCXouABwKUoeABwKQoeAFyKggcAl7I8m2Rvb6927typ/v5+eb1erVmzRnfffbeG\nhobU0NCg7u5uFRcXq7q6Wrm5uTORGQBgg2XB+3w+/fCHP1RZWZlGR0dVU1Oj66+/Xq2trVqyZInu\nvfdetbS0qLm5WVVVVTORGQBgg+UuGr/fr7KyMklSTk6OrrzySvX29qq9vV2rVq2SJFVUVOjAgQMp\nDQoAmJ5p7YPv6urS0aNH9Y1vfEMDAwPy+/2SzvwRGBwcTElAAEBybF/RaXR0VM8//7weeOAB5eTk\n2F5BOBxWOByevB0KhZSTk6NhqwU9iSd7vV4ZhnHOfad8iR+Ox2MxqM7skjIMQ9nZ2TIMw9Exncp5\n9nhO5Tw/o5M57WZMJqeTY862nHZ8cWJEXUNjcacX52WrpODf751Z5fT5spRrc93TyZlOsyWnJDU1\nNU1+HwwGFQwGE85vq+Cj0aiee+45rVy5UsuXL5d0Zqu9v79/8mt+fv6Uy04VYnR01HqlZuLJsVhM\nkUjknPt80fHEQ5oWg+rMY41EIjIMQ5FIxNExncp59nhO5Tw/o5M57WZMJqeTY862nHZ8NTCmmj98\nGnd63V0LlZ8VtZ0zGh23ve7p5Eyn2ZQzFApNaxlbu2h2796tkpIS3X333ZP3LVu2TG1tbZKktrY2\nlZeXT2vFAIDUstyC7+zs1J/+9CeVlpbq5z//uTwej+677z5VVlaqvr5era2tCgQC2rRp00zkBQDY\nZFnw1157rV577bUpp23dutXxQAAAZ/BJVgBwKQoeAFyKggcAl6LgAcClKHgAcCnbn2QFgGT4TvRI\nfd3xZyi8XNGCwMwFuoRQ8ABSq69bYztq4k7O3lwnUfApwS4aAHApCh4AXIqCBwCXouABwKUoeABw\nKQoeAFyKggcAl6LgAcClKHgAcCkKHgBcilMVALOU5TleJCmv1NF1duXNV3fvWMJ5AvOyFMhh2zET\nUPDAbGVxjhdJUu2vHV1l92mvav7704Tz1N21UIGcbEfXi+TwZxYAXIqCBwCXouABwKUoeABwKcs3\nWXfv3q0PP/xQ+fn5evbZZyVJQ0NDamhoUHd3t4qLi1VdXa3c3NyUhwUA2Ge5Bb969Wpt2bLlnPta\nWlq0ZMkSNTY2KhgMqrm5OWUBAQDJsSz4a6+9VvPmzTvnvvb2dq1atUqSVFFRoQMHDqQmHQAgaUnt\ngx8YGJDf75ck+f1+DQ4OOhoKAHDxUv5Bp3A4rHA4PHk7FAopJydHw1YLehJP9nq9MgzjnPtO+RI/\nHI/HYlBJPp9PhmEoOztbhmE4OqZTOc8ez6mc52d0MqfdjMnkdHJMcqbnd8jny1Luea+9mTSRczZo\namqa/D4YDCoYDCacP6mC9/v96u/vn/yan58fd96pQoyOjlqvxEw8ORaLKRKJnHOfLzqeeEjTYlBJ\n0WhUkUhEhmEoEok4OqZTOc8ez6mc52d0MqfdjMnkdHJMcqbndygaHb/gtTeTJnJmOsMwFAqFprWM\nrV00pmme84NdtmyZ2traJEltbW0qLy+f1koBAKlnuQXf2NioQ4cOKRKJ6OGHH1YoFFJlZaXq6+vV\n2tqqQCCgTZs2zURWAMA0WBb8xo0bp7x/69atjocBADiHT7ICgEtR8ADgUhQ8ALgUBQ8ALsUVnQCk\nVSouA9gzGlPPcOLj7y+FSwtS8ADSKhWXAewZHlfNH7i0oLv/fAHAJYyCBwCXYhcNgFnHd6JH6uuO\nP0Ne6cyFyWAUPIDZp69bYztq4k+v/fXMZclg7KIBAJei4AHApSh4AHApCh4AXIo3WQHAhtn46VgK\nHgBsmI2fjs2cPzUAAEexBQ8AcueHpyh4AJBc+eEpdtEAgEtR8ADgUhQ8ALjURe2D/+ijj/TKK6/I\nNE2tXr1alZWVTuUCAFykpLfgY7GYXnrpJW3ZskXPPfec3nvvPX355ZdOZgMAXISkC/7w4cO64oor\ndPnllysrK0u33nqrDhw44GQ2AMBFSLrg+/r6VFRUNHm7sLBQfX19joQCAFw8R4+D93g8Nmf0Kuv/\nhGQODcad5aqi/9D/vfGKuNNL8udONx4AXFI8pmmaySz4ySef6PXXX9eWLVskSS0tLZJ0wRut4XBY\n4XB48nYoFEo2KwBc0pqamia/DwaDCgaDiRcwkxSNRs0NGzaYXV1d5unTp82f/vSn5ueff2653Guv\nvZbsKmcUOZ0zGzKaJjmdRk5nJZMz6V00Xq9XDz30kLZv3y7TNHX77berpKQk2eEAAA67qH3wS5cu\nVWNjo1NZAAAO8m3btm3bTK+0uLh4pleZFHI6ZzZklMjpNHI6a7o5k36TFQCQ2TgXDQC4FAUPAC41\nYxf8mA0nJuvt7dXOnTvV398vr9erNWvW6O677053rLhisZgee+wxFRYWqqYmwYUK0mhkZEQvvPCC\nPv/8c3k8Hj388MO65ppr0h3rAvv371dra6s8Ho9KS0u1bt06ZWWl/3o4u3fv1ocffqj8/Hw9++yz\nkqShoSE1NDSou7tbxcXFqq6uVm5ubsbl3Lt3r/76178qKytL8+fP17p16zIy54Q333xT+/bt00sv\nvaS8vLw0JYyf8fe//73efvtt+Xw+3XDDDaqqqrIezOFDNac01THzX3zxxUyselpOnDhhHjlyxDRN\n0zx58qT56KOPZmTOCW+99ZbZ2Nho7tixI91R4tq5c6f5zjvvmKZpmuPj4+bw8HCaE12ot7fXXL9+\nvXn69GnTNE3z+eefN9va2tKc6ox//OMf5pEjR8yf/OQnk/f95je/MVtaWkzTNM3m5mZz79696Yo3\naaqcf/vb38xoNGqapmnu3bvX3LdvX7riTZoqp2maZk9Pj7l9+3Zz3bp1ZiQSSVO6M6bK+PHHH5u/\n+MUvzPHxcdM0TXNgYMDWWDOyi2a2nJjM7/errKxMkpSTk6Mrr7wyY8+v09vbq4MHD2rNmjXpjhLX\nyZMn1dnZqdWrV0uSfD5f2rfg4onFYhodHVU0GtWpU6dUUFCQ7kiSpGuvvVbz5s0757729natWrVK\nklRRUZERv0tT5bzuuuvk9Z6pmGuuuUa9vb3piHaOqXJK0quvvqr7778/DYkuNFXGP/7xj6qsrJTP\n55MkXXbZZbbGmpH/Qac6Mdnhw4dnYtVJ6+rq0tGjRzNyd4L07xfkyMhIuqPEdfz4cRmGoV27duno\n0aNauHChHnzwQWVnZ6c72jkKCwt1zz33aN26dZo7d66uu+46XXfddemOFdfAwID8fr+kMxslg4Px\nz+mUKVpbW3XrrbemO8aU2tvbVVRUpNLSzL2o9ldffaVDhw7pt7/9rbKzs/X9739fixYtslwubW+y\n2j4xWRqMjo7q+eef1wMPPKCcnJx0x7nAxP65srIymaYpM0OPdI3FYjpy5IjuvPNO1dXVae7cuZPn\nLMokw8PDam9v165du/Tiiy9qdHRUf/7zn9MdyzXeeOMN+Xw+3XbbbemOcoGxsTE1Nzefc46sTPx9\nikajGhkZ0RNPPKGqqirV19fbWm5GCr6wsFA9PT2Tt/v6+jLmX+DzRaNRPffcc1q5cqWWL1+e7jhT\n6uzsVHt7uzZs2KDGxkaFw2Ht3Lkz3bEuUFhYqKKiosktjRUrVujTTz9Nc6oLdXR0qLi4WHl5efJ6\nvbrpppv0z3/+M92x4vL7/erv75ck9ff3Kz8/P82J4mtra9PBgwe1cePGdEeZ0tdff62uri797Gc/\n0/r169XX16fNmzdrYGAg3dHOEQgEdOONN0qSrr76ank8HkUiEcvlZmQXzdVXX62vv/5a3d3dKigo\n0HvvvZexP/Ddu3erpKQko4+eWbt2rdauXStJOnTokN566y1t2LAhzaku5Pf7VVRUpGPHjmnBggXq\n6OjIyPMVBQIB/etf/9LY2JjmzJmjjo4OW//+zpTz/0tbtmyZ2traVFlZqba2NpWXl6cx3b+dn/Oj\njz7Sm2++qccff1xz5sxJY7JznZ2ztLRUe/bsmZy2fv161dXVpfUoGunC53L58uX6+OOPtXjxYh07\ndkzRaFSGYViOM2OfZP3oo4/08ssvT56YLBMPk+zs7FRtba1KS0vl8Xjk8Xh03333aenSpemOFtdE\nwWfqYZKfffaZXnzxRY2Pj2fMoXJTef311/X+++/L5/OprKxMP/7xjzPiMMnGxkYdOnRIkUhE+fn5\nCoVCWr58uerr69XT06NAIKBNmzZN+cZhunM2NzdrfHx8soiuueYa/ehHP8q4nBMHAUjShg0btGPH\njrQW/FQZV65cqV27dumzzz7TnDlz9IMf/ECLFy+2HItTFQCAS/FJVgBwKQoeAFyKggcAl6LgAcCl\nKHgAcCkKHgBcioIHAJei4AHApf4X5LsdQLOJiq0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f39a3795cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist([analytic_samples, gibbs_samples], bins=np.arange(16))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's important to make sure the sampler is settling down on a good estimate of the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def GibbsSamplerWithHistory(k, initial_state):\n",
    "    y_0 = initial_state\n",
    "    \n",
    "    y_current = y_0\n",
    "    x_current = 0\n",
    "    \n",
    "    y_history = []\n",
    "    x_history = []\n",
    "    \n",
    "    for _ in range(k):\n",
    "        x_history.append(x_current)\n",
    "        y_history.append(y_current)\n",
    "        \n",
    "        # implements relations above\n",
    "        x_current = x_given_y(y_current)\n",
    "        y_next = y_given_x(x_current)\n",
    "    \n",
    "        # update for next iteration\n",
    "        y_current = y_next\n",
    "    \n",
    "    return x_history, y_history"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a bunch of draws, and study the behavior of the sampler by looking at the entropy of the distribution generated if we stopped the sampler in its tracks at the kth iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "MAX_K = 10\n",
    "\n",
    "histories = []\n",
    "\n",
    "for _ in range(500):\n",
    "    x_history, y_history = GibbsSamplerWithHistory(k=MAX_K, initial_state=0.187)\n",
    "    histories.append((x_history, y_history))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dist_at_iteration(k):\n",
    "    dist = []\n",
    "    for history in histories:\n",
    "        dist.append(history[0][k])\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize(d):\n",
    "    raw = sum(d.values())\n",
    "    return {key:value/raw for key,value in d.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "\n",
    "ks = []\n",
    "entropies = []\n",
    "\n",
    "for k in range(MAX_K):\n",
    "    c = Counter(dist_at_iteration(k))\n",
    "    H = sp.stats.entropy(list(normalize(c).values()))\n",
    "    entropies.append(H)\n",
    "    ks.append(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEECAYAAADK0VhyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF2tJREFUeJzt3WtsU+cdx/HfiXMriZXE4SKSKPMImShmdCWhi0pVFkCq\nBpVgquqNTQi2Tu2AlsI2haFqoGkXtdy60JRs641p6RtvWtE69gaJizrWTqEkE3KgJYLA2ggI8QIG\nkmHHz15EeKQJl/Sc2LHz/UhVfezHfv5/QP75POf42DLGGAEAxrWMZBcAAEg+wgAAQBgAAAgDAIAI\nAwCACAMAgKRMuy8QiUS0ZcsWRaNR9ff3q6amRk8++eSgMdFoVA0NDTp9+rTcbrc2bNigiRMn2p0a\nAOAQ23sGWVlZ2rJli7Zu3apt27aptbVV7e3tg8YcOHBA+fn52rVrl5YsWaKmpqZ7eu1gMGi3vDGN\n/lIb/aWudO5N+nz9ObJMlJOTI2lgL6G/v3/I483NzZo/f74kqaamRsePH7+n1+UvLLXRX2pL5/7S\nuTfp8/Vne5lIkmKxmH7yk5/owoULeuyxxzR9+vRBj4dCIRUXF0uSMjIylJeXp6tXryo/P9+J6QEA\nNjmyZ5CRkaGtW7eqsbFRp06d0ieffHLH8VwBAwDGFsvpaxP96U9/Um5urh5//PH4fb/61a/05JNP\nqrKyUrFYTE8//bRef/31Ic8NBoODdm/8fr+TpQHAuBEIBOK3fT6ffD7fHcfbXia6cuWKMjMzNWHC\nBN24cUPHjx/X0qVLB42pqqrS4cOHVVlZqffff1+zZs0a9rWGK7izs9NuiWOW2+1WOBxOdhmjhv5S\nWzr3l869SVJJScmIP0zbDoOenh69+uqrisViMsbo4Ycf1pw5cxQIBFRRUaGqqiotWLBAr7zyitat\nWye3263nn3/e7rQAAAc5vkzkNPYMUhf9pbZ07i+de5MG9gxGim8gAwAIAwAAYQAAEGEAABBhAAAQ\nYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBh\nAAAQYQAAEGEAABBhAACQlGn3Bbq7u9XQ0KCenh5lZGRo4cKFWrx48aAxbW1t2rp1q6ZMmSJJeuih\nh/TEE0/YnRoA4BDbYeByubRy5Up5vV719fVp48aNeuCBB1RaWjpo3P3336+NGzfanQ4AUkYsZunE\nifvU0eGS19uvmTN7ZVkm2WUNy3YYFBYWqrCwUJKUm5ur0tJShUKhIWFgzNj8AwAwftx8cz53LkPl\n5RNG/c35xIn7tGRJgSIRS1lZRvv2ST7f9VGbzw7bYXCrixcv6uzZs6qsrBzy2KlTp1RXV6eioiKt\nWLFCZWVlTk4NAHeV6Dfnjg6XIhFLkhSJWOrocMnnG7XpbHEsDPr6+rRz506tWrVKubm5gx6bNm2a\ndu/erZycHLW0tGjbtm2qr68f8hrBYFDBYDC+7ff75Xa7nSpxzMnOzqa/FEZ/zohGjVpaLJ05k6Ev\nfjGmOXOMXC5rVOY6dy5j0JvzuXNZqqkZvR4rKqSsLBMPn4oKJezfTCAQiN/2+Xzy3SWFLOPA+k1/\nf79efPFFPfjgg0MOHg9n7dq1eumll5Sfn3/XsZ2dnXbLG7PcbrfC4XCyyxg19Jea/r+UkqXy8sio\nL6UEgxM+82n98qh9Wk/kXJJkjKW2tsQfMygpKRnxcxzZM2hsbFRZWdltg6Cnpyd+XKG9vV2S7ikI\nACReOi+lzJzZq337NCjoRpNlGfl818fs0tCtbIfByZMn9d5776m8vFx1dXWyLEvLly9XV1eXLMvS\nokWL9MEHH2j//v1yuVzKzs7W+vXrnagdGBcSfUZKote5vd7+QUspXm//qM118825psatcHhsHshN\nFkeWiUYTy0Spi/6ckeiljfGwlJLu/zaTtkwEjCeJPj0x0Z/UWUoZnwgDYIQSvaaeyGUUiaWU8Yow\nAEYoWZ/Ub11GAZxGGAAjlKxP6iyjYDQRBkh5iT7bJtFr6kAiEAZIeYlew2dNHemI3zNAyhtuDR/A\nyBAGSHk31/AlJWQNH0hHLBMh5XG2DWAfYYCUx9k2gH2EARyX6G/oArCPMIDjUunXnQAM4AAyHMfZ\nPUDqIQzgOM7uAVIPy0RwHN/QBVIPYQDH8Q1dIPWwTAQAIAwAAIQBAECEAQBAhAEAQIQBAECEAQBA\nDnzPoLu7Ww0NDerp6VFGRoYWLlyoxYsXDxn35ptvqrW1VTk5OVq7dq28Xq/dqQEADrEdBi6XSytX\nrpTX61VfX582btyoBx54QKWlpfExLS0tunDhgnbt2qVTp07ptdde0y9/+Uu7UwMAHGJ7maiwsDD+\nKT83N1elpaUKhUKDxjQ3N2v+/PmSpMrKSl2/fl09PT12pwYAOMTRYwYXL17U2bNnVVlZOej+UCik\n4uLi+LbH4xkSGACA5HHs2kR9fX3auXOnVq1apdzc3LuOtyxryH3BYFDBYDC+7ff75Xa7nSpxzMnO\nzqa/FEZ/qSude7spEAjEb/t8Pvnu8lOAjoRBf3+/duzYoUcffVRz584d8rjH41F3d3d8u7u7W0VF\nRUPGDVdwOBx2osQxye12018Ko7/Ulc69SQP9+f3+ET3HkWWixsZGlZWVDXsWkSRVV1fr8OHDkqSP\nP/5YeXl5KiwsdGJqAIADbO8ZnDx5Uu+9957Ky8tVV1cny7K0fPlydXV1ybIsLVq0SHPmzFFLS4ue\ne+455ebmavXq1U7UDgBwiGWMGdO/VN7Z2ZnsEkbNeNhVpb/Ulc79pXNvklRSUjLi5/ANZAAAYQAA\nIAwAAOI3kMeFWMzSiRP3qaPDJa+3XzNn9sqyxvShIgAJRhiMAydO3KclSwoUiVjKyjLat0/y+fih\negD/xzLRONDR4VIkMvCN70jEUkeHK8kVARhrCINxwOvtV1bWwLJQVpaR19uf5IoAjDUsE40DM2f2\nat8+DTpmAAC3IgzGAcsy8vmu6y7XqQIwjrFMBAAgDAAAhAEAQIQBAECEAQBAhAEAQIQBAECEAQBA\nhAEAQIQBAECEAQBAhAEAQIQBAECEAQBADl3CurGxUceOHVNBQYG2b98+5PG2tjZt3bpVU6ZMkSQ9\n9NBDeuKJJ5yYGgDgAEfCoLa2Vl//+tfV0NBw2zH333+/Nm7c6MR0AACHObJMNGPGDOXl5d1xjDHG\niakAAKMgYb90durUKdXV1amoqEgrVqxQWVlZoqYGANxFQsJg2rRp2r17t3JyctTS0qJt27apvr5+\nyLhgMKhgMBjf9vv9crvdiSgxKbKzs+kvhdFf6krn3m4KBALx2z6fT767/O5tQsIgNzc3fvvBBx/U\n66+/rqtXryo/P3/QuOEKDofDiSgxKdxuN/2lMPpLXencmzTQn9/vH9FzHDu11Bhz2+MCPT098dvt\n7e2SNCQIAADJ48ieQX19vdra2hQOh7V69Wr5/X5Fo1FZlqVFixbpgw8+0P79++VyuZSdna3169c7\nMS0AwCGWGeOn+XR2dia7hFEzHnZV6S91pXN/6dybJJWUlIz4OXwDGQBAGAAACAMAgAgDAIAIAwCA\nCAMAgAgDAIAIAwCACAMAgAgDAIAIAwCACAMAgAgDAIAIAwCACAMAgAgDAIAIAwCACAMAgAgDAIAI\nAwCACAMAgAgDAIAIAwCApEwnXqSxsVHHjh1TQUGBtm/fPuyYN998U62trcrJydHatWvl9XqdmBoA\n4ABH9gxqa2v1wgsv3PbxlpYWXbhwQbt27dLTTz+t1157zYlpAQAOcSQMZsyYoby8vNs+3tzcrPnz\n50uSKisrdf36dfX09DgxNQDAAQk5ZhAKhVRcXBzf9ng8CoVCiZgaAHAPHDlm8HlYljXkvmAwqGAw\nGN/2+/1yu92JLCuhsrOz6S+F0V/qSufebgoEAvHbPp9PPp/vjuMTEgYej0fd3d3x7e7ubhUVFQ0Z\nN1zB4XB41OtLFrfbTX8pjP5SVzr3Jg305/f7R/Qcx5aJjDEyxgz7WHV1tQ4fPixJ+vjjj5WXl6fC\nwkKnpgYA2OTInkF9fb3a2toUDoe1evVq+f1+RaNRWZalRYsWac6cOWppadFzzz2n3NxcrV692olp\nAQAOscztPs6PEZ2dnckuYdSMh11V+ktd6dxfOvcmSSUlJSN+Dt9ABgAQBgAAwgAAIMIAACDCAAAg\nwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAACDC\nAAAgwgAAIMIAACDCAAAgKdOJF2ltbdWePXtkjFFtba2WLVs26PFDhw6pqalJxcXFkqTHHntMCxYs\ncGJqAIADbIdBLBbTG2+8oc2bN6uoqEibNm3S3LlzVVpaOmjcww8/rO9973t2pwMAjALby0Tt7e2a\nOnWqJk2apMzMTM2bN0/Nzc1O1AYASBDbewahUCi+/CNJHo9H7e3tQ8b985//1IkTJzR16lStXLly\n0HMAAMnlyDGDz7Isa9B2dXW1HnnkEWVmZmr//v169dVXtXnz5iHPCwaDCgaD8W2/3y+32z0aJY4J\n2dnZ9JfC6C91pXNvNwUCgfhtn88nn893x/G2w8Dj8ejSpUvx7VAopKKiokFj8vPz47cXLlyot99+\ne9jXGq7gcDhst8Qxy+12018Ko7/Ulc69SQP9+f3+ET3H9jGD6dOn6/z58+rq6lI0GtWRI0dUXV09\naExPT0/89tGjR1VWVmZ3WgCAg2zvGWRkZOipp57SL37xCxljtGDBApWVlSkQCKiiokJVVVX629/+\npg8//FAul0v5+flas2aNE7UDABxiGWNMsou4k87OzmSXMGrGw64q/aWudO4vnXuTpJKSkhE/h28g\nAwAIAwAAYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAA\nEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAACQlOnEi7S2tmrPnj0yxqi2tlbLli0b9Hg0\nGlVDQ4NOnz4tt9utDRs2aOLEiU5MDQBwgO09g1gspjfeeEMvvPCCduzYoSNHjujTTz8dNObAgQPK\nz8/Xrl27tGTJEjU1NdmdFgDgINth0N7erqlTp2rSpEnKzMzUvHnz1NzcPGhMc3Oz5s+fL0mqqanR\n8ePH7U6b0mIxS8HgBAUCGQoGJ8gYK9klARjnbC8ThUIhFRcXx7c9Ho/a29tvOyYjI0N5eXm6evWq\n8vPz7U6fkk6cuE9LlhQoErGUlWW0b5/k811PdlkAxjFHjhl8lmXd+ZOuMWbY+4PBoILBYHzb7/fL\n7XY7WttYcO5chiKRgT+jSMTSuXNZqqlJvz6zs7PT8u/vJvpLXenc202BQCB+2+fzyefz3XG87TDw\neDy6dOlSfDsUCqmoqGjQmOLiYnV3d8vj8SgWi6m3t3fYvYLhCg6Hw3ZLHHPKyycoK8vE9wzKyyMK\nh9Nvz8Dtdqfl399N9Je60rk3aaA/v98/oufYPmYwffp0nT9/Xl1dXYpGozpy5Iiqq6sHjamqqtLh\nw4clSe+//75mzZpld9qUNnNmr/btu6zXX7+mffsua+bM3mSXBGCcs8zt1mxGoLW1VW+99ZaMMVqw\nYIGWLVumQCCgiooKVVVVKRKJ6JVXXlFHR4fcbreef/55TZ48+Z5eu7Oz0255Y9Z4+HRCf6krnftL\n594kqaSkZMTPcSQMRhNhkLroL7Wlc3/p3Jv0+cKAbyADAAgDAABhAAAQYQAAEGEAABBhAAAQYQAA\nEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQYQAAEGEAABBhAAAQ\nYQAAkJRp58lXr17Vr3/9a3V1dWny5MnasGGDJkyYMGTcN7/5TXm9XhljNHHiRNXV1dmZFgDgMFth\nsHfvXn35y1/W0qVLtXfvXr3zzjv6zne+M2Rcbm6uXnrpJTtTAQBGka1loqNHj2r+/PmSpK997Wtq\nbm4edpwxxs40AIBRZmvP4PLlyyosLJQkFRYW6sqVK8OOi0Qi2rRpk1wul5YuXaq5c+famRYA4LC7\nhsHPf/5zXb58Ob5tjJFlWfrWt751z5M0NjaqsLBQFy9e1M9+9jN94Qtf0OTJkz9fxQAAx901DH76\n05/e9rHCwkL19PTE/19QUHDbcZI0efJk+Xw+nTlzZtgwCAaDCgaD8W2/36+SkpK7NpHK3G53sksY\nVfSX2tK5v3TuTZICgUD8ts/nk8/nu+N4W8cMqqqqdOjQIUnSoUOHVF1dPWTMtWvXFI1GJUlXrlzR\nRx99pLKysmFfz+fzye/3x/+7tZl0RH+pjf5SVzr3Jg30d+t76d2CQLJ5zGDZsmV6+eWXdfDgQU2c\nOFE//OEPJUmnT5/W/v379cwzz+jTTz/V7373O2VkZMgYo2984xsqLS21My0AwGG2wiA/P3/YZaRp\n06bpmWeekSR96Utf0vbt2+1MAwAYZWP6G8j3smuTyugvtdFf6krn3qTP159l+BIAAIx7Y3rPAACQ\nGIQBAMDeAeREaGpq0ocffqjMzExNmTJFa9asGfZieKmmtbVVe/bskTFGtbW1WrZsWbJLckx3d7ca\nGhrU09OjjIwMLVy4UIsXL052WY6KxWLatGmTPB6PNm7cmOxyHHX9+nX95je/0b///W9ZlqXVq1er\nsrIy2WU55q9//asOHjwoy7JUXl6uNWvWKDNzzL8V3lZjY6OOHTumgoKC+Mk693oR0UHMGPevf/3L\n9Pf3G2OMaWpqMm+//XaSK7Kvv7/fPPvss+bixYsmEomYH//4x+aTTz5JdlmO+c9//mPOnDljjDGm\nt7fXrFu3Lq36M8aYd99919TX15sXX3wx2aU4rqGhwRw4cMAYY0w0GjXXrl1LckXO6e7uNmvXrjWR\nSMQYY8zOnTvNoUOHklyVPSdOnDBnzpwxP/rRj+L3/eEPfzB79+41xhjzzjvvmKampru+zphfJpo9\ne7YyMgbKrKysVHd3d5Irsq+9vV1Tp07VpEmTlJmZqXnz5t32In+pqLCwUF6vV9LAFWtLS0sVCoWS\nW5SDuru71dLSooULFya7FMf19vbq5MmTqq2tlSS5XK602BO/VSwWU19fn/r7+/Xf//5XRUVFyS7J\nlhkzZigvL2/Qffd6EdFbpdS+0cGDBzVv3rxkl2FbKBRScXFxfNvj8ai9vT2JFY2eixcv6uzZs2m1\nzPD73/9eK1as0PXr15NdiuMuXLggt9ut3bt36+zZs5o2bZq++93vKjs7O9mlOcLj8ejxxx/XmjVr\nlJOTo9mzZ2v27NnJLstx93oR0VuNiTC408Xwbl7i4s9//rNcLpceeeSRZJU5qizLSnYJjuvr69PO\nnTu1atUq5ebmJrscR9xcm/V6vQoGg2l3efZYLKYzZ87oqaeeUkVFhfbs2aO9e/fK7/cnuzRHXLt2\nTUePHtXu3bs1YcIE7dixQ3//+9/T9n1lJMZEGNzpYnjSwHWPWlpatHnz5gRVNLo8Ho8uXboU3w6F\nQim/q/pZ/f392rFjhx599NG0umT5yZMndfToUbW0tOjGjRvq7e1VQ0ODnn322WSX5giPx6Pi4mJV\nVFRIkmpqarR3794kV+Wc48ePa/LkycrPz5ckffWrX9VHH32UdmFwrxcRvdWYP2bQ2tqqv/zlL6qr\nq1NWVlayy3HE9OnTdf78eXV1dSkajerIkSPDXuQvlTU2NqqsrCztziL69re/rcbGRjU0NGj9+vWa\nNWtW2gSBNPAmUlxcrM7OTkkDb563u7BkKpo4caJOnTqlGzduyBij48ePp8W10owxg/ZS7+Uiop81\n5r+BvG7dOkWj0fjlZisrK/X9738/yVXZ19raqrfeekvGGC1YsCCtTi09efKktmzZovLyclmWJcuy\ntHz5cn3lK19JdmmOamtr07vvvpt2p5Z2dHTot7/9raLRaFqdzn3TH//4R/3jH/+Qy+WS1+vVD37w\ng5Q+tbS+vl5tbW0Kh8MqKCiQ3+/X3Llz9fLLL+vSpUvxi4h+9iDzZ435MAAAjL4xv0wEABh9hAEA\ngDAAABAGAAARBgAAEQYAABEGAAARBgAASf8DMErotawyBawAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f39a2e297b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(ks, entropies)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}