{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Welcome](img/welcome.jpg)\n",
    "\n",
    "\n",
    "-----\n",
    "\n",
    "![Monte Carlo Inference for Researchers in the Physical Sciences](img/title.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pV0ug9gzCfmm"
   },
   "source": [
    "## Install packages needed later in the course"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, lets make sure the necessary packages are installed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fZzYI69X_itU"
   },
   "outputs": [],
   "source": [
    "import ultranest, snowline, cmdstanpy\n",
    "print(len(dir(ultranest)) + len(dir(snowline)) + len(dir(cmdstanpy)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The answer above is the first response to put in the [questionnaire](https://johannesbuchner.github.io/PracticalInferenceForResearchersInThePhysicalSciencesCourse/preliminaries.html)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wHnbft3T2DwM"
   },
   "source": [
    "# Basic probability distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Y7Cd60IHoqeM"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VWihRRlI21sc"
   },
   "source": [
    "Useful references:\n",
    "\n",
    "* scipy: https://docs.scipy.org/doc/scipy/reference/ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "R4g2B3et6y6M"
   },
   "source": [
    "You can create probability distribution objects with \n",
    "\n",
    "* scipy.stats.uniform(minimum, width)\n",
    "* scipy.stats.norm(mean, std)\n",
    "* etc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RJT226S2-HqG"
   },
   "source": [
    "## Narrow uniform distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3fQg895423kR"
   },
   "source": [
    "Plot a distribution that has uniform probability from 1.1 to 1.6, and zero probability elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "IMTyKyAs21Y8"
   },
   "outputs": [],
   "source": [
    "x = np.linspace(1.0, 2.0, 400)\n",
    "rv_uniform = scipy.stats.uniform(1.1, 0.5)\n",
    "plt.plot(x, rv_uniform.pdf(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "khsgnGw_3eTC"
   },
   "source": [
    "What is the probability density at x=1.2?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rKp8MFK43g8b"
   },
   "outputs": [],
   "source": [
    "rv_uniform.pdf(1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hs_Gxb9N3SNe"
   },
   "source": [
    "**What is the integral over the probability density?**\n",
    "\n",
    "$\\int_{-\\infty}^{\\infty} P(x)dx=?$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "47bpmELf2wkT"
   },
   "outputs": [],
   "source": [
    "# TODO by you! (use rv_uniform)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WMwKY5_5-Mlc"
   },
   "source": [
    "## Gaussian distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "S6FVThVy76t-"
   },
   "source": [
    "Create a normal distribution object (scipy.stats.norm) centered at 1.5 with standard deviation 0.02. Call it rv. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TIflkMaE8VRd"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n",
    "rv = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2r60U_HS9Smu"
   },
   "source": [
    "Plot its probability density distribution (PDF):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8swmjBm78fbL"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AOXTSEO59U8u"
   },
   "source": [
    "Plot its cumulative probability distribution (CDF):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nFoiLc8G813I"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qNPxUTng9IpY"
   },
   "source": [
    "How much probability is below 1.48? (use rv.cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nEYTBD1L9crz"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rkPnq4g9-UDO"
   },
   "source": [
    "Does this value make sense, looking at the CDF and PDF plots?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KrjSaWfq9cQ1"
   },
   "source": [
    "How much probability is above 1.52? (use rv.cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "t2Ak1QzP9oX-"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AK2jYanG8SMB"
   },
   "source": [
    "**Combine the above to find out how much probability is between 1.44 and 1.56?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Vxv4CGnq8IrK"
   },
   "outputs": [],
   "source": [
    "# TODO by you!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oXH_k4g7-3bg"
   },
   "source": [
    "Does this match what you expect from \"3 sigma\"?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PQqMQC_z2J4J"
   },
   "source": [
    "# Generating random numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yPxsKSXz2PGf"
   },
   "source": [
    "Useful references:\n",
    "\n",
    "* numpy.random: https://docs.scipy.org/doc/numpy-1.14.0/reference/routines.random.html \n",
    "* numpy User guide: https://numpy.org/doc/stable/user/\n",
    "* matplotlib for plotting: https://matplotlib.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fq87MieBpO50"
   },
   "source": [
    "Create a grid from 1 to 1.5 with 101 points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WHkuhsC7osCv"
   },
   "outputs": [],
   "source": [
    "xgrid = np.linspace(1, 1.5, 101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "yEacJ8ro6DaN"
   },
   "outputs": [],
   "source": [
    "xgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xnhXXNOyp8AA"
   },
   "source": [
    "Set the random number generator seed to 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1uCkMnUsp6cl"
   },
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "v8l9_DTapW53"
   },
   "source": [
    "For each value, generate a random number, Gaussian distributed centered at that value.\n",
    "Use standard deviation of 0.1 for the noise introduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "agQnCKOWpVrk"
   },
   "outputs": [],
   "source": [
    "measured = # TODO by you!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "glRySWG-pp1T"
   },
   "source": [
    "Plot one against the other as a scatter plot. Assign each point a error bar of 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "l6Ny25Xlpocw"
   },
   "outputs": [],
   "source": [
    "plt.errorbar(x=xgrid, y=measured, yerr=0.1, marker='o', ls=' ')\n",
    "plt.plot(xgrid, xgrid)\n",
    "plt.xlabel('x values')\n",
    "plt.ylabel('Measured values');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Do7L9GSTqD08"
   },
   "source": [
    "**Plot a histogram of the measured values.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ghbaFFhKpuf0"
   },
   "outputs": [],
   "source": [
    "plt.hist( ????????  )\n",
    "plt.xlabel('Measured values');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Yg7H3etjqdJ7"
   },
   "source": [
    "**Compute summary statistics of this distribution: mean, median, standard deviation, 1% and 99% quantile.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4hqeD5BxqHSG"
   },
   "outputs": [],
   "source": [
    "np.quantile(measured, 0.99)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4vDJmdS-7dHY"
   },
   "source": [
    "**What is the 99% quantile (99th percentile) of this distribution?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "anm8A6Ce4xX4"
   },
   "source": [
    "# Random Walk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_VfLcTxx3aGH"
   },
   "source": [
    "Start at x=0,y=0 with seed=1. \n",
    "\n",
    "Repeat the following 10000 times: \n",
    "generate a random gaussian number, add it to x.\n",
    "generate a random gaussian number, add it to y.\n",
    "\n",
    "Plot x vs y (the sequence of 10000 points in a 2d line plot)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4g-tyYzl4J8k"
   },
   "outputs": [],
   "source": [
    "x = 0\n",
    "y = 0\n",
    "xseq = []\n",
    "yseq = []\n",
    "np.random.seed(1)\n",
    "\n",
    "\n",
    "for .... : # TODO by you!\n",
    "\n",
    "  x += np.random.normal()\n",
    "  y += np.random.normal()\n",
    "\n",
    "  # TODO by you: keep track of the x and y \n",
    "  # append to xseq and yseq\n",
    "\n",
    "  \n",
    "# visualise:\n",
    "plt.plot(xseq, yseq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AGYrlp4s5HB0"
   },
   "source": [
    "Plot histograms of x and y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KQZcf1ea4kaz"
   },
   "outputs": [],
   "source": [
    "# TODO by you!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NC5KOE1A7Jkl"
   },
   "source": [
    "The random numbers added are normal distributed with mean 0 and standard deviation 1.\n",
    "\n",
    "**Is the resulting distribution the same, narrower or wider than this?**"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Init-cloze.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}