ML-Kurs-SS2023/notebooks/02_fit_ex_3_sol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise 3: Least square fit with a 3rd order polynomial with iminuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data x,y and dy\n",
    "x = np.array([0.2 , 0.4 , 0.6 , 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2.,  2.2, 2.4, 2.6, 2.8 , 3.,  3.2 ,3.4, 3.6, 3.8,4.],dtype='d')\n",
    "dy = np.array([0.04,0.021,0.035,0.03,0.029,0.019,0.024,0.018,0.019,0.022,0.02,0.025,0.018,0.024,0.019,0.021,0.03,0.019,0.03,0.024 ], dtype='d')\n",
    "y = np.array([1.792,1.695,1.541,1.514,1.427,1.399,1.388,1.270,1.262,1.228,1.189,1.182,1.121,1.129,1.124,1.089,1.092,1.084,1.058,1.057 ], dtype='d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define fit functions - a 3rd order polynomial\n",
    "def pol3(a0, a1, a2, a3):\n",
    "    return a0 + x*a1 + a2*x**2 + a3*x**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# least-squares function = sum of data residuals squared\n",
    "def LSQ(a0, a1, a2, a3):\n",
    "    return np.sum((y - pol3(a0, a1, a2, a3)) ** 2 / dy ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import minuit2 fitting library\n",
    "from iminuit import Minuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  create instance of Minuit and use LSQ function to minimize\n",
    "LSQ.errordef = Minuit.LEAST_SQUARES\n",
    "m = Minuit(LSQ,a0=0.01, a1=0.05 ,a2=0.01 ,a3=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run migrad for minimization\n",
    "m.migrad()\n",
    "chi2 = m.fval / (len(y) - len(m.values))\n",
    "print (\"Chi2/ndof =\" , chi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run covariance  \n",
    "m.hesse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get correlation matrix\n",
    "cov = m.covariance\n",
    "print (cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# access elements of the numpy arrays\n",
    "print(cov[0, 1])\n",
    "print(cov[0, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run minos error analysis\n",
    "# The Minos algorithm uses the profile likelihood method to compute\n",
    "# (generally asymmetric) confidence intervals.\n",
    "m.minos()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a 2D contour of the function around the minimum for 2 parameters\n",
    "# and draw a 2 D contours up to 4 sigma of a1 and a2  \n",
    "m.draw_mncontour(\"a1\", \"a2\", cl=[1, 2, 3, 4])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.draw_profile(\"a2\",subtract_min=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# access fit results by parameter name and get minos asymetric errors\n",
    "print (m.merrors['a2'].lower)\n",
    "print (m.merrors['a2'].upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# more print out\n",
    "print (m.values,m.errors)\n",
    "print (m.errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access fit results\n",
    "a0_fit = m.values[\"a0\"]\n",
    "a1_fit = m.values[\"a1\"]\n",
    "a2_fit = m.values[\"a2\"]\n",
    "a3_fit = m.values[\"a3\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display fitted function \n",
    "x_plot = np.linspace( 0.1, 4.1 , 100 )\n",
    "y_fit = a0_fit + a1_fit * x_plot + a2_fit * x_plot**2 +  a3_fit * x_plot**3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.errorbar(x, y, dy , fmt=\"o\")\n",
    "plt.plot(x_plot,y_fit )           \n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n",
    "plt.title('iminuit exponential Fit')\n",
    "#plt.axis([0,30,-1.2,1.2])\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ML]",
   "language": "python",
   "name": "conda-env-ML-py"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}