Lecture16_04 - Gibbs with conditional a conjugate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gibbs with conditional a conjugate\n",
    "\n",
    "##### Keywords: conjugate prior, gibbs sampler, MCMC,  binomial, beta, beta-binomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from scipy.stats import norm, gamma\n",
    "from scipy.stats import distributions\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import time\n",
    "sns.set_style('whitegrid')\n",
    "sns.set_context('poster')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "Gibbs sampler for a model in which  one of the conditionals look like a distribution which is part of a conjugate pair. In this case we can use Bayes theorem to get the other conditional by multiplying the known conditional by a marginal which is the other part of the conjugate pair. Our example involves a $Binom$ conditional. Multiplying by a $Beta$ marginal leaves us with the other conditional as another $Beta$.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "{:.no_toc}\n",
    "* \n",
    "{: toc}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now going to take a look at a slightly more complicated case that was originally outlined in full generality by Casella and George (1992). Suppose we have a nasty looking joint distribution given as: \n",
    "\n",
    "$$p(x,y) = \\binom{16}{y} x^{y+1} (1-x)^{19-y}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Looks like a binomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For such a situation the *two* conditional distributions are not exactly obvious. Clearly we have a binomial term staring at us, so we should be looking to try and express part of the function as a binomial of the form, \n",
    "\n",
    "$$p(\\theta \\vert \\pi ) = \\binom{n}{\\theta} \\pi^{\\theta} (1-\\pi)^{n-\\theta}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It follows directly that for our example we have a binomial with $n=16$ and $\\theta =y$, \n",
    "\n",
    "$$p( y \\vert x ) = g(x) \\binom{16}{y} x^{y+1} (1-x)^{16-y} .$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The  $ x\\vert y$ conditional\n",
    "\n",
    "So, now we need the conditional for x|y, and we know from Bayes' theorem that :\n",
    "\n",
    "$$p(x \\vert y) = \\frac{p(y \\vert x)p(x)}{p(y)} $$\n",
    "\n",
    "so what we should be looking for is a conjugate prior to a Binomial distribution, which is of course a Beta distibution:\n",
    "\n",
    "$$Beta(\\alpha,\\beta) = x^{\\alpha-1}(1-x)^{\\beta-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this intuition in mind, the math is now trivial:\n",
    "\n",
    "$$p(x \\vert y) = h(y) x^{\\alpha + y - 1}(1-x)^{\\beta + n - y -1}$$\n",
    "\n",
    "which for our example question is simply:\n",
    "\n",
    "$$p(x \\vert y) \\sim Beta(y+\\alpha,n-y+\\beta)$$\n",
    "\n",
    "with $\\alpha=2$ and $\\beta=4$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The sampler\n",
    "\n",
    "With our conditionals formulated, we can move directly to our Gibbs sampler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import binom, beta\n",
    "n=16\n",
    "alph=2.\n",
    "bet=4.\n",
    "\n",
    "\n",
    "def gibbs(N=10000,thin=50):\n",
    "    x=1\n",
    "    y=1\n",
    "    samples=np.zeros((N,2))\n",
    "    for i in range(N):\n",
    "        for j in range(thin):\n",
    "            y=binom.rvs(n,x)\n",
    "            newalph=y+alph\n",
    "            newbet=n-y+bet\n",
    "            x=beta.rvs(newalph, newbet)\n",
    "          \n",
    "        samples[i,0]=x\n",
    "        samples[i,1]=y\n",
    "    return samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAIbCAYAAABrDGF/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuU1WW9P/DPngvXAREJ0UTAS4IKJGQWammmZssunmWm\nQsc4onaxTtpR0fwVrOOtND0palxSsTqdZXq0Trd1AutYR8oTlK1q1CTxDApqIggMMLf9+6M1c0Bm\nhpln75nhgddrLdfC73d/vs+z9zyz93ue+c7zFIrFYjEAAIBsVfR1BwAAgNII9QAAkDmhHgAAMifU\nAwBA5oR6AADInFAPAACZ63aoX7p0aUyZMmWn4z/84Q/jAx/4QEyaNClOP/30+Na3vlWWDgIAAJ2r\n6s6DV6xYEVdeeeVOx3/0ox/FP/3TP8WsWbPihBNOiGXLlsV1110XNTU18eEPf7hsnQUAAHbWpVDf\n0NAQixcvjttvvz0GDRoUjY2NO5y/+eabY/r06fH5z38+IiKOO+64eOGFF+Lxxx8X6gEAoId1KdQ/\n9thjsWjRopg9e3asW7cu7r333rZzf/jDH2Lt2rXx0Y9+dIeam2++ubw9BQAA2tWle+onTZoUS5cu\njenTp0ehUNjh3NNPPx0REY2NjfGxj30sjj766DjppJPiO9/5Tvl7CwAA7KRLoX7kyJFRU1PT7rl1\n69ZFRUVFfOpTn4oTTzwxvvGNb8Spp54ac+fOjR//+Mdl7SwAALCzbv2hbHuampqiubk5PvrRj8bF\nF18cEX+7p76uri7uvPPOOOOMM0ruJAAA0LGS16kfNGhQFAqFOOGEE3Y4Pm3atFi1alU0NTWV2gQA\nwC4Vi8Wk/3JqM7W9vmizL55jTn0tt5Jn6seMGRMRsdOKOE1NTVEsFqOions/NyxfvrzULgEAwB5p\n6tSp7R4vOdS/7W1vi379+sVPfvKTHTal+tnPfhYTJ07sdqiP6LizUKra2tqIiJgwYUIf94Q9nbFG\nbzDOdpQ6+/nGRUB25zZLmeEtpc2UsdYXzzFVLl+Pzia/Sw71NTU1cckll8Sdd94ZgwcPjre//e3x\nwx/+MJYvXx4LFiwo9fIAAMAuJIX6N/5k8alPfSqGDh0a3/zmN+Oee+6JsWPHxh133LHTffYAAED5\ndTvUX3rppXHppZfudHzGjBkxY8aMsnQKAADoupJXvwEAAPqWUA8AAJkT6gEAIHNCPQAAZK7kJS0B\nAHYHpaw3n0ubffUcW9vtjfb74jmm2p36aqYeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAA\nZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQ\nOaEeAAAyJ9QDAEDmhHoAAMhcVV93AADoec3NzX3dhS6rrKxMqkt9jhUV6XOcLS0tSXXFYjGprqoq\nPbo1NTUl17b2tzvXSP06FgqFpLqI9DGQ+nWsrq5OqotIHwMdMVMPAACZE+oBACBzQj0AAGROqAcA\ngMwJ9QAAkDmhHgAAMifUAwBA5oR6AADInM2nAGAvkLoRUE764jnm9LqWsnFV64ZQpVyjN6R+Pfri\n61jKJlvtMVMPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZ\nE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5rod6pcuXRpTpkzp8Pxrr70W06ZNi3nz\n5pXUMQAAoGu6FepXrFgRV155ZaePue666+K1114rqVMAAEDXdSnUNzQ0xMKFC+OCCy6IqqqqDh/3\n6KOPxn//939H//79y9ZBAACgc10K9Y899lgsWrQoZs+eHTNmzGj3MZs2bYq5c+fG7Nmzo7q6uqyd\nBAAAOtbxtPt2Jk2aFEuXLo2ampoO75W/6aab4vDDD48Pf/jDcf3115e1kwDkp1gsJtUVCoUy96Tn\nNDQ0JNVVVlYmt7n969r676amph5rs6WlJamuoiJ9LY5t27Yl1aU+x9SxGpH++nR250NPtFeu2q1b\nt3a5ppQxkCp1DKT2tZT3q1LGXXu6NKJGjhzZ6flly5bFj370o/jBD35Qlk4BAABdl/Zj4na2bt0a\nX/ziF+Ozn/1sHHjggeXoU9TW1pblOvBGW7ZsiQhjjJ5nrO0dM/WpM5/leo6ts6bPPPNMWa63u+jt\nsVPuGdOuyKmvEf/325OVK1f2Sftd1dvvH7vTTH3Jvxe59dZbY8iQIXH++edHc3Nz268Ai8ViNDc3\nl9xBAACgcyXP1C9ZsiTWrFkTkyZNajtWKBTizjvvjLvuuitplmrChAmldgva1ToejTF6mrG2d8zU\n9/U99a0z9G95y1t6rE331Hdub7mnvnWG/tBDD+1yjXvqO5cy7lasWNHhuZJD/fz583d6U/vYxz4W\nZ555Zpx77rmlXh4AANiFkkP94YcfvtOxysrKGDlyZBx55JGlXh4AANiFpN817OpXDYVCIatfnwIA\nQM66PVN/6aWXxqWXXtrpY5544onkDgEAAN3T+3/BAAAAlFXJ99QDQHv2htsw+/Xr16ftt77GqSup\ndEUpK/WkGjBgQK+3SedaV4fxtSmfcr9HmqkHAIDMCfUAAJA5oR4AADIn1AMAQOaEegAAyJxQDwAA\nmRPqAQAgc0I9AABkTqgHAIDMCfUAAJA5oR4AADIn1AMAQOaEegAAyJxQDwAAmRPqAQAgc0I9AABk\nTqgHAIDMCfUAAJA5oR4AADJX1dcdANgbNTU1JdVVVlYm1RUKhaS65ubmpLpSVFSkzTe1tLQkt1lK\nbYpSXtft+9r67/r6+l3WVVWlfeRv27Ytqa6U17RYLCbV9evXL6kudcyVInUMlPK6VldXJ7fZnbHW\nauDAgUntlSL1dU39/tidmKkHAIDMCfUAAJA5oR4AADIn1AMAQOaEegAAyJxQDwAAmRPqAQAgc0I9\nAABkLv+V9gEylMtGJ6mbXfWFUvra288zdROgN2rdNGnQoEFluV57Ujd0Ys/SG2OtHHJ5b+0JZuoB\nACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcA\ngMwJ9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZE+oBACBz3Q71S5cujSlTpuxwbNu2bXHbbbfF\naaedFsccc0ycddZZ8aMf/ahsnQQAADpW1Z0Hr1ixIq688sqdjn/pS1+KRx99ND73uc/FuHHj4tFH\nH43LL788Kioq4n3ve1/ZOgsAAOysS6G+oaEhFi9eHLfffnsMGjQoGhsb286tW7cuHnnkkbjhhhvi\n7/7u7yIi4p3vfGf87//+b9xzzz1CPdAtxWIxubZQKJSxJ7vW3Nzc4bnW59HRYyoq0u5+bGlp6dW6\n1H5GxA6fFbu7jRs3JtU1NTUl1Q0ZMiSpLmLHMdX67w0bNuyybtOmTUntpY6BmpqapLqIiDVr1iTV\njRw5Mqku9fujFNXV1Ul1lZWVyW0OHDgwqa6xsbHtNWpoaOhyXepzLEVn78udqarq1jz3bqlL36mP\nPfZYLFq0KGbPnh0zZszY4dzmzZvjvPPOi+OPP36H4+PGjYvVq1eXr6cAAEC7uhTqJ02aFEuXLo3p\n06fvNBM2evTo+NKXvhT7779/27GWlpZ47LHH4pBDDilvbwEAgJ10KdSPHDmyW79G+9rXvhbPPfdc\nXHTRRckdAwAAuqbsNxAtWLAg5s+fHxdeeGG8+93vTrpGbW1tmXsFf7Nly5aIMMZ2ZzndU99ZX7du\n3RoREU8//XRvdWe3U8rXsrel3oeb6tVXXy3LdVrvb66rq9vlY1OfY+r3VSnPMfXvMbryOuwuUl/X\nUt7nUmuLxWJs27YtIiL+/Oc/93h7fSGnvnakrKH+xhtvjMWLF8eMGTPiiiuuKOelAQCADpQl1BeL\nxbjyyivjBz/4QXzyk5+Mf/zHfyzpehMmTChHt2AnrTP0xtjuK6eZ+s5mPltn6I844oh2z1v9ZveS\n6+o3rTPTo0eP3mWd1W86ZvWbzjU2NrbN0B9++OFdrrP6TfktX768w3NleQY33nhj/OAHP4jZs2fH\nBRdcUI5LAgAAXVRyqP/jH/8Y3/zmN2PatGnx1re+NZ588sm2cxUVFTFx4sRSmwAAADpRcqj/2c9+\nFhERjz/+eDz++OM7nBs4cGCsWLGi1CYAAIBOFIq72fIEy5cvj6lTp/Z1N9hDuaee3mKs0RuMM3qL\nsbZ76Cwnp/8FFAAAsFsQ6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8A\nAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAA\nZK6qrzsA7JmKxWKvt9nS0pJU19jYmFTXv3//Ds+1Pv+OXofUNrdt25ZUV1GRNoeT2l4pbaa+NtXV\n1Ul1ERGbN29Oquvtr0dH7a9ataos12vP1q1bk+pGjhyZ3Gbq93JTU1NS3eDBg5PqIiIqKyuT6so1\nBroj9XukWCy2vZd15zVOfW1KUSgUer3N3YWZegAAyJxQDwAAmRPqAQAgc0I9AABkTqgHAIDMCfUA\nAJA5oR4AADIn1AMAQOZsPgX0iL7YACR1o5Oe2CCl9fl39Dr069cv6bqpdalK2ZQnJ8OGDevrLiSp\nra2NiIixY8f2bUfY47W+l1VViY67KzP1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNC\nPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1\nAACQOaEeAAAyJ9QDAEDmqrpbsHTp0rjiiitixYoVOxy/++6744EHHojXXnstpkyZEtdee20ccsgh\nZeso7I6am5uT6iorK3u1vVIUi8WkuqampqS66urqpLqIiG3btiXVtbS0JNVt3bq1w3Otz/+vf/1r\nu+dff/31pDaHDh2aVNdRP3aloiJ97mfz5s1JdaljYP369Ul1EREHHnhgUl1nY6Az69atS6qLiBg/\nfnzbv3c1zraX+vqkjoHhw4cn1UVEbNmyJalu4MCBSXWlvO8UCoWkuv79+ye3mSr1/byqqttxsc+k\nPsfUr+PupFvfqStWrIgrr7xyp+Pz5s2L+fPnx6xZs+K2226LjRs3xsyZM2PTpk1l6ygAANC+LoX6\nhoaGWLhwYVxwwQU7/bS2efPmuOeee+Izn/lMTJ8+PU4++eT4xje+EZs2bYoHH3ywRzoNAAD8ny6F\n+sceeywWLVoUs2fPjhkzZuxw7sknn4wtW7bEySef3HZs6NChceyxx8YvfvGL8vYWAADYSZdC/aRJ\nk2Lp0qUxffr0ne45eu655yIi4uCDD97h+OjRo2PVqlXl6SUAANChLv3lw8iRIzs8t3nz5ujXr99O\nt+UMHjzYPfUAANALSl7SslgsdvgXw6WsmgAAAHRNyWsU1dTURENDQzQ3N++wTN/mzZtjyJAhSdes\nra0ttVvQrtZl0so1xnp76azU9vpCXywrltpmT9Q1NjZGRMSaNWvaPZ+6PGnqMpGpS4yWInWp0IaG\nhqS6UpYITF1iNPU5ltLX7cfUrsbZ9lLHQOr3ZOpYjej9949S3nf6os3eVigUyv752VMsaVmCsWPH\nRrFYjNWrV+9wvK6uLsaNG1fq5QEAgF0oeab+mGOOiX79+sWSJUviwgsvjIiIDRs2xP/8z//EZz7z\nmaRrTpgwodRuQbtaZxjKNcZsPtWxvX3zqdaZ0wMOOKDd8zaf6lhfbD613377JdX1xeZT20+Y7Wqc\nbc/mUx2z+VTnqqqqyv752VP29Jn65cuXd3iu5FA/aNCgmDFjRnzta1+LQqEQY8aMia9//esxdOjQ\nOPvss0u9POzWUsN5Lu2Voi92IBw0aFCvtldTU9PhuVdeeSUiIkaMGNHu+Y6O95Tebo/esatxtr2c\nxsCwYcP6ugtkKpdw3hOSPnXf+IJdfvnlUVlZGffcc0/U19fHlClT4itf+UqnH3gAAEB5dDvUX3rp\npXHppZfucKyysjIuv/zyuPzyy8vWMQAAoGusOQkAAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEe\nAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoA\nAMicUA8AAJmr6usOQM6ampqS6orFYlJdoVBIqmtpaUmqi4ior69Pqquurk6qS31tIiI2bdqUVPf8\n888n1e2///4dntu2bVtERKxatard8y+++GJSm1VVaW/bGzduTKpLfU0jIvbbb7/k2hSvvvpqcm1l\nZWVS3WGHHZZUt2XLlqS6iIiBAwe2/Xvr1q0REfHUU0/tsi71/WrUqFFJda+//npSXUT62Bk0aFBS\nXSnvO6nvdaW8L6dKHefFYrHtNerOa5X6mUUaM/UAAJA5oR4AADIn1AMAQOaEegAAyJxQDwAAmRPq\nAQAgc0I9AABkTqgHAIDM2XwKSpC6EVBO+vXr19dd6LKampqkutTNdTpTW1sbERFjx45t93xHx6E7\nWsfZ+PHj+7gnOxsxYkRfd2G3lroRVF8oFAptG0nZUGr3ZaYeAAAyJ9QDAEDmhHoAAMicUA8AAJkT\n6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6o\nBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMhcVV93gD1TsVhMri0UCj3WZutjtn9sY2NjUnsRERUV\naT8XNzQ0JNU1Nzcn1Q0YMCCpLiJi/fr1SXUbNmxIqnv11VeT6iLSx87zzz+fVPf44493eO6MM86I\niIgFCxa0e/5tb3tbUpsDBw5MqluzZk1S3erVq5PqIiJOP/305NoUf/nLX5Jrjz766KS6N7/5zUl1\n27ZtS6qLiKisrNzpOl0Zw8OGDSu5ve6oqkqPGE1NTUl1gwcPTqpLfW8tRerr01efr+19fvZUe6Qx\nUw8AAJkT6gEAIHNlC/UtLS2xcOHCOO200+KYY46Jc845J371q1+V6/IAAEAHyhbqFy1aFP/yL/8S\nZ599dtx1110xevTomDVrVjz11FPlagIAAGhH2UL9I488Eh/84Afj4osvjne+851x8803x5ve9KZ4\n8MEHy9UEAADQjrKF+oaGhh3+6ryioiJqamqSV84AAAC6pmyhfvr06fG9730vli1bFps2bYrFixfH\nypUr48wzzyxXEwAAQDvKtk79eeedF7/61a9i5syZEfG3tUk/97nPxUknnVSuJgAAgHYUiqXsYrCd\n888/P/7yl7/EZZddFoccckg8/vjjsXDhwrjmmmvi/PPP7/J1li9fHoMGDSpHl+hDu+vmU1u3bo2I\nHTdj2l37Ws661E2yItI3gUndzKUvNoFJ3Qxs06ZNHZ7bZ599IqLjTbhSN8lJHXOpm6yVsjnb0KFD\nk2tTlLKhU+qmXtXV1Ul1Zfrobft+6coGUaVsBpWilI2H+uK9rrf19udHqW229/nZU+3Rsfr6+pg6\ndWq758ryHb58+fJYsWJF3H777XHaaadFRMSxxx4bTU1NcfPNN8dZZ52V/IZJnvriG7k7bW7/2N29\nr30tNbSk1vWF1IC97777dnhuy5YtERExevTopGvTfUOGDOnrLvS61nHWv3//Pu4JvaWvP7Ny+vza\n25Ql1K9duzYKhUJMnjx5h+NTp06NRYsWxQsvvBCHHXZYl683YcKEcnQLdlJbWxsRxhg9z1ijNxhn\n9BZjbfewfPnyDs+V5fdUY8eOjWKxuFNDv/vd76KysjJGjRpVjmYAAIB2lGWm/qijjoqTTjop5s6d\nG+vXr49DDz00fv3rX8eiRYviggsuiJqamnI0AwAAtKNsfzVz++23x2233Rbz58+PDRs2xJgxY+KL\nX/xinHPOOeVqAgAAaEfZQn2/fv3iqquuiquuuqpclwQAALogn7WfAACAdgn1AACQOaEeAAAyJ9QD\nAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8A\nAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAA\nZK6qrzvAnqmpqanX29y2bdsuH9Pc3BwREZs3b247tn79+uQ2+/fvn1TXlb62Z8OGDUl1zz//fFJd\nKV5++eWkutavUYq//vWvSXVvfvObk9vsyNixYyMi4tvf/na755944omk66b2taGhIamuX79+SXUR\nESNGjEiqGzlyZFLdxo0bk+oiIsaPH59Ul/ocq6urk+oidvyatL7Xrlu3bpd1qe9XlZWVSXVVVekR\nI7U29f2joqL35zgLhUKvt1ksFpPq+qKvdJ+ZegAAyJxQDwAAmRPqAQAgc0I9AABkTqgHAIDMCfUA\nAJA5oR4AADIn1AMAQOZsPkWPKGXTkZ5ss3UTlcGDB7cd2/7fu7vUjYeOPPLIMveEXamtrY2IiOOP\nP77d89OnT+/N7rCHeumllyIiYvjw4X3ck91D6kZZewubSO3ZzNQDAEDmhHoAAMicUA8AAJkT6gEA\nIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCA\nzJU11C9btizOOeecmDx5crznPe+JO+64I4rFYjmbAAAA3qBsoX758uVx0UUXxWGHHRYLFiyIGTNm\nxMKFC+Ouu+4qVxMAAEA7qsp1oVtvvTVOPPHEuOGGGyIi4rjjjov169fHr3/96/j0pz9drmYAAIA3\nKEuoX7duXaxYsSLuvvvuHY5ffvnl5bg8AADQibLcfvPMM89ERET//v3jE5/4REyaNCmmTZsW8+bN\nc089AAD0sLLN1BeLxZg9e3aceeaZ8Q//8A/xxBNPxN133x0DBgyIWbNmlaOZvV4pPyAVCoWkuubm\n5qS6xsbGpLqIiPr6+qS6zZs37/IxDQ0NERFRV1fXduzFF19Mai8i/fVZu3ZtUt3LL7+cVLds2bKk\nuoiI6urqpLphw4Yl1VVUpM81rFu3LqnupZdeSqqbMGFCh+dOPfXUiIi499572z0/dOjQpDZTHXjg\ngUl1BxxwQHKbo0aN6tW6rrwHdGTEiBFJdS0tLUl1AwYMSKqL2PF7pLX9rVu37rIu9Xs5VSmfWVVV\nZbs7GPYaZfmuaWpqioiIE088Ma644oqIiHj7298er732Wtx9991x4YUXditU1tbWlqNbe5y+CPWp\nbZbS19QPya7UtY7V7cNfZWVlUnsR6QE0NbSkBo8xY8Yk1UWkj51SXtdUvf1DaGfBrDW0t4b7N+rt\n1yc10JUSBFO/l1955ZWkutSvf0TEli1bkmtTlPLD6/a2bdsWERErV67c5WNTv5f7Qk593Vu0fo/I\naLuvsryrDBo0KCIiTjjhhB2OT5s2Lerr62P16tXlaAYAAGhHWWbqW2cB3zjb1Tor2t2fuDv7lfbe\nzO03nevKr95bZ+iHDx/edsztN51z+03HunL7zU9/+tN2z7v9pmNvetObkur2xttvWmfoDz300F3W\nuf2GUrTO0MtofWv58uUdnivLTP1hhx0W+++/f/zkJz/Z4fjPf/7zGDlyZBx00EHlaAYAAGhHWUJ9\noVCIyy67LB599NGYM2dOLFu2LL761a/G9773vbj00kvL0QQAANCBsv1+68Mf/nD069cvvv71r8fD\nDz8co0aNirlz58ZHPvKRcjUBAAC0o6w3rb3//e+P97///eW8JAAAsAvlWVMLAADoM0I9AABkzppR\nGemLzThSN8gpZWOd1KXetl+msiObNm2KiIjRo0e3Hdv+33uqT3ziE33dhb1O6/JvX/nKV/q4J+zJ\nWpe3LGWJTGDPYKYeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEe\nAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoA\nAMhcVV93IGfFYjGprqmpKamuoaEhqS4iva+p1q5dm1z73HPPJdVt3Lhxl48ZNWpURET8+7//e9ux\nNWvWJLUXEfHb3/42qW7w4MFJdY8//nhSXevzTrFy5cqkuv79+yfVHXTQQUl1ERGHHHJIUt0BBxyQ\nVDdmzJgOzx188MEREfGd73yn3fOp3yPTpk1LqisUCkl1o0ePTqqLiNiyZUtS3bBhw5Lq9t1336S6\niIht27Yl16ZI/f6IiGhpaWn7d+t7e3Nzc8l96khlZWVSXW9/7kSkj3PYE5ipBwCAzAn1AACQOaEe\nAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZM7mUyVI3eSiurq6V+v6wmGHHdYntbtSW1sb\nEekb+EBXtY61448/vo97QlcMGDCgr7vQZdtvBtX6OZS6QVRPshEU9C4z9QAAkDmhHgAAMifUAwBA\n5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZ\nE+oBACBzQj0AAGROqAcAgMyVPdQ3NDTEGWecEVdffXW5Lw0AALSj7KF+3rx58dxzz5X7sgAAQAfK\nGur/9Kc/xTe/+c0YPnx4OS8LAAB0oqpcF2pubo4vfOELMWvWrPjpT39arsv2iqampqS6QqGQVLdl\ny5ZerYuIeOWVV5LqnnrqqaS6V199NakuIuLPf/5zj9VdeOGFERFxzTXXtB0bNmxYUnsREUuXLk2q\nO/7445PqNmzYkFSXOsYjIvbdd9+kuvHjxyfVlTLOU9tMHQNHHnlkh+e2bt0aER1/Dw0ePDipzerq\n6qS6gQMH9mp7EREVFWnzRpWVlcltpkp9P6+qSvsYLRaLSXUR6X0F9mxlm6lfsGBBNDU1xSWXXFKu\nSwIAAF1eqXfxAAASxUlEQVRQlpn6lStXxvz58+P+++9PnrUAAADSlDxTXywW49prr42PfOQjMWnS\npHL0CQAA6IaSp9Xvv//+WLt2bSxcuDCam5t3uE+wubk56d7I2traUrvVLaXc25iipaWlV+si0u+p\nHjVqVFLdm970pqS6iIixY8cm1b3rXe/a5WNGjBgREf93b31EaffvnnHGGUl1qfdTNzc3J9WVIvX1\nGTBgQFJdKeM8tc3U59jZ/f+tz6Ojx2zbti2pzdT7qVPvb++L+7dzumc8ta/luqe+dXz19ucmex9j\nbfdXcqhfsmRJrF27Nt72tre1HSsUCvHUU0/FI488EkuXLo0DDzyw1GYAAIAOlBzq//mf/zk2b968\nw7HPf/7zMW7cuPjMZz4TI0eO7PY1J0yYUGq3usXqNx1btWpVUt3uvvrNN77xjbZjVr/pnNVvOtbZ\n6jetz6OjVWesftMxq990bvu+ts6a9vbnJnsfY233sHz58g7PlRzq27tVYsCAATFs2LBOP/AAAIDy\nKPuOshF/m0XI6Z5IAADIWY+sP/nwww/3xGUBAIB2WFQ+0u+JTFVTU9OrdRHpq9HsabdQtd4T6AdP\nelrrWEu9z589m99mA+XWI7ffAAAAvUeoBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT\n6gEAIHNCPQAAZE6oBwCAzAn1AACQOaEeAAAyJ9QDAEDmhHoAAMicUA8AAJkT6gEAIHNCPQAAZE6o\nBwCAzAn1AACQuaq+7kC5NDc3J9c2NjYm1b366qtJdRs3bkyqe/3115PqIiJ+/OMfJ9Wl9vX73/9+\nUl1ERF1dXVLdwIEDd/mYO++8MyIijj/++LZjZ599dlJ7ERGvvPJKUt3TTz+dVDd58uSkuiFDhiTV\nRUQMGzYsqS71dU0dcxERo0ePTqrr379/Ul1nr03re9KGDRvaPV9fX5/U5qBBg5LqUp9jRUX63E91\ndXVybYpCodCr7QHsTszUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAA\nMrfHbD5VWVnZ67VvfvObk9vsbW9/+9t7tb1bbrmlV9vrqtra2oiIWLduXR/3hD3diy++GBER++yz\nT7vnOzoOACnM1AMAQOaEegAAyJxQDwAAmRPqAQAgc0I9AABkTqgHAIDMCfUAAJA5oR4AADIn1AMA\nQOaEegAAyJxQDwAAmRPqAQAgc0I9AABkTqgHAIDMCfUAAJA5oR4AADIn1AMAQOaEegAAyJxQDwAA\nmasq14VaWlpi8eLF8d3vfjfWrFkTBx54YJx//vkxffr0bl+rqakpqf1Ua9asSarbuHFjUt23v/3t\npLolS5Yk1ZXiN7/5TVLdmWeemdzmn/70p6S6LVu27PIxt912W0REnHLKKW3HDjnkkKT2IiLOO++8\npLrU8TpgwICkusMPPzypLiLioIMOSqrbf//9k+oqKtLnGorFYlJdVVXZ3grbtH6Nt27d2u756urq\npOumPsfKysqkulIUCoWkutTnmNoewJ6gbJ9kd955ZyxatCg+/elPx6RJk+I3v/lN3HDDDbF169a4\n8MILy9UMAADwBmUJ9S0tLXHffffFrFmz4uKLL46IiHe84x2xbt26uOeee4R6AADoQWW5p37Tpk1x\n1llnxamnnrrD8XHjxsW6des6/PUzAABQurLM1A8dOjSuvfbanY4/+uijMWrUqOT7gAEAgF3rsdVv\nvvvd78ayZcti1qxZPdUEAAAQEYVi6jIDnfj+978fV199dZx22mltq4101fLly2PgwIHdbrOUp9HY\n2JhUl7qCybp165LqXn/99aS6UtTX1yfV7bPPPsltpt6u1ZWvR+tKLqtXr2471r9//6T2IiL23Xff\n5NoUqSvDlPIcU1dpSa3rCz2xasq2bdsiouPXfm9YqcXqNz2vddWvlM9N6A5jbfdQX18fU6dObfdc\n2Wfq77333rjqqqviPe95T9x8883lvjwAAPAGZV2c+dZbb40FCxbEWWedFddff33yrOJb3vKWbtfk\ntE790qVLk+qsU9+57qxTf9lll7Uds05956xTn2blypUREXHooYe2e9469R0zU991tbW1ERExYcKE\nPu4JezpjbfewfPnyDs+V7ZNs8eLFsWDBgvj4xz8es2fPLtdlAQCAXShLqH/llVfiq1/9ahxxxBFx\nxhlnxJNPPrnD+YkTJ3Zr9q0nZs06M2bMmF5t78Ybb+zVOv5P60zDiy++2Mc9YU/X+p5n9a/u2xtn\n3AFKVZb0/Mtf/jIaGxvjmWeeiXPPPXen88uWLYthw4aVoykAAOANyhLqzzrrrDjrrLPKcSkAAKCb\nemydegAAoHcI9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAA\nkDmhHgAAMifUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMxV9XUH2rNx48Zu19TV1SW3\nd8899yTVPfzww0l1zc3NSXVvectbkuoiIv7whz8k1dXU1CTVzZo1K6kuIuLll19Oqps4ceIuH3Po\noYdGRMTixYvbjg0aNCipvYiIo446Kqlu27ZtSXXDhw9Pqttvv/2S6iIi+vXrl1TX1NTUq+1FRFRW\nVibVFQqF5DYpv2KxmFTn6wjszczUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAA\nkDmhHgAAMrdbbj41ZMiQbtcceeSRye3dcsstvVpH36mtrY2IiBNOOKGPe7LnK2UTKfZuNpEC6D4z\n9QAAkDmhHgAAMifUAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifU\nAwBA5oR6AADInFAPAACZE+oBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5qr6ugPt+elP\nf9rtmhdffDG5vZdeeimprra2Nqnu4IMPTqo75ZRTkuoiIn72s58l1Z1xxhlJdSNGjEiqi4iorq5O\nqhs4cOAuH9PU1BQREX/961/bjg0YMCCpvYiIfv36JdVVVKT9PF0oFJLqKisrk+pK0dzcnFTXF33d\nGxSLxaS61DEHQO8yUw8AAJkra6h/4IEH4vTTT4/JkyfHueeeG7/73e/KeXkAAKAdZQv1Dz/8cMyZ\nMyc+9KEPxR133BFDhw6NWbNmxQsvvFCuJgAAgHaULdTfcccdce6558anPvWpeNe73hV33XVXDBs2\nLO67775yNQEAALSjLKH++eefjxdffDFOPvnktmNVVVVx0kknxS9+8YtyNAEAAHSgLKF+1apVUSgU\nYsyYMTscP+igg6Kuri551QUAAGDXyhLqN23aFBERgwcP3uH44MGDo6WlJerr68vRDAAA0I6yhPrW\nmfiO1jNOXZMbAADYtbJsPjVkyJCIiNi8eXMMHz687fjmzZujsrKyS5sCbW/o0KHd7kN329he6mZQ\nEydOTKpL3bColE2STjjhhKS61B/IXn/99aS6iPTNbrrS18bGxoiIWLNmTcntlVLb2xv69MUGQnv7\nZkdbtmyJiPRN6sptb/967Kl2t3HGnstY2/2VJdSPGTMmisVi1NXVxejRo9uOr169OsaOHdv9TlV1\nv1spNa1SfyBI+eGjrwwbNqxX20vdTbS3NDQ09HUX2Eu4/ZDeYJzRW4y13VdZQv3YsWPjgAMOiCVL\nlsS0adMi4m8zoj//+c93WBGnK6ZOnVqOLgEAwF6jLKE+IuKiiy6K6667LoYMGRJTpkyJb33rW7F+\n/fq44IILytUEAADQjkKxjOtN3nfffXH//ffHa6+9FuPHj4+rr746Jk2aVK7LAwAA7ShrqAcAAHqf\ntSYBACBzQj0AAGROqAcAgMwJ9QAAkDmhHgAAMifUAwBA5oR69jgPPPBAnH766TF58uQ499xz43e/\n+12nj1+xYkX8/d//fRx77LFx4oknxlVXXRWvvvpqL/WWXHV3nG1v3rx5MX78+B7sHXuS7o61devW\nxZVXXhnHHXdcHHvssfHJT34y6urqeqm35Crls/P888+PKVOmxHvf+96YN29eNDU19VJvaY9Qzx7l\n4Ycfjjlz5sSHPvShuOOOO2Lo0KExa9aseOGFF9p9/MqVK2PmzJkxZMiQuPXWW2P27NmxYsWKmDVr\nVjQ3N/dy78lFd8fZ9p555pmYP39+FAqFXugpuevuWGtqaoqZM2fGH/7wh7j++uvjpptuirq6urjo\noosELjrU3XFWV1cXF154YdTU1MS8efNi5syZsWjRorj11lt7uefsoAh7kJNPPrk4d+7ctv9vbGws\nnnLKKcXrrruu3cfPnTu3+N73vrfY1NTUduz3v/998Ygjjij+13/9V4/3lzx1d5y1am5uLp599tnF\nd7/73cXx48f3dDfZA3R3rD3wwAPFt771rcW1a9e2HautrS2eeOKJxT/+8Y893l/y1N1xNn/+/OLk\nyZOLW7dubTt26623FqdOndrjfaVjZurZYzz//PPx4osvxsknn9x2rKqqKk466aT4xS9+0W7N4Ycf\nHjNnzozKysq2Y+PGjYuIiNWrV/dsh8lSyjhrde+990Z9fX3MmDGjp7vJHiBlrC1dujROPPHE2H//\n/duOjR8/Ph577LE48sgje7zP5CdlnDU2NkZVVVX079+/7dg+++wT9fX10dDQ0ON9pn1CPXuMVatW\nRaFQiDFjxuxw/KCDDoq6urooFos71Zx33nlx/vnn73Ds0UcfjUKhEIccckiP9pc8pYyziL99cM6b\nNy+uu+66qK6u7o2ukrmUsfb000/HuHHjYt68eXHCCSfExIkT45JLLok1a9b0VrfJTMo4++AHPxiV\nlZVxyy23xIYNG+L3v/993H///XHqqadGv379eqvrvIFQzx5j06ZNERExePDgHY4PHjw4Wlpaor6+\nfpfXWLNmTXzlK1+JiRMnxjve8Y4e6Sd5Sx1n1157bZx11llxzDHH9Hgf2TOkjLV169bFQw89FL/8\n5S/jhhtuiJtvvjmeffbZuOSSS6KlpaVX+k1eUsbZ6NGj44orroh77rknjjvuuDjnnHNiv/32ixtu\nuKFX+kz7qvq6A1AurbMJHf0BYkVF5z/DrlmzJj7+8Y9HRPhjHzqUMs6+853vRF1dXcyfP79H+8ae\nJWWsNTU1RVNTUyxatChqamoi4m8zrmeffXb853/+Z7zvfe/ruQ6TpZRx9t3vfjf+3//7f3HuuefG\nGWecES+//HLcfvvtcfHFF8d9993nt5F9xEw9e4whQ4ZERMTmzZt3OL558+aorKyMgQMHdlj7zDPP\nxLnnnhv19fVx7733xkEHHdSjfSVf3R1na9eujVtuuSW+8IUvRP/+/aO5ubltxrS5ubnD23Ug5T1t\n0KBBMXny5LZAHxFx9NFHx9ChQ+OZZ57p2Q6TpZRxtnDhwjjppJNizpw5cdxxx8UHPvCBmD9/fixf\nvjz+4z/+o1f6zc6EevYYY8aMiWKxuNN6zKtXr46xY8d2WPfkk0/GjBkzorq6Ov71X/81Dj/88B7u\nKTnr7jhbtmxZ1NfXx2c/+9k46qij4qijjoovf/nLUSwW4+ijj44777yzl3pOblLe0w4++OBobGzc\n6XhTU5NlVGlXyjhbs2ZNTJ48eYdjhxxySAwbNiyeffbZnuoquyDUs8cYO3ZsHHDAAbFkyZK2Y42N\njfHzn/883vnOd7Zbs3r16rj44otj5MiR8W//9m8xevTo3uoumeruOHvPe94TDz74YDz44IPx0EMP\nxUMPPRQzZ86MQqEQDz30UHz0ox/tze6TkZT3tBNOOCFWrFgRr7zyStuxJ554Iurr62PKlCk93mfy\nkzLOxo4dG7/97W93OPb888/H+vXrfY72oco5c+bM6etOQLn069cv7rrrrmhoaIiGhoa48cYbY9Wq\nVXHTTTfF0KFDo66uLlatWhWjRo2KiIirrroqnn322bjmmmsiIuKll15q+6+ysnKnPxyCiO6NswED\nBsTIkSN3+G/lypXxy1/+MubOnWuM0anuvqcdccQR8dBDD8XSpUtjxIgR8cc//jHmzJkT48ePj899\n7nN9/GzYXXV3nA0fPjwWLFgQa9eujYEDB8Zvf/vb+OIXvxhDhw6NOXPmuKe+jxSKbuhkD3PffffF\n/fffH6+99lqMHz8+rr766pg0aVJERFx99dXxyCOPRG1tbTQ1NcVb3/rWDneOvfLKK2PmzJm92XUy\n0tVx1p7FixfHTTfd1OF52F53x1pdXV18+ctfjmXLlkVVVVWccsopcc011+xwnz28UXfH2ZIlS+Ku\nu+6KZ599NkaMGBHHH398XHbZZTF8+PC+egp7PaEeAAAy5556AADInFAPAACZE+oBACBzQj0AAGRO\nqAcAgMwJ9QAAkDmhHgAAMifUAwBA5oR6AADI3P8H15rkR4/IBEQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11404bbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out=gibbs()\n",
    "plt.hist2d(out[:,0],out[:,1], normed=True, bins=50)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}