{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Think Bayes\n", "\n", "This notebook presents example code and exercise solutions for Think Bayes.\n", "\n", "Copyright 2018 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Configure Jupyter so figures appear in the notebook\n", "%matplotlib inline\n", "\n", "# Configure Jupyter to display the assigned value after an assignment\n", "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n", "\n", "# import classes from thinkbayes2\n", "from thinkbayes2 import Pmf, Suite\n", "\n", "import thinkbayes2\n", "import thinkplot\n", "\n", "import numpy as np\n", "from scipy.special import gamma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### World Cup problem, part two\n", "\n", "> In the final match of the 2014 FIFA World Cup, Germany defeated Argentina 1-0. How much evidence does this victory provide that Germany had the better team? What is the probability that Germany would win a rematch?\n", "\n", "Scoring in games like soccer and hockey can be modeled by a Poisson process, which assumes that each team, against a given opponent, will score goals at some goal-scoring rate, $\\lambda$, and that this rate does not vary; in other words, the probability of scoring a goal is about the same at any point during the game.\n", "\n", "Based on this modeling decision, we can answer the questions by\n", "\n", "1. Defining a prior distribution for each team's goal-scoring rate against the other,\n", "2. Updating the prior based on the outcome of the game,\n", "3. Using the posterior distributions to compute the probability that Germany's goal-scoring rate is higher.\n", "4. Generating a predictive distribution for the number of goals each team would score in a rematch.\n", "\n", "I'll start with Step 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 2: Updating\n", "\n", "If goal-scoring is a Poisson process, the distribution of goals per game is Poisson with parameter $\\lambda$. To compute the distribution of $\\lambda$ we can define a new class that inherits from `thinkbayes2.Suite` and provides an appropriate `Likelihood` function:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "class Soccer2(thinkbayes2.Suite):\n", " \"\"\"Represents hypotheses about goal-scoring rates.\"\"\"\n", "\n", " def Likelihood(self, data, hypo):\n", " \"\"\"Computes the likelihood of the data under the hypothesis.\n", "\n", " hypo: goal rate in goals per game\n", " data: goals scored in a game\n", " \"\"\"\n", " # FILL THIS IN!\n", " return 1" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "from scipy.stats import poisson\n", "\n", "class Soccer2(thinkbayes2.Suite):\n", " \"\"\"Represents hypotheses about goal-scoring rates.\"\"\"\n", "\n", " def Likelihood(self, data, hypo):\n", " \"\"\"Computes the likelihood of the data under the hypothesis.\n", "\n", " hypo: goal rate in goals per game\n", " data: goals scored in a game\n", " \"\"\"\n", " return poisson.pmf(data, hypo) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Likelihood` computes the likelihood of `data` given `hypo`, where `data` is an observed number of goals, and `hypo` is a hypothetical goal-scoring rate in goals per game. We can compute the likelihood of the data by evaluating the Poisson probability mass function (PMF).\n", "\n", "Now we can get back to Step 1.\n", "\n", "### Step 1: Constructing the prior\n", "\n", "Before the game starts, what should we believe about each team's goal scoring rate against each other? We could use previous tournament results to construct the priors, but to keep things simple, I'll just use the average goal-scoring rate from all matches in the tournament, which was 2.67 goals per game (total for both teams).\n", "\n", "To construct the prior, I use a gamma distribution with a mean of 1.34 goals per game." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.310359949002256" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcVOWd7/HPr3rvptlBdhrZBEQFGhQR3Pc9ccFJ5ppcJyZ34mSbycRkJpNJbiYzyeSqycSZGyc68ZrENdEwERUVFAXZFQEB2WWnWbtpequq5/5xTldXFb0BXX2qqr/v16tfnHPq1KkfpTy/fp7znN9jzjlERETSTSjoAERERJqjBCUiImlJCUpERNKSEpSIiKQlJSgREUlLSlAiIpKWlKBERCQtKUGJiEhaUoISEZG0lBt0AB2lb9++rqysLOgwRESkDStXrjzonOvX1nlZk6DKyspYsWJF0GGIiEgbzGxHe87TEJ+IiKQlJSgREUlLSlAiIpKWlKBERCQtKUGJiEhaUoISEZG0pAR1BpxzRKPRoMMQEclKWfMcVGfac+Aory9ez1vLPyYajfKFO2dyyeRRQYclIpJVlKBOwd6KYzz23Dt8+PGuhOMPP/kGO/ceZvYNUzGzgKITEckuGuI7Bf/5/MnJqdEL81bxr0/Mo7auoZOjEhHJTkpQ7XTgcBWrNzYlp6nnlvG3913LBecMjR1b+uE2fvL4azjngghRRCSrKEG105tLNsS2LzhnKA9+4TouPG8E37n/em669LzYa6s37mLB0o1BhCgiklWUoNohGo2yYGlTgrrionNi2zk5IT7/qYu5+bKmJPXrl97jaNWJTo1RRCTbKEG1w+qNuzl0tBqAbsUFTDu37KRz7rlxKv17lwJQXVPHE39Y3JkhiohkHSWodogf3rts6ljy8nJOOqcgP48v3j0rtr9o1WZWrmtXRXkREWmGElQbKo/XsGzNtth+/PBesgvOGcqs8tGx/ceef4e6es3qExE5HUpQbXh7+SYiEa9axKhh/Rk+qHer53/utovpVlwAwMEjx3l98fqUxygiko2UoNowP25yxJWt9J4a9Sgt4u7ry2P7L735AfUN4ZTEJiKSzZSgWlFVXcsnew8DkJub0+5yRldNH0ev7sUAHKk8kXAPS0RE2kcJqhX7D1bGtgf260FxUX673pefl8utV1wQ23/xjfdpaIh0eHwiItlMCaoV++IS1IA+3U/pvdfMGEf3bkUAHDpazYJlenhXRORUKEG1Yt+huATV99QSVEF+Hrdd2dSL+sPr7xMOqxclItJeSlCt2HfwWGx7QN8ep/z+a2eMp7SkEICKI1W8s3Jzh8UmIpLtlKBaEX8P6qxT7EEBFBbkcfPlTSWQ/vT2GhWSFRFpJyWoViTcgzqNBAVwzcXjycv1Kk9s332Qj7bs7ZDYRESynRJUC+obwhw+5tXfC5nRr1e307pOaUkhl00bE9uf+/aaDolPRCTbKUG1IL731LdXN3JzT66/1143zJoY21764TYqDledUWwiIl2BElQLEof3Tn2CRLxhA3szccxgABzw6rvrzuh6IiJdgRJUC+InSAzod3r3n+LF96JeX7xeRWRFRNqgBNWCM51inqx8wrCE9aIWrth0xtcUEclmSlAt2B/3kO5ZfUrP+HqhUIjrZ50b25+7cK2mnIuItEIJqgX7kurwdYQrLzonNuX8k72H2bTjQIdcV0QkGylBNSMajXIgbqbdWadYh68lJUUFXDKlqSL6a4s+6pDriohkIyWoZhw8Wh1bpLBHaRGFBXkddu1rLh4f2160ajPVNXUddm0RkWyiBNWMfRUdO0Ei3ujh/Rk+qA8ADeEIby37uEOvLyKSLVKaoMzsOjPbaGabzezBZl4vMLNn/deXmllZ0uvDzOy4mf1NKuNMFn//qSMmSMQzs4Re1OuLP9JkCRGRZqQsQZlZDvAocD0wHrjHzMYnnXYfcMQ5Nwp4GPhx0usPA6+kKsaWdPQU82SzykdTkO8NG+7cd4SN2/Z3+GeIiGS6VPagpgGbnXNbnXP1wDPArUnn3Ao86W+/AFxpZgZgZrcBW4FOL7uwvwOKxLamuCifSyaPjO3PW6zJEiIiyVKZoAYDO+P2d/nHmj3HORcGjgF9zKwE+Bbw/dY+wMzuN7MVZraioqKiwwLfm4Ip5snih/kWv7+F4yc0WUJEJF4qE5Q1cyz5ZktL53wfeNg5d7y1D3DOPeacK3fOlffr1+80wzzpmkn3oDq+BwUwclg/ygb3BbzJEu+sVGUJEZF4qUxQu4ChcftDgD0tnWNmuUAP4DBwIfATM9sOfA34jpk9kMJYYyqP18bq5BUW5NG9W2FKPsfMuHr6uNj+G+9tSMnniIhkqlQmqOXAaDMbYWb5wGxgTtI5c4B7/e07gPnOM9M5V+acKwMeAX7knPtFCmONSZ4g4d8SS4mZ5aMSFjPcurPjhilFRDJdyhKUf0/pAeA1YD3wnHNunZn9wMxu8U97HO+e02bgG8BJU9E72yF/kUKA/r1Pb5HC9iopKmD6BWfH9tWLEhFpkpvKizvn5gJzk479Q9x2LXBnG9f4x5QE14JwOBLbzs9P6dcDePX5Giubv7NyE5+7fTr5ean/XBGRdKdKEkkikaZ5HDmh1H89E0YNik1lP1Fbz3sfbE35Z4qIZAIlqCThSFMPqjMSlJlxxUXnxPbfXKJhPhERUII6SXwPKje3c76ey6eNjc23X7d5D3vjagGKiHRVSlBJItFobDs3p3O+nt49SpgyYXhsf8HSjZ3yuSIi6UwJKkk40pSgOmOIr9GVcc9EzV+6Ibbch4hIV6UElSQ+MeR0Ug8KYPK4ofQoLQLgSOUJPtiws413iIhkNyWoJJ09SaJRbm4Ol00dE9ufr8kSItLFKUEliUTjppl3Yg8KSJjNt2ztDo5V1XTq54uIpBMlqCTRgIb4AIac1YuxIwZ4cUSjsQd4RUS6IiWoJImTJFJXh68lV140Nrb95pL1Wm1XRLosJagkidPMczr98y++YGTCarubPznQ6TGIiKQDJagk4XD8EF/n96CKCvO5eFJTAVlVlhCRrkoJKknUBduDArjqoqZnot5ZuZnauoZA4hARCZISVJKEe1AB9KAAxo44i0H+UvO1dQ0qICsiXZISVJJIQJUk4qmArIiIEtRJ4ntQQQ3xAVw2bSwhfzXf9Vv3svvA0cBiEREJghJUkvjnoEIBDfEB9OpenFhAVr0oEelilKCSxFeSCLIHBXDl9KZhvreWf6wCsiLSpShBJUmsxRdcDwpg8rhh9OpeDHgFZFet/yTQeEREOpMSVJKEJd87udRRspycUEIB2Tff0zCfiHQdSlBJgqpm3pL42Xwr1+3gSOWJAKMREek8wbfAaSbxHlTwX8+g/j0ZP3IgAFHntNquiHQZwbfAaSaoBQtbc1XcarsqICsiXUV6tMBpJKgl31sz/YKzKS7MB2DfwUrWbtoTcEQiIqmXHi1wGkmoZp6bHl9Pfl4ul8ZNlnj9vfUBRiMi0jnSowVOI5Fw0ySJUJr0oACuinsmasnqrVRV1wYYjYhI6qVPC5wmoi69Jkk0Khvcl5FD+wHefbK3ln0ccEQiIqmVPi1wmgin4SSJRposISJdSXq1wGkgEvCS762ZOWVUwmq7G7ftDzgiEZHUUYJKki7VzJtTVJjPjEkjY/uaLCEi2UwJKkk6PgcV75oZTcN8i1Zt5viJugCjERFJnfRrgQOWbpUkko0a1p+ywX0BaAhHeHu5JkuISHZKvxY4YOlWiy+ZmXHNxU29qHmLPtJkCRHJSunXAgcsmjDEl16TJBrNnDI6Nlli1/4jrN+6L+CIREQ6nhJUnGg0SmNfxEivB3XjFRflM6t8VGx/3qKPAoxGRCQ10rMFDkg4Ybn39P5qrrl4fGz7PVWWEJEslNJW2MyuM7ONZrbZzB5s5vUCM3vWf32pmZX5x6eZ2Qf+z2ozuz2VcTaKptFy7205e2i/WGWJcDjCgmVahkNEskvKEpSZ5QCPAtcD44F7zGx80mn3AUecc6OAh4Ef+8fXAuXOuQuA64BfmlluqmJtFE7jh3Sbc82Mpq/ztXfXabKEiGSVVPagpgGbnXNbnXP1wDPArUnn3Ao86W+/AFxpZuacO+GcC/vHC4FOaXnT/RmoZJdMHpWwDMfqjbsCjkhEpOOkshUeDOyM29/lH2v2HD8hHQP6AJjZhWa2DlgDfCkuYcWY2f1mtsLMVlRUVJxxwPFTzNPxGahkhQV5XH7h2Nj+q++sCzAaEZGOlcpWuLkxsuSeUIvnOOeWOucmAFOBb5tZ4UknOveYc67cOVfer1+/Mw44/iHdTOhBAVx7yYTY9oq12zlwuCrAaEREOk4qW+FdwNC4/SFA8lKwsXP8e0w9gMPxJzjn1gPVwLkpi9QXSajDlxkJanD/npw3ZgjgZfbXNeVcRLJEKlvh5cBoMxthZvnAbGBO0jlzgHv97TuA+c45578nF8DMhgNjge0pjBVImmZu6T9JotF1M5t6UW8s2UBDQ6SVs0VEMkPKEpR/z+gB4DVgPfCcc26dmf3AzG7xT3sc6GNmm4FvAI1T0S8BVpvZB8CLwF865w6mKtZG0bjl3nNy03uaebzyCcPp07MEgMrjNby3ekvAEYmInLmUTt12zs0F5iYd+4e47Vrgzmbe9xTwVCpja044HD/NPDOG+MC7X3bNjAk8/fIyAOYuXMus8jEBRyUicmYypxXuBJFo5t2DanTV9HNiEzs27TjA5h0HAo5IROTMZFYrnGKZ9hxUvJ6lxQmLGb68cE2A0YiInLnMaoVTLNMqSSS7cdbE2Pai97dwpPJEgNGIiJwZJag4iUN8mTNJotGo4f0ZO2IA4PUGX1ukB3dFJHMpQcUJZ8BaUG258dKmXtRr736kKecikrGUoOJkUjXzllw4sYzePZqmnC/+QFPORSQzKUHFSVzuPTN7ULm5OQkP7v7p7TWqci4iGUkJKk400tSQp/uCha25evo48vwHjbfurGCDloQXkQyUua1wCmRaNfOWdO9WxKzy0bH9OQtWBxiNiMjpydxWOAXiZ/Fl2nNQyW6+/PzY9vI129lz4GiA0YiInLrMboU7WCQSP0kis7+aoQN6MXn8MMCrcv6nt/TgrohklsxuhTtY/BBfyDL/q7klrhc1f+kGKo/XBBiNiMipyfxWuAMl9KByM/+rOXf0IMoG9wWgIRzh1Xf14K6IZI7Mb4U7UOI088z/asyM265o6kW98s466hvCAUYkItJ+md8Kd6BINHvuQTWafsHZCWtFvb3844AjEhFpn+xohTtIfDXzTH4OKl5ubg43XnpebP+P81cnLMwoIpKusqMV7iCRDK9m3pKrp4+juDAfgL0Vx3hv9baAIxIRaZsSVJxMr2bekuKifK6feW5s/w+vv6/yRyKS9lpNUGb267jte1MeTcASlnzPkiG+RjdeOjFW/mj77oN8sGFXwBGJiLSurVb4/Ljtr6YykHQQdZm75HtbepQWcdX0cbH9F994P8BoRETa1lYr3KXGgRJX1M2uBAVwyxXnE/L/Xus27+Hj7fsDjkhEpGVttcJDzOznZvZvcduxn84IsDMlDvFlzySJRv17lzJzyqjY/u/nrQowGhGR1uW28fo347ZXpDKQdJCtkyTi3X7VpNizUCvW7WDbroOMGNI34KhERE7WaoJyzj3ZWYGkg2xY8r0tQwf04qLzRrDkQ2+q+fOvreRv77s24KhERE7WaoIyszmtve6cu6VjwwlW/JLvOVnagwK46/ryWIJa+uE2duw5xPBBfQKOSkQkUVtDfNOBncDTwFIgO7sVvkhCNfPs/asOH9SHC88bwVI/ST33ygq+qV6UiKSZtiZJDAC+A5wL/Ay4GjjonHvbOfd2qoPrbNlWzbw1d147Jba9xO9FiYikk1ZbYedcxDn3qnPuXuAiYDPwlpn9VadE18myrZp5a0YM6cvUc8ti+8+/phl9IpJe2myFzazAzD4F/Ab4MvBz4A+pDiwIibP4sjtBQVIv6oMt6kWJSFppq9TRk8BiYDLwfefcVOfc/3bO7e6U6DpZ/BBftpU6as7IYf0onzAc8J7Ifvrl5cEGJCISp61W+M+BMXhljt4zs0r/p8rMKlMfXufqSkN8jWbfMDW2vXztdjbtUHUJEUkPbd2DCjnnSuN+uvs/pc657p0VZGfJxgUL2zJiSF8unjQytv+7P6kXJSLpoa0hvkIz+5qZ/cLM7jeztqalZ7RwOK4H1UUSFHi9qMZJ9R9+vIs1H2flCK6IZJi2WuEngXJgDXAD8H9SHlGAEh/U7ToJanD/nlx+4Tmx/d+9vEzrRYlI4Npqhcc75z7rnPslcAcwsxNiCkz8LL6ulKAA7rxuSuzv/PH2/axYtyPgiESkq2urFW5o3HDOhVMcS+ASqpl3kUkSjfr3LuXaGeNj+7+Zs5RIXG1CEZHO1uaChfEz94DzsnkWX1d7DirZp6+ZTEF+HgC79h9h/tINAUckIl1ZW7P4cpJm7uWeyiw+M7vOzDaa2WYze7CZ1wvM7Fn/9aVmVuYfv9rMVprZGv/PK073L3gqEquZd70E1bO0mNuvuiC2/8zcFdTWNbTyDhGR1ElZK2xmOcCjwPXAeOAeMxufdNp9wBHn3CjgYeDH/vGDwM3OuYnAvcBTqYozXrSL96AAbrn8PHp1LwbgaNUJ5ixYHXBEItJVpbIVngZsds5tdc7VA88AtyadcyveTEGAF4Arzcycc+875/b4x9cBhWZWkMJYgcQeVCiUvdXMW1OQn8c9NzY9vPvSm6s5UnkiwIhEpKtKZYIajLdUR6Nd/rFmz/EnYRwDkhcm+jTwvnOuLvkD/GezVpjZioqKijMOOH5SQFftQQFcPm0swwb2BqCuvoFnX9HDuyLS+VLZCjfXBUl+uKbVc8xsAt6w3xeb+wDn3GPOuXLnXHm/fv1OO9BG8Qmqq83iixcKhfjzWy6K7b+xeD3bdh0MMCIR6YpS2QrvAobG7Q8B9rR0jl+logdw2N8fArwI/A/n3JYUxgmAc46o65oP6jZn0rihTBrn/edzwBN/WKSHd0WkU6WyFV4OjDazEWaWD8wGkpeQn4M3CQK8B4HnO+ecmfUEXga+7ZxblMIYYyIJ959CWBavqNseZsbnbr+YkN+T/GjLXhZ/sDXgqESkK0lZgvLvKT0AvAasB55zzq0zsx+Y2S3+aY8DfcxsM/ANoHEq+gPAKOC7ZvaB/9M/VbFCUhWJLjpBItmQs3pxw8xzY/tPvrSYunpNOxeRzpHS4q/OubnA3KRj/xC3XQvc2cz7fgj8MJWxJevqz0C15K7rp7Bw5SYqj9dw6Gg1L775AbOvn9r2G0VEzpBaYp9m8DWvpKiAz9w0Lbb/0hsfsO9g1hUREZE0pJbYF+milczb44oLx3L2UG+WZEM4wuO/f1cTJkQk5dQS+9SDalkoFOL+Oy+JPROw6qNPWLJ6W6AxiUj2U0vsS6giYfpako0efhbXzJgQ23/iD4uoqa0PMCIRyXZqiX0Jlcxz9bU05zM3T6N7tyIADh+r5tlXVgQckYhkM7XEvq68FlR7lRQV8Pnbp8f2X357jSpMiEjKqCX2RbvwarqnYuaU0Zw7ehAAUef4xe/eIhyOBBuUiGQltcS+iJ6Dahcz44t3zSIvNweA7bsP8kctySEiKaCW2JfwoK4qSbRqUP+ezL6h6WHd515dya79RwKMSESykRKUL3G595wAI8kMN192HiP9Z6PC4Qj//vTbCcOkIiJnSgnKl1jqSD2otuTkhPjyn10WGw7duG0fL7+9NuCoRCSbKEH5Eh/UVQ+qPYYP6sOnrp4U2//tn5ayc5+G+kSkYyhB+RKrmetraa87rp7M8EHeIsgN4Qg//818zeoTkQ6hltiX8ByUhvjaLTc3h6/++ZWxob6tOyt4ft6qgKMSkWygBOVL6EFpiO+UDB/Um8/cdGFs//evreTj7fsDjEhEsoESlC+iSRJn5ObLJjJ+5EDAWyL+Z0+9qVp9InJGlKB80bjlNjRJ4tSFQiH+6rNXUFiQB8C+g5X88rl3tCyHiJw2JShfONJ0Yz9k6kGdjv69S/nSXbNi+++s3MRbyz4OMCIRyWRKUL5IJK4HpWrmp21m+Wguv3BsbP+x59/R1HMROS1qiX3xPSgtWHhm/uLTlzC4f08A6hvCPPTr16lvCAcclYhkGrXEvoQl3/Uc1BkpLMjjrz9/Nbl+QdlP9h7W/SgROWVqiX2qZt6xhg/qw32fmhHbf2vZRl5fvD7AiEQk06gl9sUP8akH1TGuvngcl01ruh/1q9+/q+ejRKTd1BL7Eob41IPqEN7aUTMpG9wX8Hqp//rEPI5V1QQcmYhkArXEvkhc/TglqI6Tn5fLN//nNZQUFQBw+Fg1P3niNRoaVK9PRFqnltgXSXhQV19LRxrQtztf+x9X0vh02Yat+/jl8ws1aUJEWqWW2Kdq5qk1efwwPnvLRbH9BUs3MmfBhwFGJCLpTi2xT9XMU+/WK85PmDTx1B/fY8W6HQFGJCLpTAnKpyXfU8/M+NJdsxg7YgDgFZV96NdvsHnHgWADE5G0pATl05LvnSMvL4dv3Xct/XqVAlBX38A/PfYKeyuOBRyZiKQbJShfVD2oTtOjtIi/+9INsZl9lcdr+OH/fVnTz0UkgRKUL6EHpUkSKTd0QC++c//15PnlkPYdrORHj72iNaREJEYtsS8al6BCGuLrFOecPSBh+vnmTw7wL796VYVlRQRQgoqJ70FpiK/zXHT+2fzFHTNj+2s37eGnT7xOOKwHeUW6OiUoX+JzUOpBdabrZk7gz26aFttf+dEOfvab+Qn3BUWk61GC8sUvWKhSR53v01dP5lNXTYrtL35/C//22wVKUiJdWEpbYjO7zsw2mtlmM3uwmdcLzOxZ//WlZlbmH+9jZgvM7LiZ/SKVMTZSNfPg/dlN07h+5rmx/YUrNvHIU/MTlkIRka4jZS2xmeUAjwLXA+OBe8xsfNJp9wFHnHOjgIeBH/vHa4HvAn+TqviSqRZf8MyM+z49g6svHhc7tmjVZh556k0lKZEuKJUt8TRgs3Nuq3OuHngGuDXpnFuBJ/3tF4Arzcycc9XOuXfxElWnCKuaeVrwluiYxbUzJsSOLX5/Cw/9+nVVQBfpYlLZEg8Gdsbt7/KPNXuOcy4MHAP6tPcDzOx+M1thZisqKirOKFj1oNKHmfGFOy9JGO5b8uE2/umxuXpOSqQLSWVL3NxUuOT1FdpzToucc48558qdc+X9+vU7peCSJcziU4IKXONw382XnRc7tubj3XzvF/9N5XFVnBDpClLZEu8ChsbtDwH2tHSOmeUCPYDDKYypRfELFoY0SSItmBn33jade25smoK+ZWcFf/+zP7L/UGWAkYlIZ0hlS7wcGG1mI8wsH5gNzEk6Zw5wr799BzDfBbSKnYb40pOZccc1k7n/zpmx7vbuA0d58KEX2bRjf6CxiUhqpawl9u8pPQC8BqwHnnPOrTOzH5jZLf5pjwN9zGwz8A0gNhXdzLYDDwGfM7NdzcwA7FAJ08yVoNLOtZdM4Oufuzr236byeA3f/fkclqzeGnBkIpIquam8uHNuLjA36dg/xG3XAne28N6yVMaWLKoeVNqbMWkkvbsX8y+/epXjJ+poCEf46RPzuOemaXzqqkmYqQKISDZRS+xTNfPMMG7kQP7567czoG93wJtR87s/LeOhJ9+gtq4h2OBEpEOpJfZFIprFlykG9e/JP3/9dsaPHBg7tvj9LXznkZc0eUIki6gl9iVWM9fXku66dyvie395U8IDvTv2HOJvfvICy9ZsDy4wEekwaol9idXM9bVkgtzcHO6/ayb/a/alsV7vidp6fvyrV3nypfe0ZIdIhlNL7IsmDPHpZnsmuWr6OH74lVvp26tb7NicBav5+5//kb0VxwKMTETOhBIUEI1GY+UrDD2om4nGlJ3FT795B1PGD48d27TjAH/9kxeYv2QDAT1eJyJnQC0xifefQrr/lLFKSwr59v3X8dmbL4z9klFX38CjT7/FT5+Yx7EqlUgSySRqjUmcwafl3jObmXH7VZP456/dxqB+PWLHl3y4ja/+87O8u2qzelMiGUIJiuRnoHT/KRuMGt6ff/3mHQlrS1VV1/Lwk2/wr0/M4/Cx6gCjE5H2UIIisYqEnoHKHoUFeXzp7kv5uy/eQJ+eJbHjSz/cxld+9CxzF67RkvIiaUytMYl1+PQMVPaZPH4Yjzx4d0Jvqqa2nsd/v4hvqeisSNpSa0xiJXP1oLJTcVE+X7r7Uv7xyzcn3JvaurOCBx96kX/77QIN+4mkGbXGJC73rh5Udps4ZjAPfesuZt8wldzcpgkxby3byAM/fIYX5q2irl41/UTSgVpjknpQegYq6+Xl5XDntVN45MG7mHpuWex4XX0DT7+8jC//76eZt+gjVaIQCZhaY0i4Ua7noLqOgf168OAXruN7f3kTQwf2jh0/UnmCXz63kK//y3O8s2KTJlKIBEStMRAOqw5fV3be2CH8n2/ewRfvmkWv7sWx43sqjvHIU2/y1R89q0QlEgC1xiQWitU9qK4pJyfENTPG8+h37+EzN11IcWF+7LXGRPVX//QM8xZ9RH1DOMBIRboOtcYkPairBNWlFeTn8amrJ/Ef3/sMd1w7haK4RLXvYCW/fG4hX/r+b3lh3iqqqmsDjFQk+6V0yfdMEVElCUnSrbiAe26Yys2Xncd/v/Uhc99ew4naegCOVdXw9MvLeOG1lcwqH82Nl05k+KA+AUcskn2UoEge4lMtPmnSmKhuu+J8Xn9vPf+94MPY81IN4QhvLtnAm0s2cM7ZA7h2xnguOv9s8vP0z0qkI+hfEslDfOpBycmKCvO55fLzuWHmuby7ajN/ensN23YdjL2+Yes+Nmzdx+O/X8Ss8tFcceE5jBjSN8CIRTKfEhSqZi7tl5ubw2XTxnLp1DFs2LqPlxeuZemH22Iz/I6fqGPuwrXMXbiWssF9uWzqGGZMHknvHiVtXFlEkilBoWrmcurMjHEjBzJu5ECOVJ7gzSUbeGPxeiqOVMXO2b77IL/efZAnX1rMhNGDmDllNNMmltG9W1GAkYtkDiXlDxCFAAARNElEQVQowMVVktCDunKqenUv5o5rJvOpqy5gzaY9zF+6gaWrt9HgV6JwwNpNe1i7aQ+/fHYh40cNZPr5I5k6cTh9enZr/eIiXZgSFKpmLh0jFApx/tghnD92CNU1dbz3wVbeWbmJdZv20PgrUNS5WLL6zxfe4eyh/Zh67nCmjB/O2UP7YqYevEgjJSgSZ/HpOSjpCCVFBVw1fRxXTR/H4WPVLFq1hfdWb2Xjtn0J523dWcHWnRU8+8oKuncr4oJzhjB53DDOHTM4oaqFSFekBEViqSP1oKSj9e5Rws2Xn8fNl5/HoaPHWbJ6G8vXbmfd5r0J5ZMqj9ewcMUmFq7YBMDQAb2YOGYwE0YNYvzIgbp3JV2OEhTJkySUoCR1+vTsxo2XTuTGSydSXVPH++t3smLtDlZv3EXl8ZqEc3fuO8LOfUeYu3At4CWsc84ewLizBzKm7CwG9O2uIUHJakpQaMl3CUZJUQGXTB7FJZNH4ZxjyycVrFr/CWs+3s3G7fsTHn+ApoT1+uL1AHTvVsToYf0ZNbwfo4b1Z+TQfvQoVS9LsocSFImTJNSDkiCYGaOG92fU8P7cdV05tXUNrN+6j3WbdrNuy142f1JxUjX1yuM1rPxoBys/2hE71qdnCWcP6cfwwX0oG9SHssF91NOSjKUEReKChboHJemgsCCPSeOGMmncUABq6xrYtOMAG7btY+O2fWzctj9WGzDeoaPVHDpazfK122PH8vNyGTqgF0MH9mbogF4MGdCLwf17clafUkL6hUzSmBIUiT0oPQcl6aiwII+JYwYzccxgAJxz7Kk4xuYdB/h4+3627Kxg++5DsWev4tU3hNmys4ItOysSjufkhBjYtwcD+/VgUP8eDOjbgwF9uzOgXw/69ixR8pLAKUEB0Uj8ku8aCpH0Z2YM7t+Twf17cunUMYBXsmv3gaNs23WQ7bsPeT97Dp00+aJRJBJl1/4j7Np/5KTXQqEQfXuW0L9PKf17d6df727061VK317d6N2zhD49SigsyEvp31FECQpVM5fskJMTYtjA3gwb2JtLpzYdP1ZVwyd7D/PJ3sPs3n+U3QeOsGvfUY5WnWjxWtFolAOHqzhwuArY0+w5JUUF9O5ZQu/uxfTq4SWtHqVF9OpRTK/SYrqXFtGztIjiwnzdA5PTogRF0pLvGuKTLNOjtIiJpU3Dg42qa+rYe+AY+w5WsqfiKPsOVrLvYCX7D1a2mrzi319dU8fOvYdbPS8nJ0SPbkV071ZE95JCSrsV0qNbId2KCyktKaB7SRHFRfmUlhRQUtT4k09urn5Z7OqUoNCS79I1lRQVxGYOJqurb+DA4eNUHK7iwKEqDh6pouLocQ4eOc6hI9Ucrqw+aRp8SyKRKIePVcfW0Wqvgvw8igvzKCkqoLgon5KifAoLvD+LCvIoLMyjqMDbbtwvzPe2CwryKMjPpTA/l8KCPPJyc9SLy0ApTVBmdh3wMyAH+JVz7l+SXi8A/h8wBTgE3O2c2+6/9m3gPiACfMU591qq4tQ0c5FEBfl53sy/Ab2afd05x9GqGo5WnuDQsWqOHKvm8LETHKuq4WjVCQ4fq6byeC1Hq2qoq284rRjq6huoq2/gSGXbvbm2GJCfn0d+Xg4F+bnk5+Yk7eeSl5dDfl4Oebk55Ofl+n/mkJvrHYv/yc0JkZubQ25uqGk/p/F4iBx/Oycn5P0Z8o+HmraVMNuWsgRlZjnAo8DVwC5guZnNcc59FHfafcAR59woM5sN/Bi428zGA7OBCcAg4A0zG+OcO3mKUgeIxE+S0IKFIm0yM3p1L6ZX9+I2F2asrWugsrqWyqoa78/jNVRV11FVXUvViVqqquuoPlFH1YlaTtTUc/xEHSdq6nCtXvXUOJoSXtWpdeRSxvBmDefm5JATMkIhI8dPZqGQ+cnMCIVCmL8fChkh884LmfeeUMgwjJwcI2TeOWYQMsNCifuhUChu2zDz3tt0jve6Wdw1Yj+A/6dh3Hbl+RTkp3aiTCp7UNOAzc65rQBm9gxwKxCfoG4F/tHffgH4hXm/VtwKPOOcqwO2mdlm/3rvpSJQTZIQSZ3CgjwKC/Lo37u03e9xzlFT20B1TR0nauuprqmnpraemtoGTtTWU1Pn/Vlb20BNnbdfVxempq6e2vowdXUN1NY3UFcfprY+TLiZ6fdBc3jDn+0dKk03N146kYL81H5GKhPUYGBn3P4u4MKWznHOhc3sGNDHP74k6b2Dk96Lmd0P3A8wbNiw0w5US76LpBczo7gon+KijmkBo9Eo9Q0R6urD1DWEqasPU18fpr4hTH3YO94QjlDv/9kQjlDf4G2Hw9HYfjgSjb0eiURjr4UjEcKRqPfT+FokSjgSIRJx3p9R5yWkaOYmpXidMUKZygTVXPjJvfaWzmnPe3HOPQY8BlBeXn7aIwJDzurJ+JEDiUQdPUq1xIFItgmFQhQWhNLq2a1oNEok4ohEvWQWjUZjSSzqEv90jftRl/Ba43406ohGvWPOeT0zh3fc+ec4/7Wo867jnFeHtPG8xpqk3vWiOL9FbdwH/Pd7L+R1wizLVCaoXcDQuP0hnPxAReM5u8wsF+gBHG7nezvMXdeVc9d15am6vIjISUKhEKEQ5KHbCi1J5ZS15cBoMxthZvl4kx7mJJ0zB7jX374DmO+cc/7x2WZWYGYjgNHAshTGKiIiaSZlPSj/ntIDwGt408yfcM6tM7MfACucc3OAx4Gn/EkQh/GSGP55z+FNqAgDX07VDD4REUlP5lxHTuYMTnl5uVuxYkXQYYiISBvMbKVzrs37KnoqVURE0pISlIiIpCUlKBERSUtKUCIikpaUoEREJC1lzSw+M6sAdpzhZfoCBzsgnFTKhBghM+JUjB0nE+JUjB3nTOMc7pzr19ZJWZOgOoKZrWjP1McgZUKMkBlxKsaOkwlxKsaO01lxaohPRETSkhKUiIikJSWoRI8FHUA7ZEKMkBlxKsaOkwlxKsaO0ylx6h6UiIikJfWgREQkLSlBiYhIWlKCAszsOjPbaGabzezBoONpjpk9YWYHzGxt0LG0xMyGmtkCM1tvZuvM7KtBx5TMzArNbJmZrfZj/H7QMbXGzHLM7H0z+1PQsTTHzLab2Roz+8DM0nI5ATPraWYvmNkG///N6UHHlMzMxvrfYeNPpZl9Lei4kpnZ1/1/N2vN7GkzK0zp53X1e1BmlgN8DFyNt5LvcuAe59xHgQaWxMxmAceB/+ecOzfoeJpjZgOBgc65VWZWCqwEbkun79LMDChxzh03szzgXeCrzrklAYfWLDP7BlAOdHfO3RR0PMnMbDtQ7pxL24dLzexJ4B3n3K/8xVOLnXNHg46rJX6btBu40Dl3psUHOoyZDcb79zLeOVfjr9k31zn361R9pnpQMA3Y7Jzb6pyrB54Bbg04ppM45xbiLeqYtpxze51zq/ztKmA9MDjYqBI5z3F/N8//Scvf0sxsCHAj8KugY8lUZtYdmIW3OCrOufp0Tk6+K4Et6ZSc4uQCRWaWCxQDe1L5YUpQXgO6M25/F2nWqGYiMysDJgFLg43kZP6w2QfAAeB151zaxeh7BPhbIBp0IK1wwDwzW2lm9wcdTDPOBiqA//KHSn9lZiVBB9WG2cDTQQeRzDm3G/gp8AmwFzjmnJuXys9UggJr5lha/kadKcysG/B74GvOucqg40nmnIs45y4AhgDTzCzthkzN7CbggHNuZdCxtGGGc24ycD3wZX8oOp3kApOB/3DOTQKqgbS8zwzgD0HeAjwfdCzJzKwX3ujSCGAQUGJmn03lZypBeT2moXH7Q0hxtzWb+fd1fg/81jn3h6DjaY0/1PMWcF3AoTRnBnCLf4/nGeAKM/tNsCGdzDm3x//zAPAi3pB5OtkF7IrrJb+Al7DS1fXAKufc/qADacZVwDbnXIVzrgH4A3BxKj9QCcqbFDHazEb4v73MBuYEHFNG8icgPA6sd849FHQ8zTGzfmbW098uwvtHtyHYqE7mnPu2c26Ic64M7//J+c65lP62eqrMrMSfDIM/bHYNkFazTJ1z+4CdZjbWP3QlkDaTdppxD2k4vOf7BLjIzIr9f+tX4t1nTpncVF48Ezjnwmb2APAakAM84ZxbF3BYJzGzp4HLgL5mtgv4nnPu8WCjOskM4M+BNf49HoDvOOfmBhhTsoHAk/5MqRDwnHMuLadwZ4CzgBe9topc4HfOuVeDDalZfwX81v8FdCvw+YDjaZaZFePNJv5i0LE0xzm31MxeAFYBYeB9UlzyqMtPMxcRkfSkIT4REUlLSlAiIpKWlKBERCQtKUGJiEhaUoISEZG0pAQlWcnMzjKz35nZVr8Mz3tmdvtpXqusM6vI++V4xnfW5/mf+TV/mvOpvu+RxuoRZvYjv8r18KRz3vCrEIicEiUoyTr+Q4QvAQudc2c756bgPew6JNjI2mZmOc65v+joCvDmae3f+9fwin+eyjV7Axf5hYxxzn0H+E/g60mnPgX85alcWwSUoCQ7XQHUO+f+b+MB59wO59y/QWxNqP/y1zF638wu94+Xmdk7ZrbK/2mzjIuZfcXMPjKzD83sGf9Yt7jrf2hmn/aP3+MfW2tmP467xnEz+4GZLQWmm9lbZlYe99o/mbd+1RIzO8s/PtLfX+6/93gzsZWZt/7Rv+M9XDnUzP7DzFZY3FpYZvYVvNpqC8xsgX/sGr/XucrMnvfrKya7A0h+MPcVYLb/IHSjOXgVEkROiRKUZKMJeA1yS74M4JybiNdwPmnewmsHgKv94qd3Az9vx2c9CExyzp0HfMk/9l28Ss8T/ePzzWwQ8GO85HkBMNXMbvPPLwHWOucudM69m3T9EmCJc+58YCHwBf/4z4CfOeem0nrtyLF4a4hN8pdv+DvnXDlwHnCpmZ3nnPu5f43LnXOXm1lf4O+Bq/zvYgXwjWauPQNvza94s4FeeGVwAHDOHQEKzKxPK3GKnEQJSrKemT3q90CW+4cuwRt2wjm3AdgBjMFbG+o/zWwNXjXp9twH+hCvjM5n8cq/gFff79HGE/wGeirwll9oMwz8Fm+dIoAIXoHd5tQDjaWYVgJl/vZ0mipe/66V+HYkLcZ4l5mtwitTM4Hm/44X+ccX+SWr7gWGN3PeQLylLIBYJe7PA38DfCbp3AN4vTSRduvytfgkK60DPt2445z7st8raFySvLklVsC7d7IfOB/vl7fa5BPM7L/w1rna45y7AW9BwVl4SyR818wm+NdPriHW0mcC1DrnIi281uCa6pFFOPV/s9VxsY/ASx5TnXNHzOzXQHNLdhveOlltDcvVJL3/buBNvAUWv2VmRc65Gv+1Qv98kXZTD0qy0Xyg0Mz+V9yx+AkAC/F/wzezMcAwYCPQA9jrnIviFb2Nv48CgHPu8865C5xzN/iTDoY65xbgLSzYE+gGzAMeaHyPP4NtKd6QWl///sw9wNtn8HdcQlMSnt3O93THS1jH/HtZ18e9VgWUxl17hpmN8uMv9r+nZOuBUXH7XwEe9pPS63hJu3HSygBgezvjFAGUoCQL+T2O2/ASwjYzWwY8CXzLP+XfgRx/KO9Z4HPOuTr/+L1mtgRvyK/65KsnyAF+41/nfbzG+SjwQ6CXPxliNd69nb3At4EFwGq8NX/+eAZ/za8B3/D/bgOBY229wTm32o9zHfAEsCju5ceAV8xsgXOuAvgc8LSZfYiXsM5p5pIv41XYx8wuAQ7FrQTwFE3DfFPw7qOFT7qCSCtUzVwkA/nPLNU455yZzQbucc7dGkAc7wI3+Ym5pXN+Bsxxzr3ZeZFJNtA9KJHMNAX4hT98dhT4nwHF8dd4Q6QtJii8GYpKTnLK1IMSEZG0pHtQIiKSlpSgREQkLSlBiYhIWlKCEhGRtKQEJSIiaen/A8MZIg+Bwr+ZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from thinkbayes2 import MakeGammaPmf\n", "\n", "xs = np.linspace(0, 8, 101)\n", "pmf = MakeGammaPmf(xs, 1.3)\n", "thinkplot.Pdf(pmf)\n", "thinkplot.decorate(xlabel='Goal-scoring rate (λ)',\n", " ylabel='PMF')\n", "pmf.Mean()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "suite = Soccer2(pmf);" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.310359949002256" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPmZmsJIQt7Aio7Mq+qOCKVdxAilZsq2i12ta9vy7ab23V1n6r37ZWW6t1ty6g1WpxxQXQYkX2HQRkjWxhC9kzy/n9MTeTmyGQAJncSXLer1deuffOvXdOBnJPnuc+9zyiqhhjjDHJxud1AMYYY0xNLEEZY4xJSpagjDHGJCVLUMYYY5KSJShjjDFJyRKUMcaYpGQJyhhjTFKyBGWMMSYpWYIyxhiTlAJeB1Bf2rVrpz169PA6DGOMMbVYuHDhblXNrW2/JpOgevTowYIFC7wOwxhjTC1EZHNd9rMuPmOMMUnJEpQxxpikZAnKGGNMUmoy96CMMeZYBYNB8vLyKCsr8zqUJiE9PZ2uXbuSkpJyVMdbgjLGGEdeXh7Z2dn06NEDEfE6nEZNVdmzZw95eXn07NnzqM5hXXzGGOMoKyujbdu2lpzqgYjQtm3bY2qNWoI6BpFIhFAo7HUYxph6ZMmp/hzrZ2ldfEdh2Ve7eObDVcz9ai+qyg3n9uK6cQO9DssYY5oUS1BHYPnGfO5+cR6b9lZvsv5lxjrW5u3nf783Bp/PGqXGmKO3c+dO7rjjDubOnUvr1q1JTU3lZz/7GRMnTvQ6tAZnV9Mj8KsaklOlGSvz+faDMygostE/xpijo6pceumlnHHGGWzYsIGFCxcybdo08vLy6nR8ONy0bjlYC6qO1mzdw0ZXcurdPpPxI7rzr7kb2bAnun3NzhKueehjXv+fC6wlZUwjN+m2xxN27tcf/kGN22fOnElqaio/+EHV6927d+eWW24hHA5z5513Mnv2bMrLy7npppu48cYbmT17Nvfeey+dOnViyZIlvPvuu4wbN44xY8Ywd+5cBg0axLXXXsuvf/1rdu3axUsvvcTIkSOZN28et99+O6WlpWRkZPDss8/Sp08fnnvuOaZPn05JSQlfffUVEydO5MEHH+Tpp59mxYoVPPTQQwA8+eSTrF69mj/96U8J+5zsKlpHz36wKrZ8fNt0Xr3rAr57bn+m/XwcZ/RqE3tt494y/vbWEi9CNMY0citXrmTo0KE1vvb000+Tk5PD/PnzmT9/Pk8++SQbN24EYN68edx///2sWhW9Tq1fv57bbruNZcuWsWbNGl5++WXmzJnDH/7wB373u98B0LdvXz799FMWL17Mfffdxy9+8YvYey1ZsoRXXnmF5cuX88orr7B161YmT57M9OnTCQaDADz77LNce+21ifw4rAVVF6FQmDlf5sfWLx5xXGw5NcXPIz86m9sem8Una/cC8I//bGLCqSfQrX1Og8dqjGk6brrpJubMmUNqairdu3dn2bJlvPbaawAUFBSwbt06UlNTGTlyZLVnjXr27MnJJ58MwIABAxg7diwiwsknn8ymTZtix0+ZMoV169YhIrHEAzB27FhycqLXr/79+7N582a6devGOeecw9tvv02/fv0IBoOx90gUS1B18O/Pv6I4qACk+WHymX0P2ue3U07jwnvfpbAiQkVYuev5ubz40/MbOlRjTD05VDdcIg0YMIDXX389tv7oo4+ye/duhg8fznHHHcdf/vIXzj+/+nVl9uzZtGjRotq2tLS02LLP54ut+3w+QqEQAHfffTdnn302b7zxBps2beKss86q8Xi/3x875vrrr+d3v/sdffv2TXjrCayLr05e/3xDbHnU8W3ITD+4bEd2Zhp3XDIgtr5iWxGvzl7TIPEZY5qGc845h7KyMh577LHYtpKSEgDOP/98HnvssVhLZ+3atRQXFx/1exUUFNClSxcAnnvuuTodM2rUKLZu3crLL7/MlVdeedTvXVeWoGqxbXcha3ZU/Sf4ztkHt54qfXNMbwZ1zYqtP/LeagpLyhManzGm6RAR3nzzTT755BN69uzJyJEjmTJlCg888ADXX389/fv3Z+jQoZx00knceOONsZbN0fjZz37GXXfdxejRo49o9N+3vvUtRo8eTevWrY/6vetKVDXhb9IQhg8fromYsPB3U7/g1XnRIZ4ds1N4/77xh90/L7+QSQ98QLnz733FqK7cNXlUvcdljKl/q1evpl+/fl6HkdQuvvhi7rjjDsaOHVun/Wv6TEVkoaoOr+1Ya0HV4uMV22PL5w/qXOv+XXOz+ebIbrH1txZ9TVFpRUJiM8aYhrJ//3569+5NRkZGnZPTsbIEdRg79haxpyTaFPIJXH3ugFqOiLrpksFkpkRrUJUElcfeXpqwGI0xpiG0atWKtWvX8s9//rPB3tMS1GF8uWVvbLlVhp+2ORl1Oi4rI5WLh3SJrb85fyslZcHDHGGMMSaeJajDWLdtf2y5XVbaYfY82C0TBpMeiLaiioPK39+xVpQxxhwJS1CHsTm/MLbcqVX6ER2bnZnGRYOr7ln9a95WysqtFWWMMXVlCeowvt5bElvu1i77iI+/dcJg0vzRVlRhRYSnZ6yot9iMMaapswR1GLsOVBWH7dmx5REfn5OVzvkDO8bW35y3hUgkUi+xGWOarjfeeAMRYc2axD/sv2nTJl5++eXY+oIFC7j11lsT/r51YQnqMPYVV3XJ9e56dA+l/ejigTiNKPKLQ7w/f2N9hGaMacKmTp3KmDFjmDZt2kGv1feUGvEJavjw4TzyyCP1+h5Hy2rxHUJRaUWs/p4AvTq3OqrzdGyTxbDuOczbVADAC7PWcuGoE+orTGNMglz/VP0/+F/pqesP/YxqUVERn332GbNmzWL8+PHcc889B02psWrVKn7zm9/w0ksv0a1bN9q1a8ewYcP4yU9+wldffcVNN91Efn4+mZmZPPnkk/Tt25drrrmGli1bsmDBAnbs2MGDDz7IZZddxp133snq1asZPHgwU6ZMYciQIfzhD3/g7bff5p577mHLli1s2LCBLVu2cPvtt8daV5deeilbt26lrKyM2267jRtuuKHePydLUIewesue2HJWqo/0tIPr79XV98cNYN7j/wVgzc5i1m7dQ+9ubY85RmNM0/Pmm28ybtw4evfuTZs2bVi0aBEQnVJjxYoV9OzZkwULFvD666+zePFiQqEQQ4cOZdiwYQDccMMNPP744/Tq1YsvvviCH/3oR8ycOROA7du3M2fOHNasWcP48eO57LLL+P3vfx9LSBAtPuu2Zs0aZs2aRWFhIX369OGHP/whKSkpPPPMM7Rp04bS0lJGjBjBpEmTaNu2fq9r1sV3CGu/3hdbbpuVekznGtGnE91bR4epK8Lj79pgCWNMzaZOncrkyZMBmDx5MlOnTgWoNqXGnDlzmDBhAhkZGWRnZ3PJJZcA0dbXf//7Xy6//HIGDx7MjTfeyPbtVdVwLr30Unw+H/3792fnzp11iueiiy4iLS2Ndu3a0b59+9hxjzzyCIMGDeKUU05h69atrFu3rt4+g0rWgjqEzTurhph3yDmyIeY1+daY4/m/t1YDMGfdHgpLysnOPLJnq4wxDedw3XCJsmfPHmbOnMmKFSsQEcLhMCLChRdeWG1KjUPVUI1EIrRq1YolS2qeNNU9jUZd67DWNPXG7Nmz+eijj/j888/JzMzkrLPOoqys7DBnOTrWgjqErbuLYstd27Y4zJ51c8UZfchOjX7cFWHl6fetFWWMqe61117j6quvZvPmzWzatImtW7fSs2dP5syZU22/MWPG8NZbb1FWVkZRURHvvPMOAC1btqRnz56xckSqytKlhy8SkJ2dTWFh4WH3iVdQUEDr1q3JzMxkzZo1zJ0794iOrytLUIews6Dqr4Hu7Y/8Gah4gYCf8wd1iq2/tXCrDTk3xlQzdepUJk6cWG3bpEmTqo2yAxgxYgTjx49n0KBBfPOb32T48OGxGXBfeuklnn76aQYNGsSAAQP497//fdj3HDhwIIFAgEGDBvHQQw/VKc5x48YRCoUYOHAgd999N6eccsoR/JR1Z9NtHMKZd71BQVk0gfz9hlMZ1a/2Sua1yS8oYdy97xF2PvI/Xj2csUO6H/N5jTH1ozFNt1FUVERWVhYlJSWcccYZPPHEEwwdOtTrsA5i023Us1AoTGF5Veumb7c29XLe3JxMBnWreuD3xVlr6+W8xpjm54YbbmDw4MEMHTqUSZMmJWVyOlY2SKIG67fvJ+K0cjICQk7WsQ+SqPSds3qz6B/Rlt6yvAPkF5SQm5NZb+c3xjQP8d1+TZG1oGrw5daqaTZaZx798081OXtQN9plRv8uCCs2WMKYJNNUbnskg2P9LBOaoERknIh8KSLrReTOGl5PE5FXnNe/EJEeca8fJyJFIvKTRMYZb8P2gthybstjewYqns/nY9zgqrmiPli6zQZLGJMk0tPT2bNnjyWpeqCq7Nmzh/T0o++BSlgXn4j4gUeBbwB5wHwRma6qq1y7XQfsU9UTRWQy8ABwhev1h4D3EhXjoWxxDTHv3PrYh5jHu/a8AUz7fDMhhb2lYT5atJnzhves9/cxxhyZrl27kpeXR35+vtehNAnp6el07dr1qI9P5D2okcB6Vd0AICLTgAmAO0FNAO5xll8D/ioioqoqIpcCG4DiBMZYox37SmPL3dtn1fv52+ZkMPi4lizYfACAlz9ZZwnKmCSQkpISq9ZgvJfILr4uwFbXep6zrcZ9VDUEFABtRaQF8HPg3sO9gYjcICILRGRBff7Fk19YEVvu1eXoqpjX5qqz+sSWl39dyK59DZ6HjTEmqSUyQUkN2+I7dg+1z73AQ6paVMPrVTuqPqGqw1V1eG5u7lGGWV0kEmF/aSi23qdr/Qwxj3f6wK7ktqgaLPH8R6tqOcIYY5qXRCaoPKCba70rsO1Q+4hIAMgB9gKjgAdFZBNwO/ALEbk5gbHGbNtTTMhJoyk+6FwPZY5q4vP5OH9Q9cESxhhjqiQyQc0HeolITxFJBSYD0+P2mQ5McZYvA2Zq1Omq2kNVewB/Bn6nqn9NYKwxa1zTbLTKCODzJe4jmnJe/2qTGc5ZvvXwBxhjTDOSsKuvc0/pZmAGsBp4VVVXish9IjLe2e1povec1gM/Bg4ait7Q8vZUFU1s06J+h5jHy83J5KQuVYMwrLKEMcZUSWglCVV9F3g3btuvXMtlwOW1nOOehAR3COXBqumUU1MS/xzz5aNPYOkr0WrDizbvp6i0gqyMxCZGY4xpDKySRJxQuGoch19qGsNRvy4ceTw5ac40HBF48WMbLGGMMWAJ6iDBUFVVB78v8QnK5/NxZv8OsfV3FuYl/D2NMaYxsAQVx92CCvgb5uO59rz+iDMCf+v+cpZvtKfYjTHGElSckKsuXsCf+BYUQM+OrTgxt6qi+T8+Wt0g72uMMcnMElQcdxdfoAG6+CpdekqP2PJn63ZT4RqsYYwxzZElqDihsPseVMN9PJNG9yIjEE2IJUHlzc/WNdh7G2NMMrIEFad6gmq4FlR6WgqnnFhVVumNuRsb7L2NMSYZWYKKE4o0/CCJSled0ze2vGZnCXn5hYfZ2xhjmjZLUHHC1RJUw7WgAIb26khnZ4JEBf5hBWSNMc2YJag47i6+hhwkUemCIVUFZD9avt1m2zXGNFuWoOK4u/hSAv4Gf//vju2HM1aCvaVhPllmBWSNMc2TJag4IXcliQbu4gNonZ3BwG4tY+tTP7HRfMaY5skSVJyIVrWgUht4kESlyaefGFtevKWAgqIyT+IwxhgvWYKK4+UovkrnDu1Oq/Toewcj8OJMqyxhjGl+LEHFCXtQ6iiez+fjbFcB2fcWWQFZY0zzYwkqjrtYbIq/4QdJVLrmvAFUpse8ggqWfLXTs1iMMcYLlqDiRNxdfAFvWlAA3Tvk0Kt9VQHZFz5e41ksxhjjBUtQcVwNKE+GmbtNdBWQ/Xz9Hisga4xpVixBxXE/qJvSgMViazJpTG8yU6oKyL4+Z62n8RhjTEOyBBXH1cNHIODtx5Oa4ufUE9vG1t+Yu8m7YIwxpoFZgorjnrAwxaNh5m5TxvaLLa/bVcLmnQUeRmOMMQ3H+ytwknGXvkv1uAUFMPCE9nRtVVVA9pkZK70NyBhjGoj3V+Ak464kEfB4kESli4d1iy3PXrXTCsgaY5oFS1Bx3NNtJEMXH8B3x/Yn1QmloDzCO19s8DYgY4xpAMlxBU4iYVcLKi0lOVpQWRmpDO/ZOrb+ypyvPIzGGGMahiWoOJEkqMVXk6vO6RNbXrW9iB17izyMxhhjEi95rsBJwtWAIjVJWlAAp/bvQoesFCA6FN4GSxhjmjpLUHGqV5JIro9n3ODOsWWbbdcY09Ql1xU4CbhH8SXLIIlK15w3oNpsux8t2uxtQMYYk0DJdQVOAq5KR6SlBLwLpAatszMY5J5t99P1HkZjjDGJZQkqjrpn1E2yLj6A757dO7a8LO8Au/YVexiNMcYkTvJdgT3mrsWXLMPM3c4c2I3cFtGWXVjhmQ9ssIQxpmmyBBXHPezA6+k2auLz+Rg3uEtsfcbSbTZYwhjTJFmCipPsXXwA13yjarDEvtIwMxZs8jQeY4xJhOS8AnskFAqjzkTrgiZNLb54bXMyGNo9J7Y+9dN1HkZjjDGJYQnKpcw1Y62Id9O918XVrsoSK7ZZZQljTNOT0AQlIuNE5EsRWS8id9bwepqIvOK8/oWI9HC2jxSRJc7XUhGZmMg4K4VdY8yTPXOPOblbtcoST72/wuOIjDGmfiXsOiwifuBR4AKgP3CliPSP2+06YJ+qngg8BDzgbF8BDFfVwcA44O8ikvCHkspdLShfcjegALhwaNVgiQ+XWWUJY0zTksiGwkhgvapuUNUKYBowIW6fCcDzzvJrwFgREVUtUdWQsz2d6Fx9CVcRqrrAJ3sXH0QHS7in4Xjzv/bgrjGm6UhkguoCbHWt5znbatzHSUgFQFsAERklIiuB5cAPXAkrRkRuEJEFIrIgPz//mAMuD1a9RZJVOapRTlY6I09oE1uf9h+bhsMY03Qk8jJcUxMkviV0yH1U9QtVHQCMAO4SkfSDdlR9QlWHq+rw3NzcYw446L4H1QhaUADXn1/Va7puVzFrtu7xMBpjjKk/iUxQeUA313pXYNuh9nHuMeUAe907qOpqoBg4KWGROoKuLj5/48hPDD6hAz3aRHO3Ijz1vlWWMMY0DYlMUPOBXiLSU0RSgcnA9Lh9pgNTnOXLgJmqqs4xAQAR6Q70ATYlMFYAKqoNM0/0u9Wfy07tGVues3Y3JWVBD6Mxxpj6kbAE5dwzuhmYAawGXlXVlSJyn4iMd3Z7GmgrIuuBHwOVQ9HHAEtFZAnwBvAjVd2dqFgrhdxdfI1hGJ/jW2f2oUVKNN6ykPLCx9aKMsY0fgkduq2q7wLvxm37lWu5DLi8huNeAF5IZGw1cQ8z9zeiJlRqip9zBnTgrSU7AHjjiy3ceNFgj6Myxphj0wjGqjUc9yAJfyNqQQHccMFJsWe3dhQG+WTJFm8DMsaYY2QJyiUUcj+o27gSVLf2OZzcJSu2/szHqz2Mxhhjjp0lKBf3g7q+RvjJfG9sv9jy8q+L2LyzwMNojDHm2DTCy3DiuLv4Ao0wQ505+Dg6t0wFovX5HntnmccRGWPM0Wt8V+EEqtaCalw9fDGXndYjtjx7db4NOTfGNFqWoFzCrmKrgcZQ66gG3z6rb7Uh5//4yIacG2Map8Z5FU6QYMg9zNzDQI5BeloKY0/qGFv/1xdbrMq5MaZRsgTlEgpVlQpsTA/qxvvBhSfHEuyuoiAfLtzkaTzGGHM0LEG5BMNVLahAY21CAZ3bZTPkuJax9Wc++tLDaIwx5uhYgnIJhataUP5GOIrP7YcXnhxbXrurmGVf7fIwGmOMOXKN+ypcz0LVhpk33hYUwLDeHTmhXQYQrXL+t3eWexyRMcYcGUtQLhWNcD6ow7nmnN6x5fmb9rNtd6GH0RhjzJGxBOUSdnXxBQKN/6O5aNTx5LaI1gMOKzz61lKPIzLGmLpr/FfheuQeZt7Yu/gAfD4f3zq1R2z941W7KCqt8C4gY4w5ApagXEIRVwuqkT6oG++qc/tXe3D3qffsXpQxpnFoGlfhelJ9FF/jb0FB9MHdcYM6x9bfmL+lWtV2Y4xJVpagXJrSKD63H108iFTnX7qgLMI/Pl7lbUDGGFMHlqBcwk2wiw+gbU4GZ/bNja1P/c8GK39kjEl6h70Ki8hzruUpCY/GY9VaUE0oQQHcdumQWPmj/OIQb3y23tuAjDGmFrVdhQe5lm9LZCDJIKxVLaiUJjDM3K1rbjajjm8dW39u5loPozHGmNrVdhXWWl5vUtyDJBrjhIW1uXXCICrvrG3dX87Hizd7Go8xxhxOoJbXu4rII4C4lmNU9daEReaB6l18TWeQRKW+3doyqGs2S/KiFSUef28lY4d09zgqY4ypWW3NhJ8CC4EFrmX3V5PiHiTR1Lr4Kt1yycDY8rr8Uj5bkedhNMYYc2iHbUGp6vMNFUgyqP6gbtNrQUG0iGy/Dpms3lkCwF/eXs7ok7p6HJUxxhzssAlKRKYf7nVVHV+/4XgrUm2QhN/DSBLr9gmDuPGJzwFYs7OEL1ZvY1S/zrUcZYwxDau2e1CnAluBqcAXQNNsVjjcxWL9TaCa+aGM6teZvh0yWeO0ov7876VMtQRljEkytd1o6Qj8AjgJeBj4BrBbVT9R1U8SHVxDc9+DSk1pui0ogFsurprQcLXTijLGmGRy2ASlqmFVfV9VpwCnAOuB2SJyS4NE18CaaiWJmow+qSu922fG1h9+a5mH0RhjzMFqvQqLSJqIfBN4EbgJeAT4V6ID84L7Qd3UJjqKz+3mi06KLa/eXmStKGNMUqmt1NHzwH+BocC9qjpCVX+jql83SHQNLOIeZt7EW1AAZwzsRq/cqmnh//imTWhojEketV2FrwJ6Ey1z9LmIHHC+CkXkQOLDa1hNcT6o2tw+vuq5qLW7Spi1xKpLGGOSQ233oHyqmu36aul8Zatqy4YKsqG48hNpTXiYudvok7oyoFOL2PrDb63wMBpjjKlSWxdfuojcLiJ/FZEbRKS2YemNWtg1BUWgGdyDqvTTSUMQp+zipr1lvPW5VTo3xnivtqvw88BwYDlwIfDHhEfkIXcLqqkPM3cbfEIHhh6XE1v/2/urbL4oY4znaktQ/VX1u6r6d+Ay4PQGiMkzkSY83UZtfnrZUConEd5+IMgrn3zpbUDGmGavtqtwsHJBVUMJjsVz7kZDir/5tKAgWun81BOq5ot64sMvqQiGPYzIGNPc1TphoXvkHjCwKY/ic7eg0ppRF1+lOy8fRsBpRe0rDfO3t5Z4G5AxplmrbRSfP27kXuBIRvGJyDgR+VJE1ovInTW8niYirzivfyEiPZzt3xCRhSKy3Pl+ztH+gEfCVYqv2XXxAXRrn8O4gR1j6698vpmCojIPIzLGNGcJuwqLiB94FLgA6A9cKSL943a7DtinqicCDwEPONt3A5eo6snAFOCFRMXpps28BQXw828NJzMl2owqDSkP/nOBxxEZY5qrRDYTRgLrVXWDqlYA04AJcftMIDpSEOA1YKyIiKouVtXKujsrgXQRSUtgrED1UXz+ZvKgbrzszDSuPK1HbH3G8p1s3lngXUDGmGYrkVfhLkSn6qiU52yrcR9nEEYB0DZun0nAYlUtj38D59msBSKyID8//5gDdieo9GbaggL44cWDaJsZ/flDCvdPm+9xRMaY5iiRCaqmCZX0SPYRkQFEu/1urOkNVPUJVR2uqsNzc3OPOtBK7gTVHGrxHUog4OcH5/eLrc/ftN+mhjfGNLhEXoXzgG6u9a5AfLns2D5OlYocYK+z3hV4A7haVb9KYJwARCKRatmzOT2oW5NJY3pxfNt0IFpI9n9fX2wP7xpjGlQiE9R8oJeI9BSRVGAyED+F/HSigyAg+iDwTFVVEWkFvAPcpaqfJTDGGPczPwL4fM23BQXRn/+XVwyLNXHz9lfw3IcrPY3JGNO8JOwq7NxTuhmYAawGXlXVlSJyn4iMd3Z7GmgrIuuBHwOVQ9FvBk4E7haRJc5X+0TFChAMV7UOfE13tvcjMrRXR07v1Sa2/vTH6ygsOehWoDHGJERCmwmq+q6q9lbVE1T1fmfbr1R1urNcpqqXq+qJqjpSVTc423+rqi1UdbDra1ciYy1ztaAsQVW5+zujSHee3i0OKvdPm+dxRMaY5qJ592O5BC1B1Sg3J5MrT+seW/9wxS5WbtrtYUTGmObCEpSjehefZSi3my4ZTPusFCBabeNXL82zARPGmISzBOUIhuwe1KEEAn7umjQoNmfUV7tLeeGjVR5HZYxp6ixBOcqrdfFZhop39uDujOrpqnb+0Vr2FZZ6GJExpqmzBOUIhl0Jyj6VGv1myqnVBkzc8+JcjyMyxjRldil2lFfYPaja5OZkcu1Zx8fWP1271ypMGGMSxhKUI+RqQfktQR3S9y8YyHGto3V7Fbhn2kLKyoOHP8gYY46CJShHyD1IwkZJHJLP5+O33xmJ3/mI8otD/O4VKyZrjKl/lqAc5SF7DqquBp7QnglDO8fW316ynUXrdngYkTGmKbIE5Qi5noMKWIaq1Z1XjKSD82xUROGXL84n5EryxhhzrCxBOSqsi++IpKb4ue/bw2OtzW0HKmz2XWNMvbIE5QhaC+qIjerXmfNOqqrh+/r8PBauta4+Y0z9sATlCFmpo6Nyz3dOoV1mAIiWQbrrhXk2qs8YUy8sQTnc80H5rQVVZ+lpKdz/3ZGxrr5dRUF+bQ/wGmPqgSUoR8hV/NTvtwR1JEb168TEYV1i6x+s2MnHizd7GJExpimwBOVwF4u1FtSRu+uKEXRtlQpEp4i/75VFVqvPGHNMLEE5whGNLQesGN8RCwT8/OF7p5HifHQF5RFu+/unNi2HMeao2ZXYYdNtHLu+3dpywzknxtaXfV3EY28v9TAiY0xjZgnKEQq7WlB++1iO1vcvGsTQ41rG1p+ZvcGGnhtjjopdiR3u6TYCNkjimDz0/dNpne4HokPPf/b8XIpKKzyOyhjT2FiCclRrQVmVGuQgAAAaNUlEQVQf3zHJyUrn91dXDT3fUxLmlsc+sftRxpgjYgnK4X5Q12+DJI7ZqH6duWp099j64q0HePiNRR5GZIxpbOxK7KieoKwFVR9umziUId2q7ke98Nlmez7KGFNnlqAcoYgNkqhvPp+Pv/zwTHJbREshRRR+NXUhefmFHkdmjGkM7ErsqDbdhg2SqDdZGan8+frRpDqfaXFQufHR2ZSUWb0+Y8zhWYJyuB/UTbEWVL0a0KMdPx0/ACH6GX9dUMHNj822QRPGmMOyK7HDuvgS6/Iz+lSr17doywHunzrPw4iMMcnOrsSOUMiKxSbaL789qtqgiX8tyGPa7DUeRmSMSWaWoBxhrWpBpVoLKiF8Ph+P/ugsOresKir7h+kr+WTJFo8jM8YkI7sSO6zUUcPITE/hiZvPIjs1+hmHFH7+0gKWb8z3ODJjTLKxK7EjolVdfCkB+1gSqWtuNn++7tTYyL6ykHLz3+fY8HNjTDV2JXZUL3VkH0uiDevdkd9OHkLl7b6C8gjX/2WWzSFljImxK7Ej4h7FF7BBEg3hvOE9ue2CvrHh5zsKg0x56GMrLGuMASxBxbiHmacE/B5G0rxc/Y0BfOe0HrH1LfvKueZPH1FWbg/yGtPcWYJyVHtQ17r4GtRPLh/O+CGdYuvrd5dy3cMfEwqFD3OUMaapsyuxw5WfCNggiQZ339WncW7/3Nj6yu3FXPfwR5akjGnGEnolFpFxIvKliKwXkTtreD1NRF5xXv9CRHo429uKyCwRKRKRvyYyxkohV9kdK3XkjQevG8NpJ7SOrS/NK+Kahz6iImhJypjmKGFXYhHxA48CFwD9gStFpH/cbtcB+1T1ROAh4AFnexlwN/CTRMUXz10WLtVaUJ7w+Xz89UdnMbJHTmzbim1FXPvQh5akjGmGEnklHgmsV9UNqloBTAMmxO0zAXjeWX4NGCsioqrFqjqHaKJqEOFqo/hskIRXfD4fj99yDqf0bBXbtnJ7MVP++IFVQDemmUlkguoCbHWt5znbatxHVUNAAdC2rm8gIjeIyAIRWZCff2yVCCKuUkdp1oLylM/n4283n12tu2/1zhKufHCGPSdlTDOSyCtxTQ8T6VHsc0iq+oSqDlfV4bm5ubUfcBjuWnw2zNx7ld19Z/ZuE9u2eV85kx/8kG27reKEMc1BIhNUHtDNtd4V2HaofUQkAOQAexMY0yFFbLqNpOPz+XjoxjO5ZHDH2LadRUG+88ePWb15t4eRGWMaQiKvxPOBXiLSU0RSgcnA9Lh9pgNTnOXLgJmqWucWVH1yDzNPTbEWVLLw+Xz8ZsporhrdPVZxYl9ZmGv/+imzlmz2ODpjTCIlLEE595RuBmYAq4FXVXWliNwnIuOd3Z4G2orIeuDHQGwouohsAv4EXCMieTWMAKxX1R7UtXtQSef/XTacW8b1wed0CpeFlJ/8YwH/+HClt4EZYxImkMiTq+q7wLtx237lWi4DLj/EsT0SGdtB7+daTrMWVFL63vkn06VNC371ymLKwxBWeOjd1azfXsA93z0Fn1UAMaZJsd9oR9j1HFSK3xJUsjp/xPE8+cPTyUmL/tdVhOmLt3PVHz6goKjBnkowxjQAS1AO960ve1A3uQ08oT1Tf3IuXVulxrat3F7MpP+dYYMnjGlC7ErscE0HZV18jUDndtm8due4ag/07i4JMeWRT3hp5ioPIzPG1BdLUA73PSh7DqpxSE9L4fFbx3LtGT1igycqIvB/b63m9sdn25QdxjRylqAc1sXXeN02cRj/d9VwslKr/t1mf7mHS+9/j+Ubj63CiDHGO3YlBkKhMOoUtRDUavE1QmOHdOe1n3+DE9tlxLbtKAxy7V8+5a//XkTEXQ3YGNMoWIICylyVskVsuvfGqmObLF69axwTh3aK1dAKKTw1eyPfeXAGeflWIsmYxsQSFBB0TYpnH0jj5vP5+PVVp/HnKSNolV71r7l6ZwmTHviQp99fZq0pYxoJux4D5a4WlM8aUE3CmYOP49//c0G1uaXKw8pfZqzj2w/OYOOO/R5GZ4ypC0tQQMg1xty6+JqOnKx0nrjtXP5n4km0SKn6d12zs4Rv/d/HPPDqPJtS3pgkZgkKKA+GYstWyLzpufyMPkyPa00FIzD1861ceO87VnTWmCRll2Mg6Kpz5LMWVJPUNieDJ247l998a1C1e1O7ioLc8fwCrv3Th9btZ0ySsQQFVFRUdfP4LT81aZeceiLv/vpiLh7Usdr9xsVbD3D5gx9z9/OfUVhS7l2AxpgYS1BYC6q5yUxP4bfXjOb5W86gd/vM2PaQwltLdnD+Pe/w0OsLrBKFMR6zBAWEXAlK7BNpNk7umcurd13AvZcPpE1m1cPZJUHl+TmbGXfP2zz5zlIbSGGMR+xyTPVh5n5rQTU7E07rxfv3XMKUMd3JdI32218W4dGP1nPu3dMtURnjAUtQVO/i89uDUM1SaoqfOyYNZ8Y9FzFhSCdcZf1iiWrsL6fz0OsLKCqt8C5QY5oRS1BARbUHdS1BNWfZmWnce/VpvP3LcZx/UntSXL8hBeURnp+zmW/86m3ufv4zduwt8i5QY5qBhE753li470HZrOEGoH3rFjxw3ens2lfMH/+1iFmrdlHh/DcpDSlvLdnBu0t3MOS4lnz//AGM6tfZ24CNaYIsQVG9iy9gGcq4VCaqfYWlPPrWUt5bso3iYLTySFhhweYDLHjic7rkpDJ+RDe+c05/sjJSazmrMaYu7GoMVITcw8w9DMQkrdbZGfzy26fw8W8u4cZzjie3RfW/7b4uqOCxj75i7N1vcdOjM/lsRZ5HkRrTdFgLCgi5ElTAah2Zw0hPS+GHlwzhxosG8eHCTTw/ay2rtxfHZmQuD8Nn6/fx2fovyG2xkLMHdOTKs/vQs2Orw57XGHMwS1BARdgqSZgj4/P5OH/E8Zw/4ng27yzgqfdXMHvVLgorqv7YyS8O8eq8PP45byvdWqcz9uROXH56bzq3y/YwcmMaD0tQQMRVzdxnfXzmCHXvkMNvpowmFArz9hcbeP3zDazaVkTlfytF2LKvnGc/3cRzn26iW+s0zujXkYmjT+CEzq29Dd6YJGYJCgi6WlABa0KZoxQI+Ll0dC8uHd2L/IISXp65mo+Xb2frvjLUmeNXgS37ynnxv5t58b+baZ+Vwojj23D+sO6MOakLPhukY0yMJSiqzwfltwuEqQe5OZncNnEYt02EjTv2M3XWl8xZs4ttB6o/5LurKMg7y3byzrKdpAeEPh1bMKZvR84b3oPuHXIOcXZjmgdLUFQfxRewLj5Tz3p2bMUvrhwFwFfb9vHP/6zj8y93sWVfOeraryykLM0rYmneeh79aD1tMvz069KSUb07cO6Q4+zelWl2LEFhpY5Mwzmhc2vuvGIkAPkFJUz/fD2frNjO2p3FlIW02r57S8POiMB9/OndNbTJ8HNihyyGntCO0QO6MKB7W+sSNE2aJSggEnF18dkwc9NAcnMyuW7cQK4bN5BIJMJ/luXxweItLN2yj20FFUSq5yv2loaZt6mAeZsKePzjr0gPCF1bpdOnc0sGn5DLqf060zXXWlmm6bAEBQRdVaqti894wefzcebg4zhz8HEAFBSV8dHiLfx39XZWfV3AjgNB4vIVZSFl/e5S1u8u5Z1lO+GNFbRIEbq0zqBn+yz6d2vN0F4drKVlGi1LUEDI9aeqPahrkkFOVjqTTu/NpNN7A9GE9enyPL74cier8vaTt68M1yNXMcVBZe2uEtbuKmHGil3w3pf4Bdq1SKGzk7h6d23NgO5t6detDYGA/+CTGJMkLEFh96BM8svJSueSU0/kklNPBCASibBsYz5zV29n6cY9bMwvJr8oSDi+mUW0ZuDOoiA7i4Is3noAFm4DomW9ctL95Gan0aVNBse1y+L4Tjn06daGEzu1suRlPGcJCgi7fquti880Bj6fj8EndGDwCR1i2yqCYZZu2MXCdTtZk7efTfnF7DhQftDgi0oRhX2lYfaVRltcsCf2mgBZqT5at0ghNzuNjq0y6Ny2BcflZtOtfTY9O+SQk5We4J/SNHeWoICwdfGZJiA1xc+IPp0Y0adTte15+YUsXr+TlVv2smlnIdv2l5JfWEHpIRIXRB8oLqyIUFhRzpZ95bDlwMHv5xey0nzkpKfQqkUq7bLTaNsynY6tM+jUugUd2rSgS7ts2man2z0wc1QsQVF9PihLUKap6ZqbTdfcbC45tfr2/IISVm3czbpt+9m48wBf7y1h14Fy9pUED5u8KlWElb0lYfaWhGFv2SH38wlkBITMVD8t0gJkpQfIyUyhZUYqrbJSaZuVTk5WGu1aptM2O4M2LTNo1zKd9LSUY/3RTSNnCYrqgyRSApagTPOQm5NZbeSgW2FJOWvz9rFhRwGbdxayfV8xuw6UsbeoggOlIYqDkYOGwR9KRKODN4qDIfKLQ3WOLyCQGhDSAj7SU3xkpPhJd5JcZqqfjLQALdJSyEoPkJWeQlZmKlnp0cSXlZlKi4zocssWaWSmBawV1wglNEGJyDjgYcAPPKWqv497PQ34BzCMaAf4Faq6yXntLuA6IAzcqqozEhVnyCYsNKaa7Mw0hvXuyLDeHWt8PRKJsDW/kLz8QrbuLmLH3mJ27i9lX1E5+0sq2F8SpKQiTElFhDo0xmoUUggFlZJgGErDQPCofx5B8Yvg9wkpfiHggxS/j4BfSPX7SPH7SAn4SA0IKX4fqQE/qQEfaQEfKc5yaoqfNOd7qt/5nuInLcV5PeAnLeAnNdUfPT7FT0rAR1qKnxS/n7TU6HEpgWiitYRZu4QlKBHxA48C3wDygPkiMl1VV7l2uw7Yp6onishk4AHgChHpD0wGBgCdgY9EpLeqhkmA6vegbJCEMbXx+Xx075BTp3qBBUVlbN9XzM49xew6UMrughL2FVWwr6icA6VBCsuClJSHKS4PUR6KUBZSguFIrMBufVAkmvDCSnlsUFQN4/QbkKCICD6i3aAi4BNxlqPfK5dFwC/i2l59v2rbfYJUntMn+Jzjo6/7XMvR8wnRWRyq3qPqPavtI0Dl6wi/mDyC7My0hH5GiWxBjQTWq+oGABGZBkwA3AlqAnCPs/wa8FcREWf7NFUtBzaKyHrnfJ8nItCwdfEZkzA5WenkZKXTt1vbOh8TiUTYV1jOngMl7CsqZ29hGYUlFRwoqaCwtIKishBFZUFKy0OUVIQorQhTHoxQFgpTEVKCoQgV4QihiBKMUOfuyIakCKpOmtSqrdW/J687JoYadYLqAmx1recBow61j6qGRKQAaOtsnxt3bJf4NxCRG4AbAI477uB+9LoKWQvKmKTi8/lom5NB25yMejlfKBSmuCxIYWmQorIKikuDFJcFKS0PUloRprg8SHkwTGl5iIpgmPJQhPKKEBWhCBWhCMFwhPJgOJrwnPVwRAmFlVBECUWi6+EIhCMRIhEIqxKOKBGFiFZ+j6aeZEyYR0oa4JGcRCaomqKP/2c51D51ORZVfQJ4AmD48OFH/U/evV0L9pdUEIlAh1YtjvY0xpgkFQj4ycnyJ9WzW6FQOJr8QuFo8gtHE18wFCEUiX4PO98jESUYChMOKyFVQqGwkxBdX+EIYVUiqoRCERQnQYaVsCpamSSd/VWjCTW2n5M1o8tVg2AiWvVaZbIFyEhN/Bi7RL5DHtDNtd4V2HaIffJEJADkAHvreGy9uf/aMYk6tTHG1CgQ8DvVOmw4/aEk8obLfKCXiPQUkVSigx6mx+0zHZjiLF8GzFRVdbZPFpE0EekJ9ALmJTBWY4wxSSZhLSjnntLNwAyiw8yfUdWVInIfsEBVpwNPAy84gyD2Ek1iOPu9SnRARQi4KVEj+IwxxiQnUW0Cd+uI3oNasGCB12EYY4yphYgsVNXhte1nY6qNMcYkJUtQxhhjkpIlKGOMMUnJEpQxxpikZAnKGGNMUmoyo/hEJB/YfIynaQfsrodwEqkxxAiNI06Lsf40hjgtxvpzrHF2V9Xc2nZqMgmqPojIgroMffRSY4gRGkecFmP9aQxxWoz1p6HitC4+Y4wxSckSlDHGmKRkCaq6J7wOoA4aQ4zQOOK0GOtPY4jTYqw/DRKn3YMyxhiTlKwFZYwxJilZgjLGGJOULEEBIjJORL4UkfUicqfX8dRERJ4RkV0issLrWA5FRLqJyCwRWS0iK0XkNq9jiici6SIyT0SWOjHe63VMhyMifhFZLCJvex1LTURkk4gsF5ElIpKU0wmISCsReU1E1jj/N0/1OqZ4ItLH+Qwrvw6IyO1exxVPRO5wfm9WiMhUEUnoFMXN/h6UiPiBtcA3iM7kOx+4UlVXeRpYHBE5AygC/qGqJ3kdT01EpBPQSVUXiUg2sBC4NJk+SxERoIWqFolICjAHuE1V53ocWo1E5MfAcKClql7sdTzxRGQTMFxVk/bhUhF5HviPqj7lTJ6aqar7vY7rUJxr0tfAKFU91uID9UZEuhD9femvqqXOnH3vqupziXpPa0HBSGC9qm5Q1QpgGjDB45gOoqqfEp3UMWmp6nZVXeQsFwKrgS7eRlWdRhU5qynOV1L+lSYiXYGLgKe8jqWxEpGWwBlEJ0dFVSuSOTk5xgJfJVNycgkAGSISADKBbYl8M0tQ0QvoVtd6Hkl2UW2MRKQHMAT4wttIDuZ0my0BdgEfqmrSxej4M/AzIOJ1IIehwAcislBEbvA6mBocD+QDzzpdpU+JSAuvg6rFZGCq10HEU9WvgT8AW4DtQIGqfpDI97QEBVLDtqT8i7qxEJEs4HXgdlU94HU88VQ1rKqDga7ASBFJui5TEbkY2KWqC72OpRajVXUocAFwk9MVnUwCwFDgMVUdAhQDSXmfGcDpghwP/NPrWOKJSGuivUs9gc5ACxH5biLf0xJUtMXUzbXelQQ3W5sy577O68BLqvovr+M5HKerZzYwzuNQajIaGO/c45kGnCMiL3ob0sFUdZvzfRfwBtEu82SSB+S5WsmvEU1YyeoCYJGq7vQ6kBqcC2xU1XxVDQL/Ak5L5BtagooOiuglIj2dv14mA9M9jqlRcgYgPA2sVtU/eR1PTUQkV0RaOcsZRH/p1ngb1cFU9S5V7aqqPYj+n5ypqgn9a/VIiUgLZzAMTrfZeUBSjTJV1R3AVhHp42waCyTNoJ0aXEkSdu85tgCniEim87s+luh95oQJJPLkjYGqhkTkZmAG4AeeUdWVHod1EBGZCpwFtBORPODXqvq0t1EdZDRwFbDcuccD8AtVfdfDmOJ1Ap53Rkr5gFdVNSmHcDcCHYA3otcqAsDLqvq+tyHV6BbgJecP0A3AtR7HUyMRySQ6mvhGr2Opiap+ISKvAYuAELCYBJc8avbDzI0xxiQn6+IzxhiTlCxBGWOMSUqWoIwxxiQlS1DGGGOSkiUoY4wxSckSlGmSRKSDiLwsIhucMjyfi8jEozxXj4asIu+U4+nfUO/nvOftzjDnIz3uz5XVI0Tkd06V6+5x+3zkVCEw5ohYgjJNjvMQ4ZvAp6p6vKoOI/qwa1dvI6udiPhV9fr6rgAvUYf7fb+daPHPIzlnG+AUp5AxqvoL4EngjrhdXwB+dCTnNgYsQZmm6RygQlUfr9ygqptV9S8QmxPqWWceo8UicrazvYeI/EdEFjlftZZxEZFbRWSViCwTkWnOtizX+ZeJyCRn+5XOthUi8oDrHEUicp+IfAGcKiKzRWS467X7JTp/1VwR6eBsP8FZn+8cW1RDbD0kOv/R34g+XNlNRB4TkQXimgtLRG4lWlttlojMcrad57Q6F4nIP536ivEuA+IfzH0PmOw8CF1pOtEKCcYcEUtQpikaQPSCfCg3AajqyUQvnM9LdOK1XcA3nOKnVwCP1OG97gSGqOpA4AfOtruJVno+2dk+U0Q6Aw8QTZ6DgREicqmzfwtghaqOUtU5cedvAcxV1UHAp8D3ne0PAw+r6ggOXzuyD9E5xIY40zf8j6oOBwYCZ4rIQFV9xDnH2ap6toi0A34JnOt8FguAH9dw7tFE5/xymwy0JloGBwBV3QekiUjbw8RpzEEsQZkmT0QedVog851NY4h2O6Gqa4DNQG+ic0M9KSLLiVaTrst9oGVEy+h8l2j5F4jW93u0cgfnAj0CmO0U2gwBLxGdpwggTLTAbk0qgMpSTAuBHs7yqVRVvH75MPFtjpuM8VsisohomZoB1PwznuJs/8wpWTUF6F7Dfp2ITmUBxCpxXwv8BPhO3L67iLbSjKmzZl+LzzRJK4FJlSuqepPTKqickrymKVYgeu9kJzCI6B9vZfE7iMizROe52qaqFxKdUPAMolMk3C0iA5zzx9cQO9R7ApSpavgQrwW1qh5ZmCP/nS12xd6TaPIYoar7ROQ5oKYpu4XoPFm1dcuVxh1/BfAx0QkWfy4iGapa6ryW7uxvTJ1ZC8o0RTOBdBH5oWubewDApzh/4YtIb+A44EsgB9iuqhGiRW/d91EAUNVrVXWwql7oDDropqqziE4s2ArIAj4Abq48xhnB9gXRLrV2zv2ZK4FPjuFnnEtVEp5cx2NaEk1YBc69rAtcrxUC2a5zjxaRE534M53PKd5q4ETX+q3AQ05S+pBo0q4ctNIR2FTHOI0BLEGZJshpcVxKNCFsFJF5wPPAz51d/gb4na68V4BrVLXc2T5FROYS7fIrPvjs1fiBF53zLCZ6cd4P/BZo7QyGWEr03s524C5gFrCU6Jw//z6GH/N24MfOz9YJKKjtAFVd6sS5EngG+Mz18hPAeyIyS1XzgWuAqSKyjGjC6lvDKd8hWmEfERkD7HHNBPACVd18w4jeRwsddAZjDsOqmRvTCDnPLJWqqorIZOBKVZ3gQRxzgIudxHyofR4Gpqvqxw0XmWkK7B6UMY3TMOCvTvfZfuB7HsXx/4h2kR4yQREdoWjJyRwxa0EZY4xJSnYPyhhjTFKyBGWMMSYpWYIyxhiTlCxBGWOMSUqWoIwxxiSl/w/JmDAB0z07PAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "germany = suite.Copy(label='Germany')\n", "argentina = suite.Copy(label='Argentina')\n", "thinkplot.Pdf(germany)\n", "thinkplot.Pdf(argentina)\n", "thinkplot.decorate(xlabel='Goal-scoring rate (λ)',\n", " ylabel='PMF')\n", "pmf.Mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "According to this prior, the goal-scoring rates are always greater than zero, with the most likely value (a priori) near 0.5. Goal scoring rates greater than 5 are considered unlikely.\n", "\n", "### Step 3: Comparing posteriors\n", "\n", "The next step is to compute the posteriors for the two teams:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "posterior mean Germany 1.1506263832709267\n", "posterior mean Argentina 0.6693817986970533\n" ] } ], "source": [ "germany = suite.Copy(label='Germany')\n", "argentina = suite.Copy(label='Argentina')\n", "germany.Update(1)\n", "argentina.Update(0)\n", "\n", "print('posterior mean Germany', germany.Mean())\n", "print('posterior mean Argentina', argentina.Mean())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Update` invokes the likelihood function for each hypothetical value of $\\lambda$ and updates the distribution accordingly.\n", "\n", "Since both teams scored fewer goals than the prior mean (1.4), we expect both posterior means to be lower. \n", "\n", "Here are the posteriors:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8XNWZ+P/Poxn1ZhV3uci94IIbNqYYm2KqIUBiQxaHDSEFQkny3UD2y27CJvmGLL84ISHZhdBCwEBMDA41gG3AgI1777JlyZYtWXJRsaQp5/fHvRrdGatrRhrJz/v1Err3zp17jwdbj845z32OGGNQSimlok1MZzdAKaWUaogGKKWUUlFJA5RSSqmopAFKKaVUVNIApZRSKippgFJKKRWVNEAppZSKShqglFJKRSUNUEoppaKSu7MbEC7Z2dlm8ODBnd0MpZRSzVi/fv1xY0zP5s7rNgFq8ODBrFu3rrOboZRSqhkikt+S83SITymlVFTSAKWUUioqaYBSSikVlbrNHJRSSrWXx+OhsLCQ6urqzm5Kt5CQkEBOTg6xsbFter8GKKWUshUWFpKamsrgwYMRkc5uTpdmjKG0tJTCwkJyc3PbdA0d4lNKKVt1dTVZWVkanMJARMjKympXb1QDVAv4/brqsFLnCg1O4dPez1IDVDMWf3GI7z2/gdfXFnZ2U5RS6pyiAaoJJadr+Gh7MV6/4b0tRyk5XdPZTVJKdXPHjh3jtttuY8iQIUyePJkZM2awdOnSzm5Wp9AA1YS1B8oC28bAyl3FndgapVR3Z4zhxhtv5JJLLiEvL4/169fzyiuvUFjYshEcn88X4RZ2LM3ia8K6vBNB+6t2H2fepP7EuTWuK9Xd3Xz//0Ts2q//7jsNHl++fDlxcXF85zv1rw8aNIjvf//7+Hw+HnroIVauXElNTQ333HMP3/72t1m5ciU/+9nP6Nu3L5s2beKdd95h7ty5XHTRRaxevZoJEyZw55138p//+Z8UFxfz0ksvMW3aNL788kseeOABzpw5Q2JiIs899xwjR47k+eefZ9myZVRVVbF//35uuukmfv3rX/PMM8+wbds2Fi1aBMDTTz/Nzp07+c1vfhOxz0l/0jai+HQ1h0qrgo5V1vhYm1fWyDuUUqp9tm/fzqRJkxp87ZlnniE9PZ21a9eydu1ann76aQ4cOADAl19+yS9+8Qt27NgBwL59+7j//vvZsmULu3bt4uWXX2bVqlU8/vjj/PKXvwRg1KhRfPLJJ2zcuJFHH32Un/zkJ4F7bdq0iVdffZWtW7fy6quvUlBQwPz581m2bBkejweA5557jjvvvDOSH4f2oBqz7kB978ntErw+K5Nvxc5iZo7I7qxmKaXOIffccw+rVq0iLi6OQYMGsWXLFpYsWQLAqVOn2Lt3L3FxcUybNi3oWaPc3FzGjRsHwNixY5kzZw4iwrhx4zh48GDg/QsXLmTv3r2ISCDwAMyZM4f09HQAxowZQ35+PgMGDGD27Nm89dZbjB49Go/HE7hHpGiAaoRzeO/mqTm8vrYQr89wsKSKAyWV5PZM7sTWKaUirbFhuEgaO3Ysr7/+emD/ySef5Pjx40yZMoWBAwfy+9//nquuuiroPStXriQ5OfjnUXx8fGA7JiYmsB8TE4PX6wXgkUce4bLLLmPp0qUcPHiQWbNmNfh+l8sVeM9dd93FL3/5S0aNGhXx3hPoEF+Djp2qH95zxwgzh2cxdUhm4PUVOzRZQikVfrNnz6a6upo//elPgWNVVdbPoquuuoo//elPgZ7Onj17qKysbPO9Tp06Rf/+/QF4/vnnW/SeCy64gIKCAl5++WUWLFjQ5nu3lAaoBjiH984bkE5SvJvLRtevrfVlXhkV1d7OaJpSqhsTEd544w0+/vhjcnNzmTZtGgsXLuSxxx7jrrvuYsyYMUyaNInzzjuPb3/724GeTVv827/9Gw8//DAzZ85sVfbfV7/6VWbOnElGRkab791SYkz3qJIwZcoUE64FCx9duiPQg7prVi7Th2VhjOHnb+4k/7h1/O7LhjBtaGZTl1FKdTE7d+5k9OjRnd2MqHbdddfx4IMPMmfOnBad39BnKiLrjTFTmnuv9qBCBA3vuYQJA3sA1m824wakB87be6y8U9qnlFKd4eTJk4wYMYLExMQWB6f20iSJEM7U8pF9U0mMcwX2h/dOCWzvP9b2sV+llOpqevTowZ49ezr0ntqDCuFzFIZNjguO30N6pVBX+7CgrIrq2u711LZSSkUTDVAhnAHK5QquxJsY56J/RiJglT7KK9FelFJKRYoGqBB+R9KIq4FS8cMcw3z7jlV0SJuUUupcpAEqhN9fvx3TwKczLGgeSgOUUkpFigaoEF5HhIpppge1v7hCFzNUSoXd0qVLERF27doV8XsdPHiQl19+ObC/bt067rvvvojftyU0QIVwxhtXzNkBKisljh5JsQBUe/wcPnGmo5qmlDpHLF68mIsuuohXXnnlrNfCvaRGaICaMmUKTzzxRFjv0VaaZh7C2SNqKECJCEN7p7Derjax71gFA7KSOqx9SqmOcdefw/Pgf0P+fFfjz6hWVFTw2WefsWLFCm644QZ++tOfnrWkxo4dO/iv//ovXnrpJQYMGEB2djaTJ0/mRz/6Efv37+eee+6hpKSEpKQknn76aUaNGsU3vvEN0tLSWLduHUePHuXXv/41t9xyCw899BA7d+5k4sSJLFy4kPPPP5/HH3+ct956i5/+9KccOnSIvLw8Dh06xAMPPBDoXd14440UFBRQXV3N/fffz9133x32z0kDVAifaTpAgTXMVxeg9hdXcNmYXh3SNqVU9/fGG28wd+5cRowYQWZmJhs2bACsJTW2bdtGbm4u69at4/XXX2fjxo14vV4mTZrE5MmTAbj77rv5n//5H4YPH86aNWv43ve+x/LlywEoKipi1apV7Nq1ixtuuIFbbrmFX/3qV4GABFbxWaddu3axYsUKysvLGTlyJN/97neJjY3l2WefJTMzkzNnzjB16lRuvvlmsrKywvpZaIAK4exBNTQHBcHzUHuPaqKEUip8Fi9ezAMPPADA/PnzWbx4Mddee23QkhqrVq1i3rx5JCZaj71cf/31gNX7+vzzz7n11lsD16upqQls33jjjcTExDBmzBiOHTvWovZce+21xMfHEx8fT69evTh27Bg5OTk88cQTgaXoCwoK2Lt3rwaoSPMGBaiGzxmQmUisS/D4DKUVtZysrKVHclwHtVAp1RGaGoaLlNLSUpYvX862bdsQEXw+HyLCNddcE7SkRmM1VP1+Pz169GDTpk0Nvu5cRqOldVgbWnpj5cqVfPjhh3zxxRckJSUxa9YsqqurW3S91tAkiRCmmSQJALcrJmg9qL2abq6UCoMlS5Zwxx13kJ+fz8GDBykoKCA3N5dVq1YFnXfRRRfxj3/8g+rqaioqKnj77bcBSEtLIzc3l7/97W+AFYQ2b97c5D1TU1MpL29dbdFTp06RkZFBUlISu3btYvXq1a16f0tpgArhrCQR01gXChjeJzWwnVesFSWUUu23ePFibrrppqBjN998c1CWHcDUqVO54YYbmDBhAl/5yleYMmVKYAXcl156iWeeeYYJEyYwduxY3nzzzSbvOX78eNxuNxMmTGDRokUtaufcuXPxer2MHz+eRx55hOnTp7fiT9lyutxGiFdXF/DBNmts9tZpOVw1vk+D523KP8kfPtgHwPA+Kfz4ulHtvrdSqnN1peU2KioqSElJoaqqiksuuYSnnnqKSZMmdXazztKe5TZ0DiqEvwVZfACDsutTyw+VVmGMQRpJqlBKqXC7++672bFjB9XV1SxcuDAqg1N7aYAK0dIhvh5JsaQmuCmv9lLj8XPsVA19eiR0RBOVUuqsYb/uSOegQgQ9qNtEj0hEzupFKaW6vu4y7REN2vtZaoAK4XxQt6FisU6Dsusz+Q4e10QJpbq6hIQESktLNUiFgTGG0tJSEhLaPrKkQ3whfM2UOnLSHpRS3UtOTg6FhYWUlJR0dlO6hYSEBHJyctr8fg1QIZzLbTQ1xAchAeq4Jkoo1dXFxsYGqjWozhfRIT4RmSsiu0Vkn4g81MDr8SLyqv36GhEZbB+PFZEXRGSriOwUkYcj2U6n4CG+poNNZnIcyfEuAKpqfZSU1zR5vlJKqZaLWIASERfwJHA1MAZYICJjQk77JnDCGDMMWAQ8Zh+/FYg3xowDJgPfrgtekdZcNXMnEWGwYx4q/7gO8ymlVLhEsgc1DdhnjMkzxtQCrwDzQs6ZB7xgby8B5og1RmaAZBFxA4lALXA6gm0N8LUwi6/OQMcwnwYopZQKn0gGqP5AgWO/0D7W4DnGGC9wCsjCClaVQBFwCHjcGFMWegMRuVtE1onIunBNavpbMcQHmiihlFKREskA1dBP99DczcbOmQb4gH5ALvBDERly1onGPGWMmWKMmdKzZ8/2thdoXRYfwKCs4FRzTU9VSqnwiGSAKgQGOPZzgCONnWMP56UDZcBtwHvGGI8xphj4DOiQ2vfOANWShLzs1DiS4uxEiRofpRW1kWqaUkqdUyIZoNYCw0UkV0TigPnAspBzlgEL7e1bgOXG6oIcAmaLJRmYDuyKYFsDHPEJdwt6UCKi81BKKRUBEQtQ9pzSvcD7wE7gNWPMdhF5VERusE97BsgSkX3AD4C6VPQngRRgG1age84YsyVSbXVqabFYJ52HUkqp8Ivog7rGmHeAd0KO/YdjuxorpTz0fRUNHe8IwUN8LQxQWVrySCmlwk1r8YVobZIEBPeg8u2KEkoppdpHA1SItgSoXmnxJMRaH2VFtZcTlZ6ItE0ppc4lGqBCBM1BtXCIT0QYkFXfiyrQeSillGo3DVAhnMVim1tuw2mgI0AdKtMApZRS7aUBKoTXEaFiWlGZ3BmgNNVcKaXaTwNUCGd+Q0vnoEBTzZVSKtw0QIVoS5IEQJ/0BNwu6/yyiloqqr1hb5tSSp1LNECFCFoPqhVDfG5XDDmZiYF97UUppVT7aIAK0Zr1oEINzNRhPqWUChcNUCGcQ3ytyeIDglLND2mihFJKtYsGqBDOYrEtfQ6qTlBFiVIteaSUUu2hAcrBGNPmJAmA/pmJgSU6ik/XUOPxhbN5Sil1TtEA5eBMMRdpebHYOvFuF33SEwLXKig7E87mKaXUOUUDlIOvDWWOQgUXjtVhPqWUaisNUA7+oASJtgWogUE1+bQHpZRSbaUByqE98091gjL5NNVcKaXaTAOUQ/BDum27xiBHgDpy4gxen7+Js5VSSjVGA5SDs5J5W3tQSfFuslPjAPD6DYdP6DCfUkq1hQYoh3AM8UHI0hs6zKeUUm2iAcqhrXX4Qg3OTg5sH9SKEkop1SYaoBzaU4fPSVPNlVKq/TRAOfjCkGYOMMjRgyos00QJpZRqCw1QDv4wPKgLkJLgJivFTpTwGY6crG5325RS6lyjAcohXEkSoMN8SinVXhqgHNqz1EYo5zBfviZKKKVUq2mAcmjPUhuhgntQGqCUUqq1NEA5hCtJAkJq8pVVaaKEUkq1kgYoh6AkiXYGqLTEWDKSYwErUaJIEyWUUqpVNEA5BPWg2jnEB6EP7GqihFJKtYYGKIfgANX+6+k8lFJKtZ0GKAfnirpuV/s/Gmcmn9bkU0qp1tEA5RDJHlRBaVXQ9ZVSSjVNA5RDuOegnIkSHp+h6KQuvaGUUi2lAcohnJUk6gzSyuZKKdUmGqAcwplmXsc5zHewRDP5lFKqpTRAOfjDtB6UU25PTTVXSqm2iGiAEpG5IrJbRPaJyEMNvB4vIq/ar68RkcGO18aLyBcisl1EtopIQiTbCuGtxVfH+SxUYekZar1aUUIppVoiYgFKRFzAk8DVwBhggYiMCTntm8AJY8wwYBHwmP1eN/BX4DvGmLHALMATqbbWicQcVEqCm95p8QB4/YbCMp2HUkqplohkD2oasM8Yk2eMqQVeAeaFnDMPeMHeXgLMEREBrgS2GGM2AxhjSo0xvgi2FQhvsVin3F71vagDOg+llFItEskA1R8ocOwX2scaPMcY4wVOAVnACMCIyPsiskFE/q2hG4jI3SKyTkTWlZSUtLvB4SwW6+Sch9IApZRSLRPJANXQT/jQJ1UbO8cNXATcbn+/SUTmnHWiMU8ZY6YYY6b07Nmzve3FH4EhPggOUHnFGqCUUqolIhmgCoEBjv0c4Ehj59jzTulAmX38Y2PMcWNMFfAOMCmCbQXAF6Yl30PlZCbhtgNe8ekaKqq9Ybu2Ukp1V5EMUGuB4SKSKyJxwHxgWcg5y4CF9vYtwHJjjAHeB8aLSJIduC4FdkSwrUDkhvji3DHkZCUG9jXdXCmlmhexAGXPKd2LFWx2Aq8ZY7aLyKMicoN92jNAlojsA34APGS/9wTwG6wgtwnYYIx5O1JtrROpIT7QeSillGotdyQvbox5B2t4znnsPxzb1cCtjbz3r1ip5h3GZ8JbLNYpt2cyK7ASOQ7oPJRSSjVLK0k4+B3P0Ea6B2WMVjZXSqmmaIByiMSDunX6pCeQFOcCoLzaS2lFbVivr5RS3Y0GKAdfBGrx1RERBus8lFJKtZgGKIdIJkmAPg+llFKtoQHKIVJp5nU0k08ppVpOA5SDP0IP6tYZ0it46Q2PTyubK6VUYzRAOUQySQKsJeB71VU29xkO6Qq7SinVKA1QDs4AFYEOFABDe6cEtvcdq4jMTZRSqhvQAOXgXG7DHYEeFMCgzASqazwYYzRAKaVUEyJaSaKr8UcwSeJU+Rne/ngrSz/ZwV5PKjECBwuOke09wSVTRtC3Z3pY76eUUl2dBiiHSDwH5fX6ePntL3nnk214vD4MECPJ+ImhtLyGl97bxNIPN3HvbZdx0eRhYbmnUkp1BzrE5xCJJIln//45by7fjMdrLQgsQIq7PnuvGjcer49Ff/mQv7z5BT7N7FNKKUB7UEGChvjCEJ/WbDnA+59tD+wP7p/NzVeeT6k/gb+tKaC8shrPSaDcKiL75vLNHCgs5cd3XUVCfGz7G6CUUl2Y9qAcgpIkXO37aI6fqOCPi1cG9qePz+Xx/3MzF04cyoi+acTFusnqkcK484YxecygwHlb9hTyxF+XazFZpdQ5TwOUgy9MPSi/38/vXvyIiqoaALJ6JPPdBbMQe15rcM+kwBBiSYWH799xObdcWb9g8JotB3jl3XVtb4BSSnUDTQYoEXnesb2wiVO7heAA1fYI9ebyzezYXwRYc04P3nE5KUnxgdfj3S4GZiUF9vOKK1lw7TSuu3R84NiS99fz2cb9bW6DUkp1dc31oCY4tu+PZEOiQVCpozZ2oTweH28u3xzYv3XuFEYP7XvWecOcD+wWW89D3TFvOhNG5gSO//6vyzlQeLxN7VBKqa6uuQB1Tk2EhCOLb82WA5RXVgOQnZESNHTn5AxQe49aAcrliuEH37gi8EyUx+tj0QsfUuvxtqktSinVlTUXoHJE5AkR+b1jO/DVEQ3sSP4wDPE5s/YunzEaVyPJFs4A5Swcm5IUz0Pfmkt8nJXFd7j4JK/pfJRS6hzUXID6P8B6YJ1j2/nVrfjaOcRXcPREYO4pRoQ500c1em56UnDh2IOO5Tdyemdwxw3TA/tvLt/MvvziVrdHKaW6siafgzLGvNBRDYkGwetBtf79H36+M7A9ddxgMtOTmzgbRvRJpfi0lem3q6ic4X1SA69dddEYPt+0n+37juA3hj8sXsl///BmYmNdrW+YUkp1Qc1l8S1r6qujGtlRnM9BtXY9qFqPlxVf7g7sX3HhmGbfM6pffUDaXVQe9JqI8N35lxIXa/0OUVBUxpIPNrSqTUop1ZU1V0liBlAALAbWYGVNd1vtWVH3i015VJ6xekO9MlOZOCqnmXfAqL71AWrfsQpqvX7i3PW/M/Ttmc7t103juaWfA7D0w41cMmU4/Xv1aFXblFKqK2puIKsP8BPgPOB3wBXAcWPMx8aYjyPduI7mb0cW3/uf7QhsX37h6MBDuU3pkRxH7/T6eai84rOX37j20nGMGNwbAJ/Pz/N2sFJKqe6uyQBljPEZY94zxiwEpgP7gJUi8v0OaV0H87VxyfeSsnJ2HzgKQExMDLMvaDw5ItRIRy9qV8gwH1hDfd+65aJA13XDjkOs357f4usrpVRX1WwqgIjEi8hXgL8C9wBPAH+PdMM6Q1uH+LbsKQxsjxvej4y0pCbODjaqX1pgO3Qeqs6QAT2Z7cgIfH7p53jt6uhKKdVdNZck8QLwOTAJ+JkxZqox5r+MMYc7pHUdrK3VzDfvrv84Jowa0Kp7OuehDhRXUtNI4Ln9ugtISogD4EjJKd7+ZFur7qOUUl1Ncz2ofwFGYJU5+kJETttf5SJyOvLN6zjGmOAsvhZGKGMMW/fUB6jxI/q36r5pibH07ZEAgNdv2H+sssHz0lMT+ercKYH9195bx8nyqlbdSymlupLm5qBijDGpjq80+yvVGJPW1Hu7GmdwihFalOQAcKiojNMVZwCrCsTg/lmtvrcz3XzXkcbj/tUXjw1k8FXXeHjt3W73rLRSSgU0N8SXICIPiMgfRORuEem2Cxy2ef7JMbw3bkROiwOb06i+9bG+oUSJOm63iztunBHY/+CLnRwpPtnq+ymlVFfQ3BDfC8AUYCtwDfD/RbxFnaStKebOBInWDu/VGdHXUZevpJLq2sYTICaPGcgYuzq63+/npbe+bNM9lVIq2jUXoMYYY75ujPlf4Bbg4g5oU6doS4q51+tj+76iwP74kc0/nNuQ1IRYcjITAWuoce+xs5+HqiMi3DGvvk7f6s157Dl4rE33VUqpaNZcgPLUbRhjuvWaD20Z4tubX0xNrfUR9cpMpU9226flRjvSzXccbjr/ZPig3syYODSw/+Ky1bpEvFKq22l2wUJn5h4wvrtm8bVliG+z8/mnNg7v1RnTvz5AbSs81ez5t183jRi7ou2O/UWs33GoXfdXSqlo01wWnyskc8/dXbP4nEN8LZ2CCkovb+PwXp2RfVOJdVk3LjpZzfHymibP79sznSsvHB3Yf+kfa7QXpZTqVtqwqETLichcEdktIvtE5KEGXo8XkVft19eIyOCQ1weKSIWI/CiS7QTw++u3W7JYYXWNhz0H69doGje8fT2oOHdMUNmj7YXNd1C/OndKYGHDQ0VlrFq/r11tUEqpaBKxACUiLuBJ4GpgDLBARELXoPgmcMIYMwxYBDwW8voi4N1ItdHJOQflbkEXavu+I/jtqDaoXxbpqYntbsN5OemB7a0tGOZLT03k+lnjAvuvvLtWSyAppbqNSPagpgH7jDF5xpha4BVgXsg587BS2QGWAHPEfpBIRG4E8oDtdAC/aV2SxJ788PWe6jgD1M4jp/H6/E2cbblh9gSSE62K6EePn2b5mt3NvEMppbqGSAao/lhrSdUptI81eI6dJXgKyBKRZODHwM+auoH98PA6EVlXUlLSrsb6W7nc+8HC44HtIQOy23XvOr3T48lOtert1Xj87Gsi3bxOcmI8N10+MbD/2nvrqPV064RLpdQ5IpIBqqGf8qGz+I2d8zNgkTGmyZ/QxpinjDFTjDFTevbs2cZmWrw+Z5JE8wHqwOH6AJWbE54AJSJBvahtLZiHArjmkvMCFdRPnK7i3U87pNOplFIRFckAVQg4S3vnAEcaO8cuo5QOlAEXAL8WkYPAA8BPROTeCLY1ZIiv6XNPV5yh9KRV1DXW7QrrCrfjggJU8/NQAPFxsdxy5eTA/t8/2BBY3VcppbqqSAaotcBwEckVkThgPrAs5JxlwEJ7+xZgubFcbIwZbIwZDPwW+KUx5g8RbGtQsVh3MxHqwOHSwPbAvpm4XOH7GEf1Sw0kaRSWneFEZW2L3nf5jFH0yrSyACuqavjHyi1ha5NSSnWGiAUoe07pXuB9YCfwmjFmu4g8KiI32Kc9gzXntA/4AXBWKnpH8bViLai8gvr5rnDNP9WJj3UxrE99bb6W9qLcbhfzr5ka2F+2fEugyrpSSnVFEX0OyhjzjjFmhDFmqDHmF/ax/zDGLLO3q40xtxpjhhljphlj8hq4xk+NMY9Hsp0QslhhMxHK2YPK7R/eAAXBw3xbC1oWoAAunjyMnN4ZANTUelj64aawt00ppTpKRANUV+L1t7xYrDODL1wJEk7jBtQHqO2HT1PrbT7dHCAmJiaoF/Xup9soPdl8JqBSSkUjDVC2lj4HVV3jCazBJMCgfplhb0vfHgn0TrOebarx+NnZTPFYp+kTchkywMpo9Hh9LPnnhrC3TymlOoIGKFtLi8XmHykN5Mr3750RKDUUTiLC+YMzAvsb8k+06r23XzctsP/hF7soKmn5MKFSSkULDVC2lhaLPVBYP/80OKf1y7u31KTB9anrm/JPBiVxNGfCyJygRQ1ffXdd2NunlFKRpgHK5iwW21QPKq/QkcGX076Hg5uS2zOZHklW76yyxsfeo40vBR/K6kVdENhftX4vBx0PFiulVFegAcrma+EQ30FHBt/g/pHrQYkIEwfV96I25p9s1ftHDenDlLGDAKs0x8tvrQ1n85RSKuI0QNmCh/gaDlBer4/8orLA/pAIZPA5TXLMQ23MP9nq9Z5uu25aoJbU+h357Nxf1OT5SikVTTRA2VqSJHG4+GRgOYusHsmkJidEtE0j+qSQFOcCoKyilvzjVa16/6B+WVw8ZXhg/0Vd1FAp1YVogLL5WvCg7gFnBfMIzj/VcbtiGD+w/pmo1g7zAXzt6qmBUky7DxzVpeGVUl2GBihb0HIbjQzxOTP4BkVw/snp/EGOdPODLU83r9MnO40rZtQvDf/Xf6wJLLSolFLRTAOUrSVJEvlFzhJHHROgzstJI9ZltafoZDWHy1pfX+/WuZMDz2sVFJWx8ss9YW2jUkpFggYoW/AQX8Pn1FWQABjQN/wVJBoSH+tiwsD6bL7V+0ubOLthPVKTmDd7QmB/8Ttrqan1hKV9SikVKRqgbM7nYBsa4qup9QTWgIoRobe9tEVHmD6svre2el9pmxId5s2eQHpqIgBlpyp5++NtYWufUkpFggYoW3O1+I4er6+H1ysrFbfb1SHtAmuYLzneut+JSg97jra+AGxCfCxfmzslsP/3DzfqchxKqagXp7BgAAAgAElEQVSmAcrm8zU9B3XYMbzXL4wr6LaE2xXDlCH1Q4pr2jDMBzBn+ij69bSyAs9U12ohWaVUVNMAZWvuQd0jxfUFV/v17NgABTDDMcy3Lu8EHl/rM/Hcbhe3X19fAundT7cHBV6llIomGqBszWXxOSuC9+2ZftbrkTa0VzJZKXEAVNX62NaKhQydLhify6ghfQCrkOyLb64OWxuVUiqcNEDZnHkHzQWofr06PkCJSFCyxBf72jbMJyLceeOFgf212w6yZXdhu9unlFLhpgHK5nU8vNrQY1DOFPPO6EEBXDC0fh5qS8Epqmq8bbrOsEG9mDVtZGD/uaWf68O7SqmoowHKFvQcVMgcVHllNeWV1QDEul1kZ6R0aNvq9MtIZGBWEgBen+HLvLJm3tG426+bRlysG4BDRWV8tHpXWNqolFLhogHK5uxAuF3BASp0/kkaKYXUES4cXj/M98mu420u/pqZnsxNl08M7L/89loqz9S0u31KKRUuGqBs/iay+ILmnzppeK/OjOFZgdJHh0qrOFBS2eZrzZs9gaweyQCcrjijK+8qpaKKBihbUwHqSCc+AxUqOd7NVMczUZ/savtKufFxsdwxb0Zg/91PtpF/pG3JF0opFW4aoGzeJh7UPex8BqqTAxTArNH1S32s2V/a5mQJgJnnD+W84f0AK0g//bdVumaUUioqaICyBZc6Cn6ts5+BCpXbM5mcTKuunsdn2pxyDlba+V23XEyM/YfemVfEJ+v2hqWdSinVHhqgbM5isW5HhDLGdPozUKFEhEtH1fei2pMsATCgTwbXzxoX2P/Lm6s1YUIp1ek0QNmcaebOKagTp6sCS1MkJ8ZHfJn3lpo+LIv4WOt/3+ETZ9hf3PZkCYBbr5pMRpqVwn6yvIrFb69tdxuVUqo9NEDZ/I2UOgp9QLczU8ydEuNcXOBIlli+vbh910uI4xs31VeYeO/Tbew+cLRd11RKqfbQAGXzNbLke1CR2CgY3nOaNaZXYHvdgTKOl7dvWG7m+UM5f/QAAAzwp1c+xuv1teuaSinVVhqgbMEr6tYHqGhLkHAamJXE6H7Wwol+Ax9sO9au64kId3/1kvrl4Y+eYOlHm9rdTqWUagsNULaWDPFFQ4p5qKvG9wlsf7r7OBXVbU85B+iVmcqCa6YG9v/2/npdkkMp1Sk0QNmC14OqP+7sQfWPwgA1tn9aIOW81utn5c72zUUBXHvpeQwdYGUJ+nx+nnx5pRaTVUp1OA1QNufP37oelN/v52hp/VLv0TbEB9aw3FXj6ntRH20vptbbvmASExPD9xZcGng2aveBo/xj5dZ2XVMppVpLA5StoQULS09W4rNXrk1LSSQhPrZT2tacqUMyyEi22lZe7eWLve0vVzS4fzY3X3l+YP+lt9ZwqKjt1dOVUqq1NEDZGlryvbisPHCsd1Zqh7eppdyuGK509KLe33oUbxuWhA91yxWTyM3JBqyhvt+/tEKz+pRSHUYDlK2hJIni0voA1TMzegMUwMUjskmKcwFQfLqmXeWP6rjdLu77+mzcbuu6eQUlvP7BxnZfVymlWiKiAUpE5orIbhHZJyIPNfB6vIi8ar++RkQG28evEJH1IrLV/j47ku2EhtPMg3pQUR6gEuJczHVk9C3bcKTdc1EAA/tmctu10wL7S95fz648fYBXKRV5EQtQIuICngSuBsYAC0RkTMhp3wROGGOGAYuAx+zjx4HrjTHjgIXAi5FqZx1/Aw/qOgNUr6y0SDeh3eaM7UVqgrVK7olKD5/sKgnLda+fNY5RQ6zg5zeG3/7lI63Vp5SKuEj2oKYB+4wxecaYWuAVYF7IOfOAF+ztJcAcERFjzEZjzBH7+HYgQUTiI9jWkCXfre8lZV1niA8gPtbFdef3Dey/vamIGk/754xiYmK4/1/mkJQQB0DJiXL++PJKXZZDKRVRkQxQ/YECx36hfazBc4wxXuAUkBVyzs3ARmPMWb+yi8jdIrJORNaVlLSvtxCUxec6ew6qVxQnSThdMqonmSlWICmv9vJhO2v01emVmcr3FswK7K/ecoAPPt8ZlmsrpVRDIhmgGqqqGvord5PniMhYrGG/bzd0A2PMU8aYKcaYKT179mzolBZzdgZcIvh8fkpPVgSO9cxIadf1O0qsK4YbJvUL7L+/5Wi7q0vUmTFxCFfNHBvYf/bvn3HwcNtX9FVKqaZEMkAVAgMc+znAkcbOERE3kA6U2fs5wFLgDmPM/gi2EwgpFhsjlJ6qDMxL9UhNIi7WHekmhM2MYVn0TrdGRKtqffx9XWHYrv2Nm2YwsK9VRd3j9fHYn9+nvLI6bNdXSqk6kQxQa4HhIpIrInHAfGBZyDnLsJIgAG4BlhtjjIj0AN4GHjbGfBbBNgb4fM71oCRk/qlr9J7quGKEW6fV/27w6e7jHChp33pRdeJi3fzwzisCBWWLy8pZ9MKHWgpJKRV2EQtQ9pzSvcD7wE7gNWPMdhF5VERusE97BsgSkX3AD4C6VPR7gWHAIyKyyf7qRQSF9qCC55+iP4Mv1ISB6YwbYJVmMgYWf3EobEkNOb0zuO/rlwX2N+8u1AUOlVJhF9HnoIwx7xhjRhhjhhpjfmEf+w9jzDJ7u9oYc6sxZpgxZpoxJs8+/nNjTLIxZqLjKzyz/Q23M2gOKkbgWFl9Db5eXWT+yUlEWDBjAG47JTGvuJJVe8I3XzR9whBuvmJSYP/vH27k800RH4lVSp1DtJIEZ9fhs4b46hMkumIPCqBXWkLQchyvrz1MZU14EiYA5l8zJbDAIcATLy5nz8H2rUmllFJ1NEDR8FIbxY4q5l3hGajGXDOxTyDtvKLay9/WhC9hIiYmhgfuuDxQ5d3j9fH/nn6Po8dPN/NOpZRqngYogpfaqCtzFNyD6roBKt7tYv70+l7Oqj3H2ZwfvgUIU5Li+fdvX0NKkpU1eLriDL/833eoqNJKE0qp9tEAxdlDfF31GajGTBqcwdQhGYH9F1YdpLzaE7br9+2ZzsPfujpQVPZw8Uke+/N71HrCN5yolDr3aIDi7Dp8zmegMtK61jNQjbn9wkGkJ1mp4afPeHnps/Bl9QGMGtKH799en9m3Y38R//3sP3V5DqVUm2mA4uylNrrL/JNTSoKbhRcPCuyvO3CCtXknwnqPiyYN4455MwL7G3Yc4rcvLtdnpJRSbaIBCvD6nQ/pBs8/dZcABTB+QA8uHpkd2P/LqoMcPRneKhDzZk/glivr08+/2LSfP77ysRaWVUq1mgYoQob4YiToGahoXweqtb56wQCyU62svmqPnyc/3BeWiudO86+ZyjWXnBfYX7FmN08uXqk9KaVUq2iAIjiLzxUj3bYHBZAY5+J7lw8j1q7YXnSymuc+ORjWHo6I8K9fmcmsaSMDx1as2c1vX1yuc1JKqRbTAEXoc1DBc1BdOcW8MQOzkviXi4Lnoz7YFt4HbEWEexZcymUX1Aepzzbs4zcvfKhBSinVIhqgaCBJoostVNgWFw7PZtbo+iVK/vZlIZvC+HwUWA/y3rNgVtASHWu2HODn//uOrsirlGqWBihCkiQwlJ7oPs9ANWX+9AEM6ZUMWAVl/3f5fvYdq2jmXa0jInzr1ou4ftb4wLGtew7z7797k+MnwnsvpVT3ogEKguZfaj3ewIqJ3eUZqMa4XTHce8UweqZaVSA8PsMT/9zLkRNnwnofEWHhjTNYcO20wLGCojIeXrRUFzxUSjVKAxTBlSSqHRUWuuvwnlNaYiwPXj2c1AQrEFfV+Fj03h5KK8I7BCci3HLlJO77+mxiYqy/dmWnKnl40Rus2rAvrPdSSnUPGqAICVA1tYHtcyFAgVX1/P6rhhMfa/11OFHp4ddv7ab4dPhXyr106gge+c41JCZYqe61Hi+LXviQF974Ap9P09CVUvU0QAGO+ETVmfoAVVel+1wwuGcy35szNLB+VGlFLb9+a3fYH+QFGD8yh//34E1Bn++yFZt59E9vUXYqPCv/KqW6Pg1QBPegnNllfbroOlBtNTYnnXuvrH9G6mSVh8fe2kVhWVXY7zWgTwaP/fArTBlbn+6+be8RHvzVa3y59WDY76eU6no0QBGcZl5ZVd9j6JN9bgUogPNy0rnPMdxXXu3lV//YxZaC8KagAyQnxvPQt+by1blTsJfhoqKqhsf+/B5PvfYp1TXhq7iulOp6NEDhfFDXUOlYx6j3ORigAEb3S+PBuSNIjLOWz6j2+Pn9P/fxwbZjYa+pJyJ87eop/PTe68lMTw4cf/+z7Tz4q9fYvDt8CywqpboWDVDUD/F5fX48dpWDuFg3GWlJndmsTjWsdwo/vm4kWfZqvMbAq6sLeOHTfGoiUAnivOH9+c2Pb2X6+NzAseKych7941v8/qUVlFeGfy5MKRXdNEBRP8RXU+tF7Keg+mSnISJNva3by8lM4t/njWZor/qezao9x/n5GzsjMi+VmpzAj/71Su697TKSE+MDx1d+uZt7f76Y9z7drpl+Sp1DNEBRP8RXU+MNzIWci/NPDUlLjOWH14xk+rDMwLGik9X8/M2dLN9RHJEhv8suGMnvfvI1ZkwcGjheUVXD00s+5Uf/vYQtOuyn1DlBAxT1Q3zVtfWT8n2yz50U8+bEuWP45qW5LLx4UCDDz+szvPz5If777d0UnQxv5Qmwqnj86M4r+PFdc+nleB7tUFEZP/vjW/z0yX+w52B4C9wqpaKLBijql9uo8QQP8al6IsLFI3vyyI1jyMlMDBzfc7SCn/19B8s2HKHWG/7ht2njBvO7n3yNBddOCyo7tXXPYR5etJRfPf2eBiqluikNUNQvWFhT4w0c63MOPaTbGv0yEvnJDaO5ekIf7Gd68foNyzYc4f8u2cYXe0vDPuwXF+vmlisn8ft/n8/sC0YR45gbXLvtIA8vWsojT7zJhh2HdOVepbqR7lsJtRV8jiSJJLsH1fsce0i3NeLcMdw8NYdpQzL5y6p8DpRY1R/KKmp55uMDfLDtGPMm92P8gPSwJppkZ6Rwz22zuPHyibz67jo+c9Tw27G/iB37i+jXM525F5/HrGkjghItlFJdj3SX3zinTJli1q1b16b3vr/lKK+uPsSGnQWkmyp6yRleefwuXC7tYDbH7zd8vKuEZRuOUF7tDXotJzORayb0ZXJuBq6Y8GdE5h8p442PNrJq/b5AL7hOfFwsM88fymUXjGT0kD7nfEamUtFERNYbY6Y0e54GKHhncxEvfXqA7fuL6GGqGJvl4slHbgtzC7u3M7U+3t9ylPe3HsXjC/47lZUSx6Wje3LRiGzSEmPDfu/isnL+sWIzy9fsbrD6RO+sNC6ZOpyZ5w9jQJ+MsN9fKdU6GqBa4e1NRTy7fA/7C46TYSqZPTKLR757bZhbeG44WVnL+1uP8fGukrOSJtwxwvmDezBjeBZj+qXhDnMP9Ux1LR+v3ct7q7ZRcPREg+fk9M5g+sQhTB07iKEDe2rPSqlOoAGqFZZtOMJT/9xF4bETZJhKbpuZy7duvTjMLTy3lFd7+Gh7MSt2FFNZc3bliZQEN1NyM5icm8Hw3ilhDVbGGPYcPMaKL3fz2Yb9VFXXNnhej9QkJo0ZyIRROYwb3p/01MQGz1NKhZcGqFZYuu4wf3xnGyUnKsg0lTw4byLXXza++TeqZtV6/aw7UMbHO0vYX9zwUhpJcS7GDUhn3IB0RvdLIz0pfMOAtR4v67bn8/nGPNZvz6fW42303IF9Mxk3oj+jhvRhVG6foNqASqnwaWmA0iw+rDTzmlr7B5cx52yR2EiIc8dw4fBsLhyeTWFZFZ/vLWVtXhknKuvniqpqfazZX8aa/WUA9MtIYHS/NIb2TmFYrxQy7XqAbbp/rJsLJw7lwolDqa7xsGHnIdZuPcjGnQVn1fc7VFTGoaIy3v54KwC9MlMZOrAXwwb2ZNjAngzun01KkmYGKtVRNEBhpZnX2FUkBH1IN1JyMpP46gVJ3Doth91F5aw/eIJN+SeDghXAkRPVHDlRzUfbiwHISI5lUHYyg7KTGJiVxIDMJDKSY1s9f5QQHxsIVn6/n735xWzYWcDWPYfZm1+M3x88Z1ZcVk5xWTlfbNofOJadkcKgvlkM7JtBTp8Mcnpn0L93j8AKwUqp8NEhPuClzw7y+JL1GCDLX8G7j98eVLVARY4xhoKyM2w5dJIdh0+TV1yJ19/838nEOBd9eyTQt0cCvdMT6J2WQO/0eLJS4gPLhLTGmepaduYdZVfeUXbmFbE3vzhQ2b4l0lMT6dsznT7Z6fTKTKV3Vio9M62vzLQk3O7Wt0mp7kqH+FrhVEU1dT8S05LjNTh1IBFhYJbVM7ru/H7UeH3sO1rB3mMV7D9WQV5JJTWes0sonan1kVdcSV4D81rJ8S6yUuLJTImjR1IsmclxpCfFkpYYa31PcJOS4A5KzEhMiGPSmIFMGjMQAK/XR/6RMvYXlLA3v5i8wuMUHC1rtJr6qfIznCo/w668ow2+np6aSEZaMpnpSfRITSIzPYm0lETSUxNJT0kkLSWB1OQEUpMSiI3VYKYURDhAichc4HeAC/izMeZXIa/HA38BJgOlwNeMMQft1x4Gvgn4gPuMMe9Hqp0ny+vnIs7lNaCiQbzbxdicdMbmWKWmfH5D0ckz5B+v4lBpFfnHqzhy4gxVtY33biprfFTWWOc3JSnORUqCm6R4F8nxbpLjrO3EOPsr1kVqr95c2L8vs2NduAROnKqg5Pgpjh0/RfHxkxQVn+Jo6almlwGpC2AHDzf/GcTFuklJiic5KZ6UxHiSEuJISoyzvifEkpAQR2J8LAnxbuLjYkmIjyU+1k1CvJu42Fji49zExbqIi7W+x7pdmk6vuqSIBSgRcQFPAlcAhcBaEVlmjNnhOO2bwAljzDARmQ88BnxNRMYA84GxQD/gQxEZYYwJ/0p5wMmK+mrcWekaoKKJK0bIyUwiJzOJmfYxYwynqjwcOVnN0VPVHLO/SsprKK2oxetr2bB1Va2vyUDXNDdINvTOpkcfMH4/Ho8XT62HWo+X2loP1TW11NR4qKmx0twFg9jfrX0gcKx+Gy9wGuR0NVAdWAIGY3CGmbr31G87Xwv+DFyuGGJdMbhcLtwuwe2OwRXjwu2KISZGcLticMUIrpgYXPZ2TEwMLpcQI4LLFWN9j7HOF8GxbZ0jMRAj1vtikMB5Is7vgiBWZREBEQK1FQPn4dyXwB82sBl0DoFzxPEp1J/rOMbZx2jomDTx2lmfv/O8Bk9rUkvfEo2/YMyeOJCE+PA/eO8UyR7UNGCfMSYPQEReAeYBzgA1D/ipvb0E+INY/yfmAa8YY2qAAyKyz77eF5Fo6Kny+gCV2UNTi6OdiNAjOY4eyXGM6R+c0FIXvEorajlRVcuJSg8nKms5fcbDqSoPp854OH3GS2WNl3BNv/oNIDG44+Jwx8UR+jSVMQaP14fH48Pj9VHrtb577S+Pz4/P68fj8+Hz+Wm0We35GWWwAt9ZWfbG/tKFIFXr/HNEny4doPoDBY79QuCCxs4xxnhF5BSQZR9fHfLe/qE3EJG7gbsBBg4c2OaGnq6sCWxn67MvXZozeDXFGENljY/yag9VNT4qa71UVvs44/FxptZHVa2XM7U+ajx+qj0+arz+wHat1299+azvzQU6EbGH21ryz83g8xm8fj8+n/0VtG3w+f34/da2328C28be9xm/tW2sr26SB6XOQZEMUA39vhf6T6Wxc1ryXowxTwFPgZXF19oG1hnaN42qWj8VZ2oY1FuX2TgXiAgpdrJEexhj8PkNHp/BYwcsrx1gPD7rNa/Pj9dv7OAS/OU3Bq+vLpgQOGZCtp3Bpm4/8B0wfvu7qV8+JvDdDmpen9WWQMDzG6v9PmewswKcv+7ezkDnB7/x220Ibpux209dG7EDozFW38xuJ472GevEwD/suoziuv/W/VkJeb2xffs2Zx9r4AXTyE5TP0RMS09soa6eQR3rjnwx7UgGqEJggGM/BzjSyDmFIuIG0oGyFr43bH6+cGbzJynVABGx5nRckIhm3ykVTpEMgWuB4SKSKyJxWEkPy0LOWQYstLdvAZYb69eKZcB8EYkXkVxgOPBlBNuqlFIqykSsB2XPKd0LvI+VZv6sMWa7iDwKrDPGLAOeAV60kyDKsIIY9nmvYSVUeIF7IpXBp5RSKjppJQmllFIdqqWVJHTJWKWUUlFJA5RSSqmopAFKKaVUVNIApZRSKippgFJKKRWVuk0Wn4iUAPntvEw2cDwMzYmkrtBG6Brt1DaGT1dop7YxfNrbzkHGmJ7NndRtAlQ4iMi6lqQ+dqau0EboGu3UNoZPV2intjF8OqqdOsSnlFIqKmmAUkopFZU0QAV7qrMb0AJdoY3QNdqpbQyfrtBObWP4dEg7dQ5KKaVUVNIelFJKqaikAUoppVRU0gAFiMhcEdktIvtE5KHObk9DRORZESkWkW2d3ZbGiMgAEVkhIjtFZLuI3N/ZbQolIgki8qWIbLbb+LPOblNTRMQlIhtF5K3ObktDROSgiGwVkU0iEpXLCYhIDxFZIiK77L+bMzq7TaFEZKT9GdZ9nRaRBzq7XaFE5EH73802EVksIgkRvd+5PgclIi5gD3AF1kq+a4EFxpgdndqwECJyCVAB/MUYc15nt6chItIX6GuM2SAiqcB64MZo+ixFRIBkY0yFiMQCq4D7jTGrO7lpDRKRHwBTgDRjzHWd3Z5QInIQmGKMidqHS0XkBeBTY8yf7cVTk4wxJzu7XY2xfyYdBi4wxrS3+EDYiEh/rH8vY4wxZ+w1+94xxjwfqXtqDwqmAfuMMXnGmFrgFWBeJ7fpLMaYT7AWdYxaxpgiY8wGe7sc2An079xWBTOWCns31v6Kyt/SRCQHuBb4c2e3pasSkTTgEqzFUTHG1EZzcLLNAfZHU3BycAOJIuIGkoAjkbyZBijrB2iBY7+QKPuh2hWJyGDgfGBN57bkbPaw2SagGPjAGBN1bbT9Fvg3wN/ZDWmCAf4pIutF5O7ObkwDhgAlwHP2UOmfRSS5sxvVjPnA4s5uRChjzGHgceAQUAScMsb8M5L31AAF0sCxqPyNuqsQkRTgdeABY8zpzm5PKGOMzxgzEcgBpolI1A2Zish1QLExZn1nt6UZM40xk4CrgXvsoeho4gYmAX8yxpwPVAJROc8MYA9B3gD8rbPbEkpEMrBGl3KBfkCyiHw9kvfUAGX1mAY49nOIcLe1O7PndV4HXjLG/L2z29MUe6hnJTC3k5vSkJnADfYczyvAbBH5a+c26WzGmCP292JgKdaQeTQpBAodveQlWAErWl0NbDDGHOvshjTgcuCAMabEGOMB/g5cGMkbaoCykiKGi0iu/dvLfGBZJ7epS7ITEJ4BdhpjftPZ7WmIiPQUkR72diLWP7pdnduqsxljHjbG5BhjBmP9nVxujInob6utJSLJdjIM9rDZlUBUZZkaY44CBSIy0j40B4iapJ0GLCAKh/dsh4DpIpJk/1ufgzXPHDHuSF68KzDGeEXkXuB9wAU8a4zZ3snNOouILAZmAdkiUgj8pzHmmc5t1VlmAv8CbLXneAB+Yox5pxPbFKov8IKdKRUDvGaMicoU7i6gN7DU+lmFG3jZGPNe5zapQd8HXrJ/Ac0D7uzk9jRIRJKwsom/3dltaYgxZo2ILAE2AF5gIxEueXTOp5krpZSKTjrEp5RSKippgFJKKRWVNEAppZSKShqglFJKRSUNUEoppaKSBijVLYlIbxF5WUTy7DI8X4jITW281uCOrCJvl+MZ01H3s+/5gJ3m3Nr3/baueoSI/NKucj0o5JwP7SoESrWKBijV7dgPEb4BfGKMGWKMmYz1sGtO57aseSLiMsbcFe4K8GJp6t/7A1jFP1tzzUxgul3IGGPMT4CngQdDTn0R+F5rrq0UaIBS3dNsoNYY8z91B4wx+caY30NgTajn7HWMNorIZfbxwSLyqYhssL+aLeMiIveJyA4R2SIir9jHUhzX3yIiN9vHF9jHtonIY45rVIjIoyKyBpghIitFZIrjtV+ItX7VahHpbR8fau+vtd9b0UDbBou1/tEfsR6uHCAifxKRdeJYC0tE7sOqrbZCRFbYx660e50bRORvdn3FULcAoQ/mvgvMtx+ErrMMq0KCUq2iAUp1R2OxfiA35h4AY8w4rB+cL4i18FoxcIVd/PRrwBMtuNdDwPnGmPHAd+xjj2BVeh5nH18uIv2Ax7CC50RgqojcaJ+fDGwzxlxgjFkVcv1kYLUxZgLwCfAt+/jvgN8ZY6bSdO3IkVhriJ1vL9/w78aYKcB44FIRGW+MecK+xmXGmMtEJBv4v8Dl9mexDvhBA9eeibXml9N8IAOrDA4AxpgTQLyIZDXRTqXOogFKdXsi8qTdA1lrH7oIa9gJY8wuIB8YgbU21NMishWrmnRL5oG2YJXR+TpW+Rew6vs9WXeC/QN6KrDSLrTpBV7CWqcIwIdVYLchtUBdKab1wGB7ewb1Fa9fbqJ9+SGLMX5VRDZglakZS8N/xun28c/sklULgUENnNcXaykLIFCJ+07gR8DtIecWY/XSlGqxc74Wn+qWtgM31+0YY+6xewV1S5I3tMQKWHMnx4AJWL+8VYeeICLPYa1zdcQYcw3WgoKXYC2R8IiIjLWvH1pDrLF7AlQbY3yNvOYx9fXIfLT+32ylo+25WMFjqjHmhIg8DzS0ZLdgrZPV3LDcmZD3fw34CGuBxR+LSKIx5oz9WoJ9vlItpj0o1R0tBxJE5LuOY84EgE+wf8MXkRHAQGA3kA4UGWP8WEVvnfMoABhj7jTGTDTGXGMnHQwwxqzAWliwB5AC/BO4t+49dgbbGqwhtWx7fmYB8HE7/oyrqQ/C81v4njSsgHXKnsu62vFaOZDquPZMERlmtz/J/pxC7QSGOfbvAxbZQekDrKBdl7TSBzjYwnYqBWiAUt2Q3eO4ESsgHBCRL4EXgB/bp/wRcNlDeSW8+mUAAAD7SURBVK8C3zDG1NjHF4rIaqwhv8qzrx7EBfzVvs5GrB/OJ4GfAxl2MsRmrLmdIuBhYAWwGWvNnzfb8cd8APiB/WfrC5xq7g3GmM12O7cDzwKfOV5+CnhXRFYYY0qAbwCLRWQLVsAa1cAl38aqsI+IXASUOlYCeJH6Yb7JWPNo3rOuoFQTtJq5Ul2Q/czSGWOMEZH5wAJjzLxOaMcq4Do7MDd2zu+AZcaYjzquZao70DkopbqmycAf7OGzk8C/dlI7fog1RNpogMLKUNTgpFpNe1BKKaWiks5BKaWUikoaoJRSSkUlDVBKKaWikgYopZRSUUkDlFJKqaj0/wOE+4wpJJGrVQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Pdf(germany)\n", "thinkplot.Pdf(argentina)\n", "thinkplot.decorate(xlabel='Goal-scoring rate (λ)',\n", " ylabel='PMF')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To answer the first question, \"How much evidence does this victory provide that Germany had the better team?\", we can compute the posterior probability that Germany had a higher goal-scoring rate:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "posterior prob Germany > Argentina 0.6983938606019376\n" ] } ], "source": [ "post_prob = germany.ProbGreater(argentina)\n", "print('posterior prob Germany > Argentina', post_prob)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the prior distributions, we would have said that Germany had a 50% chance of having the better team, or 1:1 odds. Based on the posteriors, we would say that Germany has a 70% chance. We can use the ratio of the prior and posterior odds to compute the Bayes factor, which measures the strength of the evidence." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "posterior odds Germany > Argentina 2.315582375066283\n", "Bayes factor 2.315582375066283\n" ] } ], "source": [ "prior_odds = 1\n", "post_odds = post_prob / (1 - post_prob)\n", "print('posterior odds Germany > Argentina', post_odds) \n", "k = post_odds / prior_odds\n", "print('Bayes factor', k) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Bayes factor is about 2.3, which is generally considered weak evidence.\n", "\n", "Now on to Step 4." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 4: Comparing posterior distributions\n", "\n", "**Exercise:** Write a few lines of code to \n", "\n", "1. Choose a random value of `lam` from the posterior distribution of each team.\n", "\n", "2. Choose a random number of goals for each team, conditioned on the value of `lam` you chose.\n", "\n", "3. Run that \"simulation\" many times and accumulate the distribution of wins, losses, and ties.\n", "\n", "Use the results to estimate the probability that Germany would win a rematch." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.418" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "gdr_goals = poisson.rvs(germany.Sample(1000))\n", "arg_goals = poisson.rvs(argentina.Sample(1000))\n", "np.mean(gdr_goals > arg_goals)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.329" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "np.mean(gdr_goals == arg_goals)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.253" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "np.mean(gdr_goals < arg_goals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of running simulations, you could compute the posterior predictive distributions explicitly.\n", "\n", "Write a function called `PredictiveDist` that takes the posterior distribution of $\\lambda$ and a duration (in units of games).\n", "\n", "It should loop through the hypotheses in `suite`, compute the predictive distribution of goals for each hypothesis, and assemble a \"meta-Pmf\" which is a Pmf that maps from each predictive distribution to its probability.\n", "\n", "Finally, it should use `MakeMixture` to compute the mixture of the predictive distributions." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "def PredictiveDist(suite, duration=1, label='pred'):\n", " \"\"\"Computes the distribution of goals scored in a game.\n", "\n", " returns: new Pmf (mixture of Poissons)\n", " \"\"\"\n", " metapmf = thinkbayes2.Pmf()\n", " for lam, prob in suite.Items():\n", " pred = thinkbayes2.MakePoissonPmf(lam * duration, 10)\n", " metapmf[pred] = prob\n", "\n", " mix = thinkbayes2.MakeMixture(metapmf, label=label)\n", " return mix" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "germany_pred = PredictiveDist(germany, label='germany')\n", "argentina_pred = PredictiveDist(argentina, label='argentina');" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAG4hJREFUeJzt3XucV3W97/HXB8QAJc2gffYGEiqSUGCEAdxe2mRZ+KjAvHQky6xsLOVYtrNjaVakHs/uZl664OWoO/OSWqLisSzcZYoxCkbgJaKMSc+J2OXemiTYZ//xW0w/hpEZxlkza2Zez8eDB7+11nd91+c3Im++39/6fVdkJpIkVc2g3i5AkqT2GFCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSSg2oiJgTEY9GxNqIOOMF2rwzItZExOqI+HaZ9UiS+o4oayWJiBgMPAYcBrQAy4H5mbmmrs0E4Abg0Mz8Y0S8IjN/v6N+R44cmePGjSulZklS+R544IE/ZOaojtrtUmINM4G1mbkOICKuA+YBa+rafBC4JDP/CNBROAGMGzeO5ubmEsqVJPWEiHi8M+3KnOIbDayv224p9tV7LfDaiPhpRCyLiDkl1iNJ6kPKHEFFO/vazifuAkwAZgNjgJ9ExH6Z+adtOopoApoAXvnKV3Z/pZKkyilzBNUCjK3bHgM80U6bWzJzc2b+GniUWmBtIzMXZWZjZjaOGtXhtKUkqR8ocwS1HJgQEeOB3wHHAu9q0+Z7wHzgyogYSW3Kb12JNUnSDm3evJmWlhY2bdrU26X0eUOHDmXMmDEMGTKkS+eXFlCZuSUiFgB3AoOBKzJzdUQsBJozc3Fx7M0RsQZ4Hjg9MzeWVZMkdaSlpYURI0Ywbtw4Itr7pEKdkZls3LiRlpYWxo8f36U+yhxBkZlLgCVt9p1d9zqBjxW/JKnXbdq0yXDqBhHBy1/+cjZs2NDlPlxJQpLaMJy6x4v9ORpQkqRKKnWKT5L6ugXnXNut/V181vxu7a8/M6B2wpk3rOrSeee+c3I3VyJJ29uyZQu77NJ//lrvP+9EkvqJz3/+81xzzTWMHTuWkSNHMn36dN7xjndwyimnsGHDBoYPH86ll17KxIkTOeGEE9hrr71YsWIF06ZNY8SIEfz617/mySef5LHHHuPLX/4yy5Yt44477mD06NHceuutDBkyhIULF3Lrrbfy7LPPcuCBB/LNb36TiGD27NnMmjWLpUuX8qc//YnLL7+cQw45hEMOOYSLLrqIhoYGAA466CC+/vWvM2XKlNJ+Dn4GJUkV0tzczE033cSKFSu4+eabW9cebWpq4qKLLuKBBx7gi1/8IieffHLrOY899hh33XUXX/rSlwD41a9+xe23384tt9zCu9/9bt7whjewatUqhg0bxu233w7AggULWL58Ob/4xS949tlnue2221r727JlCz/72c+44IIL+NznPgfAiSeeyJVXXtl6vb/85S+lhhMYUJJUKffccw/z5s1j2LBhjBgxgre//e1s2rSJe++9l2OOOYaGhgZOOukknnzyydZzjjnmGAYPHty6ffjhhzNkyBAmT57M888/z5w5tWVOJ0+ezG9+8xsAli5dyqxZs5g8eTI/+tGPWL16dev5Rx55JADTp09vbX/MMcdw2223sXnzZq644gpOOOGEcn8QOMUnSZXS3iOQ/vrXv7LnnnuycuXKds/Zbbfdttl+yUteAsCgQYMYMmRI6+3egwYNYsuWLWzatImTTz6Z5uZmxo4dy2c/+9ltVs7Yev7gwYPZsmULAMOHD+ewww7jlltu4YYbbuiRp0o4gpKkCjn44IO59dZb2bRpE08//TS33347w4cPZ/z48XznO98BaiH20EMPdfkaW8No5MiRPP3009x4442dOu/EE0/k1FNPZcaMGey1115dvn5nOYKSpB3o6dvCZ8yYwdy5c5k6dSp77703jY2N7LHHHlxzzTV8+MMf5pxzzmHz5s0ce+yxTJ06tUvX2HPPPfngBz/I5MmTGTduHDNmzOjUedOnT+elL30p73vf+7p03Z1V2hN1y9LY2Ji99cBCbzOX+r+HH36Y173udb1aw9NPP83uu+/On//8Z17/+tezaNEipk2b1qs1ATzxxBPMnj2bRx55hEGDOjcB197PMyIeyMzGjs51ik+SKqapqYmGhgamTZvGUUcdVYlwuvrqq5k1axbnnntup8PpxXKKT5Iq5tvf/nZvl7Cd448/nuOPP75Hr+kISpJUSQaUJKmSDChJUiUZUJKkSvImCUnaga5+veSFVPFrJ+eddx6f+tSnWrcPPPBA7r333l6sqMYRlCT1Mc8//3y39nfeeedts12FcAIDSpIq54gjjmD69Onsu+++LFq0CIDdd9+ds88+m1mzZnHfffexZMkSJk6cyMEHH8ypp57K2972NgCeeeYZ3v/+9zNjxgz2339/brnlFgCuvPJKjjzySObMmcOECRP4xCc+AcAZZ5zBs88+S0NDA8cdd1zrtQDuvvtuZs+ezdFHH83EiRM57rjjWtcKXLhwITNmzGC//fajqamp3TUEXywDSpIq5oorruCBBx6gubmZCy+8kI0bN/LMM8+w3377cf/999PY2MhJJ53EHXfcwT333MOGDRtazz333HM59NBDWb58OUuXLuX000/nmWeeAWDlypVcf/31rFq1iuuvv57169dz/vnnM2zYMFauXMk111yzXS0rVqzgggsuYM2aNaxbt46f/vSnwI4f19FdDChJqpgLL7yQqVOncsABB7B+/Xp++ctfMnjwYI466igAHnnkEV71qlcxfvx4AObP/9t6gd///vc5//zzaWhoYPbs2WzatInf/va3ALzxjW9kjz32YOjQoUyaNInHH3+8w1pmzpzJmDFjGDRoEA0NDZ16XEd38SYJSaqQu+++m7vuuov77ruP4cOHt4bM0KFDW5/5tKPptMzkpptuYp999tlm//3339/6GA3Y9lEaO9LeOR09rqO7OIKSpAp56qmneNnLXsbw4cN55JFHWLZs2XZtJk6cyLp161pHM9dff33rsbe85S1cdNFFrSG2YsWKDq85ZMgQNm/e3Okau/q4jp3lCEqSdqCnbwufM2cO3/jGN5gyZQr77LMPBxxwwHZthg0bxte+9jXmzJnDyJEjmTlzZuuxT3/603z0ox9lypQpZCbjxo3r8POhpqYmpkyZwrRp09r9HKqtrj6uY2f5uI2d4OM2pP6vCo/b6Iytj+TITE455RQmTJjAaaed1ttlbcfHbUjSAHPppZfS0NDAvvvuy1NPPcVJJ53U2yV1O6f4JKkPOu200yo5YupOjqAkqY2+9tFHVb3Yn6MBJUl1hg4dysaNGw2pFykz2bhxI0OHDu1yH6VO8UXEHOCrwGDgssw8v83xE4AvAL8rdl2cmZeVWZMk7ciYMWNoaWnZZnUGdc3QoUMZM2ZMl88vLaAiYjBwCXAY0AIsj4jFmbmmTdPrM3NBWXVI0s4YMmRI6woN6l1lTvHNBNZm5rrMfA64DphX4vUkSf1ImQE1Glhft91S7GvrqIj4eUTcGBFj2+soIpoiojkimh12S9LAUGZARTv72n7qeCswLjOnAHcBV7XXUWYuyszGzGwcNWpUN5cpSaqiMgOqBagfEY0BnqhvkJkbM/MvxealwPQS65Ek9SFlBtRyYEJEjI+IXYFjgcX1DSLi7+s25wIPl1iPJKkPKe0uvszcEhELgDup3WZ+RWaujoiFQHNmLgZOjYi5wBbg34ETyqpHktS3lPo9qMxcAixps+/sutefBD5ZZg2SpL7JlSQkSZVkQEmSKsmAkiRVkgElSaokA0qSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgElSaokA0qSVEmlBlREzImIRyNibUScsYN2R0dERkRjmfVIkvqO0gIqIgYDlwCHA5OA+RExqZ12I4BTgfvLqkWS1PeUOYKaCazNzHWZ+RxwHTCvnXafB/4F2FRiLZKkPqbMgBoNrK/bbin2tYqI/YGxmXnbjjqKiKaIaI6I5g0bNnR/pZKkyikzoKKdfdl6MGIQ8BXgnzvqKDMXZWZjZjaOGjWqG0uUJFVVmQHVAoyt2x4DPFG3PQLYD7g7In4DHAAs9kYJSRLALiX2vRyYEBHjgd8BxwLv2nowM58CRm7djoi7gY9nZnOJNQGw4Jxru3TeHq/dr5srkSS9kNJGUJm5BVgA3Ak8DNyQmasjYmFEzC3rupKk/qHMERSZuQRY0mbf2S/QdnaZtUiS+hZXkpAkVZIBJUmqJANKklRJBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgElSaokA0qSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVtMOAiohjit/H90w5kiTVdDSC+mTx+01lFyJJUr1dOji+MSKWAuMjYnHbg5k5t5yyJEkDXUcB9VZgGvCvwJfKL0eSpJodBlRmPgcsi4gDM3NDD9UkSVKHI6it9o6IRcDe9edk5pRSqpIkDXidDahrgNOBVcBfO9t5RMwBvgoMBi7LzPPbHP8QcArwPPA00JSZazrbvySp/+psQG3IzO1uktiRiBgMXAIcBrQAyyNicZsA+nZmfqNoPxf4MjBnZ64jSeqfOhtQn4mIy4AfAn/ZujMzb97BOTOBtZm5DiAirgPmAa0BlZn/Udd+NyA7WU+fteCca7t03sVnze/mSiSp2jobUO8DJgJD+NsUXwI7CqjRwPq67RZgVttGEXEK8DFgV+DQ9jqKiCagCeCVr3xlJ0uWJPVlnQ2oqZk5eSf7jnb2bTdCysxLgEsi4l3AWcB722mzCFgE0NjY2O9HWZKkzq/FtywiJu1k3y3A2LrtMcATO2h/HXDETl5DktRPdTagDgZWRsSjEfHziFgVET/v4JzlwISIGB8RuwLHAtvcaBERE+o23wr8srOFS5L6t85O8e30nXWZuSUiFgB3UrvN/IrMXB0RC4Hm4q7ABRHxJmAz8Efamd6TJA1MOwyoiBgKfAh4DbXvQF2emVs623lmLgGWtNl3dt3rj+xUtZKkAaOjKb6rgEZq4XQ4rscnSeohHU3xTdp6915EXA78rPySJEnqeAS1eeuLnZnakyTpxepoBDU1Irau9hDAsGI7gMzMl5ZanSRpwOrocRuDe6oQSZLqdfZ7UJIk9SgDSpJUSQaUJKmSOruShHrZmTes6tJ5575zZ9f4laRqcAQlSaokA0qSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKElSJZUaUBExJyIejYi1EXFGO8c/FhFrIuLnEfHDiNi7zHokSX1HaQEVEYOBS4DDgUnA/IiY1KbZCqAxM6cANwL/UlY9kqS+pcwR1ExgbWauy8zngOuAefUNMnNpZv652FwGjCmxHklSH1JmQI0G1tdttxT7XsgHgDvaOxARTRHRHBHNGzZs6MYSJUlVVWZARTv7st2GEe8GGoEvtHc8MxdlZmNmNo4aNaobS5QkVdUuJfbdAoyt2x4DPNG2UUS8CTgT+KfM/EuJ9UiS+pAyR1DLgQkRMT4idgWOBRbXN4iI/YFvAnMz8/cl1iJJ6mNKC6jM3AIsAO4EHgZuyMzVEbEwIuYWzb4A7A58JyJWRsTiF+hOkjTAlDnFR2YuAZa02Xd23es3lXl9SVLf5UoSkqRKMqAkSZVkQEmSKsmAkiRVkgElSaokA0qSVEkGlCSpkkr9HpSqZcE513bpvIvPmt/NlUhSxxxBSZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgElSaokA0qSVEkGlCSpkkoNqIiYExGPRsTaiDijneOvj4gHI2JLRBxdZi2SpL6ltICKiMHAJcDhwCRgfkRMatPst8AJwLfLqkOS1DftUmLfM4G1mbkOICKuA+YBa7Y2yMzfFMf+WmIdkqQ+qMyAGg2sr9tuAWaVeD2V5MwbVu30Oee+c3IJlUgaSMr8DCra2Zdd6iiiKSKaI6J5w4YNL7IsSVJfUGZAtQBj67bHAE90paPMXJSZjZnZOGrUqG4pTpJUbWUG1HJgQkSMj4hdgWOBxSVeT5LUj5QWUJm5BVgA3Ak8DNyQmasjYmFEzAWIiBkR0QIcA3wzIlaXVY8kqW8p8yYJMnMJsKTNvrPrXi+nNvUnSdI2XElCklRJBpQkqZJKneKTABacc22Xzrv4rPndXImkvsQRlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJPrBQlXXmDau6dN6575zc+tqHJUp9lyMoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSX4PSiqJ38GSXhwDSmpHd3xJWNKL4xSfJKmSSg2oiJgTEY9GxNqIOKOd4y+JiOuL4/dHxLgy65Ek9R2lTfFFxGDgEuAwoAVYHhGLM3NNXbMPAH/MzNdExLHA/wb+e1k1SX2BaxBKNWV+BjUTWJuZ6wAi4jpgHlAfUPOAzxavbwQujojIzCyxLkkvoK+Ho8Hcv0RZWRARRwNzMvPEYvs9wKzMXFDX5hdFm5Zi+1dFmz+06asJaCo29wEeLaXojo0E/tBhq/5nIL7vgfiewfc9kPTme947M0d11KjMEVS0s69tGnamDZm5CFjUHUW9GBHRnJmNvV1HTxuI73sgvmfwffd2HT2pL7znMm+SaAHG1m2PAZ54oTYRsQuwB/DvJdYkSeojygyo5cCEiBgfEbsCxwKL27RZDLy3eH008CM/f5IkQYlTfJm5JSIWAHcCg4ErMnN1RCwEmjNzMXA58K8RsZbayOnYsurpJr0+zdhLBuL7HojvGXzfA0nl33NpN0lIkvRiuJKEJKmSDChJUiUZUJ3Q0ZJN/VFEjI2IpRHxcESsjoiP9HZNPSkiBkfEioi4rbdr6SkRsWdE3BgRjxT/3f+xt2sqW0ScVvz5/kVEXBsRQ3u7pjJExBUR8fviu6db9+0VET+IiF8Wv7+sN2tsjwHVgbolmw4HJgHzI2JS71bVI7YA/5yZrwMOAE4ZIO97q48AD/d2ET3sq8D/zcyJwFT6+fuPiNHAqUBjZu5H7Wauqt+o1VVXAnPa7DsD+GFmTgB+WGxXigHVsdYlmzLzOWDrkk39WmY+mZkPFq//k9pfVqN7t6qeERFjgLcCl/V2LT0lIl4KvJ7anbVk5nOZ+aferapH7AIMK76HOZztv6vZL2Tmj9n+O6bzgKuK11cBR/RoUZ1gQHVsNLC+bruFAfIX9VbFKvP7A/f3biU95gLgE8Bfe7uQHvQqYAPwf4qpzcsiYrfeLqpMmfk74IvAb4Engacy8/u9W1WP+rvMfBJq/yAFXtHL9WzHgOpYp5Zj6q8iYnfgJuCjmfkfvV1P2SLibcDvM/OB3q6lh+0CTAO+npn7A89QwSmf7lR85jIPGA/8A7BbRLy7d6tSPQOqY51Zsqlfiogh1MLpmsy8ubfr6SEHAXMj4jfUpnMPjYhv9W5JPaIFaMnMraPkG6kFVn/2JuDXmbkhMzcDNwMH9nJNPen/R8TfAxS//76X69mOAdWxzizZ1O9ERFD7POLhzPxyb9fTUzLzk5k5JjPHUftv/aPM7Pf/qs7M/wesj4h9il1vZNtH4/RHvwUOiIjhxZ/3N9LPbwxpo36pufcCt/RiLe0qczXzfuGFlmzq5bJ6wkHAe4BVEbGy2PepzFzSizWpXP8DuKb4h9g64H29XE+pMvP+iLgReJDaXasr6APL/3RFRFwLzAZGRkQL8BngfOCGiPgAtbA+pvcqbJ9LHUmSKskpPklSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgGlASEino+IlcWq1d+JiOEvoq/ZW1c5j4i5O1rhvlgh/OQuXOOzEfHxHRzfLSJ+ULy+p1hLrjQd1SOVwYDSQPFsZjYUq1Y/B3yo/mDU7PT/D5m5ODPP30GTPYGdDqhO+EdgWbFczzOZuaWEa0i9yoDSQPQT4DURMa547tHXqH1Zc2xEvDki7ouIB4uR1u7Q+kywRyLiHuDIrR1FxAkRcXHx+u8i4rsR8VDx60BqX4Z8dTF6+0LR7vSIWB4RP4+Iz9X1dWbx3LG7gH1oR0S8uvji9LeAdwEPAFOL/rdb7DMiPl3U/YPieUcfL/Y3RMSyoobvbn0WUER8sKjtoYi4qb2RZkScGhFrinOv68LPX+oUA0oDSjEVdjiwqti1D3B13QKpZwFvysxpQDPwseIhdpcCbwcOAf7bC3R/IfBvmTmV2jp2q6ktuPqrYvR2ekS8GZhA7TEuDcD0iHh9REyntrTS/tQCcEZ7F8jMX2VmA7VgmglcDXyg6H+btdQiohE4qq7PxrrDVwP/MzOnFD+LzxT7b87MGcV7eBj4QDtlnAHsX5z7oXaOS93CpY40UAyrW7LpJ9TWGfwH4PHMXFbsP4DaQyl/WluajV2B+4CJ1BYV/SVAsXhsUzvXOBQ4HiAznweeaucppW8ufq0otnenFlgjgO9m5p+La3S03uMrMnNjREymFp7tORi4JTOfLfq8tfh9D2DPzPy3ot1VwHeK1/tFxDnUpiZ3p7bEV1s/p7Yk0veA73VQp9RlBpQGimeLkUerIoSeqd8F/CAz57dp10D3PWIlgP+Vmd9sc42PduYaEfENasEzpgjcCcDtEXFVZn6lnWvtrCuBIzLzoYg4gdr6bW29ldrDDecCn46Iff0MTGVwik/6m2XAQRHxGoBilevXAo8A4yPi1UW7+S9w/g+BDxfnDo7aU2r/k9roaKs7gffXfbY1uvjs6MfAOyJiWESMoDaduJ3M/BDwOeDz1J6Aensxvdc2nADuAd4eEUOL67216OMp4I8RcUjR7j3A1tHUCODJqD1q5bi2HRY3kozNzKXUHuq4daQldTtHUFIhMzcUo4ZrI+Ilxe6zMvOxiGiiNlL5A7W/+Pdrp4uPAIuK1aGfBz6cmfdFxE8j4hfAHcXnUK8D7itGcE8D787MByPiemAl8Di1acgX8k/UPkM6hL8FS3vvZ3kxVfhQ0Wcz8FRx+L3AN4qbIOpXLv80tScnP07ts6kRbGsw8K1imjCArwyQR8OrF7iaudSPRcTumfl0EUQ/Bpoy88HerkvqDEdQUv+2KCImAUOBqwwn9SWOoCRJleRNEpKkSjKgJEmVZEBJkirJgJIkVZIBJUmqpP8ChDnstLTWL1oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Hist(germany_pred, width=0.45, align='right')\n", "thinkplot.Hist(argentina_pred, width=0.45, align='left')\n", "thinkplot.decorate(xlabel='Predicted # goals',\n", " ylabel='Pmf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the predictive distributions, we can compute probabilities for the outcomes of a rematch." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior prob Germany wins rematch 0.44773969920578044\n", "Posterior prob tie 0.3270825241713642\n", "Posterior prob Argentina wins rematch 0.2251777766228554\n" ] } ], "source": [ "win = germany_pred.ProbGreater(argentina_pred)\n", "lose = germany_pred.ProbLess(argentina_pred)\n", "tie = 1 - (win + lose)\n", "\n", "print('Posterior prob Germany wins rematch', win)\n", "print('Posterior prob tie', tie)\n", "print('Posterior prob Argentina wins rematch', lose)" ] }, { "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.5" } }, "nbformat": 4, "nbformat_minor": 1 }