From 5be9d44529aeb6636d2632e2da9ec9e1f601e106 Mon Sep 17 00:00:00 2001
From: Allen Downey <downey@allendowney.com>
Date: Sun, 14 Jan 2024 21:29:59 -0500
Subject: [PATCH] Adding beta-binomial.ipynb

---
 examples/beta-binomial.ipynb | 902 +++++++++++++++++++++++++++++++++++
 1 file changed, 902 insertions(+)
 create mode 100644 examples/beta-binomial.ipynb

diff --git a/examples/beta-binomial.ipynb b/examples/beta-binomial.ipynb
new file mode 100644
index 00000000..c90cce95
--- /dev/null
+++ b/examples/beta-binomial.ipynb
@@ -0,0 +1,902 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "f07c7715",
+   "metadata": {},
+   "source": [
+    "# The All-Knowing Cube of Probability"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "94eb4081",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import empiricaldist\n",
+    "except ImportError:\n",
+    "    !pip install empiricaldist"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "e68dfce7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get utils.py\n",
+    "\n",
+    "from os.path import basename, exists\n",
+    "\n",
+    "def download(url):\n",
+    "    filename = basename(url)\n",
+    "    if not exists(filename):\n",
+    "        from urllib.request import urlretrieve\n",
+    "        local, _ = urlretrieve(url, filename)\n",
+    "        print('Downloaded ' + local)\n",
+    "    \n",
+    "download('https://github.com/AllenDowney/ThinkBayes2/raw/master/soln/utils.py')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "de8f6c49",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "949ce91c",
+   "metadata": {},
+   "source": [
+    "Suppose you run $n$ trials where the probability of success is $p$.\n",
+    "To compute the probability of $k$ successes, we can use the binomial distribution.\n",
+    "\n",
+    "For example, here's a range of value for $k$ and $n$, and a discrete grid of values for $p$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "78f531f8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ks = np.arange(101)\n",
+    "ns = np.arange(101)\n",
+    "ps = np.linspace(0, 1, 101)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d77a2b77",
+   "metadata": {},
+   "source": [
+    "We can use `meshgrid` to make a 3-D grid of $k$, $n$, and $p$, and `binom` to evaluate the binomial PMF at each point."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "9e450ec2",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(101, 101, 101)"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "\n",
+    "from scipy.stats import binom\n",
+    "\n",
+    "K, N, P = np.meshgrid(ks, ns, ps, indexing='ij')\n",
+    "cube = binom.pmf(K, N, P)\n",
+    "cube.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "95297582",
+   "metadata": {},
+   "source": [
+    "The result is the all-knowing cube of probability, so-called because it can answer all of our questions about Bernoulli trials.\n",
+    "Allow me to demonstrate."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d482aa8d",
+   "metadata": {},
+   "source": [
+    "## The binomial distribution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b7377f5f",
+   "metadata": {},
+   "source": [
+    "Suppose we are given $n$ and $p$, and we would like to know the distribution of $k$.\n",
+    "We can answer that question by selecting a vector from the cube along the $k$ axis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "ae639fe4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n = 50\n",
+    "p = 50\n",
+    "pmf_k = cube[:, n, p]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ab6c3676",
+   "metadata": {},
+   "source": [
+    "The result is a normalized PMF."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "07e338a7",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.9999999999999996"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pmf_k.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "da690210",
+   "metadata": {},
+   "source": [
+    "Here's what it looks like."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "6daa1886",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAl3UlEQVR4nO3dfVBU5+G38S8vsqsmYCKVFd+wqQ1aCagIP4xT7GSnmNA2NC0xjomUOGZspdVsh0SswnRMirHqYCINtY3JZCrV2kar0ZKh+NKmEomATW1Sk2mTQDUL2kRQjGDZ8/zRJ9tsXQ0Y4AD39Zk5U/fsfXbvc2eq15x9C7EsyxIAAIBBQu2eAAAAQF8jgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgnHC7J9Af+Xw+nT59WjfeeKNCQkLsng4AAOgCy7J0/vx5xcbGKjT02td4CKAgTp8+rXHjxtk9DQAAcB0aGxs1duzYa44hgIK48cYbJf1nASMjI22eDQAA6IrW1laNGzfO/+/4tRBAQXz0sldkZCQBBADAANOVt6/wJmgAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYJt3sCwEASt2JfwO131mbaNBMAwKfBFSAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBx+DFU4Br+98dPP2kMP44KAAMDV4AAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcWwPoNLSUsXFxcnpdCo1NVU1NTVXHfvXv/5V3/jGNxQXF6eQkBCVlJR86scEAADmsTWAduzYIY/Ho6KiItXV1SkxMVEZGRlqbm4OOv7ixYv67Gc/q7Vr18rlcvXIYwIAAPPYGkAbN27U4sWLlZubqylTpqisrEzDhg3T1q1bg46fOXOmfvzjH+u+++6Tw+HokceUpPb2drW2tgZsAABg8LItgDo6OlRbWyu32/3fyYSGyu12q7q6uk8fs7i4WFFRUf5t3Lhx1/X8AABgYLAtgM6ePavOzk7FxMQE7I+JiZHX6+3TxywoKFBLS4t/a2xsvK7nBwAAA0O43RPoDxwOx1VfUgMAAIOPbVeAoqOjFRYWpqampoD9TU1NV32Dsx2PCQAABh/bAigiIkIzZsxQVVWVf5/P51NVVZXS0tL6zWMCAIDBx9aXwDwej3JycpScnKyUlBSVlJSora1Nubm5kqSFCxdqzJgxKi4ulvSfNzm//vrr/j+fOnVKx48f1w033KDPfe5zXXpMAAAAWwNo3rx5OnPmjAoLC+X1epWUlKSKigr/m5gbGhoUGvrfi1SnT5/WtGnT/LfXr1+v9evXKz09XYcOHerSYwIAAIRYlmXZPYn+prW1VVFRUWppaVFkZKTd04GN4lbsu+b976zNDBjzztrM3p4SAOAquvPvt+0/hQEAANDXCCAAAGAcAggAABiHAAIAAMYhgAAAgHH4KQzg//vfT3xd7ye6eupxAAC9hytAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDi2B1Bpaani4uLkdDqVmpqqmpqaa47fuXOn4uPj5XQ6lZCQoP379wfcf+HCBeXl5Wns2LEaOnSopkyZorKyst48BQAAMMDYGkA7duyQx+NRUVGR6urqlJiYqIyMDDU3Nwcdf+TIEc2fP1+LFi1SfX29srKylJWVpRMnTvjHeDweVVRU6Be/+IXeeOMNLV++XHl5edqzZ09fnRYAAOjnbA2gjRs3avHixcrNzfVfqRk2bJi2bt0adPymTZs0d+5c5efna/LkyVqzZo2mT5+uzZs3+8ccOXJEOTk5mjNnjuLi4vTQQw8pMTHxE68sAQAAc9gWQB0dHaqtrZXb7f7vZEJD5Xa7VV1dHfSY6urqgPGSlJGRETB+1qxZ2rNnj06dOiXLsnTw4EG9+eab+vKXv3zVubS3t6u1tTVgAwAAg5dtAXT27Fl1dnYqJiYmYH9MTIy8Xm/QY7xe7yeOf+qppzRlyhSNHTtWERERmjt3rkpLS/XFL37xqnMpLi5WVFSUfxs3btynODMAANDf2f4m6J721FNP6ZVXXtGePXtUW1urDRs2aOnSpfr9739/1WMKCgrU0tLi3xobG/twxgAAoK+F2/XE0dHRCgsLU1NTU8D+pqYmuVyuoMe4XK5rjv/www+1cuVK7dq1S5mZmZKk2267TcePH9f69euvePnsIw6HQw6H49OeEgAAGCBsuwIUERGhGTNmqKqqyr/P5/OpqqpKaWlpQY9JS0sLGC9JlZWV/vGXL1/W5cuXFRoaeFphYWHy+Xw9fAYAAGCgsu0KkPSfj6zn5OQoOTlZKSkpKikpUVtbm3JzcyVJCxcu1JgxY1RcXCxJWrZsmdLT07VhwwZlZmZq+/btOnbsmLZs2SJJioyMVHp6uvLz8zV06FBNmDBBhw8f1vPPP6+NGzfadp7on+JW7LPlud5Zm9lnzwsACM7WAJo3b57OnDmjwsJCeb1eJSUlqaKiwv9G54aGhoCrObNmzVJ5eblWrVqllStXatKkSdq9e7emTp3qH7N9+3YVFBRowYIFev/99zVhwgQ9/vjjWrJkSZ+fHwAA6J9sDSBJysvLU15eXtD7Dh06dMW+7OxsZWdnX/XxXC6Xnn322Z6aHgAAGIQG3afAAAAAPgkBBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAME63AqiwsFAXL1703/7ggw96fEIAAAC9rVsB9Pjjj+vChQv+2xMmTNA//vGPHp8UAABAb+pWAFmWdc3bAAAAA4Ht7wEqLS1VXFycnE6nUlNTVVNTc83xO3fuVHx8vJxOpxISErR///4rxrzxxhv62te+pqioKA0fPlwzZ85UQ0NDb50CAAAYYLoVQCEhITp//rxaW1vV0tKikJAQXbhwQa2trQFbV+3YsUMej0dFRUWqq6tTYmKiMjIy1NzcHHT8kSNHNH/+fC1atEj19fXKyspSVlaWTpw44R/z97//XbNnz1Z8fLwOHTqk1157TatXr5bT6ezOqQIAgEEsvDuDLcvS5z//+YDb06ZNC7gdEhKizs7OLj3exo0btXjxYuXm5kqSysrKtG/fPm3dulUrVqy4YvymTZs0d+5c5efnS5LWrFmjyspKbd68WWVlZZKkH/zgB7rrrru0bt06/3G33HJLd04TAAAMct0KoIMHD/bYE3d0dKi2tlYFBQX+faGhoXK73aqurg56THV1tTweT8C+jIwM7d69W5Lk8/m0b98+PfLII8rIyFB9fb0mTpyogoICZWVlXXUu7e3tam9v99/uzlUsAAAw8HQrgNLT03vsic+ePavOzk7FxMQE7I+JidHf/va3oMd4vd6g471erySpublZFy5c0Nq1a/XYY4/piSeeUEVFhe655x4dPHjwqvMvLi7WD3/4wx44K/RncSv2+f/8ztpMG2cSOBfJ/vkAgGlsfxN0T/L5fJKku+++Ww8//LCSkpK0YsUKfeUrX/G/RBZMQUGBWlpa/FtjY2NfTRkAANigW1eAwsLCujSuK+8Bio6OVlhYmJqamgL2NzU1yeVyBT3G5XJdc3x0dLTCw8M1ZcqUgDGTJ0/Wyy+/fNW5OBwOORyOT5wzAAAYHLr9JugJEyYoJycn4M3P1yMiIkIzZsxQVVWV//05Pp9PVVVVysvLC3pMWlqaqqqqtHz5cv++yspKpaWl+R9z5syZOnnyZMBxb775piZMmPCp5gsAAAaPbgVQTU2NnnnmGW3atEkTJ07Ugw8+qAULFuimm266rif3eDzKyclRcnKyUlJSVFJSora2Nv+nwhYuXKgxY8aouLhYkrRs2TKlp6drw4YNyszM1Pbt23Xs2DFt2bLF/5j5+fmaN2+evvjFL+pLX/qSKioqtHfvXh06dOi65ggAAAafbr0HKDk5WU8//bTee+89eTwe7dq1S2PHjtV9992nysrKbj/5vHnztH79ehUWFiopKUnHjx9XRUWF/43ODQ0Neu+99/zjZ82apfLycm3ZskWJiYn69a9/rd27d2vq1Kn+MV//+tdVVlamdevWKSEhQT//+c/1m9/8RrNnz+72/AAAwODUrStAH3E6nbr//vt1//336+2339aiRYs0d+5cnTlzRjfffHO3HisvL++qL3kFu2qTnZ2t7Ozsaz7mgw8+qAcffLBb8wAAAOa4rgCSpH/+85967rnn9Nxzz+nixYvKz89XZGRkT84NAACgV3QrgDo6OrRr1y4988wz+uMf/6g777xTJSUluvPOO7v8CTEAAAC7dSuARo8erRtvvFE5OTn6yU9+olGjRkmS2traAsZxJQgAAPRn3QqgDz74QB988IHWrFmjxx577Ir7u/tbYAAAAHaw7bfAAAAA7NKtAJo9e7bWr1+vPXv2qKOjQ3fccYeKioo0dOjQ3pofAABAj+vW9wD96Ec/0sqVK3XDDTdozJgx2rRpk5YuXdpbcwMAAOgV3Qqg559/Xj/5yU/00ksvaffu3dq7d6+2bdvm/xFSAACAgaBbAdTQ0KC77rrLf9vtdiskJESnT5/u8YkBAAD0lm4F0L///W85nc6AfUOGDNHly5d7dFIAAAC9qdu/Bv+tb31LDofDv+/SpUtasmSJhg8f7t/3wgsv9NwMAQAAeli3AignJ+eKfffff3+PTQYAAKAvdCuAnn322d6aBwAAQJ/p1nuAAAAABgMCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHHC7Z4A0BviVuwLuP3O2kybZtJ1H5/zQJgvAAxkXAECAADGIYAAAIBx+kUAlZaWKi4uTk6nU6mpqaqpqbnm+J07dyo+Pl5Op1MJCQnav3//VccuWbJEISEhKikp6eFZAwCAgcr2ANqxY4c8Ho+KiopUV1enxMREZWRkqLm5Oej4I0eOaP78+Vq0aJHq6+uVlZWlrKwsnThx4oqxu3bt0iuvvKLY2NjePg0AADCA2B5AGzdu1OLFi5Wbm6spU6aorKxMw4YN09atW4OO37Rpk+bOnav8/HxNnjxZa9as0fTp07V58+aAcadOndJ3v/tdbdu2TUOGDLnmHNrb29Xa2hqwAQCAwcvWAOro6FBtba3cbrd/X2hoqNxut6qrq4MeU11dHTBekjIyMgLG+3w+PfDAA8rPz9cXvvCFT5xHcXGxoqKi/Nu4ceOu84wAAMBAYGsAnT17Vp2dnYqJiQnYHxMTI6/XG/QYr9f7ieOfeOIJhYeH63vf+16X5lFQUKCWlhb/1tjY2M0zAQAAA8mg+x6g2tpabdq0SXV1dQoJCenSMQ6HQw6Ho5dnBgAA+gtbrwBFR0crLCxMTU1NAfubmprkcrmCHuNyua45/o9//KOam5s1fvx4hYeHKzw8XO+++66+//3vKy4urlfOAwAADCy2BlBERIRmzJihqqoq/z6fz6eqqiqlpaUFPSYtLS1gvCRVVlb6xz/wwAN67bXXdPz4cf8WGxur/Px8vfTSS713MgAAYMCw/SUwj8ejnJwcJScnKyUlRSUlJWpra1Nubq4kaeHChRozZoyKi4slScuWLVN6ero2bNigzMxMbd++XceOHdOWLVskSSNHjtTIkSMDnmPIkCFyuVy69dZb+/bkAABAv2R7AM2bN09nzpxRYWGhvF6vkpKSVFFR4X+jc0NDg0JD/3uhatasWSovL9eqVau0cuVKTZo0Sbt379bUqVPtOgUAADDA2B5AkpSXl6e8vLyg9x06dOiKfdnZ2crOzu7y47/zzjvXOTMAADAY2f5FiAAAAH2NAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYJ9zuCQA9IW7FPv+f31mbaeNMes7Hz0kaPOcFAP0BV4AAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHH6RQCVlpYqLi5OTqdTqampqqmpueb4nTt3Kj4+Xk6nUwkJCdq/f7//vsuXL+vRRx9VQkKChg8frtjYWC1cuFCnT5/u7dMAAAADhO0BtGPHDnk8HhUVFamurk6JiYnKyMhQc3Nz0PFHjhzR/PnztWjRItXX1ysrK0tZWVk6ceKEJOnixYuqq6vT6tWrVVdXpxdeeEEnT57U1772tb48LQAA0I/ZHkAbN27U4sWLlZubqylTpqisrEzDhg3T1q1bg47ftGmT5s6dq/z8fE2ePFlr1qzR9OnTtXnzZklSVFSUKisrde+99+rWW2/V//3f/2nz5s2qra1VQ0NDX54aAADop2wNoI6ODtXW1srtdvv3hYaGyu12q7q6Ougx1dXVAeMlKSMj46rjJamlpUUhISEaMWJE0Pvb29vV2toasAEAgMHL1gA6e/asOjs7FRMTE7A/JiZGXq836DFer7db4y9duqRHH31U8+fPV2RkZNAxxcXFioqK8m/jxo27jrMBAAADhe0vgfWmy5cv695775VlWXr66aevOq6goEAtLS3+rbGxsQ9nCQAA+lq4nU8eHR2tsLAwNTU1BexvamqSy+UKeozL5erS+I/i591339WBAweuevVHkhwOhxwOx3WeBQAAGGhsvQIUERGhGTNmqKqqyr/P5/OpqqpKaWlpQY9JS0sLGC9JlZWVAeM/ip+33npLv//97zVy5MjeOQEAADAg2XoFSJI8Ho9ycnKUnJyslJQUlZSUqK2tTbm5uZKkhQsXasyYMSouLpYkLVu2TOnp6dqwYYMyMzO1fft2HTt2TFu2bJH0n/j55je/qbq6Or344ovq7Oz0vz/o5ptvVkREhD0nCgAA+g3bA2jevHk6c+aMCgsL5fV6lZSUpIqKCv8bnRsaGhQa+t8LVbNmzVJ5eblWrVqllStXatKkSdq9e7emTp0qSTp16pT27NkjSUpKSgp4roMHD2rOnDl9cl4AAKD/sj2AJCkvL095eXlB7zt06NAV+7Kzs5WdnR10fFxcnCzL6snpAQCAQWZQfwoMAAAgGAIIAAAYhwACAADGIYAAAIBxCCAAAGCcfvEpMKA74lbsC7j9ztpMm2bS9z5+7iadNwD0NK4AAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDjhdk8A+CRxK/b5//zO2kwbZ9L/fHxtJNYHALqKK0AAAMA4BBAAADAOAQQAAIxDAAEAAOMQQAAAwDgEEAAAMA4BBAAAjEMAAQAA4xBAAADAOAQQAAAwDgEEAACMQwABAADjEEAAAMA4/Bo8+hV+3fzT+/gasn4AEBxXgAAAgHEIIAAAYBwCCAAAGIcAAgAAxiGAAACAcQggAABgHD4GD1vxke3ex1cLAMCVuAIEAACMQwABAADj9IsAKi0tVVxcnJxOp1JTU1VTU3PN8Tt37lR8fLycTqcSEhK0f//+gPsty1JhYaFGjx6toUOHyu1266233urNUwAAAAOI7QG0Y8cOeTweFRUVqa6uTomJicrIyFBzc3PQ8UeOHNH8+fO1aNEi1dfXKysrS1lZWTpx4oR/zLp16/Tkk0+qrKxMR48e1fDhw5WRkaFLly711WkhiLgV+wI22If/DgBMZ3sAbdy4UYsXL1Zubq6mTJmisrIyDRs2TFu3bg06ftOmTZo7d67y8/M1efJkrVmzRtOnT9fmzZsl/efqT0lJiVatWqW7775bt912m55//nmdPn1au3fv7sMzAwAA/ZWtnwLr6OhQbW2tCgoK/PtCQ0PldrtVXV0d9Jjq6mp5PJ6AfRkZGf64efvtt+X1euV2u/33R0VFKTU1VdXV1brvvvuueMz29na1t7f7b7e0tEiSWltbr/vcTDO16KWA2yd+mBGw78QPM+RrvxgwprW1NWDf/97u6zH/63rG2H0O1zumK//9AKC/++jfbcuyPnmwZaNTp05ZkqwjR44E7M/Pz7dSUlKCHjNkyBCrvLw8YF9paak1atQoy7Is609/+pMlyTp9+nTAmOzsbOvee+8N+phFRUWWJDY2NjY2NrZBsDU2Nn5ig/A9QJIKCgoCrir5fD69//77GjlypEJCQnr8+VpbWzVu3Dg1NjYqMjKyxx8f/8E69z7WuG+wzr2PNe4bvb3OlmXp/Pnzio2N/cSxtgZQdHS0wsLC1NTUFLC/qalJLpcr6DEul+ua4z/636amJo0ePTpgTFJSUtDHdDgccjgcAftGjBjRnVO5LpGRkfwfrQ+wzr2PNe4brHPvY437Rm+uc1RUVJfG2fom6IiICM2YMUNVVVX+fT6fT1VVVUpLSwt6TFpaWsB4SaqsrPSPnzhxolwuV8CY1tZWHT169KqPCQAAzGL7S2Aej0c5OTlKTk5WSkqKSkpK1NbWptzcXEnSwoULNWbMGBUXF0uSli1bpvT0dG3YsEGZmZnavn27jh07pi1btkiSQkJCtHz5cj322GOaNGmSJk6cqNWrVys2NlZZWVl2nSYAAOhHbA+gefPm6cyZMyosLJTX61VSUpIqKioUExMjSWpoaFBo6H8vVM2aNUvl5eVatWqVVq5cqUmTJmn37t2aOnWqf8wjjzyitrY2PfTQQzp37pxmz56tiooKOZ3OPj+/YBwOh4qKiq542Q09i3Xufaxx32Cdex9r3Df60zqHWFZXPisGAAAweNj+RYgAAAB9jQACAADGIYAAAIBxCCAAAGAcAsgGpaWliouLk9PpVGpqqmpqauye0oBVXFysmTNn6sYbb9SoUaOUlZWlkydPBoy5dOmSli5dqpEjR+qGG27QN77xjSu+TBNdt3btWv/XTXyENe4Zp06d0v3336+RI0dq6NChSkhI0LFjx/z3W5alwsJCjR49WkOHDpXb7dZbb71l44wHls7OTq1evVoTJ07U0KFDdcstt2jNmjUBvxvFGnffH/7wB331q19VbGysQkJCrvjh8a6s6fvvv68FCxYoMjJSI0aM0KJFi3ThwoVenTcB1Md27Nghj8ejoqIi1dXVKTExURkZGWpubrZ7agPS4cOHtXTpUr3yyiuqrKzU5cuX9eUvf1ltbW3+MQ8//LD27t2rnTt36vDhwzp9+rTuueceG2c9cL366qv66U9/qttuuy1gP2v86X3wwQe6/fbbNWTIEP3ud7/T66+/rg0bNuimm27yj1m3bp2efPJJlZWV6ejRoxo+fLgyMjJ06dIlG2c+cDzxxBN6+umntXnzZr3xxht64okntG7dOj311FP+Maxx97W1tSkxMVGlpaVB7+/Kmi5YsEB//etfVVlZqRdffFF/+MMf9NBDD/XuxD/x18LQo1JSUqylS5f6b3d2dlqxsbFWcXGxjbMaPJqbmy1J1uHDhy3Lsqxz585ZQ4YMsXbu3Okf88Ybb1iSrOrqarumOSCdP3/emjRpklVZWWmlp6dby5YtsyyLNe4pjz76qDV79uyr3u/z+SyXy2X9+Mc/9u87d+6c5XA4rF/+8pd9McUBLzMz03rwwQcD9t1zzz3WggULLMtijXuCJGvXrl3+211Z09dff92SZL366qv+Mb/73e+skJAQ69SpU702V64A9aGOjg7V1tbK7Xb794WGhsrtdqu6utrGmQ0eLS0tkqSbb75ZklRbW6vLly8HrHl8fLzGjx/PmnfT0qVLlZmZGbCWEmvcU/bs2aPk5GRlZ2dr1KhRmjZtmn72s5/573/77bfl9XoD1jkqKkqpqamscxfNmjVLVVVVevPNNyVJf/7zn/Xyyy/rzjvvlMQa94aurGl1dbVGjBih5ORk/xi3263Q0FAdPXq01+Zm+zdBm+Ts2bPq7Oz0f8v1R2JiYvS3v/3NplkNHj6fT8uXL9ftt9/u/2Zwr9eriIiIK37cNiYmRl6v14ZZDkzbt29XXV2dXn311SvuY417xj/+8Q89/fTT8ng8WrlypV599VV973vfU0REhHJycvxrGezvD9a5a1asWKHW1lbFx8crLCxMnZ2devzxx7VgwQJJYo17QVfW1Ov1atSoUQH3h4eH6+abb+7VdSeAMGgsXbpUJ06c0Msvv2z3VAaVxsZGLVu2TJWVlf3m52QGI5/Pp+TkZP3oRz+SJE2bNk0nTpxQWVmZcnJybJ7d4PCrX/1K27ZtU3l5ub7whS/o+PHjWr58uWJjY1ljA/ESWB+Kjo5WWFjYFZ+OaWpqksvlsmlWg0NeXp5efPFFHTx4UGPHjvXvd7lc6ujo0Llz5wLGs+ZdV1tbq+bmZk2fPl3h4eEKDw/X4cOH9eSTTyo8PFwxMTGscQ8YPXq0pkyZErBv8uTJamhokCT/WvL3x/XLz8/XihUrdN999ykhIUEPPPCAHn74Yf+PbbPGPa8ra+pyua74INC///1vvf/++7267gRQH4qIiNCMGTNUVVXl3+fz+VRVVaW0tDQbZzZwWZalvLw87dq1SwcOHNDEiRMD7p8xY4aGDBkSsOYnT55UQ0MDa95Fd9xxh/7yl7/o+PHj/i05OVkLFizw/5k1/vRuv/32K77C4c0339SECRMkSRMnTpTL5QpY59bWVh09epR17qKLFy8G/Li2JIWFhcnn80lijXtDV9Y0LS1N586dU21trX/MgQMH5PP5lJqa2nuT67W3VyOo7du3Ww6Hw3ruuees119/3XrooYesESNGWF6v1+6pDUjf/va3raioKOvQoUPWe++9598uXrzoH7NkyRJr/Pjx1oEDB6xjx45ZaWlpVlpamo2zHvg+/ikwy2KNe0JNTY0VHh5uPf7449Zbb71lbdu2zRo2bJj1i1/8wj9m7dq11ogRI6zf/va31muvvWbdfffd1sSJE60PP/zQxpkPHDk5OdaYMWOsF1980Xr77betF154wYqOjrYeeeQR/xjWuPvOnz9v1dfXW/X19ZYka+PGjVZ9fb317rvvWpbVtTWdO3euNW3aNOvo0aPWyy+/bE2aNMmaP39+r86bALLBU089ZY0fP96KiIiwUlJSrFdeecXuKQ1YkoJuzz77rH/Mhx9+aH3nO9+xbrrpJmvYsGHW17/+deu9996zb9KDwP8GEGvcM/bu3WtNnTrVcjgcVnx8vLVly5aA+30+n7V69WorJibGcjgc1h133GGdPHnSptkOPK2trdayZcus8ePHW06n0/rsZz9r/eAHP7Da29v9Y1jj7jt48GDQv4dzcnIsy+ramv7rX/+y5s+fb91www1WZGSklZuba50/f75X5x1iWR/7CkwAAAAD8B4gAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIABGmDNnjpYvX273NAD0EwQQAAAwDgEEAACMQwABMNK+ffsUFRWlbdu22T0VADYIt3sCANDXysvLtWTJEpWXl+srX/mK3dMBYAOuAAEwSmlpqb7zne9o7969xA9gMK4AATDGr3/9azU3N+tPf/qTZs6cafd0ANiIK0AAjDFt2jR95jOf0datW2VZlt3TAWAjAgiAMW655RYdPHhQv/3tb/Xd737X7ukAsBEvgQEwyuc//3kdPHhQc+bMUXh4uEpKSuyeEgAbEEAAjHPrrbfqwIEDmjNnjsLCwrRhwwa7pwSgj4VYvBAOAAAMw3uAAACAcQggAABgHAIIAAAYhwACAADGIYAAAIBxCCAAAGAcAggAABiHAAIAAMYhgAAAgHEIIAAAYBwCCAAAGOf/AcOTXZUQu5f6AAAAAElFTkSuQmCC",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.bar(ks, pmf_k)\n",
+    "plt.xlabel('k')\n",
+    "plt.ylabel('PMF');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c3071e06",
+   "metadata": {},
+   "source": [
+    "Because we used `binom` to compute the cube, we should not be surprised to find that this slice from the cube is a binomial PMF.\n",
+    "But just to make sure, we can use `binom` again to confirm it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "46e816a6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pmf_binom = binom.pmf(ks, n, p/100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "5b889032",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.allclose(pmf_k, pmf_binom)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "56cd2854",
+   "metadata": {},
+   "source": [
+    "So we can think of the cube as a collection of binomial PMFs.\n",
+    "But we can also think of it as a joint distribution of $k$, $n$, and $p$, which raises the question: what do we get if we select a vector along the $n$ and $p$ axes?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5ba8700f",
+   "metadata": {},
+   "source": [
+    "## The negative binomial distribution\n",
+    "\n",
+    "Suppose we plan to run Bernoulli trials with probability $p$ until we see $k$ successes.\n",
+    "How many trials will it take?\n",
+    "\n",
+    "We can answer this question by selecting a vector from the cube along the $n$ axis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "517c99e8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "k = 25\n",
+    "p = 50\n",
+    "pmf_n = cube[k, :, p].copy()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1d7679c1",
+   "metadata": {},
+   "source": [
+    "The result is close to the answer we want, but there's something we have to fix.\n",
+    "Remember that the values in the cube come from the binomial PMF, which looks like this.\n",
+    "\n",
+    "$$Pr(k; n, p) = \\binom{n}{k} p^{k} (1-p)^{n-k}$$\n",
+    "\n",
+    "The first term is the binomial coefficient, which indicates that there are $n$ places we could find $k$ successes.\n",
+    "But if we keep running trials until we see $k$ successes, we know the last trial will be a success, which means there are only $n-1$ places we could find the other $k-1$ successes.\n",
+    "\n",
+    "So we have to adjust the values from the cube by dividing the elements by $n/k$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "a9d6f7f5",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/tmp/ipykernel_125703/780240301.py:1: RuntimeWarning: invalid value encountered in divide\n",
+      "  pmf_n /= (ns / k)\n"
+     ]
+    }
+   ],
+   "source": [
+    "pmf_n /= (ns / k)\n",
+    "pmf_n[0] = 0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ea6153c2",
+   "metadata": {},
+   "source": [
+    "Now we have to normalize the results to get a proper PMF."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "c4ff9e52",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pmf_n /= pmf_n.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7cf1fe07",
+   "metadata": {},
+   "source": [
+    "Here's what it looks like."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "97bd5171",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAG1CAYAAAARLUsBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJ+UlEQVR4nO3de3zT9b0/8Nc3SZP0fqVJWwoULJRLabmWohM3O8vEabcdhxw2GPPMn04drhsqHAV31FXc8OAE5eDmdJtMxqZMOQ7XVUU5VBDaApX7teWSXuglbUqbNt/v7480gUiB3pJPku/r+XhE6DefpO98Hw/bF5+rpCiKAiIiIiIV0YgugIiIiMjXGICIiIhIdRiAiIiISHUYgIiIiEh1GICIiIhIdRiAiIiISHUYgIiIiEh1GICIiIhIdRiAiIiISHUYgIiIiEh1hAegtWvXYsSIETAajcjJycGuXbuu2X7Tpk3IyMiA0WhEZmYm3n///SvaHDx4EHfeeSeio6MRHh6OadOmoaqqylsfgYiIiAKM0AC0ceNGFBYWYsWKFSgrK0NWVhby8/NRW1vbY/sdO3Zg3rx5uPfee1FeXo6CggIUFBSgsrLS3eb48eO46aabkJGRgY8//hj79u3Dk08+CaPR6KuPRURERH5OEnkYak5ODqZNm4Y1a9YAAGRZRmpqKh5++GE8/vjjV7SfO3cubDYbtmzZ4r42Y8YMZGdnY926dQCAe+65ByEhIfjjH//Y77pkWca5c+cQGRkJSZL6/T5ERETkO4qioKWlBcnJydBort3Ho/NRTVew2+3Ys2cPli5d6r6m0WiQl5eH0tLSHl9TWlqKwsJCj2v5+fnYvHkzAGdw+d///V88+uijyM/PR3l5OdLS0rB06VIUFBRctZaOjg50dHS4vz579izGjRvX/w9HREREwlRXV2Po0KHXbCMsANXX18PhcMBkMnlcN5lMOHToUI+vsVgsPba3WCwAgNraWrS2tuK5557DM888g5UrV2Lr1q349re/jY8++gizZs3q8X2Liorwi1/84orr1dXViIqK6s/HIyIiIh+zWq1ITU1FZGTkddsKC0DeIMsyAOCuu+7CT3/6UwBAdnY2duzYgXXr1l01AC1dutSjZ8l1A6OiohiAiIiIAkxvpq8IC0AJCQnQarWoqanxuF5TUwOz2dzja8xm8zXbJyQkQKfTXTF8NXbsWGzfvv2qtRgMBhgMhv58DCIiIgpAwlaB6fV6TJkyBSUlJe5rsiyjpKQEubm5Pb4mNzfXoz0AFBcXu9vr9XpMmzYNhw8f9mhz5MgRDB8+fJA/AREREQUqoUNghYWFWLhwIaZOnYrp06dj9erVsNlsWLRoEQBgwYIFSElJQVFREQBg8eLFmDVrFlatWoU5c+bgrbfewu7du7F+/Xr3ey5ZsgRz587FzTffjK9+9avYunUr3nvvPXz88cciPiIRERH5IaEBaO7cuairq8Py5cthsViQnZ2NrVu3uic6V1VVeSxjmzlzJjZs2IAnnngCy5YtQ3p6OjZv3owJEya423zrW9/CunXrUFRUhJ/85CcYM2YM/va3v+Gmm27y+ecjIiIi/yR0HyB/ZbVaER0djebmZk6CJiIiChB9+f0t/CgMIiIiIl9jACIiIiLVYQAiIiIi1WEAIiIiItVhACIiIiLVYQAiIiIi1WEAIiIiItUJqsNQiUjdOrocqGvp8Lim12mQGGkUVBER+SsGICIKCs0XOzHnN5/iTOPFK557bHYGHrhllICqiMhfcQiMiILC81sP4UzjRWgkwKDTwKDTQK9z/oh7ofgwDpyzCq6QiPwJe4CIKODtPtWAN3dWAQDe/I8ZyB0VDwBQFAX3/2kPPviiBj/ftBebH7zRHYqISN34k4CIApq9S8bSt/cDAO6eMtQdfgBAkiQ8U5CJ2LAQHDhvxdqPjokqk4j8DAMQEQW09Z8cx9HaVsSH67Hs9rFXPD8k0oCnCyYAANZ+dAyVZ5t9XSIR+SEGICIKWCfrbfjNh85eneXfHIfYcH2P7e6YmIw5mUnokhX8fNNedHQ5fFkmEfkhBiAiCkiKomDZ2/th75LxlfQE3JmVfM32/3XXeMSH63HI0oKXSjgURqR2DEBEFJDe3XsOpScuwBiiwbMFmZAk6Zrt4yMMeKZ7KOyVbcdRa233RZlE5KcYgIgoIG3afQYA8P9uHoVh8WG9es03MpMwaVgMHLKC9/ad92Z5ROTnGICIKOA0t3XisxMXAADfmpTSp9cWZDvb/73i7KDXRUSBgwGIiALOh4dr0CUrGGOKxIiE8D699o6JSdBqJOw704zjda1eqpCI/B0DEBEFnA8qawAAt4039fm18REG3JyeAAD4ezl7gYjUigGIiAJKe6cD247UAQDyx5v79R4F3cNmmyvOQVGUQauNiAIHAxARBZRPjtThYqcDKTGhGJ8c1a/3+Po4E8L0WlQ1tKG8umlwCySigMAAREQB5YMvnMNfXx9nuu7S96sJ0+tw2zjn8BmHwYjUiQGIiAJGl0NGySFnAOrv8JfLXd3DYFv2nUenQx5wbUQUWBiAiChg7DrVgKa2TsSGhWDaiNgBvddXbkhAfLgeF2x2bD9aP0gVElGgYAAiooDxz+7hr7yxJui0A/vxpdNqcMfEJADAZu4JRKQ6DEBEFBAURcE/v7AAAG4b4PCXi2sY7J9f1MDW0TUo70lEgYEBiIgCwv6zzTjX3I4wvRZf6d7HZ6AmpcZgeHwYLnY6UHygZlDek4gCAwMQEQUE1/DXrNFDYAzRDsp7SpLkPkWeAYhIXRiAiCggfOAe/ur77s/XcvPoIQCAz05cgCxzU0QitWAAIiK/Z2lux9HaVmgk4GtjBjcAZQ2NQWiIFhdsdhypbRnU9yYi/8UARER+r6K6EQAwxhyF6LCQQX1vvU6Dqd1L6kuPXxjU9yYi/8UARER+z3VcRXZqjFfeP3dUPABgBwMQkWowABGR36uoagLgXLXlDTNHOVeVfXbiAhycB0SkCgxAROTXHLKC/WebAQDZw2K88j0mJEch0qBDS3sXDpyzeuV7EJF/YQAiIr92pKYFbXYHIgw6jBoS4ZXvodNqMD0tDgCw4ziPxSBSAwYgIvJrFd3zfyYOjYZW07/T33vDNQ+o9ATnARGpAQMQEfk11/wfb02AdnHNA9p1soGnwxOpAAMQEfm1Ci+vAHPJMEciNiwEbXYH9p1p8ur3IiLxGICIyG+1dnS5Nyf01gRoF41GwoyR3cNgXA5PFPQYgIjIb+070wRFAVJiQpEYafT695vJ/YCIVIMBiIj8lq+Gv1xcE6F3n25Ee6fDJ9+TiMRgACIiv+WrCdAuo4ZEYEikAfYuGeXd35uIghMDEBH5JUVRLvUAeXn+j4skSe5hsFLuB0QU1BiAiMgvnW9uR21LB3QaCROSo332fXNHch4QkRowABGRX3L1/mQkRSJUr/XZ93XtB1RR3YQ2e5fPvi8R+RYDEBH5JV9PgHZJjQuFKcqALllB5VmeC0YUrBiAiMgvlVc1AgCyU2N9+n0lSULW0BgAwN7uEEZEwccvAtDatWsxYsQIGI1G5OTkYNeuXddsv2nTJmRkZMBoNCIzMxPvv/++x/M/+MEPIEmSx2P27Nne/AhENIg6HfKlE+B93AMEAFnd37OCO0ITBS3hAWjjxo0oLCzEihUrUFZWhqysLOTn56O2trbH9jt27MC8efNw7733ory8HAUFBSgoKEBlZaVHu9mzZ+P8+fPux5///GdffBwiGgSHLS1o75QRadRhZEK4z7//pO4AxB4gouAlPAC98MIL+NGPfoRFixZh3LhxWLduHcLCwvDaa6/12P7FF1/E7NmzsWTJEowdOxZPP/00Jk+ejDVr1ni0MxgMMJvN7kdsrG+70Ymo/y6f/6Px4gnwVzNhaDQkCTjTeBH1rR0+//5E5H1CA5DdbseePXuQl5fnvqbRaJCXl4fS0tIeX1NaWurRHgDy8/OvaP/xxx8jMTERY8aMwQMPPIALF66+pLWjowNWq9XjQUTifHHOOfzlmovja1HGEIwaEgEAPBiVKEgJDUD19fVwOBwwmUwe100mEywWS4+vsVgs120/e/Zs/OEPf0BJSQlWrlyJbdu24Rvf+AYcjp63ti8qKkJ0dLT7kZqaOsBPRkQDcdjiPAB1jDlSWA2u8FVR3SysBiLyHp3oArzhnnvucf89MzMTEydOxKhRo/Dxxx/j1ltvvaL90qVLUVhY6P7aarUyBBEJoigKjta0AgBGm8QFoOzUaPyt7AznAREFKaE9QAkJCdBqtaipqfG4XlNTA7PZ3ONrzGZzn9oDwMiRI5GQkIBjx471+LzBYEBUVJTHg4jEsFjb0dLRBZ1GQpqACdAurpVge880QVEUYXUQkXcIDUB6vR5TpkxBSUmJ+5osyygpKUFubm6Pr8nNzfVoDwDFxcVXbQ8AZ86cwYULF5CUlDQ4hROR1xzp7v0ZkRAOvU7cj6gMcxT0Wg2a2jpR1dAmrA4i8g7hq8AKCwvx6quv4o033sDBgwfxwAMPwGazYdGiRQCABQsWYOnSpe72ixcvxtatW7Fq1SocOnQITz31FHbv3o2HHnoIANDa2oolS5bgs88+w6lTp1BSUoK77roLN9xwA/Lz84V8RiLqvSOu+T8Ch78AQK/TYFyysze4gsNgREFH+ByguXPnoq6uDsuXL4fFYkF2dja2bt3qnuhcVVUFjeZSTps5cyY2bNiAJ554AsuWLUN6ejo2b96MCRMmAAC0Wi327duHN954A01NTUhOTsZtt92Gp59+GgaDQchnJKLeO1LjDEDppgjBlTiX4VdUN2FvdTPuyk4RXQ4RDSJJ4eD2FaxWK6Kjo9Hc3Mz5QEQ+dtfa/8Pe6ia8PH8ybs8UO2z9TvkZ/HTjXkwZHou/PTBTaC1EdH19+f0tfAiMiMhFlhUc6+4BGu0HPUCupfCVZ5vR6ZDFFkNEg4oBiIj8xtmmi7DZHdBrNRgeL24FmMuI+HBEGXXo6JLdexMRUXBgACIiv3G01hkyRg4JR4hW/I8njUa6dDAqJ0ITBRXxP2GIiLq5lsCnC14BdjnXMBg3RCQKLgxAROQ3XCvARieKn//jcvmGiEQUPBiAiMhvXFoC7089QNEAgKO1rWjt6BJcDRENFgYgIvILsqzgWK1zCEzkIahflhhlRHK0EYoC7D/Dg1GJggUDEBH5herGNrR3yjDoNBgWFya6HA/Zw2IAcBiMKJgwABGRX3BNgB41JAJajSS4Gk+ZKTEA2ANEFEwYgIjILxzxow0Qv2x895lgB89bBVdCRIOFAYiI/II/ToB2cR2KevKCDTZOhCYKCgxAROQXXENgok+B70lChAGmKAMUBThkYS8QUTBgACIi4bocMo7XOQPQaD8MQAAwLsnZC3TgHAMQUTBgACIi4U43tMHeJSM0RIuhsaGiy+mRaxjsAOcBEQUFBiAiEu5o9/yfGxIjoPGzFWAu45OdGyJ+wR4goqDAAEREwrnm//jr8BdwaQjskKUFXQ5ZcDVENFAMQEQknD8vgXcZFheGcL0W9i4ZJ+ptosshogFiACIi4S4FIP/tAdJoJIzlRGiioMEARERCdTpknOzuUUn34x4ggBOhiYIJAxARCVXd0IZOh4LQEC1SYvxzBZiLa0foL87xSAyiQMcARERCuXp/RiSEQ5L8cwWYy7gk50qwA+esUBRFcDVENBAMQEQklCsApSX41wnwPUk3OQ9qbWzrhMXaLrocIhoABiAiEurUhe4eoPhwwZVcnzFEixuGOOcpcSI0UWBjACIioU7VtwFwDoEFgnHueUAMQESBjAGIiIRyDYGNDJAA5JoIzR4gosDGAEREwrR3OnCu+SKAAOoBSuJSeKJgwABERMJUNbRBUYBIgw7x4XrR5fSKazPEqoY2WNs7BVdDRP3FAEREwgTSEniX2HA9kqONAIBD51sEV0NE/cUARETCnLosAAWScdwQkSjgMQARkTCuJfBpAReALm2ISESBiQGIiIQ5URc4myBejhOhiQIfAxARCRNImyBezrUU/mhNK+xdsuBqiKg/GICISIg2exdqrB0AAm8IbGhsKCINOtgdMo7XtYouh4j6gQGIiIRw7QAdExaCmLDAWALvIkkSxpgjAQCHLVwJRhSIGICISIhAnQDtkpHkDECHGICIAhIDEBEJ4T4FPsDm/7iMMTvnAR2ycCI0USBiACIiIU4G6B5ALmM5BEYU0BiAiEiIQN0E0WV0dwA639yO5jYeiUEUaBiAiEgI1xygQDkF/suijCFIiQkFwGEwokDEAEREPmdt70R9qx1A4PYAAUCGmROhiQIVAxAR+Zxr+CshwoAIg05wNf3HlWBEgYsBiIh8zr0CLMCOwPgyrgQjClwMQETkc65NEAPtCIwvc60EO2JpgSwrgqshor5gACIin3NvgjgksAPQiIRw6LUa2OwOnGm8KLocIuoDBiAi8rlA3wTRJUSrwQ2JEQA4DEYUaBiAiMjnAn0TxMtlcENEooDEAEREPtVos6P5onPjwECfAwRwJRhRoPKLALR27VqMGDECRqMROTk52LVr1zXbb9q0CRkZGTAajcjMzMT7779/1bb3338/JEnC6tWrB7lqIuqPk93zf8xRRoTqtYKrGTiuBCMKTMID0MaNG1FYWIgVK1agrKwMWVlZyM/PR21tbY/td+zYgXnz5uHee+9FeXk5CgoKUFBQgMrKyivavvPOO/jss8+QnJzs7Y9BRL10qj6wT4H/MtdKsJP1NrR3OgRXQ0S9JTwAvfDCC/jRj36ERYsWYdy4cVi3bh3CwsLw2muv9dj+xRdfxOzZs7FkyRKMHTsWTz/9NCZPnow1a9Z4tDt79iwefvhhvPnmmwgJCfHFRyGiXgj0M8C+bEikAbFhIZAV4Fhtq+hyiKiXhAYgu92OPXv2IC8vz31No9EgLy8PpaWlPb6mtLTUoz0A5Ofne7SXZRnf//73sWTJEowfP/66dXR0dMBqtXo8iMg7Tl1w7QEU2JsgukiShIzuYbCD5/mzgyhQCA1A9fX1cDgcMJlMHtdNJhMsFkuPr7FYLNdtv3LlSuh0OvzkJz/pVR1FRUWIjo52P1JTU/v4SYiot6oanAFoeJAEIAAYw5VgRAFH+BDYYNuzZw9efPFFvP7665AkqVevWbp0KZqbm92P6upqL1dJpF7V3QEoNS54AtBYrgQjCjhCA1BCQgK0Wi1qamo8rtfU1MBsNvf4GrPZfM32n376KWprazFs2DDodDrodDqcPn0aP/vZzzBixIge39NgMCAqKsrjQUSDr7WjCxdszlPggykAXVoJxgBEFCiEBiC9Xo8pU6agpKTEfU2WZZSUlCA3N7fH1+Tm5nq0B4Di4mJ3++9///vYt28fKioq3I/k5GQsWbIEH3zwgfc+DBFdV1X3/J/YsBBEGYNnccJoUwQkCahv7UB9a4focoioF3SiCygsLMTChQsxdepUTJ8+HatXr4bNZsOiRYsAAAsWLEBKSgqKiooAAIsXL8asWbOwatUqzJkzB2+99RZ2796N9evXAwDi4+MRHx/v8T1CQkJgNpsxZswY3344IvLgmv8zLAg2QLxcmF6H4XFhOHWhDYctLUi4wSC6JCK6DuEBaO7cuairq8Py5cthsViQnZ2NrVu3uic6V1VVQaO51FE1c+ZMbNiwAU888QSWLVuG9PR0bN68GRMmTBD1EYiol1zzf4YF0fCXS4Y5CqcutOHgeStuvCFBdDlEdB3CAxAAPPTQQ3jooYd6fO7jjz++4trdd9+Nu+++u9fvf+rUqX5WRkSDyd0DFBcquJLBN8Ycia1fWLgSjChABN0qMCLyX6eDugeoeyl8DQMQUSBgACIin7k0BBZcc4CAS3sBHalpgUNWBFdDRNfDAEREPuGQFZxpdE2CDr4eoOHx4TDoNGjvlN1DfUTkvxiAiMgnzjdfRKdDQYhWgjnKKLqcQafVSEg3RQDgjtBEgYABiIh8wtUrMjQ2DFpN73ZpDzRjTM4NERmAiPwfAxAR+UQwL4F3uTQRmoeiEvk7BiAi8okqFQSg0WaeCUYUKBiAiMgnTl8I/gDk6gE6VW9De6dDcDVEdC0MQETkE+4hsCBcAeaSGGlATFgIZAU4VtsquhwiugYGICLyCTUMgUmShDGm7nlAHAYj8msMQETkddb2TjS2dQIAUoM4AAGXNkTkjtBE/o0BiIi8rqp7/k98uB4RBr84gtBrxnAiNFFAYAAiIq9Tw/wfF9dE6CMMQER+jQGIiLxODfN/XEZ3zwGyWNvR3D3sR0T+hwGIiLwumE+B/7JIYwhSYkIBAIcs3BCRyF8xABGR17mGwIJ9ArQLJ0IT+T8GICLyOtcQ2HCVBSBOhCbyXwxARORVXQ4ZZxsvAlDHJGiAE6GJAgEDEBF51fnmdnTJCvRaDUyRRtHl+MTlQ2CKogiuhoh6wgBERF7lGv4aGhcKjUYSXI1vjEyIgE4joaW9C+ea20WXQ0Q9YAAiIq9S2/wfANDrNBg5JBwAcJgrwYj8EgMQEXmVmvYAutwYcxQAToQm8lcMQETkVa5jMNSyBN6FE6GJ/BsDEBF5lVp7gFw7QrMHiMg/MQARkVe55wDFhwuuxLdcPUDH61rR6ZAFV0NEX8YARERe03yxE80XnedhpcaFCq7Gt1JiQhGu16LToeBkvU10OUT0JQxAROQ1riMwEiL0CNPrBFfjWxqNxB2hifwYAxAReU2Vys4A+zL3SrDzXApP5G8YgIjIa9yHoMaqMwCNTWIPEJG/YgAiIq+pblTnCjCXMd0rwQ4zABH5HQYgIvKaqgbnIahqmwDtktE9BHa26aJ7MjgR+QcGICLymjMqnwMUHRaC5GjnAbBHatgLRORPGICIyCtkWcGZxu4eIJXOAQIunQzPidBE/oUBiIi8oqalHXaHDK1GQlJ3L4gaZSQ5h8EOch4QkV9hACIir3CdAZYSEwqdVr0/alw7QnMiNJF/Ue9PJSLyqupGdU+AdnFNhD5saYGiKIKrISIXBiAi8gq1HoL6ZSOHhCNEK6G1o8s9J4qIxGMAIiKvcK0AG6riCdAAEKLVYNSQCADcEJHInzAAEZFXqH0TxMuNTXINg3ElGJG/YAAiIq9Q+zlgl3MthedKMCL/wQBERIOuvdOBGmsHAPYAAZdWgnEvICL/wQBERIPONdk3XK9FbFiI4GrEcw2Bnay3ob3TIbgaIgIYgIjIC1zzf1LjwiBJkuBqxEuMNCAmLASyAhyrbRVdDhGBAYiIvKCa8388SJJ0aRiM84CI/AIDEBENOncAUvkS+Mu5NkTkPCAi/8AARESDrrrBOQdomMp3gb4ce4CI/AsDEBENOi6Bv5LrUFQGICL/wABERINKURT3EBiXwF8y2hQBSQLqWztQ39ohuhwi1fOLALR27VqMGDECRqMROTk52LVr1zXbb9q0CRkZGTAajcjMzMT777/v8fxTTz2FjIwMhIeHIzY2Fnl5edi5c6c3PwIRdWu+2ImWji4APAbjcmF6HYZ3B0KeDE8knvAAtHHjRhQWFmLFihUoKytDVlYW8vPzUVtb22P7HTt2YN68ebj33ntRXl6OgoICFBQUoLKy0t1m9OjRWLNmDfbv34/t27djxIgRuO2221BXV+erj0WkWq75PwkRBoTqtYKr8S/uHaE5EZpIuD4FoOXLl6Otrc39dWNj44ALeOGFF/CjH/0IixYtwrhx47Bu3TqEhYXhtdde67H9iy++iNmzZ2PJkiUYO3Ysnn76aUyePBlr1qxxt/n3f/935OXlYeTIkRg/fjxeeOEFWK1W7Nu3b8D1EtG1XToFnhOgv8y9Eow9QETC9SkAPfvss2htvbSJ1/Dhw3HixIl+f3O73Y49e/YgLy/vUkEaDfLy8lBaWtrja0pLSz3aA0B+fv5V29vtdqxfvx7R0dHIysrqsU1HRwesVqvHg4j65/JNEMmTa0do9gARidenAKQoyjW/7qv6+no4HA6YTCaP6yaTCRaLpcfXWCyWXrXfsmULIiIiYDQa8d///d8oLi5GQkJCj+9ZVFSE6Oho9yM1NXUAn4pI3TgB+urGJzsD0NGaVti7ZMHVEKmb8DlA3vLVr34VFRUV2LFjB2bPno3vfve7V51XtHTpUjQ3N7sf1dXVPq6WKHhUcRPEqxoaG4pIgw52h4zjdTwSg0ikPgUgSZLQ0tICq9WK5uZmSJKE1tbWfg8fJSQkQKvVoqamxuN6TU0NzGZzj68xm829ah8eHo4bbrgBM2bMwO9+9zvodDr87ne/6/E9DQYDoqKiPB5E1D+ug1A5BHYlSZI4DEbkJ/o8BDZ69GjExsYiLi4Ora2tmDRpEmJjYxEbG4uYmBjExsb2+v30ej2mTJmCkpIS9zVZllFSUoLc3NweX5Obm+vRHgCKi4uv2v7y9+3o4N4bRN7kkBWccc8B4iTonozrHgY7cI4BiEgkXV8af/TRR4NeQGFhIRYuXIipU6di+vTpWL16NWw2GxYtWgQAWLBgAVJSUlBUVAQAWLx4MWbNmoVVq1Zhzpw5eOutt7B7926sX78eAGCz2fDss8/izjvvRFJSEurr67F27VqcPXsWd99996DXT0SX1Fjb0elQoNNISIpmAOrJuO4eoAPsASISqk8BaNasWYNewNy5c1FXV4fly5fDYrEgOzsbW7dudU90rqqqgkZzqaNq5syZ2LBhA5544gksW7YM6enp2Lx5MyZMmAAA0Gq1OHToEN544w3U19cjPj4e06ZNw6efforx48cPev1EdIlr/k9KbCi0GklwNf7J3QN03gpFUSBJvE9EIkjKQJdyBSGr1Yro6Gg0NzdzPhBRH2zaXY0lf92Hm25IwJ/+I0d0OX6pvdOB8Ss+gENWsOPxryE5hj1lRIOlL7+/+zQHSKvV9upBROpUzQnQ12UM0eKGIREAOBGaSKQ+DYEpioLhw4dj4cKFmDRpkrdqIqIAVd3ACdC9MS45CodrWnDgnBW3jjVd/wVENOj6FIB27dqF3/3ud3jxxReRlpaGH/7wh5g/f36fVn4RUfA6fcEGABgeFy64Ev82LikK75Sf5URoIoH6NAQ2depUvPLKKzh//jwKCwvxzjvvYOjQobjnnntQXFzsrRqJKEBUdR+Eyl2gr+3yidBEJEa/doI2Go343ve+h5KSElRWVqK2thazZ89GQ0PDYNdHRAGizd6F+lbnXlvD4hmArsW1GeLpC21o7egSXA2ROvX7KIwzZ87gmWeewde//nUcOnQIS5Ys4YopIhVzLYGPDg1BdGiI4Gr8W1y4HuYoIwDgEHuBiIToUwCy2+3YuHEjbrvtNqSnp6OsrAyrV69GdXU1nnvuOeh0fZpSRERBpOqCMwANZ+9Pr3AYjEisPiWWpKQkREZGYuHChXj55ZeRmJgIwLn78uXYE0SkPu5DUDn/p1fGJUXhw0O1PBKDSJA+BaDGxkY0Njbi6aefxjPPPHPF865dTR0Ox6AVSESBwRWAOAG6d3goKpFYws8CI6LgcNo1BMYA1CuuIbBDlhZ0OWTotP2ekklE/dCnAHTTTTfh17/+Nd59913Y7XbceuutWLFiBUJDuekZkdpVsweoT4bHhSFMr0Wb3YGT9TakmyJFl0SkKn36J8cvf/lLLFu2DBEREUhJScGLL76IBx980Fu1EVGAcMgKznQfg8El8L2j0UjuYTBOhCbyvT4FoD/84Q94+eWX8cEHH2Dz5s1477338Oabb0KWZW/VR0QBwGJth90hI0QrISmaPcK9NTbJ2evDidBEvtenAFRVVYXbb7/d/XVeXh4kScK5c+cGvTAiChyuJfBDY8Og1UiCqwkc45KiAbAHiEiEPgWgrq4uGI1Gj2shISHo7Owc1KKIKLBUNTi3wuAS+L5x7wV0zgpFUQRXQ6QufT4N/gc/+AEMBoP7Wnt7O+6//36Eh186/PDtt98evAqJyO9dWgLP4a++GGOKhEYCLtjsqG3pgCnKeP0XEdGg6FMAWrhw4RXXvve97w1aMUQUmC4tgecp8H0Rqtdi1JAIHK1tReXZZgYgIh/qUwD6/e9/7606iCiAVXMX6H7LTInG0dpW7D/bjFvHmkSXQ6Qa3HmLiAbMNQTGc8D6bkKKcyJ05dlmwZUQqQsDEBENiLW9E41tzoUQ7AHqu4lDnQFo3xkGICJfYgAiogFxLYFPiNAjwtCnUXWCcyWYRgJqWzpQa20XXQ6RajAAEdGA8BT4gQnT6zBqSAQAYD+HwYh8hgGIiAbEPf+HAajfMrvnATEAEfkOAxARDYhrCTwPQe2/zO55QPs5D4jIZxiAiGhAuAR+4NgDROR7DEBENCCXlsBzE8T+4kRoIt9jACKifut0yDjbdBEAh8AGghOhiXyPAYiI+u18UzscsgKDToPESMP1X0BXlcn9gIh8igGIiPrt9GWnwGs0kuBqAlsmd4Qm8ikGICLqNy6BHzycCE3kWwxARNRvrl2guQJs4C6fCF3DidBEXscARET95uoB4gTogQvT63BDYvdEaM4DIvI6BiAi6jeeAj+4JnAYjMhnGICIqF8URXEPgbEHaHBwIjSR7zAAEVG/NLV1oqWjCwDnAA2WiUPZA0TkKwxARNQvJy84l8AnRRthDNEKriY4jEuK5kRoIh9hACKifjlV7wxAI3gExqAJ1Ws5EZrIRxiAiKhfTroCUAID0GDiRGgi32AAIqJ+cQWgkQxAg2pidwDae6ZJbCFEQY4BiIj65dQF9gB5w6RhsQCAiuomKIoiuBqi4MUARER9pigKTtY5A1BaAleADaaxSVHQ6zRoaut097IR0eBjACKiPqtr7YDN7oBG4hL4wabXadz7AZVXNYkthiiIMQARUZ+dqndugJgSGwqDjkvgB9vkYTEAgPLqRrGFEAUxBiAi6rOT9a0AuATeW1zzgNgDROQ9DEBE1Gcnu3uA0jgB2ismdfcAHbK0oM3eJbYYoiDFAEREfebaBJEByDuSokORFG2EQ1awjxsiEnkFAxAR9Rk3QfQ+Vy8Qh8GIvIMBiIj6RJYV9x5AaZwD5DWTUl3zgDgRmsgb/CIArV27FiNGjIDRaEROTg527dp1zfabNm1CRkYGjEYjMjMz8f7777uf6+zsxGOPPYbMzEyEh4cjOTkZCxYswLlz57z9MYhUwWJtR0eXDJ1GwtDYUNHlBC1XD1BZFTdEJPIG4QFo48aNKCwsxIoVK1BWVoasrCzk5+ejtra2x/Y7duzAvHnzcO+996K8vBwFBQUoKChAZWUlAKCtrQ1lZWV48sknUVZWhrfffhuHDx/GnXfe6cuPRRS0XPN/hsWFQacV/iMkaE1IiUaIVkJ9awfONF4UXQ5R0JEUwf+0yMnJwbRp07BmzRoAgCzLSE1NxcMPP4zHH3/8ivZz586FzWbDli1b3NdmzJiB7OxsrFu3rsfv8fnnn2P69Ok4ffo0hg0bdt2arFYroqOj0dzcjKioqH5+MqLg9KfPTuOJzZX4WkYiXvvBNNHlBLW71mzH3jPN+M28SbgzK1l0OUR+ry+/v4X+881ut2PPnj3Iy8tzX9NoNMjLy0NpaWmPryktLfVoDwD5+flXbQ8Azc3NkCQJMTExPT7f0dEBq9Xq8SCinnEFmO9c2g+I84CIBpvQAFRfXw+HwwGTyeRx3WQywWKx9Pgai8XSp/bt7e147LHHMG/evKumwaKiIkRHR7sfqamp/fg0ROrAQ1B95/J5QEQ0uIJ6AL+zsxPf/e53oSgKXnnllau2W7p0KZqbm92P6upqH1ZJFFhO1HMFmK9M7u4BOnCuGe2dDsHVEAUXnchvnpCQAK1Wi5qaGo/rNTU1MJvNPb7GbDb3qr0r/Jw+fRoffvjhNccCDQYDDAZDPz8FkXp0OWRUNzh3gR7BU+C9bmhsKBIi9KhvteOLc1ZMGR4ruiSioCG0B0iv12PKlCkoKSlxX5NlGSUlJcjNze3xNbm5uR7tAaC4uNijvSv8HD16FP/6178QHx/vnQ9ApDLnmtrR6VCg12mQHM0l8N4mSRKyuR8QkVcI7QECgMLCQixcuBBTp07F9OnTsXr1athsNixatAgAsGDBAqSkpKCoqAgAsHjxYsyaNQurVq3CnDlz8NZbb2H37t1Yv349AGf4+bd/+zeUlZVhy5YtcDgc7vlBcXFx0Ov1Yj4oURA44T4ENQwajSS4GnWYPDwG/zpYwx2hiQaZ8AA0d+5c1NXVYfny5bBYLMjOzsbWrVvdE52rqqqg0VzqqJo5cyY2bNiAJ554AsuWLUN6ejo2b96MCRMmAADOnj2Ld999FwCQnZ3t8b0++ugj3HLLLT75XETByLUCjKfA+w53hCbyDuH7APkj7gNE1LOn3v0Cr+84hf83aySWfmOs6HJUwdbRhcynPoCsAJ8tvRXmaKPokoj8VsDsA0REgYUrwHwv3KDDuGTnD/KdJy8IroYoeDAAEVGvneIp8ELkpDkXcuw82SC4EqLgwQBERL1i75JxptG5BH4kA5BP5aTFAQB2MQARDRoGICLqlerGNsgKEK7XYkgk983ypelpcZAk4FhtK+pbO0SXQxQUGICIqFdO1jmHv4bHh0OSuATel2LC9BhjigTAXiCiwcIARES94joDjIegiuEaBtt5ghOhiQYDAxAR9coJ9wRoHoEhQs5IToQmGkwMQETUK8dqnbtA35AYIbgSdZre3QN0yNKCRptdcDVEgY8BiIiuS1EUHK1pAQCkJ0YKrkadEiIM7vC56xR7gYgGigGIiK7rgs2OxrZOSBIwagh7gES5NA+IAYhooBiAiOi6jnT3/gyLC0OoXiu4GvW6NA+IE6GJBooBiIiuyzX/J53zf4Sa0d0DdOC8Fc0XOwVXQxTYGICI6LpcPUDpJs7/ESkxyoi0hHAoCrCb84CIBoQBiIiu62gNe4D8hXseEJfDEw0IAxARXdfR7iGw0ewBEi5nJDdEJBoMDEBEdE0XWjvQYLNzBZifcJ0MX3nOitaOLsHVEAUuBiAiuqYj3cNfqbFcAeYPkmNCkRoXCoescB4Q0QAwABHRNR2rdW2AyN4ff+HqBeI8IKL+YwAiomty9QBxBZj/mNG9H9COY/WCKyEKXAxARHRNR9kD5HduTk8AAOw724wGngtG1C8MQER0Ta4l8FwB5j8So4zIMEdCUYBPj9aJLocoIDEAEdFVXWjtwIXuHoZRieGCq6HLzRozBACw7QgDEFF/MAAR0VW59v9JjQtFmF4nuBq63KzRzgD0yZF6yLIiuBqiwMMARERXddR9BhiHv/zN1OFxCNNrUd/agYMWq+hyiAIOAxARXdVR9xlgnADtb/Q6DWaOcq4G4zAYUd8xABHRVbkPQWUPkF+6NAzGAETUVwxARHRVx9xngLEHyB/d3B2Adp9q5LEYRH3EAEREPWqw2VHf2r0CjGeA+aXh8eEYER+GLlnhpohEfcQAREQ9cs3/GRobinADV4D5K1cv0CfcD4ioTxiAiKhHR9wrwNj7489c84C2HamDonA5PFFvMQARUY+OdfcAcQdo/zZjZDz0Wg2qGy7i1IU20eUQBQwGICLqkesQ1BvYA+TXwg06TB0RCwDYdrhWcDVEgYMBiIh6dLSWZ4AFisuHwYiodxiAiOgKjTY76ls7ALAHKBC4zgX77EQD2jsdgqshCgwMQER0hS/OOY9WGBYXxhVgAWCMKRLmKCMudjpQevyC6HKIAgIDEBFdYf/ZZgBAZkq04EqoNyRJwm3jTQCArZUWwdUQBQYGICK6QmV3AJrAABQwZo83AwD+ecCCLocsuBoi/8cARERXYA9Q4JmeFofYsBA0tnVi18kG0eUQ+T0GICLy0NzWiaoG534yE1KiBFdDvaXTanDbOGcv0D84DEZ0XQxAROSh8pyz9yc1LhQxYXrB1VBfzM50BqAPvrBAlrkrNNG1MAARkQfX8NfElBixhVCf3TgqAZFGHWpbOlBW1Si6HCK/xgBERB72n+EE6ECl12mQN9a5GozDYETXxgBERB44ATqwzZ7gHAbbWmnh4ahE18AARERunAAd+GaNHoIwvRZnmy66wywRXYkBiIjcOAE68BlDtPjqmEQAHAYjuhYGICJy4/BXcOAwGNH1MQARkdt+7gAdFL6akQi9ToOT9TYcrmkRXQ6RXxIegNauXYsRI0bAaDQiJycHu3btumb7TZs2ISMjA0ajEZmZmXj//fc9nn/77bdx2223IT4+HpIkoaKiwovVEwWXSvYABYUIgw43pycAAP6xn8NgRD0RGoA2btyIwsJCrFixAmVlZcjKykJ+fj5qa2t7bL9jxw7MmzcP9957L8rLy1FQUICCggJUVla629hsNtx0001YuXKlrz4GUVBovtiJ0xe6J0AnMwAFutszkwAAf684y2Ewoh5IisD/M3JycjBt2jSsWbMGACDLMlJTU/Hwww/j8ccfv6L93LlzYbPZsGXLFve1GTNmIDs7G+vWrfNoe+rUKaSlpaG8vBzZ2dl9qstqtSI6OhrNzc2IiuJKGFKHHcfq8e+/3YmhsaHY/tjXRJdDA2Tr6MK0Z/+FNrsDm+7PxbQRcaJLIvK6vvz+FtYDZLfbsWfPHuTl5V0qRqNBXl4eSktLe3xNaWmpR3sAyM/Pv2r73uro6IDVavV4EKmNewfooez9CQbhBh3umOjsBfrL59WCqyHyP8ICUH19PRwOB0wmk8d1k8kEi6XnMWuLxdKn9r1VVFSE6Oho9yM1NXVA70cUiDgBOvjcPdX5s+x/95+HraNLcDVE/kX4JGh/sHTpUjQ3N7sf1dX81xKpDydAB5+pw2ORlhCONrsD7+8/L7ocIr8iLAAlJCRAq9WipqbG43pNTQ3MZnOPrzGbzX1q31sGgwFRUVEeDyI1sbZ34hQnQAcdSZLwb1OGAgA27T4juBoi/yIsAOn1ekyZMgUlJSXua7Iso6SkBLm5uT2+Jjc316M9ABQXF1+1PRH1jqv3Z2hsKGLDuQN0MPn25BRoJGDXqQacqreJLofIbwgdAissLMSrr76KN954AwcPHsQDDzwAm82GRYsWAQAWLFiApUuXutsvXrwYW7duxapVq3Do0CE89dRT2L17Nx566CF3m4aGBlRUVODAgQMAgMOHD6OiomLA84SIghmHv4JXUnQovpI+BADw1z3sBSJyERqA5s6di1//+tdYvnw5srOzUVFRga1bt7onOldVVeH8+Uvj1jNnzsSGDRuwfv16ZGVl4a9//Ss2b96MCRMmuNu8++67mDRpEubMmQMAuOeeezBp0qQrlskT0SUV1U0AOAE6WN091TkM9reyM3DI3BOICBC8D5C/4j5ApCaKomDKM/9Cg82Ov96fi6ncLybotHc6kPPLEjRf7MQbP5yOWaOHiC6JyCsCYh8gIvIPR2pa0WCzIzREi4lDY0SXQ15gDNHiruxkAMCm3VzlSgQwABGpXunxegDA1BGx0Ov4IyFYfbd7T6B/HqhBo80uuBoi8fjTjkjlSk9cAADMGBkvuBLypvHJURiXFAV7l4wNu6pEl0MkHAMQkYrJsoKdJxsAALmjGICCmSRJ+I+vpAEAXt9xCh1dDsEVEYnFAESkYgctVjS1dSJcr+USeBW4Y2IyzFFG1LV04O8V50SXQyQUAxCRin12wtn7My0tDiFa/jgIdnqdBotuHAEAePWTE+AiYFIz/sQjUrHS4875P7mc/6Ma83KGIcKgw9HaVnx8pE50OUTCMAARqZRDVrDzZHcA4vwf1YgyhmDuNOeKsFc/OSG4GiJxGICIVOrAOSta2rsQadBhXBI3/FSTRTeOgFYjYcfxC+5jUIjUhgGISKVKTzj3/5meFgcd5/+oytDYMMzJTAIA/PZT9gKROvGnHpFKuef/cPhLlX70lZEAgPf2nce5pouCqyHyPQYgIhXqcsj4/FQjAG6AqFaZQ6OROzIeDlnBbz89KbocIp9jACJSocpzVrR2dCE6NITzf1Ts/81y9gL9aedpnGUvEKkMAxCRCrmGv3LS4qDRSIKrIVFmjR6CGSPjYO+Sseqfh0WXQ+RTDEBEKuQ6/4vzf9RNkiQs/cZYAMA75WfxxTmuCCP1YAAiUplOh4zdp5w7QHP+D2WlxuCOiUlQFOC5fxwSXQ6RzzAAEanMzhMNaLM7EB+uxxhTpOhyyA8syR+DEK2ET4/W49Oj3B2a1IEBiEhl3tvrPAQzf4KZ838IADA8PhzfmzEcAFD0/iHIMs8Io+DHAESkIvYuGVu/sAAA7piYJLga8icPfy0dkQYdDpy34u97z4ouh8jrGICIVOT/jtWj+WInhkQakJPG+T90SVy4HvffMgoA8OsPjqC90yG4IiLvYgAiUhHX8NeczCRoOfxFX/LDG9NgjjLibNNF/Pe/joguh8irGICIVKK904F/HqgBwOEv6lmoXounCyYAcJ4UX1bVKLgiIu9hACJSiW1H6tDa0YWkaCMmD4sVXQ75qa+PM6EgOxmyAizZtJdDYRS0GICIVGLLvvMAnMNfXP1F1/LUneMxJNKA43U2DoVR0GIAIlKBNnsX/uUa/spKFlwN+buYMD1++a1MAM6hsHIOhVEQYgAiUoEPD9XiYqcDqXGhyBoaLbocCgCXD4X9nENhFIQYgIhUYMte5/DXHROTIUkc/qLeuXwo7Fcf8LBUCi4MQERBrqW9Ex8ergXA1V/UNzFhehR1D4X9bvtJvF12RnBFRIOHAYgoyP3rYA3sXTJGJoRjXFKU6HIowOSNM+HBrzo3SHz87f2cD0RBgwGIKMj9eVc1AGfvD4e/qD9+9vUxyBtrgr1Lxv/74x5YmttFl0Q0YAxAREFs54kL2HWyAXqtBvNyhokuhwKURiNh9T3ZGGOKRG1LB+77425OiqaAxwBEFMRe+vAYAODuqUORFB0quBoKZBEGHX67cCpiw0Kw70wzHv3rPigKT42nwMUARBSk9pxuxPZj9dBpJDzQfcgl0UCkxoXh5flToNNIeHfvOfzn5krIMkMQBSYGIKIg9dKHRwEA356cgqGxYYKroWCROyoeK78zEZIEbNhZhcf+tg8OhiAKQAxAREFo35kmfHy4DlqNhAe/eoPocijIfGfKUPz3d7OhkYBNe87gZ3+pQJdDFl0WUZ8wABEFIdfcn7uykjE8PlxwNRSMCial4KV5k6HVSNhccQ6LN1agkyGIAggDEFGQOXjeiuIDNZAk4Mfs/SEvmjMxCS/Pn4wQrYT/3Xce976xG402u+iyiHqFAYgoyKzp7v25Y2IybkiMEFwNBbv88Wb8z/enwKDT4JMjdbjjpe3YW90kuiyi62IAIgoin524gPcrned+PcTeH/KRr2WY8M6Pb8Tw+DCcbbqIu9eV4s2dp7lMnvwaAxBRkKhr6cBP/lwORQG+O3UoxpgjRZdEKjIuOQrvPnSTc8doh4z/fKcSP/vLXjRf7BRdGlGPGICIgoBDVvDTjRWobelAemIEnrpzvOiSSIWiQ0Ow/vtT8NjsDGgk4O3ys7h11cd4u+wMe4PI7zAAEQWBNR8ew/Zj9QgN0eLl+ZMRpteJLolUStO98eZb9+Vi1JBw1LfaUfiXvZi7/jMcqWkRXR6RGwMQUYDbcaweq0uOAACeKZiAdBOHvki86Wlx+Mfim/Ho7DEwhmiw62QDbn/xUzy5uRLVDW2iyyOCpLBf8gpWqxXR0dFobm5GVFSU6HKIrqq2pR23v7gd9a0d+O7UoXj+37JEl0R0hTONbfiv9w7gnwdqAABajYS7spJx/y2jMJqBnQZRX35/MwD1gAGIAsGx2hbc98c9OFFnwxhTJDY/eCNC9VrRZRFd1Y7j9Xj5o+PYfqzefS1vbCK+OzUVt4xJhF7HQQkaGAagAWIAIn/3wRcWFG6sgM3uQFK0EW/+Rw5GDuGePxQY9lY34ZWPj+ODAxa4fgPFhoXgm1nJ+NakFGSnxkCSJLFFUkBiABogBiDyVw5Zwep/HXEfdZGTFoe18ycjIcIguDKivjtW24qNn1dhc8U51LV0uK8nRRtxc/oQzBozBDeOSkB0WIjAKimQMAANEAMQ+aPyqka8UHwEnx51Dh/88MY0LL09AyFaDhtQYOtyyPi/4xfwTtkZfPBFDS52OtzPaSRg4tAYTB4Wi6zUaExKjUVqXCh7iKhHAReA1q5di1/96lewWCzIysrCSy+9hOnTp1+1/aZNm/Dkk0/i1KlTSE9Px8qVK3H77be7n1cUBStWrMCrr76KpqYm3HjjjXjllVeQnp7eq3oYgMhfdDpkbK204LX/O4nyqiYAgEGnwXPfycS3Jg0VWxyRF7R3OrDzZAO2Ha7DtiO1OF5nu6JNXLgeGeZI3JAY4XwMicDIIRFIjDRAo2EwUrOACkAbN27EggULsG7dOuTk5GD16tXYtGkTDh8+jMTExCva79ixAzfffDOKiopwxx13YMOGDVi5ciXKysowYcIEAMDKlStRVFSEN954A2lpaXjyySexf/9+HDhwAEaj8bo1MQCRKLKs4OQFGyqqmlBR3YR/HazB+eZ2AIBeq8E3s5LxwC0jcUMiV86QOpxpbMOukw3YW+38f+LAeSs6HT3/2grRSkiKDkVyjBEpMWFIjDIgIcKAhAg9EiIMiI/QIzo0BNGhIQgN0bIXKQgFVADKycnBtGnTsGbNGgCALMtITU3Fww8/jMcff/yK9nPnzoXNZsOWLVvc12bMmIHs7GysW7cOiqIgOTkZP/vZz/Dzn/8cANDc3AyTyYTXX38d99xzz3Vr8lYAsrZ3wspt4YPO5f8Huf6uQIFDViArChyyc+6O3SGjvdPR/ZBh6+hCfWsHals6UNfSgRprOw6et8La3uXx/gkRenxvxnDMzxmOIZGc60Pq1tHlwKHzLThS04Jjda04XtuKY7WtqG68CIfc+19nIVoJUcYQRBh1CNfrEG7QItygQ5heC6NOC0OIFsYQDYwhWhh0GoRoNdBrNdB3/12nlaDTSNBpNQjRSNBonF9rNBK0kgStRoLG/ScgSc4/NZIEqftPAJAkQIIEjcb5p/NruJ/DFdckXB7bXBnO8+ql633ly0wYaQgZ9Pldffn9LXS7WLvdjj179mDp0qXuaxqNBnl5eSgtLe3xNaWlpSgsLPS4lp+fj82bNwMATp48CYvFgry8PPfz0dHRyMnJQWlpaY8BqKOjAx0dlybgWa3WgXysq/rTZ6fx/NbDXnlvCh4GnQaZKdHITo3BlOGx+NrYRBh0XN5OBAAGnRZZqTHISo3xuN7lkFHT0oGzjRdxrukizjZdRF1LB+paO3ChtQP1rXY02OxovtgJh6yg06Hggs2OCza7mA9C+PEto/Do7Axh319oAKqvr4fD4YDJZPK4bjKZcOjQoR5fY7FYemxvsVjcz7uuXa3NlxUVFeEXv/hFvz5DX+g0Egzc5yIofPlfSa5/fV3+LzWNxvmvP63k/FehXqtx/4vSGKJFmF6LIREGDIl0PhIiDBg1JAIZSZGc2EzURzqtBikxoUiJCb1mO0VR0GZ3oPliJ5ovdsLW0QWb3YE215/2LncvbXunAxc7Heh0yLB3yeh0KLB3ybA75O4QJaPLoaBLdn7tUACHLHf3+sqQFUBWFCgK3D3CiuKsQcGl55wdVwpcHViu551/h/scNcX9H+Dyv355IKenfrAvj/UoPbbqu4GMIekEz9figUEAli5d6tGrZLVakZqaOujf576bR+G+m0cN+vsSEVHvSJKEcIMO4QYdkq8Tlii4Cf1nZkJCArRaLWpqajyu19TUwGw29/gas9l8zfauP/vyngaDAVFRUR4PIiIiCl5CA5Ber8eUKVNQUlLivibLMkpKSpCbm9vja3Jzcz3aA0BxcbG7fVpaGsxms0cbq9WKnTt3XvU9iYiISF2ED4EVFhZi4cKFmDp1KqZPn47Vq1fDZrNh0aJFAIAFCxYgJSUFRUVFAIDFixdj1qxZWLVqFebMmYO33noLu3fvxvr16wE4uzcfeeQRPPPMM0hPT3cvg09OTkZBQYGoj0lERER+RHgAmjt3Lurq6rB8+XJYLBZkZ2dj69at7knMVVVV0GgudVTNnDkTGzZswBNPPIFly5YhPT0dmzdvdu8BBACPPvoobDYb7rvvPjQ1NeGmm27C1q1be7UHEBEREQU/4fsA+SNuhEhERBR4+vL7m2ttiYiISHUYgIiIiEh1GICIiIhIdRiAiIiISHUYgIiIiEh1GICIiIhIdRiAiIiISHUYgIiIiEh1GICIiIhIdYQfheGPXJtjW61WwZUQERFRb7l+b/fmkAsGoB60tLQAAFJTUwVXQkRERH3V0tKC6Ojoa7bhWWA9kGUZ586dQ2RkJCRJGtT3tlqtSE1NRXV1Nc8Z8yLeZ9/gffYN3mff4H32DW/eZ0VR0NLSguTkZI+D1HvCHqAeaDQaDB061KvfIyoqiv+D+QDvs2/wPvsG77Nv8D77hrfu8/V6flw4CZqIiIhUhwGIiIiIVIcByMcMBgNWrFgBg8EgupSgxvvsG7zPvsH77Bu8z77hL/eZk6CJiIhIddgDRERERKrDAERERESqwwBEREREqsMARERERKrDAORDa9euxYgRI2A0GpGTk4Ndu3aJLimgFRUVYdq0aYiMjERiYiIKCgpw+PBhjzbt7e148MEHER8fj4iICHznO99BTU2NoIqDw3PPPQdJkvDII4+4r/E+D46zZ8/ie9/7HuLj4xEaGorMzEzs3r3b/byiKFi+fDmSkpIQGhqKvLw8HD16VGDFgcfhcODJJ59EWloaQkNDMWrUKDz99NMeZ0fxPvfPJ598gm9+85tITk6GJEnYvHmzx/O9ua8NDQ2YP38+oqKiEBMTg3vvvRetra1eqZcByEc2btyIwsJCrFixAmVlZcjKykJ+fj5qa2tFlxawtm3bhgcffBCfffYZiouL0dnZidtuuw02m83d5qc//Snee+89bNq0Cdu2bcO5c+fw7W9/W2DVge3zzz/H//zP/2DixIke13mfB66xsRE33ngjQkJC8I9//AMHDhzAqlWrEBsb627z/PPP4ze/+Q3WrVuHnTt3Ijw8HPn5+WhvbxdYeWBZuXIlXnnlFaxZswYHDx7EypUr8fzzz+Oll15yt+F97h+bzYasrCysXbu2x+d7c1/nz5+PL774AsXFxdiyZQs++eQT3Hfffd4pWCGfmD59uvLggw+6v3Y4HEpycrJSVFQksKrgUltbqwBQtm3bpiiKojQ1NSkhISHKpk2b3G0OHjyoAFBKS0tFlRmwWlpalPT0dKW4uFiZNWuWsnjxYkVReJ8Hy2OPPabcdNNNV31elmXFbDYrv/rVr9zXmpqaFIPBoPz5z3/2RYlBYc6cOcoPf/hDj2vf/va3lfnz5yuKwvs8WAAo77zzjvvr3tzXAwcOKACUzz//3N3mH//4hyJJknL27NlBr5E9QD5gt9uxZ88e5OXlua9pNBrk5eWhtLRUYGXBpbm5GQAQFxcHANizZw86Ozs97ntGRgaGDRvG+94PDz74IObMmeNxPwHe58Hy7rvvYurUqbj77ruRmJiISZMm4dVXX3U/f/LkSVgsFo/7HB0djZycHN7nPpg5cyZKSkpw5MgRAMDevXuxfft2fOMb3wDA++wtvbmvpaWliImJwdSpU91t8vLyoNFosHPnzkGviYeh+kB9fT0cDgdMJpPHdZPJhEOHDgmqKrjIsoxHHnkEN954IyZMmAAAsFgs0Ov1iImJ8WhrMplgsVgEVBm43nrrLZSVleHzzz+/4jne58Fx4sQJvPLKKygsLMSyZcvw+eef4yc/+Qn0ej0WLlzovpc9/Rzhfe69xx9/HFarFRkZGdBqtXA4HHj22Wcxf/58AOB99pLe3FeLxYLExESP53U6HeLi4rxy7xmAKCg8+OCDqKysxPbt20WXEnSqq6uxePFiFBcXw2g0ii4naMmyjKlTp+KXv/wlAGDSpEmorKzEunXrsHDhQsHVBY+//OUvePPNN7FhwwaMHz8eFRUVeOSRR5CcnMz7rDIcAvOBhIQEaLXaK1bF1NTUwGw2C6oqeDz00EPYsmULPvroIwwdOtR93Ww2w263o6mpyaM973vf7NmzB7W1tZg8eTJ0Oh10Oh22bduG3/zmN9DpdDCZTLzPgyApKQnjxo3zuDZ27FhUVVUBgPte8ufIwCxZsgSPP/447rnnHmRmZuL73/8+fvrTn6KoqAgA77O39Oa+ms3mKxYGdXV1oaGhwSv3ngHIB/R6PaZMmYKSkhL3NVmWUVJSgtzcXIGVBTZFUfDQQw/hnXfewYcffoi0tDSP56dMmYKQkBCP+3748GFUVVXxvvfBrbfeiv3796OiosL9mDp1KubPn+/+O+/zwN14441XbONw5MgRDB8+HACQlpYGs9nscZ+tVit27tzJ+9wHbW1t0Gg8f/VptVrIsgyA99lbenNfc3Nz0dTUhD179rjbfPjhh5BlGTk5OYNf1KBPq6YevfXWW4rBYFBef/115cCBA8p9992nxMTEKBaLRXRpAeuBBx5QoqOjlY8//lg5f/68+9HW1uZuc//99yvDhg1TPvzwQ2X37t1Kbm6ukpubK7Dq4HD5KjBF4X0eDLt27VJ0Op3y7LPPKkePHlXefPNNJSwsTPnTn/7kbvPcc88pMTExyt///ndl3759yl133aWkpaUpFy9eFFh5YFm4cKGSkpKibNmyRTl58qTy9ttvKwkJCcqjjz7qbsP73D8tLS1KeXm5Ul5ergBQXnjhBaW8vFw5ffq0oii9u6+zZ89WJk2apOzcuVPZvn27kp6ersybN88r9TIA+dBLL72kDBs2TNHr9cr06dOVzz77THRJAQ1Aj4/f//737jYXL15UfvzjHyuxsbFKWFiY8q1vfUs5f/68uKKDxJcDEO/z4HjvvfeUCRMmKAaDQcnIyFDWr1/v8bwsy8qTTz6pmEwmxWAwKLfeeqty+PBhQdUGJqvVqixevFgZNmyYYjQalZEjRyr/+Z//qXR0dLjb8D73z0cffdTjz+SFCxcqitK7+3rhwgVl3rx5SkREhBIVFaUsWrRIaWlp8Uq9kqJctv0lERERkQpwDhARERGpDgMQERERqQ4DEBEREakOAxARERGpDgMQERERqQ4DEBEREakOAxARERGpDgMQERERqQ4DEBEREakOAxARERGpDgMQERERqY5OdAFERL5wyy23YOLEiTAajfjtb38LvV6P+++/H0899ZTo0ohIAPYAEZFqvPHGGwgPD8fOnTvx/PPP47/+679QXFwsuiwiEoCnwRORKtxyyy1wOBz49NNP3demT5+Or33ta3juuecEVkZEIrAHiIhUY+LEiR5fJyUloba2VlA1RCQSAxARqUZISIjH15IkQZZlQdUQkUgMQERERKQ6DEBERESkOgxAREREpDpcBUZERESqwx4gIiIiUh0GICIiIlIdBiAiIiJSHQYgIiIiUh0GICIiIlIdBiAiIiJSHQYgIiIiUh0GICIiIlIdBiAiIiJSHQYgIiIiUh0GICIiIlKd/w93kUvSshxgYgAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(ks, pmf_n)\n",
+    "plt.xlabel('n')\n",
+    "plt.ylabel('PMF');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8f05f187",
+   "metadata": {},
+   "source": [
+    "This is a negative binomial distribution, which we can confirm using `nbinom`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "b6e9de88",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.9999999094998685"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from scipy.stats import nbinom\n",
+    "\n",
+    "pmf_nbinom = nbinom.pmf(ns-k, k, p/100)\n",
+    "pmf_nbinom.sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "df8e2c43",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.allclose(pmf_n, pmf_nbinom)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b35dd161",
+   "metadata": {},
+   "source": [
+    "To see why this works we can compare the binomial PMF, which is a distribution over $k$ with $n$ and $p$ as parameters:\n",
+    "\n",
+    "$$Pr(k; n, p) = \\binom{n}{k} p^{k} (1-p)^{n-k}$$\n",
+    "\n",
+    "And the negative binomial PMF, which I've written as a distribution over $n$ with $k$ and $p$ as parameters:\n",
+    "\n",
+    "$$Pr(n; k, p) = \\binom{n-1}{k-1} p^k (1-p)^{n-k}$$\n",
+    "\n",
+    "This is not the most common way to parameterize the negative binomial distribution, but it shows that the only difference is in the binomial coefficient, because we know that the last trial is a success."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "aa0033cf",
+   "metadata": {},
+   "source": [
+    "## The beta distribution\n",
+    "\n",
+    "Suppose we have 101 devices that perform Bernoulli trials with different probabilities.\n",
+    "One device has $p=0$, one has $p=0.01$, one has $p=0.02$, and so on up to one device with $p=1$.\n",
+    "\n",
+    "Now suppose we choose one of the devices so that all values of $p$ are equally likely.\n",
+    "If we run $n$ trials and see $k$ successes, what is the distribution of $p$?\n",
+    "\n",
+    "We can answer this question by selecting a vector from the cube along the $p$ axis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 65,
+   "id": "a0ff5426",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "k = 25\n",
+    "n = 50\n",
+    "k = 6\n",
+    "n = 12\n",
+    "pdf_p = cube[k, n, :].copy()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e001ecbe",
+   "metadata": {},
+   "source": [
+    "The result is not normalized."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 66,
+   "id": "59f2ee8f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "7.692307692307231"
+      ]
+     },
+     "execution_count": 66,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pdf_p.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3f2195eb",
+   "metadata": {},
+   "source": [
+    "But we can normalize it like this."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 67,
+   "id": "90da0a89",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pdf_p /= pdf_p.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5d85a72c",
+   "metadata": {},
+   "source": [
+    "And here's what it looks like."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 68,
+   "id": "7e58f780",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAGwCAYAAABSN5pGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABX8klEQVR4nO3de1xT9/0/8FcSSML9KldRvKB4R0UQq1UrFatry2ZXay86f07brTonW+ulVru2q203N3uxdXbrdfrVulrXWkeLaG+KN/Cu4F3wEi4ihIskJDm/P0JiQVBAkk8ur+fjkYf18Enyzqkkr3xuRyZJkgQiIiIispKLLoCIiIjI0TAgERERETXBgERERETUBAMSERERURMMSERERERNMCARERERNcGARERERNSEh+gCnJXJZMLly5fh5+cHmUwmuhwiIiJqBUmSUFVVhaioKMjlLfcTMSC10+XLlxETEyO6DCIiImqHoqIidO7cucWfMyC1k5+fHwDzCfb39xdcDREREbWGVqtFTEyM9XO8JQxI7WQZVvP392dAIiIicjK3mx7DSdpERERETTAgERERETXBgERERETUhEMEpFWrViE2NhZqtRrJycnYu3fvLdtv3LgR8fHxUKvVGDBgALZu3dro5y+88ALi4+Ph4+ODoKAgpKamYs+ePY3alJeX47HHHoO/vz8CAwMxc+ZMVFdXd/hrIyIiIucjPCBt2LABGRkZWLZsGfLy8jBo0CCkpaWhpKSk2fa7du3C1KlTMXPmTBw4cADp6elIT0/H0aNHrW169eqFt99+G0eOHMGPP/6I2NhYjB8/HqWlpdY2jz32GI4dO4asrCxs2bIF33//PWbPnm3z10tERESOTyZJkiSygOTkZAwbNgxvv/02APMGjDExMZg7dy4WLlx4U/spU6agpqYGW7ZssR4bPnw4EhISsHr16mafQ6vVIiAgANu2bcO4ceNw4sQJ9O3bF/v27UNiYiIAIDMzExMnTsTFixcRFRV102PodDrodLpGjxkTE4PKykquYiMiInISlkxwu89voT1Ier0eubm5SE1NtR6Ty+VITU1FTk5Os/fJyclp1B4A0tLSWmyv1+uxZs0aBAQEYNCgQdbHCAwMtIYjAEhNTYVcLr9pKM5i+fLlCAgIsN64SSQREZHrEhqQysrKYDQaER4e3uh4eHg4NBpNs/fRaDStar9lyxb4+vpCrVbj73//O7KyshAaGmp9jLCwsEbtPTw8EBwc3OLzLlq0CJWVldZbUVFRm14rEREROQ+X3Shy7NixOHjwIMrKyvDee+/h4Ycfxp49e24KRq2lUqmgUqk6uEoiIiJyREJ7kEJDQ6FQKFBcXNzoeHFxMSIiIpq9T0RERKva+/j4oGfPnhg+fDj+9a9/wcPDA//617+sj9F0ErjBYEB5eXmLz0tERETuQ2hAUiqVGDp0KLKzs63HTCYTsrOzkZKS0ux9UlJSGrUHgKysrBbb//RxLZOsU1JSUFFRgdzcXOvPt2/fDpPJhOTk5Pa+HCIiInIRwofYMjIyMH36dCQmJiIpKQkrV65ETU0NZsyYAQCYNm0aoqOjsXz5cgDAvHnzMHr0aKxYsQKTJk3C+vXrsX//fqxZswYAUFNTgz//+c944IEHEBkZibKyMqxatQqXLl3CL3/5SwBAnz59MGHCBMyaNQurV69GfX095syZg0ceeaTZFWxERETkXoQHpClTpqC0tBRLly6FRqNBQkICMjMzrROxCwsLIZff6OgaMWIE1q1bhyVLlmDx4sWIi4vD5s2b0b9/fwCAQqFAfn4+PvroI5SVlSEkJATDhg3DDz/8gH79+lkfZ+3atZgzZw7GjRsHuVyOyZMn480337Tviycih1RRq4eXUgGVh0J0KUQkiPB9kJxVa/dRICLHVlatw7cFpci/okW+pgr5Gi3KqvVQyGXo0ckHvSP8ER/hh0GdAzGiRwjk8ltfAZyIHFtrP7+F9yAREYmgravHe9+fxT9/OIfr9cabfm40SThZXI2TxdX48pD5WP9ofzybFo9RcaGQyRiUiFwZAxIRuZW6eiP+vfsCVu04jWu19QCAvpH+SO4ejD4R/ugd4Ye4cF9UXq839yhdMfcqZZ8owdFLWkx7fy9SuodgwX3xSIgJFPtiiMhmOMTWThxiI3I+Ry9V4slPcnGp4joAoEcnHzyTFo+0fuG37RG6Wq3DO9+ewSc5F6A3mgAAjwyLwYsP9ofSQ/hlLYmolVr7+c2A1E4MSETOZUdBCZ5em4davRER/mrMvzcOk4d0hoeibeHm4rVarNx2Cp/lXYQkASN7huLdx4fAT+1po8qJqCMxINkYAxKR89iwrxCLPz8Ko0nqsEDzbUEJftsQuOIj/PDhjCREBKg7qGIishWnuFgtEZEtSZKEv2edxILPjsBokvCLwdF4/1fDOqS3Z0zvMGyYnYJQXxXyNVX4+Ts7UaCp6oCqicgRMCARkUuSJAl/+vI43sg+BQCYM7YnVjw8qEPnCw3oHIDPfzsC3Tv54EplHR5avQtHL1V22OMTkTgMSETkkv69+wI+3HUeMhnw55/3xx/TettkaX5MsDc+e2oEhnYNQlWdAbM/3o/SKl2HPw8R2RcDEhG5nF1nyvDCl8cBAAsmxOOx5K42fb4gHyXe/9UwdA/1weXKOvzm37nQG0w2fU4isi0GJCJyKYVXa/H02jwYTRLSE6Lw5N3d7fK8AV6eeG96IvzUHth/4RqW/vcouAaGyHkxIBGRy6jWGTDr4/24VluPgZ0D8OrkgXbd8bpHJ1+8OXUw5DJg/b4ifJxzwW7PTUQdiwGJiFyCySQhY8NBFBRXoZOfCmueSITa0/4Xmx3bOwwL74sHALy45Th2ni6zew1EdOcYkIjIJXyy+wK+OV4MpUKOfzwxVOieRLNGdccvBkfDaJIwb/1BVNTqhdVCRO3DgERETu/itVq8lpkPAHhuUh8M6RIktB6ZTIZXfjEAcWG+KKvW4eWvTgith4jajgGJiJyaJElY/PlR1OqNGBYbhCeG23bFWmupPRUNc6CA/+RexPcnS0WXRERtwIBERE5tU94lfH+yFEoPOV6dPBByuf0mZd/O0K5B+NWIWADAok1HUKMziC2IiFqNAYmInFZplQ4vfWXe7+j3qXHo0clXcEU3++P43ogO9MKliuv46zcFossholZiQCIip/XCl8dQUVuPflH+mDXKPvsdtZWPygPLfzEAAPDhrvPIvXBNcEVE1BoMSETklL45psFXh69AIZfhtckD4alw3Lezu3t1wuQhnSFJwILPDkNnMIouiYhuw3HfUYiIWlBXb8SfGi4lMvvu7ugfHSC4ott7/md9EOqrxOmSanyw87zocojoNhiQiMjpfJJzAZcqriMyQI154+JEl9Mqgd5KLLyvDwDgnR2nuTcSkYNjQCIip1J5vR5v7zgNAJh/by8hu2W3188HRyM+wg/aOgPe/faM6HKI6BYYkIjIqbz77RlUXq9Hr3BfTB7SWXQ5baKQy7BggvkyJB/sOo9LFdcFV0RELWFAIiKncaXyOj7YeQ4AsGBCPBQOtOdRa43p3QnDuwdDbzDh71knRZdDRC1gQCIip7Ey6xR0BhOSYoNxT3yY6HLaRSaTWecifZZ3EfkareCKiKg5DEhE5BROFVdhY24RAGDhxHjIZM7Xe2SREBOISQMiIUnA65ncPJLIETEgEZFTeC2zACYJmNAvQvjFaDvCH9N6QyGXYXt+CXafvSq6HCJqggGJiBxe7oVr2HaiGAq5DM9M6C26nA7RLdQHU5NiAACvZeZDkiTBFRHRTzEgEZHDe6dhWf9DQzo75PXW2ut34+Kg8pDjQGEFdp8tF10OEf0EAxIRObQCTRWy80sgkwFPjekhupwOFeanxsOJ5l6kd7/jvkhEjoQBiYgc2j8agsPE/pHoFuojuJqON/vu7lDIZfj+ZCmOXa4UXQ4RNWBAIiKHdfFaLf576DIA4KnRrtV7ZBET7I2fDYwEAKz+7qzgaojIggGJiBzWP384B6NJwsieoRjQ2fEvSNtelvD31eHLuHC1RnA1RAQwIBGRg7parcP6fYUAgN+42NyjpvpE+mNs704wScCa79mLROQIGJCIyCF9tOs86upNGNg5ACN6hIgux+Z+M6YnAGBj7kWUVNUJroaIGJCIyOFU6wz4KOcCAOA3o3s49a7ZrTUsNghDuwZBbzDhg53nRZdD5PYYkIjI4azfW4jK6/XoHuqD8f0iRJdjFzKZDL9pmIv075wL0NbVC66IyL0xIBGRQzEYTXj/x3MAbiyBdxf3xIehV7gvqnQGbNhbJLocIrfGgEREDmXbiRJcrqxDsI8S6YOjRZdjV3K5DDPu6gYA+PeeCzCZePkRIlEYkIjIoXyy+zwA4OHEGKg9FWKLEeDBhCj4qT1w4WotvjtVKrocIrfFgEREDuN0STV2nr4KuQx4LLmL6HKE8FZ64JdDzZcf+aRhojoR2R8DEhE5jH/vNgeCe+LDERPsLbgacZ5I6QoA2FFQgqLyWsHVELknBiQicgg1OgM+y70IAJjWEBDcVbdQH4yKC4Uk3QiNRGRfDEhE5BA+P3AJVToDuoX6YGTPUNHlCDctJRYAsGF/EerqjWKLIXJDDEhEJJwkSdb5No8P7wq5Gy3tb8k98WGIDvRCRW09vmy4YC8R2Q8DEhEJt/dcOQqKq+DlqcBDQzuLLschKOQyPDbcPFH9Ew6zEdmdQwSkVatWITY2Fmq1GsnJydi7d+8t22/cuBHx8fFQq9UYMGAAtm7dav1ZfX09FixYgAEDBsDHxwdRUVGYNm0aLl9u/A0sNjYWMpms0e3VV1+1yesjolv7uCEApA+OQoCXp+BqHMeUxBgoFXIcvliJg0UVosshcivCA9KGDRuQkZGBZcuWIS8vD4MGDUJaWhpKSkqabb9r1y5MnToVM2fOxIEDB5Ceno709HQcPXoUAFBbW4u8vDw8//zzyMvLw6ZNm1BQUIAHHnjgpsd68cUXceXKFett7ty5Nn2tRHSzEm0dvj6qAQA8MTxWbDEOJsRXhZ8NjAQAfJxzXmwxRG5GJkmS0K1ak5OTMWzYMLz99tsAAJPJhJiYGMydOxcLFy68qf2UKVNQU1ODLVu2WI8NHz4cCQkJWL16dbPPsW/fPiQlJeHChQvo0sXcZR0bG4vf//73+P3vf9+qOnU6HXQ6nfXvWq0WMTExqKyshL+/f2tfLhE1sWrHafzl6wIM7RqEz34zQnQ5DudA4TX8/J1dUHnIsfe5VPawEd0hrVaLgICA235+C+1B0uv1yM3NRWpqqvWYXC5HamoqcnJymr1PTk5Oo/YAkJaW1mJ7AKisrIRMJkNgYGCj46+++ipCQkIwePBg/OUvf4HBYGjxMZYvX46AgADrLSYmphWvkIhuRZIkbNxvvubYI8P4O9WchJhA9A73g85g4mRtIjsSGpDKyspgNBoRHh7e6Hh4eDg0Gk2z99FoNG1qX1dXhwULFmDq1KmNkuLvfvc7rF+/Hjt27MCTTz6JV155Bc8++2yLtS5atAiVlZXWW1ERLyRJdKf2X7iG81dr4a1UYOKASNHlOCSZTIZfJponrm9s2CeKiGzPQ3QBtlRfX4+HH34YkiTh3XffbfSzjIwM638PHDgQSqUSTz75JJYvXw6VSnXTY6lUqmaPE1H7fbrP/EXjZwMj4aNy6bejO5I+OBqv/i8fh4oqcLK4Cr3C/USXROTyhPYghYaGQqFQoLi4uNHx4uJiRERENHufiIiIVrW3hKMLFy4gKyvrtvOEkpOTYTAYcP78+ba/ECJqsxqdAV8duQLAfGFaalmorwrj+oQBgHVIkohsS2hAUiqVGDp0KLKzs63HTCYTsrOzkZKS0ux9UlJSGrUHgKysrEbtLeHo1KlT2LZtG0JCQm5by8GDByGXyxEWFtbOV0NEbfHVkSuo1RvRPdQHQ7sGiS7H4VkuYPv5gUuoN5oEV0Pk+oT3aWdkZGD69OlITExEUlISVq5ciZqaGsyYMQMAMG3aNERHR2P58uUAgHnz5mH06NFYsWIFJk2ahPXr12P//v1Ys2YNAHM4euihh5CXl4ctW7bAaDRa5ycFBwdDqVQiJycHe/bswdixY+Hn54ecnBzMnz8fjz/+OIKC+EZNZA//2W+eT/NQYmfIZNw5+3bG9O6ETn4qlFbpsCO/BOP7Nd/LTkQdQ3hAmjJlCkpLS7F06VJoNBokJCQgMzPTOhG7sLAQcvmNjq4RI0Zg3bp1WLJkCRYvXoy4uDhs3rwZ/fv3BwBcunQJX3zxBQAgISGh0XPt2LEDY8aMgUqlwvr16/HCCy9Ap9OhW7dumD9/fqN5SURkO2dLq7H3fDnkMmDyEO6c3RoeCjl+MTga//j+LD7df5EBicjGhO+D5Kxau48CEd3s9cx8vPPtGYzt3QkfzEgSXY7TOF1ShdS/fQ+FXIacRfcgzE8tuiQip+MU+yARkfsxmiR8lmceXvslJ2e3Sc8wPwzuEgijScLmA5dEl0Pk0hiQiMiuvj9VimKtDkHentaVWdR6lhV/n+6/CA4AENkOAxIR2ZVlcnb64GioPBSCq3E+PxsYCbWnHKdLqnkBWyIbYkAiIrvR1tUj64R5H7OHhnJydnv4qT0xsb951/HPOcxGZDMMSERkN5lHNNAbTOgV7ou+kVzc0F4PDo4GAGw5fIV7IhHZCAMSEdnNfw+ZezweTIjm3kd34K4eIQj1VaK8Ro8fT5eJLofIJTEgEZFdFGvrsOvMVQDAA4OiBFfj3DwUcvxsoPkc/pfDbEQ2wYBERHbx5aHLkCQgsWsQYoK9RZfj9B5MMAekb44Xo1ZvEFwNkethQCIiu/jvwcsAbnyw051JiAlE1xBv1OqNyDpefPs7EFGbMCARkc2dKa3GkUuV8JDLMGkgA1JHkMlkeLBhqNISPomo4zAgEZHNWT7AR8WFIthHKbga1/FAgnk12/cnS1FeoxdcDZFrYUAiIpuSJAlfHDRPJE5vWJ5OHaNnmC/6R/vDYJLw1ZErosshcikMSERkU4cuVuL81Vp4eSqQ2idcdDku58FB5tDJ1WxEHYsBiYhsynJR1fH9wuGj8hBcjeu5f1AUZDJg/4VrKCqvFV0OkctgQCIimzEYTdhy2Dz0w9VrthERoMbwbiEAgC8OcbI2UUdhQCIim8k5exVl1ToE+ygxKq6T6HJcVvpgy2o2DrMRdRQGJCKymS2HzL1H9/WPgKeCbze2MqF/JDwVMpwsrsap4irR5RC5BL5jEZFN1BtN+Pq4BgAwaWCk4GpcW4CXp7WHjqvZiDoGAxIR2cSuM1dRUVuPUF8lkhvmyJDtTBxgDqFbGZCIOgQDEhHZxNaGydkT+kdAIZcJrsb13ds3nMNsRB2IAYmIOtxPh9csPRtkWxxmI+pYDEhE1OE4vCaGJYx+dZgBiehOMSARUYfj8JoYlmG2UyUcZiO6UwxIRNShOLwmDofZiDoOAxIRdSgOr4k1icNsRB2CAYmIOtRXh82Xu+DwmhipPxlmO8lhNqJ2Y0Aiog5TbzTh62PFAIBJA3jtNRECvDxxt2WYjb1IRO3GgEREHWbn6TJUXq9HqK8KSd2CRZfjtrhpJNGdY0Aiog5j+UCe0D+cw2sCcZiN6M4xIBFRh6g3mvDNcfPwGlevifXTYTb2IhG1DwMSEXWIvefKUVFbjxAfrl5zBBP6RwCAdU4YEbUNAxIRdYjMo+a9j1L7cHjNEVj+P5y4okXh1VrR5RA5HQYkIrpjJpOEbxo2h7T0XJBYQT5KJDdMlP/6mEZwNUTOhwGJiO7YwYsVKNbq4KvywIieHF5zFJawmsmARNRmDEhEdMe+bhheGxsfBpWHQnA1ZDG+rzkg5V64hhJtneBqiJwLAxIR3RFJkqxDOBP6cXjNkUQEqDG4SyAAWFcYElHrMCAR0R0pKK7C+au1UHrIMaZ3J9HlUBNp/Syr2TjMRtQWDEhEdEcsq9fujguFj8pDcDXUlCUg5Zy5iopaveBqiJwHAxIR3RFLQErj8JpD6hbqg/gIPxhMErJPlIguh8hpMCARUbtduFqDfE0VFHIZUvuEiy6HWsBhNqK2Y0AionazfOAO7x6MIB+l4GqoJZaA9N3JUtTqDYKrIXIODEhE1G4cXnMOfSL90CXYGzqDCd8VlIouh8gpMCARUbuUaOuQV1gB4MZ+O+SYZDLZT67NxmE2otZgQCKidrHsq5MQE4iIALXgauh2LL182SdKoDeYBFdD5PgYkIioXSwBaXw/Ts52BoNjAhHqq0KVzoA9566KLofI4TEgEVGbVdXVI+dMGQAOrzkLuVyG1D5hAIAs7qpNdFsOEZBWrVqF2NhYqNVqJCcnY+/evbdsv3HjRsTHx0OtVmPAgAHYunWr9Wf19fVYsGABBgwYAB8fH0RFRWHatGm4fPlyo8coLy/HY489Bn9/fwQGBmLmzJmorq62yesjcjXfnSxFvVFC91Af9AzzFV0OtdK9fc29fduOF0OSJMHVEDk24QFpw4YNyMjIwLJly5CXl4dBgwYhLS0NJSXNb2i2a9cuTJ06FTNnzsSBAweQnp6O9PR0HD16FABQW1uLvLw8PP/888jLy8OmTZtQUFCABx54oNHjPPbYYzh27BiysrKwZcsWfP/995g9e7bNXy+RK7D0QFg+cMk53NUzFN5KBS5X1uHYZa3ocogcmkwS/DUiOTkZw4YNw9tvvw0AMJlMiImJwdy5c7Fw4cKb2k+ZMgU1NTXYsmWL9djw4cORkJCA1atXN/sc+/btQ1JSEi5cuIAuXbrgxIkT6Nu3L/bt24fExEQAQGZmJiZOnIiLFy8iKirqtnVrtVoEBASgsrIS/v7+7XnpRE6p3mjC0JeyoK0z4D9PpSAxNlh0SdQGT32Si8xjGvxuXBwy7u0luhwiu2vt57fQHiS9Xo/c3FykpqZaj8nlcqSmpiInJ6fZ++Tk5DRqDwBpaWkttgeAyspKyGQyBAYGWh8jMDDQGo4AIDU1FXK5HHv27Gn2MXQ6HbRabaMbkTvae64c2joDQnyUGNwlSHQ51EaWXj/OQyK6NaEBqaysDEajEeHhjbvpw8PDodE0v1eHRqNpU/u6ujosWLAAU6dOtSZFjUaDsLCwRu08PDwQHBzc4uMsX74cAQEB1ltMTEyrXiORq7F8sI7rEwaFXCa4Gmqre+LN/99OXNGiqLxWdDlEDkv4HCRbqq+vx8MPPwxJkvDuu+/e0WMtWrQIlZWV1ltRUVEHVUnkPCRJ+sn8I65ec0ZBPkokdjX3/G07wV4kopYIDUihoaFQKBQoLm78S1pcXIyIiObffCMiIlrV3hKOLly4gKysrEbjjBERETdNAjcYDCgvL2/xeVUqFfz9/RvdiNzNsctaXKq4DrWnHCN7hoouh9rJMsz2zTEGJKKWCA1ISqUSQ4cORXZ2tvWYyWRCdnY2UlJSmr1PSkpKo/YAkJWV1ai9JRydOnUK27ZtQ0hIyE2PUVFRgdzcXOux7du3w2QyITk5uSNeGpFLsvQejYrrBC+lQnA11F6Wvav2ni9HRa1ecDVEjkn4EFtGRgbee+89fPTRRzhx4gR+85vfoKamBjNmzAAATJs2DYsWLbK2nzdvHjIzM7FixQrk5+fjhRdewP79+zFnzhwA5nD00EMPYf/+/Vi7di2MRiM0Gg00Gg30evMbQZ8+fTBhwgTMmjULe/fuxc6dOzFnzhw88sgjrVrBRuSuLAFpPJf3O7UuId7oHe4Ho0nCjoLmt1QhcnceoguYMmUKSktLsXTpUmg0GiQkJCAzM9M6EbuwsBBy+Y0cN2LECKxbtw5LlizB4sWLERcXh82bN6N///4AgEuXLuGLL74AACQkJDR6rh07dmDMmDEAgLVr12LOnDkYN24c5HI5Jk+ejDfffNP2L5jISV28VovjV7SQy4BxfRiQnN29fcNRUFyFrOPF+PngzqLLIXI4wvdBclbcB4nczYc7z+GFL48jKTYYnz7V/BA4OY9DRRV4cNVO+CgVyFt6L1QeHDIl9+AU+yARkfPIOsHds13JgOgAhPurUKM3YtcZXryWqCkGJCK6LW1dPfacLQcApDIguQTzxWtvXJuNiBpjQCKi2/quoBQGk4SeYb7oFuojuhzqIJawm32ihBevJWqCAYmIbsuyoeC4PmG3aUnOJKV7CLyVCmi0vHgtUVMMSER0S/VGE3bkm5eCp3L1mktReyqsG35yV22ixhiQiOiW9p+/Bm2dAUHenhjCi9O6HMswGwMSUWMMSER0S9kNH5xj43lxWld0T3wYZDLg6CUtrlReF10OkcNgQCKiFkmSZO1ZuJfDay4p1FeFwTGBAMyTtYnIjAGJiFp0prQG56/WQqmQY1SvTqLLIRux7IyezWE2IisGJCJqkaX3aHiPEPiqhF+ZiGzEsvnnzjNXUas3CK6GyDEwIBFRiyw9Cqlc3u/S4sJ8ERPsBb3BhB9OlYkuh8ghMCARUbPKa/TIvXANAC9O6+pkMu6qTdQUAxIRNWtHfglMEtAn0h/RgV6iyyEbswSkHQUlMJm4qzYRAxIRNSs7n8Nr7mRYbDD8VB4oq9bj4MUK0eUQCceAREQ30RmM+K6gFAB3z3YXSg85Rvc2r1TkMBsRAxIRNWPP2XLU6I3o5KfCgOgA0eWQnVjnIXG5PxEDEhHdzLJ6bVx8GOTcPdttjOndCQq5DCeLq1FUXiu6HCKhGJCIqBHz7tnmHZW5es29BHorMbSr+Xp73DSS3B0DEhE1crK4GpcqrkPlIbde6Z3ch2VSfnY+LztC7o0BiYgascw/GdEjBF5KheBqyN4svYZ7zpajWsddtcl9MSARUSPW+UccXnNL3UN9EBviDb3RhB9Oloouh0gYBiQisrparcOBogoAwDjuf+SWZDLZjYvXcpiN3BgDEhFZ7SgohSQBfSP9ERnA3bPdlSUc78gvgZG7apObYkAiIqvt3D2b0LCrttoDV2v0OMRdtclNMSAREQBAbzDh+5PmK7nfw/lHbs1TIcfoXuZdtbncn9wVAxIRAQD2njOvWurkp8JA7p7t9izDbNknOA+J3BMDEhEBuLG8/57e3D2bgDG9wiCXAfmaKly8xl21yf0wIBERJElCdsP8o3s4/4gABPkokdg1GACwnavZyA0xIBERTpdUo6j8OpTcPZt+4h4Os5EbY0AiIuu111K6h8BH5SG4GnIUltWMOWeuooa7apObYUAiIi7vp2b16OSLrpZdtU+ViS6HyK4YkIjc3LUaPXIvXAMAjI1nQKIbZDIZ7mn4N2EJ0UTuggGJyM3tKCiBSQLiI/zQOchbdDnkYFIb9sTanl8KE3fVJjfCgETk5izX20rl5pDUjGGxwfBTeaCsWofDlypFl0NkNwxIRG5MbzDh+wLzFdu5vJ+ao/SQ427uqk1uiAGJyI3tP1+OKp0BIT5KJHQOFF0OOSjuqk3uiAGJyI1ZlvePjefu2dSyMb3Nu2ofv6LF5YrrosshsgsGJCI39dPds7m8n24l2EeJIV2CANyYs0bk6hiQiNzUmdIaXLhaC6VCjpFxnUSXQw7OMkdtO+chkZtgQCJyU5YJt8ndg+HL3bPpNiyrHHeeuYpaPXfVJtfHgETkpri8n9oiLswXnYO8oDeYsPP0VdHlENkcAxKRG6qovbF79j3cPZtaQSaTWcM0l/uTO2BAInJD3xaUwmiS0DvcDzHB3D2bWscSprPzS7irNrk8BiQiN2QZXhvH1WvUBsndg+GjVKC0Soejl7mrNrm2NgWkpUuXora21vr3a9eudXhBRGRb9UYTvi1gQKK2U3koMKphxeM2bhpJLq5NAenPf/4zqqurrX/v2rUrzp492+FFEZHt7Dtfjqo6A4J9lEiICRJdDjmZG7tqcx4SubY2BSRJkm759/ZYtWoVYmNjoVarkZycjL17996y/caNGxEfHw+1Wo0BAwZg69atjX6+adMmjB8/HiEhIZDJZDh48OBNjzFmzBjIZLJGt6eeeuqOXwuRM7BcLmJM705QcPdsaqOx8WGQyYBjl7W4Usldtcl1CZ2DtGHDBmRkZGDZsmXIy8vDoEGDkJaWhpKS5rtud+3ahalTp2LmzJk4cOAA0tPTkZ6ejqNHj1rb1NTUYOTIkXjttddu+dyzZs3ClStXrLfXX3+9Q18bkSOSJMn6zf9eLu+ndgj1VWFwTCAAXpuNXFubApJMJkNVVRW0Wi0qKyshk8lQXV0NrVbb6NZaf/vb3zBr1izMmDEDffv2xerVq+Ht7Y3333+/2fZvvPEGJkyYgGeeeQZ9+vTBSy+9hCFDhuDtt9+2tnniiSewdOlSpKam3vK5vb29ERERYb35+/u3um4iZ3WmtAbnG3bPHtWLu2dT+4zjcn9yA20eYuvVqxeCgoIQHByM6upqDB48GEFBQQgKCkJgYCCCglo3p0Gv1yM3N7dRkJHL5UhNTUVOTk6z98nJybkp+KSlpbXY/lbWrl2L0NBQ9O/fH4sWLWo0+bw5Op2u3UGQyFFw92zqCNxVm9xBm94hd+zY0WFPXFZWBqPRiPDwxt384eHhyM/Pb/Y+Go2m2fYajaZNz/3oo4+ia9euiIqKwuHDh7FgwQIUFBRg06ZNLd5n+fLl+NOf/tSm5yFyNNtOWC5Oy+E1ar9e4eZdtS9eu44fTpUhrV+E6JKIOlybAtLo0aNtVYddzZ492/rfAwYMQGRkJMaNG4czZ86gR48ezd5n0aJFyMjIsP5dq9UiJibG5rUSdZRrNTd2z+byfroTll21P9x1HtknihmQyCUJm6QdGhoKhUKB4uLGY9jFxcWIiGj+ly0iIqJN7VsrOTkZAHD69OkW26hUKvj7+ze6ETmTHQUlMElAfIQfOgdx92y6M5ZeyO35pdxVm1xSmwKSQqFo1a01lEolhg4diuzsbOsxk8mE7OxspKSkNHuflJSURu0BICsrq8X2rWXZCiAyMvKOHofIkVlWHHF4jTpCUrdg+Kk8UFatw6GLFaLLIepwbRpikyQJXbt2xfTp0zF48OA7fvKMjAxMnz4diYmJSEpKwsqVK1FTU4MZM2YAAKZNm4bo6GgsX74cADBv3jyMHj0aK1aswKRJk7B+/Xrs378fa9assT5meXk5CgsLcfnyZQBAQUEBAFhXq505cwbr1q3DxIkTERISgsOHD2P+/Pm4++67MXDgwDt+TUSOSG8w4buTpQA4vEYdQ+khx929O+Grw1eQfaIEg7tw01FyMVIb7Nu3T3rqqaekwMBAafDgwdJbb70llZeXt+UhbvLWW29JXbp0kZRKpZSUlCTt3r3b+rPRo0dL06dPb9T+008/lXr16iUplUqpX79+0ldffdXo5x988IEE4KbbsmXLJEmSpMLCQunuu++WgoODJZVKJfXs2VN65plnpMrKyjbVXVlZKQFo8/2IRPj+ZInUdcEWaehLWZLRaBJdDrmITXlFUtcFW6S0v38nuhSiVmvt57dMktq+HXZdXR3+85//4IMPPsDu3btx//33Y+bMmbj33ns7NLw5Mq1Wi4CAAFRWVnI+Ejm8F744hg93ncfDiZ3x+kODRJdDLuJajR5DX86CSQJ+XDCWc9vIKbT287tdk7TVajUef/xxZGdn4+jRoygpKcGECRNQXl7e7oKJyDYkSeLyfrKJIB8lErsGA+Cu2uR62r2K7eLFi3j55Zdx7733Ij8/H8888wx7Uogc0Mnialy8dh1KDzlGxoWKLodcTGpf85y2bdxVm1xMmwKSXq/Hhg0bMH78eMTFxSEvLw8rV65EUVERXn31VXh4cGdeIkdj+eC6q0cIvJX8HaWOZbnsyO6zV1FVVy+4GqKO06Z3y8jISPj5+WH69Ol45513EBZm/uZQU1PTqB17kogchyUgjePwGtlAj06+6Bbqg3NlNfjhVBkmDuB2KeQa2tSDdO3aNRQWFuKll15C7969rddga8+12IjI9kqrdDhYVAGAy/vJdlIb/m1tO85hNnIdwq7FRkS2tz2/GJIEDIgOQGSAl+hyyEWl9gnHez+cw/aCEhiMJngohF2kgajDtCkgjRw5En/961/xxRdfQK/XY9y4cVi2bBm8vPjGS+SIshq+0d/bl8NrZDtDuwYhyNsT12rrsf/CNQzvHiK6JKI71qaY/8orr2Dx4sXw9fVFdHQ03njjDTz99NO2qo2I7sB1vRE/nCoDwOX9ZFseCjnGxpuH2bI4zEYuok0B6eOPP8Y777yDr7/+Gps3b8aXX36JtWvXwmQy2ao+ImqnH06VQmcwITrQC30i/USXQy5ufEMv5bYTxWjH/sNEDqdNAamwsBATJ060/j01NRUymcx63TMichw/HV6TyWSCqyFXNyquE5Qecly4WotTJdWiyyG6Y20KSAaDAWq1utExT09P1Ndz7wsiR2I0Sdieb97ZmPOPyB58VB64q4d57hGH2cgVtGmStiRJ+NWvfgWVSmU9VldXh6eeego+Pj7WY5s2beq4ComozQ4UXsPVGj381R5I6hYsuhxyE/f2jcCOglJkHS/G02N7ii6H6I60KSBNnz79pmOPP/54hxVDRB0jq2FzyLHxYfDkkmuyk3F9woDPgYNFFSjR1iHMX337OxE5qDYFpA8++MBWdRBRB7IMcXD1GtlTuL8ag2ICcaioAttOlODR5C6iSyJqN361JHIxZ0qrcba0Bp4KGcb07iS6HHIz9/bhxWvJNTAgEbkYy+UehncPgZ/aU3A15G7u7RsBAPjxdBlqdAbB1RC1HwMSkYvh7tkkUq9wX8QEe0FvMOGHU6WiyyFqNwYkIhdytVqH3MJrADj/iMSQyWS4t4+5FynreIngaojajwGJyIVk55dAkoB+Uf6ICuQ1EkmM1L7meUjb84thMPJKC+ScGJCIXMg3xzQAgPEN80CIREiKDUZgw8Vr952/JroconZhQCJyETU6A75vuDhtWn8Or5E4Hgo5xsWb/w1+c1wjuBqi9mFAInIR358shd5gQpdgb/QO58VpSazx/RoC0jFevJacEwMSkYv4umF4La0fL05L4t0d1wlqTzkuVVzHscta0eUQtRkDEpEL0BtMyG64OG1aP84/IvG8lAqM7mXeqNQS3omcCQMSkQvYffYqquoMCPVVYUiXINHlEAG4Eda/OcZdtcn5MCARuQDLN/R7+4ZDLufwGjmGcfHhUMhlKCiuwvmyGtHlELUJAxKRkzOZJOvu2ZaJsUSOIMDbE8O7BwPgMBs5HwYkIid38GIFSqp08FV5YESPENHlEDViGWZjQCJnw4BE5OQsHzxj48Og8lAIroaoMcumpXmFFSjR1gmuhqj1GJCInJgkSdYJsGkcXiMHFBGgxqCYQABA1glO1ibnwYBE5MROlVTjXFkNlAo5xvQOE10OUbMs4f1rrmYjJ8KAROTEvj5qHl67q2cIfFUegqshap5lmC3nTBm0dfWCqyFqHQYkIif29XHL7tncHJIcV88wX/To5IN6o4QdDRuaEjk6BiQiJ1VUXoujl7SQy4DUvpx/RI5tQn9ziM88ytVs5BwYkIiclOWDJqlbMEJ9VYKrIbq1+/pHAgB2FJSgVm8QXA3R7TEgETmprUevAAAmDogUXAnR7fWL8kdMsBfq6k34tqBUdDlEt8WAROSErlRex4HCCshknH9EzkEmk2FiQy/S1iNXBFdDdHsMSEROyDK8ltg1COH+asHVELXOfQ29ndvzS1BXbxRcDdGtMSAROaH/HTEHJMu8DiJnMKhzAKIC1KjVG/H9SQ6zkWNjQCJyMiVVddh3oRzAjZVBRM5AJpNhQkOo/x9Xs5GDY0AicjJfHyuGJAEJMYGICvQSXQ5Rm0wcYA71244XQ2fgMBs5LgYkIifzvyOW1WvsPSLnM6RLEML8VKjSGbDzdJnocohaxIBE5ESuVuuw++xVAJx/RM5JLpfhvoah4a1HOMxGjosBiciJfHO8GCYJGBAdgJhgb9HlELWLZTVb1vFi1BtNgqshah4DEpETsewfcx+H18iJDYsNRqivEpXX65Fz5qrocoiaxYBE5CQqavXWDxMOr5EzU8hlGN+wwen/jnLTSHJMDEhETuKb48UwmCTER/ihW6iP6HKI7ohlV+2vjxXDwGE2ckDCA9KqVasQGxsLtVqN5ORk7N2795btN27ciPj4eKjVagwYMABbt25t9PNNmzZh/PjxCAkJgUwmw8GDB296jLq6Ojz99NMICQmBr68vJk+ejOLi4o58WUQdbsthXnuNXMfw7sEI9lGivEaPXRxmIwckNCBt2LABGRkZWLZsGfLy8jBo0CCkpaWhpKSk2fa7du3C1KlTMXPmTBw4cADp6elIT0/H0aNHrW1qamowcuRIvPbaay0+7/z58/Hll19i48aN+O6773D58mX84he/6PDXR9RRymv01iXR9w+KElwN0Z3zUMitq9m2HL4suBqim8kkSZJEPXlycjKGDRuGt99+GwBgMpkQExODuXPnYuHChTe1nzJlCmpqarBlyxbrseHDhyMhIQGrV69u1Pb8+fPo1q0bDhw4gISEBOvxyspKdOrUCevWrcNDDz0EAMjPz0efPn2Qk5OD4cOHN1urTqeDTqez/l2r1SImJgaVlZXw9/dv9zkgao21ey7guc+Pon+0P7bMHSW6HKIOkXPmKqa+txv+ag/sX3IvlB7CBzXIDWi1WgQEBNz281vYv0a9Xo/c3FykpqbeKEYuR2pqKnJycpq9T05OTqP2AJCWltZi++bk5uaivr6+0ePEx8ejS5cut3yc5cuXIyAgwHqLiYlp9XMS3akth8zDa/cPZO8RuY6kbsEI81NBW2fAD6d4bTZyLMICUllZGYxGI8LDwxsdDw8Ph0bT/OZhGo2mTe1begylUonAwMA2Pc6iRYtQWVlpvRUVFbX6OYnuRIm2DrvPmedoTBrI+UfkOhRymXVO3ZeHOMxGjoX9ma2kUqng7+/f6EZkD18duQJJAoZ0CUTnIG4OSa7FMqcu63gx6up5bTZyHMICUmhoKBQKxU2rx4qLixER0fwmeBEREW1q39Jj6PV6VFRU3NHjENmLZfUaJ2eTKxrSJRDRgV6o0RuxI7/5BTpEIggLSEqlEkOHDkV2drb1mMlkQnZ2NlJSUpq9T0pKSqP2AJCVldVi++YMHToUnp6ejR6noKAAhYWFbXocInu4eK0WuReuQSbj8n5yTTKZDD9rGDr+kqvZyIF4iHzyjIwMTJ8+HYmJiUhKSsLKlStRU1ODGTNmAACmTZuG6OhoLF++HAAwb948jB49GitWrMCkSZOwfv167N+/H2vWrLE+Znl5OQoLC3H5svkXraCgAIC55ygiIgIBAQGYOXMmMjIyEBwcDH9/f8ydOxcpKSktrmAjEuWrht6j5G7BCPdXC66GyDbuHxSFf3x/FtknSlCtM8BXJfSjiQiA4IA0ZcoUlJaWYunSpdBoNEhISEBmZqZ1InZhYSHk8hudXCNGjMC6deuwZMkSLF68GHFxcdi8eTP69+9vbfPFF19YAxYAPPLIIwCAZcuW4YUXXgAA/P3vf4dcLsfkyZOh0+mQlpaGd955xw6vmKhtLMNrP+PqNXJh/aL8ERvijfNXa5F9ohgPJkSLLolI7D5Izqy1+ygQtde5shqM/eu3UMhl2Lt4HEJ8VaJLIrKZFd8U4K3tp5HaJwz/nD5MdDnkwhx+HyQiurUtDcue7+oZynBELs+yCOG7k6WorK0XXA0RAxKRw7oxvMbJ2eT6eoX7oVe4L+qNEr4+1vq97YhshQGJyAEdv6xFQXEVlAo50vpy+wlyDw809CJtPnhJcCVEDEhEDsnyATGuTxgCvD0FV0NkH5bJ2Tlnr+JK5XXB1ZC7Y0AicjBGk4T/NgSk9MFczUPuIybYG8NigyBJwH8Pck8kEosBicjB7DpThmKtDoHenhjbO0x0OUR29fPBnQEAmw9wmI3EYkAicjCfN3wwTBoQCaUHf0XJvUwaEAmlQo58TRWOX9aKLofcGN99iRxIrd6AzKPmFTy/GMLhNXI/Ad6euCfe3HPKydokEgMSkQPJOl6MWr0RXYK9MaRLkOhyiISwzL3778FLMJq4lzGJwYBE5EA25d2YnC2TyQRXQyTG2PhOCPDyRLFWh5wzV0WXQ26KAYnIQZRU1eGHU6UAgJ9z9Rq5MZWHwrpB6qYDFwVXQ+6KAYnIQXx56ApMEpAQE4huoT6iyyESyvIl4eujGtTqDYKrIXfEgETkID5v+KbMydlEwNCuQYgJ9kKN3ois48WiyyE3xIBE5ABOFVfh6CUtPOQy/GxglOhyiISTyWT4ecPO2p9zTyQSgAGJyAF81jA5e3SvTgj2UQquhsgxWFaz/XCqDMXaOsHVkLthQCISzGA04bM88/DaLxM7C66GyHF07+SLoV2DYDRJ1t8RInthQCIS7NuCUpRW6RDio8Q98eGiyyFyKFMSYwAAG/dfhCRxTySyHwYkIsE27C8CYF61w0uLEDU2cWAkvJUKnCurwb7z10SXQ26E78ZEApVU1WF7fgkAYMqwGMHVEDkeX5WHdU+kTxu+TBDZAwMSkUCf55kvpTC4SyDiwv1El0PkkCxfHr46fAVVdfWCqyF3wYBEJIgkSdbhtYcT2XtE1JIhXYLQvZMPrtcbseXwFdHlkJtgQCISJK/wGs6W1sDL88ZlFYjoZjKZzDpZm8NsZC8MSESCbNhnfqOfNDASfmpPwdUQObZfDOkMhVyGA4UVOFVcJboccgMMSEQCVOsM1qECDq8R3V4nPxXuiQ8DcOPLBZEtMSARCbD18BXU6o3oHuqDYbFBosshcgqWYbbPD1yC3mASXA25OgYkIgEsk7N/mRgDmUwmuBoi5zCmdyeE+alwtUaP7fm8gC3ZFgMSkZ3la7TIvXANCrkMk4dEiy6HyGl4KOSYPNR8OZ61ewoFV0OujgGJyM7+vfsCACCtXzjC/NWCqyFyLo8mdYFMZr6A7bmyGtHlkAtjQCKyo6q6enyedwkA8PjwroKrIXI+McHeGNvbPFnb8mWDyBYYkIjs6PMDl1CjN6JHJx+kdA8RXQ6RU3qi4cvFxv1FuK43Cq6GXBUDEpGdSJKET3LM33ifGN6Vk7OJ2unuXp0QE+wFbZ0BXx66LLocclEMSER2sudcOU6VVMNbqcAvGiaaElHbKeQyPJZs7kX6ePd5SJIkuCJyRQxIRHbyScN8ifTB0fDnztlEd+ThxBgoPeQ4ekmLQxcrRZdDLogBicgOSrR1+PqoBgDweDInZxPdqWAfpfUahpaha6KOxIBEZAf/t7cIBpOExK5B6BvlL7ocIpdgmaz95eHLKK/RC66GXA0DEpGN1RtNWLe3YXJ2CnuPiDpKQkwg+kf7Q28wYeN+Xp+NOhYDEpGNbTtejGKtDiE+SkzoHyG6HCKXIZPJrL1I/95zAUYTJ2tTx2FAIrKxf/14DgAwZVgMVB4KwdUQuZYHBkUjwMsTReXXkXWc12ejjsOARGRDBwqvYf+Fa/BUyDB9RKzocohcjpdSgceSuwAA/vnDWcHVkCthQCKyoX/+YO49emBQNMJ53TUim5g+IhaeChn2X7iGA4XXRJdDLoIBichGispr8b+jVwAAvx7VTXA1RK4r3F+NBwZFA7jxpYToTjEgEdnIv348B5MEjIoLRZ9ILu0nsqVZd5u/hPzv6BUUldcKroZcAQMSkQ1U1tbj04Zlx7NGdRdcDZHri4/wx6i4UJikGwsjiO4EAxKRDazbW4havRHxEX4YFRcquhwit2D5MvLp/iJU1tYLroacHQMSUQfTG0z4cJf5G+yvR3WHTCYTXBGRexgVF4r4CD/U6o1Yt7dQdDnk5BiQiDrYlsOXUazVIcxPhQcGRYkuh8htyGQy/LqhF+nDXeegN5gEV0TOzCEC0qpVqxAbGwu1Wo3k5GTs3bv3lu03btyI+Ph4qNVqDBgwAFu3bm30c0mSsHTpUkRGRsLLywupqak4depUozaxsbGQyWSNbq+++mqHvzZyL5IkYc335r1Ypo+IhdLDIX7FiNzGA4OiEOanQrFWhy8OXRZdDjkx4e/eGzZsQEZGBpYtW4a8vDwMGjQIaWlpKCkpabb9rl27MHXqVMycORMHDhxAeno60tPTcfToUWub119/HW+++SZWr16NPXv2wMfHB2lpaairq2v0WC+++CKuXLlivc2dO9emr5VcX/aJEuRrquD9k83riMh+lB5yzLjLvKLtnW9P8/Ij1G7CA9Lf/vY3zJo1CzNmzEDfvn2xevVqeHt74/3332+2/RtvvIEJEybgmWeeQZ8+ffDSSy9hyJAhePvttwGYv8GvXLkSS5YswYMPPoiBAwfi448/xuXLl7F58+ZGj+Xn54eIiAjrzcfHx9Yvl1yYJEl4I9vcUzktJRaB3krBFRG5p8eHd0GAlyfOltZgy2H2IlH7CA1Ier0eubm5SE1NtR6Ty+VITU1FTk5Os/fJyclp1B4A0tLSrO3PnTsHjUbTqE1AQACSk5NvesxXX30VISEhGDx4MP7yl7/AYDC0WKtOp4NWq210I/qpbwtKceRSJbw8FZjFjSGJhPFTe+LXI82/g29tPw0Te5GoHYQGpLKyMhiNRoSHhzc6Hh4eDo1G0+x9NBrNLdtb/rzdY/7ud7/D+vXrsWPHDjz55JN45ZVX8Oyzz7ZY6/LlyxEQEGC9xcTEtP6FksuTJAkrrb1HXRHiqxJcEZF7m35XLPzVHjhdUo2tDTvaE7WF8CE2UTIyMjBmzBgMHDgQTz31FFasWIG33noLOp2u2faLFi1CZWWl9VZUVGTnismRfXeyFIeKKqD2lFtX0RCROP5qT/y/hl6kN7NPsReJ2kxoQAoNDYVCoUBxcXGj48XFxYiIiGj2PhEREbdsb/mzLY8JAMnJyTAYDDh//nyzP1epVPD39290IwIazz16PLkrOvmx94jIEcy4qxv8VB44WVyNr481PypB1BKhAUmpVGLo0KHIzs62HjOZTMjOzkZKSkqz90lJSWnUHgCysrKs7bt164aIiIhGbbRaLfbs2dPiYwLAwYMHIZfLERYWdicvidzQj6fLcKCwAioPOWaPZu8RkaMI8PLEjLtiAQBvsBeJ2shDdAEZGRmYPn06EhMTkZSUhJUrV6KmpgYzZswAAEybNg3R0dFYvnw5AGDevHkYPXo0VqxYgUmTJmH9+vXYv38/1qxZA8C8Udjvf/97vPzyy4iLi0O3bt3w/PPPIyoqCunp6QDME7337NmDsWPHws/PDzk5OZg/fz4ef/xxBAUFCTkP5JwkScIb28y9R48md0GYn1pwRUT0U/9vZDe8v/M88jVV+OZ4MSb0b3kkgeinhAekKVOmoLS0FEuXLoVGo0FCQgIyMzOtk6wLCwshl9/o6BoxYgTWrVuHJUuWYPHixYiLi8PmzZvRv39/a5tnn30WNTU1mD17NioqKjBy5EhkZmZCrTZ/eKlUKqxfvx4vvPACdDodunXrhvnz5yMjI8O+L56c3q4zV7H/wjUoPeR4anQP0eUQUROB3kpMH9EVq3acwZvZpzC+bzjkcl7+h25PJkkS+xzbQavVIiAgAJWVlZyP5KZMJgnp7+zE4YuVmJ7SFX96sP/t70REdnetRo+Rr21Hjd6IN6cO5iWA3FxrP7/ddhUb0Z366sgVHL5YCR+lAnPuiRNdDhG1IMhHiScbenhfz8yHzmAUXBE5AwYkonbQGYx4/et8AMCTo3tw5RqRg/v1qG4I81Ph4rXr+CTnguhyyAkwIBG1w793F6Ko/DrC/FT4NXfNJnJ43koPZNzbC4B5d+3K2nrBFZGjY0AiaqPK6/V4a7t55VrGvb3grRS+1oGIWuGhoZ3RK9wXldfr8c63p0WXQw6OAYmojd759jQqausRF+aLh4Z2Fl0OEbWSh0KOhffFAwA+2HUeF6/VCq6IHBkDElEbXKq4jg92ngcALLwvHh4K/goROZOxvcMwvHsw9AYTVnxzUnQ55MD47k7UBiu+LoDeYMLw7sG4J567rhM5G5lMhsUT+wAAPj9wCUcvVQquiBwVAxJRK+VeKMemA5cAAIsn9oFMxs3miJzRwM6B1r2Qln1xjJcgoWYxIBG1Qr3RhOc+PwoAeDixMwZ2DhRbEBHdkYX3xcNbqUDuhWvYmFskuhxyQAxIRK3wYcO1nAK9PbHwvj6iyyGiOxQV6IX5qeZl/8v/l4/yGr3gisjRMCAR3cbliuv4+zbzZM7F9/VBsI9ScEVE1BF+dVcs4iP8UFFbj+VbT4guhxwMAxLRbfzpy2Oo1RsxLDaIy/qJXIinQo4//3wAAGBj7kXsPVcuuCJyJAxIRLew7Xgxvj5WDA+5DC+nD+BVwIlczNCuQZiaFAMAWLL5CPQGk+CKyFEwIBG1oFZvwLIvjgEAZo7qht4RfoIrIiJbWDAhHsE+Spwsrsa/fjwnuhxyEAxIRC1Y8c1JXKq4juhAL8wbFye6HCKykUBvpXVvpDeyT+JsabXgisgRMCARNWPn6TLrN8mX0vvxemtELm7ykGjc1TMEdfUmzN9wEPVGDrW5OwYkoiYqa+vxh08PAQAeTe6Ce+LDBVdERLYmk8nw118Ogr/aA4cuVuKt7byYrbtjQCJqYsl/j0KjrUO3UB8smcQ9j4jcRWSAF15uWNW2asdp5BVeE1wRicSARPQT/z14CV8eugyFXIa/T0ng0BqRm3lgUBQeTIiC0SRh/oaDqNEZRJdEgjAgETW4VHEdSzabLycy956eSIgJFFsQEQnx4oP9ERWgxoWrtXhpy3HR5ZAgDEhEAAxGEzI2HERVnQEJMYGYM7an6JKISJAAL0/89eFBkMmA9fuKkHlUI7okEoABiQjAa5n52HOuHN5KBf4+JQEeCv5qELmzET1CMWtUdwDAHzcewukSLv13N/wUILf334OX8N4P5iX9f3loELqF+giuiIgcwR/H90ZSbDCqdQbM/mQ/tHX1oksiO2JAIrd29FIlnv3PYQDAb8f0wKSBkYIrIiJHofSQY9VjQxAZoMbZ0hrMX38QJpMkuiyyEwYkcltXq3V48pNc6AwmjOndCX8Y31t0SUTkYDr5qbDmiUSoPOTIzi/Bym0nRZdEdsKARG6p3mjC0+vycKniOrqF+uCNRwZDwQvRElEzBnQOwPJfmPdHenP7aWQevSK4IrIHBiRyO5Ik4bnPj2D32XL4KBVY88RQBHh5ii6LiBzYL4Z0xsyR3QAAGZ8e4iaSboABidyKJEl4+asT+HT/RchlwMpHBiMu3E90WUTkBBbdF4+7e3VCrd6IGR/sw4krWtElkQ0xIJFbeTP7tPUitK9NHoh7+/I6a0TUOh4KOVY/PgRDuwah8no9nvjXXpwvqxFdFtkIAxK5jfd/PIe/N0ywXHZ/X/wyMUZwRUTkbLyVHnj/V8PQJ9IfZdU6PPbPPbhSeV10WWQDDEjkFj7dX4QXGy4ZkHFvL8y4q5vgiojIWQV4eeKTmUnoHuqDSxXX8fg/96CsWie6LOpgDEjk8j7ceQ4LPjPvdfTrkd0w9x5eRoSI7kyorwqf/DoZUQFqnCmtwcP/yEFRea3osqgDMSCRy5IkCa9n5uOFL49DkoBpKV3x3KQ+kMm4nJ+I7lx0oBfWzhqO6EAvnC2tweR3d+H4ZU7cdhUMSOSS6o0m/HHjYbzz7RkAwB/H98KfHujHcEREHapbqA8++80IxEf4oaRKhyn/yEHOmauiy6IOwIBELqdWb8Dsj/fjs7yLUMhleH3yQMy5J47hiIhsIiJAjQ1PpiCpWzCqdAZMf38vvjrMzSSdHQMSuZSTxVV48O2d2FFQCrWnHGueGIqHh3G1GhHZVoCXJz7+f0mY0C8C+oad+pdvPYF6o0l0adRODEjkEiRJwqf7i/DA2z/iVEk1wvxUWDdrOMb14T5HRGQfak8FVj02BP+vYZXsP74/iyn/yMGlCm4D4IwYkMjp1egM+MOnh/Dsfw6jrt6EUXGh2DpvFIZ0CRJdGhG5GYVchqX398Xqx4fAT+2BvMIKTHzjB2w7Xiy6NGojmSRJkuginJFWq0VAQAAqKyvh7+8vuhy3tfN0GZZsPopzZTWQy4A/jO+N34zuATkvPEtEghWV12LOujwculgJAJia1AULJvRGoLdScGXurbWf3wxI7cSAJFaJtg4vf3UCXxy6DACI8FfjzamDkdQtWHBlREQ36A0mvJaZb73EUbCPEovui8fkIZ35RU4QBiQbY0ASw2A04ZPdF7Dim5Oo1hkglwHTUmKRMb4X/NWeossjImrW3nPlWLL5CE4WVwMAErsG4aX0/ugTyc8Pe2NAsjEGJPvSG0z4LO8i3vn2NIrKzRMeB8UE4s/p/dE/OkBwdUREt1dvNOGDneewctsp1OqNkMmAif0jMeeengxKdsSAZGMMSPZRV2/Ep/uLsPrbM7hcWQcACPFRImN8L0wd1oVd1ETkdC5XXMefvzqBr47c2CtpfN9wzL0nDgM68wufrTEg2RgDkm0VaKqwcX8RPj9wCVdr9ACAMD8VZt/dHY8md4G30kNwhUREd+bEFS3e3nEaW49cgeWTeFhsEH6ZGINJAyLho+L7nC0wINkYA1LHK63S4etjGmzcX2Rd9QEAUQFqPDWmBx5OjIHaUyGwQiKijne6pBrv7DiN/x66DKPJ/JHsrVRg0oBI/GJIZyTGBsFTwV15OgoDko0xIN05k0nCkUuV2J5fgm8LShqFIg+5DOP6hOGXQ2MwuncnvjkQkcvTVNZh04GL2Lj/Is6V1ViP+6k9cHevThjbOwxjendCqK9KYJXOz6kC0qpVq/CXv/wFGo0GgwYNwltvvYWkpKQW22/cuBHPP/88zp8/j7i4OLz22muYOHGi9eeSJGHZsmV47733UFFRgbvuugvvvvsu4uLirG3Ky8sxd+5cfPnll5DL5Zg8eTLeeOMN+Pr6tqpmBqS2u1ajx8GLFThYWIGDReZb5fX6Rm36R/sjPSEa6YOj+SZARG5JkiTsv3ANG/cXIet4Ma7VNn6f7Bnmi4SYQOstPsIPHvwS2WpOE5A2bNiAadOmYfXq1UhOTsbKlSuxceNGFBQUICws7Kb2u3btwt13343ly5fjZz/7GdatW4fXXnsNeXl56N+/PwDgtddew/Lly/HRRx+hW7dueP7553HkyBEcP34carUaAHDffffhypUr+Mc//oH6+nrMmDEDw4YNw7p161pVNwPSzSRJQnmNHlcq61CsrcP5q7U4W1qNM6XVOFtag5Iq3U338VN5YFSvUIzpHYYxvTohzF8toHIiIsdkNEk4dLECO/JLsKOgBEcvaW9qo1TI0TXEGz06+aJ7Jx907+SL6EAvRAaoERGg5tSEJpwmICUnJ2PYsGF4++23AQAmkwkxMTGYO3cuFi5ceFP7KVOmoKamBlu2bLEeGz58OBISErB69WpIkoSoqCj84Q9/wB//+EcAQGVlJcLDw/Hhhx/ikUcewYkTJ9C3b1/s27cPiYmJAIDMzExMnDgRFy9eRFRU1G3rtlVAKq3SQWcwtvl+Lf1ftByXIEGSAAnmIGMe5jb/aTBKMEkSjCYJBpMEg9GEeqMEvdEIvUGCzmDEdb0RtXojrtcbUas3oKrOgIraelyr1aPyej3Ka/Qo0eqgv82FGbuH+pi/9XQJxOCYIMRH+nH4jIiola5W66w98AcKK3CoqAJVOsMt7xPg5YkwPxUCvT0R4KVEkLcnArw84aPygLdSAW+lAl5KD3h5KqD0kMNTIYPSQw6lQg6FXNb4JpNBJgNkMhlkAOQNfwcAGW78d1MtHb+dUF9Vhwe81n5+C50ir9frkZubi0WLFlmPyeVypKamIicnp9n75OTkICMjo9GxtLQ0bN68GQBw7tw5aDQapKamWn8eEBCA5ORk5OTk4JFHHkFOTg4CAwOt4QgAUlNTIZfLsWfPHvz85z+/6Xl1Oh10uhs9IFrtzSm+I/xh4yF8f7LUJo9tL6G+SkQEqNE50Bs9wnwavtWYv9lwM0ciovYL8VVhXJ9w64W4TSYJlyquW3vqz5aZ/9RU1uFKZR2u1xtReb3+pukMzuLj/5eEu3t1EvLcQgNSWVkZjEYjwsMbX3E9PDwc+fn5zd5Ho9E0216j0Vh/bjl2qzZNh+88PDwQHBxsbdPU8uXL8ac//amVr6z9lAoZVB7t61FpLqH/NNHL0JD6ZebUL2/4FmD5VqCQyyCXA54K8zcH8zcJ83+bv2Eo4KP0gJdSAX+1BwK8lQj08kSgt/kW7q9GmJ8aynbWT0REbSOXyxAT7I2YYG+M6d34Z5IkQVtngKayDmXVOlRer0dFbT0qrutRWVuPGr0BtXojanVG1NYbUac3Qm80ob7hpjeYYDBJMJkkGBtGGYwmCRLMwUySAFPDMIV5dKLheSE1qeMOXl97u546ADdZaKVFixY16rnSarWIiYnp8Of55/RhHf6YRETkfmQyGQK8zMNpveEnuhynI/SrfmhoKBQKBYqLixsdLy4uRkRERLP3iYiIuGV7y5+3a1NSUtLo5waDAeXl5S0+r0qlgr+/f6MbERERuSahAUmpVGLo0KHIzs62HjOZTMjOzkZKSkqz90lJSWnUHgCysrKs7bt164aIiIhGbbRaLfbs2WNtk5KSgoqKCuTm5lrbbN++HSaTCcnJyR32+oiIiMg5CR9iy8jIwPTp05GYmIikpCSsXLkSNTU1mDFjBgBg2rRpiI6OxvLlywEA8+bNw+jRo7FixQpMmjQJ69evx/79+7FmzRoA5i7F3//+93j55ZcRFxdnXeYfFRWF9PR0AECfPn0wYcIEzJo1C6tXr0Z9fT3mzJmDRx55pFUr2IiIiMi1CQ9IU6ZMQWlpKZYuXQqNRoOEhARkZmZaJ1kXFhZCLr/R0TVixAisW7cOS5YsweLFixEXF4fNmzdb90ACgGeffRY1NTWYPXs2KioqMHLkSGRmZlr3QAKAtWvXYs6cORg3bpx1o8g333zTfi+ciIiIHJbwfZCcFTeKJCIicj6t/fzmemwiIiKiJhiQiIiIiJpgQCIiIiJqggGJiIiIqAkGJCIiIqImGJCIiIiImmBAIiIiImqCAYmIiIioCQYkIiIioiaEX2rEWVk2INdqtYIrISIiotayfG7f7kIiDEjtVFVVBQCIiYkRXAkRERG1VVVVFQICAlr8Oa/F1k4mkwmXL1+Gn58fZDJZhz2uVqtFTEwMioqKeI03G+O5tg+eZ/vgebYPnmf7sOV5liQJVVVViIqKglze8kwj9iC1k1wuR+fOnW32+P7+/vzlsxOea/vgebYPnmf74Hm2D1ud51v1HFlwkjYRERFREwxIRERERE0wIDkYlUqFZcuWQaVSiS7F5fFc2wfPs33wPNsHz7N9OMJ55iRtIiIioibYg0RERETUBAMSERERURMMSERERERNMCARERERNcGAJMCqVasQGxsLtVqN5ORk7N2795btN27ciPj4eKjVagwYMABbt261U6XOrS3n+b333sOoUaMQFBSEoKAgpKam3vb/C93Q1n/TFuvXr4dMJkN6erptC3QRbT3PFRUVePrppxEZGQmVSoVevXrx/aMV2nqeV65cid69e8PLywsxMTGYP38+6urq7FStc/r+++9x//33IyoqCjKZDJs3b77tfb799lsMGTIEKpUKPXv2xIcffmjbIiWyq/Xr10tKpVJ6//33pWPHjkmzZs2SAgMDpeLi4mbb79y5U1IoFNLrr78uHT9+XFqyZInk6ekpHTlyxM6VO5e2nudHH31UWrVqlXTgwAHpxIkT0q9+9SspICBAunjxop0rdz5tPdcW586dk6Kjo6VRo0ZJDz74oH2KdWJtPc86nU5KTEyUJk6cKP3444/SuXPnpG+//VY6ePCgnSt3Lm09z2vXrpVUKpW0du1a6dy5c9LXX38tRUZGSvPnz7dz5c5l69at0nPPPSdt2rRJAiB9/vnnt2x/9uxZydvbW8rIyJCOHz8uvfXWW5JCoZAyMzNtViMDkp0lJSVJTz/9tPXvRqNRioqKkpYvX95s+4cffliaNGlSo2PJycnSk08+adM6nV1bz3NTBoNB8vPzkz766CNblegy2nOuDQaDNGLECOmf//ynNH36dAakVmjreX733Xel7t27S3q93l4luoS2nuenn35auueeexody8jIkO666y6b1ulKWhOQnn32Walfv36Njk2ZMkVKS0uzWV0cYrMjvV6P3NxcpKamWo/J5XKkpqYiJyen2fvk5OQ0ag8AaWlpLban9p3npmpra1FfX4/g4GBblekS2nuuX3zxRYSFhWHmzJn2KNPptec8f/HFF0hJScHTTz+N8PBw9O/fH6+88gqMRqO9ynY67TnPI0aMQG5urnUY7uzZs9i6dSsmTpxol5rdhYjPQl6s1o7KyspgNBoRHh7e6Hh4eDjy8/ObvY9Go2m2vUajsVmdzq4957mpBQsWICoq6qZfSGqsPef6xx9/xL/+9S8cPHjQDhW6hvac57Nnz2L79u147LHHsHXrVpw+fRq//e1vUV9fj2XLltmjbKfTnvP86KOPoqysDCNHjoQkSTAYDHjqqaewePFie5TsNlr6LNRqtbh+/Tq8vLw6/DnZg0TUxKuvvor169fj888/h1qtFl2OS6mqqsITTzyB9957D6GhoaLLcWkmkwlhYWFYs2YNhg4diilTpuC5557D6tWrRZfmUr799lu88soreOedd5CXl4dNmzbhq6++wksvvSS6NLpD7EGyo9DQUCgUChQXFzc6XlxcjIiIiGbvExER0ab21L7zbPHXv/4Vr776KrZt24aBAwfaskyX0NZzfebMGZw/fx7333+/9ZjJZAIAeHh4oKCgAD169LBt0U6oPf+mIyMj4enpCYVCYT3Wp08faDQa6PV6KJVKm9bsjNpznp9//nk88cQT+PWvfw0AGDBgAGpqajB79mw899xzkMvZD9ERWvos9Pf3t0nvEcAeJLtSKpUYOnQosrOzrcdMJhOys7ORkpLS7H1SUlIatQeArKysFttT+84zALz++ut46aWXkJmZicTERHuU6vTaeq7j4+Nx5MgRHDx40Hp74IEHMHbsWBw8eBAxMTH2LN9ptOff9F133YXTp09bAygAnDx5EpGRkQxHLWjPea6trb0pBFlCqcRLnXYYIZ+FNpv+Tc1av369pFKppA8//FA6fvy4NHv2bCkwMFDSaDSSJEnSE088IS1cuNDafufOnZKHh4f017/+VTpx4oS0bNkyLvNvhbae51dffVVSKpXSf/7zH+nKlSvWW1VVlaiX4DTaeq6b4iq21mnreS4sLJT8/PykOXPmSAUFBdKWLVuksLAw6eWXXxb1EpxCW8/zsmXLJD8/P+n//u//pLNnz0rffPON1KNHD+nhhx8W9RKcQlVVlXTgwAHpwIEDEgDpb3/7m3TgwAHpwoULkiRJ0sKFC6UnnnjC2t6yzP+ZZ56RTpw4Ia1atYrL/F3RW2+9JXXp0kVSKpVSUlKStHv3buvPRo8eLU2fPr1R+08//VTq1auXpFQqpX79+klfffWVnSt2Tm05z127dpUA3HRbtmyZ/Qt3Qm39N/1TDEit19bzvGvXLik5OVlSqVRS9+7dpT//+c+SwWCwc9XOpy3nub6+XnrhhRekHj16SGq1WoqJiZF++9vfSteuXbN/4U5kx44dzb7nWs7t9OnTpdGjR990n4SEBEmpVErdu3eXPvjgA5vWKJMk9gESERER/RTnIBERERE1wYBERERE1AQDEhEREVETDEhERERETTAgERERETXBgERERETUBAMSERERURMMSERERERNMCARERERNcGARERERNQEAxIRERFREx6iCyAichRjxoxB//79AQCffPIJPD098Zvf/AYvvvgiZDKZ4OqIyJ7Yg0RE9BMfffQRPDw8sHfvXrzxxhv429/+hn/+85+iyyIiO5NJkiSJLoKIyBGMGTMGJSUlOHbsmLXHaOHChfjiiy9w/PhxwdURkT2xB4mI6CeGDx/eaDgtJSUFp06dgtFoFFgVEdkbAxIRERFREwxIREQ/sWfPnkZ/3717N+Li4qBQKARVREQiMCAREf1EYWEhMjIyUFBQgP/7v//DW2+9hXnz5okui4jsjMv8iYh+Ytq0abh+/TqSkpKgUCgwb948zJ49W3RZRGRnDEhERD/h6emJlStX4t133xVdChEJxCE2IiIioiYYkIiIiIia4EaRRERERE2wB4mIiIioCQYkIiIioiYYkIiIiIiaYEAiIiIiaoIBiYiIiKgJBiQiIiKiJhiQiIiIiJpgQCIiIiJq4v8D541z58vAH1sAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(ps, pdf_p)\n",
+    "plt.xlabel('p')\n",
+    "plt.ylabel('PMF');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a6227be2",
+   "metadata": {},
+   "source": [
+    "This is a beta distribution, which we can confirm by running `beta` with a change of variables, $a = k+1$ and $b = n-k+1$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 69,
+   "id": "fb3fc2ef",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(7, 7)"
+      ]
+     },
+     "execution_count": 69,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from scipy.stats import beta\n",
+    "\n",
+    "a = k + 1\n",
+    "b = n - k + 1\n",
+    "a, b"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 70,
+   "id": "18185a97",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pdf_beta = beta.pdf(ps, a, b)\n",
+    "pdf_beta /= pdf_beta.sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 71,
+   "id": "534580f9",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 71,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.allclose(pdf_p, pdf_beta)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 72,
+   "id": "26eca5b7",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(6, 12)"
+      ]
+     },
+     "execution_count": 72,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a - 1, b + k - 1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5dd2e524",
+   "metadata": {},
+   "source": [
+    "To see why this works, let's compare the PDF of the beta distribution\n",
+    "\n",
+    "$$f(p, a, b) = \\frac{1}{B(a, b)} p^{a-1} (1-p)^{b-1} $$\n",
+    "\n",
+    "And the PMF of the binomial distribution.\n",
+    "\n",
+    "$$Pr(k; n, p) = \\binom{n}{k} p^{k} (1-p)^{n-k}$$\n",
+    "\n",
+    "With the change of variables, they are identical except for the first term, which normalizes the distributions."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7fb375ee",
+   "metadata": {},
+   "source": [
+    "## Conjugate priors\n",
+    "\n",
+    "This similarity is the reason the beta and binomial are conjugate distributions, which means they are joined together.\n",
+    "This relationship has a useful property for Bayesian statistics: if the prior distribution for $p$ is beta and the likelihood of the data is binomial, the posterior distribution is also beta.\n",
+    "\n",
+    "To see how that works, suppose the prior distribution of $p$ is beta with parameters $a$ and $b$. Here is the PDF of that distribution:\n",
+    "\n",
+    "$$p^{a-1} (1-p)^{b-1}$$\n",
+    "\n",
+    "I have omitted the normalizing factor, which we don't need because we are going to normalize the distribution after the update.\n",
+    "\n",
+    "Now suppose we see $k$ successes in $n$ trials.\n",
+    "The likelihood of this data is given by the binomial distribution, which has this PMF.\n",
+    "\n",
+    "$$p^{k} (1-p)^{n-k}$$\n",
+    "\n",
+    "Again, I have omitted the normalizing factor.\n",
+    "Now to get the unnormalized posterior, we multiply the beta prior and the binomial likelihood. The result is\n",
+    "\n",
+    "$$p^{a-1+k} (1-p)^{b-1+n-k}$$\n",
+    "\n",
+    "which we recognize as an unnormalized beta distribution with parameters $a+k$ and $b+n-k$.\n",
+    "\n",
+    "So if we observe $k$ successes in $n$ trials, we can do the update by making a beta distribution with parameters $a+k$ and $b+n-k$.\n",
+    "\n",
+    "As an example, let's start with a beta prior."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 73,
+   "id": "1441b8cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a = 2\n",
+    "b = 3\n",
+    "\n",
+    "prior = beta.pdf(ps, a, b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7e2d2437",
+   "metadata": {},
+   "source": [
+    "And suppose we see 5 successes in 10 attempts."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 74,
+   "id": "a7853483",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "k = 5\n",
+    "n = 10\n",
+    "\n",
+    "like = binom.pmf(k, n, ps)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3f235ad5",
+   "metadata": {},
+   "source": [
+    "We can compute the posterior by multiplying the prior and the likelihood, then normalizing the results."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 75,
+   "id": "7209c35d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posterior = prior * like\n",
+    "posterior /= posterior.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a3eb050d",
+   "metadata": {},
+   "source": [
+    "Or we can compute a beta distribution with the updated parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 76,
+   "id": "5dbefc32",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posterior_beta = beta.pdf(ps, a+k, b+n-k)\n",
+    "posterior_beta /= posterior_beta.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "492f22dd",
+   "metadata": {},
+   "source": [
+    "The result is the same either way."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 77,
+   "id": "31c83e74",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 77,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.allclose(posterior, posterior_beta)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f8f7ea3d",
+   "metadata": {},
+   "source": [
+    "## The all-knowing cube knows all"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "id": "fca0f595",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_beta(a, b, cube):\n",
+    "    k = a - 1\n",
+    "    n = b + k - 1\n",
+    "\n",
+    "    pdf= cube[k, n, :].copy()\n",
+    "    pdf /= pdf.sum()\n",
+    "    return pdf"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 90,
+   "id": "93710b0d",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(2, 3, 5, 10)"
+      ]
+     },
+     "execution_count": 90,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a, b, k, n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 91,
+   "id": "f4ef867d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posterior_cube = get_beta(a + k, b + n - k, cube)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 93,
+   "id": "860796a2",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 93,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.allclose(posterior_beta, posterior_cube)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c41afc17",
+   "metadata": {},
+   "source": [
+    "Think Bayes, Second Edition\n",
+    "\n",
+    "Copyright 2020 Allen B. Downey\n",
+    "\n",
+    "License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "707d9f6f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}