Week4/original/HypothesisTesting/HypothesisTesting_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing\n",
    "### Authors: \n",
    "- Christian Michelsen (Niels Bohr Institute)\n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 03-12-2019 (latest update)\n",
    "\n",
    "***\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 an unbinned KS test, and compares it to a $\\chi^2$-test and a simple comparison of means.\n",
    "\n",
    "The distributions compared are two unit Gaussians, 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",
    "***\n",
    "\n",
    "First, we import the modules we want to use:\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 an 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",
    "    \n",
    "    d = {'Entries': len(x), \n",
    "         'Mean': x.mean(),\n",
    "         'STD': x.std(ddof=1),\n",
    "        }\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:"
   ]
  },
  {
   "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:"
   ]
  },
  {
   "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 and error on mean:\n",
    "    mean_A, width_A, sdom_A = mean_std_sdom(x_A_array) \n",
    "    mean_B, width_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.05')        \n",
    "        ax_text(x_A_array, ax1, 0.04, 0.85, '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.65, '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 fit 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='Histogram, probability mu', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\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='Histogram, probability chi2', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\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='Histogram, probability Kolmogorov', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\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$? And do 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 you guess/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. 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",
    "5. Obviously, the test of the means is not sensitive the 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",
    "6. 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. :-)"
   ]
  }
 ],
 "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
}