{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Think Bayes\n",
"\n",
"This notebook presents code and exercises from Think Bayes, second edition.\n",
"\n",
"Copyright 2016 Allen B. Downey\n",
"\n",
"MIT License: https://opensource.org/licenses/MIT"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Configure Jupyter so figures appear in the notebook\n",
"%matplotlib inline\n",
"\n",
"# Configure Jupyter to display the assigned value after an assignment\n",
"%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
"\n",
"import math\n",
"import numpy as np\n",
"\n",
"from thinkbayes2 import Pmf, Cdf, Suite\n",
"import thinkplot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The flea beetle problem\n",
"\n",
"Different species of flea beetle can be distinguished by the width and angle of the aedeagus. The data below includes measurements and know species classification for 74 specimens.\n",
"\n",
"Suppose you discover a new specimen under conditions where it is equally likely to be any of the three species. You measure the aedeagus and width 140 microns and angle 15 (in multiples of 7.5 degrees). What is the probability that it belongs to each species?\n",
"\n",
"\n",
"This problem is based on [this data story on DASL](https://web.archive.org/web/20160304083805/http://lib.stat.cmu.edu/DASL/Datafiles/FleaBeetles.html)\n",
"\n",
"Datafile Name: Flea Beetles\n",
"\n",
"Datafile Subjects: Biology\n",
"\n",
"Story Names: Flea Beetles\n",
"\n",
"Reference: Lubischew, A.A. (1962) On the use of discriminant functions in taxonomy. Biometrics, 18, 455-477. Also found in: Hand, D.J., et al. (1994) A Handbook of Small Data Sets, London: Chapman & Hall, 254-255.\n",
"\n",
"Authorization: Contact Authors\n",
"\n",
"Description: Data were collected on the genus of flea beetle Chaetocnema, which contains three species: concinna (Con), heikertingeri (Hei), and heptapotamica (Hep). Measurements were made on the width and angle of the aedeagus of each beetle. The goal of the original study was to form a classification rule to distinguish the three species.\n",
"\n",
"Number of cases: 74\n",
"\n",
"Variable Names:\n",
"\n",
"Width: The maximal width of aedeagus in the forpart (in microns)\n",
"\n",
"Angle: The front angle of the aedeagus (1 unit = 7.5 degrees)\n",
"\n",
"Species: Species of flea beetle from the genus Chaetocnema\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can read the data from this file:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" Width \n",
" Angle \n",
" Species \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 150 \n",
" 15 \n",
" Con \n",
" \n",
" \n",
" 1 \n",
" 147 \n",
" 13 \n",
" Con \n",
" \n",
" \n",
" 2 \n",
" 144 \n",
" 14 \n",
" Con \n",
" \n",
" \n",
" 3 \n",
" 144 \n",
" 16 \n",
" Con \n",
" \n",
" \n",
" 4 \n",
" 153 \n",
" 13 \n",
" Con \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Width Angle Species\n",
"0 150 15 Con\n",
"1 147 13 Con\n",
"2 144 14 Con\n",
"3 144 16 Con\n",
"4 153 13 Con"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"\n",
"df = pd.read_csv('../data/flea_beetles.csv', delimiter='\\t')\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's what the distributions of width look like."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def plot_cdfs(df, col):\n",
" for name, group in df.groupby('Species'):\n",
" cdf = Cdf(group[col], label=name)\n",
" thinkplot.Cdf(cdf)\n",
" \n",
" thinkplot.decorate(xlabel=col, \n",
" ylabel='CDF',\n",
" loc='lower right')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGSxJREFUeJzt3X+QXfV53/H3g1it1IADo1VMxIJXTaVicEEyi0IMYwvjH/wokhm7BeqOITWoakMDToqLwSYaZ/AQNB4DYwoVCqM4bgwUTFFSMXLMj2EgKNZikIIEeCSM0SKDhXDAGIQk8/SPvZIuq13t3d177j137/s1szP3/Ni9D2e0++H7vd/znMhMJEkqm4OaXYAkSUMxoCRJpWRASZJKyYCSJJWSASVJKiUDSpJUSgaUJKmUDChJUikZUJKkUjq42QWMVldXV/b09DS7DEnSGD3xxBOvZub0kc5ruYDq6emhr6+v2WVIksYoIn5Wy3lO8UmSSsmAkiSVkgElSSolA0qSVEoGlCSplAoLqIi4PSJ+ERFPD3M8IuKmiNgUEesj4sNF1SJJaj1FjqBWAGcc4PiZwKzK1yLglgJrkSS1mMLug8rMRyKi5wCnLAS+kwPPnF8TEYdFxO9m5s+LqkkajYuX77vfbvnFvU2sRFf/7bN7X197zjFNrKS9ffayW/e+vufGxYW/XzM/gzoS2FK13V/Zt5+IWBQRfRHRt23btoYUJ0lqrmYGVAyxL4c6MTOXZWZvZvZOnz5idwxJ0gTQzIDqB46q2u4GtjapFklSyTQzoFYCX6is5jsZeN3PnyRJexS2SCIivgfMB7oioh/4M6ADIDNvBVYBZwGbgLeAPyyqFklS6ylyFd8FIxxP4I+Ken9JUmtrucdtSPWyev3LrHxyK+/serfZpUgagq2O1LZqDafODn9NpGbwN09tq9ZwWjB3RgOqkTSYU3wSdopotkc3b+eBn7zKzt1D3gqpGt334DruvL+Pd3buanYpdeEISlLT1RpOkw8e6v5+7dGocOqc3FH4e4ABJakEag2n02d3NaCa1tWocDrvzMbMODjFJ6lUbAZbH41o5lo0A0oTlsvIpdbmFJ8mLJeRS63N30xNWC4jl1qbU3xqCy4jl1qPIyhJUikZUJKkUnKKT1Kh7BKhsTKg1LJcRt4aRhNOdooY3kRrY1QLp/jUslxG3hpGE052ihjeaMKpUa2IiuYISi3LZeStxy4RYzeacGpUK6KiGVCaEFxGrnYyEdoY1cK5D0lSKRlQkqRSMqAkSaVkQEmSSslFEpLGxRtxVRRHUJLGxce1qyiOoCSNi49rr007doIYLwNKUt14I+7w6hVOE6VLRC2c4pOkBqhXOE2ULhG1cASl0rIZrCaqdukEMV6OoFRaNoOV2pu/2Sotm8FK7c0pPrUEm8FK7ccRlCSplBxBSTogO0WoWQodQUXEGRHxXERsiogrhzh+dEQ8FBFPRsT6iDiryHokjZ6dItQshY2gImIScDPwSaAfWBsRKzNzY9VpXwXuysxbIuJYYBXQU1RNkkbPThG1sVNE/RU5xTcP2JSZzwNExB3AQqA6oBJ4X+X1bwNbC6xH0jjZKWJ4tYZTO3WCGK8iA+pIYEvVdj/w+4POWQL8ICL+G/BbwCcKrEeSClNrOLVTJ4jxKjKghpqQHjxXcAGwIjO/GRF/APx1RHwoM99zA0xELAIWARx99NGFFCtJ9WKniPoocpFEP3BU1XY3+0/hfRG4CyAzHwemAPtNZGfmsszszcze6dOnF1SuJKlMigyotcCsiJgZEZOB84GVg855ETgdICI+yEBAbSuwJklSiygsoDJzN3ApsBp4hoHVehsi4usRsaBy2p8Cl0TEOuB7wEWZ6c0WkqRib9TNzFUMLB2v3ndN1euNwClF1iDpwLwRV2VlqyOpzXkjrsrKgJLanDfiqqzsxSdpL2/EVZk4gpIklZIBJUkqJQNKklRKBpQkqZQMKElSKRlQkqRSMqAkSaVkQEmSSsmAkiSVkp0kJKkG9z24rubHuqs+HEFJUg1qDafOyR0NqKY9GFCSVINaw+m8M3sbUE17cIpPkkbpnhsXN7uEtuAISpJUSgaUJKmUDChJUin5GZSaYvX6l1n55Fbe2fVus0uRVFKOoNQUowmnzg7/mUrtyN98NcVowmnB3BkFVyOpjJziU9Mtv9j7RtRcdokoJ0dQktreaMLJThGNY0BJanujCSc7RTSOU3ySVMUuEeVhQKkQLiOXNF5O8akQtYaTS8glDce/DipEreHkEnJJw3GKT4VzGbmksXAEJUkqJQNKklRKBpQkqZT8DEpDcpm4pGYrdAQVEWdExHMRsSkirhzmnH8fERsjYkNE/E2R9ah29Qonl5FLGqvCRlARMQm4Gfgk0A+sjYiVmbmx6pxZwFeAUzLzlxHxO0XVo9GpVzi5jFzSWBU5xTcP2JSZzwNExB3AQmBj1TmXADdn5i8BMvMXBdajMXKZuKRmKHL+5UhgS9V2f2VftdnA7Ih4LCLWRMQZQ/2giFgUEX0R0bdt27aCypUklUmRARVD7MtB2wcDs4D5wAXA8og4bL9vylyWmb2Z2Tt9+vS6FypJKp8iA6ofOKpquxvYOsQ592Xmrsz8KfAcA4ElSWpzRQbUWmBWRMyMiMnA+cDKQef8X+A0gIjoYmDK7/kCa5IktYjCAiozdwOXAquBZ4C7MnNDRHw9IhZUTlsNbI+IjcBDwBWZub2omiRJraPQG3UzcxWwatC+a6peJ/AnlS9JKsR9D64b1WPdVQ7eRSlpwqs1nDondzSgGtXKgJI04dUaTued6T1/ZWIvPklt5Z4bFze7BNXIEZQkqZQcQUkT3KObt/PAT15l5+7B98lL5eYISprgag2nyQcP1fxFah4DSprgag2n02d3NaAaqXZO8Ult5Npzjml2CVLNDjiCiogVVa8vLLwaSZIqRpriO6Hq9WVFFiJJUrWRpvhc9jNBrV7/ct0e6y5JRRgpoLoj4iYGnu205/VemfnHhVWmQtUaTp0drqOR1BwjBdQVVa/7iixEjVVrOC2YO6MB1UjS/g4YUJn5V40qRM2z/GL7j0kqnxHnbyLiwoj4cUT8uvLVFxFfaERxkqT2dcARVCWILmfgeU0/ZuCzqA8DSyOCzPxO8SVKktrRSCOo/wqcm5kPZebrmfnPmfkg8NnKMUmSCjFSQL0vM18YvLOy731FFCRJEowcUG+P8ZgkSeMy0jLzD0bE+iH2B/AvC6hHkiRg5IA6AXg/sGXQ/g8AWwupSJIkRp7i+xbwRmb+rPoLeKtyTJKkQowUUD2Zud8UX2b2AT2FVCRJEiNP8U05wLGp9SxE0uj5OHdNZCMF1NqIuCQzb6veGRFfBJ4oriyNh53K28dowslHuqvVjBRQlwP3RsTn2RdIvcBk4NwiC9PYjSac7Fbe2kYTTj7SXa1mpGaxrwAfiYjTgA9Vdv+/SjcJldRowslu5ROHj3PXRDPSCAqAzHwIeKjgWlQAO5VLalXO70iSSqmmEZQkFem+B9dx5/19vLNzV7NLUYk4gpLUdI0Kp87JHYW/h+rHgJLUdI0Kp/PO9DPZVuIUn6RSuefGxc0uQSVR6AgqIs6IiOciYlNEXHmA8z4XERkR/u+NJAkoMKAiYhJwM3AmcCxwQUQcO8R5hwJ/DPxjUbVIklpPkSOoecCmzHw+M3cCdwALhzjvz4HrgR0F1iJJajFFBtSRvPc5Uv2VfXtFxFzgqMz8uwP9oIhYFBF9EdG3bdu2+lcqSSqdIgNqqM6UexuHRcRBDDxT6k9H+kGZuSwzezOzd/r06XUsUZJUVkWu4usHjqra7ua9T+E9lIH+fg9HBMARwMqIWFB53pSGYbdySe2gyBHUWmBWRMyMiMnA+cDKPQcz8/XM7MrMnszsAdYAhlMNag0nO5VLamWF/QXLzN3ApcBq4BngrszcEBFfj4gFRb1vO6g1nOxULqmVFXqjbmauAlYN2nfNMOfOL7KWicpu5ZImKueAJEmlZKsjqcQe3bx9VI91lyYSR1BSidUaTpMPHuquDqm1GVBSidUaTqfP7mpANVJjOcUntYhrzzmm2SVIDeUISpJUSgaUJKmUDChJUikZUJKkUjKgJEmlZEBJkkrJgJIklZIBJUkqJQNKklRKBpQkqZQMKElSKdmLr4RWr3+55se6S9JEZUCVUK3h1NnhAFjld9+D67jz/j7e2bmr2aWoxfgXroRqDacFc2c0oBppfEYTTp2TOwquRq3EEVTJLb+4t9klSOMymnA670z/vWsfA0pSw9xz4+Jml6AW4hSfJKmUDChJUik5xdcELiOXpJE5gmoCl5FL0sj8C9gELiOXpJE5xddkLiOXpKEZUJLGxU4RKopTfJLGpdZwskuERssRlNREj27ezgM/eZWdu7PZpYxZreFklwiNlgElNVGt4TT54GhANeNnpwjVk1N8UhPVGk6nz+5qQDVSuTiCkkri2nOOaXYJUqkUOoKKiDMi4rmI2BQRVw5x/E8iYmNErI+IByLiA0XWI0lqHYUFVERMAm4GzgSOBS6IiGMHnfYk0JuZxwN3A9cXVY8kqbUUOYKaB2zKzOczcydwB7Cw+oTMfCgz36psrgG6C6xHktRCivwM6khgS9V2P/D7Bzj/i8D9Qx2IiEXAIoCjjz66XvVJUtPt2rWL/v5+duzY0exS6m7KlCl0d3fT0TG2e+CKDKih1sUOuWQpIv4j0At8bKjjmbkMWAbQ29vbujeMSC3IThHF6u/v59BDD6Wnp4eI1ridoBaZyfbt2+nv72fmzJlj+hlFTvH1A0dVbXcDWwefFBGfAK4GFmTmOwXWI2kM7BRRrB07djBt2rQJFU4AEcG0adPGNTIsMqDWArMiYmZETAbOB1ZWnxARc4H/xUA4/aLAWiSNkZ0iijfRwmmP8f53FTbFl5m7I+JSYDUwCbg9MzdExNeBvsxcCSwFDgH+T+U/5MXMXFBUTZLGx04RaqRCb9TNzFXAqkH7rql6/Yki31+SVJuXX36Zyy+/nLVr19LZ2UlPTw833HADs2fPblpNtjqSpDaXmZx77rnMnz+fzZs3s3HjRr7xjW/wyiuvNLUuWx0VYPX6l2t+rLsk7fHZy24t7GcfaHr2oYceoqOjg8WL950zZ84cMpMrrriC+++/n4jgq1/9Kueddx4PP/wwS5Ysoauri6effpoTTzyR7373u3X/LM2AKkCt4dTZ4QBWUvPtCZnBvv/97/PUU0+xbt06Xn31VU466SQ++tGPAvDkk0+yYcMGZsyYwSmnnMJjjz3GqaeeWte6/AtZgFrDacHcGQ2oRpLG5tFHH+WCCy5g0qRJvP/97+djH/sYa9euBWDevHl0d3dz0EEHMWfOHF544YW6v78jqIItv9ilt5Jq06xVkscddxx33333fvszh++L0NnZuff1pEmT2L17d93rMqCkNmenCH384x/nqquu4rbbbuOSSy4BYO3atRx++OHceeedXHjhhbz22ms88sgjLF26lGeffbYhdRlQUpuzU4QignvvvZfLL7+c6667jilTpuxdZv7mm29ywgknEBFcf/31HHHEEQaUNF6Pbt5e8yPV25mdIgQwY8YM7rrrrv32L126lKVLl75n3/z585k/f/7e7W9/+9uF1GRAacJqpXCafHA5Wt3YKUJl4io+TVitFE6nz+5qdhlS6TiCUlu49pxjml2CpFFyBCVJKiUDSpJUSgaUJKmUDChJanOHHHLIe7ZXrFjBpZdeesDvWblyJdddd12RZblIQiozuzyorBYsWMCCBcU+X9aAUstqhxtxGxlOdopovouX9xX2s8faF3Tbtm0sXryYF198EYAbbriBU045hRUrVtDX11fYTbpgQKmF1RpOZbkJdiwaGU52imhfb7/9NnPmzNm7/dprr+0dHV122WV86Utf4tRTT+XFF1/k05/+NM8880xD6jKg1LJqDaeJchOsXR5UlKlTp/LUU0/t3d4zOgL44Q9/yMaNG/cee+ONN/jVr37VkLoMKE0I3oiriaCMj+d59913efzxx5k6dWrD39tVfJKkYX3qU596z+dM1SOtohlQkqRh3XTTTfT19XH88cdz7LHHcuuttzbsvZ3ik6Q29+abb75n+6KLLuKiiy4CoKurizvvvHO/76k+pyiOoCRJpWRASZJKySk+aRzs9CAVx4BSabVCp4hGhZNdHtSOnOJTabVCp4hGhZNdHtSO2m4EtXr9y6x8civv7Hq32aVoBK3WKcJOD1J9tV1ANTKcOjscoNaLnSKk4hxyyCHvWWreiEawtWi7v6CNDKcFc2c05L0kaSJquxFUtTL2vZLUvq7+22cL+9ljnYUY7nEbS5YsYfPmzbz00kts2bKFL3/5y1xyySX1LLm9A0qSNPbHbaxfv541a9bw61//mrlz53L22WczY0b9Zo4KDaiIOAO4EZgELM/M6wYd7wS+A5wIbAfOy8wXiqxJkvReY33cxsKFC5k6dSpTp07ltNNO40c/+hGf+cxn6lZXYQEVEZOAm4FPAv3A2ohYmZkbq077IvDLzPxXEXE+8BfAeUXVJEllVsbFQAd63EZEHHB7vIpcJDEP2JSZz2fmTuAOYOGgcxYCf1V5fTdwetT7v1CSNGYHetzGfffdx44dO9i+fTsPP/wwJ510Ul3fu8iAOhLYUrXdX9k35DmZuRt4HZg2+AdFxKKI6IuIvm3bthVUriRpsAM9bmPevHmcffbZnHzyyXzta1+r6+dPUOxnUEONhAbfeVnLOWTmMmAZQG9vb3n73khSCxrL4zYAZs+ezbJlywqrq8iA6geOqtruBrYOc05/RBwM/DbwWoE1ubS8hZRxPn4wu0dIxSkyoNYCsyJiJvAScD7wHwadsxK4EHgc+BzwYGY6QpKkkluyZEnh71FYQGXm7oi4FFjNwDLz2zNzQ0R8HejLzJXAXwJ/HRGbGBg5nV9UPZJUVplZ9xVwZTDe8Uah90Fl5ipg1aB911S93gH8uyJrkKQymzJlCtu3b2fatGkTKqQyk+3btzNlypQx/ww7SUhSE3V3d9Pf389EXKE8ZcoUuru7x/z9BpQkNVFHRwczZ85sdhml1HbdzCVJrcGAkiSVkgElSSqlaLXbjiJiG/CzZtdRB13Aq80uYgLwOtaH17E+vI61+UBmTh/ppJYLqIkiIvoy07YW4+R1rA+vY314HevLKT5JUikZUJKkUjKgmqe4FsDtxetYH17H+vA61pGfQUmSSskRlCSplAwoSVIpGVAFiIjbI+IXEfF01b6lEfFsRKyPiHsj4rCqY1+JiE0R8VxEfLo5VZfPMNfxzyvX8KmI+EFEzKjsj4i4qXId10fEh5tXebkMdR2rjv33iMiI6Kpsex2HMcy/xyUR8VLl3+NTEXFW1TF/r8fJgCrGCuCMQfv+HvhQZh4P/AT4CkBEHMvAc7COq3zP/4yISY0rtdRWsP91XJqZx2fmHODvgD2PbzkTmFX5WgTc0qgiW8AK9r+ORMRRwCeBF6t2ex2Ht4IhriPwrcycU/laBf5e14sBVYDMfIRBj67PzB9k5u7K5hpgTw/6hcAdmflOZv4U2ATMa1ixJTbMdXyjavO3gD2rfBYC38kBa4DDIuJ3G1NpuQ11HSu+BXyZfdcQvI7DOsB1HIq/13VgQDXHfwLur7w+EthSday/sk/DiIhrI2IL8Hn2jaC8jqMQEQuAlzJz3aBDXsfRu7QyHXp7RBxe2ed1rAMDqsEi4mpgN/C/9+wa4jTX/h9AZl6dmUcxcA0vrez2OtYoIv4FcDX7wv09h4fY53Uc3i3A7wFzgJ8D36zs9zrWgQHVQBFxIfBvgc/nvhvQ+oGjqk7rBrY2urYW9TfAZyuvvY61+z1gJrAuIl5g4Fr9OCKOwOs4Kpn5Smb+JjPfBW5j3zSe17EODKgGiYgzgP8BLMjMt6oOrQTOj4jOiJjJwIfTP2pGja0gImZVbS4Anq28Xgl8obIK7WTg9cz8ecMLbAGZ+U+Z+TuZ2ZOZPQz8Mf1wZr6M13FUBn0+dy6wZ4Wfv9d14CPfCxAR3wPmA10R0Q/8GQOr9jqBv48IgDWZuTgzN0TEXcBGBqb+/igzf9OcystlmOt4VkT8a+BdBh67srhy+irgLAY+jH4L+MOGF1xSQ13HzPzLYU73Og5jmH+P8yNiDgPTdy8A/xnA3+v6sNWRJKmUnOKTJJWSASVJKiUDSpJUSgaUJKmUDChJUikZUFJBIuJbEXF51fbqiFhetf3NiLgqIu4e5vsfjojeyuurqvb3DNWZXJpoDCipOP8AfAQgIg4Cuhjobr3HR4AHMvNzNfysq0Y+RZpYDCipOI9RCSgGgulp4FcRcXhEdAIfBH65ZzQUEVMj4o5K49E7gamV/dcBUyvPG9rTw3FSRNwWERsqz8Wa2tD/MqkBDCipIJm5FdgdEUczEFSPA/8I/AHQC6wHdlZ9y38B3qo8M+xa4MTKz7kSeLvyvKHPV86dBdycmccB/8y+noTShGFAScXaM4raE1CPV23/w6BzPwp8FyAz1zMQYMP5aWY+VXn9BNBTv5KlcjCgpGLt+Rzq3zAwxbeGgRHURxgIr8Fq7T32TtXr32BfTU1ABpRUrMcYeMTKa5XHMrwGHMZASD0+6NxHGHgIIxHxIeD4qmO7IqKjAfVKpWFAScX6JwZW760ZtO/1zHx10Lm3AIdExHoGHsVe/XiGZcD6qkUS0oRnN3NJUik5gpIklZIBJUkqJQNKklRKBpQkqZQMKElSKRlQkqRSMqAkSaX0/wHLQbsU8JgsLAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot_cdfs(df, 'Width')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the distributions of angle."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGCxJREFUeJzt3X2QXXWd5/H3Nw8kEVCZJKiho52ZJaVIQdAmw05YjICKMCayOgvRqYFdgcrOMoIzi+sUymad0mJgtgRKdpiYseLDDA8LMkQniCMPw8aSMa0EhoQHAyJpUGwioqghxHz3j3uTuenc7oZOnz6/m/t+VXX1Pef++t5POul8+nfuuecXmYkkSaWZVHcASZLasaAkSUWyoCRJRbKgJElFsqAkSUWyoCRJRbKgJElFsqAkSUWyoCRJRZpSd4CXa9asWdnb21t3DEnSGH33u999JjNnjzau4wqqt7eX/v7+umNIksYoIn74UsZ5iE+SVCQLSpJUJAtKklQkC0qSVCQLSpJUpMoKKiI+HxE/iYgHhrk/IuKqiNgcEfdHxFuqyiJJ6jxVzqBWA6eMcP+7gcObH+cBf11hFklSh6nsfVCZeXdE9I4wZCnwxWysOX9PRLw6Il6XmT+qKpNUuou/+tDu2596zxtrTCLtbcFHbtp9e8Nn3lf589X5GtRhwJaW7YHmvr1ExHkR0R8R/YODgxMSTpJUrzoLKtrsy3YDM3NlZvZlZt/s2aNeHUOStB+os6AGgLkt2z3AUzVlkSQVps6CWgP8UfNsvuOA53z9SZK0S2UnSUTEtcBiYFZEDAD/E5gKkJnXAGuBU4HNwK+A/1xVFklS56nyLL5lo9yfwH+r6vklSZ3NK0lIkopkQUmSimRBSZKK1HEr6kpjse7Rrdz+yDNs39H2rXZSrW654z6uv7WfF7a/WHeUkU2a2PehOoNSV+i0cjpgSrv3sWt/1RHl1GJSTMy/TwtKXaHTyumk+bPqjqEJ1GnlNOfQV0/Ic3mIT13Hi7CqZDddubzuCMM6Z1X/hD6fMyhJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkabUHUCSut3PmMGzkw5kJ8E5q/rrjlMMZ1CSVLNd5dQppk2dmOqwoCSpZp1WTkuOmTMhz+UhPkkqyKpz+uqOUAxnUJKkIllQkqQiVVpQEXFKRDwcEZsj4mNt7n99RNwZEfdGxP0RcWqVeSRJnaOygoqIycDVwLuBI4BlEXHEkGEfB27IzGOAM4H/U1UeSVJnqXIGtRDYnJmPZeZ24Dpg6ZAxCbyyeftVwFMV5pEkdZAqz+I7DNjSsj0A/O6QMSuAb0TEnwAHAidXmEeS1EGqnEG1O7E/h2wvA1ZnZg9wKvCliNgrU0ScFxH9EdE/ODhYQVRJUmmqLKgBYG7Ldg97H8L7EHADQGZ+G5gOzBr6QJm5MjP7MrNv9uzZFcWVJJWkyoJaDxweEfMi4gAaJ0GsGTLmCeAkgIh4E42CcookSaquoDJzB3A+cBvwII2z9TZGxCcjYklz2J8B50bEfcC1wNmZOfQwoCSpC1V6qaPMXAusHbLvkpbbm4BFVWaQJHUmryQhSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSqSBSVJKpIFJUkqkgUlSSrSlLoDSFJVbrnjPq6/tZ8Xtr9Yd5SRTZpdd4IiOYOStN/qiHJqMSmi7ghFsaAk7bc6rZzmHPrqumMUxUN8krrCTVcurzvCsM5Z1V93hCI5g5IkFcmCkiQVyUN82ifrHt3K7Y88w/YdWXcUSfsZZ1DaJ51WTgdM8SwpqVNYUNonnVZOJ82fVXcMSS+Rh/g0bj71njfWHUHSfsQZlCSpSBaUJKlIlRZURJwSEQ9HxOaI+NgwY/5TRGyKiI0R8fdV5pEkdY7KXoOKiMnA1cA7gAFgfUSsycxNLWMOB/4cWJSZz0bEoVXlkSR1lipnUAuBzZn5WGZuB64Dlg4Zcy5wdWY+C5CZP6kwjySpg1RZUIcBW1q2B5r7Ws0H5kfEtyLinog4pd0DRcR5EdEfEf2Dg4MVxZUklaTKgmr3jsihb5qZAhwOLAaWAasiYq/L+Wbmyszsy8y+2bNdN0WSukGVBTUAzG3Z7gGeajPmlsx8MTN/ADxMo7AkSV2uyjfqrgcOj4h5wJPAmcAHhoz5Bxozp9URMYvGIb/HKswkqYv8jBk8O+lAdhIuadGBKptBZeYO4HzgNuBB4IbM3BgRn4yIJc1htwFbI2ITcCdwUWZurSqTpO6yq5w6xbSpvjW1VaWXOsrMtcDaIfsuabmdwJ82PyRpXHVaOS05Zk7dMYritfgkdYVV5/TVHUEvk/NJSVKRLChJUpEsKElSkSwoSVKRLChJUpEsKElSkSwoSVKRRiyoiFjdcvusytNIktQ02gzq6JbbF1QZRJKkVqMV1NDlMSRJmhCjXeqoJyKuorG2067bu2XmhytLJknqaqMV1EUtt71WvSRpwoxYUJn5hYkKIklSq1FPM4+IsyLiexHxy+ZHf0T80USEkyR1rxFnUM0iupDGek3fo/Fa1FuAyyOCzPxi9RElSd1otBnUHwOnZ+admflcZv4sM+8A3te8T5KkSoxWUK/MzMeH7mzue2UVgSRJgtEL6tdjvE+SpH0y2mnmb4qI+9vsD+C3K8gjSRIwekEdDbwG2DJk/xuApypJJEkSox/i+wzw88z8YesH8KvmfZIkVWK0gurNzL0O8WVmP9BbSSJJkhi9oKaPcN+M8QwiSVKr0QpqfUScO3RnRHwI+G41kSRJGv0kiQuBmyPig/xbIfUBBwCnVxlMktTdRrtY7NPA70XE24Ejm7v/sXk1CUmSKjPaDAqAzLwTuLPiLJIk7Tbq1cwlSaqDBSVJKpIFJUkqkgUlSSqSBSVJKlKlBRURp0TEwxGxOSI+NsK490dERkRflXkkSZ2jsoKKiMnA1cC7gSOAZRFxRJtxBwMfBv6lqiySpM5T5QxqIbA5Mx/LzO3AdcDSNuP+ArgM2FZhFklSh6myoA5jz3WkBpr7douIY4C5mfm1kR4oIs6LiP6I6B8cHBz/pJKk4lRZUNFmX+6+M2ISjTWl/my0B8rMlZnZl5l9s2fPHseIkqRSVVlQA8Dclu0e9lyF92Aa1/e7KyIeB44D1niihCQJqi2o9cDhETEvIg4AzgTW7LozM5/LzFmZ2ZuZvcA9wJLmYoiSpC5XWUFl5g7gfOA24EHghszcGBGfjIglVT2vJGn/8JKuZj5WmbkWWDtk3yXDjF1cZRZJUmfxShKSpCJZUJKkIllQkqQiWVCSpCJZUJKkIllQkqQiWVCSpCJZUJKkIllQkqQiWVCSpCJZUJKkIllQkqQiWVCSpCJZUJKkIllQkqQiWVCSpCJZUJKkIlW6oq6k/dMtd9zH9bf288L2F+uOMrJJs+tOoH3gDErSy9YR5dRiUkTdETQGFpSkl63TymnOoa+uO4bGwEN8hVr36FZuf+QZtu/IuqNII7rpyuV1RxjWOav6646gfeAMqlCdVk4HTPEQiqTxZUEVqtPK6aT5s+qOIWk/4yG+DvCp97yx7gjSHn7GDJ6ddCA7CQ+jqTLOoCS9bLvKqVNMm+p/dZ3IvzVJL1unldOSY+bUHUNj4CE+Sftk1Tl9dUfQfsoZlCSpSBaUJKlIFpQkqUgWlCSpSBaUJKlIFpQkqUiVFlREnBIRD0fE5oj4WJv7/zQiNkXE/RFxe0S8oco8kqTOUVlBRcRk4Grg3cARwLKIOGLIsHuBvsw8CrgRuKyqPJKkzlLlDGohsDkzH8vM7cB1wNLWAZl5Z2b+qrl5D9BTYR5JUgep8koShwFbWrYHgN8dYfyHgFvb3RER5wHnAbz+9a8fr3ySVLsXX3yRgYEBtm3bVneUcTd9+nR6enqYOnXqmL6+yoJqd7GutmtIRMQfAn3A29rdn5krgZUAfX19nbMOhSSNYmBggIMPPpje3l5iP1qaPjPZunUrAwMDzJs3b0yPUeUhvgFgbst2D/DU0EERcTJwMbAkM1+oMI8kFWfbtm3MnDlzvyongIhg5syZ+zQzrLKg1gOHR8S8iDgAOBNY0zogIo4B/oZGOf2kwiySVKz9rZx22dc/V2UFlZk7gPOB24AHgRsyc2NEfDIiljSHXQ4cBPzfiNgQEWuGeThJUpepdLmNzFwLrB2y75KW2ydX+fySpJfmxz/+MRdeeCHr169n2rRp9Pb2csUVVzB//vzaMnklCUnqcpnJ6aefzuLFi3n00UfZtGkTn/70p3n66adrzeWChZJUiPddcE1lj33TlcuHve/OO+9k6tSpLF/+b2MWLFhAZnLRRRdx6623EhF8/OMf54wzzuCuu+5ixYoVzJo1iwceeIC3vvWtfPnLXx7319IsKEnqcrtKZqivfOUrbNiwgfvuu49nnnmGY489lhNOOAGAe++9l40bNzJnzhwWLVrEt771LY4//vhxzeUhPklSW+vWrWPZsmVMnjyZ17zmNbztbW9j/fr1ACxcuJCenh4mTZrEggULePzxx8f9+Z1BSVIhRjoMV6U3v/nN3HjjjXvtzxz+ugjTpk3bfXvy5Mns2LFj3HM5g5KkLnfiiSfywgsv8LnPfW73vvXr13PIIYdw/fXX85vf/IbBwUHuvvtuFi5cOGG5nEFJUpeLCG6++WYuvPBCLr30UqZPn777NPPnn3+eo48+mojgsssu47WvfS0PPfTQhOSyoCRJzJkzhxtuuGGv/ZdffjmXX375HvsWL17M4sWLd29/9rOfrSSTh/gkSUWyoCRJRbKgJElFsqAkSUWyoCRJRbKgJElFsqAkqcsddNBBe2yvXr2a888/f8SvWbNmDZdeemmVsXwflCTp5VuyZAlLliwZfeA+sKAkqRDnrOqv7LFXndM3pq8bHBxk+fLlPPHEEwBcccUVLFq0iNWrV9Pf31/Zm3TBgpKkrvfrX/+aBQsW7N7+6U9/unt2dMEFF/CRj3yE448/nieeeIJ3vetdPPjggxOSy4KSpC43Y8YMNmzYsHt71+wI4Jvf/CabNm3afd/Pf/5zfvGLX0xILgtKkgox1sNwVdq5cyff/va3mTFjxoQ/t2fxSZKG9c53vnOP15laZ1pVs6AkScO66qqr6O/v56ijjuKII47gmmuumbDn9hCfJHW5559/fo/ts88+m7PPPhuAWbNmcf311+/1Na1jquIMSpJUJAtKklQkC0qSVKSuew1q3aNbuf2RZ9i+I+uOIkkaQdfNoDqtnA6YEnVHkKRadN0MqtPK6aT5s+qOoQl02/0/Zs29T/HCizvrjiLVrusKqtWn3vPGuiNIe+i0cppE5/zCp+EddNBBe5xqPhEXgn0puu4Qn1SyTiunQ3b+su4Y2o919QxKKlmJ12Xb5X0XVLcsRDe7+KsPVfbYYz1iNNxyGytWrODRRx/lySefZMuWLXz0ox/l3HPPHc/IFpQkdbuxLrdx//33c8899/DLX/6SY445htNOO405c+aMW65KCyoiTgGuBCYDqzLz0iH3TwO+CLwV2AqckZmPV5lJkrSnsS63sXTpUmbMmMGMGTN4+9vfzne+8x3e+973jluuygoqIiYDVwPvAAaA9RGxJjM3tQz7EPBsZv67iDgT+EvgjKoySVLJSjxxa6TlNiJixO19VeVJEguBzZn5WGZuB64Dlg4ZsxT4QvP2jcBJMd5/QknSmI203MYtt9zCtm3b2Lp1K3fddRfHHnvsuD53lQV1GLClZXugua/tmMzcATwHzBz6QBFxXkT0R0T/4OBgRXElSUONtNzGwoULOe200zjuuOP4xCc+Ma6vP0G1r0G1mwkNfdPESxlDZq4EVgL09fX5xgtJGkdjWW4DYP78+axcubKyXFUW1AAwt2W7B3hqmDEDETEFeBXw0wozFXmMV9ql5FPLW9105fK6I6gLVFlQ64HDI2Ie8CRwJvCBIWPWAGcB3wbeD9yRmc6QJKlwK1asqPw5KiuozNwREecDt9E4zfzzmbkxIj4J9GfmGuBvgS9FxGYaM6czq8ojSaXKzHE/A64E+zrfqPR9UJm5Flg7ZN8lLbe3AX9QZQZJKtn06dPZunUrM2fO3K9KKjPZunUr06dPH/NjeCUJSapRT08PAwMD7I9nKE+fPp2enp4xf70FJUk1mjp1KvPmzas7RpG8mrkkqUgWlCSpSBaUJKlI0WlvO4qIQeCH+/gws4BnxiFO1cw5vsw5vsw5vrop5xsyc/ZogzquoMZDRPRnZvFv2Tfn+DLn+DLn+DLn3jzEJ0kqkgUlSSpStxZUdZffHV/mHF/mHF/mHF/mHKIrX4OSJJWvW2dQkqTCWVCSpCJ1VUFFxEciYmNEPBAR10bE2C+zW6GIuKCZcWNEXFh3nlYR8fmI+ElEPNCy77ci4p8i4vvNz4fUmbGZqV3OP2h+T3dGRBGn8w6T8/KIeCgi7o+ImyPi1XVmbGZql/Mvmhk3RMQ3ImJ81/seg3Y5W+777xGRETGrjmxDsrT7fq6IiCeb388NEXFqnRmbmdp+PyPiTyLi4ebP02VVPX/XFFREHAZ8GOjLzCNprFFV3PpTEXEkcC6wEDga+P2IOLzeVHtYDZwyZN/HgNsz83Dg9uZ23Vazd84HgP8I3D3haYa3mr1z/hNwZGYeBTwC/PlEh2pjNXvnvDwzj8rMBcDXgEv2+qqJt5q9cxIRc4F3AE9MdKBhrKZNTuAzmbmg+bG2zf0TbTVDckbE24GlwFGZ+Wbgr6p68q4pqKYpwIzm8vKvYO8l6EvwJuCezPxVZu4A/hk4veZMu2Xm3TQWl2y1FPhC8/YXgPdOaKg22uXMzAcz8+GaIrU1TM5vNP/uAe4Bxr5ewTgZJufPWzYPBGo/42qYf58AnwE+SgEZYcScRRkm538FLs3MF5pjflLV83dNQWXmkzSa/gngR8BzmfmNelO19QBwQkTMjIhXAKcCc2vONJrXZOaPAJqfD605z/7kvwC31h1iOBHxqYjYAnyQMmZQe4mIJcCTmXlf3VlegvObh00/X8Kh8mHMB/5DRPxLRPxzRBxb1RN1TUE1/7KXAvOAOcCBEfGH9abaW2Y+CPwljcM8XwfuA3aM+EXaL0XExTT+7v+u7izDycyLM3MujYzn151nqOYveRdTaHkO8dfA7wALaPwS/b/rjTOsKcAhwHHARcANUdFSwF1TUMDJwA8yczAzXwS+AvxezZnaysy/zcy3ZOYJNKbX36870yiejojXATQ/Vzbl7xYRcRbw+8AHszPerPj3wPvqDtHG79D4pfS+iHicxuHS70XEa2tN1UZmPp2Zv8nMncDnaLwOXaIB4CvZ8B1gJ40LyI67biqoJ4DjIuIVzbY/CXiw5kxtRcShzc+vp/Gi/rX1JhrVGuCs5u2zgFtqzNLxIuIU4H8ASzLzV3XnGc6Qk3eWAA/VlWU4mfmvmXloZvZmZi+N/1zfkpk/rjnaXnb9ktd0Oo3D/SX6B+BEgIiYDxxAVVdhz8yu+QD+F40fogeALwHT6s40TM7/B2yicXjvpLrzDMl2LY3DDy/S+GH/EDCTxtl7329+/q1Cc57evP0C8DRwW6E5NwNbgA3Nj2sKzXlT82fpfuCrwGEl5hxy/+PArBJzNv9P+tfm93MN8LpCcx4AfLn5d/894MSqnt9LHUmSitRNh/gkSR3EgpIkFcmCkiQVyYKSJBXJgpIkFcmCkioWEac3r6L9xn14jLMj4rPjmUsqnQUlVW8ZsI4Cr54vlcyCkioUEQcBi2i8wfHM5r7FEXFXRNzYXPfp73ZdyywiTm3uWxcRV0XE19o85uyIuCki1jc/Fk3oH0qaIFPqDiDt594LfD0zH4mIn0bEW5r7jwHeTGPJl28BiyKiH/gb4ITM/EFEDHeJqytprBu0rnk5rNtoLNMi7VcsKKlay4Armreva27/I/CdzBwAiIgNQC/wPPBYZv6gOf5a4Lw2j3kycETLBaRfGREHZ+YvKvkTSDWxoKSKRMRMGhfVPDIiksYqzgmspXE9wF1+Q+Nn8aUuWTAJ+PeZ+etxjCsVx9egpOq8H/hiZr4hG1fTngv8ADh+mPEPAb8dEb3N7TOGGfcNWtZeiogF4xNXKosFJVVnGXDzkH03AR9oN7g5I/pj4OsRsY7GFdefazP0w0Bfc+XVTcDy8YsslcOrmUsFiYiDMvP55ll9VwPfz8zP1J1LqoMzKKks5zZPmtgIvIrGWX1SV3IGJUkqkjMoSVKRLChJUpEsKElSkSwoSVKRLChJUpH+PxlIAscDcZaFAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot_cdfs(df, 'Angle')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I'll group the data by species and compute summary statistics."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"grouped = df.groupby('Species')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here are the means."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" Width \n",
" Angle \n",
" \n",
" \n",
" Species \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" Con \n",
" 146.190476 \n",
" 14.095238 \n",
" \n",
" \n",
" Hei \n",
" 124.645161 \n",
" 14.290323 \n",
" \n",
" \n",
" Hep \n",
" 138.272727 \n",
" 10.090909 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Width Angle\n",
"Species \n",
"Con 146.190476 14.095238\n",
"Hei 124.645161 14.290323\n",
"Hep 138.272727 10.090909"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"means = grouped.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the standard deviations."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" Width \n",
" Angle \n",
" \n",
" \n",
" Species \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" Con \n",
" 5.626891 \n",
" 0.889087 \n",
" \n",
" \n",
" Hei \n",
" 4.622758 \n",
" 1.101319 \n",
" \n",
" \n",
" Hep \n",
" 4.142484 \n",
" 0.971454 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Width Angle\n",
"Species \n",
"Con 5.626891 0.889087\n",
"Hei 4.622758 1.101319\n",
"Hep 4.142484 0.971454"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stddevs = grouped.std()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the correlations."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Con -0.193701151757636\n",
"Hei -0.06420611481268008\n",
"Hep -0.12478515405529574\n"
]
}
],
"source": [
"for name, group in grouped:\n",
" corr = group.Width.corr(group.Angle)\n",
" print(name, corr)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Those correlations are small enough that we can get an acceptable approximation by ignoring them, but we might want to come back later and write a complete solution that takes them into account.\n",
"\n",
"### The likelihood function\n",
"\n",
"To support the likelihood function, I'll make a dictionary for each attribute that contains a `norm` object for each species."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import norm\n",
"\n",
"dist_width = {}\n",
"dist_angle = {}\n",
"for name, group in grouped:\n",
" dist_width[name] = norm(group.Width.mean(), group.Width.std())\n",
" dist_angle[name] = norm(group.Angle.mean(), group.Angle.std())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can write the likelihood function concisely."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"class Beetle(Suite):\n",
" \n",
" def Likelihood(self, data, hypo):\n",
" \"\"\"\n",
" data: sequence of width, height\n",
" hypo: name of species\n",
" \"\"\"\n",
" width, angle = data\n",
" name = hypo\n",
" \n",
" like = dist_width[name].pdf(width)\n",
" like *= dist_angle[name].pdf(angle)\n",
" return like"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hypotheses are the species names:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"dict_keys(['Con', 'Hei', 'Hep'])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"hypos = grouped.groups.keys()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll start with equal priors"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Con 0.3333333333333333\n",
"Hei 0.3333333333333333\n",
"Hep 0.3333333333333333\n"
]
}
],
"source": [
"suite = Beetle(hypos)\n",
"suite.Print()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can update with the data and print the posterior."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Con 0.9902199258865487\n",
"Hei 0.009770186966082915\n",
"Hep 9.887147368342703e-06\n"
]
}
],
"source": [
"suite.Update((140, 15))\n",
"suite.Print()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Based on these measurements, the specimen is very likely to be an example of *Chaetocnema concinna*."
]
}
],
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}