{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Binomial, Poisson and Gaussian distributions\n", "\n", "Python notebook for illustrating the Binomial, Poisson and Gaussian distributions and how they in certain limits converge towards each other (and in the end into the Gaussian). The notebook also illustrates simple fitting.\n", "\n", "## References:\n", "- Barlow: Chapter 3\n", "\n", "## Author(s), contact(s), and dates:\n", "- Author: Troels C. Petersen (NBI) and Christian Michelsen (NBI, first Python version)\n", "- Email: petersen@nbi.dk\n", "- Date: 16th of November 2020\n", "\n", "***" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np # Matlab like syntax for linear algebra and functions\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt # Plots and figures like you know them from Matlab\n", "from iminuit import Minuit # The actual fitting tool, better than scipy's\n", "import sys # Modules to see files and folders in directories\n", "\n", "from scipy import stats\n", "from scipy.stats import binom, poisson, norm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We set up the parameters for the program:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# General settings:\n", "r = np.random # Random generator\n", "r.seed(42) # Fixed order of random numbers\n", "\n", "save_plots = False\n", "verbose = True\n", "N_verbose = 10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting PDFs:\n", "\n", "First, we want to have a look at the three PDFs (Binomial, Poisson, and Gaussian), and compare them given \"the same\" parameters. Of course we can't give the same input parameters, but at least we can require that they have the same mean and widths (almost - one can not force the width of the Poisson to match that of the Binomial, once the mean is matched. The Gaussian will then have to choose between these two slightly different widths!).\n", "\n", "### Problem parameters:\n", "* The Binomial PDF needs two parameters: Number of trials (n) and probability of succes (p).\n", "* The Poisson PDF needs one parameter: The expected number (lambda - but in Python we write it \"Lambda\" or \"lamp\", as \"lambda\" is reserved!).\n", "* The Gaussian PDF needs two parameters: Mean (mu) and width (sigma)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def func_binomial_pmf(x, n, p):\n", " return binom.pmf(x, n, p)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def func_poisson_pmf(x, lamb):\n", " return poisson.pmf(x, lamb)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "def func_gaussian_pdf(x, mu, sigma) :\n", " return norm.pdf(x, mu, sigma)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Range of outcome:\n", "xmin = -0.5\n", "xmax = 13.5\n", "\n", "# Function parameters:\n", "# Binomial:\n", "n = 20\n", "p = 0.2\n", "# Poisson:\n", "Lambda = n * p # We choose this, as this is the expected number of successes.\n", "# Gaussian:\n", "mu = Lambda # Same here - the central value is n*p.\n", "sigma = np.sqrt(n*p*(1-p)) # For a Binomial process, the variance is n*p*(1-p)\n", "sigma = np.sqrt(n*p) # Alternatively, one can use the Poisson width " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "xaxis = np.linspace(xmin, xmax, 1000)\n", "yaxis_binom = func_binomial_pmf(np.floor(xaxis+0.5), n, p) # np.floor: See below\n", "yaxis_poiss = func_poisson_pmf(np.floor(xaxis+0.5), Lambda) # np.floor: See below\n", "yaxis_gauss = func_gaussian_pdf(xaxis, n*p, np.sqrt(n*p*(1-p)))\n", "\n", "# Note that the np.floor function takes the integer/rounded (towards zero) value of numbers." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAHwCAYAAABzIyTUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxN19rA8d/KnEhEgmhIiogpAzKUoAkxtFqlaFFtDa0OKN667y1uB9WqS9VtXbfavjqZa6iaam4rpUpjCiVUDFGJNAiZkHm9f5yTc09kECROxPP9fPJJsvfaaz17nYTzZA1baa0RQgghhBBCCGFgZekAhBBCCCGEEKIqkSRJCCGEEEIIIcxIkiSEEEIIIYQQZiRJEkIIIYQQQggzkiQJIYQQQgghhBlJkoQQQgghhBDCjCRJQghxhymlJiulFt3itcOUUr+UcX6jUmpoSWWVUplKKZ9babeEdporpQ4opTKUUmMros4btBeulPqjstu5nXZv9NpUJKVUlFLqhTvRVkVQSj2jlNpi6TiEEKK8JEkSQohyUErFK6WuGRONZKXU10opZ0vHdT2t9SNa6/mlnHPWWp8CUErNU0q9dxtNjQeitNYuWuvZt1EPxngmK6Vyjf2bqZQ6qpR6wiz2HVrr5rfbzs2yVLu3SynVVCm1VCl1QSmVrpSKU0r9RynlZYl4tNaLtdYPWaJtIYS4FZIkCSFE+fXSWjsDwcADwJvXF1AG98K/rQ2BI7dyoVLKppRTy4yJnDPwKrBIKVXvVgO8VymlfIHfgHNAkNa6JtAROAk8aMnYhBDibnEv/EcuhBAVSmudCGwEAsA09WmqUmoncBXwUUrVV0qtVUpdUkqdUEq9eF01DkqpZcbpavuVUq0LTyilJiqlThrPxSql+l53rTKOCqQppY4ppbqanSh1GpZSSiulfJVSLwHPAOONozbrlFKvKaVWXlf+P0qpWSXU8xMQCXxsvL6ZUspVKbXAOHJxRin1ZmGyaJyGtlMp9ZFS6hIwuRx9vBnIAJoY6+islEowiyFeKfV3pdQhYz8sU0o5mJ1/0djvl4yvQ/3r+mGUcXQlQyk1RSnVRCm1yzjqslwpZVdKuzd6bUqllFqhlPrLGO92pZS/2bl5Sqk5Sqn1xrp/U0o1MTvf3fhapymlPgZUGU1NBnZqrf+mtU4w9ud5rfUsrfVSY31uSqnvja/XZePXplEmY/92M/veNEVUKeWglFqklEpRSqUqpfYUJrPG1/qU8R5OK6WeMTtuPvXz30qps8b+3qeUCr+ureXGn6cMpdQRpVRoeftZCCEqgiRJQghxk5RS3sCjwAGzw4OBlwAX4AzwDZAA1AeeBP5pnswAjwMrAHdgCbBaKWVrPHcSCAdcgXcwjKh4ml3bDjgF1AHeBr5TSrmXN36t9VxgMTDDOHLTC1gE9FBK1TLeow0wEFhYwvVdgB3AaOP1x4H/GOP1AToBQ4DnSojZA5haVnzKoCdgB8SWUXQA0ANoDLQChhmv7wJMM573xPB6LL3u2h5ACBCGYergXAyJozeG5HdQKW3e6LUpy0agKYY+2I/hNTA3yFinG3ACYz8ppeoAKzGMXNYxxtCxjHa6GcuXxQr4GsOI4P3ANeDjct7HUAz37w3UBkYA15RSNYDZwCNaaxegAxBTSh17gDb89+d/hXmSC/TG8JrVAtbeRGxCCFEhJEkSQojyW62USgV+AX4G/ml2bp7W+ojWOg+4D8O0pgla6yytdQzwBYZEqtA+rfW3Wutc4EPAAcMbdrTWK7TW57TWBVrrZUAc0Nbs2vPALK11rvH8H0DP27kxrXUSsB3obzzUA7iotd53o2uVUtYYEqp/aK0ztNbxwL8oer/ntNb/0Vrnaa2vlVLVAGP/XsHwxvifWuvUMpqebeynS8A6DG+6wZDsfKW13q+1zgb+AbRXSjUyu/Z9rXW61voIcBjYorU+pbVOw5DMBJXUYDlem1Jprb8y9k82htGe1kopV7Mi32mto40/Q4vN7udRINbs52UW8FcZTdUxP6+UGm0c8clUSn1ujCVFa71Sa31Va52BISHrVJ77AHIxJEe+Wut8rfU+rXW68VwBEKCUctRaJxn7t6S+WGSMIU9r/S/AHjBf+/WL1nqD1jofQ6LeuqR6hBCiskiSJIQQ5ddHa11La91Qaz3qujf7Z82+rg9cMr75LHQGaFBSea11Af8ddUIpNUQpFWN8Y5uKYWSjjtm1iVprfV3d9bl984FnjV8/SwmjSKWog2HU58x1MZV4v2VYbuxfJwzT7IYopV4uo7x5onAVKNxIo755LFrrTCDluniSzb6+VsL3JW7KUY7XpkRKKWul1HTjVL10IN54yvzasu7H/OdFU3Z/pmAYQSss/7HWuhaG5MrWGI+TUur/lGFqZDqGBLmWMeG9kYXAZmCpUuqcUmqGUspWa30FQ7I8AkgyTh1sUVIFSqn/VYbNOdKM/eh6g75wUKWvZRNCiAonSZIQQlQM86TlHOCulHIxO3Y/kGj2vXfhF8qwdscLOKeUagh8DowGahvf3B6m6BqUBkop8+/vN7Z5q/EWWg20UkoFAI9RfDpYaS5iGF1oeF1M5vdbUnulB2cYjdoI9LqZ64zOmcdinAZW+7p4blo5X5vSPI1himU3DAlBo8Jqy3FtEkV/XpT59yX4Eeh3gzr/F8PITTvjxg4R18VzBXAyK39f4RfGEcx3tNZ+GKbUPYZheiVa681a6+4YkrRjGPqrCOP6owkYpkO6GfsxjfL1hRBC3BGSJAkhRAXTWp8FfgWmGRe5twKGUzTpCFFK9TP+dfxVIBvYDdTAkFBcAFBKPYdxgwgzHsBYpZStUqo/0BLYcJNhJmNYP2QedxbwLYY1ItFa6z/LU5FxStRyYKpSysWYTPwNwzqnW2LcRKAHt7aD3hLgOaVUG6WUPYZpkb8ZE6/bUZ7XpjQuGF7jFAzJxz/LLl7EesDf7OdlLGZJSwkmA+FKqQ+VUg2MsdbB8HNiHs81INW4nu3t6+qIAZ4y/oyFYlhXh7GuSKVUoHHUKR1DgpyvlKqnlOptTEqzgUwgv4T4XIA8DP1oo5SaBNQsV08IIcQdIkmSEEJUjkEYRgvOAauAt7XWW83Or8EwNekyhrU7/Yx/oY/FsJ5nF4ZEJhDYeV3dv2HYAOAihrUkT2qtU24yvi8BP+O0sdVmx+cb2yzvVLtCYzCMPpzCsGZrCfDVTdYx0LhuJhPDwv6dGDYyuCla6x+BtzBsXpCEYereUzdbTwn1lue1Kc0CDFMAEzFsRrH7Jtq9iGGt2HQMSVbTsto1bqQRhmF08qBSKsNY/hyGfgHD1DtHDD9Du4FN11XzFoZ+u4zhNVhidu4+DMl0OnAUw/q8RRjeU/yvsZ1LGNY4jSohxM0YRgmPY+iTLMo3HVMIIe4YVXRauxBCiHuZUup+DNOk7jNbjC+EEELcU2QkSQghBGBaG/U3YKkkSEIIIe5lVX6nGOPc5k+AHCBKa13ehcRCCCHKyfhvbTKG6U89LByOEEIIYVEWGUlSSn2llDqvlDp83fEeSqk/lOEp6RONh/sB32qtX8TwcDkhhBAVTGt9xfhgWH/jxhNCCCHEPctS0+3mcd1fKo275MwBHgH8gEFKKT8MC08L/8MuaZccIYQQQgghhKgwFkmStNbbMex8Y64tcML4xPMcYCmGZ0okYEiUQNZQCSGEEEIIISpZVVqT1ICiW4AmAO2A2cDHSqmewLrSLlZKvQS8BFCjRo2QFi1KfMi3EEIIIYQQQrBv376LWuu6JZ2rSklSSU/a1lrrK8BzN7pYaz0XmAsQGhqq9+7dW8HhCSGEEEIIIaoLpdSZ0s5VpelrCYC32fdeGB5IJ4QQQgghhBB3TFVKkvYATZVSjZVSdhiejr7WwjEJIYQQQggh7jGW2gL8G2AX0FwplaCUGq61zgNGA5uBo8ByrfURS8QnhBBCCCGEuHdZZE2S1npQKcc3ABtutV6lVC+gl6+v761WIYQQQghRRG5uLgkJCWRlZVk6FCHELXBwcMDLywtbW9tyX6O01pUYkmXIxg1CCCGEqCinT5/GxcWF2rVro1RJ+0wJIaoqrTUpKSlkZGTQuHHjIueUUvu01qElXVeV1iQJIYQQQlQ5WVlZkiAJcZdSSlG7du2bHgmWJEkIIYQQ4gYkQRLi7nUrv7+SJAkhhBBCVHHW1ta0adOG1q1bExwczK+//grAuXPnePLJJyu17b179zJ27Ngyy0RFRfHYY4/ddN3m8cfExLBhw3+Xpk+ePJmZM2fedJ3lFRMTQ/v27fH396dVq1YsW7bMdO706dO0a9eOpk2bMnDgQHJyciq8/Xnz5jF69OgKrbM8cW/dupWQkBACAwMJCQnhp59+Mp3bt28fgYGB+Pr6MnbsWAqX5Vy6dInu3bvTtGlTunfvzuXLlys07qqoWiVJSqleSqm5aWlplg5FCCGEEKLCODo6EhMTw8GDB5k2bRr/+Mc/AKhfvz7ffvttpbYdGhrK7NmzK6Vu8/ivT5Iqm5OTEwsWLODIkSNs2rSJV199ldTUVAAmTJjAuHHjiIuLw83NjS+//PKOxXU7yhN3nTp1WLduHb///jvz589n8ODBpnMjR45k7ty5xMXFERcXx6ZNmwCYPn06Xbt2JS4ujq5duzJ9+vQ7dk+WUq2SJK31Oq31S66urpYORQghhBCiUqSnp+Pm5gZAfHw8AQEBgGFkol+/fvTo0YOmTZsyfvx40zXffPMNgYGBBAQEMGHCBNNxZ2dnJkyYQEhICN26dSM6OprOnTvj4+PD2rWGx1WajxJFR0fToUMHgoKC6NChA3/88UeZsT766KMcOnQIgKCgIN59910A3nrrLb744gtT/Dk5OUyaNIlly5bRpk0b06hObGysKZ7SEjVnZ2feeOMNWrduTVhYGMnJyeXqx2bNmtG0aVPAkKx5eHhw4cIFtNb89NNPphGuoUOHsnr16mLXl9YXZb0OX3/9Nc2aNaNTp07s3LmzxLgmT57M4MGD6dKlC02bNuXzzz8v1/2UN+6goCDq168PgL+/P1lZWWRnZ5OUlER6ejrt27dHKcWQIUNM169Zs4ahQ4eWWW91Y5EtwIUQQggh7kbvrDtC7Ln0Cq3Tr35N3u7lX2aZa9eu0aZNG7KyskhKSioyRcpcTEwMBw4cwN7enubNmzNmzBisra2ZMGEC+/btw83NjYceeojVq1fTp08frly5QufOnXn//ffp27cvb775Jlu3biU2NpahQ4fSu3fvIvW3aNGC7du3Y2Njww8//MDrr7/OypUrS407IiKCHTt20KhRI2xsbEyJwS+//MKzzz5rKmdnZ8e7777L3r17+fjjjwFDsnDs2DG2bdtGRkYGzZs3Z+TIkcW2cb5y5QphYWFMnTqV8ePH8/nnn/Pmm2+yePFiPvjgg2Ix+fr6Fht9i46OJicnhyZNmpCSkkKtWrWwsTG8Tfby8iIxMbFYPWX1RUmvg42NDW+//Tb79u3D1dWVyMhIgoKCSuy3Q4cOsXv3bq5cuUJQUBA9e/bExcWF8PDwEssvWbIEDw+PcsVtbuXKlQQFBWFvb09iYiJeXl6mc+bXJycn4+npCYCnpyfnz58vs97qQJIkIYQQQogqrnC6HcCuXbsYMmQIhw8fLlaua9euFM6o8fPz48yZM6SkpNC5c2fq1q0LwDPPPMP27dvp06cPdnZ29OjRA4DAwEDs7e2xtbUlMDCQ+Pj4YvWnpaUxdOhQ4uLiUEqRm5tbZtzh4eHMnj2bxo0b07NnT7Zu3crVq1eJj4+nefPmJbZhrmfPntjb22Nvb4+HhwfJyclF3siDIcEqHOkKCQlh69atpvt85plnyqwfICkpicGDBzN//nysrKwo6fE4JS38L6svSnodLl68WOR1GDhwIMePHy8xpscffxxHR0ccHR2JjIwkOjqaPn36mH4GSnLhwoVyxV3oyJEjTJgwgS1btgCU+77vFZIkCSGEEEKU041GfO6E9u3bc/HixRLfFNvb25u+tra2Ji8vr8Q3v4VsbW1Nb4StrKxM11tZWZGXl1es/FtvvUVkZCSrVq0iPj6ezp07lxnrAw88wN69e/Hx8aF79+5cvHiRzz//nJCQkPLcaon3U9Y9mJcpz0hSeno6PXv25L333iMsLAwwrNlJTU0lLy8PGxsbEhISTNPTytsXpcVd3qTj+nJKKTIyMsocSWrZsmW54gZISEigb9++LFiwgCZNmgCGkaOEhIQiZQqvr1evHklJSXh6epKUlISHh0e57uNuVq3WJAkhhBBCVHfHjh0jPz+f2rVrl6t8u3bt+Pnnn7l48SL5+fl88803dOrU6ZbaTktLo0GDBoBh7c2N2NnZ4e3tzfLlywkLCyM8PJyZM2eW+GbfxcWFjIyMW4qrJM888wwxMTHFPgoTpJycHPr27cuQIUPo37+/6TqlFJGRkaZy8+fP5/HHHy9W/832Rbt27YiKiiIlJYXc3FxWrFhRatk1a9aQlZVFSkoKUVFRPPDAA7i4uJR4PzExMfj5+ZU77tTUVHr27Mm0adPo2LGj6binpycuLi7s3r0brTULFiwwXd+7d2/mz59fZr3VjSRJQgghhBBVXOGapDZt2jBw4EDmz5+PtbV1ua719PRk2rRpREZGmrYQv9U3uePHj+cf//gHHTt2JD8/v1zXhIeHU69ePZycnAgPDychIaHEJCkyMpLY2NgiGzdUpuXLl7N9+3bmzZtn6tvC6Wzvv/8+H374Ib6+vqSkpDB8+PBi199sX3h6ejJ58mTat29Pt27dCA4OLrVs27Zt6dmzJ2FhYbz11luljghdr7S4165dy6RJkwD4+OOPOXHiBFOmTDHdd+Eao08//ZQXXngBX19fmjRpwiOPPALAxIkT2bp1K02bNmXr1q1MnDixXPHczVRZQ7B3G6VUL6CXr6/vi3FxcZYORwghhBDVwNGjR2nZsqWlwxD3iMmTJ+Ps7Mzf//53S4dSrZT0e6yU2qe1Di2pfLUaSZItwIUQQgghhBC3SzZuEEIIIYQQooqYPHmypUMQVLORJCGEEEIIIYS4XZIkCSGEEEIIIYQZSZKEEEIIIYQQwowkSUIIIYQQQghhplolSUqpXkqpuWlpaZYORQghhBCiwlhbW9OmTRsCAgLo378/V69eLbN8hw4d7lBk5VcYU3x8PEuWLDEdnzdvHqNHj77h9Z07d2bv3r23HUdUVBSPPfbYDcuVN67yyM/PJygoqNR2s7OzGThwIL6+vrRr1474+HjTuWnTpuHr60vz5s3ZvHlzhcQjbqxaJUmyBbgQQgghqiNHR0diYmI4fPgwdnZ2fPbZZ2WW//XXX+9QZOVXGNP1SdK94N///neZz9r68ssvcXNz48SJE4wbN44JEyYAEBsby9KlSzly5AibNm1i1KhR5X6Ir7g91SpJEkIIIYSo7sLDwzlx4gQAH374IQEBAQQEBDBr1ixTGWdnZwCSkpKIiIgwjULt2LGD/Px8hg0bRkBAAIGBgXz00UcAxMTEEBYWRqtWrejbty+XL18GDCM4EyZMoG3btjRr1owdO3YUi2nUqFGsXbsWgL59+/L8888Dhjf/b775ZpGYJk6cyI4dO2jTpo2p7XPnztGjRw+aNm3K+PHjb9gHI0eOJDQ0FH9/f95++23T8UaNGvH666/Tvn17QkND2b9/Pw8//DBNmjQpklimp6fTt29f/Pz8GDFiBAUFBQB8/fXXNGvWjE6dOrFz505T+XXr1tGuXTuCgoLo1q0bycnJN4yxUEJCAuvXr+eFF14otcyaNWsYOnQoAE8++SQ//vgjWmvWrFnDU089hb29PY0bN8bX15fo6Ohyty1unTwnSQghhBCivDZOhL9+r9g67wuER6aXq2heXh4bN26kR48e7Nu3j6+//prffvsNrTXt2rWjU6dOBAUFmcovWbKEhx9+mDfeeIP8/HyuXr1KTEwMiYmJHD58GIDU1FQAhgwZwn/+8x86derEpEmTeOedd0yJV15eHtHR0WzYsIF33nmHH374oUhcERER7Nixg969e5OYmEhSUhIAv/zyC0899VSRstOnT2fmzJl8//33gGFaW0xMDAcOHMDe3p7mzZszZswYvL29S+2HqVOn4u7uTn5+Pl27duXQoUO0atUKAG9vb3bt2sW4ceMYNmwYO3fuJCsrC39/f0aMGAFAdHQ0sbGxNGzYkB49evDdd9/RsWNH3n77bfbt24erqyuRkZGmvnzwwQfZvXs3Sim++OILZsyYwb/+9S+2bdvGuHHjisXn5ORkGjl79dVXmTFjBhkZGaXeT2Jioul+bWxscHV1JSUlhcTERMLCwkzlvLy8SExMLLUeUXEkSRJCCCGEqOKuXbtGmzZtAMNI0vDhw/n000/p27cvNWrUAKBfv37s2LGjSJL0wAMP8Pzzz5Obm0ufPn1o06YNPj4+nDp1ijFjxtCzZ08eeugh0tLSSE1NpVOnTgAMHTqU/v37m+rp168fACEhIUXWyxQKDw9n1qxZxMbG4ufnx+XLl0lKSmLXrl3Mnj37hvfXtWtXCpdL+Pn5cebMmTKTpOXLlzN37lzy8vJISkoiNjbWlCT17t0bgMDAQDIzM3FxccHFxQUHBwdTQti2bVt8fHwAGDRoEL/88gs2NjZ07tyZunXrAjBw4ECOHz8OGEaDBg4cSFJSEjk5OTRu3BiAyMhIYmJiSo3z+++/x8PDg5CQEKKiokotp7UudkwpVepxUfkkSRJCCCGEKK9yjvhUtMI1SeZKegN9vYiICLZv38769esZPHgwr732GkOGDOHgwYNs3ryZOXPmsHz5ctO0t9LY29sDhg0k8vLyip1v0KABly9fZtOmTURERHDp0iWWL1+Os7MzLi4uN4yzsP6y2ih0+vRpZs6cyZ49e3Bzc2PYsGFkZWUVq8vKyqpIvVZWVqZ6r080Cr8vLQEZM2YMf/vb3+jduzdRUVFMnjwZ4IYjSTt37mTt2rVs2LCBrKws0tPTefbZZ1m0aFGR8l5eXpw9exYvLy/y8vJIS0vD3d3ddLxQQkIC9evXL7VvRMWRNUlCCCGEEHehiIgIVq9ezdWrV7ly5QqrVq0iPDy8SJkzZ87g4eHBiy++yPDhw9m/fz8XL16koKCAJ554gilTprB//35cXV1xc3MzrTdauHChaVSpvNq3b8+sWbOIiIggPDycmTNnFosHwMXFpcypZzeSnp5OjRo1cHV1JTk5mY0bN950HdHR0Zw+fZqCggKWLVvGgw8+SLt27YiKiiIlJYXc3FxWrFhhKp+WlkaDBg0AmD9/vul44UjS9R+FU+2mTZtGQkIC8fHxLF26lC5duhRLkMAw+lVY77fffkuXLl1QStG7d2+WLl1KdnY2p0+fJi4ujrZt2970/YqbJyNJQgghhBB3oeDgYIYNG2Z60/zCCy8UmWoHhu2uP/jgA2xtbXF2dmbBggUkJiby3HPPmTYrmDZtGmB48z9ixAiuXr2Kj48PX3/99U3FEx4ezpYtW/D19aVhw4ZcunSpxCSpVatW2NjY0Lp1a4YNG4abm9tNtdO6dWuCgoLw9/fHx8eHjh073tT1YEjoJk6cyO+//05ERAR9+/bFysqKyZMn0759ezw9PQkODjbtJDd58mT69+9PgwYNCAsL4/Tp0zfd5vUmTZpEaGgovXv3Zvjw4QwePBhfX1/c3d1ZunQpAP7+/gwYMAA/Pz9sbGyYM2cO1tbWt922uDFVnqHau01oaKiuiH30hRBCCCGOHj1a5vbNQoiqr6TfY6XUPq11aEnlq9V0O3mYrBBCCCGEEOJ2VaskSR4mK4QQQgghhLhd1SpJEkIIIYQQQojbJUmSEEIIIYQQQpiRJEkIIYQQQgghzEiSJIQQQgghhBBmJEkSQgghhKjikpOTefrpp/Hx8SEkJIT27duzatWqSm937969jB07ttLbAcNzg3744QcAZs2axdWrV03nnJ2d70gM5ZGfn09QUBCPPfZYieezs7MZOHAgvr6+tGvXjvj4eNO5adOm4evrS/Pmzdm8eXOlxNeoUSMuXrxYrrLmfW4pW7duJSQkhMDAQEJCQvjpp59KLHfp0iW6d+9O06ZN6d69O5cvXwZAa83YsWPx9fWlVatW7N+/v0LikiRJCCGEEKIK01rTp08fIiIiOHXqFPv27WPp0qUkJCRUetuhoaHMnj270tsBePfdd+nWrRtQPEmqSv7973+X+dysL7/8Ejc3N06cOMG4ceOYMGECALGxsSxdupQjR46wadMmRo0aZXpYraWY97ml1KlTh3Xr1vH7778zf/58Bg8eXGK56dOn07VrV+Li4ujatSvTp08HYOPGjcTFxREXF8fcuXMZOXJkhcQlSZIQQgghRBX2008/YWdnx4gRI0zHGjZsyJgxYwCIj48nPDyc4OBggoOD+fXXXwGIiooqMtoxevRo5s2bB8DEiRPx8/OjVatW/P3vfwdgxYoVBAQE0Lp1ayIiIorVER0dTYcOHQgKCqJDhw788ccfAMybN49+/frRo0cPmjZtyvjx44vdQ3R0NP369QNgzZo1ODo6kpOTQ1ZWFj4+PgAMGzaMb7/9ltmzZ3Pu3DkiIyOJjIw01fHGG2/QunVrwsLCSE5OLtbG5MmTGTp0KA899BCNGjXiu+++Y/z48QQGBtKjRw9yc3OBoiMte/fupXPnzuV+LRISEli/fj0vvPBCqWXWrFnD0KFDAXjyySf58ccf0VqzZs0annrqKezt7WncuDG+vr5ER0cXu37kyJGEhobi7+/P22+/bTreqFEj3n77bYKDgwkMDOTYsWMApKSk8NBDDxEUFMTLL7+M1rpYnfn5+QwbNoyAgAACAwP56KOPgP/2OcCGDRto0aIFDz74IGPHjjW97uXt13fffZcHHniAgIAAXnrppRLjKElQUBD169cHwN/fn6ysLLKzs8vs16FDh7J69WrT8SFDhqCUIiwsjNTUVJKSksrVdllsbrsGIYQQQoh7xPvR73Ps0rEKrbOFewsmtJ1Q6vkjR44QHBxc6nkPDw+2bt2Kg4MDcXFxDBo0iL1795Za/tKlS6xatYpjx46hlCI1NRUwvMndvHkzDRo0MB0rEmeLFmzfvh0bGxt++ITNUU8AACAASURBVOEHXn/9dVauXAlATEwMBw4cwN7enubNmzNmzBi8vb1N1wYHB3PgwAEAduzYQUBAAHv27CEvL4927doVaWfs2LF8+OGHbNu2jTp16gBw5coVwsLCmDp1KuPHj+fzzz/nzTffLBbjyZMn2bZtG7GxsbRv356VK1cyY8YM+vbty/r16+nTp0+p/bJt2zbGjRtX7LiTk5Mp8Xz11VeZMWMGGRkZpdaTmJhouncbGxtcXV1JSUkhMTGRsLAwUzkvLy8SExOLXT916lTc3d3Jz8+na9euHDp0iFatWgGGUZf9+/fzySefMHPmTL744gveeecdHnzwQSZNmsT69euZO3dusTpjYmJITEzk8OHDAMVe36ysLF5++WW2b99O48aNGTRoUJHz5enX0aNHM2nSJAAGDx7M999/T69evfjggw9YvHhxsZgiIiKKjVKuXLmSoKAg7O3ti5VPTk7G09MTAE9PT86fP1+sv837tbDsrapWSZJSqhfQy9fX19KhCCGEEEJUildeeYVffvkFOzs79uzZQ25uLqNHjyYmJgZra2uOHz9e5vU1a9bEwcGBF154gZ49e5pGDDp27MiwYcMYMGCAadTHXFpaGkOHDiUuLg6llGkEAaBr1664uroC4Ofnx5kzZ4q8cbWxscHX15ejR48SHR3N3/72N7Zv305+fj7h4eE3vGc7OztTnCEhIWzdurXEco888gi2trYEBgaSn59Pjx49AAgMDCyyNqgkkZGRxMTElHr++++/x8PDg5CQEKKiokotV9IIilKq1OPXW758OXPnziUvL4+kpCRiY2NNSVLh6xISEsJ3330HwPbt201f9+zZEzc3t2J1+vj4cOrUKcaMGUPPnj156KGHipw/duwYPj4+NG7cGIBBgwYVSbbK06/btm1jxowZXL16lUuXLuHv70+vXr147bXXeO2110rtr0JHjhxhwoQJbNmy5YZlzZW3X29WtUqStNbrgHWhoaEvWjoWIYQQQlQ/ZY34VBZ/f3/TiA3AnDlzuHjxIqGhoQB89NFH1KtXj4MHD1JQUICDgwNgSEwKCgpM12VlZZmOR0dH8+OPP7J06VI+/vhjfvrpJz777DN+++031q9fT5s2bYolDG+99RaRkZGsWrWK+Pj4ItPUzP/yb21tTV5eXrH7CA8PZ+PGjdja2tKtWzeGDRtGfn4+M2fOvGEf2Nramt74lla/eRxWVlZFrrGysjJdY94vhX0CNx5J2rlzJ2vXrmXDhg1kZWWRnp7Os88+y6JFi4qU9/Ly4uzZs3h5eZGXl0daWhru7u6m44USEhJM08wKnT59mpkzZ7Jnzx7c3NwYNmxYkRgL7+/6PrhRUuDm5sbBgwfZvHkzc+bMYfny5Xz11Vem8zeaGnejfs3KymLUqFHs3bsXb29vJk+ebIq7PCNJCQkJ9O3blwULFtCkSZMSY6hXrx5JSUl4enqSlJSEh4cHQLn69VbImiQhhBBCiCqsS5cuZGVl8emnn5qOmW9qkJaWhqenJ1ZWVixcuNC0GUDDhg2JjY0lOzubtLQ0fvzxRwAyMzNJS0vj0UcfZdasWaZk6OTJk7Rr1453332XOnXqFHnjWdhOgwYNAExrm25GREQEs2bNon379tStW5eUlBSOHTuGv79/sbIuLi5lTmm7HY0aNWLfvn0ARZLPwpGk6z8Kp9pNmzaNhIQE4uPjWbp0KV26dCmWIAH07t2b+fPnA/Dtt9/SpUsXlFL07t2bpUuXkp2dzenTp4mLi6Nt27ZFrk1PT6dGjRq4urqSnJzMxo0bb3g/ERERpiRk48aNpl3fzF28eJGCggKeeOIJpkyZUmwHuBYtWnDq1CnTqNCyZctu2K65woSoTp06ZGZmmtY5Abz22msl9mthgpSamkrPnj2ZNm0aHTt2LLUN836dP38+jz/+uOn4ggUL0Fqze/duXF1db3uqHVSzkSQhhBBCiOpGKcXq1asZN24cM2bMoG7dutSoUYP3338fgFGjRvHEE0+wYsUKIiMjqVGjBgDe3t4MGDCAVq1a0bRpU4KCggDIyMjg8ccfJysrC621aRH/a6+9RlxcHFprunbtSuvWrfn5559NcYwfP56hQ4fy4Ycf0qVLl5u+j3bt2pGcnGzaFKJVq1Z4eHiUOAry0ksv8cgjj+Dp6cm2bdtuuq2yvP322wwfPpx//vOfxdZD3apJkyYRGhpK7969GT58OIMHD8bX1xd3d3eWLl0KGEYEBwwYgJ+fHzY2NsyZMwdra+si9bRu3ZqgoCD8/f3x8fEpM2kwv59BgwYRHBxMp06duP/++4uVSUxM5LnnnjONoE2bNq3IeUdHRz755BN69OhBnTp1iiVvN1KrVi1efPFFAgMDadSoEQ888EC5r/344485ceIEU6ZMYcqUKQBs2bIFDw8PXnjhBUaMGEFoaCgTJ05kwIABfPnll9x///2sWLECgEcffZQNGzbg6+uLk5MTX3/99U3FXhpV3p0n7iahoaG6rAWLQgghhBDldfTo0TK3fBaiOsjMzMTZ2RmtNa+88gpNmzYtcfrh3aqk32Ol1D6tdWhJ5WW6nRBCCCGEEPe4zz//nDZt2uDv709aWhovv/yypUOyKJluJ4QQQgghxD1u3Lhx1Wrk6HbJSJIQQgghhBBCmJEkSQghhBDiBqrjGm4h7hW38vsrSZIQQgghRBkcHBxISUmRREmIu5DWmpSUFNPzw8pL1iQJIYQQQpTBy8uLhIQELly4YOlQhBC3wMHBAS8vr5u6RpIkIYQQQogy2Nra0rhxY0uHIYS4g2S6nRBCCCGEEEKYkSRJCCGEEEIIIcxUqyRJKdVLKTU3LS3N0qEIIYQQQggh7lLVKknSWq/TWr/k6upq6VCEEEIIIYQQd6lqlSQJIYQQQgghxO2SJEkIIYQQQgghzEiSJIQQQgghhBBmJEkSQgghhBBCCDOSJAkhhBBCCCGEGUmShBBCCCGEEMKMJElCCCGEEEIIYUaSJCGEEEIIIYQwI0mSEEIIIYQQQpiRJEkIIYQQQgghzEiSJIQQQgghhBBmbCwdgBDi3qC1JjM3kwtXL3Dx2kVSs1NJz0knIyeD9Jx00rP/+/W1vGvk5OeQXZBNdl422fnZhu/zs9HoYnUrFAC21rbYW9vjYO2AvY3hs521HQ7WDtSwrYGrvSs17WqaPte0r4mrnSvuDu7UdapLTbuaKKXudNcIIYQQooqRJEkIUSFOXbzMwK++52rBBbC9hLJJAZt0sE43fs5AWeWWeK3WVlDgCPmOhs/aDgpsQduDdgZt898P0wD49cmSBpUPKtf4kQdWOaCuGL63ygarq2B9FaUKSo6jwAbya0Keq/GzCzqvFuTW5j6nBqx+uRcu9s4V1mdCCCGEqJokSRJClFt+QT6JmYmcSD3BqbRTnEo9xdmMsyRkJnDx2kWo998UxgZ7HK1q42hVCwerBjhY1cLRyg1HKzccrGphr1ywUzWwtXLGBvs7NoKjtSaPbHIKMsnRho/sgnSuFVw2flwiq+Ay1wrOc63gGHlkA3Ae6LD0n9R2qI23izdeLl40rNkQ31q+NKnVBG8Xb2ys5J9UIYQQojqQ/9GFuIMKCjSXruZYOoxyuZx1ieOpR4lLPcbp9FOcST/Nnxnx5Bb8N/66jvXwcvambb2OHDxtzeV0Z+Y+9RDeNb1xs3e766euaa1Jy05j6YEYZm77lec71yRLnychM4G9yXv5/tT3prJ2VnY0cm1Ek1pN8K3lS3O35vjV9qOuU10L3oEQQgghboUkSULcQRO/O8TyvQmWDqMYZZ2BlUMi1o6Jhs8OiVjZppnOF+TUoiCnHgXZYeRne1CQXY+CnLpkFDhwyqyep9vdT2uPwDt/A5VEKUUth1q0qhtAXvoVutwXRphPbdP5q7lXOZV2ihOpJziZepITqSeIOR/DxtMbTWU8HD3wq+OHf21//GobPtd2rF1Sc0IIIYSoIiRJEuIOOpNylcZ1avB8x0YWi0HrAlJyzpJw7SiJWUdJvBZLau5fxrMKd9v61HNow30OvtSzb0I9+ybYWzvdsF6lFA/516vc4C2klpMdAKnXjQI62ToRUCeAgDoBRY5n5mRy/PJxjqQcITYlliMpR/j57M+mTSfq16hPUL0ggj2CCfIIokmtJlgp2WxUCCGEqCqqVZKklOoF9PL19bV0KEKUKCxjK93zthEQV/OOtZmvNbFk8xtZHCCLGLJIx7BxgTvWhGBPEO74Y09L7KmRZwWZZw0fbLu5xo5XfPxVQZO8Av5pY83lK/8uV3lnO2eC6wUTXC/YdOxK7hWOphzlSMoRDl44yG9Jv7H+1HoAXOxcCPIIIsgjiNB6oQTUCZD1TUIIIYQFKa2Lb6d7twsNDdV79+61dBhCFBMz5UGa61M4Ngi4ceFbpNGcIZ/dKpvdKptolUOGMvyeN9bWBGk708f9WJu2zxalK8hIxio1nq86/MDzDz1QIXVqrUnISGD/+f0cOH+A/ef3czrtNADOts6E3hdKe8/2tK/fnkY1G93167uEEEKIqkYptU9rHVrSOflTpRB3kGPBFU7VaIX/8E0VWm9mTia7knaxPWE7u5N289eVC4BhWlf3+mGEeYbR9r62shbmFln9/i2sHE5uxsUKq1MphXdNb7xrevO47+MAXMq6xJ6/9rDr3C52J+0m6mwUAPfVuI8wzzA61u9IhwYdqGl350YihRBCiHuRJElC3EFO+iopNhXznJ2zGWfZnrCdqLNR7E3eS15BHjXtatLOsx0vBr5Ie8/2eLl4yQhERXB0AyAvM6VSm3F3cOfhRg/zcKOHATibfpZdSYaE6cc/f2T1idXYKBtC6oXQybsTnb06413Tu1JjEkIIIe5FkiQJcQfV0Fc4b3trSZLWmt8v/s4Pf/7A9rPbOZl2EoDGro0Z3HIwnbw70bpua1nLUhmc3AEouHLpjjZbONI0oPkA8gvyOXTxEFFno/j57M/M2DODGXtm4OPqQyfvTnTx7kKruq1kAwghhBCiAsi7KSHuEF1QgDNXybdzKf81WnPo4iG2xG9h65mtJF1JMowk3BfCE82eoJNXJ+6veX8lRi0AcDQkSVdSz7NwV7wFA3GjDn15ol5fUt3+4sSVaE5eiWb+4QV8ffhrXGzq0My5Ay1cwqnv0Ax1EwmTvY01vVrXx9HOuhLjF0IIIe4OkiQJcYdkZ1/DQeVTcIMkSWvNwQsH2XLGkBj9deUvbKxs6FC/A6+0eYXO3p1xtXe9Q1ELwDSSlH8lhbfWHLFwMOaaGD6ssrBxPkpuzUPszVnPvtS1FOS6kpceSG56KwqyvKEcG3Q42VvzWKv6lR61EEIIUdVJkiTEHXIt4zIOgLYvedH92fSzrDu1jnUn15GQmYCtlS0d6ndgTNAYOnt3lsX6lmTnjLayZVzHOrwc0c3S0ZTiMQAyczPZlbSdqISt7EneTW7tX6jn5MlD9z9K9/sfxdulYbErk9Oz6Dn7FzKz8u500EIIIUSVJEmSEHdIzpVUAJT9f0eS0rLT2By/mXUn1xFzIQaFoq1nW0a0HkGX+7vgchNT80QlUgrl5I5TXhpOzvaWjqZMdbCnkVtfBvn1JT0nnW1/bmPj6Y0s/uNrFh77ktZ1W9O7SW8ebvSwaUTSyri5R1ZuviVDF0IIIaoMSZKEuENyMi8DoB1qsiNhB6tOrCLqbBS5Bbk0cW3Cq8Gv0tOnJ/fVuM/CkYoSObrDtTu7ccPtqmlXk8d9H+dx38c5f/U860+tZ+3JtUzZPYXp0dPp7N2Z3k16E1w3DICsvAILRyyEEEJUDZIkCXGH/JV2lvW1arI0ZS4pP6bi7uDOwOYD6dWkFy3dW8pW3VWdkzukJcC5GEtHcks8gOfcgxjm1oajGWdYl/QLG879ytYzW/Gwd+O+ur5cvDYcwzonIYQQ4t6mtNaWjqHChYaG6r1791o6DCHIL8hn57mdrDi+gu1nf6YAjX+NFjwf+iKR3pHYWttaOkRRXt8Oh8PfWjqKCpULbHdyZIWLM786OoCyIsIrgv7N+vNggwextpKd7oQQQlRfSql9WuvQEs9JkiRExbuUdYmVx1ey/Phy/rryF+4O7kRYNeSlP9aT2mcbgQGtLB2iuFmZFyBhj6WjqBxXU0jY8D9Mb/AQh53+IiUrhftq3McTTZ+gf7P+1HasbekIhRBCiApXVpIk0+1E1XLtMlxJsXQUt+yP9HgWn/6e9ee2k1OQS1id1oxvPoTO9UL586d5eOflcbWGbN99V3KuCy0etXQUlePqJbzWjuah3GZ81H8hUWejWPHHCubEzGHuobk80vgRnm35LC1rt7R0pEIIIcQdIUmSqDoK8uHfrSErzdKR3JR8IMrJkUU1Xdjr6IBjQQF9Mq/wdHoGTU7/CXvWAYaVHtnaBvsatSwarxDF2DoCYJV7DVsrW7o37E73ht05nXaaJUeXsObkGtaeXEtIvRCebfkskd6RMhVPCCFEtSZJkqg68rINCVLAk9Csh6WjuaGr+Tl8l3KARRd+IzEnlfp2rvxvnQfoWzsIVxvHYuV3nrzAv6Kz+cyham8hLe5BNg4UoFB514ocbuzamDfC3mBM8BhWxa1iydEljIsaRwPnBgxqMYgnmz1JDdsaFgpaCCGEqDySJImqIz/H8LlBCLTqb9lYynA56zJLji3hmz++IS07jWCPYP7uN5jO3p2xsSr9V+po2in266M42Mlf4EUVoxTZ2GOdd7XE0zXtajLUfyjPtHyGqLNRLDq6iJl7Z/J/h/6Pp5o/xTMtn5F1S0IIIaoVSZJElVGQl4MVMGXTSZZv2mzpcIqzuQSuP4NLNFjlwhV/SB3M/lON2L87H/ixzMuz8w3PoHG0lSRJVD3ZVg5Y52eVWcbGyoZuDbvRrWE3Dl88zFeHv+KL379gQewC+vr2Zaj/ULxcvO5QxEIIIUTlkSRJVBk5OddwAGq7OtO/qbelwzFJyzvLsazVnM35FVA0tAunuWMvarp7wU2G6VO3BrbWVpUSpxC3I1fZY1Nw7cYFjQLqBPBh5w85nXaaeUfm8W3ct6w4voKHGz3M8wHP09y9eSVGK4QQQlQuSZJElZGfa5hu17JBbSJ7+Vk4Goi7HMdnBz9jy5ktONk4MdjvWQb7Dea+GvdZOjQhKlyOlSO2NxhJKklj18a80+EdRrUexcLYhSw/vpwNpzfQ9f6ujGw9UpIlIYQQdyVJkkSVkZ+bDYCysbNoHObJUQ3bGrwY+CJD/IZQy0F2pRPVV66VA7Z5N58kFapXox5/f+DvvNjqRRYfXczC2IX8+OePdLu/GyNaj5BkSQghxF1FkiRRZeTlGJIkrG0t0v7xy8f57OBnbD2zlRq2NXip1UsM8RuCq70810hUf3nWjtjl3HqSVMjV3pVRbUbxTMtnWHR0EYtiF/HDnz/QvWF3Xm71siRLQggh7gqSJIkqoyCvcCTpzm6RfSb9DB8f+JhN8ZskORL3rHwbB+x1xT2jzNXelVfavMKzLZ81JUtbz2yle8PujA4ajY+rT4W1JYQQQlQ0SZJElVG4JulOjSSdv3qezw5+xndx32FnbceLgS8y1H+oJEfinpRv7Yi9zq7wes2TpYWxC03T8Pr49mFk65Gyxk8IIUSVJEmSqDIK1yRZVfJIUlp2Gl8f/prFRxeTp/MY0HwAL7V6iTqOdSq1XSGqsgIbRxzIRmuNUqrC63e1d2V00Giebvk0nx/6nGV/LGP9qfU83eJphgcOlz9OCCGEqFIkSRJVRkGeYSRJ2VTOSFJWXhaLjy7my8NfkpmTyaM+j/JKm1fwdqk6240LYSnamCTl5mvsbCo+SSrk7uDOhLYTeNbvWeYcmGPaPnx4wHCebvk0jjaOlda2EEIIUV6SJIkqoyDXsGi8okeSCnQBG05vYNa+WSRfTSa8QTj/E/w/soBcCDPa1gknssnKy8fOpvKf5dXAuQH/DP8nQ/2HMvvAbGbtn8WSo0sYEzyG3k16Y6XkeWJCCCEsR5IkUWUUjiRZVeAW4AfOH+CDPR/w+8Xfaenekmnh03jgvgcqrH4hqgtt64ijymHX6Yu4OjncwZY9eN73PTrWOciyk5/w1s63+PLgIgb5jqZ5rdYV1sr97k7Udbmzm8IIIYS4e0mSJKoMU5Jke/tv0BIyEpi1fxab4zfj4ejBex3fo1eTXvLXaSFKYefoAsAr83/lGncySTI3GJuaBznlsYlpMWPITQ8k+/wj6Fz32665xX0ubHo1ogJiFEIIcS+o8kmSUsoHeANw1Vo/ael4ROXRxiTJ+jZGkjJzMvn8989ZFLsIK2XFiNYjeM7/OZxsnSoqTCGqpWZedeEwzBvQmFx7NwtGEkBO/pNsTV7FFquVOLoepUu9x+nh2R9H61v7PV6w6wy/J1+r4DiFEEJUZ5WaJCmlvgIeA85rrQPMjvcA/g1YA19oraeXVofW+hQwXCn1bWXGKixP5xdu3HDzSZLWmvWn1/Ovvf/i4rWL9PLpxdjgsbK9sBDlZO1YC4B2ayMtHIlBF+BVa2tmu7myVq9kb+Jyxl1KpXfmFW52PPhBYA/+wMMVH6gQQohqqbJHkuYBHwMLCg8opayBOUB3IAHYo5RaiyFhmnbd9c9rrc9XcoyiitCm6XY3t24g7nIcU3+byr7kffjX9md25GwC6wZWRohCVF8tesIjH0B+xT8r6VbVA6YCg7LOM+3Cr7xlbc133n684fEgze1rl7uehJ3fcH/muUqLUwghRPVTqUmS1nq7UqrRdYfbAieMI0QopZYCj2utp2EYdRL3qMIkyca2fCNJmTmZfHrwUxYfXYyznTOT2k+in28/rK2sKzNMIaonexdo95KloyhRALBQF7DmxBo+2vcRA86uYlCLQbzS5hVc7FxueH1yzAF8Mv+s/ECFEEJUG5ZYxd4AOGv2fYLxWImUUrWVUp8BQUqpf5RR7iWl1F6l1N4LFy5UXLTijimcbmd9gy3AtdasP7We3qt7szB2IX18+7Cuzzr6N+svCZIQ1ZSVsqJv076s62v4XV9ydAm9VvVi3cl1aK3LvFZb22NHzh2KVAghRHVgiSSppKcUlvo/nNY6RWs9QmvdxDjaVFq5uVrrUK11aN26dSskUHGHlWMk6VTqKZ7f/DwTd0zEw8mDxY8uZnKHybg5WHKhuRDiTnG1d+XNsDf55rFvaODcgNd/eZ1hm4ZxMvVkqddoa3vsySW/oOxkSgghhChkiSQpAfA2+94LkMniAvJzyNNW2NjaFjuVk5/DpzGf8uS6Jzl++Thvhb3F4kcXy9ojIe5R/rX9WfjoQt7p8A4n007y5LonmRMzh5z84iNG2sYeG1VAbq6MJgkhhCgfSyRJe4CmSqnGSik74ClgrQXiEFWMzs8hFxtsrIoONsacj2HAugF8cvATujXsxto+axnQfIBMrRPiHmelrOjXtB9r+6zl4UYP89nBz3hy3ZPsT95ftKBxCm9OtmwDLoQQonwqNUlSSn0D7AKaK6USlFLDtdZ5wGhgM3AUWK61PlKZcYi7g8rPJRcbbK0NP5aZOZm8t/s9hmwcwtW8q8zpOocZETOo7Vj+Xa2EENWfu4M708On82m3T8nOy2bopqG8u+tdMnIyDAWMSVJedpYFoxRCCHE3qezd7QaVcnwDsKGi21NK9QJ6+fr6VnTV4k4oyCEbG2ysFdv+3MZ7v73HhasXeKblM4wJGiMPhBVClOnBBg+y6vFVzImZw6Kji4g6G8Ub7d7ApTBJypEkSQghRPlYYrpdpdFar9Nav+Tq6mrpUMStyM/lgpUt7/42kbHbxuJq78qiRxcxoe0ESZCEEOXiZOvEaw+8xpJHl+Du4M6rUa8yJ/tnLlpbkSsjSUIIIcqpWiVJ4u62mwu86O3EznPbGRs0lmWPLaNV3VaWDksIcRfyr+PPN499w6vBr3Io7yx9GnjyY9K2G24XLoQQQkAlT7cTojxSs1KZ+ttUNtmfxTcLpvb+Br86zSwdlhDiLmdrZcvwwOF4nMngm8RP+ODEJ8TkxvFm2Ju4O7hbOjwhhBBVmIwkCYv66c+f6LOmDz/8+QPPZLkxJVHR3L2ppcMSQlQjXo71WZCUzLP1HiPqbBR91/TlhzM/WDosIYQQVVi1SpKUUr2UUnPT0tIsHYq4gbTsNF7f8Tr/s+1/qOtUl6U9lzLgWg3yscHaqqTnDQshxK2xsnPABuhTK5xljy2jnlM9xkWNY8L2CaRly/8XQgghiqtWSZJs3HB3+CXxF/qt6ceG0xsY0XoESx5dQnP35lgV5JInM0CFEBXMysYBgPzcbJq6NWVxz8WMajOKLfFb6LOmDz+f/dnCEQohhKhqqlWSJKq2rLwspu6eysgfRlLTviaLey7mlTavYGttC2BIkpQkSUKIimVtZ0iSdK7hYbK2VraMbD2SJT2X4ObgxuifRjNl1xSu5cnDZoUQQhjIO9K7TUYyHF0Ld9kOTX9kXWDCuS2czLnEELc2jK3bHvsTO+HETlMZl5xkzioPC0YphKiOCpOkgtzsIsdb1m7J0p5L+c+B/zDvyDz2JO/h/fD3aVm7pSXCFEIIUYVIknS3if4/2PEvS0dRbgXAwpou/Nu9Fq75BfzfxRQ6nF7L/7N33/FRVekfx7/nzkwqkEDoRZAiXQVBBRWxUNxd1HVdF7sgzYI024oNRARUVFSkKGDvFX7qglQpIqCiiKCIiiAdifQkM+f3RxJEWkIyM3fm5vN+vfKKTGbufXZnF+eb85znSB8c8tyyktabJtEuEYDH+QJ5K0k5h56TlOBL0IAWA9S6amvdPfduXfHhFerTrI+uaXyNHEOzBQCUVISkeJO9VwqkSn2/cbuSAm3as0V3LxqmBZsW65yqZ2jQKbepbGL6EZ8/aPK3+nDVXl0axRoBeJ8vIVnSoStJB2pVtZXevvBtDVowSI8ueVRz183Vg2c+qEqplaJVJgAghngqJBljOknqVLduXbdL6CPwEQAAIABJREFUiRwblBy/lJrhdiVHNX3NdN0//37tzdmre1vdq0vrXSpjjj61LtMpI78/J0oVAigpAnntdso5ckiSpPSkdI1sO1LvrXpPD33+kC754BLd3/p+tavZLgpVAgBiiad6CUrEdLtQUHJi923bm7NXgxcMVt+ZfVUltYpe7/S6/n3CvwsMSJKUE7Ty+xj/DSC8/Im5K0m2gJAkScYY/bPeP/Vmpzd1XOnj1H9Wf90//36GOgBACRO7n7ZxeDYoGZ/bVRzWT5k/6coPr9Sb37+p6xpfp5f+9pJqp9Uu9OtzQiH5OSMJQJj9uZJ06J6kI6lZpqZe+NsL6ta0m9754R1d8X9XaHXm6ghVCACINYSkeBPKa7eLMVNWT9F/pvxHm3Zv0tPnPa0BLQYowZdwTNfICVoFfPxPEkB4BRISFLSmwHa7Q17nBNSneR+NOX+Mtu3dps5TOmvyj5MjVCUAIJbwiTTe2KDkxM5K0t6cvbp//v3676f/VYNyDfRmpzfVpnqbIl0rJ0S7HYDwS/D7tE8JMsFjC0n5WldrrTc7vanGGY1119y7dO+8e2m/AwCPIyTFm1AoZtrtVmeu1hUfXqG3f3hb3Zp204QOE1Q5tXKRr5cdDMkfw/utAMSngM9Rlvwyx7iSdKCKKRU1vv149Tixh95b9Z6u+L8r9OP2H8NYJQAglvCJNN7Y2BjcMGX1FHWe0llbdm/RM+c/oz7N+8hfzDbA3HY7VpIAhJfPMdqngEwoq1jX8Tt+9W7WW2Pa5bbfXf5/l+v9Ve+HqUoAQCyJvc0tOLqQu4Mb9gX36aGFD+ntH95W84rNNaLNiLCdI5ITCsnH4AYAEZCtgGpsXyS926vY12ot6S1/Hd2RvVJ3z7tbSz4fpYFJdZTo1uGzvoB09h1SWnV37g8AHuSpkFQSzkmyNqigHG39o/BTmsJl0+4NGvT5HVr5+3J1PuFadWnYSwr6tTFMtezJDqpsyrENewCAwvhUzdXRfi39Mi8s16sgabyk0UnSOG3Syj2b9NguqWooLJcvvFBQ+mOdVLW51KJLlG8OAN7lqZBkrZ0saXKLFi26u11LpPy8eYeytuxRh6HTo3pfX8oqJVV7VcbkaO9vV2v8dw01/v3ZYb9P+0acbg8g/B72d9fyRlX0wMVNwnZNn6Tekpr+Okt3fXqX/pPq0/A2w9W6auuw3aNAu7ZID9eRgtnRuycAlACeCkklQVZWlqxxNPSfTaNyP2utFm57W7M2P6+MhOq6pNpAZTSOXEtHqzoZEbs2gJIrwecoOxiZZZ62Ndrq1X+8qr4z++qGT25Q72a9dX2T6wt1iHax+QK534PF228FAPgrQlK8sSGF5NMVpx0X8Vvtyt6le+bdo5mbp6l9zfZ64IwHlBJIifh9ASDcAn6jpWsz9fgn30fsHq2TBym07xk98cUTem/5ArUp21sJTmT/zvSF9qm3pD179yg5oncCgJKFkBRnjA0qFIXNwaszV6vfzH765Y9fdGuLW3VNo2ui81tRAIiA+pXK6JPvNuq79X9E+E5/V6Bsun6u9KF+yuynPWuvVigrcm3EjkLqnST9uOF3ha+REABASIozxoYUivDk9hlrZuiuuXcp0Zeoce3G6dQqp0b0fgAQaeOvOSWKd/u7lmz8l26dfat2p47RkDOGqF3NdhG50+ad+5TziCNbjDOgAACHcv/AHRwTY4MRC0nWWo1dOlZ9ZvbR8WWO1+v/eJ2ABMATjDFR/WpRuYXe6PSG6pWtpwGzB2jM0jGysmG/T8BxlC0/e5IAIMwISXEmUu12e3L26PY5t+upr57S32v/XRM7TlTl1Mphvw8AlBQVUypqQocJurDOhRq9dHTuylL27rDew+8zypZfhul2ABBWnmq3KwnnJBkbkg1ztt2wa4P6zOyj77Z+p36n9FOXxl3YfwQAYZDoS9SQM4bohLInaOSSkfp1x68adc4oVSlVJSzXD/gc7ZRfCrGSBADh5KmVJGvtZGttj7S0NLdLiZjcdjtf2K63dPNSdZ7SWb/88YuePPdJdW3SlYAEAGFkjNG1ja/VU+c+pbU71qrz/3XWl5u+DMu1fQ4rSQAQCZ4KSSWBo1DY2u3eX/W+unzcRcn+ZL10wUs6u8bZYbkuAOBQZ1U/Sy///WWVTiitrv/rqnd/eLfY1/Q7RtnWJxMiJAFAOBGS4kw4VpKCoaAeWfSI7p53t5pVbKZX//6q6pb1bosiAMSK2mm19fLfXlbLSi117/x7Nfzz4coJ5RT5esYYZZmAHFaSACCsCElxxtjirSTtzt6tvjP76vnlz6tz/c4a026M0pPSw1ghAOBo0hLTNPr80bqq4VV66buX1Gdmn2INdMiRX8YSkgAgnAhJccbYoGwRQ9Km3Zt03cfXac66ObrrtLs08PSBCjiBMFcIACiI3/HrjlPv0D2n36N56+bp2o+v1cZdG4t0rRwF5DACHADCipAUZxyFitRut3LbSl3xf1fsH9BweYPLI1AdAOBYXFb/Mj113lP6dcevuuLDK7Ri24pjvkaO8cthTxIAhBUhKc4YGzrmlaS56+bq2o+vlZXV8xc8rzbV20SoOgDAsTqz2pl6vuPzMjK65qNrNGftnGN6fY7xy6HdDgDCipAUZ451cMMbK9/QzdNvVo3SNfTy315Wg3INIlgdAKAo6perr1f+/opqlaml3jN669UVrxb6tTkKsJIEAGFGSIozjgq3khSyIT26+FE98NkDal21tSZ1nKTKqZWjUCEAoCgqplTUpI6T1KZ6Gw1dOFTDPx+uYChY4OuCxi8fIQkAwoqQFGecQrTb7c3ZqwGzBmjSt5PUuX5njTp3lFIDqVGqEABQVCmBFD3e9nFd3ehqvfTdS+o7q6/25Ow56muCJkC7HQCEmadCkjGmkzFmXGZmptulRIxRUNYcud0uc1+mekzroelrpuv2lrfrrtPukt/xR7FCAEBx+Bzf/r+/Z/86W92mdtP2vduP+PwcE5DPFv2sJQDAoTwVkqy1k621PdLS0twuJWKOtpK0fud6XfPRNVq2ZZkeOfsRXd3oahljolwhACAcLm9wuR5r+5hWbF2hqz+6Wut2rjvs84LGLx8rSQAQVp4KSSWBo5DsYQY3fP/797rqo6u0afcmjW03Vu1rtXehOgBAOJ1X8zyNbz9eW/du1dUfXq2V21Ye8pygE5CfkAQAYUVIijPGhhQ6aCVp0YZFuu6j6yQrTeo4SS0rt3SnOABA2DWv1FwvdHxBjnF03cfX6fP1n//l5yEnwEoSAIQZISnOOApJB+xJmvrzVPWc1lPlU8rrxb+9qPrl6rtYHQAgEuqWrauX/vaSKqdWVq9Peunjnz7e/7OQCcjPniQACCtCUpzxKbh/T9Ir372iW2ffqkYZjfRCxxdUtVRVl6sDAERK5dTKmtRxkpqWb6rb5tymF5e/KIl2OwCIBMaexRknr91u1BejNP6b8Wpbo61GtBmhZH+y26UBACIsLTFN49qP038//a9GLBqhzbs3q5bjV0A5krUSw3oAICxYSYozVkG9W+oXjf9mvP5V7196rO1jBCQAKEESfYl6uM3D+k/9/2jitxP1fuqPCkpSiJY7AAgXQlIcyQ5m6+4K6VqYtFHXN7le97W6jzOQAKAE8jk+DTxtoHqc2ENfJK3X7RUylJW10+2yAMAzCElxYk/OHt0y8xb9r1SKOuw+Xn1P6csZSABQghlj1LtZb7Xb20hTS6Wq9+wB2p292+2yAMATCElxYEfWDvWa1kvz1s3TPVu2qU3W8W6XBACIEacHG2rw5q36bOMS9ZzWU39k/eF2SQAQ9whJMW7b3m26/n/X6+vNX2vEmUN12Y6dksPbBgDIZX0B/XPnLj3S8k4t27pMXT/uqi17trhdFgDENT5tx7ANuzbo2o+u1U+ZP2nUuaPU8bjzJUn2gHOSAAAlm3UCkqR2lU7V0+c+rTU71ujaj67Vbzt/c7kyAIhfntr1b4zpJKlT3bp13S6l2H7O/Fk9pvXQjqwdGtNujE6pdIqUtSv3h4QkAEAe60/I/Yfn2qm1E9A4v9GNabt19RvtNT4zW7WD7tYXFqffIJ3Z1+0qAJQgngpJ1trJkia3aNGiu9u1FMeq31ep29RusrJ6rsNzapTRKPcHodx/0+UfJgsAwC9lWuhlXaArG1SUJJ0saWLOTvXcuVRdMhI0rtTJqu8v5W6RxbFiivTLfEISgKjyVEjyghXbVqjH1B7yO3492/5Z1U6v/ecPbd6vA1lJAgDkyUoop4dsF115YYf9j9WXNCnzZ10/9Xpdn/W9xp097s9fuMWbTculYJbbVQAoYViSiCHfbvlW1//veiX6EzWp46S/BiRJCoVyvzO4AQCQJ+Azyg6GDnm8VlotTeo4San+VHX7Xzd9vflrF6oLAyfAQbkAoo5P2zHiq01fqdvUbiqdUFqTOk7ScWWOO/RJNr/djpUkAEAuv88oJ2QP+7MapWtoYseJSk9KV/ep3fXFxi+iXF0Y+PxSMNvtKgCUMISkGLBowyL1mNZDGckZmtRxkqqVqnb4J+btSZJDSAIA5PI5joIhK2sPH5SqlqqqiR0mqmJKRfX6pJcWrl8Y5QqLyZcghQhJAKKLkOSyBb8t0I2f3KgqqVU0scNEVU6tfOQnsycJAHCQgGMk6YirSZJUKbWSJnacqGqlqumm6Tdp3rp50Sqv+JwAe5IARB0hyUVz1s7RzdNvVo0yNTShwwRVSKlw9BfkryQx3Q4AkMfvy/13Qk7wyCFJksonl9eEDhN0fNrx6j2jt2b/Ojsa5RWfzy8F2ZMEILr4tO2SmWtmqs/MPqqTXkcT2k9QRnJGwS+ytNsBAP7Kn7eSlB06dHjDwcomldWz7Z9V/bL11XdmX03/ZXqkyys+J0C7HYCoIyS5YOaameo/u78almuoZzs8q/Sk9MK9MES7HQDgr/y+3JAULGAlKV9aYprGtR+nxuUb69bZt2r6mhgPSr4EBjcAiDpCUpTN/nX2/oA0tt1YlUkoU/gXM7gBAHCQ/Ha7wqwk5SudUFpjzh+jRuUb6dZZt2rGmhmRKq/4mG4HwAWEpCias3aO+s3qp/pl62tMuzEqnVD62C7A4AYAwEH2D24o5EpSvlIJpTTm/DFqmNFQA2YP0KxfZ0WgujCg3Q6ACwhJUTJ33Vz1ndlXddPrHvsKUr68lSTDShIAII+viCFJyltRajdGDco2UL9Z/WJzmIMvwEoSgKgjJEXBvHXz1GdGH9VNr6vx7ccrLTGtaBfKP0zW8YexOgBAPAvkT7c7hna7A5VJKKOx7ceqftn66jern+asnRPO8orPlyCFmG4HILr4tB1h89fN1y0zblHt9NrFC0iSbCgoI9FuBwDYL8GfG5I6PvGp8haVisb5t3xVx+vGT25RcP21srsbhKfAYupv1qirbx8fWABEFX/nRNCC3xbolpm36Pi04zW+XfECkiSFgkH5JAY3AAD2O7Need1yXj3tyw4W+1pZoQc0PfMBZVZ9UW3SblPVhJPDUGHxOIsS5FhWkgBEFyEpQj5f/7l6z+it48ocp/Htxxd+zPdRhII58kkyDl2SAIBcZZIC6t/uhLBdr+++F9VtajfN3/6InjrvKbWq2ips1y6KV5Ylydltc/fl8ktCAFHCp+0I+GrTV7p5xs2qUbqGnm3/rMomlQ3LdUOMAAcARFhaYprGtxuvWmm11GdmH32x8QtX6wk5gdx/YHgDgCgiJIXZ8q3LdeMnN6pCcgWNazdO5ZLKhe3aNm/jKtPtAACRlJ6UrnHtxqlSSiXdOP1GLduyzLVaQiav6SWY5VoNAEoeQlIYrfp9lXpO66lSCaX0bPtnVSGlQlivHwrmjQA3dEkCACIrIzlDz7Z/VumJ6eo5radWblvpSh37J7oy4Q5AFHkqJBljOhljxmVmZkb93mv+WKMe03oo4AT0bPtnVaVUlbDfw+a32/lYSQIARF6l1Ep6tv2zSvYnq8e0HlqduTrqNewPSbTbAYgiT4Uka+1ka22PtLTiTZE7Vut3rle3qd2UHcrW+PbjdVyZ4yJyn1Awr93OeOptAwDEsOqlq2t8+/GSpO5Tu2vtjrVRvX/I5O1JChGSAEQPn7aLafPuzeo2tZt2Zu3UuHbjVCe9TsTulb+SZHy02wEAouf4tOM1vv147QvuU7ep3bRh14ao3fvPlST2JAGIHm9+2l6/VBpSKeK3+d0x6lGxrDb7HI3bvF0Nn2kb0fuV3j/dzptvGwAgdp1Q9gSNPX+srp96vbpP7a6JHSeqfHL5iN/X7p9ux54kANHjzU/bqRWk07pF9BZ/hLLUc+s8/Zr9h0ZntNbJ1cM7pOFwdu0L6ukFW3Rc6boRvxcAAAdrXL6xRp83Wr0+6aUe03poYoeJxT4ovUAO7XYAos+bIalMVand4Ihdfnf2bt00rad+CO7UE+c9qVOrt4nYvQ6UuX2Pnpk7Q8MZ3AAAcEnzSs31xDlP6KbpN+mm6TdpXLtxSgmkROx+DG4A4Ab2JB2j7GC2+s/ur6+3fK3hZw1XmygFJEkKhawkyTEmavcEAOBgraq20og2I/TNlm/Uf3Z/ZUcwwFgfh8kCiD5C0jEI2ZAGzhuoeevm6Z7T71H7Wu2jfH9CEgAgNpxf83zd1+o+zVs3TwPnDlQwf99suNFuB8AFhWq3M8a8LWmCpI+staHIlhSbrLUa/vlwffTTR+rTvI8uPeHSqNeQt5Akn0NIAgC475J6l2j7vu16bMljKpNYRgNPGygT7l/k+Wi3AxB9hV1JekbSFZJ+MMYMM8Y0iGBNMWns12P1yopXdE2ja3R9k+tdqSGYl5JYSAIAxIquTbqqS+Muen3l6xq9dHTYr29ZSQLggkKtJFlrP5H0iTEmTdLlkqYZY36VNF7SS9ZaT//N9fqK1/X0V0/rwjoXakCLAeH/LVkh5bfbsZIEAIgl/U7pp8ysTI1ZOkbpiem6suGV4bu4jxHgAKKv0NPtjDEZkq6SdLWkLyW9LOlMSddKahuJ4mLBxz9/rAcXPqizq5+t+1vfL8e4t41rf0hiKQkAEEOMMbrn9HuUuS9Twz4fpjIJZdSpTqfwXNxJyP3OYbIAoqhQn/iNMe9I+lRSiqRO1toLrbWvW2t7SyoVyQLdNH/dfP330/+qWcVmeuTsRxTIX/J3yZ/tdoQkAEBs8Tt+DW8zXKdVPk33zLtHc9bOCct1jZ92OwDRV9hlkWettY2stQ9Za9dLkjEmUZKstS0iVp2Lvt78tfrO6qs6aXX05HlPKsmf5HZJCuWNzKDdDgAQixJ9iXri3CfUoFwD9Z/VX0s2Lin+Rfefk0S7HYDoKWxIGnKYxxaEs5BYsnr7at04/UZlJGVoTLsxKpNQxu2SJB24J8nlQgAAOILUQKpGnz9aVVKrqPeM3lr1+6riXdCX127HShKAKDrqx21jTGVjzCmSko0xzYwxzfO+2iq39c5zNu3epF6f9JLf+DWu3TiVTy7vdkn7BS3tdgCA2FcuqZzGthurJF+Sen3SSxt3bSzytUze4IZQDnuSAERPQWsSHSQ9Iqm6pJGSHs376i/prsiWFn07s3bqhk9uUOa+TD1z/jOqUaaG2yX9RSjE4AYAQHyoWqqqRp8/Wjuzd+qG6TdoR9aOIl0nPyRZBjcAiKKjhiRr7fPW2nMkXWetPeeArwutte9EqcaoyA5mq++svlq9fbUea/uYGmY0dLukQ3CYLAAgnjQo10CPtX1MP23/SX1n9lVWUYLO/pUk2u0ARM9RR4AbY66y1r4kqZYxpv/BP7fWjoxYZVFkrdW98+/VwvULNeSMIWpdrbXbJR0Wh8kCAOJNq6qtNPiMwbpr7l26e+7dGtZm2DEdp5E/3c759m1p68pIlemepDTp3HukgPsDogD8qaBzklLzvnt2zLckPfHFE5qyeop6N+uti+pe5HY5R8Q5SQCAeNSpTidt3rNZjy15TJVSK2lAiwGFf7E/WQuCjXTajg3SzqLvbYpJOXulPdukRhdLNVq6XQ2AAxw1JFlrx+Z9HxSdcqLvtRWv6bllz+nfJ/xb3Zt2d7uco/pzuh0hCQAQX7o07qKNuzZq0reTVDGloq5udHWhXuf3+XR59t364oZ2KpeaEOEqo+zHmdKLFzO5D4hBBbXbjTraz621t4S3nOiavma6hi4cqrbV2+qu0+6K+alxHCYLAIhXxhjd3vJ2bdq9SQ8velgVUiqoY62OBb4u/xeDOfmHBXpJ3n4rBQlJQKwpqN0uDKfAxaavNn2lO+bcoSblm2h4m+HyOwX9V+E+VpIAAPHM5/j00FkPqee0nrrr07uUkZShlpWP3mbmz/t3Xv4vCj0l/7NHiINygVhTULvd89EqJJp+zvxZvWf0VqWUSnrqvKeUEoiPI5/yf4lGRgIAxKskf5JGnTtK13x0jfrM6KMXLnhBdcvWPeLz968kBb0YkvJWkghJQMwpqN3ucWttX2PMZEmH/O1krb0wYpUVQ3YwpPWZew77s9/3blWfOT1lrdHg0x/Xvn3JWr/v8M+NNVt37ZMkObTbAQDiWFpimsacP0ZXfnilbpp+k17++8tHPLzd7/PwSpIv72MY7XZAzCmox+zFvO+PRLqQcFqxYYdaPTTj0B+YbKXUHCcncbN2/9JD/1r6vaTvo15fcSUFfG6XAABAsVQpVUVPnvekunzcRbfMuEXPdXhOyf7kQ56X/4vBHC+GJNrtgJhVULvdkrzvs40xCZIaKHdFaaW1NmaPvq6WnqyhlzT9y2PWhvTOumFasWOtLq0+UPUbt3KpuuJJSw6oToXUgp8IAECMa5zRWMPOGqa+M/tq4NyBeuTsRw45Q8nv5P45f1+up9BuB8SsQk0rMMb8XdIYST9KMpKON8b0tNZ+FMniiqpcaoI6n3rcXx57fMnjWrFjngacMkDXNfmPS5UBAIADnXvcubqt5W0asWiEHv/icfU/5a9n13t7T1JeZwghCYg5hR3p9qikc6y1qyTJGFNH0v9JismQdLB3f3hXzy17TpeecKmubXyt2+UAAIADXNXwKv3yxy+auGyiapSuoX+f8O/9P/P0dDtGgAMxyyn4KZKkTfkBKc9qSZsiUM9hGWMuNsaMN8a8b4xpfyyv/Wz9Zxq8YLBaVWkVF2chAQBQ0hhjdOepd+rMamfqwc8e1Px18/f/zOfz8DlJ+9vtCElArDlqSDLGXGKMuUTSt8aYD40x1xljrpU0WdKiwtzAGDPBGLPJGLPsoMc7GmNWGmNWGWPuPNo1rLXvWWu7S7pOUqF75VZvX63+M/urVlotPdr2UQXy/zICAAAxxe/49cjZj6hOeh0NmD1AP/z+Q97jHl5J2j+4IehuHQAOUdBKUqe8ryRJGyWdLamtpM2SyhbyHpMk/eVIbWOMT9LTki6Q1EjS5caYRsaYpsaYKQd9VTzgpXfnva5A2/Zu043Tb1TAF9BT5z2l0gmlC1kuAABwQ2ogVU+f97SS/cm6afpN2rJny597krwYkhgBDsSsgqbbdSnuDay1c4wxtQ56+FRJq6y1qyXJGPOapIustQ9J+sfB1zC5PXLDJH1krf3icPcxxvSQ1EOSjjvuOPWZ0Udb9mzRhA4TVK1UteL+xwAAAFFQObWynjrvKV338XW6efrNuqnBSEleXUliuh0Qqwq1J8kYk2SMuckYMzqvfW6CMWZCMe5bTdKvB/x5bd5jR9Jb0vmSLjXG9DrcE6y146y1Lay1LfYl7dNXm7/S0DOH6sQKJxajTAAAEG2NMhpp+FnDtXzrcj23coikkDdXkva327GSBMSawg5ueFFSZUkdJM2WVF3SjmLc93DTE474t5+1dpS19hRrbS9r7ZiCLp6Zlam+zfuqfa1jmvEAAABixDnHnaPbWt6mJVs+VUKFaQp5MSTtn27HShIQawobkupaa++RtMta+7ykv0tqWsBrjmatpBoH/Lm6pN+Kcb2/KJtYVl2bdA3X5QAAgAuuaniVzqnaSYnlZ2rh5k/cLif88g/Opd0OiDmFDUn568DbjTFNJKVJqlWM+y6SVM8Yc7wxJkFSZ0kfFON6f1GlVBVGfQMAEOeMMereaIBydh2v1396RN9s/sbtksLLmNx9SbTbATGnsCFpnDGmrKR7lBtmlksaXpgXGmNelbRAUn1jzFpjzPXW2hxJN0v6n6TvJL1hrf32mKs/9F6djDHj/sj8o7iXAgAAMSDJn6C9665S6UA53TLzFm3YtcHtksLLF2AlCYhBhQpJ1tpnrbW/W2tnW2trW2srWmvHFvK1l1trq1hrA9ba6tba5/Ie/9Bae4K1to619sHi/Ic44F6TrbU90tLSwnE5AADgMp9jZIOpuqrWIO3O3q0+M/toT84et8sKH8fPniQgBhV2ul2GMeZJY8wXxpglxpjHjTEZkS4OAACUbH4n96NKRkJNjWgzQt9t/U53z71b1npkkIPjp90OiEGFbbd7TdImSf+SdKmkLZJej1RRAAAAkuTz5e4xDoaszq5xtvqd0k9Tf5mqMUsLHHYbH2i3A2LSUQ+TPUA5a+0DB/x5iDHm4kgUBAAAkM+XN4gpmLdydF3j67Rq+yqNXjpatdNrq0OtDm6WV3y02wExqbArSTONMZ2NMU7e12WS/i+ShQEAAPic3JCUf5isMUb3tbpPJ1U4SXfPvVvLty53s7zio90OiElHDUnGmB3GmD8k9ZT0iqSsvK/XJPWLfHnHJn+6XWZmptulAACAMPDnhaRgMLT/sQRfgh4/53GlJ6Wr94ze2rx7s1vlFR/tdkBMOmpIstaWttaWyfvuWGv9eV+OtbZMtIosLKbbAQDgLfl7kvJXkvKVTy6vJ899UjuydqjfrH7KCma5UV7xOX4pyEoSEGsKuydJxpgLJbXJ++Msa+2UyJQEAACQK38l6Y892dq8Y99ffpYROF63Nb9Pgz+/U/fNHaL+zQe6UWKxlDN++UJBt8sAcJBChSRjzDDwDY0SAAAgAElEQVRJLSW9nPdQH2PMmdbaOyNWGQAAKPH8jiPHSKNmrNKoGasO+5yECm01Re/qrQVS9vbTolxh8UxO2KW6yXuV7HYhAP6isCtJf5N0srU2JEnGmOclfSmJkAQAACImwe9ownUt9evvRz5ANmQb6u11O/RLlcnqfGorVU9pFMUKi27Vxh3KXuJTdvY+QhIQYwrdbicpXdK2vH9m0w8AAIiKtvUrFvici/c9ocv/73JN3TJCr/39NVVKrRSFyopn7g9blL3EL9FuB8Scwo4Af0jSl8aYSXmrSEskDY1cWUXDdDsAAEqmtMQ0jTpnlHZl71L/Wf3jYpCD40hB6zACHIhBBYYkY4yRNFfS6ZLeyftqZa19LcK1HTOm2wEAUHLVLVtXQ88cqq+3fK0HFz4oa23BL3KRzxjlyMdKEhCDCgxJNvdvmPesteuttR9Ya9+31m6IQm0AAADH5Pya56t70+5654d39MbKN9wu56j8PqNs+WVYSQJiTmHb7T4zxrSMaCUAAABhcNPJN6lN9TYa9vkwLdm4xO1yjsgxRkE5MhwmC8Scwoakc5QblH40xnxtjPnGGPN1JAsDAAAoCp/j00NnPaRqpaup/6z+2rArNhtgfI5RtnwSIQmIOYUNSRdIqi3pXEmdJP0j7zsAAEDMKZNQRqPOGaV9wX3qN7Of9gX3FfyiKHOMUY78rCQBMeioIckYk2SM6SvpNkkdJa2z1v6S/xWVCgEAAIqgdnptDT1zqJZtXaZhnw9zu5xD+H1GObTbATGpoJWk5yW1kPSNcleTHo14RcXACHAAAHCgc487V92adtNb37+l91a953Y5f+EzRjnWL2MJSUCsKSgkNbLWXmWtHSvpUklnRaGmImMEOAAAONhNJ9+k0yqfpiGfDdGKbSvcLmc/x8kf3MB0OyDWFBSS9v+/1lp+zQEAAOKP3/FreJvhSktMU9+ZfZW5LzY6Tvx5gxsM5yQBMaegkHSSMeaPvK8dkk7M/2djzB/RKBAAAKC4MpIzNLLtSG3cvVF3zb1LIRtyu6Q/Bzfwe2gg5hw1JFlrfdbaMnlfpa21/gP+uUy0igQAACiukyqcpNtb3q45a+do/Nfj3S5HPid3cINjWUkCYo3f7QIAAACipXP9zlq6eame/uppNS3fVK2rtXatltyQ5JMTypK+ecu1OiKqekupbE23qwCOGSEJAACUGMYY3Xv6vVq5baXu+PQOvf6P11W1VFVXavE5Rptteu5K0tvXu1JDxNVrL135pttVAMeMkAQAAEqUlECKHj/ncXWe0lkDZg3Q8xc8rwRfQtTr8BmjScEOatL2Ul3azJ2gFlHv9JCydrtdBVAkhCQAAFDi1CxTU0POGKK+s/pq+OfDdU+re6Jeg+MYSUbbk2pI5WtH/f4Rl1haytnndhVAkRQ03S6ucJgsAAAorPNqnqcuTbroje/f0Pur3o/6/f2OkSQFQzbq944Kxy9xBhTilKdCEofJAgCAY3FLs1vUsnJLDflsiFb9viqq9/blhyTr5ZDEeHPEJ0+FJAAAgGPhd/wa0WaEUgOpGjB7gHZnR28PjWNyQ1LIqytJvoDEQbmIU4QkAABQopVPLq/hbYbrp8yfNOSzIbJRWtnZv5Lk/rm2keH4WElC3CIkAQCAEu+0KqfphpNu0OTVk/Xuqnejcs+8jKRgyKMpiXY7xDFCEgAAgKQeJ/bQ6VVO19CFQ7Vy28qI388YI59jvL0nKcjgBsQnQhIAAIAkn+PTQ2c9pNIJpXXr7Fu1K3tX5O9pjIfb7fzsSULcIiQBAADkKZ9cXiPajNCaHWs0eMHgiO9P8jlGIS+vJNFuhzhFSAIAADhAy8otdeNJN+rDnz7UWz+8FdF7+RyjnCAhCYg1hCQAAICDdD+xu1pVaaVhC4dpxbYVEbuPY8RKEhCDCEkAAAAHcYyjh856SOmJ6bp19q3ambUzIvfx+xwFvXpOEiEJccxTIckY08kYMy4zM9PtUgAAQJzLSM7Q8DbD9euOXzVowaCI7E9yjFGOZ0MS5yQhfnkqJFlrJ1tre6SlpbldCgAA8IAWlVuod7Pe+vjnj/Xm92+G/fo+Rwp5NST5AoQkxC1PhSQAAIBw69qkq1pXba0Ri0boh99/COu1fcbj5ySFciSv/ueDpxGSAAAAjsIxjh4880GlBlJ1+5zbtSdnT9iu7fMZ764kOf7c79arB0HBywhJAAAABSifXF4PnfmQVm1fpYcXPRy26/q8vidJouUOcYmQBAAAUAitq7VWl8Zd9Ob3b2raL9PCck3H8Xi7nSQFs92tAygCQhIAAEAh9W7WW00ymui++ffpt52/Fft6fqcEtNuxkoQ4REgCAAAopIAvoBFnj1DIhnTHnDuUU8wA4Bjj4XOSArnfQ0F36wCKgJAEAABwDGqUrqF7T79XX23+SqO/Gl2sa/kcL4ck9iQhfhGSAAAAjtHfav9NF9e9WM9+86wWrl9Y5Ov4SsKeJEIS4hAhCQAAoAj+e+p/VbNMTf330/9q295tRbqGt1eS8kMSgxsQfwhJAAAARZASSNEjZz+izH2ZumfePbJFWBHyeXpPUn5IYk8S4g8hCQAAoIjql6uvAS0GaM7aOXrpu5eO+fWOl1eSfLTbIX55KiQZYzoZY8ZlZma6XQoAACghLm9wuc6pcY4eW/KYVm5beUyv9TtGIfYkATHHUyHJWjvZWtsjLS3N7VIAAEAJYYzRoNaDlJ6Yrjvm3KG9OXsL/dqSsSeJkIT446mQBAAA4IaySWU15Iwh+jHzR41cMrLQr/P2OUl5ISlISEL8ISQBAACEQetqrXV1o6v16opXNWftnEK9xtsjwDknCfGLkAQAABAmfZr3Ub2y9XTPvHu0dc/WAp+f224XhcLc4ARyvxOSEIcISQAAAGGS6EvU8LOGa2fWTt07/94Cx4L7jFHI6+12hCTEIUISAABAGNUrW0/9W/TXnLVz9MbKN476XJ9jlBPy6FISIQlxjJAEAAAQZlc0uEJnVDtDDy9+WKu3rz7i83yOkVcXkjhMFvGMkAQAABBmxhgNOWOIUvwpuuPTO5QVzDrs87w9Ajx/cEO2u3UARUBIAgAAiIDyyeU1qPUgrdi2Qk99+dRhn1MiRoDTboc4REgCAACIkHOOO0f/PuHfmvTtJC1cv/CQn/sceTck+Zhuh/hFSAIAAIigW1vcqpplauquuXcpc1/mX37mcxwPn5PEniTEL0ISAABABKUEUjS8zXBt27tNgxcM/stYcJ8jD48A5zBZxC9CEgAAQIQ1ymikG0+6UVN/maqPfvpo/+M+Y5Tj2ZCUt5IUZHAD4g8hCQAAIAq6NOmikyqcpCELh2jDrg2SctvtvLuSxOAGxC9CEgAAQBT4Hb8ePPNB5YRydO+8e2WtzR3c4Nk9SfmDG9iThPhDSAIAAIiSmmVqasApA7Rg/QK9vvJ1OSXinCRWkhB//G4XAAAAUJJcVv8yzfx1ph5d/Kg6pA9XTshqzdbdbpcVdiYrWzUkQhLiEiEJAAAgiowxGtR6kC754BIt2PGUgqHr1ObhmW6XFXaJytLKJEkhBjcg/ngqJBljOknqVLduXbdLAQAAOKJKqZV09+l36/Y5t+vitt/rrAqd3S4p7J6ZsVLaJfYkIS55KiRZaydLmtyiRYvubtcCAABwNBccf4Fmrpmpab+8rB4t/6ZGGY3cLims3li0JjckrfxQ2rnR7XIio+zxUuub3a4CEeCpkAQAABBPBp4+UIs3LtZdn96l1zu9rkRfotslhY3P5+ibwElqun2NtH2N2+WEX/YeKXu31KKrlJDidjUIM0ISAACAS9IS0zT4jMG64ZMb9OQXT+rWlre6XVLY+Byje8s+pHdvPMPtUiJj/lPS1IEMpvAoRoADAAC46MxqZ+o/9f+jF5a/oEUbFrldTtg4xnj3oFyJw3I9jpAEAADgsv6n9Ff10tV199y7tTNrp9vlhIXPMd49KFc64BwoBlN4ESEJAADAZSmBFA09c6g27N6g4YuGu11OWDjGKBhyu4oIyg9JlpDkRYQkAACAGHByxZPVtUlXvbfqPc1ZO8ftcorNMaLdDnGLkAQAABAjbjjpBtVNr6tB8wcpc1+m2+UUi88xCnm63Y6Q5GWEJAAAgBiR4EvQkDOHaOverRqxaITb5RSL4/k9SfkhiXY7LyIkAQAAxJDGGY3VtUlXffDjB3Hddufz/HS7/MENrCR5ESEJAAAgxvQ6qZfqptfV/fPvj9u2O+9Pt6PdzssISQAAADEmv+1u295tcdt2l3tOkttVRBAhydMISQAAADGocUZjXd/0en3w4wea9esst8s5Zj5HCnq53c5wTpKXEZIAAABiVK8Te6le2XoavGBw3LXdOaaktNsRkryIkAQAABCjAr6AhpyR23Y3/PP4OmTWcYysp0MSgxu8jJAEAAAQwxplNFK3pt00efXkuGq78xnj7XY79iR5GiEJAAAgxvU8sadOKHuCBi2In0NmfQ4hCfGLkAQAABDj8tvutu/drmGfD3O7nEJxjJGXMxJ7kryNkAQAABAHGmY0VPcTu2vK6imauWam2+UUyPPT7diT5GmEJAAAgDjRvWl31S9bX4M/i/1pdw6HySKOEZIAAADiRMAX0JAzc9vuYv2QWZ8xCpWElSRLu50XEZIAAADiSINyDdS1aVd98OMHmrduntvlHFHuniQvhyRWkryMkAQAABBnep7YU8enHa9BCwZpV/Yut8s5LMfJHdzg2bOSGNzgaYQkAACAOJPgS9Dg1oO1YdcGjfpilNvlHJbPGEny7oQ7Bjd4GiEJAAAgDp1c8WRd0fAKvbriVX256Uu3yzmEL+9Tpmcn3NFu52mEJAAAgDh1S7NbVCW1iu6bf5/2Bfe5Xc5fOE7+ShIhCfGHkAQAABCnUgIpuq/Vffop8yeNXTrW7XL+Ir/dzvsrSexJ8iJCEgAAQBxrXa21LqpzkSYum6gV21a4Xc5+vryVJM+elWTyPkazkuRJMR+SjDENjTFjjDFvGWNucLseAACAWHNby9uUlpime+fdq5wY+dBu8laSbMjlQiKFlSRPi2hIMsZMMMZsMsYsO+jxjsaYlcaYVcaYO492DWvtd9baXpIuk9QikvUCAADEo7TENA08faC+2/adnv/2ebfLkST5cjOSd1eS2JPkaZFeSZokqeOBDxhjfJKelnSBpEaSLjfGNDLGNDXGTDnoq2Leay6UNFfS9AjXCwAAEJfa1Wyn8487X6O/Gq2fM392u5w/2+08vyeJkORFEQ1J1to5krYd9PCpklZZa1dba7MkvSbpImvtN9bafxz0tSnvOh9Ya1tLujKS9QIAAMSzgacPVKI/UffNv08hl/vcSs50O9rtvMiNPUnVJP16wJ/X5j12WMaYtsaYUcaYsZI+PMrzehhjFhtjFm/evDl81QIAAMSJ8snldXvL2/XFpi/05so3Xa3F+9PtHEmGlSSPciMkmcM8dsT/91hrZ1lrb7HW9rTWPn2U542z1raw1raoUKFCWAoFAACINxfVuUitqrTSyCUjtX7netfqcLzebiflriYRkjzJjZC0VlKNA/5cXdJvLtQBAADgOcYY3df6PllZPfDZA7IutbvlryR5tt1OIiR5mBshaZGkesaY440xCZI6S/rAhToAAAA8qVqpaurdrLc+Xfep/vfL/1ypwck/RsjDGUmOjz1JHhXpEeCvSlogqb4xZq0x5nprbY6kmyX9T9J3kt6w1n4byToAAABKmisaXKFGGY00bOEwZe7LjPr9Ha/vSZJyQ5IlJHlRpKfbXW6trWKtDVhrq1trn8t7/ENr7QnW2jrW2gfDdT9jTCdjzLjMzOj/RQAAABBLfI5P97e6X9v3bdfjXzzuwv1pt0P8cqPdLmKstZOttT3S0tLcLgUAAMB1DTMa6qqGV+mt79/SFxu/iOq9PT/dTiIkeZinQhIAAAD+6saTb1TV1KoatGCQsoJZUbsv0+0QzwhJAAAAHpYSSNHA0wdqdeZqTVg2IWr3LRnT7Rjc4FWEJAAAAI9rU72NOtbqqHFfj9NPmT9F5Z4+VpIQxwhJAAAAJcAdp96hJH+SBi8YHJWzk/IWkjw+ApyQ5FWeCklMtwMAADi88snl1f+U/lq8cbHeW/VexO9XIqbbGR8hyaM8FZKYbgcAAHBkl9S7RM0rNtcjix/R1j1bI3qvkjPdLuR2FYgAT4UkAAAAHJljHN3X6j7tztmtEYtGRPZe+StJng5JrCR5FSEJAACgBKmdXlvdm3bXhz99qHnr5kXsPvsHN3i53Y49SZ5FSAIAAChhujXtplplaumBzx7Qnpw9EbmHU2La7QhJXkRIAgAAKGESfAm6t9W9WrdznZ5Z+kxE7lEiBjc4fs5J8ihPhSSm2wEAABROy8otdUm9S/TCty9o5baVYb9+XkZS0MtzDdiT5FmeCklMtwMAACi8/qf0V5mEMhry2RCFbHjTTH67nfdXkghJXuSpkAQAAIDCS0tMU/8W/fXV5q/CfnaSj+l2iGN+twsAAACAey6qc5He/eFdjVwyUufUOEdlk8qG5bolZrrd1h+lly51u5LwM0Y6o69U6wy3K3EFIQkAAKAEM8bontPv0b8n/1uPLXlMg88YHJbrlojpdg3+Ie3YIO2O7MG8rlj/lVSuNiEJAAAAJVPdsnV1deOrNXHZRP2z3j/VrGKzYl+zREy3a3Zl7pcXDa9Voif3sScJAAAA6nViL1VJraLBCwYrO5Rd7Ov59q8kFftScIMp2futPBWSGAEOAABQNCmBFN156p1atX2VXl7+crGv5+R9yvT04AYvc/ySZSXJExgBDgAAUHTnHneu2lZvq9FLR2vDrg3FulaJGAHuZSX8oFxPhSQAAAAUz52n3SlrrYZ9PqxY1ykR0+28zHEISQAAAIAkVStVTT1P6qnpa6Zrzto5Rb7O/pUk2u3iUwk/KJeQBAAAgL+4ttG1qpNWR0MXDtWenD1Fusb+lSRCUnwiJAEAAAB/CvgCGnj6QK3buU7jvx5fpGvsn25HRopPxsfgBgAAAOBALSu31IV1LtTEbydq9fbVx/x6ptvFOQY3AAAAAIfqf0p/pfhTNGThENljHMDA4IY4x+AGAAAA4FAZyRnq07yPFm1YpCmrpxzTaxkBHufYk+QdHCYLAAAQXpeecKmalm+qkUtGamfWzkK/jul2cY7DZL2Dw2QBAADCyzGOBp42UFv3bNXopaML/bo/p9tFqjJElPGV6HY7v9sFAAAAILY1Lt9Y/zrhX3rlu1f0z7r/VL2y9Qp8TV5G0qYde7Viwx8RrtAdGamJqlA60e0yIsPxleh2O0ISAAAACtSnWR9N+2Wahi4cqgkdJsjktdMdiTFGKQk+vbxwjV5euCZKVUZXSoJPS+9rr4DPU81ZuRyflLPP7SpcQ0gCAABAgdKT0nVLs1v0wGcP6OOfP9YFx19Q4Gve6NlKv27bHYXqom/a8o1658t1yg6GPBqSSvbgBkISAAAACuVf9f6lt394W48sekRtqrdRaiD1qM9vUi1NTap5c6/42t/3SF+uU9CrgykY3AAAAAAUzOf4NPC0gdq0Z5PGLh3rdjmucpz8EecuFxIpJXxwAyEJAAAAhXZihRN1Sb1L9OLyF7V6+2q3y3FN/mAKz444L+GDGwhJAAAAOCZ9mvdRciBZQz8fKltCD4v1OR4/LNdhJQkAAAAotHJJ5dS7WW8tXL9QU3+Z6nY5rsif7hf0bEgq2YMbPBWSjDGdjDHjMjMz3S4FAADA0y474TI1KNdADy96WLuzvTnB7mh8eSHJqxmJwQ0eYq2dbK3tkZbmzSkqAAAAsSJ/iMPG3Rs17utxbpcTdfl7kjw73Y7BDQAAAMCxO7niybqwzoV6fvnz+jnzZ7fLiSqHPUmeRkgCAABAkfU7pZ+Sfcka9vmwEjXEIb/dLhRyuZBIYbodAAAAUDTlk8vrpmY3ad5v8zRjzQy3y4kaJ+9TNIMbvImQBAAAgGL5T/3/6ISyJ2jEohHam7PX7XKiwjFeb7fzS9ary2QFIyQBAACgWPyOX3eeeqd+2/WbJn470e1yomJ/SPLs4AaHlSQAAACgOFpWbqn2NdtrwjcTtH7nerfLibg/D5N1uZBIcfwMbgAAAACKa0CLAbKyGrlkpNulRJznR4AzuAEAAAAovqqlqqprk676+OePtXjDYrfLiaiSsScp6OHTco+OkAQAAICw6dKkiyqnVtawz4cp6OF2rRIRkqQSO7yBkAQAAICwSfYna0CLAVr5+0q9s+odt8uJGM/vSTJ5MaGEttwRkgAAABBWHWp2UItKLfTkF08qc1+m2+VEhPH8nqS8lSQPrwYeDSEJAAAAYWWM0Z2n3qnMrEyNWTrG7XIiIn8lyXq23c6X+52VpPhnjOlkjBmXmenN31gAAADEi/rl6uvSepfq1RWv6sftP7pdTtjl70ny/koSISnuWWsnW2t7pKWluV0KAABAiXdzs5uVEkjR8M+He27F5c/BDS4XEikMbgAAAADCr2xSWd108k1asH6BZv460+1ywurPwQ0eTUkMbgAAAAAi47L6l6lOWh09vOhh7Qvuc7ucsPH+YbIMbgAAAAAiIuAEdMepd2jtzrV6cfmLbpcTNo7XV5IY3AAAAABETquqrXRujXM17utx2rhro9vlhEXJOUyWlSQAAAAgIm5teauCoaAe/+Jxt0sJC19+SPLqXAPa7QAAAIDIqlG6hq5tfK2mrJ6ipZuXul1Ose0/TNarK0kMbgAAAAAir1vTbqqQXEEjFo2I+5Hg3j9MlpUkAAAAIOJSAinq3ay3vt78tT766SO3yymWPw+TdbmQSOEwWQAAACA6Lqp7kRqWa6jHvnhMe3L2uF1Okfnyu9E8u5KUN92OwQ0AAABAZDnG0W0tb9OGXRv0wrcvuF1OkRnPT7fLHwFOSAIAAAAirmXlljr/uPP13LLntGn3JrfLKRKf10OS4ZwkAAAAIKr6n9JfOaEcPfnlk26XUiQlZ08SK0kAAABAVNQoU0NXNrxS7696X8u3Lne7nGPm5O9JCnl0JYnBDQAAAED09Tixh9IT0/XwoofjbpR2/ghwz7bbMbgBAAAAiL7SCaV108k3afHGxZq+Zrrb5RyT/e12Xg9JtNsBAAAA0fWvE/6luul19ejiR5UVzHK7nEJz9g9ucLmQSCnhgxv8bhcAAACAksvv+HVbi9vU85OeeuW7V3Rdk+vcLqlQ8rrtvL8nad8OaW+mu7VEgj/56D+OUhlRYYzpJKlT3bp13S4FAAAAhdS6WmudVe0sjf16rC6se6HKJZVzu6QCeX5Pkj8p9/u7Pd2tI1KSj/6/MU+FJGvtZEmTW7Ro0d3tWgAAAFB4t7a4VZd8cIlGfzVad59+t9vlFMjsHwHu0ZCUUUe6eIy0Z5vblYTfms+k7z446lM8FZIAAAAQn2qn19Zl9S/T6ytfV+f6nVW3bGx3BuWvJHl1IUnGSCdf7nYVkRFIKTAkMbgBAAAAMeHGk25UaiBVDy+O/ZHg+XuSPDvdzsvyJ/cd7SlRKAMAAAAoUHpSunqd2Evzf5uvT9d96nY5R/XndDtCUtwxhCQAAADEkcsbXK6aZWpq5OKRyonh8dP7Q5JX9yR5mVPwjiNCEgAAAGJGwBdQ3+Z99WPmj3p31btul3NEf063c7kQHDva7QAAABBvzjvuPDWr2ExPf/m0dmfvdrucw9q/J4mUFH8ISQAAAIg3xhgNaDFAW/du1cRvJ7pdzmEZY2QMe5LiEnuSAAAAEI9OqnCSOtTqoOe/fV6bdm9yu5zD8hlDSIpHrCQBAAAgXvVp3kfZoWw9/dXTbpdyWI4xCobcrgLHjMENAAAAiFc1StdQ5/qd9d6q9/T979+7Xc4hHEcxf54TDoN2OwAAAMSznif2VGogVSOXjHS7lEPkriQRkuKOU3AEIiQBAAAgZqUnpatH0x6at26e5v823+1y/iJ3T5LbVeCY0W4HAACAeHd5w8tVrVQ1jVw8UsFQ0O1y9mO6XZyi3Q4AAADxLtGXqFua3aKVv6/UlNVT3C5nP5/DdLu4xHQ7AAAAeEHH4zuqSUYTjfpylPbk7HG7HEnsSYpbtNsBAADACxzjaECLAdq0e5NeWv6S2+VIkhyHPUlxyTC4AQAAAB7RonILta3RVs8te05b92x1uxw5RgqRkuIP7XYAAADwkn6n9NPenL16ZukzbpeSN92OkBR3aLcDAACAl9ROq61LT7hUb33/ln7K/MnVWowxChKS4g/T7QAAAOA1N5x0g5L8SXpsyWOu1uFzDO128YiVJAAAAHhNRnKGujbpqpm/ztTiDYtdq8PH4Ib45DC4AQAAAB50daOrVTGloh5d/KisSy1vxoh2u3hEux0AAAC8KNmfrJtPvlnLti7TtF+muVKDzxjXAhqKgXY7AAAAeFWnOp1UJ62OnvzySWWHsqN+fw6TjVOMAAcAAIBX+R2/+jTvo5//+Fnv/vBu1O/PYbJxinY7AAAAeFnbGm3VrGIzPbP0Ge3O3h3Ve3OYbJxiJQkAAABeZoxRv1P6acueLXrpu5eieu/c6XaEpLjjlZBkjEk1xiwxxvzD7VoAAAAQW5pVbKZzapyjicsm6ve9v0ftvrmHyUbtdggXt9vtjDETjDGbjDHLDnq8ozFmpTFmlTHmzkJc6g5Jb0SmSgAAAMS7Ps37aHfObo3/ZnzU7ukzYrpdPIqB6XaTJHU88AFjjE/S05IukNRI0uXGmEbGmKbGmCkHfVU0xpwvabmkjRGuFQAAAHGqTnodXVz3Yr224jWt27kuKvdkul2ccrvdzlo7R9K2gx4+VdIqa+1qa22WpNckXWSt/cZa+4+DvjZJOkfS6ZKukNTdGBMXLYIAAACIrhtOukGOcfT0l09H5X6OQ0iKS4Votyt4rSn8qkn69YA/r5V02pGebK0dKEnGmOskbbHWhg73PGNMD0k98v640xizMsDl4bwAAAWDSURBVCzVFk55SVuieD+ED+9dfOP9i2+8f/GL9y6+ef79W6IlekgPRe1+b/SK2q08/95FWc0j/cCNkGQO81iBEdxaO6mAn4+TNK6INRWLMWaxtbaFG/dG8fDexTfev/jG+xe/eO/iG+9f/OK9ix43WtfWSqpxwJ+rS/rNhToAAAAA4BBuhKRFkuoZY443xiRI6izpAxfqAAAAAIBDRHoE+KuSFkiqb4xZa4y53lqbI+lmSf+T9J2kN6y130ayjihwpc0PYcF7F994/+Ib71/84r2Lb7x/8Yv3LkoMs90BAAAA4E+M0wYAAADw/+3dW8gc5R3H8e+PRIlGJXqhaCIYQTwQ6gFbPEAvtIF4wHgrVgN6qdUWwQPeFhUspYJiKZ4CjYpExSB4CCp40VaCtkRTtYoWfUs0gnhAEbX+vZiRDjbR3ZX47G6+H3jZmbn6wcPOO7995pnRgCXpB0iyJsmrSV5Pcm3rPBpdksOTPJPk5STbklzZOpPGk2RRkr8nebR1Fo0nybIkG5O80n8HT22dSaNL8pv+vPlSkvuSLGmdSbuW5K4kO5K8NDh2UJLNSV7rPw9smVE7t4uxu7k/d25N8nCSZS0zzjNL0oSSLAJuA84CjgMuSHJc21Qaw5fAVVV1LN3Lii9z/GbOlXTrGjV7bgEer6pjgONxHGdGkuXAFcDJVbUKWET3ACZNr3uANd86di3wVFUdBTzV72v63MP/j91mYFVV/QT4F3Ddjx1qT2FJmtzPgNer6o2q+hy4H1jbOJNGVFXbq+qFfvtjuou05W1TaVRJVgDnAHe0zqLxJDkA+DlwJ0BVfV5VH7RNpTEtBvZJshjYF1/jMdWq6lng/W8dXgus77fXA+f/qKE0kp2NXVU92T8EDeBvdK/S0W5gSZrccuDtwf4CXmTPpCRHACcCz7VNojH8Abga+Kp1EI3tSOA94O7+dsk7kixtHUqjqar/AL8D3gK2Ax9W1ZNtU2kCh1TVduh+NAQObpxHk7kEeKx1iHllSZpcdnLMRwXOmCT7AQ8Cv66qj1rn0fdLci6wo6qeb51FE1kMnATcXlUnAp/grT4zo1+7shZYCRwGLE3yy7appD1Pkuvplg5saJ1lXlmSJrcAHD7YX4G3HMyUJHvRFaQNVfVQ6zwa2enAeUn+TXeb6xlJ/tw2ksawACxU1TcztxvpSpNmwy+AN6vqvar6AngIOK1xJo3v3SSHAvSfOxrn0RiSrAPOBS4s3+Wz21iSJrcFOCrJyiR70y1c3dQ4k0aUJHRrIl6uqt+3zqPRVdV1VbWiqo6g+949XVX+kj0jquod4O0kR/eHzgT+2TCSxvMWcEqSffvz6Jn44I1ZtAlY12+vAx5pmEVjSLIGuAY4r6o+bZ1nnlmSJtQvmrsceILuH8QDVbWtbSqN4XTgIrpZiH/0f2e3DiXtIX4FbEiyFTgBuKFxHo2onwHcCLwAvEh3HfGnpqH0nZLcB/wVODrJQpJLgZuA1UleA1b3+5oyuxi7W4H9gc39tcsfm4acY3GWTpIkSZL+x5kkSZIkSRqwJEmSJEnSgCVJkiRJkgYsSZIkSZI0YEmSJEmSpAFLkiRJkiQNWJIkSZIkacCSJEmaW0l+mmRrkiVJlibZlmRV61ySpOnmy2QlSXMtyW+BJcA+wEJV3dg4kiRpylmSJElzLcnewBbgM+C0qvpv40iSpCnn7XaSpHl3ELAfsD/djJIkSd/JmSRJ0lxLsgm4H1gJHFpVlzeOJEmacotbB5AkaXdJcjHwZVXdm2QR8JckZ1TV062zSZKmlzNJkiRJkjTgmiRJkiRJGrAkSZIkSdKAJUmSJEmSBixJkiRJkjRgSZIkSZKkAUuSJEmSJA1YkiRJkiRpwJIkSZIkSQNfA78DIP8zw5lfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig0, ax0 = plt.subplots(figsize=(14, 8))\n", "ax0.plot(xaxis, yaxis_binom, '-', label=f'Binomial with n={n:2d} and p={p:.3f}')\n", "ax0.plot(xaxis, yaxis_poiss, '-', label=f'Poisson with lambda={Lambda:.2f}')\n", "ax0.plot(xaxis, yaxis_gauss, '-', label=f'Gaussian with mu={mu:.2f} and sigma={sigma:.2f}')\n", "ax0.set(xlim=(xmin, xmax),\n", " title='Probability for Binomial and Gaussian', \n", " xlabel='x', \n", " ylabel='Probability')\n", "ax0.set_yscale('log')\n", "ax0.set_ylim(1e-4, 1.0)\n", "ax0.legend(loc='upper right');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Looping over processes:\n", "In the following we simulate a Binomial/Poisson process with given parameters, i.e. number of trials and probability of success. For the Poisson, these can not be specified, but the resulting expected number is naturally lambda = n * p.\n", "\n", "After having simulated the process, we fit the result with the three distributions in question, and test to what extend they match." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " With N_trials = 20 and p_success = 0.2000, the average number of successes is lambda = 4.00\n" ] } ], "source": [ "# Simulation parameters:\n", "N_experiments = 1000 # Number of simulations/experiments to perform\n", "\n", "N_trials = n # Number of trials in each experiment (taken from above!)\n", "p_success = p # Chance of succes in each trial (taken from above!)\n", "Lambda = N_trials * p_success # This is the mean and the one parameter by which the Poisson is defined!\n", "\n", "print(f\" With N_trials = {N_trials:d} and p_success = {p_success:.4f}, the average number of successes is lambda = {Lambda:.2f}\")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n_success: 6\n", "n_success: 6\n", "n_success: 6\n", "n_success: 6\n", "n_success: 4\n", "n_success: 4\n", "n_success: 3\n", "n_success: 3\n", "n_success: 6\n", "n_success: 2\n" ] } ], "source": [ "all_n_success = np.zeros(N_experiments)\n", "\n", "# Run the experiments, and fill the histogram from above:\n", "for iexp in range(N_experiments): \n", " \n", " # Simulating process defined:\n", " n_success = 0\n", " for i in range(N_trials): \n", " x = r.uniform()\n", " if (x < p_success): \n", " n_success += 1\n", "\n", " # Record result:\n", " if (verbose and iexp < N_verbose): \n", " print(f\"n_success: {n_success:4d}\")\n", " \n", " # Save Result\n", " all_n_success[iexp] = n_success" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot result:\n", "\n", "Define a histogram with the \"data\" (note and think about the binning!). Also, ask yourself what uncertainty to assign to each bin." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "counts, bin_edges = np.histogram(all_n_success, bins=N_trials+1, range=(-0.5, N_trials+0.5))\n", "bin_centers = (bin_edges[1:] + bin_edges[:-1])/2\n", "s_counts = np.sqrt(counts) # NOTE: We (naturally) assume that the bin count is Poisson distributed." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# We remove any bins, which don't have any counts in them:\n", "x = bin_centers[counts>0]\n", "y = counts[counts>0]\n", "sy = s_counts[counts>0]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAHgCAYAAABn17aGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZRtaV0f+O+PbpTixQDTDdM2YCPTAdGlDSlZ3EsS0BJF4tiQFRx6DBCDaaOAkGhGUNfUcdYiC6Pgy9LB8Da8BHFaXqRXgkJDCETrIFQjgYYW6QEGLt1DNyEBRAUbfvNHnZLqS/W95zb3nFNVz+ez1ln77Gfvc/bvnH1v3fre59nPru4OAADAaG636gIAAABWQRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBI5666gK/Feeed1xdddNGqywAAAA6oq6+++lPdff5+2w51GLrooouyvb296jIAAIADqqr+31vbZpgcAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADGlhYaiq7l1Vb62qa6vq/VX19Fn7pKo+UVXvmT0evec1z6qq66rqg1X1fYuqDQAA4NwFvvfNSX6qu99dVXdJcnVVXTXb9ivd/ct7d66qByZ5fJJvTfKNSd5cVX+7u7+0wBoBAIBBLaxnqLtv6O53z55/Lsm1SS48xUsuTfI73f2F7v5IkuuSPGRR9QEAAGNbyjVDVXVRkgcl+eNZ01Or6r1V9ZKqutus7cIkH9/zshM5dXgCAAC4zRYehqrqzklek+QZ3f3ZJM9Pcr8klyS5Iclzd3fd5+W9z/tdXlXbVbV90003LahqAADgqFtoGKqq22cnCL2yu1+bJN39ye7+Und/OckL85WhcCeS3HvPy++V5PqT37O7X9Dd6929fv755y+yfAAA4Ahb5GxyleTFSa7t7uftab9gz26PTXLN7PmVSR5fVV9fVfdNcnGSdy6qPgAAYGyLnE3uYUmekOR9VfWeWdvPJrmsqi7JzhC4jyb5sSTp7vdX1RVJPpCdmeieYiY5AABgURYWhrr7D7P/dUBvOMVrnp3k2YuqCQAAYNdSZpMDAAA4aIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEHCgTSaTVFUmk8mqSwEAjpjq7lXXcJutr6/39vb2qssAFqyqcph/VgEAq1NVV3f3+n7b9AwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDMEhN5lMUlWZTCarLgUA4FCpw3xX9/X19d7e3l51GbByVZXD/Hf5dI765wMAFqeqru7u9f226RkCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShoADbTqd3mIJAHC2CEPAGZtMJks5znQ6zcbGRpJkY2NjaYFoWZ8PAFit6u5V13Cbra+v9/b29qrLgJWrqizz73JVLe1Yq3KYfzYCAF9RVVd39/p+2/QMAWdsc3Mz3b3wx9bWVtbW1pIka2tr2draWspxNzc3V/wNAwDLoGcIjoBl9wwt03Q6zfHjx7O1tZVjx46tuhwA4JDRMwQcWrsBSBACAM42YQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShuCQ270R6bJuSAoAcFQIQ7Agk8lk4ceYTqfZ2NhIkmxsbCwtEC3jswEALJr7DMGCVNWqS1ioZf7sOMr3UQIAFst9hmAFNjc3090LfWxtbWVtbS1Jsra2lq2trYUfs7uzubm54m8XAOBrp2cIDrnpdJrjx49na2vryN6YVM8QAHBb6RmCI2w3AB3VIAQAsCjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCDrTJZHKL5VEzmUxSVUf28wHAQVbdveoabrP19fXe3t5edRmwclWVw/x3eXTOHwAsTlVd3d3r+23TMwQAAAxJGAIAAIYkDAEAAEMShgAAgCEtLAxV1b2r6q1VdW1Vvb+qnj5rv3tVXVVVH5ot7zZrr6r69aq6rqreW1UPXlRtAAAAi+wZujnJT3X3tyR5aJKnVNUDkzwzyVu6++Ikb5mtJ8n3J7l49rg8yfMXWBsAADC4hYWh7r6hu989e/65JNcmuTDJpUleNtvtZUkeM3t+aZKX9453JLlrVV2wqPoAAICxLeWaoaq6KMmDkvxxknt29w3JTmBKco/Zbhcm+fiel52YtQGncNRvSgoAsCgLD0NVdeckr0nyjO7+7Kl23aftq+5CWFWXV9V2VW3fdNNNZ6tMOLQmk0m6WxgCADhDCw1DVXX77AShV3b3a2fNn9wd/jZb3jhrP5Hk3ntefq8k15/8nt39gu5e7+71888/f3HFAwAAR9oiZ5OrJC9Ocm13P2/PpiuTPGn2/ElJXr+n/YmzWeUemuQzu8PpAAAAzrZzF/jeD0vyhCTvq6r3zNp+NslzklxRVU9O8rEkj5tte0OSRye5LslfJPmRBdYGAAAMbmFhqLv/MPtfB5QkG/vs30mesqh6AAAA9lrKbHIAAAAHjTAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCGAFZpOp7dYAgDLIwwBnGQymSzlONPpNBsbO7dd29jYWFogWtbnA4CDrnbudXo4ra+v9/b29qrLAI6Yqlu7X/TRcZh/9gPAmaiqq7t7fb9teoYATrK5uZnuXvhja2sra2trSZK1tbVsbW0t5bibm5sr/oYB4GDQMwSwQtPpNMePH8/W1laOHTu26nIA4MjRMwRwQO0GIEEIAJZPGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAKzSZTG6xBACWp7p71TXcZuvr6729vb3qMgAAgAOqqq7u7vX9tukZAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxpYWGoql5SVTdW1TV72iZV9Ymqes/s8eg9255VVddV1Qer6vsWVRcAAECy2J6hlyZ51D7tv9Ldl8web0iSqnpgkscn+dbZa/7PqjpngbUBAACDW1gY6u63J/n0nLtfmuR3uvsL3f2RJNcleciiagMAAFjFNUNPrar3zobR3W3WdmGSj+/Z58SsDQAAYCGWHYaen+R+SS5JckOS587aa599e783qKrLq2q7qrZvuummxVQJAAAceUsNQ939ye7+Und/OckL85WhcCeS3HvPrvdKcv2tvMcLunu9u9fPP//8xRYMAAAcWUsNQ1V1wZ7VxybZnWnuyiSPr6qvr6r7Jrk4yTuXWRsAADCWcxf1xlX1qiSPSHJeVZ1IspnkEVV1SXaGwH00yY8lSXe/v6quSPKBJDcneUp3f2lRtQEAAFT3vpfm3HKnqm/r7mtOu+OSra+v9/b29qrLAAAADqiqurq71/fbNu8wud+qqndW1U9U1V3PYm0AAAArMVcY6u6/m+SHszPJwXZV/XZVPXKhlQEAACzQ3BModPeHkvx8kp9J8vAkv15Vf1pV/3BRxcHZMplMUlWZTCarLgUAgANi3muGvj3JjyT5B0muSvLi7n53VX1jkml3f9Niy9yfa4Y4E1WVef68AwBwdJzqmqF5Z5P7jezcF+hnu/svdxu7+/qq+vmzUCMAAMBSzRuGHp3kL3enu66q2yW5Q3f/RXe/YmHVAQAALMi81wy9OcnanvU7ztoAAAAOpXnD0B26+893V2bP77iYkgAAABZv3jD0+ap68O5KVf2dJH95iv0BAAAOtHmvGXpGkt+tqutn6xck+V8WUxIAAMDizRWGuvtdVfWAJPdPUkn+tLv/eqGVAQAALNC8PUNJ8p1JLpq95kGze7a8fCFVAQAALNhcYaiqXpHkfknek+RLs+ZOIgwBAACH0rw9Q+tJHtjdvchiAAAAlmXe2eSuSfI/LrIQAACAZZq3Z+i8JB+oqncm+cJuY3f/4EKqAgAAWLB5w9BkkUUAAAAs27xTa7+tqr4pycXd/eaqumOScxZbGgAAwOLMdc1QVf2zJK9O8m9nTRcm+b1FFQUAALBo806g8JQkD0vy2STp7g8luceiioKzbTqd3mIJAADzhqEvdPcXd1eq6tzs3GcIbrPJZLKU40yn02xsbCRJNjY2lhaIlvX5AAC4bWqeWwdV1b9J8t+TPDHJ05L8RJIPdPfPLba8U1tfX+/t7e1VlsDXoKpWXcLCuTUXAMBqVdXV3b2+37Z5e4aemeSmJO9L8mNJ3pDk589OeYxqc3Mz3b3wx9bWVtbW1pIka2tr2draWspxNzc3V/wNAwBwKnP1DB1UeoaY13Q6zfHjx7O1tZVjx46tuhwAAJbkVD1Dc02tXVUfyT7XCHX3N3+NtcFS7AYgQQgAgF3z3nR1b5K6Q5LHJbn72S8HAABgOea6Zqi7/+uexye6+1eTfPeCawMAAFiYeYfJPXjP6u2y01N0l4VUBAAAsATzDpN77p7nNyf5aJIfOuvVAAAALMlcYai7v2vRhQAAACzTvMPk/uWptnf3885OOQAAAMtxJrPJfWeSK2fr/3OStyf5+CKKAgAAWLR5w9B5SR7c3Z9LkqqaJPnd7v7RRRUGAACwSHNNrZ3kPkm+uGf9i0kuOuvVAAAALMm8PUOvSPLOqnpdkk7y2CQvX1hVAAAACzbvbHLPrqrfT/L3Zk0/0t1/sriyAAAAFmveYXJJcsckn+3uX0tyoqruu6CaAAAAFm6uMFRVm0l+JsmzZk23T/LvFlUUAADAos3bM/TYJD+Y5PNJ0t3XJ7nLoooCAABYtHnD0Be7u7MzeUKq6k6LKwkAAGDx5g1DV1TVv01y16r6Z0nenOSFiysLzq7JZHKLJQAA1E6Hzxw7Vj0yyfcmqSRv7O6rFlnYPNbX13t7e3vVZQAAAAdUVV3d3ev7bTttz1BVnVNVb+7uq7r7X3X3Tx+EIATAwTeZTFJVemUBOJDm6hmqqiuTPKG7P7P4kuanZwjg4KuqzDsKAQDOtlP1DM1109Ukf5XkfVV1VWYzyiVJd//kWagPAABg6eYNQ/9h9gAAADgSThmGquo+3f2x7n7ZsgoCAABYhtNNoPB7u0+q6jULrgUAAGBpTheGas/zb15kIQAAAMt0ujDUt/IcAADgUDvdBArfUVWfzU4P0drseWbr3d3fsNDqAAAAFuSUYai7z1lWIQAAAMt0umFyAAAAR5IwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDACzMdDq9xRIADhJhCGAwk8lkKceZTqfZ2NhIkmxsbCwtEC3r8wFw+FV3r7qG22x9fb23t7dXXQbAoVJVqy5h4Q7zv20AnF1VdXV3r++3Tc8QwGA2NzfT3Qt/bG1tZW1tLUmytraWra2tpRx3c3Nzxd8wAIeFniEAFmY6neb48ePZ2trKsWPHVl0OAAPSMwTASuwGIEEIgINIGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQ1pYGKqql1TVjVV1zZ62u1fVVVX1odnybrP2qqpfr6rrquq9VfXgRdUFAACQLLZn6KVJHnVS2zOTvKW7L07yltl6knx/kotnj8uTPH+BdQEAACwuDHX325N8+qTmS5O8bPb8ZUkes6f95b3jHUnuWlUXLKo2AACAZV8zdM/uviFJZst7zNovTPLxPfudmLUBcIhNJpNbLAHgIDkoEyjUPm29745Vl1fVdlVt33TTTQsuC4CvxWQySXcLQwAcSMsOQ5/cHf42W944az+R5N579rtXkuv3e4PufkF3r3f3+vnnn7/QYgEAgKNr2WHoyiRPmj1/UpLX72l/4mxWuYcm+czucDoAAIBFOHdRb1xVr0ryiCTnVdWJJJtJnpPkiqp6cpKPJXncbPc3JHl0kuuS/EWSH1lUXQAAAMkCw1B3X3Yrmzb22beTPGVRtQAAAJzsoEygAAAAsFTCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIZIkk8kkVZXJZLLqUgAAYCmqu5d/0KqPJvlcki8lubm716vq7kn+7yQXJflokh/q7v92qvdZX1/v7e3txRY7kKrKKv48AADAolTV1d29vt+2VfYMfVd3X7KnsGcmeUt3X5zkLbN1ADiw9KoDHG6r7Bla7+5P7Wn7YJJHdPcNVXVBkv/U3fc/1fvoGTq79AwBnDk/OwEOtoPYM9RJ3lRVV1fV5bO2e3b3DUkyW95jRbUBAAADOHdFx31Yd19fVfdIclVV/em8L5yFp8uT5D73uc+i6gMAAI64lfQMdff1s+WNSV6X5CFJPjkbHpfZ8sZbee0Lunu9u9fPP//8ZZUMAAAcMUsPQ1V1p6q6y+7zJN+b5JokVyZ50my3JyV5/bJrAwAAxrGKYXL3TPK6qto9/m939x9U1buSXFFVT07ysSSPW0FtAADAIJYehrr7w0m+Y5/2/5pkY9n1sGM6nf7N8tixYyuuBgAAFm+V9xniNJZ134rpdJqNjZ0curGx8TfBaNHclwMAgFVayX2Gzpajfp+h2VDCI+0w//kDmE6nOX78eLa2tvSqAxxQB/E+Q8xhc3Mz3b3wx9bWVtbW1pIka2tr2draWspxNzc3V/wNA0eRXnUA5qVniCT+dxM4OvSqA7CXniFOazcACULAYadXHYB56Rnib1SV/20EOAN61QEOPj1DALAAetUBDjdhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDJPnKHc3d2RwAgFG4zxAAfA3cow3gYHOfIQBYAL3qAIebniEAAODI0jMEAABwEmEIAAAYkjAEAAAMSRgCAACGJAwBAABDEoYAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEMShgAAgCEJQwAAwJCEIQAAYEjCEAAAMCRhCAAAGJIwBAAADEkYAgAAhiQMAQAAQxKGAACAIQlDAMC+JpNJqiqTyWTVpQAsRHX3qmu4zdbX13t7e3vVZQDAkVVVOcy/KwBU1dXdvb7fNj1DAADAkIQhAABgSMIQAAAwJGEIAAAYkjAEAOxrOp3eYglw1AhDAHCILGua6+l0mo2NjSTJxsbGUgKRKbyBZTO1NgAcIlW16hIW6jD/XgIcTKbWBoAjYnNzM9298MfW1lbW1taSJGtra9na2lr4MTc3N1f87QKj0TMEAOxrOp3m+PHj2drayrFjx1ZdDsBtomcIADhjuwFIEAKOKmEIAAAYkjAEAAAMSRgCAACGJAwBAPvave+P+/8AR5XZ5AAAgCPLbHIAAAAnEYYAAIAhCUNzmkwmqSrjpgEA4IhwzdAZqKoc5u8LAABG45ohAACAkwhDc5pOp7dYAgAAh9uhDkPXX3/9Uo4znU6zsbGRJNnY2FhaIHJ9EgAALM6hvmaoqg5v8XM6zOcHAABW7cheM3TBBRekuxf+2NraytraWpJkbW0tW1tbSznu5ubmir9hADiazBILJIe8Z2iZs8lNp9McP348W1tbOXbs2FKOCQAsjlliYQxHtmdomXYDkCAEAIefiZGARBgCAA6IZQ1ZMzESsMswuTlNJpP8wi/8QjY3N/0wA4AFqKpVl7Bwh/n3LjisDJM7CyaTSbpbEAKABdnc3DQxErBUeoYAgOGYGAnGoWcIAGAPEyMBiTAEAAxod9j7UR3+7j5KMB/D5AAAjhjDAOErDtUwuap6VFV9sKquq6pnrroeAICzwdThcPAcqJ6hqjonyZ8leWSSE0neleSy7v7AfvvrGQIADgtTh8NqHKaeoYckua67P9zdX0zyO0kuXXFNAABfM1OHw8Fz0MLQhUk+vmf9xKwNAOBQW9YwsmPHjuWyyy5Lklx22WVLu2bIMDkOo4M2TO5xSb6vu390tv6EJA/p7qft2efyJJfPVu+f5INLLPG8JJ9a4vE4u5y/w8u5O9ycv8PN+Tu8nLvDzfk7e76pu8/fb8O5y67kNE4kufee9XsluX7vDt39giQvWGZRu6pq+9bGG3LwOX+Hl3N3uDl/h5vzd3g5d4eb87ccB22Y3LuSXFxV962qr0vy+CRXrrgmAADgCDpQPUPdfXNVPTXJG5Ock+Ql3f3+FZcFAAAcQQcqDCVJd78hyRtWXcetWMnwPM4a5+/wcu4ON+fvcHP+Di/n7nBz/pbgQE2gAAAAsCwH7ZohAACApRCG5lRVj6qqD1bVdVX1zFXXw3yq6t5V9daquraq3l9VT191TZy5qjqnqv6kqv79qmthflV116p6dVX96ezv4HJudsJZUVX/YvZz85qqelVV3WHVNXHrquolVXVjVV2zp+3uVXVVVX1otrzbKmtkf7dy7n5p9rPzvVX1uqq66yprPMqEoTlU1TlJfjPJ9yd5YJLLquqBq62KOd2c5Ke6+1uSPDTJU5y7Q+npSa5ddRGcsV9L8gfd/YAk3xHn8NCoqguT/GSS9e7+tuxMavT41VbFabw0yaNOantmkrd098VJ3jJb5+B5ab763F2V5Nu6+9uT/FmSZy27qFEIQ/N5SJLruvvD3f3FJL+T5NIV18QcuvuG7n737PnnsvPL2IWrrYozUVX3SvIPkrxo1bUwv6r6hiR/P8mLk6S7v9jd/321VXGGzk2yVlXnJrljTrrvHwdLd789yadPar40yctmz1+W5DFLLYq57HfuuvtN3X3zbPUd2bn3JgsgDM3nwiQf37N+In6hPnSq6qIkD0ryx6uthDP0q0n+tyRfXnUhnJFvTnJTkv9rNsTxRVV1p1UXxXy6+xNJfjnJx5LckOQz3f2m1VbFbXDP7r4h2fnPwST3WHE93Db/NMnvr7qIo0oYmk/t02YavkOkqu6c5DVJntHdn111Pcynqn4gyY3dffWqa+GMnZvkwUme390PSvL5GKJzaMyuLbk0yX2TfGOSO1XVP15tVTCeqvq57Az5f+WqazmqhKH5nEhy7z3r94rhAodGVd0+O0Hold392lXXwxl5WJIfrKqPZmd46ndX1b9bbUnM6USSE9292xP76uyEIw6H70nyke6+qbv/OslrkxxfcU2cuU9W1QVJMlveuOJ6OANV9aQkP5Dkh9u9cBZGGJrPu5JcXFX3raqvy85FpFeuuCbmUFWVnWsWru3u5626Hs5Mdz+ru+/V3Rdl5+/df+xu/zt9CHT3/5fk41V1/1nTRpIPrLAkzszHkjy0qu44+zm6ERNgHEZXJnnS7PmTkrx+hbVwBqrqUUl+JskPdvdfrLqeo0wYmsPsAranJnljdv4xuKK737/aqpjTw5I8ITs9Cu+ZPR696qJgEE9L8sqqem+SS5L86xXXw5xmPXqvTvLuJO/Lzu8LL1hpUZxSVb0qyTTJ/avqRFU9Oclzkjyyqj6U5JGzdQ6YWzl3v5HkLkmumv3u8lsrLfIIK71uAADAiPQMAQAAQxKGAACAIQlDAADAkIQhAABgSMIQAAAwJGEIYHBV1VX13D3rP11Vk7P03i+tqn90Nt7rNMd5XFVdW1VvXfSxADg6hCEAvpDkH1bVeasuZK+qOucMdn9ykp/o7u9aVD0AHD3CEAA3Z+eGmv/i5A0n9+xU1Z/Plo+oqrdV1RVV9WdV9Zyq+uGqemdVva+q7rfnbb6nqv7zbL8fmL3+nKr6pap6V1W9t6p+bM/7vrWqfjs7N/s8uZ7LZu9/TVX94qztf0/yd5P8VlX90kn7X1BVb5/dtPCaqvp7ez/H7Pk/qqqXzp7fs6peV1X/ZfY4Pmt/4qzO/1JVr5i1nV9Vr5l9hndV1cNm7Q/fc5PnP6mqu5yiju+tqmlVvbuqfreq7jxrf05VfWB2zF8+g3MJwBk4d9UFAHAg/GaS91bVvzmD13xHkm9J8ukkH07you5+SFU9PcnTkjxjtt9FSR6e5H5J3lpV/1OSJyb5THd/Z1V9fZI/qqo3zfZ/SJJv6+6P7D1YVX1jkl9M8neS/Lckb6qqx3T3/1FV353kp4PpZ4kAAANlSURBVLt7+6Qa/9ckb+zuZ896mu54ms/060ne1t2Pne1/56r61iQ/l+Rh3f2pqrr7bN9fS/Ir3f2HVXWfJG+cfR8/neQp3f1Hs3DzV0kuP7mOWU/czyf5nu7+fFX9TJJ/WVW/keSxSR7Q3V1Vdz1NzQDcRsIQAOnuz1bVy5P8ZJK/nPNl7+ruG5Kkqv6fJLth5n1J9g5Xu6K7v5zkQ1X14SQPSPK9Sb59T6/T30pycZIvJnnnyUFo5juT/Kfuvml2zFcm+ftJfu9UNSZ5SVXdPsnvdfd7TvOZvjs7QS3d/aUkn6mqJyZ5dXd/atb+6dm+35PkgVW1+9pvqKq7JPmjJM+b1ffa7j5RVV9VR1U9PMkDsxMEk+TrkkyTfDY7AepFVfUfkvz709QMwG1kmBwAu341O9fe3GlP282Z/VtRO7+xf92ebV/Y8/zLe9a/nFv+Z1ufdJxOUkme1t2XzB737e7dMPX5W6mvbqX9VnX327MTmD6R5BWzYHNyTXc4zdtUvvozJDvfy7E9n+HC7v5cdz8nyY8mWUvyjqp6wK3UUUmu2vP6B3b3k7v75uz0jr0myWOS/MGZfm4A5iMMAZDkb3o8rshOINr10ewMS0uSS5Pc/ja89eOq6naz64i+OckHszOk7MdnPSWpqr9dVXc61Zsk+eMkD6+q82ZDzS5L8rZTvaCqvinJjd39wiQvTvLg2aZPVtW3VNXtsjMkbddbkvz47LXnVNU3zNp+qKr+h1n77jC5NyV56p5jXTJb3q+739fdv5hkO8kDbqWOdyR52GzYYKrqjrPv4c5J/lZ3vyE7Qw0vOc33AsBtZJgcAHs9N3t+wU/ywiSvr6p3ZicU3Fqvzal8MDuh5Z5J/nl3/1VVvSg71xK9e9bjdFN2ekFuVXffUFXPSvLW7PSqvKG7X3+aYz8iyb+qqr9O8ueZDYFL8szsDD/7eJJrktx51v70JC+oqicn+VKSH+/uaVU9O8nbqupLSf4kyT/JzpDC36yq92bn39O3J/nnSZ5RVd81e/0Hkvx+ksefXEd331RV/yTJq2bXTSU71xB9Ljvf+R1mn/OrJrYA4Oyo7v16/gEAAI42w+QAAIAhCUMAAMCQhCEAAGBIwhAAADAkYQgAABiSMAQAAAxJGAIAAIYkDAEAAEP6/wGHYebl+OFingAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(14, 8))\n", "ax.errorbar(x, y, yerr=sy, xerr=0.5, label='Distribution of nSuccesses', fmt='.k', ecolor='k', elinewidth=1, capsize=1, capthick=1)\n", "ax.set(xlim=(xmin, xmax), ylim=(0, 1.2*np.max(y)), xlabel='Number of successes', ylabel='Frequency');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting a Binomial:\n", "\n", "Define the binomial function and plot it." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def func_binomial(x, N, n, p):\n", " return N * binom.pmf(x, n, p)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAHgCAYAAABn17aGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3zcZZn///c1M2k7PQAtLVhOFnYBgR5iCYUEpC1BREBAFOUo4LIoK7DfVbGgbjNFVFhWQB66uMoiiIgIKuLK4g8icjBRSKGcSpGDRUorPSCltOkhM9fvj8xMUpq2M+Vzmszr+XjkkWRmMvedTJvMe+7rvm5zdwEAAABAvUnFPQEAAAAAiANhCAAAAEBdIgwBAAAAqEuEIQAAAAB1iTAEAAAAoC4RhgAAAADUpUzcE3g3xo4d6xMmTIh7GgAAAAASau7cucvdfdxA19V0GJowYYK6urringYAAACAhDKzVzZ3HWVyAAAAAOoSYQgAAABAXSIMAQAAAKhLNb1nCAAAANHasGGDFi1apLVr18Y9FWAjw4YN02677aaGhoaKv4YwBAAAgIotWrRIo0aN0oQJE2RmcU8HkCS5u1asWKFFixZpzz33rPjrKJMDAABAxdauXasdd9yRIIREMTPtuOOOVa9YEoYAAABQFYIQkmhb/l0ShgAAAFBT0um0GhsbNWXKFE2dOlUdHR2SpMWLF+vjH/94qGN3dXXpoosu2uJtfv/73+u4447b6n3NmDGjfGbmMcccozfffHOzt7322mu1Zs2azV5/7rnnav78+ZKkkSNHbnXs/ubNm6d77rmn/Pndd9+tK664oqr7qFXsGQIAAEBNyWazmjdvniTpt7/9rS699FI9+OCD2mWXXXTnnXeGOnZTU5OampoCv9/+YWQg1157rc444wwNHz58k+vy+bxuuOGGbR573rx56urq0jHHHCNJOv7443X88cdv8/3VElaGAAAAULPeeustjR49WpK0cOFCTZw4UZJ000036aSTTtLRRx+tvffeW1/60pfKX3Pbbbdp0qRJmjhxombNmlW+fOTIkZo1a5YOPPBAHXnkkXr00Uc1Y8YM7bXXXrr77rslbbzq8+ijj6qlpUXvf//71dLSoueff36Lc+3u7tYpp5yiyZMn65Of/KS6u7vL102YMEHLly/X6tWrdeyxx2rKlCmaOHGibr/9dl133XVavHixZs6cqZkzZ5bnOnv2bB188MHq7OzcaJVJkr7whS9o6tSpam1t1bJlyyRtvBK1fPlyTZgwQevXr9fs2bN1++23q7GxUbfffrtuuukmXXDBBZKkV155Ra2trZo8ebJaW1v117/+VZJ09tln66KLLlJLS4v22muv0ENoWFgZAgAAwDaZ8+tnNX/xW4He5/67bKe2jxywxdt0d3ersbFRa9eu1ZIlS/S73/1uwNvNmzdPTzzxhIYOHap9991XF154odLptGbNmqW5c+dq9OjROuqoo3TXXXfpxBNP1OrVqzVjxgxdeeWV+uhHP6qvfvWruu+++zR//nydddZZm6yWvO9979NDDz2kTCaj+++/X1/+8pf185//fLPzvv766zV8+HA99dRTeuqppzR16tRNbnPvvfdql1120W9+8xtJ0sqVK7X99tvr6quv1gMPPKCxY8dKklavXq2JEyfqsssu2+Q+Vq9eralTp+pb3/qWLrvsMs2ZM0ff+c53BpzTkCFDdNlll6mrq6t8m5tuuql8/QUXXKBPfepTOuuss3TjjTfqoosu0l133SVJWrJkiR555BEtWLBAxx9/fOglimFgZQgAAAA1pVQmt2DBAt1777361Kc+JXff5Hatra3afvvtNWzYMO2///565ZVX9Nhjj2nGjBkaN26cMpmMTj/9dD300EOSeoPB0UcfLUmaNGmSpk+froaGBk2aNEkLFy7c5P5Xrlypk08+WRMnTtS//du/6dlnn93ivB966CGdccYZkqTJkydr8uTJm9xm0qRJuv/++zVr1iw9/PDD2n777Qe8r3Q6rY997GMDXpdKpfTJT35SknTGGWfokUce2eK8tqSzs1OnnXaaJOnMM8/c6L5OPPFEpVIp7b///nr99de3eYw4sTIEAACAbbK1FZwoNDc3a/ny5eVSsP6GDh1a/jidTqunp2fA0FTS0NBQ7kiWSqXKX59KpdTT07PJ7f/93/9dM2fO1C9/+UstXLhQM2bM2Op8t9bxbJ999tHcuXN1zz336NJLL9VRRx2l2bNnb3K7YcOGKZ1Ob3W8/mNmMhkVCgVJ2uZDc/vPv//Pd0s/1yRjZQgAAAA1a8GCBcrn89pxxx0ruv3BBx+sBx98UMuXL1c+n9dtt92m6dOnb9PYK1eu1K677ipp49KyzTn88MN16623SpKeeeYZPfXUU5vcZvHixRo+fLjOOOMMffGLX9Tjjz8uSRo1apRWrVpV0bwKhUJ5D89PfvITHXbYYZJ69yXNnTtXkjba47Ol+25padFPf/pTSdKtt95avq/BgpUhAAAA1JTSniGpd0Xi5ptvrniVZPz48frmN7+pmTNnyt11zDHH6IQTTtimeXzpS1/SWWedpauvvlpHHHHEVm9//vnn65xzztHkyZPV2NioadOmbXKbp59+WhdffLFSqZQaGhp0/fXXS5LOO+88ffjDH9b48eP1wAMPbHGcESNG6Nlnn9WBBx6o7bffXrfffrsk6Ytf/KI+8YlP6JZbbtlovjNnztQVV1yhxsZGXXrppRvd13XXXadPf/rTuuqqqzRu3Dj98Ic/3Or3WUusVpe0JKmpqcn7d80AAABAuJ577jntt99+cU8DGNBA/z7NbK67D9gPnTI5AAAAAHWJMAQAAACgLhGGAAAAANSl0MKQme1uZg+Y2XNm9qyZ/Wvx8pyZvWZm84pvx/T7mkvN7EUze97MPhTW3AAAAAAgzG5yPZK+4O6Pm9koSXPN7L7idde4+3/2v7GZ7S/pFEkHSNpF0v1mto+750OcIwAAAIA6FdrKkLsvcffHix+vkvScpF238CUnSPqpu69z979IelHSpv0GAQAAUFNyuZzMTLlcLu6pABuJZM+QmU2Q9H5JfypedIGZPWVmN5rZ6OJlu0p6td+XLdKWwxMAAABqQCkEBRWG0um0GhsbdcABB2jKlCm6+uqrVSgUJEldXV266KKLNvu1Cxcu1E9+8pPNXr948WJ9/OMfl9R7kOoFF1xQ1dxuuukmLV68uPz5ueeeq/nz51d1H9VYt26djjzySDU2NpbPE9qaNWvW6PTTT9ekSZM0ceJEHXbYYXr77bdDm2OShX7oqpmNlPRzSf/P3d8ys+slfU2SF99/S9KnJdkAX77JIUhmdp6k8yRpjz32CGvaAAAASKhsNqt58+ZJkpYuXarTTjtNK1eu1Jw5c9TU1KSmpgGPlJHUF4ZOO+20Ta7r6enRLrvsojvvvHOb53bTTTdp4sSJ2mWXXSRJN9xwwzbfVyWeeOIJbdiwofzzqMS3v/1t7bzzznr66aclSc8//7waGhrCmmKihboyZGYN6g1Ct7r7LyTJ3V9397y7FyT9QH2lcIsk7d7vy3eTtFjv4O7fd/cmd28aN25cmNMHAABAADo7Ozd6H6SddtpJ3//+9/Wd73xH7q7f//73Ou644yRJDz74oBobG9XY2Kj3v//9WrVqlS655BI9/PDDamxs1DXXXKObbrpJJ598sj7ykY/oqKOO0sKFCzVx4sTy/b/66qs6+uijte+++2rOnDmStMlt/vM//1O5XE533nmnurq6dPrpp6uxsVHd3d2aMWOGurq6JEm33XZbeTVm1qxZ5a8fOXKkvvKVr2jKlCk65JBD9Prrr2/yfb7xxhs68cQTNXnyZB1yyCF66qmntHTpUp1xxhmaN2+eGhsb9dJLL230NTNmzNCsWbM0bdo07bPPPnr44YclSUuWLNGuu/YVYO27774aOnToZr8vSXrxxRd15JFHasqUKZo6dWp5rP/4j//QpEmTNGXKFF1yySWSpJdeeklHH320DjzwQH3gAx/QggULJEl33HGHJk6cqClTpujwww+XJD377LOaNm2aGhsbNXnyZL3wwguSpB//+Mflyz/zmc8on88rn8/r7LPP1sSJEzVp0iRdc801lf0j2RJ3D+VNvSs9P5J07TsuH9/v439T7z4hqbdxwpOShkraU9LLktJbGuPAAw90AAAARGf+/PlbvU1bW1v5446ODs9msy7Js9msd3R0DHi7aowYMWKTy3bYYQf/29/+5g888IAfe+yx7u5+3HHH+SOPPOLu7qtWrfINGzZsdL27+w9/+EPfddddfcWKFe7u/pe//MUPOOCA8nXvec97fPny5b5mzRo/4IAD/LHHHtvoNu7uV111Vfl7mT59uj/22GPl60qfv/baa7777rv70qVLfcOGDT5z5kz/5S9/6e7ukvzuu+92d/eLL77Yv/a1r23y/V1wwQWey+Xc3b29vd2nTJni7r7J99Pf9OnT/fOf/7y7u//mN7/x1tZWd3d/4oknfNy4cX7IIYf4V77yFf/zn/+8yff+zu9r2rRp/otf/MLd3bu7u3316tV+zz33eHNzs69evdrdvfwzPOKII8r3+cc//tFnzpzp7u4TJ070RYsWubv73//+9/L39eMf/9jd3detW+dr1qzx+fPn+3HHHefr1693d/fzzz/fb775Zu/q6vIjjzyyPL/SffQ30L9PSV2+mTwR5srQoZLOlHTEO9po/4eZPW1mT0maWQxEcvdnJf1M0nxJ90r6nNNJDgAAoObMmTNHZiYzU0tLi7q7uyVJ3d3damlpKV9XWmkJQu9z3o0deuih+vznP6/rrrtOb775pjKZgXeIfPCDH9SYMWM2e92OO+6obDark046SY888sg2ze+xxx7TjBkzNG7cOGUyGZ1++ul66KGHJElDhgwpr2YdeOCBWrhw4SZf/8gjj+jMM8+UJB1xxBFasWKFVq5cudVxTzrppE3ut7GxUS+//LIuvvhivfHGGzrooIP03HPPbfY+Vq1apddee00f/ehHJUnDhg3T8OHDdf/99+ucc87R8OHDJUljxozR22+/rY6ODp188snlVZ0lS5ZI6n08zj77bP3gBz9QPt/7NL+5uVnf+MY3dOWVV+qVV15RNptVe3u75s6dq4MOOkiNjY1qb2/Xyy+/rL322ksvv/yyLrzwQt17773abrvttvr9b02Y3eQecXdz98nu3lh8u8fdz3T3ScXLj3f3Jf2+5uvu/g/uvq+7/19YcwMAAEB42trayq+8d3R0KJvNSurd69PR0VG+rq2tLZDxXn75ZaXTae20004bXX7JJZfohhtuUHd3tw455JByudY7jRgxYrP3bWabfJ7JZMoNGyRp7dq1W53jQGGtpKGhoTxOOp1WT09PRV//zrkNZOjQoQPe78iRI3XSSSfpv/7rv3TGGWfonnvu2ez3tbm5u/smcygUCtphhx00b9688lspaH3ve9/T5ZdfrldffVWNjY1asWKFTjvtNN19993KZrP60Ic+pN/97ndyd5111lnlr3/++eeVy+U0evRoPfnkk5oxY4a++93v6txzz93q9781kXSTAwAAQP3o3zWuublZ7e3tkqT29nY1NzcPeLtttWzZMn32s5/VBRdcsMkT85deekmTJk3SrFmz1NTUpAULFmjUqFFatWpVxfd/33336Y033lB3d7fuuusuHXroodp55521dOlSrVixQuvWrdP//u//lm+/ufs/+OCD9eCDD2r58uXK5/O67bbbNH369Irncfjhh+vWW2+VJP3+97/X2LFjt3ll5A9/+IP+/ve/S5LWr1+v+fPn673vfe9mv6/ttttOu+22m+666y5JvR3s1qxZo6OOOko33nij1qxZI6l3X9N2222nPffcU3fccYek3sD05JNPSup9PA4++GBddtllGjt2rF599dXyis9FF12k448/Xk899ZRaW1t15513aunSpeX7feWVV7R8+XIVCgV97GMf09e+9jU9/vjj2/T99xd6NzkAAADUt1IA6h+E3o3u7m41NjZqw4YNymQyOvPMM/X5z39+k9tde+21euCBB5ROp7X//vvrwx/+sFKplDKZjKZMmaKzzz5bo0ePHmCEPocddpjOPPNMvfjiizrttNPKnepmz56tgw8+WHvuuafe9773lW9/9tln67Of/ayy2exGDSPGjx+vb37zm5o5c6bcXcccc4xOOOGEir/nXC6nc845R5MnT9bw4cN18803V/y17/TSSy/p/PPPl7urUCjo2GOP1cc+9jGZ2Wa/r1tuuUWf+cxnNHv2bDU0NOiOO+7Q0UcfrXnz5qmpqUlDhgzRMccco2984xu69dZbdf755+vyyy/Xhg0bdMopp2jKlCm6+OKL9cILL8jd1draqilTpuiKK67Qj3/8YzU0NOg973mPZs+erTFjxujyyy/XUUcdpUKhoIaGBn33u99VNpvVOeecU169+uY3v7nNP4MS29KSXdI1NTV5qTsHAAAAwvfcc89pv/32q/rrzGyLpWJAEAb692lmc919wH7rlMkBAAAgVEEfugoEhTI5AAAAhCqXyxGEkEisDAEAAACoS4QhAAAAVIW9P0iibfl3SRgCAABAxYYNG6YVK1YQiJAo7q4VK1Zo2LBhVX0de4YAAABQsd12202LFi3SsmXL4p4KsJFhw4Zpt912q+prCEMAAACoWENDg/bcc8+4pwEEgjI5AAAAAHWJMAQAAACgLhGGAAAAANQlwhAAAACAukQYAgAAAFCXCEMAAAAA6hJhCAAAAEBdIgwBAAAAqEuEIQCJlsvlZGbK5XJxTwUAAAwy5u5xz2GbNTU1eVdXV9zTABAyM1Mt/64CAADxMbO57t400HWsDAEAAACoS4QhAAAAAHWJMAQAAACgLhGGAAAAANQlwhAAAACAukQYAgAAAFCXCEMAAAAA6hJhCAAAAEBdIgwBNS6Xy8nMlMvl4p4KAABATbFaPtW9qanJu7q64p4GEDszUy3/X96awf79AQCA8JjZXHdvGug6VoYAAAAA1CXCEAAAAIC6RBgCAAAAUJcIQwAAAADqEmEIAAAAQF0iDAEAAACoS4QhAAAAAHWJMAQAAACgLhGGAAAAANQlwhAAAACAukQYApBonZ2dG70HAAAICmEIQNVyuVwk43R2dqq1tVWS1NraGlkgiur7AwAA8TJ3j3sO26ypqcm7urringYQOzNTlP+XzSyyseJSy78bAQBAHzOb6+5NA13HyhCAqrW1tcndQ3/r6OhQNpuVJGWzWXV0dEQybltbW8w/YQAAEAVWhoBBIOqVoSh1dnaqpaVFHR0dam5ujns6AACgxrAyBKBmlQIQQQgAAAQtE/cEANSeP728Qlff92dFtRi182lX6BPfi6Z5wqTdtte/H7d/JGMBAIB4EYYAVK19wVJ1vfJ3TZswJpoBC3mlU+E3bXhlxWo9s3glYQgAgDpBGAJQtVVrN2jMiCG67bxDIhnPPtOs224Lfxnqiv9boBv/8JfQxwEAAMnAniEAVVu1tkejhg6+11IyKVO+MDgbUQAAgE0RhoAaVzqINKoDSSXp7XU9GjVs8IWhdDEMDdbOfAAAYGOEISAkuVwu9DE6OzvV2toqSWptbY0sEM3/88saOQjDUKa4L6mH1SEAAOoCYQgIyZw5c2Rmob61tLSou7tbktTd3a2WlpbQxzQzLVq6QiMHYZlcOt0bhiiVAwCgPhCGgJC0tbXJ3UN96+joUDablSRls1l1dHSEPqa7a/S48Ro5tCHmn3DwGlK9vxJZGQIAoD4QhoCQRFEm19zcrPb2dklSe3t7ZAeTpoeNGLR7hiQpnycMAQBQDwhDQI0rBaCogpC7D9oGCpl0ac9QIeaZAACAKBCGAFRlzfq8Cq7BuWeIBgoAANQVwhCAqry9rkeS6CYHAABqHmEIQFVWre0NQ6OGDb4GCpliAwX2DAEAUB8IQwCqsmrtBknSqEFYJseeIQAA6gthCEBVBnOZXLmbHGVyAADUBcIQgKq8XS6TG3xhiD1DAADUF8IQgKqsKq0MDcIyuXTp0FX2DAEAUBcIQwCqUm6gMHQwNlBgzxAAAPWEMASgKqUyuRFD0zHPJHilBgrsGQIAoD4QhgBU5e11GzR8SFqZ9OD79cGhqwAA1JfB92wGQKjeXtczKPcLSf3OGSIMAQBQFwbnMxogZj/qXKjfLVga2Xg7fTyns3/4aCRjzV/8VqRttXO5XPl96eOwlFaGNuSj2zOUy+U0Z84ctbW1hf79AQCAjZl77b4C2tTU5F1dXXFPA9jEh655SEtWdmvPsSMiGe/RRx/TtGkHRTKWJB253866sHXvyMaLypOvvqkTvvsH/c9ZTWrdb+fIxjUz1fLvYgAAkszM5rp700DXsTIEhCDvrsP2Hqv/Ov3ASMYz+4B+9SOeTL9bpQYK7BkCAKA+sGcICEHBXSmzuKeBKrFnCACA+kIYAkJQKBCGahHd5AAAqC+EISAEeffyE2vUjvKhqxE2UAAAAPEJLQyZ2e5m9oCZPWdmz5rZvxYvH2Nm95nZC8X3o4uXm5ldZ2YvmtlTZjY1rLkBYSsUxMpQDWJlCACA+hLmylCPpC+4+36SDpH0OTPbX9IlktrdfW9J7cXPJenDkvYuvp0n6foQ5waEqnfPUNyzQLUa0uwZAgCgnoQWhtx9ibs/Xvx4laTnJO0q6QRJNxdvdrOkE4sfnyDpR97rj5J2MLPxYc0PCFO+QJlcLWJlCACA+hLJniEzmyDp/ZL+JGlnd18i9QYmSTsVb7arpFf7fdmi4mVAzSm4lIooDPU/lBTvTmnPUJ49QwAA1IXQw5CZjZT0c0n/z93f2tJNB7hsk5dnzew8M+sys65ly5YFNU0gUFGWyeVyObk7YSgAac4ZAgCgroQahsysQb1B6FZ3/0Xx4tdL5W/F90uLly+StHu/L99N0uJ33qe7f9/dm9y9ady4ceFNHngX8gVXmgYKNSdDmRwAAHUlzG5yJul/JD3n7lf3u+puSWcVPz5L0q/6Xf6pYle5QyStLJXTAbWm4B5ZmRyCw6GrAADUl0yI932opDMlPW1m84qXfVnSFZJ+Zmb/JOmvkk4uXnePpGMkvShpjaRzQpwbECoOXa1NfecMEYYAAKgHoYUhd39EA+8DkqTWAW7vkj4X1nyAKHHoam1KpUxmUr5AAwUAAOpBJN3kgHrDoau1K5MybaBMDgCAukAYAkLAoau1K50y9gwBAFAnCENACCiTq10NqRR7hgAAqBOEISBg7i53yuRqVTpt7BkCAKBOEIaAgJUqrAhDtSmTMs4ZAgCgThCGgICV9puk+d9Vk9Ipo0wOAIA6wdM1IGAF730izaGrtSmTSkW6MtTZ2bnRewAAEB3CEBCwUhhKUyZXkzJp0xPznoxkrM7OTrW29h671traGlkgyuVykYwDAEDShXboKlCvSmVy7BmqTemU6elnn5VF/Ph1d3erpaUlsvEIRAAAsDIEBK7UiIwyudqUSZn2P2BisStguG8dHR3KZrOSpGw2q46OjkjGbWtri/mnDABAMhCGgID1lcnFPBFsk3QqpX3et18kYzU3N6u9vV2S1N7erubm5kjGZVUIAIBehCEgYHkaKNS0TMrUk4/unKFSAIoqCAEAgD6EISBgBfYM1bQ05wwBAFA3CENAwErPo9OsDNWkhrSVm2AAAIDBjTAEBKxcJkcWqkmsDAEAUD8IQ0DAKJOrbZlUipUhAADqBGEICFi5mxxLQzUpHXEDBQAAEB/CEBAwDl2tbRnK5AAAqBuEISBgBVpr17QMDRQAAKgbhCEgYOVucqwM1aRMKsXKEAAAdYIwBASsr0wu5olgm6RTrAwBAFAvCENAwMphiDRUkzIp0wYaKAAAUBcIQ0DAnDK5msbKEAAA9YMwBASsfOgq/7tqUibNniEAAOoFT9eAgNFau7ZlWBkCAKBuEIaAgDmHrtY0Dl0FAKB+EIaAgLEyVNs4dBUAgPpBGAICVt4zRBiqSek0YQgAgHpBGAICVu4mR5lcTWpIpdgzBABAnSAMAQHj0NXaVmqtXdr7BQAABi/CEBCwvtbapKFalCk+bqwOAQAw+BGGgICVu8mxZ6gmpdO9j1tU+4ZyudxG7wEAQHQycU8AGGxKXZlpoFCbSitDUYYhghAAAPFgZQgIWHnPEP+7alKm+MDl85TJAQAw2PF0DQgYh67Wtky5TI6DVwEAGOwokwMCxjlDta0UYp9Z/JZ2HDEk5tkEb8TQjPYcOyLuaQAAkAiEISBgfa21CUO1aOTQ3l+LZ934aMwzCc9vLjpMB+yyfdzTAAAgdoQhIGAculrbPjxxvHb49BCt7xl8ZXJ/fn2Vrvrt81q5ZkPcUwEAIBEIQ0DAOHS1tg3JpDR9n3FxTyMU22cbJEkcoQQAQC8aKAABY88Qkipd/I1f+jcKAEC9IwwBAaObHJKqFNALLA0BACCJMAQEjkNXkVSlgJ4nDAEAIIkwBASuXCbH/y4kTCmgUyYHAEAvnq4BASuXybEyhIQphSEnDAEAIIkwBASOc4aQVH1lcjFPBACAhCAMAQErhyEaKCBh6CYHAMDGCENAwDh0FUlFNzkAADZGGAIC1nfOUMwTAd6BbnIAAGyMMAQEjD1DSKryyhBlcgAASCIMAYHj0FUkVWkfG2EIAIBehCEgYBy6iqQqtXunmxwAAL0IQ0DA2DOEpErRTQ4AgI0QhoCAubtSJhkrQ0iYNN3kAADYCGEICFi+4JTIIZHS7BkCAGAjhCEgYHl3DlxFIpnRWhsAgP4IQ0DA3PvKkYAkYWUIAICNEYaAgPWWycU9C2BTdJMDAGBjhCEgYPkCZXJIplI3OVaGAADoRRgCAubuHLiKRKKbHAAAGyMMAQHLO93kkEylf5ecMwQAQC/CEBCwfEGEISRSqXyTlSEAAHoRhoCAFQquNP+zkFDplLEyBABAEU/ZgIAVKJNDgqXN6CYHAEARYQgIGHuGkGSpFN3kAAAoIQwBAestkyMMIZlSZuwZAgCgiDAEBKzgIgwhsdLGniEAAEoIQ0DA8u6iSg5JlUqxMgQAQAlhCAhYoeDlwy2BpKGbHAAAfQhDQMAKzp4hJFeKbnIAAJQRhoCA5cEZmt0AACAASURBVAuSsTKEhEqZ5KwMAQAgiTAEBK53ZSjuWQADS6dMefYMAQAgiTAEBK7g7BlCcqXoJgcAQBlhCAhYvuCUySGx0nSTAwCgLLQwZGY3mtlSM3um32U5M3vNzOYV347pd92lZvaimT1vZh8Ka15A2GiggCTr7SYX9ywAAEiGMFeGbpJ09ACXX+PujcW3eyTJzPaXdIqkA4pf819mlg5xbkBoCgVRJofEMusN7AAAIMQw5O4PSXqjwpufIOmn7r7O3f8i6UVJ08KaGxAmDl1FkqWNMjkAAEri2DN0gZk9VSyjG128bFdJr/a7zaLiZUDNKRQok0Ny0U0OAIA+UYeh6yX9g6RGSUskfat4+UDPHAf8a21m55lZl5l1LVu2LJxZAu8Ce4aQZCkzyuQAACiKNAy5++vunnf3gqQfqK8UbpGk3fvddDdJizdzH9939yZ3bxo3bly4Ewa2Qd45dBXJxcoQAAB9Ig1DZja+36cflVTqNHe3pFPMbKiZ7Slpb0mPRjk3ICiFgitNFkJCpUwiCwEA0CsT1h2b2W2SZkgaa2aLJLVJmmFmjeotgVso6TOS5O7PmtnPJM2X1CPpc+6eD2tuQJgok0OSpVKUyQEAUFJRGDKzie7+zNZv2cfdTx3g4v/Zwu2/Lunr1YwBJBGHriLJ0kaZHAAAJZWWyX3PzB41s38xsx1CnRFQ4wrunDOExEqxZwgAgLKKwpC7HybpdPU2Oegys5+Y2QdDnRlQowouyuSQWGm6yQEAUFZxAwV3f0HSVyXNkjRd0nVmtsDMTgprckBQcrmczEy5XC70sQoFDl1FcqVSNFAAAKCkojBkZpPN7BpJz0k6QtJH3H2/4sfXhDg/IBClEBRFGMrTQAEJlmLPEAAAZZV2k/uOes8F+rK7d5cudPfFZvbVUGYG1Cj2DCHJ0nSTAwCgrNIwdIyk7lK7azNLSRrm7mvc/ZbQZgfUoEKBQ1eRXHSTAwCgT6V7hu6XlO33+fDiZQDeIV9wpSM9zhioHN3kAADoU+lTtmHu/nbpk+LHw8OZElDbOHQVSZYyiSo5AAB6VRqGVpvZ1NInZnagpO4t3B6oWwXn0FUkVzplypOGAACQVPmeof8n6Q4zW1z8fLykT4YzJaC25Qs0UEBypcxUoEwOAABJFYYhd3/MzN4naV9JJmmBu28IdWZAjeLQVSQZK0MAAPSpdGVIkg6SNKH4Ne83M7n7j0KZFVDDOHQVSUY3OQAA+lQUhszsFkn/IGmepHzxYpdEGALeIc85Q0iwVMpooAAAQFGlK0NNkvZ3508osDV0k0OSpUysDAEAUFRpN7lnJL0nzIkAgwWHriLJ2DMEAECfSleGxkqab2aPSlpXutDdjw9lVqgLDzy/VH/+26rIxttu2kn67wdfCn2cDYUCh64isegmBwBAn0rDUC7MSaA+/ettT+ittT2RjTd65qf1zf9bEMlY7x0zIpJxgGqxMgQAQJ9KW2s/aGbvlbS3u99vZsMlpcOdGga7DXnXOYdO0MUf2jeS8UaOHKm333479HFSZhrWwH8PJBMrQwAA9Km0m9w/SzpP0hj1dpXbVdL3JLWGNzUMdnl3Dc2kNXxINR3et51vWBfZWEBSpcxEFgIAoFelOxs+J+lQSW9Jkru/IGmnsCaF+uDuiqrpWmdn50bvgXqVTtFNDgCAkkrD0Dp3X1/6xMwy6j1nCNhmPfmCUhF0Xevs7FRra+8iZmtra2SBKJfLRTIOUI0Ue4YAACirtGboQTP7sqSsmX1Q0r9I+nV400I9cJkuv/wyXXz0TyIbs7u7Wy0tLZGNRyBC0qTZMwQAQFmlK0OXSFom6WlJn5F0j6SvhjUpDH6l83tzbW1y91DfOjo6lM1mJUnZbFYdHR2hj+nuamtri/NHDAyIbnIAAPSptJtcQdIPim/Au1Z6YTqKMrnm5ma1t7erpaVF7e3tam5uDn1MiVUhJJOZyb33BQkOBwYA1LtKu8n9RQPsEXL3vQKfEepCaQN3OqIOCqUAFFUQApIqXQxABZfSZCEAQJ2rdM9QU7+Ph0k6Wb1ttoFtUiiW6fDCNBCtdLE4Ol/wyF6MAAAgqSraM+TuK/q9vebu10o6IuS5YRDzCMvkAPRJpUorQ+wbAgCg0jK5qf0+Tal3pWhUKDNCXSht4OaFaSBapTI5zhoCAKDyMrlv9fu4R9JCSZ8IfDaoG4VyGCINAVFKGStDAACUVNpNbmbYE0F98ULve8IQEK1ymVwh5okAAJAAlZbJfX5L17v71cFMB/WiQJkcEItSBznOGgIAoLpucgdJurv4+UckPSTp1TAmhcGv9ESMblZAtEr/59gzBABA5WForKSp7r5KkswsJ+kOdz83rIlhcOtrrU0YAqJENzkAAPpU1Fpb0h6S1vf7fL2kCYHPBnWD1tpAPGigAABAn0pXhm6R9KiZ/VKSS/qopB+FNisMeqUSnXSlcRxAIGitDQBAn0q7yX3dzP5P0geKF53j7k+ENy0MdpTJAfGgmxwAAH2qeV1+uKS33P3bkhaZ2Z4hzQl1gDI5IB6l1Vi6yQEAUGEYMrM2SbMkXVq8qEHSj8OaFAY/yuSAeKQokwMAoKzSp6IflXS8pNWS5O6LJY0Ka1IY/PrOGWJlCIhS6f+cszIEAEDFYWi99/7ldEkysxHhTQn1oPSiNHuGgGiVzxkiDAEAUHEY+pmZ/bekHczsnyXdL+kH4U0Lg52XV4aiGS+Xy230HqhXlMkBANCn0m5y/2lmH5T0lqR9Jc129/tCnRkGtdKr0umIVoZyuRxBCFDfyhDd5AAAqGBlyMzSZna/u9/n7he7+xcJQni3Sk/EKJMDohV1N7lcLicz48UIAEAibTUMuXte0hoz2z6C+aBOFCIukwPQq/QCRCHCMNT/PQAASVJRmZyktZKeNrP7VOwoJ0nuflEos8KgV3oiliYNAZEqlaYW2DMEAEDFYeg3xTcgEAUOXQViUe4mRxgCAGDLYcjM9nD3v7r7zVFNCPWhtDJEFgKiVe4mR2ttAAC2umfortIHZvbzkOeCOlIq0aFMDogW3eQAAOiztTDU/5nqXmFOBPWFMjkgHqXXH6JqoAAAQJJtLQz5Zj4G3hXK5IB4pFKUyQEAULK1BgpTzOwt9a4QZYsfq/i5u/t2oc4Og1apTI6VISBadJMDAKDPFsOQu6ejmgjqS+l5GHuGgGjRTQ4AgD5bPXQVCAOHrgLxSEV86CoAAElGGEIs+vYMkYaAKJW7yZGFAAAgDCEepTCUJgwBkSqtxlImBwAAYQgxKZ1xQgMFIFqpFGVyAACUEIYQC1prA/EorcayMgQAAGEIMSmXydFBAYhU1N3kOjs7N3oPAECSEIYQi9LzMMrkgGiVyuSiqJLr7OxUa2urJKm1tTWyQJTL5SIZBwBQ+7Z26CoQClprA/Eo/Z/7p3/+Z33yyd9GNm53d7daWloiG49ABACoBCtDiEWpRCdFGgIiVdozdP33/lvuHupbR0eHstmsJCmbzaqjoyP0Md1dbW1tcf6IAQA1hJUhxMIpkwNiUdozdPV9f9YND78c+nj7fvo/9NKjv9M/TDtClz6yVnrkgdDGGpJJ6ZpPNrIqBACoGGEIsaBMDojHmBFD9Jnpe+lvK9dGM+Duh2rRK6/o8MMODXWYNevzum/+63rmtZU6YJftQx0LADB4EIYQi3KZHCtDQKTMTJd+eL9Ix7zu1Kn69q+vCnWMv61cq/vmv658IdRhAACDDHuGEItymRxLQwACkCr+NcsXSEMAgMoRhhALyuQABClTTEMcJgsAqAZhCLHgnCEAQSp1ycuThQAAVSAMIRZ5Z88QgOBQJgcA2BaEIcTCKZMDEKC+MrmYJwIAqCmEIcSiQDc5AAEqrQyV9iMCAFAJwhBikaebHIAAlVaGetg0BACoAmEIsaBMDkCQSr9L8qwMAQCqQBhCLAo0UAAQIDNTyvpKcAEAqERoYcjMbjSzpWb2TL/LxpjZfWb2QvH96OLlZmbXmdmLZvaUmU0Na15IhtIm5zRLQwACkkml1EMYAgBUIcyVoZskHf2Oyy6R1O7ue0tqL34uSR+WtHfx7TxJ14c4LyRAaWWIhSEAQUmlaKAAAKhOaGHI3R+S9MY7Lj5B0s3Fj2+WdGK/y3/kvf4oaQczGx/W3BA/p0wOQMDSZsqzMgQAqELUe4Z2dvclklR8v1Px8l0lvdrvdouKl2GQKpXJEYaAwS2Xy230PkzpFGEIAFCdpDRQGOgZ8YB/0czsPDPrMrOuZcuWhTwthKVANzmgLuRyObk7YQgAkEhRh6HXS+VvxfdLi5cvkrR7v9vtJmnxQHfg7t939yZ3bxo3blyok0V43F1mvR2gACAI6ZTRWhsAUJWow9Ddks4qfnyWpF/1u/xTxa5yh0haWSqnw+BUcErkAAQrnTLlOXQVAFCFTFh3bGa3SZohaayZLZLUJukKST8zs3+S9FdJJxdvfo+kYyS9KGmNpHPCmheSIe+uNGEIQIDSxsoQAKA6oYUhdz91M1e1DnBbl/S5sOaC5CkUy+QAICjptHHoKgCgKklpoIA645TJAQhY2oxDVwEAVSEMIRb5gitNKzkAAUrRQAEAUCXCEGJBmRyAoGVooAAAqBJhCLGgTA5A0FI0UAAAVIkwhFhQJgcgaOkUDRQAANUhDCEWBXeRhQAEKZOigQIAoDqEIcSi4JJRJgcgQKmUqUCZHACgCoQhxMJZGQIQsLSZ8qwMAQCqQBhCLPIFV5qVIQABSlMmBwCoEmEIsaBMDkDQaKAAAKgWYQixcHel+NcHIEBpDl0FAFSJp6OIRd4pkwMQrHSKPUMAgOoQhhCLAoeuAggYDRQAANUiDCEWBXeRhQAEiZUhAEC1CEOIRaHgStNbG0CACEMAgGoRhhCLgjtlcgAClaKBAgCgSoQhxILW2gCClqG1NgCgSoQhxKJQcFElByBIaePQVQBAdQhDiEXB2TMEIFgpVoYAAFUiDCEWlMkBCFqGPUMAgCoRhhCL3gYKcc8CwGCSopscAKBKhCHEouCuNCtDAALEoasAgGoRhhCLQkG01gYQKM4ZAgBUizCEWBTcRRYCECTCEACgWoQhxIJucgCClqaBAgCgSoQhxKLglMkBCBYrQwCAahGGEAvK5AAEjQYKAIBqEYYQi0KBMjkAwUqnTAWXnFI5AECFCEOIBWVyAIJWeoGF1SEAQKUIQ4gFh64CCFo5DLEyBACoEGEIscgXXMbKEIAAlcJQoRDzRAAANYMwhFi49252BoCglH6n9JCGAAAVIgwhFgV3pfjXByBAKVaGAABV4ukoYtHbWpuVIQDBybBnCABQJcIQYlGgTA5AwEorQ5TJAQAqRRhCLOgmByBopRdYyEIAgEoRhhCL3jBEGgIQHMrkAADVIgwhFoVCX0kLAASh9DslnycMAQAqQxhCLCiTAxA0VoYAANUiDCEWlMkBCFp5ZahAGAIAVIYwhFjkKZMDELBSAwXCEACgUoQhxMIpkwMQsDQrQwCAKhGGEAvK5AAErRSGCuwZAgBUiDCEWBRchCEAgUoX/6L1sDIEAKgQYQixKBRYGQIQrHSq908aZXIAgEoRhhALWmsDCFqpgQJlcgCAShGGEIuC000OQLBSpTI5Dl0FAFSIMIRY5GmgACBgmWIaYmUIAFApwhBiQWttAEGjgQIAoFqEIcSCbnIAglb6nVIgDAEAKkQYQizyBWfPEIBAZegmBwCoEmEIkfNiPT9ZCECQUpTJAQCqRBhC5ErPUyiTAxAkGigAAKpFGELkSiUsrAwBCFKpgQJlcgCAShGGELnSq7bsGQIQpNJqM2EIAFApwhAi55TJAQgBDRQAANUiDCFyBRooAAhBqYFCnj1DAIAKEYYQuXw5DJGGAAQnnaJMDgBQHcIQIueF3veEIQBBIgwBAKpFGELkKJMDEIZ08QUWWmsDACpFGELkSmVyadIQgACVfqf05AlDAIDKEIYQudKrtkaZHIAAlcIQK0MAgEoRhhA5WmsDCAN7hgAA1crEPQEMbH1PQbf+6RWtXtcT91QCt2pt7/eUJooDCFDpBZYewhAAoEKEoYSa9+qbmvPr+XFPIzSZlGn30cPjngaAQSRTKpMjDAEAKkQYSqj1Pb39p2/750PUNGF0zLMJnknKsDQEIEDlMjn2DAEAKkQYSqieQm8YGtqQUgOhAQC2ysyUMvYMAQAqx7PshCr9Mc/QfhoAKpZOGWEIAFAxVoYSqrQBmLN4AKByKTM9/dpK3fqnV+KeSigO+8exeu+OI+KeBgAMGoShhOpbGWLxDgAqtcsOWT38wnI9/MLyuKcSimMnj9d3T5sa9zQAYNAgDCUUK0MAUL3/+9cP6K3uDXFPIxSfuvFRda/Pxz0NABhUCEMJlS82UGDPEABUblhDWsMa0nFPIxTZIWltyBfingYADCrUYCVUTz7alaFcLiczUy6Xi2Q8AEB1GtKp8rELAIBgxBKGzGyhmT1tZvPMrKt42Rgzu8/MXii+H3yH61ShvGcoHV0Y6v8eAJAsQ9Kpcgk1ACAYca4MzXT3RndvKn5+iaR2d99bUnvx87rFniEASL4oV9UzaaNMDgAClqQyuRMk3Vz8+GZJJ8Y4l9jRTQ4Aki/KVXXK5AAgeHE903ZJ/5+ZzTWz84qX7ezuSySp+H6nmOaWCKwMAQD6o0wOAIIXVze5Q919sZntJOk+M1tQ6RcWw9N5krTHHnuENb/Y0U0OANAfZXIAELxYVobcfXHx/VJJv5Q0TdLrZjZekorvl27ma7/v7k3u3jRu3Liophw5VoYAAP01pFPlTqMAgGBEHobMbISZjSp9LOkoSc9IulvSWcWbnSXpV1HPLUny+dKeIcIQAEBqSJvWszIEAIGKo0xuZ0m/NLPS+D9x93vN7DFJPzOzf5L0V0knxzC3xGBlCADQX0M6RZkcAAQs8jDk7i9LmjLA5SsktUY9n6TKF1zplKkYGkPX2dlZft/c3BzJmACAylEmBwDBo29zQvUUXF7IRzJWZ2enWlt7c2hra2s5GIWNA14BoHIZyuQAIHBxdZPDVuQLBfWsXxfZylBJd3e3WlpaIhuPQASglkW5qj6EMjkACBwrQwnVU3ANHdIgdw/9raOjQ9lsVpKUzWbV0dERybhtbW0x/5QBDEZRvcgS9ap6Qzold2l2Wy7UcQCgnph77dYfNzU1eVdXV9zTCMXsXz2jXz+5WE/MPiqS8To7O9XS0qKOjg72DAGoaVGvqEdlu4M/rtEzztZfv3WSChvWxT0dAKgZZjbX3ZsGuo6VoYTqKbjSqegenlIAIggBqHVtbW2DclX9W1ddKUm65MtfjfPHCwCDCmEoofJ554whANgGUZXJNTc3q729XZLU3t4e+otJDeneP9lf/NKsUMcBgHpCGEqonmJrbQBAckW5qp5J9/5NoIkCAASHMJRQ+UKh/IcPAIDSytD6HsIQAASFMJRQrAwBAPobUgxDPYXabXwEAElDGEqofIE9QwCAPpTJAUDwCEMJFXU3OQBAslEmBwDB49l2QrEyBADojzI5AAgeYSihot4zVGpFG1VLWgBAdSiTA4DgZeKeAAaWLxQiXRnK5XIEIQBIsFKZ3AbK5AAgMKwMJVRPnm5yAJB0Ua6ql8MQZXIAEBjCUELlC845QwCQcLlcTu4eURgqlsmxMgQAgSEMJRTd5AAA/ZVXhtgzBACB4dl2QtFNDgDQX3lliDI5AAgMYSihou4mBwBINhooAEDwCEMJFXU3OQBAsjWUzxkiDAFAUAhDCcXKEACgv1JTnfV5yuQAICiEoYRizxAAoL8hlMkBQOAIQwnVe84QDw8AoBdlcgAQPJ5tJxQrQwCA/kplchsokwOAwBCGEqqn4Epz6CoAoKihWC2wnjI5AAgMYSih6CYHAOgvlTJlUkaZHAAEiDCUUHSTAwC8UyZtlMkBQIAIQwmVL7jSRhgCAPRpSKcokwOAABGGEoo9QwCAdxqSTlEmBwABIgwlFN3kAADvlEmbNvRQJgcAQSEMJZC795bJcc4QAKCfhnRKG/KsDAFAUHi2nUCF4ot+rAwBAPobkk5pQ4GVIQAICmEogUr14HSTAwD011smx8oQAASFMJRA+eKrfqwMAQD6o0wOAIJFGEqgnmIYYmUIANBfA2VyABAowlAC5fOsDAEANtVAmRwABIowlEDllaE0Dw8AoA9lcgAQLJ5tJxB7hgAAA6FMDgCCRRhKILrJAQAGEnWZXC6Xk5kpl8tFNiYARCkT9wSwKVaGAAADaUintLJ7gx5YsDSS8aaf8i8adstvNP2Uf4lkzF1HZ7XPzqNCHwcASghDCUQ3OQDAQEaPGKLX3uzWOTc9FtmYO5+ci2y8EUPSembOh2TG3z8A0SAMJVDfyhBVjACAPrOP21+faNo90jEPPvhg/elPfwp9nF88vkg/6nxF63oKGtaQDn08AJAIQ4nUk2dlCACwqWENaTXuvkOkY65f8udIxnzy1TclSWvW5wlDACLD0kMCsWcIAFBvRgztfX129bqemGcCoJ4QhhKo3E0uTRgCAMSns7Nzo/dhGjGkdzVo9XrCEIDoEIYSiJUhAMDmRNXmurOzU62trZKk1tbW0APRcFaGAMSAMJRAdJMDAGzOnDlzZGahv7W0tKi7u1uS1N3drZaWllDH+9ARh0uSVq/Lx/njBVBnCEMJRDc5AMDmtLW1yd1Df+vo6FA2m5UkZbNZdXR0hDreE4/1dqxjZQhAlHi2nUCsDAEANieqMrnm5ma1t7dLktrb29Xc3BzqeCNLZXLrWRkCEB3CUALliw0U2DMEAIhTKQCFHYQkaXixgcIaGigAiBBhKIE4ZwgAUG9KrbXfpkwOQIQIQwlU3jNEa20AQJ0YmkkpnTKtoYECgAgRhhKoh9baAIA6Y2YaPiTNyhCASBGGEihfbqDAwwMAiE+pWUNUTRtGDMmwZwhApDJxTwCbYmUIAJAEuVwusiAkSSOGpjlnCECkWHpIoFI3ORooAADqyYihGa1mZQhAhAhDCcTKEACgHo0YkqGBAoBI1XSZ3PN/W6UZVz0QyVgr3nhDb6x4Q2N2HKMdx4wJday31va+KsbKEACgnowYmtbiN9fGPQ0AdaSmw9DwIWlN2X2HaAbbfQf95NY/6ogjTotkuPdsP0xjRgyJZCwAAJJgOA0UAESspsPQ7mOG69unvD+y8a47daq+/eurIhsPAIB6MmJoRm9TJgcgQuwZqlBnZ+dG7wEAQLBGDEmzMgQgUjUdhhYvXhzJOJ2dnWptbZUktba2RhaIomxnCgBA3IYPzWjN+rwKxUZCABC2mi6TW7JkicyibTLQ3d2tlpaWyMYjEAEA6sXIoWlJ0poNeY0cWtNPUQDUiJr+TTN+/PhIVodKK0Pd3d3KZrNqb29Xc3Nz6OMShAAA9WT4kN6nJX9buVbjRg0NdawrrrhCV155pWbNmqVLLrkk1LEkaWgmpWEN6dDHAVAdc6/dpeimpibv6uqKZKzOzk61tLSoo6MjkiAEAEC9+fWTi3XhbU/EPY1QjBya0R9mHaHthzfEPRWg7pjZXHdvGui6ml4ZilIpABGEAAAIR+t+O+nyEydqXU8h9LFeeuZxfevyf9cXvvo1/cPEqaGOtXD5at3yx1f04rJVOvC94Z5VCKA6rAxVwcxUyz8vAACSLJfLRVIiHnX5+0vL3lbrtx7UBxpe1i1fuzC0cQAMbEsrQzXdTS5KpV/O7OMBACAcc+bMkZmF/tbS0qLu7m5JfY2RwhzvH8ePkXtBv/5dR8w/YQDvRBiqUC6Xk7sThgAACElbW5vcPfS3jo4OZbNZSVI2m1VHR0e4Y/as1y47DNeB04+O+ScM4J0IQwAAIBGiesGxublZ7e3tkhRZh9jdxwzXmD32Dn0cANUhDAEAgLoTdWOk944Zrr++sSaSsQBUjjAEAADqTtR7gfcYM1yvv7VOazfkIxkvl8vJzCjvB7aC1toAAKDuRNW5rmSPHYdLkn417zXtvN2w0Mcb976DlBm9i8a97yD9/vmloY+XSaU0bc8xGpLhdXbUlsSFITM7WtK3JaUl3eDuV8Q8JQAAgHdln51HSZJm/fzp0Mda99pzev2nX5Hne3ThmSdp51O+rqG77hf6uFPSr+lXXz8v9HGAICUqDJlZWtJ3JX1Q0iJJj5nZ3e4+P96ZAQAAbLv9xm+nxTecLxs6PPSx8qtWyHvWS5K8Z72W/epKpUftGOqYOxx6mh5/zz+qe31e2SHpUMcCgpSoMCRpmqQX3f1lSTKzn0o6QRJhCAAA1LQvX/DpeA6V/e3doTeKeGzhGzr5e5069Qd/1OjhDaGOFTUz02H/OFbHTh4vs7hnE45hDWmNGpqRDdZvcAuSFoZ2lfRqv88XSTo4prkAAAAEJsrW4aeeeqpuvPFGnXrqqZF0zGt672idOm0PPbt4pVasXh/6eFFasz6vy/53vi7738H92nwmZUptSxjaxvy0rbEr6Lxm7h7sPb4LZnaypA+5+7nFz8+UNM3dL+x3m/MklQpS95X0fIRTHCtpeYTjIVg8frWLx6628fjVNh6/2sVjV9t4/ILzXncfN9AVSVsZWiRp936f7yZpcf8buPv3JX0/ykmVmFmXuzfFMTbePR6/2sVjV9t4/Gobj1/t4rGrbTx+0Uha/8PHJO1tZnua2RBJp0i6O+Y5AQAAABiEErUy5O49ZnaBpN+qt7X2je7+bMzTAgAAADAIJSoMSZK73yPpnrjnsRmxlOchMDx+tYvHrrbx+NU2Hr/axWNX23j8IpCoBgoAAAAAEJWk7RkCAAAAgEgQhipkZkeb2fNm9uL/3969x15d13Ecf74ETbmomeUULC+ZSE7R0pmUeEFn5USbNslSF60kRbQsdbbW2myYl7TpcoqGGdkIrzMvMCIoh4qhAkpKqRMcCc7yVl6QV398P792OvDjdw4Dv/x+5/XYzs75fn7fy+t7vvudc97n8/meaNnRFQAAB7JJREFUr6QL684TrZG0q6TZkpZIelLSxLozRfsk9ZP0mKR76s4SrZO0vaTpkv5a/gc3/cVOYqORdF553Vws6VZJW9edKbon6SZJKyUtbmjbQdJMSUvL/QfrzBjr1s2xu6y8di6UdIek7evM2JelGGqBpH7AtcDngeHAWEnD600VLVoNfNf2PsAhwFk5dr3SRGBJ3SGibVcD99seBuxPjmGvIWkIcA7wadv7Uv2o0Sn1pooeTAGObWq7EJhley9gVpmOzc8U1j52M4F9be8HPANc9H6H6hQphlpzMPA328/afgf4LTCm5kzRAtsrbC8oj1+n+jA2pN5U0Q5JQ4EvApPrzhKtk7QtcBhwI4Dtd2z/q95U0ab+wDaS+gMDaLruX2xebM8FXmlqHgPcXB7fDJzwvoaKlqzr2NmeYXt1mXyI6tqbsQmkGGrNEGBZw/Ry8oG615G0G3AA8HC9SaJNVwHfB9bUHSTasgewCvhlGeI4WdLAukNFa2y/CFwOvACsAF61PaPeVLEBdrK9AqovB4GP1JwnNszXgfvqDtFXpRhqjdbRlp/h60UkDQJuA861/VrdeaI1ko4DVtr+S91Zom39gQOBX9g+AHiTDNHpNcq5JWOA3YFdgIGSvlpvqojOI+liqiH/U+vO0lelGGrNcmDXhumhZLhAryFpS6pCaKrt2+vOE20ZCRwv6Xmq4alHSvp1vZGiRcuB5ba7emKnUxVH0TuMBp6zvcr2u8DtwKE1Z4r2vSRpZ4Byv7LmPNEGSacDxwGnOtfC2WRSDLVmPrCXpN0lbUV1EundNWeKFkgS1TkLS2xfWXeeaI/ti2wPtb0b1f/dH2zn2+lewPY/gGWS9i5NRwFP1Rgp2vMCcIikAeV19CjyAxi90d3A6eXx6cBdNWaJNkg6FrgAON72v+vO05elGGpBOYHtbOABqjeDabafrDdVtGgk8DWqHoXHy+0LdYeK6BATgKmSFgIjgJ/UnCdaVHr0pgMLgEVUnxeurzVUrJekW4F5wN6SlksaB0wCjpa0FDi6TMdmpptjdw0wGJhZPrtcV2vIPkzpdYuIiIiIiE6UnqGIiIiIiOhIKYYiIiIiIqIjpRiKiIiIiIiOlGIoIiIiIiI6UoqhiIiIiIjoSCmGIiI6nCRLuqJh+nxJP9pI654i6aSNsa4etnOypCWSZm/qbUVERN+RYigiIt4GviRpx7qDNJLUr43ZxwHftn3EpsoTERF9T4qhiIhYTXVBzfOa/9DcsyPpjXJ/uKQ5kqZJekbSJEmnSnpE0iJJezasZrSkP5X5jivL95N0maT5khZK+lbDemdL+g3VxT6b84wt618s6dLS9kPgs8B1ki5rmn9nSXPLRQsXS/pc436UxydJmlIe7yTpDklPlNuhpf20kvMJSbeUtg9Luq3sw3xJI0v7qIaLPD8mafB6chwjaZ6kBZJ+J2lQaZ8k6amyzcvbOJYREdGG/nUHiIiIzcK1wEJJP21jmf2BfYBXgGeBybYPljQRmACcW+bbDRgF7AnMlvRx4DTgVdsHSfoA8KCkGWX+g4F9bT/XuDFJuwCXAp8C/gnMkHSC7R9LOhI43/ajTRm/Ajxg+5LS0zSgh336OTDH9oll/kGSPglcDIy0/bKkHcq8VwM/s/1nSR8FHijPx/nAWbYfLMXNW8A3m3OUnrgfAKNtvynpAuA7kq4BTgSG2bak7XvIHBERGyjFUEREYPs1Sb8CzgH+0+Ji822vAJD0d6CrmFkENA5Xm2Z7DbBU0rPAMOAYYL+GXqftgL2Ad4BHmguh4iDgj7ZXlW1OBQ4D7lxfRuAmSVsCd9p+vId9OpKqUMP2e8Crkk4Dptt+ubS/UuYdDQyX1LXstpIGAw8CV5Z8t9teLmmtHJJGAcOpCkGArYB5wGtUBdRkSb8H7ukhc0REbKAMk4uIiC5XUZ17M7ChbTXlvULVJ/atGv72dsPjNQ3Ta/j/L9vctB0DAibYHlFuu9vuKqbe7Cafumnvlu25VAXTi8AtpbBpzrR1D6sRa+8DVM/LZxr2YYjt121PAr4BbAM8JGlYNzkEzGxYfrjtcbZXU/WO3QacANzf7n5HRERrUgxFRATwvx6PaVQFUZfnqYalAYwBttyAVZ8saYtyHtEewNNUQ8rGl54SJH1C0sD1rQR4GBglaccy1GwsMGd9C0j6GLDS9g3AjcCB5U8vSdpH0hZUQ9K6zALGl2X7Sdq2tH1Z0odKe9cwuRnA2Q3bGlHu97S9yPalwKPAsG5yPASMLMMGkTSgPA+DgO1s30s11HBED89LRERsoAyTi4iIRlfQ8AEfuAG4S9IjVEVBd7026/M0VdGyE3Cm7bckTaY6l2hB6XFaRdUL0i3bKyRdBMym6lW51/ZdPWz7cOB7kt4F3qAMgQMupBp+tgxYDAwq7ROB6yWNA94DxtueJ+kSYI6k94DHgDOohhReK2kh1fvpXOBM4FxJR5TlnwLuA05pzmF7laQzgFvLeVNQnUP0OtVzvnXZz7V+2CIiIjYO2evq+Y+IiIiIiOjbMkwuIiIiIiI6UoqhiIiIiIjoSCmGIiIiIiKiI6UYioiIiIiIjpRiKCIiIiIiOlKKoYiIiIiI6EgphiIiIiIioiOlGIqIiIiIiI70X3GmxucZ6VnyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xaxis = np.linspace(-0.5, N_trials+0.5, 1000) # This way we include all possibilties!\n", "yaxis = func_binomial(np.floor(xaxis+0.5), 1000, 20, 0.2)\n", "ax.plot(xaxis, yaxis, '-', label='Binomial distribution')\n", "ax.legend()\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting a Poisson:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def func_poisson(x, N, mu) :\n", " return N * poisson.pmf(x, mu)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAHgCAYAAABn17aGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3yU5Z3///d1zwQIBxEFracW7aoVA0SMQKJyMK5atSpaa+sJtNbDr+p23Vq03SWD7VZbrbquVr+tVdCqtdpK3Wrtajwgm7QaFFEB66FQEcpJRYRwyNyf3x9zMEiACb0Pc3g9H488Jpm5Z+5rQoD7nc91fS5nZgIAAACASuPFPQAAAAAAiANhCAAAAEBFIgwBAAAAqEiEIQAAAAAViTAEAAAAoCIRhgAAAABUpGTcA/hHDBw40AYPHhz3MAAAAAAUqdmzZ680s0FdPVbSYWjw4MFqa2uLexgAAAAAipRzbtHWHmOaHAAAAICKRBgCAAAAUJEIQwAAAAAqUkmvGQIAAEC0Nm3apMWLF2v9+vVxDwXYTK9evbT33nurqqqq4OcQhgAAAFCwxYsXq1+/fho8eLCcc3EPB5AkmZlWrVqlxYsXa9999y34eUyTAwAAQMHWr1+vXXfdlSCEouKc06677trtiiVhCAAAAN1CEEIx2pGfS8IQAAAASkoikVBtba2GDx+uESNGqKWlRZK0ZMkSffnLXw713G1tbbr88su3ecyzzz6rE088cbuvNW7cuPyemccff7w+/PDDrR578803a926dVt9/IILLtC8efMkSX379t3uuTubM2eOHn/88fzXjz76qK677rpuvUapYs0QAAAASkp1dbXmzJkjSfrjH/+oq6++Ws8995z23HNPPfzww6Geu66uTnV1dYG/bucw0pWbb75ZZ599tnr37r3FY+l0WnfeeecOn3vOnDlqa2vT8ccfL0k66aSTdNJJJ+3w65USKkMAAAAoWR999JEGDBggSVq4cKFqamokSdOmTdOpp56q4447Tvvvv7++853v5J/zwAMPaOjQoaqpqdHkyZPz9/ft21eTJ0/WoYceqqOPPlovvPCCxo0bp/3220+PPvqopM2rPi+88IIaGhp0yCGHqKGhQW+88cY2x9re3q6vfvWrGjZsmM444wy1t7fnHxs8eLBWrlyptWvX6oQTTtDw4cNVU1OjBx98ULfccouWLFmi8ePHa/z48fmxTpkyRaNGjVJra+tmVSZJ+rd/+zeNGDFCjY2NWrFihaTNK1ErV67U4MGDtXHjRk2ZMkUPPvigamtr9eCDD2ratGm69NJLJUmLFi1SY2Ojhg0bpsbGRv3tb3+TJE2aNEmXX365GhoatN9++4UeQsNCZQgAAAA7ZOr/vK55Sz4K9DWH7LmTmr508DaPaW9vV21trdavX6+lS5fq6aef7vK4OXPm6OWXX1bPnj114IEH6rLLLlMikdDkyZM1e/ZsDRgwQMccc4xmzJihU045RWvXrtW4ceP0ox/9SBMmTNC///u/68knn9S8efM0ceLELaolX/jCFzRz5kwlk0k99dRT+u53v6vf/OY3Wx337bffrt69e2vu3LmaO3euRowYscUxTzzxhPbcc0899thjkqTVq1erf//+uvHGG/XMM89o4MCBkqS1a9eqpqZG11xzzRavsXbtWo0YMUI/+clPdM0112jq1Km69dZbuxxTjx49dM0116itrS1/zLRp0/KPX3rppTr33HM1ceJE3XXXXbr88ss1Y8YMSdLSpUs1a9YsLViwQCeddFLoUxTDQGUIAAAAJSU3TW7BggV64okndO6558rMtjiusbFR/fv3V69evTRkyBAtWrRIL774osaNG6dBgwYpmUzqrLPO0syZMyVlgsFxxx0nSRo6dKjGjh2rqqoqDR06VAsXLtzi9VevXq3TTz9dNTU1+td//Ve9/vrr2xz3zJkzdfbZZ0uShg0bpmHDhm1xzNChQ/XUU09p8uTJev7559W/f/8uXyuRSOi0007r8jHP83TGGWdIks4++2zNmjVrm+PaltbWVp155pmSpHPOOWez1zrllFPkeZ6GDBmiZcuW7fA54kRlCAAAADtkexWcKNTX12vlypX5qWCd9ezZM/95IpFQR0dHl6Epp6qqKt+RzPO8/PM9z1NHR8cWx//Hf/yHxo8fr0ceeUQLFy7UuHHjtjve7XU8O+CAAzR79mw9/vjjuvrqq3XMMcdoypQpWxzXq1cvJRKJ7Z6v8zmTyaR835ekHd40t/P4O39/t/V9LWZUhgAAAFCyFixYoHQ6rV133bWg40eNGqXnnntOK1euVDqd1gMPPKCxY8fu0LlXr16tvfbaS9LmU8u2ZsyYMbrvvvskSa+99prmzp27xTFLlixR7969dfbZZ+vb3/62XnrpJUlSv379tGbNmoLG5ft+fg3P/fffryOOOEJSZl3S7NmzJWmzNT7beu2Ghgb96le/kiTdd999+dcqF1SGAAAAUFJya4akTEVi+vTpBVdJ9thjD1177bUaP368zEzHH3+8Tj755B0ax3e+8x1NnDhRN954o4466qjtHn/JJZfovPPO07Bhw1RbW6uRI0duccyrr76qK6+8Up7nqaqqSrfffrsk6cILL9QXv/hF7bHHHnrmmWe2eZ4+ffro9ddf16GHHqr+/fvrwQcflCR9+9vf1le+8hXde++9m413/Pjxuu6661RbW6urr756s9e65ZZbdP755+v666/XoEGDdPfdd2/3fZYSV6olLUmqq6uzzl0zAAAAEK758+froIMOinsYQJe6+vl0zs02sy77oTNNDgAAAEBFIgwBAAAAqEiEIQAAAAAVKbQw5Jzbxzn3jHNuvnPudefcv2TvTznn3nPOzcl+HN/pOVc7595yzr3hnDs2rLEBAAAAQJjd5Dok/ZuZveSc6ydptnPuyexjN5nZDZ0Pds4NkfRVSQdL2lPSU865A8wsHeIYAQAAAFSo0CpDZrbUzF7Kfr5G0nxJe23jKSdL+pWZbTCzv0p6S9KW/QYBAAAAIACRrBlyzg2WdIikP2fvutQ5N9c5d5dzbkD2vr0kvdvpaYu17fAEAACACpRIJFRbW6uamhqdfvrpWrdu3TaPb2hoiGhkXZs0aVJ+k9MLLrhA8+bN2+qx06ZN05IlS7b6+JQpU/TUU09JymyiunLlyoLHsXDhQt1///35r9va2nT55ZcX/PxyFHoYcs71lfQbSd8ys48k3S7p85JqJS2V9JPcoV08fYtNkJxzFzrn2pxzbStWrAhp1AAAAChW1dXVmjNnjl577TX16NFDd9xxxzaPb2lpiWhk23fnnXdqyJAhW318W2EonU7rmmuu0dFHH71D5/50GKqrq9Mtt9yyQ69VLkINQ865KmWC0H1m9ltJMrNlZpY2M1/Sz/XJVLjFkvbp9PS9JW3xk2BmPzOzOjOrGzRoUJjDBwAAQJE78sgj9dZbb0mSbrzxRtXU1KimpkY333xz/pi+fftKkpYuXaoxY8bkq0rPP/+80um0Jk2apJqaGg0dOlQ33XSTJGnOnDkaPXq0hg0bpgkTJuiDDz6QJI0bN06TJ0/WyJEjdcABB+j555/fYkxmpksvvVRDhgzRCSecoOXLl+cfGzdunNra2ro878MPP6y2tjadddZZqq2tVXt7uwYPHqxrrrlGRxxxhB566KHNqkySdP3112vkyJEaOXJk/vvw6WNy7/+qq67S888/r9raWt1000169tlndeKJJ0qS3n//fZ1yyikaNmyYRo8erblz50qSUqmUzj//fI0bN0777bdf2YWn0BooOOecpF9Imm9mN3a6fw8zW5r9coKk17KfPyrpfufcjco0UNhf0gthjQ8AAAD/oD9cJf391WBf8zNDpS9eV9ChHR0d+sMf/qDjjjtOs2fP1t13360///nPMjONGjVKY8eO1SGHHJI//v7779exxx6r733ve0qn01q3bp3mzJmj9957T6+9lrkk/fDDDyVJ5557rv77v/9bY8eO1ZQpUzR16tR8wOro6NALL7ygxx9/XFOnTs1PW8t55JFH9MYbb+jVV1/VsmXLNGTIEJ1//vmbHdPVeXfeeWfdeuutuuGGG1RXV5c/tlevXpo1a5Yk6YknntjsdXbaaSe98MILuueee/Stb31Lv//977f6/bruuut0ww035I959tln8481NTXpkEMO0YwZM/T000/r3HPP1Zw5cyRJCxYs0DPPPKM1a9bowAMP1CWXXKKqqqrt/OmUhjArQ4dLOkfSUZ9qo/1j59yrzrm5ksZL+ldJMrPXJf1a0jxJT0j6Jp3kAAAA8Gnt7e2qra1VXV2dPvvZz+rrX/+6Zs2apQkTJqhPnz7q27evTj311C2qNocddpjuvvtupVIpvfrqq+rXr5/2228/vfPOO7rsssv0xBNPaKeddtLq1av14YcfauzYsZKkiRMnaubMmfnXOfXUUyVJhx56qBYuXLjF+GbOnKmvfe1rSiQS2nPPPXXUUUdtcUxX592aM844Y6uPfe1rX8vftra2bv2bth2zZs3SOeecI0k66qijtGrVKq1evVqSdMIJJ6hnz54aOHCgdtttNy1btmyHz1NsQqsMmdksdb0O6PFtPOc/Jf1nWGMCAABAgAqs4AQtt2aoM7MtlppvYcyYMZo5c6Yee+wxnXPOObryyit17rnn6pVXXtEf//hH3Xbbbfr1r3+dnyq3NT179pSUaeTQ0dHR5TGZSVJbN2DAgC3Oe9ddd3V5bJ8+fbb6Op3Pk/s8mUzK931Jme/Lxo0btzmW3HFbe+3c+5W2/Z5LUSTd5AAAAIAwjRkzRjNmzNC6deu0du1aPfLIIzryyCM3O2bRokXabbfd9I1vfENf//rX9dJLL2nlypXyfV+nnXaavv/97+ull15S//79NWDAgHxl6d57781XiQody69+9Sul02ktXbpUzzzzzBbHdHVeSerXr5/WrFlT8LkefPDB/G19fb2kTJe52bNnS5J+97vfadOmTdt97TFjxui+++6TlJk+N3DgwG1Wq8pFmJuuAgAAAJEYMWKEJk2apJEjM725Lrjggs3WC0mZi/zrr79eVVVV6tu3r+655x699957Ou+88/KVlGuvvVaSNH36dF188cVat26d9ttvP919990Fj2XChAl6+umnNXToUB1wwAFdBqmtnXfSpEm6+OKLVV1dXdC0tw0bNmjUqFHyfV8PPPCAJOkb3/iGTj75ZI0cOVKNjY35ytKwYcOUTCY1fPhwTZo0abPvTyqV0nnnnadhw4apd+/emj59esHvt5S5QkqKxaqurs7a2triHgYAAEDFmD9/vg466KC4hwF0qaufT+fcbDOr6+p4pskBAAAAqEiEIQAAAAAViTAEAAAAoCIRhgAAAABUJMIQAAAAgIpEGAIAAECoUqmUnHNKpVJxDwXYDGEIAAAAocqFoKDCUCKRUG1trQ4++GANHz5cN954Y36/nra2Nl1++eVbfe7ChQt1//33b/XxJUuW6Mtf/rIkadq0abr00ku7NbZp06ZpyZIl+a8vuOACzZs3r1uv0R0bNmzQ0Ucfrdra2vwGrNuzbt06nXXWWRo6dKhqamp0xBFH6OOPPw5tjMWMTVcBAABQUqqrqzVnzhxJ0vLly3XmmWdq9erVmjp1qurq6lRX1+WWMpI+CUNnnnnmFo91dHRozz331MMPP7zDY5s2bZpqamq05557SpLuvPPOHX6tQrz88svatGlT/vtRiP/6r//S7rvvrldffVWS9MYbb6iqqiqsIRY1KkMAAAAIVWtr62a3Qdptt930s5/9TLfeeqvMTM8++6xOPPFESdJzzz2n2tpa1dbW6pBDDtGaNWt01VVX6fnnn1dtba1uuukmTZs2Taeffrq+9KUv6ZhjjtHChQtVU1OTf/13331Xxx13nA488EBNnTpVkrY45oYbblAqldLDDz+strY2nXXWWaqtrVV7e7vGjRuntrY2SdIDDzyQr8ZMnjw5//y+ffvqe9/7noYPH67Ro0dr2bJlW7zP999/X6eccoqGDRum0aNHa+7cuVq+fLnOPvtszZkzR7W1tXr77bc3e864ceM0efJkjRw5UgcccICef/55SdLSpUu111575Y878MAD1bNnz62+L0l66623dPTRR2v48OEaMWJE/lw//vGPNXToUA0fPlxXXXWVJOntt9/Wcccdp0MPPVRHHnmkFixYIEl66KGHVFNTo+HDh2vMmDGSpNdff10jR45UbW2thg0bpjfffFOS9Mtf/jJ//0UXXaR0Oq10Oq1JkyappqZGQ4cO1U033VTYD8m2mFnJfhx66KEGAACA6MybN2+7xzQ1NeU/b2lpserqapNk1dXV1tLS0uVx3dGnT58t7tt5553t73//uz3zzDN2wgknmJnZiSeeaLNmzTIzszVr1timTZs2e9zM7O6777a99trLVq1aZWZmf/3rX+3ggw/OP/aZz3zGVq5caevWrbODDz7YXnzxxc2OMTO7/vrr8+9l7Nix9uKLL+Yfy3393nvv2T777GPLly+3TZs22fjx4+2RRx4xMzNJ9uijj5qZ2ZVXXmnf//73t3h/l156qaVSKTMza25utuHDh5uZbfF+Ohs7dqxdccUVZmb22GOPWWNjo5mZvfzyyzZo0CAbPXq0fe9737O//OUvW7z3T7+vkSNH2m9/+1szM2tvb7e1a9fa448/bvX19bZ27Vozs/z38Kijjsq/5p/+9CcbP368mZnV1NTY4sWLzczsgw8+yL+vX/7yl2ZmtmHDBlu3bp3NmzfPTjzxRNu4caOZmV1yySU2ffp0a2trs6OPPjo/vtxrdNbVz6ekNttKnqAyBAAAgEBNnTpVzjk559TQ0KD29nZJUnt7uxoaGvKP5SotQchc827u8MMP1xVXXKFbbrlFH374oZLJrleI/PM//7N22WWXrT626667qrq6WqeeeqpmzZq1Q+N78cUXNW7cOA0aNEjJZFJnnXWWZs6cKUnq0aNHvpp16KGHauHChVs8f9asWTrnnHMkSUcddZRWrVql1atXb/e8p5566havW1tbq3feeUdXXnml3n//fR122GGaP3/+Vl9jzZo1eu+99zRhwgRJUq9evdS7d2899dRTOu+889S7d29J0i677KKPP/5YLS0tOv300/NVnaVLl0rK/HlMmjRJP//5z5VOpyVJ9fX1+uEPf6gf/ehHWrRokaqrq9Xc3KzZs2frsMMOU21trZqbm/XOO+9ov/320zvvvKPLLrtMTzzxhHbaaaftvv/tIQwBAAAgUE1NTfnfvLe0tKi6ulpSZq1PS0tL/rGmpqZAzvfOO+8okUhot9122+z+q666Snfeeafa29s1evTo/HStT+vTp89WX9s5t8XXyWQy37BBktavX7/dMXYV1nKqqqry50kkEuro6Cjo+Z8eW1d69uzZ5ev27dtXp556qn7605/q7LPP1uOPP77V97W1sZvZFmPwfV8777yz5syZk//IBa077rhDP/jBD/Tuu++qtrZWq1at0plnnqlHH31U1dXVOvbYY/X000/LzDRx4sT889944w2lUikNGDBAr7zyisaNG6fbbrtNF1xwwXbf//YQhgAAABCozl3j6uvr1dzcLElqbm5WfX19l8ftqBUrVujiiy/WpZdeusWF+dtvv62hQ4dq8uTJqqur04IFC9SvXz+tWbOm4Nd/8skn9f7776u9vV0zZszQ4Ycfrt13313Lly/XqlWrtGHDBv3+97/PH7+11x81apSee+45rVy5Uul0Wg888IDGjh1b8DjGjBmj++67T5L07LPPauDAgTtcGfm///s/ffDBB5KkjRs3at68efrc5z631fe10047ae+999aMGTMkZTrYrVu3Tsccc4zuuusurVu3TlJmXdNOO+2kfffdVw899JCkTGB65ZVXJGX+PEaNGqVrrrlGAwcO1Lvvvpuv+Fx++eU66aSTNHfuXDU2Nurhhx/W8uXL86+7aNEirVy5Ur7v67TTTtP3v/99vfTSSzv0/jujmxwAAABClQtAnYPQP6K9vV21tbXatGmTksmkzjnnHF1xxRVbHHfzzTfrmWeeUSKR0JAhQ/TFL35RnucpmUxq+PDhmjRpkgYMGLDNcx1xxBE655xz9NZbb+nMM8/Md6qbMmWKRo0apX333Vdf+MIX8sdPmjRJF198saqrqzdrGLHHHnvo2muv1fjx42VmOv7443XyyScX/J5TqZTOO+88DRs2TL1799b06dMLfu6nvf3227rkkktkZvJ9XyeccIJOO+00Oee2+r7uvfdeXXTRRZoyZYqqqqr00EMP6bjjjtOcOXNUV1enHj166Pjjj9cPf/hD3Xfffbrkkkv0gx/8QJs2bdJXv/pVDR8+XFdeeaXefPNNmZkaGxs1fPhwXXfddfrlL3+pqqoqfeYzn9GUKVO0yy676Ac/+IGOOeYY+b6vqqoq3XbbbaqurtZ5552Xr15de+21O/w9yHHbKtkVu7q6Ost15wAAAED45s+fr4MOOqjbz3PObXOqGBCErn4+nXOzzazLfutMkwNQ1Ni1HABKX9CbrgJBoTIEoOjx20QAKB47WhkCokBlCAAAAAAKQBgCAABAt1CtRzHakZ9LwhAAAAAK1qtXL61atYpAhKJiZlq1apV69erVrefRWhsAAAAF23vvvbV48WKtWLEi7qEAm+nVq5f23nvvbj2HMAQAAICCVVVVad999417GEAgmCYHAAAAoCIRhgAAAABUJMIQAAAAgIpEGAJKXCqVknOOXb0BAAC6yZVyW8S6ujpra2uLexhA7JxzZd3itNzfHwAACI9zbraZ1XX1GJUhAAAAABWJMAQAAACgIhGGAAAAAFQkwhAAAACAikQYAgAAAFCRCEMAAAAAKhJhCAAAAEBFIgwBAAAAqEiEIQAAAAAViTAEAAAAoCIRhgAUtdbW1s1uAQAAgkIYAtBtqVQqkvO0traqsbFRktTY2BhZIIrq/QEAgHg5M4t7DDusrq7O2tra4h4GEDvnnKL8u+yci+xccSnlfxsBAMAnnHOzzayuq8eoDAHotqamJplZ6B8tLS2qrq6WJFVXV6ulpSWS8zY1NcX8HQYAAFGgMgSUgagrQ1FqbW1VQ0ODWlpaVF9fH/dwAABAiaEyBKBk5QIQQQgAAAQtGfcAAJSeP7+zSjc++RdFVYza/czr9JU7ommeMHTv/vqPE4dEci4AABAvwhCAbmtesFxtiz7QyMG7RHNCP62EF37ThkWr1uq1JasJQwAAVAjCEIBuW7N+k3bp00MPXDg6kvO5i+r1wAPhl6Gu+8MC3fV/fw39PAAAoDiwZghAt61Z36F+PcvvdylJzyntl2cjCgAAsCXCEFDichuRRrUhqSR9vKFD/XqVXxhKZMNQuXbmAwAAmyMMASFJpVKhn6O1tVWNjY2SpMbGxsgC0by/vKO+ZRiGktl1SR1UhwAAqAiEISAkU6dOlXMu1I+Ghga1t7dLktrb29XQ0BD6OZ1zWrx8lfqW4TS5RCIThpgqBwBAZSAMASFpamqSmYX60dLSourqaklSdXW1WlpaQj+nmWnAoD3Ut2dVzN/h4FV5mX8SqQwBAFAZCENASKKYJldfX6/m5mZJUnNzc2QbkyZ69SnbNUOSlE4ThgAAqASEIaDE5QJQVEHIzMq2gUIykVsz5Mc8EgAAEAXCEIBuWbcxLd9UnmuGaKAAAEBFIQwB6JaPN3RIEt3kAABAySMMAeiWNeszYahfr/JroJDMNlBgzRAAAJWBMASgW9as3yRJ6leG0+RYMwQAQGUhDAHolnKeJpfvJsc0OQAAKgJhCEC3fJyfJld+YYg1QwAAVBbCEIBuWZOrDJXhNLlEbtNV1gwBAFARCEMAuiXfQKFnOTZQYM0QAACVhDAEoFty0+T69EzEPJLg5RoosGYIAIDKQBgC0C0fb9ik3j0SSibK758PNl0FAKCylN/VDIBQfbyhoyzXC0md9hkiDAEAUBHK84oGiNk9rQv19ILlkZ1vty+nNOnuFyI517wlH0XaVjuVSuVvc5+HJVcZ2pSObs1QKpXS1KlT1dTUFPr7AwAAm3Nmpfsb0Lq6Omtra4t7GMAWjr1pppaubte+A/tEcr4XXnhRI0ceFsm5JOnog3bXZY37R3a+qLzy7oc6+bb/0y8m1qnxoN0jO69zTqX8bzEAAMXMOTfbzOq6eozKEBCCtJmO2H+gfnrWoZGcz7kj9bt7uJj+R+UaKLBmCACAysCaISAEvpk85+IeBrqJNUMAAFQWwhAQAt8nDJUiuskBAFBZCENACNJm+QtrlI78pqsRNlAAAADxCS0MOef2cc4945yb75x73Tn3L9n7d3HOPemcezN7OyB7v3PO3eKce8s5N9c5NyKssQFh831RGSpBVIYAAKgsYVaGOiT9m5kdJGm0pG8654ZIukpSs5ntL6k5+7UkfVHS/tmPCyXdHuLYgFBl1gzFPQp0V1WCNUMAAFSS0MKQmS01s5eyn6+RNF/SXpJOljQ9e9h0SadkPz9Z0j2W8SdJOzvn9ghrfECY0j7T5EoRlSEAACpLJGuGnHODJR0i6c+SdjezpVImMEnaLXvYXpLe7fS0xdn7gJLjm+RFFIY6b0qKf0xuzVCaNUMAAFSE0MOQc66vpN9I+paZfbStQ7u4b4tfzzrnLnTOtTnn2lasWBHUMIFARTlNLpVKycwIQwFIsM8QAAAVJdQw5JyrUiYI3Wdmv83evSw3/S17uzx7/2JJ+3R6+t6Slnz6Nc3sZ2ZWZ2Z1gwYNCm/wwD8g7ZsSNFAoOUmmyQEAUFHC7CbnJP1C0nwzu7HTQ49Kmpj9fKKk33W6/9xsV7nRklbnptMBpcY3i2yaHILDpqsAAFSWZIivfbikcyS96pybk73vu5Kuk/Rr59zXJf1N0unZxx6XdLyktyStk3ReiGMDQsWmq6Xpk32GCEMAAFSC0MKQmc1S1+uAJKmxi+NN0jfDGg8QJTZdLU2e5+SclPZpoAAAQCWIpJscUGnYdLV0JT2nTUyTAwCgIoQ5TQ6oWHvZ33XgB4uk+W/GPRR000HeEqX9feMeBgAAiABhCAjBTcn/Vu1f3pb+EvdI0F13e/11W/rIuIcBAAAiQEt/NMgAACAASURBVBgCAmZm6qP1WtR/pD73tRu3/wQUj9afqs8rD7FmCACACkEYAgLmm+TJ1/qq/tJnhsY9HHRHv88oqTT7DAEAUCFooAAELO2bkkpLHr9rKDmJKiXlq6ODyhAAAJWAMAQEzDdTwvkyl4h7KOgur0qS5PsdkZ2ytbV1s1sAABAdwhAQMN8ylSGjMlR6vEyAfe2VlyM5XWtrqxobM9uuNTY2RhaIUqlUJOcBAKDYcbUGBCztmxLyJSpDpSeRqQwtmP+6XMT7RLW3t6uhoSGy8xGIAACgMgQEzvelhNIyjzBUcrLT5GoOOkBmFvpHS0uLqqurJUnV1dVqaWmJ5LxNTU1xfpcBACgahCEgYJlpclSGSlI2wB7wT5+P5HT19fVqbm6WJDU3N6u+vj6S81IVAgAggzAEBCxtlq0MMQu15GSnyVl6U2SnzAWgqIIQAAD4BGEICJifWzPENLnSk50mJz+6MAQAAOJDGAIC5puyYYjKUMnJVoaUjq61NgAAiA9hCAhY2vdV5dKsGSpF2WqeRbjPEAAAiA9hCAiYn05LEmuGSlF2mpwjDAEAUBEIQ0DA/OyFtGPNUOnJTpNzETZQAAAA8SEMAQHzO7JVBSpDpSf3Z0YDBQAAKgJhCAiYn7uQJgyVntyfmaXjHQcAAIgEYQgImKWpDJWs3DQ5KkMAAFQEwhAQsFwDBTn+epWcbID1aKAAAEBF4GoNCJifW3yfoDJUcrLd5GitDQBAZSAMAUHLVoYc0+RKT4LKEAAAlYQwBASMBgolLPtn5owwBABAJSAMAQHz85Uh9hkqOflNV+kmBwBAJSAMAUHLdpNzrBkqPblpckY3OQAAKgFhCAhYvoEC0+RKT26aHGuGAACoCIQhIGC+TwOFkpWbJsemqwAAVATCEBC0bFWBNUMlKLvpaoIGCgAAVATCEBAwY81Q6cpW8xLWITOLeTAAACBshCEgYPkNOwlDpScbhpJKK+0ThgAAKHeEISBo2TVDHmuGSk9umpzS6ogoDKVSqc1uAQBAdLhaAwLm56bJEYZKT7aBQlXEYYggBABAPKgMAQHLrRlSggYKJcfzZHJKurTSaabJAQBQ7ghDQMBye9R42SlXKC2+SyqptDp8P+6hAACAkDGPBwiYn92jhm5ypcn3MmHotSUfadc+PeIeTuD69Exq34F94h4GAABFgas1IGhp9hkqadnK0MS7Xoh7JKF57PIjdPCe/eMeBgAAsSMMAUFjmlxJS/booWP23VV7H1IX91AC95dla3T9H9/Q6nWb4h4KAABFgTAEBCy3zxCVodLkvCrt1S+pvYbsHvdQAte/OhPQ2UIJAIAMGigAAbPcPkNUhkqTl8xPdSw3iey/+GkjDQEAIBGGgMB90k2OwmtJSiTzUx3LjeecJMmnNAQAgCTCEBA8n25yJc2rkvzyXFOT8DJhKE0YAgBAEmEICJxRGSptiSopXZ5hKFcZYpocAAAZhCEgYM5ya4YIQyXJS+Sre+UmF4aMMAQAgCTCEBA8WmuXtoqYJhfzQAAAKBKEISBouU1XqQyVpjKeJkc3OQAANkcYAgKWmyaXSBKGSpJHNzkAACoFYQgI2Cf7DBGGSlIZhyG6yQEAsDnCEBA01gyVtjKeJpevDDFNDgAASYQhIHD5aXKJRMwjwQ7xkmXbQMHzCEMAAHRGGAIC5vy0Nlki/1t4lBgvWbattROObnIAAHRGGAKCZmml5ckjC5Wmcp4mRzc5AAA2QxgCAub8DnUoIUdlqDSV8TS5BN3kAADYDGEICJp1yOevVunyqsp3mhxrhgAA2AxXbEDAnJ9WWjRPKFmJZNlOk3OO1toAAHRGGAIC5sxXB3+1Slc5T5OjMgQAwGa4YgMC5qyDylAp86rKd9NVuskBALAZwhAQMOen5Tv+apWsRJWULs8wlOsmR2UIAIAMrtiAgFEZKnHlPE2ObnIAAGyGMAQEzBkNFEqaV74NFHIbAbPPEAAAGYQhIGAuu+kqSlSiSpJJfvktrPE8KkMAAHTGFRsQMNYMlTgvmbkt16lynqMyBABAFldsQMCc0vKZJle6cmGoTKfKJZyjmxwAAFmEISBgHmuGSluiKnNbpu21PY9ucgAA5BCGgIBlpskRhkqWV+ZhyDnWDAEAkJWMewBAufGYJlfavOyf3d9apepd4h1LCA52C5X29457GAAAFAXCEBAwZ2mlaaBQuqp3ztw+eHa84wjJQ066Y+1dkmriHgoAALEjDAEB8ywtn79apeugk6Xz/yh1bIh7JMFb9rr0x6vVo2NN3CMBAKAocMUGBMyztHzXM+5hYEclktJnR8c9inAke2Vu/XS84wAAoEgwlwcIWCYMsWYIRSi3HsoIQwAASIQhIHDOfBoooDjlwlCZdsoDAKC7CENAwDylZVSGUIxyP5dMkwMAQBJhCAicZx1KE4ZQjLzMMlFHZQgAAEkhhiHn3F3OueXOudc63Zdyzr3nnJuT/Ti+02NXO+fecs694Zw7NqxxAWHzLC2jtTaKUX7NkB/vOAAAKBJhXrFNk3RcF/ffZGa12Y/HJck5N0TSVyUdnH3OT53jV+soTZ58GY0aUYxylSGjMgQAgBRiGDKzmZLeL/DwkyX9ysw2mNlfJb0laWRYYwPC5LHpKopV7ufSpzIEAIAUz5qhS51zc7PT6AZk79tL0rudjlmcvQ8oOQmlZY7KEIoQa4YAANhM1GHodkmfl1Qraamkn2Tvd10ca129gHPuQudcm3OubcWKFeGMEvgHsGYIRYt9hgAA2EykV2xmtszM0mbmS/q5PpkKt1jSPp0O3VvSkq28xs/MrM7M6gYNGhTugIEd4MmXT2UIxSi/ZogwBACAFHEYcs7t0enLCZJyneYelfRV51xP59y+kvaX9EKUYwOCkmCfIRSr7M8lYQgAgIzQfn3tnHtA0jhJA51ziyU1SRrnnKtVZgrcQkkXSZKZve6c+7WkeZI6JH3TjP+tUZo882UeYQhFyCMMAQDQWUFhyDlXY2avbf/IT5jZ17q4+xfbOP4/Jf1nd86BErf6PWn14rhHEbikOuRTGUIxyoUhGigAACCp8MrQHc65HsrsHXS/mX0Y3pBQMe48WlrT5dKwkpaUtMHrE/cwgC1l1wyx6SoAABkFhSEzO8I5t7+k8yW1OedekHS3mT0Z6uhQ3tZ/KA05WRoxMe6RBOo7v31NGweM1IS4BwJ8WrZi6bHpKgAAkrqxZsjM3nTO/bukNkm3SDrEOeckfdfMfhvWAFHG/LQ0YLD0T42hnyqVSmnq1KlqampSKpUK9VxtXkJDEj1DPQewQ/Ld5KgMAQAgFdhNzjk3zDl3k6T5ko6S9CUzOyj7+U0hjg/lzPz8b6rDlgtAYQchSUqbKeF1tXUWEDMv808+DRQAAMgotDJ0qzL7An3XzNpzd5rZkmy1COg+S0tluDmpb6aEIwyhOHUoQRgCACCr0DB0vKT2XLtr55wnqZeZrTOze0MbHcqb+fnuVuXE9yVHGEKRMnnyCEMAAEgqfNPVpyRVd/q6d/Y+YMf42TULZVgZSvumRPm9LZSJtEtkqrIAAKDgMNTLzD7OfZH9vHc4Q0JFyC3gLsP9eHzWDKGI+S6hBA0UAACQVHgYWuucG5H7wjl3qKT2bRwPbFvuN9NlOJ3MN2OaHIqWz5ohAADyCl0z9C1JDznncjtk7iHpjHCGhIqQ+810Ga4ZSvs0UEDx8l2CNUMAAGQVuunqi865L0g6UJKTtMDMNoU6MpQ3P1cZKr/FNb6JaXIoWuY8ORGGAACQurHpqqTDJA3OPucQ55zM7J5QRoXyV85rhnwrx9l/KBO+qAwBAJBTUBhyzt0r6fOS5kj5XymaJMIQdoyVb2UozT5DKGLmPHk+DRQAAJAKrwzVSRpiZhbmYFBBcj9KZbhmiG5yKGa+S8hjmhwAAJIK7yb3mqTPhDkQVJhyXjPEpqsoYr5LyNFaGwAASYVXhgZKmuece0HShtydZnZSKKNC+ctejM18633Nb387klPuNPJU/b/nwj/XJt9n01UULaObHAAAeYWGoVSYg0AFyl6M/eH1ZXpg7oJITjlg/Pm69g/RnOtzu/SJ5DxAd5mYJgcAQE6hrbWfc859TtL+ZvaUc663pPJb7IHoZCtD9f+0m/7jrGMjOWXfvn318ccfh34ezzn1quKvB4qT7xJKEIYAAJBUeDe5b0i6UNIuynSV20vSHZIawxsaylp2zVAikVDvHt3p8L7jbNOGyM4FFCtznhy9cAAAkFR4A4VvSjpc0keSZGZvStotrEGhAkS8z1Bra+tmt0ClMipDAADkFRqGNpjZxtwXzrmkMvsMATsmH4bC7zTQ2tqqxsZMEbOxsTGyQJRKpSI5D9Ad5tFAAQCAnELnDD3nnPuupGrn3D9L+v8k/U94w0LZy4ah386YoRPP+dfITtve3q6GhobIzkcgQrHJVIY2bP9AAAAqQKG/lr9K0gpJr0q6SNLjkv49rEGh/JnfIUk67bQvy8xC/WhpaVF1dbUkqbq6Wi0tLaGf08zU1NQU57cY6JK5hJzYZwgAAKnwbnK+pJ9nP4B/mO/7mXaEEawZqq+vV3NzsxoaGtTc3Kz6+vrQzylRFUJxMpdQUr7MjM2BAQAVr9Bucn9VF2uEzGy/wEeEiuCnO5SQ5LxoGijkAlBUQQgoVrkGCr5JCbIQAKDCFbpmqK7T570kna5Mm21gh/h+ZpqOeeE3UADwCfMSSshX2jclPNIQAKCyFXQlamarOn28Z2Y3Szoq5LGhjFl2nyEXQTc5AJ24TBjy2WsIAICCp8mN6PSlp0ylqF8oI0JF8LNhKKp9hgBkmEvmK0MAAFS6QqfJ/aTT5x2SFkr6SuCjQcXIV4aYJgdEypyXXTNEGAIAoNBucuPDHggqi6Vzm65SGQIi5WUqQz7dtQEAKHia3BXbetzMbgxmOKgUuX2GWDMERMx5SjhfaSpDAAB0q5vcYZIezX79JUkzJb0bxqBQ/nLd5FyCyhAQJfNYMwQAQE6hYWigpBFmtkaSnHMpSQ+Z2QVhDQzlzc9WhkRlCIiWl1CSNUMAAEgqsLW2pM9K2tjp642SBgc+GlSO3IWYV2geBxAEcwl5tNYGAEBS4ZWheyW94Jx7RJJJmiDpntBGhbLnpzOVIY9NH4FoeUkllVY70+QAACi4m9x/Ouf+IOnI7F3nmdnL4Q0L5c5yrawclSEgUl5CnoxucgAAqPBpcpLUW9JHZvZfkhY75/YNaUyoAOwzBMQku2aIbnIAABQYhpxzTZImS7o6e1eVpF+GNSiUPz8fhugmB0Qqu2aIbnIAABReGZog6SRJayXJzJZI6hfWoFD+qAwBMcmuGTIqQwAAFByGNlrmf06TJOdcn/CGhEpgll2wQGUIiJTzEko4U5pFQwAAFByGfu2c+3+SdnbOfUPSU5J+Ht6wUPZylaGI9hlKpVKb3QKVyrLt7NPZjo4AAFSyQrvJ3eCc+2dJH0k6UNIUM3sy1JGhrOXWDHkRVYZSqRRBCNAn6/QsnY55JAAAxG+7v5Z3ziWcc0+Z2ZNmdqWZfZsghH9YvrU20+SAKLlsZciPqDKUSqXknOOXEQCAorTdMGRmaUnrnHP9IxgPKoTRTQ6IR/bvnO9HF4Y63wIAUEwK3fFyvaRXnXNPKttRTpLM7PJQRoWyl8nYkpegmxwQqWxlyFgzBABAwWHosewHEAjLTpNzrtAfQQBByFVjo5omBwBAMdvmlahz7rNm9jczmx7VgFAZctPk5Ll4BwJUmojXDAEAUMy2N0dpRu4T59xvQh4LKkm2MuQlqAwBUcp1cMz/QgIAgAq2vTDU+df2+4U5EFSW/JohjzVDQKRYMwQAQN72rkRtK58D/5DcmiHRTQ6IlEtE200OAIBitr05SsOdcx8pUyGqzn6u7NdmZjuFOjqUr2xlyLHPEBCp3D5DojIEAMC2w5CZcaWKUFh+zRA/YkCU6CYHAMAnWLCBeFg2DLFmCIiUS1RJooECAAASYQgxsdx6BY9uckCUct3kxJohAAAIQ4hJtjKUoDIERCvbzt6nMgQAAGEIMcmuGXJUhoBIOSpDAADkEYYQi9w+Q47KEBCp3EbH7DMEAABhCHHJTtFJ0E0OiJQX8aarra2tm90CAFBMCEOIhZkpbU6ec3EPBagoLpmdmppdtxem1tZWNTY2SpIaGxsjC0SpVCqS8wAASh8LNhAPP620PHlkISBSuTVDN9/0E40//iuRnbe9vV0NDQ2RnY9ABAAoBJUhxMN8mTx5pCEgUrlpcpdfdqnMLNSPlpYWVVdXS5Kqq6vV0tIS+jnNTE1NTXF+iwEAJYTKEOJhucoQYQiIkpfMbLr66Mt/0zULngn9fAee/2O9/cLT+vzIo3T1rPXSrPDO2SPp6aYzaqkKAQAKRhhCPHxfvhzT5ICI9e/dU5I0dI++sn47h3/CfQ7X4kWLNOaIw0M9zbqNaT05b5lee2+1Dt6zf6jnAgCUD8IQYmGWlk9lCIicS2QqQ19v2Ec6+JBIznnL10bov/7n+lDP8ffV6/XkvGVKh98XAgBQRlgzhFg4y1aGKA0B0XK5TVfT8Y4jYLkty9I+aQgAUDjCEOJBNzkgHl55hqFkNg2lfYt5JACAUkIYQjzMsmuGSENApHJhyMorDCWy/5akyUIAgG4gDCEerBkC4pFtrS2/I95xBIxpcgCAHUEYQjzMz4ahuAcCVJgyXTP0yTS5mAcCACgpdJNDPNhnCIhHrjL03I+lF38RySlfvqiPdPsRoZ6jp0wzerRr/sc/lPT5UM8FACgfhCHEw3yZ0U0OiFyfgdJh35A+WhLZKRd++JJqd/5sqOdwm9pV6z2tFR++Jum4UM8FACgfhCHEwplPNzkgDs5JJ9wQ6SknnOlkv7o/3JN8tFS68QuSlddaKABAuFgzhHjk9hlimhyAAOQ2k3VlthYKABCu0MKQc+4u59xy59xrne7bxTn3pHPuzeztgOz9zjl3i3PuLefcXOfciLDGhSLhZ7rJJSgNAQhCtmW4lVmXPABAuMKsDE3TlhO3r5LUbGb7S2rOfi1JX5S0f/bjQkm3hzguFINsZYjCEIBAZBtDOMIQAKAbQgtDZjZT0vufuvtkSdOzn0+XdEqn+++xjD9J2tk5t0dYY0P8PlkzRBoCEID8/klMkwMAFC7qNUO7m9lSScre7pa9fy9J73Y6bnH2PpQr82WEIaDspVKpzW5Dk68MbQr3PACAslIsDRS6uiK2Lg907kLnXJtzrm3FihUhDwthcZZWWo5uckCZS6VSMrPww1B2M1lnVIYAAIWLOgwty01/y94uz96/WNI+nY7bW1KXm2CY2c/MrM7M6gYNGhTqYBEiM/ny5KgMAQiC5yktT2LNEACgG6IOQ49Kmpj9fKKk33W6/9xsV7nRklbnptOhPDlLy7osCALAjulQggYKAIBuCW3TVefcA5LGSRronFssqUnSdZJ+7Zz7uqS/STo9e/jjko6X9JakdZLOC2tcKBa+/KKZpQmgHKSVYJocAKBbQgtDZva1rTzU2MWxJumbYY0FxceZL3OEIQDBScsjDAEAuoWrUcTCGZUhAMHylWDNEACgW7gaRTyyrbUBIChpl5BHZQgA0A1cjSIWznz5TJMDEKA0DRQAAN3E1Shi4UQ3OQDB8pWQZ4QhAEDhCEOIR3afIQAIStrRTQ4A0D1cjSIWnqXpJgcgUD5hCADQTVyNIh5UhgAEjGlyAIDu4moUsfBYMwQgYHSTAwB0F2EI8TBfvkvEPQoAZYRpcgCA7iIMIRbOfCpDAALlK0llCADQLYQhxMLJZFSGAATIZ5ocAKCbCEOIhWesGQIQLMIQAKC7CEOICZUhAMHyXUIJ0U0OAFA4whBiQWUIQNCMyhAAoJsIQ4iFk8ln01UAAWKaHACgu7gaRSycpSXCEIAA+S6phAhDAIDCcTWKWDiZfLFmCEBwzKMyBADoHsIQYuGZL3OsGQIQHKMyBADoJsIQYuHkS3STAxAg3yWVoDIEAOgGwhBikakM8eMHIEBeQh6VIQBAN3A1ilg4+fL58QMQoMw+Q37cwwAAlBCuRhELTz7d5AAEylxSSTZdBQB0A1ejiIWTMU0OQLC8BA0UAADdwtUoYuGZL+PHD0CAzEsqYUyTAwAUjqtRxMIxTQ5AwMxRGQIAdA9Xo4gFa4YABM5jnyEAQPdwNYromcmTyScMAQiSSypJNzkAQDdwNYro5eb0E4YABMi8hDxnMp/qEACgMFyNInr5MJSIdxwAykuiSpKU7tgU80AAAKWCMIToZX9ra87FPBAAZcVLSpLSacIQAKAwhCFEL1sZMipDAILkZf5N8Tex8SoAoDCEIUTPMpUhx5ohAEFymcpQR3pjzAMBAJQKrkYRvXxliB8/AMGxRCYMWQeVIQBAYbgaRfTya4b48QMQHJdfM0QYAgAUhqtRRM8sc8uaIQBByochpskBAApDGEL0smuG2GcIQKC83DQ59hkCABSGq1FEj32GAIQgN03Op7U2AKBAybgHgG1Y/5HUsSHuUQRv7YrMLZUhAEHKNlDwaaAAACgQYahYrXhD+mn9J1PKylDa6xH3EACUES8XhnwqQwCAwhCGitWapZkgNPqb0i77xj2awP3gj29rY/8xcQ8DQDnJTZPrIAwBAApDGCpW2fbTGnKy9NlR8Y4lBDP+9ykdU9Un7mEAKCPOq5LEmiEAQOFYtFGsck0GvPJsMmBm8lzcowBQTlxumhz7DAEACkQYKlZ+9j/zMm0y4JvJc6QhAMHJhSERhgAABSrPK+1ykJsm55XnTEbfRBgCECiXyFTS04QhAECBCEPFKtdFrkynyfk+lSEAwXKJzJohY80QAKBAhKFilZ8mV6ZhiDVDAALmZSvp5lMZAgAUhjBUrPxcA4UyniZHGgIQJDZdBQB0E2GoWOWnyZXnH1GaBgoAApbITpMTlSEAQIHK80q7HJT5NDlaawMImkvmWmuzZggAUBjCULGimxwAdEtu01VaawMACkUYKlZl3k0u7RtrhgAEKjdNjgYKAIBCEYaKVa4yVIbT5MxMkpgmByBQuWlyRmUIAFAgwlCx8su3MuRnshDT5AAEigYKAIDuIgwVqzKeJpf2qQwBCJ6XZJ8hAED3EIaKVRl3k/Nz0+RIQwAC5HINZwhDAIACEYaKVRl3kzOmyQEIQTLZI/MJa4YAAAUiDBWrMp4m59NAAUAIXCLz76UZYQgAUBjCULEq425y6XwYIg0BCE4imdtniE1XAQCFIQwVKz8tyUle+f0RmZ+5JQwBCFIiWxnK/zIJAIDtKL8r7XJh6bKcIicxTQ5AOBKep42WkKOBAgCgQOW3Or9c+B1lOUVO+mSaXII0BCBACc8prYT2XDFL+v0VcQ8nHDWnSYMPj3sUAFA2CEPFyk+XZSc56ZPKkGOaHIAAJTynVn+I6jcskub9Lu7hBK/9A+njZYQhAAhQeV5tlwPzy3aaHK21AYQh4Tmdv+k7uvroL+iisZ+PezjB+9l4qWND3KMAgLJCGCpS6Y5N6khLP3/6zbiHErg16zPz+ROsWAMQoNwvWDp8i3kkIUn0kNIb4x4FAJQVwlCRWrmmXYlN0g3/+5e4hxKKpOe0z4DecQ8DQBlJZtch+mUbhqpoGw4AASMMFSk/3SEnTw98Y7TqBg+IeziBc5KSlIYABCjXlCXXpKXsJHpIm1bHPQoAKCuEoWLldygtTz2rPFURGgBgu5xz8pyULtfKULKnlGbNEAAEiTBUpMxPy5eXn/YBANi+hOfKNwwxTQ4AAkcYKlZ+Wmnz2IsHALrBc06vvrda9/15UdxDCVzDB5u098YNqop7IABQRghDxSo7TS7pMUUOAAq1587Vev7NlXr+zZVxDyVw1yc/Vv9e67RL3AMBgDJCGCpSZr7SSlAZAoBu+MO/HKmP2stzKlnbrdPlMU0OAAJFGCpW2cpQL8IQABSsV1VCvarKdMPqRA8lOghDABAk5mAVq2wDhagqQ6lUSs45pVKpSM4HAOge36tS0jriHgYAlJVYwpBzbqFz7lXn3BznXFv2vl2cc086597M3pbf5jrd4afVIU/JRHRhqPMtAKC4mNdDSaMyBABBirMyNN7Mas2sLvv1VZKazWx/Sc3ZryuXRVsZAgB0X5RVdd+rUpU6pHLdVBYAYlBM0+ROljQ9+/l0SafEOJbYOT+tDiXoJgcARSzKqrp52abaNFEAgMDEdaVtkv7XOTfbOXdh9r7dzWypJGVvd4tpbMXB0kpTGQIAZFmiR+aT9MZ4BwIAZSSubnKHm9kS59xukp50zi0o9InZ8HShJH32s58Na3zxs7R885QkDAEA1LkyRBgCgKDEUhkysyXZ2+WSHpE0UtIy59wekpS9Xb6V5/7MzOrMrG7QoEFRDTlyLttAgcoQAECSlK8MMU0OAIISeRhyzvVxzvXLfS7pGEmvSXpU0sTsYRMl/S7qsRUTl22gQGUIACBJlshVhjbEOxAAKCNxTJPbXdIjzrnc+e83syeccy9K+rVz7uuS/ibp9BjGVjwsrbSSVIYAABlUhgAgcJGHITN7R9LwLu5fJakx6vEUK+enlXYJZUNj6FpbW/O39fX1kZwTANANNFAAgMDRt7lIOfPlW3RBqLExk0MbGxvzwShsbPAKAN1AGAKAwMXVTQ7b4SytDt8iqwzltLe3q6GhIbLzEYgAlLIoq+ouyTQ5AAgalaEi5Swt86pkZqF/tLS0qLq6WpJUXV2tlpaWSM7b1NQU83cZQDmK6pcsUVfVXbYy9Iuf3xHqeQCgkjgzi3sMO6yurs7a2triHkYoVv1wiP60cT+dkPp9JOdrbW1VQ0ODWlpaWDMEoKRFXVGPyrHHNeqJUS/q6Ps26qk32+MeDgCUDOfcbDOr6+oxKkNFysmXuURk58sFIIIQgFLX1NRUllX1S755mSTp7DO/Gue3FwDKCmGoSHl+Wn6EYQgAykVU0+Tq6+vVf4yxsAAAFKBJREFU3NwsSWpubg79l0leds3QaRNODvU8AFBJCENFyiktc/zxAEAxi7SqnuwpSUpvYtNVAAgKV9tFyjNfojIEAP9/e/cfJVdd33/8+Z6Z3U02CUL4kUIIipEaKEcRI5qgQonyta0V9QgltmitHotWivVHwa9a1q8/Dkj1K361pRb50i9FPVSxamsFTqRidYFEUEApIJSaQCTByI8km+zu7Of7x9xJhs2v3ZyZe+/sPB/n5MzM3Zn5vGfvye597ed9P1eZap9hSJLazTBUUpVkm5wkaadq1iY3Me51hiSpXQxDJWWbnCSpVaXPMCRJ7ebRdklVUr6ryUmSyq3ZJjcxtq3gSiRp5jAMlVSFumFIkrRDLVtAwZkhSWofw1BJ5b2AQnMp2ryWpJUkTU+zTS4ZhiSpbWpFF6DdqzABlXzDkEFIksqrlrXJGYYkqX2cGSqjlKjgOUOSVHZ5zqr31apsTzVS3TAkSe3izFAZTdQbtznODEmSpi/PWfW+ajBGzZkhSWojZ4bKKDXCkDNDkqSmvmqFMWow7kVXJaldDENlNDHeuHVmSJKU2TEzNDFWdCmSNGMYhsqo2SbnzJAkKdNXrTBKDWyTk6S2MQyVUfKcIUnS0/VVK4ylKjFhGJKkdjEMlZELKEiSJqllbXLUbZOTpHYxDJWRbXKSpEn6qxVG6SMMQ5LUNoahMtrRJufK55KkhuZqcrbJSVL7GIbKKFtNLmyTkyRlatVglBrhRVclqW0MQ2XkOUOSpEn6Ko0FFCourS1JbWMYKqM0ATgzJEnaqVIJxqOPMAxJUtsYhspox0VXPWdIkrTTeNScGZKkNvJou4xcTU6StBv1qDF/2y/gi6cXXUpnHHYc/P5niq5CUg8xDJVRtppcVA1DkqSdbqy8jMWD4zynb3bRpbTfrx+CtbfC730aKjauSMqHYaiMdqwm5+6RJO30/dpLqD7rNVzyhucVXUr7/eCzcOOHYWwrDMwtuhpJPcI/vZRQytrkDEOSpFZ91Qpj9Ymiy+iM/sHG7djWYuuQ1FMMQyU0UbdNTpK0q/5qhbGJVHQZndE3p3E7uqXYOiT1FMNQCdXrjZWCnBmSJLWqVYOxcWeGJKldDEMlNFFvnjPkzJAkaacZ3SbnzJCkAhiGSqjebJMzDEmSWvTN5Da55syQYUhSjgxDJTQx3miTq1Rtk5Mk7dQ3k9vk+myTk5Q/w1AJ7VxAwTAkSdppRrfJ9Tfb5AxDkvJjGCqh5jlDFVeTkyS1mNFtcjtmhmyTk5Qfw1AJ7VhNzpkhSVKLvNvkhoaGiAiGhoY6P9iOc4acGZKUH4+2S6h50dWqS2tLklr0VSs8MTLGTf+5IZfxTjn7ncy6+l855ex3dnzMqI9yKjgzJClXHm2XUL25tLZtcpKkFgfN6efhx0d4y1WrcxtzwZlDuY13/0CV2uhWIpfRJMkwVEppxzlDfQVXIkkqk7969XGctXRRrmO++MUv5tZbb+34ONfdvo6R2wcY3L7ZgxNJufHnTQk1V5Or2CYnSWoxq6/KCYsOzHXM0fX35TLmT9Y+zlYG6N+2xYMTSblxAYUS2nHOkG1ykqQeMWegxtY0wMS2zUWXIqmHGIZKaKK5mlzNv41JkoozPDz8tNtOmtNfZYQBJkZdQEFSfgxDJZTqzZkhw5Ak6elyWeaaRgBasWIFACtWrOh4IBocqLGVASa2G4Yk5SdS6t6Lty09opbW/OkBRZfRdilNECRuO2sNJx13TNHlSJJKJGJmrrU2sHAJ17/lEF5wWHDAeTcXXY6kGSQifpRSWrq7r3X31MO8BfDSc4uuou3W/Xorl9+xjTNmH1x0KZKkkrnoootymR1qzgyNjIwwe/ZsVq1axbJlyzo23j3rn+Shv3k9MfZ4x8aQpMm6PAwdDis+XHQVbffgfRu5Zs1tvL4yM//6J0naf3m1yS1btoxVq1axfPnyjgchgLlZm1xlbGtHx5GkVp4zVEL1iQkAaoYhSVKBmgGo00EIYLC/ykgaoFI3DEnKj2GohMbrjfO4qoYhSVKPmDNQYyuzqI2PFF2KpB5iGCqh+kQjDNWqhiFJUm8YqFXYFgP0TWyDrENCkjrNMFRC480w5MyQJKlHRATj1dmNB84OScqJYaiEmjND1Yq7R5JUnOZiDXkt2lCvDjbueOFVSTnp7tXkZihnhiRJZTA0NJRbEAJIfbOhDnzrfOgbzG3cXMxdAKd/FCrVoiuR1MIwVELN1eRcQEGS1Ev+e9YS1taPZtHGe4supb1GN8PmR2HpW+AQL6YulYlhqIScGZIk9aJNg4t576y/5dpzO7+Ud64e+C5c/TrYvMEwJJVMV4ehe3/5FKdeelMuY/1q0yY2/WoT8w+ez8Hz53d0rCe3jQPODEmSesucgSqPPL6t6DLab85hjdstG4qtQ9IuujoMDfZXef6iA/MZbNGBfOmaWzjttDfmMtxvPGMW8+f05zKWJEllMNhfY+voeNFltN/cLAxt3lhsHZJ20dVhaNH8QS47+wW5jffZlSdy2bcuzW08SZJ6yZyBGpu314suo/0GD4aoODMklZBrN0/R8PDw024lSVJ7zemvzsyZoUq1EYg2G4aksunqMPTII4/kMs7w8DArVqwAYMWKFbkFojyXM5UkqWiDAzW2jtaZyBYSmlHmLoAttslJZdPVbXLr168nIt9FBkZGRli+fHlu4xmIJEm9Yu5A4xo8W8fqzB3o6kOUXc051JkhqYS6+ifN4YcfnsvsUHNmaGRkhNmzZ7Nq1SqWLev8sp8GIUlSLxnsbxyW/PKJbRw6b6CjY1188cVccsklXHDBBVx44YUdHQtg7uChVDc90PFxJE1PpNS9U9FLly5Na9asyWWs4eFhli9fzg9/+MNcgpAkSb3mWz95hPO+fEfRZXTE0MCXeHPfKuKD6yHnrhap10XEj1JKS3f3ta6eGcpTMwAZhCRJ6owVxx7Gx157PNvHJzo+1gN3386nPvZh3vuhj7L4+BM7OtZDj23hkdXziBiB0c0wMK+j40maOsOQJEkqhU9+4mO5tIgPDw9z3nvOYXRkhM+855yOt78/sHEzn7/tGY0HmzcYhqQSsU1uioaGhvjIRz7CRRdd5Lk8kiR1QN6LIuWm2scfvf8irh74JBxzemOZ7ZmmUoOXvRfmH110JdIu9tYmZxiSJEmlMDQ0lNvMUN4LI/3uJ77G/9n2IRbPn6FNOY+vhZe/D077UNGVSLswDEmSJLXIe2Gks/5umJQS/3RufpfnyNXlL20sH37O14uuRNrF3sJQV190VZIkaX/kvTDSM+cP8otNW3MZqxALXwgP/wgmOr/4hdROhiFJktRzmu14eZ0HfNT8QR59cjvbxuq5jDc0NERE5Hee88IXwrYnYNOD+YwntckMbVyVJEnas7zOT2o66uBBAL7x44dZcMCsjo936JIXUTvoCA5d8iL+/d4NHR/vgPpiTgR44LvQP9jx8XIVVZi3oOgq1CGlO2coIl4FXAZUgStSShfv6bmeMyRJkrrBPeuf5Hcu+34uY21/+B4e/coHSfVxolpjwdkfZ2DhsR0ds8IEdw28jTmxraPjFObkd8MrP1J0FdpPXbOAQkRUgfuAVwLrgNXAypTSz3b3fMOQJEnqFv2HHEUMdH7WpP7Ur6g/9diOx9V5h1Cd19nlvA88+Y0sOzK49g8W0l+bYWdh3H8j3PcdeOctcPBziq6mc2bq0vbsPQyVrU3uJODnKaUHASLiK8AZwG7DkCRJUrf4n+/6k2KWDr/+mx1fKGL1Q5s48/JhzlpzIAcN9nV0rLwdMHEEl8ZN9H9ut8fSM0J9/mIqp36AGJxfdCnt1z93r18uWxhaCKxtebwOeHFBtUiSJLVNXucoLVu2jJUrV3LllVeycuXKXFbMW/rMg1h50lH89JEn+NWW0Y6Pl6e1owOcvfX9vKxyV9GldEQlEr//2A959nVvK7qUQpStTe5M4H+klN6WPT4HOCmldF7Lc94OvD17+Fzg3hxLPAR4bJ/PUlm5/7qX+667uf+6m/uve7nvupv7r32emVI6dHdfKNvM0DpgUcvjI4FHWp+QUvoC8IU8i2qKiDV76jdU+bn/upf7rru5/7qb+697ue+6m/svH2U7w201cExEHB0R/cDZwDcLrkmSJEnSDFSqmaGU0nhEvAu4nsbS2lemlH5acFmSJEmSZqBShSGAlNK3gW8XXcceFNKep7Zx/3Uv9113c/91N/df93LfdTf3Xw5KtYCCJEmSJOWlbOcMSZIkSVIuDENTFBGvioh7I+LnEXFh0fVoaiJiUUTcFBH3RMRPI+L8omvS9EVENSLuiIh/KboWTV1EHBgRX42I/8z+D3b+Yidqm4j4i+zn5t0R8eWImFV0TdqziLgyIjZExN0t2+ZHxI0RcX92e1CRNWr39rDvLs1+dt4ZEV+PiAOLrHEmMwxNQURUgc8DvwMcB6yMiOOKrUpTNA68N6V0LPAS4M/cd13pfOCeoovQtF0GfCeltAR4Pu7DrhERC4E/B5amlI6nsajR2cVWpX24CnjVpG0XAqtSSscAq7LHKp+r2HXf3Qgcn1J6HnAf8IG8i+oVhqGpOQn4eUrpwZTSKPAV4IyCa9IUpJTWp5Ruz+4/ReNgbGGxVWk6IuJI4PeAK4quRVMXEQcALwe+CJBSGk0pPV5sVZqmGjA7ImrAIJOu+6dySSndDGyatPkM4B+y+/8AvDbXojQlu9t3KaUbUkrj2cNbaFx7Ux1gGJqahcDalsfr8IC660TEs4AXALcWW4mm6TPAXwITRReiaXk2sBH4v1mL4xURMafoojQ1KaWHgb8GfgGsB55IKd1QbFXaDwtSSuuh8cdB4LCC69H++RPg34ouYqYyDE1N7Gaby/B1kYiYC3wNeHdK6cmi69HURMSrgQ0ppR8VXYumrQacCPxtSukFwBZs0eka2bklZwBHA0cAcyLij4qtSuo9EfFBGi3/1xRdy0xlGJqadcCilsdHYrtA14iIPhpB6JqU0nVF16NpORl4TUQ8RKM99bSI+MdiS9IUrQPWpZSaM7FfpRGO1B1eAfxXSmljSmkMuA5YXnBNmr5HI+JwgOx2Q8H1aBoi4s3Aq4E/TF4Lp2MMQ1OzGjgmIo6OiH4aJ5F+s+CaNAURETTOWbgnpfTpouvR9KSUPpBSOjKl9Cwa/+++m1Lyr9NdIKX0S2BtRDw327QC+FmBJWl6fgG8JCIGs5+jK3ABjG70TeDN2f03A98osBZNQ0S8CrgAeE1KaWvR9cxkhqEpyE5gexdwPY1fBtemlH5abFWaopOBc2jMKPw4+/e7RRcl9YjzgGsi4k7gBOATBdejKcpm9L4K3A7cReN44QuFFqW9iogvA8PAcyNiXUS8FbgYeGVE3A+8MnusktnDvvscMA+4MTt2ubzQImewcNZNkiRJUi9yZkiSJElSTzIMSZIkSepJhiFJkiRJPckwJEmSJKknGYYkSZIk9STDkCT1uIhIEfGplsfvi4ihNr33VRHxhna81z7GOTMi7omImzo9liRp5jAMSZK2A6+PiEOKLqRVRFSn8fS3Au9MKf12p+qRJM08hiFJ0jiNC2r+xeQvTJ7ZiYjN2e2pEfG9iLg2Iu6LiIsj4g8j4raIuCsiFre8zSsi4vvZ816dvb4aEZdGxOqIuDMi/rTlfW+KiC/RuNjn5HpWZu9/d0Rckm37K+ClwOURcemk5x8eETdnFy28OyJe1vo5svtviIirsvsLIuLrEfGT7N/ybPubsjp/EhFXZ9sOjYivZZ9hdUScnG0/peUiz3dExLy91HF6RAxHxO0R8U8RMTfbfnFE/Cwb86+nsS8lSdNQK7oASVIpfB64MyI+OY3XPB84FtgEPAhckVI6KSLOB84D3p0971nAKcBi4KaIeA7wJuCJlNKLImIA+EFE3JA9/yTg+JTSf7UOFhFHAJcALwR+DdwQEa9NKf2viDgNeF9Kac2kGt8IXJ9S+ng20zS4j8/0WeB7KaXXZc+fGxG/BXwQODml9FhEzM+eexnwv1NK/xERRwHXZ9+P9wF/llL6QRZutgFvn1xHNhP3IeAVKaUtEXEB8J6I+BzwOmBJSilFxIH7qFmStJ8MQ5IkUkpPRsT/A/4cGJniy1anlNYDRMQDQDPM3AW0tqtdm1KaAO6PiAeBJcDpwPNaZp2eARwDjAK3TQ5CmRcB/55S2piNeQ3wcuCf91YjcGVE9AH/nFL68T4+02k0ghoppTrwRES8CfhqSumxbPum7LmvAI6LiOZrD4iIecAPgE9n9V2XUloXEbvUERGnAMfRCIIA/cAw8CSNAHVFRPwr8C/7qFmStJ9sk5MkNX2Gxrk3c1q2jZP9rojGEXt/y9e2t9yfaHk8wdP/2JYmjZOAAM5LKZ2Q/Ts6pdQMU1v2UF/sYfsepZRuphGYHgauzoLN5Jpm7eNtgl0/AzS+L8taPsPClNJTKaWLgbcBs4FbImLJHuoI4MaW1x+XUnprSmmcxuzY14DXAt+Z7ueWJE2NYUiSBOyY8biWRiBqeohGWxrAGUDffrz1mRFRyc4jejZwL42WsndkMyVExG9GxJy9vQlwK3BKRByStZqtBL63txdExDOBDSmlvwe+CJyYfenRiDg2Iio0WtKaVgHvyF5bjYgDsm1nRcTB2fZmm9wNwLtaxjohu12cUrorpXQJsAZYsoc6bgFOztoGiYjB7PswF3hGSunbNFoNT9jH90WStJ9sk5MktfoULQf4wN8D34iI22iEgj3N2uzNvTRCywLg3JTStoi4gsa5RLdnM04bacyC7FFKaX1EfAC4icasyrdTSt/Yx9inAu+PiDFgM1kLHHAhjfaztcDdwNxs+/nAFyLirUAdeEdKaTgiPg58LyLqwB3AH9NoKfx8RNxJ4/fpzcC5wLsj4rez1/8M+Dfg7Ml1pJQ2RsQfA1/OzpuCxjlET9H4ns/KPucuC1tIktojUtrdzL8kSZIkzWy2yUmSJEnqSYYhSZIkST3JMCRJkiSpJxmGJEmSJPUkw5AkSZKknmQYkiRJktSTDEOSJEmSepJhSJIkSVJP+v8IjjsbVtwZmwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "yaxis = func_poisson(np.floor(xaxis+0.5), 1000, 4.0)\n", "ax.plot(xaxis, yaxis, '-', label='Poisson distribution')\n", "ax.legend()\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And save the figure:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "if save_plots: \n", " fig.savefig(\"BinomialPoisson.pdf\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculation of Binomial $\\chi^2$-value:\n", "\n", "In this part of the exercise, you are asked to calculate the ChiSquare value yourself, in order to ensure that you understand exactly what is going on!\n", "\n", "Above, we have (using Minuit) *fitted* the distribution, but as we know the initial values, I would like you to calculate the $\\chi^2$-value between the data and the binomial they were generated from, i.e. with NO free parameters.\n", "\n", "I suggest you use Pearson's $\\chi^2$, and require `N_obs` and/or `N_exp` > e.g. 0.1. Remember that your choice should ensure, that there is no division by zero (which is a cardinal sin in programming)!" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "N_bins = len(x) # Just to know how many bins to loop over\n", "chi2_bino = 0.0 # This you'll add to\n", "N_dof = 0 " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "for N_obs, x_i in zip(y, x):\n", " N_exp = func_binomial(x_i, 1000, 20, 0.2)\n", " if (N_obs > 0) :\n", " chi2_bino += 0.0 # Write the expression yourself!\n", "\n", "# Also calculate Ndof and Prob:\n", "Ndof_bino = 1 # Think about the number of degrees of freedom when there are no free parameters!\n", "Prob_bino = stats.chi2.sf(chi2_bino, Ndof_bino)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Just a test printout, change to match your own results! \n", "\n", "Binomial: chi2 = 9999.00 N_dof = 9999 Prob = 9999.0000\n" ] } ], "source": [ "print(f\"Just a test printout, change to match your own results! \\n\")\n", "print(f\"Binomial: chi2 = {9999:.2f} N_dof = {9999:d} Prob = {9999:6.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "\n", "# Questions:\n", "\n", "Important: Make sure you understand what process yields a Binomial, a Poisson and a Gaussian distribution. Without this knowledge, this exercise and a large fraction of the course will be lost on you!\n", "\n", "1. Plot a Binomial ($N=20$, $p=0.2$), Poisson ($\\lambda = 4$), and Gaussian ($\\mu=4$, $\\sigma=\\sqrt{4}$), i.e. same means and widths. Which distribution has the longest tail towards high values? And which one has the longest tail the other way? Does this pattern depend on the parameters (given same means and widths)? Play around with the settings (remember also to change the scale (use log) of the plot accordingly), and gain your own experience. And most importantly perhaps, in what limits do they start looking like each other?\n", "\n", "2. Producing binomially distributed numbers (using the parameters `N_experiments=1000`, `N_trials=10` and `p_success=0.2`), and plot the Binomial and Poisson distributions on top. Do they match? Quantify this by calculating the $\\chi^2$ fits of the resulting distribution with a Binomial and Poisson distribution. Do you get acceptable p-values with all of these? If not, investigate for what choice of parameters you do.\n", "\n", "3. In all of the above $\\chi^2$ fits, we have _assumed_ that the uncertainty on the count in each bin is Gaussianly distributed! Ask yourself to what extend this requirement is fulfilled? Does changing the parameters (`N_experiments`, `N_trials` and `p_success`) \"help\" fulfilling this requirement, and if so, which and how?" ] } ], "metadata": { "executable": "/usr/bin/env python", "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.6" }, "main_language": "python" }, "nbformat": 4, "nbformat_minor": 2 }