{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stave.py\n",
    "# G. Cowan / RHUL Physics / December 2020\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",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import truncexpon\n",
    "from scipy.stats import truncnorm\n",
    "from iminuit import Minuit\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 10.                                  # initial value for fit\n",
    "y = np.array([8., 9., 10., 11., 12.])     # measured values\n",
    "s = np.array([1., 1., 1., 1., 1.])        # estimates of std. dev\n",
    "v = s**2                                  # estimates of variances\n",
    "r = np.array([0.2, 0.2, 0.2, 0.2, 0.2])   # relative errors on errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def negLogL(mu):\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,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize Minuit and set up fit:\n",
    "#     initial parameter values are guesses,\n",
    "#     error values set initial step size in search algorithm,\n",
    "#     limit_param to set limits on parameters,\n",
    "#     fix_param=True to fix a parameter,\n",
    "#     errordef=0.5 means errors correspond to logL = logLmax - 0.5,\n",
    "#     pedantic=False to turn off verbose messages.\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.from_array_func(negLogL, parin, parstep, name=parname,\n",
    "    limit=parlim, fix=parfix, errordef=0.5, pedantic=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "par index, name, estimate, standard deviation:\n",
      "   0 mu          =  10.489128  +/-  0.641382\n",
      "\n",
      "free par indices, covariance, correlation coeff.:\n",
      "0 0 0.411371 1.000000\n"
     ]
    }
   ],
   "source": [
    "# do the fit, extract results\n",
    "m.migrad()                                             # minimize -logL\n",
    "MLE = m.np_values()                                    # max-likelihood estimates\n",
    "sigmaMLE = m.np_errors()                               # standard deviations\n",
    "cov = m.np_matrix(skip_fixed=True)                     # covariance matrix\n",
    "rho = m.np_matrix(skip_fixed=True, correlation=True)   # correlation coeffs.\n",
    "npar = len(m.np_values())\n",
    "nfreepar = len(cov[0])\n",
    "\n",
    "npar = len(m.np_values())\n",
    "print(r'par index, name, estimate, standard deviation:')\n",
    "for i in range(npar):\n",
    "    if not(m.is_fixed(i)):\n",
    "        print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n",
    "         \"{:.6f}\".format(MLE[i]), \" +/- \", \"{:.6f}\".format(sigmaMLE[i]))\n",
    "\n",
    "print()\n",
    "print(r'free par indices, covariance, correlation coeff.:')\n",
    "for i in range(nfreepar):\n",
    "    for j in range(nfreepar):\n",
    "        print(i, j, \"{:.6f}\".format(cov[i,j]), \"{:.6f}\".format(rho[i,j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEiCAYAAADjxEWuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5wU9f3H8dfnGgcc7ShH5ygC0hUERNQDC2rsvbcYGxr9GRNLTMSILcYeG1ExatTYFbuioChKEQSU3nvn4KhXvr8/Zs6sy97dXtvZW97Px2Mfy818Z+bzne+yn535znfGnHOIiIiESwo6ABERiU9KECIiEpEShIiIRKQEISIiESlBiIhIREoQIiISkRKEiIhEpAQhMWVmXczsQTObYGY7zMyZ2cBSyvczsy/MLM/MNprZC2bWLMptPe+vP/y1oupqtNc2zzOzGWa2y8yWmtkdZpZajuXrmdl9ZrbIX8dKM3vDzOqUsszgkLo1D5vXwszuNbOxZpbrlzm7MnWsLDM72sy+M7OdZrbGzB41s4xyLJ9mZreY2Wx/H601sw/MrHUpy3T0y+71eTOzDL+dPjSz9X6ZmytTx0SREnQAss85GLgOmAvMBPqXVNDMugHjgCXATUBD4A/AgWZ2kHNuZxTbywcuDZu2vdxRR8HMLgKeBz4AHgN6ArcBbSLEEGn5BsB4oDUwClgANAUOBWoBOyIsk+RvaztQN8Jqu+Dtu0XAdOCw8tWqapnZUOBD4HvgeqAtcAPQ3cyOdGWM3PWT7fvAIOAZvM9QQ2AA0AAoKfk/BBTg7cdwTYC/+stOA44qX60SmHNOL71i9gIygfr+vy8GHDCwhLLvAuuAzJBpOf4y10exreeBXTGqVy0/1s/Dpo/w4z0ginU8AWwCssux3auADcDD/naah82vBzQO23dnV7KuzwPjKrjsj8DPQFrItOLPwclRLP8nYDfQrxzbPMZf5s5Inze/7Vr6/872y9wci89NvL90iinOmdkI/5C3i3/KZLN/quU+M0sys2Zm9pqZbfEPj28PWz7HXz4nwrqdmY2IVV0AnHObnHNbyypnZvWAY4GXnXObQpYfB8wCzop2m/5+qm9mVoGQozUE79f+P8OmP473hXNmaQubWUPgEmCUc26Jfxol0q/d0GUy8b70/gpsiVTGObfNObcxuipULzPrCvTCq+OekFkv4cVfapv6R0vXAW8756aYWUppp978ZVLxkucjwMJIZZxzu51zq6Kvyb5DCaLmeBWoDdwKfI33S+pPwKfANuAWvF9nI8zs5KrcsJk1MrMmUbyiPo8chZ5AKjApwrzvgT7+F0ZZ0oCtQC6w2cxG+V/GVe1A//1X8Trn1uOd3jlwryV+bTCQDiw0s7eBncBOv6+mVwnLjARWA09XOOrYKmkfFQBTKXsfdQNaAjPN7Fm802rbzezHSD+AfP8HNMLbV1JO6oOoOX5wzv0WwMyewjs/fTcw0jn3V3/6aGAV8FvgnSrc9jSgXRTl/o13uqAqtPTfV0eYtwrvy7QRUNqv49XA3/HiBxgGXAb0NbNBzrndVRQr/C/eNRHmrQqZX5L9/Pd7gPnA+Xjn1G8HvjCzns65X/aFmfUGLgeOdc4VVu/BUZUpq00HlLF88T76P7zTalfgHZ3dCnzs90vNLC5sZi3w+oBucM5trSH7KK4oQdQczxT/wznnzGwS0AF4LmT6LjP7EehYxds+D+/opSxVeZhevL1IX+K7wspE5Jy7JWzSf81sDnAfcAEh+7QK1AbynXNFEebtouz9V3z0VQQc6ZzbDmBmU/F+cV+L90VY7DHgA+fcZ5WKugz+KZoGYZNrAalm1iRseq5zLr+U1ZXVptHuo3pAH+fcCj/GsXinj24Bzg0p/3d/+nNIhShB1BzLwv7O9d+XR5jetSo37Jz7pirXF6XiK5QinYdPDytTHg/jnW44klIShJnVJuyL0TkX6eig2E68L01zfm9niPQoYi2eP6Y4OfjbnGxmi/BOQRXHdg4wEOhexjqrwiHAlyXMWx/29xC8q85KUlabRruPJhQnBwDn3Aozm8Cv99EgvB82Q0pI2hIFJYiaozDSROdcpOmhx9IRLxs0s+RoN2xmTYFoyu90zuWWXSwqxachWkSY1xLvF+fm8q7UObfHzNbgXU1VmrOA0WHTSjtHERpv+JFUS0roIA1RvMy6CPPW4Z1OK3Y/8DrgzKyTP624PtlmluacC/9BUVE/svdln38EmuNdchxetjSh+2hx2LyWlH0EWtY+Cr1k+u94fXUrQ/ZR8fiZ1mbWwTm3qIzt7fOUIBJf8ZdoeMdsdjnWMZnY90HMxBvD0B94OWzeAGB6RX4Zmlk63hfU12UU/YTyXQ8/1X/vT0j/j38apgPeF3o0y7eKMK8VXp9T6N/n8uvTKcUm4o0xqZKjSOfcZuDz0Glmdj5Qyzn3eeSlShS6j74NWV8KXgf1J2UsX/yZKGkfhR7RtMb7zM6PUPZ1vNNc6RHmSQgliMS3BO/oYwi/7ri+phzriHkfhN+p+Alwrpnd4X9R4V+t0oOwX6/+JZQ7in85+5eIJjvnwgeX3Yr3uf+4jO2vJnJnaknG4XWcDmfv/WyEJAj/0sy2wAbn3AZ/e3P9/qMTzSyz+NJeMzsSb6DdkyHrPCPC9s/0p/+OvU9HxgXn3BwzmwVcbmZPhFzqej7eEVLoPkrF60vLLe6cd85tM7MPgePNrItzbq5ftgvewLlnQzZ3FXsPHMzBa59bgTlVXb9EpASR4Pwv2leA4Wbm8H5dDsH7VRvtOqqsD8IfLXyt/+cB/vvF/hfhFudc6DiCW/Auaf3azJ7E6xO4EZgNPBW26tl4o5Bz/L9bABPN7C1gHt6ptiOBE/xyr1RVneCXCwRuAp41szF4g/x64X0hveCc+yGkeH+88/p34A2kK/Z/wGfAN2Y2CqiPN8p4Ad54iuJtvRG+fTPr4f/z/fC+EjO7zf9ne//9lOLTLs65WF/+eQNecv7CzF7AS5R/wEuwb4eUa4XXpuFHpbcCR/jLP+pP+z3eOIpf6uKc+yh8wyGXYX/pnPsubN41eEfZxUfaQ/wjG4DHqvDUac0S9Eg9vUp/8b+RuOEjZJ/ymm+v8u8Aa8KmZeKNo9iG14n9Mt6gLgeMiHF9sv3tRnotiVC++Mt0O94o45fC94VfzhEyuhfvP/oLeMkhD6/P4me8QWXp1Vi/C/BOhezGu4DgTkJGDftlckra93hJbCJeh2xxfVtU9HMSsm8ivipYx+ep4Ehqf/lj8K7M2gmsxRtcWK+Ez8nzEZY/EC+R5uGNcXkX6BLFdi+mhJH7eEfaJe2n7Fj+H4mnl/k7R0RE5Fc0klpERCJSghARkYiUIEREJCIlCBERiShhLnNt0qSJy87OrtCy27dvp27dSM9aqXlUl/iTEPWYOxeA7a1b1/y6+BKiXXyVqcvUqVM3OOeaRpqXMAkiOzubKVOmVGjZcePGkZOTU7UBBUR1iT8JUQ8//nEjRtT8uvgSol18lamLmS0taZ5OMYmISERKECIiEpEShIiIRKQEISIiESlBiIhIREoQIiISkRKEiIhEtM8niIkLN/LBoj1lFxQR2cfs8wniizlreWNePvPXbgs6FBGRuLLPJ4ircjpRKxke+HRe0KGIiMSVfT5BZNZN49j2qXz80xp+XL4l6HBEROJGzBKEmd1iZpPNbKuZrTezMSHP0Y1UfpSZOTO7sbpjOzo7lcZ107j/k7nVvSkRkRojlkcQOcATwCBgKFAAfG5mmeEFzex04CBgVSwCq51iXD2kExMWbODbBRtisUkRkbgXswThnBvmnBvtnJvlnJuJ93D3psAhoeXMrB3wCHAukB+r+M4b0JaWDdK575O56DndIiLB3u67Hl6C2lw8wcxSgFeAkc652WZW6grM7HLgcoCsrCzGjRtXoUDy8vL47puvGdamiNGztvDga2Ppm1Uz74Sel5dX4f0QbxKlLolQjz5bvP65RKhLMdWlbEF+Cz4CTAcmhky7A9jonHsymhU450YBowD69evnKno/9OJ7qQ8uLGL8w1/x8Urj+jMOIzmp9AQVj3SP+/iTEPVo2BCAjIyMml8XX0K0i6+66hLIVUxm9iAwGDjNOVfoTzscuBi4NIiYAFKSk/jDUV2Yvy6Pd6atDCoMEZG4EPMEYWYPAecAQ51zi0JmDQFaAKvNrMDMCoB2wH1mtiJW8R3bozk9WtXnoc/nsaegKFabFRGJOzFNEGZW3Pk81Dk3J2z2E0AvoE/IaxXwEHBErGJMSjL+OKwrKzbv5NXJy2K1WRGRuBOzPggzexzvyqWTgc1m1tyfleecy3POrQPWhS2TD6xxzsV0gMJh+zVhQPtMHh27gNP7tqZOWs3ssBYRqYxYHkFcjXfl0lhgdcir2gfClZeZ8adjurIhbzejv1kSdDgiIoGI2U9j51y5LwlyzmVXQyhR6duuEUfu34ynxy/k/AHtaFAnNahQREQCsc/fi6k0Nw7rwrbdBTw5fmHQoYiIxJwSRCm6Nq/PKX1aMfqbxazasjPocEREYkoJogw3HN0ZBzz4mW4HLiL7FiWIMrRuVIdLBmXz5g8rmL16a9DhiIjEjBJEFK7O6UT99FTu/Sh86IaISOJSgohCgzqpXDu0E+PnrWfCfN0OXET2DUoQUbrg4Ha0blSbez6aTVGRbgcuIolPCSJKtVKSufHoLvy0aivv/RiT5xiJiARKCaIcTuzdku4t63P/J3PZXVAYdDgiItVKCaIckpKMW4/bn5VbdvLixKVBhyMiUq2UIMrpkE5NOLxzUx77YgG5O2L2RFQRkZhTgqiAm4/tytZd+TwxbkHQoYiIVBsliArYv0V9TjuwNaO/XcKKzTuCDkdEpFooQVTQDUd1xoAHP9UtOEQkMSlBVFDLhrW5dHB73p6+klkrc4MOR0SkyilBVMJVOR1pVCeNO9//Gec0eE5EEosSRCXUT0/lhqM68/3iTXzy05qgwxERqVJKEJV09kFt6JyVwd0fztHgORFJKEoQlZSSnMRtv+nGsk07+Pe3S4IOR0SkyihBVIHDOjdlaNdmPDZ2ARvydgcdjohIlVCCqCK3Hrc/O/MLeUhPnhORBKEEUUU6Ncvg/IHteGXSMuas0ZPnRKTmU4KoQtcfuR/10lO564PZuuxVRGo8JYgq1LBOGtcdsR9fz9/Al3PXBR2OiEilKEFUsQsObkeHpnUZ+cFs8guLgg5HRKTClCCqWGpyEn8+bn8Wrd/OS9/pmREiUnMpQVSDoV2bMbhTEx7+fD5bduwJOhwRkQpRgqgGZsZtx+/Ptl35PPz5/KDDERGpECWIatK1eX3O6d+WF79byry124IOR0Sk3JQgqtGNR3cho1YKt7/7ky57FZEaJ2YJwsxuMbPJZrbVzNab2Rgz6xFW5k4zm2Nm281ss5mNNbNBsYqxqjWqm8aNR3dm4qKNfDhTd3sVkZollkcQOcATwCBgKFAAfG5mmSFl5gLDgZ7AYGAx8LGZZcUwzip17oB2dGtRn5Ef/MyOPQVBhyMiErWYJQjn3DDn3Gjn3Czn3EzgAqApcEhImZecc2Odc4uccz8BNwD1gD6xirOqJScZd5zUndW5u3jiy4VBhyMiErWUALddDy9BbY4008zSgMuBrcD0Espc7pchKyuLcePGVSiQvLy8Ci8brYNbJPPUuAW0KVhJVt3qy8uxqEusJEpdEqEefbZsARKjLsVUlyg45wJ5Aa8B04DksOnHA3lAEbAS6B/N+vr27esq6ssvv6zwstFak7vTdfvLR+7S0ZOqdTuxqEusJEpdEqIehx/u3OGHJ0ZdfKqLB5jiSvheDeQqJjN7EK+P4TTnXPhj2L7EO6U0CPgYeM3MWsQ4xCqXVT+d3x+xH2PnrOOLOWuDDkdEpEwxTxBm9hBwDjDUObcofL5zbrtzboFz7jvn3G+BfOCyWMdZHS45pD0dmtblb2N+1uNJRSTuxTRBmNkjwLl4yWFOlIslAbWqL6rYSUtJYsQJ3VmycQfPfL046HBEREoVy3EQjwOX4B09bDaz5v4rw59f38xGmtkAM2trZn3N7DmgNV5/RUI4rHNTju6WxT+/WMCqLTuDDkdEpESxPIK4Gu/KpbHA6pDXjf78AqA78DYwHxgDNAYOc87NiGGc1e4vx3ejyDnu/nB20KGIiJQoZpe5OuesjPk7gFNiFE6g2mTW4crDO/LI2PmcO2ADgzo2CTokEZG96F5MAbkqpyOtG9Xmr+/+xJ4CPVhIROKPEkRA0lOT+dtJ3VmwLo9/fb3XxVwiIoFTggjQ0K5ZDOuexWNfzGf5ph1BhyMi8itKEAG7/YTuJJlx+3u6JbiIxBcliIC1bFib/zuyM1/MWccnP2mEtYjEDyWIOHDxIdl0bV6PO8b8xPbduiW4iMQHJYg4kJqcxMiTe7A6dxePjNUzrEUkPihBxIl+2ZmcfVAbnp2wmNmrtwYdjoiIEkQ8uemYrjSoncpt78yiqEgd1iISLCWIONKobho3H9uVqUs389qU5UGHIyL7OCWIOHNG39b0b5/J3R/OZv223UGHIyL7MCWIOGNm3H1KT3blF3HHmJ+CDkdE9mFKEHGoU7MMrhnaifdnrNbT50QkMEoQcerKwzvSOSuD296eRZ7GRohIAJQg4lRaShL3nNqL1Vt38Y9P5gYdjojsg5Qg4ljfdo24YGA7/j1xCdOWbQ46HBHZxyhBxLk/DutCVr10bnlrJvmFem6EiMSOEkScq5eeyp0n92DOmm2M+krPjRCR2FGCqAGO6pbFcT2b88jY+Sxanxd0OCKyj1CCqCFGnNCdWilJ3PzWTN2GQ0RiQgmihmhWP52//KYbkxZv4qXvlwYdjojsA5QgapAz+rXmsM5NufejOXpEqYhUOyWIGsTMuOfUniSZcdObM/SIUhGpVkoQNUyrhrW55biufLtwIy9PWhZ0OCKSwJQgaqBz+7dlUMfG3PPhHFZu2Rl0OCKSoJQgaiAz477TelHkHDfrVJOIVBMliBqqTWYdbj62K1/P36CHC4lItVCCqMHOH9COAe0zGfn+bFbn6lSTiFQtJYgaLCnJO9WUX1TELW/N1KkmEalSShA1XHaTutx0TFfGzV3Pq5N1qklEqk7MEoSZ3WJmk81sq5mtN7MxZtYjZH6qmd1nZjPMbLuZrTazl82sbaxirKkuOjibQR0bM/L9n1m3Q3d8FZGqEcsjiBzgCWAQMBQoAD43s0x/fh3gQOAu//0koA3wsZmlxDDOGicpybj/jN4kmfHMzN0U6l5NIlIFYpYgnHPDnHOjnXOznHMzgQuApsAh/vxc59xRzrn/OufmOucmAVcA+/svKUWrhrW5/cTuzNtcxLMTdFtwEam8IH+Z18NLUKU9Kq2+/x6xjJldDlwOkJWVxbhx4yoUSF5eXoWXjSeNnaNXpuPvH82hTu5SWter2V1MidIuiVCPPlu2AIlRl2KqS9mCTBCPANOBiZFmmlka8AAwxjm3IlIZ59woYBRAv379XE5OToUCGTduHBVdNt5s2/Mld0wq4OVFqbwz/BDSUmpukkiUdkmIejRsCEBGRkbNr4svIdrFV111iSpBmNlMoKwT28451zvK9T0IDAYGO+cKI8xPAV4CGgInRrNO8dSv5d3Q7/IXp/Lo2PncOKxL0CGJSA0V7RHEG6XMawpcCtSKZkVm9hBwNjDEObfXyXI/ObwC9ARynHMbo4xRfEd3b87pfVvzxLgFDN2/GQe2bRR0SCJSA0WVIJxzd4RPM7PawB/wOpvnAzeVtR4zewQvOeQ45+ZEmJ8KvAr08MusiSY+2dtfT+jGxIUb+cNrP/LB7wdTJ00XgolI+ZT7BLWZJfmdwwuA3wLXAn2ccx+XsdzjwCXAOcBmM2vuvzL8+SnA68BAv4wLKVO7vHHu6+qnp/KPM3qzZON2/jbm56DDEZEaqFwJwsxOBn4G7gEeAro6515w0d3j4Wq8K5fGAqtDXjf681vjjX1oCUwNK3NWeeIUz8EdG3Pl4R15dfJyPpq5OuhwRKSGibaT+hDgPuAA4DHgXufclvJsyDlnZcxfApRaRsrvhqM68+2CDdz81kx6t2lIy4Y6GBOR6ER7YvprYCfeJaXrgEvN9v4ud849WHWhSVVITU7ikbMP4LhHv+aG16bzn8sGkpykPCwiZYs2QSzDu8z15FLKOEAJIg5lN6nLiBO786c3ZvD0Vwu5OqdT0CGJSA0Q7VVM2dUch1SzM/q2Zvy89Tz46TwO6diE3m0aBh2SiMS5mjvMVsrFzLj75J40q1eL616dxvbdBUGHJCJxLqoEYWbHmtkSM2sQYV4Df97RVR+eVKUGdVJ56Kw+LNu0gxHv/RR0OCIS56I9grgGuN85lxs+w592H3BdVQYm1WNAh8YMH9KJ16eu4N3pK4MOR0TiWLQJohfweSnzvwCiug+TBO+6I/bjoOxG3PrWTBatzws6HBGJU9EmiKZAaY8qc0DjyocjsZCSnMSj5xxAWkoSw1+exq78ve6XKCISdYJYgXcUUZJegM5X1CAtGtTmwTP7MHv1Vu58X7fiEJG9RZsgPgDujHRPJDOrA/zNLyM1yJCuzbjisA785/tljPlxVdDhiEiciTZB3AU0AOab2U1mdpL/uhmY58+7u7qClOpz47AuHNi2Ibe8NZMlG7YHHY6IxJGoEoRzbh0wCJiBlwje9l93+dMOcc6tra4gpfqkJifx2LkHkpxkDH/5B/VHiMgvoh0H0QtY7pw7DmgCDMC7LXcT59xx/o32pIZq1bA2D5zRm59WbeWuD2YHHY6IxIloTzFNw0sMOOc2AyPwEsbmaopLYuzIbllcNrg9L363lPfUHyEiRJ8gwm//eRig+0YnmD8d05W+7Rpx85szmLd2W9DhiEjAdC8m+UVaShJPnHcgddJSuPLFqWzblR90SCISoGgThPNf4dMkwWTVT+ef5x7A0k07uPH1H4nuYYEikoiifR6EAS+Z2W7/73TgX2a2I7SQc+7EqgxOgjGwQ2NuObYrIz+YzdNfLeLKwzsGHZKIBCDaBPHvsL9fqupAJL78dnB7pi3bwt8/nkOvVg0Y1KlJ0CGJSIxF+8CgS6o7EIkvZsZ9p/di7tptXPvKNMZcO1jPsxbZx6iTWkqUUSuFp87vy678Qq76zw/sLtAgOpF9iRKElKpTswz+cUZvfly+hdvf/Umd1iL7ECUIKdOxPVswfEhHXp28nBcmLg06HBGJESUIicofjurCkfs342/v/8w3CzYEHY6IxIAShEQlKcl46Kw+dGxal+Ev/8DSjbrzq0iiU4KQqNVLT+VfF/YD4HcvTCFvd0HAEYlIdVKCkHJp17guj597IAvXb+f6V6dTVKROa5FEpQQh5XZIpyb85Tf78/nstTz42bygwxGRahLtSGqRX7loUDZz1mzjn18uoHPzepzYu2XQIYlIFdMRhFSImfG3k3rQPzuTG1//kSlLNgUdkohUsZglCDO7xcwmm9lWM1tvZmPMrEdYmVPN7BN/vjOznFjFJ+WXlpLE0xf0pVXD2vzuhSl6prVIgonlEUQO8ATes62HAgXA52aWGVKmLvAtcEMM45JKaFQ3jdEXHwTApc9PZvP2PQFHJCJVJWYJwjk3zDk32jk3yzk3E7gAaAocElLmRefcHcBHsYpLKi+7SV1GXdiPFZt3csVLU3XPJpEEYUHdW8fMWgCrgEOdcxPC5jUB1gNDnHPjSlnH5cDlAFlZWX1fffXVCsWSl5dHRkZGhZaNN0HW5btVBTw1YzcHt0zm8p61MAt/Um35JEq7JEI9+lx/PQATRo6s8XUplgjtUqwydRkyZMhU51y/SPOCvIrpEWA6MLGiK3DOjQJGAfTr18/l5ORUaD3jxo2josvGmyDrkgNkNJ/PPz6dx4Burbj+yM6VWl+itEtC1KNhQwAyMjJqfl18CdEuvuqqSyAJwsweBAYDg51zOh+RQIYP6cTSjTt4+PP5tGlUh9P6tg46JBGpoJgnCDN7CDgb7/TRolhvX6qXmXHXKT1ZnbuLm96cQWZGGkO6NAs6LBGpgJiOgzCzR4BzgaHOuTmx3LbETlpKEk+efyBdW9Tj6pd+YNqyzUGHJCIVEMtxEI8DlwDnAJvNrLn/yggpk2lmfYDi8RGdzKyPmTWPVZxSNeqlpzL64v40q1+LS5+fzIJ1eUGHJCLlFMsjiKuBesBYYHXI68aQMicC04Av/b//5f99ZezClKrStF4tXri0P8lJSVz03CTW5O4KOiQRKYdYjoOwEl4jQso8X1YZqVnaNa7L85ccRO7OfC56bhK5O/KDDklEoqR7MUm169GqAaMu6MviDdv57b8nsytfF66J1ARKEBITgzo14cGzejN12WaG/+cH9hQUBR2SiJRBCUJi5vheLbnzpB6MnbOO//vvdAoKlSRE4pmeByExdf7AduzcU8hdH86mVmoS/zi9N0lJlbslh4hUDyUIibnfHdaBHXsKeejzedROTWbkyT0qfd8mEal6ShASiN8f0Ykd+QU8PX4RtVOT+fNv9leSEIkzShASCDPj5mO6smtPIc9MWEydtGRuOLpL0GGJSAglCAmMmXH7Cd3ZmV/Io18soFZqMsOHdAo6LBHxKUFIoJKSjHtO7cXugiLu/2QugJKESJxQgpDAJScZD5zRGwPu/2QuBYWO647cL+iwRPZ5ShASF1KSk3jgzD4kJRkPfT6PwqIiDkgN5mmHIjXF7oJCnhq3iK5WPf9XNFBO4kZyknH/6b05s19rHv1iAW/OzyeoR+KKxLudewr53QtTeejzeczaUD23r1GCkLiSnGTce2ovzunfhvcX5XPvx3OUJETC5O7I58Lnvufr+eu577Se9M2qnpNBShASd5KSjLtO7snQNik8PX4RIz+YrSQh4lu7dRdnPj2R6cu38Ng5B3DWQW2rbVvqg5C4lJRkXNAtjbZtWvPshMVs25XP3af0JCVZv2lk37VofR4XPDuJLTv2MPri/gzer0m1bk8JQuKWN06iG/XTU3j0iwXk7sznkbMPID01OejQRGJuxootXDx6Mga8cvlAerVuWO3b1M8xiWtmxg1Hd+Gvx3fjk5/Wcunzk8nbXRB0WCIxNWH+Bs4Z9R21U5N5/cqDYwjQi14AABS0SURBVJIcQAlCaohLB7fnwTN78/3iTZz3r+/YtH1P0CGJxMT7M1ZxyfOTaJNZh7euHkSHphkx27YShNQYpx7YmqfP78vsNds48+mJrM7dGXRIItXGOcdT4xdyzcvT6NOmIf+94mCy6qfHNAYlCKlRjuyWxQuX9mdt7i5Of3Ii89ZuCzokkSqXX1jELW/N5N6P5nB8rxa8+NsBNKidGvM4lCCkxhnYoTGvXD6QPYVFnPbkt3yzYEPQIYlUma278rlk9GRenbyca4Z04tEAL8xQgpAaqUerBrwz/BBaNEjnoucm8fqU5UGHJFJpyzft4LQnvuW7RRv5++m9uHFYl0CfuKgEITVWq4a1eeOqQQzs0Jg/vjGDBz+dqwF1UmNNX76FU574hrVbd/HCpf05s1+boENSgpCarX56KqMvOYiz+rXh0S8WcMNrP7K7oHruSyNSXd6etoKznp5I7bRk3rp6EIM6Ve8AuGhpoJzUeKnJSdx7Wk/aNq7D/Z/MZeWWnTx1fl8y66YFHZpIqQoKi7jnozk8O2ExA9pn8sR5B9I4o1bQYf1CRxCSEMyM4UM68cjZfZi+fAsn/nMCP6/aGnRYIiXatH0PFz43iWcnLObiQdm8dNmAuEoOoAQhCeakPq147YqDKSh0nPrkN4z5cVXQIYns5adVuZzw2ASmLN3M/af3YsSJ3UmNw/uMxV9EIpXUp01D3rv2EHq0bMC1r0zj3o/mUFikzmuJD+/9uIrTnvyWwiLH61cczBlx0BldEiUISUjN6qXz8u8Gcu6Atjw1fiGXPj+Z3B35QYcl+7DdBYXc/u4sfv/KNHq2asCYawfTu01s7qlUUUoQkrDSUpK4+5Se3HVKD75duIGTHp/A7NXql5DYW7pxO6c9+S3/nriUywa35z+XDaRpvfjqb4gkZgnCzG4xs8lmttXM1pvZGDPrEVbGzGyEma0ys51mNs7MuscqRklM5w1oxyu/G8j2PYWc/Pg3vPz9Mo2XkJj5YMZqjn90Ass37eRfF/bjtuO7kZZSM36bxzLKHOAJYBAwFCgAPjezzJAyfwL+AFwLHASsAz4zs3oxjFMSUL/sTD667lD6t8/k1rdn8vtXp7Ntl045SfXZlV/IX96ZxfCXf6Bjsww++P1gjuqWFXRY5RKzcRDOuWGhf5vZBUAucAgwxswMuB641zn3pl/mIrwkcS7wdKxilcTUJKMW/76kP0+OX8gDn85l5oot/PPcA+nRqkHQoUmCWbAuj+tencZPq7byu0Pb88dhXWvMUUMoC+pQ28xaAKuAQ51zE8ysA7AQ6O+cmxxS7gNgg3PuogjruBy4HCArK6vvq6++WqFY8vLyyMiI3T3Wq5PqEp25mwp56sfdbMt3nNs1jSFtUvB+o1S9RGiTPtdfD8CEkSNrfF2KVUe7FDnH2GUFvDZ3D2nJcFnPWhzQrPp/h1emLkOGDJnqnOsXaV6QI6kfAaYDE/2/m/vva8PKrQVaRVqBc24UMAqgX79+Licnp0KBjBs3joouG29Ul+jkAGcM28MNr03nhZ/Xs5pG3HNqT5pUw0ClhGiTht7VNhkZGTW/Lr6qbpc1ubv44xs/8vX8DeR0acrfT+tFsxg9v6G6PmOBHPOY2YPAYOA051z4jXPCD2kswjSRSsusm8ZzFx3Ebb/Zn/Hz1nPMw1/x2c/hv09Eyvbu9JUc/dB4pizZzF2n9GD0xQfFLDlUp5gnCDN7CDgHGOqcWxQya43/3jxskWbsfVQhUiWSkozLDu3AmGsG07ReOr97YQo3vzlDz72WqGzevodrX5nGda9Op2OzDD687lDOG9Cu2k5XxlpME4SZPYLX4TzUOTcnbPZivCRxVEj5dOBQ4NuYBSn7pC7N6/Hu8EO4Kqcj/52ynGMf+UoPIpISOed4d/pKjnxwPB/NXM2NR3fm9SsOpn2TukGHVqVi1gdhZo8DFwAnA5vNrPhIIc85l+ecc2b2MPBnM5sDzANuA/KAl2MVp+y70lKSuOmYrgzt2ow/vTGD8575nrP6teHW3+wfyOMeJT6t2LyDP789i/Hz1tO7TUNeOrUn+7eoH3RY1SKWndRX++9jw6bfAYzw//13oDbwONAI+B442jmnBw9LzBzkj5l4+PP5/OvrRXw5dx13ntyDYd3Dz37KvqSwyPH8t0t44NO5ANx+QjcuPDib5ACf+FbdYjkOosy96Lxrbkfwv4QhEoj01GRuPrYrv+nZgj+9OYMrXpzKMd2b85cTutGqYe2gw5MYm7Uylz+/PZMfV+QypEtT7jy5B60b1Qk6rGqnBwaJlKJn6wa8d80hjPpqEY99MZ/xD6zn2iM6cdngDjVy4JOUz6bte/jHp3N5ZdIyGtdN49FzDuCEXi0SphO6LEoQImVITU5i+JBOnNi7JXe+/zN//3gub05dwZ0n9YibR0NK1SooLOLlSct44NN55O0u4JJB7bnuyP32ub4oJQiRKLXJrMOoC/vxxZy1jHjvZ8595nuO7pbFzcd2pUPTxBhdvK9zzjF+3nru/WgOc9ZsY1DHxow4sTuds/bN28EpQYiU09CuWQzq2IRnJyzmiS8XcPRDX3H+wHZcd8R+NNJzsGusWStzueej2XyzYCNtMmvz5HkHckyP5vvM6aRIlCBEKiA9NZnhQzpxZr82PPjZPF6YuIS3fljB8CGduPDgbGqnJQcdokRp2cYdPPDZXN6dvopGdVL56/HdOG9gW2qlqA2VIEQqoWm9Wtxzak8uHpTNPR/N5p6P5vDMhMUMz+nIOQP0JRPP1u8o4qY3ZvDmDytITjKuzunIlTkdqZ++b/UzlEYJQqQKdGlej+cv6c+kxZv4x6dzGTHmZ0Z9tYhrhu5HMz0PO64s37SDJ8Yt4LXJO0lOWsl5A9pyVU4nmjeo+fdOqmpKECJVqH/7TP57+UC+WbCRf3w6l1vfnklmujE8bTFnH9SGurX0Xy4oc9dsY9RXi3h3+kqSzMhpk8LI8w6jRQONaymJPq0iVczMGLxfEw7p1Jjx89ZzzztTufP9n3nsi/lcPCibiw7OVmd2jDjn+H7xJp4ev5Av566ndmoy5w9sx+WHdWDe9O+VHMqgBCFSTcyMnC7NYEBt6rXvzZPjFvLw5/N5avxCTu7TiosGZSfsPXyCtiu/kDE/ruLF75YyY0UujeumccNRnblgYLtfkvO8gGOsCZQgRGKgb7tGPHNRP+at3cbobxbz9rSVvDp5Of2zM7loUDZHd88iNVkjsytr6cbt/Of7Zbw2ZTlbduTTqVkGI0/uwel9W5OeqgsGyksJQiSGOmfV455Te3HTMV15bcpyXpi4lOEv/0CTjDROOaAVZ/Rrs88OyqqoHXsK+HjWGt78YQXfLtxIkhnDumdxwcBsBnbI3KfHMVSWEoRIABrWSePywzry28Ed+HLOOl6bspzR3yzhX18vpnebhpzRtzXH9mhO42p4BGoiKCxyTFq8ibd+WMGHM1ezfU8hbTJrc90R+3H2QW11RVIVUYIQCVByknFktyyO7JbFhrzdvDNtJa9PWcFt78zi9vd+YlDHxvymZwuGdW++z3dsFxQWMWnJJj6cuZqPZ61lQ95uMmql8JteLTjtwNYclJ1JUgLfejsIShAicaJJRi0uO7QDvx3cnp9Xb+XDmav5YMZqbn5rJn9+ZxYD2mcytGszcro0o2PTuvvEqZPcnflMmL+BcXPX8cWcdWzcvofaqckM7dqMY3s2Z2jXZtRJ09dYddGeFYkzZkb3lg3o3rIBNx7dhZ9Wecli7Ox1jPxgNiM/mE3bzDrkdGnKoI6N6d++MZkJcnSxp6CImSu38N2iTYyfu56pyzZTWOSon57CYZ2bclzPFuR0aaqkECPayyJxzMzo0aoBPVo14E/HdGXF5h2Mm7v+l36LFyYuBaBzVgYD2jemX3YjerVuSLvMOjXidMvm7XuYuTKXqUs3M2nxJqYt38yu/CIAuresz5WHdyCnSzMOaNOQFF3lFXNKECI1SOtGdTh/YDvOH9iOPQVFzFixhe8Xb+K7RRt584cVvPidlzDqpafQo2UDerVuQOesenRslkGHpnUDu8/QnoIilm3azoJ121mwbhszV+Yya+VWVm7ZCUCSwf4t6nNO/7YMaJ9Jv+xMmqiDPnBKECI1VFpKEv2yvS/T4UM6kV9YxLy125i1MpcZK3KZuTKX0d8sYU9h0S/LNKtXi/ZN6tKyYW1aNEj3X7VpUq8WDWqn0qB2KvXTU6L+te6cY9vuAnJ35JO7M5/NO/awOncXa3J3sTp3F6tzd7J04w6WbdpBYcg9qdo3qcsBbRty4cHt6NmqAT1aN9BN8uKQEoRIgkhNTvql7+Ksg7xp+YVFLNu0g4Xr8li4fjsL1+exZMN2Ji3exJqtu371pR0qo1YKaSlJpCYbqclJ/HP5FgD+OH4HKRPHkl9YxJ6CIvJ2F1DSvQibZKSRVT+d/VvU4/heLejQtC4dm2bQoWkGGbonVY2gVhJJYKnJSXRsmkHHCE+8KyxybMjbzaotO9m0fQ+5O/N/eW3dWUB+YZGXCAqLqFsrBQd0apREm5ZNSE1OIjU5iXrpKdRP9488aqfSqE4qLRrUpln9Whq5nACUIET2UclJRlb9dLLqRzGo7AkvwVzRK52cnN7VHJnEC10WICIiESlBiIhIREoQIiISkRKEiIhEpAQhIiIRKUGIiEhEShAiIhKREoSIiERkzpUwTr6GMbP1wNIKLt4E2FCF4QRJdYk/iVIPUF3iVWXq0s451zTSjIRJEJVhZlOcc/2CjqMqqC7xJ1HqAapLvKquuugUk4iIRKQEISIiESlBeEYFHUAVUl3iT6LUA1SXeFUtdVEfhIiIRKQjCBERiUgJQkREIlKCEBGRiBI6QZhZspndaWaLzWyX/z7SzEp8kp6ZZZuZi/A6JpaxlxBbPTN72MyWmtlOM/vWzA4qY5meZjbeL7/SzP5qZharmEuIqVz1iJc2MbPDzOw9fz86M7s4bL6Z2QgzW+XXa5yZdY9ivYeb2VT/M7rIzK6stkr8b5tVXhczyymhnboGXJdTzewTM1vvz8+Jcr3x2C7lrktl2iWhEwRwEzAc+D3QFbjO//uWKJY9BmgR8vqimmIsj2eAYcBFQE/gU+BzM2sVqbCZ1Qc+A9YCB+Hthz8CN8Qk2pKVqx4hgm6TDGAW3udoZ4T5fwL+AFyLt7/XAZ+ZWb2SVmhm7YEPgW+BA4B7gMfM7LSqDX0vVV6XEN35dTvNr4qAS1FWXeri7d+oP/dx3C7lrkuI8reLcy5hX8D7wL/Dpv0beL+UZbIBB/QLOv6wuGoDBcBJYdOnAiNLWOYqYCtQO2TabcBK/CvYakg94q5NgDzg4pC/DVgN/DmsrtuAK0pZz33A/LBpzwATa2Bdcvx2ahIv7RI2r4kfX04U64m7dqlEXSrcLol+BDEBGFJ8KGVm3YCheL8MyvKWma0zs2/M7PTqDDJKKUAysCts+k5gcAnLHAx87ZwL/SXyCdAS70s3CBWpR7F4a5NQ7YHmeEdDAPj7/StgUCnLHRy6jO8ToJ+ZpVZ1kFGqaF2KTTGz1WY21syGVFOM1S0e26Wyyt0uiZ4g7gNeBH42s3zgJ7wjiidKWSYPuBE4EzgOGAv818zOr+5gS+Oc2wZMBG4zs1bm9a+cj/dBblHCYs3xTi+FWhsyL+YqWI+4bJMwxfsz0v4ubV+X1EYpeL8Sg1DRuqzGO2o9DTgVmAuMNbPDqjzC6heP7VJRFW6XEjtrE8RZwIXAuXjJoQ/wiJktds49G2kB59wG4IGQSVPMrAneOdmXqjneslwAPAesAAqBH4BXgANLWSZ8JKSVMD2WylWPOG+TcJH2d1n7Oh7bKNL2S62Lc24u3pdPsYlmlo2X3L+q6uBiIF7bpVwq0y6JfgRxP/AP59yrzrmZzrkXgQeJrpM61PfAflUeXTk55xY65w7H68hq45zrD6QCi0tYZA17/+Jr5r+H/zqKmQrUI5K4aJMQa/z3SPu7tH1dUhsVABurJrRyq2hdIom3dopWPLZLVYqqXRI9QdTB+4UaqpDy17sP3mFaXHDObXfOrTazRnhXA71bQtGJwKFmlh4y7ShgFbCkeqMsWznqEUlctQlecluDt38B8Pf7oXhXnZRkInBk2LSjgCnOufyqDjJKFa1LJPHWTtGKx3apSlG1S6KfYhoD3Gxmi/FOMR2Ad3nYC8UFzOweoL9z7gj/74uAfGAaUAScgHdp7E2xDX1vZjYML7nNATrhHSHNBUb7839VF+Bl4HbgeTMbCXQGbgbucP7lDUEobz3ipU3MLMOPFz/+tmbWB9jknFtmZg8DfzazOcA8vCvG8vDaoXgdLwA45y70Jz0FXOMv+zRwCHAxcE5Nq4uZXY/3w+MnIA04HzgZ79x3kHXJBNoCDf0yncxsC7DGObcmUl2I33Ypd10q1S6xumQriBdQD3gY70lzO4FFwN1AekiZ54ElIX9fBPwMbMe7RHQKcH7QdfFjOxNYCOzGy/7/BBqUVBd/Wk+884y7/GVuJ6BLXCtaj3hpE/53uWD463l/vgEj/DrtAsYDPcLWMQ4YFzbtcLx+mN14v96vrIl1wesTWuD/X9sEfA0cFwd1ubiE+SNqYLuUuy6VaRfdzVVERCJK9D4IERGpICUIERGJSAlCREQiUoIQEZGIlCBERCQiJQgREYlICUJERCJSghARkYiUIEREJCIlCJEK8p/T/KSZPWBmm/znBF9nZrXM7HEz22Jmy8zsAr988bO1+4Wtx8XhA5BElCBEKuk8vEdxDgDuxbv31zt4N7jrh/eI22fMrGVgEYpUkBKESOX85Jwb4Zybj/eskQ1AvnPuEefcAuBveDe+i+ZRnSJxRQlCpHJmFP/DeXe+XAfMDJmWD2zmfw9qEqkxlCBEKif84TGuhGlJeM+ygP89uhIzS62+0EQqRwlCJHbW++8tQqb1CSIQkWgk+hPlROKGc26nmX0H3GRmC4EGwD0BhyVSIh1BiMTWpf77ZLxHWd4WYCwipdIT5UREJCIdQYiISERKECIiEpEShIiIRKQEISIiESlBiIhIREoQIiISkRKEiIhEpAQhIiIR/T/Ol1eQWtjx8wAAAABJRU5ErkJggg==\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.is_fixed('mu')):\n",
    "    plt.figure()\n",
    "    m.draw_mnprofile('mu', band=False, bound=(8.5, 11.5), bins=200)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}