|
|
{ "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()" ] }, { "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 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.8.16" } }, "nbformat": 4, "nbformat_minor": 4}
|
