{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hypothesis Testing\n", "\n", "Python notebook for illustrating the concept of Hypothesis Testing and specific test statistics; among them the very useful Kolmogorov-Smirnov test.\n", "\n", "The Kolmogorov-Smirnov test (KS-test) is a general test to evaluate if two distributions in 1D are the same. This program applies an unbinned KS test, and compares it to a $\\chi^2$-test and a simple comparison of means. The distributions compared are two unit Gaussians, where one is then modified by changing:\n", "- Mean\n", "- Width\n", "- Normalisation\n", "\n", "The sensitivity of each test is then considered for each of these changes.\n", "\n", "### References:\n", "- Barlow: p. 155-156\n", "- __[Wikipedia: Kolmogorov-Smirnov test](http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test)__\n", "- Though influenced by biostatistics, a good discussion of p-values and their distribution can be found here:\n", " [How to interpret a p-value histogram?](http://varianceexplained.org/statistics/interpreting-pvalue-histogram/)\n", "\n", "### Authors: \n", "Troels C. Petersen (Niels Bohr Institute)\n", "\n", "### Date: \n", "03-12-2020 (latest update)\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.pyplot as plt # Plots and figures like you know them from Matlab\n", "import seaborn as sns # Make the plots nicer to look at\n", "from iminuit import Minuit # The actual fitting tool, better than scipy's\n", "import sys # Module to see files and folders in directories\n", "from scipy.special import erfc\n", "from scipy import stats" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the parameters of the plot:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "r = np.random # Random generator\n", "r.seed(42) # Set a random seed (but a fixed one)\n", "\n", "save_plots = False\n", "verbose = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and define a function that calculates the mean, standard deviation and the standard deviation (i.e. uncertainty) on mean (sdom):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def mean_std_sdom(x):\n", " std = np.std(x, ddof=1)\n", " return np.mean(x), std, std / np.sqrt(len(x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up the experiment:\n", "\n", "How many experiments, and how many events in each:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "N_exp = 1000\n", "N_events_A = 100\n", "N_events_B = 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the two Gaussians to be generated (no difference to begin with!):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "dist_mean_A = 0.0\n", "dist_width_A = 1.0\n", "dist_mean_B = 0.0\n", "dist_width_B = 1.0" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 2 }, "source": [ "Define the number of bins and the range, initialize empty arrays to store the results in and make an empty figure (to be filled in later):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "N_bins = 100\n", "xmin, xmax = -5.0, 5.0\n", "\n", "all_p_mean = np.zeros(N_exp)\n", "all_p_chi2 = np.zeros(N_exp)\n", "all_p_ks = np.zeros(N_exp)\n", "\n", "# Figure for the two distributions, A and B, in the first experiment:\n", "fig1, ax1 = plt.subplots(figsize=(10, 6))\n", "plt.close(fig1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Loop over how many times we want to run the experiment, and for each calculate the p-value of the two distributions coming from the same underlying PDF (put in calculations yourself):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.70206\n", " 1: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.70206\n", " 2: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.70206\n", " 3: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.70206\n", " 4: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.36819\n", " 5: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.28194\n", " 6: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.90841\n", " 7: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.03638\n", " 8: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.90841\n", " 9: p_mean: 0.50000 p_chi2: 0.50000 p_ks: 0.28194\n", "Got to experiment number: 1000\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAf8ElEQVR4nO3de7SlZ10n+O+PqtAhBASS2EoqnhNUUIaWkC5RLuOFiw0aQNoLQRRwHOkeWkEbFwO009Jrlt3jpRG68QbiDCIkjRHSjSi3haAoCJUYx0AiMKSKVBKwEoghmpik+M0fZxd5qDrn1D5Ve5996tTns9Zetd/3vPt5fs9+zz7ru5569vtWdwcAAFhxj0UXAAAAW4mADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJCBbauqPlpV37HoOraiqvqnVfUnVfWFqvrPi65nVFXLVdVVtXPRtQAnJwEZOCFV1d6qesJh+55bVR84tN3d/1N3v+8o7ZysYex5SW5Mct/uftFaB1XVyyfvzyM3r7T1Tc79bVV1a1V9vqreXlXnLLouYPsQkAHmaAsH76UkH+t17hZVVZXkR5J8LslzNquwKT2lu09P8tVJPpvkvy64HmAbEZCBbWucZa6qR1bVnqq6pao+W1WvmBz2J5N/b57MSD6qqu5RVT9bVfuq6m+r6neq6iuGdp89+dlNVfV/HNbPy6vqkqr63aq6JclzJ31/sKpurqobqurVVXXPob2uqudX1ScmSx7+z6r62slrbqmqNx86vqrOrKo/mLT1uar606pa9W95VT26qj5SVX83+ffRk/3/T1YC74snY37Caq9P8j8neWCSFya5cKx5lb6mGeO/nozx81X1q5MAnqraUVW/XFU3VtWnknzPWv0crrtvT3JJkodO+xqAoxGQgZPFq5K8qrvvm+Rrk7x5sv/bJv/er7tP7+4PJnnu5PGdSR6U5PQkr06Sqnpokl9L8qyszF5+RZKzD+vraVkJbfdL8sYkB5P8dJIzkzwqyeOTPP+w1zwpyT9P8q1JXpzkNZM+zknysCTPnBz3oiT7k5yV5J8meVmSI2aBq+oBSd6e5L8kOSPJK5K8varO6O7nTur6xcmY37PGe/acJG9L8t8m2xescVymHOMFSb45ycOT/GCSfzHZ/+OTnz0iye4k379OP1+mqk5L8owkH5r2NQBHIyADJ7JLJzOWN1fVzVkJrmu5M8nXVdWZ3X1rd68XqJ6V5BXd/anuvjXJS7Myg7ozK+Htbd39ge6+I8m/z5EB9YPdfWl3f7G7b+vuy7r7Q919V3fvTfKbSb79sNf8Qnff0t0fTXJlkndN+v+7JH+UlfB4aBxfnWSpu+/s7j9dY5nE9yT5RHe/YdLvRUmuTvKUdcb9JZPg+QNJ3tTdd2Yl8K+5zGLKMf5f3X1zd386yR8nOW+y/weTvLK7r+3uzyX5T1OUeOnknN+S5IlJfmmacQFMQ0AGTmTf2933O/TIkTOWox9L8uAkV0+WG6w3G/rAJPuG7X1JdmZlxvaBSa499IPu/ockNx32+mvHjap68GRZxGcmyy7+Y1ZmWkefHZ7ftsr26ZPnv5Tkk0neVVWfqqqXTDmGQ+M4fLZ7LU9PcleSP5xsvzHJk6vqrNUOnnKMnxme/0PuHtOXvaer1L2a752c83+S5CeSvL+qvmqK1wEclYAMnBS6+xPd/cwkX5nkF5JcUlX3zirLE5Jcn5UvsR3yNVkJi59NckOSXYd+UFX3ysoShi/r7rDtX8/K7O3XT5Z4vCxJHeM4vtDdL+ruB2VlNvjfVtXjpxjDoXFcN2VXz8lKgP10VX0mye8lOSV3L/U43PGM8YasLCUZ65xKdx/s7rdkZYnHY6d9HcB6BGTgpFBVP1xVZ3X3F5PcPNl9MMmBJF/MylrjQy5K8tNVdW5VnZ6V2dD/1t13ZWWpwVMmX4C7Z5L/kKMHwftkZSnArVX1DUn+t+MYxwVV9XWTL7jdMhnDwVUO/cMkD66qH6qqnVX1jKx8ke0Ppujj7KysIb4gK8sgzsvKuuFfyNrLLI5njG9O8oKq2lVV90+y1qz4arVWVT0tyf2TXLWBPgHWJCADJ4snJfloVd2alS/sXdjdt0+WSPx8kj+brGX+1iS/neQNWbnCxTVJbk/yk0kyWSP8k0kuzsrM5xeS/G2Sf1yn759J8kOTY1+bu7/0diy+Psl7ktya5INJfm21az13901ZCbgvysoSkBcnuaC7b5yijx9JckV3v6u7P3PokZUv/H1TVT1sldcczxhfm+SdSf4qyeVJ3jLFa942OZe3ZOX8PWdybgCOW61zCUwAjmIyw3xzVpYWXLPoegA4fmaQATaoqp5SVadN1jD/cpK/TrJ3sVUBMCsCMsDGPS0rX4K7PitLHi5c7450AJxYLLEAAICBGWQAABjsXHQBozPPPLOXl5cXXQYAACeByy677MbuPuIGSFsqIC8vL2fPnj2LLgMAgJNAVa16505LLAAAYCAgAwDAQEAGAIDBllqDDADAiePOO+/M/v37c/vtty+6lHWdeuqp2bVrV0455ZSpjheQAQA4Jvv378997nOfLC8vp6oWXc6qujs33XRT9u/fn3PPPXeq11hiAQDAMbn99ttzxhlnbNlwnCRVlTPOOGNDs9wCMgAAx2wrh+NDNlqjgAwAAAMBGQCAmVheTqpm95j2BstvfetbU1W5+uqrZzIOARkAgJnYty/pnt1j36r3uTvSRRddlMc+9rG5+OKLZzIOARkAgBPWrbfemj/7sz/L6173OgEZAAAuvfTSPOlJT8qDH/zgPOABD8jll19+3G0KyAAAnLAuuuiiXHjhhUmSCy+8MBdddNFxt+lGIQAAnJBuuummvPe9782VV16ZqsrBgwdTVfnFX/zF47r8nBlkAABOSJdcckme/exnZ9++fdm7d2+uvfbanHvuufnABz5wXO0KyAAAzMTS0mwv87a0tH5/F110UZ7+9Kd/2b7v+77vy5ve9KbjGoclFgAnm+XlVa+dtH/HUs45uPeI/UtLyd4jdwMcYbP/Vrzvfe87Yt8LXvCC425XQAY42Ry6UOlhdlWttjsnwF1kAWbKEgsAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADADAbCwvz/ZCyMvLR+1yx44dOe+88/Lwhz88559/fv78z//8uIfhMm8AAMzGGpeRPGZTXGfyXve6V6644ookyTvf+c689KUvzfvf//7j6tYMMgAA28Itt9yS+9///sfdjhlkAABOWLfddlvOO++83H777bnhhhvy3ve+97jbFJABADhhjUssPvjBD+bZz352rrzyytRx3AbUEgsAALaFRz3qUbnxxhtz4MCB42pHQAYAYFu4+uqrc/DgwZxxxhnH1Y4lFgAAzMbS0lRXnthQe0dxaA1yknR3Xv/612fHjh3H1a2ADADAbOzdu+ldHjx4cOZtWmIBAAADARkAAAYCMgAAx6xneee8OdlojQIyAADH5NRTT81NN920pUNyd+emm27KqaeeOvVrfEkPAIBjsmvXruzfv/+4rzs8b6eeemp27do19fECMgAAx+SUU07Jueeeu+gyZs4SCwAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAIDBXANyVf10VX20qq6sqouq6tR59gcAAMdrbgG5qs5O8oIku7v7YUl2JLlwXv0BAMAszHuJxc4k96qqnUlOS3L9nPsDAIDjMreA3N3XJfnlJJ9OckOSv+vudx1+XFU9r6r2VNWeAwcOzKscAACYyjyXWNw/ydOSnJvkgUnuXVU/fPhx3f2a7t7d3bvPOuuseZUDAABTmecSiyckuaa7D3T3nUnekuTRc+wPAACO2zwD8qeTfGtVnVZVleTxSa6aY38AAHDc5rkG+S+SXJLk8iR/PenrNfPqDwAAZmHnPBvv7p9L8nPz7AMAAGbJnfQAAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAICBgAwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAyeH5eWk6sjH8rKSAPgyOxddAMCm2Lcv6T5yf9Xm1zKxBUsCIGaQAQDgywjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAICBgAwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYzDUgV9X9quqSqrq6qq6qqkfNsz8AADheO+fc/quSvKO7v7+q7pnktDn3BwAAx2VuAbmq7pvk25I8N0m6+44kd8yrPwAAmIV5LrF4UJIDSf7vqvrLqvqtqrr34QdV1fOqak9V7Tlw4MAcywHgWFyT5aRq9cfy8qLLA5i5eQbknUnOT/Lr3f2IJH+f5CWHH9Tdr+nu3d29+6yzzppjOQAci+XsS7pXf+zbt+jyAGZungF5f5L93f0Xk+1LshKYAQBgy5pbQO7uzyS5tqoeMtn1+CQfm1d/AAAwC/O+isVPJnnj5AoWn0ryo3PuDwAAjstcA3J3X5Fk9zz7AACAWXInPQAAGAjIAAAwEJABAGAw9RrkqnpskkcmubK73zW/kgAAYHHWnEGuqg8Pz388yauT3CfJz1XVETf8AACA7WC9JRanDM+fl+SJ3f0fknxXkmfNtSoAAFiQ9ZZY3KOq7p+VEF3dfSBJuvvvq+quTakOAAA22XoB+SuSXJakknRVfVV3f6aqTp/sAwCAbWfNgNzdy2v86ItJnj6XagAAYMHWvYpFVVVWrlxxdpJOcn2SD3f3NZtQGwAAbLo1A3JVfVeSX0vyiSTXTXbvSvJ1VfV8l3oDAGA7Wm8G+VVJntDde8edVXVukj9M8o1zrAsAABZivcu87Uyyf5X91+XLLwEHAADbxnozyL+d5CNVdXGSayf7zklyYZLXzbswAABYhPWuYvGfqurSJE9L8qisXNptf5JndffHNqk+AADYVOtexaK7r0py1SbVAgAAC7feGuQ1VdXLZ1wHAABsCccUkLNyhz0AANh2jikgd/fbZl0IAABsBevdKGRnkh/Lym2lH5i776T335O8rrvv3JQKAQBgE633Jb03JLk5yctz9/WQdyV5TpLfTfKMuVYGAAALsF5APr+7H3LYvv1JPlRVH59jTQAAsDDrrUH+fFX9QFV96ZiqukdVPSPJ5+dfGgAAbL71AvKFSb4/yWer6uOTWePPJPmXk58BAMC2s96d9PZmss64qs5IUt194ybVBQAAC7HunfQO6e6b5l0IAABsBcd6oxAAANiW1gzIk+sgAwDASWW9EPyhqtqf5B1J3jFZkwwAANvael/S211VS0menOSVVXV2kg8k+aMk7+/uf9ykGgEAYNOsuwa5u/d192909/cmeXSStyV5QpI/raq3b0aBAACwmaZeZ9zddyZ57+SRyYwyAABsK8d8FYvuvm6WhQAAwFbgMm8AG7C8nFTN5rG0NN/O9+9cnv0bcLzWGsPy8qIrA/iSoy6xqKqHdfeVm1EMwFa3b1/SvbU6r1q9pl1Vm1DUBq31Bm7FWoGT1jQzyL9RVR+uqudX1f3mXhEAACzQUQNydz82ybOSnJNkT1W9qaqeOPfKAABgAaZag9zdn0jys0n+9yTfnuS/VNXVVfUv51kcAABstqMG5Kr6pqr6lSRXJXlckqd09zdOnv/KnOsDAIBNNc11kF+d5LVJXtbdtx3a2d3XV9XPzq0yAABYgGkC8ncnua27DyZJVd0jyand/Q/d/Ya5VgcAAJtsmjXI70lyr2H7tMk+AADYdqYJyKd2962HNibPT5tfSQAAsDjTBOS/r6rzD21U1T9Pcts6xwMAwAlrmjXIP5Xk96rq+sn2Vyd5xvxKAgCAxTlqQO7uj1TVNyR5SJJKcnV33zn3ygAAYAGmmUFOkm9Osjw5/hFVle7+nblVBQAAC3LUgFxVb0jytUmuSHJwsruTCMgAAGw708wg707y0O7ueRcDAACLNs1VLK5M8lXzLgQAALaCaWaQz0zysar6cJJ/PLSzu586t6oAAGBBpgnIL593EQAAsFVMc5m391fVUpKv7+73VNVpSXbMvzQAANh8R12DXFU/nuSSJL852XV2kkvnWRQAACzKNF/S+zdJHpPkliTp7k8k+cp5FgUAAIsyTUD+x+6+49BGVe3MynWQAQBg25kmIL+/ql6W5F5V9cQkv5fkbfMtCwAAFmOagPySJAeS/HWSf5XkD5P87DyLAgCARZnmKhZfTPLayQMAALa1owbkqromq6w57u4HzaUiAABYoGluFLJ7eH5qkh9I8oD5lAMAAIt11DXI3X3T8Liuu1+Z5HGbUBsAAGy6aZZYnD9s3iMrM8r3mVtFAACwQNMssfjPw/O7kuxN8oNzqQYAABZsmqtYfOdmFAIAAFvBNEss/u16P+/uV8yuHAAAWKxpr2LxzUn+x2T7KUn+JMm18yoKAAAWZZqAfGaS87v7C0lSVS9P8nvd/b/OszAAAFiEaW41/TVJ7hi270iyPG0HVbWjqv6yqv5gg7UBAMCmm2YG+Q1JPlxVb83KHfWenuR3NtDHC5NcleS+Gy8PAAA21zQ3Cvn5JD+a5PNJbk7yo939H6dpvKp2JfmeJL91PEUCAMBmmWaJRZKcluSW7n5Vkv1Vde6Ur3tlkhcn+eJaB1TV86pqT1XtOXDgwJTNAlvd8nJSdeRjeXlBHS8tLaKZhdqbpVUHsX/H0qpj279j9eP3ZhMGvbRGrTuXF/N7tI6F/W4Dm6a6e/0Dqn4uK1eyeEh3P7iqHpiVL+k95iivuyDJd3f386vqO5L8THdfsN5rdu/e3Xv27NnQAICtqSpZ7c/LWvvn3vEGj597netYs+8N1jqrc7Du8TNtbPrjt+L5WWRNwLGpqsu6e/fh+6eZQX56kqcm+fsk6e7rM92tph+T5KlVtTfJxUkeV1W/O3XFAACwANME5Dt6ZZq5k6Sq7j1Nw9390u7e1d3LSS5M8t7u/uFjrhQAADbBNAH5zVX1m0nuV1U/nuQ9SV4737IAAGAxjnqZt+7+5ap6YpJbkjwkyb/v7ndvpJPufl+S9x1LgQAAsJnWDchVtSPJO7v7CUk2FIoBAOBEtO4Si+4+mOQfquorNqkeAABYqGnupHd7kr+uqndnciWLJOnuF8ytKgAAWJBpAvLbJw8AANj21gzIVfU13f3p7n79ZhYEAACLtN4a5EsPPamq39+EWgAAYOHWC8g1PH/QvAsBAICtYL2A3Gs8BwCAbWu9L+k9vKpuycpM8r0mzzPZ7u6+79yrAwCATbZmQO7uHZtZCAAAbAXr3igEAABONgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQECGbWx5Oak68rG8vOjKVrFWsWs89u9cnk2/S0urtn9NZtT+etYY87U7llcddrL627G0tKGhrXn8Iq11+tey1ti25O/2Wk6oDyicXKq7F13Dl+zevbv37Nmz6DJg26hKVvuIr7V/oX2v8YONHj+zwW3FN2kzalrFut1utKZZnecZlXMsZnZ6tth5hpNRVV3W3bsP328GGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAICBgAwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAzmFpCr6pyq+uOquqqqPlpVL5xXXwAAMCs759j2XUle1N2XV9V9klxWVe/u7o/NsU8AADguc5tB7u4buvvyyfMvJLkqydnz6g8AAGZhU9YgV9Vykkck+YtVfva8qtpTVXsOHDiwGeUAa1heTqqOfCwvb/AFVemsvv+arNXYxuzfsbRq+/t3LG201FUfa9poQ+s9llavNUurj23N4xdorfOw1mNvluY6tGuyvMFf4jXGtXONdqpy7Y6NtQWceKq759tB1elJ3p/k57v7Lesdu3v37t6zZ89c64GTSVWy2kd8VvvX/sExvGaN/et0saF+N9rO7Bo68a035Lm/HYv6BTiGQfsdgxNPVV3W3bsP3z/XGeSqOiXJ7yd549HCMQAAbAXzvIpFJXldkqu6+xXz6gcAAGZpnjPIj0nyI0keV1VXTB7fPcf+AADguM3tMm/d/YEk633VBQAAthx30gMAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAICBgAwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBAhg1aXk6qjnwsLy+6siNdk+VVi70myzNpJ0tLG65p/46lVdvam6WZdLG0evMbL3VmDZ341norNuXtWKvzGX3g9u9cXrX9/TvWGdgaNW30c5WsPrQNm9UfpbXamWVbW/EPJayiunvRNXzJ7t27e8+ePYsuA9ZVlaz2sVlr/0JtsNg1x3AMg9uS7wfbx0Y/iJvxwd1oW7Oqad7tLLImmLOquqy7dx++3wwyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGAjIAAAwEJABAGAgIAMAwEBABgCAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADAQkAEAYCAgAwDAQEAGAICBgAwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAQBgICADAMBAQAYAgIGADAAAAwEZAAAGAjIAAAwEZAAAGMw1IFfVk6rqb6rqk1X1knn2BQAAszC3gFxVO5L8apInJ3lokmdW1UPn1R8AAMzCPGeQH5nkk939qe6+I8nFSZ42x/4AAOC47Zxj22cnuXbY3p/kWw4/qKqel+R5k81bq+pv5lgTdzszyY2LLuJEVbWx/Qt0ZqpWP89rFLvmGI5hcFvw/diuTs7P80Y/iJvxwd1oWxurae3zPKuxrXf8fMfG3U7Oz/PiLK22c54BebVPQB+xo/s1SV4zxzpYRVXt6e7di66D+XKeTw7O88nBeT45OM9bwzyXWOxPcs6wvSvJ9XPsDwAAjts8A/JHknx9VZ1bVfdMcmGS/zHH/gAA4LjNbYlFd99VVT+R5J1JdiT57e7+6Lz6Y8Msazk5OM8nB+f55OA8nxyc5y2guo9YFgwAACctd9IDAICBgAwAAAMBmVTVz1RVV9WZi66F2auqX6qqq6vq/62qt1bV/RZdE7NTVU+qqr+pqk9W1UsWXQ+zV1XnVNUfV9VVVfXRqnrhomtiPqpqR1X9ZVX9waJrOdkJyCe5qjonyROTfHrRtTA3707ysO7+piQfT/LSBdfDjFTVjiS/muTJSR6a5JlV9dDFVsUc3JXkRd39jUm+Ncm/cZ63rRcmuWrRRSAgk/xKkhdnlZu4sD1097u6+67J5oeyck1ytodHJvlkd3+qu+9IcnGSpy24Jmasu2/o7ssnz7+QlQB19mKrYtaqaleS70nyW4uuBQH5pFZVT01yXXf/1aJrYdP8L0n+aNFFMDNnJ7l22N4fwWlbq6rlJI9I8heLrYQ5eGVWJqy+uOhCmO+tptkCquo9Sb5qlR/9uyQvS/Jdm1sR87Deee7u/z455t9l5b9q37iZtTFXtco+/xu0TVXV6Ul+P8lPdfcti66H2amqC5L8bXdfVlXfseh6EJC3ve5+wmr7q+qfJTk3yV9VVbLy3+6XV9Uju/szm1giM7DWeT6kqp6T5IIkj28XP99O9ic5Z9jeleT6BdXCHFXVKVkJx2/s7rcsuh5m7jFJnlpV353k1CT3rarf7e4fXnBdJy03CiFJUlV7k+zu7hsXXQuzVVVPSvKKJN/e3QcWXQ+zU1U7s/LFy8cnuS7JR5L8kLuWbi+1Movx+iSf6+6fWnQ9zNdkBvlnuvuCRddyMrMGGba/Vye5T5J3V9UVVfUbiy6I2Zh8+fInkrwzK1/cerNwvC09JsmPJHnc5DN8xWSmEZgTM8gAADAwgwwAAAMBGQAABgIyAAAMBGQAABgIyAAAMBCQAbaoqnp6VXVVfcNRjvvNqnrMYfteXlXXTS4JdnVV/XpV+ZsPMAV/LAG2rmcm+UCSC49y3Lck+dAq+3+lu89L8tAk/yzJt8+2PIDtSUAG2IKq6vSs3CDix7JOQK6qb0zy8e4+uE5z98zK7Ws/P9MiAbYpARlga/reJO/o7o8n+VxVnb/GcU9O8o41fvbTVXVFkhuyEqKvmEOdANuOgAywNT0zycWT5xdPtlfzL7J2QD60xOIrk9y7qo62VAOAJDsXXQAAX66qzkjyuCQPq6pOsiNJV9WLu7uH405Lcr/uvn699rr7zqp6R5Jvy92hG4A1mEEG2Hq+P8nvdPdSdy939zlJrkny2MOO+84kf3y0xqqqkjw6yf8380oBtiEBGWDreWaStx627/eT/NBh+9Zbf5zcvQb5yqz8j+GvzaxCgG2shv+tA+AEUlWXJ/mW7r5z0bUAbCcCMgAADCyxAACAgYAMAAADARkAAAYCMgAADARkAAAYCMgAADD4/wHwpeKhompumgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for iexp in range(N_exp):\n", " if ((iexp+1)%1000 == 0):\n", " print(f\"Got to experiment number: {iexp+1}\")\n", "\n", " # Generate data:\n", " x_A_array = r.normal(dist_mean_A, dist_width_A, N_events_A)\n", " x_B_array = r.normal(dist_mean_B, dist_width_B, N_events_B)\n", " \n", " \n", " # Test if there is a difference in the mean:\n", " # ------------------------------------------\n", " # Calculate mean and error on mean:\n", " mean_A, width_A, sdom_A = mean_std_sdom(x_A_array) \n", " mean_B, width_B, sdom_B = mean_std_sdom(x_B_array) \n", "\n", " # Consider the difference between means in terms of the uncertainty:\n", " d_mean = mean_A - mean_B\n", " # ... how many sigmas is that away?\n", "\n", " # Turn a number of sigmas into a probability (i.e. p-value):\n", " p_mean = 0.5 # Calculate yourself. HINT: \"stats.norm.cdf or stats.norm.sf may be useful!\"\n", " all_p_mean[iexp] = p_mean\n", " \n", " \n", " # Test if there is a difference with the chi2:\n", " # --------------------------------------------\n", " # Chi2 Test:\n", " p_chi2 = 0.5 # Calculate the p-value of the Chi2 between histograms of A and B yourself.\n", " all_p_chi2[iexp] = p_chi2\n", "\n", " \n", " # Test if there is a difference with the Kolmogorov-Smirnov test on arrays (i.e. unbinned):\n", " # -----------------------------------------------------------------------------------------\n", " p_ks = stats.ks_2samp(x_A_array, x_B_array)[1] # Fortunately, the K-S test is implemented in stats!\n", " all_p_ks[iexp] = p_ks\n", "\n", "\n", " # Print the results for the first 10 experiments\n", " if (verbose and iexp < 10) :\n", " print(f\"{iexp:4d}: p_mean: {p_mean:7.5f} p_chi2: {p_chi2:7.5f} p_ks: {p_ks:7.5f}\")\n", "\n", " \n", " # In case one wants to plot the distribution for visual inspection:\n", " if (iexp == 0):\n", " ax1.hist(x_A_array, N_bins, (xmin, xmax), histtype='step', label='A', color='blue')\n", " ax1.set(title='Histograms of A and B', xlabel='A / B', ylabel='Frequency / 0.05') \n", " ax1.hist(x_B_array, N_bins, (xmin, xmax), histtype='step', label='B', color='red')\n", " ax1.legend()\n", " fig1.tight_layout()\n", " \n", "fig1" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 2 }, "source": [ "## Show the distribution of fit p-values:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAPoCAYAAADDVV/dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5RlVX32++8DLVFEuTX6anNplFaDJip2FNQYYxMDeIHkQISAgCHhnFeDbzSJoiGBxEQxx2jwNV4wEC4aFEmOoGCU4AX1CNqoA7louoNIN9cWaC6CQsvv/WPPSjZF1a5N9dpVvYvvZ4wae6255trrt7vWaPphzjV3qgpJkiRJ0sbbbL4LkCRJkqSFwoAlSZIkSR0xYEmSJElSRwxYkiRJktQRA5YkSZIkdcSAJUmSJEkdMWBJ0gKX5MokL53vOsZZki8n+f1Znntakr8ecPzuJE+Z3DfJryb5wewqliTNFwOWJI2xJNcm2XtS25FJvjaxX1XPrKovz/A+S5NUkkUjKlXTqKqtquqaKdq/WlVPn9if6nctSdr0GLAkSSO3KQe3Tbk2SdL4MWBJ0gLXP/KR5PlJVia5M8nNSd7bul3cXte3KWt7JdksyXFJfpTkliRnJNm6730Pb8duTfLnk65zQpJzknwsyZ3Ake3a30iyPsmNST6QZIu+96skr0+yKsldSd6R5KntnDuTnN3ff4bPfFqSDye5sL3XV5LsMulab0iyCljV2l6Y5FtJ7mivL5z0tk9N8s12/Nwk2/W936eS3NSOXZzkmZPOXTxDLbtN8RlemmRt2z4T2Bn4TPv9vCXJ+UmOmXTO5UkOmOK9JkYoX5dkTZLbk/w/SX6lnbM+yQf6+p+Q5GNTnG8YlaQZGLAk6ZHlJOCkqno88FTg7Nb+kva6TZuy9g3gyPbz68BTgK2ADwAk2R34IHAo8CRga2DJpGvtD5wDbAN8HPg58CZgMbAXsAJ4/aRz9gGeB+wJvAU4uV1jJ+BZwCEP47MeCryjXe+7rYZ+BwAvAHZvYel84P3A9sB7gfOTbN/X/3Dg94AnAxta3wmfA5YBTwC+PcW1ZqploKp6LXAd8Kr2+/lb4HTgsIk+SZ5N73dwwYC3ekGr8zXA3wN/BuwNPBP4nSS/9nDqkiQ9lAFLksbfp9sIxPok6+kFn+ncD+yWZHFV3V1Vlwzoeyjw3qq6pqruBt4GHNxGMQ4EPlNVX6uq+4C/AGrS+d+oqk9X1QNVdW9VXVZVl1TVhqq6FvgIMPkf9O+uqjur6krgCuAL7fp30Asxzx3ujwSA86vq4qr6Gb0gsVeSnfqOv6uqbquqe4FXAKuq6sxW31nA94FX9fU/s6quqKqfAH9OL5BsDlBVp1bVXe1aJwDP7h/tG6KW2TgXWJZkWdt/LfDJ9vuYzjuq6qdV9QXgJ8BZVXVLVV0PfJWH9+crSZqCAUuSxt8BVbXNxA8PHRXqdxTwNOD7bRrcKwf0fTLwo779HwGLgCe2Y2smDlTVPcCtk85f07+T5GlJPtum0t0JvJPeiE6/m/u2751if6sB9U7WX9/dwG2t7qnqm/xZafv9o3JrJh17FL2pf5snOTHJf7bPdW3rs3iqc6ep5WFrYe1s4LAkm9Eb3TtzhtO6/POVJE3BgCVJjyBVtaqqDqE3le3dwDlJHstDR58AbgB26dvfmd7UuJuBG4EdJw4keQy9qXUPutyk/Q/RGxVa1qYovh3I7D/NjP5rhCjJVsB29D7TVPVN/qzQ+7zXT/V+7dj9wI+B36U3HXJvelMll05c9mHUMoypfken0xtpXAHc06Z2duEnwJZ9+/+jo/eVpAXPgCVJjyBJDkuyQ1U9AKxvzT8H1gEP0HvWasJZwJuS7NpCwTvpTUHbQO/Zqle1hSG2AP6SmcPS44A7gbuTPAP4nxv5WSqDv99rvyQvbvW9A7i0qtZM0/cC4GlJfjfJoiSvAXYHPtvX57AkuyfZEvgr4Jyq+nn7XD+jN4K3Jb0/p42pZTo38+DfDy1QPQD8HTOPXj0c3wVekmTnNtXxbR2+tyQtaAYsSXpk2Qe4Msnd9Ba8OLg9k3MP8DfA19uzXHsCp9L7R/vFwA+BnwLHALRnpI4BPkFvNOsu4BZ6QWM6f0JvtOcu4KPAJ2f7IZLsCNwNfG9At38Gjqc3He959EZ6plRVtwKvBP6YXlB6C/DKqvpxX7czgdOAm4BHA29s7WfQmzJ4PXAVMNVzbUPXMsC7gOPa7+dP+trPAH4J+NjUpz18VXUhvd/P5cBlPDhoSpIGSNVUMw4kSRpeG+FaT2/63w/n4HqHAc+sqilHVpKcBqytquNGXct8S3I4cHRVvXi+a5Ek9R5WliTpYUvyKuAielMD30NvNOnaubh2VXU2WjPO2nTF1zN45UhJ0hxyiqAkabb2p7dQww30vlvp4HJaxJxJ8pv0np27md4UREnSJsApgpIkSZLUEUewJEmSJKkjC/IZrMWLF9fSpUvnuwxJkiRJY+qyyy77cVXt8HDPW5ABa+nSpaxcuXK+y5AkSZI0ppL8aDbnOUVQkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSerIyAJWklOT3JLkir627ZJcmGRVe922tSfJ+5OsTnJ5kj36zjmi9V+V5IhR1StJkiRJG2uUI1inAftMajsWuKiqlgEXtX2AfYFl7edo4EPQC2TA8cALgOcDx0+EMkmSJEna1IwsYFXVxcBtk5r3B05v26cDB/S1n1E9lwDbJHkS8JvAhVV1W1XdDlzIQ0ObJEmSJG0SFs3x9Z5YVTcCVNWNSZ7Q2pcAa/r6rW1t07U/RJKj6Y1+sfPOO3dctiRpIXnRiV/k+vX3zvr8Jds8hq8f+7IOK5IkLRRzHbCmkynaakD7QxurTgZOBli+fPmUfSRJArh+/b1ce+IrZn3+0mPP77AaSdJCMterCN7cpv7RXm9p7WuBnfr67QjcMKBdkiRJkjY5cx2wzgMmVgI8Aji3r/3wtprgnsAdbSrh54GXJ9m2LW7x8tYmSZIkSZuckU0RTHIW8FJgcZK19FYDPBE4O8lRwHXAQa37BcB+wGrgHuB1AFV1W5J3AN9q/f6qqiYvnCFJkiRJm4SRBayqOmSaQyum6FvAG6Z5n1OBUzssTZIkSZJGYq6nCEqSJEnSgmXAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6si8BKwkb0pyZZIrkpyV5NFJdk1yaZJVST6ZZIvW9xfa/up2fOl81CxJkiRJM5nzgJVkCfBGYHlVPQvYHDgYeDfwvqpaBtwOHNVOOQq4vap2A97X+kmSJEnSJme+pgguAh6TZBGwJXAj8DLgnHb8dOCAtr1/26cdX5Ekc1irJEmSJA1lzgNWVV0PvAe4jl6wugO4DFhfVRtat7XAkra9BFjTzt3Q+m8/+X2THJ1kZZKV69atG+2HkCRJkqQpzMcUwW3pjUrtCjwZeCyw7xRda+KUAcf+u6Hq5KpaXlXLd9hhh67KlSRJkqShzccUwb2BH1bVuqq6H/hX4IXANm3KIMCOwA1tey2wE0A7vjVw29yWLEmSJEkzm4+AdR2wZ5It27NUK4CrgC8BB7Y+RwDntu3z2j7t+Ber6iEjWJIkSZI03+bjGaxL6S1W8W3ge62Gk4G3Am9OspreM1antFNOAbZv7W8Gjp3rmiVJkiRpGItm7tK9qjoeOH5S8zXA86fo+1PgoLmoS5IkSZI2xnwt0y5JkiRJC44BS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6MjBgJXlGkhVJtprUvs9oy5IkSZKk8TNtwEryRuBc4BjgiiT79x1+56gLkyRJkqRxs2jAsT8AnldVdydZCpyTZGlVnQRkLoqTJEmSpHEyKGBtXlV3A1TVtUleSi9k7YIBS5IkSZIeYtAzWDclec7ETgtbrwQWA7806sIkSZIkadwMCliHAzf1N1TVhqo6HHjJSKuSJEmSpDE07RTBqlqbnhcAS4ACbgC+WVVfn6sCJUmSJGlcTBuwkrwc+CCwCri+Ne8I7Jbk9VX1hTmoT5IkSZLGxqBFLk4C9q6qa/sbk+wKXAD84gjrkiRJkqSxM+gZrEXA2inarwceNZpyJEmSJGl8DRrBOhX4VpJPAGta207AwcApoy5MkiRJksbNoEUu3pXk08D+wF70vvtqLXBoVV01R/VJkiRJ0tgYNIJFVV0NXD1HtUiSJEnSWBv0DNa0kpzQcR2SJEmSNPZmFbCAyzqtQpIkSZIWgFkFrKr6TNeFSJIkSdK4G/RFw4uAo4DfAp4MFHADcC5wSlXdPycVSpIkSdKYGDSCdSbwHOAEYD/gFcBfAs8GPrYxF02yTZJzknw/ydVJ9kqyXZILk6xqr9u2vkny/iSrk1yeZI+NubYkSZIkjcqgVQT3qKqnT2pbC1yS5D828ronAf9WVQcm2QLYEng7cFFVnZjkWOBY4K3AvsCy9vMC4EPtVZIkSZI2KYNGsG5PclCS/+qTZLMkrwFun+0FkzweeAnty4qr6r6qWk/v+7ZOb91OBw5o2/sDZ1TPJcA2SZ402+tLkiRJ0qgMClgHAwcCNyf5jzZqdRPw2+3YbD0FWAf8U5LvJPnHJI8FnlhVNwK01ye0/kuANX3nr21tD5Lk6CQrk6xct27dRpQnSZIkSbMz7RTBqroWeA1Aku2BVNWPO7rmHsAxVXVpkpPoTQecTqYq7yENVScDJwMsX778IcclSZIkadSGWqa9qm7tKFxBbwRqbVVd2vbPoRe4bp6Y+tdeb+nrv1Pf+TvSW81QkiRJkjYps/2i4VmrqpuANUkmFtBYAVwFnAcc0dqOoLccPK398Laa4J7AHRNTCSVJkiRpUzLwe7CqasOIrnsM8PG2guA1wOvohb2zkxwFXAcc1PpeQG+Z+NXAPa2vJEmSJG1yBi3TfkmStcC/0VtS/dquLlpV3wWWT3FoxRR9C3hDV9eWJEmSpFEZtMjF8iS70Pseqr9PsgT4GvA54CtV9bM5qlGSJEmSxsLAZ7Cq6kdV9eGqOgB4IfAZYG/gq0nOn4sCJUmSJGlcDJoi+CBVdT/wxfZDG9GSJEmSJDWzXkWwqq7vshBJkiRJGndzvky7JEmSJC1UMwasJM+ai0IkSZIkadwNM4L14STfTPL6JNuMvCJJkiRJGlMzBqyqejFwKLATsDLJPyf5jZFXJkmSJEljZqhnsKpqFXAc8Fbg14D3J/l+kt8eZXGSJEmSNE6GeQbrl5O8D7gaeBnwqqr6xbb9vhHXJ0mSJEljY5jvwfoA8FHg7VV170RjVd2Q5LiRVSZJkiRJY2aYgLUfcG9V/RwgyWbAo6vqnqo6c6TVSZIkSdIYGeYZrH8HHtO3v2VrkyRJkiT1GSZgPbqq7p7Yadtbjq4kSZIkSRpPwwSsnyTZY2InyfOAewf0lyRJkqRHpGGewfoj4FNJbmj7TwJeM7qSJEmSJGk8zRiwqupbSZ4BPB0I8P2qun/klUmSJEnSmBlmBAvgV4Clrf9zk1BVZ4ysKkmSJEkaQzMGrCRnAk8Fvgv8vDUXYMCSJEmSpD7DjGAtB3avqhp1MZIkSZI0zoZZRfAK4H+MuhBJkiRJGnfDjGAtBq5K8k3gZxONVfXqkVUlSZIkSWNomIB1wqiLkCRJkqSFYJhl2r+SZBdgWVX9e5Itgc1HX5okSZIkjZcZn8FK8gfAOcBHWtMS4NOjLEqSJEmSxtEwi1y8AXgRcCdAVa0CnjDKoiRJkiRpHA0TsH5WVfdN7CRZRO97sCRJkiRJfYYJWF9J8nbgMUl+A/gU8JnRliVJkiRJ42eYgHUssA74HvB/AxcAx42yKEmSJEkaR8OsIvgA8NH2I0mSJEmaxowBK8kPmeKZq6p6ykgqkiRJkqQxNcwXDS/v2340cBCw3WjKkSRJkqTxNeMzWFV1a9/P9VX198DL5qA2SZIkSRorw0wR3KNvdzN6I1qPG1lFkiRJkjSmhpki+Hd92xuAa4HfGUk1kiRJkjTGhllF8NfnohBJkiRJGnfDTBF886DjVfXe7sqRJEmSpPE17CqCvwKc1/ZfBVwMrBlVUZIkSZI0joYJWIuBParqLoAkJwCfqqrfH2VhkiRJkjRuZlymHdgZuK9v/z5g6UiqkSRJkqQxNswI1pnAN5P8f0ABvwWcMdKqJEmSJGkMDbOK4N8k+Rzwq63pdVX1ndGWJUmSJEnjZ5gpggBbAndW1UnA2iS7jrAmSZIkSRpLMwasJMcDbwXe1poeBXxsYy+cZPMk30ny2ba/a5JLk6xK8skkW7T2X2j7q9vxpRt7bUmSJEkahWFGsH4LeDXwE4CqugF4XAfX/l/A1X377wbeV1XLgNuBo1r7UcDtVbUb8L7WT5IkSZI2OcMErPuqqugtcEGSx27sRZPsCLwC+Me2H+BlwDmty+nAAW17/7ZPO76i9ZckSZKkTcowAevsJB8BtknyB8C/Ax/dyOv+PfAW4IG2vz2wvqo2tP21wJK2vYT2pcbt+B2t/4MkOTrJyiQr161bt5HlSZIkSdLDN2PAqqr30Bs5+hfg6cBfVNX/nu0Fk7wSuKWqLutvnurSQxzrr/PkqlpeVct32GGH2ZYnSZIkSbM2cJn2JJsDn6+qvYELO7rmi4BXJ9kPeDTweHojWtskWdRGqXYEbmj91wI70Vu9cBGwNXBbR7VIkiRJUmcGjmBV1c+Be5Js3dUFq+ptVbVjVS0FDga+WFWHAl8CDmzdjgDObdvntX3a8S+2Z8IkSZIkaZMy4xcNAz8FvpfkQtpKggBV9caOa3kr8Ikkfw18BziltZ8CnJlkNb2Rq4M7vq4kSZIkdWKYgHV+++lcVX0Z+HLbvgZ4/hR9fgocNIrrS5IkSVKXpg1YSXauquuq6vTp+kiSJEmS/tugZ7A+PbGR5F/moBZJkiRJGmuDAlb/8uhPGXUhkiRJkjTuBgWsmmZbkiRJkjSFQYtcPDvJnfRGsh7Ttmn7VVWPH3l1kiRJkjRGpg1YVbX5XBYiSZIkSeNu4BcNS5IkSZKGZ8CSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6MucBK8lOSb6U5OokVyb5X619uyQXJlnVXrdt7Uny/iSrk1yeZI+5rlmSJEmShjEfI1gbgD+uql8E9gTekGR34FjgoqpaBlzU9gH2BZa1n6OBD819yZIkSZI0szkPWFV1Y1V9u23fBVwNLAH2B05v3U4HDmjb+wNnVM8lwDZJnjTHZUuSJEnSjOb1GawkS4HnApcCT6yqG6EXwoAntG5LgDV9p61tbZIkSZK0SZm3gJVkK+BfgD+qqjsHdZ2iraZ4v6OTrEyyct26dV2VKUmSJElDm5eAleRR9MLVx6vqX1vzzRNT/9rrLa19LbBT3+k7AjdMfs+qOrmqllfV8h122GF0xUuSJEnSNOZjFcEApwBXV9V7+w6dBxzRto8Azu1rP7ytJrgncMfEVEJJkiRJ2pQsmodrvgh4LfC9JN9tbW8HTgTOTnIUcB1wUDt2AbAfsBq4B3jd3JYrSZIkScOZ84BVVV9j6ueqAFZM0b+AN4y0KEmSJEnqwLyuIihJkiRJC4kBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSNjE7CS7JPkB0lWJzl2vuuRJEmSpMnGImAl2Rz4B2BfYHfgkCS7z29VkiRJkvRgYxGwgOcDq6vqmqq6D/gEsP881yRJkiRJD7JovgsY0hJgTd/+WuAF/R2SHA0c3XZ/luSKOapNmrAY+PF8F6FHHO+7Wcq75/f8MeY9p/ngfaf58PTZnDQuAStTtNWDdqpOBk4GSLKyqpbPRWHSBO87zQfvO8017znNB+87zYckK2dz3rhMEVwL7NS3vyNwwzzVIkmSJElTGpeA9S1gWZJdk2wBHAycN881SZIkSdKDjMUUwarakOQPgc8DmwOnVtWVA045eW4qkx7E+07zwftOc817TvPB+07zYVb3Xapq5l6SJEmSpBmNyxRBSZIkSdrkGbAkSZIkqSNjHbCS7JPkB0lWJzl2iuO/kOST7filSZbOfZVaaIa4796c5Koklye5KMku81GnFo6Z7rm+fgcmqSQuZayNNsx9l+R32t93Vyb557muUQvPEP+N3TnJl5J8p/13dr/5qFMLR5JTk9wy3Xfopuf97Z68PMkeM73n2AasJJsD/wDsC+wOHJJk90ndjgJur6rdgPcBj9yvhVQnhrzvvgMsr6pfBs4B/nZuq9RCMuQ9R5LHAW8ELp3bCrUQDXPfJVkGvA14UVU9E/ijOS9UC8qQf98dB5xdVc+lt6r0B+e2Si1ApwH7DDi+L7Cs/RwNfGimNxzbgAU8H1hdVddU1X3AJ4D9J/XZHzi9bZ8DrEgy1ZcWS8Oa8b6rqi9V1T1t9xJ639smzdYwf9cBvINemP/pXBanBWuY++4PgH+oqtsBquqWOa5RC88w910Bj2/bW+P3omojVdXFwG0DuuwPnFE9lwDbJHnSoPcc54C1BFjTt7+2tU3Zp6o2AHcA289JdVqohrnv+h0FfG6kFWmhm/GeS/JcYKeq+uxcFqYFbZi/654GPC3J15NckmTQ/wGWhjHMfXcCcFiStcAFwDFzU5oewR7uv/3G43uwpjHVSNTkNeeH6SM9HEPfU0kOA5YDvzbSirTQDbznkmxGbwr0kXNVkB4Rhvm7bhG9KTMvpTdS/9Ukz6qq9SOuTQvXMPfdIcBpVfV3SfYCzmz33QOjL0+PUA87T4zzCNZaYKe+/R156DDxf/VJsojeUPKgIUBpJsPcdyTZG/gz4NVV9bM5qk0L00z33OOAZwFfTnItsCdwngtdaCMN+9/Yc6vq/qr6IfADeoFLmq1h7rujgLMBquobwKOBxXNSnR6phvq3X79xDljfApYl2TXJFvQedDxvUp/zgCPa9oHAF8tvVtbGmfG+a9O1PkIvXPlMgjbWwHuuqu6oqsVVtbSqltJ77u/VVbVyfsrVAjHMf2M/Dfw6QJLF9KYMXjOnVWqhGea+uw5YAZDkF+kFrHVzWqUeac4DDm+rCe4J3FFVNw46YWynCFbVhiR/CHwe2Bw4taquTPJXwMqqOg84hd7Q8Wp6I1cHz1/FWgiGvO/+X2Ar4FNtTZXrqurV81a0xtqQ95zUqSHvu88DL09yFfBz4E+r6tb5q1rjbsj77o+BjyZ5E71pWkf6P8+1MZKcRW+q8+L2bN/xwKMAqurD9J712w9YDdwDvG7G9/SelCRJkqRujPMUQUmSJEnapBiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJj3hJlia5Yr7rkCSNPwOWJEmSJHXEgCVJGitttOn7SU5PcnmSc5JsOanPJ5Ps17d/WpL/q5371STfbj8vnOL9j0zygb79zyZ5adt+eZJvtHM/lWSrEX5USdIYMmBJksbR04GTq+qXgTuB1086/gngNQBJtgBWABcAtwC/UVV7tOPvH/aCSRYDxwF7t/NXAm/eyM8hSVpgDFiSpHG0pqq+3rY/Brx40vHPAS9L8gvAvsDFVXUv8Cjgo0m+B3wK2P1hXHPP1v/rSb4LHAHsshGfQZK0AC2a7wIkSZqFmrS/dQs9AH9RVecl+TLwm/RGqs5qx94E3Aw8m97/ZPzpFO+9gQf/D8hHt9cAF1bVIRtfviRpoXIES5I0jnZOslfbPgT4bFU9p/2c19o/AbwO+FXg81BxSSQAACAASURBVK1ta+DGqnoAeC2w+RTvfS3wnCSbJdkJeH5rvwR4UZLdAJJsmeRpXX8wSdJ4M2BJksbR1cARSS4HtgM+NEWfLwAvAf69qu5rbR9s510CPA34yRTnfR34IfA94D3AtwGqah1wJHBWu+4lwDO6+kCSpIUhVZNnWUiStOlKspTeiNWz5rkUSZIewhEsSZIkSeqII1iSJEmS1BFHsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IWqCRXJnnpfNcxzpJ8Ocnvz/Lc05L89YDjdyd5yuS+SX41yQ9mV/Fo6pMkDc+AJUljKMm1Sfae1HZkkq9N7FfVM6vqyzO8z9IklWTRiErVNKpqq6q6Zor2r1bV0yf2p/pdz4X++pL8aZIrktyV5IdJ/nSu65GkcWHAkiSNzKYc3Dbl2jZBAQ4HtgX2Af4wycHzW5IkbZoMWJK0QPWPfCR5fpKVSe5McnOS97ZuF7fX9W1K2F5JNktyXJIfJbklyRlJtu5738PbsVuT/Pmk65yQ5JwkH0tyJ3Bku/Y3kqxPcmOSDyTZou/9Ksnrk6xqIyTvSPLUds6dSc7u7z/DZz4tyYeTXNje6ytJdpl0rTckWQWsam0vTPKtJHe01xdOetunJvlmO35uku363u9TSW5qxy5O8sxJ5y6eoZbdpvgML02ytm2fCewMfKb9ft6S5Pwkx0w65/IkB0zzZ/LiJP9/+/Nfk+TIvsPbtve7K8mlSZ46VX1V9bdV9e2q2lBVPwDOBV401fUk6ZHOgCVJjwwnASdV1eOBpwJnt/aXtNdt2pSwbwBHtp9fB54CbAV8ACDJ7sAHgUOBJwFbA0smXWt/4BxgG+DjwM+BNwGLgb2AFcDrJ52zD/A8YE/gLcDJ7Ro7Ac8CDnkYn/VQ4B3tet9tNfQ7AHgBsHsLS+cD7we2B94LnJ9k+77+hwO/BzwZ2ND6TvgcsAx4AvDtKa41Uy0DVdVrgeuAV7Xfz98CpwOHTfRJ8mx6v4MLJp+fZOdW4/8GdgCe0+qYcAjwl/RGplYDfzNTTUkC/Cpw5cP5LJL0SGHAkqTx9ek2KrE+yXp6wWc69wO7JVlcVXdX1SUD+h4KvLeqrqmqu4G3AQe3KXUHAp+pqq9V1X3AXwA16fxvVNWnq+qBqrq3qi6rqkva6Me1wEeAX5t0zrur6s6quhK4AvhCu/4d9ALCc4f7IwHg/Kq6uKp+BvwZsFeSnfqOv6uqbquqe4FXAKuq6sxW31nA94FX9fU/s6quqKqfAH8O/E6SzQGq6tSquqtd6wTg2f2jfUPUMhvnAsuSLGv7rwU+2X4fkx0K/HtVnVVV91fVrVXVH7D+taq+WVUb6IW/5wxx/RPo/fvhn2b/ESRp4TJgSdL4OqCqtpn44aGjQv2OAp4GfL9Ng3vlgL5PBn7Ut/8jYBHwxHZszcSBqroHuHXS+Wv6d5I8Lcln21S6O4F30hvR6Xdz3/a9U+xvNaDeyfrruxu4rdU9VX2TPyttv39Ubs2kY4+iN/Vv8yQnJvnP9rmubX0WT3XuNLU8bC2snQ0clmQzeqNQZ07TfSfgPwe83U192/cww59zkj+kN6L3ilaHJGkSA5YkPQJU1aqqOoTeVLZ3A+ckeSwPHX0CuAHYpW9/Z3pT424GbgR2nDiQ5DH0ptY96HKT9j9Eb1RoWZui+HZ6iyaMyn+NECXZCtiO3meaqr7JnxV6n/f6qd6vHbsf+DHwu/SmQ+5Nb6rk0onLPoxahjHV7+h0eqNTK4B72tTOqayhNyV0oyX5PeBYYEVVre3iPSVpITJgSdIjQJLDkuxQVQ8A61vzz4F1wAP0nrWacBbwpiS7tlDwTnpT0DbQe7bqVW1hiC3oPb8zU1h6HHAncHeSZwD/cyM/S2Xw93vt1xZ22ILe80+XVtWaafpeADwtye8mWZTkNcDuwGf7+hyWZPckWwJ/BZxTVT9vn+tn9EbwtqT357QxtUznZh78+6EFqgeAv2P60SvoTfvbO8nvtM+3fZJhpgE+SJJD6X2+35hqaXlJ0n8zYEnSI8M+wJVJ7qa34MXBVfXTNsXvb4Cvt2e59gROpfeP9ouBHwI/BY4BaM9IHQN8gt5o1l3ALfSCxnT+hN5oz13AR4FPzvZDJNkRuBv43oBu/wwcT2863vPojfRMqapuBV4J/DG9oPQW4JVV9eO+bmcCp9GbTvdo4I2t/Qx6UwavB64CpnqubehaBngXcFz7/fxJX/sZwC8BH5vuxKq6DtiP3ue7jd4CF8+eRQ1/TW+k8lttNcO7k3x4Fu8jSQteqqaaeSBJ0szaCNd6etP/fjgH1zsMeGZVvW2a46cBa6vquFHXMt+SHA4cXVUvnu9aJEn/zS9ZlCQ9LEleBVxEb2rge+iNJl07F9euqmlHax5J2nTF1zN45UhJ0jxwiqAk6eHan95CDTfQ+w6og8vpEHMmyW/Se3buZnpTECVJmxCnCEqSJElSRxzBkiRJkqSOLMhnsBYvXlxLly6d7zIkSZIkjanLLrvsx1W1w8M9b0EGrKVLl7Jy5cr5LkOSJEnSmEryo9mc5xRBSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjIwtYSU5NckuSK/ratktyYZJV7XXb1p4k70+yOsnlSfboO+eI1n9VkiNGVa8kSZIkbaxRjmCdBuwzqe1Y4KKqWgZc1PYB9gWWtZ+jgQ9BL5ABxwMvAJ4PHD8RyiRJkiRpUzOygFVVFwO3TWreHzi9bZ8OHNDXfkb1XAJsk+RJwG8CF1bVbVV1O3AhDw1tkiRJkrRJWDTH13tiVd0IUFU3JnlCa18CrOnrt7a1Tdf+EEmOpjf6xc4779xx2ZKkheRFJ36R69ffO+vzl2zzGL5+7Ms6rEiStFDMdcCaTqZoqwHtD22sOhk4GWD58uVT9pEkCeD69fdy7YmvmPX5S489v8NqJEkLyVyvInhzm/pHe72lta8FdurrtyNww4B2SZIkSdrkzHXAOg+YWAnwCODcvvbD22qCewJ3tKmEnwdenmTbtrjFy1ubJEmSJG1yRjZFMMlZwEuBxUnW0lsN8ETg7CRHAdcBB7XuFwD7AauBe4DXAVTVbUneAXyr9furqpq8cIYkSZIkbRJGFrCq6pBpDq2Yom8Bb5jmfU4FTu2wNEmSJEkaibmeIihJkiRJC5YBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOzEvASvKmJFcmuSLJWUkenWTXJJcmWZXkk0m2aH1/oe2vbseXzkfNkiRJkjSTOQ9YSZYAbwSWV9WzgM2Bg4F3A++rqmXA7cBR7ZSjgNurajfgfa2fJEmSJG1y5muK4CLgMUkWAVsCNwIvA85px08HDmjb+7d92vEVSTKHtUqSJEnSUOY8YFXV9cB7gOvoBas7gMuA9VW1oXVbCyxp20uANe3cDa3/9pPfN8nRSVYmWblu3brRfghJkiRJmsJ8TBHclt6o1K7Ak4HHAvtO0bUmThlw7L8bqk6uquVVtXyHHXboqlxJkiRJGtp8TBHcG/hhVa2rqvuBfwVeCGzTpgwC7Ajc0LbXAjsBtONbA7fNbcmSJEmSNLP5CFjXAXsm2bI9S7UCuAr4EnBg63MEcG7bPq/t045/saoeMoIlSZIkSfNtPp7BupTeYhXfBr7XajgZeCvw5iSr6T1jdUo75RRg+9b+ZuDYua5ZkiRJkoaxaOYu3auq44HjJzVfAzx/ir4/BQ6ai7okSZIkaWPM1zLtkiRJkrTgGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMGLEmSJEnqiAFLkiRJkjpiwJIkSZKkjhiwJEmSJKkjBixJkiRJ6ogBS5IkSZI6YsCSJEmSpI4YsCRJkiSpIwYsSZIkSeqIAUuSJEmSOmLAkiRJkqSOGLAkSZIkqSMDA1aSZyRZkWSrSe37jLYsSZIkSRo/0wasJG8EzgWOAa5Isn/f4XeOujBJkiRJGjeLBhz7A+B5VXV3kqXAOUmWVtVJQOaiOEmSJEkaJ4MC1uZVdTdAVV2b5KX0QtYuGLAkSZIk6SEGPYN1U5LnTOy0sPVKYDHwS6MuTJIkSZLGzaCAdThwU39DVW2oqsOBl4y0KkmSJEkaQ9NOEayqtel5AbAEKOAG4JtV9fW5KlCSJEmSxsW0ASvJy4EPAquA61vzjsBuSV5fVV+Yg/okSZIkaWwMWuTiJGDvqrq2vzHJrsAFwC+OsC5JkiRJGjuDnsFaBKydov164FGjKUeSJEmSxtegEaxTgW8l+QSwprXtBBwMnDLqwiRJkiRp3Axa5OJdST4N7A/sRe+7r9YCh1bVVXNUnyRJkiSNjUEjWFTV1cDVc1SLJEmSJI21Qc9gTSvJCR3XIUmSJEljb1YBC7is0yokSZIkaQGYVcCqqs90XYgkSZIkjbtBXzS8CDgK+C3gyUABNwDnAqdU1f1zUqEkSZIkjYlBI1hnAs8BTgD2A14B/CXwbOBjG3PRJNskOSfJ95NcnWSvJNsluTDJqva6beubJO9PsjrJ5Un22JhrS5IkSdKoDFpFcI+qevqktrXAJUn+YyOvexLwb1V1YJItgC2BtwMXVdWJSY4FjgXeCuwLLGs/LwA+1F4lSZIkaZMyaATr9iQHJfmvPkk2S/Ia4PbZXjDJ44GX0L6suKruq6r19L5v6/TW7XTggLa9P3BG9VwCbJPkSbO9viRJkiSNyqCAdTBwIHBzkv9oo1Y3Ab/djs3WU4B1wD8l+U6Sf0zyWOCJVXUjQHt9Quu/BFjTd/7a1vYgSY5OsjLJynXr1m1EeZIkSZI0O9NOEayqa4HXACTZHkhV/bija+4BHFNVlyY5id50wOlkqvIe0lB1MnAywPLlyx9yXJIkSZJGbahl2qvq1o7CFfRGoNZW1aVt/xx6gevmial/7fWWvv479Z2/I73VDCVJkiRpkzLbLxqetaq6CViTZGIBjRXAVcB5wBGt7Qh6y8HT2g9vqwnuCdwxMZVQkiRJkjYlA78Hq6o2jOi6xwAfbysIXgO8jl7YOzvJUcB1wEGt7wX0lolfDdzT+kqSJEnSJmfQMu2XJFkL/Bu9JdWv7eqiVfVdYPkUh1ZM0beAN3R1bUmSJEkalUGLXCxPsgu976H6+yRLgK8BnwO+UlU/m6MaJUmSJGksDHwGq6p+VFUfrqoDgBcCnwH2/j/s3XuYXWV99//3x0Ql1EOqxJaGYFAjij5aaapY2kqDthyUaH9aoR7AUvO0Wq21/dVgeSq/p+3T2FqpPB7xgXJQQaCtoGAtgkilgsbDw0GwpJhCAkI8BEQiCH5/f+w1ujNMZnb2XrP37Mn7dV37mrXutdZe35m9rsl8ct/rXsC/JblwGAVKkiRJ0riYbojgdqrqh8ClzYumR0uSJEmS1Oh7FsGq2txmIZIkSZI07oY+TbskSZIkzVczBqwkTx9GIZIkSZI07nrpwXp/ki8keV2SxbNekSRJkiSNqRkDVlX9MvAKYBmwPslHkrxg1iuTJEmSpDHT0z1YVXUjcDzwFuB5wElJbkjym7NZnCRJkiSNk17uwXpGkhOB64FVwIuq6qnN8omzXJ8kSZIkjY1enoP1buCDwFurattEY1XdmuT4WatMkiRJksZMLwHrMGBbVT0AkOQhwG5VdU9VnTmr1UmSJEnSGOnlHqxPA4u61ndv2iRJkiRJXXoJWLtV1d0TK83y7rNXkiRJkiSNp14C1veT7D+xkuQXgG3T7C9JkiRJu6Re7sF6E3Buklub9T2Bl89eSZIkSZI0nmYMWFX1xSRPAfYFAtxQVT+c9cokSZIkacz00oMF8IvA8mb/ZyWhqs6YtaokSZIkaQzNGLCSnAk8Efgq8EDTXIABS5IkSZK69NKDtRLYr6pqtouRJEmSpHHWyyyC1wI/O9uFSJIkSdK466UHaw/ga0m+ANw70VhVR8xaVZIkSZI0hnoJWCfMdhGSJEmSNB/0Mk37Z5M8HlhRVZ9OsjuwYPZLkyRJkqTxMuM9WEleC5wHfKBpWgp8bDaLkiRJkqRx1MskF68HDgTuAqiqG4HHzWZRkiRJkjSOeglY91bVfRMrSRbSeQ6WJEmSJKlLLwHrs0neCixK8gLgXODjs1uWJEmSJI2fXgLWWmALcA3w34GLgONnsyhJkiRJGke9zCL4I+CDzUuSJEmStAMzBqwk32CKe66q6gmzUpEkSZIkjaleHjS8smt5N+BlwGNmpxxJkiRJGl8z3oNVVd/uem2uqr8HVg2hNkmSJEkaK70MEdy/a/UhdHq0HjlrFUmSJEnSmOpliODfdS3fD2wEfmtWqpEkSZKkMdbLLIK/NoxCJEmSJGnc9TJE8M3Tba+qd7ZXjiRJkiSNr15nEfxF4IJm/UXA5cAts1WUJEmSJI2jXgLWHsD+VfU9gCQnAOdW1e/OZmGSJEmSNG5mnKYd2Bu4r2v9PmD5rFQjSZIkSWOslx6sM4EvJPlnoICXAGfMalWSJEmSNIZ6mUXwr5J8EviVpuk1VfWV2S1LkiRJksZPL0MEAXYH7qqqdwGbkuwzizVJkiRJ0liaMWAleRvwFuC4pumhwIdmsyhJkiRJGke99GC9BDgC+D5AVd0KPHLQEydZkOQrST7RrO+T5KokNyb5aJKHNe0Pb9Y3NNuXD3puSZIkSZoNvQSs+6qq6ExwQZKfauncfwhc37X+duDEqloBfBc4tmk/FvhuVT0JOLHZT5IkSZLmnF4C1jlJPgAsTvJa4NPABwc5aZK9gMOB/9OsB1gFnNfscjrw4mZ5dbNOs/3gZn9JkiRJmlN6mUXwHUleANwF7Av8eVVdPOB5/x74U34y1PCxwNaqur9Z3wQsbZaXArc0tdyf5M5m/291v2GSNcAagL333nvA8iRJkiRp500bsJIsAD5VVc8HBg1VE+/5QuCOqvpSkoMmmqfYtXrY9pOGqpOBkwFWrlz5oO2SJEmSNNumDVhV9UCSe5I8uqrubOmcBwJHJDkM2A14FJ0ercVJFja9WHsBtzb7bwKW0ZkefiHwaOA7LdUiSZIkSa3p5R6sHwDXJDklyUkTr35PWFXHVdVeVbUcOBK4tKpeAXwGeGmz29HA+c3yBc06zfZLm0k3JEmSJGlOmfEeLODC5jXb3gKcneQvga8ApzTtpwBnJtlAp+fqyCHUIkmSJEk7bYcBK8neVXVzVZ2+o30GVVWXAZc1yzcBz55inx8AL5utGiRJkiSpLdMNEfzYxEKSfxxCLZIkSZI01qYLWN2z9z1htguRJEmSpHE3XcCqHSxLkiRJkqYw3SQXz0xyF52erEXNMs16VdWjZr06SZIkSRojOwxYVbVgmIVIkiRJ0rjr5TlYkiRJkqQeGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSVDD1hJliX5TJLrk1yX5A+b9sckuTjJjc3Xn27ak+SkJBuSXJ1k/2HXLEmSJEm9GEUP1v3AH1fVU4EDgNcn2Q9YC1xSVSuAS5p1gEOBFc1rDfC+4ZcsSZIkSTMbesCqqtuq6svN8veA64GlwGrg9Ga304EXN8urgTOq40pgcZI9h1y2JEmSJM1opPdgJVkOPAu4CviZqroNOiEMeFyz21Lglq7DNjVtk99rTZL1SdZv2bJlNsuWJEmSpCmNLGAleQTwj8Cbququ6Xadoq0e1FB1clWtrKqVS5YsaatMSZIkSerZSAJWkofSCVcfrqp/appvnxj613y9o2nfBCzrOnwv4NZh1SpJkiRJvRrFLIIBTgGur6p3dm26ADi6WT4aOL+r/dXNbIIHAHdODCWUJEmSpLlk4QjOeSDwKuCaJF9t2t4KrAPOSXIscDPwsmbbRcBhwAbgHuA1wy1XkiRJknoz9IBVVZ9j6vuqAA6eYv8CXj+rRUmSJElSC0Y6i6AkSZIkzScGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaMjYBK8khSb6eZEOStaOuR5IkSZImG4uAlWQB8B7gUGA/4Kgk+422KkmSJEna3lgELODZwIaquqmq7gPOBlaPuCZJkiRJ2s7CURfQo6XALV3rm4DndO+QZA2wplm9N8m1Q6pNmrAH8K1RF6Fdjtddn/L20R4/xrzmNApedxqFffs5aFwCVqZoq+1Wqk4GTgZIsr6qVg6jMGmC151GwetOw+Y1p1HwutMoJFnfz3HjMkRwE7Csa30v4NYR1SJJkiRJUxqXgPVFYEWSfZI8DDgSuGDENUmSJEnSdsZiiGBV3Z/kD4BPAQuAU6vqumkOOXk4lUnb8brTKHjdadi85jQKXncahb6uu1TVzHtJkiRJkmY0LkMEJUmSJGnOM2BJkiRJUkvGOmAlOSTJ15NsSLJ2iu0PT/LRZvtVSZYPv0rNNz1cd29O8rUkVye5JMnjR1Gn5o+Zrrmu/V6apJI4lbEG1st1l+S3mt931yX5yLBr1PzTw7+xeyf5TJKvNP/OHjaKOjV/JDk1yR07eoZuOk5qrsmrk+w/03uObcBKsgB4D3AosB9wVJL9Ju12LPDdqnoScCKw6z4WUq3o8br7CrCyqp4BnAf8zXCr1HzS4zVHkkcCbwSuGm6Fmo96ue6SrACOAw6sqqcBbxp6oZpXevx9dzxwTlU9i86s0u8dbpWah04DDplm+6HAiua1BnjfTG84tgELeDawoapuqqr7gLOB1ZP2WQ2c3iyfBxycZKqHFku9mvG6q6rPVNU9zeqVdJ7bJvWrl991AH9BJ8z/YJjFad7q5bp7LfCeqvouQFXdMeQaNf/0ct0V8Khm+dH4XFQNqKouB74zzS6rgTOq40pgcZI9p3vPcQ5YS4FbutY3NW1T7lNV9wN3Ao8dSnWar3q57rodC3xyVivSfDfjNZfkWcCyqvrEMAvTvNbL77onA09OckWSK5NM9z/AUi96ue5OAF6ZZBNwEfCG4ZSmXdjO/u03Hs/B2oGpeqImzznfyz7Szuj5mkrySmAl8LxZrUjz3bTXXJKH0BkCfcywCtIuoZffdQvpDJk5iE5P/b8leXpVbZ3l2jR/9XLdHQWcVlV/l+S5wJnNdfej2S9Pu6idzhPj3IO1CVjWtb4XD+4m/vE+SRbS6UqergtQmkkv1x1Jng/8GXBEVd07pNo0P810zT0SeDpwWZKNwAHABU50oQH1+m/s+VX1w6r6BvB1OoFL6lcv192xwDkAVfV5YDdgj6FUp11VT3/7dRvngPVFYEWSfZI8jM6NjhdM2ucC4Ohm+aXApeWTlTWYGa+7ZrjWB+iEK+9J0KCmveaq6s6q2qOqllfVcjr3/R1RVetHU67miV7+jf0Y8GsASfagM2TwpqFWqfmml+vuZuBggCRPpROwtgy1Su1qLgBe3cwmeABwZ1XdNt0BYztEsKruT/IHwKeABcCpVXVdkv8JrK+qC4BT6HQdb6DTc3Xk6CrWfNDjdfe3wCOAc5s5VW6uqiNGVrTGWo/XnNSqHq+7TwG/nuRrwAPA/1tV3x5d1Rp3PV53fwx8MMkf0RmmdYz/ea5BJDmLzlDnPZp7+94GPBSgqt5P516/w4ANwD3Aa2Z8T69JSZIkSWrHOA8RlCRJkqQ5xYAlSZIkSS0xYEmSJElSSwxYkiRJktQSA5YkSZIktcSAMDPtNAAAIABJREFUJUna5SVZnuTaUdchSRp/BixJkiRJaokBS5I0VprephuSnJ7k6iTnJdl90j4fTXJY1/ppSf6f5th/S/Ll5vVLU7z/MUne3bX+iSQHNcu/nuTzzbHnJnnELH6rkqQxZMCSJI2jfYGTq+oZwF3A6yZtPxt4OUCShwEHAxcBdwAvqKr9m+0n9XrCJHsAxwPPb45fD7x5wO9DkjTPGLAkSePolqq6oln+EPDLk7Z/EliV5OHAocDlVbUNeCjwwSTXAOcC++3EOQ9o9r8iyVeBo4HHD/A9SJLmoYWjLkCSpD7UpPVHN6EH4M+r6oIklwG/Qaen6qxm2x8BtwPPpPOfjD+Y4r3vZ/v/gNyt+Rrg4qo6avDyJUnzlT1YkqRxtHeS5zbLRwGfqKqfb14XNO1nA68BfgX4VNP2aOC2qvoR8CpgwRTvvRH4+SQPSbIMeHbTfiVwYJInASTZPcmT2/7GJEnjzYAlSRpH1wNHJ7kaeAzwvin2+VfgV4FPV9V9Tdt7m+OuBJ4MfH+K464AvgFcA7wD+DJAVW0BjgHOas57JfCUtr4hSdL8kKrJoywkSZq7kiyn02P19BGXIknSg9iDJUmSJEktsQdLkiRJklpiD5YkSZIktcSAJUmSJEktMWBJkiRJUksMWJIkSZLUEgOWJEmSJLXEgCVJkiRJLTFgSZIkSVJLDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEnSiCS5LslBo65jnCW5LMnv9nnsaUn+cprtdyd5wuR9k/xKkq/3V3F7kixPUkkWjroWSdJPGLAkaRYk2Zjk+ZPajknyuYn1qnpaVV02w/v4R/SIVNUjquqmKdr/rar2nVif6rPu1eTPNx3/O8kNSZb2X70kaVQMWJK0C5vLwW0u1zYbkgT4AHAQ8Lyq2jzaiga3q32GkgQGLEkame6ejyTPTrI+yV1Jbk/yzma3y5uvW5sha89N8pAkxyf5ryR3JDkjyaO73vfVzbZvJ/kfk85zQpLzknwoyV3AMc25P59ka5Lbkrw7ycO63q+SvC7JjUm+l+QvkjyxOeauJOd07z/D93xakvcnubh5r88mefykc70+yY3AjU3bLyX5YpI7m6+/NOltn5jkC83285M8puv9zk3yzWbb5UmeNunYPWao5UlTfA8HJdnULJ8J7A18vPl8/jTJhUneMOmYq5O8eJofzQLgNGAlcFBV3d4cN+1nPekclyX5yyT/3tTy8SSPTfLh5nP6YpLlXfvv8OeaZJ/m5/W9JJ9O8p4kH+rafkQ6Q1y3Nud9ate2jUnekuRq4PtJFiZ5arPf1ua4I5p9D2g+nwVdx7+kOVaSxpIBS5LmhncB76qqRwFPBM5p2n+1+bq4GbL2eeCY5vVrwBOARwDvBkiyH/Be4BXAnsCjgclDzVYD5wGLgQ8DDwB/BOwBPBc4GHjdpGMOAX4BOAD4U+Dk5hzLgKcDR+3E9/oK4C+a8321qaHbi4HnAPs1YelC4CTgscA7gQuTPLZr/1cDvwP8HHB/s++ETwIrgMcBX57iXDPVMq2qehVwM/Ci5vP5G+B04JUT+yR5Jp3P4KJp3urDwFOAVVX17a72Y9jBZ70DRwKvas73RODzwD8AjwGuB97W1DTTz/UjwBeabSc07znx/TwZOAt4E7Ck+b4+PilkHwUcTucaC/Bx4F/pfA5vAD6cZN+quhL4PrCq69jfbs4vSWPJgCVJs+djzf/Yb02ylU7w2ZEfAk9KskdV3d384bkjrwDeWVU3VdXdwHHAkekMx3op8PGq+lxV3Qf8OVCTjv98VX2sqn5UVduq6ktVdWVV3V9VG+kMU3vepGPeXlV3VdV1wLXAvzbnv5NOiHlWbz8SAC6sqsur6l7gz4DnJlnWtf2vq+o7VbWNzh/pN1bVmU19ZwE3AC/q2v/Mqrq2qr4P/A/gtyZ6RKrq1Kr6XnOuE4BnTuoBmqmWfpwPrEiyoll/FfDR5vPYkV8HzqmqrZPap/usp/IPVfWfXZ/Lf1bVp6vqfuBcfvI57fDnmmRv4BeBP6+q+6rqc8AFXed4OZ2f28VV9UPgHcAioLtn8aSquqX5DA+gEwzXNe93KfAJfhLKz5pYTvJI4LCmTZLGkgFLkmbPi6tq8cSLB/cKdTsWeDJwQzNc64XT7PtzwH91rf8XsBD4mWbbLRMbquoe4Nts75bulSRPTvKJZqjWXcD/otOj0+32ruVtU6w/Ypp6J+uu727gO03dU9U3+XulWe/ulbtl0raH0hn6tyDJuiT/2XxfG5t99pjq2B3UstOasHYO8MokD6ETHs6c4bAXAm9L8juT2qf7rKfS6+c03c/154DvNNfOhB1+JlX1o2b7jj6TnwNuafabfC7o9Fb9ZpKHA78JfLmqJtcmSWPDgCVJc0BV3VhVR9EZQvV24LwkP8WDe58AbgUe37W+N52hcbcDtwF7TWxIsojOMK/tTjdp/X10ei9WNEMU30pnWNds+XEPUZJH0Bm+dusO6pv8vULn++2eAGLZpG0/BL5FZ6jZauD5dIZKLp847U7U0oupPqPT6fQ+HQzc0wztnM6/0+mVe1eS3+5qn+6zHsR0P9fbgMck2b1rW/fPeLtjk6TZ3v2ZTP4MlzVhc/K5qKqv0Qlch+LwQEnzgAFLkuaAJK9MsqT5X/6JYWIPAFuAH9G5/2bCWcAfNRMRPIJOj9NHm2Fg59EZ5vVLzT0x/x8zh6VHAncBdyd5CvD7A34vlemf73VYkl9u6vsL4KqqumUH+14EPDnJbzeTJbwc2I/OELMJr0yyXxMI/idwXlU90Hxf99Lpwdudzs9pkFp25Ha2/3xoAtWPgL9j5t6riWM+S6cH5+QkL22ap/usB7HDn2vTe7QeOCHJw5I8l+2HZJ4DHJ7k4CQPBf6Yzs/533dwrqvo3Gf1p0ke2lwbLwLO7trnI8Ab6dxzeO6A35skjZQBS5LmhkOA65LcTWfCiyOr6gfNMK2/Aq5o7uU6ADiVzh/tlwPfAH5AZ+IAmnuk3kDnj9fbgO8Bd9D5A3hH/oROz8H3gA8CH+33m0iyF3A3cM00u32EzmQL36EzccYrdrRjM+HDC+n8Ef9tOhNsvLCqvtW125l0ZuD7JrAbnT/UAc6g0zOyGfgaMNV9bT3XMo2/Bo5vPp8/6Wo/A/hvwIemPuzBqupiOvc4nZbkRUzzWQ+ih5/rK+hMePJt4C/pXBP3Nsd+nc4kHv+bTk/hi+hM8jHlPWZN+xF0eqi+RedexFdX1Q1du51FZ3r6Syd9tpI0dlI11cgGSdJ80PR6bKUz/O8bQzjfK4GnVdVxO9h+GrCpqo6f7VpGLcmrgTVV9cujrmVQST4K3FBVbxt1LZI01/kAQEmaZ5qej0voDA18B53epI3DOHdV9dxbM581wxVfx/QzR85ZSX6RTq/eN+jMcLgaWDfSoiRpTDhEUJLmn9V0Jha4lc4zoI4shysMTZLfoHPv3O2M74QNPwtcRme450nA71fVV0ZakSSNCYcISpIkSVJL7MGSJEmSpJaM5B6sJKfSmb3ojqp6+qRtfwL8LbCkqr7VPF/jXXSe7H4PcExVfXm6999jjz1q+fLls1K7JEmSpPnvS1/60reqasnOHjeqSS5OA95NZwrbH0uyDHgBcHNX86F07iFYATyHzgMxnzPdmy9fvpz169e3WK4kSZKkXUmS/+rnuJEMEayqy+nMTjTZiXSexdF9Y9hq4IzquBJYnGTPIZQpSZIkSTtlztyDleQIYHNV/d9Jm5YCt3Stb2raJEmSJGlOmRPPwWqeF/JndJ618aDNU7Q9aOrDJGuANQB77713q/VJkiRJUi/mSg/WE4F9gP+bZCOwF/DlJD9Lp8dqWde+e9F5tst2qurkqlpZVSuXLNnpe9EkSZIkaWBzImBV1TVV9biqWl5Vy+mEqv2r6pvABcCr03EAcGdV3TbKeiVJkiRpKiMJWEnOAj4P7JtkU5Jjp9n9IuAmYAPwQeB1QyhRkiRJknbaSO7BqqqjZti+vGu5gNfPdk2SJEmSNKg5MURQkiRJkuYDA5YkSZIktcSAJUmSJEktMWBJkiRJUkvmxIOGJUmS5rsD113K5q3b+j5+6eJFXLF2VYsVSZoNBixJkqQh2Lx1GxvXHd738cvXXthiNZJmi0MEJUmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklgw9YCU5NckdSa7tavvbJDckuTrJPydZ3LXtuCQbknw9yW8Mu15JkiRJ6tUoerBOAw6Z1HYx8PSqegbwH8BxAEn2A44EntYc894kC4ZXqiRJkiT1bugBq6ouB74zqe1fq+r+ZvVKYK9meTVwdlXdW1XfADYAzx5asZIkSZK0E+biPVi/A3yyWV4K3NK1bVPT9iBJ1iRZn2T9li1bZrlESZIkSXqwORWwkvwZcD/w4YmmKXarqY6tqpOramVVrVyyZMlslShJkiRJO7Rw1AVMSHI08ELg4KqaCFGbgGVdu+0F3Drs2iRJkiSpF3OiByvJIcBbgCOq6p6uTRcARyZ5eJJ9gBXAF0ZRoyRJkiTNZOg9WEnOAg4C9kiyCXgbnVkDHw5cnATgyqr6vaq6Lsk5wNfoDB18fVU9MOyaJUmSJKkXQw9YVXXUFM2nTLP/XwF/NXsVSZIkSVI75sQQQUmSJEmaDwxYkiRJktQSA5YkSZIktcSAJUmSJEktmTPPwZIkSZKmcuC6S9m8dVvfxy9dvIgr1q5qsSJpxwxYkiRJmtM2b93GxnWH93388rUXtliNND2HCEqSJElSSwxYkiRJktQSA5YkSZIktcSAJUmSJEktMWBJkiRJUksMWJIkSZLUEgOWJEmSJLXEgCVJkiRJLTFgSZIkSVJLDFiSJEmS1JKFoy5AkiQN14HrLmXz1m19Hbt08SKuWLuq5Yokaf4wYEmStIvZvHUbG9cd3texy9de2HI1kjS/OERQkiRJklpiwJIkSZKklhiwJEmSJKklBixJkiRJaokBS5IkSZJaYsCSJEmSpJYYsCRJkiSpJQYsSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqyUgCVpJTk9yR5NqutsckuTjJjc3Xn27ak+SkJBuSXJ1k/1HULEmSJEkzGVUP1mnAIZPa1gKXVNUK4JJmHeBQYEXzWgO8b0g1SpIkSdJOGUnAqqrLge9Mal4NnN4snw68uKv9jOq4ElicZM/hVCpJkiRJvZtL92D9TFXdBtB8fVzTvhS4pWu/TU3bdpKsSbI+yfotW7bMerGSJEmSNNlcClg7kina6kENVSdX1cqqWrlkyZIhlCVJkiRJ25tLAev2iaF/zdc7mvZNwLKu/fYCbh1ybZIkSZI0o7kUsC4Ajm6WjwbO72p/dTOb4AHAnRNDCSVJkiRpLlk4ipMmOQs4CNgjySbgbcA64JwkxwI3Ay9rdr8IOAzYANwDvGboBUuSJElSD0YSsKrqqB1sOniKfQt4/exWJEmSJEmDm0tDBCVJkiRprBmwJEmSJKklIxkiKEnSgesuZfPWbX0fv3TxIq5Yu6rFiiRJGpwBS5I0Epu3bmPjusP7Pn752gtbrEaSpHY4RFCSJEmSWmLAkiRJkqSWGLAkSZIkqSUGLEmSJElqiQFLkiRJklpiwJIkSZKklvQVsJI8JcnBSR4xqf2QdsqSJEmSpPGz0wEryRuB84E3ANcmWd21+X+1VZgkSZIkjZt+HjT8WuAXquruJMuB85Isr6p3AWmzOEmSJEkaJ/0ErAVVdTdAVW1MchCdkPV4DFiSJEnSnHHgukvZvHVb38cvXbyIK9auarGi+a+fgPXNJD9fVV8FaHqyXgicCvy3VquTJEmS1LfNW7excd3hfR+/fO2FLVaza+hnkotXA9/sbqiq+6vq1cCvtlKVJEmSJI2hne7BqqpN6XgOsBQo4FbgC1V1RdsFSpIkSdK42OmAleTXgfcCNwKbm+a9gCcleV1V/WuL9UmSJEnS2OjnHqx3Ac+vqo3djUn2AS4CntpCXZIkSZI0dvq5B2shsGmK9s3AQwcrR5IkSZLGVz89WKcCX0xyNnBL07YMOBI4pa3CJEmSJGnc9DPJxV8n+RiwGngunWdfbQJeUVVfa7k+SZIkSRob/fRgUVXXA9e3XIskSZIkjbV+7sHaoSQntPl+kiRJkjROWg1YwJdafj9JkiRJGhutBqyq+nib7ydJkiRJ46SfBw0vBI4FXgL8HFDArcD5wClV9cNWK5QkSZKkMdHPJBdnAluBE/jJ87D2Ao4GPgS8vJXKJEmSJGnM9BOw9q+qfSe1bQKuTPIfLdQkSZIkSWOpn3uwvpvkZUl+fGyShyR5OfDd9kqTJEmSpPHST8A6EngpcHuS/2h6rb4J/GazrW9J/ijJdUmuTXJWkt2S7JPkqiQ3JvlokocNcg5JkiRJmi07PUSwqjbS3GeV5LFAqupbgxaSZCnwRmC/qtqW5Bw6ge0w4MSqOjvJ++lMsPG+Qc8nSZIkSW0baJr2qvp2G+Gqy0JgUTNT4e7AbcAq4Lxm++nAi1s8nyRJkiS1pu0HDfetqjYD7wBuphOs7qTz4OKtVXV/s9smYOlUxydZk2R9kvVbtmwZRsmSJEmStJ2dDlhN71Lrkvw0sBrYh87ztX4KOHSKXWuq46vq5KpaWVUrlyxZMhslSpIkSdK0+glLVybZBPwL8C/NPVlteD7wjaraApDkn4BfAhYnWdj0Yu1F56HGkuaQA9ddyuat2/o+funiRVyxdlWLFUmSJI1GP5NcrEzyeDq9S3/fTE7xOeCTwGer6t4+a7kZOCDJ7sA24GBgPfAZOrMWnk3nYcbn9/n+kmbJ5q3b2Lju8L6PX772wharkSRJGp2+7sGqqv+qqvdX1Yvp9DJ9nE4P1L8l6esvpaq6is5kFl8GrmlqOxl4C/DmJBuAxwKn9PP+kiRJkjTbBr6fqqp+CFzavCamW+/3vd4GvG1S803As/suUJIkSZKGpPUJK5rZACVpl+D9Z5IkqduszAgoSbsK7z+TJEnd+n4OVpKnt1mIJEmSJI27QR40/P4kX0jyuiSLW6tIkiRJksZU3wGrqn4ZeAWwDFif5CNJXtBaZZIkSZI0ZgbpwaKqbgSOpzOV+vOAk5LckOQ32yhOkiRJksbJIPdgPSPJicD1wCrgRVX11Gb5xJbqkyRJkqSxMcgsgu8GPgi8tap+PEdxVd2a5PiBK5MkSZKkMTNIwDoM2FZVDwAkeQiwW1XdU1VntlKdJEmSJI2RQe7B+jSwqGt996ZNkiRJknZJgwSs3arq7omVZnn3wUuSJEmSpPE0SMD6fpL9J1aS/AKwbZr9JUmSJGleG+QerDcB5ya5tVnfE3j54CVJkiRJ0njqO2BV1ReTPAXYFwhwQ1X9sLXKJEmSJGnMDNKDBfCLwPLmfZ6VhKo6Y+CqJEmSJGkM9R2wkpwJPBH4KvBA01yAAUuSJEnSLmmQHqyVwH5VVW0VI0mSJEnjbJBZBK8FfratQiRJkiRp3A3Sg7UH8LUkXwDunWisqiMGrkqSJEmSxtAgAeuEtoqQJEmSpPlgkGnaP5vk8cCKqvp0kt2BBe2VJkmSJEnjpe97sJK8FjgP+EDTtBT4WBtFSZIkSdI4GmSSi9cDBwJ3AVTVjcDj2ihKkiRJksbRIAHr3qq6b2IlyUI6z8GSJEmSpF3SIAHrs0neCixK8gLgXODj7ZQlSZIkSeNnkIC1FtgCXAP8d+Ai4Pg2ipIkSZKkcTTILII/Aj7YvCRJkiRpl9d3wEryDaa456qqnjBQRdKYOnDdpWzeuq3v45cuXsQVa1e1WJEkSZKGbZAHDa/sWt4NeBnwmMHKkcbX5q3b2Lju8L6PX772wharkSRJ0ij0fQ9WVX2767W5qv4eGOi/35MsTnJekhuSXJ/kuUkek+TiJDc2X396kHNIkiRJ0mwZ5EHD+3e9Vib5PeCRA9bzLuBfquopwDOB6+lMpnFJVa0ALmnWJUmSJGnOGWSI4N91Ld8PbAR+q983S/Io4FeBYwCaZ2zdl2Q1cFCz2+nAZcBb+j2PJEmSJM2WQWYR/LU2CwGeQGfa939I8kzgS8AfAj9TVbc157wtyeNaPq8kSZIktWKQWQTfPN32qnpnH7XsD7yhqq5K8i52YjhgkjXAGoC99957J08tSZIkSYMb5EHDK4HfB5Y2r98D9qNzH1Y/92JtAjZV1VXN+nl0AtftSfYEaL7eMdXBVXVyVa2sqpVLlizp4/SSJEmSNJhB7sHaA9i/qr4HkOQE4Nyq+t1+3qyqvpnkliT7VtXXgYOBrzWvo4F1zdfzB6hZkiRJkmbNIAFrb+C+rvX7gOUDVQNvAD6c5GHATcBr6PSynZPkWOBmOs/bkiRJkqQ5Z5CAdSbwhST/DBTwEuCMQYqpqq+y/QOMJxw8yPtKkiRJ0jAMMovgXyX5JPArTdNrquor7ZQlSZIkSeNnkEkuAHYH7qqqdwGbkuzTQk2SJEmSNJb6DlhJ3kbngb/HNU0PBT7URlGSJEmSNI4G6cF6CXAE8H2AqrqV/qZnlyRJkqR5YZCAdV9VFZ0JLkjyU+2UJEmSJEnjaZCAdU6SDwCLk7wW+DTwwXbKkiRJkqTxM8gsgu9I8gLgLmBf4M+r6uLWKpMkSZK0yzpw3aVs3rqt7+OXLl7EFWtXtVhRb/oKWEkWAJ+qqucDhipJkiRJrdq8dRsb1x3e9/HL117YYjW962uIYFU9ANyT5NEt1yNJkiRJY6vvIYLAD4BrklxMM5MgQFW9ceCqJEmSJGkMDRKwLmxekiRJkiT6CFhJ9q6qm6vq9NkoSJIkSZLGVT/3YH1sYiHJP7ZYiyRJkiSNtX4CVrqWn9BWIZIkSZI07voJWLWDZUmSJEnapfUzycUzk9xFpydrUbNMs15V9ajWqpMkSZKkMbLTAauqFsxGIZIkSZI07vp60LAkSZIk6cEMWJIkSZLUEgOWJEmSJLXEgCVJkiRJLTFgSZIkSVJLDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEmSJElSSwxYkiRJktQSA5YkSZIktcSAJUmSJEktMWBJkiRJUkvmXMBKsiDJV5J8olnfJ8lVSW5M8tEkDxt1jZIkSZI0lTkXsIA/BK7vWn87cGJVrQC+Cxw7kqokSZIkaQZzKmAl2Qs4HPg/zXqAVcB5zS6nAy8eTXWSJEmSNL05FbCAvwf+FPhRs/5YYGtV3d+sbwKWTnVgkjVJ1idZv2XLltmvVJIkSZImmTMBK8kLgTuq6kvdzVPsWlMdX1UnV9XKqlq5ZMmSWalRkiRJkqazcNQFdDkQOCLJYcBuwKPo9GgtTrKw6cXaC7h1hDVqlh247lI2b93W17FLFy/iirWrWq5IkiRJ6t2cCVhVdRxwHECSg4A/qapXJDkXeClwNnA0cP7IitSs27x1GxvXHd7XscvXXthyNZIkSdLOmTNDBKfxFuDNSTbQuSfrlBHXI0mSJElTmjM9WN2q6jLgsmb5JuDZo6xHkiRJknoxDj1YkiRJkjQWDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEmSJElSSwxYkiRJktQSA5YkSZIktcSAJUmSJEktMWBJkiRJUksMWJIkSZLUkoWjLkBzy4HrLmXz1m19H7908SKuWLuqxYokSZKk8WHA0nY2b93GxnWH93388rUXtliNJEmSNF4cIihJkiRJLTFgSZIkSVJLDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEmSJElSSwxYkiRJktQSn4MlSZIkad5ZunjRSJ7RasCSJEmSNO9csXbVQMfn7f0d5xBBSZIkSWqJAUuSJEmSWmLAkiRJkqSWGLAkSZIkqSVOcjGFA9ddyuat20Z2/qWLFw18U54kSZKk4TNgTWHz1m1sXHf4yM4/iukkJUmSJA3OIYKSJEmS1BIDliRJkiS1xIAlSZIkSS2ZMwErybIkn0lyfZLrkvxh0/6YJBcnubH5+tOjrlWSJEmSpjJnAhZwP/DHVfVU4ADg9Un2A9YCl1TVCuCSZl2SJEmS5pw5E7Cq6raq+nKz/D3gemApsBo4vdntdODFo6lQkiRJkqY3ZwJWtyTLgWcBVwE/U1W3QSeEAY/bwTFrkqxPsn7Lli3DKlWSJEmSfmzOBawkjwD+EXhTVd3V63FVdXJVrayqlUuWLJm9AiVJkiRpB+bUg4aTPJROuPpwVf1T03x7kj2r6rYkewJ3jK7C4Vi6eFHfDxteungRV6xd1XJFkiRJknoxZwJWkgCnANdX1Tu7Nl0AHA2sa76eP4LyhmqQgNRvMJMkSZI0uDkTsIADgVcB1yT5atP2VjrB6pwkxwI3Ay8bUX2SJEmSNK05E7Cq6nNAdrD54GHWIkmSJEn9mDMBq003fPN7Aw2VW7p4UYvVSJIkSdpVzMuA9cMHfsTGdYePugxJkiRJu5g5N027JEmSJI0rA5YkSZIktcSAJUmSJEktMWBJkiRJUksMWJIkSZLUEgOWJEmSJLXEgCVJkiRJLTFgSZIkSVJLDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEmSJElSSwxYkiRJktQSA5YkSZIktcSAJUmSJEktMWBJkiRJUksMWJIkSZLUEgOWJEmSJLXEgCVJkiRJLTFgSZIkSVJLDFiSJEmS1BIDliRJkiS1xIAlSZIkSS0xYEmSJElSSwxYkiRJktQSA5YkSZIktcSAJUmSJEktMWBJkiRJUkvGJmAlOSTJ15NsSLJ21PVIkiT4M7VJAAAHD0lEQVRJ0mRjEbCSLADeAxwK7AcclWS/0VYlSZIkSdsbi4AFPBvYUFU3VdV9wNnA6hHXJEmSJEnbSVWNuoYZJXkpcEhV/W6z/irgOVX1B137rAHWNKtPB64deqHa1e0BfGvURWiX43WnYfOa0yh43WkU9q2qR+7sQQtno5JZkCnatkuGVXUycDJAkvVVtXIYhUkTvO40Cl53GjavOY2C151GIcn6fo4blyGCm4BlXet7AbeOqBZJkiRJmtK4BKwvAiuS7JPkYcCRwAUjrkmSJEmStjMWQwSr6v4kfwB8ClgAnFpV101zyMnDqUzajtedRsHrTsPmNadR8LrTKPR13Y3FJBeSJEmSNA7GZYigJEmSJM15BixJkiRJaslYB6wkhyT5epINSdZOsf3hST7abL8qyfLhV6n5pofr7s1Jvpbk6iSXJHn8KOrU/DHTNde130uTVBKnMtbAernukvxW8/vuuiQfGXaNmn96+Dd27ySfSfKV5t/Zw0ZRp+aPJKcmuSPJlM/QTcdJzTV5dZL9Z3rPsQ1YSRYA7wEOBfYDjkqy36TdjgW+W1VPAk4E3j7cKjXf9HjdfQVYWVXPAM4D/ma4VWo+6fGaI8kjgTcCVw23Qs1HvVx3SVYAxwEHVtXTgDcNvVDNKz3+vjseOKeqnkVnVun3DrdKzUOnAYdMs/1QYEXzWgO8b6Y3HNuABTwb2FBVN1XVfcDZwOpJ+6wGTm+WzwMOTjLVQ4ulXs143VXVZ6rqnmb1SjrPbZP61cvvOoC/oBPmfzDM4jRv9XLdvRZ4T1V9F6Cq7hhyjZp/ernuCnhUs/xofC6qBlRVlwPfmWaX1cAZ1XElsDjJntO95zgHrKXALV3rm5q2KfepqvuBO4HHDqU6zVe9XHfdjgU+OasVab6b8ZpL8ixgWVV9YpiFaV7r5Xfdk4EnJ7kiyZVJpvsfYKkXvVx3JwCvTLIJuAh4w3BK0y5sZ//2G4/nYO3AVD1Rk+ec72UfaWf0fE0leSWwEnjerFak+W7aay7JQ+gMgT5mWAVpl9DL77qF/3979xJqVR3Fcfz7Sy0RQ4g7zLpB2QMpaRDaC0szcuCkoByUSrNoYtGoKGhaTSJ6SWEUaNmgLlJYEZKIQiKhPUFSLJCSIIPe5mqwN2HXi/fo3Xo9x+8HDne/99qw4Jy113/vSzNkZiFNp35LkrlV9fMpjk2Dq5e8Ww6srapnkiwAXm/z7sipD09nqROuJ/q5g/U9MPuo+Qs5tk383zZJptK0ko/XApTG00vekWQx8CiwrKr+PE2xaTCNl3PnA3OBzUn2AfOBEV90oQnq9Tv23ar6u6r2At/QFFzSyeol7+4H3gKoqm3AdGDotESns1VPv/2O1s8F1qfAZUkuSXIuzYOOI6O2GQFWtNN3AR+X/1lZEzNu3rXDtV6iKa58JkETddycq6pDVTVUVcNVNUzz3N+yqtoxOeFqQPTyHfsOcAtAkiGaIYPfntYoNWh6ybv9wCKAJFfSFFgHT2uUOtuMAPe1bxOcDxyqqgPH26FvhwhW1eEkDwKbgCnAq1X1RZIngR1VNQK8QtM63kPTubpn8iLWIOgx754CZgIb2neq7K+qZZMWtPpajzkndarHvNsELEnyJfAP8EhV/TR5Uavf9Zh3DwNrkqymGaa10pvnmogk62iGOg+1z/Y9AUwDqKoXaZ71WwrsAX4DVo17THNSkiRJkrrRz0MEJUmSJOmMYoElSZIkSR2xwJIkSZKkjlhgSZIkSVJHLLAkSZIkqSMWWJKks16S4SSfT3YckqT+Z4ElSZIkSR2xwJIk9ZW22/R1kteS7ErydpIZo7Z5M8nSo+bXJrmz3XdLkp3t5/oxjr8yyXNHzW9MsrCdXpJkW7vvhiQzT+GlSpL6kAWWJKkfXQ68XFVXA78AD4xavx64GyDJucAi4D3gR+C2qrq2Xf9srydMMgQ8Bixu998BPDTB65AkDRgLLElSP/quqra2028AN45a/z5wa5LzgDuAT6rqd2AasCbJbmADcNUJnHN+u/3WJJ8BK4CLJ3ANkqQBNHWyA5Ak6STUqPlZbdED8HhVjSTZDNxO06la165bDfwAXENzk/GPMY59mP/fgJze/g3wYVUtn3j4kqRBZQdLktSPLkqyoJ1eDmysqnntZ6Rdvh5YBdwEbGqXzQIOVNUR4F5gyhjH3gfMS3JOktnAde3y7cANSS4FSDIjyZyuL0yS1N8ssCRJ/egrYEWSXcAFwAtjbPMBcDPwUVX91S57vt1vOzAH+HWM/bYCe4HdwNPAToCqOgisBNa1590OXNHVBUmSBkOqRo+ykCTpzJVkmKZjNXeSQ5Ek6Rh2sCRJkiSpI3awJEmSJKkjdrAkSZIkqSMWWJIkSZLUEQssSZIkSeqIBZYkSZIkdcQCS5IkSZI68i/EA032Un6iOwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "N_bins = 50\n", "\n", "if (N_exp > 1):\n", " fig2, ax2 = plt.subplots(nrows=3, figsize=(12, 14))\n", " \n", " ax2[0].hist(all_p_mean, N_bins, (0, 1), histtype='step')\n", " ax2[0].set(title='Histogram, probability mu', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n", "\n", " ax2[1].hist(all_p_chi2, N_bins, (0, 1), histtype='step')\n", " ax2[1].set(title='Histogram, probability chi2', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n", " \n", " ax2[2].hist(all_p_ks, N_bins, (0, 1), histtype='step')\n", " ax2[2].set(title='Histogram, probability Kolmogorov', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n", "\n", " fig2.tight_layout()\n", "\n", "\n", " if save_plots:\n", " fig2.savefig('PvalueDists.pdf', dpi=600)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Questions:\n", "\n", "1. First run the program with one experiment (N_exp = 1) to display the two distributions A and B, when:\n", " - They are the same.\n", " - The mean of A is increased (to e.g. 0.1).\n", " - The width of A is enlarged (to e.g. 1.2).\n", " - The normalisation of A is increased.\n", "\n", "Get a feel for how much you need to change the distribution, before you can _by eye_ see that they are not the same. I.e. can you see any difference, if `mean_A` $= 0.1$? Or how about $0.2$? How do you quantify this and when do you start to doubt? And how about `width_A` $= 1.1$? Or $1.2$? Again, can you see it by eye? Finally, try to put $1050$ events into B. Is that visible? How about $1100$? And do you see an impact in the p-values?\n", "\n", "2. Could you for the test of the means have calculated how much of a change in the mean is needed for a difference to be statistically significant? Do so, and see if it somewhat matches you guess/estimate from above!\n", "\n", "\n", "3. Now run the tests 1000 times, where A and B are unit Gaussians and thus identical. How should the distributions of the test probabilities come out? And is this the case, approximately? If not, think of reasons for this, and what could be a remedy. HINT: Large statistics is always easier!\n", "\n", "\n", "4. Repeat the changes in question 1), and see which tests \"reacts\" most to these modifications. How much of a change in the mean is required for 95% of the tests (of each kind) to give a probability below 5%? How much is required for the width? And the norm?\n", "\n", "\n", "5. Possibly try to test different distributions than the Gaussian one (e.g. exponential, uniform, etc.), and see how the tests performs.\n", "\n", "\n", "NOTE: The Kolmogorov-Smirnov test has the great advantage that it can handle ANY distribution (even the Cauchy distribution - remind yourself of that one!). The reason is, that it doesn't care about any PDF, nor how far out an outlier is. It is just a matter of the difference in integrals between the two functions.\n", "\n", "\n", "## Advanced:\n", "\n", "6. Obviously, the test of the means is not sensitive the a change in the width. Make such a test yourself by calculating the widths and the uncertainty on the widths (or perhaps try the F-test!). Note that in a (unit) Gaussian the uncertainty on the width is of the same order as that of the means!\n", "\n", "\n", "## Very advanced:\n", "7. Implement in python the following tests:\n", " - Lilliefors test\n", " - Shapiro-Wilk test\n", " - Anderson-Darling test\n", " - Cramer-von-Mises test\n", " - Jarque-Bera test\n", " - Kuiper's test\n", " - Mann-Whitney-Wilcoxon test\n", " - Siegel-Tukey test\n", " \n", "and quantify under various conditions and datasets the power of each and the correlation among them. Write it up, and send it to a statistics journal. :-)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }