3_random_variables_and_distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modern Methods of Data Analysis 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Variables\n",
    "When performing observations it is essentially impossible to ban all sources of random fluctuations from the experiment (remember the thing I saied about EVERYTHING, even the measurement instrument being a quantum system and therefor susceptible to random behaviour), as a result all measurements will fluctuate. Some may only fluctuate subtily and others will fluctuate more pronounced.\n",
    "\n",
    "A variable that is subject to such random fluctuations is called a random variable. In day to day operations however observable quantities (like the voltage from your lab power supply) that can be controlled very tightly are given without an error margin. This is done because the effect of the error is assumed to be minimal and therefor it's effect on the result is neglected for the sake of simplicity.\n",
    "\n",
    "### Random Vectors\n",
    "Most often an experiment does not only measure a single observable quantity but a variety of them. So taking a single measurement of the experiment involves reading multiple instruments. The result of a measurements is therefore no longer a number but a vector $\\vec{X}$ that incorporate the readings of all instruments involved.\n",
    "All operations that can be formulated for a single random variable are simply applied to each element of the vector. This gives rise to a matrix formalism for random-variable-vectors that is far more common for real experiments than analyzing each member of the vector individually. So from now on I will include a matrix formulation of all the equations where applicable.\n",
    "\n",
    "---\n",
    "### The probability density function\n",
    "If two subsequent measurements $\\omega$ and $\\omega'$ are taken of the random variable $X$ with possible values in $S$, it is highly unlikely, that both have the same value. If many measurements are taken (into the thousands or even up to billions) regions ($x_0 \\lt \\omega \\lt x_1$) can emerge where more measurements appear (the average distance between two neighboring measurements is small) along with regiones where less appear ($x_2 \\lt \\omega \\lt x_3$) (the average distance between two neighboring measurements is comparatively large). The likelihood of a measurements appearing in an interval $\\text{d}x$ is called the probability density. This probability density can be characterized by the *probability density function* or PDF.\n",
    "\n",
    "The probability density is defined by:\n",
    "\n",
    "$$ p(x) = \\lim\\limits_{N \\rightarrow \\infty}{\\frac{n}{N}\\mathrm{d}x}$$ where $n$ is the amount of measurements in the interval $x_0 \\leq x \\leq x_1$.\n",
    "\n",
    "As the probability for measuring a value in $S$ is $P(S) = 1$ it follows that\n",
    "\n",
    "$$ \\int_{S} p(x)\\,\\mathrm{d}x = 1$$\n",
    "\n",
    "for a single random variable and \n",
    "\n",
    "$$ \\prod_{i=1}^{n} \\int_{S_v} p_i(x_i) \\mathrm{d}x_i = 1$$\n",
    "\n",
    "for an $n$ dimensional vector of random variables from the set of all possible vectors $S_v$.\n",
    "\n",
    "If $S$ is not a contiguous (as for example the results of a dice roll of an n sided die) then the above integral simplifies to a sum\n",
    "\n",
    "$$ \\sum_{x\\in S} p(x) = 1$$\n",
    "\n",
    "or a sum of summs in the case of a vector with discrete entries\n",
    "\n",
    "$$ \\sum_{x_1 \\in S_1} \\text{...} \\sum_{x_N \\in S_N} \\prod_{i=1}^N p_i(x_i) = 1$$\n",
    "\n",
    "---\n",
    "\n",
    "## Making the PDF visible with Histograms\n",
    "Using computers as we are, there is a fairly simple numerical way to approximate the PDF. This is done via a histogram. The idea of a histogram is that all measurements that fall into a specific range ($x_0 < x < x_1$) are simply counted without recording their exact value. then the amount of measurements that lie within that interval is visualized by the height of a rectangle with the base of the rectangle being the interval that the measurements lie in. This is not only done for one interval but for the entire range of measurement outcomes. A histogram looks like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counts per interval: [102  88  90 104 115 101 195  98 107]\n",
      "Edges of the intervals: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.8 0.9 1. ]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st\n",
    "\n",
    "# first lets produce some (1000) random numbers in the domain [0, 1)\n",
    "rand_nums = np.random.rand(1000)\n",
    "# lets define the edges of the intervals that the numbers should be grouped in\n",
    "bin_edges = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 0.9, 1.0]\n",
    "#now lets produce a histogram from the numbers\n",
    "hist, bins = np.histogram(rand_nums, bin_edges)\n",
    "print(\"Counts per interval: {}\".format(hist))\n",
    "print(\"Edges of the intervals: {}\".format(bins))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the bin from 0.6 to 0.8 is twice as wide as all the other bins and there also where about twice as many random numbers in that interval. From this we can deduce, that the PDF of the random numbers is likely a constant function. This sort of distribution is called a *uniform distribution*. Let's now plot the histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAARLUlEQVR4nO3df4xlZX3H8fenoDb1R0F3pBuWZcAsppS2i5lQGqvFYi0SA9oaCqmIlrpiodFoUlGTamxMtBVNTS12DQRsBEGRurHYSilKNKIOSpEfooAgu112R7BoqqX8+PaPe7a9LrPMnTl35u48+34lN3POc8655/twZz+cee75kapCktSWn5t0AZKk8TPcJalBhrskNchwl6QGGe6S1KD9J10AwJo1a2p6enrSZUjSqnLDDTf8oKqm5lu2V4T79PQ0s7Ozky5DklaVJPfsadmCwzJJDklybZJbk9yS5I1d+zOTXJ3ku93PA7v2JPlQkjuS3JTkeePriiRpFKOMuT8CvKWqjgSOBc5OciRwLnBNVW0ArunmAV4KbOhem4Dzx161JOkJLRjuVbW9qr7RTf8YuA04GDgZuLhb7WLg5d30ycDHauB64IAka8deuSRpjxZ1tkySaeBo4KvAQVW1vVt0H3BQN30wcO/QZlu7tt3fa1OS2SSzc3NziyxbkvRERg73JE8DrgDeVFU/Gl5WgxvULOomNVW1uapmqmpmamreL3slSUs0UrgneRKDYP94VX26a96xa7il+7mza98GHDK0+bquTZK0QkY5WybABcBtVfWBoUVbgDO66TOAzwy1v7o7a+ZY4MGh4RtJ0goY5Tz35wOnA99KcmPX9nbgvcDlSc4E7gFO6ZZdBZwI3AH8BHjtWCuWJC1owXCvqi8B2cPi4+dZv4Cze9YlSephr7hCVdLjrV23nvu23bvwig35pYMPYfvW70+6jCYY7tJe6r5t93LoWz876TJW1D3ve9mkS2iGd4WUpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVolAdkX5hkZ5Kbh9ouS3Jj97p717NVk0wn+enQso8sZ/GSpPmN8iSmi4C/BT62q6Gq/nDXdJLzgAeH1r+zqjaOq0BJ0uKN8oDs65JMz7csSYBTgN8Zb1mSpD76jrm/ANhRVd8dajssyTeTfDHJC/a0YZJNSWaTzM7NzfUsQ5I0rG+4nwZcOjS/HVhfVUcDbwYuSfKM+Tasqs1VNVNVM1NTUz3LkCQNW3K4J9kf+H3gsl1tVfVQVd3fTd8A3Akc0bdISdLi9DlyfzHw7arauqshyVSS/brpw4ENwF39SpQkLdYop0JeCnwFeG6SrUnO7Badys8OyQC8ELipOzXyU8BZVfXAOAuWJC1slLNlTttD+2vmabsCuKJ/WZKkPrxCVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0ymP2LkyyM8nNQ23vSrItyY3d68ShZW9LckeS25P83nIVLknas1GO3C8CTpin/YNVtbF7XQWQ5EgGz1b9lW6bv9v1wGxJ0spZMNyr6jpg1Idcnwx8oqoeqqrvAXcAx/SoT5K0BH3G3M9JclM3bHNg13YwcO/QOlu7NknSClpquJ8PPAfYCGwHzlvsGyTZlGQ2yezc3NwSy5AkzWdJ4V5VO6rq0ap6DPgo/z/0sg04ZGjVdV3bfO+xuapmqmpmampqKWVIkvZgSeGeZO3Q7CuAXWfSbAFOTfKUJIcBG4Cv9StRkrRY+y+0QpJLgeOANUm2Au8EjkuyESjgbuD1AFV1S5LLgVuBR4Czq+rR5SldkrQnC4Z7VZ02T/MFT7D+e4D39ClKktSPV6hKUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQguGe5MIkO5PcPNT210m+neSmJFcmOaBrn07y0yQ3dq+PLGfxkqT5jXLkfhFwwm5tVwNHVdWvAd8B3ja07M6q2ti9zhpPmZKkxVgw3KvqOuCB3do+X1WPdLPXA+uWoTZJ0hKNY8z9j4HPDc0fluSbSb6Y5AV72ijJpiSzSWbn5ubGUIYkaZde4Z7kHcAjwMe7pu3A+qo6GngzcEmSZ8y3bVVtrqqZqpqZmprqU4YkaTdLDvckrwFeBvxRVRVAVT1UVfd30zcAdwJHjKFOSdIiLCnck5wA/DlwUlX9ZKh9Ksl+3fThwAbgrnEUKkka3f4LrZDkUuA4YE2SrcA7GZwd8xTg6iQA13dnxrwQeHeSh4HHgLOq6oF531iStGwWDPeqOm2e5gv2sO4VwBV9i5Ik9eMVqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLctSqsXbeeJBN5rV23ftLdlxZtwVMhpb3Bfdvu5dC3fnYi+77nfS+byH6lPjxyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CXt81q8jsLz3CXt81q8jsIjd0lqkEfuWpS169Zz37Z7J12GpAWMFO5JLmTwMOydVXVU1/ZM4DJgGrgbOKWqfpjBc/f+BjgR+Anwmqr6xvhL1yRM6s9XbwEgLc6owzIXASfs1nYucE1VbQCu6eYBXsrgwdgbgE3A+f3LfGKT+jLEG0pJ2luNdOReVdclmd6t+WQGD84GuBj4AvDWrv1jVVXA9UkOSLK2qraPo+D5eDQpST+rz5j7QUOBfR9wUDd9MDA8KLu1a1u2cJeW1X5PYjDaKK0eY/lCtaoqSS1mmySbGAzbsH69wxvaiz36sH8ZatXpcyrkjiRrAbqfO7v2bcAhQ+ut69p+RlVtrqqZqpqZmprqUYYkaXd9wn0LcEY3fQbwmaH2V2fgWODB5Rxvl9SQbghspV8tGvVUyEsZfHm6JslW4J3Ae4HLk5wJ3AOc0q1+FYPTIO9gcCrka8dcs6RWOQQ2NqOeLXPaHhYdP8+6BZzdpyhJUj/efkCSGmS4S1KDDHdJapDhLkkNMtwlqUHe8ncV8ra7khZiuK9CLT41RtJ4Ge59eEMpSXspw70Pr6aTtJfyC1VJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQkq9QTfJc4LKhpsOBvwAOAF4HzHXtb6+qq5ZcoSRp0ZYc7lV1O7ARIMl+wDbgSgYPxP5gVb1/LBVKkhZtXMMyxwN3VtU9Y3o/SVIP4wr3U4FLh+bPSXJTkguTHDjfBkk2JZlNMjs3NzffKpKkJeod7kmeDJwEfLJrOh94DoMhm+3AefNtV1Wbq2qmqmampqb6liFJGjKOI/eXAt+oqh0AVbWjqh6tqseAjwLHjGEfkqRFGEe4n8bQkEyStUPLXgHcPIZ9SJIWodfDOpI8Ffhd4PVDzX+VZCNQwN27LZMkrYBe4V5V/wU8a7e203tVJEnqzStUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1qNdj9gCS3A38GHgUeKSqZpI8E7gMmGbwHNVTquqHffclSRrNuI7cX1RVG6tqpps/F7imqjYA13TzkqQVslzDMicDF3fTFwMvX6b9SJLmMY5wL+DzSW5IsqlrO6iqtnfT9wEH7b5Rkk1JZpPMzs3NjaEMSdIuvcfcgd+qqm1Jng1cneTbwwurqpLU7htV1WZgM8DMzMzjlkuSlq73kXtVbet+7gSuBI4BdiRZC9D93Nl3P5Kk0fUK9yRPTfL0XdPAS4CbgS3AGd1qZwCf6bMfSdLi9B2WOQi4Msmu97qkqv45ydeBy5OcCdwDnNJzP5KkRegV7lV1F/Dr87TfDxzf570lSUvnFaqS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoCWHe5JDklyb5NYktyR5Y9f+riTbktzYvU4cX7mSpFH0eczeI8Bbquob3UOyb0hydbfsg1X1/v7lSZKWYsnhXlXbge3d9I+T3AYcPK7CJElLN5Yx9yTTwNHAV7umc5LclOTCJAfuYZtNSWaTzM7NzY2jDElSp3e4J3kacAXwpqr6EXA+8BxgI4Mj+/Pm266qNlfVTFXNTE1N9S1DkjSkV7gneRKDYP94VX0aoKp2VNWjVfUY8FHgmP5lSpIWo8/ZMgEuAG6rqg8Mta8dWu0VwM1LL0+StBR9zpZ5PnA68K0kN3ZtbwdOS7IRKOBu4PW9KpQkLVqfs2W+BGSeRVctvRxJ0jh4haokNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYtW7gnOSHJ7UnuSHLucu1HkvR4yxLuSfYDPgy8FDiSwUOzj1yOfUmSHm+5jtyPAe6oqruq6n+ATwAnL9O+JEm7SVWN/02TVwInVNWfdPOnA79RVecMrbMJ2NTNPhe4vccu1wA/6LH9arOv9Rfs877CPi/OoVU1Nd+C/ZdeTz9VtRnYPI73SjJbVTPjeK/VYF/rL9jnfYV9Hp/lGpbZBhwyNL+ua5MkrYDlCvevAxuSHJbkycCpwJZl2pckaTfLMixTVY8kOQf4F2A/4MKqumU59tUZy/DOKrKv9Rfs877CPo/JsnyhKkmaLK9QlaQGGe6S1KBVE+4L3c4gyVOSXNYt/2qS6ZWvcrxG6PObk9ya5KYk1yQ5dBJ1jtOot61I8gdJKsmqP21ulD4nOaX7rG9JcslK1zhuI/xur09ybZJvdr/fJ06iznFJcmGSnUlu3sPyJPlQ99/jpiTP673TqtrrXwy+lL0TOBx4MvDvwJG7rfOnwEe66VOByyZd9wr0+UXAL3TTb9gX+tyt93TgOuB6YGbSda/A57wB+CZwYDf/7EnXvQJ93gy8oZs+Erh70nX37PMLgecBN+9h+YnA54AAxwJf7bvP1XLkPsrtDE4GLu6mPwUcnyQrWOO4Ldjnqrq2qn7SzV7P4HqC1WzU21b8JfA+4L9XsrhlMkqfXwd8uKp+CFBVO1e4xnEbpc8FPKOb/kXgP1awvrGrquuAB55glZOBj9XA9cABSdb22edqCfeDgXuH5rd2bfOuU1WPAA8Cz1qR6pbHKH0ediaD//OvZgv2uftz9ZCq+qeVLGwZjfI5HwEckeTLSa5PcsKKVbc8Runzu4BXJdkKXAX82cqUNjGL/fe+oIndfkDjk+RVwAzw25OuZTkl+TngA8BrJlzKStufwdDMcQz+Orsuya9W1X9OtKrldRpwUVWdl+Q3gX9IclRVPTbpwlaL1XLkPsrtDP5vnST7M/hT7v4VqW55jHQLhyQvBt4BnFRVD61QbctloT4/HTgK+EKSuxmMTW5Z5V+qjvI5bwW2VNXDVfU94DsMwn61GqXPZwKXA1TVV4CfZ3CDrVaN/ZYtqyXcR7mdwRbgjG76lcC/VfdNxSq1YJ+THA38PYNgX+3jsLBAn6vqwapaU1XTVTXN4HuGk6pqdjLljsUov9v/yOConSRrGAzT3LWSRY7ZKH3+PnA8QJJfZhDucyta5craAry6O2vmWODBqtre6x0n/S3yIr5tPpHBEcudwDu6tncz+McNgw//k8AdwNeAwydd8wr0+V+BHcCN3WvLpGte7j7vtu4XWOVny4z4OYfBcNStwLeAUydd8wr0+UjgywzOpLkReMmka+7Z30uB7cDDDP4SOxM4Czhr6DP+cPff41vj+L329gOS1KDVMiwjSVoEw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ16H8B7zotQ4VAM1IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h, b, p = plt.hist(bin_edges[:-1], bin_edges, weights=hist, edgecolor='black') # The result of this function is prevented from being printed by writing it to variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we sum up all the entries in all the intervals (that are called **bins** of a histogram) we should get back the total number of measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(hist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the visualisation is a bit misleading in terms of showing that the PDF is a uniform function. This is because for an interval twice as wide as the rest, we get twice as many counts producing an area that is four times the size it should actually be. To correct for this the height of a bin has to be divided by it's with, so that the factors don't compund but instead cancel each other out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD7CAYAAACG50QgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAQzElEQVR4nO3df4xmVX3H8fdHVrD+KL92Suju4tK42m5pG8kEMCbWukaBGpakSCAqK912U4vWiqlg/YNG/5HYSiW12K1Ql8byo9SWjaWlBDCkTZc6iEV+qExR2N2CjALbpsQf6Ld/PAcd1t2dH8/MM86e9yuZzLnnnOfec5jlM/c59z53UlVIkvrwvKUegCRpdAx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOzBj6Sa5K8niSe6fVfSTJl5Pck+Tvkxwxre39SSaTfCXJG6fVn9rqJpNcvPBTkSTNZDZn+p8CTt2r7hbghKr6ZeCrwPsBkqwHzgF+sb3mz5MckuQQ4OPAacB64NzWV5I0Qitm6lBVdyRZu1fdv0zb3AGc1cobgWur6jvA15JMAie1tsmqegggybWt7/0HOvbKlStr7dq1B+oiSdrLXXfd9c2qGttX24yhPwu/CVzXyqsY/BJ41q5WB7Bzr/qTZ9rx2rVrmZiYWIAhSlI/kjy8v7ahLuQm+QDwDPDpYfaz1z63JJlIMjE1NbVQu5UkMUToJ3k78CbgLfWjB/jsBtZM67a61e2v/sdU1daqGq+q8bGxfb47kSTN07xCP8mpwPuAM6rq6WlN24FzkhyW5HhgHfAfwOeBdUmOT3Iog4u924cbuiRprmZc009yDfBaYGWSXcAlDO7WOQy4JQnAjqr6naq6L8n1DC7QPgNcUFXfb/t5J3AzcAhwVVXdtwjzkSQdQH6SH608Pj5eXsiVpLlJcldVje+rzU/kSlJHDH1J6oihL0kdMfQlqSOGvpa1Y1cfR5Il+Tp29XFLPX1pzhbiMQzSknls905eetFnl+TYD1/6piU5rjQMz/QlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5IOYKme5LpYT3H1KZuSdABL9STXxXqKq2f6ktQRQ1+SOmLoS1JHDH1J6ogXcrUgjl19HI/t3rnUw5A0A0NfC+Jgu8NBOli5vCNJHZkx9JNcleTxJPdOqzsqyS1JHmzfj2z1SXJ5kskk9yQ5cdprNrX+DybZtDjTea6D7UMVkjSs2SzvfAr4M+DqaXUXA7dW1YeTXNy2LwJOA9a1r5OBK4CTkxwFXAKMAwXclWR7VT25UBPZF5ccJOm5ZjzTr6o7gCf2qt4IbGvlbcCZ0+qvroEdwBFJjgXeCNxSVU+0oL8FOHUhJiAtmUOe7ztJLTvzvZB7TFU92sqPAce08ipg+i0cu1rd/uql5ev73/OdpJadoS/kVlUxWLJZEEm2JJlIMjE1NbVQu5UkMf/Q/0ZbtqF9f7zV7wbWTOu3utXtr/7HVNXWqhqvqvGxsbF5Dk+StC/zDf3twLN34GwCbpxWf167i+cUYE9bBroZeEOSI9udPm9odZKkEZpxTT/JNcBrgZVJdjG4C+fDwPVJNgMPA2e37jcBpwOTwNPA+QBV9USSDwGfb/0+WFV7XxyWpH3yE98LZ8bQr6pz99O0YR99C7hgP/u5CrhqTqOT9OPaXUO9WYqL5nDwXTj3MQzScrNEdw0tpYMteJeSj2GQpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkd84NpBxkfQSjoQQ/8g89junf7dVkn75fKOJHXEM/3F0OkfuZD0k8/QXwxL+EcuXGaRdCAu70hSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkeGCv0k70lyX5J7k1yT5AVJjk9yZ5LJJNclObT1PaxtT7b2tQsxAUnS7M079JOsAn4PGK+qE4BDgHOAS4HLquplwJPA5vaSzcCTrf6y1k+SNELDLu+sAH4qyQrghcCjwOuAG1r7NuDMVt7YtmntG+JTySRppOYd+lW1G/hj4BEGYb8HuAt4qqqead12AataeRWws732mdb/6PkeX5I0d8Ms7xzJ4Oz9eOBngRcBpw47oCRbkkwkmZiamhp2d5KkaYZZ3nk98LWqmqqq7wGfAV4NHNGWewBWA7tbeTewBqC1Hw58a++dVtXWqhqvqvGxsbEhhidJ2tswof8IcEqSF7a1+Q3A/cDtwFmtzybgxlbe3rZp7bdVVQ1xfEnSHA2zpn8ngwuyXwC+1Pa1FbgIuDDJJIM1+yvbS64Ejm71FwIXDzFuSdI8DPWXs6rqEuCSvaofAk7aR99vA28e5niSpOH4iVxJ6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSNDhX6SI5LckOTLSR5I8qokRyW5JcmD7fuRrW+SXJ5kMsk9SU5cmClIkmZr2DP9jwH/XFU/D/wK8ABwMXBrVa0Dbm3bAKcB69rXFuCKIY8tSZqjeYd+ksOB1wBXAlTVd6vqKWAjsK112wac2cobgatrYAdwRJJj5z1ySdKcDXOmfzwwBfxVkruTfDLJi4BjqurR1ucx4JhWXgXsnPb6Xa1OkjQiw4T+CuBE4IqqeiXwf/xoKQeAqiqg5rLTJFuSTCSZmJqaGmJ4kqS9DRP6u4BdVXVn276BwS+Bbzy7bNO+P97adwNrpr1+dat7jqraWlXjVTU+NjY2xPAkSXubd+hX1WPAziSvaFUbgPuB7cCmVrcJuLGVtwPntbt4TgH2TFsGkiSNwIohX/8u4NNJDgUeAs5n8Ivk+iSbgYeBs1vfm4DTgUng6dZXkjRCQ4V+VX0RGN9H04Z99C3ggmGOJ0kajp/IlaSOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6MnToJzkkyd1JPtu2j09yZ5LJJNclObTVH9a2J1v72mGPLUmam4U403838MC07UuBy6rqZcCTwOZWvxl4stVf1vpJkkZoqNBPshr4deCTbTvA64AbWpdtwJmtvLFt09o3tP6SpBEZ9kz/T4H3AT9o20cDT1XVM217F7CqlVcBOwFa+57WX5I0IvMO/SRvAh6vqrsWcDwk2ZJkIsnE1NTUQu5akro3zJn+q4EzknwduJbBss7HgCOSrGh9VgO7W3k3sAagtR8OfGvvnVbV1qoar6rxsbGxIYYnSdrbvEO/qt5fVaurai1wDnBbVb0FuB04q3XbBNzYytvbNq39tqqq+R5fkjR3i3Gf/kXAhUkmGazZX9nqrwSObvUXAhcvwrElSQewYuYuM6uqzwGfa+WHgJP20efbwJsX4niSpPnxE7mS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkfmHfpJ1iS5Pcn9Se5L8u5Wf1SSW5I82L4f2eqT5PIkk0nuSXLiQk1CkjQ7w5zpPwO8t6rWA6cAFyRZD1wM3FpV64Bb2zbAacC69rUFuGKIY0uS5mHeoV9Vj1bVF1r5f4EHgFXARmBb67YNOLOVNwJX18AO4Igkx8575JKkOVuQNf0ka4FXAncCx1TVo63pMeCYVl4F7Jz2sl2tTpI0IkOHfpIXA38H/H5V/c/0tqoqoOa4vy1JJpJMTE1NDTs8SdI0Q4V+kuczCPxPV9VnWvU3nl22ad8fb/W7gTXTXr661T1HVW2tqvGqGh8bGxtmeJKkvQxz906AK4EHquqj05q2A5taeRNw47T689pdPKcAe6YtA0mSRmDFEK99NfA24EtJvtjq/hD4MHB9ks3Aw8DZre0m4HRgEngaOH+IY0uS5mHeoV9V/wpkP80b9tG/gAvmezxJ0vD8RK4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1JGRh36SU5N8JclkkotHfXxJ6tlIQz/JIcDHgdOA9cC5SdaPcgyS1LNRn+mfBExW1UNV9V3gWmDjiMcgSd0adeivAnZO297V6iRJI5CqGt3BkrOAU6vqt9r224CTq+qd0/psAba0zVcAXxnikCuBbw7x+uWotzn3Nl9wzr0YZs4vraqxfTWsmP945mU3sGba9upW90NVtRXYuhAHSzJRVeMLsa/lorc59zZfcM69WKw5j3p55/PAuiTHJzkUOAfYPuIxSFK3RnqmX1XPJHkncDNwCHBVVd03yjFIUs9GvbxDVd0E3DSiwy3IMtEy09uce5svOOdeLMqcR3ohV5K0tHwMgyR1ZNmH/kyPdUhyWJLrWvudSdaOfpQLaxZzvjDJ/UnuSXJrkpcuxTgX0mwf35HkN5JUkmV/p8ds5pzk7Pazvi/J34x6jAttFv+2j0tye5K727/v05dinAslyVVJHk9y737ak+Ty9t/jniQnDn3Qqlq2XwwuBv8X8HPAocB/Auv36vO7wCda+RzguqUe9wjm/GvAC1v5HT3MufV7CXAHsAMYX+pxj+DnvA64Gziybf/MUo97BHPeCryjldcDX1/qcQ8559cAJwL37qf9dOCfgACnAHcOe8zlfqY/m8c6bAS2tfINwIYkGeEYF9qMc66q26vq6ba5g8HnIZaz2T6+40PApcC3Rzm4RTKbOf828PGqehKgqh4f8RgX2mzmXMBPt/LhwH+PcHwLrqruAJ44QJeNwNU1sAM4IsmxwxxzuYf+bB7r8MM+VfUMsAc4eiSjWxxzfZTFZgZnCsvZjHNub3vXVNU/jnJgi2g2P+eXAy9P8m9JdiQ5dWSjWxyzmfMfAW9NsovBXYDvGs3QlsyCP7pm5LdsanSSvBUYB351qceymJI8D/go8PYlHsqorWCwxPNaBu/m7kjyS1X11JKOanGdC3yqqv4kyauAv05yQlX9YKkHtlws9zP9GR/rML1PkhUM3hJ+aySjWxyzmTNJXg98ADijqr4zorEtlpnm/BLgBOBzSb7OYO1z+zK/mDubn/MuYHtVfa+qvgZ8lcEvgeVqNnPeDFwPUFX/DryAwTNqDlaz+v99LpZ76M/msQ7bgU2tfBZwW7UrJMvUjHNO8krgLxgE/nJf54UZ5lxVe6pqZVWtraq1DK5jnFFVE0sz3AUxm3/b/8DgLJ8kKxks9zw0ykEusNnM+RFgA0CSX2AQ+lMjHeVobQfOa3fxnALsqapHh9nhsl7eqf081iHJB4GJqtoOXMngLeAkgwsm5yzdiIc3yzl/BHgx8LftmvUjVXXGkg16SLOc80FllnO+GXhDkvuB7wN/UFXL9l3sLOf8XuAvk7yHwUXdty/nk7gk1zD4xb2yXae4BHg+QFV9gsF1i9OBSeBp4Pyhj7mM/3tJkuZouS/vSJLmwNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakj/w9hJJl/Deb/ggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate the width of the bins\n",
    "bin_width = [left_edge - right_edge for left_edge, right_edge in zip(bin_edges[1:],bin_edges[:-1])]\n",
    "\n",
    "# correct the correlation between height and width\n",
    "corrected_hist = [count/width for count, width in zip(hist, bin_width)]\n",
    "\n",
    "#plot the corrected hist\n",
    "h, b, p = plt.hist(bin_edges[:-1], bin_edges, weights=corrected_hist, edgecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However you can now see, that instead of the number of measurements within any given bin, we now have a number close to the total number of events as the height on each bin so we need to come up with some sort of measure to give the number on the y-axis some meaning (as it currently needs a rather convoluted explanation).\n",
    "\n",
    "The solution to this problem is to plot the **probability density** instead of the number of counts in each bin. This can be done in the following way "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAOzUlEQVR4nO3df4hlZ33H8ffHXVNpG7W4I4TMTjbSTXFRi2GIKZYaMS2bUHb/0MouVWtZXbSNFJTSFEuU+JcttSBdq0sbUgUTo4gMdSWlNhJQ12YkGrMbIuP6Y2fVZo0x/4gm2377x72W62Rm7pmdM3N3nnm/4LLnx7PnfJ+9M5899zk/bqoKSdLW96xJFyBJ6oeBLkmNMNAlqREGuiQ1wkCXpEYY6JLUiLGBnuSOJI8leXiF9X+c5KEk30jypSS/3X+ZkqRxuhyh3wnsX2X9t4FXVdVLgfcBx3uoS5K0RjvHNaiq+5PsWWX9l0ZmTwLTXXa8a9eu2rNnxc1Kkpbx1a9+9UdVNbXcurGBvkZHgM91abhnzx7m5+d73r0ktS3Jd1da11ugJ3k1g0D/3VXaHAWOAszMzPS1a0kSPV3lkuRlwD8DB6vq8ZXaVdXxqpqtqtmpqWU/MUiSLtK6Az3JDPBp4I1V9c31lyRJuhhjh1yS3AXcAOxKsgi8B3g2QFV9GLgNeAHwoSQAF6pqdqMKliQtr8tVLofHrH8L8JbeKpIkXRTvFJWkRhjoktQIA12SGmGgS1IjDHRdsq6YniHJRF5XTHvjm7aevm/9l3rzw3Nnueqv/m0i+/7u+/9wIvuV1sMjdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXtC21+DRPn7YoaVtq8WmeHqFLUiMMdElqhIEuSY0w0CWpEZ4U1VhXTM/ww3NnJ12GpDEMdI01qasB/F5PaW0ccpGkRowN9CR3JHksycMrrE+SDyZZSPJQkmv7L/OXtXhDgCStV5chlzuBfwQ+usL6m4C9w9crgH8a/rlhWrwhQJLWa+wRelXdD/x4lSYHgY/WwEng+Umu6KtAaSJ2PNtPgNpy+jgpeiUwegnE4nDZD3rYtjQZ//O0J4K15WzqSdEkR5PMJ5k/f/78Zu5akprXR6CfA3aPzE8Plz1DVR2vqtmqmp2amuph15KkX+gj0OeANw2vdrkeeLKqHG6RpE02dgw9yV3ADcCuJIvAe4BnA1TVh4ETwM3AAvBT4E83qlhJ7fFO5P6MDfSqOjxmfQF/3ltF0nY2vLpmu/EEdD+89V+6lEzo6ppJajFYJ8Vb/yWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCB/OtUX4iFFJ4xjoW8QPz52d2FP4fBqetDU45CJJjfAIfa226RcQSLr0GehrNaEvIHDYQ9I4DrlIUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjOgV6kv1JHk2ykOTWZdbPJLkvyYNJHkpyc/+lSpJWMzbQk+wAjgE3AfuAw0n2LWn2N8A9VfVy4BDwob4LlSStrssR+nXAQlWdqaqngLuBg0vaFPDc4fTzgO/3V6IkqYsuD+e6Ehj9ZoVF4BVL2rwX+Pck7wB+Dbixl+okSZ31dVL0MHBnVU0DNwMfS/KMbSc5mmQ+yfz58+d72rUkCboF+jlg98j89HDZqCPAPQBV9WXgOcCupRuqquNVNVtVs1NTUxdXsSRpWV0C/QFgb5Krk1zG4KTn3JI23wNeA5DkxQwC3UNwSdpEYwO9qi4AtwD3Ao8wuJrlVJLbkxwYNnsX8NYkXwfuAt5cVbVRRUuSnqnTNxZV1QngxJJlt41MnwZe2W9pkqS18E5RSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRnQI9yf4kjyZZSHLrCm1en+R0klNJPt5vmZKkcXaOa5BkB3AM+H1gEXggyVxVnR5psxf4a+CVVfVEkhduVMGSpOV1OUK/DlioqjNV9RRwN3BwSZu3Aseq6gmAqnqs3zIlSeN0CfQrgbMj84vDZaOuAa5J8sUkJ5Ps76tASVI3Y4dc1rCdvcANwDRwf5KXVtVPRhslOQocBZiZmelp15Ik6HaEfg7YPTI/PVw2ahGYq6qnq+rbwDcZBPwvqarjVTVbVbNTU1MXW7MkaRldAv0BYG+Sq5NcBhwC5pa0+QyDo3OS7GIwBHOmxzolSWOMDfSqugDcAtwLPALcU1Wnktye5MCw2b3A40lOA/cBf1lVj29U0ZKkZ+o0hl5VJ4ATS5bdNjJdwDuHL0nSBHinqCQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUiE6BnmR/kkeTLCS5dZV2r01SSWb7K1GS1MXYQE+yAzgG3ATsAw4n2bdMu8uBvwC+0neRkqTxuhyhXwcsVNWZqnoKuBs4uEy79wHvB37WY32SpI66BPqVwNmR+cXhsv+X5Fpgd1V9tsfaJElrsO6TokmeBXwAeFeHtkeTzCeZP3/+/Hp3LUka0SXQzwG7R+anh8t+4XLgJcAXknwHuB6YW+7EaFUdr6rZqpqdmpq6+KolSc/QJdAfAPYmuTrJZcAhYO4XK6vqyaraVVV7qmoPcBI4UFXzG1KxJGlZYwO9qi4AtwD3Ao8A91TVqSS3Jzmw0QVKkrrZ2aVRVZ0ATixZdtsKbW9Yf1mSpLXyTlFJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqRGdAj3J/iSPJllIcusy69+Z5HSSh5J8PslV/ZcqSVrN2EBPsgM4BtwE7AMOJ9m3pNmDwGxVvQz4FPC3fRcqSVpdlyP064CFqjpTVU8BdwMHRxtU1X1V9dPh7Elgut8yJUnjdAn0K4GzI/OLw2UrOQJ8bj1FSZLWbmefG0vyBmAWeNUK648CRwFmZmb63LUkbXtdjtDPAbtH5qeHy35JkhuBdwMHqurny22oqo5X1WxVzU5NTV1MvZKkFXQJ9AeAvUmuTnIZcAiYG22Q5OXARxiE+WP9lylJGmdsoFfVBeAW4F7gEeCeqjqV5PYkB4bN/g74deCTSb6WZG6FzUmSNkinMfSqOgGcWLLstpHpG3uuS5K0Rt4pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1olOgJ9mf5NEkC0luXWb9ryT5xHD9V5Ls6btQSdLqxgZ6kh3AMeAmYB9wOMm+Jc2OAE9U1W8C/wC8v+9CJUmr63KEfh2wUFVnquop4G7g4JI2B4F/HU5/CnhNkvRXpiRpnC6BfiVwdmR+cbhs2TZVdQF4EnhBHwVKkrpJVa3eIHkdsL+q3jKcfyPwiqq6ZaTNw8M2i8P5bw3b/GjJto4CR4ezvwU8epF17wJ+NLZVW+zz9mCft4f19PmqqppabsXODn/5HLB7ZH56uGy5NotJdgLPAx5fuqGqOg4c71LxapLMV9Xserezldjn7cE+bw8b1ecuQy4PAHuTXJ3kMuAQMLekzRzwJ8Pp1wH/WeMO/SVJvRp7hF5VF5LcAtwL7ADuqKpTSW4H5qtqDvgX4GNJFoAfMwh9SdIm6jLkQlWdAE4sWXbbyPTPgD/qt7RVrXvYZguyz9uDfd4eNqTPY0+KSpK2Bm/9l6RGXNKBvh0fOdChz+9McjrJQ0k+n+SqSdTZp3F9Hmn32iSVZMtfEdGlz0leP3yvTyX5+GbX2LcOP9szSe5L8uDw5/vmSdTZlyR3JHlseFn3cuuT5IPDf4+Hkly77p1W1SX5YnAC9lvAi4DLgK8D+5a0+TPgw8PpQ8AnJl33JvT51cCvDqffvh36PGx3OXA/cBKYnXTdm/A+7wUeBH5jOP/CSde9CX0+Drx9OL0P+M6k615nn38PuBZ4eIX1NwOfAwJcD3xlvfu8lI/Qt+MjB8b2uaruq6qfDmdPMrgvYCvr8j4DvI/BM4J+tpnFbZAufX4rcKyqngCoqsc2uca+delzAc8dTj8P+P4m1te7qrqfwVV/KzkIfLQGTgLPT3LFevZ5KQf6dnzkQJc+jzrC4H/4rWxsn4cfRXdX1Wc3s7AN1OV9vga4JskXk5xMsn/TqtsYXfr8XuANSRYZXFX3js0pbWLW+vs+VqfLFnXpSfIGYBZ41aRr2UhJngV8AHjzhEvZbDsZDLvcwOBT2P1JXlpVP5loVRvrMHBnVf19kt9hcG/LS6rqfydd2FZxKR+hr+WRA6z2yIEtpEufSXIj8G7gQFX9fJNq2yjj+nw58BLgC0m+w2CscW6Lnxjt8j4vAnNV9XRVfRv4JoOA36q69PkIcA9AVX0ZeA6DZ560qtPv+1pcyoG+HR85MLbPSV4OfIRBmG/1cVUY0+eqerKqdlXVnqraw+C8wYGqmp9Mub3o8rP9GQZH5yTZxWAI5sxmFtmzLn3+HvAagCQvZhDo5ze1ys01B7xpeLXL9cCTVfWDdW1x0meCx5wlvpnBkcm3gHcPl93O4BcaBm/4J4EF4L+AF0265k3o838A/w18bfiam3TNG93nJW2/wBa/yqXj+xwGQ02ngW8AhyZd8yb0eR/wRQZXwHwN+INJ17zO/t4F/AB4msEnriPA24C3jbzHx4b/Ht/o4+faO0UlqRGX8pCLJGkNDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrxf5rwQprsI8yBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate the density for each bin\n",
    "relative_count = [count/sum(hist) for count in hist]\n",
    "\n",
    "# now do the correction for the bin width, just as before\n",
    "probability_density = [rel_cnt/width for rel_cnt, width in zip(relative_count, bin_width)]\n",
    "\n",
    "#plot a histogram approximation of the probability density\n",
    "h, b, p = plt.hist(bin_edges[:-1], bin_edges, weights=probability_density, edgecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have computed the probability density the relative amount of measurements in each bin is the area of the rectangle drawn over the interval of the corresponding bin. This makes it possible to quickly comprehend the shape of a PDF without much effort. (as this is used very often the `plt.hist` function offers the keyword argument `density` that when set to `true` converts the histogram to show the PDF by the process described above).\n",
    "As the histogram now shows the probability density, it follows that the total Area of the histogram is $P(S) = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999999"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([density*width for density, width in zip(probability_density, bin_width)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is in deed the case.\n",
    "\n",
    "And now just for fun lets look at a histogram of a gaussian (also called **normal**) distribution (you will encounter this type of distribution very often). As you will undoublty be able to see, the top edges of the histogram follow the refernce function very closely. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deXhU5fn/8fc9k4VN9ohA2ASsBNkkBlGhWBTBBVxQcatUv6VaUfuzdfva4lb7dWmpttrL0qpVqSIiIkIQUVFcEAiLQMIWUAgRNCyCEEhm5ty/P2bQSQhkQiY5s9yv65orc848Z/IJJHdOnvOc5xFVxRhjTOLyuB3AGGNM3bJCb4wxCc4KvTHGJDgr9MYYk+Cs0BtjTIJLcTtAZa1bt9bOnTu7HcMYY+LK0qVLd6hqRlWvxVyh79y5M3l5eW7HMMaYuCIim4/0mnXdGGNMgrNCb4wxCc4KvTHGJDgr9MYYk+Cs0BtjTIKLqNCLyHARWScihSJyz1HaXSYiKiLZYfvuDR23TkTOi0ZoY4wxkat2eKWIeIFngHOBrcASEZmpqgWV2h0H3A4sCtuXBYwBegLtgPdE5CRVDUTvSzDGGHM0kZzR5wCFqrpJVcuBKcCoKto9DDwGHAzbNwqYoqplqvolUBh6P2NiU8AHNnW3STCR3DDVHigK294KDAhvICKnAh1UdbaI3Fnp2M8rHdu+8icQkXHAOICOHTtGltyYaDi4B1a8CuvfgeJlULYHPCnQsiuc+FPody207eN2SmNqpdZ3xoqIB5gIjD3W91DVScAkgOzsbDudMnUv4IOFT8OCP0P5Pjg+C3pdBse1BV8pfJMPy16CxZOg+zA47/+gdTe3UxtzTCIp9MVAh7DtzNC+Q44DTgE+FBGAE4CZIjIygmONqX+7N8Nr18L2lfCT8+Gnd9E2ZxTbi5+s0Kx5AxjXP437yuaSWjCXX+ce5J2SE9i2dYtLwY05NpEU+iVAdxHpQrBIjwGuPvSiqu4BWh/aFpEPgd+pap6IHABeEZGJBC/GdgcWRy++MTVUtBheuRKcAFzxMmSNBGB7cRGd7p51WPPXgA90N096/sELo/KZuPDbYB9+8KTGmLhQ7cVYVfUD44G5wBpgqqrmi8hDobP2ox2bD0wFCoB3gFtsxI1xTdFiePlSaNgcxs3/ochXp4QWXOe7lxf853HHwHR4+zZwnDoOa0z0RNRHr6q5QG6lfROO0HZIpe1HgEeOMZ8x0bF9VbDIN8mAsbOhabsaHe7g4UH/z/l64Uzu4yVo1BrOub+OwhoTXXZnrEl8+0rg1asg/Ti4flaNi/yPhN9/UAanXg+fTIS856Ma05i6EnPz0RsTVQEfTL0O9u+AG+ZAs8NG99bcBRNh79eQe1dw6GX7/rV/T2PqkJ3Rm4Q28ZL2sGUhV7+2E2l/KiJS5aNGvClw6SQ47gR4/Rdw4Lu6CW9MlNgZvUlcWz7n9r4HmRb4KZ9dcBOdLjhy082PXViz927UEkY/Dy+MgNm/hdHP1S6rMXXIzuhNYvIdgDdvYvMe5QHfz+vmc3TIgcF3weppsG5O3XwOY6LACr1JTB//BXZ/yY0zD7CPRtF7X29qhS6ftKH3sPKbAFv/eQXNGvy4v22mTeVhYod13ZjEU7IePnkSeo/hw68m0Sma7x3wHXZj1e9lI2+mTeDJ313Og/7rgWPoCjKmDtkZvUksqjD7DkhrBMP+WC+fcqV25dXAz7jOO4+uYjN8mNhjhd4klnW58NXH8LM/BG+OqicT/ZdTSgP+kDK53j6nMZGyQm8SR8AP7z0ArbpD/1/U66feRVOe8l/KEO8XDPEsr9fPbUx1rNCbxLFiMuxYH5yawFv/l59eCgxjk3MCd6e8hk15ZmKJFXqTGMpLYf7/QYcBcLI7F0J9pPCk/zJ6eLYwOsvGOZjYYYXeJIal/4F92+GcB1ydQniWM5D1TnseGJIenArZmBhghd7EP99B+PQp6DwIOp3hahQHDxP9l5OV4YVV01zNYswhVuhN/Fv+cvBsfvCd1betB3OdbJZvC8CCx23eehMTrNCb+OYvD57NdxgAXQa7nQYAxcOjn5bBzkJYN9vtOMZYoTdxbuVrsKcoOOdMDC3v90aBH1p0Dt6hq7bevXFXRIVeRIaLyDoRKRSRe6p4/SYRWSUiK0TkExHJCu3vLCIHQvtXiMiz0f4CTBJThc/+Dif0hm5D3U5TQUCBM26F4jzY/KnbcUySq7bQi4gXeAYYAWQBVx0q5GFeUdVeqtoXeByYGPbaRlXtG3rcFK3gxrDxfdixDgbeElNn8z/oew00zgie1RvjokjO6HOAQlXdpKrlwBRgVHgDVd0bttkYsL9VTZ1pm9kREeGdB0bx9fcOaadeFZ0FRaIttSEM+BUUzoNvCtzNYpJaJHd1tAeKwra3AgMqNxKRW4A7gDTgZ2EvdRGR5cBe4Peq+nEVx44DxgF07GjTu5qj215cxNB7nmV4+l084buCdndefMS2rs8imX0jLPgzLPkXXPhXd7OYpBW1i7Gq+oyqdgXuBn4f2r0N6Kiq/Qj+EnhFRJpWcewkVc1W1eyMjPqbiMrErxu8czioqbwS+Fn1jd3UqCWcMhq+mGJLDhrXRFLoi4EOYduZoX1HMgW4GEBVy1R1Z+j5UmAjcNKxRTUmqFVD4VLvJ0wPnMVuDjtviA1hC5ScevOz4Cvl9rPbHNa1ZAuUmPoQSdfNEqC7iHQhWODHAFeHNxCR7qq6IbR5AbAhtD8D2KWqARE5EegObIpWeJOcbuiXSgPx8UJghNtRjixsgZJdwFLnfn4zbC8zh/wFDTu/cr1rySSFas/oVdUPjAfmAmuAqaqaLyIPicjIULPxIpIvIisIdtFcH9o/GFgZ2j8NuElVd0X9qzDJw3EY1z+NRc7JbNBMt9NE7EX/eXTxfMNgzyq3o5gkFNEUe6qaC+RW2jch7PntRzjuDeCN2gQ0poIvP6JbSw9/L4/xvvlK5jg5lGgzrvfO5SOnj9txTJKxO2NNfFn6AjtKHd5xctxOUiM+Ung1cDZDPF/Qjh1uxzFJxgq9iR/ffwNrZ/PiFz7KSHM7TY1NDZwNwGjvApeTmGRjhd7EjxWTwfEzaanP7STHZKtm8IlzClekfIgHm9XS1B8r9CY+OA4sfRE6D2L9zvgtkq8FziZTdnCmZ7XbUUwSsUJv4sOm+fDdZug/1u0ktTLP6c8ubcKV3vluRzFJxAq9iQ8rXoEGzaHHRW4nqZVyUpkeGMQwTx4t2Vv9AcZEgRV6E/sO7oG1s6DXaEhJdztNrb0WOJs0CXCJ9xO3o5gkYYXexL78GeA/CH2urr5tHNigmSxzuln3jak3VuhN7FvxCrQ+Cdqf6naSqJkaGMJJnmL6t7UfQVP37LvMxLadG6Hoc+hzVWwuLnKMcgMDKNNUruuT6nYUkwSs0JuYcWhBkfDHQ6N7EnCU9hf8LjYWE4mSvTTmPacfV52SCoH4vC/AxI+I5roxpj5sLy76YcZHAMHhhvTf8KnTltRf30un0P5EmfHxzcAgLmi8GDZ+ACed53Yck8DsjN7ErAGetWTKDqYFBrsdpU585PRhR6kTXJTEmDpkhd7ErMs8C/heG/Kuk+12lDrhI4Upq/2wLjc4hNSYOmKF3sSkdMo5z7uEOYEcDhL/Y+eP5OWVvuDQ0YKZbkcxCcwKvYlJQzwraCoHmOmc4XaUOrW4OAAtu8LK19yOYhKYFXoTk0Z5P6NEm7LQyXI7St3rfSV89TF8V+R2EpOgIir0IjJcRNaJSKGI3FPF6zeJyCoRWSEin4hIVthr94aOWyciNrTAVKsJpQz1LGdWYCABvG7HqXu9rwh+zJ/ubg6TsKot9CLiBZ4BRgBZwFXhhTzkFVXtpap9gceBiaFjswguJt4TGA78I/R+xhzRME8e6eLj7cBAt6PUj5ZdoH1/WG2rbpq6EckZfQ5QqKqbVLUcmAKMCm+gquHT8DUGNPR8FDBFVctU9UugMPR+xhzRKO9nFDkZLNPubkepP6dcBtu+gB2FbicxCSiSQt8eCO883BraV4GI3CIiGwme0d9Ww2PHiUieiOSVlJREmt0koJbs5UzPat52BgKJcRdsRHpeAoh135g6EbWLsar6jKp2Be4Gfl/DYyeparaqZmdkZEQrkolD53sXkSIOMwOJPdrmME3bQaczYNU0UK2+vTE1EEmhLwY6hG1nhvYdyRTg4mM81iS5kd7PWOdkslY7VN840fS8BHasg28L3E5iEkwkhX4J0F1EuohIGsGLqxXu7hCR8M7UC4ANoeczgTEiki4iXYDuwOLaxzaJqENTIcezLnQ2n0TdNodkXQzigdXWfWOiq9pCr6p+YDwwF1gDTFXVfBF5SERGhpqNF5F8EVkB3AFcHzo2H5gKFADvALeoaqAOvg6TAMacEpyyN9g/n4SaZECXnwZH31j3jYmiiGavVNVcILfSvglhz28/yrGPAI8ca0CTPK46JZUVTle2aBu3o7jnlMtg5nj4enlCLbRi3GV3xprYsKOQfm29yTN2/kh6XAieVBtTb6LKCr2JDWveAoIrLyW1hi2g2zmQ/yY4jttpTIKwQm9iQ/4MFhb52UYrt5O475RLYW8xbLVxCyY6rNAb9+3aBNtXMm2N3+0k9c+betjyiceddiUHfMrfbzr7h31tMzu6ndTEMVtK0LgvNBf7GwU+SMzFpI4s4KuwfOIhH3r+yhWnbWBi76dRPAmzfKJxh53RG/cVzIB2p7J5jw0pPGROYABt5DtOlQ3VNzamGlbojbt2bw4OJcwaVX3bJPKB05cyTeV8r/XTm9qzQm/cVRAcbWOFvqJ9NGKB05sR3kUINvrG1I4VeuOugregbZ/gnOymgtxADu1kF31kk9tRTJyzQm/c810RFOfZ2fwRvOf0p1y9nO9d5HYUE+es0Bv3rAnNjZd18dHbJanvacTHTm8r9KbWrNAb9xS8BW16QauubieJWXOcHDJlB/3b2o+qOXb23WPcsfdrKFpk3TbVmBfoj0+9jM5KdTuKiWNW6I07QjdJ0dO6bY5mD034zOkZLPQ2dbE5RlbojTsK3oLjs6B1Ei0AfoxynQF0a+mB7avcjmLilBV6U/++3w5bFlq3TYTeDfTH72jwDmJjjoEVelP/1rwNqI22idBumjL/ywDkz7DuG3NMIir0IjJcRNaJSKGI3FPF63eISIGIrBSR90WkU9hrARFZEXrMrHysSXxtMztWmJ1x/jO/oaAkgLTpUWG/ObJpa3ywa6MtHG6OSbWzV4qIF3gGOBfYCiwRkZmqGv4dtxzIVtVSEbkZeBy4MvTaAVXtG+XcJo5sLy76YYbG1uxhcPqveTpwCZ3uHl2hnc3QeGQz1vr550We0JDUnm7HMXEmkjP6HKBQVTepajkwBajQuaqq81W1NLT5OZAZ3ZgmUZznXYJXlNxAjttR4sq3+xU6nfnj3EDG1EAkhb49UBS2vTW070huBOaEbTcQkTwR+VxEquyUFZFxoTZ5JSUlEUQy8WqEZxEbnbas0w5uR4k/WaOgZC18u9btJCbORPVirIhcC2QDT4Tt7qSq2cDVwJMicthtkKo6SVWzVTU7IyMjmpFMDGnJXgZ6Csh1BgDWJ19jPS4C5MepI4yJUCSFvhgIP/3KDO2rQETOAe4DRqpq2aH9qloc+rgJ+BDoV4u8Jo4N8+bhFWWOddscm+NOgI4Dg6NvjKmBSAr9EqC7iHQRkTRgDFDhlEJE+gH/JFjkvw3b30JE0kPPWwNnAjZsIEmd71nEV04bCrRT9Y1N1bJGwbf5sMNWnjKRq7bQq6ofGA/MBdYAU1U1X0QeEpGRoWZPAE2A1ysNo+wB5InIF8B84NFKo3VMkmjO95zhybdum9rqcVHwo908ZWogosXBVTUXyK20b0LY83OOcNxnQK/aBDSJ4VzvUlLEsW6b2mrWHjoMCI6+GXyn22lMnLA7Y029ON+ziCIng1VqK0nVWtao4Lw3Oze6ncTECSv0ps41bwBnelaT6+Rg3TZR0CPUY2pj6k2ErNCbOnfRSamkSYA5gQFuR0kMzTtA+2wr9CZiVuhNnRudlUKxtmKF2kpSUZM1CratgF1fup3ExAEr9KZuHdzLeV1TQhdhrdsmag5N8Ww3T5kIWKE3dWv9O6SnCLnWbVM73tQKM31Ky84sKQ6w6IX7Kuxvm9nR7aQmBkU0vNKYY1bwFsV7HZandXM7SXwL+H6YAfSQ970zuSd1Cmfc/QLFBKcOsRlATVXsjN7UnbLvYcM83ljjQ+1bLeqCN5/BcO9il5OYWGc/fabubHgXAmVMK/C7nSQhbdE2rHY6c4F3kdtRTIyzQm/qTsFb0Ph4Pi0KuJ0kYeUGBnCqp5C27HQ7iolhVuhN3SjfDxvmQY+LcGyZ0zozxwlOKTHCum/MUVihN3Wj8D3wlf44DNDUiS+1LWucjoyw7htzFFboTd3InwGNWgeXvzN1anZgAKd51tOGXW5HMTHKCr2JPt8BWD8XelwIXhvBW9cOdd8M9y5xOYmJVVboTfQVvg++/dZtU082anvWOZmcb9035gis0JvoK3gLGraAzoPcTpI0cgMDOE3WcUITm2bCHM4KvYkufxmsmwMnXwDeVLfTJI1cZwAeUS452brKzOEiKvQiMlxE1olIoYjcU8Xrd4hIgYisFJH3RaRT2GvXi8iG0OP6aIY37mqb2bHi/CsiXJjVBMq/Z8TvJv2wz9S9DZrJBqc9l2fZL1dzuGp//YuIF3gGOBfYCiwRkZmV1n5dDmSraqmI3Aw8DlwpIi2B+4FsQIGloWN3R/sLMfVve3HRYfOvjE19lr2aR+FlL9Ep9O1l86/Uj1wnh/GdtsK+b6HJ8W7HMTEkkjP6HKBQVTepajkwBahwlU1V56tqaWjzcyAz9Pw8YJ6q7goV93nA8OhEN7EmDR/DPHnMc/rjs/ny6t2cwAC8HoG1s6pvbJJKJIW+PVAUtr01tO9IbgTm1ORYERknInkikldSUhJBJBOLBnlW0lRKeTtwuttRktJa7cC6HYHgPQzGhInqxVgRuZZgN80TNTlOVSeparaqZmdkZEQzkqlHF3kXslub8KnTy+0oSUqYtsYPX30C+3e4HcbEkEgKfTHQIWw7M7SvAhE5B7gPGKmqZTU51sS/BpRxjmcZ7wROs24bF00r8IEGrPvGVBBJoV8CdBeRLiKSBowBKqxfJiL9gH8SLPLfhr00FxgmIi1EpAUwLLTPJJizPStoIgd52xnodpSktmK7Ay26WPeNqaDaQq+qfmA8wQK9Bpiqqvki8pCIjAw1ewJoArwuIitEZGbo2F3AwwR/WSwBHgrtMwnmQu9CSrQpi5webkcxPS+BLxfAPrveZYIi+htbVXOB3Er7JoQ9P+coxz4PPH+sAU3sa8wBhnqW81pgCAG8bscxvUbDJxOhYAbk/NLtNCYG2J2xptaGepbRQHzMCli3TUxo0xMyesCqaW4nMTHCCr2ptYu8C9mmLcnTk9yOYg7pdRkUfQ7fFVXf1iQ8K/SmVpqyj596vmB2YIAtAB5Lel4a/Jg/3d0cJibYT6aplfO8eaRJgLet2ya2tOoK7U617hsDWKE3tXSh53O2OBl8oV3djmIq6zUatq+EHRvcTmJcZoXeHLNWDYUzPauZ5QwEbJbKmNPzUkDsrN5YoTfH7rKsFFLEYZbNbRObmraFzmfB6mmg6nYa4yIr9OaYXdMrlUKnHQXaqfrGxh2nXAY7C2HbF24nMS6yQm+Oze7NDO6UwvTAWVi3TQzxplZYCKbV4BvxBZQnxp5+2CIxbTM7up3W1BObfcocm1VTAZjpnOlyEFNBwHfYYjALeIKrB27hH6c+VWEIrC0IkzzsjN7UnCqsnMqCzX62qk0rHetmBgbSXnaSLevdjmJcYoXe1Ny2FbBjPZNX+txOYiIwz8mmVNO5xPuJ21GMS6zQm5pbORW8abxeYIU+HpTSgHec07jQ+znplLsdx7jACr2pmYA/OC67+zC+O+h2GBOp6YFBNJVShnqWuR3FuMAKvamZLz+E/d9C7yvdTmJq4DOnJ9u0JZd6P3Y7inGBFXpTMyunQoNm0H2Y20lMDTh4mBE4kyGeL2jFHrfjmHpmhd5ErmwfrHkbsi6G1AZupzE1ND0wiBRxGOn9zO0opp5FVOhFZLiIrBORQhG5p4rXB4vIMhHxi8joSq8FQssL/rDEoIlTa2eDrxT6jHE7iTkGGzSTlU4X675JQtUWehHxAs8AI4As4CoRyarUbAswFnilirc4oKp9Q4+RVbxu4sWK/0LzTtDB5raJV9MDg+jl+YqTxBYkSSaRnNHnAIWquklVy4EpwKjwBqr6laquBJw6yGhiwe7N8OVH0O9a8FiPX7x6OzAQn3rtrD7JRPIT2x4I//W/NbQvUg1EJE9EPheRi6tqICLjQm3ySkps5fqYtOIVQKDPVW4nMbWwk2Z86PThEu8neGyKoqRRH6dmnVQ1G7gaeFJEDluhQlUnqWq2qmZnZNgt9THHcYLdNl3PhuYd3E5jaml6YBBt5DvOOdHrdhRTTyIp9MVA+E93ZmhfRFS1OPRxE/Ah0K8G+Uws+PIj2FMU7LYxce9951R2aRP+p1+a21FMPYmk0C8BuotIFxFJA8YAEY2eEZEWIpIeet4aOBMoONawxiXLJ0OD5vCTC9xOYqKgnFSmBwYx6uQU2Gddpcmg2kKvqn5gPDAXWANMVdV8EXlIREYCiMhpIrIVuBz4p4jkhw7vAeSJyBfAfOBRVbVCH08O7A6One99hY2dTyBTAmeT5hX44lW3o5h6ENF89KqaC+RW2jch7PkSgl06lY/7DOhVy4zGTaumQaDMum0STKFm8ukWP2cuewnOuBXErswmMhsnZ45u+WRo0wva9nE7iYmyfy3zwc4NsGWh21FMHbNCb45s28rg3PN2Np+QXi/wQXpTWPqi21FMHbNCbw7TNrMjIsI/xw2g1Ke0OPvmw9YbFftTP+6V+oBeo6FgBhz4zu04pg7ZmrHmMNuLi+h591SuTb+FWYHTaXb7r2hWRTtbczQBnPpzyHseVr0OOb90O42pI3ZGb6p0sfdTGksZkwPnuB3F1KV2/eCE3sHuG1W305g6YoXeVOla73uscjqzUk90O4qpa/3HwjerYOsSt5OYOmKF3hzmzA5eTvYUMTlwLmB98Qmv95XBi7KLJ7mdxNQRK/TmMDdlp7JXGzIzMNDtKKY+pDeBvtdA/gz4/hu305g6YIXeVLR/B5dnBW+RP4DdCZvQvKk/jKA66donwPHxhws7Hja6qm1mR7eTmlqyUTemouUvk54i/LfMLsImvICPTnfPAqAc+CjwKLcMKeLVgU/hDysNNroq/tkZvflRwAeL/8X7m/xs0MNmtDAJ7sXAME6Q3Qzz5LkdxUSZFXrzozUzYW8xTy4qdzuJccGHTl+2OBlcn/Ku21FMlFmhNz9a+A9oeSKz1/vdTmJc4ODh5cC5DPCs5WTZ4nYcE0VW6E1Q0RIozoMBN2O3zSSvqYEhlGo6/5OSW31jEzes0Jugz/8B6c2g79VuJzEu2kMTXgsMYaTnU9qwy+04Jkqs0Bv4rggK3oJTrwuOqTZJ7fnAcLw4jE2Z63YUEyVW6A0sehZQGPArt5OYGFCkbZjj5HCN930ac8DtOCYKIir0IjJcRNaJSKGI3FPF64NFZJmI+EVkdKXXrheRDaHH9dEKbqKkdBfkvQCnXAbN7cYYEzTJfyFNpZQx3vluRzFRUG2hFxEv8AwwAsgCrhKRrErNtgBjgVcqHdsSuB8YAOQA94tIi9rHNlGz+F/g2w9n/T+3k5gYslK7ssg5mV+kvEOK/d0f9yL5L8wBClV1k6qWA1OAUeENVPUrVV0JOJWOPQ+Yp6q7VHU3MA8YHoXcJhrK9we7bU4aDm16up3GxJhJ/gvIlB1c0dNuoI93kRT69kBR2PbW0L5IRHSsiIwTkTwRySspKYnwrU2tLX0RDuyCs+5wO4mJQR84/VjrdOC+QengBNyOY2ohJv4oU9VJqpqtqtkZGRlux0kO/nJY+DR0OhM6DnA7jYlBioe/+y8hK8MbHJVl4lYkhb4Y6BC2nRnaF4naHGui7NBasCLCDac1gb3FnPfIu7YWrDmiOU4OBSUBWPAEOJV7Zk28iKTzbQnQXUS6ECzSY4BI76qZC/wp7ALsMODeGqc0UbG9uIhOd88iFT8Ppv2WlU4T1o3+I50qLS5isxWaQxw8/HFBGa9kFMDaWZA10u1I5hhUe0avqn5gPMGivQaYqqr5IvKQiIwEEJHTRGQrcDnwTxHJDx27C3iY4C+LJcBDoX3GRZd7P6KDp4SJ/suxFaRMdV7L90OrbvDR47aubJyK6HK6quYCuZX2TQh7voRgt0xVxz4PPF+LjCaK0inn1pQ3yXNO4kOnj9txTBxwFBj0O5hxE6zLhZMvcDuSqaGYuBhr6s9V3g9oK7v4i53Nm5rodTm07ArvP2wjcOKQFfok0jAFbkl5i88CWSx0bNy8qQFvCgydACVr4IspbqcxNWSFPomMz0kjQ/aEzuaNidChtWVPuYTFxQG2vHQTDVLE1paNI3bLW7LYv5P/HZTO+4F+LNWfuJ3GxJOwtWX/6snn1bRHuO/um3kuULGv3kZrxS47o08WHz1KkzT4k9/mmzfHbqHTkwWBXoxPeYvjKHU7jomQFfpkULIeljzHpKU+Nmqks1cYU7XH/FfRQvZxS4rdLRsvrNAng3kTILURD3xY5nYSkwDytTOv+wdzgzeXzrLN7TgmAlboE93G+bB+Dgy6g5JSu9nFRMfj/jGUkcYfUia7HcVEwAp9IvOXQe7voEUXOP3XbqcxCaSE5jzlv5Sh3uWc7VnudhxTDSv0ieyzv8HOQjj/CUht4HYak2BeDJxHodOOP6S8TBo+t+OYo7BCn6h2fwUL/gw9RkL3c91OYxKQjxQe9P+cEz3budk70+045iis0CciVci9C8QLwx91O41JYB87vZkROINbUmbQo7WVk1hl/zOJaPUbsGEunH0vNLPhlKZuPeT7OftoyL9HNrB5cGKUFfoEEL6gyAlNPOx8+RcsLPLjPes2W1DE1LldNOUh3885o0MKLHnO7TimCjYFQgI4tKAIKJp2RUMAAA16SURBVJNSJ9LEs5LfH/8nOtxV8WzeblE3dWWGcybDC//G8PcegG5DoVVXtyOZMHZGn0BGeT5lmHcpT/ivsDtgTT0Tfvn2geAsl9PHQcDvdiATxgp9gsiUb3k49QXynJN4PjDC7TgmCW3dn8KVk7dBcR4PDG1y2OyWNsuleyLquhGR4cBTgBf4t6o+Wun1dOAloD+wE7hSVb8Skc4Elx9cF2r6uareFJ3o5pBUDzyd+ndA+I3v1zj2+9u4IeBj0YWzeCPwD/7w009ZOfD/WKYnHdbMuhDrX7UVQUS8wDPACCALuEpEsio1uxHYrardgL8Cj4W9tlFV+4YeVuTrwCND0+nr2chdvnFs1ePdjmOS3AO+sXytrflb2tO0YK/bcQyRdd3kAIWquklVy4EpwKhKbUYBL4aeTwOGig3zqB/r53LnGem87D+Hd5wct9MYw/c04hbfbWTwHU+lPoMHx+1ISS+SQt8eKArb3hraV2UbVfUDe4BWode6iMhyEflIRAZV9QlEZJyI5IlIXklJSY2+gKRWsh7e+B+WbwvwR/+1bqcx5gcrtSv3+8cy2LuK/5cyze04Sa+uO3O3AR1VtR9wB/CKiDSt3EhVJ6lqtqpmZ2Rk1HGkBHFgN7w6BrxpXPxaKWWkuZ3ImAqmBM7mNf8Qbk2ZwXmeJW7HSWqRFPpioEPYdmZoX5VtRCQFaAbsVNUyVd0JoKpLgY3A4VdnTM0E/DDtBvhuC1w5mS17bPphE4uECf6xrHC68mTqM/SRQrcDJa1ICv0SoLuIdBGRNGAMUHkGo5nA9aHno4EPVFVFJCN0MRcRORHoDmyKTvQkpQqzfgMbP4AL/gKdBrqdyJgjKiONG8t/R4k2499pfyZTvnU7UlKqttCH+tzHA3MJDpWcqqr5IvKQiIwMNXsOaCUihQS7aO4J7R8MrBSRFQQv0t6kqrui/UUkovBpDcIffzqnISx/mYcXlCHZY21qAxPzdtKMX/juIpUAL6Y+RquG9j1b3yIaR6+quUBupX0Twp4fBC6v4rg3gDdqmTEp/TitwY9u9Obyv6mTecX/M57LuZFOOcEfGBuXbGLdRm3PuPI7eCntUeZd1whKd0Gjlm7HShp2Z02cuN47lz+kTmZO4DR+778BsLMiE18Waw/G+e4gK8MDky+Fg3vcjpQ0rNDHgV96Z/Fg6ovMDWRzm+9Wu/PVxK0FTh8um3oAtq+GyZcFz+xNnbOKEdOUW73TuS/1FWYFTucW3234bMJRE+dmb/DD5f+BbV/AC+fD3q/djpTwrNDHqBQPPJryL36bOo03Amdxu+8W/FbkTaLocSFc+wbs2QrPDYMdG9xOlNCs0Meig3uYfXUjxqR8yFP+S/it72YCeN1OZUx0dRkMY2eB/yD8ayhseM/tRAnLCn2s2bYSJg3h7M5e7vSN46/+y7ELryaheFN/HDLcvh9d/rSJFZt34Uy+lLvPSrfpjOuA9QXEClVY9hLk3gmNWjHkxVKKrx7idipjoi/gqzB0WIGrOMjjziQeO+dzzjn7DO7x/ZJlj13jXsYEY2f0seD7b+C1a+Ht26DTGXDTx3xWZIssm+RxgAbc6ruVB33XMdizinfS72F4NzsPjRYr9PWs8h2v1/ROZefD3Ti4+m3unHcQ79i3kCY2sZtJRsILgRGMLH+YndqUOdc0Ci5LuM+mTagt+5VZzw7d8dpTvmJC6ksM8KxlmdONO32/YuOg9nQITeRsd7uaZLVOOzKq/GGuXXwVf/BMh/XvwND7of9Y8NighGNhZ/T1rGMz4dGUSbyddh/dpJh7fTcyuvwBW8zbmDBlpDFhfhnc/Cmc0Btm3wH/GAgFM4PXs0yN2Bl9fdn1JXwykcJbm6Cej3k+MJy/+S9lL43dTmZM7Mr4CVz/NqyZCR/8EaZeB+36wVl3wMkX2Bl+hKzQ1yXHgU3zYclzsH4OeFJ5dqmPV3o/zbYfFuAyxlQpNAzzh02Ba3unMuGnSznx6+vYuMvhyUVlvFvShnWbio7yRsYKfV3YuRHyp8OKV2HXRmjUGs68HXJ+xW0T2tGptxV5Y6pVaRgmwIfAAhyGlefxy+az+fuIDRz07wkuxNP3ajjxbDvLr4IV+mhQ5Zw+mWQ33cHorFSy2wW/0T7e7OfZpeVMK9hEeeBB4EF3cxqTABw8vOPk8E55Dr1kE+etvIerfdNoufoNSvY7zFrv5611fuZt8lPq+/G4E9p3YNvWLe4Fd5EV+mPhOLBjPWxdAps/hY0f8N6l+4AGfOGcyB99A5kdOJ1tbVrB+dD2/B8PtdE0xkTPKj2RWbml/K3XmwwtX8bwBku4rO8KftGvlDJNZZnTnUV6MoucHsx44m6347rGCn01TuzUgeblX5OV4aXn8R76t/UyoL2XZg2CfYc7Sx3e3RjgnY1+CoY/TwktXE5sTPIpJ5U5zgDmOANIwU+OZy1DPcsZ4FnDbZ438aRMp/ye4+DZs4KjeE7oFXy0/gk0bg0JvlJbRIVeRIYDTwFe4N+q+mil19OBl4D+wE7gSlX9KvTavcCNQAC4TVXnRi19bQX8cGAXlO4M3pSxpyg4m953RcHn322mcOwePNIEgHL1UqiZzHK6say8O8u1G196TkC7e9g8/UI6Dbcib4zb/KTwmXMKnzmnANCU/fT3rOekxX/k3rFtYMM8WPHfHw9IOw5adoGWJwY/Nm0PTdrAcScEH03aQEq6S19NdFRb6EOLez8DnAtsBZaIyExVLQhrdiOwW1W7icgY4DHgShHJIriYeE+gHfCeiJykqtG/v993ANbPhfL9oce+is99pVC2Dw5+Fyzs+3cEnx/+FQf/c5tlQrtTefDNtXx7xt2s10w2axubKtiYOLOXxsx3+vGfDx3+9/3pALRpLPQ9wcNJrbx0a1lGt5a76NpiOSe28JDqreLsPr0ZNGwGDZpBg+ahj6FHWmNIaQCpDX/8mNoQUhpCaoPgPk9q8CKxJyX48Fba9qQEt71pwfeLskiqVg5QqKqbAERkCjAKCC/0o4AHQs+nAU9LcFzUKGCKqpYBX4YWD88BFkYnfpjy/fD69RV2BRxlXzns9yn7ypX95bDrgLKjVNlR5qHkez8lpQ47SpVv9ilb9jhs3av4nD3Auh/ep9PpA6Ie1xhTzyqN4lkbeoTz+Bxa+fbif/kG2jZyaHuc0LaJcHzjMpo3KKF5A/nh0SL0sXEaeKLV9dO+P/zyg+i8VxjRau4yE5HRwHBV/Z/Q9nXAAFUdH9ZmdajN1tD2RmAAweL/uapODu1/DpijqtMqfY5xwLjQ5k8Ir7KxozWww+0Q1bCM0WEZay/W80HiZeykqlVOlBUT/RCqOgmY5HaOoxGRPFXNdjvH0VjG6LCMtRfr+SC5MkYy100x0CFsOzO0r8o2IpICNCN4UTaSY40xxtShSAr9EqC7iHQRkTSCF1dnVmozEzjUQT4a+ECDfUIzgTEiki4iXYDuwOLoRDfGGBOJartuVNUvIuOBuQSHVz6vqvki8hCQp6ozgeeAl0MXW3cR/GVAqN1Ughdu/cAtdTLipn7EdNdSiGWMDstYe7GeD5IoY7UXY40xxsQ3m4/eGGMSnBV6Y4xJcFboa0BEHhaRlSKyQkTeFZF2bmeqTESeEJG1oZxvikhztzNVJiKXi0i+iDgiEjPD20RkuIisE5FCEbnH7TxVEZHnReTb0L0rMUdEOojIfBEpCP0f3+52pspEpIGILBaRL0IZY3ZaWRHxishyEZlVfesjs0JfM0+oam9V7QvMAia4HagK84BTVLU3sB641+U8VVkNXAoscDvIIWFTfYwAsoCrQlN4xJr/AMPdDnEUfuC3qpoFnA7cEoP/jmXAz1S1D9AXGC4ip7uc6UhuB9bU9k2s0NeAqu4N22wMxNyVbFV9V1X9oc3PCd67EFNUdY2qxtrdzz9M9aGq5cChqT5iiqouIDiyLSap6jZVXRZ6/j3BIhVTCyJr0L7QZmroEXM/yyKSCVwA/Lu272WFvoZE5BERKQKuITbP6MPdAMxxO0ScaA+Er0e3lRgrUPFGRDoD/YBF7iY5XKhLZAXwLTBPVWMuI/AkcBfg1PaNrNBXIiLvicjqKh6jAFT1PlXtAPwXGH/0d3MnY6jNfQT/jP7vkd/J3YwmcYlIE+AN4DeV/hKOCaoaCHXBZgI5InKK25nCiciFwLequjQa7xcTc93EElU9J8Km/wVygfvrME6VqssoImOBC4Gh6tKNEjX4d4wVNl1HlIhIKsEi/19Vne52nqNR1e9EZD7B6x6xdIH7TGCkiJwPNACaishkVb32WN7MzuhrQES6h22O4vBZTl0XWiTmLmCkqpa6nSeORDLVh6lGaHry54A1qjrR7TxVEZGMQ6PRRKQhwbU2YupnWVXvVdVMVe1M8Hvxg2Mt8mCFvqYeDXU/rASGEbwiHmueBo4D5oWGgT7rdqDKROQSEdkKDARmi4jrq46FLmAfmupjDTBVVfPdTXU4EXmV4HoOPxGRrSJyo9uZKjkTuA74Wej7b0XorDSWtAXmh36OlxDso6/V8MVYZ1MgGGNMgrMzemOMSXBW6I0xJsFZoTfGmARnhd4YYxKcFXpjjElwVuiNMSbBWaE3xpgE9/8BQ8H0G90jctwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate some random numbers following the 'normal' distributed\n",
    "mu = 0\n",
    "sigma = 1\n",
    "normal_rand_nums = np.random.randn(10000) * sigma + mu\n",
    "# calculate and plot the normalized histogram\n",
    "h, b, p = plt.hist(normal_rand_nums, bins=30, density=True, edgecolor='black')\n",
    "\n",
    "# plot the function as a reference\n",
    "x = np.linspace(min(b), max(b), 1000)\n",
    "l = plt.plot(x, st.norm.pdf(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Sampeling a random variable\n",
    "In the abstract an experiment is the set of circumstances needed to generate a probability density funciton (a multidimensional one when multiple readings are involved) for the observed quantities. Then the probability density function is interpreted and \"laws of nature\" are deduced from that. The only way to truly get the probability density is to perform infinite measurements, which, needless to say, is impossible. So instead a **Sample** is taken of the probability density function and these samples are then used to approximate the PDF with the help of models. A sample is nothing more than a group of measurements taken from the PDF.\n",
    "\n",
    "An experiment can have many different parameters $\\vec{\\theta}$ that have to be set before taking a sample. In these cases the PDF is dependent on these parameters $\\vec{\\theta}$. This dependence of the PDF on $\\vec{\\theta}$ is normally formulated as the \"laws\" of nature. When these laws are formulated they don't actually describe a PDF but normally a function $\\vec{y}=f(\\vec{\\theta})$. This is done because the PDF is dependent on the specific measurement setup, the underlying principle that is captured by $f$ however is not.\n",
    "\n",
    "It may be obvious, but it is of course possible to take multiple samples. Sometimes it is useful to change the parameters $\\vec{\\theta}$ before taking another sample, because the goal is to find out how the PDF is dependent on $\\vec{\\theta}$ other times another sample is generated to improve overall *statistical significance* (we will get to that at a later point).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expectation values, means and other Estimators/statistics\n",
    "\n",
    "The defining characteristic of a random variable is its probability density function. However as a rule the PDF is not known and has to be guessed at with the help of models. However there are a few indicators (or estimators/statistics) that can give you an Idea of how a PDF looks like (without using a histogram) based on a sample drawn form it. These characteristic measures are called **estimators** or **statistics**. An estimator can only be calculated from a Sample containing more than one measurement, without giving garbage results.\n",
    "\n",
    "### The mean estimator/statistic\n",
    "The most well known estimator is probably the mean-value estimator.\n",
    "\n",
    "The mean estimator is defined as:\n",
    "\n",
    "$$ mean(\\{x_i\\}) = \\langle x\\rangle = \\frac{1}{n}\\sum_{i=1}^n x_i$$\n",
    "\n",
    "where n is the size of the sample (the amount of measurements in the sample) and $x_i$ is the value of the i-th measurement. For a random vector this is defined as\n",
    "\n",
    "$$ mean(\\{\\vec{x}_i\\}) = \\langle\\vec{x}\\rangle = \\frac{1}{n} \\sum_{i=1}^n \\vec{x}_i$$\n",
    "\n",
    "So for the random numbers (both for the gaussian and the uniform distribution) we can calculate the mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the normal distribution is: -0.00891992692385321\n",
      "The mean of the uniform distribution in [0, 1): 0.4969981166748284\n"
     ]
    }
   ],
   "source": [
    "# calculate the mean for the sample drawn from the gaussian distribution\n",
    "normal_mean = 1/len(normal_rand_nums) * sum(normal_rand_nums)\n",
    "print(\"The mean of the normal distribution is: {}\".format(normal_mean))\n",
    "\n",
    "# calculate the mean for the sample drawn from the uniform distribution\n",
    "uniform_mean = 1/len(rand_nums) * sum(rand_nums)\n",
    "print(\"The mean of the uniform distribution in [0, 1): {}\".format(uniform_mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimators are random variables\n",
    "As the estimator/statistic is a function of a sample drawn from a random variable, it is itself a random variable. This means, that an estimator will vary from sample to sample. As such A sample of estimators can be generated and an estimator can be calculated from that sample (even the same estimator can be calculated on the sample of estimators). This can be done recursively without end (the only problem there is, that to generate a sample of a sample, the amount of measurements that need to be taken grows exponentially with each recursive level. This means that it is unpractical at best to calculate more than an estimator for a sample of samples).\n",
    "So lets try and calculate the mean of means for a uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.47781243525246453, 0.46234426580604626, 0.39296673714384067, 0.5303653962565584, 0.40063917442652636, 0.5806647113280826, 0.5650684736362902, 0.5373535061081534, 0.4735646481405034, 0.43183099985891504, 0.5609222811208431, 0.4151842477218578, 0.5219703454693956, 0.4890117508386031, 0.5582699498822603, 0.433752179296516, 0.503466758420625, 0.5169945672419074, 0.6330541884744576, 0.4741319349757669, 0.5992828190959582, 0.48838846418331705, 0.4774121741079808, 0.5459360264897681, 0.540823200071192, 0.532058794338502, 0.5328762554019367, 0.5064329156864403, 0.502784324013972, 0.5264897447122533, 0.3550892179189121, 0.5468931855474266, 0.4090138037599328, 0.3812258288288738, 0.5722828918779901, 0.3455364371020806, 0.5142588544934021, 0.5346841345099851, 0.5190235011153413, 0.5649426533676345, 0.4921197501425646, 0.621603956349009, 0.5961427022273881, 0.44188597069836977, 0.5970402075867833, 0.5121606209277565, 0.5966806453260758, 0.4479673014250084, 0.38174108802357126, 0.40494927536591535, 0.5547999857497148, 0.549714646074273, 0.38436465984133356, 0.5381079139627105, 0.5192885974035225, 0.4629675377544772, 0.607797001713131, 0.5175216331951755, 0.47850311954420083, 0.3870707589955148, 0.45614669296844534, 0.52620695867555, 0.46936081433161325, 0.42916615905295846, 0.5445597561178835, 0.5319797514555727, 0.5824769102219036, 0.5478093398176372, 0.39069956749062457, 0.478483021832239, 0.5401854853916289, 0.5230514956717595, 0.5724986306161716, 0.4770781399056724, 0.5661066185467976, 0.5095027992509965, 0.3805515589864631, 0.4918528640590549, 0.55006325030042, 0.4878153383213991, 0.36199388388982967, 0.4093712934188134, 0.4728729115071742, 0.5630119287715563, 0.5385040494976796, 0.5145502784166954, 0.4979104058810909, 0.6117432284358375, 0.4978320016907266, 0.5582227248313967, 0.4692079984627127, 0.5214444039669995, 0.3917695528142083, 0.5534056441184011, 0.5948205956978556, 0.5130568205704498, 0.4997616879455474, 0.5295822769414621, 0.5684663422393638, 0.4191256529093817]\n",
      "\n",
      "Mean of means: 0.5018947998535098\n"
     ]
    }
   ],
   "source": [
    "# a sample of random numbers of length n can be drawn from a uniform distribution [0,1) with np.random.rand(n)\n",
    "n = 20\n",
    "mean_sample_size = 100\n",
    "sample_of_means = [np.mean(np.random.rand(n)) for i in range(mean_sample_size)]\n",
    "print(sample_of_means)\n",
    "\n",
    "# now calculate the mean of the means\n",
    "mean_of_means = np.mean(sample_of_means)\n",
    "print(\"\\nMean of means: {}\".format(mean_of_means))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is of course also possible to draw a histogram for the sample of means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAShElEQVR4nO3df7BndX3f8ecrLNoOIXUJV1z3B2sTZIpOQHtnNROTwai4MFRM46RsO3ZttBszMIlTp9W0ndCBf7SdmE6DE7KRHUlH0eYHyVZR3LGxhI4YL2TRBSWLBLJ7BfbKUtHGJiy++8c923y5fu/e7/2e771378fnY+bM95zP+ZxzPp/9Mq97+JzvOSdVhSSpXT+w1g2QJK0sg16SGmfQS1LjDHpJapxBL0mN27DWDRjm3HPPre3bt691MyRp3bjnnnu+UVVTw9adlkG/fft2ZmZm1roZkrRuJHl0sXUO3UhS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGLRn0SbYm+eMkDyS5P8kvd+XnJDmQ5HD3uXGR7Xd3dQ4n2T3pDkiSTm2UM/oTwLur6iLg1cA1SS4C3gt8tqouAD7bLT9HknOA64BXATuA6xb7gyBJWhlLBn1VPVZV93bz3wK+AmwGrgJu6ardArx5yOZvBA5U1fGqego4AOycRMMlSaNZ1hh9ku3AK4AvAOdV1WPdqseB84Zsshk4MrB8tCsbtu89SWaSzMzNzS2nWVLzNm3ZRpI1mTZt2bbW3VdPIz8CIckPAr8PvKuqnk7y/9dVVSXp9aqqqtoL7AWYnp72tVfSgMdnj3D+ez6xJsd+9P1XrslxNTkjndEnOZP5kP9IVf1BV/xEkk3d+k3AsSGbzgJbB5a3dGWSpFUyyq9uAtwMfKWqPjCwaj9w8lc0u4E/GrL5HcBlSTZ2F2Ev68okSatklDP6nwDeCvx0koPddAXwPuANSQ4Dr++WSTKd5EMAVXUcuAH4Yjdd35VJklbJkmP0VXUXkEVWv25I/RngHQPL+4B94zZQktSPd8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhq35ItHkuwDrgSOVdXLu7KPAxd2VV4A/O+qumTIto8A3wKeBU5U1fSE2i1JGtGSQQ98GLgR+J2TBVX1T07OJ/k14Jun2P61VfWNcRsoSepnlFcJ3plk+7B13YvDfw746ck2S5I0KX3H6H8SeKKqDi+yvoDPJLknyZ6ex5IkjWGUoZtT2QXceor1r6mq2SQvBA4k+WpV3TmsYveHYA/Atm3bejZLknTS2Gf0STYA/xj4+GJ1qmq2+zwG3AbsOEXdvVU1XVXTU1NT4zZLkrRAn6Gb1wNfraqjw1YmOSvJ2SfngcuAQz2OJ0kaw5JBn+RW4PPAhUmOJnl7t+pqFgzbJHlxktu7xfOAu5LcB/wp8Mmq+vTkmi5JGsUov7rZtUj524aUfR24opt/GLi4Z/skST15Z6wkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxvV9TLGk1p1xJvPvGFpdL9q8lceO/uWqH7dFBr2kU3v2Gc5/zydW/bCPvv/KVT9mqxy6kaTGGfSS1DiDXpIaZ9BLUuMMeklq3CivEtyX5FiSQwNl/yHJbJKD3XTFItvuTPJgkoeSvHeSDZckjWaUM/oPAzuHlP96VV3STbcvXJnkDOCDwOXARcCuJBf1aawkafmWDPqquhM4Psa+dwAPVdXDVfU3wMeAq8bYjySphz5j9Ncm+VI3tLNxyPrNwJGB5aNd2VBJ9iSZSTIzNzfXo1n6frBpyzaSrPq04fl/d02OK/Ux7p2xvwncAFT3+WvAz/dpSFXtBfYCTE9PV599qX2Pzx5Zs7s1vUtU681YZ/RV9URVPVtV3wV+m/lhmoVmga0Dy1u6MknSKhor6JNsGlj8GeDQkGpfBC5I8pIkzwOuBvaPczxJ0viWHLpJcitwKXBukqPAdcClSS5hfujmEeAXurovBj5UVVdU1Ykk1wJ3AGcA+6rq/hXphSRpUUsGfVXtGlJ88yJ1vw5cMbB8O/A9P72UJK0e74yVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxi0Z9En2JTmW5NBA2X9K8tUkX0pyW5IXLLLtI0m+nORgkplJNlySNJpRzug/DOxcUHYAeHlV/Rjw58CvnGL711bVJVU1PV4TJUl9LBn0VXUncHxB2Weq6kS3eDewZQXaJkmagEmM0f888KlF1hXwmST3JNlzqp0k2ZNkJsnM3NzcBJolSYKeQZ/k3wEngI8sUuU1VfVK4HLgmiQ/tdi+qmpvVU1X1fTU1FSfZkmSBowd9EneBlwJ/LOqqmF1qmq2+zwG3AbsGPd4kqTxjBX0SXYC/wZ4U1X91SJ1zkpy9sl54DLg0LC6kqSVM8rPK28FPg9cmORokrcDNwJnAwe6n07e1NV9cZLbu03PA+5Kch/wp8Anq+rTK9ILSdKiNixVoap2DSm+eZG6Xweu6OYfBi7u1TpJUm/eGStJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNGynok+xLcizJoYGyc5IcSHK4+9y4yLa7uzqHk+yeVMMlSaMZ9Yz+w8DOBWXvBT5bVRcAn+2WnyPJOcB1wKuAHcB1i/1BkCStjJGCvqruBI4vKL4KuKWbvwV485BN3wgcqKrjVfUUcIDv/YMhSVpBfcboz6uqx7r5x4HzhtTZDBwZWD7alX2PJHuSzCSZmZub69EsSdKgiVyMraoCquc+9lbVdFVNT01NTaJZkiT6Bf0TSTYBdJ/HhtSZBbYOLG/pyiRJq6RP0O8HTv6KZjfwR0Pq3AFclmRjdxH2sq5MkrRKRv155a3A54ELkxxN8nbgfcAbkhwGXt8tk2Q6yYcAquo4cAPwxW66viuTJK2SDaNUqqpdi6x63ZC6M8A7Bpb3AfvGap0kqTfvjJWkxhn0ktQ4g16SGmfQS1LjDHpJapxBPyGbtmwjyZpMm7ZsW+vuSzqNjfTzSi3t8dkjnP+eT6zJsR99/5VrclxJ64Nn9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaN3bQJ7kwycGB6ekk71pQ59Ik3xyo86v9myxJWo6xn3VTVQ8ClwAkOQOYBW4bUvVPqsqHsUjSGpnU0M3rgK9V1aMT2p8kaUImFfRXA7cusu7Hk9yX5FNJXrbYDpLsSTKTZGZubm5CzZIk9Q76JM8D3gT87pDV9wLnV9XFwG8Af7jYfqpqb1VNV9X01NRU32ZJkjqTOKO/HLi3qp5YuKKqnq6qb3fztwNnJjl3AseUJI1oEkG/i0WGbZK8KEm6+R3d8Z6cwDElSSPq9YapJGcBbwB+YaDsnQBVdRPwFuAXk5wAvgNcXVXV55iSpOXpFfRV9X+AH15QdtPA/I3AjX2OIen71Bln0g0IrLoXbd7KY0f/ck2OvRJ8Z6yk09Ozz/ge5gnxEQiS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ474xtwRrdKt7abeJSqwz6FqzRreKt3SYutcqhG0lqnEEvSY0z6CWpcQa9JDXOoJekxvUO+iSPJPlykoNJZoasT5L/kuShJF9K8sq+x5QkjW5SP698bVV9Y5F1lwMXdNOrgN/sPiVJq2A1hm6uAn6n5t0NvCDJplU4riSJyQR9AZ9Jck+SPUPWbwaODCwf7cqeI8meJDNJZubm5ibQLK247o7ctZgkjW4SQzevqarZJC8EDiT5alXdudydVNVeYC/A9PR0TaBdWmm+vFlaF3qf0VfVbPd5DLgN2LGgyiywdWB5S1cmSVoFvYI+yVlJzj45D1wGHFpQbT/wz7tf37wa+GZVPdbnuJKk0fUdujkPuK0bM90AfLSqPp3knQBVdRNwO3AF8BDwV8C/6HlMSdIy9Ar6qnoYuHhI+U0D8wVc0+c4kqTxeWesJDXOoJekxhn0ktQ4g16SGuerBCVpocbew2zQS9JCjb2H2aEbSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqXHN3xm7aso3HZ48sXVGSvk80F/SPzx5p6tZlSepr7KGbJFuT/HGSB5Lcn+SXh9S5NMk3kxzspl/t11xJ0nL1OaM/Aby7qu7tXhB+T5IDVfXAgnp/UlWe7krSGhn7jL6qHquqe7v5bwFfATZPqmGSpMmYyK9ukmwHXgF8YcjqH09yX5JPJXnZKfaxJ8lMkpm5ublJNEuSxASCPskPAr8PvKuqnl6w+l7g/Kq6GPgN4A8X209V7a2q6aqanpqa6tssSVKnV9AnOZP5kP9IVf3BwvVV9XRVfbubvx04M8m5fY4pSVqePr+6CXAz8JWq+sAidV7U1SPJju54T457TEnS8vX51c1PAG8FvpzkYFf2b4FtAFV1E/AW4BeTnAC+A1xdVdXjmJKkZRo76KvqLuCUb8+tqhuBG8c9hiSpP591I0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY3r+3LwnUkeTPJQkvcOWf/8JB/v1n8hyfY+x5MkLV+fl4OfAXwQuBy4CNiV5KIF1d4OPFVVPwr8OvD+cY8nSRpPnzP6HcBDVfVwVf0N8DHgqgV1rgJu6eZ/D3hdklO+Z1aSNFmpqvE2TN4C7Kyqd3TLbwVeVVXXDtQ51NU52i1/ravzjSH72wPs6RYvBB4cq2HDnQt8zzHXsdb6A/ZpvbBPp6/zq2pq2IoNq92SxVTVXmDvSuw7yUxVTa/EvtdCa/0B+7Re2Kf1qc/QzSywdWB5S1c2tE6SDcDfA57scUxJ0jL1CfovAhckeUmS5wFXA/sX1NkP7O7m3wL8jxp3rEiSNJaxh26q6kSSa4E7gDOAfVV1f5LrgZmq2g/cDPzXJA8Bx5n/Y7AWVmRIaA211h+wT+uFfVqHxr4YK0laH7wzVpIaZ9BLUuPWddCP8AiGdyb5cpKDSe46eeduku1JvtOVH0xy0+q3fril+jRQ72eTVJLpgbJf6bZ7MMkbV6fFSxu3T+v5e0rytiRzA21/x8C63UkOd9PuhduulZ59enagfOGPMtbMKP/tJfm5JA8kuT/JRwfKT8vvaSxVtS4n5i8Afw34+8DzgPuAixbU+aGB+TcBn+7mtwOH1roP4/Spq3c2cCdwNzDdlV3U1X8+8JJuP2es8z6t2+8JeBtw45BtzwEe7j43dvMb13OfunXfXus+jNmnC4A/O/kdAC88nb+ncaf1fEa/5CMYqurpgcWzgNP9yvMoj5UAuIH55wb934Gyq4CPVdVfV9VfAA91+1trffp0uhq1T8O8EThQVcer6ingALBzhdq5HH36dLoapU//Evhg911QVce68tP1exrLeg76zcCRgeWjXdlzJLmme/TCfwR+aWDVS5L8WZL/meQnV7apI1uyT0leCWytqk8ud9s10qdPsE6/p87PJvlSkt9LcvLmwnX7PXWG9Qng7ySZSXJ3kjevaEtHN0qfXgq8NMn/6tq+cxnbrhvrOehHUlUfrKofAd4D/Puu+DFgW1W9AvhXwEeT/NBatXFUSX4A+ADw7rVuy6Qs0ad1+T11/juwvap+jPmzwVuWqL8enKpP59f8YwT+KfCfk/zIWjRwDBuYH765FNgF/HaSF6xpi1bAeg76UR7BMOhjwJsBuuGNJ7v5e5gfx3vpCrVzOZbq09nAy4HPJXkEeDWwv7t4udx/j9Uydp/W8fdEVT1ZVX/dLX4I+IejbrtG+vSJqprtPh8GPge8YiUbO6JR/q2PAvur6pluyPPPmQ/+0/V7Gs9aXyQYd2L+L/HDzF94PHmh5WUL6lwwMP+PmL9jF2CK7kIl8xdqZoFz1kOfFtT/HH974fJlPPdi7MOcHhdj+/Rp3X5PwKaB+Z8B7u7mzwH+gvkLfBu7+fXep43A87v5c4HDDLngfpr2aSdwy0DbjwA/fLp+T+NOp83TK5erRnsEw7VJXg88AzzF3z5356eA65M8A3wXeGdVHV/9XjzXiH1abNv7k/w34AHgBHBNVT27Kg0/hT59Yn1/T7+U5E3MfxfHmf/FClV1PMkNzD8rCuD69d4n4B8Av5Xku8yPEryvqh5Y9U4sMGKf7gAuS/IA8Czwr6v7v8jT8Xsal49AkKTGrecxeknSCAx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Lj/B/nA9XTVMbhoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h, b, p = plt.hist(sample_of_means, bins=10, edgecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though with \"only\" 100 values for the mean it is fairly difficult to tell from looking at the histogram, how the means are distributed.\n",
    "To improve the quality of the PDF we can calculate the mean more often, this means that we need to generate n times the amount of random numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of mean of means: 0.5003806264178176\n",
      "total amount of random numbers generated: 200000\n"
     ]
    }
   ],
   "source": [
    "# lets track how many numbers are generated in total\n",
    "numbers_generated = 0\n",
    "\n",
    "# definition of the mean\n",
    "def mean(data):\n",
    "    return 1/len(data) * sum(data)\n",
    "\n",
    "# define some functions that calculate the mean and mean of mean and mean of mean of means\n",
    "def generate_mean_from_uniform_sample(n):\n",
    "    global numbers_generated;\n",
    "    numbers_generated = numbers_generated + n\n",
    "    return mean(np.random.rand(n))\n",
    "\n",
    "def generate_sample_of_means(m, n):\n",
    "    return [generate_mean_from_uniform_sample(n) for i in range(m)]\n",
    "\n",
    "def generate_sample_of_mean_of_means(o, m, n):\n",
    "    return [mean(generate_sample_of_means(m, n)) for i in range(o)]\n",
    "\n",
    "mean_of_means_sample_size = 100\n",
    "mean_sample_size = 100\n",
    "sample_size = 20\n",
    "mean_of_means = generate_sample_of_mean_of_means(mean_of_means_sample_size, mean_sample_size, sample_size)\n",
    "print(\"mean of mean of means: {}\".format(mean(mean_of_means)))\n",
    "print(\"total amount of random numbers generated: {}\".format(numbers_generated))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So to generate a single mean of means of means 200,000 random numbers have to be generated. So be careful how deep you do that recursion (and check that what you do makes sense.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The variance Estimator\n",
    "besides the mean of a sample the variance is also a very common metric for distributions. The variance indicates how close the sample values are to the mean. A small variance means that the idividual measurements are close to the mean while a large variance indicates a large distance to the mean.\n",
    "\n",
    "The variance is defined as:\n",
    "\n",
    "$$ \\mathrm{var}(\\{x_i\\}) = \\frac{1}{n} \\sum_{i=0}^n (x_i-\\langle x \\rangle)^2 $$\n",
    "\n",
    "and for vectors\n",
    "\n",
    "$$ \\mathrm{var}(\\{\\vec{x}_i\\}) = \\frac{1}{n} \\sum_{i=0}^n (\\vec{x}_i-\\langle \\vec{x} \\rangle)^2 $$\n",
    "\n",
    "and can also be computed from the measurements in a sample directly (without first calculating the mean) through\n",
    "\n",
    "$$ \\mathrm{var}(\\{x_i\\}) = \\frac{1}{2N^2} \\sum_{i, j=0}^n (x_i - x_j)^2$$\n",
    "\n",
    "also for vectors\n",
    "\n",
    "$$ \\mathrm{var}(\\{\\vec{x}_i\\}) = \\frac{1}{2N^2} \\sum_{i, j=0}^n (\\vec{x}_i - \\vec{x}_j)^2$$\n",
    "\n",
    "The variance can of course also be defined if we don't have a specific sample $\\{\\vec{x}_i\\}$ but also if we have a PDF $p(x)$. In this case the we rely on the previously given definition of the mean of a random variable to calculate it's variance $\\sigma_X^2$.\n",
    "\n",
    "In the next two formulars (and in all following ones) I'll only be writing the vector formulation as the scalar case is simply the vector case wher the vector has length one. To denote a specific value in a measurement vector $\\vec{x}$ I'll write $x_i$. If we take $k$ measurements we produce $k$ measurement vectors $\\vec{x}_j$ when I wan't to denote a specific value in one of these vectors I'll write $(x_i)_j$ which means that I am refering to the i'th scalar in the j'th measurement.\n",
    "\n",
    "$$ \\mathrm{var}(\\vec{X}) = \\sum_{\\vec{x}} (\\vec{x} - \\langle \\vec{X}\\rangle)^2$$ where we sum over all possibilities for x (so $\\forall \\vec{x} \\in S$)\n",
    "\n",
    "We can of course also extend that into the continuous realm with the following expression\n",
    "\n",
    "$$ \\mathrm{var}(\\vec{X}) = \\int_{S_1} ... \\int_{S_n} (x_1 - \\langle X_1\\rangle)^2 ... (x_n - \\langle X_n\\rangle)^2 \\mathrm{d}x_1 ...\\mathrm{d}x_n$$\n",
    "\n",
    "This is essentially a volume integral over the entire subspace that yields a possible experimental outcome.\n",
    "The continuous example only really works, if the pdf is known analytically (either exact or by approximation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So lets give that a spin with the random numbers generated in the histogram example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variance calculated with mean: 0.08361814867579008\n",
      "Variance calculated without mean: 0.08361814867578825\n"
     ]
    }
   ],
   "source": [
    "# define the variance incorporating the mean\n",
    "def varm(data):\n",
    "    return 1/len(data) * sum([(val - mean(data))**2 for val in data])\n",
    "\n",
    "# define the variance without the mean\n",
    "def var(data):\n",
    "    distances = np.ravel([[(vali - valj) for valj in data] for vali in data])\n",
    "    return 1/(2*len(data)**2) * sum (distances**2)\n",
    "\n",
    "print(\"Variance calculated with mean: {}\".format(varm(rand_nums)))\n",
    "print(\"Variance calculated without mean: {}\".format(var(rand_nums)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Expectation value as generalisation of the mean\n",
    "\n",
    "For this part let's assume there is a [random variable](https://en.wikipedia.org/wiki/Random_variable) $X$ with the domain $\\mathbb{R}$ (this means that the random variable can take on any value in $\\mathbb{R}$). Let's further assume there is a function $f$ with the [domain](https://en.wikipedia.org/wiki/Function_(mathematics)) $\\mathbb{R}$ and the [codomain](https://en.wikipedia.org/wiki/Function_(mathematics)) $\\mathbb{R}$. Now I'll ask the question what is the mean of $f(X)$?\n",
    "\n",
    "To answer this question we can use the concept of an expectation value to expand the prevoius definitions of the mean and the variance. The Expectation value is (like the mean) the value that has the same probability density to it's left as it has to it's right (in the one dimensional case). The expectation value is defined (in the case of this (contiguous example) as\n",
    "\n",
    "$$\\mathrm{E}(\\vec{X}) = \\int_{S_{X}} f(\\vec{x})p(\\vec{x}) \\mathrm{d}\\vec{x} $$\n",
    "\n",
    "where $p(x)$ is the PDF for $X$ and $f(x)$ is the aforementioned function. In the case of $X$ having a discrete domain (for example $\\mathbb{N}$) the \n",
    "\n",
    "When we set $f(x) = x$ then we naturally get the mean estimator for the sample.\n",
    "\n",
    "Now we can also set that $f(x)$ to any other function of our liking (like for example a function that describes how the\n",
    "directly measured value is translated by the electronics in the experiment into the value that is actually stored in a long term storage (like a disk drive in a data center).\n",
    "\n",
    "So in short the expectation value is the mean value for an arbitrary function that is applied to the actual values transforming them in any way.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Variance as a different Expectation value\n",
    "\n",
    "If we now alter the definition of $f(x)$ to something along the line of $f(\\vec{x}) = (\\vec{x} - \\langle\\vec{x}\\rangle)^2$ then we get the Expectation value for the variance. Remember that the variance is a random variable when we look at a set of different samples (or measurement runs) and calculate the variance for each of the samples. That means that the variance itsself has a PDF and we can calculate the expectation value for any random variable which will be the variance of the underlying experiments within some margin of uncertainty.\n",
    "\n",
    "So for continuous values we get:\n",
    "\n",
    "$$ E[var(X)] = \\int_{S_X} (\\vec{x} - \\langle \\vec{x} \\rangle)^2 \\cdot p(\\vec{x}) \\mathrm{d}\\vec{x}$$\n",
    "\n",
    "and for a discrete version\n",
    "\n",
    "$$ E[var(X)] = \\sum_{\\vec{x}} (\\vec{x}_i - \\langle \\vec{x} \\rangle)^2 \\cdot p(\\vec{x}) $$\n",
    "\n",
    "This helps to show that the Expectation value can be calculated for any arbitrary estimator."
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}