{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "iminuit version: 2.8.1\n" ] } ], "source": [ "# Example of maximum-likelihood fit with iminuit version 2.\n", "# pdf is a mixture of Gaussian (signal) and exponential (background),\n", "# truncated in [xMin,xMax].\n", "# G. Cowan / RHUL Physics / August 2021\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", "import iminuit\n", "from iminuit import Minuit\n", "import matplotlib.pyplot as plt\n", "from matplotlib import container\n", "plt.rcParams[\"font.size\"] = 14\n", "print(f\"iminuit version: {iminuit.__version__}\") # should be v 2.x\n", "\n", "# 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\n", " \n", "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": 2, "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))\n", "\n", "# Initialize Minuit and set up fit:\n", "par = np.array([theta, mu, sigma, xi]) # initial values (here = 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)] # set limits\n", "m = Minuit(negLogL, par, name=parname)\n", "m.errors = parstep\n", "m.fixed = parfix\n", "m.limits = parlim\n", "m.errordef = 0.5 # errors from lnL = lnLmax - 0.5\n", "\n", "# Do the fit, get errors, extract results\n", "m.migrad() # minimize -logL\n", "MLE = m.values # max-likelihood estimates\n", "sigmaMLE = m.errors # standard deviations\n", "cov = m.covariance # covariance matrix\n", "rho = m.covariance.correlation() # correlation coeffs." ] }, { "cell_type": "code", "execution_count": 3, "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 3 -0.018131 -0.531774\n", "3 0 -0.018131 -0.531774\n", "3 3 0.415591 1.000000\n" ] } ], "source": [ "print(r\"par index, name, estimate, standard deviation:\")\n", "for i in range(m.npar):\n", " if not m.fixed[i]:\n", " print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n", " \"{:.6f}\".format(MLE[i]), \" +/- \", \"{:.6f}\".format(sigmaMLE[i]))\n", "\n", "print()\n", "print(r\"free par indices, covariance, correlation coeff.:\")\n", "for i in range(m.npar):\n", " if not(m.fixed[i]):\n", " for j in range(m.npar):\n", " if not(m.fixed[j]):\n", " print(i, j, \"{:.6f}\".format(cov[i,j]), \"{:.6f}\".format(rho[i,j]))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAF7CAYAAABhB6n0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABLDElEQVR4nO3deZgU1dXH8e+BYdh3hlXZBJFFBRlQQYJRISpuwQVcEolBgrtiEon6CkkMmERNiEAUNOIu7jEoAqIgCqKDIK6IsqjIquz7ct4/qgd7lp7pGbqnemZ+n+epp+mqW7dOdU0zZ27dutfcHREREZFEqRB2ACIiIlK2KLkQERGRhFJyISIiIgml5EJEREQSSsmFiIiIJJSSCxEREUmotLADSBUNGjTwli1bhh2GiIhIiVmwYMEGd89IdL1KLiJatmxJVlZW2GGIiIiUGDNbmYx6dVtEREREEkrJhYiIiCSUkgsRERFJKCUXIiIiklBKLkRERCShlFyIiIhIQim5EBERkYRSciEiIiIJpeRCREREEkrJhYiIiCSUkgsRERFJKCUXIiIiklBKLkRERCShlFyIiIhIQim5EBERkYRSciEikiTr169n5MiRrF+/PuxQREqUkgsRkSS56qqryMrK4pprrgk7FJESpeRCRCQJnnzySSpXrsyUKVOoVKkSzzzzTNghiZQYc/ewY0gJmZmZnpWVFXYYIiIiJcbMFrh7ZqLrVcuFiIiIJJSSCxEREUmoUJMLM7vazJab2S4zW2BmvQooW8XMJpnZYjPba2azYpRLN7M/RerdbWZfm9n1STsJEZEYxo4dS5s2bahatSqnn366nhqRciO05MLMBgBjgFFAF2AuMNXMmsfYpSKwCxgLvFJA1U8BpwNDgHbAhcDiBIUtIhKX2267jbvvvpsJEyYwf/58li1bxu9///uEHmP8+PG0atWKKlWq0LVrV+bMmVNg+dGjR9OtWzdq1apFRkYGZ599Nh9//PEh1Ttq1CjMjGuvvTbPttWrV3P55ZeTkZFBlSpV6NChA7Nnzy76iRYjrnj3GzlyJGaWY2ncuHGO/eP53Fq2bJmnHjOjX79+xT7fUs3dQ1mA+cDEXOuWAqPj2HcsMCuf9X2BzUCDosbTtWtXFxFJhPfff9/NzOfOnXtw3ZgxYzwjIyNhx3j66ac9LS3NJ0yY4J9++qlfe+21Xr16dV+5cmXMffr27ev/+c9//KOPPvLFixf7eeed540aNfLvv/++WPXOmzfPW7Zs6cccc4xfc801ObZt3LjRW7Vq5b/4xS98/vz5vmzZMn/99df9008/zTe2yy+/3EeMGJHQ841nvxEjRni7du189erVB5d169YV+XNbt25djjo++OADNzOfNGlSgfGFDcjyZPyOT0alhR4U0oF9wIW51o8DZsexf6zkYjzwOkFryLeRZOVfQI3C6lRyISKJMmDAAO/du3eOdU8++aSbWcKO0b17dx88eHCOdW3atPHhw4fHXcfWrVu9QoUK/vLLLxe53k2bNnnr1q195syZ3rt37zzJxR/+8Afv0aNH3LEUllwU93wL22/EiBHesWPHuON0z/9zy+3OO+/02rVr+/bt24tUd0lLVnIR1m2RBgS3OdbmWr8WaJy3eNxaAycBxwLnA9cS3CKZlF9hMxtiZllmlqV7oSKSCHv37uV///sf/fv3z7F+586d1K5dO0/5UaNGUaNGjQKX3M3/e/bsYcGCBfTt2zfH+r59+zJ37ty4Y926dSsHDhygbt26Ra53yJAhXHDBBZxyyin51v3SSy9x/PHHM2DAABo2bEjnzp0ZO3Zs9h+CRVLc8413v2XLltGsWTNatWrFwIEDWbZsWYHx5P7ccnN3HnroIS677DKqVatW2OmVSWkhHz/3T5nls64oKkT2v8TdNwOY2bXANDNr5O45khl3nwBMAGjePtPdwewQji4i5d6iRYvYsWMHt9xyC7feeuvB9Xv37qVLly55yg8dOpSLLrqowDqbNWuW4/2GDRvYv38/jRo1yrG+UaNGvP7663HHesMNN9C5c2dOPPHEItU7ceJEvvzySx577LGYdS9btozx48dz0003MXz4cBYtWsR1110HwLXXXsuoUaMYNWrUwfK7d+/GzLj77rsPrps6dSq9evUq9vnGs9/xxx/PpEmTOOqoo1i3bh133nknPXr04JNPPqF+/fpxfW65zZgxg+XLlzN48OCYsZV1YSUXG4D95G2laEje1oyiWA2syk4sIj6LvDYvqO612+HDtdD5UNpNRKTcW7JkCenp6SxevBiL+mvl4osvpmfPnnnK16tXj3r16hXrWJbrryF3z7MulmHDhvH222/z9ttvU7FixbjrXbJkCbfeeitz5swhPT09Zv0HDhwgMzOT0aNHA9ClSxeWLl3KuHHjuPbaa/MkVbfccgvNmjXj+ut/fLgvd1JV3PMtaL8zzjgjx7YTTjiB1q1b88gjjzBs2LA8dRX0uWWbOHEi3bp1o3PnzoXGVlaFclvE3fcAC4A+uTb1IXhqpLjeAZqaWY2odUdGXlcWtKMBT+XtNC0iUiSbN2+mQYMGtG3bljZt2tCmTRvq1KnDokWLuOCCC/KUL85tkQYNGlCxYkXWrFmTY/26devy/JWen5tuuomnnnqKN954g9atWxep3nnz5rFhwwY6depEWloaaWlpzJ49m/Hjx5OWlsbu3bsBaNKkCR06dMhRT/v27fn666+BIKnK/nzatGlDzZo186yrWrXqIZ1vcfarUaMGHTt2ZOnSpXF/brnr/u9//8uVV14ZM67yIMxxLu4FBpnZYDNrb2ZjgKbA/QBmNtrMZkbvYGYdzKwzQZ+NGmbWOfI+25PA98DDZtbRzHoSPO76nLuvKyiYOlXg5S9g255EnZ6IlEcNGjQ4eE8+2+jRoznxxBPzbUYfOnQoixYtKnDJzMw5OnN6ejpdu3ZlxowZOdbPmDGDHj16FBjfDTfcwJNPPskbb7zBUUcdVeR6zzvvPD766KM88Q0cOJBFixYdbM3o2bMnS5YsyVHPF198QYsWLQqMLz/FPd/i7Ldr1y4+//xzmjRpkmN9QZ9btEmTJlG5cmUGDhxY2GmVbcnoJRrvAlwNrAB2E7Rk/CRq2yRgRa7yKwj6VORYcpVpB0wHdgCrCJ5AqVlYLO2P6erN/+n+1EdF7WsrIvKjdevWedWqVf1Pf/qTL1++3P/+9797RkaGf/HFFwk9ztNPP+2VKlXyiRMn+qeffurXX3+9V69e3VesWHGwzH333eft2rU7+P7qq6/2mjVr+syZM3M8Nrl169Yi1Ztbfk+LvPfee56WluZ33nmnL1261J955hmvVauWjx071t2DJy6iY8hv2b179yGdbzz73XzzzT5r1ixftmyZv/vuu96vXz+vWbNmjnrj+dzc3Q8cOOBt27bN83RKKqMsPYqaikvXrl391Efdz326CFdFRCQfzz77rDdv3tyrVq3qffv2jTm2w6EaN26ct2jRwtPT0/24447z2bNn59g+YsSI7D/A3N3z/GGWveR+BLSwenPLL7lwd58yZYofc8wxXrlyZW/btq2PGTPGDxw4kCO2gpY333zzkM43nv0GDBjgTZo08UqVKnnTpk29f//+/sknn+TYP97P7Y033nDA58+fX+DnlUqSlVxoVtSIzMxMHzohiz/PgWmXwlENwo5IREQkuTQragno3x7SK8LTn4QdiYiISOml5CJKvarwsyPghc9g176woxERESmdlFzkMrAjbN4N074KOxIREZHSSclFLj0Oh8NracwLERGR4lJykUsFgwEdYd63sGJT2NGIiIiUPkou8nFhhyDJmKyOnSIiIkWm5CIfjWvAaa3gmU9gz/6woxERESldlFzEcOnRsGGnOnaKiIgUlZKLGH7SIujY+fjisCMREREpXZRcxFDBgtaLd1fB0h/CjkZERKT0UHJRgIs6QKUK8MRHYUciIiJSeii5KED9anBGG3j+U9i5N+xoRERESgclF4X4xTGwZQ/874uwIxERESkdlFwUoltTOLI+PK5bIyIiInFRclEIi3Ts/HAtfLQ27GhERERSn5KLOPQ/CqqmqfVCREQkHkou4lCrMpzbDv67JJgxVUQkHuvXr2fkyJGsX78+7FBESpSSizhddgzs3AfPfRp2JCJSWlx11VVkZWVxzTXXhB2KSIlSchGnoxtC1ybw6IdwwMOORkRS3ZNPPknlypWZMmUKlSpV4plnngk7JJESY+76TQmQmZnpWVlZBZZ5eQlc9xpMOhd+2rJk4hIREUkWM1vg7pmJrlctF0VwehvIqAaPfBh2JCIiIqlLyUURpFeEy46GWStgxaawoxEREUlNSi6K6JKjIa1C0PdCRCSWK664AjPLsfTo0SPssERKhJKLImpYPZhv5NlPYfuesKMRkVQ1aNAgWrRowUUXXcSzzz7LwoULefnllxNW/8iRI/MkL40bNy5wn7feeotzzjmHZs2aYWZMmjQp33Ljx4+nVatWVKlSha5duzJnzpwc21u2bJnn2GZGv379DpbZunUrN954Iy1atKBq1ar06NGD999//5DOubC4Ylm9ejWXX345GRkZVKlShQ4dOjB79ux8y44aNQoz49prry3y9qIcp6xTclEMgzoH8428+HnYkYhIKtq7dy+XXXYZ119/PZMnT+aCCy6gc+fONGjQIKHHadeuHatXrz64fPRRwSP9bdu2jU6dOjFmzBiqVq2ab5nJkydzww03cOutt7Jw4UJ69OjBGWecwddff32wzPvvv5/juB988AFmxkUXXXSwzODBg5k2bRqPPPIIH330EX379uW0005j1apVMeMbNGgQI0eOLHZc+dm0aRM9e/bE3XnllVf47LPPuO+++2jYsGGesu+++y4TJ07kmGOOybeugrYX5Tjlgrtrcadr164erwMH3Ps96X7aY8G/RUSiLViwwCtUqOA7d+5M2jFGjBjhHTt2LPb+1atX94cffjjP+u7du/vgwYNzrGvTpo0PHz48Zl133nmn165d27dv3+7u7jt27PCKFSv6Sy+9lKPccccd57fddlvMei6//HIfMWJEvtuKE5e7+x/+8Afv0aNHgWXc3Tdt2uStW7f2mTNneu/evf2aa64p0vZ4j5NqgCxPwu/UUFsuzOxqM1tuZrvMbIGZ9SqgbBUzm2Rmi81sr5nNKqTuk8xsn5l9nPi44fJj4YvvYd63ia5dREq7OnXqcODAAf7yl7+wcuVKDhw4kG+5UaNGUaNGjQKXgpr+ly1bRrNmzWjVqhUDBw5k2bJlhxT3nj17WLBgAX379s2xvm/fvsydOzfffdydhx56iMsuu4xq1aoBsG/fPvbv30+VKlVylK1atSpvv/12icSV7aWXXuL4449nwIABNGzYkM6dOzN27Fg81zAMQ4YM4YILLuCUU07Jt57Ctsd7nPIitOTCzAYAY4BRQBdgLjDVzJrH2KUisAsYC7xSSN11gUeBmQkLOJezj4R6VeHhRck6goiUVq1bt+aBBx7g7rvvpmXLlqSlpbFw4cI85YYOHcqiRYsKXDIz8x+C4Pjjj2fSpElMnTqViRMnsmbNGnr06MH3339f7Lg3bNjA/v37adSoUY71jRo1Ys2aNfnuM2PGDJYvX87gwYMPrqtZsyYnnngid955J6tWrWL//v08/vjjzJs3j9WrVx8slzu5euKJJ/KsmzNnTrHiyrZs2TLGjx9P69atmTZtGjfccAPDhw9n3LhxB8tMnDiRL7/8kj//+c/51lHY9niPU56khXjsYcAkd58YeX+dmZ0OXAX8IXdhd98ODAUws2OAOgXU/RDwCGDABQmM+aAqaXBJJxj3fvBYasuCohGRcuWee+7hrrvu4qabbuLkk0+mUaNGdOrUKU+5evXqUa9evWId44wzzsjx/oQTTqB169Y88sgjDBs2rFh1ZjOzHO/dPc+6bBMnTqRbt2507tw5x/rHHnuMK664gsMOO4yKFSty3HHHcfHFF/PBBx8cLDN06NAc/TRuueUWmjVrxvXXX39wXbNmzdi4cWOR48p24MABMjMzGT16NABdunRh6dKljBs3jmuvvZYlS5Zw6623MmfOHNLT0/PsX9j2eI9T3oTScmFm6UBXYHquTdOBQ3pWy8yuBhoDdx5KPfH45bHBY6lqvRCRbO+++y7Dhw9n1qxZjBo1ir59+3LsscdSsWLFPGUP9bZItBo1atCxY0eWLl1a7NgbNGhAxYoV87QGrFu3Lk+rQfb6//73v1x55ZV5th1xxBHMnj2bbdu28c033/Dee++xd+9eWrVqdbBMvXr1aNOmzcGlZs2aedZVrVq1yHFFa9KkCR06dMixrn379gc7gs6bN48NGzbQqVMn0tLSSEtLY/bs2YwfP560tDRmzZpV4Pbdu3fHdZzyJqyWiwYEtznW5lq/FjituJWa2dHACOAEd99fWEZrZkOAIQDNm8e6GxNbo+rB7ZFnPoVhJ0LtysUIWkTKlNdee40OHTrQsWPHQsvm/ss9P82aNYvruLt27eLzzz/npz/9aVzl85Oenk7Xrl2ZMWMGF1544cH1M2bM4Pzzz89TftKkSVSuXJmBAwfGrLN69epUr16djRs3Mm3aNP72t78lPa5oPXv2ZMmSJTnWffHFF7Ro0QKA8847L8+tp1/96le0bduWW2+9lWbNmtGzZ8+Y27NbMwo7TrmTjF6ihS1AU8CBXrnWjwA+j2P/scCsXOsqA58Av4haNxL4OJ6YivK0SLSP1ro3/6f7/VnF2l1EypgHH3zQK1So4HfddZd//PHH/sMPPyTlODfffLPPmjXLly1b5u+++67369fPa9as6StWrHB39/vuu8/btWuXY5+tW7f6woULfeHChV61alX/4x//6AsXLvSVK1ceLPP00097pUqVfOLEif7pp5/69ddf79WrVz9Yb7YDBw5427Zt8zzBke21117zV1991ZctW+bTp0/3Y4891rt37+579uzJEc/q1asLXHbv3h13XPmd83vvvedpaWl+5513+tKlS/2ZZ57xWrVq+dixY2N+tvk9DVLY9uIcJxWQpKdFwkou0oF9wIW51o8DZsexf37JRctIwrIvajkQta5vQXUWN7lwd7/oWfcTH3Lfu7/YVYhIGbF//34fNWqUd+jQwatUqeKA9+/fP+HHGTBggDdp0sQrVarkTZs29f79+/snn3xycPuIESM8+PvxR2+++aZH/k/MsVx++eU5yo0bN85btGjh6enpftxxx/ns2bPzHP+NN95wwOfPn59vfJMnT/bWrVt7enq6N27c2K+55hrftGlTjjLZMRa0vPnmm3HHld85u7tPmTLFjznmGK9cubK3bdvWx4wZ4wcKGEegOMlFcY6TCpKVXIQ2K6qZzQc+dPchUeu+AJ539zwdOnPtOxbo5O4nR62rBLTLVfRqoA/wc2CFu2+LVWc8s6LGMv0ruHIKjDsDzjqyWFWISBk1Y8YM+vbty+bNm6lVq1bY4YjkkKxZUcN8WuRe4DEzew94h+BJkKbA/QBmNhro7u6nZu9gZh0IWj0aADXMrDOAuy9y971AjjEtzGwdsNvdEz7WRbRTW0GL2vDgQiUXIvKjHTt28M4773DUUUcpsZByJbTkwt0nm1l94HagCUFicKa7r4wUaQIckWu3V4Ho3jHZD44X3HMzySpWgF91hpGzYcFq6NokzGhEJFU888wzvPrqq0yePDnsUERKVGi3RVLNodwWgWASsxMegl4tYPyZCQxMREQkSZJ1W0QTlyVI9XQY2AmmfgnfbAk7GhERkfAouUigX3WGCgYP5R3lV0REpNxQcpFATWvCue3g6Y9h486woxEREQmHkosEG3Ic7NwHj30UdiQiIiLhUHKRYEc1gJNbwKRFsGtf2NGIiIiUPCUXSTA0E77fCc9/FnYkIiIiJU/JRRKc0AyOaQgTFsD+A2FHIyIiUrKUXCSBWdB6sWIzTF8WdjQiIiIlS8lFkpx+BDSvDfdngcYpExGR8kTJRZJUrABXdoFFa+G978KORkREpOQouUiiCztAvaow/v2wIxERESk5Si6SqGoluKIzzFoJn6wPOxoREZGSoeQiyX55LNRMh3FqvRARkXJCyUWS1a4MvzwGXl0KX20MOxoREZHkU3JRAq7oAukV4d/Fn9FdRESk1FByUQIaVIOLO8GLn8O3mo5dRETKOCUXJeQ3XYPXCR+EG4eIiEiyKbkoIU1rQv+jgunY128POxoREZHkUXJRgq7OhL0H4KGFYUciIiKSPEouSlCrutCvLTz2EWzeFXY0IiIiyaHkooRdkwnb9qj1QkREyi4lFyWsfUYwqdnDi2Dz7rCjERERSTwlFyG44XjYsgf+o9YLEREpg5RchKBDBvzsiCC5UOuFiIiUNUouQnJ996D1YtKisCMRERFJLCUXIenUEPq2hgcXwha1XoiISBkSanJhZleb2XIz22VmC8ysVwFlq5jZJDNbbGZ7zWxWPmX6m9l0M1tvZlvNbL6ZnZPUkzgE1x8fJBYPLwo7EhERkcQJLbkwswHAGGAU0AWYC0w1s+YxdqkI7ALGAq/EKNMbeAPoF6nzVeDFgpKWMB3dEE5rFTyWulWtFyIiUkaE2XIxDJjk7hPd/TN3vw5YDVyVX2F33+7uQ919AvBtjDI3uPtd7v6eu3/p7n8EFgDnJekcDtmNxwedOh/+MOxIREREEiOU5MLM0oGuwPRcm6YDPRJ8uJrAxgTXmTBHN4I+rWHiAo3aKSIiZUNYLRcNCG5zrM21fi3QOFEHMbNrgMOAx2JsH2JmWWaWtX79+kQdtsiGnRA8OaIZU0VEpCwI+2kRz/Xe8llXLGZ2PvB34FJ3X5nvwd0nuHumu2dmZGQk4rDF0iEDzmoL/1kE3+8ILQwREZGECCu52ADsJ28rRUPytmYUWSSxeAz4pbu/fKj1lYSbToBd+2B8VtiRiIiIHJpQkgt330PQ0bJPrk19CJ4aKTYzuwh4HBjk7s8dSl0lqU096H8UPLYY1mwLOxoREZHiC/O2yL3AIDMbbGbtzWwM0BS4H8DMRpvZzOgdzKyDmXUm6LNRw8w6R95nbx8IPAEMB94ys8aRpV7JnNKhueF42O9w33thRyIiIlJ8aWEd2N0nm1l94HagCfAxcGZU/4gmwBG5dnsVaBH1PnvqL4u8DiU4p39GlmyzgZMTFHrSNK8NAzvC5E/gN12D9yIiIqWNuSek/2Spl5mZ6VlZ4Xd4WL0Vej8C5xwJd/cNOxoRESnLzGyBu2cmut6wnxaRXJrUhMuOgec/h6U/hB2NiIhI0Sm5SEHXZEK1SvD3Q+raKiIiEg4lFymofjUYchxM+woWrA47GhERkaJRcpGiBneBjGpw19ugbjEiIlKaKLlIUdXTg0dT3/sO3lgRdjQiIiLxU3KRwgZ2hJa14a/vwP4DYUcjIiISHyUXKaxSRfhdD1jyPbz4edjRiIiIxEfJRYo7sy0c0xDufTeYe0RERCTVKblIcRUMhp8Eq7bCo4vDjkZERKRwSi5KgZ6HQ+8WMPY92LQr7GhEREQKpuSilLjtJNi6B8bMDzsSERGRgim5KCXaNQieHnl0MSzbGHY0IiIisSm5KEWGnQCVK8Jd74QdiYiISGxKLkqRjOpwdWYwLPi734YdjYiISP6UXJQyg4+DpjXgz3PggIYFFxGRFKTkopSpkga/7wkfr4OXNLCWiIikICUXpdC57YKBtf46F3buDTsaERGRnJRclEIVDO74CazZBv/OCjsaERGRnJRclFLdmsE5R8L9C+CbLWFHIyIi8iMlF6XYH04KWjH+MifsSERERH6k5KIUa1oTru4GU7+Eud+EHY2IiEhAyUUpN+Q4OKwWjJwN+w6EHY2IiIiSi1KvShr8Xy9Y8j08/lHY0YiIiCi5KBN+dkQwc+q982DjzrCjERGR8k7JRRlgBiN7w7Y98Pe5YUcjIiLlnZKLMuLI+jDoWHjyY/hwTdjRiIhIeRZqcmFmV5vZcjPbZWYLzKxXAWWrmNkkM1tsZnvNbFaMcr0jde0ys2VmNjRpJ5BibjohmNzstjdhvzp3iohISEJLLsxsADAGGAV0AeYCU82seYxdKgK7gLHAKzHqbAW8GqmrCzAauM/Mzk9s9KmpZmW4vRd8tC5owRAREQlDmC0Xw4BJ7j7R3T9z9+uA1cBV+RV29+3uPtTdJwCxJhwfCnzn7tdF6pwIPAL8NhknkIrOORJ6HAZ/mwsbdoQdjYiIlEehJBdmlg50Babn2jQd6HEIVZ+YT53TgEwzq3QI9ZYaZvDnnwYTmo1+O+xoRESkPAqr5aIBwW2OtbnWrwUaH0K9jWPUmRY5Zg5mNsTMsswsa/369Ydw2NTSpl4wuNZzn8F7q8KORkREypuwnxbxXO8tn3WJqDO/9bj7BHfPdPfMjIyMQzxsarm2OzSrCbe/CXv3hx2NiIiUJ2ElFxuA/eRtpWhI3paHolgTo859wPeHUG+pU61SMPbFku/hoYVhRyMiIuVJKMmFu+8BFgB9cm3qQ/CkR3HNA07Lp84sd997CPWWSn2PCEbv/Md8WLkp7GhERKS8CPO2yL3AIDMbbGbtzWwM0BS4H8DMRpvZzOgdzKyDmXUm6D9Rw8w6R95nux84zMz+GalzMDAIuDv5p5Oa/nQypFUIxr7wQ73hJCIiEoe0sA7s7pPNrD5wO9AE+Bg4091XRoo0AY7ItdurQIuo99kN/hapc7mZnQn8g+CR1u+A6939+eScReprXANu6QH/NwteXAL9jwo7IhERKevMi/jnrJlVJmhhqAqsd/cy8ZhFZmamZ2VlhR1GUhxwOP9ZWLEJZv4C6lUNOyIREUkFZrbA3TMTXW9ct0XMrKaZXWVmbwGbgS8JWhrWmNk3ZjbRzLolOjhJjAoGd50CW3bDX+aEHY2IiJR1hSYXZnYTsAK4ApgBnAt0Bo4kGLRqBMHtlRlm9pqZtU1WsFJ87RrA0K7B2Bdvfx12NCIiUpYVelvEzJ4F/ujuBc5WEbld8mtgj7s/mLgQS0ZZvi2Sbdc+OP0J2O8w7dLgcVURESm/Qrst4u4XFpZYRMrtdvfxpTGxKC+qpMFfT4OvNwdzj4iIiCRDPLdF6pvZQ2a2xsz2mdn3ZjbXzP5uZt1LIkhJnOObweXHwqRF8L6GBhcRkSSIp0Pn48BJwF+Ay4DfAccBvYG3zWyWmbVOXoiSaLf0gGa14HevB7dKREREEime5KI3cL673+fuT7v7f4C9wECCR1IXAHPN7MgkxikJVD0d/noqLN8E98wLOxoRESlr4kkuVgH189vg7hvc/WaCETDHJDIwSa6TmsMlneDBhbBwTdjRiIhIWRJPcjEGeLiQ/hXPAr0SE5KUlFtPgsbV4bczdHtEREQSJ56nRcYCTwHzzGymmV2bz36/AMrESJ3lSc3KcNdp8OUPcO+7YUcjIiJlRVwjdLr7bcDxBAnEXQRDf38SGZ1zEzAcuDlZQUry9G4Bl3aCCQvgPT09IiIiCRD3xGXungUMNLM0ghE62wG1gQ3AG+6+ISkRStLd1gvmfAPDpsNrl0KN9LAjEhGR0qzIU667+z53z3L3JyKDZj2jxKJ0q54O9/aFb7fAnZp7REREDlGRkwspm7o1DeYeeepjmLk87GhERKQ0U3IhB910ArRvALe8Dj/sDDsaEREprZRcyEGV0+AffWHTLrj1DShkTjsREZF8FTm5MLPmZtYkxnolK6Vc+wz4bQ+Y+iVM/iTsaEREpDQqTjKwApgZY/2HZvaTQwlIwjfkOOh5OIycDV9tDDsaEREpbYqTXFwB3Bpj/QvA3w8pIgldBQtuj1RJg+tfgz37w45IRERKkyIlF2aW4e6T3P2l3Nsi60e4+/EJi05C06gG/O00+Hgd3K3JzUREpAiK2nIxV9Orlx99j4DLjoYHFsDbX4cdjYiIlBZFTS5eJUgwjoteaWY/MbN3EheWpIrbe0HbenDjNPh+R9jRiIhIaVCk5MLdbyCYXv1NM+trZp3N7DXgTUB/25ZBVSvBfafDlt1w03Q4oMdTRUSkEMUZ/vtuYBQwBXgf2Aoc4+4XJzg2SRHtM2Bkb5i9Ev6dFXY0IiKS6oraofNwM3sA+BNBYrEbeMXdNSJCGXdxJzjnyKBzp2ZPFRGRghS15WIp0AU4y917AucA/zCz2xIemaQUMxh9KrSoDddOhQ3qfyEiIjEUNbm4zN27u/sMAHd/AzgZuMrMxic6OEktNdJh/JnB8OA3TVP/CxERyV9RO3Q+l8+6D4GeBElGkZjZ1Wa23Mx2mdkCM+tVSPmjzWy2me00s1VmdoeZWa4yl5jZIjPbYWZrzOxxM2tc1Ngkfx0i/S/e+hrGvR92NCIikooSMheIu68kSDDiZmYDgDEEnUO7AHOBqWbWPEb5WsAMYC3QDbge+B0wLKpMT+Ax4BGgI3Ae0AF4okgnJAW6uBOc2w7ufRfeWhl2NCIikmoKTS7MrFU8Fbn7RgscHuexhwGT3H2iu3/m7tcBq4GrYpS/FKgGXO7uH7v788BfgWFRrRcnAt+6+z/cfbm7vwvcB2jU0AQyg7tODca/uO41+GZL2BGJiEgqiaflYp6ZPWRmJ8YqYGZ1zewq4FPg3MIqNLN0oCswPdem6UCPGLudCMxx951R66YBTYGWkffvAE3M7OxIotMAGEgw+JckULVKMOEsOHAAhr4Cu/aFHZGIiKSKeJKLl4DNwCtmtt7MXjOzh83s32b2tJktBtYBlwE3uvvYOOpsAFQkuMURbS0Qq39E4xjls7fh7vOAiwlug+wB1gMGXJ5fhWY2xMyyzCxr/fr1cYQt0VrWgX/+LJh/5PY3wNXBU0REiC+5uILg9sNhQH3gO6AO0ArYR9C/oYu793T3aUU8fu5fR5bPusLKH1xvZh2AfwF/JmgZOZ0g8Xgg38rcJ7h7prtnZmRkFDF0ATi1NdzQHZ79DJ74KOxoREQkFaTFUeYb4Hh3fznStWG4u687xONuAPaTt5WiIXlbJ7KtiVGeqH3+ALzn7tnTvi82s+3AHDO7zd2/ObSwJT83ngAfroWRs+GoBpDZNOyIREQkTPG0XNwFPG9mHxC0EFxhZr0iT28Ui7vvARYAfXJt6kPw1Eh+5gG9zKxKrvLfASsi76sRJC3Rst8bkhQVDP51OjSrCb+ZAqvUwVNEpFwrNLlw94lAJ+Bpgl/Qg4CZwEYzW2ZmL0TGmziniMe+FxhkZoPNrL2ZjSHonHk/gJmNNrOZUeWfBHYAk8ysk5n1B4YD97ofvNv/P+BcM7vKzFpHHk39F/CBu2titSSqXQUeOgd274fB/4Mde8OOSEREwhLXOBfuvsTd/0Yw/PdJQE2gO/AXYBVBC8KjRTmwu08GbgRuBxZF6j0zMmYGQBPgiKjymyPHaQpkAeOAewiSlOwykwgecb0W+Bh4LhJzoU+wyKFrUw/uOwM+/x6GaQZVEZFyy1xd/AHIzMz0rCxN+ZkIEz6Av8yBG4+Hm04IOxoREYnFzBa4e2ai642nQ6dIkVzZBb7YAP+cD0fWh35tw45IRERKUkKG/xaJZgZ/OQW6NgkmOPtgddgRiYhISVJyIUlROQ0mngWNa8Cv/wcrN4UdkYiIlBQlF5I09avBpHODjp2D/gsbdxa+j4iIlH5KLiSpWteFB8+GVVvhyimag0REpDxQciFJ160p3NMX3v8OfjtDj6iKiJR1elpESsTZRwYjd45+BxpWh//rFXT8FBGRskfJhZSY33SFNdvhoYVQvypc0y3siEREJBmUXEiJMYM7fhJ07PzbXKhbBS45OuyoREQk0ZRcSImqYHB3H9i0C257E+pUgTM1yJaISJmiDp1S4ipVhPv7QZfGcMM0eFtTyomIlClKLiQUVSvBw+dAqzrBLKrvrQo7IhERSRQlFxKa2lXg8Z9D05rBIFsLNEy4iEiZoORCQtWwOjzVHzKqweUvwYdrwo5IREQOlZILCV2jGvDU+UHnzstego/XhR2RiIgcCiUXkhKa1gwSjJrpcOmL8Mn6sCMSEZHiUnIhKePwWsEtkmppMPB5TdUuIlJaKbmQlNKiDjx7YTDA1mUvwrxvw45IRESKSsmFpJzDagUJRtOaQSfPN1eEHZGIiBSFkgtJSY2qw+TzoU09uPJ/8OrSsCMSEZF4KbmQlFW/WtDJ8+hGcM1UePTDsCMSEZF4aG4RSWm1K8OTP4drp8L/zYLvtsLvewZzlEh43GHNNli8Lnh0eM02+GEn/LALftgBO/dBzcrB0z+1KgfLEXXh2EbBUr9a2GcgIsmk5EJSXtVK8MBZMGIW/HsBrN4Gf+8D6RXDjqx8+XozTP8K5n4Li9fC+h3B+ooGGdWhXhWoWxWOaQRV0mDbHti6J5ikbsUmmPIFeKSuw2tBZlM4tx30ag5pakMVKVOUXEipkFYB7vxp0Mnzb3Nh3Xa4/6ygZUOSwz0Yb+S1r4KkYsn3wfoj6kLvFsHtqmMaQoeMIJkozLY9QSvHh2uD5c0V8OLnweis57aD89sHdYlI6WfuXnipciAzM9OzsrLCDkPi8MJn8PvXg6dKJp4NbeuFHVHZsnk3vPQ5PP0xfLohuAXVrSn0bQ19j4DmtRNznD37gwTj+c/gjeWw9wCc0Ax+2yM4nogkn5ktcPfMhNer5CKg5KJ0eX8VDH0Fdu2HMT+D01qHHVHpt2gNPLoYXlkKu/YFrQgXd4SzjoR6VZN77I074YXP4d9Zwe2WU1oGSUZHtWSIJFWykotQ73Sa2dVmttzMdpnZAjPrVUj5o81stpntNLNVZnaHmVmuMulm9qdIvbvN7Gszuz65ZyIlrVsz+N/FP07Z/q/3gmZ8KZoDDq8vg4ueg3Mnw7Sv4PyjYMpAmHoJ/PLY5CcWEPTV+HUXeGsQDO8ZzJB75pNw3VRYvz35xxeRxAqtz4WZDQDGAFcDb0dep5pZB3f/Op/ytYAZwFtAN6AdMAnYDtwTVfQp4HBgCLAUaASUwH+PUtKa1oTnLoThM+GeecH9/L+fFkzlLgXbsz+4vTThA/hqIzSrCf/XCwZ2ghrp4cVVrRJclQmXHA0TFsDED+Dtb4L+Nv3ahheXiBRNaLdFzGw+sNjdr4xatxR4zt3/kE/5q4C/Ao3cfWdk3e3AVcBh7u5m1hd4FjjC3TcUJR7dFim93OE/i2DU29CwGvzr9KBlQ/LatQ8mfwL3Z8F324LbDr/pCme2gUop+PTNF9/DzdODR17PORL+dHLQyiEiiVGmbouYWTrQFZiea9N0oEeM3U4E5mQnFhHTgKZAy8j784D3gWFm9q2ZLTWzf5lZjUTFLqnHLGhSf+HC4BfkRc/DmPmw/0DYkaWOnXvhwQ+g1yS4Y1bQ6vPoefDKxcGTGqmYWAAcWR9euAhuPhGmfgl9Hoe3VoYdlYgUJqw+Fw2AisDaXOvXAo1j7NM4RvnsbQCtgZOAY4HzgWuB0wlun+RhZkPMLMvMstav1xzfpd2xjYNflmcfCfe+C5e8CN9uCTuqcO3cG9xaOOlh+POc4DHSp/oHt5N6twgSs1RXqSJc3x1eHhj0/7j8v8E5qY+NSOoKe+ia3P89WD7rCisfvb5C5N+XuPt8d59GkGCcb2aN8lTmPsHdM909MyND3dLLgpqVg6dH7ukTDPTU53F4eFH5a8WITirunAPtGsCzF8DT50OPw0tHUpFbhwx48aLgkdg758BvZwS3eUQk9YTVoXMDsJ+8rRQNyds6kW1NjPJE7bMaWOXum6PKfBZ5bV5A3VKGmMEFHeCEw+DWN2DkbHj5C/jrqUEze1m2fQ88/lGQWKzfAT0Ph38fD93LSB+U6unw737wr/nwj/lBZ9QHzgomuhOR1BFKy4W77wEWAH1ybeoDzI2x2zygl5lVyVX+O2BF5P07QNNcfSyOjLzqTm05c1gteORc+EdfWLYxeLTxnnmwY2/YkSXe1t0w9j3o+XDQsbVd/aCl4sn+ZSexyFbB4MYT4N9nwucb4Jyng46fIpI6wnxaZADwGMEjqO8AQ4FfAx3dfaWZjQa6u/upkfK1gSXALOBOgqRhEvBHd78nUqYGQUvFu8BIoA7wAPCZu19YUDx6WqRs27AD/hhpwWhYHW4+AS7sABXDvjF4iNZvh0cWwyMfwpbd8NOWQf+E45qEHVnJ+HQ9/PIl2O9BB9WjGxa2h4hEK5MjdJrZ1cDvgSbAx8BN7v5WZNsk4GR3bxlV/mhgHNAd2AjcD/zJo07CzNoB9xF07NwIvAQMd/etBcWi5KJ8yPouuF+/cA0cVR9u7QU/aV76+iAs/SF4+uPFz4MxK/oeAdd1C+b7KG+Wb4RLXwySq4fP1dDhIkVRJpOLVKLkovxwh1e/hLveCWb67NIYhnYNfkGn8lTu+w/A7JVBn4qZy6FyxaD15dddoHXdsKML16otQYKxZhs8eDac1DzsiERKByUXSabkovzZvQ+e/iTo/PjNFmhdB4Z0hZ8fFd8snyVl7bYgzsmfwKqt0KAq/OKYYKlfLezoUsf67XDZi7BsE9x/Jpyq+WZECqXkIsmUXJRf+w4EAzTdvyAYQrxOFTirLfRvD8c1DueWyaZdwTwfryyFt78O+hScdHgwLHaf1pCeooNehW3TLvjFi8H08A+fGzwtIyKxKblIMiUX4g7vfBO0EEz7Cnbvh5a14byjglk6OzVMXgdQ9+AWzTvfwGtfBa/7DgRPvJxzJAzoCC3rJOfYZc3GnTDg+aA16omfl5/OrSLFoeQiyZRcSLStu4Nf8i98BvO+DUZmq1UZTjwsWDKbBqNdVqtUvPp37g2a7xevhfmr4N1vYfW2YFvz2tCvTTBRV6eGpa+zaSpYux0uehZ+2AWTzw8G4BKRvJRcJJmSC4ll/fYgwXjnG5j7bdDCkK1pjaAzZau6ULdKMMhT9UrBawVg255g2bIHNu8K9v1qY9B3IluDqsGAXyccBsc3g7b1lFAkwrdb4IJnYe9+ePZCdXoVyY+SiyRTciHx+nozfLQuGJjrq43B6/JNwaOQsVS0YCrz5rWDFo8j6gWv7eoHr0omkuOrjXDhs8GTNS8OgMaawlAkh2QlFynUJ16kdGheO1hyO+DB7Y7te4NRQPd7kFDUSg+ePlECUfKOqAuP/Rwueg5+9V945oJg/hkRSa5SPj6hSOqoYMHtkIbVg86XR9QN5ryoWkmJRZg6ZsD4M4MnSK6ZGtwmEZHkUnIhImVe7xYw6pRgELLb39R07SLJptsiIlIuDOwUPJ469v3gttY13cKOSKTsUnIhIuXGb08MEoy/zQ3GEDm3XdgRiZRNui0iIuWGGfz9tOCR39/NgA/XhB2RSNmk5EJEypXKafDvMyGjOlw5JZi7RUQSS8mFiJQ79avBg2fB1j0wZArs2hd2RCJli5ILESmX2mfAP/rCorXwh5l6gkQkkZRciEi5dXobGHYCvPA5PPBB2NGIlB1KLkSkXLu+ezBJ3F1vB+NgiMihU3IhIuWaGdzdB45qANdNzTkxnYgUj5ILESn3qlWCB/qBE3Tw3LE37IhESjclFyIiQIs6cN/p8PkGuOV1dfAUORRKLkREIk5uGYzi+fIX8NDCsKMRKb2UXIiIRLmmG/zsCBj1Nsz9JuxoREonJRciIlHM4N6+0LIOXDsVVm8NOyKR0kfJhYhILjXS4YGzgpE7r3oV9uwPOyKR0kXJhYhIPtrWg7/3gYVr4M9vhR2NSOmi5EJEJIZ+beHK4+DRxcEoniISn1CTCzO72syWm9kuM1tgZr0KKX+0mc02s51mtsrM7jAzi1H2JDPbZ2YfJyd6ESkPhvcMpmj/w0z4bH3Y0YiUDqElF2Y2ABgDjAK6AHOBqWbWPEb5WsAMYC3QDbge+B0wLJ+ydYFHgZlJCV5Eyo20CjD2DKhVGX7zCmzeHXZEIqkvzJaLYcAkd5/o7p+5+3XAauCqGOUvBaoBl7v7x+7+PPBXYFg+rRcPAY8A85IUu4iUIw2rw/gzYdVWGDYNDmiALZEChZJcmFk60BWYnmvTdKBHjN1OBOa4+86oddOApkDLqLqvBhoDdyYqXhGRbk3h9l7w+nIY937Y0YiktrBaLhoAFQlucURbS5AY5KdxjPLZ2zCzo4ERwKXuXujDY2Y2xMyyzCxr/XrdTBWRgg06Fs5tB/fMg7c0g6pITGE/LZK7cdHyWVdYeQA3s8rA08Bv3X15XAd3n+Dume6emZGREVfAIlJ+mcFdp8KR9eH61+DbLWFHJJKawkouNgD7ydtK0ZC8rRPZ1sQoT2SfJkAH4OHIUyL7gDuAjpH3fRMSuYiUa9Uqwf39YN8BuOqVYKAtEckplOTC3fcAC4A+uTb1IXhqJD/zgF5mViVX+e+AFcAq4Gigc9RyP/Bl5N+x6hURKZLWdYMhwhevgztmhR2NSOoJ87bIvcAgMxtsZu3NbAxB58z7AcxstJlFP0r6JLADmGRmncysPzAcuNcDeyNPkRxcgHXA7sj7bSV7eiJSlvU9IpjkbPIn8ORHYUcjklrSwjqwu082s/rA7QS3ND4GznT37G5STYAjospvNrM+wDggC9gI3EOQpIiIlLibT4CPI60X7RpA1yZhRySSGsxdD2wDZGZmelZWVthhiEgps2kXnP007N4HUy4OxsQQKS3MbIG7Zya63rCfFhERKdXqVIEJ/WDLbrhaM6iKAEouREQOWfsM+Ntp8P53mkFVBELscyEiUpac0y54emTiB9AxAwZ2CjsikfCo5UJEJEGG94RezeH2N4NWDJHySsmFiEiCpFWAcWdAs5owdAqs0gieUk4puRARSaDaVeDBc2DXfhg8BXbsDTsikZKn5EJEJMHa1oN/nQ6frYffzgA98S/ljZILEZEkOLVV0AfjlaXwr/fCjkakZOlpERGRJPlNV1jyPdz7LrSsE0zXLlIeqOVCRCRJsqdo7940uD3y/qqwIxIpGUouRESSqHIaTDgreILkyimwYlPYEYkkn5ILEZEkq1sVJp0b/PtX/w3mIxEpy5RciIiUgJZ1YOJZ8O1WGDIlmOhMpKxSciEiUkK6NYO7+8D8VXDjNNh/IOyIRJJDyYWISAk6tx3c3gte/RLumKUxMKRs0qOoIiIl7MrjYMMOuH8BNKgGN50QdkQiiaXkQkQkBMN7wvc74J/zgwTjF8eEHZFI4ii5EBEJgRncdRr8sAv+702oWwXOOjLsqEQSQ30uRERCkj2LamZTuGEaTP8q7IhEEkPJhYhIiKpWgofPgaMbwtWvwszlYUckcuiUXIiIhKxmZXjkPGjXAIa+ArNXhh2RyKFRciEikgJqV4Ynfg5t6sKV/4N3vgk7IpHiU3IhIpIi6lSBJ/oHo3n++mV4++uwIxIpHiUXIiIppF7VoAWjRW341cswTZ08pRRSciEikmIyqsPkC6BjBlz1CrzwedgRiRSNkgsRkRRUpwo8/nM4vhncNA0e/TDsiETip+RCRCRF1UiHh8+F01rB/82Cf72nuUikdAg1uTCzq81suZntMrMFZtarkPJHm9lsM9tpZqvM7A4zs6jt/c1supmtN7OtZjbfzM5J/pmIiCRHlTS4vx/8/Ci4Zx787nXYsz/sqEQKFlpyYWYDgDHAKKALMBeYambNY5SvBcwA1gLdgOuB3wHDoor1Bt4A+kXqfBV4sbCkRUQklVWqCP/oCzceD89+Cr98CTbvCjsqkdjMQ2pjM7P5wGJ3vzJq3VLgOXf/Qz7lrwL+CjRy952RdbcDVwGHeYwTMbP3gDnufnNB8WRmZnpWVlaxz0dEpCQ8/xnc8jo0rw2Tzg1eRYrLzBa4e2ai6w2l5cLM0oGuwPRcm6YDPWLsdiJBkrAzat00oCnQsoDD1QQ2xohjiJllmVnW+vXr4wldRCRU57cPOnpu2AHnTYa5GmxLUlBYt0UaABUJbnFEWws0jrFP4xjls7flYWbXAIcBj+W33d0nuHumu2dmZGTEE7eISOhOOAxeGhA8UXLpi3B/ljp6SmoJ+2mR3F8Hy2ddYeXzW4+ZnQ/8HbjU3TVSv4iUKa3rwssD4cw2MPod+M0rsGV32FGJBMJKLjYA+8nb4tCQvK0T2dbEKE/ufSKJxWPAL9395UMLVUQkNdVIh7FnwB0/CWZTPfsp+FR3eCUFhJJcuPseYAHQJ9emPgRPjeRnHtDLzKrkKv8dsCJ7hZldBDwODHL35xIVs4hIKjKDX3eBp/vDjn1w7mS4fwHsPxB2ZFKehXlb5F5gkJkNNrP2ZjaGoHPm/QBmNtrMZkaVfxLYAUwys05m1h8YDtyb/aSImQ0Enoisf8vMGkeWeiV4XiIiJa5bM3jtEvhpSxj9Ngx8Hr7eHHZUUl6Flly4+2TgRuB2YBFwEnBmVP+IJsARUeU3E7RUNAWygHHAPQRJSrahQBrwT2B11PJC0k5ERCRF1K8GD/SDe/rApxvg9Cdg8ifq7CklL7RxLlKNxrkQkbLk2y1w8wx491vocRjceQocUTfsqCTVlKlxLkREJLkOqwVP9Yc7fwofrwtaMe6eB7v2hR2ZlAdKLkREyqgKBr84Bt74JfRrC/e9B30eh5nLdKtEkkvJhYhIGZdRHf75s6Alo1IFuOJ/QYfPhWvCjkzKKiUXIiLlRI/DYdql8OeT4csfguHDr3oFluc7QYJI8Sm5EBEpRypVhF8eC7MHBbOszloJpz4GN02DL74POzopK5RciIiUQzXS4aYTYPblMOhYmPpl0B9jyBT4ULdL5BApuRARKccaVoc7esPcK+D67jDvWzhnMgx4Dl5ZCnv3hx2hlEYa5yJC41yIiMC2PfDER/Doh/Dt1iD5uKQTXNwJGtcIOzpJtGSNc6HkIkLJhYjIj/YfgFkr4NHFMHtl8Fhrr+bw86PgZ0dA1UphRyiJkKzkIi3RFYqISOlXsQKc2jpYVm6Cpz+Blz6HG6ZB9UpwRhs4px2ceBikVww7Wkk1armIUMuFiEjBDji8twpe+BxeXQpb90DN9GCytNPbQO8WQUdRKT10WyTJlFyIiMRv1z5452t47St4fTn8sDNowchsEiQZvVpA+wbB7RRJXUoukkzJhYhI8ew/AFmrYcYymLMSPo+Ml9GgajBwV/dm0L0ptK2vZCPVKLlIMjOL+UG4O2aW59/5lcmvbEGvUcfPU3dhdebeL7/YCoo997b8Yok+Tqzzjd6e3/vCziH6uLmPEY+Cjl3U2GJ9FrE+l3jONb948zv3wq5jfseKdYzcdUbXU9jPZ35yfwb5/Tu/OHMfp7DPIVa9hcURq0z0+1jxx3uNC/t5jGefeM45ntgL2z9Z+8Rbbu02ePuboCPou9/C2u3B+jpVgpaNLo2hc2M4uhHUrlzooQ4pllSSyDgL+/4XoR516BQRkdTXqAac3z5Y3OGbLUFfjfe+g6zvgtso2VrXgU4NoWMGtM+ADg2CuVCkdFNyISIiSWMGzWsHywUdgnWbd8NHa2HRGvhwbZBwvPzFj/tkVIO29YLbKG3rBUubelC/alCfpD4lFyIiUqJqV4aTmgdLtk274LMN8On64HXp9/Dcp7B9749laqVDq7rQsk7Q4tG8DjSvFWxzV+KRSpRciIhI6OpUCcbMOPGwH9e5w+ptwYRqyzbCsk3BDK4LvoOXl0B0b4P24+GwWsHSrOaPr01qQtMawUijlTQeR4lRciEiIinJDJrWDJaTW+bctmsffLsFvt4Mp94YDFH+7dZg3cI1QUtItAoW3G5pUiPoE9KoerA0jiQeDasH2+tW1RMtiaDkQkRESp0qaUE/jDb1gvd39M65fdse+G5rsKzeBqu3wnfbYM02WL4pmKBty+689VaqAA2q/bhkRJb61YI+Hw0ir/WrQd0qag2JRcmFiIiUOTXS4cj6wRLLjr2wbnuwrN8R9e/I+/U7gv4fG3bAvgP511GrcpBs1IssdatEXiP/rlsluOVTp0pQfu/+8pGQKLkQEZFyqVqloHNoyzoFlzvgQSvHhh3w/c7I645gVNLvd/74+u0WWLwWNu6CPTGmqm8zNkh86lSG2pGko3blH19rR73WSv/xfa3KULMypFVI9KeQHEouREREClDBfmx9aBNHefegVeSHXbBpZ9D/Y+MuOBcYdgJs3gWbdsPGncHrmm1B8rJpF+yN0UKSrXqlINHIdvlLkcQjPUg+aqZDjcpBYlKzcpDI1EiPrI8slSsm/8kaJRciIiIJZAbV04Pl8Fo5t91wfOz93GHnviDR2LwrGA9ky+4fX7dEvZ8f2ef7nbBiUzCJ3NY9sVtMoqVVCJKUmkmcZE7JhYiISAowC27VVKsUPMVSkHsir1Muzrl+975IorEbtu2NvO4Jli17YHtk2bonGENkblLORMmFiIhImVE5LVgaVIuv/D+SFEeoXUPM7GozW25mu8xsgZn1KqT80WY228x2mtkqM7vDcs1aZGa9I3XtMrNlZjY0uWchIiIi0UJLLsxsADAGGAV0IWidmWpmzWOUrwXMANYC3YDrgd8Bw6LKtAJejdTVBRgN3Gdm5yfvTERERCRamC0Xw4BJ7j7R3T9z9+uA1cBVMcpfClQDLnf3j939eeCvwLCo1ouhwHfufl2kzonAI8Bvk3sqIiIiki2U5MLM0oGuwPRcm6YDPWLsdiIwx913Rq2bBjQFWkaVyV3nNCDTzCodSswiIiISn7A6dDYAKhLc4oi2Fjgtxj6NgW/zKZ+9bXnk9fV8yqRFjrk6eoOZDQGGRN7uBj7O78DR3TosxsPB2evzK1vYa6y646mzoPpyrWtgZhviPa+C6kzGORR0jDgVeH7xxlbY+3g+66L8vMQTT2RdA2BD7rIFHSO/bfH8fBZWTzw/N/mtK+znkxjnF08cBcVQ1OsZ789FEY5d5HOLJ/Z49k/GPvmUy3F+YSni/xdFkdDzS3ScCaivXSLiyC3sp0U813vLZ11h5XOvj6dMsMJ9AjABwMyy3D2zwGhLMZ1f6abzK73K8rmBzq+0M7OsZNQbVp+LDcB+gpaGaA3J25qRbU2M8kTtE6vMPuD7YkUqIiIiRRJKcuHue4AFQJ9cm/oQe0yPeUAvM6uSq/x3wIqoMrlvq/QBstx976HELCIiIvEJ82mRe4FBZjbYzNqb2RiCzpn3A5jZaDObGVX+SWAHMMnMOplZf2A4cK+7Z9/yuB84zMz+GalzMDAIuDuOeCYk5rRSls6vdNP5lV5l+dxA51faJeX87MffyyXPzK4Gfg80IehMeZO7vxXZNgk42d1bRpU/GhgHdAc2EiQTf4pKLjCz3gSDjnUkaNX4q7vfXxLnIyIiIiEnFyIiIlL2lJKZ4UVERKS0KDfJhSVhHpNUYGZ/MLP3zWyLma03s/+ZWadC9mlpZp7PcnpJxR0vMxuZT5xrCtmnVFw7ADNbEeNavBKjfEpfOzP7iZm9HPnc3cwG5dpukWv6XeT6zDKzjnHUG/qcQQWdm5lVMrO/mtliM9tuZqvN7EmLMZ1B1H4nx7ieRyX9hPLGUti1m5RPnO/GUW/o1y4SR2Hnl991cDMbV0CdKXH94vk9UNLfvXKRXFgS5jFJIScD4wlGNj2F4LHb182sXhz7nk7Q3yV7eSNJMR6qJeSM8+hYBUvZtYMgxuhzO45gTJZnCtkvVa9dDYL+UzcAO/PZ/nvgZuA6gnNfB8wws5qxKrTUmTOooHOrRnDt/hJ5PRc4HHjNzOIZT6gjOa/n0gTFXBSFXTsIBimMjvPMgipMoWsHhZ9fk1zL2ZH1hX0XIfzrdzKF/x4o2e+eu5f5BZgPTMy1bikwOkb5q4AtQNWodbcDq4j0U0nVheALtB84u4AyLQl+gWWGHW8c5zMS+LgI5UvttYvEehuwCahWBq7dNmBQ1HsjGCX3tqh1VYGtwG8KqOevwNJc6x4E5qXKucUo0yFyrY4uoMzJkTINwr5ehZ0fMAmYUsR6Uu7aFeH6TQSWFFImVa9fjt8DYXz3ynzLhSVvHpNUVZOgRWpjHGVfMLN1ZvaOmV2Q5LgORetIU+ZyM3vazFoXULbUXjszM+DXwOPuvqOQ4qXl2kVrRTDI3cHvYuQ6vUXs7yKU3jmDakVe4/kuZkVupcw0s58mM6hDdFLk5+4LM5toZg0LKV8qr52Z1QAGEiQY8Ui165f790CJf/fKfHJBwfOY5B7NM1vjGOWzt6WyMcAiggHFYtlGMFPsRQTNmjOByWZ2WdKjK7r5BGOVnAFcSfD5zzWz+jHKl+Zr14fgP4EHCyhTmq5dbtmff1G+i9n75bdP9pxBKSfyR809wP/cPfecSNGyZ4I+H+hPcAtwppn9JPlRFtlrwC+BUwma17sDb5hZ5QL2KXXXLuISoDLBrNoFSdXrl/v3QIl/98KeW6QkJWMek5RiZvcCJwEnufv+WOXcfQPBf3zZsiyYGOv3wOPJjbJo3H1q9PtIB7JlwOUEA7Hlu1uu9yl/7SKuBN5390WxCpSma1eAon4XY+2T3/rQRfpYPA7UAc4pqKy7LyH4hZRtnpm1JEgg30pSiMXi7k9Hvf3IzBYAK4F+wAsF7ZrrfcpeuyhXAi+5+/qCCqXi9Svk90CJfffKQ8tFsuYxSSlm9g/gYuAUd19WjCrmA20TG1Xiufs24BNix1rqrh1ApHn5XOJvho1WKq4dwbWBon0Xs/crFXMGRRKLp4BjgFPdvTjxlYrr6e7fEcxUXVCspebaZTOzzkAmxfsuQojXr4DfAyX+3SvzyYUnbx6TlGHB0OmXEPxAfV7MajqTa0r6VBS5JkcRO9ZSde2iDAJ2A08XUi4/nSkF1w5YTvCf1cHvYuQ69SL2dxFKyZxBkXvQkwkSi5+6e4GPTBegM6XgekZazJpRcKyl4trlMoTg/4rXi7l/Z0K4foX8Hij5717YvVpLqOfsAGAPMBhoT3A/ahvQIrJ9NDAzqnztyIV4GuhEcC9tC3Bz2OeSz7mNi8R2CkGGmb3UiCqT+/wuj/wQtgfaETTh7SEYfj30c8p1fncDvQn6IhwPTImcb6m/dlExG/AFuZ5oKo3XjqCXeufIsgO4I/Lv5pHtt0SuR//I9XmaIPGrGVXHo8CjUe9bAduBf0bOe3DknM9PlXMjuMX8EsFTScfl+i5WLeDcbgTOI/hLt2PkejvQP5WuXWTb3QQd/FoSPCUxj6DlIuWvXTw/m5Ey1YDNRD1VkauOlLx+xPd7oES/eyV6ccNcgKsJstHdBC0ZP4naNglYkav80QT3zHYRZKEjSMFHGSM/yPktI2OdH8EvqE8jPzRbgCzgsrDPJcb5ZX8B9hD8x/080KEsXLuoeH8auWbd89lWqq4dPz6al3uZFNluBI8Xr45cn9lAp1x1zAJm5VrXG/gg8v1dDgxNpXPjx0eE81sGxTo3gr4yXxKMu/ADMAc4M9WuHcFji9MIxkbYQ9DXYhJweGm4dvH8bEbK/Iqgyb9pjDpS8voV8LM3MqpMiX73NLeIiIiIJFSZ73MhIiIiJUvJhYiIiCSUkgsRERFJKCUXIiIiklBKLkRERCShlFyIiIhIQim5EBERkYRSciEiIiIJpeRCREREEkrJhYiUODO70Mx2m1mLqHVjzOwrM2sUZmwicug0/LeIlDgzM+B9YKG7X2lmvyWYp6Gnuy8NNzoROVRpYQcgIuWPu7uZ3Qq8YmZfAbcRTBWtxEKkDFDLhYiExszmAt2Bs919atjxiEhiqM+FiITCzE4BjiWYCnptyOGISAKp5UJESpyZHQvMBoYB/YAa7v6zcKMSkURRciEiJSryhMhc4AF3/5OZdQIWE/S5mBVqcCKSEEouRKTEmFk94B3gLXf/TdT6yUBzdz8xtOBEJGGUXIiIiEhCqUOniIiIJJSSCxEREUkoJRciIiKSUEouREREJKGUXIiIiEhCKbkQERGRhFJyISIiIgml5EJEREQSSsmFiIiIJNT/A2YSTy+p2cyZAAAAAElFTkSuQmCC\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": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEiCAYAAADJdLQBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABCD0lEQVR4nO3dd3hUVfrA8e+bXkijt4QQQKrSe1UX2+raFVdcXVGsP8u69l3XXlaxu6usura1APaGoIKANIP0EmpCDTUJpLfz++Pe6DhMkklI5s5M3s/zzANzyznvPXMz75x7z71XjDEopZRSvhLidABKKaWaFk08SimlfEoTj1JKKZ/SxKOUUsqnNPEopZTyKU08SimlfEoTTwASkUwRmel0HEopVR+aePyUiIwSkftFJNFH9d0oIlf4oq5jJSKXisgqESkWkSwReUBEwr1Yr7mI3CYic0Vkr4gcFpGfReQ6EQn1sLyIyM0islFESux/bxIRcVvuChEx1by6NuS2u9Q5SES+F5F8ETkoIm+JSOs6rN9VRD4VkTy7HT4VkS4elsusZrvecVuuh4g8Zrdnnojss+Mb76HMuTW0lxGRDvVrlfoTkVNEZLGIFIlItog8LyLN6rC+V59HDdv8sNtyg+0YVttl7haRL0RkkIcyq/uMjIiU1a9FGleY0wGoao0C/gG8AeT6oL4bgWy7Pr8lIpdjxfgl8AJwPPA3IBm4spbVRwCPA1/Z/xYB44F/AWOAS9yWvw+4H3gb+CcwFngOSAQe9FD+/cAWt2nZtcRUZyLSC5gLZAJ32vHcBgwQkcHGmKJa1m8LzAfKgAfsybcA80WkvzFmr9sqq4An3aZtdXt/FTAJ+BCYCkQBlwGzRORaY8wrLss+Arzqtn4o8B9gozFmV03xNzQROQlrn1iC1Q4pwF+A3iLyO1PLVfb1+Dy+B/7rNm212/s7sb4DpgPP22VOBpaIyFnGmK9clr0FcE+SzbH21W9qit0xxhh9+eELuAswQKqHeZnAzAaubwMw1+ntriXGSGAf8K3b9Pvttupfy/qdgU4epv/bXv8El2ntgGLgNbdl37Cnt3WZdoW9/jAftcOndjs0d5k2zo7hFi/Wfw4oAbq4TOsGlALP1mdfAwYBzdymRQArgf1ASC3rn27Hf1c92+SN+u6/dozrgAgPn+k5Dfl52NNe9qLMEa7x2NOSgN3Az16sf51d1wRf7JN1femhNj8kIvcDj9lvt7l0m8e5LTdERBbYhwd2iMgtHsqKEJG/i0iGfbgoW0ReFpdDeCKSCXQHxrrUlemy/oMiki4iOXZdS0TkD42w6bU5EWgFvOg2/SWsP7KLalrZGLPNGJPlYdZH9r+9XKadjZXoXnBb9gV7+tme6hCROPFw2K6hiEgc1pf0u8aYQ1XTjTFzgTXAxV4UcxHwtTHml96ZMWYT1q9jj+uLSLiIxFZXoDEm3RiT7zatFPgCaAnUdhhwItZn+K4X8TcYEekBnABMteOt8g7WkYYa27O+n4eIRIlIdHXlGmMWusWDMSYHmMNv99PqTASOYCVFv6OJxz99BHxg//9WrEMWlwHrXZbpDHwO/Ih1WGAL8IyInFK1gIgI8DFwDzAL+D/gTbus2fLreZFbgD1YvZ6qum6x58UD1wKLgHuxDmuFA5+IyKm1bYiIRIpISy9fte2PA+x/l7pONMbsxzr0M+CoNbzT3v73gFtdJViHmVwtx+oZeKprNnAYKBKRr+1DMA3teKz2X+ph3hKgX03tKCLtgbY1rN9WRNq5TR8LFAL59g+cO7z4rKq0B8qBvBpiisVK5POMMdu9LLehVLdPlQPLqH2fqs/nMRGrPQvtH4S1HSJ21Z7f7qdHEZHOWD2mj0wth12doud4/JAxZpWIrMD6tfSJMSbTw2LHAeONMd8CiMjrwHbgaqwkA9Y5i9OBk40xc6pWFJEfsM6RTADeNsZ8IiKPA3uNMb85aQzkAMnGmBKX9V8AVgB/pfZjyJdw9PHs6nTGOrRTnaoE4em8yW6X+V4TkQjgdnv9+W517TXGVLoub4ypFJG9bnUVYh/qwfqCHYD1Y2ChiAwyxmyua1w1qKp3j4d5u7HOrSQBB+u5ftUyVfNXAQuADKAF8CfgCaATcENNgYo1WGEC8FktX4DnArFY59J8rbb2GHqM67t/HguBacA2rMO51wGviUgLY4z7ebTfEJGRWD8Cnqklpkvtf51oT69o4glcm6qSDoAxpkREFgNpLstcBGwEVotIS5fpS4F84CRq2TmNMRVABfzyJd0Mq6f8A94d1vkG6wS+N2o7ER8NlLknA1uxPb+ungd6A2e7Jle7rBLPq/y2LmPMNKwvkyqfiMhXWL3R+7F+4TaUqno9xVbstswxr2+M+c0hVRH5L1ZP+zoRec4Ys9FTJSISg9UmxVi99ppMtJebUctyVWWHAwlukyOBcLf9HCDPGFPTyK7a2qO2faqu7TnSdQH7B+NPwP0iMtUY47FnaI+QexfI4tcBIdW5FNiFdVjOL2niCVyezlXkYB2vrnIc1rmb/dWU4dXwWxG5CuvLoyfgOpS41mdqGGP24PnXYH0UYX25iLHPoLqIsud7TUTuAa4B/maM+cxDXZHVrFprXcaYxSKyAPhdLTGEYp23cnXI/fi+W1xUE1uU2zINvr4xxojIFOD3WD9cjko8dmKYjpXQz6jp8JmItMFqo4+r+9L1YCTVf6m67+snYvVEq1Nbe9S2Tx1re5aJyHPA68Bw4Kjr8+zzSF8BccAYY8zh6sqzh1v3AJ6q5geaX9DEE7gqqpnumhhCsEbr3FzNsjUeKwYQkUuwhrl+hnWIZR/WMfs/A3/0Yv1ojv51Wp39dg+rOlUJrB2/Hhaq0p6jhzLXFNd1WMN6nzXGPFJNXb8TkRDXP2D7eH0bD/V7koX1ZVKTZKzDLq5q+rJ0bQN37bF+ZefUUF9t60Pt21b1o6e5+wy7fd4CTgUuNsZ8X0tZl2ANpXY/xFuTlRzdi74d69zVbR6WrYlre7h/Du2pvS2O9fOAmtszGquH2RPr0PqaWsqq6l3XpT19ThOP/2qIJ/RtAQYC33vx66e6+i7GOnF/jmsvQ0T+7GUMF9Nw53iW2f8OAT5xiaUl1iHG6d5UIiKXYY2EexvrXEx1dV2F1YNc4TK9P9Yw4Z+9qCqN6nubVbI5+ku0pi/L1VjX3wzh6BFgQ4EVNX3Wxphd9jmqIR5mDwWy7V5qTaoO53ratpexzutcbYz5sJZywDosdBDrF71X7NFd37pOE5GJQKTr4Wcvue5TC13KC8M6V1fbOcxj+jxsHtvT7jnOwBoo8AdjzEL3Fd2WD8X6e1ttjKkt4TrL6fHc+vL8wrqg0+O1KVRzbQXWCe5Ml/eX2WXc5GHZMCDJ5X06sNLDch9iJbAQl2lpWCfUjRfb0Q7rUIo3r6hayorC+uOc7Tb9fns7B7hMi8E65NDSbdlzsHpsnwJhtcRdArzqoY1LgHYu05p7WP9kO6bXatqmeu4bn2P1PF0/v3F2fX9xW7YHkOI27QWqv47neZdpiYB42G++s9uws9u8p+wY7vByO7rby/+rAdrkDep/Hc9qqr+O5zyXaeF2e7ZzW9+rz6Oa/SQWa7RqHi7XQWEdrZiGdWTDq2txsHqZXre/ky/HA9BXNR8MDLN3oq+xRhJNAFrb8zLxLvGEYA2nNlgJ5CasIdXPYZ18nOCy7MtAJdbdEi4BzrKnX26v/wXWldMPYH35r8CLxNMI7XKlHc/nWD2S5+0/zjfdlqv6w7/fZdpgfj30cRXWYQnXV5pbGQ/aZbyJdVX+W/b7B9yWWwf8D+tq82uwLkgtxerNJDdCG/QBCrCuE7kBa7j8ITuOGLdlDW5fyFhJNRvrEM+t9mu7Pc01oV6B9aX8uL1dd1V97sAjbmXeZE//2UO7TgRiPWzHQ/Y6IxqgTd5w3846rDve3ocW2Pv4w1jnZebgkniBVDveN+rzeWD9QFpol3818HesH3UGq4foWubT9vRZntqzmu14296Ojr7+u6xzmzsdgL5q+HCs62a22zuTAcbZ0zPxIvHY00KxrslZYf8x5WIdynkC6OCyXFus8zh5dl2ZLvNuxzrcVmz/cU20/4iMQ+1ymf2FWALssL/A3K/yHsfRiecKe1p1ryvcyhCsL+XNdl2b7ffuvYCHsb5wc7ASznas28Z0aIztt+scYn8xFthfcu/gcjcFl+WOSjz29OOwkvdh+/U50M1tmQFYPcMd9vYfsb84L61m36upbVM9rLMF2NJA7fGGp+2sw/qnYY32LAL2Yl2kHOe2TCoeEo+3nwdWgvsG67xQqf23+C1wuofy5tbUnh6Wj7E/n+999Xd4LC+xg1ZKKaV8Qu9coJRSyqc08SillPIpTTxKKaV8ShOPUkopn9ILSGvRsmVLk5qa6nQYTUZhYaHTIQSMrCzrgvdOnTo5HEngiImJcTqEJmPZsmUHjDHut4MCNPHUKjU1lfT0dKfDaDKWLVtW+0IKgMmTJwMwdepUhyMJHAMHDnQ6hCZDRDzdTxLQQ21KKaV8THs8SgWoSZMmOR2CUvWiiUepADV0aG3PKFPKP+mhNqUCVEZGBhkZGU6HoVSdaY9HqQA1ZcoUQAcXqMCjPR6llFI+pYlHKaWUT2niaSSb9+Xz4OfrKC3328eeK6WUIzTxNJIdhwp5/cdtfL9hr9OhKKWUX9HBBY1kzHGtaBsfxfs/7eC0Pu2cDkcFoRtuuMHpEJSqF008jSQ0RLhoUEdemLOZXblFdEiMdjokFWT69u3rdAhK1YseamtEFw5KBmB6+g6HI1HBaOXKlaxcudLpMJSqM008jSi5eQyjurZkevpOKir1EeOqYb300ku89NJLToehVJ1p4mlkFw9OZlduEQs2H3A6FKWU8guaeBrZ+F5taB4bwQc/bXc6FKWU8tr9n63lu/WNMypXE08jiwwL5bz+HZi9bi8H8kucDkcppWq1IfswbyzMZMv+/EYpXxOPD1w8OJmyCsPHP+9yOhSllKrV/xZvJyIshAsHJjdK+Tqc2ge6tYljYKck3v9pO1eN7oyIOB2SCgK33Xab0yGoIJRfUs7Hy3dx5gntSIqNaJQ6tMfjIxcPTmbL/gLSs3KcDkUFie7du9O9e3enw1BB5pPlu8gvKWfisE6NVocmHh/5/fHtaBYZxvtL9Zoe1TCWLFnCkiVLnA5DBRFjDO8szqJ3+3j6Jyc2Wj2aeHwkNjKMs/q258vVuzlcXOZ0OCoIvPbaa7z22mtOh6GCyLKsHDZkH2HisE6NekpAE48PTRicTHFZJZ+t2O10KEopdZR3FmcRFxnG2f3aN2o9mnh86ISOCfRoG8cHP+nhNqWUfzmYX8JXq7M5f2BHYiIad9yZzxKPiNwvIsbtle0y/yER2SAiBSKSIyLficgItzLmeijj/VrqvcLDOkZEohprW2uIhQmDk1m9K481u/J8Xb1SSlVrWvpOSisquXRoSqPX5eseTwbQzuV1vNu8G+xpo4BtwEwRaeNWxn/dyrjGi3oL3dZpZ4wprv9m1N+5/TsSERbCNL1xqFLKT1RUGt5dmsWwtOZ0axPX6PX5+jqecmNMtqcZxph3XN+LyF+ASUA/4BuXWYXVlVEDU491GkVCTDhn9GnLx8t3cc8ZPYkKD3U6JBWg7rnnHqdDUEFi3sb97DhUxJ2n9fBJfb7u8aSJyC4R2SYi74tImqeFRCQCmAwcBla4zZ4gIgdEZK2IPCUi3qTnaBHJEpGdIvKFiPQ/ts04NhcPTuFIcTlfr9njZBgqwKWmppKamup0GCoIvLM4i1ZxkZzSq61P6vNl4lkCXAGcDlwNtAUWikiLqgVE5EwRyQeKgVuB8cYY17vUvQtcCpwIPAScD3xUS70ZwJXA2cAldtk/iki36lYQkckiki4i6fv376/TRnpjWFpzUlvE8J5e06OOwbx585g3b57TYagAt+NQId9n7GPC4GQiwnyTEnx2qM0Y87XrexFZDGwFLgeetifPwTq01hIrOU0TkeHGmD12GVNdilgtIluBJSIywBjzczX1LgIWudS7EKsX9X/ATdWsMxWYCjBo0KAGf5COiHDR4GT+OTODrfvzSWvVrKGrUE3AO+9YR6fHjBnjcCQqkL23dDsCXDKk8QcVVHFsOLUxJh9YC3RzmVZgjNlsjFlsjJkElAFX1VBMOlDhWoYX9VbY63m9TmO4YEBHQkOED3SQgVLKISXlFXzw0w5O7tmG9onRPqvXscRjD2fuAdR0oiMEiKxh/vFAaC1luNcrwAl1WacxtI6P4qQerflw2U7KKiqdDEUp1UTNXJPNwYLSRr0vmye+vI7nKREZKyKdRWQoMAOIBd4UkXgReVhEhopIiogMFJHXgY7ANHv9LiJyn4gMEpFUETkDeB9YDvzoUs93IvKYy/t/iMipIpImIv2A17ASz8u+2vbqTBiczIH8Ur5bv8/pUJRSTdA7i7Po1CKG0V1b+rReX/Z4OgLvYZ3s/wgoAYYZY7KAcqA38DGwCfgcaAGMMcasstcvBU7GGlqdATwPzAJ+Zx8+q9IF61qdKolY52vW28t3sMtd2vCbWDdjj2tFm/hIfTqpUsrnNmQf5qfMHC4dmkJIiG8f1eLLwQUTaphXCJxby/o7gLFe1JPq9v5WrBFyficsNISLBiXz0pzN7M4t8ukxVhX4HnzwQadDUAGssR/2VhO9V5vDLhqUTKWB6ek7nQ5FBZi2bdvStq1vrrtQwSW/pJyPft7ZqA97q4kmHoclN49hdLeWvLd0O+U6yEDVwaxZs5g1a5bTYagA9MnyXRSUVnCZjwcVVNHE4wf+NDyV7MPFzFq3t/aFlbLNmDGDGTNmOB2GCjCuD3vr14gPe6uJJh4/cFKP1nRMiubNhZlOh6KUCnK+ethbTTTx+IHQEOGyYZ1Ysu0QG7IPOx2OUiqIve2jh73VRBOPn7hoUDKRYSG8uTDL6VCUUkHqYH4JX/voYW810cTjJ5JiIzi7X3s+Wb6LvMIyp8NRSgWh93/a4bOHvdXEuZSnjvKn4alMS9/J9GU7uGq0xydGKPWLf/7zn06HoAJIaXklby3KZHS3lj552FtNtMfjR/p0SGBQpyTeXpxFZWWD3xRbBZnExEQSExOdDkMFiC9X72bv4RKuHNXZ6VA08fibP41IJetgIT9sbPjnAKng8vnnn/P55587HYYKAMYYXp2/ja6tmzG2Wyunw9HE429O692W1nGRvLko0+lQlJ/TxKO8tWTbIdbuPsykUZ19fl82TzTx+JmIsBD+ODSFuRn72XagwOlwlFJB4NX522geG8G5/Ts4HQqgiccv/XFICmEhwtuLdGi1UurYbDtQwHcb9jJxaApR4aFOhwNo4vFLreOjOOP4dkxftoOCknKnw1FKBbD//riN8JAQJg535r5snmji8VOXj+jEkeJyPl6+y+lQlFIBKrewlOnpO/lDv/a0jotyOpxf6HU8fmpAShK928fz1qJMLh2a4tg9lZT/ev75550OQfm5d5dup6isgkl+MITalfZ4/JSIcPmIVDbuzWfx1kNOh6P8UFRUFFFR/vMrVvmX0vJK3lyYyaiuLenZLt7pcH5DE48f+0Pf9iTGhOtdq5VH06dPZ/r06U6HofzUV6v3sPdwid/1dkATj1+LCg/l4sHJzFqXza7cIqfDUX5m9uzZzJ492+kwlB8yxvDqgq10aRXL2OOcv2DUnSYePzdxqDUS5d0lOrRaKeWdpdsOsWbXYSaNSvOLC0bdaeLxc8nNYzi5ZxveW7qD4rIKp8NRSgWAVxdsIykmnPMG+McFo+408QSAy4encqiglC9X7XE6FKWUn8s8UMC36/cycVgnv7lg1J0mngAwsmsLurSK5S29f5tSqhZVF4xe5kcXjLrT63gCQNXQ6vs+Xcvy7Tn0T0lyOiTlB6ZOnep0CMrP5BWWMS19J2f19a8LRt1pjydAnDegI80iw3hL79+mlKqGv14w6k4TT4BoFhnGBQM78sWq3ew7Uux0OMoPvP3227z99ttOh6H8RFmFdcHoyK4t6NXevy4YdaeJJ4BcPiKV8krDWwu116Ng/vz5zJ8/3+kwlJ/4avUesg8X+31vBzTxBJTOLWM5tVdb3lqUqXetVkr9ouoJo2mtYhl3XGunw6mVJp4Ac83YNA4Xl/P+TzucDkUp5Sd+ysxh9a48rhzpH08YrY0mngDTPyWJIZ2b89r8rZRVVDodjlLKD7w6fyuJMeGcP6Cj06F4RRNPALp2bBq784r5YtVup0NRDtK7UyuAzfuOMHv9XiYO7UR0hH9eMOpOr+MJQOOOa81xbZrxyg9bOadfB31WTxOlz+NRAP+au4WosFD+PDLV6VC8pj2eABQSIkwe04UN2Uf4YeN+p8NRSjlkx6FCPl2xm0uGpNCiWaTT4XhNE0+A+kPf9rSNj+KVH7Y6HYpyyKuvvsqrr77qdBjKQf/+YQuhIkwek+Z0KHWiiSdARYSFMGlUZxZtPcjKHblOh6McsHTpUpYuXep0GMoh2XnFzEjfyQWDOtI2IbDO9WniCWAThiQTFxXG1Hna61GqqfnP/K1UGMN1Y7s4HUqdaeIJYHFR4Vw6tBNfr9lD1sECp8NRSvnIwfwS3l2ynbP7tie5eYzT4dSZJp4A9+eRqYSFhPDq/G1Oh6KU8pHXf9xGcXkF158YeL0d8GHiEZH7RcS4vbJd5j8kIhtEpEBEckTkOxEZ4VbGXA9lvO9F3eeLyDoRKbH/PbcxttEJbeKjOLd/B6al7+BgfonT4SgfSkhIICEhwekwlI/lFZXx1sIsTu/Tlq6t45wOp1583ePJANq5vI53m3eDPW0UsA2YKSJt3Mr4r1sZ19RUoYgMBz4A/gf0s/+dLiJDj3Fb/MbVY9IoKa/kTX1kQpPy5JNP8uSTTzodhvKxtxdlcqSknOvHdXU6lHrzdeIpN8Zku7x+uQjFGPOOMeY7Y8xWY8xa4C9AHFaycFXoVkZeLXXeAswxxjxijFlvjHkEmGtPDwpdWzdjfK82vLUok8JSvXmoUsGqsLSc1xZs48TurejTIXB7u75OPGkisktEtonI+yLicfC5iEQAk4HDwAq32RNE5ICIrBWRp0Sktr7mcGCW27RvgBEelg1Y145NI7ewjGl689Am48UXX+TFF190OgzlQ+8u2U5OYRk3nhS4vR3w7S1zlgBXABuA1sDfgIUi0tsYcxBARM4E3gdigD3AeGPMXpcy3gWygN1Ab+AxoC8wvoZ62wJ73abttad7JCKTsRIfKSkp3m2dwwZ2as6gTkn8Z/42Jg7rRFiojhsJdqtWrXI6BOVDxWUVTJ23leFpLRjYqbnT4RwTn307GWO+NsZMM8asMsZ8C5xp13+5y2JzsA6tjQBmAtNEpJ1LGVONMd8YY1YbY94HLgZ+JyIDaqve7b14mOYa61RjzCBjzKBWrVp5u4mOu2ZsF3blFvHl6j1Oh6KUamAzlu1k35GSgO/tgIPDqY0x+cBaoJvLtAJjzGZjzGJjzCSgDLiqhmLSgQrXMjzI5ujeTWuO7gUFvJN7tKZLq1he+WErxlSbV5VSAaasopKXf9hCv+RERnRp4XQ4x8yxxCMiUUAPrENq1QkBarrz3fFAaC1lLOLoQ3HjgYVehBlQQkKEa8Z0Yd2ewyzYfMDpcJRSDeTTFbvZmVPEjSd2DYq70fvyOp6nRGSsiHS2hzLPAGKBN0UkXkQeFpGhIpIiIgNF5HWgIzDNXr+LiNwnIoNEJFVEzsA6H7Qc+NGlnu9E5DGXqp8DThKRu0Wkh4jcDZwIPOuTDfexs/u3p3VcpN48tAlo06YNbdq4X22ggk1FpeFfczfTs108J/f0/8dae8OXgws6Au8BLYH9wGJgmDEmS0RisAYLXAm0AA4CPwFjjDFVZ1BLgZOBm4FmwA7gS+ABY0yFSz1d7HkAGGMWisgE4GHgAWALcLExZkljbaiTIsNCuXJUZx7/egNrduUF9JBLVbOHHnrI6RCUD8xck83W/QW8+Mf+QdHbAR8mHmPMhBrmFQI13k3AGLMDGOtFPakeps3A6mE1CX8cmsKL32/m33O38NKltY27UEr5K2MML87ZTFqrWE7v0672FQKEjrkNQvFR4Vw+ohNfrdlDRvYRp8NRjWTKlClMmTLF6TBUI/p+wz7W7znM9eO6EhoSHL0d0MQTtK4enUZsRBjPfbfR6VBUI8nIyCAjI8PpMFQjqertdEyK5ux+7Z0Op0Fp4glSiTER/HlkKl+tzmb9nsNOh6OUqqNFWw6yfHsu14ztQniQXRAeXFujfuOqUWnERYbx7Lfa61EqkBhjeGpWBm3iI7lwYEenw2lwmniCWEJMOFeO6sw3a/eydndt91JVSvmL7zfs4+ftudx88nFEhYc6HU6D08QT5K4c1Zm4qDCe/XaT06GoBtapUyc6derkdBiqgVVWGp78JoPUFjFcOCj4ejvg2+t4lAMSosO5alQaz3y7kdU78zi+o17XEyzuvfdep0NQjeDzVbvZkH2E5yb0C7pzO1WCc6vUb/x5VCoJ0eF6rkcpP1dWUcnTszfSo20cZ50QXCPZXGniaQLio8K5enRnvtuwj5U7cp0ORzWQRx55hEceecTpMFQDmp6+k6yDhdx+andCgui6HXeaeJqIy0ekkhijvZ5gkpWVRVaWPu48WBSXVfDcdxsZkJLIST2C455s1dHE00TERYUzeUwaczL2s3x7jtPhKKXcvL0oi72HS7jjtB5Bc0+26mjiaUIuH55K89gIntERbkr5lSPFZfxr7mZGd2vJsLTAf95ObTTxNCGxkWFMHpPGvI37WZalvR6l/MWr87eRU1jGHaf2cDoUn9DE08T8aXgnWsRG6LmeINC9e3e6d+/udBjqGB3ML+HV+Vs5vU/bJnO5g17H08TERIRx7dguPPLVen7KPMTg1OZOh6Tq6bbbbnM6BNUA/j13C0VlFdx2ynFOh+IzXvV4RGS1iKyq5bWysYNVDWPisE60bBbJM7O116OUk3bnFvHW4izOG9CRrq3jnA7HZ7zt8dT0ELVWWE8OjTz2cJQvREeEcu3YNB7+cj1Lth5kaBM4mRmM/v73vwP6JNJA9sL3mzDGcMvvujkdik95lXiMMQ+4TxORaOA24DJgE3Bnw4amGtPEYZ14Zd5Wnvl2I+9PHu50OKoe9u7d63QI6hhsO1DAtPSdXDasEx2TYpwOx6fqPLhAREJEZDKwGZgE/B/Qzxgzs6GDU40nKjyU68d1YfHWQyzccsDpcJRqcp6evZGI0BBuOLGr06H4XJ0Sj4icA6wDHgOeAXoYY94yxphGiE01skuGpNAmPpJnZ1vdfaWUb6zbfZjPV+7mylGptIpremcpvB1cMFJEFgD/Az4BuhhjnjLGlDRmcKpxRYWHcsOJXVmaeYi5G/c7HY5STcZTszKIjwpj8pguTofiCG8HF8wHioCpwD7gSk+3dDDGPN1woSlfmDA4hdcXbOOxr9YzumtLwoL0NuzB6IQTTnA6BFUP6ZmH+H7DPu44rTsJ0eFOh+MIbxPPdsAA59SwjAE08QSYiLAQ7jq9B9e+8zPT0nfyx6EpToekvHTjjTc6HYKqI2MM//wmg1ZxkVwxItXpcBzj7ai21EaOQzno1N5tGZyaxNOzN/KHfu1pFqnXFSvVGL5Zm83SbYd46Jw+xEQ03b8zPa6iEBHu/X0vDuSX8MoPW5wOR3np9ttv5/bbb3c6DOWl4rIKHv5yPT3axnHJ4GSnw3GUt4MLTheRTBE56kZCIpJgzzul4cNTvtIvOZE/9G3Pf+ZvZU9ekdPhKC/k5eWRl5fndBjKS68t2MbOnCLuO7NXkz+X6u3W3wg8aYw5ai+3pz0B3NyQgSnfu/3U7lQaePKbDKdDUSqoZOcV89KczZzWuy0jurZ0OhzHeZt4TgC+rWH+90DfYw9HOSm5eQx/HpnKx8t3sWaX/pJWqqH8c+YGyisN95zR0+lQ/IK3iacVUFnDfAPoDb+CwA0ndiUpJoJHvlyvF5Uq1QB+3p7DR8t3cfXozqS0aFq3xqmOt4lnJ1avpzonALuOPRzltPiocG75XTcWbT3Id+v3OR2OqsGQIUMYMmSI02GoGlRWGh74fB2t4yK5flzTuzVOdbxNPF8CD9k3Bv0NEYkBHrSXUUHgkiEppLWK5dGv11NWUVNHVznpqquu4qqrrnI6DFWDj5bvYuWOXO46vQexepnCL7xNPI8ACcAmEblTRM62X3cBG+15jzZWkMq3wkNDuOf0nmzdX8D7S7c7HY5SASm/pJwnZm6gX3Ii5/Tr4HQ4fsWrxGOM2QeMAFZhJZiP7dcj9rSRxhi9R3sQOblna4alNeeZbzdxuLjM6XCUBzfddBM33XST02Goarw0ZzP7j5Twj7N6ERJy9C3GmjJvr+M5AdhhjDkDaAkMBYYBLY0xZxhjMhsvROUEEeFvv+9FTmEp/5qjF5X6o+LiYoqLi50OQ3mQdbCA1+Zv47wBHeifkuR0OH7H20Nty7ESDsaYHOB+rESU00hxKT/Qp0MC5/bvwOs/bmPHoUKnw1EqYDzy5XrCQoU7T+vhdCh+ydvE495PHAMcNdBABZ/bT+2OYN3GXSlVuwWbDjBr3V5uOLErbeKjnA7HLzXt+zaoWrVLiObq0Wl8umI3K3bkOh2OUn6tvKKSB79YS3LzaCaN6ux0OH7L28Rj7Jf7NK+JyP0iYtxe2S7zHxKRDSJSICI5IvKdiIyopiwRkZl2GRfUUu8VHuo1IqI/Rbx07bgutGwWySNfrtOLSv3I6NGjGT16tNNhKBfvLt3Oxr353HtGL6LCQ50Ox295O7BcgHdEpOqJo1HAf0TkNwf+jTF/qKWcDGCcy/sKt3k3ANuwDuPdCswUkW4eRszd5rZubQqB3zzqzxijZ2W91CwyjL+MP457Pl7NN2uzOa1PO6dDUsBll13mdAjKRU5BKVNmbWRElxac2ruN0+H4NW8Tz5tu79+pZ33lxphsTzOMMb8pU0T+AkwC+gHfuEwfhHVD0oGAt0O4TXX1Ku9cNKgjby7M5KEv1jPmuFZN+lkiSnny7LcbOVJcxn1n9cLTE5rVr7x9ENyfG6i+NBHZBZQCS4B7jDFb3RcSkQhgMnAYWOEyPQ54D7jGGLOvDh9utIhkAaF2eX83xiw/hu1ocsJCQ3j43D5c+PIinvt2E3frzQ4dN3nyZACmTp3qcCQqI/sI7yzZzqVDO9GjbbzT4fg9Xw4uWAJcAZwOXA20BRaKyC83FxWRM0UkHyjGOtQ23u0w28vATGPMV3WoNwO4EjgbuMQu+0cR6VbdCiIyWUTSRSR9//79dagquA1Obc6Ewcm8umAb6/ccdjocpfxCZaXhno9XExdlHZJWtfNZ4jHGfG2MmWaMWWWM+RY4067/cpfF5mAdWhsBzASmiUg7ABG5DOvRC3V65KIxZpEx5k1jzApjzHzgYmAL8H81rDPVGDPIGDOoVatWdaku6N11eg8So8O55+PVVFbqQAOl3lmSxbKsHP7++14kxUY4HU5AcGw4tTEmH1gLdHOZVmCM2WyMWWyMmQSUAVV3QTwZ6AXki0i5iJTb0z8QkQV1qLcCSHetV3kvMSaCe3/fk+Xbc3nvJ72Pm2raduUW8cTXGxjdrSXnDdD7sXnLscRjD2fuAeypYbEQINL+/71Yj1/o5/IC+CvwpzrUK3Y5NdWranBu/w6M6NKCx7/ewL4jOjhQNU3GGP728WoqDTx67vE6oKAOfJZ4ROQpERkrIp1FZCgwA4gF3hSReBF5WESGikiKiAwUkdeBjsA0AGPMLmPMGteXXfQO1wEK9vU/j7m8/4eInCoiaSLSD3gNK/G87JstDz4iwsPn9KGkrJKHv1jvdDhN1vjx4xk/frzTYTRZn63czZyM/fz11O4kN9cHvNWFL8fEdsQakdYS2A8sBoYZY7LsZ/r0xhoE0AI4CPwEjDHGrKpjPV2AHS7vE4GpWIMZ8rDuOzfGGLO0/pui0lo14/oTu/Dst5u4YGBHxhyn58J87cILL3Q6hCbrUEEpD3y+jr7JiVwxItXpcAKOzxKPMWZCDfMKgXPrUeZRfVtjTKrb+1uxRsipBnbduC58tmI3f/90Dd/cMkav1PaxqjtTR0XpTTh87eEv1nG4qIwnzj+eUH3kQZ3pvdpUvUWGhfLwuX3IOljIi99vdjqcJkefx+OMHzbu56Plu7huXBe9ZqeeNPGoYzKiizWa55V5W9i874jT4SjVqApKyrnno9V0aRXLjSd1dTqcgKWJRx2ze8/oSWxkGPd8vEZvIqqC2pRZG9mVW8Tj559AZJgeWq4vTTzqmLVoFsndp/dg6bZDTF+20+lwlGoUy7fn8N+F27hsWCcGpzZ3OpyApolHNYgLByYzODWJR79az8H8ktpXUCqAlJZXcteHq2kbH8Udp3V3OpyAp4lHNYiQEOGRc48nv7icR7/a4HQ4TcJZZ53FWWed5XQYTcLLP2whY+8RHj6nD3FR4U6HE/D03vaqwRzXJo5rxqbx0pwtXDCwI8O7tKh9JVVvmnR8Y/O+I7z4/WbOPKEdJ/fU5+w0BO3xqAb1fyd1I6V5DPd+vJrisro8q0/VVW5uLrm5uU6HEdQqKw13fbiamMhQ7v9Db6fDCRqaeFSDigoP5eFz+rD1QAFPfZPhdDhB7Y477uCOO+5wOoyg9r8lWaRn5fC33/eiZbPI2ldQXtHEoxrcmONa8afhnXh1wTZ+3HzA6XCUqpfMAwU8bt95+ny983SD0sSjGsXdp/ekS6tYbpu2ktzCUqfDUapOSssruen95YSFhvDE+SfonacbmCYe1SiiI0J5bkJ/DuSXcO8nemGpCixTZmewamceT5x/PO0To50OJ+ho4lGNpk+HBG4dfxxfrtrDJyt2OR2OUl5ZsOkAr/ywlUuGpHBan3ZOhxOUdDi1alTXju3C3Ix93PfJWgZ1aq7PLWlAF1xwgdMhBJ2D+SXcOm0FXVs3474zezkdTtDSHo9qVKEhwtMX9cMAt01bSUWlHnJrKKeccgqnnHKK02EEDWMMt89YRV5hGc9P6E90hN6LrbFo4lGNLrl5DA+e3ZulmYd4Zd4Wp8MJGtnZ2WRnZzsdRtB4c2Em32/Yx91n9KBXe33cQWPSxKN84tz+Hfj98e14etZG1uzKczqcoHDfffdx3333OR1GUFi/5zCPfr2Bk3q01ieK+oAmHuUTIsIj5/ahRbMIbn5/OUWlelcD5R+KSiv4v/eWkxgdzpMX6NBpX9DEo3wmMSaCKRf2Y8v+Ah7/er3T4SgFwINfrGPL/nyevqgfLfTuBD6hiUf51KhuLZk0qjNvLspiTsY+p8NRTdzMNXt4b+l2Jo9JY1S3lk6H02Ro4lE+d/up3eneJo47ZqzSZ/cox+zOLeLOD1dzQscEbhuvz9jxJU08yueiwkN55uJ+5BWWcfdHq/WuBvU0ceJEJk6c6HQYAami0nDLBysor6jk+Qn9iQjTr0Jf0tZWjujVPp7bT+3OrHV7ef+nHU6HE5DGjBnDmDFjnA4jIL00ZzNLtx3iwbP7kNoy1ulwmhxNPMoxk0Z1ZlTXlvzjs7Ws2JHrdDgBJzMzk8zMTKfDCDg/ZR7iue82cXa/9pynd512hCYe5ZiQEOH5S/rTOi6Sa99exr4jxU6HFFAeffRRHn30UafDCCi7cou47p1lJCdF89A5fXTotEM08ShHNY+NYOplg8gtKuX6d36mTG+poxpJUXklV72ZTkl5Ja9ePpj4qHCnQ2qyNPEox/VqH8+TF/QlPSuH15cfdjocFYQqjeGFpXlkZB/mhUv607V1M6dDatI08Si/cFbf9lw3rguzthYxa0uh0+GoIDNtXT5LdpVwzxk9Gde9tdPhNHmaeJTf+Osp3enfNoLXlh9mwwF9aqlqGD/uKGL6ugJOSo1m0qjOToej0MSj/EhoiHDL0ERaxoby5MJcDhbp/dxqMmnSJCZNmuR0GH5tS04ZLy7No0eLcCYPiNfBBH5CE4/yK80iQrhzRCLFFYZ//phLaYUONqjO0KFDGTp0qNNh+K2cogoe/zGH+KgQ7hiRSHioJh1/oYlH+Z2UhHBuGpLA5pwyXlmWp3c2qEZGRgYZGRlOh+GXSioMTyzMpaDUcPfIJBKi9KFu/kQTj/JLQztEcVGvWOZmFfP1Zh1s4MmUKVOYMmWK02H4HWMML6fnselQGTcPTSA1UYdN+xtNPMpvXdirGYPbR/LflUdYs09vJqq880lGAfO2F3NJ72YM7RDldDjKA008ym+FiHDTkATaNwtlyqJc9hXoYANVs592F/O/1fmMTI7i/J56DzZ/pYlH+bWY8BDuHJlEuYF/LsyhqKzS6ZCUn9qeV8azS/JISwrjhkEJOoLNj2niUX6vfVwYtw5NJCuvnCcW6kg3dbQDhRU8uiCH6DDhzpFJRIZp0vFnYU4HoJQ3BrSL5MbBCTy/NI8pi3K5fUQiYSFN+8vlhhtucDoEv5BTXMH9PxyioNRw/7jmtIjWEWz+zmc9HhG5X0SM2yvbZf5DIrJBRApEJEdEvhOREdWUJSIy0y7jAi/qPl9E1olIif3vuQ25bco3xnaK5uoB8aTvKeGFpXlUNPFh1n379qVv375Oh+GowyWVPPBDDjlFldw7OokuSTqCLRD4+lBbBtDO5XW827wb7GmjgG3ATBFp46Gc2wCvzjSLyHDgA+B/QD/73+kiolfeBaDTusQw8fhmLNhRzH9+Ptykr/FZuXIlK1eudDoMx+SXVvLgvEPszS/nrlGJ9GgZ4XRIyku+PtRWbozJ9jTDGPOO63sR+QswCStZfOMyfRBwMzAQ2OtFnbcAc4wxj9jvHxGRE+3pl9QtfOUPzu3RjMIyw0cbCogOE/50QlyTPJH80ksvATB16lSHI/G9wrJKHp6fw47D5dw1MonjW0c6HZKqA1/3eNJEZJeIbBOR90UkzdNCIhIBTAYOAytcpscB7wHXGGP2eVnncGCW27RvAI+H8ex6JotIuoik79+/38tqlC/9sU8zTusSw2cbC/lwfYHT4SgfKi6v5NEFOWzNKeO2YYn0b6tJJ9D4MvEsAa4ATgeuBtoCC0WkRdUCInKmiOQDxcCtwHhjjGuv5mVgpjHmqzrU25aje0Z77ekeGWOmGmMGGWMGtWrVqg5VKV8RESb1j2NcpyjeW5vPl5s0+TQFpRWGx3/MJeOAdVeCIXqBaEDy2aE2Y8zXru9FZDGwFbgceNqePAfr0FpLrOQ0TUSGG2P2iMhlQF9gUH2qd3svHqapABMiwvWDEigqN7y+4gjRYcJJnWOcDks1krJKw1OLclmzr5QbBycwMjna6ZBUPTl2HY8xJh9YC3RzmVZgjNlsjFlsjJkElAFX2bNPBnoB+SJSLiLl9vQPRGRBDVVlc3TvpjXenR9Sfi40RLh1aCJ920Tw7/TDLNpZ7HRIqhFUVBqeWZzLsj0lTB4Yz7hUTTqBzLHreEQkCuiB1cupTghQdQD3XuApt/mrgb8Cn9ZQxiJgPPCky7TxwMK6xKv8V3iocMeIRB6al8Ozi3OJGpXUJI7733bbbU6H4BMVxvD80jyW7Crhz/3iOCVNe7WBzpfX8TwlImNFpLM9lHkGEAu8KSLxIvKwiAwVkRQRGSgirwMdgWkAxphdxpg1ri+76B3GmK0u9XwnIo+5VP0ccJKI3C0iPUTkbuBE4FkfbLbykaiwEO4ZlURyQhj/XJjD2v3B/wTT7t270717d6fDaFSVxvBy+mEW7Cjm0uObcWY3vf9aMPDlobaOWCPSMoCPgBJgmDEmCygHegMfA5uAz4EWwBhjzKo61tMF6xohAIwxC4EJWOeSVgF/Ai42xiw5pq1Rfic2IoS/j2lO65hQHpp3KOgPuy1ZsoQlS4J3Ny6rMLywNI/vM4u4sFcs5/Vo5nRIqoH4cnDBhBrmFQJ1vpuAMeaoizeMMakeps3A6mGpIJcQGcJDJ7bg8R9zmLIol8v7xnHWccH5K/m1114DCMqnkOaXVto91zIu6d1M7zQdZPQmoSroxEeG8I+xzRnSIZI3Vh7hvysOU9mE73AQaLLzy7n7+4NkHCzjlqEJXNCrWZO8QDiYaeJRQSkyVLhteCJndI3hi02FPL1Y72odCDYeLOXu7w9xuLiS+8Y0Z3SKjl4LRnp3ahW0QkW4sl8crWJDeXPlEXKKDnHXyCTiIvX3lj9atLOY55fkkhQdyr2jm9MhTr+egpX+BaqgJiL84bhY/jIsgc05Zdwz5yB7C8prX1H5jDGGzzIKmLIol9TEcB47uYUmnSCnn65qEkYmR5MUFcoTP+Zwz3eHuHtUEl2bB/Yt9O+55x6nQzhmFZWG11Yc5pstRQzrEMlNQxOJDNXzOcFOezyqyejVKoJHTmpBeCj8Y+4hlu0pcTqkY5KamkpqaqrTYdRbUXklj/+Yyzdbiji7eyy3Ddek01Ro4lFNSsf4MB47qQXt40J5/MccZm8tdDqkeps3bx7z5s1zOox6OVRUwd/nHGJFdgmTB8TzpxPiCNGRa02GHmpTTU5SdCgPjmvOlMW5vLzsMBsOlDKpfzwx4YH1O+ydd6xHWI0ZM8bhSOpmRXYJL/2UR2G54e5RSQxoF/y3N1K/pYlHNUnR4SHcPTKJGevzmbGugHX7y7hpaAI99SmWjaak3PDWqiPM3FJIx7hQ7h2dRGpiYJ9nU/WjiUc1WaEhwsW94+jXNpLnluRx35xDnNszlot6NSMsRA/7NKRNh0p5fkkeu/MrOLNbDH88Pk7P5zRhmnhUk9e9RQRTxrfg9RVH+HB9ASuyS7h5aKIO6W0A5ZWGD9fnM2N9AUlRIfxjTBIntNFDa02d/mUphXXo7YbBCQxsF8m/l+Xx19kHuKJvPKekRevtWupp15FynluSy5accsakRHFV/3hiIwLrPJpqHJp4lHIxrGMUx7UI58Wf8pj682GW7Snh+kHxJEaFOh3aUR588EGnQ/Co0hhmbink7VVHiAgV/jo8keEd9RHV6leaeJRy0zw6lL+NTmLm5kLeWnWEW2cd5PpB8Qxu719fnm3buj9Y13kHiyp46ac8Vu4tpX/bCG4YlEBStP8lbeUsTTxKeRAiwhndYunTOoLnluTx+I+5jEmJ4uLezWjbzD/+bGbNmgXAKaec4nAk1rNzZm8r5P01+ZRXwuQBephSVc8//oKU8lMpCeE8fnILpq/P5/OMAn7cUcz4tBgu6Bnr+C/5GTOsR0w5mXgqKg1zs4qYti6fA4WV9G4VzrUDE2ivAzNUDXTvUKoW4aHCH/vEcVqXGKavy2f21kK+zyzkjK6xnNsjlmZN8IR5pTEs3FHMB2vz2Z1fQZekMK4dmEC/NhHay1G10sSjlJeaR4dyzcAEzu4eywdr8/k0o4BZWws5p3ssZ3SLITos+BOQMYb0PSW8tyafrLxykuPDuGNEIkPaR2rCUV7TxKNUHbVtFsbNQxM5p0cZ763J5901+Xy5qZALesYyPi2G8CC9MHLV3hLeXZPPpkNltI0N5eYhCYxMiSJUE46qI008StVTp4Rw7hqZRMbBUv63Op/XVhzhs40FXNirGSOTo4gKgh5QpTGsP1DGtLX5rNlfSovoEK4dGM+JqdF6dwdVb2L0WfQ1GjRokElPT3c6jCZj2bJlTodQL8YYVu0r5X+rj7Alp5yoUGFwh0hGp0TRt01ko3xJ5+bmApCYmNig5RpjyMwrZ8H2YhbsKOJAYSXxkSGc3yOWU7rEEBHAPbqBAwc6HUKTISLLjDGDPM3THo9SDUBE6NsmkhNaR7DuQBkLthexcGcx87cX0yxCGNExilEp0fRsGd5gt/9v6ISz+0g5C3YUs2B7EbuOVBAi0LdNBJf0jmZox8gmcQ5L+YYmHqUakIjQu1UEvVtFcGX/eFbtLWH+9mJ+yCpm1tYiWkSHMDI5itEp0XRODDumE/Kff/45AGeddVa9yzhQWMGPdrLZmluOAD1bhvP7brEM7xhFfKQmG9XwNPEo1UjCQ4SB7aIY2C6K4vJK0ndbSeirTYV8trGQ9s1COaFNJJ0SwkixX3V5JlBdE09ZpWH3kXK251mv9QdKWX+gDIAuSWFc3jeOkR2jaBGjdxpQjUsTj1I+EBUWwqiUaEalRHOktJLFO4v5cUcxc7OKKC7/9Txrq5gQUhLCSUkI+yUhtY8LI7wO54gqjWFfQYWVYA7/mmh2Hymnwq4qVCA5PowJva2BEHrBp/Il3duU8rG4iBDGp8UwPi2GSmM4UFhBlp0csvLK2ZFXzorskt8kidaxobif0998yOqt3Dxz/y/TDHCwsJLiil+TWeuYUFISwhjcPvKXnlX7ZmFBO+xb+T9NPEo5KESE1rFhtI4NY3D7X6e7HhbLyitnb34F7uNPN4ZZiSM5wfUpnoZ+bUNJibcSTHJ8GNEB9khvFfw08Sjlh8JDhE4J4XRKCGd0Ncts/K/15/vX4Yk+i0uphqCJR6kA9fzzzzsdglL1oolHqQAVFeVfzwdSylt68FepADV9+nSmT5/udBhK1ZkmHqUC1OzZs5k9e7bTYShVZ5p4lFJK+ZQmHqWUUj6liUcppZRPaeJRSinlU/o8nlqIyH4g6xiLaQkcaIBwgpW2T820fWqm7VMzp9qnkzGmlacZmnh8QETSq3sgktL2qY22T820fWrmj+2jh9qUUkr5lCYepZRSPqWJxzemOh2An9P2qZm2T820fWrmd+2j53iUUkr5lPZ4lFJK+ZQmHqWUUj6liUcppZRPaeKpIxG5XkS2iUixiCwTkeoeEImIRInIGyKySkTKRGRuNcuNtcsqFpGtInJto21AI2vo9hGRcSJiPLx6NOqGNJI6ts84EflURPaISKHdTld6WK6p7j+1tk+w7T9Q5zbqJSJzRGSvy/7xqIhEuC3n031IE08diMjFwHPAo0B/YCHwtYikVLNKKFAMvAh8WU2ZnYGv7LL6A48BL4jI+Q0bfeNrjPZx0Rto5/La1BAx+1I92mcEsBq4AOgD/BuYKiJ/dCmzKe8/tbaPi4Dff6BebVQKvAmcAnQHbgEmAQ+7lOn7fcgYoy8vX8AS4D9u0zYBj3mx7ovAXA/TnwA2uU17FVjk9Pb6SfuMAwzQ0untc7J9XJafBnyo+4/X7RM0+08DttHTrvuHE/uQ9ni8ZHdNBwKz3GbNwvrlVV/DPZT5DTBIRMKPoVyfasT2qZJuH1L5TkRObIDyfKoB2yceyHF5r/vPb7m3T5WA3n+gYdpIRLoCpwE/uEz2+T6kicd7LbEODe11m74XaHsM5batpswwu85A0Vjtswe4DjgfOA/IAL4TkTHHUKYTjrl9RORM4GR+e0Gg7j+2atonWPYfOIY2EpGFIlKM1TtaANzjMtvn+1BYYxQa5NyvuBUP0xqiTE/TA0GDto8xJgPry6LKIhFJBf4KzKtvuQ6qV/uIyEjgXeAmY8xSL8r0ND0QNGj7BOH+A/Vro4uBOKAv8CRwJ9a5nJrK9DS9QWji8d4BoIKjf1m05uhfC3WRXU2Z5cDBYyjX1xqrfTxZAkxo4DIbW73bR0RGYZ38vc8Y82+32U1+/6mlfTwJxP0HjqGNjDE77P+uE5FQ4FURedIYU44D+5AeavOSMaYUWAaMd5s1Hms0SH0tAn7nocx0Y0zZMZTrU43YPp70wzqEEjDq2z72IaGvgQeMMc96WKRJ7z9etI8n/Qiw/Qca9G8sBKvTEWq/9/0+5PQojUB6YXVXS4GrgJ5YwxrzsR54BFbX9Tu3dXph7ejvA+n2//u5zO8MFADP2mVeZddxvtPb6yftcwtwDtANa0jsY1jd//Oc3t7Gbh+sEVkFWIdG2rq8Wun+43X7BM3+U882ugy4EOgBpAEXAbuA953chxxvyEB7AdcDmUAJ1q+PMS7z3gAy3ZbPtHf037zclhkL/GyXuQ241unt9Jf2Ae4ANgNFwCFgPnCG09vpi/ax3x/VNh7asEnuP960T7DtP/Voo0vsfeMIVoJaizWwINrJfUjvTq2UUsqn9ByPUkopn9LEo5RSyqc08SillPIpTTxKKaV8ShOPUkopn9LEo5RSyqc08SjlB1weWBZIN/ZUql408SjlABGZKyIvBkq5SjUkTTxKKaV8ShOPUj4mIm9g3aLkBvvwmgFS7dl9RWSJiBSKSLqIDHBbd4SI/GDP3yUi/xaR+OrKFZFUEQkVkddEZJuIFInIJhG5Q0T07185Qnc8pXzvZqw7Av8XaGe/qm5b/xhwFzAA65b0/xMRARCR47GeFPkZ1nNVzsO6qerrtZQbgnVjyIuwbgJ5L9b9uv7ceJuoVPX0Xm1KOUBE5gJrjDE32u/HAXOA04wx39jTRmI9LTLZGLNTRN4Cyowxk1zK6QcsB9oYY/a5l1tD/Y8Dg4wx7rfDV6rR6YPglPIvq1z+v9v+tzWwExgIdBWRi12WqXpSZBdgX3WFisi1WLe77wREA+FAVgPFrFSdaOJRyr+4Pnir6nBEiMu/rwLPeFhvV3UF2onqWazHPS8EDgM3AOceY6xK1YsmHqWcUcqvT4D01s9Ab2PM5jqWOwpYYoz5ZZi1iHSpY91KNRgdXKCUMzKBIfaos5Z497f4hL3OyyLSX0S6isiZIvJKdeXaI9c2AgNE5HQR6SYif8ca/aaUIzTxKOWMp7B6J+uA/UBKbSsYY1YBY7CGXv8ArMQaBbe3lnJfAaYB7wI/2etPaZCtUKoedFSbUkopn9Iej1JKKZ/SxKOUUsqnNPEopZTyKU08SimlfEoTj1JKKZ/SxKOUUsqnNPEopZTyKU08SimlfOr/AW62CCm3Fl1VAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Make scan of lnL (for theta, if free)\n", "if not(m.fixed['theta']):\n", " plt.figure()\n", " m.draw_mnprofile('theta')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABZa0lEQVR4nO3dd3hT5dvA8e+d0b3L3nspgixBEFBEcG/c4MQt7oEDlOXEAQ5wob74c6DiniBDdkGUJUMQZEP3Tprc7x8nYMUCLbQ9Tfp8ritXmpNzTu4cQu48W1QVwzAMwzgch90BGIZhGMHBJAzDMAyjVEzCMAzDMErFJAzDMAyjVEzCMAzDMErFJAzDMAyjVFx2B1BRatSooU2aNLE7DJv4QYuAIsAH6iv2d/H7osBzviN4DQeIC3AF7t3FHruLbQ/Zj5hhhKSlS5fuVdWaJT1Xaf+bReQvoHEJT32jqmce5Jj2wESgG5AGTAJGaSkGjzRp0oSUlJQjD7iKU/VB0Qbw/o56fwffFvDtAv8u0NwD9nYAYYAbHImBW9L+mzgSQOJAYsERHbiPte4BNBv82YH7DPCnof408KeCfy/4doN/j/U8AH7AE7g5Aq9TAxw1wZEMzlqIsz4464GzPjjqIY6oSrluhmEcmohsPthzlfnzryvgLPa4LrAU+KiknUUkDvgRmBM4tjUwBcgFnqvIQKsi9aeD51fUuxQ8y6FoFWie9aTEgaspuFqAsyfiqA3OWv8kB9l3H42IlEs8JZ1FteCf5OHfA769qH+vlVT8e61tng3g34NS9O9jJRHcrcF9LOJqB+5jwdkIEVNrahhVRaUlDFXdU/yxiFwHZAEfH+SQK4AoYIiq5gMrRaQtcLeIjC9NKSOYqXrBsxgt+BE8i8D3Z+AZN7jbQuSFiPs4cB8HzsZV4otVJAJcjYBG/2wrYT9VXyChbAPfdvBtQ31/g3cN5L6D4g0cHIO62oH7mMB77QDO+uWW9AzDKBtbKpjF+h9/HfB/qvt+Jv9HD2BuIFns8z0wCmgCbKrQIG2g6oHCeWjh91AwEzQDJArcXZDIcyCsC7jbW1/MQUzECc461o3O1rbAc6qeQFXbKrRoFXhXQd77KG9bOziSUXcHxN3BSpbuDogjxpb3YRjVjV0tkv2BpsAbh9inDrD1gG27ij33n4QhIkOBoQCNGjU68OkqSdUHnnlo/udQONNqf5BYCD8FiTgNwk8K+gRRFiJh4G4H7nYIFwOB0lbROvAuRz2/gfc3tHBm4Agn6j4OwnogYd0h7HhEwu17A4YRwuxKGDcAS1R1+WH2O7DaSQ6y3dqoOhmYDNClS5cqXWWlRRvQ/M8g/3Pw7wZJgIgzkIgBENbd+uI0ABBxg/sYq2oq6goA1J8J3hWoZzF4FkLuJDT3FSAMDeuMhPeGsJPA1dJUYRlGOan0hCEitYBzgVsPs+tOrJJEcbUC97sIQqpeyP8SzZ8K3hWAC8L7IJHnQfjJJkmUgTjiIbwXEt4LAPXngGcJ6lloldiynwKeAkcdNPykQALpaaqvDOMo2FHCuBooBD44zH4LgKdEJEJVCwLb+gPbgb8qLLoKoJoPeR+juW+Cfwe4WiGxwyHibMSZbHd4IUEcMRBxMhJxMgDq2wGFc9HCuVDwLZr/MeBCw7og4X0gvA84m5vSh2GUgVRmZ6NAY/daYLaq3nDAc+OAbqraL/A4PrDvLGA00AqrW+3jqnrYbrVdunRRu8dhqBZA7rto3lvgTwN3ZyTmJgjrbb6oKpGq12r/KJwFhXOgaK31hLMhhJ+KRJwK7k5WY7xhVHMislRVu5T0XGWXMPoCLYErS3iuLtB83wNVzRSR/sDLQAqQjjX+YnzFh3l0VH2Q/xma86I1kC7sJCTmZiSsxH8Do4KJuCGsKxLWFWLvC5Q+ZqOFMyDv/9C8t8GRhIafgkScEWhDMiPUDeNAlVrCqEx2lTC0cB6aPc7q1ePugMTeb31RGVWS+nOgcA5a+BMU/hzopZYIEQMCyaOrKXkY1UpVKmGELPVtR7PGQuEP1gjlhJcgfICpeqrixBEDkWcgkWdYVYiFc9CCb6HgCzT/A3DUQiPOssbBuNqaf0+jWjMJ4yipeiD3LTTnFQAk5i6Ivs70eApCIhEQcRoScZrVUaFgJlrwJeQF2qFcLSDiXIg8F3Ee2IEvtFx99dXs3buXr776yu5Q9quKMVU39s8nEcTUsxRNPR/NGW8NsKv5rdVWYZJF0BOJRCLPxJH4GlJrHhL3OEgcmvMcuqcP/rQhaP7n/NOBzzjllFMYPHiw3WH8R1WN63DmzJnDOeecQ/361nQ4U6ZMOei+N910E3fddRfjxo2ja9euxMXFUbNmTc4++2xWrlxZbjGZhHEEVD34M0egaZeBPxdJeA1H4svWDKxGyBFHIhJ1GY7kD5AaP0L0reD7G828D93dE3/mSNS7yu4wbffrr7/SuXNnu8P4j6oa1+Hk5ORw7LHH8uKLLxIZGXnQ/VSVL7/8knPPPZdZs2Zxyy23MH/+fGbOnInL5eLUU08lLS2tXGIyVVJlpP4MNP1W8C6BqKuRmGGII9rusIxKIq7GSOwdaMxt1kDB/I8h/xM0/33U1Q6JGmSNr3HE2h1quenbty/t2rUjISGByZMn43A4GDx4ME8//TQOh/Wb888//yQjI6NSv5iralzl5YwzzuCMM84ArOq4g1myZAkFBQX06tWL77///l/Pvffee8THxzNv3jzOPvvso47JlDDKQIu2oKmXgHc5Ev8cjrjhJllUUyIOJPwEHAnPBqqsRgB+NGskuqcX/syHUO8Ku8MsN1OnTsXlcjF//nwmTpzICy+8wIcffrj/+aVLl+JwOOjYseNhzzV27FhiYmIOeZs7d26lxxWspk+fzplnnonL9d/f/9nZ2fj9fhITE8vltUwJo5TUs8QqWaBI0hTTVdbYTxxxEHUFRF4ORSvQvI+g4Cs0/xPUdSwSdTlEnhXUk0i2a9eOJ554AoBWrVrx+uuvM2PGDC677DLA+mJu1aoVMTGHn3rlpptuYtCgQYfcp3790lXvlmdcwerzzz9n1KhRJT43bNgwOnbsSI8ePcrltUzCKAXN/xzNHA7OBkjiJMTVxO6QjCpIRMB9HBJ/HBr7AORPR/P/h2YNh+yn0aiLkajLg7Kt67jjjvvX43r16rF79+79j5cuXVrqap+kpCSSkpKqXFwV4ZFHHmHMmDGH3Ofnn3+mb9++R3T+DRs2sHHjRgYMGPCf5+6++25++eUXfvnlF5zO8hlLZKqkDkPzpqKZ90FYJyT5I5MsjFIRRywSfRWS/DWS+B6EdYPcN9E9/fCn34Z6UgimQbNut/tfj0UEv9+//3FZGpbLs0rqaOL666+/6NChA0OGDKFVq1bcfPPNTJ8+nRNOOIFjjjmG9evXA3DuuefSuXNnjjnmGKZOnQrA/Pnz6dKlC16vl9TUVFq3bs327dv/8xp33nkna9asOeStW7dupXqvJZk+fTr9+vUjOvrfVeN33XUX//vf/5g5cybNmjU74vMfyJQwDkJVIfcNNOcZa22KhBfNOgtGmYkIhJ+AhJ9gDe7Mmwp5H6GFP4DrWIi+BiIGWtOXBKlNmzaRlpZGp06dSrV/eVZJHW1ca9as4aOPPqJFixYce+yxxMTEsGjRIiZOnMjEiRN58cUXefvtt0lKSiI7O5sTTjiBiy66iBNPPJG+ffsyduxY1qxZw4MPPki9evX+c/4aNWpQo0aNo34vB/P5558zZMiQf20bNmwYH3zwAbNmzaJNmzbl+nomYZRA1Ydmj4a8qRBxJhL/dFD/hzaqBnHWQ2LvQ6NvhYLP0dx30Mx7IPs5iB4CkRcH5fTrS5cuBaxf+8X7/DudTtq2bfuf/cuzSupo42rdujWtW7cGoG3btpx66qmAVdX1008/AfDCCy/w+eefA7Blyxa2bNlCy5YtGT16NB07dqRJkyZcc8015R5/Tk4OGzZsAMDv97NlyxaWL19OUlISjRo1Ys+ePSxcuJBp06btP+bWW2/lvffeY/r06SQmJrJz506A/SW3o2USxgFUPWjGMCicAdHXIzH3Von1so3QIY4oiLoMIi+xJkHMfdOafyxnIhp1KRI1GHHWtjvMUtv3xdyzZ89/bT/22GNZscK+nmKliSs8/J9aA4fDsf+xw+GgqKiIn3/+mblz57Jw4UIiIyPp0qULhYWFAOzevRuPx0NqaipFRUUl9lI6GikpKZx88sn7H48YMYIRI0YwZMgQpkyZwpdffknXrl2pXfufz8orr1gzTvTr1+9f5xoxYgQjR448+qBUNSRvnTt31iPhyxytvh0t1Z/z7hEdbxhHwl+4XH3pd6hvR2v17WinvowH1e/dYHdYIW3Tpk1a/Hviwgsv1J9//llVVefOnatnnnmmTp8+XS+66CJVVf3111/V5XLpihUrVFX1tNNO0+nTp+udd96po0aNqvT4zznnHH3qqafK/bxAih7ke9X8dC5G87+CvHcgaggSfZXd4RjViIR1wJHwojWSPOoSyP8a3XuG1UDu/d3u8KqtgQMHkp2dTceOHXn66af3N6C/+eabJCUlce655zJ27Fg++OCDSi9N9ezZc3/34cpipjcP0KINaOpF4GqDJL1n2iwMW6k/Dc19D/LeA82y1uiIvhHCTjQz5hoV6lDTm5sSBqD+XGtQnkQiCS+YZGHYThxJOGKHITVnIbEPQNFGNP0aNO1itHBuUHXJNUKHSRiAZo8F319I/PMhP221EVzEEYNEX4fUnInEjQLfXjT9OjTtMrRwvkkcRqWq9glDC2ZA/sdWj6jw7naHYxglEglDoi5Bav5gTbXu246mX42mXYl6FtsdnlFNVOuEob69aNbDVrtFzDC7wzGMw7ISx2VIzR+R2EfBtxlNuxJ/2jWmcdyocNU2YagqmvUI+HOQ+OfMokdGUBEJt6YeqfkTEvsgeFejqRfhT78F9a6zOzwjRFXbhEHhT1A4E4m9C3G3tDsawzgiIhFI9LVIzRlWKdmzEE09G3/GvWjRFrvDM0JMtUwYqoXWyFpXS4gKvqUbDeNA4ohBYm5Fas6E6Ouh4Ad07+n4s0ahvlS7wzNCRKUmDBGpKyLviMgeESkQkdUi0ucQ+zcRES3hNvCoAsl9C3xbkdiHTRdaI6SIIwFH7H1IzR8h8kLIex/deyqaMwH159gdnhHkKi1hiEgCMA8Q4EygLXA7sPsQh+0zEKhb7DbzSONQ30409zUIPw0JP/FIT2MYVZo4a+OIfwKp8TWEnWQljL2nornvoeqxOzwjSFXm5IP3AztUtXgd0KZSHpuqqjvLIwjNGgvqtwZDGUaIE1czJPEl1PMbmvMsmj3Kmv4m5m6ION2MGjfKpDKrpM4DFonIhyKyW0SWi8htUrpP7KeBY+aJyEVHGoAWzoHC75CYWxBXwyM9jWEEHQnrgCS+iyROBolAM+9E0y5CCxfZHZoRRCozYTQDbgE2AgOAF4EngVsPcUwOcC8wCDgDmAF8KCJXlrSziAwVkRQRSdmzZ8+/nlMtQLOeAGdTiL7uqN+MYQQbEUHC+yLJnyPxTwZGjV+FP/0m06PKKJVKm3xQRDxY0+aeWGzbWOB8Vf3vKisHP88rQC9VPe5Q+x04+aA/+yXInYgkvoOEl8+C6IYRzFQLIPddNPdVUG9g/ZcbEYm0OzTDRlVl8sEdwOoDtq0BGpXxPIuAMg2cUFXI/8haatUkC8MAAmM4YoYiNb6DiIGQ+wq6ZyBa8K2Zo8ooUWUmjHlA6wO2tQI2l/E8HbGST+n5NoJ/NxJ+8uH3NYxqRpy1cSQ8iyRNBUc8mjEMTR+CetfbHZpRxVRmwnge6C4iD4tICxG5GLgDeHnfDiIyTkRmFHs8REQuF5G2ItJaRO7FavOYUKZX9iy07sNM6cIwDkbCuiLJnyKxj4F3DZp6Dv6ssag/2+7QjCqi0hKGqi7B6ik1CFgJjAEeBV4ptltdoPkBhz4CpABLgEuBa1X1+TK9duFCcNQDp+kZZRiHIuJCoq9Ean4PkRdB3jvo3tPQvE9R9dsdnmGzkF9xT9WD7u4OEQNxxI+1OyzDCCrqXYFmjQLvcnB3QuJGIO5S91ExglBVafS2h2c+aA4ScZrdkRhG0BF3eyTpAyRuHBRtQlPPx5812lRTVVMhnzC04AeQGAgz04AYxpEQcSBRFwaqqS6FvPfQvQPQ/C9Nb6pqJqQThmoRFPwE4Seb9S4M4yiJIwFH/Egk+RNw1kUz70HTr0eLttodmlFJQjphUPQHaIbpTmsY5UjcxyJJHyGxj4B3KZp6Jpr7lvUDzQhpoZ0wvCuse3cHe+MwjBAj4kSiByM1voGwE9DsJ9HUQah3rd2hGRUopBOGen8HSQRnA7tDMYyQJM56SMIkJP558G9HUy9AcyaaKdRDVEgnDLwrwN3eTOFsGBVIRJDIM63SRsQANOclNPVC1LvS7tCMchbCCaMIitYjYZ3sDsQwqgVxJOFIGI8kvAL+dDT1IvyZj5klYkNI6CYMfxagEH6q3ZEYRrUiEadaK/1FXQH5H6N7+6M5r6NaaHdoxlEK7YThbASuMk1saxhGORBHPI64R5EaX0FYNzTnGXTv6WjB93aHZhyFylyitXJpLoT3N+0XQaAwv5CcjDzys/PJy84nLysfcQiRMRFExkQQHR9FbFIM7jC33aEaZSSu5kjia2jhPKsnVcbtaNRgJPYBRMy/Z7AJ3YSBImGmO62d0ndnsvG3v9j+5y7SdqSTvjODrPQcstNyyErNJjvVui/ML12PmqjYSOJrxhFfM46EWnEk1LD+jt9/H0tC7QRqNkgioVY8DkfoFqCDjYT3hLDP0OynIW8KWrQOEl5AHEl2h2aUQQgnDMCfa3cE1UZORi6/z1nNupQ/Wbd0I3/+uom0nRn7nxcR4mvEElcjltikGGo1qkGL45sSlxRLXHIsMYnRRMVGEhkbQVRsJKpKfk4B+dkF5GbmkZWaTVZqNpl7s8jYk8XuLXtZl7KRrL1ZFHl9/4nH6XJSo34SDdvUo+mxjWjavjFN2zeiUdv6hEWYUf92EHEhccNRd1s081E09SIrabgPuXimUYWEdsLQDLsjCFmqyoZfN7H4219Z8t2vrFm4Hr/Pj8MhND6mIZ0HdKD5cU1o3rEJDVrVJaFWPC53+X/cVJW87Hyy9maTsSeL9J0Z7N2Wxp6tqez5ey+bV29l+sTv8BZ6AXA4HTRu14AWxzelxfFNadmpGS07NyMiKrzcYzNKJpHng6s5mn47mnoJGn0jEnOLmb4nCITu9OYdInXxL2NwxN5tdyghQ1VZs2g9c6ctZO4nC9m1eQ8ALTs3o9vA4+l8Wocq+eXrK/KxbcNONq3Ywsbf/mLD8k1s+PUv0nakA1ZppGWnphxzYmva925H+95tiUuKtTnq0Kf+LDRrDBR8Bq42SPzTiLuN3WFVe4ea3jy0E8bckTjiHrA7lKCXk5HLD1Nm8eVr37N13Q5cbiedT+vASRd2p9sZnUisFW93iEckbWc661I2smr+WlbN+4O1SzbgKfAiIjTr0JiOfY/h+H7tad+7HVGxkXaHG7K0YAaa9Sj4M5GY2yD6BkRCu/KjKqumCSNCl8x7Hom52e5Qgtbe7WlMe/YLvp78EwV5hbTr0Yozh/bnxHO7EpMQbXd45c5T6GXdkg0s/3kVv81ayar56/AWenG6nLQ5oQVdBxxP19M70uL4pqZBvZypPw3NehwKvgX3cVZpw9XM7rCqpeqbMOa/gURfaXcoQSd9dybvj/mEryf9iM/n55TLe3HhnWfR4vimdodWqQrzC1k1fx2/zljBsp9+Z13KnwAk1o6n+1ld6HFOFzqd2p7wyKpVBRfMNP9rK3FoPhJ7L0RdhYhJzpWp+iaMhR8hkefYHUrQ8BR4mDb+Kz548jMK8wo57eqTuXz4BdRtVtvu0KqE9N2ZpHy3nEXfLGXJt8vJy84nIjqcrqcfT6/zT6D7WZ1N1VU5UN9uNOsRKJwFYSci8c8gzpp2h1VtVOOEMQ2JPMvuUILCb7NX8dx1r7Jj4y56nteVa8deQaM29e0Oq8ryerz8Nms186cvZt70xaTtzMAd7qbb6R3pM6gn3c/uTGR0hN1hBi1VhfyPrEZxR5RVRRXe2+6wqoXqmzDmT0Kih9gdSpXmK/IxdfQnTB09jTrNajPslRvodKrpF18Wfr+f1fPXMvvjBcyZtpC0HemER4bR/ezO9B/cly6ndcDpctodZlDSog1oxl1QtBairkNi7zLdbytYNU0Ykbr4l1E4Yu+1O5QqK3VHOmMue54Vc9bQf0gfbp9wHZExpkrlaPj9flb+8gezPpjH7I8XkJWaTWLtePpd0ZvThvShafvGdocYdFQL0KwnIf99cB2LJIxHXE3sDitkVZmEISJ1gSeBM4BYYCNws6rOPsQx7YGJQDcgDZgEjNLDBN6lY6wunn0rjvgnyyv8kPL32m0MP30MGbuzuOPVG+h/VR+7Qwo5Xo+Xxd/8yo/vzmLhV8vwFflo1aU5A689hVMu60l0fOj1NKtIWvADmvkw4EViH4PI881ccRWgSiQMEUkAlgG/YCWAPUAzYLuqrjnIMXHAOmAO8ATQGpgCjFTV5w71el06JurimRfgSHqzvN5CyPhj8XoePnMcDocw+quHaN21hd0hhbzMvVnMmDqX796ayaYVW4iICqfvJSdy1s0DaN2lud3hBQ317UQz7gXvYog4E4kbhThi7A4rpFSVhDEW6KOqPctwzM3AU0BtVc0PbHsEuBlocKhShpUwLsSR9MZRRh5a1ixaz4OnjSKuRizjvnuEBi3r2h1StaKqrF2yga8n/8SsD+ZRkFdI667NOfvmAfS95ETTRbcUVH2QOxnNeQmcjZHElxGXSbrlpaokjNXAd0B94GRgO/AG8PLBvvhF5F0gWVXPLLatK7AYaKaqmw72el06ROviOffgiH+iHN9FcFu/bCP39XucuORYxs9+nBr1k+0OqVrLzczlx/fm8OWr37NlzTZiE6MZeO0pnH3LAOo2NV2ZD0cLF6KZd4IWWr2oIvrbHVJIqCoJoyDw5/PAR0BHYALwoKpOPMgxPwBbVfXaYtsaAZuBE1V1wQH7DwWGAnQ+LryzGen9j73bUrm50/2ERYYxfvYT1G5s+rVXFarK77NX8/kr3zHvs8WoX+l+dmfOu+10ju9n1qQ/FPXtQDNuA+8KiL4JiRmGiOmRdjSqSsLwACmqemKxbWOB81W17UGO+QH4W1WvK7atMfAX0ENVFx7s9axxGB8ikeeW11sIWr4iH/f1e5z1yzYycfGTNG7bwO6QjIPYszWVr177gW9e/4mMPVk079iESx84jz6DTix14vB6vNVqsSnVQmt0eP40COuJJDxn1tk4CodKGJU55n4HsPqAbWuARoc4ZidQ54BttQL3uw77io6EUoYW2qaO/oQVc9cw7NWhJllUcTUbJHPN6MuYuvlV7nnzFjwFXr587YdSJYsdm3bx6Qtf88w1L/PoOU/yx+L1lRCx/UTCccSPReJGg2cJuvcC1LvC7rBCUmUmjHlYvZyKa4VVvXQwC4CTRKT4kNn+WO0ffx365Rxo7hRCdZxJaW1asZn3x37KKZf34tQrzUjZYBEWEcbAa07mjZXjeeTDw0/Rn5uZy9jLXmD6xG+JiIqgZoNk7j/1CX54Z1bFB1tFSNQgJPl/AGjqpWjO66gW2RxVaKnMhPE80F1EHhaRFiJyMXAH8PK+HURknIjMKHbM+0AeMEVEjhWRC4AHgfGHG4eBszZ45kHhd+X+RoKFz+dj/A2vEZMQxS0vXGN3OMYRcDgcpZo+/snBEyjIK+T+d27j7tdv4o5XbuCsG/uz4MsUvB5vJURaNYi7PVLjUwjvi+Y8g6YOQr1/2B1WyKi0hKGqS4DzgEHASmAM8CjwSrHd6gLNix2TiVWiqAekYCWX54Dxh31BRxI4m6B5H5bPGwhC0577ij8Wb+CWF68lvkac3eEYFWTup4tY8u1yhj59FW27t9y/vUb9ZP5YtL5CVjqsysSRhCRMRBJeBP8ONPUC/Nkvolq6teONg6vUT5Kqfg18fYjnry5h2wrgCOpSBMK6Q8FXqPqr3RTJW/7YxjsjPqTXBSdw8qWlHvpiBKG3H36fUy7vRcdTjsXptHoI+Yp8bFmzlabtG5GTkUtsYvUa3CYiEHE6hHW3JjDMfRn1LofEV/l3DbdRFiH9LSruDqA54DvocI2Q9c3kHwG44+XrTbfMELb0x9/ITsvh1Kv6/KsksW7pRv5et536LeqWmCyqS9ueOBJxJDyLxI0Fz3w0/VZUC+0OK2iFdMLAHeit611rbxw2WPzdco7r047E2gl2h2JUoCJPEeFR4STXS9z/wyA7PYd5ny1i9+Y9nH3LAMCaFLE4ESE7PYdNK7ew8fdD9TsJDRJ1ERI3BjxzA0mj4PAHGf8R2pWbrhaAEy1ai3CG3dFUmr/XbuPvP7Zx1o1m5GuoS6gVj6/Ih8/r27/tx3dmk/LDb5w25GQatamP3+//15KynkIv74/5hMXfLMNT4MXrKaJmg2TuefPmkB5hLlEXAT4061E09XJIfBlxmqlxyiKkE4ZIGOpqZs2lX43MmDoXEaH3xT3sDsWoYA1a1aXF8U15avAEzrqxP6vmr2Xpj7/T6/xuXPrQefv325c00nam88nzX/PdWzN54N3badm5Gd4CD288NJXHznmKYa/ewLG9ShxHGxIk6hJw1EAz70VTL4CECUhYiWPUjBKEdpUUgKsVFK2zO4pK9cuni+hw8jHUqGdGu4a66PhoRn3xIJ37H8d3b/9MXnY+Q5+5ipufv3r/aG+Hw7G/hDH7owXM+Xg+dZvV4q+VW1jwRQq1GtVk+NQ7GXT/ucTXDP3edBLRD0n+GCQWTRuC5v2v2rTpHK2QLmEAiKs1WvA16s+pFtMgZ6Vms3n1VvpdYQbpVSdDnxlMYX4h7nA3DoeDwvxCvnjle1p3bb5/+vqt63cw83+/kLk3m5adm+H3+Zk2/ktWzf+DYa8Opdf53f61gFbqjnQcDgnJdjBxtYDkaWjGXWjWCPD+BnEjTQ+qwwj5hIGrmXXv2wSO9vbGUglW/mINUmp3YiubIzEqW/Gp0X+btZpX7nybAVf3/SdhrN3O2sUbePC92znl8pMAaNW1BRNvf5Mdf+6kcbuGAGz/cycfPfMFaxauoyCvkON6t+OOV64PufmpxBEHiZPRnImBbrd/WFVUroZ2h1ZlVYMqqUDC8Ja4RlPI+fnDecQlx9Kuh0kY1Vm304/nzkk3cuVjF+/ftn7pRuo2r70/WagqHfq0Izczj99n/zPNW8r3v7HgyxTOuWUAd0++idTtadzR42E2r9la6e+jook4ccQOQxImgW8rmnoBWrjI7rCqrNBPGM6m4GqB5ryM+nPtjqZC5ecWMH/6YvoMOjHkfg0aZTfwmpOp2SB5f5dah9NBZEwEuVl5+Hw+RIS1KX+SvjODVsVWXUzflUHtxjU4c2h/OvQ9hgffuwOANQtCty1QIk5Gkj+1GsTTr0Xzv7I7pCop5BOGiNOaxdK/A815we5wKlRBbiGeAi+N2ta3OxSjCtnX4H3CWZ0oyC1k0dfLKMgtZOPvm3n2mpfpdsbxNG7XAFXl77XbqN2kFr4iP5tX/w2A11NEZEwEv878ZwbY7X/upCAvtAbAiasRkvwBuDuimXejOZNNY/gBQr8NA5CwTmjEOZD/ERo7PGRHPkfHRwGQm5lncyRGVdSiY1MuvudsJt7+Jh8/+wXZaTnE14xj2KtDiYgK5/X732PRN8tIrpdIZEwENx1/H70v7kF2Wg6bV2/lqhFW9dbuLXv4+YN5zPpwHqde2YdL7g+dNWfEEQ9Jb6OZ96M5z4Lvb4h7DBFTYodqkjAAxN0OLfgCNBskNLsOusOsf84ij5nS2SjZmUP7039wHxZ9vYxGbeuTVDeR2MQYVsxdw68/r+TUq/pw6QPn4ff7WfBFCi8Pe4trRl/GkCcupXUXa17QhNoJ9DyvGzUbJPPxs1/w0//NZvjUYTRt39jmd1c+RMIgfjzqbAi5k1DfZkh4CTHr64R+ldR+jsC6S/7d9sZRgbyBROEON7+GjIMLiwjjpAu707hdQ6LirG60qsrfa7bRpts/bRnte7cluW4innzP/mTh9/sJC3fT5JiGnDakLy8tGEOL45vy1aQf/zP9SDATceCIvQeJfwo8S61p0ouq35x0B6o+CcMZmPLAd/iF+oKVJ9+avjkswiQMo3T2zW5bp2ktmhzbkB0brf8fDoeDuKRYstJySNuZAVhJxeFw4Cn07l9jIzImknrN6rDo62Xk54Te/EwSeT6S9A74M9G0y1Fv6Db8l0a1qZLaX8Lw7bQ3jgq0r+1iX1uGYZRWrYY16Hdlb167+x3+WLSeBq3qsfSn38nck8XAa08BrIQhIsz6YB7PXvsKNz47mNbdWvDzh/No26NVyFaFSlgXSP4fmjYYTbsKkt5F3AcuHlo9VKMSRl1w1EVznkeL/rY7mgphEoZxNM677XReTnmK7PQcFn69lIioMB7+4C5qNkgG/ult1X9wH+545QamPPoBnzz/Fcf2bMNFd58d0ot0iasZkvR/IGFo2lWod/XhDwpB1aaEIRIGiW9Yxcr0ayHpA8SZbHdY5Wpf56/CfLOymHFkGrSsy2Mf30t2es7+dTSKvEUs/eE3ugzoiNPlRFU568b+7Nq8h63rtnPX5BtDtudhceJqAklTrYSRdiUkvoaEdbM7rEpVfUoYgLhbIomTwLcLzRpudzjlrlG7BsQkRLNiTvUY1W5UnOKLLq2Yu4avJv/IhuV/Af+UNBJrxZOVmh2SbRcHs3+shqM2mnYtWvCT3SFVqmqVMMAakyExt0Dhz6jnV7vDKVdOp5OGbeqx86/Q7QlmVL42J7SkRr0kHjxtFFPHfMLfa7cxY+pcfp25AvUr4VFhdodYqcRZF0l+H9xt0Izb0Nx3q80Av2qXMACIGgyOZDTnJbsjKXc1GiSza/Meu8MwQkhkdATDXh3KIx/exbzPFvHU4Am8etfb+Hx+rhpx8f6eVtWJOBKRxHcgvC+aPRrNuAX1Z9gdVoWrNm0YxYkjCo04G/LeR9WPSOjkzeYdmjB32kJyM3OJjo+2OxwjhHTu34HO/TuwbumfRMdHkVArnui46tvBQhzRkPAq5E1Bs59F954LCc+F9IJMofNNWUbiagZ4wB9a3Wzrt6gDwJ6taTZHYoSqVp2bU79F3UMmi9UL15G+K6PygrKJiCDR11jtGuJG065EC2baHVaFqbSEISIjRUQPuB3021pEmpSwv4rIwHIJyNXGus//vFxOV1XUalwTgE0rttgciVFdeT1eRg8azzVthjF9wrf4inyHPyjIibs9kjwdXG2teahCtOt+ZZcw1gJ1i91Ks6LRwAOOKZf0LWEdIeJ0NGciGkJrZbTu2pyaDZKZMXWO3aEY1ZQ7zM2TPzxK667NeXnYW9zS5QG2bdhhd1gVThwxSMKLgKIZw1ANve7tlZ0wilR1Z7FbaVpnUw84ptz+FSRuBDji0axHy+uUtnM6nfQ8rxvLZ67cP32DYVS2Rm3q8+T3j/LYx/ewd1saD/Qfxe6/99odVoUTVyMk/kkoWhmSSaOyE0YzEdkmIptE5AMRaVaKYz4Vkd0iMk9ELirPYMSRhETfBN7fQ2qOmPa921KY72HDr3/ZHYpRjYkIJ13YnSe/f4Ts9Bzu6TuCjb9vtjusCicR/ZHYR6FwBppxK6qhs25IZSaMRcDVwOnADUAdYL6IHGy4dQ5wLzAIOAOYAXwoIlce7AVEZKiIpIhIyp49pexaGnEG4EALQqcto1VgZtE1C0MnCRrBq2WnZjz1w6N4Crzc0WM4M9+fa3dIFU6ir0LinoDCOWj6Tag/NNaoEbsGnIhIDLAReFJVx5fymFeAXqp63OH27dKli6akpJQqFn/6HVD4AxL/FBIZ/IvBqCp39BjOzr/28Obq54lLirU7JMMgbWc6owaNZ+UvfzDw2lO45cVriIyOsDusCqX5n6KZw8HdHkmchDiS7A7psERkqaqW2DfYtm61qpoDrAJaluGwRWXcv1Qk/kkI62b1bsj7tLxPX+lEhDsn3UjG7kx+mDLL7nAMA4CkOok8M2MElz10Pt+//TO3dn2Qreu22x1WhZLIC5CECeD9A029FC0K7t6Lh0wYIpIlIjUCf2cHHpd4K+sLi0gE0AYoS/eJjmXcv3SxOKKsOabCTkSzHkI9v5X3S1S65h2aEJccy9a1of0f0gguLreLa8dczlM/PkrW3iyG9XyEVfPX2h1WhZKI/kjSFPCno2mXoN7gfb+HK2HcDmQH/r4t8Phgt0MSkWdFpI+INBWRE4BpQDTwTuD5cSIyo9j+Q0TkchFpKyKtReRe4FZgQpneYSmJRFq/BCQBzXm+Il6i0tVrXpvtG0N3wSgjeB1/SntenD+GmMRo7uv3OD++O9vukCqUhHVGkv8HONH064K2pHHIhKGq7+g/TfzZgcf/ugHvYjVgH04D4H9YYzE+BQqB7qq6r9tEXaD5Acc8AqQAS4BLgWtVtcK+zcURg8TcCJ75qGdJRb1MpWlyTENWz1/L6gXB+4vGCF31W9TlpfljOObEVjx99URevWtKSC3zeiBxtUAS3wYtRNOvRX3BN+dbqRu9RSQfmArcoap5gW0NgP8D2qhqaZJGpSlLo3dx6s9Bd3dCYu6xkkcQS92Rzt19HiNzTxYTFo6lYev6dodkGP/hK/Ix6d53+eylbxh47SncNfnG/VOohyL1LEfTh4CzEZL0TpVrCC+vRu8TgO7AchHpIiKXACuAfKDD0YdZNYgjBiQa9Qf/IKPkuok889NjqCqT73/P7nAMo0ROl5NbXriGKx+9iO/emsnzN7wW0tOJSFhHJOEVKPrLWvbVl2p3SKVW6oShqr8DXYBfgAXAe8AIVT1dVUOrotzZBPI/RwsX2h3JUavVqCaD7juXhV8uZev60J+ewQheg0cOspLG2z8z8sJnKMgLnQFvB5LwnkjiZCjaYk1Y6AuONWzKWu7rAPQBNgAeoJuIhFwnf0kYD44kNP2akFgc5cRzrNLlyrmhM2eWEXpEhCGPX8IdL1/Poq+WcVu3B1m/bKPdYVUYCe+BJL0B/h1o+vWoVv2VC0udMETkUWAO8DlW4ugMtAZWiMhJFROePcTVFEme9s/iKNnj7A7pqDRq14BGbevzzogPyUrLPvwBhmGjs28ewLjvHiYnI5fbuw/no2dCZxaGA0lYNyThBSj6A80aZXc4h1WWEsbNwNmqereqelR1LdAD+AAIuYVtrZknX4aoq6wFUvI+sDukI+Z0OnnwvTtI35XJe49/bHc4hnFYnft34PUV4+lxThdef+D/WPBl2TuwBAsJ7wvRN0H+x2j+Z3aHc0hlSRjHqeoPxTeoapGqPgj0L9+wqgYRBxI7HMJ6o1lPBHVX25admnHyZT35Ycos0nam2x2OYRxWbGIMD00dRrMOjRl//ashPSpcYu6AsO5o5qOoZ7Hd4RxUWRq9D9ptSFVDdvEFESeS8Dw4aqM5FTJmsNIMuvcc/D4/9/V7vFqshmYEv7BwNw//7y4A7uk7gs1rttocUcUQcSEJL4GzoTVZYRVdoyd0OzuXI3HEQuQZ4ElB/cHbBtC0fWNGf/UQuzfvZdQl44O+Md+oHhq1qc+zP49EVbnrpEf5fc5qu0OqEOJIQJLeAomxGsGL/rQ7pP8wCaOUJPwUoAjNnWx3KEelQ99juGn8EFbMWcOsD+fbHY5hlErjdg15cd4YEmrF80D/J5j5v1/sDqlCiLOulTTwoamD0MKqVXljEkZpuTtB5IWQOwl/9otB/et84HWnUKN+Egu+DN42GaP6qdusNi/NH0O7E1vz1OAJLPnuV7tDqhDiaoEkfwLOemj6UDT3zSrzfWMSRimJCBI3JpA0Xg7q9gyn00mLTk1J+W55yBbvjdAUkxDNE58/QNP2jXji4uf4bfYqu0OqEOKsjyR9COH90eyn0Jxn7A4JMAmjTEQcxZLGRDTvfbtDOmJDn76K+Jpx3Nfvcb55PeR6RRshLDouirHfDKdWoxo8fMZYlv4Y/MsRlEQcUUjCi4Hvmylo0d92h2QSRllZSWMUhJ9sdbUt+NHukI5Iw9b1mbhoHJ1P68ALN01m4VdL7Q7JMEotqU4iz816nPot6/LERc+xY1NozU60j4gDibkTcKC5k+wOxySMIyHiQuKfB3d7NPPeoJ3bPjo+mkc/upsWxzdh3JUvkpuZa3dIhlFqCTXjeeLzB0Dgyasm4Cn02h1ShRBnbYi6GPI/Rb0rbI3FJIwjZBUXJwBONPMhVINzHv/I6AjunHQjeVn5fPPGTLvDMYwyqd24JndPvonV89fy9JAJIbuehsTcDo5aaPqttq6jYRLGURBnHST2YfAuQXMmVJmeDGXVqnNzOp3anjcfmsqcaQvsDscwyqTPoBO54akrmf3RAsZe/gL5uVV/Er+yEkcSkvgK+DPQjNtR9dgSh0kYRyvyAog4z+o5lf1U0CaNx6bdS5sTWjDmsheY/7npblvepk6dSpMmTXA4HDRp0oSpU6faHVJIufjec7j+ySuZO20hd/QYzs6/gmO68LIQdzsk/knwLkMz70O1qNJjMAnjKImI9Y8YdSXkvYVmPYRq8C3+YvU8eZhWnZsx+tLnWWGmQi83U6dOZejQoWzevBlVZfPmzQwdOtQkjXIkIlxy/7mM/fZhdm/Zy2t3T7E7pAohkWcgsQ9Cwbdo5sOVXhVuEkY5sCYpfBSib7MaprKftDukIxIVG8mYr4dTu3ENRl8ynvTdmXaHFBIefvhh8vLy/rUtLy+Phx9+2KaIQlfn/h045+YBLPgiJWQnK5Toa602jYLP0OzRlfraJmGUExHBEXsHRF0Nee+gue/aHdIRiUuO5dGP7iE7PZcR5z1l1s8oB1u2lNyL7mDbjaNz7m0DiY6PYtQl40N31b7o2wLfNf+H5n9daS9rEkY5k9gHILwfmj0WzQ/OhV+aHdeY4e8PY8OyTdzd+zH2bk+zO6Sg1qhRozJtN45OjfrJPPh/w9j0+xbGXfEiRd7Kr+uvaCKCxN4P7uPQrMdRX+WMQzEJo5xZ06GPh7BuaOb9QbvwUq/zT2Dcd4+we8teHhow2pQ0jsKYMWOIior617aoqCjGjBljU0Shr9vpx3PLi9cw//MljLvypZDsbmuNB3satCDQnlHxHW4qLWGIyEgR0QNuOw9zTHsRmS0i+SKyTUQeExGprJiPlEiktcB7eB8067GgrZ7q0PcYHp9+P9vW7+ChgWPI3Jtld0hB6YorrmDy5Mk0btwYEaFx48ZMnjyZK664wu7QQtp5t53O9U9eyZyPF/D15NCc/kZczZDY+8AzB/Ir/sepVFY3UBEZCVwK9C222aeqJY5CEZE4YB3WOuJPYK0fPgUYqarPHe71unTpoikp9i7rqOpBM+6EwhlI4utIeG9b4zlSC79ayhMXP0e95rV56sfHSK6baHdIhlEqqsoDp41i7eINvLHqeWo2SLY7pHKn6kfTrwfvUiR5OuJqelTnE5GlqtqlpOcqu0qqSFV3FrsdasjiFUAUMERVV6rqJ8BTwN3BUMoAEAlDEp4DV2s04260aLPdIR2R7md1Zty3D7Nr8x4eOWsc+Tn5dodkGKUiItw16Ub8Pj/PD30taMdJHYqIA4kfC4ShGfegWnEDFys7YTQLVC1tEpEPRKTZIfbtAcxV1eLfTt8D9YAmFRlkeRKJRBJeBhxo2pWoJzhn1uzQ9xge/egeNv6+meFnjDVLvBpBo26z2lz35BUs+W450yd8a3c4FUKcdazxYEUr0cwRFZYYKzNhLAKuBk4HbgDqAPNF5GBlxDrAgU3/u4o99x8iMlREUkQkZc8e++ZbOZC4GiJJ74C40bTL0bxpdod0RLqdfjzDpw5j/dKN3NLlAdYu2WB3SIZRKufcMoAe53Rh0r3vsnLeH3aHUyEkot/+8RnkVUy7aaUlDFX9VlU/UtXfVfUn4KzA6w851GEHPJaDbN/3GpNVtYuqdqlZs+bRB12OxN3WWkUrrCuaNRx/1tigLB73GXQiL8wbjcvt5N5TRrJshr2zZxpGaTgcDh545zZq1E/ilWFvBeX/vVKJvhXCT0Wzn0S968v99LZ1q1XVHGAV0PIgu+zkvyWJWoH7oJz8XhyJSOIbEDUY8qagWSODcpbbFh2b8sK8MdRtVptHzhzLj+/Ntjskwzis6PhoBo8cxPplm5jzcWhOsmm1Z4wBiUJzni7389uWMEQkAmgD7DjILguAkwL77dMf2A78VbHRVRwRlzXDbfRQyP+f1e02COeeSq6byLM/j+SYnq15eshEJt37bkj2dTdCS78rTqJ5xya8cteUkF3/RRyJSMwtUDgbLZxXrueuzHEYz4pIHxFpKiInANOAaOCdwPPjRGRGsUPeB/KAKSJyrIhcADwIjNcgL0+KCBJzD0TfAvkfoRk3of7gGxgXlxTLuO8e4dxbBzJt/Jc8d/2r+HzBl/yM6sPpcnLXpBtJ35nB+2M+tTucihN1FTgbWquClmOvqcosYTQA/gesBT4FCoHuqrqvr2ldoPm+nVU1E6tEUQ9IAV4GngPGV2LMFcaae+pOJO4JKJyHpg4Kym63LreL2yZcx+CRg/hhyixGXvAMGXvMpIVG1dW6awv6XXkS0yd+y56tqXaHUyFEwqzvFt8mNGdC+Z03yH+sH1RVGLhXWlq4CM24HQBJmoK429kc0ZGZPuFbJt/3LjGJ0dz71q10O/14u0MyjBLt2LSL64+5i86ndeDxz+4nSIZ2lZk/8xHIn4Ykf4S4jyvVMVVp4J5RAgk/AUmeBhKJpl+LFv1pd0hH5LzbT+flJU+SWDuBR84ax9TRn5h2DaNKqtu0NlePuowFX6Sw+JtldodTYST2AXDUQLNGl0vPMJMwqghxNUKSpgCCpl2NFgXnGIem7Rvz0oIxnHJFL6Y89gEjL3jGzEFlVEkXDDuDxNrxfD/lZ7tDqTDiiEVi7gTvcig4+kGLJmFUIeJqiiS+DepB956P5r4blN1uwyPDeeCd27n5+atJ+W45Qzvcy7Kffrc7LMP4F6fLycmX9mLBFyns3RaabRkARJ5vTU+U/Qzqzzv8/odgEkYVI+42SI2vILw7mj3aqqLyHXJS3ypJRLhg2Jm8tHAsMQlRPDRwdMhOy2AEr/PuOB2/X5k2/iu7Q6kwIk4k7jHwbz/qsRkmYVRB4qyJJEy2ejl4f0X3not6gqMB/0AtOjZl4qJxnHBWZ14e9hZPDn7J9KIyqoy6TWvT64IT+PHd2Xg9XrvDqTAS1jWwQt/7aOHcIz6PSRhVlIggUZciyZ+BIx5NGxK0K/hFxkQy4pN7ueKRC5n1wXyubTOMb96YEbrTMxhB5bTBfchKzWb+9CV2h1KhJPZucLVAMx9F9chWITQJo4oTVzMk+SMI64xm3oc/+xlUg++XkNPp5OonLmXS8mdo0r4Rzw99jYfPGkfaznS7QzOquS4DO9KwdT2mjgntXn0i4VYDuH87HGEpwySMICCOBGsOqshLIfd1NO3SoBzkB9C4XUOe+/lxbp94Pb/9vJKhx93DjKlzTWnDsI3T6eTie89h04otbFqxxe5wKlb4yVY32/yPjuhwkzCChEgYjvgnkIQJULQFTT0XzQvOqQ1EhHNuGcArS5+mbrPaPHnVSzw4cDQ7NgblnJJGCGhzgjUHaqgnDBE3RF4IhbPQor/KfLxJGEFGIgYgNb4Ad3s068FAFVVwFqMbt23AC/NGc/vE6/lj0Xpu7HgvP/3fHLvDMqqhRm3qEx0fxW8/r7Q7lAonUYMBN5rzcpmPNQkjCImzLpI4BSIvt6qoMh8IynYNsKoDzrllAK+vGE+L45vy1OAJPDHoObb/GXxdiY3g5XQ56TKgAyk/BOeKmGUhzpoQfSUUfFnmWSVMwghSVt/qEVYjVsHnVi+qIB0dDlCrYQ2emTGCq5+4lMVfL+Patncy8fY3TRdco9K07tKCvdvSyEoNvpmjy0qirwfcaG7ZVuYzCSOIWdOk34LEPw1F69G951hVVP7gnOff6XJyxSMXMmX9BAZeewpfvvYD17YZxleTfgzp3itG1dCyczMAVv4Smku4FieOJIgYCAVflGn0t0kYIUAiz0Nqfg+R51pVVHvPQAtmHP7AKqpGvSTufG0ok397lmYdmvDizZO5o8dwVsxdY3doRgg7tlcbYhKiWfBlcA6SLSuJugQ0FwpKP8rdJIwQIY4kHPHjkKT/gSMOzbgZf+bIcl08pbI1bteQZ2aM4IF3b2fvtjTu7vMYI85/mi1/bLM7NCMEudwuEuskkJedb3colcPdGVzHoLmvoeop1SEmYYQYCeuMJH8CUddB/vto6kWod53dYR0xEeHUK3szZd0Erhl9GctnruTGDvcw5dEPKMwvtDs8I8QU5BYQER1udxiVQkSQ2HvBtxXy3i/VMSZhhCCRMBxxDyCJb4I/FU29AM15udS/IqqiiKhwLh9+AVPWT6DvpT2ZOuYThh53Dz++N5si75FNc2AYB4qKjSS/upQwAAnvCWEnojmvlKqnpUkYIUzCT0KSv4SIfmjOi+jec1DPYrvDOiqJteJ54J3beerHxwiPCufpIRO5utUdTJ/4rSlxGEfN6XLi9VSvHyASeSFoBhRtPOy+JmGEOHHWwJHwojW1CB407Ur8mcOPel58u3Xq155Jy59l1BcPklw/iZfveIshLW/ny9d+COlZR42KlZ2WQ1xyrN1hVC53W+u+aNVhdzUJo5qQ8N5Ija8h+gbI/wRNuzhol4LdR0ToflZnXvxlNM/OHEmdprV46ZbXubbtnXzz+k94Ck3iMMrGFebC5/XZHUblcjYFiUY9yw+7q0kY1YhIJI7Y+5DEt8C312oQz/8qJCb+69D3GJ6fM4oxXw8nLjmW52+cxODmtzJt/Jfk5wZvTzGj8uz7f1CQV72qNkWcENYLCmcedpohkzCqIQnvidT43Fq2MfNuNO0q1POr3WEdNRGh2+nHM3HROJ78/hEatKrHpHvfZXDz25g+8VtTVWUc0p/L/2LHxl107t/B7lAqnUScBv7d4D301CgmYVRT4qyDJL2HxD4Kvo1o2iX4029CvWvtDu2oiQid+3fg2ZkjeX7uKBq1rc/Ld7zFtW2G8dlL35CTEZwj4Y2KNW/6YhwOoc/FPewOpfKF9wUch12Nz7aEISLDRURFZOIh9mkS2OfA28DKjDVUibiR6KuQGj8hMXeBZwmaeg7+zIdQf5rd4ZWLY3u24dmZIxn77cMk1E7glTvf5tL6Q3n22lf4Y/F6u8MzqpBlM1bQsnOz6tfoDYgjFpxNoOjQsynYkjBEpDtwA/B7KQ8ZCNQtdptZQaFVS+KIQmJuRmrOgKhrIf9zdM9ANO+TkGjfEBG6DujIhAVjeSXlKU69sjezP57P7d2Hc3uP4fz8wTwzlqOa+3vtNlbPX8uJ53azOxT7uNuBd/Uhd6n0hCEi8cBU4DqgtOtzpqrqzmK34B2BVoWJI8Ea8Jc8HVzN0KyH0LQrg3qk+IFadmrGnZNu5INtk7n1pWvJTsth7OUvcFWzW/m/UdPYuz00SlZG2fwwZRZOl5OB155sdyi2EXcn8O845D52lDAmA9NUtSylhE9FZLeIzBORiw62k4gMFZEUEUnZs2fP0UdaTYm7FZL0PhI3GorWoaln4898EPVttzu0chMdF8V5t53OW2teYPSXD9KoXQPeGfEhVzS+mZEXPM2S7341M+RWE6rKwq+X0v6kNiTVSbQ7HPtEDASch9ylUhOGiNwAtAAeLeUhOcC9wCDgDGAG8KGIXFnSzqo6WVW7qGqXmjVrlkfI1ZaIA4kahNT8MVBN9RW65zT8WU+h/gy7wys3DoeDE87szFPfP8qUdS9x0d1ns2reWoafMZZrAo3kuVnBPcjROLTPXvyGv1b+TZ9BPe0OxVbirAFhJx56n8qqoxaR1sAvwEmq+kdg2yxgpareVobzvAL0UtXjDrVfly5dNCWlekxTXBnUtx3NfhEKpoPEIXHDIeI8RMTu0Mqdp9DLL58u4vOJ37J6wTqiYiPpd8VJ9Dy/G8f1aYc7zG13iEY5+XXmCh48bRQ9zu3KYx/fg8NRvTuO+rNG4ox/fKmqdinp+cpMGFcDbwPFh1E6AQX8QLSqHnbEjIgMAV5T1chD7WcSRsVQ7x9o1kjwLoOwk5D4UYiznt1hVZi1Szbw2Uvf8MuniyjM9xAdH8VJF3ZnwNV9OaZnm5BMmNXF3u1p3NzpfuKSY5iwcBxRsYf8SqkW/Bl340x8vkokjASgwQGb3wbWA2OBVVqKYETkeeBcVW12qP1Mwqg4qj7Im4rmjAcEibkVIi9FHDF2h1ZhCvIKWfbT7/zy6SLmfrKQgtxC6jarTb8rTuLky3rRqE19u0M0yqAgr5D7+j3OXyu2MHHxOBq3a2h3SFWCP30ozqTX7U8YJb74AVVSIjIO6Kaq/QKPhwBe4FesUsjZWMnlAVV9/lDnNgmj4mnRVjRrBHjmgsRA1GVI1FWIs47doVWo/Jx85n6yiJ/+bw7LZ65EVWl2XGN6X9SD3hd3p2FrkzyqMk+hl1EXP8fib5bx6Mf30Ov8E+wOqcrwZwzDmfhS0CSMKUBfVW0SeDwEeABojFWVtQ54QVX/73DnNgmj8qh3BZr7BhR8Dzgh4iwk+nrE3dLu0Crc3m2pzP1kEbM+ms/q+dYo+abtGwWSRw9T8qhiMvZk8viFz7Lylz+44+XrOfvmAXaHVKX4M0fgTHiiaiaMimQSRuXToi1o3juQPw20wEocMcMQVyO7Q6sUe7am8suni5gzbQGr5q1FVWnQqi4nnNmZ7md15thebXC5XXaHWW2tXbKB0ZeMJ21nBve9fSt9L6nevaJK4s9+GmfcAyZhGJVH/elo7puQ+y5QBJEXITG3hHxVVXF7t6Uyb/oSFn29lOUzV+L1FBEVF0mXAR3pdvrxdB3YsXr3+a9Efr+fT57/mjcfmkpyvUQe+/geWndtYXdYVZI/8xGcCWNMwjAqn/p2o7mvQt5HgAOiLkGir0Gc1auaJj8nn19nrGTBlyks/mYZaTszAGjZqSldBnSk68Djadu9pSl9VICV8/5g0r3v8sei9fQ8vxv3vHEzsYmh2znjaPnTb8OZ9LJJGIZ9tOhvNPcVyP8cUIg4E4m+AXG3tju0Suf3+9n422YWf/srKd8vZ9X8tfh9fiKiwzm2Vxs69D2WDn2PoWWnpiaBHIW/Vv3NOyM+5JdPF5FcL5Frx1xO/8F9TDfow/Cn34EzaYJJGIb91LcDzX0b8j8CzYOw3lZVVVgnu0OzTW5mLstmrGT5zBX8NmsVm1dvBSAiOpx2J7amfa+2HNurDa27tSAyOsLmaKs2VeXXmSuZ9twXLPluORHR4Vxy/3lcePdZ5tqVkmn0Nqoc9WdA3vto3rvgT4PwfkjsvYirud2h2S59Vwa/z1nD77NX8fuc1WxetRVVxeF00OL4prTp1oI2J7SkbfdW1G9Rx/xixur5NOuD+Xz71gw2/raZxNrxnHvr6Zx1U3/ia8TZHV5Q8WeOxJlQBUZ6VzaTMKo+9edB3jto7mSrV1XkRUjM7Yizlt2hVRnZ6TmsXrCOVfP+YM3Cdaxd8if5OdaSs9HxUbQ4viktOzWjZaemNOvQhAat6laLqqy929NY8EUK86YvZvnMlfiKfDTv2ITzbjudUy7vRVhEmN0hBh1VD7qnN87ai0zCMKou9aehOS9D3v8AF0ScjkRdBO4u5hf0AXw+H1tWb2XNog1sWLaR9cs28udvm/EWWsvPusNcNGrXgGbHNaZR2wY0alufRm3qU6dpraBOJHu3pbJq3lp+mxUoeQWq7uq1qMNJF5xAvyt70/TY6tF9u6JowQ9oxm046643CcOo+rRoC5r7OhR8BZoLzqZI5EUQeb41k6ZRoiJvEX//sY1NK7bw52+b2bRiMxt/30zq9n+Wm3E4HdRpWov6LetSv3kd6javTZ2mtajbtBY1G9YgOj6qSiTnIm8ROzbuYuu6Hfy5/C/WpmxgXcpG0nZY7yUyJoJje7XhuD7H0OPszjRq26BKxB0K/JkPQsEMnHVSTMIwgof686DgOzT/I2uSQ1wQ3heJugLCTjRfEKWUm5nL5jXb2LZuB1vXbWfr+h1sW7+D7Rt27q/W2iciKpzk+kkk10skqW4iyXUSSKqbSFyNOKLjIomKiyQiOoKI6PDA7Z+/nc6S11Dw+/0U5nvIz84nLyufvOx88rMLAvf5ZKfnkp2WQ+r2NHZs2s2OP3ey8689+H3/rEPSsHU9WnVtTusuLWjbvSUtOzXD6Tr0mg1G2akquqcvuI/DmTTRJAwjOGnRn2jeNCj4zGogd7VEogZD5DmImNlFj4Sqkrk3i52bdrNz0272bE0jdVsqe7enkbo9nbQd6aTtyKAg77CTRwPgdDlxuZ24wlw4nA6KPEV4C70UeX2HPxiITYymbvM61Gtem7rNatOwdX3qt6pL47b1iY6PPpq3apSSetegqecicY/jiL7cJAwjuKkWQv7X1tQjRWtAEiBqEBJ1BeKsa3d4IUdVycvOJzsth7ysfHIz8yjILaAgt9C65QXucwus5OApwuspwlfkJyzchTsijLBwN2GRYUTFRhAVF0VkTASRsRFExUYSGRtJbFIMMQlRZn2RKsCf+Qjkf47UmoPDmXTQhBG8rWBGtSISDlEXQOT54E1Bc9+B3DfQ3DfQ8N5I5MVWtZWYL5/yICJEx0URHRdldyhGBVN/JuR/YZXaHYeersYkDCOoiAiEdUXCulrTq+d/DPmfoIW3gqMGGnk+Enkx4mpid6iGERwKZwAFSNQlh921eq9HaAQ1cTXAEXsXUnMWkjAJ3B0g9y1072n4066yuglq6erRDaO60sLZ4KgFrvaH3deUMIygJ+KCiJORiJNR327I/wzN/xDNuA2cDSBqsDUoMIRXBDSMI6FFW6BwjjX2qRS9D00Jwwgp4qyFxNyI1PgBSZgAjtpo9lh0z0n4s8ag3tWEakcPwygL9Weg6TcAbiT6xlIdY0oYRkiySh0DkIgBgRUB37HWIc97Bxy1rIby8D4Q1tOUPIxqR9WDpt8Gvq1I0juIq3GpjjMJwwh54m6PJDyL+h6EwtmoZw4UfI/mTwMi0IhTkcjzAoMCzX8JI/Rp1ljwLkbin0PCSuxBWyLzv8OoNsRZA6IuRKIuRNUL3uVo/tdQ8DVa8BU4aqIRZ1rJw9XWjCg3QpLmfQz570P0DUjk2WU61iQMo1oScf/TPTduOBTOQvOnB6qtpoCrFUScblVrucxynkZoUM9iNGukVRUbc3eZjzcJw6j2RMIg4jQk4jTUnw4F36L5X6A5L0LOi6iz2f72EFPyMIKVeleg6TeCqzGSMB6Rss/JZRKGYRQjjkSIuhyJuhz17YKCH9HCHyB3krU+ubMBGt4fiTgN3McjYjoaGlWfetehadeBIxFJfOuwI7oPxrZPu4gMFxEVkYmH2a+9iMwWkXwR2SYij4n5iWdUAnHWRqKvxJH0LlJrPhI3FlwtIO//0LTL0D298WeNQj1LzABBo8rSgu/RtEtAwpDEKYizzhGfy5YShoh0B24Afj/MfnHAj8AcoCvQGpgC5ALPVWyUhvEPcSRB1EVI1EWoPwcKf0YLvoO8j9C890Di0fCeSHhva63yI1i/Q9ULFJlZeI1yoVqEZj8HeW+CuwOSMOGokgXYkDBEJB6YClwHPHaY3a8AooAhqpoPrBSRtsDdIjJezQgswwbiiIHIs5HIs1F/LnjmWNMrFM5BC74BQF3HQHgfK4G4Oxy2vlh9e62lagtno5oF4f2R2LsRR0IlvCMj1KhvG5pxP3iXWFWsscOttrqjZEeV1GRgmqrOLMW+PYC5gWSxz/dAPaDJgTuLyFARSRGRlD179pRLsIZxKOKIRiJOxxH/JFLzFyT5MyTmLpBwyH0NTbvCWj3wENSfgWaPg/zpSOwDSNxIKFqF5r5ROW/CCBmqiuZ9iu49C4pWI/HP4IgbWS7JAiq5hCEiNwAtgKtKeUgdYOsB23YVe25T8SdUdTJWQqJLly6m9GFUKhEHuI8B9zFIzM3WtNHelYgj7tAHFs4ETwqS8DQS3tfa5tuD5kxAo2866Eh09a4CzxJwNQV3Z8QRg6qaXlzVlPr2WF1mC38Ed1ck/inE1aBcX6PSEoaItAbGAiepqqcMhx74xS8H2W4YVYo44iG852H304LvwN0O3MVG3LqagjPZSggRJ/8rEagWQO7b1tTuzqaQ9x5IJMQ+hJTi9YzQov4MqzSa9x5oERL7AERdfUTdZg+nMksYPYAaWO0Q+7Y5gd4ichMQraoHrgm5E6skUVytwP0uDCPIqSp4VyHRN/y7JKH5gMD+qgQNJA1HoMH9R4g4F0fsMACrt1bua+BuhzgSrQb03LfRvPchrCMS1hPCTkBcjSr9PRoVQ/05kDcFzX3LqvaMOBOJuR1xNa2w16zMNozpQHugY7FbCvBB4O+SSh0LgJNEJKLYtv7AduCvCorTMCqP5oJ/L7ha/nu7fzf4c8HVHAhUd+0rXBf9CY7Yfy14I5EXgD8VvIGOh+pFPfPAEQuOJDRvCrq3P/6MYWa23iCnvt1WdeWeU9Ccl6wfAslf4EgYX6HJAiqxhKGqGUBG8W0ikgukqerKwONxQDdV7RfY5X1gBDBFREYDrYAHgcdNDykjJPi2gsTzT00rqD8PLdoAzhr/6ga5v2TubAreN0Bz/jmPPwOKtgCBEolmQdE6JP4ZJLyXtUm94M/6VxuHFm2Bgi9Q72prErqIsxFnzQp6s8aRskqiS9G8/4OCH4Aia0nimNsQ93GVFkdVG+ldF2i+74GqZopIf+BlrNJIOtb4i/H2hGcY5UyiwFkP9SxBwk+0tvk2g2cZhPcG9o3PcP4zqjy8L7g7oRn3WWucO2uj2c8APtjXDde3C/wZaM7LqGep1ZjuPg5xJu9vD1HvOjTjJnDUBXdray4t71qIe/jwDfVGpVB/LhR8YyWKojUgcRB1FRJ1eamnJC9PtiYMVe17wOOrS9hnBdC7kkIyjEqjqoirERrWAzxzUU8vwI9mPWFN4RB5MRCYKLEYcURD/Cg070M0/zNwJENYFyj0/pMwJCLQvdeFehahnjlIzD0QfuI/jec5z4CjHpLwHOKsjXrXo6nnQMQpEDGgEq+EUZz6c63JMAu+g8JZQCG4WiNxo6wSoCPKttiqWgnDMKqNfV/cEn0Nmp2Kpl8HEgPhvZCowYirIZr/FVo426p6KPaLUpz1kdi7IfZuVH1WScHdCnHWtZ53t0ZdLayeMpGD0KwRaObDUGM64oi3qrw8vyFxjyDO2oFjWqLhvdDCX6yJFo1Ko75t1qDNwtlQuAAoAEdNiLoYiTgT3J2qRHdpkzAMw2birIUkPA08bU146EjcP9BKvSug4AuItaaiVi2Agp+shk5nTau9I/cNKPoTiX3o3+cVJ6r+wMj0Qah3NXgWWqWHwnlWV1z3sfv3V38mSAL4dwVey28mV6wg6tsF3l9RzzLwzIOi9dYTzgbWFDQRA62xNRXQNfZomIRhGFXIvl/7+x/HPmD9ygyUHNACNO9tyHkJdbUAFDwLkbjHIfxkaxd/mjX3lfUocJcDFFkj0AH1LLUmUty/H6CF4NsEB2lEVfWAPw0ctUwiKQPVfPCus0bve5aBd5nV2QGAcAg7Hom9EML7gLNZlShJHIxJGIZRhYk4rC/2fY8dCZA4BTyLrIZyRxzEPrx/RK+qD819D9xtrbmsJNyaVyjvfZBYcLW3TqRZ4KzHv74C/Lug6G8k+vp9r/bvYIrWoakXAOGoswG4GoKzEeJsaJ3LWd+6l/gq/aVXUVTVSqhFf0DRH1aJrmgNFG0E/NZOjppW9VLUVeDuBO625TZtR2UwCcMwgow4YiHiVCTi1BKe9SAShmY9DgjqahkoOfyNxI8N9JLyg6s1eFf8e7CgZ7E1UNDd1XqdA7/0HbWRuJFo0WbrF7Lvb/AsQQ+cK0uiUWddcNSybs6aiKMWOGuBo8Y/N4kJusSimg++bYHbdtS31erO7Avcil8LRx1rBH/4AMTdFlztwFk/6N5zcSZhGEYIEYmEmJsh+gbw/mYN5JNYCD8ZcSYH9nFAWHc0/0O0cD4SfiJaONcqhUQM2L/ff87trGnNfFpsm6qCpu//AsW3DfVtt/7277GSkH8PireEM4ahjhrWFCiOZCtOiQFHDCIxVpdjCbfaWiQicF/sb/ZtCwNxYX2duf49zkQV8IF6AI91rwWgeYFb7j/3fmub7n8uC/yZ4M+y/vbttt7rv7itdgdnQwjrjDgbWYMw3W2LVQuGDpMwDCMEibggrLN1K0l4T4gajGbcbpUGEHC1QWLuLOPrCEiS1Rbitqq7Dvz9rOoHzQDfHmtUe+Cm/lTw7bVGqPt2g2602lr8OWiJEz+UjuLkn682D2Wfds4RSESx4Ii3xj44G4C7I+Ks/0/Vm7M+OGpWuYbpiiShOmBaRPYAm2146RrAXhtetyox18BcAzDXAILzGjRW1RKH+4dswrCLiKSoapfD7xm6zDUw1wDMNYDQuwamb5xhGIZRKiZhGIZhGKViEkb5m2x3AFWAuQbmGoC5BhBi18C0YRiGYRilYkoYhmEYRqmYhGEYhmGUikkYhmEYRqmYhHEIInKLiGwSkQIRWSoiJx1i3wgRmSIiv4uIV0RmHWS/PoFzFYjIRhG5qcLeQDko72sgIn1FREu4tanQN3IUyngN+orI5yKyQ0TyAtfi2hL2C6rPAZT/dagGn4V2IvKziOwq9u88Vg6YbTCoPguqam4l3IBLAC9wA9AWmADkAI0Osn808BowFJgOzCphn6ZAbuBcbQPn9gIX2v1+K/Ea9MWaq6EdUKfYzWn3+y2nazAcGA30BJoBNwNFwOXB+jmowOsQ6p+FFsDVQAegMXAOsAt4Olg/C7YHUFVvwCLg9QO2rQfGleLYiQf5snwKWH/AtjeABXa/30q8Bvu+JGrY/f4q+hoU2/8j4JNg/RxU4HWojp+F8cX/nYPts2CqpEoQKDJ2Bn444KkfgBOP4tQ9Sjjn90AXOXDhZptV4DXYJyVQXTFDRE4uh/OVu3K8BnFA8WlOg+ZzABV6HfapFp8FEWkBDARmF9scVJ8FkzBKVgNwYhUfi9uFVWQ+UnUOck5X4DWrkoq6BjuwqicuBC4A1gIzRKT3UZyzohz1NRCRs4B+/HsAVzB9DqDirkO1+CyIyHwRKcAqjfyCVV23T1B9Fsz05od24KhGKWFbeZyzpO1VRbleA1Vdi/XFsM8CEWkC3AvMOdLzVrAjugYi0hN4H7hDVReX4pwlba9KyvU6VKPPwiVALFZbxjPAA8C4w5yzpO22MwmjZHsBH//95VCL//4aKIudBzlnEZB6FOetCBV1DUqyCLi0nM9ZHo74GohIL+Ab4DFVffWAp4PpcwAVdx1KEnKfBVX9O/DnarEWz3hDRJ5R1SKC7LNgqqRKoKoeYCnQ/4Cn+gPzj+LUC4AD19XsD6SoaklLktmmAq9BSTpiVU9UKUd6DQJVKt8Cj6vqCyXsEjSfA6jQ61CSjoTQZ6EEDqwf6vtWXQqqz4Ltre5V9YZVjPQA12N1d3sRqwtd48Dz44AZBxzTDusD/wGQEvi7Y7Hn93WheyFwzusDr1E1u9BVzDW4EzgPaAkcEziHAhfY/X7L4xpg9fzJxap6KN5VtGawfg4q8DqE+mfhKuBioA1W1+JBwDbgg2D9LNgeQFW+AbcAfwGFWL8uehd7bgrw1wH7/xX4wP/rdsA+fYBlgXNuAm6y+31W5jUA7gc2APlAGjAXOMPu91le1yDw+D/vv4TrFFSfg4q4DtXgs3BZ4N84GyuxrMJq8I4M1s+Cma3WMAzDKBXThmEYhmGUikkYhmEYRqmYhGEYhmGUikkYhmEYRqmYhGEYhmGUikkYhmEYRqmYhGEYR6nYQkBVbrI4wyhPJmEYRhmJyCwRmRgs5zWM8mIShmEYhlEqJmEYRhmIyBSsqRxu3bcGNdAk8HQHEVkUWMM6RUQ6HXDsiSIyO/D8NhF5VUTiDnZeEWkiIk4ReTOwjnS+iKwXkftFxPzfNSqd+dAZRtkMw5ph9G2gbuC2b/rqccCDQCesqamniogAiEh7rJXVvsBaF+ECrIkZ3zrMeR1YE9YNwpqc7mGs+Yiuqbi3aBglM3NJGUYZicgsYKWq3hZ43Bf4GRioqt8HtvXEWl2toapuFZF3Aa+qXlfsPB2BX4Haqrr7wPMe4vWfBLqo6oHTYhtGhTILKBlG+fm92N/bA/e1gK1Y60G3EJFLiu2zb2W15sDug51URG7Cmva6MRAJuIHN5RSzYZSaSRiGUX6KL3izr+juKHb/BvB8CcdtO9gJAwnmBaxlS+cDWcCtwPlHGathlJlJGIZRdh7+WTGttJYBx6jqhjKetxewSFX3d7cVkeZlfG3DKBem0dswyu4voFugF1MNSvf/6KnAMa+JyPEi0kJEzhKRSQc7b6An1Dqgk4icLiItReRRrN5UhlHpTMIwjLJ7Fqs0sBrYAzQ63AGq+jvQG6sL7mzgN6xeVbsOc95JwEfA+8CSwPHPlcu7MIwyMr2kDMMwjFIxJQzDMAyjVEzCMAzDMErFJAzDMAyjVEzCMAzDMErFJAzDMAyjVEzCMAzDMErFJAzDMAyjVEzCMAzDMErl/wHbSqKlEt4+MQAAAABJRU5ErkJggg==\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.fixed['theta'] | m.fixed['xi']):\n", " plt.figure()\n", " m.draw_mncontour('theta', 'xi', cl=[0.683, 0.95], size=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": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEFCAYAAADOj31RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzYklEQVR4nO3deXgNZ//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/5vbt20bHyldMV9BKqRLANkAB7YDKwJvAxft82yNzc3Pj/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/AslAHa11nKGpcoGPjw8RERHs27ePAQMGyF3whDCYq6sr48eP5/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/0q8BRmmtJ2Uxi01uqxysV1YCMPG2yoIdqTuc9mm/N3RbWYzRn1La0kTqj1QpQB9ST92ZTOrpQF5p748DNtzzPVVI/cv+M7Ar7euADO//czrQpLQx+6QtI7dP3bL0er0NPA9UAqqmjaGBjmZcJ1LPZLhJ6o/KGU81c7HlbZXD9bK1bdUL6AI8Reqpg12BM8DPZtpWFvmzMTqArU3AG8BJIJnU//kbZXjvB+DkPfOfTPvLnmm6Z57GwO60MWOB/ra+XsBg4BiQCFwBwoC2Zl2ntN//a32yWG+b2lY5WS8b3Fbd07bB36QW+UFSPyB0Ntu2etxJ7mYnhBAmJceghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghRDCpKSghXgApdR/lFLHlFKJSqnflFIuRmcS+YMUtBD3oZQaQ+qz+foBwaTeIH68oaFEviH3gxYiG0qpICASqK+1jkh7bRDwX621q6HhRL4ge9BCZO99YOs/5ZzmElDGoDwin5GCFiILSqmCwHPAr/e85Qz8lfuJRH4khziEyIJSqhaphzeSgDsZ3ioI7NFa1zEkmMhXChgdQAiTepLUp0D7k/qQ1X8sALYZkkjkO1LQQmStOBCvtT76zwtKqTJAADDIqFAif5Fj0EJkLR4oqpTK+G9kCBBxz4eGQliN7EELkbWNpP77GKaUmgN0BnoB9Q1NJfIV2YMWIgta60vAS0Af4BDQAmic8ZCHENYmZ3EIIYRJyR60EEKYlBS0EEKYlBS0EEKYlBS0EEKYlBS0EEKYlBS0EEKYlBS0EEKYlBS0EEKY1P8DJLO/mVGXUlQAAAAASUVORK5CYII=\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.fixed['theta'] | m.fixed['xi']):\n", " fig, ax = plt.subplots(1,1)\n", " con = m.mncontour('theta', 'xi', cl=0.95, size=200)\n", " con = np.vstack([con, con[0]]) # close 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.12" } }, "nbformat": 4, "nbformat_minor": 4 }