Week5/Calibration/Calibration_empty.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, but outlines the considerations.\n",
    "\n",
    "**Calibration play a role in essentially all physics measurements, in the eternal competition for precision.**\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 could depend 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 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 out how well your own method works, make corrections to biases as needed, and finally find out how precise your calibrated method is.\n",
    "\n",
    "Your calibration technique should be able to do the following:\n",
    "   * Correct the offset from zero\n",
    "   * Reduce the RMS of the (residual) distribution\n",
    "   * Mitigate the tail at high (residual) values\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",
    "* Then, 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 your method using the underlying data.\n",
    "\n",
    "* Go on and look at the distribution of the bias of the relative difference as a function of the data variables (i.e. 2D histograms). Try to identify variables that are correlated with the relative differences, extract (that is fit) for the size of the effect, and apply a correction factor to eliminate these biases.\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(x))$:\n",
    "\n",
    "$R_{D} = \\frac{d - d_{true}}{d_{true}}  \\rightarrow  d_{calib} = \\frac{d}{(1 + f(x))}$\n",
    "\n",
    "Thus, define `d_calib = d / (1 + f(x))`, and continue using `d_calib` when considering other effects.\n",
    "\n",
    "<p style=\"text-align: center;\"><b>Happy calibration!</b><p>\n",
    "\n",
    "\n",
    "##  Author:\n",
    "- Troels Petersen ([email](mailto:petersen@nbi.dk))\n",
    "\n",
    "##  Date:\n",
    "12th of December 2023"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Define profile function (2D histogram to 1D plot):\n",
    "\n",
    "Profile x of the two arrays x and y with defined number of bins and range returns the x-values of the profile, the means and the standard deviations.\n",
    "\n",
    "This enables you to fit a trend in a 2D histogram (which shows any variation of your key variable of interest with other variables) with a 1D function, and thus currect of the possible effect of this other variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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 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": [
    "## Look at data:\n",
    "\n",
    "As always, it is a good idea to take a quick look at the data. First the input variables, but also what the \"raw\" (i.e. before calibration) performance of the method is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We calculate the relative difference RD, between the estimated and the known distance:\n",
    "RD = (dmeas - dknown) / dknown\n",
    "Res_Raw = np.std(RD[(-2.0 < RD) & (RD < 2.0)])    # Resolution masking out outliers\n",
    "\n",
    "# ... and then plot it:\n",
    "xminRD, xmaxRD = -1.0, 2.0 \n",
    "NbinsRD = 60\n",
    "fig_RD, ax_RD = plt.subplots(figsize=(12,8))\n",
    "ax_RD.hist(RD, bins=NbinsRD, range=(xminRD,xmaxRD), histtype='step', label=f'Resolution - Raw: {Res_Raw:5.3f}')\n",
    "ax_RD.set_xlabel('Realitive difference RD = (dmeas - dknown) / dknown', fontsize=16)\n",
    "ax_RD.set_ylabel('Frequency / 0.05', fontsize=16)\n",
    "ax_RD.legend(fontsize=20);\n",
    "\n",
    "if (SaveFigures):\n",
    "    fig_rel.savefig('Resolution_Uncalibrated.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Consider (cor)relations between variable and method performance:\n",
    "\n",
    "The way to investigate, if there is a (cor)relation between the target (distance in terms of RD) and any of the factors/variables/auxiliary observables, is to plot them against each other in 2D histograms. If there are any visible relation, then this is an effect of the ourside factors (lsig, lbkg, temp, and tsky) on the wanted observable (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 (the x \"profile\"). This graph shows any possible (average) relation in a way, which can then be fitted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibrating for \"lsig\" based on 2D plot:\n",
    "\n",
    "Well, as it seems, $R_D$ does not depend heavily on lsig. However, in the following, I'm still going to test this, and use the result to correct for and thus calibrate $R_D$. Then you will have to do this for the other variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Comments on process:\n",
    "\n",
    "The first step of the calibration (with lsig) doesn't improve anything much (it was designed not to!), except moving the mean of $R_D$ to zero, such that the estimates of your new fancy method at least were not biased. The Standard Deviation remains around 0.26, i.e. a 26% precision, and the resolution is not very Gaussian either, which is not what you hoped for. **This is what you want to improve upon!**\n",
    "\n",
    "The way forward is to consider other effects (i.e. the other variables) and correct for possible effects from these. Note that once you have calibrated for one effect/variable, **you should continue with the calibrated version of the relative difference RD**, as there might be a correlation between the effects. In the end, your calibrated measurement can (and will) be the result of several consecutive corrections (thus the order might matter slightly, though this is a secondary effect).\n",
    " \n",
    "Also note, that even if the function you find for correcting isn't perfect, it is surely still better, than doing nothing, and so it should be applied anyhow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions:\n",
    "\n",
    "1. Which of the four variables (lsig, lbkg, temp, and tsky) causes the long tail towards high values in $R_D$?\n",
    "2. Try to rank the four variables according to importance (think about what that means).\n",
    "3. Do you think that the order of the calibration makes a difference? If so, should one start or end with the most important variables?\n",
    "4. Just how good a resolution can you get? Share it with the class..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points:\n",
    "\n",
    "From this exercise you should\n",
    "1. Generally, understand the concept of calibration<br>\n",
    "   (i.e. correcting your data for effects that you can measure/determine by some mean).\n",
    "2. Be able to apply a calibration to data<br>\n",
    "   (i.e. measure/determine the size of an effect, and correct data accordingly).\n",
    "3. Have an idea of what to think of **ahead** of an experiment to ensure good calibration<br>\n",
    "   (i.e. go through the list of needed measurements and think about how to ensure their precision)."
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}