diff --git a/exercises/Chapter5_part2.ipynb b/exercises/Chapter5_part2.ipynb new file mode 100644 index 0000000..3499ea4 --- /dev/null +++ b/exercises/Chapter5_part2.ipynb @@ -0,0 +1,347 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 2 Exercises" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "from scipy import stats\n", + "import arviz as az\n", + "import pymc3 as pm\n", + "np.random.seed(123)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 6\n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'y')" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAE+CAYAAAC+4QU4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeVzVZd7/8deXw77JooiIC4oiq4gIaoqZW1ribrY7jY1NztjUNHfN3HPf4/0bvWsmx9vSmpYpG1vUxlwq11xKMsVAEBRESUHABZR9h3Ou3x8m48Iuh8OBz/Px8FF8189hOby5rut7XZpSCiGEEEII0bYsTF2AEEIIIURnJCFLCCGEEMIIJGQJIYQQQhiBhCwhhBBCCCOQkCWEEEIIYQQSsoQQQgghjKBNQpamaR9omparadrJm7a5aZr2taZpZ3/6r2tb3EsIIYQQwhxobTFPlqZpUUApsF4pFfTTtr8C+UqpVzVNexlwVUq91NS1unfvrvr373/XNQkhhBBCGFt8fPxVpVSP+vZZtsUNlFKHNE3rf9vmGcC9P/3/P4FvgCZDVv/+/YmLi2uLsoQQQgghjErTtMyG9hlzTFZPpdQlgJ/+69HQgZqm/ULTtDhN0+Ly8vKMWJIQQgghRPvoEAPflVLvKqXClVLhPXrU2+ImhBBCCGFWjBmyrmia1gvgp//mGvFeQgghhBAdSpuMyWrAF8CTwKs//Xd7ay9UU1NDdnY2lZWVbVWb6OBsbW3x9vbGysrK1KUIIYQQrdImIUvTtA1cH+TeXdO0bOBPXA9Xn2ma9nPgAjCvtdfPzs7GycmJ/v37o2laW5QsOjClFNeuXSM7OxsfHx9TlyOEEEK0Sls9XfhwA7smtMX1KysrJWB1IZqm4e7ujjwEIYQQwpx1iIHvzSEBq2uRr7cQQghzZzYhSwghhBDCnEjIEkIIIYQwAmM+XSiEEEII0WJlRVXkXShBX2PAc2A3HLrZmLqkVpGQ1cVs27aNHTt2kJuby5IlS5g8ebKpSxJCCNHFXTlfTEbyVfIulJCXVUJ5UfUt+5172OHl241evi54D3HF2d3ORJW2jHQXGtmyZctYuXLlXV1j9+7d+Pn54evry6uvvtrk8e+88w6enp4MHTqUgQMHsn79+rp9M2fO5L333uPDDz9k06ZNRq2noeNa+nqEEEJ0TiX5lex9/xSb/xJH/O5MSvIr6ePvxph5g5j122HMfSmce+b64u7lQEbyNQ5+dJpP/uso8bszUAZl6vKbJC1ZbUgphVIKC4vWZdf6ztfr9SxZsoSvv/4ab29vRowYQXR0NAEBAQ1eJykpiWXLlvHMM89w7Ngxpk2bxhNPPHHLMcuXL2fJkiUtrrG59TR0nJ+fX4tfjxBCiM6lpkrP8b2ZJO69gALCp/Vn2OS+WNveGUt6+jgTOrEvSikKr5Rz7MvzHN12jqzUAiYuDMDRteN2JUpLVgusWrWKoKAggoKCWL16NQAZGRn4+/vz7LPPEhYWRlZWFitWrMDPz4+JEyeSlpZWd/7HH39MREQEoaGhLF68GL1eX+/5Nzt27Bi+vr4MGDAAa2trFixYwPbtjU+en5ycjJ+fHwA+Pj5YW1vX7VNK8dJLLzF16lTCwsJa/Dlobj0NHdea1yOEEKLz+PF4Lp/86ShxOzLoP7Q7jyyLJDJ6QL0B62aapuHq6cDkRYGMf3wIV84XsWn5Mc6f6LhzKppfS9aul+Fyctte0zMYpjbebRUfH8+6deuIjY1FKUVkZCTjxo3D1dWVtLQ01q1bx1tvvUV8fDwbN24kISGB2tpawsLCGD58OKmpqWzatInDhw9jZWXFs88+yyeffEJUVNQt598uJyeHPn361H3s7e1NbGxso7XeCFlKKdauXcuKFSvq9q1Zs4Z9+/ZRVFREeno6zzzzDABjx46lpKTkjmutXLmSiRMntrieho5rzesRQghh/pRSHN+TydFt5+jR14nJiwLx8nVp8XU0TSPgHi96DezG3vdPsfPvyQSP682YhwZjYdGx5lg0v5BlIt999x2zZs3CwcEBgNmzZxMTE0N0dDT9+vVj5MiRAMTExDBr1izs7e0BiI6OBmD//v3Ex8czYsQIACoqKvDw8CAqKuqW82+n1J19zo1N1JmVlUVJSQnTpk0jJyeHkJAQli1bVrd/6dKlLF269I7zYmJimvFZaH49DR3X0tcjhBDC/CmD4rt/nSXpYDaDRvRkwpP+6CzvrjPN1dOBuf8RzpGtP3LiQBZWtjpGzfJto4rbhvmFrCZanIylvnBww43gdUNDoePJJ5/klVdeuWV7RkbGHeffzNvb+5YuxOzsbLy8vBo8PikpiaioKA4cOEBBQQFBQUEcOXKE0aNHN3gONL8lq7n1NHRcS1+PEEII86avMbDvwxTS43MZel8f7pnri9ZGLU46KwvGzB9EbY2e43su4NHPmYFhHm1y7bYgY7KaKSoqim3btlFeXk5ZWRlbt25l7Nix9R63detWKioqKCkp4csvvwRgwoQJbN68mdzcXADy8/PJzMxs8r4jRozg7NmznD9/nurqajZu3FjXOjZhwgRycnJuOT45OZlhw4YB4OrqyiOPPMKOHTuavE9MTAyJiYl3/Ls5YDVVT3OOa+75QgghzF91RS1frj1Benwuo2YP5J55bRewbjZ2/mB6+jiz/5+p5F8sa/Prt5aErGYKCwtj4cKFREREEBkZyaJFi+rCzO3HPfTQQ4SGhjJnzpy6IBYQEMDy5cuZPHkyISEhTJo0iUuXLjV5X0tLS9auXcuUKVPw9/dn/vz5BAYGYjAYSE9Px83N7Zbjbw5ZANOnT2fnzp13+eqbrueGadOmcfHixQaPa+p8IYQQnUN1ZS3b/i+BS2cLmbjQn7DJ/Yw2PERnZcH9vwjC0tqCXe8kU1VRa5T7tJTWWDeYKYSHh6u4uLhbtqWmpuLv72+iijqmkydP8sEHH7Bq1SpTl2I08nUXQgjzZDAodv09icxT+Ux7Jpj+Id3b5b45ZwrYvjoRm74OfKBKuVhUgZeLHb+b4sfMYb2Nck9N0+KVUuH17ZOWLDMVFBTUqQOWEEII8/X9lnQykq8xdv6gdgtYAL0Hu+Ic2Z3KjFK8r9SggJzCCn6/JZltCTlNnt/WJGQJIYQQos2cisnhxL4sgsd7E3yvd7vf/+8X80ixqmVspSX9aq7HnIoaPa/tSWvizLZnfk8XCiGEEKJDyjqdz6ENZ+gb6MaYuaaZTuFiUQVX7KGyAnJ1hn9vL6xo91qkJUsIIYQQd63gchl73j2Ji6c9kxcFYaEzTcTwcrGjVoP99jVUWNy6vb1JyBJCCCHEXamqqGXHW0lY6DQeeDYEGzvTdZT9boofdla6W7bZWen43RS/dq9FuguFEEIIcVdiNp2hOK+Cmb8Nw7l7+7cY3ezGU4Sv7UnjYqHxny5sjIQsIYQQQrRaenwuaUcvE/5A/1atRWgMM4f1Nkmoup10FwohhBCiVUoLqvjmk9N49HMifFp/U5fT4UjIEkIIIUSLKYPiwPoU9LUGJj0ViM5EA907MvmMCCGEEKLFkg5mk5VawD1zB+HS097U5XRIErK6mG3btvH0008zY8YM9u7da+pyhBBCmKFrF0s5svVH+ge7EzjWy9TldFgSsoxs2bJlrFy58q6u8dRTT+Hh4UFQUFCzjn/nnXfw9PRk6NChDBw4kPXr19ftmzlzJu+99x4ffvghmzZtalU9u3fvxs/PD19fX1599dV6j6msrCQiIoKhQ4cSGBjIn/70p0a3CyGEMA/6WgP71qVgbadj/OP+Rlv0uTOQkNWGlFIYDIamD2zh+QsXLmT37t3Nvk5SUhLLli3jxIkTbNiwgRdeeOGOY5YvX86SJUtaXKNer2fJkiXs2rWLlJQUNmzYQEpKyh3H2djYcODAAU6cOEFiYiK7d+/m6NGjDW4XQghhHhL2XuBqVinjHxuCvbO1qcvp0CRktcCqVasICgoiKCiI1atXA5CRkYG/vz/PPvssYWFhZGVlsWLFCvz8/Jg4cSJpaf9eK+njjz8mIiKC0NBQFi9ejF6vr/f820VFReHm5tbsOpOTk/Hzuz7pmo+PD9bW//4hUErx0ksvMXXqVMLCwlr8OTh27Bi+vr4MGDAAa2trFixYwPbt2+84TtM0HB0dAaipqaGmpgZN0xrcLoQQouMryqsgblcGA8N64DO0h6nL6fDMbp6svxz7C6fzT7fpNYe4DeGliJcaPSY+Pp5169YRGxuLUorIyEjGjRuHq6sraWlprFu3jrfeeov4+Hg2btxIQkICtbW1hIWFMXz4cFJTU9m0aROHDx/GysqKZ599lk8++YSoqKhbzm8LN0KWUoq1a9eyYsWKun1r1qxh3759FBUVkZ6ezjPPPAPA2LFjKSkpueNaK1euZOLEiXUf5+Tk0KdPn7qPvb29iY2NrbcOvV7P8OHDSU9PZ8mSJURGRja6XQghRMellOLQxjQsLDTGzBts6nLMgtmFLFP57rvvmDVrFg4ODgDMnj2bmJgYoqOj6devHyNHjgQgJiaGWbNmYW9//UmL6OhoAPbv3098fDwjRowAoKKiAg8PD6Kiom45/25lZWVRUlLCtGnTyMnJISQkhGXLltXtX7p0KUuXLr3jvJiYmGZdXyl1x7aGWqJ0Oh2JiYkUFhYya9YsTp48SVBQUIPbhRBCdFw/Hs/jwql8xswbhKOrjanLMQtmF7KaanEylvrCxQ03gtcN9YUOpRRPPvkkr7zyyi3bMzIy7jj/biQlJREVFcWBAwcoKCggKCiII0eOMHr06EbPa25Llre39y1dmtnZ2Xh5Nf5kiYuLC/feey+7d+++JUw1tF0IIUTHUl1Zy3efnaF7H0eC7zX9TOrmQsZkNVNUVBTbtm2jvLycsrIytm7dytixY+s9buvWrVRUVFBSUsKXX34JwIQJE9i8eTO5ubkA5Ofnk5mZeVc1TZgwgZycnFu2JScnM2zYMABcXV155JFH2LFjR5PXiomJITEx8Y5/NwcsgBEjRnD27FnOnz9PdXU1GzdurGutu1leXh6FhYXA9Va7ffv2MWTIkAa3CyGE6LiOfXGesuJqxj3ih4VMOtps8plqprCwMBYuXEhERASRkZEsWrSoLszcftxDDz1EaGgoc+bMqQtiAQEBLF++nMmTJxMSEsKkSZO4dOlSs+798MMPM2rUKNLS0vD29ub999/HYDCQnp5+x4D4m0MWwPTp09m5c+ddvPJbWVpasnbtWqZMmYK/vz/z588nMDCwbv+0adO4ePEily5dYvz48YSEhDBixAgmTZrEgw8+2OB2IYQQHVPehRKSDmYROLY3nj7dTF2OWdEa6wYzhfDwcBUXF3fLttTUVPz9/U1UUcd08uRJPvjgA1atWmXqUoxGvu5CCGFayqD4/LV4iq9W8Miykdg6WJm6pA5H07R4pVR4ffuM3pKladrzmqad0jTtpKZpGzRNszX2PbuCoKCgTh2whBBCmF7qkUtcOV/MPXN8JWC1glFDlqZpvYGlQLhSKgjQAQuMeU8hhBBC3L3qylpit5/Dc4AzgyM9TV2OWWqPMVmWgJ2maZaAPXCxHe4phBBCiLuQsPcC5cXV3DN3kEwa3UpGDVlKqRxgJXABuAQUKaXuWJVY07RfaJoWp2laXF5enjFLEkIIIUQTSgsqSfz6Ar7hHngOkMHurWXs7kJXYAbgA3gBDpqmPXb7cUqpd5VS4Uqp8B49ZJp+IYQQwpRit5/DoBSjZg40dSlmzdjdhROB80qpPKVUDbAFaHxWTCGEEEKYTN6FEk4fvczQ+/rg3N3O1OWYNWOHrAvASE3T7LXrHboTgFQj31MIIYQQraCU4vDms9g6WjF8an9Tl2P2jD0mKxbYDBwHkn+637vGvKcQQgghWicj6So5ZwqJeNAHGzuzW3mvwzH6Z1Ap9SfgT8a+jxBCCCFaT6838P2WH3HpaU/A2MbXpBXNI8vqdDHbtm3j6aefZsaMGezde8eDnkIIIbqolJiLFF4pZ/QcX3SyPmGbkM+ikS1btoyVK1e2+vysrCzGjx+Pv78/gYGBvP76602e88477+Dp6cnQoUMZOHAg69evr9s3c+ZM3nvvPT788EM2bdrUqpp2796Nn58fvr6+vPrqq40eq9frGTZs2C3rExYWFjJ37lyGDBmCv78/R44caVUdQggh2kZNlZ4fdmbgNciF/sHupi6n05CQ1YaUUhgMhjY939LSkr/97W+kpqZy9OhR3nzzTVJSUhq9TlJSEsuWLePEiRNs2LCBF1544Y5jli9fzpIlS1pco16vZ8mSJezatYuUlBQ2bNjQaD2vv/76HesPPvfcc9x///2cPn2aEydOyPqEQghhYsnfZFNRXM3IGQNk4tE2JCGrBVatWkVQUBBBQUGsXr0agIyMDPz9/Xn22WcJCwsjKyuLFStW4Ofnx8SJE0lLS6s7/+OPPyYiIoLQ0FAWL16MXq+v9/yb9erVi7CwMACcnJzw9/cnJyen0TqTk5Px8/MDwMfHB2tr67p9Sileeuklpk6dWnfdljh27Bi+vr4MGDAAa2trFixYwPbt2+s9Njs7mx07drBo0aK6bcXFxRw6dIif//znAFhbW+Pi4tLiOoQQQrSNqvIaju/JpF+QO7185f24LUnIaqb4+HjWrVtHbGwsR48e5b333iMhIQGAtLQ0nnjiCRISErh69SobN24kISGBLVu28MMPPwCQmprKpk2bOHz4MImJieh0Oj755JM7zu/Xr1+DNWRkZJCQkEBkZGSjtd4IWUop1q5dy4oVK+r2rVmzhn379rF582befvvtuu1jx44lNDT0jn/79u275do5OTn06dOn7mNvb+8GQ99vfvMb/vrXv2Jh8e9vs3PnztGjRw9+9rOfMWzYMBYtWkRZWVmjr0cIIYTxJO7Loqq8lsjoAaYupc2UJySQ8fAj1Fw07Up+Zvd85uX//V+qUk+36TVt/Ifg+Yc/NHrMd999x6xZs3BwcABg9uzZxMTEEB0dTb9+/Rg5ciQAMTExzJo1C3t7ewCio6MB2L9/P/Hx8YwYMQKAiooKPDw8iIqKuuX8hpSWljJnzhxWr16Ns7Nzg8dlZWVRUlLCtGnTyMnJISQkhGXLltXtX7p0KUuXLr3jvJiYmEbvf4NS6o5t9TUtf/XVV3h4eDB8+HC++eabuu21tbUcP36cNWvWEBkZyXPPPcerr77Kn//852bdXwghRNspL64mcX8WvsM96NHXydTl3LXaggLyVq2i8F+bsfT0pObyZay8TPekpNmFLFOpL1zccCN43VBf6FBK8eSTT/LKK6/csj0jI+OO829XU1PDnDlzePTRR5k9e3ajxyYlJREVFcWBAwcoKCggKCiII0eOMHp04xPtjx07lpKSkju2r1y5kokTJ9Z97O3tfUuXZnZ2Nl71fAMfPnyYL774gp07d1JZWUlxcTGPPfYYK1euxNvbu641bu7cuU0OnhdCCGEcx/dkoq/WEzHdx9Sl3BVlMFC0dSu5r61EX1KC289+Ro9fLcGiid+vxi9MqQ71b/jw4ep2KSkpd2xrb/Hx8So4OFiVlZWp0tJSFRgYqI4fP67Onz+vAgMD7ziuvLxcFRcXK19fX/Xaa6+pU6dOKV9fX3XlyhWllFLXrl1TGRkZd5x/O4PBoB5//HH13HPP3bHvvvvuU9nZ2bdse+WVV9QLL7xQ9/GLL76o/vCHP9zty69TU1OjfHx81Llz51RVVZUKCQlRJ0+ebPScgwcPqgceeKDu4zFjxqjTp08rpZT605/+pF588cV6z+sIX3chhOisiq9VqL8vOaj2/9O832srTqep8488qlL8hqjzDz+iKk6ntev9gTjVQKaRMVnNFBYWxsKFC4mIiCAyMpJFixYxbNiweo976KGHCA0NZc6cOYwdOxaAgIAAli9fzuTJkwkJCWHSpElcunSpyfsePnyYjz76iAMHDtSNk9q5cycGg4H09HTc3NxuOT45OfmWuqZPn87OnTvv8tX/m6WlJWvXrmXKlCn4+/szf/58AgMD6/ZPmzaNi030ga9Zs4ZHH32UkJAQEhMT+UMTXbVCCCHaXtzODBSK8Af6m7qUVtGXlnLllVc5P3s21T/+SK8VK+j38UfY+g02dWl1NNVIN5gphIeHq7i4uFu2paamymP+tzl58iQffPABq1atMnUpRiNfdyGEMI7CK+V8+j+xBI/rzdiHOk4oaQ6lFMU7dpL7l79Qe/UqLvPn4/H8b9CZ6El1TdPilVLh9e2TMVlmKigoqFMHLCGEEMbzw47z6Cw1s1sEuio9nct/Xk55bCy2QUF4v/UmdsHBpi6rQRKyhBBCiC6k4HIZZ3+4Quikvtg7Wzd9QgegLy3l6ptvkf/RR1g4OOC57E+4zJuHptOZurRGScgSQgghupC4XRnorCwYNqmvqUtpklKK4q++Ivevr13vGpw7lx7P/wbL28Yjd1QSsoQQQoguouByGWePXSF0Yl/snDp2K1bl6dNc/vNyKuLjsQ0OxvvNtdiFhJi6rBaRkCWEEEJ0ETdasUI7cCtWbUEBV9esoWDjJnTOznj+v//BZe5cNAvzmxBBQpYQQgjRBRReKefssSsMndgxx2IpvZ7Czz4jb/Xr6EtKcH34YXr8+lcme2qwLUjIEkIIIbqAuJ0Z6Cw75lis8rg4Li9fQdXp09hHRNDzP/+zQ8131VoSsoQQQohOrvBKOWeOXWbohD4dqhWr5uJFrrz2GiW7dmPZqxe9V/8fTlOm1Ls8nTmSkCWEEEJ0cnG7fmrFmtzP1KUAYKio4No/3ufaP/4BQPclS3Bf9HMs7OxMXFnbkpAlhBBCdGKFV8o5E3uZkA7QiqWUomTXLq68tpLaS5dwnjYVjxdfxMrLy6R1GYuErC5m27Zt7Nixg9zcXJYsWcLkyZNNXZIQQggjqmvFMvFYrIrkZK787ytUJCRg4+9P79f+in14vavRdBrm9zykmVm2bBkrV668q2v079+f4OBgQkNDCW/GN+Q777yDp6cnQ4cOZeDAgaxfv75u38yZM3nvvff48MMP2bRpU6vqeeqpp/Dw8CAoKKjR43bv3o2fnx++vr68+uqrTW4XQgjRtoryKjhz7AqBY3vj0M3GJDXUXLnCxZdeImPefKqzsui1/M/4bP5Xpw9YICGrTSmlMBgMRjn/4MGDJCYmcvvi2fVJSkpi2bJlnDhxgg0bNvDCCy/ccczy5ctZsmRJq+pcuHAhu3fvbvQYvV7PkiVL2LVrFykpKWzYsIGUlJQGtwshhGh7CXsz0Sxg2OT2b8UylJeT9+ab/Hj/VIp37cb96acZuHvX9TmvOvhyOG1FQlYLrFq1iqCgIIKCgli9ejUAGRkZ+Pv78+yzzxIWFkZWVhYrVqzAz8+PiRMnkpaWVnf+xx9/TEREBKGhoSxevBi9Xl/v+XcrOTkZPz8/AHx8fLC2/ncfvFKKl156ialTpxIWFtaq60dFReHWxJIGx44dw9fXlwEDBmBtbc2CBQvYvn17g9uFEEK0rdKCSlKPXMJ/tBcOLu3XiqUMBgq3bePHqdO4umYtjlFRDNi5A4/fvoDO0bHd6ugIJGQ1U3x8POvWrSM2NpajR4/y3nvvkZCQAEBaWhpPPPEECQkJXL16lY0bN5KQkMCWLVv44YcfAEhNTWXTpk0cPnyYxMREdDodn3zyyR3n9+t355MfmqYxefJkhg8fzrvvvttkrTdCllKKtWvXsmLFirp9a9asYd++fWzevJm33367bvvYsWMJDQ2949++ffta9fnKycmhT58+dR97e3uTk5PT4HYhhBBtK+HrCygDhLVjK1ZZ7DEy5s7j0su/x9LDg36ffIz366ux9vZutxo6ErMb+B7z2RmuZpW26TW793Fk7PzGJz377rvvmDVrFg4ODgDMnj2bmJgYoqOj6devHyNHjrxeX0wMs2bNwt7eHoDo6GgA9u/fT3x8PCNGjACgoqICDw8PoqKibjm/PocPH8bLy4vc3FwmTZrEkCFDiIqKqvfYrKwsSkpKmDZtGjk5OYSEhLBs2bK6/UuXLmXp0qV3nBcTE9Po628ppdQd2zRNa3C7EEKItlNeXE1KzEX8Invi3N340yJUnT9P7t/+Rum+/Vj26oXXa3/F+YEHzHIpnLZkdiHLVOoLBzfcCF431BcalFI8+eSTvPLKK7dsz8jIuOP823n99Girh4cHs2bN4tixYw2GrKSkJKKiojhw4AAFBQUEBQVx5MgRRo8e3eg9xo4dS0lJyR3bV65cycSJExs9tz7e3t63dH1mZ2fj5eXV4HYhhBBt58T+C9TWGhh+f3+j3qc2P5+rb75FwaZNWFhb0+M3z+G2cCEWtrZGva+5MLuQ1VSLk7FERUWxcOFCXn75ZZRSbN26lY8++qjR42pra/nyyy9ZvHgxEyZMYMaMGTz//PN4eHiQn59fb6i5XVlZGQaDAScnJ8rKyti7dy///d//DcCECRNYv349vXv3rjs+OTmZYcOGAeDq6sojjzzCjh07mgxZbd2SNWLECM6ePcv58+fp3bs3Gzdu5NNPP8XPz6/e7UIIIdpGZVkNyd/k4DvcA5ee9ka5h6GykvyPPuLaO+9iqKjAZd5cevzqV1h2726U+5mrrt2O1wJhYWEsXLiQiIgIIiMjWbRoUV2Yuf24hx56iNDQUObMmcPYsWMBCAgIYPny5UyePJmQkBAmTZrEpUuXmrzvlStXGDNmDEOHDiUiIoIHHniA+++/H4PBQHp6+h0D0G8OWQDTp09n586dd/nqb/Xwww8zatQo0tLS8Pb25v3336/bN23aNC5evIilpSVr165lypQp+Pv7M3/+fAIDAxvcLoQQom0kHcympkpvlFYsZTBQtH07P06bRt7fVmE/YgQDvthOr2XLJGDVQ2usG8wUwsPD1e3TFKSmpuLv72+iijqmkydP8sEHH7Bq1SpTl2I08nUXQoiWqa6oZf1/fo/XIBem/TKkTa9d+t1hcleupOr0aWwDA/H43e9wGBnZpvcwR5qmxSul6p30y+y6C8V1QUFBnTpgCSGEaLmTh3KoKq9l+NT+bXbNypQUcleupOz7I1h5e+P1t5U4T53a5Qe1N4eELCHMmL7WQFxGSXQAACAASURBVFV5LdWVtdRU6uv+a2Wjw87ZGntna2zsLeUJTiG6gJpqPYn7LtAnwI2e/Z3v+nrVWVnkvf4GxV99hc7FhZ5/+D0uCxZgYW3a9Q/NiYQsIcxETZWeq1kl5F4oIe9CCbmZJRReLqOpHn8LCw07JytcezngObAbvQZ0o+eAbtjYyY+/EJ1J6uFLVJTUED71zvkWW6L22jWuvvV3Cj77DE2nw/0Xv8B90c/ROd99cDM2gzKQWZxJyrUUUq6l8Kj/o3g5mu4JdnmXFaIDq66s5XxiHmd+yCUrNR9luJ6o7J2t6dHPiYHDemDvbI21rQ4rW0us7SyxstZRXVVLRUk1FcU1lJdUU15UxdXsUuJ3ZlwPZRq49XKgj78bfpGedO/jKK1dQpgxfa2BhL2Z9BrYjV6+Lq27RmkZ+evWkb9uHYaqKlzmzKH7kiVY9fRo42rbxu2B6tS1U5zOP01ZTRkANjobRnuN7twhS9M0F+AfQBCggKeUUkdaeh2llPwS6EI62gMZ7UmvN5CZfI0zx66QkXwVfY0BJzdbQif0oZdvNzz6Obd6iYzqylquZBRz+cciLv1YRPI32ZzYn4Wrpz2DIz0ZPOL6xIXbEnJ4bU8aFwsr8HKx43dT/Jg5rHfTNxBCmMSZY1coLahi3CN+Lf5daaiqonDjRq6+/Q76ggKcpkyhx3PPYTPAx0jVttztgSrlWgqp+am3BCo/Vz8eHPAgge6BBLgHMMBlAFYWViatuz1asl4Hdiul5mqaZg20eNIOW1tbrl27hru7uwStLkApxbVr17DtYpPZ6fUGzsReJm5nBsVXK7FzsiJgdC8GRXjiOcC5Tb73rW0t6TPEjT5Drk/9UVlWQ3p8LmeOXSZ2+zlit5/Dprc9G0uLyKEWNMgprOD3W5IBJGgJ0QEZDIrjezJx93akX5B7s89TtbUUbd9O3to3qb10CYfRo+jx/PPYBQcbsdqmmWugqo9Rp3DQNM0ZOAEMUM28UX1TONTU1JCdnU1lZaURqhQdka2tLd7e3lhZdbwfmrZ2e7jy6OfE8Kn96R/sjoWu/Z7eKb5awZljlzm44zz2esjR6TliW8t5SwNo0NvFjsMv39du9Qghmic9Ppc9751k8qJABoX3bPJ4ZTBQsvdr8t54g+pz57ANCcHj+d/gMGpUO1R7q+Z0+fm5+hHgHlD3b6DLQCwtOs5oJ1NO4TAAyAPWaZo2FIgHnlNKld1W4C+AXwD07XvnQpZWVlb4+HScZksh2oJSih+P53Fka3pduBr70GD6BZmmxda5ux3h03xY8G0KQdU6IiotmVtmw2Wdge9tazhXUNHuNQkhGqfU9Vasbh52DAxrfOyUUorSb78l7403qEpJxXrgQHqveQOniRPb5T3nRqA6de1UXaiqL1BNHzC9LlB11Baq5jJ2yLIEwoBfK6ViNU17HXgZ+K+bD1JKvQu8C9dbsoxckxAmV1ZYxbcb0jh/4ird+zjywJIQk4Wr2/V0tSOxsIIkaz0B1TpGVlkyu8yGi7ZQeKXcaMt0CCFaLisln7wLJYx/fAgWFg2/f5QdjSXv9depSEjAqk8fvP7yKs4PPoim0xmlrua2UJlDl9/dMHbIygaylVKxP328meshS4guSSlF6veXOLw5HX2tgVGzBxI6oU+7dgs25XdT/Pj9lmQqavSctNFzylpPRK0V46qt2fDnWMIm92P4/f2wtDbOm7MQovnid2fi6GqDX6RnvfvLExLIe+MNyo8cxbJnTzyXLcNlzmy0NhyK0awWKrfO1ULVXEYNWUqpy5qmZWma5qeUSgMmACnGvKcQHVXx1QoOfnya7NMFeA1yYfxjQzpkq9CNwe11Txe62vHwFD8mDejO4c3pxO3M4Myxy4x9aDD9g2WtMiFM5WJ6IRfPFjJm3iB0lrf+oVZx8hR5b7xO2aEYdG5ueLz0Eq4PL8DiLh8oMigDGcUZtwxKbypQdbQxVO3J6GsXapoWyvUpHKyBc8DPlFIFDR1f38B3IczdhVPX2Pv+KQwGxejZvgSO8UJrpGm/I8tOK+DQhjQKLpcTcE8vxswfjJWNtGoJ0d6+WnuCKxnFPLFidN3PYGVaGnlr1lC6bz+6bt1wW/Rz3B55BAsHhxZfv75AlXotlfLacuDfgSrAreMOSm8PJl27UCmVCNR7cyE6uxuDUo9uP4e7lwNTnwmmW4+O13rVEt5+rjz0xwiOfXWe43syuZhexKSnAvDo1/Fngxais8jLKiHz5DUio32wstFRmXaGq2++ScnevVg4OdF96a9xe+IJdI6OzbpefV1+9QWq6IHRXTpQtZR8doQwkurKWg6sT+XH43n4hntw3+P+nabFR2dpwaiZA+nr78bX61L4/K/xREYPYNikvmbbQieEOTm+JxMrWx2DvKvI/s3zlOzejYWjI92ffRa3J59A161bg+c2t8tPAtXdk8+YEEZQmFvOrreTKbhUxujZvoRO6tMhnhxsa739XFnwXxF88/Fpjmz9kazUfCY9FYi9sywgK4SxFF4p58f4XAbqznLxoV9iYW+P+y+fwf3JJ9G53LqkTnO7/GQMlXEYfUxWS8mYLGHu8i6U8MUbiSilmLIoiD7+bqYuyehuPDUZs/EMto5WTPtlCD36Opm6LCE6ncq0NPa9fpis2t7cc+IVPB+eifvChehcXJrd5dfVx1C1tcbGZEnIEqINXUov5Ku1J7C2t2TGc8M65NODxpSbWcyut5OpLK1hwsIAfId3zIVlhTA3lSkp5L31FtcOxfH9yP/HAPciBv1yMKm1WY12+UmgMj4JWUK0g6yUfHa+nYSjqy3Rz4Xi5Na11l68oayoit3vnOTyuSLCp/Un4kEfGaclRCtVnDhB3t//Ttk336J3tCPunicprQ5kS/hr5FleBG4NVIHdf5rYs9sACVTtxKRPFwrRFZxLyGPP+ydx9XQgemlolx6T5NDNhpnPD+PbDWnE7czgWk4pE38WgLWtvN0I0ZQbY6jOffMVVh9tx+PkRUrtNL6KsuDAUEvmnBpMQe8fmRg07t8Te0qg6rDkqyLEXUqLvcz+f6bi0c+JB381FFuHzj+LcVN0VhaMf3wI7t6OHP7XWbb/XwIP/noodo5dN3wKcbvbB6WfyjuJLu4UD8SUE5AFRQ6w78HelD04htDeQxmZ6M255BKWPvUo7l7Nm5pBmJaELCHuwo/Hc9n/YQpeg12Z9stgaa25iaZpDL2vD87utuz5xym2vHac6UuH4uxuZ+rShGh3dwSqq9fX8iuvLUdTitHpljxxVIdXdjk17s5oz88n7PHFjLS/HqaqK2tZ//fv8RnaXQKWGZHfCEK0UtbpfPZ+cIqePt14YEkIVrKWX718hvYgemkoO95Kqgta8ktCdGZ6g77etfxuPOVnq7PFz82Pmf2nE3mqml5bjkBGNlb9+tJ9+dM4R0djYX1rq++pmItUldcy/P7+JnhForUkZAnRCrmZxez6ezIuHvYSsJrBa5ALs34bxpdrEtm68jgPLBlKr4ENT5YohLloKlDZ6GwY4jaEGb4z6sZQ9bfuRenWL8hf+QE1Fy9iM3gw7n9bifOUKWiWd/5a1tcYSNx3Ae8hrvT0kZUVzImELCFaqOByGV+uOYGtoxXRS0NlDFYzdfd2ZM7vhvPFG4l8sTqBqc8E0zfQ3dRlCdFsN3f5nbp6qt5AdfNM6YHdA28ZlK4vLCT/00/J+Ohj9AUF2A0bRs8//hHH8fc2Ollx6pFLlBdVM/FnAe3yOkXbkZAlRAuU5FfyxeuJaBpEPxeKg4uNqUsyK87d7Zj94nC+XJPIzr8nM/WXwfSToCU6oBuB6kaYaqjL7+YWqoae8qu5dIn8f66n8LPPMJSX4zhuHO6/eBr74cObrkNvIGFvJj19nPH2c23z1ymMS0KWEM1UVV7Dl28kUl1Ry8zfhuHi0bUmGm0r9s7WzPjNMLavTmDn35OY9kwI/YIkaAnTae4YquYEqptVpp0h/4P3KdqxE5TCeepU3J9ehK2fX7NrOxuXS/HVSsbMG9Qpl+bq7CRkCdEMBoNi7/unKMqtIPq5UHr0kSVj7oatgxUzfjOML15PZOfbSUxdHEz/4O6mLkt0AS0NVIHugfh082n2PFRKKcqP/cC19/9B2aEYNDs7XB95GPcnn8Sqd+8W1aoMivjdmbh5OcjPh5mSkCVEMxzZ+iMXTuVz76N+9JYm+zZh62BF9HOhfPF6IrveSWbqL4LpHyK/SETbqS9QpeanUlFbAbS+hao+qqaG4j17yf/gAypTUtC5udHjuaW4LFiApWvr3jPOn7hKwaUyJv08QFZNMFMSsoRoQtrRSyR+fYHgcb0JHNuyv0RF424ErS/f+CloPSMtWqJ1mj1tgu/MNp0pXV9aRuHmf5G/fj21Fy9h3b8/nv/zP3SbEY2FbeuX1lJKEb87A+cedviGyRqg5kpClhCNuHK+mIMfp9Hbz4V75g8ydTmd0o2gtX11IrvfPcn0Xw2V1kLRqJZ0+QW6X1/LryVdfs1Rk5ND/sefUPivf2EoLcU+PBzPP/4XjveOQ7OwuOvrZ6cWkJtZwvjHhmChu/vrCdOQkCVEA8oKq9j5dhL23ayZ8nQQOnmjMxobeyumLx3K1r8lsOOtJGa+MAyPfjIfkDD+GKqWqjhxgmsffkjJ3q8BcJ4yBbeFT2IXEtKm94nblYGDiw1+kZ5tel3RviRkCVEPfY2BXe8kU12pZ+5/hMqae+3AztGa6KWhbFkZz5dvnGDWb8Nw83IwdVmiHbXnGKqWULW1lHz9Nfn/XE9FYiIWTk64Pfkkbo89ipWXV5vf71J6IRfPFjJm3iB0VvLHnTmTkCVEPQ5vSefK+WKmLg7GvbcsAdNeHF1tmPGbULa8dpwvXk9g9u+G49xd1jrsjJobqGb6zjRal1+TNRYWUvCvf1HwyafUXr6MVd++9PzD7+k2ew46R+P9ARC/OxNbRysCxrR9gBPtS0KWELc5l5hH8sFsht7XhwHDepi6nC6nWw97op8LZevfjrP99URmvxiGQzeZ9NWcNSdQDXEbwizfWXUtVO0dqG5WdfYs+R9/QtH27ajKSuxHjsTzv/8bx3FRaDrjLqGVd6GEzJPXiIwegJWNLNdl7iRkCXGT4msVHFifikc/J0bNHmjqcros996OPPjroWxfnXi96/DFMGzs5O3KHJhDC1V9lF5P6cGD5H/8CeVHj6LZ2OA8/UHcHn8CW7/B7VZH/K4MrG11BN8rTzJ3BvKuJcRP9HoDX79/CoNBMXlRIDpLGQthSp4+3Zi6OIgda5PY/U4yD/5qqHxNOhhzDVQ3qy0ooGjLFgo+3UBNTg6WXr3o8dsXcJk7t9XzW7XWtYul/JiQR/i0/tjYy5qonUHH+U4XwsSOfXGOy+eKmbwokG49ZMmcjqBvgDvjHx/C/n+mcuCjVCYuDJClRUxEb9DXLY5cX6Cys7TDz9Wvw3T5NaUi+SQFn35K8c6dqKoq7CMi8HjpP3C67z40S9PUHL8rE0sbHSH3eZvk/qLtdczvfiHa2YVT1zi+5wIBY70YFN7T1OWImwwZ1YvSgipivziHo6sto2ZKN66xtTRQ3Zg2QWfRsccQGaqqKNm9m/xPPqUyKQnN3p5us2bi+sgj2A5uvy7B+hReKSc97gqhE/vK08ydiIQs0eWVFVWx78MU3LwcGDtPJhztiIZP7UdJQSXHd2fi5GpD0Dj5S7+tNBWobnT5mVugull1ZiYFmz6jaMsW9IWFWPv40PM//5NuM2egc+oY65DG787AwtKC0El9TV2KaEMSskSXppTim49PU12pZ+bzQVham88vjq5E0zTGLRhMeWEVhzaewb6bDQNC5cnPlupsXX6NUbW1lH7zDQUbNlJ2+DDodDhNmIDrwwuwHzmyQ3U7F1+tIC32CsHjemPvLK1YnYn5/eQI0YZOH7lERvI1xswbJBNfdnAWOgsmLwpi2/8l8PX7p5j1YpjMCt+Iztrl15SaS5co/NdmCj//nNorV7Ds2ZPuv/4VLnPnYdWzY64BeHxPJpoFDJvcz9SliDYmIUt0WcXXKoj57Cxeg1wIGS/dT+bAykbHA8+G8K9Xf2DnW0nMfXkEjq4yh1ZLAlVg90AC3AI6RaC6Qen1lB46ROGmzyg9dAiUwmHMGDz/64843nuvyQayN0dpQSWpRy7hP9pLvpc7oY77nSeEESmD4sD606BgwpP+aBYdp+tANM7e2ZoHlwzl89fi2fHW9eV3rG27zlvZzYHqxtQJp/NP1zuxZ2cMVDerycmh8PPPKfx8C7VXrqDr0R33p5/GZd5crL3N4w+n43svgAHCJstYrM6o67wzCXGTk4dyyEkr4N5H/WTZFjPk3tuRKYuC2PHmCfatS2Hq4uBOGZSbMyj95pnSO0uXX2NUdTUlBw5SuHnz9bFWgMOYMfT8wx9wum88mpX5zC9VVlRFyncXGTzSU96HOikJWaLLKbxSzvdb0ukb6CZrg5mxfkHujJk/iJhNZzmy9UdGz/E1dUl3pakWqhtdfrMHzb4+KL0Tt1DVp+rsWQo/30LRF1+gz8/H0tOT7s8+i8vsWVj1Ns/Z0RO+voCh1sDwKTIWq7OSkCW6FINBsf+fqegsLRj/mH+HesJItFzwvd4UXC4n4esLuPS0N5vQfCNQ3TxT+u2BaojbkC4bqG7Ql5ZSvHMnhZ9/TuWJJLCywmn8eFzmzMZhzBijryNoTGVFVZz6NofBkZ649JTJjzsrCVmiS0k6kMXlc0VM/FmADDLtBDRNY+z8QRTnVfDtp2m49LTDa1D7LoXSFL1Bz/mi86TkpzQYqLpyC9XtlMFAeWwshVu3UrL3a1RlJda+A/F46SW6zYjG0s3N1CW2iYSvL6CvNRA+tb+pSxFG1C4hS9M0HRAH5CilHmyPewpxu+KrFcR+cY5+we4MjpBZ3TuL61M7BLL5L/Hseuck834fjrO7aca3SJdf61VfuEDRtm0UbttG7cVLWDg50W3GDFxmz8I2JKRTtTpLK1bX0WjI0jTtz8BVpdTrP328AriilHqjhfd5DkgFZFIbYRJKKb7dcAY0jXEP+3WqN2wBNvZWTPtlMJv/Es/Ot5KZ/TvjP3HYnBYq6fJrnL64mOLduynatp2K48dB03AYPRqP3/4WpwkTsLC1NXWJRpGw9wJ6vSJ8Wn9TlyKMrKl3ofeBLcDrmqZZAAuAiJbcQNM0b+ABYAXwQmuKFOJupcfncuHU9UlHndw65xt3V+fq6cCURYF8tfYEB/6ZypSng9rsiUMJVG1H1dRQevgwRdu3U7r/AKq6GusBA+jxwgt0i56OlaenqUs0qrKiKk4eysEvoicuHtKK1dk1GrKUUhmapl3TNG0Y0BNIUEpda+E9VgP/ATS4QJSmab8AfgHQt6/MFSLaVmVZDTGbzuDRz4lgmXS0U+sb6M6o2b58/3k6P+zMIOJBnxZfo7ldfjfPQ9W/W3+zXHqmvSilqExOpuiLLyneuRN9fj46Fxdc5s2j28yZ2AYFdpnW5YS9FzDoFcOlFatLaM67wj+AhYAn8EFLLq5p2oNArlIqXtO0exs6Tin1LvAuQHh4uGrJPYRoypEt6VSW1TL910Ow6IRzKYlbhU7sw7WcUn746jzuvR0YOKzhpVSkhcq4qjMzKdqxg+IvvqQ6IwPN2hrH++6jW/R0HMeMQbPuWuv0SStW19OckLUV+H+AFfBIC69/DxCtado0wBZw1jTtY6XUYy28jhCtknOmgJTDlxg2qS89+jbYmCo6EU3TuPdRPwqvlLNvXQouHva493aUQNVOanJzKdm1i6IdO6lMSgJNwz4iAvenF+E0eTI6p677c5iwR1qxuhpNqaYbjjRNexsoVEq93OobXW/JerGppwvDw8NVXFxca2/TtF0vw+Vk411fdBi1Bh2bkp/AoCxYEPxPrHS1pi5JtBM9itRqG2JOLqTWopr44NdJsSymQrv+fmenNIZgTQBWBCobArCmP5bokJbO1tBX6Ck+U0ZxahnlFypAgU1Pa7r5O+Ls74iVs3SlllU78FHizxnknsaEgXtMXU7X4RkMU1816i00TYtXSoXXt6/J7/yfBryPBOa1dWFCGNPxixEUVrox3W+zBKxOTI/iPDWkUE2KVk0KNZymmgpbRU+/fxB96tcMOvsoQ4Z8RICyIgBrfLCSQHWX9FUGStOvB6vS8+VgACsXS7qPcsHZ3xGb7l2rK7Apxy9GYFA6hveONXUpoh01NYVDAPAVsFUpdfZubqSU+gb45m6u0SaMnGhFx1CUV8Hx/4nFd3h3+j79lqnLEW2keV1+Acx2D6jr8qtItuHQp5aEuXzEqFkDTfwKzJu+tIzSb76hZM9uSr89hKquxrJXL9wWLsB52jRsAwO6zAD2lijJr+Tkfx/B/x5PXB7/xNTliHbU1NOFKcCAdqpFiDZzePNZNJ3GPXPNez27ruzmQHXq6vWn/NIK0u4YQzVn0Jzrgco9gP7O/e8cQxUFV7PKOL4nkx59nfAd3vBAeHEnQ1kZpd9+S/Gu3ZQeOoSqqsKyRw9c5s/Hedo07EKHollYmLrMDu2HHecBCH+g5U+7CvMmHeWi08k8eY3zJ64yatZAHF1lTixzcCNQ3byWX6sCVQOi5g8mP6eU/f9MwdXz+kB40TB9SQmlBw9SvGcvZd999+9gNW8ezvdPwS4sTIJVMxVeKef0kcsEj+stc/R1QRKyRKeirzEQs+kMLj3tGTqhj6nLEfVobgtVWz7lp7Oy4P7FwXz2vz+w8+9JzPv9CGwdrNrqJXUKtfn5lB44QPHXX1P2/RGoqcGyZ8/rLVaTJ2E3fLgEq1Y49uU5dJYaw2WNwi5JQpboVBL2XaAor4LpS4eis5RfCKbWnEDl7+bf6haqlnDoZsPUxcFs/dtxvv7gFA8sGdrl502ruXiRkn37KPl6H+Xx8WAwYNW7N26PPYbzlMnX1wyUYNVqV7NLORuXS9j9/bB3lgcBuiIJWaLTKMmvJH5nBgNCe9A3wN3U5XQ5tYba64HqWkqDXX7tFaga4jmgG2MfGsy3n6bxw1fniYzuWkNOlVJUnTlDyf79lO4/QOWpUwDYDPKl+zOLcZo4ERt/fxm83kZivziHtZ0lwybJSiZdlYQs0Wkc3pyOAu6ZJ4Pdja2+QHU6/zSV+kqgYwSqhgSO9SI3s5i4nRn06OvEgNAepi7JqFRNDeXxxyk5cD1Y1eTkgKZhN3QoPX77As6TJmHdv7+py+x0Lp8rIiPpKpHRA6RruguTkCU6hazT+fx4PJeI6T44u9uZupxOpb5B6bcHqiFuQ5g7eG6HC1T10TSNqAWDuZZdyr4PU5j3cjiung6mLqtN6YuKKD0UQ+nBg5TGxGAoKUGztsZh9Gjcn1mM0/jxWHbvbuoyO7XYL85h52RFyH2yXmpXJiFLmD2D3kDMprM4d7dl2GRplr8beoOec0XnmuzyM5dA1RBLKx33Lw7mX6/8wK63k5n7cjjWtub7dqiUojo9ndJvv6X020OUHz8Oej06d3ecJk/Cafx4HEaNwsKhc4XJjio7rYDs0wXcM9fXrL+vxN2Tr74weynfXaTgUhlTFwdjaWVev+xNqaVdfoHugfRz7md2gaohTm62TF4UxBerE9j/z1Tu/0WQWY1FMlRUUH7s2PVg9c231Fy8CIDNkCG4L1qE033jsQ0OloHr7UwpxZEt6Ti62hA0rrepyxEmJiFLmLWq8hpivzyP1yAXfEKl+6Mhna3Lr614+7kyarYv33+ezvE9mQy/v7+pS2qQUorq8xmUxRyi9FAM5T/8gKquRrOzw2HUKNwXL8ZxXBRWnp6mLrVLS4/PJTezhPue8Jc/+oSELGHe4nZlUllWw5h5g8yqFcKYmvuUX1cLVA0JndiH3MxiYrefw6OvM30C3ExdUh19SQnlsbGUfvcdZd8dpiY7GwBrHx9cH16Aw9go7EeEY2FjY+JKBYC+1sDRbT/i3tsBv5ESdoWELGHGivLKSTqQxZBRvejR18nU5ZhES6dN6Gxdfm1B0zTue9yf/Itl7H3/FPN+H45zd9M8PKFqa6k8eZLS77+n7PD3VCQmgl6Phb099pGRuP/8KRzGjsXaWwZTd0QnD+VQfLWSB38tc7CJ6yRkCbP1/ZYfsbC0YOSMrjHX0c2B6ka3X1p+WpcZQ2VMVjY6pv40EH73uyeZ/WIYltbG/7wppajOyKDs++8pO3KE8thjGEpKQNOwDQjAfdEiHO4ZjX1oKJq1TGbZkVVV1BK3IwPvIa707UCtocK0JGQJs5STVsC5hDwiowfg0K3zdZU0N1Dd6PKTQHX3XHraM/GpQHa+lcS3G9K47wnjTMpZc+kSZUdjKT96lLLYWGovXwbAqndvnO+/H4fRo7AfORJLV9c2v7doe9sScnhtTxoDLtUwssqKqgBnGbog6kjIEmbHYFB8t/ksjm42hE40//UJWzqGSgKV8fiEdCf8gf7E7cigp083gqLu/umwmiu5lP/wA+XHjlEeG0t1ZiYAOldX7EdG4hAZicPo0Vj16SO/nM3MtoQcfr8lGV2VnuFVtqRY1XIg5izWPWyZOUyeLBQSsoQZSjt6iatZpUz6eUC7dOm0pRuB6uan/M4UnJExVB1IxAM+5GWWELPpDN29HfEc0K1F59dcvkx5XPz1UHXsGNUZGQBYODpiHx6O6yMPYz9yJDaDBsn0CmbutT1pVNToub/SCg2Isa2lokbx2p40CVkCkJAlzExNlZ7Y7efo6ePMoPCepi6nUbWG2jsn9rypy8/e0p4hbkMkUHUwmoXGxJ8FXB+f9U4y8/4wosEuaaUUNRcuUB4Xdz1YxcVRk5UF/DtUucyfj31EBLb+Q9B08rXtTC4WVtBdrxFUrSPOppZinarbLgRIyBJm5sT+LMqKqpnydMeaOLK+FioZQ2W+bB2smPpMCJ//NY49751kxm+GobO0QNXUUHk6jYrj8ZTHH6f8gCdXkwAAIABJREFU+HH0V68CoHNxwX5EOG6PPYpdeDi2QyRUdXZe3ewYlaOnSoOjNrX/3u4iS3uJ6yRkCbNRXlzN8b2Z+AztTi9fF5PVIYGqa+ju7ci42X3ZvzGDfX/czJCLO6hITkZVXG+lsPL2xvGe0dgNC8N+eBjWvr4dKvgL4/vVkN7kZ+ZwwLaayp96fu2sdPxuip9pCxMdRpcLWVXnz1N7JRfboCB0jrKOlzmJ25lBbbWBUbMGtts9m9vlJ4HK/KmaGirPnKHixAkqTyRRceIEKiODPgNnk84E7PV9GDTXD/uwYdiFhWHVs2N3Vwvj0tcY0B/Px9LFmivdQCuqwMvFjt9N8ZPxWKJOlwtZRVu3ce3dd0HTsPEdiG1wCHYhIdiFBF8fiGplZeoSRT0Kc8s5dSiHgDFeuHoaJxy3OFB1D6SfkwQqc6SUoiYri4qkZCqTk6lITqby1ClUVRUAuu7dsQsJodvs2fQZOpSvYyw4lTmRIY8Px7mLTnwrbnXiQBZFeRVM//VQFge6m7oc0UF1uZDl/tTPsB8RTkVSEhVJSZQePEjRli0AaNbW2AwZgl1QELZBQdgGBWIzYACaZZf7NHU4R7edw+L/t3ff0XHe973n37+paEQH0UFUgr2JVCPVLVkSJZFSJFnxtVOcrJPYseNyvY6TnN27Z+/GcRR77ZM4yfWV10muJUuyRRVbcqhCdZESewEJggQBEoUE0ftg2m//GJAiJUqESAweAPN5nTNngBlgnq/0gDOf51e9LtasL5+U1/tgoKrrrqOhp+Gju/wUqGYsay3h9nZG6+oI1B0kcOAAgQMHiPT3A2D8fpIWLSLroYdIXr6M5OXL8RQVndf1d/uiIE/+7XZe+Nd9PPhXa0hO08KgiWy4f4wdLzRTvjSHMgUs+RjGWut0DedZvXq13bFjx5Qdz1pLqK0t1kVQd5DA+BVtdGQEiL0B+xfUkrRo0dmbv6YGl1ZfnjKnmvp56ns7WXNXBVfeVfGJf3+iLVRn9vFTl9/MZaNRQidOEDh0iMDBQ7H7ujoivb2xH/B48FdXk7x0CUlLl5K8dCn+6uoJtWCfPj7Apod3UVCVwT1fXY7LreUXEtUr/36Qhvc6+N3/4yoy81OcLkccZozZaa1dfaHnEr6JxhiDr6QEX0kJGevXA7E36mBzcyxwHTxE4OBBBp77NX2/eDz2Sx4P/qoqkhYswL9wAUkLFuKvna8VmuPAWss7Tx0lOd03oYVHL9ZCdW6X3+LcxSzKWaQWqhkqOjrK2JEjBOrrGas/TODwYcbq64kOD8d+wOPBX1ND2i03k7x4MUmLF+Ovrb3kzZTnzkvnxv9Syyv/foh3nmpk3YM1k/hfIzNFR/MA9VtPsfK2MgUsuaiED1kXYlwu/JWV+CsrydiwARi/Qm5pIVBXR6D+MIH6Qwy/8w79zz579vc8c+fir60lqXY+/tpa/DU1+Cor1ep1GZr3dXHyaD83fLYWX9L5f66X0uVXnl6Oy6gFYiY5+2+voYGxI0cYazjC2OHDsZXTo1EAXCkp+GtrydiwgaRFC/EvXBiXFucF1xTSeWKQvVtayC1LY8HVhZP6+jK92ajlzScaSE73sfqOcqfLkRlAIWuCjMuFb948fPPmkX7nnWcfD3d3EzhUz1hDA2OHDxNoaKBn2zZsKBT7AbcbX1kZ/poa/NXV+Guq8VVV4SsvV/i6iGgkytanG8nMT6Hmmjwaehs0y28Ws9EoofaTjB09QvDoUcaONjJ29ChjjY1nl00A8JaW4p8/n/Q77oh15S9YgLekZMpWT7/2/mq624d47eeHySpIJb88fUqOK85r2N5BR9MAN//eQnzJ+viUi0v4MVnxYEMhgs3Nsavuo0fPXn0HW1rOXnnjdse6Kaur8VdW4quowFdRjr+iAnemc2tATQdnWqi2v3aEvs3JHF79CluTXrzgGCp1+c080UCA4IkTBI81MXas8f37pubzwpQnLy92UVJdTdL8+fjnz8dfVYUr1fmlV0aHgvzyb3cQjVoe+M7qWblJuZxvbDTMY/9tG2mZfu7/9mqMS2uiSYzGZE0x4/XGWq5qzh+zEQ0ECDY1MdZ4jOCxxtiV+rFGht54A860fAHu7Gx85eWx23jrma98Hr6yMlwps2sMwJlAVdc1vrBnz0EaehoIhcI8tPtvCMzpYqC4jQdyHzg7MF1dftOfDYUItbfHwlTzcYLNzbFbUxOhkyfhnIs7b1ERvspKUteswVdVFWvxrarCnfHJ9gycSslpPu74s6Vs+vudbP7JATZ8PbYivMxe7z7TyOhAkPVfWqaAJROmlqxpwIbDhFpbGWtqInisiWBzU+zr48eJdHad97PuvFx8pWX4SkvxlpXiKyvDW1yCt6QYT17etF5x+qMC1QfHUC3KWURp03J6X/Vxz9eWU7pAU6Sno8jgIKHWVoItLYRaWgm2jt+fOEGorQ0ikbM/60pNjbXWnrl4KC/HXxn7fiZfOBzZ0cGLj9Sx6Loibvxs7bT+9yeXrqNpgF/9/Q6W3VjCdZ+Z73Q5Ms2oJWuaMx7P2Q8ebrrpvOciQ8OEThwnePw4weMnCLacIHT8BMPbthE+Z9A9xJab8BYV4S0pid2fuRUX4S0sxDN37pTtpXY5s/zGRsP8/ImtlC6ao4DlEGstkd5eQu0nCZ86Sai9nVBbW6x1qq2NUFs70fF1ps5wZ2TgLSkhecli0u+8A1/ZeAtsaSnu3NxZGUBqVufT1TLErs3HySlKY9lNJU6XJJMsGony2mP1pKb7uOqeSqfLkRlGIWuac6el4h5fn+uDooEAodZWQm1tBFtbCbW2xT4IW1sJ7N9PpK/vAy/mxpOXh7egAE9BQew+Px/P3Dy8+fl45s7FM3curqSkT1RjOBqmsa/x/UHpH2ihSvGksDBnIQ/UTqzLb89LJwgMh7h6g97Q4sGGQoS7uwl3dhLu6CDU0UG44zThjlOEOk4TPnWK0MmTZ1c/P8MkJ8cCe3ExKStWxAJ8aRm+0hK8JSW40xNzAPjVGyrpOTnMW788QlZ+CqWLsp0uSSbRvldb6WoZ4vYvLtFgd/nE9Bczg7mSkmLjV6qrL/h8dHiY0MmTsVtbO6FTJwmf6iB06hRj9fUMvfYaNhD48Oump+PJzcWTlxe7z83FnZuDJycXk5XBSd8oRzjNgcgJDgw1fGygOjPLb6JjqEYGgux5pYXqK+Yyd15ifmhfChsKEe7tJdLbS6Snh3BXN+HuLiLd3bGvu7pit9OnifT0nDcmCgCPZzxsF+BfuIC0m2/GW1CAt6gQT2Eh3sJC3FlZs7I16nIZl+HWLyziqb/fyeZHDnD/t1dr/aRZYrAnwLu/bqJ8aQ6VK/OcLkdmIIWsWcyVmvqxIcxaS3RggPDp04ROnyZ8OtayEe7sJNzVRajzNAN7dhLt6sY1Fjrvd6vHb3f6XITnJOPOyic5J5/UvEI8GZm4M5JwZ3TjytjLcHoz7vQ5uOakj9/PwZWaesEP7B2/bSYSiiZks7y1FhsIEBkcJDo0RKS/n+jgIJH+ASID/UQHBoj09Y3f+s9+He7tJTowcOEX9Xrx5OTgycnBm59P8pIlsRbLvLyz996CfNw5OVO2BMJs5EvysP5Ly/jld3fw/D/v4/5vX4E/RfugznRvPN4A1nLdQ/N1gSGXRCErgRljcGdkxG5VFeOD0i0Hu7s42NNJQ8+RWAuVtWRH57DCV8lidwnVNo/ScDqZI4ZoXz+R3l7Cfb1EevsI7N1HZGAg9qH/cZMqXC5cKSm40tJwpabiSk1lLG0uB3wbmOdtY/THrzGWnIwrNQWTlIwryY/xJ529N0l+XD4f5tyb1xvbZ9LtwXg9GI8nNgbN7Qbjwrhd4HbH3izPvGGe+7W1EI1ix++JRrHRKITD2GgUGw5DJIINh7GhUOwWDL5/PzZGdGwMOxbEBseIBgLY0VGiowGio6PYwCjRkVGiIyNEh4dj9yMjRIeGYqFqeBjC4Y89Z66UFNyZmbgyM/BkZuItLsKdlY07OwtPdvb7X+fm4snJwZWerg+HKZKem8wdf7qEZ//fPbz4SB3rv7xMW+/MYMf2dNK8r4tr7qsiPSfZ6XJkhopryDLGlAL/ARQAUeAn1tofxfOYH+eZ3W08vPkw7X2jFGUm861P17JxZbFT5ThmssdQXYiNRsdbYfrfb5EZGCQyOEB0cCh2PzxMdGiY1rYujrd0MmBqwBMhr24Tg2PdsWAyvofkTGeSknAlJWHGg6MrNRVXSgrenOxYcEobb+FLS8WdloYrbQ7ujHTc6em40t+/1wK201tRTRY3fLaWV39ez9tPHeW6BzUTbSYaGw3z5hMN5BSnsfyWi2/nJfJR4t2SFQa+aa3dZYyZA+w0xrxkrT0Y5+N+yDO72/jOpv2MhmLTytv6RvnOpv0AszpofShQdR/kcO9hxiKxQc2p3lQWZC/ggdoHWJwzPsvvE4yh+ijG5TrbSvZxzpyX1LlR/mDQz3Z/mH9c+2d8976lbFxZHGtBCgaxgUCslWh0dLy1aOz9FqRgkGgwCKHQeCtTOHYfCZ9thSISxUYjEImCjS0Ia619v7XNAi4T6zIzLjAmthaO2/N+C5jbA25XrMXs3JvHO97CNn7z+WKtbMnJsWCVlKSuuASyaF0R3e1D7NvSSubcFJbeqBmHM81bTzYw3B/k9j9ZilutkXIZ4hqyrLUngZPjXw8aYw4BxcCUh6yHNx9mNBQhN2IoCrvY548wGorw8ObDsyZkXaiF6nDP+4HqTAvVg7UPTmqguhxnzsunR30EgXf9YQIhzp4X43JhkpIgKQmt5y4zxdr7axjoCvDmEw3MyUmifGmu0yXJBDXv66J+6ymuuGOetkySyzZlY7KMMeXASuDdCzz3ReCLAGVlZXE5fntfbLuOmpCbtQEP9b4IQfP+4zPNxbr8zrRQPVj74LReKb29b5TCsKEm7ObNpBAB1/uPi8xUrvEZh09/fxcvPlLHfd9aRW7JHKfLkosIDIV49ef15BSnsWZ9hdPlyCwwJSHLGJMGPAV8zVr7oWlQ1tqfAD+B2Irv8aihKDOZtr5R2t1RDLHWrGZvlKLM6T+gcaJdfucGKqdbqCaqKDOZta0Rho1lpz983uMiM5kvycNdX17Or763g9/80z7u//Zq0rK0x+F09sYTDQSGQtz1leXaJkkmRdxDljHGSyxgPWqt3RTv432Ub326lu9s2s9JG8FiKYq46EgxfOvTtU6VdEHnBqq67joOdR86L1BNxy6/y/HVJSV0NbeyJSlIaHwSXLLXPe3Oi8ilSM30s/7Ly9j08C6e/+e93PvNVfiSNKl7OmrcdZoj2zu48u4K8krV6iiTI96zCw3wU+CQtfYH8TzWxZwZd/Xw5sN0DUapwMuD9zk7u/CTBKrp3OV3qay1uOsGcKd6OJ0DZiCxZ33K7JRbMofb/ngxL/zzPl76aR13/OlSLe0wzYwMBHntscPklc1h1e3znC5HZpF4X1KtBT4P7DfG7Bl/7K+stS/E+bgXtHFlMRtXFvPqo/Uc3XGaDcuLpuzYoWiIY33HPjJQpXpTWZi9kM/UfmbGdfldqub93XQ0DXDT5xbwp+um7lyITLXypblc/9B8Xv9FA68/dpgbP7dA65dNE9ZaXv/FYYKBMLf8wULNJpRJFe/ZhW8B0+6dpKAig4NvttN7aoTsotRJf/2LBaozLVSJFKg+yEYt7z7bSMbcZGqvKXC6HJG4W3JDCUN9Y+z87XFSMvwJuavBdHTwrXaO7e7kmvuqyClKc7ocmWUScnBAQWVsWu6ppv7LDlkTDVSzZQzVZDmys4PutmFu+6PFunKUhHHVPZWMDATZ8UIzqRk+ltygNbSc1NU6xJtPHqF0YRYrPxWfme2S2BIyZGXmp+BP9XCqsZ9FayfeTXWxMVTnzvJToPpokUiU955rIqc4jeor5jpdjsiUMcZw42drGR0M8frjDSSn+6haqX8DTggGwrz4yAH8yR4+9YeLY4sPi0yyhAxZxhgKKjI4daz/I39mIoEq0cZQTZb6d07S3znKnV9apjc2STgut4vb/ngxz/1wNy/99CBJX/VSPD/L6bISzhuPN9DXMcI9X1tJSrq2q5L4SMiQBVBQmcHxA90EhkN4ko0C1RQJhyLseKGZ/Ip0ypfmOF2OiCO8Pjfrv7ycTQ/v5IV/2c/Gr68kr0zLBkyV+q0nObztFGvuqqCkVgFX4ifhQtaZMVRH/YeBNL755F/zrnfLhwLVQ7UPnQ1UZellClSTpO6NdoZ6x7jlDxZpdpUktKRUL3d/dQVP/8MunvvRHjZ+c6UGXk+BnvZhXv/FYYprM1l9Z7nT5cgsl3Ah6x93/SM/q/sZnoiPL/A9MvsKeWidAtVUCAbC7PzPZkoWZOnqUQSYk53Ehq+vYNM/7OLZH+7hvm+uIjM/xemyZq1QMMLmRw7g9bu59QuLcWm4gsRZwoWsOyruoDa7lsU5i9nW3sE8143cs2al02UlhH1bWhkdDHHVBk1dFzkjIy+FDV9byTM/2MWzP9zNvd9cRXqutpWabDZqeeXfDtJzcpi7v7Kc1AxtcSTxl3BNNgtzFrK+cj3lGeUUVGbS0TSAjcZlu0Q5R2A4xO6XTlC+LJeCigynyxGZVrILU7nnL1YQGovw7A93M9Q75nRJs857v2micVcna3+nmrJFGg8qUyPhQta5CirTCQYi9JwcdrqUWW/3SycIBsJagFHkI+SWzOHur6xgdCjEcz/azXC/gtZkObK9gx0vNLPw2kKW31LqdDmSQBI7ZI23qHzcUg5y+UYGguzb0kLNFXPJLdHAXpGPkl+Rzl1fXs5g7xjP/GA3Q70Bp0ua8TqaBnjlPw5RWJ3BDZ+t1YQbmVIJHbIy5iaTlOZVyIqznb9tJhK2XHm3WrFELqaoJpN7vrKc4f4xnv7+Lga6Rp0uacYa6g3wwr/sIzXDxx1/shS3J6E/8sQBCf0XZ4yhoDKDU8cGnC5l1hrsCXDgzTYWXFOgWVMiE1RYncmGr61kbCTM09/fRV/HiNMlzTihsQgv/Mt+QsEId35pGclztOCoTL2EDlkQG5fV1zFCYDjkdCmz0vbfNAGwZn2Fw5WIzCz55els/MZKIuEoT39/Fz3tGjs6UZFQlN/+6z66Wga57Y8Wa/0xcYxClsZlxU3vqWHqt55kyfXFzMlOcrockRknt2QOG7++Cgw8/YNdPPliI2v/bgsVf/k8a/9uC8/sbnO6xGknEomy+ZEDtBzq5abPL6R8aa7TJUkCS/iQNbc8HeMydDSpy3CyvftcE26fmytuL3e6FJEZK7solXu/sYoQlvZNzfg7x7BAW98o39m0X0HrHNGo5ZV/O0TT3i6u+8x8Fl5b6HRJkuASPmR5/W5yilPVkjXJOk8M0rjrNCtuKdXmqyKXKTM/hV9mhul2W+4d9rFszA3AaCjCw5sPO1zd9GCt5fXHDnNkewdXb6xk2U0lTpckopAFsc2iO5oGiGpR0kmz7dlG/KkeVtxa5nQpIrPCsaFRHk8bo8kT5dOjPq4b9YCF9j7NPrTW8vZTRzn4VjtX3D5PrecybShkAYVVGYTGInS1DDpdyqzQfqSXE3U9rPr0PPzJCbdzk0hcFGUmEzLwdGqQvb4wV495WT/ipSQjsbfgsVHLO08dZe/LLSy7qUTbdsm0opAFlCzIBqDlUI/Dlcx81lq2PXOMlAwfS29Uc73IZPnWp2tJ9rqxBl5MDvFGUohFIQ+fH01K2EVLI+EoL//bQfa83MLSG0tY90CNFhuVaUUhC0hJ95FbmqaQNQmOH+jmZGM/a9ZX4PW5nS5HZNbYuLKY7963lOLMZIyB1gIv2bcWQX+IJ/92O631ifX+FQyEef7He2l4r4OrNlRy3WdqMC4FLJle1JczrnRhNntfaSE0FsHrVzi4FDZq2fbsMdJzkzSrRyQONq4sZuPK4vMe67m2lP/8H/t57kd7uHpjFStvK5v1rTkjA0F+80976Wod4ubfW8DCa4ucLknkgtSSNa50UTbRiKWtodfpUmasoztP0906xJV3V2r7CpEpkl2Yyv1/uZrKlXlsfbqR//wfBxgbDTtdVtz0dYzw1MM76T05zJ1/tlQBS6Y1fRKOK6zKwON1qcvwEkXCUbY9d4yc4lRq1uQ7XY5IQvElefj0/7aEtfdX07Svi8f/73dpmYXdh0d3nuaX393O2EiIDV9fqYVGZdpTyBrn8bopqsmk5eDse2OaCgffamegc5SrN1bh0rgIkSlnjGHFp8q477+uwu1x8dwP9/D6Y4cJBmZ+q1Y4FOH1xw6z+X8eIKswlQf/ag0FlRlOlyVyUQpZ5yhdlE3vqREGexJzps6lCgbCbH++iaKaTOYtyXG6HJGEVlCZwWf+5kqW31LKgTfbeOK/vzejh0H0dYzwq+/t5MAbbay4tYx7/+sq0nMSe9kKmTkUss5RulBLOVyKva+0MDoY4pr7qmb9gFuRmcDrc7PugRru/cYqMIZnfrCb1x6tZ2Qg6HRpE2ajloNvt/Pk325nqDfA+i8vY+3vVON262NLZg7NLjxHdlEqKRk+Wg71sGitBlNOxOhgkN0vnqByZd7ZzbZFZHooqsnkob+5km3PNrL/tTYatndwxe3zWH5zKZ5pvMTK6eMDvPF4Ax1NAxRWZ3DrFxZrk3mZkRSyzmGMoWxhNk37u4hGrcYWTcCOF5oJh6JcrVWWRaYlr9/NdQ/OZ8n1xbyzqZFtzxzjwOttXL2xivlr8qfV2lKjg0G2PXuMg2+3k5zm5ebfW8iCqwumVY0in4RC1geULsqmftspuloGmTsv3elyprX+zlEOvNHGwrWFZBWkOl2OiHyMrIJU1n9pGW2He3n7qaO8/LOD7H7xOMtuLmX+mnxHW7aCo2Hq3mpn52+bCQUiLL+llDXrK7Qtl8x4+gv+gDNb7Jw42KOQdRHv/foYLpfhyvUVTpciIhNUXJvFA3+5mobtHex+8Tiv/q96tj7dyOJ1RSy5oYS0LP+U1dJ3eoT9r7Zy6J2ThMYilC7MYt2D88ku1EWbzA4KWR9wdoudgz2svqPc6XKmrc6WQRre62DV7fNIzZy6N2URuXzGZai9qoD5V+bT1tDHvi0t7Nx8nN0vnqBieS4VK/KYtziHpDTvpB87HIrQWt9L3RttNB/oxuUy1KzOZ9nNJbqwlVlHIesCyhZls+flFoKBML4k/S/6IGtju977Uz2suq3M6XJE5BIZYyipzaKkNov+zlH2v97Kkfc6aNzdiTFQUJVB+bJc5i3OIbMg5ZJm9llr6esY4URdDycOdtPW0EckFCV5jpc1d5az+PpiUjN0oSazkxLEBZQuzGbX5hO0N/RRvkwrCn/Q8QPdtNb3su6BGvwpk3+lKyJTLyMvmXX317D2vmpOnxikeV8Xzfu72Lqpka2bGnG5DBn5KWQXppBVmErm3BQ8Xhcut8Hlid0bYKhvjIGuAINdowx0B+g7PcJIf2zpiMz8FBavK6JscQ4ltVm4vVqOQWY3hawLKKzKxON1ceJQj0LWB0QjUd556igZc5NZckPxxX9BRGYU4zLkl6eTX57OVfdUMtgToL2hl55TI/SeHKardYhjuzux9uNeBFIz/KTnJFGyIIvCqkzKFmWTnqtFRCWxxD1kGWNuB34EuIFHrLV/F+9jXi6310XR/CxatSjphxx8+yS9p0a440+XahNokQQwJzuJ2qsLz3ssHIow2B0gErZEI1Gikdi9jUJqpp852UlqpRIhziHLGOMGfgzcCrQC240xz1lrD8bzuJOhdGEWb/+qm8GegBbBGxccDfPer49RVJNJxXK18IkkKo/XrWVbRCYg3pcaVwJHrbXHrLVB4HFgQ5yPOSnO7MF3bHenw5VMHzs3H2d0MMTa+6u1fY6IiMhFxDtkFQMt53zfOv7YeYwxXzTG7DDG7OjsnB6hJqsgldzSNI7s6HC6lGlhoHuUvS+3MP+qfE2zFhERmYB4h6wLNXd8aLiktfYn1trV1trVeXl5cS5p4mpW59PRNMBA16jTpThu2zPHwMDVG6qcLkVERGRGiHfIagVKz/m+BGiP8zEnTfUVcwE4uvO0w5U4q6NpgCPbO1hxS6nGp4mIiExQvEPWdqDGGFNhjPEBDwHPxfmYkyY9N5n8ivSE7jK0UctbvzxC8hwvq26f53Q5IiIiM0ZcQ5a1Ngz8ObAZOAQ8aa2ti+cxJ1vN6ny6WoboPTXsdCmOOPzeKU4d6+eae6u0+r2IiMgnEPeFTKy1L1hr51trq6y1/0+8jzfZqq+YCwaO7Ei8LsOx0TDvbGokvyKdBR9YJ0dEREQ+nlaLu4jUTD/FNZkc3dGB/dgljmef7c83MToY5PqH5mNcWrJBRETkk1DImoDq1fn0nhqhu23I6VKmTHf7EPu2tLJoXZGWbBAREbkEClkTULUqD+MyHNmeGF2G1lrefOIIviQ3V2+odLocERGRGUkhawKS03yULsziSIJ0GTbu6qTtcC9X3VNJcprP6XJERERmJIWsCapZnc9gd4CO5gGnS4mr0FiEt391hNzSNBZf/6HF+UVERGSCFLImqGJFHi6P4egs7zLc+dtmhnrHuP4z83FpsLuIiMglU8iaIH+yh3mLcziys4NodHZ2GXa3DbH7xRPUXlVAYXWm0+WIiIjMaApZn0DNmnxG+oO0H+lzupRJF41atvyvenwpHtY+UO10OSIiIjOeQtYnUL4sF3+KhwOvtTpdyqTbt6WF080DXPeZGg12FxERmQQKWZ+A1+dm8XXFHNvTSX/nqNPlTJr+zlHefe4Y85bmULM63+lyREREZgWFrE9o2U0lGJdh75YWp0uZFNZaXnu0HuMy3PC7tRijwe4iIiKTQSHrE0rN9FOzJp9D75wkMBxyupzLVr/1JK31vVx7bxVzspOcLkdERGTWUMi6BCs+VUp4LELdm21Ol3JZhvvHePtXRymszmDxdVqHuJWWAAAJu0lEQVQTS0REZDIpZF2C3JI5lCzIYv+rrUTCUafLuSTWWt58vIFwMMpNn1ugDaBFREQmmULWJVpxaxnD/UGO7OhwupRLUr/1JI27O7ny7gqyClKdLkdERGTWUci6RGWLsskuSmXPSy0zbj/D3lPDvPF4A8W1may4tczpckRERGYlhaxLZIxh+S2ldLcN0Vrf63Q5ExYJRXnxp3V4vG4+9QeLtXWOiIhInChkXYbaKwtITvex5+UTTpcyYVufaaSrZYibf38haVl+p8sRERGZtRSyLoPb62LZjcWcqOuhu33I6XIuqnl/F3tfaWHpjSVULMt1uhwREZFZTSHrMi25vgSPz8W2Z445XcrHGu4fY8t/HCKnOJVrf6fK6XJERERmPYWsy5SU5mXNXRU07+vi2J5Op8u5oGjU8sq/HSQUiHDbHy3B43U7XZKIiMisp5A1CZbfUkpOcSpvPtFAMBB2upwPefuXR2g51Mt1D80nu0jLNYiIiEwFhaxJ4Ha7uPG/LGCob4z3ft3kdDnn2fdqK/tebWX5LaUsWlvkdDkiIiIJQyFrkhRUZrB4XRH7trTQeWLQ6XKA2ED3t55soHxZLtf+TrXT5YiIiCQUhaxJdPXGKpLSvLz2aD3RqLMLlHa2DPLiI3Xkls7htj/SelgiIiJTTSFrEiWleln3QA2njw9S94Zzm0cP9Y7x/I/34U/xsP5Ly/D6NdBdRERkqilkTbKaNfmULsxi6zONDPeNTfnxx0ZCPP/PewmOhln/5WWkZmrBUREREScoZE0yYwzX/24t0bDlpZ/VEQ5FpuzYw/1jPP39XfS0D3PbHy8mt2TOlB1bREREzqeQFQeZc1O46fMLaDvcx+b/WUckEo37MftOj7Dp4Z30dwW468+XU75UK7qLiIg4SSErTmqvKuCG351P874uXv7ZwbgOhO9sGWTTP+wiOBph49dXUrowO27HEhERkYnxOF3AbLbkhhKCYxG2bmrE63Nz0+cWYCZ5ll/7kV6e//E+fMke7vnGSrIKtNioiIjIdKCQFWerbptHKBBhxwvNeP1u1j1YgzGXH7SikSh7t7Ty7rPHSM9N4u6vrmBOdtIkVCwiIiKTQSFrClx5dwWhsQh7X2khErFcs7ESf4r3kl/vVFM/rz16mO7WIcqX5nDz7y8kOc03iRWLiIjI5VLImgLGGNbeXw0G9r7SwtGdHaxZX8GS64txeyY+LG5sNMy2Zxo58EYbqRl+7viTpVSsyJ2UljERERGZXMba+AzINsY8DNwNBIFG4A+ttX0X+73Vq1fbHTt2xKWm6aCzZZB3njpKa30vGXnJXHNvFZUr8z4yKEXCUU419nPiYA/1W08yOhhk6U0lXHVPJb4kZWQREREnGWN2WmtXX/C5OIas24At1tqwMeZ7ANbab1/s92Z7yAKw1nKiroe3nzpK78lh0nOTyMhLJjUribQsP2mZfqIRS8uhHlrrewmNRXC5DEXzM7nm3irmzkt3+j9BRERE+PiQFbemEGvti+d8uw24P17HmmmMMcxbkkPpwizqt57ieF03Q71j9LR3MzwQhPHcOycniflXFVC2KJuS2ix8yWq5EhERmSmm6lP7C8ATH/WkMeaLwBcBysrKpqgk57ncLhatK2LRuqKzj0UiUUb6g0QjlvTcJI23EhERmaEuK2QZY14GCi7w1F9ba58d/5m/BsLAox/1OtbanwA/gVh34eXUNNO53S4txSAiIjILXFbIstZ+6uOeN8b8PnAXcIuN1+AvERERkWkobt2FxpjbgW8DN1hrR+J1HBEREZHpKJ57F/4TMAd4yRizxxjzr3E8loiIiMi0Es/ZhdXxem0RERGR6S6eLVkiIiIiCUshS0RERCQOFLJERERE4kAhS0RERCQOFLJERERE4kAhS0RERCQOzHRbiN0Y0wkcj/NhcoGuOB9DPjmdl+lH52R60nmZfnROpqepOC/zrLV5F3pi2oWsqWCM2WGtXe10HXI+nZfpR+dketJ5mX50TqYnp8+LugtFRERE4kAhS0RERCQOEjVk/cTpAuSCdF6mH52T6UnnZfrROZmeHD0vCTkmS0RERCTeErUlS0RERCSuFLJERERE4iChQpYx5v8zxpw2xhxwuhaJMcaUGmNeNcYcMsbUGWP+wumaBIwxScaY94wxe8fPy//ldE0SY4xxG2N2G2N+43QtEmOMaTbG7DfG7DHG7HC6HgFjTKYx5lfGmPrxz5drHKkjkcZkGWOuB4aA/7DWLnG6HgFjTCFQaK3dZYyZA+wENlprDzpcWkIzxhgg1Vo7ZIzxAm8Bf2Gt3eZwaQnPGPMNYDWQbq29y+l6JBaygNXWWi1GOk0YY/4deNNa+4gxxgekWGv7prqOhGrJsta+AfQ4XYe8z1p70lq7a/zrQeAQUOxsVWJjhsa/9Y7fEueKbJoyxpQA64FHnK5FZLoyxqQD1wM/BbDWBp0IWJBgIUumN2NMObASeNfZSgTOdkvtAU4DL1lrdV6c90PgfweiThci57HAi8aYncaYLzpdjFAJdAI/G+9af8QYk+pEIQpZMi0YY9KAp4CvWWsHnK5HwFobsdauAEqAK40x6mJ3kDHmLuC0tXan07XIh6y11q4C7gC+PD40RZzjAVYB/2KtXQkMA3/pRCEKWeK48TE/TwGPWms3OV2PnG+8mf014HaHS0l0a4F7xsf/PA7cbIz5ubMlCYC1tn38/jTwNHClsxUlvFag9ZzW918RC11TTiFLHDU+wPqnwCFr7Q+crkdijDF5xpjM8a+TgU8B9c5Wldistd+x1pZYa8uBh4At1trPOVxWwjPGpI5P2mG8S+o2QDPYHWStPQW0GGNqxx+6BXBkMpXHiYM6xRjzC+BGINcY0wr8n9banzpbVcJbC3we2D8+/gfgr6y1LzhYk0Ah8O/GGDexi7EnrbVaMkDkw/KBp2PXi3iAx6y1/+lsSQJ8BXh0fGbhMeAPnSgioZZwEBEREZkq6i4UERERiQOFLBEREZE4UMgSERERiQOFLBEREZE4UMgSERERiQOFLBEREZE4UMgSERERiQOFLBGZtYwxa4wx+4wxSeMrc9dpD0YRmSpajFREZjVjzH8HkoBkYvuZfdfhkkQkQShkicisNr6txnYgAFxrrY04XJKIJAh1F4rIbJcNpAFziLVoiYhMCbVkicisZox5DngcqAAKrbV/7nBJIpIgPE4XICISL8aY3wPC1trHjDFu4B1jzM3W2i1O1yYis59askRERETiQGOyREREROJAIUtEREQkDhSyREREROJAIUtEREQkDhSyREREROJAIUtEREQkDhSyREREROLg/wcDYnM0WhxyDwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "x = np.arange(1,7)\n", + "y = np.random.randint(0, 10, size = 6)\n", + "\n", + "plt.figure(figsize=(10, 5))\n", + "order = [0, 1, 2, 5]\n", + "plt.plot(x, y, 'o')\n", + "for i in order:\n", + " x_n = np.linspace(x.min(), x.max(), 100)\n", + " coeffs = np.polyfit(x, y, deg=i)\n", + " ffit = np.polyval(coeffs, x_n)\n", + "\n", + " p = np.poly1d(coeffs)\n", + " yhat = p(x)\n", + " ybar = np.mean(y)\n", + " ssreg = np.sum((yhat-ybar)**2)\n", + " sstot = np.sum((y - ybar)**2)\n", + " r2 = ssreg / sstot\n", + "\n", + " plt.plot(x_n, ffit, label=f'order {i}, $R^2$= {r2:.2f}')\n", + "\n", + "plt.legend(loc=2)\n", + "plt.xlabel('x')\n", + "plt.ylabel('y', rotation=0)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected when we have a low order polynomial our model doesn't fit hte data well, with zero order polynomial we just predict the same value no matter what. With a high order polynomial the fit looks great, but this is only for data we have seen. Typically when a new data point is introduced the model will not do a good job. Additionally if there is noise in the data, the model will think the noise is meaningful, when in reality it's a mistake and means nothing." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 8" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Generate coin flips\n", + "coins = 30 \n", + "heads = 15\n", + "y_d = np.repeat([0, 1], [coins-heads, heads])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Per the advice in the chapter we will compute Bayes Factor using the SMC method" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sample initial stage: ...\n", + "Stage: 0 Beta: 0.334961 Steps: 25\n", + "100%|██████████| 2500/2500 [00:03<00:00, 788.98it/s]\n", + "Stage: 1 Beta: 1.000000 Steps: 4\n", + "100%|██████████| 2500/2500 [00:01<00:00, 1936.67it/s]\n", + "Sample initial stage: ...\n", + "Stage: 0 Beta: 0.146484 Steps: 25\n", + "100%|██████████| 2500/2500 [00:03<00:00, 700.04it/s]\n", + "Stage: 1 Beta: 1.000000 Steps: 4\n", + "100%|██████████| 2500/2500 [00:00<00:00, 3462.81it/s]\n" + ] + } + ], + "source": [ + "with pm.Model() as model_BF_0:\n", + " θ = pm.Beta('θ', 1, 1)\n", + " y = pm.Bernoulli('y', θ, observed=y_d)\n", + " trace_BF_0 = pm.sample(2500, step=pm.SMC())\n", + "\n", + "with pm.Model() as model_BF_1:\n", + " θ = pm.Beta('θ', .5, .5)\n", + " y = pm.Bernoulli('y', θ, observed=y_d)\n", + " trace_BF_1 = pm.sample(2500, step=pm.SMC())" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.5466473984070264" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model_BF_0.marginal_likelihood / model_BF_1.marginal_likelihood" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Referencing back to Figure 1.4, we can see that a $Beta(1,1)$ prior is flat, meaning all values are equally likely, but a prior $Beta(.5,.5)$ has larger values near 0 and 1. Intuitively it doesn't make much sense that our coin will be very biased heads, or very biased tails. Looking at the data, without inference, it seems like our coin is fair.\n", + "\n", + "The Bayes factor reflects this as the $Beta(1,1)$ model in the numerator yields a Bayes Factor of 1.4, indicating an \"anecotal\" preference. But in this case we're confirming that the Bayes Factor is sensitive to prior as the likelihood function of both models is the same." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 9\n", + "What does reduce sample size mean? Less samples from posterior? or less samples of data/coin flips?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 10" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.8700347100756574" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = range(0, 10)\n", + "q = stats.binom(10, 0.5)\n", + "q_pmf = q.pmf(x)\n", + "stats.entropy(q_pmf)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD4CAYAAAAejHvMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQqElEQVR4nO3dXYzld13H8c/XXQsqUaudG7stu+iqrE9Ux0Ul4kQKLNF0vShxMZBiMI3G+oTGVE2KWW98ig8XVdvIGiNoxeLFxCxWItQbU9wpRXRbG4cV2nExjC4+RJC68PVijnr4OdM9687sOZTXK5ns+f//v/+Z7yQnO++c+Z9zqrsDAAD8r8+Y9wAAALBoRDIAAAxEMgAADEQyAAAMRDIAAAz2z3uA0XXXXdcHDx6c9xgAADzDPfzww//Y3UvbHVu4SD548GDW1tbmPQYAAM9wVfWBnY653AIAAAYiGQAABjNFclUdq6rHq2q9qu7c5vjrq+rRqnpvVf1pVT136tjHq+o9k6/V3RweAAD2wiWvSa6qfUnuTvLSJBtJzlTVanc/OrXskSTL3f2Rqvq+JD+f5Dsnxz7a3S/Y5bkBAGDPzPJM8tEk6919rrufSnJfkuPTC7r7nd39kcnmQ0kO7O6YAABw9cwSydcneXJqe2OybyevS/K2qe1nV9VaVT1UVd+x3QlVdftkzdrm5uYMIwEAwN6Z5S3gapt9ve3CqlcnWU7yLVO7b+zu81X1vCTvqKq/6u73fdKddd+b5N4kWV5e3va+AQDgapnlmeSNJDdMbR9Icn5cVFU3J/mpJLd098f+e393n5/8ey7Jg0luuoJ5AQBgz80SyWeSHK6qQ1V1TZITST7pXSqq6qYk92QrkD80tf/aqnrW5PZ1SV6UZPoFfwA8Q6ysrGRlZWXeYwDsiktebtHdF6vqjiQPJNmX5FR3n62qk0nWuns1yS8keU6SP6iqJHmiu29J8vwk91TVJ7IV5D87vCsGAAAsnJk+lrq7Tyc5Pey7a+r2zTuc9+dJvupKBgQAgKvNJ+4BAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQD/D+trKxkZWVl3mOwoDw+4FPbTJFcVceq6vGqWq+qO7c5/vqqerSq3ltVf1pVz506dltV/e3k67bdHB4AAPbCJSO5qvYluTvJK5IcSfKqqjoyLHskyXJ3f3WS+5P8/OTcL0jyhiQvTHI0yRuq6trdGx8AAHbfLM8kH02y3t3nuvupJPclOT69oLvf2d0fmWw+lOTA5PbLk7y9uy9094eTvD3Jsd0ZHQAA9sYskXx9kientjcm+3byuiRvu5xzq+r2qlqrqrXNzc0ZRgIAgL0zSyTXNvt624VVr06ynOQXLufc7r63u5e7e3lpaWmGkQAAYO/MEskbSW6Y2j6Q5Py4qKpuTvJTSW7p7o9dzrkAALBIZonkM0kOV9WhqromyYkkq9MLquqmJPdkK5A/NHXogSQvq6prJy/Ye9lkHwAALKz9l1rQ3Rer6o5sxe2+JKe6+2xVnUyy1t2r2bq84jlJ/qCqkuSJ7r6luy9U1c9kK7ST5GR3X9iTnwQAAHbJJSM5Sbr7dJLTw767pm7f/DTnnkpy6v87IAAAXG0+cQ8AAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYzRXJVHauqx6tqvaru3Ob4i6vq3VV1sapuHY59vKreM/la3a3BAQBgr+y/1IKq2pfk7iQvTbKR5ExVrXb3o1PLnkjy2iQ/ts1dfLS7X7ALswIAwFVxyUhOcjTJenefS5Kqui/J8ST/E8nd/f7JsU/swYwAAHBVzXK5xfVJnpza3pjsm9Wzq2qtqh6qqu/YbkFV3T5Zs7a5uXkZdw0AALtvlkiubfb1ZXyPG7t7Ocl3JfmVqvri/3Nn3fd293J3Ly8tLV3GXQMAwO6bJZI3ktwwtX0gyflZv0F3n5/8ey7Jg0luuoz5AADgqpslks8kOVxVh6rqmiQnksz0LhVVdW1VPWty+7okL8rUtcwAALCILhnJ3X0xyR1JHkjyWJK3dPfZqjpZVbckSVV9fVVtJHllknuq6uzk9OcnWauqv0zyziQ/O7wrBgAALJxZ3t0i3X06yelh311Tt89k6zKM8bw/T/JVVzgjAABcVT5xDwAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGbgsKysrWVlZmfcYwKcQ/2/wqUgkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwGCmSK6qY1X1eFWtV9Wd2xx/cVW9u6ouVtWtw7HbqupvJ1+37dbgAACwVy4ZyVW1L8ndSV6R5EiSV1XVkWHZE0lem+R3h3O/IMkbkrwwydEkb6iqa698bAAA2DuzPJN8NMl6d5/r7qeS3Jfk+PSC7n5/d783ySeGc1+e5O3dfaG7P5zk7UmO7cLcAACwZ2aJ5OuTPDm1vTHZN4srORcAAOZilkiubfb1jPc/07lVdXtVrVXV2ubm5ox3DQAAe2OWSN5IcsPU9oEk52e8/5nO7e57u3u5u5eXlpZmvGsAANgbs0TymSSHq+pQVV2T5ESS1Rnv/4EkL6uqaycv2HvZZB8AACysS0Zyd19Mcke24vaxJG/p7rNVdbKqbkmSqvr6qtpI8sok91TV2cm5F5L8TLZC+0ySk5N9AACwsPbPsqi7Tyc5Pey7a+r2mWxdSrHduaeSnLqCGQEA4KryiXsAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMJgpkqvqWFU9XlXrVXXnNsefVVW/Pzn+rqo6ONl/sKo+WlXvmXz9xu6ODwAAu2//pRZU1b4kdyd5aZKNJGeqarW7H51a9rokH+7uL6mqE0l+Lsl3To69r7tfsMtzAwDAnpnlmeSjSda7+1x3P5XkviTHhzXHk/z25Pb9SV5SVbV7YwIAwNUzSyRfn+TJqe2Nyb5t13T3xST/kuQLJ8cOVdUjVfVnVfXN232Dqrq9qtaqam1zc/OyfgAAANhts0Tyds8I94xrPpjkxu6+Kcnrk/xuVX3u/1nYfW93L3f38tLS0gwjAQDA3pklkjeS3DC1fSDJ+Z3WVNX+JJ+X5EJ3f6y7/ylJuvvhJO9L8qVXOjQAAOylWSL5TJLDVXWoqq5JciLJ6rBmNcltk9u3JnlHd3dVLU1e+Jeqel6Sw0nO7c7oAACwNy757hbdfbGq7kjyQJJ9SU5199mqOplkrbtXk7wxye9U1XqSC9kK6SR5cZKTVXUxyceTfG93X9iLHwQAAHbLJSM5Sbr7dJLTw767pm7/R5JXbnPeW5O89QpnhLlaWVlJkjz44INznQOAZwa/Vz41+MQ9AAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZLa1srKSlZWVeY8BADzDLWpziGQAABiIZAAAGIhkAAAYiGQAABjMFMlVdayqHq+q9aq6c5vjz6qq358cf1dVHZw69hOT/Y9X1ct3b3QAANgbl4zkqtqX5O4kr0hyJMmrqurIsOx1ST7c3V+S5JeT/Nzk3CNJTiT5iiTHkvza5P4AAGBhzfJM8tEk6919rrufSnJfkuPDmuNJfnty+/4kL6mqmuy/r7s/1t1/l2R9cn8LZ1HffgQAgKuvuvvpF1TdmuRYd3/PZPs1SV7Y3XdMrfnryZqNyfb7krwwyU8neai73zTZ/8Ykb+vu+4fvcXuS25Pkxhtv/LoPfOADu/PTAQDADqrq4e5e3u7YLM8k1zb7xrLeac0s56a77+3u5e5eXlpammEkAADYO7NE8kaSG6a2DyQ5v9Oaqtqf5POSXJjxXAAAWCizRPKZJIer6lBVXZOtF+KtDmtWk9w2uX1rknf01nUcq0lOTN794lCSw0n+YndGBwCAvbH/Ugu6+2JV3ZHkgST7kpzq7rNVdTLJWnevJnljkt+pqvVsPYN8YnLu2ap6S5JHk1xM8v3d/fE9+lkAAGBXXPKFe1fb8vJyr62tzXsMAACe4a70hXsAAPBpRSQDAMBAJAMAwEAkAwDAYOFeuFdVm0nm9ZF71yX5xzl9bxabxwY78djg6Xh8sBOPjcXw3O7e9pPsFi6S56mq1nZ6hSOf3jw22InHBk/H44OdeGwsPpdbAADAQCQDAMBAJH+ye+c9AAvLY4OdeGzwdDw+2InHxoJzTTIAAAw8kwwAAAORDAAAA5GcpKqOVdXjVbVeVXfOex4WR1XdUFXvrKrHqupsVf3QvGdisVTVvqp6pKr+aN6zsDiq6vOr6v6q+pvJ/x/fOO+ZWAxV9SOT3yd/XVW/V1XPnvdMbO/TPpKral+Su5O8IsmRJK+qqiPznYoFcjHJj3b385N8Q5Lv9/hg8ENJHpv3ECycX03yx9395Um+Jh4jJKmq65P8YJLl7v7KJPuSnJjvVOzk0z6SkxxNst7d57r7qST3JTk+55lYEN39we5+9+T2v2XrF931852KRVFVB5J8W5LfnPcsLI6q+twkL07yxiTp7qe6+5/nOxULZH+Sz6qq/Uk+O8n5Oc/DDkTyVvA8ObW9ERHENqrqYJKbkrxrvpOwQH4lyY8n+cS8B2GhPC/JZpLfmlyK85tV9TnzHor56+6/T/KLSZ5I8sEk/9LdfzLfqdiJSE5qm33eF49PUlXPSfLWJD/c3f8673mYv6r69iQf6u6H5z0LC2d/kq9N8uvdfVOSf0/i9S6kqq7N1l+rDyX5oiSfU1Wvnu9U7EQkbz1zfMPU9oH40wdTquozsxXIb+7uP5z3PCyMFyW5paren63LtL61qt4035FYEBtJNrr7v//qdH+2ohluTvJ33b3Z3f+Z5A+TfNOcZ2IHIjk5k+RwVR2qqmuydQH96pxnYkFUVWXrusLHuvuX5j0Pi6O7f6K7D3T3wWz9v/GO7vaMEOnuf0jyZFV92WTXS5I8OseRWBxPJPmGqvrsye+Xl8SLOhfW/nkPMG/dfbGq7kjyQLZeZXqqu8/OeSwWx4uSvCbJX1XVeyb7frK7T89xJmDx/UCSN0+efDmX5LvnPA8LoLvfVVX3J3l3tt496ZH4eOqF5WOpAQBg4HILAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAY/BcUFwUUj0cTkwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "_, ax = plt.subplots(figsize=(12, 4))\n", + "ax.vlines(x, 0, q_pmf)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.7144085256537243" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = range(0, 10)\n", + "q = stats.binom(10, 0.25)\n", + "q_pmf = q.pmf(x)\n", + "\n", + "stats.entropy(q_pmf)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD4CAYAAAAejHvMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQsUlEQVR4nO3dX4yl913f8c+3u3VSiKCmnhu83uxGbGm2hcZl2FCihqPGSTYq8nLhiI0UZKpUq1aYAgFVBiSnWm74J2gv3NYW2QrxpyY4XIzQBjfCgRvksOM4TVgbi/GS2MMGZWEDVE0as8m3F3PSHv864znOzuw5xK+XNNrzPM/vOfMd6Wj3rWefM6e6OwAAwP/ztxY9AAAALBuRDAAAA5EMAAADkQwAAAORDAAAg4OLHmB0yy239JEjRxY9BgAAX+Eef/zxP+vule2OLV0kHzlyJOvr64seAwCAr3BV9cmdjrndAgAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhl2MZlMMplMFj0GAHADiWQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYzBXJVXWyqp6uqo2quneb4++uqier6mNV9dtV9eqZY1+oqo9Ov9b2cngAANgPB3dbUFUHktyf5M1JNpNcqKq17n5yZtkTSVa7+7NV9W+S/HSS754e+1x3v26P5wYAgH0zz5XkE0k2uvtSdz+f5KEkp2YXdPeHuvuz083Hkhza2zEBAODGmSeSb03y3Mz25nTfTt6V5AMz26+sqvWqeqyqvmu7E6rqzHTN+pUrV+YYCViUyWSSyWSy6DEAYF/tertFktpmX2+7sOqdSVaTfMfM7sPdfbmqXpPk0ar6eHc/84In634wyYNJsrq6uu1zAwDAjTLPleTNJLfNbB9KcnlcVFV3JPnxJHd29+e/tL+7L0//vJTkd5Lcfh3zAgDAvpsnki8kOVZVR6vqpiSnk7zgt1RU1e1JHshWIH96Zv/NVfWK6eNbkrwhyewb/gAAYOnsertFd1+rqnuSPJLkQJJz3X2xqs4mWe/utSQ/k+RVSX69qpLk2e6+M8lrkzxQVV/MVpD/5PBbMQAAYOnMc09yuvt8kvPDvvtmHt+xw3m/l+SbrmdAAAC40XziHgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAMRDIAAAxEMgAADEQyAAAM5orkqjpZVU9X1UZV3bvN8XdX1ZNV9bGq+u2qevXMsbur6o+mX3fv5fAAALAfdo3kqjqQ5P4kb0tyPMk7qur4sOyJJKvd/c1JHk7y09Nzvy7Je5K8PsmJJO+pqpv3bnwAANh781xJPpFko7svdffzSR5Kcmp2QXd/qLs/O918LMmh6eO3Jvlgd1/t7s8k+WCSk3szOgAA7I95IvnWJM/NbG9O9+3kXUk+8GWeCwAAC3dwjjW1zb7edmHVO5OsJvmOl3JuVZ1JciZJDh8+PMdIAACwf+a5kryZ5LaZ7UNJLo+LquqOJD+e5M7u/vxLObe7H+zu1e5eXVlZmXd2AADYF/NE8oUkx6rqaFXdlOR0krXZBVV1e5IHshXIn5459EiSt1TVzdM37L1lug8AAJbWrrdbdPe1qronW3F7IMm57r5YVWeTrHf3WpKfSfKqJL9eVUnybHff2d1Xq+onshXaSXK2u6/uy08CAAB7ZJ57ktPd55OcH/bdN/P4jhc591ySc1/ugAAAcKP5xD0AABiIZAAAGIhkAAAYiGS2NZlMMplMFj0GAMBCiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABjMFclVdbKqnq6qjaq6d5vjb6yqj1TVtaq6azj2har66PRrba8GBwCA/XJwtwVVdSDJ/UnenGQzyYWqWuvuJ2eWPZvke5P8yDZP8bnuft0ezAoAADfErpGc5ESSje6+lCRV9VCSU0n+byR39yemx764DzMCAMANNc/tFrcmeW5me3O6b16vrKr1qnqsqr7rJU0HAAALMM+V5NpmX7+E73G4uy9X1WuSPFpVH+/uZ17wDarOJDmTJIcPH34JTw0AAHtvnivJm0lum9k+lOTyvN+guy9P/7yU5HeS3L7Nmge7e7W7V1dWVuZ9agAA2BfzRPKFJMeq6mhV3ZTkdJK5fktFVd1cVa+YPr4lyRsycy8zwN9kk8kkk8lk0WMAsA92jeTuvpbkniSPJHkqyfu6+2JVna2qO5Okqr61qjaTvD3JA1V1cXr6a5OsV9X/SPKhJD85/FYMAABYOvPck5zuPp/k/LDvvpnHF7J1G8Z43u8l+abrnBEAAG4on7gHAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAg7kiuapOVtXTVbVRVfduc/yNVfWRqrpWVXcNx+6uqj+aft29V4MDAMB+2TWSq+pAkvuTvC3J8STvqKrjw7Jnk3xvkl8dzv26JO9J8vokJ5K8p6puvv6xAQBg/8xzJflEko3uvtTdzyd5KMmp2QXd/Ynu/liSLw7nvjXJB7v7and/JskHk5zcg7kBAGDfzBPJtyZ5bmZ7c7pvHnOdW1Vnqmq9qtavXLky51MDAMD+mCeSa5t9Pefzz3Vudz/Y3avdvbqysjLnUwMAwP6YJ5I3k9w2s30oyeU5n/96zgUAgIWYJ5IvJDlWVUer6qYkp5Oszfn8jyR5S1XdPH3D3lum+wAAYGntGsndfS3JPdmK26eSvK+7L1bV2aq6M0mq6lurajPJ25M8UFUXp+deTfIT2QrtC0nOTvcBAMDSOjjPou4+n+T8sO++mccXsnUrxXbnnkty7jpmBACAG8on7gEAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDsCcmk0kmk8mixwDYEyJ5yl/uAAB8iUgGAIDBXJFcVSer6umq2qiqe7c5/oqq+rXp8Q9X1ZHp/iNV9bmq+uj067/s7fgAALD3Du62oKoOJLk/yZuTbCa5UFVr3f3kzLJ3JflMd39DVZ1O8lNJvnt67Jnuft0ezw0AAPtmnivJJ5JsdPel7n4+yUNJTg1rTiX5xenjh5O8qapq78YEAIAbZ55IvjXJczPbm9N9267p7mtJ/jLJ35seO1pVT1TV71bVP9vuG1TVmapar6r1K1euvKQfAAAA9to8kbzdFeGec82nkhzu7tuTvDvJr1bV1/x/C7sf7O7V7l5dWVmZYyQAANg/80TyZpLbZrYPJbm805qqOpjka5Nc7e7Pd/efJ0l3P57kmSR//3qHBgCA/TRPJF9IcqyqjlbVTUlOJ1kb1qwluXv6+K4kj3Z3V9XK9I1/qarXJDmW5NLejA4AAPtj199u0d3XquqeJI8kOZDkXHdfrKqzSda7ey3Je5P8UlVtJLmarZBOkjcmOVtV15J8Icm/7u6r+/GDAADAXtk1kpOku88nOT/su2/m8f9O8vZtznt/kvdf54wAAHBD+cQ9AAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkANgHk8kkk8lk0WMAXyaRDAAAA5EMAAADkQwAAAORDAAAA5EMAAADkQwAAAORDAAAg7kiuapOVtXTVbVRVfduc/wVVfVr0+MfrqojM8d+dLr/6ap6696NDgAA+2PXSK6qA0nuT/K2JMeTvKOqjg/L3pXkM939DUl+PslPTc89nuR0kn+Y5GSS/zR9PgDgZcIHq/A30TxXkk8k2ejuS939fJKHkpwa1pxK8ovTxw8neVNV1XT/Q939+e7+4yQb0+cDAIClVd394guq7kpysrv/1XT7e5K8vrvvmVnzB9M1m9PtZ5K8Psm/T/JYd//ydP97k3ygux8evseZJGeS5PDhw9/yyU9+cm9+OgAA2EFVPd7dq9sdm+dKcm2zbyzrndbMc266+8HuXu3u1ZWVlTlGAgCA/TNPJG8muW1m+1CSyzutqaqDSb42ydU5zwUAgKUyTyRfSHKsqo5W1U3ZeiPe2rBmLcnd08d3JXm0t+7jWEtyevrbL44mOZbk9/dmdAAA2B8Hd1vQ3deq6p4kjyQ5kORcd1+sqrNJ1rt7Lcl7k/xSVW1k6wry6em5F6vqfUmeTHItyfd19xf26WcBAIA9sesb92601dXVXl9fX/QYAAB8hbveN+4BAMDLikgGAICBSAYAgIFIBgCAwdK9ca+qriRZ1Efu3ZLkzxb0vVluXhvsxGuDF+P1wU68NpbDq7t720+yW7pIXqSqWt/pHY68vHltsBOvDV6M1wc78dpYfm63AACAgUgGAICBSH6hBxc9AEvLa4OdeG3wYrw+2InXxpJzTzIAAAxcSQYAgIFIBgCAgUhOUlUnq+rpqtqoqnsXPQ/Lo6puq6oPVdVTVXWxqn5g0TOxXKrqQFU9UVW/uehZWB5V9Xer6uGq+sPp3x//dNEzsRyq6oem/578QVX9t6p65aJnYnsv+0iuqgNJ7k/ytiTHk7yjqo4vdiqWyLUkP9zdr03ybUm+z+uDwQ8keWrRQ7B0/mOS3+ruf5DkH8drhCRVdWuSf5tktbv/UZIDSU4vdip28rKP5CQnkmx096Xufj7JQ0lOLXgmlkR3f6q7PzJ9/D+z9Q/drYudimVRVYeS/Iskv7DoWVgeVfU1Sd6Y5L1J0t3Pd/dfLHYqlsjBJH+nqg4m+aoklxc8DzsQyVvB89zM9mZEENuoqiNJbk/y4cVOwhL5D0n+XZIvLnoQlsprklxJ8l+nt+L8QlV99aKHYvG6+0+S/GySZ5N8Kslfdvd/X+xU7EQkJ7XNPr8XjxeoqlcleX+SH+zuv1r0PCxeVX1nkk939+OLnoWlczDJP0nyn7v79iT/K4n3u5Cqujlb/1t9NMnXJ/nqqnrnYqdiJyJ568rxbTPbh+K/PphRVX87W4H8K939G4ueh6XxhiR3VtUnsnWb1j+vql9e7Egsic0km939pf91ejhb0Qx3JPnj7r7S3X+d5DeSfPuCZ2IHIjm5kORYVR2tqpuydQP92oJnYklUVWXrvsKnuvvnFj0Py6O7f7S7D3X3kWz9vfFod7siRLr7T5M8V1XfON31piRPLnAklsezSb6tqr5q+u/Lm+JNnUvr4KIHWLTuvlZV9yR5JFvvMj3X3RcXPBbL4w1JvifJx6vqo9N9P9bd5xc4E7D8vj/Jr0wvvlxK8i8XPA9LoLs/XFUPJ/lItn570hPx8dRLy8dSAwDAwO0WAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAw+D/dE/y0GU/QsAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "_, ax = plt.subplots(figsize=(12, 4))\n", + "ax.vlines(x, 0, q_pmf)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This answer intuitively makes sense, because of the lower bound of 0 on the overall distribution is less \"spread out\" with a binomial model where $p=.25$" + ] + } + ], + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}