monteCarlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo methods\n",
    "===================\n",
    "\n",
    "In this chapter we will discuss an introduction to the so-called Monte\n",
    "Carlo (MC) simulations, which use random numbers and a sequential\n",
    "description of the operations as basic concepts. Because this method\n",
    "uses principles from probability calculations and statistics, it is also\n",
    "known as the method of statistical trials. In the next sections, after a\n",
    "digression on random number generators, we will describe two of the most\n",
    "common applications of Monte Carlo methods: **integration** and\n",
    "**simulation**.\n",
    "\n",
    "```{figure} ./images/ch4/needledisplay2.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-needles\n",
    "---\n",
    "FIXME Ref: Randomly distributed needles on a strip which has the width of the\n",
    "length of one needle. 81 out of 128 needles cross the border of the\n",
    "strip in this picture. (wiki)\n",
    "```\n",
    "\n",
    "```{margin}\n",
    "George Louis Leclerc, Duke of Buffon (1707 - 1788), a French natural scientist\n",
    "```\n",
    "Probably the oldest mentioning of a Monte Carlo method, which\n",
    "illustrates all its basic elements, is known as the needle experiment by\n",
    "Buffon. The duke baffled his colleagues in 1777 with a simple method\n",
    "to get the number $\\pi$ by simply counting the number of needles thrown\n",
    "onto a strip with the same width as the length of the needles $(l)$. He\n",
    "found out that the ratio between the number of needles $(k)$ crossing\n",
    "the border of the strip (the slightly darker ones in {numref}`fig-needles`)\n",
    "and the total number of thrown needles $(n)$ is exactly $2 / \\pi$ (i.e.\n",
    "$k / n = 2 / \\pi = p$). This value is calculated analytically using the\n",
    "position dependent probability density to cross the border of the strip,\n",
    "which is an $arccos$ function, mirrored in the middle of the strip (see\n",
    "{numref}`fig-needlesmc`. The picture on the right hand side in\n",
    "{numref}`fig-needlesmc` shows the convergence to the exact value of $\\pi$ with increasing number\n",
    "of thrown needles. The dashed line in this figure is the expected error\n",
    "for the value of $\\pi$, which is calculated using the binomial\n",
    "distribution for $k$ by including the error propagation to be\n",
    "$\\frac{2n}{k^{2}} \\sqrt{np(1-p)} = 2.37 / \\sqrt{n}$.\n",
    "\n",
    "```{figure} ./images/ch4/needlemcEnglishRotatedPNG.png\n",
    "---\n",
    "width: 800px\n",
    "align: center\n",
    "name: fig-needlesmc\n",
    "---\n",
    "FIXME Ref: Probabilities for needles to cross the border of the strip depending\n",
    "on their position between the two borders (left), and the result of a MC\n",
    "simulation with its error (dashed line). (wiki)\n",
    "```\n",
    "Another example for the determination of $\\pi$ is to inscribe a circle\n",
    "in a square and drop objects randomly on it (e.g. drops of rain, see\n",
    "{numref}`fig-piMC`. From\n",
    "the ratio of the number of drops ending in the circle to the total\n",
    "number of drops in the square one finds that $\\pi = 4 * in / all$. With\n",
    "$10^6$ drops we found a value of 3.13954.\n",
    "\n",
    "```{figure} ./images/ch4/piMC.png\n",
    "---\n",
    "width: 300px\n",
    "align: center\n",
    "name: fig-piMC\n",
    "---\n",
    "FIXME: The circle inscribed in a square used to calculate $\\pi$.\n",
    "```\n",
    "\n",
    "As of today, the MC methods are preferably used in numerical mathematics\n",
    "if the formulation of the stochastic model is simpler than the\n",
    "formulation of the analytic model for the numerical solution of the\n",
    "problem. Monte Carlo methods are used in many different areas of\n",
    "research. To name just a few of them:\n",
    "\n",
    "-   Numerical problems, such as the calculation of integrals or the\n",
    "    solution of ordinary or partial differential equations.\n",
    "\n",
    "-   Quality control of products, for example the determination of the\n",
    "    lifetime of light bulbs.\n",
    "\n",
    "-   Problems from operations research, such as transport problems.\n",
    "\n",
    "-   Decision management using simulations or risk analysis in investment\n",
    "    banking.\n",
    "\n",
    "Monte Carlo methods can generally be spilt into three major steps:\n",
    "\n",
    "-   A stochastic model has to be found for the original mathematical\n",
    "    model, which describes the problem accurately enough.\n",
    "\n",
    "-   A sequence of random numbers has to be generated, which are then\n",
    "    used to simulate a realistic situation, and hence have the same\n",
    "    underlying distribution.\n",
    "\n",
    "-   Estimates for the original problem have to be found using the\n",
    "    results coming from the random numbers.\n",
    "\n",
    "Monte Carlo was eponymous for this process: the problems connected to\n",
    "gambling were motivating enough to start thinking about randomness of\n",
    "events."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Pseudo)Random Numbers Generators \n",
    "---------------------------------\n",
    "\n",
    "```{margin}\n",
    "http://www.random.org offers true random numbers generated from atmospheric noise.\n",
    "```\n",
    "\n",
    "True random generators are based on physical processes. Examples can be\n",
    "the noise level across a resistor, the time between the arrival of two\n",
    "cosmic rays or the number of radioactive decays in a fixed time\n",
    "interval. A method from the early applications of Monte Carlo technique\n",
    "in high energy physics, used the stopping azimuthal position of a\n",
    "cylinder which had been put in rotation by a motor (activated by an\n",
    "operator) and turned off by a cosmic ray recorded by a detector.\n",
    "\n",
    "The main issues with these kind of random number generators is that they\n",
    "are very slow. The solution adopted is to generate random numbers\n",
    "running an algorithm on a computer. By construction the obtained random\n",
    "sequence is *not* random (that's why they are called pseudo-random).\n",
    "We'll see in the following that the two most important parameters of any\n",
    "generator are the period length (how many numbers it can generate before\n",
    "it starts to repeat itself) and the correlation between the generated\n",
    "numbers.\n",
    "\n",
    "A simple and classic generator is the general *linear congruential\n",
    "generator*:\n",
    "\n",
    "$$\n",
    "n_{i+1}=(a\\cdot n_{i}+c)\\,\\,\\,\\mathrm{mod}\\,\\,\\,m\\hspace{8mm} u_{i}=n_{i}/m\n",
    "$$\n",
    "\n",
    "it generates uniformly distributed numbers in the interval (0\\...1\\].\n",
    "The boundary value 0 is usually not included to avoid the divide-by-zero\n",
    "danger in case the number is carelessly used in further calculations. We\n",
    "will call such a uniform probability density $U(0,1)$: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "U(0,1) =\n",
    "\\begin{cases}\n",
    " 1 \\qquad &\\text{if} \\hspace{2mm}0<u\\le1,  \\\\\n",
    " 0 &\\text{else}\n",
    "\\end{cases}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The algorithm uses three integer constants:\n",
    "the multiplicand $a$, the summand $c$ and the module $m$. Generators\n",
    "with the summand $c = 0$ are called *multiplicative linear congruential\n",
    "generators*. The initial value $n_{1}$ is also called *seed*: the choice\n",
    "of the initial value allows to steer the generation process. The\n",
    "distribution and the correlation among the first 10'000 values using the\n",
    "values $m = 2^{31}$, $a=65539$ and $c=0$ is shown in {numref}`fig-mulincongPNG`\n",
    ". It was first used in the sixties by IBM and\n",
    "became famous under the name of *RANDU*, but it had bad performance.\n",
    "\n",
    "```{figure} ./images/ch4/mulincongPNG.png\n",
    "---\n",
    "width: 600px\n",
    "align: center\n",
    "name: fig-mulincongPNG\n",
    "---\n",
    "FIXME ref: Histogram (100 bins) of the first 10'000 values generated with *RANDU*\n",
    "and correlations among three sequential values (not binned). 10'000\n",
    "values, generated with the Mersenne Twister-function, have been plotted\n",
    "for comparison in the same way on the right hand side\n",
    "graph. [wiki]\n",
    "```\n",
    "\n",
    "```{margin}\n",
    "G. Marsaglia, \\\"Random numbers fall mainly in the planes\\\", Proc.\n",
    "Natl. Acad. Sci. 61(1), 25--28 (1968)\n",
    "It is a general property of a linear generator that a sequential\n",
    "```\n",
    "$k$-tuple of random numbers lies in the $k$-dimensional space on\n",
    "$(k-1)$-dimensional hyperplanes. The maximal distance between these\n",
    "hyperplanes is an important test for linear generators (spectral test).\n",
    "The graph on the right hand side in {numref}`fig-mulincongPNG`\n",
    "compares *RANDU* (graph in the middle) to a\n",
    "far more uniform distribution from the Mersenne Twister generator.\n",
    "\n",
    "Nowadays algorithms exist with a period length of $2^{19937}$ and which\n",
    "are (for most practical purposes) uncorrelated. The random generators\n",
    "which are implemented in many computer programs are usually sufficient\n",
    "for daily use. Nevertheless, in some special cases such as lattice QCD\n",
    "calculations, far better generators are needed.\n",
    "\n",
    "Two tricks used to get random numbers with minimal correlation and\n",
    "extraordinary long period length are:\n",
    "\n",
    "-   **Combination:** Two random numbers are generated with a generator\n",
    "    each, and a new one is generated using the operations $+$, $-$ or\n",
    "    *exclusive-OR* at bit level.\n",
    "\n",
    "-   **Rearrangement:** The memory is filled with some random numbers,\n",
    "    and the result of another generator is used to determine the address\n",
    "    in the memory for the next random number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tests of Random Generators\n",
    "\n",
    "For an extensive overview see Knuth [@KnuthBook], here we will briefly\n",
    "touch upon the most important ones:\n",
    "\n",
    "-   **Test for uniform distribution.** The interval \\[0,1\\] is divided\n",
    "    into $k$ equal sub-intervals of length $1/k$. $N$ random numbers\n",
    "    $u_{i}$ are generated and it is counted how many of these numbers\n",
    "    come to lie in each of the sub-intervals. If we call the number of\n",
    "    cases in each sub-interval $N_{i}$, $i$ = 1 \\... $k$, then the sum\n",
    "    \n",
    "    $$\n",
    "    \\chi^{2}=\\sum_{i=1}^{k}\\frac{(N_{i}-N/k)^{2}}{N/k}\n",
    "    $$ \n",
    "    \n",
    "    should (for $N / k \\geq 10)$ approximately follow a $\\chi^{2}$-distribution with\n",
    "    ($k$-1) degrees of freedom. This means that on average the ratio\n",
    "    $\\chi^{2} / (k-1)$ should be 1. Similar expressions can be\n",
    "    constructed for non-uniform distributions.\n",
    "\n",
    "-   **Test for correlation.** If $n$ successive random numbers are\n",
    "    plotted as coordinates in an $n$-dimensional space, then these\n",
    "    points lie on hyperplanes, as shown above. A good generator has many\n",
    "    hyperplanes which are uniformly distributed.\n",
    "\n",
    "```{margin}\n",
    "Another definition of gap test looks for the significance of the interval between recurrence of the same digit.\n",
    "```\n",
    "-   **Gap test.** Choose two numbers $\\alpha, \\beta$ with\n",
    "    $0 \\leq \\alpha < \\beta \\leq 1$. Generate $r+1$ random numbers, which\n",
    "    are uniformly distributed in the interval \\[0,1\\]. The probability\n",
    "    that the first $r$ numbers are not included in the interval\n",
    "    $[\\alpha, \\beta]$ and the last, $r+1^{st}$ number is included in\n",
    "    the interval should be\n",
    "    $P_{r}=p(1-p)^{r}\\,\\,\\,\\rm{with}\\,\\,\\, p=\\beta -\\alpha$.\n",
    "    \n",
    "-   **Random walk test** Choose a number $0 < \\alpha <1$. Build a large\n",
    "    set of random numbers and note the number of cases $r$ in which a\n",
    "    random number is smaller than $\\alpha$. We expect this to be a\n",
    "    binomial distribution for $r$ with $p = \\alpha$. This test is very\n",
    "    sensible for large values of $r$. The test should also be made for\n",
    "    the amount of random numbers which are larger than $(1-\\alpha)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arbitrarily distributed Random Numbers\n",
    "--------------------------------------\n",
    "\n",
    "Up to now we considered random numbers generated on a constant\n",
    "distribution. More generally we will need random numbers distributed\n",
    "according to some probability density $f(x)$. For example, we might need\n",
    "random numbers following a Gaussian distribution. In this section we\n",
    "will describe the most important methods to generated arbitrarily\n",
    "distributed random numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inverse of the Cumulative Distribution\n",
    "\n",
    "A standard procedure to produce random numbers generated according to\n",
    "the distribution $f(x)$ starts with random numbers $u_{i} \\in U(0,1)$\n",
    "and transforms them using the inverse function of the cumulative\n",
    "distribution $F(x)$:\n",
    "\n",
    "$$\n",
    "f(x)\\, dx = U(0,1)\\,du\\hspace{15mm}\\int_{-\\infty}^{x}f(t)\\,dt = F(x) = u\n",
    "\\hspace{15mm} x = F^{-1}(u)\n",
    "$$ \n",
    "\n",
    "$F^{-1}$ is the inverse function of the\n",
    "cumulative distribution function $F(x)$. The method is illustrated in\n",
    "{numref}`fig:-invFtrans`. For a sequence of uniform random numbers\n",
    "$u_{i}$, the random numbers $x_{i} = F^{-1}(u_{i})$ are distributed\n",
    "according to the probability density $f(x)$. In practice one matches the\n",
    "first x-% quantile of the first distribution to the x-% quantile of the\n",
    "target distribution $f(x)$.\n",
    "\n",
    "```{figure} ./images/ch4/invFtrans.png\n",
    "---\n",
    "width: 600px\n",
    "align: center\n",
    "name: fig-invFtrans\n",
    "---\n",
    "FIXME: Generation of random numbers on a continuous distribution $f(x)$\n",
    "using the inverse of the cumulative distribution function $F(x)$.\n",
    "```\n",
    "\n",
    "This method can only be applied if the integral of the probability\n",
    "density can be expressed as an analytic function $F(x)$ and this $F(x)$\n",
    "is invertible.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Random numbers generated on an exponential distribution.\n",
    "The p.d.f. for the exponential distribution is given by the\n",
    "$f(x;\\lambda) = \\lambda e^{-\\lambda x}$ for $x\\geq 0$, and it is zero\n",
    "for $x<0$. The cumulative distribution is:\n",
    "\n",
    "$$\n",
    "u= F(x;\\lambda) = \\int _{0}^{x}\\lambda e^{-\\lambda t} dt=1-e^{-\\lambda x}\n",
    "$$\n",
    "\n",
    "Inverting the cumulative we get the sequence of exponentially\n",
    "distributed random numbers $x_{i} = - \\ln (1 - u_{i} )/ \\lambda$. Or,\n",
    "because $u_{i}$ and $1-u_{i}$ are equally distributed in the interval\n",
    "$(0...1)$ we can write $x_{i} = - \\ln(u_{i}) / \\lambda$.\n",
    "```\n",
    "\n",
    "If we have an application with very large random numbers (for example\n",
    "very long lifetimes $t \\gg \\tau = 1 / \\lambda)$, then the above method\n",
    "might not be precise enough. Very large values of $x$ are generated by\n",
    "very small values of $u$. Because floating point numbers are represented\n",
    "with finite accuracy in a computer, large values of $x$ will appear as\n",
    "discrete.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Acceptance-Rejection Method\n",
    "\n",
    "The acceptance-rejection method (also known as \"hit-or-miss\"), even\n",
    "though not very efficient, can be used to generate random numbers\n",
    "according to a given probability density $f(x)$ when the cumulative\n",
    "distribution function $F(x)$ cannot be inverted.\\\n",
    "Under the assumption that the variable $x$ is restricted to some\n",
    "interval $a < x < b$, we can determine an upper limit $c$ with\n",
    "$c \\geq \\max(f(x))$; where $\\max(f(x))$ is the maximum of $f(x)$ in the\n",
    "interval $[a,b]$. This is then fed into the following algorithm:\n",
    "\n",
    "1.  Choose $x_{i}$ uniformly from the interval \\[a,b\\]:\n",
    "    $x_{i} = a + u_{i} \\cdot (b-a)$.\n",
    "\n",
    "2.  Choose another random number $u_{j} \\in U(0,1)$.\n",
    "\n",
    "3.  If $f(x_{i}) < u_{j} \\cdot c$, then go back to step 1, otherwise\n",
    "    accept $x_{i}$ as random number.\n",
    "\n",
    "The efficiency of this method is given by the ratio between the integral\n",
    "of $f(x)$ over \\[a,b\\] and the total area $c \\cdot (b-a)$ of the space\n",
    "of all generated pairs $(u_{i}, u_{j})$. The efficiency can be increased\n",
    "if we can find a function $s(x)$ which has the *approximate* shape of\n",
    "$f(x)$, which has an invertible cdf and for all $x$ in \\[a,b\\] and\n",
    "$s(x) > f(x)$. By means of\n",
    "$$\\int_{-\\infty}^{x}s(t)dt=S(x)\\hspace{20mm}x_{i}=S^{-1}(u_{i})$$ we can\n",
    "apply the following algorithm:\n",
    "\n",
    "1.  Choose a random number $u_{i}$ and calculate\n",
    "    $x_{i} = S^{-1}(u_{i})$.\n",
    "\n",
    "2.  Choose another random number $u_{j} \\in U(0,1)$.\n",
    "\n",
    "3.  If $f(x_{i}) \\leq u_{j} \\cdot s(x_{i})$, then go back to step 1,\n",
    "    otherwise accept $x_{i}$ as random number.\n",
    "\n",
    "The random number $x_{i}$ corresponds to the $s(x)$-distribution. The\n",
    "probability that it is accepted in step 3 is so increased because we\n",
    "better target the generation of the sampling points.\\\n",
    "Choosing an envelope for the generation of the sampling points is\n",
    "particularly useful in case of steeply falling or very \"irregular\"\n",
    "distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specially distributed Random Numbers\n",
    "\n",
    "A random angle $\\phi$, distributed in $[0, 2 \\pi]$, is obtained by\n",
    "$\\phi_{i} = 2 \\pi \\cdot u_{i}$. The random unit vector is then simply:\n",
    "\n",
    "$$\n",
    "\\cos \\phi_{i} \\choose \\sin \\phi_{i}\n",
    "$$\n",
    "\n",
    "In three dimensions, an additional polar angle\n",
    "$\\theta \\in [- \\pi / 2, \\pi / 2]$ is needed (in addition to sin $\\phi$\n",
    "and cos $\\phi$). According to the differential solid angle\n",
    "\n",
    "$$\n",
    "d\\Omega = {\\rm sin}\\theta\\, d\\theta\\, d\\phi = |d\\,{\\rm cos}\\theta|\\, d\\phi\n",
    "$$\n",
    "\n",
    "we obtain $\\theta_{j} = {\\rm arcsin} (2 \\cdot u_{j} -1)$ from the\n",
    "analytical transformation. For the 3 components of the random unit\n",
    "vector we get $e_{x} = \\sin \\phi \\cdot \\cos \\theta$,\n",
    "$e_{y} = \\cos \\phi \\cdot \\cos \\theta$ and $e_{z} = \\sin \\theta$.\n",
    "\n",
    "This is one of the most commonly used distribution for random numbers.\n",
    "The simplest but only approximately correct generator for the random\n",
    "numbers $z_{i}$, which are distributed according to a Gaussian\n",
    "distribution ($1 / \\sqrt{2 \\pi} \\cdot e^{-x^{2}/2}$) in the interval\n",
    "\\[-6,6\\], is based on the central limit theorem {numref}`fig-CLT`:\n",
    "\n",
    "$$\n",
    "z_{i}=\\sum_{j=1}^{12}u_{j}-6\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random numbers in ROOT and PYTHON\n",
    "\n",
    "ROOT (TRandom3) and PYTHON (import random) use the Mersenne Twister\n",
    "algorithm wiht a period of $2^{19927}-1 \\sim 4.3 \\cdot 4.3 10^{6001}$\n",
    "equidistributed in up to 623 dimensions (for 32-bit machines). Both\n",
    "ROOT and PYTHON have predefined functions to generate numbers\n",
    "according to the most common pdf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo Integration\n",
    "-----------------------\n",
    "\n",
    "Monte Carlo techniques can be used to evaluate (definite) integrals:\n",
    "\n",
    "$$\n",
    "I = \\int_a^b f(x) dx\n",
    "$$ \n",
    "\n",
    "with $f(x)$ is any function. The\n",
    "straightforward deterministic way to evaluate the integral is to divide\n",
    "the range $[a,b]$ in n intervals and compute $I$ as: \n",
    "\n",
    "```{math}\n",
    ":label:eq:Idet\n",
    "\\frac{b-a}{n}\\sum_{i=1}^n f(x_i)\n",
    "```\n",
    "\n",
    "where $x_i = a+(i-0.5)(b-a)$. This is\n",
    "what is typically called the method of the *trapezoid*. You can evaluate\n",
    "the funciton at the central value of the bin, or at the min and max and\n",
    "approximate the function with a linear interpolation, etc \\...\n",
    "\n",
    "```{figure} ./images/ch4/trapezoid.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-trapezoid\n",
    "---\n",
    "FIXME Ref: Definite integral computed with the trapezoid method [wiki].\n",
    "```\n",
    "\n",
    "The uncertainty on the integral goes as $n^{-2}$ in one dimension.\n",
    "\n",
    "The most naive implementation of the Monte Carlo way of computing the\n",
    "integral differs in the choice of the points where we evaluate the sum:\n",
    "instead of regularly spaced points we randomly sample the axis:\n",
    "\n",
    "$$\n",
    "x_i = a+ r_i (b-a)\n",
    "$$ \n",
    "with $r_i$ random numbers uniformly distributed\n",
    "in $[0,1]$.\n",
    "\n",
    "The first consequence of this approach is that the result is non\n",
    "deterministic: different sequences of random numbers will give different\n",
    "integral values.\n",
    "\n",
    "The uncertainty on the estimate of the integral $I$ can be estimated as\n",
    "the variance of $f(x_{i})$ in the following way:\n",
    "\n",
    "$$\n",
    "V[I_{MC}]=\\sigma_{I_{MC}}^{2}=V\\Big[\\frac{b-a}{n}\\sum_{i=1}^{n}f(x_{i})\\Big]=\n",
    "\\Big(\\frac{b-a}{n}\\Big)^{2}V\\Big[\\sum_{i=1}^{n}f(x_{i})\\Big]=\\frac{(b-a)^{2}}{n}V[f(x_{i})]\n",
    "$$\n",
    "\n",
    "This equation shows us that the accuracy of the computation decreases\n",
    "with $1/\\sqrt{n}$.\n",
    "\n",
    "To compute integrals in larger dimensions with the trapezoidal method\n",
    "you need to create an d-dimensional equidistand grid of points, evaluate\n",
    "the function at those points and sum over the d-dimensional trapezoid.\n",
    "The uncertainty of the trapezoidal method in d-dimensions goes as\n",
    "$n^{-2/d}$, while the uncertainty of the Monte Carlo approach is\n",
    "$1/\\sqrt{n}$ in any dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methods to reduce the Variance\n",
    "\n",
    "The Monte Carlo integration shown above is just the most naive\n",
    "implementation of this technique. There are several ways to improve the\n",
    "numerical result (reducing the variance)\n",
    "\n",
    "Simply dividing the range of integration in two regions and generating\n",
    "half of the Monte Carlo points in each of the regions reduces the\n",
    "variance. The reason is that in this way we allow a more uniform\n",
    "sampling of the distribution.\n",
    "\n",
    "Given that we know how the integrand $f(x)$ behaves, we can sample the\n",
    "distribution with higher density of random points where the function is\n",
    "varying more rapidly. It's like in the stratification but this time the\n",
    "intervals are chosen in a more cleverer way.\n",
    "\n",
    "Because the variance of the Monte Carlo result is proportional to the\n",
    "variance of the integrand, it makes sense to transform the integral to\n",
    "get a new integral with a smaller variance than the original one. By\n",
    "introducing a function $g(x)$ we can write: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\int_{a}^{b} f(x)dx=\\int_{a}^{b} \\Big[\\frac{f(x)}{g(x)}\\Big]g(x)dx =\n",
    "\\int_{A}^{B} \\Big[\\frac{f(x)}{g(x)}\\Big]dv(x)\n",
    "\\quad \\text{with }v(x)=\\int g(x)dx \\end{aligned}\n",
    "$$ \n",
    "\n",
    "The variance of the\n",
    "new result is now proportional to the one of $f(x)/g(x)$ instead of\n",
    "being proportional to $f(x)$ alone. If $g(x)$ has been chosen carefully,\n",
    "the method of importance sampling allows us to reduce the variance of\n",
    "the Monte Carlo integration considerably. But the function $g(x)$ has to\n",
    "be integrable and invertible, and it has to describe the original\n",
    "function $f(x)$ appropriately enough.\n",
    "\n",
    "This is similar to importance sampling except that instead of dividing\n",
    "by g(x) subtract it:\n",
    "\n",
    "$$\n",
    "I = \\int f(x) dx = \\int [f(x) - g(x)] dx + \\int g(x) dx\n",
    "$$ \n",
    "\n",
    "Here $\\int g(x) dx$ must be known and $g$ is chosen such that $f-g$ has a\n",
    "smaller variance than $f$. This method does not risk the instability of\n",
    "importance sampling, nor is it necessary to invert the integral of\n",
    "$g(x)$.\n",
    "\n",
    "Up to now the Monte Carlo points were all independent. But looking at\n",
    "the variance of two general functions:\n",
    "\n",
    "$$\n",
    "V[f_1(x)+f_2(x)] = V[f_1(x)] + V[f_2(x)] + 2cov[f_1(x),f_2(x)]\n",
    "$$ \n",
    "\n",
    "and writing the integral as $I = \\int_a^b f = \\int_a^b(f_1 + f_2) dx$ we\n",
    "observe that we can reduce the covariance of $I$ by introducing a large\n",
    "negative correlation.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Suppose that we know that $f(x)$ is a monotonically\n",
    "increasing function of $x$. Then let $f_1 = 1/2 f(x)$ and\n",
    "$f_2 = 1/2 f(b-(x-a))$. Clearly the integral of $(f_1 +f_2)$ is just the\n",
    "integral of $f$. However, since $f$ is monotonically increasing, $f_1$\n",
    "and $f_2$ are negatively correlated; when $f_1$ is small, $f_2$ is large\n",
    "and vice versa. If this negative correlation is large enough,\n",
    "$V [f_1 + f_2] < V [f]$.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo Simulations\n",
    "-----------------------\n",
    "\n",
    "We will briefly describe in this section the typical applications of\n",
    "Monte Carlo techniques in high energy physics: simulation of a physics\n",
    "process, simulation of a detector, study of reconstruction programs,\n",
    "physics reach of an experiment, \"toy\" experiments.\n",
    "\n",
    "### Simulation of a physics process\n",
    "\n",
    "In particle physics, \"generators\" are programs used to simulate the\n",
    "production (and decays) of particles. The generators allow to compute\n",
    "the full kinematics of the process and to get the distributions of the\n",
    "observables, which are typically impossible to derive with analytical\n",
    "methods. The description of a proton proton collision for example\n",
    "requires the understanding of several components: the hard scatter\n",
    "(matrix elements + pdf), initial/final state radiation, parton showering\n",
    "and hadronization and decays of the generated final states. This kind of\n",
    "process can be treated only with Monte Carlo simulations. An\n",
    "introduction to Monte Carlo generators in high energy physics can be\n",
    "found in [@MCgen].\n",
    "\n",
    "### Simulation of a detector\n",
    "\n",
    "Once the particles are generated we will need to understand how they\n",
    "will be \"seen\" by a detector. This requires basically two steps, which\n",
    "in jargon are called \"simulation\" and \"digitization\". In the first the\n",
    "particles are propagated through the detector simulating their\n",
    "interactions with all sensitive and passive materials. For instance a\n",
    "high energy photon hitting a calorimeter will create an electromagnetic\n",
    "shower, while a charged particle will deposit a fraction of its energy\n",
    "in a silicon layer. These phenomena are simulated with standard programs\n",
    "like GEANT or FLUKA, MARS etc\\... which contain software libraries\n",
    "describing the details of the interactions. The digitization step\n",
    "simulates how the energy deposited in the different sensitive material\n",
    "is transformed into electric/optical signals. In this case the signals\n",
    "are a mixture of deterministic effects (e.g. amplification/shaping)\n",
    "overlapped with random effects (e.g. noise or finite accuracy of\n",
    "calibrations) described with Monte Carlo methods.\n",
    "\n",
    "### Study of reconstruction programs\n",
    "\n",
    "The electronic signals can then be analysed to \"reconstruct\" physics\n",
    "objects. Those can be e.g. the tracks of charged particles in a tracker,\n",
    "or the electromagnetic showers in a calorimeter. Even in this case,\n",
    "Monte Carlo simulations are essential tools to develop the\n",
    "reconstruction software. In this case, starting from Monte Carlo samples\n",
    "of particles, and so knowing exactly their type and their kinematics, we\n",
    "can develop the reconstruction packages to be able to reconstruct them\n",
    "with high efficiency and accuracy.\n",
    "\n",
    "### Physics reach of an experiment\n",
    "\n",
    "Simulations are not important only to decipher the results of an\n",
    "experiment once the data are collected. They are the first necessary\n",
    "step in the design of an experiment. The conception of a detector starts\n",
    "with simulations where different solutions can be tested without doing\n",
    "the actual prototypes, which would be unacceptably expensive and time\n",
    "consuming.\n",
    "\n",
    "### Toy modelling\n",
    "\n",
    "\"A good physicist knows how to build toy models.\" This quotation from an\n",
    "anonymous physicist summarize the importance of toy modeling in physics.\n",
    "Toys can be used to get a first understanding of a process picking only\n",
    "its essential characteristics (here is where the physicist has to be\n",
    "good) and avoiding all complications of real experimental life. As it\n",
    "will be seen in hypothesis testing, toys are also an ubiquitous tool\n",
    "used in statistics for hypothesis testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References\n",
    "----------\n",
    "\n",
    "Most of the material of this section comes from:\n",
    "\n",
    "-   L. Lyons [@Lyons], \"Statistics for Nuclear and Particle Physicist\":\n",
    "    Ch. 6\n",
    "\n",
    "- FIXME [@KnuthBook]\n",
    "\n",
    "- FIXME [@MCgen]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}