ch.3/Bayesian Spam Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Spam Detection\n",
    "\n",
    "Here, we'll take a look at spam detection using the bag-of-words model and the application of Bayesian detection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two hypotheses are H_0 (normal email) and H_1 (spam/phishing email).\n",
    "\n",
    "The first step is to choose a dictionary. I've chosen this fairly arbitrarily from examples I found online. In a real-world system, you would need to estimate the probabilities for each word given the two hypotheses, but here I have made them up entirely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "dictionary = ['order','report','mail','colleague','language','one','report','send','money','free']\n",
    "p_0 = [10,5,15,10,5,20,20,10,5,1]\n",
    "p_0 = p_0/np.sum(p_0)\n",
    "\n",
    "p_1 = [7,5,20,15,10,15,13,20,15,10]\n",
    "p_1 = p_1/np.sum(p_1)\n",
    "\n",
    "#Note that p_0 and p_1 are NOT the likelihood functions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate a \"message\" according to each hypothesis. Obviously this is not going to generate a realistic-looking message owing to the limited dictionary and the bag-of-words model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['one' 'report' 'one' 'one' 'order' 'money' 'one' 'colleague' 'report'\n",
      " 'send' 'mail' 'report' 'report' 'mail' 'order' 'report' 'order' 'one'\n",
      " 'report' 'report' 'order' 'money' 'report' 'mail' 'send' 'send' 'order'\n",
      " 'one' 'money' 'one' 'report' 'send' 'one' 'colleague' 'report' 'one'\n",
      " 'money' 'colleague' 'one' 'order' 'one' 'mail' 'one' 'report' 'colleague'\n",
      " 'mail' 'one' 'report' 'send' 'report']\n",
      "['one' 'one' 'send' 'colleague' 'send' 'colleague' 'language' 'order'\n",
      " 'order' 'report' 'colleague' 'report' 'one' 'report' 'money' 'one' 'mail'\n",
      " 'one' 'one' 'report' 'language' 'send' 'send' 'send' 'report' 'send'\n",
      " 'language' 'one' 'colleague' 'money' 'colleague' 'colleague' 'colleague'\n",
      " 'one' 'send' 'colleague' 'report' 'mail' 'one' 'report' 'language'\n",
      " 'colleague' 'report' 'report' 'order' 'send' 'free' 'send' 'one'\n",
      " 'language']\n"
     ]
    }
   ],
   "source": [
    "N = 50\n",
    "\n",
    "normal_message = np.random.choice(dictionary,size=N,p=p_0)\n",
    "spam_message = np.random.choice(dictionary,size=N,p=p_1)\n",
    "\n",
    "print(normal_message)\n",
    "print(spam_message)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's build a Bayes classifier. To do so, we need to choose priors and cost functions, then we compute the LLR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pi_0 = 0.9\n",
    "pi_1 = 1-pi_0\n",
    "\n",
    "C_00 = 1\n",
    "C_11 = 1\n",
    "C_10 = 2\n",
    "C_01 = 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's sample directly from the spam bag of words model using np.multinomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 3 1 1 0 1 3 1 0]\n"
     ]
    }
   ],
   "source": [
    "y = np.random.multinomial(N,p_1)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we compute the LLR and classify."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.18479008414\n",
      "0.587786664902\n"
     ]
    }
   ],
   "source": [
    "L = 0\n",
    "for i in range(0,len(dictionary)):\n",
    "    L += y[i]*np.log(p_1[i]/p_0[i])\n",
    "\n",
    "threshold = np.log(((C_10-C_00)*pi_0)/((C_01-C_11)*pi_1))\n",
    "print(L)\n",
    "print(threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did we correctly detect the spam? Let's get a better sense of the detector performance by looking at the Bayes risk, averaged over many samples and for many priors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fy4311\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "C:\\Users\\fy4311\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8VFX6x/HPkw4hgSRDh1ASeocQ\nIFFXd1fFil2xK4i6um511V1/uqvu6uo2Xdcui6JgRUXFunZ66L2EEiIthBASQuo8vz9mcLNIkgEy\nc6c879crL2buvZP55pLkyTn33HNEVTHGGGOCTZTTAYwxxpgjsQJljDEmKFmBMsYYE5SsQBljjAlK\nVqCMMcYEJStQxhhjgpIVKGOMMUHJCpQxxpigZAXKGGNMUIpxOkBzcrlc2r17d6djGGOMacSiRYv2\nqGrbpo4LqwLVvXt38vLynI5hjDGmESKy1ZfjrIvPGGNMULICZYwxJihZgTLGGBOUrEAZY4wJSlag\njDHGBCUrUMYYY4KSFShjjDFByQqUMSbgSg5Us3DLXtxudTqKCWJhdaOuMSZ4qSqLtpbw8vwC3l+x\ng+paN6f1b8/fLh1Kq3j7VWS+z74rjDF+tb+yhreXfMvL8wpYt6uMVvExXJrVlXZJ8fzjPxu44InZ\nPHt1Ft3SEp2OaoKMFShjjF+sKCzl5flbeWfpdg7W1DGwczIPXTCIc4Z0ItHbYhreLYVbpi3m3Mdn\n8/jlwzixV5PTs5kIIqr+6QMWkcnA2cBuVR14hP3jgPsBN1AL/FxVv/HuqwNWeA8tUNVzfXnPrKws\ntbn4jHFORXUtM5duZ9qCApYXltIiNppzh3TiitHpDO7S5oivKSiu4IYX89iwu4zfntmPCSf0QEQC\nnNwEkogsUtWsJo/zY4E6CSgHXmygQLUCDqiqishg4DVV7evdV66qrY72Pa1AGeOMtTv3M21+AW8t\n/payqlp6t2/FlaO7cd6wziQnxDb5+gNVtfzqtWV8uGonFwzrzJ8uGERCbHQAkhsn+Fqg/NbFp6pf\niUj3RvaX13uaCNhwHmNCSGVNHR+s3MHL8wrI21pCXEwUZw3qyBWj0hnRLeWoWkGJ8TE8ccVwHv98\nI3/7ZD35ReU8fVUWHVon+PErMMHO0WtQInI+8CDQDjir3q4EEcnD0/X3kKq+7UQ+Y8z3bSoqZ/qC\nAl5fVMi+ihp6uBL53Zn9uHBEF1IT447580ZFCbf9qBd9OiTxy1eXcs7j3/DUlSMY0S2lGdObUOK3\nLj4AbwvqvSN18R123EnAPar6Y+/zTqq6XUR6Ap8BP1LV/AZeOwmYBJCenj5i61aflhkxxhyF6lo3\nn6zexbQFW5m9sZiYKOH0AR24fFQ6Y3qmERXVvNeM1u8q44YX89ixr5L7zxvApSPTm/XzG2c5fg3K\nG6I7PhQo77GbgZGquuew7VO8n+ONpj6HXYMypnlt21vBKwsLeHVhIXvKq+jcpgWXj0rn4qwutEvy\nb/fbvopqfjp9CV9v2MM1Y7px99n9iY22uQXCgePXoJoiIplAvneQxHAgDigWkRSgQlWrRMQF5AIP\nO5XTmEhT51Y+X7ubl+dv5Yv1RQjww77tuGJUN07q3ZboZm4tNaRNyzj+fe1I/vzhWp79ejPrdpXx\nxBUjjqsb0YQWvxUoEZkOnAy4RKQQuBeIBVDVp4ALgatFpAY4CFzqLVb9gKdFxI1nKqaHVHW1v3Ia\nYzx27a/k1YXbeGVBAdtLK2mXFM9PT8nk0ux0Ordp4UimmOgofndWf/p1TObOGSs455/f8OzVWfTv\nlOxIHhNYfu3iCzTr4jPm6Ljdyjcb9zBtfgGfrNlFnVs5sZeLK0al86N+7YOqS23Ztn3cOHURpQdr\n+MvFQzhrcEenI5ljFPRdfMYY5xSXV/H6okKmzS+gYG8FqYlxTDyxB+NHptPdFZxTDg3p2oaZP83l\npqmLuGXaYtbsyOSXp/Zu9gEaJnhYgTImQqgqCzbv5eX5BXy4cifVdW6ye6Tyq9N6M3ZgB+Jjgv/G\n2HZJCUyfNJp73l7F459vZO3O/fz90qEk+XAzsAk9VqCMCXOlFTXMWFLIy/ML2Li7nOSEGK4Ync7l\n2en0ap/kdLyjFh8TzUMXDmJA52T+8O5qzn9iDs9cNYKebY968hkT5OwalDFh6NDSFtO8S1tU1boZ\n2rUNV4xK5+zBnWgRF/ytJV/MzS/mlmmLqalz88/xwzi5TzunIxkfBMV9UIFmBcpEun0V1cxY/C3T\nFxSwYXc5SfExnDesM5dld2VAp9ZOx/OLbXsrmDR1Eet27ueOsX2ZdFJPm2w2yNkgCWMihKqSt7WE\n6Ye1lh6+cDBnD+lIy7jw/jHvmtqSN28ew+2vL+fBD9ayesd+/nzhYJtsNgyE93euMWFsX0U1b3pb\nSxu9raVLsroyPjs94u4TahkXw+OXD6P/F8n85eN1300269T9W6Z5WBefMSFEVVm4pYTpC/67bPqw\n9DaMz07n7MHh31ryxX/W7OJnrywlITaKJ68cwcjuqU5HMoexa1DGhJGSA9W8ubiQ6QsKyC86QFJ8\nDOcP78xlIyOvteSLjbvLuOHFRRSWVPCHcwdy+SibbDaY2DUoY0LcofuWpi8oYNbKnd+1lh65aDBn\nWWupUZntknj7llxum76E3761glXbS7n3nAHExQTPzBimafYdbkyQOdRamraggE1FB0hKiGH8yK5c\nlp1Ov47WWvJV6xaxTL52JA9/tJanv9zEhl3lPHHlcFyt4p2OZnxkXXzGBAFVZb63tfTBCs8sD8PT\n23D5qG6cNahj2Ny35JR3ln7Lb95YTlpiHM9cncXAzuE55D5UWBefMSFg74FqZhzWWrp8VDqXZXel\nbwdrLTWXcUM709PVihun5nHRU3N4+KIhnDukk9OxTBOsQBkTYKrKvE2e1tKhOfFGdEvhLxdnWmvJ\njwZ1ac3Mn57AzS8t4rbpS1i9fT+3n94nYOtbmaNnBcqYANl7oJo3F3lG4m3ac4Bkb2tpfHY6fTqE\n3px4ocjVKp6XJ47mD++u4qkv81m7cz+PXjaM1i1sstlgZAXKGD9SVeZuKmb6gm185G0tZXVL4ZZT\nMjnTWkuOiIuJ4o/nD6Jfx2R+P3MV5/9rNs9cnUVmO5tsNthYgTLGD4rLq7z3LW1jc73W0uWj0ukd\ngjOIh6MrR3ejd/skbn5pEef/azaPjh/KD/u2dzqWqcdG8RnTjCqqa/nj+2t4Pa+Q6jo3I7unMD47\nnTMHdbS54YLUt/sOcuPUPFZt38+dY/ty4w8ynI4U9mwUnzEBtm5nGbdMW0x+UTlXjErn6jHdrbUU\nAjq3acHrN+Zw+xvLePCDtUSJcMNJPZ2OZbACZcxxU1VeXbiNe2euIikhlpcmjCI30+V0LHMUWsRF\n8+hlw1CFP85aQ1JCDJdl2/RITrMCZcxxKKus4bdvreTdZds5IdPF3y8dStskm6kgFEVHCX+/dCjl\nVbXc9dYKWiXEcPZgu1fKSX6dmEpEJovIbhFZ2cD+cSKyXESWikieiJxQb981IrLB+3GNP3MacyxW\nflvK2f/8hveXb+f20/vw4vXZVpxCXFxMFE9dOYKR3VL5+StL+XztbqcjRTR/z5w4BRjbyP7/AENU\ndShwPfAcgIikAvcCo4Bs4F4RSfFvVGN8o6pMmb2ZC56YQ1WNm1cmjeGWUzKJshs+w0KLuGieuzaL\nvh2TuOmlRczfVOx0pIjl1wKlql8BexvZX67/HUaYCBx6fDrwiaruVdUS4BMaL3TGBERpRQ03Tl3E\n799dzQm9XMz62Ylk97D1hsJNckIsL1yXTZeUFkx4IY8VhaVOR4pIjs89LyLni8ha4H08rSiAzsC2\neocVercd6fWTvN2DeUVFRf4NayLa4oISznzsaz5bu5u7z+rH89dkkZoY53Qs4ydpreJ5aeIoWreI\n5erJ89mwq8zpSBHH8QKlqm+pal/gPOB+7+Yj9ZUc8YYtVX1GVbNUNatt27b+imkimNutPP1lPpc8\nNZeoKHjj5hwmntgTEevSC3cdW7fg5YmjiImO4srn57Ntb4XTkSKK4wXqEG93YIaIuPC0mLrW290F\n2O5IMBPRisuruP6FhTz4wVpOG9Ce9356IkO7tnE6lgmg7q5Epk7IprLGzZXPz2f3/kqnI0UMRwuU\niGSK989QERkOxAHFwEfAaSKS4h0ccZp3mzEBM29TMWc+9jVz8ou5/7yB/Ovy4TapaITq2yGZKdeN\npKisiiufn0/JgWqnI0UEfw8znw7MBfqISKGITBCRm0TkJu8hFwIrRWQp8C/gUvXYi6e7b6H34z7v\nNmP8rs6tPPafDVz+7DwS42J46yc5XDW6m3XpRbhh6Sk8d00WW4oruHbKQsqrap2OFPZsLj5j6tm9\nv5Kfv7qUOfnFnDe0Ew+cP4hW8XY/u/mvT1bv4qaXFpHdPZV/XzfS5lg8Br7OxRc016CMcdpX64s4\n87GvWVxQwsMXDebvlw614mS+59T+7fnbJUOYt7mYW6ctpqbO7XSksGUFykS82jo3D3+4lqsnLyA1\nMY53bz2BS7K6WpeeadC4oZ25b9xAPl2zm1+/vgy3O3x6ooKJ/XloItr2fQe5bfoS8raWcNnIrtx7\nzgBbRND45KrR3SirrOHhD9fRKj6GB84baH/UNDMrUCZifbp6F79+Yxk1tW4evWwo44Ye8V5wYxr0\nk5MzKaus5ckv8kluEcsdY/s6HSmsWIEyEae61s2fP1zL899sZkCnZB6/fDg9XIlOxzIh6jen92H/\nwRqe/CKfpIQYfnJyptORwoYVKBNRCooruHX6YpYXlnJtTnfuOrMv8THWpWeOnYhw/7iBlFfV8vCH\n60hKiOWq0d2cjhUWrECZiPH+8h3c+eZyROCpK4czdmBHpyOZMBEVJfzl4iEcqKrlnndWkhQfw3nD\nrMv4eNkoPhP2KmvquPvtFdwybTEZ7Vrx/m0nWnEyzS42OorHLx/O6B5p/Or1ZXyyepfTkUKeFSgT\n1vKLyjnvX7N5aV4BN57Uk9dvGkPX1JZOxzJhKiE2mmevyWJg59bcMm0xczbucTpSSLMCZcLWjMWF\nnPPPb9hdVsW/rxvJXWf2IzbavuWNf7WKj+GF60bSIy2RiS/msaSgxOlIIct+Wk3Yqaiu5devL+OX\nry1jYOfWzLrtRE7p087pWCaCtGkZx9QJ2bhaxXPtvxeydud+pyOFJCtQJqys3bmfc/75DW8uLuS2\nH/Vi2sRRdGid4HQsE4HaJSfw8sRRtIiN5qrnF7BlzwGnI4UcK1AmLKgq0xcUMO7x2eyvrOWlCaP4\n5am9ibEuPeOgrqkteWliNrV1bq54bj47Sg86HSmk2E+vCXlllTXc9spS7pqxguweqcy67URyM11O\nxzIGgMx2Sbx4/Sj2H6zhyufmU1xe5XSkkGEFyoS0dTvLOPuf3zBrxQ5uP70PL1yXTdukeKdjGfM/\nBnVpzfPXjqSw5CDX/HsB+ytrnI4UEqxAmZC1Zc8BrnhuPger63hl0mhuOSWTqCibrNMEp+weqTx1\n1QjW7ihj4pQ8DlbXOR0p6FmBMiFpZ2klVz4/nzq3m2k3jGJk91SnIxnTpFP6tOMflw0lb+tebnpp\nEdW1tpZUY6xAmZCzr6KaqyfPp+RANS9cn01muySnIxnjs7MHd+LBCwbx5foifvHqUupsLakG2Vx8\nJqQcqKrl2n8vZEtxBVOuG8ngLm2cjmTMUbt0ZDpllbU88P4aWsXH8NCFg2wtqSOwAmVCRlVtHTdO\nXcTywn08eeUIcjJspJ4JXRNP7Mn+gzU89tlGWiXEcPdZ/axIHcZvBUpEJgNnA7tVdeAR9l8B3OF9\nWg7crKrLvPu2AGVAHVCrqln+ymlCQ51b+fkrS/lm4x4euWgwpw/o4HQkY47bL07tzf7KWp7/ZjOt\nW8Ry2496OR0pqPizBTUFeBx4sYH9m4EfqGqJiJwBPAOMqrf/FFW1mRYNqspvZ6zgg5U7ufusflyc\n1dXpSMY0CxHhnrP7U1ZZy98+WU+r+BiuP6GH07GCht8KlKp+JSLdG9k/p97TeUAXf2Uxoe2hD9fy\nat42fvrDTCae2NPpOMY0q6go4c8XDuJAVS33vbeapIQY+yPMK1hG8U0APqj3XIGPRWSRiExq7IUi\nMklE8kQkr6ioyK8hTeA9+UU+T3+5iatGd+OXp/Z2Oo4xfhETHcWj44dyYi8Xd7y5nA9X7nA6UlBw\nvECJyCl4CtQd9Tbnqupw4AzgFhE5qaHXq+ozqpqlqllt27b1c1oTSNMXFPDnD9dy7pBO/OHcAXYB\n2YS1+Jhonr5qBMPSU/jp9CV8td7+4Ha0QInIYOA5YJyqFh/arqrbvf/uBt4Csp1JaJzy/vId/Pat\nFZzSpy1/vWSIzRBhIkLLuBgmXzuSzHZJTJqaR96WvU5HcpRjBUpE0oEZwFWqur7e9kQRSTr0GDgN\nWOlMSuOEr9YX8fNXl5DVLYUnrhhhiwyaiNK6RSxTJ2TTqXULrpuykDU7InctKb/95IvIdGAu0EdE\nCkVkgojcJCI3eQ+5B0gDnhCRpSKS593eHvhGRJYBC4D3VfVDf+U0wWXR1hJunLqIzHZJPHfNSFrE\nRTsdyZiAc7WKZ+rEUcRFR/GPT9c3/YIw5c9RfOOb2D8RmHiE7ZuAIf7KZYLX2p37uX7KQtonx/Pi\n9dm0bhHrdCRjHNO5TQt+3K89H6zcQZ1biY7Abm7rOzFBoaC4gqueX0BCbBRTJ4yyJTOMAXIy09hf\nWcvq7ZHZzWcFyjhu937PzOQ1dW5emjCKrqktnY5kTFAYk5EGwOz8yJyzwAqUcVRpRQ1XPb+APeVV\nTLkum17tbWZyYw5pl5RAr3atmL3RCpQxAVVRXct1Uxawec8Bnr06i6FdbWZyYw6Xm+li4Za9Ebl2\n1FEXKBHpIyLP+iOMiRzVtW5unLqIpdv28dj4oeRm2szkxhxJTkYalTVulhSUOB0l4BosUCIyWEQ+\nFpGVIvKAiLQXkTeB/wCrAxfRhJs6t/KL15by9YY9PHTBYMYO7Oh0JGOC1qieaUQJzMkvbvrgMNNY\nC+pZYBpwIVAELAY2AZmq+vcAZDNhSFW5++2VvL98B787sx+XjLRJMY1pTOsWsQzq3Jo5EThQorEC\nFa+qU1R1nao+CriBO1W1MkDZTBh65KN1TF9QwE9OzuCGk2xmcmN8kZPpYknBPg5U1TodJaAaK1AJ\nIjJMRIaLyHA8iwoOrvfcmKPyzFf5PPFFPpePSuf20/s4HceYkJGTkUatW1kYYXPzNTaTxA7gb/We\n76z3XIEf+iuUCT+vLizgT7PWcvbgjtw/bqDNTG7MUcjqlkpcdBRz8os5uU87p+METIMFSlVPCWQQ\nE74+XLmDu2as4Ae92/K3S4ZG5JQtxhyPFnHRDEtvE3H3Q9l9UMavvtmwh9umL2VYegpPXjmcuBj7\nljPmWORmuli9Yz8lB6qdjhIw9tvC+M2SghImTc2jZ9tEJl8zkpZxfpub2Jiwl5uZhirM2xQ5w82t\nQBm/WL+rjOumLKRtkndm8pY2M7kxx2NwlzYkxkVH1Lx8TRYoEbnvsOfRIvKy/yKZULdtbwVXPT+f\nuOgoXpowinbJCU5HMibkxUZHkd0jNaJu2PWlBZUuIncBiEg8niXYN/g1lQlZu8s8M5NX1riZajOT\nG9OscjNdbCo6wM7SyLgd1ZcCdR0wyFuk3gU+V9Xf+zWVCUmlB2u4+vkFFJVV8e/rRtKng81Mbkxz\nOrT8RqTMKtHYXHyHbsgdBjwKXIqn5fSl3ahrDnewuo4JUxaSX1TO01eNYHh6itORjAk7/Tokk9Iy\nltkbI6Obr7FhVX897HkJ0N+73W7UNd+prnVz88uLWFxQwuOXD+fEXm2djmRMWIqKEsZkpDEnfw+q\nGvY3vNuNuua41LmVX72+jC/WFfHQBYM4c5DNTG6MP+VkuJi1Yidbiivo4Up0Oo5f+TKK72cikiwe\nz4nIYhE5zYfXTRaR3SKysoH9V4jIcu/HHBEZUm/fWBFZJyIbReTOo/uSTKCoKvfOXMm7y7Zz5xl9\nuSw73elIxoS9Q2unRcKsEr4MkrheVfcDpwHt8AyaeMiH100BxjayfzPwA1UdDNwPPAOeYezAv4Az\n8HQpjheR/j68nwmwv368npfmFXDTDzK46QcZTscxJiJ0T2tJx9YJETFQwpcCdaiT80zg36q6rN62\nBqnqV0CDU++q6hxVPbRE5Dygi/dxNrBRVTepajXwCjDOh5wmgJ77ehOPf76R8dlduWOszUxuTKCI\nCDkZLubmF+N2q9Nx/MqXArVIRD7GU6A+EpEkPGtDNacJwAfex52BbfX2FXq3HZGITBKRPBHJKyoq\nauZY5khey9vGA++v4axBHXngvEFhf6HWmGCTm5lGSUUNa3budzqKX/lSoCYAdwIjVbUCiMPTzdcs\nROQU73vccWjTEQ5r8M8EVX1GVbNUNattWxs95m8frdrJnW8u58ReLv526RCbmdwYB+RkeK5DzQ3z\nWSUauw+qr/fhUO+/Pb33P3Wj8eHpPhORwcBzwDhVPXSmC4H664B3AbY3x/uZ47Py21J+Om0JQ7q2\n4emrRhAfE+10JGMiUofWCfRsmxj2AyUaKzS/BCbx/fuhoBnugxKRdGAGcJWqrq+3ayHQS0R6AN8C\nlwGXH897meNXXevm168vo03LWJuZ3JggkJvhYsbiQmrq3MRGh+e8343dBzXJ++8x3Q8lItOBkwGX\niBQC9wKx3s/5FHAPkAY84b2GUevtqqsVkVuBj4BoYLKqrjqWDKb5PPHFRtbuLOO5q7NISYxzOo4x\nES8nI42p87ayvHAfI7qlOh3HL476z2ARORX4jaqe2thxqjq+if0TgYkN7JsFzDrabMY/Vm/fz+Of\nbeT8YZ35cf/2TscxxgCje6YhArM3FodtgWrsGtQPRWS9iJSLyEsi0l9E8vDcA/Vk4CIaJ9XUubn9\njWW0aRnHvefY7WjGBIuUxDj6d0wO6+tQjXVc/hXPNag04A089ypNVdURqjojEOGM857+Mp9V2/fz\nwHkDadPSuvaMCSa5mS6WFOzjYHWd01H8orECpar6hapWqerbQJGqPhqoYMZ563aW8eh/NnDOkE6M\nHdjB6TjGmMPkZKRRXecmb2uDcyKEtMauQbURkQvqPZf6z60VFd5qvV17yQmx/OHcAU7HMcYcwcju\nqcRECbM3FoflKgKNFagvgXMaeK54hoibMPXs15tZXljKvy4fTqqN2jMmKCXGxzAsvQ1zw3RevsaG\nmTfbbBEmtGzcXcbfP13PGQM7cNZgWz7DmGCWk+Hin59toPRgDa1bxDodp1mF591d5pjVuZXb31hO\nYlw0940b6HQcY0wTcjLScCvM3xR+0x5ZgTL/Y/I3m1lSsI/fnzuAtknxTscxxjRhWHoKCbFRzAnD\nefmsQJnvbCoq5y8fr+PU/u05d0gnp+MYY3wQFxPFyO6pYXk/lC8r6l7sXWIDEblbRGZ4J401YaTO\nrfzmjeUkxEbzx/MG2hIaxoSQ3EwXG3aXs7us0ukozcqXFtT/qWqZiJwAnA68gM0kEXZemLOFvK0l\n3HtOf9olJzgdxxhzFHLDdPkNXwrUoVuUzwKeVNV38KwJZcLElj0HePijtfywbzvOH9bg2pDGmCDV\nv1MyyQkxYdfN50uB+lZEngYuAWaJSLyPrzMhwO1WfvPmcmKjo/jT+bY6rjGhKDpKGJORFnYDJXwp\nNJfgWfpirKruA1KB2/2aygTMS/O3smDzXv7v7P50aG1de8aEqtxMF4UlBykornA6SrNpskB5l3nf\nDZzg3VQLbPBnKBMY2/ZW8NAHazmpd1suHtHF6TjGmOOQk5EGwOwwmlXCl1F89wJ3AHd5N8UCL/kz\nlPE/VeWON5cTJcJDF1jXnjGhLqNtK9olxYdVN58vXXznA+cCBwBUdTuQ5M9Qxv+mLShgTn4xvzur\nH53atHA6jjHmOIkIORlpzM3fg6o6HadZ+FKgqtXz1SqAiCT6N5Lxt2/3HeTBWWs5IdPFZSO7Oh3H\nGNNMcjJd7CmvZv2ucqejNAtfCtRr3lF8bUTkBuBT4Fn/xjL+oqrc+eZy3Ko8aF17xoSV765Dhclw\nc18GSfwFz4q6bwJ9gHtU9Z/+Dmb847W8bXy9YQ93ndGXrqktnY5jjGlGXVJa0i2tJXPCZKCEL4Mk\nbgXyVPV2Vf21qn7iyycWkckisltEVjawv6+IzBWRKhH59WH7tojIChFZKiJ5Pn0lpkk7Sg/ywHtr\nGN0zlStGdXM6jjHGD3IyXMzftJfaOrfTUY6bL118HYCFIvKaiIwV3/uEpgBjG9m/F7gN+EsD+09R\n1aGqmuXj+5lGqCq/nbGCWrfy8IVDiIqyrj1jwlFuZhplVbWs+LbU6SjHzZcuvruBXsDzwLXABhH5\nk4hkNPG6r/AUoYb271bVhUDNUSU2x+TNxd/y+boifjO2D+lp1rVnTLga09NzHSochpv7NGWRdxTf\nTu9HLZACvCEiD/splwIfi8giEZnU2IEiMklE8kQkr6ioyE9xQtuu/ZXc9+4qsruncs2Y7k7HMcb4\nUVqrePp2SAqL61C+XIO6TUQWAQ8Ds4FBqnozMAK40E+5clV1OHAGcIuInNTQgar6jKpmqWpW27Zt\n/RQndKkqv3trBVW1bv580WDr2jMmAuRkuMjbUkJlTV3TBwcxX1pQLuACVT1dVV9X1RoAVXUDZ/sj\nlPdmYFR1N/AWkO2P94kE7yzdzqdrdnP76X3o4bJb2IyJBLmZaVTVullcUOJ0lOPiyzWoe1R1q4i0\nE5H0Qx/efWuaO5CIJNZbIDEROA044khA07jdZZX8/t1VDE9vw3W5PZyOY4wJkOweqURHCXM2hvZ1\nqJimDhCRc4C/AZ3wTBrbDVgDDGjiddOBkwGXiBQC9+KZxw9VfUpEOgB5QDLgFpGfA/3xtNje8g4W\njAGmqeqHx/LFRTJV5f/eXklFdR0PXzSEaOvaMyZiJCXEMrhLa2bn7+HX9HE6zjFrskABDwCjgU9V\ndZiInAKMb+pFqtroMaq6EzjSFNr7gSE+5DKNeG/5Dj5atYs7z+hLZrtWTscxxgRYboaLJ7/Mp6yy\nhqSEWKfjHBNfrkHVqGoxECWkKBFvAAAUN0lEQVQiUar6OTDUz7nMcSgur+LemasY0qU1E0+wrj1j\nIlFOZhp1bmXB5gbv9gl6vhSofSLSCvgKeFlEHsUz1NwEqXtmrqK8spZHLh5CTLQtfmxMJBqenkJ8\nTBSzQ/g6lC+/vcYBFcAvgA+BfOAcf4Yyx+6DFTt4f/kOfvbjXvRub6uiGBOpEmKjyeqeEtL3Q/ky\niu+AqrpVtRZ4H/int8vPBJm9B6r5v3dWMrBzMpNO6ul0HGOMw3IyXKzdWcae8iqnoxyTBguUiIwW\nkS9EZIaIDPNO+roS2CUijc2xZxzyh3dXUXqwhkcuGkKsde0ZE/EOLb8xN0SnPWrst9jjwJ+A6cBn\nwERV7QCcBDwYgGzmKHy8aifvLN3Oraf0ol/HZKfjGGOCwKDOrUmKjwnZefkaK1Axqvqxqr4O7FTV\neQCqujYw0Yyv9lVU87u3V9KvYzI/OaXROXyNMREkJjqKUT1TQ/Y6VGMFqv5iIgcP2xceC96Hifve\nW03JgWoeuWiwde0ZY/5HToaLrcUVFJZUOB3lqDX222yIiOwXkTJgsPfxoeeDApTPNOGztbuYsfhb\nfnJyBgM7t3Y6jjEmyORmuoDQXH6jwQKlqtGqmqyqSaoa43186Hlo3pYcZkoP1nDXjBX0aZ/ErT/s\n5XQcY0wQ6t2+Fa5WcczZGHrdfNYfFML++P5q9pRX88jFg4mLsf9KY8z3iQhjMlzMyS/Gs7Rf6LDf\naiHqy/VFvJZXyKSTejK4Sxun4xhjglhORhq7y6rILyp3OspRsQIVgsoqa7jzzeVktmvFz35kXXvG\nmMblZniuQ4XatEdWoELQn2atZdf+Sh65aDAJsdFOxzHGBLn0tJZ0SWkRcsPNrUCFmG827GH6ggIm\nntiTYekpTscxxoSInIw05uYXU+cOnetQVqBCSHlVLXe8uZyerkR+eWpvp+MYY0JIbqaL/ZW1rN6+\n3+koPrMCFUL+/MFatpce5GHr2jPGHKUx3nn5ZodQN58VqBAxN7+YqfO2cl1OD7K6pzodxxgTYtol\nJdCrXStmh9D9UFagQkBFtadrr1taS24/vY/TcYwxISo308XCLXuprnU3fXAQsAIVAh7+cB0Feyt4\n+MLBtIizrj1jzLEZk5FGZY2bJQUlTkfxid8KlIhMFpHd3nWkjrS/r4jMFZEqEfn1YfvGisg6Edko\nInf6K2MoWLhlLy/M3cI1Y7oxqmea03GMMSFsdM80ogRmh8i8fP5sQU0BGlvYcC9wG/CX+htFJBr4\nF3AG0B8YLyL9/ZQxqB2sruM3byynS0oLfjO2r9NxjDEhrnWLWAZ1bs3cEBko4bcCpapf4SlCDe3f\nraoLgZrDdmUDG1V1k6pWA68A4/yVM5j99eN1bN5zgD9fMJjE+Bin4xhjwsCYDBdLCvZxoKrW6ShN\nCsZrUJ2BbfWeF3q3HZGITBKRPBHJKyoq8nu4QFm2bR+TZ2/m8lHp5HinyzfGmOOVm5lGrVtZsKXB\n9kPQCMYCJUfY1uCtz6r6jKpmqWpW27Zt/RgrcGrq3Nw5YwWuVvHceYZ17Rljmk9Wt1TioqOYGwLX\noYKx36gQ6FrveRdgu0NZHDH5m82s2bGfp64cTnKCLb1ljGk+LeKiGZbeJiTuhwrGFtRCoJeI9BCR\nOOAyYKbDmQKmoLiCv3+6nlP7t+f0AR2cjmOMCUO5mS5W79hPyYFqp6M0yp/DzKcDc4E+IlIoIhNE\n5CYRucm7v4OIFAK/BO72HpOsqrXArcBHwBrgNVVd5a+cwURV+d3bK4gW4b5xAxA5Um+nMcYcn5yM\nNFRh3qbg7ubzWxefqo5vYv9OPN13R9o3C5jlj1zBbOay7Xy9YQ9/OHcAHVu3cDqOMSZMDenahsS4\naGbn7+GMQR2djtOgYOzii0glB6q5793VDO3ahitHd3M6jjEmjMVGR5HdI5U5QT5QwgpUkPjTrDWU\nHqzhwQsGER1lXXvGGP/KyXCxqegAO0srnY7SICtQQWDOxj28vqiQG07qSb+OyU7HMcZEgJxM7/Ib\nQTyazwqUwypr6vjtWyvoltaSn/2ol9NxjDERol+HZFJaxgZ1N58VKIc9/tlGthRX8MfzBtkihMaY\ngImKEsZkpDEnfw+qwbkMvBUoB63bWcZTX+ZzwfDOnNDLpjMyxgRWToaLHaWVbCmucDrKEVmBcojb\nrdw1YzlJCTHcfVZETtZujHFYTkZwX4eyAuWQlxcUsLhgH3ef1Z/UxDin4xhjIlAPVyIdWycwJ0iX\n37AC5YCdpZU8/MFacjPTuGB4gxO1G2OMX4kIORku5uYX43YH33UoK1AO+P3MVVTXufnjeYNsOiNj\njKNyMtIoqahhzc79Tkf5HitQAfbxqp18uGonP/txL7q7Ep2OY4yJcLne9ebmbAy+4eZWoAKorLKG\ne95ZRd8OSdxwYk+n4xhjDB1aJ9CzbWJQXoeyAhVAf/14PbvKKnnwgkHERtupN8YEh5yMNBZs3ktN\nndvpKP/DfksGyJKCEl6Yu4WrR3djWHqK03GMMeY7uRkuDlTXsWzbPqej/A8rUAFQU+fmrhkraJ+U\nwK9P7+N0HGOM+R+je6YhQtBNe2QFKgCe/XoTa3eWcd+4ASTZEu7GmCCTkhhH/47JQXfDrhUoP9uy\n5wCPfrqBsQM6cJot4W6MCVK5mS6WFOzjYHWd01G+YwXKjw4t4R4XHcXvzx3gdBxjjGnQmIw0quvc\n5G3d63SU71iB8qO3lnzL7I3F/OaMvnRoneB0HGOMaVB291RiooTZQXQ/lBUoP9l7oJr731vN8PQ2\nXJGd7nQcY4xpVGJ8DMPS2zA3iO6H8luBEpHJIrJbRFY2sF9E5DER2Sgiy0VkeL19dSKy1Psx018Z\n/emB91dTVlnLgxcMJsqWcDfGhIAxGS5WfFtK6cEap6MA/m1BTQHGNrL/DKCX92MS8GS9fQdVdaj3\n41z/RfSPbzbsYcbib7npBxn06ZDkdBxjjPFJbkYaboV5m4Kjm89vBUpVvwIau9o2DnhRPeYBbUSk\no7/yBMrBas8S7j1cidz6w0yn4xhjjM+GpaeQEBvF3CC5H8rJa1CdgW31nhd6twEkiEieiMwTkfMa\n+yQiMsl7bF5RUZG/svrssc82ULC3gj+eP9CWcDfGhJS4mChGdk8NmvuhnCxQR7owc2hBknRVzQIu\nB/4hIhkNfRJVfUZVs1Q1q23btv7I6bM1O/bz7FebuHhEF3IybAl3Y0zoyc10sWF3ObvLKp2O4miB\nKgS61nveBdgOoKqH/t0EfAEMC3S4o1XnVu6asYLkFrH89sx+TscxxphjcmgZ+GDo5nOyQM0ErvaO\n5hsNlKrqDhFJEZF4ABFxAbnAagdz+uSleVtZum0f95zdnxRbwt0YE6IGdGpNckJMUHTzxfjrE4vI\ndOBkwCUihcC9QCyAqj4FzALOBDYCFcB13pf2A54WETeeAvqQqgZ1gdq+7yAPf7iWE3u5GDe0k9Nx\njDHmmEVHCWMy0oJi4li/FShVHd/EfgVuOcL2OcAgf+VqbqrKPe+sok7VlnA3xoSFnAwXH63aRUFx\nBelpLR3LYTNJHKePVu3k0zW7+MWPezv6H2mMMc0lN9NzHWq2w7NKWIE6Dvsra7h35ir6d0xmwgk9\nnI5jjDHNIqNtK9olxTvezWcF6jg88uE6isqqePCCQcTYEu7GmDAhIuRkpDE3fw+eqzHOsN+qx2jR\n1r28NH8r1+R0Z0jXNk7HMcaYZpWT6WJPeTXrdpU5lsEK1DGorvUs4d4xOYFfnWZLuBtjws+h+6Hm\nOLj8hhWoY/DMV/ms31XOfeMG0irebwMhjTHGMV1SWtItrSVzHBwoYQXqKG0qKuexzzZy1qCO/Lh/\ne6fjGGOM3+RkuJi/aS+1dW5H3t8K1FFQVX731kriY6K495z+Tscxxhi/yslIo6yqlhXfljry/lag\njsIbiwqZu6mYO8/oS7tkW8LdGBPevrsO5dBwcytQPtpTXsUfZ60hq1sK40faEu7GmPCX1iqevh2S\nHLsOZQXKRw+8t5oDVbU8eMEgW8LdGBMxcjJc5G0pobKmLuDvbQXKB1+uL+Ltpdu5+eRMerW3JdyN\nMZEjNzONqlo3i7eWBPy9rUA14WB1HXe/vYKebRP5yckNrptojDFhKbtHKtFR4sh1KCtQTfjHf9az\nbe9B/nT+IFvC3RgTcZISYhncpbUjE8dagWrEqu2lPPf1Zi7N6sronmlOxzHGGEfkZrhYXlhKWWVN\nQN/XClQDDi3hntIylrvO7Ot0HGOMcUxORhp1bmXB5r0BfV8rUA14Yc4WlheWcs85A2jT0pZwN8ZE\nruHdUoiPiWJ2gOflswJ1BN/uO8hfPl7HyX3acs7gjk7HMcYYRyXERpPVPSXg90NZgTqMqnLP2ytR\nhfvHDbQl3I0xBs/9UGt3lrGnvCpg72kF6jAfrNzJf9bu5len9aZrqi3hbowx8N9pj+YGcLi5XwuU\niEwWkd0isrKB/SIij4nIRhFZLiLD6+27RkQ2eD+u8WfOQ0oPepZwH9g5mWtzugfiLY0xJiQM6tya\npPiYgN4P5e8W1BRgbCP7zwB6eT8mAU8CiEgqcC8wCsgG7hWRFL8mBf784VqKy6t48PzBtoS7McbU\nExMdxaieqQG9DuXX38Kq+hXQ2LjEccCL6jEPaCMiHYHTgU9Uda+qlgCf0HihO24Lt+xl2vwCrs/t\nwaAurf35VsYYE5JyMlxsLa6gsKQiIO/ndDOhM7Ct3vNC77aGtn+PiEwSkTwRySsqKjrmIE98vpHO\nbVrwi1N7H/PnMMaYcJaT6bkONW9TYO6Hcnq98iMNkdNGtn9/o+ozwDMAWVlZRzzGF09cMYJtJRUk\n2hLuxhhzRH3aJ/H+bSfQr0NyQN7P6RZUIdC13vMuwPZGtvtNi7hoettM5cYY0yARYUCn1gFbcsjp\nAjUTuNo7mm80UKqqO4CPgNNEJMU7OOI07zZjjDERwq/9WSIyHTgZcIlIIZ6RebEAqvoUMAs4E9gI\nVADXefftFZH7gYXeT3WfqgZ2EihjjDGO8muBUtXxTexX4JYG9k0GJvsjlzHGmODndBefMcYYc0RW\noIwxxgQlK1DGGGOCkhUoY4wxQckKlDHGmKAknoF04UFEioCtx/EpXEBgV+QKTnYePOw8eNh58LDz\n8F/Hey66qWrbpg4KqwJ1vEQkT1WznM7hNDsPHnYePOw8eNh5+K9AnQvr4jPGGBOUrEAZY4wJSlag\n/tczTgcIEnYePOw8eNh58LDz8F8BORd2DcoYY0xQshaUMcaYoGQFyhhjTFCKyAIlImNFZJ2IbBSR\nO4+wP15EXvXuny8i3QOf0v98OA+/FJHVIrJcRP4jIt2cyOlvTZ2HesddJCIqImE51NiX8yAil3i/\nJ1aJyLRAZwwEH34u0kXkcxFZ4v3ZONOJnP4mIpNFZLeIrGxgv4jIY97ztFxEhjd7CFWNqA8gGsgH\negJxwDKg/2HH/AR4yvv4MuBVp3M7dB5OAVp6H98cqefBe1wS8BUwD8hyOrdD3w+9gCVAivd5O6dz\nO3QengFu9j7uD2xxOrefzsVJwHBgZQP7zwQ+AAQYDcxv7gyR2ILKBjaq6iZVrQZeAcYddsw44AXv\n4zeAH4lIYNY4Dpwmz4Oqfq6qFd6n84AuAc4YCL58PwDcDzwMVAYyXAD5ch5uAP6lqiUAqro7wBkD\nwZfzoECy93FrYHsA8wWMqn4FNLZQ7DjgRfWYB7QRkY7NmSESC1RnYFu954XebUc8RlVrgVIgLSDp\nAseX81DfBDx/LYWbJs+DiAwDuqrqe4EMFmC+fD/0BnqLyGwRmSciYwOWLnB8OQ+/B670rhI+C/hp\nYKIFnaP9HXLU/LqibpA6Ukvo8LH2vhwT6nz+GkXkSiAL+IFfEzmj0fMgIlHA34FrAxXIIb58P8Tg\n6eY7GU9r+msRGaiq+/ycLZB8OQ/jgSmq+lcRGQNM9Z4Ht//jBRW//56MxBZUIdC13vMufL+J/t0x\nIhKDpxnfWFM3FPlyHhCRHwO/A85V1aoAZQukps5DEjAQ+EJEtuDpa58ZhgMlfP25eEdVa1R1M7AO\nT8EKJ76chwnAawCqOhdIwDN5aqTx6XfI8YjEArUQ6CUiPUQkDs8giJmHHTMTuMb7+CLgM/VeFQwj\nTZ4Hb9fW03iKUzheb4AmzoOqlqqqS1W7q2p3PNfizlXVPGfi+o0vPxdv4xk4g4i48HT5bQpoSv/z\n5TwUAD8CEJF+eApUUUBTBoeZwNXe0XyjgVJV3dGcbxBxXXyqWisitwIf4RmxM1lVV4nIfUCeqs4E\nnsfTbN+Ip+V0mXOJ/cPH8/AI0Ap43TtGpEBVz3UstB/4eB7Cno/n4SPgNBFZDdQBt6tqsXOpm5+P\n5+FXwLMi8gs8XVrXhuEfsIjIdDzduS7v9bZ7gVgAVX0Kz/W3M4GNQAVwXbNnCMPzaowxJgxEYhef\nMcaYEGAFyhhjTFCyAmWMMSYoWYEyxhgTlKxAGWOMCUpWoIwxxgQlK1DGBBkReU5E+jeyf4SIrPAu\nc/BYGE5kbAxg90EZE3JEZAHwMzyzWswCHlPVcJzI10Q4a0EZ4xAR6S4ia0XkBe+Cb2+ISEsR+aKh\nuf68yxkkq+pc7+wFLwLnBTS4MQFiBcoYZ/UBnlHVwcB+PItlNqYznkk6D2n2JQ6MCRZWoIxx1jZV\nne19/BJwQhPHR8JSMMYAVqCMcdrhxaWpYlPI/65s3OxLHBgTLKxAGeOsdO+id+BZCO+bxg72LmdQ\nJiKjvaP3rgbe8XNGYxxhBcoYZ60BrhGR5UAq8KQPr7kZeA7PMgf5gI3gM2HJhpkb4xAR6Q68p6oD\nHY5iTFCyFpQxxpigZC0oY4KUiMwH4g/bfJWqrnAijzGBZgXKGGNMULIuPmOMMUHJCpQxxpigZAXK\nGGNMULICZYwxJij9P5MjMQTJeWLeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a0d05c0c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NUM_SAMPLES = 5000\n",
    "N = 10\n",
    "pi_0s = np.linspace(0,1,10)\n",
    "\n",
    "R = np.zeros(pi_0s.size)\n",
    "\n",
    "for i in range(0,pi_0s.size):\n",
    "    this_pi0 = pi_0s[i]\n",
    "    this_pi1 = 1- this_pi0\n",
    "    threshold = np.log(((C_10-C_00)*this_pi0)/((C_01-C_11)*this_pi1))\n",
    "\n",
    "    for k in range(0,NUM_SAMPLES):\n",
    "        #draw a sample from H, and then from p(y|H)\n",
    "        H = np.random.choice(2,p=[this_pi0,this_pi1])\n",
    "        if(H==0):\n",
    "            y = np.random.multinomial(N,p_0)\n",
    "            L = 0\n",
    "            for j in range(0,len(dictionary)):\n",
    "                L += y[j]*np.log(p_1[j]/p_0[j])\n",
    "            if(L<threshold):\n",
    "                R[i] += C_00/NUM_SAMPLES\n",
    "            else:\n",
    "                R[i] += C_10/NUM_SAMPLES\n",
    "        else:\n",
    "            y = np.random.multinomial(N,p_1)\n",
    "            L = 0\n",
    "            for j in range(0,len(dictionary)):\n",
    "                L += y[j]*np.log(p_1[j]/p_0[j])\n",
    "            if(L<threshold):\n",
    "                R[i] += C_01/NUM_SAMPLES\n",
    "            else:\n",
    "                R[i] += C_11/NUM_SAMPLES  \n",
    "\n",
    "plt.plot(pi_0s,R)\n",
    "\n",
    "plt.xlabel(\"pi_0\")\n",
    "plt.ylabel(\"Bayes Risk R\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's worth fiddling both with the costs and the number of words N per sample. How do the affect the Bayes risk? Which prior is the worst-case prior, corresponding with the minimax detector? How does changing the costs impact this prior?"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}