{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SMIPP 21/22 - Exercise Sheet 5\n",
    "\n",
    "## Prof. Dr. K. Reygers, Dr. R. Stamen, Dr. M. Völkl\n",
    "\n",
    "## Hand in by: Thursday, November 25th: 12:00\n",
    "### Submit the file(s) through the Übungsgruppenverwaltung\n",
    "\n",
    "\n",
    "### Names (up to two):\n",
    "### Points:      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1 Estimator of the variance (15 points)\n",
    "We want to calculate the maximum-likelihood estimate for the variance $\\sigma^2$ of a Gaussian centered at 0 and test if it is biased. Consider first the case that we know that the true mean is zero: $\\mu = 0$. (This is quite unrealistic. We will discuss further down the realistic case.)\n",
    "\n",
    "a) First of all we need to calculate the log-likelihood of a Gaussian function for $n$ measurements: $ \\ln( L(\\sigma^2))$.\n",
    "\n",
    "b) Then we need to find the value $\\sigma_0^2$ that maximizes the likelihood. It suffices to put the first derivative with respect to $\\sigma^2$ to 0:\n",
    "\n",
    "c) Now calculate the expectation value of this estimator $\\sigma_0^2$: $ E[\\sigma_0^2]$. Is the estimator biased or unbiased?\n",
    "\n",
    "d) Now assume that we don't know the true mean value $\\mu$. We replace $\\mu$ by the mean value of the measurements $\\bar{x} = \\frac{1}{n}\\sum x_i$. Write down the new estimator and calculate its expectation value. Is the new estimator now biased. What happens when the number of measurements $n$ is large?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Biased and unbiased estimators (10 points)\n",
    "\n",
    "We want to estimate numerically what is the best estimator of the width $\\sigma^2$ of a gaussian distributed set of measurements. We will consider two estimators:\n",
    "   * the population variance: $ \\sigma_P^2 = \\frac{1}{n}\\sum_{i=1}^{n}(x_i - \\mu_0)^2 $\n",
    "   * the sample variance: $ \\sigma_S^2 = \\frac{1}{n-1}\\sum_{i=1}^{n}(x_i - \\mu_0)^2 $\n",
    "   \n",
    "To test these estimators we want to generate a large number of fake experiments. In each experiment we will have $n$ measurements following a gaussian distributed function and we want to estimate its true width $\\sigma^2$.\n",
    "\n",
    "a) Generate 1000 experiments each with $n=100$ measurements. Randomly generate the measurements from a Gaussian with $\\mu=3$ and $\\sigma=5$. Plot histograms of the values of $\\sigma_P^2$ and $\\sigma_S^2$ of obtained from the 1000 experiments. Also calculate the mean values for the two estimators.\n",
    "\n",
    "b) Now do the same, but this time generate experiments with $n=4$. Do the mean of the two estimators agree?\n",
    "\n",
    "c) Now scan repeat this procedure scanning the number of measurement in each experiment $n$ between 3 and 50. Plot the mean $\\sigma_S^2$ and the mean $\\sigma_P^2$ as a function of $n$. Do you see the trend as $n$ becomes larger and larger?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3   Unbinned Maximum Likelihood Fit  (15 points)\n",
    "\n",
    "Maximum Likelihood Fits are a standard technique in particle physics to determine theory parameters from measurements (e.g. the lifetime of elementary particles). These fits are formulated as a minimisation problem which is solved by the package Minuit which is a quasi standard in particle physics. The roots of Minuit go back to the early 1970's. The python version iminuit is based on the C++ implementation Minuit2. https://iminuit.readthedocs.io/en/stable/index.html\n",
    "\n",
    "A set of measurements is assumed to follow the normalized probability density function\n",
    "\n",
    "$$ f(t) = r \\frac{1}{\\tau_1} \\exp(-t/\\tau_1) + (1-r) \\frac{1}{\\tau_2} \\exp(-t/\\tau_2)$$\n",
    "\n",
    "where $t \\ge 0$. The variable $r$ with $0 \\le r \\le 1$ defines the relative contribution of the two components with decay times $\\tau_1$ and $\\tau_2$, respectively. The data are read from the file double_exp_data.npy below.\n",
    "\n",
    "Determine the parameters $r$, $\\tau_1$, and $\\tau_2$ from an unbinned maximum likelihood fit. In order to solve this problem you need to download the file double_exp_data.npy as well and store it in the same directory as this jupyter notebook.\n",
    "\n",
    "a) Program the probability density function $f(t; \\tau_1,\\tau_2,r)$ and a function that calculates the negative log likelihood.\n",
    "\n",
    "b) Define an iminuit object and perform the unbinned maximum likelihood fit.\n",
    "\n",
    "c) Plot the total probability density function along with individual components on top of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from iminuit import Minuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data\n",
    "x = np.loadtxt(\"https://www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/double_exp_data.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcbklEQVR4nO3deXxU9b3/8dc3KwmBhBhAwi4BFFFAU9CKiIpVq6jX/lyvP+utldvaa3v1ce3D5WrVWmtrpdW297aorUur1mutIBcXpCBYEGSXfYkBQoCQhez7fO8fM5nsYWAmOSdz3s/HI4855ztnzvnMFx7vTM58z/cYay0iIhL9YpwuQEREeoYCX0TEIxT4IiIeocAXEfEIBb6IiEfEOV1AVzIyMuyoUaOcLkNEpFdZt25dobV2YNt2Vwf+qFGjWLt2rdNliIj0KsaYfR2165SOiIhHKPBFRDzClYFvjJltjJlXWlrqdCkiIlHDlYFvrX3PWjsnNTXV6VJERKKGKwNfREQiT4EvIuIRrgx8ncMXEYk8VwZ+JM7h5xytiGBFIiK9nysDP1z7iiq55NlPyCupcroUERHXiMrAX77rKAD1jbq5i4hIE1cGfrjn8CtqGwFYubeQ6rrGSJYmItJruTLwIzUO/+G/beG9TfkRqkpEpHdzZeBHUnyccboEERFXiPrAFxERPwW+iIhHuDLwdeGViEjkuTLwI/Gl7TWTMhmd0ZfqOl8EKxMR6b1cGfiRMCStD9OzMvhkV4HTpYiIuELUBj7AtNPSKatucLoMERFXiOrAHze4H6tyipwuQ0TEFaI68EekJwMw59W1VNTqk76IeJsrAz9So3QS4/xv76NtR9hXVBmJ0kREei1XBn6kplYwRlfZiog0cWXgR1K/xDinSxARcYWoD/xl988kIyWB1TnFTpciIuKoqA/8U1ISufrsTJ5YuI36Rl2EJSLeFfWBD1Db4J8Tv7iyzuFKRESc44nA/9HsM50uQUTEcZ4I/D7xsQBMe2oJPp9ueygi3uTKwO/O2TIV9yLiVa4M/EiNw28pLTkegJ2HyyO2TxGR3sSVgd8dkgOndRZ9ccjhSkREnOGZwG+SX1rtdAkiIo7wTODHxPinWYjRdAsi4lGeCfy5N07maxMGO12GiIhjPBP4U0enM0uBLyIe5pnAFxHxOk8Gvi6+EhEv8lzgv70uj+yffOx0GSIiPc6Vgd+dV9qCJlETEW9yZeB3x5W2AOnJCRHdn4hIb+LKwO8u544c4HQJIiKO8VTgt6SboYiI13g28Mc+/L7TJYiI9ChPBX7/pHjuv3x8cL2spt7BakREepanAj82xvC9i7OcLkNExBGeCvy27v7TeqdLEBHpMZ4O/E/3FDpdgohIj/F04IuIeInnA7+4so4tB7vnil4RETfxZOA/ed1EFt4zHYBzfryYq3/9qcMViYh0vzinC3DCbeeNpKa+0ekyRER6lCc/4YuIeJErA7+7Z8sUEfEiVwZ+d82W2ZFbpo5g7KCUbj+OiIjTXBn4PWnYgCT6xMc6XYaISLfzbOAnxMZwx1dHMWFIf6dLERHpEZ4N/JgYw2PXnMkpKQmUaxI1EfEAzwZ+k8y0JHKLqnTbQxGJep4P/IyURAb1S9QNUUQk6nk+8AHqGn0s2V7A/I0HufF3q5wuR0SkW3jyStu2jlXV89DfviApPpZqXYErIlFKn/CBX9002ekSRES6nQIfuG7KUGJjDLExxulSRES6jQK/AxN/9CG1DTq1IyLRRYEf0OizVNQ2AFBR20Btg0btiEh0UeCLiHiEAl9ExCMU+CIiHqHAD/jdbefyyNUTguu19TqHLyLRRYEfcMXEU7loXEZwfdpTHztYjYhI5CnwO+GzTlcgIhJZCvwWBvbr43QJIiLdpscC3xhznTHmBWPMfGPM13rquCciNSmes4Z2/20VRUScEFLgG2P+YIwpMMZsadN+hTFmpzFmjzHmga72Ya1911p7F3AHcNNJVywiIicl1NkyXwZ+A7za1GCMiQV+C1wG5AGfG2MWALHAT9u8/lvW2oLA8n8GXiciIj0opE/41trlQHGb5qnAHmttjrW2DngTuNZa+4W19uo2PwXG72fA+9ba9Z0dyxgzxxiz1hiz9ujRoyf7viLi1VW5jh5fRCSSwjmHPxQ40GI9L9DWmXuAWcD/M8Z8p7ONrLXzrLXZ1trsgQMHhlFeeIamJbH7SIVjxxcRibRwboDS0VzCnQ5mtNY+DzwfxvF6VGJ8DDsOl/HuhoNcN6Wr32MiIr1DOJ/w84DhLdaHAfnhleMun+eW8O9/2ciWg6V8/40NTpcjIhKWcAL/c2CsMWa0MSYBuBlYEImijDGzjTHzSktLI7G7k1Je458qOSUxjpV7C1mwKap+l4mIB4U6LPMNYBUw3hiTZ4y501rbAPwb8CGwHXjLWrs1EkVZa9+z1s5JTXVuTHxmWhIAozP6OlaDiEgkhXQO31p7Syfti4BFEa3IJc4bnc6mA8dateUWVjJKvwBEpJfS1ArH0eizrNtXAsDMXyxzthgRkTCEM0qn2xhjZgOzs7KyevzY/3nVGewvrmJPgX9I5rZDZWw7VNbjdYiIRJorP+E7eQ5/2mmncEP2cG6ZOoIbs4f1+PFFRLqLKwPfDUZl9OW280Y6XYaISMQo8LuQNSilXdvOw+Ws31/iQDUiIuFxZeC7YRw+QHJC+684Lv/Vcq7/r5UOVCMiEh5XBr4bxuGLiEQbVwa+m1w0rvMJ3KzVfRBFpPdQ4B/HK9+a2ulzWQ+/z8LNmnJBRHoHBX4YGn2WXYfLnS5DRCQkCvyT9F/L9jhdgojICXFl4LtllE6TxffO4MXbs1u1/fyDnQ5VIyJyclwZ+G4bpTN2cD9GZSQ7XYaISFhcGfgiIhJ5CnwREY9Q4IdozMAUXmhzHl9EpDdR4IfIGMNlEwY7XYaIyElzZeC7bZROqKy1/HLxLkoq65wuRUSkHVcGvttG6YSqrtHHc0t2s/lgKQ2NPqfLERFpxZWBHw2yHn6f+RsPOl2GiEiQAj9MXxZVcaC4CoC1ua3nyd8buE2iiIgbKPBP0oh0/4VY723K57aXVgPwzy+udrIkEZEuKfBP0OD+iQDc9JXhwbbC8lpKq+qdKklEJCQK/BP05pzz27U1Wsv2w2UOVCMiEjpXBr6bh2UOG5AEQHrfhGBb3w5uhSgi4jauDHw3D8uMj41h9UOXcsvUEcTGGAAuGt/5XbFERNzClYHvdoP79/E/9vOfz39n/UEOlVY7WZKIyHEp8CMk/1iN0yWIiHRJgR8hz3yoG6KIiLsp8EVEPEKBH4bU5IQO2/cErrD9zdI9WGvZsL+E9zbl92RpIiLtKPDD8OI3s7kxe1i79tzCSgB8Fhp9lscWbOWeNzb0dHkiIq0o8MMwNC2JsYP6AXDW0OYhpMa03s5nA49NCyIiDlDg96C31+c5XYKIeJgrA9/NV9qGo7qu0ekSRMTDXBn4br7StjOZaX2Cyx9vOxJcrqht4IuD/l9c1fWN/Otra09oojWfz3L3n9fpwi4RCZsrA783mnvjZB648nQA8kubL8I6XNa8XFJZx4dbj5B3rCrk/db7fCz64jC7jmhufREJjwI/QvomxgWnTm6pscUXtT6rL21FxDkK/G521fOfBpc/3l7gYCUi4nUK/B7UdrimiEhPUuD3oLoGHwBNZ3astdQ2aOSOiPQMBX4Pyivxj7RZsbsQgBdXfMnZj33kZEki4iEK/B5w/+XjW603fdLfcbic2sCyiEh3U+CLiHiEAr8HJCfEHnebytoGrnxuBb9cvItnPtzR4TavrMxt9dytL3zG5rxjEatTRKKbAr8HnDGk/3G3Ka2uZ/uhMp5bspvfLt3b4Tb/vWxvq+dW7i1iw34FvoiERoHfi8RoWKeIhCHO6QK8KKewgi8LK9kUOB2TW1hJQlz73715JVXBIZwVNQ09WaKIRCFXBr4xZjYwOysry+lSTsikYWlcMymT4elJnZ6WAZi/MZ/5G5vvgDXzF8v4+L4ZwfWm8J/+s6VkpvonZVuVU9hNVYuIV7jylE5vnC0T4LSBKTx/yxTuv/z0YFC3NfKUZAb162jOneblMQNTgstV9f4Lsww6nyMi4XFl4IuISOQp8F2ipKqu1Xpdmwuyahsag1MtF5TVUFUXmXP6+4oqI7IfEXE/BX6YLsjK4IZz29/I/ET97IPWY++fXbyz1fpba/OC98ad+tQSfvj25rCPufNwORc9s4zCitqw9yUi7qfAD9OEzP48c8OksPaRlhxP26nyiyvqOt44IP9Y+HfAavoroqFR8/SLeIECX0TEIxT4IiIeocDvJplpScHl+NjWQyrbfkHbr08cuS2+PG1o9PE/6/IAiOvk8tr1baZUOFxaww/f3oS1lkfe3cKlzy6joVEzcYpIMwV+N3nh9mxmjBsIwDkjBrDo+xcGn7t4/KDg8s+/cTb3zhrHsar6YFtlbfMInGEDkkM63vZDZby1No9Gn+W1z/ax92hlcAy/iAgo8LvNgL4JzBibAYAxhgmZzROoZbS48Gr0wL4kxh1/Nk0RkXAp8EVEPEKBLyLiEQr8HhYfe/wuzy+tCS7nhnglbF5gXP7RFhdRvbvhIDlHKwD/fPuLtx0JPldcWcfSnQUh7VtEooMCvxtdN2Uoc29svijrz9+exnM3T2633fSsDL538Zh27WcM6d/qy9yOpCT6Jzzdll8KwDvrDwafe3T+Vv7wjy8B+J+1B7jr1bXB515ckcPcxbtO4N2ISG+nwO9GGSmJXH9O87QLF2RlcGZm+xlAU5Pjuf/y09vdGSspPqbDefJbMqbpsePhm53NsunTxbUinqPAFxHxCAW+iIhHKPBFRDxCge8A23ZqzICWV9iG6rOcIj7dXcjWg/4vbfcUVLR6Pq+kiuq65itu9xdVsaegnKITnBJ5T0F5l+si4n49FvjGmDOMMb8zxrxtjPluTx3XjYakJnXYPiL9+NMojBuc0mr9rbV53PbSajbl+QP/bxsOtnp+6c6jPLFwa3B9xjNLmTV3eXCunlDsK6pk1tzlHCiuAmDH4TJmzV3O0XLNoy/Sm4QU+MaYPxhjCowxW9q0X2GM2WmM2WOMeaCrfVhrt1trvwPcCGSffMm931lDO75X7zkj0gC468LRnb722slDT/h4hceZW/946gPz5dcHJmOrb/CvN2qoj0ivEuon/JeBK1o2GGNigd8CVwITgFuMMROMMWcZYxa2+RkUeM01wKfAkoi9AxERCUlcKBtZa5cbY0a1aZ4K7LHW5gAYY94ErrXW/hS4upP9LAAWGGP+F3i9o22MMXOAOQAjRowIpTwREQlBSIHfiaHAgRbrecC0zjY2xswErgcSgUWdbWetnQfMA8jOztY5AxGRCAkn8Du6hLPTgLbWLgOWhXE8EREJQzijdPKA4S3WhwH54ZUjIiLdJZzA/xwYa4wZbYxJAG4GFkSiKGPMbGPMvNLS0kjsztViO7mFYWdiAnPmDEiOD/k1i7cd4cn/3d7p8//y8udMfuIjnly4DYA31uxnyXb/zJo/mr+F/MBMnE8t2hGcfbOlXy/ZzaYDx9q1i4i7hDos8w1gFTDeGJNnjLnTWtsA/BvwIbAdeMtau7Wr/YTKWvuetXZOamrHwxejxfzvXcCU4Wkn9Jqrzx7Cwnums+w/Lua+y8YBBCdYG57uH9+/4ocXc/aw0Ptu+6EyjlXV8+Kn/pk1H3znC54IhP8rq/axIXD/3I+3H+GTXUfbvf7Zxbt48/P9J/Q+RKTnhTpK55ZO2hfRxRew0rUzM/t3OstlZxLiYpgYGMffFPDTszL4+44CMlOTOFBczfD0ZIYPSGZzXvT/hSQiodPUCiIiHuHKwPfSOXwRkZ7iysD3yjl8EZGe5MrAFxGRyFPgi4h4hAJfRMQjXBn4+tJWRCTywplLp9tYa98D3svOzr7L6Vq6Q1rgKtm2Y/DjYtv//q1r8B13fxOG9GfdvpJWbdsPlZ1UbV/9aegzV//xH1/y1/XNN1y57y8bGZ3Rl3suHQvA/I0HeWrRdo6UNd8oZca4gbz6rald7vdH87ewcPMhnrnhbC45fXCw/ZF3t5CcGMuDV54Rco3h+s5r68geNYBvX3hajx1TpLu48hN+tJsyYgA7fnxFu2kVZowbeFL7u++ycax+6NJWbXmB6RBa2vzY13j9251OaApAfmlNyMf97MtiClvcKvGdDQdb3Unr092FrcIeYHkHV+q2te1QGaem9ml3u8bXPtvHKytzQ64vEj7Yepj3Nh/q0WOKdBcFvkP6xMe2a+toWp2maRO6EhNjOtxfW8nxsSc8d8+JOsELhzvV0V870DyXUE/q+SOKdA8FvoiIR7gy8PWlrYhI5Lky8HWlrYhI5Lky8EVEJPIU+CIiHqHAFxHxCAW+iIhHuDLwNUpHRCTyXBn4GqUjIhJ5rgx8ERGJPAW+i1nA57Mhb994AtseT0VNAwXl/nl1Sqrqgu019T5Kq+sBKA88AlTVNQYfS6vrsdZSVlNPR6y1+HyW/GPV+HyWoora4CRxpS32WVbdQGl1fau2jjQdr7K2gYbGriebq6lvpKC8hsbA8WvqGymprKOksi7Y1+U19ZRU1lEdeE8n4ni11tQ3UtsQ+f2KhMKVs2V61an9+zA0LYkLxw7khRVfcsO5w5k5bhAlVXW8vmY/MQb69Wn+Jxs3uB/jB/cLrvdPiu9y/zHGMHRAUki1FFXWMfUn/pkzX16ZS7/EOMprG/hgyyGeX7IbgJzCyuD28zfmA3C0vJZJj3/E1FHprMkt7nDff169n5LKOp5dvIv7LhvH3MW7+O7MMcw6YzDf+O+VTBmRxvZDZWw6cIzfLN0DQO7TV3Va66THP+LJ6yYyd/EuLj9zMD+9/uxOt73jj2v4LKeYr4wawOe5rWcYvf/y8Xzv4iymPLGYBp/lq2NOCaGnmn1ZWMnFv1jG+kcuI71vQofb3PT7VfRPiue1O7uexK6tSY9/xMv/8hVmjh90Qq8TaUmf8F1kUP8+/OOBS5gxbiC5T1/FrdNG8INZY3nsmjPZ9eSV7PjxlSQnNAf+mZmpfHjvjOB6RwGV+/RVZKb2AfyTrA0bkNwqPBfeMx2AlMQ4Jg3r/DuTW6eN4Le3noMxhrTkeP552oh220wansaQwLE2HzzW6b4Ol9awv7gKgH1F/sfymvrgXwTWhjYtdNt9FlfWsf1QeZfbNYX8xgPt69sfqKUh8Em/qKKu3TZdqaxtAOjyE/ymvFLWtvlFEyp9ypdwKfBFRDxCgS8i4hGuDHyNwxcRiTxXBr7G4YuIRJ4rA19ERCJPgS8i4hEKfBERj1Dgi4h4hAJfRMQjFPgiIh6hwBcR8QhXBr4uvBIRiTxXBr4uvBIRiTxXBr6IiESeAl9ExCMU+CIiHqHAFxHxCAW+iIhHKPBFRDxCgS8i4hEKfBERj4hzugDpHgOS46moaQCgf1I8+aU1x31N38TO/zvYwOPGA8darbdkgEOB49TU+zrd12+W7gku/3V9XnC/H2w50uoYLT22YCsAVXWNvLoql7ySanILK9l7tKLdPlftLeLXf9+NtXDt5EyKKuu4e+YY3l6XR6Ovo8r9/rL2AEkJscH1nUfKO9zujTX7Kamq47sXjWFrfhkHiqu48qwhwefnLc/h0asnsDW/jI0HjlHb4OPO6aODz1fXN3LLvM/4yT9NJCEuhtdW7ePey8bxeW4xcTExnDUslXteX883vzqKmeMHBV+Xc7SSt9fl8Y1zhvLCihxuyh5BanI8R8pquPcvG3no62cwcWjnFyu+tiqXyyacyuovi6iua2Rgv0QuPWMwFbUN/OmzfYwZmEJ63wTOHTkAn88yb0UOt503kpTA/4v5Gw8CcO3koe32eWpqH95Ys5+dh8tJiIvh+5eODb6uyYrdR4mLieH8Mad0WF9dg48XP83hzumjSYyL7XCbaLc6p4icwkqshVunjYj4/hX4UerNOefT6POH7gu3Z1NZ19Dq+UXfv5C1+4qZMKR/sO3ZGyeRf6yaFbsL+dXHuwHIGpTCnoIKJg9Pw7bIyu9eNIZhA5JYtuMoa3KLefr6s8geNYBZc5d3WVd63wSKK+sA6BMfQ029j3tnjeOXH+9qtd2Zmf3Zml8WXH95ZW5w+dH5W7s8xh1/XENtg/+9r8opAmDOjNO4/+3NXb6u5XHGD+7XaeA/9/FuDpfVcNeFp/Ho/C2s33+M3KevCj7/x3/k8tDXzwg+B7QK/Ka6fv9JDhOH9uf3y3OYPSmT///SGvolxvHSHV9h6c6jLN15tNV+f798LzX1Pq6dnMlTi3YwdnA/Lh4/iEVfHGLl3iIenb+Fd+6+oNP39sj8rdQ1Wn68cFuwLffpq9iWX8bT7+8AYOQpyXxy/8WU1zbw9Ps7mDI8jWmn+QP6B29uBFoHftM+75w+mgff+SLYfunpg4Kva9L0/r54/PIO68srqeLnH+zkaxNOJWtQSqfvI5r965/WcayqHuiewNcpnSg1OqMvWYP6ATA8PZnTT+3f6vkJmf25/fxRxMSYYNuQ1CTOHZnOjdnDg23nnZYOwMB+ia1ePzw9mbtnZjHylGQAbp46Ini8rtx+/sjg8j2XjAXg8omD220X16KuE9XQxad4gMnD0467j9mThhx3G4CuDnWcMjpU7+v8L6PjaTyJ47XV0MVOzMn/kwSF8/68oKv+jwQFvoiIR7gy8DVbpohI5Lky8DVbpohI5Lky8EVEJPIU+CIiHqHAFxHxCAW+iIhHKPBFRDzCWNu9A/3DYYw5Cuw7yZdnAIURLKc3U180U180U180i7a+GGmtHdi20dWBHw5jzFprbbbTdbiB+qKZ+qKZ+qKZV/pCp3RERDxCgS8i4hHRHPjznC7ARdQXzdQXzdQXzTzRF1F7Dl9ERFqL5k/4IiLSggJfRMQjoi7wjTFXGGN2GmP2GGMecLqe7maM+YMxpsAYs6VFW7oxZrExZnfgcUCL5x4M9M1OY0zHtx7qpYwxw40xS40x240xW40xPwi0e64/jDF9jDFrjDGbAn3xeKDdc30BYIyJNcZsMMYsDKx7sh+w1kbNDxAL7AVOAxKATcAEp+vq5vc8AzgH2NKi7efAA4HlB4CfBZYnBPokERgd6KtYp99DBPtiCHBOYLkfsCvwnj3XH/hvMZwSWI4HVgPnebEvAu/vPuB1YGFg3ZP9EG2f8KcCe6y1OdbaOuBN4FqHa+pW1trlQHGb5muBVwLLrwDXtWh/01pba639EtiDv8+igrX2kLV2fWC5HNgODMWD/WH9KgKr8YEfiwf7whgzDLgKeLFFs+f6AaLvlM5Q4ECL9bxAm9cMttYeAn8IAoMC7Z7pH2PMKGAK/k+2nuyPwGmMjUABsNha69W++BXwQ6DlDXW92A9RF/gd3WZZ406beaJ/jDEpwF+Bf7fWlnW1aQdtUdMf1tpGa+1kYBgw1RgzsYvNo7IvjDFXAwXW2nWhvqSDtl7fD02iLfDzgOEt1ocB+Q7V4qQjxpghAIHHgkB71PePMSYef9j/2Vr7TqDZs/0BYK09BiwDrsB7fXEBcI0xJhf/Kd5LjDF/wnv9AERf4H8OjDXGjDbGJAA3AwscrskJC4BvBpa/Ccxv0X6zMSbRGDMaGAuscaC+bmGMMcBLwHZr7dwWT3muP4wxA40xaYHlJGAWsAOP9YW19kFr7TBr7Sj8efB3a+1teKwfgpz+1jjSP8DX8Y/O2As87HQ9PfB+3wAOAfX4P53cCZwCLAF2Bx7TW2z/cKBvdgJXOl1/hPtiOv4/vzcDGwM/X/difwBnAxsCfbEFeDTQ7rm+aPH+ZtI8SseT/aCpFUREPCLaTumIiEgnFPgiIh6hwBcR8QgFvoiIRyjwRUQ8QoEvIuIRCnwREY/4P6HpFnvNKz4/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot data\n",
    "plt.yscale('log')\n",
    "plt.hist(x, bins=1000, density=True, histtype='step');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}