08_soccer_soln.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bite Size Bayes\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": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The World Cup problem\n",
    "\n",
    "In the 2018 FIFA World Cup final, France defeated Croatia 4 goals to 2.  Based on this outcome:\n",
    "\n",
    "1. How confident should we be that France is the better team?\n",
    "\n",
    "2. If the same teams played again, what is the chance Croatia would win?\n",
    "\n",
    "To answer these questions, we have to make some modeling decisions.\n",
    "\n",
    "* First, I'll assume that for any team against any other team there is some unknown goal-scoring rate, measured in goals per game, which I'll denote λ.\n",
    "\n",
    "* Second, I'll assume that a goal is equally likely during any minute of a game.  So, in a 90 minute game, the probability of scoring during any minute is $\\lambda/90$.\n",
    "\n",
    "* Third, I'll assume that a team never scores twice during the same minute.\n",
    "\n",
    "Of course, none of these assumptions is absolutely true in the real world, but I think they are reasonable simplifications.\n",
    "\n",
    "As George Box said, \"All models are false; some are useful.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goal-scoring\n",
    "\n",
    "Given these assumptions, we can simulate the goal-scoring process by tossing 90 imaginary coins, one for each minute of the game, with the same probability of yielding a goal.  Here's a function that does it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flip(n, p):\n",
    "    \"\"\"Flips n coins and returns the number of heads.\n",
    "    \n",
    "    n: number of coins\n",
    "    p: probability of heads\n",
    "    \"\"\"\n",
    "    return np.sum(np.random.random(n) < p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, suppose we know the goal scoring rate is 1.4; we can compute `p` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ = 1.4\n",
    "n = 90\n",
    "p = λ / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And call `flip` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip(n, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each time you call flip, you get a different number of goals.  Here's a list of 10 simulated scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 2, 2, 3, 1, 1, 2, 3, 0, 1]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[flip(n, p) for i in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we simulate many games, we expect the mean to be near 1.4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3919"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample = [flip(n, p) for i in range(10000)]\n",
    "\n",
    "np.mean(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sample contains scores from many simulated games.  The following function computes the number of times each score occurs and stores the result in a Series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_from_seq(seq):\n",
    "    \"\"\"Make a PMF from a sequence of values.\n",
    "    \n",
    "    seq: sequence\n",
    "    \n",
    "    returns: Series representing a PMF\n",
    "    \"\"\"\n",
    "    pmf = pd.Series(seq).value_counts(sort=False).sort_index()\n",
    "    pmf /= pmf.sum()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `pmf_from_seq` to compute the PMF of goals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.2453\n",
       "1    0.3470\n",
       "2    0.2446\n",
       "3    0.1143\n",
       "4    0.0356\n",
       "5    0.0097\n",
       "6    0.0028\n",
       "7    0.0005\n",
       "8    0.0002\n",
       "dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_goals = pmf_from_seq(sample)\n",
    "pmf_goals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PMF represents the distribution of goals scored, that is, each possible outcome and its probability.\n",
    "\n",
    "In this example, the PMF indicates that the probability of scoring 0 goals is about 27%, and the probability of scoring 1 goal is about 35%.\n",
    "\n",
    "Here's what the whole distribution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAETCAYAAAAs4pGmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfOUlEQVR4nO3de5hdVX3/8feHYEC5iZB6STIkQKgEqWiH2N+j5VIRggixihIURYqNWFFb1AqiCFEL4qVqxUKEeAExgtg62FhEuVR+iiQIgsFGQkQYgxgMF7kT+PSPvQKHkz0zZ8LsOUPyeT3PPLMva639nXOS8z17rb3Xlm0iIiLabdTtACIiYmxKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRT5mk0yV9ZITa6pF0r6RxZf0ySW8fibZLe9+XdPhItTeM435c0h2Sfj/Kx71Z0j6jecwmSDpR0jndjmNDs3G3A4ixTdLNwHOB1cCjwA3A14F5th8DsH3UMNp6u+0fDlTG9i3A5k8t6sePdyKwo+3DWtrffyTaHmYck4H3AdvZ/sNoHz9iXeUMIjpxoO0tgO2AU4APAmeN9EEkra9fWLYD/pjkUFmP3+f1ThJEdMz23bb7gEOAwyW9CEDSVyV9vCxvK+l7ku6StErSjyVtJOlsoAe4sHQh/bOkKZIs6UhJtwCXtGxr/RDZQdJVku6W9F1JzynH2ktSf2uMa7pUJM0EPgQcUo73i7L/8S6rEteHJf1W0h8kfV3SVmXfmjgOl3RL6R46fqDXRtJWpf7K0t6HS/v7ABcDLyhxfHWA+v8s6TZJKyS9vRx7x8HaLvt2kHSJpD+WGL8h6dkDHGOGpMWS7pF0u6TPDlCu9j0s+yZL+k6J5Y+SvjiM1/Lx97ls/ytJPynH+YWkvVpimCrpckl/knQxsO1Ar300Jwkihs32VUA/8Nc1u99X9k2g6pr6UFXFbwFuoTob2dz2qS119gR2BvYb4JBvBf4OeAFVV9cXOojxv4F/Ab5VjvfimmJvKz97A9tTdW19sa3MK4A/B14JnCBp5wEO+W/AVqWdPUvMR5TutP2BFSWOt7VXLMnsGGAfYMdSf8i211QHTqZ6bXYGJgMnDhDj54HP294S2AE4b4Byte+hqnGh7wG/BaYAE4EFpc7bGPq1fPx9ljQR+C/g48BzgPcDF0iaUMqeC1xNlRg+Boz6uFEkQcS6W0H1H7vdI8DzqfrbH7H9Yw894deJtu+z/cAA+8+2/Uvb9wEfAd5YPqyeqjcDn7W93Pa9wHHA7Lazl5NsP2D7F8AvgLUSTYnlEOA423+yfTPwGeAtHcbxRuArtpfYvh84qdO2bS+zfbHth2yvBD7L2glmjUeAHSVta/te21cOUq7uPZxBlYg+UN6vB21fUep08lq2vs+HAQttL7T9mO2LgcXAqyX1ALsDHyl/1/8AF3b4WsYISoKIdTURWFWz/VPAMuAHkpZLOraDtm4dxv7fAs9gZLocXlDaa217Y6pvzWu0XnV0P/UD6NsC42vamjiMOFr/xtblQduW9GeSFkj6naR7gHMY+LU5EtgJ+F9JiyS9ZoByA72Hk4Hf2l49wN8w1GvZ+ndtB7yhdC/dJekuqrO155e27ixfCFrbi1GWBBHDJml3qg+oK9r3lW+577O9PXAgcIykV67ZPUCTQ51hTG5Z7qH6hnsHcB/wrJa4xlF1i3Ta7gqqD6rWtlcDtw9Rr90dJab2tn7XYf3bgEkt661/71Btn0z1d/5F6To6jKrbaS22b7R9KPBnwCeBb0varKbcQO/hrUCP6geZO3ktW9+PW6nODJ/d8rOZ7VPK67F1W2w9dX9TNCsJIjomacvyrXMBcI7t62vKvEbSjpIE3EN1aeyjZfftVP3Tw3WYpOmSngXMBb5t+1Hg18Cmkg6Q9Azgw8AmLfVuB6asGWCt8U3gn8qA6OY8MWZR9w15QCWW84BPSNpC0nZUYwqdXrd/HnCEpJ3L33jCMNreArgXuKv0639goINIOkzShHJ58l1l86M15QZ6D6+i+vA+RdJmkjaV9PJSbbiv5TnAgZL2kzSutLWXpEm2f0vV3XSSpPGSXkGVqGKUJUFEJy6U9Ceqb33HU/VzHzFA2WnAD6k+tH4KfMn2ZWXfycCHS5fC+4dx/LOBr1J192wKvAeqq6qAfwDOpPpGfR/V4Ooa55fff5T085p255e2/wf4DfAg8O5hxNXq3eX4y6nOrM4t7Q/J9vepBt4vpera+WnZ9VAHbZ8EvBS4m2rQ9zuDHGomsETSvVQD1rNtP1hTrvY9LMnqQKqB9FuoXutDSp1hvZa2bwVmUQ2Ar6T6t/UBnvhMehPwMqpuzI9S3XsTo0x5YFDE2FKulPolsMlwz2YiRlLOICLGAEl/W7pTtqYaH7gwySG6LQkiYmx4B1VXy01U/f3v7G44EeliioiIAeQMIiIiaiVBRERErfVmVsVtt93WU6ZM6XYYERFPK1dfffUdtifU7VtvEsSUKVNYvHhxt8OIiHhakTTgNCbpYoqIiFpJEBERUSsJIiIiajWaICTNlLRU0rK6aZ8lHSXpeknXSrpC0vSyfYqkB8r2ayWd3mScERGxtsYGqcvUy6cBr6Ka1GuRpD7bN7QUO9f26aX8QVSTwM0s+26yvVtT8UVExOCaPIOYASwrT5h6mGqK6FmtBWzf07K6GUPP3x8REaOkyQQxkSc/QaqfmidsSXqXpJuAUynTOBdTJV1THlxe9+zjiIhoUJMJou6pVmudIdg+zfYOwAepHvgC1UNJemy/hOrhKOdK2nKtA0hzJC2WtHjlypUjGHpERDR5o1w/T3504iSqxxIOZAHw7wC2H6I8LMX21eUMYyeqp0w9zvY8YB5Ab2/v07d76tKTR66tvY8bubYiYoPW5BnEImBaeQTheGA20NdaQNK0ltUDgBvL9gllkBtJ21M94Wp5g7FGRESbxs4gbK+WdDRwETAOmG97iaS5wGLbfcDRkvaheij7ncDhpfoewFxJq6nmxj/K9qqmYo2IiLU1OheT7YXAwrZtrQ9kf+8A9S4ALmgytoiIGFzupI6IiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIio1WiCkDRT0lJJyyQdW7P/KEnXS7pW0hWSprfsO67UWyppvybjjIiItTWWICSNA04D9gemA4e2JoDiXNu72t4NOBX4bKk7HZgN7ALMBL5U2ouIiFHS5BnEDGCZ7eW2HwYWALNaC9i+p2V1M8BleRawwPZDtn8DLCvtRUTEKNm4wbYnAre2rPcDL2svJOldwDHAeOBvWupe2VZ3Yk3dOcAcgJ6enhEJOiIiKk2eQahmm9faYJ9mewfgg8CHh1l3nu1e270TJkx4SsFGRMSTNZkg+oHJLeuTgBWDlF8AvHYd60ZExAhrMkEsAqZJmippPNWgc19rAUnTWlYPAG4sy33AbEmbSJoKTAOuajDWiIho09gYhO3Vko4GLgLGAfNtL5E0F1hsuw84WtI+wCPAncDhpe4SSecBNwCrgXfZfrSpWCMiYm1NDlJjeyGwsG3bCS3L7x2k7ieATzQXXUREDCZ3UkdERK0kiIiIqJUEERERtZIgIiKiVhJERETUSoKIiIhaSRAREVErCSIiImolQURERK0kiIiIqJUEERERtZIgIiKiVhJERETUSoKIiIhaSRAREVErCSIiImolQURERK0kiIiIqJUEERERtZIgIiKiVqMJQtJMSUslLZN0bM3+YyTdIOk6ST+StF3LvkclXVt++pqMMyIi1rZxUw1LGgecBrwK6AcWSeqzfUNLsWuAXtv3S3oncCpwSNn3gO3dmoovIiIG1+QZxAxgme3lth8GFgCzWgvYvtT2/WX1SmBSg/FERMQwNJkgJgK3tqz3l20DORL4fsv6ppIWS7pS0mubCDAiIgbWWBcToJptri0oHQb0Anu2bO6xvULS9sAlkq63fVNbvTnAHICenp6RiToiIoBmzyD6gckt65OAFe2FJO0DHA8cZPuhNdttryi/lwOXAS9pr2t7nu1e270TJkwY2egjIjZwTSaIRcA0SVMljQdmA0+6GknSS4AzqJLDH1q2by1pk7K8LfByoHVwOyIiGtZYF5Pt1ZKOBi4CxgHzbS+RNBdYbLsP+BSwOXC+JIBbbB8E7AycIekxqiR2StvVT+vu0pNHpBn2Pm5k2hmrRup1gvX/tYpYTzU5BoHthcDCtm0ntCzvM0C9nwC7NhlbREQMLndSR0RErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRERG1OkoQki6QdICkJJSIiA1Epx/4/w68CbhR0imSXthgTBERMQZ0lCBs/9D2m4GXAjcDF0v6iaQjJD2jyQAjIqI7Ou4ykrQN8Dbg7cA1wOepEsbFjUQWERFd1dHzICR9B3ghcDZwoO3byq5vSVrcVHAREdE9nT4w6Mzy8J/HSdrE9kO2exuIKyIiuqzTLqaP12z76UgGEhERY8ugZxCSngdMBJ4p6SWAyq4tgWc1HFtERHTRUF1M+1ENTE8CPtuy/U/AhxqKKSIixoBBE4TtrwFfk/R62xeMUkwRETEGDDoGIemwsjhF0jHtP0M1LmmmpKWSlkk6tmb/MZJukHSdpB9J2q5l3+GSbiw/hw/7L4uIiKdkqC6mzcrvzYfbsKRxwGnAq4B+YJGkPts3tBS7Bui1fb+kdwKnAodIeg7wUaAXMHB1qXvncOOIiIh1M1QX0xnl90nr0PYMYJnt5QCSFgCzgMcThO1LW8pfCaw5Y9kPuNj2qlL3YmAm8M11iCMiItbBUFcxfWGw/bbfM8juicCtLev9wMsGKX8k8P1B6k4cLJaIiBhZQ3UxXf0U2lbNNtcWrMY6eoE9h1NX0hxgDkBPT8+6RRkREbU6uYppXfUDk1vWJwEr2gtJ2gc4HtjT9kMtdfdqq3tZTXzzgHkAvb29tcknIiLWzVBdTJ+z/Y+SLqTmG7ztgwapvgiYJmkq8DtgNtWU4a3tvwQ4A5hp+w8tuy4C/kXS1mV9X+C4of6YiIgYOUN1MZ1dfn96uA3bXi3paKoP+3HAfNtLJM0FFtvuAz5FdYXU+ZIAbrF9kO1Vkj5GlWQA5q4ZsI6IiNExVBfT1eX35ZLGU83oamCp7YeHarxM8LewbdsJLcv7DFJ3PjB/qGNEREQzOp3u+wDgdOAmqgHkqZLeYfv7g9eMiIinq06n+/4MsLftZQCSdgD+iycuS42IiPVMp9N9/2FNciiWA38YqHBERDz9DXUV0+vK4hJJC4HzqMYg3sATA8gREbEeGqqL6cCW5dt54ka2lcDWaxePiIj1xVBXMR0xWoFERMTY0ulVTJtSzZW0C7Dpmu22/66huCIioss6HaQ+G3ge1Syrl1NNffGnpoKKiIju6zRB7Gj7I8B9ZX6mA4BdmwsrIiK6rdME8Uj5fZekFwFbAVMaiSgiIsaETm+Um1cmzvsI0Ec1f9JHGosqIiK6rqMEYfvMsng5sH1z4URExFjRUReTpG0k/Zukn0u6WtLnJG3TdHAREdE9nY5BLKCaWuP1wMHAHcC3mgoqIiK6r9MxiOfY/ljL+sclvbaJgCIiYmzo9AziUkmzJW1Uft5INZtrRESsp4aarO9PVJPzCTgGOKfs2gi4F/hoo9FFRETXDDUX0xajFUhERIwtnY5BIOkgYI+yepnt7zUTUkREjAWdXuZ6CvBe4Iby896yLSIi1lOdnkG8GtjN9mMAkr4GXAMc21RgERHRXZ1exQTw7JblrTqpIGmmpKWSlklaK5lI2qPcfLda0sFt+x6VdG356RtGnBERMQI6PYM4GbhG0qVUVzTtARw3WAVJ44DTgFcB/cAiSX22b2gpdgvwNuD9NU08YHu3DuOLiIgRNmSCkCTgCuCvgN2pEsQHbf9+iKozgGW2l5d2FgCzqMYwALB9c9n32LoEHxERzRmyi8m2gf+0fZvtPtvf7SA5AEwEbm1Z7y/bOrWppMWSrsxd2xERo6/TLqYrJe1ue9Ew2lbNNg+jfo/tFZK2By6RdL3tm550AGkOMAegp6dnGE1HRMRQOh2k3psqSdwk6TpJ10u6bog6/cDklvVJwIpOA7O9ovxeDlwGvKSmzDzbvbZ7J0yY0GnTERHRgU7PIPZfh7YXAdMkTQV+B8wG3tRJxfJwovttPyRpW+DlwKnrEENERKyjoeZi2hQ4CtgRuB44y/bqThq2vVrS0cBFwDhgvu0lkuYCi233Sdod+A9ga+BASSfZ3gXYGTijDF5vBJzSdvVTREQ0bKgziK9RPY/6x1RnEdOp7qjuiO2FwMK2bSe0LC+i6npqr/cTYNdOjxMRESNvqAQx3fauAJLOAq5qPqSIiBgLhhqkfmTNQqddSxERsX4Y6gzixZLuKcsCnlnWRXWLxJaNRhcREV0z1PMgxo1WIBERMbYMZ7K+iIjYgCRBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRERG1On3kaET3XXryyLW193Ej11bEeipnEBERUSsJIiIiaiVBRERErSSIiIio1WiCkDRT0lJJyyQdW7N/D0k/l7Ra0sFt+w6XdGP5ObzJOCMiYm2NJQhJ44DTgP2B6cChkqa3FbsFeBtwblvd5wAfBV4GzAA+KmnrpmKNiIi1NXkGMQNYZnu57YeBBcCs1gK2b7Z9HfBYW939gIttr7J9J3AxMLPBWCMiok2TCWIicGvLen/Z1nTdiIgYAU0mCNVs80jWlTRH0mJJi1euXDms4CIiYnBNJoh+YHLL+iRgxUjWtT3Pdq/t3gkTJqxzoBERsbYmE8QiYJqkqZLGA7OBvg7rXgTsK2nrMji9b9kWERGjpLEEYXs1cDTVB/uvgPNsL5E0V9JBAJJ2l9QPvAE4Q9KSUncV8DGqJLMImFu2RUTEKGl0sj7bC4GFbdtOaFleRNV9VFd3PjC/yfgiImJguZM6IiJqJUFEREStJIiIiKiVBBEREbWSICIiolYSRERE1EqCiIiIWkkQERFRKwkiIiJqJUFEREStJIiIiKiVBBEREbWSICIiolYSRERE1EqCiIiIWkkQERFRKwkiIiJqJUFEREStJIiIiKiVBBEREbWSICIiolajCULSTElLJS2TdGzN/k0kfavs/5mkKWX7FEkPSLq2/JzeZJwREbG2jZtqWNI44DTgVUA/sEhSn+0bWoodCdxpe0dJs4FPAoeUfTfZ3q2p+CIiYnBNnkHMAJbZXm77YWABMKutzCzga2X528ArJanBmCIiokNNJoiJwK0t6/1lW20Z26uBu4Ftyr6pkq6RdLmkv24wzoiIqNFYFxNQdybgDsvcBvTY/qOkvwT+U9Iutu95UmVpDjAHoKenZwRCjoiINZo8g+gHJresTwJWDFRG0sbAVsAq2w/Z/iOA7auBm4Cd2g9ge57tXtu9EyZMaOBPiIjYcDWZIBYB0yRNlTQemA30tZXpAw4vywcDl9i2pAllkBtJ2wPTgOUNxhoREW0a62KyvVrS0cBFwDhgvu0lkuYCi233AWcBZ0taBqyiSiIAewBzJa0GHgWOsr2qqVgjImJtTY5BYHshsLBt2wktyw8Cb6ipdwFwQZOxRUTE4HIndURE1EqCiIiIWkkQERFRq9ExiIgNwqUnj0w7ex83Mu1EjJCcQURERK0kiIiIqJUEERERtZIgIiKiVhJERETUSoKIiIhaSRAREVErCSIiImolQURERK0kiIiIqJUEERERtZIgIiKiVhJERETUSoKIiIhaSRAREVErz4OIWB/lGRUxAnIGERERtRpNEJJmSloqaZmkY2v2byLpW2X/zyRNadl3XNm+VNJ+TcYZERFrayxBSBoHnAbsD0wHDpU0va3YkcCdtncE/hX4ZKk7HZgN7ALMBL5U2ouIiFHS5BjEDGCZ7eUAkhYAs4AbWsrMAk4sy98GvihJZfsC2w8Bv5G0rLT30wbjjYgmjdS4CGRsZJQ0mSAmAre2rPcDLxuojO3Vku4Gtinbr2yrO7H9AJLmAHPK6r2Slo5M6GwL3DF4kQ+N0KE61kFMMMpxjcWY4Gn7/iUmxua/qQ5jGnUjFdd2A+1oMkGoZps7LNNJXWzPA+YNP7TBSVpsu3ek230qElPnxmJciakzialzoxFXk4PU/cDklvVJwIqBykjaGNgKWNVh3YiIaFCTCWIRME3SVEnjqQad+9rK9AGHl+WDgUtsu2yfXa5ymgpMA65qMNaIiGjTWBdTGVM4GrgIGAfMt71E0lxgse0+4Czg7DIIvYoqiVDKnUc1oL0aeJftR5uKtcaId1uNgMTUubEYV2LqTGLqXONxqfrCHhER8WS5kzoiImolQURERK0kiIiIqJUEAUh6oaQPSvqCpM+X5Z27HddYU16nV0ravG37zC7GNEPS7mV5uqRjJL26W/HUkfT1bsfQTtIrymu1bxdjeJmkLcvyMyWdJOlCSZ+UtFWXYnqPpMlDlxw9ksZLequkfcr6myR9UdK7JD2j0WNv6IPUkj4IHAosoLr/Aqr7LmZTTfdxSrdiqyPpCNtf6cJx3wO8C/gVsBvwXtvfLft+bvulXYjpo1RzfW0MXEx1p/5lwD7ARbY/0YWY2i/lFrA3cAmA7YNGOyYASVfZnlGW/57qvfwPYF/gwm78O5e0BHhxueJxHnA/1ZQ7ryzbX9eFmO4G7gNuAr4JnG975WjH0RbTN6j+jT8LuAvYHPgO1esk24cPUv2psb1B/wC/Bp5Rs308cGO346uJ65YuHfd6YPOyPAVYTJUkAK7pYkzjyn+ce4Aty/ZnAtd1KaafA+cAewF7lt+3leU9u/jv5pqW5UXAhLK8GXB9l2L6Vevr1rbv2m69TlQ9K/tSXYa/Evhvqvu1tuhSTNeV3xsDtwPjyrqa/neeBwbBY8ALgN+2bX9+2TfqJF030C7guaMZS4txtu8FsH2zpL2Ab0vajvqpUUbDalf3x9wv6Sbb95T4HpDUlfcO6AXeCxwPfMD2tZIesH15l+JZYyNJW1N9+MnlW7Ht+ySt7lJMv2w5I/6FpF7biyXtBDzSpZhs+zHgB8APShfO/lS9DJ8GJnQhpo3KzcabUX0ZWjPjxCZAo11MSRDwj8CPJN3IE5ML9gA7Akd3KabnAvsBd7ZtF/CT0Q8HgN9L2s32tQC275X0GmA+sGuXYnpY0rNs3w/85ZqNpf+6KwmifLj8q6Tzy+/bGRv/z7YCrqb6N2RJz7P9+zKe1K0E/3bg85I+TDXp3E8l3Ur1//DtXYrpSa+F7UeoZnbok/TM7oTEWcD/Up0tHw+cL2k58FdUXeON2eDHIAAkbUQ1nfhEqn8g/cAij+7d263xnAV8xfYVNfvOtf2mLsQ0ieob++9r9r3c9v/vQkybuJoSvn37tsDzbV8/2jHVxHIA8HLboz4taickPQt4ru3fdDGGLYDtqRJpv+3buxjLTrZ/3a3jD0TSCwBsr5D0bKpxtltsNzoFURJERETUymWuERFRKwkiIiJqJUHEekGSJX2mZf39kk4coba/KungkWhriOO8QdKvJF3a4DFG5W+J9UMSRKwvHgJeVwaoxwxJ44ZR/EjgH2zv3VQ8EcORBBHri9VU8+P/U/uO9m/Nku4tv/eSdLmk8yT9WtIpkt4s6SpJ10vaoaWZfST9uJR7Tak/TtKnJC2SdJ2kd7S0e6mkc6lu5muP59DS/i8lfbJsOwF4BXC6pE+1ld9I0pckLZH0PUkL1/w9qqY+uaa0N1/SJmvaK3H9UtI8SWtdylr+3htK7J8e3ssdG4IkiFifnAa8eZjz+LyY6sa2XYG3ADu5mpLiTODdLeWmUN0NfQDVh/imVN/477a9O7A78PeqnoAI1WXTx9ue3nqwcrniJ4G/oZqyZHdJr7U9l+ru9Dfb/kBbjK8rx9+V6v6A/1fa2hT4KnCI7V2pLhN9Z6nzRdu7234R1Z3lr2mL4znA3wK72P4L4OOdvmCx4UiCiPVGuZP668B7hlFtke3byv0UN1HdQQvVN/8pLeXOs/2Y7RuB5cALqaZjeKuka4GfAdtQPR4X4KoB7i3YHbjM9krbq4FvAHsMEeMrqOYEeqzch7JmjOLPgd+0XLf/tZa29pb0M0nXUyWjXdravAd4EDhT0uuo5kGKeJIkiFjffI7qm/1mLdtWU/6tl66W8S37Wm+0e6xl/TGefAd0+w1Dprqp8t22dys/U22vSTD3DRDfuty1PFCd2u3lzOJLwMHlzOLLwKatZUpymgFcALyWar6hiCdJgoj1iu1VwHlUSWKNm3liKo5ZrNv8NW8oYwE7UN31u5TqeevvLPP1IGknSZsN1gjVmcaekrYtA9iHAkPN03QF8Ppy/OdSTQAI1fQLUyTtWNbfUtpakwzuKFNprHXVUtm+le2FVNPN7DZEDLEBGgtzxESMtM/w5Hm0vgx8V9JVwI8Y+Nv9YJZSffg+FzjK9oOSzqTqhvp5OTNZSfVtfEC2b5N0HFU3kYCFLtOmD+ICqqmdf0k1+/DPqMY+HpR0BNXcPBtTzdJ6uu2HJH2Zqpvs5rK93RZUr8mmJY61BvcjMtVGxNOApM3LBInbAFdRze+01rxYESMpZxARTw/fK5O0jQc+luQQoyFnEBERUSuD1BERUSsJIiIiaiVBRERErSSIiIiolQQRERG1kiAiIqLW/wGWk+vXdQgnuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_goals.plot.bar(color='C1', alpha=0.5)\n",
    "\n",
    "plt.xlabel('Number of goals')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Distribution of goals scored');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, if we know the goal-scoring rate, λ, we can use `flip` to simulate games and estimate the distribution of goals scored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Poisson distribution\n",
    "\n",
    "The distribution in the previous section is well approximated by the [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution).  If we know that the goal-scoring rate is λ, the probability of scoring $k$ goals is\n",
    "\n",
    "$f(k; \\lambda) = \\lambda^k exp(-\\lambda) ~/~ k!$\n",
    "\n",
    "SciPy provides a function that computes this expression; we can use it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.46596964e-01, 3.45235750e-01, 2.41665025e-01, 1.12777012e-01,\n",
       "       3.94719540e-02, 1.10521471e-02, 2.57883433e-03, 5.15766866e-04,\n",
       "       9.02592015e-05])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "\n",
    "xs = pmf_goals.index\n",
    "ys = poisson.pmf(xs, λ)\n",
    "ys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure compares these values with the results we got by simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAETCAYAAAAs4pGmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwV1Zn/8c8XZBG3KDCMAm2jwIiKkqRRM/pTHDfihmN0xESD2xAzITI4cUvUKDGj2XRwgkmIosRE0WiiHYfEYFyiBiLgCiiyBKSDIQguAYTY+vz+qGq8XKq7b7ddfZvm+369+tW1nHPqubfhPrfOqTqliMDMzKxYh3IHYGZmbZMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmVkmJwj7yCT9UNJVLdRWhaS1kjqm649LuqAl2k7b+7WkUS3VXhOOe52kNyT9pZWPu1TS0a15zDxIukbST8sdx7Zmu3IHYG2bpKVAL6AWeB+YD/wEmBQRHwBExIVNaOuCiHikvjIR8Rqw40eLetPxrgH6R8RZBe1/uiXabmIcfYH/AvaMiL+29vHNmstnEFaKkyJiJ2BP4AbgMuC2lj6IpPb6hWVPYLWTQ6Id/53bHScIK1lEvB0R1cAZwChJ+wNIukPSdelyD0kPSXpL0hpJT0rqIOlOoAL4VdqFdKmkSkkh6XxJrwGPFmwr/BDZW9Izkt6W9KCk3dJjDZNUUxhjXZeKpOHAV4Ez0uO9kO7f1GWVxnWlpGWS/irpJ5J2SffVxTFK0mtp99DX6ntvJO2S1l+Vtndl2v7RwHRgjzSOO+qpf6mk1yWtkHRBeuz+DbWd7ttb0qOSVqcx/kzSx+o5xkGSZkt6R9JKSTfWUy7zb5ju6yvpF2ksqyV9vwnv5aa/c7r9EEl/SI/zgqRhBTH0k/SEpL9Jmg70qO+9t/w4QViTRcQzQA3w/zJ2/1e6rydJ19RXkypxNvAaydnIjhHx7YI6RwCDgOPqOeTngfOAPUi6um4uIcbfAP8N3JMe78CMYuekP0cCe5F0bX2/qMxhwD8BRwFXSxpUzyH/F9glbeeINOZz0+60TwMr0jjOKa6YJrOLgaOB/mn9Rtuuqw5cT/LeDAL6AtfUE+MEYEJE7AzsDdxbT7nMv6GScaGHgGVAJdAbmJrWOYfG38tNf2dJvYH/A64DdgO+AtwvqWda9i5gDkli+AbQ6uNG5gRhzbeC5D92sfeA3Un629+LiCej8Qm/romIdRHxbj3774yIuRGxDrgK+Lf0w+qj+hxwY0QsiYi1wBXAyKKzl2sj4t2IeAF4Adgi0aSxnAFcERF/i4ilwPeAs0uM49+A2yNiXkSsB64tte2IWBQR0yNiY0SsAm5kywRT5z2gv6QeEbE2ImY2UC7rb3gQSSK6JP17bYiIp9I6pbyXhX/ns4BpETEtIj6IiOnAbOB4SRXAUOCq9HX9HvhVie+ltSAnCGuu3sCajO3fARYBv5W0RNLlJbS1vAn7lwGdaJkuhz3S9grb3o7kW3OdwquO1pM9gN4D6JzRVu8mxFH4GguXG2xb0j9Imirpz5LeAX5K/e/N+cBA4BVJsySdWE+5+v6GfYFlEVFbz2to7L0sfF17Aqen3UtvSXqL5Gxt97StN9MvBIXtWStzgrAmkzSU5APqqeJ96bfc/4qIvYCTgIslHVW3u54mGzvD6FuwXEHyDfcNYB3QrSCujiTdIqW2u4Lkg6qw7VpgZSP1ir2RxlTc1p9LrP860KdgvfD1Ntb29SSv84C06+gskm6nLUTEwog4E/gH4FvAfZJ2yChX399wOVCh7EHmUt7Lwr/HcpIzw48V/OwQETek78euRbFVZL0my5cThJVM0s7pt86pwE8j4qWMMidK6i9JwDskl8a+n+5eSdI/3VRnSdpXUjdgPHBfRLwPvAp0lXSCpE7AlUCXgnorgcq6AdYMdwPj0gHRHflwzCLrG3K90ljuBb4paSdJe5KMKZR63f69wLmSBqWv8eomtL0TsBZ4K+3Xv6S+g0g6S1LP9PLkt9LN72eUq+9v+AzJh/cNknaQ1FXSoWm1pr6XPwVOknScpI5pW8Mk9YmIZSTdTddK6izpMJJEZa3MCcJK8StJfyP51vc1kn7uc+spOwB4hORDawZwS0Q8nu67Hrgy7VL4ShOOfydwB0l3T1fgIkiuqgL+A7iV5Bv1OpLB1To/T3+vlvRsRruT07Z/D/wJ2AB8uQlxFfpyevwlJGdWd6XtNyoifk0y8P4YSdfOjHTXxhLavhb4BPA2yaDvLxo41HBgnqS1JAPWIyNiQ0a5zL9hmqxOIhlIf43kvT4jrdOk9zIilgMjSAbAV5H827qEDz+TPgscTNKN+XWSe2+slckPDDJrW9IrpeYCXZp6NmPWknwGYdYGSPrXtDtlV5LxgV85OVi5OUGYtQ1fIOlqWUzS3//F8oZj5i4mMzOrh88gzMwskxOEmZllajezKvbo0SMqKyvLHYaZ2VZlzpw5b0REz6x97SZBVFZWMnv27HKHYWa2VZFU7zQm7mIyM7NMThBmZpbJCcLMzDLlOgaRPghlAtARuDWdqbFw/4XAl0huDFoLjI6I+ZIqgZeBBWnRmaU+99jM2rf33nuPmpoaNmzImkbK6tO1a1f69OlDp06dSq6TW4JIp16eCBxDMqnXLEnVETG/oNhdEfHDtPzJJJPADU/3LY6IIXnFZ2Zbp5qaGnbaaScqKytJJpy1xkQEq1evpqamhn79+pVcL88upoOARekTpv5OMkX0iMICEfFOweoOND5/v5lt4zZs2ED37t2dHJpAEt27d2/yWVeeCaI3mz9BqoaMJ2xJ+pKkxcC3SadxTvWT9Fz64PKsZx+b2TbKyaHpmvOe5ZkgsqLZ4gwhIiZGxN7AZSQPfIHkoSQVEfFxkoej3CVp5y0OII2WNFvS7FWrVrVg6GZm9Vu6dCn777//FtsvuOAC5s+fn1Gj5Rx//PG89dZbDZaprKzkjTfe+MjHynOQuobNH53Yh+SxhPWZCvwAICI2kj4sJSLmpGcYA0meMrVJREwCJgFUVVVtvd1Tj13fcm0deUXLtWVmTXLrrbfmfoxp06blfow6eSaIWcAASf1InvY1kuQpUZtIGhARC9PVE4CF6faewJqIeF/SXiRPuFqSY6xtwk2LejFhca/GC2YYu/dKxvVv6qOUzay5amtrGTVqFM899xwDBw7kJz/5Cccffzzf/e53qaqqYscdd2Ts2LE89NBDbL/99jz44IP06tWLZcuWcd5557Fq1Sp69uzJ7bffTkVFBeeccw7bb789r7zyCsuWLeP2229nypQpzJgxg4MPPpg77rgD+HDWiB49enDKKaewfPlyNmzYwNixYxk9enSLvsbcEkRE1EoaAzxMcpnr5IiYJ2k8MDsiqoExko4meSj7m8CotPrhwHhJtSSXwF4YEWvyirWtGNffH/JmW4sFCxZw2223ceihh3Leeedxyy23bLZ/3bp1HHLIIXzzm9/k0ksv5cc//jFXXnklY8aM4fOf/zyjRo1i8uTJXHTRRTzwwAMAvPnmmzz66KNUV1dz0kkn8fTTT3PrrbcydOhQnn/+eYYM2fzCzsmTJ7Pbbrvx7rvvMnToUD7zmc/QvXv3FnuNud4HERHTgGlF2wofyD62nnr3A/fnGZuZtT83TX+VCb9b2HjBDGOPGsC4YwaWXL5v374ceuihAJx11lncfPPNm+3v3LkzJ554IgCf/OQnmT59OgAzZszgF79IHh1+9tlnc+mll26qc9JJJyGJwYMH06tXLwYPHgzAfvvtx9KlS7dIEDfffDO//OUvAVi+fDkLFy7cehKEmVlrGnfMwCZ9yH8UxVcFFa936tRp07aOHTtSW5v9BNnCel26dAGgQ4cOm5br1ovrP/744zzyyCPMmDGDbt26MWzYsBa/edBTbZiZNcNrr73GjBkzALj77rs57LDDSqr3z//8z0ydOhWAn/3sZyXXK/b222+z66670q1bN1555RVmzpzZrHYa4gRhZtYMgwYNYsqUKRxwwAGsWbOGL36xtMeI33zzzdx+++0ccMAB3HnnnUyYMKFZxx8+fDi1tbUccMABXHXVVRxyyCHNaqch7eaZ1FVVVbHVPg/Cl7malezll19m0KBB5Q5jq5T13kmaExFVWeV9BmFmZpmcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmbVhw4YNo1yX8HuqDTPburXkfUTge4kK+AzCzKyJ1q1bxwknnMCBBx7I/vvvzz333MP48eMZOnQo+++/P6NHj6buJuRhw4Yxbtw4Dj/8cAYNGsSsWbM49dRTGTBgAFdemTwjbenSpeyzzz6MGjWKAw44gNNOO43169dvcdzf/va3fOpTn+ITn/gEp59+OmvXrs31dTpBmJk10W9+8xv22GMPXnjhBebOncvw4cMZM2YMs2bNYu7cubz77rs89NBDm8p37tyZ3//+91x44YWMGDGCiRMnMnfuXO644w5Wr14NJNOHjx49mhdffJGdd955i+nD33jjDa677joeeeQRnn32WaqqqrjxxhtzfZ1OEGZmTTR48GAeeeQRLrvsMp588kl22WUXHnvsMQ4++GAGDx7Mo48+yrx58zaVP/nkkzfV22+//dh9993p0qULe+21F8uXLwe2nD78qaee2uyYM2fOZP78+Rx66KEMGTKEKVOmsGzZslxfp8cgzMyaaODAgcyZM4dp06ZxxRVXcOyxxzJx4kRmz55N3759ueaaazaberuUabwbmz48IjjmmGO4++6783pZW/AZhJlZE61YsYJu3bpx1lln8ZWvfIVnn30WgB49erB27Vruu+++JrfZ2PThhxxyCE8//TSLFi0CYP369bz66qsf8ZU0zGcQZmZN9NJLL3HJJZfQoUMHOnXqxA9+8AMeeOABBg8eTGVlJUOHDm1ym3XTh3/hC19gwIABW0wf3rNnT+644w7OPPNMNm7cCMB1113HwIH5PSDJ0323BZ7u26xk7XG676VLl3LiiScyd+7cXI/j6b7NzKxFOEGYmZVZZWVl7mcPzZFrgpA0XNICSYskXZ6x/0JJL0l6XtJTkvYt2HdFWm+BpOPyjNPMzLaUW4KQ1BGYCHwa2Bc4szABpO6KiMERMQT4NnBjWndfYCSwHzAcuCVtz8yM9jJ22pqa857leQZxELAoIpZExN+BqcCIwgIR8U7B6g5A3SsYAUyNiI0R8SdgUdqemW3junbtyurVq50kmiAiWL16NV27dm1SvTwvc+0NLC9YrwEOLi4k6UvAxUBn4F8K6s4sqts7o+5oYDRARUVFiwRtZm1bnz59qKmpYdWqVeUOZavStWtX+vTp06Q6eSYIZWzbIuVHxERgoqTPAlcCo5pQdxIwCZLLXD9StGa2VejUqRP9+vUrdxjbhDy7mGqAvgXrfYAVDZSfCpzSzLpmZtbC8kwQs4ABkvpJ6kwy6FxdWEDSgILVE4CF6XI1MFJSF0n9gAHAMznGamZmRXLrYoqIWkljgIeBjsDkiJgnaTwwOyKqgTGSjgbeA94k6V4iLXcvMB+oBb4UEe/nFauZmW0p17mYImIaMK1o29UFy2MbqPtN4Jv5RWdmZg3xndRmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy5RrgpA0XNICSYskXZ6x/2JJ8yW9KOl3kvYs2Pe+pOfTn+o84zQzsy1tl1fDkjoCE4FjgBpglqTqiJhfUOw5oCoi1kv6IvBt4Ix037sRMSSv+MzMrGF5nkEcBCyKiCUR8XdgKjCisEBEPBYR69PVmUCfHOMxM7MmyDNB9AaWF6zXpNvqcz7w64L1rpJmS5op6ZQ8AjQzs/rl1sUEKGNbZBaUzgKqgCMKNldExApJewGPSnopIhYX1RsNjAaoqKhomajNzAzI9wyiBuhbsN4HWFFcSNLRwNeAkyNiY932iFiR/l4CPA58vLhuREyKiKqIqOrZs2fLRm9mto3L8wxiFjBAUj/gz8BI4LOFBSR9HPgRMDwi/lqwfVdgfURslNQDOJRkALvF3DT9VSb8bmGz6o49agDjjhnYkuGYmbU5uSWIiKiVNAZ4GOgITI6IeZLGA7Mjohr4DrAj8HNJAK9FxMnAIOBHkj4gOcu5oejqp+Z77HoAxm0H445rbiMvAle0SDhtVvo+tYgj2/l7ZdZO5XkGQURMA6YVbbu6YPnoeur9ARicZ2xmZtawXBOEbf1uWtSLCYt7Navu2L1XMq7/yhaOyMxaixOENWhcf3/Im22rPBeTmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllKilBSLpf0gmSnFDMzLYRpX7g/4Bkqu6Fkm6QtE+OMZmZWRtQUoKIiEci4nPAJ4ClwHRJf5B0rqROeQZoZmblUXKXkaTuwDnABcBzwASShDE9l8jMzKysSprNVdIvgH2AO4GTIuL1dNc9kmbnFZyZmZVPqdN935o+/GcTSV0iYmNEVOUQl5mZlVmpXUzXZWyb0ZKBmJlZ29LgGYSkfwR6A9tL+jigdNfOQLecYzMzszJqrIvpOJKB6T7AjQXb/wZ8NaeYzMysDWgwQUTEFGCKpM9ExP2tFJOZmbUBDY5BSDorXayUdHHxT2ONSxouaYGkRZIuz9h/saT5kl6U9DtJexbsGyVpYfozqsmvzMzMPpLGuph2SH/v2NSGJXUEJgLHADXALEnVETG/oNhzQFVErJf0ReDbwBmSdgO+DlQBAcxJ677Z1DjMzKx5Guti+lH6+9pmtH0QsCgilgBImgqMADYliIh4rKD8TKDujOU4YHpErEnrTgeGA3c3Iw4zM2uGxq5iurmh/RFxUQO7ewPLC9ZrgIMbKH8+8OsG6vZuKBYzM2tZjXUxzfkIbStjW2QWTMY6qoAjmlJX0mhgNEBFRUXzojQzs0ylXMXUXDVA34L1PsCK4kKSjga+BhwRERsL6g4rqvt4RnyTgEkAVVVVmcnHzMyap7Eupv+JiP+U9CsyvsFHxMkNVJ8FDJDUD/gzMJJkyvDC9j8O/AgYHhF/Ldj1MPDfknZN148FrmjsxZiZWctprIvpzvT3d5vacETUShpD8mHfEZgcEfMkjQdmR0Q18B2SK6R+LgngtYg4OSLWSPoGSZIBGF83YG1mZq2jsS6mOenvJyR1JpnRNYAFEfH3xhpPJ/ibVrTt6oLloxuoOxmY3NgxzMwsH6VO930C8ENgMckAcj9JX4iIXzdc08zMtlalTvf9PeDIiFgEIGlv4P/48LJUMzNrZ0qd7vuvdckhtQT4a32Fzcxs69fYVUynpovzJE0D7iUZgzidDweQzcysHWqsi+mkguWVfHgj2ypg1y2Lm5lZe9HYVUzntlYgZmbWtpR6FVNXkrmS9gO61m2PiPNyisvMzMqs1EHqO4F/JJll9QmSqS/+lldQZmZWfqUmiP4RcRWwLp2f6QRgcH5hmZlZuZWaIN5Lf78laX9gF6Ayl4jMzKxNKPVGuUnpxHlXAdUk8yddlVtUZmZWdiUliIi4NV18Atgrv3DMzKytKKmLSVJ3Sf8r6VlJcyT9j6TueQdnZmblU+oYxFSSqTU+A5wGvAHck1dQZmZWfqWOQewWEd8oWL9O0il5BGRmZm1DqWcQj0kaKalD+vNvJLO5mplZO9XYZH1/I5mcT8DFwE/TXR2AtcDXc43OzMzKprG5mHZqrUDMzKxtKXUMAkknA4enq49HxEP5hGRmZm1BqZe53gCMBeanP2PTbWZm1k6VegZxPDAkIj4AkDQFeA64PK/AzMysvEq9igngYwXLu5RSQdJwSQskLZK0RTKRdHh6812tpNOK9r0v6fn0p7oJcZqZWQso9QzieuA5SY+RXNF0OHBFQxUkdQQmAscANcAsSdURMb+g2GvAOcBXMpp4NyKGlBifmZm1sEYThCQBTwGHAENJEsRlEfGXRqoeBCyKiCVpO1OBESRjGABExNJ03wfNCd7MzPLTaBdTRATwQES8HhHVEfFgCckBoDewvGC9Jt1Wqq6SZkua6bu2zcxaX6ldTDMlDY2IWU1oWxnbogn1KyJihaS9gEclvRQRizc7gDQaGA1QUVHRhKbNzKwxpQ5SH0mSJBZLelHSS5JebKRODdC3YL0PsKLUwCJiRfp7CfA48PGMMpMioioiqnr27Flq02ZmVoJSzyA+3Yy2ZwEDJPUD/gyMBD5bSsX04UTrI2KjpB7AocC3mxGDmZk1U2NzMXUFLgT6Ay8Bt0VEbSkNR0StpDHAw0BHYHJEzJM0HpgdEdWShgK/BHYFTpJ0bUTsBwwCfpQOXncAbii6+snMzHLW2BnEFJLnUT9JchaxL8kd1SWJiGnAtKJtVxcszyLpeiqu9wdgcKnHMTOzltdYgtg3IgYDSLoNeCb/kMzMrC1obJD6vbqFUruWzMysfWjsDOJASe+kywK2T9dFcovEzrlGZ2ZmZdPY8yA6tlYgZmbWtjRlsj4zM9uGOEGYmVkmJwgzM8vkBGFmZpmcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmVkmJwgzM8vkBGFmZplKfeSoWfk9dn3LtXXkFS3Xllk75QRhW6WbFvViwuJezao7du+VjDuyhQMya4ecIGyrNK7/Ssb1X1nuMMzaNY9BmJlZJicIMzPLlGuCkDRc0gJJiyRdnrH/cEnPSqqVdFrRvlGSFqY/o/KM08zMtpRbgpDUEZgIfBrYFzhT0r5FxV4DzgHuKqq7G/B14GDgIODrknbNK1YzM9tSnmcQBwGLImJJRPwdmAqMKCwQEUsj4kXgg6K6xwHTI2JNRLwJTAeG5xirmZkVyTNB9AaWF6zXpNvyrmtmZi0gzwShjG3RknUljZY0W9LsVatWNSk4MzNrWJ4JogboW7DeB1jRknUjYlJEVEVEVc+ePZsdqJmZbSnPBDELGCCpn6TOwEigusS6DwPHSto1HZw+Nt1mZmatJLcEERG1wBiSD/aXgXsjYp6k8ZJOBpA0VFINcDrwI0nz0rprgG+QJJlZwPh0m5mZtZJcp9qIiGnAtKJtVxcszyLpPsqqOxmYnGd8ZmZWP99JbWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpYp1wQhabikBZIWSbo8Y38XSfek+/8oqTLdXinpXUnPpz8/zDNOMzPb0nZ5NSypIzAROAaoAWZJqo6I+QXFzgfejIj+kkYC3wLOSPctjoghecVnZmYNy/MM4iBgUUQsiYi/A1OBEUVlRgBT0uX7gKMkKceYzMysRHkmiN7A8oL1mnRbZpmIqAXeBrqn+/pJek7SE5L+X45xmplZhty6mICsM4EosczrQEVErJb0SeABSftFxDubVZZGA6MBKioqWiBkMzOrk+cZRA3Qt2C9D7CivjKStgN2AdZExMaIWA0QEXOAxcDA4gNExKSIqIqIqp49e+bwEszMtl15JohZwABJ/SR1BkYC1UVlqoFR6fJpwKMREZJ6poPcSNoLGAAsyTFWMzMrklsXU0TUShoDPAx0BCZHxDxJ44HZEVEN3AbcKWkRsIYkiQAcDoyXVAu8D1wYEWvyitXMzLaU5xgEETENmFa07eqC5Q3A6Rn17gfuzzM2MzNrmO+kNjOzTLmeQZhtS26a/ioTfrewWXXHHjWAccdscR2GWVk5QZi1kHHHDPSHvLUrThBmH9Vj17dMO0de0TLtmLUQj0GYmVkmJwgzM8vkBGFmZpmcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmVkmJwgzM8vkBGFmZpmcIMzMLJMThJmZZXKCMDOzTJ7u26wd80OM7KNwgjBrj9JnVIzbDsYd19xGXgT8jIptmbuYzMwsU64JQtJwSQskLZJ0ecb+LpLuSff/UVJlwb4r0u0LJDX7O5CZmTVPbglCUkdgIvBpYF/gTEn7FhU7H3gzIvoDNwHfSuvuC4wE9gOGA7ek7ZmZWSvJcwziIGBRRCwBkDQVGAHMLygzArgmXb4P+L4kpdunRsRG4E+SFqXtzcgxXjPLUzouctOiXkxY3KtZTYzdeyXj+q/087tbiSIin4al04DhEXFBun42cHBEjCkoMzctU5OuLwYOJkkaMyPip+n224BfR8R9RccYDYxOV/8JWNBC4fcA3mihtlqKYypdW4zLMZWmLDF13KnHHh13+Njuzan7/rq3Xn//b2+saOmYoNXi2jMiembtyPMMQhnbirNRfWVKqUtETAImNT20hkmaHRFVLd3uR+GYStcW43JMpXFMpWuNuPIcpK4B+has9wGKs9mmMpK2A3YB1pRY18zMcpRngpgFDJDUT1JnkkHn6qIy1cCodPk04NFI+ryqgZHpVU79gAHAMznGamZmRXLrYoqIWkljgIeBjsDkiJgnaTwwOyKqgduAO9NB6DUkSYS03L0kA9q1wJci4v28Ys3Q4t1WLcAxla4txuWYSuOYSpd7XLkNUpuZ2dbNd1KbmVkmJwgzM8vkBGFmZpmcIABJ+0i6TNLNkiaky4PKHVdbk75PR0nasWj78DLGdJCkoenyvpIulnR8ueLJIukn5Y6hmKTD0vfq2DLGcLCkndPl7SVdK+lXkr4laZcyxXSRpL6Nl2w9kjpL+ryko9P1z0r6vqQvSeqU67G39UFqSZcBZwJTSe6/gOS+i5Ek033cUK7Yskg6NyJuL8NxLwK+BLwMDAHGRsSD6b5nI+ITZYjp6yRzfW0HTCe5C/9x4Gjg4Yj4ZhliKr6UW8CRwKMAEXFya8cEIOmZiDgoXf53kr/lL4FjgV+V49+5pHnAgekVj5OA9SRT7hyVbj+1DDG9DawDFgN3Az+PiFWtHUdRTD8j+TfeDXgL2BH4Bcn7pIgY1UD1jxaqkaMAAAWtSURBVCYitukf4FWgU8b2zsDCcseXEddrZTruS8CO6XIlMJskSQA8V8aYOqb/cd4Bdk63bw+8WKaYngV+CgwDjkh/v54uH1HGfzfPFSzPAnqmyzsAL5UpppcL37eifc+X630i6Vk5luQy/FXAb0ju19qpTDG9mP7eDlgJdEzXlfe/cz8wCD4A9gCWFW3fPd3X6iS9WN8uoHmznH10HSNiLUBELJU0DLhP0p5kT43SGmojuT9mvaTFEfFOGt+7ksrytwOqgLHA14BLIuJ5Se9GxBNliqdOB0m7knz4KdJvxRGxTlJtmWKaW3BG/IKkqoiYLWkg8F6ZYoqI+AD4LfDbtAvn0yS9DN8FMucsylmH9GbjHUi+DNXNONEFyLWLyQkC/hP4naSFwPJ0WwXQHxhTb6189QKOA94s2i7gD60fDgB/kTQkIp4HiIi1kk4EJgODyxTT3yV1i4j1wCfrNqb912VJEOmHy02Sfp7+Xknb+H+2CzCH5N9QSPrHiPhLOp5UrgR/ATBB0pUkE/TNkLSc5P/hBWWKabP3IiLeI5nZoVrS9uUJiduAV0jOlr8G/FzSEuAQkq7x3GzzYxAAkjqQTCfem+QfSA0wK1r37u3CeG4Dbo+IpzL23RURny1DTH1IvrH/JWPfoRHxdBli6hLJlPDF23sAu0fES60dU0YsJwCHRsRXyx1LFkndgF4R8acyxrATsBdJIq2JiJVljGVgRLxaruPXR9IeABGxQtLHSMbZXouIXKcgcoIwM7NMvszVzMwyOUGYmVkmJwhrFySFpO8VrH9F0jUt1PYd6RMScyXpdEkvS3osx2O0ymux9sEJwtqLjcCp6QB1myGpYxOKnw/8R0QcmVc8Zk3hBGHtRS3J/PjjincUf2uWtDb9PUzSE5LulfSqpBskfU7SM5JekrR3QTNHS3oyLXdiWr+jpO9ImiXpRUlfKGj3MUl3kdzMVxzPmWn7cyV9K912NXAY8ENJ3ykq30HSLZLmSXpI0rS616Nk6pPn0vYmS+pS114a11xJkyRtcSlr+nrnp7F/t2lvt20LnCCsPZkIfK6J8/gcSHJj22DgbGBgJFNS3Ap8uaBcJcnd0CeQfIh3JfnG/3ZEDAWGAv+u5AmIkFw2/bWI2LfwYOnlit8C/oVkypKhkk6JiPEkd6d/LiIuKYrx1PT4g0nuD/hU2lZX4A7gjIgYTHKZ6BfTOt+PiKERsT/JneUnFsWxG/CvwH4RcQBwXalvmG07nCCs3UjvpP4JcFETqs2KiNfT+ykWk9xBC8k3/8qCcvdGxAcRsRBYAuxDMh3D5yU9D/wR6E7yeFyAZ+q5t2Ao8HhErIqIWuBnwOGNxHgYyZxAH6T3odSNUfwT8KeC6/anFLR1pKQ/SnqJJBntV9TmO8AG4FZJp5LMg2S2GScIa2/+h+Sb/Q4F22pJ/62nXS2dC/YV3mj3QcH6B2x+B3TxDUNBclPllyNiSPrTLyLqEsy6euJrzl3L9dXJ3J6eWdwCnJaeWfwY6FpYJk1OBwH3A6eQzDdkthknCGtXImINcC9JkqizlA+n4hhB8+avOT0dC9ib5K7fBSTPW/9iOl8PkgZK2qGhRkjONI6Q1CMdwD4TaGyepqeAz6TH70UyASAk0y9USuqfrp+dtlWXDN5Ip9LY4qqldPsuETGNZLqZIY3EYNugtjBHjFlL+x6bz6P1Y+BBSc8Av6P+b/cNWUDy4dsLuDAiNki6laQb6tn0zGQVybfxekXE65KuIOkmEjAt0mnTG3A/ydTOc0lmH/4jydjHBknnkszNsx3JLK0/jIiNkn5M0k22NN1ebCeS96RrGscWg/tmnmrDbCsgacd0gsTuwDMk8zttMS+WWUvyGYTZ1uGhdJK2zsA3nBysNfgMwszMMnmQ2szMMjlBmJlZJicIMzPL5ARhZmaZnCDMzCyTE4SZmWX6/2smXrD/O+E2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmf_goals.plot.bar(color='C1', alpha=0.5, label='sample')\n",
    "plt.plot(xs, ys, '_', color='C0', markersize=15, label='binomial')\n",
    "\n",
    "plt.xlabel('Number of goals')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Distribution of goals scored')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results from the Poisson distribution are consistent with the results from simulation.\n",
    "\n",
    "In summary, given λ, we can compute the distribution of goals scored, either by simulation or using a Poisson distribution.\n",
    "\n",
    "Now let's turn it around: given the number of goals scored in a game, how can we estimate λ?\n",
    "\n",
    "First, we need a prior distribution for λ, which represents the range of possible values and their probabilities before we see the score.\n",
    "\n",
    "Then we use the data to update the distribution.  \n",
    "\n",
    "The result is a posterior distribution that represents our estimate of λ after seeing the score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The prior\n",
    "\n",
    "If you have ever seen a soccer game, you have some information about λ.  In most games, teams score a few goals each.  In rare cases, a team might score more than 5 goals, but they almost never score more than 10.\n",
    "\n",
    "I used [data from previous World Cups](https://www.statista.com/statistics/269031/goals-scored-per-game-at-the-fifa-world-cup-since-1930/) to estimate each team scores about 1.4 goals per game, on average.  So I'll set the mean of λ to be 1.4.\n",
    "\n",
    "For a good team against a bad one, we expect λ to be higher; for a bad team against a good one, we expect it to be lower.\n",
    "\n",
    "To model the distribution of goal-scoring rates, I will use a [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution), which I chose because:\n",
    "\n",
    "1. The goal scoring rate is a continuous quantity that cannot be less than 0, and the gamma distribution is appropriate for this kind of quantity.\n",
    "\n",
    "2. The gamma distribution has only one parameter, α, which is the mean.  So it's easy to construct a gamma distribution with the mean we want.\n",
    "\n",
    "3. As we'll see, the shape of the Gamma distribution is a reasonable choice, given what we know about soccer.\n",
    "\n",
    "We can use `scipy.stats.gamma` to compute the values of the gamma distribution, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gamma\n",
    "\n",
    "α = 1.4\n",
    "qs = np.linspace(0, 10, 101)\n",
    "ps = gamma.pdf(qs, α)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter, α, is the mean of the distribution.\n",
    "\n",
    "`qs` is an array of possible values for λ, from 0 to 10.\n",
    "\n",
    "`ps` is an array of probabilities, one for each value of λ.\n",
    "\n",
    "Now we can use `qs` and `ps` to make a `Pmf` that represents the prior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = pd.Series(ps, index=qs)\n",
    "prior /= prior.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I normalized the distribution so the probabilities add up to 1.\n",
    "\n",
    "We can use this function, from the previous notebook, to compute the prior mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_mean(pmf):\n",
    "    \"\"\"Compute the mean of a PMF.\n",
    "    \n",
    "    pmf: Series representing a PMF\n",
    "    \n",
    "    return: float\n",
    "    \"\"\"\n",
    "    return np.sum(pmf.index * pmf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4140818156118378"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_mean(prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior mean is a little higher than 1.4; that's because we are approximating a continuous distribution with an array of discrete values.\n",
    "\n",
    "But since the choice of prior is partly arbitrary, we don't have to be too precise.\n",
    "\n",
    "Here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxU9Z3v/9e7u6EbulmbVRpoDCA0iyzNEte4AlHUMRIxOoGoD5NJTCaTyb03k9/Emzi5mWTuTJaZxCQmGok34xI1igmRKBpBBWXfw6JsLTuINtjQdPfn90ed6hRFdXd101Wnu+rzfDzq0afO+qkSz6c+33PO9yszwznnnIuXE3YAzjnn2iZPEM455xLyBOGccy4hTxDOOecS8gThnHMuIU8QzjnnEvIE4TKOpD9KmpPC/ZukocH0zyR9o5X2O0jScUm5wfs/S7q7NfYd7C+l34vLPHlhB+BcUyTtBPoCtcAJYAHwRTM7nmh9M5uRrtjM7HPJrBd8hrvN7KVG9rUbKGqNuCR9ExhqZnfE7D9t34vLDF5BuPZippkVAROAScA/x6+giBb/m5YU2g+mMI/tXEM8Qbh2xczeBf4IjIb6Zpj/I+l14EPg/NimGUk5kv5Z0i5JByX9WlK3YFlp0Fx0l6TdwMuJjinpf0jaJ2mvpDvjlj0i6dvBdC9Jv5d0TNJRSUuC4z8KDAKeD5qQ/meiY8fMi00WH5H0lqT3JT0nqWdwrI9JqoiLZaekqyVNB74O3Bocb23Md9Wc72WOpN2SDkv6/1r8H821W54gXLsiaSDwcWB1zOy/Be4BugC74jaZG7yuAM4n0oTz47h1LgdGAtMSHG868FXgGmAYcHUj4f0jUAH0JtIk9nXAzOxvgd0EVZCZ/Vsyxw58GrgTOA+oAf6zkeND5IAvAN8BngiOd2GC1ebS9PdyCXABcBVwn6SRTR3bZRZPEK69eFbSMeA14FUiJ8CoR8xso5nVmNnpuO1uB75vZu8E1yz+CZgd9yv9m2Z2wsyqEhz3k8CvzGyDmZ0AvtlIjKeB/sBgMzttZkus6c7OGjs2wKMxx/4G8MnoRexzlMz38i0zqzKztcBaIFGicRnME4RrL24ys+5mNtjMPh93Qt3TyHbncWZVsYvIzRl9m7F97PL4CiXW/wW2A3+S9I6krzWybjLHjl++C+gA9Epiv01J5nvZHzP9Ia10Ad21H54gXCZo7Ff6XmBwzPtBRJpqDiS5/T5gYNz2iYMwqzSzfzSz84GZwFckXdXEMZqqMOKPfRo4TORurs7RBUFV0bsZ+03me3FZzhOEy3SPAf8gaYikIv7aNl+T5PZPAnMllUnqDPzvhlaUdL2koZIEfEDkttzaYPEBIm39zXVHzLHvB54ys1pgK1Ag6TpJHYjc1ZUfs90BoLSRu7rO9XtxWcAThMt0DwOPAouBHcBJ4IvJbmxmfwR+SOQOp+00cKdTYBjwEnAcWAo8YGZ/Dpb9K/DPwR1OX21G/I8CjxBp7ikAvhTE9T7weeCXwLtEKorYu5p+G/w9ImlVgv2e0/fisoN8wCDnnHOJeAXhnHMuIU8QzjnnEvIE4ZxzLiFPEM455xLKmA7CevXqZaWlpWGH4Zxz7crKlSsPm1nvRMsyJkGUlpayYsWKsMNwzrl2RVKDvQN4E5NzzrmEPEE455xLyBOEc865hDLmGoRzzrXE6dOnqaio4OTJk2GHklIFBQWUlJTQoUOHpLfxBOGcy2oVFRV06dKF0tJSIv0sZh4z48iRI1RUVDBkyJCkt/MmJudcVjt58iTFxcUZmxwAJFFcXNzsKskThHMu62VycohqyWf0BNFMZsbq1as5fPhw2KE451xK+TWIZqirq+P5559nzZo15ObmcvXVVzNlypSs+PXhnAvffffdx2WXXcbVV1+dluN5gkhSXV0dzz77LOvXr+eiiy7i8OHDLFy4kK1bt3LzzTdTVOTD9TrnUqe2tpb777+/2dvk5ua2+JjexJQEM+OZZ55h/fr1XHnllVxzzTXMnj2b66+/nj179vDcc8/hAy8551pq586djBgxgjlz5jB27FhuueUWPvzwQ0pLS7n//vu55JJL+O1vf8vcuXN56qmnAFi0aBHjx49nzJgx3HnnnZw6dQrgrG3OhVcQSTh48CAbN27ksssu49JLLwUiF3wmTpxIdXU1f/rTn9i2bRvDhw8POVLn3Ll44YUX2L9/f6vus1+/fkyfPr3J9bZs2cJDDz3ExRdfzJ133skDDzwARJ5feO211+rjg8idV3PnzmXRokUMHz6cT3/60/z0pz/ly1/+8lnbnAuvIJKwY8cOACZMmHDWssmTJ1NcXMzChQupqfHx3p1zLTNw4EAuvvhiAO644476E/ytt9561rpbtmxhyJAh9T9K58yZw+LFi+uXJ9qmJbyCSMKOHTvo2bMn3bp1O2tZbm4u06dP5ze/+Q1vvvlm/X9g51z7k8wv/VSJv9kl+r6wsPCsdZtq0k60TUt4BdGEuro6du7c2ejTh0OHDmX48OEsXryYysrKNEbnnMsUu3fvZunSpQA89thjXHLJJQ2uO2LECHbu3Mn27dsBePTRR7n88stbPSZPEE3Yu3cv1dXVTT6ePm3aNGpra88o85xzLlkjR45k3rx5jB07lqNHj/J3f/d3Da5bUFDAr371K2bNmsWYMWPIycnhc5/7XKvHlNImJknTgR8BucAvzey7ccvzgV8DE4EjwK1mtlNSKbAZ2BKsuszMWv/TJyF6/aGp0ep69uzJmDFjWLt2LVdddRUFBQVpiM45lylycnL42c9+dsa8nTt3nvH+kUceqZ++6qqrWL169Vn7id/mnGJqtT3FkZQL/ASYAZQBt0kqi1vtLuA9MxsK/AD4Xsyyt81sXPAKJTlAJEH07ds3qTa9yZMnc/r0adasWZOGyJxzLrVS2cQ0GdhuZu+YWTXwOHBj3Do3AvOC6aeAq9SGHkuuqalhz549Sfd+2L9/fwYOHMjy5cv9uQjnXNJKS0vZsGFD2GGcJZUJYgCwJ+Z9RTAv4TpmVgO8DxQHy4ZIWi3pVUmXJjqApHskrZC04tChQ60bPbBnzx5qamqa1T3upEmTOHr0aP3FI+dc25cNP+ha8hlTmSASVQLxETa0zj5gkJmNB74C/LekrmetaPagmZWbWXnv3r3POeB4O3bsQBKDBw9OepuysjKKiop46623Wj0e51zrKygo4MiRIxmdJKLjQTT32mgqL1JXAANj3pcAextYp0JSHtANOGqR/1KnAMxspaS3geHAihTGe5adO3cyYMAA8vPzk94mNzeXiRMn8uqrr3LkyBGKi4ub3sg5F5qSkhIqKipIRStEWxIdUa45UpkglgPDJA0B3gVmA5+KW2c+MAdYCtwCvGxmJqk3kURRK+l8YBjwTgpjPUt1dTXvvvsuF110UbO3LS8vZ8mSJaxYsYJp06alIDrnXGvp0KFDs5qRs0nKmpiCawr3AguJ3LL6pJltlHS/pBuC1R4CiiVtJ9KU9LVg/mXAOklriVy8/pyZHU1VrIkcPXqUuro6+vfv3+xti4qKuOCCC1i/fj11dXUpiM4551Ivpc9BmNkCYEHcvPtipk8CsxJs9zTwdCpja0r0ieguXbq0aPuxY8eyefNm3n77bYYNG9aaoTnnXFr4k9QNONcEMWzYMDp16sTatWtbMyznnEsbTxANiCaIlg4ElJuby+jRo9myZUuzBwp3zrm2wBNEAyorK+ncuTN5eS1vhRs7diw1NTVs2rSpFSNzzrn08ATRgMrKyhY3L0UNGDCAnj17sm7dulaKyjnn0scTRANaI0FI4sILL2TXrl0cO3aslSJzzrn08ATRgNZIEBBpZgK8inDOtTueIBKoq6vj+PHjrZIgunfvzqBBg9pkR1zOOdcYTxAJHD9+HGj5La7xRo0axaFDhzh48GCr7M8559LBE0QC5/oMRLyRI0cCsHHjxlbZn3POpYMniARaO0F06dKF0tJSNm3alNE9RjrnMosniARaO0FApBvww4cPezOTc67d8ASRQGVlJZKSGmY0WWVlZUjyZibnXLvhCSKByspKioqKyMlpva+nsLCQ0tJSNm7c6M1Mzrl2wRNEAq31DES8UaNGcfToUfbv39/q+3bOudbmCSKBVCWIkSNHejOTc67d8ASRQKoSROfOnTn//PP9bibnXLvgCSJOTU0NVVVVKUkQEKki3nvvPQ4cOJCS/TvnXGvxBBEnFbe4xhoxYgSSvAtw51yb5wkiTqoTRGFhIYMHD2bz5s0p2b9zzrUWTxBxUp0g4K8PzR06dChlx3DOuXPlCSJOOhLEiBEjALyZyTnXpnmCiFNZWUlubi6dOnVK2TG6dOnCoEGDPEE459o0TxBxore4SkrpcUaOHMnBgwc5cuRISo/jnHMt5QkiTqqegYgX7QLcqwjnXFvlCSJOuhJEt27dGDBggN/N5JxrszxBxElXgoBIFbFv3z6OHTuWluM551xzeIKIcerUKaqrq9OWIMrKygC8inDOtUmeIGJEx6IuKipKy/F69OhBv379PEE459qklCYISdMlbZG0XdLXEizPl/REsPxNSaVxywdJOi7pq6mMM+r06dMAdOzYMR2HAyLNTHv27Kl//sI559qKlCUISbnAT4AZQBlwm6SyuNXuAt4zs6HAD4DvxS3/AfDHVMUYr7a2FqBVBwpqSvRuJq8inHNtTSrPhJOB7Wb2jplVA48DN8atcyMwL5h+CrhKwQMIkm4C3gHSNnhCXV0dALm5uek6JL1796ZXr16eIJxzbU4qE8QAYE/M+4pgXsJ1zKwGeB8ollQI/C/gWymM7yxhVBAQqSJ27drFiRMn0npc55xrTCrPhIkeRY4fJaehdb4F/MDMjjd6AOkeSSskrWiNju/CqCAgkiDMjC1btqT1uM4515hUJogKYGDM+xJgb0PrSMoDugFHgSnAv0naCXwZ+Lqke+MPYGYPmlm5mZX37t37nAOOVhDpThD9+vWje/fu3szknGtTUpkglgPDJA2R1BGYDcyPW2c+MCeYvgV42SIuNbNSMysFfgh8x8x+nMJYgb9WEOluYpJEWVkZ77zzDlVVVWk9tnPONSRlZ8LgmsK9wEJgM/CkmW2UdL+kG4LVHiJyzWE78BXgrFth0ymsCgIiD83V1dWxdevWtB/bOecSyUvlzs1sAbAgbt59MdMngVlN7OObKQkugbAqCIDzzjuPrl27smnTJi688MK0H9855+L5k9Qxwqwgos1Mb7/9NqdOnUr78Z1zLp4niBhhVhAQaWaqra31ZibnXJvgCSJGmBUEQElJCV26dPExIpxzbYIniBhhVxCSGDlyJNu3b6e6ujqUGJxzLsoTRIywKwiINDPV1NSwbdu20GJwzjnwBHGGsLraiDVw4EAKCwu9mck5FzpPEDHC6mojVk5ODiNHjmTbtm3ezOScC5UniBhtoYIAGDVqFKdPn/ZmJudcqDxBxKirqyMnJ4egx/HQDBo0yJuZnHOh8wQRo7a2NvTqASIVTFlZGVu3bvVmJudcaMI/G7YhdXV1oV5/iDVq1Chqamr8oTnnXGg8QcRoKxUERO5mKioq8mYm51xo2sbZsI1oSxVEtJlp27Zt3jeTcy4UniBiRC9StxXezOScC1PbORu2AbW1tW2mgoBIM1OXLl3YuHFj2KE457KQJ4gYba2CiHYBvn37dk6ePBl2OM65LNN2zoZtQFurIADGjBlDbW2tj1ftnEs7TxAx2tJdTFHnnXcePXr0YMOGDWGH4pzLMm3rbBiytnQXU5QkRo8ezY4dOzh+/HjY4TjnsogniBhtsYIAGD16NGbmz0Q459Kq7Z0NQ9QWKwiAPn360KdPH29mcs6lVVIJQtLTkq6TlNEJpS1epI4aPXo0e/bs4dixY2GH4pzLEsme8H8KfArYJum7kkakMKbQtLXbXGONHj0awKsI51zaJHU2NLOXzOx2YAKwE3hR0huSPiOpQyoDTKe2XEH06NGDkpISTxDOubRJ+ueypGJgLnA3sBr4EZGE8WJKIgtBW64gIPJMxIEDBzh48GDYoTjnskCy1yCeAZYAnYGZZnaDmT1hZl8EilIZYDq15QoCIn0z5eTksHbt2rBDcc5lgWR/Lv/SzMrM7F/NbB+ApHwAMytPWXRp1tYriMLCQoYOHcr69evrx892zrlUSfZs+O0E85a2ZiBtQVt9DiLW2LFjqaysZOfOnWGH4pzLcHmNLZTUDxgAdJI0HogO1tyVSHNTRmnrTUwAw4cPJz8/n3Xr1nH++eeHHY5zLoM19XN5GvDvQAnwfeA/gtdXgK83tXNJ0yVtkbRd0tcSLM+X9ESw/E1JpcH8yZLWBK+1kv6meR+rZdp6ExNAhw4dKCsrY/PmzZw+fTrscJxzGazRs6GZzTOzK4C5ZnZFzOsGM3umsW0l5QI/AWYAZcBtksriVrsLeM/MhgI/AL4XzN8AlJvZOGA68HNJjVY7raE9VBAQaWaqrq7mL3/5S9ihOOcyWFNNTHeY2f8DSiV9JX65mX2/kc0nA9vN7J1gX48DNwKxHQrdCHwzmH4K+LEkmdmHMesUANbUB2kN7aGCABg8eDDdunVj3bp1jBkzJuxwnHMZqqmzYWHwtwjokuDVmAHAnpj3FcG8hOuYWQ3wPlAMIGmKpI3AeuBzwfIzSLpH0gpJKw4dOtREOI0zszbbF1M8SYwdO5a3336bDz74IOxwnHMZqtEKwsx+Hvz9Vgv2rQTz4iuBBtcxszeBUZJGAvMk/dHMzhhWzcweBB4EKC8vP6cqI3rbaHtIEAAXXnghS5YsYd26dVxyySVhh+Ocy0BNNTH9Z2PLzexLjSyuAAbGvC8B9jawTkVwjaEbcDTuGJslnQBGAysai+dcRBNEe2hiAiguLmbQoEGsXr2aiy++GClRrnXOuZZr6my4solXY5YDwyQNkdQRmA3Mj1tnPjAnmL4FeNnMLNgmD0DSYOACIn1ApUxtbS3QfioIgPHjx3P06FH27NnT9MrOOddMTTUxzWvpjs2sRtK9wEIgF3jYzDZKuh9YYWbzgYeARyVtJ1I5zA42vwT4mqTTQB3weTM73NJYktHeKgiAsrIyFixYwOrVqxk0aFDY4TjnMkxTTUw/NLMvS3qeBHcSmdkNjW1vZguABXHz7ouZPgnMSrDdo8CjjYfeutpjBdGxY0dGjRrFxo0bmTFjBh07dgw7JOdcBmnq2YLoSfrfUx1I2NpjBQGRZqY1a9awadMmxo0bF3Y4zrkM0tSDciuDv68S6XvpPSJNQUuDeRmjPVYQAAMHDqRnz56sXr067FCccxkm2e6+rwPeBv4T+DGwXdKMVAaWbtEE0d4qCEmMHz+e3bt3c/hwSi/TOOeyTLJnw/8ArjCzj5nZ5cAVRLrGyBjt7TmIWOPGjSMnJ4eVK5u6scw555KXbII4aGbbY96/A2TUsGbttYIAKCoqYsSIEaxdu5aamrMeOHfOuRZp9Gwo6WZJNwMbJS2QNFfSHOB5Is85ZIz2XEEATJw4kaqqKjZt2tT0ys45l4Sm7mKaGTN9ALg8mD4E9EhJRCFprxepo4YMGUKPHj1YuXIlY8eODTsc51wGaOpBuc+kK5CwtdfbXKMkMXHiRF566SUOHTpE7969ww7JOdfOJXsXU4GkL0h6QNLD0Veqg0un9l5BQORidW5url+sds61imR/Lj8K9CMywtyrRDreq0xVUGFo7xUEQGFhISNHjmTt2rU+2pxz7pwlezYcambfAE4E/TNdB2TUSDWZUEEAlJeXc/LkSdavXx92KM65di7ZBBH9OXpM0mgi3XKXpiSikGRCBQEwaNAg+vbty1tvvYVZWgbic85lqGTPhg9K6gF8g0gX3Zv46/jRGSFTKghJTJ48mQMHDrB79+6ww3HOtWNJJQgz+6WZvWdmr5rZ+WbWJzraXKbIlAoCYMyYMRQUFPDmm2+GHYpzrh1L9i6mYkn/JWmVpJWSfiipONXBpVOmVBAAHTp0YMKECfzlL3/h/fffDzsc51w7lezP5ceJdK3xCSIjvx0GnkhVUGFoz11tJDJp0iQAVqxI2SitzrkMl+zZsKeZ/YuZ7Qhe3wa6pzKwdGvvXW3E6969OxdccAErV670W16dcy2SbIJ4RdJsSTnB65PAH1IZWLplUhNT1JQpU6iqqmLdunVhh+Kca4ea6qyvUtIHwGeB/waqg9fjwD+kPrz0yaSL1FGDBw+mf//+LF261G95dc41W1MjynUxs67B3xwzywteOWbWNV1BpkMmVhCSuOiiizhy5AhbtmwJOxznXDuT9M9lSTdI+vfgdX0qgwpDXV0dkpAUdiitqqysjG7durF06dKwQ3HOtTPJ3ub6XeDviTwgtwn4+2Bexqitrc2o5qWonJwcpk6dyu7du6moqAg7HOdcO5LsGfHjwDVm9rCZPQxMD+ZljLq6uoxqXoo1fvx4CgoKvIpwzjVLc34yx97W2q21AwlbplYQAPn5+UycOJHNmzdz5MiRsMNxzrUTyZ4R/xVYLekRSfOAlcB3UhdW+mVyBQEwdepUcnNzef3118MOxTnXTjSZIBS5avsaMBV4Jnh91MweT3FsaZXJFQRAUVEREyZMYO3atd79hnMuKU2eES1yA/2zZrbPzOab2XNmtj8NsaVVplcQABdddBGAVxHOuaQk+5N5maRJKY0kZJleQQB069aNcePGsWrVKiorM2pAQOdcCiR7RryCSJJ4W9I6SeslNdl/g6TpkrZI2i7pawmW50t6Ilj+pqTSYP41Qa+x64O/VzbnQ7VENlQQAJdccgl1dXV+R5Nzrkl5Sa43o7k7lpQL/AS4BqgAlkuab2abYla7C3jPzIZKmk1kEKJbifQWO9PM9gYj2C0EBjQ3huaora3NigTRo0cPxowZw4oVK7j44ospLCwMOyTnXBvVVF9MBZK+DPwPIs8+vGtmu6KvJvY9GdhuZu+YWbT/phvj1rkRmBdMPwVcJUlmttrM9gbzNwIFkvKb8bmara6uLuObmKIuvfRSampq/FqEc65RTZ0R5wHlwHoiVcR/NGPfA4A9Me8rOLsKqF/HzGqA94H4gYg+Aaw2s1PxB5B0j6QVklYcOnSoGaGdLVsqCIBevXoxduxYli9f7tcinHMNaipBlJnZHcHworcAlzZj34k6NYrvUrTRdSSNItLs9NlEBzCzB82s3MzKe/fu3YzQzpZNFQTA5ZdfTl1dHYsXLw47FOdcG9XUGbF+pJngF35zVAADY96XAHsbWkdSHpEntI8G70uA3wGfNrO3m3nsZsumCgIi1yLGjx/PqlWreO+998IOxznXBjWVIC6U9EHwqgTGRqeDcSIasxwYJmmIpI7AbGB+3DrzgTnB9C3Ay2ZmkroTGZDon8wsLQ3l2VZBAFx22WVI8irCOZdQU+NB5AbjQUTHhMiLmW50PIig4riXyB1Im4EnzWyjpPsl3RCs9hBQLGk78BUgeivsvcBQ4BuS1gSvPufwOZuUbRUEQNeuXZk0aRJr167lXK/hOOcyT7K3ubaImS0AFsTNuy9m+iQwK8F23wa+ncrY4mVjBQGR5yJWrVrFSy+9xG233RZ2OM65NiT7zogNyMYKAqCwsJBLL72UrVu3smPHjrDDcc61IZ4gAtnQ1UZDpkyZQrdu3XjxxRd97GrnXL3sPCMmkK1NTAAdOnTgyiuvZN++faxfvz7scJxzbUR2nhETyNYmpqgxY8bQv39/Fi1axOnTp5vewDmX8TxBBLKls76GSOLaa6/lgw8+8C44nHOAJ4h62XwNIqq0tJRRo0bx+uuv+8NzzjlPEFHZXkFEXXvttUhi4cKFYYfinAuZJ4iAVxARXbt25bLLLmPLli1s27Yt7HCccyHyMyKR6gHwCiLw0Y9+lOLiYl544QVqaprbBZdzLlN4guCvCcIriIjc3FxmzJjB0aNHWbJkSdjhOOdC4mdEIs1L4BVErI985COMHTuW1157jYMHD4YdjnMuBJ4g8AqiIdOmTaOgoID58+fXf0fOuezhZ0S8gmhI586dmTZtGu+++y7Lly8POxznXJp5guCvCcIriLONGTOGoUOHsmjRIo4dOxZ2OM65NPIzIn4XU2Mkcd111yGJ5557zjvzcy6LeILAm5ia0r17d6ZPn87OnTtZtmxZ2OE459LEEwR+kToZ48aN44ILLmDRokV+V5NzWcLPiHgFkQxJzJw5k/z8fH73u9/Vf2fOuczlCQKvIJJVWFjIzJkz2b9/P4sWLQo7HOdcivkZEa8gmmPEiBGUl5ezdOlStm7dGnY4zrkU8gSBVxDNNW3aNPr168ezzz7L+++/H3Y4zrkU8TMiXkE0V15eHrfccgu1tbU8/fTTfj3CuQzlCQKvIFqiuLiYmTNnsmfPHl566aWww3HOpUBe2AG0BV5BtMzo0aPZs2cPy5Yto1+/flx44YVhh+Sca0X+kxmvIM7FtddeS2lpKc8//zx79+4NOxznXCvyMyJeQZyL3NxcZs2aRVFREU888QTHjx8POyTnXCvxBIF31neuOnfuzOzZs6mqquKxxx6juro67JCcc63Az4h4Z32toV+/ftx8883s3buXZ555xsePcC4DpDRBSJouaYuk7ZK+lmB5vqQnguVvSioN5hdLekXScUk/TmWM4E1MrWXEiBHMmDGDLVu28MILL3jPr861cylLEJJygZ8AM4Ay4DZJZXGr3QW8Z2ZDgR8A3wvmnwS+AXw1VfHF8ovUrWfy5MlMnTqV5cuX8/rrr4cdjnPuHKTyjDgZ2G5m75hZNfA4cGPcOjcC84Lpp4CrJMnMTpjZa0QSRcp5BdG6rr32WkaPHs2iRYtYsWJF2OE451oolc9BDAD2xLyvAKY0tI6Z1Uh6HygGDqcwrrN4BdG6JHHTTTdRXV3NH/7wBzp27MjYsWPDDss510ypPCMqwbz4Rulk1mn4ANI9klZIWnHo0KFmBRfLK4jWl5ubyy233EJpaSnPPvssmzdvDjsk51wzpTJBVAADY96XAPFPUtWvIykP6AYcTfYAZvagmZWbWXnv3r1bHGi0gpAS5SvXUh06dGD27NkMGDCA3/72t2zcuDHskJxzzZDKBLEcGCZpiKSOwGxgftw684E5wfQtwFNd8kcAAA/ySURBVMsWwq0vtbW15OTkeIJIgfz8fO644w4GDhzI008/zbp168IOyTmXpJQlCDOrAe4FFgKbgSfNbKOk+yXdEKz2EFAsaTvwFaD+VlhJO4HvA3MlVSS4A6rV1NXVefNSCuXn53P77bczePBgfve737Fq1aqwQ3LOJSGlnfWZ2QJgQdy8+2KmTwKzGti2NJWxxYpWEC51OnbsyKc+9SmefPJJnn/+eY4fP86ll17qVZtzbZifFYkkCK8gUi96TWLs2LG88sorLFiwwJ+4dq4N8+6+iTQxeQWRHrm5udx0000UFRXxxhtvUFlZyc0330zHjh3DDs05F8fPivg1iHSTxDXXXMOMGTPYunUrDz/8MMeOHQs7LOdcHE8QeBNTWCZPnsztt9/OsWPH+MUvfsGuXbvCDsk5F8MTBN7EFKaPfOQj3H333XTq1Il58+bxxhtveCd/zrURflbEK4iw9erVi7vvvpsRI0bw4osv8sQTT1BVVRV2WM5lPU8QeAXRFhQUFDBr1iymTZvGtm3b+PnPf87OnTvDDsu5rOZnRbyCaCskMXXqVD7zmc+Qm5vLvHnzePHFF6mpqQk7NOeykicIvIJoa0pKSvjsZz/LhAkTeOONN/jFL37Bu+++G3ZYzmUdPyviFURb1LFjR2bOnMltt91GVVUVDz30EAsXLvTxrp1LI08QeAXRlg0fPpzPf/7zTJgwgWXLlvHAAw+wefNmv9PJuTTwsyJeQbR1BQUFXH/99cydO5f8/HyefPJJfvOb33D4cFrHlXIu63iCwDvray8GDx7MZz/7WaZPn05FRQUPPPAAv//97zl+/HjYoTmXkbwvJryrjfYkJyeHKVOmMHr0aF599VVWrlzJunXrmDp1Kh/96Efp1KlT2CE6lzE8QeBNTO1RYWEhH//4x5kyZQovv/wyS5Ys4a233mLKlClMnTrVE4VzrcATBH6Ruj0rLi5m1qxZHDhwgFdffZXFixezbNkyJkyYwNSpU+nWrVvYITrXbnmCwK9BZIK+ffvyyU9+kgMHDvD666/z5ptv8tZbbzFq1CgmT57MgAEDfHAi55rJEwR+DSKT9O3bl5tvvpkrr7ySZcuWsXr1atavX0///v2ZNGkSo0aN8rEnnEuSJwi8gshE3bt3Z/r06VxxxRWsW7eO5cuXM3/+fF544QVGjRrF+PHjKSkp8arCuUZ4gsAriEyWn5/PpEmTKC8vZ/fu3axZs4YNGzawevVqevTowejRoxk9ejR9+vQJO1Tn2hxPEHgFkQ0kMXjwYAYPHsz06dPZvHkzGzZs4LXXXmPJkiX06tWLkSNHMnLkSPr16+eVhXN4gsDMMDOvILJIfn4+48aNY9y4cRw/fpxNmzaxefPm+mTRpUsXhg0bxrBhwxgyZAj5+flhh+xcKLI+QdTV1QF4BZGlioqKmDx5MpMnT+bDDz9ky5YtbNu2jQ0bNrBq1SpycnIoKSnh/PPPp7S0lAEDBpCXl/X/27gskfX/0mtrawG8gnB07tyZ8ePHM378eGpra9m9ezdvv/02O3bs4M9//jMAeXl5lJSUMHDgQAYOHEhJSYk/lOcylieIIEF4BeFi5ebmMmTIEIYMGQJAVVUVu3btqn+99tpr9T3K9uzZkwEDBnDeeefRv39/+vbtS0FBQZjhO9cqsj5BRJuYvIJwjenUqRMjRoxgxIgRAFRXV7N371727NnD3r172blzJ+vXr69fv3v37vTt25c+ffrQp08fevfuTc+ePenQoUNYH8G5Zsv6BOFNTK4lOnbsSGlpKaWlpfXzjh8/zr59+9i/fz/79+/n4MGDbN269YyxK3r06EFxcTE9e/asf/Xo0YPu3bv7tQ3X5mT9v0i/SO1aS1FRUf3dT1E1NTUcPnz4jNfRo0fZvXv3WaPjFRUV0a1bt/pX165d6dq1K126dKFLly4UFRV5EnFplfX/2ryCcKmUl5dHv3796Nev3xnzzYwTJ07w3nvv1b+OHTvGBx98wP79+9m6dSs1NTVn7a+goICioiKKioooLCykc+fO9X+jr06dOtW/OnTo4M90uBZLaYKQNB34EZAL/NLMvhu3PB/4NTAROALcamY7g2X/BNwF1AJfMrOFqYjRKwgXBkn1J/qBAweetdzMqKqq4oMPPqCysrL+dfz4cU6cOFHfnHXixAlOnTrV4HFycnIoKCg445Wfn1//6tixY/3f6HSHDh3o2LHjGX87dOhAXl4eubm5nnCySMoShKRc4CfANUAFsFzSfDPbFLPaXcB7ZjZU0mzge8CtksqA2cAo4DzgJUnDzay2teP0CsK1RZLqK4L46iNebW0tH374IR9++CFVVVV8+OGHnDx5kqqqKqqqqjh58mT969SpU1RWVnLy5Emqq6sbTS4NiSaLvLy8MxJHQ3/jXzk5OWdMR98nmo6+JDU6LemM6dh5zX1Fv//497F/s0UqK4jJwHYzewdA0uPAjUBsgrgR+GYw/RTwY0X+C9wIPG5mp4AdkrYH+1va2kF6BeHau9zc3PrrFM1lZpw+fZrq6uozXtF5p0+fTviqra2lpqYm4au6upqamhpqa2sTvurq6ur/v2uvGksejSWVZJbFH6eh5bHvhw8fznXXXdfCT9OwVCaIAcCemPcVwJSG1jGzGknvA8XB/GVx2w6IP4Cke4B7AAYNGtSiIPPz8ykrK6Nr164t2t659kxSffNSOpnZGckiftrMzpqOvmLfR6ejXebETid639gL/vqDMXZestON/U12Wfx31NDy+Pd9+/ZN8ptvnlQmiES1WPy30NA6yWyLmT0IPAhQXl5+9jechOiIZM659JHkd2S1A6lsV6kAYq++lQB7G1pHUh7QDTia5LbOOedSKJUJYjkwTNIQSR2JXHSeH7fOfGBOMH0L8LJFaqf5wGxJ+ZKGAMOAt1IYq3POuTgpq/GCawr3AguJ3Ob6sJltlHQ/sMLM5gMPAY8GF6GPEkkiBOs9SeSCdg3whVTcweScc65hSnRxpD0qLy+3FStWhB2Gc861K5JWmll5omV+b6dzzrmEPEE455xLyBOEc865hDxBOOecSyhjLlJLOgTsOodd9AIOt1I47UG2fV7wz5wt/DM3z2Az651oQcYkiHMlaUVDV/IzUbZ9XvDPnC38M7ceb2JyzjmXkCcI55xzCXmC+KsHww4gzbLt84J/5mzhn7mV+DUI55xzCXkF4ZxzLiFPEM455xLK+gQhabqkLZK2S/pa2PGkmqSBkl6RtFnSRkl/H3ZM6SIpV9JqSb8PO5Z0kNRd0lOS/hL89/5o2DGlmqR/CP5db5D0mKSCsGNqbZIelnRQ0oaYeT0lvShpW/C3R2scK6sThKRc4CfADKAMuE1SWbhRpVwN8I9mNhKYCnwhCz5z1N8Dm8MOIo1+BLxgZiOAC8nwzy5pAPAloNzMRhMZZmB2uFGlxCPA9Lh5XwMWmdkwYFHw/pxldYIAJgPbzewdM6sGHgduDDmmlDKzfWa2KpiuJHLSOGu870wjqQS4Dvhl2LGkg6SuwGVExlzBzKrN7Fi4UaVFHtApGKGyMxk4EqWZLSYyfk6sG4F5wfQ84KbWOFa2J4gBwJ6Y9xVkwckySlIpMB54M9xI0uKHwP8E6sIOJE3OBw4Bvwqa1X4pqTDsoFLJzN4F/h3YDewD3jezP4UbVdr0NbN9EPkRCPRpjZ1me4JQgnlZcd+vpCLgaeDLZvZB2PGkkqTrgYNmtjLsWNIoD5gA/NTMxgMnaKVmh7YqaHe/ERgCnAcUSroj3Kjat2xPEBXAwJj3JWRgSRpPUgciyeE3ZvZM2PGkwcXADZJ2EmlGvFLS/ws3pJSrACrMLFodPkUkYWSyq4EdZnbIzE4DzwAXhRxTuhyQ1B8g+HuwNXaa7QliOTBM0hBJHYlc0JofckwpJUlE2qU3m9n3w44nHczsn8ysxMxKifw3ftnMMvqXpZntB/ZIuiCYdRWRMd4z2W5gqqTOwb/zq8jwC/Mx5gNzguk5wHOtsdO81thJe2VmNZLuBRYSuePhYTPbGHJYqXYx8LfAeklrgnlfN7MFIcbkUuOLwG+CHz/vAJ8JOZ6UMrM3JT0FrCJyt95qMrDbDUmPAR8DekmqAP438F3gSUl3EUmUs1rlWN7VhnPOuUSyvYnJOedcAzxBOOecS8gThHPOuYQ8QTjnnEvIE4RzzrmEPEG4dk1SX0n/LekdSSslLZX0Ny3cV2lsD5mpJKlc0n+m41gxxxwn6ePpPKZr3zxBuHYreBjqWWCxmZ1vZhOJPAhXEm5kjZOUZ2YrzOxLqdh3I4vHAZ4gXNI8Qbj27Eqg2sx+Fp1hZrvM7L8AJBVI+pWk9UGHdVcE80slLZG0Kng12h2DpP6SFktaE4wzcGkwf3qw/VpJi4J5PSU9K2mdpGWSxgbzvynpQUl/An4t6WPRcSmCZQ9L+nNQCX0p5tjfCMZzeDEY3+CrCeJ7RNL3Jb0CfE/SZElvBJ/5DUkXBA/L3Q/cGnyOWyUVBsddHqyb0T0Zu+bL6iepXbs3ishTsw35AoCZjZE0AviTpOFE+qm5xsxOShoGPAaUN7KfTwELzez/BGOIdJbUG/gFcJmZ7ZDUM1j3W8BqM7tJ0pXAr4n8cgeYCFxiZlWSPhZ3jBHAFUAXYIuknxIZw+ETRHrczQs+a0MdDg4Hrjaz2mhX30FPAVcD3zGzT0i6j8hYCfcCSPoOkW5H7pTUHXhL0ktmdqKR78JlEU8QLmNI+glwCZGqYlIw/V8AZvYXSbuInEh3AT+WNA6oDeY1ZjnwcNDJ4bNmtiY4wS82sx3B/qP9819C5KSOmb0sqVhSt2DZfDOrauAYfzCzU8ApSQeBvsG+notuI+n5RmL8rZnVBtPdgHlB8jOgQwPbXEukE8NoVVIADCJ7+i9yTfAE4dqzjQQnYwAz+4KkXsCKYFai7twB/gE4QOQXeg5wsrGDmNliSZcRGXDoUUn/FzhG4q7hG+tCvrFf5qdipmuJ/L/ZUPyJxO77X4BXzOxvFBnz488NbCPgE2a2pRnHcVnEr0G49uxloEDS38XM6xwzvRi4HSBoWhoEbCHyC3ufmdUR6bgwt7GDSBpMZDyJXxDpCXcCsBS4XNKQYJ1oE1PsMT8GHD6H8TZeA2YG11KKiCSoZHQD3g2m58bMryTShBW1EPhicLEfSeNbGKfLUJ4gXLtlkZ4mbyJyot4h6S0iwy3+r2CVB4BcSeuBJ4C5QTPOA8AcScuINC811eb+MWCNpNVEKpYfmdkh4B7gGUlrg/0DfBMol7SOSA+bc87eXdKfbzmRbpzXEhnbYAXwfhKb/hvwr5Je58zk9wpQFr1ITaTS6ACsC27v/ZeWxuoyk/fm6lwbJqnIzI5L6kykOrknOqa4c6nm1yCca9selFRG5ALyPE8OLp28gnDOOZeQX4NwzjmXkCcI55xzCXmCcM45l5AnCOeccwl5gnDOOZfQ/w/XlJiI9fFLowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prior.plot(label='prior', color='gray')\n",
    "\n",
    "plt.xlabel('Goal scoring rate')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Prior distribution')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution represents our prior knowledge about goal scoring: λ is usually less than 2, occasionally as high as 6, and seldom higher than that.  And the mean is about 1.4.\n",
    "\n",
    "As usual, reasonable people could disagree about the details of the prior.  Later, we can try a different prior and see how much effect it has on the results.\n",
    "\n",
    "But first, it's time to do an update."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The update\n",
    "\n",
    "Suppose you are given the goal-scoring rate, λ, and asked to compute the probability of scoring a number of goals, $k$.  That is precisely the question we answered by computing the Poisson distribution:\n",
    "\n",
    "$ f(k; λ) = \\lambda^k \\exp(-\\lambda) ~/~ k! $\n",
    "\n",
    "For example, if λ is 1.4, the probability of scoring 4 goals in a game is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.039471954028253146"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = 4\n",
    "poisson.pmf(4, λ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose we are have an array of possible values for λ, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "λs = prior.index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the likelihood of the data for each hypothetical value of λ, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 4\n",
    "likelihood = poisson.pmf(k, λs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's all we need to do the update.  \n",
    "\n",
    "Here's the function we used in the previous notebook; it takes a prior distribution and an array of likelihoods, does a Bayesian update, and returns the total probability of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bayes_update(pmf, likelihood):\n",
    "    \"\"\"Do a Bayesian update.\n",
    "    \n",
    "    pmf: Series that represents the prior\n",
    "    likelihood: sequence of likelihoods\n",
    "    \n",
    "    returns: float probability of the data\n",
    "    \"\"\"\n",
    "    pmf *= likelihood\n",
    "    prob_data = pmf.sum()\n",
    "    pmf /= prob_data\n",
    "    return prob_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example, France scored 4 goals, so I'll make a copy of the prior and update it with the likelihoods we just computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05015532557804499"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "france = prior.copy()\n",
    "bayes_update(france, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under the prior distribution, the probability of scoring 4 goals was about 5%.\n",
    "\n",
    "Here's what the posterior distribution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxU9bn48c8zk31fWAJJSIAAIRsBAgIqiguLiluhpbWtWm1trV1v22v7u1209Xa/rd7WLlYUva51hWpFxQVBEILsCZAAgYQ1C4EEss58f3/MSRxCQhKYNTzv1+u8mDnrMwOcZ57v95zvEWMMSimlVFc2fweglFIqMGmCUEop1S1NEEoppbqlCUIppVS3NEEopZTqliYIpZRS3dIEofxORBpFZJSX9n25iFS5vd8uIpd7aN+3iMibbu+NiGR5Yt/W/rzyvYjIOBHZKCINIvJNT+9fDRyaIC4gIlIhIk3WieeIiDwmIjHnsb9M66QYcj5xGWNijDF7zmcf/ThWrjHmvbOt09fPZYx5yhgz2xNxich7InJnl/1763v5AfCeMSbWGPPQ+e5MRH4mIm3Wv6uO6QceiFP5mSaIC898Y0wMMAmYAvyXvwI538RyvtsH67E9IAPYfi4bnuVzP2cltI7pN91sKyKi55wgon9ZFyhjzAHg30AegIgMF5GlIlInIuUi8uWOdUVkqogUi8gJq/L4H2vRSuvPeutX43Rr/S+JSKmIHBOR5SKS4bYvIyJfF5EyoMxtXpb1Ol5EnhCRahHZJyL/1XFSEZHbRGS1iPxBROqAn3X9XCISKSKPW8cuwZUE3ZdXiMhV/f1c3R3bmreqSwjXiMgeEakRkd+6xf4zEfk/tzg6qxQReQC4FPiTdbw/ncP3skpEfmd97r0iMq+7v3cReQeY5Xassef7nZ+NVRk9ICKrgVPAKBG53fr30WB9V3e5rX+5iFSJyH+IyFEROSQit7stjxSR31txHrc+d6S1bJqIfCgi9SKyWTzUlHhBM8bodIFMQAVwlfU6HdevyJ9b798HHgYigEKgGrjSWrYG+IL1OgaYZr3OBAwQ4naMG4FyYDwQgqtC+dBtuQHeApKASLd5WdbrJ4BXgVhr/7uAO6xltwHtwDesfUd28xl/BXxg7T8d2AZU9fAd9OdznXFsa96qLp/tXevYI6zY77SW/Qz4P7d1TzsG8F7Hul3219fvpQ34MmAHvgYcBKSHfwenHcsD3/lpn62bY+0Hcq3tQ4FrgdGAAJfhShyTrPUvt453v7XuNdbyRGv5n619plqfdQYQbr2vtda3AVdb7wf7+/9dME9+D0AnH/5lu06OjUA9sA9XQojEdSJ1ALFu6/4SeNx6vRK4DxjUZX+nneSsef/uOLlY723Wf/AM670BruiyHwNkWf/hW4Act2V34Wov7zhZ7e/lM+4B5rq9/wo9J4j+fK4zjk33CcL92HcDK6zXp51Eux6DsySIPn4v5W7LoqxtU3r4jjqP5aHv/GdAq/XvqmMa7nas+3vZ/hXgW9bry4GmLt/9UWCa9W+pCZjQzT7+E3iyy7zlwK3+/n8XzJM2MV14bjTGJBhjMowxdxtjmoDhQJ0xpsFtvX24fpUB3AGMBXaIyHoRue4s+88AHrTK/HqgDtcvxVS3dSp72HYQEGYdu7s4zrZth+Fd1tnX04r073P15dhd19lnxXO++vK9HO54YYw5Zb3sywUInvjOAZ63/l11TAd72l5E5onIWqs5sx7Xr/5BbqvUGmPa3d6fsj7LIFwV7u5ujp8BLOz4d2ft9xJgWB9iVz3QBKHA1RyRJCKxbvNGAAcAjDFlxpjPAkOAXwMviEg0rl+pXVUCd3U5WUQaYz50W6enIYRrcDWVZLjN64yjl207HMJVEblv361+fq6+HJtujt1xojyJ65d9h5R+7Lsv38u58sR33pvO7UUkHHgR+B0w1BiTALyO60dEX2JtxtU81VUlrgrC/d9dtDHmV+cZ+wVNE4TCGFMJfAj8UkQiRKQA16/rpwBE5PMiMtgY48TVfACuJqlqwAm4X6v/V+CHIpJrbRsvIgv7GIcDeB54QERixdW5/V3g/86+5Wmet46fKCJpuNrOu9XPz9VX37eOnQ58C3jOmr8JmCkiI0QkHvhhl+2O9HQ8D30v3fLmvnsQhqvPoBpotzrT+3SpsPX3tBj4H3FdVGG3LiAIt+KdLyJzrPkRVod3mpc+xwVBE4Tq8Flc7eIHgZeBnxpj3rKWzQW2i0gj8CCwyBjTbDVlPACstsr6acaYl3H9Gn9WRE7g6iTu9oqaHnwD16/tPcAq4GlcJ4W+ug9XE8le4E3gybOs2+fP1Y/jvwpswJUQXgMeBbC+y+eALdbyf3XZ7kFggXUVUnf3Jpzv93I23tz3aaxmzG/iSkrHgM8BS/uxi+8BW4H1uJovfw3YrB85NwA/wpV8KoHvo+e48yJWZ45SSil1Gs2uSimluqUJQimlVLc0QSillOqWJgillFLdCuYBx04zaNAgk5mZ6e8wlFIqqGzYsKHGGDO4u2UDJkFkZmZSXFzs7zCUUiqoiEiPow1oE5NSSqluaYJQSinVLU0QSimlujVg+iCUUqpDW1sbVVVVNDc3+zuUgBEREUFaWhqhoaF93kYThFJqwKmqqiI2NpbMzExE+jJQ7MBmjKG2tpaqqipGjhzZ5+20iUkpNeA0NzeTnJysycEiIiQnJ/e7otIEoZQakDQ5nO5cvg9NEP1kjGHjxo3U1NT4OxSllPIq7YPoB6fTybJly9i0aRN2u52rrrqKiy66SH+pKKXO2U9+8hNmzpzJVVdd5e9QzqAJoo+cTievvPIKW7duZcaMGdTU1LB8+XJ27drFzTffTExMXx7/q5RSn3A4HNx///393sZut3spotNpE1MfGGN46aWX2Lp1K1dccQVXX301ixYt4rrrrqOyspJXX30VffCSUspdRUUF2dnZ3HrrrRQUFLBgwQJOnTpFZmYm999/P5dccgn//Oc/ue2223jhhRcAWLFiBRMnTiQ/P58vfelLtLS0AJyxja9oBdEHR48eZfv27cycOZNLL70UcHX4TJ48mdbWVt58803KysoYO3asnyNVSnX1xhtvcPjwYY/uMyUlhblz5/a63s6dO3n00Ue5+OKL+dKXvsTDDz8MuO5JWLVqVWd84Lry6rbbbmPFihWMHTuWL37xi/zlL3/h29/+9hnb+IpWEH2wd+9eACZNmnTGsqlTp5KcnMzy5ctpb2/3dWhKqQCWnp7OxRdfDMDnP//5zhP8Zz7zmTPW3blzJyNHjuz8oXnrrbeycuXKzuXdbeNtWkH0wd69e0lKSiI+Pv6MZXa7nblz5/LUU0/x0Ucfdf5jUEoFhr780veWrhewdLyPjo4+Y93emqm728bbtILohdPppKKi4qx3H2ZlZTF27FhWrlxJQ0ODD6NTSgWy/fv3s2bNGgCeeeYZLrnkkh7Xzc7OpqKigvLycgCefPJJLrvsMp/E2RNNEL04ePAgra2tvd6ePmfOHBwOx2kloVLqwjZ+/HiWLFlCQUEBdXV1fO1rX+tx3YiICB577DEWLlxIfn4+NpuNr371qz6M9kxebWISkbnAg4Ad+Icx5lddlocDTwCTgVrgM8aYChHJBEqBndaqa40xfvmmOvofentaXVJSEvn5+WzevJkrr7ySiIgIH0SnlApkNpuNv/71r6fNq6ioOO39448/3vn6yiuvZOPGjWfsp+s2vuK1CkJE7MCfgXlADvBZEcnpstodwDFjTBbwB+DXbst2G2MKrclvaXTv3r0MHTq0T+1/U6dOpa2tjU2bNvkgMqWU8i5vNjFNBcqNMXuMMa3As8ANXda5AVhivX4BuFIC6Lbk9vZ2Kisr+zz64bBhw0hPT2f9+vV6X4RSF7jMzEy2bdvm7zDOizcTRCpQ6fa+yprX7TrGmHbgOJBsLRspIhtF5H0RubS7A4jIV0SkWESKq6urPRs9UFlZSXt7e7+Gx50yZQp1dXWdHU1KKRWsvJkguqsEuv6s7mmdQ8AIY8xE4LvA0yISd8aKxvzdGFNkjCkaPHjweQfc1d69exERMjIy+rxNTk4OMTExrFu3zuPxKKWUL3kzQVQB6W7v04CDPa0jIiFAPFBnjGkxxtQCGGM2ALsBn9+mXFFRQWpqKuHh4X3exm63M3nyZMrLy6mtrfVidEop5V3eTBDrgTEiMlJEwoBFwNIu6ywFbrVeLwDeMcYYERlsdXIjIqOAMcAeL8Z6htbWVg4cONDr1UvdKSoqwmazUVxc7PnAlFLKR7yWIKw+hXuA5bguWX3eGLNdRO4Xkeut1R4FkkWkHFdT0r3W/JnAFhHZjKvz+qvGmDpvxdqduro6nE4nw4YN6/e2MTExjBs3jq1bt+J0Or0QnVIq0NntdgoLCzsnf12qej68eh+EMeZ14PUu837i9roZWNjNdi8CL3oztt503BEdGxt7TtsXFBRQWlrK7t27GTNmjCdDU0oFgcjIyLNe8t7e3k5ISGCPdqR3UvfgfBPEmDFjiIyMZPPmzZ4MSykVxB5//HEWLlzI/PnzmT17No2NjVx55ZVMmjSJ/Px8Xn31VcDV/zl+/Hi+/OUvk5uby+zZs2lqagKgvLycq666igkTJjBp0iR2794NwG9/+1umTJlCQUEBP/3pTz0Sb2CnLz/qSBDn+iAgu91OXl4eGzdupLm5We+sVspP7lu2nZKDJzy6z5zhcfx0fu5Z12lqaqKwsBCAkSNH8vLLLwOwZs0atmzZQlJSEu3t7bz88svExcVRU1PDtGnTuP56Vwt8WVkZzzzzDI888gif/vSnefHFF/n85z/PLbfcwr333stNN91Ec3MzTqez85ED69atwxjD9ddfz8qVK5k5c+Z5fU5NED1oaGggKirqvErAgoIC1q9fT0lJSbdDhSulBq6empiuvvpqkpKSANcIrj/60Y9YuXIlNpuNAwcOcOTIEcCVVDoSzOTJk6moqKChoYEDBw5w0003AXT+8HzzzTd58803mThxIgCNjY2UlZVpgvCWhoaGc25e6pCamkpSUhJbtmzRBKGUn/T2S9/X3Ifteeqpp6iurmbDhg2EhoaSmZlJc3MzwGmX19vtdpqamnococEYww9/+EPuuusuj8aqfRA98ESCEBEmTJjAvn37qK+v91BkSqmB4vjx4wwZMoTQ0FDeffdd9u3bd9b14+LiSEtL45VXXgGgpaWFU6dOMWfOHBYvXkxjYyMABw4c4OjRo+cdnyaIHngiQYCrmQlgy5Yt570vpdTAcsstt1BcXExRURFPPfUU2dnZvW7z5JNP8tBDD1FQUMCMGTM4fPgws2fP5nOf+xzTp08nPz+fBQsWeOTZNDJQBpUrKioynroxzel08vOf/5yZM2cya9as897fY489RlNTE3fffbcHolNK9aa0tJTx48f7O4yA0933IiIbjDFF3a2vFUQ3Oso0T1QQALm5uVRXV3uk5FNKKV/RBNGN870HoquOjL19+3aP7E8ppXxBE0Q3PJ0gYmNjyczMpKSkRJ8ToZSP6P+1053L96EJohueThDgGga8pqZGm5mU8oGIiAhqa2s1SViMMdTW1vb7hl29D6IbDQ0NiEifHjPaVzk5Ofz73/9m+/btDB061GP7VUqdKS0tjaqqKrzxILFgFRERQVpaWr+20QTRjYaGBmJiYrDZPFdgRUdHk5mZyfbt25k1axYB9GRVpQac0NDQfj0JUnVPm5i64al7ILrKzc2lrq6Ow4cPe3zfSinlaZoguuGtBDF+/HhERK9mUkoFBU0Q3fBWgoiKimLUqFF6NZNSKihoguiivb2dpqYmryQIcFURx44d6xyxUSmlApUmiC68cYmru+zsbESEkpISr+xfKaU8RRNEF95OENHR0WRkZFBaWuqV/SullKdogujC2wkCPrlpTq/RVkoFMk0QXfgiQXQM6avNTEqpQKYJoouGhgbsdjuRkZFeO0ZsbCwjRozQBKGUCmiaILrouMTV23c6jx8/nqNHj1JbW+vV4yil1LnSBNGFt+6B6KpjCHCtIpRSgUoTRBe+ShDx8fGkpqbq1UxKqYClCaILXyUIcFURhw4dor6+3ifHU0qp/tAE4aalpYXW1lafJYicnBwArSKUUgFJE4SbjmdRx8TE+OR4iYmJpKSkaIJQSgUkryYIEZkrIjtFpFxE7u1mebiIPGct/0hEMrssHyEijSLyPW/G2aGtrQ2AsLAwXxwOcDUzVVZWdt5/oZRSgcJrCUJE7MCfgXlADvBZEcnpstodwDFjTBbwB+DXXZb/Afi3t2LsyuFwAHj0QUG96biaSasIpVSg8eaZcCpQbozZY4xpBZ4Fbuiyzg3AEuv1C8CVYt2AICI3AnsAnz08wel0AmC32311SAYPHsygQYM0QSilAo43E0QqUOn2vsqa1+06xph24DiQLCLRwH8C93kxvjP4o4IAVxWxb98+Tp486dPjKqXU2XjzTNjdrchdn5LT0zr3AX8wxjSe9QAiXxGRYhEp9sTAd/6oIMCVIIwx7Ny506fHVUqps/FmgqgC0t3epwEHe1pHREKAeKAOuAj4jYhUAN8GfiQi93Q9gDHm78aYImNM0eDBg8874I4KwtcJIiUlhYSEBG1mUkoFFG8miPXAGBEZKSJhwCJgaZd1lgK3Wq8XAO8Yl0uNMZnGmEzgj8B/G2P+5MVYgU8qCF83MYkIOTk57Nmzh6amJp8eWymleuK1M6HVp3APsBwoBZ43xmwXkftF5HprtUdx9TmUA98FzrgU1pf8VUGA66Y5p9PJrl27fH5spZTqTog3d26MeR14vcu8n7i9bgYW9rKPn3kluG74q4IAGD58OHFxcZSUlDBhwgSfH18ppbrSO6nd+LOC6Ghm2r17Ny0tLT4/vlJKdaUJwo0/KwhwNTM5HA5tZlJKBQRNEG78WUEApKWlERsbq8+IUEoFBE0QbvxdQYgI48ePp7y8nNbWVp8fv83hZF/tSVrbnT4/tlIq8Hi1kzrY+LuCAFcz07p16ygrKyM3N9frxzt8vJnHVu9lw75jbDt4nOY2J2F2G+OHxTIhPYFbZ2QyerBvRrdVSgUWTRBu/DXUhrv09HSio6MpKSnxaoIwxvDChiru/1cJLW1O8tPi+dzUDLKGxFBRe5ItVfX8s7iKpz/az20zMvnGlWOIjwz1WjxKqcCjCcKNv4bacGez2Rg/fjybN2+mtbXVK0OPHzvZynef38S7O6uZOjKJ33yqgMxB0WesV93Qwu/f3Mmjq/fy8sYDPLhoIpeMGeTxeJRSgUn7INwEQgUBkJubS1tbG2VlZR7fd3ObgzuWrGf17lp+Nj+HZ788rdvkADA4NpxffaqAZfdcwqCYcG5/fB3LNncdLUUpNVBpgnDjdDqx2WxYI477zYgRIzqbmTzJ4TR869mNbKys56FFhdx28Uhstt4/a15qPM9/dToT0xP55rMbeXz1Xo/GpZQKTJog3DgcDr9XD+CqYHJycti1a5dHr2b6xWslLN9+hB9fm8PcvGH92jY+MpQn7pjK1eOH8rNlJSxepUlCqYHO/2fDAOJ0Ov3a/+AuNzeX9vZ2j9009+KGKh5bXcEdl4zkS5eMPKd9RITaefiWSczJHcovXivh3R1HPRKbUiowaYJwEygVBLiuZoqJifFIM9Oxk6384rUSijIS+X/XjD+vfYXYbfzhM4WMHxbHN57ZyM7D+ixtpQaqwDgbBohAqiA6mpnKysrOe2ymX7+xg4bmdh64Kb9PfQ69iQoL4R+3FhEVZueOJeupbdSxo5QaiDRBuOnopA4UnmhmKq6o49n1ldxxyUjGpcR6LLZh8ZE88sUijja08P0XtmBM14cFKqWCXeCcDQOAw+EImAoCXM1MsbGxbN++/Zy2b3M4+a9XtjE8PoJvXjnGw9HBhPQEfjgvm3d2HOXpdfs9vn+llH9pgnATaBVExxDg5eXlNDc393v7Z9btZ8fhBn56fS7R4d65J/LW6ZlcOmYQv/hXKXuqz/oIcaVUkAmcs2EACLQKAiA/Px+Hw9Hv51W3OZz87f09FGUkMjtnqJeiA5tN+O2CCYSF2PjO85tpc+hAf0oNFJog3ATSVUwdhg8fTmJiItu2bevXdq9uOsiB+ibunjXa6zf+pcRH8MBNeWyurOfvK/d49VhKKd8JrLOhnwXSVUwdRIS8vDz27t1LY2PfmnCcTsNf399Ndkoss8YN8XKELtcVDOea/BQeWlHG/tpTPjmmUsq7NEG4CcQKAiAvLw9jTJ/viXiz5AjlRxv52uXerx7c/eS6XEJswk+WbtOrmpQaAALvbOhHgVhBAAwZMoQhQ4b0qZnJGMNf3itnRFIU1+b3bziN85USH8F3rh7LezureWPbYZ8eWynleX1KECLyoohcKyIDOqEEYid1h7y8PCorK6mvrz/reh/urmVz1XHuumwUIXbf/3XdNiOT8cPiuG9ZCY0t7T4/vlLKc/p6BvkL8DmgTER+JSLZXozJbwLtMld3eXl5AL1WEU+u2UdSdBifmpTmi7DOEGK38cBNeRxpaOZ/V3h+uHKllO/06WxojHnbGHMLMAmoAN4SkQ9F5HYRGTCPGQvkCiIxMZG0tLSzJoiaxhbeLj3CzRNTiQj13+eYNCKRBZPSeGx1hXZYKxXE+vxzWUSSgduAO4GNwIO4EsZbXonMDwK5ggDXPRFHjhzh6NHuR1F9+eMDtDsNn5mS7uPIzvS9OeOw24RfvdG/+zeUUoGjr30QLwEfAFHAfGPM9caY54wx3wAGzBPtA7mCANfYTDabjc2bN5+xzBjDc8WVTByRwJihnhtz6VwNjYvgq5eN5vWth1lfUefvcJRS56CvP5f/YYzJMcb80hhzCEBEwgGMMUVei87HAr2CiI6OJisri61bt3Y+P7vDx/vrKT/ayGeK/F89dPjKzFGkxEXw83+V4HTqZa9KBZu+ng1/0c28NZ4MJBAE6n0Q7goKCmhoaKCiouK0+c+vryQqzM51E4b7J7BuRIbZ+cHccWypOs6rmw/4OxylVD+d9WwoIikiMhmIFJGJIjLJmi7H1dw0oAR6ExPA2LFjCQ8PZ8uWLZ3zTra0868tB7k2fxgxXhqU71zdWJhKXmocv39zFy3tDn+Ho5Tqh95+Ls8BfgekAf8D/N6avgv8qLedi8hcEdkpIuUicm83y8NF5Dlr+UcikmnNnyoim6xps4jc1L+PdW4CvYkJIDQ0lJycHEpLS2lrawPgta2HONnqCIjO6a5sNuEHc7KpOtbEMx/pkOBKBZOzng2NMUuMMbOA24wxs9ym640xL51tWxGxA38G5gE5wGdFJKfLancAx4wxWcAfgF9b87cBRcaYQmAu8DcR8fpP42CoIMDVzNTa2sqOHTsAWLb5IJnJUUzOSPRzZN27dMwgpo9K5n/fKdeb55QKIr01MX3eepkpIt/tOvWy76lAuTFmjzGmFXgWuKHLOjcAS6zXLwBXiogYY04ZYzrOJBGAT3o4g6GCAMjIyCA+Pp4tW7ZQd7KVD3fXck3+MJ+Ou9QfIsL3546j9mQri1ft9Xc4Sqk+6u1sGG39GQPEdjOdTSpQ6fa+yprX7TpWQjgOJAOIyEUish3YCnzVLWF0EpGviEixiBRXV1f3Es7ZGWMCdiymrkSEgoICdu/ezasbKnA4Ddf4eNyl/po0wvVcikdW7qHuZKu/w1FK9UFvTUx/s/68r7upl31393O2ayXQ4zrGmI+MMbnAFOCHIhLRTXx/N8YUGWOKBg8e3Es4Z9dx2WgwJAiACRMmYIzhpfV7GZEURe7wOH+H1KvvzRlHY2s7f3mv3N+hKKX64Kzt+iLy0NmWG2O+eZbFVYB7r2kacLCHdaqsPoZ44LS7qowxpSJyEsgDis8Wz/noSBDB0MQEkJyczKDUDLaVt3HXZSMCtnnJ3dihsdxUmMoTa/bx5UtHMSTujJyvlAogvZ0NN/Qync16YIyIjBSRMGARsLTLOkuBW63XC4B3jDHG2iYEQEQygHG4xoDyGofDdQlmsFQQAKcSx2AQJg8OjqQG8M0rx9DuNDz83m5/h6KU6sVZKwhjzJKzLe9l23YRuQdYDtiBxcaY7SJyP1BsjFkKPAo8KSLluCqHRdbmlwD3ikgb4ATuNsbUnGssfRFsFQTAlmN2YqWFk1U7oGicv8Ppk8xB0XxqUipPr9vPXZeNYlh8pL9DUkr1oLcmpj8aY74tIsvo5koiY8z1Z9veGPM68HqXeT9xe90MLOxmuyeBJ88eumcFWwVRf6qVNXvquGxYKCUlW7nmmnmEhYX5O6w++cYVY3jp4wP8+d1yfnFjvr/DUUr1oLd7CzpO0r/zdiD+FmwVxJvbj7hGbp0xhjWvfUxJSQmFhYX+DqtP0pOi+PSUdJ5bX8lXLxtNWuKAuylfqQGht6uYNlh/vo9r7KVjuJqC1ljzBoxgqyDeLDlCakIkV00aS1JSEhs3bvR3SP1yz6wsBOFP7+gVTUoFqr4O930tsBt4CPgTUC4i87wZmK91JIhgqCCa2xysLq/hyvFDsNlsTJw4kf3791NT49VuGo8anhDJoqnpvLChiso6faiQUoGor2fD3wOzjDGXG2MuA2bhGhpjwAim+yDW7qmlqc3BrOwhABQWFmKz2diwobcLywLL1y4fjU2Eh/W+CKUCUl8TxFFjjPv/4j1A9481C1LBVEG8s+MokaF2po9KBiAmJobs7Gw2b95Me3vwjHU0LN5VRfyzuIqqY1pFKBVoehuL6WYRuRnYLiKvi8htInIrsAzXfQ4DRrBUEMYY3tlxlIuzBp323OnJkyfT1NRESUmJH6Prv44q4s/v6n0RSgWa3n4uz7emCOAIcBlwOVANBObQoecoWDqpy442UnWsiSus5qUOI0eOJDExMeiamT6pIiq1ilAqwPR2o9ztvgrE34LlMtcVpa6Wva4JQkSYPHkyb7/9NtXV1Zzv2FS+9LXLR/Psukr+/O5ufnmz3hehVKDo61VMESLydRF5WEQWd0zeDs6XgqWCeHfHUXKGxZESf+Y4RoWFhdjt9qCsIj4zJZ0XNlRyoL7J3+EopSx9/bn8JJCC6wlz7+MaeK/BW0H5QzBUEPWnWineV8eV44d0uzw6OhfqdCYAACAASURBVJrx48ezefPmzqfNBYuvXT4agIff1SualAoUfT0bZhljfgyctMZnuhYYUG0BwVBBvL+rGqc5s3nJXVFREc3NzWzdutWHkZ2/4QmRfLooneeLKzmoVYRSAaGvCaLj52i9iOThGpY70ysR+UkwVBDv7awmKTqMCWkJPa4zYsQIhg4dyrp16zDGJw/i85i7Z2UB6H0RSgWIvp4N/y4iicCPcQ3RXcInz48eEAK9gjDG8EFZDZeOGYTN1vOzH0SEqVOncuTIEfbv3+/DCM9fakIkCyan8/z6Kg4d1ypCKX/rU4IwxvzDGHPMGPO+MWaUMWZIx9PmBopAryB2HG6gprGFS7IG9bpufn4+ERERfPTRRz6IzLPuvnw0TmP4iz4vQim/6+tVTMki8r8i8rGIbBCRP4pIsreD86VAryA+KHM9c/vSMb1fvhoaGsqkSZPYsWMHx48f93ZoHpWeFMWCyWk8u66Sw8eb/R2OUhe0vv5cfhbX0BqfwvXktxrgOW8F5Q+BPtTGB2U1jB0a0+3lrd2ZMmUKAMXFXntKq9d8fVaWVUVoX4RS/tTXs2GSMebnxpi91vQLoOee0iAUyENtNLc5WLe3jkuy+n7zW0JCAuPGjWPDhg1Bd8lrRxXxjFYRSvlVXxPEuyKySERs1vRp4DVvBuZrgdzEtL6ijpZ2J5eO7b3/wd1FF11EU1MTW7Zs8VJk3tNRRfz1fe2LUMpfehusr0FETgB3AU8Drdb0LPAd74fnO4HcSb2qrIYwu42LRib1a7uMjAyGDRvGmjVrgu6S1/SkKD41KY2n1+3nyAmtIpTyh96eKBdrjImz/rQZY0KsyWaMifNVkL4QyBXEyrIaJmckEhXW2xNiTycizJgxg9raWnbu3Oml6Lzn67OycDr1iial/KXPP5dF5HoR+Z01XefNoPzB6XQiIoj0fI+BP1Q3tFB66ES/m5c65OTkEB8fz5o1azwcmfeNSP6kitC+CKV8r6+Xuf4K+BauG+RKgG9Z8wYMh8MRkM1Lq8tdjxG9tB8d1O5sNhvTpk1j//79VFVVeTI0n7jnClcVoXdXK+V7fT0jXgNcbYxZbIxZDMy15g0YTqczQJuXqkmMCiV3+Lm36E2cOJGIiIigrCLSk6JYWJTOs+t0pFelfK0/P5ndL2uN93Qg/haIFYQxhjW7a5kx+uzDa/QmPDycyZMnU1paSm1trQcj9I17rsjCYPjTO1pFKOVLfT0j/hLYKCKPi8gSYAPw394Ly/cCsYKoqD3FoePNzMg6/5vWp02bht1uZ/Xq1R6IzLdSEyJZNGUE/yyupLJOnzqnlK/0miDE1Wu7CpgGvGRN040xz3o5Np8KxAqio/9hxuhz66B2FxMTw6RJk9i8eXPQDb8BcPes0dhswkMryvwdilIXjF7PiMZ1Af0rxphDxpilxphXjTGHfRCbTwViBbFmdy3D4iPITI7yyP5mzJgBEJRVxLD4SD43dQQvbTzAnupGf4ej1AWhrz+Z14rIFK9G4meBVkE4nYY1e1z9D5669DY+Pp7CwkI+/vhjGhqC74GAX5+VRZjdxh/e1ipCKV/o6xlxFq4ksVtEtojIVhHpdfwGEZkrIjtFpFxE7u1mebiIPGct/0hEMq35V1ujxm61/ryiPx/qXARaBbHjcAN1J1uZMdqzg+ZecsklOJ3OoLyiaXBsOLdfnMmyzQcpPXTC3+EoNeD1NUHMA0YBVwDzgeusP3skInbgz9a2OcBnRSSny2p3AMeMMVnAH/jkIUQ1wHxjTD5wK65nYnuVw+EIqATx4W5X/8N0DyeIxMRE8vPzKS4u5uTJkx7dty/cNXM0sREh/P7NXf4ORakBr7exmCJE5NvA93Hd+3DAGLOvY+pl31OBcmPMHmNMx/hNN3RZ5wZgifX6BeBKERFjzEZjzEFr/nYgQkTC+/G5+s3pdAZUE9Oa3bWMHBTN8IRIj+/70ksvpb29PSj7IuKjQrlr5ijeLj3Cx/uP+TscpQa03s6IS4AiYCuuSuD3/dh3KlDp9r7KmtftOsaYduA40PUn86eAjcaYlq4HEJGviEixiBRXV1f3I7QzBVIF0e5w8tHeOo83L3UYNGgQBQUFrF+/Pij7Im6/eCTJ0WH8bnnwjS+lVDDpLUHkGGM+bz1edAFwaT/23V3PatchRc+6jojk4mp2uqu7Axhj/m6MKTLGFA0efG5DUXQIpApiy4HjNLa0e+Ty1p5cdtllOJ1OVq5c6bVjeEt0eAhfn5XFh7trO5+0p5TyvN7OiJ1PmrF+4fdHFZDu9j4NONjTOiISgusO7TrrfRrwMvBFY4zXh/MMpApizW7X3c7TRvVveO/+SExMZOLEiXz88cccOxZ8TTW3TBtBWmIkv/r3DpzO4BrKXKlg0VuCmCAiJ6ypASjoeG09J+Js1gNjRGSkiIQBi4ClXdZZiqsTGlwVyjvGGCMiCbgeSPRDY4xPGsoDqYL4cHcN2SmxJMd4tduFmTNnIiJBWUWEh9j53uxxbD94gmVbuv7uUEp5Qm/Pg7Bbz4PoeCZEiNvrs44eZ1Uc9wDLgVLgeWPMdhG5X0Sut1Z7FEgWkXLgu0DHpbD3AFnAj0VkkzUNOY/P2atAqSBa2h0UVxzzavNSh7i4OKZMmcLmzZs53z4cf7h+wnDGD4vjd2/upLXd6e9wlBpwvPqT2RjzujFmrDFmtDHmAWveT4wxS63XzcaYhcaYLGPMVGPMHmv+L4wx0caYQrfpqDdjDZQKYuP+elranR6/vLUnl1xyCaGhobz99ts+OZ4n2WzCf84dR2VdE09/1NtFdUqp/vL/GTFABEoFsWZ3LTaBqf18vOi5io6O5tJLL2XXrl3s3bvXJ8f0pMvGDmb6qGQeeqecE81tvW+glOozTRCWQBlqY82eWnKHxxMfGeqzY1500UXEx8fz1ltvBd2zq0WEH10znrqTrTz8rj6aVClP8v8ZMUAEQhNTU6uDTfvrfda81CE0NJQrrriCQ4cOsXXrVp8e2xPy0+K5eVIqi1ft1eHAlfIgTRCWQGhi2rDvGK0O3/U/uMvPz2fYsGGsWLGCtrbga6r5/pxx2GzwG715TimP0QRhCYTB+tbsqcFuE6Zk+qb/wZ2IMHv2bE6cOBGUQ3AMi4/kK5eOYtnmg2zYF3z3dSgViDRBWAKhD2LN7loK0uKJCQ/xy/EzMzPJzc1l9erVQXnz3F2XjWZwbDi/eK0k6PpSlApEmiAs/q4gTra0s6XqONNH+b55yd3s2bMREZYvX+7XOM5FdHgI358zjo3763l54wF/h6NU0NMEYfF3BbG+oo52p/FL/4O7uLg4Zs6cyc6dOykrC74H8yyYlMaE9AR++e8dNLb0d3QYpZQ7TRC4qgfArxXEmt21hNqFogzf9z90NX36dJKTk3njjTdobw+uk6zNJtx3fS7VDS38rz6/WqnzogmCTxKEPyuINXtqKUxPIDLM/zfr2e125s2bR11dHR988IG/w+m3wvQEFk5OY/HqvezW51crdc40QeBqXgL/VRDHT7Wx7cBxn4y/1FejR4+moKCAVatWcfSoV0c58YofzM0mIsTOfcu0w1qpc6UJAv9XEGv31uI0cHFW4CQIgDlz5hAREcHSpUs7v6NgMTg2nO9cPZaVu6p5fethf4ejVFDSBIH/K4g1u2uJCLVRmJ7gl+P3JCoqijlz5nDgwAHWr1/v73D67YvTM8gZFsd9y7brOE1KnQNNEHySIPxVQawur2FKZhJhIYH315Gfn09WVhYrVqygvr7e3+H0S4jdxi9vzqe6sYXf6x3WSvVb4J2R/MCfVzEdbWim7GhjwDUvdRARrr32WkSEV199Neja8yekJ/DFaRk8sXYfmyuDK8Ep5W+aIPBvE1PH40Vn+Pn+h7NJSEhg7ty5VFRUsHbtWn+H02//MWccg2PC+eFLW2l3BFdfilL+pAkC/3ZSf1heS1xECLnD431+7P4oLCxk3LhxrFixIuiuaoqLCOW+63MpOXSCv63c4+9wlAoamiDwbwXx4Z4apo1Kxm4Tnx+7P0SE+fPnEx4ezssvv9z5nQWLefnDuCY/hQffLqPsSIO/w1EqKGiCwH8VRGXdKSrrmgK2/6Gr6Oho5s+fz+HDh1mxYoW/w+m3+67PIzrczvdf2ILDGVx9KUr5gyYI/FdBfLi7Bgjs/oeusrOzKSoqYs2aNezatcvf4fTL4NhwfnZ9Lpsq61m8Kvger6qUr2mCwH8VxOryWgbHhpM1JManxz1fc+bMISUlhVdeeYXjx4/7O5x+uX7CcK4aP5TfvbmT8qPa1KTU2WiCwD8VhNNp+HB3DTNGJyMS2P0PXYWEhLBgwQIcDgcvvvhiUPVHiAj/fXMe0eEhfPu5TbS261VNSvVEEwT+qSBKDp2gprGVmWMG++yYnpScnMz8+fOprKzk7bff9nc4/TIkNoJf3pzPtgMn+OPbwdVMppQvaYLAPxXEB2Wu/odLxwRHB3V38vLymDp1KmvXrmXz5s3+Dqdf5uSmsGhKOn95fzcf7an1dzhKBSRNEPingvigrJrslFiGxEX47JjeMHv2bDIzM1m2bBkHDx70dzj98uPrcshIiuK7z2/meJOO1aRUV5og8H0Fcaq1neKKY8wcG5zNS+7sdjsLFy4kJiaG5557jsbG4Hn+QnR4CH9cNJEjJ5r5wQubg24YEaW8TRMEvh+s76M9dbQ6nEHdvOQuKiqKRYsW0dTUxDPPPENra6u/Q+qzwvQE7p2XzfLtR1i8usLf4SgVUDRB4PvB+laWVRMeYmNKpv8fL+opKSkp3HzzzRw8eJCXXnopqJ4fccclI5mdM5Rfvl7Kxv3H/B2OUgHDqwlCROaKyE4RKReRe7tZHi4iz1nLPxKRTGt+soi8KyKNIvInb8YIvm9iWrmrmotGJRMR6v/Hi3pSdnY28+bNY+fOnbzxxhtB02QjIvx2wQRS4iO45+mNHDsZPBWQUt7ktQQhInbgz8A8IAf4rIjkdFntDuCYMSYL+APwa2t+M/Bj4Hveis+dLzupD9Q3sbv6JDMHSPNSV1OnTmXatGmsX7+e1atX+zucPouPCuXhWyZR3djC15/+WEd9VQrvVhBTgXJjzB5jTCvwLHBDl3VuAJZYr18ArhQRMcacNMaswpUovM6XFcSqsmqAAdFB3ZPZs2eTl5fHihUrKC4u9nc4fVaQlsAvb8rnw921PPB6qb/DUcrvQry471Sg0u19FXBRT+sYY9pF5DiQDNR4Ma4z+LKCWLmrhqFx4YwJsuE1+kNEuPHGG2ltbeW1114jLCyMgoICf4fVJ5+anMb2gydYvHovOcPiWFiU7u+QlPIbb54Ruxs/omujdF/W6fkAIl8RkWIRKa6uru5XcO58VUG0OZysLKtm5pjBQTe8Rn/Z7XYWLFhAZmYmr7zyCqWlwfOL/EfXZHNxVjL/7+VtFFfU+TscpfzGmwmiCnD/+ZUGdL2TqnMdEQkB4oE+/480xvzdGFNkjCkaPPjcm2w6Kghvn7TX762jobmdq3KGevU4gSI0NJRFixaRmprKP//5T7Zv3+7vkPokxG7jz5+bRFpiJHc+Uczu6uC5t0MpT/JmglgPjBGRkSISBiwClnZZZylwq/V6AfCO8cOlLw6HA5vN5vUE8XbpUcJCbAPm/oe+CA8P5/Of/zzp6em8+OKLbNmyxd8h9UlCVBiP3z4Vuwi3PbaO6oYWf4eklM95LUEYY9qBe4DlQCnwvDFmu4jcLyLXW6s9CiSLSDnwXaDzUlgRqQD+B7hNRKq6uQLKY5xOp9ebl4wxvFV6mItHJxMV5s2un8ATHh7OLbfcQkZGBi+//DIff/yxv0PqkxHJUTx62xSqG1q4c8l6TrW2+zskpXzKq72yxpjXjTFjjTGjjTEPWPN+YoxZar1uNsYsNMZkGWOmGmP2uG2baYxJMsbEGGPSjDEl3oqzo4LwprKjjVTWNV0wzUtdhYWF8bnPfY6srCyWLVvGypUrg+I+icL0BP73s5PYeuA4X3liA81twTO0uVLnS++kxpUgvF1BvFVyBIArsy/MBAGf9EkUFBTw7rvv8vrrrwfFHddX5wzlNwsmsKq8hnue3kib3iOhLhCaIHA1MXm7glhReoSCtHhS4oN79NbzZbfbufHGG5kxYwbFxcU8//zzQTF204LJafz8hlzeLj3Cfzy/WZ9prS4ImiDwfh9EdUMLGyvrL+jqwZ2IcPXVVzNv3jx27drF4sWLqa+v93dYvfrC9EzunZfN0s0H+Y/nN+nd1mrA0wSB95uY3t1xFGPgqpwhXjtGMJo6dSq33HIL9fX1PPLII+zbt8/fIfXqq5eN5vtzxvHKpoN889mN+shSNaBpgsD7TUxvlR5heHwEOcPivHaMYDV69GjuvPNOIiMjWbJkCR9++GHAd15/fVYW/3XteF7fepi7n9KOazVwaYLAuxVEQ3MbK3dVc3XO0AF/9/S5GjRoEHfeeSfZ2dm89dZbPPfcczQ1Nfk7rLO689JR/PzGPN4uPcoXF6/j+Cl9Ip0aeDRB4N0KYvn2I7S0O7m+MNUr+x8oIiIiWLhwIXPmzKGsrIy//e1vVFRU+Duss/rCtAweXFTIxv3HWPi3DzlYH9hJTan+0gSBdyuIVzcdYERSFJNGJHhl/wOJiDBt2jRuv/127HY7S5Ys4a233qK9PXBvULuhMJUlX5rKofpmbnp4NdsOHPd3SEp5jCYIvFdBHG1oZnV5DTcUDtfmpX5IS0vjrrvuYtKkSXz44Yc88sgjHDhwwN9h9WjG6EH882vTsYmw4K8fsnRz1yHHlApOmiDwXgXxr82HcBq4oXC4x/c90IWFhTF//nw++9nP0tTUxKOPPsry5csD9p6J7JQ4lt5zCfmp8XzzmY388t+leq+ECnqaIPBeBfHqpgPkDo8ja0isx/d9oRg7dix33303kyZNYu3atTz88MOUlpYG5JVOg2PDeerOadxy0Qj+9v4evvDoRxw94ZNnXinlFZog8E4FsbfmJJurjnOjdk6ft4iICK677jpuu+02wsPDef7553nqqaeoqfHpc6X6JCzExgM35fObBQVs3F/PvAc/4L2dR/0dllLnRBME3hms75WNBxCB+RO0eclTMjIyuOuuu5g7dy5VVVU8/PDD/Otf/6KxMfCe1/DponSWfeNiBsWEc9tj67lv2XaaWvV+CRVcLqxxp3vg6aE2nE7DK5sOMG1k8gU/9pKn2Ww2LrroIvLy8nj//ffZsGEDW7ZsYdq0aUyfPp3IyEh/h9gpa0gsr95zMb98vZTHVlfw3s5qfruggKLMJH+HplSfaAWB55uY3tt1lH21p1g0VZ9n7C3R0dFcc8013H333YwZM4YPPviABx98kHfffTegbrKLCLVz3w15PP3li2hzOFn4tzX89NVtHG/SG+tU4NMKAs93Uv/jg70Mi4/gmvxhHtun6l5ycjILFy7kyJEjvP/++6xcuZK1a9cyadIkpk2bRnx8vL9DBFyXwi7/9kx+88YOnly7j9e2Hub/XZvNjYWpegm0ClhaQeDZPoiSgyf4cHctX5yeSahdv15fGTp0KJ/+9Kf56le/yrhx4/joo4946KGHeOmll6iqqgqIq56iw0O474Y8lt5zCamJkXznuc0s+OsaNuzr82PYlfIprSDwbB/E4tV7iQy187mpIzyyP9U/Q4cO5eabb+aKK65g7dq1bNy4ka1btzJs2DCmTJlCbm4uYWFhfo0xLzWel782g+eLK/mft3bxqb+sYW5uCt+bM46sITF+jU0pd5og8FwFcbShmaWbDvKZKenER4V6IDJ1rhISEpg7dy6zZs1iy5YtrF+/nqVLl/LGG2+Qm5vLxIkTSUtL81vzjs0mLJo6gusLh/PIyr38beVulpccZn7BcL5xRRZjhuq9M8r/NEHguQri/9bup9Xh5PaLM88/KOUR4eHhTJkyhaKiIvbv38+mTZvYtm0bGzduJDExkby8PPLy8hgyxD/P6ogKC+FbV43hlmkj+McHe3liTQXLthxkTk4Kd146kskZidpHofxGEwSeqSCON7Xxf2v3cWX2EEYN1maCQCMiZGRkkJGRwdy5cyktLWXbtm2sWrWKDz74gEGDBjF+/HjGjx9PSkqKz0/Kg2LCuXdeNl+ZOYrFq/by5Np9vLH9MBPSE7h9RiZz81KICPXuc9OV6uqCTxDGGIwx511B/OGtXRw71cp3rh7rociUt4SHh1NYWEhhYSGNjY2UlJRQWlramSxiY2MZM2YMY8aMYeTIkYSHh/sstqToML43Zxx3zxrNix8fYPGqvXz7uU0kLAvlU5PS+MyUdMZq85PykQs+QTidrkdGnk8FUXLwBE+sqeCWi0aQlxoYl1WqvomJiWHq1KlMnTqVU6dOsXPnTsrKyti2bRsff/wxNpuNtLQ0Ro0aRWZmJqmpqYSEeP+/TVRYCF+YlsEtU0ewZk8tT6/bzxNrKnh01V5yhsVx08RUrpswjGHxgXNjoBp4LvgE4XC4hj841wrCGMNPl24jPjKU780e58nQlI9FRUUxceJEJk6ciMPhYP/+/ezevZu9e/fy3nvvARASEkJaWhrp6emkp6eTlpbm1bu3bTbh4qxBXJw1iJrGFv61+SAvbzrIA6+X8sDrpRSmJzA3L4Wrc4YyalC09lcoj9IEYSWIc60gXt54gPUVx/jVzfkkRPn38knlOXa7nZEjRzJy5EgAmpqa2LdvX+e0atWqznsrkpKSSE1NZfjw4QwbNoyhQ4cSEeH5IVYGxYRz28Ujue3ikeytOcnrWw+xfPthfvXvHfzq3zsYkRTF5eMGM3PMYKaOSiIuQq+kU+fngk8QHU1M51JBVNad4oHXSpmQnsCni3RYjYEsMjKS7OxssrOzAWhtbeXgwYNUVlZy8OBBKioq2Lp1a+f6CQkJDB06lCFDhjBkyBAGDx5MUlISoaGeOWmPHBTN12dl8fVZWRyob+KdHUd5f+dR/llcxRNr9mETyE+N56JRyUzOSGTSiEQGx/quL0UNDBd8gjjXJqaaxha+uHgd7U7D7xYUYLNpaX8hCQsLIzMzk8zMzM55jY2NHDp0iMOHD3P48GGOHj3Krl27TruLOzExkeTkZJKSkjqnxMREEhISzrlvIzUhki9My+AL0zJoaXfw8b561uypZc3uGh5fXcHfV+4BID0pkoLUBPLT4skbHk/2sFgGxWjSUD274BPEuXRSN7a0c/tj6zl0vImn7rxIb2pSgKvDu+Pqpw7t7e3U1NScNtXV1bF///4zno4XExNDfHx85xQXF0dcXByxsbHExsYSExPTaxIJD7EzfXQy00cnw9VjaWl3sO3ACT7ed4yNlcfYcqCe17Ye6lx/UEw441JiyBocw+ghMYweHENGchTD4iOx64+eC94FnyD6W0EcOdHMt57dSMmhEzzyxclMztChm1XPQkJCSElJISUl5bT5xhhOnjzJsWPHOqf6+npOnDjB4cOH2bVrF+3t7WfsLyIigpiYGGJiYoiOjiYqKqrzz44pMjKycwoLDWVyRiKTMxI793HsZCslh06w43ADOw+fYOfhBl78+ACNLZ8cL8xuIy0pkrTEKFITIklLjGRYfAQpcRGkxEcwJC6CmPAL/vQx4Hn1b1hE5gIPAnbgH8aYX3VZHg48AUwGaoHPGGMqrGU/BO4AHMA3jTHLvRFjXyuIdoeTJWv28Ye3dtHqcPK7hQVckT3UGyGpC4CIdJ7o09PP7L8yxtDU1MSJEydoaGjonBobGzl58mRnc9bJkydpaWnp8Tg2m42IiIjTpvDwcMLDw0kND2dkUhg3DAsnNDSRRkcIR5vh6CnD4UYHhxraOHLiFFur6jl26szhyaPC7AyODSc5Ooyk6HAGxYSREBVGYlQoiVFhxEWGEh8ZSlxkCHERocRGhBATHkKIDmIZNLyWIETEDvwZuBqoAtaLyFJjTInbancAx4wxWSKyCPg18BkRyQEWAbnAcOBtERlrjPH4I7nOVkE0trRTXFHHmj21vF1yhN3VJ7ls7GDuvyGXjORoT4eiVCcR6awIulYfXTkcDk6dOsWpU6doamri1KlTNDc309TURFNTE83NzZ1TS0sLDQ0NNDc309ra2mNyCQcyrQmgPcLGSRPKKRNGsy2SFgmjiTCaT4VRezKEKqedUw4bTQ4bDs7eNBVqg4gQiAwRwu1CRIiNiFAh3G4jIkSICLURHmIjLMRGhPVnmN31Z3iIjVDrdViInVC7EGa3E2IXQu02QjvmhdgJsdsIsdmw24UQm816L9jd5tttNuw213ybtY5NXIlVRDovG+7654XCmxXEVKDcGLMHQESeBW4A3BPEDcDPrNcvAH8S19/ADcCzxpgWYK+IlFv7W+PpIMuOnuSl5lzefqWS0NePYIwrMZxoaqOl3VVdhNqFwvQEvj9nHHNyfT8Mg1JnY7fbO/sp+ssYQ1tbG62tradNHfPa2tq6nRwOB+3t7W5TC+3t7bS1tdPU7qSx1cmpNjjV7pqaHdDsFFocQqux0eoQ2hx22o2Nk9ipx4bD2GjHNTmM4MCGAxuml4TjLYLp/FPOmMcZy5CO+e7rd92n+7qm+2XdbNc5v4cFk4aG8dg3r+1hq3PnzQSRClS6va8CLuppHWNMu4gcB5Kt+Wu7bJva9QAi8hXgKwAjRpzb8NqxURFkJIQyaFA0kZGR2ESICQ/pLIsL0uKZnJFIVJi2t6qBR0QICwvz+RDoxhgcDgdOpxOn03nGa2MMTqfTlXwchua2dprbHLS0OWhzOGlpd9Da5qDNaWhpc9DuNLQ7DW3trteOzvdOnMb13uE0ONxeO43r8cAOYzCGzvUMrtdOJxgrVqc5/TVY89zmd7zu+Hyu+afPw5rn2sPp30fnvO6WY06b0fXpJiOHeOdCGW+e9brLdV0/V0/r9GVbjDF/B/4OUFRUdE5PhCnMSmXpf56Re5RSXiQiPhmyRJ0fb/YWVQHuWibHcgAAB/lJREFUvW9pwMGe1hGRECAeqOvjtkoppbzImwliPTBGREaKSBiuTuelXdZZCtxqvV4AvGNctdZSYJGIhIvISGAMsM6LsSqllOrCazWe1adwD7Ac12Wui40x20XkfqDYGLMUeBR40uqErsOVRLDWex5Xh3Y78HVvXMGklFKqZxIID3P3hKKiIlNcXOzvMJRSKqiIyAZjTFF3y/SOFaWUUt3SBKGUUqpbmiCUUkp1SxOEUkqpbg2YTmoRqQb2nccuBgE1HgonGFxonxf0M18o9DP3T4YxZnB3CwZMgjhfIlLcU0/+QHShfV7Qz3yh0M/sOdrEpJRSqluaIJRSSnVLE8Qn/u7vAHzsQvu8oJ/5QqGf2UO0D0IppVS3tIJQSinVLU0QSimlunXBJwgRmSsiO0WkXETu9Xc83iYi6SLyroiUish2EfmWv2PyFRGxi8hGEfmXv2PxBRFJEJEXRGSH9fc93d8xeZuIfMf6d71NRJ4RkQh/x+RpIrJYRI6KyDa3eUki8paIlFl/JnriWBd0ghARO/BnYB6QA3xWRHL8G5XXtQP/YYwZD0wDvn4BfOYO3wJK/R2EDz0IvGGMyQYmMMA/u4ikAt8EiowxebgeM7DIv1F5xePA3C7z7gVWGGPGACus9+ftgk4QwFSg3BizxxjTCjwL3ODnmLzKGHPIGPOx9boB10ljwD9zVUTSgGuBf/g7Fl8QkThgJq5nrmCMaTXG1Ps3Kp8IASKtJ1RGMQCfRGmMWYnr+TnubgCWWK+XADd64lgXeoJIBSrd3ldxAZwsO4hIJjAR+Mi/kfjEH4EfAE5/B+Ijo4Bq4DGrWe0fIhLt76C8yRhzAPgdsJ//3969h0hZhXEc//5yDdkMwywpzDTQpCisLKLMNrWI6KJIBF2wCwTRhYIgCiorshsI3SwSLC2TysSMAq1UttJyLVclyn/arA3L+qMrabX9+uOcgclmp9V2ZpyZ5wPi65l5zzmj6/u855x5nwPbgR9tr6xtr6pmuO3tkG4CgUP7o9JmDxAqUdYU3/uVNBh4FbjZ9k+17k8lSTof2GH7o1r3pYpagBOBp2yfAPxKP0077KvyvPtFwGjgcOAASZfXtlf1rdkDRDdwRNGfR9CAQ9LdSRpICg6LbC+tdX+q4HTgQklfkKYRJ0t6obZdqrhuoNt2YXS4hBQwGtlUoMv2d7b/AJYCp9W4T9XyraTDAPLvO/qj0mYPEB3AGEmjJe1PWtBaXuM+VZQkkealP7U9p9b9qQbbt9seYXsU6d94le2GvrO0/Q3wlaSjc9EU0h7vjexL4FRJrfnnfAoNvjBfZDkwMx/PBF7rj0pb+qOSemX7T0k3ACtI33iYb/uTGner0k4HrgC2SOrMZXfYfrOGfQqVcSOwKN/8fA5cVeP+VJTtDyUtAT4mfVtvIw2YdkPSYqANGCapG7gbeBB4WdI1pEB5cb+0Fak2QgghlNLsU0whhBB6EQEihBBCSREgQgghlBQBIoQQQkkRIEIIIZQUASLUNUnDJb0o6XNJH0laJ2n6XtY1qjhDZiVJmiDpsWq0VdTmeEnnVbPNUN8iQIS6lR+GWga02z7K9kmkB+FG1LZn5Ulqsb3B9k2VqLvMy+OBCBChzyJAhHo2Gfjd9tOFAtvbbD8OIGmQpGclbckJ687K5aMkvSvp4/yrbDoGSYdJapfUmfcZOCOXn5vP3yTpnVw2VNIySZslfSDp+Fw+S9IzklYCCyW1FfalyK/Nl7Qmj4RuKmr7zryfw1t5f4NbS/TvOUlzJK0GHpJ0iqS1+TOvlXR0fljuXuCS/DkukXRAbrcjv7ehMxmHPdfUT1KHuncs6anZ3lwPYPs4SeOAlZLGkvLUnG17p6QxwGJgQpl6LgVW2L4/7yHSKukQYB4wyXaXpKH5vfcAG21PkzQZWEi6cwc4CZho+zdJbbu1MQ44CzgQ2CrpKdIeDjNIGXdb8mftLeHgWGCq7Z5Cqu+cKWAqMNv2DEl3kfZKuAFA0mxS2pGrJR0ErJf0tu1fy/xdhCYSASI0DElPAhNJo4qT8/HjALY/k7SNdCHdBjwhaTzQk8vK6QDm5ySHy2x35gt8u+2uXH8hP/9E0kUd26skHSxpSH5tue3femnjDdu7gF2SdgDDc12vFc6R9HqZPr5iuycfDwEW5OBnYGAv55xDSmJYGJUMAkbSPPmLwn+IABHq2SfkizGA7eslDQM25KJS6dwBbgG+Jd2h7wfsLNeI7XZJk0gbDj0v6RHgB0qnhi+XQr7cnfmuouMe0v/N3vpfSnHd9wGrbU9X2vNjTS/nCJhhe+setBOaSKxBhHq2Chgk6bqistai43bgMoA8tTQS2Eq6w95u+y9S4sIB5RqRdCRpP4l5pEy4JwLrgDMljc7vKUwxFbfZBnz/P/bbeA+4IK+lDCYFqL4YAnydj68sKv+ZNIVVsAK4MS/2I+mEvexnaFARIELdcso0OY10oe6StJ603eJt+S1zgQGStgAvAVfmaZy5wExJH5Cml/5rzr0N6JS0kTRiedT2d8C1wFJJm3L9ALOACZI2kzJszvx3dX3+fB2kNM6bSHsbbAB+7MOpDwMPSHqffwa/1cAxhUVq0khjILA5f733vr3ta2hMkc01hH2YpMG2f5HUShqdXFvYUzyESos1iBD2bc9IOoa0gLwggkOophhBhBBCKCnWIEIIIZQUASKEEEJJESBCCCGUFAEihBBCSREgQgghlPQ3dSTsvJ7p5K0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prior.plot(label='prior', color='gray')\n",
    "france.plot(label='France', color='C0')\n",
    "\n",
    "plt.xlabel('Goal scoring rate')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Posterior distribution for France')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data, `k=4`, makes us think higher values of λ are more likely and lower values are less likely.  So the posterior distribution is shifted to the right.\n",
    "\n",
    "Remember that the prior mean was about 1.4; here's the posterior mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.4140818156118378, 2.699772393342308)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_mean(prior), pmf_mean(france)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior mean is near 2.7, which is between the prior (1.4) and the data (4).\n",
    "\n",
    "Having seen 4 goals, we think France's goal-scoring rate is higher than 1.4, but we have not seen enough data to go all the way to 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Here's a prior distribution for Croatia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "croatia = prior.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Update `croatia` with the observed number of goals, `k=2`.  Plot the results and compute the posterior mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1609321178598705"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "k = 2\n",
    "λs = prior.index\n",
    "likelihood = poisson.pmf(k, λs)\n",
    "bayes_update(croatia, likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3xUZdb4v2cmvRdKIAkE6QkJPSIC0kEUFcWVFd8Ve1l3331dfX/qNnVd3WrbXeva1tdesK90RIpAkCY1CS2QQBoJSUif5/fH3OAwTDIzmTuTgM/387mf3Ln3KefemdxzzznPcx5RSqHRaDQajTOWjhZAo9FoNJ0TrSA0Go1G4xKtIDQajUbjEq0gNBqNRuMSrSA0Go1G4xKtIDQajUbjEq0gfgCISLWInOentieKyGGHzztEZKJJbc8XkcUOn5WI9DOjbaM9v92XVvoTEXlFRI6LyIZA9XuuIiJpxm8iqKNlaY1A/8bMRisIExGRAyJSa/wojhkPgygf2jPlH0ApFaWU2udLG170laGUWtlWGU+vSyn1hlJquhlyichKEbnZqf2A3ReDccA0IEUple1rYw73sdph2+q7mOcOInKtiOQY96ZIRP4jIuP81Fdn+I2ZilYQ5jNbKRUFjABGA7/uKEF8VSwd+WbWmd8KfaA3cEApVeNtRTf3I854EEUppYa2o/45iYjcDTwJPAp0B3oBzwCXt1L+B3eP3KKU0ptJG3AAmOrw+S/AZ8Z+T+AToBzIA25xKJcN5AAngGPA48bxQ4ACqo3tAuP4jcAu4DiwCOjt0JYCfgrkAvsdjvUz9mOBfwMlwEHsCsxinFsArAGeMOR8xMU1hgOvGn3vBO4FDru6B95cl6u+jWOrna7t58A+oNS4vy2yPwj8n0PZNKN8EPAHoBmoM/r7Rzvuy2rgr8Z17wcuduhrgSFTlXFuvov7dpPRf7Mhw0PG8Vuw/x7Ksf8+erb1XTq1eeoaXZxzdT/7AsuBMuP+vYFduTh+d/cA24BK4B0gzOH85cAW4/vMB2Y63LuXgCLgiNGXtZX/kWxgHVBhlP8HEOJ0zbcb13wc+Ccgxjmr8R2UGvf7p21cf6xxn69u4//1QeB94P+Ma7oZCMWuVAqN7Ukg1CgfD3xm/EaOG/spxjlPfmOXAJuNvgqABzv6meX2mdbRApxLG6c/HFOBHcDvjc9fYX97CQOGGT+yKca5dcB/GftRwBhj/4wHAHAF9gfKYOwPv18Dax3OK2AJkACEOxxr+ZH+G/gYiDba3wvcZJxbADQBPzPaDndxjX8EvjbaTwW+o3UF4c11ndE3rhXECqPvXobsNxvnHqQVBWF8XtlS1qk9T+9LI/aHuRW4A/vDQ4BI7P/wA42yPYCMVn4fztczGfvDbgT2B9PfgVVtfZdO7Z1xH93cz37YXVyhQFdgFfCk03e3AfvLTAL2l5DbjXPZ2JXGNOyeh2RgkHHuI+B54150M9q4rZV7MBIYY8iUZvTxC6dr/gyIM77jEr5XRLcDu7H/7hKM30Jr1z/TuP4zzjmUedD4Xq8wrikceBj4xriOrsBavv8fTgSuAiKM38l7wEcO7a2k7d/YRCDT6CsL+0vTFR393GrzmdbRApxLm/EPVo397eggdoUQbvygm4Foh7KPAa8a+6uAh4AuTu2d8QAA/oPx4DI+W4CTGFaEUX6yUzvKeDhYgXog3eHcbcBKY38BcMjNNe5r+Yc1Pt9K6wrCm+s6o29cKwjHvu8Elhn7D9JOBeHhfclzOBdh1E3C/lCswP7gOOMh7uZ6XgL+7PA5CvsDK62177KV+1jhsN3jxXd5BbDZ6bu7zuHzn4HnjP3ngSdctNHduHfhDsd+DKzw8H/mF8BCp+9knMPnd4H7jP3lGArL+Dzd+XfkcG4+cNRN3w/ioJCNY/nALIfPM7C7BV3VHwYcd/jc6m+slfpPurqnnWnTMQjzuUIpFaeU6q2UulMpVYv9jaxcKVXlUO4g9rcwsLsfBgC7RWSjiFzaRvu9gadEpEJEKrC7D8ShLbCbr67oAoQYfbuSo626LfR0KnOwtYJ4d12e9O1c5qAhj694cl+OtuwopU4au1HKHk+4BvvbbZGIfC4igzzst6djn0qpauzuH2++D7Ar4Dhj+2trdUWkm4i8LSJHROQEdtdKF6e2jjrsn8SutMD+kpPvou/eQDD2a2/5TT6P/Q38DERkgIh8JiJHDRke9UIGb357ZUAXD+IKzvf3tO8Eh9+YiESIyPMictCQfRUQJyJWN31g1D9fRFaISImIVGL/zThfe6dCK4jAUAgkiEi0w7Fe2P21KKVylVI/xv5P9SfgfRGJxP724UwBdvM9zmELV0qtdSjjqh7Y3RmN2P+pz5DDTd0WirA/LBzru8TL6/Kkb1z0XWjs12B/s28hyYu2PbkvraKUWqSUmobdvbQbeNGTethlP9WncW8S8e77aFM0p8+PGceylFIxwHXYXy48oQB7DMPV8XpOV1IxSqmMVtp5Fvs96m/I8IAXMnj828Pu3qzDbiW1hfM9Ou074fTf2C+BgcD5huwTjOMt8rv7rt7EHmdKVUrFAs/h+bV3CFpBBAClVAF2X+ZjIhImIlnY367fABCR60Skq1LKht1NAHaXVAlgAxzHUT8H3C8iGUbdWBG52kM5mrGb7H8QkWgR6Q3cjf1N0lPeNfqPF5EU7D5ul3h5XZ5yr9F3KvDf2AOpYA+eThCRXiISC9zvVO9Ya/35cl9EpLuIXGY83OuxuxibPbyWN4EbRGSYiIRif5ter5Q64GF9b4k25KsQkWTsAww85SXssk4REYuIJIvIIKVUEbAY+JuIxBjn+orIRW3IcAKoNiytO7yQ4V3g5yKSIiLxwH2tFVRKVQK/Bf4pIlcYb//BInKxiPy5jT7eAn4tIl1FpIvRRsvvIBqoxX7/EoDfOdVt9TfmUL9cKVUnItnAtW2U7RRoBRE4fozdZ1wILAR+p5RaYpybCewQkWrgKWCeUqrOcGX8AVhjmO9jlFILsb+Nv22Yud8BF3shx8+wv23vwz4y503gZS/qP4Td7N6P/cHwehtlPb4uL/r/GNiEXSF8jv3BhXEv38E+AmcT9kCnI08Bc8U+Se1pF+22975YsL9ZFmJ3912EPTbiFqXUMuA3wAfY3477AvM8qdtOHsIeEK/Efu8+9LSiUmoDcAP2UVGV2AddtLxp/wS7i24n9tE972O3plxxD/YHYxV2S+udVsq54kXso/a2At+6k18p9Th2Rf9r7C8lBcBd2IPqrfEI9pF324DtRj+PGOeexB5TLMUeyP7Sqa6739idwMMiUoVd8bzblvydgZbhYxqNRqPRnIa2IDQajUbjEq0gNBqNRuMSrSA0Go1G4xKtIDQajUbjknMmOVWXLl1UWlpaR4uh0Wg0ZxWbNm0qVUp1dXXunFEQaWlp5OTkdLQYGo1Gc1YhIq3OSNcuJo1Go9G4RCsIjUaj0bhEKwiNRqPRuOSciUFoNBpNY2Mjhw8fpq6urqNF6XSEhYWRkpJCcHCwx3W0gtBoNOcMhw8fJjo6mrS0NEQ6daLUgKKUoqysjMOHD9OnTx+P62kXk0ajOWeoq6sjMTFRKwcnRITExESvLSutIDQazTmFVg6uac990QrCS5RSbN68mdLS0o4WRaPRaPyKjkF4gc1m49NPP2XLli1YrVamTp3K+eefr99YNBqN1/z2t79lwoQJTJ06taNFaRWtIDzEZrPx0UcfsX37dsaOHUtpaSmLFi1i7969XHnllURFRblvRKPRaIDm5mYefvhhr+tYrR4tf20a2sXkAUopPvzwQ7Zv387kyZOZNm0a8+bN49JLL6WgoICPP/4YvfCSRqMBOHDgAIMGDeL6668nKyuLuXPncvLkSdLS0nj44YcZN24c7733HgsWLOD9998HYNmyZQwfPpzMzExuvPFG6uvrAc6oE2i0BeEBxcXF7NixgwkTJjB+/HjAHvAZOXIkDQ0NLF68mNzcXAYMGNDBkmo0mha+/PJLjh49amqbSUlJzJw50225PXv28NJLL3HhhRdy44038swzzwD2uQirV68+JR/YR14tWLCAZcuWMWDAAH7yk5/w7LPP8otf/OKMOoFGWxAesH//fgBGjBhxxrns7GwSExNZtGgRTU1NgRZNo9F0QlJTU7nwwgsBuO6660494K+55pozyu7Zs4c+ffqcesG8/vrrWbVq1anzruoECm1BeMD+/ftJSEggNjb2jHNWq5WZM2fyxhtvsH79+lM/Co1G07F48qbvL5wHrrR8joyMPKOsO/e0qzqBQlsQbrDZbBw4cKDN2Yf9+vVjwIABrFq1iqqqqgBKp9FoOiOHDh1i3bp1ALz11luMGzeu1bKDBg3iwIED5OXlAfD6669z0UUXBUROd2gF4YbCwkIaGhrcTk+fMWMGzc3Np5mGGo3mh8ngwYN57bXXyMrKory8nDvuuKPVsmFhYbzyyitcffXVZGZmYrFYuP322wMobev41cUkIjOBpwAr8C+l1B+dzocC/wZGAmXANUqpAyKSBuwC9hhFv1FKdcgda4k/uFutLiEhgczMTLZu3cqUKVMICwsLgHQajaYzYrFYeO655047duDAgdM+v/rqq6f2p0yZwubNm89ox7lOoPGbBSEiVuCfwMVAOvBjEUl3KnYTcFwp1Q94AviTw7l8pdQwY+swdbp//366d+/ukR8wOzubxsZGtmzZEgDJNBqNxr/408WUDeQppfYppRqAt4HLncpcDrxm7L8PTJFONC25qamJgoICj7Mf9ujRg9TUVDZu3KjnRWg0P1DS0tL47rvvOloMU/CngkgGChw+HzaOuSyjlGoCKoFE41wfEdksIl+JyHhXHYjIrSKSIyI5JSUl5koPFBQU0NTU5FV63NGjR1NeXn4q4OQtqrmZys8+5/jbb2M7ebJdbWg0Go0Z+FNBuLIEnF+rWytTBPRSSg0H7gbeFJGYMwoq9YJSapRSalTXrl19FtiZ/fv3IyL07t3b4zrp6elERUWxYcMGr/pSSlG1YgX7r5hD4T33cPTBh8ibPIWSZ56h+cQJb0XXaDQan/GngjgMpDp8TgEKWysjIkFALFCulKpXSpUBKKU2AflAwKcpHzhwgOTkZEJDQz2uY7VaGTlyJHl5eZSVlXlcr+jXv+bwHXdia6gn+ckn6P3Wm4QPG0bp039n/5VX0VxR0Z5L0Gg0mnbjTwWxEegvIn1EJASYB3ziVOYT4Hpjfy6wXCmlRKSrEeRGRM4D+gP7/CjrGTQ0NHDkyBG3o5dcMWrUKCwWCzk5OR6Vr/lmPZUffEj8T/6Lvp99RszMmUQMH07qc8/S69+v0XjsGIX3P6DjGhqNJqD4TUEYMYW7gEXYh6y+q5TaISIPi8hlRrGXgEQRycPuSrrPOD4B2CYiW7EHr29XSpX7S1ZXlJeXY7PZ6NGjh9d1o6KiGDhwINu3b8dms7VZVjU1cezRRwlOTqbb3XcjTuvFRmZn0/3ee6lesYLyl1/xWhaNRhNYrFYrw4YNO7V19FBVX/DrPAil1BfAF07HfuuwXwdc7aLeB8AH/pTNHS0zoqOjo9tVPysri127dpGfn0///v1bLXf8nXeo37uX5KeewtLK3In4/7qOkzk5FD/+OOHDhxMxYni7ZNJoNP4nPDy8zaHuTU1NBAWdHVmO9EzqVvBVQfTv35/w8HC2bt3aapmm48cpefrvRIwZQ/T0aa2WExF6/OERgnv2pPDee7E1NLRLJo1G0zG8+uqrXH311cyePZvp06dTXV3NlClTGDFiBJmZmXz88ceAPe45ePBgbrnlFjIyMpg+fTq1tbUA5OXlMXXqVIYOHcqIESPIz88H4C9/+QujR48mKyuL3/3ud6bKfXaosQ6gRUG0dyEgq9XKkCFD2Lx5M3V1dS5nVpf+8xls1dV0f+B+t6vSWaOjSfrd7yi4+WYqP1xI/LyOy/Co0ZwNPPTpDnYWmjsCML1nDL+bndFmmdraWoYNGwZAnz59WLhwIQDr1q1j27ZtJCQk0NTUxMKFC4mJiaG0tJQxY8Zw2WV2z3tubi5vvfUWL774Ij/60Y/44IMPuO6665g/fz733Xcfc+bMoa6uDpvNdmqpgQ0bNqCU4rLLLmPVqlVMmDDBlOvVCqIVqqqqiIiI8MkUzMrKYuPGjezcufOMVOG2+noqFy4k9tJLCPNwHYnIC8cSPnQopS88T9yVc5CQkHbLptFo/ENrLqZp06aRkJAA2Ie1P/DAA6xatQqLxcKRI0c4duwYYFcqLQpm5MiRHDhwgKqqKo4cOcKcOXMATr1wLl68mMWLFzN8uN3tXF1dTW5urlYQ/qaqqqrd7qUWkpOTSUhIYNu2bWcoiJo1a7HV1BBzySUetycidLnrLgpuuYWKhR8Rf82PfJJPozmXcfemH2gc0/W88cYblJSUsGnTJoKDg0lLS6Ourg7gtGH1VquV2traVkcwKqW4//77ue222/wis45BtIIZCkJEGDp0KAcPHqTCaR5D1aIvscTEEDlmjFdtRo670G5FPP8cSsciNJqzksrKSrp160ZwcDArVqzg4MGDbZaPiYkhJSWFjz76CID6+npOnjzJjBkzePnll6murgbgyJEjFBcXmyanVhCtYIaCALubCWDbtm2njtkaGqhavoLoKVO8dhPZrYif0lRYRMXCj3yWT6PRBJ758+eTk5PDqFGjeOONNxg0aJDbOq+//jpPP/00WVlZjB07lqNHjzJ9+nSuvfZaLrjgAjIzM5k7d66pa9LIuTL5atSoUcrTiWnusNls/P73v2fChAlMmjTJ5/ZeeeUVamtrufPOOwGoWrmSw7ffQcpzzxI9caLX7SmlODBvHs0lpfRdshixWn2WUaM5F9i1axeDBw/uaDE6La7uj4hsUkqNclVeWxAuaDHXzLAgADIyMigpKTll+lUtWowlOprIsWPb1Z6IkHjjTTQWFlLTQYuZazSacx+tIFzg6xwIZ1o09o4dO1ANDVQtW0b05ElYfBiFFD1pItbERI6/+54pMmo0Go0zWkG4wGwFER0dTVpaGjt37qR63TpsJ04QPcO3BdUlJIS4K+dQvXIljcfMC0ppNBpNC1pBuMBsBQH2NOClpaUUf/wJlshIIsdd6HObcXPnQnMzlQs/NEFCjUajOR2tIFxQVVWFiHi0zKinpKenI0DtmjVEXTTBJ/dSCyG9exMxZgwV772PcpMUUKPRaLxFKwgXVFVVERUVhcVi3u2JjIxkQEIC1spKwkeMNK3duKvn0njkCDVr15nWpkaj0YBWEC4xaw6EMy0jnWt6pbZZzhuip03DGhdHxXs6WK3RdBaOHj3KvHnz6Nu3L+np6cyaNYu9e/f63O6jjz562uex7RwJ6SlaQbjAXwqiS2kZTUFB7DGyM5qBJSSE2CuuoGrZMpqOHzetXY1G0z6UUsyZM4eJEyeSn5/Pzp07efTRR0/lWgJobm5uV9vOCmLt2rU+yeoOrSBc4C8F0bh9OyeTk9m5Z4+pq8PFXjYbmpqoWrLEtDY1Gk37WLFiBcHBwdx+++2njg0bNozm5mYmTZrEtddeS2ZmJgCPP/44Q4YMYciQITz55JOnyl9xxRWMHDmSjIwMXnjhBQDuu+++U5li58+fD3yfbbq19OG+opP1OdHU1ERtba3pCsJ28iR1e/YQfumlHD9+nGPHjpGUlGRK26GDBxOSlsaJL/5D/I90Aj+NBuDoo49Sv2u3qW2GDh5E0gMPtFnmu+++Y+RI13HGDRs28N1339GnTx82bdrEK6+8wvr161FKcf7553PRRRcxfPhwXn75ZRISEqitrWX06NFcddVV/PGPf+Qf//iHy0yxYWFhLtOHu1tGwB3agnDCH0NcAWq3fwfNzfScPAkRYefOnaa1LSLEzJrFyfXraTQxUZdGozGX7Oxs+vTpA8Dq1auZM2cOkZGRREVFceWVV/L1118D8PTTTzN06FDGjBlDQUEBubm5bbbbkj48KyuLqVOnnpY+3Be0BeGE3xTE5s0AxJ9/Pr2PHmXXrl1MnjzZtPZjLplF6TPPULVoMQn/dZ1p7Wo0Zyvu3vT9RUZGBu+//77Lc45D51tzM69cuZKlS5eybt06IiIimDhx4qlU4K3RVvpwX9AWhBP+VBAhfftijYs7NWmupKTEtPZD+/YldNAgTnzxhfvCGo3Gb0yePJn6+npefPHFU8c2btzIV199dVq5CRMm8NFHH3Hy5ElqampYuHAh48ePp7Kykvj4eCIiIti9ezfffPPNqTrBwcE0Njae0ae36cM9RSsIJ/yhIJTNRu2WLYQPt68S1ZLa10w3E0DMrFnUbt5M45Ejprar0Wg8R0RYuHAhS5YsoW/fvmRkZPDggw/Ss2fP08qNGDGCBQsWkJ2dzfnnn8/NN9/M8OHDmTlzJk1NTWRlZfGb3/yGMQ5rxtx6661kZWWdClK30J704R5di073fTpLlixh/fr1/OpXv/I5wNNC/b597Jt1CT0e+b09PQb2FOB1dXXccccdpvQB0HD4MPlTp9Htnl+SePPNprWr0Zwt6HTfbaPTfftIyxBXs5QDQO1m+6iDcGPdWLBneC0uLqasrMy0fkJSUggbmkWldjNpNBoT0ArCCX/MgajdshlLbCwhxugF+D4FuNlupthZs6jfuYsGk3yQGo3mh4tWEE74R0FsIXxoFuKQ2yk2Npbk5GR27dplal/RU6cCULV0qantajRnC+eK29xs2nNftIJwwmwFYWtooH7ffsLS0884N3jwYIqKiqioqDCtv+DkZMLS06laohWE5odHWFgYZWVlWkk4oZSirKyMsLAwr+rpeRAO1NfX09DQYKqCaNi/H5qbCe3f/4xz6enpLF26lF27dnHBBReY1mf0tKmUPPU0jcXFBHfrZlq7Gk1nJyUlhcOHD5s6hPxcISwsjJSUFK/qaAXhQMta1C35Tcygfq99BqQrBREfH09SUpL5CmKqXUFUL19O/Lx5prWr0XR2goODT81U1viOX11MIjJTRPaISJ6I3OfifKiIvGOcXy8iaU7ne4lItYjc4085W2iZgBJiwmI+LdTn5kJQEKFpaS7PDx48mIKCglPzL8wgpF8/Qnr3pmqxTt6n0Wjaj98UhIhYgX8CFwPpwI9FxNkRfxNwXCnVD3gC+JPT+SeA//hLRmdaUvCauVBQfW4uoX3SkFaUTstoJjOD1SJC9LSp1GzYQHNlpWntajSaHxb+tCCygTyl1D6lVAPwNnC5U5nLgdeM/feBKWJMQBCRK4B9wA4/yngaNmPZTqvValqb9Xl5Lt1LLXTt2pUuXbqYP5pp2jRoaqLaaXq/RqPReIo/FUQyUODw+bBxzGUZpVQTUAkkikgk8P+Ah/wo3xmYbUHYTp6ksaCgTQUBdivi4MGD1NTUmNIvQFhmJkHduuk1IjQaTbvxp4JwNRXZeexZa2UeAp5QSlW32YHIrSKSIyI5ZoxaMNuCqM/PB1wHqB0ZPHgwSin27NljSr8AYrEQPXUK1V+vxmbiCnYajeaHgz8VxGHAcfHlFKCwtTIiEgTEAuXA+cCfReQA8AvgARG5y7kDpdQLSqlRSqlRXbt29VngFgvCNAXRxggmR5KSkoiLizPdzRQ1ZQqqro6add+4L6zRaDRO+FNBbAT6i0gfEQkB5gGfOJX5BLje2J8LLFd2xiul0pRSacCTwKNKqX/4UVbgewvCLBdTfW4uEhZGsJuxxyJCeno6+/bto9bEt/3I0aOxREZSvWKFaW1qNJofDn5TEEZM4S5gEbALeFcptUNEHhaRy4xiL2GPOeQBdwNnDIUNJKZbELm5hPbti3jQXnp6Ojabjb1795rSN4CEhBA5fjxVK1egDOWn0Wg0nuLXiXJKqS+AL5yO/dZhvw642k0bD/pFOBf4w4KIHDvWo7I9e/YkJiaGnTt3MnToUFP6B4ieNJGqL7+kbscOwo2F0jUajcYTdC4mB8y0IJorKmgqLia0fz+Pyre4mfLz86mvr/e5/xYiJ0wAi4Wq5ctNa1Oj0fww0ArCATMtiPq8PMB9gNqR9PR0mpubTXUzBcXHEz5iONXLdRxCo9F4h1YQDphpQdTnejaCyZGUlBSio6NNXyMietJk6vfs0UuRajQar9AKwgFTLYjcXCxRUQQlJXlcR0QYPHgweXl5NDQ0+CxDC1GTJwFQtWJlm+Uam20cLKuhoUkHtDUajc7mehqmWhB7cwnt39/rpUvT09PZsGEDubm5ZGRk+CwHQGifPoSkpVG9YgUJ152+2PnRyjpeWbOfTQeP811hJXWNNkKsFgb3iGZoahzXj02jb1fzsttqNJqzB21BOGBmqo36/HxC+3kWoHYkNTWVyMhI091MUZMn25P3GSnNlVK8l1PAtCe+4pU1B1DAtdm9eXROJgsuTCM8xMp7OYeZ8cQqHvlsJ5W1jabKo9FoOj/agnDArFQbzSdO0Hz8OCGtpPhuC4vFwuDBg9m6dSsNDQ2mpR6PnjSR8pdfpmbNWprHTeTud7ewYk8J2X0S+PNVWaR1iTyjTklVPX9bvIeX1uxn4eYjPDVvOOP6dzFFHo1G0/nRFoQDZlkQDYfsOQpDevdqV/2MjAwaGxvJNQLdZhA+fDiW2Fgql6/gptc2sia/jAdnp/P2LWNcKgeArtGh/PGqLD69axxdokK54dUNfLrVOVuKRqM5V9EKwgGbzYbFYvE6buBM46GDAAT3ap+C6NWrl+luJgkKInLcOIqXLGProXKenjeMBRf2wWJxf61DkmN59/YLGJ4az8/f3syra/abJpdGo+m8aAXhQHNzsynxh4ZDhwAISU11U9I1FouF9PR09u7da+popsWx/Yk4WcVjg4OYOaSHV3Vjw4P5903ZTBvcnQc/3cnLq7WS0GjOdbSCcMBms5kygqnh4CGCunfHEh7e7jYyMjJoamoybdLcB5sO83hFAjaLhYnH25dWPCzYyjPzRzAjozuPfL6TFbuLTZFNo9F0TrSCcMBMCyKkne6lFlJTU4mKijLFzXS8poFHPt/JoP7JRI4cSfXK9q8yF2S18MQ1wxjcI4afvbWZPUfNW0tbo9F0LrSCcMA0C+LQQYLbGaBuocXNlJub63Nupj99uZuquib+MCeT6LtexM8AACAASURBVIkTqd+9m8bC9gebI0KC+Nf1o4gIsXLTaxspqzYvd5RGo+k8aAXhQEuQ2qc2ampoLiklpFdvn+Uxw82Uc6CctzcWcNO4PgxMiiZq0kQAn9eq7hEbzos/GUVxVT33vr8NpZwXC9RoNGc7WkE40Nzc7LMF0VBgDHH10cUEdjdTdHQ0O3bsaFf9xmYbv/7oO3rGhvHzKfacUCF9+hDcuxdVK1f6LN/Q1Djuv3gQy3cX8+aGQz63p9FoOhdaQThghgXRcNAYweSjiwm+TwGel5dHXV2d1/Xf2nCI3Uer+N1lGUSGBp1qM3riRE6u+wbbyZM+y3j9BWmM79+FRz7bxb6SNpcQ12g0ZxlaQThgigXRMgci1XcFAZCZmUlzc7PX61U3Ntt4/qt9jOodz/T07qedi5o0CdXQQM03vq9VbbEIf5k7lJAgC//z7lYam3WiP43mXEErCAfMGMXUeOgQ1i5dsEa5np3sLT179iQ+Pp7vvvvOq3ofbynkSEUtd07qe8bEv4gRI7BERZm2VnVSbBh/mDOErQUVvLBqnyltajSajkcrCAfMGMXUcND3Ia6OiAhDhgxh//79VFd75sKx2RTPfZXPoKRoJg3sdmabISFEjh9H9cqvTFur+tKsnszKTOLpZbkcKvPddaXRaDoerSAcMMOCMGMOhDNDhgxBKeXxnIjFO4+RV1zNHRPPtB5aiJ44kaaSEup2mJfO47eXZhBkEX77yXd6VJNGcw6gFYQDvloQtro6mo4eNSVA7Ui3bt3o1q2bR24mpRTPrsyjV0IEl2S2nk6jZa1qs9xMYHc1/c+0AazcU8KX3x01rV2NRtMxeKQgROQDEblERM5pheJrkLrRGOLa3iR9bTFkyBAKCgqoqKhos9za/DK2Hq7ktovOI8ja+tcVFB9P+PDhVJsw3NWRBWPTGNwjhoc+3Ul1fZOpbWs0msDi6QP/WeBaIFdE/igig/woU4fh6zDXU0n6TJgk58yQIUMA3FoRr687SEJkCFeNSHHbZtTEi6jbuZPGY8dMkRHsqTj+MGcIx6rq+Psy89KVazSawOPR01AptVQpNR8YARwAlojIWhG5QUSC/SlgIPHVgjg1B6JX+7K4tkV8fDwpKSltKojS6nqW7jrGlcOTCQt2fx3Rk+xrVVe7WavaW0b0imfuiBReWXNAB6w1mrMYj1+XRSQRWADcDGwGnsKuMJb4RbIOwHcL4iDWuDissbEmSvU9mZmZHDt2jOJi11lUF357hCab4prRnimokL59CU5JMd3NBHDPjIFYLcIfv/Ru/oZGo+k8eBqD+BD4GogAZiulLlNKvaOU+hlwzqxo73MM4tAhn5P0tUVGRgYWi4WtW7eecU4pxTs5BQzvFUf/7tEetSciRE2aRM26ddhqa02VtXtMGLdf1Jcvth9l44FyU9vWaDSBwdPX5X8ppdKVUo8ppYoARCQUQCk1ym/SBRifLYiDh/wSf2ghMjKSfv36sX379lPrZ7fw7aEK8oqruWaUd+6t6EkTUfX11KxbZ6aoANw64TySYsL4/Wc7sdn0sFeN5mzD06fhIy6Omf9E6WB8mQehGhtpLCoiJNV9cNgXsrKyqKqq4sCBA6cdf3djAREhVi4d2tOr9iJGjcISFUXV8uUmSmknPMTK/84cyLbDlXy89Yjp7Ws0Gv/S5tNQRJJEZCQQLiLDRWSEsU3E7m46p/DFxdR47BjYbAQnJ5ss1ekMGDCA0NBQtm3bdupYTX0Tn20r5JLMHkQZSfk8RUJCiJowgeoVK02bVe3IFcOSGZIcw98W76W+qdn09jUajf9w97o8A/grkAI8DvzN2O4GHnDXuIjMFJE9IpInIve5OB8qIu8Y59eLSJpxPFtEthjbVhGZ491ltQ9fXEyNR+wL8AT39O4N3luCg4NJT09n165dNDY2AvD59iJqGpo9Dk47EzV5Ms1lZdQ5KB2zsFiE/50xiMPHa3lrvU4JrtGcTbT5NFRKvaaUmgQsUEpNctguU0p92FZdEbEC/wQuBtKBH4tIulOxm4DjSql+wBPAn4zj3wGjlFLDgJnA8yLi3atxO/DJgjhid6H424IAu5upoaGB3bt3A/Dp1kLSEiMY2Tu+Xe1FTRgPQUFULTPfzQQwvn8XLjgvkb8vz9OT5zSaswh3LqbrjN00EbnbeXPTdjaQp5Tap5RqAN4GLncqcznwmrH/PjBFREQpdVIp1fIkCQMCEuH0yYIoLAQRgpOSTJbqTHr37k1sbCzbtm2jvKaBtfllzMrs0WreJXdYY2KIGD2KqhX+URAiwr0zB1JW08DLq/f7pQ+NRmM+7p6GLTmro4BoF1tbJAMFDp8PG8dcljEUQiWQCCAi54vIDmA7cLuDwjiFiNwqIjkiklNSUuJGnLZRSvmUi6nxyBGCunVDQkJ8ksMTRISsrCzy8/P5eNMBmm2KWW3kXfKE6EmTacjLp+HgQZOkPJ0RvezrUry4ah/lNQ1+6UOj0ZiLOxfT88bfh1xtbtp29TrrbAm0WkYptV4plQGMBu4XkTAX8r2glBqllBrVtWtXN+K0Tcuw0XYriMJCv8cfHBk6dChKKT7cuJ9eCRFk9Izxqb2oyZMBqFpuXvI+Z+6ZMZDqhiaeXZnntz40Go15tOnXF5Gn2zqvlPp5G6cPA45R0xSgsJUyh40YQyxw2qwqpdQuEakBhgA5bcnjCy0Kov1B6iOEDxtmpkhtkpiYSJfk3nyX18htF/Vqt3uphZCUZEIHDqR6+XISb1hgjpBODOgezZxhyfx73UFuGX8e3WLO0PkajaYT4e5puMnN1hYbgf4i0kdEQoB5wCdOZT4Brjf25wLLlVLKqBMEICK9gYHYc0D5jeZm+xDM9lgQqrmZxqNHA2pBAJyM749CGNnVnCS7UZMncXLTJpqOHzelPVf8fEp/mmyKZ1bm+60PjUZjDm1aEEqp19o676Zuk4jcBSwCrMDLSqkdIvIwkKOU+gR4CXhdRPKwWw7zjOrjgPtEpBGwAXcqpUrbK4sn+GJBNBUXQ1NTQEYwObLtuJVoqafm8G4YNdDn9qInT6Hs2eeoXvkVcXOuMEHCM0nrEslVI5J5c8MhbrvoPHrEhvulH41G4zvuRjE9afz9VEQ+cd7cNa6U+kIpNUAp1Vcp9Qfj2G8N5YBSqk4pdbVSqp9SKlsptc84/rpSKkMpNUwpNUIp9ZHvl9o2vlgQp4a4BtCCqDjZwLp95WT3CGbnzh00NPge+A0bkkFQUhJVS5eaIGHr/Gxyf2w2xT9X6FiERtOZcfe6/Lrx9698P0nOcTtn8MWCaCw0JskF0IJYvOOYPXPr2P40NjZ6vBxpW4gI0VOnUrN6NbaT/kvTnZoQwY9Gp/LOxgIOH9fpwDWazoq7UUybjL9fYc+9dBy7K2idceycwRwLwrehpt6weOcxkuPCmTpiAAkJCWzevNmUdqOnTkXV11O9erUp7bXGXZP6IQj/WK6tCI2ms+Jpuu9LgHzgaeAfQJ6IXOxPwQJNi4JorwVh7dIFS1hgRuXUNTazJq+UKYO7YbFYGD58OIcOHaK01PcwTcSokVjj4qha4l83U8+4cOZlp/L+psMUlGsrQqPpjHj6NPwbMEkpNVEpdREwCXtqjHMGX+ZBNB45EtD4wzf7yqhtbGbSoG4ADBs2DIvFwqZN7gaWuUeCgoiaNInqlStRJsQ12uKOiX2xiPCMnheh0XRKPFUQxUopx//ifYDrZc3OUnyyII4UEpwcOAWxfHcx4cFWLjgvEYCoqCgGDRrE1q1baWryPddR9LSp2KqqqNmw0ee22qJHrN2KeC/nsI5FaDSdEHejmK4UkSuBHSLyhYgsEJHrgU+xz3M4Z2ivBaFstoDOolZKsXx3MRf263LautMjR46ktrbWlGB15NixSEQEVUv9v5psixXxzxV6XoRG09lw97o829jCgGPARcBEoARoX+rQTkp7g9RNJaWoxsaAjWDKLa7m8PFaJhvupRb69OlDfHy8KW4mS1gYUePGUbVsmV/WiHDkeytCj2jSaDob7kYx3dDGdmOghAwE7R3m2lgY2DkQy3bZPXvOCkJEGDlyJIcOHcLXxIVgdzM1l5RSu2WLz225Q1sRGk3nxNNRTGEi8lMReUZEXm7Z/C1cIGmvBdGyUFBIgCyIFbuLSe8RQ1LsmSOmhg0bhtVqNcWKiJo4EQkOpmrRIp/bckeP2HCuGZ3K+5sKOFJR6/f+NBqNZ3j6uvw6kIR9hbmvsCfeq/KXUB1B+y2IwKwkB/bZ0zkHy5kyuJvL85GRkQwePJitW7eeWm2uvVijo4kcN44Tixb73c0EdisC4Bk9u1qj6TR4+jTsp5T6DVBj5Ge6BMj0n1iBp/0WxBGscXFYIiPdF/aRr/aWYFNnupccGTVqFHV1dWzfvt3n/mJmzqDp6FFqt271uS139IwL50ejUnk3p4BCbUVoNJ0CTxVEy+tohYgMwZ6WO80vEnUQvlgQgYo/rNxTQkJkCENT4lot06tXL7p3786GDRtQyreF+KImT7a7mb70v5sJ4M5J/QD0vAiNppPg6dPwBRGJB36DPUX3Tr5fP/qcwBcLIhAjmJRSfJ1byvj+XbBYWl/7QUTIzs7m2LFjHDp0yKc+v3czLQqImyk5Lpy5I1N5d+Nhiiq1FaHRdDQeKQil1L+UUseVUl8ppc5TSnVrWW3uXKE9FoRSKmAWxO6jVZRW1zOuXxe3ZTMzMwkLC2P9+vU+9xtINxPAnRP7YlOKZ/V6ERpNh+PpKKZEEfm7iHwrIptE5EkRSfS3cIGkPRZEc0UFqq4uILOov861D10d39/90qrBwcGMGDGC3bt3U1lZ6VO/gXYzpSZEMHdkCm9vKOBoZV1A+tRoNK7x9HX5beypNa7CvvJbKfCOv4TqCNqTaqNlBFNQD/9ncf06t5QB3aNcDm91xejRowHIyfFtldZTbqbFgRnNBPDTSf0MK0LHIjSajsTTp2GCUur3Sqn9xvYI0Hqk9CykPak2moqKAAhO8q+CqGtsZsP+csb1c289tBAXF8fAgQPZtGmTz0NeY2bOoKmoKGBuphYr4i1tRWg0HYqnCmKFiMwTEYux/Qj43J+CBZr2uJgai44C/l8HYuOBcuqbbIwf4D7+4Mj5559PbW0t27Zt86n/qMmTkZAQTnzxH5/a8YYWK+K5r3QsQqPpKNwl66sSkRPAbcCbQIOxvQ38j//FCxztCVI3FhUhISFYExL8JRYAq3NLCbFaOL+Pd/307t2bHj16sG7dOp+GvFqjo4maOJETX3yBMiFbrCekJkRw1YgU3txwiGMntBWh0XQE7nIxRSulYoy/FqVUkLFZlFIxgRIyELTPgigkqEcSIq0POzWDVbmljOwdT0RIkFf1RISxY8dSVlbGnj17fJIhZvalNJeVUfON7yOjPOWnk/phs+kRTRpNR+Hx67KIXCYifzW2S/0pVEdgs9kQEa8e9k1FRwnu4d8RTCVV9ewqOuG1e6mF9PR0YmNjWbdunU9yRE2YgCU6mhOffeZTO97QK/F7K0LHIjSawOPpMNc/Av+NfYLcTuC/jWPnDM3Nzd7Poi4qItjPI5jW5NmXER3vRYDaEYvFwpgxYzh06BCHDx9utxyW0FCip0+jaskSbHWBe1jfNdluRejZ1RpN4PH0iTgLmKaUelkp9TIw0zh2zmCz2bxyL6nGRpqKi/2uIFbllhAfEUxGz/Z79IYPH05YWJjPVkTspZdiq6mheuVKn9rxhtSECK4elcrbG3SmV40m0Hjzyuw4rDXWbEE6Gm8tiKbiYlCKoB5JfpNJKcW6/DLG9m07vYY7QkNDGTlyJLt27aKsrKzd7URkZxPUtSuVAXQzgd2KUCj+sVxbERpNIPH0ifgYsFlEXhWR14BNwKP+EyvweGtBNLbMgfBjDOJA2UmKKusY28/3SetjxozBarWyZs2adrchVisxs2ZR89Uqmn2coe0NyXHhzBvdi/dyCigo16vOaTSBwq2CEHvUdjUwBvjQ2C5QSr3tZ9kCircWxCkF4cc5EC3xh7F92xegdiQqKooRI0awdetWn9JvxFx6KaqxkRMBSr3Rwp2T+mKxCE8vyw1ovxrNDxm3T0RlH0D/kVKqSCn1iVLqY6XU0QDIFlC8tyCMSXJJ/nMxrcsvo0dsGGmJEaa0N3bsWACfrIiwIRmE9O1L5cKFpsjkKT1iw7k2uxcfbj7CvpLqgPat0fxQ8fSV+RsRGe1XSToY7y2IQqyxsVgizHl4O2OzKdbts8cfzJpnERsby7Bhw/j222+pqmrfgoAiQtyVc6jdsoX6fftNkctTfjqpHyFWC08s1VaERhMIPH0iTsKuJPJFZJuIbBcRt/kbRGSmiOwRkTwRuc/F+VARecc4v15E0ozj04yssduNv5O9uaj24K0F0VRYRJAf03zvPlpFeU0DY/uamzR33Lhx2Gw2n0Y0xcyeDVZrwK2IrtGh3HBhGp9uLWRX0YmA9q3R/BDxVEFcDJwHTAZmA5caf1tFRKzAP4266cCPRSTdqdhNwHGlVD/gCb5fhKgUmK2UygSux74mtl9pbm72zsV09Khf3Utr8+3xhwtMVhDx8fFkZmaSk5NDTU1Nu9oI7taNqHHjqPz4Y5QxAz1Q3DahL9FhQfxt8d6A9qvR/BBxl4spTER+AdyLfe7DEaXUwZbNTdvZQJ5Sap9SqiV/0+VOZS4HXjP23wemiIgopTYrpQqN4zuAMBEJ9eK6vMZms3kdpPbnHIh1+WX06RJJz7hw09seP348TU1NPsUiYq+8kqbiYmrWrjVRMg/6jQjmtgnnsXTXMb49dDygfWs0PzTcPRFfA0YB27FbAn/zou1koMDh82HjmMsySqkmoBJwfmW+CtislKp37kBEbhWRHBHJKSkp8UK0M/HGgmiursF24oTfRjA1NdtYv7/cdPdSC126dCErK4uNGze2OxYRPWki1ri4gLuZAG64sA+JkSH8dZFv+aU0Gk3buFMQ6Uqp64zlRecC471o21Vk1TmlaJtlRCQDu9vpNlcdKKVeUEqNUkqN6tq1fakoWvDGgmg6ah/i6q+FgrYdqaS6vsmU4a2tcdFFF2Gz2Vi1alW76ktICDGXXkrV0mUBnRMBEBkaxE8n9WNtftmplfY0Go35uHsinlppxnjD94bDQKrD5xSgsLUyIhKEfYZ2ufE5BVgI/EQp5fd0nt5YEN9PkvOPgliXb5/tPOY8/6URj4+PZ/jw4Xz77bccP94+V03clXNQDQ0Bn1kNMH9ML1Liw/njf3Zjs7U/lblGo2kddwpiqIicMLYqIKtl31gnoi02Av1FpI+IhADzgE+cynyCPQgNdgtluVJKiUgc9gWJ7ldKtd9R7gXeWBCNhf5VEGvzSxmUFE1ilF/DLkyYMAERabcVEZaeTlhGBhXvvOvTehPtITTIyj3TB7Kj8ASfbnN+79BoNGbgbj0Iq7EeRMuaEEEO+21mjzMsjruARcAu4F2l1A4ReVhELjOKvQQkikgecDfQMhT2LqAf8BsR2WJs3Xy4Trd4ZUEcLQKrlSAf3VquqG9qJufAcb+6l1qIiYlh9OjRbN26lfbGcOLmXUP93r3Ubt5isnTuuWxoTwb3iOGvi/fQ0BSY9bI1mh8S3uW39hKl1BdKqQFKqb5KqT8Yx36rlPrE2K9TSl2tlOqnlMpWSu0zjj+ilIpUSg1z2Ir9KatXMYjCIoK6dUOCvFvAxxM2H6qgvslm+vDW1hg3bhzBwcEsXbq0XfVjZ83CEhVFxTuBz7xisQj/b+ZACspreXO9u0F1Go3GW/yqIM4mvI1B+DP+YBHI9nJ50fYSGRnJ+PHj2bt3L/v3ez8z2hIZSexll3HiP1/S1M5Yhi9cNKArF5yXyNPL8zhR1+i+gkaj8RitIAy8SbXRePSo/xTEvjIyesYSGx7sl/Zdcf755xMbG8uSJUvaFUuIu+Yae7D6o4/9IF3biAgPzBpMeU0Dz6zQS5NqNGaiFYSBpy4mZbPRVFREsB/WgahtaGbLoYqAuZdaCA4OZvLkyRQVFbF9+3av64cNHED4iBFUvP12wIPVAJkpsVw5IpmXV+/X6cA1GhPRCsLAUxdTc1kZqrHRL3MgNh08TkNz4OIPjmRmZtKjRw+WLVtGY6P3rpr4edfQcPAgJ7/5xg/SuefeGQOxWODPevKcRmMaWkEYeJqsr7HQPqQy2A+J+tbtK8VqEUanBSb+4IiIMH36dE6cONGuFBzRM2ZgjY+n/P/e8IN07ukRG86t48/j062FbDqoU3BoNGagFYSBpzGI7xcK8oOCyC8jKyWWqFDzR0d5QlpaGhkZGaxZs8bryXOW0FDi5l1D9fLlNBzsmBFFt13Ul67RoTzy+c4OcXVpNOcaWkEYeGxBHPGPBVFT38S2w5VccF7g3UuOTJ8+HRFh0SLvV4yL//GPISiI8tf/zw+SuScyNIh7Zwxk86EKFm4+0iEyaDTnElpBGHhjQViiorBGR5va/8YD5TTZVIfEHxyJiYlhwoQJ7Nmzh9xc7xbmCe7WjdhZF1Px4Yc0n+iY9RrmjkhhaGocj/1nN9X13maH0Wg0jmgFgd16ADyOQfjLvRRsFUb1Dnz8wZkLLriAxMREvvzyS5qavHvIJlx/PerkSSree99P0rWNxSI8dFkGJVX1/F2vX63R+IRWEHyvIDy1IPwxB2LdvjKGpcYRHuL5okX+wmq1cvHFF1NeXs7XX3/tVd2w9HQiRo+m/I3/Q3mpXMxiWGocV49M4eU1+8nX61drNO1GKwjs7iXwwoJINteCqDzZyHdHKgOSf8lT+vbtS1ZWFqtXr6a42LssJwkLrqepsIiqJUv8JJ17/nfmIMKCrDz0qQ5YazTtRSsIPLcgmqtrsFVWmj4H4pv9ZdgUXNiv8ygIgBkzZhAWFsYnn3xy6h55QtTEiYT07k3piy922MO5a3Qo/zNtAKv2lvDF9qMdIoNGc7ajFQSeWxBNRf4ZwbQuv4ywYAvDUuNMbddXIiIimDFjBkeOHGHjxo0e1xOrlcRbb6F+5y5qvHRRmclPLuhNeo8YHvp0h87TpNG0A60g+F5BuLMgvl8oyFwFsSavlNFpCYQEdb6vIzMzk379+rFs2TIqKio8rhc7ezZBPXpQ+tzzHWZFBFktPHZlJiXV9fxNz7DWaLym8z2ROgBPRzGdmkVtYgyiuKqO3OLqTudeakFEuOSSSxARPv74Y48f9hISQuLNN1H77bec9ML6MJuhqXH8ZExv/v3NQbYWeK7gNBqNVhCA5y6mxsIiCAoiqIt5D/OW5UXHdvD8h7aIi4tj5syZHDhwgG+8yLUUd9VVWLt0oey55/0onXt+OWMgXaNCuf/D7TQ164WFNBpP0QoCz4PUjYWFBCclIR6uG+EJa/PKiAkLIqNnrGlt+oNhw4YxcOBAli1b5vGoJktYGIk3LKBm7Vpqt23zs4StExMWzEOXZbCz6ATPr9rXYXJoNGcbWkHgjQVRaPociLX7ShlzXiJWi5jartmICLNnzyY0NJSFCxeeumfuiLtmHtbYWEr+8Q8/S9g2F2f2YFZmEk8tzSX3WFWHyqLRnC1oBYEXFkRRkakjmArKT1JQXttp4w/OREZGMnv2bI4ePcqyZcs8qmONiiTx1luoWfV1h8YiAB66bAiRoVbufX8bzTY9N0KjcYdWEHhmQajGRpqOHSOop3kWxNr8UqBzxx+cGTRoEKNGjWLdunXs3bvXozrx8+cT1K0bxY8/0aGT1rpGh/LgZRlsKajg5dXeL6+q0fzQ0AoCzyyIpuJisNlMtSDW5JXRNTqUft2iTGszEMyYMYOkpCQ++ugjKisr3Za3hIXR5c47qd28meqVK/0vYBtcNrQnUwd356+L95BXrF1NGk1baAWBZxbEqSGuJs2BsNkUa/NLGds3EZHOHX9wJigoiLlz59Lc3MwHH3zgUTwi7qorCe7di5Inn0J5MSvbbESER68cQmRoEL94ZwsNTXpUk0bTGlpB4JkFYfZCQTuLTlBa3cCE/l1NaS/QJCYmMnv2bAoKCli6dKnb8hIcTNef/Zz6PXs48fnnAZCwdbpFh/HYlZl8d+QETy71zE2m0fwQ0QoCby2IJFP6/DrXHn8Y3//sCFC7YsiQIWRnZ/PNN9+wdetWt+VjZl1MaPpgih9/AlttbQAkbJ0ZGUnMG53Ks1/ls35fWYfKotF0VrSCwEMLorAIa0IClvBwU/r8OreEQUnRdIsJM6W9jmL69OmkpaXx6aefUmgo0dYQi4WkBx6gqaiIsn+9FCAJW+c3l6bTOyGCu9/dSmWtztWk0TijFQSeWxBmzYE42dBEzoHjTBhwdrqXHLFarVx99dVERUXxzjvvUF3d9voLEaNGETNrFmX/+heNRzp2WdDI0CCenDecYyfq+N/3t+q04BqNE1pB4FmyPjPnQKzfV05Ds+2sdi85EhERwbx586itreWtt96ioaGhzfLd7r0HRDj2578ESMLWGZYax30XD2LRjmO8vOZAR4uj0XQqtILAfbI+pZSx1Kg5FsSq3BJCgyyMTuv45UXNIikpiSuvvJLCwkI+/PDDNtePCO7Rg8Rbb6Fq0SJqvlkfQCldc9O4PkxP785jX+xi86HjHS2ORtNp8KuCEJGZIrJHRPJE5D4X50NF5B3j/HoRSTOOJ4rIChGpFhG/52hw52JqrqhA1daaZkGs2lvC+eclEhbc8cuLmsmgQYO4+OKL2bNnD19++WWbLpvEG28kODmZo7//PTY3Foe/ERH+MncoSbFh3PXmZo7XdKw8Gk1nwW8KQkSswD+Bi4F04Mciku5U7CbguFKqH/AE8CfjeB3wG+Aef8nniLsgdeOhQwAEp/byua8jFbXkl9Qw4RxxLzmTnZ3NmDFj2LhxI2vWrGm1nCUsjKQHH6QhP5+y554LoISuiY0I5pn5Iyipruenb36rs75qNPjXgsgG8pRS+5RSDcDbwOVOZS4HJD26eAAAHX9JREFUXjP23wemiIgopWqUUquxKwq/486CaDAUREhv3xXE6twSgHMiQN0a06dPZ8iQISxbtoycnJxWy0WNH0fs5ZdR+sKL1O3p+PkIWSlxPDYnk7X5Zfzhi10dLY5G0+H4U0EkAwUOnw8bx1yWUUo1AZVAwBMTubMgGg4eAhGCU1J87mvV3lK6x4TS/yxLr+ENIsIVV1zBgAED+Pzzz9nWRqrvbvfdhzUmhqLf/AblYYZYf3LVyBRuvLAPr6w5wHs5Be4raDTnMP5UEK7yRzg7pT0p03oHIreKSI6I5JSUlHglnCPuLYiDBCUlYQkNbXcfAI3NNlblljChf9ezLr2Gt1itVubOnUtaWhofffQRu3a5fiMPio+n+68eoG7bNsr//XqApXTNA7MGcWG/RH618DtyDpR3tDgaTYfhTwVxGEh1+JwCOM+kOlVGRIKAWMDj/0il1AtKqVFKqVFdu7bfZdNiQbT20G48eIiQXr67lzbuL6eqromp6d19butsIDg4mHnz5pGcnMx7773Hjh07XJaLmTWLqClTKHniCer2dPza0UFWC/+8dgQp8eHc/O8c8kvantuh0Zyr+FNBbAT6i0gfEQkB5gGfOJX5BLje2J8LLFcdMFupubkZi8XSqoJoKCgwRUEs3VVMSJDlnJn/4AmhoaFcd911pKam8sEHH7h0N4kIPX7/MJbYGI788pfY6gISemqTuIgQXr0hG6sIC17ZQElVfUeLpNEEHL8pCCOmcBewCNgFvKuU2iEiD4vIZUaxl4BEEckD7gZODYUVkQPA48ACETnsYgSUadhsttaHuFZV0Vxe7nOAWinFkl1HubBvIhEhQT61dbYRGhrK/Pnz6d27NwsXLuTbb789o0xQQgI9H/sjDXn5FP/5zx0g5Zn0SozgpQWjKamq5+bXNnKyoamjRdJoAopf50Eopb5QSg1QSvVVSv3BOPZbpdQnxn6dUupqpVQ/pVS2UmqfQ900pVSCUipKKZWilNrpLzlbLAhXtIxgCvbRgsgtrqagvPYH415yJiQkhGuvvZZ+/frx6aefsmrVqjPmSUSNu5CEG27g+JtvUbV8eQdJejrDUuP4+49HsP1IJbf+exN1jR0fSNdoAoWeSY1dQbRmQTSeGuLa26c+luw8BsCUQT9MBQHfxySysrJYsWIFX3zxxRkzrrv+zy8ITR9M4X33n1LOHc209O78ee5QVueVctebm2nUcyQ0PxC0gsDuYmpziCsQkprq8rynLNt1jKyUWJJiz+7srb5itVq54oorGDt2LDk5Obz77run5W6yhISQ8tRTIMLhn96FraamA6X9nrkjU/j95Rks3XWMX767Va9prflBoBUEbccgGg4dwtq1C5aIiHa3X1JVz+aCih+09eCIiDBt2jQuvvhi9u7dy8svv0xFRcWp8yGpqSQ//jfq8/MpfOBXnSbL6n9dkMZ9Fw/ik62F/PLdLXq2teacRysI2nYxNRw6SEgv39xLK3YXoxRMTe/mUzvnGtnZ2cyfP5+KigpefPFFDh48eOpc1IUX0u2Xd1O1aBFlL7zYgVKezu0X9eXeGQP5aEshP397s16yVHNOoxUEbbuYzJgDsWTXMXrGhpHeI8ands5F+vbty80330x4eDivvfYa/7+9e4+PqjwXPf575pKZJJP7jZCEJFzDVQQEFOvGCxsUEa2ytdsjglq13lr37jlHd8+x2n5s7a67PXZXW0FRtFapWhXv0ooiAhUk3BVEQsIlMVwScr/MzHP+WCshwAQhZDKQeb+fz3xmzTtrvetZXNYz73rXet8VK1a0txhSb7qJxOnT2ffb33LozTcjHOlhd144kP8zfSjvbKzgjhdMx7XRe5kEQectiGBDA/59+07pFtfaplaWbdvHlGFZvf7p6a5KT0/nlltuoaioiCVLlrBo0SIaGxut5yN+8TBx48ez9z9+Qt1xBv/rabd8pz8/v3IEf/uiktkLPuNQg5mRzuh9TIKg8xZEy67dAKfUgnh/8zc0+4NcMfroYaiMjrxeL7NmzWLq1Kl89dVXPPnkk+zcuROHx0Pu47/HU1jInrvvobGTp7Ej4YaJ+Tx23WiKy6qY9eQK9lZHdp5tw+huJkHQeQuipcy6Ju4+hT6IN9btoV9qHGP6JXe5jmghIkycOJG5c+fidDpZuHAhS5YsQWNjyZs/D0dyErtu+f5pMfJrm5mjc1h403jKq5u46olP2bTnUKRDMoxuYxIEnbcg2p+B6Ne1W1wra5v4dPt+Zo7uay4vnYTc3Fxuu+02xowZw4oVK5g/fz6Vfj/5CxYgMTGUzZlzWozZ1Oa8Aem8/INzcYhwzR9XsHj90UOOGcaZySQIjtOCKC3DmZKCM7FrnctvrS8nqDBzdPfMRBdNYmJimDFjBt/73vdobGzk6aefZunWrfR9+inE46Hsxjk0ffllpMNsV9QnkcV3nc/InCTuebGYX777hXlWwjjjmQTBcfogyspwd7H1ANblpeF9ExmYmXAq4UW1wYMHc8cddzBmzBhWrVrF/Lfewv/AA4jXS+nsG2lYvTrSIbbLSPDwwi0TuX5CP578eAc3PP0PKmsiP/CgYXSVSRAcvw+iq89AlOyvZ/3uQ1xpOqdPmdfr5fLLL2fOnDl4PB7+8skyiq+5GklOpuymm6l5991Ih9guxuXg4atG8p/XjKK4rJpLH/uEj7ZWRjosw+gSkyAIPVhfsKUFf3lFl+9ger14DyIw4yxzeam75Ofnc9tttzFt2jR21NXx6jnjaMjNZc+9/8aBp58+bZ64BviXcXm8efck0n0e5jyzmofe3Exji3lewjizRNe4050INdRGy44doEpMQUEX6lNeX7eHiYVpUT/2UndzOBxMmDCBESNG8PHHH/Ou2834YBB+/Sh16zeQ98tf4IiPj3SYAAzMTOCNuybxy3e+4JlPd/LR1n38+ppRjCtIjXRohnFCTAuC0JeYGu2JbWJHjTzp+j7aVknpgQauG39qA/wZnYuPj+eyyy7j9rvvpvqmuWw4axT1S5aw6fIZ1HQyvWkkeN1OHpo5gj9/fwKtgSCznlzJT9/YxKFG82CdcfozCYLQndSNGzbgTErq0jwQT31SQnaSl8tGZndXiEYn0tLSmHXttUz+3e8omT2bwMGDlM76F5b/9MEjBgCMtPMGpPP+jy5g9sR8nl9VysX/9TGvFe8+rS6LGcbRTIIgdB9E04aNeEeNOunnF7bsrWHF1weYfW4Bbqf54+0pWVlZTL//Pvr86Xma8/NJW7SItd+9mjcWLmT37tPjRBzvcfHQzBEsvut8clJiuXfReq7540o+Lz3hadgNo0eZMxjH9kEE6+tp3r6d2JEnf3lpwaclxLqd/Ov4U5/D2jh52SNHMvbNxSTeey9Z+/ZR8F+/Yen99zN/3jyKi4uPmHsiUkbkJPHaD87jke+OZNfBBq7+w0puf/5ztlfWRTo0wziC6aTm2BZE4+bNEAziPcn+h8raJhav28u15+SRFOfu7jCNEyQOBzm33Ur6lEvY+8ADjFu9hpqyMj7ZvIX3svswfPhwzj77bHJzcyP2hLvDIVw3vh9XjO7L/GUlPLnsa97fUsGMUX25+6KBDMoyz84YkWcSBMe2IJo2bgQgdtSok6rnT6vKaAkEmTupoDvDM7rI078/Bc8/T81bb+H81a+Y8sEH1IwezcqaGoqLi0lJSWHEiBGMGDGCzMzIzNURF+Pih5cM4vqJ/XjqkxKeW7mTNzfsZeqwPtzynULG5qeYYVqMiDEJghAtiPUbcOfm4ko98dsRDzW28qdVpVxclEn/DF84wjS6QERImjED3+TJHJj/FPLcc0zdtInARReyOT2d5cuX88knn5Cens7QoUMZOnQoffr06fGTcrrPw32XFnHrBf1ZsLyE51eV8t7mCs7KS2bueQVMG9EHrzv0pFaGES5RnyBUFVU9ogXRuHEjcWeffVL1/HbJNqoaWrh3yuDuDtHoBs6EBDL/7V5Srr+e/U88QfWrrzLqw6VMnDqVfZPOY3NtbXuySEhIYNCgQQwaNIjCwkI8Hk+PxZkaH8OPpw7hjgsH8OraPSxYXsKPFq0j+U03V4/J5dpz8hhsLj8ZPSTqE0QwaE0Z2daCaK2sxF9ejvfG2Sdcx5a9NTy3cifXT+jHiJykcIRpdBN3VibZDz1I+u23cfDZZ6l6+RXi336byRMmEHfVlezOyeGrkhI2bdrE2rVrcTgc5Obm0r9/fwoKCsjJycHlCv9/m7gYFzdMzOf68f1YueMAf/6sjOdW7uTp5SUMy07kqrNzuPysbLKTYsMeixG9oj5BBALW8AdtLYjD/Q9nndD2qspPF28iKdbNj/95SHiCNLqdOzubrPvvJ+3226l++RWqFy1i/333E5+WxgXTL2P69Ol8k5DAjh07KCkp4aOPPgLA5XKRm5tLXl4eeXl55ObmEhsbvpO0wyFMGpjOpIHp7K9r5q31e3lt3V4efucLHn7nC0bnJTNtRB+mDMuif3q86a8wupVJEHaCaGtBNG7YCC4X3mFDT2j714r3sHpnFY98dyTJcTFhi9MID1dKCum3fp+0m2+ifvlyql95heoXX6LqueeJKShg1JQpTLrkYnTgQMp27aK0tJTS0lKWL19+eO7s1FRycnLo27cv2dnZZGVl4fV2/xAr6T4PcyYVMmdSISX763lnYznvb67gkXe/5JF3v6RfahyTh2RwwaAMxvdPJdFr7qQzTo2cDg8QdYdx48bpmjVrTnq7+vp6Hn30US699FLGjx9P6dy5BA/VUPjXV791210HG7jy8U/JTY3jtR+ch8Nhfr31BoFDh6j54ANq3n7HGk48EMCVkUH8pEnW67xzCSYksHfvXnbt2sXevXvZs2cPtbW17XUkJyeTlZVFZmYmmZmZZGRkkJqaitvd/SftPdWNfPhlJR9vreTT7QdobA3gEBiZk8SE/mmMzU9hTL8UMhJ6ri/FOHOIyOeqOi7Ud6YF0eESkwaDNG3cROLl0791u/11zcxe8Bn+oPLoNaNMcuhFnElJpMyaRcqsWQSqq6lbtozaD5dSt3Qph15/HYCYAQOIGzOGUePGMv6c8cTMyqe+oYHy8nIqKiqoqKigsrKSbdu2HfEUd0pKCmlpaaSmpra/UlJSSE5O7nLfRk5yLDdMzOeGifk0+wOsLa1m5Y4DrPx6P89+upN5y3YAkJcay6icZEbmJjGibxJF2Qmk+0zSMDoX9QmiYyd1/acrCNbVEXvW6ONuU9fsZ+4zqyk/1MgLt0wwDzX1Ys7kZJKuuIKkK65AAwGatmyhfsVKGtZ+Ts1771H98ssAOHw+vMOHk1g0hIzBgxlXVITnsssIejzs37//iNfBgwcpKys75qlun89HUlJS+ysxMZHExEQSEhJISEjA5/N9axLxuJycOyCNcwekwZTBNPsDbNpTw9rSKop3VbFhTzVvbyxvXz/d52FIHx8DM3wMyPQxIMNHfloc2UmxOM2PnqgX9QmivQXR0ED5T35CTP/+JE6b2un639Q08cOXitlSXsP82WMZm2+Gbo4W4nQSO3Jk+xAsGgzS/NV2mjZtonHTRpo2baZq0V/QpsOzyLmysogpLCStXz/65OUSk5eHe/BgXFl9aIr1Un3oEFVVVVRVVVFdXU1NTQ0VFRVs27YNv99/TAxerxefz4fP5yM+Pp64uLj297ZXbGxs+yvG7WZsfgpj81Pa66iqb2FLeQ1fVtSytaKGrRW1vLp2D3XNh/cX43SQmxpLbkocOcmx5KbEkp3kpU+ilz5JXjITvfg8UX/66PXC+jcsItOAxwAn8JSqPnLU9x7gOWAscAC4VlV32t/dD9wMBIB7VPX9cMQYDAZBFdf8+firqij84x9whLgrxR8IsnBlKb9dso2WQJBHZ43ioqKscIRknCHE4cA7ZDDeIYNJvvq7AGggQOuuXTRt20bLjhJaSnbQXLKT2iVLCFRVHVmB2407I4PEjAxSMzNwpqfjSk3DmZqKs/8A/HGxNLpcNDgc1KlS29pKXUMD9fX11NXVUV5eTn19Pc3NzZ3G6HA48Hq9R7w8Hg8ej4ccj4fC1BhmZntwu1OoC7iobILKBqWiLkB5bSvf1DSwcXc1VQ3HDk8eF+MkI8FDWnwMqfEe0n0xJMfFkBLnJiUuhsRYN0mxbhJjXSR63SR4Xfg8LlxmEMszRtgShIg4gceBKcBuYLWILFbVLR1WuxmoUtWBInId8CvgWhEZBlwHDAf6An8TkcGq2u1TcgUCAQpKSpB/fEbmj/8d77Bh7d/VNftZs/MgK3cc4G9bvuHrffX80+AMfjZzOPlpp8ekNMbpRZxOYgoKQk40Fairo3X3blr3ltNaUY6/vBx/ZSX+fftoLikhsOZzAtXVEOLGkVggVoQ+Ph8Onw+nLx5HXDwSF4vExhF0uwm63fhdTvwOB36Hg1agBaGluZmW6mpagkpLMEizBqn2B2gJBmgOBAiKg6BDUIeDoFjvKSIkiYMhIqhDaHU4qSeGejw0OmNpcnhoEC+N9TEcqHezO+ikIeCgMeAgwPEvTbkd4HVBrEvwOAWvy4HXLXicDrwuwet24HE5iHE58NrvMU7r3eNy4LaXY1xO3E4hxunE5RTcTgfutjKXE5fTgcvhwOkUXA6H/Vlwdih3Ohw4HVa5w17HIVZiFZH224aPfo8W4WxBjAe2q+oOABF5CZgJdEwQM4EH7eVXgN+L9TcwE3hJVZuBEhHZbte3sruD/HrDdkZ+vo6tfQZyV1V/go9+RF2zn5rGVpr9Vv+E2ymMzkvmf04dwtThPT8Mg9E7OH0+nEVFeIuKOl1H/X4CVVUEqqsJVFfjr64mWFNLoLaGYE0Ngbo6gnX1BGtrCTY0EGxoIHCwimBDA9rUhDY2Ik1NuP1+3EBcDxyX2v8fVAREUBEUQQXr3V4OIih2mXRYPmZd6zs6fH94X4fLCVEO4AdaO37XMdhO/u92di+nIkelu2+56/Mkzg3ybXWdhPL+Q5m96Mluq69NOBNEDrCrw+fdwITO1lFVv4gcAtLs8lVHbZtz9A5E5FbgVoB+XZw7Oj7WQ2lmPktnfJ+irGQcIvg8rvZm8ajcJMbmpxAXY663GuEnLheujAxcGRmnVI8GAmhzM8HmZrSlBW17b221Xi0tqN+PtvqtzwE/+P1WWSAAgYD1XTAAgaD1fSAIGkQDQQgG0GDQKkNDL6van9VaVrW2V4Wg2pd3g6i9rBo8XG5vZ5UrgaASCAQJqFV/wP4uqLS/W9Vq+/A5QQUNKorau1crKam2N9LayrH2eLgcIKiHl9ve7Tra/5yP+nD0Kf+Iz6rHloVa/wTyxtGruPuE53J3OM96oVLp0cfV2Tonsi2qOg+YB9ZzECcbIMDo75zD6A8XM7MrGxvGaUqcTiQuDkdcT7QhjN4qnL1Fu4GOkzLnAns7W0dEXEAScPAEtzUMwzDCKJwJYjUwSEQKRSQGq9N58VHrLAZutJevAT5Uq723GLhORDwiUggMAj4LY6yGYRjGUcJ2icnuU7gLeB/rNtcFqrpZRH4GrFHVxcDTwPN2J/RBrCSCvd5fsDq0/cCd4biDyTAMw+hc1I/FZBiGEc2ONxaTeWLFMAzDCMkkCMMwDCMkkyAMwzCMkEyCMAzDMELqNZ3UIrIPKD2FKtKB/d0Uzpkg2o4XzDFHC3PMJydfVUM+ut9rEsSpEpE1nfXk90bRdrxgjjlamGPuPuYSk2EYhhGSSRCGYRhGSCZBHDYv0gH0sGg7XjDHHC3MMXcT0wdhGIZhhGRaEIZhGEZIJkEYhmEYIUV9ghCRaSKyVUS2i8h9kY4n3EQkT0SWisgXIrJZRH4Y6Zh6iog4RaRYRN6KdCw9QUSSReQVEfnS/vs+N9IxhZuI3Gv/u94kIi+KiDfSMXU3EVkgIpUisqlDWaqILBGRr+z3lO7YV1QnCBFxAo8DlwLDgO+JyLDIRhV2fuDfVXUoMBG4MwqOuc0PgS8iHUQPegx4T1WLgLPo5ccuIjnAPcA4VR2BNc3AdZGNKiyeBaYdVXYf8HdVHQT83f58yqI6QQDjge2qukNVW4CXoHfPPqqq5aq61l6uxTppHDPfd28jIrnAdOCpSMfSE0QkEbgAa84VVLVFVasjG1WPcAGx9gyVcfTCmShVdRnW/DkdzQQW2ssLgSu7Y1/RniBygF0dPu8mCk6WbUSkADgb+EdkI+kR/w/4X0Aw0oH0kP7APuAZ+7LaUyISH+mgwklV9wCPAmVAOXBIVT+IbFQ9JktVy8H6EQhkdkel0Z4gJERZVNz3KyI+4FXgR6paE+l4wklELgcqVfXzSMfSg1zAGOAPqno2UE83XXY4XdnX3WcChUBfIF5E/kdkozqzRXuC2A3kdficSy9skh5NRNxYyeEFVf1rpOPpAZOAK0RkJ9ZlxItE5E+RDSnsdgO7VbWtdfgKVsLozS4BSlR1n6q2An8FzotwTD3lGxHJBrDfK7uj0mhPEKuBQSJSKCIxWB1aiyMcU1iJiGBdl/5CVX8T6Xh6gqrer6q5qlqA9Xf8oar26l+WqloB7BKRIXbRxVhzvPdmZcBEEYmz/51fTC/vmO9gMXCjvXwj8EZ3VOrqjkrOVKrqF5G7gPex7nhYoKqbIxxWuE0CbgA2isg6u+w/VPWdCMZkhMfdwAv2j58dwNwIxxNWqvoPEXkFWIt1t14xvXDYDRF5EZgMpIvIbuCnwCPAX0TkZqxEOatb9mWG2jAMwzBCifZLTIZhGEYnTIIwDMMwQjIJwjAMwwjJJAjDMAwjJJMgDMMwjJBMgjDOaCKSJSJ/FpEdIvK5iKwUkau6WFdBxxEyw0lExonI73piXx32OVpELuvJfRpnNpMgjDOW/TDU68AyVe2vqmOxHoTLjWxkxyciLlVdo6r3hKPu43w9GjAJwjhhJkEYZ7KLgBZV/WNbgaqWqup/A4iIV0SeEZGN9oB1F9rlBSLyiYistV/HHY5BRLJFZJmIrLPnGfiOXT7N3n69iPzdLksVkddFZIOIrBKRUXb5gyIyT0Q+AJ4Tkclt81LY3y0QkY/sltA9Hfb9f+35HJbY8xv8OER8z4rIb0RkKfArERkvIivsY14hIkPsh+V+BlxrH8e1IhJv73e1vW6vHsnYOHlR/SS1ccYbjvXUbGfuBFDVkSJSBHwgIoOxxqmZoqpNIjIIeBEYd5x6/hV4X1UftucQiRORDGA+cIGqlohIqr3uQ0Cxql4pIhcBz2H9cgcYC5yvqo0iMvmofRQBFwIJwFYR+QPWHA5XY42467KPtbMBBwcDl6hqoG2ob3ukgEuAX6jq1SLyANZcCXcBiMgvsIYduUlEkoHPRORvqlp/nD8LI4qYBGH0GiLyOHA+VqviHHv5vwFU9UsRKcU6kZYCvxeR0UDALjue1cACe5DD11V1nX2CX6aqJXb9bePzn491UkdVPxSRNBFJsr9brKqNnezjbVVtBppFpBLIsut6o20bEXnzODG+rKoBezkJWGgnPwXcnWzzz1iDGLa1SrxAP6Jn/CLjW5gEYZzJNmOfjAFU9U4RSQfW2EWhhnMHuBf4BusXugNoOt5OVHWZiFyANeHQ8yLya6Ca0EPDH28I+eP9Mm/usBzA+r/ZWfyhdKz758BSVb1KrDk/PupkGwGuVtWtJ7EfI4qYPgjjTPYh4BWRH3Qoi+uwvAy4HsC+tNQP2Ir1C7tcVYNYAxc6j7cTEcnHmk9iPtZIuGOAlcA/iUihvU7bJaaO+5wM7D+F+TaWAzPsvhQfVoI6EUnAHnt5TofyWqxLWG3eB+62O/sRkbO7GKfRS5kEYZyx1Bpp8kqsE3WJiHyGNd3i/7ZXeQJwishGYBEwx76M8wRwo4iswrq89G3X3CcD60SkGKvF8piq7gNuBf4qIuvt+gEeBMaJyAasETZvPLa6Ez6+1VjDOK/HmttgDXDoBDb9T+CXIvIpRya/pcCwtk5qrJaGG9hg3977867GavROZjRXwziNiYhPVetEJA6rdXJr25zihhFupg/CME5v80RkGFYH8kKTHIyeZFoQhmEYRkimD8IwDMMIySQIwzAMIySTIAzDMIyQTIIwDMMwQjIJwjAMwwjp/wOMv+NdjdIJzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "prior.plot(label='prior', color='gray')\n",
    "france.plot(label='France', color='C0')\n",
    "croatia.plot(label='Croatia', color='C3')\n",
    "\n",
    "plt.xlabel('Goal scoring rate')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Posterior distributions for France and Croatia')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6999765866755225"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmf_mean(croatia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook, I introduced the World Cup problem, which involves estimating the goal-scoring rate of a soccer team, based on the outcome of a single game.\n",
    "\n",
    "To solve this problem, I introduced the gamma distribution, which we used to generate a prior, and the Poisson distribution, which we used to compute a likelihood.\n",
    "\n",
    "We used the prior and the likelihood to due a Bayesian update, which computes posterior distributions that estimate the goal-scoring rate for each team.\n",
    "\n",
    "[In the next notebook](https://colab.research.google.com/github/AllenDowney/BiteSizeBayes/blob/master/09_predict.ipynb) we will use those distributions to answer the questions we started with:\n",
    "\n",
    "1. What is the probability that the winning team is actually the better team?\n",
    "\n",
    "2. What is their probability of winning a rematch?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}