{ "metadata": { "name": "", "signature": "sha256:cd9f5b672c86961f2fe931c3fa564d2773377f93b082d4852ed5339e7980c4dc" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The World Cup Problem: Germany v. Argentina\n", "-------------------------------------------\n", "\n", "Allen Downey\n", "\n", "This notebook contains a solution to a problem I posed in my Bayesian statistics class:\n", "\n", "> In the final match of the 2014 FIFA World Cup, Germany defeated Argentina 1-0. How much evidence \n", "> does this victory provide that Germany had the better team? What is the probability that Germany\n", "> would win a rematch?\n", "\n", "Scoring in games like soccer and hockey can be (reasonably) well 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 is stationary; 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", "My solution uses the ThinkBayes2 framework, which is described in [Think Bayes](http://thinkbayes.com), and summarized in [this notebook](http://nbviewer.ipython.org/github/AllenDowney/ThinkBayes2/blob/master/code/framework.ipynb).\n", "\n", "I'll start with Step 2.\n", "\n", "### 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", "collapsed": false, "input": [ "# first a little house-keeping\n", "from __future__ import print_function, division\n", "% matplotlib inline" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "import thinkbayes2\n", "\n", "class Soccer(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: interarrival time in minutes\n", " \"\"\"\n", " goals = data\n", " lam = hypo\n", " like = thinkbayes2.EvalPoissonPmf(goals, lam)\n", " return like" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "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'll start with an unrealistic uniform distribution and update it with fake data until the mean matches the observed rate for a single team, 1.34 goals per game." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy\n", "import thinkplot\n", "\n", "hypos = numpy.linspace(0, 12, 201)\n", "suite = Soccer(hypos)\n", "suite.Update(0.33) # fake data chosen by trial and error to yield the observed prior mean\n", "\n", "thinkplot.Pdf(suite)\n", "suite.mean()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 3, "text": [ "1.3395835151444653" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl81PWdx/HX5AaChPtIAuEKtxLOACLBo1JssbrtWt0W\nq/ahu4rttttW3e22YC+1a2tdW0tXbVnrqlttWbAcohJFipFwQziDMQRIwhWOQO7ZP36/TH4ZZjKT\nkOT3+828n4/HPPL7/ub7m/lMi9/P7/c9fj8QEREREREREREREREREREREREREZEoNA/YBxwEHglS\n51nz/R1AlrkvCcgDtgMFwM8s9RcDJcA28zWvvYMWEZErFwscAjKAeIwGfYxfnfnAKnN7OvCR5b2u\n5t84c/8ss/xD4NvtH66IiLRGTIj3p2EkgSKgFngNuNWvzgJgmbmdB6QA/c3yRfNvAkZCOWM5ztOm\niEVEpN2ESgKpwBFLucTcF6pOmrkdi3H1UAasx+gWavQwRvfRixiJQ0REOlmoJOAN83P8z+obj6sH\nJmIkheuAHHP/88BQ873jwNNhfo+IiLSjuBDvHwXSLeV0jDP9luqkmfuszgJ/BaYAuUC55b0XgJWB\nvnz48OHewsLCECGKiIhFITAi3MqhrgTygZEYA8MJwB3ACr86K4CF5nY2UIHR/dOHpm6eLsBNGDOB\nAAZajr8N2BXoywsLC/F6vRH5+uEPf2h7DPp9+n36fZH3AoaHaNebCXUlUAcsAtZi9O+/COwFHjDf\nX4oxM2g+xgByJXCPpaFfhpFoYoCXgXfN957E6AryAp9YPk9ERDpRqCQAsNp8WS31Ky8KcNwuYFKQ\nz1wYZL+IiHSiUN1B0kFycnLsDqFD6fe5m35f9HD6XH2v2cclIiJh8Hg80Iq2PZzuIMepr29gxdrN\nHDl2CoDZ08eQNWGozVGJiLiPK68EfvfyOpa9tt5Xjo2L5ZePf43J17RqUFxEJOK09krAdWMCR0tP\n8z9/3tBsX31dPf/60//h0yMnbIpKRMSdXJcEfvPSGmpr6nzl3r26A3DhwiV+8sybaAxBRCR8rkoC\nRcXl5G7c7Sv/9j8e4KkfLCQu3hja2LOvmLytB+0KT0TEdVyVBPJ3NN1CYua00UwYM4TRI1P5/Gem\n+Pb//tX3dDUgIhImVyWB3fuKfdvWQeCvfuk639XA7r3FzZKFiIgE56oksKugKQlcPWawb7t/3xQ+\nd9NkX3n5qo87NS4REbdyTRIoP3mW0nLjmTSJiQlkDh/U7P2/uyXbt70hby8VZys7NT4RETdyTRLY\nVfCpb3tsZhpxcbHN3h+W0Z+xo4w7WtfX1bM2d3unxici4kauSQI791q6gsYNCVjnFkuX0F/XbdEA\nsYhICK5JArstSWDCmMBJ4MbrriYxMQGAwk9KKSwq7ZTYRETcyhVJwOv18klx08PIxoz0f8yxIblb\nEtdOH+0r527c0+GxiYi4mSuSwLnzl6iurgGgS5dEelzVNWjd668d79teb1lYJiIil3NFEig7UeHb\n7tenR+MNkgLKnpzp6xIqKi7nk0/LOjw+ERG3ckUSKD951rfdr2+PFusmJSUwY0qmr5z7N3UJiYgE\n44okUFredCUwIEQSAJhr6RJ6f1NBh8QkIhIJXJEEml0J9AmdBGZMySTWXEdwsPAYJ06d67DYRETc\nzHVJoH/flJD1u3VNYuK4DF85b8uBjghLRMT1XJEEyk6EPybQyDousClfSUBEJBBXJIHyE627EgBj\nllCjzdsPUVdX3+5xiYi4XThJYB6wDzgIPBKkzrPm+zuALHNfEpAHbAcKgJ9Z6vcC1gEHgLeBoC17\nfX0D5acsSSCMMQGAjMH9GNCvJwCVlVXstNx7SEREDKGSQCzwHEYiGAvcCYzxqzMfGAGMBO4Hnjf3\nVwFzgYnA1eb2LPO9RzGSQCbwrlkO6NSZ8zTUNwCQ0qMbiYnxYfws42HL2ZYuoY80LiAicplQSWAa\ncAgoAmqB14Bb/eosAJaZ23kYZ/X9zfJF828CRkI5E+CYZcAXggXQmjUC/mZMGeXb1riAiMjlQiWB\nVOCIpVxi7gtVJ83cjsXoDioD1mN0C4GRJBqX8pbRlDQuU1befLVwa0y+ZpjviWOHi0qbrTwWERGI\nC/F+uPdi9r+PQ+Nx9RjdQT2AtUAOkBugbtDv+c2vf0nxgUMAnB+VHGY4hi5JCWRNGMpm8+Hzm/IP\n8IXPTmvVZ4iIOFlubi65ubltPj5UEjgKpFvK6Rhn+i3VSTP3WZ0F/gpMxkgCZcAAoBQYCJQTxHU3\nfpHSCxsBmDXz2hDhXm7GlExLEtivJCAiESUnJ4ecnBxfecmSJa06PlR3UD7GgG8GRr/+HcAKvzor\ngIXmdjZQgdHI96Fp1k8X4CaMrqHGY+42t+8GlgcL4Pz5S77tlu4eGox1XCB/eyE1NXWt/gwRkUgV\nKgnUAYswunIKgNeBvcAD5gtgFXAYYwB5KfCguX8g8B5Gw58HrMSYCQTwBEZSOABcb5YDunCxyred\n3C0prB9lNTi1D4MG9gagqqqG3fuKQxwhIhI9QnUHAaw2X1ZL/cqLAhy3C5gU5DNPAzeG8d1cvFjt\n2+7WhiQAMG3icJYfPwVA/o5CJl09rE2fIyISaRy/YvhCZdOVQLcuiW36jKlZI3zb+dsLrzgmEZFI\n4fgkUHnpyq8EsiYMA/NBNAUHSpolFhGRaOb8JGC9EujatiuBHld1ZdTwQQB4GxrYtutwu8QmIuJ2\nzk8CljGB5K5tuxIAmDJxuG97s7qEREQAFySBmppaAGJiY8K+b1AgUyZqXEBExJ/jk0Cj5K5JLT5g\nPpRrxg4hPsGYDPXpkfJm9yQSEYlWrkkCbR0PaJSYGM+E0YN95a07NS4gIuKeJNDGmUFWzcYFth26\n4s8TEXE79ySBNq4RsJpqHRfYUYjXG+798UREIpN7kkA7XAmMGpHq+5yTp85RXHLyij9TRMTN3JME\nrnBMACA2NqbZLSM2b1eXkIhENxclgSu/EgCYck3TuED+Dk0VFZHo5qIkcOVXAtB8XGDrzsPUm88v\nFhGJRq5JAm25jXQgg9P60Nd8TGVlZRX7Dvk//0ZEJHq4Jgl0bYfZQQAej4fJlnGBfI0LiEgUc00S\naK8rAdCtpUVEGkVlEphsGRzeubeYqqqadvtsERE3cU0SaK/uIIC+va9iSHo/AOpq69i5V4+cFJHo\n5Jok0J5XAtD8FhIaFxCRaOWaJNAeK4atpjZ7voCSgIhEJ9ckgSt5oEwgWROG4Ykxfv6BwuOcPXex\nXT9fRMQNXJME2muxWKPkbkmMzUwzCl6vbi0tIlHJFUkgMTGBuLjYdv9cjQuISLQLJwnMA/YBB4FH\ngtR51nx/B5Bl7ksH1gN7gN3ANyz1FwMlwDbzNa+lANr7KqCR9T5Cm3UfIRGJQqGSQCzwHEYjPRa4\nExjjV2c+MAIYCdwPPG/urwW+BYwDsoGHgNHme17gFxgJIwtY01IQ7T0o3Gj86MEkJiYAcPTYKY6X\nnemQ7xERcapQSWAacAgowmjUXwNu9auzAFhmbucBKUB/oBTYbu6/AOwFUi3Hhf3A4PZ4oEwgCQlx\nXD1uiK+8RVcDIhJlQiWBVOCIpVxC84Y8WJ00vzoZGGf8eZZ9D2N0H72IkTiC6pKUECLMtptq7RLS\nLSREJMrEhXg/3Ocv+p/VW49LBt4AvolxRQBGl9Hj5vaPgKeB+wJ9cPGBDdSd283ixSXk5OSQk5MT\nZkjhsQ4ObzEfOenxhH2RIiJiq9zcXHJzc9t8fKgkcBRjgLdROsaZfkt10sx9APHAm8AfgeWWOuWW\n7ReAlcECGJw5m6lZI1i8+N4QobbNyGEDuap7V86dv8iZigsc/rSM4RkDOuS7RETam//J8ZIlS1p1\nfKjuoHyMAd8MIAG4A1jhV2cFsNDczgYqgDKMq4MXgQLgGb9jBlq2bwN2tRhkTMfNZI2JiWl2Qzl1\nCYlINAnVutYBi4C1GI356xgDvA+YL4BVwGGMAeSlwIPm/lnAV4C5XD4V9ElgJ8aYwByMWURBdXT3\njLVLaPM2rRcQkegRqjsIYLX5slrqV14U4LgPCZ5kFgbZH1BsbMeuabMmge27P6Gurr5DFqeJiDiN\nK1YMx8R07JVA6oBeDOzfE4Cqqhr27D8S4ggRkcjgjiTQwd1BHo+n2biAnjYmItHCHUmgAweGGzV7\n5KQWjYlIlHBFEujoMQGASZaHz+/Zf4SLl6o7/DtFROzmiiTQGWu3eqUkM3yosT6gvq6e7buLOv5L\nRURs5ookENsJ3UEAUyaqS0hEoosrkkBHzw5qNNWSBLReQESigUuSQOeEec24IcSa6wMOF5Vy6sz5\nTvleERG7uCQJdM6VQNcuiYwf1XQbJD1yUkQinTuSQCfe1bP5Iyc1LiAikc0VScDTSd1B0HxwePO2\nQ3i94d5NW0TEfVyRBDqrOwhgbGaa73GWZScqKC452WnfLSLS2VyRBDpjsVijuLjYZreQ2LTlQKd9\nt4hIZ3NFEujsJ31NnzTSt/3x1oOd+t0iIp3JFUmgsxaLNZpmuY/Qtl2fUF1d26nfLyLSWVyRBDpz\nTABg0IBeDE7rC0BNTS3b9xR16veLiHQWJYEgplm6hPLUJSQiEcodSaCTxwQAsidn+rbzNDgsIhHK\nHUmgk8cEALLGZxCfYDx9s6i4nLITFZ0eg4hIR3NJEuj8K4GkpAQmjsvwldUlJCKRSEmgBdaponlb\nlAREJPK4Igl09hTRRtMnNY0LbN5+iPr6BlviEBHpKK5IAnaMCQAMHdKPvn16AFBZWcWe/UdsiUNE\npKOE07rOA/YBB4FHgtR51nx/B5Bl7ksH1gN7gN3ANyz1ewHrgAPA20BKSwF4bOoO8ng8zbqE/rZ5\nvy1xiIh0lFBJIBZ4DiMRjAXuBMb41ZkPjABGAvcDz5v7a4FvAeOAbOAhYLT53qMYSSATeNcsBw/S\npiQAMGvaaN/2xo/32RaHiEhHCJUEpgGHgCKMRv014Fa/OguAZeZ2HsZZfX+gFNhu7r8A7AVSAxyz\nDPhCS0HYNSYAxvMF4uKNqaKHi0opLddUURGJHKFa11TA2hFeQlND3lKdNL86GRjdRHlmuT9QZm6X\nmeWgOvsGclZduySSNWGor7wpX11CIhI54kK8H+4TVfxbaetxycAbwDcxrggCfUfQ7yk+sIE/vXqC\nXflryMnJIScnJ8yQ2s/MqaPYbK4T2PjxPm6bP73TYxARCSQ3N5fc3Nw2Hx8qCRzFGOBtlI5xpt9S\nnTRzH0A88CbwR2C5pU4ZMACjy2ggUB4sgMGZs7lr4W0suHlqiFA7zsypo/jV0rcA2LrTuKtoYmK8\nbfGIiDTyPzlesmRJq44P1R2UjzHgmwEkAHcAK/zqrAAWmtvZQAVGI+8BXgQKgGcCHHO3uX03zRPE\n5UHa2B0EkDawt++uotXVNWzRA+hFJEKESgJ1wCJgLUZj/jrGAO8D5gtgFXAYYwB5KfCguX8W8BVg\nLrDNfM0z33sCuAljiuj1Zjl4kDYODDeaOXWUb1uzhEQkUoTqDgJYbb6slvqVFwU47kOCJ5nTwI1h\nfDdg7xTRRrOmjea1v3wIwKbN+/F6vbYOWIuItAf7T7HD4IQkcPXYIc0eQH/407IQR4iIOJ9LkoD9\nYcbFxTZbPawuIRGJBPa3rmGwc7GYlXX18Id5SgIi4n7OaF1DcErX+4wpo4iJNf4n27OvmBOnztkc\nkYjIlXFFEoiNdUaYPa7qStb4ptXDH+bttTEaEZEr54zWNQSPxzlhzs4e69t+f1OBjZGIiFw557Su\nLXDC7KBG181ouonq1p2HOX/hko3RiIhcGSWBVurfN4VRI4x76NXX1euGciLiaq5IAk4ZE2g0Z2ZT\nl9AHmzQuICLu5azWNQinrcydPb0pCXy05QDV1bU2RiMi0nauSAJ230DO39Ah/UhL7QPApUvV5O8o\ntDkiEZG2cUUScFp3kMfj4TrLLKENH2mWkIi4k7Na1yCc1h0EMDu7aZbQho/2Ul/fYGM0IiJt44ok\n4LQrAYDxo9Pp3as7ABVnK9m6S88YEBH3cV7rGoDTxgTAuKnd3Gsn+MrrP9xtYzQiIm3jjiTgoHUC\nVjfMbkoCuRv3qEtIRFzHJUnAmWGOH51O3z49ADh7rlKzhETEdZzZuvpxyq2k/RldQuN95fc27LIx\nGhGR1nNm6+rHgUMCPtYuofc37aGurt7GaEREWscVSSDGgbODGo0blU7/vikAnD9/ic3b1SUkIu7h\n3NbVwqljAmCsYbB2Ca3/UF1CIuIezm1dLZw6O6hR8y6hAmpr62yMRkQkfO5IAk4eFADGZKYxsH9P\nAC5cuETe1oM2RyQiEp5wksA8YB9wEHgkSJ1nzfd3AFmW/S8BZYB/H8lioATYZr7mtRikw68EPB4P\nN1x3ta+8+t1tNkYjIhK+UEkgFngOo5EeC9wJjPGrMx8YAYwE7geet7z3ewI38F7gFxgJIwtY02KQ\nDh4TaDRv7kTf9sbN+/TEMRFxhVCt6zTgEFAE1AKvAbf61VkALDO384AUYIBZ3gCcCfLZYZ/eO/1K\nAGDokP5kjhgEQG1NHe/pNhIi4gKhkkAqcMRSLjH3tbZOIA9jdB+9iJE4gnLqYjF/N89t6glb+566\nhETE+eJCvO8N83P8T9VDHfc88Li5/SPgaeC+QBWLD2zgiSd+SlJiPDk5OeTk5IQZUue7ac7V/Pql\n1TTUN7BjTxFHS0+TOqCX3WGJSATLzc0lNze3zceHSgJHgXRLOR3jTL+lOmnmvpaUW7ZfAFYGqzg4\nczbf//6/k9wtKcRH2q93z+5MyxrJR+bD599ev5177rze5qhEJJL5nxwvWbKkVceH6mfJxxjwzQAS\ngDuAFX51VgALze1soAJjRlBLBlq2b+Py2UPNg3TBmECjz97Q1CW0Zv12vN5wL6ZERDpfqCRQBywC\n1gIFwOvAXuAB8wWwCjiMMYC8FHjQcvyrwN+ATIxxg3vM/U8COzHGBOYA32oxSIevE7CaPX0MXbsm\nAlBy9CQFB/wvnEREnCNUdxDAavNltdSvvCjIsXcG2b8wyP6A3DBFtFFiYjxzZ43nr+u2ALDqna2M\nG5Ue4igREXu4onV14uMlWzLv+qYuoXXv76CqqsbGaEREgnNF6+qi3iAAsiYMJS21DwCVlVVaMyAi\njuWKJOCm7iAwbiOx4DNTfOWVb+fbGI2ISHCOb12d/CyBlsy7IYvYuFgAdu4p4pNPQ02YEhHpfI5v\nYT1u6wsy9e7ZndnTm26ztNIcKBYRcRLHJwG33DIikM/f3NQltObdbdTU6DkDIuIsjm9h3bRQzN+0\nrBEM6Gc8Z+DsuUo++KjA5ohERJpTEuhAMTExfO4zk33lFWs22xiNiMjlXJAEHB9ii265aTIe8zds\n2VGoAWIRcRTHt7BuHhMA6NenB7OzmwaI33jrIxujERFpzvEtrEsnBzXzpc/P8G2vfnebnjomIo7h\n+CTgtltGBJI1YSjDMoyHrVVX1/CWpouKiEM4voV1+5gAGGsdvrSg6Wrgz299RH19g40RiYgYHN/C\nuuk20i35zJxruKp7VwCOlZ5mk/ngGREROzk/Cbh4iqhVUlJCs8Vjf1rxNxujERExOD4JeCKgO6jR\n7bdk+35P/vZCDhdpuqiI2MvxLWykXAkADOiXwpwZY33lV978wMZoRETckAQiZEyg0Z23X+vbXvfB\nTo6XnbExGhGJdo5PApEwRdRq/OjBZE0YBkB9XT2vLd9oc0QiEs0c38K69VbSLfnKl67zba9cm0/F\n2UoboxGRaOb4JBBpVwIA0yeNZOTwQYCxeOyNlZtsjkhEopXjW9hIGxMA4+rmK19suhp4461NXLxU\nbWNEIhKtnJ8EImh2kNXcWeMZNLA3AOfPX+L/Vn9sc0QiEo3CSQLzgH3AQeCRIHWeNd/fAWRZ9r8E\nlAG7/Or3AtYBB4C3gZSgAUbQOgGr2NgY7rLMFHrlzQ1cqqqxMSIRiUahWthY4DmMRDAWuBMY41dn\nPjACGAncDzxvee/35rH+HsVIApnAu2Y5cAAROCbQaP4Nk+jXtwcAZyou8KbGBkSkk4VqYacBh4Ai\noBZ4DbjVr84CYJm5nYdxVj/ALG8AAk2Etx6zDPhCsAAiszPIkJgYz913zPWVX3lzA5UXq2yMSESi\nTagkkAocsZRLzH2treOvP0Y3Eebf/kEDjNDuoEa33DiJQQN6AXDu/EX+9/90TyER6TxxId73hvk5\n/ifs4R7XWDdo/c0bV7B4cQkAOTk55OTktOKjnS8+Po6vfXkuP33mTQBe/cuH/N3nsn13HBURaUlu\nbi65ubltPj5UEjgKpFvK6Rhn+i3VSTP3taQMo8uoFBgIlAerOGPOF1i8+N4QH+du867P4uU3PuBI\nyQkqK6t4bflG7v/qTXaHJSIu4H9yvGTJklYdH6qvJR9jwDcDSADuAFb41VkBLDS3s4EKmrp6glkB\n3G1u3w0sDxpghHcHgTH4fd9d1/vKry/fyMnT52yMSESiRagWtg5YBKwFCoDXgb3AA+YLYBVwGGMA\neSnwoOX4V4G/YcwCOgLcY+5/ArgJY4ro9WY5cIARuk7A3w2zJ/geQVlVVcN/vfyOzRGJSDRwegvr\n/e6S/+apH3zV7jg6Rd7Wg3z7339vFDwe/vDsIkYOG2hvUCLiKub91sJu2x3f1xItVwJg3FMoe8oo\no+D18p8vrMLrbc0Yu4hI6zg+CcRGwZiA1aJ75xFjLpDbsqOQjR/rWcQi0nEc38JG4q2kWzJ0SH9u\nnTfNV/71S6upq6u3MSIRiWSOTwKRfNuIYO6763q6dUsCoLjkhBaQiUiHcXwLG01jAo16piRzz5eb\nbifxwivvUlpeYWNEIhKpnJ8Eoqw7qNGXFsz0TRmtrq7hl0tX2hyRiEQi5yeBKLwSAIiLi+W7DzXd\nq+/Dj/ay4aMCGyMSkUjkgiTg+BA7zNVjh7Bg3lRf+Ze/fUtPIBORduX4FjZarwQa/ePdN5PSoxsA\nZScq+N3L62yOSEQiifOTQJSOCTTqcVVXHv76fF/5Tys2sXXnYRsjEpFI4vwkEIVTRP3dPHci0ydn\nGgWvl58+86YePiMi7cLxLWy0LRYLxOPx8Og3biM5uQsAx8vO8NyLq22OSkQigeOTQDQuFgukX58e\n/Ms/fd5XXrFmM5vyD9gYkYhEAse3sNE+JmB105xryJk13lf+2a/e5HTFBRsjEhG3c34SiOIpov48\nHg/feXABPVOSATh1+jyP/8f/Ul/fYHNkIuJWjm9ho32KqL+eKcl8/9tf9JU3bzvEy39638aIRMTN\nHJ8ENCZwuezJmSy8I8dXfuGVdzVtVETaxPEtrK4DAvv6P9zIxPFDAfA2NLD4569z6sx5m6MSEbdx\nfBLQmEBgsbExLP7eHb7VxKdOn+exH79CdXWtzZGJiJs4voVVd1BwfXtfxQ+/ewceM1Hu2VfMk88t\n1yMpRSRsjm9hNUO0ZdOyRvDw1z/rK699bxuvvPGBjRGJiJs4PglE2zOG2+LvF8zk8zdP8ZWfX/Y2\nH2zSbadFJDTHt7AaEwjN4/HwL/+0wDdQjNfLD596nW27PrE3MBFxvHBa2HnAPuAg8EiQOs+a7+8A\nssI4djFQAmwzX/OCfblH6wTCEh8fx0/+9S4GDewNQE1NLY/86GUOHj5uc2Qi4mShkkAs8BxGIz0W\nuBMY41dnPjACGAncDzwfxrFe4BcYCSMLWBM0QCWBsKX06MYzP76H3r26A1BZWcW3f/AHSo6fsjky\nEXGqUElgGnAIKAJqgdeAW/3qLACWmdt5QAowIIxjw2rdNSbQOqkDevGLx79Gt25JAJw+c55//reX\nOF52xubIRMSJQrWwqcARS7nE3BdOnUEhjn0Yo/voRYzEEZBuJd16I4YO5KkffJWEhHjAuPX0okdf\n4GjpaZsjExGniQvxfrgTzlvbUj8PPG5u/wh4GrgvUMU/vfoiu/KN3qKcnBxycnJa+VXRaeL4ofzk\nX+/isZ+8Ql1tHaXlZ3j40Rd49mf3kWaOG4iI++Xm5pKbm9vm40M13tkYg7iNA7ePAQ3Ak5Y6vwVy\nMbp7wBgIngMMDeNYgAxgJTAhwPd7/2/Nxyy4eWqAtyQcm/IP8NhP/khtTR0Affv04Jkf3UPG4H42\nRyYiHcHsPQn7xDxUd1A+xoBvBpAA3AGs8KuzAlhobmcDFUBZiGMHWo6/DdgVNEB1B12RGVMyeeoH\nC31dQydOnuUfv7eUHXuK7A1MRBwhVBKoAxYBa4EC4HVgL/CA+QJYBRzGGAReCjwY4lgwrgZ2YowJ\nzAG+FTRADQxfsWlZI/j54oUkJSUAcP78Jb75/ZfI3bjb5shExG5OP832rn53K/OuzwpdU0Lad/Ao\n31m8jDONTyPzeFh032f58hdmaQBeJEK0d3eQ7bROoP2MHpnK0qf/kbTUPsYOr5fnXljF40//iaqq\nGnuDExFbuCAJOD5EV0kd0IvfPnU/48cM9u17e/12/umR32ktgUgUcnwLq8Vi7a9nSjL/+dOvN7vp\n3IFDx7j3m7/WjedEoozjW1h1VXeMhIQ4Hnn4Nr7z0K3ExsUCcO78RR778R/5+XPL1T0kEiUcnwT0\nUJmO4/F4uG3+dJ772dfp37dp0fby1R9z37d+Q8H+Iy0cLSKRwPEtrMfj+BBd7+qxQ/jDfy4iZ9Z4\n376i4nLu/85Sfv3SGj2yUiSCOb2zxbvx433MnDrK7jiigtfrZcXafH71u79SXd3UHZQ6qDfffehW\npk4cYWN0IhKO1k4RdXwS2JS/n+zJmXbHEVWOlp7myWf/wpYdhc3258waz0P3zmPQgF42RSYioURc\nEvh420Gdgdqg8arg1y+tprKyyrc/ISGeO2+/lrtun02yebtqEXGOiEsCm7cfYso1w+2OI2qdOHWO\n3/5hLWve29Zs/1Xdu7Lw7+dw+y3ZJCbG2xSdiPiLuCSwZUchk64eZnccUW9nwac8s/Qt9h862mx/\nn95Xcefts7l13lS6mPcmEhH7RFwS2L77E64Zl2F3HALU1zew7v0dvPjKuxzze0DNVd278sXPz+D2\nW6bTMyXLtLWWAAAIQklEQVTZpghFJOKSwI49RVw9dojdcYhFbW0dK9bms+z19Zw6fb7Ze3Hxcdw4\newK33TKdcaPSdWM6kU4WcUlg975ixo1KtzsOCaC6upa/vrOVV/+84bIrA4DMEYO4bf50bpg9gW5d\nNYgs0hkiLgkU7D/CmMw0u+OQFtTXN/DOBzt5Y+WmgKuMExMTuHb6aD6Tcw3TJ40kPj7UU01FpK0i\nLgnsO1jCqBH+z7YXp9p38Ch/WZXH27k7qKm5fKVx9+5duP7aCVw3YyyTJgwjIUEJQaQ9RVwSOFB4\njJHDBoauKY5y9txFVr27lVXvbOVwUWnAOl26JDItawQzp45i5rTR9NKAssgVi7gkUPhJKcMy+tsd\nh1yBwqJS1r2/k3W5OygtD/LMAo+HzOEDyZowjEkThpE1IUPjCCJtEHFJ4HBRKUOHKAlEgoaGBnbt\nLeb9TQV8mLeXo8dOBa3riYlh9IhBTBg7hHGj0hmbmcbA/j0120gkhIhLAkVHyhmS1tfuOKSdeb1e\niktO8mHeXjZ+vI9d+4ppqG9o8ZieKcmMHZVO5rCBjBg6gKFD+pM2sLduNy5iEXFJoPjoSdIH9bY7\nDulglRer2L67iK07D7N112EOFB4HrzfkcYmJCQwd3K9ZUkhP7c2g/j01C0miUsQlgZLjp0jVXSuj\nzrnzF9lZUEzBgSMU7D9CwYGSZjeyC8UTE8PA/j1JHdiL9EFGUujXpwd9+/Sgf98e9Ol1la4gJCJ1\nRBKYBzwDxAIvAE8GqPMs8FngIvA1YFuIY3sBrwNDgCLg74GKAJ/rPVZ6moH9e4YRpkSyhoYGiktO\nsvdgCYVFZRwqKuVwUellK5bD5YmJoU+v7kZi6H0V/fr2oGePbvRMSaZnj26kmK+ePbrRtUuixiLE\nNdo7CcQC+4EbgaPAZuBOYK+lznxgkfl3OvArIDvEsU8BJ82/jwA9gUcDfL+3tPxMs0cfRorc3Fxy\ncnLsDqPDdNbvO1NxgcJPyzj0SSnFJScoOXaKkmOnKDsR6JyibRIS4unZoxvdu3chuWsSyclJnDhW\nyISJU+jeLYnk5C4kd02ke3IXunVLomtSAkmJCSQlxdMlKYEuSQkkJsa7KpHo36d7tTYJhOo0nQYc\nwjhbB3gNuJXmSWABsMzczgNSgAHA0BaOXQDMMfcvA3IJnASIiYnMS/ZI/kcInff7eqYkMyUl+bLb\njVdX13L0+GlKjp/iyNGTlJ08S/nJs5SfOMuJU+c4fSb8K4iamlrKTlQ0SyzFBzawv7imhaP8eDy+\nhJCUlEBSopEgkhLjiY+PIyEhjoT4OOLjYklIiCMuLpYEc3+8uR3vVyc+Lpa4uFhiYmKIjTVf5nac\nWY6JMffFNX/fqBNr1vEQGxNDXFwsHo/RiOjfZ/QIlQRSAet9AEowzvZD1UkFBrVwbH+gzNwuM8sB\nxcS45+xJnCMxMZ5hGf2DrjGpqanjxOlzZlI4y4mT56g4d5GKsxc4XVFJxVnjdbriQsCVz63m9XLp\nUjWXLlVf+Wd1guKDG3hnSx0xMR48Hg8xHuOvx4Nvn3V/0z7jxC3GvOqJifHgiYnxJRf/+gCNF0hN\n5eb7G3lo/n6gY5v2t/zZmzduoezi78zPDVynpe8I9r0Eib0lnjBP2jvqSjJUEgg9PcMQTnSeIJ/n\nbel7YiP0SkDslZAQR+qAXiEnHXi9Xi5V1XDmbCWVlVWcv3CJC5VVLH3+JAv+7hYuXLjEhYtVXLhQ\nxfnKS1RWVnPxUjWXqmqorq7lUlUNl6prqK2p66Rf1k684G1oIMSsXdcqK69g554iu8NwhWxgjaX8\nGEYfvtVvgS9byvswzuxbOnYfRpcRwECzHMghmpKEXnrppZdeoV+HaEdxQCGQASQA24ExfnXmA6vM\n7WzgozCObRwQBmMs4In2DFpERNrPZzFm+RzCOJsHeMB8NXrOfH8HMCnEsWBMEX0HOAC8jTGYLCIi\nIiIi0W4exjjBQS4fg3C7dGA9sAfYDXzD3nA6TCzGosGVdgfSzlKANzCmOhdgdIFGkscw/m3uAv4H\nSLQ3nCv2EsYMxF2Wfb2AdURGT0Sg3/dzjH+fO4A/Az1siOuKxGJ0H2UA8QQeh3CzAcBEczsZo7ss\nkn5fo28DrwAr7A6knS0D7jW343Dhf2AtyAAO09Twvw7cbVs07WM2kEXzRvIp4Hvm9iO4e0wy0O+7\nCWicVvkELvx9M2g+q+hRgiwkixDLgRvsDqKdpWGM+cwlsq4EemA0kpGqF8ZJSU+MBLcSY8W/22XQ\nvJFsnMEIxklZsNmJbpFB899ndRvwx5YOduIk/GCLzyJRBkYWz7M5jvb2S+C7QKTNMh8KnAB+D2wF\n/gvoamtE7es08DRQDBzDuJ/XO7ZG1DHCXqwaAe6lafZmQE5MAl67A+gkyRh9y98ELtgcS3v6HFCO\nMR4Qacu94zBmv/3G/FtJZF2lDgf+GePkZBDGv9F/sDOgTtA4tz4S/RtQgzG2E5QTk8BRjMHTRukY\nVwORJB54E+MybbnNsbS3mRj3hvoEeBW4HvhvWyNqPyXma7NZfoPmU6LdbgrwN+AUUIcxqDjT1og6\nRhnNF6uW2xhLR/kaxhouVybxcBaouZkHo1H8pd2BdII5RNaYAMAHQKa5vZjAt1Z3q2swZqx1wfh3\nugx4yNaI2kcGlw8MR9Ji1Qya/755GDO8+tgSTTsJtsgsElyL0Ve+HaPLZBvG/2mRaA6RNzvoGowr\nAddOvwvhezRNEV2GcdXqZq9ijG/UYIw13kNkLVb1/333Ykyt/5Sm9uU3tkUnIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIiIiKd4/8BBwwzGzr9F3gAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 3 }, { "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 2: Comparing posteriors\n", "\n", "The next step is to compute the posteriors for the two teams:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "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())" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "posterior mean Germany 1.16529110336\n", "posterior mean Argentina 0.677218426967\n" ] } ], "prompt_number": 4 }, { "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. Germany's posterior mean is 1.2; Argentina's is 0.7. We can plot the posteriors:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "thinkplot.Pdf(germany)\n", "thinkplot.Pdf(argentina)\n", "thinkplot.Config(xlabel='goal-scoring rate', ylabel='probability')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW9+PHPzGQhK0kIZCEhYRME2bcACrFVL1LFulTU\nFrT2KpWLtr/WXrW21bS91XqrtV5bt1plsWKrqLhQsEpEQFlk39dAFhIg+0b23x/PmZkzQ2ZyJmQy\nmcP3/XpNc/Z5JtL55tm+DwghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIQyaBRwADgMPebjmOe38\nTmCcdmwYsF33qgQe8GtJhRBCdDsbcATIBEKBHcClbtfMBj7WtqcAX7XzHCtwCkj3SymFEEJ4ZfXj\nsyejAkUe0AQsB25wu2YOsFjb3gTEAUlu11wFHAXy/VVQIYQQnvkzUPTH9cu9QDvW0TVpbtfcBvy9\ny0snhBDCEH8GijaD11m83BcGXA/8s0tKJIQQwmchfnx2Ia79CumoGoO3a9K0Y3bXAl8DZ9p7g8GD\nB7cdPXr0wksqhBAXl6PAEKMX+7NGsRUYiurMDgPmAivdrlkJzNe2s4AKoER3/nbgTU9vcPToUdra\n2kz7euyxxwJeBvl88vkuxs9n5s/W1tYGMNiXL3N/1iiagUXAatQIqFeB/cAC7fxLqBFPs1Gd3rXA\n93X3R6E6su/xYxmFEEJ0wJ+BAmCV9tJ7yW1/kYd7a4HELi+REEIIn/iz6albtba28ddVB1j4fxvY\nfOB0oIvTJbKzswNdBL+SzxfczPz5zPzZOsN9xFGwadPa21i7s4hnV+wBID4mnL/9ZAZWa7B/PCGE\n6HoWiwV8+P73d9NTt2hobGHJJ4cd++XVDRwuqmRYWlwASyXExSUhIYHy8vJAF0PoxMfHU1ZWdsHP\nMUWg+OCrE5RVN7gc27T/jAQKIbpReXm5fUSN6CG0msMFM0UfxRd7Ss47tskk/RRCCBFopggUFTUN\n5x0rOFtL4dnaAJRGCCHMJegDRVtbG9X1TY79CUOdI2o3H2x3QrcQQggfBH2gqD3XTEuraheNCAth\n0rC+jnPHTlUHqlhCCGEaQR8o9LWJmMhQMpOiHfsnTkugEEI4LV++nClTphAdHU1SUhJZWVm88MIL\ngS5Wjxf0gaKqrtGxHRMRyoB+zkBReLaW5pbWQBRLCNHDPP300/z4xz/moYceoqSkhJKSEl588UU2\nbNhAY2Njxw/QaW5u9lMpe6agDxTVdc4aRWxkKFG9QunbuxcAzS1tFJ6tC1TRhBA9RGVlJY899hgv\nvPACN910E1FRUQCMHTuWZcuWERYWRkNDAw8++CAZGRkkJydz3333ce7cOQByc3NJS0vjqaeeIiUl\nhbvvvpucnBy+853vMG/ePGJjYxk9ejSHDx/miSeeICkpiYyMDD755BNHGV577TVGjBhBbGwsgwcP\n5uWXX3acsz//mWeeISkpidTUVF5//XUAtmzZQnJyssvQ4xUrVjB27Nhu+M0pQT+PQt/0FB0ZCkBG\nUjRnKtV/4JOna8jQNUcJIQLjpvn/6LJnrVhyq0/Xf/nllzQ0NHDDDe6LbDo9/PDDHD9+nJ07dxIS\nEsIdd9zBr3/9a373u98BUFJSQnl5OSdPnqSlpYUnn3ySDz/8kJUrV/L6669z9913c/XVV7NgwQKK\niop47bXXWLBgAceOHQMgKSmJjz76iIEDB7Ju3TquvfZaJk2axLhx4xzPr6qqoqioiDVr1nDLLbdw\n4403MmnSJPr06cPq1auZNWsWAEuXLuXOO+/szK+uU0xWowgDcGl+Onm6ptvLJIToWc6ePUtiYiJW\nq/Mrb9q0acTHxxMZGcm6det45ZVXeOaZZ4iLiyM6OppHHnmE5cuXO663Wq3k5OQQGhpKr16q1WLG\njBlcffXV2Gw2brnlFkpLS3n44Yex2WzMnTuXvLw8qqqqAJg9ezYDBw503HfNNdfwxRdfOJ4fGhrK\nr371K2w2G9deey3R0dEcPHgQgPnz57Ns2TIAysrKWLNmDXfccYd/f2k6wV+j0AWKmAitRtEvxnFM\nOrSFEH369OHs2bO0trY6gsXGjRsBSE9Pp6SkhLq6OiZMmOC4p62tjdZWZx9n3759CQsLc3luv379\nHNsREREkJiY6ZkNHREQAUFNTQ2xsLKtWrSInJ4fDhw/T2tpKXV0do0ePdimjPpBFRkZSU6P+0P3u\nd7/LyJEjqaur4x//+AczZswgKSmpS343RgR9oNB3Zsfqmp7sTpRIjUKInsDX5qKuNHXqVMLDw3nv\nvfe46aabzjufmJhIREQE+/btIyUlpd1nuKfD8CU9RkNDAzfffDPLli3jhhtuwGazceONNxpOeZKW\nlkZWVhYrVqxg2bJlLFy40PB7d4Wgb3qqaqfpqX9iJDYtc2xxeT31DRfXCAUhhKu4uDgee+wxFi5c\nyDvvvEN1dTWtra3s2LGD2tparFYr99xzDz/+8Y85c0ZN1C0sLGTNmjUen+lLXqvGxkYaGxsdzV+r\nVq3y+uz2zJ8/n9///vfs2bOn3WDnT0EfKNznUQCEhdjo3yfKcTz/jKTyEOJi97Of/YxnnnmGp556\niuTkZJKTk/nhD3/IU089xbRp0/j973/PkCFDyMrKonfv3lx99dUcOnTIcX97NYqOahn2/ZiYGJ57\n7jluvfVWEhISePPNN8/rWO+ohnLTTTdx8uRJbrzxRkcfSXcJ9gUb2u7/8wZH89Iff5jFoJRYAJ5Y\nvoOv9qvEgD+9ZRQzRrVfnRRCdA2LxSLZY/1s6NChvPTSS3zjG98wdL2n/ya+rkcR/DWKdjqzAVIS\nIh3bRaUyl0IIEdxWrFiBxWIxHCS6UtB3ZrsEisj2A0VxWX23lkkIIbpSdnY2Bw4cYOnSpQF5/6AP\nFE1aio6wECvhoTbH8ZQ++hqF9FEIIYJXbm5uQN8/6Jue7GIiw1w6g1J1NYpTUqMQQohOM02giNU1\nOwEkxIQTFqI+XlVdIzW60VFCCCGM83egmAUcAA4DD3m45jnt/E5gnO54HPA2sB/YB2R5eyN9RzaA\n1Woh2aVWIR3aQgjRGf4MFDbgeVSwGAHcDlzqds1sYAgwFLgX0CeG/xPwsXbPaFTA8CjGrUYBrh3a\nEiiEEKJz/BkoJgNHgDygCVgOuKdunAMs1rY3oWoRSUBv4Argb9q5ZqDS25vZZ2Xr6Tu0T8kQWSGE\n6BR/Bor+QL5uv0A71tE1acBA4AzwGrANeAWIxAv3pieQDm0hRHB44oknuOeeewJdDI/8OTzW6BRN\n99mBbahyjQcWAVuAZ4GHgV+533wgV40rfvdEPP1tN5Gdne04l5IQ4diWIbJCCFBzEnbt2kVxcfF5\n2WC7Q25uLvPmzSM/3/k38iOPPOL397yQIbb+DBSFQLpuPx1VY/B2TZp2zKJdu0U7/jYqUJxnePY8\nAG79xhCyZw5yOZciNQohhE5eXh6bN29mwIABrFy5kltuuaXd6/TpyM0gOzvb5Y/onJwcn+73529i\nK6qTOhMIA+YCK92uWQnM17azgAqgBChGNUldop27Ctjr7c2s7WQtSYgNJ9TmHCJbJ1lkhbioLVmy\nhKuuuop58+axePFix/G77rqL++67j9mzZxMdHU1ubi7btm1j3LhxxMbGcuuttzJ37lx++ctfOu75\n8MMPGTt2LPHx8UyfPp3du3c7zmVmZvL0008zZswY4uLiuO2222hoaKC2tpZrr72WoqIiYmJiiI2N\n5dSpUzz++OPMm6f+6M3Ly8NqtbJkyRIyMjLo27evY5U9gM2bNzN16lTi4+NJTU3l/vvvp6nJv8P/\n/VmjaEY1Ha1GjYB6FTVyaYF2/iXUqKbZqE7vWuD7uvvvB95ABZmjbufOY2knUtisVvrFR1B4VjU7\nlZTXMzA55rzrhBD+d8NjvqXV9ub9nGs6dd+SJUvIyclh8uTJ5OTkcObMGfr27QvAm2++yapVq5g6\ndSpVVVWMGjWKBx98kIULF7Jy5Upuu+02HnpIjfLfvn07P/jBD/jwww+ZOHEiS5cuZc6cORw6dIjQ\n0FAsFgv//Oc/Wb16NeHh4UyfPp3XX3+dBQsW8K9//Yvvfe97Lk1P7WWO3bBhA4cOHeLgwYNMnjyZ\nm2++mWHDhhESEsKf/vQnJk6cSH5+Ptdeey1/+ctf+NGPftSp34kR/q5brQKGoYbAPqEde0l72S3S\nzo9BdVzb7QQmacdvooNRT57SICbFO/spSsql+UmIi9X69espLCxkzpw5DB06lBEjRvDGG284zn/7\n299m6tSpAOzYsYOWlhbuv/9+xyJDkydPdlz78ssvs2DBAiZNmoTFYmH+/PmEh4fz1VdfOa554IEH\nSE5OJj4+nuuvv54dO3YA7a9j0d6xxx57jPDwcEaPHs2YMWMc948fP57JkydjtVrJyMjg3nvv5fPP\nP++aX5IHpmmEs3rI5Z7sEihkiKwQF6vFixdzzTXXEBOjWhW+853vuDQ/paWlObaLioro3991kGZ6\nurM79cSJEzz99NPEx8c7XgUFBRQVFTmuSU5OdmxHREQ4ljU1Sn9/ZGQktbWqZeTQoUNcd911pKSk\n0Lt3bx599FFKS0t9eravgj4poJ21vU4KICnOGSiKpUYhRMB0trmoK9TX1/OPf/yD1tZWx1KnDQ0N\nVFZWsmvXrvMWIUpJSaGwsNDlGSdPnmTIkCEADBgwgEcffZSf//znPpelvWYmX5ZVve+++5gwYQJv\nvfUWUVFRPPvss7zzzjs+l8MXpqlRePpFJyVI05MQF7v33nuPkJAQ9u/fz86dO9m5cyf79+/n8ssv\nZ8mSJeddP23aNGw2G88//zzNzc28//77bNmyxXH+nnvu4cUXX2Tz5s20tbVRW1vLRx99ZKjWkJSU\nRGlpKVVVVY5jviz4VFNTQ0xMDJGRkRw4cIAXXnih45sukGkChYcKBcnxziGyEiiEuDgtWbKEu+++\nm7S0NPr160e/fv1ISkpi0aJFvPHGG7S0tLj8sRkaGsqKFSt49dVXiY+P54033uC6665zzLuYMGEC\nr7zyCosWLSIhIYGhQ4eyZMkSj3+w6mssw4cP5/bbb2fQoEEkJCRw6tSp82o03moYf/jDH/j73/9O\nbGws9957L7fddptPNZLOCPqlUOf8ajUAC751KbMnp593QV1DM7f/7jMAQm1W3vrFN7CZaHy0ED2F\n2ZdCnTJlCgsXLuTOO+8MdFEMk6VQ3dg8VCkiw0PoHaX+CmhqaaW8urE7iyWECFLr1q2juLiY5uZm\nFi9ezJ49e5g1a1agixUQpunM9iYpPoLKWhUgTpXVkdi7V4BLJITo6Q4ePMitt95KbW0tgwcP5u23\n3yYpKSnQxQoI0wQKT6OeQAWKQwVqGkZJeT2jBnZXqYQQweqee+7p0Yn6upNpmp68NbalSIe2EEJ0\nmmkCRUc1CjuZSyGEEL4xTaDwNjosSWZnCyFEp5mnj8JLpEiWfE9C+F18fLzfx/ML38THx3fJcy6K\nQJEQG06IzUJzSxsVtY3UNzQTEW6ajy5Ej1BWVhboIgg/uSianmxWK0m6Dm3ppxBCCONMEyi8dWaD\npBsXQojOMk2g6KhpVNKNCyFE55gmUHjrowBJNy6EEJ1lmkDR0VgLSTcuhBCdY5pA0VEfhaQbF0KI\nzjFNoOho/LZ7Z3Zrq3nTIQshRFcyTaDooEJBZHgIsZHOdONl1Q3dUCohhAh+JgoUHc8ITU7Qd2jL\nyCchhDDCNIHCSOYAmUshhBC+83egmAUcAA4DD3m45jnt/E5gnO54HrAL2A5s7uiNjOSY0acbLy6T\nQCGEEEb4M+GRDXgeuAooBLYAK4H9umtmA0OAocAU4AUgSzvXBmQDhhLIdDTqCSTduBBCdIY/axST\ngSOomkETsBy4we2aOcBibXsTEAfo1xo0vvi3z01P0kchhBBG+DNQ9AfydfsF2jGj17QB/wa2Ah2u\nR2gxEFMk3bgQQvjOn01PRicqePqGvxwoAvoCn6D6Or5wv+hA7lIA/lz9GTfPmUV2drbHN5J040KI\ni1Fubi65ubmdvt+f35KFQLpuPx1VY/B2TZp2DFSQADgDvItqyjovUAzPngfAA/dOYWj/3l4LZLNa\n6RcXQVGpanYqqagnMynGwEcRQojglZ2d7fJHdE5Ojk/3+7PpaSuqkzoTCAPmojqz9VYC87XtLKAC\nKAEiAfs3eBRwDbDb25sZmUcBMkRWCCF85c8aRTOwCFiNGgH1KmrE0wLt/EvAx6iRT0eAWuD72rlk\nYIWujG8Aa7y9mdEVGFXOp1JAAoUQQhjh7wb6VdpL7yW3/UXt3HcMGOvLGxkZHguuHdoyRFYIITpm\nnpnZBkfS6tONF5fJEFkhhOiIaQKF1eAnkXTjQgjhG/MEik52Zku6cSGE8M40gcJoZ7akGxdCCN+Y\nJlAYrVGApBsXQghfmCZQGMkeaydzKYQQwjjTBAqDo2MByfkkhBC+ME2g8K1GoVuXQgKFEEJ4ZZpA\nYXTCHUiNQgghfGGaQOFDhUIChRBC+MA0gcKXUU/2dOMA5TUN1Dc0+6tYQggR9EwUKIxfa083bldS\nIbUKIYTwxDSBwpfObJAhskIIYZRpAoUvndkgOZ+EEMIo0wQK38KEa41ChsgKIYRnpgkUvnRmg1ug\nkHTjQgjhkWkCha9VChkiK4QQxpgmUNh87aNI0PVRVEi6cSGE8MQ0gcLHlifXdOPNrZTVSLpxIYRo\nj5FAsQL4lsFrA8bXPgpwSzcu/RRCCNEuI1/+LwDfBY4ATwLD/FqiTvJ1HgXIXAohhDDCSKD4BLgD\nGA/kAZ8CG4HvA6F+K5kPfO2fsJMObSGE6JjR5qQ+wF3AfwLbgOeACaggEnCdqEwAkm5cCCGMMBIo\n3gXWA5HA9cAcYDmwCIjp4N5ZwAHgMPCQh2ue087vBMa5nbMB24EPvL1JZ5qdQJqehBDCiBAD17wC\nfOx2LBxoQNUqPLEBzwNXAYXAFmAlsF93zWxgCDAUmILqD8nSnf8RsI8OAlJnOrJBmp6EEMIIIzWK\n/2nn2JcG7puM6gDPA5pQtZAb3K6ZAyzWtjcBcUCStp+GCiR/pYPpdJ0NFH3c0o2fa5R040II4c5b\noEhB1RgiUB3ZE7Sf2ahmqI70B/J1+wXaMaPX/BH4GdDa0RtZOzlw1z3duPRTCCHE+bw1Pf0HcCfq\ni/tp3fFq4OcGnm10qrN7dcACXAecRvVPZHu7+UDuUsJCrDzesI7s7Gyys71efp7khEiKStUciqLS\nOjKTOup2EUKI4JKbm0tubm6n7/cWKF7XXjcD73Ti2YVAum4/HVVj8HZNmnbsZlSz1GygFxALLAHm\nu7/J8Ox5xESG8vhDV3aiiJCWGMXWg2ewWKDwbG2nniGEED2Z+x/ROTk5Pt3vLVDMA5YCmcBPdMct\nqNrCMx08eyuqkzoTKALmAre7XbMSNXpqOaoTuwIoRtVY7LWWmcCDtBMknAXqXB/FgcNnWb/uOMcL\nKrFY4K0PDzBhQDyDMuM79TwhhDAjb6379n6IGA+vjjSjgsBq1Milt1AjnhZoL1CjqY6hOr1fAhZ6\neJbXZqzOTLjbuqOIR3+7llMFleoN2iCvuJpH/2ct+w6e8fl5QghhVp2cqtZjtM351WriY8J5/cGZ\nhm+qqm7gxz9fTUXlOZpp45h23AoMBnqFh/KbR7IZMijBH2UWQoiA0uaeGf7+99b09H9ezrUBDxh9\nE3/zdXjsy4u3UVF5DoDE2AhC+0ZQe66ZU8XVtLS00dDQzHMvb+YPv7masFCbP4oshBBBw1vT09eo\nfoavPbx6DF9ank4WVLJxs3NE7sL/nMig1FjCw2ykpsTQpgWGgqIqVnyw39NjhBDiotHRqKeg4EuN\n4tN1xx3bUyb0Z+LYVNYfK+NQQSVhoTamTBvA1s/zAHjngwPMmJZBarIMmRVCXLy81Sj+pP38oJ3X\nSj+XyycWg1WKpuYWctefcOxfc+VgANISnfMH+6bEMnxoIgAtLa388/19XVhSIYQIPt5qFEu0n0+3\nc65HrRtqtD6xdfspqrWV7BL7RDJ6ZD8AUhOjHNcUltYyb+5oHv3tZwCs23iSW+ZcSv+U2C4tsxBC\nBIuO+igAclG5ncqBUtRaFJ/7t1i+sRqsUeSuz3NsX3l5Jjab+vhp+kBxto5LL0lkzGXJALS1tfHP\n96SvQghx8TKSJelbqHkOz6GywR5FzZjuMYzEiebmVnbvP+3YnzFtgGM7OSHCMRfjTGU9DY0tzL1x\nhOP8+k0nOVsqS6UKIS5ORgLFM8CVqBnSM1G5l/7oxzL5zMh6FEfzyjl3TmWH7dsnyqWDOizE5lib\noq0NisrqGD40kZHD+wLQ2trG6s+O+qHkQgjR8xkJFFWoGoXdMe1Yj2Fk0NPufSWO7ctG9D0vuPTv\n42x+Kjijcj7Nvnqo49iatcdobGy5wJIKIUTw8RYobtZeW1GpNu7SXh9qx3oMI8Nj9+xzpuUYNaLf\neef799UFCi054OTxqST2USOiqmsaWP/VyQstqhBCBB1vgeJ6VLrvXqiU3/ampzPasR6jo0DR2NTC\ngcNnHfuXDT8/UKS5jXwCsNmszPrmEMfxT3KPnXefEEKYnbfhsXd1VyEuVEcVikNHSmlsUs1GKUkx\njlqCnn4uRdFZZ8f1N2dk8uY7e2hpaeXgkVLyC6tI7y9DZYUQFw8jfRQRqCywfwFeA/6mvXqMjobH\nutQmRvRt9xqXPoqztbS2qqkivWN7MWlcquPcZ7qZ3UIIcTEwEiiWotaxnoWaU5EO1PixTD7rqIfi\n+IkKx/ZQDxlhY6PCiI0MA6ChqYXSqnOOc9+cOdCxnbvhBM3NHa7OKoQQpmEkUAwBfokKDotRcyim\n+LNQvuqoRpF3stKxnTkgzuN1+uanQt28ibGXJZGgDZ+trDrH9l3FnS2qEEIEHSOBolH7WQmMAuKA\n9ttvAsTbPIq6+iZOlagKkNVqYUD/3h6v1afyKNAti2qzWZk5PcOxv+7LEwghxMXCSKB4BUgAfoFK\nBrgPeMqfhfKVtwrFifxK7Kmp0lJjCQvzvL6EfuRT/mnX1rUZU50zubdsP0VdfVPnCiuEEEHGaKAo\nQ+V3GoiqTbzoz0L5yluNIu+ks38iM91zsxNAum4uxckztS7nMtLjGJCmaiONjc1s3lbYmaIKIUTQ\nMRIoElGr3W0HtqHSj/fxZ6F85a1GkZevCxQZ3gNFRr9ox/bJ0zW0tbkmydXXKtZ/mY8QQlwMjASK\n5agJdzcBt6Am3L3lz0L5ytuEu7wTzo7sgR0EisTevYjqpaaW1NQ3UVrd4HJ+epYzUOzYU+xYTlUI\nIczMSKBIBn4DHEflefotarhsj+Fp4aKWllZO6GsU6Z47skE1YQ1wq1XoJfWNYvglalGj1tY2lyVV\nhRDCrIwEijXA7dq1VmCudqzH8FSfKD5d45iRHR8XQe/YjjOPuASKkvOni1yhq1Ws/0oChRDC/LwF\nihqgGrgHeAM1TLYReBO41/9FM85T09Mp3Rd9Wqqxda/1/RQnTp8fKKZPSXcseHTg8FmK27lGCCHM\nxFugiAZitJcVlRcqRNs29q2rZnMfAA4DD3m45jnt/E5gnHasF7AJ2IEajvuEtzexevgUp4qdX+LJ\nSdHtX+TGW9MTQGxMOGNGOlveJKOsEMLsjDQ9AdyAWjv7D6isskbYUCvizQJGoJqvLnW7ZjZq5vdQ\nVC3lBe34OdRiSWOB0dr25Z7eyNPwWP1f+ylJvtco8s/U0tJ6frqOK3SjnzZsKjD0XCGECFZGAsWT\nwAPAXmC/tu31L3zNZNSCR3lAE2r01A1u18xBpQUBVYOIw9lRbs+hEYYKOmWe3shT01NRsT5QGKtR\nxEaFER8TDqicTyXl9eddM2l8KmGhauLeifwKCop61DpOQgjRpYyumX0NKmPsq6gawnUG7usP6Ht7\nC7RjHV2Tpm3bUE1PJcBaVBNUuzyNji0uqXZsGw0U4FqryGunQzsyIpTxY1Ic+xs2Sae2EMK8vK1H\nYdeG+ku/VNuPw54To+P7jHD/mrff14JqeuoNrEat1Z3rfvOB3KXU74umbl8S2dnZZGdnA9DU3MJp\nx7oSFpL7GQ8Umckx7DiqPm5ecTXTRpw/Gnja5DS+2qqanTZuLmDujSMNP18IIbpTbm4uubm5nb7f\nSKB4AjUjey3qS30m8LCB+wpRKcnt0lE1Bm/XpGnH9CqBj4CJtBMohmfPY8aoFH56yyiX46fP1Dpm\nVif2ifCa48ndQF1/xvHi6navmTg2lbCwEBobm8kvrORkQaUjxYcQQvQk+j+iAXJycny6v6OmJyvQ\nCkwF3gXe0baXG3j2VlQndSaqn2EuKqmg3kpgvradBVSgmpoSUTUXUAsnXY1KIdJ+Idv5FJ3pn7Ab\nmNxxoOjVK4QJuuYnmXwnhDCrjgJFK/DfQBHwPuqL/ZTBZzejVsZbjepfeAvVGb5AewF8jJrtfQR4\nCVioHU8BPkP1UWwCPgA+9fRGlnam3BWX+D7iyS6tbyShIepXc6byHNV1je1ed3mWszK0YVPBebmh\nhBDCDIw0PX0CPIj6otenVPU4CklnlfbSe8ltf1E79+0Gxht4PtB+UkD9ZDujcyjsbFYrGf2iOaKN\nZjpeXM3oQefnQRw/Opnw8BAaGpopPFXFifxKrwsjCSFEMDIy6uk24L+AdcDXuleP0V6upyJdk1Gq\nj4ECVIe2nafmp/DwEJf1tGX0kxDCjIwEikuBP6NmTm9HpRwf4c9C+aq9eRSndSvU+VqjABjkEig8\np+mYPkXX/LQ5X5qfhBCmYyRQLEEFiz+hZlqP0I71GO4Vira2NkpLnRPlEvtE4it9h/axYs8T6saN\nSiaiVyig+kWOnajweK0QQgQjI4FiJPAD1PDYz4D/1I71GO4pPCqrGmhqVlljIyNCiYwI9fmZ+qan\ngjO1NGhZaN2FhdmYPF6an4QQ5mUkUGxDDYm1y6Kn9VG41SjOltU5tjtTmwCIDA9xrKHd0tpGnod+\nCoBpuuanjZuk+UkIYS5GAsVEYANwApW3aaN2bDewy28l84F7H8XZUl2gSOhcoAAYnBrr2D7iJZ/T\nmMuSHLUaHZihAAAcCElEQVSW02drOXLMyIAwIYQIDkaGx87yeykukNWtk+Js2YX1T9gNSY3l811q\n2oi3QBEWamPKxP6s/SIPUHMqhg7uUcuKCyFEpxmpUeR18Ao49zFPZ8/qm54iOv3cIQZrFADTJ+ua\nn7ZI85MQwjyMrkfRo7nPo+iKPgpQI5/szVoFZ2qpb2j2eO3okUlER4Wp9y+t4+CRUo/XCiFEMDFF\noHAfHttVfRQR4SGk9VUd2q1tbRzz0qEdEmJlygRnFvWNm2VBIyGEOZgkULhGijOlXVOjANfmp6Md\nND/pcz9t3JxPa6s0Pwkhgp/pAkVzcyvlFee0PQt9EjrfRwGugeJQYaXXay+7tB8x0Wp1vLLyeg4c\nPntB7y2EED2BKQKFvkJRWl6Pfe2j+N69CA0xvg5Fe4b2d64xcTDfe6Cw2axMnZTm2JfU40IIMzBJ\noHBGitIu6si2G5gcQ5iWcvx0RT1l1Q1er5822RkovtxcSEtL6wWXQQghAskUgULfmX2mi4bG2oWG\nWF2anw4WeM/lNHJ4X3rH9gKgvLKe/Yek+UkIEdzMESis/qtRAAxLd64x4Wvzk+R+EkIEO1MECv2Y\nJ2dHNiTEXXiNAmBYur6fouPssC7NT1sKaG6W5ichRPAyR6DQ1SjKdIEirnevLnn+sDRnjeJIURXN\nHfQ7XHpJIvFakKqqbmDn3pIuKYcQQgSCKQKFfnhseaUzz1N8XNcEioSYcPppX/yNza0cO+V54h2o\n5qcZUwc49j/fcKJLyiGEEIFgkkDh3K6odNYo4ruo6Qlcm58OGGh+mjk9w7G9+etC6uqbuqwsQgjR\nnUwRKPTDY/V9FF1VowC4LCPesb0nr7zD6zMHxJGhdYI3NrXwpaT0EEIEKVMFirr6Jhq0xH2hITai\nIn1f2c6TEbpAsf9khaH0HDOmOZuf1n0pzU9CiOBkikBhb3oqr3Dtn3BfIvVCpPeNIjZSZYetqmsk\n/0xth/eofgpVht37zrgkKxRCiGDRHYFiFnAAOAw85OGa57TzO4Fx2rF01Drde4E9wAOe3sA+j6Ki\n0jlruiubnUDVWkZmOmsV+0523PzUJyGS0SP7aXttrPvyZJeWSQghuoO/A4UNeB4VLEYAtwOXul0z\nGxgCDAXuBV7QjjcB/w8YiVqn+7/auRdwjnrS1yi6amis3ogBzmGye090HCjAtVP78w0nZEEjIUTQ\n8XegmAwcQa2E1wQsB25wu2YOsFjb3gTEAUlAMbBDO14D7AdSvb2Za0d21414srss07VD28iXftbE\n/oSFqRVn8wsrOXai4xFTQgjRk/g7UPQH9DksCrRjHV2T5nZNJqpJalN7b2LTmp7KK/0z4skuIyma\nqF6qg7y8uoGCsx33U0T0CiVLt6DRZ+uOd3m5hBDCn0L8/Hyj7Szuvc76+6KBt4EfoWoWLg7kLmVJ\ncV/Wp8RytqofkAj4J1DYrFZGDYznq/2nAdh5rIz0vtEd3vfNmQMdo54+33CC+XNHEx7u71+9EEIo\nubm55Obmdvp+f39bFaI6pe3SUTUGb9ekaccAQoF3gGXAe+29wfDseXz/xsu4cmwqj//+c4rOqnQZ\n/uijABgzqI8jUOw6Vsp1UwZ0cIfKKJvUL5qS0zXU1Tfx5dYCsqdn+qV8QgjhLjs7m+zsbMd+Tk6O\nT/f7u+lpK6qTOhMIA+YCK92uWQnM17azgAqgBFXLeBXYBzzr9V20+oi/ZmXrjRmU4Njek1dOS2vH\nCf+sVgtXzxzo2P9krTQ/CSGCh78DRTOwCFiN+sJ/C9UpvUB7AXwMHEN1er8ELNSOTwe+B1wJbNde\ns9p7E+eoJ//2UQCk9omkr1ZbqT3XzJFC7+to22VfkekYxrv/0BkKTxm7TwghAq075lGsAoahhsA+\noR17SXvZLdLOjwG2acfWa+Ubi+rIHgf8q703sFosNDW3UF1jn0dhoXdMeFd+BgeLxcKYQX0c+zuO\nlhq6LyEugoljnYO2/p0rtQohRHAwx8xsK1TqJtvF9Q7HZvPfRxsz2Nn8tN1goAC4KtvZ/LR2fR5N\nzS1dWi4hhPAHUwQKi8VCRZWz2cm+FKm/jBvcx9HcdTC/kuq6RmP3jUomId65TsWWbUV+K6MQQnQV\nkwQKqKpy1ih6x/qn2ckuJjLMkXa8ta2Nrw8bq1XYbFa+OcNZq1jz2TG/lE8IIbqSKQKFzWKhqsYZ\nKGL91D+hN2FoomP768NnDd931cyBjmSFu/aVcLLA+xrcQggRaKYIFBaLxbVG0c2BYvvRs4aGyQL0\nTYwia6JzpvZHaw53edmEEKIrmSRQQGW1rkbh56YngIHJMSRoAam6romD+cZrBt+6Zqhj+/MNJ3Sj\ntYQQoucxSaBwrVF0R9OTxWJh4iV9Hfv22dpGXHpJIoO0hZAam1r4ZK30VQghei5TBAqb1eJao+iG\nQAGQdWk/x/ZX+08bTiFusVhcahWrPj1Kc7OxpishhOhupggUFrp31JPd6IEJRPVS6bJKKuo5dqra\n8L3Ts9Idw3hLy+rY9HVhB3cIIURgmCNQBKhGERpiZdKwzjU/hYXauObKQY79D1cf6tKyCSFEVzFH\noACqq7u/RgGuzU8b95X4tILdf3xzsGMG+cEjpew9cKbLyyeEEBfKFIGipbWVuvomQGVqjYoM67b3\nHj84kfBQGwAFZ2vJKzlvyQyPEuIimDnNuVTq2yv3d3n5hBDiQpkiUNTXNTm2Y6LDHVlau0N4mI3J\nw53NT5/vOuXT/TddP9wxAW/nnmIO+5A7SgghuoMpAkWdLlB0V/+EXvboFMf2F3uKaW013vyUmhzD\n5VnOdZukViGE6GnMESjq9YGi+5qd7MYO7kOs1tx1tvIce0+U+3T/zddf6tjesr2I4ycqurR8Qghx\nIUwRKGpqndlbY2P8mzm2PSE2K9NHJjn2fW1+GpDWm6yJaY79dz6QWoUQoucwRaDQNz1154gnvZmj\nkx3b6/eUUN/Q7NP9N89x1io2bi4g3+DKeUII4W+mCBS1tYHtowAYnh5HWmIUAPWNzWzYW+LT/YMz\n4xnv6Oto4+9v7+7iEgohROeYIlC4Nj11fx8FqLQcV493ZoX9ZJvvM61vu3mkY3vT14XsOyjzKoQQ\ngWeKQFGrCxT+Xt3Om+wxKYTY1FDXA/kVnPBhTgXAkIEJXDF1gGN/yfJdPk3gE0IIfzBFoHCtUQSm\n6QkgLjqcycOcM7U/3nzS52fcccsoQkLUf5ZDR0v5aqvkgBJCBJYpAoW+jyJQndl2syc750Ss3XmK\nat3QXSOS+kYx+ypnZtml/9hFU3NLl5VPCCF8ZYpAUR2AhICeXJYZT2ZSDAANTS38uxN9FTfPGe5I\nQ1JcUsNqWVtbCBFA3REoZgEHgMPAQx6ueU47vxMYpzv+N6AE8DwEqA3q6u1NTxZiogPTmW1nsVi4\nLsvZz/Dx5nzDy6TaxUSHuwyXXf7OHsoq6rusjEII4Qt/Bwob8DwqWIwAbgcudbtmNjAEGArcC7yg\nO/eadq9H+i/hmOgwRzbWQJoxKtkxU/t0RT3rdhf7/IzZVw0hRauZ1NU38bdlO7q0jEIIYZS/v1Un\nA0eAPKAJWA7c4HbNHGCxtr0JiAPss9e+ALzmw2hpcY4KCnSzk114qI3rdbWKd7447lP+J4CwMBsL\n7hrv2N+4OZ+vd/o241sIIbqCvwNFfyBft1+gHfP1Go9aWnteoAC4dnI6keFq9bv8M7V8dcD4okZ2\no0cmMXO6Mw35y69v49w532Z8CyHEhQrx8/ON/hntnhfc8J/fh79YRlhtI1Ys9O09E7jScOH8KSYi\nlFmT0lixPg+At3KPkTW8n88p0O+6fQzbdhZTXdPAmdJa3np3L3fePsYPJRZCmFVubi65ubmdvt/f\ngaIQSNftp6NqDN6uSdOOGZIx5XZiS+sIwcK48YM7XVB/uGFqBh9vzudcYwt5JdV8saeYmbqU5Eb0\nju3FnbeN5vm/bgFg5b8OMWl8KiN0S7AKIYQ32dnZZGdnO/ZzcnJ8ut/fTU9bUZ3UmUAYMBdY6XbN\nSmC+tp0FVKBGOhnS0tLqqI70pKYnUBPw5mQ5m47+vvYoTc2+jYACuPKKTC7Tllxta2vj2Rc3UVPT\n2MFdQgjRNfwdKJqBRcBqYB/wFrAfWKC9AD4GjqE6vV8CFurufxPYCFyC6sf4vvsbuPRRBHiyXXu+\nPT2DmMhQAIrL6li1Jb+DO85nsVh44N7JREdpa16U1vHi619Leg8hRLfojrGkq4BhqCGwT2jHXtJe\ndou082OAbbrjtwOpQDiqeeo194e3tjj/Qu/dw2oUAFG9Qrnl8oGO/eW5R6moafByR/sS+0Ry390T\nHfsbN+fz2bq8riiiEEJ4FfhJBxfIZXhsD6xRAHxrygD6aynIa881s/TTI516ztRJaVx95SDH/itL\nt3OyoLJLyiiEEJ4Ef6BobXP2UUT3zEARGmLlB7OGOfb/va2QPXllnXrW3XeMJS01FoDGxmaeeHYD\nVdW+11CEEMKo4A8U+qanHlqjAJgwNJEpw50jlf68ch8NTb4n+wsPD+Gn/5VFuDZHo+R0DU89t1ES\nBwoh/Cb4A4WuMzsmQIsWGXXvty4lqpf6gi8qreONTjZBZaTH8eMfTsE+/WTfwTO8/Po26dwWQvhF\n0AcK2tTXZWREKKEhtkCXxqvE2F7cefUljv33vzzB9iNnO/WsKRP6871bRzn2P113nHc/PHDBZRRC\nCHfBHyg0PbUj2901E/ozYWiiY/9P7+6lvBOjoABu/NYwlxQfy/65mw9XH7rgMgohhJ5pAkVPHBrb\nHovFwv3fHkm81vFeXtPAU2/t7NREPIvFwn13T2Skru/jb2/s4F+dbNISQoj2mCJQWLD0uFnZ3sRH\nh/OjG0di0YZr7TtZwav/OtipPoawUBs//8nlDL/EWUt5efE21qyVxY6EEF0j6ANFT03f0ZFxQxKZ\np1vydNWWfN7beKJTz4roFcovfnoFQwf1cRx78bWt/PP9fdLBLYS4YEEfKOyCpY9C76bpmVwxKtmx\n//qaQ6zdWdSpZ0VGhPLLn13BoMx4x7E339nDX17dSnMnmrWEEMLONIEivnevQBfBZxaLhftvGMnI\nDOeX+3Pv7uXzXZ1boCg6Koych2cyekSS49in647zuz+ulySCQohOM0+giIsIdBE6JTzUxiO3jyVT\nW/a0ta2NZ1fsYc3X7tnYjYmKDOPRBy8n+/JMx7Edu4v56S8/4eCR0q4oshDiIhP0gcLeRxEfF3w1\nCruYiFB+fecEMpKiARUs/rxyH8tzj/q8hCpAaIiN+++ZxK3fHuk4dqa0ll/8z1re/ehAp54phLh4\nBX2gsAvGpie93lFh/Hr+BAanxDqOvbn2KE8s30FNfZPPz7NYLNx200h+dv80oiLVjPWWllaWvrWL\nXz2RS0FRVZeVXQhhbkEfKJw1iuBsetKLiw7nt9+fyJhBCY5jmw+e4acvb+LYqc59sU+dlMYffnM1\nlwx2jojad/AMP/nFGpav2Etjo+SIEkJ459sCzj1P2+h5yxkVEcayl24MdFm6THNLK0v/fdhluGxo\niJVbZwzi29MzCOtEqpLm5laWr9jLex+7Nj0l9YvmthtHcsXUAT6v5y2ECE4WNYnL8P/hg/2boW30\nvOXMSOnN//1+VqDL0uU27C3h+ff3UtfQ7DiW3jeKH153KZdlJni507O8kxW88LevOXzMtWN7QFpv\nvnvLKCaMTZGAIYTJXZSBYs7wfvzm51cGuix+UXi2lmfe2c0Rtz6FqSOSmDtzEAOTY3x+ZktLK6s/\nO8ryFXupqXUdNts/JZbr/mMo2dMzHKnMhRDmctEFijHzljM/K4OfLMwKdFn8pqW1lVVbClj27yPU\nNza7nJs2IombrxjIkNRYD3d7VlvXyPsfH+KD1YdoaHB9bkx0ONmXZ5A9PZOBGXEXVH4hRM9yUQaK\nH80axvfvGBvosvjd2cpz/G31QTbsLTnv3JDUWK6ZkMaMUclE+FgTKK84x3sfH+DfucepP3f+CKvM\nAXFcMXUAk8en0j/F94AkhOhZLspAkXPbWG6YPazjq03i2KkqluceY9OB0+edCw+1MWFoItNGJDHh\nkkQifQgadfVNfPr5cT5ac5jTZ2vbvaZ/SiwTx6YwemQSl16SSK9e0jwlRLC56ALF2HnL+dMPs5gx\nLaPjq03m2Kkq3v/yBBv2lrSbpjzUZmVYem9GD0xg9KAEBqfGGhox1dLSyu59p1m7Po9NXxfR6Nbc\nZWezWRk6KIGhgxMYOiiBwZnxJCdF2/8RCiF6qJ4WKGYBzwI24K/A79u55jngWqAOuAvY7sO9bWPn\nLWfJw1cyakS/ri15EKmua2TtzlOs+bqA/DPt1wQAQmwWMvrFMDAlhsEpsQxOjSG9b7TXWkddfROb\nthayeVsh23eXeAwadlGRYQzOjCe9fyypKTGkpsSQlhJDQnyEBBAheoieFChswEHgKqAQ2ALcDuzX\nXTMbWKT9nAL8CcgyeC9A27h5y1n55LWkdaIzt6fLzc0lOzvbp3tOnq7hy30lfLn/NMeLqw3dExcV\nRkqfSJLjI0lJiCCxdy/iosOJjw4jPiac2MhQbFYrjY0t7Nxbwq69Jezed5qTBZWGyxUeHkLfPpEk\nxEfQt08kfRIiyD+xi29+40ri4yKIiQ4jNiacsLCevZytLzrz3y+YmPnzmfmzge+Bwp8NzJOBI0Ce\ntr8cuAHXL/s5wGJtexMQByQDAw3c6xDMeZ686cw/1gH9ohnQL5q52YMpq25g9/Eydh0vY09eOcVl\nde3eU1HbSEVtI/tPVrR73ma1EBsVRnx0GFG9QokMD2HQpP4MmdSfqspzVFaco7S0jpKSGs6da8KK\n+hdoAcd2U0MTdUWVnCxSwcUC7N+5gi+3u/4TDAsLISY6TAWO6HAiI0MJD7cRGRFKeHgIEb1C6BUe\nQkSE+hkeFkJoqJWQECuhoTZCQ7RtbT9E2w/TtrtzjojZv2zM/PnM/Nk6w5+Boj+Qr9svQNUaOrqm\nP5Bq4F4ArFYLkRGhF1xYM0qICWfm6BRmjk4BoKa+iePF1RwtquLoqWqOF1dRXFZPU4v39SpaWtso\nr26gvNrL2t5WsKZEE9bcSkNjM01NrTQ2tdDU2EpjcwutLecnIiwDDqOO27++LY1NUNaEpazWeczt\nPk9f9e0dP+9ei1oR0WJRO1b7MXVCnbM6r7Ffp35oxzw83KI/YIGjmw7wZen7nq+zuJ5ybZlzPen5\nuu5nL//BLw/xddWHgS2Mn3T1Zwv0f7ML5c9AYTRF6QX9CsPCbNL2bVB0RCijBiYwaqBzVndLaytl\nVQ0UldZRXF7PqbI6ymtUUCivaaSippGqOuNrWai/4MPOO97a2kZzc6t6taifFWE2ekWE0NLSRktL\nG62trfg9sW2b9j9tjh2/qW5oprDynF/fI5Aq6ho57mF0XLAz82frabKAf+n2HwEecrvmReA23f4B\nIMngvaCap9rkJS95yUtePr2O0EOEAEeBTCAM2AFc6nbNbOBjbTsL+MqHe4UQQpjAtajRS0dQtQKA\nBdrL7nnt/E5gfAf3CiGEEEIIIUTXmIXq0zhM+/0XwSwdWAvsBfYADwS2OH5hQ02u/CDQBfGDOOBt\n1HDufahmVTN5BPVvczfwdyA8sMW5YH8DSlCfxy4B+AQ4BKxB/TcNVu19vv9F/fvcCawAegegXH5n\nQzVJZQKhmK8PIxmwZzmMRjXBmenzAfwEeANYGeiC+MFi4G5tOwRz/Z8wEziGMzi8BdwZsNJ0jSuA\ncbh+kT4F/Le2/RDwZHcXqgu19/muxrnC6ZME9+fzaCquo6Ie1l5m9R7wzUAXogulAf8GrsR8NYre\nqC9Ss0pA/eESjwqCH6AyKAS7TFy/SO0jMEH94XaguwvUxTJx/Xx6NwLLvN0crGtme5qoZ0aZqL8G\nNgW4HF3pj8DPAO8z/YLTQOAM8BqwDXgFiAxoibpWGfA0cBIoAipQQd9sklDNNWg/k7xcG+zuxjn6\ntF3BGijaAl2AbhKNauv+EVAT4LJ0leuA06j+CTPOlAxBjd77i/azFnPVdgcDP0b9AZOK+jf63UAW\nqBvY5x6Y0aNAI6qvyaNgDRSFqA5fu3RUrcJMQoF3UFXC9wJclq40DZXj6zjwJvANYElAS9S1CrTX\nFm3/bVyHfQe7icBGoBRoRnWETgtoifyjBNXkBJCC+uPGbO5CzWUzbaA3+4Q8C+rL84+BLoifzcR8\nfRQA64BLtO3HaT9FfrAagxqJF4H6d7oY+K+AlqhrZHJ+Z7Z9NOXDBH9nbyaun28WauRaYkBK043M\nPCHvclT7/Q5UE8121H9Ys5mJOUc9jUHVKMw69PC/cQ6PXYyq/QazN1H9LY2ovs/vozrt/405hse6\nf767UdMKTuD8fvlLwEonhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQojs9Dvy0m94rFfhnN72X\nXQZweze/pzCZYE3hIURX6a4cPiGoSU/f8dOzPRkI3OGH9xQXEQkUIlj9EpX6+QtUQjN7rWAsau11\n+6xo+4zae4DNqNnub6NSUHgTBXykXb8buFU7PgnYoB3fpF3XC5UtdhcqY2y2du1dqJnnn6IWwclA\npb+wn1sBrELN/tWn+fgBKuvAJlT22f9rp3yPA0uB9ajZ0Rmo1CFfa6+p2nVPotYj2I5KLmlFLVqz\nGfU7ureD34MQQgSlSagvvjBU9tJDqIWQQH1ZX6Ft5+DMl5Wgu/83wCJt+zHab3q6GXhZtx+rvd9R\nYIJ2LBq1iNZPgb9qx4ahUiOEo4JBPs5glYkz385d2rNitGvzUKnyU1EJE+NQNYV1wHPtlO9xVJoQ\n+wJCEbrtoTiTErrn07oXlTEU7fotWrmE8EhqFCIYTUdl1G1EpV+3fxHGovIqfaHtLwZmaNujtOO7\nUNkyR3TwHrtQq4A9icq9VYUKAqdQf7GjvXeLVh77wi8HUYHiElSz1ieoNRva8ylQDTSglkzNBCYD\nn2v3NKP6NNpLx96Gqq00aPthqGC1C/gHziSZ7vdeA8xHBdqvUAF0iIfyCQF4b9sUoqdqw/UL0NO6\nFvrjr6PSm+9GLd2Z7XZtGvCh9uwXULWJccC3gN+ivtTf9VImT2Wo9XJPg267BfX/R/c+E29rdtTp\ntv8fKojNQ9Vyznm5bxEqgAlhiNQoRDDaAFyPajqJRn2Zg/qrvxxVAwD1pZmrbUcDxahMp9/D+YVs\n/yIuQPVvjEMFiRTUl+0bwB+04we14xO1e2JQX8pf4MzpfwkwANV/4uvCTG2opqCZOJuebsZYh3us\n9vlA1Rhs2na1Vk671cBCnH8kXoK5VuATfiA1ChGMtqKaXXahFpjZDVRq5+4EXkR9+R1FpYwG1fm9\nCbVM6SZU4ADPq5eNQnX6tgJNwA+1n3NRncsRqL/or0KlaH5BK0+zVoYmD89u0/1s732LgN+hOpvL\nUAGnysPvQX//X1ALXc1HrSdvXxFxJ6q2sgPV4f4cqolrGyqQnUatmSyEEKYTpf2MRP0VPjaAZelq\n9s8WggqINwSwLEIIEbTeQHXI7se5EplZ/C/Oz/ZsgMsihBBCCCGEEEIIIYQQQgghhBBCCCGEEEII\nIYQQZvL/AXZDPGlMfW16AAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 5 }, { "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", "collapsed": false, "input": [ "post_prob = germany > argentina\n", "print('posterior prob Germany > Argentina', post_prob)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "posterior prob Germany > Argentina 0.704322830196\n" ] } ], "prompt_number": 6 }, { "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", "collapsed": false, "input": [ "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) " ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "posterior odds Germany > Argentina 2.38206700457\n", "Bayes factor 2.38206700457\n" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Bayes factor is 2.4, which is generally considered weak evidence.\n", "\n", "Now on to Step 4.\n", "\n", "### Step 4: The predictive distributions\n", "\n", "If we knew the actual goal scoring rate, $\\lambda$, we could predict how many goals each team work score in a rematch. The distribution of goals would be Poisson with parameter $\\lambda t$.\n", "\n", "We don't actually know $\\lambda$, but we can use the posterior distribution of $\\lambda$ to generate a predictive distribution for the number of additional goals." ] }, { "cell_type": "code", "collapsed": false, "input": [ "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, 15)\n", " metapmf[pred] = prob\n", "\n", " mix = thinkbayes2.MakeMixture(metapmf, label=label)\n", " return mix\n", "\n", "germany_pred = PredictiveDist(germany, label='germany')\n", "argentina_pred = PredictiveDist(argentina, label='argentina')\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "`PredictiveDist` takes the posterior distribution of $\\lambda$ and a duration (in units of games).\n", "\n", "It loops through the hypotheses in `suite`, computes the predictive distribution of goals for each hypothesis, and assembles a \"meta-Pmf\" which is a Pmf that maps from each predictive distribution to its probability.\n", "\n", "Finally, it uses `MakeMixture` to compute the mixture of the distributions. Here's what the predictive distributions look like." ] }, { "cell_type": "code", "collapsed": false, "input": [ "thinkplot.Hist(germany_pred, width=0.45, align='right')\n", "thinkplot.Hist(argentina_pred, width=0.45, align='left')\n", "thinkplot.Config(xlabel='predicted # goals', ylabel='probability', xlim=[-0.5, 7])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEPCAYAAACk43iMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHQNJREFUeJzt3Xt4VeWd6PFvCAiogMEgKJeiQovUEbSKF6SmtnbgCFqY\nKleB0SqnilqxPaDWMRzntE570LbDGaVUi+BMsVjHWqtiOzReOipQbtKCAgpy8YKCCNgLkZw/3pVk\nE0Kygll7Z2d9P8+znqy191rv/iXK+u33st4XJEmSJEmSJEmSJEmSJEmSpAYZDKwF1gFTD3FOCbAc\nWA2UZSUqSVLiCoH1QE+gFbACOKXGOccAfwS6RcfF2QpOkgQtEix7ACEJbAT2AfOBS2ucMwb4BbAl\nOn4vwXgkSTUkmQS6ApszjrdEr2XqDXQEfgcsBa5IMB5JUg0tEyy7IsY5rYAzgC8CRwIvAi8R+hAk\nSQlLMglsBbpnHHenutmn0mZCE9Cfo+05oB81kkC/fv0qVq5cmVykktQ8bQB65erDW0YB9ASOoPaO\n4T7AbwmdyEcCrwB9aymrIgl33HFHIuUmyZiTl2/xVlTkX8z5Fm9FRX7GTIwWmSRrAuXAZGAh4SZ/\nP7AGmBS9P4swfPRpYBWwH5gN/CnBmCRJGZJMAgBPRVumWTWO/2+0SZKyLMnRQU1eSUlJrkNoMGNO\nXr7FC/kXc77FC/kZcxwFuQ4gpqh5S5IUV0FBAdRzn0+6OUhSM9SxY0d27tyZ6zAUKSoqYseOHYd1\nrTUBSQ1WUFCA/yabjkP994hTE0h1n4AkpV2qm4PuXr2+UcubcmrOnsmQpMNiTUCSUswkIEkplurm\nIEmNZ9qi1xIt/64LP51o+WllTUCSalFeXp7rELLCJCCp2Vm2bBmnn3467du35/LLL2fkyJHcfvvt\nADzxxBP079+foqIiBg4cyCuvvFJ1Xc+ePfne977HaaedRrt27diwYQMtWrRgzpw59OjRg2OPPZb7\n7ruPJUuWcNppp1FUVMT1119fdf2GDRu48MILKS4uplOnTowbN45du3YdUP6MGTPo168fxxxzDKNG\njeKvf/0rAKeeeipPPPFE1bn79u2juLiYpGdQNglIalb+9re/MXz4cK688kp27tzJ6NGjeeyxxygo\nKGD58uVcddVVzJ49mx07djBp0iQuueQS9u3bV3X9/Pnzeeqpp/jggw8oLCwEYPHixaxfv5758+dz\n44038p3vfIdFixbxxz/+kZ///Oc899xzVdffdtttvPXWW6xZs4bNmzdTWlpa9V5BQQELFixg4cKF\nvPHGG6xatYo5c+YAMGHCBB566KGqc5988km6du1Kv379Ev17mQQkNSsvvfQSH3/8Mddffz2FhYUM\nHz6cAQMGUFFRwezZs5k0aRJnnXUWBQUFjB8/ntatW/PSSy8B4SZ9ww030LVrV1q3bl1V5u23384R\nRxzBRRddRLt27RgzZgzFxcWccMIJDBo0iOXLlwNw8skn88UvfpFWrVpRXFzMTTfdxLPPPntAfDfc\ncANdunShqKiIYcOGsWLFCgDGjh3Lr3/9a/bs2QPAvHnzuOKK5BdbNAlIala2bdtG164HrmTbvXtY\n32rTpk3MmDGDoqKiqm3Lli1s27btoHMzde7cuWq/bdu2Bx1X3rjfeecdRo0aRbdu3ejQoQNXXHEF\n77///gFldenSpdZrTzjhBAYOHMgjjzzCBx98wNNPP83YsWMP988Qm0lAUrNy/PHHs3Xr1gNee/PN\nN4Fwg7/tttvYuXNn1bZnzx5GjhxZdW401UKDVF5z6623UlhYyOrVq9m1axfz5s1j//799V5XqbJJ\naMGCBZx33nkcf/zxDY6loUwCkpqV8847j8LCQmbOnEl5eTm//OUvWbJkCQUFBVx99dXcd999LF68\nmIqKCvbu3XtAE8zhyJyzZ8+ePRx11FG0b9+erVu38v3vfz/2tQDDhw9n2bJl/OhHP2L8+PGHHVND\n+JyApEbRVMbxt2rVikcffZSvfe1r3HLLLQwZMoShQ4fSunVrPve5zzF79mwmT57MunXraNu2LYMG\nDapzrYD6agaZ799xxx2MHz+eDh060Lt3b8aNG8cPfvCDOq/NvL5NmzaMGDGChx9+mBEjRsT/pT+B\nVM8i6txB0uHJt1lEzz77bK699lomTJiQ61Dqdeedd7Ju3Trmzp0b+xpnEZWkDM899xxvv/025eXl\nPPjgg6xevZrBgwfnOqx67dixgwceeIBrrrkma59pEpDU7Lz66qtVD4Tdc889PPLIIweM6GmKZs+e\nTY8ePRgyZAjnn39+1j7X5qBGZHOQ0iLfmoOaO5uDJEmHxSQgSSlmEpCkFDMJSFKKmQQkKcVMApKU\nJc8//zx9+vTJdRgHcNoISY2isYdc15SPQ7BbtGjB+vXrOemkkwAYNGgQa9euzXFUB0q6JjAYWAus\nA6bW8n4JsAtYHm3fTjgeSary8ccfJ/4ZTf15iiSTQCEwk5AI+gKjgVNqOe9Z4PRo++cE45GUEnfd\ndRe9evWiffv2fPazn+Wxxx4DYM6cOQwcOJApU6ZQXFzM9OnT2bFjB8OGDaNDhw4MGDCAb3/72wwa\nNKiqrLVr13LRRRdx7LHH0qdPHxYsWFD13sSJE7nuuusYOnQo7du355xzzuH1118H4POf/zwA/fr1\no127dixYsICysrID1iuoa7nJDz74gKFDh3LcccfRsWNHhg0bdtAU2Y0hySQwAFgPbAT2AfOBS2s5\nL1+eWpaUJ3r16sULL7zAhx9+yB133MG4ceN4++23gbBU5Mknn8y7777LrbfeyrXXXku7du145513\nePDBB5k7d27VzJ579+7loosuYty4cWzfvp358+dz7bXXsmbNmqrPevjhhyktLWXnzp306tWL2267\nDaBqyclVq1axe/duLrvssoPirGu5yf3793PVVVfx5ptv8uabb9K2bVsmT57c6H+rJJNAV2BzxvGW\n6LVMFcB5wErgSUKNQZI+ka9+9atVK3hdfvnl9O7dm8WLFwNhBa/rrruOFi1aVE07PX36dNq0acMp\np5zChAkTqppwnnjiCU488UQmTJhAixYt6N+/PyNGjDigNjBixAjOPPNMCgsLGTt2bNVykXEdarnJ\njh07Mnz4cNq0acPRRx/NrbfeetBSlY0hyY7hOA1hy4DuwEfAEOAxoGlMSi4pb82dO5d77rmHjRs3\nAmGxl/fee4/CwsIDmmO2b99OeXn5Aa9169atan/Tpk28/PLLFBUVVb1WXl5eteBLQUHBIZeajKvm\ncpOVS11+9NFH3HTTTSxcuJCdO3dW/R4VFRWHtfrZoSSZBLYSbvCVuhNqA5l2Z+w/Bfwb0BHYUbOw\n0tLSqv2SkpI6F4GQlF6bNm3immuuYdGiRZx77rkUFBRw+umnV327z7yBdurUiZYtW7J582Z69+4N\nwObN1Q0YPXr04IILLuCZZ57J7i8BzJgxg9dee43Fixdz3HHHsWLFCs4444w6k0BZWRllZWUN+pwk\nk8BSoDfQE9gGjCR0DmfqDLxLqDUMIPQPHJQA4MAkIEmHsnfvXgoKCiguLmb//v3MnTuX1atXAweP\n1CksLGTEiBGUlpbyk5/8hE2bNjFv3jw+9alPAXDxxRczbdo0Hnrooap1iFesWEG7du3o06dPvSN/\nOnfuzIYNG6qGiDbEnj17aNu2LR06dGDHjh1Mnz693mtqfkGOc02SSaAcmAwsJIwUuh9YA0yK3p8F\nfBX4enTuR8CoBOORlKCmMo6/b9++3HzzzZx77rm0aNGC8ePHc/7551ct5VjzW/TMmTOZOHEiXbp0\noU+fPowePZqlS5cC0K5dO5555hmmTJnClClT2L9/P/379+fuu+8GDl4esvK1SqWlpUyYMIE///nP\nzJ49m06dOtXZlJNZ3je+8Q3GjBlDcXExXbt2ZcqUKTz++OON8jc64DMbvcRkuJ6A1IQ05/UEpk6d\nyrvvvstPf/rTXIcSm+sJSNJhevXVV1m1ahUVFRUsXryYBx54gOHDh+c6rKxx2ghJqbZ7925Gjx7N\ntm3b6Ny5M9/85je55JJLch1W1pgEJKXamWeeybp163IdRs7YHCRJKWYSkKQUMwlIUorZJyCpwYqK\nihp16gJ9MpnTWjSUSUBSg+3YUeuD/cpDNgdJUoqZBCQpxUwCkpRiJgFJSjGTgCSlmElAklLMJCBJ\nKWYSkKQUMwlIUoqZBCQpxUwCkpRiJgFJSjGTgCSlmElAklLMJCBJKWYSkKQUMwlIUoqZBCQpxUwC\nkpRiJgFJSjGTgCSlWNJJYDCwFlgHTK3jvLOAcmBEwvFIkjIkmQQKgZmERNAXGA2ccojz/gV4GihI\nMB5JUg1JJoEBwHpgI7APmA9cWst51wOPANsTjEWSVIskk0BXYHPG8ZbotZrnXArcGx1XJBiPJKmG\nlgmWHeeG/gNgWnRuAXU0B5WWllbtl5SUUFJS8smik6RmpqysjLKysgZdk2Qb/DlAKaFPAOAWYD+h\n/b/S6xkxFAMfAVcDj9coq6KiovErCXevXt+o5U05tVejlidJn0RBQQHUc59PsiawFOgN9AS2ASMJ\nncOZTsrY/ynwKw5OAJKkhCSZBMqBycBCwgig+4E1wKTo/VkJfrYkKYZ8GZJZUVFRwbRFrzVqoccd\n17j94jYHSWpK4jQH+cSwJKWYSUCSUswkIEkpZhKQpBQzCUhSipkEJCnFTAKSlGImAUlKsThJ4FHg\n4pjnSpLySJwb+73AWMLaAHcBn0k0IklS1sRJAr8BxgBnEBaI+S/gv4F/BFolFpkkKXFxm3iOBSYC\nXwOWAT8CPkdIEJKkPBVnFtH/BPoA84BhwFvR6/OBPyQUlyQpC+IkgdnAkzVeaw38lVAbkCTlqTjN\nQf+nltdebOxAJEnZV1dN4HjgBKAtoVO4gLAWcHvgyORDkyQlra4k8PfABKArMCPj9d3ArUkGJUnK\njrqSwJxo+wfgF9kIRpKUXXUlgSsII4J6AlMyXq9sFro7ubAkSdlQVxKobPdvR7jpVyqocSxJylN1\nJYFZ0c/SLMQhScqBupLAv9bxXgVwQyPHIknKsrqSwB8IN/uCWt6zOUiSmoH6RgdJkpqxupLAD4Eb\ngV/V8l4FcEkiEUmSsqauJDA3+jmjlvdsDpKkZqC+PgGAMsKEcX2A/cCrwN+SDUuSlA1xZhG9GLgP\neD06PgmYxMEzi0qS8kycWUTvBr4AXBBtJcA9McsfDKwF1gFTa3n/UmAlsJxQ87gwZrmSpEYQpybw\nIWF94UqvR6/VpxCYCXwJ2AosAR4H1mSc81vgl9H+3xEWsOkVo2xJUiOoKwn8Q/RzKaHp5+fR8WXR\na/UZQEgeG6Pj+YRv/plJYG/G/tHAezHKlSQ1krqSwDCqRwG9S2gKAtgOtIlRdldgc8bxFuDsWs77\nCvBdwvoFX45Rbl6Ztui1Ri3vrgs/3ajlSUq3upLAxE9YdtxhpI9F2yDCrKWfqe2k0tJSXnjjfQB6\n9BtAj/615RNJSq+ysjLKysoadE2cPoG2wFVA32i/8uZ+ZT3XbQW6Zxx3J9QGDuX5KJ5jgfdrvlla\nWspfGvlbtSQ1JyUlJZSUlFQdT58+vd5r4owOmgd0Joz0KSPczPfEuG4p0JuwHsERwEhCx3Cmk6me\nm+iM6OdBCUCSlIw4NYFewFcJnboPAv8BvBDjunJgMrCQMFLofkKn8KTo/VmEzufxwD5CYhnVgNgl\nSZ9QnCRQ+XTwLsIwzreBTjHLfyraMs3K2P9etEmSciBOEpgNdAS+TWjOORq4PcmgJEnZETcJADwL\nnJhgLJKkLIvTMVxMWGVsObCMMMX0sUkGJUnKjjhJYD7hYbERhA7i7cDDSQYlScqOOM1BXYA7M47/\nmTDcU5KU5+LUBJ4BRkfntiAkgGeSDEqSlB111QT2UP108DcID41BSAR7gZsTjEuSlAV1JYGjsxaF\nJCkn4vQJQHha+POEmsGz1L74vCQpz8RJAncBZwH/Tpjn5wbgPOCWBOPSIdy9en39JzXAlFNdw0dK\ns7hrDPcHPo6O5wArMAlIUt6LMzqoAjgm4/gY4q8VIElqwuLUBL5LeFL4d4TmoAuAaUkGJUnKjvqS\nQAtgP3AuoV+ggpAA3ko4LklSFtSXBPYD/4swTcQvkw9HkpRNcfoEfgN8k7CiWMeMTZKU5+L0CYwi\nNANdV+N1p5WWpDwXJwmcQkgA5xOah14A7k0yKElSdsRJAnOBDwnrCBQAY6LXLkswLklSFsRJAp8F\n+mYcLwL+lEw4kqRsitMxvIwwRLTSOcAfkglHkpRNcWoCZwK/BzYTOoh7AK8Cr0THpyUWnSQpUXGS\nwODEo5Ak5UScJLAx6SAkSbkRp09AktRMmQQkKcVMApKUYiYBSUqxbCSBwcBaYB0wtZb3xwIrgVWE\noagOOZWkLIm70PzhKgRmAl8CtgJLgMeBNRnnvE5YxH4XIWH8mPBAmiQpYUnXBAYA6wnDTPcB84FL\na5zzIiEBALwMdEs4JklSJOkk0JXwpHGlLdFrh3IV8GSiEUmSqiTdHNSQBem/AFwJDKztzdLSUl54\n430AevQbQI/+Z3/y6CSpGSkrK6OsrKxB1ySdBLYSViSr1J1QG6jpNGA2oU9gZ20FlZaW8pdFrzV6\ngJLUXJSUlFBSUlJ1PH369HqvSbo5aCnQG+gJHAGMJHQMZ+oBPAqMI/QfSJKyJOmaQDkwGVhIGCl0\nP2Fk0KTo/VnAPwFFVK9Wto/QoSxJSljSSQDgqWjLNCtj/2vRJknKMp8YlqQUMwlIUoqZBCQpxUwC\nkpRiJgFJSjGTgCSlmElAklIsG88JKI9MS2Bqjrsu/HSjlympcVgTkKQUMwlIUoqZBCQpxUwCkpRi\nJgFJSjGTgCSlmElAklLMJCBJKWYSkKQUMwlIUoqZBCQpxUwCkpRiJgFJSjGTgCSlmElAklLMJCBJ\nKWYSkKQUMwlIUoqZBCQpxUwCkpRi2UgCg4G1wDpgai3v9wFeBP4C3JyFeCRJkZYJl18IzAS+BGwF\nlgCPA2syznkfuB74SsKxSJJqSLomMABYD2wE9gHzgUtrnLMdWBq9L0nKoqSTQFdgc8bxlug1SVIT\nkHRzUEXC5SsP3L16faOWN+XUXo1anpRmSSeBrUD3jOPuhNpAg5WWlvLCG+8D0KPfAHr0P/uTRydJ\nzUhZWRllZWUNuibpJLAU6A30BLYBI4HRhzi3oK6CSktL+cui1xo1OElqTkpKSigpKak6nj59er3X\nJJ0EyoHJwELCSKH7CSODJkXvzwK6EEYNtQf2AzcCfYE9CccmSamXdBIAeCraMs3K2H+bA5uMJElZ\n4hPDkpRiJgFJSjGTgCSlWDb6BKRETWvkUWN3XfjpRi1PasqsCUhSipkEJCnFTAKSlGImAUlKMZOA\nJKWYSUCSUswkIEkpZhKQpBQzCUhSipkEJCnFTAKSlGImAUlKMZOAJKWYSUCSUswkIEkp5noCUg13\nr17fqOVNObXXQa+5BoKaCmsCkpRiJgFJSjGTgCSlmElAklLMJCBJKWYSkKQUMwlIUor5nICkevlc\nQ/OVdBIYDPwAKAR+AvxLLef8CBgCfARMBJYnHJPU7GTjATc1T0k2BxUCMwmJoC8wGjilxjn/A+gF\n9AauAe5NMJ6DrF/ycjY/rlEYc/LyLV7Iv5jLyspyHUKD5WPMcSSZBAYA64GNwD5gPnBpjXMuAR6M\n9l8GjgE6JxjTATbk2T8cMOZsyLd4If9izscbaj7GHEeSzUFdgc0Zx1uAs2Oc0w14J8G4JOXYi+/u\nSLwJq7H7Mdo0amlNR5JJoCLmeQWHeZ0kZU1jJ65sTCyYa+cAT2cc3wJMrXHOfcCojOO11N4ctIKQ\nHNzc3Nzc4m+NW91qoJbABqAncAThRl5bx/CT0f45wEvZCk6SlLwhwKuEbHRL9NqkaKs0M3p/JXBG\nVqOTJEmSlD6DCX0Q6zi4r6IpeoAwauqVXAcSU3fgd8AfgdXADbkNJ5Y2hKHKK4A/Ad/NbTixFRIe\nsvxVrgOJaSOwihDz4tyGEtsxwCPAGsL/G+fkNpw6fYbwt63cdpEf//6yqpDQBNUTaEXt/RVNzSDg\ndPInCXQB+kf7RxOaBZv63xjgyOhnS0If1fk5jCWuKcC/A4/nOpCY3gA65jqIBnoQuDLabwl0yGEs\nDdECeIvwpeyQJ6RRnAfZmprngZ25DqIB3iYkV4A9hG9QJ+QunNg+in4eQfiysCOHscTRjTDA4icc\nPNy6KcunWDsQvoQ9EB2XE75d54MvEQbobD7UCWlNArU9pNY1R7GkQU9CLSYfHmttQUhe7xCas/6U\n23DqdQ/wLWB/rgNpgArgt8BS4OocxxLHicB24KfAMmA21TXGpm4U8B91nZDWJFCR6wBS5GhCW+qN\nhBpBU7ef0IzVDfg8UJLTaOo2FHiX0O6bT9+sBxK+FAwBriN8y27KWhJGLv5b9HMvMC2nEcVzBDAM\nWFDXSWlNAls5sI2sO6E2oMbVCvgF8BDwWI5jaahdwK+BM3MdSB3OI8y/9QbwM+BCYG5OI4rnrejn\nduA/Cc2zTdmWaFsSHT9CfgxnHwL8gfB3Vg1xHmRrinqSPx3DBYQb0j25DqQBigmjQADaAs8BX8xd\nOA1yAfkxOuhIoF20fxTwe+DLuQsntueAykUQSql9WvymZj4wIddBNGW1PcjWlP0M2Ab8ldCf8Y+5\nDade5xOaVlZQPVRtcE4jqt/fEdp8VxCGMH4rt+E0yAXkx+igEwl/3xWEocP58G8PoB+hJrASeJSm\nPzroKOA9qhOuJEmSJEmSJEmSJEmSJEmSJDUXJVQ/RDWMuqcN7wB8/TA+oxS4uY73jwJ+E+0/T/JP\n6tcXj1QlrdNGKP8dzv+7v6LuJz2LgGsPo9z65qI6F3gxKn8vyU/25txYis0koKamJ2Gxn4cIM3gu\nIEzhAGHq77sI86FcRphu4L+j458TvnFDeDJ5TfT68IyyJwL/Gu13JsxbU/n06rlR2ScTnm6uTBbf\nIix8spLwDbvSbYQnzp8nLOJRm8qy5gFjCLNm9iM8ldyplvNvj3735wkzP1Z+m+9PWNug8mnVyqkt\nro5iW0GYz6YtB7uBsLDPSsJT55LUpPUkfFM+Nzq+n+qb4RvAN6P9YuBZqm98Uwk30TbAm4QbMMDD\nVE+nMJHqJPAw1asttQDaA5/iwLmZvgzMyjjnV4QZLz9HmFaiDeGx/HWEhV0O5QlCLeCfCNOV1OYs\nQsI4gjDz6msZZa6ieqbN6VTPx5S5MMudwORo/46Ma7cSJvIj+h2lA1gTUFO0mdB8AqFGkLm618PR\nz3OAvoSawHJgPNCD8K38DcIEgZXX1zbN8heAe6P9/cCHtZz35WhbTqhVfAboHcXzKPAXYDchydQ1\nlfNxhAWB+hFu6LUZSJhp9W+EKbcr+zHaE/oqno+OHyRMcQ1hrqPnozLHEv4eNa0i1CrGAh/XEaNS\nyiSgpiizTbugxvHejP3fEOalPx34LLUvUFLXzTnOHPzfzfiMT1O9ulTmtYcq515CzaI3IZH8PaFW\ncGMt51bELDPz9TmEPozTCDWEtrWcdzHw/whTHy8hrJYmVTEJqCnqQfVC3mOo/hac6WXCt+fKZp+j\nCDfbtYQmpZOi10cf4jP+i+qRQIWEb9y7OXDWxYWEdWUr+xq6EtrynwO+QnVz0FBq74z9OuHm/L+j\n839NSCY/rOXc3xNGL7UmNAddHL3+IaEWUVkbugIoi/aPJizj2QoYlxFDQcbPHtH50wg1isrfRQLC\nvPpSU/MqYcWpBwidmpXNNpk32u2ENv6fEW6cEDpr1wHXEG64HxESyFEZ11eWcSPwY+AqQjPJ/yQk\nlt8Tvr0/SehnOIXqpqndhJvtckKz1ErCyl6L6/hdLiCsqzCI6pt3bZYSmpVWEZa2fIXqdWwnAPcR\n5uLfQPU04rdHMW+Pfh5d4/csJHRKdyAkhB8SkookNVk9yZ+FcxpbZbI6ktB00z+HsSglrAmoKUrr\nOPcfEzp32xDa+1fkNBpJkiRJkiRJkiRJkiRJkiRJ+ej/A/OIHUm+x5IfAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the predictive distributions, we can compute probabilities for the outcomes of a rematch." ] }, { "cell_type": "code", "collapsed": false, "input": [ "win = germany_pred > argentina_pred\n", "lose = germany_pred < argentina_pred\n", "tie = 1 - (win + lose)\n", "\n", "print('posterior prob Germany wins rematch', win)\n", "print('posterior prob Argentina wins rematch', lose)\n", "print('posterior prob tie', tie)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "posterior prob Germany wins rematch 0.450511695667\n", "posterior prob Argentina wins rematch 0.225540779355\n", "posterior prob tie 0.323947524978\n" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the score is tied after 90 minutes the teams play an additional 30 minute period. We can use `PredictiveDist` again to compute the distribution of scores after 1/3 of a game:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "germany_pred_overtime = PredictiveDist(germany, 1/3, label='germany')\n", "argentina_pred_overtime = PredictiveDist(argentina, 1/3, label='argentina')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "win = germany_pred_overtime > argentina_pred_overtime\n", "lose = germany_pred_overtime < argentina_pred_overtime\n", "tie = 1 - (win + lose)\n", "\n", "print('posterior prob Germany wins rematch', win)\n", "print('posterior prob Argentina wins rematch', lose)\n", "print('posterior prob tie', tie)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "posterior prob Germany wins rematch 0.256156067229\n", "posterior prob Argentina wins rematch 0.13931727161\n", "posterior prob tie 0.604526661162\n" ] } ], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "print('Total prob Germany wins a rematch', 0.45 + 0.32 * 0.26)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Total prob Germany wins a rematch 0.5332\n" ] } ], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "print('Total prob Argentina wins', 0.23 + 0.32 * 0.14)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Total prob Argentina wins 0.2748\n" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "print('Prob of draw after overtime', 0.32 * 0.60)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Prob of draw after overtime 0.192\n" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 } ], "metadata": {} } ] }