probability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probability\n",
    "\n",
    "I was promised a book on data analysis and the first chapter is about probability, why?\n",
    "\n",
    "The concept of probability, and how to work with it, is a pre-requisite to perform data analysis. Probabilities are the mathematical tool that allows us to make quantitative statements about data.\n",
    "\n",
    "<!-- Xcheck [test label](testLabel) and {eq}`my_label` and link to another [`.md`](./unfolding.md) and to {figure} -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomness\n",
    "\n",
    "Probability and randomness come as very closed ideas. \n",
    "Let's begin by trying to understand what we mean by randomness.\n",
    "\n",
    "```{margin}\n",
    "A good outline concerning randomness in classical systems can be found in J. Fords article \"How random is a coin toss?\" FIXME [Ref. 1](#refs). The author makes the analogy between quantum mechanics stemming from the finiteness of the Planck's constant, special relativity stemming from the finiteness of the speed of light and the complexity theory stemming from dropping the assumption of infinite precision.\n",
    "```\n",
    "\n",
    "The classical example for randomness comes from tossing a coin where the outcome is head or tail. Because it is a classical system, its outcome can in principle be predicted by evaluating the equations of motion. So how can the aspect of randomness arise from a, at least in principle, deterministic system? The issue is that in order to predict precisely the evolution of a physical system it would be necessary to \"prepare\" it in a configuration with infinitely precise initial conditions, which is in practice not possible. On this case, probability comes in as a convenient tool to predict the outcome of the experiment (head/tail).\n",
    "\n",
    "A different situation is encountered in quantum mechanics when we study the physics at the atomic and sub-atomic scale. Here the governing laws are \"intrinsically\" probabilistic.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability\n",
    "\n",
    "There are several definitions of probability. We will start from the\n",
    "axiomatic definition and then introduce the frequency interpretation.\n",
    "After having defined the conditional probability we will consider the\n",
    "subjective interpretation arising from the Bayes theorem. Later in the\n",
    "course we will address the effect on statistical inference coming from\n",
    "the frequentists/Bayesian interpretation of probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Axiomatic definition\n",
    "\n",
    "We beging from the axiomatic definition of probability (also known as\n",
    "the Kolmogorov axioms). Let $S$ be the set of possible outcomes of an\n",
    "experiment. For every possible outcome $E$ the probability $p(E)$\n",
    "fulfills the following axioms:\n",
    "\n",
    "-   $P(S) = 1$\n",
    "\n",
    "-   $P(E) \\ge 0$, E $\\in$ S\n",
    "\n",
    "-   $P(\\cup E_i) = \\sum P(E_i)$ for any set of disjoint $E_i$ and $E_j$\n",
    "    (i.e. when $E_i$ and $E_j$ are mutually exclusive)\n",
    "\n",
    "It follows immediately that:\n",
    "\n",
    "-   $P(E) = 1-P(E^*)$ where $S = E \\cup E^*$ and $E$ and $E^*$ are\n",
    "    disjoint\n",
    "\n",
    "-   $P(E) \\le 1$\n",
    "\n",
    "-   $P(\\emptyset) = 0$ where $\\emptyset$ is the null set\n",
    "\n",
    "From these axioms we implement working interpretations of probability:\n",
    "the frequency limit and the Bayesian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability as frequency limit\n",
    "\n",
    "The most popular definition of probability is based on the limit of\n",
    "relative frequencies. Assume we conduct an experiment which has a\n",
    "certain number of outcomes (events). Suppose we prepare N identical\n",
    "experiments and find that the outcome $E_i$ occurs $N_i$ times. We\n",
    "assign the probability P$(E_i)$ to the outcome $E_i$, defined by its\n",
    "relative frequency of occurrence:  \n",
    "$P(E_i)=\\lim_{N\\rightarrow \\infty} \\frac{N_i}{N}.$\n",
    "\n",
    "It's easy to verify\n",
    "that this definition satisfies the axioms given above. This definition\n",
    "is also called the objective posterior probability, because the\n",
    "probability is defined a posteriori, i.e. _after the outcomes of the\n",
    "experiment are known_.\n",
    "The \"frequency limit\" approach is very useful in practice, however:\n",
    "\n",
    "-   the limit does not exist in a strict mathematical sense. This is\n",
    "    because there is no deterministic rule linking the outcome of\n",
    "    experiment $j$ with the outcome of experiment $j+1$.\n",
    "\n",
    "-   how do we prepare $N$ identical experiments? Is it sufficient if\n",
    "    they are very similar? (e.g. back to the coin toss: at each toss the\n",
    "    coin has some abrasion and the $(j+1)^{th}$ toss is not identical to\n",
    "    the $j^{th}$).\n",
    "\n",
    "-   nobody can conduct infinitely many experiments. When does the series\n",
    "    converge to the limit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional probability\n",
    "\n",
    "Let A and B be two different events. The probability for A to happen is\n",
    "$P(A)$ and correspondingly $P(B)$ is the probability for B to happen.\n",
    "The probability that either A or B happens is given by:  \n",
    "\n",
    "$$\n",
    "P(A\\,\\cup\\, B)=P(A)+P(B)-P(A\\,\\cap\\,B)\n",
    "$$\n",
    "\n",
    "\n",
    "where $P(A\\, \\cap \\, B)$ denotes the probability that A and B occur together.\n",
    "If A and B are\n",
    "mutually exclusive, then $P(A\\, \\cap \\, B)=0$.\n",
    "\n",
    "Let $P(A|B)$ be the **conditional probability** that event A occurs,\n",
    "given that event B has already occured. Then\n",
    "$P(A\\, \\cap \\, B)=  P(A|B) \\cdot P(B)$: the probability that A and B\n",
    "happen it's the probability that A happens given B, multiplied by the\n",
    "probability that B happens. If the two events are independent, then\n",
    "$P(A|B)=P(A)$, i.e. the occurrence of A does not depend on B, and so\n",
    "$P(A\\, \\cap\\, B)=P(A)\\cdot P(B)$. The conditional probability $P(A|B)$\n",
    "can be written as:  \n",
    "\n",
    "$$\n",
    "P(A|B)=\\frac{P(A\\, \\cap\\, B)}{P(B)}.\n",
    "$$\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Suppose you draw a card from a deck of 52 cards and it's red. \n",
    "What is the probability that it is a diamond ?\n",
    "\n",
    "$\\mbox{P(diamond|red)=}\\frac{\\mbox{P(diamond}\\, \\cap\\, \\mbox{red)}}{\\mbox{P(red)}}$ \n",
    "\n",
    "$\\mbox{P(diamond}\\, \\cap\\, \\mbox{red)}$ = number of red diamonds divided by the total number of cards = 13/52 \n",
    "\n",
    "$\\mbox{P(red)}$ = number of red cards divided by the total number of cards = 26/52 \n",
    "\n",
    "$\\mbox{P(diamond|red)} = \\frac{13/52}{26/52} = 13/26  = 0.5$ \n",
    "\n",
    "which makes sense: 50% of the red cards are diamonds \n",
    "\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Suppose you draw a card from a deck of 52 cards and it's red. \n",
    "\n",
    "What is the probability that it is a queen ?\n",
    "\n",
    "$\\mbox{P(queen|red)}=\\frac{\\mbox{P(queen}\\, \\cap\\, \\mbox{red)}}{\\mbox{P(red)}} $\n",
    "\n",
    "$\\mbox{P(queen}\\, \\cap\\, \\mbox{red)}$ = number of red queens divided by the total number of cards = 2 / 52\n",
    "\n",
    "$\\mbox{P(red)}$ = number of red cards divided by the total number of cards = 26/52 \n",
    "\n",
    "$\\mbox{P(queen|red)}= \\frac{2/52}{26/52} = 1/13$\n",
    "\n",
    "which makes sense: only 2 cards are queen out of the 26 red cards.\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Suppose you draw a card from a deck of 52 cards and it's red. \n",
    "\n",
    "What is the probability that it is a queen of diamonds ?\n",
    "\n",
    "$\\mbox{P(queen of diamonds|red)=}\\frac{\\mbox{P(queen of diamonds}\\, \\cap\\, \\mbox{red)}}{\\mbox{P(red)}} $\n",
    "\n",
    "$\\mbox{P(queen of diamonds}\\, \\cap\\, \\mbox{red)}$ = number of queens of diamonds divided by the total number of cards = 1 / 52\n",
    "\n",
    "$\\mbox{P(red)}$ = number of red cards divided by the total number of cards = 26/52 \n",
    "\n",
    "$\\mbox{P(queen|red)}= \\frac{1/52}{26/52} = 1/26$\n",
    "\n",
    "which makes sense: only 1 cards are queen of diamonds out of the 26 red cards.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes' Theorem\n",
    "\n",
    "The Bayes' theorem formalize the relation between the conditional probability P(A|B) and P(B|A):\n",
    "```{margin}\n",
    "Thomas Bayes, British clergyman \\[1702 - 1761\\]. The so-called \"Bayes' Theorem\" is named after him, but it has been independently re-discovered by Pierre-Simon Laplace \\[1749 -- 1827\\].\n",
    "```\n",
    "$$\n",
    "P(A|B)=\\frac{P(B|A)\\cdot P(A)}{P(B)}\n",
    "$$\n",
    "It can be proven in one line by writing $P(A \\cap B)$ in two different ways:\n",
    "\n",
    "$P(A\\, \\cap\\, B)= P(A|B) \\cdot P(B) =  P(B|A) \\cdot P(A) = P(B\\, \\cap\\, A)$\n",
    "\n",
    "and hence    \n",
    "\n",
    "$P(A|B)=\\frac{P(B|A)\\cdot P(A)}{P(B)}.$  \n",
    "\n",
    "In the general case of $n$-classes of events with the properties $A_{i}$, the theorem\n",
    "generalizes to:  \n",
    "\n",
    "$$\n",
    "P(A_i|B)=\\frac{P(B|A_i)P(A_i)}{\\sum_i P(B|A_i)\\cdot P(A_i)}.\n",
    "$$\n",
    "\n",
    "The theorem states that the probability that A happens given B is equal to\n",
    "the probability that B happens given A (note that A and B are inverted)\n",
    "times your **prior knowledge** about A and divided by a **normalization\n",
    "factor** P(B). The normalization P(B), the probability for B to happen\n",
    "given that A or not-A happens, in all practical cases is expresses as:\n",
    "$P(B) = P(B|A)P(A)+P(B|\\mbox{not-A})P(\\mbox{not-A})$.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Consider a topic coming from virology. We assume that 0.1%\n",
    "of all swans of a certain colony are infected by an influenza virus. A\n",
    "specially developed test for influenza shall, in case of an infected\n",
    "bird, have a detection efficiency of 98%. Unfortunately, the probability\n",
    "of error (false positive) is 3%, which means that the test indicates an\n",
    "infection in 3% of all cases where the bird is not infected. \n",
    "\n",
    "We ask now for the probability $P$ that, after having had a positive test, a swan\n",
    "is actually infected by the influenza virus.\n",
    "\n",
    "According to the given information, we have: \n",
    "\n",
    "$P({\\rm influenza})=0.001$ (0.1% of all swans are infected) \n",
    "\n",
    "$P({\\rm non-influenza})=1-0.001=0.999$. \n",
    "\n",
    "The probabilities for the test response are:\n",
    "\n",
    "$P(+|{\\rm influenza})=0.98$ and\n",
    "\n",
    "$P(-|{\\rm influenza})=1-0.98=0.02$. \n",
    "\n",
    "Where $P(+|{\\rm influenza})$  \\($P(-|{\\rm influenza})$\\) denote the positive \\(negative\\) response\n",
    "under the condition that the swan is actually infected. \n",
    "\n",
    "Furthermore, the probability of a wrong result of the test is given by\n",
    "\n",
    "$P(+|{\\rm non-influenza})=0.03$ and\n",
    "\n",
    "$P(-|{\\rm non-influenza})=1-0.03=0.97$. \n",
    "\n",
    "Therefore, the probability $P({\\rm influenza}|+)$ that a swan is infected, given the test was positive, is\n",
    "\n",
    "$\\begin{aligned}\n",
    "P({\\rm influenza}|+)\n",
    "&=&\\frac{P(+|{\\rm influenza})\\cdot P({\\rm influenza})}{P(+|{\\rm influenza})\\cdot P({\\rm influenza})+P(+|{\\rm non-influenza})\\cdot P({\\rm non-influenza})}\\\\\n",
    "&=&\\frac{0.98\\times 0.001}{0.98\\times 0.001+0.03\\times 0.999}\\\\\n",
    "&=&3\\%.\n",
    "\\end{aligned}$\n",
    "\n",
    "This means that an infection is present in only 3% of all cases in which\n",
    "the test was positive! \n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subjective probability or the Bayesian interpretation\n",
    "\n",
    "In Bayesian inference the probability is interpreted as a subjective\n",
    "\"**degree of belief**\\\" which can be modified by observations (more\n",
    "data). This is the strength of the Bayesian theorem: it provides a\n",
    "quantitative way to update the initial knowledge (prior belief) about a\n",
    "proposition when new data becomes available. Try to compute as an\n",
    "exercise how the $P({\\rm influenza}|+)$ of the previous example changes\n",
    "if you get two consecutive positive tests.\n",
    "\n",
    "A typical application of Bayes theorem in physics is to take the Bayes' theorem and rewrite\n",
    "it interpreting P(A) as the probability (\"belief\") that a theory is\n",
    "correct before doing the experiment; P(B$|$A) = P(result $|$ theory) is\n",
    "the probability of getting the result if the theory is true; P(B) = P\n",
    "(result) is the probability of getting the result irrespective of\n",
    "whether the theory is true or not, and P(A$|$B) = P (theory $|$ result)\n",
    "is our belief in the theory after having obtained the result.  \n",
    "\n",
    "$$\n",
    "P(\\mbox{theory}|\\mbox{result})=\\frac{P(\\mbox{result}|\\mbox{theory}) \\cdot P(\\mbox{theory})}{P(\\mbox{result})}.\n",
    "$$  \n",
    "\n",
    "It's important to notice the *inversion of the logic*: for\n",
    "$P(\\mbox{theory}|\\mbox{result})$ you have collected the data and you are\n",
    "evaluating the probability of the theory to be right; for\n",
    "$P(\\mbox{result}|\\mbox{theory})$ (called the **likelihood**, we will\n",
    "come back to the reason of this name later in these notes when talking\n",
    "about fits) you are estimating the probability to obtain such a data\n",
    "distribution given a certain theory.\\\n",
    "\\\n",
    "The fundamental difference between the frequentist approach and the\n",
    "Bayesian approach relies in the interpretation. The frequentist's\n",
    "probability is interpreted as a *State of Nature*, whereas the Bayesian\n",
    "probability is a *State of Knowledge* inducing inevitably some\n",
    "subjectivity. Thus the probability of an event $P(E)$ depends on the\n",
    "information which is accessible to the observer. The function $P(E)$ is\n",
    "therefore not purely an intrinsic function of the event, it rather\n",
    "depends on the knowledge and information possessed by the observer. A\n",
    "question like \"what is the probability that SuperSymmetry is a true\n",
    "symmetry of nature?\\\" has no meaning in frequentist inference: it is or\n",
    "it is not. There is no probability associated to it. We can instead\n",
    "associate a probability in Bayesian interference, interpreting the\n",
    "probability as degree of belief.\\\n",
    "\\\n",
    "There is no a priori way to assign the \"prior\" assumption $P(A)$ (in the\n",
    "previous example $P(\\mbox{theory})$): the assignment of the prior probability is\n",
    "subjective. The usual prescription is to assume complete ignorance about\n",
    "the prior P(A) and take all values of A as equiprobable.\\\n",
    "There are objections to this postulate:\n",
    "\n",
    "-   if we are completely ignorant about P(A), how do we know it is a\n",
    "    constant?\n",
    "\n",
    "-   a different choice of P(A) would give a different P(A$|$B)\n",
    "\n",
    "-   if we are ignorant about P(A), we are also ignorant about P($A^2$)\n",
    "    or P(1/A), etc\\... taking any of these (P($A^2$) or P(1/A), etc\\...)\n",
    "    as constant would imply a different P(A), giving a different\n",
    "    posterior probability.\n",
    "\n",
    "These objections are usually answered by the assertion (supported by\n",
    "experience) that P(A$|$B) usually converges to about the same value\n",
    "after several experiments irrespective of the initial choice of prior\n",
    "P(A).\n",
    "\n",
    "The distinction between the frequentist and Bayesian approaches to\n",
    "statistical inference reaches its climax when addressing the problem of\n",
    "setting confidence/credible intervals in Sec. [`Confidence Intervals`](./confidenceIntervals.md). The following is an example to give an idea of the problem.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Let's take the measurement of the mass of the electron as\n",
    "an example to understand the difference between the frequentist and\n",
    "Bayesian interpretation of confidence intervals. You measure the mass of\n",
    "the electron to be $520 \\pm 10~keV/c^2$, i.e. you measured $520~keV/c^2$\n",
    "with an apparatus with a resolution of $10~keV/c^2$. \n",
    "It is tempting to conclude that \"the mass of the electron is between 510 and 530 $keV/c^2$ with 68%\n",
    "probability\". This is not the frequentist's meaning of probability. In\n",
    "the frequentist interpretation, the statement that the electron has a\n",
    "certain mass with a certain probability is nonsense. The electron has a\n",
    "definite mass, the problem is that we do not know what the value is. It\n",
    "sounds much more like a Bayesian statement: with a resolution, or\n",
    "\"error\\\", of $\\sigma = 10~keV/c^2$, the probability that we will measure\n",
    "a mass m when the true value is $m_e$ is  \n",
    "\n",
    "$P(m | m_e) \\propto e^{-(m - m_e)^2/2\\sigma^2}$\n",
    "\n",
    "this is the likelihood\n",
    "term. Then by Bayes' theorem, the probability that the true mass has the\n",
    "value $m_e$ after we have measured a value m is  \n",
    "\n",
    "$P(m_e | m) = \\frac{P(m | m_e)P_{prior}(m_e)}{P(m)}$\n",
    "\n",
    "$\\propto P (m | m_e) \\;\\; \\mbox{assuming} \\;\\;\\;P_{prior}(m_e) = const$\n",
    "\n",
    "$\\propto e^{-(m-m_e)^2/2\\sigma^2}$ \n",
    "\n",
    "What we typically state is that\n",
    "\"the 68% credible interval for the mass of the electron is between 510 and $530~keV/c^2$. Note the use of \"credible level\" instead of \"probability\".\n",
    "\n",
    "This and other subtleties will be discussed further when discussing confidence intervals.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability density function\n",
    "\n",
    "We define **random variable** any function of the data. The **event\n",
    "space** is the set of all possible values of a random variable. A random\n",
    "variable which can take any value between two arbitrarily given points\n",
    "in the event space is called a **continuous** variable; conversely, if\n",
    "the variable can only take certain values it is called a **discrete**\n",
    "variable. In the same manner, data described by discrete or continuous\n",
    "variables are called discrete data or continuous data respectively. The\n",
    "distribution $f(x)$ of a random variable $x$ is called **probability\n",
    "density function (p.d.f.)**. $f(x')dx'$ is the probability to find $x$\n",
    "in the interval between $x'$ and $x'+dx'$ and it is normalized\n",
    "$\\int_{-\\infty}^{+\\infty}f(x')dx'=1$ (the probability to find $x$\n",
    "anywhere in its event space is 1). Note that $f(x)$ is *not* a\n",
    "probability but $f(x)dx$ is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative distribution function\n",
    "\n",
    "Let $x$ be a one-dimensional continuous random variable distributed\n",
    "according to $f(x)$. The **cumulative distribution function (cdf)**\n",
    "$F(x')$ gives the probability that the random variable $x$ will be found\n",
    "to have a value less than or equal to $x'$:  \n",
    "\n",
    "$$\n",
    "F(x') = \\int_{-\\infty}^{x'} f(x)dx\n",
    "$$ \n",
    "\n",
    "It follows trivially (with some abuse of notation) that\n",
    "$F(-\\infty)=0$ and $F(+\\infty)=1$. The function $F$ is a monotonously\n",
    "(but not necessarily strictly monotonously) rising function of $x$. \n",
    "\n",
    "The probability density function $f(x)$ is then simply $f(x)=dF(x)/dx$. \n",
    "\n",
    "The function $F$ is dimensionless whereas the function $f$ has dimension\n",
    "$1/x$. The probability to observe the random variable $x$ between two\n",
    "values $x_1$ and $x_2$ can be written in terms of the cdf as:\n",
    "\n",
    "$$\n",
    "P(x_1 \\le x \\le x_2)=\\int_{x_1}^{x_2} f(x')dx'=F(x_2)-F(x_1)\n",
    "$$\n",
    "\n",
    "<!-- \n",
    "\n",
    "The relationship between $f$ and $F$ is depicted in\n",
    "\n",
    "```{figure} ./tmp/Section1Bilder/FWHM.png\n",
    "---\n",
    "name: directive-fig\n",
    "class: bg-primary mb-1\n",
    "height: 150px\n",
    "width: 200px\n",
    "align: center\n",
    "---\n",
    "A density function $f(x)$ as well as its cumulative function $F(x)$.\n",
    "```\n",
    "-->\n",
    "\n",
    "The following is an example to show what is the cumulative of a pdf (chosen to be a gaussian - see [next chapter](./probabilityDistributions.html#gaussian-or-normal-distribution) for more details on the gaussian pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# import relevant libraries\n",
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "x= np.arange(-5,5,0.01) # range of the x axis\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('pdf(x)')\n",
    "ax.plot(x, norm.pdf(x),\"black\") # plot the gaussian pdf (norm.pdf) as a function of x\n",
    "ax.set_ylim(0,0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, norm.pdf(x), 'black')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('pdf(x)')\n",
    "x=np.arange(-5,1,0.01) # cumulative plotted from -5 to 1\n",
    "ax.set_ylim(0,0.5)\n",
    "ax.fill_between(x,norm.pdf(x),alpha=0.5, color='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('cdf(x)')\n",
    "x= np.arange(-5,5,0.01) # range of the x axis\n",
    "ax.plot(x, norm.cdf(x), 'black')\n",
    "ax.set_ylim(0,1.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean, Median and Mode\n",
    "\n",
    "The **arithmetic mean** $\\bar{x}$ (or simply mean value) of a set of $N$\n",
    "numbers $X_{i}$: \n",
    "\n",
    "$\\bar{x}=\\frac{1}{N}\\sum_{i=1}^N X_i$\n",
    "\n",
    "The mean of a function of $x$, ($\\bar{f}$) can be calculated analogously:\n",
    "\n",
    "$\\bar{f}=\\frac{1}{N}\\sum_{i=1}^N f(X_i).$\n",
    "\n",
    "If the $N$ data points are\n",
    "classified in $m$ intervals (i.e. bins of a histogram), and if $n_{i}$ stands for the number of entries in the\n",
    "interval (bin) $i$, then:\n",
    "```{margin}\n",
    "See also the weighted mean in [Measurements uncertainties](./errors.html#weighted-mean-sec-weigthedmean)\n",
    "```\n",
    "\n",
    "$\\bar{x}=\\frac{1}{N}\\sum_{i=1}^{m}n_iX_i.$\n",
    "\n",
    "The **median** of a random variable $x$ divides a frequency distribution\n",
    "into two equally sized halves:\n",
    "\n",
    "$\\int_{-\\infty}^{x_{median}}f(x')dx'=\\int_{x_{median}}^{+\\infty}f(x')dx'=0.5.$\n",
    "\n",
    "The **mode** corresponds to the value of $x$ where the probability\n",
    "density $f(x)$ has a maximum. The mode is not necessarily unique: if a\n",
    "distribution has two maxima, we call it *bimodal*, if it has several\n",
    "maxima, we call it *multimodal*. When only one maximum is present the\n",
    "mode is also called **most probable value**.\n",
    "\n",
    "<!--\n",
    "A sketch to illustrate the meaning of the mean, median, mode is shown in\n",
    "Fig. [1.2](#BarlowComic){reference-type=\"ref\" reference=\"BarlowComic\"}.\\\n",
    "\n",
    "![[\\[BarlowComic\\]]{#BarlowComic label=\"BarlowComic\"} The distribution\n",
    "of the monthly income of Americans around the year 1950. This pictorial\n",
    "representation explains well the differences between mean, mode and\n",
    "median. Tricky question: which of the three describes the most important\n",
    "property of the distribution? Ask yourself what is for you the most\n",
    "important property, and what is the message you want to convey\n",
    "highlighting that. (Taken\n",
    "from [@Barlow])](Section1Bilder/BarlowComic){#BarlowComic\n",
    "width=\"0.5\\\\linewidth\"}\n",
    "\n",
    "-->\n",
    "\n",
    "A rough relation between mode, median and mean (true for unimodal and\n",
    "\"not very skewed\" distributions) is given by\n",
    "\n",
    "$\\mbox{ Mean - Mode}=3\\times\\mbox{ (Mean - Median)}.$\n",
    "\n",
    "So by knowing\n",
    "two out of the three, the third one can be estimated easily by this\n",
    "formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "mu,sigma = 0,5\n",
    "X = np.random.normal(mu, sigma, 11)\n",
    "\n",
    "print (\"array = \", X)\n",
    "print (\"size of the array = \", np.size(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute the mean\n",
    "sum = 0\n",
    "for x in X:\n",
    "    sum += x\n",
    "print (\"mean = \", sum/np.size(X))\n",
    "\n",
    "# or using the numpy mean function\n",
    "print (\"mean = \", np.mean(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the median\n",
    "np.sort(X)\n",
    "print (\"sorted array = \", X)\n",
    "np.median(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the mode of this array:\n",
    "X = [ 1,2,2,4,5,5,5,6,6,7,8,9]\n",
    "\n",
    "from scipy.stats import mode\n",
    "# this will give you back the mode and the number of times it appears\n",
    "\n",
    "print (\"mode\", mode(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the distribution to see visually the mode\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "n, bins, patches = plt.hist(X, facecolor='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **geometric mean** $\\mu_g$ is defined as:\n",
    "\n",
    "$\\mu_g = \\sqrt[N]{x1\\cdot x2\\cdot \\dots \\cdot x_N}.$\n",
    "\n",
    "It is used to\n",
    "characterize the mean of a geometric sequence\n",
    "$(a, ar, ar^2, ar^3, ...)$. The geometric interpretation of the\n",
    "geometric mean of two numbers, a and b, is the length of one side of a\n",
    "square whose area is equal to the area of a rectangle with sides of\n",
    "lengths a and b.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "A population of bacteria grows from 2000 to 9000 in 3 days. \n",
    "What is the daily grow (assuming a constant rate r) ?\\\n",
    "$1^{st}$ day: $n_1$ = 2000 + 2000 r\\\n",
    "$2^{nd}$ day: $n_2$ = n1 + n1 r = 2000 $(1+r)^2$\\\n",
    "$3^{rd}$ day: $n_3$ = n2 + n2 r = 2000 $(1+r)^3$ = 9000\\\n",
    "$\\Rightarrow  1+r = \\sqrt[3]{4.5} \\Rightarrow r = 65.1\\%$\\\n",
    "```\n",
    "\n",
    "The function \n",
    "\n",
    "$f(x) = a \\cdot (1+f)^x$\n",
    "\n",
    "is called geometrical or exponential growth when $x$ is discrete or continuous respectively.\\\n",
    "It covers a particularly important role in finance where it describes\n",
    "the compound interest.\\\n",
    "\\\n",
    "The **harmonic mean** $H$ is defined as:\n",
    "\n",
    "$\\frac{1}{H}=\\frac{1}{N}\\sum_{i} \\frac{1}{X_{i}}.$\n",
    "\n",
    " It is characterize\n",
    "the mean value of a harmonic sequence\n",
    "\n",
    "$\\frac{1}{a}\\;, \\;\\frac{1}{a+d}\\;, \\;\\frac{1}{a+2d}\\;,...\\;,\\; \\frac{1}{a+kd}\\;,...$\n",
    "\n",
    "```{admonition} Example: \n",
    ":class: tip\n",
    "A car travels at 80 km/h for the half of the trip and 100\n",
    "km/h for the second half. What's his average speed? 2/(1/80+1/100) = 89\n",
    "km/h (averaging over periods of time not distances)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantiles\n",
    "\n",
    "Quantiles are values taken from the inverse of the\n",
    "cdf of a random variable. When they are taken at regular intervals \n",
    "they get special names. If the set of data is split into:\n",
    "- two equally sized parts, then the value in the middle is the median. \n",
    "- four equally sized parts, then the four values are called *quartiles* Q1, Q2, Q3 and Q4. The value of Q2 coincides to the median. \n",
    "- ten parts we call them *deciles*\n",
    "- hundred parts we call them *percentile*.\n",
    "\n",
    "The quantiles allow to describe any data distribution without knowing what is its underlying pdf. For example, when a baby is born the pediatrician will compare her/his weight or height with respect to a reference population of babies of the same age, in terms of percentiles. When a baby is in the 20-th weight percentile, it means that 20% of the babies in the reference sample will be lighter.\n",
    "\n",
    "The following is an example of deciles on a gaussian distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.figure(figsize=[5,15])\n",
    "\n",
    "# cumulative of the gaussian pdf\n",
    "plt.subplot(211)\n",
    "x = np.linspace(-5,5,100)\n",
    "plt.plot(x, norm.cdf(x),'k-', alpha=0.6)\n",
    "plt.xlabel(r'x [units]')\n",
    "plt.ylabel(r'cdf(x)')\n",
    "\n",
    "# draw vertical lines\n",
    "for i in range(0,11):\n",
    "    step = i*0.1\n",
    "    plt.plot([-5,5], [step, step], 'k--', lw='0.2')\n",
    "    plt.plot([norm.ppf(step),norm.ppf(step)], [step, 0], 'k--', lw='0.2')\n",
    "\n",
    "# gaussian pdf \n",
    "plt.subplot(212)\n",
    "plt.subplots_adjust(top=0.5)\n",
    "x = np.linspace(-5,5,100)\n",
    "plt.plot(x, norm.pdf(x),'k-', alpha=0.6)\n",
    "plt.xlabel(r'x [units]')\n",
    "plt.ylabel(r'pdf(x)')\n",
    "\n",
    "# write out the corresponding numerical values\n",
    "s =[]\n",
    "p =[]\n",
    "for i in range(0,11):\n",
    "    step = i*0.1\n",
    "    s.append(\"{:.1f}\".format(step))\n",
    "    p.append(\"{:.3f}\".format(norm.ppf(step)))\n",
    "    plt.plot([norm.ppf(step),norm.ppf(step)], [norm.pdf(norm.ppf(step)), 0], 'k--', lw='0.2')\n",
    "\n",
    "print (s)\n",
    "print (p)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expectation value\n",
    "\n",
    "The **expectation value** of a random variable $x$ (or first moment) is defined as:\n",
    "\n",
    "$<x>=\\int_{-\\infty}^{\\infty} x' f(x') dx'$\n",
    "\n",
    "and for discrete variables\n",
    "$r$ as: \n",
    "\n",
    "$<r>=\\sum r_i P(r_i).$\n",
    "\n",
    " where $f(x)dx$ is the probability that\n",
    "the value x is found to be in the interval $(x, x+dx)$ and in the\n",
    "discrete case $P(r_i)$ is the probability that the value $r_i$ occurs.\\\n",
    "\\\n",
    "The expectation value reduces to the arithmetic mean when the\n",
    "probability for any value to occur ($f(x)$ in the continuous case or\n",
    "$P(r)$ in the discrete case) is constant. Take e.g. the discrete case\n",
    "and set the probability for any event $r_i$ $(i = 1,...,N )$ to occurr\n",
    "to be the same. Because $P(r_i)$ is normalized $\\sum_{i=1}^N P(r_i) = 1$\n",
    "then for each $i$, $P(r_i) = 1/N$:\n",
    "\n",
    "$< r > = \\sum r_i P(r_i) = \\frac{1}{N}\\sum r_i$\n",
    "\n",
    "More generally, the expectation value for any *symmetric* probability\n",
    "distribution will coincide with the arithmetic mean. In case of\n",
    "*asymmetric* distributions, the expectation value will be pulled towards\n",
    "the side of the tail (more details on this when we will see some\n",
    "examples of probability distibution functions).\\\n",
    "\\\n",
    "The expectation value of a function $h(x)$ is defined by\n",
    "$<h>=\\int h(x')f(x')dx'$.\\\n",
    "\\\n",
    "The expectation value is a linear operator, i.e.\n",
    "\n",
    "$<a\\cdot g(x)+b\\cdot h(x)> = a  <g(x)>+\\,b  <h(x)>,$\n",
    "\n",
    "but in general $<fg>\\ne<f><g>$. The equality is true only if $f$ and $g$ are\n",
    "independent.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "In quantum mechanics, to compute the expectation value of A,\n",
    "the eigenvalues (outcomes of a measurement) are weighted on their\n",
    "probability to occur:\n",
    "$\\langle A \\rangle_\\psi = \\sum_j a_j |\\langle\\psi|\\phi_j\\rangle|^2$\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variance and Standard Deviation\n",
    "\n",
    "The expectation values of $x^{n}$ and of $(x-<x>)^{n}$ are called the\n",
    "$n^{th}$ *algebraic* moment $\\mu_{n}$ and the $n^{th}$ *central* moment\n",
    "$\\mu'_{n}$, respectively. Central refers to the fact that the\n",
    "expectation value of x is subtracted from each value of x (which in case\n",
    "of a symmetric distribution effectively centers it at zero).\\\n",
    "\\\n",
    "The first algebraic moment $\\mu_{1}$ is equal to the expectation value\n",
    "$<x>$ and it is usually just called $\\mu$. The first central moment is\n",
    "zero by definition. \n",
    "\n",
    "The second central moment is a measure of the width\n",
    "of a probability distribution and is called **variance** $V(x)$. Its\n",
    "square root is called the **standard deviation** $\\sigma$:\n",
    "\n",
    "$\\begin{aligned}\n",
    "V(x) &=& <(x-<x>)^2>              = <x^2 + <x>^2 -2x<x> >= \\\\ \n",
    "     &=&  <x^2> + <x>^2 - 2<x><x> = <x^2>-<x>^2          =\\sigma^2\\end{aligned}$\n",
    "     \n",
    "     \n",
    "where we just used the linearity of the expectation value.\n",
    "\n",
    "We will see when discussing the characteristic function that any a pdf\n",
    "is uniquely described by its moments (see\n",
    "[Characteristic function](./probabilityDistributions.html#characteristic-function-sec-characteristic))\n",
    "\n",
    "\n",
    "It is important to notice that quantities like the variance or the\n",
    "standard deviation are defined using expectation values, and they can\n",
    "only be determined if the \"true\" underlying probability density of the\n",
    "sampling distribution is known. When analysing data we need to distinguish the distribution of the\n",
    "collected data, the (known) **sampling distribution**, from the (unknown) **parent\n",
    "distribution** that we are effectively sampling with our measurements.\n",
    "One of the goals of data analysis (called parameter estimation) is to\n",
    "estimate the parameters of the parent distribution from the properties\n",
    "of the collected data sample. \n",
    "\n",
    "When the the mean $<x>$ of the parent distribution is not known, we define the \n",
    "*sample variance*, commonly called $s^{2}$, as: \n",
    "\n",
    "$\\begin{aligned}\n",
    "\\label{samplevariance}\n",
    "  s^2&=& \\frac{1}{N-1}\\sum_i(x_i-\\bar{x})^2 \\\\\\end{aligned}$.\n",
    "\n",
    "To compute this expression, one has to loop on the data a first time to compute the sample mean, and then a second time to compute the sample variance. A little algebra, transforms the previous expression into \n",
    "\n",
    "$\\begin{aligned}\n",
    "     s^2 &=& \\frac{1}{N-1}\\left(\\sum_i x_i^2-\\frac{1}{N}\\left(\\sum_i x_i\\right)^2\\right)                                                                   \\end{aligned}$\n",
    "     \n",
    "that can be computed by looping only once on the events. This is particularly useful when the number of data points is huge.\n",
    "\n",
    "The value of $s^{2}$ can be understood as the best estimate of the\n",
    "\"true\" variance of the parent distribution. The origin of the factor\n",
    "$\\frac{1}{N-1}$ instead of the usual $\\frac{1}{N}$ will become clear\n",
    "when we will discuss the [theory of estimators](./likelihood.html#properties-of-the-estimators-sec-propestimator)\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "The following table collects some of the quantities defined above,\n",
    "for the Maxwell distribution. The probability density for the magnitude\n",
    "of the velocity $v$ of molecules in an ideal gas at temperature $T$ is\n",
    "given by:\n",
    "\n",
    "$$\n",
    "f(v)=N\\cdot (m/2\\pi k_{B}T)^{\\frac{3}{2}}\\exp(-mv^2/2k_{B}T)\\cdot 4\\pi v^2.\n",
    "$$\n",
    "\n",
    "Here, $m$ is the mass of the molecule and $k_{B}$ is the Boltzmann\n",
    "constant.\n",
    "\n",
    "|  Quantity                         |              Value            |\n",
    "| ----------------------------------|------------------------------ |\n",
    "| Mode (most probable value) $v_m$  |      $(2kT/m)^{1/2}$          |\n",
    "|  Mean $<v>$                       |       $(8kT/\\pi m)^{1/2}$     |\n",
    "|  Median                           |   $v_{median}= 1.098\\cdot v_m$|\n",
    "|  RMS-velocity $v_{rms}$           |         $(3kT/m)^{1/2}$       |\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full Width at Half Maximum\n",
    "\n",
    "Another way to characterize the spread of the data is to compute the\n",
    "Full Width at Half Maximum (FWHM). Given a\n",
    "distribution P(x), it can be computed as the difference between the two\n",
    "values of x at which P(x) is equal to half of its maximum. For a\n",
    "gaussian distribution the relation between FWHM and the standard\n",
    "deviation is: \n",
    "\n",
    "$\\mbox{FWHM} = 2\\sqrt{2\\ln2}\\sigma \\sim 2.355\\sigma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Gaussian mu = 0, sigma = 1\n",
    "max = norm.pdf(0)\n",
    "halfMax = max/2.\n",
    "\n",
    "# Constant array at helfMax\n",
    "f=np.empty(10000)\n",
    "f.fill(halfMax)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "x= np.arange(-5,5,0.001) \n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('pdf(x)')\n",
    "ax.plot(x, norm.pdf(x),'k') # Gaussian\n",
    "ax.plot(x, f,'k') # Constant\n",
    "\n",
    "# Find the intersection between the gaussian and the constant \n",
    "# by finding where the sign of the difference of the two functions changes\n",
    "idx = np.argwhere(np.diff(np.sign(norm.pdf(x) - f))).flatten() \n",
    "plt.plot(x[idx], f[idx], 'ro')\n",
    "ax.set_ylim(0,0.5)\n",
    "plt.show()\n",
    "\n",
    "print (\"Half Maximum = \", halfMax)\n",
    "print (\"Intersections: \", x[idx])\n",
    "print (\"FWHM = \", x[idx[1]]-x[idx[0]])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher Moments\n",
    "\n",
    "The third moment is called **skewness** $\\gamma_1$ and it is often\n",
    "defined as \n",
    "\n",
    "$\\gamma_1=\\mu'_3/\\sigma^3=\\frac{1}{\\sigma^3}<(x-<x>)^3>=\n",
    "\\frac{1}{\\sigma^3}(<x^3>-3<x><x^2>+2<x>^3).$\n",
    " \n",
    "\n",
    "The quantity $\\gamma_1$ is\n",
    "dimensionless and characterizes the skew. It is zero for symmetric\n",
    "distributions, and positive(negative) for asymmetric distributions with\n",
    "a tail to the right(left).\n",
    "\n",
    "Another definition of skewness often used is the *Pearson's skew* given\n",
    "by $(mean-mode)/\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import skewnorm\n",
    "\n",
    "fig, (axl, axr) = plt.subplots(ncols=2, figsize=(15, 6))\n",
    "\n",
    "x = np.linspace(-5,5,1000)\n",
    "yl_ = skewnorm.pdf(x, a=0)\n",
    "yl = skewnorm.pdf(x, a=-1)\n",
    "yr_ = skewnorm.pdf(x, a=0)\n",
    "yr = skewnorm.pdf(x, a=1)\n",
    "\n",
    "axl.plot(x, yl_, 'k-', alpha=0.6, label='skewness=0')\n",
    "axl.plot(x, yl, 'r', linewidth=2, label='skewness=-1')\n",
    "axr.plot(x, yr_, 'k-', alpha=0.6, label='skewness=0')\n",
    "axr.plot(x, yr, 'r', linewidth=2, label='skewness=1')\n",
    "axl.legend()\n",
    "axr.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **kurtosis** $\\gamma_2$ is defined as\n",
    "$\\gamma_{2}=\\mu'_4/\\sigma^4-3$, which is a (dimensionless) measure of\n",
    "how the distribution behaves in the tails compared to its maximum. The\n",
    "factor $-3$ is subtracted to obtain a kurtosis of zero for a Gaussian\n",
    "distribution. A positive $\\gamma_{2}$ means that the distribution has a\n",
    "larger maximum and larger tails than a Gaussian distribution with the\n",
    "same values for mean and variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm, kurtosis\n",
    "\n",
    "sigmas = [0.2, 0.5, 1, 1.5, 2]\n",
    "size = 10000\n",
    "\n",
    "fig, ax = plt.subplots(figsize=[10,6])\n",
    "for sigma in sigmas:\n",
    "    data = norm.rvs(size=size, loc=0, scale=sigma)\n",
    "    k = kurtosis(data)\n",
    "    x = np.linspace(-5, 5, size)\n",
    "    y = norm.pdf(x, loc=0, scale=sigma)\n",
    "    ax.plot(x, y, label=\"k={:.2f}\".format(k))\n",
    "    ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some useful scipy functions to compute moments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "mu,sigma = 0,5\n",
    "X = np.random.normal(mu, sigma, 101)\n",
    "#print (X)\n",
    "\n",
    "from scipy.stats import moment\n",
    "for i in range(1,5):\n",
    "    print (i,\"-moment\", moment(X, moment=i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import skew\n",
    "skew(X)\n",
    "\n",
    "# left tail\n",
    "Xleft = np.array([0,1,1,1,3,3,3,3,3,4])\n",
    "print (\"left\", skew(Xleft))\n",
    "\n",
    "# right tail\n",
    "Xright = np.array([0,1,1,1,3,4,5,6])\n",
    "print (\"right\", skew(Xright))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import kurtosis\n",
    "kurtosis(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful Inequalities\n",
    "\n",
    "Two useful inequalities can be used to estimate the upper limits\n",
    "probabilities if the underlying distribution is not known. Both\n",
    "inequalities are given without proof.\\\n",
    "\\\n",
    "Let $x$ be a positive random variable. Then it holds:\n",
    "\n",
    "$P(x\\geq a) \\leq \\frac{<x>}{a}.$\n",
    "\n",
    "This inequality provides an upper limit for the probability \n",
    "of random events located in the tails of the distribution.\n",
    "\n",
    "Let $x$ be a positive random variable. Then it holds:\n",
    "\n",
    "$P(\\,|~x~-<x>|\\geq k)\\leq \\frac{\\sigma^2}{k^2}$ \n",
    "\n",
    "So for example the probability for a result to deviate for more than three standard\n",
    "deviations from the expectation value, is always smaller than 1/9,\n",
    "independent of the underlying probability distribution. The inequality\n",
    "is true in general, but it gives quite a weak limit (for a Gaussian\n",
    "distribution, the probability to lie outside three standard deviations\n",
    "is about 0.0027, see next chapter). It is nevertheless useful for\n",
    "qualitative considerations, if nothing at all is known about the\n",
    "underlying distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More than one dimension\n",
    "\n",
    "When doing measurements we are often performing a sampling of a\n",
    "multi-dimensional dimensional pdf. For example in the case of a 2D\n",
    "pdf, the probability to observe $X_1$ in $[x_1, x_1+dx_1]$ and $X_2$ in\n",
    "$[x_2, x_2+dx_2]$ is: \n",
    "\n",
    "$\\begin{aligned}\n",
    "P(x_1< X_1 < x_1+dx_1, x_2 < X_2 < x_2+dx_2) &=& f(x_1,x_2)dx_1dx_2\\\\\n",
    "P(a< X_1 < b, c < X_2 < d) &=& \\int_a^b dx_1 \\int_c^d dx_2 f(x_1,x_2)\\\\\\end{aligned}$\n",
    "\n",
    "For multi-dimensional pdfs we can define the concept of marginal and\n",
    "conditional pdf. To get the idea let's take a 2D example:\n",
    "\n",
    "-   **Marginal** pdf: it's the pdf describing $x_1$ independently of the\n",
    "    value of $x_2$; the 2D probability can be \"marginalised / projected\"\n",
    "    by integrating over one variable \n",
    "\n",
    "    $\\begin{aligned}\n",
    "    f_1(x_1) &=& \\int_{-\\infty}^{+\\infty}f(x_1,x_2)dx_2 \\\\\n",
    "    f_2(x_2) &=& \\int_{-\\infty}^{+\\infty}f(x_1,x_2)dx_1 \\\\\n",
    "    \\end{aligned}$\n",
    "\n",
    "-   **Conditional** pdf: it's the pdf of $x_2$ for a fixed value of\n",
    "    $x_1$: $f(x_2|x_1)$\n",
    "\n",
    "    $\\begin{equation}\n",
    "    P(a<X_2<b|X_1 = x1) = \\int_a^bf(x_2|x_1)dx_2\n",
    "    \\end{equation}$\n",
    "\n",
    "In plain english: the marginal pdf ignores the values of the other\n",
    "variable(s), the conditional pdf fix the value of the other variable(s).\n",
    "\n",
    "FIXME Normalization Coditional probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import axes3d\n",
    "\n",
    "# multivariate gaussian mu = 0; sigma = 1\n",
    "def gg(x, y):\n",
    "    return 1./(2*math.pi) * np.exp(-0.5*(x*x +y*y))\n",
    "\n",
    "# gaussian mu = 1; sigma \n",
    "# def g(x, sigma):\n",
    "#     return 1./math.sqrt(2*math.pi)/sigma * np.exp(-0.5*(x*x))\n",
    "\n",
    "x = np.linspace(-5, 5, 1000)\n",
    "y = np.linspace(-5, 5, 1000)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "G = gg(X, Y)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax1 = plt.axes(projection='3d')\n",
    "ax1.plot_surface(X, Y, G,cmap='Greys', linewidth=0)\n",
    "ax1.set_xlabel('X')\n",
    "ax1.set_ylabel('Y')\n",
    "ax1.set_zlabel('Z')\n",
    "\n",
    "# Conditional distribution\n",
    "# p(X=x, y=1)\n",
    "xline = x\n",
    "yline = -1.+0.*y\n",
    "zline = gg(xline,yline)\n",
    "\n",
    "ax1.plot3D(xline, yline, zline, 'red', label='p(X=x,Y)',linewidth=3.0)\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax2 = fig.gca(projection='3d')\n",
    "#ax1 = plt.axes(projection='3d')\n",
    "ax2.plot_surface(X, Y, G, cmap='Greys')\n",
    "cset = ax2.contourf(X, Y, G, zdir='x', offset=-5, cmap='Greys')\n",
    "cset = ax2.contourf(X, Y, G, zdir='y', offset=5, cmap='Greys')\n",
    "ax2.set_xlim(-5, 5)\n",
    "ax2.set_ylim(-5, 5)\n",
    "ax2.set_xlabel('X')\n",
    "ax2.set_ylabel('Y')\n",
    "ax2.set_zlabel('Z')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformation of variables\n",
    "\n",
    "Let's consider a 2-dimensional event space as an example (easily\n",
    "generalizable to the N-dimensional case) and call the two random\n",
    "variables $(x,y)$. How does the pdf $f(x,y)$ transforms under a change\n",
    "of variable to $(X,Y)$? \n",
    "\n",
    "```{image} ./images/ch1/transformation.png\n",
    ":width: 400px\n",
    ":align: center\n",
    "```\n",
    "\n",
    "Consider a small interval $A$ around a point $(x,y)$ that we will\n",
    "transform to a small interval $B$ around the point $(X,Y)$ such that\n",
    "$P(A)=P(B)$:\n",
    "\n",
    "$$\n",
    "P\\left[(X,Y) \\in B\\right] = P\\left[(x,y) \\in A \\right] = \\int_A f(x,y) dx dy.\n",
    "$$\n",
    "\n",
    "\n",
    "Assume that the transformation \n",
    "\n",
    "$$\n",
    "X = u_x(x,y)\\\\\n",
    "Y = u_y(x,y)\n",
    "$$\n",
    "\n",
    "is bijective so that the inverse exist\n",
    "(together with the first derivative): \n",
    "\n",
    "$$\n",
    "x = w_x(X,Y)\\\\\n",
    "y = w_y(X,Y).\n",
    "$$\n",
    "\n",
    "Then \n",
    "\n",
    "$$\n",
    "P(A) = P(B)\\\\\n",
    "\\int_A f(x,y) dx dy &=& \\int_B f(w_x(X,Y), w_y(X,Y)) |J| dX dY\n",
    "$$\n",
    "\n",
    "where J is the Jacobian determinant: \n",
    "\n",
    "$$\n",
    "J = \n",
    "\\left| \\begin{array}{cc}\n",
    "\\frac{\\partial w_x}{\\partial x} \\frac{\\partial w_x}{\\partial y}\\\\\n",
    "\\frac{\\partial w_y}{\\partial x} \\frac{\\partial w_y}{\\partial y}\n",
    "\\end{array}\\right|\n",
    "$$\n",
    "\n",
    "So the p.d.f. in (X,Y) is the p.d.f. in (x,y) times the Jacobian: \n",
    "\n",
    "$$\n",
    "g(X,Y) = f(x,y) |J| = f(w_x(X,Y), w_y(X,Y)) |J|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Covariance and Correlation\n",
    "\n",
    "Let $f(x_1,x_2) dx_1 dx_2$ be the joint probability to observe\n",
    "$x_1\\in [x_1,x_1+dx_1]$ and $x_2\\in [x_2, x_2+dx_2]$. The two variables\n",
    "$x_1$ and $x_2$ are **independent** if and only if they fulfill the\n",
    "following relation: \n",
    "\n",
    "$f(x_1,x_2)=f(x_1)\\cdot f(x_2).$\n",
    "\n",
    "If this is the\n",
    "case, the two variables are said to be **uncorrelated**. If the above\n",
    "condition is not fulfilled, then the variables are dependent, i.e.\n",
    "**correlated**.\n",
    "\n",
    "The **covariance** $cov(x_{1},x_{2})$ between two variables is defined\n",
    "as \n",
    "\n",
    "$cov(x_1,x_2)=<(x_1-<x_1>)\\cdot (x_2-<x_2>)>=<x_1x_2>-<x_1><x_2>.$\n",
    "\n",
    "If two variables are independent then $<x_1x_2>=<x_1><x_2>$ and so the\n",
    "covariance is zero.\\\n",
    "Knowing the covariance between two variables, we can write the general\n",
    "expression of the variance of their sum:\n",
    "\n",
    "$V(x_1+x_2)=V(x_1)+V(x_2)+2\\cdot cov(x_1,x_2).$ \n",
    "\n",
    "Some more properties of the covariance are:\n",
    "($a\\in \\mathbb{R}$)\n",
    "\n",
    "-   $cov(x,a) = 0$\n",
    "\n",
    "-   $cov(x,x) = V(x)$\n",
    "\n",
    "-   $cov(x,y) = cov(y,x)$\n",
    "\n",
    "-   $cov(ax,by) = ab~cov(x,y)$\n",
    "\n",
    "-   $cov(x+a,y+b) = cov(x,y)$ is translation invariant (shift origin)\n",
    "\n",
    "-   $cov(x,y)$ has units !\n",
    "\n",
    "The **correlation** coefficient between $x_{1}$ and $x_{2}$ is defined\n",
    "as: \n",
    "\n",
    "```{math}\n",
    ":label: rho\n",
    " \\rho_{x_1x_2} = \\frac{cov(x_1,x_2)}{\\sqrt{V(x_1)V(x_2)}}\n",
    "```\n",
    "\n",
    "The correlation coefficient is the covariance normalized by the square root\n",
    "of the product of the variances, to get a value between +1 and -1.\n",
    "\n",
    "If two variables are independent, given that their covariance is zero,\n",
    "also $\\rho_{x_1 x_2}=0$. The inverse is not necessarily true. This means\n",
    "that $\\rho_{x_1 x_2}=0$ can hold but $x_{1}$ and $x_{2}$ can nevertheless\n",
    "be dependent, as illustrated by the examples in {numref}`fig-wikiCorr`.\n",
    "\n",
    "```{figure} ./images/ch1/wikiCorr.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-wikiCorr\n",
    "---\n",
    "Several sets of (x, y) points, with the Pearson correlation coefficient of x and y for each set. \n",
    "[wikipedia](https://en.wikipedia.org/wiki/Correlation)\n",
    "```\n",
    "\n",
    "\n",
    "Note in particular the second raw where the correlation coefficient \n",
    "is not defined for the horizontal line because one of the variables \n",
    "has null variance; and the third raw, where the correlation coefficient is always\n",
    "zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a sample of size $n$\n",
    "($(x_{1},y_{1}),(x_{2},y_{2}),\\ldots, (x_{n},y_{n})$), the *sample\n",
    "covariance* or empirical covariance $s_{xy}$, which is the best estimate\n",
    "for the (true) covariance $s_{xy}$ is given by:\n",
    "\n",
    "$$\n",
    "s_{xy}=\\frac{1}{n-1}\\sum_i(x_i-\\bar{x})(y_i-\\bar{y}).\n",
    "$$\n",
    "\n",
    "and the empirical correlation $r_{xy}$ also called\n",
    "Pearson-correlation-coefficient gives the best estimate for the (true)\n",
    "correlation coefficient $\\rho_{xy}$:\n",
    "\n",
    "$$\n",
    "r_{xy}=\\frac{s_{xy}}{\\sqrt{s_x} \\sqrt{s_y}}.\n",
    "$$ \n",
    "\n",
    "Here, the sample standard deviations are labeled with $s_{x}$ and $s_{y}$,\n",
    "respectively.\n",
    "\n",
    "A word of caution: \"correlation is not causation\". The fact that one can\n",
    "find a correlation between two observables does not necessarily means\n",
    "that there is a causality relation between the two (see {numref}`fig-correlationCausation`)\n",
    "\n",
    "```{figure} ./images/ch1/chocoNobel.jpg\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-correlationCausation\n",
    "---\n",
    "\"Chocolate consumption, cognitive function, and Nobel laureates.\n",
    "*N. Engl. J. Med. 2012 Oct 18; 367(16):1562-4.*\n",
    "```\n",
    "\n",
    "\n",
    "## References\n",
    "\n",
    "Most of the material of this section is beautifully exposed in:\n",
    "\n",
    "-   R. Feynman [@Feynman], \"Feynamn lectures on physics\": Ch.6\n",
    "\n",
    "-   W. Metzger [@Metzger], \"Statistical Methods in Data Analysis\": Ch.2\n",
    "\n",
    "-   L. Lyons [@Lyons], \"Statistics for Nuclear and Particle Physicist\":\n",
    "    Ch. 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problems\n",
    "\n",
    "1. Take as starting point the example in the Bayes' theorem section and suppose that a swan is tested positive twice in two subsequent tests. What would then be the probability that the swan is actually infected by the influenza virus."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}