Week1/original/CentralLimit/CentralLimit_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Central Limit Theorem (CLT)\n",
    "\n",
    "The aim of this notebook is to illustrate the Central Limit Theorem (CLT) through concrete examples.\n",
    "\n",
    "When you add random numbers from different distributions, but with similar variance (or standard deviation), together and plot the distribution of these sums, you end up with a Gaussian distribution, as dictated by the CLT. \n",
    "The example also illustrates how widths (and therefore uncertainties) are added in quadrature, as one has to divide the sum by the square root of the number of random numbers that went into the sum in order to get a Gaussian of unit width (when using random numbers of unit width, i.e. RMSE $= \\sigma = 1$).\n",
    "\n",
    "\n",
    "## References:\n",
    "- **R. Barlow**: page 49 (and page 45 for Uniform distribution)\n",
    "- **G. Cowan**: page 33\n",
    "- __[Wikipedia: \"Central limit theorem\"](http://en.wikipedia.org/wiki/Central_limit_theorem)__\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   8th of November 2021"
   ]
  },
  {
   "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                                             # Modules to see files and folders in directories"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we set the random seed for the random number generator (RNG). This ensures reproducability (the same results every time the notebook is restarted). "
   ]
  },
  {
   "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 - more on that later.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here we set the parameters for the experiement. We are going to play around with these more.\n",
    "For now we choose to take __10 uniform numbers__ and consider their sum __1000 times__. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_experiments = 1000                  # Number of sums produced\n",
    "N_uniform     = 10                    # Number of uniform numbers used in sum\n",
    "N_exponential = 0                     # Number of exponential numbers used in sum\n",
    "N_cauchy      = 0                     # Number of cauchy numbers used in sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we set up some global constants and bool flags about the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "verbose = True                        # Print some numbers or not?\n",
    "N_verbose = 10                        # If so, how many?\n",
    "save_plots = False                    # Save the plots produced to file(s)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Initial functions\n",
    "We load the external functions from the `AppStat2021/External_Functions` directory. The first line adds the relevant path to the places it will go look for functions, and the second line imports the function wanted. If you get an error here, check that you have copied \"External Functions\" to your directory, and that you refer to the correct path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../../External_Functions')\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax   # Useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop over process:\n",
    "\n",
    "Having loaded everything that we need, we can start the actual program. We start out by initializing a counter to count how many of the produced sums that fall outside some range ($\\pm 3\\sigma$) and some zero-filled numpy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N3_sigma = 0        # Counter for the number of produced sums, that fall outside +-3 sigma\n",
    "\n",
    "x_uniform = np.zeros((N_uniform, N_experiments))\n",
    "x_exponential = np.zeros((N_exponential, N_experiments))\n",
    "x_cauchy = np.zeros((N_cauchy, N_experiments))\n",
    "x_sum = np.zeros((N_experiments))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note on code:\n",
    "\n",
    "We loop over \"experiments\" to make the code more explicit. This could in principle be avoided (which typically also makes the code faster), but for now we keep it, also to be able to comment on what goes on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iexp in range(N_experiments): \n",
    "\n",
    "    if ((iexp+1) % 500 == 0): \n",
    "        print(\"At iexp : \", iexp+1)         # Show progress!\n",
    "    sum_value = 0.0                         # sum_value is the number we are going to add random numbers to!\n",
    "                                            # According to the CLT, it should be Gaussianly distributed.\n",
    " \n",
    "    # Generating uniform numbers (with mean 0, and RMS of 1):\n",
    "    x_uniform[:, iexp] = (np.random.uniform(size=N_uniform) - 0.5) * np.sqrt(12)   # ...sqrt(12). Why? Possibly see references above!\n",
    "    sum_value += x_uniform[:, iexp].sum()               # We add the numbers in \"sum_value\"\n",
    "\n",
    "    # Generating exponential numbers (with mean 0, and RMS of 1):\n",
    "    x_exponential[:, iexp] = np.random.exponential(size=N_exponential) - 1   # Exponential starting at -1. Why?\n",
    "    sum_value += x_exponential[:, iexp].sum()\n",
    "\n",
    "    # Generating numbers according to a Cauchy distribution (1 / (1 + x^2)):\n",
    "    x_cauchy[:, iexp] = np.random.standard_cauchy(size=N_cauchy)\n",
    "    sum_value += x_cauchy[:, iexp].sum()\n",
    "\n",
    "    N_total = N_uniform + N_exponential + N_cauchy\n",
    "    sum_value = sum_value / np.sqrt(N_total)              # Ask yourself, why I divide by sqrt(N)?\n",
    "    x_sum[iexp] = sum_value\n",
    "\n",
    "    # Count the number of times, that the sum is not within 3 sigma of a unit Gaussian:\n",
    "    if not (-3.0 < sum_value < 3.0):\n",
    "        N3_sigma += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have created three 2-dimensional arrays of shape (N_pdf, N_experiments), e.g. (10, 1000). We now flatten the arrays to get 1D arrays in order to see their distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_uniform = x_uniform.flatten()\n",
    "x_exponential = x_exponential.flatten()\n",
    "x_cauchy = x_cauchy.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the input distributions:\n",
    "\n",
    "It is always important to visualize out data to see if our code produced the expected results, find outliers and just generally get a better understand of it.\n",
    "\n",
    "We first define the number of bins and the ranges of the different distributions. Note that it is important to think about this - to many/few bins or poor choice of range may ruin the plot, and one should try to get a \"regular\" bin widths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_bins = 100\n",
    "x_ranges = [(-2.5, 2.5), (-1.5, 8.5), (-5.0, 5.0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now loop over our three distributions and plot them (if they are not empty). In each subplot we plot histograms of the distributions with the given number of bins and ranges as defined above and with the mean, standard deviation and truncated standard deviation for each distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(12,6))\n",
    "\n",
    "x_all = [x_uniform, x_exponential, x_cauchy]\n",
    "titles = ['Uniform', 'Exponential', 'Cauchy']\n",
    "\n",
    "# Here, we simultaneously loop \"smartly\" (using zip) over axis, x, titles, and ranges:\n",
    "for ax_i, x, title, x_range in zip(ax, x_all, titles, x_ranges):    \n",
    "    if len(x) > 0 :\n",
    "        ax_i.hist(x, bins=N_bins, range=x_range, histtype='step')\n",
    "        ymax = ax_i.get_ylim()[1]*1.2                              # Making a little space above the bins\n",
    "        ax_i.set(title=title, ylim=(0, ymax), xlim=x_range)\n",
    "        d = {'Entries': len(x),\n",
    "             'Mean': x.mean(),\n",
    "             'Std': x.std(ddof=1),                 # The truncated Std. includes only values inside range!\n",
    "             'Std_truncated': x[(x_range[0]<x) & (x<x_range[1])].std(ddof=1),\n",
    "        }\n",
    "        text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "        add_text_to_ax(0.1, 0.97, text, ax_i, fontsize=12)\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "if save_plots:\n",
    "    fig.savefig('Plot_CentralLimit_Input.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Outside $\\pm 3 \\sigma$\n",
    "\n",
    "Below we write how many of the experiments had a result outside the range [-3,3], i.e. beyond $\\pm 3 \\sigma$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Number of experiments beyond 3 sigma / total: {N3_sigma:4d} / {N_experiments:d} = {N3_sigma/N_experiments:6.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- __Question__: *Do the numbers above make sense? What did you expect it to be?*\n",
    "***\n",
    "## Draw output plots with corresponding fits to the screen:\n",
    "\n",
    "Having visualized the individual distributions above, we now visualize the sum of the different distributions, `x_sum`. We first plot a histogram of it in the range from `xmin` to `xmax`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "xmin, xmax = -5, 5\n",
    "\n",
    "fig2, ax2 = plt.subplots(figsize=(16, 7)) \n",
    "hist2 = ax2.hist(x_sum, bins=N_bins, range=(xmin, xmax), histtype='step', linewidth=2)\n",
    "ax2.set(xlabel='Sum', ylabel='Frequency / 0.1', title='Histogram of x_sum (sum of input values)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we overlay the plot with a unit Gaussian (i.e. not fitted to the data), to see if the resulting sums actually distribute themselves Gaussianly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define your PDF / model \n",
    "def gauss_pdf(x, mu, sigma):\n",
    "    \"\"\"Normalized Gaussian\"\"\"\n",
    "    return 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-(x - mu) ** 2 / 2. / sigma ** 2)\n",
    "\n",
    "def gauss_extended(x, N, mu, sigma):\n",
    "    \"\"\"Non-normalized Gaussian\"\"\"\n",
    "    return N * gauss_pdf(x, mu, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "N_scale = (xmax-xmin) / N_bins            # The scale factor between histogram and the fit. Takes e.g. bin width into account.\n",
    "x_gauss = np.linspace(xmin, xmax, 1000)   # Create the x-axis for the plot of the fitted function\n",
    "y_gauss = N_scale*gauss_extended(x_gauss, len(x_sum), 0, 1)                   # Unit Gaussian\n",
    "ax2.plot(x_gauss, y_gauss, '-', color='blue', label='Unit Gauss (no fit)') \n",
    "ax2.legend(loc='upper right', fontsize=20)\n",
    "fig2.tight_layout()\n",
    "\n",
    "# Note, we refer to the old plot \"ax2\". Had we done another plot in between, \n",
    "# we would not have been able to plot on top of the old figure with the matlab syntax.\n",
    "if save_plots :\n",
    "    fig2.savefig('Histogram.pdf', dpi=600)\n",
    "fig2        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "First make sure that you understand what the Central Limit Theorem (CLT) states! Then, acquaint yourself with the program. Make sure that you read through it, as many\n",
    "of these features will be used onwards. Do you understand why the uniform distribution needs to go from $\\pm \\sqrt 3$ in order to give a distribution with a width of $1$ (i.e. unit) and why you subtract $1$ from the exponential distribution (and how this works at all)?\n",
    "\n",
    "# Questions:\n",
    "\n",
    "1. What is the mean and RMSE (i.e. Standard Deviation) of the input distributions?\n",
    "\n",
    "2. Why is there a $\\frac{1}{\\sqrt N}$ at the bottom of cell [7] in the line `sum_value = sum_value / np.sqrt(Ntotal)` (when summing up the various contributions to sum)?\n",
    "    Hint: Assume that I always wanted to compare the distribution of sums with a UNIT Gaussian.\n",
    "\n",
    "3. Using a sum of 10 uniform random numbers with mean 0 and width 1, what is the expected \n",
    "    width of the resulting distribution according to CLT? What is the probability of\n",
    "    obtaining a number beyond 3 sigma, i.e. how many numbers did you get beyond 3 sigma?\n",
    "    What would you expect from a true Gaussian distribution?\n",
    "    And what about the same question for 3.5 sigma? And 4.0 sigma?\n",
    "    Put additional counters into the code, to test any effects in the tails (which have little statistics as it is).\n",
    "    Increase the number of experiments run to (much) more than 1000...\n",
    "\n",
    "4. Now try to add 10 exponential (i.e. set `N_exponential=10` and rerun the program). Does that give something Gaussian? What about 1000?\n",
    "    Then try to add 10 cauchy numbers (i.e. set `N_cauchy=10` and rerun the program). Does that give something Gaussian? What about 1000?\n",
    "    If not Gaussian, why do the Cauchy distribution \"ruin\" the Gaussian distribution?\n",
    "    And is this in conflict with the Central Limit Theorem?\n",
    "\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    "5. If one used a trunkated mean throwing away the top and bottom e.g. 10%,\n",
    "    will the truncated mean of 1000 Cauchy numbers then converge to a Gaussian (possibly not with unit width)?\n",
    "\n",
    "6. What is the lowest number of uniform random numbers needs to be added, before the probability\n",
    "    for the sum to follow a Gaussian distribution is greater than 1% (on average)\n",
    "    when considering 1000 sums (i.e. `N_experiments=1000`)? Here, a $\\chi^2$ fit is needed.\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.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}