From 87246a74ad5659771135bfc8329929403a94a42c Mon Sep 17 00:00:00 2001
From: Mauro Donega <mauro.donega@cern.ch>
Date: Tue, 30 Mar 2021 19:03:26 +0200
Subject: [PATCH] cleaned

---
 notebooks/ubinnedLikelihood.ipynb | 55 +++++++++++++++++++++++--------
 1 file changed, 41 insertions(+), 14 deletions(-)

diff --git a/notebooks/ubinnedLikelihood.ipynb b/notebooks/ubinnedLikelihood.ipynb
index 454051a..835fd69 100644
--- a/notebooks/ubinnedLikelihood.ipynb
+++ b/notebooks/ubinnedLikelihood.ipynb
@@ -4,7 +4,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Unbinned Maximum Likelihood fit"
+    "# Unbinned Maximum Likelihood fit\n",
+    "Implement by hand an unbinned maximum likelihood fit"
    ]
   },
   {
@@ -64,12 +65,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
+    "# this is my model for the data\n",
     "def likelihood_point(x, position, width):\n",
-    "    return 1.0/3/np.sqrt(2*np.pi*width**2)*np.exp(-0.5*((x-position)/(width))**2.0) + 2.0/3/np.sqrt(2*np.pi*5**2)*np.exp(-0.5*((x-91)/(5))**2.0)"
+    "    gauss1 = 1./np.sqrt(2*np.pi*width**2)*np.exp(-0.5*((x-position)/(width))**2.0)\n",
+    "    gauss2 = 1./np.sqrt(2*np.pi*5**2)*np.exp(-0.5*((x-91)/(5))**2.0)\n",
+    "    f = 1./3.\n",
+    "    return f * gauss1 + (1-f)* gauss2 \n",
+    "# 1.0/3/np.sqrt(2*np.pi*width**2)*np.exp(-0.5*((x-position)/(width))**2.0) + 2.0/3/np.sqrt(2*np.pi*5**2)*np.exp(-0.5*((x-91)/(5))**2.0)"
    ]
   },
   {
@@ -84,7 +90,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -104,16 +110,29 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 41,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "position: 117.72333147980623\n",
-      "negative error: [3.31211666]\n",
-      "positive error: [3.39091994]\n"
+      "Optimization terminated successfully. \n",
+      "\n",
+      "Fit results:\n",
+      "     fun: 291.4301328503986\n",
+      "     jac: array([-3.81469727e-06])\n",
+      " message: 'Optimization terminated successfully.'\n",
+      "    nfev: 44\n",
+      "     nit: 4\n",
+      "    njev: 22\n",
+      "  status: 0\n",
+      " success: True\n",
+      "       x: array([117.72327492])\n",
+      "\n",
+      "Uncertainty scan:\n",
+      "negative error: [3.3120601]\n",
+      "positive error: [3.3909765]\n"
      ]
     },
     {
@@ -131,13 +150,20 @@
    ],
    "source": [
     "solution = minimize(nll_1d, [100.0], method='CG')\n",
+    "print (solution.message,\"\\n\")\n",
+    "print (\"Fit results:\")\n",
+    "print (solution)\n",
+    "\n",
     "min_pos = solution.x[0]\n",
     "min0 = solution.fun\n",
+    "\n",
+    "# scan numberically the likelihood to get the uncertainties\n",
     "scan_points = np.linspace(110.0,126.0,50)\n",
     "plt.plot(scan_points, [nll_1d([x]) - min0 for x in scan_points])\n",
     "\n",
     "nll_1sigma = lambda x: nll_1d([x]) - min0 - 0.5\n",
-    "print(\"position:\", min_pos)\n",
+    "\n",
+    "print(\"\\nUncertainty scan:\")\n",
     "print(\"negative error:\", min_pos - fsolve(nll_1sigma, min_pos-0.5))\n",
     "print(\"positive error:\", fsolve(nll_1sigma, min_pos+0.5) - min_pos)\n",
     "plt.show()"
@@ -159,7 +185,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 42,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -179,20 +205,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 46,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "position: 118.31548192622421 width: 13.629783202046086\n"
+      "position = 118.3155 \n",
+      "width = 13.6297\n"
      ]
     }
    ],
    "source": [
     "solution = minimize(nll, [120.0, 10], method='CG')\n",
-    "print(\"position:\", solution.x[0], \"width:\", solution.x[1])"
+    "print(\"position = {:.4f}\".format(solution.x[0]), \"\\nwidth = {:.4f}\".format(solution.x[1]))"
    ]
   },
   {
@@ -204,7 +231,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 48,
    "metadata": {},
    "outputs": [
     {
-- 
GitLab