Skip to content

Instantly share code, notes, and snippets.

@AllenDowney
Last active November 1, 2018 14:01
Show Gist options
  • Save AllenDowney/1506edbddaf7f1280e80458e977e41e4 to your computer and use it in GitHub Desktop.
Save AllenDowney/1506edbddaf7f1280e80458e977e41e4 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Think Bayes\n",
"\n",
"This notebook presents code and exercises from Think Bayes, second edition.\n",
"\n",
"Copyright 2018 Allen B. Downey\n",
"\n",
"MIT License: https://opensource.org/licenses/MIT"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Configure Jupyter so figures appear in the notebook\n",
"%matplotlib inline\n",
"\n",
"# Configure Jupyter to display the assigned value after an assignment\n",
"%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
"\n",
"import math\n",
"import numpy as np\n",
"\n",
"from thinkbayes2 import Pmf, Cdf, Suite\n",
"import thinkplot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Game of Ur problem\n",
"\n",
"In the Royal Game of Ur, players advance tokens along a track with 14 spaces. To determine how many spaces to advance, a player rolls 4 dice with 4 sides. Two corners on each die are marked; the other two are not. The total number of marked corners -- which is 0, 1, 2, 3, or 4 -- is the number of spaces to advance.\n",
"\n",
"For example, if the total on your first roll is 2, you could advance a token to space 2. If you roll a 3 on the next roll, you could advance the same token to space 5.\n",
"\n",
"Suppose you have a token on space 13. How many rolls did it take to get there?\n",
"\n",
"Hint: you might want to start by computing the distribution of k given n, where k is the number of the space and n is the number of rolls.\n",
"\n",
"Then think about the prior distribution of n."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's a Pmf that represents one of the 4-sided dice."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Pmf({0: 0.5, 1: 0.5})"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"die = Pmf([0, 1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's the outcome of a single roll."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Pmf({0: 0.0625, 1: 0.25, 2: 0.375, 3: 0.25, 4: 0.0625})"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"roll = sum([die]*4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I'll start with a simulation, which helps in two ways: it makes modeling assumptions explicit and it provides an estimate of the answer.\n",
"\n",
"The following function simulates playing the game over and over; after every roll, it yields the number of rolls and the total so far. When it gets past the 14th space, it starts over."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def roll_until(iters):\n",
" \"\"\"Generates observations of the game.\n",
" \n",
" iters: number of observations\n",
" \n",
" yields: number of rolls, total\n",
" \"\"\"\n",
" for i in range(iters):\n",
" total = 0\n",
" for n in range(1, 1000):\n",
" total += roll.Random()\n",
" if total > 14:\n",
" break\n",
" yield(n, total)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now I'll the simulation many times and, every time the token is observed on space 13, record the number of rolls it took to get there."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"pmf_sim = Pmf()\n",
"for n, k in roll_until(1000000):\n",
" if k == 13:\n",
" pmf_sim[n] += 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the distribution of the number of rolls:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"501019"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf_sim.Normalize()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 0.01727080210530938\n",
"5 0.14913606070827654\n",
"6 0.29612849013710063\n",
"7 0.27942253686985924\n",
"8 0.16166652362485256\n",
"9 0.06677391476171562\n",
"10 0.021727718908863738\n",
"11 0.006071626026158689\n",
"12 0.001429087519634984\n",
"13 0.00031136543723890716\n",
"14 5.788203640979684e-05\n",
"15 3.991864579985989e-06\n"
]
}
],
"source": [
"pmf_sim.Print()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHQVJREFUeJzt3XuUVeWd5vHvQ3FT0QSFTlpLLFS8oCBqiSYmWknUkI4LdAZHGJOhZ+zF5IImccXRGBdJTHQR7WSMrRMlhsZJi2iMOtUZjEEj2okoN8sLKC0SRyowE4TEtMYb8Js/9i5yOJy6n815q3g+a9Wqs/d+3/f89imop/Y++7xbEYGZmVlqBtS6ADMzs0ocUGZmliQHlJmZJckBZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSBta6gGoZMWJENDQ01LoMMzNrx8qVK1+LiJFdbd9vAqqhoYEVK1bUugwzM2uHpP/TnfY+xWdmZklyQJmZWZIcUGZmlqRC34OSNAn4AVAH3B4Rc8q2fw74IrAdeAOYGRFr8m1fAy7Ot10aEQ8VWauZ7T3ee+89Wltbefvtt2tdSr80dOhQ6uvrGTRoUK/GKSygJNUBtwBnA63AcknNbQGUWxARt+btJwPfByZJGgtMA44DDgYelnRURGwvql4z23u0tray//7709DQgKRal9OvRARbtmyhtbWV0aNH92qsIk/xTQTWRcT6iHgXWAhMKW0QEX8qWdwPaLt74hRgYUS8ExG/Bdbl45mZ9drbb7/NQQcd5HAqgCQOOuigqhydFnmK7xBgQ8lyK3BqeSNJXwQuAwYDHy/p+2RZ30Mq9J0JzAQYNWpUVYo2s72Dw6k41XptizyCqlThbveXj4hbIuII4Arg6m72nRsRjRHROHJklz/7ZWZmfUCRR1CtwKEly/XAxg7aLwR+2MO+ZmY9Nus7d1V1vJuvnt6ldtdeey0LFiygrq6OAQMGcNttt/GjH/2Iyy67jLFjx/a6jrYJDEaMGNFum+uuu46rrrpq5/KHP/xhnnjiiV4/dzUUGVDLgTGSRgO/I7vo4T+WNpA0JiJeyhc/DbQ9bgYWSPo+2UUSY4BlBdba71X7P2AlXf1PaWawdOlSfv7zn7Nq1SqGDBnCa6+9xrvvvsvtt9++R+soD6hUwgkKPMUXEduAWcBDwAvAPRGxWtI1+RV7ALMkrZbUQvY+1Iy872rgHmAN8Avgi76Cz8z6k02bNjFixAiGDBkCwIgRIzj44INpamraOW3bsGHDuOKKKzj55JM566yzWLZsGU1NTRx++OE0NzcDMH/+fGbNmrVz3HPPPZclS5bs9nznnXceJ598Mscddxxz584F4Morr+Stt95iwoQJXHTRRTufE7Kr8S6//HKOP/54xo0bx9133w3AkiVLaGpqYurUqRxzzDFcdNFFROz2DkxVFPpB3YhYFBFHRcQREXFtvm52RDTnj78UEcdFxISI+FgeTG19r837HR0RDxZZp5nZnnbOOeewYcMGjjrqKL7whS/w2GOP7dbmzTffpKmpiZUrV7L//vtz9dVXs3jxYu6//35mz57dreebN28eK1euZMWKFdx0001s2bKFOXPmsM8++9DS0sKdd965S/v77ruPlpYWnnnmGR5++GEuv/xyNm3aBMDTTz/NjTfeyJo1a1i/fj2/+c1vev5CdMAzSZiZ1cCwYcNYuXIlc+fOZeTIkVx44YXMnz9/lzaDBw9m0qRJAIwbN44zzzyTQYMGMW7cOF555ZVuPd9NN93ECSecwGmnncaGDRt46aWXOmz/61//munTp1NXV8cHPvABzjzzTJYvXw7AxIkTqa+vZ8CAAUyYMKHbtXRVv5nN3Mysr6mrq6OpqYmmpibGjRvHHXfcscv2QYMG7bxke8CAATtPBw4YMIBt27YBMHDgQHbs2LGzT6XPHy1ZsoSHH36YpUuXsu+++9LU1NTp55Q6Om3XVkfbPrTVUm0+gjIzq4G1a9fuchTT0tLCYYcd1u1xGhoaaGlpYceOHWzYsIFly3a/nuz1119n+PDh7Lvvvrz44os8+eRfPmY6aNAg3nvvvd36nHHGGdx9991s376dzZs38/jjjzNx4p6dL8FHUGa216vFFahvvPEGl1xyCX/84x8ZOHAgRx55JHPnzmXq1KndGuf0009n9OjRjBs3juOPP56TTjpptzaTJk3i1ltvZfz48Rx99NGcdtppO7fNnDmT8ePHc9JJJ+3yPtT555/P0qVLOeGEE5DE9ddfzwc/+EFefPHFnu90N6moqy/2tMbGxvANC9vny8zN/uKFF17g2GOPrXUZ/Vql11jSyoho7OoYPsVnZmZJckCZmVmSHFBmtlfqL29vpKhar60Dysz2OkOHDmXLli0OqQK03Q9q6NChvR7LV/FZIXxRhqWsvr6e1tZWNm/eXOtS+qW2O+r2lgPKzPY6gwYN6vXdXq14PsVnZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZklyQJmZWZIcUGZmliQHlJmZJckBZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZkkqNKAkTZK0VtI6SVdW2H6ZpDWSnpX0iKTDSrZtl9SSfzUXWaeZmaWnsFu+S6oDbgHOBlqB5ZKaI2JNSbOngcaI+LOkzwPXAxfm296KiAlF1WdmZmkr8ghqIrAuItZHxLvAQmBKaYOIeDQi/pwvPgnUF1iPmZn1IUUG1CHAhpLl1nxdey4GHixZHipphaQnJZ1XqYOkmXmbFZs3b+59xWZmlozCTvEBqrAuKjaUPgM0AmeWrB4VERslHQ78StJzEfHyLoNFzAXmAjQ2NlYc28zM+qYij6BagUNLluuBjeWNJJ0FfB2YHBHvtK2PiI359/XAEuDEAms1M7PEFBlQy4ExkkZLGgxMA3a5Gk/SicBtZOH0+5L1wyUNyR+PAE4HSi+uMDOzfq6wU3wRsU3SLOAhoA6YFxGrJV0DrIiIZuAGYBjwU0kAr0bEZOBY4DZJO8hCdE7Z1X9mZtbPFfkeFBGxCFhUtm52yeOz2un3BDCuyNrMzCxtnknCzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0tSoQElaZKktZLWSbqywvbLJK2R9KykRyQdVrJthqSX8q8ZRdZpZmbpKSygJNUBtwCfAsYC0yWNLWv2NNAYEeOBe4Hr874HAt8ATgUmAt+QNLyoWs3MLD1FHkFNBNZFxPqIeBdYCEwpbRARj0bEn/PFJ4H6/PEngcURsTUi/gAsBiYVWKuZmSWmyIA6BNhQstyar2vPxcCDPexrZmb9zMACx1aFdVGxofQZoBE4szt9Jc0EZgKMGjWqZ1WamVmSijyCagUOLVmuBzaWN5J0FvB1YHJEvNOdvhExNyIaI6Jx5MiRVSvczMxqr8iAWg6MkTRa0mBgGtBc2kDSicBtZOH0+5JNDwHnSBqeXxxxTr7OzMz2EoWd4ouIbZJmkQVLHTAvIlZLugZYERHNwA3AMOCnkgBejYjJEbFV0rfJQg7gmojYWlStZmaWniLfgyIiFgGLytbNLnl8Vgd95wHziqvOzMxS5pkkzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkdBpSk+SWPZxRejZmZWa6zI6gTSh5/qchCzMzMSg3sZHvskSrMemHWd+4q/Dluvnp64c9hZrvqLKDqJd0EqOTxThFxaWGVmZnZXq2zgLq85PGKIgsxMzMr1WFARcQde6oQMzOzUh0GlKTmjrZHxOTqlmNmZpbp7BTfh4ANwF3AU2TvRZmZmRWus8vMPwhcBRwP/AA4G3gtIh6LiMc6G1zSJElrJa2TdGWF7WdIWiVpm6SpZdu2S2rJvzo8kjMzs/6nw4CKiO0R8YuImAGcBqwDlki6pLOBJdUBtwCfAsYC0yWNLWv2KvC3wIIKQ7wVERPyL59KNDPby3R2ig9JQ4BPA9OBBuAm4L4ujD0RWBcR6/NxFgJTgDVtDSLilXzbjm7WbWZm/VxnF0ncQXZ670HgWxHxfDfGPoTs/as2rcCp3eg/VNIKYBswJyIeqFDfTGAmwKhRo7oxtJmZpa6zI6jPAm8CRwFfktQ2s4SAiIgDOuhb6YKK7sxMMSoiNko6HPiVpOci4uVdBouYC8wFaGxs7HOzXhQ9A4JnPzCzvqyzz0H1ZrbzVuDQkuV6YGNXO0fExvz7eklLgBOBlzvsZGZm/UZns5kPlfRlSTdLmimp0/esSiwHxkgaLWkwMA3o0tV4kobn730haQRwOiXvXZmZWf/X2RHSHUAj8BzwN8D3ujpwRGwDZgEPAS8A90TEaknXSJoMIOkUSa3ABcBtklbn3Y8FVkh6BniU7D0oB5SZ2V6ksyOisRExDkDSj4Fl3Rk8IhYBi8rWzS55vJzs1F95vyeAcd15LjMz6186O4J6r+1BfkRkZma2R3R2BHWCpD/ljwXsky935So+MzOzHuvsKr66PVWImZlZqd5cRm5mZlYYB5SZmSXJAWVmZklyQJmZWZIcUGZmliQHlJmZJckBZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZklyQJmZWZIcUGZmliQHlJmZJckBZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpakQgNK0iRJayWtk3Rlhe1nSFolaZukqWXbZkh6Kf+aUWSdZmaWnsICSlIdcAvwKWAsMF3S2LJmrwJ/Cywo63sg8A3gVGAi8A1Jw4uq1czM0lPkEdREYF1ErI+Id4GFwJTSBhHxSkQ8C+wo6/tJYHFEbI2IPwCLgUkF1mpmZokpMqAOATaULLfm66rWV9JMSSskrdi8eXOPCzUzs/QUGVCqsC6q2Tci5kZEY0Q0jhw5slvFmZlZ2ooMqFbg0JLlemDjHuhrZmb9QJEBtRwYI2m0pMHANKC5i30fAs6RNDy/OOKcfJ2Zme0lCguoiNgGzCILlheAeyJitaRrJE0GkHSKpFbgAuA2SavzvluBb5OF3HLgmnydmZntJQYWOXhELAIWla2bXfJ4Odnpu0p95wHziqzPzMzS5ZkkzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJI0sNYFmPU1s75zV+HPcfPV0wt/DrPU+QjKzMyS5IAyM7MkOaDMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS1KhASVpkqS1ktZJurLC9iGS7s63PyWpIV/fIOktSS35161F1mlmZukpbCYJSXXALcDZQCuwXFJzRKwpaXYx8IeIOFLSNOC7wIX5tpcjYkJR9ZmZWdqKPIKaCKyLiPUR8S6wEJhS1mYKcEf++F7gE5JUYE1mZtZHFBlQhwAbSpZb83UV20TENuB14KB822hJT0t6TNJHC6zTzMwSVORksZWOhKKLbTYBoyJii6STgQckHRcRf9qlszQTmAkwatSoKpRsZmapKPIIqhU4tGS5HtjYXhtJA4H3AVsj4p2I2AIQESuBl4Gjyp8gIuZGRGNENI4cObKAXTAzs1opMqCWA2MkjZY0GJgGNJe1aQZm5I+nAr+KiJA0Mr/IAkmHA2OA9QXWamZmiSnsFF9EbJM0C3gIqAPmRcRqSdcAKyKiGfgx8BNJ64CtZCEGcAZwjaRtwHbgcxGxtahazcwsPYXesDAiFgGLytbNLnn8NnBBhX4/A35WZG1mZpY2zyRhZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZklyQJmZWZIcUGZmliQHlJmZJckBZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpYkB5SZmSWp0BsWmlnvzfrOXYU/x81XTy/8Ocy6y0dQZmaWJAeUmZklyaf4SvhUiplZOnwEZWZmSXJAmZlZkhxQZmaWJAeUmZklyQFlZmZJckCZmVmSHFBmZpYkfw7KzHbyZwEtJT6CMjOzJBUaUJImSVoraZ2kKytsHyLp7nz7U5IaSrZ9LV+/VtIni6zTzMzSU1hASaoDbgE+BYwFpksaW9bsYuAPEXEk8N+B7+Z9xwLTgOOAScD/yMczM7O9RJHvQU0E1kXEegBJC4EpwJqSNlOAb+aP7wVulqR8/cKIeAf4raR1+XhLC6zXzPYwv+dlHVFEFDOwNBWYFBF/ly9/Fjg1ImaVtHk+b9OaL78MnEoWWk9GxD/l638MPBgR95Y9x0xgZr54NLC2kJ1p3wjgtT38nHtSf94/71vf1Z/3r7/v234RMbKrHYo8glKFdeVp2F6brvQlIuYCc7tfWnVIWhERjbV6/qL15/3zvvVd/Xn/9oJ9a+hOnyIvkmgFDi1Zrgc2ttdG0kDgfcDWLvY1M7N+rMiAWg6MkTRa0mCyix6ay9o0AzPyx1OBX0V2zrEZmJZf5TcaGAMsK7BWMzNLTGGn+CJim6RZwENAHTAvIlZLugZYERHNwI+Bn+QXQWwlCzHydveQXVCxDfhiRGwvqtZeqNnpxT2kP++f963v6s/7530rUdhFEmZmZr3hmSTMzCxJDigzM0uSA6oXJNVJelrSz2tdSzVJer+keyW9KOkFSR+qdU3VIukrklZLel7SXZKG1rqm3pA0T9Lv888Utq07UNJiSS/l34fXssaeamffbsj/XT4r6X5J769ljb1Raf9Ktn1VUkgaUYvaequ9fZN0ST593WpJ13c2jgOqd74EvFDrIgrwA+AXEXEMcAL9ZB8lHQJcCjRGxPFkF+9Mq21VvTafbDqwUlcCj0TEGOCRfLkvms/u+7YYOD4ixgP/CnxtTxdVRfPZff+QdChwNvDqni6oiuZTtm+SPkY2S9D4iDgO+PvOBnFA9ZCkeuDTwO21rqWaJB0AnEF2hSUR8W5E/LG2VVXVQGCf/HN3+9LHP18XEY+TXQFbagpwR/74DuC8PVpUlVTat4j4ZURsyxefJPuMZJ/Uzs8OsnlJ/xsVJifoK9rZt88Dc/Ip7IiI33c2jgOq524k+0e0o9aFVNnhwGbgH/PTl7dL2q/WRVVDRPyO7K+2V4FNwOsR8cvaVlWID0TEJoD8+1/VuJ6i/BfgwVoXUU2SJgO/i4hnal1LAY4CPprfueIxSad01sEB1QOSzgV+HxEra11LAQYCJwE/jIgTgTfpu6eIdpG/FzMFGA0cDOwn6TO1rcp6QtLXyT4jeWeta6kWSfsCXwdm17qWggwEhgOnAZcD9+STg7fLAdUzpwOTJb0CLAQ+LumfaltS1bQCrRHxVL58L1lg9QdnAb+NiM0R8R5wH/DhGtdUhP8n6a8B8u+dnkrpSyTNAM4FLor+9UHOI8j+eHom/91SD6yS9MGaVlU9rcB9kVlGdvapw4tAHFA9EBFfi4j6fOLDaWRTNPWLv8Qj4v8CGyQdna/6BLveIqUvexU4TdK++V9un6CfXABSpnQKsRnA/6phLVUlaRJwBTA5Iv5c63qqKSKei4i/ioiG/HdLK3BS/n+yP3gA+DiApKOAwXQyc7sDyiq5BLhT0rPABOC6GtdTFflR4b3AKuA5sn//fXpqGUl3kd0n7WhJrZIuBuYAZ0t6iexqsDm1rLGn2tm3m4H9gcWSWiTdWtMie6Gd/esX2tm3ecDh+aXnC4EZnR0Be6ojMzNLko+gzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDiirmnz25e+VLH9V0jerNPZ8SVOrMVYnz3NBPoP7o70Y4438e0Olmao76HdVT5+zG8/RJKnih5M72payDmbO/nY+63mLpF9KOrhWNVrPOKCsmt4B/l1qtwiQVNeN5hcDX4iIj1VpvO4oPKCAJtqfPaOjbSmbT4VZwYEbImJ8REwAfk7/nUKo33JAWTVtI/vg61fKN5QfAZUcZTTlE0feI+lfJc2RdJGkZZKek3REyTBnSfqXvN25ef+6/B5By/O/lv9rybiPSlpA9qHc8nqm5+M/L+m7+brZwEeAWyXdUNZ+t/EkXZb3f17Slzt6YSQdl+9TS17nmLLtc8hmWW+RdGdXx5d0cf56LJH0I0k35+tHSvpZ/rosl3S6pAbgc8BX8uf5aMk4u22TdJikR/J6H5E0qsLzn5m3b1E2ufD++Wv1uLL7Na2RdKukAXn7H0paoex+QN8qGecUSU9IeiZ/nfZv72dbrr1ZwSPiTyWL+9GHZwffa0WEv/xVlS/gDeAA4BXgfcBXgW/m2+YDU0vb5t+bgD8Cfw0MAX4HfCvf9iXgxpL+vyD7o2oM2TQwQ4GZwNV5myHACrL5zJrIJrodXaHOg8mmPRpJNoHlr4Dz8m1LyO4XVd5nl/GAk8mCaj9gGLAaOLFs3xqA5/PH/0A2dxxkU7zsU+n1K3nc7vhl+/EKcCAwCPgX4OZ82wLgI/njUcAL+eNvAl9t5+e3yzbgn8k+7Q/ZzOEPVOjzz8Dp+eNh+evZBLxNNjN+Hdk9nKbmbQ7Mv9flr/X4/PVYD5ySbzsgH6fiz7ad2ne+1mXrrwU2AM8DI2v9f8Rf3fvyEZRVVWR/tf5PshsDdtXyiNgU2X1iXgbaboHxHNkvnjb3RMSOiHiJ7BfaMcA5wH+S1AI8BRxEFmAAyyLitxWe7xRgSWSTxrbNiH1GF+osHe8jwP0R8WZEvEE28exH2+/KUuAqSVcAh0XEW508V1fGnwg8FhFbI5v89qcl284Cbs5fl2bgAEn7d2EfS32ILOgAfpLXVO43wPclXQq8P/5yr6ZlEbE+IrYDd5X0/Q+SVgFPA8cBY4GjgU0RsRyyf0P5OB39bLskIr4eEYeS/Yxndaev1Z4DyopwI9l7OaX3kdpG/u9Nksj+am7zTsnjHSXLO8j+km5TfoomAAGXRMSE/Gt0/OUeT2+2U1+HU/x3oHS8bo0REQuAycBbwEOSPt5Jl66M31GbAcCHSl6XQyLi37pYbnt2O0UWEXOAvwP2AZ6UdEw7bUPSaLKj6k9Edkfc/012FKxKY9Pxz7a7FgD/vod9rUYcUFZ1EbEVuIcspNq8QnbaCrJ7Mg3qwdAXSBqQvy91OLAWeAj4vKRBkM2SrM5vsPgUcKakEfkFD9OBx7pZy+PAecpmRt8POJ/sFFtFkg4H1kfETWRHNOMrNHuvbT+6OP6yfD+GK7tDcOkv4F9ScsQgaUL+8N/IJlutpHzbE2Sz9QNcBPy6wn4dEdks3N8lOwXXFlATJY3O33u6MO97AFnIvy7pA8Cn8rYvAgcrv4Fd/v7TQHr2sy2trfRoa3L+PNaHOKCsKN9j13u9/Ijsl+ky4FTaP7rpyFqyIHkQ+FxEvA3cTnY7kFXKLjO+jV2PunYT2V1mvwY8CjwDrIqIbt2SIiJWkb0vtows8G6PiKc76HIh8Hx+uuoYstOg5eYCz0q6syvjR3aH4Ovy7Q+TvQ6v55svBRrziwvWkF0AAdl7RueXXyTRzrZLgf+sbFb7z5K9J1juy/lFHM+QHR223eF2Kdks6s8DvyU7XfkM2am91WQzW/8m349389fnH/JxFpMdWXXpZ6v2ZwWfk9f2LNnpwkr1W8I8m7lZHyZpWES8kR9x3A/Mi4j7a1xTE9nFFufWsg7r+3wEZda3fTM/Kms7UnmgxvWYVY2PoMzMLEk+gjIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJP1/IORnQlZsL9AAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"thinkplot.Hist(pmf_sim, label='Simulation')\n",
"thinkplot.decorate(xlabel='Number of rolls to get to space 13',\n",
" ylabel='PMF')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bayes\n",
"\n",
"Now let's think about a Bayesian solution. It is straight forward to compute the likelihood function, which is the probability of being on space 13 after a hypothetical `n` rolls.\n",
"\n",
"`pmf_n` is the distribution of spaces after `n` rolls.\n",
"\n",
"`pmf_13` is the probability of being on space 13 after `n` rolls."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 0.008544921875\n",
"5 0.0739288330078125\n",
"6 0.14878177642822266\n",
"7 0.13948291540145874\n",
"8 0.08087921887636185\n",
"9 0.033626414369791746\n",
"10 0.010944152454612777\n",
"11 0.002951056014353526\n",
"12 0.0006854188303009323\n",
"13 0.00014100133496341982\n",
"14 2.6227807875534026e-05\n"
]
},
{
"data": {
"text/plain": [
"0.4999919364007537"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf_13 = Pmf()\n",
"for n in range(4, 15):\n",
" pmf_n = sum([roll]*n)\n",
" pmf_13[n] = pmf_n[13]\n",
" \n",
"pmf_13.Print()\n",
"pmf_13.Total()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total probability of the data is very close to 1/2, but it's not obvious (to me) why.\n",
"\n",
"Nevertheless, `pmf_13` is the probability of the data for each hypothetical values of `n`, so it is the likelihood function.\n",
"\n",
"### The prior\n",
"\n",
"Now we need to think about a prior distribution on the number of rolls. This is not easy to reason about, so let's start by assuming that it is uniform, and see where that gets us.\n",
"\n",
"If the prior is uniform, the posterior equals the likelihood function, normalized."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 0.017090119365747274\n",
"5 0.1478600505840099\n",
"6 0.2975683518003199\n",
"7 0.2789703298127999\n",
"8 0.16176104650522904\n",
"9 0.0672539133567936\n",
"10 0.021888657911956433\n",
"11 0.005902207214774348\n",
"12 0.0013708597687294604\n",
"13 0.00028200721791321923\n",
"14 5.245646172683854e-05\n"
]
}
],
"source": [
"posterior = pmf_13.Copy()\n",
"posterior.Normalize()\n",
"posterior.Print()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That sure looks similar to what we got by simulation. Let's compare them."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4FeW59/HvTUAQ8YCG3SoJBneDFgwgRrRKgW3R4hZB+qqAtaJY87aAh03BarVI2dXLIm0tYiuoiG1FVESKCrWgorWikmAURBFElAi7TUHpi4KScO8/1iTvIuScNcmT8PtcV67M4ZmZe8Lhl5k18zzm7oiIiISmVVMXICIiUhkFlIiIBEkBJSIiQVJAiYhIkBRQIiISJAWUiIgESQElIiJBUkCJiEiQFFAiIhKk1k1dQKqkp6d7VlZWU5chIiJVKCgo+Ke7d6pt+xYTUFlZWeTn5zd1GSIiUgUz+7Au7XWLT0REgqSAEhGRICmgREQkSC3mMygRqb+9e/dSVFTEnj17mroUaQHatWtHRkYGbdq0adB+Yg0oMxsM/AZIA+539zsqrP8BMA4oBXYBee6+Llp3E3BVtO5ad382zlpFDmZFRUUcfvjhZGVlYWZNXY40Y+7O9u3bKSoqomvXrg3aV2y3+MwsDbgHOA/oDowys+4Vms1z9xx37w1MA34VbdsdGAn0AAYDv432JyIx2LNnD8ccc4zCSRrMzDjmmGNScjUe52dQfYGN7r7J3b8E5gPDkhu4+7+SZg8Dyob3HQbMd/cv3P0DYGO0PxGJicJJUiVVf5fivMXXGdiSNF8EnF6xkZmNAyYAhwBnJ237aoVtO1eybR6QB9ClS5eUFC0iImGIM6Aqi1A/YIH7PcA9ZnYpcAswug7bzgZmA+Tm5h6wXiJbnoTN86C0ET8AT2sHWZdC5vDGO6Y0a2bGhAkT+OUvfwnA9OnT2bVrF1OmTGm0Gq644gqGDBnCRRddxPe//30mTJhA9+4VP5movc2bNzNkyBDWrl1b5fL8/Hx+//vfM2PGDKZMmUKHDh2YOHFirfaf3H7y5Mn079+fQYMGlXdckJ6eXu/aKxPXfqsS5y2+IiAzaT4D2FpN+/nAhfXcVqrT2OEEieNtnte4x5RmrW3btixcuJB//vOf9dq+pKQkpfXcf//9DQqn2srNzWXGjBkN3s/UqVMZNGhQCioKR5xXUKuAbDPrCnxM4qGHS5MbmFm2u2+IZs8HyqYXA/PM7FfAcUA28HqMtbZsUTit2fBx7IfKyU66E9vYoSipseKC+PY98KkqV7Vu3Zq8vDx+/etfc9ttt+237sMPP2TMmDEUFxfTqVMnHnzwQbp06cIVV1zB0UcfzRtvvEGfPn04/PDD+eCDD9i2bRvvvfcev/rVr3j11VdZunQpnTt35qmnnqJNmzZMnTqVp556it27d3PmmWcya9asAz43GThwINOnT2fr1q1MnjwZgN27d/Pll1/ywQcfUFBQwIQJE9i1axfp6enMnTuXY489loKCAsaMGUP79u3p169fjT+SFStWMH36dJ5++un9lt93330sXLiQhQsXsnXrVsaNG0dxcTHt27fnvvvu46STTtqvffLVH8Ddd9/NU089xd69e3n88cc56aST2LFjB2PGjGHTpk20b9+e2bNn07NnzyqXb9++nVGjRlFcXEzfvn1xT9yo+uyzz7jkkksoKiqitLSUn/70p4wYMaLGc62r2ALK3UvMbDzwLInHzOe4+9tmNhXId/fFwHgzGwTsBT4hcXuPqN1jwDqgBBjn7qVx1XowmfX3SbHuf+bVo+L9D05atHHjxtGzZ09uuOGG/ZaPHz+eyy+/nNGjRzNnzhyuvfZaFi1aBMB7773H8uXLSUtLY8qUKbz//vu88MILrFu3jm984xs88cQTTJs2jeHDh/PMM89w4YUXMn78+PLQ+d73vsfTTz/NBRdU/vd26NChDB06FIBLLrmEAQMGsHfvXq655hr+9Kc/0alTJx599FFuvvlm5syZw5VXXsndd9/NgAEDmDSpfv/eZs6cyV/+8hcWLVpE27ZtycvL49577yU7O5vXXnuNsWPH8vzzz1e7j/T0dFavXs1vf/tbpk+fzv3338+tt97KKaecwqJFi3j++ee5/PLLKSwsrHL5z372M/r168fkyZN55plnmD17NgB//vOfOe6443jmmWcA2LlzZ73Osyaxvgfl7kuAJRWWTU6avq6abW8DbqtqvYi0PEcccQSXX345M2bM4NBDDy1fvnLlShYuXAgkAiU5wC6++GLS0v7/WyjnnXcebdq0IScnh9LSUgYPHgxATk4OmzdvBuCFF15g2rRpfP755+zYsYMePXpUGVBlpk2bxqGHHsq4ceNYu3Yta9eu5ZxzzgGgtLSUY489lp07d/Lpp58yYMCA8lqXLl1ap5/BH/7wBzIyMli0aBFt2rRh165dvPLKK1x88cXlbb744osa9/Od73wHgFNPPbX8Z/fyyy/zxBNPAHD22Wezfft2du7cWeXyl156qXzb888/n44dOwKJn+XEiRP58Y9/zJAhQ/jmN79Zp3OsLfUkISL7q+Y2XGO4/vrr6dOnD1deeWWVbZJvxx122GH7rWvbti0ArVq1ok2bNuVtW7VqRUlJCXv27GHs2LHk5+eTmZnJlClTanxn57nnnuPxxx/npZdeAhIvo/bo0YOVK1fu1+7TTz9t8CPWJ598MoWFheUvuu7bt4+jjjqKwsLCOu2n7OeQlpZW/vlc2S26ZGZW5fLk78m6detGQUEBS5Ys4aabbuLcc88tvyJNJfXFJyJBOfroo7nkkkt44IEHypedeeaZzJ8/H4CHH364Vp/tVKUsjNLT09m1axcLFiyotv2HH37I2LFjeeyxx8qv6k488USKi4vLA2rv3r28/fbbHHXUURx55JG8/PLL5bXW1SmnnMKsWbMYOnQoW7du5YgjjqBr1648/vjjQCJk3nzzzTrvF6B///7lNa1YsYL09HSOOOKIWi1funQpn3zyCQBbt26lffv2XHbZZUycOJHVq1fXq56a6ApKRILzox/9iJkzZ5bPz5gxgzFjxnDnnXeWPyRRX0cddRRXX301OTk5ZGVlcdppp1Xbfu7cuWzfvp3hwxOvTBx33HEsWbKEBQsWcO2117Jz505KSkq4/vrr6dGjBw8++GD5QxLf/va361Vjv379mD59Oueffz7Lli3j4Ycf5oc//CE///nP2bt3LyNHjqRXr1513u+UKVO48sor6dmzJ+3bt+ehhx6qdvmtt97KqFGj6NOnDwMGDCh/33TNmjVMmjSp/Cr1d7/7Xb3OsyZW2aVdc5Sbm+sasLAK0UMLazZ8HP9DErdUeEiiiW8XSe288847fP3rX2/qMqQFqezvlJkVuHtubfehW3wiIhIkBZSIiARJASUiIkFSQImISJAUUCIiEiQ9Zi6xWbPhY9jQJ9Zj5GR3Vs/pIi2UAkpSL61d43YUW9ZzugIqZcb//JGU7m/mLaNq1e62225j3rx5pKWl0apVK2bNmsV9993X4GEvytRmuIjbb7+dn/zkJ+XzZ555Jq+88kqDjy11p1t8knpZlyZCqjGp5/Rmb+XKlTz99NOsXr2at956i+XLl5OZmdlow16Uuf322/ebVzg1HV1BSeplDofM4cx6MbW/hVdmZrbGnGoptm3bRnp6enkfcmVXOWXDXuTm5tKhQwfGjRvH8uXL6dixI7fffjs33HADH330EXfddRdDhw5l7ty55Ofnl/dEMWTIECZOnMjAgQP3O96FF17Ili1b2LNnD9dddx15eXnceOON7N69m969e9OjRw8efvhhOnTowK5du3B3brjhBpYuXYqZccsttzBixAhWrFjBlClTSE9PZ+3atZx66qn88Y9/TNmw5wczXUGJSBDOPfdctmzZQrdu3Rg7diwvvvjiAW0+++wzBg4cSEFBAYcffji33HILy5Yt48knn6xzZ6Vz5syhoKCA/Px8ZsyYwfbt27njjjs49NBDKSwsPKAfvYULF1JYWMibb77J8uXLmTRpEtu2bQPgjTfe4K677mLdunVs2rSJv/3tb/X/QUg5BZSIBKFDhw4UFBQwe/ZsOnXqxIgRI5g7d+5+bQ455JD9hs8YMGBA+dAaZUNp1NaMGTPo1asXZ5xxBlu2bGHDhg3Vtn/55ZcZNWoUaWlpfOUrX2HAgAGsWrUKgL59+5KRkUGrVq3o3bt3nWuRyukWn4gEIy0tjYEDBzJw4EBycnLKOy0tU3H4jOShNcqGlGjdujX79u0r36ayoTRWrFjB8uXLWblyJe3bt2fgwIE1DrlRXb+lZXWUnUOqh58/WOkKSkSCsH79+v2uYgoLCzn++OPrvJ+srCwKCwvZt28fW7Zs4fXXXz+gzc6dO+nYsSPt27fn3Xff5dVXXy1f16ZNG/bu3XvANv379+fRRx+ltLSU4uJiXnrpJfr27Vvn+qT2dAUlIgeo7WPhqbRr1y6uueYaPv30U1q3bs3XvvY1Zs+ezUUXXVSn/Zx11ll07dqVnJwcTj75ZPr0OfBdvMGDB3PvvffSs2dPTjzxRM4444zydXl5efTs2ZM+ffrs9znU8OHDWblyJb169cLMmDZtGl/96ld5991363/SUi0Nt3EwaOzhNiKpfpem0uP1S3qKT0N71JuG25BU03AbIiLSYimgREQkSAooEQGqf0pNpC5S9XdJASUitGvXju3btyukpMHcne3bt9OuXcO7O9NTfCJCRkYGRUVFFBcXN3Up0gK0a9eOjIyMBu9HASUitGnThq5duzZ1GSL70S0+EREJkgJKRESCpIASEZEgxRpQZjbYzNab2UYzu7GS9RPMbJ2ZvWVmz5nZ8UnrSs2sMPpaHGedIiISntgekjCzNOAe4BygCFhlZovdfV1SszeAXHf/3Mx+CEwDRkTrdrt777jqExGRsMV5BdUX2Ojum9z9S2A+MCy5gbu/4O6fR7OvAg1/LlFERFqEOAOqM7Alab4oWlaVq4ClSfPtzCzfzF41swsr28DM8qI2+Xp/Q0SkZYnzPSirZFmlr6mb2WVALjAgaXEXd99qZicAz5vZGnd/f7+duc8GZkOiN/PUlC0iIiGI8wqqCMhMms8AtlZsZGaDgJuBoe7+Rdlyd98afd8ErABOibFWEREJTJwBtQrINrOuZnYIMBLY72k8MzsFmEUinP6RtLyjmbWNptOBs4DkhytERKSFi+0Wn7uXmNl44FkgDZjj7m+b2VQg390XA3cCHYDHzQzgI3cfCnwdmGVm+0iE6B0Vnv4TEZEWLta++Nx9CbCkwrLJSdODqtjuFSAnztpERCRs6klCRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQlSrAFlZoPNbL2ZbTSzGytZP8HM1pnZW2b2nJkdn7RutJltiL5Gx1mniIiEJ7aAMrM04B7gPKA7MMrMuldo9gaQ6+49gQXAtGjbo4FbgdOBvsCtZtYxrlpFRCQ8cV5B9QU2uvsmd/8SmA8MS27g7i+4++fR7KtARjT9bWCZu+9w90+AZcDgGGsVEZHAxBlQnYEtSfNF0bKqXAUsrcu2ZpZnZvlmll9cXNzAckVEJCRxBpRVsswrbWh2GZAL3FmXbd19trvnuntup06d6l2oiIiEJ86AKgIyk+YzgK0VG5nZIOBmYKi7f1GXbUVEpOWKM6BWAdlm1tXMDgFGAouTG5jZKcAsEuH0j6RVzwLnmlnH6OGIc6NlIiJykGgd147dvcTMxpMIljRgjru/bWZTgXx3X0zill4H4HEzA/jI3Ye6+w4z+28SIQcw1d13xFWriIiEJ7aAAnD3JcCSCssmJ00PqmbbOcCc+KoTEZGQqScJEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEjVBpSZzU2aHh17NSIiIpGarqB6JU1fF2chIiIiyVrXsN4bpQqRVFhxQeMcJ60dZF0KmcMb53giB6maAirDzGYAljRdzt2vja0ykdpIawele1iz4ePYD5WT3TkxUboHNs9TQInErKaAmpQ0nR9nISL1knVpIiwaW+mexj+myEGm2oBy94caqxCReskcDpnDmfXiI7EfaubVoxrvNqKIVB9QZra4uvXuPjS15YiIiCTUdIvvG8AW4BHgNRKfRYmIiMSupsfMvwr8BDgZ+A1wDvBPd3/R3V+saedmNtjM1pvZRjO7sZL1/c1stZmVmNlFFdaVmllh9FXtlZyIiLQ81QaUu5e6+5/dfTRwBrARWGFm19S0YzNLA+4BzgO6A6PMrHuFZh8BVwCVfcq92917R1+6lSgicpCp6RYfZtYWOB8YBWQBM4CFtdh3X2Cju2+K9jMfGAasK2vg7pujdfvqWLeIiLRwNT0k8RCJ23tLgZ+5+9o67Lszic+vyhQBp9dh+3Zmlg+UAHe4+6JK6ssD8gC6dOlSh12LiEjoarqC+h7wGdANuM7MynqWMMDd/Yhqtq3sgYq69EzRxd23mtkJwPNmtsbd399vZ+6zgdkAubm5zafXiy1PJt7daewXTEVEmpGa3oNqSG/nRUBm0nwGsLW2G7v71uj7JjNbAZwCvF/tRs1FFE6Nba+3afRjiojUV029mbczs+vNbKaZ5ZlZjZ9ZJVkFZJtZVzM7BBgJ1OppPDPrGH32hZmlA2eR9NlVs9cUvRCktSN/11mNf1wRkXqqKXAeAvYCfwX+E+hBLXs1d/cSMxsPPAukAXPc/W0zmwrku/tiMzsNeBLoCFxgZj9z9x7A14FZ0cMTrUh8BtVyAirJrL9PqrlRA828ehQAbzVCbwsiIqlSU0B1d/ccADN7AHi9Ljt39yXAkgrLJidNryJx66/idq8AOXU5loiItCw1fca0t2zC3UtirkVERKRcTVdQvczsX9G0AYdG87V5ik9ERKTeanqKL62xChEREUnWkMfIRUREYqOAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREghRrQJnZYDNbb2YbzezGStb3N7PVZlZiZhdVWDfazDZEX6PjrFNERMITW0CZWRpwD3Ae0B0YZWbdKzT7CLgCmFdh26OBW4HTgb7ArWbWMa5aRUQkPHFeQfUFNrr7Jnf/EpgPDEtu4O6b3f0tYF+Fbb8NLHP3He7+CbAMGBxjrSIiEpg4A6ozsCVpvihaFve2IiLSAsQZUFbJMk/ltmaWZ2b5ZpZfXFxcp+JERCRscQZUEZCZNJ8BbE3ltu4+291z3T23U6dO9S5URETCE2dArQKyzayrmR0CjAQW13LbZ4Fzzaxj9HDEudEyERE5SMQWUO5eAownESzvAI+5+9tmNtXMhgKY2WlmVgRcDMwys7ejbXcA/00i5FYBU6NlIiJykGgd587dfQmwpMKyyUnTq0jcvqts2znAnDjrExGRcKknCRERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCFOuQ7yIt2ooL4j9GWjvIuhQyh8d/LJHAKKBE6iKtHZTuYc2Gj2M/VE52ZyjdA5vnKaDkoKRbfCJ1kXVpIqQaU+mexj2eSCB0BSVSF5nDIXM4s158JPZDzcyeF/sxREKmKygREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREghRrQJnZYDNbb2YbzezGSta3NbNHo/WvmVlWtDzLzHabWWH0dW+cdYqISHhiew/KzNKAe4BzgCJglZktdvd1Sc2uAj5x96+Z2UjgF8CIaN377t47rvpERCRscV5B9QU2uvsmd/8SmA8Mq9BmGPBQNL0A+JaZWYw1iYhIMxFnQHUGtiTNF0XLKm3j7iXATuCYaF1XM3vDzF40s2/GWKeIiAQozq6OKrsS8lq22QZ0cfftZnYqsMjMerj7v/bb2CwPyAPo0qVLCkoWEZFQxHkFVQRkJs1nAFuramNmrYEjgR3u/oW7bwdw9wLgfaBbxQO4+2x3z3X33E6dOsVwCiIi0lTiDKhVQLaZdTWzQ4CRwOIKbRYDo6Ppi4Dn3d3NrFP0kAVmdgKQDWyKsVYREQlMbLf43L3EzMYDzwJpwBx3f9vMpgL57r4YeAD4g5ltBHaQCDGA/sBUMysBSoEfuPuOuGoVEZHwxDrchrsvAZZUWDY5aXoPcHEl2z0BPBFnbSIiEjb1JCEiIkFSQImISJAUUCIiEiQFlIiIBEkBJSIiQVJAiYhIkBRQIiISJAWUiIgESQElIiJBUkCJiEiQFFAiIhIkBZSIiARJASUiIkFSQImISJAUUCIiEqRYx4MSkRRZcUHjHCetHWRdCpnDG+d4ItVQQImEKq0dlO4BYM2Gj2M9VE5258RE6R7YPE8BJUHQLT6RUGVdmgipxhaFokhT0xVUmS1PJn5z1D9OCUXm8PIrmVkvPxLroWZeParxbiOK1JICqkwUTnHfSoGk2ylN8duxiEgzoVt8ZRr7yqnsw2gREamUrqAqMevvk2Ld/8yrRyXNxXvrRkSkudIVlIiIBEkBJSIiQVJAiYhIkBRQIiISJD0kISIHaox3otStktRAASUiCVHXSo32LqC6VZIa6BafiCQ0RddK6rlFqhHrFZSZDQZ+A6QB97v7HRXWtwV+D5wKbAdGuPvmaN1NwFVAKXCtuz8bZ60iB72oa6VZL8b/bt7M7HmxH0Oav9gCyszSgHuAc4AiYJWZLXb3dUnNrgI+cfevmdlI4BfACDPrDowEegDHAcvNrJu7l8ZVr4g0EQ0lIlWI8wqqL7DR3TcBmNl8YBiQHFDDgCnR9AJgpplZtHy+u38BfGBmG6P9rYyxXhFpLE01lMj7cxJfLdXAp5q6gpQyd49nx2YXAYPd/fvR/PeA0919fFKbtVGbomj+feB0EqH1qrv/MVr+ALDU3RdUOEYekBfNngisj+VkapYO/LOJjh03nVvzpHNrnlr6uR3m7p1qu0GcV1BWybKKaVhVm9psi7vPBmbXvbTUMrN8d89t6jrioHNrnnRuzdNBcG5Zddkmzqf4ioDMpPkMYGtVbcysNXAksKOW24qISAsWZ0CtArLNrKuZHULioYfFFdosBkZH0xcBz3vinuNiYKSZtTWzrkA28HqMtYqISGBiu8Xn7iVmNh54lsRj5nPc/W0zmwrku/ti4AHgD9FDEDtIhBhRu8dIPFBRAowL/Am+Jr/NGCOdW/Okc2uedG5JYntIQkREpCHUk4SIiARJASUiIkFSQDWQmaWZ2Rtm9nRT15JKZnaUmS0ws3fN7B0z+0ZT15QqZvZfZva2ma01s0fMrJE7oEstM5tjZv+I3issW3a0mS0zsw3R945NWWN9VXFud0Z/L98ysyfN7KimrLG+Kju3pHUTzczNLL0pamuoqs7NzK4xs/XRv79pNe1HAdVw1wHvNHURMfgN8Gd3PwnoRQs5RzPrDFwL5Lr7ySQe4BnZtFU12FxgcIVlNwLPuXs28Fw03xzN5cBzWwac7O49gfeAmxq7qBSZy4Hnhpllkugi7qPGLiiF5lLh3MzsP0j0EtTT3XsA02vaiQKqAcwsAzgfuL+pa0klMzsC6E/iKUvc/Ut3/7Rpq0qp1sCh0bt37Wnm79i5+0sknoJNNgx4KJp+CLiwUYtKkcrOzd3/4u4l0eyrJN6TbHaq+HMD+DVwA5V0TtBcVHFuPwTuiLqww93/UdN+FFANcxeJv0j7mrqQFDsBKAYejG5f3m9mhzV1Uang7h+T+M3tI2AbsNPd/9K0VcXiK+6+DSD6/m9NXE9cxgBLm7qIVDGzocDH7v5mU9cSg27AN83sNTN70cxOq2kDBVQ9mdkQ4B/uXtDUtcSgNdAH+J27nwJ8RvO9RbSf6LOYYUBXEj3lH2ZmlzVtVVIfZnYzifckH27qWlLBzNoDNwOTm7qWmLQGOgJnAJOAx6LOwaukgKq/s4ChZrYZmA+cbWZ/bNqSUqYIKHL316L5BSQCqyUYBHzg7sXuvhdYCJzZxDXF4e9mdixA9L3G2ynNiZmNBoYA3/WW8zLnv5P4xenN6P+VDGC1mX21SatKnSJgoSe8TuLOU7UPgSig6sndb3L3jKjzw5EkumlqEb+Ju/v/AFvM7MRo0bfYf5iU5uwj4Awzax/99vYtWsgDIBUkdyM2GvhTE9aSUtFAqD8Ghrr7501dT6q4+xp3/zd3z4r+XykC+kT/HluCRcDZAGbWDTiEGnpuV0BJVa4BHjazt4DewO1NXE9KRFeFC4DVwBoS/waadfcyZvYIibHSTjSzIjO7CrgDOMfMNpB4IuyO6vYRqirObSZwOLDMzArN7N4mLbKeqji3FqGKc5sDnBA9ej4fGF3T1a+6OhIRkSDpCkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAkhYv6hX6l0nzE81sSor2PdfMLkrFvmo4zsVRr/IvNGAfu6LvWZX1oF3NdmZmz0d9NFbXbrqZnV3f+kQqUkDJweAL4DuhDV1gZml1aH4VMNbd/yNF+6uL/wTedPd/1dDublpIl1gSBgWUHAxKSLyM+18VV1S8Akq6yhgYdWj5mJm9Z2Z3mNl3zex1M1tjZv+etJtBZvbXqN2QaPu0aNyiVdG4Rf83ab8vmNk8Ei8KV6xnVLT/tWb2i2jZZKAfcK9oy1zsAAAC30lEQVSZ3Vmh/QH7M7MJ0fZrzez66n4wZtYjOqfCqM7sSpp9l6gniujq6x0zuy8a0+cvZnYogLt/CBzTgrrmkSamgJKDxT3Ad83syDps04vEeF85wPeAbu7el8TwKtcktcsCBpAYeuVeSwyAeBWJntJPA04DrjazrlH7vsDN7t49+WBmdhzwCxLdwfQGTjOzC919KpBPot+5SZXUWb4/MzsVuBI4nUSnnFeb2SnVnOMPgN+4e28gl0T3OhWdBSR3ipwN3BON6fMp8H+S1q2O2os0mAJKDgrR7anfkxissLZWufu2aPya94GyYTnWkAilMo+5+z533wBsAk4CzgUuN7NC4DXgGBL/sQO87u4fVHK804AVUUe2Zb10969Fncn76wc86e6fufsuEp3hfrOabVcCPzGzHwPHu/vuStoc7e7/L2n+A3cvjKYL2P9n8Q8SvcSLNJgCSg4md5G4skke26qE6N9B1HnsIUnrvkia3pc0v4/E0AFlKvYX5oAB17h77+ira9K4U59VUV+1Qw9UI3l/ddqHu88DhgK7gWereMihxMyS/69I/rmUsv/Pol20L5EGU0DJQcPddwCPkQipMpuBU6PpYUCbeuz6YjNrFX0udQKwHngW+KGZtYFE781W86CPrwEDzCw9euBhFPBiHWt5Cbgw6q39MGA48NeqGpvZCcAmd59Bogf0npU0Wx+dV210A2r9hKBIdVrX3ESkRfklMD5p/j7gT2b2OvAcVV/dVGc9iSD5CvADd99jZveTuPW1OroyK6aGYdfdfZuZ3QS8QOJKaIm712mYDHdfbWZzgdejRfe7+xvVbDICuMzM9gL/A0ytpM0zwEBgY3XHjsL4ayQ+LxNpMPVmLiLVssSAh79393NqaDecxPhFP22cyqSl0y0+EamWu28D7qvpRV0Sd2R+WUMbkVrTFZSIiARJV1AiIhIkBZSIiARJASUiIkFSQImISJAUUCIiEqT/Ba3tD7M1DMW6AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"thinkplot.Hist(pmf_sim, label='Simulation')\n",
"thinkplot.Pmf(posterior, color='orange', label='Normalized likelihoods')\n",
"thinkplot.decorate(xlabel='Number of rolls (n)',\n",
" ylabel='PMF')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the posterior distribution based on a uniform prior matches the simulation, it seems like the uniform prior must be correct. But it is not obvious (to me) why."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment