{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bayesian Zig Zag\n", "\n", "Developing probabilistic models using grid methods and MCMC.\n", "\n", "Thanks to Chris Fonnesback for his help with this notebook, and to Colin Carroll, who added features to pymc3 to support some of these examples.\n", "\n", "Copyright 2018 Allen Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n", "\n", "import numpy as np\n", "import pymc3 as pm\n", "\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulating hockey\n", "\n", "I'll model hockey as a Poisson process, where each team has some long-term average scoring rate, lambda, in goals per game.\n", "\n", "For the first example, I'll assume that lambda is somehow known to be 2.7. Since regulation play (as opposed to overtime) is 60 minutes, we can compute the goal scoring rate per minute." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.045000000000000005, 0.0020250000000000003)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lam_per_game = 2.7\n", "min_per_game = 60\n", "lam_per_min = lam_per_game / min_per_game\n", "lam_per_min, lam_per_min**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we assume that a goal is equally likely during any minute of the game, and we ignore the possibility of scoring more than one goal in the same minute, we can simulate a game by generating one random value each minute." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def simulate_game(p, n=60):\n", " goals = np.random.choice([0, 1], n, p=[1-p, p])\n", " return np.sum(goals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And simulate 10 games." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 6, 2, 3, 0, 1, 3, 2, 2, 1]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "size = 10\n", "sample = [simulate_game(lam_per_min) for i in range(size)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we simulate 1000 games, we can see what the distribution looks like. The average of this sample should be close to lam_per_game." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2.706, 2.7)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "size = 1000\n", "sample_sim = [simulate_game(lam_per_min) for i in range(size)]\n", "np.mean(sample_sim), lam_per_game" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PMFs\n", "\n", "To visualize distributions, I'll start with a probability mass function (PMF), which I'll implement using a Counter.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from collections import Counter\n", "\n", "class Pmf(Counter):\n", " \n", " def normalize(self):\n", " \"\"\"Normalizes the PMF so the probabilities add to 1.\"\"\"\n", " total = sum(self.values())\n", " for key in self:\n", " self[key] /= total\n", " \n", " def sorted_items(self):\n", " \"\"\"Returns the outcomes and their probabilities.\"\"\"\n", " return zip(*sorted(self.items()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are some functions for plotting PMFs." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "plot_options = dict(linewidth=3, alpha=0.6)\n", "\n", "def underride(options):\n", " \"\"\"Add key-value pairs to d only if key is not in d.\n", "\n", " options: dictionary\n", " \"\"\"\n", "\n", " for key, val in plot_options.items():\n", " options.setdefault(key, val)\n", " return options\n", "\n", "def plot(xs, ys, **options):\n", " \"\"\"Line plot with plot_options.\"\"\"\n", " plt.plot(xs, ys, **underride(options))\n", "\n", "def bar(xs, ys, **options):\n", " \"\"\"Bar plot with plot_options.\"\"\"\n", " plt.bar(xs, ys, **underride(options))\n", "\n", "def plot_pmf(sample, **options):\n", " \"\"\"Compute and plot a PMF.\"\"\"\n", " pmf = Pmf(sample)\n", " pmf.normalize()\n", " xs, ps = pmf.sorted_items()\n", " bar(xs, ps, **options)\n", " \n", "def decorate_pmf_goals():\n", " \"\"\"Decorate the axes.\"\"\"\n", " plt.xlabel('Number of goals')\n", " plt.ylabel('PMF')\n", " plt.title('Distribution of goals scored')\n", " legend()\n", " \n", "def legend(**options):\n", " \"\"\"Draw a legend only if there are labeled items.\n", " \"\"\"\n", " ax = plt.gca()\n", " handles, labels = ax.get_legend_handles_labels()\n", " if len(labels):\n", " plt.legend(**options)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the results from the simulation look like." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHaBJREFUeJzt3Xu8lWWd9/HPV46ShwgYGzkI5iGRLahIlqY0MkajqfnoCKMONfowNYM1GaSYmVlPUk096ZNlpIyVlprhK2pQB0epzANsFAI8JBLBDiuE1CRPwO/54742s9isva+9Yd97bff+vl8vXqx1H677t9aC9V336boUEZiZmbVkj1oXYGZmnZ/DwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYbtN0vWSPt1ObQ2T9JKkHun5QkkXtkfbqb27JE1pr/basN3PS3pO0u87eLtrJE3oyG2WQdKVkm6udR3dWc9aF2Cdm6Q1wH7AFmAr8DjwXWB2RGwDiIgPt6GtCyPi3uaWiYi1wF67V/X27V0JHBQR51W0/772aLuNdQwFPgEcEBF/7Ojtm7UH71lYa7w/IvYGDgBmAZcAN7b3RiR11R8vBwAbHRSFLvw5d2kOC2u1iHghIuYB5wBTJI0CkHSTpM+nxwMl/VTS85I2SfqFpD0kfQ8YBvwkHWb6pKThkkLSBZLWAvdVTKv8QnmbpEWSXpD0Y0lvSdsaL6mhssbGwy6SJgKXAeek7S1L87cf1kp1XS7pt5L+KOm7kvZN8xrrmCJpbTqE9Knm3htJ+6b1N6T2Lk/tTwAWAPunOm5qZv1PSnpW0npJF6ZtH9RS22ne2yTdJ2ljqvEWSW9uZhvjJNVLelHSHyR9tZnlqn6Gad5QSXNTLRslfb0N7+X2zzlNP1bSg2k7yySNr6hhhKSfSfqzpAXAwObee+sYDgtrs4hYBDQA764y+xNp3iCKw1eXFavE+cBair2UvSLiSxXrnAgcBry3mU3+I/BPwP4Uh8OubUWNdwNfAG5L2xtdZbEPpj/vAQ6kOPz19SbLHA8cCpwEXCHpsGY2+f+AfVM7J6aaP5QOub0PWJ/q+GDTFVOwXQxMAA5K62fbblwduJrivTkMGApc2UyN1wDXRMQ+wNuA25tZrupnqOI80k+B3wLDgcHArWmdD5J/L7d/zpIGA/8JfB54CzAd+JGkQWnZ7wNLKELic0CHn2eyHTksbFetp/hP3tTrwF9THJ9/PSJ+EfkOyK6MiM0R8XIz878XESsiYjPwaeDv0xfX7joX+GpErI6Il4CZwKQmezWfjYiXI2IZsAzYKXRSLecAMyPizxGxBvgKcH4r6/h74D8iYmVE/AX4bGvbjohVEbEgIl6NiA3AV9k5bBq9DhwkaWBEvBQRD7ewXLXPcBxFKM1In9crEfFAWqc172Xl53weMD8i5kfEtohYANQDfydpGHAM8On0un4O/KSV76WVxGFhu2owsKnK9C8Dq4D/krRa0qWtaGtdG+b/FuhF+xyW2D+1V9l2T4pf040qr176C9VPvg8Eeldpa3Ab6qh8jZWPW2xb0l9JulXS7yS9CNxM8+/NBcAhwJOSFks6tZnlmvsMhwK/jYgtzbyG3HtZ+boOAM5Oh6Cel/Q8xV7cX6e2/pR+HFS2ZzXksLA2k3QMxZfVA03npV+/n4iIA4H3AxdLOqlxdjNN5vY8hlY8Hkbxy/c5YDPQr6KuHhSHTlrb7nqKL63KtrcAf8is19Rzqaambf2ules/CwypeF75enNtX03xOo9Ih5fOozg0tZOIeDoiJgN/BXwRuEPSm6os19xnuA4YpuonqFvzXlZ+Huso9hjfXPHnTRExK70f/ZvUNqzaa7KO47CwVpO0T/o1eitwc0Qsr7LMqZIOkiTgRYrLbbem2X+gOJ7dVudJGimpH3AVcEdEbAV+DfSVdIqkXsDlQJ+K9f4ADG88OVvFD4CPp5Ope/E/5ziq/XJuVqrlduD/SNpb0gEU5yBae1/A7cCHJB2WXuMVbWh7b+Al4Pl0HmBGcxuRdJ6kQemS5+fT5K1VlmvuM1xE8UU+S9KbJPWVdFxara3v5c3A+yW9V1KP1NZ4SUMi4rcUh6Q+K6m3pOMpQstqyGFhrfETSX+m+DX4KYrj4h9qZtmDgXspvsAeAr4REQvTvKuBy9Nhh+lt2P73gJsoDgn1BT4KxdVZwL8AN1D80t5McWK20Q/T3xslPVql3Tmp7Z8DvwFeAS5qQ12VLkrbX02xx/X91H5WRNxFcdL+forDPw+lWa+2ou3PAkcBL1CcMJ7bwqYmAislvURxsntSRLxSZbmqn2EKrvdTnIRfS/Fen5PWadN7GRHrgNMpTp5voPi3NYP/+U76B+AdFIc6P0Nxb4/VkDz4kVnnkq64WgH0aetejllZvGdh1glI+kA65NKf4nzCTxwU1pk4LMw6h3+mOBzzDMX5gY/UthyzHfkwlJmZZXnPwszMsrpMh14DBw6M4cOH17oMM7M3lCVLljwXEYNyy3WZsBg+fDj19fW1LsPM7A1FUqvujvdhKDMzy3JYmJlZlsPCzMyyusw5CzN743v99ddpaGjglVeq9UJiu6Nv374MGTKEXr167dL6pYZFGtTlGqAHcEPqUbJy/sXAhRS9U24A/il1IoakrUBjR3VrI+K0Mms1s9praGhg7733Zvjw4RT9GFp7iAg2btxIQ0MDI0aM2KU2SjsMlbqLvo5ilLCRwGRJI5ss9hgwNiKOAO4AKkdPezkixqQ/DgqzbuCVV15hwIABDop2JokBAwbs1h5bmecsxgGr0shZr1F0a3165QIRcX8aGQzgYXbs09/MuiEHRTl2930tMywGs+PIWA20PHLYBcBdFc/7qhhc/mFJZ1RbQdLUtEz9hg0bdr9iMzOrqsxzFtVirGpHVJLOA8ay49jBwyJivaQDgfskLY+IZ3ZoLGI2MBtg7Nix7uTKrIuZOXen8bV2y9Vn1u3SehdeeCEXX3wxI0c2PZLedo03EA8c2PzIwF/4whe47LLLtj9/17vexYMPPrjb294dZYZFAzsODzmEYujFHUiaQDGgzokR0TjYCxGxPv29WtJC4EiKHjmtHbX3f8ZqdvU/qFlnccMNN3To9pqGRa2DAso9DLUYODgNs9gbmATMq1xA0pHAt4DTIuKPFdP7S+qTHg8EjgMeL7FWMzMANm/ezCmnnMLo0aMZNWoUt912G+PHj9/endBee+3FJZdcwtFHH82ECRNYtGgR48eP58ADD2TevOIr7qabbmLatGnb2zz11FNZuHDhTts644wzOProozn88MOZPXs2AJdeeikvv/wyY8aM4dxzz92+TSiuapoxYwajRo2irq6O2267DYCFCxcyfvx4zjrrLN7+9rdz7rnn0t49ipcWFmnglmnAPcATwO0RsVLSVZIar276MrAX8ENJSyU1hslhQL2kZRRDTc6KCIeFmZXu7rvvZv/992fZsmWsWLGCiRMn7jB/8+bNjB8/niVLlrD33ntz+eWXs2DBAu68806uuOKKZlqtbs6cOSxZsoT6+nquvfZaNm7cyKxZs9hzzz1ZunQpt9xyyw7Lz507l6VLl7Js2TLuvfdeZsyYwbPPPgvAY489xte+9jUef/xxVq9ezS9/+cvdeyOaKPU+i4iYD8xvMq1yMPoJzaz3IOBjF2bW4erq6pg+fTqXXHIJp556Ku9+97t3mN+7d+/tAVJXV0efPn3o1asXdXV1rFmzpk3buvbaa7nzzjsBWLduHU8//TQDBgxodvkHHniAyZMn06NHD/bbbz9OPPFEFi9ezD777MO4ceMYMqS4oHTMmDGsWbOG448/vk31tMR3cJuZVTjkkENYsmQJ8+fPZ+bMmZx88sk7zO/Vq9f2y1D32GMP+vTps/3xli3FSLg9e/Zk27Zt29epdn/DwoULuffee3nooYfo168f48ePz94H0dKhpcY6AHr06LG9lvbivqHMzCqsX7+efv36cd555zF9+nQeffTRNrcxfPhwli5dyrZt21i3bh2LFi3aaZkXXniB/v37069fP5588kkefvjh7fN69erF66+/vtM6J5xwArfddhtbt25lw4YN/PznP2fcuHFtrm9XeM/CzDqtWlxJt3z5cmbMmMEee+xBr169+OY3v8n06dPb1MZxxx3HiBEjqKurY9SoURx11FE7LTNx4kSuv/56jjjiCA499FCOPfbY7fOmTp3KEUccwVFHHbXDeYsPfOADPPTQQ4wePRpJfOlLX+Ktb30rTz755K6/4FbqMmNwjx07Njz4Udv50lnrTJ544gkOO+ywWpfRZVV7fyUtiYixuXV9GMrMzLIcFmZmluWwMLNOpascGu9sdvd9dViYWafRt29fNm7c6MBoZ43jWfTt23eX2/DVUGbWaQwZMoSGhgbci3T7axwpb1c5LMys0+jVq9cuj+Rm5fJhKDMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVmWw8LMzLIcFmZmluWwMDOzrJ61LsC6r5lzl5e+javPrCt9G2bdgfcszMwsq9SwkDRR0lOSVkm6tMr8iyU9LulXkv5b0gEV86ZIejr9mVJmnWZm1rLSwkJSD+A64H3ASGCypJFNFnsMGBsRRwB3AF9K674F+AzwDmAc8BlJ/cuq1czMWlbmnsU4YFVErI6I14BbgdMrF4iI+yPiL+npw8CQ9Pi9wIKI2BQRfwIWABNLrNXMzFpQZlgMBtZVPG9I05pzAXBXW9aVNFVSvaT6DRs27Ga5ZmbWnDLDQlWmRdUFpfOAscCX27JuRMyOiLERMXbQoEG7XKiZmbWszLBoAIZWPB8CrG+6kKQJwKeA0yLi1basa2ZmHaPMsFgMHCxphKTewCRgXuUCko4EvkURFH+smHUPcLKk/unE9slpmpmZ1UBpN+VFxBZJ0yi+5HsAcyJipaSrgPqImEdx2Gkv4IeSANZGxGkRsUnS5ygCB+CqiNhUVq1mZtayUu/gjoj5wPwm066oeDyhhXXnAHPKq87MzFrLd3CbmVmWw8LMzLLckWAn4A71zKyz856FmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZZ7nbVuyT39mrWN9yzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7OsUsNC0kRJT0laJenSKvNPkPSopC2Szmoyb6ukpenPvDLrNDOzlpU2+JGkHsB1wN8CDcBiSfMi4vGKxdYCHwSmV2ni5YgYU1Z9ZmbWemWOlDcOWBURqwEk3QqcDmwPi4hYk+ZtK7EOMzPbTWUehhoMrKt43pCmtVZfSfWSHpZ0RvuWZmZmbVHmnoWqTIs2rD8sItZLOhC4T9LyiHhmhw1IU4GpAMOGDdv1Ss3MrEVl7lk0AEMrng8B1rd25YhYn/5eDSwEjqyyzOyIGBsRYwcNGrR71ZqZWbPKDIvFwMGSRkjqDUwCWnVVk6T+kvqkxwOB46g412FmZh2rtLCIiC3ANOAe4Ang9ohYKekqSacBSDpGUgNwNvAtSSvT6ocB9ZKWAfcDs5pcRWVmZh2ozHMWRMR8YH6TaVdUPF5McXiq6XoPAnVl1mZmZq3nO7jNzCzLYWFmZlkOCzMzy3JYmJlZVothIemmisdTSq/GzMw6pdyexeiKxx8rsxAzM+u8cmHRlu45zMysi8rdZzFE0rUU/Tw1Pt4uIj5aWmVmZtZp5MJiRsXj+jILMTOzzqvFsIiI73RUIWZm1nm1GBa54Uwj4rT2LcfMzDqj3GGod1IMYPQD4BGqj1FhZmZdXC4s3koxhvZk4B+A/wR+EBErW1zLzMy6lBYvnY2IrRFxd0RMAY4FVgELJV3UIdWZmVmnkO2iPA1CdArF3sVw4FpgbrllmZlZZ5I7wf0dYBRwF/DZiFjRIVWZmVmnktuzOB/YDBwCfExS4x3dAiIi9imzODMz6xxy91m4V1ozM8sehuoLfBg4CPgVMCeNrW1mZt1Ibs/hO8BYYDnwd8BXSq/IzMw6ndw5i5ERUQcg6UZgUfklmZlZZ5Pbs3i98YEPP5mZdV+5PYvRkl5MjwXsmZ77aigzs24kdzVUj44qxMzMOi9fGmtmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZpYaFpImSnpK0StKlVeafIOlRSVskndVk3hRJT6c/U8qs08zMWlZaWEjqAVwHvA8YCUyWNLLJYmuBDwLfb7LuW4DPAO8AxgGfkdS/rFrNzKxlZe5ZjANWRcTqiHgNuBU4vXKBiFgTEb8CtjVZ973AgojYFBF/AhYAE0us1czMWlBmWAwG1lU8b0jT2m1dSVMl1Uuq37Bhwy4XamZmLSszLFRlWrTnuhExOyLGRsTYQYMGtak4MzNrvTLDogEYWvF8CLC+A9Y1M7N2VmZYLAYOljRCUm9gEjCvleveA5wsqX86sX1ymmZmZjVQWlikYVinUXzJPwHcHhErJV0l6TQAScdIagDOBr4laWVadxPwOYrAWQxclaaZmVkN5IZV3S0RMR+Y32TaFRWPF1McYqq27hxgTpn1mZlZ6/gObjMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8sq9Q5uM9vZzLnLS9/G1WfWlb4N6168Z2FmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVmWw8LMzLIcFmZmluWwMDOzLIeFmZll9SyzcUkTgWuAHsANETGryfw+wHeBo4GNwDkRsUbScOAJ4Km06MMR8eEya505d3mZzQNw9Zl1pW/DzKwMpYWFpB7AdcDfAg3AYknzIuLxisUuAP4UEQdJmgR8ETgnzXsmIsaUVZ+ZmbVemYehxgGrImJ1RLwG3Aqc3mSZ04HvpMd3ACdJUok1mZnZLigzLAYD6yqeN6RpVZeJiC3AC8CANG+EpMck/UzSu6ttQNJUSfWS6jds2NC+1ZuZ2XZlhkW1PYRo5TLPAsMi4kjgYuD7kvbZacGI2RExNiLGDho0aLcLNjOz6soMiwZgaMXzIcD65paR1BPYF9gUEa9GxEaAiFgCPAMcUmKtZmbWgjLDYjFwsKQRknoDk4B5TZaZB0xJj88C7ouIkDQonSBH0oHAwcDqEms1M7MWlHY1VERskTQNuIfi0tk5EbFS0lVAfUTMA24EvidpFbCJIlAATgCukrQF2Ap8OCI2lVWrmZm1rNT7LCJiPjC/ybQrKh6/ApxdZb0fAT8qszYzM2s938FtZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLKrW7DzPrXDx8sO0q71mYmVmWw8LMzLIcFmZmluWwMDOzLIeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8vy4Edm1iE88NIbm/cszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMskq9GkrSROAaoAdwQ0TMajK/D/Bd4GhgI3BORKxJ82YCFwBbgY9GxD1l1mpmXZevxNp9pe1ZSOoBXAe8DxgJTJY0ssliFwB/ioiDgP8LfDGtOxKYBBwOTAS+kdozM7MaKHPPYhywKiJWA0i6FTgdeLximdOBK9PjO4CvS1KafmtEvAr8RtKq1N5DJdZrZtbuuspejSKinIals4CJEXFhen4+8I6ImFaxzIq0TEN6/gzwDooAeTgibk7TbwTuiog7mmxjKjA1PT0UeKqUF1PdQOC5DtxeZ+HX3b1019cN3ee1HxARg3ILlblnoSrTmiZTc8u0Zl0iYjYwu+2l7T5J9RExthbbriW/7u6lu75u6N6vvZoyr4ZqAIZWPB8CrG9uGUk9gX2BTa1c18zMOkiZYbEYOFjSCEm9KU5Yz2uyzDxgSnp8FnBfFMfF5gGTJPWRNAI4GFhUYq1mZtaC0g5DRcQWSdOAeygunZ0TESslXQXUR8Q84Ebge+kE9iaKQCEtdzvFyfAtwL9GxNayat1FNTn81Qn4dXcv3fV1Q/d+7Tsp7QS3mZl1Hb6D28zMshwWZmaW5bDYBZImSnpK0ipJl9a6no4gaaik+yU9IWmlpI/VuqaOJKmHpMck/bTWtXQUSW+WdIekJ9Pn/s5a19QRJH08/RtfIekHkvrWuqbOwGHRRq3sxqQr2gJ8IiIOA44F/rWbvO5GHwOeqHURHewa4O6IeDswmm7w+iUNBj4KjI2IURQX50yqbVWdg8Oi7bZ3YxIRrwGN3Zh0aRHxbEQ8mh7/meKLY3Btq+oYkoYApwA31LqWjiJpH+AEiisWiYjXIuL52lbVYXoCe6Z7v/rhe7wAh8WuGAysq3jeQDf50mwkaThwJPBIbSvpMF8DPglsq3UhHehAYAPwH+nw2w2S3lTrosoWEb8D/h1YCzwLvBAR/1XbqjoHh0Xbtaorkq5K0l7Aj4B/i4gXa11P2SSdCvwxIpbUupYO1hM4CvhmRBwJbAa6/Pk5Sf0pjhSMAPYH3iTpvNpW1Tk4LNqu23ZFIqkXRVDcEhFza11PBzkOOE3SGopDjn8j6ebaltQhGoCGiGjce7yDIjy6ugnAbyJiQ0S8DswF3lXjmjoFh0XbtaYbky4ndR1/I/BERHy11vV0lIiYGRFDImI4xWd9X0R0+V+aEfF7YJ2kQ9Okk9hxeIGuai1wrKR+6d/8SXSDE/utUepIeV1Rc92Y1LisjnAccD6wXNLSNO2yiJhfw5qsXBcBt6QfRauBD9W4ntJFxCOS7gAepbgC8DHc7Qfg7j7MzKwVfBjKzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhXYKkkPSViufTJV3ZTm3fJOms9mgrs52zU++u95e4jQ55Ldb1OCysq3gVOFPSwFoXUin1UtxaFwD/EhHvKases13lsLCuYgvFzVMfbzqj6a9pSS+lv8dL+pmk2yX9WtIsSedKWiRpuaS3VTQzQdIv0nKnpvV7SPqypMWSfiXpnyvavV/S94HlVeqZnNpfIemLadoVwPHA9ZK+3GT5PSR9I42x8FNJ8xtfj6STUkd/yyXNkdSnsb1U1wpJs9PdyE3rmCXp8VT7v7ft7bbuxmFhXcl1wLmS9m3DOqMpxqqoo7hD/ZCIGEfRHflFFcsNB06k6Kr8+jQgzgUUvZIeAxwD/G9JI9Ly44BPRcQOY35I2h/4IvA3wBjgGElnRMRVQD1wbkTMaFLjmWn7dcCFwDtTW32Bm4BzIqKOokeGj6R1vh4Rx6QxGfYETm1Sx1uADwCHR8QRwOdb+4ZZ9+SwsC4j9YL7XYrBa1prcRqr41XgGaCxO+rlFF/QjW6PiG0R8TRF1xdvB04G/jF1f/IIMAA4OC2/KCJ+U2V7xwALU0d1W4BbKMaNaMnxwA/T9n8PNJ7TOJSi07tfp+ffqWjrPZIekbScIpgOb9Lmi8ArwA2SzgT+kqnBujmHhXU1X6P4xV859sIW0r/1dDimd8W8Vyseb6t4vo0d+05r2i9OUHRXf1FEjEl/RlSMfbC5mfqqdXGf09w6VaenPY5vAGelPY5vAzsMDZqCahxFL8JnAHfvQl3WjTgsrEuJiE3A7RSB0WgNcHR6fDrQaxeaPjudO3gbxcBAT1F0JvmR1HU7kg5pxQBBjwAnShqYTn5PBn6WWecB4H+l7e8HjE/TnwSGSzooPT8/tdUYDM+l8Ud2uvopTd83dQT5bxSHxMya5V5nrSv6CjCt4vm3gR9LWgT8N83/6m/JUxRfxPsBH46IVyTdQHGo6tG0x7KB4ld6syLiWUkzKQ4lCZgfET/ObPtHFF1lrwB+TRE4L6QaPgT8MA0Buhi4PiJelfRtikNpa9L0pvameE/6pjp2ujDArJJ7nTV7A5C0V0S8JGkAsAg4Lp2/MOsQ3rMwe2P4qaQ3U5xv+ZyDwjqa9yzMzCzLJ7jNzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyy/j8F2VDr6ed0lQAAAABJRU5ErkJggg==\n", "text/plain": [ "