{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7046e3dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iminuit version: 2.16.0\n"
     ]
    }
   ],
   "source": [
    "# stave.py\n",
    "# Computes Student's t average, with the number of degrees of\n",
    "# freedom nu related to the relative uncertainty r in the\n",
    "# standard deviations of the measurements, nu = 1/(2*r**2), see\n",
    "# G. Cowan, Eur. Phys. J. C (2019) 79 :133, arXiv:1809.05778.\n",
    "# Uses iminuit version 2.x (not compatible with v 1.x).\n",
    "# G. Cowan / RHUL Physics / December 2022\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import truncexpon\n",
    "from scipy.stats import truncnorm\n",
    "import iminuit\n",
    "from iminuit import Minuit\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "print(f\"iminuit version: {iminuit.__version__}\")      #  should be v 2.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aa82f19d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set input values\n",
    "mu = 10.                                  # initial value for fit\n",
    "y = np.array([17., 19., 15., 3.])         # measured values\n",
    "s = np.array([1.5, 1.5, 1.5, 1.5])        # estimates of std. dev\n",
    "v = s**2                                  # estimates of variances\n",
    "r = np.array([0.2, 0.2, 0.2, 0.2])        # relative errors on errors\n",
    "x = np.array([1., 2., 3., 4.])            # arbitrary x coord. for measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "150af0a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function to be minimized\n",
    "class NegLogL:\n",
    "\n",
    "    def __init__(self, y, s, r):\n",
    "        self.setData(y, s, r)\n",
    "        \n",
    "    def setData(self, y, s, r):\n",
    "        self.data = y, s, r\n",
    "\n",
    "    def __call__(self, mu):\n",
    "        y, s, r = self.data\n",
    "        v = s ** 2\n",
    "        lnf = -0.5*(1. + 1./(2.*r**2))*np.log(1. + 2.*(r*(y-mu))**2/v)\n",
    "        return -np.sum(lnf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "add98882",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize Minuit and set up fit:\n",
    "negLogL = NegLogL(y, s, r)          # instantiate function to be minimized\n",
    "parin     = np.array([mu])          # initial values\n",
    "parname = ['mu']\n",
    "parstep = np.array([0.5])           # initial setp sizes\n",
    "parfix  = [False]                   # change these to fix/free parameters\n",
    "parlim  = [(None, None)]\n",
    "m = Minuit(negLogL, parin, name=parname)\n",
    "m.errors = parstep\n",
    "m.fixed = parfix\n",
    "m.limits = parlim\n",
    "m.errordef = 0.5                           # errors from lnL = lnLmax - 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7f629d63",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "par index, name, estimate, standard deviation:\n",
      "   0 mu          =  16.205276  +/-  0.987898\n"
     ]
    }
   ],
   "source": [
    "# do the fit, extract results\n",
    "m.migrad()                                        # minimize -logL\n",
    "MLE = m.values                                    # max-likelihood estimates\n",
    "sigmaMLE = m.errors                               # standard deviations\n",
    "cov = m.covariance                                # covariance matrix\n",
    "rho = m.covariance.correlation()                  # correlation coeffs.\n",
    "muHat = MLE[0]\n",
    "sigma_muHat = sigmaMLE[0]\n",
    "\n",
    "print(r\"par index, name, estimate, standard deviation:\")\n",
    "for i in range(m.npar):\n",
    "    if not m.fixed[i]:\n",
    "        print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n",
    "         \"{:.6f}\".format(MLE[i]), \" +/- \", \"{:.6f}\".format(sigmaMLE[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "07d791e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwqklEQVR4nO3dd3RVVfr/8feTRhIggZBQEwgtoUroxYYURxB0RhQYRcWG3cEy8xsYG+qgM5ZBZ8aCBawICiogNlBE6S10QgmEhF6TQHqyf3/cy3dFJj252bc8r7XO4rbkfE4uebLv3vvsI8YYlFJKeR8/2wGUUkq5hhZ4pZTyUlrglVLKS2mBV0opL6UFXimlvJQWeKWU8lJa4JVSyktpgVfWiMhoEVkhIlkisrSE5/1F5DkROSQimSKyUUQalPK9XhKR3c7X7RSRW6qZLVZEfnJm2ykiQ8p4bQMReV9Ejjm3py94PkFEfhGRdBFJE5EnL3j+QRHZJyIZIrJORC6pTnZfOB5VQcYY3XSzsgFDgNHAk8DSEp5/DvgRaAUI0AUILuV7TQE64Gi09AVOAwOqkW0l8AoQAowCzgBRpbx2BvAZEArEAnuB24o9vx34O+APtAUOA9c4n+sLnAN6Oo/xXuA44F+BjE8DT3vL8ehW85v1ALrV0hsN+4E/A5udv4DvAk2Ab4BMYDHQ0PnagUBaCV8/xEXZ7rywwAMNgbNA2yp+z/nAo1X82jggF6hf7LFfgHtKef0JoHex+5OBX4rdzwI6Fbv/GTDJeXsMsKbYc3UBAzSrQM4KFXhPOR7dan7TLhrfMgoYiuMXfiSO4j4ZiMTR8n2oKt9URP4qImdK26qYtStQAFwvIkdEZJeI3F/BPCFAb2BbFffdGUg2xmQWe2yT8/FSd3vB7S7F7k8DbhGRQBGJB/rj+IMKjvfAX0T6iog/cDuQCBypYvaSeNvxqAoKsB1A1ap/G2OOAojIL8AxY8xG5/0vgMFV+abGmBeAF2ospUM0EI7jj1FroD2wRER2GWN+KOdr38RRwL6r4r7rAekXPJYOtCjl9d8CfxWRW3F8KrodR/fGeQuBD4DHcHRrPGOMWet8LhOYC/yKo5CeAYYZZ/O3hnjb8agK0ha8bzla7HZ2Cffr1W6cMmU7/33GGJNtjNkMfAoML+uLRORFHK3N0aUVFRH5RkTOOrebSnjJWSDsgsfCcBSvkjzkzLsb+AqYBaQ59xWBo2A+AwQDMcDvROQ+59feiaOAdgaCgHHAQhFpXkr2hcU+Gf0VRyE+/2lpYSn53PZ4lGtpgVclOUexFpvzo3ZUaS8WkcnFCub/bFXMsNn5b4VbfiIyBRgGXGmMySjtdcaYYcaYes7t4xJesg1oIyL1iz3WjVK6fIwxp4wxNxljmhpjOuP4vVrjfLoNUGiM+cAYU2CMSeO3f6i6AQuMMbuMMUXGmG9xDFoOKGVfI4wxDYwxDXB8anrh/H1jzIhSDtltj0e5lhZ4VZJdQLCIXC0igcDjQJ3SXmyMmVqsYP7PVtrXOadBBuPoKvQTkWDn/jDG7MUxEPg3EakjIh1xDOCV2EoVkUnAjcBQY8zJKh73+ePZhaPf+Clnpj8AF+Hoeihp321FpJHzeIYBE3DMAALHz1JE5EYR8RORps7j2OR8fi1wtYi0EYfzYyRbq3MM3nw8qhJsj/LqVjsbF8yCAT6i2AwMHB+tFxe7Px5Hy+sYjr7W33x9DWUaj6OFXnybWez5Fji6A84CycDdxZ67CdhW7L7BMVPkbLFtcjWyxQJLcXRVJF3ws7sUOFvs/mjgEI7ZJYnA7y74XoNwFL50HIONbwOhzucER3fHARxdJjuAmyuY8WkqPk3S7Y9Ht5rfxPmmKKWU8jLaRaOUUl5KC7xSSnkpLfBKKeWltMArpZSXcqszWSMjI01sbKztGEr9RlJSEgDx8fGWkyj1v9avX3/CGFPieSpuVeBjY2NZt26d7RhK/cbAgQMBWLp0qdUcSpVERFJKe067aJRSyku5VQteKXf0+OOP246gVJVogVeqHEOGlHrxI6XcmnbRKFWOxMREEhMTbcdQqtK0Ba9UOSZOnAjoIKvyPNqCV0opL6UFXimlvJQWeKWUsuinnceYsXwfeQVFNf69tcArpZRFry/dw/sr9hPgJ+W/uJJ0kFWpckydOtV2BOWldh3NZO3+00wa1gE/LfBK1b4BA/Ryoso1Pll9gCB/P67vGe2S769dNEqVY8WKFaxYscJ2DOVlsvIKmLshjWFdm9KoXqmXPK4WbcErVY7JkycDOg9e1ayFmw6TmVPATX1buWwf2oJXSikLPl5zgPaN69E7tqHL9qEFXimlatnWg+lsSj3DjX1bIlLzg6vnaYFXSqla9vHqAwQH+nFdd9cMrp6nBV4ppWpRZk4+XyUeZORFzQkPDXTpvnSQValyTJs2zXYE5UW+SjxEVl4hN/Vz3eDqeVrglSpHQkKC7QjKSxhj+Hj1ATo3D6NbdLjL96ddNEqVY/HixSxevNh2DOUFNhw4w47DGS4fXD1PW/BKleO5554D9MpOqvo+XLmf+nUC+H1Ci1rZn7bglVKqFpw4m8uiLUcY1TOaunVqp22tBV4ppWrB7LWp5BUWMa4WBlfP0wKvlFIuVlhk+HhVChe3a0S7xvVqbb9a4JVSysWW7DjKofQcbq7F1jvoIKtS5XrrrbdsR1Ae7sNVKTQLD2ZIxya1ul8t8EqVIz4+3nYE5cGSj5/ll90neHRoHAH+tdtpol00SpVjwYIFLFiwwHYM5aE+XJVCoL8wtk/LWt+3tuCVKsfLL78MwMiRIy0nUZ4mK6+Az9enMaxLM6Lqu+aiHmXRFrxSSrnIlxsPkZlTwC39a3dw9TyXF3gR8ReRjSKy0NX7Ukopd2GM4YOV++nYLIyerVx3UY+y1EYL/k/AjlrYj1JKuY11KafZeSSTW/q3qpV1Z0ri0gIvItHA1cA7rtyPUkq5mw9WplA/OIBrE5pby+DqQdZpwF+A+qW9QEQmABMAWras/VFmpcrz4Ycf2o6gPMyxzBy+3XqYm/vFEhpkby6Ly1rwIjICOGaMWV/W64wx040xvYwxvaKiolwVR6kqi4mJISYmxnYM5UE+WX2A/ELDuH52G62u7KK5GLhGRPYDnwKDROQjF+5PKZeYPXs2s2fPth1DeYjcgkI+WnWAK+KjaBNVe+vOlMRlBd4YM8kYE22MiQXGAj8aY8a5an9Kucobb7zBG2+8YTuG8hBfbz7MibO5jL+4te0oOg9eKaVqijGGGcv30zaqLpe1j7Qdp3YKvDFmqTFmRG3sSymlbNlw4DRbDqYz/uLW1qZGFqcteKWUqiHvLd9PWHAAo3rUziX5yqMFXimlasChM9l8u/UIY/u0tDo1sjj3SKGUG/v8889tR1Ae4MNVKRhjrK07UxIt8EqVIzLS/mCZcm/ZeYXMWnOAKzs1JbphqO04/8crumiOZuRw+lye7RjKS82cOZOZM2fajqHc2JeJBzmTlc9tF8fajvIbHl/g07PzGfLyz0xbvMt2FOWltMCrshhjmLncsWpkn9YRtuP8hscX+PCQQEYmNOfj1QdIPn7WdhyllI9ZufckSUczue3iWLeYGlmcxxd4gIlD2hMU4Mc/v02yHUUp5WPeW76fRnWDuKabvVUjS+MVBb5x/WDuvqwt3247wrr9p2zHUUr5iJST51iy8yg39m1JcKC/7Tj/wysKPMBdl7Wmcf06TF20A2OM7ThKKR8wY/l+/EUY1899pkYW5zUFPjQogEeGxrHhwBm+3XrEdhzlRRYtWsSiRYtsx1BuJj0rnznrUrkmoTlNwoJtxymR1xR4gBt6xRDXpB7/+HYneQVFtuMoLxEaGkpoqPvMbVbu4ZM1B8jKK+TOS9rYjlIqryrw/n7CpGEd2X8yi09Wp9iOo7zE66+/zuuvv247hnIjeQVFzFyxj0vaRdKpeZjtOKXyqgIPMDA+igFtG/Hqkt1k5OTbjqO8wJw5c5gzZ47tGMqNLNh0iKMZudx1mfu23sELC7yIMHl4R05n5fPG0r224yilvIwxhrd/SSa+SX23WPO9LF5X4AG6tAjnD91b8N6v+zh0Jtt2HKWUF1m+5yQ7j2Ryx6XuseZ7WbyywAM8emUcBnjpez35SSlVc6b/kkxU/Tpcm+B+JzZdyGsLfHTDUG4bEMsXGw+y7VC67ThKKS+QdCSTZbuOM35ALHUC3O/Epgt5bYEHuO+KdoSHBPL8op168pOqsqVLl7J06VLbMZQbeOeXZEIC/bmpb0vbUSrEqwt8eEggDw5qz697TvDzruO24yilPNixjBy+SjzE6F7RNAgNsh2nQry6wAPc3K8VLSNCeX7RTgqLtBWvKu+ll17ipZdesh1DWfb+yv3kFxVx+yWtbUepMK8v8EEBfvzlqniSjmYyd32a7TjKAy1cuJCFCxfajqEsysor4KNVB/hdp6a0alTXdpwK8/oCD3B112YkxDTgpe+TOJdbYDuOUsrDfL4+jfTsfO66zHNa7+AjBV5EeGJEJ45l5vLmz3ryk1Kq4goKi5i+LJmerRrSs5V7XbGpPD5R4AF6tmrINd2aM31ZMgf15CelVAUt3HyYtNPZ3Ht5W9tRKs1nCjzA/xvWAYB/fLPTchLlSUJCQggJCbEdQ1lgjOGNpXuJb1KfQR0a245TaT5V4Fs0CGHCZW2Yv+kQ61NO246jPMQ333zDN998YzuGsuCnpGMkHc3knoFt8PNz72UJSuJTBR7gnsvb0rh+HZ5duJ0inTaplCrDG0v30qJBCCMucv9lCUricwW+bp0A/vy7eBJTz7Bg8yHbcZQHePbZZ3n22Wdtx1C1bO3+U6zdf5oJl7Uh0N8zS6Vnpq6mUT2i6dIijBe+2Ul2XqHtOMrNLVmyhCVLltiOoWrZG0v3ElE3iNG9YmxHqTKfLPB+fsKTIzpzOD2H6cuSbcdRSrmZHYcz+HHnMW4bEEtIkPsvKlYanyzwAH1aRzC8a1Pe/HkvR9JzbMdRSrmRt37eS90gf27pH2s7SrX4bIEHmDSsI4VFhn9+p9MmlVIOqaeyWLD5MDf2bUl4aKDtONXi0wU+JiKU2y9pzbwNB9mcdsZ2HOWmGjVqRKNGjWzHULXk7V+S8RO44xL3vt5qRbiswItIsIisEZFNIrJNRKa4al/Vcf8VbYmsF8QzC7brmvGqRHPnzmXu3Lm2Y6hacOJsLrPXpjKqRzRNw4Ntx6k2V7bgc4FBxphuQAJwlYj0c+H+qqR+cCCPXRnPupTTLNx82HYcpZRFM5bvI6+wiAmXeX7rHVxY4I3DWefdQOfmlk3kG3rF0Ll5GFMX7SArT1ebVL81adIkJk2aZDuGcrH0rHzeX5HC8C7NaBNVz3acGuHSPngR8ReRROAY8IMxZnUJr5kgIutEZN3x43auuuTvJ0y5xjFt8o2lutqk+q2VK1eycuVK2zGUi81YsY+zuQU8MKid7Sg1xqUF3hhTaIxJAKKBPiLSpYTXTDfG9DLG9IqKinJlnDL1io3g9wnNeWtZMgdOZlnLoZSqfZk5+bz36z6u7NSEjs3CbMepMbUyi8YYcwZYClxVG/urqknDOxLgJzz79XbbUZRSteiDlSlk5BTw0OD2tqPUKFfOookSkQbO2yHAEMCtJ5w3CQvmwUHt+WH7Ub1It1I+4mxuAW//ksygDo3p0iLcdpwa5coWfDPgJxHZDKzF0Qfv9he2vP2SWFpH1mXKgm3kFRTZjqPcQHR0NNHR0bZjKBf5aFUKZ7LyedCL+t7PC3DVNzbGbAa6u+r7u0qdAH+eHNGJ22auZeaKfUy4zPOu4qJq1kcffWQ7gnKR7LxC3l6WzKXtI+nesqHtODXOp89kLc0VHRozuENjXl28m2MZuk6NUt7q49UpnDyXx5+8rO/9PC3wpXhiRCfyCw0vfOvWwwaqFkycOJGJEyfajqFqWE5+IW8tS2ZA20b0ivWsi2lXlBb4UsRG1uXOSx3r1KxPOWU7jrIoMTGRxMRE2zFUDZu9NpXjmbk8OMg7W++gBb5M91/RjqZhwTw9fzuFenk/pbxGbkEhbyzdS5/YCPq18c7WO2iBL1PdOgFMGt6BLQfTmbMu1XYcpVQN+Xx9GkcycnhwcDtEPO9i2hWlBb4c13RrTp/YCF78LokzWXm24yilqim3oJDXf9pL95YNuKRdpO04LqUFvhwiwpRrO5Oenc8/v0uyHUdZEBcXR1xcnO0YqobMXpvKwTPZPDI0zqtb7+DCefDepGOzMG4bEMu7y/dxQ89or5wvq0o3ffp02xFUDcnJL+Q/P+6hT+sIr2+9g7bgK2zi0Dia1A/m8S+3UlCoZ7gq5Yk+WpXCscxcHvWB1jtoga+wenUCeHJkJ7YdyuCjVSm246haNGHCBCZMmGA7hqqmc7kFvL50L5e2j6RvG9+4BKMW+EoY1qUpl8VF8fL3u/QMVx+ya9cudu3aZTuGqqaZK/Zz6lwejwz1nfEULfCVICI8c01ncguLeO7rHbbjKKUqKD07n7d+3suQjo19agxNC3wlxUbW5d7L2zJ/0yF+3X3CdhylVAW8++s+MnIKeNiHWu+gBb5K7h3YllaNQnnyq63kFhTajqOUKsPpc3m89+s+hndtSufm3rXee3m0wFdBcKA/U67pTPKJc7y9LNl2HOViCQkJJCQk2I6hquitZcmcyyvg4SG+1XoHnQdfZQPjGzO8a1P+/eMerunWgpaNQm1HUi4ybdo02xFUFR3PzOX9Ffu5tltz2jepbztOrdMWfDU8OaIzAX7CU/O3YowuRqaUu3lj6V7yCov4kw+23kELfLU0DQ/m4aFx/JR0nO+2HbUdR7nIuHHjGDdunO0YqpLSTmfx0aoURvVoQevIurbjWFFmF42ItCzreWPMgZqN43luHRDL5+vTeHr+Ni5u14j6wYG2I6kalpaWZjuCqoJXftiFCEz00dY7lN+C/xpY6Pz362L3VwP7XBvNMwT6+/H8dV05mpnDS7oYmVJuYfuhDL7YeJDxF8fSvEGI7TjWlFngjTFdjTEXOf/tCowElgNngYm1kM8jdG/ZkFv7x/LBqhTWp5y2HUcpn/ePb3cSFhzIfZe3sx3Fqgr1wYtIexGZCXwDrAc6GWP+7cpgnuax38XTNCyYyfO2kFegi5EpZcuKPSf4eddx7r+iLeGhvt1lWmaBF5EuIjILmAssBroYY94xxuTXSjoPUq9OAM9e24Wko5lMX7bXdhxVg/r370///v1tx1AVUFRkeOHbnbRoEMIt/WNtx7GuvHnwm4BUHH3vfYA+xZfYNMY85LponmdIpyYM79qU137cw/CuzWgTVc92JFUDnn/+edsRVAUt2nqYzWnpvHxDN4ID/W3Hsa68An8HoBO8K+HpkZ35ZfcJJn+xhVl39fOJNaeVcgd5BUW8+F0SHZrW5/fdW9iO4xbKK/CfAvWNMceLPygijYEMl6XyYI3Dgpk0rCOTv9jCZ+vSGN07xnYkVU2jRo0CYO7cuZaTqLLMWnOAlJNZzLitN/5+2rCC8gdZXwMuLeHxocC/aj6OdxjbO4Y+sRH8fdEOjmfm2o6jqunkyZOcPHnSdgxVhrO5Bby2ZDf92kQwMC7Kdhy3UV6Bv8QYM+/CB40xHwOXuSaS5/PzE6Ze14XsvEKmLNhmO45SXm/6smROnstj0rCO2i1aTHkFvqyflC5zUIZ2jevzwKB2LNx8mO+2HbEdRymvdSwjh3d+Sebqrs3oFtPAdhy3Ul6RPiYifS58UER6A8dLeL0q5t6BbenYLIzHv9zKmaw823GU8kovfpdEQaHhL1fF247idsor8H8G5ojI0yIy0rlNAeY4n1NlCPT348XrL+LUuTyeXaiX+PNUgwcPZvDgwbZjqBJsPZjO5xvSGH9xLK0a+eaCYmUpcxaNMWaNswV/PzDe+fA2oK8x5piLs3mFLi3Cuffytvznpz2M6NaMK+Ib246kKumJJ56wHUGVwBjDMwu30zA0iAcG+faSBKUp70zWlsaYY8aYp4wxo5zbk1rcK+fBwe1o37gek+dtISNHTwJWqiZ8t+0Ia/ad4pGhcYTpKq4lKq+L5svzN0REJwFXUZ0Af168oRtHM3J4fpF21XiaYcOGMWzYMNsxVDG5BYVMXbSTuCb1GKvnmpSqMrNo2lTmG4tIjIj8JCI7RGSbiPyp8vG8R0JMA+66tA2z1qTy6+4TtuOoSsjOziY7O9t2DFXM+yv2c+BUFo9f3YkAf53QV5ryfjKmlNsVUQA8aozpCPQD7heRTpX8Hl7l4aFxtImsy1/nbeZcboHtOEp5pJNnc/n3kj1cER/FZXpSU5nKK/DdRCRDRDKBi5y3M0QkU0TKXKrAGHPYGLPBeTsT2AH49AIRwYH+/PP6izh4Jpt/fLvTdhylPNLLP+wiK7+Qv13d0XYUt1feBT/8jTFhxpj6xpgA5+3z98MquhMRiQW647gS1IXPTRCRdSKy7vhx759a3ys2gvEDYvlgZQrL92hXjVKVsfVgOrPWHOCW/q1o17i+7Thuz+WdVyJSD8d68hONMf/T6jfGTDfG9DLG9IqK8o2PW3/5XQfaRNblz59t0lk1HmDEiBGMGDHCdgyfV1RkePKrrTSqG8TDQ333OquV4dICLyKBOIr7xyWtaeOrQoL8eWVMAkczc5kyf7vtOKocjz32GI899pjtGD5v3saDbDhwhv93VQedFllBLivw4ljx511ghzHmFVftx1MlxDTg/oFtmbshTdeqUaocGTn5vPDNTrq3bMCoHtG243gMV7bgLwZuBgaJSKJzG+7C/XmcBwa1p3PzMCbP28KJs7qssLsaOHAgAwcOtB3Dp726eDcnz+XyzDVd8NO13ivMZQXeGPOrMUaMMRcZYxKc2yJX7c8TBQX48a8xCWTmFjB53haM0YtnKXWhXUczmbliP3/s05Ku0eG243gUPUPAsrgm9Xnsyji+336UeRsO2o6jlFsxxvDUV9uoVyeAP1+pq0VWlhZ4N3DHJW3oExvB0/O3cfCMnjGp1HnzNx1iZfJJHrsyjoZ1g2zH8Tha4N2Av5/w0g3dKDKGP3+2iaIi7apRKj07n2cX7uCi6HBu7NvKdhyPpAXeTbRsFMoTIzqxYu9J3vk12XYcVczo0aMZPXq07Rg+58XvdnLqXC5T/9BVL6JdRWWuB69q15jeMfyUdIwXv0tiQNtIurTQASV3cN9999mO4HM2HjjNx6sPcNuA1vp7UA3agncjIsIL111Eo7p1eGjWRrLydEEyd5CVlUVWVpbtGD6joLCIyV9spWlYMI9cqWesVocWeDfTsG4Qr4zpxr6T53hmgZ7l6g6GDx/O8OF6CkdtmbF8PzsOZ/DUyM7Uq6OdDNWhBd4NDWgbyT2Xt+XTtal8s+Ww7ThK1ZqDZ7J55YddDO7QmN91bmI7jsfTAu+mHhkaR7focP46bwuHdOqk8gGOOe9bAZhybWccq52o6tAC76YC/f14dWx38guLeHh2IoU6dVJ5uYWbD7N4xzEeGRpHdMNQ23G8ghZ4NxYbWZcp13Rm9b5TvPnzXttxlHKZU+fyeHr+NrrFNOD2S1rbjuM1dATDzV3fM5qfdx3nlR920a9NBD1bRdiO5HPGjx9vO4LXm7JgGxk5+fxz1EU6570GaQvezYkIf/9DV1o0COHBTzZy+lye7Ug+Z/z48VrkXWjJjqN8lXiI+69oR3xTvUpTTdIC7wHCQwL57409OHE2j0fmJOpSBrXsxIkTnDihl1d0hYycfP72xVbim9TnvoHtbMfxOlrgPUTX6HAeH9GRn5KO89YyXcqgNl1//fVcf/31tmN4pecX7eRYZg7/vP4iggK0HNU0/Yl6kJv7teLqrs146fsk1u4/ZTuOUtWyYs8JZq05wJ2XtqFbTAPbcbySFngPIiK8MKorMQ1DeOCTDZzUq0ApD5WRk8+fP99Mm8i6PDxElyNwFS3wHqZ+cCD/vakHp7PyeXiOLi2sPNOzC7ZzOD2bl0d3IyTI33Ycr6UF3gN1bh7OUyM7sWzXcd7Q+fHKw/yw/SifrU/jvoHt6N6yoe04Xk3nwXuoG/u0ZHXyKV7+PonuMQ0Y0C7SdiSvde+999qO4DVOns1l0rzNdGoWxkOD29uO4/W0wHsoEWHqdV3ZcTiDB2ZtZP4DF+vp3S4yZswY2xG8gjGGv32xlYzsAj6+M0FnzdQC/Ql7sHp1Anjr5p7kFxRxz0fryckvtB3JK6WmppKammo7hsf7MvEg3247wiNXxukJTbVEC7yHaxNVj2ljE9h6MIPJX2zBGB10rWk333wzN998s+0YHi3tdBZPfrWNXq0actelbWzH8Rla4L3A4I5NmDikPfM2HOSDlSm24yj1GwWFRfzp00Qw8K8xCbrWTC3SAu8lHhrUniEdG/Pswu2s2acnQSn38dqPe1ifcprn/tCFmAgdJ6pNWuC9hJ+f8MqYBFpGhHLfx+s5nK4XCVH2rU4+yX9+3M2oHtFcm9DCdhyfowXei4QFB/LWzT3Jzivkno82kFugg67KnjNZeUycnUirRnWZcm1n23F8khZ4L9O+SX1eHt2NTalnmDRPB11rwqOPPsqjjz5qO4ZHMcbw17lbOHE2l9fGdteLZ1uiP3UvdFWXZjw8JI5/Ld5F26h63H+FLsNaHSNHjrQdweN8suYA3247wt+Gd6RrdLjtOD5LC7yXemhwO5JPnOXF75JoHVmX4V2b2Y7ksZKSkgCIj4+3nMQzbElLZ8r87VweF8Udevk9q7TAeykR4R+jLiL1VBaPzEmkRYMQXZK1iu6++24Ali5dajeIB0jPyufej9cTWS+IaWMS8NMpkVZpH7wXCw70Z/otvYisV4c7P1jHoTM6s0a5TlGR4ZE5iRzNyOG/N/WgYd0g25F8nhZ4LxdZrw7v3tqb7LxC7nx/HedyC2xHUl7qjZ/3smTnMR6/upOuEukmtMD7gPim9fnPjd3ZeSSDP32aSKGuIa9q2Iq9J3j5+yRGdmvOLf1b2Y6jnFxW4EXkPRE5JiJbXbUPVXED4xvz1MjOLN5xlOe+3q7TJ1WNOZyezUOzNtI6si7PX9cVEe13dxeuHGSdCfwH+MCF+1CVcOuAWPafPMeM5ftpFh7MhMva2o7kER5//HHbEdxWTn4hd3+4nuy8Qmbd1U/nu7sZl70bxphlIhLrqu+vquaJqztxLDOXqYt2ElW/Dn/oHm07ktsbMmSI7QhuyXEy02Y2p6Uz/eaetG+iSwC7G/1z62P8/IRXRnfj5Nlc/vzZZiLr1eHS9lG2Y7m1xMREABISEqzmcDfTlyXzZeIhHh0ax5Wdm9qOo0pgfZBVRCaIyDoRWXf8+HHbcXxCnQDH9Ml2jetxz4fr2Xow3XYktzZx4kQmTpxoO4Zb+SnpGC98u5OruzbjgUF6prS7sl7gjTHTjTG9jDG9oqK0JVlbwoIDmXlbHxqEBjF+xlr2nThnO5LyEHuPn+WhWRvp2DSMF2+4SAdV3Zj1Aq/saRoezPu396HIGMa9s1pPhFLlOnUujztmriXI34/pt/QkNEh7ed2ZK6dJzgJWAvEikiYid7hqX6rq2jWuxwe39yEjO59x76zmeGau7UjKTeXkF3Ln+2s5nJ7D9Ft66kXePYDLCrwx5o/GmGbGmEBjTLQx5l1X7UtVT5cW4cy4rTeH03O4+d3VpGfl246k3ExRkeHh2YlsTD3DtDEJ9GwVYTuSqgDtolEA9IqNYPotPUk+fo5bZ6zhrC5p8H+mTp3K1KlTbcewauqiHXyz1bH87zBdmdRjaIFX/+fS9lG89sfubDmYzh0z15KVp0UeYMCAAQwYMMB2DGtmLt/HO7/uY/yAWF3+18NogVe/cVWXprwyuhtr95/ithla5AFWrFjBihUrbMew4psth5mycDtDOzXhiRGddMaMh9EhcPU/zl8c+eHZidw2Yy0zbuvt07MlJk+eDPjeevA/7zrOQ59upHtMA14dm4C/ru3ucbQFr0p0bUIL/jUmQVvyPmrd/lPc/eE62jWuz4zb+vj0H3hPpgVelerCIq9ryfuGbYfSuW3mWpqHh/DhHX0IDwm0HUlVkRZ4VabzRX5dymnGvbuaM1l5tiMpF0o+fpZb3l1D/ToBfHhnXyLr1bEdSVWDFnhVrmsTWvD6TT3YdjCDsdNXcSwzx3Yk5QL7T5zjpndWIwIf3dmXFg1CbEdS1aQFXlXI7zo35b3xvTlwKosb3lxJ6qks25FqzbRp05g2bZrtGC6VfPwsY6avJLegiA/v6EubqHq2I6kaoAVeVdgl7SP58I6+nD6Xxw1vrmTPsUzbkWpFQkKCVy8VvOfYWcZOX0VBoWHWXf3o2CzMdiRVQ7TAq0rp2aohs+/uT0GR4fo3V7Ju/ynbkVxu8eLFLF682HYMl9h9NJOx01dRZODTCf2Ib6oX7fAmWuBVpXVsFsbce/vTMDSIG99ZzdebD9uO5FLPPfcczz33nO0YNW7nEceYip84irtekcn7aIFXVdKqUV3m3juAri3Cuf+TDby9LFkv5O1B1uw7xQ1vriTQ349PJ/SjXWPtc/dGWuBVlUXUDeLjO/tydddm/H3RDp6av43CIi3y7u67bUcY9+5qourX4fN7++uAqhfT09NUtQQH+vPvP3anRcMQpi9LJuVkFq+N7U54qJ4c444+XXOAyV9s4aLoBrw3vjcRdYNsR1IupC14VW1+fsLk4R15/rqurNh7gmv/+yu7j/rGDBtPYYzh30t289d5W7gsLopP7uqrxd0HiDv1m/bq1cusW7fOdgxVDev2n+KejzaQnVfAv8YkcGXnprYjVVtSUhIA8fHxlpNUTU5+IZPnbWHexoNc170F/7j+IgL9tW3nLURkvTGmV0nP6busalSv2AgWPHgxbRvXY8KH65m2eBdFHt4vHx8f77HF/VhmDn98exXzNh7k0aFxvDy6mxZ3H6LvtKpxzcJDmHN3f67r0YJpi3dz64w1Hn2t1wULFrBgwQLbMSpt68F0rv3PcnYezuTNcT14cHB7Xc/dx2gXjXIZYwyz16by1PxthIUEMm1MAhe3i7Qdq9IGDhwIeNZ68PM3HeIvn28iIjSIt2/tRefm4bYjKRfRLhplhYgwtk9LvnrgYsJDAhn37mpe+WGXTqV0oZz8QiZ/sYWHZm2kS/NwvnrgEi3uPkwLvHK5Dk3DmP/AxYzqEc1rS3Yz+q2V7DtxznYsr7P3+Fl+/9/lfLL6AHdf3oZZE/oRVV+X+/VlWuBVrQgNCuClG7oxbUwCu49mMuzVZcxYvs/jB2DdxZcbDzLy379yNCOHGeN7M2lYRx1MVXqik6pdv+/egn5tGvHXeZuZsmA73207wovXdyMmItR2NI906lweT361lYWbD9M7tiGv/bE7zcJ1HXfloIOsygpjDHPWpfLswh0UGcMjQ+O4dUCsW7Y6U1NTAYiJibGc5Le+2XKYx7/cSkZOPg8Oas99A9sS4IY/P+VaZQ2yaoFXVqWdzuKJL7fyU9Jx4pvU59nfd6FP6wjbsdzaybO5PDl/G19vPkyXFmG8dEM3OjTVNdx9lRZ45daMMXy//SjPLNjOwTPZXNejBZOGdXSbAcLZs2cDMGbMGKs5CosMn6w5wMvfJ3Eut4A/DW7P3Ze3dctPPar2aIFXHiErr4D//LiHt39Jpk6APxMua8Mdl7Smbh27Q0XuMA9+dfJJnl6wnR2HM+jXJoIp13TRi3MooOwCr4Osym2EBgXwl6s6MKpnNC9+m8QrP+zig5Up/GlwO8b2aemTLdXUU1n887skFmw6RIsGIbx+Uw+GdWmqZ6SqCtECr9xO26h6vHlzTzYcOM0L3+zkia+28e6v+7jvinb8PqEFQQHeX+jTTmfx35/28Nm6NPz9hIcGt+fey9sSEuRvO5ryIFrgldvq0bIhsyf0Y2nScV78Lom/fL6ZV77fxZ2XtmZsn5bUs9x14wqH07P57097mL02FUG4sW9L7hvYjqbhwbajKQ/kfb8hyquICFd0aMzA+CiW7T7Bm0v38tzXO3htyW7G9WvFH/u09Io59FvS0nlv+T4Wbj4EwOheMdx/RTuaN9A57arqdJBVeZzE1DO8uXQv328/QpGBS9tHMqZ3DEM7NaFOQM13YZw4cQKAyMiaXSgtK6+ARVuO8OmaA6xLOU3dIH9G947hjktaE93Q8/9oqdqhs2iUVzp0JpvP1qUxZ10qB89kE1E3iKu7NmNYl6b0aR3hlif9FBQWsWbfKRZsPsSCTYc5m1tAbKNQxvVrxejeMYQF66UOVeVogVderbDI8OueE8xee4Afdx4jJ7+IBqGBDO3YhCs7N6Vvm4hqFc6ZM2cCMH78+Cp9/dncAlbuPcmPO4/x/bYjnDyXR0igP8O7NmN0r2j6tI7QWTGqyqwVeBG5CngV8AfeMca8UNbrtcCr6srOK+TnXcf4dusRluw4RmZuAX4CnZuH07d1BH3bNKJnq4aVuh5pZefBp2fns/HAaTaknGbN/lOsTzlNfqEhNMifQR0ac3XXZgyMb6wzYlSNsDIPXkT8gf8CQ4E0YK2IzDfGbHfVPpUKCfLnqi7NuKpLM/IKili3/xSr9p1iVfJJPliVwju/7gOgSVgdOjYLI65JfWIiQmkVEUrzBsE0qluH8JBA/PxKb1EbYzibW8DJs3kcPJPNwdPZpJw6R9KRTJKOZpJ6KhsAP4FOzcO445I2XBYXSa9WET4xxVO5D1fOoukD7DHGJAOIyKfAtYAWeFUrggL8GNAukgHOq0jl5BeSmHqGLWnp7DicwfbDGazYe5K8gqLffJ2/n1A/OIDgAH/qBPqxOfUMBrj4hR/JyS8kPTufgguWOfb3E9pE1qVbdAPG9m5J95gGdItpYP0sXOXbXPm/rwWQWux+GtD3wheJyARgAkDLli1dGEf5uuBAf/q1aUS/No3+77GiIsOxzFwOnMricHo2J8/mcepcHhk5+eTmF5FTUMieOgEI0L9tI+oE+NEgNJAGIUE0rBtEiwYhRDcMoWl4sE+eaavcmysLfEmfcf+nw98YMx2YDo4+eBfmUep/+PkJTcODyzyRaNOb9QB46YZutRVLqRrhygKfBhRfQDsaOOTC/SnlEosWLbIdQakqceVnyrVAexFpLSJBwFhgvgv3p5RLhIaGEhqqJx4pz+OyFrwxpkBEHgC+wzFN8j1jzDZX7U8pV3n99dcBuO+++ywnUapyXDrEb4xZBOjnW+XR5syZA2iBV55Hh/2VUspLaYFXSikvpQVeKaW8lBZ4pZTyUm61mqSIHAdSgEjghOU4NcmbjkePxX150/HosVRcK2NMVElPuFWBP09E1pW2Opon8qbj0WNxX950PHosNUO7aJRSyktpgVdKKS/lrgV+uu0ANcybjkePxX150/HosdQAt+yDV0opVX3u2oJXSilVTVrglVLKS1kv8CLynogcE5GtxR57VkQ2i0iiiHwvIs1tZqyoko6l2HOPiYgRkUgb2aqilPfmaRE56HxvEkVkuM2MFVXaeyMiD4pIkohsE5F/2spXWaW8N7OLvS/7RSTRYsQKK+VYEkRklfNY1olIH5sZK6qUY+kmIitFZIuILBCRsFoLZIyxugGXAT2ArcUeCyt2+yHgTds5q3oszsdjcCybnAJE2s5ZzffmaeAx29lq6FiuABYDdZz3G9vOWZ3jueD5l4EnbeesxnvzPTDMeXs4sNR2zmocy1rgcuft24FnayuP9Ra8MWYZcOqCxzKK3a1LCZf6c0clHYvTv4C/4CHHcV4Zx+NxSjmWe4EXjDG5ztccq/VgVVTWeyMiAowGZtVqqCoq5VgMcL6lG46HXA2ulGOJB5Y5b/8AjKqtPNYLfGlE5O8ikgrcBDxpO09Vicg1wEFjzCbbWWrQA84utPdEpKHtMNUQB1wqIqtF5GcR6W07UA25FDhqjNltO0g1TARedNaAl4BJduNUy1bgGuftG/jtpUxdym0LvDHmb8aYGOBj4AHbeapCREKBv+HBf6BK8AbQFkgADuPoCvBUAUBDoB/wZ2COs/Xr6f6Ih7Tey3Av8LCzBjwMvGs5T3XcDtwvIuuB+kBebe3YbQt8MZ9Qix9palhboDWwSUT247jw+AYRaWo1VTUYY44aYwqNMUXA24BHDH6VIg2YZxzWAEU4FobyWCISAFwHzLadpZpuBeY5b3+GB/8/M8bsNMZcaYzpieMP797a2rdbFngRaV/s7jXATltZqsMYs8UY09gYE2uMicVRUHoYY45YjlZlItKs2N0/4Pj46am+BAYBiEgcEITnr2A4BNhpjEmzHaSaDgGXO28PAjy2u0lEGjv/9QMeB96srX279JqsFSEis4CBQKSIpAFPAcNFJB5HiyoFuMdewoor6ViMMR770bKU92agiCTgGATbD9xtK19llHIs7wHvOae05QG3GudUB3dXxv+1sXhY90wp781dwKvOTyQ5wAR7CSuulGOpJyL3O18yD5hRa3k85P+zUkqpSnLLLhqllFLVpwVeKaW8lBZ4pZTyUlrglVLKS2mBV0opL6UFXimlvJQWeKWU8lJa4JVyEpFYEdkpIu+IyFYR+VhEhojIchHZLSJ9nOvhP1bsa7aKSKzF2EqVSgu8Ur/VDngVuAjoANwIXAI8Bky2mEupStMCr9Rv7XOuIVQEbAOWOJcv2ALEWk2mVCVpgVfqt3KL3S4qdr8Ix9pNBfz29ya4lnIpVWla4JWqnP04LsmGiPTAsRy0Um5JC7xSlTMXiHBe0PpeYJfdOEqVTleTVEopL6UteKWU8lJa4JVSyktpgVdKKS+lBV4ppbyUFnillPJSWuCVUspLaYFXSikv9f8BRtzdPALe1z0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make scan of -lnL\n",
    "if not m.fixed[\"mu\"]:\n",
    "    plt.figure()\n",
    "    m.draw_mnprofile('mu', band=False, bound=(muHat-3.*sigma_muHat, muHat+3.*sigma_muHat), size=200)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "28a27e06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEFCAYAAAAYKqc0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgnUlEQVR4nO3dfXgV9Z338feXhxAS6A0pATFIojxZpICYaoGrQi21ukippV3kThVbadZKBamuxXqtYGuw3boKK2KvsCDWjUC1Sl26iqgQpYIaIApsRFAIiuHBtVYwtyDJ9/7jDGkeIUFyJmfyeV3XuTLzm9/MfHMIn0xm5vzG3B0REYmWNmEXICIip5/CXUQkghTuIiIRpHAXEYkghbuISAQp3EVEIqhd2AU0Vbdu3TwrKyvsMkREQrdx48YP3D29vmUJF+5ZWVkUFRWFXYaISOjMrLShZTotIyISQQp3EZEIUriLiERQXMLdzM4yszVmVmJm28xsetA+28z2mllx8PqHeNQjIhJ18bqgegy42d03mVlnYKOZrQ6W3efu98SpDhGRViEu4e7uZUBZMH3IzEqAjHjsW0SkNYr7OXczywLOB14Jmn5qZm+Y2WIz69rAOrlmVmRmRQcPHoxXqSIiCSuu4W5mnYA/Aje5+8fAg0AfYCixI/t/q289d89392x3z05Pr/d+fRERqSZu4W5m7YkFe4G7PwHg7vvdvcLdK4GFwIXxqkdEJMridbeMAYuAEne/t1p7z2rdrgS2xqMeEZGoi9fdMiOBq4EtZlYctP0CmGRmQwEHdgP/FKd6REQiLV53y6wDrJ5F/x2P/YuItDb6hKqISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBcRiSCFu4hIBCncRUQiSOEuIhJBCncRkQhSuIuIRJDCXUQkghTuIiIRpHAXEYkghbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwl4QzevRoRo8eHXYZIi1aXMLdzM4yszVmVmJm28xsetCeZmarzWxH8LVrPOoREYm6eB25HwNudvcvAV8FpprZQGAm8Ly79wOeD+ZFRORziku4u3uZu28Kpg8BJUAGMB54OOj2MPCdeNQjIhJ1cT/nbmZZwPnAK0APdy+D2C8AoHu86xERiaK4hruZdQL+CNzk7h83Yb1cMysys6KDBw82X4EiIhERt3A3s/bEgr3A3Z8ImvebWc9geU/gQH3runu+u2e7e3Z6enp8ChYRSWDxulvGgEVAibvfW23RU8DkYHoy8Kd41CMiEnXt4rSfkcDVwBYzKw7afgH8GviDmV0H7AG+H6d6REQiLS7h7u7rAGtg8TfiUYOISGuiT6hKQikoKGDDhg0UFhaSlZVFQUFB2CWJtEgKd0kYBQUF5ObmcuTIEQBKS0vJzc1VwIvUQ+EuCeP222+nvLy8Rlt5eTm33357SBWJtFwKd0kYe/bsaVK7SGumcJeE0bt37ya1i7RmCndJGHl5eaSkpNRoS0lJIS8vL6SKRFouhbskjJycHPLz8+nQoQMAmZmZ5Ofnk5OTE3JlIi1PvD7EJKfJ8YdUrF27NtQ6wpKTk8PChQuB1vseiDSGjtxFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBcRiSDd5y4JR/e3i5ycjtxFRCJI4Z5A9KAKEWmshDst88lnsHZ32FXE33MrCrjnFzUfVHHdj3MpOQhjvqOxVURaI+vQObXBZe4ez1o+twGDs/3x1UVhlxF338rOouy90jrtPXtlsqpod/wLEpHQDTm7647K8r/2r2+ZTsskiH17638gRUPtItK6KdwTxBkZ9T+QoqF2EWndFO4JYtpteSR3rPmgiuSOKUy7TQ+qEJG6FO4JYuyEHGbdk09SUuxBFT17ZTLrnnzGTtDFVBGpK+HulmnNxk7I4Y//GXtQxeIn14ZbjITmR1eOBvQzICeWcHfLZJ9jXnRX2FWIiISvzZQuultGRKQ1SbjTMoc6XcCWb7S++9yP05/kop8B+buuDS7RkbuISAQp3EVEIkjhLiISQQp3EZEISrgLqq2dLqKJSGPE5cjdzBab2QEz21qtbbaZ7TWz4uD1D/GoRSSR/fmPBbyxcQNF6wv5VnYWf/6jxvSX+sXryH0JMB/4fa32+9z9nqZsqI3B//6/01WWSOI4Pqb/0aOxMf3L3itl9i25HD6qMf1brcqKioYWxe0TqmaWBax090HB/GzgcFPDPTs724uKWu997tJ6ZWVlUVpad0z/zMxMdu/eHf+CJHRmttHds+tbFvYF1Z+a2RvBaZsG78Y3s1wzKzKzooMHD8azPpEWY8+e+sfub6hdWrcww/1BoA8wFCgD/q2hju6e7+7Z7p6dnp4ep/JEWpbevesfu7+hdmndQgt3d9/v7hXuXgksBC4MqxaRRJCXl0dKSs0x/VNSUsjL05j+Uldo4W5mPavNXglsbaiviEBOTg75+fl06BAb0z8zM5P8/HxycnQxVeqKy90yZrYUGA10M7P3gFnAaDMbCjiwG/ineNQikshycnJYuDA2pv/atWvDLUZatLiEu7tPqqd5UTz2LSLSGoV9t4yIiDQDhbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIL0sA6RBKMPL0lj6MhdRCSCThruZvacmQ2JRzEiInJ6NObI/VbgPjN7qNZgXyIi0kKdNNzdfZO7XwKsBJ4xs1lm1rH5SxMRkVPVqHPuZmbAdmIP2LgR2GFmVzdnYSIicuoac859HbAXuA/IAK4lNnzvhWaW35zFiYjIqWnMrZDXA9u87pO0bzSzkmaoSUREPqeThru7n+gJSWNPYy0iInKafK773N39ndNViIiInD76EJOISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhEUl3A3s8VmdsDMtlZrSzOz1Wa2I/jaNR61iIi0BvE6cl8CXFarbSbwvLv3A54P5kVE5DSIS7i7+4vAh7WaxwMPB9MPA9+JRy0iIq1BmOfce7h7GUDwtXuItYiIREpCXFA1s1wzKzKzooMHD4ZdjohIixdmuO83s54AwdcDDXV093x3z3b37PT09LgVKCKSqMIM96eAycH0ZOBPIdYiIhIp8boVcimwHhhgZu+Z2XXAr4FvmtkO4JvBvIiInAbt4rETd5/UwKJvxGP/IiKtTUJcUBURkaZRuIuIRJDCXUQkghTuIiIRpHAXEYkghbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBcRiSCFu4hIBCncRUQiSOEuIhJBCndpMR5//HHMrMnrderUiSVLlpz+gkQSmMJdWp3Zs2czaNCgsMsQaVYKdzmtjh07hruHXYZIq6dwr0dBQQFZWVm0adOGrKwsCgoKmnV/o0eP5ic/+Qk333wzaWlppKenM2/ePI4cOcLUqVPp0qULvXv35pFHHqmx3t69e7nqqqvo2rUrXbt2ZezYsezYsaNq+dtvv8348eM544wzSE1NZdiwYaxcubLGNp544gkGDx5Mx44dSUtLY9SoUezfvx+o/wh3yZIldOrUqWr+eJ8lS5bQp08fOnTowCeffMLf/vY3cnNz6d69O507d2bUqFEUFRXV2Nbvf/97MjMzSUlJ4Yorrqja74ns3LmT0aNHk5yczIABA+p8PwAzZ85kwIABdOzYkaysLG699VY+/fTTqvrvvPNOtm3bhplhZlWndO69914GDx5MamoqGRkZTJkyhY8++uikNYm0RAr3WgoKCsjNzaW0tBR3p7S0lNzc3GYP+IKCAjp37swrr7zCzJkzuemmm/jOd75D//79KSoqYvLkyUyZMoX3338fgPLycr7+9a+TnJxMYWEh69evp2fPnowZM4by8nIADh8+zOWXX87q1at5/fXXmTBhAt/97nd58803Adi3bx9XXXUVkydPpqSkhBdffJGrr766ybXv2rWLRx99lMcee4zXX3+dDh06MHbsWPbu3cvKlSvZvHkzF198MZdccgllZWUAvPLKK1x77bXk5uZSXFzMuHHjuOOOO064n8rKSq688koqKytZv349ixcvZvbs2Rw5cqRGv9TUVBYvXkxJSQkLFixg2bJl5OXlATBx4kRuvvlmBgwYQFlZGWVlZUycOBGANm3aMHfuXLZt28ajjz7Kq6++yo033tjk90OkRXD3hHpdcMEF3pwyMzMdqPPKzMxstn2OGjXKv/rVr1bNV1ZWerdu3XzcuHFVbUePHvX27dv7Y4895u7uixYt8r59+3plZWVVn2PHjnlaWpovX768wX1ddNFF/qtf/crd3Tdu3OiA7969u96+s2bN8vPOO69G20MPPeSpqak1+rRr18737dtX1fb88897amqql5eX11h3yJAh/pvf/Mbd3SdNmuRjxoypsfy6667z2I9k/VatWuVt2rTx0tLSqraXXnrJAX/ooYcaXO/BBx/0Pn36nPD7qs/TTz/tSUlJXlFRcdK+ImEAiryBrGwX2m+VFmrPnj1Naj9dBg8eXDVtZnTv3p0vf/nLVW3t27ena9euHDhwAICNGzeya9cuOnfuXGM75eXlvP322wB88skn3HnnnaxcuZKysjI+++wzPv3006p9DRkyhDFjxjBo0CAuvfRSxowZw/e+9z3S09ObVHuvXr3o0aNH1fzGjRspLy+vs51PP/20qraSkhLGjRtXY/nw4cNZtGhRg/spKSkhIyOD3r17V7VddNFFtGlT8w/Qxx9/nLlz57Jz504OHz5MRUUFFRUVJ/0+XnjhBe6++25KSkr429/+RkVFBUePHmXfvn2ceeaZJ11fpCVRuNfSu3dvSktL621vTu3bt68xb2b1tlVWVgKxUxRDhw5l2bJldbaVlpYGwC233MIzzzzDPffcQ79+/UhJSeGaa67h6NGjALRt25Znn32WDRs28Oyzz7Jo0SJuu+02CgsLGTJkCG3atKlzcfSzzz6rs7/U1NQa85WVlfTo0YOXXnqpTt8vfOELAKd00bUx62zYsIGrrrqKWbNmcd9999GlSxeeeuopbrnllhOuV1paytixY/nxj3/ML3/5S774xS+yadMmJk2aVPV+iSQShXsteXl55ObmVp23BkhJSak6Z9tSDBs2jKVLl9KtWze6dOlSb59169ZxzTXXMGHCBODvR879+/ev6mNmDB8+nOHDh3PHHXdw3nnnsXz5coYMGUJ6ejr79+/H3avuPy8uLm5Ubfv376dNmzacc8459fYZOHAgGzZsqNFWe76+dfbu3cu7777LWWedBcCrr75a9QsP4C9/+QsZGRn8y7/8S1Vb7V/WSUlJdY7ki4qKOHr0KPfddx9t27YFqPdirUiiCP2CqpntNrMtZlZsZkUnX6N55eTkkJ+fT2ZmJmZGZmYm+fn55OTkhF1aDTk5OfTo0YPx48dTWFjIrl27ePHFF7n55pur7pjp378/Tz75JJs2bWLLli384Ac/qLprBGJhetddd/Haa6+xZ88ennrqKd59910GDhwIxO7i+fDDD5kzZw5vv/02ixYt4vHHHz9pbWPGjGHkyJGMHz+ep59+ml27drF+/XpmzZpVdTQ/bdo0nnvuOe6++2527NjBwoULefLJJ0+63XPPPZdrrrmG4uJi1q9fz4wZM2jX7u/HKP3792fv3r0UFBTwzjvv8OCDD7J06dIa28nKyqK0tJRNmzbxwQcfcOTIEfr160dlZSVz585l165dLF26lLlz5zbq30KkRWroZHy8XsBuoFtj+zf3BdUwjBo1yqdOnVqj7bzzzvNZs2bVaOvRo4fff//9VfP79u3za6+91tPT0z0pKcmzsrL8hz/8oR88eNDd3Xfv3u3f+MY3PCUlxTMyMvy3v/2tjx071idPnuzu7v/zP//jl112mXfv3t2TkpK8T58+VRc8j/vd737nvXv39pSUFJ84caLPnTu3zgXV+i5Ofvzxxz5t2jTPyMjw9u3be69evXzixIm+c+fOqj6LFy/2s846y5OTk/2yyy7z+++//4QXVN3dt2/f7hdffLEnJSV53759/U9/+pOnpqbWuKA6c+ZM79atm6empvqVV17pCxYsqLHdTz/91CdMmOBdunSpcTF23rx5fuaZZ3pycrJfcsklvnz5cgd8165dJ6xJJCyc4IKqecgfODGz3UC2u3/QmP7Z2dle+35pEZHWyMw2unt2fctCPy1D7FbDZ81so5nl1tfBzHLNrMjMig4ePBjn8kREEk9LCPeR7j4MuByYamYX1+7g7vnunu3u2U29TU9EpDUKPdzd/f3g6wHgSeDCcCsSEUl8oYa7maWaWefj08ClwNYwaxIRiYKw73PvATwZ3EPdDnjU3Z8JtyQRaelGjx4NwNq1a0OtoyULNdzd/R1gSJg1iIhEUejn3EVE5PRTuIuIRJDCXUQkghTuIiIRpHBvgQ4ePMjs2bPRp3FF5FQp3Fugn/zkJxQVFTF16tSwSxGRBKVwb2EeffRROnTowMqVK2nfvj1/+MMfwi5JRBJQ6KNCNpVGhRQRfYgppqWPCilx9OKLL/Ltb3+bjIwMzIwlS5bU26+srIzJkyeTnp5OcnIyAwcOpLCwsN6+d999N1/5ylf4whe+QHp6OuPGjWPr1pqjSDR2v6fDggULOPvss0lOTuaCCy6o93F/tR06dIibbrqJzMxMOnbsyIgRI3jttdea3EeaX0FBARs2bKCwsJCsrCwKCgrCLqlFUri3MocPH2bQoEHMmzePjh071tvno48+YuTIkbg7f/7znykpKeH++++ne/fu9fZfu3YtN9xwAy+//DIvvPAC7dq1Y8yYMXz44YdN2u+JXHvttcyePfuk/ZYvX8706dP5xS9+webNmxkxYgSXX375SR9wPmXKFFatWsXDDz/Mli1bqh4Yvnfv3ib1keZVUFBAbm4uR44cAWKPUMzNzVXA16ehp3i01FcUn8Tk7v7yyy+7mfmhQ4eq2j788EMHvLi4uFn2WfsJRsfddtttPmLEiFPe7qFDh7xNmzb+1FNPNWm/JzJ58uQ6T6aqz4UXXuhTpkyp0da3b1+fOXNmg+uUl5d727ZtfcWKFTXahw0b5rfffnuj+0jzy8zMdGLPgKjxyszMDLu0UHCCJzHpyL2FKC4upl+/fnTq1KmqbfPmzSQlJVU90/S4OXPm0KlTpxO+GnMqoiErVqzgoosuYuLEiXTv3p2hQ4cyf/78449FPKlDhw5RWVlJ165dT7mGU3H06FE2btzIpZdeWqP90ksv5eWXX25wvWPHjlFRUUFycnKN9o4dO7Ju3bpG95Hm19BfYCf7y6w1CntUyOb3qIWz3//btAvVxcXFDBs2rEbb5s2bGThwIO3bt6/Rfv311/OP//iPJ9xeRkZGk/Zf3TvvvMOCBQuYMWMGM2fOpLi4mBtvvBGAn/70pyddf/r06QwdOpThw4efcg1z5sxhzpw5VfNHjhzBzLjnnnuq2p5++mm+9rWvVc1/8MEHVFRU0KNHjxrb6tGjB88991yD++rcuTPDhw/nrrvuYtCgQZxxxhksXbqU9evX07dv30b3kebXu3dvSktL622XmqIf7gmiuLiYCRMm1GjbtGkTQ4cOrdM3LS2NtLS0ZqulsrKS7Oxs7r77bgDOP/98duzYwQMPPHDScP/Zz37GunXrWLduHW3btj3lGmr/Avv5z39ORkYG06ZNq2pr6BdYMIR0FXev01bbI488wo9+9CN69epF27ZtGTZsGJMmTWLTpk1N6iPNKy8vj9zcXMrLy6vaUlJSyMvLC7Gqlin64d7EI+gwVFZWsnXr1jo/oEVFRdxwww11+tc+qq1P7aPapujZs2edU0Ff+tKXmDdv3gnXmzFjBsuWLWPNmjWcc845p7Tv42r/AuvcuTNpaWknPEru1q0bbdu2Zd++fTXaDxw4UOdovrY+ffpQWFjIJ598wscff0zPnj2ZOHEiZ599dpP6SPPKyckB4LrrruPIkSNkZmaSl5dX1S5/F/1wTwDbt2+nvLycM888s6pty5YtvPXWW/UeuTf3aZmRI0eyffv2Gm1vvfUWmZmZDa4zffp0li1bxtq1azn33HNPed+fR1JSEhdccAGrV6/m+9//flX76tWr6/xV1JDU1FRSU1P561//yqpVq/jXf/3XU+ojzScnJ4eFCxcCus/9RBTuLUBxcTEA8+fPZ8aMGezevZubbroJoOqWr+o+z2mZw4cPs3PnTiD2F8OePXsoLi4mLS2t6rzljBkzGDFiBHl5eUycOJHNmzfz7//+7zX+Wpg/fz7z58/nzTffZOrUqTzyyCOsWLGCrl27Vh05H7+429j91q7z8OHDVfO//vWvAWoclaelpZGUlFRjvZ/97GdcffXVXHjhhYwcOZLf/e53vP/++1x//fX11n7cqlWrqKys5Nxzz2Xnzp388z//MwMGDOCHP/xhk/qItBgN3UbTUl9RvBXy1ltv9W9+85s+btw4T0pK8kGDBvmKFSu8a9eufsUVV5zWfa1Zs6beW8kmT55co9/KlSt98ODB3qFDB+/Xr5/PmzfPKysrq5bPmjXLYz8+Xu/2gBq3LjZ2v7W3f6LXmjVr6l33gQce8MzMTE9KSvJhw4Z5YWFhvduubvny5X7OOed4UlKSn3HGGT516lT/6KOPmtxH4mPUqFE+atSosMsIHSe4FVLDD7QA3/rWtxg2bFjVBUwROTENPxCj4QdauNdff53BgweHXYaIRIjCPWT79+9n//79CncROa10QTVkPXr0aPQnP0VEGktH7iIiEaQjdxFJOK39Qmpj6MhdRCSCFO4iIhGkcBcRiSCFu4hIBCncRUQiSOEuIhJBCncRkQhSuIuIRFDo4W5ml5nZdjPbaWYzw65HRCQKQg13M2sLPABcDgwEJpnZwBOvJSIiJxP2kfuFwE53f8fdjwLLgPEh1yQikvDCDvcM4N1q8+8FbSIi8jmEPXCY1dNWZ/xbM8sFcoPZw2a2vXafVqYb8EHYRYSstb8Hrf37B70HAA0+tT7scH8POKvafC/g/dqd3D0fyI9XUS2dmRU19Git1qK1vwet/fsHvQcnE/ZpmdeAfmZ2tpklAVcBT4Vck4hIwgv1yN3dj5nZT4FVQFtgsbtvC7MmEZEoCPu0DO7+38B/h11HgtEpKr0Hrf37B70HJ2R6fqeISPSEfc5dRESagcI9gZjZYjM7YGZbw64lDGZ2lpmtMbMSM9tmZtPDrinezCzZzF41s9eD9+DOsGsKg5m1NbPNZrYy7FpaKoV7YlkCXBZ2ESE6Btzs7l8CvgpMbYXDVRwBLnH3IcBQ4DIz+2q4JYViOlASdhEtmcI9gbj7i8CHYdcRFncvc/dNwfQhYv+5W9Unmj3mcDDbPni1qgtnZtYLGAv8R9i1tGQKd0lIZpYFnA+8EnIpcReckigGDgCr3b21vQdzgVuBypDraNEU7pJwzKwT8EfgJnf/OOx64s3dK9x9KLFPdF9oZoNCLiluzOwK4IC7bwy7lpZO4S4JxczaEwv2And/Iux6wuTuHwFraV3XYUYC3zaz3cRGkb3EzP4z3JJaJoW7JAwzM2ARUOLu94ZdTxjMLN3MugTTHYExwJuhFhVH7n6bu/dy9yxiw5W84O4/CLmsFknhnkDMbCmwHhhgZu+Z2XVh1xRnI4GriR2tFQevfwi7qDjrCawxszeIjc202t11O6DUoU+oiohEkI7cRUQiSOEuIhJBCncRkQhSuIuIRJDCXUQkghTuIs3IzEab2Yg47zOrtY4cKn+ncJdWw8zCePLYaCCu4f55hfQ+yWmmcJdmERw9vmlm/2FmW82swMzGmNlfzGyHmV0Y9EsNxql/LRife3y19V8ys03Ba0TQ3tPMXgw+wLTVzL4WtB+utu/vmdmSYHqJmd1rZmuA35hZHzN7xsw2Bts/t1q/B4Px4t8xs1FBXSXHtxX0u9TM1gc1PRaMc4OZ7TazO4P2LWZ2bjC42fXAjKDer9V6j2YH+1gb7HNate99a7V+t5jZ7GB6rZndF7wHJWb2FTN7InhP76q2+XZm9rCZvWFmj5tZSrD+BWZWGHz/q8ysZ7XtzjGzQmLD6Uqic3e99DrtLyCL2PjrXyZ2ELERWAwYMB5YEfSbA/wgmO4CvAWkAilActDeDygKpm8Gbg+m2wKdg+nD1fb9PWBJML0EWAm0DeafB/oF0xcR+/j68X7LqtX3ca3ahwLdgBeB1GCdnwN3BNO7gRuD6RuA/wimZwO3NPAezQZeBjoE2/5fYkP4ZgFbq/W7BZgdTK8FfhNMTwfeJ/ap1Q7Ae8AXg/UdGBn0Wxxso32wv/SgfSKxh9If3+6CsH9u9Dp9L/35Jc1pl7tvATCzbcDz7u5mtoVYAAFcSmwgqFuC+WSgN7HQmm9mQ4EKoH+w/DVgcTCA2Ap3L25EHY+5e0VwlD0CeCw2TA0QC8Xj/qtafftr1Z5FbBTGgcBfgvWTiA0Hcdzxgcw2At9tRF0Af3b3I8ARMzsA9GjEOk8FX7cA29y9LKjzHeAs4CPgXXf/S9DvP4FpwDPAIGB1UH9boKzadpc3smZJAAp3aU5Hqk1XVpuv5O8/ewZMcPft1VcMTkPsB4YQO3r+FGIPLDGzi4k9rOERM/utu/+emg+sSK5VxyfB1zbARx4bLvdE9VavtXq9FcTGcpl0kvUraPz/rer7Ob7eMWqeMq39/ZysTqj7AA8n9l5vc/fhDdTySQPtkoB0zl3Ctgq4MRjxETM7P2j/P0CZu1cSGyysbbA8k9h43guJjRA5LOi/38y+ZGZtgCvr25HHxn7fZWbfD7ZlZjakCbVuAEaaWd9g/RQz63+SdQ4BnZuwD4j9UutuZl80sw7AFU1cH6C3mR0P8UnAOmA7kH683czam9l5p7BtSQAKdwnbr4idC34juIj4q6B9ATDZzDYQOyVz/KhyNFBsZpuBCcC8oH0msXPrL1DzVENtOcB1ZvY6sI3Y+fVGcfeDwLXAUouNyrgBOPckq/0XcGV9F1RPsJ/PgF8Se8rUSk5tSN8SYu/fG0Aa8KC7HyV2PeI3wfdfTILdySONp1EhRUQiSEfuIiIRpHAXEYkghbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIL+Px6ICLQu1AVmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot fit\n",
    "fig, ax = plt.subplots(1,1)\n",
    "plt.errorbar(x, y, yerr=s, xerr=0, color='black', fmt='o', label='measured data')\n",
    "plt.xlabel(r'measurement number')\n",
    "plt.ylabel(r'$y$', labelpad=5)\n",
    "xMin = 1.e-4\n",
    "xMax = 5. - 1.e-4\n",
    "yMin = 0.\n",
    "yMax = 26.\n",
    "plt.xlim(xMin, xMax)\n",
    "plt.ylim(yMin, yMax)\n",
    "xPlot = np.array([xMin, xMax])\n",
    "fit = np.array([muHat, muHat])\n",
    "plotLabel = r'$\\hat{\\mu} = ' + '{:.2f}'.format(muHat) + '\\pm' + '{:.2f}'.format(sigma_muHat) + '$'\n",
    "plt.plot(xPlot, fit, 'orange', linewidth=2, label=plotLabel)\n",
    "fitUp = np.array([muHat+sigma_muHat, muHat+sigma_muHat])\n",
    "fitLo = np.array([muHat-sigma_muHat, muHat-sigma_muHat])\n",
    "plt.fill_between(xPlot, fitLo, fitUp, color='dodgerblue', alpha=0.2)\n",
    "plt.subplots_adjust(left=0.15, right=0.9, top=0.9, bottom=0.15)\n",
    "handles, labels = ax.get_legend_handles_labels()  # Tweak legend\n",
    "handles = [handles[1], handles[0]]\n",
    "labels = [labels[1], labels[0]]\n",
    "handles = [handles[0][0], handles[1]]    # turn off error bar for data in legend\n",
    "plt.legend(handles, labels, loc='lower left', fontsize=14, frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1197fcb",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}