{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bayesian Statistics Made Simple\n", "\n", "Code and exercises from my workshop on Bayesian statistics in Python.\n", "\n", "Copyright 2020 Allen Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Euro problem\n", "\n", "Here's a problem from David MacKay's book, [*Information Theory, Inference, and Learning Algorithms*](http://www.inference.org.uk/mackay/itila/p0.html), which is the book where I first learned about Bayesian statistics. MacKay writes:\n", "\n", "> A statistical statement appeared in The Guardian on\n", "Friday January 4, 2002:\n", ">\n", "> >\"When spun on edge 250 times, a Belgian one-euro coin came\n", "up heads 140 times and tails 110. ‘It looks very suspicious\n", "to me’, said Barry Blight, a statistics lecturer at the London\n", "School of Economics. ‘If the coin were unbiased the chance of\n", "getting a result as extreme as that would be less than 7%’.\"\n", ">\n", "> But [asks MacKay] do these data give evidence that the coin is biased rather than fair?\n", "\n", "To answer this question, we have to make some modeling choices.\n", "\n", "First, let's assume that if you spin a coin on edge, there is some probability that it will land heads up. I'll call that probability $x$.\n", "\n", "Second, let's assume that $x$ varies from one coin to the next, depending on how the coin is balanced and maybe some other factors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With these assumptions we can formulate MacKay's question as an inference problem: given the data --- 140 heads and 110 tails --- what do we think $x$ is for this coin?\n", "\n", "This formulation is similar to the 101 Bowls problem we saw in the previous notebook.\n", "\n", "But in the 101 Bowls problem, we are told that we choose a bowl at random, which implies that all bowls have the same prior probability.\n", "\n", "For the Euro problem, we have to think harder about the prior. What values of $x$ do you think are reasonable?\n", "\n", "It seems likely that many coins are \"fair\", meaning that the probability of heads is close to 50%. Do you think there are coins where $x$ is 75%? How about 90%?\n", "\n", "To be honest, I don't really know. To get started, I will assume that all values of $x$, from 0% to 100%, are equally likely. Then we'll come back and try another prior.\n", "\n", "Here's a uniform prior from 0 to 100." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "xs = np.linspace(0, 1, num=101)\n", "p = 1/101\n", "prior = pd.Series(p, index=xs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's what it looks like." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbBklEQVR4nO3dfbRddX3n8fcHELGKBkxwKA8GakBjlZpekJmpFW2xBB+Co7akdqBoF2UE265Wx9ihVWo7g+20ViqFBQpCnYK09SEusQxFLZYS5QYxEBGJKUoKlVA7WB9p8Dt/7B28PZx770n4nXu9yfu11l1n7/17OL9fbtb53P1w9k5VIUlSC3vM9wAkSbsOQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCrSD5gk30hy+HyPQ9oZhoo0B5LcleTbfWB8NcmlSZ4wrG5VPaGqNs/1GKUWDBVp7ry0qp4ArACOBs6eWphkr0fT+aNtL7VgqEhzrKr+EfgY8KNJKsmZSe4E7gTotz2tX35SksuTbE3y5SRnJ9mjL/vFJDckeUeSrwFvnacpSQ/zLxtpjiU5BDgR+ADwYuAk4LnAt4dU/xPgScDhwJOB/wvcC7ynL38ucCVwAPCYcY5bGkW895c0fknuAhYD24AHgI8CvwF8C/ipqvr4lLoFLAP+oS9/TlV9vi/7ZWB1VR2X5BeB36mqQ+dwKtKM3FOR5s5JVfU3UzckAbh7mvqLgb2BL0/Z9mXgoCnr07WV5oXnVKT5N93hgvuBfwOeOmXbocA/jtBWmheGivQDqqoeAq4Cfi/JvkmeCvw68L75HZk0PUNF+sH2euCbwGbg74A/By6Z1xFJM/BEvSSpGfdUJEnNGCqSpGYMFUlSM4aKJKmZ3frLj4sXL66lS5fO9zAkaUFZv379/VW1ZFjZbh0qS5cuZXJycr6HIUkLSpIvT1fm4S9JUjOGiiSpGUNFktSMoSJJamasoZLkhCR3JNmUZM2Q8iQ5ry/fkGTFlLJLktyX5LaBNvsnuTbJnf3rfgPlh/bPAX/D+GYmSRpmbKGSZE/gfGAlsBxYnWT5QLWVdA8jWgacDlwwpey9wAlDul4DXFdVy4Dr+vWp3kH3qFZJ0hwb557KMcCmqtpcVQ/SPfJ01UCdVcDl1VkHLEpyIEBVXQ98bUi/q4DL+uXL6B7FCkCSk+ju5rqx4TwkSSMaZ6gcxL9/Kt0W/v0T60atM+gpVXUvQP96AECSxwNvAs6ZqXGS05NMJpncunXrrJOQJI1unKGSIdsG77M/Sp1RnQO8o6q+MVOlqrqoqiaqamLJkqFfCJUk7aRxfqN+C3DIlPWDgXt2os6gryY5sKru7Q+V3ddvfy7wyiS/DywCvpfkO1X1rp2dgCRpx4xzT+UmYFmSw5LsDZwMrB2osxY4pb8K7Fjgge2HtmawFji1Xz4V+DBAVT2vqpZW1VLgj4H/aaBI0twaW6hU1TbgLOAa4HbgqqramOSMJGf01a6mO7G+CbgYeN329kmuAG4EjkyyJclr+6JzgeOT3Akc369Lkn4A7NaPE56YmChvKClJOybJ+qqaGFbmN+olSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNTPWUElyQpI7kmxKsmZIeZKc15dvSLJiStklSe5LcttAm/2TXJvkzv51v3778UnWJ7m1f33hOOcmSXqksYVKkj2B84GVwHJgdZLlA9VWAsv6n9OBC6aUvRc4YUjXa4DrqmoZcF2/DnA/8NKqehZwKvBnbWYiSRrVOPdUjgE2VdXmqnoQuBJYNVBnFXB5ddYBi5IcCFBV1wNfG9LvKuCyfvky4KS+/mer6p5++0ZgnySPbTkhSdLMxhkqBwF3T1nf0m/b0TqDnlJV9wL0rwcMqfMK4LNV9d3BgiSnJ5lMMrl169ZZ3kqStCPGGSoZsq12os6OvWnyTODtwC8PK6+qi6pqoqomlixZ8mjeSpI0YJyhsgU4ZMr6wcA9O1Fn0Fe3HyLrX+/bXpDkYOCDwClV9aWdHLckaSeNM1RuApYlOSzJ3sDJwNqBOmuBU/qrwI4FHth+aGsGa+lOxNO/fhggySLgo8Cbq+qGRnOQJO2AsYVKVW0DzgKuAW4HrqqqjUnOSHJGX+1qYDOwCbgYeN329kmuAG4EjkyyJclr+6JzgeOT3Akc36/Tv9fTgN9Kckv/M+x8iyRpTFL1qE5hLGgTExM1OTk538OQpAUlyfqqmhhW5jfqJUnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqZqRQSfKSJAaQJGlGowbFycCdSX4/yTNG7TzJCUnuSLIpyZoh5UlyXl++IcmKKWWXJLkvyW0DbfZPcm2SO/vX/aaUvbnv644kPzPqOCVJbYwUKlX1C8BzgC8Blya5McnpSfadrk2SPYHzgZXAcmB1kuUD1VYCy/qf04ELppS9FzhhSNdrgOuqahlwXb9O3/fJwDP7dn/aj0GSNEdGPqRVVV8H/gq4EjgQeDlwc5LXT9PkGGBTVW2uqgf7dqsG6qwCLq/OOmBRkgP797se+NqQflcBl/XLlwEnTdl+ZVV9t6r+AdjUj0GSNEdGPafysiQfBD4OPAY4pqpWAkcBb5im2UHA3VPWt/TbdrTOoKdU1b0A/esBO9JXv4c1mWRy69ats7yVJGlH7DVivVcC7+j3Hh5WVd9K8ppp2mTIttqJOqMaqa+qugi4CGBiYmJn30uSNMSoh7/uHQyUJG8HqKrrpmmzBThkyvrBwD07UWfQV7cfIutf73sUfUmSGho1VI4fsm3lLG1uApYlOSzJ3nQn0dcO1FkLnNJfBXYs8MD2Q1szWAuc2i+fCnx4yvaTkzw2yWF0J/8/M0tfkqSGZjz8leS/Aa8DfiTJhilF+wI3zNS2qrYlOQu4BtgTuKSqNiY5oy+/ELgaOJHupPq3gNOmvPcVwHHA4iRbgLdU1XuAc4GrkrwW+Arwqr6/jUmuAj4PbAPOrKqHRvpXkCQ1karpTyskeRKwH/C/6C/d7f1rVQ27MmtBmZiYqMnJyfkehiQtKEnWV9XEsLLZTtRXVd2V5Mwhne6/KwSLJKmd2ULlz4GXAOvprqSaeoVVAYePaVySpAVoxlCpqpf0r4fNzXAkSQvZbCfqV8xUXlU3tx2OJGkhm+3w1x/OUFbACxuORZK0wM12+OsFczUQSdLCN9vhrxdW1ceT/Jdh5VX1gfEMS5K0EM12+Ov5dDeRfOmQsgIMFUnSw2Y7/PWW/vW0mepJkgSj3/r+yf0TGm9Osj7JO5M8edyDkyQtLKPeUPJKYCvwCrrb4G8F3j+uQUmSFqZRn6eyf1W9bcr67yY5aQzjkSQtYKPuqXwiyclJ9uh/fhb46DgHJklaeGa7pPhf+f49v34deF9ftAfwDeAtYx3dD7BzPrKRz9/z9fkehiTtlOU//ETe8tJnNu93tqu/9m3+jpKkXdao51RIsh/d0xT32b5t8BHDu5NxJLwkLXQjhUqSXwJ+le6577cAxwI34r2/JElTjHqi/leBo4Ev9/cDew7dZcWSJD1s1FD5TlV9ByDJY6vqC8CR4xuWJGkhGvWcypYki4APAdcm+RfgnnENSpK0MI0UKlX18n7xrUk+ATwJ+OuxjUqStCDtyNVfK4CfoPveyg1V9eDYRiVJWpBGvaHkbwOXAU8GFgOXJjl7nAOTJC08o+6prAaeM+Vk/bnAzcDvjmtgkqSFZ9Srv+5iypcegccCX2o+GknSgjbbvb/+hO4cyneBjUmu7dePB/5u/MOTJC0ksx3+muxf1wMfnLL9k2MZjSRpQZvthpKXbV9OsjdwRL96R1X92zgHJklaeEa999dxdFd/3UV3G/xDkpy6O99QUpL0SKNe/fWHwIuq6g6AJEcAVwA/Pq6BSZIWnlGv/nrM9kABqKovAo8Zz5AkSQvVqKGyPsl7khzX/1xMd/J+RklOSHJHkk1J1gwpT5Lz+vIN/bf2Z2yb5KgkNya5NclHkjyx3/6YJJf1229P8uYR5yZJamTUUDkD2Aj8Ct1t8D/fb5tWkj2B84GVwHJgdZLlA9VW0j34axlwOnDBCG3fDaypqmfRXZH2xn77q4DH9tt/HPjlJEtHnJ8kqYFZz6kk2QNYX1U/CvzRDvR9DLCpqjb3/VwJrKILpO1WAZdXVQHrkixKciCwdIa2RwLbLxC4FrgG+C267888PslewOOABwEfIi9Jc2jWPZWq+h7wuSSH7mDfBwF3T1nf0m8bpc5MbW8DXtYvvwo4pF/+S+CbwL3AV4D/XVVfGxxUktOTTCaZ3LrV54xJUkujHv46kO4b9dclWbv9Z5Y2GbKtRqwzU9vXAGcmWQ/sS7dHAt2e0UPADwOHAb+R5PBHdFJ1UVVNVNXEkiVLZpmCJGlHjHpJ8Tk70fcWvr8XAd3z7Qcf7DVdnb2na9s/dfJF8PClzS/u6/w88Nf9lzLvS3IDMAFs3omxS5J2wox7Kkn2SfJrdIeZnk73HJW/3f4zS983AcuSHNZ/G/9kYHDvZi1wSn8V2LHAA1V170xtkxzQv+4BnA1c2Pf1FeCFfV+PB44FvjDCv4EkqZHZDn9dRvfX/q10V2L94agdV9U24Cy6E+m3A1dV1cYkZyTZfuXY1XR7EpuAi4HXzdS2b7M6yRfpAuMe4NJ++/nAE+jOudwEXFpVG0YdryTp0Ut34dU0hcmt/SW69FdVfaaqVkzbYIGZmJioycnJ2StKkh6WZH1VTQwrm21P5eGbRvZ7D5IkTWu2E/VHJdn+XY8Aj+vXA1RVPXGso5MkLSiz3fp+z7kaiCRp4Rv1eyqSJM3KUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc0YKpKkZgwVSVIzhookqRlDRZLUjKEiSWrGUJEkNWOoSJKaMVQkSc2MNVSSnJDkjiSbkqwZUp4k5/XlG5KsmK1tkqOS3Jjk1iQfSfLEKWXP7ss29uX7jHN+kqR/b2yhkmRP4HxgJbAcWJ1k+UC1lcCy/ud04IIR2r4bWFNVzwI+CLyxb7MX8D7gjKp6JnAc8G/jmp8k6ZHGuadyDLCpqjZX1YPAlcCqgTqrgMursw5YlOTAWdoeCVzfL18LvKJffhGwoao+B1BV/1xVD41rcpKkRxpnqBwE3D1lfUu/bZQ6M7W9DXhZv/wq4JB++QigklyT5OYk/33YoJKcnmQyyeTWrVt3cEqSpJmMM1QyZFuNWGemtq8BzkyyHtgXeLDfvhfwE8Cr+9eXJ/mpR3RSdVFVTVTVxJIlS2afhSRpZHuNse8tfH8vAuBg4J4R6+w9Xduq+gLdoS6SHAG8eEpff1tV9/dlVwMrgOsazEWSNIJx7qncBCxLcliSvYGTgbUDddYCp/RXgR0LPFBV987UNskB/esewNnAhX1f1wDPTvJD/Un75wOfH+P8JEkDxranUlXbkpxF92G/J3BJVW1MckZffiFwNXAisAn4FnDaTG37rlcnObNf/gBwad/mX5L8EV0gFXB1VX10XPOTJD1SqgZPc+w+JiYmanJycr6HIUkLSpL1VTUxrMxv1EuSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFUlSM4aKJKkZQ0WS1IyhIklqZqyhkuSEJHck2ZRkzZDyJDmvL9+QZMVsbZMcleTGJLcm+UiSJw70eWiSbyR5wzjnJkl6pLGFSpI9gfOBlcByYHWS5QPVVgLL+p/TgQtGaPtuYE1VPQv4IPDGgT7fAXys+YQkSbMa557KMcCmqtpcVQ8CVwKrBuqsAi6vzjpgUZIDZ2l7JHB9v3wt8IrtnSU5CdgMbBzTnCRJMxhnqBwE3D1lfUu/bZQ6M7W9DXhZv/wq4BCAJI8H3gScM9OgkpyeZDLJ5NatW0eejCRpduMMlQzZViPWmanta4Azk6wH9gUe7LefA7yjqr4x06Cq6qKqmqiqiSVLlsxUVZK0g/YaY99b6PciegcD94xYZ+/p2lbVF4AXASQ5AnhxX+e5wCuT/D6wCPheku9U1btaTEaSNLtxhspNwLIkhwH/CJwM/PxAnbXAWUmupAuFB6rq3iRbp2ub5ICqui/JHsDZwIUAVfW87Z0meSvwDQNFkubW2EKlqrYlOQu4BtgTuKSqNiY5oy+/ELgaOBHYBHwLOG2mtn3Xq5Oc2S9/ALh0XHOQJO2YVA2e5th9TExM1OTk5HwPQ5IWlCTrq2piWJnfqJckNWOoSJKaMVQkSc3s1udU+qvMvvwoulgM3N9oOAvB7jZfcM67C+e8Y55aVUO/6Ldbh8qjlWRyupNVu6Ldbb7gnHcXzrkdD39JkpoxVCRJzRgqj85F8z2AOba7zRec8+7COTfiORVJUjPuqUiSmjFUJEnNGCqzSHJCkjuSbEqyZkh5kpzXl29IsmI+xtnSCHN+dT/XDUn+PslR8zHOlmab85R6Ryd5KMkr53J84zDKnJMcl+SWJBuT/O1cj7G1Ef5vPynJR5J8rp/zafMxzlaSXJLkviS3TVPe/vOrqvyZ5ofuDslfAg6ne8bL54DlA3VOBD5G92CxY4FPz/e452DO/wnYr19euTvMeUq9j9PdXfuV8z3uOfg9LwI+Dxzarx8w3+Oegzn/JvD2fnkJ8DVg7/ke+6OY808CK4Dbpilv/vnlnsrMjgE2VdXmqnoQuBJYNVBnFXB5ddYBi5IcONcDbWjWOVfV31fVv/Sr6+georaQjfJ7Bng98FfAfXM5uDEZZc4/D3ygqr4CUFULfd6jzLmAfZMEeAJdqGyb22G2U1XX081hOs0/vwyVmR0E3D1lfUu/bUfrLCQ7Op/X0v2ls5DNOuckBwEvp38o3C5glN/zEcB+ST6ZZH2SU+ZsdOMxypzfBTyD7kmztwK/WlXfm5vhzYvmn1/jfPLjriBDtg1egz1KnYVk5PkkeQFdqPzEWEc0fqPM+Y+BN1XVQ90fsQveKHPeC/hx4KeAxwE3JllXVV8c9+DGZJQ5/wxwC/BC4EeAa5N8qqq+PuaxzZfmn1+Gysy2AIdMWT+Y7i+YHa2zkIw0nyTPBt4NrKyqf56jsY3LKHOeAK7sA2UxcGKSbVX1oTkZYXuj/t++v6q+CXwzyfXAUcBCDZVR5nwacG51Jxw2JfkH4OnAZ+ZmiHOu+eeXh79mdhOwLMlhSfYGTgbWDtRZC5zSX0VxLPBAVd071wNtaNY5JzmU7lHO/3UB/9U61axzrqrDqmppVS0F/hJ43QIOFBjt//aHgecl2SvJDwHPBW6f43G2NMqcv0K3Z0aSpwBHApvndJRzq/nnl3sqM6iqbUnOAq6hu3LkkqramOSMvvxCuiuBTgQ2Ad+i+0tnwRpxzr8NPBn40/4v9221gO/wOuKcdymjzLmqbk/y18AG4HvAu6tq6KWpC8GIv+e3Ae9NcivdoaE3VdWCvSV+kiuA44DFSbYAbwEeA+P7/PI2LZKkZjz8JUlqxlCRJDVjqEiSmjFUJEnNGCqSpGYMFe2y+rsJ35LktiR/0X/XokW/VydZlGTpDHd//WSSsV1mneS983Gn5CSvSnJ7kk/M9XtrYTBUtCv7dlX9WFX9KPAgcEaLTqvqxKr6fy36WoBeS/fFzxfM90D0g8lQ0e7iU8DTkuyf5EP9syPW9bebIcnz+72aW5J8Nsm+SQ5Mcv2UvZ3n9XXvSrK473evJJf1/f3lsL2hJC9KcmOSm/s9picMlD8jyWemrC9NsqFf/u0kN/Xvf1GG3Hhs6niSTCT5ZL/8+HTP07ipn9Oqfvszk3ymn9eGJMuG9Lk6ya39+759+1jo7vN2YZI/GKj/8iR/038z+8AkX0zyH0b83WgXYqhol5dkL7rnvtwKnAN8tqqeTffsjMv7am8AzqyqHwOeB3yb7tbv1/TbjqK70eCgI4GL+v6+Drxu4L0XA2cDP11VK4BJ4Nen1qmq24G9kxzeb/o54Kp++V1VdXS/t/U44CU7MPX/AXy8qo4GXgD8QZLH0+2xvbOf1wTd/Z+mjvmHgbfT3VTxx4Cjk5xUVb/Tj//VVfXGgTl8EPgn4EzgYuAtVfVPOzBW7SIMFe3KHpfkFroPwq8A76H7S/vPAKrq48CTkzwJuAH4oyS/Aiyqqm1094o6LclbgWdV1b8OeY+7q+qGfvl9PPKOzccCy4Eb+rGcCjx1SD9XAT/bL/8c8P5++QVJPt3fNuSFwDNHnz4vAtb07/tJYB/gUOBG4DeTvAl4alV9e6Dd0cAnq2pr/+/wf+ge9jSb1wNvBr5bVVfswDi1C/HeX9qVfbv/a/xhww4fAVVV5yb5KN19kNYl+emquj7JTwIvBv4syR9U1eWDbWdZD3BtVa2eZazvB/4iyQf68dyZZB/gT4GJqrq7D7d9hrTdxvf/QJxaHuAVVXXHQP3bk3y6n9c1SX6pD9ip7XbGQXT3CHtKkj128eeQaBruqWh3cz3wauiev053a/evJ/mRqrq1qt5Ot2fz9CRPBe6rqovp9nKGPb/70CT/sV9eDfzdQPk64D8neVr/nj+U5IjBTqrqS8BDwG/x/b2U7QFxf38eZrqrve6ie+4JwCumbL8GeP32IE3ynP71cGBzVZ1Hd5faZw/092ng+UkWJ9mzn9eMz6fvDzFeSnfI8HYGDvFp92GoaHfzVmCiPxF+Lt3hKIBf609Kf47ufMrH6O7uekuSz9J9WL9zSH+3A6f2/e0PXDC1sKq2Ar8IXNHXWUf3fI5h3g/8Av35lP4Ks4vpzgV9iO5w3DDnAO9M8im6YNrubXR3pN3QX/r8tn77zwG39YfFns73zyttH/O9dIexPkH3HPebq+rD07z3dr8JfKqqPkUXKL+U5BmztNEuyLsUS5KacU9FktSMoSJJasZQkSQ1Y6hIkpoxVCRJzRgqkqRmDBVJUjP/H/qma1Z28egKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "prior.plot()\n", "\n", "plt.xlabel('Possible values of x')\n", "plt.ylabel('Probability')\n", "plt.title('Prior');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the likelihoods for heads and tails:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "likelihood_heads = xs\n", "likelihood_tails = 1 - xs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we toss the coin twice and get one heads and one tails.\n", "\n", "We can compute the posterior probability for each value of $x$ like this:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "posterior = prior * likelihood_heads * likelihood_tails\n", "posterior /= posterior.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's what the posterior distribution looks like." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5VUlEQVR4nO3deXhU1fnA8e+bjRAIhJAAAQIJEDaRNSwiiLvghrsgCioWqVvVLmptbe3qWit14Qd1ww1Rq6JiERUElC0IhB1CCCQQSNgCJIRs7++PubRpDMkEMrkzyft5nnlm5t5z7rwnhHlz7j33HFFVjDHGGG8FuR2AMcaYwGKJwxhjTI1Y4jDGGFMjljiMMcbUiCUOY4wxNWKJwxhjTI1Y4jAmQIjIOBH50u04jBG7j8OYyolIBtAaKAXygTnAvap69BSO9Xugi6reXJsxGuMG63EYU7UrVLUp0B8YCPzGjSBEJOQ06oqI2P91U2vsl8kYL6jqLuALoJeIXCki60XkkIgsEJEeJ8qJyEMisktEjojIZhG5QERGAr8GbhSRoyKyxinbXEReEZFsp86fRCTY2XeriHwnIs+JyAHg9862xeU+a6iIrBCRPOd5aLl9C0TkzyLyHVAAdKqTH5RpECxxGOMFEYkHLgWOAO8C9wOxeE5ffSoiYSLSDbgHGKiqkcAlQIaq/hv4C/CeqjZV1T7OYd8ASoAuQD/gYuCOch87GEgHWgF/rhBPNPA5MAVoCfwN+FxEWpYrdgswCYgEdtTCj8EYwBKHMdX5WEQOAYuBb4ENwOeqOk9Vi4FngMbAUDzXQhoBPUUkVFUzVHVbZQcVkdbAKOB+Vc1X1RzgOWBMuWK7VfUfqlqiqscqHOIyYKuqvunsfxfYBFxRrszrqrre2V98mj8HY/7DEocxVbtKVaNUtaOq3gW0pdxf76paBmQC7VQ1DU9P5PdAjojMFJG2JzluRyAUyHZOeR0C/g9P7+KEzCri+p84HDuAdl7WN+aUWeIwpmZ24/nSBzwXnoF4YBeAqr6jqsOcMgo86RStOHwxEzgOxDiJKUpVm6nqGeXKVDXk8X/icHQ4EYcX9Y05ZZY4jKmZWcBlzkXvUODneBLA9yLSTUTOF5FGQCFwDM/pK4C9QMKJ0U2qmg18CTwrIs1EJEhEOovICC/jmAN0FZGbRCRERG4EegKf1VpLjTkJSxzG1ICqbgZuBv4B7MNzTeEKVS3Cc33jCWf7HjynnX7tVH3fed4vIj84r8cDYXiumxwEPgDivIxjP3A5nsS1H/gVcLmq7jud9hnjDbsB0BhjTI1Yj8MYY0yNWOIwxhhTI5Y4jDHG1IglDmOMMTVyyhOnBZKYmBhNSEhwOwxjjAkoK1eu3KeqsRW3N4jEkZCQQEpKitthGGNMQBGRSuc4s1NVxhhjasQShzHGmBqxxGGMMaZGLHEYY4ypEUscxhhjasQShzHGmBqxxGGMMaZGGsR9HMa4paxMyT16nN2HjnGwoIiD+cXkHSumqLSMopIySsqUsGAhLCSI8NBgoiLCaBERSssmjWjXojHNG4e63QRjfsQShzG1ZE9eIaszD7Fpz2E2Zh8mLecoWQePcbyk7JSP2bxxKB2iI+jaOpIecZH0jGtG7/gomjay/7rGPfbbZ8wp2nu4kG8357IkfT8rMg6QdfAYACKQ2LIJXVtHckGP1sS3aEzbqMa0bNqI6IgwmjUOoVFIMGEhQQQJlJQpRSVlHCsu5VBBMQcLisg9cpysgwVkHjhGxv58Fm7N5cMfsgAIEujZthnJHaMZnhTDWZ1bEhFm/5VN3bHfNmNqYMveI3yWms1XG/ayIfswADFNGzEosQW3n51Ivw5RdG/TjMZhwV4fMzRYCA0OokmjEGKaNjppuX1Hj7N+92FW7jhISsYB3luRyevfZxAWEsTgxGhG9mrDyDPa0LKKYxhTGxrECoDJyclqc1WZU7X3cCEfrMzik9W72LL3KEECyQnRnNetFed2i6V7m0hEpM7jOl5SyortB1mwOYdvNuWQvi+f4CBhaOeWXNO/HaN6xREe6n0CM6YiEVmpqsk/2m6Jw5gfKytT5m/O4Z1lO5m/OYcyhYEJLbiiT1tG9mpDq8hwt0P8H6rKpj1H+Cx1N5+uyWbngQIiw0O4qm87bjmrI11bR7odoglAljgscRgvFBSV8OHKLF79LoPt+/KJjWzE9QPac0NyPAkxTdwOzytlZcrS7fuZtSKTOev2UFRSxvCkGCYOS2RE11hXekcmMFnisMRhqnCksJgZS3bwz0XpHCwopk98FBOHJTKqVxtCgwP3dqeD+UW8s3wnb3yfQc6R4/SIa8Z953fhkjPaEBRkCcRUzRKHJQ5TiYKiEl77LoNpC9PJO1bMed1iufu8Lgzo2KJe/WVeVFLG7DW7eWl+Gun78unWOpIHLurKJWe0rlftNLXLEoclDlNOSWkZH6zM4m/ztpBz5DgXdG/FfRck0Sc+yu3QfKq0TPksdTdTvt7Kttx8BnRswa8v7c6AjtFuh2b8kCUOSxzGsWTbfn43ex1b9h6lf4cofn1pD5ITGtYXZ0lpGe+vzOI5J3FedmYcj17Wg7ZRjd0OzfiRkyUOn568FZGRIrJZRNJE5OFK9ouITHH2p4pI/3L7XhWRHBFZd5Jj/0JEVERifNkGU3/sySvk3ndXMXb6UgqKSnl5XH8+/OnQBpc0AEKCgxg7qAMLfnkuD1zYla827uWCZ7/lpQVpHC8pdTs84+d8ljhEJBh4ERgF9ATGikjPCsVGAUnOYxLwcrl9rwMjT3LseOAiYGftRm3qo7Iy5c2lO7jwb98yd/0efnZBEl89OIJRZ8Y1+PP7EWEh/OxCz89jeFIMT/17M5dNWczKHQfdDs34MV/2OAYBaaqarqpFwExgdIUyo4EZ6rEUiBKROABVXQgcOMmxnwN+BdT/82zmtGzfl8+Y6Uv57cfr6BsfxbwHzuGBi7rajXEVxEdHMG18Mq/dOpCC4yVcN/V7Hv90PQVFJW6HZvyQLxNHOyCz3PssZ1tNy/wPEbkS2KWqa6opN0lEUkQkJTc31/uoTb2gqsxYksGo5xeyMfswT13bmzcnDqJjy8C4F8Mt53VvxZcPjuCWIR157bsMLn1+ET/stN6H+V++TByVnQOo2EPwpsx/C4tEAI8Cj1X34ao6TVWTVTU5Nja2uuKmHsk5Ushtr6/gsU/WMzixJV89OIIbBsY3+NNS3mraKIQ/jO7FzElDKC5Vrp+6hOfmbaG49NRn+TX1iy8TRxYQX+59e2D3KZQprzOQCKwRkQyn/A8i0ua0ozX1wrdbchn590Us2bafP4w+g9dvG0jrZv41PUigGNKpJV/cP5zRfdry/NdbueH/lrDr0DG3wzJ+wJeJYwWQJCKJIhIGjAFmVygzGxjvjK4aAuSpavbJDqiqa1W1laomqGoCnsTTX1X3+KgNJkCUlJbx9NxNTHh1ObFNG/H5fcMYf1aC9TJOU7PwUP52Y1/+MbYfW/ce5bIpi/hm0163wzIu81niUNUS4B5gLrARmKWq60VksohMdorNAdKBNGA6cNeJ+iLyLrAE6CYiWSIy0VexmsC27+hxxv1zGS/O38aNyfF8fPfZdGllk/rVpiv6tOXTe4fRtnljbn89hSf/vYnSMhub0lDZDYAmoK3NymPSmykcLCjiz1edybUD2rsdUr1WWFzK459u4N3lOxnRNZYpY/rRPMKWt62vXLkB0Bhf+mhVFtdN/Z4gET6YPNSSRh0IDw3mr9ecyV+uPpPvt+1j9IuL2bL3iNthmTpmicMEnLIy5em5m3jgvTX06xDF7HvOple75m6H1aDcNLgDMycNIb+olGtf+p6FW2zIe0NiicMElMLiUu6buYoX529jzMB43pw42JZKdcmAjtHMvuds2kdHcNvrK3h72Q63QzJ1xBKHCRgH8ou4afpSPkvN5pFR3fnrNWcG9FoZ9UFc88a8P/kszkmK4dGP1vGXORsps4vm9Z79rzMBIetgAddN/Z51uw/z0rj+3Dmisw219RNNG4UwfXwytwzpyLSF6fz8/TV2s2A9F+J2AMZUZ8veI4x/ZTn5RSW8efsgBndq6XZIpoKQ4CD+MPoMWjdrxDNfbuFAfhEv39yfiDD7iqmPrMdh/NqqnQe5fuoSylSZdedZljT8mIhwz/lJ/PWaM1m0NZebpi8jr6DY7bCMD1jiMH5rWfp+bv7nMpo3DuXDnw6lR1wzt0MyXhg7qAMvjRvAht2HGTt9KfuPHnc7JFPLLHEYv7Roay4TXltOm+bhzLrzLOKjI9wOydTAyF5tmDZ+ANtyjzJm2lJyDhe6HZKpRZY4jN+ZvzmHia+nkNCyCe/deRZtmtskhYHo3G6teO22gew6dIwbpy1lryWPesMSh/ErC7fkcuebK0lq3ZSZk4YQY/doBLShnWN4c+Igcg4XMnb6UnKOWPKoDyxxGL+xeOs+fjIjhc6xTXlr4mCiIsLcDsnUggEdo3n99kHsySvkpunLyD1i1zwCnSUO4xeWpe/njhkrSIxpwtt3DKZFE0sa9cnAhGheu3Uguw4eY9w/l3Iwv8jtkMxpsMRhXJeadYiJb6TQLqoxb90xmGhLGvXS4E4teeXWZDL2FzDhteUcKbShuoHKEodx1da9R5jw6nKiIkJ5+w67plHfDe0cw8vj+rNh92HueCOFwuJSt0Myp8ASh3FN5oECxv1zGSHBQbw1cbCNnmogLujRmmdv6MPyjAPc9fYPNj1JALLEYVxxIL+I8a8up7C4lDcnDiIhponbIZk6NLpvO/44uhffbMrhkX+tpSEsKFef2EQyps4VFJVw++sr2H3oGG/dMZjubeyO8Ibo5iEdyTlynClfb6VNs3B+cUk3t0MyXvJpj0NERorIZhFJE5GHK9kvIjLF2Z8qIv3L7XtVRHJEZF2FOk+LyCan/EciEuXLNpjaVVJaxr3vrCI16xDPj+nHwIRot0MyLnrgwiTGDIznhflpvLnU1vMIFD5LHCISDLwIjAJ6AmNFpGeFYqOAJOcxCXi53L7XgZGVHHoe0EtVewNbgEdqN3LjK6rK72av5+tNOfxhdC9G9mrjdkjGZSLCn67qxYU9WvHYJ+v4asNet0MyXvBlj2MQkKaq6apaBMwERlcoMxqYoR5LgSgRiQNQ1YXAgYoHVdUvVbXEebsUsIWmA8Qri7fz9rKdTB7RmZuHdHQ7HOMnQoKDmDK2H73aNue+matYtyvP7ZBMNXyZONoBmeXeZznbalqmKrcDX1S2Q0QmiUiKiKTk5tp6yG6bu34Pf56zkUvPbMOv7Fy2qSAiLIRXJiQT1TiUiW+sYE+eTU3iz3yZOCpbnq3i0AlvylR+cJFHgRLg7cr2q+o0VU1W1eTY2FhvDml8ZN2uPO6fuZre7aP42w19CQqylfvMj7VqFs4rtw4k/3gpt7++goKikuorGVf4MnFkAfHl3rcHdp9CmR8RkQnA5cA4tXF8fi3nSCE/mZFCdJMw/jk+mfDQYLdDMn6sR1wz/nFTPzbtOczPZ62x9cv9lC8TxwogSUQSRSQMGAPMrlBmNjDeGV01BMhT1eyqDioiI4GHgCtVtcAXgZvacbyklMlvruRQQTHTxg8gNtLuCjfVO69bKx4Z1YMv1u3hH9+kuR2OqYTPEodzAfseYC6wEZilqutFZLKITHaKzQHSgTRgOnDXifoi8i6wBOgmIlkiMtHZ9QIQCcwTkdUiMtVXbTCnTlX5zUfr+GHnIZ65vg9ntG3udkgmgNwxPJFr+rfjua+28O91Vf4taVzg0xsAVXUOnuRQftvUcq8VuPskdceeZHuX2ozR+MYb32fw/sos7ju/C5f1jnM7HBNgRIS/XH0m6bn5PDhrDYkxTenWJtLtsIzDphwxtW5FxgH+9PlGLuzRmvsv7Op2OCZAhYcG83+3DKBJoxAmv7WSwzabrt+wxGFqVc7hQu56+wfioyP42419bASVOS2tm4Xz0rj+ZB4o4MH37GK5v7DEYWpNUUkZd739A0cLS5h68wCahYe6HZKpBwYmRPPoZT34auNeXlpgF8v9gSUOU2v++sVGUnYc5Knretv5aFOrbh2awOi+bXl23hYWbbUbet1micPUin+vy+a17zK4dWgCV/Rp63Y4pp4REf56zZkktWrK/TNXs/ew3VnuJksc5rTt2J/PL99PpU98FL++tIfb4Zh6KiIshJfG9aegqJR731lFiS0A5RpLHOa0FBaXcvc7PyACL4ztR1iI/UoZ3+nSKpK/XNOL5RkH+Nu8LW6H02DZ/3JzWv46ZyPrdh3m2Rv6Eh8d4XY4pgG4ul97xgyM56UF2/h2i13vcIMlDnPK5m3YyxtLdnD72Ylc1LO12+GYBuT3V55Bt9aR/HzWanKPHHc7nAbHEoc5JXvyCvnVB2s4o20zHhpl06SbuhUeGsyUsf04UljCz9+3+zvqmiUOU2OlZcoD763meEkZ/xjbj0YhNuOtqXvd2kTy28t7snBLLq8s3u52OA2KJQ5TY1O/3caS9P08fuUZdIpt6nY4pgEbN7gDl5zRmqfmbmJtlq0cWFcscZgaWZuVx3PztnB57ziuG2Cr9hp3iQhPXtublk0acf97qzhWVOp2SA2CJQ7jtWNFpdz/3ipimjbiz1ediYjNQ2XcFxURxjPX92Fbbj5PfLHR7XAaBEscxmtPfLGRbbn5PHtDH5pH2DxUxn8MS4rhtrMTeGPJDhZsznE7nHrPEofxyoLNOf8Zent2lxi3wzHmRx4a2Z2kVk355QepHMgvcjuces0Sh6lWXkExD32YStfWTfnVSBt6a/xTeGgwfx/Tl0MFRfz2k3Vuh1Ov+TRxiMhIEdksImki8nAl+0VEpjj7U0Wkf7l9r4pIjoisq1AnWkTmichW57mFL9tg4Pefrmff0SKevb4v4aE29Nb4rzPaNue+85P4PDWbz1J3ux1OveWzxCEiwcCLwCigJzBWRHpWKDYKSHIek4CXy+17HRhZyaEfBr5W1STga+e98ZG56/fw0apd3H1eF85sb+uGG//303M707t9c3778Tq7q9xHfNnjGASkqWq6qhYBM4HRFcqMBmaox1IgSkTiAFR1IXCgkuOOBt5wXr8BXOWL4A0cyC/i0Y/W0jOuGfecZ0u9m8AQEhzEs9f3Ib+olEc/Wouq3VVe23yZONoBmeXeZznbalqmotaqmg3gPLeqrJCITBKRFBFJyc21idBOxWOfrCPvWDHP3tDHZr01ASWpdSQ/v6grX27Yyyer7ZRVbfPlt0Flg/wrpn5vypwSVZ2mqsmqmhwbG1sbh2xQvly/h89Ss7n3/CR6xDVzOxxjauyO4Z3o1yGKxz9dz76jdsqqNvkycWQB8eXetwcqpn5vylS098TpLOfZBm3Xsrxjxfzm43X0iGvGT8/t7HY4xpyS4CDhqWt7k3+8lN/NXu92OPWKLxPHCiBJRBJFJAwYA8yuUGY2MN4ZXTUEyDtxGqoKs4EJzusJwCe1GbSBP3++gf35RTx9XW9Cg+0UlQlcSa0juff8Lnyems3c9XvcDqfe8Nm3gqqWAPcAc4GNwCxVXS8ik0VkslNsDpAOpAHTgbtO1BeRd4ElQDcRyRKRic6uJ4CLRGQrcJHz3tSSRVtzmZWSxaRzOtGrnY2iMoFv8rmd6RHXjN98vI68gmK3w6kXpCGMOEhOTtaUlBS3w/B7BUUlXPzcQsKCg5jzs+F2z4apN9btymP0i99x/YD2PHFtb7fDCRgislJVkytut/MQ5j+e/2orWQeP8ddrzrSkYeqVXu2aM3FYIjNXZLIsfb/b4QQ8SxwGgPW78/jn4u2MGRjP4E4t3Q7HmFp3/4VJtG/RmEc+WsvxEpt+/XRY4jCUlimP/GstLSLCeGRUD7fDMcYnIsJC+NNVvUjPzeflBdvcDiegeZU4RORyEbEkU0+98X0GqVl5PHZFT5su3dRr53ZrxZV92vLS/G2k5RxxO5yA5W0yGANsFZGnRMT+JK1H9uQV8uyXmxnRNZYrese5HY4xPvfby3vSOCyY33y8zqYjOUVeJQ5VvRnoB2wDXhORJc6UHpE+jc743B8/30BJmfLH0b1sRT/TIMRGNuJXI7uxNP0AH6/e5XY4Acnr00+qehj4EM9khXHA1cAPInKvj2IzPrZwSy6fp2Zz93ld6NAywu1wjKkzYwd2oE98FH/+fCN5x+zejpry9hrHlSLyEfANEAoMUtVRQB/gFz6Mz/hIYXEpj32yjsSYJtw5opPb4RhTp4KChD9f1YsD+UU8M3ez2+EEHG97HNcBz6lqb1V9WlVzAFS1ALjdZ9EZn/m/b9PJ2F/AH0afQaMQu2fDNDy92jVn/FkJvLVsB6lZh9wOJ6B4mziynfUx/kNEngRQ1a9rPSrjU5kHCnhpQRqX9Y5jeJLNHGwargcv7kpM00b89pP1lJXZhXJveZs4Lqpk26jaDMTUnT9+toEgEX5zmQ2QMw1bs/BQHhnVnTWZh/hgZZbb4QSMKhOHiPxURNYC3Z01wU88tgOpdROiqU3fbsnlyw17ufeCLsQ1b+x2OMa47up+7RjQsQVP/nuTXSj3UnU9jneAK/BMXX5FuccAZ4iuCSBFJWU8Pns9iTFNmDgs0e1wjPELIsLjV57BgYIinpu3xe1wAkJ1iUNVNQO4GzhS7oGIRPs2NFPbXv1uO+n78nnsip52QdyYcnq1a85Ngzrw5tIdbNpz2O1w/J43PQ6AlUCK87yy3HsTIHIOFzLl661c2KMV53WrdJl2Yxq0X1zcjcjwEB6fvcHuKK9GlYlDVS93nhNVtZPzfOJhg/8DyFNzN1NcWsZvLuvpdijG+KUWTcJ48KKuLEnfz9z1e90Ox6+FVLVTRPpXtV9Vf6jdcIwvpGZ5RozceU4nEmKauB2OMX7rpkEdeGvpDv4yZyPndY+1U7onUd2pqmereDxT3cFFZKSIbBaRNBF5uJL9IiJTnP2p5RPVyeqKSF8RWSoiq0UkRUQGedfUhklV+cOnG4hpGsY953dxOxxj/FpIcBC/vbwnOw8U8OriDLfD8VtV9jhU9bxTPbCIBAMv4rkHJAtYISKzVXVDuWKjgCTnMRh4GRhcTd2ngMdV9QsRudR5f+6pxlnffZqaTcqOgzxxzZlEhtuU6cZUZ3hSLBf2aMUL32zl2gHtaBUZ7nZIfqe6+zjOd56vqexRzbEHAWmqmq6qRXgmRxxdocxoYIZ6LAWiRCSumroKNHNeNwd2e9nWBqewuJQnv9hEz7hmXJ8c73Y4xgSMRy/rSVFpGc/OteG5lamyxwGMwDOx4RWV7FPgX1XUbQdklnufhadXUV2ZdtXUvR+YKyLP4El8Q6tsQQP2yuLt7Dp0jGeu70NwkE2Zboy3EmOaMP6sBF79bju3np1Aj7hm1VdqQKo7VfU75/m2Uzh2Zd9UFce4naxMVXV/Cjygqh+KyA3AK8CFP/pwkUnAJIAOHTp4G3O9se/ocV5esI0Le7TmrM62hrgxNXXv+V34YGUWf5mzkRm3D7L1asrxdlr1ls5F7B9EZKWIPC8i1X0bZQHlz4+058enlU5Wpqq6E/hvT+d9PKe1fkRVp6lqsqomx8Y2vIn8npu3hcLiUh65tLvboRgTkKIiwrjvgiQWbd3Hgi25bofjV7yd5HAmkAtci2eK9VzgvWrqrACSRCRRRMLwLD87u0KZ2cB4Z3TVECBPVbOrqbsbzyk0gPOBrV62ocHYuvcIM1dkMm5wBzrHNnU7HGMC1i1DOtKxZQR/+XwjJaVlbofjN7xNHNGq+kdV3e48/gREVVVBVUuAe4C5wEZglqquF5HJIjLZKTYHSAfSgOnAXVXVder8BHhWRNYAf8E5HWX+669fbCIiNJj7LkhyOxRjAlpYSBAPj+zO1pyjzEqx2XNPqO7i+AnzRWQMMMt5fx3weXWVVHUOnuRQftvUcq8VzzxYXtV1ti8GBngZd4OzZNt+vtmUw0Mju9OyaSO3wzEm4I3s1Ybkji147qstXNWvLRFh3n5t1l/VDcc9IiKHgTvxzFtV5DxmAg/4PjxTE6rKE19sJK55OLedneB2OMbUCyLCI5d2J/fIcf65aLvb4fiF6uaqilTVZs5zkKqGOI8gVbXxaX7m87XZrMnK48GLuhIealMlGFNbBnSM5pIzWvN/325j39HjbofjOm+vcSAiLURkkIicc+Lhy8BMzRSXlvH03M10ax3JNf3bux2OMfXOLy/pTmFJGS98k+Z2KK7zdjjuHcBCPBerH3eef++7sExNvbt8Jzv2F/DQqG52s58xPtClVVNuSI7n7WU72LE/3+1wXOVtj+NnwEBghzN/VT88Q3KNH8g/XsKUr7cyODHa1towxoceuDCJkKAgnvmyYU9F4m3iKFTVQgARaaSqm4BuvgvL1MSri7ez72gRD43qbne3GuNDrZqFc/uwBD5ds5v1u/PcDsc13iaOLBGJAj4G5onIJ9jkgn7hYH4R0xamc1HP1vTv0MLtcIyp9yad05nmjUN5Zu5mt0NxjVeJQ1WvVtVDqvp74Ld45oe6yodxGS9N/XYbR4tK+MXF1gE0pi40bxzK5BGdmb85lxUZB9wOxxU1GVXVX0TuA3oDWc5058ZFe/IKef37DK7u245ubSLdDseYBuPWoQm0imzEU//e1CDXJ/d2VNVjwBtASyAGeE1EfuPLwEz1pnyzlTJVHrioq9uhGNOgNA7zTOmzIuMgCzY3vHFC3vY4xgIDVfV3zlTrQ4BxvgvLVGfH/nxmrchkzMAOxEdHuB2OMQ3OjQPj6dgygqfnbqasrGH1OrxNHBlA+fUTGwHbaj0a47Xnv95KcJBwr60jbowrQoOD+NkFSWzIPszc9XvcDqdOVTdX1T9EZApwHFgvIq+LyGvAOuBoXQRofiwt5ygfr9rF+LM60qqZrYdsjFtG921H59gmPPfVFkobUK+jumkeU5znlcBH5bYv8Ek0xit//2oL4aHBTB7R2e1QjGnQgoOEBy7qyj3vrOKz1N2M7tvO7ZDqRHVLx75x4rWzoNKJq7CbVbXYl4GZym3ac5jPUrO5+7zONm26MX7g0l5xdG+Txt+/2splZ8YREuz1YNWA5e2oqnPxrLT3IvASsMUmOXTHc/O2EBkewqTh1tswxh8EBQkPXtSV7fvy+WjVLrfDqRPepsZngYtVdYSqngNcAjznu7BMZdbtymPu+r3cMawTzSNC3Q7HGOO4qGdrerdvzpRvtlLcAJaY9TZxhKrqf+6vV9UtgH1z1bHnv95Ks/AQbhuW4HYoxphyRIT7L0wi88AxPvqh/vc6vE0cK0XkFRE513lMx3PBvEoiMlJENotImog8XMl+EZEpzv5UEenvTV0RudfZt15EnvKyDQFt3a485m3Yyx3DO9Es3HK2Mf7mvG6t6N2+Of+YX/97Hd4mjsnAeuA+PFOsb3C2nZSIBOO5JjIK6AmMFZGeFYqNApKcxyTg5erqish5wGigt6qeATzjZRsC2onexq22JKwxfqkh9TqqXXVdRIKAlaraC/hbDY49CEhT1XTnODPxfOFvKFdmNDBDPZO9LBWRKBGJAxKqqPtT4AlVPQ6gqjk1iCkgnehtPHhRV+ttGOPHyvc6ru7fjtB6OsKq2lapahmwRkQ61PDY7YDMcu+znG3elKmqbldguIgsE5FvRWRgZR8uIpNEJEVEUnJzA3suGettGBMYGkqvw9t0GIfnzvGvRWT2iUc1dSpbUajirZUnK1NV3RCgBZ75sn4JzJJKVi9S1WmqmqyqybGxsdWE6r/W7/b0NiYOs2sbxgSC8r2Oknp6raPaU1WOx0/h2FlAfLn37fnx4k8nKxNWRd0s4F/O6a3lIlKGZ8bewO5WnMQL36QR2ch6G8YEChHh3vOT+MmMFD5ZvZtrB7R3O6RaV91cVeEicj9wPdAd+E5Vvz3xqObYK4AkEUl07jofA1TspcwGxjujq4YAeaqaXU3dj4Hznfi64kky+7xrbmDZsvcIX6zbw61nJ9C8sfU2jAkUF/ZoRY+4Zrw4P61ezmFV3amqN4BkYC2eEU7PentgVS0B7gHmAhuBWaq6XkQmi8iJEVlzgHQgDZgO3FVVXafOq0AnEVkHzAQmaD1dSeXF+Wk0CQvm9rMT3Q7FGFMDnl5HF9L35TNnbbbb4dQ6qeo7V0TWquqZzusQYLmq9j9pBT+VnJysKSkp1Rf0I9v35XPBswv4yTmdeGRUD7fDMcbUUFmZcvHfFxIswhc/G05QUGWXbv2biKxU1eSK26vrcfxnIkOnF2DqyIvz0wgLCeKOYZ3cDsUYcwqCgoR7zuvC5r1H+HLDXrfDqVXVJY4+InLYeRwBep94LSKH6yLAhijzQAEfrdrF2EEdiI20GXCNCVSX944joWUEL8zfWq/WJq8ycahqsKo2cx6RqhpS7nWzugqyoZm2MJ0ggUnnWG/DmEAWEhzEXed2Yd2uwyzcWn/G8NTP2xoDWM6RQt5LyeTa/u2Ja97Y7XCMMafpqn7tiGsezovz09wOpdZY4vAzry7OoKS0jDttdT9j6oWwkCB+MrwTy7cfICXjgNvh1ApLHH4kr6CYt5bu4LLebUmMaeJ2OMaYWjJmUDzRTcJ4acE2t0OpFZY4/MiMJRkcPV7CXedab8OY+iQiLITbz07gm005bNgd+OOKLHH4iYKiEl77PoPzu3vuODXG1C+3nJVA00YhvLQg8K91WOLwE7NWZHIgv8h6G8bUU80bh3LzkI7MWZvNjv35bodzWixx+IHi0jKmL9rOwIQWJCdEux2OMcZHbj87gZCgIKYvSnc7lNNiicMPfJ6aza5Dx7jzHOttGFOftWoWzjX92/F+Shb7jh53O5xTZonDZarK1G+3kdSqKed3b+V2OMYYH/vJOZ0oKi3jje8z3A7llFnicNm3W3LZtOcIk87pFJCToBljaqZzbFMu7tmaGUt2kH88MKcAtMThsqnfbiOueTij+1ZcVdcYU19NHtGZvGPFvLt8p9uhnBJLHC5ak3mIpekHmDgskbAQ+6cwpqHo16EFgxOjeWXxdooDcHlZ+7Zy0bRF6USGh3DjwPjqCxtj6pU7R3QiO6+Qz1MDb6EnSxwuyTxQwBdrs7lpcAciw21ZWGMamnO7tqJLq6ZMW5gecFOu+zRxiMhIEdksImki8nAl+0VEpjj7U0Wkfw3q/kJEVERifNkGX3ll8XaCRLhtqC0La0xDFBQk/GR4IhuyD/P9tv1uh1MjPkscIhIMvIhnrfKewFgR6Vmh2CggyXlMAl72pq6IxAMXAQF5ZelQQRGzUjK5sm9b2jQPdzscY4xLRvdtR0zTRkxbGFg3BPqyxzEISFPVdFUtAmYCoyuUGQ3MUI+lQJSIxHlR9zngV0Bg9e8cby/bSUFRKT8Zbgs1GdOQhYcGc+vQjny7JZfNe464HY7XfJk42gGZ5d5nOdu8KXPSuiJyJbBLVddU9eEiMklEUkQkJTc399Ra4APHS0p5/fsMhifF2GSGxhjGDe5I49DggJqGxJeJo7K72Sr2EE5WptLtIhIBPAo8Vt2Hq+o0VU1W1eTY2Nhqg60rs1fvJvfIcVsW1hgDQIsmYdyQ3J5PVu8i53Ch2+F4xZeJIwsoP860PbDbyzIn294ZSATWiEiGs/0HEWlTq5H7iKryyuLtdG8TybAuAXlN3xjjA7cPS6SkTHlz6Q63Q/GKLxPHCiBJRBJFJAwYA8yuUGY2MN4ZXTUEyFPV7JPVVdW1qtpKVRNUNQFPgumvqnt82I5a8/22/Wzac4Tbz05ExKYXMcZ4dGzZhAt7tOatpTsoLC51O5xq+SxxqGoJcA8wF9gIzFLV9SIyWUQmO8XmAOlAGjAduKuqur6Kta68sng7MU3DuLJvW7dDMcb4mYnDEjlYUMy/ftjldijVCvHlwVV1Dp7kUH7b1HKvFbjb27qVlEk4/Sjrxrbco3yzKYf7L0wiPDTY7XCMMX5mcGI0vdo149XvtjNmYLxfT3pqd47XkVcXbycsJIibh3R0OxRjjB8SESYOSyQt5yjfbvWfkaCVscRRBw7mF/HhD1lc1bctMU0buR2OMcZPXXZmW1pFNuLVxdvdDqVKljjqwLsrdlJYXMbtw2x6EWPMyYWFBDFhaAKLtu5jy17/vSHQEoePFZeW8eaSHQzt3JLubeyGP2NM1cYO6kCjkCBe+y7D7VBOyhKHj81dv4fsvEJuO9t6G8aY6kU3CeOqvu34aFUWhwqK3A6nUpY4fOz17zLoEB1h64kbY7x227AECovLmLkis/rCLrDE4UNrs/JI2XGQCUMTCPbjoXXGGP/SvU0zzurUkhnfZ1DihysEWuLwode+206TsGCuT27vdijGmABz29kJ7M4r5MsNe90O5UcscfhIzpFCPk3dzXUD2tPMVvgzxtTQBT1aEx/dmNe+87+huZY4fOTdZZkUlyoThia4HYoxJgAFBwkTzkpgRcZB1u3Kczuc/2GJwweKS8t4e9kORnSNpVNsU7fDMcYEqOsHxNM4NJg3l/jXrLmWOHxg7vo95Bw5zoShNr2IMebUNY8I5ap+7fh49S4O5vvP0FxLHD4w4/sddIiOYERXG4JrjDk9E4Z25HhJGbNS/GdoriWOWrZh92GWZxxg/FkdbQiuMea0dW/TjMGJ0by5dAelZRUXUXWHJY5aNmNJBuGhQVw/IL76wsYY44UJQxPIOniM+Zty3A4FsMRRqw4VFPHx6l1c3a8dzSNsCK4xpnZc1LM1bZqF88aSDLdDASxx1KoPVmZRWFzGLUMS3A7FGFOPhAYHMW5wBxZt3ce23KNuh+PbxCEiI0Vks4ikicjDlewXEZni7E8Vkf7V1RWRp0Vkk1P+IxGJ8mUbvFVWpry1dAcDE1rQs63NgmuMqV03DoonNFh4e+lOt0PxXeIQkWDgRWAU0BMYKyI9KxQbBSQ5j0nAy17UnQf0UtXewBbgEV+1oSYWp+0jY3+BrfBnjPGJVpHhjOwVx/srMykoKnE1Fl/2OAYBaaqarqpFwExgdIUyo4EZ6rEUiBKRuKrqquqXqnrip7YU8IuJoGYs2UFM0zBG9mrjdijGmHrqliEdOVJYwqdrdrsahy8TRzug/MDjLGebN2W8qQtwO/DFaUd6mrIOFvDNpr3cODCeRiHBbodjjKmnBia0oFvrSGYs2YGqe0NzfZk4KruJoWJLT1am2roi8ihQArxd6YeLTBKRFBFJyc317cLv7y73nHMcO6iDTz/HGNOwiQg3n9WR9bsPszrzkGtx+DJxZAHlb2ZoD1TsX52sTJV1RWQCcDkwTk+SdlV1mqomq2pybGzsKTeiOsdLSnlvRSbnd29N+xYRPvscY4wBuLpfO5qEBfPmUvfmr/Jl4lgBJIlIooiEAWOA2RXKzAbGO6OrhgB5qppdVV0RGQk8BFypqgU+jN8r/163h31Hi7jlLLsobozxvaaNQrimf3s+S83mgEvzV/kscTgXsO8B5gIbgVmqul5EJovIZKfYHCAdSAOmA3dVVdep8wIQCcwTkdUiMtVXbfDG28t20iE6guFdYtwMwxjTgNw8pCNFJWV8uDLLlc8P8eXBVXUOnuRQftvUcq8VuNvbus72LrUc5inbuvcIy7cf4OFR3QmyeamMMXWkW5tIkju24J3lO7ljeCIidfv9Y3eOn4Z3lu8kNFi4foBfjAg2xjQg44Z0YPu+fJZs21/nn22J4xQdKyrlw5VZjOwVR8umjdwOxxjTwIzqFUdURChvL6v7O8ktcZyiz1J3c7iwhHGDbQiuMabuhYcGc13/9s7CcYV1+tmWOE7RO8t30jm2CYMTo90OxRjTQI0d3IGSMuX9lLq9SG6J4xSs353Hqp2HuGlwxzq/KGWMMSd0jm3K0M4teXf5zjpd5MkSxyl4d/lOwkKCuLZ/ZbOgGGNM3blpcAeyDh5j4VbfzpBRniWOGiooKuHjVbu5/Mw4oiLC3A7HGNPAXdyzDS2bhDFzed1dJLfEUUOfpWZz9HgJY2xeKmOMHwgLCeK6Ae35amMOOYfr5iK5JY4aete5KD4woYXboRhjDAA3DoyntEx5v47uJLfEUQOb9hxm1c5DjB3UwS6KG2P8RqfYpgzpFM3MFTspq4OL5JY4amDm8kzCgoO4pr/dKW6M8S9jB3Ug88Axvtu2z+efZYnDS4XFpfzrhyxG9mpDdBO7KG6M8S+XnNGGqIhQZi7PrL7wabLE4aU5a7M5XFhiizUZY/xSeGgw1zp3kuceOe7Tz7LE4aWZKzJJaBnBkE52p7gxxj+NHRRPSZny0SrfXiS3xOGF9NyjLN9+gBsH2kVxY4z/6tIqkgEdWzBzRaZP1yS3xOGFWSlZBAcJ1w6wO8WNMf7txuR40nPzWbnjoM8+wxJHNYpLy/hgZRbnd29Fq8hwt8MxxpgqXdY7jiZhwcxc4buL5JY4qjF/Uw77jh7nxuR4t0MxxphqNWkUwhV92vJ5ajZHCot98hk+TRwiMlJENotImog8XMl+EZEpzv5UEelfXV0RiRaReSKy1Xn26S3c763IpFVkI87tFuvLjzHGmFpz48B4jhWX8umabJ8c32eJQ0SCgReBUUBPYKyI9KxQbBSQ5DwmAS97Ufdh4GtVTQK+dt77xJ68QuZvzuH65PaEBFvnzBgTGPrGR9GtdSTvpfjmdJUvvw0HAWmqmq6qRcBMYHSFMqOBGeqxFIgSkbhq6o4G3nBevwFc5asGfPhDFmUKN9hpKmNMABERbhgYz5rMQ2zac7jWj+/LxNEOKJ/uspxt3pSpqm5rVc0GcJ5bVfbhIjJJRFJEJCU399TmqY+NbMSNyfF0bNnklOobY4xbru7XjnO6xlJcUvvDckNq/Yj/VdkNDxVbcLIy3tStkqpOA6YBJCcnn9JP7obkeOttGGMCUnSTMGbcPsgnx/ZljyMLKP+t2x7Y7WWZqurudU5n4Tzn1GLMxhhjquHLxLECSBKRRBEJA8YAsyuUmQ2Md0ZXDQHynNNPVdWdDUxwXk8APvFhG4wxxlTgs1NVqloiIvcAc4Fg4FVVXS8ik539U4E5wKVAGlAA3FZVXefQTwCzRGQisBO43ldtMMYY82Piy/lM/EVycrKmpKS4HYYxxgQUEVmpqskVt9vNCcYYY2rEEocxxpgascRhjDGmRixxGGOMqZEGcXFcRHKBHadYPQbw/erv/sXa3DBYmxuG02lzR1X90QyvDSJxnA4RSalsVEF9Zm1uGKzNDYMv2mynqowxxtSIJQ5jjDE1YomjetPcDsAF1uaGwdrcMNR6m+0ahzHGmBqxHocxxpgascRhjDGmRixxOERkpIhsFpE0EfnROubO1O9TnP2pItLfjThrkxdtHue0NVVEvheRPm7EWZuqa3O5cgNFpFRErqvL+GqbN+0VkXNFZLWIrBeRb+s6xtrmxe91cxH5VETWOG2+zY04a5OIvCoiOSKy7iT7a/f7S1Ub/APP1O3bgE5AGLAG6FmhzKXAF3hWJxwCLHM77jpo81CghfN6VENoc7ly3+CZ9v86t+P28b9xFLAB6OC8b+V23HXQ5l8DTzqvY4EDQJjbsZ9mu88B+gPrTrK/Vr+/rMfhMQhIU9V0VS0CZgKjK5QZDcxQj6VA1ImVCANUtW1W1e9V9aDzdimelRgDmTf/zgD3Ah8S+KtLetPem4B/qepOAFVtCG1WIFJEBGiKJ3GU1G2YtUtVF+Jpx8nU6veXJQ6PdkBmufdZzraalgkkNW3PRDx/sQSyatssIu2Aq4GpdRiXr3jzb9wVaCEiC0RkpYiMr7PofMObNr8A9MCzHPVa4GeqWlY34bmmVr+/fLYCYICRSrZVHKfsTZlA4nV7ROQ8PIljmE8j8j1v2vx34CFVLfX8QRrQvGlvCDAAuABoDCwRkaWqusXXwfmIN22+BFgNnA90BuaJyCJVPezj2NxUq99fljg8soD4cu/b4/lrpKZlAolX7RGR3sA/gVGqur+OYvMVb9qcDMx0kkYMcKmIlKjqx3USYe3y9vd6n6rmA/kishDoAwRq4vCmzbcBT6jn5H+aiGwHugPL6yZEV9Tq95edqvJYASSJSKKIhAFjgNkVyswGxjujE4YAeaqaXdeB1qJq2ywiHYB/AbcE8F+g5VXbZlVNVNUEVU0APgDuCtCkAd79Xn8CDBeREBGJAAYDG+s4ztrkTZt34ulhISKtgW5Aep1GWfdq9fvLehyAqpaIyD3AXDyjMl5V1fUiMtnZPxXPCJtLgTSgAM9fLQHLyzY/BrQEXnL+Ai/RAJ5Z1Ms21xvetFdVN4rIv4FUoAz4p6pWOqQzEHj5b/xH4HURWYvnFM5DqhrQU62LyLvAuUCMiGQBvwNCwTffXzbliDHGmBqxU1XGGGNqxBKHMcaYGrHEYYwxpkYscRhjjKkRSxzGGGNqxBKHCWjODLarRWSdiLzv3ItQG8edIyJRIpJQxYyjC0TEZ8OTReR1N2bnFZHrRWSjiMyv6882gcEShwl0x1S1r6r2AoqAybVxUFW9VFUP1caxAtBEPDc+nud2IMY/WeIw9ckioIuIRIvIx866A0udaVMQkRFO72S1iKwSkUgRiRORheV6LcOdshkiEuMcN0RE3nCO90FlvRoRuVhElojID07Pp2mF/T1EZHm59wkikuq8fkxEVjifP00qmSSrfDwikiwiC5zXTcSzFsMKp02jne1niMhyp12pIpJUyTHHisha53OfPBELnjnJporI0xXKXy0iXzl3H8eJyBYRaePlv42pRyxxmHpBRELwrBmyFngcWKWqvfGsvTDDKfYL4G5V7QsMB47hmVZ8rrOtD57J7yrqBkxzjncYuKvCZ8cAvwEuVNX+QArwYPkyqroRCBORTs6mG4FZzusXVHWg02tqDFxeg6Y/CnyjqgOB84CnRaQJnp7X8067kvHMVVQ+5rbAk3gm+usLDBSRq1T1D07841T1lxXa8BGwB7gbmA78TlX31CBWU09Y4jCBrrGIrMbzZbcTeAXPX8xvAqjqN0BLEWkOfAf8TUTuA6JUtQTP3Ea3icjvgTNV9Ugln5Gpqt85r9/ix7MEDwF6At85sUwAOlZynFnADc7rG4H3nNfnicgyZwqM84EzvG8+FwMPO5+7AAgHOgBLgF+LyENAR1U9VqHeQGCBquY6P4e38SwGVJ17gUeA46r6bg3iNPWIzVVlAt0x56/q/6jsVA+gqvqEiHyOZ86epSJyoaouFJFzgMuAN0XkaVWdUbFuNe8FmKeqY6uJ9T3gfRH5lxPPVhEJB14CklU100lg4ZXULeG/f+iV3y/Ataq6uUL5jSKyzGnXXBG5w0mi5eudinZ45rRqLSJBDWAdC1MJ63GY+mghMA4862njmTb8sIh0VtW1qvoknh5KdxHpCOSo6nQ8vZXK1mLuICJnOa/HAosr7F8KnC0iXZzPjBCRrhUPoqrbgFLgt/y3t3EiCexzroucbBRVBp51MwCuLbd9LnDviWQpIv2c505AuqpOwTMzau8Kx1sGjBCRGBEJdtpV5XrjzunA1/Cc3ttIhdNxpuGwxGHqo98Dyc7F5yfwnDoCuN+5ELwGz/WNL/DMKLpaRFbh+UJ+vpLjbQQmOMeLBl4uv1NVc4FbgXedMkvxrO9QmfeAm3Gubzgjt6bjuTbzMZ5TZ5V5HHheRBbhST4n/BHPLKipzrDhPzrbbwTWOaewuvPf6zwnYs7Gc8ppPp51uX9Q1U9O8tkn/BpYpKqL8CSNO0SkRzV1TD1ks+MaY4ypEetxGGOMqRFLHMYYY2rEEocxxpgascRhjDGmRixxGGOMqRFLHMYYY2rEEocxxpga+X8ELT8jjAc6ogAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior.plot()\n", "\n", "plt.xlabel('Possible values of x')\n", "plt.ylabel('Probability')\n", "plt.title('Posterior')\n", "\n", "posterior.idxmax()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1\n", "\n", "Go back and run the update again for the following outcomes:\n", "\n", "* Two heads, one tails.\n", "\n", "* 7 heads, 3 tails.\n", "\n", "* 70 heads, 30 tails.\n", "\n", "* 140 heads, 110 tails." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A better prior\n", "\n", "Remember that this result is based on a uniform prior, which assumes that any value of $x$ from 0 to 100 is equally likely.\n", "\n", "Given what we know about coins, that's probably not true. I can believe that if you spin a lop-sided coin on edge, it might be somewhat more likely to land on heads or tails. \n", "\n", "But unless the coin is heavily weighted on one side, I would be surprised if $x$ were greater than 60% or less than 40%.\n", "\n", "Of course, I could be wrong, but in general I would expect to find $x$ closer to 50%, and I would be surprised to find it near 0% or 100%.\n", "\n", "We can represent that prior belief with a triangle-shaped prior.\n", "Here's an array that ramps up from 0 to 49 and ramps down from 50 to 0." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "ramp_up = np.arange(50)\n", "ramp_down = np.arange(50, -1, -1)\n", "\n", "ps = np.append(ramp_up, ramp_down)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I'll put it in a `Series` and normalize it so it adds up to 1." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "triangle = pd.Series(ps, xs)\n", "triangle /= triangle.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the triangle prior looks like." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3q0lEQVR4nO3debxV8/7H8de7OUKlVJqRIUPJkQwhRGVIhsoPGXK7XeLiDrjuJcO913i5rrFMxSVJyJW6GTNFpyQq0YBSdC4pQ8qpz++P7zpspzPsU3vvdfben+fjsR9n77W+a63Pt+F89lrruz5fmRnOOedcKtSIOwDnnHO5w5OKc865lPGk4pxzLmU8qTjnnEsZTyrOOedSxpOKc865lPGk4lwSJHWXND8Dx2knySTVSvexouPNkXRoJo7l8oP8ORWXryR9m/BxC2AtsD76/Gsz+3cMMbUDFgO1zaw408d3bnNl5NuQc9WRmTUoeS/pY+AcM3u+dDtJtXLtF/zm9ikX/0xcavjlL+dKkXSopKWSLpH0OfBAybKENpdKWijpG0lzJfVLWHempNck3SRppaTFknonrG8vaWq07fOS7pD0cDmxbCPpPknLJX0m6VpJNctpO1zSOEmPRfueKalTwvqPoz7NBr6TVCtadkS0vq6kWyUti163Sqpb3p/JZv4xuxzlScW5sjUHGgNtgSFlrF8IdAe2Aa4CHpbUImH9fsB8oAlwA3CfJEXrHgHeBrYFhgOnVxDHKKAY2AnYGzgSOKeC9n2Bx6PYHwGeklQ7Yf0pwNFAwzLONC4HugGdgU5AV+DPCesr+zNxzpOKc+XYAFxpZmvNbE3plWb2uJktM7MNZvYY8BHhl3CJT8xspJmtJySGFkAzSW2AfYErzGydmb0GTCgrAEnNgN7AhWb2nZmtAG4BBlYQ9wwzG2dmPwL/AOoREkWJ28xsSVl9Ak4FrjazFWZWREiWiQmvwj8T58DvqThXniIz+6G8lZIGARcD7aJFDQhnJSU+L3ljZt9HJyklbb4ys+8T2i4BWpdxmLZAbWD5zyc51Ijal+endWa2Ibpkt31Z68uwPfBJwudPSm1b4Z+Jc+BJxbnylDssUlJbYCRwOPCmma2XNAtQedskWA40lrRFQmIpK6FASABrgSZVuCn+074k1QBaAcsS1lc03HMZIZHNiT63qcK2zgF++cu5TbEl4RdsEYCks4A9ktnQzD4BCoHhkupI2h84tpy2y4H/AjdL2lpSDUk7SjqkgkPsI+mE6DmXCwlJaVqS/XoU+LOkppKaAFcAZQ4gcK48nlScqyIzmwvcDLwJfAHsCbxehV2cCuwPfAlcCzxG+OVflkFAHWAusBIYR7g/U56ngQFR29OBE6L7K8m4lpDwZgPvATOjZc4lzR9+dC5mkh4DPjCzKzdzP8OBnczstJQE5twm8DMV5zJM0r7RZawaknoRhgE/FXNYzqWE36h3LvOaA+MJz6ksBX5jZu/EG5JzqeGXv5xzzqWMX/5yzjmXMnl9+atJkybWrl27uMNwzrmsMmPGjP+ZWdOy1uV1UmnXrh2FhYVxh+Gcc1lF0iflrfPLX84551LGk4pzzrmU8aTinHMuZTypOOecSxlPKs4551ImrUlFUi9J8yUtkHRpGesl6bZo/WxJXaLlrSW9JGmepDmSfpuwTWNJUyR9FP1slLDusmhf8yUdlc6+Oeec21jakko0j/YdhJnrOgKnSOpYqllvoEP0GgLcFS0vBn5nZrsRZq07L2HbS4EXzKwD8EL0mWj9QGB3oBdwZ3lzeTvnnEuPdJ6pdAUWmNkiM1sHjCEUzkvUFxhtwTSgoaQWZrbczGYCmNk3wDygZcI2o6L3o4DjE5aPiaY6XQws4JfTuzqXHcxg9lhYtTTuSJyrsnQmlZb8curSpfycGJJuI6kdsDfwVrSoWTR5UckkRttV4XhIGiKpUFJhUVFRVfrjXGa8PQLG/woePgnWfRd3NM5VSTqTSllTq5auXllhG0kNgCeAC81sdQqOh5mNMLMCMyto2rTMKgPOxWfJdJh8ObToBEUfwDMXhjMX57JEOpPKUn4593bpubIrbCOpNiGh/NvMxie0+UJSi6hNC2BFFY7nXPX13Zfw+JmwdQsY9DT0+BO8NxYK74s7MueSls6kMh3oIKm9pDqEm+gTSrWZAAyKRoF1A1aZ2XJJAu4D5pnZP8rY5ozo/RmE6VNLlg+UVFdSe8LN/7dT3y3n0mDDehh/DnxXBP0fgvqNoPvvYaeeMOky+GxG3BE6l5S0JRUzKwaGAZMJN9rHmtkcSUMlDY2aTQQWEW6qjwTOjZYfSJhf+zBJs6JXn2jddUBPSR8BPaPPmNkcYCxhLu9JwHlmtj5d/XMupV65Hha+CH1ugO07h2U1asAJI6BBcxh7Bnz/VawhOpeMvJ6kq6CgwLxKsYvdR8/Dv0+CTqfA8XeCSt0e/GwG3N8L2h8C/zc2JBvnYiRphpkVlLXO/3U6F6evPw2XvZrtDkffvHFCAWi5D/T6OyyYAq/elPkYnasCTyrOxaV4bbistWE99B8NdbYov23BYNizP7z0t3CZzLlqypOKc3GZ/CdYNhP63gHb7lhxWwmOvRWa7gpPnOMPRrpqy5OKc3GY/ThMvxcOOB86HpfcNnW2hAEP/XyGU7wuvTE6twk8qTiXaSvmwTMXQJsD4PDhVdu2SQfoezt8Vgj//XNawnNuc3hScS6T1n4DYwdBnQZw8gNQs1bV97F7P+h2Lrx9D7w3LvUxOrcZPKk4lylmMOF8+HIBnHQ/bNV80/fV82povR9MuACK5qcuRuc2kycV5zLlrXtgzpNw+BXQvvvm7atmbTj5QahdHx47HdZ+m5IQndtcnlScy4Qlb8N/L4dd+sCBF6Zmn1tvDyfdB19+BM/81gtPumrBk4pz6fbd/0KhyG1awfF3lf2A46ba4dBQePL9cWE0mXMx86TiXDptWA9PDA6Jpf9oqN8w9cc46HfQ4ahQeHKplx1y8fKk4lw6vfx3WPQyHH1TmCMlHWrUgH53h5L5Y88IJfSdi4knFefS5cP/wtQbofNp0GVQeo+1ReNwJvTdijBr5AYv0O3i4UnFuXRY+Un45d5sz3CWkgnb7w29b4CFL4Rk5lwMPKk4l2rFa+HxM8A2QP9RYdhvpuxzJnT6P3j5OljwfOaO61zEk4pzqTbpUlj2ThjpVVmhyFSTQgn97TrCE7+Cr5dk9vgu76U1qUjqJWm+pAWSLi1jvSTdFq2fLalLwrr7Ja2Q9H6pbR5LmA3yY0mzouXtJK1JWHd3OvvmXJnefQwK74cDLoDdjoknhjpbhPsr638MQ5m98KTLoLQlFUk1gTuA3kBH4BRJHUs1602YS74DMAS4K2Hdg0Cv0vs1swFm1tnMOgNPAOMTVi8sWWdmQ0tv61xafTE3PITY9kA4/Mp4Y2myU5hF8rPC8NClcxmSzjOVrsACM1tkZuuAMUDfUm36AqMtmAY0lNQCwMymAuVOyi1JQH/g0bRE71xV/LAaxp4OdbcKdb02pVBkqnU8DvYfBm+P8MKTLmPSmVRaAokXdJdGy6rapjzdgS/M7KOEZe0lvSPpFUllFleSNERSoaTCoqKiJA/lXAXMYMIw+GpxqMe1OYUiU+2I4dBm/1B4csUHcUfj8kA6k0pZtShKFydKpk15TuGXZynLgTZmtjdwMfCIpK032rnZCDMrMLOCpk2bJnko5yow7U6Y+zQccSW0OzDuaH6pZm046YFwn2Xs6aH0vnNplM6kshRonfC5FbBsE9psRFIt4ATgsZJlZrbWzL6M3s8AFgI7b1LkziXr02kw5QrY9Zhwc7462rpFuCT35YJwxuKFJ10apTOpTAc6SGovqQ4wEJhQqs0EYFA0CqwbsMrMliex7yOAD8zsp4m6JTWNBgcgaQfCzf9FqeiIc2X6tigqFNk6zDOfykKRqdb+YDjsLzBnfLjH4lyapC2pmFkxMAyYDMwDxprZHElDJZWMzJpI+MW/ABgJnFuyvaRHgTeBXSQtlTQ4YfcD2fgG/cHAbEnvAuOAoWZW7o1+5zbLhvXwxNmwZmX6CkWm2oEXws69YfLlsGR63NG4HCXL41PhgoICKyz0qq5uE7xwNbx6czhD2fu0uKNJ3pqVcM8hsKEYfj0VtmwSd0QuC0maYWYFZa3zJ+qdq6oPJ4eE0mVQdiUUgPqNYMBDoRT/E4O98KRLOU8qzlXFyo9Docjme0HvLC3a2KIT9LkxlOR/+bq4o3E5xpOKc8n68QcYG5Ww7z8aateLN57N0WUQdD4Vpt4AH02JOxqXQzypOJesSZfA8neh3z3QuH3c0WweCfrcBM32CGdeX38ad0QuR3hScS4Zsx6FGQ/CQRfBLr3jjiY1SgpPblgfZowsXht3RC4HeFJxrjJfzIH/XATtukOPP8cdTWptu2MoPLlsJkz+U9zRuBzgScW5ivywCh47HeptU30KRababsfCAefD9Hth9ti4o3FZzpOKc+Uxg6fPCyO+Tn4AGmwXd0Tpc/hwaHNAKN3/xdy4o3FZzJOKc+V58w6Y9wz0vAraHhB3NOlVs1ZInHUahBFuXnjSbSJPKs6V5ZM3Q6HI3Y4Nc5Lkg62ah0t8Xy2Ep4d54Um3STypOFfatytCochGbat/ochUa98dDr8C5j4Fb/mM3K7qPKk4l2h9MYw7O9yg7/9QuEGfbw68EHY5Gv77Z/j0rbijcVnGk4pziV66Fj5+FY65BZrvEXc08ZDCMONtWocztm99hlSXPE8qzpX4YCK8dgvscyZ0PiXuaOJVv2F4MHLNV1540lWJJxXnIMwv/+TQUGyx1/VxR1M9tNgrlHJZ/Aq8/Pe4o3FZwpOKcyWFIqXsLxSZal1OD+X9p94YSv47V4m0JhVJvSTNl7RA0qVlrJek26L1syV1SVh3v6QVkt4vtc1wSZ9JmhW9+iSsuyza13xJR6Wzby6HPPcH+Hw2nDACGrWLO5rqp89N0HxPGD8EVn4SdzSumktbUonmi78D6A10BE6R1LFUs96EueQ7AEOAuxLWPQj0Kmf3t5hZ5+g1MTpeR8I0w7tH291ZMme9c+V6598wczR0/x3s7N9DylS7fjiDMwtndD/+EHdErhpL55lKV2CBmS0ys3XAGKBvqTZ9gdEWTAMaSmoBYGZTgarMMd8XGGNma81sMWHe+66b3QuXuz5/D569GNofDD0ujzua6q3xDtDvLlg+CyZtdNHBuZ+kM6m0BJYkfF4aLatqm7IMiy6X3S+pUVX2JWmIpEJJhUVFPlQyb/2wKnzrrt8ITrwfavhJbaV2PRoO/C3MeADeHRN3NK6aSmdSKesx5NJ1H5JpU9pdwI5AZ2A5cHNV9mVmI8yswMwKmjZtWsmhXE4yg6fODRNTnfwgNPB/B0k77ApoexA8c2GYEsC5UtKZVJYCrRM+twKWbUKbXzCzL8xsvZltAEby8yWuKu/L5ak3/gUf/Ad6Xg1tusUdTXapWSvUB6u3dZgS4IfVcUfkqpl0JpXpQAdJ7SXVIdxEn1CqzQRgUDQKrBuwysyWV7TTknsukX5AyeiwCcBASXUltSfc/H87FR1xOeTj1+H54dCxL3Q7N+5ostNWzeCkB8KUAE+f54Un3S+kLamYWTEwDJgMzAPGmtkcSUMlDY2aTQQWEW6qjwR++l8u6VHgTWAXSUslDY5W3SDpPUmzgR7ARdHx5gBjgbnAJOA8M/PHgN3Pvvkcxp0V5pc/7vb8KhSZau0OhCOGw7wJMO3OuKNx1Ygsj79lFBQUWGFhYdxhuExYXwyj+8JnM+BXL0Cz3eOOKPuZwWOnwYeT4Mxn/VJiHpE0w8wKylrnT9S7/PDi1fDJa3DsPz2hpEpJ4cmGbbzwpPuJJxWX+z54Fl7/JxScDZ0GxB1Nbqm3TVR4ciU8cbYXnnSeVFyO+2oRPPkb2H5v6HVd3NHkpuZ7wtH/gMVT4aW/xh2Ni5knFZe7flzzc6HIk0dBrbpxR5S79j4V9j4dXr0Z5k+KOxoXI08qLndN/H0oxXLCyDA1sEuvPjdC873gySFhuLHLS55UXG6a+RC88zAc/AfY+ci4o8kPJYUnwQtP5jFPKi73LJ8dzlJ2OBQOvSzuaPJL4/bQ7x5Y/i5MuiTuaFwMPKm43LLmaxh7OtRvDCfe54Ui47BLbzjoIpjxIMx6NO5oXIZ5UnG5o6RQ5Kql0H8UbNkk7ojyV48/Q7vu8J+L4PP3K2/vcoYnFZc7Xv8nzH8WjrwWWvtUOrH6qfDkNuHM8YdVcUfkMsSTissNH78GL1wNu/eD/YZW3t6lX4Pt4OQHwhTET53rhSfzhCcVl/2++RwePyvMTnjcv7xQZHXS9gDoeVWYauDN2+OOxmWAJxWX3db/GBLKum9hwENQd6u4I3Kl7T8MdjsWplwJn7wRdzQuzTypuOz2wlXw6RuhUOR2u8UdjSuLBH3vCA+gPn4WfPNF3BG5NPKk4rLXvGfCLI4Fg2Gv/nFH4ypSbxvo/1C4Yf/E4DAVgctJaU0qknpJmi9pgaRLy1gvSbdF62dL6pKw7n5JKyS9X2qbGyV9ELV/UlLDaHk7SWskzYped6ezby5mXy4MN3+37wK9/h53NC4ZzfeAY26Bj1+Fl66NOxqXJmlLKpJqAncAvYGOwCmSOpZq1psw7W8HYAhwV8K6B4FeZex6CrCHme0FfAgkPjK90Mw6Ry8fApSr1n0f5kevUTM8j+KFIrNH51NgnzPhtVvgg4lxR+PSIJ1nKl2BBWa2yMzWAWOAvqXa9AVGWzANaFgyB72ZTQW+Kr1TM/tvNFUxwDSgVdp64KofM3j2d7BiLpxwb5ggymWXXtdDi07w5FD4anHc0bgUS2dSaQksSfi8NFpW1TYVORt4LuFze0nvSHpFUveyNpA0RFKhpMKiIp+pLuvMHAXvPgKH/BE6HBF3NG5T1K4XCk9K4cHIH9fEHZFLoaSSiqRjJFU1AZX1sEDpp5+SaVNeTJcDxcC/o0XLgTZmtjdwMfCIpK032rnZCDMrMLOCpk2bJnMoV10smwUT/wg7HgaHeLHCrNaoHZwwIkxNMPEPcUfjUijZRDEQ+EjSDZKSHbe5FGid8LkVsGwT2mxE0hnAMcCpZuExXTNba2ZfRu9nAAuBnZOM1VV3a1aGb7VbNg2XvbxQZPbb+Sjo/jt4J5qmwOWEpJKKmZ0G7E34Rf2ApDejy0gVPWk2Heggqb2kOoTENKFUmwnAoGgUWDdglZktrygWSb2AS4DjzOz7hOVNo8EBSNqBcPN/UTL9c9Xchg3h+vvq5VGhyG3jjsilSo/Lof3B4T7Z8tlxR+NSIOlLWma2GniCcMO9BdAPmCnp/HLaFwPDgMnAPGCsmc2RNFRSycisiYRf/AuAkcC5JdtLehR4E9hF0lJJg6NVtwNbAVNKDR0+GJgt6V1gHDDUzDa60e+y0Ou3wIeT4Ki/QauCuKNxqVSjJpx4P9RvFCb2WvN13BG5zSRLosibpOOAs4AdgYeAUWa2QtIWwDwzy8q5WgsKCqywsDDuMFxFFk+F0X1DocgT7/O6Xrnq07fgwT6wcy8Y8LD/PVdzkmaYWZnf8JI9UzkJuMXM9jKzG81sBUB0+ensFMXp3C+tXgbjzoZtO8Cxt/kvmlzWZj/oeU0oPPnGbXFH4zZDskllefTcyE8kXQ9gZi+kPCrnfioU+X1UKLJB3BG5dOv2G+h4PDw/PExl4LJSskmlZxnLeqcyEOd+4fnhsGQaHHcbNN0l7mhcJkhh6oLGO0SFJz+POyK3CSpMKpJ+I+k9YNeo1lbJazHgQzVcesx9Osy90XUI7HlS3NG4TKq3dSg8ue7bcOnTC09mncrOVB4BjgWejn6WvPaJhhk7l1r/WwBPnQctC+DIv8YdjYtDs45wzK3wyethagOXVSpLKmZmHwPnAd8kvJDUOL2hubyz7vswrLRmbTj5QahVJ+6IXFw6DYCCs8NN+3n/iTsaVwW1Kln/COHJ9RmE8imJw28M2CFNcbl8YwbPXhwKRZ72BDRsXfk2Lrf1ug6WvROmONhuN9h2x7gjckmo8EzFzI6JfrY3sx2inyUvTygudWY8CO8+CodeCjsdHnc0rjqoVRdOHgU1asDYM7zwZJao8EwlcdKsspjZzNSG4/LSZzPhuT/CjofDwX+MOxpXnTRqC/1GwCMnw7O/h+PviDsiV4nKLn/dXME6Aw5LYSwuH33/VfgW2qAZnDAyfCt1LtHOR4YvG1NvCA9JdhkUd0SuAhUmFTPrkalAXB4qKRT5zXI4e7IXinTlO/RSWDo9nK206BRerlqq7DmVw6KfJ5T1ykyILme9djN8NDnMMd9qn7ijcdVZjZpw4r2wZZMwlfSalXFH5MpR2bWGQ6Kfx5bxOiaNcblct+hleOlvsOfJsO85cUfjssGWTcJQ89WfwZO/CWe6rtqp7PLXldHPszITjssLqz6DcYNDochjbvVCkS55rbuGh2InXQKv3wrdL447IldKstMJbyvpNkkzJc2Q9E9JfgHcVV3xOnj8zDA81AtFuk2x36/DVAgvXhOmRnDVSrJDbcYARcCJhDL4RcBj6QrK5bApV8DSt71QpNt0JYUnt90p1AdbXeFksS7Dkk0qjc3sGjNbHL2uBRpWtpGkXpLmS1og6dIy1is6A1oQFarskrDufkkrJL1fapvGkqZI+ij62Shh3WXRvuZLOirJvrlMeX88vHUX7DfUC0W6zVN3K+g/GtZ9B+POClMluGoh2aTykqSBkmpEr/7AsxVtEM0XfwehRH5H4BRJHUs1602YS74DMAS4K2Hdg0CvMnZ9KfCCmXUAXog+E+17ILB7tN2dJXPWu2qg6EOYcD606homY3Juc223W5i87dM3w1QJrlqobEjxN5JWA78m1AFbF73GABdVsu+uwAIzW2RmJdv0LdWmLzDagmlAQ0ktAKJJwcqaY74vMCp6Pwo4PmH5GDNba2aLCfPed60kRpcJ674LhSJr1fVCkS619opGD755O8ydEHc0jsprf21lZltHP2uYWa3oVcPMtq5k3y2BJQmfl0bLqtqmtGZmtjyKbzmwXVX2JWmIpEJJhUVFRZUcym02M3jmQij6IDxnsE1lf73OVdFRf4OW+8DT58GXC+OOJu8lXRNDUiNJXSUdXPKqbJMyltkmtElWUvsysxFmVmBmBU2bNt3EQ7mkFd4H742FHn+CHb2qj0uDkjPgGjXDg5Hrvo87oryW7JDic4CpwGTgqujn8Eo2Wwok1i9vBSzbhDalfVFyiSz6uWIz9uXS6bMZMOky2KkndP993NG4XNawDZxwb5g64dnfhTNkF4tkz1R+C+wLfBLVA9ubMKy4ItOBDpLaS6pDuIle+qLnBGBQNAqsG7Cq5NJWBSYAZ0TvzyDMSlmyfKCkupLaE27+v51E31w6/FQosjmcMMILRbr063AEHHIJvPsIzBxVeXuXFsn+T//BzH4AkFTXzD4AKnzIwMyKgWGEs5p5wFgzmyNpqKShUbOJwCLCTfWRwLkl20t6FHgT2EXSUkmDo1XXAT0lfQT0jD5jZnOAscBcYBJwnpmtT7J/LpU2bIDxQ+DbL6D/g7CFTxLqMuSQP4bLrBP/CMtmxR1NXpIlcZoo6UngLOBCQrn7lUBtM+uT1ujSrKCgwAoLC+MOI/e8cgO89Fc4+h+w7+DK2zuXSt99CfccHM6Ofz0V6jeqfBtXJZJmmFlBWeuSOlMxs35m9rWZDQf+AtzHz0N5nfvZwhdDoci9ojnGncu0LbeF/qPCk/ZPDvXCkxlWldFfXSRdAOwFLI2ePXHuZ6s+gyfOgaa7wjG3eKFIF59WBWGo8YeT4PVb4o4mryQ7+usKwoOG2wJNgAck/TmdgbksU7wOHj8DiteGQpF1tow7Ipfvuv4K9jgRXrzWC09mULJnKqcA+5rZlVE5/G7AqekLy2WdKX8JM/P1vQOadIg7GufCmfKxt4UpFsadDav9CYNMSDapfAzUS/hcF/BHV13w/hPw1t3Q7VzY/fi4o3HuZ3UbhDPndd/D4154MhMqq/31L0m3AWuBOZIelPQA8D7wbSYCdNVc0Ycw4QJovR/0vDruaJzbWNNdwlQLS6bBlCvjjibnVTjzI1Ay3nYG8GTC8pfTEo3LLmu/hbGnQ616cNIDULN23BE5V7Y9T4Ilb8G0O8LskX5GnTaVTSf802Op0VPxO0cf55uZn0fmMzP4z4Xwvw/h9Ce9UKSr/o78K3w2E54eBs32gCY7xR1RTkp29NehwEeE+VHuBD5MoqCky2XT74X3Hg+FInc4NO5onKtcrTqh8GTN2uEMe913cUeUk5K9UX8zcKSZHWJmBwNHAT74O18tjQpFdjgKDvpd3NE4l7yGrcMUDCvmwX8u8sKTaZBsUqltZvNLPpjZh4BfQM9H330ZnkfZqgX0u9sLRbrss9PhcOilMPsxmPFA3NHknMpu1JeYIek+4KHo86mEm/cun2xYD+N/FQpFnj3ZC0W67HXwH8NzVc9dAi06Q8sucUeUM5L9mjkUmANcQCiDPzda5vLJ1Bth4QvQ+3r/T+iyW40acMJIaNAsTNHwfVkzl7tNUWlSkVQDmGFm/zCzE6LikreY2doMxOeqiwXPw8vXwV4DYZ+z4o7Guc23RWM4eRR8szxM1eCFJ1Oi0qRiZhuAdyW1yUA8rjr6ekkoFLldRzjmH14o0uWOVvtA7+tgwRR49ea4o8kJyd5TaUF4ov5t4KdxeGZ2XFqictVH8dpwY359MfQf7YUiXe4pGAyfvhXmAGpVADv2iDuirJZsUrlqU3YuqRfwT6AmcK+ZXVdqvaL1fYDvgTPNbGZF20p6jJ9nnWwIfG1mnSW1I8wwWTJKbZqZ+X2fzTX58jDXfP+H/GExl5skOPZW+Pw9eGIw/PpVf5h3M1SYVCTVI9yQ3wl4D7gvmia4UpJqEh6W7AksBaZLmmBmcxOa9SbMJd8B2A+4C9ivom3NbEDCMW4GViXsb6GZdU4mPpeE98bB9JGw/zDo6CelLofV2TIUnhxxKDx+Jpz5bHhY0lVZZfdURgEFhITSm/AQZLK6AgvMbFE0odcYoG+pNn2B0RZMAxpKapHMttFZTn/g0SrE5JK14gOYcD602R+OGB53NM6lX5MO0Pd2WPo2TLki7miyVmVJpaOZnWZm9wAnAd2rsO+WwJKEz0ujZcm0SWbb7sAXZvZRwrL2kt6R9IqkMmOVNERSoaTCoqKi5HuTT9Z+E8pY1GnghSJdftm9H+z3G3jrLnh/fNzRZKXKkspPRSOTveyVoKwhQqVrIpTXJpltT+GXZynLgTZmtjdwMfCIpK032onZCDMrMLOCpk2blht83jILpey/XAAn3Q9bt4g7Iucyq+fV0KprOFMv+jDuaLJOZUmlk6TV0esbYK+S95JWV7LtUqB1wudWQOmp18prU+G2kmoBJwCPlSwzs7Vm9mX0fgZhErGdcVXz9giYMx4O+wu0r8qJqXM5oqTwZK26XnhyE1SYVMysppltHb22MrNaCe83OgsoZTrQQVL7qGz+QGBCqTYTgEEKugGrzGx5EtseAXxgZktLFkhqGt3gR9IOhJv/iyr9E3A/WzI9jPbauTcceGHc0TgXn21awon3QdF8eOZCLzxZBckOKa4yMyuWNAyYTBgWfL+ZzZE0NFp/NzCRMJx4AWFI8VkVbZuw+4FsfIP+YOBqScXAemComXnthWR997/wPMrW20O/u7xQpHM79oAel8NL10Kb/WDfc+KOKCvI8jgDFxQUWGFhYeUNc92G9fDwifDJG3DOFGjRKe6InKseNmyARwfAopfh7EnQcp+4I6oWJM0ws4Ky1vnXUQevXA+LXoI+N3pCcS5RjRrQ7x5o0NwLTybJk0q+++h5eOUG6HwqdBkUdzTOVT9bNIb+o8KUD+N/5YUnK+FJJZ99/SmMPwea7Q59bvJCkc6Vp2UX6HVdqNY99ca4o6nWPKnkq+K14XR+w/qoUOQWcUfkXPVWcDbsNQBe/jsseCHuaKotTyr5avKfYNlMOP5O2HbHuKNxrvqT4JhboOmuYSqIr5dUvk0e8qSSj2aPhen3wgHnw27Hxh2Nc9mjzpYw4GFY/2MoPFm8Lu6Iqh1PKvlmxTx45rfQ9kA4fHjc0TiXfZrsFApPflYI/7087miqHU8q+eSH1fDYaVGhyPuhZtqefXUut+1+PHQ7L5Q1em9c3NFUK55U8oVZKJD31WI4+QHYqnncETmX3XpeBa27hQKsKz6IO5pqw5NKvph2F8x9Cg6/AtodFHc0zmW/mrXDF7Q6W8DYQbD227gjqhY8qeSDT6fBlL/ALkfDgb+NOxrncsfW24fCk19+BM9c4IUn8aSS+74tCqNUtmkVhg/7A47OpdYOh4TCk+8/AW+PjDua2HlSyWUb1sMTg0O9ov6joX7DuCNyLjcddDHs3Cs8/7VketzRxMqTSi576W+w+BU4+iYvFOlcOtWoAf3uDjOlPn4mfPdl3BHFxpNKrvpwMrx6E3Q+zQtFOpcJ9RuFKwLfrQg19TasjzuiWKQ1qUjqJWm+pAWSLi1jvSTdFq2fLalLZdtKGi7pM0mzolefhHWXRe3nSzoqnX2r1lZ+AuOHQPM9w1mKcy4ztt8bet8AC18M1b/zUNqefoum9r0D6EmYc366pAlmNjehWW/CtL8dgP2Au4D9ktj2FjP7xW9LSR0JM0LuDmwPPC9pZzPLr68LP/4QhjeahW9NtevHHZFz+WWfM2HJ22Geolb7Qocj4o4oo9J5ptIVWGBmi8xsHTAG6FuqTV9gtAXTgIaSWiS5bWl9gTFmttbMFhOmKO6ayg5lhUmXwvJZYUrgxjvEHY1z+UeCo2+G7TqGy2B5VngynUmlJZD4p7k0WpZMm8q2HRZdLrtfUqMqHA9JQyQVSiosKiqqSn+qv3fHwIwHwrMoux4ddzTO5a86W4QrBeuL4fEzwlQTeSKdSaWsByJKPxlUXpuKtr0L2BHoDCwHbq7C8TCzEWZWYGYFTZs2LWOTLPXFHHjmQmh7EBx2RdzROOea7BSeDftsBkzOn8KT6UwqS4HWCZ9bAcuSbFPutmb2hZmtN7MNwEh+vsSVzPFy0w+r4bHTod7WXijSueqk43Gw/zCYPhJmPx53NBmRzqQyHeggqb2kOoSb6BNKtZkADIpGgXUDVpnZ8oq2je65lOgHvJ+wr4GS6kpqT7j5/3a6OldtmMHT58HKj+HkB2GrZnFH5JxLdMRwaLN/KOOSB4Un05ZUzKwYGAZMBuYBY81sjqShkoZGzSYCiwg31UcC51a0bbTNDZLekzQb6AFcFG0zBxgLzAUmAeflxcivaXfCvAnhH27bA+KOxjlXWs3acNIDYcqJsafD2m/ijiitZHlcAK2goMAKCwvjDmPTfToNHjw6lIcY8LDX9XKuOlv8Kow+DjoeHy5TZ/H/V0kzzKygrHX+RH22KikU2bCNF4p0Lhu07w6H/QXmjA+Te+UoTyrZaMN6eOJsWLMyDFust03cETnnknHghbBz7zAabElu3vL1pJKNXvorLJ4KR/8jlGJxzmWHGjXCg8lbbx8Vnvxf3BGlnCeVbDN/Erx6cygSufepcUfjnKuq+o1gwEMhoTwxOOcKT3pSySYrP4Ynh0DzvaD3jXFH45zbVC06QZ8bYdHL8PJ1cUeTUp5UskVJoUiICkXWizce59zm6TIIOp8KU2+Aj6bEHU3KeFLJFs/9EZa/C/3ugcbt447GObe5JOhzEzTbA8b/Cr7+NO6IUsKTSjaY9QjMHAUHXQS79I47GudcqpQUntywPlyJyIHCk55UqrvP34f/XATtukOPP8cdjXMu1bbdMTxrtuwdmHRZ3NFsNk8q1dkPq0JZh3oNvVCkc7lst2PhgAug8D6YPTbuaDaLJ5Xq6qdCkZ+EQpENtos7IudcOh1+JbQ9EJ75LXwxt/L21ZQnlerqzdth3jPQ82pou3/c0Tjn0q1mrXBFoqTw5A+r445ok3hSqY4+fh2mXAm7HQf7nxd3NM65TNmqOZz8AHy1GCYMC1cssownlermmy9g3FnQqB30vcMLRTqXb9odBEdcCXOfhml3xR1NlXlSqU7WF8O4s8Np74CHwkyOzrn8c8AFsOsxMOUvYYqLLOJJpTp58Rr45DU49lZotnvc0Tjn4iKFKxXbtA6FJ78tijuipKU1qUjqJWm+pAWSLi1jvSTdFq2fLalLZdtKulHSB1H7JyU1jJa3k7RG0qzodXc6+5ZyH0yE12+Ffc6CTgPjjsY5F7f6DcMVizUrs6rwZNqSiqSawB1Ab6AjcIqkjqWa9SbMJd8BGALclcS2U4A9zGwv4EMg8WmhhWbWOXoNJVt8tQieHAotOkOv3Cou55zbDM33hKNvhsWvwEt/izuapKTzTKUrsMDMFpnZOmAM0LdUm77AaAumAQ0ltahoWzP7bzSHPcA0oFUa+5B+P64J5RkkLxTpnNvY3qfB3qfDqzfBh5PjjqZS6UwqLYElCZ+XRsuSaZPMtgBnA88lfG4v6R1Jr0jqXlZQkoZIKpRUWFRUDa5TTvwDfP4enDACGrWNOxrnXHXU58Zw1jJ+SJgCoxpLZ1Ipayxs6UHX5bWpdFtJlwPFwL+jRcuBNma2N3Ax8IikjYZPmdkIMysws4KmTZtW0oU0e+dheOch6P572PmoeGNxzlVftetD/4fCcytjzwhTYVRT6UwqS4HWCZ9bAcuSbFPhtpLOAI4BTjULTweZ2Voz+zJ6PwNYCOyckp6kw/LZ8OzvoP0h0ONPcUfjnKvuGreHfnfD8lkwaaNxT9VGOpPKdKCDpPaS6gADgQml2kwABkWjwLoBq8xseUXbSuoFXAIcZ2bfl+xIUtPoBj+SdiDc/F+Uxv5tujVfh/so9RvDifdBjZpxR+Scywa79oEDL4QZD8C7Y+KOpkxpK3trZsWShgGTgZrA/WY2R9LQaP3dwESgD7AA+B44q6Jto13fDtQFpig8bT4tGul1MHC1pGJgPTDUzL5KV/82WUmhyFVL4MyJ0CDmS3DOuexy2F/gsxnwzIXhPks1e6ZNloW1ZVKloKDACgsLM3vQ1/8JU66Ao/4O+5+b2WM753LDN1/APQdDnS1hyMsZr74haYaZFZS1zp+oz6SPX4fnr4KOx0O338QdjXMuW23VLBSeXPkxPH1utSo86UklU775PBSKbNwejvuXF4p0zm2etgfAEcPDFBlv3hF3ND/xpJIJ64th3GBY+00YFuiFIp1zqXDA+VHhySvgkzfijgbwpJIZL1wVCkUecys0K12pxjnnNpEU5rdv1BYePyvca4mZJ5V0m/cfeOM2KBgMnQbEHY1zLtfU2yZcAflhVSg8ub648m3SyJNKOn25EJ76DWzfBXr9Pe5onHO5qvkecMw/4ONX4aVrYw3Fk0q6/LgmlFOoURP6j4JadeOOyDmXyzr/H3Q5A167JUylERNPKulgFkqwfPE+nDASGraJOyLnXD7ofQO06BSm0vhqcSwheFJJh5mjYda/4eA/QIeecUfjnMsXteuFKTREKAUVQ+FJTyqptmxWKGe/Qw84tPoWfXPO5ahG7aDfCPh8Njz3h4wf3pNKKq1ZGb4dbNkETrzXC0U65+KxSy846OJw1eSdf1fePoU8qaTKhg3w5G9g9Wdw8oMhsTjnXFx6XA7tusOzF4eJADPEk0qqvH4rfPgcHPlXaN017micc/muZi046X6o1zBcQflhVUYO60klFRZPhRevgd37wX6/jjsa55wLGmwXrpys/ASeykzhSU8qm2v1chh3Nmy7kxeKdM5VP233h55Xwwf/gTf+lfbDeVLZHOt/DJWH130XhvHV3SruiJxzbmP7nwe7HQfPDw9TcKRRWpOKpF6S5ktaIGmj8bXRNMK3RetnS+pS2baSGkuaIumj6GejhHWXRe3nSzoqnX0Dwl/Qp2/CsbfBdrul/XDOObdJJOh7RxhuPC69hSfTllSi+eLvAHoDHYFTJJUu0dubMJd8B2AIcFcS214KvGBmHYAXos9E6wcCuwO9gDtL5qxPi7kT4M3bYd9fwV4np+0wzjmXEvW2hgEPwQ+rwyX7NBWeTOeZSldggZktMrN1wBigb6k2fYHRFkwDGkpqUcm2fYFR0ftRwPEJy8eY2VozW0yY9z49w7C+XBjmmW+5Dxz117QcwjnnUq7Z7nDsrWEqjhevScsh0plUWgJLEj4vjZYl06aibZuZ2XKA6Od2VTgekoZIKpRUWFRUVKUO/aRGzZBQTvZCkc65LNNpIOw3NG01CWulZa9BWcOgSo9nK69NMttuyvEwsxHACICCgoJNG1/XqB0MemqTNnXOudj1vj5tu07nmcpSoHXC51bAsiTbVLTtF9ElMqKfK6pwPOecc2mUzqQyHeggqb2kOoSb6BNKtZkADIpGgXUDVkWXtCradgJwRvT+DODphOUDJdWV1J5w8//tdHXOOefcxtJ2+cvMiiUNAyYDNYH7zWyOpKHR+ruBiUAfwk3174GzKto22vV1wFhJg4FPgZOjbeZIGgvMBYqB88xsfbr655xzbmOyDDy2X10VFBRYYWFh3GE451xWkTTDzArKWudP1DvnnEsZTyrOOedSxpOKc865lPGk4pxzLmXy+ka9pCLgk83YRRPgfykKJxvkW3/B+5wvvM9V09bMmpa1Iq+TyuaSVFjeCIhclG/9Be9zvvA+p45f/nLOOZcynlScc86ljCeVzTMi7gAyLN/6C97nfOF9ThG/p+Kccy5l/EzFOedcynhScc45lzKeVCohqZek+ZIWSLq0jPWSdFu0frakLnHEmUpJ9PnUqK+zJb0hqVMccaZSZX1OaLevpPWSTspkfOmQTJ8lHSpplqQ5kl7JdIyplsS/7W0kPSPp3ajPZ8URZ6pIul/SCknvl7M+9b+/zMxf5bwIZfcXAjsAdYB3gY6l2vQBniPMPNkNeCvuuDPQ5wOARtH73vnQ54R2LxKmbDgp7rgz8PfckDCVRJvo83Zxx52BPv8JuD563xT4CqgTd+yb0eeDgS7A++WsT/nvLz9TqVhXYIGZLTKzdcAYoG+pNn2B0RZMAxqWzEyZpSrts5m9YWYro4/TCLNsZrNk/p4Bzgee4OfZRrNZMn3+P2C8mX0KYGbZ3u9k+mzAVpIENCAkleLMhpk6ZjaV0IfypPz3lyeVirUEliR8Xhotq2qbbFLV/gwmfNPJZpX2WVJLoB9wdwbjSqdk/p53BhpJelnSDEmDMhZdeiTT59uB3QhTkb8H/NbMNmQmvFik/PdX2mZ+zBEqY1npMdjJtMkmSfdHUg9CUjkorRGlXzJ9vhW4xMzWhy+xWS+ZPtcC9gEOB+oDb0qaZmYfpju4NEmmz0cBs4DDgB2BKZJeNbPVaY4tLin//eVJpWJLgdYJn1sRvsFUtU02Sao/kvYC7gV6m9mXGYotXZLpcwEwJkooTYA+korN7KmMRJh6yf7b/p+ZfQd8J2kq0AnI1qSSTJ/PAq6zcMNhgaTFwK7A25kJMeNS/vvLL39VbDrQQVJ7SXWAgcCEUm0mAIOiURTdgFVmtjzTgaZQpX2W1AYYD5yexd9aE1XaZzNrb2btzKwdMA44N4sTCiT3b/tpoLukWpK2APYD5mU4zlRKps+fEs7MkNQM2AVYlNEoMyvlv7/8TKUCZlYsaRgwmTBy5H4zmyNpaLT+bsJIoD7AAuB7wjedrJVkn68AtgXujL65F1sWV3hNss85JZk+m9k8SZOA2cAG4F4zK3NoajZI8u/5GuBBSe8RLg1dYmZZWxJf0qPAoUATSUuBK4HakL7fX16mxTnnXMr45S/nnHMp40nFOedcynhScc45lzKeVJxzzqWMJxXnnHMp40nF5ayomvAsSe9Lejx61iIV+50oqaGkdhVUf31ZUtqGWUt6MI5KyZJOljRP0kuZPrbLDp5UXC5bY2adzWwPYB0wNBU7NbM+ZvZ1KvaVhQYTHvzsEXcgrnrypOLyxavATpIaS3oqmjtiWlRuBkmHRGc1syS9I2krSS0kTU042+ketf1YUpNov7UkjYr2N66ssyFJR0p6U9LM6IypQan1u0l6O+FzO0mzo/dXSJoeHX+Eyig8lhiPpAJJL0fvt1SYT2N61Ke+0fLdJb0d9Wu2pA5l7PMUSe9Fx72+JBZCnbe7Jd1Yqn0/Sc9HT2a3kPShpOZJ/t24HOJJxeU8SbUI8768B1wFvGNmexHmzhgdNfs9cJ6ZdQa6A2sIpd8nR8s6EQoNlrYLMCLa32rg3FLHbgL8GTjCzLoAhcDFiW3MbB5QR9IO0aIBwNjo/e1mtm90tlUfOKYKXb8ceNHM9gV6ADdK2pJwxvbPqF8FhPpPiTFvD1xPKKrYGdhX0vFmdnUU/6lm9odSfXgS+Bw4DxgJXGlmn1chVpcjPKm4XFZf0izCL8JPgfsI37QfAjCzF4FtJW0DvA78Q9IFQEMzKybUijpL0nBgTzP7poxjLDGz16P3D7NxxeZuQEfg9SiWM4C2ZexnLNA/ej8AeCx630PSW1HZkMOA3ZPvPkcCl0bHfRmoB7QB3gT+JOkSoK2ZrSm13b7Ay2ZWFP05/Jsw2VNlzgcuA9aa2aNViNPlEK/95XLZmujb+E/KunwEmJldJ+lZQh2kaZKOMLOpkg4GjgYeknSjmY0uvW0lnwVMMbNTKon1MeBxSeOjeD6SVA+4EygwsyVRcqtXxrbF/PwFMXG9gBPNbH6p9vMkvRX1a7Kkc6IEm7jdpmhJqBHWTFKNHJ+HxJXDz1RcvpkKnAph/nVCaffVknY0s/fM7HrCmc2uktoCK8xsJOEsp6z5u9tI2j96fwrwWqn104ADJe0UHXMLSTuX3omZLQTWA3/h57OUkgTxv+g+THmjvT4mzHsCcGLC8snA+SWJVNLe0c8dgEVmdhuhSu1epfb3FnCIpCaSakb9qnB++ugS4wOES4bzKHWJz+UPTyou3wwHCqIb4dcRLkcBXBjdlH6XcD/lOUJ111mS3iH8sv5nGfubB5wR7a8xcFfiSjMrAs4EHo3aTCPMz1GWx4DTiO6nRCPMRhLuBT1FuBxXlquAf0p6lZCYSlxDqEg7Oxr6fE20fADwfnRZbFd+vq9UEvNywmWslwjzuM80s6fLOXaJPwGvmtmrhIRyjqTdKtnG5SCvUuyccy5l/EzFOedcynhScc45lzKeVJxzzqWMJxXnnHMp40nFOedcynhScc45lzKeVJxzzqXM/wPFBsvkBQUYPgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "triangle.plot(color='C1')\n", "\n", "plt.xlabel('Possible values of x')\n", "plt.ylabel('Probability')\n", "plt.title('Triangle prior');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's update it with the data." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "posterior1 = prior * likelihood_heads**140 * likelihood_tails**110\n", "posterior1 /= posterior1.sum()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "posterior2 = triangle * likelihood_heads**140 * likelihood_tails**110\n", "posterior2 /= posterior2.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And plot the results, along with the posterior based on a uniform prior." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7N0lEQVR4nO3deXyU9bnw/881k4SQnSzsq4oIiiJGpZW6Va0Lahdb9dQq9pzHl3vt8pz6dFHaes5Te1rPY1srP6sWba1atafFSrV1oUqrVVAEEVCgCBHIvocsk7l+f9zfoUOYJJNk7plJcr1fL15k7vW6J5P7mu/3vr/XLaqKMcYY01Mg1QEYY4xJT5YgjDHGxGQJwhhjTEyWIIwxxsRkCcIYY0xMliCMMcbEZAnC+EJEPi8if0ri/u4QkRoR2ZesffpFRFaIyB1J2tdOETkrGftKBRH5hojc736eKSIqIhmpjmu4sAQxArg/8v0i0iIilSLyCxHJG8L2lonIr4YSk6o+oqrnDGUb8RKRacBXgXmqOlFElorImgTvY5KIrBSRPe4kM7OX5YpFpLrn/kVkgYisE5E29/+CRMaXSiJyjIg85xL0IQOrRORGEVkrIh0isiLG/I+LyBb33rwkIjP62NdqEfm3eGNT1f9U1biXNwezBDFyXKiqecBC4ETgW6kKZCjf0MQz0M/lDKBWVasGu98eMcSKPww8C3ymn9XvBDb32F4W8HvgV8A44CHg9276SNAF/Ab4117m7wHuAB7sOUNESoHfAt8GioG1wOP+hGkGTFXt3zD/B+wEzop6/V/AH9zPFwGbgAZgNTA3armvAx8CzcBW4OPAuUAn3h99C/C2W7YQeADY69a5Awi6eUuBvwL/DdS5eUuBNVH7+ijwBtDo/v9o1LzVwH+4bewHjohxjLcC212s7wKfctPPcuuEXbyPA+1At3vd4JYbA/wQ2AVUAsuBsW7e6UCFez/2Ab/s473OABSYGWPeR4BXgat7HPs57j2TqGm7gHN72ccK4B7gGXe8fwcOj5p/FPBn915vBT4XNe8C4C2gCdgNLOux7S8AHwC1wDejPzvASXgn6Cb3Ht01wM/hEYD2Mf8OYEWPadcAf4t6net+n0fFWP8/3O+13f1uf+qm3+2OtQlYB3wsap1lwK/czzPd7y4j6nO7w73H/wA+n+q/5XT7Zy2IEcZ1t5wPvCUiRwKPArcAZcAq4GkRyRKROcCNwImqmg98Atipqs8C/wk8rqp5qnqc2/RDQAjvJHA83kkvuul+Mt4f23i8P+TomIrxTnY/BkqAu4BnRKQkarEv4J0s8vFOYD1tBz6Gl6i+A/xKRCap6vPAecAeF++lwLXAq+51kVv/TuBIYIE7hinAbVHbn4j3DXaGi2NARCSId1K/Ee8kFO1oYIO6s5KzwU3vzeXuOMcB23DvqYjk4iWHX+O915cDPxORyLZagSuBIrxkcZ2IfNKtOw+4F++9noz3u5gatc+7gbtVtQA4HK9V4LejgbcjL1S1Fe93fch7o6rfBF4BbnS/2xvdrDfwfq/FeO/LEyKS3ddO3fv4Y+A89/n/KLB+qAcz0liCGDl+JyINwBrgL3gn+UuBZ1T1z6rahfcNeizeH0M33rfqeSKSqao7VXV7rA2LyAS8k/AtqtqqXlfOfwOXRS22R1V/oqohVd3fYxMXAO+r6i/d/EeBLcCFUcusUNVNbn5XzxhU9QlV3aOqYVV9HHgf7xtvv0REgP8FfFlV61S12b0/0fGHgdtVtSNG/PG4Gfi7qq6LMS8Pr+UUrREvGfbmt6r6uqqGgEfwToAAS/AS+S/ce/Um8BRwCYCqrlbVje592oD3BeE0t+4leC3Ll1W1A69bJxy1zy7gCBEpVdUWVX0tzmMfisG8NwdR1V+paq17P36E97meE8eqYeAYERmrqntVdVPcUY8SliBGjk+qapGqzlDV691JbjJR38ZVNYzXFJ+iqtvwWhbLgCoReUxEJvey7RlAJrBXRBpcIvr/8L7BRuzuI7aD4nA+wPsWH8/6iMiVIrI+av/HAKV9rROlDMgB1kWt/6ybHlGtqu1xbq9nbJPxEsQ3e1mkBSjoMa0Ar2ujN9F3Y7XhnUjB+12cHDkOdyyfx2sBISInuwu91SLSiNeairxPk4l6n9239dqo/fwrXitri4i8ISJL+ogvUQbz3hxERL4qIptFpNG9H4X089lwxx5pbe4VkWdE5KgBRT4KWIIY2fbgnVCAA9+kp+H1h6Oqv1bVxW4ZxeuGgUO7SHYDHUCpS0JFqlqgqtHdAH2VBT4oDmd6JI7+1nd3tfwcr/umxHUbvQNIL6v03FYNXr/20VHxF6p3UT+e+PtzEjAJeNfdZns3cJKI7HNdT5uAY937H3Gsmz5Qu4G/RB1Hketuuc7N/zWwEpimqoV411oi+92L9/sHQERy8LqZAFDV91X1crzEfyfwpOuK8dMmINKNGen6OZze35uDfk8i8jG8a0efA8a5z0YjvX82/rkh1edU9Wy8390WvM+YiWIJYmT7DXCBu40wE+9W0A7gbyIyR0TOFJExeBf99uN1O4F3gXJm5G4iVd0L/An4kYgUiEhARA4XkdMO2WNsq4AjReRfRCRDRC4F5gF/iHP9XLwTQzWAiFyN14LoTSUwNXKXkGs5/Rz4bxEZ77YxRUQ+Eef+cetk43VfAIyJ6uf+I94F0AXu3214F4oXqGo33kX4buBmERkjIpG+8xcHsn/nD3jv5RdEJNP9O1FE5rr5+UCdqraLyEnAv0St+ySwREQWu/fmu0SdA0TkChEpc+9Xg5vc7ebtFJGlsQJyd55lA1nudbb7XEXmZ7j5QSDo5kfuFPsfvG6ez7hlbsO7XrOll+OvBA6Lep2Pd22sGsgQkds4tEUSK+YJInKRS0gdeC2Z7n5WG3UsQYxgqroVuAL4Cd636AvxboftxDvRfd9N34f3rfEbbtUn3P+1IvKm+/lKvBPAu0A93slmUpxx1OL1nX8Vr0vj34ElqloT5/rvAj/Cu0OoEpiPd8dTb17E+wa6T0Qi+/g63sXe10SkCXie+Pqpo+3HO5GA941zv4uvQ1X3Rf7hfYPtcj/j3u9P4r2HDcAX8boEOwe4f9z1k3Pwrp/swfvd3ck/E9f1wHdFpBnvZPubqHU3ATfgtTL24v0eK6I2fy6wSURa8FpBl7lEk4XX0ujtmsQM915EvvXvx7u7KuJbbtqteJ/H/W4aqlqNd+vwf7h4Tubga0M93Q1cIiL1IvJj4Dm8BP0eXrdlO/10VzoBvM/jHry7wU7De+9MFDn4xgpjjDmYiCwGbnDdT2YUsQRhjDEmJutiMsYYE5MlCGOMMTFZgjDGGBPTiCp7W1paqjNnzkx1GMYYM2ysW7euRlXLYs0bUQli5syZrF27NtVhGGPMsCEisWqfAdbFZIwxpheWIIwxxsRkCcIYY0xMI+oahDEmPXR1dVFRUUF7+6AK5BofZGdnM3XqVDIzM+NexxKEMSbhKioqyM/PZ+bMmRxcxNakgqpSW1tLRUUFs2bNins962IyxiRce3s7JSUllhzShIhQUlIy4BadJQhjjC8sOaSXwfw+LEEYM0zUf/g+W1/8ZarDMKOIJQhjhgNVGh+5mjkv30htTWWqo0l7O3fu5JhjDn6m1LJly/jhD3/Y53pr167l5ptvBqCjo4OzzjqLBQsW8Pjjj/sW60DcdtttPP/880nbn12kNmYY6H53JTPbNgLwzusvcdr5fT1TxwxWeXk55eXlALz11lt0dXWxfv36uNfv7u4mGAz6Elt3dzff/e53B7zOUOKxFoQx6S7USdezt7E9PImwCvXvvZrqiIa9008/na9//eucdNJJHHnkkbzyyisArF69miVLllBVVcUVV1zB+vXrWbBgAdu3b+eFF17g+OOPZ/78+Xzxi1+ko6MD8Er8fPe732Xx4sU88cQTzJw5k2984xt85CMfoby8nDfffJNPfOITHH744SxfvvyQWHbu3MlRRx3FVVddxbHHHssll1xCW1tbzG0vXbqUJ598EiDueIbCWhDGpLu1D5DdvJP/DH+dO3J/Q1Hd27R1hsjJGh5/vt95ehPv7mlK6DbnTS7g9guPHtI2QqEQr7/+OqtWreI73/nOQV0348eP5/777+eHP/whf/jDH2hvb+f000/nhRde4Mgjj+TKK6/k3nvv5ZZbbgG8MQZr1qwB4NZbb2XatGm8+uqrfPnLX2bp0qX89a9/pb29naOPPpprr732kFi2bt3KAw88wCmnnMIXv/hFfvazn/G1r33tkG0/++yzgHeX2NKlS+OKZyisBWFMOttfj/7lTt4IHEfXrI8TmFrOsfI+L2+tTnVkaa23O3aip3/6058G4IQTTmDnzp19bm/r1q3MmjWLI488EoCrrrqKl19++cD8Sy+99KDlL7roIgDmz5/PySefTH5+PmVlZWRnZ9PQ0HDI9qdNm8Ypp5wCwBVXXHHQyb3ntgcTz2ANj68gxoxWr98P+xu4rePf+fzREykJnkLG9id48+23OHf+pFRHF5ehftMfjJKSEurr6w+aVldXd9AgsTFjxgAQDAYJhUJ9bq+/RzPn5uYe9Dqy7UAgcODnyOtY++qZ0KJf99z2YOIZLGtBGJPOKjfSMHYam3UGZ8+bQMa0EwFo2f4qoe5wioNLX3l5eUyaNIkXXngB8JLDs88+y+LFiwe1vaOOOoqdO3eybds2AH75y19y2mmnJSzeXbt28eqr3rWlRx99tN84/Y4nwhKEMemsfic7QmUsmFbEhIJsGD+XUDCH2V1bWftBff/rj2IPP/wwd9xxBwsWLODMM8/k9ttv5/DDDx/UtrKzs/nFL37BZz/7WebPn08gEIh5LWGw5s6dy0MPPcSxxx5LXV0d1113XUrjiZD+mirDSXl5udoDg8xIEv6/M3iktZzms+7k+tOPAKD7wfN554N9rDzxV3x7ybwURxjb5s2bmTt3bqrDGBZ27tzJkiVLeOedd3zfV6zfi4isU9XyWMv72oIQkXNFZKuIbBORW2PMP0pEXhWRDhH5WtT0aSLykohsFpFNIvIlP+M0Ji3tbyDQ0cAuHc858yYemBycdiJHywes3rSr375oY4bCtwQhIkHgHuA8YB5wuYj0/LpTB9wM9BzeGAK+qqpzgUXADTHWNWZka/CeBNk4ZjJHjM/75/Sp5WQQoqBhC/VtXSkKziTKzJkzk9J6GAw/WxAnAdtUdYeqdgKPARdHL6CqVar6BtDVY/peVX3T/dwMbAam+BirMemn3ksQoYLpB0+f4vUGLAhsY0/D/mRHZUYRPxPEFGB31OsKBnGSF5GZwPHA33uZf42IrBWRtdXVdm+4GUHqdwIgxTMPnl4wic7cyRwf2MbeRnsgj/GPnwki1kiVAXWYikge8BRwi6rGHIqpqveparmqlpeVlQ0iTGPSk9Z/QKPmMq740M91eMoJLJBt7G20FoTxj58JogKYFvV6KrAn3pVFJBMvOTyiqr9NcGzGpL1Q7Q52aRmTi8YeMi9r6kKmB6qpqa1NQWRmtPAzQbwBzBaRWSKSBVwGrIxnRfGGET4AbFbVu3yM0Zi0Fa77gN06nslF2YfMCxRNBaCttiLZYQ0LtbW1LFiwgAULFjBx4kSmTJly4HVnZycAK1eu5Pvf/37C9x1dUC+R9uzZwyWXXJLw7fbFt1IbqhoSkRuB54Ag8KCqbhKRa9385SIyEVgLFABhEbkF746nY4EvABtFZL3b5DdUdZVf8RqTVsJhMpt3sUvn8tEYLQjyvdteQw1xN8pHlZKSkgNlupctW0ZeXt6B4nfgFeq76KKLDtRMSnehUIjJkycPKPGEQiEyMoZ2ive1FpM7oa/qMW151M/78LqeelpD7GsYxowOLfsIhLuo0DImFcZKEF4dJmnZm+TAhq+lS5dSXFzMW2+9xcKFC5k/fz5r167lpz/9KU8//TR33HEHnZ2dlJSU8MgjjzBhwgSWLVvGrl272LFjB7t27eKWW2458ECh733vezzyyCNMmzaN0tJSTjjhhIOSEMC6dev4yle+QktLC6WlpaxYsYJJkyYdEld2djabNm2isrKSu+66iyVLlrBixQqeeeYZ2tvbaW1t5cEHHzwwoK69vZ3rrruOtWvXkpGRwV133cUZZ5xxyDovvvjikN4zK9ZnTDpydzDtCUykJDfr0PmuBZG9v4rusBIMpPH3qT/eCvs2JnabE+fDeQPvHnrvvfd4/vnnCQaDrFix4sD0xYsX89prryEi3H///fzgBz/gRz/6EQBbtmzhpZdeorm5mTlz5nDdddfx9ttv89RTT/HWW28RCoVYuHAhJ5xwwkH76urq4qabbuL3v/89ZWVlPP7443zzm9/kwQcfPCSunTt38pe//IXt27dzxhlnHKix9Oqrr7JhwwaKi4sPqjh7zz33ALBx40a2bNnCOeecw3vvvXfIOkNlCcKYdOTGQHTkTSUQ6+Q/Jp/OYC5loXpqWjq8Ok2mX5/97GdjPmGtoqKCSy+9lL1799LZ2XlQ1dcLLriAMWPGMGbMGMaPH09lZSVr1qzh4osvZuxYr3V34YUXHrLNrVu38s4773D22WcD3tPderYeIj73uc8RCASYPXs2hx12GFu2bAHg7LPPjnmiX7NmDTfddBPgFe6bMWPGgQTR2zqDYQnCmHTU8AFhhEDR9F4X6coZz/jOevY07E/vBDGIb/p+6a0M9k033cRXvvIVLrroIlavXs2yZcsOzIsu1x0pDR5PiRNV5eijjz5QpbUvvZX77i3evvafqFLfYNVcjUlP9Tupppjx4wp7XyZ/EhOk3gbLJUBjYyNTpnjjeB966KF+l1+8eDFPP/007e3ttLS08MwzzxyyzJw5c6iurj6QILq6uti0aVPM7T3xxBOEw2G2b9/Ojh07mDNnTp/7P/XUU3nkkUcAr9ts165d/a4zGNaCMCYNaf1OPgiXxbzFNSKzaDITKt7nbSu3MWTLli3js5/9LFOmTGHRokX84x//6HP5E088kYsuuojjjjuOGTNmUF5eTmHhwck8KyuLJ598kptvvpnGxkZCoRC33HILRx996AOU5syZw2mnnUZlZSXLly8nO7vvFuH111/Ptddey/z588nIyGDFihUHtXQSxcp9G5OGun94FP/TOJuuC+/h8pNidzPpn75N51/v5c4T/sJtFyX/qW19GQ3lvltaWsjLy6OtrY1TTz2V++67j4ULFw54O0uXLmXJkiVJGeMw0HLf1oIwJt10tRNo2cfu8Ec5vrD3b5JSMJkx0kVTfSWQXgliNLjmmmt49913aW9v56qrrhpUckh3liCMSTeNuxGUXTqeJbEGyUW4W107bbBcSvz6179OyHaib7dNN3aR2ph0425x3a1lTOozQbhbJpvSc7DcSOq+HgkG8/uwBGFMuqn3LpDWj5lM3pg+GvmuBTGmvYrOUDgZkcUtOzub2tpaSxJpQlWpra3t9+J3T9bFZEy6afiALskkqyD2oKoD8rwEMZ4GKpvamVack4Tg4jN16lQqKiqwZ7Skj+zsbKZOjVXZqHeWIIxJN017qZZSJo3rZ8BTZjZdWUVMCHmD5dIpQWRmZh40GtkMT9bFZEy6aa2iKlzQ5xiIiHDeRCZKnQ2WM76wBGFMmgm3VFPZnR+7imsPwcLJjJd69tiT5YwPLEEYk2bCLVXUaCFT+rqDyckonMykQAN7G6wFYRLPEoQx6STcTXB/HTUUMKmPQXIH5E+klAb2NbT4H5sZdSxBGJNO2moRlBotjPks6kPkTyRImLb6Kv9jM6OOJQhj0kmrd1toHQVMjKsF4d0KG26y0dQm8SxBGJNOWryWQFd2KZnBOP483ViJsR3V7O/s9jMyMwpZgjAmnbTWeP/nlsW3vGtBTBTvyXLGJJIlCGPSSavXgpD88fEtnzseRZgg9dS2dvoYmBmNLEEYk05aq+kig7F5cT5TOJhBaGwp46mn1loQJsF8TRAicq6IbBWRbSJya4z5R4nIqyLSISJfG8i6xoxIrdXUagEl+fEXVdO8idaCML7wLUGISBC4BzgPmAdcLiLzeixWB9wM/HAQ6xoz4nQ3V1GjBZTkZcW9TqBwMhOknjpLECbB/GxBnARsU9UdqtoJPAZcHL2Aqlap6htA10DXNWYk8hJEIaW58T9fOFgwiQnSYF1MJuH8TBBTgN1RryvctISuKyLXiMhaEVlrpYXNsNdSTS0FFOfG34KQgsmUSiP1LW0+BmZGIz8ThMSYFu/TQ+JeV1XvU9VyVS0vK4vz1kBj0pEqwfaaAXcxkefd8dTdVOlTYGa08jNBVADTol5PBeId7jmUdY0ZnjqaCXZ3eF1MefF3MZFTAkA4MobCmATxM0G8AcwWkVkikgVcBqxMwrrGDE+uzEatDqyLidxSAKS11o+ozCjm2xPlVDUkIjcCzwFB4EFV3SQi17r5y0VkIrAWKADCInILME9Vm2Kt61esxqQFlyCaguPIyQrGv16OlyCCHXV+RGVGMV8fOaqqq4BVPaYtj/p5H173UVzrGjOiuQQRGluKSKzLcL1wXUz53Y20dYbIybInCZvEsJHUxqQLV6hP4q3DFDF2HGECFEsTtS02FsIkjiUIY9KFu8gczB9ggggECI0pophmG01tEsoShDHporWaRvIoys8d8KrhscUUSzN1rTZYziSOJQhj0oS2DrzMxgE5pZRIEzXWxWQSyBKEMWmiu7mKai2kZCC3uDoZ+WWMo9nqMZmEsgRhTJrQSKG+AdRhigjmlVIiliBMYlmCMCZNSFsNNVo4qC4myS2lSFqoa273ITIzWlmCMCYdhDrJ6GykVgsGVmYjIqeEIGH2N9toapM4liCMSQdt3i2uNRQOrMxGhBtNHW6xisYmcSxBGJMO3CC5Addhisj1RlPTZi0IkziWIIxJB26QXFtmCdmZA6jDFOHKbQT3W4IwiWMJwph00Oq1IMLuRD9grospP+zVYzImESxBGJMOXKE+cscPbn2XWMbRbPWYTMJYgjAmHbRU0U4WufmFg1s/M5tQRq6NhTAJZQnCmHTQWkMdhZTmD+IWV6c7u9ir6Gr1mEyCWIIwJg1oaw3V4UHewRSRW+pVdLUuJpMgliCMSQPdLdXUad6gymxEBPNKKZYm62IyCWMJwpg0oG211DHISq5OpB6TPRPCJIolCGPSQGB/7eDLbDiSW0qxWBeTSRxLEMakWmcbwdB+6jV/aNcgckrIppOW5oaEhWZGN0sQxqSaK49RS/6Qupgig+W6W2oSEZUxliCMSTlXqK9OCyjOGVoLAkBbrdyGSQxfE4SInCsiW0Vkm4jcGmO+iMiP3fwNIrIwat6XRWSTiLwjIo+KSLafsRqTMu6E3jVmHBnBIfxJ5notCKvHZBLFtwQhIkHgHuA8YB5wuYjM67HYecBs9+8a4F637hTgZqBcVY8BgsBlfsVqTEq5LqZB12GKcOvndTdYPSaTEH62IE4CtqnqDlXtBB4DLu6xzMXAw+p5DSgSkUluXgYwVkQygBxgj4+xGpM6rospmFs2tO1E6jHZnUwmQfxMEFOA3VGvK9y0fpdR1Q+BHwK7gL1Ao6r+KdZOROQaEVkrImurq+1hKWYYaq0hRJDsvHFD2052IeFAJiXSTH2bJQgzdHElCBFZIiIDTSYSY5rGs4yIjMNrXcwCJgO5InJFrJ2o6n2qWq6q5WVlQ/wGZkwqtNXSQD7FQ6jDBIAIoTHjKKbJBsuZhIj3pH8Z8L6I/EBE5sa5TgUwLer1VA7tJuptmbOAf6hqtap2Ab8FPhrnfo0ZVrS1hppwPiVDGQMRkVNCsTRTbwnCJEBcCUJVrwCOB7YDvxCRV13XTn4fq70BzBaRWSKShZdkVvZYZiVwpbubaRFeV9JevK6lRSKSIyICfBzYPLBDM2Z48OowDXGQnBPI80ZTWz0mkwhxdxupahPwFN7F5knAp4A3ReSmXpYPATcCz+Gd3H+jqptE5FoRudYttgrYAWwDfg5c79b9O/Ak8Caw0cV534CPzphhINzq1WFKRIII5pVRYgX7TIJkxLOQiFwEXA0cDvwSOElVq0QkB+/k/5NY66nqKrwkED1tedTPCtzQy7q3A7fHE58xw5lXh+kwDktAgpCcEoqlxRKESYi4EgRwCfDfqvpy9ERVbRORLyY+LGNGie4QwY5G6klMFxO5pRTSQkNL29C3ZUa9eLuY9vZMDiJyJ4CqvpDwqIwZLfbXIyi1WjCkZ0Ec4MZChJqtHpMZungTxNkxpp2XyECMGZWi6jCNy80c+vYO1GOyBGGGrs8uJhG5Du/C8eEisiFqVj7wVz8DM2ZUcGU22jOLGJMRHPr2XD0maa8b+rbMqNffNYhfA38E/i8QXWyvWVXtE2jMULlv+uGc4sRsz5X8HtNRT6g7PLTif2bU6+/To6q6E+9Oo+aof4hIgj7RxoxirotJ3Il9yFwLoliaqG/rSsw2zagVTwtiCbAOr0xGdGkMBQ7zKS5jRoc2ryGemZ+gBDHW+95WIk3Ut3VSNtTyHWZU6zNBqOoS9/+s5IRjzCjTWkMTuRTm5SRme8EMurIKKQ65iq4TErNZMzr1d5F6YV/zVfXNxIZjzOiibTWuzEbivumHx5ZQvL/JKrqaIeuvi+lHfcxT4MwExmLMqNPdUkOtJqhQnyO5pZTUNbPNRlObIeqvi+mMZAVizGgUbqn2nkWdwASRkV9Gseyxiq5myPrrYjpTVV8UkU/Hmq+qv/UnLGNGibZa6nQu4/MSlyACuaWUBKyiqxm6/rqYTgNeBC6MMU/xntNgjBkMVTLa66gjn7kJbEGQW8o4mqlraU/cNs2o1F8X0+3u/6uTE44xo0hHM4FwF7UJ7mIip5QgYTpabCyrGZp4HzlaIiI/FpE3RWSdiNwtIiV+B2fMiOYGydVrfmIK9UW4wXLdLVaPyQxNvOPwHwOqgc/glf6uBh73KyhjRgU3SK4lWMTYrATUYYpwZTuCrs6TMYMV7/MgilX1e1Gv7xCRT/oQjzGjh6vDFBqb4Ko1rmxHRnstqor31F5jBi7eFsRLInKZiATcv88Bz/gZmDEjnutiCuQmuLfWdTEVaCOtnd2J3bYZVfq7zbWZf9Zg+grwKzcrALRgjwQ1ZvBcF1Agryyx23UtiGKaqWvpJG9MvB0Fxhysv7uY8pMViDGjTmsNHWSRm1eQ2O1mZhPKyKUk1ERtawfTSxJU58mMOnF/tRCRccBsIDsyredjSI0xA9BWS53mU5KX+Iqr3dnFjOtotnpMZkjiShAi8m/Al4CpwHpgEfAqVovJmEHrbqmmVvMZl8gxEBG5pZQ0NlHZYgnCDF68F6m/BJwIfODqMx2Pd6trn0TkXBHZKiLbROTWGPPFja/YJiIboqvHikiRiDwpIltEZLOIfCTOWI0ZFrqbKqnRwoQW6osI5pVSLNaCMEMTb4JoV9V2ABEZo6pbgDl9rSAiQeAe4DxgHnC5iMzrsdh5eN1Ws4FrgHuj5t0NPKuqRwHHAZvjjNWY4aG1mhoKE1rqOyKYV0qJNFFr9ZjMEMSbICpEpAj4HfBnEfk9sKefdU4CtqnqDlXtxBtsd3GPZS4GHlbPa0CRiEwSkQLgVOABAFXtVNWGOGM1Jv2pkrG/mhotTGyZDUdyXQuipSPh2zajR1zXIFT1U+7HZSLyElAIPNvPalOA3VGvK4CT41hmChDC68L6hYgch/fI0y+pamvPnYjINXitD6ZPnx7P4RiTeu2NBMJdVGuBL11M5JQyhi5aWxoTv20zasTbgkBEForIzcCxQIVrFfS5SoxpGucyGcBC4F5VPR5oBQ65hgGgqveparmqlpeVJfh+cmP80updwqvRQooTWOr7gAP1mPq9VGhMr+It1ncb8BBQApTifbP/Vj+rVQDTol5P5dBuqd6WqcBLQn9305/ESxjGjAwtVQDUB8aR78dANjdYTlqtHpMZvHhbEJcDJ6rq7a4E+CLg8/2s8wYwW0RmiUgWcBmwsscyK4Er3d1Mi4BGVd2rqvuA3SISuRD+ceDdOGM1Jv21VALQlV3qT60k14KQdiv5bQYv3q8uO/EGyEWeQDIG2N7XCqoaEpEbgeeAIPCgqm4SkWvd/OXAKuB8YBvQBkQ/d+Im4BGXXHb0mGfM8Oa6mDR3vD/bz/HqO+V01dPVHSYzGHdvsjEH9FeL6Sd41wQ6gE0i8mf3+mxgTX8bV9VVeEkgetryqJ8VuKGXddcD5f3tw5hhqaWKbgJkF5T6s/3cSD2mJupbOxlfkN3PCsYcqr8WxFr3/zrgf6Kmr/YlGmNGi9Yq6imgtGCsP9vPyqM7kEmxNFPTYgnCDE5/xfoeivzsunqOdC+3qmqXn4EZM5JpSxXVWkhZfuIHyQEgQnd2CSWdTdTYWAgzSPHWYjod7y6mnXi3pk4TkausWJ8xg9PdXEV1uIAyHwr1RWhOCcXNzVQ3W4IwgxPvReofAeeo6lYAETkSeBQ4wa/AjBnJtKWKamYyvsC/BJGRV0ZJ1W7etxaEGaR4b23IjCQHAFV9D8j0JyRjRjhVAm1emQ0/WxDB/DJKxFoQZvDibUGsE5EHgF+615/Hu3BtjBmojmaC3R1egvDrGgRAjlewzxKEGax4E8S1eLej3ox3DeJl4Gd+BWXMiOZGUddoob93F+WWkEM7DU3N/u3DjGj9JggRCQDrVPUY4C7/QzJmhGv1EkRTsJjcrKB/+3HlNjqbrR6TGZx+r0Goahh4W0SsVKoxieBaEN25Zf6U2YjI9YpXiktIxgxUvF1Mk/BGUr+OV1kVAFW9yJeojBnJXJmNQJ5PZTYi8icCMLajho5QN2MyfGytmBEp3gTxHV+jMGY0iZTZKPS5PH3eBADGSwO1LZ1MLvJp1LYZsfqrxZSNd4H6CGAj8ICqhpIRmDEjVmsVDeRTWpDj735cgphAPdXNHZYgzID1dw3iIbyCeRvxnh/9I98jMmaE626uosrnUdQAZGQRyi5mvDTYra5mUPrrYpqnqvMB3DiI1/0PyZiRrbup0v8xEI7mTWB8awPVNpraDEJ/LYgDBfmsa8mYxNCWKmoo9LXMRkSwYDLjpd5aEGZQ+mtBHCciTe5nAca614L3OIcCX6MzZqRRJWN/NTU6nyPy/C/BHSiYyMTAW5YgzKD0V+7b7oszJpE6Wwh2tyeti4m8CZTQSE3Tfv/3ZUYcew6hMckUKbNBISV5Wf7vL38iGXTTbqOpzSBYgjAmmdwguY4xpcl5TrQbLCfNe/3flxlxLEEYk0yuBRHO8XmQXESelyAyWqvwHgFvTPwsQRiTTK4uUqDA5zIbEfneYLnCcB2tnd3J2acZMSxBGJNMLdWEEbILk5QgXAtiPDZYzgycrwlCRM4Vka0isk1Ebo0xX0Tkx27+BhFZ2GN+UETeEpE/+BmnMcmiLVU0aB4lBbnJ2WFmNl1ZhUywsRBmEHxLECISBO7BK9ExD7hcROb1WOw8YLb7dw1wb4/5XwI2+xWjMckWaqqk2udHjfYUzp1g5TbMoPjZgjgJ2KaqO1S1E3gMuLjHMhcDD6vnNaBIRCYBiMhU4ALgfh9jNCapuhs/pEqLkjMGwpGCiW40dXvS9mlGBj8TxBRgd9TrCjct3mX+H/DvQLivnYjINSKyVkTWVlfbvd4mvQWaPuRDLWV8vv+jqCMyCycxQRqoaelM2j7NyOBngoj1qKye99nFXEZElgBVqrquv52o6n2qWq6q5WVlSbp10JjBCHWQ1V7NHi1NbgsifyJl0kB1k7UgzMD4mSAqgGlRr6cCe+Jc5hTgIhHZidc1daaI/Mq/UI1JgqYPAdhDSVITBPmTyCJEW5O1sM3A+Jkg3gBmi8gsEckCLgNW9lhmJXClu5tpEdCoqntV9f+o6lRVnenWe1FVr/AxVmP811gBQHWgjILseB/mmADuwUHhpn3J26cZEXz7lKpqSERuBJ4DgsCDqrpJRK5185cDq4DzgW1AG3C1X/EYk3IuQezPmYxIrN5Vn7hyG8HWyuTt04wIvn6NUdVVeEkgetryqJ8VuKGfbawGVvsQnjHJ1eDdj6EFPe/V8JlrQYxpryIcVgKBJCYnM6zZSGpjkqVxN7WMY2JxYXL361oQpdpAw/6ufhY25p8sQRiTJNpYwe5wMVPHjU3ujrNy6crIY7zUU2VjIcwAWIIwJklC9bv5UEuYNi4n6fvuzh3PeKlnT4M9OMjEzxKEMcmgSqCpgj1ayrTiJLcggEDBJMZLAxX1liBM/CxBGJMMbbUEu9vZoyVMTUELIrNoMhOkgd11bUnftxm+LEEYkwyN3h1MH1LK5KLkldmIkLwJXoKotQRh4mcJwphkcGMgOnImMyYjmPz9509iDJ3U1dloahM/SxDGJINLEMFx0/pZ0CfuVteuhp7VbozpnSUIY5KhsYJ2sigqnpia/RdO9f7r3EdTu42FMPGxBGFMEoTrd/FhuISpJUl6klxPxYcBMF0q7UK1iZslCGOSoKtuFx9qKdOSPUguIreM7owcZkolu+vsVlcTH0sQxiSBNFXwoZam5BZXLwCB4lnMkEoq6q0FYeJjCcIYv3W1k9Vewx4tSckguYhAyeEcFqi0wXImbpYgjPGbe1DQPillYkHyx0BESPFhTJEqKmqbUxaDGV4sQRjjNzdIriN3MhnBFP7JFc8iixAddRWpi8EMK5YgjPGbGwMRKErRGIgIdydTZuM/8B7FYkzfLEEY4zeXIHJKpqc2DpcgJnbvpa61M7WxmGHBEoQxPgvV76JSi5hUkuQHBfWUP5nuQJY3FsIuVJs4WIIwxmedtR+4Mt8pusU1IhAgVDiDmXarq4mTJQhjfBao284HOj6lt7hGBEsOZ4bss8FyJi6WIIzxU3sj2W17eS88NSVPkuspo/RwZgaq2F3XmupQzDBgCcIYP1VvBWBHYDqleWNSHAxQPIuxdNBc+2GqIzHDgK8JQkTOFZGtIrJNRG6NMV9E5Mdu/gYRWeimTxORl0Rks4hsEpEv+RmnMb6p2gxAS8ERBAKS4mCA4lkABOp2pDgQMxz4liBEJAjcA5wHzAMuF5F5PRY7D5jt/l0D3Oumh4CvqupcYBFwQ4x1jUl/1VtoJ4vMklmpjsTjbnUd27yLcNjGQpi++dmCOAnYpqo7VLUTeAy4uMcyFwMPq+c1oEhEJqnqXlV9E0BVm4HNwBQfYzXGF92V77ItPIV5U4pSHYqncDphyWAqe6lq7kh1NCbN+ZkgpgC7o15XcOhJvt9lRGQmcDzw91g7EZFrRGStiKytrrbHKZr0Eq7czFadwjGTUzwGIiKYQUfuZGZIJR/U2oVq0zc/E0SsDteebdo+lxGRPOAp4BZVbYq1E1W9T1XLVbW8rKxs0MEak3D7G8hsq+T98FSOmZImCQKvqusMqeTdvTH/pIw5wM8EUQFEF5+ZCvR8IG6vy4hIJl5yeERVf+tjnMb4o3oLALszZzI1VQ8KimHM+COYFahkY0VDqkMxac7PBPEGMFtEZolIFnAZsLLHMiuBK93dTIuARlXdKyICPABsVtW7fIzRGP9UvQtAcMJcvI90mig+jHza+GC3VXU1fcvwa8OqGhKRG4HngCDwoKpuEpFr3fzlwCrgfGAb0AZc7VY/BfgCsFFE1rtp31DVVX7Fa0yidVdupl3HMHH67FSHcjB3JxN122jtCJE7xrfTgBnmfP1kuBP6qh7Tlkf9rMANMdZbQ+zrE8YMG+0fbuJ9ncLR6XIHU8TE+QDMlx28u7eJE2cWpzggk65sJLUxPgnUbkm7C9QAFE6hO38KJwTeY2NFY6qjMWnMEoQxfmirY2xHLTsD05lVkpvqaA4RnH4yJwa3sfFDSxCmd5YgjPGDK7HRVTInPUps9DR9EROpYd/ubamOxKQxSxDG+CDsEsTYqcekOJJeTDsJgJL6t2npCKU4GJOuLEEY44OmXRto0rFMn5FmdzBFTDiG7uBYFsp7vLvHBsyZ2CxBGOODrr3vsk2ncMzUolSHElswk+7JC1kYeI8NNmDO9MIShDGJ1tlGUd3bbGQ2h5el3wXqiKyZizgm8AFbd1emOhSTpixBGJNoO1aTqZ1sH7eYjGAa/4lNO5kMuuna/WaqIzFpKo0/vcYMTx3vrqJZx5J75KmpDqVvU08EYFLzBprbu1IcjElHliCMSaRwmO6tz/KX8LFcsGBGqqPpW04xrQWHsVDeY4MNmDMxWIIwJpH2rieno5oNYxdx9OSCVEfTr6yZH6E88D7PbOhZaNkYSxDGJFTrO88QVqHw2AvSq4JrLzJnLmKcNPPOhnV0hLpTHY5JM5YgjEmg9k3P8KbO5uPlc1MdSnxmfgyAj3X9jZe2VKU4GJNuLEEYkyhNeyhp2sz6sYuYMyE/1dHEp3gW4cPO5MrM5/ndul2pjsakGUsQxiRI04Y/AJA17/xh0b0UETj5GiZQR8b7q6hv7Ux1OCaNWIIwJkEa1j/N7nAZH120ONWhDMzsc+jMm8bn5U/8YePeVEdj0oglCGMSIFzxFlNrXuFvY0/liOHSvRQRCJK56H/xkeC7rHt9TaqjMWnEEoQxQxUOU/vETdRqAVmnfS3V0QyKLPwCocAYyqueYkd1S6rDMWnCEoQxQ1T3txWUNW7kqZJr+ORHhsndSz3lFNM199N8OriG7z35KqHucKojMmnAEoQxQ6Bt9QRfXMZbeiQXfuHLw+ridE9jF19PjnRwyZ4f8JPnt6Q6HJMGLEEYM1iqbHv0f5PX3cTexXcwZVz6Vm6Ny6Rj4RP/yQXB1ylecztr3qtOdUQmxSxBGDMIodZ63r/nM8ze/QTP5X2K8846J9UhJcZHbqDr5Bu4Kvgn3nrsdrseMcr5miBE5FwR2Soi20Tk1hjzRUR+7OZvEJGF8a5rTEqoUrFhNTU/WsSs6pf4Xdm1nHL98mHdtdRT5ifuoGn2J7kp/Ajrf3IZP3/yaZqs2uuolOHXhkUkCNwDnA1UAG+IyEpVfTdqsfOA2e7fycC9wMlxrmuMrzQcprm+hrrKnTTt+wcd773IlKq/MDW8l72U8Nppv+KTZ16Q6jATLxCg4NKf0/bHiSx5cwVZ77zCmneOp27SqeRNO5opsxcwYdI08nPGEgyMnMRoDuVbggBOArap6g4AEXkMuBiIPslfDDysqgq8JiJFIjIJmBnHugnz/vdOIFM7/Ni0SbHeT1/q5iui6v1PmICGyaSLMdrBGDopkG4iNVk7NJN3sxewbfpS5p69lMXjJybhCFIkI4ucC++Ej3+dyhfu4Zj1v6Bo792wF3jdW2S/ZtEmY+kik26CqAQII4Cg4v0fiybrGEaRtmAh877514Rv188EMQXYHfW6Aq+V0N8yU+JcFwARuQa4BmD69OmDCrQxdyaBsJUYGLl6O1G56SIoAZAAKgE0mIVmZBMOjkXySskomsrYkqlMn3cyxxcUJS/sdJBTzIQLvw1LvoW2VFK1/W2q/rGBruZawu3N0NEM4S4IdyMaAg2DS7ixWXrwQyjTn9LyfiaIWH+VPT8dvS0Tz7reRNX7gPsAysvLB/XpK//KU4NZzZjRQwTJn8iEBROZsOATqY7GJImfCaICmBb1eirQ86kkvS2TFce6xhhjfOTnXUxvALNFZJaIZAGXASt7LLMSuNLdzbQIaFTVvXGua4wxxke+tSBUNSQiNwLPAUHgQVXdJCLXuvnLgVXA+cA2oA24uq91/YrVGGPMocS7gWhkKC8v17Vr16Y6DGOMGTZEZJ2qlseaZyOpjTHGxGQJwhhjTEyWIIwxxsRkCcIYY0xMI+oitYhUAx8McvVSoCaB4QwHdswj32g7XrBjHqgZqloWa8aIShBDISJre7uSP1LZMY98o+14wY45kayLyRhjTEyWIIwxxsRkCeKf7kt1AClgxzzyjbbjBTvmhLFrEMYYY2KyFoQxxpiYLEEYY4yJaVQlCBE5V0S2isg2Ebk1xnwRkR+7+RtEZGEq4kykOI758+5YN4jI30TkuFTEmUj9HXPUcieKSLeIXJLM+PwQzzGLyOkisl5ENonIX5IdY6LF8dkuFJGnReRtd8xXpyLORBGRB0WkSkTe6WV+4s9fqjoq/uGVDd8OHIb3QKK3gXk9ljkf+CPeE+0WAX9PddxJOOaPAuPcz+eNhmOOWu5FvJLzl6Q67iT8novwnuk+3b0en+q4k3DM3wDudD+XAXVAVqpjH8IxnwosBN7pZX7Cz1+jqQVxErBNVXeoaifwGHBxj2UuBh5Wz2tAkYhMSnagCdTvMavq31S13r18De/pfcNZPL9ngJuAp4CqZAbnk3iO+V+A36rqLgBVHe7HHc8xK5AvIgLk4SWIUHLDTBxVfRnvGHqT8PPXaEoQU4DdUa8r3LSBLjOcDPR4/hXvG8hw1u8xi8gU4FPA8iTG5ad4fs9HAuNEZLWIrBORK5MWnT/iOeafAnPxHle8EfiSqoaTE15KJPz85eczqdONxJjW8x7feJYZTuI+HhE5Ay9BLPY1Iv/Fc8z/D/i6qnZ7Xy6HvXiOOQM4Afg4MBZ4VUReU9X3/A7OJ/Ec8yeA9cCZwOHAn0XkFVVt8jm2VEn4+Ws0JYgKYFrU66l43ywGusxwEtfxiMixwP3Aeapam6TY/BLPMZcDj7nkUAqcLyIhVf1dUiJMvHg/2zWq2gq0isjLwHHAcE0Q8Rzz1cD31eug3yYi/wCOAl5PTohJl/Dz12jqYnoDmC0is0QkC7gMWNljmZXAle5ugEVAo6ruTXagCdTvMYvIdOC3wBeG8bfJaP0es6rOUtWZqjoTeBK4fhgnB4jvs/174GMikiEiOcDJwOYkx5lI8RzzLrwWEyIyAZgD7EhqlMmV8PPXqGlBqGpIRG4EnsO7A+JBVd0kIte6+cvx7mg5H9gGtOF9Axm24jzm24AS4GfuG3VIh3ElzDiPeUSJ55hVdbOIPAtsAMLA/aoa83bJ4SDO3/P3gBUishGv++Xrqjpsy4CLyKPA6UCpiFQAtwOZ4N/5y0ptGGOMiWk0dTEZY4wZAEsQxhhjYrIEYYwxJiZLEMYYY2KyBGGMMSYmSxBmWHBVV9eLyDsi8oS7lz8R210lIkUiMrOPKpmrRcS3W39FZEUqKsqKyGdFZLOIvJTsfZvhwRKEGS72q+oCVT0G6ASuTcRGVfV8VW1IxLaGoX/FGyR4RqoDMenJEoQZjl4BjhCRYhH5nat9/5orGYKInOZaG+tF5C0RyReRSSLyclQr5GNu2Z0iUuq2myEiD7ntPRmrlSIi54jIqyLypmvJ5PWYP1dEXo96PVNENrifbxORN9z+75MYhaCi4xGRchFZ7X7OFe95AG+4Y7rYTT9aRF53x7VBRGbH2OblIrLR7ffOSCx4dbeWi8h/9Vj+UyLyvBuRO0lE3hORiXH+bswIYgnCDCsikoH33IqNwHeAt1T1WLza/w+7xb4G3KCqC4CPAfvxyl0/56Ydh1fErac5wH1ue03A9T32XQp8CzhLVRcCa4GvRC+jqpuBLBE5zE26FPiN+/mnqnqiawWNBZYM4NC/CbyoqicCZwD/JSK5eC2pu91xlePV44mOeTJwJ17BugXAiSLySVX9rov/86r6v3scw/8A+4AbgJ8Dt6vqvgHEakYISxBmuBgrIuvxTmq7gAfwvgH/EkBVXwRKRKQQ+Ctwl4jcDBSpagivds/VIrIMmK+qzTH2sVtV/+p+/hWHVrZdBMwD/upiuQqYEWM7vwE+536+FHjc/XyGiPzdlX44Ezg6/sPnHOBWt9/VQDYwHXgV+IaIfB2Yoar7e6x3IrBaVavd+/AI3oNn+nMT8H+ADlV9dABxmhFk1NRiMsPefvct+YBYXTSAqur3ReQZvLo0r4nIWar6soicClwA/FJE/ktVH+65bj+vBfizql7eT6yPA0+IyG9dPO+LSDbwM6BcVXe7RJUdY90Q//ziFj1fgM+o6tYey28Wkb+743pORP7NJcvo9QZjCl7NpgkiEhjhz1EwvbAWhBnOXgY+D97zlvHKWTeJyOGqulFV78RrcRwlIjOAKlX9OV7rI9bzeqeLyEfcz5cDa3rMfw04RUSOcPvMEZEje25EVbcD3cC3+WfrIXKyr3HXLXq7a2kn3nMbAD4TNf054KZIUhSR493/hwE7VPXHeNU8j+2xvb8Dp4lIqYgE3XH1+Txq1433C7xuuc306EYzo4clCDOcLQPK3UXg7+N1+QDc4i7Ivo13/eGPeFUw14vIW3gn3rtjbG8zcJXbXjFwb/RMVa0GlgKPumVew3u+QCyPA1fgrj+4O6V+jnft5Hd4XV6xfAe4W0RewUsyEd/Dq9y5wd2O+z03/VLgHdf1dBT/vA4TiXkvXlfRS3jPbX5TVX/fy74jvgG8oqqv4CWHfxORuf2sY0Ygq+ZqjDEmJmtBGGOMickShDHGmJgsQRhjjInJEoQxxpiYLEEYY4yJyRKEMcaYmCxBGGOMien/B5uTr8oDeqieAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior1.plot(label='Uniform prior')\n", "posterior2.plot(label='Triangle prior')\n", "\n", "plt.xlabel('Possible values of x')\n", "plt.ylabel('Probability')\n", "plt.title('Posterior after 140 heads, 110 tails')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The posterior distributions are almost identical because, in this case, we have enough data to \"swamp the prior\"; that is, the posteriors depend strongly on the data and only weakly on the priors.\n", "\n", "This is good news, because it suggests that we can use data to resolve arguments. Suppose two people disagree about the correct prior. If neither can persuade the other, they might have to agree to disagree.\n", "\n", "But if they get new data, and each of them does a Bayesian update, they will usually find their beliefs converging.\n", "\n", "And with enough data, the remaining difference can be so small that it makes no difference in practice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summarizing the posterior distribution\n", "\n", "The posterior distribution contains all of the information we have about the value of $x$. But sometimes we want to summarize this information.\n", "\n", "We have already seen one way to summarize a posterior distribution, the Maximum Aposteori Probability, or MAP:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.56" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "posterior1.idxmax()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`idxmax` returns the value of $x$ with the highest probability.\n", "\n", "In this example, we get the same MAP with the triangle prior:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.56" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "posterior2.idxmax()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to summarize the posterior distribution is the posterior mean.\n", "\n", "Given a set of values, $x_i$, and the corresponding probabilities, $p_i$, the mean of the distribution is:\n", "\n", "$\\sum_i x_i p_i$\n", "\n", "The following function takes a Pmf and computes its mean. Note that this function only works correctly if the Pmf is normalized." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def pmf_mean(pmf):\n", " \"\"\"Compute the mean of a PMF.\n", " \n", " pmf: Series representing a PMF\n", " \n", " return: float\n", " \"\"\"\n", " return np.sum(pmf.index * pmf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's the posterior mean based on the uniform prior:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5595238095238096" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pmf_mean(posterior1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's the posterior mean with the triangle prior:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5574349943859507" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pmf_mean(posterior2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The posterior means are not identical, but they are close enough that the difference probably doesn't matter.\n", "\n", "In this example, the posterior mean is very close to the MAP. That's true when the posterior distribution is symmetric, but it is not always true." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If someone asks what we think $x$ is, the MAP or the posterior mean might be a good answer." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Credible intervals\n", "\n", "Another way to summarize a posterior distribution is a credible interval, which is a range of quantities whose probabilities add up to a given total.\n", "\n", "The following function takes a `Series` as a parameter and a probability, `prob`, and return an interval that contains the given probability.\n", "\n", "If you are interested, it computes the cumulative distribution function (CDF) and then uses interpolation to estimate percentiles." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "from scipy.interpolate import interp1d\n", "\n", "def credible_interval(pmf, prob):\n", " \"\"\"Compute the mean of a PMF.\n", " \n", " pmf: Series representing a PMF\n", " prob: probability of the interval\n", " \n", " return: pair of float\n", " \"\"\"\n", " # make the CDF\n", " xs = pmf.index\n", " ys = pmf.cumsum()\n", " \n", " # compute the probabilities\n", " p = (1-prob)/2\n", " ps = [p, 1-p]\n", " \n", " # interpolate the inverse CDF\n", " options = dict(bounds_error=False,\n", " fill_value=(xs[0], xs[-1]), \n", " assume_sorted=True)\n", " interp = interp1d(ys, xs, **options)\n", " return interp(ps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's the 90% credible interval for `posterior1`." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.50259405, 0.60603433])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "credible_interval(posterior1, 0.9)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And for `posterior2`." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.50114608, 0.60373672])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "credible_interval(posterior2, 0.9)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The credible interval for `posterior2` is slightly narrower." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }