Week5/original/Calibration/Calibration_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibration Exercice\n",
    "\n",
    "This exercise is about calibration, that is to make a certain type of method or measurement more precise by taking into account effects that otherwise blurs the precision of the method by correcting for these effects. The below example is constructed/invented, but outlines the considerations.\n",
    "\n",
    "## Description\n",
    "  You're developping a new method for measuring the distance to stars, and want to\n",
    "  calibrate and thus improve this method, such that the precision obtained is unbiased\n",
    "  and has a minimal variance. You know that the method depends on several factors, such as:\n",
    "   * Amount of signal light from the star `lsig`\n",
    "   * Amount of background light in the surrounding sky `lbkg`\n",
    "   * Temperature of star `temp`\n",
    "   * Transparency of sky `tsky`\n",
    "\n",
    "In order to determine the influence of these factors, and how much you need to correct for each of them, you consider 10.000 stars with known distances (measured by another method, e.g. triangulation). From these, you can find how well your own method works, make corrections to biases as needed, and finally find out how precise your calibrated method is. Happy calibration.\n",
    "\n",
    "\n",
    "## Your Task\n",
    "\n",
    "* As always look at the data and get a feel for each of the variables. A good idea might be to plot them all to know what range to expect them in.\n",
    "\n",
    "* First, consider the raw distribution of \"relative differences\" ($R_D$) between the observed and actual distance: $R_{D} = \\frac{(D_{obs} - D_{known})}{D_{known}}$. You'll notice that this distribution is far from an ideal (and narrow) gaussian, which is why you'll need to calibrate the underlying data. Your calibration technique should be able to do the following:\n",
    "\n",
    "    - Reduce the RMS of the distribution\n",
    "    - Correct the offset from zero\n",
    "    - Mitigate the tail at high values\n",
    "\n",
    "* Secondly, look at the distribution of the bias and relative precision as a function of the data variables. Try to identify variables that are correlated with the relative differences, and apply a correction factor to eliminate these biases.\n",
    "\n",
    "\n",
    "Note that if you are on average say 50% too high compared to the true values, then you need to correct by 50%, i.e. divide by (1 + 0.5), and in general, if your measurement is $f(x)$ off, where $f(x)$ describes the offset, then you need to divide by $(1 + f(T))$:\n",
    "\n",
    "$R_{D} = \\frac{d - d_{true}}{d_{true}}  \\rightarrow  d_{calib} = \\frac{d}{(1 + f(x))}$\n",
    "\n",
    "\n",
    "Thus, define `d_calib = d / (1 + f(x))`, and continue using `d_calib` when considering other effects.\n",
    "\n",
    "\n",
    "##  Author: \n",
    "- Troels Petersen ([email](mailto:petersen@nbi.dk))\n",
    "\n",
    "##  Date:   \n",
    "10th of December 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from iminuit import Minuit\n",
    "from scipy import stats\n",
    "import os, sys                                         # Modules to see files and folders in directories\n",
    "from os.path import dirname as parent_folder"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set som plotting standards:\n",
    "font = {'family' : 'serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 18}\n",
    "mpl.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "verbose = True\n",
    "Nverbose = 10\n",
    "SaveFigures = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Profile x of the two arrays x and y with defined number of bins and range\n",
    "# returns the x-values of the profile, the means and the standard deviations.\n",
    "def profile_x(x, y, bins=(50, 50), xyrange=[(0, 50), (-1,1)]):\n",
    "    \n",
    "    H, xedges, yedges = np.histogram2d(x, y, bins=bins, range=xyrange)\n",
    "    x_center = 0.5*(xedges[1:] + xedges[:-1])\n",
    "    y_center = 0.5*(yedges[1:] + yedges[:-1])\n",
    "    \n",
    "    wsums = H.sum(1)\n",
    "    \n",
    "    mask = wsums > 0\n",
    "    \n",
    "    mean = (H*y_center).sum(1)[mask] / wsums[mask]\n",
    "    mean_squared = (H*y_center**2).sum(1)[mask] / wsums[mask]\n",
    "    std = np.sqrt( mean_squared - mean**2 ) / np.sqrt(wsums[mask]) \n",
    "\n",
    "    return x_center[mask], mean, std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading the data (text) file:\n",
    "\n",
    "There are multiple ways we can retrieve the information from a text file, but here we use the `loadtxt` function of the `numpy` package. This is a very powerful function that allows you to parse an entire textfile in a single line, _provided that the text file is formatted properly_. For more advanced parsing option, you can also check out `numpy.genfromtxt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "filename = \"data_calib.txt\"\n",
    "dknown, dmeas, lsig, lbkg, temp, tsky = np.loadtxt(filename, unpack=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that this worked:\n",
    "if (verbose) :\n",
    "    for i in range(Nverbose) :\n",
    "        print(\"  Distance (known) = {:8.2f}    Distance (measured) = {:6.3f} \".format(dknown[i], dmeas[i]))\n",
    "print(\"The TOTAL number of entries read is: \", len(dmeas))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make Histograms and vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ioff() # We use to hide the empty plot on this cell\n",
    "fig_rel, ax_rel = plt.subplots(figsize=(15, 8));\n",
    "ax_rel.set_title('Hist relative resolution');\n",
    "Nbins = 200\n",
    "xmin, xmax = -2.0, 2.0\n",
    "binwidth = (xmax-xmin) / Nbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "fig_lsig2D, ax_lsig2D = plt.subplots(figsize=(12,8));\n",
    "ax_lsig2D.set_title('Hist lsig 2D')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop over data and make plots: (TO DO YOURSELF)\n",
    "----\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The initial relative resolution:\n",
    "distrel = (dmeas - dknown) / dknown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply calibration as found from the 2D plot below:\n",
    "dmeas_calib = dmeas      # So here is so far no calibration - that is what you have to change!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Once you have made the calibration, the (improved) relative resolution becomes:\n",
    "distrel_calib = (dmeas_calib - dknown) / dknown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Determine the resolution (quantified as the RMS) before and after the calibration:\n",
    "mask_distrel       = (xmin < distrel) & (distrel < xmax)\n",
    "mask_distrel_calib = (xmin < distrel_calib) & (distrel_calib < xmax)\n",
    "print(\"  The initial and final resolutions are:  {:6.3f}  and  {:6.3f} \\n\".format(\n",
    "        distrel[mask_distrel].std(ddof=1), distrel_calib[mask_distrel_calib].std(ddof=1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Compare the raw and calibrated data on a plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_rel.hist(distrel, bins=Nbins, range=(xmin, xmax), histtype='step', label='Raw')\n",
    "hist_rel_calib = ax_rel.hist(distrel_calib, bins=Nbins, range=(xmin, xmax), histtype='step', label='Calibration')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting attributes\n",
    "ax_rel.set_xlim(xmin, xmax)\n",
    "ax_rel.set_xlabel('Realitive precision (dmeas - dknown) / dknown',fontsize=18)\n",
    "ax_rel.set_ylabel('Frequency / 0.02', fontsize=18)\n",
    "ax_rel.legend(loc='best', fontsize=18)\n",
    "\n",
    "d = {'Entries':\"{:d}\".format(len(distrel[mask_distrel])), \n",
    "     'Mean':\"{:.3f}\".format(distrel[mask_distrel].mean()), \n",
    "     'STD Dev.':\"{:.3f}\".format(distrel[mask_distrel].std(ddof=1))}\n",
    "\n",
    "ax_rel.text(0.04, 0.95, nice_string_output(d), family='serif', fontsize=18, \n",
    "        transform=ax_rel.transAxes, verticalalignment='top')\n",
    "fig_rel.tight_layout()\n",
    "plt.ion()\n",
    "fig_rel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if (SaveFigures):\n",
    "    fig_rel.savefig('UncalibratedCalibrated.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "##  Calibration using 2D histograms:\n",
    "\n",
    "The way to investigate, if there is a (cor)relation between the target (relative distance) and any of the factors/variables/auxiliary observables, is to plot them against each other. If there are any visible relation, then this is an effect of the ourside factors (lsig, lbkg, temp, and tsky) on the wanted observable (relative distance), which we want to correct for.\n",
    "\n",
    "In the following, we produce this 2D histogram along with a 1D graph on top showing the mean and uncertainty for each bin in x. This graph shows any possible (average) relation in a way, which can be easily fitted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histograms for making calibration:\n",
    "ax_lsig2D.hist2d(lsig, distrel, bins=50, range=[(0, 50), (-1,1)], cmin=1, alpha=0.5)\n",
    "ax_lsig2D.set_xlabel('Signal light from star (lsig)');\n",
    "ax_lsig2D.set_ylabel('Realitive precision (dmeas - dknown) / dknown');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This calculates the 1D \"profile\" of the 2D histogram:\n",
    "x_center_lsig2D, mean_lsig2D, std_lsig2D = profile_x(lsig, distrel, bins=(50, 50), xyrange=[(0, 50), (-1,1)])\n",
    "x_binwidth_lsig2D = x_center_lsig2D[1] - x_center_lsig2D[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "ax_lsig2D.errorbar(x_center_lsig2D, mean_lsig2D, xerr=x_binwidth_lsig2D/2, yerr=std_lsig2D,  fmt='r.', \n",
    "                   ecolor='r', elinewidth=1, capsize=1, capthick=1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit of the effect from \"lsig\" in the 1D graph:\n",
    "def simple_fit(x, p0, p1):\n",
    "    return p0 + p1*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2_object_simple = Chi2Regression(simple_fit, x_center_lsig2D, mean_lsig2D, std_lsig2D)\n",
    "minuit_simple = Minuit(chi2_object_simple, pedantic=False, p0=0.0, p1=0.0)\n",
    "minuit_simple.migrad() # fit\n",
    "p0, p1 = minuit_simple.args\n",
    "print(\"Simple fit result:\")\n",
    "for name in minuit_simple.parameters:\n",
    "    print(\"Fit value: {0} = {1:.5f} +/- {2:.5f}\".format(name, minuit_simple.values[name], minuit_simple.errors[name]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_fit = np.linspace(0, 50, 1000)\n",
    "y_fit_simple = simple_fit(x_fit, *minuit_simple.args)\n",
    "ax_lsig2D.plot(x_fit, y_fit_simple, '--k')\n",
    "\n",
    "d = {'Entries':\"{:d}\".format(len(lsig)), \n",
    "     'Mean x' :\"{:.3f}\".format(lbkg.mean()), \n",
    "     'Mean y' :\"{:.5f}\".format(distrel.mean()), \n",
    "     'RMS x'  :\"{:.3f}\".format(lbkg.std(ddof=1)), \n",
    "     'RMS y'  :\"{:.5f}\".format(distrel.std(ddof=1))}\n",
    "\n",
    "ax_lsig2D.text(0.03, 0.2, nice_string_output(d), family='monospace', \n",
    "        transform=ax_lsig2D.transAxes, fontsize=16, verticalalignment='top')\n",
    "fig_lsig2D.tight_layout()\n",
    "fig_lsig2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your Task (reminder)\n",
    "\n",
    " As always look at the data and get a feel for each of the variables. A good idea might\n",
    " be to plot them all to know what range to expect them in.\n",
    "\n",
    " Next, consider the distribution of \"relative differences\" (RD) between the observed\n",
    " and actual distance: RD = (dist_obs - dist_known) / dist_known\n",
    " \n",
    " The Standard Deviation is 0.27, i.e. a 27% precision, and it is neither centered at zero\n",
    " (as it should, not to be biased), nor very Gaussian. There is also a long tail towards very\n",
    " high values. This is what you want to improve upon!\n",
    "\n",
    " Finally, there is the distribution of the bias and relative precision as a function of\n",
    " the signal luminosity (lsig). As you can see, the response does not depend on lsig,\n",
    " and so there seems to be no (varying) bias from this variable. Check the other three\n",
    " variables, and if you find some bias, try to correct for it, such that you get the\n",
    " relative difference to be the most narrow Gaussian centered around 0.\n",
    "\n",
    " Note that if you are on average say 50% too high compared to the true values, then\n",
    " you need to correct by 50%, i.e. divide by (1 + 0.5), and in general, if your\n",
    " measurement is f(x) off, where f(x) describes the offset, then you need to divide\n",
    " by (1 + f(T)). Why:\n",
    "\n",
    "   RD = (d - d_true) / d_true  =>  d_true = d / (1 + f(x))\n",
    "\n",
    " Thus define d_calib = d / (1 + f(x)), and continue using d_calib when considering\n",
    " other effects. In the end, your calibrated measurement can be the result of several\n",
    " consecutive corrections (note that the order might matter!).\n",
    "\n",
    "\n",
    "Questions:\n",
    "----------\n",
    "1. What corrections do you apply for each of the variables, and how much do each of them improve on the result?\n",
    "\n",
    "2. What is the final resolution you obtain?\n",
    "\n",
    "3. Do you think that further improvements are possible in this example (constructed by me)? How about in real and more complicated data?\n",
    "\n",
    "Advanced Questions:\n",
    "-------------------\n",
    "4. Try at the end to figure out, which variables the final resolution depends on.\n",
    "     Can you give an estimate of the uncertainty for each single star?"
   ]
  },
  {
   "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
}