{ "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": { "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": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "# Get utils.py\n", "\n", "import os\n", "\n", "if not os.path.exists('utils.py'):\n", " !wget https://github.com/AllenDowney/ThinkBayes2/raw/master/code/soln/utils.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Distributions\n", "\n", "In statistics a **distribution** is a set of values and their\n", "corresponding probabilities.\n", "\n", "For example, if you toss a coin, there are two possible outcomes with\n", "approximately equal probabilities.\n", "\n", "If you roll a six-sided die, the set of possible values is the numbers 1\n", "to 6, and the probability associated with each value 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", "\n", "This library provides a class called `Pmf`, which represents a\n", "**probability mass function**.\n", "\n", "`empiricaldist` is available from the Python Package Index (PyPI). You\n", "can [download it here](https://pypi.org/project/empiricaldist/) or\n", "install it with `pip`. For more details, see the Preface.\n", "\n", "To use `Pmf` you can import it like this:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from empiricaldist import Pmf" ] }, { "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": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
heads0.5
tails0.5
\n", "
" ], "text/plain": [ "heads 0.5\n", "tails 0.5\n", "dtype: float64" ] }, "execution_count": 4, "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": [ "The two outcomes have the same probability, $0.5$.\n", "\n", "This example makes a `Pmf` that represents the distribution of outcomes\n", "of a six-sided die:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
11
21
31
41
51
61
\n", "
" ], "text/plain": [ "1 1\n", "2 1\n", "3 1\n", "4 1\n", "5 1\n", "6 1\n", "dtype: int64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die = Pmf()\n", "for x in [1,2,3,4,5,6]:\n", " die[x] = 1\n", " \n", "die" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Pmf` creates an empty `Pmf` with no values. The `for` loop adds the\n", "values $1$ through $6$, each with \"probability\" $1$.\n", "\n", "In this `Pmf`, the probabilities don't add up to 1, so they are not\n", "really probabilities. We can use `normalize` to make them add up to 1." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die.normalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The return value from `normalize` is the sum of the probabilities before normalizing.\n", "\n", "Now we can see that the total is 1 (at least within floating-point error)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
10.166667
20.166667
30.166667
40.166667
50.166667
60.166667
\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", "dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.9999999999999999" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way make a `Pmf` is to provide a sequence of values.\n", "\n", "In this example, every value appears once, so they all have the same probability." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
10.166667
20.166667
30.166667
40.166667
50.166667
60.166667
\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", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die = Pmf.from_seq([1,2,3,4,5,6])\n", "die" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More generally, values can appear more than once, as in this example:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
M0.090909
i0.363636
p0.181818
s0.363636
\n", "
" ], "text/plain": [ "M 0.090909\n", "i 0.363636\n", "p 0.181818\n", "s 0.363636\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "letters = Pmf.from_seq(list('Mississippi'))\n", "letters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "The `Pmf` class inherits from a Pandas `Series`, so anything you can do\n", "with a `Series`, you can also do with a `Pmf`.\n", "\n", "For example, you can use the bracket operator to look up a value and\n", "returns the corresponding probability." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.36363636363636365" ] }, "execution_count": 11, "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 value that's not in the\n", "distribution, you get a `KeyError`.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "KeyError\n" ] } ], "source": [ "try:\n", " letters['t']\n", "except KeyError as e:\n", " print('KeyError')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also call a `Pmf` as if it were a function, with a value in parentheses." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.36363636363636365" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "letters('s')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the value is in the distribution the results are the same. \n", "But if the value is not in the distribution, the result is $0$, not an error." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "letters('t')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With parentheses, you can also provide a sequence of values and get a sequence of probabilities." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.16666667, 0.16666667, 0. ])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "die([1,4,7])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As these examples shows, the values in a `Pmf` can be integers or strings. In general, they can be any type that can be stored in the index of a Pandas `Series`.\n", "\n", "If you are familiar with Pandas, that will help you work with `Pmf` objects. But I will explain what you need to know as we go along." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The cookie problem\n", "\n", "In this section I'll use a `Pmf` to solve the cookie problem from Section XX.\n", "Here's the statement of the problem again:\n", "\n", "> Suppose there are two bowls of cookies.\n", "> Bowl 1 contains 30 vanilla cookies and 10 chocolate cookies.\n", "> Bowl 2 contains 20 of each.\n", ">\n", "> Now suppose you choose one of the bowls at random and, without\n", "> looking, select a cookie at random. The cookie is vanilla. What is\n", "> 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": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
Bowl 10.5
Bowl 20.5
\n", "
" ], "text/plain": [ "Bowl 1 0.5\n", "Bowl 2 0.5\n", "dtype: float64" ] }, "execution_count": 16, "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": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
Bowl 10.375
Bowl 20.250
\n", "
" ], "text/plain": [ "Bowl 1 0.375\n", "Bowl 2 0.250\n", "dtype: float64" ] }, "execution_count": 17, "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.\n", "We can use `normalize` to compute the posterior probabilities:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.625" ] }, "execution_count": 18, "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": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
Bowl 10.6
Bowl 20.4
\n", "
" ], "text/plain": [ "Bowl 1 0.6\n", "Bowl 2 0.4\n", "dtype: float64" ] }, "execution_count": 19, "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": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.6" ] }, "execution_count": 20, "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": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
Bowl 10.692308
Bowl 20.307692
\n", "
" ], "text/plain": [ "Bowl 1 0.692308\n", "Bowl 2 0.307692\n", "dtype: float64" ] }, "execution_count": 21, "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", "Here's the update." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
Bowl 10.529412
Bowl 20.470588
\n", "
" ], "text/plain": [ "Bowl 1 0.529412\n", "Bowl 2 0.470588\n", "dtype: float64" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "likelihood_chocolate = [0.25, 0.5]\n", "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 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 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 represent 101 hypotheses, numbered from 0 to 100." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "hypos = np.arange(101)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is a NumPy array, which we can use to make the prior distribution:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "101" ] }, "execution_count": 24, "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 values.\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", "\n", "Since all hypotheses have the same prior probability, this distribution is **uniform**.\n", "\n", "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": 25, "metadata": {}, "outputs": [], "source": [ "likelihood_vanilla = hypos/100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can compute the posterior distribution in the usual way:\n" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5000000000000001" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "posterior1 = prior * likelihood_vanilla\n", "posterior1.normalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following figure shows the prior distribution and the posterior distribution after one vanilla cookie." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "from utils import decorate, savefig\n", "\n", "def decorate_bowls(title):\n", " decorate(xlabel='Bowl #',\n", " ylabel='PMF',\n", " title=title)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1fnH8c9DWJWtIIoSNgVENhUj4FJ3ERVBW6xIK+BSaiu1ttpfcVfUqrVqtVgVRUSroFKBaFVUcK2IBLES9gBBwr7vW5Ln98e9seOQjZDJTDLf9+s1r9zl3HPPvTPJk3vvmeeYuyMiIpJoqsW7ASIiIoVRgBIRkYSkACUiIglJAUpERBKSApSIiCQkBSgREUlIClBSpZjZdjM7uoL2dZqZLQr3eWlF7LOyMbOfm9n7EfNuZm3C6RfN7P74te6HzCzbzM4rYt0cMzurgpuU9BSgpEzCX+Zd4R/nNWY22szqHkR9rcI/XtUPpl3uXtfdlxxMHQdgODAi3OdEM/vYzK6roH1XCu7+irv3jHc7Dpa7d3T3j+PdjmSjACUH4xJ3rwt0BU4G7ohXQw42sJVx+5bAnIPZb1QbUsqrLpGqQAFKDpq7rwDeBToBmNlRZpZuZhvNLMvMfllQ1sy6mVmGmW0Nr7weC1d9Gv7cHF6VnRKWv8bM5pnZJjObbGYtI+pyM7vBzBYBiyKWFdxCamBmL5nZOjNbZmZ3mFm1cN1gM/uPmT1uZhuBe6KPK2zrNDPbbGarzGyEmdUM1y0GjgbeCtv7IPBjYEQ4PyIs197MPgjPxQIz+1lE/S+a2dNm9o6Z7QDOLqQNxZ3Le8zs9fAYt4W3odKitv1XePxLzezGwt4/M+thZqsjA6SZXWZm35Z0HiLO+fXh7c5NZvaUmVnEef68sP1GteFHZvZ22NZN4XRqMeWbm9mbYfkNEee7Wvg+LzOzteG5aRCxXZ/wPG0Or3iPK6L+9uE56x/Of3/7L9zHMDNbHO77dTNrVNIxShm4u156HfALyAbOC6ebE1xJ3BfOfwL8A6gNnACsA84N100Drgqn6wI9wulWgAPVI/ZxKZAFHAdUJ7hC+yJivQMfAI2AOhHL2oTTLwGTgHph/QuBa8N1g4Fc4Ldh3XUKOcaTgB7h+lbAPOCmws5BOP8xcF3E/KHAcuDqsI6uwHqgY7j+RWALcBrBP4u1C2lDcefyHmA3cBGQAjwIfBmuqwbMBO4CahIE0yXABUW8n4uB8yPm3wCGlfI8OPA20BBoEbaxV8R5/jyqbJuI478/nG4M/BQ4JHy/3gAmFtHWFOC/wOPhOa4NnB6uu4bgM3M0wefrTeDlcF07YAdwPlAD+L+wbM3I9zN8n74Dehfxeb8J+BJIBWoBzwJj4/07WRVfcW+AXpXzFf7Cbgc2A8vCP6J1CIJVHlAvouyDwIvh9KfAvcBhUfW1Yv8A9S5hQAnnqwE7gZbhvAPnRNXjQJvwj9geoEPEul8BH4fTg4HvDvCYbwImRJ2D4gLUFcBnUXU8C9wdTr8IvFTM/ko6l/cAH0as6wDsCqe7Rx8fcCswuoh93Q+8EE7XC/+QtyzleXDCABHOv87/gttgShGgCtnHCcCmItadQhAEqxeybgrwm4j5Y4F9BMH1TuD1qM/TCuCsiPfzXiAHOLuQz3tBgJpH+E9COH9kwT4q8ncwGV66xScH41J3b+juLd39N+6+CzgK2Oju2yLKLQOahdPXEvwnO9/MZphZ72Lqbwk8Ed6O2QxsBCyiLgiuUApzGMGVw7Ii2lHctgCYWbvwVtNqM9sK/Dmst7RaAt0L2h8ew8+BpqVsQ0nnEmB1xPROoLYFz9NaAkdF7fs24Igi9vUq8BMzqwX8BPja3ZdBqc9DdDsOqMOMmR1iZs+Gt+a2Evwj09AKfy7XHFjm7rmFrDuK/d/z6gTH/YN17p5PcP4jz+f1BFfpHxXT3JbAhIjzOo/gH4mizq2UkQKUlLeVQCMzqxexrAXBf6q4+yJ3vxI4HHgYGG9mhxL8Zx1tOfCrMAgWvOq4+xcRZYpKx7+e4L/alhHLvm9HCdsWeBqYD7R19/oEf+CtmPLR9S0HPolqf113/3Up21DsuSzBcmBp1L7ruftFhTbcfS7BH+8LgQEEAavAgZ6HsriZ4Gqne7iPM8Llhe1nOdDCCu/YspL93/NcYE30uvA5WXN+eD6vD+t+vJi2LgcujDq3tT14FivlSAFKypW7Lwe+AB40s9pm1oXgqukVADP7hZk1Cf973Rxulkdwyyaf4NlBgWeAW82sY7htAzO7vJTtyCO41fSAmdWzoHPFH4B/HsDh1AO2AtvNrD3w6xLKr4lq/9tAOzO7ysxqhK+Ti3owX8gxFHsuS/AVsNXM/mRmdcwsxcw6mdnJxWzzKnAjQXB4I2L5gZ6HsqgH7CLoJNMIuLuYsl8Bq4CHzOzQ8NycFq4bC/zezFpb8LWHPwOvhVdbrwMXm9m5ZlaDICjuITjHBbYBvYAzzOyhIvb/DMHnqiWAmTUxs75lOWgpngKUxMKVBM+UVgITCJ65fBCu6wXMMbPtwBNAf3ff7e47gQeA/4S3Tnq4+wSCq6xx4W2fTIL/8EvrtwTPUpYAnxP8AX7hALa/heBqYhvwHPBaCeWfAPqFvdCeDG/N9QT6E5yL1QTHU+sA2lDcuSxSGKAvIXiWs5TgivJ5oEExm40FzgKmuvv6iOUHeh7K4m8EzzDXE3RAeK+oghHH1oagM0MOwfM+CN7flwluES4l6ETy23C7BcAvgL+H+7mE4KsSe6Pq30zQkeJCM7uvkCY8AaQD75vZtrC93Q/4iKVE5q4BC0VEJPHoCkpERBKSApSIiCQkBSgREUlIClAiIpKQDirBZmV32GGHeatWreLdDBGRpDZz5sz17t4kenlSB6hWrVqRkZER72aIiCQ1M1tW2HLd4hMRkYSkACUiIglJAUpERBJSTJ9BmVkvgrQgKcDz7v5Q1PpaBGP2nARsAK5w92wzOx94iCAb9V7gj+4+NdzmJII0/XWAd4DfubuH+bteI0gLkw38zN03HWib9+3bR05ODrt37z7wAxZq165NamoqNWrUiHdTRKSSi1mACtPkP0WQ0yoHmGFm6WHW5ALXEoz50iYcufJhgpxa6wlyZK00s07AZP6XEv9pYAhB/qt3CHK7vQsMA6a4+0NmNiyc/9OBtjsnJ4d69erRqlUrwkFBpZTcnQ0bNpCTk0Pr1q3j3RwRqeRieYuvG5Dl7kvCZIzjgOiMv32BMeH0eOBcMzN3n+XuK8PlcwjGuKllZkcC9d19mgdJBF8iGHU1uq4xEcsPyO7du2ncuLGCUxmYGY0bN9bVp4iUi1gGqGb8cDC2HH44MNgPyoTp8LcQDP0c6afALHffE5bPKaLOI9x9VVjXKoLxhvZjZkPMLMPMMtatW1dowxWcyk7nTkTKSywDVGF/qaJTpxdbJhwH6GGCobpLW2ex3H2ku6e5e1qTJvt9L0xEREpp1948Rn2+lNy8/JjUH8sAlUMwWmWBVIIxbQotE46O2YBgWG/MLJVg/JuB7r44onxqEXWuCW8BEv5cW25HkqDuuusuPvzww3g3Q0SS0EcL1tLzb59w39tz+SxrfckblEEsA9QMoG04smVNgkHb0qPKpAODwul+BAOluZk1BP4N3Oru/ykoHN6622ZmPcLhmgcCkwqpa1DE8iopLy+P4cOHc9555x3QNiIiB2PN1t385pWZXD16BjVTqjH2lz04+9hCn6gctJgFqPCZ0lCCHnjzgNfdfY6ZDTezPmGxUUBjM8siGI57WLh8KMFomXea2Tfhq+AM/JpgZNAsYDFBDz4IuqWfb2aLCHoOFjVcc8LLzs6mffv2DBo0iC5dutCvXz927txJq1atGD58OKeffjpvvPEGgwcPZvz48QBMmTKFE088kc6dO3PNNdewZ88egP22EREpi7x858X/LOXcRz9hyry13NKzHe/+7gxOOSa620D5ien3oNz9HYKu4JHL7oqY3g1cXsh29wP3F1FnBtCpkOUbgHMPssk/8N5777F69eryrJKmTZvSq1evEsstWLCAUaNGcdppp3HNNdfwj3/8Awi+Z/T5559/3z4Ieh4OHjyYKVOm0K5dOwYOHMjTTz/NTTfdtN82IiIH6tuczdw+IZPZK7ZwRrsm3Ne3Iy0bHxrz/SqTRIJq3rw5p512GgC/+MUvvg8wV1xxxX5lFyxYQOvWrWnXrh0AgwYN4tNPP/1+fWHbiIiUZOvufdw9KZO+T/2H1Vt38/crT2TM1SdXSHCCJM9mXpLSXOnESnR37YL5Qw/d/4MRfCWsaIVtIyJSFHfnndmrufetOazbvoeBPVpy8wXHUr92xWaI0RVUgvruu++YNm0aAGPHjuX0008vsmz79u3Jzs4mKysLgJdffpkzzzyzQtopIlXLdxt2Mnj0DG549WsOr1+LSTecxr19O1V4cAIFqIR13HHHMWbMGLp06cLGjRv59a9/XWTZ2rVrM3r0aC6//HI6d+5MtWrVuP766yuwtSJS2e3Nzeepj7I4//FPmLlsE3df0oFJN5xOl9SGcWuTbvElqGrVqvHMM8/8YFl2dvYP5l988cXvp88991xmzZq1Xz3R24iIRJu+ZAO3T8wka+12LuzUlLsv6UjTBrXj3SwFKBGRZLVxx17+/M48xs/MIfVHdXhhcBrntD8i3s36ngJUAmrVqhWZmZnxboaIVFH5+c74mTn8+d15bN+dy6/POoYbz2lLnZop8W7aDyhAiYgkkYVrtnH7hNnMyN7Eya1+xAOXdabdEfXi3axCKUCJiCSBXXvzeHLqIp77dAl1a1fnLz/tQr+TUqlWLXFHIFCAEhGp4j5asJa7JmWyfOMufto1ldsuak/jurXi3awSKUCJiFRRq7fsZvjbc3hn9mraHF6XcUN60OPo2OXOK28KUFXUxIkTadeuHR06dDig7dLT05k7dy7Dhg0rubCIJKS8fOeladk8+v5C9uXlc0vPdgw54xhqVq9cX31VgKqiJk6cSO/evQ8oQOXm5tKnTx/69OlTcuGIbapX18dIJFF8m7OZ2ybMJnPF1gpN7BoL+suSgLKzs+nVqxfdu3dn1qxZtGvXjpdeeolp06Zxyy23kJuby8knn8zTTz9NrVq1GDZsGOnp6VSvXp2ePXvyk5/8hPT0dD755BPuv/9+/vWvfwFwww03sG7dOg455BCee+452rdvz+DBg2nUqBGzZs2ia9eudO7cmYyMDEaMGMGyZcu45pprWLduHU2aNGH06NG0aNFiv20effTROJ8xEdm6ex+PTl7AS18uo0ndWowYcCIXdz5yv7yelYkCVDHufWsOc1duLdc6OxxVn7sv6VhiuejhNh577DGeffbZ/YbUGDhwIBMmTGD+/PmYGZs3b6Zhw4b06dOH3r17069fPyDINPHMM8/Qtm1bpk+fzm9+8xumTp0KwMKFC/nwww9JSUn5QXaKoUOHMnDgQAYNGsQLL7zAjTfeyMSJE/fbRkTix915+9tVDH97LuvjmNg1FirXDckkEj3cxpQpUwodUqN+/frUrl2b6667jjfffJNDDjlkv7q2b9/OF198weWXX84JJ5zAr371K1atWvX9+ssvv7zQQDNt2jQGDBgAwFVXXfWDMaWK2kZEKs6yDTsYNHoGvx07iyPinNg1FmJ6BWVmvYAngBTgeXd/KGp9LeAl4CRgA3CFu2ebWWNgPHAy8KK7Dw3L1wM+i6giFfinu99kZoOBR4AV4boR7v78wbS/NFc6sVLay/Lq1avz1VdfMWXKFMaNG8eIESO+vzIqkJ+fT8OGDfnmm28KraO0w3FEtklDeIjEz57cPEZ+soQRH2VRI6Uad1/SgYGntCIlgb/TVBYxu4IysxTgKeBCoANwpZlFP7G/Ftjk7m2Ax4GHw+W7gTuBWyILu/s2dz+h4AUsA96MKPJaxPqDCk7xFj3cxnnnnVfokBrbt29ny5YtXHTRRfztb3/7PgjVq1ePbdu2AVC/fn1at279/ZDv7s5///vfEttw6qmnMm7cOABeeeWVYof8EJGK8eWSDVz0xGc8+sFCzjvuCD78w5lcfVrrKhecILa3+LoBWe6+xN33AuOAvlFl+gJjwunxwLlmZu6+w90/JwhUhTKztsDh/PCKqsqIHm7j97//faFDamzbto3evXvTpUsXzjzzTB5//HEA+vfvzyOPPMKJJ57I4sWLeeWVVxg1ahTHH388HTt2ZNKkSSW24cknn2T06NF06dKFl19+mSeeeCLWhy0iRdiwfQ83v/5f+o/8kr15+YwefDJP/bxrQmQdjxUraTTWMlds1g/o5e7XhfNXAd0LbteFyzLDMjnh/OKwzPpwfjCQFrlNxLZ3AfXd/ZaIsg8C64CFwO/dfXlxbUxLS/OMjIwfLJs3bx7HHXdcmY65vGRnZ9O7d+9KmzA2Ec6hSFWRn++8MXM5D747nx17chlyxtEMPTvxErseDDOb6e5p0ctj+QyqsOvN6GhYmjJF6Q9cFTH/FjDW3feY2fUEV2bn7NcosyHAEIAWLVqUclciIhVvweogsWvGsk10a9WIBy7rRNsETewaC7EMUDlA84j5VGBlEWVyzKw60ADYWFLFZnY8UN3dZxYsc/cNEUWe43/Ps37A3UcCIyG4gir5MCqehtsQSW479+by5JQsnv9sCfVqV+cv/bpw+Umplfo7TWURywA1A2hrZq0Jetb1BwZElUkHBgHTgH7AVC/dPccrgbGRC8zsSHcv6DvdB5hX1oa7e9J9EMpLrG4ZiySLqfPXcOfEOazYvIufpaUy7MLjaHRozXg3Ky5iFqDcPdfMhgKTCbqZv+Duc8xsOJDh7unAKOBlM8siuHLqX7C9mWUD9YGaZnYp0NPd54arfwZcFLXLG82sD5Ab1jW4LO2uXbs2GzZsoHHjxgpSB8jd2bBhA7VrV92HtiKxsmrLLu5Nn8t7c4LErq8N6UH3SpTYNRZi1kmiMiisk8S+ffvIyclh9+4iOxBKMWrXrk1qaio1alSNLwqKxFpuXj5jpi3jsfcXkJvv3HhuW37546MrXWLXgxGPThKVUo0aNWjdunW8myEiSeC/y4PErnNWbuXMdk24r28nWjTePxtMslKAEhGpYFt37+OvkxfwcpjY9akBXbmoc1M9VoiiACUiUkEiE7tu2L6HQae04uae7ahXRXLnlTcFKBGRCpC9fgd3Tsrks0Xr6dysAaMGpdEltWG8m5XQFKBERGJoT24ez4aJXWumVOOeSzpwVRVM7BoLClAiIjEybfEG7pg4m8XrdnBxlyO5q3cHjqivr2GUlgKUiEg527B9D39+Zz7/+jqH5o3q8OLVJ3PWsYfHu1mVjgKUiEg5yc93Xs8IErvu3JvL0LPbMPScNtSuUXUSu1YkBSgRkXIwf/VWbp+Qycxlm+jWuhEPXJpciV1jQQFKROQg7NybyxNTFjHqs6XUq12dR/p1oV8SJnaNBQUoEZEymjJvDXdNChK7Xn5SKrdelLyJXWNBAUpE5ABFJnZtq8SuMaMAJSJSSpGJXfPc+b9ex3Ld6cmV2LUiKUCJiJTCN8s3c3uY2PWsY4PErs0bKbFrLClAiYgUY8uuILHrP6cv4/B6tXj6513p1UmJXSuCApSISCHcnbe+XcV9YWLXwae24g/nK7FrRVKAEhGJEpnYtUtqA14YdDKdUxvEu1lJJ6ZP9sysl5ktMLMsMxtWyPpaZvZauH66mbUKlzc2s4/MbLuZjYja5uOwzm/C1+HF1SUiUlp7cvN4csoiev7tU2Z9t5l7+3Rkwm9OU3CKk5hdQZlZCvAUcD6QA8wws3R3nxtR7Fpgk7u3MbP+wMPAFcBu4E6gU/iK9nN3z4haVlRdIiIl+mLxeu6YmMkSJXZNGLG8guoGZLn7EnffC4wD+kaV6QuMCafHA+eambn7Dnf/nCBQlVahdZW9+SKSDNZv38MfXvuGAc9NJzfPefHqk3lqQFcFpwQQy2dQzYDlEfM5QPeiyrh7rpltARoD60uoe7SZ5QH/Au53dy9tXWY2BBgC0KJFizIclohUBfn5zmsZy3koTOz623PacMPZSuyaSGIZoAq7evEylIn2c3dfYWb1CALUVcBLpa3L3UcCIwHS0tJK2peIVEGRiV27t27EA5d1os3hSuyaaGIZoHKA5hHzqcDKIsrkmFl1oAGwsbhK3X1F+HObmb1KcCvxpbLUJSLJZefeXJ74cBHPf76UBnVqKLFrgotlgJoBtDWz1sAKoD8wIKpMOjAImAb0A6aGt+sKFQaehu6+3sxqAL2BD8tSl4gklw/nruHu9CCx6xVpzRl2YXt+pMSuCS1mASp8DjQUmAykAC+4+xwzGw5kuHs6MAp42cyyCK52+hdsb2bZQH2gppldCvQElgGTw+CUQhCcngs3KbIuEUleKzfv4p70Obw/dw1tD6/L6786hW6tG8W7WVIKlswXGWlpaZ6REd1bXUSqgty8fF78IpvHPlhIvju/O7cd157eWoldE5CZzXT3tOjlyiQhIlXOrO82cduETOat2so57Q/n3j4dldi1ElKAEpEqY8uufTwyeT6vTP+OI+rV5plfdOWCjkrsWlkpQIlIpefupP93Jfe9PY+NO/Zw9amt+UPPdtStpT9xlZnePRGp1CITux6f2oAXrz6ZTs2UO68qUIASkUppT24ez3y8hKc+zqJWSjXu69uRAd1bklJNt/OqCgUoEal0vsgKE7uu38Elxx/FnRcfx+HKnVflKECJSKWxfvseHvj3PCbMWkHLxocw5ppunNmuSbybJTGiACUiCS8/3xk3YzkPvTuPXfvylNg1SShAiUhCm7dqK3dMDBK79ji6Efdf2pk2h9eNd7OkAihAiUhCikzs2rBODR772fFcdmIzfacpiShAiUjC+WDuGu4JE7te2a05f+rVnoaHKLFrslGAEpGEEZnY9dgj6jH++lNIa6XErslKAUpE4i46seuferXnuh+3pkaKErsmMwUoEYmrr7/bxO1hYtdz2x/OPUrsKiEFKBGJiy079/Hw5PmM/aogsetJXNDxCHWCkO8pQIlIhXJ3Jn2zkvv/PZeNO/ZyzWmt+f35Suwq+4vpDV4z62VmC8wsy8yGFbK+lpm9Fq6fbmatwuWNzewjM9tuZiMiyh9iZv82s/lmNsfMHopYN9jM1pnZN+Hrulgem4gcuCXrtnPVqK+46bVvaNawDulDT+fO3h0UnKRQMftUmFkK8BRwPpADzDCzdHefG1HsWmCTu7cxs/7Aw8AVwG7gTqBT+Ir0V3f/yMxqAlPM7EJ3fzdc95q7D43VMYlI2ezel8cznyzmHx8tplZ1JXaV0onlvy3dgCx3XwJgZuOAvkBkgOoL3BNOjwdGmJm5+w7gczNrE1mhu+8EPgqn95rZ10BqDI9BRA7S54vWc+ekTJau30Gf44/ijt7HcXg9JXaVksUyQDUDlkfM5wDdiyrj7rlmtgVoDKwvqXIzawhcAjwRsfinZnYGsBD4vbsvL2S7IcAQgBYtWpT6YETkwKzbtocH/j2Xid+spFXjQ3j52m78uK0Su0rpxTJAFXbt7mUos3/FZtWBscCTBVdowFvAWHffY2bXA2OAc/ar3H0kMBIgLS2txH2JyIHJz3fGzviOh9+dz659edx4Tht+o8SuUgaxDFA5QPOI+VRgZRFlcsKg0wDYWIq6RwKL3P1vBQvcfUPE+ucInmeJSAWau3Irt0+czazvNiuxqxy0WAaoGUBbM2sNrAD6AwOiyqQDg4BpQD9gqrsXe1VjZvcTBLLropYf6e6rwtk+wLyDPgIRKZUde3L524cLeeE/2UrsKuUmZgEqfKY0FJgMpAAvuPscMxsOZLh7OjAKeNnMsgiunPoXbG9m2UB9oKaZXQr0BLYCtwPzga/DD/8Id38euNHM+gC5YV2DY3VsIvI/789ZzT3pc1i5ZbcSu0q5shIuWKq0tLQ0z8jIiHczRCqlFZt3cfekOXw4bw3tm9bjgcs6cVJLJXaVA2dmM909LXq5vh0nIgdkX14+o/+zlMc/WATArRe255rTldhVyp8ClIiU2tffbeK2N2czf/U2zjsuSOya+iMldpXYUIASkRJFJnZtWr82z151Ej07KLGrxJYClIgUyd2Z+M0KHvj3PDbt3Me1YWLXQ5U7TyqAPmUiUqjF67Zz58RMvli8gROaN2TMNZ3oeFSDeDdLkogClIj8wO59efzj48U88/FiatWoxv2XdmJAtxZUU2JXqWAKUCLyvcjErn1POIo7Lu5Ak3q14t0sSVIKUCLCum17uP/fc5kUJnb957XdOb3tYfFuliQ5BSiRJJaf77z61Xc8/N589uzL58Zz2/Kbs45RYldJCApQIklq7sqt3DZhNt8s38ypxzTmvks7cUwTJXaVxKEAJZJkduzJ5fEPFjL6iyCx6+NXHM+lJyixqyQeBSiRJDI5TOy6astuBnRvwZ8uaE+DQ2rEu1kihVKAEkkCOZt2ck/63O8Tu44Y0JWTWv4o3s0SKVaxAcrMXnT3weH0IHcfUyGtEpFyEZ3Y9baL2nP1aUrsKpVDSVdQx0dM/45gGHURqQRmLtvI7RMyldhVKq2SAlTyDhYlUklt3rmXh99bwNivvuOoBrUZedVJ9OzYNN7NEjlgJQWoVDN7ErCI6e+5+43FbWxmvYAnCEbUfd7dH4paXwt4CTgJ2ABc4e7ZZtYYGA+cDLzo7kMjtjkJeBGoA7wD/M7d3cwaAa8BrYBs4GfuvqmE4xOpMtydCbOCxK6bd+3jutOV2FUqt5I+uX+MmD6goWfNLAV4CjgfyAFmmFm6u8+NKHYtsMnd25hZf+Bh4ApgN3An0Cl8RXoaGAJ8SRCgegHvAsOAKe7+kJkNC+f/dCBtFqmsFq/bzh0TMpm2JEjs+tJlSuwqlV+xAeogO0V0A7LcfQmAmY0D+gKRAaovcE84PR4YYWbm7juAz82sTWSFZnYkUN/dp4XzLwGXEgSovsBZYdExwMcoQEkVp8SuUpWV1Isvvbj17t6nmNXNgOUR8zlA96LKuHuumW0BGgPri6kzJ6rOZuH0Ee6+KqxrlZkdXlgFZjaE4AqMFi1aFNP84r333nusXr26zNuLHKzFOyciRq0AABlFSURBVGrw9up6bNqXQuf6u7ng8O3kzl/FS/Pj3TJJJk2bNqVXr14xqbukW3ynEASQscB0gmdRpVVY2ehOF6UpczDl9y/sPhIYCZCWlqZOIFLpbMs1Jq+tS+bW2jSqkctVzTdzzKH74t0skXJXUoBqSvAM6UpgAPBvYKy7zylF3TlA84j5VGBlEWVyzKw60ADYWEKdqUXUucbMjgyvno4E1paijWUWq/8YRIqSFyZ2HRkmdr3pvGO4/kwldpWqq9hv67l7nru/5+6DgB5AFvCxmf22FHXPANqaWWszqwn0B6JvGaYDg8LpfsBUdy/yqia8hbfNzHpYkDhsIDCpkLoGRSwXqfTmrNzCT57+gjsnZtIltQHv3fRjbjqvnYKTVGkl9j8Nu4JfTHAV1Qp4EnizpO3CZ0pDgckE3cxfcPc5ZjYcyHD3dGAU8LKZZRFcOfWP2G82UB+oaWaXAj3DHoC/5n/dzN8NXwAPAa+b2bXAd8DlJbVRJNFtL0js+p+lNDq0Jk/0P4E+xx+lxK6SFKyYCxbMbAxBN+93gXHunllRDasIaWlpnpFxQL3nRSqEuzN5zhrufWsOq7fu5spuSuwqVZeZzXT3tOjlJV1BXQXsANoBvzOzgmhmgLt7/fJtpogEiV3n8OG8tbRvWo+nft6Vri2U2FWST0nfg1JGSZEKsi8vn1GfL+WJD5XYVQRK/h5UbeB6oA3wLcFzpNyKaJhIMolM7Hp+hyO4p09HmjWsE+9micRVSbf4xgD7gM+Ai4COBFnNRaQcbN65l4fenc+4GcuV2FUkSkkBqoO7dwYws1HAV7FvkkjV5+68+fUKHnhnHlt27WPIGUfzu3PbKrGrSISSfhu+/3p62G08xs0Rqfqy1m7njomz+XLJRk5s0ZAHLu1Mh6PU30gkWokDFprZ1nDagDrhvHrxiRyg3fvy+MdHWTz9yWLq1Ejhz5d1pv/JzZXYVaQIJfXi09fURcrBpwvXceekTJZt2MmlJxzF7Rd3oEm9WvFulkhC0w1vkRhau2039709j7f+u5KjDzuUV6/rzqltDot3s0QqBQUokRjIy3denb6Mv0xeECZ2bcuvzzqGWtV1U0KktBSgRMpZ5oot3D4xk/8u38xpbRpz/6WdaX3YofFulkilowAlUk6278nlsfcX8uIXSuwqUh4UoEQOUnRi1wHdWvB/SuwqctAUoEQOwvKNQWLXKfPXctyR9ZXYVaQcKUCJlMG+vHye/2wpT0xZSDUz7rj4OAaf2orqSuwqUm4UoEQOUEb2Rm6bMJuFa7YrsatIDMX03z0z62VmC8wsy8yGFbK+lpm9Fq6fbmatItbdGi5fYGYXhMuONbNvIl5bzeymcN09ZrYiYt1FsTw2ST6bduxl2L++pd8z09i+O5fnBqbx3MA0BSeRGInZFZSZpQBPAecDOcAMM0sPh20vcC2wyd3bmFl/4GHgCjPrQDD8e0fgKOBDM2vn7guAEyLqXwFMiKjvcXf/a6yOSZKTu/Ovr1fwZyV2FalQsfwN6wZkufsSADMbB/QFIgNUX+CecHo8MMKCPrl9CYaY3wMsNbOssL5pEdueCyx292UxPAZJcllrt3HHxEy+XLKRri0a8sBlnTnuSKWgFKkIsQxQzYDlEfM5QPeiyoTZ0rcAjcPlX0Zt2yxq2/7A2KhlQ81sIJAB3Ozum6IbZWZDgCEALVq0OJDjkSSye18eI6Zm8eynSuwqEi+xfAZV2G+yl7JMsduaWU2gD/BGxPqngWMIbgGuAh4trFHuPtLd09w9rUmTJkW3XpLWJwvX0fPxTxnxURaXdDmKqbecxYDuLRScRCpYLK+gcoDmEfOpwMoiyuSYWXWgAbCxFNteCHzt7msKFkROm9lzwNvlcAySRNZu3c3wt+fy9rerlNhVJAHEMkDNANqaWWuCzgz9gQFRZdKBQQTPlvoBU93dzSwdeNXMHiPoJNGWH47meyVRt/fM7Eh3XxXOXgZklvPxSBWVl++8Mn0Zj7y3gD15+fzh/Hb86syjldhVJM5iFqDCZ0pDgclACvCCu88xs+FAhrunA6OAl8NOEBsJghhhudcJOlTkAje4ex6AmR1C0DPwV1G7/IuZnUBwKzC7kPUi+8lcsYXbJszm25wt/LjtYQzv20mJXUUShLlHPxZKHmlpaZ6RkRHvZkgcbNu9j8c+WMiYL7JpdGgt7rqkA5d0OVKJXUXiwMxmunta9HJ9kUOSirvzbuZq7n1rDmu37eHn3Vvwxwva06COEruKJBoFKEkayzfu5K5JmXy0YB0djqzPM784iROV2FUkYSlASZW3Ly+f5z5bwpNTFimxq0glogAlVdqM7I3cHiZ27dWxKXdd0oGjlDtPpFJQgJIqadOOvTz07nxey1hOs4Z1GDUojXOPOyLezRKRA6AAJVWKuzN+Zg5/fmce23bn8qszg8Suh9TUR12kstFvrVQZWWu3cduETL5aGiR2/fNPOtO+qRK7ilRWClBS6UUmdj2kZnUe+klnfpamxK4ilZ0ClFRqHy9Yy12T5vDdxp385MRm3HbxcRxWt1a8myUi5UABSiqlNWFi139/u4qjmxzKq7/szqnHKLGrSFWiACWVSl6+8/K0bP76/kL2KrGrSJWmACWVRnRi1/v6dqKVEruKVFkKUJLwtu3ex6PvL+SlaUFi1yf6n0Cf449SYleRKk4BShKWu/Ne5mruCRO7/qJ7S2654FgldhVJEgpQkpCiE7s+e1UaJzRvGO9miUgFUoCShLI3N0js+vepi0gx487eHRh0SksldhVJQjH9rTezXma2wMyyzGxYIetrmdlr4frpZtYqYt2t4fIFZnZBxPJsM5ttZt+YWUbE8kZm9oGZLQp/ahyFSuarpRu5+MnPeGTyAs5s14QPbz6Ta09vreAkkqRidgVlZinAUwTDs+cAM8ws3d3nRhS7Ftjk7m3MrD/wMHCFmXUgGP69I3AU8KGZtSsY9h04293XR+1yGDDF3R8Kg+Ew4E+xOj4pPxt37OWhd+fxekaOEruKyPdieYuvG5Dl7ksAzGwc0BeIDFB9gXvC6fHACAu6ZvUFxrn7HmCpmWWF9U0rZn99gbPC6THAxyhAJTR3542ZOTyoxK4iUohY/iVoBiyPmM8BuhdVxt1zzWwL0Dhc/mXUts3CaQfeNzMHnnX3keHyI9x9VVjXKjM7vLBGmdkQYAhAixYtynhocrAWrdnG7RMy+Sp7I2ktf8QDl3Xm2Kb14t0sEUkgsQxQhX1JxUtZprhtT3P3lWEA+sDM5rv7p6VtVBjQRgKkpaVFt0dibNfePP4+dREjP13CobWU2FVEihbLAJUDNI+YTwVWFlEmx8yqAw2AjcVt6+4FP9ea2QSCW3+fAmvM7Mjw6ulIYG35H5IcjI8WrOWuSZks37iLn3ZN5baL2tNYiV1FpAix7B41A2hrZq3NrCZBp4f0qDLpwKBwuh8w1d09XN4/7OXXGmgLfGVmh5pZPQAzOxToCWQWUtcgYFKMjksO0Jqtu7nhla+5evQMaqRUY+wve/Doz45XcBKRYsXsCip8pjQUmAykAC+4+xwzGw5kuHs6MAp4OewEsZEgiBGWe52gQ0UucIO755nZEcCEMMVNdeBVd38v3OVDwOtmdi3wHXB5rI5NSicyseu+vHxu6dmOX56hxK4iUjoWXLAkp7S0NM/IyCi5oBywb3M2c/uETGavCBK73n9pJ1o2VmJXEdmfmc1097To5erPK+Vq6+59PBYmdm1ctxZ/v/JEenc5UoldReSAKUBJuXB33pm9mnvfmsO67Xu4qkeQ2LV+bSV2FZGyUYCSg/bdhp3clZ7JxwvW0fGo+jw3MI3jldhVRA6SApSUWUFi1yenLKJGSjXu6t2BgUrsKiLlRAFKymT6kg3cPjGTrLXbubBTU+6+pCNNG9SOd7NEpApRgJIDsnHHXv78zjzGz8wh9Ud1eGFwGue0V2JXESl/ClBSKvn5zviZOTz4bpDY9ddnHcON57SlTk19p0lEYkMBSkq0cM027lBiVxGpYApQUqTIxK51a1fnLz/tQr+TUpXYVUQqhAKUFOqj+Wu5c1ImOZuU2FVE4kMBSn5g9ZbdDH97Du/MXk2bw+sybkgPehzdON7NEpEkpAAlQJDYdcwX2Tz6/gJy851berZjyBnHULO6vtMkIvGhACV8m7OZ2ybMJnPFVs5o14T7+nZUYlcRiTsFqCS2dfc+Hp28gJe+XEaTurUYMeBELu6sxK4ikhgUoJKQu/Pv2asY/tZc1m3fw8AeLblZiV1FJMEoQCWZZRt2cNekOXyycB2dmtXn+UFpdElVYlcRSTwxfQJuZr3MbIGZZZnZsELW1zKz18L1082sVcS6W8PlC8zsgnBZczP7yMzmmdkcM/tdRPl7zGyFmX0Tvi6K5bFVNntz8xkxdRE9H/+Umcs2cfclHZh0w+kKTiKSsGJ2BWVmKcBTwPlADjDDzNLdfW5EsWuBTe7exsz6Aw8DV5hZB4Lh3zsCRwEfmlk7guHfb3b3r82sHjDTzD6IqPNxd/9rrI6psvpyyQbuUGJXEalkYnmLrxuQ5e5LAMxsHNAXiAxQfYF7wunxwAgLntD3Bca5+x5gqZllAd3cfRqwCsDdt5nZPKBZVJ0S2rB9Dw++O//7xK6jB5/M2e0Pj3ezRERKJZYBqhmwPGI+B+heVBl3zzWzLUDjcPmXUds2i9wwvB14IjA9YvFQMxsIZBBcaW066KOohPLznTdmLufBd+ezXYldRaSSimWAKqyvspeyTLHbmlld4F/ATe6+NVz8NHBfWO4+4FHgmv0aZTYEGALQokWL4o+gElqweht3TJzNjOxNdGvViPsv60S7I5TYVUQqn1gGqBygecR8KrCyiDI5ZlYdaABsLG5bM6tBEJxecfc3Cwq4+5qCaTN7Dni7sEa5+0hgJEBaWlp0wKy0du3N44kpi3j+MyV2FZGqIZYBagbQ1sxaAysIOj0MiCqTDgwCpgH9gKnu7maWDrxqZo8RdJJoC3wVPp8aBcxz98ciKzKzI919VTh7GZAZo+NKOFPnr+GuSXPI2bSLy09K5daLjqPRoTXj3SwRkYMSswAVPlMaCkwGUoAX3H2OmQ0HMtw9nSDYvBx2gthIEMQIy71O0PkhF7jB3fPM7HTgKmC2mX0T7uo2d38H+IuZnUBwiy8b+FWsji1RrNqyi3vT5/LenCCx62tDetBdiV1FpIow9ypzl+uApaWleUZGRrybccBy8/IZM20Zj4WJXW88ty2//PHRSuwqIpWSmc1097To5cokUcl8s3wzt0+YzZyVWzmzXRPu69uJFo0PiXezRETKnQJUJbF19z7+OnkBL4eJXZ8a0JWLOjdVYlcRqbIUoBKcu/P2t6sY/vZcNmzfw6BTWnFzz3bUU2JXEaniFKAS2LINO7hz0hw+XbiOzs0aMEqJXUUkiShAJaA9uXmM/GQJf/8oi5op1bjnkg5cdUorUvSdJhFJIgpQCWba4g3cMXE2i9ft4OLOR3LXJR04or4Su4pI8lGAShAbtu/hgXfm8ebXK2jeqA6jrz6Zs49VYlcRSV4KUHEWmdh1x55cbjj7GIaercSuIiIKUHE0f/VW7piQScayTXRr3YgHLu1EWyV2FREBFKDiYufeXJ6YsohRny2lXu3qPNIvSOyq7zSJiPyPAlQFmzIvSOy6YrMSu4qIFEcBqoJEJnZtq8SuIiIlUoCKsejErn+84FgldhURKQUFqBiKTOx61rFNGN5HiV1FREpLASoGtuwKErv+c/oyDq9Xi3/8vCsXdlJiVxGRA6EAVY7cnbe+XcV9YWLXwae24g/nK7GriEhZKECVk+z1O7hzUiafLVpPl9QGvDDoZDqnNoh3s0REKq2YPqk3s15mtsDMssxsWCHra5nZa+H66WbWKmLdreHyBWZ2QUl1mlnrsI5FYZ0V0nd7T24ef5+yiJ5/+5RZ323m3j4dmfCb0xScREQOUswClJmlAE8BFwIdgCvNrENUsWuBTe7eBngceDjctgPQH+gI9AL+YWYpJdT5MPC4u7cFNoV1x9QXi9dz4ROf8egHCzm/wxFMuflMBp2qrOMiIuUhlldQ3YAsd1/i7nuBcUDfqDJ9gTHh9HjgXAt6EvQFxrn7HndfCmSF9RVaZ7jNOWEdhHVeGsNj45HJ8xnw3HRy85wXrz6ZpwZ0VdZxEZFyFMtnUM2A5RHzOUD3osq4e66ZbQEah8u/jNq2WThdWJ2Ngc3unltI+R8wsyHAEIAWLVoc2BFF6Na6MYYx9Jw21K6hxK4iIuUtlgGqsPtcXsoyRS0v7IqvuPL7L3QfCYwESEtLK7RMaZzZrglntmtS1s1FRKQEsbzFlwM0j5hPBVYWVcbMqgMNgI3FbFvU8vVAw7COovYlIiKVSCwD1Aygbdi7riZBp4f0qDLpwKBwuh8w1d09XN4/7OXXGmgLfFVUneE2H4V1ENY5KYbHJiIiMRazW3zhM6WhwGQgBXjB3eeY2XAgw93TgVHAy2aWRXDl1D/cdo6ZvQ7MBXKBG9w9D6CwOsNd/gkYZ2b3A7PCukVEpJKy4OIjOaWlpXlGRka8myEiktTMbKa7p0UvV0ptERFJSApQIiKSkBSgREQkISlAiYhIQkrqThJmtg5YdhBVHEbwHaxkluznINmPH3QOdPwHf/wt3X2/zAdJHaAOlpllFNbzJJkk+zlI9uMHnQMdf+yOX7f4REQkISlAiYhIQlKAOjgj492ABJDs5yDZjx90DnT8MaJnUCIikpB0BSUiIglJAUpERBKSAlQZmVkvM1tgZllmNize7Yk1M2tuZh+Z2Twzm2NmvwuXNzKzD8xsUfjzR/FuayyZWYqZzTKzt8P51mY2PTz+18JhYKosM2toZuPNbH74WTglmT4DZvb78POfaWZjzax2Vf8MmNkLZrbWzDIjlhX6nlvgyfDv4rdm1vVg9q0AVQZmlgI8BVwIdACuNLMO8W1VzOUCN7v7cUAP4IbwmIcBU9y9LTAlnK/KfgfMi5h/GHg8PP5NwLVxaVXFeQJ4z93bA8cTnIuk+AyYWTPgRiDN3TsRDPnTn6r/GXgR6BW1rKj3/EKC8fvaAkOApw9mxwpQZdMNyHL3Je6+FxgH9I1zm2LK3Ve5+9fh9DaCP0zNCI57TFhsDHBpfFoYe2aWClwMPB/OG3AOMD4sUtWPvz5wBuFYa+6+1903k0SfAYIx9OqEo3cfAqyiin8G3P1TgvH6IhX1nvcFXvLAlwQjnR9Z1n0rQJVNM2B5xHxOuCwpmFkr4ERgOnCEu6+CIIgBh8evZTH3N+D/gPxwvjGw2d1zw/mq/jk4GlgHjA5vcz5vZoeSJJ8Bd18B/BX4jiAwbQFmklyfgQJFvefl+rdRAapsrJBlSdFf38zqAv8CbnL3rfFuT0Uxs97AWnefGbm4kKJV+XNQHegKPO3uJwI7qKK38woTPmfpC7QGjgIOJbilFa0qfwZKUq6/EwpQZZMDNI+YTwVWxqktFcbMahAEp1fc/c1w8ZqCS/jw59p4tS/GTgP6mFk2wS3dcwiuqBqGt3ug6n8OcoAcd58ezo8nCFjJ8hk4D1jq7uvcfR/wJnAqyfUZKFDUe16ufxsVoMpmBtA27L1Tk+BBaXqc2xRT4fOWUcA8d38sYlU6MCicHgRMqui2VQR3v9XdU929FcH7PdXdfw58BPQLi1XZ4wdw99XAcjM7Nlx0LjCXJPkMENza62Fmh4S/DwXHnzSfgQhFvefpwMCwN18PYEvBrcCyUCaJMjKziwj+g04BXnD3B+LcpJgys9OBz4DZ/O8ZzG0Ez6FeB1oQ/AJf7u7RD1SrFDM7C7jF3Xub2dEEV1SNgFnAL9x9TzzbF0tmdgJBJ5GawBLgaoJ/dJPiM2Bm9wJXEPRqnQVcR/CMpcp+BsxsLHAWwbAaa4C7gYkU8p6HgXsEQa+/ncDV7p5R5n0rQImISCLSLT4REUlIClAiIpKQFKBERCQhKUCJiEhCUoASEZGEVL3kIiISK2aWR9B134A8YKi7f3GQdW5397rFrH8QmAw0BNq7+0MHsz+RWNEVlEh87XL3E9z9eOBW4MEK2Gd3gu+vnUnw3TaRhKQAJZI46hMM11Awrs4j4bhDs83sinD5P8ysTzg9wcxeCKevNbP7i6s8rO9b4GRgGsGXTJ82s7tieEwiZaZbfCLxVcfMvgFqA0cS5PgD+AlwAsGYS4cBM8zsU+BT4McEKWWahdsAnE6QzaBI7v5HM3sDuAr4A/Cxu59WvocjUn50BSUSXwW3+NoTpId5KUwXczow1t3z3H0N8AnBlc9nwI/DwSLn8r+knacApXl2dSLwDdA+3F4kYekKSiRBuPs0MzsMaELhwxbg7ivCYR96EVxNNQJ+BmwPB5IsVJhD70WC7NLrCQbbs/Dq7RR331WexyJSHnQFJZIgzKw9QfLhDQTB5wozSzGzJgQj2X4VFp0G3BSW+Qy4hRI6O7j7N+5+ArAQ6ABMBS4Ir94UnCQh6QpKJL4KnkFBcNU0yN3zzGwCwW27/xIM+PZ/4XAXEASjnu6eZWbLCK6iSuyNFwa6Te6eb2bt3V23+CShKZu5iIgkJN3iExGRhKQAJSIiCUkBSkREEpIClIiIJCQFKBERSUgKUCIikpAUoEREJCH9P+2nLqvs2zlaAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "prior.plot(label='prior', color='gray')\n", "posterior1.plot(label='posterior')\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 contain the most vanilla cookies.\n", "In between, the shape of the posterior distribution is a line because the 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": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.6699999999999999" ] }, "execution_count": 29, "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": 30, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwV5dn/8c9FQgj7EkBkBwGRHQyLVetuUVHQasUFcEWr1trl6aP+Wlvtpk9rXapVUUBwAXdLF6XuuCAQBGSHsElYwxYSIJDl+v0xE01jNiAn5yTn+369ziuz3DNzzZzJuc49c5+5zd0RERGJNXWiHYCIiEhplKBERCQmKUGJiEhMUoISEZGYpAQlIiIxSQlKRERikhKU1EhmlmNmXatpWyeb2epwm6OqY5s1jZm9ZWbjwuFrzOyTYvPczLpFL7pvmNnpZpZRxrxTzWxldcckZVOCknKZ2XozOxB+OG8zs8lm1ugo1tc5/MBKPJq43L2Ru689mnUchvuAx8JtvmlmH5rZDVW5gUisszq5+3nuPiXacRwNd//Y3Y+PdhzyDSUoqYwL3b0RMAgYDPwyWoEcbWI7wuU7AUuPZrslYkioqnWJ1GZKUFJp7r4JeAvoA2Bmbc1shpntMrN0M7uxqKyZDTGzNDPbG9a8/hLOmhX+3RPWyk4Ky19nZsvNbLeZzTSzTsXW5WZ2q5mtBlYXm9YtHG5qZlPNLNPMNpjZL82sTjjvGjP71MweMrNdwG9K7lcY62wz22NmW8zsMTNLCuetAboC/wjj/SNwKvBYOP5YWK6nmb0THouVZvaDYut/1syeMLN/m9k+4IwS2/99yXWa2b1m9tdwfl0z22dm/xeO1zezXDNrHo5fZGZLw/g/NLMTSnv/zOxJM/tziWl/N7OfhsN3mtkaM8s2s2VmdnGxcteY2Sdm9ufwPVpnZucVm1+pGqCZXWBmC8LzYqOZfev9KFF+pJktDMuvMbPh4fTyzr16ZvawmW0OXw+bWb0y1n97uK/trcTlv3Abr4Xn1Tozu73YvLLOb6lK7q6XXmW+gPXA2eFwB4KaxG/D8Y+AvwHJwAAgEzgrnDcbGBMONwKGhcOdAQcSi21jFJAOnAAkEtTQPis234F3gBZA/WLTuoXDU4G/A43D9a8Crg/nXQPkAz8K112/lH08ERgWzu8MLAfuKO0YhOMfAjcUG28IbASuDdcxCNgB9A7nPwtkAScTfClMLiWGkus8E1gcDn8HWAPMKTZvUTjcA9gHnAPUBX4RHsukUrbx3TBOC8ebAweAtuH4ZUDbMMbLw/UeW+w45gE3AgnAD4HNxdb1dfxh2U9KvH9F79XpQN9wG/2AbcCoMs69IeFxOycs3w7oWYlz7z7gc6A10Ar4jG/O2dOBjHD4V8AXQKtS5tUB5gP3AEkEX1LWAt8r7/zWq4o/f6IdgF6x/SL4cM4B9gAbwg+F+gTJqgBoXKzsH4Fnw+FZwL1AyxLr68y3E9RbhAklHK8D7Ac6heMOnFliPQ50Cz8sDwK9is27CfgwHL4G+Oow9/kO4I0Sx6C8BHU58HGJdTwF/DocfhaYWsE2S66zPpALpAB3AncDGeGH4b3Ao2G5XwEvlzh2m4DTS9mGAV8B3w3HbwTeLyemhcDIYscxvdi8BuF70KZk/JSToErZxsPAQ2XMe6q0eZU499YA5xeb9z1gfTh8enh8/gJ8AjQtVu50vklQQ0ueN8BdwOTyzm+9qvalS3xSGaPcvZm7d3L3W9z9AME37V3unl2s3AaCb7kA1xN8u19hZvPMbEQ56+8EPBJeotoD7CL4MG1XrMzGMpZtSfANd0MZcZS3LABm1sPM/mlmW81sL/CHcL2V1QkYWhR/uA9XAW0qG0NJ4TFOA04jqPl8RFATODmc9lFYtC3F9t3dC8NtFd//onkOTAeuCCddCbxQNN/MxoaX04r2oQ//fRy2FlvX/nDwsBrMmNlQM/sgvGyWBdxM2ce6A0GyKamic++/jkk43LbYeDNgPPBHd88qY9udgLYl3tO7gWPC+YdzfssRUoKSI7UZaGFmjYtN60jw7RR3X+3uVxBcZnkAeNXMGhJ8my5pI3BTmASLXvXd/bNiZcp67P4OgktPnYpN+zqOCpYt8gSwAuju7k0IPoisnPIl17cR+KhE/I3c/YeHEUNp8z8iuJw3EJgXjn+P4NJX0b28zRTbdzMzgg/2TZRuGnCpBff4hgKvhct1Ap4GbgNS3L0ZsITyj8OReBGYAXRw96bAk+VsYyNwXCnTyz33KHFMwnmbi43vBkYAk83s5HK2va7Ee9rY3c+Hcs9vqUJKUHJE3H0jwTf6P5pZspn1I/hW+QKAmV1tZq3Cb/R7wsUKCO4VFBJc0y/yJHCXmfUOl21qZpdVMo4C4GXg92bWOPyg/Snw/GHsTmNgL5BjZj0J7q+UZ1uJ+P8J9DCzMWGDhrpmNrisxgqVXCcECWkssMzdDxFeRiP44MwMy7wMXGBmZ5lZXeBnBJc8P6MU7r6A4D14Bpjp7kXvTdGXh0wAM7uWsDFMFWtMUPvJNbMhBLW4skwErg33rY6ZtTOznhWdewRJ+Jdm1srMWhLcR/qv88HdPySo5b5hZkNL2fZcYK+Z/W/YKCXBzPqY2WAo9/yWKqQEJUfjCoJ7SpuBNwjuubwTzhsOLDWzHOARYLS754aXhn4PfBpeOhnm7m8QfAudHl5iWwKcR+X9iOCG/lqC+wovApMOY/mfE3xQZhPUIl6qoPwjBLWQ3Wb2aHip6VxgNMGx2EqwP6W2HKvMOsNpnxHciyqqLS0juC9VNI67rwSuBv5KUJu8kOBnAYfK2dY04GyC41S0nmXAgwQ3/7cRNGT49DDir6xbgPvMLJsgcbxcVkF3n0vQ8OQhgsYSH/FNzai8c+93BJdHvwQWEzSE+F0p638nXP8MMzuxxLwCgmM5AFhHcGyfAZqGRUo9vyt7EKRyilrgiIiIxBTVoEREJCYpQYmISExSghIRkZikBCUiIjHpqB68WVO0bNnSO3fuHO0wRESkFPPnz9/h7q1KTo+LBNW5c2fS0tKiHYaIiJTCzDaUNl2X+EREJCYpQYmISEyKaIIys+EW9I2TbmZ3ljK/npm9FM6fY2adw+lDwodWLjSzRfbf/dKUu04REakdInYPyoJeQx8n6MslA5hnZjPCR6oUuR7Y7e7dzGw0weNhLid41E2qu+eb2bHAIjP7B8GzwipaZ6Xk5eWRkZFBbq6eTnK4kpOTad++PXXr1o12KCJSi0WykcQQgv5j1gKY2XRgJMHzxIqM5JseTl8l6FHUij3KH4IOyYqex1SZdVZKRkYGjRs3pnPnzgQPgJbKcHd27txJRkYGXbp0iXY4IlKLRfISXzv+uw+cDL7dR83XZdw9n+CBkCnwdb8xSwke9nhzOL8y6yRcfnzYJXNaZmbmt+bn5uaSkpKi5HSYzIyUlBTVPEUk4iKZoEr75C/5ZNoyy7j7HHfvDQwm6IohuZLrJFx+grununtqq1bfal4fbFzJ6YjouIlIdYhkgsog6DitSHv+u9Ow/ypjZokEj7LfVbyAuy8n6EqhTyXXKSIi1WDDzn3MWBS5j+BIJqh5QHcz62JmSQR95cwoUWYGMC4cvhR43909XCYRvu7p83hgfSXXGVfefPNNli077FtwzJgxg/vvvz8CEYlIPMjOzeP6KWnc8/clZO3Pi8g2IpagwntGtwEzgeXAy+6+1MzuM7OLwmITgRQzSyfoBbWo2fgpBC33FhJ0RnaLu+8oa52R2oea4EgSVH5+PhdddBF33ln5Vvr5+fmHG5qI1FIFhc7t0xawfsc+/nbVIJo2iEyL3rjosDA1NdVLPupo+fLlnHDC4fTIXfXWr1/P8OHDGTp0KAsWLKBHjx5MnTqV2bNn8/Of/5z8/HwGDx7ME088Qb169bjzzjuZMWMGiYmJnHvuuVxyySWMGDGCpk2b0rRpU1577TUAbr31VjIzM2nQoAFPP/00PXv25JprrqFFixYsWLCAQYMG0bdvX9LS0njsscfYsGED1113HZmZmbRq1YrJkyfTsWPHby3z4IMPfh17LBw/EYmOP/57OU/NWstvR/VhzLBOFS9QATOb7+6pJafHxbP4KnLvP5aybPPeKl1nr7ZN+PWFvSsst3LlSiZOnMjJJ5/Mddddx1/+8heeeuop3nvvPXr06MHYsWN54oknGDt2LG+88QYrVqzAzNizZw/NmjXjoosuYsSIEVx66aUAnHXWWTz55JN0796dOXPmcMstt/D+++8DsGrVKt59910SEhJ49tlnv47htttuY+zYsYwbN45JkyZx++238+abb35rGRGRV+dn8NSstYwZ1qlKklN59KijKOvQoQMnn3wyAFdffTXvvfceXbp0oUePHgCMGzeOWbNm0aRJE5KTk7nhhht4/fXXadCgwbfWlZOTw2effcZll13GgAEDuOmmm9iyZcvX8y+77LJSE83s2bO58sorARgzZgyffPJJhcuISPyZv2EXd7++mJO6pnDPhb0ivj3VoKBSNZ1IqWyT7cTERObOnct7773H9OnTeeyxx76uGRUpLCykWbNmLFy4sNR1NGzY8LBjquwyIlK7Zezez/ip82nbLJknrh5E3YTI129Ug4qyr776itmzZwMwbdo0zj77bNavX096ejoAzz33HKeddho5OTlkZWVx/vnn8/DDD3+dhBo3bkx2djYATZo0oUuXLrzyyitA8NSHRYsWVRjDd77zHaZPnw7ACy+8wCmnnFLl+ykiNVfOwXxumJLGoYJCJl4zmGYNkqplu0pQUXbCCScwZcoU+vXrx65du/jJT37C5MmTueyyy+jbty916tTh5ptvJjs7mxEjRtCvXz9OO+00HnroIQBGjx7Nn/70JwYOHMiaNWt44YUXmDhxIv3796d37978/e9/rzCGRx99lMmTJ9OvXz+ee+45HnnkkUjvtojUEAWFzh3TF7B6ew5/u2oQx7VqVG3bViu+KFq/fj0jRoxgyZIlUY3jSMTC8RORyPvDv5czYdZa7hvZm7EndY7INspqxacalIiIlGr63K+YMGst407qFLHkVB4lqCjq3Llzjaw9iUjt99maHfzyzSV8t0crfjUi8i32ShPXCSoeLm9Ggo6bSO22bsc+fvj8F3Rp2ZDHrhxIYjW02CtN3Cao5ORkdu7cqQ/bw1TUH1RycnK0QxGRCNi97xDXPTuPhDrGxHGDaZIcvY5J4/Z3UO3btycjI4PS+oqS8hX1qCsitcuh/EJuen4+m/YcYNqNQ+mY8u0HAlSnuE1QdevWVY+wIiIhd+eu1xczd90uHhk9gBM7tYh2SPF7iU9ERL7x+AfpvPZFBj89pwcjB5TaUXm1U4ISEYlzf1+4iT//ZxWXDGzHj87sFu1wvqYEJSISx+at38X/vPIlQ7u04I/f71vp54NWByUoEZE4tW7HPm6cmkb75vV5asyJ1EuMrZ4LlKBEROLQrn2HuHbyXOqYMfna6nsA7OGI21Z8IiLxKjevgBunprE5K5dpNw6lU0psdqujGpSISBwpLHR+9vIi5m/YzcOXx0Zz8rIoQYmIxJH7317BvxZv4e7ze3J+32OjHU65lKBEROLE1NnrmTBrLWOGdeLGU7tGO5wKKUGJiMSB/yzdym9mLOXsE1rz6wt7xVRz8rIoQYmI1HILvtrN7dMX0LddUx69InpPJz9cNSNKERE5Iut37OP6KWm0bpzMxGsG0yCp5jTeVoISEamlduQcZNzkubg7z147mJaN6kU7pMMS0QRlZsPNbKWZpZvZnaXMr2dmL4Xz55hZ53D6OWY238wWh3/PLLbMh+E6F4av1pHcBxGRmmj/oXyuf3Ye2/bmMvGawXRt1SjaIR22iNX1zCwBeBw4B8gA5pnZDHdfVqzY9cBud+9mZqOBB4DLgR3Ahe6+2cz6ADOB4o/Xvcrd0yIVu4hITZZfUMhtLy5g8aYsnhqTyqCOzaMd0hGJZA1qCJDu7mvd/RAwHRhZosxIYEo4/CpwlpmZuy9w983h9KVAspnVrLqpiEgUuDu/fHMJ76/Yzm9H9eGcXsdEO6QjFskE1Q7YWGw8g/+uBf1XGXfPB7KAlBJlvg8scPeDxaZNDi/v/crKaCtpZuPNLM3M0tRrrojEi4feXc30eRu57YxuXDW0U7TDOSqRTFClJQ4/nDJm1pvgst9NxeZf5e59gVPD15jSNu7uE9w91d1TW7VqdViBi4jURM9/voFH31vND1Lb87Nze0Q7nKMWyQSVAXQoNt4e2FxWGTNLBJoCu8Lx9sAbwFh3X1O0gLtvCv9mAy8SXEoUEYlrby/Zyj1/X8JZPVvzh4tjq1+nIxXJBDUP6G5mXcwsCRgNzChRZgYwLhy+FHjf3d3MmgH/Au5y90+LCptZopm1DIfrAiOAJRHcBxGRmDdn7U5un76A/h2a8diVg2rMD3ErErG9CO8p3UbQAm858LK7LzWz+8zsorDYRCDFzNKBnwJFTdFvA7oBvyrRnLweMNPMvgQWApuApyO1DyIisW75lr3cMDWNDs3rM2ncYOonxVang0fD3EveFqp9UlNTPS1NrdJFpHbZuGs/lzzxGQlmvHbLd2jXrH60QzoiZjbf3VNLTq85z7wQEZGv7cw5yNhJczmUX8grN59UY5NTeWrHhUoRkTiSczCfaybPY0vWASZdk0qPYxpHO6SIUA1KRKQGyc0rYPzUNJZt2cvTY0+M6R5xj5ZqUCIiNURBoXPH9IV8tmYnf76sH2f2rLlPiagMJSgRkRogeITRYt5eupV7RvTi4oHtox1SxClBiYjUAA+8vZJpc4NHGF13Spdoh1MtlKBERGLcEx+u4cmP1jBmWKda8QijylKCEhGJYdPmfsUDb6/gov5tufei3rXiEUaVpQQlIhKj/rFoM3e/sZgzjm/Fgz/oT5068ZOcQAlKRCQmvb9iGz95aSGDO7fgb1edSN1a8ny9wxF/eywiEuM+X7uTHz7/BScc24SJ41Jr1fP1DocSlIhIDFm0cQ83TEmjQ4sGTLluCI2T60Y7pKhRghIRiRErtu5l7KS5NG9Yl+evH0qLhknRDimqlKBERGLA2swcrn5mLvXrJvDiDcNo0zQ52iFFnRKUiEiUZezez9XPzMHdef6GoXRo0SDaIcUEJSgRkSjampXLlU/PIedgPs9dP5RurRtFO6SYoaeZi4hEyY6cg1z1zOfs2neI528YSq+2TaIdUkxRDUpEJAr27D/E1c/MYfOeXCZdM5gBHZpFO6SYoxqUiEg1yzqQx5iJc1m7Yx+Txg1mSJfa26fT0VANSkSkGgW94c5lxda9PHn1IE7p3jLaIcUs1aBERKrJ/kP5XDt5Ll9mZPH4lYNqfYeDR0s1KBGRanDgUAE3TElj/obdPDJ6AMP7tIl2SDFPNSgRkQjLzSvgxqlpzF67kwcv68+Ifm2jHVKNoBqUiEgE5eYVcNNz8/l0zQ7+7/v9uGRQ7e+qvaooQYmIRMjB/AJueeELPlqVyf2X9OWy1A7RDqlGiWiCMrPhZrbSzNLN7M5S5tczs5fC+XPMrHM4/Rwzm29mi8O/ZxZb5sRwerqZPWrx1L2kiNQYB/MLuPWFL3h/xXb+cHFfLh/cMdoh1TgRS1BmlgA8DpwH9AKuMLNeJYpdD+x2927AQ8AD4fQdwIXu3hcYBzxXbJkngPFA9/A1PFL7ICJyJA7lF3LrC1/w7vLt/HZUH64cquR0JCJZgxoCpLv7Wnc/BEwHRpYoMxKYEg6/CpxlZubuC9x9czh9KZAc1raOBZq4+2x3d2AqMCqC+yAiclgO5Rdyywvzv05OY4Z1inZINVYkE1Q7YGOx8YxwWqll3D0fyAJSSpT5PrDA3Q+G5TMqWCcAZjbezNLMLC0zM/OId0JEpLKK7jm9u3w7vx3ZW8npKEUyQZV2b8gPp4yZ9Sa47HfTYawzmOg+wd1T3T21VatWlQhXROTIHcwv4Jbnv+Dd5du4b2RvxpzUOdoh1XiRTFAZQPEmK+2BzWWVMbNEoCmwKxxvD7wBjHX3NcXKF2+jWdo6RUSqVW5eATc/N5/3Vmznd6P6MFbJqUpEMkHNA7qbWRczSwJGAzNKlJlB0AgC4FLgfXd3M2sG/Au4y90/LSrs7luAbDMbFrbeGwv8PYL7ICJSrqLfOX2wMpM/XNyXq3VZr8pELEGF95RuA2YCy4GX3X2pmd1nZheFxSYCKWaWDvwUKGqKfhvQDfiVmS0MX63DeT8EngHSgTXAW5HaBxGR8hQ9vmjW6uB3TmqtV7UsaAxXu6WmpnpaWlq0wxCRWmT/oXyufzaNz9ft5E+X9ufSE/WEiCNlZvPdPbXkdD2LT0TkMOUczOe6yfNI27CLh34wgFEDS21MLEdJCUpE5DBkHcjjmrDLjEdGD+TC/nrwa6QoQYmIVNLufYcYM2kOK7dm8/iVg9RlRoQpQYmIVEJm9kHGTJzD2h37mDAmlTN6tq54ITkqSlAiIhXYknWAq56Zw5Y9uUy+ZjAnd1M37dVBCUpEpBxf7dzPlc98Ttb+PKZeP4TBnVtEO6S4oQQlIlKG9O05XPXM5xzML+SFG4fSr32zaIcUV5SgRERKsWRTFmMnzaWOGdPHD6NnmybRDinuKEGJiJSQtn4X1z47jybJdXn+hqF0adkw2iHFJSUoEZFiPl6dyfip8zm2aTLP3TCUds3qRzukuKUEJSISenvJFm6ftpDjWjdi6nVDaNW4XrRDimtKUCIiwMvzNnLn618ysGNzJo0bTNMGdaMdUtxTghKRuPf0rLX8/t/LOa1HK564ehANkvTRGAv0LohI3HJ3/jRzJX/7cA0X9DuWh34wgKTESHaTJ4dDCUpE4lJ+QSG/fHMJ0+dt5MqhHfntyD4k1LFohyXFKEGJSNzJzSvgx9MXMHPpNn50Zjd+ek4Pgk66JZYoQYlIXNmbm8f4qWl8vnYX94zoxXWndIl2SFIGJSgRiRvb9+YybvI8Vm/L5uHL1dFgrFOCEpG4sG7HPsZOmsPOnENMvGYwp/VoFe2QpAJKUCJS6y3auIfrnp2HAy/eOIwBHfTQ15pACUpEarUPVm7nlue/IKVRElOuG8JxrRpFOySpJCUoEam1XknbyJ2vL6Znm8ZMvnYwrRsnRzskOQzl/iLNzJ4tNjwu4tGIiFQBd+fR91bzP69+yXeOS+Glm05ScqqBKvrJdP9iwz+OZCAiIlUhv6CQu15fzF/eWcUlA9sxcdxgGtXTxaKaqKIE5UezcjMbbmYrzSzdzO4sZX49M3spnD/HzDqH01PM7AMzyzGzx0os82G4zoXhq/XRxCgitce+g/ncODWN6fM2ctsZ3XjwB/316KIarKKvFe3N7FHAig1/zd1vL2tBM0sAHgfOATKAeWY2w92XFSt2PbDb3buZ2WjgAeByIBf4FdAnfJV0lbunVRC7iMSR7XtzuW7KPJZt3svvL+7DVUM7RTskOUoVJaj/KTZ8uAlhCJDu7msBzGw6MBIonqBGAr8Jh18FHjMzc/d9wCdm1u0wtykicWjVtmyunTyP3fsPMXHcYM7oqQsrtUG5CcrdpxzFutsBG4uNZwBDyyrj7vlmlgWkADsqWPdkMysAXgN+5+5HdSlSRGquT9N3cPNz86mflMDLN51En3ZNox2SVJFyE5SZzShvvrtfVN7ipS1yBGVKusrdN5lZY4IENQaY+q2Nm40HxgN07NixglWKSE308ryN3P3GYrq2asjka4eoe/ZapqJLfCcR1HCmAXMoPaGUJQPoUGy8PbC5jDIZZpYINAV2lbdSd98U/s02sxcJLiV+K0G5+wRgAkBqaqpqWCK1SGGh8+f/BP04ndq9JY9fNYgmyeoBt7apqHlLG+BugoYKjxA0eNjh7h+5+0cVLDsP6G5mXcwsCRgNlKyRzQCKfl91KfB+eZfrzCzRzFqGw3WBEcCSCuIQkVokN6+AH01bwN8+XMOVQzsy6ZrBSk61VEX3oAqAt4G3zawecAXwoZnd5+5/rWDZfDO7DZgJJACT3H2pmd0HpLn7DGAi8JyZpRPUnEYXLW9m64EmQJKZjQLOBTYAM8PklAC8Czx9BPstIjXQ9uxcbpw6ny8z9nD3+T258dSu6sepFrOK2heEiekCguTUmaDWM6noUltNkJqa6mlpapUuUpMt27yXG6bMY/f+PB4ePYDv9W4T7ZCkipjZfHdPLTm9okYSUwgu770F3OvuupwmItXu3WXb+PH0BTROrssrN6ulXryoqJHEGGAf0AP4sZkVVbcMcHdvEsngRCS+uTtPzVrLA2+voE/bpjwzLpVjmuiZevGiontQekaIiETFwfwC7np9Ma9/sYkL+h3Lny/tT/2khGiHJdWookt8ycDNQDfgS4J7T/nVEZiIxK/t2bnc/Nx8vvhqD3ec3Z0fn9VdjSHiUEWX+KYAecDHwPlAb/RUcxGJoCWbsrhxahq79x/i8SsHcUG/Y6MdkkRJRQmql7v3BTCzicDcyIckIvHqn19u5uevLKJFgyRevfk7agwR5ypKUHlFA+HvmiIcjojEo4JC58HwyRAndmrOk1efSKvG9aIdlkRZRQmqv5ntDYcNqB+OqxWfiFSJrAN53DF9AR+szOSKIR34zUW9qZeoxhBScSs+nSUiEjHp27MZP3U+X+3az+9G9eHqYerDSb6hfpBFJCreXrKVn728kPpJCbxww1CGdk2JdkgSY5SgRKRaFRQ6D7+7ir++n07/Ds148upBHNtU3WTItylBiUi12bP/ED+evpCPVmVyeWoH7h3Zm+S6upMgpVOCEpFqsWRTFjc/P5/tew/y+4v7cOWQjvrxrZRLCUpEIu7V+Rn8vzcW06JhEi/ffBIDOjSLdkhSAyhBiUjE5OYVcO8/ljFt7lec1DWFv145kJaN9PsmqRwlKBGJiI279vPDF+azZNNebjn9OH56Tg8SE/T8aak8JSgRqXLvLd/GT19eRKE7T49N5Zxex0Q7JKmBlKBEpMrkFxTyp/+s5KmP1tK7bRP+dtUgOqU0jHZYUkMpQYlIldialcvt0xcwd90urhzakXtG9FITcjkqSlAictQ+WpXJT15aSG5eAQ9fPoBRA9tFOySpBZSgROSI5RcU8uA7q3jiw5q9BB8AABW0SURBVDX0bNOYx64cRLfWjaIdltQSSlAickQ27TnAj6ctIG3Dbq4Y0oFfX6inQkjVUoISkcM2c+lWfvHqlxQUOo+MHsDIAbqkJ1VPCUpEKi03r4A//ns5U2ZvoG+7pvz1ioF0bqlWehIZSlAiUimrt2Xzo2kLWLE1m+tO7sL/nne8OhaUiIroz7rNbLiZrTSzdDO7s5T59czspXD+HDPrHE5PMbMPzCzHzB4rscyJZrY4XOZR09MmRSLK3XlhzgYufOwTMrMPMvnawdxzYS8lJ4m4iCUoM0sAHgfOA3oBV5hZrxLFrgd2u3s34CHggXB6LvAr4OelrPoJYDzQPXwNr/roRQRg175DjH9uPv/vjSUM7tyCt+44lTOObx3tsCRORLIGNQRId/e17n4ImA6MLFFmJDAlHH4VOMvMzN33ufsnBInqa2Z2LNDE3We7uwNTgVER3AeRuDVrVSbfe3gWH63M5JcXnMCUa4fQunFytMOSOBLJe1DtgI3FxjOAoWWVcfd8M8sCUoAd5awzo8Q6S20+ZGbjCWpadOzY8XBjF4lbuXkFPPD2CiZ/up7urRsx5doh9GrbJNphSRyKZIIq7d6QH0GZIyrv7hOACQCpqanlrVNEQks2ZfGTlxayensO407qxF3nn6DfNknURDJBZQAdio23BzaXUSbDzBKBpsCuCtbZvoJ1ishhKih0npq1hofeWUXzBklMuW4Ip/VoFe2wJM5FMkHNA7qbWRdgEzAauLJEmRnAOGA2cCnwfnhvqVTuvsXMss1sGDAHGAv8NRLBi8SLDTv38dOXFzF/w27O69OGP1zcl+YNk6IdlkjkElR4T+k2YCaQAExy96Vmdh+Q5u4zgInAc2aWTlBzGl20vJmtB5oASWY2CjjX3ZcBPwSeBeoDb4UvETlM7s6Lc7/i9/9aTkId4+HLBzByQFv0yw2JFVZOhaXWSE1N9bS0tGiHIRIzNu85wP++9iUfr97BKd1a8qfL+nFs0/rRDkvilJnNd/fUktP1JAmROOLuvJKWwW//uYwCd347qg9XD+2oWpPEJCUokTixec8B7np9MR+tymRIlxb8+dL+dExpEO2wRMqkBCVSy7k7L6dt5Hf/XE5+ofPrC3sx7qTO1KmjWpPENiUokVps46793P3GYj5evYNhXVvwwPf70SlFTx+XmkEJSqQWKix0psxez59mrsSA347szVVDO6nWJDWKEpRILbN6WzZ3vr6Y+Rt2c1qPVvzhkr60a6YWelLzKEGJ1BIH8wv42wdr+NuH6TSql8hfftCfiwe2Uws9qbGUoERqgbnrdnH3G4tJ357DqAFt+dWIXqQ0qhftsESOihKUSA22Z/8h/vjvFbyUtpH2zevz7LWDOV39NUktoQQlUgO5O28s2MTv/7WcPQfyuOm0rtxxVg/qJ+nJ41J7KEGJ1DDp23P45ZuL+XztLgZ0aMbzl/TlhGPVX5PUPkpQIjXE/kP5PP5BOhNmraVBUiJ/uLgvowd3UNNxqbWUoERinLszc+lW7vvHMjZn5XLJoHbcff4JtFQjCKnllKBEYtiazBzu/ccyZq3KpGebxjxyxUAGd24R7bBEqoUSlEgMyjmYz1/fW82kT9eRnJjAPSN6MfakTiQm1Il2aCLVRglKJIYUFgat8x54ewXbsw/yg9T2/GJ4T13Ok7ikBCUSI774ajf3/mMZizbuoX+HZkwYm8qADs2iHZZI1ChBiUTZ5j0H+L+3V/Dmws20blyPv/ygP6MGtFPrPIl7SlAiUbLvYD5PfrSGCbPW4sCtZxzHLad3o2E9/VuKgBKUSLUrKHReSdvIg++sIjP7IBf1b8svhh9P++bq3VakOCUokWri7ny4KpP7/72ClduyObFTc54acyKDOjaPdmgiMUkJSqQaLNq4h/vfWsHstTvpnNKAJ68exPd6t1FXGCLlUIISiaC1mTk8+J9V/GvxFlIaJnHvRb25YkhHkhL1eyaRiihBiUTAlqwDPPreal5Oy6BeYh1uP6s747/blUZqACFSafpvEalCO3MO8tSstUz5bD2F7owZ1olbz+hGq8b6oa3I4YpogjKz4cAjQALwjLvfX2J+PWAqcCKwE7jc3deH8+4CrgcKgNvdfWY4fT2QHU7Pd/fUSO6DSGVkHcjjmY/XMumTdRzIK2DUwHb85OwedGihlnkiRypiCcrMEoDHgXOADGCemc1w92XFil0P7Hb3bmY2GngAuNzMegGjgd5AW+BdM+vh7gXhcme4+45IxS5SWXtz85j8yXomfrKWvbn5XNDvWH5ydne6tW4c7dBEarxI1qCGAOnuvhbAzKYDI4HiCWok8Jtw+FXgMQuaNY0Eprv7QWCdmaWH65sdwXhFKi07N48pn63n6Y/XkXUgj3N7HcOPz+5O77ZNox2aSK0RyQTVDthYbDwDGFpWGXfPN7MsICWc/nmJZduFww78x8wceMrdJ5S2cTMbD4wH6Nix49HtiUgoa38ekz9bx6RP1rE3N5+zT2jNHWf3oE87JSaRqhbJBFXaDzy8kmXKW/Zkd99sZq2Bd8xshbvP+lbhIHFNAEhNTS25XZHDsjPnIJM+XcfUzzaQfTCfc3sdw+1ndVdiEomgSCaoDKBDsfH2wOYyymSYWSLQFNhV3rLuXvR3u5m9QXDp71sJSqQqbM3KZcKstUyb+xW5+QUM792G287spkt5ItUgkglqHtDdzLoAmwgaPVxZoswMYBzBvaVLgffd3c1sBvCimf2FoJFEd2CumTUE6rh7djh8LnBfBPdB4tSazBwmfLSW1xdkUOgwckBbbjn9ODV+EKlGEUtQ4T2l24CZBM3MJ7n7UjO7D0hz9xnAROC5sBHELoIkRljuZYIGFfnAre5eYGbHAG+Ej4dJBF5097cjtQ8Sf774ajcTPlrLzGVbSUqowxVDOnLjqV3VXFwkCsy99t+eSU1N9bS0tGiHITGqsNB5b8V2Jsxaw7z1u2mSnMiYkzpx7cld1JOtSDUws/ml/aZVT5KQuLX/UD6vzc9g0qfrWbdjH+2a1eeeEb24fHAH9ckkEgP0XyhxZ9OeA0ydvZ7pczeSdSCP/h2a8egVAzm/TxsSE/QQV5FYoQQlccHdmbtuF89+tp6ZS7cC8L3ebbjh1C4M6thc3V6IxCAlKKnV9h3M582Fm3hu9gZWbM2maf263PjdrowZ1kk92IrEOCUoqZVWbcvmxTlf8doXGWTn5tPr2Cbcf0lfRg5oR/2khGiHJyKVoAQltUZuXgFvLdnCtDkbmbt+F0kJdTivbxvGDOvEiZ10GU+kplGCkhpvxda9TJ+7kde/yGBvbj6dUhpw13k9ufTE9qSombhIjaUEJTVS1oE8ZizazCtpG/kyI4ukhDp8r08brhjcgWFdU6hTR7UlkZpOCUpqjPyCQj5evYNXv8jgnWXbOJRfSM82jblnRC9GDWxHi4ZJ0Q5RRKqQEpTENHdn6ea9vLFgEzMWbSYz+yDNGtRl9OAOXHpie/q2a6p7SyK1lBKUxKT1O/bxzy838+bCzaRvz6FugnH68a35/qD2nNGzFfUS1RJPpLZTgpKYsWnPAd5avIV/LNrMoowsAAZ3bs7vL+7DBX2PpVkDXcITiSdKUBJVG3ftZ+bSrfxr8RYWfLUHgD7tmnD3+T0Z0a8tbZvVj3KEIhItSlBSrdyd9O05/GfZNt5asoUlm/YC0LttE34x/Hgu6HssnVIaRjlKEYkFSlAScfkFhSzYuId3lm3jnWXbWLdjHwADOzbjrvN6MrxPGyUlEfkWJSiJiN37DjFrdSYfrNjOh6sy2bM/j7oJxknHteS6U7pwzgnH0KZpcrTDFJEYpgQlVaKg0Fm8KYtZqzL5cOV2Fm7cQ6FDi4ZJnNXzGM46oTWndm9J4+S60Q5VRGoIJSg5Iu7Ohp37+WzNTj5N38En6TvIOpCHGfRr15TbzuzO6ce3on/7ZiToqQ4icgSUoKTSMnbvZ87aXXy+diefrdnJpj0HAGjTJJlzex3DqT1acfJxKXr+nYhUCSUoKVVhobN2Rw5z1+0mbf0u5q7fRcbuICE1b1CXoV1SuPm0rnynW0u6tmyopzmISJVTghIA9ubmsTgjiwVf7Wb+ht188dUesg7kAdCyURKDO7fghlO6MOy4FHq0bqyHsYpIxClBxaEDhwpYtiWLLzOyWLwpi0Ub97Amc9/X87u3bsR5fdowqGNzUjs3p4tqSCISBUpQtZi7s23vQVZs3cuKrdks27yXpZuzWLdjH4UelGnVuB792jVl1IB2DOjYjH7tmtG0gVraiUj0KUHVAoWFzta9uazN3Mfq7dmkb89h9fYcVm7N/voyHUC7ZvU54dgmXNCvLX3bNaVf+6Yc00S/RRKR2KQEVUPkFxSyJSuXjbv3s3HXftbv3M+GnftYv2M/63bs40BewddlmyQn0uOYxpzf91h6tmnM8W0a07NNYz1sVURqlIgmKDMbDjwCJADPuPv9JebXA6YCJwI7gcvdfX047y7geqAAuN3dZ1ZmnTXRgUMF7Mg5yPbsg2zfm8vW8LVlTy5bsg6weU8wXlB0XQ6om2B0aN6AjikNGNY1ha6tGtK1ZUO6HdOIVo3q6Z6RiNR4EUtQZpYAPA6cA2QA88xshrsvK1bsemC3u3czs9HAA8DlZtYLGA30BtoC75pZj3CZitZZbdydg/mFHMwv5FB+Ibl5Bew/VMD+Q/nsP1RAdm4+OQfzyc7NI+vAN689+/PYte8Qu/cfYmfOIXIO5n9r3XUTjDZNk2nbtD5DurSgXbP6dGhRnw7NG9ChRQPaNquvH8CKSK0WyRrUECDd3dcCmNl0YCRQPJmMBH4TDr8KPGbBV/+RwHR3PwisM7P0cH1UYp1V6q7XF/PWki0UFjpOcL8nP3wVr9FURqN6iTStX5fmDevSvEESnVIa0LxBEq0a1wtejepxTJNkjmlSj+YNktSUW0TiWiQTVDtgY7HxDGBoWWXcPd/MsoCUcPrnJZZtFw5XtE4AzGw8MB6gY8eOR7YHwKCOzUhKMMwMMzCMuglGYoKRWKcOSYl1qBe+kusm0CApkQZJCdRPSqBxciKN69WlYb0EmtSvS92EOkcch4hIvIlkgirt63/JKkdZZcqaXtonfKnVGHefAEwASE1NPbyqTjGXpXbgstQOR7q4iIgcoUh+pc8Ain+ytwc2l1XGzBKBpsCucpatzDpFRKQWiGSCmgd0N7MuZpZE0OhhRokyM4Bx4fClwPvu7uH00WZWz8y6AN2BuZVcp4iI1AIRu8QX3lO6DZhJ0CR8krsvNbP7gDR3nwFMBJ4LG0HsIkg4hOVeJmj8kA/c6u4FAKWtM1L7ICIi0WNBhaV2S01N9bS0tGiHISIipTCz+e6eWnK6mpWJiEhMUoISEZGYpAQlIiIxSQlKRERiUlw0kjCzTGDDUayiJbCjisKpyXQcdAyK6DgEdByq5hh0cvdWJSfGRYI6WmaWVloLk3ij46BjUETHIaDjENljoEt8IiISk5SgREQkJilBVc6EaAcQI3QcdAyK6DgEdBwieAx0D0pERGKSalAiIhKTlKBERCQmKUFVwMyGm9lKM0s3szujHU91MLMOZvaBmS03s6Vm9uNwegsze8fMVod/m0c71upgZglmtsDM/hmOdzGzOeFxeCns+qXWMrNmZvaqma0Iz4mT4vFcMLOfhP8PS8xsmpklx8O5YGaTzGy7mS0pNq3U998Cj4afl1+a2aCj2bYSVDnMLAF4HDgP6AVcYWa9ohtVtcgHfubuJwDDgFvD/b4TeM/duwPvhePx4MfA8mLjDwAPhcdhN3B9VKKqPo8Ab7t7T6A/wbGIq3PBzNoBtwOp7t6HoLuf0cTHufAsMLzEtLLe//MI+u/rDowHnjiaDStBlW8IkO7ua939EDAdGBnlmCLO3be4+xfhcDbBB1I7gn2fEhabAoyKToTVx8zaAxcAz4TjBpwJvBoWqdXHwcyaAN8l6LsNdz/k7nuIw3OBoP+8+mHv3w2ALcTBueDuswj66yuurPd/JDDVA58Dzczs2CPdthJU+doBG4uNZ4TT4oaZdQYGAnOAY9x9CwRJDGgdvciqzcPAL4DCcDwF2OPu+eF4bT8nugKZwOTwMuczZtaQODsX3H0T8GfgK4LElAXMJ77OheLKev+r9DNTCap8Vsq0uGmXb2aNgNeAO9x9b7TjqW5mNgLY7u7zi08upWhtPicSgUHAE+4+ENhHLb+cV5rwHstIoAvQFmhIcDmrpNp8LlRGlf5/KEGVLwPoUGy8PbA5SrFUKzOrS5CcXnD318PJ24qq6+Hf7dGKr5qcDFxkZusJLu+eSVCjahZe5oHaf05kABnuPiccf5UgYcXbuXA2sM7dM909D3gd+A7xdS4UV9b7X6WfmUpQ5ZsHdA9b6iQR3BSdEeWYIi68zzIRWO7ufyk2awYwLhweB/y9umOrTu5+l7u3d/fOBO/9++5+FfABcGlYrFYfB3ffCmw0s+PDSWcBy4izc4Hg0t4wM2sQ/n8UHYe4ORdKKOv9nwGMDVvzDQOyii4FHgk9SaICZnY+wbfmBGCSu/8+yiFFnJmdAnwMLOabey93E9yHehnoSPAPe5m7l7x5WiuZ2enAz919hJl1JahRtQAWAFe7+8FoxhdJZjaAoJFIErAWuJbgy21cnQtmdi9wOUEr1wXADQT3V2r1uWBm04DTCbrV2Ab8GniTUt7/MHk/RtDqbz9wrbunHfG2laBERCQW6RKfiIjEJCUoERGJSUpQIiISk5SgREQkJilBiYhITEqsuIiIVBUzKyBovm9AAXCbu392lOvMcfdG5cz/IzATaAb0dPf7j2Z7ItVFNSiR6nXA3Qe4e3/gLuCP1bDNoQS/YTuN4PdtIjWCEpRI9DQh6KKhqB+dP4V9DS02s8vD6X8zs4vC4TfMbFI4fL2Z/a68lYfr+xIYDMwm+GHpE2Z2TwT3SaTK6BKfSPWqb2YLgWTgWILn+wFcAgwg6G+pJTDPzGYBs4BTCR4h0y5cBuAUgicYlMnd/8fMXgHGAD8FPnT3k6t2d0QiRzUokepVdImvJ8HjYKaGj4c5BZjm7gXuvg34iKDm8zFwathh5DK+eUjnSUBl7l0NBBYCPcPlRWoM1aBEosTdZ5tZS6AVpXdTgLtvCrt6GE5Qm2oB/ADICTuTLFX4/LxnCZ4mvYOggz0La28nufuBqtwXkUhQDUokSsysJ8FDiHcSJJ/LzSzBzFoR9GI7Nyw6G7gjLPMx8HMqaOzg7gvdfQCwCugFvA98L6y9KTlJjaAalEj1KroHBUGtaZy7F5jZGwSX7RYRdPD2i7CrCwiS0bnunm5mGwhqURW2xgsT3W53LzSznu6uS3xSo+hp5iIiEpN0iU9ERGKSEpSIiMQkJSgREYlJSlAiIhKTlKBERCQmKUGJiEhMUoISEZGY9P8B5QdQyXFoxB8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior2.plot(label='posterior')\n", "decorate_bowls('Posterior after two vanilla cookies')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "%TODO improve this transition...\n", "\n", "Because the likelihood function is a line, the posterior after two cookies is a parabola.\n", "At this point the high-numbered bowls are the most likely because they contain the most vanilla cookies, and the low-numbered bowls have been all but eliminated.\n", "\n", "Now suppose we draw again and get a chocolate cookie.\n", "Here's the update:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.2462686567164179" ] }, "execution_count": 31, "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": 32, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3gVZfbA8e9JJ5QgIdQACZBIb0ZAwYq6qBQLKNhQVHQXu66r+9t1d9V1dS2oC2JFikixx7IWQMSCQOgdQg81JJQQSD+/P+7gXmNCCtzMzc35PM99mPLOO2fmDjl33nlnRlQVY4wxxt8EuR2AMcYYUxJLUMYYY/ySJShjjDF+yRKUMcYYv2QJyhhjjF+yBGWMMcYvWYIy1YKIHBGR1lW0rj4istFZ5xVVsU63eO9XEZkoIk86w+eLSJq70f2aiMwVkduqcH1xIqIiElJV6zS/ZgnK/IqIbBWRY84frr0i8raI1DmJ+k7Jf3JVraOqm0+mjgp4HBjrrPNjX/xhFJHnnCSYJSLrROSmU1l/eVXFfhWRa0TkJxE5KiJzfbkut/hjQg8ElqBMSQaqah2gB3Am8Be3AjnZxFbJ5VsBq09mvcViCC5hcjYwEIgCRgAvicjZp2qdfiYTeBF42u1ATPViCcqUSlV3Av8FOgGISDMRSRaRTBFJFZHbj5cVkZ4ikiIih50zrxecWfOcfw86Z2VnOeVHishaETkgIl+JSCuvulRERovIRmCj17S2znCUiEwWkXQR2SYifxGRIGfezSLyo4iMEZFM4O/Ft8uJdb6IHBSR3SIyVkTCnHmbgNbAp068/wLOAcY642Odcu1E5BtnX6wXkWu86p8oIuNF5AsRyQYuKGHf/k1V16lqkaouAL4Hzirpe3D20wCv8RAR2S8iPZzx90Rkj4gcEpF5ItKxWCzjRORz52xtgYi0Kbav25a03mIxPCIim5w61ojIlWUt47Wts1R1JrCrPOVFZLCILHOOpU0i0t9rdivn+80Ska9FpKHXcoNEZLXzvc4VkfZe81qIyIfOMZPh9T0GOcfPNhHZ5xxXUaXEdYvzXWSJyGYRucOZXhvP/5NmzjFyxPm/EuS13zJEZKaINCjvfjOAqtrHPr98gK3ARc5wCzxnEk84498BrwARQDcgHejnzJsP3OgM1wF6O8NxgAIhXuu4AkgF2gMheM7QfvKar8A3QAOglte0ts7wZOAToK5T/wbgVmfezUABcLdTd60StvEMoLczPw5YC9xX0j5wxucCt3mN1wZ2ALc4dfQA9gMdnfkTgUNAHzw/AiPK2Oe1gN1A/1LmPwZM9Rq/HFjnNT7S2RfheM5UlnnNm4jnDKanE+tUYHqxfd3Wq+yTzvD5QJpXuaFAM2d7rsVzBti0gsfWbcDcMsr0dPbdxc66mgPtvL6HTUCis8/mAk878xKdmC4GQoGHnWMsDAgGlgNjnO8uAujrte9S8fwoqQN8CEwp6dh19nsbQIDzgKNAj5L2lzPtPuBnINb5bl4Dprn9f7w6fVwPwD7+9XH+OB8BDgLb8CSkWniSVSFQ16vsv4CJzvA84B9Aw2L1/eo/uTPtvzgJxRkPcv6zt3LGFbiwWD0KtHX+2OQCHbzm3XH8Dx+eBLW9gtt8H/BRsX1wogR1LfB9sTpeA/7mDE8EJldg/ZOALwEpZX5bIAuIdManAo+VUra+s6+ivGJ502v+Zfw6uZUrQZWwnmXA4Aru5/IkqNeAMaXMmwv8xWv8D8CXzvBfgZnFjqmdznachefHVEgJdc4G/uA1fjqQz/9+vGhJyzllPwbuLW1/4fnh089rvOnxuiuy32ryx5r4TEmuUNX6qtpKVf+gqsfw/HrOVNUsr3Lb8PzCBbgVz6/YdSKyyLtJqgSt8FxzOSgiB/H8whevusBzhlKShnh+FW8rJY4TLQuAiCSKyGdOs9hh4Cmn3vJqBfQ6Hr+zDdcDTcobg1csz+JpQr1Gnb9ixalqKp4/dgNFJBIYBLzrLB8sIk87zUiH8SRXim3PHq/ho3jOFCpERG5ymt2Ob28nKrbPyqsFnrOk0pS2Lc3wOiZUtQjPd9DcqXObqhaUUN+vlnOGQ4DGxQuKyKUi8rPTrHsQT7I/0T5oBXzktc/W4vmR95u6Tcms+6Qpr11AAxGp65WkWuL5lYqqbgSGO9eCrgLeF5FoPL9Ai9sB/FNVp55gfaU9Zn8/nl+hrYA1xeMoY9njxgNLgeGqmiUi9wFDKhDLDuA7Vb24Asv8hoj8A7gUOE9VD5dRfBowHM+ZwRonaQFcBwwGLsKTnKKAA3gS/ikhnuuDbwD9gPmqWigiy07lOrzswNOMVlG7gM7HR0RE8CSmnXjOuFuKSEgJSWoXnmPpuJZ4moj34mmaO15fOPABcBPwiarmi8jH/G8flHacj1TVHyuxPQbrJGHKSVV3AD8B/xKRCBHpguesaSqAiNwgIjHOL9eDzmKFeJpWivC08R/3KvDo8Yv54un0MLSccRQCM4F/ikhd54/nA8A7FdicusBh4IiItAN+X0b5vcXi/wxIFJEbRSTU+ZzpfVG+LCLyKJ7kcrGqZpRjkenAJU6s7xbbllwgA4jEczZ4qtXG8wc4HTydBXA6zjjjx28liCtpYecsLwLPD+Ig5/gJLWVdbwG3iEg/p5NBc+c7KstM4HJnuVDgQTz75SdgIZ5rfE+LSG1n/X2c5aYB94tIvHhup3gKmFFCIgvDcx0pHSgQkUvxfB/H7QWii3WweBXPcdrK2Q8xIjK4HNtiHJagTEUMx9Muvwv4CM81l2+cef2B1SJyBHgJGKaqOap6FPgn8KPT1NFbVT8CngGmO81Sq/CcSZTX3XguiG8GfsDzB3tCBZZ/CE9yyMJzZjCjjPIvAUPE0+PwZecM8hJgGJ59sQfP9oRXIIan8PxaP35D8BER+XNphVV1N56OKGcXi3cynmapnXjOKH+uQAzloqprgOed9e/Fc6bifVbQwiuGktwIHMNz5nqOM/xGKetaiKfzyRg8nSW+49dnOKXFuB64AfgPnrPsgXhul8hzftQMxHMtbzuQhuc6IniOmyl4rqFuAXLwHF/F688C7sGTCA/gOX6Sveavw5PsNjvHeTM8x00y8LWIZOH5bnqVtS3mf6SUZm9jjCkXEfkLkK6qr7kdiwkslqCMMcb4JWviM8YY45csQRljjPFLlqCMMcb4pRp9H1TDhg01Li7O7TCMMaZGW7x48X5VjSk+vUYnqLi4OFJSUtwOwxhjajQR2VbSdGviM8YY45csQRljjPFLlqCMMcb4pRp9DcoYYyoiPz+ftLQ0cnJy3A6lWoqIiCA2NpbQ0NIexfhrlqCMMaac0tLSqFu3LnFxcXgemG7KS1XJyMggLS2N+Pj4ci1jTXzGGFNOOTk5REdHW3KqBBEhOjq6QmeflqCMMaYCLDlVXkX3nTXxGWNck19YxP4juaRn5bLvcC5H8wspLCqisAhCgoSoyFDq1wqlYZ1wmtWvRXCQJYeaxKcJSkT643knSjDwpqo+XWx+OJ732ZyB54Vr16rqVudNrO8DZwITVfUup3xd4HuvKmKBd1T1PhG5GXiW/72TZqyqvumzjTPGVEhOfiHLdhwkZWsma/dksWFPFpv3Z1NYVL43KoSFBBEfXZu2jerQvWV9ereOpn3Tepa0KuHjjz8mMTGRDh06VGi55ORk1qxZwyOPPOKjyH7NZwlKRIKBccDFeF4QtkhEkp2Xnx13K3BAVduKyDA8L327Fs9Lw/6K562dv7y503lpWDevdSwGPvSqb8bxZGaMcd+OzKN8uWoPc9btY8n2A+QWFAHQokEtTm9cj4s7NCb2tEhi6oYTUzecOuHBBAcFESxCXmERh47lc+hYHulZuWxOz2ZT+hGWpx3k85W7AagbEcL5pzdiQJemnJcYQ0RosJubW218/PHHDBgwoEIJqqCggEGDBjFo0KAKLRMSUvk048szqJ5AqqpuBhCR6cBgPG/9PG4w8Hdn+H1grIiIqmYDP4hI29IqF5EEoBG/PqMyxrgs40guHy7ZSfLyXazceQiAdk3qcn2vVpzVJpqecQ2IiixfN+PS7D50jAWbM/lp035mrd3Hp8t3USc8hMs6N+Hms+Pp0KzeqdgUv7R161b69+9Pr169WLp0KYmJiUyePJn58+fz0EMPUVBQwJlnnsn48eMJDw/nkUceITk5mZCQEC655BKuuuoqkpOT+e6773jyySf54IMPABg9ejTp6elERkbyxhtv0K5dO26++WYaNGjA0qVL6dGjB507dyYlJYWxY8eybds2Ro4cSXp6OjExMbz99tu0bNnyN8s8//zzld5WXyao5sAOr/E0fvu641/KqGqBiBwCovG8srksw/GcMXm3D1wtIucCG4D7VXVH8YVEZBQwCqBly5bl3BRjzImoKj9vzmTqgm18tXoP+YVK1xb1efTSdlzaqSktoyNP6fqaRtXiiu7NuaJ7c/ILi5i/KYPPVuzi0+W7mZmSRq/4BtzaN56LOzT2WaeGf3y6mjW7Dp/SOjs0q8ffBnYss9z69et566236NOnDyNHjuSFF17gtddeY/bs2SQmJnLTTTcxfvx4brrpJj766CPWrVuHiHDw4EHq16/PoEGDGDBgAEOGDAGgX79+vPrqqyQkJLBgwQL+8Ic/MGfOHAA2bNjArFmzCA4OZuLEib/EcNddd3HTTTcxYsQIJkyYwD333MPHH3/8m2VOhi8TVElHRfHG5vKUKc0w4Eav8U+BaaqaKyJ3ApOAC39TuerrwOsASUlJ9jphY05CYZHy9eo9jP9uEyvSDhFVK5QberdieM+WJDauWyUxhAYHcW5iDOcmxvB/l3VgRsp2Jv20jVFTFtM1NoqH+7ejT9uGVRJLVWnRogV9+vQB4IYbbuCJJ54gPj6exMREAEaMGMG4ceO46667iIiI4LbbbuPyyy9nwIABv6nryJEj/PTTTwwdOvSXabm5ub8MDx06tMREM3/+fD780HOF5cYbb+Thhx8uc5mK8mWCSgNaeI3HArtKKZMmIiFAFJBZVsUi0hUIUdXFx6epaoZXkTfwXM8yxvhAUZHy+crdjJm1gc3p2cRFR/LUlZ25qkdzV68DRUWGMurcNozsE89HS3fy4qyNXP/mAvq2bcjfBnYg4RQmzfKc6fhKec8KQ0JCWLhwIbNnz2b69OmMHTv2lzOj44qKiqhfvz7Lli0rsY7atWtXOKbyLlMWX94HtQhIEJF4EQnDc8aTXKxMMjDCGR4CzCnWZFea4cA07wki0tRrdBCwtlJRG2NKparM25DOoHE/cPe0pYQFBzHuuh7MfvB8ruvV0m86KYQEBzE0qQWzHzyPvw7owMqdh7js5e954ev15OQXuh3eSdu+fTvz588HYNq0aVx00UVs3bqV1NRUAKZMmcJ5553HkSNHOHToEJdddhkvvvjiL0mobt26ZGVlAVCvXj3i4+N57733AM93vHz58jJjOPvss5k+fToAU6dOpW/fvqd8O32WoFS1ALgL+ApPspipqqtF5HEROd4N5C0gWkRSgQeAX/ouishW4AXgZhFJExHv7ibXUCxBAfeIyGoRWQ7cA9zsg80ypsbauj+bkRMXcdOEhRw8ms+Ya7vy+T3ncHmXpn7b1TsiNJhb+8Yz+8HzuLxzU16ek8plL33P4m1lNtT4tfbt2zNp0iS6dOlCZmYm999/P2+//TZDhw6lc+fOBAUFceedd5KVlcWAAQPo0qUL5513HmPGjAFg2LBhPPvss3Tv3p1NmzYxdepU3nrrLbp27UrHjh355JNPyozh5Zdf5u2336ZLly5MmTKFl1566ZRvp5TvhCUwJSUlqb2w0JgTO5ZXyPi5qbz63WbCQoK4t18CN53divAQ/zhbqojvN6bz6Icr2X0ohwcuTuTO89pUKLmuXbuW9u3b+zDCsm3dupUBAwawatUqV+OorJL2oYgsVtWk4mXtSRLGmFL9vDmDP32wgm0ZR7miWzP+fFl7GtWLcDusSjsnIYYv7j2HP3+4kme/Ws+PqfsZc203GlfjbQpk9iw+Y8xvZOcW8Ngnqxj2+s+owru39+LFYd2rdXI6rl5EKP8Z3p1nru7Mku0HGPifH1iRdtDtsMotLi6u2p49VZQlKGPMryzZfoBLX/qeKT9v45Y+cXx53zmc3SawummLCNee2ZKPR/chNDiIoa/O57MVxTsZl6wmXxY5WRXdd5agjDGA556msXM2MvTV+RQWKTNGncXfBnYkMixwrwS0a1KPT+7qQ6fmUdz17lL+M3vjCf+IRkREkJGRYUmqEo6/Dyoiovxn4YF75Bljym3f4Rzumb6UnzdnMrBrM568ohNRtU7ucUTVRcM64bx7ey8e/WAlz3+zgcyjefz18g4EldB5IjY2lrS0NNLT012ItPo7/kbd8rIEZUwNt3BLJqPfXcKRnAKeHdKFIWfE1rh3HoWHBPP8NV2pHxnGhB+3cCSngKev7vKbHn6hoaHlfhusOXmWoIypoVSVt37Ywr/+u45WDSJ559ZenN6kah5P5I9EhL8OaE/diBBemr2R7LwCXhrWndBguxLiFktQxtRAOfmFPPLBCj5etov+HZvw7NAu1I2oGU16JyIi3H9xInUjQnjy87WEBC1nzLXd/PZG5EBnCcqYGmZfVg6jJi9m2Y6DPHRJIqMvaFvjmvTKcts5rSkoUp7+7zoiw4L511WdbR+5wBKUMTXI6l2HuG1SCgeP5vPqDWfQv1MTt0PyW3ee14bs3AL+MyeVyLAQ/jqgvSWpKmYJypgaYt6GdH7/zmLq1QrlvTvPolPzKLdD8nsPXJzIkdwCJvy4heg6YYy+oNR3qBofsARlTA3wXsoOHv1wJW0b1WHSyJ72aJ9yEhEeG9CBzOw8nv1qPa2iIxnQpZnbYdUYlqCMCWCqytg5qTz/zQb6tm3I+Bt6WGeIChIRnrm6C7sOHuOBmctpGlWLM1qd5nZYNYL1nzQmQBUVKU98tpbnv9nAld2bM+HmMy05VVJEaDCv3ZhEs6gIbp+cwraMbLdDqhEsQRkTgAoKi3j4gxVM+HELN58dx/NDuxIWYv/dT0aD2mG8fUtPilS5Y8pijuYVuB1SwLMj1pgAk1dQxF3vLuX9xWncf1EifxtY8mN7TMXFN6zNy8O6s35vFo9+uNKeyedjlqCMCSC5BYX8/p3FfLl6D38b2IF7L0qwrtGn2LmJMTx4cSKfLNvF5Pnb3A4noFmCMiZA5OQXMmryYmav28dTV3bmlj72zDhf+cP5bbmofSOe+GxNtX99vD+zBGVMADiWV8htk1KYtzGdf1/dhet6tXQ7pIAWFCQ8f003Yk+rxeipSzmQned2SAHJpwlKRPqLyHoRSRWRR0qYHy4iM5z5C0QkzpkeLSLfisgRERlbbJm5Tp3LnE+jE9VlTKDLyS9k1JQUfty0n+eGdOWaM1u4HVKNEFUrlLHX9SAjO9euR/mIzxKUiAQD44BLgQ7AcBHpUKzYrcABVW0LjAGecabnAH8FHiql+utVtZvz2VdGXcYErOPXnH5I3c+zQ7py9Rnlf9eOOXmdmkfx0CWn8+XqPcxYtMPtcAKOL8+gegKpqrpZVfOA6cDgYmUGA5Oc4feBfiIiqpqtqj/gSVTlVWJdlQ/fGP+WV1DE6KlL+XZ9Ov+6sjNDLDm54vZzWtOnbTT/+HQNm9KPuB1OQPFlgmoOeP+kSHOmlVhGVQuAQ0B0Oep+22ne+6tXEqpsXcZUO4VFyv0zlzFr7V6euKITw3raNSe3BAUJL1zTjYjQIO6dvpT8wiK3QwoYvkxQJZ29FG+kLU+Z4q5X1c7AOc7nxorUJSKjRCRFRFLstc2mOlJV/u+jlXy+Yjf/d1l7buzdyu2QarzG9SL411VdWLXzMOPnbnI7nIDhywSVBnhfrY0FdpVWRkRCgCjghH02VXWn828W8C6epsRy16Wqr6tqkqomxcTEVHCTjHGXqvLPz9cyfdEO7rmwLbef29rtkIyjf6cmDOzajP/M2ci6PYfdDicg+DJBLQISRCReRMKAYUBysTLJwAhneAgwR0/QFUZEQkSkoTMcCgwAVlWmLmOqo3HfpvLmD57HF91/caLb4Zhi/jGoI/UiQvnjeysosKa+k+azBOVcB7oL+ApYC8xU1dUi8riIDHKKvQVEi0gq8ADwS1d0EdkKvADcLCJpTg/AcOArEVkBLAN2Am+UVZcxgWDawu089/UGrurenMcGdLAnRPihBrXDeHxwJ1buPMTr3292O5xqT2rySUZSUpKmpKS4HYYxZfpq9R5+/85izk2M4Y2bkggNtnvs/dnv31nM7LX7+O9959Ampo7b4fg9EVmsqknFp9tRboyfW7A5g7unLaVLbH1eub6HJadq4PHBnYgIDeKxT1bZDbwnwY50Y/zYxr1Z3D45hdjTajHh5jOJDLN3jFYHMXXD+WP/dvyYmkHy8uJ9w0x5WYIyxk/tO5zDzW8vIjw0mEm39KRB7TC3QzIVcF3PlnSJjeLJz9dyOCff7XCqJUtQxvih7NwCbpm4iANH85gw4kxaNIh0OyRTQcFBwpNXdGL/kVxe+HqD2+FUS5agjPEzBYVF3PXuEtbtyWLcdT3oHBvldkimkrrE1ueGXq2YPH8rq3YecjucascSlDF+RFX5x6dr+HZ9Ok8M7sQF7Rq5HZI5SQ/97nROiwzj8U/XWIeJCrIEZYwfmfjTVqb8vI07zm1t73QKEFG1QnngkkQWbs3ky1V73A6nWrEEZYyfmL12L098toZLOjTmT/3buR2OOYWuTWpBuyZ1eeq/a8nJL3Q7nGrDEpQxfmDNrsPcPW0pHZtF8eKwbgQF2VMiAklIcBB/ubwDOzKPMfGnrW6HU21YgjLGZelZudw+OYV6EaG8OSLJ7nUKUH0TGtKvXSPGzkklPSvX7XCqBUtQxrgot6CQO99ZTEZ2Lm/clETjehFuh2R86M+Xtycnv5Axs6zbeXlYgjLGJZ73Oq1i8bYDPD+0m3UnrwHaxNTh+l4tmbFoB1v2Z7sdjt+zBGWMS976YQvvL07j3n4JXN6lqdvhmCpy14UJhIcE8cI3dhZVFktQxrhg3oZ0nvpiLf07NuHefgluh2OqUEzdcEb2iefT5bvs5t0yWIIypopty8jm7mlLSWhUl+ev6Wo99mqgUee1pn5kKM99vd7tUPyaJShjqlB2bgG3T/a8g+z1m86gdrj12KuJ6kWE8vvz2jB3fTo/b85wOxy/ZQnKmCqiqjw4czmp+44w7roetIqu7XZIxkUjzo6jcb1wnv1qvT0CqRSWoIypIq/M3cSXq/fw6KXt6ZvQ0O1wjMsiQoO5+8IEFm87wA+p+90Oxy9ZgjKmCszbkM5zX69nYNdm3HZOvNvhGD8xNCmWplERvDRro51FlcCnCUpE+ovIehFJFZFHSpgfLiIznPkLRCTOmR4tIt+KyBERGetVPlJEPheRdSKyWkSe9pp3s4iki8gy53ObL7fNmPLakXmUu6ct5fTGdXnm6s6IWKcI4xEeEszvz29DyrYD/LTJrkUV57MEJSLBwDjgUqADMFxEOhQrditwQFXbAmOAZ5zpOcBfgYdKqPo5VW0HdAf6iMilXvNmqGo35/PmKdwcYyrlWF4hd0xZjKry2o1n2GOMzG9ck9SCxvXC7SyqBL48g+oJpKrqZlXNA6YDg4uVGQxMcobfB/qJiKhqtqr+gCdR/UJVj6rqt85wHrAEiPXhNhhTaarKXz5exdo9h3lpWHfrFGFKFBEazO/Pa8PCrZn8vDnT7XD8ii8TVHNgh9d4mjOtxDKqWgAcAqLLU7mI1AcGArO9Jl8tIitE5H0RaVHZwI05Fd5duJ0PlqRx94UJ9uJBc0LDerakUd1wXpptT5fw5ssEVVJDe/Hz1/KU+W3FIiHANOBlVd3sTP4UiFPVLsAs/ndmVnzZUSKSIiIp6enpZa3KmEpZvuMg/0hew7mJMfakCFOmiNBg7jyvDT9vziRlq51FHefLBJUGeJ/FxAK7SivjJJ0ooDzfzuvARlV98fgEVc1Q1ePPsH8DOKOkBVX1dVVNUtWkmJiYcm2IMRVxIDuPP0xdQkzdcF66thvB9qQIUw7DerbgtMhQXv1uk9uh+A1fJqhFQIKIxItIGDAMSC5WJhkY4QwPAeZoGVcJReRJPInsvmLTvZ+2OQhYexKxG1MpRUXKfTOWkZ6Vy/gbenBa7TC3QzLVRGRYCCPOjmPW2n1s2Jvldjh+wWcJyrmmdBfwFZ5kMVNVV4vI4yIyyCn2FhAtIqnAA8AvXdFFZCvwAnCziKSJSAcRiQX+D0+vwCXFupPf43Q9Xw7cA9zsq20zpjTjvk3luw3pPDawA11i67sdjqlmRpwVR63QYF77bnPZhWsAqcndGpOSkjQlJcXtMEyA+DF1Pze+tYBBXZsx5tpudr+TqZS/J6/mnZ+3Me/hC2hWv5bb4VQJEVmsqknFp9uTJIw5BfYezuHe6UtpHVOHf15pN+OayrvtnHgUz/vCajpLUMacpILCIu5+dynZuYWMv76HPaHcnJTY0yIZ1LUZ0xZu5+DRPLfDcZUlKGNO0gvfbGDh1kyeuqoTCY3ruh2OCQB3nNeao3mFTF2w3e1QXGUJypiTMHf9Pl6Zu4lhZ7bgyu72UBNzarRrUo++bRsyef5W8guL3A7HNZagjKmk3YeOcf+MZbRrUpe/D+rodjgmwNzaN569h3P5YuVut0NxjSUoYyqhoLCIe6YtJa+giFeu70FEaLDbIZkAc15iDK0b1mbCD1tq7ENkLUEZUwljZm1g0dYDPHVVZ1rH1HE7HBOAgoKEW/rEsTztEEu2H3Q7HFdYgjKmguZtSP/lutPgbsWff2zMqXNVj1jqRYQw4cea2eXcEpQxFbDvcA73z1hGQqM6/G2gXXcyvlU7PIThPVvy5ao97Dx4zO1wqpwlKGPKqbBIuXf6Mo7mFTLuuh7UCrPrTsb3bjo7DoDJ87e6GYYrLEEZU06vfJvK/M0Z/GNwR7vfyVSZ5vVrcXH7xryXkkZOfqHb4VQpS1DGlMPCLZmMmbWBK7o1Y+gZdr+TqVo3ntWKzOw8/ruqZnU5twRlTBkOZOdx7/SltGwQyZP2nD3jgrPbRNO6YW2mzN/mdj+iapIAACAASURBVChVyhKUMSegqvzx/RXsP5LL2Ot6UMees2dcICJc37sVS7YfZPWuQ26HU2UsQRlzApN+2sqstXt59NL2dGoe5XY4pgYb0iOWiNAg3vm55jyfzxKUMaVYvesQT32xjn7tGnFLnzi3wzE1XFRkKIO6NuPjpTs5nJPvdjhVwhKUMSU4mlfA3dOWclrtUJ4d2tWuOxm/cGPvOI7lF/Lh4jS3Q6kSlqCMKcHfPlnNlv3ZjLm2Gw1qh7kdjjEAdI6NomtsFO8s2F4jns9nCcqYYpKX7+K9xWmMPr8tZ7dp6HY4xvzK8J4tSd13pEY8n8+nCUpE+ovIehFJFZFHSpgfLiIznPkLRCTOmR4tIt+KyBERGVtsmTNEZKWzzMvitL2ISAMR+UZENjr/nubLbTOBaUfmUf7vw5X0aFmf+y5KcDscY35jQNdmRIYFM2NR4HeW8FmCEpFgYBxwKdABGC4iHYoVuxU4oKptgTHAM870HOCvwEMlVD0eGAUkOJ/+zvRHgNmqmgDMdsaNKbf8wiLumb4UgJeGdSck2BoYjP+pEx7CwC7N+HT5brICvLOEL/8H9gRSVXWzquYB04HBxcoMBiY5w+8D/UREVDVbVX/Ak6h+ISJNgXqqOl89DbCTgStKqGuS13RjyuWlWRtZuv0gT13VmRYNIt0Ox5hSXduzBcfyC/lsRWA/WcKXCao5sMNrPM2ZVmIZVS0ADgHRZdTp3X3Fu87GqrrbqWs30KjSkZsaZ/6mDMbNTeWapFgGdm3mdjjGnFD3FvVJbFyH6Yt2lF24GvNlgiqpX27xbiflKXMy5X9bgcgoEUkRkZT09PSKLGoC1IHsPO6fsYz46Nr26nZTLYgI157ZkuU7DrJ292G3w/EZXyaoNKCF13gssKu0MiISAkQBmWXU6f2kTu869zpNgMebAveVVIGqvq6qSaqaFBMTU85NMYFKVXn4gxVkZOfy8vDuRIbZo4xM9XBl9+aEBQcxI4DPonyZoBYBCSISLyJhwDAguViZZGCEMzwEmKMn6NzvNN1liUhvp/feTcAnJdQ1wmu6MaWaumA736zZy5/6t7NHGZlqpUHtMC7p2JiPlu4ktyAwX8PhswTlXFO6C/gKWAvMVNXVIvK4iAxyir0FRItIKvAAXj3vRGQr8AJws4ikefUA/D3wJpAKbAL+60x/GrhYRDYCFzvjxpRqw94snvhsDeckNGRkn3i3wzGmwoYmteDQsXzmrC2xwajak5pwN3JpkpKSNCUlxe0wjAty8gu5YtyP7D+Syxf3nkOjuhFuh2RMhRUWKWc/PZvOzaN4c8SZbodTaSKyWFWTik8/4RmUiEz0Gh5xgqLGVCtP/3cd6/Zk8eyQrpacTLUVHCRc0b05c9ens/9IrtvhnHJlNfF19Rq+15eBGFNV5qzby8SftnJLnzguaGd3I5jq7eoesRQUKcnLivdBq/7KSlA1t/3PBKR9WTn88b0VtGtSlz/1b+d2OMactMTGdencPIoPlgTeE87L6lMbKyIv47n/6PjwL1T1Hp9FZswpVlSkPPTeCo7kFjB9VG8iQoPdDsmYU+LqHs35+6drWLfnMO2a1HM7nFOmrDOoPwKLgRSvYe+PMdXGhB+3MG9DOn8Z0IGExnXdDseYU2Zg12aEBAkfLtnpdiin1AnPoFR10onmG1NdrNp5iGe+XMfFHRpzQ6+WbodjzCkVXSecC9o14qOlO3n4d6cHzIOOT5igRKT4jbW/oqqDTjTfGH9wNK+Ae6cvpUHtMJ65uou9HdcEpKt7xPLNmr38kLqf808PjM4/ZV2DOgvPw1ynAQso+Vl4xvi1Jz5by+b92bxzay97O64JWBe0i6FuRAjJy3cFTIIq6zywCfBnoBPwEp4nNOxX1e9U9TtfB2fMyfpy1R6mLdzOqHNb06etvR3XBK7wkGAu69SUr1btISc/MB59dMIEpaqFqvqlqo4AeuN5vNBcEbm7SqIz5iTsPnSMRz5cQefmUTx48eluh2OMzw3u1ozsvEJmB8ijj8q8kua8lv0q4B1gNPAy8KGvAzPmZBQWKQ/MWE5eQREvDetGWEhgXDQ25kR6tY6mUd1wPlkWGL35yuokMQlP895/gX+o6qoqicqYk/TavE3M35zBv4d0oXVMHbfDMaZKBAcJA7s2Y8r8bRw6mk9UZKjbIZ2Usn5W3ggk4nnM0XwROex8skQkcN+SZaq1ZTsO8sLXG7i8S1OGnhFb9gLGBJDB3ZqRV1jEl6ur/+vgy7oGFaSqdb0+9ZxPXVUNnNuVTcA4klvAPdOW0rheBE9d0dm6lJsap3PzKOIb1iZ5efV/Nl9ZTzOPEJH7RGSs86p0e92o8WuPfbKKtANHeXFYt2rfvGFMZYh4mvl+2pTBvsM5bodzUspq4psEJAErgcuA530ekTGV9MmynXy4ZCd3X5jAmXEN3A7HGNcM7tYMVfh0RfVu5isrQXVQ1RtU9TU8r2Q/pwpiMqbCtmcc5f8+WsUZrU7j7gvbuh2OMa5qE1OH9k3r8fmK6t3MV1aCyj8+4LzC3Ri/k19YxL0zliLAi9d2C5jnkBlzMgZ0acqS7QfZdfCY26FUWpkvLPTuuQd0sV58xt+8NGsjS7cf5KmrOtOiQaTb4RjjFy7r3BSAL1ZW32a+snrxBRfruRdSkV58ItJfRNaLSKqIPFLC/HARmeHMXyAicV7zHnWmrxeR3znTTheRZV6fwyJynzPv7yKy02veZRXdGab6+WnTfsbNTeWapFgGdm3mdjjG+I34hrXp0LRe4CaokyEiwcA44FKgAzBcRDoUK3YrcEBV2wJjgGecZTsAw4COQH/gFREJVtX1qtpNVbsBZwBHgY+86htzfL6qfuGrbTP+4UB2Hg/MWE58dG3+Pqij2+EY43cur+bNfL5srO8JpKrqZlXNA6YDg4uVGYynpyDA+0A/8dy4MhiYrqq5qroFzzMAexZbth+wSVW3+WwLjN9SVR7+YAWZ2Xm8PLw7kWF2B4QxxVX3Zj5fJqjmeF7VcVyaM63EMk4njENAdDmXHYbnNSDe7hKRFSIyQUROO7nwjT975+dtfLNmLw/3P51OzaPcDscYv1Tdm/l8maBKuoVfy1nmhMuKSBgwCHjPa/54oA3QDdhNKfdsOTccp4hISnp6eunRG7+1bs9hnvh8LeefHsPIPvFuh2OMX6vOzXy+TFBpQAuv8VigeKf8X8o4T6mIAjLLseylwBJV3Xt8gqrudV4PUgS8wW+bBI+Xe11Vk1Q1KSYmplIbZtxzLK+Qu99dSlStUJ4b2pWgIHuUkTEncnk1bubzZYJaBCSISLxzxjMMKP4K+WRghDM8BJijqupMH+b08osHEoCFXssNp1jznog09Rq9ErAnrwegxz9bQ2r6EcZc042GdcLdDscYvxdXjZv5fJagnGtKdwFfAWuBmaq6WkQeF5FBTrG3gGgRSQUeAB5xll0NzATWAF8Co1W1EEBEIvG82bf4O6n+LSIrRWQFcAFwv6+2zbjj8xW7f3k7bt8EezuuMeXVv1MTlmw/WO2ezSeeE5aaKSkpSVNSUtwOw5TDjsyjXPby97SJqcN7d55FqD0twphyW78ni9+9OI8nr+jEDb1buR3Ob4jIYlVNKj7d/pcbv5dfWMTd05aCwn+Gd7fkZEwFJTauQ3zD2ny1eo/boVSI/U83fu+FbzawbMdBnr66iz3KyJhKEBEu6diY+ZsyOHQ0v+wF/IQlKOPX5m1IZ/zcTQzv2ZLLuzQtewFjTIn6d2xCQZEyZ/3esgv7CUtQxm/tO5zDAzOXkdi4Do8NKP6ULGNMRXSNrU/jeuF8uar6NPNZgjJ+qbBIuW/GMo7kFjDuuh7UCgt2OyRjqrWgIOGSDk34bkM6x/IK3Q6nXCxBGb/0yrep/LQpg8cHdSKhcV23wzEmIPTv1ISc/CK+21A9nqJjCcr4nQWbMxgzawNXdGvG0KRYt8MxJmD0jG9AVK1Qvq4mvfksQRm/knEkl3umL6VVdG2evLIznofbG2NOhdDgIPq1b8SstXvJLyxyO5wyWYIyfqOoSLl/5nIOHM1n7HXdqRNur9Aw5lS7pENjDucUkLL1gNuhlMkSlPEb47/bxLwN6fxtYAc6NrNXaBjjC+ckxBAWHMSstf7f3dwSlPELC7dk8vzX6xnYtRnX9WzpdjjGBKza4SGc1SaaWWv34u+PurMEZVy3/0gud09bQssGkTx1ZSe77mSMj13UoTHbMo6yKf2I26GckCUo46rCIuW+6cs4cDSfV64/g7oRoW6HZEzAu6h9IwC+WbPP5UhOzBKUcdXLszfyQ+p+Hh/UkQ7N6rkdjjE1QtOoWnRqXo/Zfn4dyhKUcc33G9N5ec5GrurRnGvPbFH2AsaYU6Zfu8Ys3n6AjCO5bodSKktQxhW7Dx3jvunLSGhUhyevsOtOxlS1izs0RhXmrPPfZj5LUKbK5RUUMXrqEnLyC3nl+jOIDLP7nYypah2b1aNJvQhmr7UEZcwvnvpiLUu2H+SZIV1o26iO2+EYUyOJCP3aN2LexnRy8v3z4bGWoEyVSl6+i4k/bWVkn3gGdGnmdjjG1GgXtW/M0bxCFm7JdDuUElmCMlVm494sHvlgBUmtTuPRy9q5HY4xNd5ZbaIJDwny2+tQPk1QItJfRNaLSKqIPFLC/HARmeHMXyAicV7zHnWmrxeR33lN3yoiK0VkmYikeE1vICLfiMhG59/TfLltpmIO5+Rzx5TFRIYFM/a6HoQG228jY9wWERrM2W2imbu+hiUoEQkGxgGXAh2A4SJS/LWotwIHVLUtMAZ4xlm2AzAM6Aj0B15x6jvuAlXtpqpJXtMeAWaragIw2xk3fqCoSHlo5nK2ZR5l7HU9aBIV4XZIxhjHBe0asTXjKFv2Z7sdym/48mdsTyBVVTerah4wHRhcrMxgYJIz/D7QTzz9jQcD01U1V1W3AKlOfSfiXdck4IpTsA3mFBj/3Sa+XrOXP1/Wnt6to90Oxxjj5fxEz1MlvvXDZj5fJqjmwA6v8TRnWollVLUAOAREl7GsAl+LyGIRGeVVprGq7nbq2g00KikoERklIikikpKeXj3eKlmdzduQznPOQ2BH9olzOxxjTDEtoyNpE1Obb/2wmc+XCaqkOy+LPzq3tDInWraPqvbA03Q4WkTOrUhQqvq6qiapalJMTExFFjUVtC0jm7unLeX0xnV55mp7+aAx/uqC0xuxYHMm2bkFbofyK75MUGmA9/NrYoFdpZURkRAgCsg80bKqevzffcBH/K/pb6+INHXqagr438+BGiQ7t4BRkxcD8PqNSXYzrjF+7IJ2jcgrLOKnTRluh/IrvkxQi4AEEYkXkTA8nR6Si5VJBkY4w0OAOep5QUkyMMzp5RcPJAALRaS2iNQFEJHawCXAqhLqGgF84qPtMmVQVf74/nI27sti7HXdaRkd6XZIxpgTODOuAbXDgv2umc9nP2tVtUBE7gK+AoKBCaq6WkQeB1JUNRl4C5giIql4zpyGOcuuFpGZwBqgABitqoUi0hj4yGkqCgHeVdUvnVU+DcwUkVuB7cBQX22bObFX5m7ii5V7+PNl7TgnwZpRjfF3YSFB9E1oyNx1+1BVv2mO92m7i6p+AXxRbNpjXsM5lJJIVPWfwD+LTdsMdC2lfAbQ7yRDNifpmzV7ee7r9Qzq2ozbz2ntdjjGmHK64PRGfLV6L+v3ZtGuiX+8+sbuljSnzIa9Wdw3fSmdmkXx7yFd/OZXmDGmbOef7un4/N16/+ndbAnKnBIHsvO4bVIKkeEhvHFTEhGhwWUvZIzxG02iIkhsXId5Gy1BmQCSX1jE6HeXsOdwDq/feIY9KcKYaurchBgWbTnA0Tz/6G5uCcqcFFXlsU9W89OmDJ6+qjPdW9ojEI2prs5NjCGvsIgFm/3j6eaWoMxJeeuHLUxbuJ3RF7Thqh6xbodjjDkJPeMbEB4SxHcb/KOZzxKUqbRZa/byzy/WcmmnJjx48eluh2OMOUkRocH0ah3tN9ehLEGZSlm96xD3OD32XrimG0FB1mPPmEBwbkJDNqdnk3bgqNuhWIIyFbfr4DFGTlxEVK1Q3hyRRK0w67FnTKA4L9Fzc/28DftdjsQSlKmgrJx8Rk5cxNHcQt6+5Uwa17Mee8YEkraN6tA0KoLv/aCZzxKUKbf8wiL+MHUJqfuOMP6GM/zmbnNjzKkjIpybEMMPqfspKCxyNRZLUKZcVJVHP1zJ9xv389SVnemb0NDtkIwxPnJuYgxZOQUsTzvoahyWoEy5PPf1et5fnMa9/RK45swWZS9gjKm2+rZtSJDAdy5fh7IEZco0Zf5Wxn27ieE9W3DfRQluh2OM8bGoyFC6xNbnx1RLUMaP/Xflbh5LXs1F7RvxxOBO9gBYY2qIvm0bsmzHQbJy8l2LwRKUKdUPG/dz7/RldG9Rn/8M70FIsB0uxtQUfdo2pLBIXX3skf3FMSVauv0Ao6ak0DqmNm/f3NPudTKmhunRqj4RoUH84GIznyUo8xsb9mZxy8RFNKwTzuSRPYmKDHU7JGNMFQsPCaZnfLSr16EsQZlf2bo/mxveXEBYcBDv3NqLRnYjrjE1Vt+20Wzcd4S9h3NcWb8lKPOLtANHuf7NBRQUKe/c1ouW0ZFuh2SMcVGftp77Hd06i/JpghKR/iKyXkRSReSREuaHi8gMZ/4CEYnzmveoM329iPzOmdZCRL4VkbUislpE7vUq/3cR2Skiy5zPZb7ctkCz51AO172xgKycfCaP7Eli47puh2SMcVn7JvVoUDvMtetQIb6qWESCgXHAxUAasEhEklV1jVexW4EDqtpWRIYBzwDXikgHYBjQEWgGzBKRRKAAeFBVl4hIXWCxiHzjVecYVX3OV9sUqPYdzuG6N38mMzuPKbf2pFPzKLdDMsb4gaAg4aw2nutQqlrlt5n48gyqJ5CqqptVNQ+YDgwuVmYwMMkZfh/oJ549MBiYrqq5qroFSAV6qupuVV0CoKpZwFqguQ+3IeDtO5zDsDd+Zs+hHCbcfKa9EdcY8yt92zZk7+FcNqUfqfJ1+zJBNQd2eI2n8dtk8ksZVS0ADgHR5VnWaQ7sDizwmnyXiKwQkQkiUuJfWhEZJSIpIpKSnu7+03rdtPdwDsNe/5m9h3KYNLInPeMbuB2SMcbP9HWuQ/2wseqb+XyZoEo6F9RyljnhsiJSB/gAuE9VDzuTxwNtgG7AbuD5koJS1ddVNUlVk2JiYk68BQFs96FjnuR02JOczoyz5GSM+a0WDSJp2SCSH1IzqnzdvkxQaYD3U0VjgV2llRGRECAKyDzRsiISiic5TVXVD48XUNW9qlqoqkXAG3iaGE0JtmVkM/TV+ezPymXyrT1JsuRkjDmBPm2jWbAlg8Ki4ucYvuXLBLUISBCReBEJw9PpIblYmWRghDM8BJijqupMH+b08osHEoCFzvWpt4C1qvqCd0Ui0tRr9Epg1SnfogCwcW8WQ1+dT3ZuAe/e3pszWllyMsacWO/W0WTlFLB616EqXa/PevGpaoGI3AV8BQQDE1R1tYg8DqSoajKeZDNFRFLxnDkNc5ZdLSIzgTV4eu6NVtVCEekL3AisFJFlzqr+rKpfAP8WkW54mgK3Anf4atuqqxVpBxkxYSEhwUHMuOMs60pujCmXs1pHAzB/UwZdYutX2XrFc8JSMyUlJWlKSorbYVSJeRvSufOdxZwWGcbU23oR17C22yEZY6qRfs/PpUWDSCbecuqvnojIYlVNKj7dniRRA3y0NI2RExfRKro2H/3hbEtOxpgKO6tNNIu2ZJJfha+BtwQVwFSV8XM3cf+M5STFncaMO3rbs/WMMZVydpuGZOcVsiKt6q5DWYIKUHkFRTzywUqe+XIdl3dpyqSRPakXYU8lN8ZUTm/nOtTPm6uuu7klqAB06Gg+IyYsZEbKDu6+sC3/Gdad8BB7n5MxpvIa1A6jXZO6zN9UdQnKZ734jDtS92UxavJi0g4c44VrunJVj1i3QzLGBIjeraOZvmg7uQWFVfKj186gAsiXq3YzeOyPHM4pYOrtvSw5GWNOqbPaRJOTX8TyHVVzHcoSVAAoKCzi31+u4853lpDQuC6f3d3XHl1kjDnlesdHI0KVNfNZgqrm9hzK4bo3F/DK3E0M79mCGXf0pkmU9dQzxpx6UZGhdGxWj/mbq+bBsXYNqhqbs24vD85cTm5BkV1vMsZUibNaRzNp/jZy8guJCPXtdSg7g6qGjuYV8LdPVjFyYgpNomrx6d19LTkZY6pEr/ho8gqKWL7joM/XZWdQ1czibZk8OHM5WzOOMrJPPA/3P93nv2KMMea4M+MaIAILtmTSy7k3ylcsQVUT2bkFvPDNBt7+cQvN6tdi2u29OauNbw8OY4wpLioylHZN6rFgSwaeF034jiWoauCr1Xv4e/Jqdh/K4fpeLXn0svbUCbevzhjjjl7xDZi+aDt5BUWEhfjuSpFdg/JjqfuyGDlxEXdMWUxUrVA++P3Z/PPKzpacjDGu6hXfgJz8Ilbu9O39UPaXzg/tP5LLi7M2MG3hDiJDg/nzZe24pU88ocH2e8IY476e8Z77LBduyeSMVqf5bD2WoPzI/iO5vPH9Zt6Zv43cgiJu6NWSe/olEF0n3O3QjDHmF9F1wmnbqA4LtmTw+/Pb+Gw9lqD8wI7Mo0z8aStTF2wjr6CIAV2ace9FCbSJqeN2aMYYU6Je8Q34ZNkuCouU4CDxyTosQbmkqEj5cdN+Jv20jdnr9hIkwuBuzRh9QVtLTMYYv9czvgFTF2xnza7DdI6N8sk6LEFVIVVlze7DJC/fxafLdrHrUA7RtcMYfX5bruvVkmb1a7kdojHGlEuveM9tLgu2ZFTPBCUi/YGXgGDgTVV9utj8cGAycAaQAVyrqludeY8CtwKFwD2q+tWJ6hSReGA60ABYAtyoqnm+3L7yyMrJZ/6mDOZtTGfehv1szzxKSJBwbmIMf7q0Hf07NbF3NRljqp0mURG0io5kwZZMbjuntU/W4bMEJSLBwDjgYiANWCQiyaq6xqvYrcABVW0rIsOAZ4BrRaQDMAzoCDQDZolIorNMaXU+A4xR1eki8qpT93hfbZ+3vIIisnLy2XUwh50Hj7Ij8xirdx1ixc5DbNmfjSpEhgVzVutoRp3bmss6N6VB7bCqCM0YY3ymV3wDvl6zl6IiJcgH16F8eQbVE0hV1c0AIjIdGAx4J6jBwN+d4feBsSIizvTpqpoLbBGRVKc+SqpTRNYCFwLXOWUmOfX6LEE9+9U6Js/fxrG8QgqK9Dfzm9SLoHNsFFd0a05S3GkktWrg0xvajDGmqvWMj2ZmShob9mXRrkm9U16/LxNUc2CH13ga0Ku0MqpaICKHgGhn+s/Flm3uDJdUZzRwUFULSij/KyIyChgF0LJly4ptkZdOzaK4ukcskWHBRIYFUzs8hKZRtYg9zfOpH2lnSMaYwHZWm2iu69XSZ/do+jJBlXS+V/xUo7QypU0vaS+cqPxvJ6q+DrwOkJSUVGKZ8ri0c1Mu7dy0sosbY0y117x+LZ66srPP6vdlm1Ma0MJrPBbYVVoZEQkBooDMEyxb2vT9QH2njtLWZYwxphrxZYJaBCSISLyIhOHp9JBcrEwyMMIZHgLMUVV1pg8TkXCnd14CsLC0Op1lvnXqwKnzEx9umzHGGB/zWROfc03pLuArPF3CJ6jqahF5HEhR1WTgLWCK0wkiE0/CwSk3E0+HigJgtKoWApRUp7PKPwHTReRJYKlTtzHGmGpKPCcfNVNSUpKmpKS4HYYxxtRoIrJYVZOKT7d+z8YYY/ySJShjjDF+yRKUMcYYv2QJyhhjjF+q0Z0kRCQd2HYSVTTEcw9WTVbT90FN336wfWDbf/Lb30pVY4pPrNEJ6mSJSEpJPU9qkpq+D2r69oPtA9t+322/NfEZY4zxS5agjDHG+CVLUCfndbcD8AM1fR/U9O0H2we2/T5i16CMMcb4JTuDMsYY45csQRljjPFLlqAqSUT6i8h6EUkVkUfcjsfXRKSFiHwrImtFZLWI3OtMbyAi34jIRuff09yO1ZdEJFhElorIZ854vIgscLZ/hvMamIAlIvVF5H0RWeccC2fVpGNARO53jv9VIjJNRCIC/RgQkQkisk9EVnlNK/E7F4+Xnb+LK0Skx8ms2xJUJYhIMDAOuBToAAwXkQ7uRuVzBcCDqtoe6A2Mdrb5EWC2qiYAs53xQHYvsNZr/BlgjLP9B4BbXYmq6rwEfKmq7YCuePZFjTgGRKQ5cA+QpKqd8LzyZxiBfwxMBPoXm1bad34pnvf3JQCjgPEns2JLUJXTE0hV1c2qmgdMBwa7HJNPqepuVV3iDGfh+cPUHM92T3KKTQKucCdC3xORWOBy4E1nXIALgfedIoG+/fWAc3Hetaaqeap6kBp0DOB5h14t5+3dkcBuAvwYUNV5eN7X562073wwMFk9fsbzpvOmlV23JajKaQ7s8BpPc6bVCCISB3QHFgCNVXU3eJIY0Mi9yHzuReBhoMgZjwYOqmqBMx7ox0FrIB1422nmfFNEalNDjgFV3Qk8B2zHk5gOAYupWcfAcaV956f0b6MlqMqREqbViP76IlIH+AC4T1UPux1PVRGRAcA+VV3sPbmEooF8HIQAPYDxqtodyCZAm/NK4lxnGQzEA82A2niatIoL5GOgLKf0/4QlqMpJA1p4jccCu1yKpcqISCie5DRVVT90Ju89fgrv/LvPrfh8rA8wSES24mnSvRDPGVV9p7kHAv84SAPSVHWBM/4+noRVU46Bi4AtqpquqvnAh8DZ1Kxj4LjSvvNT+rfRElTlLAISnN47YXgulCa7HJNPOddb3gLWquoLXrOSgRHO8Ajgk6qOrSqo6qOqGquqcXi+7zmqej3wLTDEKRaw2w+gqnuAHSJyujOpH7CGGnIM4Gna6y0ikc7/h+Pb///t3T9oZFUUx/HvjyyygiwStngmOAAAAtpJREFUkkIEC6thG2OxuMEEU+kiYiGsaXQl7JZbyBIFt1AEJcK2ol1YbAJaROxsRBMwoEKyBixEhC22CCjLQjBYjD+Le4eI5I8kM85l8/tUw8yb9+Yxbzhz7r3vnBNzDfzDft/5F8CluprvPHCvNxR4FKkkcUSSnqf8gx4BFm2/P+SPNFCSpoBVYJPdOZjrlHmoT4HHKD/gi7b/PaF6X5E0A8zbfkHS45SMahRYB16x/ecwP98gSZqgLBJ5APgVmKP80T0R14Ckd4FZyqrWdeAKZY7lvr0GJC0BM5S2GlvAO8Dn7PGd18D9IWXV3x/AnO0fjnzsBKiIiGhRhvgiIqJJCVAREdGkBKiIiGhSAlRERDQpASoiIpp06vBNImJQJHUpS/cFdIGrtr895j63bT90wOsLwJfAw0DH9gfHOV7EoCSDihiuHdsTtp8A3gIW/odjPkW5f+0Zyr1tEU1KgIpoxxlKu4ZeX50bte/QpqTZ+vxHkl6sj5clLdbHlyW9d9DO6/5+BM4Ba5SbTD+W9PYAzyniyDLEFzFcD0raAE4Dj1Bq/AG8BExQei6NAd9LWgFWgGlKSZlH63sApijVDPZl+w1JnwGvAteAr20/3d/TieifZFARw9Ub4utQysN8UsvFTAFLtru2t4BvKJnPKjBdm0X+xG7Rzkngv8xdPQlsAJ36/ohmJYOKaITtNUljwDh7ty3A9p3a9uECJZsaBV4GtmsjyT3VGno3KdWlf6M021PN3iZt7/TzXCL6IRlURCMkdSjFh3+nBJ9ZSSOSximdbL+rm64Br9dtVoF5DlnsYHvD9gTwM3AW+Ap4rmZvCU7RpGRQEcPVm4OCkjW9ZrsraZkybHeL0vDtzdruAkowetb2L5JuU7KoQ1fj1UB31/Zfkjq2M8QXTUs184iIaFKG+CIiokkJUBER0aQEqIiIaFICVERENCkBKiIimpQAFRERTUqAioiIJv0NPr2F7A3fOpsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior3.plot(label='posterior')\n", "decorate_bowls('Posterior after 2 vanilla, 1 chocolate')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now Bowl 100 has been eliminated because it contains no chocolare 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 posteori probability\", where \"a posteori\" is unnecessary Latin for \"posterior\".\n", "\n", "To compute the MAP, we can use the `Series` method `idxmax`:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "67" ] }, "execution_count": 33, "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": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "67" ] }, "execution_count": 34, "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": { "tags": [ "remove-cell" ] }, "source": [ "Here's a figure for the book." ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAGoCAYAAACdRPr5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hVVfbw8e9KI/TeA4TeCWLoKEWKIoINsY2gKCqiODrOTxx0UOcdu6POCIqgIBZQRh1soIKhSJFgQapSQgi9hkAIpKz3j3OCMaQnN+cmWZ/nuU/uPXefc9c9IYtT9t5LVBVjjPGFAK8DMMaUXpZgjDE+YwnGGOMzlmCMMT5jCcYY4zOWYIwxPmMJxpQZIhIjIgO9jqMssQRjsiUiE0QkWkTOiMisLN6/RES2iEiiiHwrIk0yvHediKx034vK5XP6iUiaiJx0H3tE5PGi/0amuFmCMTnZC/wDeDPzGyJSC/gIeBSoAUQD8zI0OQq8BDyd189S1UqqWgnoA4wVkSsLEbvxA5ZgTLZU9SNV/QQ4ksXbVwMbVfVDVU0CpgARItLGXfcbVf0AJ0nl93N3AiuBdunLRKSXiKwVkXj3Zy93eX8R+SVDu29E5PsMr1dklahEpJt7dHZCRA6IyIv5jdPkzhKMKaj2wM/pL1T1FLDdXV4oItIS6A2sdl/XAD4HXgFqAi8Cn4tITWAV0EJEaolIENABCBORyiJSHrgQWJ7Fx7wMvKyqVYDmwAeFjduczxKMKahKQHymZfFA5QJur4GIHBeRE8CvwBpghfve5cBvqjpHVVNU9X1gC3CFe/QUDVwMRALr3fV6Az3c9bI6AkvGTUyqelJVVxcwbpMDSzCmoE4CVTItqwIkFHB7e1W1mntEUQ04Dcx232sA7MrUfhfQ0H2+FOiHk2SWAlFAX/exNJvPGwu0Ara4p1zDChi3yYElGFNQG4GI9BciUhHnVGNjYTesqvHAe8AV7qK9QJNMzRoDe9znmRPMUnJJMKr6m6reANQBngHmu9/BFCFLMCZbIhIkIqFAIBAoIqHudQ6Aj4EOInKN2+YxYL2qbnHXDXSXBwEB7rrBefzcSsD1/J6svgBaiciNbkyjcC4Af+a+vxJoDXQDvlfVjTgJqTuwLJvPuFlEaqtqGnDcXZyapx1j8k5V7WGPLB84d4Y002NKhvcH4lwLOY1zWhKe4b0xWaw7K5vP6Qek4Zx2ncS5a/U50CJDmz7AOpzrPOuAPpm2sQr4NsPr+cDmTG1igIHu83eAg+7nbQSu9Hp/l8aHuDvbGGOKnJ0iGWN8xhKMMcZnLMEYY3zGEowxxmeCcm9S8tWqVUvDw8O9DsOYUmvdunWHVbV25uU+TTAicinOmI9AYIaqPp3p/XLA2zjjRY4Ao1Q1RkQG4YzCDQHOAg+p6hJ3nQuBWUB5nP4REzWXW2Hh4eFER0cX5VczxmQgIpl7WgM+PEUSkUDgVeAynE5RN4hIu0zNxgLHVLUF8C+cHpUAh3HGmXQERgNzMqwzDRgHtHQfl/rqOxhjCseX12C6AdtUdYeqngXmAiMytRnB7+NN5gOXiIio6o+qmj7MfyMQKiLlRKQ+UEVVV7lHLW8DNmeIMX7KlwmmIbA7w+s4fh+cdl4bVU3B6aVZM1Oba4AfVfWM2z4ul20CICLj3Pk+og8dOlTgL2FMWfRj7DHGvR3Noo37C7UdX16DkSyWZb5WkmMbEWmPc9o0OB/bdBaqTgemA0RGRp7XJjk5mbi4OJKSkrJa3WQhNDSUsLAwgoPzNKTIlDCqyvLfDjM1ahurdxylavlgBratW6ht+jLBxAGNMrwO4/zZzdLbxLmD6KriTLWIiIThDKi7RVW3Z2gflss28xZcXByVK1cmPDwckazylslIVTly5AhxcXE0bdrU63BMEUpNUxZu2M+0pdvYsOcEdauU429D23JD98ZUKle4FOHLBLMWaCkiTXGG1V8P3JipzQKci7irgGuBJaqqIlINZ7DbJFX9Lr2xqu4TkQQR6YEzIdEtwL8LElxSUpIll3wQEWrWrImdbpYeZ1JS+eiHPUxftoOdh0/RtFZFnrmmI1de0JByQYFF8hk+SzCqmiIiE4BFOLep31TVjSLyBBCtqguAmcAcEdmGc+Ryvbv6BKAF8KiIPOouG6yqB4G7+f029Zfuo0AsueSP7a/S4eSZFN5bs4sZy3dyMOEMHRtWZepNXRjSvh6BAUX7O/ZpPxhV/QKnr0rGZY9leJ4EjMxivX/gzGaf1TajceZdNcbkw5GTZ5i1MobZK2M4kZRCr+Y1eeG6CPq0qOWz/zzKRE/eku6xxx7j4osvZuBAqxlm8i/uWCJvLNvBvOjdnElJY0i7etzVrzmdG1Xz+WdbgvFzqampPPHEE/leJzCwaM6hTcn124EEpkVt538/7yVA4MrODbmzb3Na1KlUbDHYYEcPxcTE0KZNG0aPHk2nTp249tprSUxMJDw8nCeeeII+ffrw4YcfMmbMGObPnw/A4sWLueCCC+jYsSO33XYbZ86cAThvHVN2rdt1jNtnRzPoX8v4csN+RvcMZ+lD/XluZESxJhewIxgAFi5cyP79hetQlFm9evW49NLcRzFs3bqVmTNn0rt3b2677TamTp0KOH1OVqxYcS4+cO58jRkzhsWLF9OqVStuueUWpk2bxv3333/eOqZsUVWW/nqIqVHb+X7nUapVCOb+gS0Z3TOc6hVDPIvLjmA81qhRI3r37g3AzTfffC5BjBo16ry2W7dupWnTprRq1QqA0aNHs2zZ73NaZ7WOKd1SUtNY8PNehr6ygjFvrWX30UQeG9aOlQ8P4P6BrTxNLmBHMAB5OtLwlcxX79NfV6x4fgWN3OZPzmodUzolJTt9WF5ftp1dRxJpVrsiz17biSs7NyQkyH+OGyzBeCw2NpZVq1bRs2dP3n//ffr06cOPP/6YZds2bdoQExPDtm3baNGiBXPmzKFv377FHLHxUkJSMu+uiWXmip0cSjhDRFhVJt18IYPb1SWgiPuwFAX/SXVlVNu2bZk9ezadOnXi6NGj3H333dm2DQ0N5a233mLkyJF07NiRgIAA7rrrrmKM1njl8MkzPLdoC72eXsLTX26hdd3KvHd7dz65pzeXdqjnl8kF7AjGcwEBAbz22mt/WBYTE/OH17NmzTr3/JJLLsnyCCfzOqZ02H00kTeW72De2t2cTU3j0vb1uLtfczqF+b4PS1GwBGOMH9qy/wSvRW3n0/X7CBC4pksY4y5uRrPaxXububAswXgoPDycDRs2eB2G8SPRMUeZGrWdJVsOUiEkkFt7hXP7Rc2oVzXU69AKxBKMMR5TVaK2HmJa1Ha+jzlK9QrB/HlgK0b3akK1Ct7eZi4sSzDGeCQlNY3Pf9nHtKjtbNmfQMNq5ZlyRTuu69qICiGl40+zdHwLY0qQpORUPlwXx/Rl29l99DQt61TihZERDO/cgODA0nVj11/LltTEmQS8KzBLVSdkWCcKqA+cdhelzxNjjF87kZTMO6t38eaKnRw+eZbOjarx6OXtGNjWP/uwFAV/LVuSBDwK/CWbzd+kqp3dR4lNLrfddht16tShQwffT2+zd+9err32WgCioqIYNmwY4NwCnzBhQk6rmkI6mJDEMwu30PupJTy7cCvtGlTl/Tt68PH4Xgxu7799WIqCL49gzpUtARCR9LIlmzK0GQFMcZ/PB/7jli05BawQkRY+jM9zY8aMYcKECdxyyy0+/6wGDRqcG5FtiseuI6eYvmwHH66LIzk1jaEd63N33+Z0aFjV69CKTUkoW5KVt0TkJxF5VErwPI4XX3wxNWrUyPb9+Ph4wsPDSUtLAyAxMZFGjRqRnJzMG2+8QdeuXYmIiOCaa64hMTERcJLWfffdR69evWjWrNm5pBITE5PrkdKnn35K9+7dueCCCxg4cCAHDhwoom9atmzed4L73v+R/s9H8WF0HNd0aciSB/vx6o1dylRyAT8vW5KNm1R1j4hUBv4L/AnnOs4fNywyDqcCJI0bN85xg49/upFNe0/k8rH5065BFf5+RftCbaNq1apERESwdOlS+vfvz6effsqQIUMIDg7m6quv5o477gBg8uTJzJw5k3vvvReAffv2sWLFCrZs2cLw4cPPnRrlpk+fPqxevRoRYcaMGTz77LO88MILhfoOZYWqsjbmGNOitvHt1kNUDAnkjouacVufptStUjL7sBQFvy1bkh1V3eP+TBCR93BOxc5LMLnVRSopRo0axbx58+jfvz9z585l/PjxAGzYsIHJkydz/PhxTp48yZAhQ86tc+WVVxIQEEC7du3ydRQSFxfHqFGj2LdvH2fPnrXyJHmQlqZ8u/UgU6O2s27XMWpUDOEvg1vxpx7hVK1g9aP8smxJdht0k1A1VT0sIsHAMOCbwgZa2CMNXxo+fDiTJk3i6NGjrFu3jgEDBgDOqdAnn3xCREQEs2bNIioq6tw65cqVO/c8tykeMrr33nt54IEHGD58OFFRUUyZMqWovkapk5KaxmfrnT4sWw84fVgeH96e6yIbUT7EpitN569lSxCRGKAKECIiV+JUd9wFLHKTSyBOcnnDV9/BH1SqVIlu3boxceJEhg0bdm6u3YSEBOrXr09ycjLvvvsuDRtmWUE3X+Lj489tZ/bs2bm0LptOn03lw3W7mb5sB3HHTtOqbiX+NSqCYZ1KXx+WouCXZUvc98Kz2eyFRRWf12644QaioqI4fPgwYWFhPP7444wdO/a8dqNGjWLkyJF/OEp58skn6d69O02aNKFjx44kJCQUOp4pU6YwcuRIGjZsSI8ePdi5c2eht1laxJ9OZs6qGN76LoYjp87SpXE1plzRngFt6pTq28yFJfk5hC6pIiMjNTo6+g/LNm/eTNu2bT2KqOQqa/vt4IkkZq7YybtrYjl5JoV+rWszvl8LuoZXt0J0GYjIOlWNzLzchgoYk4WYw6d4fdkO/rsujpS0NC7v1IC7+zanXYMqXodWoliCMSaDDXvieW3pdr74ZR9BgQFcGxnGnRc3o0lNm++4IMp0glFVO8zNh9J6Oq2qrNnpzMOy7NdDVCoXxB0XN2Ns76bUKcN9WIpCmU0woaGhHDlyhJo1a1qSyQNV5ciRI4SGlp4/uLQ0ZfGWg0yN2saPscepVSmEh4a05uYeTaha3vqwFIUym2DCwsKIi4vj0KFDXodSYoSGhhIWFuZ1GIWWnJrGgp/28trS7fx28CRh1cvz5Ij2jIxsRGiw9WEpSmU2wQQHB1tP1TLm9NlU5q2N5Y3lO9lz/DSt61bmpVGdGdapPkHWh8UnymyCMWVHfGIyb6+K4a2VMRw9dZau4dV5YoTTh8VOj33LEowptQ6cSGLG8h28tyaWU2dTGdCmDnf3a07X8OxHsJuiZQnGlDo7Dp1k+rIdfPTDHlLS0rgiogF39W1O2/rWh6W4WYIxpcaGPfFMi9rOFxv2ERwYwHVdwxh3UXMa16zgdWhlVo4JRkRmqeoY9/loVbURcMavqCqrdhxhWtR2lv92mMrlgrirb3Nu692U2pXL5b4B41O5HcFEZHg+EbAEY/xCWpry1aYDTFu6nZ93H6dWpXL836VtuKlHY6qEWh8Wf5FbgimdXTdNiXU2JY3//bSH15ZuZ/uhUzSuUYF/XNmBay8Msz4sfii3BBMmIq/gTG2Z/vwcVb3PZ5EZk0Hi2RTmfr+bGct3sDc+ibb1q/DKDRcwtEM968Pix3JLMA9leB6dbats+Kgu0oXALKA8zlwzE3OaBc+UbMcTzzJrZQyzVsZwPDGZbuE1+H9XdaRf69rWh6UEyDHBFOaiboa6SINw5t5dKyILVDVj2ZJzdZFE5Hqcukij+L0uUgf3kdE0nMm8V+MkmEuBLwsap/FP++PdPizfx5J4NpWBbetwV9/mRFoflhIlt7tIC3J6X1WH5/B2kddFEpH6QBVVXeW+fhu4kiJIMAsXLmT//v2F3YwppMNnAvnuaHl+jg9FgY5VztC7YSJ19RAbojaywesAy5h69epx6aWXFnj93E6ReuLULXofWEPWZUayk1VdpO7ZtXHn8E2vi3Q4h23GZdpmlpPR5qdsifHentNBrDhSgc0nQwgSuLBaEr1qJFI9JM3r0Ewh5JZg6uGc4tyAUxHgc+B9Vd2Yh237oi5Sntvnt2xJYbK0KRhVZeX2I0yN2sZ3u45QJTSICf3DGdM7nFqVrA9LaZDbNZhUYCGw0L0gewMQJSJPqOq/c9m2L+oixbnbyWmbxs+lpilfbdzPtKXbWR8XT53K5XhkaBtu6NaYytaHpVTJdaiAm1gux0ku4cArwEd52HaR10VS1X0ikiAiPXBO2W4Bckt0xk+cTUnjkx/38Nqy7ew4dIomNSvwz6s6cs2FDSkXZH1YSqPcLvLOxrmL8yXwuKrm+RqbL+oiuXeg7ub329RfYneQ/N6pMym8/30sM5bvZP+JJNrVr8J/bryAyzrUJ9BKfpRqOZYtEZE04FSGRemNBVBVLRHDU7MqW2J87+ips8xeGcPsVU4flu5NazC+fwsublnL+rCUMgUqW6Kq1kXS5Nue46eZsXwHc7/fzenkVAa2rcv4/s3p0ri616GZYpbbKVIocBfQAliPc5qTUhyBmZJn28EEXlu6g09+3APA8M5OLaGWdSt7HJnxSm4XeWcDycByYCjQHmdUtTHn/LT7ONOitvHVpgOUCwrg5h5NuP2ipoRVt3lYyrrcEkw7Ve0IICIzge99H5IpCVSVFdsOM/Xb7azacYSq5YO5t38LxvRuSo2KIV6HZ/xEbgkmOf2Je1fIx+EYf5eapizcsJ9pS7exYc8J6lYpx9+GtuWG7o2pVM4mSDR/lOuEUyJywn0uQHn3dYm6i2QK70xKKh//sIfXl+1g5+FTNKtVkaev7shVXawPi8lebneR7F9OGXfyTArvrdnFjOU7OZhwho4NqzL1pi4MaV/P+rCYXNkxrcnSkZNnmLUyhtkrYziRlEKv5jV54boI+rSwPiwm7yzBmD+IO5bIjOU7mbs2lqTkNIa0r8vd/VrQuVE1r0MzJZAlGAPArwcSeC1qO//7eS8CXHVBQ+7s24wWdawPiyk4SzBl3A+xx5gWtZ2vNx2gfHAgo3uGc/tFTWlQrbzXoZlSwBJMGaSqLPvtMFO/3caanUepViGYiZe0ZEyvcKpbHxZThCzBlCGpacoXv+xjWtR2Nu07Qf2qoUy+vC03dGtMRevDYnzA/lWVAUnJqXz0wx6mL9tOzJFEmtWuyLPXduLKzg0JCbLxrMZ3fJpgClq2xH1vEk7VgVTgPlVd5C6PARLc5SlZDRE3joSkZN5dE8vMFTs5lHCGTmFVmXZTFwZbHxZTTHyWYApTtkRE2uFMPtUeaAB8IyKt3Ck8AfqranYTg5d5h0+e4a3vdvL2ql0kJKXQp0UtXhrVmV7Na1ofFlOsfHkEU+CyJe7yuap6BtjpznjXDWdqTZON3UcTeWP5Duat3c3Z1DQubV+Pu/s1p1OY9WEx3vBlgilM2ZKGOIXVMq6bXp5Ega9ERIHX3eoB5ylLZUu27k9gWtQ2Pl2/jwCBqy8IY1zfZjSvXcnr0EwZ58sEU5iyJTmt21tV94pIHeBrEdmiqsvOa5zPsiUl0bpdR5n67XYWbzlIhZBAbu0Vzu0XNaNe1VCvQzMG8G2CKUzZkmzXVdX0nwdF5GOcU6fzEkxppapEbT3EtKjtfB9zlOoVgnlgUCtu6dmEahWsD4vxL75MMAUuW+KWrH1PRF7EucjbEvheRCoCAaqa4D4fDDzhw+/gN1JS0/jc7cOyZX8CDaqG8vcr2jGqayMqhFhvA+OffPYvszBlS9x2H+BcEE4B7lHVVBGpC3zs3gkJAt5T1YW++g7+ICk5lfnr4pi+bAexRxNpUacSz4+MYHhEA+vDYvxejmVLSouSWLbkRFIy7652+rAcPnmGiEbVGN+vOYPa1iXA+rAYP1OgsiWm+B1KOMOb3+3knVW7SDiTwkUta3F3v870bGZ9WEzJYwnGT8QeSWT68u18EB1HcmoaQzvU5+5+zenQsKrXoRlTYJZgPLZ53wmmRW3ns/V7CQoI4OouDbmzb3Oa1qrodWjGFJolGI98v/Mo06K28e3WQ1QMCeT2i5oxtk9T6laxPiym9LAEU4xUlSVbDjItajvRu45Ro2IIDw5qxS09w6laIdjr8IwpcpZgikFKahqfrXf6sGw9kEDDauWZckU7RnVtTPkQK9xgSi9LMD6UlJzKh9G7eX3ZDuKOnaZV3Uq8MDKC4Z0bEBxofVhM6WcJxgfiTyfzzupdvLliJ0dOnaVL42pMuaI9A9rUsT4spkyxBFOEDiYkMXPFTt5dHcvJMyn0a12bu/s2p1vTGtaHxZRJlmCKQMzhU0xfvoP56+JISU3j8k4NuKtvM9o3sD4spmyzBFMIG/fGMy1qO1/8so+ggACujQxj3EXNCLc+LMYAlmDyTVVZs/Mo06K2s/TXQ1QqF8QdFzdjbO+m1LE+LMb8gSWYPEpLUxZvOcjUqG38GHucmhVDeGhIa27u0YSq5a0PizFZsQSTi+TUNBb8tJfXlm7nt4MnCatenidHtGdkZCNCg60PizE5sQSTjdNnU5m3NpY3lu9kz/HTtKlXmZev78zlHesTZH1YjMmTklgXKcdtFlZ8YjJvr4rhrZUxHD11lsgm1XlihNOHxW41G5M/JaoukrtObtsskAMn0vuw7OLU2VT6t67N+P4t6Bpeo7CbNqbMKml1kcjDNvMtOuYoN76xhpS0NK6IaMBdfZvTtn6VwmzSGEPJrIuU2zaB/NVF6hRWjVv7hHNTtyY0rlkhx7bGmLzz5dVKX9RFyss2nYWq01U1UlUja9eunWOgIUEBTLqsrSUXY4qYLxNMfuoikce6SHnZpjHGT/gywZyriyQiITgXbRdkapNeFwky1EVyl18vIuXcukotge/zuE1jjJ8oUXWRALLaZm6xrFu37rCI7MpD2LWAw/n9rsXEYis4f47Pn2ODvMfXJKuFZaIuUl6JSHRWtV38gcVWcP4cnz/HBoWPz7qkGmN8xhKMMcZnLMH80XSvA8iBxVZw/hyfP8cGhYzPrsEYY3zGjmCMMT5jCcYY4zOWYHCmgBCRrSKyTUQe9jiWRiLyrYhsFpGNIjLRXV5DRL4Wkd/cn9U9jjNQRH4Ukc/c101FZI0b3zy3I6QXcVUTkfkissXdhz39ad+JyJ/d3+sGEXlfREK92nci8qaIHBSRDRmWZbmvxPGK+zeyXkS65OUzynyCyTCtxGVAO+AGd7oIr6QAD6pqW6AHcI8bz8PAYlVtCSx2X3tpIrA5w+tngH+58R3DmYrDCy8DC1W1DRCBE6Nf7DsRaQjcB0SqageczqLp05R4se9mAZdmWpbdvroMp0d9S5xBxNPy9AmqWqYfQE9gUYbXk4BJXseVIZ7/4cx/sxWo7y6rD2z1MKYw9x/fAOAznEGoh4GgrPZpMcZVBdiJe/Miw3K/2Hf8PntADZxe9J8BQ7zcd0A4sCG3fQW8DtyQVbucHmX+CIasp5VomE3bYiUi4cAFwBqgrqruA3B/1vEuMl4C/gqkua9rAsdVNcV97dU+bAYcAt5yT99miEhF/GTfqeoe4HkgFtgHxAPr8I99ly67fVWgvxNLMPmYAqI4iUgl4L/A/ap6wut40onIMOCgqq7LuDiLpl7swyCgCzBNVS8ATuH9qeQ57vWMEUBTnJkaK+KcemTm+b+/LBTod2wJxg+ngBCRYJzk8q6qfuQuPiAi9d336wMHPQqvNzBcRGKAuTinSS8B1dwpN8C7fRgHxKnqGvf1fJyE4y/7biCwU1UPqWoy8BHQC//Yd+my21cF+juxBONnU0C4U4bOBDar6osZ3so4tcVonGszxU5VJ6lqmKqG4+yrJap6E/AtzpQbnsWnqvuB3SLS2l10Cc6IfL/YdzinRj1EpIL7e06Pz/N9l0F2+2oBcIt7N6kHEJ9+KpUjLy52+dsDGAr8CmwH/uZxLH1wDj3XAz+5j6E41zkWA7+5P2v4wX7rB3zmPm+GM2fPNuBDoJxHMXUGot399wlQ3Z/2HfA4sAXYAMwBynm174D3ca4FJeMcoYzNbl/hnCK96v6N/IJzJyzXz7ChAsYYn7FTJGOMz1iCMcb4jCUYY4zPWIIxxviMJRhjjM9YgjHG+IwlGGOMz1iCMcb4jCUYY4zPWIIxxviMJRhjjM9YgjHG+ExQ7k2KnohcijN3aiAwQ1WfzvR+OeBt4ELgCDBKVWNEpBu/F4ISYIqqfpzb59WqVUvDw8OL8BsYYzJat27dYVWtnXl5sSeYDJNsD8IZIr5WRBao6qYMzcYCx1S1hYikT4o8CmeIe6SqpriT4fwsIp/q79MNZik8PJzo6GiffB9jDIjIrqyWe3GK1A3Ypqo7VPUszqxoIzK1GQHMdp/PBy4REVHVxAzJJBT/nFrQGOPyIsHkZfLgc23chBKPMxEOItJdRDbiTHpzV3ZHLyIyTkSiRST60KFDRfwVjCkDEg5AIeeL8iLB5GXy4GzbqOoaVW0PdAUmiUhoVh+iqtNVNVJVI2vXPu/U0BiTk/2/wGu9YdlzhdqMFxd58zJ5cHqbOHcy5KrA0YwNVHWziJwCOuBMkZgvycnJxMXFkZSUlN9Vy6zQ0FDCwsIIDg72OhTjS7Gr4d3roFwlaHdloTblRYI5N8k2sAdn4ugbM7VJn3h4Fc5kyEtUVd11drsXeZsArYGYggQRFxdH5cqVCQ8Px5l/2eREVTly5AhxcXE0bdrU63CMr2z7Bub9CSrXh1s+gWqNC7W5Yj9Fcq+ZTAAW4ZT1/EBVN4rIEyIy3G02E6gpItuAB/i9tk0fnDtHPwEfA+NV9XBB4khKSqJmzZqWXPJIRKhZs6Yd8ZVmG/4L710PNZvDbQsLnVzAo34wqvoF8EWmZY9leJ4EjMxivTk4M7EXCUsu+WP7qxRbOwM+/ws07gk3zoXQqkWyWU8SjDHGT6jC0mcg6ilodRmMfAuCyxfZ5m2ogEd2795N//79adu2Le3bt+fll1/26eft3buXa691antFRUUxbNgwAGbNmsWECRN8+tnGT6Wlwhd/cZJLxI0w6p0iTS5gRzCeCQoK4oUXXtHduU0AACAASURBVKBLly4kJCRw4YUXMmjQINq1a+eTz2vQoAHz58/3ybZNCZScBB+Pg03/g173waAnwAenwHYE45H69evTpUsXACpXrkzbtm3Zs2fPH9rEx8cTHh5OWloaAImJiTRq1Ijk5GTeeOMNunbtSkREBNdccw2JiYkAjBkzhvvuu49evXrRrFmzc0klJiaGDh065BjTp59+Svfu3bngggsYOHAgBw4cKOqvbfzB6ePwzjVOchn8Dxj8pE+SC9gRjOPLh52ORUWpXke47Onc2+H88f/444907979D8urVq1KREQES5cupX///nz66acMGTKE4OBgrr76au644w4AJk+ezMyZM7n33nsB2LdvHytWrGDLli0MHz783KlRbvr06cPq1asREWbMmMGzzz7LCy+8kI8vbfzeib3wzrVw+Fe4egZ0Ou9eSpGyBOOxkydPcs011/DSSy9RpUqV894fNWoU8+bNo3///sydO5fx48cDsGHDBiZPnszx48c5efIkQ4YMObfOlVdeSUBAAO3atcvXUUhcXByjRo1i3759nD171vq7lDYHNzvJJSkebvoQmvf3+UdagoE8H2kUteTkZK655hpuuukmrr766izbDB8+nEmTJnH06FHWrVvHgAEDAOdU6JNPPiEiIoJZs2YRFRV1bp1y5cqde56f2uP33nsvDzzwAMOHDycqKoopU6YU6HsZP7RzOcy9CYJD4dYvoH6nYvlYuwbjEVVl7NixtG3blgceeCDbdpUqVaJbt25MnDiRYcOGERgYCEBCQgL169cnOTmZd999t0hiio+Pp2FDZ9zp7Nmzc2ltSoxf5sM7V0PlejD262JLLmAJxjPfffcdc+bMYcmSJXTu3JnOnTvzxRdfZNl21KhRvPPOO4waNercsieffJLu3bszaNAg2rRpUyQxTZkyhZEjR3LRRRdRq1atItmm8ZAqLH8B/jsWwrrC2EVQvUmxhiD5OYQuqSIjIzXzhFObN2+mbdu2HkVUctl+KyFSk+HzB+GH2dBxJIx4FYLK5b5eAYnIOlWNzLzcrsEYU9okxcMHo2HHt3DRg9B/MgR4c7JiCcaY0uR4LLw3yrkNPfw/0OVPnoZTphOMqtoAvnwoC6fTJVpcNLx/PaSchZv/C836eR1R2b3IGxoaypEjR+yPJo/S54MJDc1yAkHjtV/mw6zLIaQi3P6NXyQXKMNHMGFhYcTFxWHz9eZd+ox2xo9kHA3duJczYLFiTa+jOqfMJpjg4GDrqWpKtrOJ8L/xsPFjZzT0FS/59E5RQZTZBGNMiRa/B+beCPt+dkZC97rPZwMWC8MSjDElze61MO8mOHsKbngfWl/mdUTZKrMXeY0pkX58F2YNdSaGuv0bv04uYEcwxpQMqcnw1aOwZppzh+jat6BCDa+jypUlGGP83anD8OEYiFkOPcbDoCchsGT86ZaMKI0pq/b+6NQpOnkQrnodIq73OqJ8sQRjjL/66T349H6oWNupU9Swi9cR5ZsnF3lF5FIR2Soi20Tk4SzeLyci89z314hIuLt8kIisE5Ff3J8Dijt2Y3wu5Sx89gB8cjc07g53Li2RyQU8OIIRkUDgVWAQTg3qtSKyQFU3ZWg2Fjimqi1E5HrgGWAUcBi4QlX3ikgHnOqQDYv3GxjjQ/FxzkjoPdHQeyIMeKzEXG/JiheRdwO2qeoOABGZC4wAMiaYEcAU9/l84D8iIqr6Y4Y2G4FQESmnqmd8H7YxPrZ9Cfz3ducI5rq3od0IryMqNC9OkRoCuzO8juP8o5Bzbdxa1vFA5gEW1wA/ZpdcRGSciESLSLSNNzJ+LS0Nop6BOVc711vGfVsqkgt4cwSTVX/mzEOac2wjIu1xTpsGZ/chqjodmA7OjHb5D9OYYnDqMHw0DrYvhk6jYNi/nBHRpYQXCSYOaJThdRiwN5s2cSISBFQFjgKISBjwMXCLqm73fbjG+EjsavjwVkg84iSWC2/1y/FEheHFKdJaoKWINBWREOB6YEGmNguA0e7za4ElqqoiUg34HJikqt8VW8TGFKW0NFjxL3hrqDP6+favIfK2UpdcwIMjGFVNEZEJOHeAAoE3VXWjiDwBRKvqAmAmMEdEtuEcuaT3LpoAtAAeFZFH3WWDVfVg8X4LYwro5CH4+E7nlKj9VXDFKxB6fsG90qLMVhUwptjtiHKutyTFw5B/lqqjFqsqYIxXUpPh2386p0W1WsKfPoa67b2OqlhYgjHGl47ugP/e4XSc63ILXPp0qbpLlBtLMMb4gir8PBe+eAgkwJleoUPW9cdLswLfRRKRWRmej86hqTFly+ljTrnWT+6Ceh3h7u/KZHKBwh3BRGR4PhGwaunG7FjqDFI8eQAueQx63w8BgV5H5ZnCJJjSf/vJmLxKToLFT8DqV6FmCxj7dYkdAV2UCpNgwkTkFZxu/enPz1HV+woVmTElxd4f4eO74NAW6Hq7M8t/GbqQm5PCJJiHMjy3Tiam7ElNhuUvwLLnnEGKN38ELS7xOiq/UuAEo6p2zcWUXQc2Okct+9c7gxQvewbKV/c6Kr9T4AQjIpnHD/2Bqg4v6LaN8VupyfDdS870CuWrOaVa217hdVR+qzCnSD1x5mx5H1hD1lMsGFN67P8F/nePU02x/dUw9Hm/qgPtjwqTYOrhTHt5A3Ajzijn91V1Y1EEZozfSDkDy56HFS86p0HXzYF2doCeF4W5BpMKLAQWikg5nEQTJSJPqOq/iypAYzwVuwYW3AuHt0Kn6+HSp0pEwTN/UaihAm5iuRwnuYQDrwAfFT4sYzyWFO/0a1k7E6qGwU3zoeUgr6MqcQpzkXc20AH4EnhcVTcUWVTGeEUVtnwGX/wVEvZB97tgwN+gXGWvIyuRCnME8yfgFNAKmCgi6T17BVBVLb2z6JjS6XgsfPl/sPULqNvRuUMUdqHXUZVohbkG40nRNmOKXMpZWD0Vlj7jvB70pFMDugTXI/IXhTlFCgXuwpnCcj3O1JcpRRWYMcVi53L4/EHnIm7roXDZs1CtUe7rmTwpTIqeDSQDy4GhQHucUdXG+L8T++CrybBhPlRrDDfMg9aXeh1VqVOYBNNOVTsCiMhM4PuiCckYH0o5A6unOeOHUpOh7/9Bnz9DcHmvIyuVCpNgktOfuJUCiiAcY3xEFX5dBIsegaPbndOhIf8PajTzOrJSrVATTonICfe5AOXd13YXyfiXg1ucxLJ9MdRsaX1ailGB7wSpaqCqVnEflVU1KMPzHJOLiFwqIltFZJuIPJzF++VEZJ77/hoRCXeX1xSRb0XkpIj8p6CxmzLi1GH4/C8wrRfERTulQsavsuRSjIr9PpyIBAKv4oxjigPWisgCVd2UodlY4JiqthCR63HqUI8CkoBHcTr4dSjeyE2JkXwa1rwGy1+Es6cg8lbo94gNTPSAFzf6uwHbVHUHgIjMBUYAGRPMCGCK+3w+8B8REVU9BawQkRbFGK8pKdLS4JcPYPGTcCIOWg6BwU9C7dZeR1ZmeZFgGuJM85AuDuieXRv3AnI8UBM4nNcPEZFxwDiAxo0bFyZe4+9UYdti+GYKHPgF6kfAVdOg6cVeR1bmeZFgsrrdlHkC8by0yZGqTgemg1M6Nj/rmhJk9/fwzeOwawVUawLXzHTmagmwjub+wIsEEwdk7CoZBuzNpk2ciAQBVYGjxROeKRH2/wJL/h/8+qUzH+7Q56HLaAgK8Toyk4EXCWYt0FJEmgJ7gOtxJqzKaAEwGlgFXAssUVU7CjHOLeeop2DTJxBaFQY86ox4LlfJ68hMFoo9wbjXVCYAi4BAnDFMG0XkCSBaVRcAM4E5IrIN58jl+vT1RSQGqAKEiMiVwOBMd6BMaXRoqzMYccNHTkmQix+CnvfYRNt+TsrCgUFkZKRGR1tllRLpwCanW//GjyG4AnQfBz3vtVvOfkZE1qlqZOblNh7d+Kc9Pzg1h7Z8BiGVnPFCPe+BirW8jszkgyUY4z9UIWa500Fux7fONZa+/+dcY7F5cEskSzDGe2mpzpHKipdg7w9QsQ4MnAKRYyHUhrSVZJZgjHfOnoKf3oNVr8KxnVC9KQz7F0TcCMGhXkdnioAlGFP84uPg+zdg3SxIOg4NI2HQ49BmGAQEeh2dKUKWYEzxUIXda5xBiJsWAOoklJ4ToHHmkSKmtLAEY3zr7Cn4ZT6sfcPpfRtaFXrcDd3GQfUmXkdnfMwSjPGNg1tg3Vvw0/twJh7qtHeur3Qa5XSUM2WCJRhTdM6egk3/gx/ehthVEBgC7UZA5G3QuCfYtKpljiUYUziqsGcd/DjH6cZ/5gTUaO7UFup8o3WMK+MswZiCiY+D9R/Az+/D4V8hqLxztNLlFmjSy45WDGAJxuTH6eOweYGTWGJWAAqNesAVLztzsFinOJOJJRiTszMn4deFzunPtq8h9axT6qPfw9DpOiv7YXJkCcacL+mEU0No0yew7RtISYLKDaDrHdDxWmhwgZ0CmTyxBGMcCfudI5XNn8HOpc6RSqV6zjWVdlc6d4FsGkqTT5ZgyipV2L8efv3KmXZyzzpnefVwpxNc2+EQ1tWSiikUSzBlSeJRZxqEbUucU5+T+53lDbpA/8nQ+jKo295Of0yRsQRTmp1NdMb/7FwKO6Jg70+AOt31mw+AFoOcKoeV6ngdqSmlLMGUJmcSnDIeu1bCru+ccqlpyRAQ5Jzu9JvkJJaGXWzUsikWlmBKKlVnDpW4aCep7F4DBzaApoEEOsXHetwNTftC4x42677xhCWYkkAVjsc6F2X3/gR7f3Rmfjt9zHk/pBI0vNCZab9xDwjrZgnF+AVLMP7m9DGnRMfBzc7jwEanHGpSvPO+BEKdds5cKmGRTmKp085OeYxf8iTBiMilwMs4dZFmqOrTmd4vB7wNXAgcAUapaoz73iRgLJAK3Keqi4ox9KJx5iQc3wXHYuDoDjiyDQ5vgyO/wckDv7cLqQR12jrd8Ot1dE576rSDkAqehW5MfhR7ghGRQOBVYBBOidi1IrIgU/G0scAxVW0hItcDzwCjRKQdThG29kAD4BsRaaWqqcX7LbKRluYcgZw6CCcPOskiYR+c2OsMDjyxxznVSTzyx/XK14CaLZw7OjVbOkmldmuo2tj6oZgSzYsjmG7ANlXdASAic4ERQMYEMwKY4j6fD/xHRMRdPldVzwA73cqP3XBKzBZcUjzEfAeoc5E0LRXSUiA12enRmnIGUk5D8mk4e9I5Ajl70lkvKd5JKolH4fRRZ/3MQipBlYZQrZFzFFKtCVRrDDWaOhNdW0kOU0p5kWAaArszvI4DMk/Keq6NW2o2HqjpLl+dad2GWX2IiIwDxgE0btw454iO7oS5N+Qt+qDyzgXUkEpOf5LQqs5pS/nqztwnFWo6xdgr1XG62leuZ6OMTZnlRYLJqpto5vq12bXJy7rOQtXpwHRwSsfmGFGtVjAuCiTA+YiAIOeiaWCwMytbUHkICnFKl9rFVGPyzIsEEwc0yvA6DNibTZs4EQkCqgJH87hu/oVUcEYIG2OKlBdXENcCLUWkqYiE4Fy0XZCpzQJgtPv8WmCJqqq7/HoRKSciTYGWwPfFFLcxJp+K/QjGvaYyAViEc5v6TVXdKCJPANGqugCYCcxxL+IexUlCuO0+wLkgnALc4zd3kIwx5xHnwKB0E5FDwK48NK0FHPZxOAVlsRWcP8fnz7FB3uNroqq1My8sEwkmr0QkWlUjvY4jKxZbwflzfP4cGxQ+PuvFZYzxGUswxhifsQTzR9O9DiAHFlvB+XN8/hwbFDI+uwZjjPEZO4IxxviMJRhjjM9YgsGZn0ZEtorINhF52ONYGonItyKyWUQ2ishEd3kNEflaRH5zf1b3OM5AEflRRD5zXzcVkTVufPPcXtpexFVNROaLyBZ3H/b0p30nIn92f68bROR9EQn1at+JyJsiclBENmRYluW+Escr7t/IehHpkpfPKPMJJsP8NJcB7YAb3HlnvJICPKiqbYEewD1uPA8Di1W1JbDYfe2licDmDK+fAf7lxncMZ04fL7wMLFTVNkAETox+se9EpCFwHxCpqh1werKnz3fkxb6bBVyaaVl2++oynKE5LXFmKZiWp09Q1TL9AHoCizK8ngRM8jquDPH8D2dyrq1AfXdZfWCrhzGFuf/4BgCf4YxyPwwEZbVPizGuKsBO3JsXGZb7xb7j92lIauAM0/kMGOLlvgPCgQ257SvgdeCGrNrl9CjzRzBkPT9NlnPMFDcRCQcuANYAdVV1H4D708tiRi8BfwXSZ9eqCRxX1RT3tVf7sBlwCHjLPX2bISIV8ZN9p6p7gOeBWGAfEA+swz/2Xbrs9lWB/k4sweRjjpniJCKVgP8C96vqCa/jSSciw4CDqrou4+IsmnqxD4OALsA0Vb0AOIX3p5LnuNczRgBNcaZ8rYhz6pGZ5//+slCg37ElGF/NMVMIIhKMk1zeVdWP3MUHRKS++3594KBH4fUGhotIDDAX5zTpJaCaO3cPeLcP44A4VV3jvp6Pk3D8Zd8NBHaq6iFVTQY+AnrhH/suXXb7qkB/J5Zg8jY/TbFx5x6eCWxW1RczvJVxjpzRONdmip2qTlLVMFUNx9lXS1T1JuBbnLl7PItPVfcDu0WktbvoEpypPfxi3+GcGvUQkQru7zk9Ps/3XQbZ7asFwC3u3aQeQHz6qVSOvLjY5W8PYCjwK7Ad+JvHsfTBOfRcD/zkPobiXOdYDPzm/qzhB/utH/CZ+7wZzuRf24APgXIexdQZiHb33ydAdX/ad8DjwBZgAzAHKOfVvgPex7kWlIxzhDI2u32Fc4r0qvs38gvOnbBcP8OGChhjfMZOkYwxPmMJxhjjM5ZgjDE+YwnGGOMzlmCMMT5jCcYY4zOWYIwxPmMJxhjjM5ZgjDE+YwnGGOMzlmCMMT5jCcYY4zNBuTcp+WrVqqXh4eFeh2FMqbVu3brDqlo78/IykWDCw8OJjo72OgxjSi0R2ZXVcjtFMsb4TJk4gjEmJ6rKkaQjHE06SvyZeJLTklFVQoNCqRRciToV6lCtXDWcSehMfliCMWVKmqbx27HfiD4QzcbDG9l6bCuxJ2JJSk3Kcb3yQeVpUa0FrWu0pkudLnSr1426FesWU9QllyUYU+qdST3DirgVLNm9hOVxyzl25hgAtcvXpnWN1vSo34OGlRpSq3wtqparSkhgCIKQlJpEwtkEDpw6wJ6Te/jt2G8silnE/F/nA9C+ZnsGNRnEFc2voE4FL6vI+K8yMWVmZGSk2kXesmfjkY18uPVDvor5ioTkBKqEVOHisIvp1aAXkXUjqV+pfr63maZpbD26lZV7V/LNrm/YcGQDgRJIv0b9GNN+DJ3rdPbBN/F/IrJOVSPPW24JxpQmKWkpfBP7DW9vfJtfDv9C+aDyDGoyiMubXk63+t0ICijag/bYE7HM/20+H//2McfPHKdrva5M7DKRiNoRRfo5/s4SjCWYUi0lLYXPd3zO6+tfZ3fCbhpXbsyNbW/kiuZXUCWkis8/PzE5kfm/zufNDW9yJOkIQ8KH8FDkQ2XmOo0lGEswpZKqsiR2CS/98BIxJ2JoU6MNd3a6k/6N+hMYEFjs8SQmJ/LWxreYtWEWgQGB3N/lfq5rfR0BUrp7hFiCsQRT6mw+spmnv3+aHw7+QLOqzbivy30MaDTAL24n7z6xmydWP8Hqfavp1aAX/+j9D2pXOK+ja6lhCcYSTKmRcDaBl394mQ+2fkD10OpMuGACV7W4qsivrxSWqvLhrx/y3NrnKB9Unuf6Pkf3+t29Dssnskswpfu4zZQ6i2MXM+KTEXz464fc0OYGPr3qU0a2Gul3yQVARLiu9XXMGzaP6qHVufPrO3l749uUhf/U0/nfb8WYLBxPOs5T3z/FFzu/oHX11rwy4BU61OrgdVh50qxaM967/D3+tuJvPBf9HLEJsUzqNsmTa0TFzRKM8Xsr965k8orJHEs6xvjO47m94+0EBwR7HVa+VAyuyIv9XuTlH17mzQ1vcjDxIM/1fY5ygeW8Ds2n7BTJ+K3k1GSeX/s8d359J5VDKvPe5e9xd8TdJS65pAuQAP584Z+Z1G0SUbujmLB4AqdTTnsdlk/5NMGIyKUislVEtonIw1m8X05E5rnvrxGRcHd5TRH5VkROish/Mq0T5W7zJ/dhfbRLoT0n9zB64Whmb5rNqNajmDtsLm1rtvU6rCJxY9sbebL3k6zZt4bx34wv1UnGZ6dIIhIIvAoMAuKAtSKyQFU3ZWg2Fjimqi1E5HrgGWAUkAQ8CnRwH5ndpKp2W6iUWha3jEnLJ6GqvNjvRQY1GeR1SEVuRIsRBAUE8ciKR/jzt3/mlQGvEBIY4nVYRc6XRzDdgG2qukNVzwJzgRGZ2owAZrvP5wOXiIio6ilVXYGTaEwZkaZpTPt5GhMWT6B+xfrMGzavVCaXdJc3u5y/9/w73+39joeXP0xqWqrXIRU5XyaYhsDuDK/j3GVZtlHVFCAeqJmHbb/lnh49Ktn0qhKRcSISLSLRhw4dyn/0plidSj7FA1EPMPWnqQxrNow5Q+fQqEojr8PyuatbXs1DkQ/x9a6veT76ea/DKXK+vIuU1R9+5g4AeWmT2U2qukdEKgP/Bf4EvH3eRlSnA9PB6WiXe7jGK3tP7mXCkgnsOL6Dv3b9Kze3vdkveuMWl1va38K+U/t4Z/M7NKzUkJvb3ex1SEXGlwkmDsj4X1AYsDebNnEiEgRUBY7mtFFV3eP+TBCR93BOxc5LMKZk+OngT0z8diLJqclMHTiVXg16eR2SJ/4S+Rf2ndrHs2ufJbxqOH0a9vE6pCLhy1OktUBLEWkqIiHA9cCCTG0WAKPd59cCSzSHbo4iEiQitdznwcAwYEORR26KxVcxXzF20VgqBVfincvfKbPJBSAwIJB/9vknLau35K/L/krsiVivQyoSPksw7jWVCcAiYDPwgapuFJEnRGS422wmUFNEtgEPAOduZYtIDPAiMEZE4kSkHVAOWCQi64GfgD3AG776DsZ3Zm+czYNLH6RdzXa8M/QdmlVt5nVInqsQXIGX+79MgAQw8duJJCYneh1SodlgR1Os0jSN56OfZ86mOQxqMoh/9vknoUGhXoflV1buXcldX9/FiBYjeLL3k16Hkyc22NF4Ljk1mYeXPcycTXO4qe1NPN/3eUsuWejVoBd3dLqDT7Z9woLtma8qlCyWYEyxSExO5N4l9/JlzJfc3+V+/q/r/5X6SZgK4+6Iu4msG8k/Vv+jRF+Psd+w8bn4M/GM+3ocq/at4oleTzC249gydRu6IIICgnjqoqcICghi0opJpKSleB1SgViCMT51+PRhxi4ay6Yjm3ih7wtc1fIqr0MqMepVrMfk7pNZf2g9M36Z4XU4BWIJxvjM/lP7uXXhrcQmxPKfS/7DwCYDvQ6pxBnabCiXNb2M139+na1Ht3odTr5ZgjE+EZcQx5iFYzh0+hCvDXytTPdxKaxHuj1ClXJVeGzlYyXuVMkSjClysSdiuXXRrSScTWDG4Bl0qdvF65BKtGqh1Xik+yNsOrKJtzeVrE7rlmBMkYqJj+HWhbeSlJLEm0PeLDHTWvq7wU0Gc0njS5j601R2J+zOfQU/YQnGFJkd8Tu4bdFtpGgKM4fMpHWN1l6HVGqICA93e5hACeSpNU+VmInDLcGYIrEjfgdjF40lVVOZOXgmraq38jqkUqdexXqM7zye5XuWsyR2idfh5IklGFNoO+N3MnbRWFSVN4e8SYvqLbwOqdS6se2NtKreiqfXPl0iptq0BGMKZdeJXYxdNJY0TWPmkJk0r9bc65BKteCAYCZ1m8T+U/uZtWGW1+HkyhKMKbDdCbsZu2gsKWkpzBg8w5JLMYmsF8ngJoN5c8Ob7D+13+twcmQJxhTI3pN7GbtoLEmpSbwx+A1aVm/pdUhlygORD6AoL6570etQcmQJxuTbgVMHGLtoLCeTTzJ90HS7W+SBhpUacku7W/hy55dsPLzR63CyZQnG5Mvh04e5/avbOXbmGK8PfJ12Ndt5HVKZdVuH26herjovrnvRb29bW4IxeXY86Th3fHUHBxIPMPWSqXSs3dHrkMq0SiGVuDPiTr7f/z0r9qzwOpwslcTKjheKyC/uOq9kV7bEFK0TZ08w7utxxJ6I5d8D/m3d//3Eda2uI6xSGC/98BJpmuZ1OOfxWYLJUNnxMqAdcIM7r25G5yo7Av/CqewIv1d2/EsWm54GjANauo9Liz56k1FiciL3fHMPvx3/jX/1/xfd63f3OiTjCg4M5p4L7uHXY7/y1a6vvA7nPDkmGBGZleH56ByaZqXIKzuKSH2giqqucqsPvA1cmc+4TD6cST3DfUvuY/3h9Tx78bNcHHax1yGZTC4Lv4zmVZsz9aepflcdMrcjmIgMzyfmc9u+qOzY0N1OTtsErLJjUUhOS+bBqAf5fv/3/KP3P0p1GdeSLDAgkPGdx7Mzfidf7PzC63D+ILcEU5hL076o7Jjn9qo6XVUjVTWydu3aOWzSZCU1LZVHlj/C0rilTO4xmSuaX+F1SCYHA5sMpE2NNrz282t+NWdMbgkmzL2Q+u8Mz889clk3P5UdyWNlxzh3Ozlt0xRSmqbx+KrHWRizkAcvfJDrWl/ndUgmFwESwJ2d7iQ2IZaFMQu9Duec3ErHPpTheX4LC52r7IhTIO164MZMbdIrO64iD5UdVXWfiCSISA9gDXAL8O98xmVyoKo8t/Y5Pt72MXdF3MWYDmO8Dsnk0YDGA2hRrQVvrH+DoU2H+kXVhhwTjKrOzun9XNZNEZH0yo6BwJvplR2BaFVdgFPZcY5b2fEoThICzlV2rAKEiMiVwGBV3QTcDcwCygNfug9TRF796VXe2fwON7e9mfER470Ox+RDgAQwrtM4/rrsr3yz6xsGhw/2OqScKzuKSI5Vn1R1eE7v+wur7Jg3szbM4oV1L3B1y6uZ0nOKlRYpgVLTfOkCSwAAFmBJREFUUrnyf1dSPqg884bNK7bfYXaVHXM7ReqJc5fnfZxTEvsXV0p9+OuHvLDuBYaED+GxHo9ZcimhAgMCubXDrfx95d9ZvW81PRv09DSe3E7S6gGPAB2Al4FBwGFVXaqqS30dnCken+34jCdXPcnFYRfzVJ+nCAwI9DokUwjDmg2jdvnavLnhTa9DyTnBqGqqqi5U1dFAD2AbECUi9xZLdMbnlsQuYfKKyUTWi+SFvi8QHBjsdUimkEICQ7i53c2s3reaTUc2eRpLrpeZ3fFCVwPvAPcArwAf+Tow43sr967kL0v/Qvua7fn3gH9bIfpSZGSrkVQMruj5rHe5DRWYDawEugCPq2pXVX1SVfcUS3TGZ3448AP3f3s/Tas2ZerAqVQMruh1SKYIVQ6pzLUtr+WrXV95OutdbkcwfwJa4QwTWCUiJ9xHgoic8H14xhc2HtnIPYvvoW6Furw+6HWqlqvqdUjGB25seyOK8t6W9zyLIbdrMAGqWjnDo4r7qKyqVYorSFN0fjv2G3d9fRdVy1XljcFvUKt8La9DMj7SoFIDLml8CfN/nU/i/2/v3KOqrNI//nm4k4IppqGU0kSClqB5y0uJZlialZc5lr+EWVrLLqbjbdIZW07TWlP5059Y5mSlJDrpZGlqM2aaVup4HXREKCUlQEURFbnKAfbvj/eF0Lgp5+UcdH/WOuucs89+937efc551r5+H3uBU2yobYjkIyKTReRd8/BgbcvaGhcmNSeV5zY/h5ebFx8M+oDbm9zubJM0FjO241hyi3NZ/1ONW9oso7Yh0sdAN+Aw8Bgwz3KLNJaQkZvB+M3jUSg+iPqAO/zvqP0iTaMn/LZw7g24l7//8HenyGrW5mA6KqX+Ryn1PsZZoX4NYJPGwWTmZzJ+83gKSwpZMmgJdzW7y9kmaRoIEeHpsKc5kXOCPZl7Grz+2hyMvfyFqdeiaWRkFWQxfvN4ci7n6AgANylR7aNo7t2cVT+savC6axWcqrxyBHTWq0iNh3OF5xi3eRxnC86y+OHFdGrZydkmaZyAt7s3T4U8xbb0bZzOO92gdde2iuR+1cqRh15FahxcKLrAc5uf43Tead4b+B4RrSKcbZLGify2w29RSvHp0U8btF7nC0ZoHM7FoouM3zye9Nx03hn4Dt1u/9UhV81NRtumbXkw6EHWpqzFXmav/QIHoR3MDUbO5Rye//p5UnNSWRi5kF6BvZxtksZFGHXPKM4VnuPb9IY7p6wdzA1EzuUcntv8HD9d/InYAbH0btvb2SZpXIg+bfvQ+pbWrDm6psHq1A7mBqE86mK5c+nbtq+zTdK4GB5uHgwPGc6uU7vIyM2o/QIH4JKRHc3PZprpP4pIVKX0VDOy40ER0TJ1wPmi84zbPE47F02tDA8Zjojw+bGGEURwyciOZr7RQCeMyI3vmeWVE6mUiqhKou9m41zhOcZ9NY6fL/3MOwPf0c5FUyO3N7mdPm368EXKFw0SpM3KHsx1R3Y001cppS4rpU5gCF31sNDWRklmfia/2/Q7Tuad5L2B79G7jZ5z0dTOUyFPcbbwLLtO7bK8LisdTH0iO9Z0rQI2i8gBEXm+uspv9MiO6bnpxGyKIaswi/cHvU+PQO1/NXWjf1B/mns3Z13KOsvrstLB1CeyY03X9lFKdcUYer0kIlUGS76RIzsev3icmH/FkGfP48NHPqRLqy7ONknTiPB092TIXUPYlr6Ni0UXLa3LSgdTn8iO1V6rlCp/Pgus5SYbOiWeSyR6UzSlqpRlUcu4t+W9zjZJ0wh58u4nsZfZ+fLEl5bWY6WDqYjsKCJeGJO2V4tSlEd2hCsjO64HRpurTMFACLBXRJqIiB+AiDQBHgESLbwHl2LP6T2M+2ocTTybsPzR5YQ0D3G2SZpGSocWHQhrEWa5ToxlDsacUymP7JgM/KM8sqOIlAds+wgIMCM7TgFeNa89AvwDSAI2AS8ppUqB1sAOETkE7AW+VEq5TiBeC9mUuokXtrxAm6ZtWP7ocu70v9PZJmkaOY//5nGSspP46eJPltVRY2THG4XGHtlxZfJK3tr7Fl1adWHhgIVaQ1fjEM4VnuPhTx8mplMMk++fXK+yqovsqHfyujBlqoy5++by5t43ibwjUgt0axxKS9+W9G7Tm43HN1KmyiypQzsYF6WwpJBp305jedJyngl9hvn95+u4RRqH8/hvHudMwRn2Ze6zpHwt4u2CZBVk8co3r3Ak+wjTuk1jbMexDo8VbbfbycjIoKioyKHlahoX7VQ7YjvFUpZZRvLF5Frz+/j4EBQUhKdn3SKAagfjYhzJPsKkbyZxqfgSsZGxRN4ZaUk9GRkZ+Pn50b59ex3o/ianWW4zLhVfokOLDrhJ9YMapRTZ2dlkZGQQHBxcp7L1EMmF+PL4l0T/Kxo3cWP5o8stcy4ARUVFBAQEaOeioZl3M8pUGXnFeTXmExECAgKuqderezAugL3Mzvz981mRvIKurboyv/98AnwDLK9XOxcNQBPPJri7uZNTnIO/d81KuNf6m9EOxsmcyT/DjO9m8J+z/2FM2BimdpuKp1vdxrcajSMQEZp5NePC5QuUlpXi7uZe+0V1RA+RnMjOkzsZtWEUyeeTebPfm7za49Wbxrmkp6cTGRlJWFgYnTp1IjY21tL6Tp06xciRIwHYvn07Q4cOBSAuLo6XX375ussdPHgwt956a0V5daWyDVZQl/IPHjzIP//5T8AYJimlyC3Odagd2sE4geLSYubum8uELRMI8A1g1dBVDLlriLPNalA8PDyYN28eycnJ7N69m0WLFpGUlGRZfW3atGHNGsdLRU6fPp34+HiHl9sQVHYwvh6+eLp5klOc49A69BCpgfnx/I/M2jGLoxeOMrrDaKZ2m+r0/S1v7X2LH87/4NAyQ1uE8ocef6j288DAQAIDAwHw8/MjLCyMkydP0rHjL5pkOTk5hIeHc/z4cdzc3CgoKKBDhw4cP36cuLg4lixZQnFxMXfffTfx8fHccsstxMTE4O/vz/79+8nMzOTtt99m5MiRpKamMnToUBITqz+6tmHDBt544w2Ki4sJCAhg5cqVtG7dusb7HDhwINu3b68xT0pKChMmTCArKwt3d3c+/dQIHZKXl8fIkSNJTEzk/vvvZ8WKFYgIW7duZdq0aZSUlNC9e3cWL16Mt7c3+/btY9KkSeTn5+Pt7c3WrVvx9PTkhRdeYP/+/Xh4eDB//nwiI69cHNi7dy+TJ0+msLAQX19fli1bRnBwMK+99hqFhYXs2LGDmTNn0j2yO7+f+HvSjqZRWlLKnDlzeOKJqyWcrg3dg2kg7KV2/nbob4z+cjTZhdm8O+Bd/tjrj053Lq5AamoqCQkJ9OzZ84r0Zs2aER4ezrffGir4GzZsICoqCk9PT4YPH86+ffs4dOgQYWFhfPTRRxXXnT59mh07drBx40ZeffVXSq3V0rdvX3bv3k1CQgKjR4/m7bffdsj9jRkzhpdeeolDhw6xa9euCseakJDAggULSEpK4vjx4+zcuZOioiJiYmJYvXo1hw8fpqSkhMWLF1NcXIzNZiM2NpZDhw6xZcsWfH19WbRoEQCHDx/mk08+ITo6+lerPKGhoXz33XckJCTw+uuvM2vWLLy8vHj99dex2WwcPHgQm83G4nmL6dG3B19//zXbtm1j+vTp5Ofn1+vedQ+mATiUdYg///vPHLtwjMHtBzOr5yya+zR3tlkV1NTTsJq8vDxGjBjBggUL8Pf/9QqGzWZj9erVREZGsmrVKl588UUAEhMT+dOf/sTFixfJy8sjKqpCtpknn3wSNzc3OnbsyJkzZ+psS0ZGBjabjdOnT1NcXFznvR41kZuby8mTJ3nqqacAY6NaOT169CAoKAiAiIgIUlNT8fPzIzg4mHvuuQeA6OhoFi1axMCBAwkMDKR79+4AFW21Y8cOJk6cCBiOpF27dhw9evQKG3JycoiOjubYsWOICHZ71XGRtm3dxqX8SyxfvBwvNy+KiopIS0sjLCzsuu9fOxgLuVB0gdj/xPLZsc9odUsrFkYutHRvS2PDbrczYsQIxowZw/Dhw6vMM2zYMGbOnMn58+c5cOAAAwYMACAmJoZ169YRHh5OXFzcFcMUb2/vitfXcph34sSJTJkyhWHDhrF9+3bmzJlzXfdVmZrqr2ynu7s7JSUl1eZXSlW5RFyX+5s9ezaRkZGsXbuW1NRU+vfvX20dS/++lBZ3tOCeFvfg4VZ/96CHSBZQXFpMXGIcQz4fwhcpXxDTKYb1T67XzqUSSinGjRtHWFgYU6ZMqTZf06ZN6dGjB5MmTWLo0KG4uxtLqLm5uQQGBmK321m5cqVDbMrJyaFtW0OZ9eOPP65I37t3L2PHjr2uMv39/QkKCmLdOkOe8vLlyxQUFFSbPzQ0lNTUVFJSUgCIj4/noYceIjQ0lFOnTrFvn3FmKDc3l5KSEh588MGK+z969ChpaWl06NCh2vuKi4urSPfz8yM395dVo6ioKFZ8sIIyVUZucS4JCQnXdc+V0Q7GgdjL7Kw5uoYha4cw78A8urTuwmfDPmNqt6k08WzibPNcip07dxIfH88333xDREQEERERFSsaV2Oz2VixYgU2m60i7S9/+Qs9e/Zk0KBBhIaGOsSmOXPmMGrUKPr160fLli0r0tPS0vD19a3ymn79+jFq1Ci2bt1KUFAQX3311a/yxMfHs3DhQjp37kzv3r3JzMys1gYfHx+WLVvGqFGjuO+++3Bzc2PChAl4eXmxevVqJk6cSHh4OIMGDaKoqIgXX3yR0tJS7rvvPmw2G3FxcVf0jABmzJjBzJkz6dOnD6Wlv0QSiIyMJCkpiYiICFavXs3s2bMpKyljxEMj6NW1F7Nnz77WJvwVWg/GARTYC1ibspa4I3Fk5mfSuWVnJnad6NJhW5OTk+s1tr6ZmD59Os8++yydO3d2tikNQmZ+JueLztOheYcqN91V9dupTg9Gz8HUg5QLKXx27DPWpawjz55H11Zdmd1rNv3a9tPb8G8g5s6d62wTGhR/L3+yC7PJs+fVW3/IUgcjIoOBWMAd+FAp9eZVn3sDy4H7gWzAppRKNT+biRGYrRR4RSn1VV3KtJr0S+lsSdvCptRNJGUn4eHmwSPtHuHp0KeJaBXRkKZoNJbg6+GLh5sHl4ovua6DqRTZcRBGlIB9IrJeKVV5u2ZFZEcRGY0R2dF2VWTHNsAWEbnHvKa2Mh2GvcxO+qV0ErMTOXj2IHtO7yEtNw2AjgEdmdF9Bo8FP9YgBxOtoLqVCc3NjYjg5+VHzuUcylTZFRIO1zqlYmUPpiKyI4CIlEd2rOwMngDmmK/XAO9eHdkROGGKgpeHJ6mtzGvmVN4pliYu5XLpZYpKijhfdJ6zBWfJyM2gRJUAxonTbq278UzYM/S/oz9tm14dQ65x4ePjQ3Z2tpZs0FSJv5c/F4oukFecV3HCulwPpvJentqw0sFUFZ2xZ3V5lFIlIlI5suPuq64t/0fXViZgRHYEnge4886aFfjz7fl8/fPXeLl74ePuQ3Of5oQ0D+Hhdg8T3CyYji06Etws2KGnTJ1NUFAQGRkZ3IhRLzX1RylFVkEWBZ4FNPVqWpFermhXV6x0MFZEdqxqWb3KPptSagmwBIxVpOrNhJDmIXxr+7amLDccnp6eDtmpqrlxCSkNwdO9fqf7G1tkx7qUqdFoHEB9nQs0ssiOdSxTo9G4CJYNkcw5lfLIju7A0vLIjsB+pdR6jMiO8eYk7nkMh4GZrzyyYwm/RHakqjKtugeNRlM/boqdvCKSBfxch6wtgXMWm3O9aNuuH1e2z5Vtg7rb104pddvViTeFg6krIrK/qu3OroC27fpxZftc2Taov336sKNGo7EM7WA0Go1laAdzJUucbUANaNuuH1e2z5Vtg3rap+dgNBqNZegejEajsQztYDQajWVoB4OhMSMiP4pIiojUPc6FNbbcISLbRCRZRI6IyCQzvYWIfC0ix8xnp4YlEBF3EUkQkY3m+2AR2WPat9rcae0Mu24VkTUi8oPZhg+4UtuJyO/N7zVRRD4RER9ntZ2ILBWRsyKSWCmtyrYSg4Xmf+S/ItK1LnXc9A6mkm7No0BH4GlTj8ZZlABTlVJhQC/gJdOeV4GtSqkQYKv53plMApIrvX8L+D/TvgsYWj/OIBbYpJQKBcIxbHSJthORtsArQDel1L0Yu9HLdZCc0XZxwOCr0qprq0cxjuyEYKgULK5TDUqpm/oBPAB8Ven9TGCms+2qZM8XGAJbPwKBZlog8KMTbQoyf3wDgI0Yp9/PAR5VtWkD2uUPnMBcvKiU7hJtxy/yJC0wjulsBKKc2XZAeyCxtrYC3geeripfTY+bvgdD1bo1LqEmJSLtgS7AHqC1Uuo0gPncynmWsQCYAZSZ7wOAi0qZ6lzOa8O7gCxgmTl8+1BEmuAibaeUOgn8L5AGnAZygAO4RtuVU11bXdf/RDuYuunWNDgi0hT4DJislLrkbHvKEZGhwFml1IHKyVVkdUYbegBdgcVKqS5APs4fSlZgzmc8AQRjSME2wRh6XI3Tf39VcF3fsXYwLqgxIyKeGM5lpVLqczP5jIgEmp8HAmedZF4fYJiIpAKrMIZJC4BbTU0fcF4bZgAZSqk95vs1GA7HVdruYeCEUipLKWUHPgd64xptV051bXVd/xPtYFxMY8bUJP4ISFZKza/0UWXtnGiMuZkGRyk1UykVpJRqj9FW3yilxgDbMDR9nGafUioTSBeR8tCGAzEkP1yi7TCGRr1E5Bbzey63z+ltV4nq2mo9MNZcTeoF5JQPpWrEGZNdrvYAHgOOAj8Bf3SyLX0xup7/BQ6aj8cw5jm2AsfM5xYu0G79gY3m67swRMFSgE8BbyfZFAHsN9tvHdDcldoO+DPwA5AIxAPezmo74BOMuSA7Rg9lXHVthTFEWmT+Rw5jrITVWoc+KqDRaCxDD5E0Go1laAej0WgsQzsYjUZjGdrBaDQay9AORqPRWIaVkR01NyEiUoqxjClAKfCyUmpXPcvMU0o1reHzv2KEsrkVCFVKvVmf+jSOQ/dgNI6mUCkVoZQKxzg4+tcGqLMnxnmth4DvG6A+TR3RDkZjJf4Y8gPleiJzTR2UwyJiM9PfE5Fh5uu1IrLUfD1ORN6oqXCzvP8C3YF/A+OBxSLymoX3pLkG9BBJ42h8ReQg4INx3H+AmT4cY5dtOEYwr30i8h3wHdAPYyt6W/MaMHY0r6qpIqXUdBH5FHgWmAJsV0r1ceztaOqD7sFoHE35ECkUQ8xouXnupi/wiVKqVCl1BvgWo+fxPdDPFNVK4pfDdg8AdZm76YJxnCLUvF7jQugejMYylFL/FpGWwG1UfdwfpdRJU8ZgMEZvpgXwWyBPKZVbXdkiEoGhyBaEIdh0i5EsB4EHlFKFjrwXzfWhezAayxCRUAxZyGwM52EztXxvAx7EOOAHxvzJZDPP98A0apmsVUodVEpFYBxS7Qh8A0SZvSftXFwE3YPROJryORgwei3RSqlSEVmLMew5hHFafIYy5BXAcCaPKKVSRORnjF5MratBpqO6oJQqE5FQpZQeIrkY+jS1RqOxDD1E0mg0lqEdjEajsQztYDQajWVoB6PRaCxDOxiNRmMZ2sFoNBrL0A5Go9FYxv8DKaAreFpnFowAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.figure(figsize=(4, 6))\n", "plt.subplot(311)\n", "prior.plot(label='prior', color='gray')\n", "posterior1.plot(label='1 vanilla', color='C0')\n", "plt.ylabel('PMF')\n", "plt.title('101 Bowls')\n", "plt.legend()\n", "\n", "plt.subplot(312)\n", "posterior2.plot(label='2 vanilla', color='C1')\n", "plt.ylabel('PMF')\n", "plt.legend()\n", "\n", "plt.subplot(313)\n", "posterior3.plot(label='2 vanilla, 1 chocolate', color='C2')\n", "decorate_bowls('')\n", "\n", "savefig('fig02-01')" ] }, { "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 again:\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 again using a `Pmf`.\n", "I'll use integers to represent the hypotheses:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "hypos = [6, 8, 12]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And I can make the prior distribution like this:\n" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.333333
80.333333
120.333333
\n", "
" ], "text/plain": [ "6 0.333333\n", "8 0.333333\n", "12 0.333333\n", "dtype: float64" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prior = Pmf(1/3, hypos)\n", "prior" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As in the previous example, the prior probability gets broadcast across the hypotheses.\n", "\n", "Now we can compute the likelihood of the data:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "likelihood1 = 1/6, 1/8, 1/12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And use it to compute the posterior distribution.\n" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.444444
80.333333
120.222222
\n", "
" ], "text/plain": [ "6 0.444444\n", "8 0.333333\n", "12 0.222222\n", "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": {}, "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", "Now we can do the update in the usual way:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.000000
80.692308
120.307692
\n", "
" ], "text/plain": [ "6 0.000000\n", "8 0.692308\n", "12 0.307692\n", "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%.\n", "\n", "One note about the `Pmf` class: if you multiply a `Pmf` by a sequence, the result is a `Pmf`:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "empiricaldist.empiricaldist.Pmf" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(prior * likelihood1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you do it the other way around, the result is a `Series`:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "pandas.core.series.Series" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(likelihood1 * prior)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We usually want the posterior distribution to be a `Pmf`, so I usually put the prior first. But even if we do it the other way, we can use `Pmf` to convert the result to a `Pmf`:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.055556
80.041667
120.027778
\n", "
" ], "text/plain": [ "6 0.055556\n", "8 0.041667\n", "12 0.027778\n", "dtype: float64" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Pmf(likelihood1 * prior)" ] }, { "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": 45, "metadata": {}, "outputs": [], "source": [ "def update_dice(pmf, data):\n", " \"\"\"Update a pmf based on new data.\n", " \n", " pmf: Pmf of possible dice and their probabilities\n", " data: integer outcome\n", " \"\"\"\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 `qs` from the `Pmf`, which is the index of the `Series`; in this example, it represents the hypotheses.\n", "\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 die.\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": 46, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.333333
80.333333
120.333333
\n", "
" ], "text/plain": [ "6 0.333333\n", "8 0.333333\n", "12 0.333333\n", "dtype: float64" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pmf = prior.copy()\n", "pmf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And use `update_dice` to do the updates." ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.000000
80.692308
120.307692
\n", "
" ], "text/plain": [ "6 0.000000\n", "8 0.692308\n", "12 0.307692\n", "dtype: float64" ] }, "execution_count": 47, "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.\n" ] }, { "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", "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 as we see 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 \\...\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?" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
60.000000
80.835052
120.164948
\n", "
" ], "text/plain": [ "6 0.000000\n", "8 0.835052\n", "12 0.164948\n", "dtype: float64" ] }, "execution_count": 48, "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?" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
40.000000
60.000000
80.391304
120.347826
200.260870
\n", "
" ], "text/plain": [ "4 0.000000\n", "6 0.000000\n", "8 0.391304\n", "12 0.347826\n", "20 0.260870\n", "dtype: float64" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\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": 50, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
BlackWhite0.6
RedGreenBlue0.4
\n", "
" ], "text/plain": [ "BlackWhite 0.6\n", "RedGreenBlue 0.4\n", "dtype: float64" ] }, "execution_count": 50, "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": 51, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "0.30000000000000004" ] }, "execution_count": 51, "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": 52, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
priorlikelihoodunnormposterior
identical0.3333331.00.3333330.5
fraternal0.6666670.50.3333330.5
\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": 52, "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": 54, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
identical0.333333
fraternal0.666667
\n", "
" ], "text/plain": [ "identical 0.333333\n", "fraternal 0.666667\n", "dtype: float64" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "# Here's a solution using a Pmf\n", "\n", "hypos = ['identical', 'fraternal']\n", "prior = Pmf([1/3, 2/3], hypos)\n", "prior" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
probs
identical0.5
fraternal0.5
\n", "
" ], "text/plain": [ "identical 0.5\n", "fraternal 0.5\n", "dtype: float64" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "likelihood = 1, 1/2\n", "posterior = prior * likelihood\n", "posterior.normalize()\n", "posterior" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }