Week1/original/ChiSquareTest/ChiSquareTest_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ChiSquare Test:\n",
    "This Python macro illustrates the use of ChiSquare as a goodness-of-fit measure,\n",
    "how the ChiSquare distribution comes about, how to use the ChiSquare value (along with\n",
    "number of degrees of freedom, Ndof) as a test, and that it actually works!\n",
    "\n",
    "The program generates a certain number of datasets, each consisting of 9 points along\n",
    "a line. These are then fitted, and the result and the Chi2 of the fit is recorded along\n",
    "with the probability of the fit.\n",
    "\n",
    "## References:\n",
    "- Barlow: Chapter 6\n",
    "- Cowan: Chapter 2.7, Chapter 7\n",
    "- Bevington: Chapter 6\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   8th of November 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial Question/Task:\n",
    "Make sure you've read the relevant references and that you understand not only what the ChiSquare is, but also that it follows the ChiSquare distribution, and that the probability of obtaining such a ChiSquare or worse (assuming that the hypothesis/PDF is the correct one!) can be calculated from it. Try to go to the \"External Functions\" and find the Chi2Regression function, and then read through the code to see, that it does as you would expect (note: It is essentially the last line, in which the Chi2 calculation is done).\n",
    "\n",
    "Executive Summary: The ChiSquare method is used both for getting the best fit, and for testing the quality/probability of this fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "r = np.random                         # Random generator\n",
    "r.seed(1)                            # Set a random seed (but a fixed one - more on that later.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nexp = 1\n",
    "Npoints = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters used in this problem (constant, linear coef., quadratic coef., and uncertainty in y) for generating a line (parabola) with N points and given uncertainty on points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha0 = 3.6\n",
    "alpha1 = 0.3\n",
    "alpha2 = 0.0\n",
    "sigmay = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating and fitting data:\n",
    "In the following, we generate data points following a simple line. The linear fit (and only this) can be done analytically, as discussed in Barlow chapter 6.2 and outlined here: http://www.nbi.dk/~petersen/Teaching/Stat2019/StraightLineFit.pdf .\n",
    "The numerical fit of the line (and any other function) is done iteratively by Minuit. The code is slightly shorter, but a lot slower, as the code will typically test many (50-5000) possible combination of fit parameters, before it finds the best values for the fitting parameters!\n",
    "\n",
    "The whole thing is contained in a loop, allowing us to repeat the generation and fitting, so that we can see the result on more than just a case-by-case. For this reason, we save the result of the fit (in pre-defined arrays)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arrays for storing fit results:\n",
    "array_alpha0 = np.zeros(Nexp)\n",
    "array_alpha1 = np.zeros(Nexp)\n",
    "array_Chi2 = np.zeros(Nexp)\n",
    "array_Prob = np.zeros(Nexp)\n",
    "\n",
    "# Loop, repeating the data generation and fit:\n",
    "for iexp in range(Nexp) : \n",
    "\n",
    "    # Generating data by filling values into (x,y) and associated uncertainties (here chosen to be 0 for x):\n",
    "    x = np.arange(Npoints)+1\n",
    "    ex = np.zeros_like(x)\n",
    "    y = alpha0 + alpha1*x + r.normal(0, sigmay, Npoints)   # + alpha2*x**2\n",
    "    ey = sigmay*np.ones_like(x)\n",
    "\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Numerical fit:\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Define a fit function:\n",
    "    def fit_function(x, alpha0, alpha1):\n",
    "        return alpha0 + alpha1*x\n",
    "\n",
    "    # Now we define a ChiSquare to be minimised (using probfit), where we set various settings and starting parameters:\n",
    "    chi2_object = Chi2Regression(fit_function, x, y, ey) \n",
    "    minuit = Minuit(chi2_object, pedantic=False, alpha0=3.0, alpha1=0.0, print_level=0)  \n",
    "    minuit.migrad();                                                  # Perform the actual fit\n",
    "    minuit_output = [minuit.get_fmin(), minuit.get_param_states()]    # Save the output parameters in case needed\n",
    "    \n",
    "    # Here we extract the fitting parameters and their errors\n",
    "    alpha0_fit = minuit.values['alpha0']\n",
    "    alpha1_fit = minuit.values['alpha1']\n",
    "    sigma_alpha0_fit = minuit.errors['alpha0']\n",
    "    sigma_alpha1_fit = minuit.errors['alpha1']\n",
    "    \n",
    "    Nvar = 2                     # Number of variables (alpha0 and alpha1)\n",
    "    Ndof_fit = Npoints - Nvar    # Number of degrees of freedom = Number of data points - Number of variables\n",
    "    \n",
    "    # In Minuit, you can just ask the fit function for the Chi2:\n",
    "    Chi2_fit = minuit.fval                          # The chi2 value\n",
    "    Prob_fit = stats.chi2.sf(Chi2_fit, Ndof_fit)    # The chi2 probability given N degrees of freedom\n",
    "    \n",
    "    # Fill the arrays with fit results (to produce plots of these at the end):\n",
    "    array_alpha0[iexp] = alpha0_fit\n",
    "    array_alpha1[iexp] = alpha1_fit\n",
    "    array_Chi2[iexp] = Chi2_fit\n",
    "    array_Prob[iexp] = Prob_fit\n",
    "    \n",
    "    # Let us see what the fit gives for the first couple of data sets:\n",
    "    if (iexp < 10) :\n",
    "        print(f\"  Fit: a0={alpha0_fit:6.3f}+-{sigma_alpha0_fit:5.3f}  a1={alpha1_fit:5.3f}+-{sigma_alpha1_fit:5.3f}  p={Prob_fit:6.4f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit the data and plot the result:\n",
    "Below we plot the latest fit (i.e. the last one from above), so that we can actually see what is going on (and potentially wrong!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.errorbar(x, y, ey, fmt='ro', ecolor='k', elinewidth=1, capsize=2, capthick=1)\n",
    "ax.plot(x, fit_function(x, *minuit.args), '-r')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = {'alpha0':   [alpha0_fit, sigma_alpha0_fit],\n",
    "     'alpha1':   [alpha1_fit, sigma_alpha1_fit],\n",
    "     'Chi2':     Chi2_fit,\n",
    "     'ndf':      Ndof_fit,\n",
    "     'Prob':     Prob_fit,\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.95, text, ax, fontsize=14)\n",
    "fig.tight_layout()\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (save_plots) :\n",
    "    fig.savefig(\"FitToLine.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Many experiments\n",
    "In case we have more than one \"experiment\" (i.e. fit), we would like to see how the fit results and the Chi2 values distribute themselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if (Nexp > 1) :\n",
    "    \n",
    "    fig2, ax2 = plt.subplots(nrows=2, ncols=2, figsize=(14,10))\n",
    "    ax2 = ax2.flatten()\n",
    "\n",
    "    # Distribution of alpha0 values:\n",
    "    ax2[0].hist(array_alpha0, bins=80, range=(2.8, 4.4), histtype='step')\n",
    "    ax2[0].set_title('Histogram of alpha0')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_alpha0), array_alpha0.mean(), array_alpha0.std(ddof=1))\n",
    "    ax2[0].text(0.62, 0.95, string, family='monospace', transform=ax2[0].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution of alpha1 values:\n",
    "    ax2[1].hist(array_alpha1, bins=60, range=(-0.05, 0.25), histtype='step')\n",
    "    ax2[1].set_title('Histogram of alpha1')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_alpha1), array_alpha1.mean(), array_alpha1.std(ddof=1))\n",
    "    ax2[1].text(0.05, 0.95, string, family='monospace', transform=ax2[1].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution (normed) of Chi2 values, along with Chi2-distribution (for Ndof=7) drawn on top:\n",
    "    ax2[2].hist(array_Chi2, bins=25, range=(0.0, 25.0), histtype='step', density=True)\n",
    "    ax2[2].set_title('Histogram of Chi2')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_Chi2), array_Chi2.mean(), array_Chi2.std(ddof=1))\n",
    "    ax2[2].text(0.62, 0.95, string, family='monospace', transform=ax2[2].transAxes, fontsize=12, verticalalignment='top')\n",
    "    x_axis = np.linspace(0, 20, 1000)\n",
    "    y_axis = stats.chi2.pdf(x_axis, df=7)\n",
    "    ax2[2].plot(x_axis, y_axis, 'r-')\n",
    "    ax2[2].text(0.45, 0.65, \"Note: This is not a fit!\", transform=ax2[2].transAxes, fontsize=16, verticalalignment='top')\n",
    "\n",
    "    # Distribution of Chi2 fit probabilities:\n",
    "    ax2[3].hist(array_Prob, bins=20, range=(0.0, 1.0), histtype='step')\n",
    "    ax2[3].set_title('Histogram of Prob')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_Prob), array_Prob.mean(), array_Prob.std(ddof=1))\n",
    "    ax2[3].text(0.3, 0.2, string, family='monospace', transform=ax2[3].transAxes, fontsize=12, verticalalignment='top')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    if (save_plots) :\n",
    "        fig2.savefig('FitResultDistributions.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remark on $\\chi^2$ rule of thumb:\n",
    "A general rule is that the \"reduced chisquare\" (defined as $\\chi^2$/Ndof) should be around 1. That is true, but the rule definitely has its limitations, since \"around\" depends on the size of Ndof! This is examplified below, considering two cases with the same reduced chisquare of 1.5, but with **very** different conclusions.\n",
    "\n",
    "The reason for this can be found in the mean and width of the $\\chi^2$ distribution (with Ndof = k), which are k and sqrt(2k). From the central limit theorem, the $\\chi^2$ should approach a Gaussian for \"large\" values of k (k > 50 is probably enough), but the **relative** width (i.e. width/mean) becomes smaller and smaller as k grows, and a reduced chisquare of e.g. 1.5 becomes more and more unlikely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Numbers to illustrate the point, that using the ratio Chi2/Ndof (\"reduced chisquare\") is only a rule of thumb:\n",
    "print(\"\\n  Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:\")\n",
    "print(\"  Low number of degrees of freedom:   Prob(  3.0,  2) = %9.7f\"%(stats.chi2.sf(  3.0,   2)))\n",
    "print(\"  High number of degrees of freedom:  Prob(300.0,200) = %9.7f\"%(stats.chi2.sf(300.0, 200)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "\n",
    "\n",
    "# Questions:\n",
    "\n",
    " 1) Run the code such that it does exactly one fit (initial setting), and take a look at the line fitted.\n",
    "    Does this look reasonable, and are you satisfied with the fit probability? What is the chance that the\n",
    "    input for the data could actually be from a flat distribution, i.e. without the slope?\n",
    "    [There is a correct (slower) way of doing this, and another quick but approximate way of testing this].\n",
    "\n",
    " 2) Does the fit reproduce the input numbers well (include the uncertainties in your\n",
    "    answer)?\n",
    "\n",
    " 3) Now increase the number of experiments to e.g. 1000, and rerun the macro. Figure\n",
    "    out what you see in the window split in 2-by-2, and go through each of these to\n",
    "    see, if you understand every feature of the distributions shown, and if you are\n",
    "    happy with them! Is this what you expected? Are they also as wide as you would expect?\n",
    "    Possibly try to make some fits of these output distributions to check, if you're satisfied!\n",
    "    <br>\n",
    "    This somehow makes this the \"long\" question without any leading questions.\n",
    "    But you should try to be statistically minded enough and have an idea of what to\n",
    "    look for, at least to some degree. We will discuss this further in class.\n",
    "\n",
    " 4) Investigate if the distributions of probabilities is flat, or if it has some\n",
    "    slope to it. Do you understand why the distributions of probabilities is flat?\n",
    "    (This is typically conceptually hard to begin with, but don't worry!)\n",
    "\n",
    " 5) Find the line of code where the random points for the line are generated, and add a\n",
    "    quadratic term in x with the coefficient alpha2 = -0.04. Run the program again with Nexp=1,\n",
    "    and start by looking at the single fit in the graph. Can you see this change?\n",
    "    Now run 1000 experiments again, and see what has changed in the distributions\n",
    "    in the 2-by-2 window when including such a term. Are these changes as you would\n",
    "    have expected them?\n",
    "    <br>\n",
    "    NOTE: The reason why the quadratic term is more clear in the latter case is NOT\n",
    "    that it is a better method. It is simply because there is 1000 times more statistics!\n",
    "\n",
    " 6) With the quadratic coefficient from question 5), the linear fit of course does\n",
    "    not do very well. What changes are needed for the fit to be good again? Make\n",
    "    these changes, and see that the condition in question 4) is again met.\n",
    "\n",
    "### (Semi-)Advanced questions:\n",
    "\n",
    " 7) Try to generate data according to more advanced functions, and then fit the data with these.\n",
    "    Do you still obtain a flat distribution in probability? Do you also manage to fit these with\n",
    "    simpler/other functions and still maintain the flat distribution?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}