From 7a56cdda45dff79e4bd6224462f577b0d0d599fa Mon Sep 17 00:00:00 2001
From: Mauro Donega <mauro.donega@cern.ch>
Date: Sat, 9 May 2020 00:34:48 +0200
Subject: [PATCH] binned likleihood expo fit
---
binnedLikelihood.ipynb | 1395 ++++++++++++++++++++++++++++++++++++++++
1 file changed, 1395 insertions(+)
create mode 100644 binnedLikelihood.ipynb
diff --git a/binnedLikelihood.ipynb b/binnedLikelihood.ipynb
new file mode 100644
index 0000000..33ebab1
--- /dev/null
+++ b/binnedLikelihood.ipynb
@@ -0,0 +1,1395 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Binned Likelihood\n",
+ "\n",
+ "\n",
+ "In this notebook we will be using probfit together with iminuit:\n",
+ "\n",
+ "probfit:\n",
+ "https://probfit.readthedocs.io/en/latest/\n",
+ "\n",
+ "iMinuit:\n",
+ "https://iminuit.readthedocs.io/en/latest/index.html#\n",
+ "\n",
+ "\n",
+ "1) fit an exponential and compare with the results from formula (average) \n",
+ "\n",
+ "2) same with gaussian pdf\n",
+ "\n",
+ "3) same with poisson\n",
+ "\n",
+ "4) Extended likelihood\n",
+ "\n",
+ "5) Binned likelihood\n",
+ "\n",
+ "6) model = flat background + gaussian \n",
+ " parameters = gaussian fraction / mu / sigma / const \n",
+ " generate toy data \n",
+ " fit different sets of parameters \n",
+ " show how uncertainties change with different parameters fixed \n",
+ " constrain parameters with gaussian constrains in the likelihood "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 1) fit an exponential"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "%matplotlib inline\n",
+ "import matplotlib.pyplot as plt\n",
+ "from math import exp, pi, sqrt\n",
+ "import probfit\n",
+ "from probfit import BinnedLH\n",
+ "from iminuit import Minuit, describe\n",
+ "from scipy.stats import norm, chi2, lognorm\n",
+ "import scipy.stats"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Generate data\n",
+ "# set the seed to always get the same samples\n",
+ "np.random.seed(seed=123456)\n",
+ "\n",
+ "# Generate a toy dataset on an exponential distribution (background)\n",
+ "# pdf = lambda * exp(-lambda * x) ; scale = 1/lambda\n",
+ "data = scipy.stats.expon.rvs(loc= 100, scale = 25, size=10000)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 40,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFBCAYAAADzMv2/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWfElEQVR4nO3df7RlZX3f8fdH5IcLdAw4yXKBkwGhGEITNBOqJlpC0zr+GDAWLSSrsS5WpiRFYtOEwoqxJKtpos0Pi0HIhBCIEQiiLQyOslwi0hpLAEUDDOiIGkatxHQ5/kiKgt/+cfaVw829Z/aee/fd58x5v9Y66+zznH32+c4e5sOzz7P3s1NVSJLaedLQBUjSLDE0JakDQ1OSOjA0JakDQ1OSOjA0JakDQ1OSOjA0JamDqQnNJD+Q5LIk1yf5+aHrkaSl9BqaSa5I8nCSexa1b07yQJJdSS4AqKqdVXUO8BpgU591SdK+Sp+XUSZ5MfAN4E+r6sSm7QDgU8A/B3YDdwBnVdV9SU4DLgD+oKqu3tv2n/GMZ9TGjRv7Kl/SnLrrrru+UlXrl3rvyX1+cVXdlmTjouaTgV1V9SBAkmuB04H7qupG4MYk7wWWDM0kW4GtABs2bODOO+/sqXpJ8yrJ55d7r9fQXMaRwENjr3cD/yTJKcCrgIOBHct9uKq2AdsANm3a5GwjktbUEKGZJdqqqm4Fbl3bUiSpmyFGz3cDzxp7fRTwxQHqkKTOhgjNO4Djkhyd5CDgTODGLhtIsiXJtj179vRSoCQtp+9Tjq4BPgocn2R3krOr6lHgXOBmYCdwXVXd22W7VbW9qrauW7du9YuWpAn6Hj0/a5n2HUwY7JGkaTU1VwRJ0iwwNCWpg5kMTQeCJA1lJkPTgSBJQ5nJ0FyJiy4augJJs2zuQlOSVsLQlKQOZjI0HQiSNJSZDE0HgiQNZSZDc6UuuuiJA0IODklqay5DU5L21RDzaU4Ne5iSuprJnqYDQZKGMpOh2cdA0OLfOSVpKTMZmpI0FENTkjowNCWpA0NTkjowNCWpg5kMTU85kjSUmQxNrz2XNJS5viJoKUudq+n5m5IWzGRPU5KGYmi24NVCkhYYmpLUgaEpSR0YmpLUgaEpSR3MZGh6crukocxkaHpyu6ShzGRoStJQDE1J6sDQ7MCT3CUZmpLUgaEpSR0YmvvIQ3VpPhma+8CwlOaXoSlJHcxkaHpFkKShzGRoekWQpKHMZGhOEweEpPliaEpSB4amJHVgaEpSB4bmKvL3TWn/Z2hKUgeG5iqxhynNB0NTkjowNHvgb5vS/svQ7JHBKe1/DE1J6sDQlKQODE1J6mAmQ9Op4SQNZSZDc5amhnMkXdq/zGRoStJQDE1J6sDQlKQODE1J6sDQXCMOCEn7B0NTkjowNCWpA0NzQB6yS7PH0FxjBqU02wxNSerA0JSkDgzNgXiILs0mQ1OSOjA0p4CDQ9LsMDQlqQNDc8rY65Smm6E5RQxLafoZmpLUgaEpSR0YmpLUgaEpSR0sG5pJnpbkt5K8I8lPL3rv7X0Uk+SVSf4oyQ1J/kUf3yFJKzGpp/knQIB3A2cmeXeSg5v3nt/2C5JckeThJPcsat+c5IEku5JcAFBV/6Oqfg74N8C/6vIHkaS1MCk0n11VFzRBdhrwMeCWJEd0/I4rgc3jDUkOAC4BXgqcAJyV5ISxVd7YvD+3PF9Tmk5PnvDewUmeVFXfAaiq30yyG7gNOKztF1TVbUk2Lmo+GdhVVQ8CJLkWOD3JTuC3gfdV1ceW2l6SrcBWgA0bNrQtQ5JWxaSe5nbg1PGGqroK+A/At1b4vUcCD4293t20vR74SeCMJOcs9cGq2lZVm6pq0/r161dYhiR1s2xPs6rOX6b9/cBxK/zeLL3puhi4eIXb3q+MH6J7uC4Nb6hTjnYDzxp7fRTwxYFqkaTWhgrNO4Djkhyd5CDgTODGth9OsiXJtj179vRWoCQtpffQTHIN8FHg+CS7k5xdVY8C5wI3AzuB66rq3rbbrKrtVbV13bp1/RQ9pRxRl4Y3afT8u5I8p6ruX3ju8gVVddYy7TuAHV22JUlDa9vTvHrRswZkb1MaTtfD86VGvSVpbszkhB0OBEkaykyG5rwOBEkaXtfQrF6qkKQZ0TY0s+hZA/P0I2kYbUPzRYueJWkutQrNqvrG+PPQHAiSNBQHgmach+nS2prJ0JSkobQKzSRPSXJ838VI0rTba2gm2QLcDby/eX1SktYzEknS/qRNT/MiRren+CpAVd0NbOyvpL1zIEjSUNqE5qNVNVXp5ECQpKG0mRrunua+5wckOQ44D/iLfsuSpOnUpqf5euAHgUcYTQ23B3hDn0WpO089ktZGm57mjwBvqqpfXWhI8jxG90GXpLnSpqd5M3BLku8ba7u8p3okaaq1Cc0HgP8K3JrkhU2bE3dImkttQrOq6ibgNOAPkpzLwFPEecrR8vxdU+pXm9AMQFV9mtEsRy8GfqjPovbGU44kDWWvA0FV9dyx5W8Cr0myodeqJGlKLRuaSc6vqrckuXiZVc7rqSatkoVDdQ/ZpdUzqae5s3m+ay0K0eoxJKX+LBuaVbW9eb5qoS3Jk4DDqupra1CbJE2dNrMcXZ3kaUkOBe4DHkjyK/2XJknTp83o+QlNz/KVwA5gA/Cve61KkqZUm9A8MMmBjELzhqr6Nt7KV9KcahOafwh8DjgUuC3J9wOD/qbpye3dLTWhh5N8SN3tNTSr6uKqOrKqXlZVBfw18BP9lzaxJk9u72A8GA1KaWXazHL0BE1wPtpDLZI09bwbpSR1YGhKUgdtztN8dZKnNstvTPKeZhJiSZo7bXqav1ZVX0/y48BLgKuAS/stS5KmU5vQfKx5fjlwaVXdABzUX0mSNL3ahOYXkvwh8BpgR5KDW35OkvY7bcLvNYzuE7S5qr4KHA547bmkudTm5Pa/A24AvtlMPnwgcH/fhU3iFUGShtJm9Pz1wJeBDwDvbR439VzXRF4RJGkobQ7PfxE4vqp+sKr+cfMY9B5B6oeXV0p71yY0HwI8Dp5jXq8uPa7NtecPMrrn+XuBRxYaq+r3eqtKg/G+QtJkbULzr5vHQXh+pqQ51+YWvr++FoVobdmTlPbNpFv4vrWq3pBkO0vM1F5Vp/VamdaMASq1N6mn+Y7m+XfWohBJmgWTbuF7V/P84SQHAc9h1ON8oKq+tUb1SdJU2etvmkleDlwGfAYIcHSSf1tV7+u7OEmaNm1Gz38X+Imq2gWQ5NmMrgoyNCXNnTYntz+8EJiNB4GHe6pHM8CT3TXPJo2ev6pZvDfJDuA6Rr9pvhq4Yw1qk6SpM6mnuaV5HMJowo5/CpwC/A3wPb1Xpqlnb1PzaNLo+evWshBJmgWd73s+DZJsAbYce+yxQ5cyV+xZSjN62wrn01w7BqX0RMuGZpIXJMlaFiNJ027S4flrgUuSfAp4P/D+qvo/a1OWhmYPU1rapIGgcwCSPAd4KXBlknXAhxiF6Eeq6rHlPi9J+6M2N1a7v6p+v6o2A6cC/4vRuZq3912cJE2bTqPnVfX3wI7mIUlzZyZHzyVpKIamVp3Xpmt/1ua+54cmeVKz/I+SnJbkwP5Lk6Tp06aneRtwSJIjgQ8CrwOu7LMoSZpWbUIzVfV3wKuAt1XVTwEn9FuWJE2nVqGZ5AXAzzCafBhm9Jp1SVqpNqH5BuBC4L9X1b1JjmF0grv0DwZ9ugwAOWCkWdTmvucfBj6c5NDm9YPAeX0XJknTqM3o+QuS3AfsbF7/cJK3916ZJE2hNofnbwVeAvwtQFV9Anhxn0VJ0rRqdXJ7VT20qMmJOiTNpTaj4A8leSFQSQ5i9Hvmzn7LkqTp1KaneQ7w74Ajgd3ASc1rSZo7bUbPv8LoHE1JmnuT7nt+flW9JcnbGN3v/AmqytOOJM2dSYfnC79b3gnctcRjVSU5JskfJ7l+tbetYXkCu/Ynk253sT3JAcCJVfUr+7LxJFcArwAerqoTx9o3A/8NOAC4vKp+uzlp/mxDU9I0mzgQ1NwD6EdWsP0rgc3jDU0QX8LovkMnAGclcQIQSTOhzSlHH09yI/Au4JsLjVX1nr19sKpuS7JxUfPJwK6mZ0mSa4HTgfvaFJxkK7AVYMOGDW0+oimycKjuIbtmVZtTjg5ndDXQqcCW5vGKFXznkcD4yfK7gSOTHJHkMuC5SS5c7sNVta2qNlXVpvXr16+gDEnqrk1P8/Kq+sh4Q5IfW8F3Zom2qqq/ZXROqCRNrTY9zbe1bGtrN/CssddHAV9cwfYkac1MOk/zBcALgfVJfmnsracxGvXeV3cAxyU5GvgCcCbw0102kGQLsOXYY49dQRnq277OsylNs0k9zYOAwxgF61PHHl8Dzmiz8STXAB8Fjk+yO8nZVfUocC5wM6NzQa+rqnu7FF1V26tq67p167p8TJJWbNJ5mguTD19ZVZ/fl41X1VnLtO8AduzLNiVpSG0Ggg5Osg3YOL5+VZ3aV1GSNK3ahOa7gMuAy3EeTUlzrk1oPlpVl/ZeSQcOBE0fB3o0L9qccrQ9yS8keWaSwxcevVc2gQNBkobSpqf52uZ5fNKOAo5Z/XIkabq1mYT46LUoRJJmwbKH50nOH1t+9aL3/kufRUnStJr0m+aZY8uLJ9DYzICSbEmybc+ePUOWoZ4svpLIQSZNk0mhmWWWl3q9phwIkjSUSaFZyywv9VqS5sKkgaAfTvI1Rr3KpzTLNK8P6b0ySZpCk649X8lMRpK0X2pzcrskqTGToeno+XybNKLuSLv6NpOh6ei5pKHMZGhK0lAMTUnqwNCUpA4MTUnqYCZD09FzSUOZydB09FzSUGYyNCVpKIamJHVgaEpSB4amJHVgaEpSB4amJHVgaEpSB23uez51kmwBthx77LFDl6JVNGlaN6d807SYyZ6mJ7dLGspMhqYkDcXQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQOvCJIg1vJ1T4Ln+26jfH121yJ5BVJWjCTPU2vCJI0lJkMTUkaiqEpSR0YmpLUgaEpSR0YmpLUgaEpSR0YmpLUgaEpSR0YmpLUgaEpSR0YmpLUgaEpSR0YmpLUgaEpSR04n6YGsdrzU7adE3Op9Ra3LbetSXNrtpl3c7nv12yZyZ6m82lKGspMhqYkDcXQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQODE1J6sDQlKQOpuZulEkOBd4OfAu4tareOXBJkvQP9NrTTHJFkoeT3LOofXOSB5LsSnJB0/wq4Pqq+jngtD7rkqR91ffh+ZXA5vGGJAcAlwAvBU4AzkpyAnAU8FCz2mM91yVJ+6TXw/Oqui3JxkXNJwO7qupBgCTXAqcDuxkF591MCPMkW4GtABs2bFj9ojUzLrpoddqXWq/Nugvr7O379ratNuu13dY02x/+DDDMQNCRPN6jhFFYHgm8B/iXSS4Fti/34araVlWbqmrT+vXr+61UkhYZYiAoS7RVVX0TeN1aFyNJXQzR09wNPGvs9VHAFweoQ5I6GyI07wCOS3J0koOAM4Ebu2wgyZYk2/bs2dNLgZK0nL5POboG+ChwfJLdSc6uqkeBc4GbgZ3AdVV1b5ftVtX2qtq6bt261S9akiboe/T8rGXadwA7+vxuSeqDl1FKUgeGpiR1MJOh6UCQpKHMZGg6ECRpKDMZmpI0FENTkjpIVQ1dwz5L8jfA5zt+7BnAV3ooZ7VZ5+qblVqtc3XtS53fX1VLTm4x06G5L5LcWVWbhq5jb6xz9c1Krda5ula7Tg/PJakDQ1OSOpjH0Nw2dAEtWefqm5VarXN1rWqdc/ebpiStxDz2NCVpnxmaktTBfhWaS90yOMnhST6Q5NPN8/c07UlycXMb4U8med4U1HpRki8kubt5vGzsvQubWh9I8pI1rPNZST6UZGeSe5P8YtM+Vft1Qp1TtU+THJLkL5N8oqnz15v2o5Pc3uzPP28m6CbJwc3rXc37Gweu88oknx3bnyc17UP/ezogyceT3NS87m9/VtV+8wBeDDwPuGes7S3ABc3yBcCbm+WXAe9jdM+i5wO3T0GtFwG/vMS6JwCfAA4GjgY+AxywRnU+E3hes/xU4FNNPVO1XyfUOVX7tNkvhzXLBwK3N/vpOuDMpv0y4Oeb5V8ALmuWzwT+fI3253J1XgmcscT6Q/97+iXgauCm5nVv+3O/6mlW1W3A/13UfDpwVbN8FfDKsfY/rZH/DTw9yTPXptJla13O6cC1VfVIVX0W2MXoVsi9q6ovVdXHmuWvM5pt/0imbL9OqHM5g+zTZr98o3l5YPMo4FTg+qZ98f5c2M/XA/8syVI3J1yrOpcz2L+nJEcBLwcub16HHvfnfhWay/i+qvoSjP5hAd/btC93K+Ghndsc3lyxcMjLlNTaHMo8l1GvY2r366I6Ycr2aXMoeTfwMPABRr3cr9boVjCLa/lunc37e4Ajhqizqhb25282+/P3kxy8uM7GWv69vxU4H/hO8/oIetyf8xCay1nyVsJrXsUTXQo8GzgJ+BLwu0374LUmOQx4N/CGqvrapFWXaFuzWpeoc+r2aVU9VlUnMboT68nAD0yoZWrqTHIicCHwHOBHgcOB/zhknUleATxcVXeNN0+oZcV1zkNofnnhMKF5frhpn7pbCVfVl5v/UL8D/BGPHy4OWmuSAxkF0Tur6j1N89Tt16XqnNZ92tT2VeBWRr8BPj3Jwj27xmv5bp3N++to/7POate5ufkZpKrqEeBPGH5//hhwWpLPAdcyOix/Kz3uz3kIzRuB1zbLrwVuGGv/2WbU7/nAnoXDzaEs+g3op4CFkfUbgTObkb+jgeOAv1yjmgL8MbCzqn5v7K2p2q/L1Tlt+zTJ+iRPb5afAvwko99fPwSc0ay2eH8u7OczgFuqGcUYoM77x/5HGUa/E47vzzX/e6+qC6vqqKrayGhg55aq+hn63J9rOcLV9wO4htEh2LcZ/R/lbEa/V3wQ+HTzfHg9Pjp4CaPfk/4K2DQFtb6jqeWTzV/uM8fW/9Wm1geAl65hnT/O6PDlk8DdzeNl07ZfJ9Q5VfsU+CHg40099wBvatqPYRTau4B3AQc37Yc0r3c17x8zcJ23NPvzHuDPeHyEfdB/T00Np/D46Hlv+9PLKCWpg3k4PJekVWNoSlIHhqYkdWBoSlIHhqYkdWBoamYl2Zjk75tL/fZ1G5uSXNwsn5LkhXtZ/0VJ7svY7FSaL4amZt1nanSp3z6pqjur6rzm5SnAxNCsqv/J6PxPzSlDU1MpyY82k0IckuTQZk7HE/fymY154vykv5zkomb51iRvbuaI/FSSFzXtpyS5qZnk4xzg3zfzRL4oyauT3JPRnJK39faH1Ux58t5XkdZeVd2R5EbgPwNPAf6sqlZ6SPzkqjo5o4mI/xOjSwMXvu9zSS4DvlFVvwOQ5K+Al1TVFxYuKZQMTU2z3wDuAP4fcN5e1m1jYbKRu4CNLdb/CHBlkuvGPqs55+G5ptnhwGGMZmI/pMX6j/LE/6YXf+aR5vkxWnQYquoc4I2MZsW5O8mazGOp6WZoapptA34NeCfw5hbrfxn43iRHNJPjvqLj932dUUADkOTZVXV7Vb0J+ApPnPpMc8rDc02lJD8LPFpVVyc5APiLJKdW1S3Lfaaqvp3kNxjN2P5Z4P6OX7sduD7J6cDrGQ0KHcdoBp8PMrqnkOacsxxpZjUj3jdV1cRR9f3lezUdPDzXLHsMWLeSk9u7ak5V2s7ocF1zyJ6mJHVgT1OSOjA0JakDQ1OSOjA0JakDQ1OSOvj/o4kEMf2/CQsAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<Figure size 360x360 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.figure(figsize=[5,5])\n",
+ "plt.subplot(111)\n",
+ "plt.hist(data, bins=150, range=[100,400], color='blue', alpha=0.5)\n",
+ "plt.xlabel(r'x [units]')\n",
+ "plt.ylabel(r'Entries / bins size = 2')\n",
+ "plt.yscale('log', nonposy='clip')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 41,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "['loc', 'scale']"
+ ]
+ },
+ "execution_count": 41,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "def exp_func(x, loc, scale):\n",
+ " return scipy.stats.expon.pdf(x, loc, scale)\n",
+ "\n",
+ "blh = BinnedLH(exp_func, data, bins=150, bound=(100,400))\n",
+ "describe(blh)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 42,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAE1CAYAAACflJmiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxUVZrw8d+ppEIlkEBYBCTI0kEWEZGETQUDsolbO0O3ONrdLrhN2zrTo6O+bfc79rz9mVftmXntblulXbC7GUFtbTd2hUZkEXEDZZXNBAQJgexkO+8ft27lplKpVFLLvbfq+X4++Vh1c+vepy7ycM+5zzlHaa0RQggBHrsDEEIIp5CEKIQQfpIQhRDCTxKiEEL4SUIUQgg/SYhCCOEnCVEIIfwkIQohhF/CEqJSaqRS6mml1KtKqbsSdV4hhIhUVAlRKfW8Uuq4UmpH0PY5SqndSql9SqkHAbTWO7XWdwLfBwqjOa8QQsRDtHeIi4A51g1KqTTgSeByYBRwvVJqlP93VwMbgHejPK8QQsRcVAlRa70eOBm0eQKwT2u9X2tdBywBrvHv/6bW+iLghmjOK4QQ8ZAeh2MOAL62vC8GJiqlioC/A7oAy9r6sFLqduB2gMzMzIKBAwdGfOKmpiY8nuYcn15VRWZJCU1dulA1aFBHvkNcBcfpZG6J1S1xgntidUuc0PFY9+zZc0Jr3afVL7TWUf0Ag4EdlvffA561vP8B8NvOHLugoEB3xNq1a1tuqKnROjNTa9C6pKRDx4qnVnE6mFtidUucWrsnVrfEqXXHYwU+0iFyTjzSfzFgva3LA47E4Tzt8/ngssuM1ytW2BKCEMI94pEQtwLDlFJDlFIZwHzgzY4cQCl1lVJq4enTp6OP5vLLjf8ua7OVLoQQQPRlNy8Bm4DhSqlipdStWusG4G5gJbATeFlr/UVHjqu1fktrfXv37t2jCc9gJsTVq6G+PvrjCSGSVlQPVbTW17exfRlhHpwk1JAhMHIk7NwJGzfCpZfaHZEQwqHc8QgpWtJsFkJEwJEJsbN9iNV1jTy5dh/bDpW1/MXcucZ/ly+PUYRCiGTkyITYmT7EbYfK2P9tJY+v3M0Nz25umRQvuQS6dYPt2+Hrr9s+iBAipTkyIXbG5v2lmOsH1jc0sXl/afMvu3SBGTOM13KXKIRoQ9IkxElDe+FRCg+glCI3K6PlDtKPKIRoRzyG7kVNKXUVcFV+fn7EnykYlMuhHT409TQ0aX759hcM75fNz/+6nfLaBn4/dTJjAN59F86cMe4ahRDCwpF3iJ2tQ2xs0oFmc119E7986wt2Hq2guKyG7y8rpnrEKKishA0bYh+0EML1HJkQO6trl3R8Xg8eoAn4rPh0IEGeqW9i64iJxhtpNgshQkiqhJiVkcbiBZO4eFhvDv/XvBa/08DCbiOMNwl+sPK73/2O/Px8lFKcOHGiOSatueeee8jPz2fMmDF8/PHHgd+9+OKLDBs2jGHDhvHiiy/GNb5bbrmFs846i9GjR8f1PEI4nSMTYjRjmQsG5fJPM84FIE0ZP6aP+g/nTNduxqiVAwdiFW67Lr74YtasWcOgoCnIli9fzt69e9m7dy8LFy7krruMlRVOnjzJI488wpYtW/jwww955JFHKCtrWVt58OBBioqKYhLfTTfdxAqZ/EIIZybEaMcyFwzKxedN46ezhvPLa0Zzau3zHHnuHzm46F7+PGy4sdPy5Tz22GOcf/75XHDBBTz44IMx/AYtXXjhhQwePLjV9jfeeIMf/vCHKKWYNGkSp06d4ujRo6xcuZKZM2fSs2dPcnNzmTlzZlwT1tSpU+nZs2fcji+EWzjyKXMseBT8eFo+f/nLXxieUcqNL61ieA9YMPsi5gKfLlrEX9PT2bJlC1lZWZw8GTzxNyxevJjHH3+81fb8/HxeffXVqGMsKSnBOgFuXl4eJSUlbW4HuPbaazlw4AB1dXUcPnyYsWPHAnDvvfdy8803Rx2TEKksaROiacOGDdxx8w+5ZYZxZ3hpURFbX3uNv33yCTc/8QRZWVkAIe+QbrjhBm64IX6rHRjzVLaklGpzO8Drr78OGE3mm266iXXr1sUtPiFSjSObzLHUKrlkZcHQoeiGBtSePWE/u3jxYsaOHdvqZ968eWE/F6m8vDy+tgwlLC4u5uyzz25zuxAivhyZEGM5QezUqVNZunQpjY2NrPl4L++seo+cwqnMAh577k98sMtoioZqMt9www18+umnrX5i0VwGuPrqq/njH/+I1prNmzfTvXt3+vfvz+zZs1m1ahVlZWWUlZWxatUqZs+e3eKzgwcPlrtDIWLMkQkxlhPEXnvttYwZM4ZzR41m7pyZpE26kd96RzEH+G5jPdMumUzfISO47+f/J/rA2/Cb3/yGvLw8iouLGTNmTKBfcu7cuQwdOpT8/Hxuu+02fv/73wNG8/3nP/8548ePZ/z48fziF78INOmvvfbakHetL7zwQqfju/7665k8eTK7d+8mLy+P5557LvovLYQLJW0fYmVlJWD0vT3++OMMnnsHj6/cDcCnTY2c8nXjsZoKPrzx1xzsOYCNXg/bDpVRMCg35rHcc8893HPPPYH35p2dUoonn3wy5GduueUWbrnlllbbzT7EWHrppZdifkwh3MiRd4jxYEz+YLxu9KSxfsg4AIr2bwNCzJAjhEg5SXuHGKxgUC6v3HkRm/eXkpuVwZnqWbBzPdP2f8SiwqvxpnuYNLSX3WEKIWyUMneIYCTFH0/L5x8mnsP3Hr4NgCklX/DQpeeweMGkFs3lgwcPOmIo25w5c+jRowdXXnlli+0HDhxg4sSJDBs2jOuuu466urq4xfDKK69w3nnn4fF4+Oijj+J2HiHsllIJsYWzzoLx4/HUneGOpsNx6TuMhfvvv58//elPrbY/8MAD/PM//zN79+4lNzc35IOQoqIiDh48GHUMo0eP5rXXXmPq1KlRH0sIJ3NkQoy27KaqqoorrriCCy64gNGjR7N06VIAtm7dykUXXcQFF1zAhAkTqJg+nYPAlNtuY9y4cYwbN46NGze2Ol5jYyP3338/48ePZ8yYMTzzzDNRfLuOueyyy8jOzm6xTWvNe++9F6iH/NGPfsRf//rXuMUwcuRIhg8fHrfjC+EUjuxD1Fq/BbxVWFh4W2c+v2LFCs4++2zeeecdAE6fPk1dXR3XXXcdS5cuZfz48ZSXl5O5fTtnPfooq71evnhtDW++/zEL7vwxX37+SYvjPffcc3Tv3p2tW7dy5swZLr74YmbNmsWQIUNa7DdlyhQqKipaxfPrX/+aGeYSBjFQWlpKjx49SE83/visQ/teeOEFnnjiCQD27dvH3LlzycjIYMiQIXF5Qi1EMnFkQozW+eefz3333ccDDzzAlVdeyZQpU9i+fTv9+/dn/PjxAOTk5MCkSVT17Mndhw/z+rjzqfN2oeHkEbYdKsP6eGXVqlV8/vnngYLs06dPs3fv3lYJ8f3330/I9ws3tO/mm28OjGkuKipi0aJFISeWEEK0lpQJ8dxzz2Xbtm0sW7aMhx56iFmzZvHd7343kDQC0tL47wED6HvyJPdfMIfnx1/D4V9fy+b9pVwxpPnSaK357W9/22q0SLBE3SH27t2bU6dO0dDQQHp6ugztEyJGkjIhHjlyhJ49e3LjjTfSrVs3Fi1axIMPPsiRI0fYunUr48ePp6KigszMTE4PGEDe9u1MOrCN32Zmg24iNyuDP27aT019IwCzZ8/mqaeeYvr06Xi9Xvbs2cOAAQPo2rVri/Mm6g5RKcW0adN49dVXmT9/Pi+++CLXXHNNq/1kaJ8QHZOUCXH79u3cf//9eDwevF4vTz31FBkZGSxdupSf/OQn1NTUkJmZyZo1a/jHRx7h71es4JXDnzN+QgGrMnz8r9e303D6GN+W1bDtUBkLFizg4MGDjBs3Dq01ffr0ietDDKspU6awa9cuKisrA8PqZs+ezaOPPsr8+fN5+OGHufDCC7n11luBln2IVtH0Ib7++uv85Cc/4dtvv+WKK65g7NixrFy5MqrvJYQTqVD9UU5RWFioO1L3tm7dus7NIn3RRbBpE/sW/pmZX/UIrMPiUfAvs4bz42mRr/4XiU7HaQO3xOqWOME9sbolTuh4rEqpbVrrwuDtjiy7STj/ms1n3nwL6z8PHqVk9IoQKUQSIsDcuQAM+3gDvnRjsft0j6J3twzuXfIJ2w6Vhf+8ECIpODIhxnI+xIhceCH07UvGkRL+MrU7/zLbWIvleMUZistquOHZzZIUhUgBjkyIsZwPMSIeT6DZfN6nH/DjafmUVdfR5G8/y0w4QqQGRyZEW/gTormI/aShvfB5PaQp8KZ7eOWjr7nk0ffkTlGIJCYJ0TRzJqSlwQcfwOnTFAzKZfGCSfx01nBumjyYQ6XV0nwWIslJQjTl5hrlNw0NsGYNYEwXNmloL55Zvz/w9LlOms9CJC1JiFZBzWaAzftLW5TiKODIqRq5SxQiCUlCtPKX37B8OfgL1q1LD6QpaNKweMth5j21kQm/Wi2JUYgkIgnRaswYOPtsOHoUPvsMaF564P7Zw5k/4ZzA3aIGjlfUSZ+iEElEEqKVUs13iZZms7n0wN+NywvcLZqkJEeI5CEJMViIfkRTwaBcRvTLpk92BmZelMWphEgeSTnbTVRmzID0dNi0CcrKjKfPFsvuNdYV2XaojDv//BEq1DGEEK7kyDvEhA/ds8rJgSlToKkJVq0Ku2tpZZ30IwqRRByZEBM+dC9YmGazafP+UhnaJ0SScWRCtJ35YGXFCuNOMQRrOY70IwqRHCQhhjJqFJxzDhw/Dh9/HHIXazlOvxyfTBMmRBKQhBhKG+U3wcyhfYdPyjhnIZKBJMS2RNCPCG33Jc59Yr3MjiOEy0hCbMv06ZCRAR9+CCdOtLmbtS+xScOLGw/wP1sOs+ubCrlrFMJlJCG2pVs3uPRSY0xzmBXmzL7Ef5h4DmAM5/vFGzvkCbQQLiQJMZwIm80Fg3IZ0CMzMM65oal5fhx5Ai2Ee0hCDMd8sLJyJTQ2ht3V2nQ2KaBfjo/HVuyKT3xCiJiShBjOuefC0KFQWgpbt4bd1RznnONrHg3pUVBV10CJzJ8ohCtIQgwnwvIb07J7p/LCzROa50/0KE5U1AUersx9Yj3XPbMpjgELIaIhCbE9EfYjmsw7xbzcTL5XODDQr1hb38Se45WUnKqhui5881sIYQ9JiO0pKgKfD7Ztg2PHIvrIsnunsuGB6a3mT2xo1BSX1XDgRJU0oYVwIEmI7cnKgmnTjNcrVnToo9b5E6201lKKI4QDSUKMRAebzVbL7p3KTRcNCdqqpBRHCAeShBgJMyGuWmUsU9pB1pKcdI/i7B4+Cgblhv+QECLhEpYQlVLfVUr9QSn1hlJqVqLOGxP5+UYJzqlTsHlzhz9unRln6R2TOVVTL0+bhXCgqBKiUup5pdRxpdSOoO1zlFK7lVL7lFIPAmit/6q1vg24CbgumvPaogPlN6GYC1UVDMqlqUlLbaIQDhTtHeIiYI51g1IqDXgSuBwYBVyvlBpl2eVh/+/dJYp+RKtth8qorW+U2kQhHCiqhKi1Xg+cDNo8Adintd6vta4DlgDXKMOjwHKtdehZV51s6lTjifNnn0FJSacPs3l/aYvaxF3fVHDgRKXcLQrhAEpr3f5e4Q6g1GDgba31aP/7ecAcrfUC//sfABOBPcCPgK3Ap1rrp9s43u3A7QB9+/YtWLJkScSxVFZW0q1bt05/l/aM/tnP6L1xI7vuu49vrriiU8eormukqrKSb2pabldKMbR3V7Iy0mIQaezE+5rGilviBPfE6pY4oeOxTps2bZvWujB4ezyWIQ21MqfWWv8G+E17H9ZaLwQWAhQWFuqioqKIT7xu3To6sn+H3XgjbNzIiP37GRHFeZ59+W3+sNdLeW3LJ9Z5uYoND3T+uPEQ92saI26JE9wTq1vihNjFGo+nzMXAQMv7POBIHM6TeGY/4urVUF/f6cPkn9WNF26e0OpfjrqGRulTFMJG8UiIW4FhSqkhSqkMYD7wZkcOYOu6zOEMHmwsQFVRAR98ENWhCgblMrJ/NulpitysdBTG5LI7j0qfohB2ibbs5iVgEzBcKVWslLpVa90A3A2sBHYCL2utv+jIcW1flzmcKMtvrLJ9XgrOyWXBlO8EHrRojMQoSw8IkXjRPmW+XmvdX2vt1Vrnaa2f829fprU+V2v9Ha31r2ITqkPEqPwGYOkdk1l6x+SQk8vK0gNCJJ4M3euoSy4x1lv54gs4fDgmh7ROAmHmRQ2883lydL0K4RaOTIiO7UMEYyW+mTON18uXx+ywy+6dytafzWRkf2MuxXN6ZlFe2yDNZiESyJEJ0dF9iBDTZnOwbJ+XHF86h09WU1xWw/ee3ihJUYgEcWRCdDwzIb77Lpw5E9NDL71jMleMOTuwjGmThr98XBzTcwghQnNkQnR0kxkgLw/GjIGqKnj//ZgfftLQXi1qFF/dVix3iUIkgCMTouObzBDT8ptgBYNyW8yy3dgoT5yFSARHJkRXiGM/IsBTNxYGSnFksXshEkMSYmdNngzdu8Pu3bB/f8wPb129b2jvrrLYvRAJIAmxs7xemOWf+DuG5TdW2T4vA3pkku3zxuX4QoiWHJkQHf9QxRTnZrM5kqWitl5m2BYiARyZEF3xUAVgjn+y8LVroaYm/L6dtO1QGbu+qQjMsC1JUYj4cWRCdI3+/WHcOCMZ/u1vcTnF5v2lgZrE2vom7l3ySVzOI4SQhBi9OJbfQMslTBXGnIlylyhEfEhCjFac+xGtEz+AMTXY957eKBPJChEHjkyIrnmoAjBxIvTsCV99BXv3xuUU2T4vXdLTAnMmNmlaLT8ghIieIxOiax6qAKSlwezZxus4Pm1+Yv6FgaazR0G6R8mTZyFizJEJ0XXi3GyG5qZzepqie2Y6h0qr5cmzEDEWj1X3Us/s2aCU8aS5qgq6do3bqRoaNWXVzc3l2vomrv/DZi4c2COwbekdkwECfYzmeyFEeHKHGAtnnQXjxxtTga1dG7fThOo3NJvPQojoSUKMlTiX3wAt+hEVcFZ2BiP6ZQNQcqqGitqWS6PKCBchOkYSYqxY+xG1Dr9vJ1knfBjZP5shvbtRXddIdV0jxWU1LZYwlREuQnScIxOiq8puTIWF0KcPHDoEu+I3M411woeK2noOlVYHfmddwvS1j4sDI1xkBT8hIuPIhOiqshuTx9M8tjmOzWZzwgdzqYFQ96L1DU1okPkUheggRyZE10pA+Y1V8HrO5ktvuoe/H5cXaF4vXjCJgkG5CYlJCDeTsptYmjXLuFN8/32oqIDs7LiezuxT3PdtFfl9jFKf8toGnph/IQWDcsn2ecn2eSUZChEhuUOMpV69YNIkqK83VuRLAGNYnyeQ/Ab0yAyZAK97ZpOMfRaiHZIQYy0B5TfBRvXPCVl8bfY1tkeSpRAGaTLH2uWXw8MPN5ffqPgWTVsTXqTJL9J9hUg1cocYa2PHQr9+UFICO3bYEkKoOz6zSDu4eFsI0cyRCdGVdYgmjyfhT5vDue6ZTcx9Yn2gSHvXNxWSFIVogyMToivrEK0clBDBePJsFmmHmktRhvgJYXBkQnS9mTONeRI/+AAccJeb40tvcxkCGeInRDNJiPHQowdcfDE0NsLq1Qk/ffAdX7bP22oZAjP5WRexkiF+ItVJQowXG8pvoPUd37HyWkpOGUukWpchMJOfdbSLN93DO58fkRIckbIkIcaL2Y+4fDk0NSXstNY7vrr6psDM2ru+qSDdo1qNb7bOoLN4wSSyfV5AahNFapKEGC/nnw8DBsA338BnnyXstNY7Po9HtViYqqFJhxzfnO3zcrqmnsdWxG+WHiHcQBJivChlS7PZesf3y2tGt1iYKseXHnZ4nxCpThJiPNlUfmMmvX+YeE4gOb5y50Usu3dqi/2kgFuIliQhxtNll4HXC5s3w8mTtoQQ6o7wy6PlLRLh0jsmM6p/DhW19YEHMubs29V1jXaELYQtJCHGU04OTJliPFRZtcruaIDm5Besorae4rKawAMZc/btAyeqpDZRpAxJiPFmQz9ipLPcmMxSnVCr+mmtpTZRpAxHJkRXj2UOZvYjrliR0PIbU7jkaPYXWtdfAcjKSAvMvq2UarX8gJTkiGTlyITo+rHMViNHwqBB8O23sG2b3dEEWPsLX/no60AC9CgY3CuLkf2NhzFDeneVJ9IiZTgyISYVm8pvwgnuL2xs0vTOziAvN5MR/bJbzL6dlZFmb7BCJJAkxERw0Ow3ofoLveke+nTrEljeVIhUJQkxEaZPh4wM2LrVaDrbyDq0D4xibeuQvWD7T1QF+gs72ncofY3CbSQhJkLXrlBUZCwpsHKlraG0GNqnIC+3ZY1iJE+oZf5EkawkISaKQ/oRrUP7zP7CUNpKjNaHMTJ/okg2khATxexHXLnSmCfRRuYDk870F1pn35b5E0WykVX3EmXYMPjOd+Crr8jZtcsY1ucgkRZym7NvN+nmKcSESBZyh5golvKbnlu22BxM5JqadKC/sKK2nvLaBs7pmdVqCrFQpK9RuI0kxETyN5t7bd5scyCG9h6gbDtURm19I8VlNVy/cBM7jxp9h4dPVpPjS6dgUG6bT5JlrRbhRtJkTqSiIvD5yN6715g4tl8/W8KItHm8eX9p85IDjbrFZLOhxj0Hfza4r1FGvAinkzvERMrMNGoSwRjb7HCThvYKDOnzpqkWw/tyfOH/LQ1eq0X6GoUbSEJMNIeU30SiYFAuPm8aebmZvHT75MD45rbKdazN5+C1WuTuULiBNJkTzSy/WbUKGhog3dl/BB6PCkwwa45xjpS5vyRD4RZyh5hoQ4dSPXCgsYD9JhnWJoSTSEK0QenEicYLFzSbQzEfylz3zKZAac3cJ9bz5dFymyMTIjoJS4hKqaFKqeeUUq8m6pxOddLlCdFkHca365sKGq2zRgjhQlElRKXU80qp40qpHUHb5yildiul9imlHgTQWu/XWt8azfmSxakxYyArCz7/HEpK7A6n06zD+Mx1n4Vws2jvEBcBc6wblFJpwJPA5cAo4Hql1Kgoz5NUdEYGzJhhvFm+3N5gIhC8Sp/JHMYHRimOAhmZIlwtqoSotV4PBK+vOQHY578jrAOWANdEc56k5JLym6G9uwZW6Qse2ZLt8wZKa87pmcWZhqYWI1M6utiVEHZTWkfXzFFKDQbe1lqP9r+fB8zRWi/wv/8BMBH438CvgJnAs1rr/2jjeLcDtwP07du3YMmSJRHHUllZSbdu3Tr9XRKlsrKSXlVVTJ4/n4bMTD544w2015kzVVdWVnK81rgNHNq7a2D7/hNVLfarb2iirtFYREsBfXN89MnuktA43fBnD+6J1S1xQsdjnTZt2jatdWHw9ngUwakQ27TWuhS4s70Pa60XAgsBCgsLdVFRUcQnXrduHR3Z3y7r1q1j8pVXwr//O+lffMGl6enGsD4HWrduHX8qMRLb0nnNd3uPPbGe8toGcnzpZPu8HCuv5WBpNQA+r4fFCyaErD80m95L75jc4nUs4nTDnz24J1a3xAmxizUeT5mLgYGW93nAkTicx/1c0mwOZp24Ydc3FRwrr+XwSSMZKuAXV54nxdjCleKRELcCw5RSQ5RSGcB84M2OHCCp1mUOx0GLT3WEdeKGJg0nq+oC7zVQVl1nW2xCRCPaspuXgE3AcKVUsVLqVq11A3A3sBLYCbystf6iI8dNqnWZw7n4YsjOhi+/hEOH7I6mTcHzGgavy9Kza0aL9zKRg3CraJ8yX6+17q+19mqt87TWz/m3L9Nan6u1/o7W+lexCTUJZWTAzJnGa4eW31TXNbaa1zB4XZa+OT5G9MsmI93DiH7ZUTWXw63UJ6v4iXhz5NC9lGkyg+ObzVVnQq+hErwuS7bPS5d0j6zrLFzNkQkxZZrM0JwQ330XzpyxN5YQunZJl3kNRcpwZEJMKQMGwAUXQHU1rF9vdzStZGWkhZ3XsKPF19b+SFlzRTiNJEQncHj5jdk8jraUxlquY12jRdZcEU7hyISYUn2I4Ph+xFhpsc6KZY0WWd9ZOIUjE2JK9SECTJ4M3bvDnj3w1Vd2R9Npo/rnhG0+t1hnxbJGi/RNCqdwZEJMOenpMHu28dqh5TexYC3Xsa7RImuuCKeQhOgUKdJstvZHhuqbDPegRR7CiHiThOgUc/zTSq5dCzU19sZik3CL28vC9yIRHLnkm1LqKuCq/Px8u0NJnH79oKAAtm2Ddeua7xgdKrivsL3Sm3AjTMzflZyqaXNxe1n4XiSCI+8QU+6hisnh5TfxluNr/vc5+EGLLHwvEsGRCTFlWfsRo5y4142yfV6yMtJCPmiRhe9FIkhCdJIJE6BnT9i/H/butTsaW6R5VJtF4LEqEBeiLZIQnSQtrfnhioOazbI2ikgVjkyIKTdSxSpFym+EcCJHJsSUfagCRoG2UvC3v0FVVfv7CyFixpFlNymtTx+jL3HLFnjvPbjqKrsjitp1z2ziy6PlgeVMO/pZ6NhCVMElPuE+G8uFroT7OfIOMeWlePmNEHaRhOhESVx+Y31AIw9rhNNIQnSiggKj6Xz4MOzcaXc0MdHYpNsdo3ysvJYzDU1U1NbbEKEQkhCdyeNJqqfNFbX1VNc1tjtG+WBpNXUNTez6pkLGKgtbODIhpnTZjSmJEmJ5bUPgdfBksNYxyqYmjUwYK2zhyISY0mU3plmzjDvFDRugvNzuaKIS6Rhlk6ztLOziyIQoMIbwTZ4M9fXGinwuFukY5cG9smKytrMQnSV1iE42dy588IHRbL72WrujiUp7Y5TN9Zz75vjafPIsT6RFvMkdopOZ/YjLlydd+Y0QTiR3iE42diz07w8lJbB9O4wZY3dEEYnH6I+2jmluv2t48/vgUTGRjlwJF3d730lGvCQHuUN0MqWS4mlzRW19p+sLzRrFjnzWrHk8Vl7b7mdlnRZhJQnR6VyeEM06w87UF1prFHd9UxFRUrTWPB4srQ772eB1WqQgXDgyIUodosXMmcY8iRs3wqlTdkfTYdY6w47WFwZ/1lrP2JZQ+7T12eB1WiI5vkhujkyIUodo0b07XJ3lnAYAAA+XSURBVHIJNDbC6tV2R9Nh1jrDjtYXBn/WWs/YllD7tPXZ4HVaIjm+SG6OTIgiiItnvzHrDDtTX2itURzRLztQmhOOteZxcK+ssJ8NXqclkuOL5CYJ0Q2s5TdNTfbG0gnZPi9d0j2dSjjmOiod+axZ89g3x9fuZ2WdFmElCdENRo+GvDw4dgw+/dTuaIRIWpIQ3UApVzebhXALSYhu4fLyGyHcQBKiW1x2GXi9xlorpTI1lhDxIAnRLbKzYepU46HKqlV2RyNEUpKE6CbSjyhEXElCdBOzH3HFCkeX38RjfHBbxzS3V9c1Bt4Hj5u2rtkSLq5wcbf3nWRMdHKQhOgmI0bA4MFw4gR89JHd0YQUPD44FgmiorY+5DGt5zpwoor/2XK4xbjpitr6Fp81xzaHiitc3O19p3h8Z2EPR45VUkpdBVyVn59vdyjOYpbf/P73RrN5wgS7I2oleHywOXZ5VP+cDk38ap2yq7y2ocUx713yCQN6ZFJyqsYy1lnzf5fvDDn2OXjNltp6Iy5rMXZw3OY5lt4xOeR3KhiUG4hx6rl9Qv4++LuEmxosGaYPS4bv4Mg7RBnLHIbDy2+CxwfHYm2UHF96yDHH1rHHCujZNaPV2GfrZ02hxlSHG9fc3neKx3cW9nBkQhRhTJsGXboYTebjx+2OppXg8cGxGBKX7fOGHHNsHbfs86bRN8fXYty0uTRBJGu2hBvX3N53isd3FvaQhOg2XbtCUZGxpMDKlXZHE1I8xge3dUxz3LLHf4sWaty0+dm+Ob6wY6rDxd3ed5Ix0clBEqIbSfmNEHEhCdGNzH7ElSuhQSY1FSJWJCG60bBhkJ8PZWXw4Yd2RyNE0pCE6FbSbBYi5iQhupXDy2+EcCNJiG516aWQmQmffAJHj9odjRBJQRKiW2VmwvTpxusVK+yNRYgkIQnRzaQfUYiYkoToZmY/4qpVUC+LrAsRLUmIbjZkiDEDTnk5bNrU/v5CiLAcOduN6IC5c2HXLqPZPHWq3dGEFM3sJ+ZnrbPftCd4Zp1wn3XzzCwi9uQO0e2k/EaImJGE6HZTphgTPmzfDsXFdkcjhKslLCEqpboqpV5USv1BKXVDos6b9Lp0gRkzjNfLl9sbixAuF1VCVEo9r5Q6rpTaEbR9jlJqt1Jqn1LqQf/mvwNe1VrfBlwdzXlFECm/ESImor1DXATMsW5QSqUBTwKXA6OA65VSo4A84Gv/bo1RnldYmf2Ia9ZAXZ29scRIJAs+HSuvDSwo1dSkKTlV02JxqVAa/fu1t+5JJItGRbqwVCyP5WTJ8B2iSoha6/XAyaDNE4B9Wuv9Wus6YAlwDVCMkRSjPq8IMnAgjB4NlZWwYYPd0UQt0gWfDpZWU9fQxM6jFdTUN1JcVhNYXCqUitp6qusa210MKpJFo4L3aeucnTmWGxNKMnwHAKW1bn+vcAdQajDwttZ6tP/9PGCO1nqB//0PgInAA8DvgFpgg9Z6cRvHux24HaBv374FS5YsiTiWyspKunXr1unvkijxiHPoM89wzpIlfP397/PVXXfF7LidiXX/iSojpt5dO3XObyvO8E15LWCsldI3x0fFGWPex+wu6YHfWfXNhGM1xuuMNA/D+2W3iqe+oYm6xqYWx+2T3aXVftZzBJ/f/E7BMXrTPHjTPa2+c6jvkqnqW1zTUPtY47JLR/7s7f4OHf3/dNq0adu01oXB2+NRh6hCbNNa6yrg5vY+rLVeCCwEKCws1EVFRRGfeN26dXRkf7vELc4lSxi4fTsDY3jszsT6lLn62rzO1fhtO1TGA09vpEmDz+th8YIJPLZiFwD/OmdE4HcmBfz0/Ab+c7uxoNSIftksm99ck2nGU1Fbz5dHK4Dm41qn/Df3s54j+PzmdwqOcWjvrmT7vK2+c6jvUnHgsxbXNNQ+TliKoCN/9nZ/h1j9nYpH07UYGGh5nwccicN5hNXFF0N2NuzcCQcP2h1NVMIt2mT9nblo1Mj+2WR6jcWmzMWlQrEuShVuMahIFo0KtyhVtMdyQjLsqGT4DhCfhLgVGKaUGqKUygDmA2925ABKqauUUgtPnz4dh/CSlNcLs2YZr5Og/CaSBZ+si0Z5/ItNtZWYTOaiVO39hY1k0ahIF5aK5bGcLBm+Q7RlNy8Bm4DhSqlipdStWusG4G5gJbATeFlr/UVHjivrMneSlN8IEZWo+hC11te3sX0ZIH8rE22OvwLq3XehthZ8PnvjEcJlHFn+Ik3mTjr7bBg7FmpqYP16u6MRwnUcmRClyRwFaTYL0WmOTIgiCjL7jRCdJgkx2UyaBD16wN69sG+f3dEI4SqOTIjShxiF9HSYPdt4nQTlN0IkkiMTovQhRkn6EYXoFEcmRBEl8w5x7VqorrY3FiFcRBJiMurbFwoL4cwZWLfO7miEcA1JiMlKms1CdJgjE6I8VIkBa/lNlFO8CZEqHLkMqdb6LeCtwsLC2+yOxbXGj4deveDAAdizB4YPT9ip47G0Z1vHNJccff7Vd9rcLxZLkrb3nWQ50+TgyDtEEQNpac1jm6XZLEREJCEmM+lHFKJDJCEms1mzQCljoofKSrujEcLxHJkQ5aFKjPTuDRMnGivxvfee3dEI4XiOTIgyUiWGpNksRMQcmRBFDEn5jRARk4SY7MaNg7POgq+/hi+/tDsaIRxNEmKy83hkjkQhIiQJMRVIP6IQEZGEmApmzjTuFDdsgPJyu6MRwrEcmRCl7CbGcnPhoougoQHWrLE7GiEcy5EJUcpu4kCazUK0y5EJUcSB+WBl+XIpvxGiDZIQU8UFF0D//nDkCHz+ud3RCOFIkhBThVLSbBaiHZIQU4kkRCHCkoSYSmbMMJYp3bgRysrsjkYIx5GEmEpycuCSS6CpCVavtjsaIRzHkQlR6hDjSJrNQrTJkQlR6hDjyFp+09RkbyxCOIwjE6KIo/POg4ED4fhx+OQTu6MRwlEkIaYaKb8Rok2SEFORJEQhQpKEmIqmT4eMDNiyBU6csDsaIRxDEmIq6tYNpk41xjSvWmV3NEI4hiTEVCXNZiFakYSYqszymxUroLHR3liEcAhJiKlq+HAYMgRKS+Gjj+yORghHkISYqqT8RohWJCGmMkmIQrTgyIQoY5kTpKgIunQxmszHjtkdjRC2c2RClLHMCZKVBdOmGa9XrrQ3FiEcwJEJUSSQNJuFCJCEmOrM8puVK41lSoVIYZIQU11+PgwbBqdOGUP5hEhhkhCFNJuF8JOEKCQhCuEnCVEYEz1kZcGnnxrrNguRoiQhCvD5jCnBwBjbLESKkoQoDNJsFkISovAzy29Wr4b6entjEcImkhCFYfBgGDkSysuNheyFSEGSEEUzaTaLFCcJUTSThChSnCRE0eySS4z1VnbsgK+/tjsaIRJOEqJolpEBM2YYr5cvtzcWIWwgCVG0JM1mkcISlhCVUkOVUs8ppV5N1DlFJ5jlN2vWwJkz9sYiRIJFlBCVUs8rpY4rpXYEbZ+jlNqtlNqnlHow3DG01vu11rdGE6xIgLw8OP98qKqCDRvsjkaIhIr0DnERMMe6QSmVBjwJXA6MAq5XSo1SSp2vlHo76OesmEYt4kuazSJFRZQQtdbrgZNBmycA+/x3fnXAEuAarfV2rfWVQT/HYxy3iCdJiCJFKa11ZDsqNRh4W2s92v9+HjBHa73A//4HwESt9d1tfL4X8CtgJvCs1vo/2tjvduB2/9vhwO5IvwzQGzjRgf3t4pY4wT2xuiVOcE+sbokTOh7rIK11n+CN6VEEoEJsazO7aq1LgTvbO6jWeiGwsFMBKfWR1rqwM59NJLfECe6J1S1xgntidUucELtYo3nKXAwMtLzPA2QyPSGEa0WTELcCw5RSQ5RSGcB84M3YhCWEEIkXadnNS8AmYLhSqlgpdavWugG4G1gJ7ARe1lp/Eb9QI9KpprYN3BInuCdWt8QJ7onVLXFCjGKN+KGKEEIkOxm6J4QQfq5JiKFGyyileiqlViul9vr/m+vfrpRSv/GPoPlcKTXOAbH+m1KqRCn1qf9nruV3D/lj3a2Ump3AOAcqpdYqpXYqpb5QSt3r3+6o6xomTideU59S6kOl1Gf+WB/xbx+ilNriv6ZL/f3uKKW6+N/v8/9+sM1xLlJKHbBc07H+7bb+nfLHkKaU+kQp9bb/feyvqdbaFT/AVGAcsMOy7THgQf/rB4FH/a/nAssxSoMmAVscEOu/AfeF2HcU8BnQBRgCfAWkJSjO/sA4/+tsYI8/Hkdd1zBxOvGaKqCb/7UX2OK/Vi8D8/3bnwbu8r/+R+Bp/+v5wFKb41wEzAuxv61/p/wx/BT4H4x6aOJxTV1zh6hDj5a5BnjR//pF4LuW7X/Uhs1AD6VU/8RE2masbbkGWKK1PqO1PgDswxgFFHda66Na64/9ryswHo4NwGHXNUycbbHzmmqtdaX/rdf/o4HpgDmxSfA1Na/1q8BlSqlQNb6JirMttv6dUkrlAVcAz/rfK+JwTV2TENvQV2t9FIy/NIA5ZnoAYJ3htJjwf4ES5W5/c+N5sxmKQ2L1NysuxLhTcOx1DYoTHHhN/U27T4HjwGqMO9RT2qjMCI4nEKv/96eBXnbEqbU2r+mv/Nf0v5VSXYLj9Ev0n/3/A/4VaPK/70UcrqnbE2JbOjSKJkGeAr4DjAWOAv/p3257rEqpbsBfgH/SWpeH2zXEtoTFGiJOR15TrXWj1nosxmCFCcDIMPHYFmtwnEqp0cBDwAhgPNATeMDuOJVSVwLHtdbbrJvDxNPpWN2eEI+Zt+3+/5qTSDhuFI3W+pj/f8Am4A80N+FsjVUp5cVIMou11q/5NzvuuoaK06nX1KS1PgWsw+hz66GUMofKWuMJxOr/fXci726JdZxz/N0TWmt9BngBZ1zTi4GrlVIHMSaRmY5xxxjza+r2hPgm8CP/6x8Bb1i2/9D/ZGwScNpsAtolqL/lWsB8Av0mMN//ZGwIMAz4MEExKeA5YKfW+r8sv3LUdW0rTode0z5KqR7+15nADIw+z7XAPP9uwdfUvNbzgPe0/2mADXHusvxDqDD65KzX1Ja/U1rrh7TWeVrrwRgPSd7TWt9APK5pop8UdfYHeAmjWVSP8S/ArRj9Au8Ce/3/7ambn6A9idF3sx0odECsf/LH8rn/D6y/Zf+f+WPdDVyewDgvwWhKfA586v+Z67TrGiZOJ17TMcAn/ph2AL/wbx+KkZT3Aa8AXfzbff73+/y/H2pznO/5r+kO4M80P4m29e+UJe4imp8yx/yaykgVIYTwc3uTWQghYkYSohBC+ElCFEIIP0mIQgjhJwlRCCH8JCEKIYSfJEQhhPCThCiEEH7/H6Xlbj11tNXKAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "<table>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td/>\n",
+ "<th title=\"Variable name\">\n",
+ "Name\n",
+ "</th>\n",
+ "<th title=\"Value of parameter\">\n",
+ "Value\n",
+ "</th>\n",
+ "<th title=\"Hesse error\">\n",
+ "Hesse Error\n",
+ "</th>\n",
+ "<th title=\"Minos lower error\">\n",
+ "Minos Error-\n",
+ "</th>\n",
+ "<th title=\"Minos upper error\">\n",
+ "Minos Error+\n",
+ "</th>\n",
+ "<th title=\"Lower limit of the parameter\">\n",
+ "Limit-\n",
+ "</th>\n",
+ "<th title=\"Upper limit of the parameter\">\n",
+ "Limit+\n",
+ "</th>\n",
+ "<th title=\"Is the parameter fixed in the fit\">\n",
+ "Fixed\n",
+ "</th>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#FFFFFF;\">\n",
+ "<td>\n",
+ "0\n",
+ "</td>\n",
+ "<td>\n",
+ "loc\n",
+ "</td>\n",
+ "<td>\n",
+ "100.0\n",
+ "</td>\n",
+ "<td>\n",
+ "1.0\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td>\n",
+ "1\n",
+ "</td>\n",
+ "<td>\n",
+ "scale\n",
+ "</td>\n",
+ "<td>\n",
+ "10.0\n",
+ "</td>\n",
+ "<td>\n",
+ "1.0\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "</table>\n"
+ ],
+ "text/plain": [
+ "-------------------------------------------------------------------------------------------\n",
+ "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
+ "-------------------------------------------------------------------------------------------\n",
+ "| 0 | loc | 100.0 | 1.0 | | | | | |\n",
+ "| 1 | scale | 10.0 | 1.0 | | | | | |\n",
+ "-------------------------------------------------------------------------------------------"
+ ]
+ },
+ "execution_count": 42,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "m = Minuit(blh, \n",
+ " loc=100, scale= 10,\n",
+ " errordef=0.5, #remember up is 0.5 for likelihood and 1 for chi^2\n",
+ " pedantic=False)\n",
+ "\n",
+ "# Show() is the same thing as draw(). But show the figure immediately.\n",
+ "# For all parameters and return vars:\n",
+ "# https://probfit.readthedocs.io/en/latest/api.html#probfit.costfunc.UnbinnedLH.draw\n",
+ "plt.figure(figsize=[5,5])\n",
+ "plt.yscale('log', nonposy='clip')\n",
+ "plt.ylim([0.1,1000])\n",
+ "blh.show(m, print_par=True)\n",
+ "m.get_param_states()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 43,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/mauro/anaconda3/envs/fits/lib/python3.7/site-packages/ipykernel_launcher.py:1: LogWarning: x is really small return 0\n",
+ " \"\"\"Entry point for launching an IPython kernel.\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAE1CAYAAACflJmiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxN9R/H8dd39sGMxtgGNYyxrzViUsmSslcoNElZUhQttn4qqaQdCS0opBnSQoiSpGRfI0KWjH0ZZozZ5/v74y7urGbMcs6983k+HvfRveeee87nHrw795zvorTWCCGEADejCxBCCLOQQBRCCCsJRCGEsJJAFEIIKwlEIYSwkkAUQggrCUQhhLCSQBRCCKtiC0SlVD2l1MdKqUVKqaeKa79CCJFXBQpEpdRspdQZpdTuTMs7KKX+UUodVEqNAdBa79VaPwk8BDQryH6FEKIoFPQM8Qugg+MCpZQ7MA3oCNQH+iil6lvf6wb8AfxSwP0KIUShK1Agaq3XAhcyLW4OHNRaH9JaJwNRwH3W9ZdorVsCEQXZrxBCFAWPIthmVeCYw+tooIVSqjXQHfAGluf0YaXUE8ATAL6+vmE33nhjnnecnp6Om5v57xM5S53gPLU6S53gPLU6S52Q/1r3799/TmtdIcsbWusCPYDqwG6H1w8CMx1e9wWmXs+2w8LCdH78+uuv+VrfKM5Sp9bOU6uz1Km189TqLHVqnf9agS06m8wpiviPBhxP66oBJ4pgP0IIUaiKIhA3A7WUUjWUUl5Ab2BJfjaglOqqlPr00qVLRVCeEEJkr6DNbiKB9UAdpVS0UmqA1joVeBpYCewFFmqt9+Rnu1rrH7TWT5QtW7Yg5QkhRL4U6KaK1rpPDsuXk8uNk6KiUlKKe5dCCBfiHLeQ8uKrr2jRty8cOWJ0JUIIJ2XKQMz3NUSt4euv8Tl9mpNderD137NFW6AQwiWZMhDzfQ1RKXaMe4/LAeUI2rONDY8/y9ajMUVbpBDC5ZgyEK/HuhjNT089RzqKJ/+I4uj3K4wuSQjhZFwmEMNDAjnZoDEfh/fEXadz78QRECNniUKIvDNlIF5PO8Sw4ACCyvow6Y4ItgfVofTpE8Q88hidJv/GHW+vlp/QQohrMmUgXm87xLR0TYq7B8O6jSTOy5eA5UtovPIbomMSiJi5QUJRCJErUwbi9Srt7YGPpxvHb6jMy/cMAeCVXz6l5vljJKWk8822aIMrFEKYmUsFYikvd+YPDOf2WuX5cOU0vm3QhlIpSXy45F08U1NYtDW6SM8S+/fvT8WKFWnYsGGG5RcuXKB9+/bUqlWL9u3bExcXB1gG1hg2bBihoaE0btyYbdu2ZbvdDh060KRJExo0aMCTTz5JWlpanmuaOHEioaGh1KlTh5UrV2a7TkREBHXq1KFhw4b079+fFGsD9zVr1tClSxeaNm1K06ZNee211+yfWbFiBXXq1CE0NJS33nrLvvzOO++0r1+lShXuv//+PNcqhOGyG/HB6AfQFfg0NDQ0XyNY2Ea82HLkglaePrrxcwv1kRsqaw36s2b36ZAxS/VHqw/ka5v58dtvv+mtW7fqBg0aZFg+cuRIPXHiRK211hMnTtS9e/fWWmu9bNky3aFDB52enq7Xr1+vmzdvnu12L126pLXWOj09XXfv3l1HRkZmWSc4ODjLsj179ujGjRvrxMREfejQIR0SEqJTU1OzrLds2TKdnp6u09PTde/evfX06dO11pbjGR4enmX91NRUHRISov/991+dlJSkGzdurPfs2ZNlve7du+s5c+Zk+50KmyuPzGIUZ6lTa3OPdlNguoB9mcOCA/DxdOeJrjez972PeQHF5C2LiZnRn9g9v9nXe+edd2jUqBFNmjRhzJgxBa67VatWlCtXLsvyxYsX069fPwD69evHunXr7MsfffRRlFKEh4dz8eJFTp48meXz/v7+AKSmppKcnIxSKk/1LF68mN69e+Pt7U2NGjUIDQ1l06ZNWdbr1KkTSimUUjRv3pzo6NwvLWzatInQ0FBCQkLw8vKid+/eLF68OMM6cXFxrF69Ws4QhVMxZSAWBjcFQ9uEEn9DKpuDq7MT2OWexoy3xnHy5El+/PFHvv/+ezZu3MjOnTsZNWpUlm3Mnz/f/vPP8dGzZ8981XL69GmCgoIACAoKIsbaHOj48eM4DoBbrVo1jh8/nu027r33XipWrIifn599/xMmTLDXdOLECfvzoUOH5nv7ACkpKcybN48OHa7OCvH333/TpEkTOnbsyJ49e/K83e+++4527drZw1wIZ1AUI2abyh9//MFjL/0P96++otqvv3JXUBCbN27kt99/5/HHH6dUqVIA2Z7ZRUREEBFRdLMdWM7cM8rp7G/lypUkJiYSERHB6tWrad++PWPHjmXs2LEAVK9enR07dlz39gGGDBlCq1atuPPOOwG45ZZbiIqKomPHjixfvpz777+fAwcO5Gm7kZGRDBw4MMd9CWFGLnuGaKO1Bjc3mDcPAgPh5En44Qe01tf86VlYZ4iVKlWy/xQ+efIkAQEBgOXM6tixq7MtREdHU6VKlRy34+PjQ7du3bL8PM1JfrY/fvx4zp49ywcffGBf5u/vj6+vL2D5WZ2SksK5c+euud3z58+zadMmOnfunKc6hTALUwZiYQ4Q26pVKxYsWEBa5cps/N/rrAXC5sxl+4Fknn1tEuv2WX7qXbiQea4syxnijh07sjwWLVqUrxq6devGnDlzAJgzZw4tW7a0L587dy5aazZs2EDZsmXtP61tLl++bA/T1NRUli9fTt26dbPs40g2o/x069aNqKgokpKSOHz4MAcOHKB58+ZZ1ps5cyYrV64kMjIyw7wUp06dsp8Nbtq0ifT0dAIDA7n11ls5cOAAhw8fJjk5maioKLp162b/3Ndff02XLl3w8fHJ13ESwnDZ3Wkxy6Mgc6qULl1aa225MztixAgdUruu9qwQrIdUv1lr0AfLVdMVbo/QnoE36orV6+jHhzyXr31lp3fv3rpy5craw8NDV61aVc+cOVNrrfW5c+d027ZtdWhoqG7btq1evHixvbYhQ4bokJAQ3bBhQ71582b7tpo0aaK11vrUqVO6WbNmulGjRrp+/fr66aef1ikpKVprrd944w3dpEmTLI8hQ4bYt/PGG2/okJAQXbt2bb18+XL78o4dO+rjx49rrbV2d3fXISEh9s+PHz9ea6311KlTdXBwsG7cuLFu0aKFXrdunf3zy5Yt07Vq1dIhISH6jTfeyHAc7rrrLv3jjz8W+HjmhyvfETWKs9SpdeHdZTY89HJ7FOYkUx+tPqCDRy/VtZ//Ru8rf5PWoCMb36ODRy/VwaOX6jovLddbjlzI1/6ulyv/RTOKs9SptfPU6ix1au3izW6KQnhIIG4Kkjy9GdZtFEnunvTe9ROd9v0BQEpqOhsOnTe4SiGEkVz+LrNNWHAAXz/Zkg2HzhNQqhEr3U/Q7bM3mbhiKjuq1OZCYGXCQwKNLlMIYaASc4YIllAc2iaUh1vcRLdP3oCuXSmbFM+iP6Yz/7FmhAUH2Nc9cuRIli54hemDDz6gfv36NG7cmHbt2nH06FH7e+7u7vY72o43KxwlJSXRq1cvQkNDadGihf2mSuY7425ublma4+Qkp206OnbsGG3atKFevXo0aNCAKVOm2N979dVXqVq1qn3fy5dfnVYnpy6EU6ZMoWHDhjRo0IDJkyfnqU4hikx2v6PN8ijyierPntU6KMhyKfXVVzO8dfjw4Sxd8ArLr7/+qlevXq3j4+O11lpPnz5dP/TQQ/b3bTeEcjNt2jQ9ePBgrbXWkZGRGT5vs2vXLl2jRo0syw8fPqzvuuuuPG0z8zE9ceKE3rp1q9Za69jYWF2rVi17t71x48bpd999N8t2c+pC+Ndff+kGDRro+Ph4nZKSotu1a6f3799/ze+eHVe+3mUUZ6lTaxe/hljQZjfx8fF07tyZJk2a0LBhQxYsWADA5s2badmyJU2aNKF58+bEeXtz5N13uRO45dVXuaV2bf78888s20tLS2PkyJHceuutNG7cmE8++aQgXw+ANm3a2BuFh4eHX7O7XGaO3QF79uzJL7/8YrlL5iAyMpI+fbKdGPG6txkUFMQtt9wCgJ+fH/Xq1cu194ttu9l1Idy7dy/h4eGUKlUKDw8P7rrrLr777rs81ytEYTNlIOoC9mVesWIFVapUYefOnezevZsOHTqQnJxMr169mDJlCjt37mTVqlX4+vpS8YEH+PmFF9gGzDp/iYGDnsyyvVmzZlG2bFk2b97M5s2b+eyzzzh8+HCW9RxHenF8rFq1Ktd6Z82aRceOHe2vExMTadasGeHh4Xz//ffZfsax+5yHhwdly5bl/PmMN4UWLFiQIRAfeOABmjZtSqdOndiyZYu9vs8//zzHbcbGxuZY95EjR9i+fTstWrSwL/voo49o3Lgx/fv3v2YXxYYNG7J27VrOnz/PlStXWL58eYYG30IUN5e8qdKoUSNGjBjB6NGj6dKlC3feeSd//fUXQUFB3HrrrcDVARPi4+Ppc/AoBzy8KH3hDIdizrH1yAUcb6/89NNP7Nq1y94g+9KlSxw4cIAaNWpk2O/vv/+e71q//PJLtmzZwm+/XR104r///qNKlSocOnSItm3b0qhRI2rWrJnhc5nP3CBj97mNGzdSqlSpDNdBbWdfR44c4bHHHmPNmjXX3GZOLl++TI8ePZg8ebL9WD711FO8/PLLKKV4+eWXeeGFF5g9e3aOtdarV4/Ro0fTvn17ypQpQ5MmTfDwcMm/ksJJmPIMsaBq167N1q1badSoES+++CKvvfYaWmffVW/SpEkkeZUl6PEPWePpQ6pO59L0jD+JtdZMnTrV3lPl8OHD3HPPPVm2ld8zxFWrVjFhwgSWLFmCt7e3fbmtG1xISAitW7dm+/btWT7r2H0uNTWVS5cuZeiPHRUVla+fyzltM7vBGVJSUujRowcRERF0797dvrxSpUq4u7vj5ubGoEGD7CPr5NbVb8CAAWzbto21a9dSrlw5atWqla+ahShMLhmIJ06coFSpUjzyyCOMGDGCbdu2UbduXU6cOMHmzZsBy/BUtn/0jWsHEx1YjV71WpEGtJr6Ot8vWktCimUg1nvvvZcZM2bYB07dv38/8fHxWfb7+++/Z9vV7+67786y7vbt2xk8eDBLliyhYsWK9uUxMTEkJSUBcO7cOdatW0f9+vWzfN6xO+CiRYto27atPfDT09P5+uuv6d27d7bHp3r16lnODq+1TRutNQMGDKBevXo8//zzGd5zHLrsu+++s5+d5taF8MyZM4DlrPjbb7/Nd4gLUZhc8vfJX3/9xciRI3Fzc8PT05MZM2bg5eXFggULeOaZZ0hISMDX15dVq1YxZMgQevToQWn3xegGYfjsdscjMYF6b7/McWDr0RgGDhzIkSNHuOWWW9BaU6FChRyv7eXVyJEjuXz5Mg8++CAAN910E0uWLGHv3r0MHjwYNzc30tPTGTNmjD0QX3nlFZo1a0a3bt0YMGAAffv2JTQ0lHLlyhEVFWXf9tq1a6lWrRohISEZ9vnAAw9ke+1z+PDhPP7449lu87///uPEiRMMHDiQ5cuXs27dOubNm0ejRo1o2rQpAG+++SadOnVi1KhR7NixA6UU1atXt998atCgAQ899BD169fHw8ODadOm4e7uDkCPHj04f/48np6eTJs2zT7whRCGyO7Ws1keRd7sJpMtRy7ohs8u1EfLVtIa9KfNHyiSEbZduTmDUZylTq2dp1ZnqVNrF292Y5QNh84T512K4V1HkqrcGLTpO+45lrdGzUII5yeB6CA8JBAfTzd2Vq3LlDstA8OWG/oEO7fuN7gyIURxMGUgFuZ4iPkRFhzA/IHhvHBvHapMHM/6mxoReDmGi70i2HpYBn4QwtWZMhB1ARtmF4Stv/OFpDSe6/wCMT5+3PXvFq68P6nYaxFCFC9TBqIZhIcEcjGwIi92GgZAi4/f4bHhnxbpvM5CCGNJIObA9vO50dOPsa1TL7zSUnjpy9cZOH2NhKIQLkoCMRdhwQGEhwQSUfdB9gfeROiFaEat/EQGkhXCRUkgXsOGQ+dJ8PRhWLeRJLl70mfnSgKWL5azRCFckATiNdimHthXsQZvtekPQOep43jm7cU0n/CzBKMQLkQC8RpsUw+MvLcOyU8N4Zeat1I2KZ5JS9/n3KUEImZukFAUwkVIIOaBrSlO97AbGd35Wc6UDqDFsd0MXb9QJqcSwoVIIOZDWHAAFUOqMf7BUQAMXxdJ85P7ZHIqIVyES452U5SWD28FtOKU9wkqfzKV9xe/y+kXegIySosQzs6UZ4hGdd3LjxMvjGVX5VpUuXiaE737sfXIBaNLEkIUkCkD0ciue3m1PjqOYV1HEO/pQ+e/f+PijM+MLkkIUUCmDERnEB4SyH+BVRnX3jIpVeuPXof9MiqOEM5MAvE62Zrj1HhhCL80aYP7lXjie/aC5GSjSxNCXCcJxAIICw4gvGZ5nm/zJMfKVqL0Xzs4NWyE0WUJIa6TBGIBbTh0nkvepRnedQSpyo3Kn0yFn3+m05S13PH2amm0LYQTkUAsIFvXvm1V6/Hh7ZYZ4xL6RHDmUDTRMdKTRQhnIoFYQLZriQ+3uIlptz3Exhsb4nv+LG8vmwxaS08WIZyIBGIhCAsOoOoNvqS5ufNslxe46FOGdv9upt+2pXh6uElPFiGchARiIbH9dD7pX4HRHSyjbP/v19m0jD/OOyv2GVydECIvJBALSVhwAHUr++Hv48HKOi35qkkHvK2jbJ87EyPXEYVwAtKXuRAtH96KrUdjePDjP3m97UBaRO+m5tn/6P/tVCKSFSHlS+Pn48lTdYyuVAiRHTlDLGS2M8XASgEs/98kktw9iNixgrt2/8H+M5c5fjGBK8lpRpcphMiGBGIRWD68FX+MbkvLHu14u/XjALy94kPKx5wlOiaBw+fi5Se0ECYkgViEwoID2NAlgnW1b+WGxMtMXvoebulpaK2lKY4QJiSBWMSWP3sX/0ycytnSNxB+bDdPbVgEKGmKI4QJSSAWgyZhtXmhy/MAPPfHfBqePERYsAwoK4TZFFsgKqXuV0p9ppRarJS6p7j2awZhwQEMf/tptj84AA+dTtg7b/HY5J+NLksIkUmBAlEpNVspdUYptTvT8g5KqX+UUgeVUmMAtNbfa60HAY8BvQqyX2cUFhzAzV9Oh1tuwe/saR6e87aMsi2EyRT0DPELoIPjAqWUOzAN6AjUB/oopeo7rPKS9f2Sx8uL3e99TIq3N/fs+IWFw9+k05S19PpkvdGVCSEoYCBqrdcCmU9zmgMHtdaHtNbJQBRwn7J4G/hRa72tIPt1Zr8RwG+PDgLg5RXTSdj7D4fPXZZmOEKYgNJaF2wDSlUHlmqtG1pf9wQ6aK0HWl/3BVoA+4F+wGZgh9b64xy29wTwBEClSpXCoqKi8lzL5cuXKVOmzHV/l+JwJTmN+Lg4mr77HrU2ruN0jVAWjZuI9vQipHxpSnm5G11iBs5wTMF56gTnqdVZ6oT819qmTZutWutmWd7QWhfoAVQHdju8fhCY6fC6LzD1erYdFham8+PXX3/N1/pG+WzBD7rl6EX6mH9FrUFPb9FTB49eqm9/6xejS8vCWY6ps9SptfPU6ix1ap3/WoEtOpvMKYq7zNHAjQ6vqwEnimA/Tiu0Yhk+fKotz3YdQZpyY/DGb2h5ZAfJqWlyTVEIAxVFIG4GaimlaiilvIDewJL8bMAZ5mUuqLDgAK7c2oKP7uiDG5pJyz4g9fRZ9p6Mk2uKQhikoM1uIoH1QB2lVLRSaoDWOhV4GlgJ7AUWaq335Ge72gnmZS4Mfj6ebOg9mBONmlHp8gXe/nEKWmvOxCXL1ANCGKBAw39prfvksHw5sLwg2y4JFgy+DYBdTWdTuu1ttD+4ib7blzHvli72qQekR4sQxUe67plA45aN+KjPKABeWj2LOmePoIFlu+TSqxDFyZSBWBKuIWY29vNXYOBAvNNSmLHsPUL93IlNTJWfzUIUI1MGYkm5hpjF5Mkcr3QTIaeP8OiiqUTHJPDgx39KKApRTEwZiCVW6dJU/fF70jw9eWT7cu7Zv550Dd9siza6MiFKBFMGYkn8yWx3882cGD0OgLd//JBKcedYtDVazhKFKAamDMQS+5PZ6sbxL/JnrWYEJMYxaekH6JQUGWFbiGJgykAs8dzcKP3Vl5wrdQMt/9vFU5u/lRG2hSgGEogm1aRZHab0/R8Aw3/7km9nfGNwRUK4PglEE9t/8x0sbdcL9/Q0npk1DmJjjS5JCJdmykAs0TdVHCwYfBtdls3hYLVaVD53gvOPDTK6JCFcmikDsaTfVHG09dQVBnd8gSue3gR+t5DDkz8xuiQhXJYpA1FcteHQef4tV43x7Z4AoNLo5+Dffw2uSgjXJIFocuEhgbgpWND4HpbXuZ1SyQnE93gIUlKMLk0IlyOBaHJhwQHUrexHBX9vXuzwDNH+FSi9cxsLOj0uA8kKUchMGYhyUyUjPx9PvD3cueRThue6vECacuPBX76iwb4tRpcmhEsxZSDKTZWMFgy+jSm9b8ZNweYbGzKtZS/ctOap2ePZsf2g0eUJ4TJMGYgiK9tPZw93xdx2j7C5an0qxJ7nfO++MuG9EIVEAtHJpKZpziVpnu06gljv0rTbv4ElT71Cr0/W2x82mV8LIXIngehEYhNT7c+Pl63Ii/c+DcCLqz7jxuPSFEeIgpJAdCK264gACtjcvB0rW3TCJzWZpz59haTYyxnWj0tM4fjFBBk6TIg8kkB0IrbriNUCfKkX5EeN8mWY3GUo/5arSs3Th+keNcU+henWozHsOxVHdEyCzOAnRB6ZMhCl2U3O/Hw8qXqDL34+nsQlprDvsmZY15Eku3nw6LZlNN72OxEzN/DttmjSteUzthn8hBC5M2UgSrObnC0YfJv90blxFTSwp3Io79z1KADv/DiFgJizaLD/vPb0cJPxFIXIA1MGosgbW7c+gFm33s/a6jdTLiGW95dNokfTKvaf1/MHhsv8zkLkgQSiE7NdU/TycKNelbJ80v9lYkrfQMsjOwhbONP+81rCUIi8kUB0cpZufW74+XiSWrEyHz/2kuWNl16i5uG/7etJm0Qhrk0C0QXUD/JnweDbANjRqCUMHw6pqby5aCILHm54zc9LWAph4WF0AaJgbEGY4XnSLbBmDezcCU8/DXPm2NexBZ/j54QQFnKG6IJ6fbGN5+8fDb6+MHcufPWVvZF2XKKMoyhETkwZiNIOseCOB1WHyZMBuNJ/EPH7DhAdk8C+U3ESikLkwJSBKO0QC8mgQWy8uTWlkq4wefG7eKSlkq4z9okG6eInhI0pA1EUEqX45JExnL2hAjef/Ifh6yJRQHJqmj38pIufEFdJILogxzO++NL+TBvwKulKMXTDQlr8t4szccn28Ntw6Lx08RPCSgLRxWQ+4zsdm8iqCnVZ0L4vbloz6Yf3uSEh1h5+jr1dPD3cWLbrhDTBESWWBKKLcTzjS05J5+j5K0THJPBK0x7sDq5P0OXzvP3jh3i6K8JDAjOMoDN/YDh+Pp6AtE0UJZMEootxPONzc1NYs5EU5c6rvcYS71Oaew9sYGWpf+xd+vx8PLmUkMI7K/YZU7QQJiGB6GIcz/heu6/h1XBUcKVKNWZGjAIg+I2x8PffuWxJiJJHAtEF2QZ1eLjFTfZw/PrJliwf3oo/b23PmvBOkJDA0Xu68chHazJ8Vhpwi5JMAtHFZTfizejWAzlZoRrBx/8l4ttpLBh8G/WD/C0DzlpvyOw9Gcfhc5e5kpxmYPVCFC8JxBJmweDbqF69Mh8OHE+quwcdf10ES5cSl5hCdEyC/YaMBs7EJXP4XLy0TRQlhgSiC7KNqJ2bQ8H1iLpvMAAp/R7j/MGjWXqwAGitpW2iKDFMGYjSl7nw5BaOkXf0ZHPoLXheOM97P3yA0ukAlPJyx3ovBqVUlukHpEmOcFWmDETpy1z04hJT2HsmniH3PssFX3/uPLqDQZu+w01B9cBS1Auy3IypUb60jLgtSgxTBqIoWo7XC8+WKceozs8CMGrtXLomH8fPx9N+M6aUl7vB1QpRfCQQSxhb1z7H64V/1A1n8Z3d8UhPY8y81/FOvGJghUIYRwKxhHHs2gfg7+PB/IHhfP3QMxytWpOgs9H0X/CB/f1D5+Lt1wvze+1QrjUKZyOBWMJk6NqnoFqApY1iiqc3Uwa8Bj4+tF6/nJabf85xGzJ+onBVEogljGPXvrqV/eyDOQAcr1IDJk0CYPjX77PgnspZPu/YeFvGTxSuRgKxBLLdMHEMQ7vBg+H++yE2FiIiUGkZe6rEJqbK+InCZcmsewLINAvfzJmweTOsX0/TKpGsaD/M/pa/jwduCtK1ZfzEzG0UhXBmcoYosgoMhC+/BKVo8u1Cym/byNajMcQlphCbmMpN5UrZx0/MrY2iXGsUzkbOEEuwXLv3tW7NySHPETTtA8ZGTuC+0kGc8yqDxnIzpm5lP8KCA3Kc59nWvCddQ8TMDdcMTyHMQM4QS6C89HUG+P6+gZyqWZsqced4bdmHaG25eJjdzH2ZyVwtwhlJIIocNa9dmZ+GPk+cly8d9/9Jn50rAcsZor9P7j8uMs/VItcahTOQQBQ5CgsOILlqVT7o/hwA41Z/xu0pZ7I017FxbIidea4W+bksnIEEosiVm5vi79ZdoG9ffFKSmLBwIuXc0/P02ewGpxXCzCQQRd5Mm8apClWpHn2Ah7+bYXQ1QhQJCUSRN35+fDhgPKlu7nRavZCmf/1Jr0/W25vWdJqylr9PxhpdpRAFUmyBqJQKUUrNUkotKq59isL1b/X6LLjvCQCGzH0DjzOn7N349p2KI81x1AghnFCBAlEpNVspdUYptTvT8g5KqX+UUgeVUmMAtNaHtNYDCrI/Ybwf2kdA27aUjbvIsHlvotMt1xPTNaRKIAonV9AzxC+ADo4LlFLuwDSgI1Af6KOUql/A/QgD/X0y1n73WLu5wbx5xJYuS4sDWxi4+XvA0hRHgfRMEU6tQIGotV4LXMi0uDlw0HpGmAxEAfcVZD/COCHlS1M/yB9waD3+JecAABkYSURBVNBdpQoz+o0FLKNst4s/xk3lSpGUmp5hFJy8NgAXwiyUrffBdW9AqerAUq11Q+vrnkAHrfVA6+u+QAtgHDABaA/M1FpPzGF7TwBPAFSqVCksKioqz7VcvnyZMmXKXPd3KS7OUidYaj2TaGlhHVK+tH35oXPxhM+cQf2Vy7gUVIVFb07iiqc3YDlTrOTvQwU/72Kt05mOqTPU6ix1Qv5rbdOmzVatdbPMy4uiL7PKZpnWWp8HnrzWh7XWnwKfAjRr1ky3bt06zztes2YN+VnfKM5SJ1hqnXfcEmwLel4923tnylq+vGMos3f9TfWTh6k8/QuGth0KgI+nG/MHNs+2/aFj3+ec+kFfb53OdEydoVZnqRMKr9aiuMscDdzo8LoacKII9iMMYhu44dDldJ7s9AJJHl503vwjnff+jgJe6dJAGmMLp1QUgbgZqKWUqqGU8gJ6A0vyswGZl9ncHAdu2Fe+Ou/dMwiAiSs/ouql08RcSTawOiGuX0Gb3UQC64E6SqlopdQArXUq8DSwEtgLLNRa78nPdmVeZnPJPK5h5nlZfr6rOz/VCsc/KZ7JP7xH+E3y5yacU0HvMvfRWgdprT211tW01rOsy5drrWtrrWtqrScUTqnCCFeS07LMoZJ5XpZKZX2Z3e9FTpcJpNnxvYTNm3bd+8ttpj6ZxU8UNVN23ZOfzOYRn5T9HCqZ52XR5Svwv+4jSVcKXn8dfv/dqJKFuG6mDET5yWwepb098jyu4aYaTVlyzyOQng4RERAjDbSFczFlIArzKOXlnuu4hpkbXy/sNgiaN4djx+CJJyBTO1fH65Ey54owGwlEcU35Gdcwzd0DvvoK/Pxg0SKYPdv+nq25TnRMAn0+Xc/ekzK/szAXUwaiXEN0cjVrwvTplufDhsG+fUCmeVbSNLZzR5lzRZiFKQNRriE6p/pB/ld/Pj/yiOU64pUr0KcPJCVlnGfFXdm7NMmcK8IsTBmIwkVMnw4hIbBjB7z4YobmOpFP3Ea9IJlzRZiLBKIoOv7+luuJHh4waRKsWJHhemR21yZzu9EiN2FEUZNAFEWrRQt47TXL8379KBubebS4qxxvumS+0ZLbe0IUlqIY7abAlFJdga6hoaFGlyJykHmEmlxHrBk1it1zFtHwn208NecN3h76Xoa3bb1Pjl9MyNII3Hb2mN3E9/IzWxQ2U54hyk0VF+PuzkePjSOutD8379lAx9ULs13N3+fq/58z32iRie9FcTBlIArXExNQgY/7/g+AiO+mw/btWdbx8/GklJd7tjdaZOJ7URwkEEWx2dK0FT+1egCPtFTo0wfvpIQs67i7qRwbgcvE96KoSSCKayrMuVHm9hzGsaAa8M8/9Fs4uVC2KURhMWUgSk8V15Xi5c2Uga+Btzft1v1Ai62rjS5JCDtTBqLcVHFtx6rWhPcsd5qfmP82gRdOGVyREBamDEThWnp9sp6/T8ZmXDh0KFsa30GZK3E8M3s8bmmpOX42v4PC2j6Tl8/KoLPCkQSiMIZSfPzo/7hQtjz1Du7kgRVzja5ICHM2zBauy/HmzMwXOkLTKNLbt6fn0tnsrtOMv0tXN644UeLJGaIoFmnpOvt+yO3a8U2b3rjpdIbMfAWvy7HEJaYYU6Qo8SQQRZGLS0zhSnJajn2UXwp7iJ2Va1H54hnGL5vKvpOx0ldZGMKUgSjNblxLbOLVGyaZB4PdcOg8SW6eDOs2kstevnTZ9zs9dq2SAWOFIUwZiNLsxrXkpY/y0YAqvNL+SQBeXfUJdyFniKL4mTIQhWvJax/lba268EODuyidkkjDEU9CUpKBVYuSSO4yi2JxrT7Ktvmdvxk4lq5TjsK2bTB2LNTqYV+vsLoPCpETOUMUppLgWwYiI8HdHd5/nyZ7NhhdkihBJBBFoStw74/wcBg/HoAhc97AP/ZCjtvMvDy7XjF57bmS2/sF+axwHhKIosjFJaaQlJqev/aFY8bAXXdxQ+wFHp/1OnEJyXn+qK3N4+nYRI5fTMh1vzJPi3AkgSiKlG0ulOTUdPadist78Li7s2viVC76lKHlPxsJXzo/T4Hq2ObxyPkrRMcksO9UXLafzTxPizQIF6YMRGmH6Doc50JJ1+SrfeHvib6M7jgMgNFrPqfioX3X/Ixjm0ebdJ398szztGS3jihZTBmI0g7RdTjOheKmyNdcKOEhgfxcpyXzm3bAOy2V16ImWCa+z4Vjm0cbN5X98szztGS3jihZTBmIwnXY2hl6ebhRt7Jfvob/t312dvdn+K9SMDed+Q+eey7Xzzi2eaweWIpqAb7Urexnb9aT3fZt7SOzW0eULBKIosj5+Xji7eF2XYHj5+NJ+YoBfDToNVI8POHTT+Hbb3P9jK3NYyV/H6re4JvrfmWeFuFIAlE4haPVavFl96GWFwMHwrFjxhYkXJIEonAaK9o8CJ07Q0wMPPIIpKUZXZJwMRKIwnkoBZ9/DpUrw9q1MHGi0RUJFyOBKJxLhQow1zrdwKuvUuvQX8bWI1yKBKJwPu3bw4gRkJbGsFmv4ptw2eiKhIuQQBTOacIECAuj4vmTDJr/DmhtdEXCBUggikJXFP2Ds2zTywsiI0nw8uX2LasIWL7Cvl7mftO2z9r6NudUV251X+s7SZ9o1yCBKApV5v7BhREQcYkp2W5zq1d5Xrl7MAD1p37ID1//lqHfdFxiSobP2vo2Z1dXbnVf6zsVxXcWxjBlXyWlVFega2hoqNGliHzK3D/Y1ne5fpB/jgO8ZrfccSit2MTUDNscHrWdqjf4cvxiAtEN29Hq0Fa67vudkOGDcI94h3R3zwz9l9Mz/ZpOTLHU5dgYO3Pdtn0sGHxbtt8pLDjAXmOr2hWyfT/zd8ltgNu8rGN2rvAdTHmGKH2ZnVfm/sH56bucE38fj2z7HPv7eIBSjL13KLHlK9Dg5AFG/D4PuNp/2fGzNtn1qc6tX/O1vlNRfGdhDFMGonBemfsHF0aXOD8fz2z7HNv6LfsHVWDNsJGkubkzeOO3tP5vh73/suNnqweWyrFPdW79mq/1nYriOwtjSCCKQlcU/YNz2qat3/K5evVZ1Lk/AO/+8D5VUi5n+Wwlf59c+1TnVve1vpP0iXYNEojCZXzX8VH2hjahwuUYnpo7QZriiHyTQBQuQ7u5M7X/OGJ9yhD21590WLPI6JKEk5FAFC7lfLnKjOv6LAAR30zjpuiDBlcknIkEonA5q+rfwao7uuGVmsywWePwSk40uiThJCQQhUua++BwjlcO5saTh+m7aKrR5QgnIYEoXFKSty9TBownxcOTe9Z+x63bfzO6JOEEJBCFyzp6Y22+emAIAIO/nEi5mDMGVyTMTgJRuLTlbR9ie4Pb8IuP5enZ43FLl1G2Rc4kEIVrU4rp/cZy0b8cDQ5sp/+6r42uSJiYBKJwebH+5Zje7yUAhv46l1qHdhtckTArU452I1xLQUY/sX3WcfSba8k8sk6vT9azs0E4S+/uTZdVUTwz+1X430NQtqxTj8wiCp+cIYoSI/K+J/m7ck0qnTsBQ4ZI1z6RhQSiKDFSPb0Y3WMMiV4+8NVX8OWXRpckTKbYAlEpVVopNUcp9ZlSKqK49iuEoyPlb+TzXs9bXgwZAgela5+4qkCBqJSarZQ6o5TanWl5B6XUP0qpg0qpMdbF3YFFWutBQLeC7FeIgljTsjM8+CBcvgwPPwzJyUaXJEyioGeIXwAdHBcopdyBaUBHoD7QRylVH6gGHLOuJo3BRK7yMuHT6dhE+4RS6ema4xcTMkwulZ20dM3xS4nseOkduOkm2LwZXnklX/vPzzqFvS0zc4XvUKBA1FqvBS5kWtwcOKi1PqS1TgaigPuAaCyhWOD9CteW1wmfjpy/QnJqOntPxpGQkkZ0TIJ9cqnsxCWmcCXZsl7vr/ey74OPwc0N3nkHfvklT/vPaZ2c9nk923LGQHGF7wCgdAHvtCmlqgNLtdYNra97Ah201gOtr/sCLYDRwEdAIvCH1np+Dtt7AngCoFKlSmFRUVF5ruXy5cuUKVPmur9LcXGWOuH6aj10Lh6AkPKlr2ufZ+OSOBVrGaFGAZX8fYhLskwY5eftYX/PUSVfOJ1gee7l7kadyn5Z6klJTSc5LT3Ddpt9G0WNL74gKTCQRe98SJJ/2Qz7yLx/23fKXKOnuxueHm5ZvnN238VXpWQ4ptmtU8HPO/8HrpDl58/e6O+Q37+nbdq02aq1bpZ5eVG0Q1TZLNNa63jg8Wt9WGv9KfApQLNmzXTr1q3zvOM1a9aQn/WN4ix1wvXVOsM2+1rP62vjt/VoDKM//pN0DT6ebswf2Jx3VuwDYFSHuvb3bBTwfKNU3v/LMqFU3cp+LO/dKks9cYkp/H0yDri63Rod28GBA3ivW0edKTN496m3GdWxXo77t32nzDWGlC+Nn49nlu+c3XeJO7wzwzHNbh0zTEWQnz97o79DYf2bKoqfrtHAjQ6vqwEnimA/wkXlNmmT43u2SaPqBfnh6+lOtQBf++RS2bFNSpVhux4eMH8+lC1Ls11/cM9v3+Zp0qjcJqXK63fJzzpm5wrfAYomEDcDtZRSNZRSXkBvYEl+NqCU6qqU+vTSpUtFUJ5wBnmZ8Mlx0ig362RTOQWTjW1SqgzbDQ6Gzz4DsIyd+NdfeZo0Kq8TSxXmtszMFb5DQZvdRALrgTpKqWil1ACtdSrwNLAS2Ass1Frvyc92ZV5mUewefJDVt3fFKzUZ+vTBMznJ6IqEAQp6l7mP1jpIa+2pta6mtZ5lXb5ca11ba11Taz2hcEoVomh98dCznKh0E+zZQ99vZJTtksiUzV/kJ7MwQpK3Lx8OGA+entz727c027HW6JJEMTNlIMpPZmGUwzfVgbfeAuDJeW/C8eMGVySKkykDUQhDPfssO+q3wC8+Fvr2hTTpWFVSSCAKkZmbG9Mfe5mLfgHw66/w7rtGVySKiSkDUa4hCqNd8i/HDOso27z8MjUP/21sQaJYmDIQ5RqiMIMdDW+DZ5+F1FSGzXoF34R4o0sSRcyUgSiEabz1FjRtSuVzJ+gf9Z7R1YgiJoEoRG68vSEykkQvH1ptXCmjbLs4CUQhrqVuXeY89Kzl+ZAhcOiQsfWIImPKQJSbKsJsVt/elQ03t4a4OOjTB1JyH4hWOCdTTkOqtf4B+KFZs2aDjK5F5F9RTO2Z0zZtU47OXrQsx/UyT0ma330ALHiyJfT6Fpo0gU2bYNw4CO6ap88K52HKM0QhTCkgwDJUmJsbvPUWDfZtMboiUcgkEIXIjzvvhJdeAq15+vPXKHNZLuu4EglEIfLr5ZehZUvKXTpn6e8sE967DFMGotxUEaZmHWU73rcMt+78HT7+2OiKRCExZSBKTxVhetWr81nEKMvz55+H3btzX184BVMGohDOYH2zu/m1ZWdITLQ0xUlIMLokUUASiEIUwOcPPQe1alnOEEeNMrocUUASiEIUQJJPKYiMBE9P+Ogj+OEHo0sSBSCBKERBhYXBm29anj/+OJyQWXedlQSiEIXh+efhnnvg/Hl49FFITze6InEdTBmI0uxGOB03N5gzBypUgF9+gfdkqDBnZMpAlGY3wilVrgyff255PnYsbN5sbD0i30wZiEI4rc6dYdgwSE21NMWJizO6IpEPEohCFLa337aMivPvv/D000ZXI/JBAlGIwubjY2mK4+sLc+dy+6afjK5I5JEEohBFoV49mDwZgIFfvUuFc9IUxxlIIApRVAYNgu7dKZUYz7BZ42SUbScggShEUVEKPvuMcwEVqX14D4wfb3RF4hpMGYjSDlG4jHLl+OjxcaQrZenNsmaN0RWJXJgyEKUdonAle2vfzHcd+1kGkn3kETxiY40uSeTAlIEohKv5pnN/CA+H48ep8957Msq2SUkgClEM0tw94KuvwN+fCr//Dp9+anRJIhsSiEIUlxo1rk438Nxz8PffxtYjspBAFKI49enDqXvvtYyu3aePZbRtYRoSiEIUswPDhkFoKOzaBaNHG12OcCCBKEQxSytVynI90cMDPvwQli0zuiRhJYEohBFuvRUmTLA8f+wxbrh0zth6BCCBKIRxRoyAu++Gc+cY+vnrKBll23ASiEIYxc0N5s6F8uVpvG8znVdFGV1RiSeBKISRgoJg9mwA+iz+GLZuNbigks2UgSh9mUWJ0rUrK1r3xCPNOsr25ctGV1RimTIQpS+zKGm+7DGUo1VrwoED8MwzRpdTYpkyEIUoaVI8vZky4DXLaNtffAFRcj3RCBKIQpjE8So1YNIky4vBg+HIEUPrKYkkEIUwk8GD4f77ITYWHn7YMnufKDYSiEKYiVIwcyZUrQrr18NrrxldUYkigSiE2QQGwrx5lnCcMAHWrjW6ohJDAlEIM2rTBl58EdLTISICYmKMrqhEkEAUwqxefRVatIDoaMsMfjLKdpGTQBTCrDw9LaPi+PnBN99Yri2KIiWBKISZhYTAjBmW58OHw969xtbj4iQQhTC7iAjo2/fqKNtJSUZX5LIkEIVwBtOmWc4Wd+6EMWOMrsZlSSAK4Qz8/CAy0jLK9uTJ8OOPRlfkkiQQhXAWzZvD669bnvfrB6dOGVuPC5JAFMKZjBoFbdvC2bOWUJRRtguVBKIQzsQ2ynZgIPz0k+Xnsyg0xRaISqkQpdQspdSi4tqnEC6palWYNcvyfMwY2LbN2HpcSJ4CUSk1Wyl1Rim1O9PyDkqpf5RSB5VSud760lof0loPKEixQgir++6DIUMgJUVG2S5EeT1D/ALo4LhAKeUOTAM6AvWBPkqp+kqpRkqppZkeFQu1aiEEvPceNGwI+/dbGm2LAstTIGqt1wIXMi1uDhy0nvklA1HAfVrrv7TWXTI9zhRy3UIIX19LUxwfH8tEVQsXGl2R01M6jx3GlVLVgaVa64bW1z2BDlrrgdbXfYEWWuunc/h8IDABaA/M1FpPzGG9J4AnrC/rAP/k9csA5QFnmPHbWeoE56nVWeoE56nVWeqE/NcarLWukHmhRwEKUNksyzFdtdbngSevtVGt9afAp9dVkFJbtNbNruezxclZ6gTnqdVZ6gTnqdVZ6oTCq7Ugd5mjgRsdXlcDThSsHCGEME5BAnEzUEspVUMp5QX0BpYUTllCCFH88trsJhJYD9RRSkUrpQZorVOBp4GVwF5godZ6T9GVmifX9VPbAM5SJzhPrc5SJzhPrc5SJxRSrXm+qSKEEK5Ouu4JIYSV0wRidr1llFLllFI/K6UOWP8bYF2ulFIfWnvQ7FJK3WKCWl9VSh1XSu2wPjo5vPeitdZ/lFL3FmOdNyqlflVK7VVK7VFKDbcuN9VxzaVOMx5TH6XUJqXUTmut463LayilNlqP6QLrdXeUUt7W1wet71c3uM4vlFKHHY5pU+tyQ/9NWWtwV0ptV0ottb4u/GOqtXaKB9AKuAXY7bDsHWCM9fkY4G3r807Aj1iaBoUDG01Q66vAiGzWrQ/sBLyBGsC/gHsx1RkE3GJ97gfst9ZjquOaS51mPKYKKGN97glstB6rhUBv6/KPgaesz4cAH1uf9wYWGFznF0DPbNY39N+UtYbnga+wtIemKI6p05wh6ux7y9wHzLE+nwPc77B8rrbYANyglAoqnkpzrDUn9wFRWuskrfVh4CCWXkBFTmt9Umu9zfo8DsvNsaqY7LjmUmdOjDymWmtt61jsaX1ooC1gG9gk8zG1HetFQDulVHZtfIurzpwY+m9KKVUN6AzMtL5WFMExdZpAzEElrfVJsPyjAWx9pqsCxxzWiyb3f0DF5Wnrz43Ztp+hmKRW68+Km7GcKZj2uGaqE0x4TK0/7XYAZ4CfsZyhXtSWlhmZ67HXan3/EhBoRJ1aa9sxnWA9ppOUUt6Z67Qq7j/7ycAowDYAZCBFcEydPRBzkq9eNMVkBlATaAqcBN63Lje8VqVUGeAb4FmtdWxuq2azrNhqzaZOUx5TrXWa1ropls4KzYF6udRjWK2Z61RKNQReBOoCtwLlgNFG16mU6gKc0VpvdVycSz3XXauzB+Jp22m79b+2QSRM14tGa33a+hcwHfiMqz/hDK1VKeWJJWTma62/tS423XHNrk6zHlMbrfVFYA2Wa243KKVsXWUd67HXan2/LHm/3FLYdXawXp7QWusk4HPMcUxvB7oppY5gGUSmLZYzxkI/ps4eiEuAftbn/YDFDssftd4ZCwcu2X4CGiXT9ZYHANsd6CVAb+udsRpALWBTMdWkgFnAXq31Bw5vmeq45lSnSY9pBaXUDdbnvsDdWK55/gr0tK6W+ZjajnVPYLW23g0woM59Dv8jVFiuyTkeU0P+TWmtX9RaV9NaV8dyk2S11jqCojimxX2n6HofQCSWn0UpWP4PMADLdYFfgAPW/5bTV++gTcNy7eYvoJkJap1nrWWX9Q8syGH9sdZa/wE6FmOdd2D5KbEL2GF9dDLbcc2lTjMe08bAdmtNu4FXrMtDsITyQeBrwNu63Mf6+qD1/RCD61xtPaa7gS+5eifa0H9TDnW35upd5kI/ptJTRQghrJz9J7MQQhQaCUQhhLCSQBRCCCsJRCGEsJJAFEIIKwlEIYSwkkAUQggrCUQhhLD6PzslFG56NHOfAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "<table>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td/>\n",
+ "<th title=\"Variable name\">\n",
+ "Name\n",
+ "</th>\n",
+ "<th title=\"Value of parameter\">\n",
+ "Value\n",
+ "</th>\n",
+ "<th title=\"Hesse error\">\n",
+ "Hesse Error\n",
+ "</th>\n",
+ "<th title=\"Minos lower error\">\n",
+ "Minos Error-\n",
+ "</th>\n",
+ "<th title=\"Minos upper error\">\n",
+ "Minos Error+\n",
+ "</th>\n",
+ "<th title=\"Lower limit of the parameter\">\n",
+ "Limit-\n",
+ "</th>\n",
+ "<th title=\"Upper limit of the parameter\">\n",
+ "Limit+\n",
+ "</th>\n",
+ "<th title=\"Is the parameter fixed in the fit\">\n",
+ "Fixed\n",
+ "</th>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#FFFFFF;\">\n",
+ "<td>\n",
+ "0\n",
+ "</td>\n",
+ "<td>\n",
+ "loc\n",
+ "</td>\n",
+ "<td>\n",
+ "100.25\n",
+ "</td>\n",
+ "<td>\n",
+ "0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td>\n",
+ "1\n",
+ "</td>\n",
+ "<td>\n",
+ "scale\n",
+ "</td>\n",
+ "<td>\n",
+ "25.07\n",
+ "</td>\n",
+ "<td>\n",
+ "0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "</table>\n"
+ ],
+ "text/plain": [
+ "-------------------------------------------------------------------------------------------\n",
+ "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
+ "-------------------------------------------------------------------------------------------\n",
+ "| 0 | loc | 100.25 | 0.25 | | | | | |\n",
+ "| 1 | scale | 25.07 | 0.25 | | | | | |\n",
+ "-------------------------------------------------------------------------------------------"
+ ]
+ },
+ "execution_count": 43,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "m.migrad()\n",
+ "plt.figure(figsize=[5,5])\n",
+ "plt.yscale('log', nonposy='clip')\n",
+ "plt.ylim([0.1,1000])\n",
+ "blh.show(m)\n",
+ "m.get_param_states()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 44,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/mauro/anaconda3/envs/fits/lib/python3.7/site-packages/ipykernel_launcher.py:1: LogWarning: x is really small return 0\n",
+ " \"\"\"Entry point for launching an IPython kernel.\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "<table>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td/>\n",
+ "<th title=\"Variable name\">\n",
+ "Name\n",
+ "</th>\n",
+ "<th title=\"Value of parameter\">\n",
+ "Value\n",
+ "</th>\n",
+ "<th title=\"Hesse error\">\n",
+ "Hesse Error\n",
+ "</th>\n",
+ "<th title=\"Minos lower error\">\n",
+ "Minos Error-\n",
+ "</th>\n",
+ "<th title=\"Minos upper error\">\n",
+ "Minos Error+\n",
+ "</th>\n",
+ "<th title=\"Lower limit of the parameter\">\n",
+ "Limit-\n",
+ "</th>\n",
+ "<th title=\"Upper limit of the parameter\">\n",
+ "Limit+\n",
+ "</th>\n",
+ "<th title=\"Is the parameter fixed in the fit\">\n",
+ "Fixed\n",
+ "</th>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#FFFFFF;\">\n",
+ "<td>\n",
+ "0\n",
+ "</td>\n",
+ "<td>\n",
+ "loc\n",
+ "</td>\n",
+ "<td>\n",
+ " 100.25\n",
+ "</td>\n",
+ "<td>\n",
+ " 0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "-0.25\n",
+ "</td>\n",
+ "<td>\n",
+ " 0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "<tr style=\"background-color:#F4F4F4;\">\n",
+ "<td>\n",
+ "1\n",
+ "</td>\n",
+ "<td>\n",
+ "scale\n",
+ "</td>\n",
+ "<td>\n",
+ " 25.07\n",
+ "</td>\n",
+ "<td>\n",
+ " 0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "-0.25\n",
+ "</td>\n",
+ "<td>\n",
+ " 0.25\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "<td>\n",
+ "\n",
+ "</td>\n",
+ "</tr>\n",
+ "</table>\n"
+ ],
+ "text/plain": [
+ "-------------------------------------------------------------------------------------------\n",
+ "| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
+ "-------------------------------------------------------------------------------------------\n",
+ "| 0 | loc | 100.25 | 0.25 | -0.25 | 0.25 | | | |\n",
+ "| 1 | scale | 25.07 | 0.25 | -0.25 | 0.25 | | | |\n",
+ "-------------------------------------------------------------------------------------------"
+ ]
+ },
+ "execution_count": 44,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "m.minos()\n",
+ "m.get_param_states()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 2) same with gaussian pdf\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def line(x, a, b):\n",
+ " return a + x * b\n",
+ "\n",
+ "data_x = np.linspace(0, 1, 10)\n",
+ "# precomputed random numbers from a normal distribution\n",
+ "offsets = np.array([-0.49783783, -0.33041722, -1.71800806, 1.60229399, 1.36682387,\n",
+ " -1.15424221, -0.91425267, -0.03395604, -1.27611719, -0.7004073 ])\n",
+ "data_y = line(data_x, 1, 2) + 0.1 * offsets # generate some data points with random offsets\n",
+ "plt.plot(data_x, data_y, \"o\")\n",
+ "plt.xlim(-0.1, 1.1);\n",
+ "def least_squares(a, b):\n",
+ " yvar = 0.01\n",
+ " return sum((data_y - line(data_x, a, b)) ** 2 / yvar)\n",
+ "\n",
+ "m = Minuit(least_squares, pedantic=False)\n",
+ "m.migrad() # finds minimum of least_squares function\n",
+ "m.hesse() # computes errors \n",
+ "\n",
+ "# draw data and fitted line\n",
+ "plt.plot(data_x, data_y, \"o\")\n",
+ "plt.plot(data_x, line(data_x, *m.values.values()))\n",
+ "\n",
+ "# print parameter values and uncertainty estimates\n",
+ "for p in m.parameters:\n",
+ " print(\"{} = {} +/- {}\".format(p, m.values[p], m.errors[p]))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Minuit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# A complete instance of Minuit contains the starting point of each parameter (<name>) and \n",
+ "# the step size (<err_name>) for the dimension corresponding to that parameter.\n",
+ "#\n",
+ "# To correctly compute the uncertainties, it needs to know if you are doing a max likelihood \n",
+ "# or a least squares fit:\n",
+ "# errordef = 0.5 for negative log-likelihood functions\n",
+ "# errordef = 1 for least-squares functions\n",
+ "#\n",
+ "# You can specify the limits on the parameters to be fit (\"limit_varname\")\n",
+ "# lower limit: use limit_<name> = (<value>, None) or (<value>, float(\"infinity\"))\n",
+ "# upper limit: use limit_<name> = (None, <value>) or (-float(\"infinity\"), <value>)\n",
+ "# two-sided limit: use limit_<name> = (<min_value>, <max_value>)\n",
+ "#\n",
+ "# You can fix some of the parameters (ignore that dimension in the fit) setting <fix_name> = True/False\n",
+ "# fix_a=True,\n",
+ "\n",
+ "m = Minuit(least_squares, \n",
+ " a=5, b=5,\n",
+ " fix_a = True,\n",
+ " error_a=0.1, error_b=0.1,\n",
+ " limit_a=(0, None), limit_b=(0, 10),\n",
+ " errordef=1)\n",
+ "\n",
+ "m.get_param_states()\n",
+ "\n",
+ "\n",
+ "# Once Minuit is constructed you can still fix/release parameters as:\n",
+ "m.fixed[\"a\"] = False\n",
+ "m.fixed[\"b\"] = True\n",
+ "\n",
+ "# Trick to run over all parameters:\n",
+ "for key in m.fixed:\n",
+ " m.fixed[key] = False\n",
+ "m.migrad()\n",
+ "\n",
+ "# To change the value of a fixed parameter (or reset it to a different value) you can access it as:\n",
+ "m.values[\"a\"] = 0.5"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# MIGRAD"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Migrad performs Variable-Metric Minimization. It combines a steepest-descends algorithm along with \n",
+ "# line search strategy. Migrad is very popular in high energy physics field because of its robustness.\n",
+ "m.migrad()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pprint import pprint\n",
+ "\n",
+ "# To understand the results of the fit there are two dict-like objects:\n",
+ "# The first one is\n",
+ "pprint (m.get_fmin())\n",
+ "\n",
+ "# The most important one here is is_valid. If this is false, the fit did not converge and the result is useless !\n",
+ "#\n",
+ "# When is_valid = False it can be that the fit function is not analytical everywhere in the parameter space \n",
+ "# or does not have a local minimum (the minimum may be at infinity, the extremum may be a saddle point \n",
+ "# or maximum). Indicators for this are:\n",
+ "# is_above_max_edm=True\n",
+ "# hesse_failed=True\n",
+ "# has_posdef_covar=False\n",
+ "# has_made_posdef_covar=True\n",
+ "#\n",
+ "# Migrad reached the call limit before the convergence so that has_reached_call_limit=True. \n",
+ "# The used number of function calls is nfcn, and the call limit can be changed with the keyword argument ncall \n",
+ "# in the method Minuit.migrad\n",
+ "#\n",
+ "# Migrad detects converge by a small edm value, the estimated distance to minimum. This is the difference between\n",
+ "# the current minimum value of the minimized function and the next prediction based on a local quadratic \n",
+ "# approximation of the function (something that Migrad computes as part of its algorithm). If the fit did not \n",
+ "# converge, is_above_max_edm is true.\n",
+ "#\n",
+ "# To have a reliable uncertainty determination, you should make sure that:\n",
+ "# has_covariance = True\n",
+ "# has_accurate_covar = True\n",
+ "# has_posdef_covar = True\n",
+ "# has_made_posdef_covar = False \n",
+ "# hesse_failed = False\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The second one is a list of dict-like objects which contain information about the fitted parameters:\n",
+ "pprint(m.get_param_states())\n",
+ "\n",
+ "# Important fields are:\n",
+ "# index: parameter index.\n",
+ "# name: parameter name.\n",
+ "# value: value of the parameter at the minimum.\n",
+ "# error: uncertainty estimate for the parameter value."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Hesse"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The Hesse algorithm numerically computes the matrix of second derivatives at the function minimum \n",
+ "# (called the Hessian matrix) and inverts it.\n",
+ "# Pros:\n",
+ "# (Comparably) fast computation.\n",
+ "# Provides covariance matrix for error propagation.\n",
+ "# Cons:\n",
+ "# Wrong if function does not look like a hyperparabola around the minimum\n",
+ "\n",
+ "m.hesse()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Covariance - colored table\n",
+ "m.matrix()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# to access the matrix\n",
+ "cov = m.matrix()\n",
+ "print (cov)\n",
+ "\n",
+ "# or as a numpy array\n",
+ "npcov = m.np_matrix()\n",
+ "print (npcov)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Correlation - colored table\n",
+ "m.matrix(correlation=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# to access the matrix\n",
+ "corr = m.matrix(correlation=True)\n",
+ "print (corr)\n",
+ "\n",
+ "# or as a numpy array\n",
+ "npcov = m.np_matrix(correlation=True)\n",
+ "print (npcov)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Minos"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Minos implements the so-called profile likelihood method, where the neighborhood around the function minimum \n",
+ "# is scanned until the contour is found where the function increase by the value of errordef. \n",
+ "# Pros:\n",
+ "# Good for functions which are not very close to a hyper-parabola around the minimum\n",
+ "# Produces pretty confidence regions for scientific plots\n",
+ "# Cons:\n",
+ "# Takes really long time\n",
+ "# Result is difficult to error-propagate, since it cannot be described by a covariance matrix\n",
+ "#\n",
+ "# The results contain information as:\n",
+ "# At Limit: Whether Minos hit a parameter limit before the finishing the contour.\n",
+ "# Max FCN: Whether Minos reached the maximum number of allowed calls before finishing the contour.\n",
+ "# New Min: Whether Minos discovered a deeper local minimum in the neighborhood of the current one."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.minos()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.get_param_states()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Other ways to access the fits results\n",
+ "pprint(m.np_values())\n",
+ "pprint(m.np_errors())\n",
+ "pprint(m.np_merrors()) # The layout returned by Minuit.np_merrors() follows the convention \n",
+ " # [abs(delta_down), delta_up] that is used by matplotlib.pyplot.errorbar.\n",
+ "pprint(m.np_covariance())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Plot contours"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.draw_mncontour('a','b',nsigma=4) # nsigma=4 says: draw four contours from sigma=1 to 4"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.draw_profile('a');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.draw_mnprofile('a');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from __future__ import print_function\n",
+ "import numpy as np\n",
+ "%matplotlib inline\n",
+ "import matplotlib.pyplot as plt\n",
+ "from scipy.optimize import curve_fit, minimize, fsolve\n",
+ "from scipy.stats import norm, chi2, lognorm\n",
+ "import scipy.stats"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Generate data\n",
+ "# set the seed to always get the same samples\n",
+ "np.random.seed(seed=123456)\n",
+ "\n",
+ "# Generate a toy dataset on an exponential distribution (background)\n",
+ "# pdf = lambda * exp(-lambda * x) ; scale = 1/lambda\n",
+ "b = scipy.stats.expon.rvs(loc= 100, scale = 25, size=10000)\n",
+ "\n",
+ "# Generate a toy dataset on an gaussian distribution (signal)\n",
+ "#mu = 125, sigma = 1\n",
+ "s = scipy.stats.norm.rvs(loc=125, scale=1, size=100)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#Sum the two datasets\n",
+ "d = np.concatenate((s,b))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# histograma of the data\n",
+ "plt.figure(figsize=[15,5])\n",
+ "plt.subplot(131)\n",
+ "plt.hist(b,bins=80, range=[100,180])\n",
+ "plt.xlabel(r'x')\n",
+ "plt.ylabel(r'# entries / bin size = 2')\n",
+ "plt.subplot(132)\n",
+ "plt.hist(s,bins=80, range=[100,180])\n",
+ "plt.xlabel(r'x')\n",
+ "plt.ylabel(r'# entries / bin size = 2')\n",
+ "plt.subplot(133)\n",
+ "plt.hist(d,bins=80, range=[100,180])\n",
+ "plt.xlabel(r'x')\n",
+ "plt.ylabel(r'# entries / bin size = 2')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Build a pdf summing an exponential and gaussian component with a relative fraction f:\n",
+ "# f = #gauss / (#gauss+#exp) ; (1-f) = #exp / (#gauss+#exp)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def likelihood_point(x, f, b_slope, s_mu, s_sigma):\n",
+ " return f*(1./np.sqrt(2*np.pi*s_sigma**2)*np.exp(-0.5*((x-s_mu)/(s_sigma))**2.0)) + (1-f)*(1./b_slope)*np.exp(-x/b_slope)\n",
+ "\n",
+ "def nll(f, b_slope, s_mu, s_sigma):\n",
+ " return -np.sum([np.log(likelihood_point(x, f, b_slope, s_mu, s_sigma)) for x in d]) # VECTORIZE THIS !!!\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "f = 0.1\n",
+ "b_slope = 20\n",
+ "s_mu = 125\n",
+ "s_sigma = 1\n",
+ "\n",
+ "x = 125\n",
+ "print (likelihood_point(x, f, b_slope, s_mu, s_sigma))\n",
+ "\n",
+ "print (nll(f, b_slope, s_mu, s_sigma))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "f = 0.01\n",
+ "b_slope = 25\n",
+ "s_mu = 125\n",
+ "s_sigma = 1\n",
+ "\n",
+ "from iminuit import Minuit\n",
+ "\n",
+ "# def f(x, y, z):\n",
+ "# return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2\n",
+ "m = Minuit(nll, pedantic=False)\n",
+ "\n",
+ "m.migrad() # run optimiser\n",
+ "print(m.values) # {'x': 2,'y': 3,'z': 4}\n",
+ "\n",
+ "#m.hesse() # run covariance estimator\n",
+ "#print(m.errors) # {'x': 1,'y': 1,'z': 1}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "measurements = np.array([97.8621, 114.105, 87.7593, 93.2134, 86.6624, 87.4629, 79.7712, \\\n",
+ "91.5024, 87.7737, 89.6926, 133.506, 91.4124, 94.4401, 97.3968, \\\n",
+ "108.424, 103.197, 88.2166, 142.217, 89.0393, 102.438, 95.7987, \\\n",
+ "94.5177, 96.8171, 90.903, 132.463, 92.3394, 84.1451, 87.3447, \\\n",
+ "92.2861, 84.4213, 124.017, 90.4941, 95.7992, 92.3484, 95.9813, \\\n",
+ "88.0641, 101.002, 97.7268, 137.379, 96.213, 140.795, 99.9332, \\\n",
+ "130.087, 108.839, 90.0145, 100.313, 87.5952, 92.995, 114.457, \\\n",
+ "90.7526, 112.181, 117.857, 95.2804, 115.922, 117.043, 104.317, \\\n",
+ "126.728, 87.8592, 89.9614, 100.377, 107.38, 88.8426, 93.3224, \\\n",
+ "138.947, 102.288, 123.431, 114.334, 88.5134, 124.7, 87.7316, 84.7141, \\\n",
+ "91.1646, 87.891, 121.257, 92.9314])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1-D Maximum likelihood fit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We have a set of measurements which are distributed according to the sum of two Gaussians (e.g. this can be signal and background).\n",
+ "\n",
+ "$\\rho = \\frac{1}{3}\\frac{1}{\\sqrt{2\\pi \\sigma^2}} e^{-\\frac{1}{2}\\left(\\frac{x-p}{\\sigma}\\right)^2} + \\frac{2}{3}\\frac{1}{\\sqrt{2\\pi \\sigma_b^2}} e^{-\\frac{1}{2}\\left(\\frac{x-p_b}{\\sigma_b}\\right)^2}$ \n",
+ "\n",
+ "where for one of the two peaks the parameters are known already\n",
+ "\n",
+ "$p_b = 91.0$ \n",
+ "$\\sigma_b = 5.0$ \n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "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)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we assume the width of the peak we want to fit is already known: $\\sigma = 15.0$.\n",
+ "Perform a 1-D Maximum Likelihood fit for the position of the peak $p$.\n",
+ "\n",
+ "Complete the functions below which return the likelihood and negative log likelihood (NLL)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def likelihood_1d(params):\n",
+ " return np.prod([likelihood_point(x, params[0], 15.0) for x in measurements])\n",
+ "\n",
+ "def nll_1d(params):\n",
+ " return -np.log(likelihood_1d(params))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Minimize the NLL and give the best-fit result, including asymetric errors and plot the NLL."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "solution = minimize(nll_1d, [100.0], method='CG')\n",
+ "print ('Results:')\n",
+ "print (solution)\n",
+ "min_pos = solution.x[0]\n",
+ "min0 = solution.fun\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",
+ "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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# histogram of the data\n",
+ "plt.figure()\n",
+ "plt.hist(measurements,bins=100, range=[50,150])\n",
+ "plt.xlabel(r'x')\n",
+ "plt.ylabel(r'# entries / bin size = 2')\n",
+ "\n",
+ "# overlapped with the projection of the fit\n",
+ "x = np.linspace(measurements.min(), measurements.max(), 100)\n",
+ "plt.plot(x, len(measurements)*likelihood_point(x, min_pos, 15), \"r-\") # Plot of the data and the fit\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2-D Likelihood fit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we perform the 2-D Maximum Likelihood fit, fitting for both $\\sigma$ and $p$ at the same time."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def likelihood(params):\n",
+ " return np.prod([likelihood_point(x, params[0], params[1]) for x in measurements])\n",
+ "\n",
+ "def nll(params):\n",
+ " return -np.log(likelihood(params))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Minimize the NLL and find the best-fit result."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "solution = minimize(nll, [120.0, 10], method='CG')\n",
+ "print(\"position:\", solution.x[0], \"width:\", solution.x[1])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Create a 2D contour plot of the 1, 2 and 3 $\\sigma$ contours of the NLL and plot the best-fit solution."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scanA = np.linspace(110.0,130.0,50)\n",
+ "scanB = np.linspace(5,20,50)\n",
+ "minValue = nll(solution.x)\n",
+ "Z = [[nll([a,b]) - minValue for b in scanB] for a in scanA]\n",
+ "p1 = plt.contour(scanB, scanA, Z, [0.01,0.5, 2.0, 4.5])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compute numerically the error matrix of the NLL for the 2-D fit."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from scipy.misc import derivative\n",
+ "\n",
+ "# compute the error matrix\n",
+ "A = np.linalg.inv([\n",
+ " [\n",
+ " derivative(lambda x: nll([x, solution.x[1]]), solution.x[0], n=2),\n",
+ " derivative(lambda y: derivative(lambda x: nll([x, y]), solution.x[0]), solution.x[1])\n",
+ " ],\n",
+ " [\n",
+ " derivative(lambda x: derivative(lambda y: nll([x, y]), solution.x[1]), solution.x[0]),\n",
+ " derivative(lambda y: nll([solution.x[0], y]), solution.x[1], n=2)\n",
+ " ]\n",
+ "])\n",
+ "print(A, \"\\nsigma(position):\", np.sqrt(A[0,0]), \"sigma(width):\", np.sqrt(A[1,1]))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Binned ML fit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "With the same data as above, we now perform a binned ML fit and compare with the unbinned fit.\n",
+ "First, create a histogram of the data using np.histogram."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nBins = 10\n",
+ "histoMax = 170\n",
+ "histoMin = 70\n",
+ "binWidth = (histoMax-histoMin)/nBins\n",
+ "h0 = np.histogram(measurements, bins=nBins, range=(histoMin, histoMax))\n",
+ "print(h0[0])\n",
+ "print(h0[1])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compute the binned NLL:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def nll_binned(params):\n",
+ " # params is a list of [position, sigma]\n",
+ " expected = [likelihood_point(x+binWidth/2, params[0], params[1])*(binWidth/2)*sum(h0[0]) for x in h0[1]]\n",
+ " return sum([-np.log(expected[i]**h0[0][i]) for i in range(nBins)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Minimize the binned NLL:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "solution_binned=minimize(nll_binned, [120.0, 10], method='CG')\n",
+ "print(solution_binned)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Make a contour plot of the 1,2, and 3 $\\sigma$ contours for the binned NLL and overlay it with the unbinned contours."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scanA = np.linspace(110.0,130.0,50)\n",
+ "scanB = np.linspace(5,20,50)\n",
+ "Z_binned = [[nll_binned([a,b]) - solution_binned.fun for b in scanB] for a in scanA]\n",
+ "\n",
+ "fig1, ax2 = plt.subplots(constrained_layout=True)\n",
+ "\n",
+ "p1 = ax2.contour(scanB, scanA, Z, [0.01,0.5, 2.0, 4.5])\n",
+ "p2 = ax2.contour(scanB, scanA, Z_binned, [0.01,0.5, 2.0, 4.5], linestyles=\"dotted\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat the same for 50 bins:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scanA = np.linspace(110.0,130.0,50)\n",
+ "scanB = np.linspace(5,20,50)\n",
+ "Z_binned = [[nll_binned([a,b]) - solution_binned.fun for b in scanB] for a in scanA]\n",
+ "\n",
+ "fig1, ax2 = plt.subplots(constrained_layout=True)\n",
+ "\n",
+ "p1 = ax2.contour(scanB, scanA, Z, [0.01,0.5, 2.0, 4.5])\n",
+ "p2 = ax2.contour(scanB, scanA, Z_binned, [0.01,0.5, 2.0, 4.5], linestyles=\"dotted\")\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.7"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
--
GitLab