{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bayesian Statistics Made Simple\n", "\n", "Code and exercises from my workshop on Bayesian statistics in Python.\n", "\n", "Copyright 2020 Allen Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The cookie problem\n", "\n", "> Suppose you have two bowls of cookies. Bowl 1 contains 30 vanilla and 10 chocolate cookies. Bowl 2 contains 20 vanilla of each.\n", ">\n", "> You choose one of the bowls at random and, without looking into the bowl, choose one of the cookies at random. It turns out to be a vanilla cookie.\n", ">\n", "> What is the chance that you chose Bowl 1?\n", "\n", "Assume that there was an equal chance of choosing either bowl and an equal chance of choosing any cookie in the bowl.\n", "\n", "Here are the hypotheses and prior probabilities." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "hypos = 'Bowl 1', 'Bowl 2'\n", "probs = 1/2, 1/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute the answer, I'll use a Pandas `Series` to represent the hypotheses and their probabilities." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.5\n", "Bowl 2 0.5\n", "dtype: float64" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prior = pd.Series(probs, hypos)\n", "prior" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This `Series` represents a probability mass function (PMF)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we compute the likelihood of the data under each hypothesis.\n", "\n", "* The chance of getting a vanilla cookie from Bowl 1 is 3/4.\n", "\n", "* The chance of getting a vanilla cookie from Bowl 2 is 1/2." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "likelihood = 3/4, 1/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step is to multiply the priors by the likelihoods:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.375\n", "Bowl 2 0.250\n", "dtype: float64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "unnorm = prior * likelihood\n", "unnorm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is called `unnorm` because it is an \"unnormalized posterior\".\n", "\n", "To compute the posteriors, we have to divide through by $P(D)$, which is the sum of the unnormalized posteriors." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.625" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prob_data = unnorm.sum()\n", "prob_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we get 5/8, which is what we got by computing $P(D)$ directly.\n", "\n", "Now we divide by `prob_data` to get the normalized posteriors:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.6\n", "Bowl 2 0.4\n", "dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "posterior = unnorm / prob_data\n", "posterior" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The posterior probability for Bowl 1 is 0.6, which is what we got using Bayes's Theorem explicitly.\n", "\n", "As a bonus, we also get the posterior probability of Bowl 2, which is 0.4.\n", "\n", "The posterior probabilities add up to 1, which they should, because the hypotheses are \"complementary\"; that is, either one of them is true or the other, but not both.\n", "\n", "When we add up the unnormalized posteriors and divide through, we force the posteriors to add up to 1. This process is called \"normalization\", which is why the total probability of the data is also called the \"[normalizing constant](https://en.wikipedia.org/wiki/Normalizing_constant#Bayes'_theorem)\"\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1\n", "\n", "Suppose we put the first cookie back, stir, draw another cookie from the same bowl, and it's a chocolate cookie. What is the probability we drew both cookies from Bowl 1?\n", "\n", "Hint: The prior for the second update is the posterior from the first update." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.6\n", "Bowl 2 0.4\n", "dtype: float64" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prior2 = posterior\n", "prior2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now \n", "\n", "1. Compute the likelihood of the data under each hypothesis,\n", "2. Multiply the new prior by the likelihoods.\n", "3. Divide through by the total probability of the data." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "likelihood2 = 1/4, 1/2" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.15\n", "Bowl 2 0.20\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "unnorm2 = prior2 * likelihood2\n", "unnorm2" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.35" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "prob_data2 = unnorm2.sum()\n", "prob_data2" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Bowl 1 0.428571\n", "Bowl 2 0.571429\n", "dtype: float64" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "posterior2 = unnorm2 / prob_data2\n", "posterior2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 101 Bowls\n", "\n", "Suppose instead of 2 bowls there are 101 bowls:\n", "\n", "* Bowl 0 contains no vanilla cookies,\n", "\n", "* Bowl 1 contains 1% vanilla cookies,\n", "\n", "* Bowl 2 contains 2% vanilla cookies,\n", "\n", "and so on, up to\n", "\n", "* Bowl 99 contains 99% vanilla cookies, and\n", "\n", "* Bowl 100 contains all vanilla cookies.\n", "\n", "As in the previous problem, there are only two kinds of cookies, vanilla and chocolate. So Bowl 0 is all chocolate cookies, Bowl 1 is 99% chocolate, and so on.\n", "\n", "Suppose we choose a bowl at random, choose a cookie at random, and it turns out to be vanilla. What is the probability that the cookie came from Bowl $x$, for each value of $x$?\n", "\n", "To represent the prior, I'll use a Pandas `Series` with 101 equally spaced quantities from 0 to 1." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.00 0.009901\n", "0.01 0.009901\n", "0.02 0.009901\n", "0.03 0.009901\n", "0.04 0.009901\n", "dtype: float64" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xs = np.linspace(0, 1, num=101)\n", "prob = 1/101\n", "\n", "prior = pd.Series(prob, xs)\n", "prior.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what it looks like." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbgElEQVR4nO3de7RkZX3m8e/DzQsCDXbjIlxsTBqSNoriEZmMTogZDI3RxqgRkgmgZrWMkJvR2JlFYoi5oCuJCRFhkHBLJjJEQduRyDDthYgQOY3cGkTaFqSFJW10YfBGwN/8sffBsqhTp06z6xyK/n7WqlW79vu+u963Guo5+1LvTlUhSVIXdljsDkiSnjgMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBXpcSbJA0metdj9kLaFoSItgCR3JvluGxhfS3J+kqcNqltVT6uqzQvdR6kLhoq0cF5RVU8DDgVeCJzaW5hkp8ey8cfaXuqCoSItsKr6KvDPwE8nqSQnJ7kDuAOgXfcT7fIeSS5KsjXJXUlOTbJDW3ZikquTvCfJN4A/WqQhSY/wLxtpgSXZHzgauBR4OXAM8CLguwOq/y2wB/As4OnA/wXuBf6uLX8RcDGwN7DzOPstjSLO/SWNX5I7gaXAQ8D9wMeA3wW+A/x8VX2ip24BK4Avt+XPr6pb27I3AcdV1RFJTgT+uKoOWMChSEO5pyItnGOq6v/1rkgCcPcs9ZcCuwB39ay7C9i35/VsbaVF4TkVafHNdrjg68B/AM/sWXcA8NUR2kqLwlCRHqeq6mHgEuBPk+yW5JnAW4B/WNyeSbMzVKTHt98Avg1sBj4D/CNw3qL2SBrCE/WSpM64pyJJ6oyhIknqjKEiSeqMoSJJ6sx2/ePHpUuX1vLlyxe7G5I0UTZs2PD1qlo2qGy7DpXly5czPT292N2QpImS5K7Zyjz8JUnqjKEiSeqMoSJJ6oyhIknqzFhDJclRSW5PsinJ2gHlSXJGW35TkkN7ys5Lcl+SW/ra7JXkyiR3tM979pUf0N4H/K3jG5kkaZCxhUqSHYEzgVXASuC4JCv7qq2iuRnRCmANcFZP2QXAUQM2vRZYX1UrgPXt617voblVqyRpgY1zT+UwYFNVba6qB2luebq6r85q4KJqXAssSbIPQFVdBXxjwHZXAxe2yxfS3IoVgCTH0MzmurHDcUiSRjTOUNmXH70r3RZ+9I51o9bp94yquhegfd4bIMmuwNuB04Y1TrImyXSS6a1bt845CEnS6MYZKhmwrn+e/VHqjOo04D1V9cCwSlV1TlVNVdXUsmUDfxAqSdpG4/xF/RZg/57X+wH3bEOdfl9Lsk9V3dseKruvXf8i4DVJ3g0sAX6Q5HtV9d5tHYAkaX7GuadyHbAiyYFJdgGOBdb11VkHHN9eBXY4cP/Moa0h1gEntMsnAB8BqKqXVNXyqloO/DXwZwaKJC2ssYVKVT0EnAJcAdwGXFJVG5OclOSkttrlNCfWNwHvB9480z7JB4BrgIOTbEnyxrbodODIJHcAR7avJUmPA9v17YSnpqbKCSUlaX6SbKiqqUFl/qJektQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFktSZsYZKkqOS3J5kU5K1A8qT5Iy2/KYkh/aUnZfkviS39LXZK8mVSe5on/ds1x+ZZEOSm9vnl45zbJKkRxtbqCTZETgTWAWsBI5LsrKv2ipgRftYA5zVU3YBcNSATa8F1lfVCmB9+xrg68Arquo5wAnA33czEknSqMa5p3IYsKmqNlfVg8DFwOq+OquBi6pxLbAkyT4AVXUV8I0B210NXNguXwgc09b/fFXd067fCDw5yZO6HJAkabhxhsq+wN09r7e06+Zbp98zqupegPZ57wF1Xg18vqq+31+QZE2S6STTW7duneOtJEnzMc5QyYB1tQ115vemybOBdwFvGlReVedU1VRVTS1btuyxvJUkqc84Q2ULsH/P6/2Ae7ahTr+vzRwia5/vmylIsh9wGXB8VX1pG/stSdpG4wyV64AVSQ5MsgtwLLCur8464Pj2KrDDgftnDm0NsY7mRDzt80cAkiwBPgb8flVd3dEYJEnzMLZQqaqHgFOAK4DbgEuqamOSk5Kc1Fa7HNgMbALeD7x5pn2SDwDXAAcn2ZLkjW3R6cCRSe4Ajmxf077XTwB/kOSG9jHofIskaUxS9ZhOYUy0qampmp6eXuxuSNJESbKhqqYGlfmLeklSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmdGCpUkv5jEAJIkDTVqUBwL3JHk3Ul+atSNJzkqye1JNiVZO6A8Sc5oy29KcmhP2XlJ7ktyS1+bvZJcmeSO9nnPnrLfb7d1e5JfGLWfkqRujBQqVfXfgOcDXwLOT3JNkjVJdputTZIdgTOBVcBK4LgkK/uqrQJWtI81wFk9ZRcARw3Y9FpgfVWtANa3r2m3fSzw7Lbd+9o+SJIWyMiHtKrqW8CHgIuBfYBXAdcn+Y1ZmhwGbKqqzVX1YNtudV+d1cBF1bgWWJJkn/b9rgK+MWC7q4EL2+ULgWN61l9cVd+vqi8Dm9o+SJIWyKjnVF6Z5DLgE8DOwGFVtQo4BHjrLM32Be7ueb2lXTffOv2eUVX3ArTPe89nW+0e1nSS6a1bt87xVpKk+dhpxHqvAd7T7j08oqq+k+QNs7TJgHW1DXVGNdK2quoc4ByAqampbX0vSdIAox7+urc/UJK8C6Cq1s/SZguwf8/r/YB7tqFOv6/NHCJrn+97DNuSJHVo1FA5csC6VXO0uQ5YkeTAJLvQnERf11dnHXB8exXY4cD9M4e2hlgHnNAunwB8pGf9sUmelORAmpP/n5tjW5KkDg09/JXkvwNvBn48yU09RbsBVw9rW1UPJTkFuALYETivqjYmOaktPxu4HDia5qT6d4DX97z3B4AjgKVJtgDvqKq/A04HLknyRuArwGvb7W1McglwK/AQcHJVPTzSpyBJ6kSqZj+tkGQPYE/gz2kv3W39e1UNujJrokxNTdX09PRid0OSJkqSDVU1NahsrhP1VVV3Jjl5wEb3eiIEiySpO3OFyj8CvwhsoLmSqvcKqwKeNaZ+SZIm0NBQqapfbJ8PXJjuSJIm2Vwn6g8dVl5V13fbHUnSJJvr8NdfDikr4KUd9kWSNOHmOvz1cwvVEUnS5Jvr8NdLq+oTSX5pUHlVXTqebkmSJtFch79+lmYSyVcMKCvAUJEkPWKuw1/vaJ9fP6yeJEkw+tT3T2/v0Hh9kg1J/ibJ08fdOUnSZBl1QsmLga3Aq2mmwd8K/O9xdUqSNJlGvZ/KXlX1zp7Xf5LkmDH0R5I0wUbdU/lkkmOT7NA+fhn42Dg7JkmaPHNdUvzv/HDOr7cA/9AW7QA8ALxjrL17HDvtoxu59Z5vLXY3JGmbrPyx3XnHK57d+Xbnuvprt87fUZL0hDXqORWS7ElzN8Unz6zrv8Xw9mQcCS9Jk26kUEny68Bv0dz3/QbgcOAanPtLktRj1BP1vwW8ELirnQ/s+TSXFUuS9IhRQ+V7VfU9gCRPqqovAAePr1uSpEk06jmVLUmWAB8GrkzyTeCecXVKkjSZRgqVqnpVu/hHST4J7AF8fGy9kiRNpPlc/XUo8GKa361cXVUPjq1XkqSJNOqEkn8IXAg8HVgKnJ/k1HF2TJI0eUbdUzkOeH7PyfrTgeuBPxlXxyRJk2fUq7/upOdHj8CTgC913htJ0kSba+6vv6U5h/J9YGOSK9vXRwKfGX/3JEmTZK7DX9Pt8wbgsp71nxpLbyRJE22uCSUvnFlOsgtwUPvy9qr6j3F2TJI0eUad++sImqu/7qSZBn//JCdszxNKSpIebdSrv/4SeFlV3Q6Q5CDgA8ALxtUxSdLkGfXqr51nAgWgqr4I7DyeLkmSJtWoobIhyd8lOaJ9vJ/m5P1QSY5KcnuSTUnWDihPkjPa8pvaX+0PbZvkkCTXJLk5yUeT7N6u3znJhe3625L8/ohjkyR1ZNRQOQnYCPwmzTT4t7brZpVkR+BMYBWwEjguycq+aqtobvy1AlgDnDVC23OBtVX1HJor0t7Wrn8t8KR2/QuANyVZPuL4JEkdmPOcSpIdgA1V9dPAX81j24cBm6pqc7udi4HVNIE0YzVwUVUVcG2SJUn2AZYPaXswMHOBwJXAFcAf0Px+ZtckOwFPAR4EvIm8JC2gOfdUquoHwI1JDpjntvcF7u55vaVdN0qdYW1vAV7ZLr8W2L9d/iDwbeBe4CvAX1TVN/o7lWRNkukk01u3ep8xSerSqIe/9qH5Rf36JOtmHnO0yYB1NWKdYW3fAJycZAOwG80eCTR7Rg8DPwYcCPxukmc9aiNV51TVVFVNLVu2bI4hSJLmY9RLik/bhm1v4Yd7EdDc377/xl6z1dlltrbtXSdfBo9c2vzyts6vAB9vf5R5X5KrgSlg8zb0XZK0DYbuqSR5cpLfpjnM9JM091H59Mxjjm1fB6xIcmD7a/xjgf69m3XA8e1VYIcD91fVvcPaJtm7fd4BOBU4u93WV4CXttvaFTgc+MIIn4EkqSNzHf66kOav/ZtprsT6y1E3XFUPAafQnEi/DbikqjYmOSnJzJVjl9PsSWwC3g+8eVjbts1xSb5IExj3AOe3688EnkZzzuU64PyqumnU/kqSHrs0F17NUpjc3F6iS3tV1eeq6tBZG0yYqampmp6enruiJOkRSTZU1dSgsrn2VB6ZNLLde5AkaVZznag/JMnMbz0CPKV9HaCqavex9k6SNFHmmvp+x4XqiCRp8o36OxVJkuZkqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjpjqEiSOmOoSJI6Y6hIkjoz1lBJclSS25NsSrJ2QHmSnNGW35Tk0LnaJjkkyTVJbk7y0SS795Q9ty3b2JY/eZzjkyT9qLGFSpIdgTOBVcBK4LgkK/uqrQJWtI81wFkjtD0XWFtVzwEuA97WttkJ+AfgpKp6NnAE8B/jGp8k6dHGuadyGLCpqjZX1YPAxcDqvjqrgYuqcS2wJMk+c7Q9GLiqXb4SeHW7/DLgpqq6EaCq/q2qHh7X4CRJjzbOUNkXuLvn9ZZ23Sh1hrW9BXhlu/xaYP92+SCgklyR5PokvzeoU0nWJJlOMr1169Z5DkmSNMw4QyUD1tWIdYa1fQNwcpINwG7Ag+36nYAXA7/aPr8qyc8/aiNV51TVVFVNLVu2bO5RSJJGttMYt72FH+5FAOwH3DNinV1ma1tVX6A51EWSg4CX92zr01X19bbscuBQYH0HY5EkjWCceyrXASuSHJhkF+BYYF1fnXXA8e1VYIcD91fVvcPaJtm7fd4BOBU4u93WFcBzkzy1PWn/s8CtYxyfJKnP2PZUquqhJKfQfNnvCJxXVRuTnNSWnw1cDhwNbAK+A7x+WNt208clObldvhQ4v23zzSR/RRNIBVxeVR8b1/gkSY+Wqv7THNuPqampmp6eXuxuSNJESbKhqqYGlfmLeklSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZwwVSVJnDBVJUmcMFUlSZ8YaKkmOSnJ7kk1J1g4oT5Iz2vKbkhw6V9skhyS5JsnNST6aZPe+bR6Q5IEkbx3n2CRJjza2UEmyI3AmsApYCRyXZGVftVXAivaxBjhrhLbnAmur6jnAZcDb+rb5HuCfOx+QJGlO49xTOQzYVFWbq+pB4GJgdV+d1cBF1bgWWJJknznaHgxc1S5fCbx6ZmNJjgE2AxvHNCZJ0hDjDJV9gbt7Xm9p141SZ1jbW4BXtsuvBfYHSLIr8HbgtGGdSrImyXSS6a1bt448GEnS3MYZKhmwrkasM6ztG4CTk2wAdgMebNefBrynqh4Y1qmqOqeqpqpqatmyZcOqSpLmaacxbnsL7V5Eaz/gnhHr7DJb26r6AvAygCQHAS9v67wIeE2SdwNLgB8k+V5VvbeLwUiS5jbOULkOWJHkQOCrwLHAr/TVWQeckuRimlC4v6ruTbJ1trZJ9q6q+5LsAJwKnA1QVS+Z2WiSPwIeMFAkaWGNLVSq6qEkpwBXADsC51XVxiQnteVnA5cDRwObgO8Arx/Wtt30cUlObpcvBc4f1xgkSfOTqv7THNuPqampmp6eXuxuSNJESbKhqqYGlfmLeklSZwwVSVJnDBVJUme263Mq7VVmdz2GTSwFvt5RdybB9jZecMzbC8c8P8+sqoE/9NuuQ+WxSjI928mqJ6LtbbzgmLcXjrk7Hv6SJHXGUJEkdcZQeWzOWewOLLDtbbzgmLcXjrkjnlORJHXGPRVJUmcMFUlSZwyVOSQ5KsntSTYlWTugPEnOaMtvSnLoYvSzSyOM+Vfbsd6U5LNJDlmMfnZprjH31HthkoeTvGYh+zcOo4w5yRFJbkiyMcmnF7qPXRvhv+09knw0yY3tmF+/GP3sSpLzktyX5JZZyrv//qoqH7M8aGZI/hLwLJp7vNwIrOyrczTwzzQ3Fjsc+NfF7vcCjPlngD3b5VXbw5h76n2CZnbt1yx2vxfg33kJcCtwQPt678Xu9wKM+X8A72qXlwHfAHZZ7L4/hjH/F+BQ4JZZyjv//nJPZbjDgE1VtbmqHgQuBlb31VkNXFSNa4ElSfZZ6I52aM4xV9Vnq+qb7ctraW6iNslG+XcG+A3gQ8B9C9m5MRllzL8CXFpVXwGoqkkf9yhjLmC3JAGeRhMqDy1sN7tTVVfRjGE2nX9/GSrD7Qvc3fN6S7tuvnUmyXzH80aav3Qm2ZxjTrIv8Cram8I9AYzy73wQsGeSTyXZkOT4BevdeIwy5vcCP0Vzp9mbgd+qqh8sTPcWReffX+O88+MTQQas678Ge5Q6k2Tk8ST5OZpQefFYezR+o4z5r4G3V9XDzR+xE2+UMe8EvAD4eeApwDVJrq2qL467c2Myyph/AbgBeCnw48CVSf6lqr415r4tls6/vwyV4bYA+/e83o/mL5j51pkkI40nyXOBc4FVVfVvC9S3cRllzFPAxW2gLAWOTvJQVX14QXrYvVH/2/56VX0b+HaSq4BDgEkNlVHG/Hrg9GpOOGxK8mXgJ4HPLUwXF1zn318e/hruOmBFkgOT7AIcC6zrq7MOOL69iuJw4P6qunehO9qhOcec5ACaWzn/2gT/1dprzjFX1YFVtbyqlgMfBN48wYECo/23/RHgJUl2SvJU4EXAbQvczy6NMuav0OyZkeQZwMHA5gXt5cLq/PvLPZUhquqhJKcAV9BcOXJeVW1MclJbfjbNlUBHA5uA79D8pTOxRhzzHwJPB97X/uX+UE3wDK8jjvkJZZQxV9VtST4O3AT8ADi3qgZemjoJRvx3fidwQZKbaQ4Nvb2qJnZK/CQfAI4AlibZArwD2BnG9/3lNC2SpM54+EuS1BlDRZLUGUNFktQZQ0WS1BlDRZLUGUNFj1vtbMA39DyWP8btPS/J0T2vXzlsRuIuJPnNJLcl+V9jfI/Lkyxplx9on5fPNjPtOCW5YNAMzknOTbJyofujhefvVPR49t2qet6ggnbCv8xzXqbn0fwy/nKAqlrHo3/81rU308w68OVxvUFVHT13rcVVVb++2H3QwnBPRROj/ev7tiTvA64H9k9yVpLp9t4Xp/XUfWGae73cmORzSfYA/hh4XbvX87okJyZ5b1v/mUnWt/eUWN/OGjDzl/cZ7bY2D/orvK33liS3tI/fbtedTTPN+rokv9NX/1+TPLvn9aeSvCDJYe17fb59PrgtPzHJpUk+nuSOJO/uaXtnkqVzfG7/kuT69vEzs9Q7vh3/jUn+fo7PZeD6vu29s/38dmjHN9Wuf1mSa9q+/FOSp7XrT09ya7vNv5htPHqcW+z5/n34mO0BPEwzud8NwGXAcppfdh/eU2ev9nlH4FPAc2nulbEZeGFbtjvNXvmJwHt72j7yGvgocEK7/Abgw+3yBcA/0fwBtpJm6vT+fr6AZkbbXWmmS98IPL8tuxNYOqDN7wCntcv7AF/s7Wu7/F+BD/X0dTOwB/Bk4C5g//73AB5on5fT3kMDeCrw5HZ5BTA9oD/PBm7v2c5ec3wuwz6v1wDvBv4nP/yB9ado9hKXAlcBu7br304zQ8Ne7fvP1F+y2P/9+di2h4e/9Hj2I4e/2nMqd1Vz34cZv5xkDU1o7EPzxV/AvVV1HUC1M8xm+OzC/wn4pXb572m+FGd8uJrDbLemmQ+q34uBy6qZeJEklwIvAT4/5P0uAa6kmTbjl2mCC5rQuDDJinYcO/e0WV9V97fvcSvwTH502vLZ7Ay8N8nzaIL6oAF1Xgp8sNopSapq5h4cs30uwz6vP6C52dOaAe9zOM2/0dXtv8cuwDXAt4DvAecm+Rjwf0YYlx6HDBVNmm/PLCQ5EHgrzR7JN5NcQPNXfHjstx/obf/9nuVByTTvufCr6qtJ/i3NbM+vA97UFr0T+GRVvaoN0U/N0o+HGf3/398BvkYzw/AONF/e/Ub9zGar07v+OuAFSfbqCafe97myqo57VAeSw2gmczwWOIUm6DRhPKeiSbY7Tcjc3+5BrGrXfwH4sSQvBEiyW5KdgH8HdptlW5+l+TID+FXgM/Pox1XAMUmemmRXmpt5/csI7S4Gfg/Yo6pubtftAXy1XT5xHn0YZg+aPbcfAL9Gc6iw33qavb6nAyTZq10/2+cy7PP6OHA68LEk/Z/3tcB/TvIT7fs8NclB7XmVParqcuC3aS6q0ARyT0UTq6puTPJ5mnMYm4Gr2/UPJnkd8LdJngJ8l+b8xCeBtUluAP68b3O/CZyX5G3AVuYxW2tVXd/uJc3cc+Pcqhp26GvGB4G/odk7mfFumsNfbwE+MWof5vA+4ENJXkvzGXy7v0I1s/X+KfDpJA/THLo7kdk/l6GfV1X9Uxso69JzGXdVbU1yIvCBJE9qV59KE/gfSTKzp/kjFzZocjhLsSSpMx7+kiR1xlCRJHXGUJEkdcZQkSR1xlCRJHXGUJEkdcZQkSR15v8DgqkRX09lQdcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "prior.plot()\n", "\n", "plt.xlabel('Fraction of vanilla cookies')\n", "plt.ylabel('Probability')\n", "plt.title('Prior');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a prior, we need to compute likelihoods.\n", "\n", "Here are the likelihoods for a vanilla cookie:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "likelihood_vanilla = xs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And for a chocolate cookie." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "likelihood_chocolate = 1 - xs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute unnormalized posteriors, we multiply the priors and the likelihoods." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "unnorm = prior * likelihood_vanilla" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To normalize, we divide through by the total probability of the data." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "posterior = unnorm / unnorm.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the posterior looks like." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA30ElEQVR4nO3dd3zUZbbH8c+X3kGkSAtFQKkihKJrr4Aooq6irth2kbt6d93duwKWFdsuupa1c1Gx7LqiAiIqLmJBbChBMdRIKEIgUqXXJOf+8fvlOsSQDDCTySTn/XrllZl5nuf3O0/EOfMrcx6ZGc4551wsVEh0AM4558oOTyrOOedixpOKc865mPGk4pxzLmY8qTjnnIsZTyrOOedixpOKcwdJ0pWS3kt0HKWVpHclXR0+vkbSpxFtJqlt4qJz8eZJxZVqklZI2iVpu6S1kp6XVOswtjdK0r8OJyYze9nMzjmcbZRlZtbPzF5MdBwuMTypuGRwvpnVAroDPYHbExWIpEqHMVaS/P85V6b5P3CXNMxsNfAu0BlA0gWSFkjaLGmGpA75fSUNl7Ra0jZJGZLOlNQXuBW4LDzy+TbsW1fSc5KywzH3SqoYtl0j6TNJj0jaBIwq5JTOiZJmS9oS/j4xom2GpPskfQbsBNoUN89i5rVC0v9ISg/396qkahHtAyTNDcd+LqnrAfYxRtKDBV57U9Ifw8cjJC0N/34LJQ2K6HeNpE8lPSjpR0nLJfUrMOdfRzHP8yR9I2mrpFWSRhU3xpV+nlRc0pDUAugPfCOpPfAKcDPQEJgKvCWpiqRjgJuAnmZWGzgXWGFm/wH+CrxqZrXM7Lhw0y8COUBb4HjgHCDyTbE3sAxoBNxXIKb6wDvAY8CRwMPAO5KOjOh2FTAUqA18X8wcDziviG6XAn2B1kBX4JpwbHdgHHBDGMv/AlMkVS1kV/8mSK4Kxx4Rznt82L4UOBmoC9wF/EtSkwJ/kwygAfAA8Fz+tg7CDmAIUA84D/gvSRce5DZcKeNJxSWDyZI2A58CHxMkhsuAd8xsupntAx4EqgMnArlAVaCjpMpmtsLMlha2YUmNgX7AzWa2w8zWAY8AgyO6rTGzx80sx8x2FdjEecASM/tn2P4KsBg4P6LPC2a2IGzfV8xci5pXvsfMbI2ZbQLeArqFr/8G+F8z+9LMcsPrGnuAPoXs5xPACBIHwCXAF2a2BsDMXg/3kWdmrwJLgF4R4783s2fMLJcgKTcBGhczt/2Y2QwzmxfuI50gmZ56MNtwpY8nFZcMLjSzembW0sx+G76xNyXiU7+Z5QGrgGZmlknwSX8UsE7SeElND7DtlkBlIDs8ZbSZ4BN+o4g+q4qIbb84Qt8DzaIcX+T2IucV0eeHiMc7gfwbF1oCf8qfRziXFuE292NBJdnxwOXhS1cAL+e3SxoScRptM8EpxwaFxWBmO8OHB3UDhaTekj6StF7SFmBYgX24JORJxSWrNQRvokBwEZzgDXQ1gJn928xOCvsYcH/YtWBZ7lUEn+YbhImrnpnVMbNOEX2KKuW9XxyhlPw4ohhf5PYKzqsYq4D7IuZRz8xqhEdPhXkFuERSS4LTWRPDfbYEniE4hXikmdUD5gMHe3qrOP8GpgAtzKwuMCYO+3AlzJOKS1avAeeFF+ArA38iSA6fSzpG0hnhtYTdwC6CU2IAa4FW+XdhmVk28B7wkKQ6kipIOlpStKdhpgLtJV0hqZKky4COwNsHGqDgtuYZBzuvKGJ5BhgWHgFIUs3wYnjtwjqb2TfAeuBZYJqZbQ6bahIkwvVhvNcS3hwRY7WBTWa2W1IvgqMll+Q8qbikZGYZwK+Ax4ENBNcwzjezvQTXU0aHr/9AcCrr1nDo6+HvjZK+Dh8PAaoAC4EfgQkE1wiiiWMjMIDgzX8jcAswwMw2FDGsBfDZIcyruFjSCK6rPBHOI5PwIn4RXgHOIjhqyN/OQuAh4AuCJNzlQPEept8Cd0vaBvyFIKG6JCdfpMu5kiVpLnBmmJCcK1M8qTjnnIsZP/3lnHMuZjypOOecixlPKs4552LmkIvjlQUNGjSwVq1aJToM55xLKnPmzNlgZg0LayvXSaVVq1akpaUlOgznnEsqkg5Yw85PfznnnIsZTyrOOedixpOKc865mIlrUpHUV8ECSZmSRhTSLkmPhe3p4XoQSGoRVi9dFC5W9PuIMfUlTZe0JPx9RETbyHBbGZLOjefcnHPO/VzckoqClfOeJFiroiNwuaSOBbr1A9qFP0OBp8PXc4A/mVkHgrUgbowYOwL4wMzaAR+EzwnbBwOdCBYweiqMwTnnXAmJ55FKLyDTzJaFxfDGAwML9BkIvGSBWUA9SU3MLNvMvgYws23AIn5aT2IgwaJAhL8vjHh9vJntMbPlBMX0IhcVcs45F2fxTCrN2H9xoiz2X2goqj6SWhEs8fpl+FLjsFx5ftny/MWUotkfkoZKSpOUtn79+oOZj3POuWLEM6kUtthOweqVRfaRVItg4aCbzWxrDPaHmY01s1QzS23YsNDv7jjnXJllZrw6eyXvL1wbl+3HM6lkEawbka85wap2UfUJFyiaCLxsZpMi+qyV1CTs0wRYdxD7c865cmvlxp1c+eyXDJ84j8lzo1lM9ODFM6nMBtpJai2pCsFF9CkF+kwBhoR3gfUBtphZdriE6nPAIjN7uJAxV4ePrwbejHh9sKSqkloTXPz/KvbTcs655JKbZzz7yTLO+cfHpGdt4a+DuvDY4OPjsq+4lWkxsxxJNwHTgIrAODNbIGlY2D6GYCnW/gQX1XcC14bDfwFcBcwLFzQCuNXMphKs6PeapOuBlcAvw+0tkPQawep9OcCNZpa/hKxzzpVL363dxi0T0pm7ajNnHNuI+wZ1pknd6nHbX7lepCs1NdW89pdzrizam5PHUzMyefKjTGpVrcSoCzpxwXFNCU4EHR5Jc8wstbC2cl1Q0jnnyqJvV23mlgnpZKzdxgXHNeXO8ztyZK2qJbJvTyrOOVdG7Nqby8PTM3ju0+U0ql2NZ4ekclbHxiUagycV55wrA75YupERk9L5fuNOLu/VgpH9O1CnWuUSj8OTinPOJbGtu/cx+t3F/PvLlbQ8sgb//k1vTjy6QcLi8aTinHNJ6sPFa7l10nzWbdvNb05uzR/PPobqVRJb8tCTinPOJZmN2/dw99sLeXPuGo5pXJsxV/WgW4t6iQ4L8KTinHNJw8yY8u0a7nprIdt27+Pms9rx29PaUqVS6Vkay5OKc84lgewtu7hj8nzeX7SO41rU44GLu3LMUbUTHdbPeFJxzrlSLC/PGD97FX+buoh9eXncfl4Hrv1FaypWOPwvMcaDJxXnnCulVmzYwYhJ6cxatokTjz6S0Rd1JeXIGokOq0ieVJxzrpTJyc1j3GfLeei976hSsQKjL+rCZT1bxKTESrx5UnHOuVJkUfZWhk9MJz1rC2d1aMy9F3bmqLrVEh1W1DypOOdcKbAnJ5cnP1rKUx9lUrd6ZR6//HgGdG2SFEcnkTypOOdcgn2z8keGT0znu7XbGXR8M+4Y0JH6NaskOqxD4knFOecSZOfeHB567zvGfbaco+pU4/lrenL6sY0SHdZh8aTinHMJ8HnmBkZMmsfKTTv5VZ8Uhvc9ltoJKAAZa3H9GqakvpIyJGVKGlFIuyQ9FranS+oe0TZO0jpJ8wuMeVXS3PBnRf7KkJJaSdoV0TYmnnNzzrlDsWXXPkZMTOeKZ7+kguDVoX2498IuZSKhQByPVCRVBJ4EzgaygNmSppjZwohu/QjWkm8H9AaeDn8DvAA8AbwUuV0zuyxiHw8BWyKal5pZt5hOxDnnYuS9BT9w++T5bNyxlxtObcMfzmpPtcqJLQAZa/E8/dULyDSzZQCSxgMDCdaQzzcQeMmCNY1nSaonqYmZZZvZTEmtDrRxBbdEXAqcEbcZOOdcDGzYvodRUxbwdno2HZrU4bmre9Kled1EhxUX8UwqzYBVEc+z+OkopKg+zYDsKLZ/MrDWzJZEvNZa0jfAVuB2M/uk4CBJQ4GhACkpKVHsxjnnDo2ZMXnuau56ayE79+Typ7PbM+y0o6lcsfQUgIy1eCaVwm6utkPocyCXA69EPM8GUsxso6QewGRJncxs634bNxsLjAVITU2Ndl/OOXdQ1mzexW1vzOOjjPV0T6nHA5d0pW2j0lcAMtbimVSygBYRz5sDaw6hz89IqgRcBPTIf83M9gB7wsdzJC0F2gNphxK8c84dirw84+WvVjJ66iLyDO48vyNDTmhVagtAxlo8k8psoJ2k1sBqYDBwRYE+U4CbwustvYEtZhbNqa+zgMVmlpX/gqSGwCYzy5XUhuDi/7IYzMM556KybP12Rkycx1crNnFS2wb87aIutKhfugtAxlrckoqZ5Ui6CZgGVATGmdkCScPC9jHAVKA/kAnsBK7NHy/pFeA0oIGkLOBOM3subB7M/qe+AE4B7paUA+QCw8xsU7zm55xz+XJy83j20+U8Mv07qlaqwAOXdOWXPZonXYmVWFBw41X5lJqaamlpfnbMOXfoFq7Zyi0Tv2X+6q2c26kx9wzsTKM6yVMA8lBImmNmqYW1+TfqnXPuEOzJyeXxDzIZ8/FS6tWowtNXdqdflyaJDivhPKk459xBmvP9Jm6ZkM7S9Tu4qHsz/jKgI/VqJGcByFjzpOKcc1HasSeHv0/L4MUvVtC0bnVevK4Xp7ZvmOiwShVPKs45F4VPlqxn5KR5rN68i6v6tOSWvsdSq6q/hRbkfxHnnCvClp37uPedhbw+J4s2DWvy2g0n0LNV/USHVWp5UnHOuQP4z/xs7nhzAZt27OW3px3N785sV+YKQMaaJxXnnCtg3bbd3PnmAt6d/wMdm9Th+Wt60rlZ2SwAGWueVJxzLmRmTJiTxb3vLGLXvlz+fO4xDD2lTZkuABlrnlSccw5YtWknt74xj0+WbCC15RGMvrgrbRvVSnRYSceTinOuXMvLM176YgUPTMtAwF0XdOKqPi2pUE4KQMaaJxXnXLmVuW47Iyamk/b9j5zSviF/HdSZ5keUrwKQseZJxTlX7uzLzWPszGU8+v4SqlepyIO/PI6LuzcrlwUgY82TinOuXJm/egu3TEhnYfZWzuvShFEXdKJh7aqJDqvM8KTinCsXdu/L5dEPljB25jLq16zCmF/1oG/noxIdVpnjScU5V+bNXrGJ4RPSWbZhB7/s0Zzbz+tI3RqVEx1WmeRJxTlXZm3fk8MD/1nMS198T/MjqvPP63txcjsvABlPcf1Gj6S+kjIkZUoaUUi7JD0WtqdL6h7RNk7SOknzC4wZJWm1pLnhT/+ItpHhtjIknRvPuTnnSrcZGes495GZ/HPW91z7i1ZMu/kUTyglIG5HKpIqAk8CZwNZwGxJU8xsYUS3fgRrybcjWKP+6fA3wAvAE8BLhWz+ETN7sMD+OhIsM9wJaAq8L6m9meXGbFLOuVLvxx17ueedhUz6ejVHN6zJhGEn0KOlF4AsKfE8/dULyDSzZQCSxgMDgcikMhB4yYI1jWdJqiepiZllm9lMSa0OYn8DgfFmtgdYLikzjOGLWEzGOVe6mRlT5/3AnVPms3nnPv77jLbcdEZbqlbyApAlKZ5JpRmwKuJ5Fj8dhRTVpxmQXcy2b5I0BEgD/mRmP4bjZhWyrf1IGgoMBUhJSSl+Fs65Um/d1t3c8eZ8pi1YS5dmdXnput50bFon0WGVS/G8plLYt4jsEPoU9DRwNNCNIPk8dDDbMrOxZpZqZqkNG/r5VeeSmZnx2uxVnPXwx8zIWM/Ifsfyxm9P9ISSQPE8UskCWkQ8bw6sOYQ++zGztfmPJT0DvH2o23LOJa9Vm3YyctI8Ps3cQK/W9Rl9URfaNPQCkIkWzyOV2UA7Sa0lVSG4iD6lQJ8pwJDwLrA+wBYzK/LUl6QmEU8HAfl3h00BBkuqKqk1wcX/r2IxEedc6ZGbZ4z7dDnnPDKTuas2c++FnRn/mz6eUEqJuB2pmFmOpJuAaUBFYJyZLZA0LGwfA0wF+gOZwE7g2vzxkl4BTgMaSMoC7jSz54AHJHUjOLW1Argh3N4CSa8R3AiQA9zod345V7YsWbuNWyam883KzZx+TEPuG9SFpvWqJzosF0HBjVflU2pqqqWlpSU6DOdcMfbl5vH0jKU88WEmNatW5M7zOzGwW1MvAJkgkuaYWWphbf6NeudcqTYvawt/nvAti3/YxoCuQQHIBrW8AGRp5UnFOVcq7d6Xyz/eX8IznyzjyJpVGHtVD87p5AUgSztPKs65UufLZRsZMWkeyzfsYHDPFozs34G61b0AZDLwpOKcKzW27d7H/f9ZzL9mraRF/eq8/Ove/KJtg0SH5Q6CJxXnXKnw0eJ13PrGPH7YupvrT2rNn85pT40q/haVbPy/mHMuoTbt2Mvdby1g8tw1tGtUi4n/dSLdU45IdFjuEHlScc4lhJnxdno2o6YsYMuuffzuzHbcePrRXgAyyXlScc6VuLVbd3PbG/N5f9Faujavy8u/6c2xR3m9rrLAk4pzrsSYGa/OXsV9UxexLzePW/sfy3W/aE2linFdL9CVIE8qzrkSsXLjTkZMSufzpRvp06Y+oy/qSqsGNRMdlosxTyrOubjKzTOe/2w5D76XQeUKFbhvUGcu75lChQpeYqUs8qTinIubjB+CApDfrtrMmcc24t5BnWlS1wtAlmWeVJxzMbc3J4+nZmTy5EeZ1K5WmUcHd+OC47wAZHngScU5F1PfrtrMLRPSyVi7jYHdmvKXAR050gtAlhueVJxzMbFrby4PT8/guU+X06h2NZ4dkspZHRsnOixXwjypOOcO2+dLNzBy0jy+37iTK3qnMKLfsdSp5gUgy6O43hwuqa+kDEmZkkYU0i5Jj4Xt6ZK6R7SNk7RO0vwCY/4uaXHY/w1J9cLXW0naJWlu+DMmnnNzzsHW3fsYOWkeVzzzJQCv/KYPfx3UxRNKORa3IxVJFYEngbOBLGC2pClmtjCiWz+CteTbAb2Bp8PfAC8ATwAvFdj0dGBkuFzx/cBIYHjYttTMusV+Ns65gt5fuJbbJs9j/bY9DD2lDX84qz3Vq3iJlfIunqe/egGZZrYMQNJ4YCDBGvL5BgIvWbCm8SxJ9SQ1MbNsM5spqVXBjZrZexFPZwGXxG0Gzrmf2bh9D3e9tZAp367h2KNqM/aqVI5rUS/RYblSIp5JpRmwKuJ5Fj8dhRTVpxmQHeU+rgNejXjeWtI3wFbgdjP7pOAASUOBoQApKSlR7sY5Z2ZM+XYNo6YsYPueHP54dnuGnXo0VSp5iRX3k3gmlcJuSLdD6FP4xqXbgBzg5fClbCDFzDZK6gFMltTJzLbut3GzscBYgNTU1Kj25Vx5l71lF7e9MZ8PF6+jW4t6PHBJV9o3rp3osFwpFFVSkTQAmGpmeQex7SygRcTz5sCaQ+hTWDxXAwOAM8NTZ5jZHmBP+HiOpKVAeyDtIGJ2zkXIyzNemb2Sv01dTG6ecceAjlxzYisqeokVdwDRHrcOBpZIekBShyjHzAbaSWotqUq4jSkF+kwBhoR3gfUBtphZkae+JPUluDB/gZntjHi9YXhzAJLaEFz8XxZlrM65ApZv2MHlz8zitjfm07V5XabdfArXn9TaE4orUlRHKmb2K0l1gMuB5yUZ8DzwipltO8CYHEk3AdOAisA4M1sgaVjYPgaYCvQHMoGdwLX54yW9ApwGNJCUBdxpZs8R3BFWFZgelnyYZWbDgFOAuyXlALnAMDPbdFB/DeccObl5PPfpch6e/h1VKlXg/ou7cGlqCy+x4qKi8OxRdJ2lBsCvgJuBRUBb4DEzezwu0cVZamqqpaX52THn8i3K3srwiemkZ23h7I6NuffCzjSuUy3RYblSRtIcM0strC3aayoXEBxFHA38E+hlZusk1SBILkmZVJxzgT05uTz50VKe+iiTutUr8/jlxzOgaxM/OnEHLdq7vy4BHjGzmZEvmtlOSdfFPiznXEn5euWPDJ+QzpJ127no+GbcMaAjR9SskuiwXJKKNqlkF0woku43s+Fm9kEc4nLOxdnOvTk8OO07nv98OUfVqcbz1/Tk9GMbJTosl+SiTSpn81MplHz9CnnNOZcEPsvcwIhJ6azatIur+rTklr7HUNvrdbkYKDKpSPov4LfA0ZLSI5pqA5/FMzDnXOxt2bWP+95ZyGtpWbRuUJNXh/ahd5sjEx2WK0OKO1L5N/Au8DcgssrwNr9d17nkMm3BD9wxeT4bd+xl2KlHc/NZ7ahW2QtAutgqLqmYma2QdGPBBkn1PbE4V/qt37aHUVMW8M68bDo0qcNzV/ekS/O6iQ7LlVHRHKkMAOYQ1OSKvL/QgDZxiss5d5jMjDe+Wc3dby9k555c/nzuMQw9pQ2VK3oBSBc/RSYVMxsQ/m5dMuE452Jh9eZd3PbGPGZkrKd7SlAAsm0jLwDp4q+4C/Xdi2o3s69jG45z7nDk5Rkvf/k9o99djAGjzu/IkBNaUcHrdbkSUtzpr4eKaDPgjBjG4pw7DEvXb2fkxHl8tWITJ7drwF8HdaFF/RqJDsuVM8Wd/jq9pAJxzh2anNw8nvlkOY+8/x3VKlXg75d05ZIezb3EikuI4k5/nWFmH0q6qLB2M5sUn7Ccc9FYsGYLwyemM3/1Vvp2Ooq7L+xEo9peANIlTnGnv04FPgTOL6TNAE8qziXA7n25PPFhJmM+Xkq9GlV4+sru9OvSJNFhOVfs6a87w9/XFtXPOVdy0lZsYvjEdJau38HF3Ztzx4AO1KvhBSBd6RBt6fsjgTuBkwiOUD4F7jazjXGMzTkXYceeHP4+LYMXv1hB07rVefG6XpzavmGiw3JuP9F+C2o8sB64mKAM/nrg1eIGSeorKUNSpqQRhbRL0mNhe3rkLcySxklaJ2l+gTH1JU2XtCT8fURE28hwWxmSzo1ybs6VejO/W885j8zkxS9WMKRPS6b94RRPKK5Uijap1Deze8xsefhzL1CvqAHhevFPElQz7ghcLqljgW79CNaSbwcMBZ6OaHsB6FvIpkcAH5hZO+CD8DnhtgcDncJxT+WvWe9cstq8cy//8/q3DBn3FVUrV+D1G07groGdqVU12gLjzpWsaJPKR5IGS6oQ/lwKvFPMmF5AppktM7O9BEc7Awv0GQi8ZIFZQD1JTQDC9VsKqy02EHgxfPwicGHE6+PNbI+ZLSdY975XlPNzrtR5d142Zz08kze+Wc2Npx/N1N+dTGqr+okOy7kiFXdL8TZ+qvn1R+BfYVMFYDvBdZYDaQasinieBfSOok8zILuI7TY2s2wAM8uWlL+qUDNgViHb2o+koQRHRaSkpBSxG+cSY9223dz55gLenf8DnZrW4cXretKpqReAdMmhuLu/DqdYUGHfvLJD6BPL/WFmY4GxAKmpqYe6L+dizsyY+PVq7nl7Ibv2eQFIl5yiPjEbXhBvB/z/N6sKLjFcQBbQIuJ5c2DNIfQpaK2kJuFRShNg3WFsy7lSIevHndz6xnxmfree1JZHMPrirrRtVCvRYTl30KL6CCTp18BMYBpwV/h7VDHDZgPtJLWWVIXgIvqUAn2mAEPCu8D6AFvyT20VYQpwdfj4auDNiNcHS6oqqTVBAvyq2Mk5l0B5ecaLn6/gnEdmMmfFJu66oBOv3XCCJxSXtKI9Uvk90BOYZWanSzqWILkckJnlSLqJIAFVBMaZ2QJJw8L2McBUoD/BRfWdwP9/yVLSK8BpQANJWcCdZvYcMBp4TdL1wErgl+H2Fkh6DVgI5AA3mllulPNzrsRlrtvOiInppH3/I6e0b8hfB3Wm+RFeANIlN5kVf1lB0mwz6ylpLtDbzPZImmtm3eIdYDylpqZaWlpaosNw5cy+3DzGzlzGo+8voXqVitwxoCMXd2/mBSBd0pA0x8xSC2uL9kglS1I9YDIwXdKP+PUK5w7a/NVbuGVCOguzt9K/y1HcdUFnGtaumuiwnIuZqJKKmQ0KH46S9BFQF/hP3KJyrozZvS+XRz9YwtiZy6hfswpjftWDvp2PSnRYzsXcwdz91Z2fan99Fn6h0TlXjNkrNjF8QjrLNuzg0tTm3Na/I3VrVE50WM7FRbQFJf9CcEE8v9T985JeD8u1OOcKsX1PDn//z2JemvU9zepV55/X9+Lkdl6vy5Vt0R6pXA4cb2a7ASSNBr4GPKk4V4gZGeu47Y35rNmyi2tObMX/nHMMNb1elysHov1XvoLgS4+7w+dVgaXxCMi5ZPbjjr3c885CJn29mraNajFh2In0aHlE8QOdKyOKq/31OME1lD3AAknTw+dnE6yp4pwjKLHy7vwf+Mub89m8cx//fUZbbjqjLVUreaFsV74Ud6SS/yWOOcAbEa/PiEs0ziWhdVt3c8eb85m2YC1dmtXlpet607FpnUSH5VxCFFdQMr/EPGGplfbh0wwz2xfPwJwr7cyM19OyuPedhezJyWNkv2O5/qTWVPICkK4ci/bur9MI1i5ZQVANuIWkq4spKOlcmbVq005GTprHp5kb6NW6PqMv6kKbhl6vy7loL9Q/BJxjZhkAktoDrwA94hWYc6VRblgA8u/TMqhYQdx7YWeu6JVChQpeYsU5iD6pVM5PKABm9p0k//aWK1eWrN3GLRPT+WblZk47piF/HdSFpvWqJzos50qVaJPKHEnPAf8Mn19JcPHeuTJvb04eYz5eyhMfZlKzakX+cVk3BnZr6gUgnStEtEllGHAj8DuCayozgafiFZRzpUV61mZumZDO4h+2MaBrE0Zd0IkGtbwApHMHUmxSkVQBmGNmnYGH4x+Sc4m3e18uj0z/jmc+WUbD2lV5ZkgqZ3dsnOiwnCv1ik0qZpYn6VtJKWa2siSCci6RZi3byIiJ6azYuJPLe7VgRL8O1K3ulxCdi0a0p7+aEHyj/itgR/6LZnZBUYMk9QUeJVj58VkzG12gXWF7f4KVH68xs6+LGivpVeCYcBP1gM1m1k1SK2ARkH9DwSwzGxbl/Jxj2+59jH53MS9/uZKU+jX49697c2LbBokOy7mkEm1SKXLp4MJIqgg8SVDSJQuYLWmKmS2M6NaPYC35dkBv4Gmgd1FjzeyyiH08BGyJ2N7SZF+N0iXGh4vXctsb81m7dTfXn9SaP53TnhpVvACkcweruNpf1Qgu0rcF5gHPmVlOlNvuBWSa2bJwW+OBgQRryOcbCLxkwZrGsyTVk9QEaFXc2PAo51LgjCjjce5nNu3Yy91vLWDy3DW0b1yLp648keNTvACkc4equI9iLwL7gE8Ijio6Ar+PctvNgFURz7MIjkaK69MsyrEnA2vNbEnEa60lfQNsBW43s08KBiVpKDAUICUlJcqpuLLGzHgrPZtRUxawbfc+fn9mO248vS1VKnmJFecOR3FJpaOZdQEIv6fy1UFsu7Cb+C3KPtGMvZzgW/35soEUM9soqQcwWVInM9u630bMxgJjAVJTUwtu05UDP2zZze2T5/P+orUc17wu91/Sm2OP8gKQzsVCcUnl/4tGmlnOQX7ZKwtoEfG8ObAmyj5VihorqRJwERFlYsxsD0GJfsxsjqSlBAUw03CO4Ohk/OxV/PWdRezLy+O2/h247qTWVPQSK87FTHFJ5ThJ+Z/0BVQPnwswMyvq491soJ2k1sBqYDBwRYE+U4CbwmsmvYEtZpYtaX0xY88CFptZVv4LkhoCm8wsV1Ibgov/y4qZnysnvt+4gxET5/HFso2c0OZIRl/chZZH1kx0WM6VOcWVvj/kFYbCI5ubgGkEtwWPM7MFkoaF7WOAqQS3E2cS3FJ8bVFjIzY/mP1PfQGcAtwtKQfIBYaZ2aZDjd+VDbl5xvOfLefB9zKoXKECf7uoC4N7tvASK87FiYIbr8qn1NRUS0vzs2NlVcYPQQHIb1dt5qwOjbj3wi4cVbdaosNyLulJmmNmqYW1+Y34rszZm5PHUzMyefKjTGpXq8xjlx/P+V2b+NGJcyXAk4orU+au2szwCelkrN3GwG5NufP8TtSvWSXRYTlXbnhScWXCrr25PPReBuM+W06j2tV47upUzuzgBSCdK2meVFzS+3zpBkZMnMfKTTu5oncKI/sdS+1qXgDSuUTwpOKS1pZd+/jb1EWMn72KVkfWYPzQPvRpc2Siw3KuXPOk4pLS9IVruX3yPNZv28MNp7Th5rPaU73KId8B75yLEU8qLqls3L6HUW8t5K1v13DsUbV5ZkgqXZvXS3RYzrmQJxWXFMyMN+eu4a63FrB9Tw5/OKs9/3Xa0V4A0rlSxpOKK/XWbN7F7ZPn8+HidRyfUo/7L+5K+8a1Ex2Wc64QnlRcqZWXZ/z7q5WMfncxuXnGXwZ05OoTW3kBSOdKMU8qrlRavmEHwyem89XyTZzUtgF/u6gLLerXSHRYzrlieFJxpUpObh7PfrqcR6Z/R5VKFbj/4i5cmuoFIJ1LFp5UXKmxKHsrwyemk561hXM6NuaeCzvTuI4XgHQumXhScQm3JyeXJz/M5KkZS6lXozJPXtGd/l2O8qMT55KQJxWXUHO+/5HhE9PJXLedi7o3447zOnKEF4B0Lml5UnEJsXNvDn+flsELn6+gad3qvHBtT047plGiw3LOHaa4fnNMUl9JGZIyJY0opF2SHgvb0yV1L26spFGSVkuaG/70j2gbGfbPkHRuPOfmDt2nSzZwziMzef6zFVzVpyXT/nCKJxTnyoi4HalIqgg8CZwNZAGzJU0xs4UR3foRrCXfjmCN+qeB3lGMfcTMHiywv44Eywx3ApoC70tqb2a58ZqjOzhbdu3jvncW8lpaFq0b1OS1G06gV+v6iQ7LORdD8Tz91QvINLNlAJLGAwOByKQyEHjJgjWNZ0mqJ6kJ0CqKsQUNBMab2R5guaTMMIYvYjstdyimLfiBOybPZ+OOvQw79WhuPqsd1Sp7AUjnypp4JpVmwKqI51kERyPF9WkWxdibJA0B0oA/mdmP4ZhZhWxrP5KGAkMBUlJSDmI67lCs37aHUW8t4J30bDo0qcNzV/ekS/O6iQ7LORcn8bymUtj9oBZln6LGPg0cDXQDsoGHDmJ/mNlYM0s1s9SGDRsWMsTFgpkx6esszn7kY6YvWMufzz2GKTf9whOKc2VcPI9UsoAWEc+bA2ui7FPlQGPNbG3+i5KeAd4+iP25ErB68y5unTSPj79bT4+WR3D/xV1o28gLQDpXHsTzSGU20E5Sa0lVCC6iTynQZwowJLwLrA+wxcyyixobXnPJNwiYH7GtwZKqSmpNcPH/q3hNzv1cXp7x0hcrOOfhj5m9YhOjzu/I6zec4AnFuXIkbkcqZpYj6SZgGlARGGdmCyQNC9vHAFOB/kAmsBO4tqix4aYfkNSN4NTWCuCGcMwCSa8RXMzPAW70O79KztL12xkxMZ3ZK37k5HYN+OsgLwDpXHmk4Mar8ik1NdXS0tISHUZSy8nN45lPlvPI+99RvXJFbj+vA5f0aO4lVpwrwyTNMbPUwtr8G/XukC1Ys4XhE9OZv3orfTsdxd0XdqJRbS8A6Vx55knFHbTd+3J5/MMljPl4GUfUqMJTV3anf5cmxQ90zpV5nlTcQZnz/SZumZDO0vU7uLh7c+4Y0IF6NbwApHMu4EnFRWXHnqAA5ItfBAUgX7yuF6e29+/5OOf250nFFWvmd+sZOWkea7bs4uoTWvHnc4+hZlX/p+Oc+zl/Z3AHtHnnXu55exETv86iTcOavH7DCaS28gKQzrkD86TiCvXuvGzueHMBP+7cy42nH81/n+EFIJ1zxfOk4vazbutu/vLmAv6z4Ac6Na3Di9f1pFNTr9flnIuOJxUHBAUgJ8zJ4p63F7I7J49b+h7D0JPbUKliXNdxc86VMZ5UHKs27eTWN+bxyZIN9Gx1BKMv7srRDWslOiznXBLypFKO5ReAfGBaBgLuGdiJK3u3pEIFL7HinDs0nlTKqcx12xk+MZ053//Iqe0bct+gzjQ/wgtAOucOjyeVcmZfbh5jZy7j0feXUKNqRR6+9DgGHd/MC0A652LCk0o5Mn/1Fv48IZ1F2Vs5r0sTRl3QiYa1qyY6LOdcGeJJpRzYvS+XRz9YwtiZy6hfswpjftWDvp2PSnRYzrkyyJNKGffV8k2MmJjOsg07uDS1Obf170jdGpUTHZZzroyK65cQJPWVlCEpU9KIQtol6bGwPV1S9+LGSvq7pMVh/zck1QtfbyVpl6S54c+YeM6ttNu+J4c7Js/n0v/9gr25efzr+t48cMlxnlCcc3EVtyMVSRWBJ4GzgSxgtqQpZrYwols/grXk2wG9gaeB3sWMnQ6MDJccvh8YCQwPt7fUzLrFa07JYkbGOm6dNI/srbu57het+Z9z21Ojih+UOufiL57vNL2ATDNbBiBpPDCQYA35fAOBlyxY03iWpHqSmgCtDjTWzN6LGD8LuCSOc0gqP+7Yyz1vL2TSN6tp26gWE4adSI+WRyQ6LOdcORLPpNIMWBXxPIvgaKS4Ps2iHAtwHfBqxPPWkr4BtgK3m9knBQdIGgoMBUhJSYlqIqWdmTF13g/cOWU+m3fu43dntOXGM9pStZIXgHTOlax4JpXCvvhgUfYpdqyk24Ac4OXwpWwgxcw2SuoBTJbUycy27rcRs7HAWIDU1NSC8SSdtVt3c8fk+by3cC1dmtXln9f3pkOTOokOyzlXTsUzqWQBLSKeNwfWRNmnSlFjJV0NDADODE+dYWZ7gD3h4zmSlgLtgbRYTKa0MTNeS1vFve8sYm9OHiP7Hcv1J7X2ApDOuYSKZ1KZDbST1BpYDQwGrijQZwpwU3jNpDewxcyyJa0/0FhJfQkuzJ9qZjvzNySpIbDJzHIltSG4+L8sjvNLmJUbdzLyjXQ+y9xIr9b1uf/irrRuUDPRYTnnXPySSnh31k3ANKAiMM7MFkgaFraPAaYC/YFMYCdwbVFjw00/AVQFpoelRWaZ2TDgFOBuSTlALjDMzDbFa36JkJtnvPD5Ch6clkHFCuK+QZ25vGeKF4B0zpUaCs8elUupqamWlpYcZ8e+W7uNWyakM3fVZs44thH3DepMk7rVEx2Wc64ckjTHzFILa/MvL5Rye3PyGPPxUh7/cAm1qlbi0cHduOC4pl4A0jlXKnlSKcXSszZzy4R0Fv+wjfOPa8qo8ztyZC0vAOmcK708qZRCu/bm8o/3v+OZT5bRsHZVnhmSytkdGyc6LOecK5YnlVJm1rKNjJiYzoqNO7m8VwtG9u9AnWper8s5lxw8qZQS23bvY/S7i3n5y5Wk1K/Bv3/dmxPbNkh0WM45d1A8qZQCHy5ey21vzGft1t38+qTW/OmcY6hexUusOOeSjyeVBNq4fQ93v72QN+euoX3jWjx15Ykcn+IFIJ1zycuTSgKYGVO+XcNdby1k2+59/P7Mdtx4eluqVPISK8655OZJpYRlb9nFHZPn8/6idRzXoh4PXNyVY46qneiwnHMuJjyplJC8PGP87FX8beoi9uXlcft5Hbj2F62p6CVWnHNliCeVErBiww5GTEpn1rJNnNDmSEZf3IWWR3oBSOdc2eNJJY5y84xxny7noekZVK5Qgb8O6sLlvVp4iRXnXJnlSSVOMn7Yxi0T0/l21WbO6tCYey/szFF1qyU6LOeciytPKjG2NyePp2Zk8uRHmdSpVpnHLz+eAV2b+NGJc65c8KQSQ3NXbWb4hHQy1m5j0PHNuGNAR+rXrJLosJxzrsR4UomBXXtzeXh6Bs99upzGdaox7ppUzjjWC0A658qfuH7bTlJfSRmSMiWNKKRdkh4L29MldS9urKT6kqZLWhL+PiKibWTYP0PSufGcW77Pl26g76MzeeaT5VzeK4X3/nCKJxTnXLkVt6QiqSLwJNAP6AhcLqljgW79CNaSbwcMBZ6OYuwI4AMzawd8ED4nbB8MdAL6Ak+F24mLLbv2MXJSOlc88yUCxg/tw32DulDbKwo758qxeJ7+6gVkmtkyAEnjgYHAwog+A4GXLFjTeJakepKaAK2KGDsQOC0c/yIwAxgevj7ezPYAyyVlhjF8EeuJpWdt5jcvpbF+2x5uOKUNN5/V3gtAOucc8U0qzYBVEc+zgN5R9GlWzNjGZpYNYGbZkhpFbGtWIdvaj6ShBEdFpKSkHMR0fpJSvwbtG9fmmSGpdG1e75C24ZxzZVE8k0ph99BalH2iGXso+8PMxgJjAVJTU4vbZqHq1ajCP68vmB+dc87F80J9FtAi4nlzYE2UfYoauzY8RUb4e91B7M8551wcxTOpzAbaSWotqQrBRfQpBfpMAYaEd4H1AbaEp7aKGjsFuDp8fDXwZsTrgyVVldSa4OL/V/GanHPOuZ+L2+kvM8uRdBMwDagIjDOzBZKGhe1jgKlAfyAT2AlcW9TYcNOjgdckXQ+sBH4Zjlkg6TWCi/k5wI1mlhuv+TnnnPs5BTdelU+pqamWlpaW6DCccy6pSJpjZqmFtflSg84552LGk4pzzrmY8aTinHMuZjypOOeci5lyfaFe0nrg+8PYRANgQ4zCSQblbb7gcy4vfM4Hp6WZNSysoVwnlcMlKe1Ad0CUReVtvuBzLi98zrHjp7+cc87FjCcV55xzMeNJ5fCMTXQAJay8zRd8zuWFzzlG/JqKc865mPEjFeecczHjScU551zMeFIphqS+kjIkZUoaUUi7JD0WtqdL6p6IOGMpijlfGc41XdLnko5LRJyxVNycI/r1lJQr6ZKSjC8eopmzpNMkzZW0QNLHJR1jrEXxb7uupLckfRvO+dpExBkrksZJWidp/gHaY//+ZWb+c4AfgrL7S4E2QBXgW6BjgT79gXcJVp7sA3yZ6LhLYM4nAkeEj/uVhzlH9PuQYMmGSxIddwn8d65HsJRESvi8UaLjLoE53wrcHz5uCGwCqiQ69sOY8ylAd2D+Adpj/v7lRypF6wVkmtkyM9sLjAcGFugzEHjJArOAevkrUyapYudsZp+b2Y/h01kEq2wms2j+OwP8NzCRn1YbTWbRzPkKYJKZrQQws2SfdzRzNqC2JAG1CJJKTsmGGTtmNpNgDgcS8/cvTypFawasinieFb52sH2SycHO53qCTzrJrNg5S2oGDALGlGBc8RTNf+f2wBGSZkiaI2lIiUUXH9HM+QmgA8FS5POA35tZXsmElxAxf/+K28qPZYQKea3gPdjR9EkmUc9H0ukESeWkuEYUf9HM+R/AcDPLDT7EJr1o5lwJ6AGcCVQHvpA0y8y+i3dwcRLNnM8F5gJnAEcD0yV9YmZb4xxbosT8/cuTStGygBYRz5sTfII52D7JJKr5SOoKPAv0M7ONJRRbvEQz51RgfJhQGgD9JeWY2eQSiTD2ov23vcHMdgA7JM0EjgOSNalEM+drgdEWXHDIlLQcOBb4qmRCLHExf//y019Fmw20k9RaUhVgMDClQJ8pwJDwLoo+wBYzyy7pQGOo2DlLSgEmAVcl8afWSMXO2cxam1krM2sFTAB+m8QJBaL7t/0mcLKkSpJqAL2BRSUcZyxFM+eVBEdmSGoMHAMsK9EoS1bM37/8SKUIZpYj6SZgGsGdI+PMbIGkYWH7GII7gfoDmcBOgk86SSvKOf8FOBJ4KvzknmNJXOE1yjmXKdHM2cwWSfoPkA7kAc+aWaG3piaDKP873wO8IGkewamh4WaWtCXxJb0CnAY0kJQF3AlUhvi9f3mZFuecczHjp7+cc87FjCcV55xzMeNJxTnnXMx4UnHOORcznlScc87FjCcVV2qF1YDnRvy0OsztdZPUP+L5BUVVJI4FSb+TtEjSy3Hcx1RJ9cLH28PfrQ5UmTaeJL1QWAVnSc9K6ljS8biS599TcaXZLjPrVlhDWPBPB1mXqRvBN+OnApjZFH7+5bdY+y1B1YHl8dqBmfUvvldimdmvEx2DKxl+pOKSRvjpe5Gkp4CvgRaSnpaUFq59cVdE354K1nr5VtJXkuoCdwOXhUc9l0m6RtITYf+Wkj4I15T4IKwakP/J+7FwW8sK+xQe9vujpPnhz83ha2MIyqxPkfSHAv2/lNQp4vkMST0k9Qr39U34+5iw/RpJkyT9R9ISSQ9EjF0hqUExf7dPJH0d/px4gH5Dwvl/K+mfxfxdCn29wPbuCf9+FcL5pYavnyPpizCW1yXVCl8fLWlhuM0HDzQfV8olut6///jPgX6AXILifnOBN4BWBN/s7hPRp374uyIwA+hKsFbGMqBn2FaH4Kj8GuCJiLH//xx4C7g6fHwdMDl8/ALwOsEHsI4EpdMLxtmDoKJtTYJy6QuA48O2FUCDQsb8AbgrfNwE+C4y1vDxWcDEiFiXAXWBasD3QIuC+wC2h79bEa6hAdQAqoWP2wFphcTTCciI2E79Yv4uRf29LgEeAP6Xn75gPYPgKLEBMBOoGb4+nKBCQ/1w//n96yX635//HNqPn/5ypdl+p7/CayrfW7DuQ75LJQ0lSBpNCN74Dcg2s9kAFlaYVdHVhU8ALgof/5PgTTHfZAtOsy1UUA+qoJOANywovIikScDJwDdF7O81YDpB2YxLCRIXBEnjRUntwnlUjhjzgZltCfexEGjJ/mXLD6Qy8ISkbgSJun0hfc4AJlhYksTM8tfgONDfpai/1x0Eiz0NLWQ/fQj+G30W/veoAnwBbAV2A89Kegd4O4p5uVLIk4pLNjvyH0hqDfwPwRHJj5JeIPgULw5/+YHI8XsiHheWmQ66Fr6ZrZa0UUG158uAG8Kme4CPzGxQmERnHCCOXKL///cPwFqCCsMVCN68C4r2b3agPpGvzwZ6SKofkZwi9zPdzC7/WQBSL4JijoOBmwgSnUsyfk3FJbM6BElmS3gE0S98fTHQVFJPAEm1JVUCtgG1D7CtzwnezACuBD49iDhmAhdKqiGpJsFiXp9EMW48cAtQ18zmha/VBVaHj685iBiKUpfgyC0PuIrgVGFBHxAc9R0JIKl++PqB/i5F/b3+A4wG3pFU8O89C/iFpLbhfmpIah9eV6lrZlOBmwluqnBJyI9UXNIys28lfUNwDWMZ8Fn4+l5JlwGPS6oO7CK4PvERMELSXOBvBTb3O2CcpD8D6zmIaq1m9nV4lJS/5sazZlbUqa98E4BHCY5O8j1AcPrrj8CH0cZQjKeAiZJ+SfA32FGwgwXVeu8DPpaUS3Dq7hoO/Hcp8u9lZq+HCWWKIm7jNrP1kq4BXpFUNXz5doKE/6ak/CPN/W5scMnDqxQ755yLGT/95ZxzLmY8qTjnnIsZTyrOOedixpOKc865mPGk4pxzLmY8qTjnnIsZTyrOOedi5v8AJeJ16Yvka4cAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior.plot()\n", "\n", "plt.xlabel('Fraction of vanilla cookies')\n", "plt.ylabel('Probability')\n", "plt.title('Posterior, one vanilla');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we put the first cookie back, stir the bowl, draw from the same bowl again, and get a vanilla cookie again.\n", "\n", "What's are the posterior probabilities now?\n", "\n", "We can do another update, using the posterior from the first draw as the prior for the second draw." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "prior2 = posterior\n", "unnorm2 = prior2 * likelihood_vanilla\n", "posterior2 = unnorm2 / unnorm2.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's what it looks like." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAykUlEQVR4nO3dd3xW5f3/8debsPfGsDcKKisMba2oVRG1aOvABVItWouj1bb+utTafmttrbtaqghqBXGCSrUUZxGEsEFWREaYYYURQtbn98c5aW9jSO5A7tx3ks/z8bgf9xnXOedz3YT7c5/rnHNdMjOcc865aNWIdwDOOecqF08czjnnysQTh3POuTLxxOGcc65MPHE455wrE08czjnnysQTh3NHIekaSf+KdxyJSFJHSQclJYXzH0q6MZy+XtJ/4huhiyVPHC4hSNog6XD4ZbRD0nOSGh7H/u6V9OLxxGRm/zCz845nH0cjqbMkk1QzFvuPNTPbZGYNzSw/3rG4iueJwyWSi82sITAAGAT8Kl6BHM8XugL+f8tVWf7H7RKOmW0B/gmcDCDpO5JWStoXNomcVFhW0s8lbZF0QNIaSedIGg78ArgyPINZGpZtIulZSdvCbX4X0dRyvaQ5kh6WtAe4t2iTi6TTJS2QlBm+nx6x7kNJv5c0B8gCupZSzY/D931hjKdJ2ihpYLi/a8Mzkt7h/I2S3gyn60h6RNLW8PWIpDpFDxCW2yfp5IhlrcIzu9aSmkl6W1KGpL3hdPsidbo//FwOSPqXpJbhuqjPmCQ9KmmzpP2SFko6o7RtXGLzxOESjqQOwAhgsaSewBTgDqAVMBN4S1JtSb2A8cAgM2sEnA9sMLN3gf8DXg6bU/qGu54M5AHdgf7AecCNEYceAqwHWgO/LxJTc+Ad4DGgBfAX4B1JLSKKXQeMAxoBG0up5rfC96ZhjHOBj4BhEevXA2dGzH8UTv8SGAr0A/oCgynm7MzMjgCvA1dFLL4C+MjMdhL8/38O6AR0BA4DTxTZzdXAWILPpDZwVyn1Ks6CMNbmwEvAK5LqHsN+XILwxOESyZuS9gH/IfiS/D/gSuAdM5tlZrnAn4F6wOlAPlAH6C2plpltMLMvituxpDbABcAdZnYo/OJ8GBgVUWyrmT1uZnlmdrjILi4E1pnZC+H6KcBq4OKIMpPMbGW4PvcY6v8R/0sUZwB/iJg/k/8ljmuA35rZTjPLAO4jSFrFeYmvJo6rw2WY2W4ze83MsszsAEGyPLPI9s+Z2drw85hGkADKxMxeDI+VZ2YPEfyb9Srrflzi8MThEsklZtbUzDqZ2S3hl1VbIn69m1kBsBloZ2ZpBGci9wI7JU2V1PYo++4E1AK2hc03+4C/EfySLrS5hNi+EkdoI9Auyu2j8RFwhqQTgCTgZeAbkjoDTYAlR4llY7isOO8D9SQNkdSJ4Iv/DQBJ9SX9LWwi20/QfNa0sPkutD1iOgso8w0Lku6UtCps4tsX1qVlWffjEocnDpfothJ86QPBhWegA7AFwMxeMrNvhmUM+GNYtGi3z5uBI0DLMDk1NbPGZtYnokxJXUV/JY5Qx8I4oti+qK+VDRNhFnAb8HF4FrCdoPnrP2HSLC6WjuGyrx8k2GYawVnH1cDb4X4B7iT45T/EzBrzv+YzlaEeJQqvZ/ycoImsmZk1BTLL8xiu4nnicIluGnBheNG7FsGX3RHgU0m9JJ0dXhjOJmijL7w9dAfQufDuJjPbBvwLeEhSY0k1JHWTVLRp5mhmAj0lXS2ppqQrgd7A20fbQMEtwR8eZXUGUMDXL6J/RHDdprBZ6sMi8xBc8/lVeKG7JfAboKRbj18iaPK7Jpwu1IjgM9sXXsO5p4R9HKtGBNeVMoCakn4DNI7BcVwF8sThEpqZrQGuBR4HdhFcU7jYzHII2sofCJdvJ2h2+kW46Svh+25Ji8Lp0QQXeD8H9gKvAslRxrEbuIggce0GfgZcZGa7StisAzDnKPvLIrimMCdsOhsarvqI4Mv246PMA/wOSAWWAcuBReGyo8X+GXCIoDnrnxGrHiG4XrQLmAe8W0JdjtV74THXEjSpZXP8TXouzuQDOTkXG5KWAOeESce5KsMTh3POuTLxpirnnHNl4onDOedcmXjicM45VyaVsmfOsmrZsqV17tw53mE451ylsnDhwl1m1qro8mqRODp37kxqamq8w3DOuUpFUrF9rnlTlXPOuTLxxOGcc65MYpo4JA0Px0hIk3R3Mesl6bFw/TJJA8LldSXNl7RUwTgM90Vs01zSLEnrwvdmsayDc865r4pZ4gh72HySoCvr3sBVhYPSRLgA6BG+xgFPhcuPAGeH4yj0A4ZHdMlwNzDbzHoAs8N555xzFSSWZxyDgTQzWx/2KzQVGFmkzEjgeQvMI+jSOTmcPxiWqRW+LGKbyeH0ZOCSGNbBOedcEbFMHO34amdm6Xx17IISy0hKCvv62QnMCjtqA2gT9nRa2ONp5HgK/yVpnKRUSakZGRnHWxfnnHOhWCaO4vrbL9ox1lHLmFm+mfUD2gODI8dNjoaZTTCzFDNLadXqa7chO+ecO0axTBzpBN1KF2rP1webKbWMme0jGJNgeLhoh6RkgPB9Z7lF7JxzVUR2bj73zljJnkM55b7vWCaOBUAPSV0k1SYY23lGkTIzgNHh3VVDgUwz2xYOUNMUQFI94NsE4zsXbjMmnB4DTI9hHZxzrtIxM3715gomz93Aii2Z5b7/mD05bmZ5ksYTDOSSBEw0s5WSbg7XP00wqtoIoHDIzLHh5snA5PDOrBrANDMrHGntAWCapBuATcDlsaqDc85VRi8v2MyrC9O57ZwefKtn+TfVV4vxOFJSUsy7HHHOVQfL0zP53tOfMqRLcyaNHUxSjWMf3l3SQjNLKbrcnxx3zrkqYl9WDj/8x0JaNqjNo6P6H1fSKEm16OTQOeequoIC4/apS9ixP5tpN51G8wa1Y3YsP+Nwzrkq4NHZ6/hobQb3XNyH/h1j2xOTJw7nnKvk3l+9g0dnr+N7A9pzzZCOMT+eJw7nnKvENu3O4o6pS+id3JjfX3oyUmyua0TyxOGcc5XU4Zx8xr2QiiSevnYgdWslVchx/eK4c85VQmbG3a8vY82OAzx3/SA6tqhfYcf2Mw7nnKuEJn26gelLtnLnuT0Z1qvYvl5jxhOHc85VMp+t383v31nFub3bcMuw7hV+fE8czjlXiWzdd5hb/rGIjs3r89AVfakRo4f8SuKJwznnKons3Hx++OJCjuQVMGF0Co3r1opLHH5x3DnnKgEz49dvrmBpeiZ/u24g3Vs3jFssfsbhnHOVwIvzNvLKwnRuO7s75/c5Ia6xeOJwzrkE99n63dz31uecfWJr7vh2z3iH44nDOecS2ZbCi+Et6vPIqH5xuRhelCcO55xLUIdz8rnphVRy8gr4exwvhhflF8edcy4BFT4ZvnLrfp4ZnUK3VvG7GF6Un3E451wCmvDxeqYv2cpd5/XinJPaxDucr/DE4ZxzCeaDNTt54N3VXHhqMrcM6xbvcL7GE4dzziWQ9RkHuW3KYk46oTF/uuzUCukmvaw8cTjnXILIPJzLjZNTqZVUgwmjB1K/dmJehvbE4ZxzCSAvv4Bbpyxm894snr52IO2bVVw36WWVmOnMOeeqmT/8czUfr83gge+ewuAuzeMdTon8jMM55+JsWupmnv3Pl1x/emdGDY79mOHHyxOHc87F0fwv9/DLN5ZzRo+W/OrCk+IdTlRimjgkDZe0RlKapLuLWS9Jj4Xrl0kaEC7vIOkDSaskrZR0e8Q290raImlJ+BoRyzo451ysbN6Txc0vLqRDs/o8cdUAaiZVjt/yMbvGISkJeBI4F0gHFkiaYWafRxS7AOgRvoYAT4XvecCdZrZIUiNgoaRZEds+bGZ/jlXszjkXaweP5HHj5FTy8gt4ZkwKTeonRnci0YhlehsMpJnZejPLAaYCI4uUGQk8b4F5QFNJyWa2zcwWAZjZAWAV0C6GsTrnXIXJLzBun7KYtIyD/PWagXRNoO5EohHLxNEO2Bwxn87Xv/xLLSOpM9Af+Cxi8fiwaWuipGbFHVzSOEmpklIzMjKOsQrOOVf+/jBzFbNX7+Te7/Thmz1axjucMotl4ijucUcrSxlJDYHXgDvMbH+4+CmgG9AP2AY8VNzBzWyCmaWYWUqrVq3KGLpzzsXGlPmbeCa8g+q6oZ3iHc4xiWXiSAc6RMy3B7ZGW0ZSLYKk8Q8ze72wgJntMLN8MysA/k7QJOaccwnv07Rd/PrNFZzZs1WluYOqOLFMHAuAHpK6SKoNjAJmFCkzAxgd3l01FMg0s20KOmd5FlhlZn+J3EBScsTspcCK2FXBOefKxxcZB7n5xYV0admAx6/uX2nuoCpOzO6qMrM8SeOB94AkYKKZrZR0c7j+aWAmMAJIA7KAseHm3wCuA5ZLWhIu+4WZzQQelNSPoElrA3BTrOrgnHPlYc+hHL4/aQG1a9Zg4vWDEmZApmMls6KXHaqelJQUS01NjXcYzrlq6EhePtc+8xlL0zOZOm4oAzoWez9PQpK00MxSii6vvOdKzjmX4MyMu19bzoINe3no8r6VKmmUxBOHc87FyKOz1/HG4i3ceW5PLu7bNt7hlBtPHM45FwOvL0rnkX+v47KB7Rl/dvd4h1OuPHE451w5m7d+Nz9/bRmndW3B/116SkKO4nc8PHE451w5Stt5kJteWEjH5vV5+tqB1K5Z9b5mq16NnHMuTjIOHOH65+ZTK0lMGju4UnVcWBY+AqBzzpWDwzn53Ph8KrsOHuHlcafRoXniDv16vDxxOOfcccovMG6buphl6fv427UD6duhabxDiilvqnLOueNgZtz31kpmfb6Dey7qzXl9Toh3SDHnicM5547DhI/X8/zcjfzgjC5c/40u8Q6nQnjicM65YzRj6Vb+8M/VXHhqMv/vgsrb221ZeeJwzrljMPeL3dw1bSmDOzfnocv7UqNG1XpWoySeOJxzroxWb9/PuBdS6dSiPn8fnULdWknxDqlCeeJwzrky2LrvMNdPXED92klM+n7VfVajJJ44nHMuSplZuVz/3HwOHclj0tjBtGtaL94hxYU/x+Gcc1HIzs3nxucXsGFXFpO+P4iTkhvHO6S48cThnHOlyMsv4NYpi0nduJfHr+rP6d1axjukuPKmKuecK4GZ8evp/3vA76JTq864GsfKE4dzzpXg4VlrmTJ/E7cM61ZtHvArjScO55w7ismfbuCx99O4IqU9Pz2/V7zDSRieOJxzrhhvL9vKvW+t5NzebarkYEzHwxOHc84V8cm6DH788hJSOjXj8av6UzPJvyoj+afhnHMRFm/ay00vLKRbq4Y8M2ZQtXsqPBqeOJxzLrR2xwHGTlpAq0Z1eP6GwTSpV/2eCo9GTBOHpOGS1khKk3R3Mesl6bFw/TJJA8LlHSR9IGmVpJWSbo/YprmkWZLWhe/NYlkH51z1sHlPFqOfnU+tpBq88P0htG5UN94hJayYJQ5JScCTwAVAb+AqSb2LFLsA6BG+xgFPhcvzgDvN7CRgKPCjiG3vBmabWQ9gdjjvnHPHbOeBbK599jMO5+bzwg2D6dii6g77Wh5iecYxGEgzs/VmlgNMBUYWKTMSeN4C84CmkpLNbJuZLQIwswPAKqBdxDaTw+nJwCUxrINzrorLzMpl9LPz2bn/CM+NHcSJJ1TfrkSiFcvE0Q7YHDGfzv++/KMuI6kz0B/4LFzUxsy2AYTvrYs7uKRxklIlpWZkZBxrHZxzVdihI3mMnTSf9RmHmDB6IAM6est3NGKZOIq76dnKUkZSQ+A14A4z21+Wg5vZBDNLMbOUVq1alWVT51w1kJ2bz7gXUlmyeR+PXdWPM3r490S0Ypk40oEOEfPtga3RlpFUiyBp/MPMXo8os0NSclgmGdhZznE756q43LDTwjlpu/nTZX0ZfnJyvEOqVGKZOBYAPSR1kVQbGAXMKFJmBjA6vLtqKJBpZtsUPKL5LLDKzP5SzDZjwukxwPTYVcE5V9UUFBg/fWUpsz7fwW9H9uF7A9vHO6RKJ2bdqptZnqTxwHtAEjDRzFZKujlc/zQwExgBpAFZwNhw828A1wHLJS0Jl/3CzGYCDwDTJN0AbAIuj1UdnHNVi5nxyzdX8OaSrfz0/F6MPq1zvEOqlGRW9LJD1ZOSkmKpqanxDsM5F0dmxv1vr2LinC8Zf1Z37vJOC0slaaGZpRRd7k+OO+eqhYf+tZaJc75k7Dc6c+d5PeMdTqXmicM5V+U9PnsdT3yQxlWDO/Cbi3p7T7fHyROHc65Km/DxFzw0ay3f7d+O313i3aOXB08czrkqa/KnG/i/mau58NRkHrzsVJJqeNIoD544nHNV0kufbeKeGcFATI9c2c/H1ChHUX2Ski6S5J+6c65SmJa6mV+8sZyzerXiiav7U8uTRrmK9tMcBayT9KCkk2IZkHPOHY/XF6Xz89eWcUaPljx17UDq1PSBmMpbVInDzK4l6GjwC+A5SXPDTgQbxTQ655wrg+lLtnDXK0s5rWsL/j46xUfvi5Goz9/CTgZfI+gePRm4FFgk6dYYxeacc1GbsXQrP355CYO7NOeZMZ40YinaaxzfkfQG8D5QCxhsZhcAfYG7Yhifc86V6q2lW7lj6mIGdW7OxOsHUb92zHpTckTfV9VlwMNm9nHkQjPLkvT98g/LOeei89bSrdzx8hJSOjfnubGeNCpCtE1V24omDUl/BDCz2eUelXPORWH6ki3cPnUxAzs24zk/06gw0SaOc4tZdkF5BuKcc2UxfckWfhxxptGgjieNilLiJy3ph8AtQDdJyyJWNQLmxDIw55w7mjcWp3PntKUM7uLXNOKhtE/7JeCfwB+AuyOWHzCzPTGLyjnnjuKV1M387LVlnNa1Bc+MSfGkEQelfeJmZhsk/ajoCknNPXk45yrSlPmb+H+vL+eMHi2ZcF0K9Wr7LbfxEM0Zx0XAQsCAyB7CDOgao7icc+4rXpi7gV9PX8mwXq14+tqB/pxGHJWYOMzsovC9S8WE45xzX/fMJ+v53Tur+PZJbXjymv7ejUiclXZxfEBJ681sUfmG45xzX/XkB2n86b01XHhKMo+M6ucdFiaA0pqqHiphnQFnl2Mszjn3X2bGX2at5fH307i0fzv+dNmp3jV6giitqeqsigrEOecKmRn3v72KiXO+ZNSgDvz+0lN8EKYEUlpT1dlm9r6k7xa33sxej01YzrnqKr/A+OUby5m6YDNjv9HZxwhPQKU1VZ1J0LHhxcWsM8ATh3Ou3OTmF3DntKXMWLqVH53VjbvO6+VJIwGV1lR1T/g+tmLCcc5VV9m5+Yx/aRH/XrWTnw3vxS3Dusc7JHcU0Xar3kLSY5IWSVoo6VFJLaLYbrikNZLSJN1dzHqF+02TtCzyLi5JEyXtlLSiyDb3StoiaUn4GhFNHZxzievgkTzGPreAf6/ayf0j+3jSSHDR3qIwFcgAvkfQxXoG8HJJG0hKAp4k6AyxN3CVpN5Fil0A9Ahf44CnItZNAoYfZfcPm1m/8DUzyjo45xLQ3kM5XPP3eczfsIe/XNGX607rHO+QXCmiTRzNzex+M/syfP0OaFrKNoOBNDNbb2Y5BMlnZJEyI4HnLTAPaCopGSDsxt27NHGuCtuemc0Vf5vLqu0HePragXx3QPt4h+SiEG3i+EDSKEk1wtcVwDulbNMO2Bwxnx4uK2uZ4owPm7YmSmoWRXnnXILZsOsQlz39Kdsys5k8djDn9m4T75BclEpMHJIOSNoP3ETQb1VO+JoK/LiUfRd3K4QdQ5mingK6Af2AbRzlIUVJ4ySlSkrNyMgoZZfOuYq0Yksmlz39KVk5+bz0gyGc1q3US6YugZSYOMyskZk1Dt9rmFnN8FXDzBqXsu90oEPEfHtg6zGUKRrTDjPLN7MC4O8ETWLFlZtgZilmltKqVatSQnXOVZS5X+xm1IR51KmZxCs3n8ap7ZvGOyRXRlE/vy+pmaTBkr5V+CplkwVAD0ldJNUGRgEzipSZAYwO764aCmSa2bZS4kiOmL0UWHG0ss65xPLuiu2MmTif5CZ1efWHp9GtVcN4h+SOQVQjoEi6Ebid4IxgCTAUmEsJfVWZWZ6k8cB7QBIw0cxWSro5XP80MBMYAaQBWcB/nxeRNAUYBrSUlA7cY2bPAg9K6kfQpLWBoBnNOZfgXpy3kd9MX0HfDk2ZOGYQzRrUjndI7hjJrLRLCiBpOTAImGdm/SSdCNxnZlfGOsDykJKSYqmpqfEOw7lqycx45N/reHT2Os4+sTVPXN3fR+2rJCQtNLOUosuj/dfLNrNsSUiqY2arJfUq5xidc1VMXn4Bv56+kinzN3HZwPb84buneLfoVUC0iSNdUlPgTWCWpL2UchHbOVe9Hc7J59YpQRcitwzrxk/P936nqoqoEoeZXRpO3ivpA6AJ8G7MonLOVWp7DuXw/UkLWJq+j/tH9vGnwauYqBsaw36kvklwUXpO+DS4c859xYZdh7j+uflsy8zmqWsGMvzkE+Idkitn0XZy+BtgMtACaAk8J+lXsQzMOVf5LN60l+8+9SmZh3N56QdDPGlUUdGecVwF9DezbABJDwCLgN/FKjDnXOXy3srt3D51MW0a12XS2MF0adkg3iG5GIk2cWwA6gLZ4Xwd4ItYBOScq1zMjIlzNvC7dz6nb/umPDMmhZYN68Q7LBdDpQ0d+zjBNY0jwEpJs8L5c4H/xD4851wiy8sv4P63P2fy3I0M73MCD1/Zj3q1k+Idloux0s44Cp+aWwi8EbH8w5hE45yrNA4eyeP2KYuZvXon477VlbuHn0iNGn67bXVQ2tCxkwunw/6meoaza8wsN5aBOecS19Z9h7lhciprdxzg/ktO5rqhneIdkqtA0fZVNYzgrqoNBF2hd5A0JhxsyTlXjSxPz+SGyQs4nJPPxOsHcWZP7326uon24vhDwHlmtgZAUk9gCjAwVoE55xLPzOXb+Mm0JbRoUIcXfjiEXic0indILg6iTRy1CpMGgJmtlVQrRjE55xKMmfHE+2k8NGstAzo2ZcJov3OqOos2cSyU9CzwQjh/DcEFc+dcFZedm8/PXl3GjKVbuaRfWx743qnUreV3TlVn0SaOm4EfAbcRXOP4GPhrrIJyziWG7ZnZjHshlWXpmdx1Xk9+dFZ376jQlZ44JNUAFprZycBfYh+Scy4RLN28jx88n8rBI3lMuG4g5/Xx7kNcoNS+qsKxvZdK6lgB8TjnEsDri9K5/G9zqV2zBq/fcronDfcV0TZVJRM8OT4fOFS40My+E5OonHNxkZdfwB/fXc3fP/mSoV2b89drBtLch3h1RUSbOO6LaRTOubjbl5XDrVMW88m6XYw5rRO/uqi3j9bnilVaX1V1CS6MdweWA8+aWV5FBOacqzirtu1n3Aup7Mg8wh+/dwpXDvKWaXd0pZ1xTAZygU+AC4DewO2xDso5V3FmLN3Kz19dRuN6NZl601AGdGwW75BcgistcfQ2s1MAwuc45sc+JOdcRcjNL+CBf67m2f98SUqnZvz12gG0blQ33mG5SqC0xPHfjgzNLM/v33auasg4cITxLy3isy/3MOa0Tvzywt7UrunXM1x0SkscfSXtD6cF1AvnBZiZNY5pdM65crdw4x5+9I/F7Ducw8NX9uXS/u3jHZKrZErrVt37FXCuijAzJn26gd+/s4q2Tevx2g9Pp0/bJvEOy1VCMT03lTRc0hpJaZLuLma9JD0Wrl8maUDEuomSdkpaUWSb5pJmSVoXvvuVPOdKcfBIHrdOWcx9b33OsF6teevWb3rScMcsZolDUhLwJP+7G+sqSb2LFLsA6BG+xgFPRaybBAwvZtd3A7PNrAcwO5x3zh3F6u37+c7j/2Hm8m38bHgvJlw3kCb1vHNrd+xiecYxGEgzs/VmlgNMBUYWKTMSeN4C84CmkpIBwkGi9hSz35EEtwkTvl8Si+CdqwpeSd3MJU/O4cCRPF76wVBuGdbdh3d1xy3aJ8ePRTtgc8R8OjAkijLtgG0l7LeNmW0DMLNtkloXV0jSOIKzGDp29IeZXPWSlZPHr99cyWuL0jm9WwseHdWfVo18/AxXPmKZOIr7WWPHUOaYmNkEYAJASkpKuezTucpgzfYD/OilRXyRcZDbzunB7ef0IMnPMlw5imXiSAc6RMy3B7YeQ5midkhKDs82koGdxx2pc1WAmfHygs3c+9ZKGtapxYs3DOEb3VvGOyxXBcXyGscCoIekLpJqA6OAGUXKzABGh3dXDQUyC5uhSjADGBNOjwGml2fQzlVG+7NzGT9lMXe/vpyUTs2Zefs3PWm4mInZGUf4pPl44D0gCZhoZisl3RyufxqYCYwA0oAsYGzh9pKmAMOAlpLSgXvM7FngAWCapBuATcDlsaqDc5XBok17uX3qYrbuy+Znw3tx87e6+QVwF1Myq/rN/ykpKZaamhrvMJwrV/kFxlMfpvHwv9dxQuO6PHZVfwZ28seaXPmRtNDMUoouj+U1DudcjGzdd5ifTFvCvPV7uLhvW35/6ck0ruvPZriK4YnDuUrmraVb+eUby8kvMP58eV++N6Ad3gGpq0ieOJyrJPZn53Lv9JW8vngL/To05dFR/ejUokG8w3LVkCcO5yqBeet3c+e0pWzLPMxt5/Tg1rO7+7CuLm48cTiXwI7k5fOXf61lwifr6dS8Pq/cfLpfAHdx54nDuQS1YksmP5m2hLU7DnL1kI78csRJNKjj/2Vd/PlfoXMJJje/gKc+/ILHZq+jeYPaPDd2EGf1KrZLNufiwhOHcwlkzfYD3PnKElZs2c/Ifm257zt9aFq/drzDcu4rPHE4lwDy8gv428freeTfa2lctxZPXTOAC05JjndYzhXLE4dzcbZq235+9uoylm/J5MJTk/ntd/rQoqF3ge4SlycO5+IkJ6+AJz9I48kP0mhavxZ/vWYAI/wsw1UCnjici4NFm/Zy92vLWLvjIJf0a8s9F/ehWQO/luEqB08czlWgQ0fy+PO/1jDp0w2c0Lguz45J4ZyT2sQ7LOfKxBOHcxVk9qod/Gb6SrZmHua6oZ342fATaejPZbhKyP9qnYuxHfuzue+tlcxcvp2ebRry6s2nMbBT83iH5dwx88ThXIzkFxgvzN3An/+1lpz8An56fi9+cEZXatf0PqZc5eaJw7kYWLp5H796cwXLt2RyRo+W3D/yZDq39J5sXdXgicO5crQvK4cH31vDlPmbaNmwDo9f1Z+LTk328TJcleKJw7lyUFBgvLJwM398dw2Zh3MZe3oXfnxuDxr5qHyuCvLE4dxxWrJ5H/dMX8HS9ExSOjXjtyNPpnfbxvEOy7mY8cTh3DHaeSCbP727hlcWptO6UR0eubIfI/u19WYpV+V54nCujI7k5fPcnA08PnsdOfkF3PStrtx6Tg9/JsNVG/6X7lyUzIz3Vm7nD/9czcbdWXz7pNb88sLedPG7pVw144nDuSgsT8/k/nc+Z/6Xe+jZpiHPf38w3+rZKt5hORcXMX0SSdJwSWskpUm6u5j1kvRYuH6ZpAGlbSvpXklbJC0JXyNiWQdXvW3ek8XtUxdz8RP/4YudB/n9pScz87YzPGm4ai1mZxySkoAngXOBdGCBpBlm9nlEsQuAHuFrCPAUMCSKbR82sz/HKnbn9mXl8NSHX/DcpxsQcMuwbtw8rBuN/fZa52LaVDUYSDOz9QCSpgIjgcjEMRJ43swMmCepqaRkoHMU2zpX7rJz85n06Qb++kEaB47kcWn/dtx1Xi/aNq0X79CcSxixTBztgM0R8+kEZxWllWkXxbbjJY0GUoE7zWxv0YNLGgeMA+jYseMxVsFVF7n5Bby6MJ1H/72O7fuzOfvE1vz0/F6clOzPYzhXVCwTR3E3s1uUZUra9ing/nD+fuAh4PtfK2w2AZgAkJKSUvS4zgHBE99vL9/Gw7PW8uWuQwzo2JRHRvVjaNcW8Q7NuYQVy8SRDnSImG8PbI2yTO2jbWtmOwoXSvo78Hb5heyqi8Jbax+etY41Ow5w4gmNeHZMCmef2Nof4HOuFLFMHAuAHpK6AFuAUcDVRcrMIGh2mkrQFJVpZtskZRxtW0nJZrYt3P5SYEUM6+CqGDPj36t28si/17Jy6366tmrA41f158JTkqlRwxOGc9GIWeIwszxJ44H3gCRgopmtlHRzuP5pYCYwAkgDsoCxJW0b7vpBSf0Imqo2ADfFqg6u6jAzZn2+g0dnr2Pl1v10alGfP1/el0v6taVmko+P4VxZKLihqWpLSUmx1NTUeIfh4qCgwPjniu08/v46Vm8/QKcW9bn17B6eMJyLgqSFZpZSdLk/Oe6qpNz8AqYv2crTH31B2s6DdG3VgIcu78tITxjOHTdPHK5KycrJY9qCzfz9ky/Zsu8wJ57QiMev6s+IU5JJ8msYzpULTxyuSthzKIfJn27g+bkb2JuVy8BOzbj/kj6c1cvvknKuvHnicJXa+oyDTJzzJa8uTCc7t4Bvn9SGm8/sSkrn5vEOzbkqyxOHq3TMjLlf7GbinA3MXr2DWkk1uLRfO37wrS50b90o3uE5V+V54nCVRnZuPtOXbOG5ORtYvf0AzRvU5taze3Dd0E60alQn3uE5V2144nAJb/OeLF6ct5GXUzezLyuXE09oxIOXncp3+ralbq2keIfnXLXjicMlpPwC48M1O3lx3kY+XJtBDYnz+7ThuqGdGdq1uV/wdi6OPHG4hLI9M5tpqZt5ecFmtuw7TKtGdRh/VneuHtKR5CbetblzicATh4u73PwCPli9k2mpm3l/9U4KDE7v1oJfjDiJ8/q0oZY/sOdcQvHE4eJm7Y4DvLYwndcWbWHXwSO0bFiHm87sxqhBHejUokG8w3POHYUnDlehdh88wltLt/Laoi0s35JJzRri7BNbc0VKB87s1crPLpyrBDxxuJjLyslj1uc7eHPxFj5et4v8AqN3cmN+fVFvRvZrS8uGfiutc5WJJw4XE9m5+Xy8NoMZS7cye9VODufmk9ykLj84oyuX9G/LiSf4kKzOVVaeOFy5yc7N56O1Gcxcvo3Zq3Zy8EgezerX4rsD2nFx37YM7tzcB0tyrgrwxOGOS+bhXD5cs5N3V2zno7UZZOXk07R+LS48JZkRpyZzercWft3CuSrGE4crs817spi9agezVu3gs/V7yCswWjWqw6X92zH85BMY2tWThXNVmScOV6ojefks3LCXD9dm8P7qnaTtPAhA99YNufGMrpzbuzX9OzTzZijnqglPHO5rzIwvMg4xJ20XH6/NYO763WTl5FMrSQzp0oKrBnfk7BNb06WlP2vhXHXkicMBQfPT3PW7mbd+N5+m7Wb7/mwAOrWoz/cGtOfMnq0Y2q0FDev4n4xz1Z1/C1RDhWcUCzbsYf6XwWvLvsMAtGhQm6FdW/DNHi35ZveWdGheP87ROucSjSeOamB/di7L0zNZsnkfCzfuZdGmvezLygWgZcPaDO7SnB+c0YXTurWkZ5uG3vOsc65EnjiqmAPZuazefoDl6Zks35LJsvR9rN91CLNgfbdWDTi/9wkM7NSMgZ2b0bVlA08Uzrky8cRRSeXlF7BpTxZrth9g9fYDrNl+gFXb97Nxd9Z/y7RpXIdT2jVhZL929OvQlL7tm9Kkfq04Ru2cqwo8cSQwM2PPoRw27sliw65DfLnrEOszDpG28yBf7jpETn4BABJ0al6fPm0bc/nA9vRu25iT2zahdeO6ca6Bc64qimnikDQceBRIAp4xsweKrFe4fgSQBVxvZotK2lZSc+BloDOwAbjCzPbGsh6xkpdfwK6DOWzLPMz2zGy2ZmazZe9h0vdmkb73MJv3ZHHgSN5/yyfVEO2b1aN7q4YMO7EV3Vs1pNcJjejRuhH1avsQqs65ihGzxCEpCXgSOBdIBxZImmFmn0cUuwDoEb6GAE8BQ0rZ9m5gtpk9IOnucP7nsapHacyMnPwCDufkk5WTz6EjeezPzuNAdi77s/PIzMphX1Yue7Ny2XPoCLsP5bDrYA4ZB7LZfSjnv9ceCtWvnUS7pvVo16wegzo3o2OLBnRqXp8urRrQoVl9atf0J7Kdc/EVyzOOwUCama0HkDQVGAlEJo6RwPNmZsA8SU0lJROcTRxt25HAsHD7ycCHxChxPDZ7HTOWbqXADDMoMCM3r4DcAiM3v4AjuQUcycunwErfV4PaSTRvWJsWDerQtkld+nVoQqtGdWndqA7JTeqS3KQeyU3q0rR+Lb9Y7ZxLaLFMHO2AzRHz6QRnFaWVaVfKtm3MbBuAmW2T1Lq4g0saB4wD6Nix4zFVoHWjOvRq0wgJJFFDUCupBrWSRK2kGtSpWYO6tZKoU7MG9WvXpEGdJOrXrkmjuoWvWjStX4sm9WpRp6Y3JTnnqoZYJo7ifjYX/W1+tDLRbFsiM5sATABISUkp07aFRg3uyKjBx5Z0nHOuqoplg3k60CFivj2wNcoyJW27I2zOInzfWY4xO+ecK0UsE8cCoIekLpJqA6OAGUXKzABGKzAUyAyboUradgYwJpweA0yPYR2cc84VEbOmKjPLkzQeeI/gltqJZrZS0s3h+qeBmQS34qYR3I47tqRtw10/AEyTdAOwCbg8VnVwzjn3dbKi94NWQSkpKZaamhrvMJxzrlKRtNDMUoou94cCnHPOlYknDuecc2XiicM551yZeOJwzjlXJtXi4rikDGDjMW7eEthVjuFUBl7n6sHrXD0cT507mVmrogurReI4HpJSi7uroCrzOlcPXufqIRZ19qYq55xzZeKJwznnXJl44ijdhHgHEAde5+rB61w9lHud/RqHc865MvEzDuecc2XiicM551yZeOIISRouaY2ktHAs86LrJemxcP0ySQPiEWd5iqLO14R1XSbpU0l94xFneSqtzhHlBknKl3RZRcZX3qKpr6RhkpZIWinpo4qOsbxF8XfdRNJbkpaGdR4bjzjLk6SJknZKWnGU9eX7/WVm1f5F0HX7F0BXoDawFOhdpMwI4J8EoxMOBT6Ld9wVUOfTgWbh9AXVoc4R5d4n6Pb/snjHHeN/46bA50DHcL51vOOugDr/AvhjON0K2APUjnfsx1nvbwEDgBVHWV+u319+xhEYDKSZ2XozywGmAiOLlBkJPG+BeUDTwpEIK6lS62xmn5rZ3nB2HsFIjJVZNP/OALcCr1H5R5eMpr5XA6+b2SYAM6sOdTagkSQBDQkSR17Fhlm+zOxjgnocTbl+f3niCLQDNkfMp4fLylqmMilrfW4g+MVSmZVaZ0ntgEuBpyswrliJ5t+4J9BM0oeSFkoaXWHRxUY0dX4COIlgOOrlwO1mVlAx4cVNuX5/xWwEwEpGxSwrep9yNGUqk6jrI+ksgsTxzZhGFHvR1PkR4Odmlh/8IK3UoqlvTWAgcA5QD5graZ6ZrY11cDESTZ3PB5YAZwPdgFmSPjGz/TGOLZ7K9fvLE0cgHegQMd+e4NdIWctUJlHVR9KpwDPABWa2u4Jii5Vo6pwCTA2TRktghKQ8M3uzQiIsX9H+Xe8ys0PAIUkfA32Bypo4oqnzWOABCxr/0yR9CZwIzK+YEOOiXL+/vKkqsADoIamLpNrAKGBGkTIzgNHh3QlDgUwz21bRgZajUussqSPwOnBdJf4FGqnUOptZFzPrbGadgVeBWypp0oDo/q6nA2dIqimpPjAEWFXBcZanaOq8ieAMC0ltgF7A+gqNsuKV6/eXn3EAZpYnaTzwHsFdGRPNbKWkm8P1TxPcYTMCSAOyCH61VFpR1vk3QAvgr+Ev8DyrxD2LRlnnKiOa+prZKknvAsuAAuAZMyv2ls7KIMp/4/uBSZKWEzTh/NzMKnVX65KmAMOAlpLSgXuAWhCb7y/vcsQ551yZeFOVc865MvHE4Zxzrkw8cTjnnCsTTxzOOefKxBOHc865MvHE4eIu7IV2ScSr83Hur5+kERHz3ympJ9zyIOk2Sask/SOGx5gpqWk4fTB873y0HlFjSdKk4noOlvSMpN4VHY+rWP4ch0sEh82sX3Erwo7oVMa+hPoRPAE+E8DMZvD1h8DK2y0ET9d/GasDmNmI0kvFl5ndGO8YXOz5GYdLOOGv6FWS/gosAjpIekpSajh+wn0RZQcpGCtkqaT5kpoAvwWuDM9erpR0vaQnwvKdJM0OxySYHT4dX/gL+rFwX+uL+zUdlvuJpBXh645w2dME3XjPkPTjIuU/k9QnYv5DSQMlDQ6PtTh87xWuv17S65LelbRO0oMR226Q1LKUz+0TSYvC1+lHKTc6rP9SSS+U8rkUu7zI/u4PP78aYf1SwuXnSZobxvKKpIbh8gckfR7u889Hq49LYPHuR95f/gLyCTqdWwK8AXQmeIp5aESZ5uF7EvAhcCrBeAvrgUHhusYEZ9HXA09EbPvfeeAtYEw4/X3gzXB6EvAKwY+p3gRdcxeNcyBBb6oNCLrjXgn0D9dtAFoWs82PgfvC6WRgbWSs4fS3gdciYl0PNAHqAhuBDkWPARwM3zsTjsEA1AfqhtM9gNRi4ukDrInYT/NSPpeSPq/LgAeBv/G/h4k/JDjbawl8DDQIl/+coCeC5uHxC8s3jfffn7/K/vKmKpcIvtJUFV7j2GjBuAGFrpA0jiAxJBN8uRuwzcwWAFjYu6lK7tX2NOC74fQLBF98hd60oEnscwV9GBX1TeANCzoERNLrwBnA4hKONw2YRdAFxBUEyQmCxDBZUo+wHrUitpltZpnhMT4HOvHVLrGPphbwhKR+BMm4ZzFlzgZetbCLDTMrHMPhaJ9LSZ/XrwkGBBpXzHGGEvwbzQn/PWoDc4H9QDbwjKR3gLejqJdLMJ44XKI6VDghqQtwF8GZxV5Jkwh+jYvj79o+cvsjEdPFZZ8y97NuZlsk7VbQy/CVwE3hqvuBD8zs0jBRfniUOPKJ/v/pj4EdBL3b1iD4gi4q2s/saGUily8ABkpqHpGAIo8zy8yu+loA0mCCTgZHAeMJkpmrRPwah6sMGhMkkszwTOCCcPlqoK2kQQCSGkmqCRwAGh1lX58SfGEBXAP8pwxxfAxcIqm+pAYEAz59EsV2U4GfAU3MbHm4rAmwJZy+vgwxlKQJwRlYAXAdQbNeUbMJzt5aAEhqHi4/2udS0uf1LvAA8I6kop/3POAbkrqHx6kvqWd4naOJmc0E7iC4kcFVMn7G4RKemS2VtJjgmsJ6YE64PEfSlcDjkuoBhwmuF3wA3C1pCfCHIru7DZgo6adABmXoJdTMFoVnO4XjNjxjZiU1UxV6FXiU4Cyj0IMETVU/IRjfvDz8FXhN0uUEn8GhogUs6Cn298BHkvIJmtmu5+ifS4mfl5m9EiaNGYq4BdrMMiRdD0yRVCdc/CuCpD5dUuEZ41duJnCVg/eO65xzrky8qco551yZeOJwzjlXJp44nHPOlYknDuecc2XiicM551yZeOJwzjlXJp44nHPOlcn/B0m7+ierT+ekAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior2.plot()\n", "\n", "plt.xlabel('Fraction of vanilla cookies')\n", "plt.ylabel('Probability')\n", "plt.title('Posterior, two vanilla');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2\n", "\n", "Suppose we put the second cookie back, stir the bowl, draw from the same bowl again, and get a chocolate cookie.\n", "\n", "Now what are the posterior probabilities?" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "prior3 = posterior2\n", "unnorm3 = prior3 * likelihood_chocolate\n", "posterior3 = unnorm3 / unnorm3.sum()" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABBW0lEQVR4nO3dd3hUZfbA8e9JJYGQkEIvCb1JDUVB7IpYsCLWRWVZrKuuq+6urrqrv3Vd3bXLWlhlVRA7KooVKy2hd0IPBAi9pZDk/P64N+sYUiY4k5tJzud57pOZW887M5kz933vfV9RVYwxxphACPM6AGOMMXWHJRVjjDEBY0nFGGNMwFhSMcYYEzCWVIwxxgSMJRVjjDEBY0nFeEpErhSRz7yOozYSkbYiclBEwt3nM0VkrPt4jIh8722EgSEiJ4tIdg0f8xUReagmj1lfWFKpR0Rkg4jkuV9U20XkPyLS6Bfs7wERee2XxKSqr6vqmb9kHxURkVQRURGJCMb+g01VN6lqI1Ut9jqW+sw3mZuqWVKpf85T1UZAP2AAcK9XgfySL3tx2OfXmFrG/inrKVXdAnwC9AQQkfNFZJmI7HV/mXUrXVdE7haRLSJyQERWichpIjIc+CNwmXvms8hdN15EXhaRHHebh3yqb8aIyA8i8i8R2Q08ULYaR0ROEJF5IrLP/XuCz7KZIvKwiPwAHAbaV1HMb92/e90YjxeRjSLS393fVe6ZTHf3+VgRed99HC0iT4jIVnd6QkSiyx7AXW+viPT0mZfinhE2FZEmIvKRiOSKyB73cesyZfqr+7ocEJHPRCTZXeb3mZaIPCkim0Vkv4hkisiJVW3js203N4697mfgfJ9lr4jIsyLysRvfHBHp4LO8q4h8LiK73c/GqEqOk+ieHW91X4v3yyz/nYjscD871/rMjxeRSe5ruFFE7vX9QSEivxaRFW58y0WkX1XlKnPcCt8jEXkYOBF4xv0MPVPdctc7qmpTPZmADcDp7uM2wDLgr0Bn4BBwBhAJ3AVkAVFAF2Az0NLdLhXo4D5+AHitzDHeB/4NNASaAnOB37jLxgBFwC1ABBDjzvveXZ4I7AGudpdf7j5PcpfPBDYBPdzlkVWUNxVQIMJn3iTgd+7jF4C1wA0+y253H/8FmO2WIQX4EfhrBceZCDzs8/wm4FP3cRJwMRALxAFvAe/7rDvTjaGz+3rMBB4pL3532Vif1/J7n/1c5R4rAvgdsA1o4MdnItJ9r//ovt+nAgeALu7yV4DdwEB3368DU9xlDd3PxrXusn7ATqBHBcf6GHgTaOIe9yR3/snu5+Iv7vwROD8amvi8Lx+4r18qsBq43l12KbAF56xbgI5AOz/L9VA13qOxPs+rVe76NnkegE01+GY7SeUgsBfYCDznfpHdB0z1WS/M/Uc92f0n3QGcTpkvccokFaAZUADE+My7HPjafTwG2FRmH2P4KalcDcwts3wWMMZ9PBP4SzXKm8rRSeV6YJr7eAUwlp++JDcC/dzHa4ERPtudBWyo4DinA+t8nv8AXFPBun2APT7PZwL3+jy/kZ8S0s/ip5KkUs5x9gC9/XiNTsRJQGE+8yYDD7iPXwFe8lk2AljpPr4M+K7M/v4N3F/OcVoAJbiJosyyk4G8Mu/TDmAwEO5+prr7LPsNMNN9PAP47TGW66FqvEe+ScXvctfHKSQbMM0vcoGqfuE7Q0Ra4nyhAqCqJSKyGWilqjNF5DacBNJDRGYAd6jq1nL2XfoLMUdESueF4fyqK7W57EY+fhaHayPQys/t/fEN8JiINMf5wnoTuF9EUoF4YGEFsWx055XnKyBGRAbhfJH1Ad4DEJFY4F/AcJxf6ABxIhKuPzXAb/PZ12Gg2hdPiMjvcBJkS5xE1BhI9mPTlsBmVS3xmVf2Na8ovnbAIBHZ67M8AvhvOcdpA+xW1T0VxLFLVYvKOU4yzplG2feiNL42OD8AyvKnXIDf75Gv6pS73rE2FQOwFecfBXAawXH+WbcAqOobqjrUXUeBv7urlu3iejPOr8pkVU1wp8aq2sNnncq6xf5ZHK62pXH4sX1ZR62rqlk4X1i3At+q6gGcL81xOL/8S7+EysbS1p139EGcbabinJVdAXzk7hecqqguwCBVbQwMc+fLUTs6Rm77yd3AKJwzgQRgn5/H2Aq0kZ9f9FD2Na/IZuAbn/c6QZ2r1W6oYN1EEUnwY7++dgJHOPq9KI1vM9Ch7EZUr1xVvUflfc79LXe9Y0nFgPOFeI44DfCROP9kBcCPItJFRE51G6nzcaopSn+9bQdSS/9xVTUH+Ax4XEQai0iYiHQQkZP8jGM60FlErhCRCBG5DOgOfFTRBuJc1jyzgsW5OFUuZRv0vwFudv+CU73h+xycqpJ73Ub3ZODPQGWXT7+BUy1ypfu4VBzOa7ZXRBKB+yvZx7GKw2mTyAUiROTPOGcqwP/uA6koGc/BaU+7S0QiReRk4Dxgih/H/Qjn/bra3TZSRAaIz0UepdzPxifAc27DeKSIDDtqj0dvV4zz+XxYROJEpB1wBz+9Fy8Bd4pIf3F0dNepTrmqeo+28/PPkN/lro8sqRhUdRVOQ+/TOL8Mz8O59LgQiAYecedvw2m4/qO76Vvu310iMt99fA1OdcVynHr9t3Hq0/2JYxdwLk5S24VzwcC5qrqzks3a4LRhlLe/w8DDwA/uFUCD3UXf4HyRfFvBc4CHgAxgMbAEmO/Oqyj20i+xljhfnqWewGm32onT8P9pJWU5VjPcY67GqeLJ5+fVhG1w2qaO4r7H5wNnuzE+h9MetLKqg7pnY2cCo3HODLbhnMUedZWc62qcs46VOG0mt1V1DNctOK/tOuB7nKQ90Y3hLZz3+A2chvj3gcRqlusJKn+PngQuca8Me+oYyl2viNvIZExIEpGFwGluQjLlEJGXgLdUdYbXsZi6z5KKMcaYgLHqL2OMMQFjScUYY0zAWFIxxhgTMPX65sfk5GRNTU31OgxjjAkpmZmZO1U1pbxl9TqppKamkpGR4XUYxhgTUkSkbM8X/2PVX8YYYwLGkooxxpiAsaRijDEmYCypGGOMCRhLKsYYYwLGkooxxpiAsaRijDEmYOr1fSrGmMDJP1LMxl2H2b4/n+3789l1qJCi4hKKShRViI+JpEnDSBIbRtM+uSGtEmIICwvYWGWmlghqUhGR4ThjEYTjjHP9SJnl4i4fgTMa3xhVne8um4gztsYOVe3ps82bOKO0ASQAe1W1jzsc7ApglbtstqqOD1LRjKnXVJWNuw7zfdZO5m3YzYqc/azNPURxif+9njeIDKNT0zgGpCZyQockBrZPpHGDyCBGbWpC0JKKiIQDzwJnANnAPBGZpqrLfVY7G+jkToOA592/AK8AzwCTfPerqpf5HONxnGFTS61V1T4BLYgxBoCSEmXB5j18uCiHz5dvZ8vePACaxkXTs1U8Z3ZvTufmcbSIb0CzuAYkNYoiKiKMcHHORvbnH2HP4SPkHihgXe5BVm8/yIqc/bw+ZyMTf1hPeJgwpGMyF/ZtyZndm9Mw2ipSQlEw37WBQJaqrgMQkSnASJwRAUuNBCapM6jLbBFJEJEWqpqjqt+6Zx/lcs9yRgGnBq0Exhi27cvnjTkbeWf+FrbszSMqIoyTOqcw/qT2DOmYTFpyQ0SqrsZKiI0iITaKtOSGDExL/N/8/CPFLNi0l29W5/Lhoq3c/uYiYiKXclG/Vlw/NI32KY2CWTwTYMFMKq34+ZCm2fx0FlLZOq2AHD/2fyKwXVXX+MxLE5EFwH7gXlX9ruxGIjIOGAfQtm1bPw5jTP2UuXEPL323js+Wb6dElWGdUvjdmZ05o3sz4gJYTdUgMpzjOyRxfIck7jqrC5mb9vB2RjZvZWbzxtxNnNGtGbee1omereIDdkwTPMFMKuX9dClb4erPOhW5HJjs8zwHaKuqu0SkP/C+iPRQ1f0/27nqC8ALAOnp6TbspTFlzFm3i6e/yuL7rJ3Ex0QydmgaVw5qR9uk2KAfOyxMGJCayIDURO48qwuTZm1g0qyNnPfM91zUtzV3ntWZFvExQY/DHLtgJpVsoI3P89bA1mNY5ygiEgFcBPQvnaeqBUCB+zhTRNYCnQHrhtgYP6zI2c//TV/Bd2t2ktwomj+N6MaVg9sSG+VN20ZKXDS/O7MLY09sz3NfZ/GfHzbw8ZKt3HJqJ34zrD0R4XZHRG0UzE/LPKCTiKQBW4DRwBVl1pkG3Oy2twwC9qmqP1VfpwMrVTW7dIaIpAC7VbVYRNrjNP6vC0A5jKnTdhzI5/EZq3krczNxDSK595xuXDmoHTFR4V6HBjiXIv9hRDeuGtyOhz9ewT9mrGLGsm3845LedGke53V4poygJRVVLRKRm4EZOJcUT1TVZSIy3l0+AZiOczlxFs4lxdeWbi8ik4GTgWQRyQbuV9WX3cWj+XnVF8Aw4C8iUgQUA+NVdXewymdMqCspUd6Yu4m/f7qS/CPFXDskjVtO7UhCbJTXoZWrTWIsE67uz8eLc/jzB0s59+nv+P1ZXfj1ie39ulDA1AxxLryqn9LT09UG6TL10ZrtB7jn3SVkbtzD8e2TeOjCnnQIoausdh0s4E/vLeXTZds4s3sz/nFpb+Jj7B6XmiIimaqaXt4yq5Q0ph4pKVFe+m4d5zz9PetyD/LYpb1549eDQiqhACQ1iub5q/px37nd+WrlDs57+nuWb91f9YYm6CypGFNPbNmbx5UvzeGhj1cwrFMKn99xEpf0bx2yVUciwvVD03jzN4MpKCrm0gk/8u3qXK/DqvcsqRhTD3yxfDsjnvyOxdl7efTiXrx4TX+SG0V7HVZA9G+XyAc3DaVNYizXvTKPqRmbq97IBI0lFWPqsCPFJfxt+grGTsqgdZMYPr71REYNaBOyZycVaR7fgLfGH+/cQPn2Yp79OsvrkOot61zHmDpq58ECbnxtPnM37OaqwW2595zuNIisHZcJB0Ncg0gmjhnA799axD9mrOJIcQm3nd7Z67DqHUsqxtRBS7fsY9ykDHYfLuTJ0X0Y2aeV1yHViMjwMB4f1YeI8DCe+GINJQq3n96pzp2Z1WaWVIypYz5avJU731pEk9go3h5/Qr3rMys8THj04l6ECTz15RoEuP0MO2OpKZZUjKkjVJXnv1nLo5+uIr1dE56/qj8pcXWjMb66wsKERy7qRYnCk1+uIalRFNccn+p1WPWCJRVj6oCi4hLu+2AZk+du4vzeLfnHpb2Ijqi77Sf+cBLLcew9XMj905aR1DCac3q18DqsOs+u/jImxOUVFvPrSRlMnruJG0/uwBOX9an3CaVURHgYT1/ej/5tm3D7mwv5ce1Or0Oq8yypGBPC9uUd4eqX5zBzdS4PXdCTu4Z3tXHfy4iJCuelX6XTLimWG16bz/qdh7wOqU6zpGJMiNpxIJ/RL8xmUfZenrm8H1cNbud1SLVWQmwUE8cMIEzg15MyOJB/xOuQ6ixLKsaEoK178xg1YRYbdh7i5V8NsLYCP7RJjOXZK/uxfuchfjtlIcUl9bcz3WCypGJMiNm8+zCXvTCLXQcLeW3sIIZ1TvE6pJBxQodk7j/P6YTy8c9WeR1OnWRJxZgQsnHXIUa/MJt9h4/w+q8H0b9dE69DCjlXD27HZelteG7mWr5etcPrcOocSyrGhIiNuw5x2b9nc6iwiDd+PZherRO8DikkiQgPjuxB1+Zx/G7qIrbty/c6pDrFkooxISB7z2GueHEO+UXFvDF2cL27Sz7QGkSG88wV/cg/UsytUxZQVFzidUh1hiUVY2q5nH15XPHiHA7kH+G16wfRvWVjr0OqEzo2bcRDF/Rk7vrdPPXlGq/DqTMsqRhTi+08WMCVL85h96FCJl0/yM5QAuyifq25uF9rnvk6i8yNe7wOp04IalIRkeEiskpEskTknnKWi4g85S5fLCL9fJZNFJEdIrK0zDYPiMgWEVnoTiN8lv3B3dcqETkrmGUzJtj25x/hmpfnsnVfHv+5dgB92iR4HVKd9MD53WkRH8Odby3icGGR1+GEvKAlFREJB54Fzga6A5eLSPcyq50NdHKnccDzPsteAYZXsPt/qWofd5ruHq87MBro4W73nBuDMSEnr7CYsa9ksGbHASZc1Z8BqYleh1RnxTWI5B+X9mL9zkM8+qldZvxLBfNMZSCQparrVLUQmAKMLLPOSGCSOmYDCSLSAkBVvwV2V+N4I4EpqlqgquuBLDcGY0LKkeISbnw9k3kbd/Ovy/pwcpemXodU553QIZkxJ6Tyyo8b+CHL+gf7JYKZVFoBvoNFZ7vzqrtOeW52q8smikjphfp+7UtExolIhohk5Obm+nEoY2qOqnLPO0v4elUuD19wHOf2aul1SPXG3cO70j65IXe9vZhDBVYNdqyCmVTK69WubL8I/qxT1vNAB6APkAM8Xp19qeoLqpququkpKXYnsqldHvtsFe/Mz+a20ztxxaC2XodTr8REhfP3S3qxZW8eT3yx2utwQlYwk0o20MbneWtg6zGs8zOqul1Vi1W1BHiRn6q4qr0vY2qTV3/cwLNfr+XygW347WmdvA6nXhqQmsjlA9sw8YcNLN2yz+twQlIwk8o8oJOIpIlIFE4j+rQy60wDrnGvAhsM7FPVnMp2Wtrm4roQKL06bBowWkSiRSQNp/F/biAKYkywfbZsGw98uIzTuzXjryN72pjqHrpneDeaxEbxx/eWWKeTxyBoSUVVi4CbgRnACmCqqi4TkfEiMt5dbTqwDqdR/UXgxtLtRWQyMAvoIiLZInK9u+hREVkiIouBU4Db3eMtA6YCy4FPgZtUtThY5TMmUBZn7+W3UxbSq1U8T1/el4hwu33MS/Gxkfz5vO4szt7Hqz9u8DqckCOq9TcTp6ena0ZGhtdhmHose89hLnj2R6Ijwnj/piH1dkz52kZVGfOfeWRs2M3Xvz+ZpnENvA6pVhGRTFVNL2+Z/SQyxiP7849w3SvzKCgq5pVrB1hCqUVEhAfO70FhcYndu1JNllSM8UBxiXLr5AWsyz3EhKv606lZnNchmTLSkhty3dA03s7MZuHmvV6HEzIsqRjjgf+bvoKZq3J5cGQPhnRM9jocU4FbTu1ESlw0D0xbRok12vvFkooxNWzy3E28/P16rh2SypWDbFz52qxRdAR3D+/Kws17eXfBFq/DCQmWVIypQXPW7eK+95dyUucU/jSim9fhGD9c1LcVvdsk8PdPV1qHk36wpGJMDdmyN48bX59P26RYnr7CLh0OFWFhwp/P7UbugQJe/m691+HUevapNqYG5BUWM25SBoVFJbx4TTqNG0R6HZKphv7tEjmzezP+/e06dh0s8DqcWs2SijFBpqrc9c5ilufs56nL+9IhpZHXIZljcNfwLhwuLOKZr7O8DqVWs6RiTJC9+N06Ply0ld+f1YVTulo39qGqY9M4RqW34bXZG9m067DX4dRallSMCaIfs3byyCcrGXFcc244qYPX4Zhf6LbTOxMeJjz2md0QWRFLKsYEyZa9edw8eQEdUhrx6CW9rZPIOqB5fAOuHZLGtEVbWbltv9fh1EqWVIwJgvwjxdzwWiZHikqYcHV/GkVHeB2SCZDfDGtPo+gInvxijdeh1EqWVIwJggc/XMbi7H08Pqq3NczXMQmxUVw7JJVPlm5j+VY7WynLkooxAfZ2ZjaT527mhpM7cGaP5l6HY4Jg7ND2xEVH8OSXNkJkWZZUjAmgFTn7+dN7Szi+fRK/O6Oz1+GYIImPjeS6oWnMWLadZVtthEhfllSMCZD9+Ue44bVM4mMiecoG26rzrhuaRlyDCJ6wtpWfsU+9MQGgqtzzzmI278nj2Sv72dgo9UB8TCRjh7bn8+XbWZFjbSulLKkYEwCv/riB6Uu2cddZXRiQmuh1OKaGjDkhlYZR4Uz4Zq3XodQallSM+YUWbd7Lw9NXcFrXpvz6xPZeh2NqUHxsJFcObseHi7baXfauoCYVERkuIqtEJEtE7ilnuYjIU+7yxSLSz2fZRBHZISJLy2zzDxFZ6a7/nogkuPNTRSRPRBa604Rgls0YgH15R7jpjfk0jWvA46N6ExZmNzjWN9cPTSMiLIx/f2tnKxDEpCIi4cCzwNlAd+ByEeleZrWzgU7uNA543mfZK8Dwcnb9OdBTVXsBq4E/+Cxbq6p93Gl8QApiTAVUlbvfXsy2ffk8fUVfEmKjvA7JeKBZ4wZc3L81b2Vms+NAvtfheC6YZyoDgSxVXaeqhcAUYGSZdUYCk9QxG0gQkRYAqvotsLvsTlX1M1UtHSlnNtA6aCUwphL/nb2RT5dt4+7hXenXtonX4RgPjT+pPUXFJbz8vY23Esyk0grY7PM8251X3XUqcx3wic/zNBFZICLfiMiJ5W0gIuNEJENEMnJzc6txKGN+snTLPh76aAWndm3K9UPTvA7HeKxdUkPO6dWS12dvYn/+Ea/D8VQwk0p5lct6DOuUv3ORPwFFwOvurBygrar2Be4A3hCRxkftXPUFVU1X1fSUlBR/DmXMzxwsKOKWyQtIbBjFY5daO4px/GZYew4WFDF13uaqV67DgplUsoE2Ps9bA1uPYZ2jiMivgHOBK1VVAVS1QFV3uY8zgbWA3dJsAu7P7y9l465DPDm6D4kNrR3FOHq2imdQWiL/+WEDRcUlXofjmWAmlXlAJxFJE5EoYDQwrcw604Br3KvABgP7VDWnsp2KyHDgbuB8VT3sMz/FvTgAEWmP0/i/LnDFMQbenZ/Nuwu2cOtpnRjUPsnrcEwtc/3QNLbszWPGsu1eh+KZoCUVtzH9ZmAGsAKYqqrLRGS8iJRemTUd54s/C3gRuLF0exGZDMwCuohItohc7y56BogDPi9z6fAwYLGILALeBsar6lEN/cYcq/U7D3Hf+0sZmJrIzad09DocUwud1q0ZqUmxvPR9/f09K27tUb2Unp6uGRkZXodhQkBhUQkXP/8jm3Yf5pPfnkjLhBivQzK11Ks/buD+act454YT6N+ubl4VKCKZqppe3jK7o94YPzz22SqWbNnH3y/uZQnFVOqS/q1p3CCCifX08mJLKsZU4bs1ubzw7TquGNSW4T1tfBRTuYbREVw+qC2fLM1hy948r8OpcZZUjKnEroMF/G7qIjo2bcR955TtEMKY8l09uB0Ab8zZ6HEkNc+SijEVUFXufmcxew8f4anRfYmJCvc6JBMiWjeJ5dSuzZgydzMFRcVeh1OjLKkYU4HX5mziixU7uOfsrnRvedR9tMZU6prj27HrUCGfLNnmdSg1ypKKMeXI2nGAhz9ezrDOKYw5IdXrcEwIGtoxmbTkhvx3dv2qArOkYkwZhUUl/HbKQmKjInjskl7WDYs5JmFhwlWD25G5cU+9GsfekooxZTz++SqWbd3P3y/uRdPGDbwOx4SwS/q1pkFkGP+dVX/OViypGONj1tpdvPDtOi4f2JYzujfzOhwT4uJjI7mgTyveX7iFfYfrR+/FllSMce07fIQ7pi4kLakh953bzetwTB1x1eB25B8p4YNFW7wOpUb4lVRE5FwRsQRk6rT7PlhK7oEC/nVZH2KjIrwOx9QRPVvF07NVYybP3Ux96BbL30QxGlgjIo+KiP2EM3XOBwu3MG3RVn57Wid6t0nwOhxTx1w2oC0rcvazOLvuN9j7lVRU9SqgL84YJf8RkVnuCIpxQY3OmBqwZW8e976/lP7tmnDDyR28DsfUQSP7tCQmMpwp8zZ5HUrQ+V2lpar7gXdwxppvAVwIzBeRW4IUmzFBV1Ki3Dl1ESUlyr9G9SEi3Gp5TeA1bhDJOb1aMG3hVg4VFHkdTlD526Zyvoi8B3wFRAIDVfVsoDdwZxDjMyaoJv6wnlnrdnH/eT1omxTrdTimDrt8YBsOFRbz0eIqB7cNaf7+LLsE+Jeq9lLVf6jqDgB35MXrghadMUG0atsBHp2xijO6N+PS9NZeh2PquH5tm9CxaSMmz63bY9j7m1RyVPVb3xki8ncAVf0y4FEZE2QFRcXc9uZCGjeI4G8XHYeI3TVvgktEGD2gDQs372XVtgNehxM0/iaVM8qZd3YgAzGmJj3xxRpW5OznkYt6kdwo2utwTD1xYd9WRIQJ78zP9jqUoKk0qYjIDSKyBOgqIot9pvXA4poJ0ZjAytiwm39/s5bRA9pwut01b2pQUqNoTunalPcWbKGouMTrcIKiqjOVN4DzgA/cv6VTf/cy40qJyHARWSUiWSJyTznLRUSecpcvFpF+PssmisgOEVlaZptEEflcRNa4f5v4LPuDu69VInJWVfGZ+udgQRF3TF1EqyYx3HuuDbplat7F/VqTe6CA79bs9DqUoKgqqaiqbgBuAg74TIhIYmUbikg48CxONVl34HIRKftffDbQyZ3GAc/7LHsFGF7Oru8BvlTVTsCX7nPcfY8GerjbPefGYMz/PPzxcjbvOcw/R/WhUbTdNW9q3qldm9IkNpK362gVmD9nKgCZQIb7N9PneWUGAlmquk5VC3HubxlZZp2RwCR1zAYSRKQFgHthwO5y9jsSeNV9/Cpwgc/8KapaoKrrgSw3BmMA+GrldibP3cxvhnVgQGqlv4mMCZqoiDBG9mnF58u318lOJitNKqp6rvs3TVXbu39Lp/ZV7LsV4HvtXLY7r7rrlNVMVXPcuHKAptXZl9sTQIaIZOTm5lZxKFNX7D5UyF1vL6Fr8zhuP6OT1+GYeu7ifq0pLCrhoyV1756VSs//fds4yqOq8yvbvLxNjmEdf/m1L1V9AXgBID09ve737mZQVe59fwn78gr57/UDiY6wWlHjrZ6tGtO5WSPezszmykHtvA4noKqqVH68kmUKnFrJ8mygjc/z1kDZtOzPOmVtF5EWqprjVpXt+AX7MvXAtEVbmb5kG3cP70q3FjbWvPGeiHBxv9b87ZOVrMs9SPuURl6HFDBVVX+dUslUWUIBmAd0EpE0EYnCaUSfVmadacA17lVgg4F9pVVblZgG/Mp9/CucK9NK548WkWgRScNp/J9bxb5MHZezL4/73M4ixw2rqsbWmJpzQd9WiMAHC+vWb9+qqr9OVdWvROSi8par6rsVbauqRSJyMzADCAcmquoyERnvLp8ATAdG4DSqHwau9Tn2ZOBkIFlEsoH7VfVl4BFgqohcD2wCLnX3t0xEpgLLgSLgJlUt9uM1MHWUqnLX24spKlH+Oao34TbWvKlFmjVuwPHtk5i2aCu3nd6pzvTqUFX110k4nUieV84yBSpMKgCqOh0ncfjOm+DzWHEuVy5v28srmL8LOK2CZQ8DD1cWk6k/Xpu9ke/W7OShC3rSLqmh1+EYc5SRfVpy9ztLWLJlH71aJ3gdTkBUmlRU9X7377WVrWdMbbN+5yH+b/pKhnVO4cpBbb0Ox5hyDe/ZgvveX8b7C7bWmaTib9f3Se6d7/NFJFNEnhSRpGAHZ8yxKC5Rfjd1IZHhwqMX96oz1Qqm7omPieSUril8uHgrxSV142JUfzuUnALkAhfjdIOfC7wZrKCM+SX+/e1a5m/ay18v6Enz+AZeh2NMpUb2aUXugQJmr9vldSgB4W9SSVTVv6rqend6CEgIYlzGHJPlW/fzr89XM+K45pzfu6XX4RhTpVO7NqVRdAQfLNzidSgB4W9S+VpERotImDuNAj4OZmDGVFdBUTF3TF1IfEwUD11gY6SY0NAgMpyzejTnk6XbyD8S+hesVtX1/QER2Q/8BqcfsEJ3mgLcHvzwjPHfk1+sYeW2A/z94uNIbBjldTjG+G1kn5YcyC9i5qrQ7zqqqpsf41S1sfs3TFUj3ClMVe3WZFNrZG7cw4Rv1jIqvTWndbMxUkxoOaFDEk1iI5m+pKp7v2s/v/v+dsct6QT8r+Wz7BDDxnjhcGERd761iBbxMdxnY6SYEBQRHsbwni34YOEW8o8U0yAydPun8/eS4rHAtzh3xz/o/n0geGEZ479HPlnJ+p2H+MelvYhrEOl1OMYck3OOa8HhwuKQrwLzt6H+t8AAYKOqngL0xbms2BhPfb9mJ5NmbeTaIamc0CHZ63CMOWaD2yeS2DCKj0O8CszfpJKvqvkAIhKtqiuBLsELy5iq7cs7wu/fXkSHlIbcPbyr1+EY84s4VWDN+XLF9pC+CszfpJItIgnA+8DnIvIB1q288diD05ax40AB/xzVJ6TroI0p9VMV2I6qV66l/GqoV9UL3YcPiMjXQDzwadCiMqYKny7N4d0FW7j11I70bpPgdTjGBMSgtESSGkbx0eIchvds4XU4x6Q6V3/1A4bi9E78gzvuvDE1LvdAAX98byk9WzXm5lNtaGBTd5RWgb07fwt5hcXERIXeGbi/V3/9GXgVSAKSgf+IyL3BDMyY8qgqf3h3CQcLivjnqD5ERfhbg2tMaDjnuBbkHQndKjB//yMvBwao6v1ud/iDgSuDF5Yx5Xs7M5svVmznrrO60LlZnNfhGBNwA9MSaRIbyYxl27wO5Zj4m1Q24HPTIxANrA14NMZUYvPuwzz44XIGpiVy3ZA0r8MxJigiwsM4vVszvly5g8KiEq/Dqbaq+v56WkSeAgqAZSLyioj8B1gKHKyJAI0BKClR7nxrEarK45f2JsyGBjZ12PCezTmQX8SsEOwOv6ozlQwgE3gP+CPwNTAT+BPwSVU7F5HhIrJKRLJE5J5ylos7+FeWiCx2LwaodFsReVNEFrrTBhFZ6M5PFZE8n2UTyh7PhK6JP6xnzvrd3H9eD9okxnodjjFBNaRjMg2jwvl0aehVgVU1nPCrpY9FJAro7D5dpapHKttWRMKBZ4EzgGxgnohMU9XlPqudjdOfWCdgEPA8MKiybVX1Mp9jPA7s89nfWlXtU1lcJvSs3n6AR2es4vRuTbk0vbXX4RgTdA0iwzm5a1M+X76dhy7oSXgInZn7e/XXycAanC/654DVIjKsis0GAlmqus69/HgKMLLMOiOBSeqYDSSISAt/thVnsIxRwGR/ymBCU2FRCXdMXUhcdAR/u8iGBjb1x1k9mrPzYAELNu3xOpRq8beh/nHgTFU9SVWHAWcB/6pim1bAZp/n2e48f9bxZ9sTge2qusZnXpqILBCRb0TkxCriMyHgqS/XsHTLfh6+8DhS4qK9DseYGnNKlxSiwsNCrgrM36QSqaqrSp+o6mqgqu5gy/tJqX6u48+2l/Pzs5QcoK2q9gXuAN4QkaPGfBGRcSKSISIZubnWJ2ZtlrlxN8/NzOKS/q0Z3rO51+EYU6PiGkQypGMSM5ZvQ7Xs11/t5W9SyRSRl0XkZHd6EacBvzLZQBuf5605ur+witapdFsRiQAuAt4snaeqBaq6y32ciXPJc2fKUNUXVDVdVdNTUlKqKILxyqGCIu6Y6oyRcv95NkaKqZ+G92zO5t15LM/Z73UofvM3qYwHlgG34nSDv9ydV5l5QCcRSXMb+UcD08qsMw24xr0KbDCwT1Vz/Nj2dGClqmaXzhCRFLeBHxFpj9P4v87P8pla5qGPV7Bp92H+Oaq3jZFi6q3TuzVDBD5fvt3rUPxWZd9fIhIGZKpqT+Cf/u5YVYtE5GacAb3CgYmqukxExrvLJwDTgRFAFnAYuLaybX12P5qjG+iHAX8RkSKgGBivqrv9jdfUHl8s387kuZsYN6w9g9oneR2OMZ5JahRNv7ZN+HLFDm47/aiKl1qpyqSiqiUiskhE2qrqpursXFWn4yQO33kTfB4rcJO/2/osG1POvHeAd6oTn6l9dh4s4J53F9O1eRy/OzM0/omMCabTuzXj75+uJGdfHi3iY7wOp0r+Vn+1wLmj/ksRmVY6BTMwU/+oKve8s4T9eUU8MboP0RGh10OrMYF2eremAHy5IjQ6mPS36/sHgxqFMcCb8zbzxYrt3HtON7o2P+rCPWPqpY5NG9EuKZYvV2znqsHtvA6nSpUmFRFpgNMg3xFYArysqkU1EZipXzbsPMRfPlrOCR2SrLNIY3yICKd1bcZrczZyuLCI2Ci/h8HyRFXVX68C6TgJ5WycmyCNCagjxSXc9uZCIsKEx6yzSGOOcnr3phQWlfDdmp1eh1KlqlJed1U9DkBEXgbmBj8kU988/VUWCzfv5Zkr+tIyofY3RBpT0wakJhLXIIIvlm/nrB61+0bgqs5U/tdppFV7mWDI3LibZ75aw0X9WnFur5Zeh2NMrRQZHsYpXZry1codFJfU7rvrq0oqvUVkvzsdAHqVPhaR0LnF09RKB/KPcNubC2mZEMOD5/fwOhxjarXTuzdj16FCFm6u3R1MVtX1vV3TaYLm/mnL2LInj6m/Od7umjemCid1TiE8TPh6ZS792yV6HU6F/L1PxZiAmrZoK+/O38LNp3YiPbX2/oMYU1vEx0TSv10Tvl5Vu+9XsaRialz2nsP86b0l9GubwK2ndvQ6HGNCxildmrJs636278/3OpQKWVIxNaq4RLn9zYWowpOj+xIRbh9BY/x1SlenZ/WZtfhsxf6jTY165qss5m3Yw18vsLHmjamuLs3iaBHfgK9X1t6xoCypmBozb8NunvxyNRf1bcWFfW2seWOqS0Q4uUtTvs/aSWFRidfhlMuSiqkR+w4f4bYpC2mTGMtfLujpdTjGhKxTuqRwsKCIjI21c2QPSyom6FSVP7y3mO3783lqdF8aRdfuvouMqc2GdEwmMlyYuap2VoFZUjFBN3nuZqYv2cadZ3Whd5sEr8MxJqQ1jI5gUFoSX6+snY31llRMUK3adoAHP1zGiZ2SGXdie6/DMaZOOLlLCmt2HGTz7sNeh3IUSyomaPIKi7n5jfnENYjkn6P6WO/DxgTIKV2dgbtmrq59VWCWVEzQPPjhMrJyD/LEZX1IiYv2Ohxj6oz2yQ1p3SSGb+tbUhGR4SKySkSyROSecpaLiDzlLl8sIv2q2lZEHhCRLSKy0J1G+Cz7g7v+KhE5K5hlM5X7YOEWpszbzA0ndWBop2SvwzGmThERhnVOYdbaXRwprl2XFgctqYhIOPAszuBe3YHLRaR7mdXOBjq50zjgeT+3/Zeq9nGn6e423YHRQA9gOPCcux9Tw9bmHuSP7y4hvV0T7jijs9fhGFMnDeuUzMGCIhZs2ut1KD8TzDOVgUCWqq5T1UJgCjCyzDojgUnqmA0kiEgLP7ctayQwRVULVHU9kOXux9Sg/CPF3PT6fKIiwnj6CuuGxZhgOaFjMuFhUuuqwIL5H98K2OzzPNud5886VW17s1tdNlFEmlTjeIjIOBHJEJGM3Nza9WbUBQ9+uJyV2w7wz8v60CLeRnE0JlgaN4ikb5sEvl1Tu77HgplUyrvUp+yQZRWtU9m2zwMdgD5ADvB4NY6Hqr6gqumqmp6SklLOJuZYfbBwC5PnbuKGkztwSpemXodjTJ13YqcUlmzZx+5DhV6H8j/BTCrZQBuf562BrX6uU+G2qrpdVYtVtQR4kZ+quPw5ngmSNdsPcM87SxiYmmjtKMbUkGGdk1GF77N2eh3K/wQzqcwDOolImohE4TSiTyuzzjTgGvcqsMHAPlXNqWxbt82l1IXAUp99jRaRaBFJw2n8nxuswpmfHCooYvxrmTSMDufpK/oSae0oxtSIXq0TiI+JrFXtKkHrhElVi0TkZmAGEA5MVNVlIjLeXT4BmA6MwGlUPwxcW9m27q4fFZE+OFVbG4DfuNssE5GpwHKgCLhJVYuDVT7jUFX+8O4S1u88xGtjB9GscQOvQzKm3ggPE4Z2TOa7NbmoKiLe32Ac1J793Mt9p5eZN8HnsQI3+butO//qSo73MPDwscZrqm/SrI1MW7SVO8/szAkd7H4UY2rasM7JfLwkh9XbD9KleZzX4dgd9ebYZW7czV8/Ws5pXZty48k2LLAxXjixk3PBUW2pArOkYo7JjgP53PDafFo1ieGfl1m/XsZ4pWVCDO1TGvLD2trRWG9JxVTbkeISbn5jAfvzjzDhqv7Ex0R6HZIx9drQjsnMWbe7VowGaUnFVNvDH69g7vrd/O2i4+jWorHX4RhT7w3pmEzekWIWbt7rdSiWVEz1vJ2ZzSs/buC6IWk2zrwxtcTg9kmESe24X8WSivHbos17+eN7SzihQxJ/HNHV63CMMa74mEiOa53AD5ZUTKjIPVDA+NcySWkUzTNX9LOOIo2pZYZ2TGLh5r0cyD/iaRz2zWCqVFBUzPjXMtlzuJB/X92fxIZRXodkjCljSIdkikuUuet3exqHJRVTKVXlT+8tJXPjHh6/tA89W8V7HZIxphz92jUhOiLM83YVSyqmUi99t563M7P57WmdOKdXi6o3MMZ4okFkOAPTEvkxa5encVhSMRX6auV2/u+TFYw4rjm/Pa2T1+EYY6pwQodkVm0/wI4D+Z7FYEnFlGv51v3c8sYCerRszGOX9rY75o0JAUM7Ov3veXm2YknFHGX7/nyuf3UecQ0ieflXA4iNCmq/o8aYAOnesjHxMZH86GGXLfZtYX7mcGERY1/NYF/eEd4af7x1ZW9MCAkPEwalJTJrnZ2pmFqgqLiEW95YwLKt+3j68r70aGlXehkTao7vkMTm3Xlk7znsyfEtqRjAuXT4z9OW8eXKHTxwfg9O69bM65CMMcdgcPskAGav8+Z+FUsqBoDnZq7ljTmbGH9SB645PtXrcIwxx6hLsziaxEYya603VWCWVAxvZ2bzjxmruKBPS+46q4vX4RhjfoGwMGFw+yRmr9uFM7huDR+/xo9oapUvlm/n7ncWM6RjEo9eYpcOG1MXDG6fxJa9eWTvyavxYwc1qYjIcBFZJSJZInJPOctFRJ5yly8WkX5VbSsi/xCRle7674lIgjs/VUTyRGShO00IZtnqgrnrd3PTG/Pp0bIx/746nagI+41hTF1wfAenXcWLKrCgfYuISDjwLHA20B24XES6l1ntbKCTO40Dnvdj28+BnqraC1gN/MFnf2tVtY87jQ9OyeqG5Vv3c/2r82jVJIb/jBlAo2i7utyYuqJT00YkNYzy5NLiYP40HQhkqeo6VS0EpgAjy6wzEpikjtlAgoi0qGxbVf1MVYvc7WcDNlJUNWXtOMjVL8+hUXQE/71+EEmNor0OyRgTQCLetasEM6m0Ajb7PM925/mzjj/bAlwHfOLzPE1EFojINyJyYnlBicg4EckQkYzc3Fz/SlKHbNx1iCtfmo2I8PrYQbRKiPE6JGNMEAzukETOvnw27qrZ+1WCmVTKa/EtmzIrWqfKbUXkT0AR8Lo7Kwdoq6p9gTuAN0TkqAHUVfUFVU1X1fSUlJQqilC3bN2bxxUvzqGgqITXxw6ifUojr0MyxgTJ8e79KjVdBRbMpJINtPF53hrY6uc6lW4rIr8CzgWuVPfcTlULVHWX+zgTWAt0DkhJ6oCte/MY/cJs9ucd4b/XDaJL8zivQzLGBFGHlIakxEUzuw4llXlAJxFJE5EoYDQwrcw604Br3KvABgP7VDWnsm1FZDhwN3C+qv7vvE5EUtwGfkSkPU7j/7ogli9kbHETyp5DhUy6fiDHtbbuV4yp60SEgWmJzF2/u0bbVYKWVNzG9JuBGcAKYKqqLhOR8SJSemXWdJwv/izgReDGyrZ1t3kGiAM+L3Pp8DBgsYgsAt4Gxquqt+Nq1gJOQpnFnsOF/HfsIPq2beJ1SMaYGjIoLZGcffk1er9KUK8jVdXpOInDd94En8cK3OTvtu78jhWs/w7wzi+Jt65Zv/MQV700h/35R3jt+kH0bpPgdUjGmBo0MC0RcO5Ja5MYWyPHtLvd6qgVOfu5dMIs8o4UM/nXgy2hGFMPdW4aR3xMJHPX11yljd3xVgct2LSHMf+ZR0xkOK+NHUzHpnaVlzH1UViYMCA1kTnra66x3s5U6pivVm7nihfnEB8TyVvjj7eEYkw9NygtkQ27DrN9f82MW29JpQ6ZMncTv56UScemjXjnhhNqrA7VGFN7+bar1ARLKnVASYnyz89Xc8+7SxjaMZkp4waTEmddrxhjoEfLxsRGhddYUrE2lRCXV1jMnW8t4uMlOYxKb83DFx5HZLj9VjDGOCLCw+jfromdqZiqbduXz6h/z2L60hz+OKIrf7+4lyUUY8xRBrdPYtX2A+w5VBj0Y9k3UIias24X5z79PetyD/LSNemMG9YBERtgyxhztNJ2lXkbgn+2YkklxKgqL3y7litemkPjBhG8d9MQTuvWzOuwjDG1WK/W8URFhDGnBqrArE0lhOw5VMg97y5mxrLtnN2zOY9e0ou4BpFeh2WMqeWiI8Lp0zqBjBo4U7GkEiJ+zNrJHVMXsetQAfee043rh6ZZdZcxxm/pqU144dt1HC4sIjYqeF/9Vv1Vy+UfKeb/pq/gypfnEBsdzns3DmHsie0toRhjqmVAaiJFJcrCzXuDehw7U6nF5qzbxT3vLmH9zkNcMagt957TLai/MIwxdVc/t4fyzA17OKFDctCOY99QtdCeQ4U89tkqXp+ziTaJMbw+dhBDOgbvQ2CMqfviYyPp0iyOeRv3BPU4llRqkaLiEt6Yu4nHP1vNwYIirhuSxp1ndbazE2NMQKSnNmHawq0UlyjhYcGpQrdvq1pAVZmxbBv//Hw1q7cf5IQOSdx/Xg8b8tcYE1DpqU14fc4mVm07QPeWjYNyDEsqHiopUb5etYN/fbGapVv20z6lIc9f2Y/hPZtbQ7wxJuDS2zk3QWZs3G1JpS7JKyzmnfnZTPxhPetyD9EmMYbHLu3NBX1aEmHdrBhjgqR1kxiaNY4mY8Merjk+NSjHsKRSQ1SV+Zv28u78bD5ctJX9+UUc1yqeJy7rwzm9WlifXcaYoBMR0lMTg3oTZFCTiogMB54EwoGXVPWRMsvFXT4COAyMUdX5lW0rIonAm0AqsAEYpap73GV/AK4HioFbVXVGMMtXlYKiYuas283Xq3bw5YodbNp9mAaRYQzv0ZwrBrVjQGoTq+YyxtSoAe2a8PHiHLbszaNVQkzA9x+0pCIi4cCzwBlANjBPRKap6nKf1c4GOrnTIOB5YFAV294DfKmqj4jIPe7zu0WkOzAa6AG0BL4Qkc6qWhysMpZSVXIPFrB5dx6bdh9iSfZ+FmXvZemWfRQUlRAdEcbxHZK45dSOnH1cCxpF2wmiMcYb6aluu8qG3bTq0yrg+w/mt9tAIEtV1wGIyBRgJOCbVEYCk1RVgdkikiAiLXDOQiradiRwsrv9q8BM4G53/hRVLQDWi0iWG8OsQBdsRc5+bpm8gLzCYg4XFnGooJjC4pL/LW8QGUbPlvFcNbgdQzomcXz7ZGKiwgMdhjHGVFvX5nE0jAonY8MeRoZYUmkFbPZ5no1zNlLVOq2q2LaZquYAqGqOiDT12dfscvb1MyIyDhgH0LZt22oU5ycNoyLo3KwRMZERxEaFExsdTsv4GNokxtCmSSxpyQ2twd0YUytFhIcxakAbWjcJznDjwUwq5TUWqJ/r+LPtsRwPVX0BeAEgPT29qn2Wq21SLM9d2f9YNjXGGM/df16PoO07mD+ns4E2Ps9bA1v9XKeybbe7VWS4f3dU43jGGGOCKJhJZR7QSUTSRCQKpxF9Wpl1pgHXiGMwsM+t2qps22nAr9zHvwI+8Jk/WkSiRSQNp/F/brAKZ4wx5mhBq/5S1SIRuRmYgXNZ8ERVXSYi493lE4DpOJcTZ+FcUnxtZdu6u34EmCoi1wObgEvdbZaJyFScxvwi4KaauPLLGGPMT8S58Kp+Sk9P14yMDK/DMMaYkCIimaqaXt4yu0TJGGNMwFhSMcYYEzCWVIwxxgSMJRVjjDEBU68b6kUkF9j4C3aRDOwMUDihoL6VF6zM9YWVuXraqWpKeQvqdVL5pUQko6IrIOqi+lZesDLXF1bmwLHqL2OMMQFjScUYY0zAWFL5ZV7wOoAaVt/KC1bm+sLKHCDWpmKMMSZg7EzFGGNMwFhSMcYYEzCWVKogIsNFZJWIZInIPeUsFxF5yl2+WET6eRFnIPlR5ivdsi4WkR9FpLcXcQZSVWX2WW+AiBSLyCU1GV8w+FNmETlZRBaKyDIR+aamYww0Pz7b8SLyoYgscst8rRdxBoqITBSRHSKytILlgf/+UlWbKphwut1fC7QHooBFQPcy64wAPsEZeXIwMMfruGugzCcATdzHZ9eHMvus9xXOkA2XeB13DbzPCThDSbR1nzf1Ou4aKPMfgb+7j1OA3UCU17H/gjIPA/oBSytYHvDvLztTqdxAIEtV16lqITAFGFlmnZHAJHXMBhJKR6YMUVWWWVV/VNU97tPZOKNshjJ/3meAW4B3+Gm00VDmT5mvAN5V1U0Aqhrq5fanzArEiYgAjXCSSlHNhhk4qvotThkqEvDvL0sqlWsFbPZ5nu3Oq+46oaS65bke55dOKKuyzCLSCrgQmFCDcQWTP+9zZ6CJiMwUkUwRuabGogsOf8r8DNANZyjyJcBvVbWkZsLzRMC/v4I28mMdIeXMK3sNtj/rhBK/yyMip+AklaFBjSj4/CnzE8Ddqlrs/IgNef6UOQLoD5wGxACzRGS2qq4OdnBB4k+ZzwIWAqcCHYDPReQ7Vd0f5Ni8EvDvL0sqlcsG2vg8b43zC6a664QSv8ojIr2Al4CzVXVXDcUWLP6UOR2Y4iaUZGCEiBSp6vs1EmHg+fvZ3qmqh4BDIvIt0BsI1aTiT5mvBR5Rp8EhS0TWA12BuTUTYo0L+PeXVX9Vbh7QSUTSRCQKGA1MK7PONOAa9yqKwcA+Vc2p6UADqMoyi0hb4F3g6hD+1eqryjKrapqqpqpqKvA2cGMIJxTw77P9AXCiiESISCwwCFhRw3EGkj9l3oRzZoaINAO6AOtqNMqaFfDvLztTqYSqFonIzcAMnCtHJqrqMhEZ7y6fgHMl0AggCziM80snZPlZ5j8DScBz7i/3Ig3hHl79LHOd4k+ZVXWFiHwKLAZKgJdUtdxLU0OBn+/zX4FXRGQJTtXQ3aoasl3ii8hk4GQgWUSygfuBSAje95d102KMMSZgrPrLGGNMwFhSMcYYEzCWVIwxxgSMJRVjjDEBY0nFGGNMwFhSMbWW2xvwQp8p9Rfur4+IjPB5fn5lPRIHgojcKiIrROT1IB5juogkuI8Pun9TK+qZNphE5JXyenAWkZdEpHtNx2Nqnt2nYmqzPFXtU94Ct8M/qWa/TH1w7oyfDqCq0zj65rdAuxGn14H1wTqAqo6oei1vqepYr2MwNcPOVEzIcH99rxCR54D5QBsReV5EMtyxLx70WXeAOGO9LBKRuSISD/wFuMw967lMRMaIyDPu+u1E5Et3TIkv3V4DSn95P+Xua115v8Ld9e4QkaXudJs7bwJON+vTROT2MuvPEZEePs9nikh/ERnoHmuB+7eLu3yMiLwrIp+KyBoRedRn2w0iklzF6/adiMx3pxMqWO8at/yLROS/Vbwu5c4vs7+/uq9fmFu+dHf+mSIyy43lLRFp5M5/RESWu/t8rKLymFrO6/7+bbKpogkoxuncbyHwHpCKc2f3YJ91Et2/4cBMoBfOWBnrgAHussY4Z+VjgGd8tv3fc+BD4Ffu4+uA993HrwBv4fwA647TdXrZOPvj9GjbEKe79GVAX3fZBiC5nG1uBx50H7cAVvvG6j4+HXjHJ9Z1QDzQANgItCl7DOCg+zcVdwwNIBZo4D7uBGSUE08PYJXPfhKreF0qe70uAR4F/s1PN1jPxDlLTAa+BRq68+/G6aEh0T1+6foJXn/+bDq2yaq/TG32s+ovt01lozrjPpQaJSLjcJJGC5wvfgVyVHUegLo9zErlvQsfD1zkPv4vzpdiqffVqWZbLk5/UGUNBd5Tp+NFRORd4ERgQSXHmwp8jtNtxiicxAVO0nhVRDq55Yj02eZLVd3nHmM50I6fd1tekUjgGRHpg5OoO5ezzqnA2+p2SaKqpWNwVPS6VPZ63Ycz2NO4co4zGOc9+sF9P6KAWcB+IB94SUQ+Bj7yo1ymFrKkYkLNodIHIpIG3IlzRrJHRF7B+RUv/PLhB3y3L/B5XF5mqnZf+Kq6RUR2idPb82XAb9xFfwW+VtUL3SQ6s4I4ivH///d2YDtOD8NhOF/eZfn7mlW0ju/8eUB/EUn0SU6+x/lcVS8/KgCRgTidOY4GbsZJdCbEWJuKCWWNcZLMPvcM4mx3/kqgpYgMABCROBGJAA4AcRXs60ecLzOAK4HvqxHHt8AFIhIrIg1xBvP6zo/tpgB3AfGqusSdFw9scR+PqUYMlYnHOXMrAa7GqSos60ucs74kABFJdOdX9LpU9np9CjwCfCwiZV/v2cAQEenoHidWRDq77SrxqjoduA3nogoTguxMxYQsVV0kIgtw2jDWAT+48wtF5DLgaRGJAfJw2ie+Bu4RkYXA38rs7lZgooj8HsilGr21qup89yypdMyNl1S1sqqvUm8DT+KcnZR6FKf66w7gK39jqMJzwDsicinOa3Co7Arq9Nb7MPCNiBTjVN2NoeLXpdLXS1XfchPKNPG5jFtVc0VkDDBZRKLd2ffiJPwPRKT0TPNnFzaY0GG9FBtjjAkYq/4yxhgTMJZUjDHGBIwlFWOMMQFjScUYY0zAWFIxxhgTMJZUjDHGBIwlFWOMMQHz/xL+xmh+eeKKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "posterior3.plot()\n", "\n", "plt.xlabel('Fraction of vanilla cookies')\n", "plt.ylabel('Probability')\n", "plt.title('Posterior, two vanilla, one chocolate');" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.67" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "posterior3.idxmax()" ] } ], "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.8.5" } }, "nbformat": 4, "nbformat_minor": 1 }