{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------------------------------\n", "\n", "\n", "Copyright 2016 Allen Downey\n", "\n", "MIT License: http://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "import thinkbayes2\n", "import thinkplot\n", "\n", "import numpy as np\n", "from scipy import stats\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's an update that takes the prior probability of being sick `prior`=$x$, and the likelihoods of the data, `prob_given_sick`=$P(fever|sick)$ and `prob_given_not` =$P(fever|not sick)$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def update(prior, prob_given_sick, prob_given_not):\n", " suite = thinkbayes2.Suite()\n", " suite['sick'] = prior * prob_given_sick\n", " suite['not sick'] = (1-prior) * prob_given_not\n", " suite.Normalize()\n", " return suite['sick']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we start with $x=0.1$ and update with the assumption that fever is more likely if you're sick, the posterior goes up to $x\\prime = 0.25$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.24999999999999997" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prior = 0.1\n", "prob_given_sick = 0.9\n", "prob_given_not = 0.3\n", "\n", "post = update(prior, prob_given_sick, prob_given_not)\n", "post" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now suppose we don't know $s =$ `prob_given_sick`=$P(fever|sick)$ and $t = $ `prob_given_not` =$P(fever|not sick)$, but we think they are uniformly distributed and independent." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.5, 0.5)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dist_s = thinkbayes2.Beta(1, 1)\n", "dist_t = thinkbayes2.Beta(1, 1)\n", "dist_s.Mean(), dist_t.Mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compute the distribute of $x\\prime$ by drawing samples from the distributions of $s$ and $t$ and computing the posterior for each." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n = 1000\n", "ss = dist_s.Sample(n)\n", "ts = dist_t.Sample(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just checking that the samples have the right distributions:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{'xscale': 'linear', 'yscale': 'linear'}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXCXvYFBCQVcCtWtGqdfnW1qhfFauW1qWK\n1oWKYtXWpVq1rYX61Yo7Wq0VF9xatSr9ubRW65JarVaQnSRAEghkI+tknyQz9/z+mJBklpCBZO5k\nZt7PxyOPR+acM5OP95G8vZx77rnGWouIiCSntHgXICIisaOQFxFJYgp5EZEkppAXEUliCnkRkSSm\nkBcRSWLdhrwx5hljzA5jzNpdjHnUGLPZGLPaGHNE75YoIiJ7Kpoz+aXA6V11GmPOAGZYaw8A5gN/\n7KXaRESkh7oNeWvtp0D1LobMBl5oG/tfYKQxZlzvlCciIj3RG3PyE4HtnV4XtbWJiEic6cKriEgS\n698Ln1EETO70elJbWxhjjDbKERHZA9ZasyfvizbkTdtXJG8B1wKvGmOOAzzW2h1dfZA2RAtYuHAh\nCxcujHcZfYKORQcdiw6pdizuWraOd/6bj7fFF9a3ZvF5e/y53Ya8MebPQAYw2hizDVgADASstXaJ\ntfbvxpjvGmNygQZg7h5XIyKSAr7Mq2LJx/nslT4ArOWDFXldjp3ulLOmBz+r25C31l4UxZjrelCD\niEhKaGz28bMXVwPgqWtkxYbyLseOsE3sY+t57eGr6PfoT/b4Z/bGnLzsgYyMjHiX0GfoWHTQseiQ\nbMciq6iWh97dREV1HVuLq7ocN86p4fGbzmL65DGkpfV8bYxxc47cGGM1Jy8iqaTU4+XXr68HILeg\nDE99U3tfum1hjK0DIA3L72//IZPH7x32GcaYmF94FRGR3VRU1cSCZRsAKNxRHRTw051yDLDfxDE8\ncMu5GLNHGd4thbyISAy89FkBmdmBOfdGbwulFbUADLdexto6JuwzkgdvPZ+BA2Ibwwp5EZFe9sy/\ntvD55koAcvJLqG9qAWCqU0l/HE469iCuu+gkV2pRyIuI9KI7/5rFtspGAJpbWtsDfobTsZLGrYAH\nhbyISI80t/p5e1UJpTVeVhd42tsd67BuczF72UZG2wYgcEfpa4vnu1qfQl5EZA+0+h1+snRlxL6q\nmgbyCysY2SngR40cylN3XuJmiYBCXkRkt9V7fdzw0uqIfWWVtWwrrWasU8twmtvb/7jgYrfKC6Jd\nKEVEdkNFXXOXAT99qJ+W4q1Md8qDAn7p3ZfRr1984lY3Q4mIRGl9YQ2L/7E5rP3+Cw5l3q+fD2tP\nHzyQF+/9cY9/bk9uhlLIi4hEoai6iQVvbAhqa2n1UZ+zusstet945Ope+dm641VEJIayi2p58N1N\nQW1bCisYWb01YsD/6OxjmX3y4e4U1w2FvIhIN0IDvqm5tcuA//P9VzBo4AB3CouCQl5EpAuOY7nq\n2a+C2pqaW2jatDYs4J/73eUMHzrYveKipJAXEYlgS3kDd7+Z3f66samZrPxSZrRtLNZZb829x4JC\nXkQkRG1Ta1DAb9xaSl1DM/s6nrCxfTngQSEvIhIk9Ay+sLSauoZmpjiVDMAJGtvXAx4U8iIiAPgd\nS1F1U1DAbyuppKyqPmhzsb2Gp/P7X11I+pCB8ShztynkRSSl+fwOf/woP2hzMbCs2LANCDzcY6fr\nLzmZ7xx9oMsV9oxCXkRS1qbSOu57Z2NQW0l5DUVlgcCf3uki634TxyRcwINCXkRSUGOzj5+9GL7/\nTENTc3vAT+sU8Il4Br+TtjUQkZTi8ztcHbJF8M4HbO/reEinNajvoVvPZ+qE0W6WGEbbGoiIRKGw\nqpGFy7KC20qraaqrY4qtScjVM91RyItISqiqbwkL+PpGL0PK8xgaMrZfvzReffBK94qLIYW8iCS9\ntds9PPpeblCbz+enNm8jQzq1ffuoA7jh0lPcLS7GFPIiktR+svQrWv2h1wItddmrggIeSLqAB4W8\niCSxa59fGRbwdQ1eWvPXJ9T+Mz2hkBeRpGOt5cpnvgprb/B48G3PDQv4pXdf5k5hcaBnvIpIUrDW\nsrrAw0dZZRED/obT96dle25Y+/P3zGXEsNCJm+ShdfIikvB+//5m1myr6bJ/zvGT+cOSvwa17T0i\nnafuvARj9mj5uat6sk5eZ/IiktA+WL9jlwF/34UzGTvAF9b+9P9dmhAB31OakxeRhLVyazWvfLE9\nrP3rk0aQPiCNmfv058pfPRvW/8SCi90or09QyItIQiqr9fKHD/KC2g6dNIJzZo7mpntfA+C9CO87\nduY0xo4a7kKFfYPm5EUk4Xhb/Fz3wqqgtr3SBzCNaj76b06X7xsyeCAv3HM5aWmJNVPdkzl5hbyI\nJIxI4Q5QW9eIb2tW2NLInSbsM5LbrjqDiWP3im2BMRLzDcqMMbOAxQQu1D5jrb03pH8E8BIwBegH\nPGitfW5PChIRicRxbMSAz84rYVxjUcSAT9YbnHZHt2fyxpg0YBNwClAMLAcutNbmdBpzOzDCWnu7\nMWYMsBEYZ631hXyWzuRFZLeFn8Fb6hubydmyI+jBHjtd+N1vcv7pR7lZYkzF+kz+GGCztbag7Ye9\nAswGOk98WWDnlYzhQGVowIuI7Alrg8/gHcdhZfZ29nMqmEHwSePXD5jAwmvPTomlkdGKJuQnAp3X\nKBUSCP7OHgPeMsYUA8OAC3qnPBFJZZ6GFm5+eW1Q28rs7UxzKkgLCfhbfnwaxx0+3c3yEkJvLaE8\nHVhlrT3ZGDMD+KcxZqa1tj504MKFC9u/z8jIICMjo5dKEJFkYK3l8Q/yQh6sHZBfWM4UpzIo4A+a\nNp6b557KqJGhu8InrszMTDIzM3vls6KZkz8OWGitndX2+jbAdr74aox5B7jHWvtZ2+sPgVuttStC\nPktz8iLSpZc+KyAzuzxiX0urj9KcbIbT3N72wqK5DB0yyK3y4ibWc/LLgf2NMVOBEuBCYE7ImALg\nf4HPjDHjgAOB/D0pSERST3Orn2ufD185A1Bb30Tu1lKm2iqGdzqDv+6ik1Ii4Huq25C31vqNMdcB\n79OxhDLbGDM/0G2XAHcBzxljdk6e/cJaWxWzqkUkaXyZV8WSjyOfE7aWFOCvKGdahL6Tjj0otoUl\nCd0MJSJx8++N5Tz/74Kw9v2H+Vn+eeQze4A/3XcFgwcNiGVpfUrMb4YSEeltpR5vhIC3+PKzWN7Q\nFDZ+8r6juPv62Zqi2U0KeRFxjc/vsDy/mmf+tSWsL32AoWjligjvgovPOpZzTv1GrMtLSgp5EXHF\n0k+28Nmmyi56bcSAP3C/cdxz4w9iW1iSU8iLSEzVe3089O4mtlU2Ru6vb6BlS3bY1gRPLLg4pbYE\njhWFvIjEzK7WvR84fhj/+ee/6IcNC/jXHr4q4bYD7qsU8iLS6/yOZf6z4Q/TBpgyOp07vv81zrvh\nSfpF6H/mrksV8L1IIS8iverL/CqWfBR53fstZx7EAeOGct4NT4b1nXfakcw5M3RbLOkphbyI9JqV\nW6sjBvysmeM575hJfLluK+cv+lNY/6/mf5cjD5niRokpRzdDiUivqGls5ed/XhPW/tBFM/n1w8sY\nPGgAedvD5+d1gbV7evyfiMTV05n5fJEbvJPJlNHpTMHD25lru3gX/N/PZnPIjH1jXV7C0x2vIhIX\n1lqufCbyBdYDBtbz+vtdB7wezecOhbyI7LGuAv7So0fzyNL3wtp/PvdUpk4YnbAP1E5Emq4RkT0y\n7+nwO1S/MXUv9qWOP73z36D2GZP34b6bz3WrtKSj6RoRcdXv398c1vazEyfw28feijheAR8/OpMX\nkahFmoN3rENJdhbpfm/E92juved0Ji8iMbetspE7/5oV1NbobaFs82ZGEh7wRx4yhV/N/65b5UkX\nFPIi0q1I8+/bSqportzBGIL3fh+WPojnfnc5xuzRiaf0MoW8iETk8zvc8vJa6ry+sL4VGwqY7FQx\nHH9Q++N3XMT4MSPcKlGioJAXkTB+x3L10pVh7S2tPtZuKmKGE37nqube+yZdeBWRIF3d4NTc0kr1\nxg0MIvjM/pYfn8Zxh093q7yU1JMLr9rPU0TaWWu55eXwu1R/cuIkGjeuCQv42686QwHfx+lMXkQA\n+PfG8qAHa5eU11BS7mGyv5L+OGHjNT3jHi2hFJEeuevNLLaWdzyer6XVR1GZhwmOJ2LAv/bwVW6W\nJz2gkBdJYV9tqeaJD/PaX/v9DtW1jWwtrmSsU8sQWoPG9+uXxqsPXqnlkQlE0zUiKcRaS1ZRLcvz\nq/l0U0XnHqpqGskvrGBv28AoG/zQ7ftvPpfJ40cxYECkB/ZJrGm6RkS69Zs31lNcHXxnakNTM9n5\npe2vJzrVDA65uPrqg1fSv7/CPVEp5EWSXFV9C794JXzFzJaiCio9DQAMsj4m2mo6nyoeN3Ma1196\nigI+wSnkRZJcaMBba/kqaxuDrI+9aGFv2xC0lvqoQ6byy/lnuFukxIxCXiSJZRXVBr0uq6rFX5TP\njJApGQM8v2guQ4cMcrE6cYNCXiQJ1Xlb+SK3ile/2A6Az+dn9cZCpjvlDAgZe/D08dx9/ffdL1Jc\noZAXSRLeFj93vLGe6oaOZY87p2YAhtgWQpdnPLHgYsaOGu5ileI2hbxIEnjty+28t3ZHUFvngAfY\n19a0f3/czGnccsXprtUn8aOQF0lg1lqWfrKV/2yubG8rr66joLiq/fVw62UfW9d+Fv+ba87i8IMm\nuVypxItCXiRBeVv8XPfCqqC2sqo6tpVUMdmpYmDIXu/P3zOXYem6sJpqFPIiCeKvK4rYUt4QtmJm\np6qaBgqLK5hhK8P6HvnlBQr4FBVVyBtjZgGLCWxN/Iy19t4IYzKAh4EBQLm19qRerFMkpf3xwzxW\nbKmO2Nfc0kp5VT2NFTuYZuvb279+wASuvegk9tl7mPaaSWHdhrwxJg14DDgFKAaWG2PetNbmdBoz\nEngcOM1aW2SMGROrgkVSzY4ab8SAr/TUs6UocNY+zSlnaKe+mQdOYsG1Z7lUofRl0ZzJHwNsttYW\nABhjXgFmAzmdxlwEvGGtLQKw1laEfYqI7Baf34n4CL69+7WwbnUOA/ExHSdsWWS/fmkKeGkXTchP\nBLZ3el1IIPg7OxAYYIz5GBgGPGqtfbF3ShRJPb9/fzNrtnUseWz1+dhWUo3P72dsXWHQWftOe49I\n5/gjpnPFuSe4V6j0eb114bU/cCRwMjAU+NwY87m1NreXPl8kJUQ6e/fUNZK7rZw06zA1wkXV7xx9\nANdcmKFtgCWiaEK+CJjS6fWktrbOCoEKa60X8BpjPgEOB8JCfuHChe3fZ2RkkJGRsXsViySpSBdX\n6xq85G4rZ5pTHvZA5gn7jOT3v57jXoHimszMTDIzM3vls7p9aIgxph+wkcCF1xLgS2COtTa705iD\ngd8Ds4BBwH+BC6y1WSGfpYeGiITwNLRwc4SHZxfuqMaWFTKM5qD2a+dkMGPKWKZOGOVWiRJnMX1o\niLXWb4y5DnifjiWU2caY+YFuu8Ram2OMeQ9YC/iBJaEBLyLBymq9fLShjA82lIX1rdhQwFQn+AHa\nacbw2uL5bpYoSUCP/xNx0ebSOu59Z+Mux+Suz2a0bQhrf+ORq2NVlvRxevyfSAJ49Yvt/HP9jl2M\nsFSsW8logk+EzjvtSOacGbqgTSQ6CnkRF3hb/F0GvDFw5uH7suqLr6gOCfjfXnc2Xz9gohslSpJS\nyIvE2N9Xl7BsRfCCtEED0njo4sMZ1L8fdQ1eLv/lc2Hv+801ZyngpccU8iIx5Dg2LOABHr/sSCDw\nxKZIAb/gmrOYqe2ApRco5EVi5LNNFSz9ZGtY++OXfQOfz8+jf/qYz1aG3y94+rcOVcBLr1HIi8SA\nt8UfMeCfnnc0BcVV3HTvXyK+7y8PXUW/fqG3PYnsOYW8SAyEPswD4LZZ0zj/xiU4jhPhHfDUnZco\n4KXXKeRFepHfsfxnc/gmrEcPb+DW+1+L+J5FN/2A/aeM1Z7vEhO6GUqkh8rrmrn91XWMTB9ATWNr\nUF+rz09d9qqw7YB30g1OEo2e3AylkBfpAWstVz7zVVi733EoLK2msaqS8Tb8cX0P3Xo+UyeMdqNE\nSQK641UkDiI9SBsCwb8qeztjnHrG0hTU94NTjuBH3zvOrRJFwnYvFZEoNLdGDvgbT5tB1fqVzHDK\nGUlT0B/YL644XQEvrtN0jchusNay6O0c8sqCNxDz+x1qN2+A1uaw9/z4nG9x5omHuVWiJCFN14i4\nJNL8e6vPR2326oj/LH75gXkMHKA/M4kf/faJROG5T7by6abwpZE7KmtJK94c8Q/p9cXztSxS4k4h\nL7ILjmNZ9E4O+SHTMw1NzXhyc0inNew9D9/2Q6bsq6c2Sd+gkBfZhaueDZ6esdbyVdY2xjs1DI0Q\n8E8suJixo4a7VZ5ItxTyIl1Y9HZOWNv6DXnMsJ6w9qMOmcov55/hRlkiu0Wra0RCvL+ulL/8tzCo\nzef3U561jsH4gtoP3X8Ct847naFDBrlZoqQYra4R6SV3vL6eEo83qM3b3ErB5nzGhgT8s3ddxsjh\nQ9wsT2S3KeRFgOX5VTz5UX5QW3Orj3WbipjgVIcF/E2Xn6qAl4SgkJeUFrr3TKvPz/bSKjy1jaQ7\nXmbYurD36AxeEolCXlJSYVUjC5dlARa/Y/H5/BSVeaiqaQRgtFPPXiH7zsyYvA93XT9bNzdJQtFv\nq6SUL/OqWPJxYFpma1EFFZ7g9e/ptpnxtjZsa2Dd2CSJSiEvSa+yvplbX1nX/nrnWvdQE5xqhnSa\ne7/q/G/zzcP2Y9TIoa7UKRILCnlJauW1zdz+l0DAO9ZhZdb2sDFTnUr60/FIvktnH89ZJx6mR/FJ\nUlDIS1IKvaDqdxxWZXcEfLptZqytox+Wa+acyCnHfS0eZYrEnEJekk5owDc0NZOdXwrAFKeK/vjb\n59zvvekc9p86Ng5VirhDIS9Jp3PA19Q3sbmgjEG2lYnW0x7uz/3ucoYPHRyfAkVcpJCXpDLv6RXt\n31d66tlSVMkEx8OQTpuJ/em+Kxg8aEA8yhNxnUJeksKX+VUs6XTHamFpNaWVtUxxKhnQ6aLqG49c\nHY/yROJGIS8J78aXVlPn7Vj62OrzU1pZy3inRgEvKU8hLwmroq6Z215dF9TmWIc1GwsZaRsZSgsA\nB08fz10/mx2PEkXiTlsNS8JZu93Do+/lBrU1eVvYkFdCum1hnK0Jet6qzuAl0WmrYUkJb68q5s2v\nittfN7f6qK5pxNvcQoWngbFOLcNpDnrPyw/Mc7tMkT5FIS99mrWWTzdV8Py/C9rbqmsbyNse/FDt\nyU4VA/EHtT1712XaTExSXlR/AcaYWcBiIA14xlp7bxfjvgn8B7jAWrus16qUlOPzO1y9dGVIq2V1\nTiE+f+Bi6lDbzLgIm4ktvfsyRgzTVsAiEEXIG2PSgMeAU4BiYLkx5k1rbU6EcYuA92JRqKSGirpm\n7v/bRirrW9rbrLVs3LqD+sbAVMxI28QYW9/en5aWxqnHf42J4/bizBMPc71mkb4smjP5Y4DN1toC\nAGPMK8BsIPQpxz8FXge+2asVSkoo9Xj59evrw9pzt5XhqQvs6z7YtjKh012rMw+cxG+uOVNbAIvs\nQjQhPxHovHVfIYHgb2eMmQB831p7kjEmqE+kO45jgwLecRxWdtpMLNJ8+2O/nsO++4x0rUaRRNVb\nV6UWA7d2eq1TK4lK6LNV128uwtviw2CZ5FSHhfvP557K/xwxw+0yRRJWNCFfBEzp9HpSW1tnRwOv\nmMC/m8cAZxhjWq21b4V+2MKFC9u/z8jIICMjYzdLlmTgbfFz3Qur2l933ilyhlMeNl7bAUsqyczM\nJDMzs1c+q9uboYwx/YCNBC68lgBfAnOstdldjF8KvB1pdY1uhhJrLU9+lM+KLdXtbZ66RnK3BYJ9\nP6eCfnT8jmhaRiTGN0NZa/3GmOuA9+lYQpltjJkf6LZLQt+yJ4VIcsvbUc+yFUVsLKkDoMrTQFFZ\nNf379aPB28IQ28IEW9M+/u7rv8/B08fHq1yRpKFtDSSm6r0+bnt1Ld7Wjo3CcvJLqG/qWCIZevb+\n0r0/Zsjgga7WKdKXaVsD6XNCn87U1spXWdvY+f/5SNsQaJ8Zkd6lkJeY+EnY3aqwYsM20qzDGNvA\nMLykEVjrfsZ3vs6kcXsxYexe7hcqkuQU8tKrFr2dQ+6O+pBWy4oN29jX8ZDe6QlNF511DOeeeqS7\nBYqkGM3JS6/p/Oi9zorXrQl6/B7APTf+gAP3G+dGWSIJT3PyElclniYWvR26y0XA3i3VVIcE/EO3\n/pCpE0a5UZpIytOZvPTIQ+9uIquoNqit1ednb28F27dsDxv/xIKLGTtquFvliSQFncmLqzwNLXyY\nVca7a0rD+rYUVtC/uoT6kFUzAN/42mQFvIjLFPKyW/6zuYJn/7U1Yt+KDQVhF1d3unT28cw++fAY\nVycioRTyErWiqqaIAV9d20BDQR4zaAnru/z7/8PpJxyiJzSJxIn+8iQqO2q8LFi2Iahte0kVvooS\nhuNlaMhuFuedfhRzvqtHC4jEm0JedslxLNc8txKfszPELQXFVZRX10fcLRLgdzd8n4Omad8Zkb5A\nIS9dCt4OOHBD006RAn765H24/+ZzXapORKKhkJeI1hfWsPgfmwFwrMPKrMByyIHWxwTrCRp77qlH\nMuvbhzJq5FDX6xSRXVPIS5jtlY3tAV9T18TmbWUAjHdqGNrp4uq1czI4+biD41KjiERHIS9hfvvX\nLCCwJHKn6U55+zMdH7/jIsaPGRGHykRkdynkJci8p1fg9zsUl3VMyXSef39h0VyGDhkUj9JEZA8o\n5AUAn9/h6qUr8dQ2kru9I9SntwX8kMEDeXHRXAKP8RWRRKGQT3EtPofnPtnKF3kVlFfVs72049mr\n45waDDBp3N488ssL4lekiOwxhXwKyy6u5cG/byI7r4QGb/DdqhOcaobg4/gjZvDzy/83ThWKSE8p\n5FNQU4ufn7atf29qbg0L+GlOBWlYPYpPJAko5FPMnX/NYltlIxDYEnhDbnF73wTHwxBa+e53vs6P\nz/lWvEoUkV6kkE8RtU2t3PSnNe2vK6rr2Fpc1f56hlPOzAMnseDas+JRnojEiEI+ybX4HK55ruOh\n2g1NzWTnd+wDP9WppD8OAL+55kzX6xOR2FLIJ7HVBR4e+2du++v8wnKqagJTNWOcOkbibe97+YF5\nWh4pkoQU8knGWsvWikae/CiPirrABdUqTwP5RRXtY/ZzKujXtjXwsPRBPH/P3LjUKiKxp5BPEo3N\nPh59P5fcHfVB7Zu27qC2IXDG3nlrAoAXF/2Y9CEDXaxSRNymkE9wjmO56tmvIvZ56hqpbfCSZh32\ns5XtAf+DU47g4rOP1fSMSApQyCcoay33vrMx7Mx9p7oGL7nbypnoVDMIHwY4aNp4fvajk7W5mEgK\nMdba7kf11g8zxrr585KRtZZ3Vpfw5lfFEftbfX4qstcxGF9Yn25uEklMxhistXv0T2+dySeQv68p\nYdnyorD2nWved15QHRzhvQp4kdSkkE8AZbVefvmX9RH7VudsZ3Srhxk0R+w/6pCp3HjZKbEsT0T6\nME3X9HGZ2WW89Nm2oDZrLdZainKyGeZv6vK9Ty78EWP2HhbrEkUkxjRdk6Q+21QREvCWDZtLSGuu\nZ5ytJTS+jzxkCj+58EQ9a1VE2ink+yCf3+GhdzexqXTnyhmLYy0rs7YH3cjU2c9+dDInfvNAdwsV\nkT5PId8HXb00sNdMc0srWXkl+J1AqHd+DF9n9998LtMn7+NafSKSOBTyfcy8p1cAkF9YQVVNAwCD\nbSsTrSdo3H0/P5cZUxTsIrJradEMMsbMMsbkGGM2GWNujdB/kTFmTdvXp8aYw3q/1OT37poSGpqa\nWbGhoD3gJziesIBffPsFCngRiUq3q2uMMWnAJuAUoBhYDlxorc3pNOY4INtaW2OMmQUstNYeF+Gz\ntLqmCy//axP3/r91QW2TnSoG4gc0JSOSymK9uuYYYLO1tqDth70CzAbaQ95a+0Wn8V8AE/ekmFTU\n4nM4/+F/UdDpAR4Ao5x6BuLn/FlHceEZ34xTdSKS6KIJ+YnA9k6vCwkEf1fmAe/2pKhU8WV+Fbe/\n8CWVbVMznR2yz0AeuOUiBg8aEIfKRCRZ9OqFV2PMScBc4ISuxixcuLD9+4yMDDIyMnqzhITgbfFx\n9qIPqW3w0twSvMfMMOsl84EL6d+/X5yqE5F4y8zMJDMzs1c+K5o5+eMIzLHPant9G2CttfeGjJsJ\nvAHMstbmdfFZKT8n3+Rt4fjb347Y9519DY/c8gNtASwiQXoyJx/N6prlwP7GmKnGmIHAhcBbIQVM\nIRDwl3QV8AJVNQ0RA36QbeV7Mwbx6C/OUcCLSK/qdrrGWus3xlwHvE/gfwrPWGuzjTHzA912CXAH\nMAr4gwmkVKu1dlfz9innmt+9zn/KwwP8kTmHMfPgSew9Ij0OVYlIstMGZTHmqW3kgjteYUfayLC+\nZbedxvRxw+NQlYgkkp5M1yjkY2hFdjHzlnwese/tX53O5DHaIVJEuqddKPsYn99h3lPLWb2xMKxv\nyr6jWPbzE+nfL6qbjUVEekQh38v+lV3OH/+5iewtpWF955+wP788Z6YuroqIaxTyvcTn8/O3VcUs\n+WATRWWesP7P7zmbIYMHxqEyEUllCvkechyH829cgpf+FKXtHdY/2ali6cI5CngRiQuFfA9Yaznn\nxqfYZkbjmOA59n7WYT9byeuL52t6RkTiRqtrdlNFXTMLlm1g5JAB/P3zjUDkAL/wkKHcduUsd4sT\nkaSkJZQumvf0ChqbmsnKD7+wevC0cRx/0DjmnzSdwQO194yI9A4toXRBflk9v122gey8Elp8/rD+\nIw6axO2zD+GgfXVzk4j0HQr5bvgdy6V/+JwNeSVhfZOcKgzw5B0XMGHMCPeLExHphkJ+F8o9DZz6\n239E7BtiWxg/YhBP3XmJLqyKSJ+lkA/R3Orn89xKnvoojw25xWH9A6yPqWl1/O76szlwv3FxqFBE\nJHopH/IJ4hqEAAAH0UlEQVR+x7K6wEO918eLnxUE2vxOxICf5lTw2oPz9EAPEUkYKR3yG0vquP9v\nG9tfW2spLvNQUlEbNvbfd5/F8PRBbpYnItJjKRvy855egd9x8NQ24m1ppbK6IWzVzBDbwhhbz9uP\nXBmnKkVEeialQt5ay4PvbmJtQTVrNxWF9Q+3XlrpxzhbS38cAF5fPN/tMkVEek3KhHy918cNL62m\npLwmbAOxNOswzVYGtV11/rc57VuHaOWMiCS0pA95ay0fZZXx8ufbydlSSn1jc1D/ZKeKgfj51pH7\nM3HsXsw64VBGDh8Sp2pFRHpXUof8vKdXtH+/YkNBUN8+Th0j8PLCorkMHaILqiKSnJIy5BuafVz/\n4mogsBxyVc72oP7pTjnfOHgyd/zkzHiUJyLimqQL+Y+yyvjzf7YBgamazgG/l21ktG3g5rmncfwR\n0+NVooiIa5Im5K21vPBpAf/eWAFAQXEl5dX1AEx1KttXyyy89mwOO3Bi3OoUEXFTUoR8qcfLr19f\nj7WWdZuLaGntWO++n1NBPwLbG7/28FWkpekB2iKSOhI+5B/4+0ZyiutwHIeV2cFz7zOc8vbvn7rz\nEgW8iKSchA352qZWbvrTGiD84uq+jod0WgFNz4hIaku4kP/NG+sprva2v96QW0xTcyDQDZbpTmBO\n/s6ffo9D958QlxpFRPqKhAn5LeUN3P1mNs2tPkrLa9ovqu7U+exda99FRAISIuR33tS0paiCSk9D\nWP80p5yds+2vPHAlAwZoK2AREUiAkH/0vc3UN3rJ2bIjqH2cU8MQWttXztx70znsP3VsPEoUEemz\n+nTIV9W38ObneVTWdJy9D7Nexto6DHDe6Udx5ne+zohh2mtGRCSSPhvyz3+4kcf/kR20x/vO/Wbm\n//A7nPatQ+JYnYhIYuhzIV9U1cQVT3xKaYSnM40Z4OPP98/X9r8iIlEy1lr3fpgxdlc/b9Ff1/HK\nJ5si9l1+1Ghu+FFGjCoTEem7jDFYa/fo7LZPnMkXVTXx85e+Cru4CjDStPDkdRkcPH3fOFQmIpLY\n4h7yuaV1XPjAB/j8TljfDSdN4vLvHRuHqkREkkNUIW+MmQUsBtKAZ6y190YY8yhwBtAAXG6tXd3V\n55XVeln0dg6e2kZWbioJ6x/t1PP2ojmkDxkY5X+GiIhE0u2OXcaYNOAx4HTgUGCOMebgkDFnADOs\ntQcA84E/dvV5WYUezr3vAz5akRcx4Kc75Xz4yGVJH/CZmZnxLqHP0LHooGPRQceid0SzLeMxwGZr\nbYG1thV4BZgdMmY28AKAtfa/wEhjzLhIH3bRgx9S19AcqYt3bjuZZY9cHW3tCU2/wB10LDroWHTQ\nsegd0YT8RKDzHr6FbW27GlMUYUxEk5wqzj04nVUPncOkcXtH8xYREYlSXC+8PnXFURx96FStexcR\niZFu18kbY44DFlprZ7W9vg2wnS++GmP+CHxsrX217XUOcKK1dkfIZ7m3KF9EJInEcp38cmB/Y8xU\noAS4EJgTMuYt4Frg1bb/KXhCA74nRYqIyJ7pNuSttX5jzHXA+3Qsocw2xswPdNsl1tq/G2O+a4zJ\nJbCEcm5syxYRkWi4uq2BiIi4KyZPtjbGzDLG5BhjNhljbu1izKPGmM3GmNXGmCNiUUdf0N2xMMZc\nZIxZ0/b1qTHmsHjU6YZofi/axn3TGNNqjDnHzfrcFOXfSIYxZpUxZr0x5mO3a3RLFH8jI4wxb7Vl\nxTpjzOVxKDPmjDHPGGN2GGPW7mLM7uemtbZXvwj8jyMXmAoMAFYDB4eMOQP4W9v3xwJf9HYdfeEr\nymNxHDCy7ftZqXwsOo37EHgHOCfedcfx92IksAGY2PZ6TLzrjuOxuB24Z+dxACqB/vGuPQbH4gTg\nCGBtF/17lJuxOJPv1ZunEly3x8Ja+4W1tqbt5RdEeX9BAorm9wLgp8DrQJmbxbksmmNxEfCGtbYI\nwFpb4XKNbonmWFhgeNv3w4FKa63PxRpdYa39FKjexZA9ys1YhHxMb55KMNEci87mAe/GtKL46fZY\nGGMmAN+31j4BJPNKrGh+Lw4ERhljPjbGLDfGXOJade6K5lg8BhxijCkG1gDXu1RbX7NHuRn3XSgl\nwBhzEoFVSSfEu5Y4Wgx0npNN5qDvTn/gSOBkYCjwuTHmc2ttbnzLiovTgVXW2pONMTOAfxpjZlpr\n6+NdWCKIRcgXAVM6vZ7U1hY6ZnI3Y5JBNMcCY8xMYAkwy1q7q3+uJbJojsXRwCsmcAv0GOAMY0yr\ntfYtl2p0SzTHohCosNZ6Aa8x5hPgcALz18kkmmMxF7gHwFqbZ4zZAhwMrHClwr5jj3IzFtM17TdP\nGWMGErh5KvSP9C3gUmi/ozbizVNJoNtjYYyZArwBXGKtzYtDjW7p9lhYa6e3fU0jMC9/TRIGPET3\nN/ImcIIxpp8xJp3AhbZsl+t0QzTHogD4X4C2OegDgXxXq3SPoet/we5Rbvb6mbzVzVPtojkWwB3A\nKOAPbWewrdbaY+JXdWxEeSyC3uJ6kS6J8m8kxxjzHrAW8ANLrLVZcSw7JqL8vbgLeK7T0sJfWGur\n4lRyzBhj/gxkAKONMduABcBAepibuhlKRCSJxeRmKBER6RsU8iIiSUwhLyKSxBTyIiJJTCEvIpLE\nFPIiIklMIS8iksQU8iIiSez/A8cnjD4vy1iUAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(thinkbayes2.Cdf(ss))\n", "thinkplot.Cdf(thinkbayes2.Cdf(ts))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now computing the posteriors:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "posts = [update(prior, s, t) for s, t in zip(ss, ts)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the distribution of values for $x\\prime$ looks like:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{'xscale': 'linear', 'yscale': 'linear'}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHJlJREFUeJzt3Xl0VOXdB/DvLxsQQkJYAkIgYAgJi4Js8gracYOAtW5Y\nxbdatSp91fe0tm9dWn0bezxV39O6tGoVDopoEa2oIFrFhVHZUYGw7wQISNiFQMIk83v/yBjmTpYZ\nJjP3mbnz/ZzDOfPceXLn6z3hy/WZO3dEVUFERM6UZDoAERFFD0ueiMjBWPJERA7GkicicjCWPBGR\ng7HkiYgcLGjJi8hUEdknIqXNzPmbiGwWkZUiMjiyEYmIKFyhnMm/AmBsU0+KyDgA+apaAGASgBcj\nlI2IiFooaMmr6gIAh5uZchWA6b65SwFkiUiXyMQjIqKWiMSafHcAu/zG5b5tRERkGN94JSJysJQI\n7KMcQA+/ca5vWwMiwhvlEBGFQVUlnJ8LteTF96cxcwDcA+BNERkJ4Iiq7mtqR7whWp2SkhKUlJSY\njhETeCxOc+Kx8Hq9WLF+F7bs3I9FK7Zi977m3uI7bd2Sueg/8sdRTmePwt5dsafiCIp6d8VZnbOC\nzvd6FV71on9+N2RmtMbAgvBXwIOWvIjMAOAC0FFEdgL4I4A0AKqqk1X1QxEZLyJbAFQCuC3sNEQU\ns2prvZi3aB1qaryW7VWnPFi+egdyOmY2+JljlSexZvMeuyI2KTOjDVJTGl+dVgUOHa3Eef16oHWr\ntEbnVJ/yQCC4bsx5EAn9hDo7Mx2dO7QLK3OkBC15Vb0phDn3RiYOEdlNVbF3/1Hs2HMQ4vsf9qPH\nTmLrrv3ISG+Fg0crsfDbLUH3s3XX/mhHRddOmUgSQVJSEgYX9Wh0Tq3Xi5yO7XDB4HwkJQmyM9PP\nqJidJhJr8hQGl8tlOkLM4LE4LZrHYvWmcuwoP4iFK7agc4d2EBHU1nqxZNW2qL1moMLeXdEpOwOj\nzsvHoMJctG6V2uRct7uIvxsRIHaukYuIck2e6MypKjaXVeDw9yfqt9XUelG6cTcy0lsBAE5UncJX\n32xBr24drT8LxYZt30U0z7gLByIp6fTZsaqi8uQpDOnXs9F37/K6dUSPrtkRzZBIRCTsN15Z8kQx\n4mTVKSwt3Y7Kk6eweOU2dGjfFkki+OqbzbZlyM5MR99edZ9lrKr2QBUYVJQLr1fRpnUqRg/pg3Zt\nW9uWh+q0pOS5XENko+MnquH11r1xuWL9Liwr3Y5NZRU4dLTStgzDB/ZCamoyhvTridSUZABA5w4Z\n6NurS0KvXTsVS54ogo6fqMb+Q8cAAPsPH8eu7w5hxtxlyOnQDhW+7S01qDAXaamn/+ru3X8UFw0v\nqH/TFAAKeze8s0iHrLYhXb5HzsLlGqJmHDxyHEtLt2Pd1u+wYdtedO/S3vL8D5cHpqYkw1NTG7HX\nHXlubwDA4H49kJaagiQRDB2Qh/Q2jV/iR87GNXmiFqqq9mDfwe/rxx5PLR546p2ovd4Pb5YeP1GN\ni88vRJ8eOTi7RycU5OVwyYQa4Jo8UYhOVp2C1+9E45SnFnc8Mj3ir3NW5yykpaZgT8UR9OmZg47Z\nbfHT4mHIbNuab1ySrVjy5CherxflFUehqvh4wVocPFKJ5Wt2oFvnLOzZfzTs/SYlJeG8oh7ol98V\n+T06W55LS01B79yO9Wfg/uvlRKZxuYbi1payCjz6wlwU5OVg1cbdEdlnZkYbZLVrAwDweGrw3YHv\n8doTt3MtnIzimjwljFOeGny0YC1efW9xi/aT3vp0aZ+oOgUAmPKnm9Ehq22L9ksUDVyTJ8cr23MQ\nv3nyX2f0M7ldsrF732GMv2ggCnt3RV63jsjKaI3MjDZRSkkUe1jyFFOOHDuBJSu349t1O7FjzwF4\narz4/vjJoD/3h0njkZqSjN65ndCmVSqSk/l9OEQAS54MqKmpxaayChw/UY1N27/DstU7UF5x5Iz3\nU5CXg8fvu4aXHBI1gyVPUXX02Els3lmBbbv24+s1ZS2+HW166zRMf+I2FjtRiFjyFBFbyiqwbPUO\nzPrkWwiAjLatcayyKuz9JScnod/ZXTFsQC/k9+yMtm3S0D2nPVJ891ohotDw6hpqkW/WluHPk/8d\n9s8X5OXg+IlqDC7qgeHn9MLZuZ2Qkd6KZ+pEfnh1DdnK6/XiowVr8fGCdSF/XycADCzohvaZ6Rh1\nXh8MH5jHIieyAUueQvbGh8uxrHQ7du491OScnmd1wDl9u2PchQOR3joNKSlJaNumlY0picgfS56a\ntXbLHkx+66ugZ+zJyUl466m7bEpFRKFiyVMDXq8XM+Yuw7ufrQw6d/TQPrhoaAGGDsizIRkRnSmW\nPNU7VlmFW38/Lei8yy/oh0tHFqEgr+EXUxBRbGHJJ7h5C9fhpbe+DDovLTUFf7z7xyjsza+II4on\nLPkEs3f/USxeuQ3/nLs0pPlX/Ogc3H7tqCinIqJoYckniMqT1bj7TzNw/ER1SPN5R0YiZ2DJO9y6\nrXvxyN9mB503fGAv/ObWy/iFF0QOw7/RDlVx6Bj+69F/Njtn1JA+uO7yIeh5VjbX2YkciiXvILW1\nXkyfvQRzvyhtdt6zv78BuV2ybUpFRCax5B3i98+8h43bv2vy+e457fHgncXoltPexlREZBpLPs4t\nX7MDT0z5qNk5/3r6LiQl8Us0iBIRSz5Ord+6Fw8384bqT4uH4afFQ7nWTpTgWPJxRlUx7d3FTa67\nP37fNejbi59EJaI6LPk48p/3T0VVtafR5y45vwg/v/o/kJHOOz4S0Wks+TjwxofL8fbH3zT5/Kxn\nf2ljGiKKJyz5GPfo83NRuml3o89NvGIEJowZYnMiIoonLPkYVbpxNx59YW6jz5XccyXO6dvd5kRE\nFI9CKnkRKQbwDIAkAFNV9cmA5zMBvA6gJ4BkAH9V1WmRjep8qzeV4+///BwHj1Q2+vx/DM7H/9x2\nuc2piCieBf0ibxFJArAJwKUA9gBYDuBGVd3gN+chAJmq+pCIdAKwEUAXVa0J2Be/yLsJR46dwC8e\nnt7k8/dMdOGSkUU2JiKiWBHtL/IeAWCzqpb5XmwmgKsAbPCbowDa+R63A3AwsOCpae9+ugKvv9/4\nrX8L8nLw+H3X8Hp3IgpLKCXfHcAuv/Fu1BW/v+cAzBGRPQAyANwQmXjOpqqYPnsJ5sxf1eC5v94/\nAb26dzKQioicJFJvvI4FsEJVLxGRfACfiMi5qno8cGJJSUn9Y5fLBZfLFaEI8WfCr19qdPtbT92F\n5GTehoAoUbndbrjd7ojsK5Q1+ZEASlS12Dd+EID6v/kqInMBPK6qC33jzwA8oKpfB+yLa/I+j734\nAVas39VgO695J6JA0V6TXw6gj4jkAdgL4EYAEwPmlAG4DMBCEekCoC+AbeEESgRvz/u2QcHfM9GF\ni88vNJSIiJwqaMmraq2I3AtgHk5fQrleRCbVPa2TATwGYJqI/HBDlftV9VDUUsexo8dO4o0Pllm2\n3XfLZRg9tI+hRETkZEGXayL6Ygm+XHP02Enc/vCrlm0dstpiyp9uNpSIiOJBS5ZrWPI2OVl1Cj97\n4OUG27kGT0TBtKTkeQmHDVS10YJ/8693GkhDRImE966JMlXFs6993mA7v62JiOzAko8SVcWMucvw\nzqcrGjz3xl/uYMETkS1Y8lFQVe3Bf94/tdHnbhw/HGmpPOxEZA+eTkZBUwXvGlGI68cOtTkNESUy\nnlJG2PMz3A223X7tKIy/aCBvMkZEtmPJR9DxE9X4fOkGy7aXH/s5stq1MZSIiBIdl2siZO2WPfj5\nQ69Ytt169QUseCIyimfyETDt3UV4313aYPuVF59rIA0R0Wk8k4+AxgqeH3QioljAM/kWmvKvrxps\n460KiChW8Ey+BWpqavHRgrWWbSx4IoolLPkWuOG3UyzjOydcaCgJEVHjWPJh+u3/vd1gW/GFAwwk\nISJqGks+DF99vRk7yg9Ytk1+9GeG0hARNY0lf4ZqamrxzGufWbb9790/Rsf2GYYSERE1jSV/hu59\nbKZl3L5dOgYV5hpKQ0TUPJb8Gdp/+JhlPPWxWwwlISIKjiV/Bt75xHpv+Cd/c62hJEREoWHJn4F/\nzl1qGffJyzGUhIgoNCz5EG3asc8y/snFgwwlISIKHUs+BKqKh55+17LtlqtGGkpDRBQ6lnwIFq7Y\nahlnZ6bzC0CIKC6w5EPw9KufWsZT/nSzoSRERGeGJR/EnPmrLOOh/fN4Fk9EcYMlH8Sr7y22jO//\nxRhDSYiIzhxLvhmvzVliGU+8YgRSUpINpSEiOnMs+SaoKt77bKVl24QxQwylISIKD0u+CS+/s9Ay\n/tmV5xtKQkQUPpZ8Ez78co1lfM1l5xlKQkQUPpZ8Iw4drbSMf3vb5YaSEBG1DEu+EV99s8UyvmBw\nvqEkREQtw5JvxPTZi4NPIiKKAyz5AKpqGV84tMBQEiKilgup5EWkWEQ2iMgmEXmgiTkuEVkhImtE\nZH5kY9pn0cptlvHEK4YbSkJE1HIpwSaISBKA5wBcCmAPgOUiMltVN/jNyQLwPIAxqlouIp2iFTja\nnpr2iWXcpWOmoSRERC0Xypn8CACbVbVMVT0AZgK4KmDOTQBmqWo5AKjqgcjGtEfg3SZ7dY/bf6uI\niACEVvLdAezyG+/2bfPXF0AHEZkvIstFJC5v0xh4Fv/wL8cbSkJEFBlBl2vOYD9DAFwCoC2AxSKy\nWFW3NP9jsePIsROW8cCCbsjOTDeUhogoMkIp+XIAPf3Gub5t/nYDOKCqVQCqRORLAIMANCj5kpKS\n+sculwsul+vMEkfJvY/NtIx/d/tYQ0mIKNG53W643e6I7EsCLxlsMEEkGcBG1L3xuhfAMgATVXW9\n35wiAH8HUAygFYClAG5Q1XUB+9Jgr2fKdb960TKe9ewvDSUhIrISEahqWF9kEfRMXlVrReReAPNQ\nt4Y/VVXXi8ikuqd1sqpuEJGPAZQCqAUwObDgY1ngPzx/mMS1eCJyhqBn8hF9sRg9k5/rLsUr7y6q\nH/MsnohiSUvO5BP+E6+VJ6stBU9E5CQJX/IPPzvbMr792lGGkhARRV5Cl7yqYufeQ5Zt4y8aaCgN\nEVHkJXTJPzfDbRlPf+I2iIS17EVEFJMSuuTdyzZaxm3btDKUhIgoOhK25KuqPZbxHRNGG0pCRBQ9\nCVvyT0371DIuHj3AUBIiouhJ2JLfvLPCMuZaPBE5UcKW/PfHT9Y/vnPChQaTEBFFT8KWvL8BBd1M\nRyAiioqELPk1m6030eyek2UoCRFRdCVkyS9aYf0e16SkhDwMRJQAErLdPll0+gaZ+T06G0xCRBRd\nCVnyXr87YZ4/qLfBJERE0ZVwJX/KU2MZjxx0tqEkRETRl3AlP3+p9VYG3TrzTVcicq6EK/mv15ZZ\nxvwQFBE5WcKV/LfrdtY/HtCH18cTkbMlVMkHfvXgmFH9DSUhIrJHQpX8/sPHLeMR5/QyE4SIyCYJ\nVfLPz5hvGaelphhKQkRkj4Qq+TWb99Q/Tk1JNpiEiMgeCVXy/m69+gLTEYiIoi5hSt7r9VrGl44s\nMpSEiMg+CVPy/ks1AJCayuUaInK+hCn5Fet3mY5ARGS7hCn5OfNX1T/O7ZJtMAkRkX0SouQPHrFe\nH3/Fj84xlISIyF4JUfJ/e/1zy3g4PwRFRAkiIUp+/yHrmXx2ZrqhJERE9kqIkt938Pv6x3ddf6HB\nJERE9nJ8yS9audUyHjYwz1ASIiL7Ob7kX5610DLu2D7DUBIiIvs5vuT97y7ct1cXc0GIiAxwfMkf\nOXai/vFPLh5kMAkRkf0cXfI1NbWWcX7PzoaSEBGZEVLJi0ixiGwQkU0i8kAz84aLiEdEro1cxPCt\n3brXMs7p0M5QEiIiM4KWvIgkAXgOwFgAAwBMFJEGt3D0zXsCwMeRDhmul9780nQEIiKjQjmTHwFg\ns6qWqaoHwEwAVzUy778BvA2gIoL5WsT/+vhO2byqhogSTygl3x2A/y0cd/u21RORbgCuVtV/AJDI\nxQtfba31/vH33XKZoSREROZE6o3XZwD4r9UbL/pJJa9bxoW9efkkESWeUL7JuhxAT79xrm+bv2EA\nZoqIAOgEYJyIeFR1TuDOSkpK6h+7XC64XK4zjByaw9+fsIzrohERxT632w232x2RfYn6f1qosQki\nyQA2ArgUwF4AywBMVNX1Tcx/BcD7qvpOI89psNeLhOpTHtz0u6n14z//+moU9u4a9dclIooGEYGq\nhnWmGvRMXlVrReReAPNQt7wzVVXXi8ikuqd1cuCPhBMkkj78co1lzIInokQVynINVPUjAIUB215q\nYu7tEcjVIq+/v9R0BCKimOC4T7wGLgf1OKuDoSREROY5ruSnvbvYMn7ivmsMJSEiMs9xJT/3i1LL\nuHWrVENJiIjMc1TJBy7V3DBumKEkRESxwVElv333Acv4+rFDDSUhIooNjir53/1llmXMD0ARUaJz\nVMn7y0hvZToCEZFxjin5xSu3WcZ//8ONhpIQEcUOx5T8X16ZZxlnZrQxlISIKHY4puT9uUYUBp9E\nRJQAHFHypRt3W8Z33/gjQ0mIiGKLI0r+0RfmWsbJyY74zyIiarG4b8MtZdZvGxzaP89QEiKi2BP3\nJT/l7QWW8f2/GGMoCRFR7In7ki/bc7D+cft26UhJSTaYhogotsR9yXtqausf/2LCKINJiIhiT1yX\nfOXJast4QH43Q0mIiGJTXJf8ivW7LOOsdvwAFBGRv7gu+adf/dR0BCKimBa3JX/KU2MZ988/y1AS\nIqLYFbclv6P8oGX8yH9dYSgJEVHsituSf+eTFZZxWmqKoSRERLErbkt++Zod9Y9zu2SbC0JEFMPi\nsuQ/+GK1Zcyv+SMialxclvzL7yy0jEcP7WMoCRFRbIu7kq+q9ljG5/XrYSgJEVHsi7uSf3X2Ysv4\nD5PGG0pCRBT74q7k5y1cZxmLiKEkRESxL+5K3t81lw42HYGIKKbFVcmrqmV8ycgiQ0mIiOJDXJX8\nvoPHLOOzOmcZSkJEFB/iquT//eUay5jr8UREzYurkp/7RWn94zat0wwmISKKD3FT8qs3lVvGYy7o\nZygJEVH8iJuSL3n+fcv45p+MNJSEiCh+xE3J+zu3by7X44mIQhBSyYtIsYhsEJFNIvJAI8/fJCKr\nfH8WiMg5kQwZeCuDe25yRXL3RESOFbTkRSQJwHMAxgIYAGCiiAReoL4NwEWqOgjAYwCmRDLkvEXW\nT7l2ys6I5O6JiBwrlDP5EQA2q2qZqnoAzARwlf8EVV2iqkd9wyUAukcy5PylGyO5OyKihBFKyXcH\nsMtvvBvNl/gdAP7dklCBdu49VP/4omEFkdw1EZGjRfQ780TkYgC3ARjd1JySkpL6xy6XCy6Xq9l9\nllccsYyvuCiiy/1ERDHH7XbD7XZHZF8SeD+YBhNERgIoUdVi3/hBAKqqTwbMOxfALADFqrq1iX1p\nsNcLdN2vXrSM335mEq+sIaKEIiJQ1bCKL5TlmuUA+ohInoikAbgRwJyAAD1RV/A3N1Xw4TjlqWmw\njQVPRBS6oMs1qlorIvcCmIe6fxSmqup6EZlU97ROBvAIgA4AXpC6Fvao6oiWhlvwzRbL+I2/3NHS\nXRIRJZSQ1uRV9SMAhQHbXvJ7fCeAOyMbDZj9+SrLOC01om8hEBE5Xsx+4rWq2oPd+w7Xjwt7dzWY\nhogoPsVsyf95svUqzGsvP89QEiKi+BWzJb9h+3eW8bABeYaSEBHFr5gt+dpab/3j39x6ucEkRETx\nKyZLvuKQ9Wv+BhflGkpCRBTfYrLkA+9V07ZNK0NJiIjiW0yW/FsffV3/ODOjjcEkRETxLeZK/rU5\nSyzjH7t4rxoionDFXMm/99lKy/i6y4cYSkJEFP9iquQD7zg5YexQQ0mIiJwhpkp+8ltfWsYTxw83\nlISIyBliquTXbN5T/zi9dZrBJEREzhAzJV9TU2sZ/8/tYwwlISJyjpgp+a/XllnGA/LPMpSEiMg5\nYqbkn5n+mWWckpJsKAkRkXPERMlXnqyGx2+5pj/P4omIIiImSv6WB1+xjO+56WJDSYiInMV4yVdV\nexps69op00ASIiLnMV7ygW+4TvvzrWaCEBE5kPGSf/rVTy3jdm1bG0pCROQ8Rkv+wOHjlnFBXo6h\nJEREzmS05P1vKQwAv79rnKEkRETOZLTkP1uywTLmveOJiCLL+Jr8D650nWs6AhGR4xgr+Q++WG0Z\n/+zK8w0lISJyLmMl//I7Cy1j3saAiCjyjJS8qlrGY0cNMBGDiMjxjJT8xwvWWcZ3Xj/aRAwiIscz\nUvJT3v7KMhYREzGIiBzP9pL3er2W8aUji+yOQESUMGwv+V3fWb+s+67rL7Q7AhFRwrC95H/75FuW\nMa+qISKKHttL3v+6ms7Z7ex+eSKihGL0E6+/vuVSky9PROR4xkp+aP88FJ3d1dTLExElhJBKXkSK\nRWSDiGwSkQeamPM3EdksIitFZHCwffLaeCKi6Ata8iKSBOA5AGMBDAAwUUSKAuaMA5CvqgUAJgF4\nMdh+O2VnhBXYKdxut+kIMYPH4jQei9N4LCIjlDP5EQA2q2qZqnoAzARwVcCcqwBMBwBVXQogS0S6\nNLfTRP8AFH+BT+OxOI3H4jQei8gIpeS7A9jlN97t29bcnPJG5hARkc2MvPH6zEM3mHhZIqKEI4F3\nhGwwQWQkgBJVLfaNHwSgqvqk35wXAcxX1Td94w0AfqSq+wL21fyLERFRo1Q1rDXulBDmLAfQR0Ty\nAOwFcCOAiQFz5gC4B8Cbvn8UjgQWfEtCEhFReIKWvKrWisi9AOahbnlnqqquF5FJdU/rZFX9UETG\ni8gWAJUAbotubCIiCkXQ5RoiIopfUXnjNRofnopXwY6FiNwkIqt8fxaIyDkmctohlN8L37zhIuIR\nkWvtzGenEP+OuERkhYisEZH5dme0Swh/RzJFZI6vK1aLyK0GYkadiEwVkX0iUtrMnDPvTVWN6B/U\n/cOxBUAegFQAKwEUBcwZB+AD3+PzASyJdI5Y+BPisRgJIMv3uDiRj4XfvM8AzAVwrencBn8vsgCs\nBdDdN+5kOrfBY/EQgMd/OA4ADgJIMZ09CsdiNIDBAEqbeD6s3ozGmXxUPjwVp4IeC1VdoqpHfcMl\ncO7nC0L5vQCA/wbwNoAKO8PZLJRjcROAWapaDgCqesDmjHYJ5VgogB9uWdsOwEFVrbExoy1UdQGA\nw81MCas3o1Hy/PDUaaEcC393APh3VBOZE/RYiEg3AFer6j8AOPlKrFB+L/oC6CAi80VkuYjcbFs6\ne4VyLJ4D0F9E9gBYBeBXNmWLNWH1ZiiXUJINRORi1F2VlMh3bnsGgP+arJOLPpgUAEMAXAKgLYDF\nIrJYVbeYjWXEWAArVPUSEckH8ImInKuqx00HiwfRKPlyAD39xrm+bYFzegSZ4wShHAuIyLkAJgMo\nVtXm/nctnoVyLIYBmCl1NzbqBGCciHhUdY5NGe0SyrHYDeCAqlYBqBKRLwEMQt36tZOEcixuA/A4\nAKjqVhHZDqAIwNe2JIwdYfVmNJZr6j88JSJpqPvwVOBf0jkAbgHqP1Hb6IenHCDosRCRngBmAbhZ\nVbcayGiXoMdCVc/2/emNunX5ux1Y8EBof0dmAxgtIskiko66N9rW25zTDqEcizIAlwGAbw26L4Bt\ntqa0j6Dp/4MNqzcjfiav/PBUvVCOBYBHAHQA8ILvDNajqiPMpY6OEI+F5UdsD2mTEP+ObBCRjwGU\nAqgFMFlV1xmMHRUh/l48BmCa36WF96vqIUORo0ZEZgBwAegoIjsB/BFAGlrYm/wwFBGRgxn9jlci\nIoouljwRkYOx5ImIHIwlT0TkYCx5IiIHY8kTETkYS56IyMFY8kREDvb/f0JjbjLfbgoAAAAASUVO\nRK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cdf = thinkbayes2.Cdf(posts)\n", "thinkplot.Cdf(cdf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's the mean:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.15586200108737797" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cdf.Mean()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "This result implies that if our prior probability for $x$ is 0.1, and then we learn that the patient has a fever, we should be uncertain about $x\\prime$, and this distribution describes that uncertainty. It says that the fever probably has little predictive power, but might have quite a lot.\n", "\n", "The mean of this distribution is a little higher than the prior, which suggests that our priors for $s$ and $t$ are not neutral with respect to updating $x$. It's surprising that the effect is not symmetric, because our beliefs about $s$ and $t$ are symmetric. But then again, we just computed an arithmetic mean on a set of probabilities, which is a bogus kind of thing to do. So maybe we deserve what we got." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just for fun, what would we have to believe about $s$ and $t$ to make them neutral with respect to the posterior mean of $x$?" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{'xscale': 'linear', 'yscale': 'linear'}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdcVGe+x/HPQ+9VRQUUFbFr7MaKvcT0aDRt0930srtJ\n7t27d912k2w2m6bZdE0zarpGo2LB3kvsBcQCKIr0zgzP/QMycAAFcZiBmd/79fL1cn7nMPPLCXw5\nPuec51Faa4QQQjgmF3s3IIQQovFIyAshhAOTkBdCCAcmIS+EEA5MQl4IIRyYhLwQQjiwOkNeKfWx\nUipNKbX/Cvu8rZQ6oZTap5S6zrotCiGEaKj6nMnPAyZebqNSajLQSWvdGZgFvGel3oQQQlyjOkNe\na70JyLzCLjcDn1Xsux0IVEqFWac9IYQQ18IaY/LhwNkqr1MqakIIIexMLrwKIYQDc7PCe6QAkVVe\nR1TUalBKyUQ5QgjRAFpr1ZCvq2/Iq4o/tVkCPAEsUkoNAbK01mmXeyOZEK3c7NmzmT17tr3baBLk\nWFSSY1HJ2Y7F0ZPn+WzJNo4lnSdd+ZKtfADw8fJg6ys3Nfh96wx5pdQCIBYIVUqdAf4MeABaa/2B\n1nq5UmqKUioByAceaHA3QgjhhPYeOcvf31sGQBbeZCtvAMJbBXHHsI5sfaXh711nyGut76rHPk82\nvAUhhHBOZWVlfPLdFn7eeBCATOVDhvLF19uDjhEt6NchlN8Mj+LBa/gMa4zJiwaIjY21dwtNhhyL\nSnIsKjn6sVix8RA/rNnHxcxcALLxJkP5Et4qkDYtgwjwduOh2A64uDRoKN5C2XKMXCmlZUxeCOHM\nSkpNvPzBCvYfT7bUinDjnGswMVFh+Pl4Eezrzp9v7YGfV/l5uFKqwRdeJeSFEMJGTpxO4/V5qy1n\n7wDFuEJYFGEtg1BKEeTjzt/v6ImXh6tln2sJeRmuEUIIG1iz7QjvLdpIWVmZpebm7Y1vZDQenp7l\nr10Uj4/rZAj4ayUhL4QQjUhrzeIVu1m8YpehPnZ4L44X+5BXZAbAzVXxxLhoOrbys+rnS8gLIUQj\nKSk18fq8OHYdOm2ptQz254l7x/LZjjTyikqB8jP430/pQnSYdQMeJOSFEKJRXMjI5V+frCLx7EVL\nrV2bEP7w8GTeiDtJZn6ppX7/yKhGCXiQkBdCCKvKLyzmi6XbWbX5sKE+qFcUz943lg/XnzYE/F1D\n2zEkOrTR+pGQF0IIK9lx4BSvfrSiRn3GlIFMm9ifRdvOsu90lqU+uU9rxnRv1ag9ScgLIcQ1Kisr\n4/3FG1m99YihHhzgw6w7RzKwZxQHzmYTd7ByWq+hnUO5fWBEo/cmIS+EENcgIzuff89fzZGT5wz1\nicN68Mi04SilOJKaw9y4BMu2jq18uXd4e5v0JyEvhBANoLXm27i9fL1yNyaT2VKPbB3MHx6aSHir\nIAAS0vJ4a8UJTGXlD4J6urvw0KgOuLvaZjkPCXkhhLhKyWmZPPN/i2rUbxl7HXdPHYSLS3mAn88q\n4p8/HaUi3/F0d+HFqV0JC/SyWa8S8kIIUU+lpWbe+Gw1O/YnGequri689PAk+nVvZ6nlFJby1soT\nloD3cnfhhaldaRfqY8uWJeSFEKI+Tqde4rVPVnHuYrah3rdbJC8+NAl398qpCIpKzLz601Eu5hZb\nak9N6GzzgAcJeSGEuCKtNd+v3seCn7ZTfXrFl5+7lZioMEMtr8jEX74/VONhpy5t/G3QbU0S8kII\ncRkJpy8w96t4zpzLsNTc3Fy5fXxfbhvXFzc340RixaVmXlt2zBDwtw0MZ3hMC5v1XJ2EvBBCVJOb\nX8Qn321mw64ThnpEWDDP3jeWDhE1Q1trzdzViaRkFlpqw2JCmdKnTaP3eyUS8kIIUcWabUf48OtN\nlFa5LRJg1MAYHrtzlGHs/VcmcxkfxidxOCXHUrttQDhTrrNvwIOEvBBCAOUPNX3641Y27U4w1Lt2\nbM1Dtw2jY2TLWr9Oa81H8UnsTsq01IZEhzC5T+tG7be+JOSFEE7vl2PJvD4vjvzCYkP97qmDuWVs\nH8t977VZvD2ZXVUCvn+HYO4fEYVS17Y2q7VIyAshnNavC3p8s3I3ZVWWJu3bLZLHZowiNOjK0/+u\n3H/eMB9N3/ZB/HZMxyYT8CAhL4RwUikXspg9ZykZ2fmWmrubK7+9cySjBsbUGdRbTqTz9Y7Kxbjb\nhfrwyOimFfAgIS+EcDJaa75ZtYeFy3ca6t5eHvzzd7fRtmLOmSvZeTKDT9afsrwOC/Tk91Ni8HCz\nzXw0V0NCXgjhNC5m5PLKRys5lZJuqA/qFcXv7h9f47732uxOyuT9tSctrwO83fj9lC74eDbNOG2a\nXQkhhBVprVmx6RCf/biNklKTpR4a5Mvt4/sxcXiPer1PQloe/1mTaHnt6+nKCzd0JdjXw+o9W4uE\nvBDCoSUlp/P+4g2cOH3BUO8R3ZYXH56Ir7dnvd7nfFaRYU54Hw9X/nhzN1oF2G5GyYaQkBdCOKxv\nVu3hq2U7DLWqqzXVV0ZeCS8vPUJ+cfkDUm4uimcmdm7yAQ8S8kIIB6S1Zs6CeOJ3HDPURw/uwgO3\nDq332TvAxdxiXl9+zBLwrhUB3ynsyrdXNhUS8kIIh5KVW8Bf5v5kmFQssnUwt47ry6iBMVf1XqmZ\nhby27Bi5RZXj+I+M7kC38ACr9dvYJOSFEA5j96HTvPbJKsO8M907teFPj92Ah/vVxV1hiZlXfjpK\nQXHle/1mRHsGdAixWr+2ICEvhGj2SkvNvPn5Grb9ctJQH9K7A8/fPx7Xq1xP1WQu47VlxywB76LK\n54Qf2tl+UwY3lIS8EKJZO5WSzlufrzUMz3i4u/HEzFiG94++6vczmcv4YN1JzlwqsNQm92nTLAMe\nJOSFEM1UXkEx7y/ewNa9iYYVm9q1CeFPj91ASKDvVb9niamMN1cc5/j5PEttUMcQbunf1god24eE\nvBCiWdFa8/PGg3z87eYa224dex133zi4QfPHFJWY+d1Xv1BcWmapjezagnuHtW9y89FcDQl5IUSz\nsWVfIguX7STlQpah7uPlwXO/GUe/7u0a9L75xSZeX37cEPC9IwObfcBDPUNeKTUJeBNwAT7WWr9a\nbXsA8AXQDnAFXtdaz7duq0IIZ5Wemcfr8+M4fiqtxrY7JvRjxpSBDQ7jrPwS/vrDYXIKK2+T7BEe\nwJPjo5t9wAMorauvP15tB6VcgOPAWCAV2AnM0FofrbLPfwEBWuv/Ukq1AI4BYVprU7X30nV9nhBC\nVLVu+zHmLFhXoz5xWA9m3jAQf9+GP3VqLtP87YfDJGdUrss6pU9rbhsY0eD3bAxKKbTWDfqNU58z\n+UHACa316YoPWwjcDBytso8G/Cv+7g9cqh7wQghxNbTWvPHZGjbvMS7HFxEWzEuPTKJNy8Br/oyv\ntp4xBPztA8OZbOeFt62tPiEfDpyt8jqZ8uCvag6wRCmVCvgBd1qnPSGEMzqSeI7/eftHQ83Px5P7\nbh7C2CHdrPIZP+xKIf7IRcvrcT1aOVzAg/UuvE4E9mqtxyilOgFxSqneWuu86jvOnj3b8vfY2Fhi\nY2Ot1IIQorkrKTXx7lfr2bj7hKEe3iqIV56/DR/va5/SV2vNom1nWX2oclbKPu0CmT448prf21ri\n4+OJj4+3ynvVZ0x+CDBbaz2p4vVLgK568VUp9RPwstZ6c8XrNcCLWutd1d5LxuSFELU6lJDK21+s\nJT3TeG7o7+vFR3+9t14LetTHR/En2ZZQ+eBURIg3L03tipeHdd6/MTT2mPxOIFop1R44B8wAZlbb\n5zQwDtislAoDYoCTCCFEHUpLzcz7fgsrNx8y1KPCW3D/LdfTKybcKp+jtWbx9mRDwIcFePL85Jgm\nHfDXqs6Q11qblVJPAquovIXyiFJqVvlm/QHwd2C+Ump/xZe9oLXOuMxbCiEEUH72/tbna7iUlW+o\n3zl5ANMnDbDqZy3de464g5W3YLYL9eGlG7s2yXVZranO4RqrfpgM1wghgKLiUj79cSurNh821EOD\nfPn9AxOIiQqz6uet3H+er3ckW15Hh/nx/OSmufB2bRp7uEYIIazmUEIqr8+PIzu30FCfNqk/d4zv\nZ7Wx91/tSsowBHxUS59mFfDXSkJeCGET2bmFfPjNJrbuSzTUI8KCeWzGKLp2bG31z0xIy+OjdUmW\n12GBnjw70XkCHiTkhRCNzGwuY+XmQzUmFFPAPTcN4eYxfRpl+oCCYhMfrjuJqax8iNjHw5XfTe6C\nn5dzxZ5z/dcKIWwq7VIO/56/moQzFwz1TpEteXzmKKLCG2eO9tyiUv72wxEy8kqA8oW3X7qxKyF+\n136ffXMjIS+EaBTrth/jvcUbMFVZis/D3Y3pk/pz67i+jfa5BcUmXll61BLwAHcPa0fbYO9G+8ym\nTEJeCGFVOXmFfPB1zbH3kQM68+i0EXh7Ne7Z9PyNp0jLLra8vndYe0Z0admon9mUScgLIaxCa82y\n9Qf4avlOiopLLXU/H0+euGs0g3pFNXoPcQfT2HOqcq75W/q3ZVQ35w14kJAXQlhBfmExr3y4gsOJ\n5wz1gT2jeOKu2GuaDri+Nh9PZ9G2yrkUe0cGMrVv8122z1ok5IUQ12TnwVPM+XIdeQWVQyTubq7M\nmj6S0YO72KSHzcfTmbfhlOV1u1AfHh3d0Saf3dRJyAshGiQ7t5B/zVtV4+x9WL9onpg5Ck8Pd5v0\nsSMxwxDwYYGePDups0PPR3M1JOSFEFctOS2Tv/1nWY0ZI++aOojbx/ezWR97T2XywbrKuRADvN14\nfnIMAd62+QXTHEjICyHqTWvNohW7+C5uL2Zz5aLXUeEteHTacLp0sP5Tq5frY/kv5/l+V4ql1tLf\nk99NiSHUz9MmPTQXEvJCiHpJz8zjrc/XGIZnFPDo9JFMGNbdZn2YzGXM33jKMGWwt4crT0+MpoW/\nBHx1EvJCiCsqKi7ly5+2s2bbMYpLKm+N9HB3478fnWy1+d7ro8RUxpsrjnP8fOUwUbCvO89OjKFN\nkHM+7FQXCXkhxGVt3pvIv+fH1agP6NGep+8dg6+37c6cawv4XpGBPDa2k1NNOHa1JOSFEDUkJafz\n9hdrOXPOuPaPAp65dywjBnS2aT9aa+bEJRgCfnzPMKYPjmiUyc0ciYS8EMIiN7+Ilz9cwbGk84a6\nh7sbI/pH8+i0EVaf770uWmvmrk7kcEqOpTa+Zxh3Dmk6C283ZRLyQggANuw6zlufr61R79qxNX94\ncAJB/j4276nUXMbcuAQOJlcG/KCOIUwfHGHzXporCXkhnFxGdj5vf7GWA8dTamz7y5M30rOz7S6s\nVlViKuOfy45y6mKBpdYjIoAHR0XJEM1VkJAXwkmZTGa+W72XH9b8YrhrJsjfh6mxvZg8oidenvZ5\nqEhrzfwNpwwBP7RzKPcNb4+bq1xkvRoS8kI4Ga0136zaw5K1v1BQVGLY1r97e566Z7RNJhS7nBJT\nGfM2JLHzZKalNrJrC+4d1l7O4BtAQl4IJ7Ln8BnmLFhXYxHtkEBfZt05kgE92tups3KmijH4QynG\nMXgJ+IaTkBfCCeQXFjPv+y2s236sxrY7JvbnptG9bXrPe23KyjTvrk40BPyQ6BDuHyFj8NdCQl4I\nB6a1Zu32o3z24zbDVMBgu5Wa6kNrzfvrTrL/bLalNjymBb8ZIWfw10pCXggHtW77Mb6N28O5i9mG\neuf2rXjqnjGEtwqyU2dGJaYy3l51gqOpuZba0M6hEvBWIiEvhIPRWvPvT1ezZa9xjVV3N1cemTac\nMYO7NpnwNJfpGrdJ9osK4oGRMkRjLRLyQjiQTbsT+Pi7zeTkGS+sjh3Sld/ccr3dx92r0lozNy7B\nEPBDokN4aFQHCXgrkpAXwgEUFZfyyXebWbPtqKHu6eHOK8/fSrs2IXbq7PI+33zaMAbfKzKQh2Nl\nyT5rk5AXohnTWhO/4zhzFqyrse3mMX2458bBuLg0rYeHtNYs3p7MhqPpltqQ6BAeHNnBjl05Lgl5\nIZqpjOx8/u+Dn0lKTjfUo8Jb8D+/nUJwgO3nmqlLWZnmvbWJ7DmVZanFtPbj/hFRuLjIEE1jkJAX\nopkxmcwsWbefxSt2UWoyG7ZNHNaDB28bavOZIutDa81H8UmGgO/YypfnJsfIVAWNSEJeiGbkzLkM\n/j0/jrPnMw31pjAdQV0Wb09mx8nK+el7RQbyxLhOEvCNTEJeiGbg14ea3l+80bCAto+XB0/ePZrB\nvZvueLbWmh92pxJ3MM1S6xcVxGNjO8ldNDYgIS9EE1dcUsrcr9azeU+CoX7HxP5Mm9CvSQ7N/Epr\nzXtrT7I7qfJfHhEh3jw6uqMEvI1IyAvRhO04cIpXP1phqPn5ePLiw5Po3qmNnbqqH601C7aeMQR8\nZIg3z0zsLEM0NlSvkFdKTQLeBFyAj7XWr9ayTyzwBuAOXNRaj7Zin0I4lfzCYr5Yup24zYcN9U6R\nLfnzE1Ob1ENNtSksMfPWyhMkpFWuydozIoDHx0XLots2prTWV95BKRfgODAWSAV2AjO01ker7BMI\nbAEmaK1TlFIttNbptbyXruvzhHB263ce5+0vjMvw+Xh5MH5otyZ533t1F3KK+M/qRM5mVD51Gx3m\nx3OTOuPp3nSHlpoypRRa6waNb9XnTH4QcEJrfbriwxYCNwNVH627C/hWa50CUFvACyGuLCk5nfcX\nb+DE6QuGurubK2+8NJ0WwX526qz+MvJK+J+vD1JW5VxuUMcQHhwVJUM0dlKfkA8HzlZ5nUx58FcV\nA7grpdYBfsDbWuvPrdOiEI7NbC7jw282ErfliKHu6urCmMFdePj24U364uqvMvJKePWno4aAH9m1\nBfcNj7JbT8J6F17dgH7AGMAX2KqU2qq1Trjylwnh3BJOX+Ddhes5nXrJUB89uAu/ufn6Jn3fe1W/\nnMninVXGH/fpgyOY0Ku1nToSv6pPyKcA7aq8jqioVZUMpGuti4AipdQGoA9QI+Rnz55t+XtsbCyx\nsbFX17EQDmLT7gTe+Gy1oRbdrhW/vXMkHSJa2Kmrq3fmUgH/WWOc1njm9ZGM7RFmp46av/j4eOLj\n463yXvW58OoKHKP8wus5YAcwU2t9pMo+XYF3gEmAJ7AduFNrfbjae8mFV+H0Lmbk8vnS7TXue584\nrAePTBverO4fz8ov4R9LjpCZX2qp3Te8PSO7trRjV46nUS+8aq3NSqkngVVU3kJ5RCk1q3yz/kBr\nfVQptRLYD5iBD6oHvBDOzmQyszR+P4tX7Kak1GSp+3h58Mx9Y+2+iPbVOnUxn7dXnSCnsPy/xdVF\n8czEznQPD7BzZ6KqOs/krfphciYvnNSew2eY991mUqstxRcV3oL/fnQSoUFN/86ZqnYlZfDhuiTM\nVa6yPhLbgcHRoXbsynE19i2UQogGSrmQxf++vYSs3AJDvUWwH49MG9Hszt611qw6kMbXO5ItNTcX\nxf0joyTgmygJeSEaQVlZGR99s5mVmw/V2DbzhkHcNLo3Hu7N68evrEzz6aZTbD5eeSdQkI87z0+O\noW2wtx07E1fSvL7LhGgGDhxPYfbcpTXq7m6u/PnxqXRr4nPO1Kag2MTbqxIM0xS0CvDkhRu6EOTr\nYcfORF0k5IWwkjPnMnjlwxWkXcqpse2pu0cTO6iLHbq6dimZhfz75+NkF1TeQRMd5sdzkzvj2Qwe\n0nJ2EvJCWMGqzYf56NtNhrneAfp2i+QPD07A08PdTp1dm9TMQv617Bi5RZV3Aw2PacE9w9rJNAXN\nhIS8ENfgQkYucxes4+CJVEM90N+bfzxzC21aBtqps2t3JDWHOXEJFJeW/+Jyc1FMHxLJmO6t7NyZ\nuBoS8kI00Lrtx3h/8QbDOqsuSvHbGSMZO6SbHTu7duuPXOSrrWcwVdwi6eaieFrugW+WJOSFuEpZ\nuQXM+34Lm3Ybn1gd2DOKWXeOJDjAx06dXTutNQu3nWXNocqZMD3dXXhmYmdiWvvbsTPRUBLyQlyF\nzXsTmbsgnuKSyouQ/r5ePDZjVJNeZ7W+ft5/3hDwvp6u/OGGLkSENN9fXM5OQl6IejCby/hi6XaW\nrPvFUB/YM4on7optNrNFXsm6wxf4bmfl3IMxrf14fFw0fl4SE82Z/N8Tog7JaZm888U6Es5UOcP1\n9uTR6SMY3i/ajp1Zz86TGSzYesbyuqW/J89MlJWcHIGEvBBXsHbbUT74eqPh4mqbloH89ambCAn0\ntWNn1nMkNYeP45P4dVqpAG83XpzaRQLeQUjIC1GL9Mw8/vnxShLPXjTUh/WL5vEZo/DybJ73vVd3\n7Fwub608YbmLJtjXnZdu7CpPsToQCXkhqjmSeI7Z7/6EqcrZe4CfNy8+NJGuHR1npaPEtDzeXZ2A\nyVwe8L6erjw3KYZQP087dyasSUJeiAplZWV8vXIPi1fsMtQH9GjPb2eMata3RlZ3PquIt1aeoKCk\n/BeZt4crL07tKhONOSAJeSGofXjG08Odp+4ezfXXdbRjZ9Z3IaeI15YfswQ8wGNjO0nAOygJeeH0\n1u88zrsL1xuGZ8JCA3jx4Ym0b+tYc6Sn5xbzz5+OWSYbc1HIak4OTkJeOK2i4lL+s2h9jSdXB/WK\n4vnfjMfdwe4uScko5O1VJ8iqMpvkw7Ed6RHRfOfXEXWTkBdO6cy5DP7x/nLSMyvnR2+ua63Wx8kL\neby27Bil5srl+p4Y14m+UcF27ErYgoS8cDrrth/jvcUbDMMzXTq05g8PTnCoi6u/2n82i3dXJ1ru\nonFzUfxmRJQEvJOQkBdOQ2vNu1+tZ+32o4b6tEn9uXPSAJRq0DrJTdrm4+nM23DK8lopeGFqFzq2\nal4Lh4uGk5AXDk9rzZptR/lxzT5SL2Zb6mGhATx19+hmuRxfXYpKzHyyIYk9p7IsNS93F56fHCMB\n72Qk5IVDS0pO5/evfVOj3iO6LS88NBE/H8d78Cc1s5A5cQlcyCm21IJ83HluUgzhIXKbpLORkBcO\nqaTUxD/eX15jxSaAsUO68tiMUQ45PLPzZAafrE8yXGAd1DGEmUMj8fdyjKkYxNWRkBcO51jSef74\n5g/oWrb9+8VpDnfvO5QPSf2wO5Vl+85ZakrBrQPCmdLH8YajRP1JyAuHUVZWxnuLNrBm29Ea2158\neBKDekXZvikb0Fozb8Mptpy4ZKmF+Hnw0KgOdGkjqzk5Owl54RDOns/knx+tMFxYBZg4rAcP3jYU\nNzfHerDpV3lFJv61/BjJGYWWWnSYH4+P60SAtwzPCAl54QB+3niQT77dTJmuHKAJDfLl9w9MICYq\nzI6dNa6TF/L4eH0SadmVF1iDfd353ZQY3F1d7NiZaEok5EWzdSkrj/8sXM/eI2ctNQXcPPY67p46\nCBcXxww6rTUr9p/n2ypL9QH0CA/gqQnRuEnAiyok5EWzU1pqZtWWw3y9cje5+UWWuqeHO3+cNZke\n0W3t2F3jOpdVyKcbT5OQlmeo3z4wnMlygVXUQkJeNCuZOQW88elqDiUYb40c1i+ae24cTKsQx7zQ\nqLVm9aELfLcz2XB7ZIifB3cPbUefdkF27E40ZRLyotnYvj+J9xZtICev8iJjgJ83v71zJIN7d7Bj\nZ43LZC5j/sZTbEvIMNRHdmnBzKHtZPxdXJGEvGjy0i7l8L/vLDHMGAkwckBnZk0f6TDrrdamuNTM\nGytOGIZngn3deXR0Rzq3dsx/tQjrkpAXTVrclsO8t2iDoebj5cEj04YzckCMnbqyDZO5jHfiEgwB\n3z08gCfHR+PhJmfvon4k5EWTVFZWxj8/XsXOg6cMdX9fL954abpDTglcVV6RiT99c5DcIpOlNq5H\nK+4cEumQ0zGIxiMhL5oUrTXLNxxk+YYDnE/PsdQ9Pdx58eGJ9OkSYcfubCOnsJT/W3LEEPBje7Ri\nxvXt7NiVaK7qFfJKqUnAm4AL8LHW+tXL7DcQ2ALcqbX+zmpdCqdQUmrij2/9yMkqi2kDRIW34I+z\nJhMS6Gunzmwnp7CUV5YeJT23xFLrGRHAjCGRduxKNGd1hrxSygWYA4wFUoGdSqkftdZHa9nvFWBl\nYzQqHNuFjFz+8d5yktMyDfXxQ7vxyB0jcHWCO0gu5hTzxorjhimCZ14fydgejvvUrmh89TmTHwSc\n0FqfBlBKLQRuBqrPAvUU8A0w0KodCoemtWbJuv0sWLbDsByfv68X//fsLbRt5Rz3fx9MzubNFScM\ntTuHSMCLa1efkA8HzlZ5nUx58FsopdoCt2itRyulDNuEuJxzF7N57pXFlFYJd4A7Jw9g+qQBdurK\n9o6m5vDOqgRD7a6h7RjTvZWdOhKOxFoXXt8EXqzyWi7/iytavfUIH32zyRDwoUG+zJo+kv492tux\nM9vRWvPz/vN8V20OmjsGRUjAC6upT8inAFUv60dU1KoaACxU5fd2tQAmK6VKtdZLqr/Z7NmzLX+P\njY0lNjb2KlsWzVlBYQlvf7G2xq2RbVoG8u8Xp+Hh7jw3fH219SxrD1+wvHZ3VbwwtSsdWjr+BWZx\nZfHx8cTHx1vlvZTWta2fU2UHpVyBY5RfeD0H7ABmaq2PXGb/ecDS2u6uUUrpuj5POK49h8/wzpfr\nDNMSADwxM5YxQ7raqSv7+HZnMj//ct7yOtTPg+cmxdA6yMuOXYmmSimF1rpBIyR1njZprc1KqSeB\nVVTeQnlEKTWrfLP+oPqXNKQR4bguZuSyYNkONuwyXljs2rE1T941mjYtA+3UmX0s2ZNqCPge4QE8\nPq4Tnu6OubCJsK86z+St+mFyJu9UtNb8FH+AL3/abhh7d3V14YmZsYwa6NjTElSntWbhtrOsOVQ5\nRBMd5sfzk2NkmgJxRY16Ji9EQ5jNZcyeu5TDiecM9d4xEfx2xkjCQgPs1Jl9aK15d00ie09lWWot\n/D14blJnCXjRqCTkhdWlXcph7oJ4Q8CHBvny2ztH0a+78z2aX1BsYk5cAsfPV0401tLfkxendpEh\nGtHoJOSF1ZSVlbE0/gCf/bjVUG/fNpSXn7sFTw/HnRL4cs5cKuD9NYmkVXmKtWdEAI+Pk5kkhW1I\nyAuryMxOHsaTAAAXIElEQVQp4M3PVnPwhHHFplEDY3jyrliHXW/1crILSll54DxxB9OoehlqZNcW\n3D20Pa4u8iiJsA0JeXFNCotK+Gr5TlZsOoTZXGapBwf48MBtwxjWt5Mdu7OPIyk5vP7z8Rp1WYdV\n2IOEvGiwzJwCnv7HQgqKSgz1CcO688CtQ53qwSYov7j6w+5Ulu0zXmwOC/Dk/pFRspKTsAvn+ikU\nVrPr0Gne+HQ1RcWlhvofHpzAkD4d7dSV/WTll/DumkROXsi31DzdXRjVtSW3DgiXdViF3UjIi6uS\nnVvI/B+21HiwaVi/aJ6+ezRubs53t8i+01nMiTNOMObp7sKfb+1OqwB5glXYl4S8qLcDx1N4fX4c\nuflFhvpDtw9jyshedurKfi7mFPPNzmR2JxnnwB/ZtQV3DIzAx1N+vIT9yXehqFNhUQlzv1rPtn2J\nhjkrotu14tn7xjrdtARaa7acuMS8DadqbLtzSCTje8oc8KLpkJAXl1VUXMqabUeZ//0WyqrcB+jp\n4c79t1zPhGHd7didfaRkFjJvQxKnLhYY6hEh3swa05E2Qd526kyI2knIixq01sRtOcL7izfUuv1f\nf7jdaVZsqmpXUgbzNpyiuLTyVlFfT1eGxbRg2qAIymfaFqJpkZAXBiaTmXcWrGPTbuOFRF9vT6ZN\n7M/U2F5OF2bmMs3i7caJxQCGRIcwY0g7/Lzkx0g0XfLdKSxOpaTz7lfrSTx70VLz9/Wib7dIHrp9\nOH4+nnbszj72nc5iwdYzZORVPgvg7eHKk+Oj6dJG7nsXTZ+EvMBkMvPJd1tYufmQoT5yQGcemzHK\n6R5q+tWSPaks2WOcpqGlvye/vyGGUD/n+4Unmifn/OkVFvuOnuVv/1lWo37T6D7cd/MQpxuagfIH\nm15eepRLecYneUd3b8m0QZEysZhoViTknVR+YTHfx+3l+zX7DHUfLw9mTR/J8P7RdurMvg6czWb+\nxlNkF1Q+yRvi58HTE6KJCPGxY2dCNIyEvJPRWvP1yt18vWK34bZIKF/Q478fnYy7k85x/vMv5/h2\np3GN+nahPvzhhi54ezjnMRHNn4S8E0lKTuftL9Zy5lyGod62ZSBP3zuGzu2d8yGenMJS3luTaFjU\nw81Fcd+I9gzt3MKOnQlx7STkncSKjYf48JuNhpqLUkyfPIBbx17nlHPOAKw9fIEFW84Yai39PXlu\ncmeZd0Y4BAl5B3cqJZ0FP+1k9+HThvrw/tE8fPtw/H2dM8gKik18tyuF+CMXDfW+UUE8OLKDDM8I\nhyEh78DWbDvCu1+tN9QiwoJ5fOYounRobaeu7G/nyQw+33SaghKzoX7HoAgm9Xbe4yIck4S8AzKZ\nzHy+ZDs/rd9vqA/qFcXT94zB28vDTp3Z16W8Yr7ensyuarNGRoR48+jojrQNlnlnhOORkHcwJpOZ\nv7z7E4cTK1cn8vPx5Pn7x9OnS4QdO7MfrTUr9p9nyZ5USs2VdxR5e7hy24BwYru1dMrnAYRzkJB3\nIMeSzjPny3WkXsy21Nq2DOTvz9xCoL9znqXmFZl4Z9UJEqus2AQwuFMIM66PxN/L3U6dCWEbEvLN\nnNaaFZsO8dE3m2psG9KnI8/dN9Zp75zZeuISH69PqlF/bGwn+ncItkNHQtiehHwzdvTkeeYuMJ65\n/+r28f2YecNApxyGKCg28eWWM2xPND4P0KWNP89M7CzTEginIiHfDOXmFzHv+y2s33m8xjY3N1f+\n+9HJTjv+fig5mw/WnSS/uPLOGTcXxaQ+rbmlf7gdOxPCPiTkm5ll6w+wcPlOCoqMk2eNH9qN6ZMG\nEBLoa6fO7EtrzfJfzvPD7hSqztbQOzKQ+4a3J8jXOe8oEkJCvpkoKTXxzpfr2LI30VDvENGCZ+4d\nS2Rr5x1j/uV0Fou2n+VCTrGlphSM6xHG9MGyYpNwbhLyzUBWbgGz5/7E2Wpzztw9dTC3je9rp67s\nLz23mI/ik0hIyzPUg33deWpCZ9qFyqyRQkjIN3G7Dp3m9XlxlJSaLDUvT3f+8sSNRLdvZcfO7Cc9\nt5gfdqewLSGjxrb2LXx44YYueDrpTJpCVCch30SVlJr441s/cvKscW6VqaN6c/+t1zvlEESJqYxv\ndyaz/shFTGXGaZKDfd25fWAEQ6JD7dSdEE2ThHwTlJSczhufriblQpah/uRdoxk9uIudurIfrTX7\nz2azYMuZGqs1tW/hw5Q+begXFeSUv/iEqIuEfBOSkZ3P96v3snzDwRrb/vrUTfSIbmuHruxrW8Il\nlu5JJa3KRdVfPT0xmt6RQXboSojmQ0K+ibjcWqszpgzkjgn9nO4stcRUxqLtZ1lfbSpgpeCmfm25\noU8bXFyc65gI0RD1Cnml1CTgTcAF+Fhr/Wq17XcBL1a8zAUe01ofsGajjqqsrIzvVu9j4bIdhnpY\naABP3zOGrh2da+pbrTXrj17ki81namwLD/bm2UmdCZZ73oWoN6WrrfNZYwelXIDjwFggFdgJzNBa\nH62yzxDgiNY6u+IXwmyt9ZBa3kvX9XnOZM/hM8z/fkuNsfeH7xjO+Ou7OdWcM1prdiZl8sHakzW2\ntQrw5P6RUcS09rdDZ0LYn1IKrXWD/ulanzP5QcAJrfXpig9bCNwMWEJea72tyv7bAHl+/AoKCkv4\n5PvNrNt+zFAPDfLlhQcnOt2tkQeTs3lzxYkadTcXxahuLZkxJNLphquEsJb6hHw4cLbK62TKg/9y\nHgZ+vpamHNnabUeZ+1W8oaaAG0b15u4bB+Hh7jyXSc5eKmBOXEKNO2YAOrf247GxnQjwlqmAhbgW\nVk0UpdRo4AFg+OX2mT17tuXvsbGxxMbGWrOFJuvMuQze+nwtp1LSDfXWLQJ46ZHJTjUtwZlLBfz8\nyzl2nsyssc3NRfGnW7oTHuKc898LARAfH098fLxV3qs+Y/JDKB9jn1Tx+iVA13LxtTfwLTBJa51Y\n852cc0w+r6CYxSt28fOGg5RV+W93dXVh3JBuPHzHMFxcnGPq2+JSM0v3nmPF/vM1tgX7uvPsxBgJ\ndyFq0dhj8juBaKVUe+AcMAOYWa2BdpQH/L2XC3hno7Vm0YpdfBe3F7O5zFJXwNB+0Tx02zCnWq1p\n47GLLNmTSmZ+qaHeJsiL6YMj6RkRIOPuQjSCOkNea21WSj0JrKLyFsojSqlZ5Zv1B8CfgBDgXVX+\nk1qqtb7SuL1Dy8kr5M9zlnKm2oRika2DefKu0U5zYbXYZGbRtrPsP5NNVoEx3MMCPLm+cyg3XNdG\nwl2IRlTncI1VP8wJhmu27jvJu1/F15jv/eYxfbj3piFOEWgFxSY+23SaA8nZFJeW1dh+64BwJvdu\nLQ8zCVFPjT1cI+qhtNTM7HeXcvSkcbw5dlAXHr59GN5ejv8Aj8lcxmebTrM14RK1/S5XCv52e09a\nB3nZvjkhnJSEvBXsOnSa9xdtICM731JTwDP3jmXEgM72a8yGTl3M58N1J2udY2Zcj1ZM7tOGQB+5\nHVIIW5OQvwZFxaXMWRDP1n01rzXP/d+7CAsNsH1TNpaSUcirPx2loMRcY9utA8IZ16OVzO0uhB1J\nyDfQ7kOnmbMgnpy8QkP9rqmDuHXsdQ5/W+SZSwW8vyax1jP32G4tmTYoQsJdiCZAQv4qpWfm8f7i\nDew5bJxAq0NECx65YzhdOjj2hGKpmYV8syOZ/Weza93+Xzd2pVOYn427EkJcjoT8VVi8Yhffr95n\nWIrPRSlm3jCIW8dd59B3zqRlF/HtzmT2nMqqdftT46Pp017mdheiqZGQrwetNf/7zhIOJ54z1Ht2\nbsvjM2Mdeuy9sMTM3NUJHE3NrbHNz8uNm/u1JbZbS4f+BSdEcyYhX4ek5HRe+WgF6Zl5hvrz949n\nWN9Oduqq8eVX3Ou+O6nm/DLRYX6M6tqS6zvLeqpCNHUS8lfwy7FkXv1oJcUllU9rBvh58z+zptCp\nXUs7dtZ4svJLWLjtLLtqCXeAWWM6MrBjiI27EkI0lIR8LbTWfLNqDwuX7zTUb4ztzX03D3HIO2eK\nS82sOpjGj7tTa93eMyKAZyZ2lmEZIZoZCflq0i7l8OZnazh+Ks1S83B34+l7xnD9dR3t2FnjKDGV\nsXRvKmsOXaDEVHMKgo6tfHl6Qmf8vORbRYjmSH5yK2itidtyhPcXbzDUI8KCefHhibRt5Vh3jvwa\n7j//UnPaX4A7BkUwpnsrPNwc718tQjgTCXkgK7eA9xdtYMeBU4Z6326RvPDQRIdarelSXjGbj19i\n8/H0Gisyebm7MLVvW4bHtJAzdyEchNP/JG/cdYL3v95IYZVZI0MCfbnnxsGMGhhjx86s60JOEUv2\npLItIaPW7QM6BPObEVF4e8hTqkI4EqcN+S37Elm56RAHTxgvNA7sGcXT94zBx7v5zxppMpexPTGD\nnSczOJicU2O7m4tiat82TOrdGjdXGZYRwhE5XcinXshi/vdb2X34tKHu5+PJI3eMYHj/aDt1Zj0m\ncxmvLTtG4oX8WreHBXgS08af2wdGyLCMEA7OaX7CtdZs2XeSOV+uM0xLANC/e3ueumc0/r7Ne57z\n5IwCvth8hoS0vFq3R7X0YUqfNvRtHyS3QgrhJJwi5NMz8/j3p6s5lmS8k2Ro306MHdKV67pG2qmz\na2cyl7HvTBaLtyeTUe1C6q+6tvVn2qAI2rfwtXF3Qgh7c+iQ11qzYdcJPvluM3kFlVPiuru5Mmv6\nSEYP7mLH7hpOa82+01msPJBGWnYRuUWmWvcb3zOMWweEy22QQjgxhw357NxC3vlyLXuPnDXURw7o\nzPRJA2jTMtBOnTVcamYhS/akXnbKAYDWQV4M7BDMuJ5h+Ho67P9eIUQ9OWQKLFt/gAXLdlBUXDnn\njLubK8/fP55BvaLs11gDJaTlsXzfucvO4Q4wpnsrRnVtSXiItw07E0I0dQ4V8hnZ+by3cEONO2cG\n9Yri0ekjCQ7wsVNnV+/MpQI2H09nzaELl90n2Ned6YMj6ds+SG6BFELUymFCfseBU/xn4foay/E9\nd9+4ZnNbZFGJmR/3pLLnVGaNp1F/FRnizZDoUMb1DMPVRe6QEUJcWbMP+dz8It76fE2Nsff+3dvz\n+F2jCPJv2mfvZWWafWey+GF3CqmZRZfdz9Pdhd9P6UKHlnKHjBCi/pptyJtMZr6N28tP8fspqDIl\ngaeHO7Omj2jSUxKUmss4nJLD0r2pnLpYcNn92oX60KWNP5P7tCbA292GHQohHEWzDPmk5HRe/vBn\nLmUZn+js0qE1z943llYh/nbq7PKKS83sOJnBpxtP17nvDde1YVCnEMKD5SKqEOLaNKuQ11rz5dLt\nfL9mn6HeVM/ei0vN7EzKJO5AGimZhVfct3t4AJN6tyamtZ9cRBVCWE2zCfnc/CJenx/HgeMphvrE\nYT144NahuLs3jdkTi01m4g6kceJ8HodSak4KVlXXtv6M6d6KXpGBuEuwCyEaQbMI+XMXs5k9d6lh\nMe22LQP5/YMTad/W/uuNaq3ZeTKTnScz2Hs664r7RrX04a7r29GxlZ+NuhNCOLMmHfK/Tir2weIN\nhmkJJg7rwSPThtt1kq28IhM7TmawYMuZK+4X6OPOsM6hDI4OpW2Ql0wMJoSwqSYb8geOpzDv+y2c\nTr1kqSngwduHMWVkL5v3o7UmNauINYfS2HA0/Yr7erq74Ofpxszr23Fde8daNlAI0bw0uZC/lJXH\n21+srbGYh6urCy8+NJH+PdrbrJeUjEKOpOZw7FxuncMwAMNiQhnauQUxrf3kjF0I0SQ0mZDXWrP7\n8Ble/uDnGtvatQnhL0/eSIBf491SmFdkIiEtj8QLeexIzLjsE6fVDegQzJTr2hAZ4i3BLoRocppE\nyKdn5vHOlzXP3gH++tRN9Ihua/XPzCsysSspg03H06/4QFJtpg2KoEdEABEhTftpWiGEsHvIb913\nkre/WGtYrcnD3Y3bJ/Tj9vF9rXJ2rLUmp9DE/rNZ7EjMIKuglHNZl59CoLqY1n4M6hRCv6hgefJU\nCNGs1CvklVKTgDcBF+BjrfWrtezzNjAZyAfu11rvq75PVSaTmY+/28yqzYcN9dhBXZg5ZSAtgq/+\nFsPiUjNpOcWczy4iJaOQghITSRfzr+pMfWTXFrQJ8sbP040BHYPl/nUhRLNWZ8grpVyAOcBYIBXY\nqZT6UWt9tMo+k4FOWuvOSqnBwHvAkNrez2Qys2rLYb5euccwY6Sfjyez7hzJ0Os6XbGf4lIz6Xkl\npGQWknQhnz2nMnFzUaTlFF/x62oT7OtOv6hgurcNILq1n00X2YiPjyc2NtZmn9eUybGoJMeikhwL\n66hPqg0CTmitTwMopRYCNwNHq+xzM/AZgNZ6u1IqUCkVprVOq/5mz7y8iPPpxidBO0W25Ln7x+Ph\n5cnJC3mkZRdTVGpme2IGl/KKySooxcfDlfxicwP/Myv1igyke3gA/aOCCfHzuOb3ayj5Bq4kx6KS\nHItKciysoz4hHw5Uncc3mfLgv9I+KRW1GiG/9xJoFYgZF0wubnSNDCUv2J8//XD8ik1cbcAH+7rT\nLtSH69oH0aWNPy39PeXuFyGE07H5hdcC5QlA69AAWrcMwM3VlRKTbtB7+Xi4Eh7iTffwALq08Seq\nha8sWi2EEFUora8csEqpIcBsrfWkitcvAbrqxVel1HvAOq31oorXR4FR1YdrlFINS3MhhHByWusG\nDUXU50x+JxCtlGoPnANmADOr7bMEeAJYVPFLIau28fiGNimEEKJh6gx5rbVZKfUksIrKWyiPKKVm\nlW/WH2itlyulpiilEii/hfKBxm1bCCFEfdQ5XCOEEKL5apSrlEqpSUqpo0qp40qpFy+zz9tKqRNK\nqX1Kqesao4+moK5joZS6Syn1S8WfTUop20+xaSP1+b6o2G+gUqpUKXWbLfuzpXr+jMQqpfYqpQ4q\npdbZukdbqcfPSIBSaklFVhxQSt1vhzYbnVLqY6VUmlJq/xX2ufrc1Fpb9Q/lvzgSgPaAO7AP6Fpt\nn8nAsoq/Dwa2WbuPpvCnnsdiCBBY8fdJznwsquy3BvgJuM3efdvx+yIQOASEV7xuYe++7Xgs/gt4\n+dfjAFwC3OzdeyMci+HAdcD+y2xvUG42xpm85eEprXUp8OvDU1UZHp4CApVSYY3Qi73VeSy01tu0\n1tkVL7dR/nyBI6rP9wXAU8A3wAVbNmdj9TkWdwHfaq1TALTWV17EoPmqz7HQgH/F3/2BS1prEw5G\na70JyLzCLg3KzcYI+doenqoeXJd7eMrR1OdYVPUwUHOuZcdQ57FQSrUFbtFa/4fyNWIcVX2+L2KA\nEKXUOqXUTqXUvTbrzrbqcyzmAN2VUqnAL8AzNuqtqWlQbtp9FkpRTik1mvK7kobbuxc7ehOoOibr\nyEFfFzegHzAG8AW2KqW2aq0T7NuWXUwE9mqtxyilOgFxSqneWuu8ur5QNE7IpwDtqryOqKhV3yey\njn0cQX2OBUqp3sAHwCSt9ZX+udac1edYDAAWqvL5J1oAk5VSpVrrJTbq0VbqcyySgXStdRFQpJTa\nAPShfPzakdTnWDwAvAygtU5USiUBXYFdNumw6WhQbjbGcI3l4SmllAflD09V/yFdAtwHlidqa314\nygHUeSyUUu2Ab4F7tdaJdujRVuo8FlrrjhV/OlA+Lv+4AwY81O9n5EdguFLKVSnlQ/mFtiM27tMW\n6nMsTgPjACrGoGOAkzbt0nYUl/8XbINy0+pn8loenrKoz7EA/gSEAO9WnMGWaq2rTwDX7NXzWBi+\nxOZN2kg9f0aOKqVWAvsBM/CB1vrwFd62Warn98XfgflVbi18QWudYaeWG41SagEQC4Qqpc4AfwY8\nuMbclIehhBDCgcmUjUII4cAk5IUQwoFJyAshhAOTkBdCCAcmIS+EEA5MQl4IIRyYhLwQQjgwCXkh\nhHBg/w8HVF7m3N08yQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dist_s = thinkbayes2.Beta(1, 1)\n", "dist_t = thinkbayes2.Beta(1.75, 1)\n", "n = 10000\n", "ss = dist_s.Sample(n)\n", "ts = dist_t.Sample(n)\n", "thinkplot.Cdf(thinkbayes2.Cdf(ss))\n", "thinkplot.Cdf(thinkbayes2.Cdf(ts))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.097733895049857883" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "posts = [update(prior, s, t) for s, t in zip(ss, ts)]\n", "np.array(posts).mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now here's a version that simulates worlds where $x$ is known and $s$ and $t$ are drawn from uniform distributions. For each $s$-$t$ pair, we generate one patient with a fever and compute the probability that they are sick." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def prob_sick(x, s, t):\n", " return x * s / (x * s + (1-x) * t)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "dist_s = thinkbayes2.Beta(1, 1)\n", "dist_t = thinkbayes2.Beta(1, 1)\n", "n = 10000\n", "ss = dist_s.Sample(n)\n", "ts = dist_t.Sample(n)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = 0.1\n", "probs = [prob_sick(x, s, t) for s, t in zip(ss, ts)]" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.15306534465160232" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.array(probs).mean()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cohort = np.random.random(len(probs)) < probs" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.1439" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cohort.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "April 6, 2016\n", "\n", "Suppose \n", "* t is known to be 0.5\n", "* x is known to be 0.1\n", "* s is equally likely to be 0.2 or 0.8, but we don't know which\n", "\n", "If we take the average value of s and compute p = p(sick|fever), we would consider p to be the known quantity 0.1 " ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.1" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s = 0.5\n", "t = 0.5\n", "x = 0.1\n", "p = x * s / (x * s + (1-x) * t)\n", "p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we propagate the uncertainty about s through the calculation, we consider p to be either p1 or p2, but we don't know which." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.06493506493506493, 0.06493506493506493)" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 0.8\n", "p1 = x * s / (x * s + (1-x) * t)\n", "p1, 5/77" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.21739130434782605, 0.21739130434782608)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 0.2\n", "p2 = x * s / (x * s + (1-x) * t)\n", "p2, 5/23" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we were asked to make a prediction about a single patient, we would average the two possible values of p." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.1411631846414455" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(p1 + p2) / 2" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(-0.47000362924573613, 0.91629073187415466)" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def logodds(p):\n", " return np.log(p / (1-p))\n", "\n", "logodds(p1) - logodds(0.1), logodds(p2) - logodds(0.1)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.39156220293917254" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_mix = (p1 + p2) / 2\n", "logodds(p_mix) - logodds(0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So let's simulate a series of patients by drawing a random value of y, computing p, and then tossing coins with probability p." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import random\n", "\n", "def simulate_patient():\n", " s = 0.5\n", " t = random.choice([0.2, 0.8])\n", " x = 0.1\n", " p = x * s / (x * s + (1-x) * t)\n", " return random.random() < p\n", "\n", "simulate_patient()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this simulation, the fraction of patients with fever who turn out to be sick is close to 0.141" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.1402" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "patients = [simulate_patient() for _ in range(10000)]\n", "sum(patients) / len(patients)" ] }, { "cell_type": "code", "execution_count": 153, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = 0.1\n", "s = 0.5\n", "t1 = 0.2\n", "t2 = 0.8" ] }, { "cell_type": "code", "execution_count": 154, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pandas as pd\n", "d1 = dict(feverbad=0.5, fevergood=0.5)\n", "d2 = dict(sick=x, notsick=1-x)\n", "d3 = dict(fever='t', notfever='1-t')\n", "\n", "iterables = [d1.keys(), d2.keys(), d3.keys()]\n", "\n", "index = pd.MultiIndex.from_product(iterables, names=['first', 'second', 'third'])\n", "df = pd.DataFrame(np.zeros(8), index=index, columns=['prob'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "code", "execution_count": 155, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
prob
firstsecondthird
fevergoodnotsickfever0.360
notfever0.090
sickfever0.025
notfever0.025
feverbadnotsickfever0.090
notfever0.360
sickfever0.025
notfever0.025
\n", "
" ], "text/plain": [ " prob\n", "first second third \n", "fevergood notsick fever 0.360\n", " notfever 0.090\n", " sick fever 0.025\n", " notfever 0.025\n", "feverbad notsick fever 0.090\n", " notfever 0.360\n", " sick fever 0.025\n", " notfever 0.025" ] }, "execution_count": 155, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t_map = dict(fevergood=t2, feverbad=t1)\n", "\n", "for label1, p1 in d1.items():\n", " t = t_map[label1]\n", " for label2, p2 in d2.items():\n", " for label3, p3 in d3.items():\n", " if label2 == 'sick':\n", " p = p1 * p2 * 0.5\n", " else:\n", " p = p1 * p2 * eval(p3)\n", "\n", " df.prob[label1, label2, label3] = p\n", " \n", "df " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If there are two kinds of people, some more fever prone than others, and we don't know which kind of patient we're dealing with, P(sick | fever)" ] }, { "cell_type": "code", "execution_count": 157, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.09999999999999998" ] }, "execution_count": 157, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.prob[:, 'sick', 'fever'].sum() / df.prob[:, :, 'fever'].sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we know fevergood, P(sick|fever) is" ] }, { "cell_type": "code", "execution_count": 158, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.064935064935064929" ] }, "execution_count": 158, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_sick_fevergood = df.prob['fevergood', 'sick', 'fever'].sum() / df.prob['fevergood', :, 'fever'].sum()\n", "p_sick_fevergood" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we know feverbad, P(sick|fever) is" ] }, { "cell_type": "code", "execution_count": 149, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.21739130434782605" ] }, "execution_count": 149, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_sick_feverbad = df.prob['feverbad', 'sick', 'fever'].sum() / df.prob['feverbad', :, 'fever'].sum()\n", "p_sick_feverbad" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we think there's a 50-50 chance of feverbad" ] }, { "cell_type": "code", "execution_count": 159, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.1411631846414455" ] }, "execution_count": 159, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(p_sick_fevergood + p_sick_feverbad) / 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If fevergood, here's the fraction of all patients with fever:" ] }, { "cell_type": "code", "execution_count": 143, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.38500000000000006" ] }, "execution_count": 143, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_fever_fevergood = df.prob['fevergood', :, 'fever'].sum()\n", "p_fever_fevergood" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If fever bad, here's the fraction of all patients with fever" ] }, { "cell_type": "code", "execution_count": 144, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.11500000000000002" ] }, "execution_count": 144, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_fever_feverbad = df.prob['feverbad', :, 'fever'].sum()\n", "p_fever_feverbad" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So if we started out thinking there's a 50-50 chance of feverbad, we should now think feverbad is less likely" ] }, { "cell_type": "code", "execution_count": 146, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.22999999999999998" ] }, "execution_count": 146, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_feverbad_fever = p_fever_feverbad / (p_fever_feverbad + p_fever_fevergood)\n", "p_feverbad_fever" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And if we compute the weighted sum of the two possible worlds" ] }, { "cell_type": "code", "execution_count": 150, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.099999999999999978" ] }, "execution_count": 150, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p_feverbad_fever * p_sick_feverbad + (1-p_feverbad_fever) * p_sick_fevergood" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }