{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Think Bayes: Chapter 3\n", "\n", "This notebook presents example code and exercise solutions for Think Bayes.\n", "\n", "Copyright 2016 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "% matplotlib inline\n", "\n", "import thinkplot\n", "from thinkbayes2 import Hist, Pmf, Suite, Cdf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Dice problem\n", "\n", "Suppose I have a box of dice that contains a 4-sided die, a 6-sided\n", "die, an 8-sided die, a 12-sided die, and a 20-sided die.\n", "\n", "I select a die from the box at random, roll it, and get a 6.\n", "What is the probability that I rolled each die?\n", "\n", "The `Dice` class inherits `Update` and provides `Likelihood`" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class Dice(Suite):\n", " def Likelihood(self, data, hypo):\n", " if hypo < data:\n", " return 0\n", " else:\n", " return 1/hypo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the update looks like:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4 0.0\n", "6 0.3921568627450981\n", "8 0.29411764705882354\n", "12 0.19607843137254904\n", "20 0.11764705882352944\n" ] } ], "source": [ "suite = Dice([4, 6, 8, 12, 20])\n", "suite.Update(6)\n", "suite.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's what it looks like after more data:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4 0.0\n", "6 0.0\n", "8 0.9432484536722124\n", "12 0.0552061280612909\n", "20 0.001545418266496554\n" ] } ], "source": [ "for roll in [6, 8, 7, 7, 5, 4]:\n", " suite.Update(roll)\n", " \n", "suite.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The train problem\n", "\n", "The Train problem has the same likelihood as the Dice problem." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class Train(Suite):\n", " def Likelihood(self, data, hypo):\n", " if hypo < data:\n", " return 0\n", " else:\n", " return 1/hypo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But there are many more hypotheses" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.0028222671142652746" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hypos = range(1, 1001)\n", "suite = Train(hypos)\n", "suite.Update(60)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the posterior looks like" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl0XeV97vHvz5ony7Isz7ItD2DLDDYoxgxJG0zAJClOCmlMmpS29NKsC23a5CYXmrvSlru4DbdpSO4qSUMgKaElhjgpuIRACEOAlNjIzJ6wjLEtbGPZkmXN4+/+cbaPjmXJ2jo60tE5ej5raZ293/3urffVtvVo73cP5u6IiIhMSnYDRERkfFAgiIgIoEAQEZGAAkFERAAFgoiIBBQIIiICKBBERCSgQBAREUCBICIigcxkN2A4pk2b5gsWLEh2M0REUsbWrVuPuntZmLopFQgLFiyguro62c0QEUkZZrYvbF2dMhIRESBkIJjZWjPbZWY1ZnbrAMtzzOyhYPlmM1sQs+y2oHyXmV0VUz7FzDaa2U4z22FmFyeiQyIiEp8hA8HMMoC7gauBSuB6M6vsV+1GoMHdFwN3AXcG61YC64HlwFrgO8H2AL4NPOHuS4HzgR0j746IiMQrzBHCKqDG3d9x905gA7CuX511wP3B9EZgjZlZUL7B3TvcfS9QA6wys8nAh4D7ANy9092Pj7w7IiISrzCBMAc4EDNfG5QNWMfdu4FGoPQM6y4E6oAfmtmrZnavmRXE1QMREUmIMIFgA5T1f6vOYHUGK88ELgC+6+4rgRbgtLEJADO7ycyqzay6rq4uRHNFRCQeYQKhFiiPmZ8LHBysjpllAsVA/RnWrQVq3X1zUL6RSECcxt3vcfcqd68qKwt1Ke2oOXz0BJ1d3Ultg4jIaAkTCC8DS8yswsyyiQwSb+pXZxNwQzB9HfCMR97NuQlYH1yFVAEsAba4+2HggJmdHayzBtg+wr6Mql/+Zjs3/+8H+e+3P0h7R1eymyMiknBDBkIwJnAL8CSRK4EedvdtZna7mV0TVLsPKDWzGuCLBKd/3H0b8DCRX/ZPADe7e0+wzl8A/25mbwArgP+TuG4l3vcefh6AhhOt/MfTryW5NSIiiRfqTmV3fxx4vF/Z12Km24FPDbLuHcAdA5S/BlQNp7HjRe3hhmQ3QUQk4XSnchwaTrQmuwkiIgmnQIjDcQWCiKQhBUIc6htbkt0EEZGEUyCElJOdFZ3u6u45Q00RkdSkQAipZHJespsgIjKqFAghFRflnzLf09ObpJaIiIwOBUJImRmn/qgam9uS1BIRkdGhQIhT/XENLItIelEghBR5Ekefel16KiJpRoEQJ92LICLpRoEQp2O6F0FE0owCIU4NCgQRSTMKhDgdP6GrjEQkvSgQQuo3pszR483JaYiIyChRIMTpmAJBRNKMAiEk7/ca6ebWDr05TUTSigJhBOoadJQgIulDgTACOm0kIulEgRBS/0FlgKM6QhCRNKJAGAGdMhKRdKJAGAEdIYhIOlEgjMDRhqZkN0FEJGEUCCOgIwQRSScKhJD6P/4aImMIA5WLiKQiBcII9PT0crxJzzQSkfQQKhDMbK2Z7TKzGjO7dYDlOWb2ULB8s5ktiFl2W1C+y8yuiil/18zeNLPXzKw6EZ0ZTYMdCWgcQUTSxZCBYGYZwN3A1UAlcL2ZVfardiPQ4O6LgbuAO4N1K4H1wHJgLfCdYHsnfdjdV7h71Yh7MoaKCnKj00fqNY4gIukhzBHCKqDG3d9x905gA7CuX511wP3B9EZgjZlZUL7B3TvcfS9QE2wvpc0onRydfv/oiSS2REQkccIEwhzgQMx8bVA2YB137wYagdIh1nXgl2a21cxuGn7Tk2fGtJhAOKZAEJH0kBmijg1Q1v+E+mB1zrTupe5+0MymA0+Z2U53f/60bx4Ji5sA5s2bF6K5o2+WAkFE0lCYI4RaoDxmfi5wcLA6ZpYJFAP1Z1rX3U9+HgH+g0FOJbn7Pe5e5e5VZWVlIZo7OmLHlGdOK45OH65TIIhIeggTCC8DS8yswsyyiQwSb+pXZxNwQzB9HfCMRy7L2QSsD65CqgCWAFvMrMDMigDMrAC4Enhr5N0ZG2VTC6OHPkcbmuju7klqe0REEmHIU0bu3m1mtwBPAhnAD9x9m5ndDlS7+ybgPuABM6shcmSwPlh3m5k9DGwHuoGb3b3HzGYA/xEZdyYTeNDdnxiF/o2K7KxMpk4p4NjxFpzIDWqzyoqHXE9EZDwLM4aAuz8OPN6v7Gsx0+3ApwZZ9w7gjn5l7wDnD7exydT/PoQZpZM5drwFiIwjKBBEJNXpTuU4mJ16pZHGEUQkHSgQ4nTKvQi60khE0oACIU6zYq40UiCISDpQIMTBMKaXFkXnD9U1JrE1IiKJoUAIqf+deLGDyIfqGvUYbBFJeQqEOBUV5DK5MA+Aru4evV9ZRFKeAiGkgY4AZk/vO0p47/3jY9kcEZGEUyCMwJzpU6LTB48oEEQktSkQ4hDcYc2cGSXRMh0hiEiqUyCENNCY8SmnjI40jGFrREQST4EQh+AA4ZRTRjpCEJFUp0AYgRmlk8nIiPwIG0600tbemeQWiYjET4EwAhkZk065Y1lHCSKSyhQIcTg5qAz9xxEUCCKSuhQIIQ12J3L5zKnR6f2H6seqOSIiCadAGKF5s/sCYd/BY0lsiYjIyCgQRmj+7NLotI4QRCSVKRBCGuyU0axpfVcaHTveQnNrx1g2S0QkYRQIcYgZUyYzM+OUcQSdNhKRVKVASIB5s/oeYaHTRiKSqhQIcYi97BROHUfQEYKIpCoFQgLMmxV76ameaSQiqUmBENKZ3oc2v9+lp3p7moikIgVCAkwtLqAwPweA9o4uvT1NRFJSqEAws7VmtsvMaszs1gGW55jZQ8HyzWa2IGbZbUH5LjO7qt96GWb2qpk9NtKOjLoz/NVvZiyY0zeOsGd/3Vi0SEQkoYYMBDPLAO4GrgYqgevNrLJftRuBBndfDNwF3BmsWwmsB5YDa4HvBNs76QvAjpF2YjxYVF4WnX7ngAJBRFJPmCOEVUCNu7/j7p3ABmBdvzrrgPuD6Y3AGotcirMO2ODuHe6+F6gJtoeZzQU+Btw78m6MNTutZNG86dHpGh0hiEgKChMIc4ADMfO1QdmAddy9G2gESodY91vAV4DeYbc6CYYaJ449QthzoE4DyyKScsIEwul/Dp9+0c1gdQYsN7OPA0fcfeuQ39zsJjOrNrPqurrx8Ze3DdCrGaVFFORFBpZb2jp4/1jTGLdKRGRkwgRCLVAeMz8XODhYHTPLBIqB+jOseylwjZm9S+QU1OVm9m8DfXN3v8fdq9y9qqysbKAq44KZnXaUICKSSsIEwsvAEjOrMLNsIoPEm/rV2QTcEExfBzzjkXMmm4D1wVVIFcASYIu73+buc919QbC9Z9z9swnoT1ItnhcTCPuPJLElIiLDlzlUBXfvNrNbgCeBDOAH7r7NzG4Hqt19E3Af8ICZ1RA5MlgfrLvNzB4GtgPdwM3u3jNKfRkz/R9dcdJCHSGISAobMhAA3P1x4PF+ZV+LmW4HPjXIuncAd5xh288Bz4VpRzKFGSQ+5QjhwFHcfdDwEBEZb3SncgJNKymkuCgPgLb2Tj3XSERSigIhpDBHCGbG0oqZ0fmd7xwazSaJiCSUAiHBli6MCYS9h5PYEhGR4VEgxOFMwwKnHiEoEEQkdSgQEmzh3GlkZUYe13Skvon6xpYkt0hEJBwFQhzOdOVQZmYGS+b3Pddoh44SRCRFKBBCGs6TiZYtnBWd3qVxBBFJEQqEURA7sKwjBBFJFQqEUXB2xYzoU/32HqijubUjqe0REQlDgRDScB5nXZCXE30/ggNv7X5vlFolIpI4CoQ4hHkYxXln9b0y4s23FQgiMv4pEEbJuQoEEUkxCoRRsnThTDKD+xHeO3KcY8ebk9wiEZEzUyDEIcwTTLOzMllaMSM6/9bu/u8UEhEZXxQIIcXziuRzz5obnX59V20CWyMikngKhFG04uy+QHh1xwF6e3uT2BoRkTNTIIQ0nMtOT1o0ryz6foQTzW3U7Ndb1ERk/FIgjCIz48LK+dH56rf2JbE1IiJnpkCIw3Bei3nh8nnR6eptCgQRGb8UCCHFM6gMcP7Zc8nIiPyY9x08xtEGXX4qIuOTAiEOwzhAIC83m+WLZkfnt+ooQUTGKQXCGKg6J2YcQYEgIuOUAmEMXLi8LxBe31VLS5uefioi448CIQ7DGVQGmDltMgvLywDo6ellyxvvjkKrRERGRoEQkg/rnWmnu2TFwuj0f722Z6TNERFJuFCBYGZrzWyXmdWY2a0DLM8xs4eC5ZvNbEHMstuC8l1mdlVQlmtmW8zsdTPbZmZ/n6gOjVeXrFwUnX5tZ61emiMi486QgWBmGcDdwNVAJXC9mVX2q3Yj0ODui4G7gDuDdSuB9cByYC3wnWB7HcDl7n4+sAJYa2arE9Ol0RHPncqxZpROZnHw0pze3l62vLE3Ec0SEUmYMEcIq4Aad3/H3TuBDcC6fnXWAfcH0xuBNRY50b4O2ODuHe6+F6gBVnnEyQvys4Kvkf3GTQGxRwm/eVWnjURkfAkTCHOAAzHztUHZgHXcvRtoBErPtK6ZZZjZa8AR4Cl33zzQNzezm8ys2syq6+rGx7OAhjek3Cd2HOGNXbU0nGhNTINERBIgTCAM9Puv/1/zg9UZdF1373H3FcBcYJWZnTPQN3f3e9y9yt2rysrKQjR3/CqbWkTlolkA9Lrz65ffTnKLRET6hAmEWqA8Zn4u0P9tL9E6ZpYJFAP1YdZ19+PAc0TGGNLe5RctjU4/u3nXiMcmREQSJUwgvAwsMbMKM8smMki8qV+dTcANwfR1wDMe+U23CVgfXIVUASwBtphZmZlNATCzPOAKYOfIuzN6Yn9vD/c+hFgXr1hITnYWALXvN7B735GRNk1EJCGGDIRgTOAW4ElgB/Cwu28zs9vN7Jqg2n1AqZnVAF8Ebg3W3QY8DGwHngBudvceYBbwrJm9QSRwnnL3xxLbtfEpNyeLS1b2jSU8s3lc56CITCCZYSq5++PA4/3KvhYz3Q58apB17wDu6Ff2BrByuI1NF5dftJRnN+8C4MVX9vDHn7iE3JysJLdKRCY63akcUiLP9S9bOJNZZcUAtLV38nz17oRtW0QkXgqEJDAzrrp0eXT+8Rfe0uCyiCSdAiEOIxhTjrp89dlkZ0XO2B04VM/2PYdGvlERkRFQICRJQV4Ov7vqrOj8479+M4mtERFRIMRlJJedxlp7Wd+9eFvefJcj9U0J2a6ISDwUCCGNxin++bOncs6SyOs1e935z2dfT/w3EREJSYGQZJ9Y03f17VP/tYPGprYktkZEJjIFQkgjfUHOYFYsncuCOdMA6Oru4ecaSxCRJFEgJJmZce2VfUcJj7/wlt65LCJJoUCIQ6IGlU9afV4Fs2NuVPvFC9sSun0RkTAUCCGN5n1jkyZN4pNX9B0lPPr0a3rFpoiMOQXCOPGhqiXRx1m0tnfyyK9eTXKLRGSiUSDEIcFnjADIzMzg+o+tis4/9us3OXa8+QxriIgklgJhHLlkxUIq5vZdcfSTJ7cmuUUiMpEoEOJgcb9VeYjtmvG5a1ZH559+aSf7Dh4ble8lItKfAiGksXoa6XlnzeG8s+YCkbuX7934op6EKiJjQoEwzpgZf3rtpUyaFNk12/cc4jev7Elyq0RkIlAghDSWf6WXzyzhYx/qe/Dd/Y++RHtH15h9fxGZmBQI49QfrK2iuCgPgPrGFn7885eT3CIRSXcKhDiMxmWn/eXnZXPDuouj8z//9Rvs2nt49L+xiExYCoRx7ENVS1ixtBwAB+5+8Dk6u7qT2iYRSV8KhHHMzPj8pz9ETnYWAO8dOc7Dv6hOcqtEJF0pEEKKHVNO9MPtzqRsahF/FHNvwiNPv8Zbu98bs+8vIhOHAiEFXHVZZfTNag5860dPc6JZL9IRkcQKFQhmttbMdplZjZndOsDyHDN7KFi+2cwWxCy7LSjfZWZXBWXlZvasme0ws21m9oVEdSgdmRlf+NwaigpyAWg40crdDz6nG9ZEJKGGDAQzywDuBq4GKoHrzayyX7UbgQZ3XwzcBdwZrFsJrAeWA2uB7wTb6wa+5O7LgNXAzQNsc1xJ9i/fqcUF3PKHH47OV2/bx2PP6e1qIpI4YY4QVgE17v6Ou3cCG4B1/eqsA+4PpjcCayxyon0dsMHdO9x9L1ADrHL3Q+7+CoC7NwE7gDkj7056q1o+n4//znnR+R89+hJvvq3xBBFJjDCBMAc4EDNfy+m/vKN13L0baARKw6wbnF5aCWwO3+zkGstB5f4++3sXsXjedCDyrKNv/PCXvH/sRNLaIyLpI0wgDPTbr//5k8HqnHFdMysEfgr8lbsP+FvNzG4ys2ozq66rqwvR3PSWlZXB//yzq5hSlA9Ac2sHX//+E3q0hYiMWJhAqAXKY+bnAgcHq2NmmUAxUH+mdc0si0gY/Lu7/2ywb+7u97h7lbtXlZWVhWhu+ptaXMBXbrySjIzI7tt/qJ5v/PCXdHf3JLllIpLKwgTCy8ASM6sws2wig8Sb+tXZBNwQTF8HPOORUdhNwPrgKqQKYAmwJRhfuA/Y4e7fTERHRlvsIVESzxhFnV0xkz//gw9G51/dcYDvPvR80ge/RSR1DRkIwZjALcCTRAZ/H3b3bWZ2u5ldE1S7Dyg1sxrgi8CtwbrbgIeB7cATwM3u3gNcCnwOuNzMXgu+PprgvqW9NauXcd1VF0bnn9uyiwcf25LEFolIKssMU8ndHwce71f2tZjpduBTg6x7B3BHv7IXGXh8QYZp/dVV1B9v4ZnNOwH42a9eJT8vm09esTLJLRORVKM7lUMar6diTj7v6MLK+dGyf/vPzTzy9GtJbJWIpCIFQhrIyJjEl/7kCpYvnh0te2DTb3n0mdeT2CoRSTUKhJBOebjdODzblZOdxd/cdPUpofCjR1/i4Seqx+3RjYiMLwqENJKbEwmFykWzomUP/aKa7//kRXp7e5PYMhFJBQqEOIyHy04Hk5uTxVf//KOcf/bcaNmTv9nGP/3wKb1cR0TOSIGQhk4eKVx24eJo2W/f2Mvf3f0Yx5tak9gyERnPFAhpKjMzg7/63JpTHoa3a+9hvvKNn7Jnvx4BIiKnUyCEFDswm8yH2w2HmfHHn7yYP/7EJdFh8GPHW/jqtx/hherdSW2biIw/CoQ0Z2b83ofP46uf/xj5udkAdHX38K0HnuZfHvq1xhVEJEqBMEGsXFbO17/0+8yZPiVa9tR/7eAr//QzDhxuSGLLRGS8UCCElA7X8s+ZPoWvf/H3uWTlomjZgUP1fPkfN/LEC9vSoo8iEj8FwgSTn5fNF2+4gs9/+kNkZWYAkVNI39/4An9393/qZTsiE5gCIQ4pMqY8KDPjI5dUcueXrqV8Zkm0/K3dB/nrr/9ERwsiE5QCYQKbP3sq//d/XMsn16yIXoXU0dnF9ze+wN986xH21h5NavtEZGwpEOKQKpedhpGdlclnr1nNP3zxk8yd0Xe08Pa77/Plf9zIPQ+/QFNLexJbKCJjRYEQUrqfQFkyfwbf+PJ1XPuRC6Kv5nQij734izs28IsX3tIrOkXSnAJBorKyMvjMx1dx161/wIqlfa/Cbmpp596NL/KFf3iIF1+p0fiCSJpSIIQ0kX4Jzpk+hf/1+Y/ylRuvoqykKFp++OgJ7rr/V3z5Gz/lle37J9TPRGQiCPUKTZl4zIyLzqtg5bJyfvHCNn76y1doaesAYG/tUe743uMsKi/juqsu5APnzE+rcRWRiUqBEIfx+IKc0ZKdlcm6y89nzeqlPPKrV3ns12/SFYwl7DlQx533PsG8WVO59soLuPj8hdHxBxFJPQoECaUwP4fPXrOatR88h0eefo2nXtoRHWTef6ieu+7/FQ+UFPLRD53LmtVLKczPSXKLRWS49OecDMu0kkL+7LrL+O7XPsM1Hz6fnOys6LKjDc386NGX+G9fe4B7Hn6B2vf1jCSRVKIjhJBOeafyxDljNKipxQXc8ImL+eQVK/j582/x5IvbovcrdHZ18+RvtvHkb7ZRuWgWV1y8jNXnV5wSHiIy/igQZEQmF+Zx/Uc/wLUfWckLW3fz2HNvsv9QfXT59j2H2L7nEPduzOaDFy7hI5cso2LutCS2WEQGo0CQhMjOymTN6mVcftFS3tp9kJ//+k22bttHb3Bo1dreGT1qmDdrKpdduJjLLljMjNLJSW65iJwUKhDMbC3wbSADuNfdv95veQ7wI+BC4BjwaXd/N1h2G3Aj0AP8pbs/GZT/APg4cMTdz0lIb0aTrrkPxcw496w5nHvWHOobW3h2yy6e+e1ODh/te4rq/kP1PPjYFh58bAtL5k/ngxcu4ZKViyiZnJ/ElovIkIFgZhnA3cBHgFrgZTPb5O7bY6rdCDS4+2IzWw/cCXzazCqB9cByYDbwKzM7y917gH8F/plIkEgamlpcwLUfuYDfv2Il22oO8tRLO9j8+t7oZasAu/cdYfe+I/zwZ7/h7IUzWXVuBavOXcCssuIktlxkYgpzhLAKqHH3dwDMbAOwDogNhHXA3wXTG4F/tsidSuuADe7eAew1s5pgey+5+/NmtiARnRgLqfhO5fHCzDhnyRzOWTKHtvZOtrz5Li++UsNrO2vp7e0FIs9N2vnOYXa+c5gfPfoS5TNLouGwaF6ZfuYiYyBMIMwBDsTM1wIXDVbH3bvNrBEoDcp/22/dOcNpoJndBNwEMG/evOGsKuNQXm42v/OBs/idD5zFieY2fvv6Xl7Yupsdew6d8gDBA4cbOHC4gZ8+9QrFRXmcf/ZcViwt5/ylc5lSpFNLIqMhTCAM9KdZ/xPqg9UJs+4Zufs9wD0AVVVV4+JEvv5YTYzJhXlceWklV15aSWNTG9Xb3mXLG+/y+q7aU04rNTa18Xz1bp6v3g3A/NmlrFxWzvlnz2XpwplkZ+naCJFECPM/qRYoj5mfCxwcpE6tmWUCxUB9yHVFKC7KY83qZaxZvYz2ji5e31XLljffZeu2fae9j2HfwWPsO3iMR55+jYyMSSyZP53li2ZTuXgWSytmkpuj+x1E4hEmEF4GlphZBfAekUHiz/Srswm4AXgJuA54xt3dzDYBD5rZN4kMKi8BtiSq8ZKecnOyuOi8Ci46rwJ3Z8/+Ol7bVcvrOw+wc+/70XEHgJ6e3ujYw0+fgklmLCwvY/niWSxdOIuzFkzXKSaRkIYMhGBM4BbgSSKXnf7A3beZ2e1AtbtvAu4DHggGjeuJhAZBvYeJDEB3AzcHVxhhZj8GfheYZma1wN+6+30J72GCxJ6r0gDn2DEzFs+fzuL507nuygtoa+/kzd0HeX3nAd7YVcvBusZT6ve6U7P/CDX7j/DoM68DUFZSxJIF0zlr/gyWzJ/OwvJpOs0kMgBLpWfaV1VVeXV1dVK+97Vf+Jfo9MZv/blCYZw43tQauRu65hDbag6ecpf0YCZNmkTFnFIWz5tOxdxSKuZMY97sqQoJSUtmttXdq8LU1f8ASWlTivK5ZMUiLlmxCIi83W3HO4fZsecQb+97nz37604ZoAbo7e1lz4E69hyoi5ZNMmPOzBIq5pSycG4ZFXNLWTBnmp7aKhOKAkHSSlFBLqvOXcCqcxcA0N3dw/5D9bz97hHe3vc+u999/7TTTBA51XTgUD0HDtVHr2YCKJ1SwLxZUymfOZXymSXMnVlC+cwS8nKzx6pLImNGgSBpLTMzg4XlZSwsL2PtB5cDkaOImv11vFNbx97aY+ytrTvl0Rqxjh1v4djxFl7dceCU8mklhZTPLIkGxezpU5hVVszkwlydTpSUpUAIof84i/7Dp7aiglxWLitn5bK+K6Lb2jt5971jvFN7NPp54HA9PT29A27jaEMzRxuaTwuK/NxsZpUVM2t6MbPLpjC7rDg6X5Cn008yvikQRIjcQb1s0SyWLZoVLevu7uHQ0RMcOFzPgUORO6drD9fz3pHGUy59jdXa3nna+MRJkwvzmFVWzIzSIqZPLWJ6aRHTp05memkRpcUFZGZmjFr/RMJQIIgMIjMzIzgtVAIr+spPBsX+Q/XUBo/YOFTXyKG6Rjo6uwbd3onmNk40t7Fr7+HTlhkwraSIsqmFTC+dHAmMIDRKpxRSWlxAVpYCQ0aXAkFkmE4JihjuTsOJ1mg4HKpr5OCR45Hpoyei76AeiAN1DU3UNTSxfc+hAetMLsyjdEoB06YUMrW4gNKSmOkpkS+9lU5GQoEgkiBmxtTiAqYWF7B88exTlrk7dQ3NHK5rpK6hifePNnGkPvg6doKGE61Dbv/kEcbe2qOD1inMz6F0SiFTi/MpLspn6uR8ppz8KspjyuR8SoryycvN0liYnEaBEMIpj75OYjskdZlZ9DTQQDq7uqlraKauvokjxyIhcaShmSPHTlDf2EJDY2v07XNn0tzaQXNrB/sOHjtjveysTEqCoCgpyqO4KJ8pk/MomRwJkskFuRQV5lJcmEdBXrbCY4JQIIiMA9lZmcyZPoU506cMuLynp5fjTa3Ry2CPHW+OfDa2cLShmfpgerDB7v46u7p5/9gJ3j828OW2sSZNmkRRQQ6TC/MoLsylqCCPyQW5TC7KjXwW5lFcmMfkwlyKCiJlGiBPTQoEkRSQkTEpMrg8pXDQOu7O8aY2jjU009DUyvETrTScaOX4iTYam1qpP9FX1v/u7TPp7e2lsamNxqY2DgxdHYg8oLAwP4fC/FwK87MpzMuhsCCXovwc8vNyKCo4uSyHwvwcCvJzKMrPITdHp7KSSYEQQgo97kkmMDOjZHL+kO+mdnfa2ruioXG8qY2GxpZIWDS10dTcTmNzG00tkc/2jsGvnBpMe0cX7R1dHG1oHtZ6kyZNioREXjaFBZHAyMvNpiAvm4LcbHKD6fzcbPKDz4K8bPLzcsjPzSYvJ4uMjEnDbq9EKBBEJhgzi/wyzcse9BRVrM6ubk40t9PU0s6JlnZONLVFPpvbIsHRHJQHn03NbcN7C1aM3t7e6OA5AzxiJIyc7Czyc7MoyMshL/rZFyR5uVnRz9ycbHKzM8kLwiSn3/REO1pRIISQDLxRAAAICElEQVSg9ynLRJadlcm0kkKmlQx+uiqWu9Pa3klTSwctrR00t0UGuptb2mlqDcpaO2hubY8Ogre0ddDU0kFnV/eI29vR2UVHZ1eoK7fOxIDsIFxycyJfeTlZ5OVkk5OTGZ3Ozc2KhErOyZCJfOVkZZKbk0l2VmZkPjuTnKzMcT2+okAQkYQyMwrycuJ6VEdXV080QFpaO2hqbae1rZO29i5a2jtoa+ukpb2TlrZO2oLP1vZOWts6aG3vor29M+6jk/6cvnBJpIyMSUFYBCGR3RcWuTl987nBZ3Z2JrnZmcyePoULKkf3vfIKBBEZN7KyMijJGnocZDAnx0daTwmNjkigtHVEy9o7uyIB0hH5hX9yur2ji7aOLtraO4c18D4cPT29tPZEgmw4LjqvQoEw7uiUkci4FTs+Mq1k6Ppn0tvbS1tMSLS3B5+dp063tZ+s00l7R3ckWDpPhk03nV3dkfLOLjo6uuI+ghmLFzgpEEREBjBp0qS4T30Nxt3p7u6NhENnNx1d3XR0dNHe2U1HZyRMOruC6SBUOoN6i8rLEtaOwSgQQoi97FQHCCISLzMjKyuDrKwMigqS3ZrT6YJdEREBFAgiIhJQIITgCbuQTURk/FIgiIgIoEAYNt2pLCLpKlQgmNlaM9tlZjVmdusAy3PM7KFg+WYzWxCz7LagfJeZXRV2myIiMraGDAQzywDuBq4GKoHrzayyX7UbgQZ3XwzcBdwZrFsJrAeWA2uB75hZRshtiojIGApzH8IqoMbd3wEwsw3AOmB7TJ11wN8F0xuBf7bIuZV1wAZ37wD2mllNsD1CbDMhnnhhG/sP1Y9oG73e99IRnTASkXQVJhDmwCnvxagFLhqsjrt3m1kjUBqU/7bfunOC6aG2mRBbt+/jle37R2PTIiJpJcwYwkB/FPe/DnOwOsMtP/2bm91kZtVmVl1XV3fGho6FykWzkt0EEZFREeYIoRYoj5mfCxwcpE6tmWUCxUD9EOsOtU0A3P0e4B6AqqqqYd8QcNVlyxP2hMD83GwuOq8iIdsSERlvwgTCy8ASM6sA3iMySPyZfnU2ATcALwHXAc+4u5vZJuBBM/smMBtYAmwhcoQw1DYTomr5/NHYrIhI2hkyEIIxgVuAJ4EM4Afuvs3Mbgeq3X0TcB/wQDBoXE/kFzxBvYeJDBZ3Aze7ew/AQNtMfPdERCQs8xR6g3xVVZVXV1cnuxkiIinDzLa6e1WYurpTWUREAAWCiIgEFAgiIgIoEEREJKBAEBERIMWuMjKzOmBfHKtOA44muDnjnfo8MajPE8NI+jzf3cvCVEypQIiXmVWHvewqXajPE4P6PDGMVZ91ykhERAAFgoiIBCZKINyT7AYkgfo8MajPE8OY9HlCjCGIiMjQJsoRgoiIDCHtA8HM1prZLjOrMbNbk92eRDGzcjN71sx2mNk2M/tCUD7VzJ4ys93BZ0lQbmb2/4KfwxtmdkFyexCf4J3cr5rZY8F8hZltDvr7kJllB+U5wXxNsHxBMts9EmY2xcw2mtnOYH9fnM772cz+Ovg3/ZaZ/djMctNxP5vZD8zsiJm9FVM27P1qZjcE9Xeb2Q0jaVNaB4KZZQB3A1cDlcD1ZlaZ3FYlTDfwJXdfBqwGbg76divwtLsvAZ4O5iHyM1gSfN0EfHfsm5wQXwB2xMzfCdwV9LcBuDEovxFocPfFwF1BvVT1beAJd18KnE+k/2m5n81sDvCXQJW7n0Pk8fjrSc/9/K/A2n5lw9qvZjYV+FsiryBeBfztyRCJi7un7RdwMfBkzPxtwG3Jbtco9fVR4CPALmBWUDYL2BVMfw+4PqZ+tF6qfBF5s97TwOXAY0RetHQUyOy/v4m8a+PiYDozqGfJ7kMcfZ4M7O3f9nTdz/S9n31qsN8eA65K1/0MLADeine/AtcD34spP6XecL/S+giBvn9cJ9UGZWklOExeCWwGZrj7IYDgc3pQLR1+Ft8CvgL0BvOlwHF37w7mY/sU7W+wvDGon2oWAnXAD4NTZfeaWQFpup/d/T3gG8B+4BCR/baV9N/PJw13vyZ0f6d7INgAZWl1WZWZFQI/Bf7K3U+cqeoAZSnzszCzjwNH3H1rbPEAVT3EslSSCVwAfNfdVwIt9J1GGEhK9zs43bEOqCDy2t0CIqdL+ku3/TyUwfqZ0P6neyDUAuUx83OBg0lqS8KZWRaRMPh3d/9ZUPy+mc0Kls8CjgTlqf6zuBS4xszeBTYQOW30LWCKmZ18FWxsn6L9DZYXE3m9a6qpBWrdfXMwv5FIQKTrfr4C2Ovude7eBfwMuIT0388nDXe/JnR/p3sgvAwsCa5QyCYyOLUpyW1KCDMzIu+y3uHu34xZtAk4eaXBDUTGFk6W/1FwtcJqoPHkoWkqcPfb3H2uuy8gsh+fcfc/BJ4Frguq9e/vyZ/DdUH9lPvL0d0PAwfM7OygaA2Rd5Sn5X4mcqpotZnlB//GT/Y3rfdzjOHu1yeBK82sJDi6ujIoi0+yB1XGYNDmo8DbwB7gq8luTwL7dRmRQ8M3gNeCr48SOX/6NLA7+Jwa1DciV1ztAd4kchVH0vsRZ99/F3gsmF4IbAFqgJ8AOUF5bjBfEyxfmOx2j6C/K4DqYF8/ApSk834G/h7YCbwFPADkpON+Bn5MZJyki8hf+jfGs1+BPw36XwP8yUjapDuVRUQESP9TRiIiEpICQUREAAWCiIgEFAgiIgIoEEREJKBAEBERQIEgIiIBBYKIiADw/wEE84Un9E6tagAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Pdf(suite)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's how we can compute the posterior mean" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "333.41989326371095" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def mean(suite):\n", " total = 0\n", " for hypo, prob in suite.Items():\n", " total += hypo * prob\n", " return total\n", "\n", "mean(suite)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or we can just use the method" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "333.41989326371095" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "suite.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sensitivity to the prior\n", "\n", "Here's a function that solves the train problem for different priors and data" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def MakePosterior(high, dataset, constructor=Train):\n", " \"\"\"Solves the train problem.\n", " \n", " high: int maximum number of trains\n", " dataset: sequence of observed train numbers\n", " constructor: function used to construct the Train object\n", " \n", " returns: Train object representing the posterior suite\n", " \"\"\"\n", " hypos = range(1, high+1)\n", " suite = constructor(hypos)\n", "\n", " for data in dataset:\n", " suite.Update(data)\n", "\n", " return suite" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's run it with the same dataset and several uniform priors" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "500 151.84958795903822\n", "1000 164.30558642273363\n", "2000 171.33818109150937\n" ] } ], "source": [ "dataset = [30, 60, 90]\n", "\n", "for high in [500, 1000, 2000]:\n", " suite = MakePosterior(high, dataset)\n", " print(high, suite.mean())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results are quite sensitive to the prior, even with several observations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Power law prior\n", "\n", "Now let's try it with a power law prior." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class Train2(Train):\n", "\n", " def __init__(self, hypos, alpha=1.0):\n", " Pmf.__init__(self)\n", " for hypo in hypos:\n", " self[hypo] = hypo**(-alpha)\n", " self.Normalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what a power law prior looks like, compared to a uniform prior" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X2QHPV95/H3Z2Z2Vw/oCWnBsh6QCHIcObZFWDB5MEns4IjYRuSCYygSQ4orXVyh4jsnueC7i3MhSVVcyR0XVygH/AiOMSbkwaoER+cYnLucDdbKYIEgmLXA0kpgCfSAQNLuzs73/uie3d7Z2Z1eaVcrTX9eVVPT/etf93QzYj77+/36QRGBmZlZabZ3wMzMzgwOBDMzAxwIZmaWciCYmRngQDAzs5QDwczMAAeCmZmlHAhmZgbkDARJGyU9I6lP0q1Nln9Y0lOSdkj6mqQLMstulPRs+roxU36JpCfSbX5ckqbnkMzM7GSo1ZXKksrAd4ErgX5gG3B9RDyVqfOzwKMRcUzSB4GfiYj3SzoX6AV6gAC2A5dExCFJ3wI+BDwCPAh8PCK+Mtm+LFu2LNasWXNyR2pmVlDbt29/KSK6W9Wr5NjWZUBfROwCkHQfsAkYCYSIeDhT/xHgV9Lpnwe+GhEH03W/CmyU9HVgYUR8My2/B7gGmDQQ1qxZQ29vb45dNjOzOknfz1MvT5fRCmBPZr4/LZvIzYz+sE+07op0Ou82zcxshuVpITTr22/azyTpV0i6h366xbpT2eZmYDPA6tWrW+2rmZmdpDwthH5gVWZ+JbCvsZKknwP+K3B1RAy0WLc/nZ50mwARcVdE9ERET3d3yy4wMzM7SXkCYRuwTtJaSZ3AdcCWbAVJFwN3koTB/syircC7JC2RtAR4F7A1Il4Ajkq6PD276APAl6fheMzM7CS17DKKiKqkW0h+3MvAZyJip6TbgN6I2AL8KXAO8Nfp2aO7I+LqiDgo6Q9JQgXgtvoAM/BB4HPAXJIxh0kHlM3MbGa1PO30TNLT0xM+y8jMbGokbY+Inlb18gwqn9UOvjrIgaMDVIdrLD2ni9ctnjPbu2RmdkZq+0D4xrMv8ffbk/Hqd29Yzi/2+OxWM7Nm2v5eRh3l0UMcrNZmcU/MzM5sbR8IlfLoJQ9Dww4EM7OJtH0gdFZGD7E6fPYMoJuZnW5tHwiVUiYQam4hmJlNpO0DoaMy2mXkMQQzs4m1fyCU3GVkZpZH+wdCZgzBg8pmZhNr+0AYe5aRWwhmZhNp+0DIXofgFoKZ2cQKEAi+DsHMLI8CBIIHlc3M8ihUILiFYGY2sbYPBA8qm5nl0/aB0OkWgplZLm0fCNkWgscQzMwmlisQJG2U9IykPkm3Nll+haRvS6pKujZT/rOSHs+8Tki6Jl32OUnPZZZtmL7DGlUuCaWZMFwLhmsOBTOzZlo+IEdSGbgDuBLoB7ZJ2hIRT2Wq7QZuAn47u25EPAxsSLdzLtAH/O9Mld+JiAdO5QBakUSlpJHxg+pwjXKpPJMfaWZ2VsrTQrgM6IuIXRExCNwHbMpWiIjnI2IHMFkn/bXAVyLi2Env7Ukae/sKtxDMzJrJEwgrgD2Z+f60bKquA77YUPbHknZIul1S10lsMxcPLJuZtZYnENSkbEp/ZktaDrwZ2Jop/gjwRuBS4FzgdydYd7OkXkm9Bw4cmMrHjvBT08zMWssTCP3Aqsz8SmDfFD/nl4G/i4ihekFEvBCJAeCzJF1T40TEXRHRExE93d3dU/zYxNiL09xlZGbWTJ5A2Aask7RWUidJ18+WKX7O9TR0F6WtBiQJuAZ4corbzG3s7SvcQjAza6ZlIEREFbiFpLvnaeD+iNgp6TZJVwNIulRSP/A+4E5JO+vrS1pD0sL4l4ZNf0HSE8ATwDLgj079cJrr8NXKZmYttTztFCAiHgQebCj7aGZ6G0lXUrN1n6fJIHREvGMqO3oqKh5UNjNrqe2vVAbfAtvMLI+CBEKmhVB1l5GZWTOFCwQPKpuZNVeIQMhehzDoQDAza6oQgZC9dYXveGpm1lwhAsG3rjAza60QgeBnIpiZtVaIQMgOKnsMwcysuUIEglsIZmatFSIQOjyGYGbWUkECwVcqm5m1UpBA8O2vzcxaKV4gVN1CMDNrpiCBkBlUrrmFYGbWTCECoeIWgplZS4UIBA8qm5m1VohA6Mzey8hdRmZmTeUKBEkbJT0jqU/SrU2WXyHp25Kqkq5tWDYs6fH0tSVTvlbSo5KelfSl9HnNM8JdRmZmrbUMBEll4A7gKmA9cL2k9Q3VdgM3Afc22cTxiNiQvq7OlH8MuD0i1gGHgJtPYv9zqZR8+2szs1bytBAuA/oiYldEDAL3AZuyFSLi+YjYAeT6tZUk4B3AA2nR3cA1ufd6inz7azOz1vIEwgpgT2a+Py3La46kXkmPSKr/6C8FDkdE9SS3OSUdJQ8qm5m1UslRR03KpvJn9uqI2CfpQuAhSU8Ar+TdpqTNwGaA1atXT+FjR7mFYGbWWp4WQj+wKjO/EtiX9wMiYl/6vgv4OnAx8BKwWFI9kCbcZkTcFRE9EdHT3d2d92PH8O2vzcxayxMI24B16VlBncB1wJYW6wAgaYmkrnR6GfCTwFMREcDDQP2MpBuBL0915/PqaLj9dfLxZmaW1TIQ0n7+W4CtwNPA/RGxU9Jtkq4GkHSppH7gfcCdknamq/8I0CvpOyQB8CcR8VS67HeBD0vqIxlT+PR0HliWpDFnGvlaBDOz8fKMIRARDwIPNpR9NDO9jaTbp3G9bwBvnmCbu0jOYDotOiolqoPDQHItQrYbyczMCnKlMuAWgplZC8UJhMw4wqCvVjYzG6cwgdBZ9qmnZmaTKUwgZK9F8MVpZmbjFSYQKr4FtpnZpAoTCH6uspnZ5AoaCG4hmJk1KlAgjL1a2czMxipQILiFYGY2mcIEggeVzcwmV5hAGNNCqLrLyMysUSEDoVpzC8HMrFFhAsG3rjAzm1xhAmHMrSt8czszs3EKEwhjbl3hFoKZ2TiFCYTs7a99pbKZ2XiFCQTf3M7MbHK5AkHSRknPSOqTdGuT5VdI+rakqqRrM+UbJH1T0k5JOyS9P7Psc5Kek/R4+towPYfUnK9UNjObXMtHaEoqA3cAVwL9wDZJWzLPRgbYDdwE/HbD6seAD0TEs5JeD2yXtDUiDqfLfyciHjjVg8jDVyqbmU0uzzOVLwP60mcgI+k+YBMwEggR8Xy6bMwvbUR8NzO9T9J+oBs4zGnmK5XNzCaXp8toBbAnM9+flk2JpMuATuB7meI/TruSbpfUNcF6myX1Suo9cODAVD92hG9/bWY2uTyBoCZlU/pFlbQc+DzwaxFR//P8I8AbgUuBc4HfbbZuRNwVET0R0dPd3T2Vjx3DXUZmZpPLEwj9wKrM/EpgX94PkLQQ+Efgv0XEI/XyiHghEgPAZ0m6pmaMB5XNzCaXJxC2AeskrZXUCVwHbMmz8bT+3wH3RMRfNyxbnr4LuAZ4cio7PlWVTAvBt64wMxuvZSBERBW4BdgKPA3cHxE7Jd0m6WoASZdK6gfeB9wpaWe6+i8DVwA3NTm99AuSngCeAJYBfzStR9ag011GZmaTynOWERHxIPBgQ9lHM9PbSLqSGtf7K+CvJtjmO6a0p6eoozLaZTTsexmZmY1TmCuVKyW3EMzMJlOYQMi2EDyGYGY2XnECoeTbX5uZTaYwgTDmSmW3EMzMxilMIJRLQmkm1MIDy2ZmjQoTCJLGPlfZA8tmZmMUJhBg7NXKgw4EM7MxChYI2cdousvIzCyrsIFQrbmFYGaWVahAGPtMBLcQzMyyChUIY7uM3EIwM8sqWCBkboHt007NzMYoViBU3EIwM5tIoQKhUvJpp2ZmEylUIGRbCH5qmpnZWMUKBN8C28xsQsUKhEr2tFMHgplZVq5AkLRR0jOS+iTd2mT5FZK+Lakq6dqGZTdKejZ93Zgpv0TSE+k2P54+W3lGjb2XkbuMzMyyWgaCpDJwB3AVsB64XtL6hmq7gZuAexvWPRf4feBtwGXA70taki7+BLAZWJe+Np70UeTU4ecqm5lNKE8L4TKgLyJ2RcQgcB+wKVshIp6PiB1A46/szwNfjYiDEXEI+CqwUdJyYGFEfDMiArgHuOZUD6aVDl+pbGY2oTyBsALYk5nvT8vymGjdFel0y21K2iypV1LvgQMHcn5scxW3EMzMJpQnEJr17ef983qidXNvMyLuioieiOjp7u7O+bHNjb2XkQPBzCwrTyD0A6sy8yuBfTm3P9G6/en0yWzzpHV6UNnMbEJ5AmEbsE7SWkmdwHXAlpzb3wq8S9KSdDD5XcDWiHgBOCrp8vTsog8AXz6J/Z+S7KDyoG9dYWY2RstAiIgqcAvJj/vTwP0RsVPSbZKuBpB0qaR+4H3AnZJ2puseBP6QJFS2AbelZQAfBD4F9AHfA74yrUfWxJyO0cM9Pjg80x9nZnZWqeSpFBEPAg82lH00M72NsV1A2XqfAT7TpLwX+NGp7OypWjivY2T6yPGh0/nRZmZnvEJdqbxw7mggvOJAMDMbo1CBsGhMIFRncU/MzM48hQqE+V3lkVtgHx8c9sCymVlGoQJBEgvmjg6buNvIzGxUoQIBPI5gZjaRwgWCxxHMzJorXCBkWwhHjrmFYGZWV7hAWORrEczMmipcICz0oLKZWVMFDAS3EMzMmil0ILziMQQzsxGFCwSfZWRm1lzhAiE7huAuIzOzUYULhLmd5ZFnKw9WawwM+TbYZmZQwECQ5IFlM7MmChcI0HAtggeWzcyAggbCQg8sm5mNkysQJG2U9IykPkm3NlneJelL6fJHJa1Jy2+Q9HjmVZO0IV329XSb9WXnTeeBTcYDy2Zm47UMBEll4A7gKmA9cL2k9Q3VbgYORcRFwO3AxwAi4gsRsSEiNgC/CjwfEY9n1ruhvjwi9k/D8eTiaxHMzMbL00K4DOiLiF0RMQjcB2xqqLMJuDudfgB4pyQ11Lke+OKp7Ox0WeRbYJuZjZMnEFYAezLz/WlZ0zoRUQWOAEsb6ryf8YHw2bS76PeaBAgAkjZL6pXUe+DAgRy725rPMjIzGy9PIDT7oY6p1JH0NuBYRDyZWX5DRLwZeHv6+tVmHx4Rd0VET0T0dHd359jd1nzHUzOz8fIEQj+wKjO/Etg3UR1JFWARcDCz/DoaWgcRsTd9PwrcS9I1dVr49hVmZuPlCYRtwDpJayV1kvy4b2moswW4MZ2+FngoIgJAUgl4H8nYA2lZRdKydLoDeA/wJKfJmFtgHxsi3VUzs0KrtKoQEVVJtwBbgTLwmYjYKek2oDcitgCfBj4vqY+kZXBdZhNXAP0RsStT1gVsTcOgDPwz8MlpOaIcujrKdHWUGBiqUa0FxweHmdfV8j+FmVlby/UrGBEPAg82lH00M32CpBXQbN2vA5c3lL0GXDLFfZ1Wi+Z2sH9oAEjGERwIZlZ0hbxSGfxsZTOzRgUOhMw4wgkPLJuZFTgQfLWymVlWYQPB1yKYmY1V3EDw7SvMzMYobiBkWggvHR2cxT0xMzszFDYQViyZOzK95+AxX5xmZoVX2EBYek4n8zrLABwbGObQa+42MrNiK2wgSGLluWNbCWZmRVbYQABYtXTeyPSelx0IZlZshQ6EsS2E47O4J2Zms6/QgbDq3NEWQr9bCGZWcIUOhNcvmUspfbTP/qMDDAwNz+4OmZnNokIHQmelxPmL5gAQAf3uNjKzAit0IACszg4s+0wjMyuwwgeCB5bNzBK5AkHSRknPSOqTdGuT5V2SvpQuf1TSmrR8jaTjkh5PX3+ZWecSSU+k63xckqbroKbCp56amSVaBoKkMnAHcBWwHrhe0vqGajcDhyLiIuB24GOZZd+LiA3p69cz5Z8ANgPr0tfGkz+Mk5ftMtp76LhvYWFmhZWnhXAZ0BcRuyJiELgP2NRQZxNwdzr9APDOyf7il7QcWBgR34zkF/ge4Jop7/00WDi3Y+RhOQNDNfa/MjAbu2FmNuvyBMIKYE9mvj8ta1onIqrAEWBpumytpMck/Yukt2fq97fY5mmTvR5ht7uNzKyg8gRCs7/0G/tVJqrzArA6Ii4GPgzcK2lhzm0mG5Y2S+qV1HvgwIEcuzt1K5eODiz71FMzK6o8gdAPrMrMrwT2TVRHUgVYBByMiIGIeBkgIrYD3wPekNZf2WKbpOvdFRE9EdHT3d2dY3enbnWmhfD9l16bkc8wMzvT5QmEbcA6SWsldQLXAVsa6mwBbkynrwUeioiQ1J0OSiPpQpLB410R8QJwVNLl6VjDB4AvT8PxnJS1580fmf7ui68yWK3N1q6Ymc2aloGQjgncAmwFngbuj4idkm6TdHVa7dPAUkl9JF1D9VNTrwB2SPoOyWDzr0fEwXTZB4FPAX0kLYevTNMxTdl5C+dw/qIuAAarNZ7e+8ps7YqZ2azR2XSaZU9PT/T29s7Ith/4Vj//tONFAH7qDcu46Yo1M/I5Zmanm6TtEdHTql7hr1Su23DB4pHp7+w+TK129gSlmdl0cCCkLuyez4I5yfUIR09Uee6AB5fNrFgcCKlSSbx19Wgr4fHdh2dxb8zMTj8HQka22+ix7zsQzKxYHAgZP7JiAR3l5Jq5Fw+f4MXDJ2Z5j8zMTh8HQkZXpcybVi4amXe3kZkViQOhwVtXjwZC73MHJ6lpZtZeHAgN3rp6MZX0QcvPHzjGsy8eneU9MjM7PRwIDRbO7eDyi5aOzNcvVjMza3cOhCZ+/i3nj0x/Z/cR9h3yHVDNrP05EJpYvnjumFNQtz7hVoKZtT8HwgQ2vuV1I9OP9h3k0GuDs7g3ZmYzz4EwgYvOP4eLzj8HgGot+OqTP5jlPTIzm1kOhElclWklPPzUfl+oZmZtzYEwibesXsQFy5KnqQ0NB/f86/OcTbcLNzObCgfCJCRx49vXkF6WwHdffJV/+beZea6zmdlscyC0sHrpPK566/KR+Qe+1c/BVz3AbGbtJ1cgSNoo6RlJfZJubbK8S9KX0uWPSlqTll8pabukJ9L3d2TW+Xq6zcfT13nTdVDT7d0blo88YvPEUI27//V5P0DHzNpOy0CQVAbuAK4C1gPXS1rfUO1m4FBEXATcDnwsLX8JeG9EvBm4Efh8w3o3RMSG9LX/FI5jRnVWStz09jUj8zv7X+H+R/fM3g6Zmc2APC2Ey4C+iNgVEYPAfcCmhjqbgLvT6QeAd0pSRDwWEfvS8p3AHEld07Hjp9u61y3gqreOnnX0zzv387WdPhXVzNpHnkBYAWT/HO5Py5rWiYgqcARY2lDnl4DHImIgU/bZtLvo9yRpSns+C/5dzwouWbtkZP6+R/bwuB+kY2ZtIk8gNPuhbuxAn7SOpDeRdCP9h8zyG9KupLenr19t+uHSZkm9knoPHJjdM3wkcfNPr+XC8+YDEAF/+bXv+TbZZtYW8gRCP7AqM78S2DdRHUkVYBFwMJ1fCfwd8IGI+F59hYjYm74fBe4l6ZoaJyLuioieiOjp7u7Oc0wzqrNS4pYrL6J7QdLzVa0Fdz60i4eeOmOHQMzMcskTCNuAdZLWSuoErgO2NNTZQjJoDHAt8FBEhKTFwD8CH4mI/1evLKkiaVk63QG8B3jy1A7l9Fk4t4Pf+oU3jJx5FAH3fmM3D3yrn2GffWRmZ6mWgZCOCdwCbAWeBu6PiJ2SbpN0dVrt08BSSX3Ah4H6qam3ABcBv9dwemkXsFXSDuBxYC/wyek8sJm2bEEXt773jSPdR5A8O+FP//EZXjo6MMmaZmZnJp1Nt2Lo6emJ3t7e2d6NMQaGhrnzoV3s2HNkpGxeZ5kbfmI1l/3QuZwFY+Vm1uYkbY+Inlb1fKXyKerqKHPLlRfxiz0rRm5xcWxwmE9+/Tn+x1e+y14/XMfMzhIOhGlQKol3b1jOre99I8sWdI6U/9u+o/zB3+7k3m/s5rCfp2BmZzh3GU2z44PDfHn7Xh56aj/Z8eVKWbz9h5dx1VuWc+45nRNvwMxsmuXtMnIgzJC9B49z7zd388wLR8eUlwQbLljMz64/jzcuX+AxBjObcQ6EM0BEsGPPEf7hsRd47sBr45afv6iLH79oKW/7oaV0Lzwr7+hhZmcBB8IZJCJ4au8rPPidF8e1GOouPG8+F1+wmIsvWMLrFs85zXtoZu3MgXCG2nvoOA8/tZ9v9r3MwFCtaZ3zF3XxphWLeNOKhfzw8gXM6Syf5r00s3biQDjDDQwN853dR3ik72We7D/CRBc4lwSrl81j3fkLWPe6c7iwez6L53tQ2szycyCcRY6eGGLH7iM8/v3D7Nz7CoPV5i2HuiXzO1jbPZ8Lls1n5blzuWDpPBbN6/AAtZk15UA4Sw1Wa/T94FV27j3C03uPsufgMfJ8RfO6yqxYMpfXL5nL6xfP4XWL5nD+ojksPafTQWFWcHkDoXI6dma2/dKH/nK2d+GkDSNO0MEJdXCCCoPqoNb0buPNiaAjhulgmArDdESNSjpdoUaJmMLWzGw2/c2f//qMbr8QgXA2KxPMZ5D5kVzpHAGDlBmgg0Gl75SpqflF54EYVIXB+lfd8OtfIqjEMGVqVEjCohzJdDnz8iXtZu3PgXCWEdDFMF0MjzyCKIDhKDFImUEqDKrMEBWGKDM8QVDU1dLAGPchDUpRo0xkAiKdj8x02uIou+VhdlYqRCDMdDPrTHZsoMqBowPsf2WAA0cHePnVQV46OsDBVwd5+dXBlgPYp6Kro8S8zjLzOivM7Swn011l5nQk03M6y8ztKDM3fe/qKDGno5y+SsztLNNZKXkMxOw0KUQgFNm8rgoXdFW4YNn8ccsiglcHqhx+bYhDxwY59OoQh48NcuTYEIePDfHKiSGOHBvi6PEq1ZN48M/AUI2BoRqHXhs66f2XkqfUdVWSsKhPd2WnK6VkOi3rrJToKCt9L42UdZZLdKTL6uUd5WS+XJKDxwrPgVBgklgwp4MFczpYtXTehPUigmODwxw9XuXI8SFePVFNXgPJ+9ETSdlrA8McG0jKjw0O5zo7qpWI0WB55Xj11Dc4AYmRcOgol6iURaWchEilLCqlseUdZVEpjS6rlEvpexIu4+ZLpfRdlEfKRuuWS6KsdJmy64lS+u7AspnmQLCWJDG/q8L8rkru22pEBCeGaryWhsPx9PXaQJUTQ8n0scFhTgwNc2KwxvGhZHpgqDayfKBam9EurbH7m5zyO1gFGD4tnzlVEpQ1GhDZ96SckWAppWFTUvKeTDMyn11nXJmSW7rXyyVGtpXdh/r2SpntZNcpabSOMsuV3Q8l/74EY7Ypja9f3y8JxNhtj5lPt2lTlysQJG0E/hwoA5+KiD9pWN4F3ANcArwMvD8ink+XfQS4meT/st+MiK15tmlnN0nJ2MAp3najVgsGqjUGhoY5MVRjcDgznQbGQHU4M11jqJrUq5cNVmsMDQdDaVl9eb1sqFqb8ErxM0kEVCOgFvjpGpNLQmNsQNQDg8yyUikJIzUEkQCahBNAKQ2gejkwGmyZ9ZUJO2WWJ/Nq2D/S7WjkvV63vh9zO8v80qUrZ/S/W8tAkFQG7gCuBPqBbZK2RMRTmWo3A4ci4iJJ1wEfA94vaT1wHfAm4PXAP0t6Q7pOq22aUSpNT7C0MlwLqsNJUFSHg+pwpNO1zHQwXEtCpDocVGv196Qs2cbY8uF0friWzNe3UX9Va7W0TmN5UKsFw5FsYzhiWrrgiiIiOfuuNvIf7ez/j7dgTmX2AwG4DOiLiF0Aku4DNgHZH+9NwH9Ppx8A/kJJm20TcF9EDADPSepLt0eObZqdNkm3SpmujjP3RoK1TFBUa0EtRgOkPl2rkYRIQ736dAQj9evT1VptTHktGAmjWm10vr5stN74ZWOm089oLKtvMxrrZbYZtUhOp87sazC6rXr9elnAaNnZ/9vfVKn+jN4ZlCcQVgB7MvP9wNsmqhMRVUlHgKVp+SMN665Ip1tt08wySiXReRp+FM520RgYaZo0hk80hFKybmYZo8E1tny0rB5aMDaM6uEWAWRaK/Xt1tclRtfPfh5ALR0+q9evlM+MQGi2F40ZPFGdicqbXS3VNNclbQY2A6xevXrivTQzI9OX70sjpyzPHQn6gVWZ+ZXAvonqSKoAi4CDk6ybZ5sARMRdEdETET3d3d05dtfMzE5GnkDYBqyTtFZSJ8kg8ZaGOluAG9Ppa4GHImn3bAGuk9QlaS2wDvhWzm2amdlp1LLLKB0TuAXYSnKK6GciYqek24DeiNgCfBr4fDpofJDkB5603v0kg8VV4DciYhig2Tan//DMzCwvPw/BzKzN5X0egu9qbGZmgAPBzMxSDgQzMwPOsjEESQeA709hlWXASzO0O2eqIh4zFPO4i3jMUMzjPtVjviAiWp63f1YFwlRJ6s0zkNJOinjMUMzjLuIxQzGP+3Qds7uMzMwMcCCYmVmq3QPhrtnegVlQxGOGYh53EY8Zinncp+WY23oMwczM8mv3FoKZmeXUloEgaaOkZyT1Sbp1tvdnpkhaJelhSU9L2inpQ2n5uZK+KunZ9H3JbO/rdJNUlvSYpH9I59dKejQ95i+lN01sK5IWS3pA0r+l3/mPt/t3Lek/pf+2n5T0RUlz2vG7lvQZSfslPZkpa/rdKvHx9Pdth6Qfm679aLtAyDzy8ypgPXB9+ijPdlQFfisifgS4HPiN9FhvBb4WEeuAr6Xz7eZDwNOZ+Y8Bt6fHfIjksa7t5s+Bf4qINwJvJTn+tv2uJa0AfhPoiYgfJbkRZv0Rve32XX8O2NhQNtF3exXJnaPXkTwr5hPTtRNtFwhkHvkZEYNA/fGcbSciXoiIb6fTR0l+IFaQHO/dabW7gWtmZw9nhqSVwLuBT6XzAt5B8vhWaM9jXghcQXJnYSJiMCIO0+bfNckdmeemz1mZB7xAG37XEfF/SO4UnTXRd7sJuCcSjwCLJS2fjv1ox0Bo9sjPFRPUbRuS1gAXA48C50fEC5CEBnDe7O3ZjPhfwH8G0ocMshQ4HBHVdL4dv/MLgQPAZ9Ousk9Jmk8bf9cRsRf4M2A3SRAcAbbT/t913UTf7Yz9xrVjIOR55GdbkXQO8DfAf4yIV2Z7f2aSpPcA+yNie7a4SdV2+84rwI8Bn4jgN8Z9AAABp0lEQVSIi4HXaKPuoWbSPvNNwFrg9cB8ku6SRu32XbcyY//e2zEQcj+esx1I6iAJgy9ExN+mxT+oNyHT9/2ztX8z4CeBqyU9T9Id+A6SFsPitFsB2vM77wf6I+LRdP4BkoBo5+/654DnIuJARAwBfwv8BO3/XddN9N3O2G9cOwZCYR7Pmfadfxp4OiL+Z2ZR9pGmNwJfPt37NlMi4iMRsTIi1pB8tw9FxA3AwySPb4U2O2aAiHgR2CPph9Oid5I8ibBtv2uSrqLLJc1L/63Xj7mtv+uMib7bLcAH0rONLgeO1LuWTlVbXpgm6RdI/mqsP57zj2d5l2aEpJ8C/i/wBKP96f+FZBzhfmA1yf9U74uIxgGrs56knwF+OyLeI+lCkhbDucBjwK9ExMBs7t90k7SBZCC9E9gF/BrJH3Vt+11L+gPg/SRn1D0G/HuS/vK2+q4lfRH4GZK7mv4A+H3g72ny3abh+BckZyUdA34tIqblUZJtGQhmZjZ17dhlZGZmJ8GBYGZmgAPBzMxSDgQzMwMcCGZmlnIgmJkZ4EAwM7OUA8HMzAD4/8ZIoA2ZH7x0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "high = 100\n", "hypos = range(1, high+1)\n", "suite1 = Train(hypos)\n", "suite2 = Train2(hypos)\n", "thinkplot.Pdf(suite1)\n", "thinkplot.Pdf(suite2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's see what the posteriors look like after observing one train." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8HNWV6PHf6U27N1nebSRjO8ZmsYPYMYGwBAjgTIJfTEIwE2bIxkuA4Q2QhEyGl2RgBkLChJAhQEKYTIBA8jAMGQIGAiSALQMBDBjLtrAFXrTY2lut7j7vjyrJ3a2W1JZV7m7pfD+f/qjq1q3bt9S2Tt9bt+4VVcUYY4wZab5sV8AYY8zoZAHGGGOMJyzAGGOM8YQFGGOMMZ6wAGOMMcYTFmCMMcZ4wgKMMcYYT1iAMcYY4wkLMMYYYzwRyHYFsmny5MlaWVmZ7WoYY0xeWb9+faOqVgyVb0wHmMrKSmpqarJdDWOMySsi8n4m+ayLzBhjjCcswBhjjPGEBRhjjDGeGNP3YIwxo0NPTw/19fWEw+FsV2VUKSwsZNasWQSDwWGdbwHGGJP36uvrKSsro7KyEhHJdnVGBVWlqamJ+vp6qqqqhlWGp11kInK2iGwUkVoRuS7N8QIRedA9/oqIVLrp5SLyrIi0i8hPEvKXicjrCa9GEfmRe+xSEWlIOPZ3Xl6bMSZ3hMNhysvLLbiMIBGhvLz8gFqFnrVgRMQP3AGcCdQD60Rktaq+nZDtMmCPqs4TkZXAzcBngTBwA3C4+wJAVduAJQnvsR74XUJ5D6rqFR5d0ojricXZ1dbNzPGF9h/DmANk/4dG3oH+Tr1swRwL1KrqFlWNAA8Ay1PyLAfuc7cfBk4XEVHVDlV9ESfQpCUi84EpwAsjX3Xv9cTi3PbsFu584X1+tbY+29UxxpgR52WAmQlsT9ivd9PS5lHVKNAClGdY/kU4LRZNSPuMiLwhIg+LyOzhVfvg+GBvmJauKADv7e4gEo1nuUbGmIOppqaGr3/96wB0d3dzxhlnsGTJEh588MEs12zkeHmTP13bSoeRZyArgS8k7D8G/EZVu0Xkyzgto4/3q5TI5cDlAHPmzMnwrUZeXHXQfWPM6FZdXU11dTUAr732Gj09Pbz++usZnx+LxfD7/V5Vb0R42YKpBxJbEbOADwfKIyIBYDzQPFTBInIUEFDV9b1pqtqkqt3u7s+Bo9Odq6p3qWq1qlZXVAw5lY4xxmSkrq6Oww/vu2XMLbfcwne/+11OPfVUrr32Wo499lgWLFjACy84vfrPPfcc5513Hrt37+biiy/m9ddfZ8mSJWzevJk1a9awdOlSjjjiCL74xS/S3e38aausrOTGG2/k5JNP5re//S2nnnoqV111FaeccgqHHXYY69at49Of/jTz58/n29/+dlZ+D4m8bMGsA+aLSBXwAU6L43MpeVYDq4CXgAuBZ1K6vAZyEfCbxAQRma6qO9zdC4B3DqDuxpg89Zkf/NGzsh/55lnDOi8ajbJ27VqeeOIJ/vmf/5mnn36679iUKVO4++67ueWWW3j88ccJh8OceuqprFmzhgULFnDJJZdw5513cuWVVwLOsykvvvgiAD/72c8IhUI8//zz/PjHP2b58uWsX7+eSZMmceihh3LVVVdRXp7pXYeR51kLxr2ncgXwJM4f+4dUdYOI3CgiF7jZ7gHKRaQWuBroG8osInXAD4FLRaReRBYlFP+/SAkwwNdFZIOI/BX4OnCpB5dljDH77dOf/jQARx99NHV1dYPm3bhxI1VVVSxYsACAVatW8fzzz/cd/+xnP5uU/4ILnD+nRxxxBIsXL2b69OkUFBQwd+5ctm/fTjZ5+qClqj4BPJGS9p2E7TCwYoBzKwcpd26atOuB64db12yzWzDG5LdAIEA8vm+wTuLzIwUFBQD4/X6i0eig5QzViVNSUpK031u2z+fr2+7dH+q9vGZP8ueIuuZODptWlu1qGJP3htuNdaCmTp3K7t27aWpqorS0lMcff5yzzz57v8tZuHAhdXV11NbWMm/ePO6//34+9rGPeVBj71mAyRFhG6ZsTF4LBoN85zvf4bjjjqOqqoqFCxcOq5zCwkJ+8YtfsGLFCqLRKMcccwxf/vKXR7i2B4dkdk99dKqurtZsLTi2pbGDe17a1z964dLpLJ01Pit1MSbfvfPOOxx22GHZrsaolO53KyLrVbV6qHNtun5jjDGesABjjDHGExZgckRHdyzbVTDGmBFlASZHvLurLdtVMMaYEWUBJkeUFNiAPmPM6GIBxhhjjCcswBhjzCjRO4FmrrAAY4wxeSrbU8EMxQJMjtiww27yG5Ov6urqWLhwIatWreLII4/kwgsvpLOzEyDt1Ptr167tmwDz0UcfpaioiEgkQjgcZu5cZ6rFzZs3c/bZZ3P00UezbNky3n33XQAuvfRSrr76ak477TSuvfbaAeu0du1aTjzxRJYuXcqJJ57Ixo0bATj33HN54403AFi6dCk33ngjADfccAN33333iP5e7M5yjhjDEyoYM6K+9di7npX9/fMHnv5l48aN3HPPPZx00kl88Ytf5Kc//SlXXHEFl156ab+p96+44gpee+01AF544QUOP/xw1q1bRzQa5bjjjgPg8ssv52c/+xnz58/nlVde4atf/SrPPPMMAO+99x5PP/30oAuOLVy4kOeff55AIMDTTz/NN7/5TR555BFOOeUUXnjhBSorKwkEAvz5z38G4MUXX+Tiiy8eqV8VYAEmZwR86Rb3NMbki9mzZ3PSSScBcPHFF3P77bdz5pln9pt6/4477uDKK69k3rx5vPPOO6xdu5arr76a559/nlgsxrJly2hvb+cvf/kLK1bsm2y+d9ExgBUrVgy5mmVLSwurVq1i06ZNiAg9PT0ALFu2jNtvv52qqio++clP8tRTT9HZ2UldXR0f+chHRvR3YgHGGGNGgIj02x9srsdly5bxhz/8gWAwyBlnnMGll15KLBbjlltuIR6PM2HChAGXUE6dsj+dG264gdNOO43f//731NXVceqppwJwzDHHUFNTw9y5cznzzDNpbGzk5z//OUcfnXYR4ANiAcYYM6oM1o3lpW3btvHSSy9xwgkn8Jvf/IaTTz550Kn3TznlFC655BIuueQSKioqaGpqYufOnSxevBgRoaqqit/+9resWLECVeWNN97gqKOOyrg+LS0tzJw5E4Bf/vKXfemhUIjZs2fz0EMPccMNN9DQ0MA111zDNddcM6K/D7Cb/DkjGrebMMbks8MOO4z77ruPI488kubmZr7yla8kTb1/xBFH4PP5+qbeP+6449i1axennHIKAEceeSRHHnlkX0vo17/+Nffccw9HHXUUixcv5tFHH92v+vzjP/4j119/PSeddBKxWPJUVMuWLWPq1KkUFxezbNky6uvrWbZs2Qj8FpLZdP05Ml0/wD+ds4BQwGK+Mfsr29P119XVcd555/HWW29lrQ5esen6R4lte7qyXQVjjBkxFmByyNhtSxqT3yorK0dl6+VAeRpgRORsEdkoIrUicl2a4wUi8qB7/BURqXTTy0XkWRFpF5GfpJzznFvm6+5rymBlGWPGhrHc3e+VA/2dehZgRMQP3AGcAywCLhKRRSnZLgP2qOo84DbgZjc9DNwADDSs4fOqusR97R6iLGPMKFdYWEhTU5MFmRGkqjQ1NVFYWDjsMrwcpnwsUKuqWwBE5AFgOfB2Qp7lwHfd7YeBn4iIqGoH8KKIzNuP9xuorLz5F9fQ3s38iqHHtxtjks2aNYv6+noaGhqyXZVRpbCwkFmzZg37fC8DzEwgcZhUPXDcQHlUNSoiLUA50DhE2b8QkRjwCPA9N4hkVJaIXA5cDjBnzpxhXJZ3drR0D53JGNNPMBikqqoq29UwKby8B5Nu7pPU1kQmeVJ9XlWPAJa5ry/sT1mqepeqVqtqdUVFxRBvdXCF/DZdjDFm9PAywNQDsxP2ZwEfDpRHRALAeKB5sEJV9QP3ZxvwXzhdccMqyxhjjHe8DDDrgPkiUiUiIWAlsDolz2pglbt9IfDMYPdMRCQgIpPd7SBwHtA7NnC/yjLGGOMtz+7BuPdBrgCeBPzAvaq6QURuBGpUdTVwD3C/iNTitDZW9p4vInXAOCAkIp8CzgLeB550g4sfeBr4uXvKgGXliw072zn/iGzXwhhjRoank12q6hPAEylp30nYDgMrUs9zj1UOUGzaKT8HKytfdEZye3U6Y4zZH/Ykfw4pLbDJrY0xo4cFGGOMMZ6wAGOMMcYTFmBySEtX1Ka6MMaMGhZgckxjRyTbVTDGmBFhASbHRGPWgjHGjA4WYIwxxnjCAkyOsfaLMWa0sACTY7Y0dma7CsYYMyIswOSYSCye7SoYY8yIsABjjDHGExZgjDHGeMICTI55e2dbtqtgjDEjwgJMjolE7R6MMWZ0sACTYwqD/mxXwRhjRoQFmBwQ9Eu2q2CMMSPOAkyO+WBvONtVMMaYEWEBJgf4fcktGLsPY4wZDSzA5IBpZQVJ++FoLEs1McaYkeNpgBGRs0Vko4jUish1aY4XiMiD7vFXRKTSTS8XkWdFpF1EfpKQv1hE/ltE3hWRDSJyU8KxS0WkQURed19/5+W1jSQRKCmwm/vGmNHFswAjIn7gDuAcYBFwkYgsSsl2GbBHVecBtwE3u+lh4AbgmjRF36KqC4GlwEkick7CsQdVdYn7unsEL8dziR+ErTlmjBkNvGzBHAvUquoWVY0ADwDLU/IsB+5ztx8GThcRUdUOVX0RJ9D0UdVOVX3W3Y4ArwKzPLyGgyaWEFW2NtmEl8aY/OdlgJkJbE/Yr3fT0uZR1SjQApRnUriITADOB9YkJH9GRN4QkYdFZPYA510uIjUiUtPQ0JDZlRwEnZF9N/ZtwktjzGjgZYBJ93BHaudPJnn6FywSAH4D3K6qW9zkx4BKVT0SeJp9LaPkwlXvUtVqVa2uqKgY6q0OmmMOmZDtKhhjzIjyMsDUA4mtiFnAhwPlcYPGeKA5g7LvAjap6o96E1S1SVW73d2fA0cPs95Zt7stku0qGGPMAfMywKwD5otIlYiEgJXA6pQ8q4FV7vaFwDOqg9/iFpHv4QSiK1PSpyfsXgC8cwB1P+h6ErrF3m+2ezDGmPwX8KpgVY2KyBXAk4AfuFdVN4jIjUCNqq4G7gHuF5FanJbLyt7zRaQOGAeERORTwFlAK/At4F3gVREB+Ik7YuzrInIBEHXLutSra/PChKJg33ZZgWcfizHGHDSe/iVT1SeAJ1LSvpOwHQZWDHBu5QDFpp24S1WvB64fVkVzwOyJRdmugjHGjCh7kj8Hbdzdke0qGGPMAbMAkyNKQslP8g9xK8oYY3KeBZgcMWtCYdK+hRdjTL6zAJMjRASxZWGMMaOIBZgc1dzRk+0qGGPMAbEAk0MSb7vsauseOKMxxuQBCzA55JBJNlTZGDN6WIDJIYkjyer3dmWxJsYYc+AswOSQ9si+lSztHowxJt9ZgMkhiUOVQwH7aIwx+c3+iuWQaWUFfdt1tuiYMSbPWYDJIYkPVzZ3WheZMSa/ZRRgROQREfmkiFhA8tCchAkvi0P2qzbG5LdM/4rdCXwO2CQiN4nIQg/rNGaVJkzTn7iEsjHG5KOMAoyqPq2qnwc+CtQBT4nIX0Tkb0UkOPjZZrjau6PZroIxxgxbxv0wIlKOs4jX3wGvAT/GCThPeVKzMagomPxxNHXY0snGmPyV0YJjIvI7YCFwP3C+qu5wDz0oIjVeVW6sERGmlIXY3eYElrhNqWyMyWOZrmh5t7s6ZR8RKVDVblWt9qBeY1ZhwvMvG3e1U1VenMXaGGPM8GXaRfa9NGkvjWRFjKMxsVvMpu83xuSxQQOMiEwTkaOBIhFZKiIfdV+nAkN+tRaRs0Vko4jUish1aY4XiMiD7vFXRKTSTS8XkWdFpF1EfpJyztEi8qZ7zu0izioqIjJJRJ4SkU3uz4kZ/xZyyAlVk/q2W7vsJr8xJn8N1YL5BHALMAv4IXCr+7oa+OZgJ4qIH7gDOAdYBFwkIotSsl0G7FHVecBtwM1uehi4AbgmTdF3ApcD893X2W76dcAaVZ0PrHH3847ft6/Z8tcPWrNYE2OMOTCDBhhVvU9VTwMuVdXTEl4XqOrvhij7WKBWVbeoagR4AFiekmc5cJ+7/TBwuoiIqnao6os4gaaPiEwHxqnqS+osWv8r4FNpyrovIT2vlCU8C1NRGspiTYwx5sAMepNfRC5W1f8EKkXk6tTjqvrDQU6fCWxP2K8Hjhsoj6pGRaQFKAcaBymzPqXMme721N7Rbaq6Q0SmDHBNl+O0gJgzZ84g1c+OxAkvG9ptmLIxJn8N1UVW4v4sBcrSvAaT7hZ16sDbTPIcSP7+mVXvUtVqVa2uqKjYn1MPCkm5wq6eWPqMxhiT4wZtwajqf7g//3kYZdcDsxP2ZwEfDpCnXkQCwHigeYgyZw1Q5i4Rme62XqYDu4dR56ybXJLcLdYZiVEU9A+Q2xhjctdQXWS3D3ZcVb8+yOF1wHwRqQI+AFbizGeWaDWwCmfI84XAM+69lYHeb4eItInI8cArwCXAv6eUdZP789HB6p6rRIQJxQH2djojyDq6Y5SXDHGSMcbkoKEetFw/3ILdeypXAE8CfuBeVd0gIjcCNaq6GrgHuF9EanFaLit7zxeROmAcEBKRTwFnqerbwFeAXwJFwB/cFziB5SERuQzYBqwYbt2zrTNhZcv3mzuZM6lokNzGGJObhuoiu2+w40Nxn/5/IiXtOwnbYQYIBKpaOUB6DXB4mvQm4PQDqG7OKA0FaI4668HYbDHGmHw11IOWP3J/PiYiq1NfB6eKY8/CaaV92y9t3ZPFmhhjzPAN1UV2v/vzFq8rYvaJJcxyWV5iqyEYY/LTUF1k692ffxKREM6MygpsdB+eNB44bFoZr9TtBWBrU1eWa2OMMcOT6XT9nwR+BmzGeRalSkS+pKp/GPxMMxwF/uSey2gsTsBvSygbY/JLptP13wqcpqq1ACJyKPDf7BvBZUbQrImFSfvRuBKwR2GMMXkm06/Fu3uDi2sLefogYz7wiRAK7Hukv7mzJ4u1McaY4RnqQctPu5sbROQJ4CGcezArcB6kNB6JRPfd6N/V1s2M8YWD5DbGmNwzVAvmfPdVCOwCPgacCjQAebneSr44JOHhyvebOrNYE2OMGZ6hRpH97cGqiElWkLB0ckvYFh4zxuSfTEeRFeIsDrYYpzUDgKp+0aN6jXmzJxbx3u4OALY121BlY0z+yfQm//3ANJwVLv+EM4txm1eVMjAlYbGxcDSexZoYY8zwZBpg5qnqDUCHOz/ZJ4EjvKuWqSwvznYVjDHmgGQaYHrHye4VkcNx1m2p9KRGBki+BwOwx4YqG2PyTKYB5i4RmQjcgLPuytvAzZ7VyhBMeXK/ob07SzUxxpjhyegmv6re7W7+CZjrXXVMotkTC9m+JwzA2zvaWDCldIgzjDEmd2TUghGRchH5dxF5VUTWi8iPRKTc68qNdYkLj9ltfmNMvsm0i+wBnKlhPoOztHEj8KBXlTKOo+dM6Ntev60lizUxxpj9l+lkl5NU9f8m7H/PXcbYeGhCUaYfjzHG5J5MWzDPishKEfG5r/+FM5uy8dCciclDlXti1lFmjMkfQy2Z3CYircCXgP8CIu7rAeCqoQoXkbNFZKOI1IrIdWmOF4jIg+7xV0SkMuHY9W76RhH5hJv2ERF5PeHVKiJXuse+KyIfJBw7N/NfQ24qK0ieo7+h3dZ4M8bkj6HmIisbbsEi4gfuAM4E6oF1IrJaVd9OyHYZsEdV54nISpyhz58VkUXASpypaWYAT4vIAlXdCCxJKP8D4PcJ5d2mqqNmeefURcZqGzpsVmVjTN7IeJlEEblARG5xX+dlcMqxQK2qbnGXV34AWJ6SZzlwn7v9MHC6iIib/oCqdqvqVqDWLS/R6cBmVX0/02vIR7JvWRgaO6wFY4zJH5kOU74J+AbOA5ZvA99w0wYzE9iesF/vpqXNo6pRoAUoz/DclcBvUtKuEJE3RORe98HQvHdC5b7L2PChTf9mjMkfmbZgzgXOVNV7VfVe4Gw3bTCSJk0zzDPouSISAi4Afptw/E7gUJwutB04yzz3r5TI5SJSIyI1DQ0NA9c+R0woDvZt26SXxph8knEXGTAhYXt8BvnrgdkJ+7OADwfKIyIBt9zmDM49B3hVVXf1JqjqLlWNqWoc+Dn9u9R6892lqtWqWl1RUZHBZWTX3JRJL7ujsQFyGmNMbsk0wPwL8JqI/FJE7gPWAz8Y4px1wHwRqXJbHCtx5jFLtBpY5W5fCDyjquqmr3RHmVUB84G1CeddREr3mIhMT9j9G+CtDK8tp00bV5C0v7fTFh8zxuSHIZ/kc2+6vwgcDxyD0311raruHOw8VY2KyBXAk4AfuFdVN4jIjUCNqq4G7gHuF5FanJbLSvfcDSLyEM79nijwNVWNufUpxhmZ9qWUt/xXEVmC05VWl+Z4XhJJ7i3c0tTJ1JSgY4wxuWjIAKOqKiL/T1WPpn8LZKhznwCeSEn7TsJ2GFgxwLnfB76fJr0TZyBAavoX9qdu+WR8UYCWLqflUrNtLydUjYrxC8aYUS7TLrKXReQYT2tiBnTIpKK+7Z2tNm2/MSY/ZBpgTsMJMpvdYcBvisgbXlbM7HPkjHFJ+85tKmOMyW2ZzqZ4jqe1MIOaV1GStB+OxikK+gfIbYwxuWHQACMihcCXgXnAm8A97gOR5iBKXd3y3Z3tLJ2dyUhxY4zJnqG6yO4DqnGCyzkM8PCi8V7At2802Yad9kS/MSb3DdVFtkhVjwAQkXtIfhbFHESzJxaytakLgHd2tme5NsYYM7ShWjA9vRvWNZZdx1Xa0GRjTH4ZqgVzlLseDDgPWBa5+4LziMy4gU81I+nQyck3+vd09jAxYZ4yY4zJNYO2YFTVr6rj3FeZqgYSti24HETFoeRRY6/Vt2SpJsYYk5n9mezSZFlJQpB53QKMMSbHWYDJuswfmlw8fd8Co00dPYPkNMaY7LMAkyWqsHlnC69ubuLt7XszOufImcm9kjZ1vzEml1mAyZL6xnb2tEeIq/Lu9j10dQ89SC9xTjKA93Z3eFU9Y4w5YBZgsqQ7ZXXKP6zfPkDOfXwpU/f/eUvziNbJGGNGkgWYHPHr5zZllG9+wrxk2/eEvaqOMcYcMAswOSSTWZKPq5yQtN8Tiw+Q0xhjsssCTA6pqW0cMs/CqaVJ++/usmljjDG5yQJMDrn190MvsZO6hPIz7w0dlIwxJhsswOSQnmgso26yj0zddx9md1vEyyoZY8yweRpgRORsEdkoIrUicl2a4wUi8qB7/BURqUw4dr2bvlFEPpGQXueuqPm6iNQkpE8SkadEZJP7My9nh1y3qWHIPGctrEja7+qx52GMMbnHswAjIn7gDpx1ZBYBF4nIopRslwF7VHUecBtws3vuImAlsBg4G/ipW16v01R1iapWJ6RdB6xR1fnAGnc/79z88OtD5plaVpC0v2ajdZMZY3KPly2YY4FaVd2iqhHgAWB5Sp7lOIuaATwMnC7OTYblwAOq2q2qW4Fat7zBJJZ1H/CpEbiGrIgOMTIs9T7MS1v3eFkdY4wZFi8DzEwg8enBejctbR53vZkWoHyIcxX4o4isF5HLE/JMVdUdblk7gCkjdB2eG1ccStr/z2eHfibmrMOSu8li8cznNDPGmIPBywAjadJS/woOlGewc09S1Y/idL19TURO2a9KiVwuIjUiUtPQMPT9joMhFEz+GB5b+/6Q5xyf8jxMzbbM5jMzxpiDxcsAUw/MTtifBXw4UB4RCQDjgebBzlXV3p+7gd+zr+tsl4hMd8uaDuxOVylVvUtVq1W1uqKiIl2WrPjquYuT9rfubB0gp6MgkLw+zOo3d414nYwx5kB4GWDWAfNFpEpEQjg37Ven5FkNrHK3LwSeUWec7mpgpTvKrAqYD6wVkRIRKQMQkRLgLOCtNGWtAh716Lo88fGjZiTtX3Pvy0Oes3RW8uzK8QyGOBtjzMHiWYBx76lcATwJvAM8pKobRORGEbnAzXYPUC4itcDVuCO/VHUD8BDwNvA/wNdUNQZMBV4Ukb8Ca4H/VtX/ccu6CThTRDYBZ7r7eUNEmDEpeVnktq7Bn3E5Z3HybaZ171s3mTEmd0gmD/aNVtXV1VpTUzN0Rg/88c0d3Pp0LQBTS4P86u+Pp6s7ysW3PtOXx+/z8dB1Zwxazrceezdp//vnLxz5yhpjTAIRWZ/ymEha9iR/DikqCCTtx+Jxuod4iLJ6zviUc8buFwZjTG6xAJNj7vzqsqT9L9z67KD5zzt8atL+/7yddmyDMcYcdBZgcsyUCcmrVsbicTrCPQPmD/qTP8K/2EOXxpgcYQEmB/00pRVzyQ8Hb8VccERyK2ZXa/eI18kYY/aXBZgcNDWlFQPw/u62AfMfe0jyQ5e3/2nriNfJGGP2lwWYHPWrq09L2r/67pcGzCsiTChOHiBgK10aY7LNAkyOKikMMnty8uqV//bIXwfM/6WTDkna/8919Z7UyxhjMmUBJofd9vcnJO2/vHEXe9rT318ZVxhM2q9t6Mxo8TJjjPGKBZgcJiJce+GSpLS/u/1PA+b/uxPnJO0/+U5uTOZpjBmbLMDkuGMX9F914PJ/fz5t3qry4qT9FzY3e1InY4zJhAWYPPDb685M2m9qC/Pnd3amzbti6fSk/T/Yg5fGmCyxAJMHfD7h5kuPS0r74e/foLE13C/vklnJU8e8uLnZ7sUYY7LCAkyemDdjPMfMT16/5ks/eZ5ImrnKPrMkuRVzz0vbPK2bMcakYwEmj1y3Ymm/tIv+bQ3xlAkuPzo7uRWztanLnosxxhx0FmDyzMPXn9kvbcVNT/ULMqkjyr77xHue1ssYY1JZgMkzIsJv/s/p/dJX3PRU0r2W1BFlAJsaOjytmzHGJLIAk4dCQT/3XXVav/QL/+UpeqL7usK+/Yn5Scd/+fJ2W1bZGHPQWIDJU6VFQe79xqn90lf+69M0tzlP+xeF/Bw+vSzp+A2PbzwY1TPGGAsw+Wx8SajfpJgAf//vf+p7Tuai6pn9jtds2+t53YwxxtMAIyJni8jxCAC+AAAayElEQVRGEakVkevSHC8QkQfd46+ISGXCsevd9I0i8gk3bbaIPCsi74jIBhH5RkL+74rIByLyuvs618tryxUlhcG092R++Ps3WHHTUwD8nzMOTTr2+7/upCMSPSj1M8aMXZ4FGBHxA3cA5wCLgItEZFFKtsuAPao6D7gNuNk9dxGwElgMnA381C0vCvyDqh4GHA98LaXM21R1ift6wqtryzWhoD/t6LJ4XPnMD/6IH+W4yuQ1Y37wZK09gGmM8ZSXLZhjgVpV3aKqEeABYHlKnuXAfe72w8DpIiJu+gOq2q2qW4Fa4FhV3aGqrwKoahvwDtC/D2gMEhEe+eZZzJ02rt+xS297jtfe/KBf+rftfowxxkNeBpiZwPaE/Xr6B4O+PKoaBVqA8kzOdbvTlgKvJCRfISJviMi9IjLxwC8h//zbF4/nxs9X90t/o66Jd97cRldK19i3Hnv3YFXNGDPGeBlgJE1aap/MQHkGPVdESoFHgCtVtdVNvhM4FFgC7ABuTVspkctFpEZEahoaRud09osPmcRD157RL12ArRs/pKa2Ial77EfPbjmItTPGjBVeBph6YHbC/izgw4HyiEgAGA80D3auiARxgsuvVfV3vRlUdZeqxlQ1Dvwcp4uuH1W9S1WrVbW6oqIiXZZRwe/38cg3z+KK8w5PSvcBRd3drN/cSE1tA6A0tEf4yfNbs1JPY8zo5WWAWQfMF5EqEQnh3LRfnZJnNbDK3b4QeEadr9argZXuKLMqYD6w1r0/cw/wjqr+MLEgEUmc4fFvgLdG/Iry0GlHzug3AMAfj1MQiQBQU+sEmg/3hvnuE3ZPxhgzcgJeFayqURG5AngS8AP3quoGEbkRqFHV1TjB4n4RqcVpuax0z90gIg8Bb+OMHPuaqsZE5GTgC8CbIvK6+1bfdEeM/auILMHpSqsDvuTVteWb3gEA7V09rLrtWQCCsRja00Mk6Cy1vH5zIwD/0NXDrZ85fMCyjDEmUzKWh6pWV1drTU1NVt77j2/u4NanawGYWhrkV39//EF77+2N7Vx5118A6PH76Q6F+uU5fU4Z//iZJf3SjTFGRNarav/RRCnsSf4xaPbkUh755lnc/qWTCMZiFHZ398uzZlsbn/qXp/jMD/5I7YctWailMSbfedZFZnLfzPISHvnmWXR1R/ncD5+ls7Aw6XhXYSGiyrW/fKUv/1XLj6AqzbM2xhiTygKMoaggwO+vP5OeWJwv3LuOPe2RvmMqQntREcXhMB80dXDNvS8DEAr4+cq5izh50TR8vnSjyo0xY50FGNMn6PfxwN8fx7cee5dYLM6mHS20h50HM3tbNyVdXQgQicb48eo3+fHqNwFYOncyF582n8qpZQMVb4wZYyzAmH6+f/5Cahs6+MXLzmQKkZ4YW3a10h6O0lFURCAWozASSTrntS2NvLalsW//o4dOZvnxlSyeMxFndLkxZqyxAGPSmldRwvfO+wjffnwjoaCfhbOcmXeisTjbG9tpavMTjEYp6OlJe/6rmxt5dfO+gFNcEODc6jmcduQMpk4osqBjzBhgAcYMSET4/vkLaWjv5kfPOk/6B/w+qqaOo2oqgLKnPcLupna0dfDlmDu7ozz85y08/Od909KUFAb55DFzqJ5XwSFTSgn4bVCjMaOJBRgzpIrSAr5//kJe3NzEH95OnL9NmFhawMTSAqAcVeWU2WX89yvv82Hz4AEHoCPcw0MvbOahFzb3pfl8wlGV5SycPYHDD5lE5ZRSCkP2z9SYfGT/c03GTj60nJMPLeexN3fycl3/VTFFhBfq2xk3s5wLTl3AafPLqdvVxh/Wb+e5N1OnoUsvHtd+93N6jS8JcdjsiSyaPZHKKaXMmVJKWVH/h0SNMbnBAozZb+cfMY3zj5jGn7c088SG3WnzPLepiec2NQFwxkem8ZVzFxHw+9zRaa385Z2dvFu/l807WtOen05LR4SX393Fy+/u6nesMBRgZnkxlVPKOGRKGXMqSpk2sYhJpQX4revNmKywAGOG7aS5kzhp7iSaOyPcumbgKf+f3tjI0xudFsmCKSWcs2gKXzxzYVKetq4I733QwqubG3mrrpn6pqG72BKFI1E272gdNGCNLwkxZXwRUyYUMW1CEdMnlTB1QhGTxxUyqazA7gEZM8IswJgDNqk4xPfPdwLGs+/tCybpvLe7g/d271sa4KOzx7Ps0ElUlIY4el4FR89LXkIhGotTt6uNLTvbePeDvby7fQ+79nYNq54tHRFaOiJsGmTqm4DfRzQWp7QoyOLZEykI+SkrCjKzvITyskJKi4JMLC1gfHGQgqDfRsMZMwgLMGZEnbZgMqctmEwsrry0tTllUEB/r25v4dXt+/7gV5SGOKFqIofPKKMkFCDg9zFvxnjmzRjPWR+dlXSuqtLc3s2HTZ1sa2hnW0M7m3e0Ut/UQU80Nqz6R2NxANq7enjlvfTdf+kcOn0cHzR1sqSqnK5IlGkTi5lUVkDQ72NGeQkBnzCprIDigkDfy4KTGe0swBhP+H3SNyhAVdm4u4M/vL2bxvbIoOc1tEdY/eYuVr+57z7L7ImFLJhSyoIpJUxL6MoSEcrLCikvK+SIyklpy+uJxtnd0sWuvV3s3tvFjuZOdu3tZOeeLrY3to/Y9fZ2zb280an3X7c2ZXxuYShAwC+UFQYZVxIi6PcxeZwzc8KcilI6wj3MqSglFlfKxxUSCvgoDDlBqijkpygUsOl6TE6yAGM8JyIsnFrKwqmlAESicWq272XNxkbCPfEhz9++J8z2PWHWpHS9TRtXwOwJhcyaWMTsCUVMKgkSTLmPEgz4mFlewszykkHfIxaL0xbuobmtm4aWMM1tYdrDUep2t9HUGsbnEzZ92EI8PvLLW4QjznQ87V097NjTOawy/D4fsXicwlCACSUhdu7pZNGciXT3xOgIR/nIrAk0toSZUV7M+OIQrZ0Rpk0spqQwQDgSY/L4QgoCfhQYVxQkGPARDPgoCgUIBXwE/c6+tbrM/rAAYw66UMDHiVWTOLHKaXXE4sqmhg5eqdvDe7szv7m/s7Wbna3drNvW/55KWWGAitIQU8sKmFpWwOTSEBOLg5QVBPCn+bbv9/uYUFLAhJIC5mYwW7Sq0t0To6Wzh85wDw2tYbq6o+ztiLBjTycFAT9v1DVTMb6Qt95vxu8TOrujGV/b/orFnUAdjkTZ6Qast7ft6Tu+0w1cG7Y1H/B7lRQG6Qg7MzjMnuw8ILt1VyvzZ4ynrMhZwK6+sYPDKycR9PuIq9LW2cOcKaUE/T6i8TiCUO62Rv1+IRqNM6G0gIBPCPh9iDiTsAZ8PvxuWijgI+B3tv0+sWCXByzAmKzz+5JbOODcC3l/TxebGzqpbezgg73h/SqzLRylLRxlS+PALYJQQJhWVoDfJ0wdV8jkkhAFAR/jiwKMKwxQEgpQGPThS/OHTEQoDAXch0CL9nsJg2gsTmd3lM7uKG1dPYQjMVo6uolE4zS2hYlG40RjypZdrVSML+K1zY3MLC9hW0M7XZEY5WUF7NzTSSjgJzLM+03D1RtcgKRuxtTBE8/89YOk/d7uw5FWWhQk4HOCT2NrmOKCAFMmFDnByS80tIQJuS3ZgE/wuwHs/V1tLJozEX9fEBMaWsNMGV/kBjfpC2572ruZNrEYv0/wuel+nxCJxikrCvaV25uuCoUhPz5x8/r3neOTfUFytAdKCzAmJwX8Pg6dXMKhk0s4i30jyyLRODtaw9Q1d7G1sZMPWrrojAzdzZZOJKps2+MErq1NQ49MKwg438bnlhfTEo4yd3Ix3dE408cVEosrk0uCiAhlhQGCfqEo6Kcg4COQ5o9IwO9jXHGIccUhpk0cVvX7qCo90ThdkZjbJdZDT8zZb3eDwY7mTopCASLRGLU7WqkYV0hPLM6bdc3MnlxKLB5n3aYGDp0+jngctu5qxSdCaVGQ1s7B75tlW3tX8nx4nd1R6na19cu3M0334/4Oh/eC4AQtn0/6BqeMKw71BaPeQNTRHaW1M0Ll1LK+Y76EoNXSEaGjO8pc97gkHHPKjrNrbyfzpo8n6PexbPE05s0Y7+m1WYAxeSUU8HHIpGIOmVTMx+aVJx2LxZXmzgiN7RF2tXWzuy3CjtYwje0RRuLWSXfUCWQb3W68na29K4EOveJnYcBHOBpnfFGAaFzpicWpnFTM7vZuZowrpKQgQFNHhEMmFRH0+WiPRJnqtq6iMWV8kdO1J0BxyO98s/b53C4lIeATxhUHEQkBRQd+sWnEYnF63FdXd4xoPE5PNN73Bz4ad7oN97ZHKAj66InFae3sob2rh+LCAD3ROB82dyLiTH7aE43T3NbN7pYuZpaXuEEyyqYPWzhkShmxWJxoXPsCQ3FBgJ6YDnuEYK5SlFhciSV8TxosqKcLnoma2wZv7fcOSFnz1w+4++uneDoVk6cBRkTOBn4M+IG7VfWmlOMFwK+Ao4Em4LOqWuceux64DIgBX1fVJwcrU0SqgAeAScCrwBdUNbe/epkR5fcJFaUFVJQWcNi09OvSxOJKazjK3q4eOiMxGtsj7G7vJhZX9nb1UL83jE9I+s8+EsJucGrp2ncfpvd+097Ofd1Mg3XpZcrvg8Kgn4BPaOmKEvQLk0tDBNx7T9v3hPnI1BL87jfb3mufPaHI/cbrdAHu7exh6jgnyPnFmSeuJ6YEfEJx0I/Ph/PtOBTEL840P+WlhVSUl+IT95j7EyDgd79NJ6T7fOBD9pUlgggIDNp1pKpEonGiMacrMeYGN+cPdZxYXIm6rbiY++2iN280FqepLUxpYbAvXyyubGto7xu915seicap3eEGvLgSiznl90TjvL61iUWzJxJXJa7q1iPO5h2tTCorpDDoJ669wUOJRGO0d/VQEPQTV+f31XvvLBu6IlEaW7uZNTkPA4yI+IE7gDOBemCdiKxW1bcTsl0G7FHVeSKyErgZ+KyILAJWAouBGcDTIrLAPWegMm8GblPVB0TkZ27Zd3p1fSY/+X3CxOIgE4uDGeWPxuJ09cTp7InRFo4SjStNHRGiMaU13ENPXNnT2YNPYEdLN2WFAXa2djOhOMDeTu9u6g8mFoeO7n3f8ntiyo6W7qQ8G3f17xpKzQPw1o7Bvy17yQl0TuDpiTlBojjkjGTrDVJOMHK2AZo6epg2riApSPUGTOkrU/CJnzd2djBjXCEFQR8+vzBp2gTUDXABIChQjHDM1PF9QbJ+bxdVk4oRYPFh0/vKTXy/j7OvXuLWEaCxPUJFaSghiDonCU7AdIKOE3DiqqDipiuqEI3F8PmE9u4YRQEfqrjH1GkVR+ME/b595wAaVzcYK53dPZQUBnjqtQ/oCPfgPwjBzcsWzLFArapuARCRB4DlQGKAWQ58191+GPiJOF9blgMPqGo3sFVEat3ySFemiLwDfBz4nJvnPrfcEQ8wm3a08vM/DTwtSqb2dqVfR8XkloDfR5nfR1lhgKllBcMqIxZXuqNxuqNxIrE4Xe437dZwFL8IkVicxvYIRSE/XT0xdrd1M7EoSHskxs7WbipKQ0Rjzki76eMLUIWG9m5icSgJ+emJx4lER374dLbFFVCIse/aMrnftq/rcmh1zfs/K8TutjzuGIlE6QgECRf48MXjqHr778bLADMT2J6wXw8cN1AeVY2KSAtQ7qa/nHLuTHc7XZnlwF5VjabJP6J2t4b5axa/1Zn84/cJxSE/xSG/p+/T+0026nbJ9MSUaDxOPA4xt6umIxLD73aLxeJKJBanMxKjIOBz0lQJ98RpDUcpK3C6cnrP3d3WTUHA3zfYIa5KPA4dkSi72rqZMb6QuDr1cMqC1q4e2rpjVJSG3HNI+tnb0goFnJFXsbj7zXv0xcucE/f52NMVZbaH7+FlgEnXgZr6z2agPAOlp5uNcLD8/SslcjlwOcCcOXPSZTno5lWUDp3JmCGICEG/EPQ2jh0UmhCE1A1ymrDf2wWUuB+NKSJOcEo85uxrX/BTha6emNOd5LyZm989Tuo5Tlmt4SilBf6+4KcKcZw3U/aV7Ta8kvZ3tXZTURZKOldJeC/3fUk4P7Fevb+T+r3hvsEffecklLdve19Zvef2br/1fuKzUPnbgqmHpOA4C0hdFKQ3T72IBIDxQPMQ56ZLbwQmiEjAbcWkey8AVPUu4C6A6urq/f7tVlaUcP7hU/f3tAHNmFDEeUtmjFh5xowGIs7AAn/a747mQMTDEVrDTmfP9AnejDjs5WWAWQfMd0d3fYBz0/5zKXlWA6uAl4ALgWdUVUVkNfBfIvJDnJv884G1OC2VfmW65zzrlvGAW+ajXlzUzEklXHHmgqEzGmNMDrpk2dyD9l6eBRj3nsoVwJM4Q4rvVdUNInIjUKOqq4F7gPvdm/jNOAEDN99DOAMCosDXVDUGkK5M9y2vBR4Qke8Br7llG2OMyRLxehRBLquurtaamppsV8MYY/KKiKxX1eqh8tkSfsYYYzxhAcYYY4wnLMAYY4zxhAUYY4wxnrAAY4wxxhNjehSZiDQA7w/z9Mk4D3iOJXbNY4Nd89hwINd8iKpWDJVpTAeYAyEiNZkM0xtN7JrHBrvmseFgXLN1kRljjPGEBRhjjDGesAAzfHdluwJZYNc8Ntg1jw2eX7PdgzHGGOMJa8EYY4zxhAWYYRCRs0Vko4jUish12a7PSBCR2SLyrIi8IyIbROQbbvokEXlKRDa5Pye66SIit7u/gzdE5KPZvYLhExG/iLwmIo+7+1Ui8op7zQ+KSMhNL3D3a93jldms93CJyAQReVhE3nU/7xNG++csIle5/67fEpHfiEjhaPucReReEdktIm8lpO335yoiq9z8m0Rk1YHUyQLMfhIRP3AHcA6wCLhIRBZlt1YjIgr8g6oeBhwPfM29ruuANao6H1jj7oNz/fPd1+XAnQe/yiPmG8A7Cfs3A7e517wHuMxNvwzYo6rzgNvcfPnox8D/qOpC4Cicax+1n7OIzAS+DlSr6uE4S32sZPR9zr8Ezk5J26/PVUQmAf+EsxT9scA/9QalYXGW9bRXpi/gBODJhP3rgeuzXS8PrvNR4ExgIzDdTZsObHS3/wO4KCF/X758euGsfroG+DjwOM6ido1AIPXzxlmH6AR3O+Dmk2xfw35e7zhga2q9R/PnDMwEtgOT3M/tceATo/FzBiqBt4b7uQIXAf+RkJ6Ub39f1oLZf73/WHvVu2mjhtslsBR4BZiqqjsA3J9T3Gyj5ffwI+Afgbi7Xw7sVWfpbUi+rr5rdo+3uPnzyVygAfiF2y14t4iUMIo/Z1X9ALgF2AbswPnc1jO6P+de+/u5jujnbQFm/6VbJHzUDMUTkVLgEeBKVW0dLGuatLz6PYjIecBuVV2fmJwmq2ZwLF8EgI8Cd6rqUqCDfd0m6eT9NbtdPMuBKpwl2EtwuohSjabPeSgDXeOIXrsFmP1XD8xO2J8FfJiluowoEQniBJdfq+rv3ORdIjLdPT4d2O2mj4bfw0nABSJSBzyA0032I2CCiPQuJ554XX3X7B4fj7PUdz6pB+pV9RV3/2GcgDOaP+czgK2q2qCqPcDvgBMZ3Z9zr/39XEf087YAs//WAfPdESghnJuFq7NcpwMmIgLcA7yjqj9MOLQa6B1Jsgrn3kxv+iXuaJTjgZbepni+UNXrVXWWqlbifI7PqOrngWeBC91sqdfc+7u40M2fV99sVXUnsF1EPuImnQ68zSj+nHG6xo4XkWL333nvNY/azznB/n6uTwJnichEt+V3lps2PNm+KZWPL+Bc4D1gM/CtbNdnhK7pZJym8BvA6+7rXJy+5zXAJvfnJDe/4Iym2wy8iTNCJ+vXcQDXfyrwuLs9F1gL1AK/BQrc9EJ3v9Y9Pjfb9R7mtS4BatzP+v8BE0f75wz8M/Au8BZwP1Aw2j5n4Dc495h6cFoilw3ncwW+6F57LfC3B1Ine5LfGGOMJ6yLzBhjjCcswBhjjPGEBRhjjDGesABjjDHGExZgjDHGeMICjBlzRERF5NaE/WtE5LsjVPYvReTCoXMe8PuscGdCfjYlvVJEPjfMMv8yMrUzxmEBxoxF3cCnRWRytiuSyJ2pO1OXAV9V1dNS0iuBtAEm4an1tFT1xP14f2OGZAHGjEVRnOVir0o9kNoCEZF29+epIvInEXlIRN4TkZtE5PMislZE3hSRQxOKOUNEXnDzneee7xeRfxORde76G19KKPdZEfkvnAfeUutzkVv+WyJys5v2HZwHY38mIv+WcspNwDIReV2cNVAuFZHfishjwB9FpFRE1ojIq265ywe41udk35oxv3afgMe97rfda7hlv3/zZkwZ9BuNMaPYHcAbIvKv+3HOUcBhOPNSbQHuVtVjxVmc7X8DV7r5KoGPAYcCz4rIPOASnOk4jhGRAuDPIvJHN/+xwOGqujXxzURkBs5aJEfjrFfyRxH5lKreKCIfB65R1ZqUOl7npvcGtktxpqI/UlWb3VbM36hqq9uCe1lEVmv/J66XAotx5qH6M3CSiLwN/A2wUFVVRCbsx+/OjEHWgjFjkjozRf8KZyGqTK1T1R2q2o0zxUZvgHgTJ6j0ekhV46q6CScQLcSZ0+kSEXkdZxmEcpzFngDWpgYX1zHAc+pM0hgFfg2csh/17fWUqvZO1ijAD0TkDeBpnKnYp6Y5Z62q1qtqHGfaoEqgFQgDd4vIp4HOYdTFjCEWYMxY9iOcexklCWlR3P8XbrdQKOFYd8J2PGE/TnJvQGproHca9P+tqkvcV5Wq9gaojgHql27q9OFILP/zQAVwtKouAXbhzL2VKvFaYzgLc0VxWluPAJ8C/meE6mdGKQswZsxyv9U/xL6lcgHqcLqkwFlDJDiMoleIiM+9LzMXZ7XAJ4GvuEsiICILxFnoazCvAB8TkcnuAICLgD8NcU4bUDbI8fE4a+D0iMhpwCEZXA9unUuB8ar6BE534JJMzzVjk92DMWPdrcAVCfs/Bx4VkbU4s88O1LoYzEacQDAV+LKqhkXkbpxuplfdllEDTitgQKq6Q0Sux5lWXoAnVPXRwc7BmSE5KiJ/xVmjfU/K8V8Dj4lIDU7X17v7cV1lOL+bQrc+/QZJGJPIZlM2xhjjCesiM8YY4wkLMMYYYzxhAcYYY4wnLMAYY4zxhAUYY4wxnrAAY4wxxhMWYIwxxnjCAowxxhhP/H+K7En5IjS6FwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dataset = [60]\n", "high = 1000\n", "\n", "thinkplot.PrePlot(num=2)\n", "\n", "constructors = [Train, Train2]\n", "labels = ['uniform', 'power law']\n", "\n", "for constructor, label in zip(constructors, labels):\n", " suite = MakePosterior(high, dataset, constructor)\n", " suite.label = label\n", " thinkplot.Pmf(suite)\n", "\n", "thinkplot.Config(xlabel='Number of trains',\n", " ylabel='Probability')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The power law gives less prior probability to high values, which yields lower posterior means, and less sensitivity to the upper bound." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "500 130.70846986256004\n", "1000 133.2752313750312\n", "2000 133.99746308073065\n" ] } ], "source": [ "dataset = [30, 60, 90]\n", "\n", "for high in [500, 1000, 2000]:\n", " suite = MakePosterior(high, dataset, Train2)\n", " print(high, suite.mean())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Credible intervals\n", "\n", "To compute credible intervals, we can use the `Percentile` method on the posterior." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(69, 869)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hypos = range(1, 1001)\n", "suite = Train(hypos)\n", "suite.Update(60)\n", "\n", "suite.Percentile(5), suite.Percentile(95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have to compute more than a few percentiles, it is more efficient to compute a CDF.\n", "\n", "Also, a CDF can be a better way to visualize distributions." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYFNXZ9/HvLSKisoiAyOag4gIEN1xxX3F5XEFBTUQR1MSVqI/GxO194pK4oTEG3BUFAY2gIaIiLlGUxV0UBRVlUVAQkSAyzP3+UTVj90xPT83Q1TU9/ftc11zT59Sp4q4p4J46p+occ3dEREQA1ks6ABERqT+UFEREpIKSgoiIVFBSEBGRCkoKIiJSQUlBREQqKCmIiEgFJQUREamgpCAiIhXWTzqA2mrdurWXlJQkHYaISEGZOXPmt+7epqZ2BZcUSkpKmDFjRtJhiIgUFDObF6Wduo9ERKSCkoKIiFRQUhARkQpKCiIiUkFJQUREKsSWFMzsfjNbbGYfVLPdzOwOM5tjZu+Z2S5xxSIiItHEeafwINAny/YjgK7h1xDg7hhjEREpWB98uoA33v2M9z9ZwE+r18T6Z8X2noK7v2JmJVmaHAs87MF6oG+YWUsz28LdF8UVk4hIIVi6fCX/fuUDnnzh7Srbhv3hZDpuvmlsf3aSL691AL5KKc8P66okBTMbQnA3QefOnfMSnIhIPrg7Xyz4jieef5up78ytsX2zjTaMNZ4kk4JlqPNMDd19BDACoFevXhnbiIgUAnfn3dnzGTfpLT76rHYdI8cetCMtmjWNKbJAkklhPtAppdwRWJhQLCIisSgtXcsb733O2GdnMv+bZbXe/6j9f8XRB/SkbatmMURXVZJJYQJwnpmNBvYAlms8QUQK3c9rSnntrbmMmjiN775fWat9WzRrSr/Dd+XA3bdjwyaNY4owu9iSgpmNAg4AWpvZfOBqoDGAu/8DmAgcCcwB/gucEVcsIiJxWbNmLVOmzWbcczNrnQS27tSGEw7dmd16lNCoUf14bSzOp48G1LDdgd/F9eeLiMRhzZq1/OetOYx8+k2+X/HfWu37q207cHKfXmy/VTvMMg2rJq/gps4WEcmnsrIypr77OY9PnM6Cxd/Xat89e3ah7+G7UtJhs3qbBCpTUhARSeHufPDpQh4aP5XP539bq3333nlr+h2+K523aBVTdPFTUhCRojdv4VJG/Wsa0z/4olb77dmzCycdsRtbti/cJFCZkoKIFJ3lK1Yx7rmZTHwl49Rs1dq125b0P7IXXTq2LpjuoNpSUhCRBq+0dC0vTf+E+598ndU/R587qOuWbfn1MXvSbestGmwSqExJQUQaHHdn7pdLeGj8VGbNjf76U7vWzRlw5O7stdNW9eYR0XxTUhCRBmHlqtX88/m3+efkd2q136+P2ZM++3RP7GWx+kZJQUQKkrvz/icLuGfsqyxcsjzyfvvvti2nHLU7rTfdJMboCpeSgogUjJWrVjP22Zk8/dJ7kffp0rE1px+7Fz26ti+acYF1oaQgIvWWuzP782+4Z9x/+GJB9HcGTj16D47avwdNNlCXUG0pKYhIvVJaupanX3qPkU+/GXmfHl3bc+YJ+zSo9wWSoqQgIon7dtmP3P/ka7z53ueR9xl43N4csW931l+/UYyRFR8lBRFJxKy5i7hz5IssXroiUvvu27Rn0Im92bL9ZjFHVtyUFEQkL9ydF9/8mL+PejnyPn0P24XjD9lZj4vmkZKCiMRm9c9reOK5t3ni+bcitW+28Yacf+qB7NKts54USoiSgojk1A8/ruKh8W/w0rTZkdr33LYjQ07aly3atIg5MolCSUFE1tmSpSsYMfZV3pr1ZaT2R+7Xg1OP3kPdQvWQkoKI1MmiJcv522NT+PizryO1H3Ribw7v3b1o5xQqFEoKIhLZgsXfM+zhycz9akmNbTdu2oTzTj2Q3XpsqfGBAqKkICJZLVqynNseeiFSImjbqhkX/vpgtt+qXR4ikzgoKYhIFd8u+5HbHn4hUtdQp3abcvHph+j9gQZCSUFEAFix8ifueuylSEtSdt6iFZeceRgd2raMPS7JLyUFkSK2+uc1PPDP13n+9Y9qbNup3aZccuZhdNx80zxEJklRUhApMu7OhCnv8fD4qTW2bde6OZcNOlxdQ0VESUGkSMz4cB43jPh3je022agJfxhyBNt10WBxMVJSEGnA5n65hMtueSJS28sH92G3HiXxBiT1npKCSAPz85pSBlxyb6S2Z/Xdhz77dNd7BFJBSUGkgRgx5lUmvfZhje0O792dgcfvxQaN9c9fqtLfCpECNm/hdwy9aWyktvdc92tatdg45oik0CkpiBQYd+fMPz7MDz+uqrHtBacdxP67bZuHqKShiDUpmFkfYBjQCLjX3W+stL0z8BDQMmxzubtPjDMmkUI17f0vuOneZ2ts13XLttxw8fEaJ5A6iS0pmFkj4C7gUGA+MN3MJrj7rJRmfwTGuPvdZtYNmAiUxBWTSKEpKyuj38UjIrVV95DkQpx3CrsDc9z9MwAzGw0cC6QmBQeah59bAAtjjEekYEz/4AtuvKfmu4JTj96DEw7dOQ8RSbGIMyl0AL5KKc8H9qjU5hrgOTM7H9gYOCTGeETqNXfnjCsfYsXKn2psO/rmwTRu3CgPUUmxiTMpZOrQ9ErlAcCD7n6Lme0FPGJmPdy9LO1AZkOAIQCdO3eOJViRpHz19TIuuuHxGttd9OuD2bdX1zxEJMUszqQwH+iUUu5I1e6hQUAfAHefamYbAq2BxamN3H0EMAKgV69elROLSEGK+l7B47cMZv31dVcg+RFnUpgOdDWzLsACoD9wSqU2XwIHAw+a2Q7AhkDNK3mIFKjS0rWc/Pt7amw36MTeHLnfr/IQkUi62JKCu5ea2XnAJILHTe939w/N7DpghrtPAH4P3GNmFxN0LQ10d90JSIMT9SWzh288g42bNslDRCKZ1ZgUzGwG8ADwmLsvq83Bw3cOJlaquyrl8yygd22OKVJIxk6ayeiJ07O22a1HCZcP7pOniESyi3Kn0B84g+A9g/IE8Zx+oxep3m+ve4xvvvsha5vrzj+G7tu0z1NEItHUmBTcfQ5wpZn9CTgauB8oM7P7gWHuvjTmGEUKwpo1a+l/Sc3jBaNuPkuT0Um9Felvppn1JLhbOBJ4AngU2Ad4EdgptuhECsC3y37k7GtGZm2zXZd2XH/RcXmKSKTuoowpzAS+B+4jmJtodbjpTTPTeIAUrY8/+5orhz2Vtc2Qfvty+D7d8xSRyLqLcqfQr3yqinJm1sXdP3f3E2KKS6Teenn6J9wx8sWsbW6/4mQ6tdMC91J4oiSFccAuGep2zX04IvXXuOfeYtS/pmVtM/KmM2m64QZ5ikgk96pNCma2PdAdaGFmqXcEzQleMhMpCvc98R8mvvJB1jbjbj9bU1VLg5DtTmE7gqeNWgL/k1K/AhgcZ1Ai9cGdj07hpWmzq93eetNNGH7NaXmMSCR+1SYFdx8PjDezvdx9ah5jEklUTclg5x068cdzjspjRCL5k6376DJ3/wtwipkNqLzd3S+INTKRPBs+5hWee21WtdsP2WsHzu2/fx4jEsm/bN1HH4XfZ+QjEJGkjJo4nXGTZla7/ZgDd+T04/bKY0QiycnWffR0+P2h/IUjkj+T3/iIv496udrtxx+8E6cds2ceIxJJXrbuo6epuihOBXc/JpaIRGL2/icLuOaup6vdfvT+PTnjhL3zGJFI/ZGt++jmvEUhkgfffPcDv73usWq3995lG4aerhVhpbhl6z6q/r5apID8tHoNp152X7Xbu27ZlhuH6uV8EcjefTTG3U8ys/dJ70YywN29Z+zRiawDd6fvRcOr3W7AuGHn5C8gkQKQrfvowvD70fkIRCSXbhjxb2Z8OK/a7WNvG8J6662Xx4hECkO27qNF4fd5ZtYO2J3gjmG6u3+dp/hEauX1d+ZyywPPV7v9kRvPZKOmmptIpDpRps4+C7iKYO0EA+40s+vc/f64gxOJavmKVZz5x+qfnr71f/uxZfvN8hiRSGGKMkvqpcDO7v4dgJltBrxOsAKbSKJqGjcYdGJvjtzvV3mMSKSwRUkK8wkmwSu3AvgqnnBEohs54Q3+OfmdjNs6br4pw/5wcp4jEil82Z4+Ghp+XECwytp4gjGFY4Hsk8qLxGjx0hWce+2j1W4fc+sQGjXSILJIXWS7U2gWfp8bfpUbH184ItWrqato2B9OpuPmWu1MZF1ke/ro2nwGIpLNs69+yD3jXs24TdNSiOROlKeP2gCXEazCVrHimrsfFGNcIgCs/nkNp1xa/dvIWvFMJLeiDDQ/CjxO8BLbOcDpwJI4gxIB+PPwibw168uM2+64sj8d2rbMc0QiDV+UpLCZu99nZheG8yG9bGaaF0lis3T5SgZf9UjGbbv1KOHywX3yHJFI8YiSFNaE3xeZ2VHAQqBjfCFJMTv7mpF8u+zHjNtG3zyYxo0b5TkikeISJSn8n5m1AH4P3Ak0By6ONSopOtkeMz23//4cstcOeY5IpDjVmBTc/Znw43LgwHjDkWJ0zjWPsmTZiozbNJAskl9Rnj7aChgG7AWUAVOBi939s5hjkwbu+xX/ZdAfH8647U/nHsVO23fKc0QiEuW1z8eAMUA7oD0wFhgV5eBm1sfMZpvZHDO7vJo2J5nZLDP70MyqXxZLGpS/3jep2oQw7vazlRBEEhJlTMHcPfVRkJFmdl6NO5k1Au4CDiWYP2m6mU1w91kpbboCVwC93X2ZmbWtXfhSaH5eU8qAS+7NuO3Ks49kl26d8xyRiKTKNvdRq/DjlPC3/NEEcx+dDPwrwrF3B+aUdzOZ2WiCeZNmpbQZDNzl7ssA3H1xrc9ACsaUN2fzt8emZNymsQOR+iHbncJMgiRQ/i/17JRtDvy/Go7dgfTZVOcDe1Rqsy2Amb0GNAKucfdnKx/IzIYAQwA6d9ZvkoUm25xFg/vuS599u+c5IhGpTra5j7qs47Ez/drnlcrrA12BAwjefXjVzHq4+/eVYhkBjADo1atX5WNIPZbtRbRRN5/FBo2j9GCKSL5EefqoMXAusF9Y9RIw3N3XVLtTYD6QOlrYkeDFt8pt3giP9bmZzSZIEtNrDl3qu0effpMnX3i7Sn2Pru259rxjEohIRGoS5de0u4HGwN/D8q/DurNq2G860NXMuhCsydAfOKVSm6eAAcCDZtaaoDtJj7oWuGzdRZqzSKR+i5IUdnP3HVPKL5rZuzXt5O6l4VNKkwjGC+539w/N7DpghrtPCLcdZmazgLXApeXLfkphyrZWsgaTReq/KElhrZlt7e5zoeJltrVRDu7uE4GJlequSvnswNDwSwrc5Dc+4u+jqs6V+D8H9GTg8VrvQKQQREkKlxI8lvoZweDxlsAZsUYlBaffxSMoKyurUn/31afStlWzDHuISH2UNSmY2XrAKoLB3+0IksLH7r46D7FJASgrK6PfxSMyblN3kUjhyTrNhbuXAbe4+2p3f8/d31VCkHLfLvsxY0LouW1Hnhh2jhKCSAGK0n30nJmdCDwZjgGI8OqMT7n9kclV6q87/xi6b9M+gYhEJBeiJIWhwMZAqZn9RNCF5O7ePNbIpN760x3jmTV3UZX6kTedSdMNN0ggIhHJlSjrKWiUUCqceOE/MtZr/ECkYah2TMHM2prZ7Wb2jJldb2a6Myhia9eWZUwIJR1aa/xApAHJNtD8MLCSYAnOZsAdeYlI6p2Vq1Zz0tCqA8rn9t+fWy7rm0BEIhKXbN1H7dz9yvDzJDN7Kx8BSf2yYPH3XPDn0VXqb7/iZDq12zSBiEQkTtmSgpnZpvwy22mj1LK7L407OEnWW7O+5M/DJ1ap14CySMOVLSm0IFhTIbWzuPxuwYGt4gpKkjf+xXd5ePzUKvVjbxvCeutFWcVVRApRtvUUSvIYh9Qjwx6ZzCszPq1S/8SwcxKIRkTySSucSJrf/2UcXyz4tkq9EoJIcVBSkAqZJrXruW1Hrv7d0QlFJCL5pqQgQOaX0o4/eCdOO2bPBKIRkaRESgpmtg/Q1d0fMLM2wCbu/nm8oUm+ZEoIZ5+0H4f17pZANCKSpChrNF8N9CKYOvsBgqU5RwK94w1N8iFTQrj0zMPYc0c9XCZSjKLcKRwP7Ez4OKq7LzQzzYfUAGRKCFf/9mh6btcxgWhEpD6IkhR+dnc3Mwcws41jjknyIFNC+POFx7H9Vu0SiEZE6osobyGNMbPhQEszGwy8ANwTb1gSp0wJ4fqLlBBEJNrU2Teb2aHADwTjCle5+/OxRyax6HfR8Cp1f77wOLbrooQgItEGmi8GxioRFL6LbxxDWaXF86497390hyAiFaJ0HzUnmCX1VTP7nZltHndQknt3jHyRLxelz2F4+eA+9OjaIaGIRKQ+qjEpuPu17t4d+B3QHnjZzF6IPTLJmacmv8PL0z9Jq/vdgAPYrUdJMgGJSL1Vm+kuFwNfA98BbeMJR3Jt5ofzeGTCG2l1/frsykF7bp9QRCJSn9WYFMzsXDN7CZgMtAYGu3vPuAOTdbdg8fdcP+LfaXV79OxC/yN2SygiEanvorynsCVwkbu/E3cwkjsrV62usmLaJhs14bJBhycUkYgUgmqTgpk1d/cfgL+E5Vap27XyWv1VWrqW31z+QJX6h244I4FoRKSQZLtTeAw4mmD1NSd9BTatvFZPuTsn/77qu4VaD0FEosi28trR4fcu+QtH1lXfDC+njb1tSAKRiEghijLQPDlKXTX79jGz2WY2x8wuz9Kur5m5mfWKclzJ7O7RL1epe+TGM7WmsohElm1MYUNgI6C1mW3KL91HzQneV8jKzBoBdwGHAvOB6WY2wd1nVWrXDLgAeLNOZyAAvDt7Pi9M/Sit7tb/7cdGTTdIKCIRKUTZfoU8m2A8Yfvwe/nXeIL/7GuyOzDH3T9z95+B0cCxGdr9P4LB7J9qEbekWLHyJ677+zNpdWf13Yct22+WUEQiUqiqTQruPiwcT7jE3bdy9y7h147u/rcIx+4AfJVSnh/WVTCznYFO7p7+P5pE5u4M/MODaXXtWjfniH17JBOQiBS0KLOk3mlmPYBuwIYp9Q/XsKtlqKuYjc3M1gNuAwbWFIOZDQGGAHTu3Lmm5kUl08DyXX86JYFIRKQhiDLQfDVwZ/h1IEFXzzERjj0f6JRS7ggsTCk3A3oAL5nZF8CewIRMg83uPsLde7l7rzZt2kT4o4vDfU/8p0rdmFv1pJGI1F2Ux1L6AgcDX7v7GcCOQJMI+00HuppZFzPbAOgPTCjf6O7L3b21u5e4ewnwBnCMu8+o7UkUo3kLlzLxlQ/S6kZcexqNGulJIxGpuyj/g6xy9zKg1MyaE0yMV+OLa+5eCpwHTAI+Asa4+4dmdp2ZRbnTkGqUlZUx9KYxaXWD++7LZi03SSgiEWkoosx9NMPMWhIswTkT+BGYFuXg7j4RmFip7qpq2h4Q5ZgC/S4eUaWuz77dE4hERBqaKAPNvw0//sPMngWau/t78YYl1XnoqalV6sbdfnYCkYhIQ5Tt5bVdsm1z97fiCUmq8+2yH5kw5d20ugevH4hZpge9RERqL9udwi1ZtjlwUI5jkSzcnbOvGZlWd/ZJ+9Fs4w2r2UNEpPayTYh3YD4DkezOuPKhKnWH9e6WQCQi0pDVOKZgZr/JVB/h5TXJkdmff82KlemzgGjmUxGJQ5Snj1LXbtyQ4J2FtwAlhTxwd/5w+1Npdbdc1lczn4pILKI8fXR+atnMWgCPxBaRpKk8jcVWndpQ0qF1QtGISENXl183/wt0zXUgUtX7nyyoUvfXS05MIBIRKRZRxhSe5peJ7NYjmBhvTPV7SC64O9fc9XRa3d+v0kR3IhKvKGMKN6d8LgXmufv8mOKR0OlXPJhW3mGrLdh8s+bJBCMiRSPKmMLLAOG8R+uHn1u5+9KYYytaCxZ/z8pVq9Pq/u/CTOsTiYjkVpTuoyEEq6OtAsoI1klwIkyKJ3VzwZ9Hp5XvuLJ/QpGISLGJ0n10KdDd3b+NOxiB4WNeSStv3LQJHdq2TCgaESk2UZ4+mkvwxJHEbPXPa3jutVlpdQ/dMDCZYESkKEW5U7gCeN3M3gQqOrrd/YLYoipSp1x6X1p56MBDNdmdiORVlKQwHHgReJ9gTEFi8Om8b6rU9d556wQiEZFiFiUplLr70NgjKXKX3/rPtPKD1w9MJhARKWpRxhSmmNkQM9vCzFqVf8UeWREZ82z6stTbdG6rKbFFJBFR7hTKX6O9IqVOj6TmSFlZGY//Oz0p3Dj0+ISiEZFiF+XltS75CKRYDb1pbFr57JP20+CyiCRG6ykk6KfVa/jq62VpdVo4R0SSpPUUEnTqZemPoN409ISEIhERCWg9hYR89/2PVeq22bJtApGIiPxC6ykkZMjVI9PK9/1fxl46EZG80noKCViw+PsqdS2bbZRAJCIi6bSeQgIqz4I66uazEopERCRdtUnBzLYBNi9fTyGlfl8za+Luc2OPrgGat/C7tHLzTZqyQeMouVlEJH7ZxhRuB1ZkqF8VbpM6qPxewohrTksoEhGRqrIlhRJ3f69ypbvPAEpii6gBm7cwfbG69m1a0Lhxo4SiERGpKltSyDb5TtNcB1IMht6UPj5/6/+elFAkIiKZZUsK081scOVKMxsEzIwvpIZpydL0nrg2mzbTXYKI1DvZRjgvAv5pZqfySxLoBWwARJqxzcz6AMOARsC97n5jpe1DgbMInmpaApzp7vNqdQYF4pxrH00r36l1l0WkHqo2Kbj7N8DeZnYg0COs/pe7vxjlwGbWCLgLOBSYT3DnMcHdU9ebfBvo5e7/NbNzgb8AJ9fhPOq1FSt/Sisb6C5BROqlKNNcTAGm1OHYuwNz3P0zADMbDRwLVCSF8Njl3gAa5KM451d6L+GhG89IKBIRkezqMs1FVB2Ar1LK88O66gwC/p1pQ7jIzwwzm7FkyZIchhi/0tK1Ve4UNm7aJKFoRESyizMpZFoUwDPUYWanEYxX/DXTdncf4e693L1XmzZtchhi/O54NP0mS2MJIlKfxfkq7XygU0q5I7CwciMzOwS4Etjf3VfHGE8iXntrTlq5fduWCUUiIlKzOO8UpgNdzayLmW0A9AcmpDYws52B4cAx7r44xlgS8fL0T9LKV559ZEKRiIhEE1tScPdS4DxgEvARMMbdPzSz68zsmLDZX4FNgLFm9o6ZTajmcAXpjpHpD2rt0q1zQpGIiEQT60xs7j4RmFip7qqUz4fE+ecn6etvf0grH7X/rxKKREQkuji7j4raRTc8nlYeeNxeCUUiIhKdkkIMysrKWFO6tqLcqNF6rLeeftQiUv/pf6oYPPjU1LTyXX8ckFAkIiK1o6QQg3+9/H5auU2rZglFIiJSO0oKObZoyfK08sDj9k4oEhGR2lNSyLGLbkxfM+HoA/TUkYgUDiWFHHJ3SlMGmJtv0hSzTLN9iIjUT0oKOfTStPQ3mG/6/QkJRSIiUjdKCjn0t8fSJ79rqwFmESkwSgo58vOa0rRyn326JxSJiEjdKSnkyPAxr6aVf3PsnglFIiJSd0oKOfLStNlp5SYbNE4oEhGRulNSyIGly1emlQed2DuhSERE1o2SQg4Mf/yVtPIR+/ZIKBIRkXWjpJADMz6cl1bWuwkiUqiUFNZR5a6jS888LKFIRETWnZLCOqrcdbRHzy4JRSIisu6UFNaRuo5EpCFRUlgHK1etTitfNujwhCIREckNJYV1cO+4/6SVd/9VSTKBiIjkiJLCOnhlxqdpZXUdiUihU1Koo9U/r0krD+m3b0KRiIjkjpJCHT1TacnNg/fcPqFIRERyR0mhjh57Zlpaef31GyUUiYhI7igp1IG7p5WP3E/TWohIw6CkUAez5i5KK/c9bNeEIhERyS0lhTp48KmpaeUWzZomFImISG4pKdTBZ18tqfi8dac2CUYiIpJbSgq19NPq9EdRBx6/d0KRiIjknpJCLT390ntp5W5bb5FQJCIiuRdrUjCzPmY228zmmNnlGbY3MbPHw+1vmllJnPHkwuiJ05MOQUQkNrElBTNrBNwFHAF0AwaYWbdKzQYBy9x9G+A24Ka44omDVlgTkYZm/RiPvTswx90/AzCz0cCxwKyUNscC14SfxwF/MzPzyi8CrKPFS1fw1AvvrPNxlixbkVY+9qAd1/mYIiL1SZxJoQPwVUp5PrBHdW3cvdTMlgObAd/mMpAfVqxi0msf5vKQALRp1SznxxQRSVKcYwqZpgytfAcQpQ1mNsTMZpjZjCVLlmTYJf+OO3inpEMQEcm5OO8U5gOdUsodgYXVtJlvZusDLYCllQ/k7iOAEQC9evWqdddS61abcFbffWq7W7U2a7kJu3brnLPjiYjUF3EmhelAVzPrAiwA+gOnVGozATgdmAr0BV7M9XgCQMtmG2lQWEQkgtiSQjhGcB4wCWgE3O/uH5rZdcAMd58A3Ac8YmZzCO4Q+scVj4iI1CzOOwXcfSIwsVLdVSmffwL6xRmDiIhEpzeaRUSkgpKCiIhUUFIQEZEKSgoiIlJBSUFERCpYDK8FxMrMlgDz6rh7a3I8hUYB0DkXB51zcViXc97S3WtcFazgksK6MLMZ7t4r6TjySedcHHTOxSEf56zuIxERqaCkICIiFYotKYxIOoAE6JyLg865OMR+zkU1piAiItkV252CiIhkURRJwcz6mNlsM5tjZpcnHU+umFknM5tiZh+Z2YdmdmFY38rMnjezT8Pvm4b1ZmZ3hD+H98xsl2TPoO7MrJGZvW1mz4TlLmb2ZnjOj5vZBmF9k7A8J9xekmTcdWVmLc1snJl9HF7vvRr6dTazi8O/1x+Y2Sgz27ChXWczu9/MFpvZByl1tb6uZnZ62P5TMzt9XWJq8EnBzBoBdwFHAN2AAWbWLdmocqYU+L277wDsCfwuPLfLgcnu3hWYHJYh+Bl0Db+GAHfnP+ScuRD4KKV8E3BbeM7LgEFh/SBgmbtvA9wWtitEw4Bn3X17YEeCc2+w19nMOgAXAL3cvQfB9Pv9aXjX+UGgT6W6Wl1XM2sFXE2w3PHuwNXliaRO3L1BfwF7AZNSylcAVyQdV0znOh44FJgNbBHyYC8fAAAGGElEQVTWbQHMDj8PBwaktK9oV0hfBKv4TQYOAp4hWNb1W2D9ytecYD2PvcLP64ftLOlzqOX5Ngc+rxx3Q77O/LJ+e6vwuj0DHN4QrzNQAnxQ1+sKDACGp9SntavtV4O/U+CXv1zl5od1DUp4u7wz8CawubsvAgi/tw2bNZSfxe3AZUBZWN4M+N7dS8Ny6nlVnHO4fXnYvpBsBSwBHgi7zO41s41pwNfZ3RcANwNfAosIrttMGvZ1Llfb65rT610MScEy1DWoR67MbBPgCeAid/8hW9MMdQX1szCzo4HF7j4ztTpDU4+wrVCsD+wC3O3uOwMr+aVLIZOCP+ew++NYoAvQHtiYoPuksoZ0nWtS3Tnm9NyLISnMBzqllDsCCxOKJefMrDFBQnjU3Z8Mq78xsy3C7VsAi8P6hvCz6A0cY2ZfAKMJupBuB1qaWflKgqnnVXHO4fYWBEu/FpL5wHx3fzMsjyNIEg35Oh8CfO7uS9x9DfAksDcN+zqXq+11zen1LoakMB3oGj61sAHBYNWEhGPKCTMzgnWuP3L3W1M2TQDKn0A4nWCsobz+N+FTDHsCy8tvUwuFu1/h7h3dvYTgWr7o7qcCU4C+YbPK51z+s+gbti+o3yDd/WvgKzPbLqw6GJhFA77OBN1Ge5rZRuHf8/JzbrDXOUVtr+sk4DAz2zS8wzosrKubpAdZ8jSQcyTwCTAXuDLpeHJ4XvsQ3Ca+B7wTfh1J0Jc6Gfg0/N4qbG8ET2LNBd4neLIj8fNYh/M/AHgm/LwVMA2YA4wFmoT1G4blOeH2rZKOu47nuhMwI7zWTwGbNvTrDFwLfAx8ADwCNGlo1xkYRTBmsobgN/5BdbmuwJnhuc8BzliXmPRGs4iIVCiG7iMREYlISUFERCooKYiISAUlBRERqaCkICIiFZQUpCCYmZvZLSnlS8zsmhwd+0Ez61tzy3X+c/qFM5xOqVRfYman1PGYr+cmOpGAkoIUitXACWbWOulAUoWz8EY1CPitux9Yqb4EyJgUUt7ezcjd967Fny9SIyUFKRSlBEsRXlx5Q+Xf9M3sx/D7AWb2spmNMbNPzOxGMzvVzKaZ2ftmtnXKYQ4xs1fDdkeH+zcys7+a2fRw/vqzU447xcweI3iJqHI8A8Ljf2BmN4V1VxG8bPgPM/trpV1uBPY1s3csWENgoJmNNbOngefMbBMzm2xmb4XHPbaac33Jfllz4dHwTWDC854VnsPNtf7JS1HJ+luISD1zF/Cemf2lFvvsCOxAMA/OZ8C97r67BQsSnQ9cFLYrAfYHtgammNk2wG8IphLYzcyaAK+Z2XNh+92BHu7+eeofZmbtCeby35Vgvv/nzOw4d7/OzA4CLnH3GZVivDysL09GAwmmhe7p7kvDu4Xj3f2H8E7pDTOb4FXfPN0Z6E4w781rQG8zmwUcD2zv7m5mLWvxs5MipDsFKRgezAD7MMHiK1FNd/dF7r6aYHqA8v/U3ydIBOXGuHuZu39KkDy2J5hD5jdm9g7BlOSbESxwAjCtckII7Qa85MFEbqXAo8B+tYi33PPuXj6hmwHXm9l7wAsE0yJvnmGfae4+393LCKY8KQF+AH4C7jWzE4D/1iEWKSJKClJobifom984pa6U8O9y2GWyQcq21Smfy1LKZaTfKVf+rbt8SuLz3X2n8KuLu5cnlZXVxJdpGuO6SD3+qUAbYFd33wn4hmCun8pSz3UtwWI0pQR3NU8AxwHP5ig+aaCUFKSghL89j+GXZRgBviDoroFgDv7GdTh0PzNbLxxn2IpgVatJwLnh9OSY2bYWLG6TzZvA/mbWOhyEHgC8XMM+K4BmWba3IFhDYo2ZHQhsGeF8CGPeBGjh7hMJusp2irqvFCeNKUghugU4L6V8DzDezKYRzCpZ3W/x2cwm+M97c+Acd//JzO4l6IJ5K7wDWULw23a13H2RmV1BMMWzARPdfXy2fQhmPi01s3cJ1uxdVmn7o8DTZjaDoFvo41qcVzOCn82GYTxVBupFUmmWVBERqaDuIxERqaCkICIiFZQURESkgpKCiIhUUFIQEZEKSgoiIlJBSUFERCooKYiISIX/D9Gl1rF1Hqk5AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cdf = Cdf(suite)\n", "thinkplot.Cdf(cdf)\n", "thinkplot.Config(xlabel='Number of trains',\n", " ylabel='Cumulative Probability',\n", " legend=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Cdf` also provides `Percentile`" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(69, 869)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cdf.Percentile(5), cdf.Percentile(95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise:** To write a likelihood function for the locomotive problem, we had\n", "to answer this question: \"If the railroad has `N` locomotives, what\n", "is the probability that we see number 60?\"\n", "\n", "The answer depends on what sampling process we use when we observe the\n", "locomotive. In the book, I resolved the ambiguity by specifying\n", "that there is only one train-operating company (or only one that we\n", "care about).\n", "\n", "But suppose instead that there are many companies with different\n", "numbers of trains. And suppose that you are equally likely to see any\n", "train operated by any company.\n", "In that case, the likelihood function is different because you\n", "are more likely to see a train operated by a large company.\n", "\n", "As an exercise, implement the likelihood function for this variation\n", "of the locomotive problem, and compare the results." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solution\n", "\n", "# Suppose Company A has N trains and all other companies have M.\n", "# The chance that we would observe one of Company A's trains is $N/(N+M)$.\n", "# Given that we observe one of Company A's trains, the chance that we\n", "# observe number 60 is $1/N$ for $N \\ge 60$.\n", "\n", "# The product of these probabilities is $1/(N+M)$, which is just the\n", "# probability of observing any given train.\n", "\n", "# If N<>M, this converges to $1/N$, which is what we saw in the previous\n", "# solution.\n", "\n", "# More generally, if M is unknown, we would need a prior distribution for\n", "# M, then we can do a two-dimensional update, and then extract the posterior\n", "# distribution for N.\n", "\n", "# We'll see how to do that soon." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }