Week4/HypothesisTesting/HypothesisTesting_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing\n",
    "\n",
    "Python notebook for illustrating the concept of Hypothesis Testing and specific test statistics; among them the very useful Kolmogorov-Smirnov test.\n",
    "\n",
    "The Kolmogorov-Smirnov test (KS test) is a general test to evaluate if two distributions in 1D are the same. This program applies three hypothesis tests to determine, if two distributions are the same:\n",
    "1. A simple comparison of means.\n",
    "2. A $\\chi^2$-test between histograms.\n",
    "3. An unbinned KS test\n",
    "\n",
    "The distributions compared are two unit Gaussians (A and B), where one is then modified by changing:\n",
    "- Mean\n",
    "- Width\n",
    "- Normalisation\n",
    "\n",
    "The sensitivity of each test is then considered for each of these changes.\n",
    "\n",
    "### References:\n",
    "- Barlow: p. 155-156\n",
    "- __[Wikipedia: Kolmogorov-Smirnov test](http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test)__\n",
    "- Though influenced by biostatistics, a good discussion of p-values and their distribution can be found here:\n",
    "  [How to interpret a p-value histogram?](http://varianceexplained.org/statistics/interpreting-pvalue-histogram/)\n",
    "\n",
    "### Authors: \n",
    "Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "07-12-2022 (latest update)\n",
    "\n",
    "***"
   ]
  },
  {
   "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.special import erfc\n",
    "from scipy import stats\n",
    "\n",
    "sys.path.append('../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression, BinnedLH, UnbinnedLH\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters of the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.random             # Random generator\n",
    "r.seed(42)                # Set a random seed (but a fixed one)\n",
    "\n",
    "save_plots = False\n",
    "verbose = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The small function below is just a simple helper function that takes a 1D-array input along with an axis, position, and color arguments, and plot the number of entries, the mean, and the standard deviation on the axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ax_text(x, ax, posx, posy, color='k'):\n",
    "    d = {'Entries': len(x), \n",
    "         'Mean': x.mean(),\n",
    "         'STD': x.std(ddof=1),\n",
    "        }    \n",
    "    add_text_to_ax(posx, posy, nice_string_output(d), ax, fontsize=12, color=color)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and finally a function that calculates the mean, standard deviation and the standard deviation (i.e. uncertainty) on mean (sdom):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_std_sdom(x):\n",
    "    std = np.std(x, ddof=1)\n",
    "    return np.mean(x), std, std / np.sqrt(len(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the experiment:\n",
    "\n",
    "How many experiments, and how many events in each:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_exp = 1\n",
    "N_events_A = 100\n",
    "N_events_B = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the two Gaussians to be generated (no difference to begin with!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_mean_A  =  0.0\n",
    "dist_width_A =  1.0\n",
    "dist_mean_B  =  0.0\n",
    "dist_width_B =  1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "Define the number of bins and the range, initialize empty arrays to store the results in and make an empty figure (to be filled in later):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "N_bins = 100\n",
    "xmin, xmax = -5.0, 5.0\n",
    "\n",
    "all_p_mean = np.zeros(N_exp)\n",
    "all_p_chi2 = np.zeros(N_exp)\n",
    "all_p_ks   = np.zeros(N_exp)\n",
    "\n",
    "# Figure for the two distributions, A and B, in the first experiment:\n",
    "fig1, ax1 = plt.subplots(figsize=(10, 6))\n",
    "plt.close(fig1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop over how many times we want to run the experiment, and for each calculate the p-value of the two distributions coming from the same underlying PDF (put in calculations yourself):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iexp in range(N_exp):\n",
    "    if ((iexp+1)%1000 == 0):\n",
    "        print(f\"Got to experiment number: {iexp+1}\")\n",
    "\n",
    "    # Generate data:\n",
    "    x_A_array = r.normal(dist_mean_A, dist_width_A, N_events_A)\n",
    "    x_B_array = r.normal(dist_mean_B, dist_width_B, N_events_B)\n",
    "    \n",
    "    \n",
    "    # Test if there is a difference in the mean:\n",
    "    # ------------------------------------------\n",
    "    # Calculate mean, standard deviation, and error on mean (sdom):\n",
    "    mean_A, std_A, sdom_A = mean_std_sdom(x_A_array)\n",
    "    mean_B, std_B, sdom_B = mean_std_sdom(x_B_array)\n",
    "\n",
    "    # Consider the difference between means in terms of the uncertainty:\n",
    "    d_mean = mean_A - mean_B\n",
    "    # ... how many sigmas is that away?\n",
    "\n",
    "    # Turn a number of sigmas into a probability (i.e. p-value):\n",
    "    p_mean    = 0.5      # Calculate yourself. HINT: \"stats.norm.cdf or stats.norm.sf may be useful!\"\n",
    "    all_p_mean[iexp] = p_mean\n",
    "    \n",
    "    \n",
    "    # Test if there is a difference with the chi2:\n",
    "    # --------------------------------------------\n",
    "    # Chi2 Test:\n",
    "    p_chi2 = 0.5         # Calculate the p-value of the Chi2 between histograms of A and B yourself.\n",
    "    all_p_chi2[iexp] = p_chi2\n",
    "\n",
    "    \n",
    "    # Test if there is a difference with the Kolmogorov-Smirnov test on arrays (i.e. unbinned):\n",
    "    # -----------------------------------------------------------------------------------------\n",
    "    p_ks = stats.ks_2samp(x_A_array, x_B_array)[1]    # Fortunately, the K-S test is implemented in stats!\n",
    "    all_p_ks[iexp] = p_ks\n",
    "\n",
    "\n",
    "    # Print the results for the first 10 experiments\n",
    "    if (verbose and iexp < 10) :\n",
    "      print(f\"{iexp:4d}:  p_mean: {p_mean:7.5f}   p_chi2: {p_chi2:7.5f}   p_ks: {p_ks:7.5f}\")\n",
    "\n",
    "    \n",
    "    # In case one wants to plot the distribution for visual inspection:\n",
    "    if (iexp == 0):\n",
    "        \n",
    "        ax1.hist(x_A_array, N_bins, (xmin, xmax), histtype='step', label='A', color='blue')\n",
    "        ax1.set(title='Histograms of A and B', xlabel='A / B', ylabel='Frequency / 0.1')        \n",
    "        ax_text(x_A_array, ax1, 0.04, 0.95, 'blue')\n",
    "\n",
    "        ax1.hist(x_B_array, N_bins, (xmin, xmax), histtype='step', label='B', color='red')\n",
    "        ax_text(x_B_array, ax1, 0.04, 0.75, 'red')\n",
    "        \n",
    "        ax1.legend()\n",
    "        fig1.tight_layout()\n",
    "\n",
    "        \n",
    "fig1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Show the distribution of hypothesis test p-values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_bins = 50\n",
    "\n",
    "if (N_exp > 1):\n",
    "    fig2, ax2 = plt.subplots(nrows=3, figsize=(12, 14))\n",
    "    \n",
    "    ax2[0].hist(all_p_mean, N_bins, (0, 1), histtype='step')\n",
    "    ax2[0].set(title='Probability of means', xlabel='p-value', ylabel='Frequency / 0.02')\n",
    "    ax_text(all_p_mean, ax2[0], 0.04, 0.25)\n",
    "    \n",
    "\n",
    "    ax2[1].hist(all_p_chi2, N_bins, (0, 1), histtype='step')\n",
    "    ax2[1].set(title='Probability of chi2', xlabel='p-value', ylabel='Frequency / 0.02')\n",
    "    ax_text(all_p_chi2, ax2[1], 0.04, 0.25)\n",
    "    \n",
    "    ax2[2].hist(all_p_ks, N_bins, (0, 1), histtype='step')\n",
    "    ax2[2].set(title='Probability of Kolmogorov', xlabel='p-value', ylabel='Frequency / 0.02')\n",
    "    ax_text(all_p_ks, ax2[2], 0.04, 0.25)\n",
    "\n",
    "    fig2.tight_layout()\n",
    "\n",
    "\n",
    "    if save_plots:\n",
    "        fig2.savefig('PvalueDists.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions:\n",
    "\n",
    "1. First run the program with one experiment (N_exp = 1) to display the two distributions A and B, when:\n",
    "    - They are the same.\n",
    "    - The mean of A is increased (to e.g. 0.1).\n",
    "    - The width of A is enlarged (to e.g. 1.2).\n",
    "    - The normalisation of A is increased.\n",
    "\n",
    "Get a feel for how much you need to change the distribution, before you can _by eye_ see that they are not the same. I.e. can you see any difference, if `mean_A` $= 0.1$? Or how about $0.2$? How do you quantify this and when do you start to doubt? And how about `width_A` $= 1.1$? Or $1.2$? Again, can you see it by eye? Finally, try to put $1050$ events into B. Is that visible? How about $1100$?<br>\n",
    "When you have an idea about when you can see effects by eye, now check if you see an impact in the p-values?\n",
    "\n",
    "2. Could you for the test of the means have calculated how much of a change in the mean is needed for a difference to be statistically significant? Do so, and see if it somewhat matches your calculation/estimate from above!\n",
    "\n",
    "\n",
    "3. Now run the tests 1000 times, where A and B are unit Gaussians and thus identical. How should the distributions of the test probabilities come out? And is this the case, approximately? If not, think of reasons for this, and what could be a remedy. HINT: Large statistics is always easier!\n",
    "\n",
    "\n",
    "4. Repeat the changes in question 1), and see which tests \"reacts\" most to these modifications. How much of a change in the mean is required for 95% of the tests (of each kind) to give a probability below 5%? How much is required for the width? And the norm?\n",
    "\n",
    "\n",
    "5. Possibly try to test different distributions than the Gaussian one (e.g. exponential, uniform, etc.), and see how the tests performs.\n",
    "\n",
    "\n",
    "NOTE: The Kolmogorov-Smirnov test has the great advantage that it can handle ANY distribution (even the Cauchy distribution - remind yourself of that one!). The reason is, that it doesn't care about any PDF, nor how far out an outlier is. It is just a matter of the difference in integrals between the two functions.\n",
    "\n",
    "\n",
    "## Advanced:\n",
    "\n",
    "6. Obviously, the test of the means is not sensitive to a change in the width. Make such a test yourself by calculating the widths and the uncertainty on the widths (or perhaps try the F-test!). Note that in a (unit) Gaussian the uncertainty on the width is of the same order as that of the means!\n",
    "\n",
    "\n",
    "## Very advanced:\n",
    "7. Implement in python the following tests:\n",
    "     - Lilliefors test\n",
    "     - Shapiro-Wilk test\n",
    "     - Anderson-Darling test\n",
    "     - Cramer-von-Mises test\n",
    "     - Jarque-Bera test\n",
    "     - Kuiper's test\n",
    "     - Mann-Whitney-Wilcoxon test\n",
    "     - Siegel-Tukey test\n",
    "     \n",
    "and quantify under various conditions and datasets the power of each and the correlation among them. Write it up, and send it to a statistics journal. :-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points:\n",
    "\n",
    "The exercise should give you a good idea about the process of hypothesis testing, namely that you calculate a p-value, which under the null hypothesis is uniformly distributed.\n",
    "\n",
    "You should be aware that in this exercise, we \"play Gods\" and repeat a perfect experiment over and over, to see the changes in the p-value distribution, when we change one of the input distributions. However, when you have real data, you will only get a single p-value! And you should be in a position, where you trust (and have perhaps tested using simulation) that the p-value really is uniform for the null hypothesis, and thus points to something \"non-null\", if very low.\n",
    "\n",
    "The three tests (means, chi2, and KS) are some of the very standard approaches, but there are certainly also others."
   ]
  }
 ],
 "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.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}