{ "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", "# truncated in [xMin,xMax].\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 scipy.stats import chi2\n", "from iminuit import Minuit\n", "import matplotlib.pyplot as plt\n", "from matplotlib import container\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, par):\n", " theta = par[0]\n", " mu = par[1]\n", " sigma = par[2]\n", " xi = par[3]\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(par):\n", " pdf = f(xData, par)\n", " return -np.sum(np.log(pdf))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Initialize Minuit and set up fit:\n", "# initial parameter values are guesses,\n", "# error values set initial step size in search algorithm,\n", "# limit_param to set limits on parameters (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", "par = np.array([theta, mu, sigma, xi]) # initial values (here equal true values)\n", "parname = ['theta', 'mu', 'sigma', 'xi']\n", "parstep = np.array([0.1, 1., 1., 1.]) # initial setp sizes\n", "parfix = [False, True, True, False] # change these to fix/free parameters\n", "parlim = [(0.,1), (None, None), (0., None), (0., None)]\n", "m = Minuit.from_array_func(negLogL, par, parstep, name=parname,\n", " limit=parlim, fix=parfix, errordef=0.5, pedantic=False)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# do the fit, get errors, extract results\n", "m.migrad() # minimize -logL\n", "MLE = m.np_values() # max-likelihood estimates\n", "sigmaMLE = m.np_errors() # standard deviations\n", "cov = m.np_matrix(skip_fixed=True) # covariance matrix\n", "rho = m.np_matrix(skip_fixed=True, correlation=True) # correlation coeffs.\n", "npar = len(m.np_values())\n", "nfreepar = len(cov[0])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "par index, name, estimate, standard deviation:\n", " 0 theta = 0.204551 +/- 0.052736\n", " 3 xi = 5.107878 +/- 0.644563\n", "\n", "free par indices, covariance, correlation coeff.:\n", "0 0 0.002797 1.000000\n", "0 1 -0.018131 -0.531774\n", "1 0 -0.018131 -0.531774\n", "1 1 0.415591 1.000000\n" ] } ], "source": [ "npar = len(m.np_values())\n", "print(r'par index, name, estimate, standard deviation:')\n", "for i in range(npar):\n", " if not(m.is_fixed(i)):\n", " print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n", " \"{:.6f}\".format(MLE[i]), \" +/- \", \"{:.6f}\".format(sigmaMLE[i]))\n", "\n", "print()\n", "print(r'free par indices, covariance, correlation coeff.:')\n", "for i in range(nfreepar):\n", " for j in range(nfreepar):\n", " print(i, j, \"{:.6f}\".format(cov[i,j]), \"{:.6f}\".format(rho[i,j]))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAF7CAYAAABhB6n0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5gV5dnH8e8Ny8LS21KVjkgRQZYiSDBRiIotWMCSaAwSEEXFJGJ5hSQETCwJEYiCRuxij0EREAVBmruCgCiiFBWpCkgvu/f7x+yS7Y1zdrb8Ptc113JmnnnmnmPh3qeauyMiIiISKeXCDkBERERKFyUXIiIiElFKLkRERCSilFyIiIhIRCm5EBERkYhSciEiIiIRFRN2AMVF3bp1vVmzZmGHISIiUmSSkpJ2unt8pOtVcpGqWbNmJCYmhh2GiIhIkTGzTdGoV90iIiIiElFKLkRERCSilFyIiIhIRCm5EBERkYhSciEiIiIRpeRCREREIkrJhYiIiESUkgsRERGJKCUXIiIiElFKLkRERCSilFyIiIhIRCm5EBERkYhSciEiIiIRpeRCREREIkrJhYiIiESUkgsRkSjZsWMHY8aMYceOHWGHIlKklFyIiETJsGHDSExMZPjw4WGHIlKklFyIiETB888/T8WKFZkxYwYVKlTgpZdeCjskkSJj7h52DMVCQkKCJyYmhh2GiIhIkTGzJHdPiHS9arkQERGRiFJyISIiIhEVanJhZjeZ2QYzO2RmSWbWO5eylcxsmpmtNLOjZjYvh3KxZvan1HoPm9nXZjYiai8hIpKDiRMn0qpVK+Li4jjvvPM0a0TKjNCSCzMbCEwAxgGdgUXATDNrksMt5YFDwETgrVyqfgE4DxgCtAGuAFZGKGwRkXy55557ePDBB5kyZQpLly5l/fr1/OEPf4joMyZPnkzz5s2pVKkSXbp0YcGCBXneM378eLp27Ur16tWJj4/noosuYvXq1YWue9y4cZgZN998c5ZrW7Zs4brrriM+Pp5KlSrRrl075s+fX/AXLURcBblnzJgxmFmGo0GDBhnK5Od7a9asWZZ6zIz+/fuf0DuXSO4eygEsBaZmOrcOGJ+PeycC87I53w/YA9QtaDxdunRxEZFI+Oijj9zMfNGiRcfPTZgwwePj4yP2jBdffNFjYmJ8ypQpvmbNGr/55pu9SpUqvmnTplzv69evn//73//2VatW+cqVK/3SSy/1+vXr+/fff1/guhcvXuzNmjXzjh07+vDhwzNc27Vrlzdv3tx/+ctf+tKlS339+vX+7rvv+po1a3KM7brrrvPRo0dH9J3zc8/o0aO9TZs2vmXLluPH9u3bC/y9bd++PUMdH3/8sZuZT5s2Lcf4wgYkejT+jo9GpXk+FGKBY8AVmc5PAubn4/6ckovJwLsErSHfpiYr/wSq5lWnkgsRiZSBAwd6nz59Mpx7/vnn3cwi9oxu3br54MGDM5xr1aqVjxo1qkD17N2718uVK+dvvvlmgerevXu3t2jRwufOnet9+vTJklzcdddd3rNnzwLFkldyUZh3zs89o0eP9vbt2xco1uy+t8zGjh3rNWrU8P379xeo7qIUreQirG6RugTdHNsynd8GNMhaPN9aAGcBpwOXATcTdJFMy66wmQ0xs0QzS1RfqIhEwtGjR/nvf//LgAEDMpw/ePAgNWrUyFJ+3LhxVK1aNdcjczP+kSNHSEpKol+/fhnO9+vXj0WLFhUo3r1795KSkkKtWrUKVPeQIUO4/PLL+dnPfpZtvW+88Qbdu3dn4MCB1KtXj06dOjFx4sS0XwQLrDDvXJB71q9fT+PGjWnevDmDBg1i/fr1ucaT+XvLzN154oknuPbaa6lcuXJer1fqxIT8/Mz/llk25wqiXOr9V7v7HgAzuxmYZWb13T1DMuPuU4ApAE3aJrg7mJ3A00WkzFuxYgUHDhzgzjvv5O677z5+/ujRo3Tu3DlL+aFDh3LllVfmWmfjxo0zfN65cyfJycnUr18/w/n69evz7rvvFijeW2+9lU6dOnHmmWfmu+6pU6fy5Zdf8swzz+RY7/r165k8eTK33347o0aNYsWKFdxyyy0Ax8dnjBs3jnHjxh2/5/Dhw5gZDz744PFzM2fOpHfv3oV65/ze0717d6ZNm8app57K9u3bGTt2LD179uTTTz+lTp06+freMpszZw4bNmxg8ODBOX1FpVpYycVOIJmsrRT1yNqaURBbgM1piUWqz1J/Nsmt7m374ZNt0OlE2k1EpMxbu3YtsbGxrFy5Ekv328pVV11Fr169spSvXbs2tWvXLtSzLNNvQ+6e5VxuRo4cycKFC1m4cCHly5fPV91r167l7rvvZsGCBcTGxuZYd0pKCgkJCYwfPx6Azp07s27dOiZNmnQ8ucicWN155500btyYESP+N8Evc2JVmHfO657zzz8/w/UePXrQokULnnrqKUaOHJmlvty+tzRTp06la9eudOrUKdfYSqtQukXc/QiQBPTNdKkvwayRwvoQaGRmVdOdOyX156bcbjTghawDpkVECmTPnj3UrVuX1q1b06pVK1q1akXNmjVZsWIFl19+eZbyhekWqVu3LuXLl2fr1q0Zzm/fvj3Lb+k5uf3223nhhRd47733aNGiRb7rXrx4MTt37qRDhw7ExMQQExPD/PnzmTx5MjExMRw+fBiAhg0b0q5duwx1tG3blq+//vr459q1ax//jlq1akW1atWynIuLiyv0Oxf2e6patSrt27dn3bp1+f7eMtf/n//8hxtvvDHHZ5R2Ya5z8TBwvZkNNrO2ZjYBaAQ8CmBm481sbvobzKydmXUiGLNR1cw6pX5O8zzwPfCkmbU3s14E011fcfftuQVTsxK8+QXsOxK5FxSRsqdu3brH++PTjB8/njPPPDPbJvShQ4eyYsWKXI+EhIyrM8fGxtKlSxfmzJmT4fycOXPo2bNnnjHeeuutPP/887z33nuceuqpBar70ksvZdWqVVniGzRoECtWrDjemtGrVy/Wrl2boY4vvviCpk2b5hlfdgrzzoX9ng4dOsTnn39Ow4YNM5zP7XtLb9q0aVSsWJFBgwbl9VqlVzRGieb3AG4CNgKHCVoyfpLu2jRgY6byGwnGVGQ4MpVpA8wGDgCbCWagVMsrlrYdu3iTf7i/sKqgY21FRP5n+/btHhcX53/60598w4YN/sADD3h8fLx/8cUXEX3Oiy++6BUqVPCpU6f6mjVrfMSIEV6lShXfuHHj8TKPPPKIt2nTJsN9N910k1erVs3nzp2bYdrk3r17C1R3etnNFlm2bJnHxMT42LFjfd26df7SSy959erVfeLEicfL7N27N0MM2R2HDx8+oXfOzz133HGHz5s3z9evX+9Llizx/v37e7Vq1TKUyc/35u6ekpLirVu3zjJDpbiiNE1FLY5Hly5d/Jyn3S95sQD/VEREsvHyyy97kyZNPC4uzvv165fr2g4nYtKkSd60aVOPjY31M844w+fPn5/h+ujRo9N+ATsuu1/QgCxTQPOqO73skgt39xkzZnjHjh29YsWK3rp1a58wYYKnpKRkiS+34/333z/hd87rnoEDB3rDhg29QoUK3qhRIx8wYIB/+umnhfre3nvvPQd86dKlOX5fxUm0kgvtipoqISHBh05J5M8LYNY1cGrdsCMSERGJLu2KWgQGtIXY8vDip2FHIiIiUnIpuUindhz8vCW89hkcOhZ2NCIiIiWTkotMBrWHPYdh1ldhRyIiIlIyKbnIpOfJcHJ1rXkhIiJSWEouMilnMLA9LP4WNu4OOxoREZGSR8lFNq5oFyQZ0zWwU0REpMCUXGSjQVU4tzm89CkcSQ47GhERkZJFyUUOrjkNdh7UwE4REZGCUnKRg580DQZ2Prsy7EhERERKFiUXOShnQevFks2w7oewoxERESk5lFzk4sp2UKEcPLcq7EhERERKDiUXuahTGc5vBa+ugYNHw45GRESkZFBykYdfdoQfj8B/vwg7EhERkZJByUUeujaCU+rAs+oaERERyRclF3mw1IGdn2yDVdvCjkZERKT4U3KRDwNOhbgYtV6IiIjkh5KLfKheES5pA/9ZG+yYKiKSHzt27GDMmDHs2LEj7FBEipSSi3y6tiMcPAavrAk7EhEpKYYNG0ZiYiLDhw8POxSRIqXkIp9OqwddGsLTn0CKhx2NiBR3zz//PBUrVmTGjBlUqFCBl156KeyQRIqMuetvSoCEhARPTEzMtcyba+GWd2DaJfDTZkUSloiISNSYWZK7J0S6XrVcFMB5rSC+Mjz1SdiRiIiIFF9KLgogtjxcexrM2wgbd4cdjYiISPGk5KKArj4NYsoFYy9ERHJyww03YGYZjp49e4YdlkiRUHJRQPWqBPuNvLwG9h8JOxoRKa6uv/56mjZtypVXXsnLL7/M8uXLefPNNyNW/5gxY7IkLw0aNMj1ng8++ICLL76Yxo0bY2ZMmzYt23KTJ0+mefPmVKpUiS5durBgwYIM15s1a5bl2WZG//79j5fZu3cvt912G02bNiUuLo6ePXvy0UcfndA75xVXTrZs2cJ1111HfHw8lSpVol27dsyfPz/bsuPGjcPMuPnmmwt8vSDPKe2UXBTC9Z2C/UZe/zzsSESkODp69CjXXnstI0aMYPr06Vx++eV06tSJunXrRvQ5bdq0YcuWLcePVatyX+lv3759dOjQgQkTJhAXF5dtmenTp3Prrbdy9913s3z5cnr27Mn555/P119/fbzMRx99lOG5H3/8MWbGlVdeebzM4MGDmTVrFk899RSrVq2iX79+nHvuuWzevDnH+K6//nrGjBlT6Liys3v3bnr16oW789Zbb/HZZ5/xyCOPUK9evSxllyxZwtSpU+nYsWO2deV2vSDPKRPcXYc7Xbp08fxKSXHv/7z7uc8EfxYRSS8pKcnLlSvnBw8ejNozRo8e7e3bty/0/VWqVPEnn3wyy/lu3br54MGDM5xr1aqVjxo1Kse6xo4d6zVq1PD9+/e7u/uBAwe8fPny/sYbb2Qod8YZZ/g999yTYz3XXXedjx49OttrhYnL3f2uu+7ynj175lrG3X337t3eokULnzt3rvfp08eHDx9eoOv5fU5xAyR6FP5ODbXlwsxuMrMNZnbIzJLMrHcuZSuZ2TQzW2lmR81sXh51n2Vmx8xsdeTjhutOhy++h8XfRrp2ESnpatasSUpKCn/5y1/YtGkTKSkp2ZYbN24cVatWzfXIrel//fr1NG7cmObNmzNo0CDWr19/QnEfOXKEpKQk+vXrl+F8v379WLRoUbb3uDtPPPEE1157LZUrVwbg2LFjJCcnU6lSpQxl4+LiWLhwYZHEleaNN96ge/fuDBw4kHr16tGpUycmTpyIZ1qGYciQIVx++eX87Gc/y7aevK7n9zllRWjJhZkNBCYA44DOwCJgppk1yeGW8sAhYCLwVh511wKeBuZGLOBMLjoFasfBkyui9QQRKalatGjBY489xoMPPkizZs2IiYlh+fLlWcoNHTqUFStW5HokJGS/BEH37t2ZNm0aM2fOZOrUqWzdupWePXvy/fffFzrunTt3kpycTP369TOcr1+/Plu3bs32njlz5rBhwwYGDx58/Fy1atU488wzGTt2LJs3byY5OZlnn32WxYsXs2XLluPlMidXzz33XJZzCxYsKFRcadavX8/kyZNp0aIFs2bN4tZbb2XUqFFMmjTpeJmpU6fy5Zdf8uc//znbOvK6nt/nlCUxIT57JDDN3aemfr7FzM4DhgF3ZS7s7vuBoQBm1hGomUvdTwBPAQZcHsmg01SKgas7wKSPgmmpzXKLRkTKlIceeoj777+f22+/nbPPPpv69evToUOHLOVq165N7dq1C/WM888/P8PnHj160KJFC5566ilGjhxZqDrTmFmGz+6e5VyaqVOn0rVrVzp16pTh/DPPPMMNN9zASSedRPny5TnjjDO46qqr+Pjjj4+XGTp0aIZxGnfeeSeNGzdmxIgRx881btyYXbt2FTiuNCkpKSQkJDB+/HgAOnfuzLp165g0aRI333wza9eu5e6772bBggXExsZmuT+v6/l9TlkTSsuFmcUCXYDZmS7NBk5orpaZ3QQ0AMaeSD358avTg2mpar0QkTRLlixh1KhRzJs3j3HjxtGvXz9OP/10ypcvn6XsiXaLpFe1alXat2/PunXrCh173bp1KV++fJbWgO3bt2dpNUg7/5///Icbb7wxy7WWLVsyf/589u3bxzfffMOyZcs4evQozZs3P16mdu3atGrV6vhRrVq1LOfi4uIKHFd6DRs2pF27dhnOtW3b9vhA0MWLF7Nz5046dOhATEwMMTExzJ8/n8mTJxMTE8O8efNyvX748OF8PaesCavloi5BN8e2TOe3AecWtlIzOw0YDfRw9+S8MlozGwIMAWjSJKfemJzVrxJ0j7y0BkaeCTUqFiJoESlV3nnnHdq1a0f79u3zLJv5N/fsNG7cOF/PPXToEJ9//jk//elP81U+O7GxsXTp0oU5c+ZwxRVXHD8/Z84cLrvssizlp02bRsWKFRk0aFCOdVapUoUqVaqwa9cuZs2axd/+9reox5Ver169WLt2bYZzX3zxBU2bNgXg0ksvzdL19Otf/5rWrVtz991307hxY3r16pXj9bTWjLyeU+ZEY5RoXgfQCHCgd6bzo4HP83H/RGBepnMVgU+BX6Y7NwZYnZ+YCjJbJL1V29yb/MP90cRC3S4ipczjjz/u5cqV8/vvv99Xr17tP/zwQ1Sec8cdd/i8efN8/fr1vmTJEu/fv79Xq1bNN27c6O7ujzzyiLdp0ybDPXv37vXly5f78uXLPS4uzv/4xz/68uXLfdOmTcfLvPjii16hQgWfOnWqr1mzxkeMGOFVqlQ5Xm+alJQUb926dZYZHGneeecdf/vtt339+vU+e/ZsP/30071bt25+5MiRDPFs2bIl1+Pw4cP5jiu7d162bJnHxMT42LFjfd26df7SSy959erVfeLEiTl+t9nNBsnremGeUxwQpdkiYSUXscAx4IpM5ycB8/Nxf3bJRbPUhOVYuiMl3bl+udVZ2OTC3f3Kl93PfML9aHKhqxCRUiI5OdnHjRvn7dq180qVKjngAwYMiPhzBg4c6A0bNvQKFSp4o0aNfMCAAf7pp58evz569GgPfn/8n/fff99T/5+Y4bjuuusylJs0aZI3bdrUY2Nj/YwzzvD58+dnef57773ngC9dujTb+KZPn+4tWrTw2NhYb9CggQ8fPtx3796doUxajLkd77//fr7jyu6d3d1nzJjhHTt29IoVK3rr1q19woQJnpLLOgKFSS4K85ziIFrJRWi7oprZUuATdx+S7twXwKvunmVAZ6Z7JwId3P3sdOcqAG0yFb0J6Av8Atjo7vtyqjM/u6LmZPZXcOMMmHQ+XHhKoaoQkVJqzpw59OvXjz179lC9evWwwxHJIFq7ooY5W+Rh4BkzWwZ8SDATpBHwKICZjQe6ufs5aTeYWTuCVo+6QFUz6wTg7ivc/SiQYU0LM9sOHHb3iK91kd45zaFpDXh8uZILEfmfAwcO8OGHH3LqqacqsZAyJbTkwt2nm1kd4F6gIUFicIG7b0ot0hBomem2t4H0o2PSJo7nPnIzysqXg193gjHzIWkLdGkYZjQiUly89NJLvP3220yfPj3sUESKVGjdIsXNiXSLQLCJWY8noHdTmHxBBAMTERGJkmh1i2jjsgipEguDOsDML+GbH8OORkREJDxKLiLo152gnMETWVf5FRERKTOUXERQo2pwSRt4cTXsOhh2NCIiIuFQchFhQ86Ag8fgmVVhRyIiIhIOJRcRdmpdOLspTFsBh46FHY2IiEjRU3IRBUMT4PuD8OpnYUciIiJS9JRcREGPxtCxHkxJguSUsKMREREpWkouosAsaL3YuAdmrw87GhERkaKl5CJKzmsJTWrAo4mgdcpERKQsUXIRJeXLwY2dYcU2WPZd2NGIiIgUHSUXUXRFO6gdB5M/CjsSERGRoqPkIoriKsANnWDeJvh0R9jRiIiIFA0lF1H2q9OhWixMUuuFiIiUEUouoqxGRfhVR3h7HXy1K+xoREREok/JRRG4oTPElod/FX5HdxERkRJDyUURqFsZruoAr38O32o7dhERKeWUXBSR33YJfk75ONw4REREok3JRRFpVA0GnBpsx75jf9jRiIiIRI+SiyJ0UwIcTYEnlocdiYiISPQouShCzWtB/9bwzCrYcyjsaERERKJDyUURG54A+46o9UJEREovJRdFrG18sKnZkytgz+GwoxEREYk8JRchuLU7/HgE/q3WCxERKYWUXISgXTz8vGWQXKj1QkREShslFyEZ0S1ovZi2IuxIREREIkvJRUg61IN+LeDx5fCjWi9ERKQUCTW5MLObzGyDmR0ysyQz651L2UpmNs3MVprZUTObl02ZAWY228x2mNleM1tqZhdH9SVOwIjuQWLxpFovRESkFAktuTCzgcAEYBzQGVgEzDSzJjncUh44BEwE3sqhTB/gPaB/ap1vA6/nlrSE6bR6cG7zYFrqXrVeiIhIKRFmy8VIYJq7T3X3z9z9FmALMCy7wu6+392HuvsU4Nscytzq7ve7+zJ3/9Ld/wgkAZdG6yVO1G3dg0GdT34SdiQiIiKREUpyYWaxQBdgdqZLs4GeEX5cNWBXhOuMmNPqQ98WMDVJq3aKiEjpEFbLRV2Cbo5tmc5vAxpE6iFmNhw4CXgmh+tDzCzRzBJ37NgRqccW2MgewcwR7ZgqIiKlQdizRTzTZ8vmXKGY2WXAA8A17r4p24e7T3H3BHdPiI+Pj8RjC6VdPFzYGv69Ar4/EFoYIiIiERFWcrETSCZrK0U9srZmFFhqYvEM8Ct3f/NE6ysKt/eAQ8dgcmLYkYiIiJyYUJILdz9CMNCyb6ZLfQlmjRSamV0JPAtc7+6vnEhdRalVbRhwKjyzErbuCzsaERGRwguzW+Rh4HozG2xmbc1sAtAIeBTAzMab2dz0N5hZOzPrRDBmo6qZdUr9nHZ9EPAcMAr4wMwapB61i+qlTsSt3SHZ4ZFlYUciIiJSeDFhPdjdp5tZHeBeoCGwGrgg3fiIhkDLTLe9DTRN9zlt6y9L/TmU4J3+kXqkmQ+cHbHgo6RJDRjUHqZ/Cr/tEnwWEREpacw9IuMnS7yEhARPTAx/wMOWvdDnKbj4FHiwX9jRiIhIaWZmSe6eEOl6w54tIpk0rAbXdoRXP4d1P4QdjYiISMEpuSiGhidA5QrwwAkNbRUREQmHkotiqE5lGHIGzPoKkraEHY2IiEjBKLkopgZ3hvjKcP9C0LAYEREpSZRcFFNVYoOpqcu+g/c2hh2NiIhI/im5KMYGtYdmNeCvH0JyStjRiIiI5I+Si2KsQnn4fU9Y+z28/nnY0YiIiOSPkoti7oLW0LEePLwk2HtERESkuFNyUcyVMxh1FmzeC0+vDDsaERGRvCm5KAF6nQx9msLEZbD7UNjRiIiI5E7JRQlxz1mw9whMWBp2JCIiIrlTclFCtKkbzB55eiWs3xV2NCIiIjlTclGCjOwBFcvD/R+GHYmIiEjOlFyUIPFV4KaEYFnwJd+GHY2IiEj2lFyUMIPPgEZV4c8LIEXLgouISDGk5KKEqRQDf+gFq7fDG1pYS0REiiElFyXQJW2ChbX+uggOHg07GhERkYyUXJRA5Qzu+wls3Qf/Sgw7GhERkYyUXJRQXRvDxafAo0nwzY9hRyMiIvI/Si5KsLvOClox/rIg7EhERET+R8lFCdaoGtzUFWZ+CYu+CTsaERGRgJKLEm7IGXBSdRgzH46lhB2NiIiIkosSr1IM/F9vWPs9PLsq7GhERESUXJQKP28Z7Jz68GLYdTDsaEREpKxTclEKmMGYPrDvCDywKOxoRESkrFNyUUqcUgeuPx2eXw2fbA07GhERKctCTS7M7CYz22Bmh8wsycx651K2kplNM7OVZnbUzOblUK5Pal2HzGy9mQ2N2gsUM7f3CDY3u+d9SNbgThERCUloyYWZDQQmAOOAzsAiYKaZNcnhlvLAIWAi8FYOdTYH3k6tqzMwHnjEzC6LbPTFU7WKcG9vWLU9aMEQEREJQ5gtFyOBae4+1d0/c/dbgC3AsOwKu/t+dx/q7lOAnDYcHwp85+63pNY5FXgK+F00XqA4uvgU6HkS/G0R7DwQdjQiIlIWhZJcmFks0AWYnenSbKDnCVR9ZjZ1zgISzKzCCdRbYpjBn38abGg2fmHY0YiISFkUVstFXYJujm2Zzm8DGpxAvQ1yqDMm9ZkZmNkQM0s0s8QdO3acwGOLl1a1g8W1XvkMlm0OOxoRESlrwp4t4pk+WzbnIlFndudx9ynunuDuCfHx8Sf42OLl5m7QuBrc+z4cTQ47GhERKUvCSi52AslkbaWoR9aWh4LYmkOdx4DvT6DeEqdyhWDti7XfwxPLw45GRETKklCSC3c/AiQBfTNd6ksw06OwFgPnZlNnorsfPYF6S6R+LYPVO/++FDbtDjsaEREpK8LsFnkYuN7MBptZWzObADQCHgUws/FmNjf9DWbWzsw6EYyfqGpmnVI/p3kUOMnM/pFa52DgeuDBonih4uhPZ0NMuWDtCz/RDicREZF8iAnrwe4+3czqAPcCDYHVwAXuvim1SEOgZabb3gaapvuc1uBvqXVuMLMLgL8TTGn9Dhjh7q9G5y2KvwZV4c6e8H/z4PW1MODUsCMSEZHSzryAv86aWUWCFoY4YIe7l4ppFgkJCZ6YmBh2GFGR4nDZy7BxN8z9JdSOCzsiEREpDswsyd0TIl1vvrpFzKyamQ0zsw+APcCXBC0NW83sGzObamZdIx2cREY5g/t/Bj8ehr8sCDsaEREp7fJMLszsdmAjcAMwB7gE6AScQrBo1WiC7pU5ZvaOmbWOWrRSaG3qwtAuwdoXC78OOxoRESnN8uwWMbOXgT+6e667VaR2l/wGOOLuj0cuxKJRmrtF0hw6Buc9B8kOs64JpquKiEjZFVq3iLtfkVdikVrusLtPLomJRVlRKQb+ei58vSfYe0RERCQa8tMtUsfMnjCzrWZ2zMy+N7NFZvaAmXUriiAlcro3hutOh2kr4CMtDS4iIlGQnwGdz5+oQksAACAASURBVAJnAX8BrgV+D5wB9AEWmtk8M2sRvRAl0u7sCY2rw+/fDbpKREREIik/yUUf4DJ3f8TdX3T3fwNHgUEEU1KTgEVmdkoU45QIqhILfz0HNuyGhxaHHY2IiJQ2+UkuNgN1srvg7jvd/Q6CFTAnRDIwia6zmsDVHeDx5bB8a9jRiIhIaZKf5GIC8GQe4yteBnpHJiQpKnefBQ2qwO/mqHtEREQiJz+zRSYCLwCLzWyumd2czX2/BErFSp1lSbWKcP+58OUP8PCSsKMREZHSIl8rdLr7PUB3ggTifoKlvz9NXZ1zNzAKuCNqUUrU9GkK13SAKUmwTLNHREQkAvK9cZm7JwKDzCyGYIXONkANYCfwnrvvjE6IEm339IYF38DI2fDONVA1NuyIRESkJCvwluvufszdE939udRFs15SYlGyVYmFh/vBtz/CWO09IiIiJ6jAyYWUTl0bBXuPvLAa5m4IOxoRESnJlFzIcbf3gLZ14c534YeDYUcjIiIllZILOa5iDPy9H+w+BHe/B3nsaSciIpKtAicXZtbEzBrmcF7JSgnXNh5+1xNmfgnTPw07GhERKYkKkwxsBObmcP4TM/vJiQQk4RtyBvQ6GcbMh692hR2NiIiUNIVJLm4A7s7h/GvAAycUkYSunAXdI5ViYMQ7cCQ57IhERKQkKVByYWbx7j7N3d/IfC31/Gh37x658CQs9avC386F1dvhQW1uJiIiBVDQlotF2l697OjXEq49DR5LgoVfhx2NiIiUFAVNLt4mSDDOSH/SzH5iZh9GLiwpLu7tDa1rw22z4PsDYUcjIiIlQYGSC3e/lWB79ffNrJ+ZdTKzd4D3Af1uWwrFVYBHzoMfD8PtsyFF01NFRCQPhVn++0FgHDAD+AjYC3R096siHJsUE23jYUwfmL8J/pUYdjQiIlLcFXRA58lm9hjwJ4LE4jDwlrtrRYRS7qoOcPEpweBO7Z4qIiK5KWjLxTqgM3Chu/cCLgb+bmb3RDwyKVbMYPw50LQG3DwTdmr8hYiI5KCgycW17t7N3ecAuPt7wNnAMDObHOngpHipGguTLwiWB799lsZfiIhI9go6oPOVbM59AvQiSDIKxMxuMrMNZnbIzJLMrHce5U8zs/lmdtDMNpvZfWZmmcpcbWYrzOyAmW01s2fNrEFBY5PstUsdf/HB1zDpo7CjERGR4igie4G4+yaCBCPfzGwgMIFgcGhnYBEw08ya5FC+OjAH2AZ0BUYAvwdGpivTC3gGeApoD1wKtAOeK9gbSW6u6gCXtIGHl8AHm8KORkREips8kwsza56fitx9lwVOzuezRwLT3H2qu3/m7rcAW4BhOZS/BqgMXOfuq939VeCvwMh0rRdnAt+6+9/dfYO7LwEeAbRqaASZwf3nBOtf3PIOfPNj2BGJiEhxkp+Wi8Vm9oSZnZlTATOrZWbDgDXAJXlVaGaxQBdgdqZLs4GeOdx2JrDA3Q+mOzcLaAQ0S/38IdDQzC5KTXTqAoMIFv+SCKpcAaZcCCkpMPQtOHQs7IhERKS4yE9y8QawB3jLzHaY2Ttm9qSZ/cvMXjSzlcB24FrgNnefmI866wLlCbo40tsG5DQ+okEO5dOu4e6LgasIukGOADsAA67LrkIzG2JmiWaWuGPHjnyELek1qwn/+Hmw/8i974FrgKeIiJC/5OIGgu6Hk4A6wHdATaA5cIxgfENnd+/l7rMK+PzMfx1ZNufyKn/8vJm1A/4J/JmgZeQ8gsTjsWwrc5/i7gnunhAfH1/A0AXgnBZwazd4+TN4blXY0YiISHEQk48y3wDd3f3N1KENo9x9+wk+dyeQTNZWinpkbZ1IszWH8qS75y5gmbunbfu+0sz2AwvM7B53/+bEwpbs3NYDPtkGY+bDqXUhoVHYEYmISJjy03JxP/CqmX1M0EJwg5n1Tp29USjufgRIAvpmutSXYNZIdhYDvc2sUqby3wEbUz9XJkha0kv7bEhUlDP453nQuBr8dgZs1gBPEZEyLc/kwt2nAh2AFwn+gr4emAvsMrP1ZvZa6noTFxfw2Q8D15vZYDNra2YTCAZnPgpgZuPNbG668s8DB4BpZtbBzAYAo4CH3Y/39v8XuMTMhplZi9Spqf8EPnZ3bawWRTUqwRMXw+FkGPxfOHA07IhERCQs+Vrnwt3XuvvfCJb/PguoBnQD/gJsJmhBeLogD3b36cBtwL3AitR6L0hdMwOgIdAyXfk9qc9pBCQCk4CHCJKUtDLTCKa43gysBl5JjTnPGSxy4lrVhkfOh8+/h5HaQVVEpMwy1xB/ABISEjwxUVt+RsKUj+EvC+C27nB7j7CjERGRnJhZkrsnRLre/AzoFCmQGzvDFzvhH0vhlDrQv3XYEYmISFGKyPLfIumZwV9+Bl0aBhucfbwl7IhERKQoKbmQqKgYA1MvhAZV4Tf/hU27w45IRESKipILiZo6lWHaJcHAzuv/A7sO5n2PiIiUfEouJKpa1ILHL4LNe+HGGdqDRESkLFByIVHXtRE81A8++g5+N0dTVEVESjvNFpEicdEpwcqd4z+EelXg/3oHAz9FRKT0UXIhRea3XWDrfnhiOdSJg+Fdw45IRESiQcmFFBkzuO8nwcDOvy2CWpXg6tPCjkpERCJNyYUUqXIGD/aF3YfgnvehZiW4QItsiYiUKhrQKUWuQnl4tD90bgC3zoKF2lJORKRUUXIhoYirAE9eDM1rBruoLtscdkQiIhIpSi4kNDUqwbO/gEbVgkW2krRMuIhIqaDkQkJVrwq8MADiK8N1b8AnW8OOSERETpSSCwld/arwwmXB4M5r34DV28OOSEREToSSCykWGlULEoxqsXDN6/DpjrAjEhGRwlJyIcXGydWDLpLKMTDoVW3VLiJSUim5kGKlaU14+Ypgga1rX4fF34YdkYiIFJSSCyl2TqoeJBiNqgWDPN/fGHZEIiJSEEoupFiqXwWmXwatasON/4W314UdkYiI5JeSCym26lQOBnmeVh+Gz4SnPwk7IhERyQ/tLSLFWo2K8Pwv4OaZ8H/z4Lu98IdewR4lEh532LoPVm4Ppg5v3Qc/HIQfDsEPB+DgMahWMZj9U71icLSsBafXD446lcN+AxGJJiUXUuzFVYDHLoTR8+BfSbBlHzzQF2LLhx1Z2fL1Hpj9FSz6FlZugx0HgvPlDeKrQO1KUCsOOtaHSjGw7wjsPRJsUrdxN8z4Ajy1rpOrQ0IjuKQN9G4CMWpDFSlVlFxIiRBTDsb+NBjk+bdFsH0/PHph0LIh0eEerDfyzldBUrH2++B8y1rQp2nQXdWxHrSLD5KJvOw7ErRyfLItON7fCK9/HqzOekkbuKxtUJeIlHzm7nmXKgMSEhI8MTEx7DAkH177DP7wbjCrZOpF0Lp22BGVLnsOwxufw4urYc3OoAuqayPo1wL6tYQmNSLznCPJQYLx6mfw3gY4mgI9GsPvegbPE5HoM7Mkd0+IeL1KLgJKLkqWjzbD0LfgUDJM+Dmc2yLsiEq+FVvh6ZXw1jo4dCxoRbiqPVx4CtSOi+6zdx2E1z6HfyUG3S0/axYkGe3VkiESVdFKLkLt6TSzm8xsg5kdMrMkM+udR/nTzGy+mR00s81mdp+ZWaYysWb2p9R6D5vZ12Y2IrpvIkWta2P471X/27L9n8uCZnwpmBSHd9fDla/AJdNh1ldw2akwYxDMvBp+dXr0EwsIxmr8pjN8cD2M6hXskHvB83DLTNixP/rPF5HICm3MhZkNBCYANwELU3/ONLN27v51NuWrA3OAD4CuQBtgGrAfeChd0ReAk4EhwDqgPlAE/3uUotaoGrxyBYyaCw8tDvrzHzg32MpdcnckOehemvIxfLULGleD/+sNgzpA1djw4qpcAYYlwNWnwZQkmPoxLPwmGG/Tv3V4cYlIwYTWLWJmS4GV7n5junPrgFfc/a5syg8D/grUd/eDqefuBYYBJ7m7m1k/4GWgpbvvLEg86hYpudzh3ytg3EKoVxn+eV7QsiFZHToG0z+FRxPhu31Bt8Nvu8AFraBCMZx988X3cMfsYMrrxafAn84OWjlEJDJKVbeImcUCXYDZmS7NBnrmcNuZwIK0xCLVLKAR0Cz186XAR8BIM/vWzNaZ2T/NrGrEgpdixyxoUn/tiuAvyCtfhQlLITkl7MiKj4NH4fGPofc0uG9e0Orz9KXw1lXBTI3imFgAnFIHXrsS7jgTZn4JfZ+FDzaFHZWI5CWsMRd1gfLAtkzntwENcrinQQ7l064BtADOAk4HLgNuBs4j6D7JwsyGmFmimSXu2KE9vku60xsEf1ledAo8vASufh2+/THsqMJ18GjQtXDWk/DnBcE00hcGBN1JfZoGiVlxV6E8jOgGbw4Kxn9c95/gnTTGRqT4Cnvpmsz/e7BszuVVPv35cql/vtrdl7r7LIIE4zIzq5+lMvcp7p7g7gnx8RqWXhpUqxjMHnmob7DQU99n4ckVZa8VI31SMXYBtKkLL18OL14GPU8uGUlFZu3i4fUrgymxYxfA7+YE3TwiUvyENaBzJ5BM1laKemRtnUizNYfypLtnC7DZ3fekK/NZ6s8mudQtpYgZXN4OepwEd78HY+bDm1/AX88JmtlLs/1H4NlVQWKx4wD0Ohn+1R26lZIxKFVi4V/94Z9L4e9Lg8Goj10YbHQnIsVHKC0X7n4ESAL6ZrrUF1iUw22Lgd5mVilT+e+AjamfPwQaZRpjcUrqT/XUljEnVYenLoG/94P1u4KpjQ8thgNHw44s8vYehonLoNeTwcDWNnWClornB5SexCJNOYPbesC/LoDPd8LFLwYDP0Wk+AhztshA4BmCKagfAkOB3wDt3X2TmY0Hurn7OanlawBrgXnAWIKkYRrwR3d/KLVMVYKWiiXAGKAm8BjwmbtfkVs8mi1Suu08AH9MbcGoVwXu6AFXtIPyYXcMnqAd++GplfDUJ/DjYfhps2B8whkNw46saKzZAb96A5I9GKB6Wr287xGR/ymVK3Sa2U3AH4CGwGrgdnf/IPXaNOBsd2+WrvxpwCSgG7ALeBT4k6d7CTNrAzxCMLBzF/AGMMrd9+YWi5KLsiHxu6C/fvlWOLUO3N0bftKk5I1BWPdDMPvj9c+DNSv6tYRbugb7fZQ1G3bBNa8HydWTl2jpcJGCKJXJRXGi5KLscIe3v4T7Pwx2+uzcAIZ2Cf6CLs5buSenwPxNwZiKuRugYvmg9eU3naFFrbCjC9fmH4MEY+s+ePwiOKtJ2BGJlAxKLqJMyUXZc/gYvPhpMPjxmx+hRU0Y0gV+cWr+dvksKtv2BXFO/xQ274W6cfDLjsFRp3LY0RUfO/bDta/D+t3w6AVwjvabEcmTkosoU3JRdh1LCRZoejQpWEK8ZiW4sDUMaAtnNAiny2T3oWCfj7fWwcKvgzEFZ50cLIvdtwXEFtNFr8K2+xD88vVge/gnLwlmy4hIzpRcRJmSC3GHD78JWghmfQWHk6FZDbj01GCXzg71ojcA1D3oovnwG3jnq+DnsZRgxsvFp8DA9tCsZnSeXdrsOggDXw1ao577RdkZ3CpSGEouokzJhaS393Dwl/xrn8Hib4OV2apXhDNPCo6ERsFql5UrFK7+g0eD5vuV22DpZljyLWzZF1xrUgP6two26upQr+QNNi0Otu2HK1+GHw7B9MuCBbhEJCslF1Gm5EJysmN/kGB8+A0s+jZoYUjTqGowmLJ5LahVKVjkqUqF4Gc5YN+R4PjxCOw5FNz71a5g7ESaunHBgl89ToLujaF1bSUUkfDtj3D5y3A0GV6+QoNeRbKj5CLKlFxIfn29B1ZtDxbm+mpX8HPD7mAqZE7KW7CVeZMaQYtHy9rBzzZ1gp9KJqLjq11wxcvBzJrXB0IDbWEokkG0kotiNCZepGRoUiM4MkvxoLtj/9FgFdBkDxKK6rHB7BMlEEWvZS145hdw5Svw6//AS5cH+8+ISHSV8PUJRYqPchZ0h9SrEgy+bFkr2PMiroISizC1j4fJFwQzSIbPDLpJRCS6lFyISKnXpymM+1mwCNm972u7dpFoU7eIiJQJgzoE01MnfhR0aw3vGnZEIqWXkgsRKTN+d2aQYPxtUbCGyCVtwo5IpHRSt4iIlBlm8MC5wZTf38+BT7aGHZFI6aTkQkTKlIox8K8LIL4K3Dgj2LtFRCJLyYWIlDl1KsPjF8LeIzBkBhw6FnZEIqWLkgsRKZPaxsPf+8GKbXDXXM0gEYkkJRciUmad1wpG9oDXPofHPg47GpHSQ8mFiJRpI7oFm8TdvzBYB0NETpySCxEp08zgwb5wal24ZWbGjelEpHCUXIhImVe5AjzWH5xggOeBo2FHJFKyKbkQEQGa1oRHzoPPd8Kd72qAp8iJUHIhIpLq7GbBKp5vfgFPLA87GpGSS8mFiEg6w7vCz1vCuIWw6JuwoxEpmZRciIikYwYP94NmNeHmmbBlb9gRiZQ8Si5ERDKpGguPXRis3DnsbTiSHHZEIiWLkgsRkWy0rg0P9IXlW+HPH4QdjUjJouRCRCQH/VvDjWfA0yuDVTxFJH9CTS7M7CYz22Bmh8wsycx651H+NDObb2YHzWyzmd1nZpZD2bPM7JiZrY5O9CJSFozqFWzRftdc+GxH2NGIlAyhJRdmNhCYAIwDOgOLgJlm1iSH8tWBOcA2oCswAvg9MDKbsrWAp4G5UQleRMqMmHIw8XyoXhF++xbsORx2RCLFX5gtFyOBae4+1d0/c/dbgC3AsBzKXwNUBq5z99Xu/irwV2BkNq0XTwBPAYujFLuIlCH1qsDkC2DzXhg5C1K0wJZIrkJJLswsFugCzM50aTbQM4fbzgQWuPvBdOdmAY2AZunqvgloAIyNVLwiIl0bwb294d0NMOmjsKMRKd7CarmoC5Qn6OJIbxtBYpCdBjmUT7uGmZ0GjAaucfc8J4+Z2RAzSzSzxB071JkqIrm7/nS4pA08tBg+0A6qIjkKe7ZI5sZFy+ZcXuUB3MwqAi8Cv3P3Dfl6uPsUd09w94T4+Ph8BSwiZZcZ3H8OnFIHRrwD3/4YdkQixVNYycVOIJmsrRT1yNo6kWZrDuVJvach0A54MnWWyDHgPqB96ud+EYlcRMq0yhXg0f5wLAWGvRUstCUiGYWSXLj7ESAJ6JvpUl+CWSPZWQz0NrNKmcp/B2wENgOnAZ3SHY8CX6b+Oad6RUQKpEWtYInwldvhvnlhRyNS/ITZLfIwcL2ZDTaztmY2gWBw5qMAZjbezNJPJX0eOABMM7MOZjYAGAU87IGjqbNIjh/AduBw6ud9Rft6IlKa9WsZbHI2/VN4flXY0YgULzFhPdjdp5tZHeBegi6N1cAF7p42TKoh0DJd+T1m1heYBCQCu4CHCJIUEZEid0cPWJ3aetGmLnRpGHZEIsWDuWvCNkBCQoInJiaGHYaIlDC7D8FFL8LhYzDjqmBNDJGSwsyS3D0h0vWGPVtERKREq1kJpvSHHw/DTdpBVQRQciEicsLaxsPfzoWPvtMOqiIQ4pgLEZHS5OI2weyRqR9D+3gY1CHsiETCo5YLEZEIGdULejeBe98PWjFEyiolFyIiERJTDiadD42rwdAZsFkreEoZpeRCRCSCalSCxy+GQ8kweAYcOBp2RCJFT8mFiEiEta4N/zwPPtsBv5sDmvEvZY2SCxGRKDineTAG46118M9lYUcjUrQ0W0REJEp+2wXWfg8PL4FmNYPt2kXKArVciIhESdoW7d0aBd0jH20OOyKRoqHkQkQkiirGwJQLgxkkN86AjbvDjkgk+pRciIhEWa04mHZJ8Odf/yfYj0SkNFNyISJSBJrVhKkXwrd7YciMYKMzkdJKyYWISBHp2hge7AtLN8NtsyA5JeyIRKJDyYWISBG6pA3c2xve/hLum6c1MKR00lRUEZEiduMZsPMAPJoEdSvD7T3CjkgkspRciIiEYFQv+P4A/GNpkGD8smPYEYlEjpILEZEQmMH958IPh+D/3odaleDCU8KOSiQyNOZCRCQkabuoJjSCW2fB7K/CjkgkMpRciIiEKK4CPHkxnFYPbnob5m4IOyKRE6fkQkQkZNUqwlOXQpu6MPQtmL8p7IhEToySCxGRYqBGRXjuF9CqFtz4X/jwm7AjEik8JRciIsVEzUrw3IBgNc/fvAkLvw47IpHCUXIhIlKM1I4LWjCa1oBfvwmzNMhTSiAlFyIixUx8FZh+ObSPh2FvwWufhx2RSMEouRARKYZqVoJnfwHdG8Pts+DpT8KOSCT/lFyIiBRTVWPhyUvg3Obwf/Pgn8u0F4mUDKEmF2Z2k5ltMLNDZpZkZr3zKH+amc03s4NmttnM7jMzS3d9gJnNNrMdZrbXzJaa2cXRfxMRkeioFAOP9odfnAoPLYbfvwtHksOOSiR3oSUXZjYQmACMAzoDi4CZZtYkh/LVgTnANqArMAL4PTAyXbE+wHtA/9Q63wZezytpEREpziqUh7/3g9u6w8tr4FdvwJ5DYUclkjPzkNrYzGwpsNLdb0x3bh3wirvflU35YcBfgfrufjD13L3AMOAkz+FFzGwZsMDd78gtnoSEBE9MTCz0+4iIFIVXP4M734UmNWDaJcFPkcIysyR3T4h0vaG0XJhZLNAFmJ3p0mygZw63nUmQJBxMd24W0AholsvjqgG7cohjiJklmlnijh078hO6iEioLmsbDPTceQAunQ6LtNiWFENhdYvUBcoTdHGktw1okMM9DXIon3YtCzMbDpwEPJPddXef4u4J7p4QHx+fn7hFRELX4yR4Y2Awo+Sa1+HRRA30lOIl7Nkimf9zsGzO5VU+u/OY2WXAA8A17q6V+kWkVGlRC94cBBe0gvEfwm/fgh8Phx2VSCCs5GInkEzWFod6ZG2dSLM1h/Jkvic1sXgG+JW7v3lioYqIFE9VY2Hi+XDfT4LdVC96Adaoh1eKgVCSC3c/AiQBfTNd6kswayQ7i4HeZlYpU/nvgI1pJ8zsSuBZ4Hp3fyVSMYuIFEdm8JvO8OIAOHAMLpkOjyZBckrYkUlZFma3yMPA9WY22MzamtkEgsGZjwKY2Xgzm5uu/PPAAWCamXUwswHAKODhtJkiZjYIeC71/Adm1iD1qF2E7yUiUuS6NoZ3roafNoPxC2HQq/D1nrCjkrIqtOTC3acDtwH3AiuAs4AL0o2PaAi0TFd+D0FLRSMgEZgEPESQpKQZCsQA/wC2pDtei+a7iIgUB3Uqw2P94aG+sGYnnPccTP9Ugz2l6IW2zkVxo3UuRKQ0+fZHuGMOLPkWep4EY38GLWuFHZUUN6VqnQsREYmuk6rDCwNg7E9h9fagFePBxXDoWNiRSVmg5EJEpJQqZ/DLjvDer6B/a3hkGfR9FuauV1eJRJeSCxGRUi6+Cvzj50FLRoVycMN/gwGfy7eGHZmUVkouRETKiJ4nw6xr4M9nw5c/BMuHD3sLNmS7QYJI4Sm5EBEpQyqUh1+dDvOvD3ZZnbcJznkGbp8FX3wfdnRSWii5EBEpg6rGwu09YP51cP3pMPPLYDzGkBnwibpL5AQpuRARKcPqVYH7+sCiG2BEN1j8LVw8HQa+Am+tg6PJYUcoJZHWuUildS5ERGDfEXhuFTz9CXy7N0g+ru4AV3WABlXDjk4iLVrrXCi5SKXkQkTkf5JTYN5GeHolzN8UTGvt3QR+cSr8vCXEVQg7QomEaCUXMZGuUERESr7y5eCcFsGxaTe8+Cm88TncOguqVIDzW8HFbeDMkyC2fNjRSnGjlotUarkQEcldisOyzfDa5/D2Oth7BKrFBpulndcK+jQNBopKyaFukShTciEikn+HjsGHX8M7X8G7G+CHg0ELRkLDIMno3RTa1g26U6T4UnIRZUouREQKJzkFErfAnPWwYBN8nrpeRt24YOGubo2hWyNoXUfJRnGj5CLKzCzHL8LdMbMsf86uTHZlc/uZ7vlZ6s6rzsz3ZRdbbrFnvpZdLOmfk9P7pr+e3ee83iH9czM/Iz9ye3ZBY8vpu8jpe8nPu2YXb3bvntc/x+yeldMzMteZvp68/v3MTubvILs/Zxdn5ufk9T3kVG9eceRUJv3nnOLP7z/jvP59zM89+Xnn/MSe1/3Ruie/5bbtg4XfBANBl3wL2/YH52tWClo2OjeATg3gtPpQo2KejzqhWIqTSMaZ13//BahHAzpFRKT4q18VLmsbHO7wzY/BWI1l30Hid0E3SpoWNaFDPWgfD23joV3dYC8UKdmUXIiISNSYQZMawXF5u+DcnsOwahus2AqfbAsSjje/+N898ZWhde2gG6V17eBoVRvqxAX1SfGn5EJERIpUjYpwVpPgSLP7EHy2E9bsCH6u+x5eWQP7j/6vTPVYaF7r/9u7+xg7qvuM49/Hsb1m2TUGXN4CzhLRJgS7ISkloQVsoKh5EaWCtmnSqLgqEIJEkjopTUtFUaQEIREISokozh9OhBITKRFKSAkJjoAoUIRpIkpfCI1NI7ANmNisl32zyekf59zs+HLnvthz78ysn480Gt+ZM7Pn7O/O3p/PPTMHxpbFHo8Vy2DF0rgvBCceVeLkwszMSrdsSXxmxlknzm0LAbZPxAnVtuyCLbvjDK5PbINvPw3Z0QanfglOXBqXN47OrY8fhRNG4pNGF/l5HAPj5MLMzCpJghNG47JmbP990/vguXH4xStwwSfiI8qf2xO3/WRH7AnJWqD4dcvxI3FMyLGHx+W4lHgcc3jcf+RhvqOlCE4uzMysdpYsjOMwTjkqvr5+9f77J2Zh2564bJ+A7Xtg2wTsmICtu+MEbeMzrz/vogWwfHhu+Y20HD0cx3wsT+ujh+HIJe4NyePkwszM5p2RxfBbR8clz+ReePHVuLw0mfl3ev3SZBz/sXMS9v2q9TmWDsVk46i0HLkkrdO/j1wSv/JZtiSW3/vaoZGQOLkwM7ND0vCiODh0bFn7cr8KsZdj5yS8PJXWk/GppC9Pza2fG4cnX4Bd7gbNHwAACxBJREFU0zCbM1X9Kf8cE59lQ3BESjqOGJpbH5FZL10893rpEIwOwcIFxf8e+sHJhZmZWRsLNNf7cEoX5UOIvSK/nIbdU3H8x65puBhY9254ZRp2z8CuqbjeMRGTl93TsDenh6Th8EUx0Wi47J6UeCyOycfoYhgZionJ6FBMZEYWp+1pGXpD/++scXJhZmZWIAkOXxyXk5buv+/j78o/LgSY2hcTjVem4/NAxmfm1uOZ14+lY16egmd3x0nk9szm95hkLVwQk5TRPk4y5+TCzMysAqT4Vc3wongXSzufT+t7P7j/9pl9KdGYgYm9aT0bl/FZeDUte2bjM0Qe6UtLnFyYmZnNG0ML47J8uLvyt/apHqUODZF0taStkqYlPSHpnA7lV0l6SNKUpOclXa+mWYskrU7nmpa0RdJV/W2FmZmZZZWWXEj6AHAb8DngHcTemfskrcgpvxT4AfAC8LvAx4C/BdZlypwM/Gs61zuAG4EvSrq0fy0xMzOzrDJ7LtYBG0II60MI/x1CuAbYDnw0p/xfAMPAZSGEp0II3wRuAtZlei+uAraFEK5J51wPfAX4VH+bYmZmZg2lJBeSFgO/A3y/adf3gd/LOews4EchhKnMtvuBE4CxTJnmc94PnCFp0cHU2czMzLpT1oDO5cAbiF9xZL0A/EHOMccBz7Uo39i3Na0faFFmYfqZ27M7JF0JXJlezgBPtfrB2WEdyrk5uLG9VdlO67xzd3POdudr2rZc0s5u29XunP1oQ7uf0aW27eu2bp1ed/O77uX90k190rblwM7msu1+Rqt93bw/O52nm/dNq22d3p/ktK+berSrQ6/x7PZ90cPP7rlt3dS9m+P7cUyLcvu1ryw9/r3oRaHtK7qeBZzvLUXUo1nZd4uEptdqsa1T+ebt3ZSJG0K4E7gTQNLmEMIZbWtbY25fvbl99TWf2wZuX91J2tyP85Y15mIn8BqxpyHrGF7fm9GwI6c8mWPyyuwDXj6gmpqZmVlPSkkuQgizwBPAhU27LiT/mR6PAudIWtJUfhvwbKZM89cqFwKbQwh7D6bOZmZm1p0y7xa5BVgr6XJJp0q6jTg48w4ASTdK2pQp/zVgEtggaaWkS4BPA7eEEBpfedwBnCjpC+mclwNrgZu7qM+dxTSrsty+enP76ms+tw3cvrrrS/s097k8eJKuBq4FjicOpvybEMLDad8GYE0IYSxTfhVwO3AmsIuYTHwmk1wgaTXxoWOnEXs1bgoh3DGI9piZmVnJyYWZmZnNPzWZGd7MzMzq4pBJLtSHeUyqQNLfS3pc0riklyR9R9LKDseMSQotlvcMqt7dknRDi3ru6HBMLWIHIOnZnFh8N6d8pWMn6VxJ306/9yBpbdN+pZhuS/F5UNJpXZy39DmD2rVN0iJJN0l6UtKrkrZL+ppypjPIHLcmJ55v7XuDXl+XTrHb0KKe/9bFeUuPXapHp/a1ikOQdHubc1Yift18Dgz62jskkgv1YR6TClkDfIn4ZNPzibfdPiDpqC6OfQ9xvEtj+WGf6niwnmb/eq7KK1iz2EGsY7Zt7yQ+k+UbHY6rauxGiOOnPg5Mtdh/LfBJ4Bpi218EfiBpNO+Eqs6cQe3aNkyM3WfT+mLgJOB7krp5ntBp7B/PZwqqcy86xQ7iQwqz9XxfuxNWKHbQuX3HNy0Xpe2drkUoP35r6Pw5MNhrL4Qw7xfgMWB907ZngBtzyn8UGAcOy2z7R+B50jiVqi7EC+g14KI2ZcaIH2BnlF3fLtpzA/BUD+VrG7tU1+uA3cDwPIjdBLA281rEp+Rel9l2GLAH+Eib89wEPNO07cvAo1VpW06Zt6VYrWpTZk0qs7zseHVqH7ABuLfH81Qudj3Ebz3wdIcyVY3ffp8DZVx7877nQv2bx6SqRok9Uru6KPstSS9K+rGkP+lzvQ7Gm1NX5lZJGyW9uU3Z2sZOkoC/Bu4KIUx2KF6X2GWdTHzI3a+vxRSnh8m/FqG+cwYtTetursXN6auUTZLO62elDtLZ6X33M0nrJR3ToXwtYydpBPhzYoLRjarFr/lzYODX3rxPLmg/j0nz0zwbjssp39hXZbcBPyU+UCzPBHGm2D8jdmtuAu6W9OH+V69njxGfVfJe4Ari7/8RSUfnlK9z7C4k/hH4cpsydYpds8bvv5drsXFcq2MacwZVTvpPzeeB74QQmudEymrMBH0pcAnxK8BNks7tfy179j3gL4ELiN3rZwI/lDTU5pjaxS75EDBEnFW7narGr/lzYODXXtlziwxSP+YxqRRJtwBnA2eHEF7LKxdC2En8w9ewWXFirGuBu/pby96EEO7Lvk4DyLYAlxEfxNbysKbXlY9dcgXweAjhp3kF6hS7Nnq9FvOOabW9dGmMxV3AMuCP2pUNITxN/EBqeFTSGDGBfLhPVTwgIYSNmZf/IekJ4P+A9wPfando0+vKxi7jCuCeEMJL7QpVMX4dPgcGdu0dCj0X/ZrHpFIk3Qp8EDg/hLDlAE7xGPCbxdaqeCGECeA/ya9r7WIHkLqXL6b7btisWsSOGBvo7VpsHFeLOYNSYvF14LeBC0IIB1K/WsQzhLCNOFN1u7rWJnYNkk4HzuDArkUoMX5tPgcGfu3N++Qi9G8ek8pQfHT6h4hvqP85wNOcTtOU9FWUYvJW8utaq9hlrAVmgI0dyrVSi9gBW4l/rH59LaY4nUP+tQg1mTMofQd9NzGxOC+E0PaW6TZqEc/UY/ZG2te1FrFrciXxb8UDB3h8KfHr8Dkw+Guv7FGtAxo5+wFgFrgcOJX4fdQE8Ka0/0ZgU6b8ESkQG4GVxO/SxoFPlt2WFm27PdXtfGKG2VhGMmWa23dZehOeCryF2IU3S3z8eultamrfzcBq4liEdwH3pvbWPnaZOgv4GU13NNUxdsRR6qenZRK4Pv17Rdr/dykel6T4bCQmfqOZc3wV+Grm9cnAq8AXUrsvT22+tCptI37FfA/xrqR3Nl2Lh7Vp2yeAPyb+T/e0FO8AXFKl2KV9NxMH+I0R75J4lNhzUfnYdfPeTGWGgVfI3FXRdI5Kxo/uPgcGeu0NNLhlLsDVxGx0htiTcW5m3wbg2abyq4jfmU0Ts9B/ooK3MqY3cqvlhrz2ET+g/iu9acaBzcCHy25LTvsaF8As8Q/3N4G3zYfYZep7XorZmS321Sp2zN2a17xsSPtFvL14e4rPQ8DKpnM8CDzYtG018O/p+t0KXFWltjF3i3CrZW1e24hjZf6X+NyFXwI/At5XtdgRb1u8n/hshFniWIsNwEl1iF03781U5q+IXf4n5JyjkvFr8967IVNmoNee5xYxMzOzQs37MRdmZmY2WE4uzMzMrFBOLszMzKxQTi7MzMysUE4uzMzMrFBOLszMzKxQTi7MzMysUE4uzMzMrFBOLszMzKxQTi7MbOAk/amkGUlvymy7TdLPJR1bZt3M7OD58d9mNnCSBDwO/CSEcIWkTxHnafj9EMIz5dbOzA7WwrIrYGaHnhBCkPQPwHcl/Ry4jjhVtBMLs3nAPRdmVhpJjwBnAheFEO4ruz5mVgyPuTCzUkg6H3g7cSroF0qujpkVyD0XZjZwkt4OPASsA94PjIQQ/rDcWplZUZxcmNlApTtEHgH+JYTwGUkrgSeJYy4eLLVyZlYIJxdmNjCSjgJ+DDwcQvhIZvvdwIoQwlmlVc7MCuPkwszMzArlAZ1mZmZWKCcXZmZmVignF2ZmZlYoJxdmZmZWKCcXZmZmVignF2ZmZlYoJxdmZmZWKCcXZmZmVignF2ZmZlao/wflkk8vIzttMAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot fitted pdf\n", "yMin = 0.\n", "yMax = f(0., MLE)*1.1\n", "fig = plt.figure(figsize=(8,6))\n", "xCurve = np.linspace(xMin, xMax, 100)\n", "yCurve = f(xCurve, MLE)\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)$')\n", "plt.figtext(0.6, 0.8, r'$\\hat{\\theta} = $' + f'{MLE[0]:.4f}' +\n", " r'$\\pm$' + f'{sigmaMLE[0]:.4f}')\n", "plt.figtext(0.6, 0.72, r'$\\hat{\\xi} = $' + f'{MLE[3]:.4f}' +\n", " r'$\\pm$' + f'{sigmaMLE[3]:.4f}')\n", "plt.xlim(xMin, xMax)\n", "plt.ylim(yMin, yMax)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEiCAYAAADJdLQBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3gc1dX48e9R726yZMmyrOLebdxtbIHpIQFCAiaBQOiEFNIbv7xOAiG8IQRC8oYQSGghBAiE0EwxFti4YdybXFRsy5LcZMvq7f7+mBFZr1fSquzMSjqf59nH3in3nrk72rMzc+eOGGNQSimlnBLidgBKKaX6Fk08SimlHKWJRymllKM08SillHKUJh6llFKO0sSjlFLKUZp4eiARKRSRpW7HoZRSnaGJJ0iJyHwRWSIi/R2q7+sicoMTdXWViHxZRLaISK2IFInIz0Uk3I/1BorId0UkV0TKRKRCRDaIyB0iEupjeRGRb4nIbhGps//9poiI13I3iIhp5TWiO7fdo87pIvK+iFSKyDEReVpEkjqw/ggReVVETtrt8KqIZPtYrrCV7XrWa7kxInKf3Z4nReSwHd/5PsrMbaO9jIgM7VyrdJ6IXCAia0SkRkRKReT3IhLXgfX9+jza2OZ7vJabYcew1S7zkIi8LiLTfZTZ2mdkRKShcy0SWGFuB6BaNR/4H+BJ4IQD9X0dKLXrC1oicj1WjG8AjwATgbuBYcCN7aw+F/g18Kb9bw1wPvB/wALgGq/lfwYsAZ4B/hdYCDwM9Ad+4aP8JcA+r2ml7W1TR4nIOCAXKAR+aMfzXWCaiMwwxtS0s/4QYAXQAPzcnnwXsEJEphpjyrxW2QL8xmtavtf7m4GbgH8BjwFRwHXAOyJyuzHmzx7L3gs87rV+KPAXYLcxprit+LubiJyLtU+sxWqHdOA7wHgROc+0c5d9Jz6P94G/eU3b6vX+h1jfAS8Cv7fLvBVYKyKfNca86bHsXYB3khyIta++3VbsrjHG6CsIX8CPAANk+JhXCCzt5vp2Ablub3c7MUYCh4H3vKYvsdtqajvrZwLDfUz/k73+JI9pKUAt8ITXsk/a04d4TLvBXn+2Q+3wqt0OAz2m5dgx3OXH+g8DdUC2x7SRQD3wUGf2NWA6EOc1LQLYDBwBQtpZ/2I7/h91sk2e7Oz+a8e4A4jw8Zle3p2fhz3tUT/KnOsZjz1tAHAI2ODH+nfYdS12Yp/s6EtPtQUhEVkC3Ge/LfA4bM7xWm6miKy0Tw8cEJG7fJQVISL/T0Ty7NNFpSLyqHicwhORQmA0sNCjrkKP9X8hIutFpNyua62IfC4wW9+mc4DBwB+8pv8R64/sqrZWNsYUGGOKfMx62f53nMe0y7AS3SNeyz5iT7/MVx0iEu/rtF13EZF4rC/p54wxx1umG2NygW3A1X4UcxXwljHm06MzY8werF/HPtcXkXARiW2tQGPMemNMpde0euB1IBFo7zTgtVif4XN+xN9tRGQMMAl4zI63xbNYZxrabM/Ofh4iEiUi0a2Va4xZ5RUPxphyYDmn76etuRY4hZUUg44mnuD0MvBP+//fxjplcR2w02OZTOA14COs0wL7gN+JyAUtC9jXIl4BfgK8A3wDeMou613573WRu4ASrKOelrpaklgCcDuwGvgp1mmtcODfInJhexsiIpEikujnq739cZr97zrPicaYI1infqadsYZ/Uu1/j3rVVYd1msnTRqwjA191vQtUADUi8pZ9Cqa7TcRq/3U+5q0FprTVjiKSCgxpY/0hIpLiNX0hUA1U2j9wfuDHZ9UiFWgETrYRUyxWIv/QGLPfz3K7S2v7VCPwCe3vU535PK7Fas9q+wdhe6eIPaVy+n56BhHJxDpietm0c9rVLXqNJwgZY7aIyCasX0v/NsYU+lhsFHC+MeY9ABH5K7AfuAUryYB1zeJiYJExZnnLiiLyAdY1ksXAM8aYf4vIr4EyY8xpF42BcmCYMabOY/1HgE3A92j/HPI1nHk+uzWZWKd2WtOSIHxdNznkMd9vIhIBfN9ef4VXXWXGmGbP5Y0xzSJS5lVXNfapHqwv2GlYPwZWich0Y8zejsbVhpZ6S3zMO4R1bWUAcKyT67cs0zJ/C7ASyAMGAV8B7geGA3e2FahYnRUWA/9p5wvwCiAW61qa09prj1ldXN/781gFvAAUYJ3OvQN4QkQGGWO8r6OdRkTmYf0I+F07MX3Z/teN9vSLJp6ea09L0gEwxtSJyBogy2OZq4DdwFYRSfSYvg6oBM6lnZ3TGNMENMGnX9JxWEfKH+DfaZ23sS7g+6O9C/HRQIN3MrDV2vM76vfAeOAyz+Rql1Xne5XT6zLGvID1ZdLi3yLyJtbR6BKsX7jdpaVeX7HVei3T5fWNMaedUhWRv2Edad8hIg8bY3b7qkREYrDapBbrqL0t19rLvdTOci1lhwP9vCZHAuFe+znASWNMWz272muP9vapjrbnPM8F7B+MHwNLROQxY4zPI0O7h9xzQBH/7RDSmi8DxVin5YKSJp6ey9e1inKs89UtRmFduznSShl+db8VkZuxvjzGAp5didt9poYxpgTfvwY7owbry0WMfQXVQ5Q9328i8hPgNuBuY8x/fNQV2cqq7dZljFkjIiuB89qJIRTrupWn497n973iopXYoryW6fb1jTFGRH4LfAbrh8sZicdODC9iJfRL2jp9JiLJWG30Smtfuj7Mo/UvVe99/RysI9HWtNce7e1TXW3PBhF5GPgrMAc44/48+zrSm0A8sMAYU9FaeXZ36zHAA638QAsKmnh6rqZWpnsmhhCs3jrfamXZNs8VA4jINVjdXP+DdYrlMNY5+68CX/Jj/WjO/HXamiP2EVZrWhJYCv89LdQilTO7MrcV1x1Y3XofMsbc20pd54lIiOcfsH2+PtlH/b4UYX2ZtGUY1mkXT219WXq2gbdUrF/Z5W3U19760P62tfzoGeg9w26fp4ELgauNMe+3U9Y1WF2pvU/xtmUzZx5Ffx/r2tV3fSzbFs/28P4cUmm/Lbr6eUDb7RmNdYQ5FuvU+rZ2ymo5uu5IezpOE0/w6o4n9O0DzgLe9+PXT2v1XY114f5yz6MMEfmqnzFcTfdd4/nE/ncm8G+PWBKxTjG+6E8lInIdVk+4Z7CuxbRW181YR5CbPKZPxeomvMGPqrJo/WizRSlnfom29WW5Fev+m5mc2QNsFrCprc/aGFNsX6Oa6WP2LKDUPkptS8vpXF/b9ijWdZ1bjDH/aqccsE4LHcP6Re8Xu3fXe57TRORaINLz9LOfPPepVR7lhWFdq2vvGmaXPg+bz/a0jxxfwuoo8DljzCrvFb2WD8X6e9tqjGkv4brL7f7c+vL9wrqh0+e9KbRybwXWBe5Cj/fX2WV808eyYcAAj/frgc0+lvsXVgIL8ZiWhXVB3fixHSlYp1L8eUW1U1YU1h/nu17Tl9jbOc1jWgzWKYdEr2UvxzpiexUIayfuOuBxH21cB6R4TBvoY/1FdkxPtLVNndw3XsM68vT8/HLs+r7jtewYIN1r2iO0fh/P7z2m9QfEx36zzG7DTK95D9gx/MDP7RhtL/9/3dAmT9L5+3i20vp9PJ/3mBZut2eK1/p+fR6t7CexWL1VT+JxHxTW2YoXsM5s+HUvDtZRpt/t7+bL9QD01coHA7PtnegtrJ5Ei4Eke14h/iWeEKzu1AYrgXwTq0v1w1gXHxd7LPso0Iw1WsI1wGft6dfb67+Odef0z7G+/DfhR+IJQLvcaMfzGtYRye/tP86nvJZr+cNf4jFtBv899XEz1mkJz1eWVxm/sMt4Cuuu/Kft9z/3Wm4H8Hesu81vw7ohtR7raGZYANpgAlCFdZ/InVjd5Y/bccR4LWvw+kLGSqqlWKd4vm2/9tvTPBPqDVhfyr+2t+tHLZ87cK9Xmd+0p2/w0a7XArE+tuOX9jpzu6FNnvTezg6se769D6209/F7sK7LLMcj8QIZdrxPdubzwPqBtMou/xbg/2H9qDNYR4ieZT5oT3/HV3u2sh3P2NuR5vTfZYfb3O0A9NXGh2PdN7Pf3pkMkGNPL8SPxGNPC8W6J2eT/cd0AutUzv3AUI/lhmBdxzlp11XoMe/7WKfbau0/rmvtPyLjUrtcZ38h1gEH7C8w77u8czgz8dxgT2vtdYNXGYL1pbzXrmuv/d77KOAerC/ccqyEsx9r2Jihgdh+u86Z9hdjlf0l9yweoyl4LHdG4rGnj8JK3hX26zVgpNcy07CODA/Y23/K/uL8civ7Xlttm+FjnX3Avm5qjyd9bWcH1r8Iq7dnDVCGdZNyvNcyGfhIPP5+HlgJ7m2s60L19t/ie8DFPsrLbas9fSwfY38+7zv1d9iVl9hBK6WUUo7QkQuUUko5ShOPUkopR2niUUop5ShNPEoppRylN5C2IzEx0WRkZHSpjKqqKmJjWx1Rvs/zbJ9Dp/wZEKBviWyOpC7kzKHAEg9Y404eHTbI6ZCCSmvt40tqfIfHke3x3Pr++eSTT44aY7yHgwI08bQrIyOD9evXd6mM3NxccnJyuiegXsizfZbkLnE1lmA0unI0eXF5Z0y/4a4nAXjyoRucDSjItNY+vizJWRLYYIKQW98/IuJrPElAT7UppZRymCYepZRSjtLEo5RSylGaeJRSSjlKE49SSilHaeJRSinlKE08SimlHKWJJ0D2Hq7kF6/toL4xaB97rpRSrtAbSAPkwPFq/vpRATMzBxDldjBKKRVE9IgnQBaMGsyQhCie//iA26EopVRQ0cQTIKEhwlXT0/hg9xGO1ejpNqWUaqGJJ4C+OH0YACuKG12ORCmlgocmngAaNjCG+SMSWXGwkaZmfcS4UkqBJp6Au3rGMI7VGlbuPep2KEopFRQ08QTY+eOSiQ+Hf3683+1QlFLKb//z6jaW7SwLSNmaeAIsMiyUuUPDeHdHGUcr/XtYlVJKuWlXaQVPrS5i35HKgJSviccBC9LCaWgyvLKh2O1QlFKqXX9fs5+IsBC+eNawgJSviccBQ+NCOGv4AJ7/eD/GaCcDpVTwqqxr5OUNB7l0UgoDYiMCUocmHodcPWMY+45Usb6o3O1QlFKqVf/eWExVfRPXzh4esDo08TjkMxNTiIsM4/l1OpKBUio4GWN4dk0R41MTmDqsf8Dq0cTjkNjIMD47OZU3th6iorbB7XCUUuoMnxSVs6v0FNfOHo6IBKweTTwOWjxjGLUNzfxn0yG3Q1FKqTM8u6aI+MgwLpuSGtB6NPE4aFJaP8YMieefOnCoUirIHKus482tpVx5VhoxEYF9cIFjiUdEloiI8XqVesz/pYjsEpEqESkXkWUiMterjFwfZTzfTr03+FjHiIjjTysQERbPGMbW4pNsKz7pdPVKKdWqF9YfpL6pmS/PSg94XU4f8eQBKR6viV7z7rSnzQcKgKUikuxVxt+8yrjNj3qrvdZJMcbUdn4zOu+KqWlEhIXwwno96lFKBYemZsPf1xYxO2sgI5PjA16f04mn0RhT6vE60jLDGPOsMWaZMSbfGLMd+A4QD0zxKqPaqwx/Dh2M1zql7a8SGP1iwrlkwhBe2VhMbUOTW2EopdSnPtx9hIPlNQHtQu3J6cSTJSLFIlIgIs+LSJavhUQkArgVqAA2ec1eLCJHRWS7iDwgIv6k52gRKRKRgyLyuohM7eJ2dMnVM9I5VdvIW9tK3AxDKaUAeGZNEYPjI7lg3BBH6nPy0ddrgRuAXUAScDewSkTGG2OOAYjIpcDzQAxQApxvjPEcpe45oAg4BIwH7gMmA+e3UW8ecCOwGesI6lvARyIy2Rizx9cKInIrVuIjOTmZ3NzcTmzuf1VWVp5WhjGG5Bjh0Xe2MuDk3i6V3Rt4ts/oytHuBhOEIpsifbZLdFMMoG3WWvv40tW/5Z7I+/vH25HqZpbvquHS7HBWrfzQkZgcSzzGmLc834vIGiAfuB540J68HOvUWiJwC/CCiMwxxpTYZTzmUcRWEckH1orINGPMhlbqXQ2s9qh3FdZR1DeAb7ayzmPAYwDTp083OTk5HdtYL7m5uXiXcb3s5X+X5pE+fjpZg+O6VH5P59k+S3KXuBpLMBpdOZq8uLwzps8JrQbwOa8vaa19fLkm55oARxN8fH3/ePrfpbsQ2cePv3g2qf2jHYnJte7UxphKYDsw0mNalTFmrzFmjTHmJqABuLmNYtYDTZ5l+FFvk72e3+sEwhempREaIvxTOxkopVxS19jEPz8+wKKxyY4lHXAx8djdmcdgnVJrTQgQ2cb8iUBoO2V41yvApI6sEwhJCVGcOyaJf31ykIamZjdDUUr1UUu3lXKsqt6xTgUtnLyP5wERWSgimSIyC3gJiAWeEpEEEblHRGaJSLqInCUifwXSgBfs9bNF5GciMl1EMkTkEqzrQRuBjzzqWSYi93m8/x8RuVBEskRkCvAEVuJ51Kltb83iGcM4WlnPsp2H3Q5FKdUHPbumiOGDYjh7RKKj9Tp5xJMG/APrYv/LQB0w2xhTBDRidRZ4BdgDvAYMAhYYY7bY69cDi4C37TJ+D7wDnGefPmuRjXWvTov+WNdrdtrLD7XLXReAbeyQhaMGk5wQqU8nVUo5bldpBR8XlvPlWemEhARuXDZfnOxcsLiNedXAFe2sfwBY6Ec9GV7vvw18278onRUWGsJV04fxx+V7OXSixtFzrEqpvu3ZNUUBfdhbW3SsNpddNX0YzQZeXH/Q7VCUUn1EZV0jr2woDujD3tqiicdlwwbGcPbIRP6xbj+N2slAKeWAloe9Xedwp4IWmniCwFfmZFBaUcs7O8raX1gppbrA82FvUwL4sLe2aOIJAueOSSJtQDRPrSp0OxSlVC/n1MPe2qKJJwiEhgjXzR7O2oLj7CqtcDscpVQv9oxDD3triyaeIHHV9GFEhoXw1Koit0NRSvVSxyrreMuhh721RRNPkBgQG8FlU1L598ZiTlY3uB2OUqoXev7jA4497K0tmniCyFfmZFDT0MSLn+j4bUqp7lXf2MzTqws5e2SiIw97a4smniAyYWg/pg8fwDNrimhuNm6Ho5TqRd7YeoiyijpunJ/pdiiaeILNV+ZmUHSsmg92H2l/YaWU8oMxhsdXFDAiKY6FIwe7HY4mnmBz0fghJMVH8tTqQrdDUUr1EmsLjrP9UAU3zc90fFw2XzTxBJmIsBC+NCud3LwjFBytcjscpVQv8PiKAgbGRnDF1KFuhwJo4glKX5qZTliI8Mxq7VqtlOqa0qpmlu0q49pZ6USFh7odDqCJJyglJURxycQUXvzkAFV1jW6Ho5Tqwd4taiA8JIRr57gzLpsvmniC1PVzh3OqtpFXNha7HYpSqoc6UV3PiuJGPjcllaT4KLfD+ZQmniA1LX0A41MTeHp1IcZo12qlVMc9t24/9U1wUxB0ofakiSdIiQjXz81gd1kla/KPux2OUqqHqW9s5qlVhYwfFMLYlAS3wzmNJp4g9rnJqfSPCddRq5VSHfbm1hLKKuq4ICPc7VDOoIkniEWFh3L1jGG8s6OU4hM1boejlOohjDE8vjKf7MGxTEwMjp5snjTxBLlrZ1k9UZ5bq12rlVL+WVdwnG3FFdw0P4sQl5650xZNPEFu2MAYFo1N5h/rDlDb0OR2OEqpHuDxlQUMiAnn89OC44ZRb5p4eoDr52RwvKqeN7aUuB2KUirIFRyt4r2dZVw7e3jQ3DDqTRNPDzBvxCCyB8fy9OpCt0NRSgW5v31UQHhICNcF0Q2j3jTx9AAtXas3HzzJxv3lboejlApSJ6sbeHH9QT47ObhuGPWmiaeH+Py0NOIiw3hax29TSrXiuXX7qWloCrobRr1p4ukh4iLD+MJZaby+5RCHT9W6HY5SKsg0NFk3jM4bMYhxqcF1w6g3TTw9yPVzM2hsNjy9So96lFKne3NrCaUVtUF/tAOaeHqUzMRYLhw3hKdXF+qo1UqpTxlj+MuKfLIGx5IzKsntcNqliaeHuW1hFhW1jTz/8QG3Q1FKBYmWG0ZvnBccTxhtjyaeHmZq+gBmZg7kiRX5NDQ1ux2OUioIPLGygP4x4Vw5Lc3tUPyiiacHun1hFodO1vL6lkNuh6KUctnew6d4d2cZ184aTnREcN4w6k0TTw+UMyqJUclx/PmDfH1Wj1J93P8t30dUWChfnZfhdih+08TTA4WECLcuyGZX6Sk+2H3E7XCUUi7Zf6yaVzcf4pqZ6QyKi3Q7HL9p4umhPjc5lSEJUfz5g3y3Q1FKueTRD/cRKsKtC7LcDqVDNPH0UBFhIdw0P5PV+cfYfOCE2+EopRxWerKWl9Yf5AvT0xjSL3iHx/FFE08PtnjmMOKjwnjsQz3qUaqv+cuKfJqM4Y6F2W6H0mGaeHqw+KhwvjxrOG9tK6HoWJXb4SilHHKsso6/ry3issmpDBsY43Y4HaaJp4f76rwMwkJCeHxFgduhKKUc8tePCqhrbOZr5/S8ox1wMPGIyBIRMV6vUo/5vxSRXSJSJSLlIrJMROZ6lZHro4zn/aj7ShHZISJ19r9XBGIb3ZCcEMUVU4fywvoDHKusczscpVSAnaxp4OlVRVw8YQgjkuLdDqdTnD7iyQNSPF4TvebdaU+bDxQAS0Uk2auMv3mVcVtbFYrIHOCfwN+BKfa/L4rIrK5uTLC4ZUEWdY3NPKWPTFCq13tmdSGn6hr5Ws4It0PpNKcTT6MxptTj9elNKMaYZ40xy4wx+caY7cB3gHisZOGp2quMk+3UeRew3BhzrzFmpzHmXiDXnt4rjEiK4/xxyTy9upDqeh08VKneqqqukSdWFnDO6MFMGNrP7XA6zenEkyUixSJSICLPi4jPzuciEgHcClQAm7xmLxaRoyKyXUQeEJH2jjXnAO94TXsbmOtj2R7r9oVZnKhu4AUdPFSpXusf6/ZTXt3A18/tuUc7AGEO1rUWuAHYBSQBdwOrRGS8MeYYgIhcCjwPxAAlwPnGmDKPMp4DioBDwHjgPmAycH4b9Q4ByrymldnTfRKRW7ESH8nJyeTm5vq1ga2prKzschn+GNk/hEfe3cmwukJCe8AItS0822d05Wh3gwlCkU2RPtslusnqzdTX26y19vHFib/DQKlvMvzhwxrGDgzhVMEWcv3sT+TU909HOJZ4jDFveb4XkTVAPnA98KA9eTnWqbVE4BbgBRGZY4wpsct4zKOIrSKSD6wVkWnGmA1tVe/1XnxM84z1MeAxgOnTp5ucnJx2tq5tubm5dLUMfzQklXHL0+upHDiKy6YMDXh93cWzfZbkLnE1lmA0unI0eXF5Z0yfE1oN4HNeX9Ja+/hyTc41AY4mcJ5dU8SJum388bqZzBuR6Pd6Tn3/dIRr3amNMZXAdmCkx7QqY8xeY8waY8xNQANwcxvFrAeaPMvwoZQzj26SOPMoqMdbNCaJ7MGxOnioUr1MQ1Mzf8rdx5Rh/ZmbPcjtcLrMtcQjIlHAGKxTaq0JAdoa+W4iENpOGas581Tc+cAqP8LsUUJChNsWZLOjpIKVe4+6HY5Sqpu8uukQxSdq+Po5IxDpOafRW+PkfTwPiMhCEcm0uzK/BMQCT4lIgojcIyKzRCRdRM4Skb8CacAL9vrZIvIzEZkuIhkicgnW9aCNwEce9SwTkfs8qn4YOFdEfiwiY0Tkx8A5wEPObLmzLpuaSlJ8pA4eqlQv0dRs+L/cvYxNSWDR2OB/rLU/nDziSQP+gXW/zstAHTDbGFMENGJ1FngF2AO8BgwCFhhjttjr1wOLsHqk5QG/x+qtdp4xpsmjnmys+3sAMMasAhZjXUvaAnwFuNoYszYwm+muyLBQbpyfycq9R9lW3F5Pc6VUsFu6rZT8I1XceU52rzjaAWc7FyxuY1410OZoAsaYA8BCP+rJ8DHtJawjrD7hS7PS+cP7e/lT7j7++OVpboejlOokYwx/WL6XrMGxXDwhpf0Veggdq60XSogK5/q5w3lzWwl5pafcDkcp1Unv7zrMzpIKvpYzokfdItEeTTy91C1nZxEbEcbDy3a7HYpSqhNajnbSBkRz2ZRUt8PpVpp4eqn+MRF8dV4Gb24tZWdJhdvhKKU6aPW+Y2zcf4LbFmYTHtq7vqp719ao09w8P4v4yDAeek+PepTqSYwxPPBOHskJkXzxrDS3w+l2mnh6sX4x4dw4P5O3t5ex/ZD2cFOqp3h/12E27D/BtxaNIio81O1wup0mnl7uxvmZxEeF8dB7e9wORSnlh+Zmw2/eziNjUAxfnN77jnZAE0+v1y86nJvnZ/HujjK2HtSjHqWC3WtbDrGr9BTfPn9Ur7u206J3bpU6zVfnZ9AvOlyv9SgV5Bqamnnw3d2MGRLPZyf1rp5snjTx9AEJUeHccnYmy3YdZvOBE26Ho5RqxYvrD1J0rJrvXziakF503443TTx9xPVzM+gfo0c9SgWr2oYmHl62m2np/Tl3TO8Yk601mnj6iPiocG5dkMXyvCNs3F/udjhKKS/PrC6irKKOH1w0pteMydYaTTx9yPVzMhgYG8HvtIebUkHlVG0D/5e7l7NHJjI7q+c/b6c9mnj6kNjIMG5dkMWHu4/wSZEe9SgVLB5fUUB5dQM/uHCM26E4QhNPH/OVOcMZFBuh13qUChLHKut4fEU+F08YwsS0fm6H4whNPH1MTEQYty/MZsWeo3xceNztcJTq8/6Uu4+ahia+e8Eot0NxjF+JR0S2isiWdl6bAx2s6h7Xzh5OYlwkv3tXj3qUctOhEzU8vaaIK6elMSIp3u1wHOPvg+DaeojaYOBGILLr4SgnREeEcvvCLO55Yydr848xqw9czFQqGD3y/h4w8K3zRrodiqP8SjzGmJ97TxORaOC7wHVYj6v+YfeGpgLp2tnD+fOH+fzuvd08f+sct8NRqs8pOFrFC+sPct3s4aQNiHE7HEd1+BqPiISIyK3AXuAm4BvAFGPM0u4OTgVOVHgoX8vJZk3+cVbtO+p2OEr1OQ++u5vIsBDuPGeE26E4rkOJR0QuB3YA9wG/A8YYY542xphABKcC65qZ6SQnRPLQu3vQj1Ap5+w4VMFrmw9x47xMBsf3vasU/nYumCciK4G/A/8Gso0xDxhj6gIanQqoqPBQ7jxnBOsKj5O7+4jb4SjVZzzwTh79oou6/50AAB+RSURBVMO5ZUGW26G4wt/OBSuAGuAx4DBwo68hHYwxD3ZfaMoJi2ek89eVBdz35k7OHpFIWC8dhl2pYLG+8Djv7zrMDy8aQ7/ocLfDcYW/iWc/YIDL21jGAJp4epiIsBB+dPEYbn92Ay+sP8iXZqW7HZJSvZYxhv99O4/B8ZHcMDfD7XBc42+vtowAx6FcdOH4IczIGMCD7+bxuSmpxEX6+3tEKdURb28vZV3Bce65fALREb3vkdb+0vMqChHhp58Zx9HKev78wT63w1GqV6ptaOKeN3YyZkg8i2cMczscV/nbueBiESkUkTMGEhKRfva8C7o/POWUKcP687nJqfxlRT4lJ2vcDkepXueJlQUcLK/hZ5eO6/PXUv3d+q8DvzHGnPSeYU+7H/hWdwamnPeDi0bTbOA3b+e5HYpSvUrpyVr+uHwvF40fwtwRiW6H4zp/E88k4L025r8PTO56OMpNaQNiuHFeJi9vKGZb8Rm/MZRSnfS/S3fR2Gz4ySVj3Q4lKPibeAYDzW3MN4AO+NULfO2cbAbGRnDPGzv0plKlusGG/eW8vLGYW87OJH1Q3xoapzX+Jp6DWEc9rZkEFHc9HOW2hKhwvn3eSNbkH2fZzsNuh6NUj9bcbPj5aztIio/kazl9b2ic1vibeN4AfmkPDHoaEYkBfmEvo3qBxTPTyRocy6/e2klDU1sHukqptry8sZjNB07wo4vHEKu3KXzK38RzL9AP2CMiPxSRy+zXj4Dd9rxfBSpI5azw0BB+cvFY8o9U8Y91+90OR6keqbKukfuX7mLKsP5cPmWo2+EEFb8SjzHmMDAX2IKVYF6xX/fa0+YZY8oCFaRy3qKxSczJGsRD7+2horbB7XCU6nH+uHwvR07V8T+fHUdIyJlDjPVl/t7HMwk4YIy5BEgEZgGzgURjzCXGmMLAhajcYN1UOpby6nr+uHyv2+Eo1aMUHaviiRUFfH7aUKamD3A7nKDj76m2jVgJB2NMObAEKxGVByguFQQmDO3H56em8bePCjlwvNrtcJTqMe59YydhocIPLxrjdihByd/E432cuAA4o6OB6n2+d+EoQkRvKlXKXyv3HOWdHWXcec4IkhOi3A4nKPXtcRtUu1L6RXPL2Vn8Z/MhNh044XY4SgW1xqZmfvH6doYNjOam+ZluhxO0/E08xn55T/ObiCwREeP1KvWY/0sR2SUiVSJSLiLLRGRuK2WJiCy1y/hCO/Xe4KNeIyL6U8RPty3MJjEuknte15tKlWrLc+v2s7uskp9eMo6o8L47+nR7/O1YLsCzItLyxNEo4C8ictqJf2PM59opJw/I8Xjf5DXvTqAA6zTet4GlIjLSR4+573qt255qINsr1toOrN+nxUWG8d0LRvHjl7eydFspF09McTskpYJOeVU9v31nN3OzB3Hh+GS3wwlq/iaep7zeP9vJ+hqNMaW+ZhhjTitTRL4D3ARMAd72mD4da0DSswB/u3Cb1upV/rlq+jCeWlXIPW/sZOHowcRE6M1wSnl66L3dnKpt4GefHYevJzSr//L3QXBf7ab6skSkGKgH1gI/Mcbkey8kIhHArUAFsMljejzwD+A2Y8zhDny40SJSBITa5f0/Y8zGLm1JHxMaItxz+QS+8OhqHnpvjw52qJSHvNJTPLt2P1+eNZwxQxLcDifoOfmzdS1wA7ALSALuBlaJyHhjzDEAEbkUeB6IAUqA871Osz0KLDXGvNmBevOAG4HNQDzW0dJHIjLZGLPH1woicitW4iM5OZnc3NwOVHemysrKLpcRLBamhfH4inzSGg+RntA957A922d05ehuKbM3iWyK9Nku0U3WgJN9vc1aax9fAvF32GwMv1pbS1SoYVbMkaD7Ww/G7x/HEo8x5i3P9yKyBsgHrgcetCcvxzq1lgjcArwgInOMMSUich3Woxemd7De1cBqj3pXYR31fAP4ZivrPAY8BjB9+nSTk5PTkSrPkJubS1fLCBZTZtaz6Lcf8MrBKP51+9xuuSPbs32W5C7pcnm9zejK0eTFndmdfU6odYnV17y+pLX28eWanGu6vf6nVxey98R2fvvFyVx6Vlq3l99Vwfj941p3amNMJbAdGOkxrcoYs9cYs8YYcxPQANxsz14EjAMqRaRRRBrt6f8UkZUdqLcJWO9Zr/Jf/5gI7r50LBv3n+A5HcdN9XHFJ2q4/61dnD0ykc9P0/HY/OVa4rG7M4/BOqXWmhAg0v7/T7EevzDF4wXwPeArHahX7HLaqle14fIpQ5k3YhD3L93F4VPaOVD1TcYY7n5lK80GfnXFRO1Q0AGOJR4ReUBEFopIpojMAl4CYoGnRCRBRO4RkVkiki4iZ4nIX4E04AUAY0yxMWab58su+oBnBwX7/p/7PN7/j4hcKCJZIjIFeAIr8Tzq0Kb3OiLCLy+bQF1jM798fafb4Sjliv9sPsTyvCN878LRDBuoD3jrCCePeNKweqTlAS8DdcBsY0wR0AiMxxrxeg/wGtYTTRcYY7Z0sJ5swPNGk/5Y12t2Au8AQ+1y13V+U1TW4DjuzBnBa5sP8eHuI26Ho5SjjlfV8/PXdjB5WH9umJvhdjg9jpOdCxa3Ma8auKITZZ5xbGuMyfB6/22sm1FVN7s9J4tXNxdz97+38c63F+id2qrPuOf1HVTUNHD/lRMJ1UcedJiO1aY6LTIslHsun8D+49U88r7PnulK9Tof7D7CyxuLuSMnW+/Z6SRNPKpL5mZbvXke+zCfPWWn3A5HqYCqqmvkJy9vJXtwLF8/d4Tb4fRYmnhUl/30krHERobx01e20dysg4iq3uu37+ym+EQNv75yEpFhemq5szTxqC4bFBfJTy4ey7rC47z0yUG3w1EqIDbuL+dvqwq4bvZwZmQMdDucHk0Tj+oWX5yexsyMgfzqrZ0cq6xrfwWlepD6xmZ+9K+tDEmI4gcX9e0hirqDJh7VLUSEe6+YQFVdI/e+qff2qN7l0Q/2kVd2insun0B8VLjb4fR4mnhUtxmZHM+tC7J4eUMxq/YddTscpbrF3sOn+MP7e7l0UgqLxupzdrqDJh7Vrb5x7kjSB8Zw9yvbqG3oyLP6lAo+zc2GH/1rKzGRoSz53Hi3w+k1NPGobhUVHsq9V0wg/2gVv3m7b4+arHq+v68tYn1ROXd/ZhyJcZHtr6D8oolHdbuzRw7m+jnDeWJlASv36Ck31TMVHq3i1/bI01fqyNPdShOPCogfXTyW7MGxfPfFTZyornc7HKU6pL6xmW8+v5Gw0BDuv3KSjjzdzTTxqICIjgjl4cVTOVZZz09f2YYxemOp6jl++24eWw6e5P4rJ5LaP9rtcHodTTwqYCYM7cd3LhjFG1tLeGVjsdvhKOWXlXuO8ucP8rlmZjoXTUhpfwXVYZp4VEDdtiCbmRkD+dmr2zlwvNrtcJRq07HKOr79wiZGJMXxs0vHuR1Or6WJRwVUaIjw26smA/CdFzbRpGO5qSBljOH7L23hZHUDv188legIHYstUDTxqIAbNjCGX1w2no8Ly/nzh/vcDkcpn55aVcj7uw7z40vGMC5VH3cQSJp4lCOumDqUz0xK4cF3drOt+KTb4Sh1mp0lFfzqrV2cOyZJnyjqAE08yhEiwr2XTyAxLpJvPb+Rmnod1UAFh5r6Jr7xj430jw7nN1/QrtNO0MSjHNM/JoLfXjWZfUequO8tHUhUBYdfvL6DfUcqefCqKQzS0QkcoYlHOWreiERump/J06uLWJ532O1wVB+3dFsJ/1i3n1sXZDF/ZKLb4fQZmniU475/4WhGJ8fzg5e26LN7lGsOnajhh//ayqS0fnz3fH3GjpM08SjHRYWH8tDiKZysbuDHL2/VUQ2U45qaDXf9cxONTc38fvFUIsL0q9BJ2trKFWNTEvjBRaN5Z0cZHxxsdDsc1cf8cfle1hUc5xeXTSAjMdbtcPocTTzKNTfOy+TskYk8u7OeTQdOuB2O6iM+LjzOw8v2cNmUVD6vo067QhOPck1IiPD7xVPpHync9sx6Dp+qdTsk1csVn6jhjmc/YdiAaH55+QTtOu0STTzKVQNiI/jm1Egqahq549kNNDXrF4EKjIbGEG5+aj11jc08fv0MEqLC3Q6pz9LEo1yXnhDKb744iU+Kylm7M93tcFQvZAys3JpJXmkFj1wzlRFJcW6H1Kdp4lFB4dJJqdyRk83uA0nk7R/sdjiql9m0N5WisoH85JKx5IxOcjucPk8Tjwoa37tgNEMTT7B2Zzpl5fqLVHWPgpKBbN43lBFDj3DT/Ey3w1Fo4lFBJDREWDA5n7joepZvHEFVrZ6DV11z9GQMK7dmktT/FHPGF2lngiChiUcFlcjwJs6dtofGphCWbxxBY5N+UajOqa4N5/0NI4mKaOCcqXsJDdEblYOFJh4VdPrH1XL2pHyOnoxj9fYMdGAD1VGNTcL7G0dQ3xjKoml7iI7Um5SDiSYeFZSGJ59gcnYx+w4lsmu/XgxW/jMGVm3L5OjJOM6elM/AhBq3Q1JeNPGooDVlxCGGJZWzblc6Jcfi3Q5H9RDbCoaQXzKIqSMPMjxZR8QIRpp4VNASgbMn5ZMQU0vupmwqayLcDkkFuQOH+/HJ7jQyhxxjUlaJ2+GoVmjiUUEtIqyZc6ftodkI728YQUOj7rLKt/JT0XywOZtBCdXMm1iIdmALXvpXrIJev9g6Fk7Op7wyhvc3jNSebuoMVTURvPfJSMLDrF6RYaHNboek2qCJR/UIaYNPMn9iASXHE/hgUzbNOqabsp2qF97+ePSnPdhioxrcDkm1w7HEIyJLRMR4vUo95v9SRHaJSJWIlIvIMhGZ20pZIiJL7TK+4EfdV4rIDhGps/+9oju3TTkjO/UYs8cVcuDIAFZszaRZu1n3ebX1YTyxLYbqunDOP2sPif2q3Q5J+cHpI548IMXjNdFr3p32tPlAAbBURJJ9lPNdoMmfCkVkDvBP4O/AFPvfF0VkVie3QbloTPoRzhp1gIKSQazZPlzv8enD6hpCeefjURyrDWHRtD0kDah0OyTlpzCH62s0xpT6mmGMedbzvYh8B7gJK1m87TF9OvAt4CygzI867wKWG2Putd/fKyLn2NOv6fAWKNdNzCqlvjGUrfmphIc1MX30Qb2Q3MfUN4bw7vpRnKiM5ivjqmHQKbdDUh3gdOLJEpFioB5YC/zEGJPvvZCIRAC3AhXAJo/p8cA/gNuMMYf9HHdpDvCI17S3ga+3toKI3GrXT3JyMrm5uf7U06rKysoul9GbebbP6MrRfq0zKgVia+pZU5hCavMAzk2vD2CE7opsivTZLtFNMYD/bdZb1DfB37bHcPxUKF8eU8OkhDDq/GyDvvh3GIzfP04mnrXADcAuIAm4G1glIuONMccARORS4HkgBigBzjfGeB7VPAosNca82YF6h3DmkVGZPd0nY8xjwGMA06dPNzk5OR2o7ky5ubl0tYzezLN9luQu8Xu90ZPgiGTy7v5ETsaUMS7jcGACdNnoytHkxeWdMX1OqHU9w9e83qqxSVi2YSSlFaEsmJxPaMpx6lppH1+uyel7JzmC8fvHscRjjHnL872IrAHygeuBB+3Jy7FOrSUCtwAviMgcY0yJiFwHTAamd6Z6r/fiY5rqYURg3oQCGppCWbdrOOFhzYxMO+p2WCpAmpqF3E0jKDnWj/kT88lMOe52SKqTXOtObYypBLYDIz2mVRlj9hpj1hhjbgIagJvt2YuAcUCliDSKSMuof/8UkZVtVFXKmUc3Sfh3fUgFuZAQWDh5H6mDTrJqWwaFpQPcDkkFQHMzfLg5i4NH+jNnXCEjhh5zOyTVBa4lHhGJAsZgnVJrTQgQaf//p8AkrCOilhfA94CvtFHGauB8r2nnA6s6GLIKUqEhhnOm7mVw/0r7yynB7ZBUN2o2sHJrFkVlA5k5Zj+j04+4HZLqIifv43lARBaKSKbdlfklIBZ4SkQSROQeEZklIukicpaI/BVIA14AMMYUG2O2eb7sog94dlCw7/+5z6Pqh4FzReTHIjJGRH4MnAM85MR2K2eEhzWz6Kw99I+vYfnGEZQe10FFewNjYPX2DPJLBjFt1AHGZeiJit7AySOeNKweaXnAy0AdMNsYUwQ0AuOBV4A9wGvAIGCBMWZLB+vJxrpHCABjzCpgMda1pC1YR0dXG2PWdmlrVNCJDG/igum7iYuu5531o/S0Ww/X1Cys3JrJnoODmZxdzKQsn3diqB7Iyc4Fi9uYVw10eDQBY8wZ/amNMRk+pr2EdYSlermoiEYunrWT9zeMJHdTNjPGHGC8/kruceoaQu0j1wSmjjjIpGwdabo30bHaVK8TFdHEBTPyGJ5czse70lm3c5iOcNCDnKqO4M01YzlcHsfZk/YxeUSJ3iDcyzh9A6lSjggLNSycso+PdzWwo2gIVbURnD0pn7BQzUDB7MiJWJZtGEmzES6YsZshA3VEgt5IE4/qtUIEZo7ZT1xUHR/npVOzPpxFU/cQGeHXMH/KYUWlA/hwSxYxkfWcd9Ye+sXVuh2SChA91aZ6NREYn1nGwsl7OXoiljfWjuVUtT7JNJgYA9sLklm+KZuB8dV8Zs5OTTq9nCYe1SdkppRzwYw8auvCeWPNOI6ejHE7JIV1Y+janel8nJfO8ORyLpy5i6iIxvZXVD2aJh7VZwwZWMkls3cSGtLM0nVjOHikn9sh9WkNjSG8v3Eku/YnMyGzhJwp+/QaXB+hiUf1Kf3javnM7J0kxNaybMNIdh9IdDukPqm6Npy31o2h+Eg/Zo8r1Edb9DGaeFSfExPVwMUzd5EyqIJV2zNZuTWD+kb9U3BK8dEEXl89joqqKBadtYcxOgROn6O92lSfFB7WzHnT9rB5Xwpb9qVSejyesycVkKxPsQyYxqYQ1uelsWt/Mv3jajhv+m4Gxte4HZZygSYe1WeFhBimjjzE0MQKPtySydK1Y5iYXcKU7EOEhOi1hu509GQsH27OpKI6mnHDS5k26qBez+nDNPGoPi9pQCWXzdvO2p3pbNmXSvGRfiyYlK9dertBc7NYR5X5qcRE1nPhjF2k6GOq+zxNPEphnXqbP7GQtMEnWb09g/+sGseMMQcYPeyIXvTupJOVUazYmsnRk3FkpR5l1tj9RIbrzbtKE49Sp8kYUk5S/0pWbs1kzY4MDh7pz7wJBURH6r0l/jIGdu1PYn1eGmGhhpwpe8kYUu52WCqIaOJRyktMVAPnT9/NzqJk1u9O49WPJjBvQgHDkk66HVrQq6oN56OtmRw61o+hiSeYN6GQmKgGt8NSQUYTj1I+iMC4jDJSBlXw4ZYslm0YRVbKMaaOLCY+ps7t8IJOU7Ow+8BgNu4dSnOzMHtcoZ6mVK3SxKNUGwbE13DpnB1s3pvK9sIhFJQOYFTaUSZnH9Jf8lhD3uw7lMimvalU1UaSPKCCeRMKSYjV5Kxap4lHqXaEhhimjSpmdPphtuxLZffBRPYWD2Ls8MNMzCzpk6NdGwOFpQPZuDeViqpoBiVUMXd8IamJFXqUo9qliUcpP8VGNTBnfBHjM0rZtHco2wqGkHdgMBMySxk3vIzwsGa3Qww4Y+DAkf5s3DOU8lMx9I+r5pype0hPOqEJR/lNE49SHZQQW8eCyflMzCphw56hbNyTxs6iZCZlHWJ0+hFCe+HNp8ZAybEENuwZytGTccTH1LJg0j4yUo4ToglHdZAmHqU6aUB8DYum7eVweSwb9qSxbtdwthcOYXL2ITJTjveKIyBjoPR4PJv3pVJ6PIHYqDrmji9gxNBjOrqD6jRNPEp1UdKAKi6ckUfJsQQ+2Z3Gqu2ZrNuVzrCkE2SlHCM1saJHHQUZA8dPxVBQMpCCkoFU1UYSFdHAzDFFvfaITjlLE49S3UAEUhMrSBm0g7LyOApKBlFYOpCCkkFEhDeSkXycrNTjJA84FbTXQk5WRVJQMoiCkoGcrIpGpJmhiRVMG1VMelJ5rziCU8FBE49S3UjEeuDckIGVzBy7n0NHEygoGUR+ySB2H0wiJrKezJTjZKUcY2BCtetJqKomgoJS68jmWEUsYEgecIpxw8sYPqRcnwaqAkITj1IBEhpiGJZ0kmFJJ2loDOHAkf4UHBrIzqIkthcOISGmhtTECgbE1TAgvob+8dVEBPCooqlZqKiKovxUNOWV0Rwuj6esPB6AQQlVzBi9n4yU48Tq/UkqwDTxKOWA8LBmslKOk5VynLr6UIrKBlBQMoi9xYk0NoV+ulxsVB0D4mvsZFTNgPgaGjt4ScUYqKyJ/DTBlJ+K5kRlNCerojDGeuCdSDMD4mqYOuIgmSnH9YZP5ShNPEo5LDKiiVHDjjJq2FE7SURwojKa8lMxnyaL4qMJnyaJ18QQFz0RkdMz0GdOxgLwyooJp02vqo04LZnFRdcyIK6G9KQT9LePrhJia7WTgHKNJh6lXCQC8TH1xMfUnzYIqedpMTmeRlFj9Rnrhodap+UGeDzF0xhITTz539N3cTXaKUAFHU08SgWh0BBjnXKLr2F0wmDy4vadsUz/OCvh5Ew5c55SwSzE7QCUUkr1LZp4lFJKOUoTj1JKKUdp4lFKKeUoTTxKKaUcpYlHKaWUozTxKKWUcpQmHqWUUo4SY3TYjLaIyBGgqIvFJAJHuyGc3krbp23aPm3T9mmbW+0z3Bgz2NcMTTwOEJH1xpjpbscRrLR92qbt0zZtn7YFY/voqTallFKO0sSjlFLKUZp4nPGY2wEEOW2ftmn7tE3bp21B1z56jUcppZSj9IhHKaWUozTxKKWUcpQmHqWUUo7SxNNBIvI1ESkQkVoR+UREzm5j2SgReVJEtohIg4jktrLcQrusWhHJF5HbA7YBAdbd7SMiOSJifLzGBHRDAqSD7ZMjIq+KSImIVNvtdKOP5frq/tNu+/S2/Qc63EbjRGS5iJR57B+/EpEIr+Uc3Yc08XSAiFwNPAz8CpgKrALeEpH0VlYJBWqBPwBvtFJmJvCmXdZU4D7gERG5snujD7xAtI+H8UCKx2tPd8TspE60z1xgK/AFYALwJ+AxEfmSR5l9ef9pt3089Pj9BzrVRvXAU8AFwGjgLuAm4B6PMp3fh4wx+vLzBawF/uI1bQ9wnx/r/gHI9TH9fmCP17THgdVub2+QtE8OYIBEt7fPzfbxWP4F4F+6//jdPr1m/+nGNnrQc/9wYx/SIx4/2YemZwHveM16B+uXV2fN8VHm28B0EQnvQrmOCmD7tFhvn1JZJiLndEN5jurG9kkAyj3e6/5zOu/2adGj9x/onjYSkRHARcAHHpMd34c08fgvEevUUJnX9DJgSBfKHdJKmWF2nT1FoNqnBLgDuBL4PJAHLBORBV0o0w1dbh8RuRRYxOk3BOr+Y2ulfXrL/gNdaCMRWSUitVhHRyuBn3jMdnwfCgtEob2c9x234mNad5Tpa3pP0K3tY4zJw/qyaLFaRDKA7wEfdrZcF3WqfURkHvAc8E1jzDo/yvQ1vSfo1vbphfsPdK6NrgbigcnAb4AfYl3LaatMX9O7hSYe/x0Fmjjzl0USZ/5a6IjSVspsBI51oVynBap9fFkLLO7mMgOt0+0jIvOxLv7+zBjzJ6/ZfX7/aad9fOmJ+w90oY2MMQfs/+4QkVDgcRH5jTGmERf2IT3V5idjTD3wCXC+16zzsXqDdNZq4DwfZa43xjR0oVxHBbB9fJmCdQqlx+hs+9inhN4Cfm6MecjHIn16//GjfXzpcfsPdOvfWAjWQUeo/d75fcjtXho96YV1uFoP3AyMxerWWIn1wCOwDl2Xea0zDmtHfx5Yb/9/isf8TKAKeMgu82a7jivd3t4gaZ+7gMuBkVhdYu/DOvz/vNvbG+j2weqRVYV1amSIx2uw7j9+t0+v2X862UbXAV8ExgBZwFVAMfC8m/uQ6w3Z017A14BCoA7r18cCj3lPAoVeyxfaO/ppL69lFgIb7DILgNvd3s5gaR/gB8BeoAY4DqwALnF7O51oH/v9GW3jow375P7jT/v0tv2nE210jb1vnMJKUNuxOhZEu7kP6ejUSimlHKXXeJRSSjlKE49SSilHaeJRSinlKE08SimlHKWJRymllKM08SillHKUJh6lgoDHA8t60sCeSnWKJh6lXCAiuSLyh55SrlLdSROPUkopR2niUcphIvIk1hAld9qn1wyQYc+eLCJrRaRaRNaLyDSvdeeKyAf2/GIR+ZOIJLRWrohkiEioiDwhIgUiUiMie0TkByKif//KFbrjKeW8b2GNCPw3IMV+tQxbfx/wI2Aa1pD0fxcRARCRiVhPivwP1nNVPo81qOpf2yk3BGtgyKuwBoH8KdZ4XV8N4DYq1Sodq00pF4hILrDNGPN1+30OsBy4yBjztj1tHtbTIocZYw6KyNNAgzHmJo9ypgAbgWRjzGHvctuo/9fAdGOM93D4SgWcPghOqeCyxeP/h+x/k4CDwFnACBG52mOZlidFZgOHWytURG7HGu5+OBANhANF3RSzUh2iiUep4OL54K2W0xEhHv8+DvzOx3rFrRVoJ6qHsB73vAqoAO4EruhqsEp1hiYepdxRz3+fAOmvDcB4Y8zeDpY7H1hrjPm0m7WIZHewbqW6jXYuUModhcBMu9dZIv79Ld5vr/OoiEwVkREicqmI/Lm1cu2ea7uBaSJysYiMFJH/h9X7TSlXaOJRyh0PYB2d7ACOAOntrWCM2QIswOp6/QGwGasXXFk75f4ZeAF4DvjYXv+33bIVSnWC9mpTSinlKD3iUUop5ShNPEoppRyliUcppZSjNPEopZRylCYepZRSjtLEo5RSylGaeJRSSjlKE49SSilH/X8zRQ+JQdvOCAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Make scan of lnL (for theta, if free)\n", "if not(m.is_fixed('theta')):\n", " plt.figure()\n", " m.draw_mnprofile('theta')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAESCAYAAADuVeJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3xP1x/H8dfJkGVLEDMhBLESidh779amRatWbWq2aM0qalSLqKIo1VapXXuvxN4rYhMxEhmyzu+PhJ+2EkHy/SbffJ6Px33w/d5z7/ncjrxz7jhXaa0RQgghXsfM2AUIIYRIGyQwhBBCJIkEhhBCiCSRwBBCCJEkEhhCCCGSRAJDCCFEkkhgCCGESBKDBYZS6kullP7Xcvc129RXSh1QSoUopR4opdYopYoaqmYhhBD/Z+gRxgXA8aWlVEINlVLOwBpgD+AO1AFsgA0pX6YQQoh/szBwf9Fa60RHFS8pB1gCI7TWMQBKqUnAdqWUvdb6QWIb29vbaycnp3cqNq2IjYklOiqGmOgYoqNiiI2OISYmltjoWGJiYoiNiX3xOTYmhpjoWGJjYnmXZ/zNzc0wtzDH3NIccwtzLOL/NLc0x+Kl7+MWOfMpRFrh5+f3QGvt8Kp1hg6MQkqpW0AkcAgYqbW+mkBbXyAK+EQp9SNgC3QGjrwuLACcnJzw9fVNprJTh5BHT7l6MgD/U9fxPxmA/+nrXDt9g/CnEa9sn8HakoxZ7bBzsMMuiy2ZstmRKXtGMma1I1O2jGTMFve9tZ011nZWWNtZYZPRGms7a6xsMxAbE0t4SAShwWFxfz4JI+ThU4KDQnjyIJjgoBAe33/Co3tPeHz/CcH3nvKqqWbMLczJYp+JrDmzkDVnZrLlzopDPnty5s+BQ357HPLnIGcBezJly5jS/wiFEK+hlApIcJ2h5pJSSjUEMgHngZzAF0AxwE1rHZTANlWB3wB74k6fHQMaaq3vJ9C+O9AdoECBAuUCAhI87lQv/Gk4Z/Zf5JLfVS4du8olv6vc9f//YWfKnhHnUgVwLlmA3M45yZYrK9lyZSFbrixkcchMxqx2ZLDOYNCaY6JjeBwYzKN7j3kSGMzj+8E8CQzm0f0ncZ8D44Ll4Z3HPLj1kJjomH9snz13Vpzij8mpZAGcSxWgYIl8WNtaGfQ4hEjPlFJ+WmvPV64z1uSDSqmMwFXga631t69YnxvYDawGlhMXNmPjV9fSWscmtn9PT0+dlkYYWmuungzAd/MJfP8+zuk954iOivuBmqdwLlw8ClHEoxCFyxTEuXRBcjhmQyll5KrfXkxMDI/uPSHwRhCBNx5w71og187ewP/UdQLO3CAyIgoApRSOhXNRqHRBXMo64+LhjIu7Mzkcsxn5CIQwTYkFhqFPSb2gtX6qlDoDFEmgSW8gVGs99PkXSqkPgBtAJWBvyleZssJDIziy8RgH1/nh9/cJHt59DIBzqQK8168RHnXL4OpV2CRP1Zibm2OfJzv2ebJT3Puf/wnExMRw5+p9/E9d59rp6/ifvs7VEwHsXXXoRZtsubLg4u5M4bLOuHoVxq2SK9lyZTX0YQiRrhgtMJRS1sSdktqRQBNbIOZf3z3/nGavooY+CeXguqPsWXUQ303HeRYeSeYcmfCoWxrPemUoV68M9nmyG7tMozI3NydfEUfyFXGk6vveL74PDQ7j6okALh/z5/Jxfy4f9efo1lMvTm3lKZwLt8rFKFHRFbfKrhQskQ8zszT7n4oQqY4hr2FMBdYC14m7hjEKqAaU0loHxN8BVV5rXTu+fS1gK/AV8Atxp6QmAiWA4lrr0MT6S02npGJiYjiy8Tjr52/Bb/MJoiKjyZEnG1Xe86ZqywqUrFoMc3NzY5eZJkVGRHLpqD9n91/gzP7znNl3gceBwUDcdZ7S1UtQprobZWq44VQyvwSIEK+RKq5hKKVWEBcQ9kAgcBAYpbU+G79+EVBDa+300jbtgKFAUSA8fpthz7dJTGoIjEf3n7BpwXbW+2zhXkAg2R2zUat9Faq29KaYdxH54ZUCtNbcvnKX03vPc2r3OU7sOvPiZoHMOTJRunoJ3GuVwr12SfIVzZOmrwMJkRJSRWAYmjED49blO6z8Zg1bft5FVGQ0ZWuVpFmv+lRs5omFpdHOAqZb9wICObnrLMd3nub49tPcvx53V7Z93uyUrVUSjzpxpwPlGogQEhgG43/6Oiu+/pOdK/ZhbmlB/S41aNGvEQWL5zNoHSJhz0cgx7ad5viOUxzffponD0IAKOLhjFcDd7waulPcuwjmFnKaUKQ/Ehgp7M7Veyz+8le2L9uLtZ0VTXvW4/2BTeTWzzQgNjaWK8evcXjjMY5sOsa5AxeJjdVkyp6R8o3cqdjUC8/6ZbDLbGvsUoUwCAmMFBIcFMLPX65kvc8WzMzNeK9fI9oMbU7m7JlStF+RckIePeXolpMcXO/H4Q3HCA4KwcLSnNI13KjSojyVWpSXXwSESZPASGaxsbFs+mkHC0Ys4+njUBp2rc0Ho1pinzdHivQnjCMmOoazBy5ycK0v+/86ws2Ld1BKUbxCESq3KE/VVhVwdM5l7DKFSFYSGMno5qU7TOkym7MHLlKqanH6zu6Kc6mCyd6PSF201lw/d5O9qw6zb/UhLh31B6CoZ2Gqt65ocuHRpUsXHjx4wLp164xdygupsSZTlFhgyH2dSaS1ZtPCHfTyGMKN87cYuqgP03Z+JWGRTiilKFgiPx2/aMkPvt/w85XZdJv8AUrB/GFL6VS4D328h7NyyhruXnvlVGcmr1atWnTq1MnYZfxHaq3rdXbv3k2zZs3ImzcvSikWLVqUYNuePXsycOBAJk2ahJeXF5kzZ8bBwYGmTZty+vTpZKtJAiMJnj4OZUL76Uzr+gOuXi7MOzGNup2qyz386Zijcy7aDGnO7ENf8/OV2Xzy9QdoHRceHxbqTd8KI1jz/SaCH4YYu1SDOXbsGOXKlTN2Gf+RWut6nadPn1KyZElmzpyJjY1Ngu201qxdu5bmzZuzc+dOPv30U/bv38/27duxsLCgTp06PHz4MFlqkocCXuP03nNM+mAWQbcf8fGEDrQZ2kyeyhb/4Oici7ZDm9N2aHPu+N9j928H2b58D7P7LmDe4MVUaOZJvU418KxfJk0+h1OjRg1KlChB1qxZ8fHxwczMjE6dOvHNN9+8ePj0ypUrPH782KA/mFNrXcmlUaNGNGrUCIg7HZeQI0eOEBERQZUqVdi8efM/1i1ZsoQsWbKwb98+mjZt+s41yQgjEdt/2cNntb7CwtKcGXvH0X7EexIWIlHPw2PesanMOfoNTXvV5+TOM4xq9jXt8/dk7qBFXD7ub+wy39iyZcuwsLBg//79zJ49mxkzZvDrr7++WO/n54eZmRlly5Z97b4mTpxIxowZE1327Nlj8LrSqtWrV9O4cWMsLP77y0hISAixsbFky5Y8d/alvV93DOSP6euYO3gxpauX4Ks/h5Ixq52xSxJpjEtZZ1zKOvPJ5I4c2XScLT/vYs33m/hjxnpc3J1p8HEtanWokiZmIy5RogRjx8a9XaBo0aLMnz+fbdu20b59eyDuB3PRokXJmPH1x9KzZ0/atGmTaJu8efMavK60as2aNYwbN+6V6/r370/ZsmWpWLFisvQlgfEvWmt+HL6MlVPWULWlN8OX9DP4i4iEabHMYEmlZl5UauZFcFAI25fvZdNP25nddwE+Q36massKNOpWh1JVi6fa62KlS5f+x+c8efJw//7/L+77+fkl+bRP9uzZyZ49eWZkTs66UsIXX3zBhAkTEm2zY8cOatSo8Vb7v3z5MlevXqV+/fr/WTdo0CD27t3L3r17k+3MiJySeklsbCzfdpvLyilraNqzHp+vGChhIZJV5hyZaNGnIXOPTuEH38nU71KTA2t9GVxjDF3dBvLH9HWp8kK5paXlPz4rpYiN/f87zN7kwnJynpJ6l7quXbtGmTJl6Ny5M0WLFqVXr16sXr0ab29v3NzcuHTpEgDNmzenXLlyuLm5sWzZMgD279+Pp6cnUVFRBAUF4erqyu3bt//Tx4ABAzh37lyiS/ny5ZN0rK+yevVqateujZ3dP8+ADBw4kOXLl7N9+3YKFSr01vv/NxlhxNNaM+vTH9n003Y6ft6SzmPbptrf9oRpKOJRiCI/FKLblA/ZtfIAG+ZvYe7gxfz0+S/UbFeFpp/Wx9WzsLHLfC1/f38ePnyIh4dHkton5ympd63r3LlzrFy5EhcXF0qWLEnGjBk5dOgQs2fPZvbs2cycOZOFCxeSPXt2QkJC8Pb2plWrVlSqVIkaNWowceJEzp07x/Dhw8mTJ89/9m9vb4+9vf07H0tC1qxZQ+fOnf/xXf/+/VmxYgU7d+6kWLFiydqfBAZxYTFn4CLW+2yh3bAWEhbCoGzsrGnwUU0afFSTKyeusXbO32xbtpvNi3ZQ1LMwTXvVp0bbSqn23eZ+fn5A3G/7L9/zb25uTvHixf/TPjlPSb1rXa6urri6ugJQvHhx6tSpA8Sd6tq6dSsAM2bMYM2aNQBcv36d69evU6RIEcaPH0/ZsmVxcnLio48+Svb6nz59yuXLl4G4sx/Xr1/n+PHjZM+enQIFChAYGMjBgwf5/fffX2zTu3dvlixZwurVq8mWLRt3794FeDFye1fp/pSU1pqfRv7Cn7M28F6/Rnw8sYOEhTCawmWcGDC3OytuzqP3rI+JCI1gWtcf6JC/B3MHL+b2lbvGLvE/nv9grly5MqVKlXqxvG4UkRrqsrL6fwibmZm9+GxmZkZ0dDQ7duxgz549HDx4kBMnTlCsWDGePXsGwP3794mMjCQoKIjo6Ohkr9/X1xd3d3fc3d0JDw9nzJgxuLu7M3r0aADWrl2Ll5cXuXL9f4aBH374gZCQEGrXro2jo+OLZerUqclTlNbaJJdy5crppNi0cLuuo1rp6T3m6djY2CRtI4ShxMbG6uM7TuuxbabpehZtdF2z1npk4wn60IajOiYmxtjlpWn+/v765Z8TLVu21Dt27NBaa71nzx7duHFjvXr1at2qVSuttdbHjh3TFhYW+tSpU1prrevVq6dXr16tBwwYoMeNG2fw+ps1a6YnT56c7PsFfHUCP1fT9SmpgHM3md1nAWVquNH3+64yshCpjlKKMjXiXjH74FYQ6322st5nC583nkgel9w061Wfel1qpIlbc9OiBg0aMGfOHMqWLUuJEiVeXEBfsGAB2bNnp3nz5tSrVw8vLy+aN29OqVKlDFZb5cqVX9w+bCjpdvLBZ+HP6FthJA/vPGLu8anY50n5c6pCJIeoyCj2/HGIv37YxJl9F7C2s6JJj3q0HNRE/jsW70wmH3wFnyFL8D91naGL+8r/ZCJNscxgSa32VZixZzxz/L6hcovyrJqxjk6FejOjpw93rt4zdonCRKXLEcbpvecYWG007/VrxKczkv/uBiEM7faVu6yc8hd/L9pBTEwsNdtVpt3w93Byy2/s0kQaI+/DeElUZBS9PIYS/jSCH09/i03GhGeBFCKteXD7IX98u4518/4mIvQZlVt40X7E+7h6uRi7NJFGyCmpl6z85i8Czt6k3/efSFgIk2OfJzs9pnZi2bU5fDi6NSd3naWP9wiG1R/HiZ1nMNVfEIVhpKsRRsDZG/TyGErF5l6M+nWQkSoTwnDCQsJZN/dvfv92LY/uPaFEJVc+GNUKz3pl5K5A8UqpYoShlPpSKaX/tST6FJKKM0ApdV4p9UwpdUcp9fXb9B8TE8PUrnOwyWRDn1kfv91BCJHG2Gayoc2Q5iy5+j19Z3/Cg5tBjGw4gX6VPufwxmMy4hBvxNCnpC4Aji8tr7tpeRrwKTAMKA40Ana/Tcerpq/n/KFL9J71MdlyZX2bXQiRZlnZWNHs0/osujiLAXO78+juYz5vPJF+FUdyaMNRCQ6RJAY7JaWU+hJopbUumcT2rsBpoLTW+tyb9vfyKan71wPp4tofrwZl+XLVEBmKi3QvKjKKLT/vZvnEP7h7LZCinoX5YFQrKjQpJ/9/pHOp4pRUvEJKqVtKKX+l1AqlVGLz7jYHrgINlFJXlVLXlFKLlVI5E9pAKdVdKeWrlPINDAx88f3+v3yJehZFt8kfyP8MQhD3LEejT2qz8MIsBs3vScjDp4xuPplPPYex989D/5giXIjnDBkYh4AuQEOgG5Ab2K+UypFA+0JAQaBd/HYfAsWAtUqpV9attfbRWntqrT0dHBxefH9s2ylyO+ckX9H/Tj8sRHpmYWlBw661+encDIYs7E1YSDhftZxKL4+h7PnjoJyqEv9gsMDQWm/UWq/UWp/UWm8FmsT33zmBTcwAK+BDrfVurfUe4kKjPOCV1H5jomM4vuM0HrUNN8eLEGmNhaUF9TrX4KezMxj2c1+inkUxtvU0BlQdxblDl4xdnkgljPYchtb6KXAGKJJAkztAtNb64kvfXQKigQJJ7efUnnOEBYdTrl6Zt65ViPTC3MKcOh9UY/7pbxno05M7V+7Sr+JIJn0wk/vXA1+/A2HSjBYYSilr4k4x3UmgyT7AQin18ivHChH30qeApPaz89f9WNtZ4dXQ/a1rFSK9MTc3p9EntVl08Ts6jHyfvasO8VGx/iwY+QuhwWHGLk8YiSGfw5iqlKqulHJWSnkDvwN2wOL49ZOUUtte2mQrcBT4SSnlrpRyB34i7lpIwtPQviQ6Kpo9fxykYjNPbOysk/V4hEgPbDPZ8NH49iw8P5Mq73uz4us/6ezSh9WzNxIdlfwvDRKpmyFHGPmA5cQ9i7EKeAZU0Fo/Hy04Ai9GE1rrWOKuc9wn7tmLzcBNoHn8utc6vuMMwUEh1GhbOdkOQoj0KGcBB0Ys7c/sw1/jVLIA3/f7iU9KDmLPqkNyYTwdMempQXrVHchv09ay5slirGxS5/uQhUhrtNYcWn+U+cOWcP3cLUpVLU6v6V0o4pHYXfIirUhNz2EY1NmDF3Fxd5KwECIZKaWo0KQcPiem0X9Od26cv0Vvr+FM6/oDD+8+MnZ5IgWZbmBouHjkCsUrFDV2JUKYJHMLc5r0qMuii7NoObAJW5fu5iPX/qyYvJrIiEhjlydSgMkGRnhoBBFhzyhZuZixSxHCpNllsaPH1E7MPz2d0jVKsGDEMrqWGMCOFfvk+oaJMdlrGAVzO2u3xxX47f4C7DLbGrscIdKNo1tPMu+zn7l6MoDiFYrQY2pn3Cq5GrsskUTp8hrG00ehuNcpJWEhhIF51CnND36TGbzgU+4FPGBAlS8Y3+5bHtx+aOzSxDsy2cCIioymynvexi5DiHTJ3NycBh/VZNHFWXQa04YDf/nyidtANszfKqep0jCTDQwA+7zZjV2CEOmajZ01H45pjc/Jabi4OzO9xzyG1P6KW5cTmuBBpGYmHRihT2QKAyFSg7wujkzZNoaB83pw6ehVupcezC8TVxEVGWXs0sQbMOnAePpYAkOI1EIpRaNudVhwdgbejT1Y+MVyenkM5dSeN34/mjASkw4MuRdciNTHPk92Rv/2GePXDudZ2DMGVR/N1I9/IDgoxNilidcw6cCwyyJ3SAmRWnk3Lsf809NpN6wFW5fupmuJAez8VZ7dSM1MOjAyZcto7BKEEImwtrWi66SOzPGbTC4nBya0n8GY977hwa0gY5cmXsGkA0NGGEKkDc6lCjJz/wS6T+nE0S0n6eo2kHXztsi7xVMZkw4MSysLY5cghEgic3NzWg9uyrwTUylarhAze/kwsNpo/E8l+X1pIoWZdGAoM5M+PCFMUl4XR77ZOoahi/pw88JtepUbxo/DlxIR9szYpaV7Jv0TVSljVyCEeBtKKep2qs7C8zOp+2E1fv1mDd1KDeLswYvGLi1dM/HAkMQQIi3LnCMTgxd8yrSdX4HWDKw6iuWT/pRrG0Zi0oEhhDANpauVYO6xKVRrVYGfPv+FYfXGyWSGRmDSgXHzosxXI4SpsMtix8hfBjD4x16cP3iJ7qUHs2vlfmOXla6YbGCYmSnO7L9g7DKEEMlIKUWDj2sx5+g35C3iyPh205nQYQbBD+UpcUMw2cCwtrPikt8VY5chhEgB+YrmYcaecXQZ1449vx+kW6nBHNl0zNhlmTyTDQxLK0vuX39g7DKEECnE3MKcjp+35LuDE8mUzY6RjSYyo6cP4U/DjV2ayTJYYCilvlRK6X8td5O4bRGlVIhS6mlS+7PIYMHDu4+JfCbTJwthyop4FOIH38m0GtSUDfO30qPsEE7vlRlwU4KhRxgXAMeXllKv20AplQFYAex+k44yWFsCcO309TcuUgiRtmSwzkCPqZ2YuuNLtNYMqj6GpeN+l4kMk5mhAyNaa333pSUwCdtMBk4Cv71JR9Z21gCcP3T5zasUQqRJpauVYN7xqdTqWIXFY37l2LZTxi7JpBg6MAoppW4ppfyVUiuUUoUSa6yUagw0Afq9aUeWGSzI6pCZS0evvm2tQog0yDaTDcMW92XKtjF41Cn92van9pzj4Do/A1SW9hkyMA4BXYCGQDcgN7BfKZXjVY2VUo7AfOBDrXWS7plTSnVXSvkqpXwDAwPJ5eQg0yQLkQ4ppShbs+Rr20U+i+LO1Xt8220OK6esMUBlaZvBAkNrvVFrvVJrfVJrvZW4kYMZ0DmBTZYCc7TWB9+gDx+ttafW2tPBwYEcebLz4JY8DSqEeLUMVpbU6lCF/MXyEhYid1e9jtFuq9VaPwXOAEUSaFILGKOUilZKRQMLALv4z92T0kfhMk4EnLnJ/Rtye60Q4r9iYmKY3WcB2XJlocvYdgByoTwRRgsMpZQ1UAxIaP6OUkDZl5bRQHj835N0Abxup+pordny8653L1gIYXJ+m7qWgHM3GbKwNxAXIC9PWirh8U+GfA5jqlKqulLKWSnlDfwO2AGL49dPUkpte95ea3365QW4BcTGf36UlD4dC+WiTA03ti3bkwJHJIRIi25cuAXA3j8Psf2XPQz06YmVjRUxMTGYm5sbubrUzZCvpMsHLAfsgUDgIFBBa/38dVqOQOHk7rRiU0/mDl7M/RsPyJnfPrl3L4RIQyLCnjGz13wyZrMj8EYQvWd+TIFieYmNjcXc3JzY2FjMzMwIuvOIgLM3ObjWlyz2men4RUtjl54qGPKidzutdR6tdQatdV6tdUut9dmX1nfRWjslsv0irXXGN+3XvXbcs4FyP7YQwtrWinFrh/Ms7Bl3/e9Tpobbi3UxMTEvwmLRF8vZ8/sB8hZx5OzBC3zz0WwjVp16mOxcUs85lcxP5hyZOL33vLFLEUKkAjZ21kza+AV1O1Xn78U7CXn0FDMzsxcjjAUjl5EjT3ZaDW5K894NGLqoD9a21nI9g3QQGGZmZriWd+H8oUvGLkUIkYr0nNYZV6/CbFm8i3sBcZNOrJu7BWsbK+p1qUFeF0cAfp28mvvXA+UNnhj2GobRFPcugu+m49zxv4ejcy5jlyOESCUKlsiPtZ31i7nnHtwKwrW8C9kdswGwb/VhLvpdZaBPDwA2zN/KkwchWFia0/qzZkar21hMfoQBUP+jmljbWfFd7x9lWCmE+IdcBR3Ilisrkc+i8D99naKehbG2teLS0ausnr2R5n0aojWsnLKGP7/bQIHiednx6z5+mbjK2KUbXLoIjJz57flofHuObDrOjhX7jF2OECIVymBliaNzLr7r/SPbf9nD+LbfUrNtZYp4OOO7+Tg3zt/i8+UDqdyiPJ98/QE3LtwiJibG2GUbVLoIDIBmvetTqExBlk9aJaMMIcQrfTrjIzzrlyXkUSidvmxLo251uH35LucPX6JOp+o4ueUnNDiMY9tOUai0U7p7biNdXMMAMDc3p0mPesz6dD6Xjl6laLlkf+RDCGECOox8/x+fl4z9jYpNPSlTPe4W3AuHLxP6JIyKTcsRHhrB8e2nuXDkMgVL5Kdmu8rGKNlg0s0IA6BG20pYZrBgx3I5LSWEeL2nj0PJmjMLbYY0B+DErjMc2XQcR+eclKjoyoR20zm41he7zLb8+s1qfv92rZErTlnpZoQBkClbRmwz2xARGmHsUoQQaYBNJmuio6IZ22YaziULcP3cTZzcCtD6s2Z889FstNYM9OkJgEP+HAScvWnkilNWugoMgLCQCGwz2Ri7DCFEGmBubs64NcNZNv4PIiMiaT/ifQqVLshfP2zm4pEr/Hh6+ou2QbcfERMdg9baZJ/ZSFeBobUmU/aM7Fl1iPcGNMY+T3ZjlySESAM6ftHyRRAE3XnEkU3HGPZz3xfr/U8FsHLqX3y5aojJhgWks2sYSinG/PEZj+8/YXi9cTx5EGzskoQQacTzILDLYotdFlvyFol7Ejw0OIzRzSfTdmhzinsXMem7MNNVYACUqFCUcX8N587Ve3zR9Ot0dx+1EOLdxETHcP/GAxaMWMamhTuY0H4GFZp48n7/xgAywjA1ZWq4MdCnJ+cPXWLbUnlXhhAi6ewy2zL571FERkTxJDCYuh9Wp/esjwHTf+GSMtUD9PT01L6+vgmuj42NpV/FkTy8+5iF52diZWNlwOqEEKbGVC52K6X8tNaer1qXLkcYEDeL7ccTOhB4I4j9axIOFiGESApTCIvXSbeBAVCmphu2mWw4tfvs6xsLIcQbMMXro+k6MMzNzXGr7MrRbadM8l+uEMI4/pi+jgFVRnHz0h1jl5Ks0nVgANTtVINbl+6w8pu/jF2KEMJE2OfLwc0Lt+nlPoR187aYzMXwdB8YNdpWonqbiiwe8ysX/a4YuxwhhAmo3roi809No0RlV2b28uGLppMIDgoxdlnvLN0HhlKKfj90I2vOzHzXZ4HJ/CYghDAu+7w5mLTxc3rP/JhjW08xvcc8Y5f0ztJ9YABkzp6J9iPe5/yhS5zZf8HY5QghTISZmRkt+jbkwzFt2LvqEAfX+Rm7pHcigRGvXpcaZMqekWXjfyc2NtbY5QghTEirwU1wcsvPd31+TNNTEhksMJRSXyql9L+Wu4m0r6GUWqOUuqOUClNKnVRKfZxS9dnYWfPBqFb4bj7BbDk1JYRIRpYZLBn0Y6+4eezqjyfk0VNjl/RWDD3CuAA4vrSUSqRtJeAU0AooCcwBfJRSHVKquPf6NaLt0Oasnfs38z77WUJDCJFsinsX4djKX50AACAASURBVMtVQwg4c4ORjSYSGhxm7JLemKEDI1prffelJTChhlrriVrrL7TW+7TWV7XWc4BVQMuUKk4pRddJHWnRpyF/TF/H34t3plRXQoh0yKuBO1/8OohLflcZXm8cj+49NnZJb8TQgVFIKXVLKeWvlFqhlCr0httnBh4ltFIp1V0p5auU8g0MTDCLEqWUoteMLpSo5Mr8oUsIfpj2b4UTQqQelZp7Mfr3wfifuk6/iiMJOJd23tJnyMA4BHQBGgLdgNzAfqVUjqRsrJRqAtQGfBJqo7X20Vp7aq09HRwc3rpQMzMz+v/QjZBHoSz8fPlb70cIIV6lUjMvpu38imfhkfSv9DnHd5w2dklJYrDA0Fpv1Fqv1Fqf1FpvBZrE99/5ddsqpSoDvwD9tNaHU7hUAAqVLkitDlXYtXK/XMsQQiQ7Vy8XZh2YiH3e7Ixq9jWXj/sbu6TXMtpttVrrp8AZoEhi7ZRSVYCNwOj46xgG41apGCGPQrl77b4huxVCpBO5nXIyectoMmXLyOjmk3l4N8Ez7qmC0QJDKWUNFAMSnJ1LKVWNuLD4Sms9w1C1PefqVRiAFZP+JDoq2tDdCyHSgRyO2Ri7ZhghQU8Z894Uwp+GG7ukBCUaGEqpYKWUffzfQ+I/v3J5XUdKqalKqepKKWellDfwO2AHLI5fP0kpte2l9jWIC4u5wDKlVO745e0vTrwhF3dnWg9uyoYftzGiwfg0/cCNECL1cnF3ZvjSflz0vcLnTSal2tBI9I17SqnOwAqt9bP4vydIa7040Y6UWgFUA+yBQOAgMEprfTZ+/SKghtba6aXPr+oz4HmbxLzujXtvYsuSXUzvPo8cjlkZv34kBYvnS5b9CiHEy3b+uo9JHWfiVqUYE9aNwCajjcFrSOyNe0l+RatS6n2t9apXfK+AoVrrye9WZvJKzsAAuHDkMqOafY2VrRWzD00ii33mZNu3EEI89zw0vBq6M+6v4QZ/k19yvaJ1mVLqR6WU7Us7zgfsAAa+Y42pnquXC1+tHkbQ7UeMa/OtXNMQQqSIGm0r02NaZw6tP8r2X/Yau5x/eJPA8AYqAMeVUp5KqbbETd0RDpRJieJSm+LeRRg0vycndp5h6bjfjV2OEMJENe/TgOIVivDDgIU8Dnxi7HJeSHJgaK1PAp7AXuAAsAQYo7VuqLW+l0L1pTp1PqhG1VYVWP3dRsJCUueFKSFE2mZubs6g+b0ICw5jYoeZREVGGbsk4M1vqy0DVAcuA5FAeaVUpmSvKpVrO6Q5oU/CWD9vi7FLEUKYKCe3/Aya34tj204xs+f8VPEAcZIDQyk1CtgNrCEuOMoBrsAppVTVlCkvdXL1csGzfhkWfrGcI5uPG7scIYSJqtupOh+MasXmRTtYNGqF0UPjTUYYvYCmWutBWutIrfUFoCKwAtiaItWlYiOW9adAiXyMafENfltOGLsckQYsW7YMJycnzMzMcHJyYtmyZcYuSaQBnb5sQ4OPavLLxFVM+mAmEWHPjFeM1jpJC2CfyLpqSd2PoZZy5crplPbkQbDuXmawbmTTXp/YdSbF+xNp19KlS7Wtra0GXiy2trZ66dKlxi5NpAGxsbH6l4mrdF2z1rqX51B9/8aDFOsL8NUJ/FxN8nMYaU1yP4eRkMeBTxhUfQwP7zzi211jKVS6YIr3KdIeJycnAgIC/vN9wYIFuXbtmuELEmnSgbW+fP3BLGwyWbPw/MwUebAvuZ7DEK+Q1SELX2/6HJuM1oxoMJ47/unmhjHxBq5fv/5G3wvxKhWbejLmj88Iuv2IfauPGLx/CYxkkLOAA5M2fUFkRBRftZxKZESksUsSqUyBAgXe6HshElK2VklyFXRg27LdBu9bAiOZOLnlZ9jPfbly/BpzByU6rZZIhyZMmICtre0/vrO1tWXChAlGqkikVWZmZtT5oBpHt5zk9N5zhu3boL2ZuApNytHms2asnfs3G+anuxvHRCI6duyIj48PBQsWRClFwYIF8fHxoWPHjsYuTaRBrYc0I7dzTiZ2nEnIo6cG61cueiez6KhoRjX7Gt/NJ+g+pROtBzc1eA1CCNN34chl+lf+ggpNyzH6t8GYmSXP7/9y0duALCwt+Gr1MKq1rojPkJ/56fNfjP6wjRDC9Lh6udBt8gfs+/Mw3/f7ySA/ZyxSvId0KIOVJSN/6Y9dZluWT/oTC0sLOn3ZxthlCSFMzPsDGhN0+yG/TVuLXRZbPp7QIUX7k8BIIebm5gz06UFsTCxLxv5GZvtMtOjT0NhlCSFMiFKKbt98SOiTMJZP+pPcTjlp1K1OivUnp6RSkFKKgT49qNTcix/6L2TLkl3GLkkIYWKUUvSb042yNd34cfjSFH2VtARGCjO3MGfkL/0pU9ONKV2+Z72PzHArhEhe5ubm9J7VldDgcBaNWpFi/UhgGICVjRXj1w7Hq2FZZvT04Y/p64xdkhDCxDi55ad57was99nKBd8rKdKHBIaBWNlY8eWqIVRt6c3cwYtZO2ezsUsSQpiYzl+1IVvurMzsOY+Y6Jhk378EhgFZZrDk8+UD8W7swey+Cziw1vDPiQghTJddFjs+nd6FS0f9WTN7U7LvXwLDwMwtzPl8+QAKuzszsf0MTuw6Y+yShBAmpFrripRv5M6CkcvwP/XfGZLfhcECQyn1pVJK/2u5+5ptSimldimlwpVSt5RSo5VSylA1pxSbjDZMWDcCh/w5GFZ3nFwIF0IkG6UUny34lIxZ7RjXdjrhoRHJtm9DjzAuAI4vLaUSaqiUygxsAe4BXkA/YAgwKOXLTHnZcmVl1oGJeNQpxYyePszq/SPRUdHGLksIYQKy5crKsCX9uHnhNvOHLEm2/Ro6MKK11ndfWgITadsRsAU6a61Pa63/ACYDg0xhlAHE/QawdjitBzdl7ZzNjG7xjXFfvyiEMBketUvRtFc9Nvy4jfs3HiTLPg0dGIXiTy35K6VWKKUKJdK2IrBHax3+0nebgTyAU0oWaUjm5uZ0n9KJgfN64LvpOCMbTSA0OMzYZQkhTECbIc0B+H3a2mTZnyED4xDQBWgIdANyA/uVUjkSaJ+buNNRL7v30rr/UEp1V0r5KqV8AwMTG7ykPo261WHEsv6c3X+RoXXG8uj+E2OXJIRI43IVdKBWhypsmL+VwJtB77w/gwWG1nqj1nql1vqk1nor0CS+/86JbfavzyqB75/34aO19tRaezo4OLx70QZWs11lxvzxGddOX6ev94hkv8NBCJH+dPqyDVpr5g9792sZRrutVmv9FDgDFEmgyV3+O5LIGf+nyb44u2JTT77dNZaoyGj6V/6CQ+v9jF2SECINy+2UkzZDmrNj+T5O7j77TvsyWmAopayBYsCdBJocAKrGt3uuLnAbuJay1RmXq5cL3x+eRL6ijoxqNplNP203dklCiDSs7bAWOOTL8c7v5zHkcxhTlVLVlVLOSilv4HfADlgcv36SUmrbS5v8AoQBi5RSJZVS7wPDgW91OngjkX3eHEzbNRaPuqWZ9skceVZDCPHWrG2taDusBWf2XeDkrrcfZRhyhJEPWE7csxirgGdABa318xP1jkDh54211k+IG1HkAXyB74FpwLcGrNmobOysGbt6KOUbuTOjpw+rZqyXt/cJId5Kg49rkj13VpaM/e2tf44Y8qJ3O611Hq11Bq11Xq11S6312ZfWd9FaO/1rm1Na62paa2uttaPW+qv0MLp4WQbrDIz5YwiVW3gxZ9AixrX91qAvfRdCmAYrGyvaj3yfEzvPsPv3g2+1D5lLKg3IYGXJ6N8/o+ukjuxffYSe7kM4vfecscsSQqQxTXvWw8XdmTkDFxIWEv76Df5FAiONMDMzo92wFszYOw4LS3MG1xjDb1P/klNUQogkM7cwp98P3Qi6/Yg/Z2544+0lMNKYYuWL8IPfN1R+3xufoUuY8tH3REZEGrssIUQaUdy7CKWqFmf37wfeeFsJjDTILrMtX6wYSKcxbdjy8y4+q/UlD24/NHZZQog0osp73lw9GcCtywk91fBqEhhplJmZGR+Oac3o3wbjf/I6vdyHcGTzcWOXJYRIA6q09EYpxd+Ldr7RdhIYaVzVlhWYfXgSWXNlYWTDCcwftlSmSRdCJCpnfnsqNfdk3bwtPAtP+gzZEhgmoGCJ/Mw+NInG3eqwcsoaBtf8UiYvFEIk6r3+jQkOCmHb0j1J3kYCw0RY2VgxYF4PRv4ygCvH/OlXYQQBZ28YuywhRCpVuloJCpbIx9Zlu5O8jQSGianZrjLTdn5FZEQU/Sp9Ltc1hBCvpJSiWquKnN5znqA7j5K0jQSGCXL1cuG7gxPJVdCBkQ0nMKPHPEKfhBq7LCFEKlOtdUW01uxcsS9J7SUwTFTOAg7MOjCRVoOasnHBNrq6DWTf6sPGLksIkYo4ueWnZJVirJq5Pkk3y0hgmDBrWyt6TO3ErAMTyWKfmS/fn8L4djIXlRDi/9oMac796w/Y/dvrH+STwEgHXL1c+P7I13QZ1469qw7To+xnnNojc1EJIcC7sQd5Cudi+/K9r20rgZFOWFha0PHzlszcNx7LDBZ8VnMMi0avkGc2hEjnzMzMcPEoxI0Lt1/f1gD1iFTE1cuFOUenUPuDaiwb/wd9K4zk8nF/Y5clhDCi/EXzcNf/PlGRUYm2k8BIh2wz2TB0UR9G/zaYoNsP6VN+BAu/WE7ks8T/YxFCmKZCZQoSGxPLRd+ribaTwEjHqraswI9nplOrQxV+mbiKXh5DZLQhRDpUtmZJlFIc3XIy0XYSGOlc5uyZGLqoDxPWjyT0SRj9Koxkzfeb5D0bQqQjmXNkoki5QhzbcSrRdhIYAoDyDd2Ze2wK7nVKMbvvAr5qNVVuvxUiHSniUYhrp64n2kYCQ7yQ1SEL4/4aTvcpnTi41o+Pi/Vnvc8WYmJijF2aECKFFSiWl5BHic8IIYEh/sHMzIzWg5vy3cGJ5HPNw4yePnxabhhHtyU+VBVCpG0F3fK9to0EhnilIh6F+HbXWEatHERYcBjD6o5ldIvJ3L1239ilCSFSgFvlYlhaWSbaRgJDJOj5bJYLzs6g68QOHNt2im4lB/Hb1L+IiZbTVEKYEmtbK8rUKJFoG6MFhlJqpFJKK6Vmv6ZdfaXUAaVUiFLqgVJqjVKqqKHqFJDBOgPthr/HgjPTKVu7JD5Dl9C7/HAuHLls7NKEEMnIya1AouuNEhhKqQpANyDRm36VUs7AGmAP4A7UAWyADSldo/ivnAUcGLt6GKN//4zH95/Qt8JIvu02l0f3Hhu7NCFEMgh/GpHoeoMHhlIqC7AM6Aq87q0d5QBLYITW+rLW+jgwCSislLJP2UrFqyilqPq+NwvOTOf9/o34e/FOuhTtx8opa+RJcSHSuLCQsETXG2OE4QP8rrXenoS2vkAU8IlSylwplQnoDBzRWj/4d2OlVHellK9SyjcwMDB5qxb/YJfFjp7fdmH+qWmUqlac+cOW0q3kQA5tOGrs0oQQb8kig0Wi6w0aGEqpboALMCop7bXW14C6wFfAM+AJUApokkB7H621p9ba08HBIVlqFonL75qX8WtHMHHj51hksOCLJpOY0GGGnKYSIg2yy2yb6HqDBYZSyhWYCHTUWkcmcZvcwALgZ8ALqAGEACuVUnKHVyriVb8sc49NofNXbdm36hBdSwxg08IdMsWIEGmIta1VousN+UO3ImAPnFZKRSulooHqwKfxn19VaW8gVGs9VGt9TGu9G/ggfrtKBqtcJIllBks+GNWKucenUtAtP9O6/kC/Sp9zZPNxCQ4h0oCwkPBE1xsyMFYTdzqp7EuLL7Ai/u+vGnXYAv++4f/5ZxlhpFIFiuVl2s6vGDS/Jw/vPGJkwwn0ryzBIURq97r54wz2Q1dr/VhrffrlBQgFHsZ/1kqpSUqpbS9tth7wUEqNUUoVUUp5AAuBG4CfoWoXb87MzIyGXWuz6OIsBsztTtDt/weH35YTxi5PCPEKz8ISv1qQ2n5LdwQKP/8QfydVB6A5cAzYTNxdUw201onPkiVSBcsMljTuXpdFF2cxcF4PHt55zPD64xlWfxyXjib+shYhhGHZZUn8orcy1VMEnp6e2tfX19hliH+JfBbFujl/s2zCHwQHhVCrQxW6jG2HY6Fcxi5NiHRvzsBFfDrjIz+tteer1qe2EYYwcRmsLHl/QGN+vvwd7Ue8x74/D/Nx8f7M+nQ+92/859EaIYQBPX0i05uLVMguix0fT+jAoouzqP9RLTYu2EZnlz4SHEIYSeSzKPavPpJoGwkMYVT2eXMwYG53Fl387j/B8eBWkLHLEyLd8N18nKePZYQh0oBcBR3+ExxdivZjwchfCH3NMFkI8e4Orz+KbWabRNtIYIhU5Xlw/HR+JpXfK8+Kr/+kk0tfVs1YL5MbCpGCTu45R6mqxRNtI4EhUiVH51yMWNqfH3wnU7isE3MGLeLjYv35Y/o6GXEIkcwe3X/CjfO3KFlFAkOkYUU8CjH571FM3Pg59vmyM3fwYtrn78nsvgu4efG2scsTwiTsXXUIgHJ1SyfaTp7DEGnKRb8rrP5uIztX7CMqMhqvhu68378x5eqWRill7PKESJP6eA8nMiKKecenYmZmJs9hCNNQtFxhhi7qw7KAOXQa04bLR68yosF4eroPYcuSXURHRRu7RCHSlGtnbnDhyBUafFTrtb90SWCINClbrqx8OKY1S6/N4bOfPiU2JpZvOs+mU+E+/DZtLaHBib85TAgRZ/+auGcvarav/Nq2EhgiTctgZUn9LjXxOTmN8etGkMclNz5DfqZDgZ7MHbyYO1fvGbtEIVI1383HKeLhTLZcWV/bVgJDmASlFN6NPJi6/UtmH/6a8o08WP3dRjoX6cuo5l/jt+WETK0uxL88uveYswcuUq5e2SS1l4vewmQ9uBXEunlbWO+zlcf3n5C/WF6a925Avc7VscmY+ANKQpi62NhYRjaayMldZ5lz9BsKFs8HgFJKLnqL9Mc+bw66jG3HsoA5DF3cB9tM1szuu4D2+Xsyf9hSAm/K1CMi/Vo+8U/8/j5B75kfvQiL15ERhkhXzh64wB8z1rP3j4MoMzNqtKtEq0FNcSnrbOzShDCYE7vOMLT2V9RoV5nhS/r94+6oxEYYEhgiXbrjf48/Z25g00/bCX8aQZkabtTuWJXKLcqTOUcmY5cnRIoJefSUHmU/I4N1Bub4Tf7P6VkJDCES8PRxKBvmb2XdvC3cuXoPcwtz3GuXpEbbylRrXREbO2tjlyhEstFaM7HDDPb8cYiZ+8bj6uXynzYSGEK8htaaS0evsvu3A+z+/SB3rt7DNrMNtdpXoVG3OhTxKGTsEoV4Zxt+3Mb07nP5aHx7Oox8/5VtJDCEeANaa87sO8+GH7exa+V+IiOiKFSmIA0+qkWtDlXIYp/Z2CUK8cZO7DrD8HrjKFOzJBPWj8Dc3PyV7SQwhHhLIY+esv2XvWxetINLflexsDTHu0k5aneshncjdzJYZzB2iUK81p2r9+jjPYIs9pmYdWAiGbPaJdhWAkOIZOB/KoDNC3ewY8U+Ht59jF0WW6q1qkitDlUoVa14gr+xCWFMd/zvMbz+eEKCQvju0CTyujgm2l4CQ4hkFBMdw7Htp9m6dBf7/jxMROgzcuTJRo22lanVoQpFPArJzLkiVbh8zJ+RjSYQHRnN+HUjKFHR9bXbpMrAUEqNBCYA32ut+yTSTgH9gZ6AM/AQWKy1Hp7Y/iUwhCFEhD3j4Fpfti/fy5GNx4iOiiFfUUdqta9KzfaVyVc0j7FLFOnU0a0n+arlVOyy2jJp4+cULJE/SdulusBQSlUAlgPBwJ7XBMa3QBNgCHAKyAI4aq03JNaHBIYwtOCHIexddZgdy/dwYudZtNY4lypA5RblqfK+N4VKF5SRh0hxWmv++mEzcwYuokDxvEzcMBL7vDmSvH2qCgylVBbgKNANGA2cTigwlFKuwGmgtNb63Jv0I4EhjOnBrSB2/3aQvasPcXrPebTWOBbKRaXmXlR5rzzFKxaVax4i2T0Lf8bMXvPZ8vMuvBt7MHxJv0QvcL9KaguMX4FrWuthSqmdJB4YQ4GuwFygL3FzX+0Chmit7yfWjwSGSC0e3X/Cgb982fvnIY5vO0VUZDRZc2ahYlNPKjX3wr12SaxsrIxdpkjj7vjf46uWU7l6IoAPx7Sm4xctMTN78+kCU01gKKW6EXctoqLWOjIJgTEX6AKcIO6UlAamxq+uqLWO/Vf77kB3gAIFCpQLCAhIicMQ4q2FBodxZOMx9q0+zOENxwgLCcfKJgMedUtTsakn3o09yJ47m7HLFGlIVGQUf87cwJKxv2FhacHwJX3xblzurfeXKgIj/vTSXqCq1vp8/Hc7STwwfIg7deWqtb4Y/11R4AJQQWt9KKH+ZIQhUrvIZ1Gc3HWWg2t9ObDWl/vXHwDg6lWYCk08qdCkHIXLOsl1D5GgY9tPMbvvAq6fu0WFpuXoPfNjcjvlfKd9ppbA6AIsBGJe+tqcuFFDLGCntX72r22+AkZqrS1f+k4BkUAHrfVvCfUngSHSEq01V08GcHCdH4fW+3H+0GW01mTPnRWPuqXxqBO35HCU0YeAwJtB+Axdws4V+8jtnJNPZ3xExaav/Bn/xlJLYGQF/j3p+kLgEjAROKP/VYxSqh6wGXDRWl+J/64wcBnw1lofTqg/CQyRlj26/4TDG47it+UER7ec5MmDEACcSuanbM2SlKnhRunqJcicXWbWTU/uXrvPr5NXs3nhDlCKtkOb0254i2S9BpYqAuOVnf/rlJRSahJQXmtdO/6zGXAEeAoMiN9sBmAFVPr3NYyXSWAIUxEbG8vVEwH4bTnJ0a0nOLPvAs/CI1FK4Vy6AGWqu1G2VklKVyvxxnfEiLTh5sXbLP/6T7Yu2Y25uRn1utSk3fAW73z66VXSUmAsAmporZ1eauMIzAIaAOHAFmCQ1vpeYvuWwBCmKioyigtHrnBixxlO7DrDmX3niYyIwsxM4eJRiLI14gKkRMWi2GWRAEnLLh/3Z+WUNez6dT8WGSxo3L0urT9rhkO+pD9X8aZSbWCkJAkMkV5EPovi/KFLHN9+muM7TnPu4EWio2JQSuFUMj9ulVxxq1wMt0qu5HbOKRfRU7mwkHB2rtjHhh+3cuHIFaztrGj2aQNaDWpCtlxZU7x/CQwh0pHw0AjOHbjImX0XOHPgAucOXCQsJByAbLmy4FrehWJeRXAt74KrV2EyZcto5IpFVGQUfn+fZOfKfexffYTwpxE4lcxPo0/qUPuDqga9ViWBIUQ6FhMTQ8CZm5zZd55zhy9x/tBlbpy/9WJ9vqKOFPUsjEtZZ1w8nHFxd5YQMYCwkHBO7jrLvj8PsW/1YUIehZIpmx1V3q9Ag661KO5dxCijQQkMIcQ/hD4J5YLvVS4cvsz5w5e4dPQqgTeCXqzP7eRAYXdnXMo6U9AtPwVL5COvS27MLWQ6k7cVEx3DhSOX429eOMm5g5eIiY7BNpMNlVp4UaNNJTzqlsYyg+Xrd5aCJDCEEK/15EEwl4/5c+moP1eOx/1569KdF+stLM3J55qHgiXyUbBEXIgUKJ6PPC65yWBl3B9yqVHok1DOHbrM2f0XOHvwIucOXiQsOBylFC4ezpSrUxqPuqVxq1wsVf3zk8AQQryV8NAIrp+7xfWzNwk4e4OAczcJOHOTu/73ef6zw8xMkcspJ/lc85CviCN5Cucmj0tu8hTORS4nB6P/xmwIwUEhBJy9ScDZm1w+epUzBy4QcOYmWmvMzBROJQtQvEJR3GuVxL12KTLnSL3PzyQWGBaGLkYIkXbY2Fnj6lkYV8/C//g+IuwZN87filsu3ObmxdvcuHCbU7vPEhH6/wkblFJkd8yKQ357HPLnIGe+HDjkt8c+Xw5y5MmGfd7s5MiTLU2ESlRkFPeuBXLz4h1uXrzNzQu3uXHxNtfP3eLx/Scv2tllsaVExaJUa1URt0quuJZ3wS6zrRErTz4ywhBCJButNY/vP+HW5bvcuXKP21fuEngjiMCbD+L+vBFERNiz/2yX1SEz2fNkI1O2jNhlscUuqy12mW2xzWSDtZ01Nhmtsbazivvzpb/HfW+NpZUFllaWcUsGixcXi7XWREVGExkeybPwSCIjInkWFklYcBihweGEPYn7M/RJGGHBYYQFh8etCwknLDic4KAQQoJCCA56+uJOs+cy58hEvqKOFCieL/40Xdxiny/HW80Sm1rICEMIYRBKKbLlykq2XFkpWbnYf9ZrrQl5+JQHtx4SdPshD27FLUG3HhJ09xGhj8O4c/UeoU/CCH0SRnhIOLGxb/5LrWUGC8zMzYiMiOJNfim2yWiNbWYbbDPbYpfZhqwOmSlQLC+Zsmckc45M5CroQL6ijuQrmidVn1ZKKSY7wlBKBQKGmN/cHnhggH5So/R87JC+j1+O3XQV1Fo7vGqFyQaGoSilfBMavpm69HzskL6PX449fR572j3RJoQQwqAkMIQQQiSJBMa78zF2AUaUno8d0vfxy7GnQ3INQwghRJLICEMIIUSSSGAIIYRIEgkMIYQQSSKB8RKl1KdKKX+lVIRSyk8pVTWRttZKqUVKqZNKqaj4182+ql0GpdTY+P0+U0pdV0r1S7GDeAcpdPwdlFLHlVJhSqm7SqmlSqncKXYQb+kNj72GUmqNUupO/HGdVEp9/Ip21eP3FaGUuqqU6pmyR/F2kvvYlVLvK6X+VkoFKqVClFKHlFLNUv5I3k5K/Lt/qX0VpVS0Uup0ylRvWBIY8ZRSbYGZwETAHdgPbFRKFUhgE3MgApgNrE9k18uJex95d8AVaA2cTKayk01KHL9SqjKwBFgMuAEtgBLAsmQt/h29xbFXAk4BrYCSwBzARynV4aV9OgMb4vfl/r/27j1UijKM4/j3Z0UUJRRmGmb3oiIyi6ASE8roCtUBu0EX6o/KLhASlVRGCJTaOwAABdBJREFUQYlJYoVKFyUqoiv4R5FQHakQTYukC2SQoSezC1FYkRZPf7xzZB336OzuzM5mvw8MuzvzzrPznD173pl33vO+wMPA45L6qsqjHVXkDpwJvAtckMV8E3hjR3+I61JR/oOx9wOeA96p4NDrERFeUk+x5cBTuXVrgIcL7PsE0N9k/TnAr8CIuvOrKf9pwLe5ddcBm+rOt6zcG8q/DLzW8HomsCZX5mlgWd35Vp37EGVWALPrzreb+QOvA/cDM4DP6s61jMVXGKRmI+BkYElu0xLSGUW7LgY+Au6QtF7SGklzJfXU/JcV5v8hMFrSRUpGAJeTzjh7Qom5Dwd+aXh9WpOYbwOnSOqJsbwrzL2ZfQuU6aoq85d0MzAKeKiTY+w1rjCSEaQmlo259RtJH3q7DgcmACcCfcAtpOapRR3ErEIl+UfEMuAKUhPUZuBHQMA17casQMe5S7oQOItt/6Fr1BAxd8/esxdUlXu+zFRgDKl5spdUkr+kE0hXFldFxD/lHGpvcIWxrfx/MarJulYMy/a/MiKWR8TbpEqjT9KBHcStSqn5SzoOmAs8SDqTO5f0RVzQbswKtZV7dp/mReC2iFhRIGaz9XWrIvfBMn3ALNIfz26MHt2O0vKXtCfwEjAtIr4p+0Dr5vkwkp+Af9j+rGIk2599tGIDMBARvzas+zJ7HNth7DJVlf/dwIqImJW9Xi3pd+B9SdMjYl0HscvSdu6SJpCa1+6LiHm5zd8PEfNv4Oe2j7ZcVeU+WKaPdFVxdUQs7vxwS1dF/qNJHTsWSlqYrRuWdtHfwPkRkW8C+8/wFQYQEZuBVcDk3KbJpF4T7foQOCh3z+Lo7LFnzrYqzH9v0hey0eBr0QPazV3SROAt4IGImNOkyDLg7CYxV0bElvaPuDwV5o6kKcDzwLUR8Wo5R1yuivIfAE4AxjUs84Gvs+edfJ/qV/dd915ZgMtI7ew3AMeSutptIk0mAqlb5Du5fY4j/RK8BKzMno9r2L4PsA54hdSt9AzgM+CVuvPtUv7XAluAm0j3c84gdQJYVXe+neQOTAJ+JzW1jGpYDmgoc1hWZk4W84bsPfrqzrcLuV+efe6358rsX3e+3ci/yXvMYBfpJVX7AfTSAtwMrAX+Ip15TGzYtghYmyu/ltTWuc2SK3MMqdfFH6SzjyeBfevOtYv53wp8nuW/gdTmO6buXDvJPXu9Xd5Nfj5nAh9nMb8Bbqw7z27kDvQPUaa/7ly79dnn4s9gF6kwPFqtmZkV4nsYZmZWiCsMMzMrxBWGmZkV4grDzMwKcYVhZmaFuMIwM7NCXGGYdSibVCey0XjNdlmuMMxaJKlf0hP/lbhmZXGFYWZmhbjCMGuBpEWkIT+mZs1QARyabT4xm7/6D0krJY3P7Xu6pKXZ9gFJ8yQNHyqupEMl7SbpmWzO6T+zSbjulOTvrnWdf+nMWnM7aSTahaShrEeTBpiENFDdXcB40hDmL0gSbJ1UZwmwmDSh1qWkwRqf3UncYaQxyKaQBsebDtxDmurWrKs8lpRZiyT1kwaTuyV7PQl4Dzg30iRZg5PrfAAcHBHrJT0HbImI6xvijAM+AQ6MiB/ycXfw/o8Ap0REfvh0s0p5AiWz8qxueP5d9jgSWE+acfBISZc1lBmcE+QI4Iehgkq6kTT89iHAXsAe9NB8Kvb/4QrDrDyNEyMNXroPa3h8GnisyX4DQwXMKpg5wDTS5Du/AVOBSzo9WLNWucIwa91mYLcW9/kYOD4ivm4x7gRgeURs7W4r6YgW39usFL7pbda6tcCpWS+mERT7Hs3M9pkv6SRJR0q6UNKCoeJmPaG+AsZLOk/SUZLuJfWmMus6VxhmrXuUdDXwBfAjMHZnO0TEamAiqQvuUuBTUq+qjTuJuwB4mTRT4UfZ/rNLycKsRe4lZWZmhfgKw8zMCnGFYWZmhbjCMDOzQlxhmJlZIa4wzMysEFcYZmZWiCsMMzMrxBWGmZkV8i9anwt0QyN1KAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Make a contour plot of lnL = lnLmax - 1/2 (here for theta and xi)\n", "if not(m.is_fixed('theta') | m.is_fixed('xi')):\n", " plt.figure()\n", " m.draw_mncontour('theta', 'xi', nsigma=1, numpoints=200);\n", " plt.plot(MLE[0], MLE[3], marker='o', linestyle='None', color='black', label=r'$(\\hat{\\theta}, \\hat{\\xi})$')\n", " plt.figtext(0.6, 0.8, r'$\\ln L = \\ln L_{\\rm max} - 1/2$')\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEFCAYAAADOj31RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXgNZ//H8fedkMW+JZGURCTSWhohIfaldloaa1Gttij1VHctHpQWraofpaWqtLW1PGqnsZMQErHE1toitRNLlWyW+/dH0lyJJgTnZOYk39d1zSXOmdzzGcPHZM4sSmuNEEII87EzOoAQQoisSUELIYRJSUELIYRJSUELIYRJSUELIYRJFTA6gCWVKVNGV6hQwegYQgiRY9HR0fFaa5es3stTBV2hQgV27dpldAwhhMgxpVRcdu/JIQ4hhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDCpPHWanXg0Wmtu3rzJlStXKFSoEKVKlcLOTv7vFsJoUtB5UEJCAhcuXOD8+fNcuHCBCxcuEB8fz5UrV7hy5QqXL19O//qfycHBgRIlSpCUlMT169cpXbo0bm5uuLm54erqmv71vb93cXGhYMGCRq+yEHmS6QpaKXUS8MrirdVa63a5HMe0UlJS2Lt3L9u2bSMyMpLTp0+nF3JKSgply5bFzc0t/VcXFxc8PDyoVq0apUqVyjSVLFkSJyen9LFv3brFpUuX0sv94sWL6V/HxMRkei0+Pp5ixYplKnB3d3cqVKiQaSpevLiBf1pC2CbTFTRQC7DP8Ht3IBpYaEwcc4iPj2f79u3p0+7du/Hx8aF+/fq0bt0ab2/v9EIuVqwYSqlHXlbBggXx8PDAw8PjgfPevXuXy5cvZyrxM2fOcOzYMdavX09cXByxsbG4urpSu3bt9KlGjRoUKlTokTMKkR8osz9RRSk1DPgA8NBaJ9xv3qCgIJ2XLvWOjY1l0aJFLFq0iCNHjlCnTh3q1atHvXr1CA4OplixYkZHzJG7d+9y5MgRIiMj2blzJ5GRkRw6dIgnn3wyU2lXrlwZe3v7Bw8oRB6ilIrWWgdl+Z6ZC1ql7gYeJ/Xwxn+ymacf0A/A09MzMC4u28vabcKff/7Jzz//zMKFC/nzzz/p2LEjXbp0oXHjxhQoYMYfeB5NUlIS+/btIzIyMn06e/YsgYGBNG/enFatWlGzZk0pbJHn2XJBtwRCgRpa670Pmt9W96CTkpJYunQps2bNIjo6ms6dO9O1a9c8V8oPcvXqVXbs2MG6desIDQ3lwoULtGjRglatWtGyZcscHXIR4n4uXbrE119/zcCBA3FxyfIGclb53vux5YJeBHhprWvnZH5bK+i4uDgmTpzI3LlzqVmzJq+++iohISGZPrDLz06dOsXatWsJDQ1l/fr1lCtXjlatWtG2bVsaNmyYr/7zEpbRuXNnkpKSKFSoEAsXPtzHWo/zvfdzv4I27cmuSilXoAPwndFZLO348eP06dOHGjVq4OjoSHR0NOvWraN79+5SzhmUL1+e1157jYULF3Lx4kVmzJhB4cKF+eCDDyhbtiwvvfQSixcv5saNG0ZHFTZg/vz5ODo6snLlSgoWLPhQJfs43/s4TLsHrZQaDAwn9cPBv3PyPWbfgz5//jzDhw9nyZIlDBw4kEGDBlG6dGmjY9mkU6dOsXz5cpYtW8aOHTto2LAhHTp0oH379pQtW9boeELkmM3tQad9ONgH+Dmn5WxmCQkJfPLJJ1SrVo0SJUpw9OhRRo0aJeX8GMqXL8/AgQNZu3Ytp06dolevXmzatInKlSvTrFkzZs+ezV9//WV0TCEeiykLGmgCVCIPHN4IDQ3l6aefZv/+/URFRfHFF19QsmRJo2PlKcWLF+eFF15gwYIFnD17lgEDBrB8+XI8PT3p0qULS5cuJTk52eiYQjw0Uxa01nqT1lpprSONzvKoLly4QI8ePRgwYABff/01CxcuxNvb2+hYeZ6zszOdO3dmyZIlnDx5kpYtW/J///d/eHh48MYbb7B7926jIworu3r1Km5ubhw/fjzT61OnTsXX1xdnZ2dat27NpUuXcjzmq6++ilIq01SvXj0g9cPDiRMnWnQd0mmt88wUGBiojXbnzh393XffaRcXFz148GB98+ZNoyMJrXVcXJweNWqU9vT01DVr1tTTpk3T165dMzqWsIL3339f9+7dO9NrQ4cO1V5eXnrDhg163759ulKlSv+a5362bNmivby8dNeuXfWiRYv0nj179KVLl7TWWsfExOiSJUs+8t8nYJfOptMML1VLTkYX9LVr13SLFi107dq19d69ew3NIrJ2+/Zt/dtvv+nOnTvr4sWL6969e+vw8HB99+5do6MJC7h586YuUaKEDgsLS38tKipKK6X09u3b01+bPHmydnFxydGYKSkpunz58vrLL7/Mdp7AwEA9derUR8p8v4I25SEOW3Tu3DkaN26Mn58f27dvp3r16kZHElmwt7enVatW6ZfPV6lShVdffZWqVasyceJErly5YnTEPGHevHlUqFABOzs7KlSowLx583JluatXr8bOzo769eunvzZhwgQaNWpE3bp1019zcXEhPj4+R2Pu37+fM2fO8MYbb2Q7T/v27VmwYMGjB8+GFLQFHD16lPr169O5c2emTJkilyfbCFdXVz744AN+//13pk+fzp49e/Dx8aF///4cOnTI6Hg2a968efTr14+4uDi01sTFxdGvX79cKemwsDACAwPTbxZ269YtVqxYQceOHTPNl5iYmOM7LJYoUYK7d+8yZswY4uLiuHv37r/mqV27NpGRkSQmJj7+SmSU3a61LU5GHOKIiorSZcuW1d99912uL1tY3vnz5/XHH3+s3dzcdMuWLfWqVav0nTt3jI5lU7y8vDTwr8nLy8vqy+7QoYN+6aWX0n8fGRmpAe3k5KQLFy6cPjk4OOjg4OAcj/vtt99qJycnDWillN69e3em9/ft26cBfezYsYfOjBzisI5169bRpk0bpk+fTp8+fYyOIyzAzc2NkSNHEhcXR48ePRg2bBiVK1fmm2++kSsWc+jPP/98qNctKTExMdPVuH/88QcODg7ExMSwd+/e9Mnf3z/TYZD7+fLLLxk2bBjvvPMOoaGh7NmzB39//0zzODs7py/forJrbluccnMPev78+drV1VVv3bo115Ypct/du3f1li1bdEhIiC5durT+6KOP9Pnz542OZWpG7kH36NFDd+nSJf33U6dO1R4eHpnmuXTpki5QoECmDw2zExERoQsUKKAPHDhw3/l27NihgUf6u4HsQVvW5MmTGTx4MOvXr6dhw4ZGxxFWpJSiUaNG/Prrr0RGRnL9+nUqV67Mm2++ia3f2tZaxowZ86+HMRQqVIgxY8ZYfdk1atTI9PlBmTJl+PvvvzMdNx43bhx169bN9KFhdn777TeqVKlC1apV7zvfgQMH8PDwwM3N7dHDZyW75rbFKTf2oGfNmqV9fX31yZMnrb4sYU7nzp3TH374oS5VqpR+6aWX9KFDh4yOZDpz587VXl5eWimlvby89Ny5c3NluTExMdrOzk7Hx8drrbW+ePGidnZ21qNHj9axsbH6iy++0C4uLvrIkSM5Gm/mzJnazs5Of/bZZ/rAgQP6ypUrWc738ssv61dfffWRMiPnQVvGvn37dJkyZfTBgwetuhxhG65evao//fRT7erqqkNCQnRkZKTRkYTWuk6dOpnOSV60aJH29PTUzs7OumXLlln+hzp79mwN6NjY2Eyv37lzR48dO1ZXqVIl/UPCjh07ZponMTFRFytWTEdERDxSXiloC7h+/br28/PTP/30k9WWIWzTzZs39eTJk3X58uV18+bN9caNG+XCFwOtWbNG+/n56du3b+f4e0aMGKGrVKmib926dd/51q5dqwH9119/pb82depU3aJFi0fOe7+ClmPQOTRgwAAaNWpEr169jI4iTKZQoUIMGjSIY8eO0aNHD/r370+TJk0IDw83Olq+1Lp1awYOHMjp06dz/D2rV69m6tSp930IREJCAtu2beOpp57K9DzQggULMmXKlMfKnK3smtsWJ2vtQS9cuFD7+fnJfTVEjty6dUvPnj1be3l56TZt2ujo6GijIwkLmD17tq5Vq5bet2+fRcflPnvQpr1h/6Owxg37L1y4QPXq1Vm2bBnBwcEWHVvkbcnJycycOZMxY8bQoEEDxo4di6+vr9GxhMnY3A37zUJrzeuvv85rr70m5SwemqOjIwMHDuTo0aMEBARQp04dBg0a9FC3uRT5mxT0fSxbtowjR44wYsQIo6MIG1a4cGGGDh3K4cOHAahcuTJjx44lISHB4GTC7KSgs5GUlMR7773HV199haOjo9FxRB7g4uLCV199xY4dO9i7dy9PPvkks2fP5s6dO0ZHEyYlBZ2NiRMn4u/vT/PmzY2OIvIYX19fFi5cyKJFi5g1axYBAQGsWbOGvPR5kLAM+ZAwC6dPnyYgIIDIyEgqVqxogWRCZE1rzfLly/nwww954oknGD9+PIGBgUbHErlIPiR8SB9++CH9+/eXchZWp5SiQ4cOHDhwgK5du/Lcc8/Rs2dPTp48aXQ0YQJS0PfYtm0bW7duZciQIUZHEflIgQIFeP311zly5Ah+fn4EBgYyZMgQbt68aXQ0YSAp6Ay01rz11lt8/vnnFC5c2Og4Ih8qUqQII0eO5MCBA5w6dYqqVauybNkyOT6dT0lBZ3DixAnOnj1L9+7djY4i8jl3d3fmzp3L7Nmz+eijj2jfvj2xsbFGxxK5TAo6g82bN9O0adP055kJYbSmTZuyb98+6tWrR61atRg7dizJyclGxxK5RAo6g82bN9OkSROjYwiRiYODA0OGDCEqKoqIiAiqV6/Oxo0bjY5lSn///Tdvv/02Xl5eODs7U69ePaKiotLf7927N0qpTFOdOnUyjfHuu+9SqlQpypcv/68H3a5YsYIGDRrk3iGn7G7SYYvT49ws6fbt29rNzU0fP378kccQIjcsXbpUe3p66h49euhz584ZHcdUunbtqp966im9adMmffToUT1y5EhdrFgxffr0aa116o31mzdvrs+dO5c+Xb58Of37ly9frt3c3HRUVJSeP3++dnJy0pcuXdJap95y2MfHx+L3g0duN/pgO3fuxMXFRU6tE6bXoUMHDh06hKenJ08//TRTp06VqxFJfWDr4sWL+eyzz2jSpAm+vr58/PHH+Pr6Mm3atPT5HB0dKVu2bPpUqlSp9PcOHz5MkyZNCAoKonv37hQrViz92P/QoUN58cUXqVKlSq6tkxR0miVLlhASEmJ0DCFypHDhwowbN44tW7bwv//9j9q1a7Nv3z6jYxnq9u3b3LlzJ9NTvSH1idsZ780dHh6Oq6srfn5+9O3bl4sXL6a/V716dXbt2sXVq1eJjo4mMTERX19fduzYwaZNmxg6dGiurQ8ghzj+4efnJ/ftFTbp7t27etasWdrFxUWPHDlSJycnW3wZZPGUbmtPj6Ju3bq6QYMG+vTp0/r27dt6zpw52s7OTvv5+WmttV6wYIFetmyZjomJ0cuXL9f+/v66atWqOikpKX2MkSNHah8fH12tWjX966+/6pSUFO3v76/DwsL09OnT9VNPPaVr1qypt23bZqk/W3nk1f1cvnxZFy1a9KEekSOE2Zw+fVq3a9dO+/v76927dxsdxxDHjh3TjRo10oC2t7fXtWrV0j179tSVK1fOcv4zZ87oAgUK6MWLF2c75ieffKJff/11HRMTo11dXfXZs2d1aGio9vDwsMh/hvcraDnEAezevZsaNWpgb29vdBQhHtkTTzzBihUreO+992jVqhXDhw/Pd6fk+fj4sGXLFm7cuMGpU6eIjIzk1q1beHt7Zzm/h4cH5cqV4+jRo1m+f+TIEWbNmsXnn3/Opk2baNSoEe7u7rRs2ZKUlBT++OMPa66OFDTArl27CArK8l4lQtgUpRQvvfQSe/fuZd++fQQFBWHppwzZgsKFC+Pu7s7Vq1cJDQ2lQ4cOWc4XHx/PmTNncHd3/9d7Wqc+sGPChAkUL16cu3fvcuvWrfT3bt26ZfUPZ6WggfXr19OwYUOjYwhhMR4eHixbtoyPPvqIdu3aMXDgQC5fvmx0LKsLDQ1lzZo1xMbGsm7dOpo2bcqTTz7JK6+8wo0bN3j//feJiIjg5MmTbN68meeeew5XV9csTxD4/vvvKVGiBB07dgSgQYMGbNy4kfDwcKZNm0bBggV58sknrbtC2R37sMXpUY5BX716VRctWlTfuHHjob9XCFsQHx+vBw4cqF1cXPSUKVP0rVu3jI5kNb/88ouuWLGidnBw0GXLltUDBw7U165d01prnZCQoFu2bKldXFx0wYIFtaenp3755Zf1n3/++a9xzp8/r728vNLPn/7H2LFjdZkyZbS3t7des2aNRTIjD43N3vz581mwYAErVqywUiohzGH//v28/fbbXLp0iUWLFll/70/kiNwP+j62b99O06ZNjY4hhNU9/fTTrF+/njfffJNGjRqxbt06oyOJB8j3BV2oUCFSUlKMjiFErlBK0bdvXxYtWkSvXr2YPHkyeemn6Lwm3xd0yZIluXr1qtExhMhVjRo1IiIigjlz5tCuXTvOnTtndCSRBSnokiW5du2a0TGEyHXe3t5EREQQFBREjRo1+PXXX42OJO6R7ws6JSUFBwcHo2MIYYiCBQsyevRoli5dyuDBg+nTp488ZstE8n1B//XXX5QoUcLoGEIYqk6dOuzZs4eUlBSCgoKIiYkxOpJACpr169cTHx9vdAwhDFe0aFF++uknhg4dSrNmzfj666/lA0SD5fuCjo6OlkfcC5FBr1692L59O7NmzSIkJIQrV64YHSnfyvcFPWzYMPz9/Y2OIYSpVKpUie3bt+Pj40NAQABbt241OlK+ZMqCVkq5K6V+VEpdUkolKaUOKaUaW2NZbm5uXLhwwRpDC2HTHB0d+fLLL5k+fTrdunXj448/5vbt20bHyldMV9BKqRLANkAB7YDKwJvAxft936Nyc3Pj/Pnz1hhaiDyhbdu27N69m/DwcJ555hlOnTpldKR8w3QFDQwGzmmtX9JaR2qtY7XWG7TWh62xMCcnJ5KSkqwxtBB5hru7O6GhobRp04ZatWrJU8VziRkL+nlgp1LqF6XURaXUXqXUf5RSyhoLO336NOXKlbPG0ELkKfb29gwZMoR58+bRo0cPuUw8F5ixoCsCbwAngFbAZOAzYGBWMyul+imldimldl26dOmhF3bq1CnKly//GHGFyF+aNWtGREQEs2bN4pVXXpGfQK3IjAVtB+zWWg/RWu/RWs8GviKbgtZaz9BaB2mtg1xcXB56YVLQQjw8b29vtm/fTmJiIo0bN+bMmTNGR8qTzFjQ54BD97x2GPC0xsIcHByIi4uzxtBC5GmFCxfm559/pmPHjtSuXZuwsDCjI+U5ZizobcC9dxL3A6zSoh999BHff/89x44ds8bwQuRpSik+/PBDZs6cSbdu3XjvvfdITEw0OlaeYcaC/j+gjlJqmFLKVynVBRgEfG2NhT3xxBMMGDCAqVOnWmN4IfKFNm3aEBMTw9mzZwkICCAiIsLoSHmC6Qpaax1F6pkcXYEDwBhgOPCNtZbZrFkzIiMjrTW8EPlCmTJlWLBgAePGjSMkJISFCxcaHcnmFTA6QFa01quAVbm1vMDAQPbu3cudO3ewt7fPrcUKkSd17NgRX19fWrZsiaOjIx06dDA6ks0y3R60EYoWLYqrqyuxsbFGRxEiT/D392fVqlX07duXJUuWGB3HZklBp6lSpQqHDt178ogQ4lEFBgayZs0a3nzzTcaPHy8XtTwCKeg0UtBCWF5gYCA7duxgwYIF9OnTRx7Q/JCkoNNUrVqVgwcPGh1DiDynXLlyhIWFcenSJVq3bi33l34IUtBpZA9aCOspUqQIS5YsoUaNGtStW5ejR48aHckmSEGnqVKlCseOHZNbKQphJfb29nz55Ze89957NGzYkC1bthgdyfSkoNMULVqUd999l0GDBhkdRYg8rV+/fsydO5cuXbrwww8/GB3H1KSgM/jwww/ZvXs3+/btMzqKEHla8+bN2bJlC5988glDhgzh7t27RkcyJSnoDJycnKhevbqcDy1ELqhcuTI7d+4kPDycrl27kpCQYHQk05GCvoeHhwfnzp0zOoYQ+UKZMmVYv349zs7ONG7cWB4/dw8p6Hu4u7tLQQuRixwdHfnpp59o27YtzzzzjDzEOQMp6HtUqFCBzZs3c+vWLaOjCJFvKKUYNWoU3bp1o1mzZly8aJVnRNscKeh79OjRgyJFitC/f3+5NFWIXDZixAg6duxI06ZNOXv2rNFxDCcFfY+CBQuyaNEiYmJi+OYbq93hVAiRBaUUo0eP5sUXX6Rhw4b5/gN7KegsFC5cmJkzZ/Lpp5/KJ8tCGGDIkCG8++67NGrUKF/fq10KOhvVq1enTp06zJ071+goQuRLAwcOZPLkyTz77LN89dVX+fKQoxT0fTRt2pSYmBijYwiRb3Xs2JGIiAh+/PFHunTpws2bN42OlKukoO/Dz8+PyMjIfPk/txBm4ePjw/bt2ylQoABvvPGG0XFylRT0fTRv3py7d+/y/fffGx1FiHzN0dGR77//nqioKH788Uej4+QaKej7KFCgAD/88ANDhgyRe9gKYbDChQuzcOFC3n//fTZv3mx0nFwhBf0A1apVo3HjxvJcNSFMoFq1aixatIguXbrki5KWgs6BF154genTp5OUlGR0FCHyvSZNmqSX9Lp164yOY1VS0DkQEhKCr68vISEhUtJCmECTJk349ddfefHFF1m0aJHRcazmgQWtlApXSpXKjTBmZW9vz5w5cyhWrBgvvvii0XGEEEDDhg1Zu3Ytb7/9Nt9++63RcawiJ3vQ9YByGV9QSlVXSql7Z1RKFbFUMLMpUKAAc+bMITo6mm3bthkdRwhB6gVlW7duZcyYMcyfP9/oOBaX00Mcvv98kVbM0UC1jDMopRYBfymlDiqlqlouonk4ODgwbNgwxo4da3QUIUQaHx8f1qxZwzvvvJPnjknntKA7Z/i6fNr3uf3zglKqBNARaAUsAvLsicMtWrSQqwuFMJmqVauyePFievbsmad+ws1pQddRSr2llHIA+gBJQOMM73sASVrr9cDnwDTLxjQPNzc3Ll68KLdCFMJkGjRowNy5cwkJCWH79u1Gx7GInBT0bOBl4F3gJjAMeAvop5R6Km2eZ4ETAFrrRK11nr3Ux8nJiVGjRhEcHExUVJTRcYQQGbRs2ZI5c+bw/PPPs3fvXqPjPLYCD5pBa/0agFLKB3gauKa1jlVKFQNilFIHST0e/V+rJjWRjz76iMqVK9OuXTu+//57nnvuOaMjCSHStGrViq+++oquXbsSHR1N0aJFjY70yNTj3AhIKVWH1OPOp4DZ2uC7CgUFBeldu3bl2vIiIyNp164dmzZtolq1ag/+BiFErunXrx83b95k7ty5ZHHSmWkopaK11kFZvfdYF6porXdorUdprWcZXc5GqF27Nl9++SXPP/88iYmJRscRQmQwadIkDh06xLhx44yO8sjkSsLH9NJLL+Hh4cGGDRuMjiKEyKBQoUKsWrWKGTNm8NNPPxkd55FIQVvAs88+m6cvNxXCVnl4eLB69WoGDx7MrFmzjI7z0KSgLaB3795EREQwYcIEo6MIIe5RpUoVtmzZwrhx4xg8eDB37twxOlKOSUFbgKurKxs3bmTatGlMmTLF6DhCiHs8+eST7Nixg8jISHr16mV0nByTgraQcuXKsXHjRj7//HMWL15sdBwhxD1Kly5NaGgoERERNnMhixS0BXl5ebF8+XL69++frx8VL4RZOTo6MnjwYD799FOjo+SIFLSF1axZk2nTptGrVy859U4IE3rllVc4fvw4P/zwg9FRHkgK2go6d+6Mv78///1vvrm4Ugib4eTkxJIlS/jggw/IzQvbHoUUtJV8/fXXrFq1ihEjRpAPr+HJE+bNm0eFChWws7OjQoUKzJs3z+hIwkKqVKnC9OnT6dixI3/++afRcbL1wHtxiEfj6urK1q1badWqFTdu3ODLL7809eWmIrN58+bRr18/EhISAIiLi6Nfv34A9OzZ08howkI6derEqVOnaNGiBWFhYbi6uhod6V9kD9qKXF1d2bRpExs3bmTixIlGxxEPYdiwYenl/I+EhASGDRtmUCJhDW+//TZdu3alTZs2pnzeqBS0lZUoUYLly5czceJEli1bZnQckUPZ/dhr5h+HxaMZPXo03t7ejBgxwugo/yIFnQs8PT1ZunQpffv2ZeXKlUbHETng6en5UK8L26WUYtq0acyZM8d050ebrqCVUh8rpfQ903mjcz2uWrVqsXLlSl577TXZk7YBY8aMoVChQpleK1SoEGPGjDEokbAmFxcXpk6dSu/evf91aMtIpivoNH8A7hmmp42NYxm1a9dm1apV9OnTR55raHI9e/ZkxowZeHl5oZTCy8uLGTNmyAeEeVinTp0ICgoy1ecMj3XDfmtQSn0MdNZaP/Qd8HP7hv2Pat68eXz88cfs3LmTUqVKGR1HCJHmypUr+Pv78/3339OqVatcWabVbthvRRWVUmeUUrFKqZ+VUhWNDmRJPXv2pGvXrtSpU4fDhw8bHUcIkaZUqVLMnz+fl19+mZMnTxodx5QFvRPoDbQB+gJlge1KqdJZzayU6qeU2qWU2nXp0qXcS/mYxowZw5AhQ2jcuDGrV682Oo4QIk2jRo346KOP6NSpE7dv3zY0i+kOcdxLKVWE1CeGf6a1vu/JxLZyiCOjiIgIQkJCmDZtGiEhIUbHEUIAWmuaNWvGCy+8kH6BkrXY4iGOdFrrG8BBoJLRWayhbt26rF69mv79+8vZHUKYhFKK8ePHM2rUKG7evGlYDtMXtFLKCXgKOGd0FmupWbMmK1eupH///nzzzTdGxxFCAEFBQTRt2pTRo0cblsF0Ba2UmqCUaqyU8lZKBQP/AwoDPxoczapq1apFeHg4U6ZM4a233rKpx/IIkVdNnDiRH374gd27dxuyfNMVNFAOWEDqudC/AslAHa11nKGpcoGPjw8RERHs2xrAjXAAABNgSURBVLePAQMGyF3whDCYq6sr48eP5/XXX+fu3bu5vnzTFbTW+gWttYfW2kFr/YTWupPW+pDRuXJLiRIlWLFiBXv37mXo0KFS0kIYrFevXiilWLhwYa4v23QFLaBo0aKsXr2aNWvW0K1bN65evWp0JCHyLTs7O8aPH8/QoUNz/SlJUtAmVaZMGXbs2IGHhwfVq1dn8+bNRkcSIt9q0qQJtWvXZvjw4bm6XCloE3NycmLSpEnMmDGD7t27M23aNKMjCZFvTZkyhXnz5uXqHe+koG1A69atCQ8PZ9KkSXzwwQeGfFghRH7n4uLC559/nqt70VLQNuKfMzyioqJo3bo1Fy5cMDqSEPlOt27diImJybX7dEhB25BSpUqxfv16goODqVGjBmvXrjU6khD5iqOjIz169ODbb7/NleVJQduYAgUK8MknnzB37lx69+4tx6WFyGXvvPMOM2bM4MqVK1ZflhS0jXrmmWcICwtjwoQJfPrpp3K+tBC5pEKFCnTs2JFJkyZZfVlS0DbMx8eH8PBwFi9eTPv27YmLy/MXWwphCu+99x7fffcdt27dsupypKBtnLu7Ozt37qROnToEBgbyxRdfWP0vjRD53VNPPUXFihVZs2aNVZcjBZ0HODg4MGzYMHbu3Mm6deuoX78+sbGxRscSIk8bMGAAEyZMsOrhRSnoPMTHx4fQ0FC6d+9OcHAwv/76q9GRhMizunfvTnx8PL/99pvVliEFnccopXjnnXdYuXIl77//Pr179+by5ctGxxIiz7G3t2fUqFF89tlnVluGFHQeVbt2bWJiYihRogTVqlXj559/ljM9hLCwli1bEh0dbbX7t0tB52FFihRh0qRJLF26lDFjxvDCCy9w48YNo2MJkWcUL14cNzc3jhw5YpXxpaDzgeDgYKKioihSpAjBwcFW+8skRH7UpEkTq13VKwWdTzg5OTFz5kzeeust6tWrx5gxYwx9GKYQecXzzz/P0qVLrTK2FHQ+opSiX79+REREEBMTg5+fH9999x23b982OpoQNqt58+bs3LmTpKQki48tBZ0PVapUiV9++YWlS5cyf/58AgICcvUet0LkJc7Ozvj6+nLw4EGLjy0FnY/VqlWLjRs3MnLkSDp37sx//vMfrl+/bnQsIWxOjRo1iIyMtPi4UtD5nFKKLl26cPDgQZKSkqhatSo//vijPBRAiIfQvn17Fi9ebPFxpaAFACVLlmTmzJksXLiQGTNmEBgYyPr1642OJYRNaNu2LdHR0RZ/kIYUtMikbt26hIeH89///pcBAwbw/PPPc+7cOaNjCWFqzs7OlC1b1uJX7UpBi39RStGpUycOHjyIv78/AQEBzJs3T65EFOI+rl27RokSJSw6phS0yJaDgwOjR49m9erVjBs3jtatW7N161YpaiGykJycjKOjo0XHlIIWDxQYGEh0dDRdunShT58+1K9fn+XLl8sHiUJkULJkSa5du2bRMaWgRY44OjrSp08fDh8+zLvvvsuoUaPw9/dn6dKlskctBKn35ZCCFoayt7enc+fO7Nq1i/HjxzNixAjq169PeHi40dGEMNSNGzcoUqSIRceUghaPRClF27Zt2bNnDwMGDODFF1+kXbt2REREGB1NCENcvXqVkiVLWnRMKWjxWOzt7enVqxe///47zz77LD169OCZZ55h/fr1cuhD5CvOzs4Wvx+HFLSwCCcnJwYMGMCRI0d45ZVXGDRoEHXq1GHz5s1GRxMiV5QqVYr4+HiLjikFLSyqYMGC9OrViwMHDvDuu+/y6quv0qFDB/744w+jowlhNVevXiU2NpYKFSpYdFwpaGEVdnZ2dOvWjUOHDtGgQQPq169Pnz59iI6ONjqaEBb3yy+/0KpVK0qVKmXRcaWghVU5OTnxwQcf8Pvvv+Pt7U2nTp0IDAzk22+/5e+//zY6nhAWsWLFCrp27WrxcaWgRa4oU6YMw4YN48SJE4wdO5a1a9fi6enJ66+/zuHDh42OJ8Qju3XrFuHh4TRt2tTiY0tBi1xlZ2dHq1atWLx4MYcOHcLd3Z0mTZrQrl07Nm7cKGd+CJvz22+/UaVKFUqXLm3xsaWghWHc3d35+OOPOXnyJCEhIfznP/+hRo0afPfdd/L0cWEzvv76awYMGGCVsaWgheGcnZ3p06cPBw4c4LPPPmPVqlV4enoyYMAA9u7da3Q8IbJ1/vx5du7caZXjzyAFLUzEzs6O1q1bs3TpUvbv34+HhwcdOnQgODiYOXPmkJKSYnREITIJDQ2lWbNmODk5WWV8KWhhSk888QTDhw/nxIkTDB8+nB9//JGKFSsyfvx4i9+QRohHkZKSwrRp02jfvr3VliEFLUzN3t6eZ599lvXr17Ny5Ur2799PxYoVGThwIOvXr5e9amGYN998E3d3d1588UWrLUMKWtiMgIAA5syZQ0xMDO7u7gwfPhxXV1c6d+7M//73P5KTk42OKPKJ6dOnEx4ezk8//YSdnfVqVOWl05qCgoL0rl27jI4hctHFixdZtWoVc+fOZd++fXTr1o3evXsTFBSEUsroeCIPCgsLo1OnTmzfvh1fX9/HHk8pFa21DsrqPdmDFjbN1dWVV155hQ0bNhAdHU3ZsmXp3r07VatWZdy4cRw/ftzoiCIPiY2NpVu3bvz0008WKecHkYIWeYaXlxfDhw/n6NGjfPvtt/z555/Uq1ePwMBAPvvsMylr8VguXLhAy5YtGTZsGK1bt86VZZq+oJVSQ5VSWik11egswjYopWjYsCHTpk3jzJkzfPHFF8TFxaWX9eeff86JEyeMjilsyMWLF2nTpg09e/Zk4MCBubZcUxe0UqoO0BeIMTqLsE0FChTgmWeeyVTWJ0+epG7dugQGBjJu3DhiYmLkEnORrV27dhEUFMSzzz7LyJEjc3XZpi1opVRxYB7wGnDV4DgiD8hY1mfPnmXChAmcPXuWkJAQypcvT9++fVmyZIncZU8AoLVm9uzZtGnThkmTJjF69Ohc/+DZtGdxKKV+AU5qrT9USm0GDmit/5PFfP2AfgCenp6BcXFxuRtU2DytNUePHmXNmjWsXr2a7du3U7t2bdq0aUPbtm2pXLmynBGSz1y9epUBAwawf/9+Fi5cSNWqVa22LJs7i0Mp1RfwBYY/aF6t9QytdZDWOsjFxcX64USeo5TCz8+Pt956i9DQUM6fP8/bb7/N8ePHadOmDd7e3rzxxhusWLGCmzdvGh1XWNmmTZsICAjAxcWFXbt2WbWcH8R0e9BKqSeBcKCh1vr3tNc2k80edEZyHrSwNK01hw8fZvXq1axevZqoqCjq1atH69atadSoEdWrV6dAgQJGxxQWsH//foYNG8a+ffv45ptvaNeuXa4s93570GYs6N7AbOBOhpftAQ3cBQprrbO8ZEwKWljb9evX2bBhA6GhoYSFhXH69Gnq1KlDw4YNadiwIbVr18bZ2dnomOIhxMbGMnLkSEJDQxkyZAj9+/e32s2PsmJrBV0CKHfPy7OBo8BY4KDOJrQUtMht8fHxbNu2jbCwMMLCwjh48CABAQHphV2/fn2KFy9udExxD60127ZtY/LkyWzYsIE333yT9957j2LFiuV6Fpsq6KzIIQ5hK27cuMGOHTvSCzsqKgofH5/0sg4ODqZChQryoaNBkpOT+eWXX5g8eTLXr19n0KBBvPzyy4YU8z+koIUwSEpKCtHR0YSFhREREcHOnTu5ffs2wcHBBAcHU7NmTQICAnB3d5fStpI7d+4QHh7OggULWLx4MTVq1OCtt96iTZs2Vr3RUU7ZfEHnlBS0MDutNadPnyYyMpKdO3eyZ88e9uzZg52dHQEBAdSoUYOAgAACAgLw8/PD3t7e6Mg26ebNm2zYsIGVK1eycuVK3NzceOGFF+jatSve3t5Gx8tECloIE9Nac/bsWfbu3ZtpOnv2LNWqVSMgIAB/f3/8/Pzw9fXF09NTivseN27cYOfOnYSFhREeHk5kZCS1atXiueeeo127dlSqVMnoiNmSghbCBl2/fp2YmBj27t1LTEwMx44d4+jRo1y6dAlvb298fX2pVKlS+q+VKlWifPnyeb68k5OTOXjwILt372b37t1ERUVx+PBhAgICaNCgAQ0aNKBhw4Y28+GsFLQQeUhCQgInTpzg6NGj6aWdsbzLlStHuXLlKF++fPqvGb8uXbq06Y933717l/Pnz3PixIlM04EDB/j999/x9fWlZs2a1KhRg8DAQIKCgnL11DhLkoIWIp9ITEzk1KlTnD59OttfExIScHNzo2zZsri5uaVPxYsXp0iRIv+aihYtSuHChXF0dMTJyQlHR0ccHR0fuKd++/ZtEhMTSUxMJCEhgcTERG7evMn169e5fv06f/31F3/99Vemr8+ePcuJEyc4efIkxYoVo2LFiumTt7c3VatW5emnn85T55pLQQsh0iUkJHDhwgXOnz/P+fPn07/++++/uXHjRrZTcnIyycnJJCUlkZycjJ2dXabCdnBwIDk5Ob2Q79y5g7OzM4UKFcLZ2Tn962LFilG8eHGKFy/+r689PDzSy7hw4cJG/1HlivsVtFyjKkQ+U6hQIby9vR/rbAatNbdv304v6+TkZFJSUnB0dEwvZAcHB9MfSjE7KWghxENTSlGwYEEKFixI0aJFjY6TZxl/lrYQQogsSUELIYRJSUELIYRJSUELIYRJSUELIYRJSUELIYRJSUELIYRJ5akrCZVSlwAzPta7DBBvdAgryIvrJetkO/LKenlprbN84nWeKmizUkrtyu5STluWF9dL1sl25NX1ykgOcQghhElJQQshhElJQeeOGUYHsJK8uF6yTrYjr65XOjkGLYQQJiV70EIIYVJS0EIIYVJS0EIIYVJS0A9JKfWGUipWKZWklIpWSjW8z7xOSqkflFIxSqlbSqnN2czXOG2sJKXUCaVUf6utQDYsvV5KqSZKKZ3F9JRVVyRzhodZpyZKqWVKqXNKqYS0dXs1i/lsbVs9cL1scFtVUUptUkpdyLAdxiqlHO6Zz/Bt9di01jLlcAK6AbeAvkBlYApwA/DMZv7CwHSgH7AU2JzFPN7AzbSxKqeNfQvoZOPr1QTQQBWgbIbJ3qTrNBT4FKgPVAQGALeBHja+rXKyXra2rXyB3kB1wAtoD1wAxptpW1nkz8boALY0ATuB7+557SgwLgffOzWbIvscOHrPazOBCBtfr3/+0ZextW2VYf6FwOK8sq3us155YVtNzLgdzLCtLDHJIY4cSvvxKRBYe89ba4F6jzF03SzGDAWClFIFH2PcHLHiev1jV9qP1xuUUk0tMN4DWXCdigFXM/w+r2yre9frHza5rZRSvkBrYEuGlw3dVpYiBZ1zZQB7Un+UyugCqT8OPqqy2YxZIG2Z1mat9TpH6o/TnYCOwB/ABqVUo8cYM6cee52UUs8Czch8MYTNb6ts1ssmt5VSartSKonUve1wUg/n/MPobWUR8lTvh3fvlT0qi9csMWZWr1uTRddLa/0Hqf/Q/xGhlKoAvA9sfdRxHzbGPb/P0ToppeoD84FBWuvIHIyZ1evWZNH1suFt1Q0oSuqx6C+AD4FxDxgzq9dNSwo65+KBO/z7f3VX/v0/9cM4n82Yt4HLjzFuTllrvbKyE3jBwmNm5ZHXSSnVAFgNjNBaT7vnbZvdVg9Yr6yYfltprU+lfXlIKWUPzFRKfaG1vo3x28oi5BBHDmmtU4BooMU9b7UAtj/G0BFA8yzG3KW1vvUY4+aIFdcrKwGk/jhtVY+6Tmk/0q8BRmmtJ2Uxi01uqxysV1ZMva2yYEfqDqd92u8N3VYWY/SnlLY0kfojVQrQh9RTdyaTejqQV9r744AN93xPFVL/sv8M7Er7OiDD+/+cDjQpbcw+acvI7VO3LL1ebwPPA5WAqmljaKCjGdeJ1DMZbpL6o3LGU81cbHlb5XC9bG1b9QK6AE+ReupgV+AM8LOZtpVF/myMDmBrE/AGcBJIJvV//kYZ3vsBOHnP/CfT/rJnmu6ZpzGwO23MWKC/ra8XMBg4BiQCV4AwoK1Z1ynt9/9anyzW26a2VU7Wywa3Vfe0bfA3qUV+kNQPCJ3Ntq0ed5K72QkhhEnJMWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghhDApKWghHkAp9R+l1DGlVKJS6jellIvRmUT+IAUtxH0opcaQ+my+fkAwqTeIH29oKJFvyP2ghciGUioIiATqa60j0l4bBPxXa+1qaDiRL8getBDZex/Y+k85p7kElDEoj8hnpKCFyIJSqiDwHPDrPW85A3/lfiKRH8khDiGyoJSqRerhjSTgToa3CgJ7tNZ1DAkm8pUCRgcQwqSeJPUp0P6kPmT1HwuAbYYkEvmOFLQQWSsOxGutj/7zglKqDBAADDIslchX5Bi0EFmLB4oqpTL+GxkCRNzzoaEQViN70EJkbSOp/z6GKaXmAJ2BXkB9Q1OJfEX2oIXIgtb6EvAS0Ac4BLQAGmc85CGEtclZHEIIYVKyBy2EECYlBS2EECYlBS2EECYlBS2EECYlBS2EECYlBS2EECYlBS2EECYlBS2EECb1/ySzv5lP0UADAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Confidence region from lnL = lnLmax - Q/2 (here for theta and xi)\n", "if not(m.is_fixed('theta') | m.is_fixed('xi')):\n", " CL = 0.95\n", " Q = chi2.ppf(CL, nfreepar) # lnL = lnLmax - Q/2\n", " sig = np.sqrt(Q) # number of sigmas of contour\n", " fig, ax = plt.subplots(1,1)\n", " contour = m.mncontour('theta', 'xi', sigma=sig, numpoints=200)[2]\n", " contour.append(contour[0]) # close the contour\n", " con = np.array(contour)\n", " plt.plot(MLE[0], MLE[3], marker='o', linestyle='None', color='black', label=r'$(\\hat{\\theta}, \\hat{\\xi})$')\n", " plt.plot(con[:,0], con[:,1], color='black', linewidth=1, label=r'95%')\n", " plt.xlabel(r'$\\theta$', labelpad=10)\n", " plt.ylabel(r'$\\xi$', labelpad=10)\n", " plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2)\n", " handles, labels = ax.get_legend_handles_labels()\n", " plt.legend(handles, labels, loc='upper right', fontsize=14, frameon=False)\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 }