{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import seaborn as sns\n",
    "from IPython.core.display import HTML\n",
    "import pymc3 as pm\n",
    "import theano\n",
    "import theano.tensor as tt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Random variables and probability distributions\n",
    "\n",
    "![](images/random-variable-1.svg)\n",
    "\n",
    "* Random != unpredictable, but only means we are uncertain of its exact value\n",
    "* We express our uncertainty of the possible values of x in an probability distribution\n",
    "    * Discrete probability distributions\n",
    "    * Continuous probability distributions\n",
    "    * Probability functions $P(x)$ map $x$ to a probability of x being equal to that value.\n",
    "    * A probability distribution spans all possible values of $x$, which means the total probability == 1. $$\\int{P(x)} dx = 1$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<h2>Discrete probability distributions</h2>\n",
       "\n",
       "<img src=images/dieprob.png width=600px />\n",
       "<ul>\n",
       "<li>Maps discrete values x to a probability <b>mass</b> Pmf(x)</li>\n",
       "<li>Examples:</li>\n",
       "    <ul>\n",
       "        <li>Dice</li>\n",
       "        <li>Head / Tails</li>\n",
       "        <li>Having a flat tire tommorow</li>\n",
       "    </ul>\n",
       "</ul>\n",
       "\n",
       "What is P(x <= 5)\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTML(\"\"\"\n",
    "<h2>Discrete probability distributions</h2>\n",
    "\n",
    "<img src=images/dieprob.png width=600px />\n",
    "<ul>\n",
    "<li>Maps discrete values x to a probability <b>mass</b> Pmf(x)</li>\n",
    "<li>Examples:</li>\n",
    "    <ul>\n",
    "        <li>Dice</li>\n",
    "        <li>Head / Tails</li>\n",
    "        <li>Having a flat tire tommorow</li>\n",
    "    </ul>\n",
    "</ul>\n",
    "\n",
    "What is P(x <= 5)\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Continuous probability distributions\n",
    "\n",
    "![](images/gauss.png)\n",
    "* Maps continuous values to a probability **density** Pdf(x)\n",
    "* Examples:\n",
    "    * Length of people\n",
    "    * Housing prices in Amsterdam\n",
    "    * Duration of waiting time at dentist\n",
    "    \n",
    "How de we get probability mass outputs?\n",
    "\n",
    "Try in python:\n",
    "\n",
    "```python\n",
    "from scipy import stats\n",
    "norm = stats.norm()\n",
    "```\n",
    "\n",
    "$\\mu = 0$\n",
    "\n",
    "$\\sigma = 1$\n",
    "\n",
    "What is $P(x > 0)$, $P(x > 0.02)$, $P(x < 0.69)$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15700 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x > 0.02) 0.49202168628309806\n",
      "P(x < 0.69) 0.7549029063256906\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.496"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(-3, 3)\n",
    "\n",
    "norm = stats.norm(0, 1)\n",
    "\n",
    "plt.plot(x, norm.pdf(x))\n",
    "plt.plot(x, norm.cdf(x))\n",
    "plt.show()\n",
    "\n",
    "\n",
    "print('P(x > 0.02)', 1 - norm.cdf(0.02))\n",
    "print('P(x < 0.69)', norm.cdf(0.69))\n",
    "\n",
    "# monte carlo\n",
    "n = 10000\n",
    "(norm.rvs(n) > 0.02).sum() / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Two dimensional discrete distribution\n",
    "\n",
    "* Two processes $x$ and $y$\n",
    "* Only discrete values for $x$ and $y$ possible\n",
    "\n",
    "![](images/2d-discrete.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Two dimensional discrete distribution\n",
    "\n",
    "Properties:\n",
    "$$P(x, y) \\ge 0$$\n",
    "$$\\sum_{x=0}^1\\sum_{y=0}^1 P(x, y) = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Two dimensional continuous distribution\n",
    "\n",
    "* Two processes $x$ and $y$\n",
    "* An infinite number of values for $x$ and $y$ possible\n",
    "\n",
    "![](images/2d-continuous.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Two dimensional continuous distribution\n",
    "\n",
    "Properties:\n",
    "$$P(x, y) \\ge 0$$\n",
    "$$\\int_{x=-\\infty}^ \\infty \\int_{y=-\\infty}^\\infty P(x, y) dx dy = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Marginal probability: discrete\n",
    "\n",
    "* Reducing a the dimensionality of a probability distribution. **ND -> (N - 1)D**\n",
    "\n",
    "![](images/marginal1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Marginal probability: discrete\n",
    "\n",
    "Answers previous sheet:\n",
    "$$P(x=1) = P(x=1,y=1) + P(x=1, y=0) = 0.4$$\n",
    "\n",
    "$$P(x) = \\sum_{i=1}^{N_y} P(x, y=y_i)$$\n",
    "\n",
    "Sum the probabilities in the margin\n",
    "\n",
    "![](images/marginal2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Marginal probability: continuous\n",
    "\n",
    "$$P(x) = \\int_{y = -\\infty}^\\infty P(x, y) dy$$\n",
    "\n",
    "Can be approximated by sampling:\n",
    "\n",
    "![](images/marginal3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conditional probability: discrete\n",
    "\n",
    "![](images/conditional_prob.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conditional probability: continuous\n",
    "\n",
    "![](images/conditional_prob_continuous.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' formula\n",
    "\n",
    "Jointly distributed events can be written as a multiplication of a marginal and conditional distribution.\n",
    "$$P(A,B) = P(B) \\cdot P(A|B)$$\n",
    "\n",
    "And can be written for both **\n",
    "$$P(A,B) = P(A) \\cdot P(B|A)$$\n",
    "\n",
    "$$ P(B) \\cdot P(A|B) = P(A) \\cdot P(B|A)$$\n",
    "\n",
    "divide by $P(B)$\n",
    "\n",
    "$$ P(A|B) = \\frac{P(B|A) \\cdot P(A)}{P(B)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise: Conditional probability, breast cancer mammogram\n",
    "\n",
    "* **C** Individual has cancer\n",
    "* **NC** Individual does NOT have cancer\n",
    "* **+** Individual has been positively tested\n",
    "* **-** Individual has been negatively tested\n",
    "\n",
    "Prior probabilities:\n",
    "\n",
    "$P(C) = 0.01$,&nbsp;&nbsp;&nbsp; $P(+|C)=0.9$,&nbsp;&nbsp;&nbsp; $P(+|NC)=0.08$\n",
    "\n",
    "Posterior probability:\n",
    "\n",
    "What is $P(C|+)$?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise: Conditional probability, breast cancer mammogram\n",
    "\n",
    "\n",
    "$$P(C|+) = \\frac{P(+|C) \\cdot P(C)}{P(+)}$$\n",
    "\n",
    "Working out marginal probability $P(+) = \\sum_x P(+, x)$ \n",
    "\n",
    "$= P(+, C) + P(+, NC)$\n",
    "\n",
    "$= P(+|C) \\cdot P(C) + P(+|NC) \\cdot P(NC)^{**}$\n",
    "\n",
    "\n",
    "$$P(C|+) = \\frac{0.9 \\cdot 0.01}{(0.9 \\cdot 0.01) + (0.08 \\cdot 0.99)} \\cdot 100\\% = 11\\%$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A Simple model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 12.,  54., 118., 178., 234., 187., 133.,  66.,  16.,   2.]),\n",
       " array([2.48806141, 3.01859083, 3.54912026, 4.07964968, 4.6101791 ,\n",
       "        5.14070852, 5.67123794, 6.20176736, 6.73229678, 7.2628262 ,\n",
       "        7.79335562]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 5316 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(2)\n",
    "# Central Limit theorem states: that the means of a subset of samples of any distribution will approximate a gaussian distribution given enough samples.\n",
    "\n",
    "# Uniform distribution between 0, 10\n",
    "uniform = stats.uniform(0, 10)\n",
    "\n",
    "# Sample 10 samples 1000 times and calculate the mean.\n",
    "means = np.array([uniform.rvs(10).mean() for _ in range(1000)])\n",
    "\n",
    "plt.figure(figsize=(14, 4))\n",
    "plt.hist(means)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise: Fit a gaussian distribution\n",
    "\n",
    "Fit a gaussian distribution on the `means` array.\n",
    "\n",
    "Hint: A gaussian has two parameters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f04a8ee9be0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 17196 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu = means.mean()\n",
    "sd = means.std()\n",
    "\n",
    "m = stats.norm(mu, sd)\n",
    "\n",
    "x = np.linspace(0, 10)\n",
    "plt.figure(figsize=(14, 4))\n",
    "plt.plot(x, m.pdf(x))\n",
    "sns.distplot(means, kde=False, norm_hist=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What did we just do?\n",
    "\n",
    "The gaussian we fitted is actually a model that explains the data we've seen.\n",
    "\n",
    "The models has two parameters $\\theta$, the mean $\\mu$ and the standard deviation $\\sigma$.\n",
    "\n",
    "We fitted the model on a frequentist way. Our model would have a really bad fit, if we only drew 5 times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cumulative distribution\n",
    "\n",
    "* A cummulative distribution shows the integral of the probability distribution.\n",
    "* $Cdf(a)$ = $P(x \\le a)$\n",
    "\n",
    "```python\n",
    "plt.plot(x, m.pdf(x))\n",
    "plt.plot(x, m.cdf(x))\n",
    "```\n",
    "![](images/cdf.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise\n",
    "\n",
    "Given our model **m**\n",
    "\n",
    "* What is $P(x \\le 4|\\theta)$\n",
    "* What is $P(x > 8|\\theta)$\n",
    "* What is $P(5 < x < 6|\\theta)$\n",
    "* What is the probability of our data, given our models parameters? $P(x|\\theta)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exersice: answers I\n",
    "\n",
    "Using the cumulative density function:\n",
    "```python\n",
    "print(f\"\"\"\n",
    "P(x <= 4) = {m.cdf(4)}\n",
    "P(x > 8) = {1 - m.cdf(8)}\n",
    "P(5 < x < 6) = {m.cdf(6) - m.cdf(5)}\n",
    "\"\"\")\n",
    "```\n",
    "```text\n",
    "P(x <= 4) = 0.1495479734480526\n",
    "P(x > 8) = 0.0003465131086946016\n",
    "P(5 < x < 6) = 0.3527724808375644\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exersice: answers II\n",
    "\n",
    "$$P(x|\\theta) = \\prod_{i=0}^{N}{P(x_i|\\theta)}$$\n",
    "\n",
    "Where $N$ is the number of data points. \n",
    "\n",
    "The following code returns `0.0` due to nummerical underflow. $P(x)^N$ leads to very small numbers.\n",
    "```python\n",
    "np.prod([m.pdf(x) for x in means])  \n",
    "```\n",
    "\n",
    "### Log probability\n",
    "A trick used is converting to log space. \n",
    "Hint: $ln(x \\cdot y) = ln(x) + ln(y)$\n",
    "\n",
    "```python\n",
    "np.sum([np.log(m.pdf(x)) for x in means])\n",
    "```\n",
    "\n",
    "```text\n",
    "-1316.6759133098722\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Likelihood\n",
    "\n",
    "A distribution is plotted by varying the data $x$ given constant $\\theta$ parameters.\n",
    "\n",
    "$$P(x|\\theta)$$\n",
    "\n",
    "We could also vary the parameters $\\theta$ given constant (observed) data $x$.\n",
    "\n",
    "This is called a likelihood function.\n",
    "\n",
    "$$\\mathcal{L}(\\theta|x)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Equivalence principle\n",
    "\"The probability of data conditional on a parameter value equals the likelihood of that parameter value given the same data\"\n",
    "\n",
    "$$P(x|\\theta) = \\mathcal{L}(\\theta|x)$$\n",
    "\n",
    "Note:\n",
    "\n",
    "The likelihood function is not a probability distribution.\n",
    "\n",
    "$$\\int_{x=-\\infty}^ \\infty P(x) dx \\neq 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' model definition\n",
    "\n",
    "\\begin{equation}\n",
    "P(\\boldsymbol{\\theta}\\ |\\ D) =  \\frac{P(D\\ |\\ \\boldsymbol{\\theta})\\cdot P(\\boldsymbol{\\theta})}{P(D)}\n",
    "\\end{equation}\n",
    "\n",
    "Where:\n",
    "* $P(\\boldsymbol{\\theta}\\ |\\ D)$ is the __posterior__: the distribution of the parameter(s) after taking into account the observed data.\n",
    "* $P(D\\ |\\ \\boldsymbol{\\theta})$ is the __likelihood__: the probability of observing data D, given model parameters $\\boldsymbol{\\theta}$.\n",
    "* $P(\\boldsymbol{\\theta})$ represents the __prior__ (initial belief) about the distibution of the model parameters $\\boldsymbol{\\theta}$.\n",
    "* $P(D)$ is the __evidence__ of observing $D$. This is a marginal probability. Since this factor is the same for all possible model parameters being considered, it does not influence the _relative_ probabilities of different model parameters.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' model definition\n",
    "\n",
    "\\begin{equation}\n",
    "P(\\boldsymbol{\\theta}\\ |\\ D) =  \\frac{P(D\\ |\\ \\boldsymbol{\\theta})\\cdot P(\\boldsymbol{\\theta})}{P(D)}\n",
    "\\end{equation}\n",
    "\n",
    "Note that:\n",
    "The prior and the likelihood are functions, where the evidence is just a constant (Needed for normalization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example\n",
    "\n",
    "I observe two coin flips [H, H]\n",
    "\n",
    "Is this a fair coin?\n",
    "\n",
    "We can model a coin with a Bernoulli distribution. \n",
    "\n",
    "$$ x \\sim Bernoulli(p) $$\n",
    "\n",
    "![](images/bernoulli.png)\n",
    "\n",
    "So we are interested in the parameter $p$ of the distribution.\n",
    "\n",
    "Let's restrict our model to 3 parameters.\n",
    "\n",
    "| $\\theta$          | $\\mathcal{L}$ | prior | unnormalized posterior | posterior            |\n",
    "|------------|----------|-------|------------------------|----------------------|\n",
    "| 0.3          | 0.09      | 0.33  | 0.33 * 0.09 = 0.03     | 0.03 / 0.326 = 0.09 |\n",
    "| 0.5          | 0.25      | 0.33  | 0.33 * 0.25 = 0.083     | 0.083 / 0.326 = 0.26 |\n",
    "| 0.8          | 0.64      | 0.33  | 0.33 * 0.64 = 0.213     | 0.213 / 0.326 = 0.65  |\n",
    "| marginal p | -        | 1     | 0.326                  | 1                    |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solution in code\n",
    "\n",
    "```python\n",
    "prior = np.ones(3) / 3\n",
    "likelihood = stats.bernoulli(np.array([.3, .5, .8])).pmf(1)**2\n",
    "unnormalized_posterior = prior * likelihood\n",
    "normalized_posterior = unnormalized_posterior / np.sum(unnormalized_posterior)\n",
    "print(normalized_posterior)\n",
    "```\n",
    "\n",
    "```text\n",
    "array([0.09183673, 0.25510204, 0.65306122])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Practicallity\n",
    "\n",
    "In the previous example we restricted the models parameters. In practical models we can have thousands of parameters, each having an infinite domain. \n",
    "\n",
    "Computing $P(D)$ is therefore often impractical / impossible. \n",
    "\n",
    "As $P(D)$ is only a normalization constant, algorithms can still approximate a posterior. This is due to posterior being proportional to the prior * likelihood.\n",
    "\n",
    "$$ P(\\theta|D) \\propto P(D|\\theta) \\cdot P(\\theta)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PyMC3\n",
    "\n",
    "* Approximates posterior probability by sampling (for now it's magic ;) )\n",
    "\n",
    "![](images/PyMC3_banner.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Same example in PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 14000/14000 [00:02<00:00, 5751.58draws/s]\n"
     ]
    }
   ],
   "source": [
    "coin_flips = [1, 1]\n",
    "with pm.Model():\n",
    "    p = pm.Uniform('p', 0, 1)\n",
    "    \n",
    "    likelihood = pm.Bernoulli('observed', p=p, observed=coin_flips)\n",
    "    \n",
    "    trace = pm.sample(draws=3000)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7f049d161630>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f049d18f5f8>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 54112 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(p < 0.3): 0.09008333333333333\n",
      "P(p < 0.5): 0.24941666666666668\n",
      "P(p < 0.8): 0.6446666666666667\n"
     ]
    }
   ],
   "source": [
    "for p in [.3, .5, .8]:\n",
    "    print(f'P(p < {p}):', (trace['p'] < p).sum() / trace['p'].shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Same example in PyMC3\n",
    "\n",
    "* Do we really have a uniform prior when we have a coin flip?\n",
    "* Update your prior with a Beta distribution (Tip: Google beta!)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f049d0b2080>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 14080 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 1)\n",
    "plt.plot(x, stats.beta(4, 4).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 14000/14000 [00:03<00:00, 4603.02draws/s]\n"
     ]
    }
   ],
   "source": [
    "coin_flips = [1, 1]\n",
    "with pm.Model():\n",
    "    p = pm.Beta('p', 4, 4)\n",
    "    \n",
    "    likelihood = pm.Bernoulli('observed', p=p, observed=coin_flips)\n",
    "    \n",
    "    trace = pm.sample(draws=3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7f048d4c3630>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f048d4e2048>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 57484 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Who will be elected?\n",
    "\n",
    "There will be a new referendum about a possible 'Nexit'.\n",
    "\n",
    "People have to vote $0$ for staying and $1$ for leaving.\n",
    "\n",
    "Maurice de Hond has done a survey and gives a probability of \n",
    "\n",
    "\\begin{align}\n",
    "f(x)= \n",
    "\\begin{cases}\n",
    "    0.8,& \\text{staying }  \\\\\n",
    "    0.2,              & \\text{leaving}\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "Meanwhile the first votes are coming in. What is the probability of a 'Nexit' occurring?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 1, 1, 1, 1, 0, 0, 1])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "votes = np.load('data/votes.npy')\n",
    "votes[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f048b80eba8>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15624 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# translate Maurice's probabilities to a prior believe\n",
    "\n",
    "x = np.linspace(0, 1)\n",
    "plt.plot(x, stats.beta(1.3, 3).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5357.62draws/s]\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of data points: 2\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 48236 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5601.20draws/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of data points: 10\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 46532 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5444.02draws/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of data points: 50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 44044 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5619.20draws/s]\n",
      "The acceptance probability does not match the target. It is 0.8897717351599277, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.8937936259844127, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of data points: -1\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 44576 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in [2, 10, 50, -1]:\n",
    "    with pm.Model():\n",
    "        p = pm.Beta('p', 1.3, 3)\n",
    "        pm.Bernoulli('likelihood', p=p, observed=votes[:n])\n",
    "\n",
    "        trace = pm.sample(10000)\n",
    "\n",
    "    print('Number of data points:', n)\n",
    "    pm.traceplot(trace)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('data/votes.npy', stats.bernoulli(0.68).rvs(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A/B/C testing\n",
    "\n",
    "* Three different colors are tested for the Xomnia website, and randomly assigned to visitors (not uniformly)\n",
    "* Should we switch colors?\n",
    "* Assume Xomnia goes bankrupt if we don't get at least 35 calls. We expect 1000 visitors next month, what is the probability that we don't make it? Which color should we use?\n",
    "\n",
    "color | number of visitors | number of calls\n",
    "--- | --- | ---\n",
    "blue | 2000 | 35\n",
    "red | 500 | 8\n",
    "green | 1000 | 18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5808.91draws/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 43852 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Simple model for one color\n",
    "\n",
    "with pm.Model():\n",
    "    \n",
    "    # prior\n",
    "    p = pm.Uniform('p', lower=0, upper=1)\n",
    "    conversions = pm.Binomial('conversions', 2000, p, observed = [35])\n",
    "    \n",
    "    samples = pm.sample(10000)\n",
    "    \n",
    "pm.traceplot(samples);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining parameters\n",
    "\n",
    "* You can easily combine multiple parameters and observations in the same model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains: 100%|██████████| 42000/42000 [00:07<00:00, 5600.55draws/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 64304 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with pm.Model():\n",
    "    \n",
    "    p = pm.Uniform('p', lower=0, upper=1, shape=3)\n",
    "    conversions = pm.Binomial('conv', [2000, 500, 1000], p, observed = [35, 8, 18])\n",
    "    \n",
    "    trace = pm.sample(10000)\n",
    "    \n",
    "pm.traceplot(trace)\n",
    "\n",
    "blue = trace['p'][:, 0]\n",
    "red = trace['p'][:, 1]\n",
    "green = trace['p'][:, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f048b8700f0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 26648 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "sns.distplot(blue, color='blue', hist=False)\n",
    "sns.distplot(red, color='red', hist=False)\n",
    "sns.distplot(green, color='green', hist=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior samples\n",
    "\n",
    "* `pymc3` does not generate a real distribution, but samples from the posterior distribution.\n",
    "* With those samples you can \n",
    "  * Understand the data, draw conclusions\n",
    "  * Generate predictions\n",
    "\n",
    "## Exercise\n",
    "1. What is the probability that the conversion rate of blue is higher than 2%?\n",
    "2. Which color should we choose?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2357"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(blue > 0.02).sum() / blue.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6184 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "calls = stats.binom(n=1000, p=blue).rvs()\n",
    "fig = plt.hist(calls, bins=calls.max() - calls.min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[16, 23, 17, ..., 19, 24, 19],\n",
       "       [13, 19, 23, ...,  8, 28, 17],\n",
       "       [17, 20, 20, ..., 21, 10, 29]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calls = stats.binom(n=1000, p=np.stack([blue, red, green])).rvs()\n",
    "calls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 40)"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 14364 characters]\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16, 8))\n",
    "plt.title('Calls histogram per color')\n",
    "plt.hist(calls.T, (calls.max(1) - calls.min(1)).min(), color='brg', ec='black')\n",
    "plt.xlabel('Calls')\n",
    "plt.xlim(0, 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00165  0.019875 0.009775]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 3 artists>"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 14032 characters]\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean = (calls > 35).mean(1)\n",
    "print(mean)\n",
    "\n",
    "plt.figure(figsize=(16, 8))\n",
    "plt.title('Probability of evading bankrupty')\n",
    "plt.bar(np.array(list('brg')), mean, color='brg', ec='black')\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
