leastSquares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter Estimation - Least Squares\n",
    "====================================\n",
    "\n",
    "Together with the maximum likelihood method, the method of least squares\n",
    "(LS or LSQ) is very often used in parameters estimation. In this chapter\n",
    "we will give a description of the method together with a few examples\n",
    "and show its relation with the ML method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Least Squares Method\n",
    "------------------------\n",
    "\n",
    "Let's take a set of independent Gaussian random variables $y_i$, with\n",
    "$i=1,\\ldots,N$ and let's assume that each $y_i$ is distributed around an\n",
    "unknown mean $\\mu_i$ with variance $\\sigma^2_i$, where the mean is\n",
    "predicted by a function $f(x_i;a)$ {numref}`fig-chi2Sketch`.\n",
    "In the typical application the $y_i$ are\n",
    "the (independent) measurements with uncertainty $\\sigma_i$ and\n",
    "$f(x_i;a)$ is the functional form of the \"model\" for which you are\n",
    "interested in estimating the value of some parameters (in this case\n",
    "$a$).\n",
    "\n",
    "```{figure} ./images/ch7/chi2Sketch.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-chi2Sketch\n",
    "---\n",
    "FIXME: Sketch to illustrate the notation.\n",
    "```\n",
    "\n",
    "If the data $\\{y_i\\}$ are Gaussian distributed around the mean\n",
    "$f(x_i;a)$ then the sum: \n",
    "\n",
    "```{math}\n",
    ":label: eq-ChisquareEquation\n",
    "\\chi^{2} = \\sum_{i=1}^{N}\\left(\\frac{y_{i} -f(x_{i};a)}{\\sigma_{i}} \\right)^{2}\n",
    "```\n",
    "\n",
    "obeys a $\\chi^{2}$-distribution with $(N-1)$ degrees of freedom (the\n",
    "number of measurements, minus the number of fitted parameters). In the\n",
    "general case of $p$-parameters to be fitted\n",
    "($f(x_i;a)\\to f(x_i;\\vec{a})$ ) the number of degrees of freedom will be\n",
    "$N-p$.\n",
    "\n",
    "To find the best estimate for the parameter $a$ we proceed in a way\n",
    "similar to what we have done in the ML-method: we look for the value of\n",
    "the parameter $a$ which minimize the $\\chi^2$. If you interpret the\n",
    "$\\chi^2$ as a distance, the LS method corresponds to minimize the\n",
    "distance between the measured data and the considered model. This boils\n",
    "down to calculating the minimum of the $\\chi^2$ as:\n",
    "\n",
    "$$\n",
    "\\frac{d\\chi^2}{da} = 2 \\sum_{i=1}^N \\frac{d}{da}f(x_i;a) \\cdot \\frac{y_i - f(x_i ; a)}{\\sigma_i^2} = 0.\n",
    "$$\n",
    "\n",
    "In case of $p$ parameters $a_1,\\ldots,a_p$ and $f(x;\\vec{a})$, the idea\n",
    "is the same, just the minimization will have to be performed\n",
    "simultaneously in $p$ dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connection to the Likelihood Function\n",
    "\n",
    "The simplest way to see the relation between the LS and the ML methods\n",
    "is to take a set of data $(x_{i},y_{i})$ for which the $x_{i}$ are known\n",
    "precisely and the $y_{i}$ are known with uncertainties $\\sigma_{i}$.\n",
    "Under the assumption that the $y_{i}$ are Gaussian distributed (for\n",
    "instance when coming from the CLT), the probability to observe $y_i$\n",
    "given the prediction $f(x_i;a)$ is:\n",
    "\n",
    "$$\n",
    "p(y_i|a)=\\frac{1}{\\sigma_i\\sqrt{2\\pi}}e^{-(y_i-f(x_i|a))^2/2\\sigma_i^2}\n",
    "$$\n",
    "\n",
    "From this we can build the likelihood function for the complete data set\n",
    "as: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L(a,\\vec{y})&=&\\prod_ip(y_i|a) \\\\\n",
    "\\ln L(a,\\vec{y})&=&-\\frac{1}{2}\\sum_{i} \\left(\\frac{y_i-f(x_i|a)}{\\sigma_i}\\right)^2\n",
    "-\\sum_{i} \\ln\\sigma_i\\sqrt{2\\pi},\\end{aligned}\n",
    "$$ \n",
    "\n",
    "where only the first\n",
    "term depends on $a$. To minimize the negative log-likelihood as a\n",
    "function of the parameter $a$ we will have to minimize:\n",
    "\n",
    "$$\n",
    "\\chi^2=\\sum_{i} \\left(\\frac{y_i-f(x_i|a)}{\\sigma_i}\\right)^2\n",
    "$$ \n",
    "\n",
    "which corresponds to the $\\chi^2$ procedure shown in the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitting a Straight Line\n",
    "-----------------------\n",
    "\n",
    "As a first example of the application of the LS method, we take a set of\n",
    "$N$ independent measurements $(x_i,y_i)$ where we assume that the model\n",
    "is linear, and in particular that $f(x) = mx$ (i.e. a straight line\n",
    "passing through the origin). The quantity which has to be minimized is\n",
    "then: \n",
    "\n",
    "$$\n",
    "\\label{}\n",
    "\\chi^2=\\sum_i\\left(\\frac{y_i-mx_i}{\\sigma_i}\\right)^2\n",
    "$$ \n",
    "\n",
    "We furthermore\n",
    "assume that all the uncertainties are the same:\n",
    "$\\sigma_{i} = \\sigma \\, \\, \\forall i$. Differentiating with respect to\n",
    "$m$ and equating to zero to obtain the best estimator $\\hat{m}$, we\n",
    "have: \n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\frac{\\partial \\chi^2}{\\partial m}&=&-\\frac{2}{\\sigma^2}\\sum_i(x_iy_i-mx_i^2)\\\\\n",
    "\\sum_i(x_iy_i-mx_i^2) &=&0 \\\\\n",
    "\\sum_i x_iy_i&=&m\\sum_i x_i^2 \\\\\n",
    "\\hat{m}&=&\\sum_i\\frac{x_i y_i}{N\\overline{x^2}}=\\frac{\\overline{xy}}{\\overline{x^2}} \\qquad \\left(\\sum x_i^2 = N \\overline{x^2}\\right)\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "The variance of $\\hat{m}$ can be determined by error propagation to be:\n",
    "\n",
    "$$\n",
    "V(\\hat{m})=\\sum_i\\left(\\frac{x_i}{N\\overline{x^2}}\\right)^2\\sigma^2=\\frac{\\sigma^2}{N\\overline{x^2}}.\n",
    "$$\n",
    "\n",
    "When the straight line does not pass through the origin,\n",
    "$f(x_{i};m,b) = mx_{i} + b$ the solution of the LS method is:\n",
    "\n",
    "```{math}\n",
    ":label: eq-ls_slope-intercept\n",
    "\\begin{aligned}\n",
    "\\hat{m} &=& \\frac{\\overline{xy}-\\bar{x}\\bar{y}}{\\overline{x^2}-\\bar{x}^2} \\\\\n",
    "\\hat{b} &=& \\bar{y}-\\hat{m}\\bar{x} \n",
    "\\end{aligned}\n",
    "```\n",
    "\n",
    "The variances are given by: \n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\sigma_m^2 &=& V(\\hat{m})=\\frac{\\sigma^2}{N(\\overline{x^2}-\\bar{x}^2)} \\\\\n",
    "\\sigma_b^2 &=& V(\\hat{b})=\\frac{\\sigma^2\\overline{x^2}}{N(\\overline{x^2}-\\bar{x}^2)} \\\\\n",
    "cov(\\hat{m},\\hat{b}) &=& -\\frac{\\sigma^2\\bar{x}}{N(\\overline{x^2}-\\bar{x}^2)} \\\\\\end{aligned}\n",
    "$$\n",
    "\n",
    "The $\\chi^{2}$ for the best fit is:\n",
    "\n",
    "$$\n",
    "\\chi^2 =\\frac{V(y)}{\\sigma^2}(1-\\rho^2(x,y))\n",
    "$$\n",
    "\n",
    "Note that $V(y)$ is not the same as $\\sigma^{2}$! \n",
    "\n",
    "$V(y) = \\bar{y^{2}} - \\bar{y}^{2}$ is the\n",
    "variance of the whole data sample, whereas $\\sigma$ describes the\n",
    "standard deviation of a single measurement around its true value (here\n",
    "we assumed $\\sigma_i = \\sigma \\; \\forall i$).\n",
    "\n",
    "If the errors $\\sigma_{i}$ are not all the same, we have to minimize the\n",
    "following expression: \n",
    "\n",
    "$$\n",
    "\\sum_i\\frac{(y_i-mx_i-b)^2}{\\sigma_i^2}.\n",
    "$$ \n",
    "\n",
    "The solution of this minimization can again be obtained by the equations\n",
    "given above, if the two means $\\bar{x}$ and $\\bar{y}$ are replaced by\n",
    "the weighted means. Furthermore the normalization is no longer given\n",
    "simply given by $N$ but it is given by $\\sum_{i} 1/\\sigma_{i}^{2}$:\n",
    "\n",
    "$$\n",
    "\\frac{\\sum_i y_i}{N}\\to \\frac{\\sum_iy_i/\\sigma_i^2}{\\sum_i 1/\\sigma_i^2}\n",
    "$$\n",
    "\n",
    "And the quantity $\\sigma^{2}$ has to be replaced in the expressions for\n",
    "the variance by \n",
    "\n",
    "$$\n",
    "\\sigma^2\\to \\frac{N}{\\sum_i 1/\\sigma_i^2}.\n",
    "$$\n",
    "\n",
    "For the sake of completeness we give here the full solution in a form\n",
    "which is particularly useful when plugged into a program. We define:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "S_1&=&\\sum_i w_i=\\sum_i 1/\\sigma_i^2 \\\\\n",
    "S_x &=&\\sum w_i x_i \\quad S_y=\\sum_i w_iy_i \\\\\n",
    "S_{xx}&=&\\sum_iw_ix_i^2 \\quad S_{xy}=\\sum_iw_ix_iy_i \\quad S_{yy}=\\sum_iw_iy_i^2\\end{aligned}\n",
    "$$\n",
    "\n",
    "With this notation we can rewrite the linear system of equation as\n",
    "\n",
    "$$\\left(\\begin{array}{cc}\n",
    " S_1  & S_x \\\\\n",
    "S_x & S_{xx} \\\\\n",
    "\\end{array}\\right)\\cdot {b \\choose m} ={S_y \\choose S_{xy}}.\n",
    "$$\n",
    "\n",
    "With the expression of the determinant $D$: \n",
    "$$\n",
    "\\begin{aligned}\n",
    "D&=&S_1S_{xx}-S_x^2 \\\\\\end{aligned}\n",
    "$$ \n",
    "\n",
    "we can compute: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{b}&=&(S_{xx}S_y-S_xS_{xy})/D \\\\\n",
    "\\hat{m}&=&(-S_{x}S_y+S_1S_{xy})/D\\end{aligned}\n",
    "$$\n",
    "\n",
    "Once the slope $\\hat{m}$ and the intercept $\\hat{b}$ are obtained, we\n",
    "can compute the uncertainty for an arbitrary interpolated (or\n",
    "extrapolated) point $y$ for a given $x$. For such a given $x$ the\n",
    "predicted value for $y$ is just $y = \\hat{m} x + \\hat{b}$ and the error\n",
    "for the interpolated value $y$ is then: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "V(y)&=&V(\\hat{m}x+\\hat{b})=V(\\hat{m}x)+V(\\hat{b})+2\\cdot cov(\\hat{m}x,\\hat{b})\\\\\n",
    "V(y)&=&x^2V(\\hat{m})+V(\\hat{b})+2x\\cdot cov(\\hat{m},\\hat{b}) \\\\\n",
    "V(y)&=&\\frac{\\sigma^2(x-\\bar{x})^2}{N(\\overline{x^2}-\\bar{x}^2)}+\\sigma^2/N\\end{aligned}\n",
    "$$\n",
    "\n",
    "We will see in Ch. [MVA](mva.html#definitions) how this will be used as a simple learning\n",
    "algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Considering the systematic Errors\n",
    "\n",
    "As an example we consider a straight line fit where all the measurements\n",
    "$y_{i}$ have a common statistic error $\\sigma$ and a common systematic\n",
    "error $S$. We know from the discussion about systematic errors in\n",
    "Ch. [Uncertainties](errors.html#statistical-vs-systematic-uncertainties)\n",
    "that the covariance matrix $cov(y_{i},y_{j})$ can\n",
    "be written as $cov(y_{i},y_{j}) = \\delta_{ij}\\sigma^{2} +S^{2}$. The\n",
    "estimators for the slope and the intercept, $\\hat{m}$ and $\\hat{b}$,\n",
    "respectively, are again given by {eq}`eq-ls_slope-intercept`. \n",
    "\n",
    "The complete formula for the variances reads\n",
    "therefore as follows: \n",
    "\n",
    "```{math}\n",
    ":label: eq-summand\n",
    "\\begin{aligned}\n",
    "V(\\hat{m})&=&\\frac{1}{N^2(\\overline{x^2}-\\bar{x}^2)^2}\\sum_{i,j}(x_i-\\bar{x})(x_j-\\bar{x})\\cdot cov(y_i,y_j) \\\\\n",
    "&=&\\frac{1}{N^2(\\overline{x^2}-\\bar{x}^2)^2}\\left(\\sum_i(x_i-\\bar{x})^2\\sigma^2+\n",
    "\\sum_{i,j}(x_i-\\bar{x})(x_j-\\bar{x})S^2\\right)\\\\\n",
    "&=&\\frac{1}{N^2(\\overline{x^2}-\\bar{x}^2)^2}\\left(\\sum_i(x_i-\\bar{x})^2\\sigma^2\\right)\n",
    "\\end{aligned}\n",
    "```\n",
    "\n",
    "The second summand in Eq. {eq}`eq-summand` vanishes, because $1/N \\sum x_{i} = \\bar{x}$.\n",
    "\n",
    "The variance for $\\hat{b}$ is:\n",
    "\n",
    "$$\n",
    "V(\\hat{b})=\\frac{1}{N^2(\\overline{x^2}-\\bar{x}^2)^2}\\sum_{i,j}(\\overline{x^2}-\\bar{x}x_i)\n",
    "(\\overline{x^2}-\\bar{x}x_j)\\cdot cov(y_i,y_j)\n",
    "$$ \n",
    "\n",
    "The sum $\\sum_i (\\overline{x^2}-\\bar{x}x_i)=N(\\overline{x^2}-\\bar{x}^2)$ does\n",
    "not vanish in this equation, thus an additional term appears which is\n",
    "just $S^{2}$; a common systematic error only influences the variance of\n",
    "the intercept, but it does not change the variance of the slope, as we\n",
    "would have naively expected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Straight Line Fit with Errors on both Variables\n",
    "\n",
    "Now we allow for both variables $x_{i}$ and $y_{i}$ to have errors\n",
    "$\\sigma_{x_{i}}$ and $\\sigma_{y_{i}}$ (see {numref}`fig-chi2UncertaintyY` and {numref}`fig-chi2UncertaintyXY`). \n",
    "This means that the sum of the squared residuals of the error ellipsis of the\n",
    "straight line has to be minimized, i.e.:\n",
    "\n",
    "$$\n",
    "S(m,b)=\\sum_i\\frac{(y_i-mx_i-b)^2}{\\sigma^2_{y_i}+m^2\\sigma^2_{x_i}}\n",
    "$$\n",
    "\n",
    "```{figure} ./images/ch7/chi2UncertaintyY.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-chi2UncertaintyY\n",
    "---\n",
    "FIXME: LS fit with uncertainty only on \"y\" or both on \"x\" and \"y\"\n",
    "```\n",
    "\n",
    "```{figure} ./images/ch7/chi2UncertaintyXY.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-chi2UncertaintyXY\n",
    "---\n",
    "FIXME: LS fit with uncertainty only on \"y\" (left) or both on \"x\"\n",
    "and \"y\".\n",
    "```\n",
    "\n",
    "As usual, the two equations $\\partial S / \\partial m = 0$ and\n",
    "$\\partial S / \\partial b = 0$ have to be solved.\\\n",
    "The condition $\\partial S / \\partial b = 0$ leads to\n",
    "\n",
    "$$\n",
    "\\hat{b}=\\frac{\\sum y_i/\\kappa_i-\\hat{m}\\sum x_i/\\kappa_i}{\\sum 1/\\kappa_i}\n",
    "$$\n",
    "\n",
    "where $\\kappa_i=\\sigma^2_{y_i}+m^2\\sigma^2_{x_i}$. For $\\hat{m}$, if the\n",
    "errors are all the same on $x$ (i.e. $\\sigma_{x_{i}} = \\sigma_{x}$) and\n",
    "on $y$ (i.e. $\\sigma_{y_{i}} = \\sigma_{y}$), then the solution for the\n",
    "straight line fit is given by: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{m}&=&\\frac{\\sigma_x}{\\sigma_y}(A\\pm\\sqrt{A^2+1}) \\\\\n",
    "A&=&\\frac{\\sigma_x^2V(y)-\\sigma_y^2V(x)}{2\\sigma_x\\sigma_y\\cdot cov(x,y)} \\\\\n",
    "\\bar{y}&=&\\hat{m}\\bar{x}+\\hat{b}\n",
    "\\end{aligned}\n",
    "$$ \n",
    "\n",
    "Other particular cases can be found making assumption on the uncertainties, but to solve the\n",
    "generic problem numerical techniques are typically used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix Notation and the uncertainty on the fitted parameters\n",
    "------------------------------------------------------------\n",
    "\n",
    "Let ${\\bf a}$ be a vector of $n$ parameters $\\{a_{i}\\}, i=1,\\ldots,n$.\n",
    "The $N$ measured data points can be represented by a vector\n",
    "${\\bf y} = \\{y_{i}\\},~ i=1,\\ldots,N$ and the function to be fitted as\n",
    "$f(x_{i};{\\bf a})$. The $\\chi^2$ expression in matrix form becomes:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\chi^2&=&\\sum_i\\sum_j[y_i-f(x_i;{\\bf a})]V_{ij}^{-1}[y_j-f(x_j;{\\bf a})]\\\\\n",
    "      &=&({\\bf y}-{\\bf f})^T{\\bf V}^{-1}({\\bf y}-{\\bf f})={\\bf r^T}{\\bf V}^{-1}{\\bf r}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Here, ${\\bf r} = {\\bf y} - {\\bf f}$ is the vector of residuals and\n",
    "$\\bf{V}$ is the covariance matrix. By differentiating $\\chi^{2}$ w.r.t\n",
    "each $a_{i}$ and equating each of them to zero yields $n$ equations,\n",
    "which have to be solved simultaneously in order to get the estimator\n",
    "${\\bf \\hat{a}}$. Many mathematical packages (e.g. MatLab, Octave,\n",
    "etc\\...) offer a way to solve matrix problems; in slang moving from the\n",
    "single components equations to the matrix form is called\n",
    "\"vectorization\". Because they use highly optimized algorithms, it is\n",
    "always preferable to work with the vectorized version of the problem.\n",
    "\n",
    "Let $f(x;{\\bf a})$ be a linear function of the parameters ${\\bf a}$:\n",
    "$f(x;{\\bf a})=\\sum_r c_r(x) a_r$. Where the linearity is on the\n",
    "${\\bf a}$, the $c_r(x)$ can be any function of $x$. Written in matrix\n",
    "form, it reads:\n",
    "\n",
    "$$\n",
    "\\chi^2=({\\bf y}-{\\bf Ca})^T{\\bf V}^{-1}({\\bf y}-{\\bf Ca}).\n",
    "$$\n",
    "\n",
    "If we have $N$ data points and $n$ coefficients $(n \\le N)$, then ${\\bf y}$\n",
    "and ${\\bf a}$ are column vectors with dimension $N$ and $n$,\n",
    "respectively. The covariance matrix ${\\bf V}$ has dimension $N \\times N$\n",
    "and the matrix ${\\bf C}$ has dimensions $N \\times n$.\n",
    "\n",
    "Minimizing the $\\chi^2$, the equations $\\partial \\chi^{2} / \\partial {\\bf a} = 0$ are:\n",
    "\n",
    "$$\n",
    "\\partial \\chi^{2} / \\partial {\\bf a} = -2 ( \\bf{C}^T \\bf{V}^{-1} {\\bf \\bf{y}} - \\bf{C}^T \\bf{V}^{-1}\\bf{C} {\\bf a}) = 0 \n",
    "$$\n",
    "\n",
    "The solution for the estimator $\\hat{ {\\bf a}}$ is then:\n",
    "\n",
    "$$\n",
    "{\\bf \\hat{a}}=({\\bf C^T}{\\bf V}^{-1}{\\bf C})^{-1}{\\bf C^T}{\\bf V}^{-1}{\\bf y} := \\bf{B}{\\bf y}\n",
    "$$\n",
    "\n",
    "which means that the solution ${\\bf a}$ are linear functions of the\n",
    "measurements ${\\bf y}$.\n",
    "\n",
    "Using error propagation we can find the covariance matrix for the\n",
    "${\\bf \\hat{a}}$:\n",
    "\n",
    "$$\n",
    "\\bf{U} = \\bf{B} \\bf{V} \\bf{B}^T = [{\\bf C^TV}^{-1}{\\bf C}]^{-1}\n",
    "$$\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Let's go back to the fit of a straight line of the form\n",
    "$f(x) = mx + b$ to $N$ data points, which have independent and common\n",
    "errors, such that ${\\bf V} = \\sigma^{2} I$, i.e.\n",
    "$V_{ij} = \\sigma^{2} \\delta_{ij}$. In matrix notation we then get:\n",
    "\n",
    "$$\n",
    "\\begin{array}{cc}\n",
    "f_1 & = b+ mx_1 \\\\\n",
    "f_2 & = b+ mx_2 \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "f_N & = b+ mx_N  \\\\\n",
    "\\end{array} \n",
    "\\qquad\n",
    "{\\bf C} = \\left( \\begin{array}{cc}\n",
    "1 & x_1 \\\\\n",
    "1 & x_2 \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "1 & x_N \\\\\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\bf \\hat{a}} = \\sigma^2({\\bf C^TC})^{-1}\\frac{1}{\\sigma^2}{\\bf C^Ty}\n",
    "$$\n",
    "\n",
    "which can be explicitly written as: \n",
    "\n",
    "$$\n",
    "\\label{matrix_ls}\n",
    "{\\bf \\hat{a}}=\\left(\\begin{array}{c}\n",
    "\\hat{b} \\\\\n",
    "\\hat{m}\n",
    "\\end{array} \\right)\n",
    "=\n",
    "\\left(\\begin{array}{cc}\n",
    "\\sum_i 1 & \\sum_i x_i \\\\\n",
    "\\sum_i x_i & \\sum_i x_i^2 \\\\\n",
    "\\end{array} \\right)^{-1}\n",
    "\\left(\\begin{array}{c}\n",
    "\\sum_i y_i \\\\\n",
    "\\sum_i x_iy_i\n",
    "\\end{array} \\right)\n",
    "$$ \n",
    "\n",
    "The inversion of the $2 \\times 2$-matrix gives:\n",
    "\n",
    "$$\n",
    "\\frac{1}{N(\\overline{x^2}-\\bar{x}^2)}\n",
    "\\left(\\begin{array}{cc}\n",
    "\\overline{x^2}& -\\bar{x} \\\\\n",
    "-\\bar{x} & 1 \\\\\n",
    "\\end{array} \\right)\n",
    "$$ \n",
    "\n",
    "which finally leads to:\n",
    "\n",
    "$$\n",
    "{\\bf \\hat{a}}=\\left(\\begin{array}{c}\n",
    "\\hat{b} \\\\\n",
    "\\hat{m}\n",
    "\\end{array} \\right)\n",
    "=\n",
    "\\frac{1}{N(\\overline{x^2}-\\bar{x}^2)}\n",
    "\\left(\\begin{array}{cc}\n",
    "\\overline{x^2}& -\\bar{x} \\\\\n",
    "-\\bar{x} & 1 \\\\\n",
    "\\end{array} \\right)\n",
    "\\left(\\begin{array}{c}\n",
    "\\sum_i y_i \\\\\n",
    "\\sum_i x_iy_i\n",
    "\\end{array} \\right).\n",
    "$$ \n",
    "\n",
    "Which corresponds to the expressions for\n",
    "$\\hat{m}$ and $\\hat{b}$ which we extracted\n",
    "in {eq}`eq-ls_slope-intercept`.\n",
    "\n",
    "The variance for the estimator\n",
    "$\\hat{ {\\bf a}}$ is: \n",
    "\n",
    "$$\n",
    "{\\bf V}({\\bf\\hat{a}})=\\left(\\begin{array}{cc}\n",
    "V(b)& cov(b,m) \\\\\n",
    "cov(b,m) & V(m) \\\\\n",
    "\\end{array} \\right)\n",
    "=\n",
    "\\frac{\\sigma^2}{N(\\overline{x^2}-\\bar{x}^2)}\\left(\\begin{array}{cc}\n",
    "\\overline{x^2}& -\\bar{x} \\\\\n",
    "-\\bar{x} & 1 \\\\\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "In a second example we fit the parabola\n",
    "$f(x) = a_{0} + a_{1} x + a_{2} x^{2}$ to $N$ data points. Again we\n",
    "assume the errors being independent and equal for all data points. The\n",
    "matrix ${\\bf C}$ is now given by: \n",
    "\n",
    "$$\n",
    "{\\bf C}=\\left( \\begin{array}{ccc}\n",
    "1 & x_1 & x_1^2\\\\\n",
    "1 & x_2 & x_2^2\\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "1 & x_N & x_N^2 \\\\\n",
    "\\end{array} \\right)\n",
    "$$ \n",
    "\n",
    "Going through the same steps as for the linear\n",
    "case we obtain: \n",
    "\n",
    "$$\n",
    "{\\bf \\hat{a}}=\n",
    "\\left( \\begin{array}{c}\n",
    "\\hat{a_0} \\\\\n",
    "\\hat{a_1} \\\\\n",
    "\\hat{a_2}\\\\\n",
    "\\end{array} \\right)\n",
    "=\\left( \\begin{array}{ccc}\n",
    "\\sum_i 1 & \\sum_i x_i & \\sum_i x_i^2 \\\\\n",
    "\\sum_i x_i & \\sum x_i^2 & \\sum_i x_i^3 \\\\\n",
    "\\sum_i x_i^2 & \\sum x_i^3 & \\sum_i x_i^4 \\\\\n",
    "\\end{array} \\right)^{-1}\n",
    "\\left( \\begin{array}{c}\n",
    "\\sum_i y_i  \\\\\n",
    "\\sum_i x_iy_i \\\\\n",
    "\\sum x_i^2y_i\\\\\n",
    "\\end{array} \n",
    "\\right)\n",
    "$$ \n",
    "\n",
    "The generalization of this method to higher-order\n",
    "polynomials follows the same pattern.\n",
    "```\n",
    "\n",
    "Another way to write the (inverse) of the covariance matrix is:\n",
    "\n",
    "$$\n",
    "(U^{-1})_{ij} = \\frac{1}{2}\\left[ \\frac{\\partial^2 \\chi^2}{\\partial a_i \\partial a_j }  \\right]_{\\bf a = \\hat{a}}\n",
    "$$\n",
    "\n",
    "which is the expression of the RCF bound in the case that the\n",
    "measurements ${\\bf y}$ are Gaussian distributed, and where we use the\n",
    "relation with the log likelihood $- \\ln L = \\chi^2/2$ (a part from an\n",
    "overall constant).\n",
    "\n",
    "Again in case of a function $f$ linear in ${\\bf a}$, the $\\chi^{2}$\n",
    "becomes quadratic in ${\\bf a}$:\n",
    "\n",
    "$$\n",
    "\\chi^2({\\bf a})=\\chi^2({\\bf \\hat{a}})+\\frac{1}{2}\\sum_{i,j=1}\\left[\\frac{\\partial^2\\chi^2}{\\partial a_i\\partial a_j}\\right]_{a=\\hat{a}}(a_i-\\hat{a}_i)(a_j-\\hat{a}_j)\n",
    "$$\n",
    "\n",
    "Using the expression of the variance we just found, the equation above\n",
    "corresponds to contours in the parameter space defined by\n",
    "$\\hat{a}_i\\pm \\hat{\\sigma_i}$ and therefore giving the $\\pm 1 \\sigma$\n",
    "deviations from the estimators:\n",
    "\n",
    "$$\n",
    "\\chi^2({\\bf a}\\pm \\hat{\\sigma})=\\chi^2({\\bf \\hat{a}})+1=\\chi^2_{min}+1.\n",
    "$$\n",
    "\n",
    "Hence the $\\chi^{2}$-function is a parabola with a minimum at $\\hat{a}$\n",
    "and the errors $\\sigma$ for the estimators are determined by\n",
    "$\\chi^{2}_{min}+1$. In general, if the function $f$ is not linear in the\n",
    "parameters, the contour will not be an ellipse, but it will still define\n",
    "a confidence region reflecting the statistical uncertainty on the fitted\n",
    "parameters. The precise construction of the confidence region will be\n",
    "developed in Sec. [Confidence Limits](confidenceIntervals.html#confidence-belt-neyman-frequentist-construction-sec-neyman). It is important to notice that the\n",
    "confidence level of the region defined by the contour, depends on the\n",
    "number of parameters fitted: 6.83% for one parameter, 39.4% for two,\n",
    "19.9% for three, etc\\... (see again Sec. [Confidence Limits](confidenceIntervals.html#confidence-belt-neyman-frequentist-construction-sec-neyman))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binned $\\chi^{2}$ fit\n",
    "---------------------\n",
    "\n",
    "In this section we will apply the LS method to binned data. In the case\n",
    "of $f$ being a probability distribution function (a p.d.f. instead of\n",
    "any generic function), we can interpret the value of $f$ integrated over\n",
    "a given range (\"bin\"), as the expected number of events in that bin\n",
    "$f_i = E[y_i]$: \n",
    "$$\n",
    "\\label{eq:LS}\n",
    "f_i({\\bf a}) = n\\int_{x_i^{min}}^{x_i^{max}} f(x;{\\bf a}) dx = n p_i({\\bf a})\n",
    "$$\n",
    "\n",
    "where $x_i^{min}$ and $x_i^{max}$ are the bin boundaries, the\n",
    "$p_i({\\bf a})$ is the probability to have an event in the bin $i$, given\n",
    "the parameters ${\\bf a}$ and n is the overall normalization.\\\n",
    "The fit proceeds as before minimizing the $\\chi^2$ w.r.t. the parameters\n",
    "${\\bf a}$:\n",
    "\n",
    "$$\n",
    "\\chi^2(y_i|{\\bf a})=\\sum_i \\frac{(y_i - f_i({\\bf a}))^2}{\\sigma_i^2}\n",
    "$$\n",
    "\n",
    "where here $\\sigma_i$ is the variance of the *expected* number of\n",
    "entries in bin $i$. If the number of entries in bin $i$ is small\n",
    "compared to the total number of entries in the histogram then we can\n",
    "assume that they are Poisson distributed and the variance is equal to\n",
    "the mean $\\sigma_i^2 = f_i({\\bf a}) = n p_i({\\bf a})$.\\\n",
    "Often, instead of using the variance of the expected number of entries\n",
    "in bin $i$, the variance of the *observed* number of entries in bin $i$\n",
    "is used, leading to:\n",
    "\n",
    "$$\n",
    "\\chi^2(y_i|{\\bf a})=\\sum_i \\frac{(y_i - f_i({\\bf a}))^2}{y_i}.\n",
    "$$\n",
    "\n",
    "This new expression is called the **modified LS** method. Even if\n",
    "computationally easier to implement (contraty to $\\sigma_i^2$ the\n",
    "observed $y_i$, obviously, does not depend on $\\bf{a}$), it brings in\n",
    "the problem of the variance estimation for bins which are poorly\n",
    "populated or have no entries at all. As a rule of thumb, try to have at\n",
    "least 5 entries per bin. Two situations are rather typical: small\n",
    "statistics in the tails of a distribution, or the whole histogram is\n",
    "sparsely populated. In the first case, try to rebin, for the latter just\n",
    "move to an unbinned ML fit to use the full information you have in your\n",
    "data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use of the $\\chi^2$ to test the goodness of fit \n",
    "-----------------------------------------------\n",
    "\n",
    "If the $\\chi^{2}$ is large after minimization, then the function is\n",
    "probably badly chosen (i.e. not correctly describing the data);\n",
    "intuitively, the $\\chi^{2}$ should be small if the function describes\n",
    "the data. On the other hand, the $\\chi^{2}$ can be small if you have too\n",
    "many degrees of freedom (fit n-points with a polynomial of degree n) or,\n",
    "even with a bad model, when the uncertainties are overestimated\n",
    "($\\sigma_i^2$ sits at the denominator). You can always get very small\n",
    "$\\chi^{2}$ if you assume large enough uncertainties! If the errors are\n",
    "too large the whole method loses its predictive power.\n",
    "\n",
    "We have already encountered in [Chi2 distribution](probabilityDistributions.html#chi-2-distribution) \n",
    "the $\\chi^{2}$-distribution:\n",
    "\n",
    "$$\n",
    "f(\\chi^2;n)=\\frac{2^{-n/2}}{\\Gamma(n/2)}\\chi^{n-2}e^{-\\chi^2/2}.\n",
    "$$ \n",
    "\n",
    "The distribution depends on $n$, the number of degrees of freedom, which is\n",
    "the number of data points minus the number of parameters of the model.\n",
    "Because the $\\chi^{2}$ distribution has expectation value $n$ and\n",
    "variance $2n$, we expect the $\\chi^{2}$ divided by the number of degrees\n",
    "of freedom to be approximately one: $\\chi^{2} / n \\approx 1$. If the\n",
    "$\\chi^{2} / n$ is much larger than one, the data are not properly\n",
    "described by the model. Let's introduce the definition of $p-$value to\n",
    "quantify the level of agreement of the model with data:\n",
    "\n",
    "$$\n",
    "p=\\int_{\\chi^2}^\\infty f(x';n) dx' \\label{eq:pvalue}\n",
    "$$\n",
    "\n",
    "where $f(x';n)$\n",
    "is the $\\chi^2$ distribution for $n$ degrees of freedom. Values can be\n",
    "computed in FIXME ROOT using TMath::ChisquareQuantile(Double_t p, Double_t ndf).\n",
    "\n",
    "The $p-$value for a given $(\\chi^2; n)$ is a measurement of the\n",
    "\"goodness of fit\"; when repeating the experiment several times it gives\n",
    "the probability, under the hypothesis $f$, of obtaining a result as\n",
    "incompatible with $f$ or worse (i.e. $\\chi^2$ equal or larger) than the\n",
    "one actually observed. The threshold on the $p-$value used to reject the\n",
    "model is subjective; typical values used are of a few percent.\n",
    "{numref}`fig-chi2pdg`.\n",
    "maps the relation between the $\\chi^2$, the number of degrees of freedom\n",
    "and the $p-$value. In particular is shown the example where, for $n=4$,\n",
    "a $\\chi^2>6$ will be observed in 20% of the cases.\n",
    "\n",
    "```{figure} ./images/ch7/chi2pdg.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-chi2pdg\n",
    "---\n",
    "FIXME: One minus the\n",
    "$\\chi^2$ cumulative distribution, $1 - F (\\chi^2; n)$, for n degrees of\n",
    "freedom. This gives the $p-$value for the $\\chi^2$ goodness-of-fit\n",
    "test.\n",
    "```\n",
    "\n",
    "We will come back in [Hypothesis Testing](hypothesisTesting.html#hypotheses-and-tests-statistics) to different quantitative ways to\n",
    "evaluate the compatibility of the data with a model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References\n",
    "----------\n",
    "\n",
    "FIXME Most of the material of this section comes from:\n",
    "\n",
    "-   G. Cowan [@CowanBook], \"Statistical Data Analysis\": Ch. 7"
   ]
  }
 ],
 "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
}