|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "code", |
| 5 | + "execution_count": null, |
| 6 | + "metadata": {}, |
| 7 | + "outputs": [], |
| 8 | + "source": [ |
| 9 | + "import pandas as pd\n", |
| 10 | + "import janitor as jn\n", |
| 11 | + "import pymc3 as pm\n", |
| 12 | + "import matplotlib.pyplot as plt\n", |
| 13 | + "import seaborn as sns\n", |
| 14 | + "import numpy as np\n", |
| 15 | + "from utils import ECDF\n", |
| 16 | + "\n", |
| 17 | + "%load_ext autoreload\n", |
| 18 | + "%autoreload 2\n", |
| 19 | + "%matplotlib inline\n", |
| 20 | + "%config InlineBackend.figure_format = 'retina'" |
| 21 | + ] |
| 22 | + }, |
| 23 | + { |
| 24 | + "cell_type": "markdown", |
| 25 | + "metadata": {}, |
| 26 | + "source": [ |
| 27 | + "# Darwin's Finches\n", |
| 28 | + "\n", |
| 29 | + "A research group has taken measurements of the descendants of the finches that Charles Darwin observed when he postulated the theory of evolution.\n", |
| 30 | + "\n", |
| 31 | + "We will be using Bayesian methods to analyze this data, specifically answering the question of how quantitatively different two species of birds' beaks are.\n", |
| 32 | + "\n", |
| 33 | + "## Data Credits\n", |
| 34 | + "\n", |
| 35 | + "The Darwin's finches datasets come from the paper, [40 years of evolution. Darwin's finches on Daphne Major Island][data]. \n", |
| 36 | + "\n", |
| 37 | + "One row of data has been added for pedagogical purposes.\n", |
| 38 | + "\n", |
| 39 | + "[data]: (https://datadryad.org/resource/doi:10.5061/dryad.g6g3h). \n", |
| 40 | + "\n", |
| 41 | + "Let's get started and load the data." |
| 42 | + ] |
| 43 | + }, |
| 44 | + { |
| 45 | + "cell_type": "code", |
| 46 | + "execution_count": null, |
| 47 | + "metadata": {}, |
| 48 | + "outputs": [], |
| 49 | + "source": [ |
| 50 | + "from data import load_finches_2012\n", |
| 51 | + "df = load_finches_2012()" |
| 52 | + ] |
| 53 | + }, |
| 54 | + { |
| 55 | + "cell_type": "markdown", |
| 56 | + "metadata": {}, |
| 57 | + "source": [ |
| 58 | + "### Exercise\n", |
| 59 | + "\n", |
| 60 | + "View a random sample of the data to get a feel for the structure of the dataset." |
| 61 | + ] |
| 62 | + }, |
| 63 | + { |
| 64 | + "cell_type": "code", |
| 65 | + "execution_count": null, |
| 66 | + "metadata": {}, |
| 67 | + "outputs": [], |
| 68 | + "source": [ |
| 69 | + "# Your code below\n" |
| 70 | + ] |
| 71 | + }, |
| 72 | + { |
| 73 | + "cell_type": "markdown", |
| 74 | + "metadata": {}, |
| 75 | + "source": [ |
| 76 | + "**Note:** I have added one row of data, simulating the discovery of an \"unknown\" species of finch for which beak measurements have been taken.\n", |
| 77 | + "\n", |
| 78 | + "For pedagogical brevity, we will analyze only beak depth during the class. However, I would encourage you to perform a similar analysis for beak length as well." |
| 79 | + ] |
| 80 | + }, |
| 81 | + { |
| 82 | + "cell_type": "code", |
| 83 | + "execution_count": null, |
| 84 | + "metadata": {}, |
| 85 | + "outputs": [], |
| 86 | + "source": [ |
| 87 | + "# These are filters that we can use later on.\n", |
| 88 | + "fortis_filter = df['species'] == 'fortis'\n", |
| 89 | + "scandens_filter = df['species'] == 'scandens'\n", |
| 90 | + "unknown_filter = df['species'] == 'unknown'" |
| 91 | + ] |
| 92 | + }, |
| 93 | + { |
| 94 | + "cell_type": "markdown", |
| 95 | + "metadata": {}, |
| 96 | + "source": [ |
| 97 | + "### Exercise\n", |
| 98 | + "\n", |
| 99 | + "Recreate the estimation model for finch beak depths. A few things to note:\n", |
| 100 | + "\n", |
| 101 | + "- Practice using numpy-like fancy indexing.\n", |
| 102 | + "- Difference of means & effect size are optional." |
| 103 | + ] |
| 104 | + }, |
| 105 | + { |
| 106 | + "cell_type": "code", |
| 107 | + "execution_count": null, |
| 108 | + "metadata": {}, |
| 109 | + "outputs": [], |
| 110 | + "source": [ |
| 111 | + "with pm.Model() as beak_depth_model:\n", |
| 112 | + "\n", |
| 113 | + " # Your model defined here.\n" |
| 114 | + ] |
| 115 | + }, |
| 116 | + { |
| 117 | + "cell_type": "markdown", |
| 118 | + "metadata": {}, |
| 119 | + "source": [ |
| 120 | + "### Exercise\n", |
| 121 | + "\n", |
| 122 | + "Perform MCMC sampling to estimate the posterior distribution of each parameter." |
| 123 | + ] |
| 124 | + }, |
| 125 | + { |
| 126 | + "cell_type": "code", |
| 127 | + "execution_count": null, |
| 128 | + "metadata": {}, |
| 129 | + "outputs": [], |
| 130 | + "source": [ |
| 131 | + "# Your code below.\n" |
| 132 | + ] |
| 133 | + }, |
| 134 | + { |
| 135 | + "cell_type": "markdown", |
| 136 | + "metadata": {}, |
| 137 | + "source": [ |
| 138 | + "### Exercise\n", |
| 139 | + "\n", |
| 140 | + "Diagnose whether the sampling has converged or not using trace plots." |
| 141 | + ] |
| 142 | + }, |
| 143 | + { |
| 144 | + "cell_type": "code", |
| 145 | + "execution_count": null, |
| 146 | + "metadata": {}, |
| 147 | + "outputs": [], |
| 148 | + "source": [ |
| 149 | + "# Your code below.\n" |
| 150 | + ] |
| 151 | + }, |
| 152 | + { |
| 153 | + "cell_type": "markdown", |
| 154 | + "metadata": {}, |
| 155 | + "source": [ |
| 156 | + "### Exercise\n", |
| 157 | + "\n", |
| 158 | + "Visualize the posterior distribution over the parameters using the forest plot." |
| 159 | + ] |
| 160 | + }, |
| 161 | + { |
| 162 | + "cell_type": "code", |
| 163 | + "execution_count": null, |
| 164 | + "metadata": {}, |
| 165 | + "outputs": [], |
| 166 | + "source": [ |
| 167 | + "# Your code below.\n" |
| 168 | + ] |
| 169 | + }, |
| 170 | + { |
| 171 | + "cell_type": "markdown", |
| 172 | + "metadata": {}, |
| 173 | + "source": [ |
| 174 | + "### Exercise\n", |
| 175 | + "\n", |
| 176 | + "Visualize the posterior distribution of the means. " |
| 177 | + ] |
| 178 | + }, |
| 179 | + { |
| 180 | + "cell_type": "code", |
| 181 | + "execution_count": null, |
| 182 | + "metadata": {}, |
| 183 | + "outputs": [], |
| 184 | + "source": [ |
| 185 | + "# Your code below.\n" |
| 186 | + ] |
| 187 | + }, |
| 188 | + { |
| 189 | + "cell_type": "markdown", |
| 190 | + "metadata": {}, |
| 191 | + "source": [ |
| 192 | + "### Discuss\n", |
| 193 | + "\n", |
| 194 | + "- Is the posterior distribution of beaks for the unknown species reasonable?" |
| 195 | + ] |
| 196 | + }, |
| 197 | + { |
| 198 | + "cell_type": "markdown", |
| 199 | + "metadata": {}, |
| 200 | + "source": [ |
| 201 | + "### Exericse\n", |
| 202 | + "\n", |
| 203 | + "Perform a posterior predictive check to visually diagnose whether the model describes the data generating process well or not." |
| 204 | + ] |
| 205 | + }, |
| 206 | + { |
| 207 | + "cell_type": "code", |
| 208 | + "execution_count": null, |
| 209 | + "metadata": {}, |
| 210 | + "outputs": [], |
| 211 | + "source": [ |
| 212 | + "samples = pm.sample_ppc(trace, model=beak_depth_model, samples=2000)" |
| 213 | + ] |
| 214 | + }, |
| 215 | + { |
| 216 | + "cell_type": "markdown", |
| 217 | + "metadata": {}, |
| 218 | + "source": [ |
| 219 | + "Hint: Each column in the samples (key: \"likelihood\") corresponds to simulated measurements of each finch in the dataset. We can use fancy indexing along the columns (axis 1) to select out simulated measurements for each category, and then flatten the resultant array to get the full estimated distribution of values for each class." |
| 220 | + ] |
| 221 | + }, |
| 222 | + { |
| 223 | + "cell_type": "code", |
| 224 | + "execution_count": null, |
| 225 | + "metadata": {}, |
| 226 | + "outputs": [], |
| 227 | + "source": [ |
| 228 | + "fig = plt.figure()\n", |
| 229 | + "ax_fortis = fig.add_subplot(2, 1, 1)\n", |
| 230 | + "ax_scandens = fig.add_subplot(2, 1, 2, sharex=ax_fortis)\n", |
| 231 | + "\n", |
| 232 | + "# Extract just the fortis samples.\n", |
| 233 | + "\n", |
| 234 | + "# Compute the ECDF for the fortis samples.\n", |
| 235 | + "\n", |
| 236 | + "ax_fortis.plot(x_s, y_s, label='samples')\n", |
| 237 | + "\n", |
| 238 | + "# Extract just the fortis measurements.\n", |
| 239 | + "\n", |
| 240 | + "# Compute the ECDF for the fortis measurements.\n", |
| 241 | + "\n", |
| 242 | + "ax_fortis.plot(x, y, label='data')\n", |
| 243 | + "\n", |
| 244 | + "ax_fortis.legend()\n", |
| 245 | + "ax_fortis.set_title('fortis')\n", |
| 246 | + "\n", |
| 247 | + "# Extract just the scandens samples.\n", |
| 248 | + "\n", |
| 249 | + "# Compute the ECDF for the scandens samples\n", |
| 250 | + "\n", |
| 251 | + "ax_scandens.plot(x_s, y_s, label='samples')\n", |
| 252 | + "\n", |
| 253 | + "# Extract just the scandens measurements.\n", |
| 254 | + "\n", |
| 255 | + "# Compute the ECDF for the scanens measurements.\n", |
| 256 | + "\n", |
| 257 | + "\n", |
| 258 | + "ax_scandens.plot(x, y, label='data')\n", |
| 259 | + "ax_scandens.legend()\n", |
| 260 | + "ax_scandens.set_title('scandens')\n", |
| 261 | + "ax_scandens.set_xlabel('beak length')\n", |
| 262 | + "\n", |
| 263 | + "sns.despine()\n", |
| 264 | + "plt.tight_layout()" |
| 265 | + ] |
| 266 | + }, |
| 267 | + { |
| 268 | + "cell_type": "markdown", |
| 269 | + "metadata": {}, |
| 270 | + "source": [ |
| 271 | + "## Summary\n", |
| 272 | + "\n", |
| 273 | + "1. NumPy-like fancy indexing lets us write models in a concise fashion.\n", |
| 274 | + "1. Posterior estimates can show up as being \"unreasonable\", \"absurd\", or at the minimum, counter-intuitive, if we do not impose the right set of assumptions on the model.\n" |
| 275 | + ] |
| 276 | + }, |
| 277 | + { |
| 278 | + "cell_type": "code", |
| 279 | + "execution_count": null, |
| 280 | + "metadata": {}, |
| 281 | + "outputs": [], |
| 282 | + "source": [] |
| 283 | + } |
| 284 | + ], |
| 285 | + "metadata": { |
| 286 | + "kernelspec": { |
| 287 | + "display_name": "bayesian-modelling-tutorial", |
| 288 | + "language": "python", |
| 289 | + "name": "bayesian-modelling-tutorial" |
| 290 | + }, |
| 291 | + "language_info": { |
| 292 | + "codemirror_mode": { |
| 293 | + "name": "ipython", |
| 294 | + "version": 3 |
| 295 | + }, |
| 296 | + "file_extension": ".py", |
| 297 | + "mimetype": "text/x-python", |
| 298 | + "name": "python", |
| 299 | + "nbconvert_exporter": "python", |
| 300 | + "pygments_lexer": "ipython3", |
| 301 | + "version": "3.6.6" |
| 302 | + } |
| 303 | + }, |
| 304 | + "nbformat": 4, |
| 305 | + "nbformat_minor": 2 |
| 306 | +} |
0 commit comments