{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evidence of gluten sensitivity\n",
    "------------------------------\n",
    "\n",
    "This notebook contains an exploration of results from this paper:\n",
    "\n",
    "http://onlinelibrary.wiley.com/doi/10.1111/apt.13372/epdf\n",
    "\n",
    "which reports results of a study of gluten sensitivity.\n",
    "\n",
    "Copyright 2015 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",
    "from scipy import stats\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following class defines the function that computes the likelihood of the data given the hypothetical number of subjects who are gluten sensitive.\n",
    "\n",
    "I assume that a subject who is gluten sensitive has a 95% chance of successfully identifying gluten flour based on resumption of symptoms, and that a subject who is not gluten sensitive has only a 40% chance of identifying gluten flour (and a 20% chance of detecting no difference).\n",
    "\n",
    "The function works by computing the PMF of the number of gluten identifications conditioned on gs, and then selecting the actual number of identifications from the PMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Gluten(thinkbayes2.Suite):\n",
    "    \n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"\n",
    "        hypo: int hypothetical # who are gluten sensitive\n",
    "        data: (int, int) # who correctly identify gluten, others\n",
    "        \"\"\"\n",
    "        gs = hypo\n",
    "        yes, no = data\n",
    "        n = yes + no\n",
    "        ngs = n - gs\n",
    "        \n",
    "        pmf1 = thinkbayes2.MakeBinomialPmf(gs, 0.95)\n",
    "        pmf2 = thinkbayes2.MakeBinomialPmf(ngs, 0.1)\n",
    "        pmf = pmf1 + pmf2\n",
    "        return pmf[yes]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior is uniform from 0 to 35; that is, there are equally likely to be any number of GS subjects in the study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEACAYAAACpoOGTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFWdJREFUeJzt3X2MVWeBx/HvlRd33VFx0y6VgWSo0Fi6q9K1MIa23Lq1\nUkxgN6ZSElIXky2xImw1CjQaZptNLKaJLNttIes0oqYlDd02ky4vtgl3o4mdlpaX8jKVQckCFTRu\nUWGbFMLsH88zzOH2Gc6Z6YWZs/P9JDf3nOc855zn3gvnN8/znDsDkiRJkiRJkiRJkiRJkiRJUinN\nAbqAQ8CKfuqsi9v3ANNj2SRgB7Af2Acsy9T/OPBzYC/QAbw/lrcAbwG74uPRBr0GSVKDjQK6CRfu\nMcBu4Pq6OnOBLXF5JvBiXL4G+ERcbgJeBz4a118GbonLi4EH43IL8FqjGi9Junw+BWzLrK+Mj6z1\nwILMehcwPnGsZ4G/icunMuWTCL0MMCAkadh4T872ZuBoZv1YLMurM7GuTgth6Kkzru8H5sfluwgh\n0WsyYXipBtyc0z5J0mWSFxA9BY9TucR+TcBmYDlwOpZ9CbgP2Bm3vx3L3yCExXTga8AT9M1PSJKu\noNE5249z8U/3kwg9hEvVmRjLIMxbPA38mDDE1Ot14LNx+Trgc3H5bfrC4lXgMDA1Ll/wkY98pOfw\n4cM5TZck1TkMTGnUwUbHA7YAY8mfpG6lb5K6AvwQ+F7iuFfH5/fEOn8f168iTIwDXEsIo3GJ/XvK\nbPXq1UPdhHelzO0vc9t7emz/UCt7+yk+KgTk9yDOAUuB7fHC3Q4cBJbE7RsI4TCXcLfTGcJdSQCz\ngEWEW1l3xbJVhEnvhcBXYtnTwA/i8q2EO5rOAufjebIT2pKkKyQvIAC2xkfWhrr1pYn9fkb/cxzr\n4qPef8SHJGmI5U1S6zKoVqtD3YR3pcztL3PbwfYPtbK3f6Dq7z4qizicJkkqqlKpwACu+/YgJElJ\nBoQkKcmAkCQlGRCSpCQDQpKUZEBIkpIMCElSkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRA\nSJKSDAhJUpIBIUlKMiAkSUkGhCQpyYCQJCUZEJKkJANCkpRUJCDmAF3AIWBFP3XWxe17gOmxbBKw\nA9gP7AOWZep/HPg5sBfoAN6f2bYqHqsLuKPIi5AkNV5eQIwCHiGExDRgIXB9XZ25wBRgKnAv8Fgs\nPwvcD9wAtAJfAT4at30f+CbwMeAZ4BuxfBqwID7PAR4t0EZJ0mWQd/GdAXQDRwgX/E3A/Lo684CN\ncbkTGAeMB04Au2P5aeAg0BzXpwI/jcsvAJ+Py/OBJ+O5jsRzzyj+ciRJjZIXEM3A0cz6Mfou8peq\nM7GuTgth6Kkzru+nL2juIgxHAUyI+1/qfJKkK2B0zvaegsepXGK/JmAzsJzQkwD4EmHe4tuEOYi3\nB9qGtra2C8vVapVqtVqwqZI0MtRqNWq12qD3r7+w12sF2gjzARAmkM8DazJ11gM1wvAThMnl2cBJ\nYAzwHLAVWNvPOa4DfgTMBFbGsofi8zZgNX09j149PT1Fs0uSBFCpVCD/un9B3hDTTsJ8QQswljCB\n3FFXpwO4Jy63AqcI4VAB2oEDvDMcrs6c/1v0TWx3AHfHc02O536p6IuRJDVO3hDTOWApsJ1wR1M7\nYbJ5Sdy+AdhCuJOpGzgDLI7bZgGLCLey7oplqwi9goWEu5oAngZ+EJcPAE/F53PAfRQf5pIkNVDh\nrsYw4xCTJA1Qo4eYJEkjlAEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRASJKSDAhJUpIBIUlK\nMiAkSUkGhCQpyYCQJCUZEJKkJANCkpRkQEiSkgwISVKSASFJSjIgJElJBoQkKcmAkCQlGRCSpCQD\nQpKUVCQg5gBdwCFgRT911sXte4DpsWwSsAPYD+wDlmXqzwBeAnYBLwM3xfIW4K1Yvgt4tNjLkCQ1\n2uic7aOAR4DbgeOEi3kHcDBTZy4wBZgKzAQeA1qBs8D9wG6gCXgF+AkhbL4LfBvYDtwZ12+Lx+um\nL2QkSUMkrwcxg3DBPkK44G8C5tfVmQdsjMudwDhgPHCCEA4Apwmh0hzXfw18MC6PI4SPJGkYyetB\nNANHM+vHCL2EvDoTgZOZshZCr6Azrq8EfgY8TAipT2XqTiYML/0e+FasJ0m6wvICoqfgcSqX2K8J\n2AwsJ/QkANoJcxLPAHcBjwOfAd4gzF28CdwIPAvcAPyx/oRtbW0XlqvVKtVqtWBTJWlkqNVq1Gq1\nQe9ff2Gv1wq0ESaqAVYB54E1mTrrgRph+AnCHMNsQg9iDPAcsBVYm9nnD8AHMm04Rd+QU9YO4OvA\nq3XlPT09RbNLkgRQqVQg/7p/Qd4cxE7C5HMLMBZYQJikzuoA7onLrYSL/cnYiHbgABeHA4R5jdlx\n+dPAL+LyVYSJcYBr47l/WeiVSJIaKm+I6RywlHC30SjCBf8gsCRu3wBsIdzJ1A2cARbHbbOARcBe\nwpwChB7INuBe4N+A9xJua703br8VeJAwIX4+nufUYF+cJGnwCnc1hhmHmCRpgBo9xCRJGqEMCElS\nkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRASJKSDAhJUpIBIUlKMiAkSUkGhCQpyYCQJCUZ\nEJKkJANCkpRkQEiSkgwISVKSASFJSjIgJElJBoQkKcmAkCQlGRCSpKQiATEH6AIOASv6qbMubt8D\nTI9lk4AdwH5gH7AsU38G8BKwC3gZuCmzbVU8VhdwR5EXIUlqvErO9lHA68DtwHHCxXwhcDBTZy6w\nND7PBP4FaAWuiY/dQBPwCjCfcOGvAd8BtgN3At8EbgOmAU8QAqMZeAG4Djhf166enp6eAb5USRrZ\nKpUK5F/3L8jrQcwAuoEjwFlgE+EinzUP2BiXO4FxwHjgBCEcAE4TQqU5rv8a+GBcHkcIH+Kxn4zn\nOhLPPaPoi5EkNc7onO3NwNHM+jFCLyGvzkTgZKashTD01BnXVwI/Ax4mhNSnYvkE4MW6YzUjSbri\n8gKi6DhOfZclu18TsBlYTuhJALQT5iSeAe4CHgc+M5A2tLW1XViuVqtUq9WCTZWkkaFWq1Gr1Qa9\nf95YVCvQRpiohjCBfB5Yk6mznjCnsCmudwGzCT2IMcBzwFZgbWafPwAfyLThFGHIaWUseyg+bwNW\n09fz6OUchCQNUKPnIHYCUwlDRGOBBUBHXZ0O4J643Eq42J+MjWgHDnBxOECYW5gdlz8N/CJzrLvj\nuSbHc79U9MVIkhonb4jpHOEOpe2EO5raCZPNS+L2DcAWwh1M3cAZYHHcNgtYBOwl3M4KoQeyDbgX\n+DfgvcBbcR1CmDwVn88B91F8mEuS1ECFuxrDjENMkjRAjR5ikiSNUAaEJCnJgJAkJRkQkqQkA0KS\nlGRASJKSDAhJUpIBIUlKyvsm9bA163MPDHUTJOn/NXsQkqQkA0KSlOTvYpKkEcLfxSRJaggDQpKU\nZEBIkpIMCElSkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRASJKSigTEHKALOASs6KfOurh9\nDzA9lk0CdgD7gX3Askz9TcCu+PhVfAZoAd7KbHu02MuQJDVa3t+DGAU8AtwOHAdeBjqAg5k6c4Ep\nwFRgJvAY0AqcBe4HdgNNwCvA83HfuzP7Pwycyqx30xcykqQhkteDmEG4YB8hXPA3AfPr6swDNsbl\nTmAcMB44QQgHgNOEYJhQt28F+ALw5MCbLkm6nPICohk4mlk/Fsvy6kysq9NC6BV01pXfApwEDmfK\nJhOGl2rAzTntkyRdJnlDTEX/6EL97xfP7tcEbAaWE3oSWQuBJzLrbxDmLt4EbgSeBW4A/lh/wra2\ntgvL1WqVarVasKmSNDLUajVqtdqg98/7wxGtQBthohpgFXAeWJOps57w0/6muN4FzCb0DMYAzwFb\ngbV1xx5N6G3cSAiGlB3A14FX68r9g0GSNECN/oNBOwmTzy3AWGABYZI6qwO4Jy63EiacT8ZGtAMH\neGc4QJj4PsjF4XAVYWIc4Np47l/mvwxJUqPlDTGdA5YC2wkX7nbCRX1J3L4B2EK4k6kbOAMsjttm\nAYuAvfTdxvoAoTcBIWzqJ6dvBR4kTIifj+c5hSTpivNvUkvSCOHfpJYkNYQBIUlKMiAkSUkGhCQp\nyYCQJCUZEJKkJANCkpRkQEiSkgwISVKSASFJSjIgJElJBoQkKcmAkCQlGRCSpCQDQpKUZEBIkpIM\nCElSkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlFQkIOYAXcAhYEU/ddbF7XuA6bFsErAD2A/s\nA5Zl6m8CdsXHr+Jzr1XxWF3AHUVehCSp8UbnbB8FPALcDhwHXgY6gIOZOnOBKcBUYCbwGNAKnAXu\nB3YDTcArwPNx37sz+z8MnIrL04AF8bkZeAG4Djg/mBcnSRq8vB7EDKAbOEK44G8C5tfVmQdsjMud\nwDhgPHCCEA4ApwnBMKFu3wrwBeDJuD4/Lp+N5+yObZAkXWF5AdEMHM2sH4tleXUm1tVpIQw9ddaV\n3wKcBA7H9Qlx/0udT5J0BeQNMfUUPE7lEvs1AZuB5YSeRNZC4InBtKGtre3CcrVapVqtFmimJI0c\ntVqNWq026P3rL+z1WoE2wkQ1hAnk88CaTJ31QI0w/ARhcnk2oWcwBngO2AqsrTv2aEIP4UbgjVi2\nMj4/FJ+3Aat5Z8+jp6enaHZJkgAqlQrkX/cvyBti2kmYfG4BxhImkDvq6nQA98TlVsKE88nYiHbg\nAO8MBwgT3wfpC4feY90dzzU5nvulQq9EktRQeUNM54ClwHbCHU3thIv6krh9A7CFcCdTN3AGWBy3\nzQIWAXvpu431AUJvAkLY9E5O9zoAPBWfzwH3UXyYS5LUQIW7GsOMQ0ySNECNHmKSJI1QBoQkKcmA\nkCQlGRCSpCQDQpKUZEBIkpIMCElSkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRASJKSDAhJ\nUpIBIUlKMiAkSUkGhCQpyYCQJCUZEJKkJANCkpRkQEiSkgwISVJSkYCYA3QBh4AV/dRZF7fvAabH\nsknADmA/sA9YVrfPV4GDcduaWNYCvAXsio9HC7RPknQZjM7ZPgp4BLgdOA68DHQQLuy95gJTgKnA\nTOAxoBU4C9wP7AaagFeA5+O+twHzgI/FeldnjtdNX8hIkoZIXg9iBuGCfYRwId8EzK+rMw/YGJc7\ngXHAeOAEIRwAThOCYUJc/zLwnXhMgN8OqvWSpMsmLyCagaOZ9WOxLK/OxLo6LYReQWdcnwrcCrwI\n1IBPZupOJgwv1YCbc9onSbpM8oaYegoep3KJ/ZqAzcByQk+i97wfIgxF3QQ8BVwLvEGYu3gTuBF4\nFrgB+GP9Cdva2i4sV6tVqtVqwaZK0shQq9Wo1WqD3r/+wl6vFWgjTFQDrALO0zepDLCe8NP+prje\nBcwGTgJjgOeArcDazD5bgYeA/4rr3YT5i9/VnX8H8HXg1brynp6eotklSQKoVCqQf92/IG+IaSdh\nOKgFGAssIExSZ3UA98TlVuAUIRwqQDtwgIvDAULP4NNx+bp47N8BVxEmxiH0KKYCvyz6YiRJjZM3\nxHQOWApsJ1y42wmTzUvi9g3AFsKdTN3AGWBx3DYLWATsJcwpADxA6D08Hh+vAW/TFzC3Ag8SJq/P\nx/OcGuyLkyQNXuGuxjDjEJMkDVCjh5gkSSOUASFJSjIgJElJBoQkKcmAkCQlGRCSpCQDQpKUZEBI\nkpIMCElSkgEhSUoyICRJSQaEJCnJgJAkJRkQkqQkA0KSlGRASJKSDAhJUpIBIUlKMiAkSUkGhCQp\nyYCQJCUZEJKkpCIBMQfoAg4BK/qpsy5u3wNMj2WTgB3AfmAfsKxun68CB+O2NZnyVfFYXcAdBdon\nSboM8gJiFPAIISSmAQuB6+vqzAWmAFOBe4HHYvlZ4H7gBqAV+Epm39uAecDHgL8EHo7l04AF8XkO\n8GiBNpZOrVYb6ia8K2Vuf5nbDrZ/qJW9/QOVd/GdAXQDRwgX/E3A/Lo684CNcbkTGAeMB04Au2P5\naUJvYUJc/zLwnXhMgN/G5/nAk7H8SDz3jOIvpxzK/o+szO0vc9vB9g+1srd/oPICohk4mlk/Fsvy\n6kysq9NCGHrqjOtTgVuBF4Ea8MlYPiHuf6nzSZKugNE523sKHqdyif2agM3AckJPove8HyIMPd0E\nPAVc+y7bIEm6glqBbZn1Vbxzono9cHdmvYswxAQwBtgO/GPdPluB2Zn1buAqYGV89NoGzEy0q5sQ\nHD58+PDho/ijmwYaDRwmDBGNJcwppCapt8TlVsKwEYRexQ+B7yWOuwT4p7h8HfDfcXlaPMdYYHI8\nd33vRJI0TNwJvE5InlWxbEl89Hokbt8D3BjLbgbOEy74u+LjzrhtDPAj4DXgFaCaOdYD8VhdwGcb\n+kokSZIkjSxFvrQ3nB0B9hJ6Uy8NbVMKeRw4Sejp9fpz4HngF8BPCLc1D1ep9rcR7o7r7dXOufLN\nKqy/L5uW5TPor/1tDP/P4E8Id13uBg4QbsuH8rz3/bW/jeH/3g/KKMLQUwthiCo1HzLc/YrwD6ws\nbiHcnpy9wH4X+GZcXgE8dKUbNQCp9q8GvjY0zRmwa4BPxOUmwlDv9ZTnM+iv/WX5DN4Xn0cT5lZv\npjzvPaTbP6D3vkzfUi7ypb0yKNOk+0+BN+vKsl+M3Aj87RVt0cCk2g/l+QxSXzZtpjyfQX/th3J8\nBv8bn8cSfkB9k/K895BuPwzgvS9TQBT50t5w1wO8AOwE/mGI2zJY4wnDNsTn8ZeoO1x9lXBDRTvD\nd4igXgt9XzYt42fQQmh/712OZfgM3kMIuJP0DZWV6b1PtR/K8d4P2OeBf8+sLwL+dYjaMlgfjs9X\nEz64W4awLUW1cPEQTf1P5P9z5ZoyKC1c3P6/IPwEVQH+mfCfZLhrItzt1/vTatk+gybCD0W97S/b\nZ/BBQrDdRvnee+hrf5UBvvdl6kEcJ0x69ZrExb+Wowx+HZ9/CzxDOX/P1EnC2DKEwPvNELZlMH5D\n35eGvs/w/wzGAE8Tbgt/NpaV6TPobf+P6Wt/2T6D3wP/Cfw15Xrve/W2/5MM8L0vU0DsJPwOpxbC\nmNoCoGMoGzRA7wPeH5f/jPCrzF/rv/qw1QF8MS5/kb7/9GXx4czy3zG8P4MK4Se8A8DaTHlZPoP+\n2l+Gz+Aq+oZf/hT4DOGun7K89/21/5pMneH63g9a6kt7ZTGZMKy0m3DLXxna/yTwBvA2Yf5nMeEu\nrBcY/rf5wTvb/yXCt/v3EsZgn2V4jyGnvmw6h/J8Bv19WbYMn8FfAa8S2r4X+EYsL8t731/7y/De\nS5IkSZIkSZIkSZIkSZIkSZIkSZLUGP8H7SnxvybDzrkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc4a5e160d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = 12, 23\n",
    "prior = Gluten(range(0, 35+1))\n",
    "thinkplot.Pdf(prior)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Here's the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.032678923853482846"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior = prior.Copy()\n",
    "posterior.Update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8G3ed//GXfMVX4iRO4ttxmqTN1aRHmiY9wFzdcmxb\njt1SYLkKdIEUfl0oZReWJrCPLeePH8cWSilQyi4t+4CyBVp6sBh6pGnS5m7S3M1px3ZiJ74iW9bv\nj+9IGimybNkezYz8fj4eenhmNCN/Ijnz0fcGERERERERERERERERERERERERERERGaFrgV3AHuCO\nJM+/F9gCbAWeBZbanjtoHd8EvOBolCIi4opcYC/QAOQDm4GFCeesAsqs7WuB523PHQCmOxuiiIgM\nJ8fB116BSRQHgX7gQeD6hHPWAZ3W9nqgNuH5gIPxiYjICDiZKGqAw7b9I9axodwMPGrbDwNPARuB\nj457dCIiMiJ5Dr52OI1zXwd8GLjSduxK4DgwE3gS09bx9LhFJyIiI+JkojgK1Nn26zClikRLgXsx\nbRSnbMePWz9bgYcxVVlxiWLu3Lnhffv2jVe8IiITxT5gnttBgElC+zCN2QUkb8yux7RjrEw4XgxM\ntrZLMD2irknyO8J+duedd7odwpgofnf5OX4/xx4O+z9+0qvxcbREMQCsBh7H9IC6D9gJ3GI9fw/w\nJWAa8APrWD+m5FAJ/MYW438CTzgYq4iIDMHJRAHwmPWwu8e2/RHrkWg/cJFTQYmIyMg52etJhtHY\n2Oh2CGOi+N3l5/j9HDv4P/50+X2cglXdJiIiIxUIBCCN+79KFCIikpIShYiIpKREISIiKSlRiIhI\nSkoUIiKSkhKFiIikpEQhIiIpKVGIiEhKShQiIpKSEoWIiKSkRCEiIikpUYiISEpKFCIikpIShYiI\npKREISIiKSlRiIhISkoUIiKSkhKFiIikpEQhIiIpKVGIiEhKShQiIpKSEoWIiKSkRCEiIikpUYiI\nSEpKFCIikpIShYiIpKREISIiKSlRiIhISkoUIiKSkhKFiIikpEQhIiIpKVHIqAwOhgmHw26HISIZ\noEQhadu5u40PrX6E277wBPsOnnI7HBFxmBKFpO3Xj+zkTNdZDh3p5Atf+V+anj3odkgi4iCnE8W1\nwC5gD3BHkuffC2wBtgLPAkvTuFZc0D8QYvuu1uh+sD/Ed+95gft+sYmBgUEXIxMRpziZKHKB72Nu\n+IuAm4CFCefsB16DSRBfAX6UxrXigj37ThIMDpxz/A9P7GHt1/9CR2efC1GJiJOcTBQrgL3AQaAf\neBC4PuGcdUCntb0eqE3jWnHB1h0t0e2rVtazcnltdH/HrlY+t+Yp9u4/6UZoIuIQJxNFDXDYtn/E\nOjaUm4FHR3mtZMjWHSei25dfWsPtt67iPe9aAgQAaGvv4Qv/9mf+/PRBdwIUkXGX5+Brp9N38nXA\nh4Er0712zZo10e3GxkYaGxvT+LWSjp7efnbva4/uX7hoFoFAgHddt4g5s6fx7bufp6e3n/6BEN+7\n9wX2HTzJB2+6iLw89ZkQcVNTUxNNTU2jvj4wfqGcYyWwBtPOAPDPwCDwtYTzlgK/sc7bm+a1YfXl\nz5wNm45x17efAeC82dP45lfeFPf88eYzfPU7z3H4aGf02KILZvLZ1auYWlaY0VhFZGiBQADSuP87\n+VVvIzAfaAAKgBuBRxLOqcckifcRSxIjvVYyzN4+ceGiWec8X1U5ma/e+XpWXRZrt3j5lVbuvm9j\nRuITEWc4mSgGgNXA48DLwEPATuAW6wHwJWAa8ANgE/DCMNeKi7a9HGufWLq4Iuk5RYX5fHb1Kt73\ndxdGj7245TjBYMjx+ETEGU5WPWWCqp4y5FRHHzd/yhTq8vJy+PndN1BYmLqJ6+OffZSWE10AfOsr\n1zBn9lTH4xSR4Xmp6kmyiL3a6YJ55cMmCYC6minRbXu7hYj4ixKFjMhIqp0SxSeK0+Mek4hkhhKF\nDCscDrP15ViJYjSJ4tARlShE/EqJQoZ1rPkMbe09ABQX5TNvzrQRXVdfUxbdPnxMJQoRv1KikGHZ\nq50WL5hJbu7I/mxqq6dEGs1obulWzycRn1KikGHZp+0YabUTQEFBLhWzSqy9MEdUqhDxJSUKSSkU\nGhxVQ3ZEXbWtnUI9n0R8SYlCUtp/8BTdPUEAppUVUVs9Oa3r62tt7RTq+STiS0oUklJcaWLJrGib\nw0ipi6yI/ylRSEqjbZ+IqKtRiULE75QoZEjBYIidu9ui+0uTTAQ4nJqqydFSSMuJbs6ePXd1PBHx\nNiUKGdLOPW30D5gurTVVUyifXpz2axQU5FJVUWrthTWeQsSHlChkSPb5nZYtSb80EaF2ChF/U6KQ\nIdnbJy5clH77RITaKUT8TYlCkurqCrLvwCnATEm8eMHMUb9Wfa2tRKE5n0R8R4lCktq28wSRpcvn\nzZlOaUnBqF9LVU8i/qZEIUnZ2yeWLh59+wSYJVJzckzPpxNt3fT29Y/p9UQks5QoJKmxTNuRqCA/\nl6qK2IjuIypViPiKEoWco7Wtm2PNZwAoKMjjgnnlY37NuHYKJQoRX1GikHPYezstPH8GBQW5Y35N\ntVOI+JcShZwjvtppbO0TEeoiK+JfShQSZ7TLng4nbllUTTcu4itKFBLn8NHTdHT2AVBaUsCc+qnj\n8rpVlaXRlfHa2nvo6VXPJxG/UKKQOPZusRcumhXt1jpW+Xm5VFeWRvfV80nEP5QoJM5YpxVPpbZa\n7RQifqREIVEDA4Ns3xVLFMvGOVGonULEn5QoJOrAq6fo6zPrRcwsL6FiVsm4vr7GUoj4kxKFRO1/\ntSO6fcH88rSXPR2OusiK+JMShUQdsCWKObPHp7eTXVVFKXl55k+u/WQP3T3Bcf8dIjL+lCgk6uAh\nW6IYp26xdnl5OVRXxuZ8UqlCxB+UKASAUGiQg4djDcwNDiQKgPpaVT+J+I0ShQBwvKWLYNA0ZE+b\nWsTUskJHfo/mfBLxHyUKAZyvdoqI6yKr1e5EfEGJQoD4Hk9ONGRHqEQh4j9KFAI43+MponJWKfl5\nZtryUx29dHWr55OI1zmdKK4FdgF7gDuSPL8AWAf0AZ9JeO4gsBXYBLzgXIgSDoc58Oqp6L6TVU+5\nuTnUVKvnk4ifOJkocoHvY5LFIuAmYGHCOe3ArcA3k1wfBhqBi4EVjkUpnOrs4/SZswAUFuZRMat0\nmCvGRu0UIv7iZKJYAezFlAz6gQeB6xPOaQU2Ws8nM75DgyUpe7VTQ/3UcZsxdij2RKFZZEW8z8lE\nUQMctu0fsY6NVBh4CpNIPjqOcUkCe7XTeQ62T0TU26by0OSAIt6X5+Brh8d4/ZXAcWAm8CSmrePp\nsQYl5zrwqvMD7ezU80nEX5xMFEeBOtt+HaZUMVLHrZ+twMOYqqxzEsWaNWui242NjTQ2NqYZphw4\nFCtRZCJRzJpZQkF+LsH+EB2dfZzpOsvk0kmO/16RiaqpqYmmpqZRX+9kotgIzAcagGPAjZgG7WQS\nK8WLMY3hZ4AS4BpgbbIL7YlC0tfT209zSxdgeiTZq4WckpubQ231FPZbVV6Hjpxm8YKZjv9ekYkq\n8Uv02rVJb6dDcrKNYgBYDTwOvAw8BOwEbrEeAJWYdozbgC8Ch4BS6/jTwGZgPfB74AkHY52wXrXN\n71RTNZmCgtyM/N746ie1U4h4mZMlCoDHrIfdPbbtZuKrpyK6gIucCkpi4sZPZKAhO0LtFCL+oZHZ\nE5y9a+x5s6dl7PfaFzHSWAoRb1OimOAOHIofQ5EpdVoWVcQ3lCgmsIGBQQ4fid2kG+qdb8iOmDWj\nhIICU/N5+sxZOk/3Zex3i0h6lCgmsKPHT9M/EAJgRnlxRruo5uQEqNWcTyK+oEQxgR3I0BoUQ9Gc\nTyL+oEQxgcVPLZ65huwILYsq4g9KFBNYptagGIq6yIr4gxLFBBUOh+OXP3UlUcSXKMLhsU4PJiJO\nUKKYoNrae6Kry5UUFzCzvDjjMcwsL2bSJNPz6UzXWTo6z2Y8BhEZnhLFBBW/BkUZgUDml/7IyQlo\nKg8RH0iVKH5m2/6Aw3FIhu13uSE7Qu0UIt6XKlEss23/H6cDkcxyu30iorY6liiOt5xxLQ4RGZqq\nniaogy6PoYioqoitz33cmu5cRLwl1eyxtcB3MWtF1Ni2waxe9ylnQxOndHUFOdHWDUBeXk7ct/pM\nq66Mjc4+1qxEIeJFqRLF7ZiEEABeTHhO/Rh97ODhWGmivqaMvDz3CpYVM0ui261t3QwMDLoaj4ic\nK1Wi+FmmgpDMcnugnd2kSXlMn1bEyVO9DA6GOdHWHVfKEBH3pUoUvyNWokgUBq5zJCJx3P5XM7tG\n9nCqKko5eaoXgOaWLiUKEY9JlShWAkeAX2KWI4X4NgrxqYOHYuMV3C5RAFRVTGbHrlZADdoiXpQq\nUVQBbwJush5/wCSNHRmISxwS7A9x5Jh9DQoPJIpKe88ndZEV8ZpUrYYDmPWu348pXewF/gKszkBc\n4pAjR08TCg0CUDGrlOKifJcjgspZsUTR3NLtYiQikkyqEgVAIfBW4N1AA/Ad4GGHYxIHxa+R7X5p\nAuLHUhxTiULEc1IligeAxcCjwJeBbRmJSBx14FCsIdsL7RMQX6JobetRF1kRj0n1v/G9wPnAp4Hn\ngDO2hybl8akDr8Yasr3QPgFQWJjHtKlFAIRCg7S197gckYjYpUoUOUApMDnJw72hvDJqg4Nhz0zd\nkaja1qB9rFnVTyJekipRFAG3Ad8HbmH49gzxuBOt3fT29QMwZfIkpk8rcjmimLgG7RPqIiviJakS\nxf3ApcB24C3AtzISkTjmwCH7GhRTXVmDYihxkwNqzicRT0lVSlgIXGht/xjY4Hw44iQv9niKqLKN\nxtagOxFvGW4cRbJt8akDHpu6w05VTyLelapEsRTTwymiyLYfRg3avnPgkDdWtUvGXvXU0tpNKDRI\nbq66yIp4Qar/ibnE93TKQ72efKvzdF904r2Cgry4XkZeUFiYx7SyWBfZVnWRFfEMfWWbIOzdYmfX\nlnny23rcnE9q0BbxDO/dLcQR+z20BsVQ4pdF1VgKEa9Qopgg4gbaeTRR2Bu01fNJxDuUKCaIvfu9\n2+Mpwl711KxEIeIZShQTQOfpvmhVTl5eji9KFM0nNN24iFcoUUwAu/eejG7PbZhGQX6ui9EMLa5E\ncaIrum6GiLjL6URxLbAL2APckeT5BcA6oA/4TJrXygi9src9un3+vHIXI0mtqDBfXWRFPMjJRJGL\nmVDwWmARZjnVhQnntAO3At8cxbUyQq/sbYtuX+DhRAFQUVES3VY7hYg3OJkoVmCWTz0I9AMPAtcn\nnNMKbLSeT/daGYGBgUH22BqyF8yf4WI0w6uu0JxPIl7jZKKoAQ7b9o9Yx5y+VmwOHu4gGDRTdc0o\nL/bU1OLJVNpLFJrzScQTnEwUYZeuFZvdtvYJr5cmAKrsJQqNzhbxBCcXIzoK1Nn26zAlg3G9ds2a\nNdHtxsZGGhsb04kx6+3a44+G7Ij40dlKFCLjoampiaamplFf72Si2AjMBxqAY8CNmEbpZBJX0Bnx\ntfZEIeeylygumDvdxUhGpjJuFtkuzSIrMg4Sv0SvXbs2reudTBQDwGrgcUwvpvuAnZhlVQHuASox\nCyJNAQaBT2N6OXUNca2k4WRHLyfazMC1goI8z00tnkxxUT5Tywrp6OxjYGCQtpO9VMwsGf5CEXGM\n0+tgP2Y97O6xbTcTX8U03LWShlds1U5zG6aRl+ePb+aVFaV0dPYBpkFbiULEXf64c8ioxFU7zfd+\n+0RE/PrZmkVWxG1KFFnMPiJ7gQ8asiPUoC3iLUoUWap/IMS+A7GBdn7o8RRRpUF3Ip6iRJGlDhzs\noH8gBEDFrFKmlhW6HNHI2Xs+aRoPEfcpUWQpv1Y7QXzVk2aRFXGfEkWWsg+081NDNpgusmVTTAlo\nYGCQ9lO9LkckMrEpUWQp+4yx58/1V6IANWiLeIkSRRZqa+/hpPUtfNKkPGbXlbkcUfrUTiHiHUoU\nWWjXHntpYrovp8CIH0uhRCHiJv/dQWRY9obsC+Z5f8bYZOKrnjToTsRNShRZyL5Gtt8asiMqE3o+\niYh7lCiyTDAYYv+rtoF2PpgxNpm4LrIt3QwOaokSEbcoUWSZvQdORscd1FRNYXLpJJcjGp2S4gKm\nTDax9w+EaD/Z43JEIhOXEkWWiW+f8Ge1U4S6yIp4gxJFlnnFZyvapVKpRCHiCUoUWSQcDrN7n32N\nbH8nCvvkgBpLIeIeJYos0tLaHV3wp7gon9rqKS5HNDaqehLxBiWKLGKvdpo/t5ycnMSlyP1FiULE\nG5Qoskg2NWRDfBtFy4kudZEVcYkSRRaxL33q9/YJgNKSgmj33mB/KDp/lYhklhJFlujt6+fAoQ5r\nL8B8nw60S6TqJxH3KVFkib37TxEOm6qZ+toplBQXuBzR+LAnimPNmvNJxA1KFFnCvv5ENrRPRCS2\nU4hI5ilRZAl7j6dsShSqehJxnxJFFgiHw/E9nrKgITtCiULEfUoUWeBY8xm6uoOA6SlUXTl5mCv8\nI2509gnNIiviBiWKLJBY7RQI+HugnV1pqa2LbHBAXWRFXKBEkQVe2Zed1U4RlbO0iJGIm5QoskB8\nicKfS5+monYKEXcpUfhcd0+QQ0dOAxAIBJh33jSXIxp/mm5cxF1KFD5nShOmgXdO/VSKCvPdDcgB\n8cuiKlGIZJoShc+t23Akur3wguyrdgKortTobBE3KVH4WLA/FJcorry8zsVonGOvelIXWZHMU6Lw\nsc1bm+np7Qdg1oySrBqRbTe5dFJcF9kjx067HJHIxKJE4WNPP38oun3lyrqsGj+RaMnCmdHtzdua\nXYxEZOJRovCp3r5+Nmw6Ht2/emW9i9E4b+niiuj2lh0tLkYiMvE4nSiuBXYBe4A7hjjnu9bzW4CL\nbccPAluBTcALzoXoTxteOkYwOABAXc0UZteVuRyRs5YtiSWKHbvaCPaHXIxGZGJxMlHkAt/HJItF\nwE3AwoRz3gLMA+YDHwN+YHsuDDRikscKB+P0pWeePxzdvmplfVZXO4EZnV1hjdAOBgfiBhmKiLOc\nTBQrgL2YkkE/8CBwfcI51wH3W9vrgalAhe357L77jdKZrrNsstXTX5WlvZ0SLbNXP21X9ZNIpjiZ\nKGqAw7b9I9axkZ4TBp4CNgIfdShGX1q34Sih0CAA886bTlUWzRabir36Se0UIpmT5+Brj7Sz+1Cl\nhquAY8BM4ElMW8fTiSetWbMmut3Y2EhjY2M6MfrSM7beTtneiG134aJZBAIBwuEw+w6c4kzX2Wi3\nWREZWlNTE01NTaO+3slEcRSw14nUYUoMqc6ptY6BSRIArcDDmKqslIliIjh5qpftO1utvUDWDrJL\nprSkgHlzprNnv5m2ZOuOExPq3y8yWolfoteuXZvW9U5WPW3ENFI3AAXAjcAjCec8Arzf2l4JdAAt\nQDEQqU8pAa4BtjkYq288u/4wkcLa4gUzmD6tyN2AMsxe/bRV1U8iGeFkohgAVgOPAy8DDwE7gVus\nB8CjwH5Mo/c9wCes45WY0sNmTCP374EnHIzVNxJ7O000F10Y36AdDms6DxGnOVn1BPCY9bC7J2F/\ndZLr9gMXORKRjzW3dFnVLpCbm8MVK2pdjijzzp9bTmFhHn19A5xo66a5pWvCNOaLuEUjs33kmfWx\nRuyLllRMyIbcvLwcliyYFd1X7ycR5ylR+MjT6yZ2tVPE0iUaTyGSSUoUPvHq4Q4OH+0EoKAgjxWX\nVrsckXvsA++2vXwiOqZERJyhROET9tLE8ouqsnIlu5GqrZ4c7e3V09vPnv0nXY5IJLspUfhAOBzm\nWVv7xNWrJm61E5i1wZep+kkkY5QofGD3vpO0tHYDUFyUz8VLK12OyH3LNO24SMYoUfjAM+tipYmV\ny2spyM91MRpvsK9PsWffyehKfyIy/pQoPC4UGuS5F2Izn1y9SlNWAEwtK6Shfipg3qMd0WlNRGS8\nKVF43I5drZzq7AWgbEohSxbOGuaKiSO++knLo4o4RYnC4+xTdlx5eS25ufrIIuIbtE+4GIlIdtNd\nx8P6B0I8vzFW7XTl5RO7t1OiRRfMJD/PtNccPX6atvYelyMSyU5KFB62eWsLXd1BAGbNKGHB/HKX\nI/KWgoJcFl4wI7qvbrIizlCi8LCnbQsUXbmyLuvXxR4NrXon4jwlCo/as6/dWnvCmEgr2aXD3qC9\ndUcLg4OadlxkvClReFD/QIj/uG9jdK2FpYsrmF1X5nJU3tRQP5Upk80suqfPnOXgoQ6XIxLJPkoU\nHvTw71/h0JHYBID/+KFLVe00hJycQNzgO1U/iYw/JQqPOXSkk//+n5ej++991xIqZ5W6GJH3xY2n\nUIO2yLhTovCQUGiQu+/bGJ02e/555bzlTfNcjsr77A3au3a3EQyGXIxGJPsoUXjIY3/ax+59saVO\nP/mR5RpgNwIzyoupqZoCQLA/xM7dms5DZDzpLuQRLa3d/OJX26L77/zbhdTXqgF7pJYtsS2Pquon\nkXGlROEB4XCYH/70RYLBAQDqasp453ULXI7KX5Ytjk29rkQhMr6UKDyg6ZlX2bI9MqldgE/evDw6\nNYWMzOKFM8nJMT3DDhzqoKOzz+WIRLKHEoXLTnX08dP/2hzdf9vfzOf8eZqqI13FRflx79u2lzVJ\noMh4UaJw2Y8f2BQ3n9N73rnE5Yj8a+miWO+nX/9uJ319Ay5GI5I9lChctP7Fo6zbEJum4+M3L6ew\nMM/FiPzttVfUR6vsDh3p5O6fxEa3i8joKVG4pKs7yI9+9lJ0//VXz4kbOCbpq6qczMc+cEl0/5nn\nD/H7x/e4GJFIdlCicMkDD22Nrlw3rayID960zOWIssMbXjuHa143N7p//4Nb2L5T7RUiY6FEkWGD\ng2Ee/sMunmzaHz32kfdfTGlpgYtRZZcPv+8i5p9nGrYHB8N86z+e16JGImOgRJFBJzt6+fI3/soD\nD22NHlu5vJZVl9W6GFX2KcjP5fZbV1E2pRCAztN9fON76wj2a2oPkdFQosiQjZuPcdu/PMFW2+ym\n888r5x8/dKmLUWWvGeXFfOaTK6NjK/bsb+cnv9g8zFUikozf564Oe71XSzAY4oGHtvKHJ+2NqgHe\n8bYFvPsdi8nLU6520u/+uDtunMonb76MN7x2josRibjPWrZgxPd/9cV00OGjp/n2D56PW0xn2tQi\nPn3Lirg1FMQ5b/ub+ezZf5JnrGVlf3T/S8yuK2PeedNdjkzEP1SicEA4HObJpgP85D83R+dvAlh+\ncTWrP3JZdEU2yYy+vgE+/+U/RReDmlFezDfWvjHahiEy0aRbolCiGGcdnX3c+/OXWLfhSPRYfl4u\nH7hpKW9+4zytVOeS481nuP3Op+jp7QfMKO5/vf1qTeMuE5ISReYD4OjxM2zYdIwNm47xyp72uNHA\ndTVTuO3jK2mon+pilAKmQ8G//99novuvv3oO7/jbBVRXTnYxKpHM81qiuBb4f0Au8GPga0nO+S7w\nZqAH+CCwKY1rXUkUodAgO3e3sWHTMTZuOs7xljNJz7vmdXP50HuWMWmSmoK84sHf7OBXv90Rd6yh\nfiqrLjPdlGurp7gUmUjmeClR5AKvAG8EjgIbgJuAnbZz3gKstn5eDnwHWDnCayEDiWJgYJCTHb20\ntffQ0trN1u0tvLjleHQiv3MFWDC/nLe/dQGXXVKd8rWbmppobGwc95gzxY/xDw6G+fp3n+OFl47S\n2ryDmZWL456vqylj1WW1XLGilrqaKZ6uKvTj+x/h59jB//F7qdfTCmAvcNDafxC4nvib/XXA/db2\nemAqUAnMGcG1aRscDBMMhjgbHKDvbIizZwfMIxiiu6eftvYe2k/20NreQ1t7D20nezl5qnfYieUm\nTcrj4gsrWX5RFZcsq2Jq2cgaSf3+x+bH+HNyAtx+6yqeXX+YtWv/SH7eUvoHYgPxDh/t5PDRTn71\n2x3UVE1h+UVVlE8vorS0gNIS85hs23azjcOP73+En2MH/8efLicTRQ1w2LZ/BFNqGO6cGqB6BNcC\n8Pm1f2IwHGYwFCY0GGbQeoQGBxkMhRkMhzl7NsTZYCiuB9JYTZ9WxGUXV3PZxdUsWTiLggItNOQX\nubk5vOaK2Vy9qp7P3XEdL205zroNR3hxS3Pc38jR46c5evx0ytcqLsqnpLiAgoJccnIC5OYGyAkE\nyLH9zM0JkJMTIBAIREsokYJK9Kfty91ICzFPrzsU1+biJ36OHfwff7qcTBQjrRMaU9l+9772sVw+\nItPKiphRXsyM8iLqasq47JJqzps91dPVEjIyxUX5XLWynqtW1tPXN8BLW03S2Lj5OGfPDv/Foqe3\nP9qTKtOONZ9h4+ZjrvzusfJz7OD/+L1kJfBH2/4/A3cknPND4N22/V1AxQivBVM9FdZDDz300COt\nx148Ig/YBzQABcBmYGHCOW8BHrW2VwLPp3GtiIhkgTdjei/txZQKAG6xHhHft57fAlwyzLUiIiIi\nIiLj41pMm8YekrdfeN1BYCtmgOEL7oYyrJ8ALcA227HpwJPAbuAJTNdmr0oW/xpMb7pN1uPazIc1\nYnXAn4EdwHbgU9Zxv3wGQ8W/Bn98BoWY7vubgZeBu6zjfnj/h4p9Df5478ckF1Ml1QDk4882jAOY\nPzQ/uBq4mPgb7deBz1nbdwBfzXRQaUgW/53AP7kTTtoqgYus7VJMlexC/PMZDBW/nz6DYutnHqYt\n9Sr88/4niz2t996vM6LZB/P1ExuQ5zd+6V/7NHAq4Zh9sOT9wA0ZjSg9yeIH/7z/zZgvQwBdmIGn\nNfjnMxgqfvDPZxBZS7cA80X1FP55/5PFDmm8935NFEMN1POTMPAUsBH4qMuxjEYFpjoH66cfF9i4\nFdOJ4j68WW2QTAOmdLQef34GDZj4Iz0c/fIZ5GCSXQuxajS/vP/JYgf/vPej9k7gXtv++4DvuRTL\naFVZP2diPsSrXYxlJBqIr7pJ/IZ+MnOhjEoD8fHPwnyjCgD/hvnP4nWlwIvEvrn67TMoxXwxisTv\nx8+gDJOZUNEBAAAGDElEQVTkXof/3v9I7I2k+d77tURxFNNAFlGHKVX4yXHrZyvwMKY6zU9aMHXP\nYJLeCRdjGY0TxAYf/Rjvv//5wK+BB4DfWsf89BlE4v8Fsfj99hkAdAJ/AC7FX+8/xGJfTprvvV8T\nxUZgPrEBeTcCj7gZUJqKgcgiCCXANcR/2/WDR4APWNsfIPaf3y+qbNtvx9vvfwDzje9lzNT7EX75\nDIaK3y+fwQxiVTNFwJswPYX88P4PFXul7Rwvv/dj5ucBeXMw1U2bMd0FvR7/L4FjQBDTNvQhTI+t\np/B218CIxPg/DPwc0z15C+Y/uFfrl8H0UhnE/L3YuzP65TNIFv+b8c9ncCHwEib+rcDt1nE/vP9D\nxe6X915ERERERERERERERERERERERERERETET+7CTBdwA/D5Ub7GGuAz4xSPUxqB3w1zThnwcedD\nGRe3AP9gbX+Q+IFx9+K/mZpljPw6Mlv8YQVmbpnXAn8d5WuExy+cIWXi/8E04BMZ+D3j4R7MVCFg\nRhxX2577KGb2VxGRMfk6ZsTnacwo3NPW/hcTzssF9lvbU4EQZhQvmMQyDzNv/n2YWS/3YWa8jPgn\nzNQD24BPDxHL3cAGzAj4NbbjBzHrB7yImQLmGuA5a/9XmKlVEl1GbLGpbxCb9qCRWIliDfEloG3A\nbMxU+D3WtV+znrsds2jVFltsDZgb8Y+smB/HLD6T6O+s194M/MU6lmvFFXnNj9niawL+23rtX9he\n56uY2US3YD43+7/hncAZzAJhL1lxNGHmObrFdj6YkkdkYs73YWa33QT8EH0hFZEhLAe+g1ks5ZkU\n5z0GLALehrnB/QswiVgCWQM8i5lUrhxow9wQL8XctIswN/XtxBbHsZtm/czFJJsl1v4B4LPW9gzM\nzbbI2r8D+Nckr7UduNzavsv6/RCfKO7k3ERRj0kW9vl0rsF8cwdzI/0dZgbhBswaK0ut5x4C3psk\nlq3EqoSmWD8/BnzB2p6ESZANVnwdmJJBAJMQr8S8n7tsrxl5HfuiNn8mfi37yP4MzOqSEY8CV2Cq\npR7BvN9gEvU/IL6mTC9OidzIF5K6quJp4DWYm+RdmBLFcsxNDkzV0+8xN892zKyXldZ5vwF6gW5r\nO9lU7TdiSgkvAYsxSSniIevnSuv4c5hvwe/H3NztpmKmyV5v7f8X6S26k3juNdZjkxXfBZgSFJgk\nFklCL2Ju9omexSyW8xFMMo685vut13weMxfRPMx7+AJmvqswphQyG5M8+jAltrdj3suRxA4mYe/H\nJM5yYAHm/XsD5rPfaMXxeszcZuJjecOfIpKWZcDPgFrMzaQYc6N5CfONsy/h/L9i6u6rgC9hqmMa\niW/TCNq2Q5i/2zDxN7AA57ZnzMF8u1+OmWL5p8RX43Tbtp8E3jPsvy7+9yUzQPwXsGTVRhF3YaqY\n7BqAs7b9ELGSjt3HMW1Ab8Ukk0ut46sx/xa7xiSvmW/9XIG5ub/LuvYNSX7XUO1EDwJ/jymV/MZ2\n/H5MyVCyhEoUMt62YFYw240pTfwv5pvuJZybJMB8070Cc9M6a11/C6kbv8OYksgNxKqebrCO2U3B\nJIPTmNkx3zzE663HVMXMtfZLMNPY23Vg6usj8/a/e4jXOkisquYSYt+mzxCbWh5M28OHibWF1GAW\nsRqpuZj37k7MmiZ11mt+gtgXwPOJrZecTAmmpPQYpqppmXU8sqBNJO4p514KmHVUbgBuwiQNgD9h\nkk7k3zKdc0tn4jMqUYgTZhJb7WsB8fXgiYLAIWJLY/4VU11kr89P9o12E6bk8oK1fy8mydhtsc7b\nhZlefKi2klZMY+wvMXX7YOr69yScd7P1ewYxbRqdSWL8Nab6ZzsmAb1iHW/HVBdtw9Tn34FJpOus\n589gGoHDnPvvTfbv/zommQUwU11vwVRXNWBKbwFMNd3bU7zmZOB/MKWeAHCb7bnI+T/DNEj3YBK6\nXQdmjYmFmKomMNWMX8RMu52DqTL8BOYzFhHJevaeUJ8Hvu1WICIi4k1/jymhbMP0Uip3NxwRERER\nERERERERERERERERERERERERSdP/B8j5v3/66oMAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc475574b90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "thinkplot.PrePlot(1)\n",
    "thinkplot.Pdf(posterior)\n",
    "thinkplot.Config(xlabel='# who are gluten sensitive', \n",
    "                 ylabel='PMF', legend=False)\n",
    "#thinkplot.Save(root='gluten1', formats=['png'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the posterior credible intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 13)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior.CredibleInterval(90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 13)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior.CredibleInterval(95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare to a null hypothesis in which none of the respondents are gluten sensitive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.3956946672325594e-05"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior2 = Gluten([0])\n",
    "prior2.Update(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bayes factor in favor of the null hypothesis is about 8, which is moderate evidence.  If you started thinking the probability of the null was 50%, you should now believe it is about 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.461538461538462"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.11 / 0.013"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8947368421052632"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "8.5/9.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}