Week3/original/CoffeeUsage/CoffeeUsage_empty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coffee Usage\n",
    "\n",
    "Python notebook for analysing the coffee usage in the NBI HEP group, and determining the length of the Christmas vacation from a fit to it. The exercise lies in writing a stepwise fit function, and making it work on real data.\n",
    "\n",
    "For a period in 2009-2010, the usage of the old coffey machine in the NBI HEP group was (somewhat irregularly) monitored. Below is the count of total number of cups of coffey ever brewed at given dates. We will mainly focus on the data from around Christmas...\n",
    "\n",
    "\n",
    "- 28479   4/11-2009     NOTE: This day, we in the following define as day 0! \n",
    "- 28674  13/11-2009 \n",
    "- 28777  18/11-2009 \n",
    "- 28964  25/11-2009 \n",
    "- 29041  27/11-2009 \n",
    "- 29374  10/12-2009 \n",
    "- ~29650   8/ 1-2010\n",
    "- 30001  29/ 1-2010 (?)\n",
    "- 30221   8/ 2-2010\n",
    "- 30498  21/ 2-2010\n",
    "- 32412  17/ 5-2010\n",
    "- 33676  11/ 8-2010\n",
    "- 34008   9/ 9-2010`\n",
    "\n",
    "### Author: \n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 29-11-2021 (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                                             # Modules to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the external functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": [
    "and define the parameters of the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting functions:\n",
    "\n",
    "First we define a simple linear function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_coffee_simple(x, p0, p1): \n",
    "    return p0 + p1*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now YOU should add an extra parameter, which is the length of the vacation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "def func_coffee_advanced(x, p0, p1, p2):\n",
    "\n",
    "    # Here goes your function in terms of p0, p1, and p2 (vacation length)...\n",
    "    return p0 + p1*x + p2*x*x/100.0       # Just a silly function!\n",
    "    \n",
    "# Numpy-way of \"vectorizing\" the function above:\n",
    "func_coffee_advanced_vec = np.vectorize(func_coffee_advanced)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why we are using the last line is to allow `func_coffee_advanced` to take Numpy arrays as inputs, so we do not have to loop over entire lists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coffee Usage\n",
    "\n",
    "Data set (So small that we will not use a seperate file for it!). We define 4th of November 2009 to be day 0, and count from there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "days = np.array([     0,     9,    14,    21,    23,    36,    65,    76,    86,    99,   194,   280,   309 ])\n",
    "cups = np.array([ 28479, 28674, 28777, 28964, 29041, 29374, 29650, 30001, 30221, 30498, 32412, 33676, 34008 ])\n",
    "edays = np.zeros_like(days)\n",
    "ecups = np.ones_like(cups)*30.0   # I estimate the uncertainty to be 30 cups, but perhaps you disagree?\n",
    "                                  # Certainly, it is NOT the square root, as it is not at all a Poisson number!\n",
    "\n",
    "# Print the number of days:\n",
    "for day, cup in zip(days, cups):\n",
    "    print(f\"days: {day:3.0f}    cups: {cup:5.0f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit and plot graph on screen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Define the relevant range:\n",
    "xmin, xmax = 60, 105\n",
    "mask_range = (xmin < days) & (days < xmax)\n",
    "\n",
    "# Fit the data points in the range (xmin, xmax)\n",
    "chi2_object_coffee = Chi2Regression(func_coffee_simple, days[mask_range], cups[mask_range], ecups[mask_range])\n",
    "minuit_coffee = Minuit(chi2_object_coffee, p0=28000.0, p1=20.0)\n",
    "minuit_coffee.migrad();\n",
    "p0, p1 = minuit_coffee.values[:]\n",
    "\n",
    "print(\"Simple fit\")\n",
    "for name in minuit_coffee.parameters:\n",
    "    print(\"Fit value: {0} = {1:.5f} +/- {2:.5f}\".format(name, minuit_coffee.values[name], minuit_coffee.errors[name]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data along with the fit on top of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data:\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.errorbar(days, cups, xerr=edays, yerr=ecups, fmt='o', label='Data')\n",
    "ax.set(xlim=(-20, 400), xlabel='Days', ylabel='Cups', title='Coffe usage as function of time')\n",
    "\n",
    "# Plot the fit on top:\n",
    "x_fit = np.linspace(xmin, xmax, 1000)\n",
    "y_fit = func_coffee_simple(x_fit, *minuit_coffee.values[:])\n",
    "ax.plot(x_fit, y_fit, '-', label='Simple Fit')\n",
    "\n",
    "# Get the ChiSquare probability:\n",
    "chi2_lin = minuit_coffee.fval\n",
    "ndof_lin = len(days[mask_range]) - len(minuit_coffee.values[:])\n",
    "chi2_prob_lin = stats.chi2.sf(chi2_lin, ndof_lin)\n",
    "\n",
    "# Include fit results in the plot:\n",
    "d = {'Chi2': chi2_lin,\n",
    "     'Ndof': ndof_lin,\n",
    "     'Prob': chi2_prob_lin,\n",
    "     'p0': [minuit_coffee.values['p0'], minuit_coffee.errors['p0']],\n",
    "     'p1': [minuit_coffee.values['p1'], minuit_coffee.errors['p1']],\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.97, text, ax, fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced fit:\n",
    "Try to fit the advanced function in the range (-10, 110):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions:\n",
    "\n",
    "First, only consider the first 100 days of data.\n",
    "\n",
    "0. From the ChiSquare of the linear fit to four points (65-100 days), is the estimated uncertainty of 30 cups reasonable?\n",
    "\n",
    "1. Assuming the error of 30 cups, do the numbers follow the hypothesis of constant use? Quantify this, and find out how large the error has to be (roughly), for this hypothesis to be credible.\n",
    "\n",
    "2. Does taking into account Christmas vacation looks like it will improve the above hypothesis? Can you actually fit the length of vacation? Try to rewrite the function \"func_coffee_advanced\", such that it includes \"if\", dividing the function into two linear functions with the same slope, that has an interval with no usage (i.e. the Christmas vacation). What length do you obtain? And how much better is the model?\n",
    "\n",
    "Now consider the whole data set.\n",
    "\n",
    "3. The total number of cups of coffey ever brewed was 36716, after which the old coffey machine was decommissioned. From the above data, estimate when this happened (including error!). And when would you estimate that the coffey machine was commissioned originally?"
   ]
  }
 ],
 "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
}