{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Distributions"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"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": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:12.094262Z",
"iopub.status.busy": "2021-04-16T19:35:12.093730Z",
"iopub.status.idle": "2021-04-16T19:35:12.095822Z",
"shell.execute_reply": "2021-04-16T19:35:12.095402Z"
},
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# If we're running on Colab, install empiricaldist\n",
"# https://pypi.org/project/empiricaldist/\n",
"\n",
"import sys\n",
"IN_COLAB = 'google.colab' in sys.modules\n",
"\n",
"if IN_COLAB:\n",
" !pip install empiricaldist"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:12.099615Z",
"iopub.status.busy": "2021-04-16T19:35:12.098997Z",
"iopub.status.idle": "2021-04-16T19:35:12.100895Z",
"shell.execute_reply": "2021-04-16T19:35:12.101335Z"
},
"tags": [
"remove-cell"
]
},
"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,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:12.104987Z",
"iopub.status.busy": "2021-04-16T19:35:12.104287Z",
"iopub.status.idle": "2021-04-16T19:35:13.067675Z",
"shell.execute_reply": "2021-04-16T19:35:13.068108Z"
},
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"from utils import set_pyplot_params\n",
"set_pyplot_params()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the previous chapter we used Bayes's Theorem to solve a cookie problem; then we solved it again using a Bayes table.\n",
"In this chapter, at the risk of testing your patience, we will solve it one more time using a `Pmf` object, which represents a \"probability mass function\".\n",
"I'll explain what that means, and why it is useful for Bayesian statistics.\n",
"\n",
"We'll use `Pmf` objects to solve some more challenging problems and take one more step toward Bayesian statistics.\n",
"But we'll start with distributions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Distributions\n",
"\n",
"In statistics a **distribution** is a set of possible outcomes and their corresponding probabilities.\n",
"For example, if you toss a coin, there are two possible outcomes with\n",
"approximately equal probability.\n",
"If you roll a six-sided die, the set of possible outcomes is the numbers 1 to 6, and the probability associated with each outcome is 1/6.\n",
"\n",
"To represent distributions, we'll use a library called `empiricaldist`.\n",
"An \"empirical\" distribution is based on data, as opposed to a\n",
"theoretical distribution.\n",
"We'll use this library throughout the book. I'll introduce the basic features in this chapter and we'll see additional features later."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Probability Mass Functions\n",
"\n",
"If the outcomes in a distribution are discrete, we can describe the distribution with a **probability mass function**, or PMF, which is a function that maps from each possible outcome to its probability.\n",
"\n",
"`empiricaldist` provides a class called `Pmf` that represents a\n",
"probability mass function.\n",
"To use `Pmf` you can import it like this:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.071512Z",
"iopub.status.busy": "2021-04-16T19:35:13.071040Z",
"iopub.status.idle": "2021-04-16T19:35:13.073278Z",
"shell.execute_reply": "2021-04-16T19:35:13.072869Z"
}
},
"outputs": [],
"source": [
"from empiricaldist import Pmf"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"If that doesn't work, you might have to install `empiricaldist`; try running\n",
"\n",
"```\n",
"!pip install empiricaldist\n",
"```\n",
"\n",
"in a code cell or\n",
"\n",
"```\n",
"pip install empiricaldist\n",
"```\n",
"\n",
"in a terminal window."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following example makes a `Pmf` that represents the outcome of a\n",
"coin toss."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.085000Z",
"iopub.status.busy": "2021-04-16T19:35:13.081870Z",
"iopub.status.idle": "2021-04-16T19:35:13.088193Z",
"shell.execute_reply": "2021-04-16T19:35:13.087667Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
heads
\n",
"
0.5
\n",
"
\n",
"
\n",
"
tails
\n",
"
0.5
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"heads 0.5\n",
"tails 0.5\n",
"Name: , dtype: float64"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coin = Pmf()\n",
"coin['heads'] = 1/2\n",
"coin['tails'] = 1/2\n",
"coin"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`Pmf` creates an empty `Pmf` with no outcomes.\n",
"Then we can add new outcomes using the bracket operator.\n",
"In this example, the two outcomes are represented with strings, and they have the same probability, 0.5."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also make a `Pmf` from a sequence of possible outcomes.\n",
"\n",
"The following example uses `Pmf.from_seq` to make a `Pmf` that represents a six-sided die."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.092152Z",
"iopub.status.busy": "2021-04-16T19:35:13.091700Z",
"iopub.status.idle": "2021-04-16T19:35:13.101915Z",
"shell.execute_reply": "2021-04-16T19:35:13.102297Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
0.166667
\n",
"
\n",
"
\n",
"
2
\n",
"
0.166667
\n",
"
\n",
"
\n",
"
3
\n",
"
0.166667
\n",
"
\n",
"
\n",
"
4
\n",
"
0.166667
\n",
"
\n",
"
\n",
"
5
\n",
"
0.166667
\n",
"
\n",
"
\n",
"
6
\n",
"
0.166667
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"1 0.166667\n",
"2 0.166667\n",
"3 0.166667\n",
"4 0.166667\n",
"5 0.166667\n",
"6 0.166667\n",
"Name: , dtype: float64"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"die = Pmf.from_seq([1,2,3,4,5,6])\n",
"die"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, all outcomes in the sequence appear once, so they all have the same probability, $1/6$.\n",
"\n",
"More generally, outcomes can appear more than once, as in the following example:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.109370Z",
"iopub.status.busy": "2021-04-16T19:35:13.108697Z",
"iopub.status.idle": "2021-04-16T19:35:13.111096Z",
"shell.execute_reply": "2021-04-16T19:35:13.111446Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
M
\n",
"
0.090909
\n",
"
\n",
"
\n",
"
i
\n",
"
0.363636
\n",
"
\n",
"
\n",
"
p
\n",
"
0.181818
\n",
"
\n",
"
\n",
"
s
\n",
"
0.363636
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"M 0.090909\n",
"i 0.363636\n",
"p 0.181818\n",
"s 0.363636\n",
"Name: , dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"letters = Pmf.from_seq(list('Mississippi'))\n",
"letters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The letter `M` appears once out of 11 characters, so its probability is $1/11$.\n",
"The letter `i` appears 4 times, so its probability is $4/11$.\n",
"\n",
"Since the letters in a string are not outcomes of a random process, I'll use the more general term \"quantities\" for the letters in the `Pmf`.\n",
"\n",
"The `Pmf` class inherits from a Pandas `Series`, so anything you can do with a `Series`, you can also do with a `Pmf`.\n",
"\n",
"For example, you can use the bracket operator to look up a quantity and get the corresponding probability."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.116296Z",
"iopub.status.busy": "2021-04-16T19:35:13.115451Z",
"iopub.status.idle": "2021-04-16T19:35:13.118339Z",
"shell.execute_reply": "2021-04-16T19:35:13.118724Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.36363636363636365"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"letters['s']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the word \"Mississippi\", about 36% of the letters are \"s\".\n",
"\n",
"However, if you ask for the probability of a quantity that's not in the distribution, you get a `KeyError`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.122337Z",
"iopub.status.busy": "2021-04-16T19:35:13.121841Z",
"iopub.status.idle": "2021-04-16T19:35:13.124479Z",
"shell.execute_reply": "2021-04-16T19:35:13.124086Z"
},
"tags": [
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"try:\n",
" letters['t']\n",
"except KeyError as e:\n",
" print(type(e))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also call a `Pmf` as if it were a function, with a letter in parentheses."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.128416Z",
"iopub.status.busy": "2021-04-16T19:35:13.127840Z",
"iopub.status.idle": "2021-04-16T19:35:13.130899Z",
"shell.execute_reply": "2021-04-16T19:35:13.130488Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.36363636363636365"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"letters('s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the quantity is in the distribution the results are the same. \n",
"But if it is not in the distribution, the result is `0`, not an error."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.134469Z",
"iopub.status.busy": "2021-04-16T19:35:13.133893Z",
"iopub.status.idle": "2021-04-16T19:35:13.136395Z",
"shell.execute_reply": "2021-04-16T19:35:13.136752Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"letters('t')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With parentheses, you can also provide a sequence of quantities and get a sequence of probabilities."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.140773Z",
"iopub.status.busy": "2021-04-16T19:35:13.140250Z",
"iopub.status.idle": "2021-04-16T19:35:13.142982Z",
"shell.execute_reply": "2021-04-16T19:35:13.142616Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0.16666667, 0.16666667, 0. ])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"die([1,4,7])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The quantities in a `Pmf` can be strings, numbers, or any other type that can be stored in the index of a Pandas `Series`.\n",
"If you are familiar with Pandas, that will help you work with `Pmf` objects. \n",
"But I will explain what you need to know as we go along."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Cookie Problem Revisited\n",
"\n",
"In this section I'll use a `Pmf` to solve the cookie problem from <<_TheCookieProblem>>.\n",
"Here's the statement of the problem again:\n",
"\n",
"> Suppose there are two bowls of cookies.\n",
">\n",
"> * Bowl 1 contains 30 vanilla cookies and 10 chocolate cookies. \n",
">\n",
"> * Bowl 2 contains 20 vanilla cookies and 20 chocolate cookies.\n",
">\n",
"> Now suppose you choose one of the bowls at random and, without looking, choose a cookie at random. If the cookie is vanilla, what is the probability that it came from Bowl 1?\n",
"\n",
"Here's a `Pmf` that represents the two hypotheses and their prior probabilities:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.149760Z",
"iopub.status.busy": "2021-04-16T19:35:13.149078Z",
"iopub.status.idle": "2021-04-16T19:35:13.152305Z",
"shell.execute_reply": "2021-04-16T19:35:13.151707Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Bowl 1
\n",
"
0.5
\n",
"
\n",
"
\n",
"
Bowl 2
\n",
"
0.5
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"Bowl 1 0.5\n",
"Bowl 2 0.5\n",
"Name: , dtype: float64"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prior = Pmf.from_seq(['Bowl 1', 'Bowl 2'])\n",
"prior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This distribution, which contains the prior probability for each hypothesis, is called (wait for it) the **prior distribution**.\n",
"\n",
"To update the distribution based on new data (the vanilla cookie),\n",
"we multiply the priors by the likelihoods. The likelihood\n",
"of drawing a vanilla cookie from Bowl 1 is `3/4`. The likelihood\n",
"for Bowl 2 is `1/2`."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.158260Z",
"iopub.status.busy": "2021-04-16T19:35:13.157620Z",
"iopub.status.idle": "2021-04-16T19:35:13.160416Z",
"shell.execute_reply": "2021-04-16T19:35:13.160045Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Bowl 1
\n",
"
0.375
\n",
"
\n",
"
\n",
"
Bowl 2
\n",
"
0.250
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"Bowl 1 0.375\n",
"Bowl 2 0.250\n",
"Name: , dtype: float64"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likelihood_vanilla = [0.75, 0.5]\n",
"posterior = prior * likelihood_vanilla\n",
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is the unnormalized posteriors; that is, they don't add up to 1.\n",
"To make them add up to 1, we can use `normalize`, which is a method provided by `Pmf`."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.164762Z",
"iopub.status.busy": "2021-04-16T19:35:13.164137Z",
"iopub.status.idle": "2021-04-16T19:35:13.167103Z",
"shell.execute_reply": "2021-04-16T19:35:13.166662Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.625"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The return value from `normalize` is the total probability of the data, which is $5/8$.\n",
"\n",
"`posterior`, which contains the posterior probability for each hypothesis, is called (wait now) the **posterior distribution**."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.172660Z",
"iopub.status.busy": "2021-04-16T19:35:13.171984Z",
"iopub.status.idle": "2021-04-16T19:35:13.174476Z",
"shell.execute_reply": "2021-04-16T19:35:13.174842Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Bowl 1
\n",
"
0.6
\n",
"
\n",
"
\n",
"
Bowl 2
\n",
"
0.4
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"Bowl 1 0.6\n",
"Bowl 2 0.4\n",
"Name: , dtype: float64"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the posterior distribution we can select the posterior probability for Bowl 1:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.178618Z",
"iopub.status.busy": "2021-04-16T19:35:13.177927Z",
"iopub.status.idle": "2021-04-16T19:35:13.180255Z",
"shell.execute_reply": "2021-04-16T19:35:13.180610Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.6"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior('Bowl 1')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the answer is 0.6.\n",
"\n",
"One benefit of using `Pmf` objects is that it is easy to do successive updates with more data.\n",
"For example, suppose you put the first cookie back (so the contents of the bowls don't change) and draw again from the same bowl.\n",
"If the second cookie is also vanilla, we can do a second update like this:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.187902Z",
"iopub.status.busy": "2021-04-16T19:35:13.187276Z",
"iopub.status.idle": "2021-04-16T19:35:13.189938Z",
"shell.execute_reply": "2021-04-16T19:35:13.190290Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Bowl 1
\n",
"
0.692308
\n",
"
\n",
"
\n",
"
Bowl 2
\n",
"
0.307692
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"Bowl 1 0.692308\n",
"Bowl 2 0.307692\n",
"Name: , dtype: float64"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior *= likelihood_vanilla\n",
"posterior.normalize()\n",
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the posterior probability for Bowl 1 is almost 70%.\n",
"But suppose we do the same thing again and get a chocolate cookie.\n",
"\n",
"Here are the likelihoods for the new data:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.193226Z",
"iopub.status.busy": "2021-04-16T19:35:13.192782Z",
"iopub.status.idle": "2021-04-16T19:35:13.194564Z",
"shell.execute_reply": "2021-04-16T19:35:13.194987Z"
}
},
"outputs": [],
"source": [
"likelihood_chocolate = [0.25, 0.5]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's the update."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.202639Z",
"iopub.status.busy": "2021-04-16T19:35:13.202008Z",
"iopub.status.idle": "2021-04-16T19:35:13.204858Z",
"shell.execute_reply": "2021-04-16T19:35:13.205318Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Bowl 1
\n",
"
0.529412
\n",
"
\n",
"
\n",
"
Bowl 2
\n",
"
0.470588
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"Bowl 1 0.529412\n",
"Bowl 2 0.470588\n",
"Name: , dtype: float64"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior *= likelihood_chocolate\n",
"posterior.normalize()\n",
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the posterior probability for Bowl 1 is about 53%.\n",
"After two vanilla cookies and one chocolate, the posterior probabilities are close to 50/50."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 101 Bowls\n",
"\n",
"Next let's solve a cookie problem with 101 bowls:\n",
"\n",
"* Bowl 0 contains 0% 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 version, 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 solve this problem, I'll use `np.arange` to make an array that represents 101 hypotheses, numbered from 0 to 100."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.209286Z",
"iopub.status.busy": "2021-04-16T19:35:13.208453Z",
"iopub.status.idle": "2021-04-16T19:35:13.210850Z",
"shell.execute_reply": "2021-04-16T19:35:13.211365Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"hypos = np.arange(101)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use this array to make the prior distribution:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.217907Z",
"iopub.status.busy": "2021-04-16T19:35:13.217166Z",
"iopub.status.idle": "2021-04-16T19:35:13.220630Z",
"shell.execute_reply": "2021-04-16T19:35:13.219991Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"101"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prior = Pmf(1, hypos)\n",
"prior.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As this example shows, we can initialize a `Pmf` with two parameters.\n",
"The first parameter is the prior probability; the second parameter is a sequence of quantities.\n",
"\n",
"In this example, the probabilities are all the same, so we only have to provide one of them; it gets \"broadcast\" across the hypotheses.\n",
"Since all hypotheses have the same prior probability, this distribution is **uniform**.\n",
"\n",
"Here are the first few hypotheses and their probabilities."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.228558Z",
"iopub.status.busy": "2021-04-16T19:35:13.227628Z",
"iopub.status.idle": "2021-04-16T19:35:13.232003Z",
"shell.execute_reply": "2021-04-16T19:35:13.231443Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.009901
\n",
"
\n",
"
\n",
"
1
\n",
"
0.009901
\n",
"
\n",
"
\n",
"
2
\n",
"
0.009901
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"0 0.009901\n",
"1 0.009901\n",
"2 0.009901\n",
"Name: , dtype: float64"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prior.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The likelihood of the data is the fraction of vanilla cookies in each bowl, which we can calculate using `hypos`:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.237433Z",
"iopub.status.busy": "2021-04-16T19:35:13.236599Z",
"iopub.status.idle": "2021-04-16T19:35:13.239986Z",
"shell.execute_reply": "2021-04-16T19:35:13.240390Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0. , 0.01, 0.02, 0.03, 0.04])"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likelihood_vanilla = hypos/100\n",
"likelihood_vanilla[:5]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can compute the posterior distribution in the usual way:\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.247286Z",
"iopub.status.busy": "2021-04-16T19:35:13.246696Z",
"iopub.status.idle": "2021-04-16T19:35:13.249668Z",
"shell.execute_reply": "2021-04-16T19:35:13.249258Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"prior.plot(label='prior', color='C5')\n",
"posterior1.plot(label='posterior', color='C4')\n",
"decorate_bowls('Posterior after one vanilla cookie')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior probability of Bowl 0 is 0 because it contains no vanilla cookies.\n",
"The posterior probability of Bowl 100 is the highest because it contains the most vanilla cookies.\n",
"In between, the shape of the posterior distribution is a line because the likelihoods are proportional to the bowl numbers.\n",
"\n",
"Now suppose we put the cookie back, draw again from the same bowl, and get another vanilla cookie.\n",
"Here's the update after the second cookie:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.493905Z",
"iopub.status.busy": "2021-04-16T19:35:13.493399Z",
"iopub.status.idle": "2021-04-16T19:35:13.497012Z",
"shell.execute_reply": "2021-04-16T19:35:13.497425Z"
},
"tags": [
"hide-output"
]
},
"outputs": [
{
"data": {
"text/plain": [
"0.6699999999999999"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior2 = posterior1 * likelihood_vanilla\n",
"posterior2.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's what the posterior distribution looks like."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.520648Z",
"iopub.status.busy": "2021-04-16T19:35:13.513566Z",
"iopub.status.idle": "2021-04-16T19:35:13.640728Z",
"shell.execute_reply": "2021-04-16T19:35:13.640354Z"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA37UlEQVR4nO3dd5wV5dn/8c+1BZbewaUJUqQozaWJilgiIko09gKWxJKYxDzmSUjiEzWahF8S46OJD2oCitiisRGjsSCCBZSliBSBFSlLXXqHLdfvj5ldzi7bKGfP2T3f9+t1Xjtzz9wz18yePdfOzH3u29wdERGReJMU6wBERERKowQlIiJxSQlKRETikhKUiIjEJSUoERGJS0pQIiISl5SgpNoxs+vM7N0q3N+DZrbZzDZU1T6rk5K/DzNzM+scTj9tZg/GLrrizGylmZ1XxrJFZnZ21UYk5VGCknKFf9D7zGy3mW00s6fMrP4xbO8+M3v2WGJy9+fc/VvHso3KMrN2wN1AD3c/wcxuNLOPj/M+jvs2q1JV/j6iyd17uvuHsY5DDlGCksq42N3rA/2A/sA9sQrEzFKOoa6Z2ZG+508Etrj7pqPdb4kYjjp+kUSjBCWV5u5rgbeBUwDM7JLwtsh2M/vQzLoXrmtmPzeztWa2y8yWmtm5ZjYc+CVwVXhF9kW4biMzm2Bm68M6D5pZcrjsRjP7xMweNrOtwH0lrzjM7HQzm21mO8Kfp0cs+9DMfmtmnwB7gZNKHpeZjTWzr8NYF5vZpWH5ecB7QOsw3n8AjwODw/nt4Xq1zexPZrY6vMp83MzqhMvONrPs8HxsAJ4qse/uJbdpZh3Dn0nhOn83s00RdZ41s7vC6dZmNsXMtppZlpl9r7TfnZkNMrMNhec1LLvUzBaE0wPMbGa43/Vm9lczqxWxrpvZ7Wa23My2mdljZmYRv6MKrwDNrImZvWlmOeE23jSztuWs387MXg3X32Jmfw3Lk8zsHjNbZWabzOwZM2sUUa/M92WJ7Xczs2/M7Opwvuj2X7iPwvfFFjN7ycyahsvSwt/BlnAfs82sVUXHL0fB3fXSq8wXsBI4L5xuBywCHgC6AnuA84FU4GdAFlALOBlYA7QO63UAOoXT9wHPltjH68ATQD2gJfA5cFu47EYgD/ghkALUCcs+Dpc3BbYBN4TLrwnnm4XLPwRWAz3D5amlHOMVQGuCf9iuCo8rPVx2NpAdsW7RviPK/heYEsbSAPgX8PuI+nnA/wNqA3VK2X9p21wNnBZOLwVWAN0jlvUNp6cD/wekAX2AHODcMn6XXwPnR8y/DIwNp08DBoXnqAOwBLgrYl0H3gQaA+3D/QwvLf5w3c7h9NPAg+F0M+A7QN3wPL0MvF5GrMnAF8DD4fsiDTgjXHYzwXvtJKA+8CowOVxW5vsy8v1McDdgNTCyjPf6XcAsoG34e3sCeCFcdlv4O64bxnka0DDWf6s18RXzAPSK71f4R7sb2A6sCj8M6wD/A7wUsV4SsJbgA7kzsCn8IEgtsb37iEhQQCvgABEf3ARJZlo4fSOwusQ2ij4QCRLT5yWWzwRuDKc/BH5zhMc8HxgVTp9NOQkKsPADsVNE2WDgm4j6B4G0cvZXbJth2WTgv4ATCBLUH4DbgY7h7yKJ4B+GfKBBRL3fA0+XsZ8HgYnhdIMw7hPLWPcu4LWIeSdMEOH8SxxKbiXPSakJqpR99AG2lbFsMEESTCll2VTg+xHzJwO5BMm1zPdlxPv5fiAbGFbKe70wQS0hItED6RH7uBn4FOgVq7/LRHnpfrhUxrfd/f3IAjNrTZCwAHD3AjNbA7Rx9w/DW1D3AT3N7B3gv9x9XSnbPpHgP9314R0jCD5U1kSss6ZkpQjF4gitAtpUsj5mNpogGXQIi+oDzcurE6EFwX/ScyLiN4L/rAvluPv+Sm6v0HTgEoIP0hkEifYGYD/wUXi+WwNb3X1XRL1VQEYZ23we+NTM7gAuA+a6+yoAM+sK/DmsW5fgg3hOifqRrRj3EpynSjOzugRXRMOBJmFxAzNLdvf8Equ3A1a5e14pmyr5O18Vxtuq5LLI92XE+rcD0919Wjnhngi8ZmYFEWX54T4mh/G9aGaNgWeBX7l7bjnbk6OgZ1BytNYR/BEDQQMEgj/atQDu/ry7nxGu4wS3uAinI60huIJq7u6Nw1dDd+8ZsU55Xe4XiyPUvjCOiuqb2YnA34A7CW4LNgYWEiSZ0pTc1mZgH9AzIv5GHjQqqUz8ZS2fDpxJcAU2HfgYGAIMDechOPamZtYgol7JYz+0E/fFBB/eFwLXEiSsQuOBr4Au7t6Q4FlhWefgaN1NcLUzMNzHWWF5aftZA7S30huVlPydtye4jbqx5LKS78vQ7eG2Hy4n1jXAhRG/08bunubua909193vd/cewOnASGB0OduSo6QEJUfrJeAiCxo/pBJ8+Bwg+A/9ZDM7x8xqE/zHv4/gv08IPkQ6FDYAcPf1wLvAQ2bWMHw43cnMhlYyjreArmZ2rZmlmNlVQA+C5yWVUY8gQeQAmNlNhI1AyrARaFvYgMDdCwgS3MNm1jLcRhszu6CS+z9sm+F2lxOct+uBGe6+M1zvO4QJyt3XENxq+n344L4XcAvwXDn7eh74EUFyeDmivAGwE9htZt2AO44g/spqQHBM28MGB/eWs+7nwHpgnJnVC49vSLjsBeAnYWOS+sDvgH+EV1tlvi8jtr2L4CruLDMbV8b+Hwd+G/4Dg5m1MLNR4fQwMzs1bHCyk+DWX8krQDkOlKDkqLj7UoIPz78QXEVcTNAc/SDBQ+VxYfkGgoYPvwyrFn4objGzueH0aILGFYsJGjj8k+Cef2Xi2ELwH+zdwBaCh+Ij3X1zJesvBh4ieG61ETgV+KScKh8QNBTZYGaF+/g5wYP4WWa2E3if4EqhskrbJgSJaIu7r46YN2BexDrXENyaXAe8Btzr7u+Vs68XCK7KPihxjn5KcFW1iyDh/uMI4q+s/yV4frmZoAHCf8paMbzldzHB88zVBLc6rwoXTyS4zTYD+Ibgn6AfhvXKe19Gbn87QUOKC83sgVJCeISg4cu7ZrYrjHdguOwEgvfoToJnVdMJbvPJcWbuGrBQRETij66gREQkLilBiYhIXFKCEhGRuKQEJSIicSkhvqjbvHlz79ChQ6zDEBGRUsyZM2ezu7coWZ4QCapDhw5kZmbGOgwRESmFmZXsDQbQLT4REYlTSlAiIhKXopqgzGy4BWMBZZnZ2FKWm5k9Gi5fYGb9wvI0M/vczL4Ix3W5P6JOUzN7z4Jxad4zsyYltysiItVf1J5Bhf1UPUbQnUg2MNvMpoRdyxS6EOgSvgYSdFg5kKDvrHPcfXfYn9bHZva2u88CxgJT3X1cmPTGEnQ1c0Ryc3PJzs5m//4j7WRaANLS0mjbti2pqamxDkVEaqhoNpIYAGS5+woAM3sRGEXQ31qhUcAzHvS3NMvMGptZetiB6O5wndTw5RF1zg6nJxEMQ3DECSo7O5sGDRrQoUMHIoZJkEpwd7Zs2UJ2djYdO3aMdTgiUkNF8xZfG4qPw5NN8TFZyl3HzJLNbD7BwHfvuftn4TqtwgRW2BN2y9J2bma3mlmmmWXm5OQctnz//v00a9ZMyekomBnNmjXT1aeIRFU0E1Rpn/wle6Ytcx13z3f3PgRDLg8ws/KGQDh8I+5PunuGu2e0aHFY8/pg50pOR03nTkRWLclhwcelthA/LqJ5iy+bYKCwQm0JhgQ4onXcfbuZfUgwfstCYGPhbUAzSye4whIRkSq0PWcPk387nT079rNx1XbOu7bXcf/HNZpXULOBLuGgYrWAqwnGV4k0BRgdtuYbBOwIE0+LcChlzKwOcB7BaJ+FdcaE02OAN6J4DHHv9ddfZ/HixRWvWMKUKVMYN66ssdpERMp2YH8ukx8MkhPArH8vY9e2fcd9P1G7gnL3PDO7E3gHSAYmuvsiM7s9XP44wWioIwgGe9sL3BRWTwcmhS0Bk4CX3L1whNRxwEtmdgvBQGZXROsYqoPXX3+dkSNH0qNHj0rXycvL45JLLuGSSy45ojopKQnR8YiIlMPd+efDM1m3YisASclJXPfLs2jYtO5x31dUP3Hc/S2CJBRZ9njEtAM/KKXeAqBvGdvcApx7fCONjZUrVzJ8+HAGDhzIvHnz6Nq1K8888wwzZ87kpz/9KXl5efTv35/x48dTu3Ztxo4dy5QpU0hJSeFb3/oWl112GVOmTGH69Ok8+OCDvPLKKwD84Ac/ICcnh7p16/K3v/2Nbt26ceONN9K0aVPmzZtHv379OPXUU8nMzOSvf/0rq1at4uabbyYnJ4cWLVrw1FNP0b59+8PqPPTQQzE+YyISa1Nf+JKFn64umh91R39OOqVVVPalf4mBX1wcvdGaf/+v68tdvnTpUiZMmMCQIUO4+eab+fOf/8wTTzzB1KlT6dq1K6NHj2b8+PGMHj2a1157ja+++gozY/v27TRu3JhLLrmEkSNHcvnllwNw7rnn8vjjj9OlSxc+++wzvv/97/PBBx8AsGzZMt5//32Sk5N5+umni2K48847GT16NGPGjGHixIn86Ec/4vXXXz+sjogktgUfr2LqCwuK5k+/uBsDLugStf2pq6MYa9euHUOGDAHg+uuvZ+rUqXTs2JGuXbsCMGbMGGbMmEHDhg1JS0vju9/9Lq+++ip16x5+Ob17924+/fRTrrjiCvr06cNtt93G+vXri5ZfccUVpSaamTNncu211wJwww038PHHH1dYR0QSS/byLfzz4ZlF8537pDPiln5R3aeuoGKssq1eUlJS+Pzzz5k6dSovvvgif/3rX4uujAoVFBTQuHFj5s+fX+o26tWrd8QxVbaOiNRcO7bsZfKD08k9mAdA89YNuebnZ5CcHN1rHCUoKr4NF02rV69m5syZDB48mBdeeIHzzjuPJ554gqysLDp37szkyZMZOnQou3fvZu/evYwYMYJBgwbRuXNnABo0aMCuXbsAaNiwIR07duTll1/miiuuwN1ZsGABvXv3LjeG008/nRdffJEbbriB5557jjPOOCPqxy0i1cPBA3k8+9vp7Ny6F4C0erUY/euzqVu/dtT3rVt8Mda9e3cmTZpEr1692Lp1Kz/5yU946qmnuOKKKzj11FNJSkri9ttvZ9euXYwcOZJevXoxdOhQHn74YQCuvvpq/vjHP9K3b1++/vprnnvuOSZMmEDv3r3p2bMnb7xRcSv8Rx99lKeeeopevXoxefJkHnnkkWgftohUA+7OP/93JtnLtwCQlJTEtWPPpEWbhlWyfwsa0tVsGRkZXnLAwiVLltC9e/cYRRRYuXIlI0eOZOHChTGN42jFwzkUkeh577kv+ODFL4vmR90xgEEjuh73/ZjZHHfPKFmuKygRETnMvA+/KZacBl10clSSU3mUoGKoQ4cO1fbqSURqrlVLcnjlkVlF8137tWbk906r8jgSOkElwu3NaNG5E6mZtm7czeTfTic/Lx+Alu0acfXPot9irzQJm6DS0tLYsmWLPmiPQuF4UGlpabEORUSOo317DjLp/mlFfezVa5jGmHuHUaderZjEk7DNzNu2bUt2djaljRUlFSscUVdEaoa8vHyeH/cRm9bsACA5JZkb7hlK01b1YxZTwiao1NRUjQYrIkJwV+RfT2SSNT+i55mfDObE7qWPpVdVEvYWn4iIBD5+fQmf/2d50fx51/Wm91kdYhdQSAlKRCSBLfx0NW8/Na9ovs/ZHTnnqiMawDxqlKBERBLUmmWbeemhT4sai3Xo0ZLv/GjQcR8Z92gpQYmIJKBtm3Yz6TcfFusA9oZ7hpKSGj+jFyhBiYgkmL27D/D0fYeak9dtUJsx955N3QbR7wD2SChBiYgkkLzcfJ773Yxizcmv/9VQmreumg5gj4QSlIhIgnB3Xv3LLFZ8ubGo7IqfDKZjz5YxjKpsSlAiIgni/ecXMG/aN0XzF4zuExfNycuiBCUikgBmv5tVrHfy/hd0YejlPWMYUcWUoEREarilc9bx+mOfF813Pa01o+7oHzfNycuiBCUiUoOt/Xorz4/7iIKCAgBan9SUa39+Zkx6Jz9S8R+hiIgclW2bdjPp/mkc3J8LQJOW9Rlz79nUrpMa48gqRwlKRKQG2rvrABN//QG7tu0DIK1eLW68bxgNm9aNcWSVF9UEZWbDzWypmWWZ2dhSlpuZPRouX2Bm/cLydmY2zcyWmNkiM/txRJ37zGytmc0PXyOieQwiItVN7sE8Jv3mQzav3QkcGjqjZbtGMY7syERtuA0zSwYeA84HsoHZZjbF3RdHrHYh0CV8DQTGhz/zgLvdfa6ZNQDmmNl7EXUfdvc/RSt2EZHqqqCggBf/+AmrvwrGujMzrrz7dE46pVWMIzty0byCGgBkufsKdz8IvAiMKrHOKOAZD8wCGptZuruvd/e5AO6+C1gCtIlirCIi1V7huE6LZ60pKhtxSz96nXFiDKM6etFMUG2ANRHz2RyeZCpcx8w6AH2BzyKK7wxvCU40syal7dzMbjWzTDPL1Ki5IpIIpr20kFlvLSuaP/PSHpwxqnsMIzo20UxQpTWw9yNZx8zqA68Ad7n7zrB4PNAJ6AOsBx4qbefu/qS7Z7h7RosWsR0VUkQk2ma/m8V7z35RNN/7rA5ceFPfGEZ07KKZoLKBdhHzbYF1lV3HzFIJktNz7v5q4QruvtHd8929APgbwa1EEZGE9dXstcW+iNu5TzqX/2Rw3H8RtyLRTFCzgS5m1tHMagFXA1NKrDMFGB225hsE7HD39Rac1QnAEnf/c2QFM0uPmL0UWBi9QxARiW+rluQc9kXc635xJikp8TOu09GKWis+d88zszuBd4BkYKK7LzKz28PljwNvASOALGAvcFNYfQhwA/Clmc0Py37p7m8BfzCzPgS3AlcCt0XrGERE4tnG1duLDTrYpFV9brxvGGl1a8U4suPDCof6rckyMjI8MzMz1mGIiBw323P2MP6/32Hnlr0A1GuUxh1/vIBm6Q1iHNmRM7M57p5Rslw9SYiIVDN7du5n4q8/KEpOtdJSufG+YdUyOZVHCUpEpBo5sD+XSfd/SE72oRFxb7hnKG07N4txZMefEpSISDWRl5vPs7+dwZplm4FDvUR07n1CjCOLDiUoEZFqoKCggJf+/ClZ89cXlY26o3+17SWiMpSgRETinLvzxvjZfPnxqqKy867rzcALu8YwquhTghIRiXPvTv6Cz/+zvGj+9Iu7cc5Vp8QwoqqhBCUiEsc+em0xH758qD+CvsM6MvJ7p1X7XiIqQwlKRCROZb6XxVsT5xbNd+vfhst+NCghkhMoQYmIxKUFH6/i1b8cGsSh4ymtuHZszejCqLKUoERE4szSOet46aFPKezpp/VJTbnhnqGk1opa73RxSQlKRCSOfLNoE8/9bgb5efkAtGjbiJt+cw516tWM/vWOhBKUiEicyF6+hUn3TzvU+WvL+tz8wDnUb5QW48hiQwlKRCQObFi1nafu/YAD+3IBqN+4Drc8eC6Nm9eLcWSxowQlIhJjm9ftZMI9U9m76wAAdRvU5pYHz61xnb8eKSUoEZEY2rZpNxPumcru7fsAqF0nlZvuP4cTTmwc28DigBKUiEiM7Ny6lwn3TGV7zh4AUmulMObeYbTtUvN6Jj8aSlAiIjGwe8d+/v6rqWxZvws4NGxGx54tYxxZ/FCCEhGpYnt3HWDi/0wtGtMpKTmJ635xJl36psc4sviiBCUiUoX27TnIU/d+wPpvtgHBmE5X/XQI3Qe0jXFk8UcJSkSkihzYl8tT935A9vItQJCcvvPjQTV6TKdjoQQlIlIFDuzP5en7prFm6eaism//YACnndsphlHFNyUoEZEoO3ggj8kPTGfl4k1FZZfc3p8BF3SJYVTxTwlKRCSKcg/mMfnB6Xy9YENR2YhbTmPwRSfHMKrqQQlKRCRKcg/m8exvZ5A1f31R2QVj+nLmt7vHMKrqQwlKRCQK8nLzee73H7Fs7rqism/d0IezL+8Zw6iql6gmKDMbbmZLzSzLzMaWstzM7NFw+QIz6xeWtzOzaWa2xMwWmdmPI+o0NbP3zGx5+LNJNI9BRORI5eXm8+zvZrA0c21R2XnX9mLYlafEMKrqJ2oJysySgceAC4EewDVm1qPEahcCXcLXrcD4sDwPuNvduwODgB9E1B0LTHX3LsDUcF5EJC4EV07Fk9Owq07lnKtPjWFU1VM0r6AGAFnuvsLdDwIvAqNKrDMKeMYDs4DGZpbu7uvdfS6Au+8ClgBtIupMCqcnAd+O4jGIiFRaXm4+z4/7iK9mRySnK0/h/Ot6YWYxjKx6imaCagOsiZjP5lCSqfQ6ZtYB6At8Fha1cvf1AOHPUjuuMrNbzSzTzDJzcnKO9hhERCql8MppyefZRWXDrjyF86/vreR0lKKZoEr7jfiRrGNm9YFXgLvcfeeR7Nzdn3T3DHfPaNGixZFUFRE5IrkH83j2dzOKXTkNvbynktMximaCygbaRcy3BdZVdh0zSyVITs+5+6sR62w0s/RwnXRgEyIiMVKYnIo9c7ryFC4Y3UfJ6RhFM0HNBrqYWUczqwVcDUwpsc4UYHTYmm8QsMPd11vwW50ALHH3P5dSZ0w4PQZ4I3qHICJStoMH8njmgeksm3Pof+9hV52qK6fjJCVaG3b3PDO7E3gHSAYmuvsiM7s9XP448BYwAsgC9gI3hdWHADcAX5rZ/LDsl+7+FjAOeMnMbgFWA1dE6xhERMpyYH8ukx8o3kPEOVefynnXqkHE8WLuJR8L1TwZGRmemZkZ6zBEpIY4sC+Xp++fxspFh54wnHddb85VU/KjYmZz3D2jZHnUrqBERGqifXsO8vR901j91aHWwReM6aseIqJACUpEpJL27jpQbDwnCDp+Vd960aEEJSJSCbt37GfCPVPZsHJbUdnFt2Zw+sXdYhhVzaYEJSJSgZ1b9zLhnqlsWrMDCEbCHfX9AQwcrvGcokkJSkSkHNs27WbCPVPZsn4XECSny+8aTL9zTopxZDWfEpSISBly1u5kwj1T2bF5DwBJyUlcdffp9DqzQ2wDSxBKUCIipVi/chsT/+cDdm/fB0BySjLXjj2DHgPbVVBTjhclKBGREtYs28xT905j3+4DAKTWTuGGe4bSpU96jCNLLEpQIiIRVny5kUm/+ZCD+3MBqF0nlTH3DqNjz1IHTpAoUoISEQl9NXstz4/7iNyDeQDUbVCbm35zDm07N4txZIlJCUpEBJg//RtefngmBfkFADRsWpdbHjyXlu0axTiyxKUEJSIJb9Zby5jy+GwK+yZt0qo+tzx4Ls1OaBDjyBKbEpSIJCx3Z9pLC3nv2S+Kylqd2Jibf3MODZvWjWFkAkpQIpKg3J1//30On0z5qqis3cnNufHeYdRtUDuGkUkhJSgRSTh5efm8+ugs5k37pqisc590rv/VWdROS41hZBJJCUpEEsrBA3k8P+6jYkO0nzKkPVfdPYSU1OQYRiYlKUGJSMLYu+sAT98/jTVLNxeVDRjehVF39CcpKSmGkUlplKBEJCFsz9nDxF9/QE72jqKyYVeewvnX99YQ7XFKCUpEaryNq7fz1L3Tijp9NTNGfu80jeUU55SgRKRGW7FwI5MfnM7+PQeBoNPXK/9rsHokrwaUoESkxlr46Wr+8adPyMvNB6BWWirX/+osdfpaTShBiUiNNPPfS/nXE5lFvUPUb1yHm+4fRuuTmsY4MqmscputmNnTEdNjoh6NiMgxcnfefmpusa6LmrduyB1/ukDJqZqpqF1l74jpH0czEBGRY5WXm88//vQJM15dXFTW7uTm3P7Hb9G0Vf0YRiZHo6IE5ceycTMbbmZLzSzLzMaWstzM7NFw+QIz6xexbKKZbTKzhSXq3Gdma81sfvgacSwxikjNsHf3ASb++gO+mLGyqKz7wLZ897fnUa9hWuwCk6NW0TOotmb2KGAR00Xc/UdlVTSzZOAx4HwgG5htZlPcfXHEahcCXcLXQGB8+BPgaeCvwDOlbP5hd/9TBbGLSILYunE3T983rdh3nAaN6MrFt2XoC7jVWEUJ6r8jpjOPcNsDgCx3XwFgZi8Co4DIBDUKeMaDG8WzzKyxmaW7+3p3n2FmHY5wnyKSYLKXb+Hp+6exZ8f+orILxvRl6Hd66Au41Vy5CcrdJx3DttsAayLmszl0dVTeOm2A9RVs+04zG02QNO92923HEKeIVFOLZq3hH3/8pGgEXH3HqWYpN0GZ2ZTylrv7JeVVL63KUaxT0njggXC9B4CHgJsP27nZrcCtAO3bt69gkyJSnbg7H7/xFW9PnFvUUq9ug9rccM9QOvRoGePo5Hip6BbfYIIrnBeAzyg9oZQlG2gXMd8WWHcU6xTj7hsLp83sb8CbZaz3JPAkQEZGxjE19hCR+JGfX8CbT2Yy661lRWXN0hsw5t5htGjTMIaRyfFW0dPDE4BfAqcAjxA0eNjs7tPdfXoFdWcDXcyso5nVAq4GSl6RTQFGh635BgE73L3c23tmFvkV8EuBhWWtKyI1y749B3nmNx8WS04ndm/BHX+6QMmpBqroGVQ+8B/gP2ZWG7gG+NDMfuPuf6mgbp6Z3Qm8AyQDE919kZndHi5/HHgLGAFkAXuBmwrrm9kLwNlAczPLBu519wnAH8ysD8EtvpXAbUd60CJS/WzduJtJ909j05pDLfV6n9WB7/x4EKm11ClOTWSF92/LXCFITBcRJKcOBFc9E919bXn14klGRoZnZh5pI0QRiRerluQw+bfTi7XUG3bVqZx/XS+11KsBzGyOu2eULK+okcQkgtt7bwP3u7tup4lIlZr7wQpe/ctn5OcFHb4mpyTznR8Pou/ZHWMcmURbRdfFNwB7gK7Aj82s8HLLAHd33fQVkahwd955Zj7T/7moqKxeozRu+NVQTuzeIoaRSVWp6BmUvoItIlXuwL5c/vHQJyz5LLuorNWJjRn9P2erT70EUtEtvjTgdqAzsIDg2VNeVQQmIolpy4ZdTH5gOhtXby8q69a/DVf9dAhpdWvFLjCpchXd4psE5AIfEbS264l6NReRKPl6wQaeH/cRe3cdKCo789IeDL+xj/rUS0AVJage7n4qgJlNAD6PfkgikmjcnZn/Xsa//z6HgvwCIGgMcdkPB9LvnJNiHJ3ESkUJKrdwIvxeU5TDEZFEk5ebzxvjZ5P5XlZRWYMmdbj+V0Npf3LzGEYmsVZRguptZjvDaQPqhPNqxScix2zn1r08+7sZrFm6uaisbZdmXPfLs2jcvF4MI5N4UFErvuSqCkREEsuqJTk8P+4jdm7dW1TWd1hHLr1zoHqGEKDiKygRkePus/8s519PZBZ9+dbMGHFLP4Zc0k09Q0gRJSgRqTK5B/P41xOZzH730POmug1qc/XPzqBLn/RyakoiUoISkSqxPWcPz/1+BtnLtxSVtT6pKdf/6iyatNSXb+VwSlAiEnXL56/nxT98XOz7Tb2HduCyHw6iVm19DEnp9M4Qkahxdz58eRHvPftF0ci3SclJXHRLPwaPPFnPm6RcSlAiEhV7dx/g5T9/ylezD43M06BJHa4de6aGZZdKUYISkeNu7ddbee73M9i2cXdRWYeeLbnmZ2fQsGndGEYm1YkSlIgcN+7O5+9k8eaTmeTl5heVn3lpD741ujcpKfpqpVSeEpSIHBcH9ufy+mOfM//Db4rKatdJ5fK7BnPK6e1jGJlUV0pQInLMNqzazvPjPiIne0dR2QkdmnDt2DNp0UY9osnRUYISkaPm7sx5/2umPJ5J7sFDQ8VlnN+ZS27PUJdFckz07hGRo3JgXy5vjP+cedMO3dJLrZ3CqDv6c9q5nWIYmdQUSlAicsTWfr2VF/7fR2xZv6uorGW7Rlz3i7No2a5RDCOTmkQJSkQqzd2Z+eZS3po4r6ijV4DTzuvExbdlUDstNYbRSU2jBCUilbJ7x35eeWRmsS/e1kpL5ds/GEDfszvGMDKpqZSgRKRCWV9s4OU/f1ps7KbWnZpyzc/OoHlrtdKT6EiK5sbNbLiZLTWzLDMbW8pyM7NHw+ULzKxfxLKJZrbJzBaWqNPUzN4zs+XhzybRPAaRRJaXl8/bT89j4v9MLZachozqzh1/vEDJSaIqagnKzJKBx4ALgR7ANWbWo8RqFwJdwtetwPiIZU8Dw0vZ9Fhgqrt3AaaG8yJynG3K3sH4n77DjFcWFXX0Wq9RGjfedw4jv3saKanqFUKiK5q3+AYAWe6+AsDMXgRGAYsj1hkFPOPBu3+WmTU2s3R3X+/uM8ysQynbHQWcHU5PAj4Efh6dQxBJPIXdFf37b3OKfbepS990rvjJ6TRoUieG0UkiiWaCagOsiZjPBgZWYp02wPpyttvK3dcDuPt6Myu1W2Qzu5Xgqoz27dXNikhl7Nq+j1cfnVWsIURySjLDb+yj4dilykUzQZX2TvajWOeouPuTwJMAGRkZx2WbIjXZ4s/W8OpfPmPPjv1FZa3aN+aq/x5Cegc96pWqF80ElQ20i5hvC6w7inVK2lh4G9DM0oFNxxypSALbv/cg//77XDLfyypWPuSSblwwpo+6K5KYieY7bzbQxcw6AmuBq4FrS6wzBbgzfD41ENhRePuuHFOAMcC48OcbxzVqkQSyYuFG/vnwTLZtOjRuU8Omdbn8J4Pp0ic9hpGJRDFBuXuemd0JvAMkAxPdfZGZ3R4ufxx4CxgBZAF7gZsK65vZCwSNIZqbWTZwr7tPIEhML5nZLcBq4IpoHYNITXXwQB7vTv6CT6d8VdRCD6DXmScy6o4B1G1QO4bRiQQs8s1ZU2VkZHhmZmaswxCJC6u+yuGf/zuTzWt3FpXVqV+bUXf0p/dZHWIXmCQsM5vj7hkly3VzWSRB5B7M4/3nv+SjVxcXu2rqelprLvvhIBo101DsEl+UoEQSwKqvcnjlkVnFBhSsXSeVi757Ghnnd1LzcYlLSlAiNdjBA3m8/9wCPn59SbGrpk69TuA7Px5Ek5b1YxidSPmUoERqqBVfbuTVv8wqNmZTrbRURtzSjwEXdNZVk8Q9JSiRGmbfnoP85+l5fP6f5cXKO/dJ57IfDtRVk1QbSlAiNciimauZ8nhmsZ7H0+rV4qJb+nHaeXrWJNWLEpRIDbBz617+9UQmCz9dXay8+8C2fPv7A2jYVC30pPpRghKpxgoKCvj8P1m8M2k++/ceLCqv37gOF9+WwalD2uuqSaotJSiRamr9N9t47bHPWLN0c7HyjPM7c+HNfalbX71BSPWmBCVSzRzYl8t7zy1g5r+WUlBQUFTevE1DLv3BQE46tVUMoxM5fpSgRKoJd2fRzDX868lMdm451AgiOSWZYVf2ZOjlPTXKrdQoSlAi1UDO2p28+WQmy+YWH42mU68TuOSO/rRs2yhGkYlEjxKUSBw7sD+X6S8vYsarS8jPyy8qr9cojYu+exp9hnZQIwipsZSgROKQu/PlJ6t5a8JcdmzeU1RuZgwc0ZXzr++lRhBS4ylBicSZDau28+aTmXy9YEOx8nYnN2fUHQNo06lpjCITqVpKUCJxYu+uA7z37Bd89vbyYh271muUxvAxfdQThCQcJSiRGMvLy+fzt5fz/vNfsm/3gaLypKQkBl98Mudecyp16tWKYYQisaEEJRIj7s7SzHX8e8KcYqPbQtCx68jvnUar9o1jE5xIHFCCEomBdSu28vZT88iav75YebP0Boy4pR/dB7TV7TxJeEpQIlVo++Y9vPfsF8z74Jtiz5nS6tZi2FWncPrFJ+vLtiIhJSiRKrBvz0FmvLKYT974ityDeUXlZsaA4V0477pe1G+UFsMIReKPEpRIFOXl5jPrrWVM+8dC9u46UGxZt/5tGH5jXz1nEimDEpRIFBQUFDBv2je8/9wCtufsKbas9UlNGXFLPzr1OiFG0YlUD0pQIsdRYYeu7z37BZvW7Ci2rEmr+lwwug+9zjxRDSBEKkEJSuQ4cHeWzV3Pu5Pns+7rrcWW1WuYxrCrTmHghV3UAELkCEQ1QZnZcOARIBn4u7uPK7HcwuUjgL3Aje4+t7y6ZnYf8D0gJ9zML939rWgeh0hZ3J2vF2zk/ee+YNWSnGLLatdJ5czLenDGqG7UrpMaowhFqq+oJSgzSwYeA84HsoHZZjbF3RdHrHYh0CV8DQTGAwMrUfdhd/9TtGIXqYwVX27k/ecX8M3CjcXKU2ulMOiirgy9vAf1GqplnsjRiuYV1AAgy91XAJjZi8AoIDJBjQKe8eALIbPMrLGZpQMdKlFXpMq5Oyu+3MjUF748LDElpyTT/4LODLuyJw2b1o1RhCI1RzQTVBtgTcR8NsFVUkXrtKlE3TvNbDSQCdzt7ttK7tzMbgVuBWjfvv1RHoJIwN1ZPm89H7z45WG38pKSk8g4vxNnX9GTJi3rxyhCkZonmgmqtGZKXsl1yqs7HnggnH8AeAi4+bCV3Z8EngTIyMgouV+RSnF3Fs9aw7SXFrE2a0uxZUnJSfQ75yTOufoUJSaRKIhmgsoG2kXMtwXWVXKdWmXVdfei+ypm9jfgzeMXskggLy+fBTNWMf2fiw5rLp6cksxp552kKyaRKItmgpoNdDGzjsBa4Grg2hLrTCG4XfciwS28He6+3sxyyqprZunuXtjD5qXAwigegySYA/tzyXz3az56bUmxkWwhaPzQ/4LOnHlZdxo3rxejCEUSR9QSlLvnmdmdwDsETcUnuvsiM7s9XP448BZBE/MsgmbmN5VXN9z0H8ysD8EtvpXAbdE6Bkkcu7btY+abS5n11vJiYzJB0Fx80EVdGTKqGw0a14lRhCKJxyJ7VK6pMjIyPDMzM9ZhSBzasGo7n7zxFfOmfUN+Xn6xZfUapTHkkm4MuqirBgwUiSIzm+PuGSXL1ZOEJJzCXh8+eWMJy+etP2x50xMacOal3TntvJNIraU/EZFY0V+fJIwD+3OZO3UFn/5r6WEj2AK069qcMy/rTs/B7UhKSopBhCISSQlKaryctTuZ9e9lzJ26gv17DxZbZmb0HNyOM77dnfbdmqsTV5E4ogQlNVJ+fgFffZ7NrLeWHzasOgQj2J52ficGj+xKsxMaxCBCEamIEpTUKNs37yHz3a+Z/U4WO7fuPWx58zYNOf3ik+l3zknqwFUkzilBSbWXl5fPsjnrmP1OFksz11GyZaqZ0W1AGwaPPJnOvU/QbTyRakIJSqqtnLU7mfP+18yduoJd2/Ydtrx+4zr0/1Yn+l/QWT0+iFRDSlBSrezbc5AvP17F3KkrDuu0tVCn3icwYHgXegxqS0qKBggUqa6UoCTu5eXlkzV/A/OnfcPiWdnkHsw7bJ0GTepw2nmdyDi/E83S1ehBpCZQgpK45O6sWbaFL6av5IsZK9mzY/9h6yQlJ9F9QBsyzu9Ml37pJCfru0siNYkSlMQNd2fDqu18MX0lC2asYtum3aWu1/qkpvQ9pyN9zu5I/UYasVakplKCkpgqTEoLP1nNgo9WldrDA0DDpnXpc3YH+gzrSHqHJlUcpYjEghKUVDl3Z23WVhZ+uppFn65h87rSk1JavVr0HNyOvsM60vGUlup+SCTBKEFJlcjLy+ebLzexeNYaFn+Wzc4th3+JFqBWWirdB7Sh11kn0rVfa1JS1QpPJFEpQUnU7N6xn6WZa/lq9lqWz13PgX25pa5XKy2Vbv1bc+oZJ3JyRmv1IC4igBKUHEcFBQVkL9/KsjnrWJq5lrVZWw/r1aFQ3Qa16TagDaec3p7OfU5QUhKRw+hTQY7Jtk27yZq/geXz1pM1f8Nho9FGatKyPj0GtaXHoHac2KOFmoWLSLmUoOSI7Ny6l28WbuLrBRv5+osNbN2wq8x1k5KSaN+tOSf3b0OPgW1p0bah+sETkUpTgpIyuTvbc/awctEmVi7OYcWXG8tsBl6oYdO6dO57AidntKFznxOoW792FUUrIjWNEpQUycvLZ8M321m1JIdVS3JY/dVmdmzeU26d1NopdDylJZ37pNO1Xzot2zXSVZKIHBdKUAnK3dm8bhdrs7aQvWwLq5duZv2KbeTl5pdbLyU1mXYnN+ekU1vRqfcJtDu5mTpkFZGoUIJKAHm5+Wxas4P132xj/YptrFuxjXVfby2z2XekWmmptO/WnA49W9LxlJa069pMLe5EpErok6YGyc8vYOuG3WxctZ1Na3awcdV2NqzczuZ1uyjIL6jUNpq0rE/77s1p360F7bs1J71jE7W2E5GYUIKqZgoKCti5ZR9b1u8KX7vJyd7B5rW72LphN/l55d+ii1SvURptuzSjTeemtOvanLZdm6nzVRGJG0pQccTd2b83l11b97Fjy152bA5e2zftYXvOHrZt3M32nL1HlIQKNWlVn/SOTWjdqSnpHRvTulNTGjWrqwYNIhK3opqgzGw48AiQDPzd3ceVWG7h8hHAXuBGd59bXl0zawr8A+gArASudPdt0TyOo+HuHNyfx/69uRzYm8u+3QfD1wH27DzA3p3Bzz079rN7e/DatW0/B/dX/FyoPA2b1aVV+0a0bN+YFm0bkt6xCa3aN6J2ndTjdGQiIlUjagnKzJKBx4DzgWxgtplNcffFEatdCHQJXwOB8cDACuqOBaa6+zgzGxvO/zwax7B98x7+/bc5ABQUeNBtjwfPegrynfy8AvJy88nPDX7mHsjn4IE8cg/kk3sgr8xufo5VvUZpNEtvQLP0+jRLb0DzNg3DVwNqpykRiUjNEM0rqAFAlruvADCzF4FRQGSCGgU848En+Swza2xm6QRXR2XVHQWcHdafBHxIlBLUwX15LPx0dTQ2XabUWik0bFaHhs3q0rBpHRq3rEfjFsGrSat6NG5ZT0lIRBJCNBNUG2BNxHw2wVVSReu0qaBuK3dfD+Du682sZWk7N7NbgVsB2rdvf3RHcIyPZ2qlpZJWN5VadVKoU79W8KpXizoNalO/UW3qNUqjbsPa1G+cRoMmdajfOI20uql6LiQiQnQTVGmfsiXveZW1TmXqlsvdnwSeBMjIyDiqe22NmtXlmp+fSVJSEI4lGWZBH3PJKUkkJRspqckkpyaRWiuZlNRkaqWlUCsthdTayRpgT0TkGEQzQWUD7SLm2wLrKrlOrXLqbjSz9PDqKR3YdFyjjlC7Tiq9zjgxWpsXEZFyRPNf/NlAFzPraGa1gKuBKSXWmQKMtsAgYEd4+668ulOAMeH0GOCNKB6DiIjESNSuoNw9z8zuBN4haCo+0d0Xmdnt4fLHgbcImphnETQzv6m8uuGmxwEvmdktwGrgimgdg4iIxI5Fqyl0PMnIyPDMzMxYhyEiIqUwsznunlGyXE/xRUQkLilBiYhIXFKCEhGRuKQEJSIicSkhGkmYWQ6w6hg20RzYfJzCqc50HnQOCuk8BHQejs85ONHdW5QsTIgEdazMLLO0FiaJRudB56CQzkNA5yG650C3+EREJC4pQYmISFxSgqqcJ2MdQJzQedA5KKTzENB5iOI50DMoERGJS7qCEhGRuKQEJSIicUkJqgJmNtzMlppZlpmNjXU8VcHM2pnZNDNbYmaLzOzHYXlTM3vPzJaHP5vEOtaqYGbJZjbPzN4M5xPqPJhZYzP7p5l9Fb4nBifaOQAws5+Efw8LzewFM0tLhPNgZhPNbJOZLYwoK/O4zewX4eflUjO74Fj2rQRVDjNLBh4DLgR6ANeYWY/YRlUl8oC73b07MAj4QXjcY4Gp7t4FmBrOJ4IfA0si5hPtPDwC/MfduwG9Cc5FQp0DM2sD/AjIcPdTCIYBuprEOA9PA8NLlJV63OHnxNVAz7DO/4Wfo0dFCap8A4Asd1/h7geBF4FRMY4p6tx9vbvPDad3EXwgtSE49knhapOAb8ckwCpkZm2Bi4C/RxQnzHkws4bAWcAEAHc/6O7bSaBzECEFqGNmKUBdglG+a/x5cPcZwNYSxWUd9yjgRXc/4O7fEIz1N+Bo960EVb42wJqI+eywLGGYWQegL/AZ0Coc8ZjwZ8sYhlZV/hf4GVAQUZZI5+EkIAd4KrzN+Xczq0dinQPcfS3wJ4JBUtcTjP79Lgl2HiKUddzH9TNTCap8VkpZwrTLN7P6wCvAXe6+M9bxVDUzGwlscvc5sY4lhlKAfsB4d+8L7KFm3sYqV/iMZRTQEWgN1DOz62MbVVw6rp+ZSlDlywbaRcy3Jbisr/HMLJUgOT3n7q+GxRvNLD1cng5silV8VWQIcImZrSS4vXuOmT1LYp2HbCDb3T8L5/9JkLAS6RwAnAd84+457p4LvAqcTuKdh0JlHfdx/cxUgirfbKCLmXU0s1oED/+mxDimqDMzI3jmsMTd/xyxaAowJpweA7xR1bFVJXf/hbu3dfcOBL/7D9z9ehLoPLj7BmCNmZ0cFp0LLCaBzkFoNTDIzOqGfx/nEjybTbTzUKis454CXG1mtc2sI9AF+Pxod6KeJCpgZiMInkMkAxPd/bexjSj6zOwM4CPgSw49e/klwXOol4D2BH+wV7h7yYenNZKZnQ381N1HmlkzEug8mFkfgkYitYAVwE0E/9wmzDkAMLP7gasIWrnOA74L1KeGnwczewE4m2BYjY3AvcDrlHHcZvYr4GaC83SXu7991PtWghIRkXikW3wiIhKXlKBERCQuKUGJiEhcUoISEZG4pAQlIiJxKSXWAYgkGjPLJ2jCb0A+cKe7f3qM29zt7vXLWf574B2gMdDN3ccdy/5EqoKuoESq3j537+PuvYFfAL+vgn0OJPge21CC77iJxD0lKJHYaghsg6AHDzP7Yzje0JdmdlVY/n9mdkk4/ZqZTQynbzGzB8vbeLi9BUB/YCbBl0vHm9mvo3hMIseFbvGJVL06ZjYfSAPSgXPC8suAPgRjLjUHZpvZDGAGcCZBNzJtwjoAZxD0EVgmd/9vM3sZuAH4L+BDdx9yPA9GJFp0BSVS9Qpv8XUjGNTtmbB/tzOAF9w93903AtMJrnw+As4MB4NbzKGOOgcDlXl21ReYD3QL64tUC7qCEokhd59pZs2BFpQ+VAHuvjYc7mE4wdVUU+BKYHc4oGSpwj70niboUXozwSB7Fl69DXb3fcfvSESOP11BicSQmXUj6Ih4C0HyucrMks2sBcFItoU9Qc8E7grX+Qj4KRU0dnD3+e7eB1gG9AA+AC4Ir96UnCTu6QpKpOoVPoOC4KppjLvnm9lrBLftviAY5O1n4XAXECSjb7l7lpmtIriKqrA1Xpjotrl7gZl1c3fd4pNqQ72Zi4hIXNItPhERiUtKUCIiEpeUoEREJC4pQYmISFxSghIRkbikBCUiInFJCUpEROLS/wehQtVxQB7ljQAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"posterior2.plot(label='posterior', color='C4')\n",
"decorate_bowls('Posterior after two vanilla cookies')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After two vanilla cookies, the high-numbered bowls have the highest posterior probabilities because they contain the most vanilla cookies; the low-numbered bowls have the lowest probabilities.\n",
"\n",
"But suppose we draw again and get a chocolate cookie.\n",
"Here's the update:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.645875Z",
"iopub.status.busy": "2021-04-16T19:35:13.645126Z",
"iopub.status.idle": "2021-04-16T19:35:13.650340Z",
"shell.execute_reply": "2021-04-16T19:35:13.649962Z"
},
"tags": [
"hide-output"
]
},
"outputs": [
{
"data": {
"text/plain": [
"0.2462686567164179"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likelihood_chocolate = 1 - hypos/100\n",
"\n",
"posterior3 = posterior2 * likelihood_chocolate\n",
"posterior3.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's the posterior distribution."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.684232Z",
"iopub.status.busy": "2021-04-16T19:35:13.671117Z",
"iopub.status.idle": "2021-04-16T19:35:13.798448Z",
"shell.execute_reply": "2021-04-16T19:35:13.797816Z"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABE0ElEQVR4nO3dd3gc1dX48e+RZEm2XOQiV9mWbEvuXbg3isEGYwPBYBKwCSHEL5AGhAB505MXfqEkIRAI1XTbdCehG3DvTe623OUqy72on98fM1qtZFla2VrNSjqf59nHe2fmzp5Zr/bs3Llzr6gqxhhjTKgJ8zoAY4wxpjSWoIwxxoQkS1DGGGNCkiUoY4wxIckSlDHGmJBkCcoYY0xIsgRlQp6IfE9EvqjC1/uTiBwWkQNV9ZpeEJHnReTX7vNRIpLut26niFzhXXTFicjtIjK/il/zWxG5sypf0xRnCcoU434xnRWRUyJyUEReFZH6F7G/34nImxcTk6q+papXXsw+AiUibYH7gW6q2jIYX4wico2IzBeRYyJyQEReFJEGlfkagVDVqar6x2C+hoi0EpFZIrJPRFREEoL5el4JtYReU1iCMqW5VlXrA/2AS4D/9SoQEYm4iLoiIhX9jLcHMlX10IW+bokYSou/EfAnoDXQFYgHHq+M1wtBBcBnwHe8DsRUP5agzHmp6l7gU6AHgIiMF5H17i//b0Wka+G2IvJLEdkrIidFZLOIXC4iY4BHgJvdM7I17raNRORlEdnv1vmTiIS7624XkQUi8lcROQL8ruRZjIgMEZFlInLc/XeI37pvReTPIrIAOAN0KHlcIvKQiGxzY90gIte7y68AvgRau/HOAJ4HBrvlY+52USLyhIjsds8ynxeRuu66USKS7r4fB4BXS3lf31bVz1T1jKoeBV4Ehpb2f+DG+l6JZX8Xkafd598XkY3usWwXkR/5bVcYy/0icsh9v7/vt36aiPyptNct8XoDRGSR+/++X0SeEZHI8uq5x3pQVf8JLAtkexFpKyIfiEiGiGSKyDMl1j8hIkdFZIeIjPVb3to9UzsiImki8kO/deEi8ojf//kK90y5zM9SidftKCJfuzEdFpG3RCTWXfcG0A74t/s5edBdPkhEFrrv2xoRGRXIe2D8qKo97OF7ADuBK9znbYH1wB+BZOA0MBqoAzwIpAGRQGdgD9DarZcAdHSf/w54s8RrfAT8C4gBmgNLgR+5624H8oAfAxFAXXfZfHd9E+AocJu7/ha33NRd/y2wG+jurq9TyjFOxDl7CQNudo+rlbtuFJDut63vtf2W/Q2Y5cbSAPg38Khf/Tzg/wFRQN0A3vO/AdPPs649TqJt6JbDgf3AILd8DdAREGCku22/ErH8wf0/u9pd39hdPw3403mO2/9z0B8Y5L6fCcBG4GcV/FxFAAoklLFNOLAG+Kv72YgGhvn9P+QCP3S3+x9gHyDu+jnAP906fYAM4HJ33S+AtTifUwF6A00D/Czd6T7vhPPZjwLigLnA30p7v9xyGyDTfc/D3LqZQJzXf+PV6eF5APYIrYf7h3YKOAbscv/o6wK/Bmb6bRcG7HW/2DoBh4ArKJEQKJGggBZANn5f3O4Xwzfu89uB3SX2cTtFCeo2YGmJ9YuA293n3wJ/qOAxrwYmuM9HUUaCcr/gTuMmYHfZYGCHX/0cIDrA1x7tfikml7HNfGCy3/bbytj2I+CnfrGcBSL81h+iKLlNI4AEVcpr/Az4sILvcSAJajBOYokoZd3tQJpfuZ67v5Y4P6TygQZ+6x8FprnPNxf+/5bYZyCfpTvPE+t1wKrzvV/AL4E3StT5HJhSkfettj8uuH3f1GjXqepX/gtEpDVOwgJAVQtEZA/QRlW/FZGf4SSj7iLyOXCfqu4rZd/tcX7N7xeRwmVhOGdghfaUrOSnWByuXTi/WAOpj4hMBu7DORsAqA80K6uOnzicL8cVfvELzq/6QhmqmlXejkRkEPA2cKOqbilj07dxkvjrwHfdcuE+xgK/xTnDDXNjW+tXN1NV8/zKZ3CON2Aikgw8BaS4+48AVlRkHwFqC+wqEa8/X69KVT3jvv/1cc6GjqjqSb9td7nxFu53Wyn7C+SzBICINAeeBobjnDWH4fywOJ/2wEQRudZvWR3gmzLqmBLsGpQJ1D6cPzrA6YCA84e/F3zXVYa52yhOExfuc397cM6gmqlqrPtoqKrd/bYpa4j9YnG42hXGUV59EWmPc83nXpymnFhgHU6SKU3JfR3GOSvp7hd/I3U6lQQSf2EcfXGaCe9Q1dnlbP4uMEpE4oHrcROUiEQB7wNPAC3cY/mkjGO5UM8Bm4AkVW2Ic12xsl8DnM9GO6l4x5h9QBMp3hPS/zOxB6cZtLR65X2WCj2K8//ay30PbqX4e1Da5/wNv89IrKrGqOpjAR2RASxBmcDNBK4Rp/NDHZyu2NnAQhHpLCKXuV+YWThf4PluvYNAgri96VR1P/AF8KSINBSRMPcC9MgA4/gESBaR74pIhIjcDHQD/hNg/RicL5MMcDoZ4HYCOY+DQHxhpwBVLcBJcH91f1UjIm1E5KoAXx8R6YHTs+3Hqvrv8rZX1Qyc5qZXcZoSN7qrInGuiWQAee7ZVDC64zcATgCnRKQLzvUfH3E6pvzufJVFJNqNEyDKLZdmKc71tcdEJEZEokWk1M4j/lR1D7AQeNSt0wv4AfCWu8lLwB9FJEkcvUSkKRX7LDXAbfoWkTY417X8HaR4h5w3gWtF5Cq3k0a022klvrzjMUUsQZmAqOpmnF+N/8A5i7gWpzt6Ds6Xz2Pu8gM4HR8ecau+6/6bKSIr3eeTcb5cN+A0k7wHtAowjkxgHE6CzMTprDFOVQ8HWH8D8CTOtYaDQE9gQRlVvsbpKHJARApf45c4HUQWi8gJ4CucC/CBuh+nqfBlt9fXKRFZX06dt3Gu8fma99wmrZ/g/Hg4itP8N6sCcQTqAXffJ3GS84wS69tS9nt4FufLHZwzsbOlbaSq+Tifq044HV3ScTqxBOIWnCbbfcCHwG9V9Ut33VM479EXOIn2ZZxroBX5LP0e57aL48B/gQ9KrH8U+F+3x94DbtKcgPN3kIFzRvUL7Du3Qgp7wBhjTIW5ZwTvqupgr2MxNY8lKGOMMSHJTjeNMcaEJEtQxhhjQpIlKGOMMSGpVt+o26xZM01ISPA6DGOMqdVWrFhxWFXjSi6v1QkqISGB5cuXex2GMcbUaiJSckQPwJr4jDHGhChLUMYYY0KSJShjjDEhqVZfgzLGmIrIzc0lPT2drKxyB6s3pYiOjiY+Pp46deoEtL0lKGOMCVB6ejoNGjQgISEBv+lWTABUlczMTNLT00lMTAyojjXxGWNMgLKysmjatKklpwsgIjRt2rRCZ592BmWMCQnZWbnknM2joEApyC9AwoS69SOJjIoIqYQQSrFUNxV97yxBGWOqTG5OHvu2HeXArmMc2HmMQ7uPcfzwGU4ezSInK7fUOuER4cQ0jKJp6wbExTckLr4hbZOb0SapCRER4aXWMTVDUBOUiIwB/o4zHfZLJWeTdGdl/TtwNc5U1Ler6kp33Ss4c7UcUtUefnVmUDT3TixwTFX7iEgCsBHY7K5brKpTg3RoxpgAqCp7046wecU+tqceYPemw+Tl5pdf0U9+Xj4njpzhxJEz7Fh30Le8TmQE7bo0o1PfVvQc1o6mLRuUsRfj76OPPiI5OZlu3bpVqN6sWbPYsGEDDz30UJAiKy5oCUpEwoFngdE4E48tE5FZ7oRxhcYCSe5jIM7U0gPdddOAZ4DX/ferqr4JzETkSZwJxAptU9U+lXogxpgKUVX2bMkkdd4u1i/czbGM0wHVC48IJzqmDuHhYUiYUJBfwJmTOeTnlZ7QcnPy2JZ6gG2pB/j8tVXEJzWl1/D29L+iI/UaRJVaxzg++ugjxo0bV6EElZeXx/jx4xk/fnyF6kREXHiaCeYZ1AAgTVW3A4jIdJwZJv0T1ATgdXUmpVosIrEi0kpV96vqXPesqFTu2ddNwGVBOwJjTMDOnMpm9Tc7WPp5Ggd3HTvvds1aN6RNUhNato+lZUJjmrSsT/3G0dSNiTznGoWqkpuTz8mjZzm89yQZ6cfZv+MouzZkkLn/ZLFt07dmkr41ky/fTKXPpYkMHd+ZFu1ig3Ck3tq5cydjxoxh4MCBrFq1iuTkZF5//XUWLVrEAw88QF5eHpdccgnPPfccUVFRPPTQQ8yaNYuIiAiuvPJKbrjhBmbNmsWcOXP405/+xPvvvw/APffcQ0ZGBvXq1ePFF1+kS5cu3H777TRp0oRVq1bRr18/evbsyfLly3nmmWfYtWsXd9xxBxkZGcTFxfHqq6/Srl27c+o8+eSTF3yswUxQbXCmOS6UTtHZUVnbtAH2B7D/4cBBVd3qtyxRRFbhTOv8v6o6r2QlEbkLuAugXbt2AbyMMaYsRw+dYt6HG1n+xTZyc/LOWV+3fhRdB7YhqW8rEnu0oFHTegHvW0SIjIqgacsGNG3ZgM79W/vWHTt8mm1rDrBuwW62rjrgO9PKzclj2edbWfb5VroPbsuVt/WhedtGF3+gJTx87ZuVvs9Cj/771jLXb968mZdffpmhQ4dyxx138NRTT/Gvf/2L2bNnk5yczOTJk3nuueeYPHkyH374IZs2bUJEOHbsGLGxsYwfP55x48Zx4403AnD55Zfz/PPPk5SUxJIlS7j77rv5+uuvAdiyZQtfffUV4eHhTJs2zRfDvffey+TJk5kyZQqvvPIKP/nJT/joo4/OqXMxgpmgSuuuUXL63kC2OZ9bgHf8yvuBdqqaKSL9gY9EpLuqnii2c9UXgBcAUlJSbDphYy5Qxt4TfDNjLWvm7KKgoKDYujpREfQa3p7eIxJI7Nk8KJ0ZYpvF0P/yjvS/vCNnT+ewfuFuFv1nC/u2H/Fts37RHjYsTqfvZYmMvrU3sc1iKj0OL7Rt25ahQ4cCcOutt/LHP/6RxMREkpOTAZgyZQrPPvss9957L9HR0dx5551cc801jBs37px9nTp1ioULFzJx4kTfsuzsbN/ziRMnlppoFi1axAcffADAbbfdxoMPPlhunYoKZoJKB9r6leOBfRewzTlEJAK4AehfuExVs4Fs9/kKEdkGJAM2XLkxleh45hlmv7OWFV9uOycxtUpszMCxSfQemUB0vcgqi6luTCQpozvR/4qO7NyQwYKPN7J+kdM4o6qsnL2ddQv2MPrWXgy5tjNhYdX7FtBAu2tHRESwdOlSZs+ezfTp03nmmWd8Z0aFCgoKiI2NZfXq1aXuIyYmsKTuH1OgdcoTzAS1DEgSkURgLzAJ+G6JbWYB97rXpwYCx1U1kOa9K4BNqppeuEBE4oAjqpovIh1wOl5sr4TjMMbg3Kc05931zP9o0zlNeR17tWTkxO506t3S0/uERITE7s1J7N6cvduO8MXrq9my0vnNm5OVy39fWsHqb3dyw48H0rpDk4t6rfKa4YJp9+7dLFq0iMGDB/POO+9wxRVX8K9//Yu0tDQ6derEG2+8wciRIzl16hRnzpzh6quvZtCgQXTq1AmABg0acPKkcw2vYcOGJCYm8u677zJx4kRUldTUVHr37l1mDEOGDGH69OncdtttvPXWWwwbNqzSjzNoPyNUNQ+4F/gcp/v3TFVdLyJTRaSw+/cnOEkkDXgRuLuwvoi8AywCOotIuoj8wG/3kyjevAcwAkgVkTXAe8BUVT2CMeaiqCrrFu7mb3f/h29mriuWnDr2bsn/PDGGO/98BUl9WoXUTaxtOjbh+7+/jDv/fEWxzhJ70zJ59r7PmPvBBpz+WdVP165dee211+jVqxdHjhzh5z//Oa+++ioTJ06kZ8+ehIWFMXXqVE6ePMm4cePo1asXI0eO5K9//SsAkyZN4vHHH6dv375s27aNt956i5dffpnevXvTvXt3Pv7443JjePrpp3n11Vfp1asXb7zxBn//+98r/Tiluv4HVYaUlBS1CQuNOb9jGaf58NklbFlRvOW9TaemXDWlD0l9WnkUWcXk5eUz74ONfD19bbH7sJL7tebGnw+mQWzdgPazceNGunbtGqwwA7Jz507GjRvHunXrPI3jQpX2HorIClVNKbmtjSRhjDmHqrL08zQ+eXllsREeYhpFM+b2vvS/vENInS2VJyIinEtv6kHPYe2Y+dRC9mw+DMCWlft4+sefcOsjI2jf9ZwZx43HqveVQmNMpTt2+DQv/+9sPnp2iS85iQiDrk7mvuevJeWKjtUqOflr1rohdz02mpE3dvctO3XsLC8+8hWrvqkel6wTEhKq7dlTRdkZlDHGZ+2CXXzwjyVknc7xLYuLb8R3fjqI9l1qxhlGREQ4Y6b0pWPvlsx4fAGnT2SRn5fPzKcWcmjPCa68rXeZCVhVq22C9lpFLynZGZQxhuysXD74x2LefmyeLzmJCCNv7M6P/z62xiQnf0l9WnH3U2OKdaD49t11zHhiAXnnGV4pOjqazMzMatu5wkuF80FFR0cHXMc6SVgnCVPLZew9wVuPzi02PFHj5vW5+YGhteK6TNaZHKY/voDNy/f6lnUdEM8tvxxGncjijUw2o+7FOd+MuufrJGEJyhKUqcXWL9rNe39bTNaZoia93iMSmHD3AOrGVN2Ntl4rKChg1vPLWfLpFt+yjr1actuvRxIVHdj05ObCnS9BWROfMbWQqvLFG6t58//m+pJTRJ1wvvOTQdz8wNBalZwAwsLCmPA/lxTrPLEt9QCv/uZrss8zT5UJPktQxtQy2Vm5vPHnOXwzs6gnWOPm9Zn6+FWkjO5UazsAiAhjpvTlqil9fct2bczgzT/PLXUQXBN8lqCMqUWOHjrF87/4go1LfKOEkdyvNff+bSxtOl7c0D81xagbu3PNnb5hPklbvZ93/jL/vB0nTPBYgjKmlkhPy+Sf93/OgZ1HfcuGXdeVKb8dZRP8lTBsQldG31o0Ft3GJem897dF1nuvitl9UMbUAptX7OPtx+b5brwNjwjnursvIWV0J48jC12X3tSD7LN5zH1/PQBr5uykSYv6XHlbH28Dq0UsQRlTwy3/Mo0Pn1nqmxqjbv0obv3VCDr0aOFxZKHNuSbVh+wzub7efd/MXEeTlvUtsVcRa+IzpoZSVb59dx3vP73Yl5waN6/P1L9cackpQCLCtT9KoXNKG9+yD59dStqaAx5GVXtYgjKmBlJVPpu2is9fX+1b1rpDE6Y+fmVQpj+vycLDw7jlwWG0SmwMQEF+AW89OpeMvSfKqWkuliUoY2qYgoICPnx2KXM/2OBb1rFXS+56bDQNm9TzMLLqK6puHSb/ZpTv/cs6ncNb/zfX7pEKMktQxtQg+fkFzHxyIcs+3+pb1nVgPFN+O4qoujYiwsWIbRbD5N+M8g1/dHD3MT54erH17AsiS1DG1BB5efnMeGIBa+bu9C3rd1kHvvfwiHPGlDMXpk3HJky4+xJfOXXeLhb+e7OHEdVslqCMqQHy8vKZ/pf5rJ2/y7ds0NXJ3PizwYSH2595Zep/eUcGjk32lT95ZSU7NxzyMKKayz65xlRzeXn5vP3YPNYv2uNbNuTaLoyfekmtHbYo2Mb9sD9tk5sBTqeJGU8s4KzfHFqmcgQ1QYnIGBHZLCJpIvJQKetFRJ5216eKSD+/da+IyCERWVeizu9EZK+IrHYfV/ute9jd12YRuSqYx2ZMKMh3vxz9hy4afn03xv2wvyWnIIqoE853Hx5O3frOCBzHMk7z0bNL7HpUJQtaghKRcOBZYCzQDbhFRLqV2GwskOQ+7gKe81s3DRhznt3/VVX7uI9P3NfrBkwCurv1/unGYEyNVFBQwLtPLWTdgt2+ZcOv78bY7/e15FQFYpvFcMOPB/rKqfN2sfLr6jFtfHURzDOoAUCaqm5X1RxgOjChxDYTgNfVsRiIFZFWAKo6FzhSgdebAExX1WxV3QGkuTEYU+OoKu8/vbhYh4ih47tYcqpiPYa045Iri0aVmPX8cjL3n/QwopolmAmqDbDHr5zuLqvoNqW5120SfEVEGl/kvoypVlSV/7y4gpWzi36tD7o6mWvutGY9L1zzw/40a9MQgJysXGY+tdA3coe5OMFMUKX9pZRsoA1km5KeAzoCfYD9wJMV2ZeI3CUiy0VkeUZGRjkvZUzo+ertVBb+e5OvfMmVnaxDhIeioutw8wNDCY9wrijs3pRhXc8rSTATVDrQ1q8cD+y7gG2KUdWDqpqvqgXAixQ14wW0L1V9QVVTVDUlLi4uoAMxJlTM+3ADX09f6yv3Gt6e6+4ZYMnJY/GdmnLpTUWz8X7x+hpr6qsEwUxQy4AkEUkUkUicDgyzSmwzC5js9uYbBBxX1f1l7bTwGpXreqCwl98sYJKIRIlIIk7Hi6WVcSDGhIIVs7fxySsrfeXOKW2YeN8QwsLsbpFQMHJid1omOFcccnPyeN9GmbhoQftkq2oecC/wObARmKmq60VkqohMdTf7BNiO06HhReDuwvoi8g6wCOgsIuki8gN31V9EZK2IpAKXAj93X289MBPYAHwG3KOqNgWmqRE2LdvLB08v8ZUTe7Tgew8PJyLCOqqGioiIcG782WDfD4Yd6w6y+JMtHkdVvUltzvApKSm6fPlyr8Mwpky7Nx/mpUe+IjcnD3BGJb/z/66gbkykx5GZ0nzxxmq+mek07ERG1+Hnz40jtlmMx1GFNhFZoaopJZdb24AxIexQ+nFe+/03vuTUuEV9bv/dpZacQthlk3r6pjTJycrlk5dXllPDnI8lKGNC1MmjZ5n22284czIbgJhG0dzxh8to0Liux5GZskTUCWfC3UW3YK6dv4stK8vs+2XOwxKUMSEoOyuX1/7wLUcPnQKcpqLbf3spzVo39DgyE4gOPVrQ77IOvvKs55f5zoJN4CxBGRNi8vMLmP6X+exNywScacdv+eUw4pOaehyZqYixd/Qj2m2Kzdx/kjnvbSinhinJEpQxIURV+c8Ly9m0bK9v2YS7B9AlxQZFqW7qN4pmzJQ+vvK3764n84DdG1URlqCMCSEL/725WNfkURN7MHBMkocRmYtxyVWdfGe++Xn5fPbqKo8jql4sQRkTIjYt28t/X1rhK/ca3p4rb+vtYUTmYoWFhXHtj4pm4F23cDfb1x30MKLqxRKUMSFg/46jvPOX+b6RB9p2bsaNPxtsQxjVAO06N6P3iARf+ZOXV9oIEwGyBGWMx04ePctrf/iWnKxcABo3r89t/zuSOpERHkdmKsuY2/v6/j/3pmXavFEBsgRljIdyc/J4489zOH74NADR9SKZ8ttRNIi1e51qkti4GIZd39VX/uL1NWS7P0jM+VmCMsYjqspHzy5lz+bDgNOdfNKDw2jRLtbbwExQjLyxGw2b1APgxJEzzP9oUzk1jCUoYzwy/6ONxZp6rv5BPzr3b+1hRCaYoqLrMPrWXr7yvA82cPpElocRhT5LUMZ4YPOKfXzq1+U4ZXQnho7v4mFEpir0vayDb5y+7LO5fPvueo8jCm2WoIypYhl7TzD98aIeewndmjPhf2xG3NogPDyM0bcW3Tqw+L9bOOZefzTnsgRlTBXKOpPDm3+eQ9bpHAAaNYvhe4+MIKKOzetUW3Qf3NZ3825ebj6z315bTo3ayxKUMVVEVZn51EIO7TkOQJ3ICG791QjqN4r2ODJTlUSEMbf39ZVXfLWNQ+nHPYwodFmCMqaKfD19LRuXpPvKN/xkIPGdbADY2qhjr5Z06tMKcH64zH471eOIQpMlKGOqwKZle/nK70to+PXd6DMy0cOIjNeumtzH93zt/N2+M2tTxBKUMUF2eN8JZjy5wFfu2LslV/mNcm1qp/ikpnS5xBmlXlX5erpdiyopqAlKRMaIyGYRSRORh0pZLyLytLs+VUT6+a17RUQOici6EnUeF5FN7vYfikisuzxBRM6KyGr38Xwwj82YQGRn5fLm/831dYqIjYvhlgeHER5uvw0NXH5L0X1RqfN22bWoEoL2VyIi4cCzwFigG3CLiHQrsdlYIMl93AU857duGjCmlF1/CfRQ1V7AFuBhv3XbVLWP+5haKQdizAVSVT78xxIO7joGOFOBf++REcQ0tE4RxhGf1JRk9+ZsVeXbmevKqVG7BPNn3AAgTVW3q2oOMB2YUGKbCcDr6lgMxIpIKwBVnQscKblTVf1CVQvnTl4MxAftCIy5CIv+s5k1c3f6ytfdM8A6RZhz+J9Frf52Jxl7T3gYTWgJZoJqA+zxK6e7yyq6TVnuAD71KyeKyCoRmSMiwysSrDGVadfGDP778kpfeeDYZPpf3tHDiEyoate5GUl9i3r02VlUkWAmqNJuiy85CUog25S+c5FfAXnAW+6i/UA7Ve0L3Ae8LSINS6l3l4gsF5HlGRkZgbyUMRVy8thZ3n5sHgX5BYDTjDPuh/09jsqEsmJnUXN2cvTQKQ+jCR3BTFDpQFu/cjyw7wK2OYeITAHGAd9Td7wYVc1W1Uz3+QpgG5Bcsq6qvqCqKaqaEhcXV4HDMaZ8+fkFTP/LfE4cOQNAvQZRfO9hGynClK191zg69GwBQEF+gY107gpmgloGJIlIoohEApOAWSW2mQVMdnvzDQKOq+r+snYqImOAXwLjVfWM3/I4t2MGItIBp+OFzQpmqtTst1PZvtaZ0ltEuPmBocTGxXgclakORnynu+/5si/SOHMy28NoQkPQEpTbkeFe4HNgIzBTVdeLyFQRKexh9wlOEkkDXgTuLqwvIu8Ai4DOIpIuIj9wVz0DNAC+LNGdfASQKiJrgPeAqap6TicLY4Jl07K9fON3/eDyW3qS3M+mzzCBSe7XilaJjQHIzc5j0X82exyR96RwROXaKCUlRZcvX+51GKYGOHroFP/46aecPeX86k3q24rbf3cpYWF2v5MJ3Oo5O5jxhHNTd70GUTz4ynVERdfxOKrgE5EVqppScrn99RhzkfLy8nnn/833JadGzWK46f6hlpxMhfUc1p7GzesDcOZkNiu+3OZxRN6yvyBjLtJnr65izxZn2vaw8DBu+eUwG6HcXJDw8DCGXd/VV57/0Sby3d6gtZElKGMuwrqFu1kwq6jH1Zjb+9K+i/UONRcuZXRH32gjRw+dYv2iPeXUqLksQRlzgTIPnOT9pxf7yl0HxjNsgk3bbi5OZFQEA69O8pUXzqq9Xc4tQRlzAfJy85n+l/m+QWAbN6/PjT8bbNO2m0ox6OpkwiOce+d2bczwNSHXNpagjLkAn01bRfrWTADCI8K55ZfDqFc/yuOoTE3RoHFdeo9o7ysvqKVnUZagjKmg9Yv3FPvCGPv9vrRNbuZhRKYmGjK+qLl47fzdHM88U8bWNZMlKGMq4OihU7z3t0W+cteB8Qy5trOHEZmaqk3HJiT2KBr+aMknWzyOqOpZgjImQHl5+bzjd90pNi6GG39q151M8AwdX/TjZ8mnW8nJzitj65rHEpQxAfrqrVT2bHbvdwoLY9KDw6jXwK47meDpOjC+2I27qX7zi9UGlqCMCcDWVfuZ8956X3n0bb3tficTdGFhYQweVzQpw6L/bqE2DU9nCcqYcpw8epYZTy7wlZP7tWbkd7p5GJGpTfpf0ZE6kREA7Nt2xNd7tDawBGVMGVSVmU8t5PTxLMDp/jvxviF23clUmXoNoug5vJ2vvLgWdZawBGVMGea8t5601c4UZSLCTfcPtXH2TJUbdHVRM9/aebtrzVxRlqCMOY9dmzL48s1UX3nUxO506t3Sw4hMbRWf1JTWHZsAkJuTx4qvasco55agjCnF2dM5zHh8AQUFzkjS7bvGcfl3e3kclamtRITB1xSdRS35dGut6CxhCcqYElSVD59ZwtFDpwCIjonk5geGEh5ufy7GO71GJBAdEwlA5v6TpK0+4HFEwWd/ccaUsPzLbaydv8tXvuHHA333ohjjlcioCPpf3sFXXvrZVg+jqRqWoIzxc2jPcf79r+W+8oAxSfQc2r6MGsZUnQFjiqbh2LAknVNu79KaKqgJSkTGiMhmEUkTkYdKWS8i8rS7PlVE+vmte0VEDonIuhJ1mojIlyKy1f23sd+6h919bRaRq4J5bKbmyc3JY/rj88nNcYaTad62Edfc2d/jqIwp0rxtI9p3dW4QL8gvYNXX2z2OKLiClqBEJBx4FhgLdANuEZGSdzeOBZLcx13Ac37rpgFjStn1Q8BsVU0CZrtl3H1PArq79f7pxmBMQD6btpr9O44CEFEnnEkPDiMyKsLjqIwpLuXKTr7nSz9Pq9GdJYJ5BjUASFPV7aqaA0wHJpTYZgLwujoWA7Ei0gpAVecCR0rZ7wTgNff5a8B1fsunq2q2qu4A0twYjCnXpmV7Wfhvvyk07uhHq4TGZdQwxhs9h7Ujqm4dAA7vPcGujRkeRxQ8wUxQbYA9fuV0d1lFtympharuB3D/bX4R+zKGE0fOFJ9CY0B8sS69xoSSqOg69B6Z4Csv+yLNu2CCLJgJqrSxYEqeiwayTWW+HiJyl4gsF5HlGRk195eHCYyq8u5fF3H6hHOxuWGTenznp4NsKCMT0i7xa+ZbO383Z90pYGqaYCaodKCtXzke2HcB25R0sLAZ0P33UEX2paovqGqKqqbExdlo1LXd/I82lhjKaAgxDW0oIxPa2nRqQku3CTo3O4/UebvKqVE9BTNBLQOSRCRRRCJxOjDMKrHNLGCy25tvEHC8sPmuDLOAKe7zKcDHfssniUiUiCTidLxYWhkHYmqm9LRMPn99ja884jvd6NjLhjIyoU9EuOTKjr7y8hrazBe0BKWqecC9wOfARmCmqq4XkakiMtXd7BNgO06HhheBuwvri8g7wCKgs4iki8gP3FWPAaNFZCsw2i2jquuBmcAG4DPgHlXND9bxmeotOyuXGY8vID/P+YjEJzXliu/ZUEam+ugzKpHwCKejcvrWTA7tOe5xRJVPanIXxfKkpKTo8uXLy9/Q1DjvP72Y5V86vzojo+vwk6evpmmrBh5HZUzFvPXoXNYt3A3AyBu7M2ZKX48jujAiskJVU0ouL/MMSkSm+T2fUsamxlQbqfN3+ZITwPipKZacTLXU97JE3/NV3+zwDW5cU5TXxNfb7/lPgxmIMVXhWMZpPnxmia/ce0QC/S7rUEYNY0JXcv/WxLjzk53IPMO21IMeR1S5yktQtbf9z9Q4BQUFzHhyAVlul9zGzesz4e4B1qXcVFsREeH08bsnqqYNfVRegop3x8r7h99z36MqAjSmsnz77np2rnfuSggLC+PmXwylrjt9gTHVVV+/FoB1C/eQfTbXw2gqV3kDjf3C77n1JjDV1q5NGcx+e62vfNktPWnfxe6DM9Vf6w6NadE+loO7jpGbnce6hbvpf3nH8itWA2UmKFV9raz1xlQHWWdymPFE0ey4Cd2ac+lN3T2OypjKISL0u6wDn766EoCVs7fXjgQlIiVvrC1GVcdXbjjGVL5Zzy/j6MGi2XFvun8IYWE2FZqpOfqMSuCzaatQVbavPcixw6eJbRbjdVgXrbwmvsE4A7C+Ayyh9PHujAlZq77dwapvdvjK199rs+Oamqdhk3p07N3SN2xX6txdjLih5OxG1U95PyNbAo8APYC/44zccFhV56jqnGAHZ8zFyDxwko//WTTaVb/LO9BrmM2Oa2qmPqMSfM/XzNnpWRyVqcwEpar5qvqZqk4BBuEMSfStiPy4SqIz5gLl5eUz44kFvh5NTVs1YPyPLvE4KmOCp9ugtkTUcYY+2rf9SI0Y+qjchnh38NUbgDeBe4CngQ+CHZgxF+ObGevYs/kwAGHhYUz6xTDfJG/G1ER1YyLpcknRFHhr5u70LphKUt5QR68BC4F+wO9V9RJV/aOq7q2S6Iy5ANvXHeSbGet85dG39iY+qamHERlTNfwnMlwzZ2e1nw6+vDOo24BknGGOFonICfdxUkROBD88YyrmzKlsZj650PeH2bFXS0Z+p/pfLDYmEJ1TWhNdz7n5PHP/SdK3Znoc0cUp7xpUmKo28Hs0dB8NVLVhVQVpTCBUlQ+fWcLxw6cBqNcgion3DbGhjEytUScygu5DiuZtXTO3ek9kWF4TX7SI/ExEnnGnSi+vW7oxnlnx1TbWLdjtK9/wk0E0alrPw4iMqXr+zXypc3dW6xHOy2view1IAdYCVwNPBj0iYy5Axt4TzPpX0WhcA8cm031Q2zJqGFMzdezVgvqxdQE4efQsO9Yd8jiiC1deguqmqreq6r+AG4HhVRCTMRWSl5fP9Mfnk5udB0Dzto24+gf9PI7KGG+EhYXRa3jR/X5r5+8uY+vQVl6C8g2L607hbkzI+eL1NezbdgSAiDrhTHpwGJFR1hptaq8eQ9v5nq9buJv8/OrZzFfuhIX+PfeAXtaLz4SSLSv3Me/DDb7y2O/3pVVCYw8jMsZ7Cd3iaNjEuf56+niWb5qZ6qa8XnzhJXruRVSkF5+IjBGRzSKSJiIPlbJe3Lml0kQkVUT6lVdXRGaIyGr3sVNEVrvLE0TkrN+65yv0Tphq5+Sxs8x8aqGv3OWSNgwe19nDiIwJDSJS7CwqdV717M0XtCGdRSQceBYYC3QDbhGRkjekjAWS3MddwHPl1VXVm1W1j6r2Ad6n+KgW2wrXqerUYB2b8Z6q8t7fFnH6eBYADRrX5Ts/HWxdyo1x+Seo9Yv2VMtmvmDOOTAASFPV7aqaA0wHJpTYZgLwujoWA7Ei0iqQuuJ8E92EM9K6qWUWzNrElhX7fOWJ9w2hfqNoDyMyJrSUbOarjr35gpmg2uBM1VEo3V0WyDaB1B0OHFTVrX7LEkVklYjMERHrcVhD7d12hM+mrfaVR9zQjaQ+rbwLyJgQVLKZb+386tfMF8wEVVpbS8mBoc63TSB1b6H42dN+oJ2q9gXuA94WkXOuk7k3HC8XkeUZGRnnDd6EpuysXKb/ZT75efkAxCc1ZfRtvT2OypjQ1HOYXzPfwurXzBfMBJUO+N8pGQ/sC3CbMuu6I1rcAMwoXKaq2aqa6T5fAWzDGUewGFV9QVVTVDUlLi7uAg7LeGnW88s4vM/pQBoZXYdJvxhGRES4x1EZE5rad/Vr5jtR/Zr5gpmglgFJIpIoIpHAJKDkFPKzgMlub75BwHFV3R9A3SuATaqaXrhAROLczhWISAecjhfbg3VwpuqtnrODlbOL/kuvu/sSmrZq4GFExoS2ks186xdWr5t2g5ag3Bt77wU+BzYCM1V1vYhMFZHCHnaf4CSRNOBF4O6y6vrtfhLndo4YAaSKyBrgPWCqqh4JysGZKpe5/yQfPuM3O+5lHeh7aQcPIzKmevAfPHbD4vRqNQVHUG+3V9VPcJKQ/7Ln/Z4rziSIAdX1W3d7Kcvex+l2bmqYvLx83vnLfHKy/GbHnWqz4xoTiIRuccQ0jOb0iSxOHDnD7s2Had+lelzeCGYTnzGV4vPXVrM3zZnXJjzCGcrIZsc1JjBhYWF0HRjvK29YtKeMrUOLJSgT0jYt28v8jzb6ymNu70N8J5sd15iK6D64qJlv/aI91aaZzxKUCVnHDp/m3b8WH8po6PguHkZkTPXUsXcLIqOdVofM/Sc5sOuYtwEFyBKUCUn5+QXMfHIhZ05mA9CwST0bysiYC1QnMoLOKa195Q2L08vYOnRYgjIh6evpa9mx7iDgdJW9+RdDbSgjYy5CD//efNXkOpQlKBNytqUe4JsZ63zly7/biw49WngYkTHVX3L/1oS7N7Xv236EIwdPeRxR+SxBmZBy8thZZjyxwHcRt2Ovllx6U3ePozKm+ouuF0mnPi195fXV4CzKEpQJGQUFznWnk0fPAhDTKJqb7h9CWJh9TI2pDP69+TYstgRlTMC+fXc9aav3A851p5vuG+IbR8wYc/G6DGjj62i0a0OGrxNSqLIEZULC9nUH+eqtVF951MTuJPdrXUYNY0xFNYitS3yycx+hqrJ5+V6PIyqbJSjjuVPHs5jxeNF1p4Tuzbn8u708jsqYmqmb36gSG5dagjLmvAoKCpj++HxOHDkDQEzDaCb9Yhjh4fbRNCYYugwoSlBbVuwjz51bLRTZt4Dx1NfT17FtzQHAue408b4hNGpq152MCZYW7RrRuEV9ALLP5ob0HFGWoIxntq7az9fT1/rKo27qQef+dt3JmGASkeLNfEtCd1QJS1DGE8czz5xzv9MV3+3pcVTG1A7+zXwbl4TuHFGWoEyVy8vL553/N4/TJ7IAaNC4Ljc/MNTudzKmiiR0jyM6JhKAYxmnQ3bwWPtGMFXus2mr2bUxA3Dmqpn04DAaNK7rcVTG1B4REeEk92vlK4dqM58lKFOlUufvYsHHRfM7XTm5t42zZ4wHulaD61CWoEyVOZR+nPf/vthX7jaoLSNu6OZhRMbUXsn9W/ua1dO3ZnLy2FmPIzqXJShTJbLP5vLW/80lJysXgKatGnDjz2x+J2O8Uq9+FO26NvOVt67c72E0pQtqghKRMSKyWUTSROShUtaLiDztrk8VkX7l1RWR34nIXhFZ7T6u9lv3sLv9ZhG5KpjHZgKnqrz390Uc2nMccCZP+97DI6jrXqQ1xnijc0ob3/NQHPYoaAlKRMKBZ4GxQDfgFhEp2Z4zFkhyH3cBzwVY96+q2sd9fOLW6QZMAroDY4B/uvsxHpv/0UbWLdjtK19/7wBaJTb2MCJjDFBslt0tK/eTn1/gYTTnCuYZ1AAgTVW3q2oOMB2YUGKbCcDr6lgMxIpIqwDrljQBmK6q2aq6A0hz92M8tC31AJ++uspXHnRNZ/pe2sHDiIwxhVq2j6VRsxgAsk7nsHvTYY8jKi6YCaoN4D/hSLq7LJBtyqt7r9sk+IqIFP4UD+T1EJG7RGS5iCzPyMioyPGYCjp66BTv/L/5vpsA23eN45o7+5VTyxhTVUSk2FlUqDXzBTNBlXb1u+Ttyufbpqy6zwEdgT7AfuDJCrweqvqCqqaoakpcXFwpVUxlyMnO480/z/XdjFs/ti7ffWg4ERHW6mpMKCmeoPZ5GMm5gpmg0oG2fuV4oOTRn2+b89ZV1YOqmq+qBcCLFDXjBfJ6pgqoKh8+s4R9248AEBYexvceHm6TDxoTgjr2bkm4+8PxwM6jHDt82uOIigQzQS0DkkQkUUQicTowzCqxzSxgstubbxBwXFX3l1XXvUZV6Hpgnd++JolIlIgk4nS8WBqsgzPnN//jTaz+doevfO1dKSR0a+5hRMaY84mKrkNij6K/zy0rQud3fdASlKrmAfcCnwMbgZmqul5EporIVHezT4DtOB0aXgTuLquuW+cvIrJWRFKBS4Gfu3XWAzOBDcBnwD2qGroTndRQW1fv59NXVvrKl1zZiYFjkzyMyBhTnlBt5pNQHcW2KqSkpOjy5cu9DqPGOLzvBM/e9xlZp3MAaNu5GXc9OpqIOnbdyZhQlrH3BE9NdRq4IqPr8Ou3b6zSv1sRWaGqKSWX20gSplKcPZ3D63+c40tODZvU49ZHRlhyMqYaaNa6AU1bNQAgJyvXN5iz1yxBmYtWUFDAjCcWkJFeNFLEbf870jpFGFNNiAhJfYsu729dFRrDHlmCMhfts2mri90/ccNPBhKf1NTDiIwxFZXkN/1GqIzLZwnKXJSln29l3ocbfOWRN3anz8hEDyMyxlyIDj1bEBbupIR924+ExOjmlqDMBduWeoCPn1vmK3cb1JarJvfxLiBjzAWLrhdJ+65FgxekhUAznyUoc0Ey9p7grUfnUeAOLtm6QxNuun+ITZ9hTDXmfx1qSwg081mCMhV26ngWr/3+G86eygacHnu3/XokUdF1PI7MGHMxil2HWrUfr29DsgRlKiQ3J483/jSHzP0ngaIee7HuiMjGmOqrTccmxDSKBuD08Sz2bT/qaTyWoEzAVJV3/7qI3ZuceyREhEkPDrUee8bUECJCpz4tfWWvu5tbgjIB+2zaKtbO3+UrX3Nnf7oNbFtGDWNMdZPcr2jYo60rvR32yBKUCciCWZuY+0FRd/Ih13Zh6PguHkZkjAkG/zOoXRsPk52V61kslqBMuVLn7+K/L63wlbsOjLeJB42poRo2qUfLBGce2Py8fLanHvQsFktQpkzbUg8w88mFxWbFveXBYYSF2UfHmJoq2a8337Y1BzyLw75lzHnt236EN/88l/w8Z9aSuPhGTP71KOpERngcmTEmmDr1KUpQaastQZkQk7H3BK/85muyzhSNTv79319KvQZRHkdmjAm29t3ifLPsHtx9jBNHzngShyUoc45jh0/zyq9nc/p4FgDRMZHc/vtLady8vseRGWOqQmRUBO27NvOVvWrmswRlijl1PItXf/M1xzJOA1AnKoLbf3cprdyLpsaY2qFYM58lKOO1M6eyeeXXszm0x5nXKTwinFsfGUH7LnHl1DTG1DQdexd1N09bfcCTYY8sQRkAss7kMO2337B/hzO0iYhw0/1Dit20Z4ypPdp0akJ0TCQAJzLPkLH3RJXHENQEJSJjRGSziKSJyEOlrBcRedpdnyoi/cqrKyKPi8gmd/sPRSTWXZ4gImdFZLX7eD6Yx1aTZGfl8trvv2XPlsO+Zd/56SB6DWvvYVTGGC+Fh4fRoWcLX9mL3nxBS1AiEg48C4wFugG3iEi3EpuNBZLcx13AcwHU/RLooaq9gC3Aw37726aqfdzH1OAcWc1SmJx2bjjkW3bd3QPof3lHD6MyxoQC/1El0lZX/bh8wTyDGgCkqep2Vc0BpgMTSmwzAXhdHYuBWBFpVVZdVf1CVfPc+ouB+CAeQ42WnZXLtN99w451RXeKX3NnfwaOTfYwKmNMqPC/DrVj7SHy3fnfqkowE1QbYI9fOd1dFsg2gdQFuAP41K+cKCKrRGSOiAwvLSgRuUtElovI8oyMjMCOpAYqTE471xedOV19Rz+GTejqYVTGmFAS16YhDZvWA5zr1HvTjlTp6wczQZU2tWrJbiDn26bcuiLyKyAPeMtdtB9op6p9gfuAt0Wk4Tk7UX1BVVNUNSUurnb2Tjt7OodXfv118eT0g/4Mv75kC6wxpjYrOf1GVTfzBTNBpQP+czHEAyXHbj/fNmXWFZEpwDjge+r2fVTVbFXNdJ+vALYB1lZVwqnjWbz0yFe+OZ3AadYbfp2dORljztXJr5mvqgeODWaCWgYkiUiiiEQCk4BZJbaZBUx2e/MNAo6r6v6y6orIGOCXwHhV9Y2/ISJxbucKRKQDTseL7UE8vmrneOYZXnjoS/ZtLzpNv/auFGvWM8acV6JfT77dmw6Tl5tfZa8dtFE/VTVPRO4FPgfCgVdUdb2ITHXXPw98AlwNpAFngO+XVdfd9TNAFPCliAAsdnvsjQD+ICJ5QD4wVVWrtsE0hGXsPcGrv/mao4dOAc6p+w0/GUTKFdZbzxhzfrHNYmjWuiGH950gNyeP3ZsP06FHi/IrVoKgDkutqp/gJCH/Zc/7PVfgnkDruss7nWf794H3Lybemmr35sO89vtvOHMyG4Cw8DBufmCo3edkjAlIh14tOLzPuVF3e+rBKktQNpJEDbdp2V5eeuQrX3KqExXB5F+PsuRkjAmY/w2729dW3XUom9inBlv8yRZmPb/MN4ZWTMNopvx2FG2Tm5VT0xhjinQocR0qJzuPyKjgpw87g6qBCgoK+PcLy/n4uaW+5NS4eX1+9JcrLTkZYyqsQeO6NG/bCHCmgd+1sWruIbUEVcOcPZ3DG3+aw8J/b/Iti09qyv88cRVxbc65LcwYYwLifxa1o4qa+SxB1SCH9hznn/d/xqZle33Legxtxw8fHU2DxnU9jMwYU9116FWUoLZV0f1Qdg2qhli/aDczn1pETlaub9moiT248rbeuN3xjTHmgvmfQaVvzSQ7K5eo6DpBfU1LUNVcXl4+X76xhrkfbPAtqxMZwXd+OojeIxK8C8wYU6PENIymZUJjDuw8SkF+ATvXZ9C5f3Dni7MEVY0dOXiK6Y/PZ8/monmcGreoz22/GkmrRJui3RhTuTr2asGBnc6kptvXHrQEZUqXOm8nHz67lKzTOb5lyf1bc/P9Q6nXIMrDyIwxNVWHXi1YMMvpgLU9NfgTGFqCqmbOnMzm4+eWkjpvl29ZWHgYV03uw/Dru9r1JmNM0CR0b46IoKrs23aU7LO5RNUN3nUoS1DVyIYle/jo2aWcPHrWt6xx8/pMenAY7Trb/U3GmOCqVz+KFu1jnetQBQXs3nSYpL6tgvZ6lqCqgWOHT/OfF5azftGeYstTRnfi6h/0o25MpEeRGWNqmw49i65D7Vh30BJUbZWXl8+i/2zhq7dSi3Ufb9C4Ljf8eBBdLiltkmFjjAmehO5xvoEAdqw7VM7WF8cSVAhSVdYv2sNn01aRuf9ksXUpozsx9vt9rSOEMcYTiX4jme/ZkkluTh51IoOTSixBhRBVZfvag3z55ppzxrpq3rYR190zkMTuzT2KzhhjoH6jaOLiG5GRfpz8vHz2bM4sdhNvZbIEFQJUlS0r9/PNjLXnJKa69aO49OYeDB6XTEREuEcRGmNMkcQezclIPw4416EsQdVA2Vm5rPp6B4v+s5lDe44XWxceEc6ga5K57OYe1pxnjAkpiT2as/SzrUBwr0NZgqpiqsqOdYdYM3cnqXN3kXUmp9j68IhwUkZ3ZOSN3WjcvL5HURpjzPkl+F1q2L3pMHl5+UFp4bEEVQWys3LZsfYQW1buY/2iPZzIPHPONpHRdUgZ3ZHhN3QltlmMB1EaY0xgYpvF0KRlA44cOEluTh57tx6hfde4Sn+doCYoERkD/B0IB15S1cdKrBd3/dXAGeB2VV1ZVl0RaQLMABKAncBNqnrUXfcw8AMgH/iJqn4ezOMrTfbZXI4cPMWBHUfZm3aE9K2ZpG89Qn5efqnbN2vTkMHXJNP3sg52P5MxptpI7NGcIwecXsY71h+qXglKRMKBZ4HRQDqwTERmqeoGv83GAknuYyDwHDCwnLoPAbNV9TERecgt/1JEugGTgO5Aa+ArEUlW1dIzw0U4fSKL/760gpysPHKz88nJyiP7bC7HD5/hzMnscuvHNIym5/D29BmZQLsuzWx4ImNMtZPYozkrvtoGOB0lRt3YvdJfI5hnUAOANFXdDiAi04EJgH+CmgC8rs685ItFJFZEWuGcHZ2v7gRglFv/NeBb4Jfu8umqmg3sEJE0N4ZFlX1gBQXKqm92VKhOy4TGJPdrRVK/1iR0j7MeecaYas3/fqhdGzLIzy8gPLxy58ANZoJqA/iPzZOOc5ZU3jZtyqnbQlX3A6jqfhEpvFrXBlhcyr6KEZG7gLsA2rVrV4HDKRIZff63LTwinNi4esTFN6RNp6a06dSE+KSmNqOtMaZGadw8hkbNYjh++DTZZ3M5uOsYrTs0qdTXCGaCKq3dSgPcJpC6F/J6qOoLwAsAKSkp5e2zVJFREdx03xDqREW4j3AioyJo2LQuDRrXtSY7Y0yNJyJcNqkH4RFhdOjZgti4yu/cFcwElQ609SvHA/sC3CayjLoHRaSVe/bUCijshB/I61UKEaHvpR2CsWtjjKk2BlyVFNT9V26DYXHLgCQRSRSRSJwODLNKbDMLmCyOQcBxt/murLqzgCnu8ynAx37LJ4lIlIgk4nS8WBqsgzPGGBNcQTuDUtU8EbkX+Bynq/grqrpeRKa6658HPsHpYp6G0838+2XVdXf9GDBTRH4A7AYmunXWi8hMnI4UecA9wejBZ4wxpmqI04GudkpJSdHly5d7HYYxxtRqIrJCVVNKLg9mE58xxhhzwSxBGWOMCUmWoIwxxoQkS1DGGGNCUq3uJCEiGcCui9hFM+BwJYVTXdX296C2Hz/Ye2DHf/HH315VzxlttlYnqIslIstL63lSm9T296C2Hz/Ye2DHH7zjtyY+Y4wxIckSlDHGmJBkCerivOB1ACGgtr8Htf34wd4DO/4gsWtQxhhjQpKdQRljjAlJlqCMMcaEJEtQF0hExojIZhFJE5GHvI4n2ESkrYh8IyIbRWS9iPzUXd5ERL4Uka3uv429jjWYRCRcRFaJyH/ccm07/lgReU9ENrmfhcG16T0QkZ+7n/91IvKOiETX9OMXkVdE5JCIrPNbdt5jFpGH3e/FzSJy1cW8tiWoCyAi4cCzwFigG3CLiHTzNqqgywPuV9WuwCDgHveYHwJmq2oSMNst12Q/BTb6lWvb8f8d+ExVuwC9cd6LWvEeiEgb4CdAiqr2wJkKaBI1//inAWNKLCv1mN3vhElAd7fOP93vywtiCerCDADSVHW7quYA04EJHscUVKq6X1VXus9P4nwxtcE57tfczV4DrvMkwCogIvHANcBLfotr0/E3BEYALwOoao6qHqMWvQc4c+jVFZEIoB7OrN01+vhVdS5wpMTi8x3zBGC6qmar6g6cuf4GXOhrW4K6MG2APX7ldHdZrSAiCUBfYAnQwp0FGfff5h6GFmx/Ax4ECvyW1abj7wBkAK+6zZwviUgMteQ9UNW9wBM4E6Xux5kB/AtqyfGXcL5jrtTvRktQF0ZKWVYr+uuLSH3gfeBnqnrC63iqioiMAw6p6gqvY/FQBNAPeE5V+wKnqXnNWeflXmeZACQCrYEYEbnV26hCTqV+N1qCujDpQFu/cjzOqX6NJiJ1cJLTW6r6gbv4oIi0cte3Ag55FV+QDQXGi8hOnCbdy0TkTWrP8YPzuU9X1SVu+T2chFVb3oMrgB2qmqGqucAHwBBqz/H7O98xV+p3oyWoC7MMSBKRRBGJxLkoOMvjmIJKRATn2sNGVX3Kb9UsYIr7fArwcVXHVhVU9WFVjVfVBJz/769V9VZqyfEDqOoBYI+IdHYXXQ5soPa8B7uBQSJSz/17uBznWmxtOX5/5zvmWcAkEYkSkUQgCVh6oS9iI0lcIBG5GueaRDjwiqr+2duIgktEhgHzgLUUXYN5BOc61EygHc4f8ERVLXlBtUYRkVHAA6o6TkSaUouOX0T64HQSiQS2A9/H+aFbK94DEfk9cDNOr9ZVwJ1AfWrw8YvIO8AonGk1DgK/BT7iPMcsIr8C7sB5j36mqp9e8GtbgjLGGBOKrInPGGNMSLIEZYwxJiRZgjLGGBOSLEEZY4wJSZagjDHGhKQIrwMwpjYTkXycrvsC5AP3qurCi9znKVWtX8b6R4HPgVigi6o+djGvZ0yw2BmUMd46q6p9VLU38DDwaBW85kCc+9dG4tzbZkxIsgRlTOhoCBwFZ+QOEXncnXdorYjc7C7/p4iMd59/KCKvuM9/ICJ/Kmvn7v5SgUuARTg3mT4nIr8J4jEZc8Gsic8Yb9UVkdVANNAKuMxdfgPQB2fOpWbAMhGZC8wFhuMMKdPGrQMwDGeMwPNS1V+IyLvAbcB9wLeqOrQyD8aYymRnUMZ4q7CJrwvOBG+vu+O8DQPeUdV8VT0IzME585kHDHcnhttA0aCdg4FArl31BVYDXdz6xoQsO4MyJkSo6iIRaQbEUfq0BajqXnfahzE4Z1NNgJuAU+5EkqVyx9CbhjO69GGcyfbEPXsbrKpnK+9IjKkcdgZlTIgQkS44gw9n4iSfm0UkXETicGayLRwVehHwM3ebecADlNPZQVVXq2ofYAvQDfgauMo9e7PkZEKSnUEZ463Ca1DgnDVNUdV8EfkQp9luDc6Ebw+6012Ak4yuVNU0EdmFcxZVbm88N9EdVdUCEemiqtbEZ0KajWZujDEmJFkTnzHGmJBkCcoYY0xIsgRljDEmJFmCMsYYE5IsQRljjAlJlqCMMcaEJEtQxhhjQtL/BwUf+ta8q8jmAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"posterior3.plot(label='posterior', color='C4')\n",
"decorate_bowls('Posterior after 2 vanilla, 1 chocolate')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now Bowl 100 has been eliminated because it contains no chocolate cookies.\n",
"But the high-numbered bowls are still more likely than the low-numbered bowls, because we have seen more vanilla cookies than chocolate.\n",
"\n",
"In fact, the peak of the posterior distribution is at Bowl 67, which corresponds to the fraction of vanilla cookies in the data we've observed, $2/3$.\n",
"\n",
"The quantity with the highest posterior probability is called the **MAP**, which stands for \"maximum a posteriori probability\", where \"a posteriori\" is unnecessary Latin for \"posterior\".\n",
"\n",
"To compute the MAP, we can use the `Series` method `idxmax`:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.802724Z",
"iopub.status.busy": "2021-04-16T19:35:13.802249Z",
"iopub.status.idle": "2021-04-16T19:35:13.806979Z",
"shell.execute_reply": "2021-04-16T19:35:13.806610Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"67"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior3.idxmax()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or `Pmf` provides a more memorable name for the same thing:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.811076Z",
"iopub.status.busy": "2021-04-16T19:35:13.810496Z",
"iopub.status.idle": "2021-04-16T19:35:13.815123Z",
"shell.execute_reply": "2021-04-16T19:35:13.815492Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"67"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior3.max_prob()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you might suspect, this example isn't really about bowls; it's about estimating proportions.\n",
"Imagine that you have one bowl of cookies.\n",
"You don't know what fraction of cookies are vanilla, but you think it is equally likely to be any fraction from 0 to 1.\n",
"If you draw three cookies and two are vanilla, what proportion of cookies in the bowl do you think are vanilla?\n",
"The posterior distribution we just computed is the answer to that question.\n",
"\n",
"We'll come back to estimating proportions in the next chapter.\n",
"But first let's use a `Pmf` to solve the dice problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Dice Problem\n",
"\n",
"In the previous chapter we solved the dice problem using a Bayes table.\n",
"Here's the statement of the problem:\n",
"\n",
"> Suppose I have a box with a 6-sided die, an 8-sided die, and a 12-sided die.\n",
"> I choose one of the dice at random, roll it, and report that the outcome is a 1.\n",
"> What is the probability that I chose the 6-sided die?\n",
"\n",
"Let's solve it using a `Pmf`.\n",
"I'll use integers to represent the hypotheses:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.818821Z",
"iopub.status.busy": "2021-04-16T19:35:13.818241Z",
"iopub.status.idle": "2021-04-16T19:35:13.820082Z",
"shell.execute_reply": "2021-04-16T19:35:13.820472Z"
}
},
"outputs": [],
"source": [
"hypos = [6, 8, 12]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can make the prior distribution like this:\n"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.826805Z",
"iopub.status.busy": "2021-04-16T19:35:13.826279Z",
"iopub.status.idle": "2021-04-16T19:35:13.828958Z",
"shell.execute_reply": "2021-04-16T19:35:13.828592Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
"6 0.444444\n",
"8 0.333333\n",
"12 0.222222\n",
"Name: , dtype: float64"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior = prior * likelihood1\n",
"posterior.normalize()\n",
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior probability for the 6-sided die is $4/9$.\n",
"\n",
"Now suppose I roll the same die again and get a 7.\n",
"Here are the likelihoods:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.859245Z",
"iopub.status.busy": "2021-04-16T19:35:13.858594Z",
"iopub.status.idle": "2021-04-16T19:35:13.860853Z",
"shell.execute_reply": "2021-04-16T19:35:13.860402Z"
}
},
"outputs": [],
"source": [
"likelihood2 = 0, 1/8, 1/12"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The likelihood for the 6-sided die is 0 because it is not possible to get a 7 on a 6-sided die.\n",
"The other two likelihoods are the same as in the previous update.\n",
"\n",
"Here's the update:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.867636Z",
"iopub.status.busy": "2021-04-16T19:35:13.867119Z",
"iopub.status.idle": "2021-04-16T19:35:13.870075Z",
"shell.execute_reply": "2021-04-16T19:35:13.869606Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
6
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
8
\n",
"
0.692308
\n",
"
\n",
"
\n",
"
12
\n",
"
0.307692
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"6 0.000000\n",
"8 0.692308\n",
"12 0.307692\n",
"Name: , dtype: float64"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior *= likelihood2\n",
"posterior.normalize()\n",
"posterior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After rolling a 1 and a 7, the posterior probability of the 8-sided die is about 69%."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Updating Dice\n",
"\n",
"The following function is a more general version of the update in the previous section:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.874004Z",
"iopub.status.busy": "2021-04-16T19:35:13.873571Z",
"iopub.status.idle": "2021-04-16T19:35:13.875415Z",
"shell.execute_reply": "2021-04-16T19:35:13.875784Z"
}
},
"outputs": [],
"source": [
"def update_dice(pmf, data):\n",
" \"\"\"Update pmf based on new data.\"\"\"\n",
" hypos = pmf.qs\n",
" likelihood = 1 / hypos\n",
" impossible = (data > hypos)\n",
" likelihood[impossible] = 0\n",
" pmf *= likelihood\n",
" pmf.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first parameter is a `Pmf` that represents the possible dice and their probabilities.\n",
"The second parameter is the outcome of rolling a die.\n",
"\n",
"The first line selects quantities from the `Pmf` which represent the hypotheses.\n",
"Since the hypotheses are integers, we can use them to compute the likelihoods.\n",
"In general, if there are `n` sides on the die, the probability of any possible outcome is `1/n`.\n",
"\n",
"However, we have to check for impossible outcomes!\n",
"If the outcome exceeds the hypothetical number of sides on the die, the probability of that outcome is 0.\n",
"\n",
"`impossible` is a Boolean `Series` that is `True` for each impossible outcome.\n",
"I use it as an index into `likelihood` to set the corresponding probabilities to 0.\n",
"\n",
"Finally, I multiply `pmf` by the likelihoods and normalize.\n",
"\n",
"Here's how we can use this function to compute the updates in the previous section.\n",
"I start with a fresh copy of the prior distribution:\n"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.881390Z",
"iopub.status.busy": "2021-04-16T19:35:13.880863Z",
"iopub.status.idle": "2021-04-16T19:35:13.883304Z",
"shell.execute_reply": "2021-04-16T19:35:13.883671Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
"6 0.000000\n",
"8 0.692308\n",
"12 0.307692\n",
"Name: , dtype: float64"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"update_dice(pmf, 1)\n",
"update_dice(pmf, 7)\n",
"pmf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is the same. We will see a version of this function in the next chapter."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"This chapter introduces the `empiricaldist` module, which provides `Pmf`, which we use to represent a set of hypotheses and their probabilities.\n",
"\n",
"`empiricaldist` is based on Pandas; the `Pmf` class inherits from the Pandas `Series` class and provides additional features specific to probability mass functions.\n",
"We'll use `Pmf` and other classes from `empiricaldist` throughout the book because they simplify the code and make it more readable.\n",
"But we could do the same things directly with Pandas.\n",
"\n",
"We use a `Pmf` to solve the cookie problem and the dice problem, which we saw in the previous chapter.\n",
"With a `Pmf` it is easy to perform sequential updates with multiple pieces of data.\n",
"\n",
"We also solved a more general version of the cookie problem, with 101 bowls rather than two.\n",
"Then we computed the MAP, which is the quantity with the highest posterior probability.\n",
"\n",
"In the next chapter, I'll introduce the Euro problem, and we will use the binomial distribution.\n",
"And, at last, we will make the leap from using Bayes's Theorem to doing Bayesian statistics.\n",
"\n",
"But first you might want to work on the exercises."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:** Suppose I have a box with a 6-sided die, an 8-sided die, and a 12-sided die.\n",
"I choose one of the dice at random, roll it four times, and get 1, 3, 5, and 7.\n",
"What is the probability that I chose the 8-sided die?\n",
"\n",
"You can use the `update_dice` function or do the update yourself."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.904691Z",
"iopub.status.busy": "2021-04-16T19:35:13.904111Z",
"iopub.status.idle": "2021-04-16T19:35:13.907043Z",
"shell.execute_reply": "2021-04-16T19:35:13.906676Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
6
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
8
\n",
"
0.835052
\n",
"
\n",
"
\n",
"
12
\n",
"
0.164948
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"6 0.000000\n",
"8 0.835052\n",
"12 0.164948\n",
"Name: , dtype: float64"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"pmf = prior.copy()\n",
"for data in [1, 3, 5, 7]:\n",
" update_dice(pmf, data)\n",
" \n",
"pmf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:** In the previous version of the dice problem, the prior probabilities are the same because the box contains one of each die.\n",
"But suppose the box contains 1 die that is 4-sided, 2 dice that are 6-sided, 3 dice that are 8-sided, 4 dice that are 12-sided, and 5 dice that are 20-sided.\n",
"I choose a die, roll it, and get a 7.\n",
"What is the probability that I chose an 8-sided die?\n",
"\n",
"Hint: To make the prior distribution, call `Pmf` with two parameters."
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.914938Z",
"iopub.status.busy": "2021-04-16T19:35:13.914323Z",
"iopub.status.idle": "2021-04-16T19:35:13.916995Z",
"shell.execute_reply": "2021-04-16T19:35:13.917367Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
4
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
6
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
8
\n",
"
0.391304
\n",
"
\n",
"
\n",
"
12
\n",
"
0.347826
\n",
"
\n",
"
\n",
"
20
\n",
"
0.260870
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"4 0.000000\n",
"6 0.000000\n",
"8 0.391304\n",
"12 0.347826\n",
"20 0.260870\n",
"Name: , dtype: float64"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"# Notice that I don't bother to normalize the prior.\n",
"# The `Pmf` gets normalized during the update, so we\n",
"# don't have to normalize it before.\n",
"\n",
"ps = [1,2,3,4,5]\n",
"qs = [4,6,8,12,20]\n",
"pmf = Pmf(ps, qs)\n",
"update_dice(pmf, 7)\n",
"pmf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:** Suppose I have two sock drawers.\n",
"One contains equal numbers of black and white socks.\n",
"The other contains equal numbers of red, green, and blue socks.\n",
"Suppose I choose a drawer at random, choose two socks at random, and I tell you that I got a matching pair.\n",
"What is the probability that the socks are white?\n",
"\n",
"For simplicity, let's assume that there are so many socks in both drawers that removing one sock makes a negligible change to the proportions."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.925514Z",
"iopub.status.busy": "2021-04-16T19:35:13.924848Z",
"iopub.status.idle": "2021-04-16T19:35:13.927865Z",
"shell.execute_reply": "2021-04-16T19:35:13.927501Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
BlackWhite
\n",
"
0.6
\n",
"
\n",
"
\n",
"
RedGreenBlue
\n",
"
0.4
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"BlackWhite 0.6\n",
"RedGreenBlue 0.4\n",
"Name: , dtype: float64"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"# In the BlackWhite drawer, the probability of getting a match is 1/2\n",
"# In the RedGreenBlue drawer, the probability of a match is 1/3\n",
"\n",
"hypos = ['BlackWhite', 'RedGreenBlue']\n",
"prior = Pmf(1/2, hypos)\n",
"likelihood = 1/2, 1/3\n",
"posterior = prior * likelihood\n",
"posterior.normalize()\n",
"posterior"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.931373Z",
"iopub.status.busy": "2021-04-16T19:35:13.930918Z",
"iopub.status.idle": "2021-04-16T19:35:13.933340Z",
"shell.execute_reply": "2021-04-16T19:35:13.933685Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.30000000000000004"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"# If I drew from the BlackWhite drawer, the probability the\n",
"# socks are white is 1/2\n",
"\n",
"posterior['BlackWhite'] / 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:** Here's a problem from [Bayesian Data Analysis](http://www.stat.columbia.edu/~gelman/book/):\n",
"\n",
"> Elvis Presley had a twin brother (who died at birth). What is the probability that Elvis was an identical twin?\n",
"\n",
"Hint: In 1935, about 2/3 of twins were fraternal and 1/3 were identical."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.942870Z",
"iopub.status.busy": "2021-04-16T19:35:13.942005Z",
"iopub.status.idle": "2021-04-16T19:35:13.950507Z",
"shell.execute_reply": "2021-04-16T19:35:13.950906Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
prior
\n",
"
likelihood
\n",
"
unnorm
\n",
"
posterior
\n",
"
\n",
" \n",
" \n",
"
\n",
"
identical
\n",
"
0.333333
\n",
"
1.0
\n",
"
0.333333
\n",
"
0.5
\n",
"
\n",
"
\n",
"
fraternal
\n",
"
0.666667
\n",
"
0.5
\n",
"
0.333333
\n",
"
0.5
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" prior likelihood unnorm posterior\n",
"identical 0.333333 1.0 0.333333 0.5\n",
"fraternal 0.666667 0.5 0.333333 0.5"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Solution\n",
"\n",
"# The trick to this question is to notice that Elvis's twin was a brother.\n",
"# If they were identical twins, it is certain they would be the same sex.\n",
"# If they were fraternal twins, the likelihood is only 50%.\n",
"\n",
"# Here's a solution using a Bayes table\n",
"\n",
"import pandas as pd\n",
"\n",
"table = pd.DataFrame(index=['identical', 'fraternal'])\n",
"table['prior'] = 1/3, 2/3\n",
"table['likelihood'] = 1, 1/2\n",
"\n",
"table['unnorm'] = table['prior'] * table['likelihood']\n",
"prob_data = table['unnorm'].sum()\n",
"\n",
"table['posterior'] = table['unnorm'] / prob_data\n",
"table"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-16T19:35:13.957464Z",
"iopub.status.busy": "2021-04-16T19:35:13.956748Z",
"iopub.status.idle": "2021-04-16T19:35:13.959501Z",
"shell.execute_reply": "2021-04-16T19:35:13.959878Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"