Royal Game of Ur ABC

ur_abc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "from tqdm import trange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import StrMethodFormatter\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set()\n",
    "pct_formatter = StrMethodFormatter('{x:.1%}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use [approximate Bayesian computation](https://en.wikipedia.org/wiki/Approximate_Bayesian_computation) (ABC) to solve the [Game of Ur](https://www.allendowney.com/blog/2018/10/21/the-game-of-ur-problem/) problem that Allen Downey [posed on Twitter](https://twitter.com/AllenDowney/status/1054053993566081026) recently.\n",
    "\n",
    "Our observation is that the piece is on the thirteenth square, and we want to infer the posterior distribution of the number of turns it took to get there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "OBS = 13"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with all Bayesian inference, ABC requires both a prior distribution on the number of rolls and a data generation process (likelihood).  Since it must take at least four rolls to reach spot 13, we place a $U(4, 40)$ prior on the number of rolls.  Note that while the number of rolls needed to reach spot 13 is theoretically unbounded, it would be extremely unlucky not to pass spot 13 within 40 rolls, so we use this truncated uniform prior for computational simplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_rolls = int(np.ceil(OBS / 4))\n",
    "prior = stats.randint(min_rolls, 10 * min_rolls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each turn results in moving a number of spaces drawn from a $\\textrm{Bin}\\left(4, \\frac{1}{2}\\right)$ distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dice = stats.binom(4, 0.5)\n",
    "\n",
    "def simulate_spot(rolls, dice):\n",
    "    return dice.rvs(size=rolls).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest version of ABC generates posterior samples by the following algorithm.\n",
    "\n",
    "1. Generate a sample from the prior.\n",
    "2. Simulate the data generating process given the prior sample.\n",
    "3. If the simulated data matches the observed data, save the sample.  If not, repeat the above process.\n",
    "\n",
    "The following function implements this algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_abc_sample(obs, prior, dice):\n",
    "    while True:\n",
    "        rolls = prior.rvs()\n",
    "        spot = simulate_spot(rolls, dice)\n",
    "        \n",
    "        if spot == obs:\n",
    "            return rolls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now generate 5,000 samples from the posterior using ABC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_SAMPLES = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5000/5000 [01:09<00:00, 72.00it/s] \n"
     ]
    }
   ],
   "source": [
    "post_sample_counts = Counter(generate_abc_sample(OBS, prior, dice) for _ in trange(N_SAMPLES))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "post_df = (pd.DataFrame.from_dict(post_sample_counts, orient='index')\n",
    "             .rename(columns={0: 'samples'})\n",
    "             .sort_index())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>samples</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>716</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1440</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>841</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   samples\n",
       "4       79\n",
       "5      716\n",
       "6     1458\n",
       "7     1440\n",
       "8      841"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "post_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we visualize the posterior distribution of these samples.  The fact that the largest observed sample is much smaller than the point at which we truncated our prior (40) is circumstantial evidence that truncation has not had too large an effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFYCAYAAACyKp7WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XtYVPWi//HPAKKhJEIzoGWhPHnD\nvKWdFM17pjx2tUBump1OHje2tUyFTC1FRd3W9lKaV9IKSjvlPjvFrCx3oW61bWX6VJaGiFwUCQJ1\nxPn94a95IsWZ9pkbrvfreXo2a61Z3/ksdPuZdZm1TDabzSYAAGAIft4OAAAAPIfiBwDAQCh+AAAM\nhOIHAMBAKH4AAAyE4gcAwEACvB3AE0pKKlwyTrNmQSorq3LJWK7ki7nI5BwyOc8Xc5HJOWRynqty\nmc3BdS5jj/8PCAjw93aEK/LFXGRyDpmc54u5yOQcMjnPE7nctsdfXV2tqVOn6tSpUzp37pzGjRun\ndu3aafLkyaqpqZHZbNaCBQsUGBhYa705c+bowIEDMplMSk9PV6dOnZSVlaUtW7aoa9eumjJliiRp\n8+bNKi0t1ZgxY9y1CQAAXHPctsf/8ccfq2PHjtqwYYNeeuklzZs3T4sXL1ZCQoLeeOMN3XLLLdq4\ncWOtdfbs2aNjx44pJydHGRkZysjIkCRt2bJF2dnZOnz4sKqqqnTu3Dlt2rRJSUlJ7ooPAMA1yW3F\nP2zYMD3++OOSpMLCQoWHh2v37t0aOHCgJKl///7Ky8urtU5eXp4GDRokSYqKilJ5ebkqKyvVoEED\nSVJoaKgqKiqUlZWlxMTEy44WAACAq3P7Of74+HhNmjRJ6enpqq6utpd1WFiYSkpKar22tLRUzZo1\ns0+HhoaqpKRENptNVqtVxcXF8vPz0/79+xUUFKS0tDStW7fO3ZsAAMA1w+1X9WdnZ+vQoUN65pln\n9NvnATnzbKBfXzNy5EilpKQoNjZWK1asUGpqqhYtWqRVq1YpLS1NJ0+eVERERJ3jNGsW5LILJq52\npaQ3+WIuMjmHTM7zxVxkcg6ZnOfuXG4r/q+//lphYWFq3ry52rdvr5qaGjVu3Fhnz55Vo0aNVFRU\nJIvFUmsdi8Wi0tJS+3RxcbHMZrNiY2MVGxuro0eP6vDhw+rYsaOsVqv8/PwUERGhgoKCqxa/q76y\nYTYHu+yrga7ki7nI5BwyOc8Xc5HJOWRynqtyeeXrfHv37tWaNWskXTqEX1VVpV69eik3N1eStG3b\nNvXp06fWOjExMfblBw8elMViUZMmTezLly5dqvHjx0uSrFarbDabCgsLL/sAAQAArsxtxR8fH6/T\np08rISFB//Vf/6Xp06dr/Pjxevfdd5WQkKAzZ87o/vvvlyRNnDhRZ8+eVbdu3RQdHa34+HjNnj1b\nM2bMsI+3d+9eRUZGKjw8XJI0fPhwxcfHy9/fXy1btnTXZgAAcE0x2Zw52V7PuepwzrV+aMiVyOQc\nMjnPF3ORyTlkcl69PtQPAAB8D8UPAICBUPwAABiIIZ7OB9caM+8jl4yzZuoAl4wDAHAee/wAABgI\nxQ8AgIFQ/AAAGAjn+HFN4LoDAHAOe/wAABgIxQ8AgIFQ/AAAGAjFDwCAgVD8AAAYCMUPAICBUPwA\nABgIxQ8AgIFQ/AAAGAjFDwCAgVD8AAAYCMUPAICBUPwAABgIxQ8AgIFQ/AAAGAjFDwCAgVD8AAAY\nCMUPAICBUPwAABgIxQ8AgIFQ/AAAGAjFDwCAgVD8AAAYCMUPAICBUPwAABgIxQ8AgIFQ/AAAGAjF\nDwCAgVD8AAAYCMUPAICBUPwAABgIxQ8AgIEEuHPw+fPna9++fbpw4YKeeOIJffTRRzp48KBCQkIk\nSY899pj69etXa505c+bowIEDMplMSk9PV6dOnZSVlaUtW7aoa9eumjJliiRp8+bNKi0t1ZgxY9y5\nCQAAXFPcVvy7du3Sd999p5ycHJWVlemBBx7QnXfeqaeeekr9+/e/4jp79uzRsWPHlJOToyNHjig9\nPV05OTnasmWLsrOz9eijj6qqqkr+/v7atGmTVq5c6a74AABck9xW/D169FCnTp0kSddff72qq6tV\nU1Nz1XXy8vI0aNAgSVJUVJTKy8tVWVmpBg0aSJJCQ0NVUVGh9957T4mJiQoMDHRXfAAArkluO8fv\n7++voKAgSdLGjRt11113yd/fXxs2bFBKSoomTpyo06dP11qntLRUzZo1s0+HhoaqpKRENptNVqtV\nxcXF8vPz0/79+xUUFKS0tDStW7fOXZsAAMA1x63n+CVp+/bt2rhxo9asWaOvv/5aISEhat++vV59\n9VUtXbpU06dPr3Ndm80mSRo5cqRSUlIUGxurFStWKDU1VYsWLdKqVauUlpamkydPKiIios5xmjUL\nUkCAv0u2x2wOdsk4ruarua7GFzN7IxO/B+f5Yi4yOYdMznN3LrcW/86dO7V8+XKtWrVKwcHB6tmz\np33ZgAEDNHPmzFqvt1gsKi0ttU8XFxfLbDYrNjZWsbGxOnr0qA4fPqyOHTvKarXKz89PERERKigo\nuGrxl5VVuWR7zOZglZRUuGQsV/LVXI74YmZPZ/LFPztfzCT5Zi4yOYdMznNVrqt9eHDbof6KigrN\nnz9fK1assF/FP378eOXn50uSdu/erVtvvbXWOjExMcrNzZUkHTx4UBaLRU2aNLEvX7p0qcaPHy9J\nslqtstlsKiwslMVicddmAABwTXHbHv/777+vsrIyTZgwwT7vwQcf1IQJE3TdddcpKChIc+fOlSRN\nnDhRc+fOVbdu3RQdHa34+HiZTCbNmDHDvu7evXsVGRmp8PBwSdLw4cMVHx+v1q1bq2XLlu7aDAAA\nriluK/64uDjFxcVdNv+BBx64bN6LL75o/3nSpElXHK979+7q3r27fToxMVGJiYkuSAoAgHFw5z4A\nAAyE4gcAwEAofgAADITiBwDAQCh+AAAMhOIHAMBAKH4AAAyE4gcAwEAofgAADITiBwDAQCh+AAAM\nhOIHAMBAKH4AAAyE4gcAwEAofgAADITiBwDAQCh+AAAMhOIHAMBAKH4AAAyE4gcAwEAofgAADITi\nBwDAQCh+AAAMhOIHAMBAKH4AAAyE4gcAwEAofgAADITiBwDAQCh+AAAMhOIHAMBAKH4AAAyE4gcA\nwEAofgAADITiBwDAQCh+AAAMhOIHAMBAKH4AAAyE4gcAwEAofgAADITiBwDAQALcOfj8+fO1b98+\nXbhwQU888YRuu+02TZ48WTU1NTKbzVqwYIECAwNrrTNnzhwdOHBAJpNJ6enp6tSpk7KysrRlyxZ1\n7dpVU6ZMkSRt3rxZpaWlGjNmjDs3AQCAa4rb9vh37dql7777Tjk5OVq1apXmzJmjxYsXKyEhQW+8\n8YZuueUWbdy4sdY6e/bs0bFjx5STk6OMjAxlZGRIkrZs2aLs7GwdPnxYVVVVOnfunDZt2qSkpCR3\nxQcA4JrksPgnTpyozz///A8P3KNHD/31r3+VJF1//fWqrq7W7t27NXDgQElS//79lZeXV2udvLw8\nDRo0SJIUFRWl8vJyVVZWqkGDBpKk0NBQVVRUKCsrS4mJiZcdLQAAAFfnsPgHDx6s7OxsxcbGatmy\nZTp58qRTA/v7+ysoKEiStHHjRt11112qrq62l3VYWJhKSkpqrVNaWqpmzZrZp0NDQ1VSUiKbzSar\n1ari4mL5+flp//79CgoKUlpamtatW+fstgIAYHgOz/EPGzZMw4YNU1VVlT7++GM99dRTaty4sR59\n9FH16tXL4Rts375dGzdu1Jo1a3T33Xfb59tsNofr/vqakSNHKiUlRbGxsVqxYoVSU1O1aNEirVq1\nSmlpaTp58qQiIiLqHKdZsyAFBPg7fD9nmM3BLhnH1Xw119X4YmZvZOL34DxfzEUm55DJee7O5dTF\nfdXV1dq2bZveffddXbx4Uf3799drr72m3bt3a+LEiXWut3PnTi1fvlyrVq1ScHCwgoKCdPbsWTVq\n1EhFRUWyWCy1Xm+xWFRaWmqfLi4ultlsVmxsrGJjY3X06FEdPnxYHTt2lNVqlZ+fnyIiIlRQUHDV\n4i8rq3JmMx0ym4NVUlLhkrFcyVdzOeKLmT2dyRf/7Hwxk+SbucjkHDI5z1W5rvbhweGh/rS0NA0Z\nMkQHDhzQ1KlTlZ2drYSEBL3yyivauXNnnetVVFRo/vz5WrFihUJCQiRJvXr1Um5uriRp27Zt6tOn\nT611YmJi7MsPHjwoi8WiJk2a2JcvXbpU48ePlyRZrVbZbDYVFhZe9gECAABcmcM9/vbt2+u5556z\nn6+XpH/961/q0qWLZs+eXed677//vsrKyjRhwgT7vHnz5mnatGnKyclRixYtdP/990u6dAHh3Llz\n1a1bN0VHRys+Pl4mk0kzZsywr7t3715FRkYqPDxckjR8+HDFx8erdevWatmy5R/fcgAADKjO4v/5\n55915swZ/f3vf1f//v116tQpSZf2tKdMmaLc3Fx16NChzoHj4uIUFxd32fy1a9deNu/FF1+0/zxp\n0qQrjte9e3d1797dPp2YmKjExMQ63x8AAFyuzuL/4osvlJWVpUOHDmnUqFH2+X5+furdu7dHwgEA\nANeqs/j79u2rvn376s0339TIkSM9mQkAALhJncW/adMmPfTQQyoqKrLfiOe3/vznP7s1GAAAcL06\ni9/P79IF/wEBbr2dPwAA8KA6W/2+++7TxYsXNW7cOE/mAQAAblRn8Xfo0EEmk+my+TabTSaTSYcO\nHXJrMAAA4Hp1Fv/hw4c9mQMAAHiAw4v7rnRhn8TFfQAA1EcOL+7z93fNw20AAID31Vn8DzzwgCQp\nNTVVP//8s44ePSqTyaRWrVrVun8+AACoPxx+V2/dunV65ZVXFBkZqYsXL+r48eNKTU3ldrkAANRD\nDot/06ZN2r59u4KDLz3ir7y8XElJSRQ/AAD1kMPH8oaHh9tLX5KaNm2qm2++2a2hAACAe9S5x79x\n40ZJUosWLTR27Fj16tVLfn5+2rVrl/3RuAAAoH6ps/j37dtn/7lZs2b2G/YEBwerurra/ckAAIDL\n1Vn8c+fOrXOl1157zS1hAACAezm8uO/QoUNavny5ysrKJEnnz5/XyZMnlZKS4vZwAADAtRxe3Pf8\n88/r7rvvVnl5ucaMGaPIyEjNnz/fE9kAAICLOSz+Ro0aKTY2VsHBwerXr58yMjK0evVqT2QDAAAu\n5rD4z507p2+//VYNGzbUnj17VF5eroKCAk9kAwAALubwHP+kSZP0008/6cknn9TkyZN16tQpPf74\n457IBgAAXMxh8d9+++32n3Nzc90aBgAAuJfDQ/3//Oc/9dBDD6lLly7q2rWr4uLitH//fk9kAwAA\nLuZwj/+FF15Qenq6unXrJpvNpn379mnmzJnavHmzJ/IBAAAXclj8YWFh6tmzp306JiZGLVq0cGso\nAADgHnUWf35+viTptttu05o1a+z36s/Ly1OHDh08FhAAALhOncU/atQomUwm2Ww2SdKGDRvsy0wm\nk5588kn3pwPqsTHzPnLZWGumDnDZWACMrc7i/+gj1/2jBQAAfIPDc/zFxcV66aWX9NVXX8lkMqlL\nly6aMGGCQkNDPZEPAAC4kMOv802fPl3R0dFatGiRFi5cqNatWys9Pd0T2QAAgIs53OOvrq5WYmKi\nfbpNmzacBgAAoJ5yuMdfXV2t4uJi+/TJkyd1/vx5t4YCAADu4XCPf9y4cXrwwQdlNptls9l0+vRp\nZWRkeCIbAABwMYfF37dvX23fvl1Hjx6VJLVq1UoNGzZ0dy4AAOAGDg/1p6SkqFGjRmrXrp3atWtH\n6QMAUI853ONv3769/vrXv6pr165q0KCBff5vb+MLAADqB4fFf+jQIUnS3r177fNMJhPFDwBAPeSw\n+NevX++JHAAAwAMcnuPfs2ePHnzwQXXu3FldunRRXFyc/vWvf3kiGwAAcDGHe/xz5szRlClTdPvt\nt8tms2nv3r2aOXOm3n33XU/kAwAALuRwjz8kJEQ9e/ZUYGCgGjZsqJiYGIWHhzs1+LfffqtBgwbZ\nn+w3depUDR8+XMnJyUpOTtaOHTsuW2fOnDmKi4tTfHy8vvzyS0lSVlaW4uPjlZmZaX/d5s2btWbN\nGqdyAACASxzu8Xfu3Fnr1q1T7969dfHiRe3atUtRUVHKz8+XJLVs2fKK61VVVWnWrFmXXQT41FNP\nqX///ldcZ8+ePTp27JhycnJ05MgRpaenKycnR1u2bFF2drYeffRRVVVVyd/fX5s2bdLKlSv/6PYC\nAGBoDov/b3/7myTptddeqzV/69atMplM+vDDD6+4XmBgoFauXPmHyjkvL0+DBg2SJEVFRam8vFyV\nlZX2rxGGhoaqoqJC7733nhITExUYGOj02AAAwIni/3cfyBMQEKCAgMuH37Bhg9auXauwsDA999xz\ntR7vW1paqujoaPt0aGioSkpKZLPZZLVaVVxcLD8/P+3fv18dOnRQWlqa2rZtq9GjR181S7NmQQoI\n8P+3tuP3zOZgl4zjar6a62p8MbMvZpI8n4vfg/PI5BwyOc/duRwWvyvdd999CgkJUfv27fXqq69q\n6dKlmj59ep2vt9lskqSRI0cqJSVFsbGxWrFihVJTU7Vo0SKtWrVKaWlpOnnypCIiIuocp6ysyiX5\nzeZglZRUuGQsV/LVXI74YmZfzCR5Npev/n3yxVxkcg6ZnOeqXFf78ODw4j5X6tmzp9q3by9JGjBg\ngL799ttayy0Wi0pLS+3TxcXFMpvNio2N1ZtvvqnevXvr7Nmz6tixo6xWq/z8/BQREaGCggJPbgYA\nAPWWw+L/7SN5/6/Gjx9vvyhw9+7duvXWW2stj4mJUW5uriTp4MGDslgsatKkiX350qVLNX78eEmS\n1WqVzWZTYWGhLBaLyzICAHAtc3iof9KkSZdd2OeMr7/+WpmZmSooKFBAQIByc3OVlJSkCRMm6Lrr\nrlNQUJDmzp0rSZo4caLmzp2rbt26KTo6WvHx8TKZTJoxY4Z9vL179yoyMtL+VcLhw4crPj5erVu3\nrvObBQAAoDaHxR8ZGanJkydf9pCeESNGXHW9jh07XvF2v0OGDLls3osvvmj/edKkSVccr3v37ure\nvbt9OjExUYmJiY7iAwCA33BY/FarVf7+/vab6fzKUfEDAADf47D4fz0cf+bMGZlMJjVt2tTtoXDJ\nmHn/3lcpr2TN1AEuGwsAUH85LP79+/dr8uTJ+uWXX2Sz2RQSEqIFCxbotttu80Q+AADgQg6L/y9/\n+YtefvlltWnTRpL0zTffKCMjQ6+//rrbwwEAANdy+HU+Pz8/e+lLUocOHeTv75q74AEAAM9yqvi3\nbdumyspKVVZW6v3336f4AQCopxwe6n/++ec1a9YsPfvsszKZTOrSpYuef/55T2QDAAAu5tT3+Fev\nXu2JLAAAwM3qLP7Zs2dr2rRpSkhIkMlkumw5F/cBAFD/1Fn8v96gZ8KECR4LAwAA3KvO4m/Xrp0k\n6YMPPtCzzz7rsUAAAMB9HF7V7+/vr7y8PJ07d04XL160/wcAAOofhxf3vf3228rKypLNZpPJZLL/\n76FDhzyRDwAAuJDD4t+3b58ncgAAAA9weKi/vLxcmZmZeuaZZyRJH330kU6fPu32YAAAwPUcFv+0\nadPUvHlz5efnS5LOnz+vKVOmuD0YAABwPYfFf/r0aaWkpKhBgwaSpHvuuUdnz551ezAAAOB6Dotf\nkqxWq/0mPqWlpaqqqnJrKAAA4B4OL+5LTEzUiBEjVFJSorFjx+qrr77ie/0AANRTDot/2LBh6tat\nm7744gsFBgbqhRde0PXXX++JbAAAwMUcHup/7LHHFBERoaFDh2rgwIGyWCxKTEz0RDYAAOBide7x\nb968WcuWLdOJEyfUr18/+3yr1aobbrjBE9kAAICL1Vn89957r2JjY/Xss89q/Pjx9vl+fn6yWCwe\nCQcAAFzrqof6/f39NXXqVFVVVenGG2/UDz/8oHfeeYcb+AAAUE85PMc/efJkFRcX6+jRo5o3b55C\nQkK4qh8AgHrKYfFXV1crJiZGW7duVVJSkhITE2W1Wj2RDQAAuJhTxX/69Gnl5uaqX79+stlsKi8v\n90Q2AADgYg6Lf/jw4br77rt15513qnnz5lq2bJn+4z/+wxPZAACAizm8gc+oUaM0atSoWtPBwcFu\nDQUAANzD4R7/kSNHlJKSom7duun222/XhAkTdOzYMU9kAwAALuaw+GfNmqUxY8boH//4hz799FPF\nx8dr5syZHogGAABczWHx22w29evXT0FBQWrcuLEGDx6smpoaT2QDAAAu5rD4rVarDh48aJ/+8ssv\nKX4AAOophxf3TZkyRU8//bT9bn1ms1mZmZluDwYAAFzPYfF37txZW7duVUVFhUwmk5o0aeKJXAAA\nwA3qLP7Kykq9/PLL+uGHH9SjRw+NGjVKAQEOPycAAAAfVuc5/l+v3I+Li9P333+vpUuXeioTAABw\nkzp34QsKCrRw4UJJ0l133aXRo0d7KhMAAHCTOvf4f3tY39/f3yNhAACAe9VZ/CaT6arTAACg/qnz\nUP8XX3yhfv362adPnTplfzqfyWTSjh07HA7+7bffaty4cRo9erSSkpJUWFioyZMnq6amRmazWQsW\nLFBgYGCtdebMmaMDBw7IZDIpPT1dnTp1UlZWlrZs2aKuXbtqypQpkqTNmzertLRUY8aM+fe2HAAA\nA6qz+Ldu3fp/GriqqkqzZs1Sz5497fMWL16shIQEDR06VIsWLdLGjRuVkJBgX75nzx4dO3ZMOTk5\nOnLkiNLT05WTk6MtW7YoOztbjz76qKqqquTv769NmzZp5cqV/6eMAAAYTZ2H+m+88car/udIYGCg\nVq5cKYvFYp+3e/duDRw4UJLUv39/5eXl1VonLy9PgwYNkiRFRUWpvLxclZWVatCggSQpNDRUFRUV\nysrKUmJi4mVHCwAAwNW57Yv5AQEBl33vv7q62l7WYWFhKikpqbW8tLRU0dHR9unQ0FCVlJTIZrPJ\narWquLhYfn5+2r9/vzp06KC0tDS1bdvW4TcOmjULUkCAay5QNJvr5yOJfTE3mZzn6Vz8HpxHJueQ\nyXnuzuW1O/LYbDanXzNy5EilpKQoNjZWK1asUGpqqhYtWqRVq1YpLS1NJ0+eVERERJ3jlJVVuSSz\n2RyskpIKl4zlab6Ym0zO82QuX/177ou5yOQcMjnPVbmu9uHB4UN6XCkoKEhnz56VJBUVFdU6DSBJ\nFotFpaWl9uni4mKZzWbFxsbqzTffVO/evXX27Fl17NhRVqtVfn5+ioiIUEFBgSc3AwCAesujxd+r\nVy/l5uZKkrZt26Y+ffrUWh4TE2NffvDgQVksllrPBli6dKnGjx8v6dJTA202mwoLCy/7AAEAAK7M\nbYf6v/76a2VmZqqgoEABAQHKzc3VwoULNXXqVOXk5KhFixa6//77JUkTJ07U3Llz1a1bN0VHRys+\nPl4mk0kzZsywj7d3715FRkYqPDxckjR8+HDFx8erdevWatmypbs2AwCAa4rbir9jx45av379ZfPX\nrl172bwXX3zR/vOkSZOuOF737t3VvXt3+3RiYqISExNdkBQAAOPw6KF+AADgXRQ/AAAGQvEDAGAg\nFD8AAAZC8QMAYCAUPwAABkLxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAABkLxAwBgIBQ/\nAAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAABkLxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAA\nBkLxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAABhLg7QAAPGfMvI9cNtaaqQNcNhYAz2GP\nHwAAA6H4AQAwEIofAAADofgBADAQih8AAAOh+AEAMBCKHwAAA6H4AQAwEI8W/+7du3XnnXcqOTlZ\nycnJmjVrVq3ln3/+uUaMGKG4uDgtW7ZMkvTjjz8qPj5eycnJKisrkyRVVFRo9OjRunjxoifjAwBQ\n73n8zn133HGHFi9efMVls2fP1urVqxUeHq6kpCQNGTJE77zzjp555hnl5+dr69atGjlypFasWKEn\nnnhCfn4csAAA4I/wmebMz89X06ZN1bx5c/n5+alv377Ky8vTzz//LLPZLLPZrPLychUUFCg/P189\ne/b0dmQAAOodjxf/999/r7Fjx2rkyJH67LPP7PNLSkoUGhpqnw4NDVVJSYkiIiL0008/6ejRo7rx\nxhu1ZMkSjR49WtOnT9f06dN15swZT28CAAD1lkcP9UdGRio1NVVDhw5Vfn6+UlJStG3bNgUGBta5\nzsMPP6z09HQ1btxYKSkpCg4O1u7duzV06FBJUnZ2tsaOHXvV923WLEgBAf4u2QazOdgl43iaL+Ym\nk/N8MZc3MvF7cA6ZnOOLmST35/Jo8YeHh2vYsGGSpJtvvlk33HCDioqK1LJlS1ksFpWWltpfW1RU\nJIvFovDwcK1evVqS9Kc//UkZGRlauHChhg0bposXL+rvf/+7w/ctK6tySX6zOVglJRUuGcvTfDE3\nmZzni7k8nckX//9HJueQyXmuynW1Dw8ePdS/efNme4mXlJTo1KlTCg8PlyTddNNNqqys1PHjx3Xh\nwgV9/PHHiomJsa+7fft29ejRQyEhIQoLC9OJEydUWFgoi8XiyU0AAKBe8+ge/4ABAzRp0iR9+OGH\nslqtmjlzpv73f/9XwcHBGjx4sGbOnKmnn35akjRs2DC1atVKknThwgVt3LhRS5YskSQ9+OCDmjJl\niiRpwYIFntwEAADqNY8Wf5MmTbR8+fI6l/fo0UM5OTmXzQ8ICKi13i233KLs7Gy3ZAQA4FrmM1/n\nAwAA7kfxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAABkLxAwBgIBQ/AAAGQvEDAGAgFD8A\nAAZC8QMAYCAUPwAABkLxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMAYCAUPwAABkLxAwBgIBQ/AAAG\nQvEDAGAgAd4OAMDYxsz7yGVjrZk6wGVjAdcq9vgBADAQih8AAAOh+AEAMBCKHwAAA6H4AQAwEK7q\n//+4shgAYATs8QMAYCAUPwBSF6j0AAANiklEQVQABkLxAwBgIBQ/AAAGQvEDAGAgFD8AAAZC8QMA\nYCAUPwAABkLxAwBgIBQ/AAAGwi17AeAKXHUbb27hDV/j8T3+OXPmKC4uTvHx8fryyy9rLfv88881\nYsQIxcXFadmyZZKkH3/8UfHx8UpOTlZZWZkkqaKiQqNHj9bFixc9HR8AgHrNo8W/Z88eHTt2TDk5\nOcrIyFBGRkat5bNnz9aSJUv05ptv6rPPPtP333+vt99+W88884weeughbd26VZK0YsUKPfHEE/Lz\n40wFAAB/hEebMy8vT4MGDZIkRUVFqby8XJWVlZKk/Px8NW3aVM2bN5efn5/69u2rvLw8/fzzzzKb\nzTKbzSovL1dBQYHy8/PVs2dPT0YHAOCa4NFz/KWlpYqOjrZPh4aGqqSkRE2aNFFJSYlCQ0NrLcvP\nz1dERIR++uknHTt2TDfeeKOWLFmi0aNHa/r06ZKkp556SiEhIZ7cDADwCl+87sAXH2nui5l8iclm\ns9k89WbPPfec+vbta9/rHzlypObMmaNWrVpp//79Wr16tf3c/ttvv638/HwlJiYqPT1djRs3VkpK\ninJzcxUWFqbOnTtLkg4cOKCxY8d6ahMAAKjXPHqo32KxqLS01D5dXFwss9l8xWVFRUWyWCwKDw/X\n6tWrtXjxYq1du1Z/+tOfdPz4cd14441q3ry5jh8/7slNAACgXvNo8cfExCg3N1eSdPDgQVksFjVp\n0kSSdNNNN6myslLHjx/XhQsX9PHHHysmJsa+7vbt29WjRw+FhIQoLCxMJ06cUGFhoSwWiyc3AQCA\nes2jh/olaeHChdq7d69MJpNmzJihb775RsHBwRo8eLD++c9/auHChZKku+++W4899pgk6cKFC0pN\nTdWSJUvUoEEDHTt2TFOmTJEkLViwQC1btvTkJgAAUG95vPgBAID38EV4AAAMhOIHAMBAKP4/4OzZ\nsxo0aJDeeecdb0fR7t27deeddyo5OVnJycmaNWuWtyPZbd68Wffee68efPBB7dixw9tx9Pbbb9t/\nT8nJyeratau3I+mXX35RamqqkpOTFR8fr507d3o7ki5evKjnnnvOfovsI0eOeDXPt99+q0GDBmnD\nhg2SpMLCQiUnJyshIUF//vOfdf78ea9nkqTXXntN0dHR+uWXXzye50qZCgsLNXr0aCUlJWn06NEq\nKSnxeqYvvvhCI0eOVHJysh577DGdPn3a45mulOtXO3fuVNu2bX0i09SpUzV8+HD7v1fu+DeUh/T8\nAa+88oqaNm3q7Rh2d9xxhxYvXuztGLWUlZVp2bJl2rRpk6qqqrRkyRL169fPq5kefvhhPfzww5Iu\n3TZ6y5YtXs0jSf/zP/+jVq1a6emnn1ZRUZFGjRplvyW1t3z44YeqqKhQdna2fvrpJ2VkZGjFihVe\nyVJVVaVZs2bVukPn4sWLlZCQoKFDh2rRokXauHGjEhISvJrp3Xff1alTp7z27aIrZXrppZf0yCOP\naNiwYXr99de1du1aTZ482auZ1q5dq/nz56tly5ZaunSp3nrrLY/ff+VKuSTp3LlzevXVV+1fLfeF\nTE899ZT69+/vtvdlj99JR44c0ffff+/1EvN1eXl56tmzp5o0aSKLxeJTRyIkadmyZRo3bpy3Y6hZ\ns2Y6c+aMJOnnn39Ws2bNvJxIOnr0qDp16iRJuvnmm3XixAnV1NR4JUtgYKBWrlxZq1B3796tgQMH\nSpL69++vvLw8r2caNGiQJk6cKJPJ5NEsV8s0Y8YMDRkyRFLtv2fezLR48WK1bNlSNptNRUVFioiI\n8GimunJJ0vLly5WQkKDAwECfyeRuFL+TMjMzNXXqVG/HqOX777/X2LFjNXLkSH322WfejiNJOn78\nuM6ePauxY8cqISHB4/84X82XX36p5s2be+WT/e/FxsbqxIkTGjx4sJKSkuxfT/WmNm3a6B//+Idq\namr0ww8/KD8/3/5ETE8LCAhQo0aNas2rrq62/+McFhbm8UPYV8r0631IvOVKmYKCguTv76+amhq9\n8cYbGj58uNczSdKnn36qe+65R6Wlpbr33ns9mqmuXD/++KMOHz6soUOHejxPXZkkacOGDUpJSdHE\niRPdclqE4nfCu+++qy5duvjU/QIiIyOVmpqqV155RZmZmXr22We9cs7zSs6cOaOlS5dq3rx5SktL\nk698Y3Tjxo164IEHvB1DkvTee++pRYsW+uCDD5SVlaUXXnjB25HUt29f3XbbbUpMTFRWVpZat27t\nM392v+eruXxFTU2NJk+erDvvvNNnHmh21113aevWrWrdurVeffVVb8eRJM2dO1dpaWnejlHLfffd\np0mTJum1115T+/bttXTpUpe/B+f4nbBjxw7l5+drx44dOnnypAIDAxUREaFevXp5LVN4eLiGDRsm\n6dJh2RtuuEFFRUVe/3ASFhamrl27KiAgQDfffLMaN26s06dPKywszKu5pEuHiqdNm+btGJKk/fv3\nq3fv3pKkdu3aqbi4WDU1NfL39/dqrokTJ9p/HjRokE/8uf0qKChIZ8+eVaNGjey39MaVpaWl6ZZb\nblFqaqq3o0iSPvjgAw0ePFgmk0lDhgzRkiVLvB1JRUVF+uGHHzRp0iRJl24hn5SUdNmFf5722w9q\nAwYM0MyZM13+HuzxO+Gll17Spk2b9NZbb+nhhx/WuHHjvFr60qUr51evXi1JKikp0alTpxQeHu7V\nTJLUu3dv7dq1SxcvXlRZWZmqqqp84vx1UVGRGjdu7JXzeFdyyy236MCBA5KkgoICNW7c2Oulf/jw\nYfvez6effqoOHTrIz893/ono1auX/Zbf27ZtU58+fbycyDdt3rxZDRo00JNPPuntKHZLlizRoUOH\nJF16sFqrVq28nOjSztP27dv11ltv6a233pLFYvF66UvS+PHjlZ+fL+nSzsqtt97q8vdgj7+eGjBg\ngCZNmqQPP/xQVqtVM2fO9IlSCw8P15AhQ/TII49IkqZNm+YT5fH7xz57W1xcnNLT05WUlKQLFy64\n5VP9H9WmTRvZbDaNGDFCDRs2tN8+2xu+/vprZWZmqqCgQAEBAcrNzdXChQs1depU5eTkqEWLFrr/\n/vu9nqlXr176/PPPVVJSoscff1xdunTx6BX0V8p06tQpNWzYUMnJyZKkqKgoj/79ulKm2bNn6/nn\nn5e/v78aNWqk+fPneyzP1XItWbLEq491v1KmpKQkTZgwQdddd52CgoI0d+5cl78vt+wFAMBAvL8r\nBgAAPIbiBwDAQCh+AAAMhOIHAMBAKH4AAAyE4geuEcePH1fbtm21efPmWvMHDBjgkvHbtm2rCxcu\nuGSsuuTm5mrgwIF6++23/631BwwYoGPHjumdd96x35gFQG0UP3ANiYyM1LJly1RZWentKP+WTz75\nRI899pj9aYoAXI8b+ADXEIvFot69e+vll1++7EYy77zzjj7//HP7jXmSk5P13//93/L399fy5csV\nERGhr776Sp07d1bbtm31wQcf6MyZM1q5cqX9aWrLly/Xrl279MsvvygzM1Nt2rTR4cOHlZmZqQsX\nLshqtWr69Onq0KGDkpOT1a5dOx06dEhZWVm17ky4Y8cOLVu2TI0aNdJ1112nWbNm6YsvvtAnn3yi\nffv2yd/fX3FxcfbX/36snTt3XrZ+XXeuXLhwoXbt2qXAwECFh4crMzPTJ252BXgLe/zANebRRx/V\nJ598oh9++MHpdb788ktNmTJFmzZt0t/+9jddf/31Wr9+vaKjo7V161b766KiorRhwwYlJCTYHx7y\nzDPP6Pnnn9f69es1c+bMWs9DCAoK0oYNG2qVfnV1taZNm6YlS5Zo/fr1uuuuu/TSSy/pnnvuUZ8+\nffSf//mftUr/92OdP3/+iutfSXl5uV5//XXl5OTojTfe0ODBg1VaWur07wW4FlH8wDUmMDBQkydP\nVkZGhtPrREVFKSQkRA0bNlRISIi6du0q6dItmH972iAmJkaS1K1bN3333Xc6deqUfvzxRz377LNK\nTk5WRkaGKisrdfHiRfvrfu/o0aMKCwuzH0W444479NVXXznM+OtYf2T9pk2bqk+fPkpKStKaNWvU\nrVs3tWjRwtlfC3BN4lA/cA3q27ev3nzzTX3wwQf2eSaTqdZrrFar/effPyDot9O/vav3r89dsNls\nMplMCgwMVIMGDbR+/for5mjQoMFl836f49exHPl1rD+6/uLFi3XkyBF98sknSkpK0pIlS9S+fXuH\n7wdcq9jjB65R6enp+stf/qLz589Lkpo0aaKTJ09Kkk6dOqXvvvvuD4+Zl5cn6dJjhdu0aaPg4GDd\ndNNN+uSTTyRJP/74o8Pnh0dGRurUqVM6ceKEfczOnTs7neGPrJ+fn69169YpKipKY8aM0eDBg3X4\n8GGn3wu4FrHHD1yjbr75Zg0ZMkTLly+XdOkw/erVq/XII48oKirKfjjfWf7+/vruu++UnZ2tsrIy\nLViwQJKUmZmp2bNn69VXX9WFCxc0derUq47TqFEjZWRkaOLEiQoMDFRQUNAfOi3xR9YPDw/XN998\noxEjRqhx48Zq2rSpzzyjHvAWns4HAICBcKgfAAADofgBADAQih8AAAOh+AEAMBCKHwAAA6H4AQAw\nEIofAAADofgBADCQ/wciAKU2ZyPjdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x396 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = (post_df.div(N_SAMPLES)\n",
    "             .plot(kind='bar', legend=False, rot=0))\n",
    "\n",
    "ax.set_xlabel(\"Number of rolls\");\n",
    "\n",
    "ax.yaxis.set_major_formatter(pct_formatter);\n",
    "ax.set_ylabel(\"Posterior probability\");"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}