Exercise 3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise 3: Bayesian computation with Stan and Farmer Jöns\n",
    "============\n",
    "### *Rasmus Bååth (adapted for Python by Christophe Carvenius)*\n",
    "Here follows a number of data analytic questions. Use [Stan](http://mc-stan.org/) and Python to build models that probe these questions. The Stan documentation can be found here: http://mc-stan.org/documentation/ . You can find the answers to the exercise questions here: https://github.com/rasmusab/bayesianprobabilitiesworkshop/blob/master/Exercise%203%20-%20solutions.ipynb\n",
    "\n",
    "## 1. Getting started\n",
    "Below is a code scaffold you can use. Right now the scaffold contains a simple model for two binomial rates, but this should be replaced with a model that matches the relevant questions.\n",
    "\n",
    "**→ Read through the code to see if you can figure out what does what and then run it to make sure it works. It should print out some statistics and some pretty graphs.**\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import pystan # install with pip install pystan\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# The Stan model as a string.\n",
    "model_string = \"\"\"\n",
    "data {\n",
    "  # Number of data points\n",
    "  int n1;\n",
    "  int n2;\n",
    "  # Number of successes\n",
    "  int y1[n1];\n",
    "  int y2[n2];\n",
    "}\n",
    "\n",
    "parameters {\n",
    "  real<lower=0, upper=1> theta1;\n",
    "  real<lower=0, upper=1> theta2;\n",
    "}\n",
    "\n",
    "model {  \n",
    "  theta1 ~ beta(1, 1);\n",
    "  theta2 ~ beta(1, 1);\n",
    "  y1 ~ bernoulli(theta1);\n",
    "  y2 ~ bernoulli(theta2); \n",
    "}\n",
    "\n",
    "generated quantities {\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "y1 = [0, 1, 0, 0, 0, 0, 1, 0, 0, 0]\n",
    "y2 = [0, 0, 1, 1, 1, 0, 1, 1, 1, 0]\n",
    "data_list = {'y1' : y1, 'y2' : y2, 'n1' : len(y1), 'n2' : len(y2)}\n",
    "\n",
    "# Compiling and producing posterior samples from the model.\n",
    "stan_samples = pystan.stan(model_code = model_string, data = data_list)\n",
    "\n",
    "# Plotting and summarizing the posterior distribution\n",
    "print(stan_samples)\n",
    "print(stan_samples.plot())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint 1:** `bernoulli` is the Bernoulli distribution which is the special case of the binomial distribution when there is just one trial `(bernoulli(x) === binomial(1, x))`. That is, if the data is coded as 4 successes out of 6 `(x = 4; n = 6)` it would be most convenient to use a binomial distribution. If the data is coded like `[1, 1, 1, 1, 0, 0]` it would be more convenient to use a Bernoulli distribution. The result would in any case be the same.\n",
    "\n",
    "**Hint 2:** Stan has quite a lot of different built in data types and two that sounds the same, but aren’t, are vectors and arrays. Vectors are simple, they are lists of real numbers and `vector[4] v;` would define a vector or length 4. Arrays are more general in that they can contain other data types, for example `int a[4]` would define an array of integers of length 4. Note the different placement of the `[]`-brackets compared to defining a vector.\n",
    "\n",
    "**Hint 3:** When defining parameters it’s important to properly define the *support*, that is, for what values the parameter has a defined meaning. For example, the support of a mean is on the whole real line (-Inf to Inf) so that can simply be declared by `real mu;`. A *standard deviation*, on the other hand, can’t be below 0.0, which could be written like this: `real<lower=0> sigma`. Finaly, a rate has to be between 0 and 1 which would be written like `real<lower=0, upper=1> theta;`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Manipulating samples\n",
    "To inspect and manipulate samples from individual parameters it is useful to convert the Stan “object” into a simple dataframe which gets one column per parameter: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta1</th>\n",
       "      <th>theta2</th>\n",
       "      <th>lp__</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.338945</td>\n",
       "      <td>0.673828</td>\n",
       "      <td>-15.336131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.125290</td>\n",
       "      <td>0.825161</td>\n",
       "      <td>-17.500812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.110733</td>\n",
       "      <td>0.596067</td>\n",
       "      <td>-15.812474</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.536350</td>\n",
       "      <td>0.324687</td>\n",
       "      <td>-18.623693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.174577</td>\n",
       "      <td>0.282452</td>\n",
       "      <td>-17.472201</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     theta1    theta2       lp__\n",
       "0  0.338945  0.673828 -15.336131\n",
       "1  0.125290  0.825161 -17.500812\n",
       "2  0.110733  0.596067 -15.812474\n",
       "3  0.536350  0.324687 -18.623693\n",
       "4  0.174577  0.282452 -17.472201"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = pd.DataFrame(stan_samples.extract())\n",
    "s.head()\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is useful as you can, for example, plot and compare the individual parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.96499999999999997"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The probability that the rate theta1 is smaller than theta2\n",
    "np.mean(s.theta1 < s.theta2) # Can also be written as sum(s.theta1 < s.theta2) / len(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x114d779e8>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       "  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Created with matplotlib (http://matplotlib.org/) -->\n",
       "<svg height=\"256pt\" version=\"1.1\" viewBox=\"0 0 386 256\" width=\"386pt\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       " <defs>\n",
       "  <style type=\"text/css\">\n",
       "*{stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:100000;}\n",
       "  </style>\n",
       " </defs>\n",
       " <g id=\"figure_1\">\n",
       "  <g id=\"patch_1\">\n",
       "   <path d=\"M 0 256.117188 \n",
       "L 386.601563 256.117188 \n",
       "L 386.601563 0 \n",
       "L 0 0 \n",
       "L 0 256.117188 \n",
       "z\n",
       "\" style=\"fill:none;\"/>\n",
       "  </g>\n",
       "  <g id=\"axes_1\">\n",
       "   <g id=\"patch_2\">\n",
       "    <path d=\"M 36.65 235.239062 \n",
       "L 371.45 235.239062 \n",
       "L 371.45 12.039062 \n",
       "L 36.65 12.039062 \n",
       "z\n",
       "\" style=\"fill:#ffffff;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_3\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 49.803084 235.239062 \n",
       "L 78.149484 235.239062 \n",
       "L 78.149484 233.676662 \n",
       "L 49.803084 233.676662 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_4\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 78.149484 235.239062 \n",
       "L 106.495885 235.239062 \n",
       "L 106.495885 227.873462 \n",
       "L 78.149484 227.873462 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_5\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 106.495885 235.239062 \n",
       "L 134.842285 235.239062 \n",
       "L 134.842285 208.455062 \n",
       "L 106.495885 208.455062 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_6\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 134.842285 235.239062 \n",
       "L 163.188685 235.239062 \n",
       "L 163.188685 153.771062 \n",
       "L 134.842285 153.771062 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_7\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 163.188685 235.239062 \n",
       "L 191.535086 235.239062 \n",
       "L 191.535086 67.169463 \n",
       "L 163.188685 67.169463 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_8\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 191.535086 235.239062 \n",
       "L 219.881486 235.239062 \n",
       "L 219.881486 27.663062 \n",
       "L 191.535086 27.663062 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_9\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 219.881486 235.239062 \n",
       "L 248.227886 235.239062 \n",
       "L 248.227886 16.503062 \n",
       "L 219.881486 16.503062 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_10\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 248.227886 235.239062 \n",
       "L 276.574286 235.239062 \n",
       "L 276.574286 107.568663 \n",
       "L 248.227886 107.568663 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_11\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 276.574286 235.239062 \n",
       "L 304.920687 235.239062 \n",
       "L 304.920687 190.375862 \n",
       "L 276.574286 190.375862 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_12\">\n",
       "    <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 304.920687 235.239062 \n",
       "L 333.267087 235.239062 \n",
       "L 333.267087 226.534262 \n",
       "L 304.920687 226.534262 \n",
       "z\n",
       "\" style=\"fill:#0000ff;stroke:#000000;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_13\">\n",
       "    <path d=\"M 36.65 12.039062 \n",
       "L 371.45 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_14\">\n",
       "    <path d=\"M 36.65 235.239062 \n",
       "L 36.65 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_15\">\n",
       "    <path d=\"M 371.45 235.239062 \n",
       "L 371.45 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_16\">\n",
       "    <path d=\"M 36.65 235.239062 \n",
       "L 371.45 235.239062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
       "   </g>\n",
       "   <g id=\"matplotlib.axis_1\">\n",
       "    <g id=\"xtick_1\">\n",
       "     <g id=\"line2d_1\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 235.239062 \n",
       "L 36.65 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_2\">\n",
       "      <defs>\n",
       "       <path d=\"M 0 0 \n",
       "L 0 -4 \n",
       "\" id=\"m17367a8d9f\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_3\">\n",
       "      <defs>\n",
       "       <path d=\"M 0 0 \n",
       "L 0 4 \n",
       "\" id=\"m730a5358fc\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_1\">\n",
       "      <!-- −0.4 -->\n",
       "      <defs>\n",
       "       <path d=\"M 31.78125 66.40625 \n",
       "Q 24.171875 66.40625 20.328125 58.90625 \n",
       "Q 16.5 51.421875 16.5 36.375 \n",
       "Q 16.5 21.390625 20.328125 13.890625 \n",
       "Q 24.171875 6.390625 31.78125 6.390625 \n",
       "Q 39.453125 6.390625 43.28125 13.890625 \n",
       "Q 47.125 21.390625 47.125 36.375 \n",
       "Q 47.125 51.421875 43.28125 58.90625 \n",
       "Q 39.453125 66.40625 31.78125 66.40625 \n",
       "M 31.78125 74.21875 \n",
       "Q 44.046875 74.21875 50.515625 64.515625 \n",
       "Q 56.984375 54.828125 56.984375 36.375 \n",
       "Q 56.984375 17.96875 50.515625 8.265625 \n",
       "Q 44.046875 -1.421875 31.78125 -1.421875 \n",
       "Q 19.53125 -1.421875 13.0625 8.265625 \n",
       "Q 6.59375 17.96875 6.59375 36.375 \n",
       "Q 6.59375 54.828125 13.0625 64.515625 \n",
       "Q 19.53125 74.21875 31.78125 74.21875 \n",
       "\" id=\"BitstreamVeraSans-Roman-30\"/>\n",
       "       <path d=\"M 37.796875 64.3125 \n",
       "L 12.890625 25.390625 \n",
       "L 37.796875 25.390625 \n",
       "z\n",
       "M 35.203125 72.90625 \n",
       "L 47.609375 72.90625 \n",
       "L 47.609375 25.390625 \n",
       "L 58.015625 25.390625 \n",
       "L 58.015625 17.1875 \n",
       "L 47.609375 17.1875 \n",
       "L 47.609375 0 \n",
       "L 37.796875 0 \n",
       "L 37.796875 17.1875 \n",
       "L 4.890625 17.1875 \n",
       "L 4.890625 26.703125 \n",
       "z\n",
       "\" id=\"BitstreamVeraSans-Roman-34\"/>\n",
       "       <path d=\"M 10.59375 35.5 \n",
       "L 73.1875 35.5 \n",
       "L 73.1875 27.203125 \n",
       "L 10.59375 27.203125 \n",
       "z\n",
       "\" id=\"BitstreamVeraSans-Roman-2212\"/>\n",
       "       <path d=\"M 10.6875 12.40625 \n",
       "L 21 12.40625 \n",
       "L 21 0 \n",
       "L 10.6875 0 \n",
       "z\n",
       "\" id=\"BitstreamVeraSans-Roman-2e\"/>\n",
       "      </defs>\n",
       "      <g transform=\"translate(24.50859375 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-2212\"/>\n",
       "       <use x=\"83.7890625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"147.412109375\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"179.19921875\" xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_2\">\n",
       "     <g id=\"line2d_4\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 84.478571 235.239062 \n",
       "L 84.478571 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_5\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"84.4785714286\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_6\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"84.4785714286\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_2\">\n",
       "      <!-- −0.2 -->\n",
       "      <defs>\n",
       "       <path d=\"M 19.1875 8.296875 \n",
       "L 53.609375 8.296875 \n",
       "L 53.609375 0 \n",
       "L 7.328125 0 \n",
       "L 7.328125 8.296875 \n",
       "Q 12.9375 14.109375 22.625 23.890625 \n",
       "Q 32.328125 33.6875 34.8125 36.53125 \n",
       "Q 39.546875 41.84375 41.421875 45.53125 \n",
       "Q 43.3125 49.21875 43.3125 52.78125 \n",
       "Q 43.3125 58.59375 39.234375 62.25 \n",
       "Q 35.15625 65.921875 28.609375 65.921875 \n",
       "Q 23.96875 65.921875 18.8125 64.3125 \n",
       "Q 13.671875 62.703125 7.8125 59.421875 \n",
       "L 7.8125 69.390625 \n",
       "Q 13.765625 71.78125 18.9375 73 \n",
       "Q 24.125 74.21875 28.421875 74.21875 \n",
       "Q 39.75 74.21875 46.484375 68.546875 \n",
       "Q 53.21875 62.890625 53.21875 53.421875 \n",
       "Q 53.21875 48.921875 51.53125 44.890625 \n",
       "Q 49.859375 40.875 45.40625 35.40625 \n",
       "Q 44.1875 33.984375 37.640625 27.21875 \n",
       "Q 31.109375 20.453125 19.1875 8.296875 \n",
       "\" id=\"BitstreamVeraSans-Roman-32\"/>\n",
       "      </defs>\n",
       "      <g transform=\"translate(72.3371651786 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-2212\"/>\n",
       "       <use x=\"83.7890625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"147.412109375\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"179.19921875\" xlink:href=\"#BitstreamVeraSans-Roman-32\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_3\">\n",
       "     <g id=\"line2d_7\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 132.307143 235.239062 \n",
       "L 132.307143 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_8\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"132.307142857\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_9\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"132.307142857\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_3\">\n",
       "      <!-- 0.0 -->\n",
       "      <g transform=\"translate(124.355580357 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_4\">\n",
       "     <g id=\"line2d_10\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 180.135714 235.239062 \n",
       "L 180.135714 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_11\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"180.135714286\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_12\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"180.135714286\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_4\">\n",
       "      <!-- 0.2 -->\n",
       "      <g transform=\"translate(172.184151786 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-32\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_5\">\n",
       "     <g id=\"line2d_13\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 227.964286 235.239062 \n",
       "L 227.964286 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_14\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"227.964285714\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_15\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"227.964285714\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_5\">\n",
       "      <!-- 0.4 -->\n",
       "      <g transform=\"translate(220.012723214 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_6\">\n",
       "     <g id=\"line2d_16\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 275.792857 235.239062 \n",
       "L 275.792857 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_17\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"275.792857143\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_18\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"275.792857143\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_6\">\n",
       "      <!-- 0.6 -->\n",
       "      <defs>\n",
       "       <path d=\"M 33.015625 40.375 \n",
       "Q 26.375 40.375 22.484375 35.828125 \n",
       "Q 18.609375 31.296875 18.609375 23.390625 \n",
       "Q 18.609375 15.53125 22.484375 10.953125 \n",
       "Q 26.375 6.390625 33.015625 6.390625 \n",
       "Q 39.65625 6.390625 43.53125 10.953125 \n",
       "Q 47.40625 15.53125 47.40625 23.390625 \n",
       "Q 47.40625 31.296875 43.53125 35.828125 \n",
       "Q 39.65625 40.375 33.015625 40.375 \n",
       "M 52.59375 71.296875 \n",
       "L 52.59375 62.3125 \n",
       "Q 48.875 64.0625 45.09375 64.984375 \n",
       "Q 41.3125 65.921875 37.59375 65.921875 \n",
       "Q 27.828125 65.921875 22.671875 59.328125 \n",
       "Q 17.53125 52.734375 16.796875 39.40625 \n",
       "Q 19.671875 43.65625 24.015625 45.921875 \n",
       "Q 28.375 48.1875 33.59375 48.1875 \n",
       "Q 44.578125 48.1875 50.953125 41.515625 \n",
       "Q 57.328125 34.859375 57.328125 23.390625 \n",
       "Q 57.328125 12.15625 50.6875 5.359375 \n",
       "Q 44.046875 -1.421875 33.015625 -1.421875 \n",
       "Q 20.359375 -1.421875 13.671875 8.265625 \n",
       "Q 6.984375 17.96875 6.984375 36.375 \n",
       "Q 6.984375 53.65625 15.1875 63.9375 \n",
       "Q 23.390625 74.21875 37.203125 74.21875 \n",
       "Q 40.921875 74.21875 44.703125 73.484375 \n",
       "Q 48.484375 72.75 52.59375 71.296875 \n",
       "\" id=\"BitstreamVeraSans-Roman-36\"/>\n",
       "      </defs>\n",
       "      <g transform=\"translate(267.841294643 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-36\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_7\">\n",
       "     <g id=\"line2d_19\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 323.621429 235.239062 \n",
       "L 323.621429 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_20\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"323.621428571\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_21\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"323.621428571\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_7\">\n",
       "      <!-- 0.8 -->\n",
       "      <defs>\n",
       "       <path d=\"M 31.78125 34.625 \n",
       "Q 24.75 34.625 20.71875 30.859375 \n",
       "Q 16.703125 27.09375 16.703125 20.515625 \n",
       "Q 16.703125 13.921875 20.71875 10.15625 \n",
       "Q 24.75 6.390625 31.78125 6.390625 \n",
       "Q 38.8125 6.390625 42.859375 10.171875 \n",
       "Q 46.921875 13.96875 46.921875 20.515625 \n",
       "Q 46.921875 27.09375 42.890625 30.859375 \n",
       "Q 38.875 34.625 31.78125 34.625 \n",
       "M 21.921875 38.8125 \n",
       "Q 15.578125 40.375 12.03125 44.71875 \n",
       "Q 8.5 49.078125 8.5 55.328125 \n",
       "Q 8.5 64.0625 14.71875 69.140625 \n",
       "Q 20.953125 74.21875 31.78125 74.21875 \n",
       "Q 42.671875 74.21875 48.875 69.140625 \n",
       "Q 55.078125 64.0625 55.078125 55.328125 \n",
       "Q 55.078125 49.078125 51.53125 44.71875 \n",
       "Q 48 40.375 41.703125 38.8125 \n",
       "Q 48.828125 37.15625 52.796875 32.3125 \n",
       "Q 56.78125 27.484375 56.78125 20.515625 \n",
       "Q 56.78125 9.90625 50.3125 4.234375 \n",
       "Q 43.84375 -1.421875 31.78125 -1.421875 \n",
       "Q 19.734375 -1.421875 13.25 4.234375 \n",
       "Q 6.78125 9.90625 6.78125 20.515625 \n",
       "Q 6.78125 27.484375 10.78125 32.3125 \n",
       "Q 14.796875 37.15625 21.921875 38.8125 \n",
       "M 18.3125 54.390625 \n",
       "Q 18.3125 48.734375 21.84375 45.5625 \n",
       "Q 25.390625 42.390625 31.78125 42.390625 \n",
       "Q 38.140625 42.390625 41.71875 45.5625 \n",
       "Q 45.3125 48.734375 45.3125 54.390625 \n",
       "Q 45.3125 60.0625 41.71875 63.234375 \n",
       "Q 38.140625 66.40625 31.78125 66.40625 \n",
       "Q 25.390625 66.40625 21.84375 63.234375 \n",
       "Q 18.3125 60.0625 18.3125 54.390625 \n",
       "\" id=\"BitstreamVeraSans-Roman-38\"/>\n",
       "      </defs>\n",
       "      <g transform=\"translate(315.669866071 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-38\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_8\">\n",
       "     <g id=\"line2d_22\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 371.45 235.239062 \n",
       "L 371.45 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_23\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m17367a8d9f\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_24\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m730a5358fc\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_8\">\n",
       "      <!-- 1.0 -->\n",
       "      <defs>\n",
       "       <path d=\"M 12.40625 8.296875 \n",
       "L 28.515625 8.296875 \n",
       "L 28.515625 63.921875 \n",
       "L 10.984375 60.40625 \n",
       "L 10.984375 69.390625 \n",
       "L 28.421875 72.90625 \n",
       "L 38.28125 72.90625 \n",
       "L 38.28125 8.296875 \n",
       "L 54.390625 8.296875 \n",
       "L 54.390625 0 \n",
       "L 12.40625 0 \n",
       "z\n",
       "\" id=\"BitstreamVeraSans-Roman-31\"/>\n",
       "      </defs>\n",
       "      <g transform=\"translate(363.4984375 246.8375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-31\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
       "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"matplotlib.axis_2\">\n",
       "    <g id=\"ytick_1\">\n",
       "     <g id=\"line2d_25\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 235.239062 \n",
       "L 371.45 235.239062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_26\">\n",
       "      <defs>\n",
       "       <path d=\"M 0 0 \n",
       "L 4 0 \n",
       "\" id=\"maaf83d69fe\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_27\">\n",
       "      <defs>\n",
       "       <path d=\"M 0 0 \n",
       "L -4 0 \n",
       "\" id=\"m3160b66a81\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"235.2390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_9\">\n",
       "      <!-- 0 -->\n",
       "      <g transform=\"translate(26.2875 237.9984375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_2\">\n",
       "     <g id=\"line2d_28\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 190.599063 \n",
       "L 371.45 190.599063 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_29\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"190.5990625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_30\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"190.5990625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_10\">\n",
       "      <!-- 200 -->\n",
       "      <g transform=\"translate(13.5625 193.3584375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-32\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_3\">\n",
       "     <g id=\"line2d_31\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 145.959063 \n",
       "L 371.45 145.959063 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_32\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"145.9590625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_33\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"145.9590625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_11\">\n",
       "      <!-- 400 -->\n",
       "      <g transform=\"translate(13.5625 148.7184375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_4\">\n",
       "     <g id=\"line2d_34\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 101.319063 \n",
       "L 371.45 101.319063 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_35\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"101.3190625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_36\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"101.3190625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_12\">\n",
       "      <!-- 600 -->\n",
       "      <g transform=\"translate(13.5625 104.0784375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-36\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_5\">\n",
       "     <g id=\"line2d_37\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 56.679063 \n",
       "L 371.45 56.679063 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_38\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"56.6790625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_39\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"56.6790625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_13\">\n",
       "      <!-- 800 -->\n",
       "      <g transform=\"translate(13.5625 59.4384375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-38\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_6\">\n",
       "     <g id=\"line2d_40\">\n",
       "      <path clip-path=\"url(#p8bbfefa1ad)\" d=\"M 36.65 12.039062 \n",
       "L 371.45 12.039062 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-dasharray:1.000000,3.000000;stroke-dashoffset:0.0;stroke-width:0.5;\"/>\n",
       "     </g>\n",
       "     <g id=\"line2d_41\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"36.65\" xlink:href=\"#maaf83d69fe\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"line2d_42\">\n",
       "      <g>\n",
       "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"371.45\" xlink:href=\"#m3160b66a81\" y=\"12.0390625\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_14\">\n",
       "      <!-- 1000 -->\n",
       "      <g transform=\"translate(7.2 14.7984375)scale(0.1 -0.1)\">\n",
       "       <use xlink:href=\"#BitstreamVeraSans-Roman-31\"/>\n",
       "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "  </g>\n",
       " </g>\n",
       " <defs>\n",
       "  <clipPath id=\"p8bbfefa1ad\">\n",
       "   <rect height=\"223.2\" width=\"334.8\" x=\"36.65\" y=\"12.0390625\"/>\n",
       "  </clipPath>\n",
       " </defs>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115f31c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting distribution of the difference between theta1 and theta2\n",
    "(s.theta2 - s.theta1).hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Calculate the probability that the difference between the two underlying rates is smaller than 0.2.**\n",
    "\n",
    "*Hint:* `abs(x - y)` calculates the absolute difference between x and y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Cows and disease\n",
    "\n",
    "Farmer Jöns has a huge number of cows. Earlier this year he ran an experiment where he gave 10 cows medicine A and 10 medicine B and then measured whether they got sick `0` or not `1` during the summer season. Here is the resulting data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cowA = [0, 1, 0, 0, 0, 0, 1, 0, 0, 0]\n",
    "cowB = [0, 0, 1, 1, 1, 0, 1, 1, 1, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Jöns now wants to know: How effective are the drugs? What is the evidence that medicine A is better or worse than medicine B?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Cows and milk\n",
    "Farmer Jöns has a huge number of cows. Earlier this year he ran an experiment where he gave 10 cows a special diet that he had heard could make them produce more milk. He recorded the number of liters of milk from these “diet” cows and from 15 “normal” cows during one month. This is the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "diet_milk = [651, 679, 374, 601, 401, 609, 767, 709, 704, 679]\n",
    "normal_milk = [798, 1139, 529, 609, 553, 743, 151, 544, 488, 555, 257, 692, 678, 675, 538]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Jöns now wants to know: Was the diet any good, does it results in better milk production?**\n",
    "\n",
    "**Hint 1:** To model this you might find it useful to use the Normal distribution which is called `normal` in Stan. A statement using `normal` could look like:\n",
    " \n",
    "`for(i in 1:n ) {`\n",
    "\n",
    "`  y[i] ~ normal(mu, sigma)`\n",
    "\n",
    "`}`\n",
    "\n",
    "Where `mu` is the mean and `sigma` is the standard deviation and `y` is a vector of length `n`. Since Stan is partly vectorized the above could also be written without the loop like `y ~ normal(mu, sigma)`.\n",
    "\n",
    "**Hint 2:** You will have to put priors on `mu` and `sigma` and here there are many options. A lazy but often OK shortcut is to just use `uniform` distributions that are wide enough to include all thinkable values of the parameters. If you want to be extra sloppy you can actually skip putting any priors at all in which case Stan will use uniform(-Infinity, Infinity), but it’s good style to use explicit priors.\n",
    "\n",
    "## Bonus questions \n",
    "If you have made it this far, great! Below are a couple of bonus questions. How far can you reach?\n",
    "## 5. Cows and Mutant Cows\n",
    "Farmer Jöns has a huge number of cows. Due to a recent radioactive leak in a nearby power plant he fears that some of them have become mutant cows. Jöns is interested in measuring the effectiveness of a diet on normal cows, but not on mutant cows (that might produce excessive amounts of milk, or nearly no milk at all!). The following data set contains the amount of milk for cows on a diet and cows on normal diet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diet_milk = [651, 679, 374, 601, 4000, 401, 609, 767, 3890, 704, 679]\n",
    "normal_milk = [798, 1139, 529, 609, 553, 743, 3, 151, 544, 488, 15, 257, 692, 678, 675, 538]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the data points might come from mutant cows (also known as outliers).\n",
    "\n",
    "**→ Jöns now wants to know: Was the diet any good, does it results in better milk production for non-mutant cows?**\n",
    "\n",
    "**Hint:** Basically we have an outlier problem. A conventional trick in this situation is to supplement the normal distribution for a distribution with wider tails that is more sensitive to the central values and disregards the far away values (this is a little bit like trimming away some amount of the data on the left and on the right). A good choice for such a distribution is the t-distribution which is like the normal but with a third parameter called the “degrees of freedom”. The lower the “degrees of freedom” the wider the tails and when this parameter is larger than about 50 the t-distribution is practically the same as the normal. A good choice for the problem with the mutant cows would be to use a t distribution with around 3 degrees of freedom:\n",
    "\n",
    "`y ~ student_t(3, mu, sigma);`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, you could also estimate the “degrees of freedom” as a free parameter, but that might be overkill in this case…\n",
    "\n",
    "## 6. Chickens and diet\n",
    "Farmer Jöns has a huge number of cows. He also has chickens. He tries different diets on them too with the hope that they will produce more eggs. Below is the number of eggs produced in one week by chickens on a diet and chickens eating normal chicken stuff:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "diet_eggs = [6, 4, 2, 3, 4, 3, 0, 4, 0, 6, 3]\n",
    "normal_eggs = [4, 2, 1, 1, 2, 1, 2, 1, 3, 2, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Jöns now wants to know: Was the diet any good, does it result in the chickens producing more eggs?**\n",
    "\n",
    "**Hint:** The `poisson` distribution is a discrete distribution that is often a reasonable choice when one wants to model count data (like, for example, counts of eggs). The `poisson` has one parameter `λλ` which stands for the mean count. In Stan you would use the Poisson like this:\n",
    "\n",
    "`y ~ poisson(lambda);`\n",
    "\n",
    "where y would be a single integer or an integer array of length n ( defined like int y[n];) and lambda a real number bounded at 0.0 (`real<lower=0> lambda;`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Cows and milk in a different data format\n",
    "It's often common to have all data in a data frame. Inspect the following data frame d:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = pd.DataFrame({'milk' : [651, 679, 374, 601, 401, 609, 767, 709, 704, 679, 798, 1139,\n",
    "                            529, 609, 553, 743, 151, 544, 488, 555, 257, 692, 678, 675, 538],\n",
    "                 'group' : [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, \n",
    "                            2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at `d` you should see that it contains the same data as in exercise (4) but coded with one cow per row (The mutant cows were perhaps just a dream…). The diet group is coded as a 1 and the normal group is coded as a 2. This data could be read into Stan by using the following data list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_list = dict(y = d.milk, x = d.group, n = len(d.milk), \n",
    "                  n_groups = max(d.group))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Modify the model from (4) to work with this data format instead.**\n",
    "\n",
    "**Hint:** In your Stan code you can loop over the group variable and use it to pick out the parameters belonging to that group like this:\n",
    "\n",
    "`for(i in 1:n) {`\n",
    "\n",
    "`  y[i] ~ normal( mu[x[i]], sigma[x[i]] ) `\n",
    "\n",
    "`}`\n",
    "\n",
    "Where `mu` and `sigma` now are 2-length vectors. This is also known as indexception: You use an index (`i`) to pick out an index (`x[i]`) to pick out a value (`mu[x[i]]`). As indexing is vectorised in Stan this can actually be shortened to just:\n",
    "\n",
    "`y ~ normal( mu[x], sigma[x] );`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Cows and more diets\n",
    "Farmer Jöns has a huge number of cows. He also has a huge number of different diets he wants to try. In addition to the diet he already tried, he tries another diet (let’s call it diet 2) on 10 more cows. Copy-n-paste the following into R and inspect the resulting data frame `d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "d = pd.DataFrame({'milk' : [651, 679, 374, 601, 401, 609, 767, 709, 704, 679, 798, 1139, 529,\n",
    "                            609, 553, 743, 151, 544, 488, 555, 257, 692, 678, 675, 538, 1061,\n",
    "                            721, 595, 784, 877, 562, 800, 684, 741, 516],\n",
    "                  'group' : [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
    "                             2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It contains the same data as in the last exercise but with 10 added rows for diet 2 which is coded as group = 3.\n",
    "\n",
    "**→ Now Jöns now wants to know: Which diet seems best, if any? How much more milk should he be expecting to produce using the best diet compared to the others?**\n",
    "\n",
    "**Hint:** If you looped or used vectorization in a smart way you should be able to use the same model as in Question 7.\n",
    "\n",
    "## 9. Cows and sunshine\n",
    "Farmer Jöns has a huge number of cows. He is wondering whether the amount of time a cow spends outside in the sunshine affects how much milk she produces. To test this he makes a controlled experiment where he picks out 20 cows and assigns each a number of hours she should spend outside each day. The experiment runs for a month and Jöns records the number of liters of milk each cow produces. Copy-n-paste the following into R and inspect the resulting data frame `d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "d = pd.DataFrame({'milk' : [685, 691, 476, 1151, 879, 725, 1190, 1107, 809, 539,\n",
    "                            298, 805, 820, 498, 1026, 1217, 1177, 684, 1061, 834],\n",
    "                 'hours' : [3, 7, 6, 10, 6, 5, 10, 11, 9, 3, 6, 6, 3, 5, 8, 11, \n",
    "                            12, 9, 5, 5]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**→ Using this data on hours of sunshine and resulting liters of milk Jöns wants to know: Does sunshine affect milk production positively or negatively?**\n",
    "\n",
    "**Hint 1:** A model probing the question above requires quite small changes to the model you developed in Question 8.\n",
    "\n",
    "**Hint 2:** You do remember the equation for the (linear regression) line? If not, here it is: `mu <- beta0 + beta1 * x;`\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}