{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Gibbs Sampling" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(6, 8)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import random\n", "def roll_a_die():\n", " return random.choice([1,2,3,4,5,6])\n", "\n", "def direct_sample():\n", " d1 = roll_a_die()\n", " d2 = roll_a_die()\n", " return d1, d1 + d2\n", "\n", "direct_sample()" ] }, { "cell_type": "code", "execution_count": 156, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(3, 6)" ] }, "execution_count": 156, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def random_y_given_x(x):\n", " \"\"\"equally likely to be x + 1, x + 2, ... , x + 6\"\"\"\n", " return x + roll_a_die()\n", "\n", "def random_x_given_y(y):\n", " if y <= 7:\n", " # if the total is 7 or less, the first die is equally likely to be\n", " # 1, 2, ..., (total - 1)\n", " return random.randrange(1, y)\n", " else:\n", " # if the total is 7 or more, the first die is equally likely to be\n", " # (total - 6), (total - 5), ..., 6\n", " return random.randrange(y - 6, 7)\n", "\n", "def gibbs_sample(num_iters=100):\n", " x, y = 1, 2 # doesn't really matter\n", " for _ in range(num_iters):\n", " x = random_x_given_y(y)\n", " y = random_y_given_x(x)\n", " return x, y\n", "\n", "gibbs_sample()" ] }, { "cell_type": "code", "execution_count": 154, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 154, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFndJREFUeJzt3X+MXWWdx/H3pz+hgK1s3Wb6I2mzVk01S2mbilujSEVr\nJRazLmld2aKYslnQdtcEwf2jEJcNG3+gRtM4UpaitdAghIYQpZYqq1mK01qxtCBdqNvpDowVRFBp\nbfvdP+4ZuAz3zr1zz5l7zznzeSUnvee5Z57znEzy7TPf8/xQRGBmZp01ptMNMDMzB2Mzs1xwMDYz\nywEHYzOzHHAwNjPLAQdjM7McaBiMJd0iqV/SvqqyL0h6TNIjku6WNGVkm2lmVm7N9IxvBZYNKtsO\nvC0i/hr4FXBtxu0yMxtVGgbjiHgQeHZQ2f0RcSI5fQiYOQJtMzMbNcZlUMcngDvqfSlpDbCmcjZ+\nIUzN4JZmViYL5/S9pmz3UxyNiDekqfeNUvyxiev64AcRMTgD0FapgrGkfwVOAJvrXRMR3UB35frp\n8XJcNjNL9Nxw/WvK9FF+nbbePwJXNHHddTnoJbYcjCVdBlwELA0vcGFmOSSy+fO/HVpqp6RlwNXA\nuyOa+ivAzKztxgCnd7oRTWoYjCVtAc4HpkrqBdZTGT0xEdguCeChiPjHEWynmZWYPrq+RulrUxfD\nrhcYn7qW9mgYjCNiVY3ijSPQFjOzTJU+TWFmVgSl6hmbmRWVe8ZmZjngnrGZWQ6UajSFmVlRuWds\nZpYTRQlyRWmnmdmwuWdsZpYDHk1hZpYDRXqB522XzKy0BtIUjY6G9UinSXpY0i8kPSrp+qR8jqRd\nkg5KukPShKR8YnJ+MPl+dqN7OBibWWkNpCkaHU04BlwQEecA84Flks4D/gO4KSLeCDwHXJ5cfznw\nXFJ+U3LdkJymMCu445NrL6gz4flai+/kU8yvsZ7x3vT1ZvUCL1km+MXkdKBDHcAFwEeT8k3AdcAG\nYEXyGeBO4OuSNNRyw+4Zm1lpDaNnPFVST9Xxml0wJI2VtBfop7IP6P8Av6vagq4XmJF8ngEcBki+\nfx74i6Ha6p6xmZXWMHrGRyNi0VAXRMRJYL6kKcDdwFvStq+ag7GZlZbIfjRFRPxO0k7gHcAUSeOS\n3u9M4Ehy2RFgFtAraRwwGfjtUPU6TWFmpSVg/LjGR8N6pDckPWIknQ5cCBwAdgIfSS5bDdyTfN6W\nnJN8/0Cj7encMzYruCK9qKtHe0dopw/BuGai3ImGV3QBmySNpdKJ3RoR90raD9wu6d+An/PKxhsb\ngW9LOgg8C6xsdAMHYzMrLQnGj01fT0Q8Apxbo/xJYHGN8peAvxvOPRyMzay0mu4Z50BBmmlmNnwS\njJ/Y6VY0x8HYzMqrQCsFFaSZZmYtcDA2M8uJgkS5gjTTzKwFAjIYTdEODsZmVl5OU5iZ5YAAj6Yw\nM+sw94zNzHLAwdjMLCf8As/MrMPcMzYzywEHYzOzHPBoCjOzHHDP2Cyf4ooauxB/s9iLs8cnay/C\nrpuL81zraywkn35peRyMzcxyoUDToRvugSfpFkn9kvZVlZ0tabukJ5J/Xz+yzTQza8FAz7jRkQPN\nbEh6K7BsUNk1wI6ImAvsSM7NzPJl4AVeoyMHGgbjiHiQyoZ61VYAm5LPm4CLM26XmVl6BeoZt9qM\naRHRl3x+GphW70JJa4A1lbPJLd7OLBtFf1lXS5Fe1NVzPSOzO/SoeoEXESEphvi+G+gGkKbXvc7M\nbEQUJBg3kzOu5RlJXQDJv/3ZNcnMLCMDoykaHTnQajDeBqxOPq8G7smmOWZmGSpQzriZoW1bgP8G\n3iypV9LlwI3AhZKeAN6bnJuZ5UtGoykkzZK0U9J+SY9KWpuUXyfpiKS9ybG86meulXRQ0uOS3t/o\nHg3/T4iIVXW+Wtr4EczMOii7F3gngM9ExB5JZwG7JW1PvrspIr74qttK84CVwFuB6cAPJb0pIk7W\nu0GraQozs/zLKE0REX0RsSf5/AJwAJgxxI+sAG6PiGMR8RRwEFg81D0cjM2svJoPxlMl9VQda+pW\nKc0GzgV2JUVXSXokma08MBt5BnC46sd6GTp4OxibWck1N5riaEQsqjq6a1Ul6Uzge8C6iPg9sAH4\nK2A+0Ad8qdVm5uQ9opnZCMhw0oek8VQC8eaIuAsgIp6p+v5bwL3J6RFgVtWPz0zK6nLP2MzKK7vR\nFAI2Agci4stV5V1Vl30YGFhQbRuwUtJESXOAucDDQ93DPWMzK6/sesZLgEuBX0ram5R9DlglaT4Q\nwCHgCoCIeFTSVmA/lZEYVw41koLMmmlmlkcZBeOI+ElS22D3DfEzNwA3NHsPB2MzK68CLS7vYGxm\n5TWaVm0zM8stAad1uhHNcTA2s/JymsKK7vmJtRf2nnys2AuZ76yxYPl7ai5sXhzfrrMI+6UFeq54\nd41du3+cQcVOU5iZ5URBolxBmmlm1gKnKczMcqBAaQpFtG9busoeeHUXQzIzq3L97ohYlKaGRW9R\n9GxsfJ3eSep7pVWQ/zPMzFpQoJ5xQZppZtYCB2MzsxxwMDYzywmPpjAz6zD3jM3McmBgcfkCcDA2\ns/Jyz9jMLAccjM3McsDB2MwsH8KjKczMOivGwHEvLm9m1lkhODF2TBNXnhrxtjTiYGxmpRUSJ8c1\nE+aOj3hbGnEwNrNSOzm2GEljB2MzK61AnCzIfOhmkilmZoUUiBOMbXg0ImmWpJ2S9kt6VNLapPxs\nSdslPZH8+/qkXJK+JumgpEckLWh0DwdjMyutQBxnYsOjCSeAz0TEPOA84EpJ84BrgB0RMRfYkZwD\nfACYmxxrgA2NbuA0RQbiiho7236zODvz1hIX1N5xWA8U/LlW1fhdbSn2Mx0dW/t3NfVkcZ4r1tb4\nvXw1g3ozSlNERB/Ql3x+QdIBYAawAjg/uWwT8CPgs0n5bVHZSukhSVMkdSX11JSqZyzpn5Mu+z5J\nWyQVZESfmY0WJxnb8ACmSuqpOuruDydpNnAusAuYVhVgnwamJZ9nAIerfqw3Kaur5Z6xpBnAp4F5\nEfEnSVuBlcCtrdZpZpalgZxxE442sweepDOB7wHrIuL3kl65V0RIanlT0bRpinHA6ZL+DEwC/i9l\nfWZmmamkKbLJxkoaTyUQb46Iu5LiZwbSD5K6gP6k/Agwq+rHZyZl9etPszt08kbxBuBPwP0R8fc1\nrlnDy1tCT14I61q+n5mNJul3h5636PTY3DOn4XULdGDIe6nSBd4EPBsR66rKvwD8NiJulHQNcHZE\nXC3pg8BVwHLg7cDXImLxUG1oOWecDOFYAcwBpgNnSPrY4OsiojsiFlUedFKrtzMzG7aATIa2AUuA\nS4ELJO1NjuXAjcCFkp4A3pucA9wHPAkcBL4F/FOjG6Tpv78XeCoifgMg6S7gb4DvpKjTzCxD2aQp\nIuInVBbkrGVpjesDuHI490jTyv8FzpM0iUqaYinQk6I+M7NMFWkGXsvBOCJ2SboT2ENlQPTPge6s\nGmZmloXSB2OAiFgPFGdkuZmNKqOiZ2xmlneBOFaQ7aEdjM2stNwzNjPLAQdjM7OcaHIcccc5GJtZ\naWU5HXqkFaOVZmYtcJrCzCwHKqMpJnS6GU1xMDaz0nKawswsJ5ymMDPrMOeMzcxywMHYzCwHPB26\njoWT+uiZ9+pdYNVT/HWGfsVrd7Z9U8HXT4oldXaH/mmxn+v45Nc+14Tni/1M8e46v6sfF+e5av9e\n0tfrnrGZWU44GJuZddgwdofuOAdjMyutIo0zTrU79LBvpunx8kbRZmZDSr879PRFXfHJnk80vO7z\n+vfU90qrGP9lmJm1IBDHPR3azKyznDM2M8uBIuWMi9FKM7MWeWibmVmHFWnSx5hON8DMbKQM5Iwb\nHc2QdIukfkn7qsquk3RE0t7kWF713bWSDkp6XNL7G9XvnrGZlVZlNEVma1PcCnwduG1Q+U0R8cXq\nAknzgJXAW4HpwA8lvSkiTtar3D1jMyutgTRFo6OpuiIeBJ5t8tYrgNsj4lhEPAUcBBYP9QMOxmZW\nak0G46mSeqqO4cxOu0rSI0ka4/VJ2QzgcNU1vUlZXU5TmFlpDWOc8dEWZ+BtAD4PRPLvl4DGU/5q\ncDA2s9Ia6XHGEfHMwGdJ3wLuTU6PALOqLp2ZlNXlNIWZldbAdOhGR6skdVWdfhgYGGmxDVgpaaKk\nOcBc4OGh6nLP2MxKK8vp0JK2AOdTyS/3AuuB8yXNp5KmOARcARARj0raCuwHTgBXDjWSAhyMzazk\nskpTRMSqGsUbh7j+BuCGZut3MDaz0irSDDwHYzMrLQfjOhZO7aPnbwdtSPrN4myaWE8pNyS9rM4m\nl7cW/LlWvfa5tKXgz/S+Or+r+4v9XFnxEppmZh12ijFZToceUamGtkmaIulOSY9JOiDpHVk1zMws\nC1lNhx5paXvGXwW+HxEfkTQBmJRBm8zMMjEqcsaSJgPvAi4DiIjjwPFsmmVmll5QnJxxy7tDJwOd\nu6kMaj4H2A2sjYg/DLpuDS9vCT15IaxL0VwzGz3S7w59+qJ5Madnc8PrDmhBx3eHTpMzHgcsADZE\nxLnAH4BrBl8UEd0RsajyoM5imFn7ZLmE5khLkzPuBXojYldyfic1grGZWacE4liKtSfaqeVgHBFP\nSzos6c0R8TiwlErKwswsF0bT7tCfAjYnIymeBD6evklmZtnJSxqikVTBOCL2Ah1NepuZ1TMqhraZ\nmeVdIE6ecjA2M+uoOCWOvVSM6dAOxmZWWhHi5An3jM3MOitwMDYz67QIceLPDsZmZh0mTp0sRpgr\nRivNzFoRgNMUZmYddkrwUjHCXDFaaWbWqhOdbkBzHIzNrLwqCxoXgoOxmZWXg3FtC2f00fPpQbtD\nf7b4O9jG9TV2HF5f7Of6RY0drwHOKfqu10tq/K5+WuxnsiEE8OdsqpJ0C3AR0B8Rb0vKzgbuAGYD\nh4BLIuI5SaKyLd1y4I/AZRGxZ6j6U21IamaWawEca+Jozq3AskFl1wA7ImIusINX1nT/ADA3OdYA\nGxpV7mBsZuU1kKZodDRTVcSDwLODilcAm5LPm4CLq8pvi4qHgCmSuoaq3zljMyuv5nPGUyX1VJ13\nR0R3Ez83LSL6ks9PA9OSzzOAw1XX9SZlfdThYGxm5dV8MD6adkPSiAhJre3wTJuD8e4jXeiza9p5\ny7Yo+su6Wor+oq4ev6wbZUZ+NMUzkroioi9JQ/Qn5UeAWVXXzUzK6nLO2MzKLaOccR3bgNXJ59XA\nPVXl/6CK84Dnq9IZNTlNYWbldQp4KZuqJG0BzqeSX+4F1gM3AlslXQ78Grgkufw+KsPaDlIZ2tZw\nf1AHYzMrrwzTFBGxqs5XS2tcG8CVw6nfwdjMyssz8MzMcsDB2MwsJxyMzcw6zD1jM7McOAX8qdON\naI6DsZmVVwAnO92I5jgYm1m5OU1hZtZhzhmbmeWAg7GZWQ5kOB16pDkYm1m5uWdsZtZhTlOYmeVA\nhhuSjrT27g49vY+eKwbtDl2Chdnj7TV2HN5V7OeK5bV3h9Z9xX4uG2VG0zhjSWOBHuBIRFyUvklm\nZhkZZWmKtcAB4HUZ1GVmlp2gMNOhU227JGkm8EHg5myaY2aWoYE0RaMjB9L2jL8CXA2cVe8CSWuA\nZBfSyWj9upS3zJ+i54drcW7YSqFAaYqWe8aSLgL6I2L3UNdFRHdELKpsgz2p1duZmQ3fQDAeuQ1J\nM5OmZ7wE+JCk5cBpwOskfSciPpZN08zMUirQ0LaWe8YRcW1EzIyI2cBK4AEHYjPLnVGSMzYzy6/R\ntjZFRPwI+FEWdZmZZaZAaQr3jM2svEbTDDwzs1zLaLSEpEPAC1TC+4mIWCTpbOAOYDZwCLgkIp5r\npf5Ukz7MzHIt+6Ft74mI+ZWhugBcA+yIiLnAjuS8JQ7GZlZeAy/wGh2tWwFsSj5vAi5utSIHYzMr\nr+Z7xlMl9VQda+rUdr+k3VXfT4uIvuTz08C0VpvqnLGZlVtzaYijVamHet4ZEUck/SWwXdJj1V9G\nREiKFlvpYGxmJZbh0LaIOJL82y/pbmAx8Iykrojok9QF9Ldav9MUZlZeGa3aJukMSWcNfAbeB+wD\ntgGrk8tWA/e02lT3jM2svLJbtW0acLckqMTN70bE9yX9DNgq6XLg18Alrd7AwdjMyusUmSwuHxFP\nAufUKP8tsDT9HRyMzazsPAPPzCwHWh7f0F5tDcYLxvbx0Jmv3nV4wvPF31EiVtXYHXpL8Z/LzNrH\noynMzHLAwdjMLAecMzazEstoOEUbtDUY7znZxYTna035Ljbnh83yqjiry7tnbGYllt2sj5HmYGxm\nJeaesZlZDjgYm5nlQOAXeGZmHeecsZlZDjhNYWaWA+4Zm5nlgHvGZmY54J6xmVkOeDq0mVkOOE1h\nZpYTTlOYmXWYe8ZmZjngYGxmlgMeTWFmlgMeTWFmlgNOU9Q0X338eMKrd1KefMy7ZJjZSClOmqLl\nDUklzZK0U9J+SY9KWptlw8zM0hvoGTc6GpO0TNLjkg5KuibrlqbpGZ8APhMReySdBeyWtD0i9mfU\nNjOzlLLpGUsaC3wDuBDoBX4maVuW8a7lYBwRfUBf8vkFSQeAGYCDsZnlRGYv8BYDByPiSQBJtwMr\nyDDeZZIzljQbOBfYVeO7NcDAltDHJh9j36uvuH7wjxTRVOBopxuRsTI+E5Tzucr4TABvTl9F3w/g\nuqlNXHiapJ6q8+6I6K46nwEcrjrvBd6evn2vSB2MJZ0JfA9YFxG/H/x98kDdybU9EbEo7T3zpozP\nVcZngnI+VxmfCSrPlbaOiFiWRVvaoeUXeACSxlMJxJsj4q5smmRmljtHgFlV5zOTssykGU0hYCNw\nICK+nF2TzMxy52fAXElzJE0AVgLbsrxBmp7xEuBS4AJJe5NjeYOf6W7wfVGV8bnK+ExQzucq4zNB\njp4rIk4AVwE/AA4AWyPi0SzvoYjIsj4zM2tBqpyxmZllw8HYzCwH2hKMJd0iqV/SvsZXF0NZp4NL\nOk3Sw5J+kTxXKQaCQ2UWlaSfS7q3023JiqRDkn6ZvLNJPRQsDyRNkXSnpMckHZD0jk63qR3akjOW\n9C7gReC2iHjbiN+wDSR1AV3V08GBi4s+HTwZJXNGRLyYDF38CbA2Ih7qcNNSk/QvwCLgdRFxUafb\nkwVJh4BFEVGaSR+SNgH/FRE3JyMXJkXE7zrdrpHWlp5xRDwIPNuOe7VLRPRFxJ7k8wtU3rDO6Gyr\n0ouKF5PT8clR+Le8kmYCHwRu7nRbrD5Jk4F3URk2S0QcHw2BGJwzzsRQ08GLKPlzfi/QD2yPiDI8\n11eAq6ksVlAmAdwvaXey9EDRzQF+A/xnklK6WdIZnW5UOzgYp9RoOngRRcTJiJhPZZbRYkmFTi1J\nugjoj4jdnW7LCHhnRCwAPgBcmaQEi2wcsADYEBHnAn8AMl+uMo8cjFMo+3Tw5M/DnUBh5vfXsQT4\nUJJfvZ3KRKXvdLZJ2YiII8m//cDdVFYXK7JeoLfqr7E7qQTn0nMwblFZp4NLeoOkKcnn06ms3/pY\nZ1uVTkRcGxEzI2I2lWmsD0TExzrcrNQknZG8PCb5U/59MHhVxGKJiKeBw5IGVmxbyihZlrct2y5J\n2gKcD0yV1Ausj4iN7bj3CBqYDv7LJL8K8LmIuK+DbcpCF7ApWUx7DJVpn6UZClYy04C7K/0CxgHf\njYjvd7ZJmfgUsDkZSfEk8PEOt6ctPB3azCwHnKYwM8sBB2MzsxxwMDYzywEHYzOzHHAwNjPLAQdj\nM7MccDA2M8uB/weM5VAIUh+uMQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "x = []\n", "y = []\n", "\n", "for _ in range(10000):\n", " i, j = gibbs_sample()\n", " x.append(i)\n", " y.append(j)\n", " \n", "# plt.scatter(x, y)\n", "plt.hist2d(x, y, bins=50, cmap=plt.cm.jet)\n", "plt.colorbar()" ] }, { "cell_type": "code", "execution_count": 155, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "defaultdict(.>,\n", " {(1, 2): [256, 320],\n", " (1, 3): [288, 262],\n", " (1, 4): [260, 270],\n", " (1, 5): [306, 289],\n", " (1, 6): [277, 272],\n", " (1, 7): [260, 256],\n", " (2, 3): [291, 269],\n", " (2, 4): [259, 273],\n", " (2, 5): [293, 283],\n", " (2, 6): [292, 288],\n", " (2, 7): [289, 266],\n", " (2, 8): [264, 301],\n", " (3, 4): [260, 271],\n", " (3, 5): [264, 270],\n", " (3, 6): [279, 297],\n", " (3, 7): [274, 258],\n", " (3, 8): [292, 301],\n", " (3, 9): [277, 275],\n", " (4, 5): [296, 274],\n", " (4, 6): [283, 264],\n", " (4, 7): [258, 290],\n", " (4, 8): [257, 278],\n", " (4, 9): [259, 290],\n", " (4, 10): [290, 294],\n", " (5, 6): [269, 253],\n", " (5, 7): [276, 280],\n", " (5, 8): [267, 289],\n", " (5, 9): [286, 257],\n", " (5, 10): [289, 256],\n", " (5, 11): [312, 241],\n", " (6, 7): [276, 273],\n", " (6, 8): [306, 286],\n", " (6, 9): [262, 276],\n", " (6, 10): [291, 264],\n", " (6, 11): [283, 311],\n", " (6, 12): [259, 303]})" ] }, "execution_count": 155, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from collections import defaultdict\n", "def compare_distributions(num_samples=10000):\n", " counts = defaultdict(lambda: [0, 0])\n", " for _ in range(num_samples):\n", " counts[gibbs_sample()][0] += 1\n", " counts[direct_sample()][1] += 1\n", " return counts\n", "\n", "compare_distributions()" ] } ], "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.6.1" } }, "nbformat": 4, "nbformat_minor": 0 }