{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Example of maximum-likelihood fit with iminuit.\n", "# pdf is a mixture of Gaussian (signal) and exponential (background),\n", "# both truncated in [xMin,xMax].\n", "# Here parameters other than signal fraction theta\n", "# and exponential mean xi are treated as constant.\n", "# G. Cowan / RHUL Physics / November 2020\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": [ "# define pdf and generate data\n", "np.random.seed(seed=1234567) # fix random seed\n", "theta = 0.2 # fraction of signal\n", "mu = 10. # mean of Gaussian\n", "sigma = 2. # std. dev. of Gaussian\n", "xi = 5. # mean of exponential\n", "xMin = 0.\n", "xMax = 20.\n", "\n", "def f(x, theta, xi):\n", " fs = stats.truncnorm.pdf(x, a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n", " fb = stats.truncexpon.pdf(x, b=(xMax-xMin)/xi, loc=xMin, scale=xi)\n", " return theta*fs + (1-theta)*fb" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "numVal = 200\n", "xData = np.empty([numVal])\n", "for i in range (numVal):\n", " r = np.random.uniform();\n", " if r < theta:\n", " xData[i] = stats.truncnorm.rvs(a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n", " else:\n", " xData[i] = stats.truncexpon.rvs(b=(xMax-xMin)/xi, loc=xMin, scale=xi)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Function to be minimized is negative log-likelihood\n", "def negLogL(theta, xi):\n", " pdf = f(xData, theta, xi)\n", " return -np.sum(np.log(pdf))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Create instance of 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 (needed here to keep pdf>0).\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", "\n", "m = Minuit(negLogL, theta=0.5, xi=5., error_theta=0.1, error_xi=0.2, limit_theta=(0., 1.), \n", " limit_xi=(0, None), fix_xi=False, errordef=0.5, pedantic=False)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# do the fit, get errors, extract result as arrays\n", "m.migrad() # minimize -logL\n", "MLE = m.np_values() # parameter estimates\n", "sigmaMLE = m.np_errors() # their standard deviations\n", "cov = m.np_matrix() # covariance matrix\n", "rho = m.np_matrix(correlation=True) # correlation coeffs." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "parameter index, estimate and standard deviation:\n", "0 0.2046388058666193 0.052734483201438365\n", "1 5.1085667979611316 0.6447893594778562\n", "\n", "parameter indices, covariance and correlation coeff.:\n", "0 0 0.0027969098228253996 1.0\n", "0 1 -0.018132494205288098 -0.5316588752771374\n", "1 0 -0.018132494205288098 -0.5316588752771374\n", "1 1 0.41588235455534234 1.0\n" ] } ], "source": [ "# rename for convenience below\n", "thetaHat = MLE[0]\n", "sigmaThetaHat = sigmaMLE[0]\n", "xiHat = MLE[1]\n", "sigmaXiHat = sigmaMLE[1]\n", "\n", "# line below gives: AttributeError: 'iminuit._libiminuit.Minuit' object has no attribute 'nfit'\n", "#npar = m.nfit\n", "npar = len(m.np_values())\n", "print(r'parameter index, estimate and standard deviation:')\n", "for i in range(npar):\n", " print(i, MLE[i], sigmaMLE[i])\n", "\n", "print()\n", "print(r'parameter indices, covariance and correlation coeff.:')\n", "nfreepar = len(cov[0])\n", "for i in range(nfreepar):\n", " for j in range(nfreepar):\n", " print(i, j, cov[i,j], rho[i,j])" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAF7CAYAAAAE6luLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5xU1f3/8dcHkN6bgEpTVAQpumIQUSxoYhdiS1E0apTEaEw1+osxiTFGQyTFEMtXbInEiIliQ7FGUNxVxIoNsFBcVJC21M/vjzMbxmFmd2Z3ds/M7vv5eNzH7Nx77rmfu3eH+XDuueeYuyMiIiJS35rEDkBEREQaJyUhIiIiEoWSEBEREYlCSYiIiIhEoSREREREolASIiIiIlE0ix1Asenatav37ds3dhgiIiL1pqysbIW7d8t3vUpCctS3b19KS0tjhyEiIlJvzGxxXdSr2zEiIiIShZIQERERiUJJiIiIiEShJERERESiUBIiIiIiUSgJERERkSiUhIiIiEgUSkJEREQkCiUhIiIiEoWSEBEREYlCSYiIiIhEoSREREREolASIiIiIlEoCREREZEolISIiIhIFEpCREREJAolISIiIhKFkhARERGJQkmIiIiIRKEkRERERKKInoSY2UQzW2hmFWZWZmajqyjb0symmtl8M9tkZk+mKTPVzDzNsjapzIQMZVrW0WmKiIhIiqhJiJmdAkwGfgMMB2YDD5lZ7wy7NAUqgD8DD2QocyHQM2V5D/hnSrl1qeXcvaLGJyMiIiI5aRb5+BcDU939xsT7C8zsy8D5wCWphd19LXAegJkNATqmKbMKWFX53sxGAf2Bb25f1Jfl4yREREQkd9FaQsysObAvMDNl00zggDwe6hzgNXefnbK+lZktNrMPzWyGmQ3P4zFFRESkGjFvx3Ql3F5ZnrJ+OdAjHwcwsw7AScCNKZsWAGcBxwOnEW7xPGtmAzLUc66ZlZpZaXl5eT5CExERafSid0wFPOW9pVlXU98gJDq3f+GA7nPc/VZ3n+fuzwCnAO8CF6QN0P0Gdy9x95Ju3brlKTQREZHGLWYSsgLYwvatHt3ZvnWkps4B7nH3T6sq5O5bgFIgbUuIiIiI5F+0JMTdNwJlwNiUTWMJT8nUipmNAIay/a2YdGUNGAIsre1xRUREJDuxn46ZBNxuZnOBZwlPvvQCpgCY2VXACHc/rHIHM9sLaE7oU9LWzIYBuPu8lLrPBd4Gnko9qJldDjyX2N4e+B4hCTk/nycnIiIimUVNQtx9mpl1AS4jjNXxKnCUuy9OFOkJ7Jqy24NAn6T3LyVerXKFmbUDTgV+6e7p+pd0BG4g3ApalajjIHefW7szEhERkWxZ+u9oyaSkpMRLS0tjhyEiIlJvzKzM3UvyXW8hPB0jIiIijZCSEBEREYlCSYiIiIhEoSREREREolASIiIiIlEoCREREZEolISIiIhIFEpCREREJAolISIiIhKFkhARERGJQkmIiIiIRKEkRERERKJQEiIiIiJRKAkRERGRKJSEiIiISBRKQkRERCQKJSEiIiIShZIQERERiUJJiIiIiEShJERERESiUBKSo4UrwT12FCIiIsVPSUiO1myEsqWxoxARESl+SkJy1MRg6suxoxARESl+SkJy1LkVPPg2LFsTOxIREZHipiQkR11awVaHO16JHYmIiEhxUxKSo+ZN4fD+8PdXoGJz7GhERESKl5KQGpgwFD5ZDw+8HTsSERGR4qUkpAZG7QK7dYZb5ulxXRERkZpSElIDZqE15JWP9biuiIhITSkJqaFxe0L75nCrHtcVERGpESUhNdSmOZw0CB58B5brcV0REZGcKQmphTOGwJatcNv82JGIiIgUn+hJiJlNNLOFZlZhZmVmNrqKsi3NbKqZzTezTWb2ZJoyY8zM0yx7ppQbb2avm9mGxOuJucbepyMcsSvc+Qqs35Tr3iIiIo1b1CTEzE4BJgO/AYYDs4GHzKx3hl2aAhXAn4EHqql+ENAzafnfA7VmNhKYBtwJDEu83m1m++d6DmcPh88qYPqbue4pIiLSuMVuCbkYmOruN7r7G+5+AbAUOD9dYXdf6+7nufsNwIfV1P2xuy9LWrYkbbsIeMLdr0wc90rgycT6nOzXC4Z0h5tfCiOpioiISHaiJSFm1hzYF5iZsmkmcEAeDlFqZkvNbJaZHZKybWSa4z5Sk+Oawdn7wLufwROLahipiIhIIxSzJaQr4fbK8pT1y4Eetai3siVlPDAOWADMMrODksr0yOW4ZnaumZWaWWl5efl224/aDXq2hZterEXUIiIijUzs2zEAqTcxLM267CtzX+DuU9y9zN3nuPtE4GHghzU9rrvf4O4l7l7SrVu37bbv0DQMXjb7Q3ht+xxFRERE0oiZhKwAtrB960N3tm+lqK3ngQFJ75fl+7inDYbWO4S+ISIiIlK9aEmIu28EyoCxKZvGEp6SyadhhNs0lebk+7gdWsLJe8F9CzR4mYiISDZi346ZBEwws7PNbKCZTQZ6AVMAzOwqM5uVvIOZ7WVmwwh9Stqa2bDE+8rtF5nZCWY2wMwGmdlVwAmEx3orTQYONbNLzGxPM7sEOAS4rjYnc9Yw2LwVbtXgZSIiItVqFvPg7j7NzLoAlxHG8ngVOMrdFyeK9AR2TdntQaBP0vvKGyCWeG0OXAvsBKwHXgOOdvcHk44728xOBX4NXAG8C5zi7s/X5nwqBy+7Yz58d79we0ZERETSM9dc9DkpKSnx0tLSjNtLl8D4u+EXB8OZwzIWExERKRpmVubuJfmuN/btmAanpFcYwOzGF2HTlurLi4iINFZKQurAefvCR6thxtvVlxUREWmslITUgUP7wYDO8Lcy0N0uERGR9JSE1IEmFlpD3lgBTy2uvryIiEhjpCSkjhy3RxjK/a9lsSMREREpTEpC6kjzpvCt4fDchzBvWexoRERECo+SkDp02mBo3wKmqDVERERkO0pC6lDb5vDNIfDwO7Dws9jRiIiIFBYlIXXszKHh1oxaQ0RERL5ISUgd69YGThkE97wBS1bHjkZERKRwKAmpB+ftC04YN0REREQCJSH1YKf2cOKe8I9XoXxt7GhEREQKg5KQejKxBDZthZtfqr6siIhIY6AkpJ707wTHDIDb5sPKitjRiIiIxKckpB59Zz9YuwmmzosdiYiISHxKQurRnl1hbH/4v3mwZmPsaEREROJSElLPvrsfrNoAd8yPHYmIiEhcSkLq2bAeMLo33PgirN8UOxoREZF4lIRE8L0RsGI93PFK7EhERETiURISwYidYNQuMKVUrSEiItJ4KQmJ5KL91RoiIiKNm5KQSJJbQ9apNURERBohJSERfV+tISIi0ogpCYlov53gwF3gb2oNERGRRkhJSGQXfSm0htyucUNERKSRURIS2X69wrghfytTa4iIiDQuSkIKwIX7wyfr4baXY0ciIiJSf5SEFID9esHBfeCvZbB6Q+xoRERE6oeSkALxo5GwsgJueil2JCIiIvVDSUiB2HtH+MpuIQn5dH3saEREROqekpACcvGXYO1G+Gtp7EhERETqnpKQArJ7Fxg3EG59GZaviR2NiIhI3YqehJjZRDNbaGYVZlZmZqOrKNvSzKaa2Xwz22RmT6YpM87MZppZuZmtNrPnzey4lDITzMzTLC3r4BRzctH+sMXhj3NjRyIiIlK3oiYhZnYKMBn4DTAcmA08ZGa9M+zSFKgA/gw8kKHMwcDjwNGJOh8E7k2T3KwDeiYv7l5R87PJj94d4NRBcNdr8P6q2NGIiIjUndgtIRcDU939Rnd/w90vAJYC56cr7O5r3f08d78B+DBDmQvd/bfuPtfd33H3K4Ay4ITti/qy5CWP51Ur3xsBTQ2uez52JCIiInUnWhJiZs2BfYGZKZtmAgfk+XDtgM9S1rUys8Vm9qGZzTCz4Xk+Zo3t2BbOGArT34AFK2JHIyIiUjditoR0JdxeWZ6yfjnQI18HMbPvADsDtyetXgCcBRwPnEa4xfOsmQ3IUMe5ZlZqZqXl5eX5Cq1KE0ugXXP47bP1cjgREZF6F/t2DICnvLc062rEzMYD1wBfd/fF/zug+xx3v9Xd57n7M8ApwLvABWkDdL/B3UvcvaRbt275CK1anVrBxP3g8UXwXNobTyIiIsUtZhKyAtjC9q0e3dm+dSRniQTkduB0d7+vqrLuvgUoBdK2hMRy5jDo2Ta0hnhe0jIREZHCES0JcfeNhA6jY1M2jSU8JVNjZnYycAcwwd3/lUV5A4YQOsUWjJbNwgBmLy2DB9+JHY2IiEh+xb4dMwmYYGZnm9lAM5sM9AKmAJjZVWY2K3kHM9vLzIYR+pS0NbNhifeV208F7gR+CjxtZj0SS+ekMpeb2ZFm1j+x782EJGRKHZ9vzsYPhD26wDWzYdOW2NGIiIjkT7OYB3f3aWbWBbiMMFbHq8BRSf03egK7puz2INAn6X3llG+WeD2PcF7XJZZKTwFjEj93BG4g3ApalajjIHcvuCHCmjaBn4yCs+6Df7wKpw+NHZGIiEh+mKuzQU5KSkq8tLR+J3dxh1PvgXc+hacmQNvm9Xp4ERFp5MyszN1L8l1v7NsxkgUzuORAWLEebiiLHY2IiEh+KAkpEsN6wLG7w99ehKWrY0cjIiJSe0pCishPRoVbM9fMiR2JiIhI7SkJKSK7tIdvDYd73oD5tR5JRUREJC4lIUVmYgl0bQW/eloDmImISHFTElJk2rWAH4yEuUvgIQ1gJiIiRUxJSBE6eVAYwOyqZ2HD5tjRiIiI1IySkCLUrAlcNhreXwVTX44djYiISM0oCSlSB/WBQ/rCn+bCinWxoxEREcmdkpAidtloWL8Zfler6f5ERETiUBJSxHbrDGcOg3++Bi8vix2NiIhIbpSEFLkLR0DX1nD5U7BVj+yKiEgRURJS5Nq1gJ+OgpeWwfQ3YkcjIiKSPSUhDcC4gTC8B/z2Wfh8Q+xoREREsqMkpAFoYnDFweEpmT/OjR2NiIhIdpSENBBDe4RBzG6ZB29/GjsaERGR6ikJaUB+fAC0bgaXP6l5ZUREpPApCWlAuraGHx0Az34A970VOxoREZGq5ZyEmFkLM+tnZnuZWbe6CEpq7ut7w5DuYZZddVIVEZFCllUSYmbtzOx8M3saWAW8A7wKLDOzD8zsRjPbry4Dlew0bQJXHho6qf5+TuxoREREMqs2CTGz7wOLgLOAR4HjgWHA7sBI4HKgGfComT1sZgPqLFrJypAd4ZtD4Lb58MrHsaMRERFJr1kWZQ4ADnb3VzNsnwv8n5mdB3wLOBh4O0/xSQ398AB46B247HG495TwGK+IiEghqbYlxN1PqiIBSS63wd2vd/eb8hOa1EaHFnDpaJi3HP5R7dUTERGpf9ncjvmvmXWuj2Akv07YA0buHEZSLV8bOxoREZEvyqZj6gHAzskrzGyomW3XwG9mbfMVmNSeWeikWrEZrng6djQiIiJflO0jurtV/pBIPsqAwckFzOxuYJWZvWZmg/IXotTGrp3gu/vB/W/BrIWxoxEREdkm2yTkq0k/75LYb8fKFWbWERgHHAncDdycrwCl9s4vgd27hE6qazbGjkZERCTINgn5kpldaGbNgbOBCsJTMJV6ARXu/hhwNfDX/IYptdG8Kfz2MFi6RmOHiIhI4cgmCbkFOAO4GFgLXApcCJxrZnsmyhwDvAfg7uvd/dY6iFVqYd+eYeyQW+bBvGWxoxEREclinBB3/xaAme0K7A2sdPeFZtYemG9mrxH6h1xWp5FKrf34AJj5HvzkMZhxGuzQNHZEIiLSmGU9d4y7b3b3l9x9YeL974GDgH8D3wZ+VzchSr60awG/PgTe/ASmlMWORkREGrtsRkzNyN2fA57LUyxSD8b2h6MHwB/nwpG7hg6rIiIiMeQ8i26+mdlEM1toZhVmVmZmo6so29LMpprZfDPbZGZPZih3cKKuCjN7LzGkfGqZ8Wb2upltSLyemMfTKmi/HANtm8MPH4XNW2NHIyIijVXUJMTMTgEmA78BhgOzgYfMrHeGXZoSnsz5M/BAhjr7AQ8m6hoOXAX8yczGJ5UZCUwD7iRMxncncLeZ7Z+H0yp4XVvDr8bAy8vhphdjRyMiIo2VuXt+KgqJw4funvX/rc3seWC+u5+TtO5t4F/ufkk1+/4ZGOzuY1LWXw2Mc/cBSetuAga5+8jE+2lAZ3cfm1TmMaDc3U+r6rglJSVeWlqa7SkWLHc47wF4YhE88DUYoIH5RUQkAzMrc/eSfNebz5aQRcDLZnZQNoUTY47sC8xM2TSTMFR8TY1MU+cjQImZ7VBNmdoct6iYhU6qrXaAHz0KW3RbRkRE6lk+k5CzgOnANVmW70q4vbI8Zf1yoEct4uiRoc5miWNWVSbtcc3sXDMrNbPS8vLyWoRWWLq1gSsOhpeWwc3zYkcjIiKNTd6SEHef6u6Xu3uu/SpS7wdZmnU5h5OmztT1WR/X3W9w9xJ3L+nWrVstQyssx+8BR/SHa2fD25/GjkZERBqTapMQM+tiZjeb2TIz22xmn5jZbDO7xsxG1OLYK4AtbN/60J3tWylysSxDnZuBT6opU5vjFqXKmXZb7wDffwQ2bYkdkYiINBbZtITcARwIXAl8A/gRsA9h7pj/mtmTZtY/1wO7+0bCbLxjUzaNJTzZUlNzgMPT1Fnq7puSyuT7uEWrexv4zaHwysfw5xdiRyMiIo1FNknIwcB4d/+Tu9/l7v8HbAJOJUxcVwbMNrPda3D8ScAEMzvbzAaa2eREnVMAzOwqM5uVvIOZ7WVmwwj9O9qa2bDE+0pTgJ3N7LpEnWcDE4Brk8pMBg41s0vMbE8zuwQ4BLiuBufQIBw1AMbtCX+aq7llRESkfmQzYupHQNpxNd19BfADM1tK+GL/Si4Hd/dpZtaFMO9MT+BV4Ch3X5wo0hPYNWW3B4E+Se9fSrxaos6FZnYU8AfgfGAJ8D13vyfpuLPN7FTg18AVwLvAKe7+fC7xNzRXjIHnPgy3ZR78WnhyRkREpK5UO06ImX2XMIPuqe4+N7FuNTDU3d9LvO8DvObubes43ugayjghmTz7AXxtOkwYGpISERGRaOOEuPufgX8Ac8xsViIpSd3vm0DDeXa1ERu1C3xrGEx9GZ5ZXH15ERGRmsrqEV13vxTYn5Bo/BZoBbxmZh+Y2Urgp8AP6ixKqVc/HhVGUL34Ufh0fexoRESkocp6nBB3L3X3U4GOwAjgbMK8LOcCfd19et2EKPWtZTP445dhZUUYTTVPI/uLiIh8QTYdU7/A3TcDpYlFGqi9usElo+CKp+H2+XD60NgRiYhIQxN1Fl0pbGcOgzF94NfPwIIVsaMREZGGRkmIZGQG146Fds3hgoehYnPsiEREpCFREiJV6tYGfn8ELPgErnwmdjQiItKQ5JyEmFlvM+uZYb2SmgZoTN/w2O5t8+GRd2NHIyIiDUVNkoZFwKwM6182s4NqE5AUpp+Mgr27ww8fhfdXxY5GREQagpokIWcBP8uwfjpwTa0ikoLUohlcfxTg8N2HYKNm2xURkVrKKQkxs27uPtXd/526LbH+cnffP3/hSSHp3QF+NxZeXg6//W/saEREpNjl2hIy28z610kkUhS+slt4dPfmeeofIiIitZNrEvIgIRHZJ3mlmR1kZs/mLywpZJeMgiHqHyIiIrWUUxLi7hcC1wJPmNkRZjbMzB4GngDer4sApfC0aAZ/SfQPmfigxg8REZGaybljqrtfC/wGmAG8AKwGhrj7aXmOTQpY7w7whyPhlY/h50/GjkZERIpRrh1TdzGzvwG/JCQgG4AH3P21ughOCtvh/eGCETDtNfjHq7GjERGRYpNrS8jbwHDgGHcfBRwH/MHMLs17ZFIUvr8/HNwntIa8vCx2NCIiUkxyTUK+4e4j3P1RAHd/HBgDnG9m1+c7OCl8TZvA5COhe2s47wH4dH3siEREpFjk2jH1X2nWvQyMIiQj0gh1agVTjoFP1sMFD8HmrbEjEhGRYpCXuV7cfTEhEZFGau/u8OtD4L8fwFUayExERLJQbRJiZv2yqcjdP7Ngl9qHJcXo5EFhILObXoK7X48djYiIFLpsWkLmmNnNZjYyUwEz62Rm5wOvA8fnLTopOpeNhgN3gZ89DmVLY0cjIiKFrFkWZf4NrAMeMLMtQBmwFKgAOgF7AQOBucBF7v5IHcUqRaBZkzCQ2XF3wbdnwP2nQs92saMSEZFClE1LyFnA1cDOQBdgCdAR6AdsBm4Fhrv7KCUgAtCxJdx0LKzfDOfM0IiqIiKSXjYtIR8A+7v7fWYG8FN3/7huw5Jit3uX8Oju2feHOWb++GVoYrGjEhGRQpJNS8hvgXvM7EXAgbPMbLSZta/b0KTYHd4ffjoK7n8Lfj8ndjQiIlJoqm0JcfcbzexpQofTYcAEwrDtTc1sMTCvcnH3++owVilC394XFq2CP78AfTqEJ2hEREQgy3FC3H2Bu/+OMGz7gUA7YARwJfARMBa4ra6ClOJlBr8aA6N7wyWPw7MfxI5IREQKhbl77BiKSklJiZeWlsYOo+h8vgHG/ROWr4V7T4bdOseOSEREsmVmZe5eku968zJiqkh12reAW46H5k1gwn+gfG3siEREJDYlIVJvdmkPNx8HK9bBhPtgzcbYEYmISExKQqReDesB1x8Fb5SHWXc3bokdkYiIxBI9CTGziWa20MwqzKzMzEZXU35vM3vKzNab2Udm9nNLDGCS2D7VzDzNsjapzIQMZVrW5blKcGg/uPpweOb9MIbIVnVLEhFplLIZrKzOmNkpwGRgIvDfxOtDZraXu7+fpnx74FHgaWA/YA9gKrAW+H2i2IXAT1N2fTaxT7J1wK7JK9y9ohanIzk4aS/4eC38bjZ0bxPmnBERkcYlahICXAxMdfcbE+8vMLMvA+cDl6Qp/3WgNXCGu68HXjWzgcDFZjbJg1XAqsodzGwU0B/4Zkpd7u7L8nw+koOJJSERufFF6Noazts3dkQiIlKfot2OMbPmwL7AzJRNM4EDMuw2EngmkYBUegToBfTNsM85wGvuPjtlfSszW2xmH5rZDDMbXkWs55pZqZmVlpeXZyomOTKDnx8Ex+4OV/0X7nwldkQiIlKfYvYJ6Qo0BZanrF8O9MiwT48M5Su3fYGZdQBOAm5M2bSAMDHf8cBphBmBnzWzAekO6u43uHuJu5d069YtQ2hSE02bwKQj4NC+cOnjMP3N2BGJiEh9id4xlTAfTTJLs6668unWA3yDkOjc/oUK3Oe4+63uPs/dnwFOAd4FLsg6asmb5k3hr0fDl3aGH86Eh9+JHZGIiNSHmEnICmAL27dgdGf71o5KyzKUJ8M+5wD3uPunVQXi7luAUiBtS4jUvZbN4KZjYciOcMHD8PTi2BGJiEhdi5aEuPtGoIww70yysUBq/41Kc4DRKY/SjgWWAIuSC5rZCGAo29+K2U7iEd8hwNJsYpe60bY53Ho87NYJzpkBcz6MHZGIiNSl2LdjJgETzOxsMxtoZpMJnUynAJjZVWY2K6n83wmP1k41s8FmNo7wOO4k334SnHMJE+49lXpQM7vczI40s/5mNgy4mZCETMn3CUpuOrSE208Mo6tO+A/M1oR3IiINVtQkxN2nARcBlwHzCDP0HuXulY3xPUkayyPx+O1YQqJSCvyFMD7IpOR6zawdcCpwU5rkBKAjcAPwBuFpnJ2Ag9x9bt5OTmqsa2u4azz07gBn3qeZd0VEGirNopsjzaJbf1asg69Nh0Ur4f+OgwN7x45IRKRx0iy60uh0bQ3/GAf9OsJZ98Ez6qwqItKgKAmRgtalNfx9HPTvBGfdD4+8GzsiERHJFyUhUvC6JPqIDOoG5z8A09+IHZGIiOSDkhApCh1bwp0nwv47wfdnwq0vx45IRERqS0mIFI02zeGW4+GI/vDzJ+FPc0H9qkVEipeSECkqLZvB9UfBiXvCtXNCMrJla+yoRESkJprFDkAkVzs0DZPedW8DfyuDJavhT1+B1jvEjkxERHKhlhApSk0MfnYg/HIMPL4ITr0HytdGDkpERHKiJESK2hlD4YZjYMEncOI/4Z0qpyoUEZFCoiREit7Y/jBtPKzfBCdOgycWxY5IRESyoT4h0iAM6wH/OTXMvnvmf+Ano+C8fcEsdmSN28oKeL08LG9+EobiX1kRls8qYOOWMHtyu8qlBezaKYwJM7g7DOgc+gCJSMOkJEQajJ3bw/ST4EePwW+fDV98vzscWqnDar0pXwtPvw9PLoKypfDR6m3burWGHm3DmC+7tA8zJjdvAms2wZqNYVlZAf98HdZtCvs0bwpDusPRA+CoAWF/EWk4NIFdjjSBXeFzh7+Wwu9mw8BuMOUo6NMxdlQN19ufwn0Lwm2wVz4O67q1hpE7w6DuMKhruA5dW2dX35atsHAlvFYOr34M/30fXl8Rtu3XKyQk4/YMSYyI1I+6msBOSUiOlIQUj8cXwoWPhKTkt4fBMbvHjqjh+Hgt3P8WTH8zJApNDPbtCWP6wJi+sFe3sC5f3v0MHngbHngr3NZp2xy+uTecvU/2yY2I1JySkAKhJKS4fPg5fPcheGkZfH0w/PzgMOCZ5G6rw9OL4fb54bHorQ57dw+tEsfuDt3a1E8cr5WHlq4Zb4XbNacNDv1/erarn+OLNEZKQgqEkpDis2lLGF11Shns2QX+fFTo8CjZqeynccd8WLwq3Go5eS84cWDc3+O7n8FfX4B7F0CzJnDhCDhnH3VkFakLSkIKhJKQ4vXkojD53dqN8IORcPZwaKqH1DN67zO4+SX41xtQsTn0xzh9CHx5t9ACUSg++Bx+9TQ88i7s3gV+cwjst1PsqEQaFiUhBUJJSHH7eC1c+jjMfA/26QnXjg2PhErgDnOXwI0vwmPvhWTjhD3hzKGhc2khe+y9MJfQR6vhlEFw6YHqvCqSL0pCCoSSkOLnDv9eAJc/Gf6H/8MD4KxhoUm/sdqyFR5+N8zF8/Jy6NwqdPz85pD66+uRD+s2wXXPw00vhj4i138FhvaIHZVI8VMSUiCUhDQcy9fCz2bBYwthYFe4Ygzs38ia8Ss2w92vh5aPxaugX8dwm+qrexV3B96XlsF3Hky0fIuq9D8AAB4mSURBVI2GCUM1cJ1IbSgJKRBKQhoWd3jwHbjymdCMf9zu4UuroQ+KtWJdeMrltvnw6XoYtiOcVwJH9G84/WRWVsAPH4VH34Mjd4VrxkKHFrGjEilOSkIKhJKQhmn9Jri+NNyOaNokPPL5reFhPIqG5O1Pw62Ke9+EDVvgsH7w7X1hRK+G2VLgDje9FEbQ3akd3HK8+gCJ1ISSkAKhJKRhe38V/PqZ8KRFp5bhC/qModC6iId+37I1jGZ623x4ajG0aApfHQjf2qfxfCGXLoFzZ8DmrXDTsTCikd12E6ktJSEFQklI4/DyMpj0HDy5GLq2CrcqThtcXC0jn66Haa/BHa+EQdt2bANfT3Q27dwqdnT17/1VMOE/4ZHe34+F4/aIHZFI8VASUiCUhDQupUvgD8/Bfz8ICcj4gWGsjN0KdLCzTVtCa8c9b4QOtxu3hDlcTh8CY/trIK+VFaFF5PmP4CcHwPklDfM2lEi+KQkpEEpCGqeXl8Gt88N8KRu3wIG7hJaRQ/vFv1Wz1WHeMrjvrTCR3CfroUsrOH6PEOPuXeLGV2g2bA4zLf9nQXhq5hcHKxERqU5dJSFF/BCeSP0Z2gMm9QgDYN31WhjC/DsPhcdYD+kbZnY9tC+0qafbNes3hdaZR98LE/WVrwsDi43tD+P3hIP6qNUjkxbN4Lojw8R3N78UksorD83vhHsikh0lISI56NIavrNfeHrmhSVhZteH3glLi6YwdMcwZPiIXmFE1vZ5eiT0k3Xw4rJwe6hsKcxfHp5uadc8zFp7eD84pJ8eQc1WE4P/NxpaNoW/lIbbWFcf3nAeTxYpFrodkyPdjpFUW7ZC6dLQKjH3ozC1/RYPX3T9OkKfDonXjtC7Q0gcWu+wbXGHNZvCnDZrNsLnG8LAYYtXwsKVsGglLFkTjrVDExjcHfbtGVpgRuxUWPO4FBt3+OPc0An5+D1g0hGNe+RckUx0O0akQDVtEkZarRxtde3GMGLnC0tgwSchmXjuozCkeC46tYS+HWH/nWGPLlDSE/besbhHMi00ZnDh/iG5u3p2aBH501eUiIjUF/1zJpJnbZrDgb3DUskdPl4HH6yCtZtCQrJ+U/jZEvu0bQ5tdoB2LaB3e02+Vp8m7hdalH71DLR4NLSIqI+ISN2LnoSY2UTgR0BP4DXgInd/poryewN/BkYAnwJ/A37liftKZjYGeCLNrgPd/c2kesYDvwJ2Bd4FLnX3e/NxTiKpzMI4HTsW0WRwjc3Z+8D6zXDtnJAM/voQPTUjUteiJiFmdgowGZgI/Dfx+pCZ7eXu76cp3x54FHga2A/YA5gKrAV+n1J8ECFJqVSeVM9IYBpwOTAdGAfcbWaj3P35vJyciBSd7+4Xbqf9tSz01/nZgUpEROpS7JaQi4Gp7n5j4v0FZvZl4HzgkjTlvw60Bs5w9/XAq2Y2ELjYzCb5F3vZfuzuKzIc9yLgCXe/MvH+SjM7JLH+tFqek4gUKTP4yahwm+yGF8Mtsgv3jx2VSMMVrfuVmTUH9gVmpmyaCRyQYbeRwDOJBKTSI0AvoG9K2VIzW2pmsxIJRmo9qcd9pIrjikgjYQZXjAnz60x6Dm6ZFzsikYYrZh/wrkBTYHnK+uVAjwz79MhQvnIbwFJCS8p4wm2WBcAsMzsoi3rSHtfMzjWzUjMrLS8vT1dERBqQJhbGDTlyV7jiKZjxVuyIRBqmQngQLXWgEkuzrrry/1vv7gvcfYq7l7n7HHefCDwM/LCmx3X3G9y9xN1LunXrVkVoItJQNGsCf/xyGJPl+zNh9gexIxJpeGImISuALWzf+tCd7VspKi3LUJ4q9gF4HhiQRT1V1SEijUzLZnDzcWHAuXNnwBtqCBXJq2hJiLtvBMqAsSmbxgKzM+w2BxhtZi1Tyi8BFlVxuGGE2zTJ9eRyXBFppDq2hNtOCGO5nP4f+PDz2BGJNByxb8dMAiaY2dlmNtDMJhM6mU4BMLOrzGxWUvm/A+uAqWY22MzGAT8FJiWNE3KRmZ1gZgPMbJCZXQWcQBhbpNJk4FAzu8TM9jSzS4BDgOvq+oRFpPj0age3HQ8Vm+D0f8OqitgRiTQMUZMQd59GeCz2MmAecCBwlLsvThTpSRhMrLL8KkKLRS+gFPgLYXyQSUnVNgeuBeYDzyTqPNrdpyfVMxs4FTgjUe504BSNESIimezRFW48Fj74HM59ADZsjh2RSPHTBHY50gR2Io3bvW/CRY/AuD3D8O4azEwaA01gJyJSAE7cE95fFcYQ6d0Bvv+l2BGJFC8lISIiOfreiDAZ4XXPwy7t4at7xY5IpDgpCRERyZEZ/OYw+Gg1/HRW6Lh6wC6xoxIpPrGfjhERKUrNm8KUY6BPRzjvAXj3s9gRiRQfJSEiIjXUoQXcclwYXfXM/8Cn66vfR0S2URIiIlILvTvADcfAsjXw7Rl6dFckF0pCRERqqaQXXDMW5i4JfUQ08oFIdtQxVUQkD47fAxatDI/u9usUnqARkaopCRERyZPvjYCFK+H3c6BfRzh299gRiRQ23Y4REckTM7j6MNivF/xgJpQtrX4fkcZMSYiISB61aBY6qu7YFs65P4yuKiLpKQkREcmzzq3Co7ubtsJZ98HnG2JHJFKYlISIiNSB3TrDlKNDH5GJD8KmLbEjEik8SkJEROrIqF3gykPgmffh50/q0V2RVHo6RkSkDp06GBavgutLoU8HOC/vk6GLFC8lISIidexHB8AHn8NVz8IuHeDoAbEjEikMSkJEROpYE4Nrx8LS1fD9R2DHNmGUVZHGTn1CRETqQctmcOOx0LMtnDMjjK4q0tgpCRERqSedW8HU40MH1QmadVdESYiISH3q1wluOhaWrA5jiKzbFDsikXiUhIiI1LOSXvCnr8DLy+G7D8HmrbEjEolDSYiISARH7gq/GgOzFsLPHtcYItI46ekYEZFIvjEElq2FP82FHm3g4pGxIxKpX0pCREQi+sGXYPkamDwXurcJiYlIY6EkREQkIjO46jD4ZD1c9gS0awHH7xE7KpH6oT4hIiKRNWsC1x8F++8EF88M/UREGgMlISIiBaBls/Do7l5d4fwH4LkPY0ckUveUhIiIFIh2LeDWE8L8Mt+6H+Yvjx2RSN1SEiIiUkA6t4I7T4SOLeH0f8Mb5bEjEqk7SkJERApMj7bw9xPDLZqv3QtvrogdkUjdUBIiIlKA+nSEu8bDDk3gtOmwQImINEBKQkREClTfjjBtfHh65rTp8NYnsSMSya/oSYiZTTSzhWZWYWZlZja6mvJ7m9lTZrbezD4ys5+bmSVtH2dmM82s3MxWm9nzZnZcSh0TzMzTLC3r6jxFRGqiX6fQItK0CZx2jxIRaViiJiFmdgowGfgNMByYDTxkZr0zlG8PPAosB/YDvgf8CLg4qdjBwOPA0Yk6HwTuTZPcrAN6Ji/uXpGfMxMRyZ9dE4mIGZz8L3hFT81IAxG7JeRiYKq73+jub7j7BcBS4PwM5b8OtAbOcPdX3f0e4Grg4srWEHe/0N1/6+5z3f0dd78CKANOSKnL3X1Z8lInZygikge7doJ/nQRtdoBTp8PzH8WOSKT2oiUhZtYc2BeYmbJpJnBAht1GAs+4+/qkdY8AvYC+VRyuHfBZyrpWZrbYzD40sxlmNjzr4EVEIujbMSQiO7aBb94Lj2tkVSlyMVtCugJNCbdWki0HemTYp0eG8pXbtmNm3wF2Bm5PWr0AOAs4HjgNqACeNbMBGeo418xKzay0vFwP7YtIPD3bwd1fhQFd4JwZcP9bsSMSqbnYt2MAPOW9pVlXXfl06zGz8cA1wNfdffH/KnCf4+63uvs8d38GOAV4F7gg7QHdb3D3Encv6datW9VnIyJSx7q0hn+Mg316wAUPwY0vglf1r6ZIgYqZhKwAtrB9C0Z3tm/tqLQsQ3lS90kkILcDp7v7fVUF4u5bgFIgbUuIiEihad8Cbj8Rvrwb/PoZ+H9PwuatsaMSyU20JMTdNxI6jI5N2TSW8JRMOnOA0SmP0o4FlgCLKleY2cnAHcAEd/9XdbEkOrUOIXSKFREpCi2bhdl3v70v3D4fzrkf1m6MHZVI9mLfjpkETDCzs81soJlNJnQynQJgZleZ2ayk8n8nPFo71cwGm9k44KfAJPfQGGlmpwJ3JtY/bWY9EkvnykrM7HIzO9LM+pvZMOBmQhIype5PWUQkf5oY/OxAuPIQeGoxnPQvWLYmdlQi2YmahLj7NOAi4DJgHnAgcFRS/42ewK5J5VcRWj56EW6f/AX4PSGZqXQe0Ay4jtCyUblMTyrTEbgBeIPwNM5OwEHuPje/ZygiUj++MQT+7zhYtBKO+QfM1SO8UgTM1ZspJyUlJV5aWho7DBGRtN76BL49A97/HC4bDROGhkHORGrDzMrcvSTf9ca+HSMiInm0exf4z6lwaF/4xVNw4SOwflPsqETSUxIiItLAtG8BfzsGfjQS7lsAJ0zTnDNSmJSEiIg0QE0MvjsCbj0ByteFfiK3vqzxRKSwKAkREWnADu4Dj3wdRu4MP38SzrwPytfGjkokUBIiItLAdWsDU4+HK8bA7A/gyDvhkXdjRyWiJEREpFEwC0/KzDgtTIB37ozwFM1yjSkiESkJERFpRHbvAvedCj8+AJ5YBIfdDnfMh63qKyIRKAkREWlkdmgK39kv9BXZuztc+gScdDe8kmnWLpE6oiRERKSR6tcJ/j4OrjkcFq6EY++CH8zUsO9Sf5SEiIg0YmZw8iB48owwEd59b8GYW+EPz8E6DXImdUxJiIiI0L4FXHIgzPomHNYPrnseDrwFppRqZl6pO0pCRETkf3p3gL8cBfeeDIO7w1XPwoFT4foXYI2SEckzJSEiIrKdfXrCbSeEZGTv7nD1bBh1C/zuWVi6OnZ00lAoCRERkYwqk5F/nwL77wR/LQstIxc8BC8ujR2dFLtmsQMQEZHCN7wH3HAMvL8qzEEz7bXQiXVwdzhpLzhhD+jYMnaUUmzMNZtRTkpKSry0tDR2GCIiUa3ZCPe8EZKR18qheVM4oj98dS84cJcwFok0HGZW5u4l+a5XLSEiIpKzts3hjKFhea0c7n4d/v0mzHg7tIgc0R+OHgCjlJBIFdQSkiO1hIiIpLdhMzy1GB54Gx5bGFpLOrSAQ/vBIX3DjL66ZVOc1BIiIiIFrUUzOGLXsFRshv++HxKSJxbBvW9CEwt9S8b0hQN2hiE7hts40nipJSRHZlbtL8zdMbOst6WuS36fen3MLGP51P0q16Xbr7pjZqo/UyzJx0l3zsn7JK9LF2+mc0itI139Vanq2OliSXf+qeuTfw/Jx0kXe3XnmhpruvPKFHO6mFLPJdPvqaq/qXTxVvf3nRp/pp9T46/u76yquLO5Ntn8rtIdM9PfeaZrXN3fYzb7VHfOye+rir26/Wu6TzbnmVx28xZn/sfwxMKQkMz/OGxr1QxKesGXdoIRO4VHgVvtkFW1NZJL3DHlM87qPv851KOWEBERKT5Nm4QWkOE94OKR8Ol6eP4jeO7DsFwzJ5Rr1gQGdk2U7QmDu8GuncL+0jApCRERkXrVuRV8ZbewQEhKXlwKLy6Dl5aGp25umx+2tWwWEpNB3cLr7l3Cor4lDYOSEBERiapzKzi8f1gAtmyFdz4NT928+nF4/c8CuOOVbft0bwO7d4b+ncKya6cwK/BO7ULfEykOSkJERKSgNG0Ce3QNy7iBYZ07LFkNCz6Btz6Ftz+Btz6B6W9+cU6b5k1hl/ZhDpzeHaBPB9i5fUhOKuvJ0KVJIlASIiIiBc8MdmoflkP7bVvvDuXr4L3PwrJ41bblhSXbT7o3eAr0agc920KPttAr8dqjLezYJrx2bKlEpb4oCRERkaJlFm7NdG8DX9r5i9vc4bMK+Ohz+HA1HHVRGGJ+yeowCd/r5SGBSdW8aaivW+ttr91aQ9fW0CXx2q01dGkVBm1TwlJzSkJERKRBMgv9TTq3gr13DOt+cfAXy2zcAsvXwvI1216XrYWP10L5Wli0EuZ+FJKZdFo0DfV3SSQlnVqG951aQeeW4bVTy23rO7YM46lIoF+FiIg0WpV9SHZpX3W5jVvCUzwr1m1bPlkflk8TP3+6HhauhM/Ww+qNmetq1SwkJR1bQofK1xbh5w4tQpn73wo/t2+x7bVdi4Y3uJuSEBERkWo0b7qt70g2NmwOrSefrU+8VsDKipCorKyAVRWwckPY/s6nsGpDWLdhS9j/uw+lr7dVs5CQVCYl7ZqHJKVd823v21a+b75tv4Wfhe1tm4fWm0K5haQkREREJM9aNMstaalUsRlaXQSPfiMkJp9XJBKUDaF15fMN25bVG0JC8/6q8PPqjduSmFRjbtv2c7Mm0GaHkJC0aZ70c8q6Njts+7muKAkREREpEC0T38q7d6nZ/hs2hyeCVm8MicmQi8L6PxwJaxKJytqNsGZTKLem8v1GWLZm27Z1m2Dz1vycU1WUhIiIiDQQLZqFpUvrL64ft2du9biHfjDrNoWkpPdF+YsxWfQR+c1sopktNLMKMyszs9HVlN/bzJ4ys/Vm9pGZ/dxSZqYys4MTdVWY2Xtmdl6aesab2etmtiHxemK+z01ERKQYmYVkplOr6jvt1kbUJMTMTgEmA78BhgOzgYfMrHeG8u2BR4HlwH7A94AfARcnlekHPJioazhwFfAnMxufVGYkMA24ExiWeL3bzPbP8ymKiIhIBhZzWmMzex6Y7+7nJK17G/iXu1+Spvz5wNXAju6+PrHuMuB8YGd3dzO7Ghjn7gOS9rsJGOTuIxPvpwGd3X1sUpnHgHJ3P62amKv9hVU31Xm6KcWrmkY95fhZTTVf1dTj2RwzlynWq5ouPfUcMk1Fn805pNaRrv6qVHXsdLGkO//qpotPfl/VtOvpzjU11kxTrqeLOV1MqedS3VTt1f09ZPo9ZDp+dT+nxl/d31lVcWdzbbL5XaU7Zqa/80zXuLq/x2z2qe6ck99XFXt1+9d0n1ymms/ntPS1VUixVCWfcVb3+c+hnjJ3L8lHTMmitYSYWXNgX2BmyqaZwAEZdhsJPFOZgCQ8AvQC+iaVSa3zEaDEzHaopkym44qIiEiexeyY2hVoSri1kmw5cHiGfXoAH6YpX7ltYeL1sTRlmiWOuTRRJt1xe6Q7qJmdC5ybeLsBeDVDfJXlc9qWui75fS7l0+2Xqa4qjtnVzFZUV3822+riHLKtvwoZzy+X2DLFWt0+mbZl8zvKZpuZdQVWpCuXS31VXYds68n2by/TcdJs+8K5VbVPLsfL5n0uf5/Z/D1m2Cfjtauu/lxiyKXeXI5ZlUzXL4Yc/73IRV7PL99x5qG+PfIRR6pCeDomtY3I0qyrrnzq+pqWSXtcd78BuAHAzErrokmqUOj8iltDPr+GfG6g8yt2jeH86qLemB1TVwBb2L71oTvbt1JUWpahPEn7ZCqzGfikmjKZjisiIiJ5Fi0JcfeNQBkwNmXTWMKTLenMAUabWcuU8kuARUllUm/njAVK3X1TUplcjisiIiJ5FnuckEnABDM728wGmtlkQifTKQBmdpWZzUoq/3dgHTDVzAab2Tjgp8Ak39b1dwqws5ldl6jzbGACcG1SPZOBQ83sEjPb08wuAQ4Brssi5htqfrpFQedX3Bry+TXkcwOdX7HT+dVA1Ed0AcxsIvBjoCehw+f33f3pxLapwBh375tUfm/gL8AI4DNC0vHLpCQEMzsY+AMwiNBKcrW7T0k57leBXwP9gXeBS919et2cpYiIiKSKnoSIiIhI4xT7doyIiIg0UkpCUlgdzGVTCBL9X14ws8/NrNzM7jezwdXs09fMPM3y5fqKO1tm9os0cS6rZp+iuHYAZrYow7V4IEP5gr12ZnaQmd2X+J27mU1I2W6J67kkcW2eNLNBWdRb7ZxR9aGq8zOzHczsajObb2ZrzWypmf3dMkxVkbTfmAzXM8dpyWovi+s3NU2cz2VRb8Ffv8T2dNfBzewvVdRZENcvm++B+v78KQlJYnUwl00BGQNcTxgV9lDCI8uPmVnnLPb9MqHPTuXyeB3FWFsL+GKce2cqWGTXDkKMyee2D2Fcm39Ws18hXru2hP5fFwLr02z/MfAD4ALCeX8MPGpm7TJVaFnMGVWPqjq/1oRrd2Xi9XhgF+BhM8tm3KZBfPF6vp2nmHNR3fWDMGBkcpxHVVVhEV0/+OJ59QSOTayv7rMI8a/fGKr/Hqjfz5+7a0kswPPAjSnr3gauylD+fOBzoFXSusuAj0j0tynUhfBB2wIcW0WZvoQvupLY8WZxPr8AXs2hfNFeu0SslwIrgdbFfO2ANcCEpPdGGNX40qR1rYDVwLerqOdq4O2UdTcBcwrp/DKU2StxrfauosyYRJmusa9ZdecHTAVm5FhPMV+/G4EF1ZQp1Ov3he+BGJ8/tYQkWN3NZVOo2hFawj7Loux0M/vYzJ618FRRoeqfaEJdaGZ3mVn/KsoW7bUzMwO+Bdzh7uuqKV4s165SP8JAgv/7HCau0dNUPbdTNnNGFarKidKz+SyWJm7hzDKzQ+oyqFo6MPF395aZ3Whm3aspX5TXz8zaAqcSEpFsFNr1S/0eqPfPn5KQbaqayybtnDJknoOmclshmwzMIwzclska4IfAyYTm1FnANDP7Rt2Hl7PnCePBfAU4h/D7n21mXTKUL+ZrN5bwj8VNVZQppmuXrPJ3n8vnsHK/dPtUzhlVkBL/+fk9cL+7p86LlWwpofVuPDCOcOtxlpkdVPdR5uxh4HTgMEKz/gjgcTNrUcU+RXn9gK8BLYBbqylXqNcv9Xug3j9/hTB3TKGpi7lsCoqZTQIOBA509y2Zyrn7CsI/kJVKLUyQ9mPgjrqNMjfu/lDy+0RHuPeAMwiD4qXdLeV9wV+7hHOAF9x9XqYCxXTtMsj1c5hpn3TrC0KiD8gdQEfguKrKuvsCwhdXpTlm1peQaD5dRyHWiLvflfT2FTMrAxYDRwNVjcVUVNcv4Rzg3+5eXlWhQrx+1XwP1NvnTy0h29TVXDYFxcz+AJwGHOru79WgiueBAfmNKv/cfQ3wGpljLbprB5Bo1j6e7Jt/kxXDtat8oinXuZ2ymTOqYCQSkH8AQ4DD3L0mMRbD9cTdlxBmP68q1qK6fgBmNgwooWafRYh4/ar4Hqj3z5+SkASvu7lsCoaFYfG/RvjDe7OG1QwjNC0WtMQ12ZPMsRbVtUsyAdgA3FVNuXSK4dotJPyD9r/PYeIajabquZ2ymTOqICTukU8jJCCHuHuVj5JXoRiuJ4kWuJ2oOtaiuX5JziX8W/FYDfePcv2q+R6o/89f7N65hbQApwAbgbOBgYT7ZWuAPontVwGzksp3SFywu4DBhHt9nwM/iH0uac7tL4nYDiVkrJVL26Qyqed3RuKPdSCwB6HpcCNhaP3o55RyftcCBxP6SuwPzEicb9Ffu6SYDXiLlCe4iu3aEXrkD0ss64CfJ37undj+k8S1GJe4NncRksN2SXXcBtyW9L4fsJYw/9PAxGd4IzC+kM6PcAv834SnsPZJ+Sy2quL8LgJOIPzPeVDiejswrsDOr23isziS0MF7DOEL6sOGcP2SyrQGVpH0FElKHQV5/cjue6BeP3/1enGLYQEmErLbDYSWkYOStk0FFqWU35twT6+CkNVeTgE+4pn4g0+3/CLT+RG+yF5P/HF9DpQC34h9LhnOr/KDspHwD/w9wF4N4dolxXtI4pqNSLOtaK4d2x5XTF2mJrYb4ZHrpYlr8xQwOKWOJ4EnU9YdDLyY+OwuBM4rtPNj26PT6ZYJmc6P0JfnHcK4FZ8CzwBHFeD5tSI8FfFx4rO4OLF+l4Zw/ZLKnEm41dArQx0Fef2q+Nv7RVKZev38ae4YERERiUJ9QkRERCQKJSEiIiIShZIQERERiUJJiIiIiEShJERERESiUBIiIiIiUSgJERERkSiUhIiIiEgUSkJEREQkCiUhIlKwzOwkM9tgZn2S1k02s3fNbMeYsYlI7WnYdhEpWGZmwAvAS+5+jpn9kDAPxyh3fztudCJSW81iByAikom7u5n9DHjAzN4FLiVMQa4ERKQBUEuIiBQ8M5sNjACOdfeHYscjIvmhPiEiUtDM7FBgKGGK8eWRwxGRPFJLiIgULDMbCjwFXAwcDbR19yPjRiUi+aIkREQKUuKJmNnA39z9l2Y2GJhP6BPyZNTgRCQvlISISMExs87As8DT7v7/27VDIoChGAqCUVIfVVRtVZnSL6AzR3bhQ4EH8hz7OzPX7t7ZccBvRAgAkPCYCgAkRAgAkBAhAEBChAAACRECACRECACQECEAQEKEAAAJEQIAJD7fKN40Mho04QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot fitted pdf\n", "yMin = 0.\n", "yMax = f(0., thetaHat, xiHat)*1.2\n", "fig = plt.figure(figsize=(8,6))\n", "xCurve = np.linspace(xMin, xMax, 100)\n", "yCurve = f(xCurve, thetaHat, xiHat)\n", "plt.plot(xCurve, yCurve, color='dodgerblue')\n", "\n", "# Plot data as tick marks\n", "tick_height = 0.05*(yMax - yMin)\n", "xvals = [xData, xData]\n", "yvals = [np.zeros_like(xData), tick_height * np.ones_like(xData)]\n", "plt.plot(xvals, yvals, color='black', linewidth=1)\n", "plt.xlabel(r'$x$')\n", "plt.ylabel(r'$f(x; \\theta, \\xi)$')\n", "plt.xlim(xMin, xMax)\n", "plt.ylim(yMin, yMax)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAESCAYAAADuVeJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3yN5//H8dcngxAxE8QOIfaMvdUopdqiVAcdVimlNVt0oVpq1Ld2S9Votaq1qvZesfeovXcSSci6fn+c8EtVIiE5d3LyeT4e5xHn3Nd935+rJe9z3eO6xRiDUkop9ThOVheglFIqddDAUEoplSAaGEoppRJEA0MppVSCaGAopZRKEA0MpZRSCaKBoZRSKkHsFhgi8omImIdelx+zThMR2SIiwSJyXUT+EJFi9qpZKaXU/7P3COMo4B3rVSauhiLiA/wBbAAqAA2BDMDS5C9TKaXUw1zsvL9IY0y8o4pYKgGuwEBjTBSAiIwAVouIpzHmenwre3p6mkKFCj1VsalFdFQ0kRFRREVGERkRRXRkFFFR0URHRhMVFUV0VPSD99FRUURFRhMdFc3T3OPv7OyEs4szzq5OOLu44BLz09nVGRcXZ5xdnW3LXZxxdtEjn0qlFjt37rxujPF61DJ7B0ZhEbkAhAPbgEHGmJNxtA0AIoB3RGQakBHoAOx4XFgAFCpUiICAgCQqO2UIvnWHk/vOcGr/WU7tO8OpA2c5feAcYXfuPrJ9OjdXMmV1x93LHfcsGfHI5k6mbO54ZMuER7ZMZMpm+9zN3Q039/Qxr///c3RUNGHBdwkJCrX9DAwl+OYdgm4EE3QjmMDrQdy+GsitK4HcvhpI0JU7PGqqGWcXZ7J4epA1Zxay5sxMttxZ8crnSc78OfDK74lX/hzkLOCJR7ZMyf2fUCn1GCJyJs5l9ppLSkSaAh7AESAn8DFQHChljLkRxzq1gfmAJ7bDZ7uBpsaYq3G07wx0BihQoEClM2fi7HeKF3YnjIObj3F850mO7z7J8Z0nuXzq/7vtkT0TPmUK4FO6ALl9cpItV1ay5cpCttxZyeLpQaas7qRzS2fXmqMiowi8HvQgQG5fDSLwWhC3rgYSeC2I29dsn9+4eIvrF24SHRX9r/Wz585KoZg+FSpdAJ8yBShYMh9uGdPbtR9KpWUistMY4//IZVZNPigimYCTwJfGmG8esTw3sB5YCMzFFjafxSxuYIyJfnid2Pz9/U1qGmEYYzi57wwBy/cS8PceDmw4TGREFAB5iuTCt2JhilYsTJFyBfEpW5Ac3tkQEYurfnJRUVHcunybq+ducO3cDa6cvsrpQ+c4feAcZw6eI/xuBAAigneRXBQuWxDf8j74VvTBt4IPObyzWdwDpRxTfIFh70NSDxhj7ojIQaBoHE26AyHGmH73PxCR14BzQA1gY/JXmbzCQu6yY9luti7eyc6/93Lz8m0AfMoU4MWezajYqBx+lYs45KEaZ2dnPPPmwDNvDqj272VRUVFcOnmVU/vPcvrAWU4dOMvJvWfYuGDbgzbZcmXBt4IPRcr74Fe5CKVq+JEtV1Y790KptMWywBARN2yHpNbE0SQjEPXQZ/ffp9qzqCGBIWxdvIsNC7YS8Nce7oWFkzmHBxUblcW/cTkqNS6HZ57sVpdpKWdnZ/IV9SZfUW9qv1T1wechQaGc3HuGE7tPcWLPKU7sOsWulfuJivz/kVipmsUpWd2PUjX9KFgyH05OqfavilIpjj3PYYwCFgFnsZ3DGAzUAcoYY87EXAFVxRjzTEz7BsBK4FNgDrZDUsOBkkAJY0xIfPtLSYekoqKi2LFsD0umrmDn8r1EhEeSI082ar1YldqtqlG6dnGcnZ2tLjNVCr8bzvFdpzi0+SgHNx/h4Kaj3L4WBNjO85StW5JydUtRrl4pCpXOrwGi1GOkiHMYIjIPW0B4AteArcBgY8yhmOUzgHrGmEKx1mkH9AOKAWEx6/S/v058UkJg3LoayF/TV7NkygqunLlGdu9sNHilFrVbVaV41aL6yysZGGO4+M9lDmw8wr71h9i39iCXT18DIHMOD8rWLUmFBmWo8Exp8hXLk6rPAymVHFJEYNiblYFx4cQlfvnqD1b8uI6I8EjKNyjN892aUP15f1xcLTsKmGZdPn2VfesOsWftAfasPsC1c7aL8jzzZqd8g9JUbGg7HKjnQJTSwLCbUwfOMnfEAtb9vBlnVxeadKzHCz2bUbBEPrvWoeJ2fwSye9UBdq/ez941Bwi8HgxA0UqFqdykPJWbVqBE1aI4u+hhQpX2aGAks0snrzDzk59ZPXsjbu7padG1MS/1bq6XfqYC0dHR/LPnNNuX7WbHX7s5vOUY0dEGj+yZqNKsAtVbVMa/STncM2e0ulSl7EIDI5kE3Qjmx09+YcmUFTi7OPNiz2a06fs8mbN7JOt+VfIJvnWHXSv2sXXJTrYv3U3QjWBcXJ0pV780NV+oQo2WlfWLgHJoGhhJLDo6mr++X8P0gbO5czuEpm8/w2tDWqf5y2EdTVRkFIe2HGProgA2/7mD88cuISKUqFaUWi9WpVarqnj75LK6TKWSlAZGEjp//BJfd5zAoS3HKFO7BO9NeBufMgWTfD8qZTHGcPbweTYu2M6mhds4vusUAMX8i1C3TXVqt66m4aEcggZGEjDGsHzGWv7Xczqu6V3p9k1HGr5eRy/LTKMunbrChl+3sv7XLRzd8Q8AfpWLUKd1deq0qU7uQjktrlCpJ6OB8ZTu3A5hbNfJrPtlC+XqlaL/j+/hlS9HkmxbpX6XTl1h/XxbeBwLsIVH8Sq+NHy9LvVfqanntFSqooHxFA5sPMyI18Zz4+ItOn7WljZ9n9e7slWc7ofH6rkbOLn3DK7pXKj2vD+N36iHf5Nyeh+OSvE0MJ7Q6jkb+Krj/8hV0JOBs3tRvEpc8yQq9V8n9pxixcx1rJq9nsDrwWTNmYVn2teiUYd6FClXyOrylHokDYwn8NuYxUz6YCZl65bk09/7kSmrexJWp9KSiPAIdizbw4pZ69i6KIDIiCh8K/jw7FsNaNC+lkPORqxSLw2MRDDGMG3AbH75+g9qt6rKgFk97f4gIuW4gm4Es3ruRv76fjX/7DlNOjdXareuxnOdGlG6VnG9iEJZTgMjgaKjoxnTeTJ/fb+aFl0b0/3bt/R8hUo2x3edZNm0Vayas4HQoDDyF8/Lc50a0qhDXT1RriyjgZEAxhjGdZvKkikrePWjVnT4rK1+21N2ERZyl3W/bGHp1BUc3nqcdG6u1G9XixbvNsHPv4jV5ak0RgPjMYwxTOw9g9/HL6Vd/xd4a3h7DQtliX/2nmbRxL9ZNXs9d0PuUcy/CC26NaFe2xr6bHNlFxoY8TDG8P2gOcwbuZCXej1H1286aFgoy4UEhrBi1noWTVzO2cMX8MjmTuOO9Xn+3SbkKZLb6vKUA9PAiMfyGWsY9dZ3NO/SiJ7fddKwUCmKMYZ96w7x58TlbFywDRNtqNy0PC27N8W/STl9CJdKchoYcThz+Dw9Kg+geFVfvvx7sJ7gVina9Qs3WDJlJUumrODWlUDy+Obm+W5NaNyxnl6aq5KMBsYj3Au7x3vVBnHr8m0m7v5aZ5pVqUZEeAQbftvGn9/9xcFNR3FzT0/zLo1p1ae5/j1WT00D4xG+7TGNP79bzvClg6j8bAU7VqZU0jmx+xS/frOINXM34uziTOOO9WnbryXehXXmXPVk4guMNHkA9MDGw/z53XJe7NlMw0Klar4VfBgwqyc/HB1P4471+XvGGjr69eTL18dz+uA5q8tTDibNjTAiwiPoVrEfYXfuMu3AN2TIlMGC6pRKHtcv3uS3bxazePLf3A25R80XKvPKwJfwq+xrdWkqldARRiy/fPUnZw6dp+d3nTQslMPxzJOdLqPeYPbpibw+pA371h2iR9WB9G/yOXvXHsRRvyAq+0hTI4wzh87RrWI/arxQmY/n9bGoMqXsJzQ4jMWT/ubXbxZx60ogJWv48drg1vg3LqeXkKtHShEjDBH5RETMQ6/Lj1lHROR9ETkiIvdE5JKIfPkk+4+KimLU2xPJ4JGB7uPffrJOKJXKZPTIwMt9WzLr5P94b8I7XDt3nUFNh9GzxkdsX7ZbRxwqUex9SOoo4B3rVeYx7UcD7wL9gRJAM2D9k+x4wZglHNl2nB7fvk22nFmeZBNKpVrpM6Tn+XebMPP4t7w/qTM3L93io+eG07P6ILYt3aXBoRLEboekROQToLUxpnQC2/sBB4CyxpjDid1f7ENSV89eo6NfL6o0Lc/Q3/rqUFyleRHhEayYuY65IxZw+fQ1ivkX4bXBranWvJL++0jjUsQhqRiFReSCiJwSkXkiUjieti2Bk8CzInJSRE6LyEwRyRnXCiLSWUQCRCTg2rVrDz7f/GcAEfcieGfk6/qPQSnANZ0rzTo15Iej4+kztStBN4IZ0nIk7/r3Z+Pv24iOjra6RJUC2TMwtgEdgaZAJyA3sFlEcsTRvjBQEGgXs97rQHFgkYg8sm5jzBRjjL8xxt/Ly+vB57tW7sO7cC7yFfVOoq4o5RhcXF1o+vYz/HBkHB9+/y6hwWF82moU3Sr2Y8NvW/VQlfoXuwWGMWaZMeYXY8w+Y8xKoHnM/jvEU1t64HVjzHpjzAZsoVEFqJzQ/UZFRrF3zUEqPvO40yVKpV0uri406Vif7w+Npf+P7xFxL4LP2ozm/dqDObztuNXlqRTCsvswjDF3gINA0TiaXAIijTHHYn12HIgECiR0P/s3HCY0OIxKjcs9ca1KpRXOLs40fK0OUw98Q+8pXbn0z2V6Vh/EiNfGcfXstcdvQDk0ywJDRNywHWK6FEeTTYCLiMR+5FhhwAU4k9D9rP15M27u6ancVKcAUSqhnJ2dafbOM8w49i3tB73ExgXbeLN4L6YPmkNIUKjV5SmL2PM+jFEiUldEfESkKvAr4A7MjFk+QkRWxVplJbAL+F5EKohIBeB7bOdCEvQovciISDb8tpUaLSvr08qUegIZPTLw5hev8MORcdR6qSrzvvydDr49WDhhGZERkVaXp+zMniOMfMBcbPdiLADuAdWMMfdHC97Ag9GEMSYa23mOq9juvVgOnAdaxix7rN2rDxB0I5h6bWsmWSeUSotyFvBi4E+9mLD9SwqVLsD/en7PO6X7sGHBNj0xnoY49NQg3Rr1Zv7oRfwROJP0GXSEoVRSMMawbckupvafxdnDFyhTuwTdxnSkaMX4rpJXqUVKug/Drg5tPYZvhUIaFkolIRGhWvNKTNk7ml4TO3PuyAW6Vx7A6Le/4+blW1aXp5KR4waGgWM7/qFEtWJWV6KUQ3J2caZ5l0bMODaeVr2bs/Kn9bzp14t5IxcSfjfc6vJUMnDYwAgLucvd0HuUrlnc6lKUcmjuWdzpMuoNph4YQ9l6JZk+cDZvl3yfNfM26fkNB+Ow5zAK5vYxpW5XY/7V6bhnzmh1OUqlGbtW7mPyhz9yct8ZSlQrSpdRHShVw8/qslQCpclzGHduhVCxUVkNC6XsrGLDsny3cyQfTH+XK2eu836tj/mi3Tdcv3jT6tLUU3LYwIgIj6TmC1WsLkOpNMnZ2Zln36zPjGPjeWPoy2z5M4B3SvVm6dSVepgqFXPYwADwzJvd6hKUStMyuLvx+tA2TN47Gt8KPozpMpm+z3zKhRNxTfCgUjKHDoyQQJ3CQKmUIF9Rb75aOYTek7twfNdJOpf9gDnDFxARHmF1aSoRHDow7tzWwFAqpXBycqJZp4ZMPzSWqs9V5IeP59KtYj/2b0j089GURRw6MPRacKVSHs882Rky/0O+WDSAuyH36FN3CKPe+o6gG8FWl6Yew6EDwz2LXiGlVEpV9blKTD3wDW37tWTlT+t5u+T7rP1Z791IyRw6MDyyZbK6BKVUPDK4u/HOl68xcedIchb0YtgrYxn64ldcv3DD6tLUIzh0YLhn1RGGUqmBT5mCjN88jM5fv8GuFft4u1RvFk9eoc8WT2EcOjBc07taXYJSKoGcXZxp80ELJu8dRbFKhRnXbQq96wzh1P4EPy9NJTOHDgwRsboEpVQi5fX15quVQ+n7Q3fOH71It0r9mTbgJ+6G3rO6tDTPwQPD6gqUUk9CRGjcoR4/HBlHo9fr8PNXf9CpTB8ObT1mdWlpmmMHhpNDd08ph5c5hwcfTH+XUWs+wUQbetcezNwRv+u5DYs49m9UvTxPKYdQrm4pJu3+mjqtq/H9R3Po3/hznczQAg4dGOeP6Xw1SjmKTFndGTTnfT6Y1o0jW4/TuewHrPtls9VlpSkOGxhOTsLBzUetLkMplYREhGffasDEXV+Rt6g3X7Qbw/BXxxJ0U+8StweHDQw39/Qc33XS6jKUUskgX7E8jN3wOR0/b8f6+VvpXPYDdizfY3VZDs9hA8M1vStXz163ugylVDJxdnHm1Y9a8e3W4bbDVU2HMa7bFMLuhFldmsOyW2CIyCciYh56XU7gukVFJFhE7iR0fy7pXLh56Rbh93T6ZKUcWdGKhfkuYCSt+7RgyZSVdCnflwMbdQbc5GDvEcZRwDvWq8zjVhCRdMA8YH1idpTOzXaX9+kDZxNdpFIqdUnnlo4uo96wXX5rDH3qDuWnz3/ViQyTmL0DI9IYcznW61oC1hkJ7APmJ2ZHbu5uABzZdiLxVSqlUqWydUoyec8oGrxai5lDf2b3qv1Wl+RQ7B0YhUXkgoicEpF5IlI4vsYi8hzQHOiZ2B25pnMhq1dmPfGtVBqT0SMD/We+x9erhlKxYdnHtt+/4TBbF++0Q2Wpnz0DYxvQEWgKdAJyA5tFJMejGouINzAVeN0Yk6Br5kSks4gEiEjAtWvXyFXIS6dJVioNEhHK1y/92Hbh9yK4dPIK33SayC9f/2GHylI3uwWGMWaZMeYXY8w+Y8xKbCMHJ6BDHKv8BEw0xmxNxD6mGGP8jTH+Xl5e5MiTnesX9G5QpdSjpUvvSoP2tchfPC+hwXp11eNYdlmtMeYOcBAoGkeTBsBQEYkUkUhgOuAe875zQvZRuGxBzhw8z9VzenmtUuq/oqKimNBjOtlyZaHjZ+0A9ER5PCwLDBFxA4oDcc3fUQYoH+s1BAiL+XOCToA37lAPYwwrflz39AUrpRzO/FGLOHP4PH1/6A7YAiT2YxE0PP7NnvdhjBKRuiLiIyJVgV8Bd2BmzPIRIrLqfntjzIHYL+ACEB3z/lZC9uldOBfl6pVi1ewNydAjpVRqdO7oBQA2/r6N1XM20HtKV9JnSE9UVBTOzs4WV5eyudhxX/mAuYAncA3YClQzxtx/nJY3UCSpd1q9hT+TPpjJ1XPXyZnfM6k3r5RKRe6G3mNct6lkyubOtXM36D7uLQoUz0t0dDTOzs5ER0fj5OTEjUu3OHPoPFsXBZDFMzOvftzK6tJTBHue9G5njMljjElnjMlrjGlljDkUa3lHY0yheNafYYzJlNj9VnjGdm+gXo+tlHLLmJ7PFw3gXug9Lp+6Srl6pR4si4qKehAWMz6ey4Zft5C3qDeHth7lqzcnWFh1yuGwc0ndV6h0fjLn8ODAxiNWl6KUSgEyuLsxYtnHNHqjLn/PXEvwrTs4OTk9GGFMHzSbHHmy0/qDFrTs/iz9ZvTALaObns8gDQSGk5MTflV8ObLtuNWlKKVSkK6jO+BXuQgrZq7jyhnbpBOLJ63ALUN6GnesR15fbwB+HrmQq2ev/etkeFplz3MYlilRtSgBf+3h0qkrePvksrocpVQKUbBkftzc3R7MPXf9wg38qviS3TsbAJsWbufYzpP0ntIFgKVTVxJ4PRgXV2fafPi8ZXVbxeFHGABN3qyPm3t6vu0+TYeVSql/yVXQi2y5shJ+L4JTB85SzL8Ibhltz9NZOGEZLXs0xRj45es/+P3bpRQokZc1P29izvAFVpdud2kiMHLm9+TNL15hx197WDNvk9XlKKVSoHTpXfH2ycW33aexes4Gvmj7DfXb1qRoRR8Clu/h3JELfDS3NzVfqMI7X77GuaMXiIqKsrpsu0oTgQHwfPcmFC5XkLkjFugoQyn1SO+OfRP/JuUJvhXCG5+0pVmnhlw8cZkj24/T8I26FCqVn5CgUHav2k/hsoXS3H0baeIcBoCzszPNuzRm/LtTOb7rJMUqJfktH0opB9B+0Ev/ej/rs/lUb+FPubq2S3CPbj9BSGAo1VtUIizkLntWH+DojhMULJmf+u1qWlGy3aSZEQZAvbY1cE3nwpq5elhKKfV4d26HkDVnFl7u2xKAvesOsuOvPXj75KRkdT+GtRvD1kUBuGfOyM9fLeTXbxZZXHHySjMjDACPbJnImDkDd0PuWl2KUioVyODhRmREJJ+9PBqf0gU4e/g8hUoVoM2Hz/PVmxMwxtB7SlcAvPLn4Myh8xZXnLzSVGAAhAbfJaNHBqvLUEqlAs7Oznz+xwBmf/Eb4XfDeWXgSxQuW5A/v1vOsR3/MO3AmAdtb1y8RVRkFMYYh71nI00FhjEGj+yZ2LBgGy++/xyeebJbXZJSKhV49eNWD4LgxqVb7PhrN/1/fO/B8lP7z/DLqD/5ZEFfhw0LSGPnMESEob99yK0rtxnQ+HMCrwdZXZJSKpW4HwTuWTLiniUjeYva7gQPCQplSMuRtO3XkhJVizr0VZhpKjAASlYrxud/DuDSySt83OJLoqOjrS5JKZWKREVGcfXcdaYPnM1fP6xh2Ctjqdbcn5d6PQegIwxHU75+ad6f3IUj247rszKUUoninjkjI/8eTPjdCAKvBdHo9bp0H/8W4PgPXBJH7aC/v78JCAiIc3l0dDTvVRvErSu3mXF0POnc0tmxOqWUo3GUk90istMY4/+oZWlyhAG2WWzfGvYK187dYNPCHVaXo5RK5RwhLB4nzQYG2A5NZfTIwP71hx7fWCmlEsER55lK04Hh7OJMqZp+7Fq13yH/5yqlrPHbmMW8X2sw549fsrqUJJWmAwOg0Rv1uHD8EvO//tPqUpRSDsIzXw7OH71Itwp9WTx5hcOcDE/zgVGvbQ3qvlydGUN+5viuk1aXo5RyAHXbVGfq/tGUrOnHuG5T+LjFCIJuBFtd1lNL84EhIvT8rhNZc2ZmvD5gSSmVRDzz5mDEso/oPu4tdq/cz5guk60u6aml+cAAyJzdg1cGvsSRbcc5uPmo1eUopRyEk5MTL7zXlNeHvszGBdvYunin1SU9FQ2MGI071sMjeyZmf/Gr3v2tlEpSrT9oTqFS+fm2x7RUPSWR3QJDRD4REfPQ63I87euJyB8icklEQkVkn4i8lVz1ZXB347WPWxOwfC8TekzXQ1NKqSTjms6VPtO6cftqIAOafEHwrTtWl/RE7D3COAp4x3qViadtDWA/0BooDUwEpohI++Qq7sVezXi5b0sWTfqbyR/+qKGhlEoyJaoW5ZMFfTlz8ByDmg0nJCjU6pISzd6BEWmMuRzrdS2uhsaY4caYj40xm4wxJ40xE4EFQKvkKk5EeOfLV2nZ/Vl+G7OYFT+uS65dKaXSoMrPVuDjn/twfOdJBjT+nFtXbltdUqLYOzAKi8gFETklIvNEpHAi188M3IproYh0FpEAEQm4di3OLIqXiPDuuDcpWb0YU/vNIuhm6r8UTimVctRoWZkhv37Aqf1n6Vl9EGcOp56n9NkzMLYBHYGmQCcgN7BZRHIkZGURaQ48A0yJq40xZooxxt8Y4+/l5fXEhTo5OdHzu04E3bzDjI/nPfF2lFLqUWo8X5nRaz/lXlg4vWp8xJ41B6wuKUHsFhjGmGXGmF+MMfuMMSuB5jH77/C4dUWkJjAH6GmM2Z7MpQJQpFwhGrSvxdqfN+m5DKVUkvOr7Mv4LcPxzJudwc9/yYk9p6wu6bEsu6zWGHMHOAgUja+diNQClgFDYs5j2E2pGsUJvhXC5dNX7blbpVQakbtQTkauGIJHtkwMaTmSm5fjPOKeIlgWGCLiBhQH4pydS0TqYAuLT40xY+1V231+lYsAMG/E70RGRNp790qpNCCHdzY++6M/wTfuMPTFrwm7E2Z1SXGKNzBEJEhEPGP+HBzz/pGvx+1IREaJSF0R8RGRqsCvgDswM2b5CBFZFat9PWxhMQmYLSK5Y15PfnIikXwr+NDmgxYsnbaKgc9+kapvuFFKpVy+FXwY8FNPjgX8w0fNR6TY0Ij3iXsi0gGYZ4y5F/PnOBljZsa7I5F5QB3AE7gGbAUGG2MOxSyfAdQzxhSK9f5R+zxzv018HvfEvcRYMWsdYzpPJod3Vr5YMoiCJfIlyXaVUiq2tT9vYsSr4yhVqzjDFg8kQ6YMdq8hvifuJfgRrSLykjFmwSM+F6CfMWbk05WZtJIyMACO7jjB4Oe/JH3G9EzYNoIsnpmTbNtKKXXf/dCo3LQCn/85wO5P8kuqR7TOFpFpIpIx1obzAWuA3k9ZY4rnV9mXTxf258bFW3z+8jd6TkMplSzqta1Jl9Ed2LZkF6vnbLS6nH9JTGBUBaoBe0TEX0TaYpu6IwwolxzFpTQlqhalz9Su7F17kJ8+/9XqcpRSDqplj2cpUa0o373/A7evBVpdzgMJDgxjzD7AH9gIbAFmAUONMU2NMVeSqb4Up+FrdajduhoLv11GaHDKPDGllErdnJ2d6TO1G6FBoQxvP46I8AirSwISf1ltOaAucAIIB6qIiEeSV5XCte3bkpDAUJZMXmF1KUopB1WoVH76TO3G7lX7Gdd1aoq4gTjBgSEig4H1wB/YgqMS4AfsF5HayVNeyuRX2ZdKjcvxw8dz2bF8j9XlKKUcVKM36vLa4NYsn7GGGYPnWR4aiRlhdANaGGP6GGPCjTFHgerAPGBlslSXgg2a04sCJfMx9IWv2LVyn9XlKKUc1BufvMyzb9ZnzvAFjHhtHHdD71lWS2ICo6wx5u/YHxhjIo0xA4BGSVtWypc5uwdfrRhCfr88DH7+S/ZvOGx1SUopByQi9JnWjbeGtWftvM30qTuEa+dvWFJLYk56X49n2fqkKSd1yZzDg5ErBpOrUE4GP/8lJ/edsbokpZQDEhFeGfginy7sx4Vjl3iv2kBL7gbXZ3o/paxeWfjyr4/IkMmNgU2HcelUmrlgTCllZ9Vb+DP0tw+5cfEWmxbusPv+NTCSQLRnCf0AACAASURBVM4CXoz462PCw8L5tNUowu+GW12SUspBlW9QmlwFvVg12/4HdjQwkkihUvnp/+N7/LPnNJP6xDutllJKPTEnJycavlaHXSv2cWCjfc+damAkoWrNK/Hyh8+zaNLfLJ2a5i4cU0rZSZu+z5PbJyfDXx1H8K07dtuvBkYSe3PYK/g3KceYLpOZP3qR1eUopRyQe+aMDJrzPjcv3Wb0OxOJjo62y341MJKYi6sLny7sT5021ZnS90e+/2iO5TfbKKUcj19lXzqNfI1Nv2/nfz2/t8vvGZdk30MalC69K4Pm9MI9c0bmjvgdF1cX3vjkZavLUko5mJfef44bF28yf/Qi3LNk5K1h7ZN1fxoYycTZ2ZneU7oQHRXNrM/mk9nTgxd6NLW6LKWUAxEROn31OiGBocwd8Tu5C+WkWaeGybY/PSSVjESE3lO6UP15f77r9QMrZq2zuiSllIMREXpO7ET5+qWYNuCnZH2UtAZGMnN2ceajue9Trl5Jvu74P5ZM0RlulVJJy9nZme7j3yYkKIwZg+cl2340MOwgfYb0fLF4IJWblmds1yn8Nmax1SUppRxMoVL5adn9WZZMWcnRgH+SZR8aGHaSPkN6PlnQl9qtqjLpg5ksmvT341dSSqlE6PDpy2TLnZVxXScTFRmV5NvXwLAj13SufDS3N1Wfq8iEHtPYunin1SUppRyIexZ33h3TkeO7TvHHhL+SfPsaGHZ2/5xGkQo+DGs3hr3rDlpdklLKgdRpU50qzSowfdBsTu1P2hm07RYYIvKJiJiHXpcfs04ZEVknImEickFEhoiI2Kvm5JIhUwaGLR6IV/4c9G/0uZ4IV0olGRHhw+nvkimrO5+3HUNYyN0k27a9RxhHAe9YrzJxNRSRzMAK4ApQGegJ9AX6JH+ZyS9brqyM3zKcig3LMLbrFMZ3n0ZkRKTVZSmlHEC2XFnpP6sn549eZGrfWUm2XXsHRqQx5nKs17V42r4KZAQ6GGMOGGN+A0YCfRxhlAHYvgEsGkCbD1qwaOJyhrzwFffCrHv8olLKcVR8pgwtujVm6bRVXD0X5/PvEsXegVE45tDSKRGZJyKF42lbHdhgjIn9WKnlQB6gUHIWaU/Ozs50/voNek/uQsBfexjUbDghQaFWl6WUcgAv920JwK9JNBGqPQNjG9ARaAp0AnIDm0UkRxztc2M7HBXblVjL/kNEOotIgIgEXLsW3+Al5WnWqSEDZ/fi4Kaj9G/0GbevBVpdklIqlctV0IsG7WuxdOrKJHkOuN0CwxizzBjzizFmnzFmJdA8Zv8d4lvtofcSx+f39zHFGONvjPH38vJ6+qLtrH67mgz97UNO7T/Le1UHcurAWatLUkqlcm988jLGGKb2f/pzGZZdVmuMuQMcBIrG0eQy/x1J5Iz56bAPzq7ewp9v1n1G+L1IetX4iG1L9F4NpdSTy10oJy/3bcmauZvYt/7QU23LssAQETegOHApjiZbgNox7e5rBFwETidvddbyq+zL/7aPIF8xbwY/P5K/vl9tdUlKqVSsbf8X8MqX46mfz2PP+zBGiUhdEfERkarAr4A7MDNm+QgRWRVrlTlAKDBDREqLyEvAAOAbkwaeSOSZNwej131GxUZlGf3ORL1XQyn1xNwypuflfi05uOko+9Y9+SjDniOMfMBcbPdiLADuAdWMMfdvRfQGitxvbIwJxDaiyAMEAP8DRgPf2LFmS2Vwd+Ozhf2o0qwCY7tOYcHYJfr0PqXUE2n6dgOy5crCrM/mP/HvEXue9G5njMljjElnjMlrjGlljDkUa3lHY0yhh9bZb4ypY4xxM8Z4G2M+TQuji9jSuaVj6G99qflCZSb2mcHnbb+x60PflVKOIX2G9LT/qBV71x5k/a9bn2gbOpdUKpAuvStDfv2Qt0e8yuaFO+haoS8HNh62uiylVCrTomtjfCv4MLH3D4QGhz1+hYdoYKQSTk5OtOv/AmM3fo6LqzMf1BvK/FF/6iEqpVSCObs40/O7Tty4eIvfxy1N9PoaGKlM8SpF+W7nV9R8qSpT+s3i6zf/R/jdcKvLUkqlEiWqFqVM7RKs/3VLotfVwEiF3DNn5ON5vXlj6Mus+HEdHzb4hBuXblldllIqlaj1YlVO7jvDhRNx3dXwaBoYqZSTkxOvD23DkPkfcGrfWbqW/5Ady/dYXZZSKhWo1aoqIsLfM9Ymaj0NjFSudqtqTNg+gqy5sjCo6TCmDfhJp0lXSsUrZ35ParT0Z/HkFYmaIVsDwwEULJmfCdtG8Fynhvz81R982OATbl3VyQuVUnF7sddzBN0IZtVPGxK8jgaGg0ifIT3vT+7CoDnvc2LXKXpWG8iZQ+esLksplUKVrVOSgiXzsXL2+gSvo4HhYOq3q8notZ8SfjeCnjU+0vMaSqlHEhHqtK7OgQ1HEnzRjAaGA/Kr7Mu3W4eTq6AXg5oOY2yXyYQEhlhdllIqhanTpjrGGNbO25Sg9hoYDipnAS/GbR5G6z4tWDZ9FW+X6s3mP3ZYXZZSKgUpVCo/pWsVZ8G4JQm6WEYDw4FlcHejy6g3GL9lOFk8MzP0xa8Y9soYnYtKKfXAy31bcvXsddbPf/yNfBoYaYBfZV/+t+NLOn7ejg2/baNL+Q/Zv0HnolJKQdXnKpKnSC5Wz9342LYaGGmEi6sLr37UinGbvsA1nQsf1h/KjCHz9J4NpdI4JycnfCsW5tzRi49va4d6VAriV9mXibu+5pnX6jD7i994r9ogTuw5ZXVZSikL5S+Wh8unrhIRHhFvOw2MNCijRwb6zejB4F/6cP3CTXpUGcgPH88l/F78f1mUUo6pcLmCREdFcyzgZLztNDDSsDqtqzP94BgatK/FnOEL6Faxr442lEqDytcvjYiwa8W+eNtpYKRxmXN40G9GD4YtGURIYCg9qw3ij//9pc/ZUCoNyZzDg6KVCrN7zf5422lgKACqNK3ApN1fU6FhGSa8N51PW4/Sy2+VSkOKVizM6f1n422jgaEeyOqVhc//HEDnr99g66KdvFW8F0umrCAqKsrq0pRSyaxA8bwE34p/RggNDPUvTk5OtPmgBd9uHU4+vzyM7TqFdyv1Z9eq+IeqSqnUrWCpfI9to4GhHqloxcJ8s+4zPv65D6FBofRv9BlDXhjJ5dNXrS5NKZUMStUsjmt613jbaGCoOIkIddtUZ/qhsbw9vD27V+2nU+k+zB/1J1GRephKKUfiljE95eqVjLeNZYEhIoNExIjIhMe0ayIiW0QkWESui8gfIlLMXnUqSOeWjnYDXmT6wTGUf6Y0U/rNonuVARzdccLq0pRSSahQqQLxLrckMESkGtAJiPeiXxHxAf4ANgAVgIZABmBpcteo/itnAS8+W9ifIb9+yO2rgbxXbRDfdJrErSu3rS5NKZUEwu7cjXe53QNDRLIAs4G3gcc9taMS4AoMNMacMMbsAUYARUTEM3krVY8iItR+qSrTD47hpV7N+HvmWjoW68kvX/+hd4orlcqFBofGu9yKEcYU4FdjzOoEtA0AIoB3RMRZRDyADsAOY8z1hxuLSGcRCRCRgGvXriVt1epf3LO40/WbjkzdP5oydUowtf9PdCrdm21Ld1ldmlLqCbmkc4l3uV0DQ0Q6Ab7A4IS0N8acBhoBnwL3gECgDNA8jvZTjDH+xhh/Ly+vJKlZxS+/X16+WDSQ4cs+wiWdCx83H8Gw9mP1MJVSqZB75ozxLrdbYIiIHzAceNUYE57AdXID04EfgcpAPSAY+EVE9AqvFKRyk/JM2v01HT5ty6YF23i75Pv89cManWJEqVTELWP6eJfb85dudcATOCAikSISCdQF3o15/6hKuwMhxph+xpjdxpj1wGsx69WwW+UqQVzTufLa4NZM2jOKgqXyM/rt7+hZ4yN2LN+jwaFUKhAaHBbvcnsGxkJsh5PKx3oFAPNi/vyoUUdG4OEL/u+/1xFGClWgeF5Gr/2UPlO7cvPSLQY1HUavmhocSqV0j5s/zm6/dI0xt40xB2K/gBDgZsx7IyIjRGRVrNWWABVFZKiIFBWRisAPwDlgp71qV4nn5ORE07efYcax8bw/qTM3Lv5/cOxcsdfq8pRSj3AvNP6zBSntW7o3UOT+m5grqdoDLYHdwHJsV009a4yJf5YslSK4pnPluc6NmHFsPL0nd+HmpdsMaPIF/Zt8zvFd8T+sRSllX+5Z4j/pLY56iMDf398EBARYXYZ6SPi9CBZP/JvZw34j6EYwDdrXouPn7fD2yWV1aUqleRN7z+DdsW/uNMb4P2p5ShthKAeXLr0rL73/HD+e+JZXBr7Ipt+381bxXox/dypXz/3n1hqllB3dCdTpzVUK5J7FnbeGtWfGsfE0ebMBy6avooNvDw0OpSwSfi+CzQt3xNtGA0NZyjNvDt6f1JkZx779T3Bcv3DD6vKUSjMClu/hzm0dYahUIFdBr/8ER8diPZk+aA4hjxkmK6We3vYlu8iYOUO8bTQwVIpyPzi+PzKOmi9WYd6Xv/OG73ssGLtEJzdUKhnt23CYMrVLxNtGA0OlSN4+uRj4Uy++CxhJkfKFmNhnBm8V78VvYxbriEOpJHbraiDnjlygdC0NDJWKFa1YmJF/D2b4so/wzJedSR/M5JX8XZnw3nTOH7todXlKOYSNC7YBUKlR2Xjb6X0YKlU5tvMfFn67jLXzNhERHknlphV4qddzVGpUFhGxujylUqUeVQcQfjeCyXtG4eTkpPdhKMdQrFIR+s3owewzE3lj6Muc2HWSgc9+QdcKfVkxax2REZFWl6hUqnL64DmO7viHZ99s8NgvXRoYKlXKlisrrw9tw0+nJ/Lh9+8SHRXNVx0m8EaRHswfvYiQoPifHKaUstn8h+3ei/qv1HxsWw0MlaqlS+9Kk471mbJvNF8sHkge39xM6fsj7Qt0ZdIHM7l08orVJSqVogUs30PRij5ky5X1sW01MJRDEBGqNqvIqNWfMGH7l1RpVpGF3y6jQ9H3GNzyS3au2KtTqyv1kFtXbnNoyzEqNS6foPZ60ls5rOsXbrB48gqWTFnJ7auB5C+el5bdn6Vxh7pkyBT/DUpKObro6GgGNRvOvnWHmLjrKwqWyAeAiOhJb5X2eObNQcfP2jH7zET6zexBRg83Jrw3nVfyd2Vq/5+4dl6nHlFp19zhv7Pz7710H/fmg7B4HB1hqDTl0Jaj/DZ2CRt/24o4OVGvXQ1a92mBb3kfq0tTym72rjtIv2c+pV67mgyY1fNfV0fFN8LQwFBp0qVTV/h93FL++n41YXfuUq5eKZ55tTY1X6hC5hweVpenVLIJvnWHLuU/JJ1bOibuHPmfw7MaGErF4c7tEJZOXcniySu4dPIKzi7OVHimNPXa1qROm+pkcHezukSlkowxhuHtx7Lht22M2/QFfpV9/9NGA0OpxzDGcHzXSdbP38L6X7dy6eQVMmbOQINXatGsU0OKVixsdYlKPbWl01YxpvMk3vziFdoPeumRbTQwlEoEYwwHNx1h6bRVrPtlM+F3IyhcriDPvtmABu1rkcUzs9UlKpVoe9cdZEDjzylXvzTDlgzE2dn5ke00MJR6QsG37rB6zkaWz1jD8Z0ncXF1pmrzSjzzah2qNqtAOrd0Vpeo1GNdOnmFHlUHksXTg/FbhpMpq3ucbTUwlEoCp/afYfkPa1gzbxM3L9/GPUtG6rSuToP2tShTp0Sc39iUstKlU1cY0OQLgm8E8+22EeT19Y63vQaGUkkoKjKK3asPsPKndWz6fTt3Q+6RI0826rWtSYP2tShasbDOnKtShBO7TzGo2TAiwyP5YvFASlb3e+w6KTIwRGQQMAz4nzGmRzztBOgFdAV8gJvATGPMgPi2r4Gh7OFu6D22Lgpg9dyN7Fi2m8iIKPIV86bBK7Wp/0pN8hXLY3WJKo3atXIfn7YahXvWjIxY9hEFS+ZP0HopLjBEpBowFwgCNjwmML4BmgN9gf1AFsDbGLM0vn1oYCh7C7oZzMYF21kzdwN71x7CGINPmQLUfKEKtV6qSuGyBXXkoZKdMYY/v1vOxN4zKFAiL8OXDsIzb44Er5+iAkNEsgC7gE7AEOBAXIEhIn7AAaCsMeZwYvajgaGsdP3CDdbP38rGhds4sOEIxhi8C+eiRsvK1HqxCiWqF9NzHirJ3Qu7x7huU1nx4zqqPleRAbN6xnuC+1FSWmD8DJw2xvQXkbXEHxj9gLeBScB72Oa+Wgf0NcZcjW8/Ghgqpbh1NZAtfwaw8fdt7Fm1n4jwSLLmzEL1Fv7UaFmZCs+UJn2G9FaXqVK5S6eu8GmrUZzce4bXh7bh1Y9b4eSU+OkCU0xgiEgnbOciqhtjwhMQGJOAjsBebIekDDAqZnF1Y0z0Q+07A50BChQoUOnMmTPJ0Q2lnlhIUCg7lu1m08LtbF+6m9DgMNJnSEfFRmWp3sKfqs9VJHvubFaXqVKRiPAIfh+3lFmfzcfF1YUBs96j6nOVnnh7KSIwYg4vbQRqG2OOxHy2lvgDYwq2Q1d+xphjMZ8VA44C1Ywx2+Lan44wVEoXfi+CfesOsXVRAFsWBXD17HUA/CoXoVpzf6o1r0SR8oX0vIeK0+7V+5nw3nTOHr5AtRaV6D7uLXIXyvlU20wpgdER+AGIivWxM7ZRQzTgboy599A6nwKDjDGusT4TIBxob4yZH9f+NDBUamKM4eS+M2xdvJNtS3ZyZNsJjDFkz52Vio3KUrGh7ZXDW0cfCq6dv8GUfrNYO28TuX1y8u7YN6ne4pG/4xMtpQRGVuDhSdd/AI4Dw4GD5qFiRKQxsBzwNcb8E/NZEeAEUNUYsz2u/WlgqNTs1pXbbF+2m50r9rJrxT4CrwcDUKh0fsrXL025eqUoW7ckmbPrzLppyeXTV/l55EKW/7AGRGjbryXtBryQpOfAUkRgPHLnDx2SEpERQBVjzDMx752AHcAd4P2Y1cYC6YEaD5/DiE0DQzmK6OhoTu49w84V+9i1ci8HNx3lXlg4IoJP2QKUq1uK8g1KU7ZOyURfEaNSh/PHLjL3y99ZOWs9zs5ONO5Yn3YDXnjqw0+PkpoCYwZQzxhTKFYbb2A88CwQBqwA+hhjrsS3bQ0M5agiwiM4uuMf9q45yN51Bzm46QjhdyNwchJ8KxamfD1bgJSsXgz3LBogqdmJPaf45es/WPfzZlzSufBc50a0+fB5vPIl/L6KxEqxgZGcNDBUWhF+L4Ij246zZ/UB9qw5wOGtx4iMiEJEKFQ6P6Vq+FGqZnFK1fAjt09OPYmewoUGh7F23iaWTlvJ0R3/4OaenufffZbWfZqTLVfWZN+/BoZSaUhYyF0ObznGwU1HObjlKIe3HCM0OAyAbLmy4FfFl+KVi+JXxRe/ykXwyJbJ4opVRHgEO//ex9pfNrF54Q7C7tylUOn8NHunIc+8Vtuu56o0MJRKw6Kiojhz8DwHNx3h8PbjHNl2gnNHLjxYnq+YN8X8i+Bb3gffij74VvDRELGD0OAw9q07xKbft7Fp4XaCb4Xgkc2dWi9V49m3G1CialFLRoMaGEqpfwkJDOFowEmObj/Bke3HOb7rJNfO3XiwPHchL4pU8MG3vA8FS+WnYMl85PXNjbOLTmfypKIiozi640TMxQv7OLz1OFGRUWT0yECNFypT7+UaVGxUFtd0ro/fWDLSwFBKPVbg9SBO7D7F8V2n+GeP7eeF45ceLHdxdSafXx4KlsxHwZK2EClQIh95fHOTLr21v+RSopDAEA5vO8GhzUc5tPUYh7ceIzQoDBHBt6IPlRqWpWKjspSqWTxF/ffTwFBKPZGwkLucPXyBs4fOc+bQOc4cPs+Zg+e5fOoq9393ODkJuQrlJJ9fHvIV9SZPkdzk8c1NniK5yFXIy/JvzPYQdCOYM4fOc+bQeU7sOsnBLUc5c/A8xhicnIRCpQtQoloxKjQoTYVnypA5R8q9fya+wHCxdzFKqdQjg7sbfv5F8PMv8q/P74be49yRC7bX0YucP3aRc0cvsn/9Ie6G/P+EDSJCdu+seOX3xCt/DnLmy4FXfk888+UgR55seObNTo482VJFqESER3Dl9DXOH7vE+WMXOX/0IueOXeTs4Qvcvhr4oJ17loyUrF6MOq2rU6qGH35VfHHPnNHCypOOjjCUUknGGMPtq4FcOHGZS/9c4eI/l7l27gbXzl+3/Tx3g7uh9/6zXlavzGTPkw2PbJlwz5IR96wZcc+ckYweGXBzdyNDJjfc3NPbfsb6s+1zN1zTu+Ca3tX2Sufy4GSxMYaI8EjCw8K5FxZO+N1w7oWGExoUSkhQGKGBtp8hgaGEBoUSGhRmWxYcRmhQGEE3ggm+EUzQjTsPrjS7L3MOD/IV86ZAiXwxh+lsL898OZ5oltiUQkcYSim7EBGy5cpKtlxZKV2z+H+WG2MIvnmH6xducuPiTa5fsL1uXLjJjcu3CLkdyqWTVwgJDCUkMJSw4DCioxP/pdY1nQtOzk6E340gMV+KM3pkIIOHGxkzZ8Q9cwayemWmQPG8eGTPROYcHuQq6EW+Yt7kK5YnRR9WSi4OO8IQkWuAPeY39wSu22E/KVFa7juk7f5r3x1XQWOM16MWOGxg2IuIBMQ1fHN0abnvkLb7r31Pm31PvQfalFJK2ZUGhlJKqQTRwHh6U6wuwEJpue+QtvuvfU+D9ByGUkqpBNERhlJKqQTRwFBKKZUgGhhKKaUSRAMjFhF5V0ROichdEdkpIrXjaesmIjNEZJ+IRMQ8bvZR7dKJyGcx270nImdFpGeydeIpJFP/24vIHhEJFZHLIvKTiOROtk48oUT2vZ6I/CEil2L6tU9E3npEu7ox27orIidFpGvy9uLJJHXfReQlEflbRK6JSLCIbBOR55O/J08mOf7fx2pfS0QiReRA8lRvXxoYMUSkLTAOGA5UADYDy0SkQByrOAN3gQnAkng2PRfb88g7A35AG2BfEpWdZJKj/yJSE5gFzARKAS8AJYHZSVr8U3qCvtcA9gOtgdLARGCKiLSPtU0fYGnMtioAI4BvRaRVcvXjSSRH34G6wGrguZhtLgV+j+8XsVWSqf/3t50N+BFYlQylW8MYoy/blWLbgKkPfXYcGJGAdScAax/xeWMgEPC0un8W9f9D4MxDn70J3LG6v0nV91jtfwF+i/V+JHD8oTbTgC1W9ze5+x5Hm+3AaKv7a8/+AwuAocAnwAGr+5oULx1hYDtsBFQC/n5o0d/YvlE8qReAHUAfETkvIsdFZLyIpKjnXyZj/zcB3iLSQmw8gXbYvnGmCEnY98zArVjvqz9im8sBfxFJEXN5J2PfH8UjAW3sKjn7LyLvArmBL56mxpRGA8PGE9shlisPfX4F2//0J1UYqAWUA1oBPbAdnprxFNtMDsnSf2PMFuAVbIegwoFrgAAdnnSbyeCp+y4izYFn+PcNXbnj2KZLzD5TguTq+8NtugP5sB2eTEmSpf8iUgbbyOJVY0xU0pSaMmhg/NvDdzHKIz5LDKeY9dsbY7YZY5ZjC41WIpLrKbabXJK0/yJSEhgPfI7tm9yz2P4hTn7SbSajJ+p7zHmaOUBPY8z2BGzzUZ9bLTn6fr9NK+BrbL887TF79JNIsv6LSHpgHvChMeZUUhdqNX0ehs11IIr/fqvIyX+/fSTGJeCCMSYw1meHY34WeMptJ6Xk6v9AYLsx5uuY9/tEJATYICIfGWPOPcW2k8oT911EamE7vDbEGDPxocWX49hmJHDjiatNWsnV9/ttWmEbVbxhjPnz6ctNcsnRf29sF3b8ICI/xHzmZFtFIoFmxpiHD4GlGjrCAIwx4cBOoNFDixphu2riSW0C8jx0zqJYzM8U820rGfufEds/yNjuvxdSgCftu4jUAZYBnxpjxj6iyRag4SO2GWCMiXjyipNOMvYdEXkZ+AnoaIz5NWkqTlrJ1P8LQBmgfKzXJOBEzJ+f5t+T9aw+655SXkBbbMfZ3wFKYLvU7g62h4mA7bLIVQ+tUxLbX4J5QEDMn8vHWp4JOAfMx3ZZaU3gADDf6v7aqf8dgQigG7bzOTWxXQSw0+r+Pk3fgXpACLZDLbljvbxitfGJaTM2ZpvvxOyjldX9tUPf28X8f+/1UJvsVvfXHv1/xD4+wUGukrK8gJT0At4FTgP3sH3zqBNr2Qzg9EPtT2M71vmv10Nt/LBddRGK7dvH/wAPq/tqx/6/BxyM6f8lbMd881nd16fpe8z7//T7Ef996gK7YrZ5CuhqdT/t0XdgbRxt1lrdV3v9v39o+5/gIIGhs9UqpZRKED2HoZRSKkE0MJRSSiWIBoZSSqkE0cBQSimVIBoYSimlEkQDQymlVIJoYCj1lGIeqmNiZuNVymFpYCiVSCKyVkQmpJbtKpVUNDCUUkoliAaGUokgIjOwTfnRPeYwlAEKxSwuF/P86lARCRCRig+tW0NE1sUsvyAiE0Ukc1zbFZFCIuIsItNjnjkdFvMQrn4iov92ld3pXzqlEqcXtplof8A2lbU3tgkmwTZR3QCgIrYpzGeLiMCDh+r8DfyJ7YFaL2GbrPH7x2zXCdscZC9jmxzvI2AQtkfdKmVXOpeUUokkImuxTSbXI+Z9PWAN8KyxPSTr/sN1NgL5jTHnReRHIMIY83as7ZQHdgO5jDFXH95uPPv/EvA3xjw8fbpSyUofoKRU0tkX688XY37mBM5je+Kgr4i0jdXm/jNBigBX49qoiHTFNv12QSAD4EoKep6KSjs0MJRKOrEfjHR/6O4U6+c0YMwj1rsQ1wZjAmYs8CG2h+8EAd2BF5+2WKUSSwNDqcQLB5wTuc4uoJQx5kQit1sL2GaMeXC5rYgUSeS+lUoSetJbqcQ7DVSJuYrJk4T9OxoZs84kEakgIr4i0lxEJse13ZgroY4BFUWkqYgUFZHB2K6mUsruNDCUSrxR2EYDh4BrQIHHrfB/7d2xDYQwEETRqYwOKYAW6IKY6K4LqjCByW+QzLuPnwAAAE9JREFULnwvt6WNviyt5DHGN8mSuYJ7JPlkblVdP+7dkuyZPxWez/n1L1PAS7akAKh4YQBQEQwAKoIBQEUwAKgIBgAVwQCgIhgAVAQDgMoNzTwoHSrGFUcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Make a contour plot of lnL = lnLmax - 1/2\n", "if nfreepar == 2:\n", " plt.figure()\n", " m.draw_mncontour('theta','xi', nsigma=1, numpoints=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 }