Week1/original/CentralLimit/CentralLimit_empty.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:   12th of November 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "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 later on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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, like $\\pi$, and the bool flags about the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "pi = 3.14159265358979323846264338328  # Selfexplanatory!!!\n",
    "pi = np.pi                            # Another way of doing it - surely better!\n",
    "\n",
    "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 `AppStat2020/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": 5,
   "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": 6,
   "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": [
    "### Problem:\n",
    "\n",
    "Fill the above numpy arrays with numbers according to the mentioned distributions, ensuring that the uniform and exponential have mean=0 and SD=1. Let the sum be the sum of all these numbers (possibly a mix), and consider the distribution of the sums. Do they look Gaussian?\n",
    "In order to compare the distribution of sums with a UNIT Gaussian, you might want to divide the sum values by sqrt(N), where N is the number of input values. As a part of checking, count how many values are outside 3 sigma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw 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": "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 RMS 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 of 10 Cauchy numbers (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. How few/many 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 using 1000 sums (i.e. `N_experiments=1000`)? Here, a $\\chi^2$ fit is needed.\n",
    "   "
   ]
  }
 ],
 "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
}