{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bayes Theorem Problems\n", "\n", "This notebook presents code and exercises from Think Bayes, second edition.\n", "\n", "Copyright 2016 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "% matplotlib inline\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "import numpy as np\n", "\n", "from thinkbayes2 import Hist, Pmf, Cdf, Suite, Beta\n", "import thinkplot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The sock problem\n", "\n", "Yuzhong Huang\n", "\n", "There are two drawers of socks. The first drawer has 40 white socks and 10 black socks; the second drawer has 20 white socks and 30 black socks. We randomly get 2 socks from a drawer, and it turns out to be a pair(same color) but we don't know the color of these socks. What is the chance that we picked the first drawer.\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "drawer 1 0.566666666667\n", "drawer 2 0.433333333333\n" ] } ], "source": [ "# Solution\n", "\n", "pmf = Pmf(['drawer 1', 'drawer 2'])\n", "pmf['drawer 1'] *= (40/50)**2 + (10/50)**2\n", "pmf['drawer 2'] *= (30/50)**2 + (20/50)**2\n", "pmf.Normalize()\n", "pmf.Print()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "drawer 1 0.568965517241\n", "drawer 2 0.431034482759\n" ] } ], "source": [ "# Solution\n", "\n", "pmf = Pmf(['drawer 1', 'drawer 2'])\n", "pmf['drawer 1'] *= (40/50)*(39/49) + (10/50)*(9/49)\n", "pmf['drawer 2'] *= (30/50)*(29/49) + (20/50)*(19/49)\n", "pmf.Normalize()\n", "pmf.Print()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "class Socks(Suite):\n", " \n", " def Likelihood(self, data, hypo):\n", " \"\"\"Probability of data under hypo.\n", " \n", " data: 'pair' or 'no pair'\n", " hypo: tuple, number of (white, black) socks\n", " \"\"\"\n", " white, black = hypo\n", " total = white + black\n", " like = white/total*(white-1)/(total-1) + black/total*(black-1)/(total-1)\n", " if data == 'pair':\n", " return like\n", " else:\n", " return 1-like " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEphJREFUeJzt3H+sX3V9x/HnqzQsm2SdG4JJS0FpFWxUaBZWh9G7mUkL\nCXVxYTRmimaxGRDJIONXZqgmi7JlmSBBNAKCkVVmolTDoBK4Lo5ZSKDStaU/plRoaHURXMBs/Hrv\nj+8p/fKl/dzvt7fcS3ufj+Sk3/M57/P5fs4nt/d1zznf801VIUnS/sya7gFIkl7fDApJUpNBIUlq\nMigkSU0GhSSpyaCQJDUNFRRJliZ5NMnWJJftp+baJNuSrE9yStc2L8m9STYm2ZDkU331f5bkP5O8\nmGTxQF9XdH1tTvLByRygJGlyJgyKJLOA64AzgEXAiiQnDdQsA06sqoXASuCGbtMLwMVVtQh4D3BB\n374bgD8FfjDQ18nAOcDJwDLg+iQ5sMOTJE3WMGcUpwHbqmpHVT0PrAaWD9QsB24FqKp1wJwkx1bV\nrqpa37U/A2wG5nbrW6pqGzAYAsuB1VX1QlU9BmzrxiBJmgbDBMVc4PG+9Se6tlbNzsGaJCcApwDr\nRny/V/UlSZo6U3IzO8lRwLeAi7ozC0nSIWL2EDU7gfl96/O6tsGa4/ZVk2Q2vZD4elXdMeT77bOv\nfkn8kipJOgBVNdJ932HOKB4EFiQ5PsmRwLnAmoGaNcBHAZIsAZ6uqt3dtpuATVV1TeM9+ge9Bjg3\nyZFJ3gIsAB7Y105V5VLFVVddNe1jeL0szoVz4Vy0lwMx4RlFVb2Y5EJgLb1gubGqNidZ2dtcX6mq\nO5OcmWQ78CxwHkCS04GPABuSPAwUcGVV3ZXkQ8AXgaOB7yVZX1XLqmpTktuBTcDzwPl1oEcnSZq0\nYS49UVV3AW8faPvywPqF+9jv34Ej9tPnd4Dv7Gfb54DPDTM2SdJryyezDwNjY2PTPYTXDediL+di\nL+dicnKoXtVJ4hUpSRpREuo1uJktSZrBDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRk\nUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaF\nJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiS\nmgwKSVLTUEGRZGmSR5NsTXLZfmquTbItyfokp3Rt85Lcm2Rjkg1JPtVX/8Yka5NsSXJ3kjld+/FJ\nfp3koW65/mAcqCTpwEwYFElmAdcBZwCLgBVJThqoWQacWFULgZXADd2mF4CLq2oR8B7ggr59Lwfu\nqaq3A/cCV/R1ub2qFnfL+Qd+eJKkyRrmjOI0YFtV7aiq54HVwPKBmuXArQBVtQ6Yk+TYqtpVVeu7\n9meAzcDcvn1u6V7fAnyor78cyMFIkg6+YYJiLvB43/oT7P1lv7+anYM1SU4ATgF+1DUdU1W7Aapq\nF3BMX/kJ3WWn+5K8d4gxSpJeI7On4k2SHAV8C7ioqp7dT1l1/z4JzK+qp5IsBr6T5B3dGckrrFq1\n6uXXY2NjjI2NHdRxS9Khbnx8nPHx8Un1kapqFyRLgFVVtbRbvxyoqrq6r+YG4L6q+ma3/ijw/qra\nnWQ28D3gX6vqmr59NgNjXc2bu/1P3sf73wdcUlUPDbTXRGOXJL1SEqpqpMv7w1x6ehBY0H0a6Ujg\nXGDNQM0a4KPdIJYAT++5rATcBGzqD4m+fc7rXn8MuKPb/+juBjpJ3gosAH4yykFJkg6eCS89VdWL\nSS4E1tILlhuranOSlb3N9ZWqujPJmUm2A8/SBUCS04GPABuSPEzv8tKVVXUXcDVwe5JPADuAc7q3\nfB/w2STPAS8BK6vq6YN4zJKkEUx46en1yktPkjS61+rSkyRpBjMoJElNBoUkqcmgkCQ1GRSSpCaD\nQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigk\nSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLU\nZFBIkpoMCklSk0EhSWoyKCRJTQaFJKlpqKBIsjTJo0m2JrlsPzXXJtmWZH2SU7q2eUnuTbIxyYYk\nn+qrf2OStUm2JLk7yZy+bVd0fW1O8sHJHqQk6cBNGBRJZgHXAWcAi4AVSU4aqFkGnFhVC4GVwA3d\npheAi6tqEfAe4IK+fS8H7qmqtwP3Ald0fb0DOAc4GVgGXJ8kkzpKSdIBG+aM4jRgW1XtqKrngdXA\n8oGa5cCtAFW1DpiT5Niq2lVV67v2Z4DNwNy+fW7pXt8CfKh7fTawuqpeqKrHgG3dGCRJ02CYoJgL\nPN63/gR7f9nvr2bnYE2SE4BTgB91TcdU1W6AqtoFHDNsX5KkqTMlN7OTHAV8C7ioqp7dT1lNxVgk\nSaOZPUTNTmB+3/q8rm2w5rh91SSZTS8kvl5Vd/TV7O4uT+1O8mbg5xP1NWjVqlUvvx4bG2NsbGyI\nw5GkmWN8fJzx8fFJ9ZGq9h/ySY4AtgAfAJ4EHgBWVNXmvpozgQuq6qwkS4AvVNWSbtutwH9X1cUD\n/V4N/LKqru4+SfXGqrq8u5n9DeAP6F1y+j6wsAYGmmSwSZI0gSRU1UgfEJrwjKKqXkxyIbCW3qWq\nG6tqc5KVvc31laq6M8mZSbYDzwLndQM6HfgIsCHJw/QuL11ZVXcBVwO3J/kEsIPeJ52oqk1Jbgc2\nAc8D55sIkjR9JjyjeL3yjEKSRncgZxQ+mS1JajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoM\nCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQ\nJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUNHu6BzAZx3/g\niukegiQd9jyjkCQ1GRSSpKZU1XSP4YAkqUN17JI0XZJQVRllH88oJElNBoUkqWmooEiyNMmjSbYm\nuWw/Ndcm2ZZkfZJT+9pvTLI7ySMD9e9Kcn+SHye5I8lRXfvxSX6d5KFuuX4yByhJmpwJgyLJLOA6\n4AxgEbAiyUkDNcuAE6tqIbAS+FLf5pu7fQd9Fbi0qt4NfBu4tG/b9qpa3C3nj3JAkqSDa5gzitOA\nbVW1o6qeB1YDywdqlgO3AlTVOmBOkmO79R8CT+2j34XdNoB7gA/3bRvpRosk6bUzTFDMBR7vW3+i\na2vV7NxHzaCNSc7uXp8DzOvbdkJ32em+JO8dYoySpNfIdN7M/gRwQZIHgTcAz3XtTwLzq2oxcAlw\n2577F5KkqTfMV3jsBOb3rc/r2gZrjpug5hWqaivdvYskC4Gzuvbn6EKjqh5K8l/A24CHBvtYtWrV\ny6/HxsYYGxsb4nAkaeYYHx9nfHx8Un1M+MBdkiOALcAH6P21/wCwoqo299WcCVxQVWclWQJ8oaqW\n9G0/AfhuVb2zr+1NVfWL7mb5zcB9VfW1JEcDv6yql5K8FfgB8M6qenpgXD5wJ0kjek0euKuqF4EL\ngbXARmB1VW1OsjLJJ7uaO4GfJtkOfBl4+ZNKSW4D7gfeluRnST7ebVqRZAuwCdhZVV/r2t8HPJLk\nIeB2YOVgSEiSpo5f4SFJM4hf4SFJOugMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQ\nSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUk\nqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKa\nDApJUpNBIUlqMigkSU1DBUWSpUkeTbI1yWX7qbk2ybYk65Oc2td+Y5LdSR4ZqH9XkvuT/DjJHUmO\n6tt2RdfX5iQfPNCDkyRN3oRBkWQWcB1wBrAIWJHkpIGaZcCJVbUQWAl8qW/zzd2+g74KXFpV7wa+\nDVza9fUO4BzgZGAZcH2SjHhckqSDZJgzitOAbVW1o6qeB1YDywdqlgO3AlTVOmBOkmO79R8CT+2j\n34XdNoB7gA93r88GVlfVC1X1GLCtG4MkaRoMExRzgcf71p/o2lo1O/dRM2hjkrO71+cA8ybRlyTp\nNTJ7Gt/7E8AXk3waWAM8N2oHq1atevn12NgYY2NjB2tsknRYGB8fZ3x8fFJ9DBMUO4H5fevzurbB\nmuMmqHmFqtpKd+8iyULgrFH76g8KSdKrDf4R/ZnPfGbkPoa59PQgsCDJ8UmOBM6ldwbQbw3wUYAk\nS4Cnq2p33/Z0y96G5E3dv7OAvwVu6Ovr3CRHJnkLsAB4YKSjkiQdNBMGRVW9CFwIrAU20rvRvDnJ\nyiSf7GruBH6aZDvwZeD8PfsnuQ24H3hbkp8l+Xi3aUWSLcAmYGdVfa3raxNwe9d+J3B+VdVBOVpJ\n0shyqP4OTmJ+SNKIklBVIz1y4JPZkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZ\nFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0Eh\nSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKk\nJoNCktRkUEiSmoYKiiRLkzyaZGuSy/ZTc22SbUnWJzm1r/3GJLuTPDJQ/+4k/5Hk4SQPJPn9rv34\nJL9O8lC3XD+ZA5QkTc6EQZFkFnAdcAawCFiR5KSBmmXAiVW1EFgJfKlv883dvoP+Hriqqk4FrgL+\noW/b9qpa3C3nj3JAM9H4+Ph0D+F1w7nYy7nYy7mYnGHOKE4DtlXVjqp6HlgNLB+oWQ7cClBV64A5\nSY7t1n8IPLWPfl8C5nSvfwfY2bctQx+B/E/Qx7nYy7nYy7mYnNlD1MwFHu9bf4JeeLRqdnZtuxv9\n/jVwd5J/pBcMf9i37YQkDwG/Aj7dhY0kaRpM583svwIuqqr59ELjpq79SWB+VS0GLgFuS3LUNI1R\nklRVzQVYAtzVt345cNlAzQ3An/etPwoc27d+PPDIwD5PD6z/aj/vfx+weB/t5eLi4uIy+jLR7/3B\nZZhLTw8CC5IcT++v/XOBFQM1a4ALgG8mWUIvBPovO4VX33fYmeT9VfWDJB8AtgIkORr4ZVW9lOSt\nwALgJ4ODqirvY0jSFJgwKKrqxSQXAmvpXaq6sao2J1nZ21xfqao7k5yZZDvwLPDxPfsnuQ0YA34v\nyc/ofdLpZuCTwDVJjgD+t1sHeB/w2STP0bvhvbKqnj5YByxJGk26yziSJO3TIflk9jAPAB6u9vUA\nY5I3JlmbZEuSu5PMafVxuEgyL8m9STYm2ZDkU137jJuPJL+RZF33AOuGJFd17TNuLqD3/Ff3wO6a\nbn1GzgNAkseS/HjPw81d20jzccgFxTAPAB7mbubVDzBeDtxTVW8H7gWumPJRTY8XgIurahHwHuCC\n7mdhxs1HVf0f8EfdA6ynAMuSnMYMnIvORcCmvvWZOg/Qu4Q/VlWnVtWeRxtGmo9DLigY7gHAw9Z+\nHmBcDtzSvb4F+NCUDmqaVNWuqlrfvX4G2AzMY+bOx6+7l79B7/5jMQPnIsk84Ezgq33NM24e+oRX\n/64faT4OxaDY1wOAc6dpLK8Xx+z5lFlV7QKOmebxTLkkJ9D7S/pH9D6aPePmo7vc8jCwC/h+VT3I\nzJyLfwL+hl5Q7jET52GPAr6f5MEkf9m1jTQfw3w8VoeeGfUJhe6BzG/Re4DzmSSDxz8j5qOqXgJO\nTfLbwLeTLOLVx35Yz0WSs4DdVbU+yVij9LCehwGnV9WTSd4ErE2yhRF/Lg7FM4qdwPy+9Xm88nui\nZqLde75bK8mbgZ9P83imTJLZ9ELi61V1R9c8Y+cDoKr+BxgHljLz5uJ04OwkPwH+GfjjJF8Hds2w\neXhZVT3Z/fsL4Dv0Lt+P9HNxKAbFyw8AJjmS3gOAa6Z5TFNt8AHGNcB53euPAXcM7nAYuwnYVFXX\n9LXNuPlIcvSeT64k+U3gT+jds5lRc1FVV1bV/Kp6K73fDfdW1V8A32UGzcMeSX5rz1cgJXkD8EFg\nAyP+XBySz1EkWQpcw94HAD8/zUOaMv0PMNL70sWr6P2V8C/AccAO4JyZ8JBiktOBf6P3g7/n6wmu\nBB4AbmcGzUeSd9K7KTmrW75ZVX+X5HeZYXOxR5L3A5dU1dkzdR6SvAX4Nr3/G7OBb1TV50edj0My\nKCRJU+dQvPQkSZpCBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWr6f+lYk2M3ReBAAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "n = 50\n", "socks = Socks()\n", "for white in range(n+1):\n", " socks[white, n-white] = 1\n", "socks.Normalize()\n", "thinkplot.Pdf(socks)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX+x/H3dxICBqRLkaqgoFhAVkBhNYgF8KeAKFUF\nFEFXRFFZEAusFdYGLCqiKE0Fu7CiFDG6NkAEQekoRYRQAoYAMcKc3x+5mRQSMglJJsl8Xs8zj3PO\nnDvznWvIJ7ece805h4iIiC/UBYiISNGgQBAREUCBICIiHgWCiIgACgQREfEoEEREBAgyEMysvZmt\nNbP1ZjYsmzHjzWyDma0ws6ZeX2kzW2xmy81slZmNTDe+kpnNN7N1ZjbPzCrkz1cSEZG8yDEQzMwH\nTACuApoAPc2scaYxHYAGzrkzgIHARADn3J9AW+dcM6Ap0MHMWniLDQcWOucaAYuAB/LnK4mISF4E\ns4XQAtjgnNvinPsLmAl0yjSmEzANwDm3GKhgZtW99iFvTGkgEnDplpnqPZ8KdM7rlxARkRMXTCDU\nArala//m9R1vzPbUMWbmM7PlwE5ggXNuqTemmnMuDsA5txOolvvyRUQkvxT4QWXnnN/bZVQbaGlm\nZ2c3tKBrERGR7EUGMWY7UDddu7bXl3lMneONcc4lmNnnQHtgNRBnZtWdc3FmVgPYldWHm5mCQkQk\nD5xzlpvxwWwhLAUamlk9M4sCegCzM42ZDdwMYGatgP3eL/qqqWcPmdlJwBXA2nTL9PWe9wE+yq4A\n55wezjFy5MiQ11BUHloXWhdaF9k/Dh/+M7gEyCTHLQTn3FEzGwTMJyVAJjvn1pjZwJSX3STn3Fwz\n62hmG4GDQD9v8ZrAVO9MJR8wyzk313ttDPC2md0CbAG6ZVeD3+/H59OUCRGRnGzcEkfXeybladlg\ndhnhnPsUaJSp7+VM7UFZLLcKuCCb94wHLg/m84c//wH/vq9rMENFRMJWUlIy/R6axv6EQzkPzkKx\n+LN71tzveX/B8lCXEXIxMTGhLqHI0LpIo3WRJtzXxeCn3mbr7/F5Xt6cK9rHbM3M1b1sONEnRfHf\nlwbRoM4poS5JRKTIee39r/nXC/8NtLcuGo0rgIPKRcKhw8nc8uBUkpOPhLoUEZEi5ce123jqlU8D\n7bNOr5mn9yk2gQCwefte7hk9K9RliIgUGYmHkrht5IzAH8vly5Xh9SdvztN7FYtA6Na+eeD5x1/8\nxJQPvw1hNSIiRceAkW8QtycBAJ/PGDu8GzVPqZin9yoWgTDmvutodFr1QPuJl+eyan3muXEiIuHl\n+akL+fqHjYF2/65taHfRWXl+v2IRCD6fj6lP9eXkcmUASE4+Qv9HppF4KCnElYmIhMbXyzfynzc+\nD7SbN6nLAwPan9B7FotAAKh5SkWeH3YDPl/KQfOduxO4fdQbIa5KRKTw7dmXyKDH3uLoUT8AVSuV\n47XH+5zwBN5iEwgAV1x8Nrd0bR1o/2/ZRp55fX4IKxIRKVx+v58+D7xO/B8pk88iIyN4aWRvKpaP\nPuH3LlaBAPDggA5ccHbatfZeeDOWBd+sDl1BIiKFaOgz7/HTht8D7XtvbkeLc+vny3sXu0Dw+Xy8\n/kQfqlU+GQC/33HPU2+z5fe9Ia5MRKRgTZ/9He/O+yHQvuLis7izd9t8e/9iFwgAFctH8+pjN1E6\nKuVSTImH/uSmYa+TlJQc4spERArG8jVbefSljwPt02pX5cWHe+XrZxTLQAA4v3EdHr3r2kB7y+97\nueOxN0NYkYhIwYjfn0j/h6cHJp+dXLY000f3IyoqqOuTBq3YBgJAj44X0qPjhYH2ou/W8fzUhSGs\nSEQkf6UcRJ7Cnn2JAPgifPznwZ7UqVk53z+rWAcCwFNDOnNeo9qB9vg3PufzxetCWJGISP554PkP\nWZluIu7dN7albctGx1ki74p9IPh8PqaP7kfVSuUA8B/1c9cTb7FtR94vASsiUhS8+fESZs5dGmhf\n1qoR99wc1G1k8qTYBwKkHWRO3Z924OCf3DjsNR1kFpFia9nqrYycMCfQrndqFV7K54PImZWIQABo\ndlZdHrnj6kB78/a99HtoGn6/P4RViYjk3s49CfR/KO1y/+WiSzN9TD/KlIkq0M8tMYEAcNO1rTIc\nZP5m+SZGpbthhIhIUZeUlEzvoZMDM5EjInyMH9GdeqdWKfDPLlGBACkHmS88t16gPfXDb3nz4yUh\nrEhEJHgD//UGG7fuCrT/eetVJ3QF09wocYHg8/mY8mRfatdIux74yP/MZsmqzaErSkQkCE9MnEvs\nkvWBdufLm3J790sK7fNLXCAAlIsuw5tP9+fksqUBSP7rKLc9PI3tcftCXJmISNbenb+MV979X6B9\nfuM6PD/shkKtoUQGAqQckZ84sjeRkREA7D9wmF5DJ+vMIxEpcpav2cqI5z/EuZT2qdUqMGNMvxO+\nnHVuldhAAGjT/AweHNgh0N68fS+3PDxdZx6JSJGxc08Ctzw4lT+9M4rKRkcxfXQ/ypc7qdBrKdGB\nAHDLda3p3vFvgfbXP2zk4f/MDmFFIiIpkpKS6TX01QxnFP1nRA8a1quew5IFo8QHAsDoIV0ynHk0\nY/ZiJr3zZQgrEpFw5/f7uXnEFDZt3R3oG3rLlYV2RlFWwiIQUs88qntq2sWgnnplHnO//CmEVYlI\nOLvv3++y+MdfA+2uVzbjjh6XhrCiMAkESDnz6O3nbqNyhZTbzPmP+hkyehbL12wNcWUiEm6en7qQ\n9xcsD7Qvano6zwy9PoQVpQgqEMysvZmtNbP1ZjYsmzHjzWyDma0ws6ZeX20zW2RmP5vZKjMbnG78\nSDP7zcx+8B7t8+crZa/mKRWZ8mRfTipTCoCkP4/Qd8RUXQhPRArNO/OWMW76okC7Yd1qTHmiT6Gf\nUZSVHCswMx8wAbgKaAL0NLPGmcZ0ABo4584ABgITvZeOAPc655oAFwF3Zlr2OefcBd7j0xP/Ojk7\nv3Edxo/oTkREylffn3CIHve/QuKhpML4eBEJY18v38iI5z/AeeeXVqt8MrOeu63Ar1EUrGAiqQWw\nwTm3xTn3FzAT6JRpTCdgGoBzbjFQwcyqO+d2OudWeP2JwBqgVrrl7ES/QF5c2boJD93eMdD+bed+\net73qk5HFZECs2nbbgaOnEHyX0eBlAvWvfH0rYFL9xcFwQRCLWBbuvZvZPylntWY7ZnHmFl9oCmw\nOF33IG8X06tmViHImvPFLde1pm+XiwLtleu3M2DkjMIsQUTCRPz+RHre/yoHDv4JQKlSEbw8qjdn\n1g/N6aXZKZSdVmZWDngXuNvbUgB4ETjdOdcU2Ak8Vxi1pPevQddy+UVpe7AWfLOGB8d9WNhliEgJ\nduhwMtcPmUTcngQAzODxwZ1o0/yMEFd2rGDu0LwdqJuuXdvryzymTlZjzCySlDCY7pz7KHWAc253\nuvGvAHPIxqhRowLPY2JiiImJCaLs4Lw8qjed75rIKu8WdTNmL6ZqxXIM6VNwdyUSkfBw5MgRetz/\naoa5BoN6t81wmf78EhsbS2xs7Am9h6Ue3Mh2gFkEsA5oB+wAlgA9nXNr0o3pCNzpnLvazFoBY51z\nrbzXpgF7nHP3ZnrfGs65nd7zIcCFzrljbgdkZi6nGk9U4qEkOg6cwJbf93qfCY8N7sRN17Yq0M8V\nkZLL7/fTd8RUvliadvXS66+6gGf/WTgXrDMznHO5Ok6b4y4j59xRYBAwH/gZmOmcW2NmA81sgDdm\nLvCrmW0EXgbu8ApqDfQGLjOz5ZlOL/23ma00sxXApcCQ3BSen8pFl+HdcQOoVvlkAJyDkRPm8HHs\nylCVJCLF3H3/fjdDGFzWqhFP3981hBXlLMcthFArjC2EVJu27abLXS/xx4HDAERFRTLlyT60btaw\nUD5fREqGJybOZdI7aZeyvuDsurw3bmChzjUokC2EcNKgzilMH90vMHEtOfkItz0yndWbdoS4MhEp\nLibO+jJDGJxRrxpvPX1rkZh4lpOiX2EhO79xHSb960aiSqXcR+HgoZSrEWo2s4jk5J15yxjzatoc\n29o1KvLu2IFFZuJZThQIWbjkb2fy9NCu+LzZzPv+OETXe15mV3xCiCsTkaLq0//9xPDnPsDvT9nF\nXbVSOd4dO5CK5aNDXFnwFAjZ6NyuGQ+lu7lO3J4Eutz1EvH7E4+zlIiEo9gl67jryVkcOZIyC/nk\nsqV565n+1DylYg5LFi0KhOO4tWsbBt94WaD92879dBk8kYTEwyGsSkSKkm9W/MLAUTNI9u54dlKZ\nUrz+ZN8iNws5GAqEHNzX7woG3PD3QHvz9r10uWuiLoYnIixbvZVbH5pK0p8pYVA6KpLJj9/MhefU\nD21heaRACMKDt3fkxmtbBtobt+7ihnsmkZSUHMKqRCSUftq4nT7DX+PQ4ZTfA1GlInjpkV7F+jR1\nBUKQnri7M12vbBZor960g273vRrYTBSR8LF+cxy97p8cuFhdZISPcSN6hPT2l/lBgZALzwy9nqsv\nPSfQ/nHtNnoPm6zLZouEkS2/76X7vZMCE1h9ET6eHno9HS85J4cliz4FQi74fD4mPNSTy1o1CvQt\nWbmZXkMVCiLhYMvve7lu8ETi/zgEgM9nPH7XtVx3RbMcliweFAi55PP5mPzYzbS+IG0/4bcrfqHH\n/a9y5Ih2H4mUVKmXttmzL+XUczMYMaAjva9pmcOSxYcCIQ98Ph8zxvSj5fmnBfoW//irQkGkhNq0\nbTfX3z2RvfsPAilhMKx/e267oU2IK8tfCoQ88vl8zHymPxc1PT3Qt3TVFrrd+4oONIuUIOs3x9H1\n7rTdRGYw/LYO3NHj0hBXlv8UCCfA5/Px5tO38vfmabuPlv28le73KRRESoL1m+O4Ycgk9gXCwHhw\n4NXc3v2SEFdWMBQIJ8jn8zFtdD8uvfDMQN8Pq7dy/RDNUxApzlZv2sH197zM/oS0A8ijBl1T4nYT\npadAyAc+n48pT/ahbcu0s49+XLuNGxQKIsXSTxu3Zzy11AuDvp0vCnFlBUuBkE98Ph+vPX4zl1/U\nONC3cv12rrnzxcBfGCJS9H334y90GzKJhMSUy9P4fMbjgzvRp1PJDgPQHdPynd/vZ+CoN5j/9epA\nX52alXl//ECqVS4fwspEJCeffbuGfzz2ZuDaRL4IH0/c3YleV7cIcWW5l5c7pikQCoDf72fImHf4\ncOGKQF/1KuV5Z+wA6p1aJYSViUh23l+wnKHPvBe4hHVkZATP/rMrndsVz0lnCoQiZuSE2Uz54NtA\nu1KFaN58uj9nN6gZwqpEJLMpH37Lv178L/6jKVccKFM6khcf7lWsr02kQCiCnp+6kLHTPgu0Ty5b\nmtef7FtsL48rUtKMn/4Zz01dSOqvmbLRUUx5sh8tzq0f0rpOlAKhiJr83lc8PnFu4NZ6ZUqXYuLI\n3hnOShKRwvfoi/9l8ntfB9oVy0fzxtO3cE7DWiGsKn8oEIqwd+YtY/hzH2TYP/nE3Z3o0fHCEFcm\nEn78fj93PTGT/8auCvRVq3wyM5+7jQZ1TglhZflHgVDEzf/6ZwY9PpM/vVnMZnBX78u4r98VIa5M\nJHwkJSXTe9hrfP/TlkBf7RoVeXfswGJ3D+TjUSAUA9/9+Av9H54WuLEGQOfLm/L8sBvw+TQtRKQg\n7YpPoPu9r/DLtj2BvrNOr8nMZ/tTsXx0CCvLfwqEYmLTtt30vO9V4vYmBPpaNT2NqU/0pUyZqBBW\nJlJyrd8cR6+hk9kdfyDQ9/fmDXnt8T5ERUWGsLKCoUAoRvbsS6T7va+wceuuQN+Z9asz69n+VK5Y\nLoSViZQ8Xy/fyIBHZpB4KG3LvHvHvzF6SJcSu2WuQChmkpKSuXnEFBb/+Gugr3rV8rz1TP8Sc2BL\nJNTembeMB57/gL/+Sjmhw+czhtzcjsE3tQtxZQVLgVAMZTWruWx0FP8pATfsFgm1xyZ+zOR3vwrM\nMYgqFcFT93bh+iubh7awQpCXQAhqW8nM2pvZWjNbb2bDshkz3sw2mNkKM2vq9dU2s0Vm9rOZrTKz\nwenGVzKz+Wa2zszmmVmF3BReUvh8PsY90J27bmyLef/rDh5K5raRM3hp5hehLU6kmEpKSuamYa/x\n6jtpYVAuujRTnuobFmGQVzluIZiZD1gPtAN+B5YCPZxza9ON6QAMcs5dbWYtgXHOuVZmVgOo4Zxb\nYWblgGVAJ+fcWjMbA+x1zv3bC5lKzrnhWXx+id5CSO+decsYMfbDDDfX6XTZ+Yx9oFuJ3c8pkt92\n7N5Pz/sn8+tvaWcSnVqtAtNH96NhveohrKxwFdQWQgtgg3Nui3PuL2Am0CnTmE7ANADn3GKggplV\nd87tdM6t8PoTgTVArXTLTPWeTwU656bwkuiGq5oz89nbqFwh7fS3jxb9yLW6hLZIUL778Rfa3zY+\nQxhccHZdFky+J6zCIK+CCYRawLZ07d9I+6We3ZjtmceYWX2gKfCd11XNORcH4JzbCVQLtuiSrPnZ\ndZn36t2cWT/th3fV+u1ccetYVm/aEcLKRIq26bO/48Zhr7Hfu6kNQLf2zXlv3EDKRZcJYWXFR6Gc\nfOvtLnoXuNs5dzCbYdnuFxo1alTgeUxMDDExMflZXpFTrXJ5Pn5pEP947E0WfLMGgF3xB7hu8Es8\neU8XrruieF6OV6Qg+P1+hj//AbPmfh/oi4yMYMSA9tzateTe7jKz2NhYYmNjT+g9gjmG0AoY5Zxr\n77WHA845NybdmInA5865WV57LXCpcy7OzCKB/wKfOOfGpVtmDRDjjanhLX/MaTXhdAwhK8+8Pp8X\n3owNXBgPoEfHC3lqSGcdV5Cwt3NPAv1GTMmw9XxyuTJMfKQXbZqfEcLKQq+gjiEsBRqaWT0ziwJ6\nALMzjZkN3OwV0QrYn7o7CHgNWJ0+DNIt09d73gf4KDeFh4v7+13J+Ad7cFKZUoG+mXOX8n93vMCu\n+ITjLClSsn35/Xqu6p9xV2q9U6vwycS7wj4M8iqoeQhm1h4YR0qATHbOjTazgaRsKUzyxkwA2gMH\ngb7OueVm1hr4ElhFyi4hB4xwzn1qZpWBt4E6wBagm3NufxafHdZbCKk2bomj74PT2LYjPtBXsXw0\nLzzUQz/8EnaefX0BL7wVy1HvhjYAV7VpwvgHuunyLx5NTCvhkpKSGfTEzMBxBYCICB939W7LkD6X\nh7AykcKReCiJ2x6ZwTfLNwX6SpWK4J+3XsmAGy4JYWVFjwIhTEyc9SVPvzY/cG8FgIubNeDlUb0p\nX+6kEFYmUnCWr9nK7f96g52703aVVq1UjkmP3kTzs+uGsLKiSYEQRpas2swd/3qDPfsSA31VKpZl\n7APduORvZ4awMpH85ff7GTd9ERPejM3wR1DzJnV57fE+Je6y1flFgRBm9iccot+DU/lh9dZAn89n\n3HRtK0bd+X86C0mKvV3xCQx4ZAbL16RNc/L5jFu6tubBAR30M34cCoQw5Pf7eW7qQl566wuOpDvA\n1ui06rzy6E3UO7VKCKsTybu5X/7EsGffIyExKdBXuUI040Z011ZwEBQIYWzZ6q3849GM+1ejT4ri\n4TuuptfVLUJYmUjuJCcfYfjz7/Pe/OUZ+i9u1oCXHumlXURBUiCEuUOHk7ln9NvM++rnDP1tWzZi\n7PBu+ockRd6y1Vu5+8lZGU6vjoqKZOgtV+gsolxSIAgAsz75nlEvzOHQ4eRAX8WTT+LRwdfS6bKm\nIaxMJGtHjhzh0ZfmMmPO4gxzC+rXqsIrj96U4dpeEhwFggRs2xHPgJEzjrkgnrYWpKjJaqvAzOh6\nZTPG3NuFyMiSd7/jwqBAkAxST9d7ceYXGe6xoK0FKQqy2yqoXqU8z/yzqw4cnyAFgmRp45Y4Bj0+\nizW/HLu18Nw/r6dyxXIhqkzC1dKfNjNk9DvHbBV0anc+Y4Z00eUn8oECQbKV3dZC2egoBt/YjgE3\ntNE53VLgEg8lMfy5D/j4i1UZruCrrYL8p0CQHGW3tXBm/eo8M7Qr5zeuE6LKpKR7Y85ixkyexx/p\nbmCjrYKCo0CQoPj9fl566wsmvBWb4UwkX4SPLu3O5/HBnYk+Sf84JX+s3xzHfWPeYeX67Rn669So\nxOh7u+hqvQVEgSC5sis+gaFPv0fskvUZ+itViObBgR254armIapMSoKkpGQee3kuM+d+n+EaRKWj\nIrm1a2uG3nKldlMWIAWC5Mln367hwfEfsWPXHxn6G59eg8cGd6LFufVDU5gUS36/nykffsu46YvY\nn3Aow2stzz+NZ4deT52alUNUXfhQIEieHTlyhCcnfcr02d+R/FfaX3NmcMnfzuSJuzvpH7Hk6PPF\n6xj1whw2b9+bob9qpXKMGnQN18ScF6LKwo8CQU7Ypm27GTH2A75b8WuG/qhSEVx/5QU8eHtHykWX\nCVF1UlSt3xzHiLEfsHTVlgz9UVGRdGvfnIcHdtRB40KmQJB88+X36xk5YQ6/bNuTob98uTLcdn0b\n/tHzUs0gFXbs3s8TL3/C3C9/yjC5zMxo2/JMnrynMzVPqRjCCsOXAkHyld/vZ/qcxYydupD4PzLu\nC65cIZrbbvg7A25oo2AIQ7viE1KC4ItVGXYxAjRpeCqP3d1JdzELMQWCFIikpGSenrKAGbO/I+nP\nIxleq1qpXCAYdMZIybdnXyJPTprLnNhVGSY4AlSvWp4RA9rTuV2zEFUn6SkQpEDtik/gqUmf8t/Y\nlcf8VVit8skM6HYJt3a9WMFQAu3Zl8i/X5vHR5+tOOaPgioVy9L/+jbc3v0S/b8vQhQIUih27kng\nyUmfMPfLVfyVKRgqV4ime4cLGdQ7RgefS4D1m+N4dsoCFn239pg/AipXiOaWrq25o/sl2m1YBCkQ\npFClHlD85H8/Z5h4BHBSmVJ0vOQc7ut7BbWqVwpRhZJXXy3bwNjpn/H9T1vJ/O+vYvlo+nW5iH/0\niCEqSkFQVCkQJCS27YjnqVc+YcG3a4/ZrxwR4aN1swbc0+dyHWQs4vx+P29/uoxX3vmKjVt3HfN6\n5QrR9L6mJYN6xugU0mJAgSAhFb8/kbHTF/H+wuUcSHdj9FSn16lKt/Z/o0+ni3StpCJk2454JrwZ\ny6df/XzMzGJIuWvZLV1bc9M1LXWMoBhRIEiRkJSUzGsffMOUD78lbk/CMa+XKV2KSy88kzt6XEKz\ns7TVEAp+v5/3Fy5n2keLWbX+twyXok51fuM6DOp1KVe2bhKCCuVEKRCkSPH7/Xz42Y+8+u5X/Lzx\n9yzHnF6nKte2PZ+brm1F1Uq6UU9B+2njdqZ99B3zvlqd5dZAVFQkbS5owL19ruDcM2uFoELJLwoE\nKbLWb45j4qwvmP/1ag4c/POY130RPs5pUJNrLjuPXle30BlK+WjTtt1M++hbFn67ht927s9yTJ2a\nlel6RTNu7dqa8uVOKuQKpSAUWCCYWXtgLOADJjvnxmQxZjzQATgI9HPOLff6JwP/B8Q5585LN34k\ncBuQevRqhHPu0yzeV4FQghw5coS35n7PWx8vzXarITIygmZn1abTZU3pcnlThUMebNsRz1tzl/LJ\n/3465vIjqVK3Bvpf34bWzRoWcoVS0AokEMzMB6wH2gG/A0uBHs65tenGdAAGOeeuNrOWwDjnXCvv\ntTZAIjAti0A44Jx7LofPVyCUUOs3xzH5va/5fMm6LI81AERG+Gh0Wg0u+dsZXHdFM86sX72Qqywe\n/H4/X36/gdmf/8jilb9muyVgZjQ6rTpXX3IOfbtcrK2BEqygAqEVMNI518FrDwdc+q0EM5sIfO6c\nm+W11wAxzrk4r10PmJNFICQ6557N4fMVCGFgyarNzJizmC+/X8++P47dt52qepXytDivPle1aULb\nFmeG9dbDjt37mfvlTyz8di0/rtvGwUPJ2Y49rXZV2rdpQp/OrXSxuTCRl0AIZlZJLWBbuvZvQIsc\nxmz3+uJyeO9BZnYT8D1wn3PujxzGSwnV4tz6tDi3Pn6/ny+WbuCtuUtZsurXY8Ihbm8Ccz5fyZzP\nV+KL8FGneiXOa1SLNhc05MqLz6JyxZJ7YHr95jgWfL2a71b+yupNO9izL/G44+ueWpm2LRpx87Ut\naVhPW1aSs1BOM3wReNQ558zsceA54NasBo4aNSrwPCYmhpiYmMKoT0LA5/PRtmUj2rZsBMDyNVv5\nYOEKvv5hE7/8tjvD6ZH+o362/L6XLb/vZc7nKxn2LNSsVoEz6lajScNTueDsOrQ49zQqlo8O1dfJ\ns2074vlmxS+sWLOVtb/GsXHrLhKymNuRXpnSkZxzRi3atWpMl8ubaksgzMTGxhIbG3tC7xHsLqNR\nzrn2XjuYXUZrgUuPt8so02dk+7p2GUmq+P2JfPDZChZ+u5a1v+w45pLc2alcIZp6p1bhjHrVaFC3\nGmfUPYWzG9akepXyIZ1o5ff72bh1N2t+2cHGLbvZtG0XG7fuZtvOfRw6nP3un1RmRq3qFWh2Vl06\nXnIuV17cWNcUkoCCOoYQAawj5aDyDmAJ0NM5tybdmI7And5B5VbA2NSDyt7r9Un5hX9uur4azrmd\n3vMhwIXOuV5ZfL4CQbL06/Y9LPhmNd+u+IWfN+7I9sB0dqKiIqlUPprqlU+metXyVCofTZWK5Til\ncjmqVSlPjarlOfWUCpxctgzRZUrl+MvW7/eTnHyEpOQj7P3jIDv3/MGO3X+wOz6R3fEH2LP/IPsT\nDrFzzx/sjj/A/gOHs5wQlp3ICB/1alXh/Ea1+XvzM7j8osY6KCzZKujTTseRdtrpaDMbSMqWwiRv\nzASgPWmnnf7g9b8JxABVSDmmMNI597qZTQOaAn5gMzAwdYsi02crECQou+ITWPTdOn5cu421v8ax\ndUd8jvvZc8MMInw+IiJSHj4zjvr9HD3q56jfZbhj2ImKjPBR45TynFarKmc3qMmF59bn7xc01DWE\nJGiamCaSSULiYRav/JVlP29lw5Y4duxOYFf8Afb9cZAj+fgLPK/KlC5FlYplqV6lPLWqV6RJw1Np\ncd5pnH/mqdr9IydEgSASJL/fz7ad+1izaQfrN8cRF3+A+P0H2ZdwiD8OHObAwSQOHPqTw0nJHDnq\nD/qv/9Rjmr+CAAAJcklEQVStiFKlIih7UmnKRZemfNkyVCwf7e2SKkvtGpVoVL86jU+vqct1SIFR\nIIgUkPTHBw7/mcyhpL84etRP6ahITipTijJRpYI6ziBSWBQIIiIC5C0QdHFzEREBFAgiIuJRIIiI\nCKBAEBERjwJBREQABYKIiHgUCCIiAigQRETEo0AQERFAgSAiIh4FgoiIAAoEERHxKBBERARQIIiI\niEeBICIigAJBREQ8CgQREQEUCCIi4lEgiIgIoEAQERGPAkFERAAFgoiIeBQIIiICKBBERMSjQBAR\nESDIQDCz9ma21szWm9mwbMaMN7MNZrbCzJql659sZnFmtjLT+EpmNt/M1pnZPDOrcGJfRURETkSO\ngWBmPmACcBXQBOhpZo0zjekANHDOnQEMBF5K9/Lr3rKZDQcWOucaAYuAB/L0DUREJF8Es4XQAtjg\nnNvinPsLmAl0yjSmEzANwDm3GKhgZtW99lfAvizetxMw1Xs+Feic+/JFRCS/BBMItYBt6dq/eX3H\nG7M9izGZVXPOxQE453YC1YKoRURECkhkqAtIx2X3wqhRowLPY2JiiImJKYRyRESKj9jYWGJjY0/o\nPcy5bH8PpwwwawWMcs6199rDAeecG5NuzETgc+fcLK+9Frg0dQvAzOoBc5xz56VbZg0Q45yLM7Ma\n3vJnZfH5LqcaRUQkIzPDOWe5WSaYXUZLgYZmVs/MooAewOxMY2YDN3tFtAL2p4ZBam3eI/Myfb3n\nfYCPclO4iIjkrxwDwTl3FBgEzAd+BmY659aY2UAzG+CNmQv8amYbgZeBf6Qub2ZvAt8AZ5rZVjPr\n5700BrjCzNYB7YDR+fi9REQkl3LcZRRq2mUkIpJ7BbXLSEREwoACQUREAAWCiIh4FAgiIgIoEERE\nxKNAEBERQIEgIiIeBYKIiAAKBBER8SgQREQEUCCIiIhHgSAiIoACQUREPAoEEREBFAgiIuJRIIiI\nCKBAEBERjwJBREQABYKIiHgUCCIiAigQRETEo0AQERFAgSAiIh4FgoiIAAoEERHxKBBERARQIIiI\niEeBICIiQJCBYGbtzWytma03s2HZjBlvZhvMbIWZNc1pWTMbaWa/mdkP3qP9iX8dERHJq8icBpiZ\nD5gAtAN+B5aa2UfOubXpxnQAGjjnzjCzlsBEoFUQyz7nnHsuf7+SiIjkRTBbCC2ADc65Lc65v4CZ\nQKdMYzoB0wCcc4uBCmZWPYhl7US/gIiI5I9gAqEWsC1d+zevL5gxOS07yNvF9KqZVQi6ahERyXcF\ndVA5mL/8XwROd841BXYC2nUkIhJCOR5DALYDddO1a3t9mcfUyWJMVHbLOud2p+t/BZiTXQGjRo0K\nPI+JiSEmJiaIskVEwkdsbCyxsbEn9B7mnDv+ALMIYB0pB4Z3AEuAns65NenGdATudM5dbWatgLHO\nuVbHW9bMajjndnrLDwEudM71yuLzXU41iohIRmaGcy5Xx2lz3EJwzh01s0HAfFJ2MU32fqEPTHnZ\nTXLOzTWzjma2ETgI9Dvest5b/9s7PdUPbAYG5qZwERHJXzluIYSathBERHIvL1sImqksIiKAAkFE\nRDwKBBERARQIIiLiUSCIiAigQBAREY8CQUREAAWCiIh4FAgiIgIoEERExKNAEBERQIEgIiIeBYKI\niAAKBBER8SgQREQEUCCIiIhHgSAiIoACQUREPAoEEREBFAgiIuJRIIiICKBAEBERjwJBREQABYKI\niHgUCCIiAigQRETEo0AQERFAgSAiIp6gAsHM2pvZWjNbb2bDshkz3sw2mNkKM2ua07JmVsnM5pvZ\nOjObZ2YVTvzriIhIXuUYCGbmAyYAVwFNgJ5m1jjTmA5AA+fcGcBAYGIQyw4HFjrnGgGLgAfy5RuV\nYLGxsaEuocjQukijdZFG6+LEBLOF0ALY4Jzb4pz7C5gJdMo0phMwDcA5txioYGbVc1i2EzDVez4V\n6HxC3yQM6Ic9jdZFGq2LNFoXJyaYQKgFbEvX/s3rC2bM8Zat7pyLA3DO7QSqBV+2iIjkt4I6qGx5\nWMblexUiIhI859xxH0Ar4NN07eHAsExjJgLd07XXAtWPtyywhpStBIAawJpsPt/poYceeuiR+0dO\nv98zPyLJ2VKgoZnVA3YAPYCemcbMBu4EZplZK2C/cy7OzPYcZ9nZQF9gDNAH+CirD3fO5WVrQ0RE\ncinHQHDOHTWzQcB8UnYxTXbOrTGzgSkvu0nOublm1tHMNgIHgX7HW9Z76zHA22Z2C7AF6Jbv305E\nRIJm3m4ZEREJc0V2pnIwk+FKMjObbGZxZrYyXV/YTeYzs9pmtsjMfjazVWY22OsPx3VR2swWm9ly\nb12M9PrDbl2kMjOfmf1gZrO9dliuCzPbbGY/ej8bS7y+XK+LIhkIwUyGCwOvk/L90wvHyXxHgHud\nc02Ai4A7vZ+FsFsXzrk/gbbOuWZAU6CDmbUgDNdFOncDq9O1w3Vd+IEY51wz51wLry/X66JIBgLB\nTYYr0ZxzXwH7MnWH3WQ+59xO59wK73kiKWen1SYM1wWAc+6Q97Q0KccAHWG6LsysNtAReDVdd1iu\nC1JO9c/8+zzX66KoBkIwk+HCUbVwnsxnZvVJ+cv4O8J0YqO3i2Q5sBNY4JxbSpiuC+B5YCgpoZgq\nXNeFAxaY2VIz6+/15XpdBHPaqRRdYXNGgJmVA94F7nbOJZpZ5u8eFuvCOecHmplZeeADM2vCsd+9\nxK8LM7saiHPOrTCzmOMMLfHrwtPaObfDzE4B5pvZOvLwc1FUtxC2A3XTtWt7feEuzrtGFGZWA9gV\n4noKhZlFkhIG051zqfNVwnJdpHLOJQCxQHvCc120Bq41s1+At4DLzGw6sDMM1wXOuR3ef3cDH5Ky\n2z3XPxdFNRACk+HMLIqUCW2zQ1xTKBgZLwOSOpkPjjOZrwR6DVjtnBuXri/s1oWZVU09U8TMTgKu\nIOWYStitC+fcCOdcXefc6aT8fljknLsJmEOYrQszi/a2oDGzssCVwCry8HNRZOchmFl7YBxpE9pG\nh7ikQmVmbwIxQBUgDhhJSvK/A9TBm8znnNsfqhoLg5m1Br4k5Qc8dUr+CGAJ8DbhtS7OJeXgoM97\nzHLOPWFmlQmzdZGemV0K3OecuzYc14WZnQZ8QMq/jUjgDefc6LysiyIbCCIiUriK6i4jEREpZAoE\nEREBFAgiIuJRIIiICKBAEBERjwJBREQABYKIiHgUCCIiAsD/A88YLwXMU/sQAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "socks.Update('pair')\n", "thinkplot.Pdf(socks)\n", "thinkplot.Config(ylim=[0,0.03])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Chess-playing twins\n", "\n", "Allen Downey\n", "\n", "Two identical twins are members of my chess club, but they never show up on the same day; in fact, they strictly alternate the days they show up. I can't tell them apart except that one is a better player than the other: Avery beats me 60% of the time and I beat Blake 70% of the time. If I play one twin on Monday and win, and the other twin on Tuesday and lose, which twin did I play on which day?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AB 0.222222222222\n", "BA 0.777777777778\n" ] } ], "source": [ "# Solution\n", "\n", "pmf = Pmf(['AB', 'BA'])\n", "pmf['AB'] = 0.4 * 0.3\n", "pmf['BA'] = 0.7 * 0.6\n", "pmf.Normalize()\n", "pmf.Print()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "class Chess(Suite):\n", " \n", " prob_I_beat = dict(A=0.4, B=0.7)\n", " \n", " def Likelihood(self, data, hypo):\n", " \"\"\"Probability of data under hypo.\n", " \n", " data: sequence of 'W' and 'L'\n", " hypo: sequence of 'A' and 'B'\n", " \"\"\"\n", " total = 1\n", " for outcome, twin in zip(data, hypo):\n", " like = self.prob_I_beat[twin]\n", " if outcome == 'W':\n", " total *= like\n", " else:\n", " total *= 1-like\n", " return total" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AB 0.222222222222\n", "BA 0.777777777778\n" ] } ], "source": [ "# Solution\n", "\n", "chess = Chess(['AB', 'BA'])\n", "chess.Update('WL')\n", "chess.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1984\n", "\n", "by Katerina Zoltan\n", "\n", "The place: Airstrip One. The reason: thoughtcrime. The time: ???\n", "\n", "John's parents were taken by the Thought Police and erased from all records. John is being initiated into the Youth League and must pass a test. He is asked whether his parents are good comrades. It is not clear what John's admission officer knows:\n", "\n", "1. He may know that John's parents have been erased and that John did not give them away.\n", "2. He may know only that John's parents have been erased.\n", "3. He may not know that John's parents have been erased.\n", "\n", "It is a well known fact that children who have parents that are 'good comrades' have twice the chances of passing the test. However, if the admission officer knows that their parents committed thoughtcrime (but not that they have been erased yet), a child that gave his parents away has three times the chances of getting in than a child who did not give them away.\n", "\n", "And if the admission officer knows the specifics of the arrest, a child that denies that the records are false and their parents existed has a 1/3 chance of getting in, while one who pretends that his parents never existed has a 2/3 chance. Lying to an admission officer that knows the parents have been erased will ensure that the child does not get in. Telling an admission officer that your parents do not exist when he does not know this will give you a 1/3 chance of getting in.\n", "\n", "There is a 60% chance the admission officer knows nothing, a 25% chance that he knows the parents have been erased, and a 15% chance that the officer knows all of the details. John says that he never had parents and is admitted into the Youth League. What did his admission officer know?" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "everything 0.15\n", "nothing 0.6\n", "something 0.25\n" ] } ], "source": [ "# Solution\n", "\n", "officer = {'everything':0.15, 'something':0.25, 'nothing':0.6}\n", "\n", "class ThoughtPolice(Suite):\n", "\n", " def Likelihood(self, data, hypo):\n", " if data == 'gave away':\n", " if hypo == 'everything':\n", " return 0\n", " elif hypo == 'something':\n", " return 1\n", " else:\n", " return 3\n", " elif data == 'none':\n", " if hypo == 'everything':\n", " return 2\n", " elif hypo == 'something':\n", " return 1\n", " else:\n", " return 1\n", " else: # data == 'good comrades'\n", " if hypo == 'everything':\n", " return 0\n", " elif hypo == 'something':\n", " return 0\n", " else:\n", " return 2\n", " \n", "pmf = ThoughtPolice(officer)\n", "pmf.Print()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "everything 0.260869565217\n", "nothing 0.521739130435\n", "something 0.217391304348\n" ] } ], "source": [ "# Solution\n", "\n", "pmf.Update('none')\n", "pmf.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Where Am I? - The Robot Localization Problem\n", "\n", "by Kathryn Hite\n", "\n", "Bayes's Theorem proves to be extremely useful when building mobile robots that need to know where they are within an environment at any given time. Because of the error in motion and sensor systems, a robot's knowledge of its location in the world is based on probabilities. Let's look at a simplified example that could feasibly be scaled up to create a working localization model.\n", "\n", "**Part A:** We have a robot that exists within a very simple environement. The map for this environment is a row of 6 grid cells that are colored either green or red and each labeled $x_1$, $x_2$, etc. In real life, a larger form of this grid environment could make up what is known as an occupancy grid, or a map of the world with places that the robot can go represented as green cells and obstacles as red cells.\n", "\n", "|G|R|R|G|G|G|\n", "|-|-|-|-|-|-|\n", "|$x_1$|$x_2$|$x_3$|$x_4$|$x_5$|$x_6$|\n", "\n", "The robot has a sensor that can detect color with an 80% chance of being accurate.\n", "\n", "Given that the robot gets dropped in the environment and senses **red**, what is the probability of it being in each of the six locations?" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0.0833333333333\n", "1 0.333333333333\n", "2 0.333333333333\n", "3 0.0833333333333\n", "4 0.0833333333333\n", "5 0.0833333333333\n" ] } ], "source": [ "# Solution\n", "\n", "colors = 'GRRGGG'\n", "locs = range(len(colors))\n", "data = 'R'\n", "\n", "pmf = Pmf(locs)\n", "for hypo in pmf:\n", " if colors[hypo] == data:\n", " pmf[hypo] *= 0.8\n", " else:\n", " pmf[hypo] *= 0.2\n", "pmf.Normalize()\n", "pmf.Print()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solution\n", "\n", "class Robot(Suite):\n", " \n", " colors = 'GRRGGG'\n", " \n", " def Likelihood(self, data, hypo):\n", " \"\"\"\n", " \n", " data: 'R' or 'G'\n", " hypo: index of starting location\n", " \"\"\"\n", " if self.colors[hypo] == data:\n", " return 0.8\n", " else:\n", " return 0.2" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0.0833333333333\n", "1 0.333333333333\n", "2 0.333333333333\n", "3 0.0833333333333\n", "4 0.0833333333333\n", "5 0.0833333333333\n" ] } ], "source": [ "# Solution\n", "\n", "robot = Robot(locs)\n", "robot.Update('R')\n", "robot.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Part B:** This becomes an extremely useful tool as we begin to move around the map. Let's try to get a more accurate knowledge of where the robot falls in the world by telling it to move forward one cell.\n", "\n", "The robot moves forward one cell from its previous position and the sensor reads **green**, again with an 80% accuracy rate. Update the probability of the robot having started in each location." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "class Robot2(Suite):\n", " \n", " colors = 'GRRGGG'\n", " \n", " def Likelihood(self, data, hypo):\n", " \"\"\"\n", " \n", " data: tuple (offset, 'R' or 'G')\n", " hypo: index of starting location\n", " \"\"\"\n", " offset, color = data\n", " index = (hypo + offset) % len(self.colors)\n", " if self.colors[index] == color:\n", " return 0.8\n", " else:\n", " return 0.2" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0.0833333333333\n", "1 0.333333333333\n", "2 0.333333333333\n", "3 0.0833333333333\n", "4 0.0833333333333\n", "5 0.0833333333333\n" ] } ], "source": [ "# Solution\n", "\n", "robot = Robot2(locs)\n", "robot.Update((0, 'R'))\n", "robot.Print()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0.030303030303\n", "1 0.121212121212\n", "2 0.484848484848\n", "3 0.121212121212\n", "4 0.121212121212\n", "5 0.121212121212\n" ] } ], "source": [ "# Solution\n", "\n", "robot.Update((1, 'G'))\n", "robot.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Red Dice problems\n", "\n", "Suppose I have a six-sided die that is red on 2 sides and blue on 4 sides, and another die that's the other way around, red on 4 sides and blue on 2.\n", "\n", "I choose a die at random and roll it, and I tell you it came up red. What is the probability that I rolled the second die (red on 4 sides)?" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Blue 2/3\n", "Red 1/3\n" ] } ], "source": [ "# Solution\n", "\n", "from fractions import Fraction\n", "\n", "d1 = Pmf({'Red':Fraction(2), 'Blue':Fraction(4)}, label='d1 (bluish) ')\n", "d1.Print()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Blue 1/3\n", "Red 2/3\n" ] } ], "source": [ "# Solution\n", "\n", "d2 = Pmf({'Red':Fraction(4), 'Blue':Fraction(2)}, label='d2 (reddish)')\n", "d2.Print()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "d2 (reddish) 1/2\n", "d1 (bluish) 1/2\n" ] } ], "source": [ "# Solution\n", "\n", "dice = Pmf({d1:Fraction(1), d2:Fraction(1)})\n", "dice.Print()" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "class Dice(Suite):\n", " def Likelihood(self, data, hypo):\n", " \"\"\"\n", " data: 'Red' or 'Blue'\n", " hypo: a Die object\n", " \"\"\"\n", " return hypo[data]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "d2 (reddish) 1/2\n", "d1 (bluish) 1/2\n" ] } ], "source": [ "# Solution\n", "\n", "prior = Dice({d1:Fraction(1), d2:Fraction(1)})\n", "prior.Print()" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "d2 (reddish) 2/3\n", "d1 (bluish) 1/3\n" ] } ], "source": [ "# Solution\n", "\n", "posterior = prior.Copy()\n", "posterior.Update('Red')\n", "posterior.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scenario B\n", "\n", "Suppose I roll the same die again. What is the probability I get red?" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Blue 4/9\n", "Red 5/9\n" ] } ], "source": [ "# Solution\n", "\n", "from thinkbayes2 import MakeMixture\n", "\n", "predictive = MakeMixture(posterior)\n", "predictive.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scenario A\n", "\n", "Instead of rolling the same die, suppose I choosing a die at random and roll it. What is the probability that I get red?" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Blue 1/2\n", "Red 1/2\n" ] } ], "source": [ "# Solution\n", "\n", "from thinkbayes2 import MakeMixture\n", "\n", "predictive = MakeMixture(prior)\n", "predictive.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scenario C\n", "\n", "Now let's run a different experiment. Suppose I choose a die and roll it. If the outcome is red, I report the outcome. Otherwise I choose a die again and roll again, and repeat until I get red.\n", "\n", "What is the probability that the last die I rolled is the reddish one?" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solution\n", "\n", "# On each roll, there are four possible results, with these probabilities:\n", "\n", "# d1, red 1/2 * 1/3\n", "# d1, blue 1/2 * 2/3\n", "# d2, red 1/2 * 2/3\n", "# d2, blue 1/2 * 1/3\n", "\n", "#On the last roll, I tell you that the outcome is red, so we are left with two possibilities:\n", "\n", "# d1, red 1/2 * 1/3\n", "# d2, red 1/2 * 2/3\n", "\n", "# The likelihood ratio is 2 to 1, so we can use that to update the prior:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "d2 (reddish) 2/3\n", "d1 (bluish) 1/3\n" ] } ], "source": [ "# Solution\n", "\n", "posterior = prior.Copy()\n", "posterior[d1] *= 1\n", "posterior[d2] *= 2\n", "posterior.Normalize()\n", "posterior.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scenario D\n", "\n", "Finally, suppose I choose a die and roll it over and over until I get red, then report the outcome. What is the probability that the die I rolled is the reddish one?" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "d2 (reddish) 1/2\n", "d1 (bluish) 1/2\n" ] } ], "source": [ "# Solution\n", "\n", "# In this case, the likelihood of the data is the same regardless of\n", "# which die I rolled, so the posterior is the same as the prior.\n", "\n", "posterior = prior.Copy()\n", "posterior.Print()" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solution\n", "\n", "#In summary, each of the four scenarios yields a different pair of posterior\n", "# and predictive distributions.\n", "\n", "# Scenario Posterior probability of d2 Predictive probability of red\n", "# A 2/3 1/2\n", "# B 2/3 5/9\n", "# C 2/3 1\n", "# D 1/2 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The bus problem\n", "\n", "Allen Downey\n", "\n", "Two buses routes run past my house, headed for Arlington and Billerica. In theory, the Arlington bus runs every 20 minutes and the Billerica bus every 30 minutes, but by the time they get to me, the time between buses is well-modeled by exponential distributions with means 20 and 30.\n", "\n", "Part 1: Suppose I see a bus outside my house, but I can't read the destination. What is the probability that it is an Arlington bus?\n", "\n", "Part 2: Suppose I see a bus go by, but I don't see the destination, and 3 minutes later I see another bus. What is the probability that the second bus is going to Arlington?" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "def generate_times(lam, n=10):\n", " gaps = np.random.exponential(lam, n)\n", " times = np.cumsum(gaps)\n", " for time in times:\n", " yield time" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "15.1381021981\n", "35.6595368992\n", "40.2555876685\n", "44.5338062274\n", "58.7218133396\n", "82.9148030062\n", "113.755194822\n", "126.646534094\n", "138.89916925\n", "145.582801415\n" ] } ], "source": [ "# Solution\n", "\n", "for time in generate_times(20, 10):\n", " print(time)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "def generate_buses(names, lams, n):\n", " buses = [generate_times(lam, n) for lam in lams]\n", " times = [next(bus) for bus in buses]\n", "\n", " while True:\n", " i = np.argmin(times) \n", " yield(names[i], times[i])\n", " times[i] = next(buses[i])" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "('A', 27.901039669324383)" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "next(generate_buses('AB', [20, 30], 10))" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Solution\n", "\n", "res = []\n", "for bus, time in generate_buses('AB', [20, 30], 1000):\n", " res.append((bus, time))" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solution\n", "\n", "buses, times = zip(*res)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.6045949214026602" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "hist = Hist(buses)\n", "hist['A'] / hist.Total()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }