mva.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multivariate Analysis Methods\n",
    "=============================\n",
    "\n",
    "In this chapter we will describe some multi-variate analysis (MVA)\n",
    "methods that are very frequently used in particle physics. To put things\n",
    "in context the methods we will discuss belongs to the much wider area of\n",
    "Artificial Intelligence (AI). AI studies the systems that perceive the\n",
    "environment. Within this huge field, machine learning (ML) is the\n",
    "technology of getting computers to act without being explicitly\n",
    "programmed to do so. Typically we then distinguish supervised learning\n",
    "from unsupervised learning. In supervised learning you instruct an\n",
    "algorithm by examples: data are presented to the algorithm with a tag\n",
    "and the algorithm learns how to associate data to the different tags by\n",
    "analysing their characteristics. In unsupervised learning (a.k.a.\n",
    "clustering algorithms) the algorithm will discover by itself the\n",
    "different tags (populations) in data by looking at their\n",
    "characteristics.\\\n",
    "In today's HEP, supervised learning is by far the most used type of\n",
    "learning. These algorithms are used to solve two classes of problems:\n",
    "classification and regression. In a classification problem the goal is\n",
    "to subdivide the elements of a dataset into a discrete set of classes\n",
    "depending on their characteristics. Typical examples are\n",
    "signal/background separation or particle identification: photon/jets,\n",
    "etc\\... A regression problem is conceptually identical but the classes\n",
    "instead of belonging to a discrete set are represented by a continuum.\n",
    "Most commonly is used to improve energy measurements resolution. The\n",
    "reconstructed energy of a jet is affected by several detector effects\n",
    "which depend on the position of the jet in the detector, its shape, the\n",
    "reconstructed energy, etc\\... The regression algorithm assigns the jet\n",
    "to a continuous value, its \"regressed energy\". It does it by looking at\n",
    "the examples he has been trained on and recognizing which (true) energy\n",
    "is closer to the case at hand.\\\n",
    "The key point of these algorithms is that they learn from a training set\n",
    "of data and then they use the acquired knowledge to take decisions on\n",
    "data they have never seen before.\\\n",
    "MVA techniques are coded in several packages. The most used in HEP is\n",
    "`TMVA` [@TMVA] that comes with `ROOT`. Other packages often used are\n",
    "`scikit` [@scikit] in python and `R` [@R].\n",
    "\n",
    "Definitions\n",
    "-----------\n",
    "\n",
    "In this section we will set the language that we will use to study the\n",
    "MVAs. We will build on the concept of test statistics developed in the\n",
    "previous chapters and so we conveniently use the language of statistics\n",
    "(we could have chosen to use the ML - computer science language, it's\n",
    "just a different naming, the concepts are the same).\\\n",
    "Let's take as a working example a classification problem with only two\n",
    "classes (we will address the regression case later). In the language of\n",
    "statistics, implementing a classifier means to choose a decision\n",
    "boundary that allows to separate the two classes. To do so we\n",
    "characterize each event with a number of variables\n",
    "$\\vec{x} = (x_1,x_2, \\ldots, x_n)$ and define the decision boundary as a\n",
    "hyper-surface in this $n$-dimensional space $y(\\vec{x})=c$. The function\n",
    "$y(\\vec{x})$ is our test statistic which compresses the full information\n",
    "contained in the $n$ variables into one number.\\\n",
    "\\\n",
    "**Example**  As an example imagine you have to separate tracks into muon\n",
    "and electron candidates. The variables $\\vec{x}$ could be $p_T$, $\\eta$,\n",
    "$\\phi$, some PID on the track, number of hits in the muon chamber,\n",
    "presence of an electromagnetic cluster in the same direction of the\n",
    "track, etc\\... The muon candidates will preferentially populate a\n",
    "certain region of this multidimensional space (e.g. large number of hits\n",
    "in the muon chambers) while the electrons a different one (energy\n",
    "deposits in the electromagnetic calorimeter). Formally the distribution\n",
    "of $\\vec{x}$ will follow some n-dimensional joint p.d.f. depending on\n",
    "the hypothesis (muon/electron): pdf$(\\vec{x}|H_i), \\; i= \\mu, e$.\\\n",
    "\\\n",
    "The choice of the boundary depends on several factors. It depends on the\n",
    "variables you choose (physics driven), on the type of classifier you\n",
    "want to use, on computational issues (CPU, memory) and typically it\n",
    "boils down to finding a compromise between performance and complexity of\n",
    "the classifier.\\\n",
    "\\\n",
    "**Example**  Looking at\n",
    "Fig. [1.1](#fig:decisionBoundary){reference-type=\"ref\"\n",
    "reference=\"fig:decisionBoundary\"} we have two variables and two classes\n",
    "\"red\" and \"blue\" distributed according to different pdfs. The different\n",
    "pictures show different choices of decision boundary: (left) just a cut\n",
    "on each of the variables , (middle) a linear combination of the two\n",
    "variables (right) a non linear combination of the two\n",
    "variables [@CowanMVA].\\\n",
    "\\\n",
    "\n",
    "![Different choices of decision\n",
    "boundary.[\\[fig:decisionBoundary\\]]{#fig:decisionBoundary\n",
    "label=\"fig:decisionBoundary\"}](Section10Bilder/decisionBoundary.png){#fig:decisionBoundary\n",
    "width=\"100%\"}\n",
    "\n",
    "Build the decision boundary\n",
    "---------------------------\n",
    "\n",
    "For the present discussion on classification, signal and background are\n",
    "any two generic classes (photons/jets, jets/b-jets, cats/dogs), for\n",
    "simplicity we will keep calling them signal and background throughout\n",
    "this section.\\\n",
    "We've seen in\n",
    "Sec.[\\[SectionNeymanPearson\\]](#SectionNeymanPearson){reference-type=\"ref\"\n",
    "reference=\"SectionNeymanPearson\"} that the optimal way to set the\n",
    "decision boundary having the highest background rejection for a given\n",
    "signal efficiency for a simple hypothesis is to use the likelihood\n",
    "ratio: $$\\label{eq:LRNP}\n",
    " y(\\vec{x}) = \\frac{p(\\vec{x}|s)}{p(\\vec{x}|b)}.$$ Unfortunately, in any\n",
    "non trivial practical case, we don't have an analytic expressions for\n",
    "the pdf of the two classes and typically we recur to Monte Carlo\n",
    "techniques to model the pdfs. In this section we will see a few methods\n",
    "to build the pdf from Monte Carlo and see why in practice they are not\n",
    "used for problems with a large number of dimensions.\n",
    "\n",
    "### Histogramming\n",
    "\n",
    "The easiest way to build the pdf for signal and background from a Monte\n",
    "Carlo sample is to fill two $n-$dimensional histograms. This method has\n",
    "the advantage of being trivial to implement, and being computationally\n",
    "very light: the information is binned into the histogram once and then\n",
    "the original dataset can be discarded. The drawback of this approach is\n",
    "the generic problem of all histograms: need to choose a proper binning\n",
    "(not to coarse, not too fine) and the unavoidable discontinuities at the\n",
    "bin boundaries.\n",
    "\n",
    "### Kernel Density Estimators (KDE) and K-Nearest Neighbors (KNN)\n",
    "\n",
    "With these methods we build the decision boundary \"event-by-event\"\n",
    "judging if an event of the n-dimensional space of the variables is to be\n",
    "assigned to the signal region or to the background region. Intuitively,\n",
    "given a point $\\vec{x}$ in the n-dimensional space of the variables we\n",
    "count how many signal and background events are present in a local\n",
    "region around $\\vec{x}$ and then we apply a simple majority vote: that\n",
    "portion of the space is assigned to signal or background depending on\n",
    "which of the two has the largest population. The use of the concept of\n",
    "locality requires the definition of a metric in the $n$-dimensional\n",
    "space. The choice of the metric requires special consideration because\n",
    "the different dimensions might have different units and scales\n",
    "($E\\in[10,500]$ GeV, $\\eta\\in[-2.5,2.5]$, $PID\\in[0,10]$ MeV/cm,\n",
    "etc\\...). Often the variables are rescaled or remapped to have\n",
    "comparable numerical values in all dimensions.\\\n",
    "Both Kernel Density Estimators (KDE) and K-Nearest Neighbors (KNN) are\n",
    "built around the estimator: $$p(\\vec{x}) = \\frac{K}{NV}$$ where\n",
    "$\\vec{x}$ is the point of the space we're sampling, N is the total\n",
    "number of events in the sample and $K$ is the number of events in the\n",
    "hyper-volume $V$.\\\n",
    "To optimize the classifier performance you can:\n",
    "\n",
    "-   fix $K$ and determine $V$ $\\to$ KNN\n",
    "\n",
    "-   fix $V$ and determine $K$ $\\to$ KDE\n",
    "\n",
    "For the KNN you fix the number of events $K$ (a.k.a. smoothing\n",
    "parameter) you want to have in your region of the space and you increase\n",
    "the volume $V$ until it contains such number. Then you perform the\n",
    "majority vote: you count how many events are of type signal, how many\n",
    "are of type background and the class that has the largest number defines\n",
    "whether that portion of space is to be assigned to signal or\n",
    "background.\\\n",
    "For the KDE you fix the hyper-volume $V$ and change the number of events\n",
    "K. The \"shape\" for the hyper-volume $V$ is typically chosen to be a\n",
    "gaussian in $n-$dimensions. This allows to have a smooth function for\n",
    "the model. (An hyper-cube would do but it would introduce\n",
    "discontinuities at the edges). In practice you place a gaussian kernel\n",
    "of standard deviation $h$ centred about each data point, then for a\n",
    "given $\\vec{x}$ you add up the contribution from all the gaussians and\n",
    "you normalize (divide by N):\n",
    "$$p(\\vec{x}) =\\frac{1}{N}\\sum_{i=1}^{N} \\frac{1}{\\sqrt{2\\pi h^2}}\\exp\\left(\\frac{-|\\vec{x} - \\vec{x_i}|^2}{2h^2}\\right)$$\n",
    "\n",
    "### Training, bias and variance\n",
    "\n",
    "The classification methods shown so far use a dataset with labelled data\n",
    "(signal or background) to build the decision boundary. We will then use\n",
    "this boundary to classify new data. This procedure is called \"supervised\n",
    "learning\"; the building of the boundary is called **training** step\n",
    "while the use of the boundary is usually called **application** (or\n",
    "simply classification) step. The important point to notice is that the\n",
    "classifier is tuned on a sample of labelled data representative of the\n",
    "parent distribution of the classes, but it will be applied on new data\n",
    "coming from a separate sampling. Clearly, even if the training sample\n",
    "and the new data come from a sampling of the same distribution, they\n",
    "will be affected by statistical fluctuations. The training phase can\n",
    "reduce the mis-classification rate to zero on the training sample\n",
    "itself, but it will have poor performance when applied to new data.\\\n",
    "The training has to be optimized minimizing two competing effects, the\n",
    "\"bias-variance\" tradeoff:\n",
    "\n",
    "-   *bias* (mis-classification): the fraction of signal events ending up\n",
    "    in the background region (or viceversa)\n",
    "\n",
    "-   *variance* (over-training) limitation in the generalization of the\n",
    "    classifier to different data samples\n",
    "\n",
    "In general all methods have one or more parameters that can be tuned to\n",
    "find the optimal classifier. The KNN and the KDE methods have a\n",
    "\"smoothing parameter\": $K$ for the KNN method and $h$ for the KDE. By\n",
    "tuning this parameter we can optimize our classifier. Having a large\n",
    "value of the smoothing parameter will produce well populated regions,\n",
    "reducing the statistical fluctuation, but generally increasing the bias;\n",
    "viceversa you can reduce the value of the smoothing parameter bringing\n",
    "to zero the bias but increasing the variance, i.e. making the training\n",
    "very sensitive to the specific sample used for training and limiting its\n",
    "generalization.\n",
    "\n",
    "### Curse of dimensionality and learning algorithms\n",
    "\n",
    "All methods described above (histogramming, KDE and KNN) become\n",
    "unmanageable when using a large number of input variables (dimensions).\n",
    "This problem goes under the name of \"the curse of dimensionality\". To\n",
    "understand the issue, suppose your data are uniformly distributed in a\n",
    "$D$-dimensional unit cube. The methods above will all try to catch a\n",
    "fraction $r$ of events in a small portion of space. Let's suppose this\n",
    "portion of space has the shape of a hyper-cube with side $s$ (its volume\n",
    "being $s^D$).\\\n",
    "The fraction of the volume taken by this cube is $r = s^D /1$ (the total\n",
    "space is a unit cube) and its side is simply $s = r^{1/D}$. If you want\n",
    "a sampling fraction $r=0.001$ (one per mill of the whole space) you will\n",
    "need different sizes $s$ depending on the number of dimensions. For\n",
    "$D=1$ $s = 0.001$, for $D=3$ $s = 0.1$ which is 10% of the whole space,\n",
    "for $D=30$ $s=0.8$ which is 80% of the whole space. To properly build\n",
    "the estimator, you will need to have a very large training set to\n",
    "adequately populate all corners of the $D$-dimensional space\\\n",
    "Learning algorithms will provide better ways to sample the space.\n",
    "\n",
    "Fisher discriminant\n",
    "-------------------\n",
    "\n",
    "The Fisher discriminant approximates the likelihood ratio in\n",
    "Eq. [\\[eq:LRNP\\]](#eq:LRNP){reference-type=\"ref\" reference=\"eq:LRNP\"}\n",
    "with a linear combination of the input dataset.\n",
    "$$y(\\vec{x}) = \\frac{p(\\vec{x}|s)}{p(\\vec{x}|b)} = \\sum_{i=1}^N w_ix_i = \\vec{w}^T\\vec{x}$$\n",
    "The decision boundary will be a constant in the one dimensional case, a\n",
    "straight line in a 2-dimensional case and in general a $n-1$ hyperplane\n",
    "in an $n$-dimensional problem.\\\n",
    "In the case of a Fisher discriminant the training phase consists in\n",
    "finding the best set of weights $\\vec{w}$ which maximize the signal to\n",
    "background (Fisher) separation. The separation is defined as:\n",
    "$$J(\\vec{w}) = \\frac{(\\tau_s - \\tau_b)^2}{\\Sigma_s^2 + \\Sigma_b^2}$$\n",
    "where $\\tau_s, \\tau_b$ and $\\Sigma_s, \\Sigma_b$ are the mean and width\n",
    "(covariance) of the signal and background (see\n",
    "Fig. [1.2](#fig:Fisher1){reference-type=\"ref\" reference=\"fig:Fisher1\"}).\n",
    "\n",
    "![Example to fix the notation in a mono-dimensional\n",
    "case.[\\[fig:Fisher1\\]]{#fig:Fisher1\n",
    "label=\"fig:Fisher1\"}](Section10Bilder/Fisher1.png){#fig:Fisher1\n",
    "width=\"60%\"}\n",
    "\n",
    "The training translates into writing the means and the covariance\n",
    "matrices as a linear combination of the variables and then find the\n",
    "weights that maximize the separation $J(\\vec{w})$. It is useful to\n",
    "notice that the numerator of $J$ is the separation between the classes\n",
    "(the distance between the means) and the denominator is the separation\n",
    "within each class (the spread of the variables).\\\n",
    "The means and covariances can be written as: $$\\begin{aligned}\n",
    "(\\mu_k)_i &=& \\int x_i p(\\vec{x}|H_k) d\\vec{x} \\\\\n",
    "(V_k)_{ij} &=& \\int (x-\\mu_k)_i(x-\\mu_k)_j p(\\vec{x}|H_k) d\\vec{x}\\end{aligned}$$\n",
    "where k = \"signal\" or \"background\" and $i,j=1,...,n$ are the components\n",
    "of $\\vec{x}$. From this we can compute the mean and variance of\n",
    "$y(\\vec{x})$ $$\\begin{aligned}\n",
    "\\mu_k &=& \\int y(\\vec{x}) p(\\vec{x}|H_k) d\\vec{x} = \\vec{w}^T \\vec{\\mu}_k\\\\\n",
    "\\Sigma_k^2 &=& \\int (y(\\vec{x})-\\tau_k)^2  p(\\vec{x}|H_k) d\\vec{x} = \\vec{w}^T V_k \\vec{w}\\end{aligned}$$\n",
    "The numerator of $J(\\vec{w})$, i.e. the separation between classes, then\n",
    "becomes $$\\begin{aligned}\n",
    "        (\\tau_s - \\tau_b)^2 &=& \\sum_{i,j=1}^N  w_i w_j (\\mu_s-\\mu_b)_i(\\mu_s-\\mu_b)_j\\\\\n",
    "                           &=& \\sum_{i,j=1}^N w_i w_j B_{ij} = \\vec{w}^TB\\vec{w}\\end{aligned}$$\n",
    "and the denominator\n",
    "$$\\Sigma_s^2 + \\Sigma_b^2 = \\sum_{i,j}^N w_i w_j (V_s + V_b)_{ij} = \\vec{w}^T W \\vec{w}.$$\n",
    "At this point we just need to maximize the ratio\n",
    "$$J(\\vec{w}) = \\frac{\\vec{w}^TB\\vec{w}}{\\vec{w}^T W \\vec{w}}$$ by\n",
    "solving $\\partial J  / \\partial w_i = 0$. With the obtained weights we\n",
    "can then build the Fisher discriminant $y(\\vec{x}) = \\vec{w}^T \\vec{x}$.\n",
    "\n",
    "![Example in 2-dimensions (i.e. 2 variables). Signal is represented by\n",
    "orange circles, background with blue ones. The red straight line\n",
    "represents the Fisher discriminant which would be a linear combination\n",
    "of the variables $x_1$ and $x_2$. The orange and blue distributions\n",
    "visually show the separation between the projections of the two\n",
    "populations.[\\[fig:Fisher2\\]]{#fig:Fisher2\n",
    "label=\"fig:Fisher2\"}](Section10Bilder/Fisher2.png){#fig:Fisher2\n",
    "width=\"40%\"}\n",
    "\n",
    "The Fisher discriminant provides a linear decision boundary (see\n",
    "Fig. [1.3](#fig:Fisher2){reference-type=\"ref\" reference=\"fig:Fisher2\"}).\n",
    "If the classes you want to separate show some particular non linear\n",
    "structure like the one in Fig. [1.4](#fig:Fisher3){reference-type=\"ref\"\n",
    "reference=\"fig:Fisher3\"} then you can still use a Fisher discriminant\n",
    "after a suitable remapping of the variables. Very often the non linear\n",
    "structure (especially in a multidimensional space) is not evident and\n",
    "for this you can use other MVA tools like BDT and ANN (see next\n",
    "chapters).\n",
    "\n",
    "![The two classes show a clear symmetry: remapping them allows to still\n",
    "use a linear discriminant.[\\[fig:Fisher3\\]]{#fig:Fisher3\n",
    "label=\"fig:Fisher3\"}](Section10Bilder/Fisher3.png){#fig:Fisher3\n",
    "width=\"80%\"}\n",
    "\n",
    "Decision trees\n",
    "--------------\n",
    "\n",
    "Because of the importance gained in recent years we will describe in\n",
    "some detail how to apply decision trees to solve classification and\n",
    "regression problems.\n",
    "\n",
    "### Classification\n",
    "\n",
    "A decision tree is a (binary) tree used to partition the variables'\n",
    "space into rectangles and then by majority vote assign each rectangle to\n",
    "a class. To understand how it works we take the example in\n",
    "Fig. [1.5](#fig:BDT1){reference-type=\"ref\" reference=\"fig:BDT1\"}.\n",
    "\n",
    "![Classification example.[\\[fig:BDT1\\]]{#fig:BDT1\n",
    "label=\"fig:BDT1\"}](Section10Bilder/BDT1.png){#fig:BDT1 width=\"100%\"}\n",
    "\n",
    "Let's consider a two dimensional space ($x^1$, $x^2$) and consider the\n",
    "usual two classes \"signal\" and \"background\" here shown in red and blue\n",
    "respectively as in Fig. [1.5](#fig:BDT1){reference-type=\"ref\"\n",
    "reference=\"fig:BDT1\"}-(a). In this example the variables can assume\n",
    "continuous values, but in reality can be any discrete or even non\n",
    "numerical value (an example can be a loose/tight selection criterion).\n",
    "To grow a decision tree means to place cuts (binary choices of the type\n",
    "pass/fail) to reach the minimal signal/background misclassification at\n",
    "each step.\\\n",
    "The first step is to decide which variable to cut on: we choose the\n",
    "variable that provides the greatest increase in the separation between\n",
    "the two classes in the two daughters node relative to the parent. To\n",
    "quantify the separation we typically use as a metric the \"Gini index\"\n",
    "defined as\n",
    "$$\\mbox{Gini}=P(1-P) \\qquad \\mbox{with} \\qquad P=\\frac{\\sum_{\\mbox{signal}}w_i}{\\sum_{\\mbox{signal}}w_i + \\sum_{\\mbox{background}}w_i}$$\n",
    "where the $w_i$ are the number of signal/background events in each node\n",
    "(or in weights in case of weighted events). The minimal\n",
    "misclassification (maximal separation) is reached when $P=0$ or $P=1$\n",
    "(selecting all signal events is equivalent to select all background\n",
    "events), while you get maximal misclassification (random guess, minimal\n",
    "separation) when $P=0.5$. Applied to the example in\n",
    "Fig. [1.5](#fig:BDT1){reference-type=\"ref\" reference=\"fig:BDT1\"} we will\n",
    "scan the two dimensions and look for the cut that maximizes the\n",
    "separation. In this case we select the variable $x^2$ and a cut value of\n",
    "1.5 (b). We then repeat the procedure and scan again the two dimensions\n",
    "to find the cut that minimizes the misclassification, this time it lead\n",
    "a value of 2.0 on $x^1$. Every time we repeat the procedure we have to\n",
    "choose both the dimension and the value of the cut. It can happen as in\n",
    "(d) that the same dimension is selected twice in a raw this time cutting\n",
    "at 1.5. The example then continues with (e) select $x^1$ and cut at 2.1\n",
    "and (f) again selecting $x^1$ and cutting at 1.6. The procedure will\n",
    "continue until we reach a minimum number of points in each of the\n",
    "rectangles. As we will see later in this section the minimum number of\n",
    "points gives a handle to limit the over-training. From the set of\n",
    "rectangles in Fig. [1.5](#fig:BDT1){reference-type=\"ref\"\n",
    "reference=\"fig:BDT1\"} we can build a binary decision tree as in\n",
    "Fig. [1.6](#fig:BDT2){reference-type=\"ref\" reference=\"fig:BDT2\"}\n",
    "\n",
    "![Decision tree corresponding to the classification example in\n",
    "Fig. [1.5](#fig:BDT1){reference-type=\"ref\"\n",
    "reference=\"fig:BDT1\"}.[\\[fig:BDT2\\]]{#fig:BDT2\n",
    "label=\"fig:BDT2\"}](Section10Bilder/BDT2.png){#fig:BDT2 width=\"50%\"}\n",
    "\n",
    "The last layer of the binary tree are called \"leaves\" (pictured as a\n",
    "circle) and they contain a certain number of signal and background\n",
    "events. In Fig. [1.6](#fig:BDT2){reference-type=\"ref\"\n",
    "reference=\"fig:BDT2\"} they are shown as colored numbers: red for signal\n",
    "and blue for background. Now we have to choose how to assign each leaf\n",
    "to either signal or background. We do this by a majority vote. The class\n",
    "with the largest population defines if that leaf (or rectangle in the\n",
    "previous schematics) is associated with one class or the other. Again in\n",
    "the same figure the leaves are colored in red and blue according to the\n",
    "chosen class. The procedure is equivalent to writing a piece-wise\n",
    "constant function over the plane.\\\n",
    "The described procedure is generally called \"training\": we use labelled\n",
    "data (i.e. we know what is signal and what is background, as we would\n",
    "have with a Monte Carlo sample) and we build the tree. This is the point\n",
    "where the algorithm learns how to classify the data. Once this is done,\n",
    "we can apply the decision tree to a new dataset (never seen before by\n",
    "the algorithm) and use it to classify new elements. For example a point\n",
    "$(1.7, 1.8)$ in the previous case, would be classified as signal because\n",
    "it will land in the leaf/rectangle with signal=3 and background=1.\\\n",
    "\n",
    "### Regression\n",
    "\n",
    "A regression problem is conceptually equivalent to a classification one\n",
    "but the target variable, instead of being discrete is represented by a\n",
    "continuous value. In the example in\n",
    "Fig. [1.7](#fig:BDTregr1){reference-type=\"ref\"\n",
    "reference=\"fig:BDTregr1\"}-(a) we consider a training sample on the\n",
    "$\\mathbb{R}$ plane $(x,y)$.\n",
    "\n",
    "![Regression example.[\\[fig:BDTregr1\\]]{#fig:BDTregr1\n",
    "label=\"fig:BDTregr1\"}](Section10Bilder/BDTregr1.png){#fig:BDTregr1\n",
    "width=\"100%\"}\n",
    "\n",
    "The idea is again the same as the one used for the classification: given\n",
    "a set of variables assign the correct category. Here the only difference\n",
    "is that the category are defined on a continuous sample. In this example\n",
    "the variable is $x$ while $y$ represent the continuous category\n",
    "(target). The strategy remains to set cuts such that the obtained\n",
    "partitions contain elements with similar characteristics (in\n",
    "classification we were grouping signal or background events here we\n",
    "group events with similar values of $y$). The first cut is placed at\n",
    "$x=2.5$ (b) (intuitively, the left points are on average below the right\n",
    "ones). The we continue setting the cuts values at $x=2.0$ (c), $x=3.0$\n",
    "(d), $x=1.0$ (e), $x=4.1$ (f). The procedure stops when we reach a\n",
    "minimum number of points in each of the regions. As in the\n",
    "classification case the minimum number of points gives an handle to\n",
    "limit the over-training. The sequence of cuts is in practice encoded in\n",
    "a binary tree as shown in Fig. [1.8](#fig:BDTregr2){reference-type=\"ref\"\n",
    "reference=\"fig:BDTregr2\"}. Once the regions/leaves are defined we can\n",
    "fit a simple model in each of them. By far the most used fit model is a\n",
    "simple constant ending up as in the classification case with a\n",
    "piece-wise constant function over the variables' space. The values shown\n",
    "in the leaves of Fig. [1.8](#fig:BDTregr2){reference-type=\"ref\"\n",
    "reference=\"fig:BDTregr2\"} represent the results of the fit to a constant\n",
    "shown as red lines in (g). The model used to fit the points in the final\n",
    "leaves can be more involved than a simple constant. If you know that the\n",
    "the points in the leaves will be distributed according to some\n",
    "principle, you can try to fit them with a parametric description of that\n",
    "distribution. The advantage of using a more complex model instead of a\n",
    "simple constant is that you will extract more information from the data\n",
    "(e.g. instead of just getting the mean value in a region you can extract\n",
    "also the width of the distribution or other parameters). To apply the\n",
    "regression tree we just need to take the particular point in the\n",
    "variables' space we want to regress and, as in the case of the\n",
    "classifier, pass it through the decision tree to obtain the regressed\n",
    "value.\n",
    "\n",
    "![Decision tree corresponding to the regression example in\n",
    "Fig. [1.7](#fig:BDTregr1){reference-type=\"ref\"\n",
    "reference=\"fig:BDTregr1\"}.[\\[fig:BDTregr2\\]]{#fig:BDTregr2\n",
    "label=\"fig:BDTregr2\"}](Section10Bilder/BDTregr2.png){#fig:BDTregr2\n",
    "width=\"50%\"}\n",
    "\n",
    "### Stabilizing the decision trees\n",
    "\n",
    "Decision trees as described in the previous sections, cannot be used\n",
    "because are too sensitive to the particular statistical fluctuations of\n",
    "the training sample or in other words they have large variance. The\n",
    "situation changes when we apply aggregation techniques to stabilize the\n",
    "algorithm. The two most commonly used in HEP are \"bagging\" and\n",
    "\"boosting\". In the following we will see how they work when applied to\n",
    "decision trees, but the aggregation methods are completely general and\n",
    "can be applied to any classification/regression MVA.\n",
    "\n",
    "#### Bagging\n",
    "\n",
    "Bagging stands for Bootstrapping AGGregation. The intuition is to use a\n",
    "resampling technique, the bootstrapping, shown in\n",
    "Sec. [\\[sec:resampling\\]](#sec:resampling){reference-type=\"ref\"\n",
    "reference=\"sec:resampling\"} to smooth out the statistical fluctuations\n",
    "of the specific training sample. We will see how it works by using as an\n",
    "example the regression problem in\n",
    "Fig. [1.9](#fig:bagging1){reference-type=\"ref\"\n",
    "reference=\"fig:bagging1\"}.\n",
    "\n",
    "![Decision tree example for bagging.[\\[fig:bagging1\\]]{#fig:bagging1\n",
    "label=\"fig:bagging1\"}](Section10Bilder/bagging1.png){#fig:bagging1\n",
    "width=\"80%\"}\n",
    "\n",
    "Suppose the red curve represents the truth mapping between the\n",
    "variables' space and $y$; $y = f(x)$. The training sample, represented\n",
    "by the grey dots is just one sampling from the truth distribution. The\n",
    "Bagging technique takes N datasets independently drawn with repetition\n",
    "from the initial training sample, for each it builds a tree and it\n",
    "aggregates the resulting trees (computes the average response on the\n",
    "different trees). $$\\begin{aligned}\n",
    "       &D^1 = &\\{ (x_1^{(1)},y_1^{(1)}), (x_2^{(1)},y_2^{(1)}), \\ldots, (x_n^{(1)},y_n^{(1)})  \\}\\\\\n",
    "       &D^2= &\\{ (x_1^{(2)},y_1^{(2)}), (x_2^{(2)},y_2^{(2)}), \\ldots, (x_n^{(2)},y_n^{(2)})  \\}\\\\\n",
    "       &\\ldots& \\\\\n",
    "       &D^N = &\\{ (x_1^{(N)},y_1^{(N)}), (x_2^{(N)},y_2^{(N)}), \\ldots, (x_n^{(N)},y_n^{(N)})  \\}\\end{aligned}$$\n",
    "For each of the resampled datasets $D^i$, given a value of $x$ we will\n",
    "get a regressed value of $y(i)$. Suppose that the estimator of $Y$ is\n",
    "unbiased: $E[Y] = y = f(x)$, the variance (square distance from the true\n",
    "value) is $E[(Y-y)^2] = \\sigma^2(Y)$. If we define\n",
    "$Z =  \\frac{1}{N}\\sum_{i=1}^Ny(i)$, its expectation value is\n",
    "$E[Z] = \\frac{1}{N}\\sum_{i=1}^N y = y$ and its variance is\n",
    "$E[(Z-y)^2] = E[(Z-E[Z])^2] = \\sigma^2(Z) = \\sigma^2(\\frac{1}{N}\\sum_{i=1}^Ny(i)) = \\left( \\frac{1}{N} \\right)^2 \\sigma^2 (\\sum_{i=1}^Ny(i)) = \\frac{1}{N}\\sigma^2(y)$.\n",
    "The larger the number of resampling the smaller the variance.\\\n",
    "The same technique can be applied to classification problems. The only\n",
    "difference is that the classes instead of being defined in a continuous\n",
    "dataset are defined on a discrete one $c_i \\, , \\, i=1, \\ldots, n$ and\n",
    "so the numerical average over the different trees becomes a majority\n",
    "vote.\n",
    "\n",
    "#### Boosting\n",
    "\n",
    "The intuition behind the boosting aggregation technique is to train\n",
    "several trees sequentially, each learning from the errors of the\n",
    "previous ones [@ihler]. Let's take as an example the classification\n",
    "problem in Fig. [1.10](#fig:boost1){reference-type=\"ref\"\n",
    "reference=\"fig:boost1\"}.\n",
    "\n",
    "![Decision tree example for boosting.[\\[fig:boost1\\]]{#fig:boost1\n",
    "label=\"fig:boost1\"}](Section10Bilder/boosting1.png){#fig:boost1\n",
    "width=\"100%\"}\n",
    "\n",
    "We have two variables $(x^1, x^2)$ and two classes: signal in red and\n",
    "background in blue (see Fig. [1.10](#fig:boost1){reference-type=\"ref\"\n",
    "reference=\"fig:boost1\"} (a)). We assign a numerical value to signal and\n",
    "background (e.g. s = +1; b = -1), the reason for this will become clear\n",
    "when we will discuss the details of the boosting algorithm. In (b) we\n",
    "have trained a single level decision tree (which just corresponds to a\n",
    "single cut on the $x^2$ variable). The tree assigns the points above the\n",
    "cut to signal and the ones below to background. It correctly classifies\n",
    "4 signal events while it mis-classifies one background event as signal\n",
    "and it correctly classifies 5 background events but it mis-classifies 3\n",
    "signal events as background. The idea here is to focus on the\n",
    "misclassified events, so we grow a second tree but this time we apply a\n",
    "weight to the events. We increase the weight of the misclassified and\n",
    "reduce the weight to the correctly classified as shown in (c); the\n",
    "larger/smaller events' weights are pictured with larger/smaller fonts.\n",
    "Then we train a new tree like in (d). This time because of the different\n",
    "weights it will obviously differ from the previous one and it will find\n",
    "it to be more discriminating to set a cut on the $x^1$ variable. This\n",
    "tree correctly classifies all the signal events below the cuts,\n",
    "correctly classifies six background events and mis-classifies 3 events\n",
    "above the cut. Repeating the procedure we increase the weights of the\n",
    "misclassified events and decrease the ones of the events we classified\n",
    "correctly as in (e). Then we train a new tree as in (f) and again repeat\n",
    "the procedure.\\\n",
    "As a last step (see Fig. [1.11](#fig:boost2){reference-type=\"ref\"\n",
    "reference=\"fig:boost2\"}) we sum the trees we obtained with some weights\n",
    "which we will derive below. The regions that sum to a positive value\n",
    "will be assigned to signal while the ones obtaining a negative value to\n",
    "background. The results is that the final classifier is more complex and\n",
    "performant than any of the intermediate ones.\n",
    "\n",
    "![Composition of the different trees.[\\[fig:boost2\\]]{#fig:boost2\n",
    "label=\"fig:boost2\"}](Section10Bilder/boosting2.png){#fig:boost2\n",
    "width=\"100%\"}\n",
    "\n",
    "\\\n",
    "Back on how to set the weights. There are several different algorithms\n",
    "to compute the boost weights: AdaBoost, GradientBoost, etc\\... Here we\n",
    "will describe the adaptive boost (AdaBoost). The algorithm is described\n",
    "by the pseudo-code in Fig. [1.12](#fig:boost3){reference-type=\"ref\"\n",
    "reference=\"fig:boost3\"}.\n",
    "\n",
    "![AdaBoost pseudo-code.[\\[fig:boost3\\]]{#fig:boost3\n",
    "label=\"fig:boost3\"}](Section10Bilder/boosting3.png){#fig:boost3\n",
    "width=\"100%\"}\n",
    "\n",
    "The algorithm starts from a uniform set of weights and it evolves them\n",
    "based on the mis-classifications. It starts by training a classifier (in\n",
    "our case a decision tree, but it could be any) based on the training\n",
    "sample (X,Y) and the initial set of weights (w). Then it applies the\n",
    "classifier to the dataset (X) obtaining a prediction $\\hat{Y}$. At this\n",
    "point it checks how many errors it made. It does it by producing a\n",
    "vector of mistakes (Y==$\\hat{Y}$ is a vector with \"1\" where the\n",
    "predicted class is correct and \"0\" when wrong) and taking the scalar\n",
    "product with the vector of weights. The result is a number \"e\" which\n",
    "represents the error rate. With this it computes $\\alpha(i)$ which is a\n",
    "smooth decreasing function of the error rate, for an error rate greater\n",
    "than 50% $\\alpha$ is negative, while it is positive for error rates\n",
    "below 50%). With this coefficient it computes a new vector of weights by\n",
    "multiplying the original one by\n",
    "$w_j \\to w_j \\exp(-\\alpha_i\\cdot(Y_j\\cdot\\hat{Y}_j)$. At the very\n",
    "beginning, we assigned a numerical value of \"+1\" to the signal and \"-1\"\n",
    "to the background. Because of this $Y\\cdot\\hat{Y}$ is a vector of \"+1\"\n",
    "and \"-1\". If the predicted and the true class coincide they will have\n",
    "equal signs and so the product will give a \"+1\", while if the prediction\n",
    "does not match the true value the signs will be opposite and the product\n",
    "give a \"-1\". When multiplied by $\\alpha_i$, the weight of the correct\n",
    "matches will be decreased, the weight of the wrong ones increased.\n",
    "Finally the vector of weights is normalized. This procedure will be\n",
    "repeated Nboost times which is a parameter that can be optimized to get\n",
    "the best performance from the classifier. The final classifier is the\n",
    "weighted sum of all the prediction.\\\n",
    "\n",
    "### Comments on BDT\n",
    "\n",
    "When using any MVA technique is extremely important to check for\n",
    "over/under-traning. Overtraining happens when the complexity of the\n",
    "classifier allows it to learn about the specific statistical fluctuation\n",
    "(noise) of the training dataset. This can be due to a traning sample too\n",
    "small for the number of variables used or a poorly chosen set of\n",
    "parameters in the algorithm training. The result of overtraining is to\n",
    "obtain an artificially good result in the classification/regression\n",
    "because the algorithms learns too many irrelevant detail of the traning\n",
    "sample at hand. Fig. [1.13](#fig:overtraining){reference-type=\"ref\"\n",
    "reference=\"fig:overtraining\"} shows a classifier trained using different\n",
    "parameters on the same dataset.\n",
    "\n",
    "![Examples of a learning algorithm trained in different\n",
    "ways.[\\[fig:overtraining\\]]{#fig:overtraining\n",
    "label=\"fig:overtraining\"} [@CowanMVA]](Section10Bilder/overtraining.pdf){#fig:overtraining\n",
    "width=\"100%\"}\n",
    "\n",
    "The figure on the left shows an example of overtraning. The algorithm is\n",
    "able to get a perfect separation between the two classes on the traning\n",
    "sample but it will have poor performance on a different dataset that\n",
    "will unavoidably have different statistical fluctuations. The middle\n",
    "figure shows another case of poor traning (under-training), where the\n",
    "algorithm complexity was not allowed to grow sufficiently. The non\n",
    "linearities of the training dataset are not model by the algorithm. The\n",
    "right picture shows an algorithm correctly trained.\\\n",
    "A common way to check for over/under training is to split the training\n",
    "sample in two. On the first segment of data the \"training-set\", we train\n",
    "the algorithm, on the second one the \"test-set\" (which is statistically\n",
    "independent from the former) we test its performance. A properly trained\n",
    "algorithm will show small bias and variance on both the training and the\n",
    "test samples.\\\n",
    "An important characteristics of BDTs is that adding correlated variables\n",
    "will have little or no effect on the performance of the algorithm. The\n",
    "reason for this is that the best variable to cut on will be chosen at\n",
    "each step (based for instance on the Gini-index): any useless variable\n",
    "will just never be used.\n",
    "\n",
    "Artificial Neural Networks {#SectionNN}\n",
    "--------------------------\n",
    "\n",
    "Artificial Neural networks (ANN or NN for short) had a large expansions\n",
    "in the '80s '90s. Then the interest in these algorithms diminished, but\n",
    "they are now back as the state of the art technology in machine learning\n",
    "(especially for classification problems) due to the enormous growth in\n",
    "power of modern computers (in particular through the use of GPUs). In\n",
    "HEP, BDT are presently dominating the scene, but it is conceivable,\n",
    "following the machine learning expansion, a return to neural networks in\n",
    "the near future.\\\n",
    "As the name suggests, neural networks were initially inspired by the\n",
    "goal of producing artificial systems to simulate the functioning of a\n",
    "brain. In reality the structure of an artificial neural network is only\n",
    "loosely inspired by nature. The modelling of a basic unit (a neuron) is\n",
    "shown in Fig. [1.14](#fig:neuron){reference-type=\"ref\"\n",
    "reference=\"fig:neuron\"}.\n",
    "\n",
    "![Sketch of a brain cell, the neuron.[\\[fig:neuron\\]]{#fig:neuron\n",
    "label=\"fig:neuron\"} [@neuron]](Section10Bilder/neuron.png){#fig:neuron\n",
    "width=\"50%\"}\n",
    "\n",
    "A neuron receives signals through the dendrites (input wires), process\n",
    "them in the body of the cell (the computational unit) and outputs a\n",
    "signal through the axion. The axion is then connected to one of more\n",
    "other neurons building a network. This basic idea is replicated through\n",
    "software in a neural network [@Ng].\\\n",
    "We will model a neuron as in\n",
    "Fig. [1.15](#fig:logisticUnit){reference-type=\"ref\"\n",
    "reference=\"fig:logisticUnit\"}. A number of inputs are fed to a\n",
    "computational unit which provides an output.\n",
    "\n",
    "![Representation of a computational\n",
    "unit.[\\[fig:logisticUnit\\]]{#fig:logisticUnit\n",
    "label=\"fig:logisticUnit\"}](Section10Bilder/logisticunit.png){#fig:logisticUnit\n",
    "width=\"50%\"}\n",
    "\n",
    "On the left are the inputs (in this case $(x_1, x_2, x_3$)), in the\n",
    "centre the computational unit and on the right the output which is a\n",
    "function of the inputs. Conventionally the first input is $x_1$ is fixed\n",
    "to \"1\" and it is called the bias neuron (this will simplify the\n",
    "vectorization of the method). The vector $\\theta$ contains all the\n",
    "parameters describing the neural network (called \"weights\"). The output\n",
    "function (also called activation function) is typically given by a\n",
    "sigmoid (but any well behaved turn-on function would do):\n",
    "$$g(z) = \\frac{1}{1+e^{-z}}$$ where\n",
    "$z = \\theta_{1,1}~x_1 + \\theta_{1,2}~x_2 + \\theta_{1,3}~x_3$ This simple\n",
    "neural network with just 1 layer is called single layer perceptron.\\\n",
    "More complex networks can be put together as the one in\n",
    "Fig. [1.16](#fig:NN){reference-type=\"ref\" reference=\"fig:NN\"}.\n",
    "\n",
    "![Example of a neural network with one hidden layer.[\\[fig:NN\\]]{#fig:NN\n",
    "label=\"fig:NN\"}](Section10Bilder/NN.png){#fig:NN width=\"80%\"}\n",
    "\n",
    "The first layer is called the input layer (Layer 1), the last layer is\n",
    "the output layer (Layer 3) and the layer in the midlle (which generally\n",
    "might be more than one) is called hidden layer. The $a_i^{(j)}$ is the\n",
    "activation unit $i$ in layer $j$. We can translate this schematics into\n",
    "its corresponding mathematical expression: $$\\begin{aligned}\n",
    "a_1^{(2)} & = & g\\left(\\theta_{1,0}^{(1)}~x_0 + \\theta_{1,1}^{(1)}~x_1 + \\theta_{1,2}^{(1)}~x_2 + \\theta_{1,3}^{(1)}~x_3\\right)\\\\\n",
    "a_2^{(2)} & = & g\\left(\\theta_{2,0}^{(1)}~x_0 + \\theta_{2,1}^{(1)}~x_1 + \\theta_{2,2}^{(1)}~x_2 + \\theta_{2,3}^{(1)}~x_3\\right)\\\\\n",
    "a_3^{(2)} & = & g\\left(\\theta_{3,0}^{(1)}~x_0 + \\theta_{3,1}^{(1)}~x_1 + \\theta_{3,2}^{(1)}~x_2 + \\theta_{3,3}^{(1)}~x_3\\right)\\\\\\end{aligned}$$\n",
    "where $x_0  = 1$ a bias neuron.\\\n",
    "Similarly, the output layer can be written as:\n",
    "$$h_\\theta(x) = a_1^{(3)} =  g\\left(\\theta_{1,0}^{(2)}~a_0^{2} + \\theta_{1,1}^{(2)}~a_1^{2} + \\theta_{1,2}^{(2)}~a_2^{2} + \\theta_{1,3}^{(2)}~a_3^{2}\\right)\\\\$$\n",
    "where $a_0^{2} = 1$ is again a bias neuron. The process of going from\n",
    "the input to the output layer, performing all the computations is called\n",
    "**forward propagation**.\\\n",
    "It is worth noting that each of the nodes of the neural network only\n",
    "performs a *simple action* on the inputs; the complexity of the final\n",
    "result is given by the composition of all the simple actions, i.e. by\n",
    "the network architechture. Another important point is that at each\n",
    "hidden layer, it's the network itself which will decide what the inputs\n",
    "$a_i^{(j)}$ will be. The assignment of the weights is obtained by the\n",
    "traning step.\\\n",
    "The optimization of the weights $\\theta$ of the neural network is done\n",
    "by minimizing a cost function based on the error between the predicted\n",
    "classification and the true one. You can think of it as:\n",
    "$$\\mbox{cost}(i) \\sim (h_\\theta(x^{(i)}) - y^{(i)}) ^2$$ basically how\n",
    "close is the classification $h_\\theta(x^{(i)})$ to the true value\n",
    "$y^{(i)}$.\\\n",
    "The actual definition of the cost function and the algorithm to\n",
    "efficiently solve the minimization problem is beyond the scope of this\n",
    "notes (see [@Ng]). The most used algorithm is called \"back-propagation\"\n",
    "and generally speaking it finds the optimal weights by minimizing the\n",
    "classification error at each layer, starting from the output layer and\n",
    "proceding backwards to the first hidden layer.\n",
    "\n",
    "MVA examples\n",
    "------------\n",
    "\n",
    "**MiniBooNE** [@miniboone] was the first experiment that published an\n",
    "analysis based on a BDT selection. The experiment searches for neutrino\n",
    "oscillation and need to separate events generated by electrons, muons or\n",
    "neutral pions. The identification process is based on Cerenkov radiation\n",
    "in a tank with its inner surface covered with photomultiplier (PMT). The\n",
    "classifier is based on a BDT and the chosen inputs are the number of\n",
    "photomultiplier hits, the energy of the candidate and the radius of the\n",
    "reconstructed Cerenkov rings. Notice that in this case, instead of using\n",
    "as input variables the output of the PMTs, the information is\n",
    "pre-processed and \"high-level\" variables are instead used.\n",
    "\n",
    "![Events from MiniBooNe: (left) an electron candidate, (middle) a muon\n",
    "candidate, (right) a $\\pi^0$\n",
    "candidate.[\\[fig:miniboone\\]]{#fig:miniboone\n",
    "label=\"fig:miniboone\"}](Section10Bilder/miniboone.png){#fig:miniboone\n",
    "width=\"100%\"}\n",
    "\n",
    "\\\n",
    "Another example of a classifier based on a BDT is the **photon\n",
    "identification** in the CMS $H\\to \\gamma \\gamma$ analysis [@hig13001]. A\n",
    "high energy photon is reconstructed with the electromagnetic calorimeter\n",
    "where it develops a narrow electromagnetic shower. A jet faking a photon\n",
    "is typically a jet where a light neutral meson ($\\pi^0$ or $\\eta$) gets\n",
    "most of the momentum of the jet. The neutral meson decays to a photons\n",
    "pair but because of its boost (the analysis searches for events with\n",
    "high energy photons) the two photons are collimated and their showers\n",
    "tend to overlap in the electromagnetic calorimeter. A jet faking a\n",
    "photon will then appear as a shower in the electromagnetic calorimeter\n",
    "with a shower shape which on average will be broader than a true photon.\n",
    "The classifier used in CMS to separate photons from jets faking photons\n",
    "is based on a BDT trained on several input variables which describe the\n",
    "shape of the shower. The algorithm is trained on a simulated (MC) set of\n",
    "photons and jets. The reason why we use simulations instead of collision\n",
    "data is that by simulating ourself the datasets, we know how to assign\n",
    "unabiguously each candidate to its correct category (photon/jet).\\\n",
    "\\\n",
    "Typically the output of a classifier will not be a binary value\n",
    "(\"signal/background\"),but a continuous value on which we will set a cut\n",
    "to separate the two classes. In\n",
    "Fig. [1.19](#fig:classifier){reference-type=\"ref\"\n",
    "reference=\"fig:classifier\"} we see the output of four different\n",
    "classifiers applied to the same toy dataset. An easy way to compare the\n",
    "performance of a classifier is the **\"receiving operator curve\" (ROC)**\n",
    "show in Fig. [1.19](#fig:classifier){reference-type=\"ref\"\n",
    "reference=\"fig:classifier\"} (bottom). This curve shows the background\n",
    "rejection (which is simply $1-\\epsilon_{bkg}$) vs. the signal\n",
    "efficiency. The best classifier will show a curve with a sharp edge in\n",
    "the top right of the plot (both maximal background rejection and maximal\n",
    "signal efficiency). A classifier which assigns randomly events to signal\n",
    "and background will result in a ROC curve which is just the diagonal\n",
    "$1-\\epsilon_{bkg} = 1 - \\epsilon_{sig}$ (i.e.\n",
    "$\\epsilon_{bkg} = \\epsilon_{sig}$). A standard figure of merit for to\n",
    "compare classifiers performance is the area under the ROC curve (AUC).\n",
    "In HEP the performance are often compared by fixing the efficiency at a\n",
    "given value (e.g. 90%) and then comparing the background rejections.\\\n",
    "\n",
    "![Example of a classifier output (left) and the corresponding ROC\n",
    "curves [@TMVA].[\\[fig:classifier\\]]{#fig:classifier\n",
    "label=\"fig:classifier\"}](Section10Bilder/classification1.png \"fig:\"){#fig:classifier\n",
    "width=\"100%\"} ![Example of a classifier output (left) and the\n",
    "corresponding ROC curves [@TMVA].[\\[fig:classifier\\]]{#fig:classifier\n",
    "label=\"fig:classifier\"}](Section10Bilder/classification2.png \"fig:\"){#fig:classifier\n",
    "width=\"60%\"}\n",
    "\n",
    "\\\n",
    "A typical regression problem encountered in HEP comes with **energy\n",
    "corrections**. Let's take again as an example the energy corrections\n",
    "applied on electrons in CMS [@egm11001]. An electron will deposit most\n",
    "of its energy in the electromagnetic calorimeter. However some of it\n",
    "will be emitted by bremstrahlung in the material in front of the\n",
    "calorimeter, some in the non-instrumented gaps of the calorimeter and\n",
    "some might even not be correctly collected by the clustering\n",
    "algorithm[^1]. The regression problem consists in assigning an energy\n",
    "value (which is continuous, that's why is a regression and not a\n",
    "classification) to the electron which is the closest to its generated\n",
    "energy. The algorithm is trained on a simulated (MC) sample of\n",
    "electrons. The reason to use simulation instead of data is to know\n",
    "precisely the true value of energy which is the target of the\n",
    "regression. The input variables chosen for the BDT are: the energy sum\n",
    "obtained by the clustering algorithm, several shower shape variables and\n",
    "the position of the electromagnetic shower in the detector. To show the\n",
    "effect of the regression, CMS used a sample of Z-bosons decaying to\n",
    "$e^+ e^-$. The better the estimation of the energy, the better the\n",
    "energy resolution, the narrower is the invariance mass peak is going to\n",
    "be. In Fig. [1.20](#fig:regression){reference-type=\"ref\"\n",
    "reference=\"fig:regression\"} the blue curve is the invariant mass of the\n",
    "$Z\\to e^+e^-$ using a very simple clustering algorithm where only the\n",
    "energy collected in a 5x5 matrix around the position of the electron is\n",
    "used; in red are the same events, this time using the CMS clustering\n",
    "algorithm (called \"supercluster\" in the legend): by better collecting\n",
    "the energy of the particle, the resolution on the invariant mass\n",
    "improves; the histogram in black shows the final improvement obtained\n",
    "from the energy regression described above on the electrons.\n",
    "\n",
    "![Effect of the energy regression on the $Z\\to e^+ e^-$ invariant\n",
    "mass.[\\[fig:regression\\]]{#fig:regression\n",
    "label=\"fig:regression\"}](Section10Bilder/regression.png){#fig:regression\n",
    "width=\"50%\"}\n",
    "\n",
    "\\\n",
    "\\\n",
    "MVAs are very commonly used also outside HEP. All major companies use\n",
    "them: searching suggestions, translate text, vision, voice recognition,\n",
    "tuning of websites advertisements, look and feel, suggestions for\n",
    "purchasing, etc\\... A simple non HEP example is the face recognition\n",
    "used in digital cameras. Faces are decomposed in basic features: left\n",
    "right symmetry, two darker shapes in the top half of the picture (the\n",
    "eyes), around the middle a vertical darker shape (nose) and a horizontal\n",
    "shape in the bottom half of the picture (the mouth). The algorithm is\n",
    "trained by showing a large sample of signal (faces) and background\n",
    "(non-faces). In Fig. [1.21](#fig:nonHEP){reference-type=\"ref\"\n",
    "reference=\"fig:nonHEP\"} you can see a few examples of the classifier\n",
    "output applied on different pictures.\n",
    "\n",
    "![Example of the output of a face recognition classifier: (top left) a\n",
    "face recognized as such; (top right) a cloud recognized as a face;\n",
    "(bottom left) pine leaves not recognized a faces; (bottom right) a face\n",
    "with a particular make up made to trick the face recognition\n",
    "algorithm.[\\[fig:nonHEP\\]]{#fig:nonHEP\n",
    "label=\"fig:nonHEP\"}](Section10Bilder/nonHEP.png){#fig:nonHEP\n",
    "width=\"70%\"}\n",
    "\n",
    "References\n",
    "----------\n",
    "\n",
    "-   Cowan academic training lectures on MVA 2008\n",
    "\n",
    "-   Machine learning classes - youtube\n",
    "\n",
    "-   Alexander Ihler, University of California Irvine on \"Boosting\" -\n",
    "    youtube\n",
    "\n",
    "-   TMVA [@TMVA]: tmva.sourceforge.net\n",
    "\n",
    "-   Andrew Ng, Stanford: \"Machine learning\" at Coursera\n",
    "\n",
    "[^1]: The CMS calorimeter is build from thousands of $PbWO_4$ crystals\n",
    "    each of which is readout by photo-sensitive detector. The light\n",
    "    detected by these sensors is a function of the energy of the\n",
    "    inpinging particles. A clustering algorithm groups crystals with\n",
    "    some significant energy deposit and creates a particle candidate.\n"
   ]
  }
 ],
 "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
}