{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to hypothesis testing\n", "> A Summary of lecture \"Statistical Thinking in Python (Part 2)\", via datacamp\n", "\n", "- toc: true \n", "- badges: true\n", "- comments: true\n", "- author: Chanseok Kang\n", "- categories: [Python, Datacamp, Data_Science, Statistics]\n", "- image: images/frog-swarmplot.png" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "sns.set()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def ecdf(data):\n", " \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n", " # Number of data points: n\n", " n = len(data)\n", "\n", " # x-data for the ECDF: x\n", " x = np.sort(data)\n", "\n", " # y-data for the ECDF: y\n", " y = np.arange(1, n + 1) / n\n", "\n", " return x, y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Formulating and simulating a hypothesis\n", "- Hypothesis testing\n", " - Assesment of how reasonable the observed data are assuming a hypothesis is true\n", "- Null hypothesis ($H_0$)\n", " - Another name for the hypothesis you are testing\n", "- Permutation\n", " - Random reordering of entries in an array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generating a permutation sample\n", "Permutation sampling is a great way to simulate the hypothesis that two variables have identical probability distributions. This is often a hypothesis you want to test, so in this exercise, you will write a function to generate a permutation sample from two data sets.\n", "\n", "Remember, a permutation sample of two arrays having respectively ```n1``` and ```n2``` entries is constructed by concatenating the arrays together, scrambling the contents of the concatenated array, and then taking the first ```n1``` entries as the permutation sample of the first array and the last ```n2``` entries as the permutation sample of the second array." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def permutation_sample(data1, data2):\n", " \"\"\"Generate a permutation sample from two data sets.\"\"\"\n", " \n", " # Concatenate the data sets: data\n", " data = np.concatenate((data1, data2))\n", " \n", " # Permute the concatenated array: permuted_data\n", " permuted_data = np.random.permutation(data)\n", " \n", " # Split the permuted array into two: perm_sample_1, perm_sample_2\n", " perm_sample_1 = permuted_data[:len(data1)]\n", " perm_sample_2 = permuted_data[len(data1):]\n", " \n", " return perm_sample_1, perm_sample_2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualizing permutation sampling\n", "To help see how permutation sampling works, in this exercise you will generate permutation samples and look at them graphically.\n", "\n", "We will use the Sheffield Weather Station data again, this time considering the monthly rainfall in June (a dry month) and November (a wet month). We expect these might be differently distributed, so we will take permutation samples to see how their ECDFs would look if they were identically distributed." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "rain_june = np.array([ 66.2, 39.7, 76.4, 26.5, 11.2, 61.8, 6.1, 48.4, 89.2,\n", " 104. , 34. , 60.6, 57.1, 79.1, 90.9, 32.3, 63.8, 78.2,\n", " 27.5, 43.4, 30.1, 17.3, 77.5, 44.9, 92.2, 39.6, 79.4,\n", " 66.1, 53.5, 98.5, 20.8, 55.5, 39.6, 56. , 65.1, 14.8,\n", " 13.2, 88.1, 8.4, 32.1, 19.6, 40.4, 2.2, 77.5, 105.4,\n", " 77.2, 38. , 27.1, 111.8, 17.2, 26.7, 23.3, 77.2, 87.2,\n", " 27.7, 50.6, 60.3, 15.1, 6. , 29.4, 39.3, 56.3, 80.4,\n", " 85.3, 68.4, 72.5, 13.3, 28.4, 14.7, 37.4, 49.5, 57.2,\n", " 85.9, 82.1, 31.8, 126.6, 30.7, 41.4, 33.9, 13.5, 99.1,\n", " 70.2, 91.8, 61.3, 13.7, 54.9, 62.5, 24.2, 69.4, 83.1,\n", " 44. , 48.5, 11.9, 16.6, 66.4, 90. , 34.9, 132.8, 33.4,\n", " 225. , 7.6, 40.9, 76.5, 48. , 140. , 55.9, 54.1, 46.4,\n", " 68.6, 52.2, 108.3, 14.6, 11.3, 29.8, 130.9, 152.4, 61. ,\n", " 46.6, 43.9, 30.9, 111.1, 68.5, 42.2, 9.8, 285.6, 56.7,\n", " 168.2, 41.2, 47.8, 166.6, 37.8, 45.4, 43.2])\n", "\n", "rain_november = np.array([ 83.6, 30.9, 62.2, 37. , 41. , 160.2, 18.2, 122.4, 71.3,\n", " 44.2, 49.1, 37.6, 114.5, 28.8, 82.5, 71.9, 50.7, 67.7,\n", " 112. , 63.6, 42.8, 57.2, 99.1, 86.4, 84.4, 38.1, 17.7,\n", " 102.2, 101.3, 58. , 82. , 101.4, 81.4, 100.1, 54.6, 39.6,\n", " 57.5, 29.2, 48.8, 37.3, 115.4, 55.6, 62. , 95. , 84.2,\n", " 118.1, 153.2, 83.4, 104.7, 59. , 46.4, 50. , 147.6, 76.8,\n", " 59.9, 101.8, 136.6, 173. , 92.5, 37. , 59.8, 142.1, 9.9,\n", " 158.2, 72.6, 28. , 112.9, 119.3, 199.2, 50.7, 44. , 170.7,\n", " 67.2, 21.4, 61.3, 15.6, 106. , 116.2, 42.3, 38.5, 132.5,\n", " 40.8, 147.5, 93.9, 71.4, 87.3, 163.7, 141.4, 62.6, 84.9,\n", " 28.8, 121.1, 28.6, 32.4, 112. , 50. , 96.9, 81.8, 70.4,\n", " 117.5, 41.2, 124.9, 78.2, 93. , 53.5, 50.5, 42.6, 47.9,\n", " 73.1, 129.1, 56.9, 103.3, 60.5, 134.3, 93.1, 49.5, 48.2,\n", " 167.9, 27. , 111.1, 55.4, 36.2, 57.4, 66.8, 58.3, 60. ,\n", " 161.6, 112.7, 37.4, 110.6, 56.6, 95.8, 126.8])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZydVZno+99a77jHql1zpTInEDAkgGQiSBC9DRwFvU4gIthtN6ft092c5tq0Q+vpP0Tl49D29XTftrW79TpAM9mKfU/beERFTAhJwBhopkwklUrN057eca37x65UKAJRhEpCZX0/n/qk3ryV2uvdqdrPXs961/MIrbXGMAzDOO3Jkz0AwzAM49RgAoJhGIYBmIBgGIZhTDEBwTAMwwBMQDAMwzCmmIBgGIZhACcgIFQqFa688kp6e3uPOffkk0/yzne+k8svv5y//Mu/JEmS2R6OYRiG8RLs2fzmO3fu5BOf+AT79+9/0fO33HILt956K+eddx4f//jHueuuu3jf+973G3//sbEqSjW2UbS25hkZqbwawz5lmWucG+b6Nc7164PX7jVKKSiVci95flYDwl133cVf/dVf8Rd/8RfHnDt06BBBEHDeeecB8M53vpMvf/nLLysgKKWnA8KR47nOXOPc8IqvUWvQmiRufNiuxLJApwphScJYEgTg+42PE+GhB+HnD0re9raUlatP3P9hbbhGbbhGti1Lti17wh73VP45rX37+/DD++GKy8i+/+2/8b+b1YDw6U9/+iXPDQ4O0t7ePn3c3t7OwMDAbA7HMOYGrSEMSSJF39OTaNdDa01bC9huIxDseq6IsBoZ4TVrZj8oPPQgvPfaHEki+H++AnfcDm/YNLuPCY1g8IuvbEYAGrjoQxtPaFA4FdW+/X3m/V+/i4Uivf92+vjGbxwUZjUgHI9SCiHE9LHWesbxb6K1NT/juL298KqM7VRmrnFuaG8vHHmTjxCNj2Mc+YIjhwg0AqEVqu5QLSs86riOTRxpfMtGuh5JmLJvt8+zuyUrzhJceqlNc/PsXs+jOxKSRKC0IE40j+7I8I53zf7Ly/BISM/eXSyYeJqDTSvIiYtpO0E/P6fqz+nIz3+GhUKiAUX+5z+j9eb3/0b/9qQFhK6uLoaGhqaPh4eH6ejoeFnfY2SkMj1ta28vMDRUflXHeKox1zg3tLcXGBgoE9SnXuKlwPPF0aCgNSiFDkK00iSVgFi4JEpi5Xy0hrG9Y4weLNP/t//CwtoeDrefw46O/8KewRKenXLPQ82kSmJZiu7uSd582ct7s/VyrVwWYIlmtJZYKFYuG2doaPZzVbV//yEb/+1TWKQswKLvTUvRrb95iuS3dSr/nNYuvoTiv34dUKRIKhdfgpoaq5TimDfSz3fSAkJPTw+e57Fjxw4uuOACvv/977Np0wmYYxrGSaY1BHVNVA6xLLBt0K6HsMR0OkgnKeFwmVi4PPfLIUShicnRkEWr26hVFPu2DuJv+wXX7PwCFgm/eLaPG8RHSLARaJQWgCBNBb94IOXNl83ur/o5bYe4edGd7NVnc5bzLOe0bQCWzepjAuR2/AKLdOrdcEpuxy/QLyNnPhdl3/92+vjGqbeG8GJuvPFGbrrpJlatWsUXvvAFPvGJT1CpVFi5ciU33HDDiR6OYbw6tEaljRdiaQmknHFqRmpIa0ApSBKiWCJ9EGigcTKONGEgCCcigiAiqgSUWjMMVQPKoxFRLaVp905e9+jtbGMND3IJ+1lIqi00Ej2VLtAoLKlYe04FmN2ckS0VK5p7ObctRFQmsKWa1cc7In3zm9B3/CNKpWhpkb75TWZzFY2gwG8RGE9IQHjggQemP//a1742/flZZ53FPffccyKGYBivCq00WjXSPEJOpWG0Jq2FTExMvej7HsWmRlDQGuq1RrCQliCTbbzoJz9+EPf+/6B69vmo/2MjlOYBEMfw1KNVJvprHHp4N60dkoln+5l8XTc7Nmv2281syOzkj3/0p2xNX89l/G8iXCwSJI0XYUso3rbgIWp2jrPbD/OGN14w689L8cwuFq+bT1wN6ThvIcUzu2b9MQHklVcw/I+3Y/34gUYwuPKKE/K4c9VJSxkZxmuNVppwMpw+9ooeQgq00tTrEKUWrpWC0ijVCAhpopkcaqSG0hTcHo/4pw+x+I+vQSQRyrJ5NvkimQXXYOctymMJT+8YIh2vcWhnL7lzm4mHJtj1zBL+esslpNrieyxkrb6An3MJES4pNhp4Q8cvWbTMZom/nwvWQ2lpljBoxfKdWX9unJYCZ938TuLRMvNWzGM8nfWHnCavvAJ95RVmZvAqMAHBMI5jOt2DRicpOlVIx0IlCq009vZHsH7+c6JV60mXraGSaLI20ykjrTRRLQatQEi0cuGnP0MkEVIrSBO8Xz6GTt8DQBqlpEGMEookiAjCkH2jHfzk0FkkqpESirH5iXgTb9Q/wSUiRGChuGjh02za5FAZrOC3LEb4Nr5n42RnPyBAIygc+eAUXXA1js8EBMN4CVpD+uAWvM0PEl6wAb1+LfXBCbBdhCXJPr2T5mvejogi8o7LwU/8C6PLX88CS0NrFhCkYcKen+ymXovIZF065q+ESy5B3eagk5gUi12cwbmJxAEyvsYZ7iMemaBeKXLnw+fzi75zSLQFCMTUuoC9sgtZa+VP99zK451v4pyOAd5yfRdOzqPjfJ/28xeiUo2TdU5YQDBe+0xAMIyXILdsoe39b0fEEcp2Ofy1uwm7zyTjS7QG6+c/R0QRQqWIOMJ68MeMZpZjVVJKnT34BYfyQI3Duw7jNecY3zNCeWAJTedcwH9c/Emadm3hcPFMRp1FLBuLybX6yChgSVuFJ+qdfG3ocmJlo2ncMSRIOaezl6uWP8Lqi7qJ2v+U8x47yH95vcDJLWT+RUuxsy6WZ2N55lfbePnMT41x2nrRjWFpCmmKtWMHmS/c1njB1wqZxGQf/hlDVy4njAWW1ATr3oByXEQckVoO+9tXE04ExK5GJY0F3jTRBBMJkaigJlLSRBPXE4bblvJUVx6VaLJxim019tNICU8PdvKdHedOBQNJY0uawiblyuXbOKtjhFzLQvz2HNmuIpnWLNJ3cYoZMxswXhETEIzT0tTt/tM8D4RKYWgI67HHaP6vH0BEEWiFFhJl2Uyes47RgRA79JCOhX/uGrb++Tco7nqYHeFZ9I3msCd307zCRbIAAGkJJgYnSfo1tmzckupnJZl4nKZ6P0qlzO9oIt/c+FV87Nkit9xzOVHc2JUsUUhS3rz8SV5f+BUbL86TbT+TeWt6EFKSn9+G11rAznnYGRMMjFfGBATjtHJkVnCkIoRlgdi8BW/rz0kv3EC87EycRx9FxI2ZgRKSifMvYv/b/wjOXIM3VKHQYlFPHJIwRZ93LkOrXk/lp/vpiMZwsgVKHbIx0wCSIKGjp4DTmSMeqJIECY4DZ27sYavu5vHD3ZQ6FJbVGM+WR3ziVDZKVKBYf84Yv3/hNta+uUTgnEfTvCLassl1ZrAsKEgJUr50+QvDeBlMQDBe0+xtW3E2P0S88Q0ka9fPOPeCUkAAhEFjHwFCNNI6W35Bzx+8G5FEaMflwOf/ieqScynYLiKJSaXNT868jr7HNEvUE/Tvq9N63lKyRZe2jk4OPjlCZajMnp8dJh+NYbmQX5NFvG1VY3w27Nol2PdwiSW5Gr9jgxaSbY85fPQHbyZRFrc/prln/STrNsGFGyIcyyVGYIuUP7h0F2d1V7DblyICAZ6HZUss10JaJgIYry4TEIzXLHvbVprf9TaII3Bcxu+9bzooHEkJHfnT80CljXIRtg1RPaU+PEnHf9wPU7MB4hhn+3ZGrv8jHvpvX6F9z2M8NNjNqFpAvGc/i3o0EyPgxx5p7FIfC5AT/Xhhih2OYdsROBnCKCatRwA8tl3zrb6rSbTFQxMpl24/yLJzJU/XV5BoG4UkThWbt9is2wRrzov4h0/tZcfTzZzTeYiLLy2RuvPIziuRTRROzsX2TDAwZofZy2G8Zvl33QFhgEhTiCOczQ8dqQt3zIeUkMaK6lhIdSQgnAjQqaK+/iK0ZaNlY51g/Jz1VCKX/vnnsuvSP+Q5fwlZOYYKYoZqNmmQ4lXHUFFMGiu0AuVmp8pWKHwXbNuZnpk8+kSOZKqkRKItHn0ih+3bXLCygmMpLKFwbM3GDTEAwrF5/Tl1/vCdvZx/bow1rwsrlwGlsByJ45tgYMweM0MwXlOOpIhUqQX/jm+D1mgAyybc8AbqdRqlJBLNxHgjPRQEGlVUDO6vMrjjWQSaUouFXSxQ9hfy3O98mAW1/VTPPp/B7FLKTw7Ru6tMoSNDvW+U0ckUOTlAbgT6hxwOPjuBNwor1xTR9QQ1MoGvInqr3RyaXIJ7dsRFzY2a/J0LnanbRjUaQedCB8uRXPrWLP8of8WOp5vYdCmsvagVAMu1KC5uR4UxOZXi5wQi74Hb2PswXS7DMGaBCQjGa8aMFJEQoFTjpVYIKu+5juEz1kO9UVdIJYrKgUkyBZvJ3SOotiwHdw7DRA1RKFCfCGltFQgt2J1ZyujrNlAdrVJ/NsApFQjLCQu6Jyh5MdkWh2hUQSbFVRGlXISq9BOMtFPoaSbfkuFXz+b41q53kqSSn/yLZvn/Ocn6i2G83Kg+qpFIoRgv2+hU4eQ9Ln1fFxeX63gdTQi7saos0GQKNrrJQ6QJwnWmVr5NIDBmnwkIxinvyKzA6j3YyPenKXrq7hotBNpxqbz9msZisdbU6xoLTRpEhCohqkVoWSBNIapH2EkZ3WkhXYd61SaNNU4aktQ0SRaaXIVOFeMVAXFIVtdRWmHZWWxrEjupo4SFbYP0XZTrs314KVHa2E2cTK0JrL8YNl6Y4DqaJFXYlmbjhcl0JzMVp0jfRXrPu11UNPoiCJ2CFCYYGCeUCQjGKW3GrMCywGoUcsNxKX/qNsTAAIPL1lLtOpv6YIDlu4z3h8g4pPLkXrycw6FHh5gsK2ojIRnfx8tAplSkGtvU6hG9Yx7V/TFoSVPBIa2HtIgh2iPB4Fg/k0mGdLzxPZuLkCtKsm0FmhY2UR0PeHSny8OHV0yNWGNZcNHFjXLW6y4UfOO23Wx9LMv682usu7ADYUm89uJ0/+MjAQJovPh73q9ppWYYs8MEBOOUZW/bSvbzn4UoRCiFBurv/wBq/oLGbaYXrCWerBMORfhSodOUKFVkc4JkIsUr+GSaChQWKHpWtNLn1cnIgPaeLNK1SP0cTb5DxxnddC/wSdOInpV5bFcwbjXhq5TWRSVaVnZRfXaAntXzGWtv48HxM9h4ZoCddSmds5A9D1koLQGBEJpr3hOw7sKjL+Qb3+TzhiskKjraQeyYQPB8JhAYJ4kJCMYp4UhaiLdeDsvPOToziMJGO0kpwXEJr752+tbSNEyYPDTBoZ2DYNkEbjOhzDCycx9OVCV4Zi9+d4HaSILs7KH3iQnccIzenZL5q9qQ2YBgLGBs29NwUGJ5Du0liFTKwa37sT2H8b0BSTDAc2Ot3HfPSv7XrjNIteQrd2ru7K6wuCfg3FUpjq1JEnAczTvemR59gz/1wv+iswHDOMWYgGCcdDPSQn/9Oex77msEhzhqzAykJNr0Rmq3fGw6GKhEMbp3jAPbD9G7N8ApFJiUIUU5jO7rRYmAcn9K2F5EiCoFu057bpLIdohiRWUypUWUKeiANnWIrpYWgvGQ8r4MdnMTlbJDvqcV0R5Qmd/F/3ziXcS9U+kqBHECWx62WfGukAsvEPzzFw7y2O4SF1+UcP4aC63ldEB4yfSQYZxiTEAwTrrp/QRaQxRN7zzGcdE0Np09PxhAIyCE5YigqlFhQuRGhElIlATUgxhXKpRS+G5CJCR+IYNWmko5wHMFtgWkCicNECol1TZCNzaoWZ6HdAUFFbFLLuCBXRuJlTV9+6gQGseBjRsTpATLEaw5L2Ljm2pYrgXCnpHxMYHAeK0wAcE4IV6qxIS9beuM/QTCtqe/ZvzexkwhvvAiknPPhyhCyandvaHiue297HtoD7t3a3I9RSyvji0HCPf1Ipqy1GKBGE3R5YhwYBAvGCcXBDjKoiB8rCCmVq6TBBGVwREcN8fjQ/N4YlcXnUpRVSW+8vQF07WFpNTYluaaayKueW/K2jUQjHuNRjl+Ea/oNeoKSWGWAIzXJBMQjFl3vBITzuaHIE2n9xPwe783fS5Zu57kgrUwPg59fSgFkxShuZmJgZiRwwHVCoTVmOzACLE1gl6YJ9fkI3s6aW7VNPeUqAw/Q7T3KcLBiNYFLdiFPFldg1IzdrEAxTYKHT6P9/Xw2Z9eSpxKbHk2lzljU13KBFJoLtoQ8ZEPV1l3kTXVEk3gN/vH9lg2jNcoM481Zt2LlZg4Yjo1ZFng+XDDDRypP6GVRsUpOoxQCqJEkkYJNglaQ1DX1CoRajLE0ikiTNBpQKIdYiFxpMKPR5FoaF/S6CgwMYmUGpnPIADGx9lXX8g9uzfx4KFziBOBUoIkFWRzEsfWWFLjupq/uLnKurXpjDuAhBRI2+wgNuYGM0MwZtULU0JYjZTQETNSQxvfQGnDBugdbhSlCzQkCerwCGH/GMpyGBGt5K0caZBQqTr096Z44Si1XovW7ixevpPEydFc8hjfV6Y+FpJWqsRDA1i2RLS2oKVNZl4bKtE8uq/I/9xzAamykFJh25CmGtuGq65KeMtVo/zns1ku3pSyfn1jfcDkg4y5ygQEY1b5d90BcTydEqpfe90xZaqTteuP/p1u1B9KtUSFESIKCWWOuOTiZm0yaZZcQaIjRabTo2d5jvrBHIUzOliwMEvLqsVElYBiRuMHI5TOmMdIU4auC5YyWrVp7nBxwipN84o89iub27d0TzevB7jmbRUWL5NsWBNy/joHy7V4i51OjdRMqI25zQQEY9Z43/w6/rf/36Ozg6l9BC9qqnmBSjW14RppLWRs/whRoFCD/UzEOfLtOSJXY0moD9Xo+9EvmXy2H0E/lYEq/rrFVONDRKOTtLp16nueo9I7Si3M451foG/vKE9sr/FkfwdnrvL57D8vJorE9BqBY2uueXfIhouOXIBjZgPGacUEBGNW2Nu2Uvjoh2csGL/Y7ACYblqgNdRHakweGieYCBg/MIb0MlQnbSbKMbUkQKV1CsUuhp4eIKfK1FHUyKFwGasKMuPDOFFMNQ6wpMSd10y5V6GrNQarTXzuvvOIE4l8WKPU0WBw0doqt/zJBBsuyTYWjM1uYeM0ZAKCMSv8u+6AJJm6cx+wrKOzgxd2t59qHqClRRqlxLWIILZIwxRVHqE2EoPIkOIQhxG2pbEcgUwTQBEjcYRGK0HdcWmuTSJd0L6PsCwsT1AsOWx+qpsobtxCigJpgdCNmcFHPlxj3fmqMR5pUkPG6ckEBONVN50q4mgwKN/2RZI16xq9hqPo6Ltvz5v+PA0TDvemHNw1SN/+Ovn6EBlLEYzXGYnyMCZocgMm91tkrQhX+hyiToRNqnO0Do9TaBOIKCRbyhG5LeDlOGy18MB3F/HTnW1TO401tgOfvrXG2Ihg49o669ZoEJ4JBsZpzQQE41X1oqmi93+A8PrfbfSyVKoRELLZxudaN8pYux6V8YiB3cMMHYyoH+wnr4ZIZYrULo6XxXI1aaaZmlsiHO1nwrJQtGMxgMRDME5zaw6aFtPS5qJsh155Bn99+0qiRE5NRBoF6K69JuQDv6sag9YOaNukiYzTnnk7ZLyqXjJVdKSXpdVoBEOSNP6cegFOUkFQU9THA8LyJGktIk0kcRwhkwAbSTm0SbFpKgKWg5aNCqgaF00NJyOwWztwbHB9ye6Jefzzvy2cCgZi6uE0nqe5+prk6KCPpIlMMDBOc2aGYLxqjpsqqtchCBovuq6Lth2UsEAJdKro2xtwaEcffQ88Rd8je2mpHCZp1mSboZ7vZGK4DgXJeOTTO2Sj/BakyiM5SJ2IJgKyXidpcxtOUudpdTa3/ONZREnjRV7Kxt6Ca68JufqahLXr9Ml7ogzjFGUCgvGqeMlU0Q2/d3TdwLYhjht9jqspYZiC45COjDG4ZT/R0wcY392PE9cJs80EnkYWM2SzFl7NolioMJLpoKmrSH2oTFTI0oRNkYRWYljSQnFJN015xV13d08Fg8Zc5bzVCZ+6NWLtOnVynyjDOIWZlJHxqvi1dxUdKUWRKNJYkaQCkgQZBah6QBLBZA2qlRQ7rRPXAixbg22jJGjbI9ESkabYNmgBtpTkSJEoLK+JXCZDrsnh8YMlntrjzxjfqtXKBAPD+DXMDMF4xabLU/CCVNHa9Y1gEMeNrQajFbSCehgTBGWSWoSTtXGG+yg/1Uu4f4zhIYEXZxFocnEWO3Wp1yKEI4gSF68ksVVAIZNgS4kkwQFaWmD+msX0jjVz863ziKKjawaOw8w1A8MwXpQJCMYrMt3m8sjs4PmpIpi+k0h7PlpEaFuShGDrOsKTOJ6FzPh4i7pxXZ8lyTiJlcEGci0u85e41Oxmim4WuucjLEGp3SWUIctWZRlpXU5aUXSe2cLS9d188RtthJFA60a56k2bUm65JWTtWjM7MIxfZ1YDwg9+8AP+/u//niRJ+MAHPsB111034/wTTzzB//gf/4M4junu7ubzn/88xWJxNodkvErsbVvx7rqDzB3fgSQ+2ubS9Y6mipSCeh1dD1CTVSrPjVIfCyiPBIyXIYlSMq05CgQ8t32cwWcHmdhXo8AYQngEzQrR56H8MYaal2KPx/hxmSY1ClFMfTxG7h4kTjNEHRnu/Y9mbr/Tm7q9tLGIbIKBYfzmZm0NYWBggC996UvcfvvtfO973+POO+9k9+7dM77m05/+NDfddBP33XcfS5Ys4Z/+6Z9mazjGq+hIf4PMN78OUYhQCqbaXE73OtC6EQzKFcLJkPqBAUYPVQl7DzPWN8rkoQkGByT9fbD7oCYINZ6VkJEBvkjJZqtgJYRugbE0D1KQteuUvBr5Fh/fk/jzW5Gv68FZOZ+h7Jl89m9aSFNo7DWAa6+NTTAwjJdh1gLC5s2b2bBhA83NzWSzWS6//HJ++MMfzvgapRTVahWAer2O7/sv9q2MU8h0iigKGy0vmWps43qNNpdr1jVmC6lCpRqlBaiUNEkJKxHRWJ1aFVSQoJOYVEvCqqIWgECQtRTCUihpkaaCNE5wiMhlLJSTQUsa1VMluAWfXFORvrCLuzavIFVH7yqyLLj66vhkPlWG8ZozaymjwcFB2tvbp487Ojr41a9+NeNrPvrRj/LBD36Qz3zmM2QyGe66666X9RitrfkZx+3thd9+wK8RJ/Uav/pV+JM/aWwqm9phLGwbPvhBuOEGSuvXN9JDWhAFCuHGKM8lnKgxeLDO+C8PUhmeQDkF0mILOgzwMyFhLcSrVhkZFThxDRBYZY2d80E7OJk8+YymqxvaW9roWtKKdCzCJ0O29rXyf//4dcTpkUVksCzB3/0dvOUtuZP3XP0ac/1nda5fH8zNa5y1gKCUQjxv56fWesZxEAT85V/+Jd/4xjdYvXo1X//61/nIRz7CV7/61d/4MUZGKijVeJfa3l5gaKj86l3AKehkXqO9bSvNf/zHMxaPo01vbMwKplJE+sAQYTlCIYkqERkRUB2rM/bkYXY/PUF0aJKhaoa8mGAigIJQpIWEiXGFKuWx6gpPtTJaSbEcSbbJpbnDpXVBgtNUo7UwxrDqIAqzxJMKr1Pz2OC8qWDQmBmcf37Cpz7VWDcYGjopT9WvNdd/Vuf69cFr9xqlFMe8kZ5xfrYeuKuri6Hn/UYODQ3R0dExffzMM8/geR6rV68G4JprruGRRx6ZreEYr9CM3scAlnU0GACooykiKwnRUURYT4gm64AmTiRhvUIURwTKhlThegLhOCRaENcCHJ0SCQudCrI5hWU7BMLGdiDf7JI6LlpI8gWB7dsI22J4bOaP8KpVyqwbGMZvadYCwsaNG9myZQujo6PU63Xuv/9+Nm3aNH1+0aJF9Pf3s3fvXgB+/OMfs2rVqtkajvEKqVLL0UY3QO1DfzKzy1kYEVUiovEqtbE6enQEvW8vE/95kIG94wT7R+kbcaiP1xkd0+SDOlQrJNUytC4gIE8tkEQ1jaUj6jpHvihYuKKV5s4CIpfD8R3c9hKR8Mm2ZBiOWtiyLTs1osZdRWbdwDB+e7OWMurs7OTmm2/mhhtuII5j3v3ud7N69WpuvPFGbrrpJlatWsVnP/tZ/uzP/gytNa2trXzmM5+ZreEYr5CzaycwlZiREt3UdLSvgW7MDLTn42dCgqps7Eq2XGRrjlJRMTyakM8WycVVYtuitLAVv5jBz7tYVitNpZT6QEw8kiWXqZNpa6J1QZHO83rwix7V1KdlfoaCl8PPWfg+3PcDZ6pGXqOC6XXXmbuKDOOVmNV9CFdddRVXXXXVjL/72te+Nv35JZdcwiWXXDKbQzBeBccUrbMd4gsvahSsE6KxC7mmqAxWqe8+RHSoj9qhUQhDhpIiwslwaNthxveOIRSoQhN6sEaUyVNc3EqYhlT3DzCxb4Jy5NDaEpKpZCm1VhnsC2iVDqmjSYSD7UqyU03NJiYaWx1AozWsWpUe5yoMw/h1zE5l47hetGjde99Hcs7qRsE6y0JhUa+kTPSWGXmin3r/MLJSIck2E9g5xPg4uhKRtSLGlYVOFZOhIpMT2FmXTHUEippyS5Ym4eO1ZSm2uTTPb0G0F8gv6iCXF3hNHrYnkBK2bZN85Svu1CgbM4SxMVOayzBeCRMQjON6scXk6Z3IQqDjBCUat5kmtTp2FBDFgrSmkHZCSkqa6EbvgkRgoRE6ILFKKOmigwCpYzQay7WxLYHr+3jNGbTn4fo2jiuwXAvHk4ip1/zNm+3pTWhH9h1s3GjqFRnGK2ECgnFcxywm/+Efk1ywFspltNLUa5rYkpQf383gjgNMPraXidgiCjS0CATjyDBG9Q0wRkKEhVUvkB+v4LdY+JaDm0jGaxZhonC72uheFrH8rDx+d5HMohJ+s8/jTzpsedhm48aEtWsVpZKaLlEB8KEPRWb9wDBeIRMQjOOasZgMyInxRitM1yXVFpNjCcngIEM79lA7NMxkJaBuFbAl6LAGsWZS5Yj9AlZk0S4VmXbI9sxYlH8AACAASURBVPgUC5Bp9pgYSMkv7qRz/gT+8oWsWA3Ni0rkmmyUl2HXfzq85z1Z4hgcx+Xee2uMjUmkBKUaReyamk7q02QYc4JJuhrH9cK+Yvp5fYdVqkmilHSySlyN0HFCWo9Jg4g4jhHCIdWSOE1xrRBH2qSWRgsJnoeSgJ/FciwcHUMuQ6bJpdjiIH0XJW0e2+nwxS96RBGkqSCOG+mijRsTXBcsS+O6Jl1kGK8GM0MwjksXGtVnjwSG5IwVACjXp15LGXt2kLEHH2Pkl308tzegEgskIW3WOMwLiCKHJHRJ6z5SVcgo0EoiMjkmlIMUTWS6LVqbPFLhsnBVCdHp8NgTLv/fvzvceadLkjTuJpKy0dvgSNro3ntr08HBpIsM45UzAcF4Sfa2rWS/8rcA03cYybFRtO1QDwQTB8YJH97K5DMjDA3XGYtz2ITYBR+3PkbdK+DmBC1xgi7lQaWIjEfa1YXT3kRuXjct831KBZ/uFW0E1YTuFSV+tiXmdz+YJQyZWid48d4Ga9cq1q6NTtrzYxhzjUkZGS/pxe4wijdsRCNQUYKemCCoxlTrKeFEjEKiiUjqEcq1kanG90DbGYSwyNgKrR2wbWzfR+VsZDaD51tYQuHlHbyczXe/60wFg6Ndz1zX9DYwjNlmZgjGSzrmDqPf/QOS9Reilaa8f5DK/sNUnjjI+AEYSnIcpolWqnQlo2RkQC1WaJHF8wWR4xBqm7TQQnenQ9PyFryWDJ3zHNq6O5AZD6UsHn5E8t3vOjOa3Fx3XczVV5tdyIYx20xAMKbZ27bibH6IeOMbSNauR46NNkpcT3VD06USeB46TLGCGlaiiDq7kAWJ2AfLPEUmzlPMeHQv84iKHeTmdZJogUptSMqI9gW0tmiazukk3+ziFVy8nE2SCrIWPPwwx5Sj+Pznw5P91BjGacEEBGNmO8w0Acdl/N77iDe+AVwPHUfgOETnXYAKItJyQO8PttL3eD9PbJ1kkBLD5ChQpovn6GScibANMV+SpDbVSFKLfUhjvPIglttOLkwRlkQ4NnHSuHsIYHLSlKMwjJPFBITT3JF2mIQBaD21XhDhbH6I+n//MON3fw/nJz8mWrOealsPPNfH+O4h+p8coG/Co4JLDQ+FIsMA7QQokSPNe3R22YwXmrAqIU6hgFtqoqsdSisX0tzl4eZdhGwEg0wGHnlE8uUvHxmZKUdhGCeaCQinOf+uOyAMEFNrBVoIcNzG7ABIXr+GZOlykjAlPjCAkC5BOSaupqRhQoJFHQePGpIEEAhLI2wbZfu4voXEJhQelm/h5B2ka2O7EmlLbBviGNIUtm415SgM42QyAeE0Nl3F9MjCse0QvO96gqvfO90bmclJ0uFRhveMMf5kL6NDdaSCPXtt+sqCfTRhkVKkD5eYmIBK0oqfZIjcJsJYQL6N1u4W5s13yLVlKS3KYGesRqooAdcFz4NqFVOOwjBOIhMQTlMvWsX0fe+neuttjS8IQ5ASXa1R7pvg4M7DTDzRRy1w0DplMi2QFlyy5YQe9uEh8HGx6ER0FwkWLaXndctoL2SQpTy5Zpem7kbrPseTJKkgk2nMDDwPHn1U8g//MLN6qSlHYRgnlknQnqZesoqp1o3SFI236ugkJSoH1CcjklpAVI2phYoQ0ahMSkBCHpeQDAkxDkkmj+u5aA3KlWRasuSb7Onv7biN/QVp2uhrICXcddeRZjdg0kWGcXKYGcJp6kVbYq5Z17jN54h8ntT2GH6uQm3PMPt3jVCWzdi5ElW/hAjr5Bmkm35cNDVy2LkW/JYeSmfMo+OsNlo6PRwfZMYj0B7CkgShmE4TSQnbt0vuuMOZfljLgttuC0y6yDBOMBMQTlNybBSEaCwmCzHdElNbNlpIhGgs7abSIc0V0R0lrO6Qrp4OaMkiygWanJDDP8nS1tFGXLVwetrofN0S/GUrWHaOTevCPLn2LJZOSG0PqSwsq7HP4EiaaPNmm95eMb2YLAS8//0xN9xgZgeGcaKZgHAasrdtxX50x9EZgtao5hZULSDoH0dojchl8dps4kMD7Pn3XfT9Yjf9Kgd7PNzlLvXyIMnIfqphH5P7Y2rkUZPNqLyiPTdGsLyFJEqJgxQrI7EcSRIdKVIH3/62zcc+5qMU2HZjVgAa1xVcfXV8cp8gwzhNmYBwmpmx74CpGzylRI4OE4zViSZrWFJjS0Facxl8apDB3joV5VKnSJpalGxFT9BLwR8jE1bopxsPl2ThItpWdDH/zCJ+ewkrlwHHQbsWUgo8r7GMsH275GMf86d3JCeJ5vrrY+bP17z1rR7Ll5tUkWGcDCYgnGZedN+B6xFtuKjRDjOISLRA+ilaaaJ6QhgoFIoEjedKHKmRKkIKkEgSJBmZwXUllitxih6eL9EatLQQ8kiRusbH3Xc708EANFIyXauovd1jaOgkPkGGcRozAeE0Ym/biv+db83cd3Dd9QTveS/B8pUM7+wneXo/VsbH1yVUVzuHxrMMHA4YAQSHSCObzrJFYXw/FRImiahTQBbaWHxmiTPWttI8r4hf9HCLPn5GHOmnA8A3v2nz7W8fWUDWZgHZME4hJiCcRvy77oAknr7VNPqdy6l87kuoSo2x3aNM7B3GmgjwvQwqUVRHAvbsHKZKDkkOj73Mc+p058PG98h0Ua1Dy9KFFJd1smpNka6z2vFzFnbOw8822lwesW2b5KMf9Z+3gKzNArJhnELMPoTTiBwcnHGcdnSA1qSxIgkSSBNqtZSJ/ipRkJJGKWmYAhJBjAPYloQkRCPQtoWNhe3aODmPfMnFcxstMoU1MxjA8/caHC1NYRaQDePUYWYIpwnvm1/Hvf/fAY6mi95zLWmsqPWNMvxUP4e372Fgf5VCu09QbGPZfE10YJCAYTIM4VGjLQhxxAJ6KTBa9iljs7inxJmrm2ld3IT2M7hZ+5hU0bZtku98x6SKDONUZgLCaeDFy1RcT3XVOpLhSQafGqZ+sI+gb5wDo3myTguHn/TJdYzjWTGLqGBRo0hIoegyIT1GM8txzuwiM5ayYkMXy9e0kuks4GQs/Iz4tbODyy5LTKrIME4xJiDMcfa2rWQ//9ljylTU33MtSaRIg5i4FhBWIsbrENdjnLRKFEuUtEEn2IAFSCKkBVrYCNfD81ySAuRLPn5GopAIeWwwABgcFDOOOzr0sV9kGMZJZQLCHDa95yBqdKvXQoBlUf7sFymvXMfwwRrxwRGGnhyg9kwvw70p0YTLcxmXs85OaVtaRGZLZHkCiCgxSVjtwpEu3tIucj3NtGVc5q1qx2txIZPB82emiqCRLvrxj4/8qDXaYpq1A8M49ZiAMIfN2HMgJdGmN1L9848RnLeOpJziyQThWXitRSo93TQlNQqxhbe0yJmrHRwHul/fysRIBxMDeSKrGffCVcw/bzldnYvofl0bSjoU52XJFi00gh07JFu22GzcmEyvD2zebM+4s+i660x/ZMM4FZmAMEe92J6D6p9/jPLKdVQP16lMxIw+eYj67oMcemAXB/bFjFckVbIsGtzHkBOxbF03k3tHCQYGiahSTyOqj49glbrJdgrsnE+iLVIs6oHgu9+1+eQnG+UoXNfl3ntrrF2r2LgxwXFcQOM4ZnZgGKcqExDmqBfuOahd+juEr19PbTih3DdJWImYeHwPYd9hxgYj6oFFhCQFLE9jV4YJ+8fwwhHq+BTRjNFJZn4HxcVtLFhZom2eR4qF6wm2bZN88pNHy1FEkWbzZpu1axtNbu69t8bmzTNnDoZhnFpMQJijXrjnQLV1NPoTJIqkHqFqdaKJKmGsCGqaNNFoEnwCYtVKpC1EEmPJFAdNjMZC4BcKJL6Lk3XIFS2qNYHW8P3vH1uO4vn9DNauVaxdG53Q58AwjJfHBIQ5yN62FfdH/wE0Zgfadhi98lpspUknK0w8sZdq3xgHNz/H+HBKpVZnHJcaAh9FODpJpnkxXkcT9sgQmjIVLKRfpHTGApasXci8MwpIqxEMbr/d5u67zR4Dw3itm9Wdyj/4wQ94y1vewmWXXcZ3vvOdY87v3buX66+/nre97W38/u//PhMTE7M5nNPG89NFANU3Xs7E69ZTm4gZ3bmXYP8BBv6zj8nDIcp1CHBwyNKEpoCktU2TZItUhup4TS7+vE4yXjddG5fTccFiula2o6VFFMEvfyn51Kd808/AMOaAWQsIAwMDfOlLX+L222/ne9/7HnfeeSe7d++ePq+15o/+6I+48cYbue+++zj77LP56le/OlvDOa288A7/uLWDNNGoeoiKYqJQUR+vksZVYmWhcWgsDSs8KijHQ7kultRo28HybZwc2MUMblMG27PQutHX4JFHjt5BZMpRGMZr26wFhM2bN7Nhwwaam5vJZrNcfvnl/PCHP5w+/8QTT5DNZtm0aRMAH/rQh7juuutmazinlWTVucDRwDC26FzKQ3XKh8YY3z/GwR0jTOwtU8ZDTlYJcRCktDBGkXFyk4fw66MUF7fQvPJMorYFqIVLcZafSc+KJlxPEIYQBNDcrKbaLzce7UMfikyqyDBeo2ZtDWFwcJD29vbp446ODn71q19NHx84cIC2tjY+/vGP8+STT7J06VI++clPvqzHaG3Nzzhuby+8skG/BrzkNW7ZAj/9KbzxjRBV0VIilEJLSUs8jmxzcKMCwcJuRhYP0rq0QPp4wrLXCdpGXYqZBDEmSd0CKy/I0LphGZ0LWrHfdynF3VX8jCS3fAFLzmrCy0iamhpdzg4cODoEKWHePI/2dm92rnEOmevXONevD+bmNc5aQFBKIZ63ZVVrPeM4SRIeeeQRvv3tb7Nq1Sr+5m/+httuu43bbrvtN36MkZEKSjXemba3FxgaKr96F3AKeqlrtLdtpfkdV0IcgeNS/vTnyLsexBHKcng06GLwzh1krZDH//d+xnc8R00HKODxwxABdSpkUGhc9qUF9KIzaB6r8ux/DtB7MAUpWNjUit9v4biCKILHH5d84xvZo+Ow4dxzawwN/fYzhNP5/3GumOvXB6/da5RSHPNG+vlmLSB0dXWxffv26eOhoSE6Ojqmj9vb21m0aBGrVq0C4Morr+Smm26areHMaf5dd0A01aMgCrF+tZO9X70P56c/YVetg/G0A/HcHqq2ZHw8InJcgsijwH58mihQxyEhIk9mno9achbZtgw6VTilPAuaHCqTIX7BJUkFMgXXhX/915m3mr75zWaPgWG8ls3aGsLGjRvZsmULo6Oj1Ot17r///un1AoDzzz+f0dFRnnrqKQAeeOABVq5cOVvDmdNerExc+Zz1PPeO/8ZQbhm6PkEcpZRDQVTT2FGKRYBGEAMuKRpJbLl4GZdCcwYn4yItgWVJUqVxfRs3Y2FZTC0iw/CwKVhnGHPJrM0QOjs7ufnmm7nhhhuI45h3v/vdrF69mhtvvJGbbrqJVatW8Xd/93d84hOfoF6v09XVxec+97nZGs6cpgvFxp9Tx9E55zI5mrD/0cMM79yHVRnGtyEodRKMxVQJcakxRkyJKjEVUtqwch6qdT6FJe3k2zP4TR7dy4r0HUqxXAstLWy7sVZw9902999vCtYZxlwyqxvTrrrqKq666qoZf/e1r31t+vNzzz2Xe+65ZzaHMOfZ27aS/crfAkz3OggPjTDSMUm8bx9utYyQinIlYcJxsV2w6i4FyihyuDgE5Gk6o4Pms3oozW+ic3mJ5vlNjRlCIUvTfIlrK2SmkSr65S8lH/uYbwrWGcYcY1povsY5mx86ptdB+fyLcFREVA0ZL8cEtRCdJNQjUPU6LpoAlxhJjgQLCa6DZdvkm7NkShmk06heKi2BbWmSVCCkwLZhyxaz98Aw5qLjBoR3vOMd058/+OCDsz4Y4+VTpZajFU2B8n/9E/o7V3LgsQEOPDrGUO8gQ/vHONQb4B3YRxTVCKkQEJMSUccDIVm4rMT8FS20LG8j31FA5jJgWcTSw8m5OHmPYlOj+U2pZPYeGMZcdNyAoPXRRcIvfelLsz4Y4+Vzdu0EmC5ToUfGqR8awhrYS2WggkAQkiXFJqNqlHzJQhTtOPiU8HtKtL2+iwVre5i/tofSim78rmZENgNC4PmCXEGSyx/thLZrlzX1aI1HLZdf0BHHMIzXpOOuIbxwH4Fx6nlhVdM00SRBjAojdKWOIgU0KS44MRmhUQgUNhaCjB2SLc3HyXu4rsDOOAjHBiEQgBCgNTz6qGTbNpuLLjI1igxjrvqNF5XFC/siGiedvW0r7o/vB6aqmlo2vRe/i8nhANk3gsMwKVUqFAmA+VYTftGiWgnwSZGEOK2vo/vcDtqWNuMt6EA2ZZFW4/9aCPC8RgvMa6/NEsfgOC633hrgug5xbBreGMZcctyAMDk5yY9+9CO01pTLZe6///4Z5y+77LJZHZxxfM7mhyBJpu8uGn/7+wjPW0vp0d2ECzOUDuXoP6zI4lHwc7S/cSEFRzC0oIAulcgM9zH/8nPoWdVCYWkrVkseYcsZPZG3b5d84QseUQRKNRaRx8Yk//qvpuGNYcw1xw0I8+bN45vf/CYA3d3dfOtb35o+J4QwAeEkExMToFRjaVdr+lvOpPfRXgYfeo7+//UsE0EvMYoKLaigzuKdisliN5XhCGv/QaxCjCiPQs1Fel1IW7Jtu5x+oQd417uyU8EApGzMCI4EAdPwxjDmluMGhOcHAOPUcsz+AwT2QB9j9JH27iMIIlIgJgu4WNjUPB+nKUPWtrBVheKiZlqXtdM0rwnhe2zbbvGudx1NDV1zTUwcN2YGUmo2bUq55ZbQzAgMY476tWsI1WqVf/u3f+OZZ57B931WrFjBFVdcgeu6J2J8xktoNMFJpvcfaMuit3MlSaRxSVCEgEIh0FhYSFzHQipF6tj42QK5JgfbB+m7CKsxMziSGjpyE4HjADRmBiYYGMbcdtzbTp977jne+ta3cv/99+N5jZLG99xzD1dccQWHDh06IQM0ZrK3bSV3y5/hf6cxezuyGW3zm/87jx7IMbbjcfTuA9Top4JimDxD+ABE2XlYrkPWDsk1O7QuLtG8vJtMdwlhW0xloACNUrBqVcq999b46Ecj7r23ZoKBYcxxx50hfPnLX+bmm2/m7W9/+4y/v/vuu/nCF75g9iacaFu20Pyut0EYgNbTi8ljV17DVvsSvIF9TEzWGR6dRJKlSjOQx8clb+fxStDSU2S+r8mf2UNpQZH8ghZkzmfbdouvfOXIrK9RjmJsTLJ2bWLWCgzjNHHcGcIzzzxzTDAAeM973sO+fftmbVDGS/jmNyEMEFPpHC0EeD5jb34HcT0iKAeEYSNJlKAQ2MRoXFzwXHwHfCcFG9yih1fwkI4FQnDXXTNLWVsW0wvLhmGcHo47Q7As6yXPmX0JJ5a9bSv88z8fLVNhOwTXXU/9Xe9lRC+kMrQFeWiQZHiCfiRDeIyQQyLIUKc1n6Gtp0jXihJdzUUKS5uxcx7Ctti2TXLHHc7UIzWCwW23BSZFZBinmd94p7Jx8tjbtpL9/Gdn7DmoX3c91c//DXGoCLYO0dbtgyoxZrmkw1k60xCFR4mA0rmLOWNViXnrlrLo7AJZL0UU8o061kKwebM9PTsQQvP+98fccIOZHRjG6ea4AaG/v59bb731Rc8NDAzMyoCMmbxvfp3CRz/c6EqjNVpKcD3Cq69Fa6iUNU9uGWTPfTvw6Cee6nNQpos6GomF8+TY/9/efYdZVZ2LH//uesqcc6bP0KuAAjYMilixYAEs0aigklx9NPnFds29JthibozRGBN/P8xN1LSrEXPxaqIhNwE0xjpAxgpoUEA6TO9z2i7r98eZOTBDUcphmOH9PE8e3Ods9l6LnWfeWWvt9b60l4QJ2j5BLYWmdHAcCAZB07osJiuVWUwWQhx+9hgQrr766t1+N2vWrAPeGNGVWbksEww6RgZoGunTzyR+x524E0/Cc6F6s0PNijUEqSeFjSKNBhhEKSCOIkbBiEJCw/IpKrUwwoFMUQPTBNve7WKyEOLws8eAcPPNN+/0WTqdlj0IB0n3WgeaaWaDAWSSzjkpD6/Vw8XrSGSXBAKARh46CcCKRYiEdSxbywQDpah8z6bi3SCbt+gdowNZTBbicLfHgJBOp7n33ns555xzOPfccwG45ZZbKCoq4v7778c0c1pw7bDXvdYBt9+eDQYATlqxZW2CzR+1YtNEHu00oJHGwMOlAYsoeRRHfIYODxCM2uD7PP1CPnO+G8P3MwMF0wTXVei6LCYLcTjb49zA3LlzaWtrY8KECdnPvv/979Pc3Mxjjz2W88Yd7rrXOqClJfud68Jnq102L9+M2bKVCD4KE4Moiv4EKaSYGAOPKmDoWJvyMcUYIZvKlXnM+W4M183sSHZdmDnT4c4707z0UlwWk4U4jO0xILz22mv85Cc/obi4OPtZeXk5Dz/8MK+88krOG3e4617rYEfpNKTaPVQySZo4aXwyj9NHoQhgoII2Zp5NJKhjhzKjuYqldpfyl7qeSV99221S9UyIw90eA4JlWQSDwZ0+j0Qiso6QY4Gnf4u9+K8A2X0H6trZuC6kUtBQr1j3STsblm5Fow6bKtpw2UIIFws3nE8w4DM4P06NPpifvzicyo8iFJYZUv5SCLFLe1wE0HWdtrY2IpFIl8/b2tpwXZlayJXs20WdC8qaRmLWtXDcydStbSWVgg1rXVKbNmPXrCdCgmaG4mFjEaOAECPGG1iBUtqj/bjt/52K4+hYT8KVVzro+vYMpvn5Pd1bIcShYo8jhOnTp3PPPfcQj8ezn8Xjce655x6phZBD3d8uwjBIXjGTdBrSSR/luHgpB5VO4qQUHjoeCvDRMDA0m2BegHBQ4+/rxpBK63i+htNR2My2wTAUti1vFAkhttvjCOGrX/0q9913H6eccgqjRo3C933Wrl3LjBkzuOmmmw5WGw873d8uin/jZpLHn0TjJp+Nq1pwUopUUytebRLbaSFFggSFxAmh4TBqlEmsxKImUcLClSOzU0SGkVkvuOIKR6qdCSF28rlTRvfffz/f+MY3+Oijj9B1nWOOOYaysrKD1b7D0o5vFyldR+Xn43mguQ4xK0VSt8nzfFLDbLyjC1i3YgA6cWLYWIQYc8YA8k8YzYo3B+P5Op0bzmbOdLIBQDKYCiG622NA2Lp1KwMGDGDgwIEMHDiwy3dvvPEGp59+ek4bdzgKPP1bgs88BWxfTHYmn4qTcHlnUR0bPtqGBhSEfZreX0vLilUYbKGd0cQJEwXqNyTJH9tGQaANv2N0ICkphBCfZ49rCDtOC91yyy1dvpNaCAde98VkNI3EzKtJTziRqk9bqPtoI04iTXNNKzXb2gk1bMQjDRTQRoggLiEUdgkU9Quytq6k48qZnQwrVuw+e60QQuxxhNBZRhFg06ZNu/1OHBjdy2JiGKSumImT9GhrdkinFH7KwW1L4bo+npvCQeFhYGKgYxJAkRcOYYUs6uolAAghvrgvnP66eypsSY19YJmVy3Yqi9n60E9IHjuRDe81suK9VtZ/uBm9pY5AUMOwPJo31JMmRR0x6tHoTxNl6AwYW84f3xrE35Z2vlOqMM3MgrIQQuzOFx4hiNzKjA6c7OggPfUCktf+C03b0qz4oJXarUmsms3YKk2eF8c0NbRYPskWD4s8+uFRToDBRwb5tHEA3/2//bM7kjVNcfXVjrxRJITYoz0GBN/3aW5uRimF53nZ/wbwPFmgPJC6p6nwyspQCpIJhRN38BvbSDZ52MWg0mk03UUlErhYKHR0bEK0Ywd0FlT0w3U1dsxgKqMDIcTn2WNA+PTTT5k0aVI2CJx00vZMmzJldOCYlcuw/7YY2P5mUeqKmbiOYsv6NHUbU6z/pJ0yaqA+iUInGC1DmQ6uAxoeJu0UkKAmfzIvvtGv48pSDlMI8cXtMSCsWrXqYLXjsGZVvLVTeUx34kk4rT7K9ykv12jr7+O3lBIr1QgYHv3HFZBsiVL2eitWsolQgcFRkwfypjYB18uMDqQcphBib+S0NNaCBQu48MILmTp1KvPmzdvtea+99hpnnXVWLptySNM6alhmtgwo3PHH4rs+ba2K9yta+McfVrFqyce4DWtp+GQN6Y83EP9wM5vf3kZrso4WHFJNcbwt24gF26UcphBin+Sswk11dTWPPvoof/jDH7Btm6uuuoqTTjqJI444ost5dXV1/OhHP8pVMw55ZuUywo//DCA7QtBra0g0p6nZ4FD/WTXBplrCJNBwcAnSQgijRZEO2OS1OURoJ9+GgtEDWNs8pOPKmfWDzN4DGSEIIT5fzkYIFRUVTJo0iYKCAsLhMOeddx4LFy7c6bx77rlnl6U6Dxe73Htw4sn4moHn+niJJH5zAp8kPhoKUGhoukYoaqHhY2AQjEE4EqC+fed05UII8UXkbIRQU1NDaWlp9risrIzly5d3Oefpp59m7NixHHvssft0j+Lirmm5S0uj+3SdHrNkCTzbde9B6sGfYk+cjNEeQEs5pD/ehGIbARqpw8VAEUTDsUdgaa0E0BmmNVA+qJj3uJDXlsbo3JlsWRpf/7pNaWnvql3R657jPujrfezr/YO+2cecBQTf97u8iaSU6nL86aefsnjxYv7rv/6LqqqqfbpHfX0bfiZZD6WlUWprW/ev0QdZ5IlfEXScbInMlsnnsOW0y2hf1UJ1rcHq1z7DrPkMF40wJq2ECBDExsbKS1JQahLVobR0EFuKj+M7T4+j821gTVPMmuVwxBEpamt7rIt7rTc+x73V1/vY1/sHvbePuq7t9It0l+9zdeN+/fpRu8NPotra2i5ZUhcuXEhtbS2XXXYZN954IzU1NcyaNStXzTkkdd/2ly4oRilFKqloa0qTSCXILDU7KHx0PDQ0FD6GbRMIRQmGNew8k+WNo/A82XsghNh3OQsIkydPZsmSJTQ0NJBIJFi8eHGX7Ki33norixYt4qWXXuLJJ5+krKyMZ5998eHy0wAAIABJREFUNlfNOSSpaCzzZ8dx05DxbPwkQc3GdmqXrmb925tpIkmcFpoAjxQeHiYeRcVB7JCiMASFQwvwQ7GO60hpTCHEvslZQCgvL+f2229n9uzZXHLJJUyfPp1jjjmGG264gRUrVuTqtr1G97eL0DTc2iaaNtTS+GkVyY2bMJKteAQwCKPIoxCDIgKUDTYI5SuGjQpTPH44m8tO4zevjeu4soamIaUxhRB7LWdrCAAzZsxgxowZXT775S9/udN5gwYN4tVXX81lUw453d8uUrpB1ZAJ+IkUbQmTeLuHlUqiMPBQmTeLAA1FMGpg2gZYIQLRNO+sKd5pukhKYwoh9lZOA4LYNbNyGcHfPwNsf7to083fY0VzOX5dLTXb4qj2NImEhoMCdDwMfCCAw4AhRZSPDhEssrH6lVIaKpDpIiHEfpOA0AO6p6povXI2LRfPIrpoG07+AELuVsL9o1hlcfw3WtCxcTAZRJTBZw+i7NQxDDhhAEZZKbqbpmljAZoGSmXSVch0kRBiX0hA6AF+YVGXVBXtQ4+ianUzyxd+SqI+SbK6FdPS8eNxmqnCAHQMmkhQvq4a7ag8DFWGaSiwAhSWaqgdSmUWFsroQAix9yQg9ABrxYcA2fUDe+UHfFj/JVrrkljxOLqTwHMUIVrwieMTJUA7GgGKR0QxygtQxWVQVASWxYqVVseVJV2FEGLf5TS5ndi17vsPPA8SzSlIuaRbfDySmDj4HcvIGhDAw8YkaBnouoYWDkEgAIZBTY2kIhdC7D8ZIRxkZuWybDGczsBQWziS1RX1tNVtI0gDOhY+4OPi4GMSpx2f0XaCwmHF5E+agFFcAJrG00+bLF7c+RilVKYQYt9JQDiIzMplFFx2EaSSQGd2U53mz2qgtgwLnzx00rgUYmCSYivFFNOKTz5HTTKJTTya4LD+aKZBZaXOnDlBKZUphDggZMroILIq3oJ0Cq2jAp3SNHzLZk3JsThtLi5JXAwsFB4pdAJ4gIWODeiWiRkJYAQzawYVFWY2GEi6CiHE/pIRwkHU5e0ioP30c/lowuVUvjMUhzfR8fBwUTgU4eDhECCIi0GUNINGDSE2fii6ZQCZt4k63y4C2X8ghNg/EhAOIr2xATQNTSmUppEaP4H6USdRtLWRViOM8jwMFP0G6/QrAwojBD/YQv+jBhLML6JgynEQ3J7KOvM2EYCGrsv+AyHE/pGAcBCYlcuwKt7KlMpUKrv/oNW1ePOP63nv5SqirMEnRIh2qjcZxDalSRfEaGtKE1lXg1nejmkdh2dlHtnTT5s880zn66aZxWRJVyGE2B8SEHIsu5DspEHTdhgh6KQ21LDxoyYKqEJDx8UjhI/CJmVbRIvCBNvjFI8fQVGJTeiIwbRb5i4Xk2fOlMVkIcT+kYCQY8Hnfg+pZCYI6DoYBkoplGmxpf8xtCxMYpFCJ7M0rKPIwyVgauihEEa+jWnqBGIh3lmdz+L/DrB5i47vgywmCyEOJAkIOWRWLiM473fbp4lMi9YHHkLbuo2qvMH8o7KYFB8QoAGXFBYtpDEpIIVu52OrFFpRmNIRBWzNO57rZhWTToNhgGmC6yp0HR56KCmjAyHEfpOAkEOZFNdONkVF+uxzSM24lFRtC8v/dwv//GcdxSTwMdEwiAERLIJYFAztR+HYfmimQdGJY/nV/4whlQLQcF3FBRe4TJjgM3myK8FACHFASEDIoc4dyZ28kjJQCtfXaW6Ik2yMo/A6Kh64WGgoQMMnGNPRQyECYQ0jHEDpXdNTlJUpbrstfRB7I4To6yQg5Ejg6d9iL/4rQHa6KPWVK1GNjbSta+LtRa1saq2njAQmbVj4pEigo1NqG5T2C2GUhSgcECE8uIzCss7XTTN7Do4+2uuRfgkh+i4JCDlgVi4jOuffwPOyNQ8SV1+Lc9wJJFZ8ysoP22jfUEsZaQxMbMLYJAhRSDkJ8o8qp+ToQcTGDkJF8lm5sYAnfh2mo9gmmqZobJRN5kKIA0t+quSAVfHW9mAAYBikrpiJ8nziLR7NdUmSOLgoFD7gYXe8Z6STRyBfJxAJYIQD6AGLJf8IdLxiClIiUwiRKxIQcsAvLNr+ZhEQ/8bNuF86kVRrmhXvtVL7aRNhqrFowsMggYGLIo8mBsSSDJk4CmvoQDwrDLF8ikroSFGRISkqhBC5IFNGOdClAI6uo/LzwfdJxT1SVhQG5WMti5KHTgM1lBAijMbRR5cz8uQyCk8/AXtQOVbIxA+GdkhRAbqOpKgQQuSEBIQDLPD0bwk+8xSwfTHZOfkUSCZJ1MdZ+t8f0/rZWnxqqMVGoVNPAosE7tYULTUxCjzQ2lMY4QDvfWDx7LNW9vqSokIIkSsyZXQAdV9MRtNIXDULd/wx+Mk0q99txN9QQ6DjVVMIYOKSh02EOPbQAYRGjcCI5BEoimBFAvzxRQvXhc4F5bPPln0HQojckIBwAO1yMfkrV4Hn4aU9Eq0JEp6DRwrw8dEw0PCBECZ2xMTIz8PXTdB1DEtH61Yds6ysewFOIYQ4MGTK6ADqks0UiH/9Jtwjx+Jt2kJjnc/WNWnSVOMTx8WlFRcThygthGknNPoIwkeNQg/aBGIBNF0jGu28WmbLmuw/EELkigSEA8SsXEb48Z8BZPceqFAYf1sVDau2sr42yLY1tQSANBHSeJTgYxCgDEXpyAIKxw8lMrCIYMxG0zUqK3Uef7xzQ5qGpiH7D4QQOSM/XQ6QXU0XORNOINnq0NoGSUfHVwofhYHCwUfrmDSy0LELgwRK8zuSV4CmazuVyJQFZSFELskI4QDZae/BtV8jXjKANRVVbNmo2NDoUNsexcQniQvYKFKE8Cgt1smbfAqhfoXYliKYn5ku6l4i8/bbkQVlIUTOSEA4QLqXx/SDebS5QZr1GN6QIuw8g5FqI1vWhSjFopFtlDKAgpjDkJnHM/j84wmVxwiGFbqRmS76858tNA2UypTILCjo6V4KIfoyCQgHgFm5DPO9d7uUx2xq8Vn9yqdUVjhs89K4roH//qeYNFKHhkaQVqoJtERoao9SkDBYV+nz3sd5FJUZ3HNvkHQ6s0NZ1xW2DWee2cMdFUL0aRIQ9lO2RGYqCXQuKOs0r29kZbNJbb1LKNJOKu5QwFoa0Ahh0oZGhCCxiEGkPMqaLVFu/nYJjpOptOn74PuZkcHpp3vccUeKk0/Oo7a2Z/srhOi7ZFF5P1kVb0E6hdaRbEhpGr5lsSU2AjeRAF0n4ZiYbhKFA1iAl1lIxiRYEsKIRXl/uUU6DZ6n4XmZFBWGkRkZ3HFHStYOhBA5JyOE/aQ1N4PvZxeTq48+nX/0O5uKjwza/WraakIYRQ60NFFMPXEsbExa8BiDouCIIZQc1Y/kKqujTrJCqUwCu/x8pCKaEOKgyekIYcGCBVx44YVMnTqVefPm7fT9K6+8wsUXX8xFF13EN7/5TZqbm3PZnAOu+94DNI21DGFFSzGp1hThgCIUTJOntVESTROnkHxsivDoj6J40iCKz5nI+vpCnno23HHVzH6D/Hy47TbJaiqEOHhyFhCqq6t59NFHefbZZ3nxxReZP38+a9asyX7f1tbG9773PZ588kn+9Kc/MWbMGB577LFcNScnuu89ULrB+tgReCmFkwDfSeG1p1DKJdVmYKGw8DEJYmEQilpY+VHeXR7sst9A6h0IIXpCzgJCRUUFkyZNoqCggHA4zHnnncfChQuz3zuOw3333Ud5eTkAY8aMYdu2bblqTk50T1Xx0fGX8d5nNtXvb6Exqaja2ELcS9Jc1UpzeyutuOi0kqKVQhopPWkSkSMG0tqud9lvIPUOhBA9IWcBoaamhtLS0uxxWVkZ1dXV2ePCwkLOPfdcAJLJJE8++STnnHNOrppzwO00XYRGU007nuOgMIhjE0YDdAaabcQ0lyLyKMSmnALKi0qI9M9n0zabp57qmp5C6h0IIXpCzhaVfd9H2yFVp1Kqy3Gn1tZWbrrpJo488kguvfTSvbpHcXGky3FpaXTfGrsvPqyks65lZ6qKzZFRGO2g4QABfDLvFFmWhZVvEUso7HgUK98g1j9A/6ElPDU/ryO9dYZpakybFqC0NLDL2x7UPvYQ6WPv19f7B32zjzkLCP369eOdd97JHtfW1lJWVtblnJqaGq6//nomTZrEXXfdtdf3qK9vw/cz0yylpVFqa1v3r9F7IWDnEd2hruU/jzyX97eFiTc0spECArikacemGUNFCRWGiYZMEvE4njKwysp4Z0Me8+dvv4ZhwIMPJjniCHeX+w0Odh97gvSx9+vr/YPe20dd13b6RXpHOQsIkydP5rHHHqOhoYFQKMTixYu5//77s997nsc3vvENLrjgAr75zW/mqhk5YVYuI/Dnl3ZIVaGTDuaTd+QI2pbXUeIpjEQ7/YeWoLU7DDurP+X5IdLxCFpzE0ZhjJKxQ/jT6mC2+I2mKa65xmH2bFlMFkL0jJwFhPLycm6//XZmz56N4zhcfvnlHHPMMdxwww3ceuutVFVV8fHHH+N5HosWLQJg/PjxPPDAA7lq0gERePq32apoKIXSdXzDYtnGEKtr/0kSaCBIf2oxNrSSJoqz0qAmz8BNp1BOGrNflIJRira40WXvgdQ6EEL0pJxuTJsxYwYzZszo8tkvf/lLAI4++mhWrVqVy9sfcNkSma6brXnQPmESr0fP4IPX81DUA9CPenSSaASJammK7DpajRKs0mK0aAGBoeV8oo/lyd/suPdASa0DIUSPkp3KX4BZuQyr4i2MzZvA97vUPNhy2f9h/RtJNH8DPkmSRNBRRFGEcLGjNmkjimmaWGGNVDgMsRiLKop2qJUsew+EED1PAsIemJXLCDz3e0K/nweem1n1NS2U64CuU3v3Q6xsLmZ15YfUsAUfhU0zGqDh46PjuSa2qeEXxrDHjiAdG8a7Gwfy0qJYx10yweChh5Ky90AI0aMkIOxGlyymSmVHBYlrvoo/aDCpE09hkxrE1l8sIVWXwCYEJIEwYRIU4VMyPEqLyiM88SjU4BGUHTeM2jX5zP1VJLszWRaThRCHCgkIu2BWLiP84wezWUwVmfUCLJvkV2aSPn4iXtojuWQbbe0+Pgl8khgofHwMTIKk0HUTzDDBoEEyGkKP5PHKq4GdpoquuMLp2Q4LIQQSEHaSHRmkU5ksproOpkly5jUkvjKTtlHHkqpP4yZdWmrbaahqxaUeDYMkKcAiQgthAB0iUZvU4CMpOnIY6zfbLFhgddxJpoqEEIcWCQjdZOsb+D5K00iffibt/34nzoSJKF/htTh4SsdJ+fhWiPxR/WlfU0OSJgoI4aMzYkAZeYNLKT9xOHmlUcLjyigYbPO/L5uy70AIcciSgNBNl/oGShE/bwZNQ8ajauKgaTQ1Q92GVlKtKTa8von1f/0Yi9VoBGkHHKJUba0l5ISIDW3DLygkGg5gWVBU5Mu+AyHEIUsCwg66J6xTmkZ6Yw2bPqglUhbGTXgkVAC3oZb2Rpdta7YRopU0AXQUBgY+CkUpwUHl5B89ktjRwykutwiGYOVKo+NOmfWDFSsMQEYIQohDgwSEDtmF5B3qG2AYtE84mXTCxU+5uCkPTxn4jo/rA6kkChcX0NFJoxHAQOFiRsIECiPY+WGsgM4u8voJIcQhRQICHQvJl07PLCSTGRko3WDllfewxRhD7Yef0VIdxAyYBMJBEv/8jFRjkviGBppxiGPh4JAiRhmNGEXDGDOpP4UTjiCvPIJhZXYgR6Odiewyf8qUkRDiUCIBAQg+9/vMQjKZH9VtRxzNe9PvYLU2DH19G15DOwX9A6TrmnCUi9PShlOfwGlvJUUJ0IyJh47FYAoonjiCQaeNpnhwGN3OjA4qK3Uef3zHugeSqkIIcWiRn0iA6nbcMnAM24qOwk/6aKk0XjJFuj2FSqVw2jxSro6hGei+hqKdTFw1yMNARyMYCxAqDmNY26eKKipMKZMphDikyQgBcI8+FtgeGNYli1nz6iridU2Uj+2P/+kG6lKNBLDwC2LUrouTbEvTDDhouMQxsAhTR3DoQMZOG094aCno2+NtR7VNpEymEOJQJSMEwFz+AdBZChPsdauw6qsx6jZjbVtDS2uShBNjS51Lot0jNigfVRCgKBChHxplBAiiMXzwQIq/NI7SUUUYtpG9/s7TRVImUwhx6JERAqDXVHc59tMOXiJBolWnJaHhx1MEaSaRVniYeBqYho5p++R5BklXI4RGMOQTiAWwokEqK3UqKkwmT3apqDA79h/IdJEQ4tB12AcEs3IZgVdfATKTOb5msMwZQ/PqbfikadocprUxj+jGZiK2RWGJRUtdG3khm+SgclTtRqw6j8F2G/bAkYw4eywfr4twxZVhHAcsy+YHP0hi25BOK3Rd0lUIIQ5Nh31ACD73e3AcNMBH45Php9IeHI2KK9w2KBxXRFFdO4NP7k80WU/wqCKG6AU0xQ1KNtQRSFnUbvEZOUYjMGQI9aqcnz0aIJ0G38+MCBobdV54IZ4dMUgwEEIcig7rgGBWLiM473eZUpiAj85rn0Vp5y1cLBKUsOr1GP2iilC0loaAzrDRoJrq8VpcEluTJJIJ6qsV9e4QKt8bx8vvZgrf+D7ousKyyAaBiRPTPd1lIYTYrcM6IATn/x5cJ7v/4FOGUIdLCkWKGGl08kocVDSf6BEx2h2TcGEQM1JEoL4VI96MbgXZ4pfx8BvTSHt6x5tEGrquOP10jzvuSMmIQAjRKxzWbxmpbjsQWgiQwidT78zAxSRiOwQiBk4oih4KYJpgWRArsjCDJp5h81nLgI5gkHlPSdMUto0EAyFEr3JYjxDcUBTYvv8gMzpoZj2DcMknTSHHlueR309RGmwhXGxjR23idWlwFSUji1C6yZBGC/XO9iudf77LzTfLPgMhRO9y2AYEY9lSIr/6OUDHgjIE0PGI0cowSggTJEC0v8awkRb9xpegO2mMfJtwQTGJljShQDGx8jw2VPTruGpm8qmsTEkwEEL0OodtQAj897NonptdP1BorKUYBx0TG9CxDYNQXhAzYBAIAIaBr+v4jkcoZmNaBr5poxmSylQI0fsdlgHBWLaU0HPPAh17D9B4njG8TxmNhLFpI0oehYPyKDmiH/llGpg2esQmNrKUVJuDbmhoho4dsTluos68/+m8mmQxFUL0TodlQLDefhPNdbNTRcs4ij8znnzKKcKh35EljJnYn0B5GYNPHEKs2MaKBPAdDyNgkpcXQPkKTdfQdI2VH3XWSc68XSRZTIUQvdFh95NLe7sCtbQSVKZMpgasR6GjswWHWmxUnUdjnUbCt9ADmf8pBZppoBk6mq6hm5k/n37a5JlnOgOCwjQlLYUQonc6rEYI2tsVFF15SbYQjgZ4QAATH58ICh2N8LEDKZs4koIRJYRKIgRKgmgamWBgbI+hlZU6c+YEs2mtNU0xc6YjC8pCiF7psBohBP77WbR0Cr1jrt8HPHTWUIxBAB2DCD52OEyoJIYdttF0jXc/tHnsF2Heea9r/NxV0rorrnAOdreEEOKAOGxGCMavfk3e/8yjczuah8brDGcJQ/iECBotRNEpKR9E2Ylj2Jrqx8d/jzJwncm994WyiepeeCGeHQFMnuxi27YkrRNC9AmHRUBQby0l/5470Hwvu5D8FiOYxyl4+BRgYONx5MnDKRo3jLZof26dMwDHAU3L5CXqTFRXUWFmcxJNnOhL0johRJ/R5wOC++rbBB9+EM3fvufAR+dNhmHgY+ATwCCqhckrimDHwvzl1RipFCiVeWtI1zPpKDoT1e1IktYJIfqKPh0Q3FffpnTWRei+k11AVmg8zfGspZAIPkF8ii2LAeMLaR5wPAs+PY4/LIxmy12aJvzwh0kaG3UZBQgh+rQ+HRC8//wVhr89m+laynmKSTQRo4AEw44dju34DD1jKFv0sXzrV6eQSmvZjKWdbw3Nni2vkQoh+r6cvmW0YMECLrzwQqZOncq8efN2+v6f//wnX/7ylznvvPO4++67cd39/8Hb8MxC6q75Ng3PLCTdmujy3SYKeZdT+AcX0cxxlJba5PeLUdQvxodVg0k7WpeMpYGAvDUkhDh85CwgVFdX8+ijj/Lss8/y4osvMn/+fNasWdPlnDvuuIPvfve7LFq0CKUUzz333H7ds+GZhQz/1tWMXvwEw791NbWRYbjo+ICDzjNcxN+4jXe5lGf1b5M45nyGXT6Z/NO+xKkX5GHbYBiZ1NWzZztd3igSQoi+LmdTRhUVFUyaNImCggIAzjvvPBYuXMjNN98MwJYtW0gmkxx33HEAfPnLX2bu3LnMmjVrn+/pL3wVAxej48VSL1HHc/2uZ6C+ln9uLaL5yOn4q0wUOh46a7WxTJvugK5z1lEaL/SXN4aEEIevnAWEmpoaSktLs8dlZWUsX758t9+XlpZSXV29X/fUzz8Lb/GvAQ8Pg8BZ59H8jzYanbG40TZOO6qJv6zxcX2wbDh1ioYZ3P5PIG8MCSEOZzkLCL7vo2nb00Irpbocf973X0RxcaTL8Zjbv0JNYZjkn14meNG5nPC1aYxc10jTukYi/aPkFQWYcF2ct94NMeUck8mn5O1j73pOaWm0p5uQc9LH3q+v9w/6Zh9zFhD69evHO++8kz2ura2lrKysy/e1tbXZ47q6ui7ffxGNje34fiYNRXFxhPr6NowZZ5A34wwA6uvbIGYRPTZz3Tgw/mQYf3IKSFFfv4+d6yGdfezLpI+9X1/vH/TePuq6RmHh7n8RzllAmDx5Mo899hgNDQ2EQiEWL17M/fffn/1+4MCBBAIB3n33XU444QReeuklTj/99L26R/eOdR8x9EXSx76hr/exr/cP+mYfNaWU+vzT9s2CBQt44okncByHyy+/nBtuuIEbbriBW2+9laOPPppVq1Zxzz330NbWxrhx43jwwQexbTtXzRFCCLEHOQ0IQggheo/DKv21EEKI3ZOAIIQQApCAIIQQooMEBCGEEIAEBCGEEB0kIAghhAAkIAghhOjQJwLC59Vd6I2uvfZapk2bxsUXX8zFF1/Mhx9+2Gf62dbWxvTp09m8eTOQyYw7Y8YMpk6dyqOPPpo9Lxf1Mg6W7n288847mTp1avZ5vvzyy8Du+34o+9nPfsa0adOYNm0aDz/8MND3nuGu+tiXnuFuqV6uqqpKTZkyRTU2Nqr29nY1Y8YMtXr16p5u1n7xfV+deuqpynGc7Gd9pZ8ffPCBmj59uho3bpzatGmTSiQS6owzzlAbN25UjuOo6667Tr322mtKKaWmTZum3n//faWUUnfeeaeaN29eTzb9C+veR6WUmj59uqquru5y3p76fqh6++231ZVXXqlSqZRKp9Nq9uzZasGCBX3qGe6qj4sXL+4zz3BPev0IYce6C+FwOFt3oTf77LPPALjuuuu46KKLeOaZZ/pMP5977jnuu+++bCLD5cuXM3ToUAYPHoxpmsyYMYOFCxfusl5Gb+lv9z4mEgm2bt3KXXfdxYwZM5g7dy6+7++274ey0tJS5syZg23bWJbFyJEjWb9+fZ96hrvq49atW/vMM9yTXl9T+fPqLvRGLS0tnHzyydx77704jsPs2bO54IIL+kQ/H3jggS7Hu3p+1dXVOamXcbB072NdXR2TJk3ivvvuIxqN8vWvf53nn3+ecDi8y74fykaNGpX97/Xr1/PXv/6Va665pk89w131cd68efzjH//oE89wT3r9COFA1FU41Bx//PE8/PDDRKNRioqKuPzyy5k7d26f6yfs/vn1pec6ePBg/vM//5OysjJCoRDXXnstr7/+eq/u4+rVq7nuuuv49re/zeDBg/vkM9yxjyNGjOhzz3BXen1A6F5XoXvdhd7onXfeYcmSJdljpRQDBw7sc/2E3T+/A1Ev41DxySefsGjRouyxUgrTNHvt/3ffffddvva1r/Fv//ZvXHrppX3yGXbvY197hrvT6wPC5MmTWbJkCQ0NDSQSCRYvXrzXdRUONa2trTz88MOkUina2tr44x//yI9//OM+10+AY489lnXr1rFhwwY8z+PPf/4zp59+epd6GcA+1cs4VCil+OEPf0hzczOO4zB//nzOPffc3fb9ULZt2zZuuukmHnnkEaZNmwb0vWe4qz72pWe4J71+DaG8vJzbb7+d2bNnZ+suHHPMMT3drP0yZcoUPvzwQy655BJ832fWrFmccMIJfa6fAIFAgIceeohbbrmFVCrFGWecwfnnnw/AI4880qVexuzZs3u4tfvmyCOP5MYbb2TmzJm4rsvUqVOZPn06wG77fqj69a9/TSqV4qGHHsp+dtVVV/WpZ7i7PvaVZ7gnUg9BCCEE0AemjIQQQhwYEhCEEEIAEhCEEEJ0kIAghBACkIAghBCigwQE0Sf97Gc/45VXXgFgzpw5/PrXv97leWPGjKGhoeGA3/9vf/sbP/jBDz73vF/84heceeaZ3Hnnnbs9Z/PmzRx//PEAPPbYY3z/+9/f5Xme5/H1r3+durq6fWv0HqxcuZJ77733gF9XHFokIIg+admyZT2aavnss8/mnnvu+dzznn/+eR555BEefPDB/b7nb37zG0488URKSkr2+1rdjR8/Htd1+fvf/37Ary0OHb1+Y5ronZYtW8ZPf/pT+vfvz7p16wiFQtx444387ne/Y926dUydOpW77roLgPnz5/O73/0OXdcpKSnh3nvvZfjw4cyZM4dIJMInn3xCVVUVY8aM4Uc/+hEvvvgiK1eu5OGHH8YwDADef/99rrrqKurq6hg1ahQ/+clPCIfD2fb8y7/8CxdccAFXXHEFAD//+c9pamrKtqHT+PHjOfvss1m1ahWPPPIIn3zyCfPnz8dxHJqbm7nhhhuYNWsWf/jDH1i0aBFPPPEE1157Lccddxzvvfce27Zt4+STT+b+++/nW9/6FtXV1dx9993cdtttDBgwgB//+Mek02lqa2uZPHkyP/zhD7/Qv2eqosoQAAAEkElEQVQikeCpp55iwYIFQGYksXHjRqqrq6mtrWXcuHGcdNJJvPjii2zevJk77riD6dOnf+HzAK688kq+973vMWXKlP17+OLQ1QMpt4VQS5cuVUcddZT66KOPlFJKXX/99dkc9PX19WrcuHGqqqpKVVRUqHPOOUfV19crpZR64YUX1AUXXKB831ff+c53uuStv+SSS9Tzzz+vlFLqmmuuUX/961+VUkp95zvfUZdffrmKx+PKdV116aWXqj/+8Y9KKaVGjx6t6uvr1csvv6wuu+wypZRSnuepKVOmqLVr1+7U7tGjR2f/bltbm7riiitUQ0ODUkqp999/Xx133HHZdt54443Zttx6663K8zzV2tqqTj31VLVkyRKllFJTpkxRy5cvV0opdfvtt6ulS5dmr33SSSepFStWqE2bNmWvO3fuXPUf//EfO7Xr1VdfVddcc032eO7cuWrKlCmqpaVFJRIJNXHiRPXggw8qpZR6+eWX1dSpU/fqvE7HH3+82rhx4+c8XdFbyQhB9JhBgwYxduxYAIYMGUI0GsW2bYqKisjLy6O5uZk333yTCy+8kKKiIiCTU/+BBx7IViI77bTTsG0bgNGjR9Pc3LzLe51zzjmEQiEgk964+7rBlClTeOCBB1i1ahXV1dUMGjSIESNG7PJaX/rSlwDIy8vj8ccf5/XXX2f9+vWsWrWKeDy+y78zZcoUdF0nEokwdOjQXbbzoYce4o033uDxxx/ns88+I5VKEY/HKSgo2OO/I2RqaAwZMqTLZ5MnTyYajQKZtMynnXYakPm3bmpq2uvzIPPM1q1bx+DBgz+3TaL3kTUE0WM6f5B3Ms2dfz/xfX+nz5RS2fWBYDCY/VzTNNRuMrHseO1dnWcYBldeeSXPP/88L7zwAlddddVu29051VRVVcUll1zCli1bOOGEE/jXf/3X3f6dL9LOa665htdff50RI0Zw0003UVZWttv+dNeZbnpHX+Tfd2/O6/yucxpO9D0SEMQh7bTTTuMvf/lL9jf6F154gYKCAoYOHbrHv2cYxl4vKn/lK1/hlVde4aOPPuLcc8/93PNXrlxJUVER3/zmNzn11FOzC66e5+3VfSFTFGnFihX8+7//O1OnTqWqqoqNGzfuMiDuyvDhw9m0adNe33dvKKXYunUrw4cPz+l9RM+RKSNxSDvllFP42te+xle/+lV836eoqIgnnngCXd/z7zJnnXUWP/3pT3Ec5wvfq7i4mPHjxzNy5Egsy/pCbXv++ec5//zz0TSNE088kaKiIjZs2PCF79kpFotx4403cumllxIOhykvL2fChAls2LDhC03PTJ48mbvvvpuWlhZisdhe3/+LWLFiBUOGDGHAgAE5ub7oeZLtVIgODQ0NXH755cybN4/+/fv3dHP22uOPP45hGNxwww05uf6cOXM4//zzOfPMM3NyfdHzZMpICOC5557jwgsv5Prrr++VwQDguuuuY+nSpV0qeB0oK1euRNM0CQZ9nIwQhBBCADJCEEII0UECghBCCEACghBCiA4SEIQQQgASEIQQQnSQgCCEEAKA/w/a1aeIimsCSQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for _ in range(50):\n", " # Generate permutation samples\n", " perm_sample_1, perm_sample_2 = permutation_sample(rain_june, rain_november)\n", " \n", " # Compute ECDFs\n", " x_1, y_1 = ecdf(perm_sample_1)\n", " x_2, y_2 = ecdf(perm_sample_2)\n", " \n", " # Plot ECDFs of permutation sample\n", " _ = plt.plot(x_1, y_1, marker='.', linestyle='none', color='red', alpha=0.02)\n", " _ = plt.plot(x_2, y_2, marker='.', linestyle='none', color='blue', alpha=0.02)\n", " \n", "# Create and plot ECDFs from original data\n", "x_1, y_1 = ecdf(rain_june)\n", "x_2, y_2 = ecdf(rain_november)\n", "_ = plt.plot(x_1, y_1, marker='.', linestyle='none', color='red')\n", "_ = plt.plot(x_2, y_2, marker='.', linestyle='none', color='blue')\n", "\n", "# Label axes, set margin, and show plot\n", "plt.margins(0.02)\n", "_ = plt.xlabel('monthly rainfall (mm)')\n", "_ = plt.ylabel('ECDF')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test statistics and p-values\n", "- Test Statistics\n", " - A single number that can be computed from observed data and from data you simulate under the null hypothesis\n", " - It serves as a basis of comparison between the two\n", "- p-value\n", " - The probability of obtaining a value of your test statistics that is at least as extreme as what was observed, under the assumption the null hypothesis is true\n", " - **NOT** the probability that the null hypothesis is true\n", "- Statistical significance\n", " - Determined by the smallness of a p-value\n", "- Null Hypothesis Significance Testing (NHST)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generating permutation replicates\n", "A permutation replicate is a single value of a statistic computed from a permutation sample." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def draw_perm_reps(data_1, data_2, func, size=1):\n", " \"\"\"Generate multiple permutation replicates.\"\"\"\n", " \n", " # Initialize array of replicates: perm_replicates\n", " perm_replicates = np.empty(size)\n", " \n", " for i in range(size):\n", " # Generate permutation sample\n", " perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)\n", " \n", " # Compute the test statistics\n", " perm_replicates[i] = func(perm_sample_1, perm_sample_2)\n", " \n", " return perm_replicates" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Look before you leap: EDA before hypothesis testing\n", "Kleinteich and Gorb (Sci. Rep., 4, 5225, 2014) performed an interesting experiment with South American horned frogs. They held a plate connected to a force transducer, along with a bait fly, in front of them. They then measured the impact force and adhesive force of the frog's tongue when it struck the target.\n", "\n", "Frog A is an adult and Frog B is a juvenile. The researchers measured the impact force of 20 strikes for each frog. In the next exercise, we will test the hypothesis that the two frogs have the same distribution of impact forces. But, remember, it is important to do EDA first! Let's make a bee swarm plot for the data. They are stored in a Pandas data frame, ```df```, where column ```ID``` is the identity of the frog and column ```impact_force``` is the impact force in Newtons (N).\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
dateIDtrial numberimpact force (mN)impact time (ms)impact force / body weightadhesive force (mN)time frog pulls on target (ms)adhesive force / body weightadhesive impulse (N-s)total contact area (mm2)contact area without mucus (mm2)contact area with mucus / contact area without mucuscontact pressure (Pa)adhesive strength (Pa)
02013_02_26I31205461.95-7858841.27-0.290387700.823117-2030
12013_02_26I42527444.08-9832481.59-0.181101940.0724923-9695
22013_03_01I11745342.82-8502111.37-0.15783790.0521020-10239
32013_03_01I21556412.51-45510250.74-0.1703301580.524718-1381
42013_03_01I3493360.80-9744991.57-0.4232452160.122012-3975
\n", "
" ], "text/plain": [ " date ID trial number impact force (mN) impact time (ms) \\\n", "0 2013_02_26 I 3 1205 46 \n", "1 2013_02_26 I 4 2527 44 \n", "2 2013_03_01 I 1 1745 34 \n", "3 2013_03_01 I 2 1556 41 \n", "4 2013_03_01 I 3 493 36 \n", "\n", " impact force / body weight adhesive force (mN) \\\n", "0 1.95 -785 \n", "1 4.08 -983 \n", "2 2.82 -850 \n", "3 2.51 -455 \n", "4 0.80 -974 \n", "\n", " time frog pulls on target (ms) adhesive force / body weight \\\n", "0 884 1.27 \n", "1 248 1.59 \n", "2 211 1.37 \n", "3 1025 0.74 \n", "4 499 1.57 \n", "\n", " adhesive impulse (N-s) total contact area (mm2) \\\n", "0 -0.290 387 \n", "1 -0.181 101 \n", "2 -0.157 83 \n", "3 -0.170 330 \n", "4 -0.423 245 \n", "\n", " contact area without mucus (mm2) \\\n", "0 70 \n", "1 94 \n", "2 79 \n", "3 158 \n", "4 216 \n", "\n", " contact area with mucus / contact area without mucus \\\n", "0 0.82 \n", "1 0.07 \n", "2 0.05 \n", "3 0.52 \n", "4 0.12 \n", "\n", " contact pressure (Pa) adhesive strength (Pa) \n", "0 3117 -2030 \n", "1 24923 -9695 \n", "2 21020 -10239 \n", "3 4718 -1381 \n", "4 2012 -3975 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv('./dataset/frog_tongue.csv', skiprows=14)\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEJCAYAAABohnsfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd2BV9f3/8ee5M3vBTQJZ7CEj7KUGEQkgIIhoFSpUv+1XrbWt1omKu1i1otaitf3VLu1XihIEIaAVGQbZslcgQ5KQvZO7z++PwMUYIBdIcu5N3o9/9HzuuTfvRMyL85mKqqoqQgghhBd0WhcghBDCf0hoCCGE8JqEhhBCCK9JaAghhPCahIYQQgivSWgIIYTwmoSGEEIIrxm0LqC1lZfX4nbLUhQhhPCGTqcQGRl8wdfbfWi43aqEhhBCtBDpnhJCCOE1CQ0hhBBek9AQQgjhNQkNIYQQXmv3A+G+LOd0NWu35WB3uLluaFcG9+ysdUlCCHFREhoaqaq187sPd2O1uwDYm1nC4z8eRu/4CI0rE0KIC5PuKY3sP1nqCQwAFdhxpEi7goQQwgsSGhrpFBbQpK3zedqEEMKXSGhopF9SJFcPjPVc94oL59rkrhpWJIQQzVPa+3GvpaU1Pr0ivLC8DpvdRWJMqNalCCEEOp1Cp04hF3xdBsI1FhMZpHUJQgjhNemeEkII4TUJDY2dLqsj53S11mUIIYRXpHtKI6qq8pfVh9l68DQAPePCeOi2IQSa5T+JEMJ3teqTxttvv820adOYNm0ar7zyCgBPPPEEqampzJw5k5kzZ/L5558DkJGRwYwZM0hNTWXJkiWezzh8+DCzZ89m8uTJPPnkkzidztYsuc0cya3wBAbAibwqNu3N17AiIYRoXquFRkZGBlu2bGHFihWkpaVx8OBBPv/8cw4cOMC//vUvVq5cycqVK5k0aRJWq5WFCxeydOlS1qxZw4EDB9i4cSMAjzzyCIsWLWLdunWoqsqyZctaq+Q2VVZlbdJWep42IYTwJa0WGhaLhccffxyTyYTRaKRnz57k5+eTn5/PwoULmTFjBm+99RZut5t9+/aRlJREQkICBoOBGTNmkJ6eTl5eHlarlSFDhgAwe/Zs0tPTW6vkNjWoZycCzXrPtaLAqH4xGlYkhBDNa7UO9N69e3v+PTs7m7Vr1/LBBx+wfft2nnnmGUJDQ7nnnntYvnw5QUFBWCwWz/3R0dEUFhZSVFTUqN1isVBYWHhJdVxsvrGWLMDL91/LJxsysdqdTB3XjeESGkIIH9fqo67Hjx/nnnvu4dFHH6VHjx788Y9/9Lx25513kpaWxuTJk1EUxdOuqiqKouB2u8/bfil8eXFfqEnHgsl9PNfFxTKLSgihreYW97XqQPiuXbv4yU9+wm9+8xtuvvlmjh49yrp16zyvq6qKwWAgNjaW4uJiT3txcTHR0dFN2ktKSoiOjm7NkoUQQlxEq4VGQUEB999/P6+99hrTpk0DGkLit7/9LZWVlTgcDj766CMmTZpEcnIyWVlZ5OTk4HK5WL16NSkpKcTFxWE2m9m1axcAK1euJCUlpbVKFkII0YxW23vqxRdf5OOPPyYxMdHTdvvtt+N2u/nggw9wOp2kpqby8MMPA7B161YWL16MzWZj/PjxPPHEEyiKwpEjR3jqqaeoqalhwIABLF68GJPJ5HUdvtw9JYQQvqa57inZsFAIIYSHpmMaQggh2hcJDSGEEF6T0BBCCOE1CQ0hhBBek9AQQgjhNQkNIYQQXpPQEEII4TUJDSGEEF6TY+I0VGd1svXgaax2J2MHxBIVFqB1SUIIcVGyIlwjdoeLZ9/fwemyOgCCzAYW/WQE0ZFBGlcmhOjIZEW4j9p3otQTGAB1Nieb9xVoWJEQQjRPQkMj5zsW5FLPChFCiLYmoaGRwT07kxB97hEwNMjI+OSuGlYkhBDNkzENDdkcLnYeKcLmcDGiXzRhQd5v+S6EEK1Btkb34dAQQghfIwPhQgghWoyEhhBCCK9JaAghhPCahIYQQgivSWgIIYTwmoSGEEIIr0loCCGE8JqEhhBCCK9JaAghhPCahIYQQgivSWgIIYTwmoSGxk7kVXIouwyX2611KUII0Sw57lUjblXlD8v3sfdEKQBdOwfz+LxhhAQaNa5MCCEurFWfNN5++22mTZvGtGnTeOWVVwDIyMhgxowZpKamsmTJEs+9hw8fZvbs2UyePJknn3wSp9MJQH5+PvPmzWPKlCncd9991NbWtmbJbeZQVpknMADyS2rZ+G2ehhUJIUTzWi00MjIy2LJlCytWrCAtLY2DBw+yevVqFi5cyNKlS1mzZg0HDhxg48aNADzyyCMsWrSIdevWoaoqy5YtA+C5555j7ty5pKenM3DgQJYuXdpaJbep6jqHV21CCOFLWi00LBYLjz/+OCaTCaPRSM+ePcnOziYpKYmEhAQMBgMzZswgPT2dvLw8rFYrQ4YMAWD27Nmkp6fjcDjYsWMHkydPbtTeHgzu1alRV5RepzB2QKyGFQkhRPNabUyjd+/enn/Pzs5m7dq1/PjHP8ZisXjao6OjKSwspKioqFG7xWKhsLCQ8vJyQkJCMBgMjdrbg+AAI0/NH87nO09hs7tIGdKVpNhQrcsSQoiLavWB8OPHj3PPPffw6KOPotfryc7O9rymqiqKouB2u1EUpUn72X9+3w+vm3OxE6i0ZrGEMqBPjNZlCCGE11o1NHbt2sUvf/lLFi5cyLRp09i+fTvFxcWe14uLi4mOjiY2NrZRe0lJCdHR0URFRVFdXY3L5UKv13vuvxRy3KsQQnhPs+NeCwoKuP/++3nttdeYNm0aAMnJyWRlZZGTk4PL5WL16tWkpKQQFxeH2Wxm165dAKxcuZKUlBSMRiMjRoxgzZo1AKSlpZGSktJaJWtOVVV2Hinioy+Ps+d4cfNvEEKINqaoqtoqfw1/8cUX+fjjj0lMTPS03X777XTr1o3Fixdjs9kYP348TzzxBIqicOTIEZ566ilqamoYMGAAixcvxmQykZeXx+OPP05paSldunTh9ddfJzw83Os6/OlJ4z9fZbL2m1zP9eyUHkwf1027goQQHU5zTxqtFhq+wldCQ1VVcgqrCTQbiIkMOu/rP399EzaHy9MWFmzijQeuacsyhRAdXHOhcdExDbvdzkcffcT69evJyspCr9fTo0cPpkyZws0334zJZGrxgtujWquD1/7vW3JOVwOQktyVn0ztR+apSj7bmo3N4WLCsHgMegXb95ZqGPWyy4sQwrdcMDS2b9/OCy+8wPDhw5k/fz7x8fEYDAZOnTrF5s2bueWWW1i4cCFjx45ty3r90pe78zyBAbBpbz5De3fmnZUHsDsa9pw6klvB9cPi+HJ3w6pwBbjp6m4aVCuEEBd2wdDYsGED//73vwkJafyY0rt3byZMmEBNTQ1vv/22hIYXyqqsTdr2nSj1BMZZOkXh6QUjOJlfRZ+ECBKifXe6sBCiY5IxjTZwOLuM1/7vW85WERxg4GczBvDGf/Y2um/uDb25YURC2xcohBBnXPaYRlpa2kU/eNasWZdfVQfTv1sUv7hlEJv3FhBoNnDjmETiLCFcPyyODXvyUFUY0D2Ka5O7al2qEEJc1AWfNO69997zvmHr1q0oisK3337bqoW1FF940riY8mobDqeL6PPMqBJCiLbWYlNuS0pKeOyxxygqKuL3v/89ffr0abEiW5Ovh4YQQviSFlkRvnHjRmbOnEm3bt34+OOP/SYwhBBCtKxm12n87ne/Y+3atbz00ktMmDChreoSQgjhgy4YGpmZmTz00ENYLBY+/fRTOnfu3JZ1CSGE8EEXHNNITk5GVVXGjh173u3I33333VYvriXImIYQQnjvsqfcPvvss61RT4eXvi2Xr77NI9BkYOa13RnSS57ghBD+Qxb3taGdR4pYmnbAc63XKfzu3rFEhQVoWJUQQpxzRRsWAqxZs4Y333yTqqqqRu1bt2698uo6mIPZZY2uXW6VI7nljBvYRaOKhBDi0jQbGq+++ipPPfVUo3MxxOVJjGl6BnhitJwLLoTwH82GRlxcHBMnTmyLWtq9awd34UReJd8cLMRo1DGkVyfeXrEfm8PF9UPjmHF1d61LFEKIi2p2TGP58uWcOHGClJQUDIZzGTNy5MhWL64l+NKYxllWu5PSSiuL/rqd7//0fz5rICP6XdoZ6EII0ZKueExj27ZtbNq0iS1btjRqX7Vq1ZVX10EFmAyczK/ih3F9JLdcQkMI4dOaDY1Dhw6xadMmzGZzW9TTYXTvEuZVmxBC+JJm957q3LkzTqezLWrpUOKjQ5h7Q2+CAwwY9DpuGB7P2IGxWpclhBAX1eyTRkxMDDNnzmTcuHGNzgR/6qmnWrWwjuCGEQlcPzweVVXR6+Q8cCGE72s2NBITE2W6bSvSKQqcZ5sWcWlUtxNbxoc4jmegBIYTMOZ2DN2Gal2WEO3OBWdPlZWVERUVddE3l5aW0qlTp1YprKX44uwp0fLse9di2/bRuQa9keB5r6MLkHUwQlyKyz5PY+HChbz//vtUVlY2ea2mpoa//OUvPP744y1TpRBXyHX62A8aHLiLszWpRYj27ILdU0uXLuWvf/0r06dPp3v37iQlJeF2u8nNzSUrK4v58+ezdOnStqxViAvSxfSGnD3nGvQGdJ2TtCtIiHaq2cV9VquVb775hpMnT6IoCt27d28yKO7L/K17qqC0lrTNWZRVWRnVP4YbRsSfd2t60ZjqcmL7+h+eMQ3z2Dswdh+hdVlC+J0WOyPcX/lCaLjdKss2ZLJ5Xz7BAUbmXNeTUf1jOJxTzuqMbGwOFxOGxjH6qhgee3cr5dU2z3vnT+7LdUPjNKzef7jK83Ce2I4SFI6x99UoRllbJMSluuIV4eLKbdqbz/od3wFQb3Px51WH6BwWwJJle3G63ACczK+i1upsFBgAu48VS2h4wVWYSd2ql8HdsKbIcWwLQTOflqc0IVpYqy8OqKmpYfr06Zw6dQqAJ554gtTUVGbOnMnMmTP5/PPPAcjIyGDGjBmkpqayZMkSz/sPHz7M7NmzmTx5Mk8++aRfLjQ8fqrxZAKXW+XrAwWewDjrVHENel3jX3IxkUGtXl97YD/0pScwANxFJ3EXZmpYkRDtU6uGxt69e7njjjvIzs72tB04cIB//etfrFy5kpUrVzJp0iSsVisLFy5k6dKlrFmzhgMHDrBx40YAHnnkERYtWsS6detQVZVly5a1Zsmtoldc4+1BdIrCVd2aTmfuFhvKrRN6YdA3/GdJjAlh+jgZzPWGotM3bTxfmxDiijQbGrW1tTz33HMsWLCAiooKFi1aRG1trVcfvmzZMp555hmioxs24auvryc/P5+FCxcyY8YM3nrrLdxuN/v27SMpKYmEhAQMBgMzZswgPT2dvLw8rFYrQ4YMAWD27Nmkp6dfwberjfFD4pg4PB6TUUdkqJn/md6f4X2jmTI60fNkMbR3Z64d3IXUkQkseeBqXvrZaJ69axThIdIv7w3jwFQwnjsBUR8/EH10Dw0rEqJ9anZM48UXXyQ6OprS0lLMZjM1NTUsWrSI3//+981++EsvvdTouqSkhDFjxvDMM88QGhrKPffcw/LlywkKCsJisXjui46OprCwkKKiokbtFouFwsLCS/n+fIJOpzBvUh/mTerTqP22Cb24cUwSDqebyNBz4RAcYCQ4wNjWZfo1facEgm9bjDN7F0pQBIYkWQ0uRGtoNjQOHz7M4sWL2bhxI4GBgbz22mtMnz79sr5YQkICf/zjHz3Xd955J2lpaUyePLnRgKWqqiiKgtvtPm/7pbjYLABfYGn+FuEtSyjOzmHozIHoZOaUEK2i2dDQ/WAjPZfL1aTNW0ePHiU7O5vJkycDDSFgMBiIjY2luLjYc19xcTHR0dFN2ktKSjxdXd7yhSm3ovWp9nrqv/gjrlMHwBiAeeQcTANv0LosIfzOZW8jctbIkSN59dVXsVqtbN68mV/84heMHj36sopRVZXf/va3VFZW4nA4+Oijj5g0aRLJyclkZWWRk5ODy+Vi9erVpKSkEBcXh9lsZteuXQCsXLmSlJSUy/raon2z713TEBgADiu2rR/irinVtigh2qFmQ+Phhx8mKCiI0NBQlixZQr9+/Xj00Ucv64v169eP//3f/+WOO+5g2rRp9O/fn+nTp2M2m3n55Zd54IEHuPHGG+nRowdTpkwB4LXXXmPx4sVMmTKFuro65s+ff1lfW7Rv7vL8xg2qG3dFgTbFCNGOebUifMeOHYwcOZKKigp27tzJDTf4z2O/dE91DPbDX2Hb/LdzDaYgQub+HsUUqFlNQvijK+6eWrJkCW+99RbQsA/Ve++9JxsVCp9j7Dce06hb0UV2RR93FUE3/kYCQ4hW0OyTxvTp01mxYgVGY8MUULvdzuzZs1m9enWbFHil5ElDCCG8d8VPGg6HwxMYAEajUfbzEUKIDqrZKbfDhg3jN7/5DXPmzEFRFNLS0khOTm6L2oQQQviYZrun6urqeOutt8jIyMBgMDB27Fh+8YtfEBjoH/3F0j0lhBDeu+Kt0d955x051lUIIQTgxZjGV1991QZliO/LK6nldFmd1mUIIUQTzT5pxMfHc/fddzNs2DCCg4M97XfddVerFtYROZxu/vDJPg6cLANgRF8L984ciE4nEw+EEL6h2dCIiIgAIC8vr9WL6WjcqsrhnHJsdheDekSx7VCRJzAAdh4tZm9mCUP7yLaGQgjf0GxoLF68GGgIDafTSVKSHArUEtxuldf+bw9HcisAsEQEMOw84VBSaW3r0oQQ4oKaDY2cnBx+/vOfU1RUhNvtJjIykj/96U/07NmzLeprtw5klXoCA6C4worD6UavU3Cdme1lNOgY0ruzViUKIUQTzQ6EP//88/z0pz9lx44d7Nq1i/vuu4/nnnuuLWprdxxONyfyK6mpd1BnbXrWucmg58HbkhnSqzPD+1h45PahWCL8Y2qzEKJjaPZJo7S0lJtvvtlzfcstt/C3v/2tNWtql3JOV7PkP3upqrVjNOi4/fpehIeYqKyxA2DQ6xg3MJb46JDznh8uhBC+oNnQcLlcVFRUeAbEy8rKmnmHOJ/lX2VSVdsQEA6nm483nmTRT0aweV8BVruLawd3IT7at08ZFEKIZkPjxz/+MT/60Y+YOnUqiqKwZs0aFixY0Ba1tSslVbZG13U2J/W2hi4qmVArhPAXF9xGZO/evZ49prZu3cqWLVtwu91ce+21jBs3rk2LvBK+so3IJ5tOsjoj23PdOz6c4op6Ks50TxkNOp5eMIJ4izxtCCG009w2IhcMjZtvvpkVK1awYMEC/v73v7daga3NV0LD7VZZtyOXAyfLiLME07VTEP9Yd6zRPVNGJXLb9b00qtC/qaob+66VOI5noASFYx51K4YufbUuSwi/c9l7TzmdTu6++24OHTrEvffe2+T1d999t2Uq7CB0OoWpo5OYOrphncu+EyVN7gkMaLa3UFyA49AG7LtXAqBWF1OfvoSQea+jmII0rkyI9uWCv6X+/Oc/880335CVlcXkyZPbsqYOYWD3TvRLjGi0uG98cleNq/JfrryDjRscVlyFJzAkDNKmICHaqWa3Rt+2bRujR49uq3panK90T52PW1U5nF2OzdGwjYjRoNe6JL9l25WGfVfauQZFT/Dc19AFR2pXlBB+6LLHNNoLXw4N0XJUhw3rhj/hzN4D5iDMo2/D1G+81mUJ4XckNCQ0OhTVXg8GI4pOxoeEuBxXfEb4+TKlsrLyyqoSopUopkAJDCFaUbOhMXv27CZt8+bNa5VihBBC+LYL/pVswYIF7N+/H6vVyrBhwzztbrebQYNkRooQQnREFxzTqKmpoaKigoULF3rO1AAwGAxYLBZ0umYfUnyCjGkIIYT3LntMIyQkhPj4eJYuXcrq1auJi4sD4C9/+QtWqxwMJIQQHVGzjwtPPPEEFRUNC9DCwsJQFIWnn3661QsTQgjhe5oNjezsbB577DEAQkNDWbhwIcePH2/1woQQQvieZkPD6XRSU1Pjua6trT3vNNzzqampYfr06Zw6dQqAjIwMZsyYQWpqKkuWLPHcd/jwYWbPns3kyZN58skncTobtgzPz89n3rx5TJkyhfvuu4/a2tpL+uaEEEK0rGZDY9asWdx66628+eabvPXWW9x+++3nnYb7Q3v37uWOO+4gOzsbAKvVysKFC1m6dClr1qzhwIEDbNy4EYBHHnmERYsWsW7dOlRVZdmyZQA899xzzJ07l/T0dAYOHMjSpUuv4FsVQghxpZoNjXvuuYeHH36Y6upq6urqePjhh7nrrrua/eBly5bxzDPPEB0dDcC+fftISkoiISEBg8HAjBkzSE9PJy8vD6vVypAhQ4CGdSHp6ek4HA527Njh2SzxbLsQQgjteLV0duLEiUycOBFoWCGenZ1Nt27dLvqel156qdF1UVERFovFcx0dHU1hYWGTdovFQmFhIeXl5YSEhGAwGBq1X6qLTR0TQghxaZoNjX//+9+8+uqr1NfXe9qioqL4+uuvL+kLud1uFOXcwaaqqqIoygXbz/7z+3547Q1Zp9FxqC4nrtPHUALD0UfFaV2OEH7psg9hOuvPf/4z77//Pu+88w6//vWv2bBhA6dPn77kQmJjYykuLvZcFxcXEx0d3aS9pKSE6OhooqKiqK6uxuVyodfrPfcLcT7umlLqPv0tak0pAMarrifgmvkaVyVE+9PsmEZERATJycn079+f0tJS7rvvPnbs2HHJXyg5OZmsrCxycnJwuVysXr2alJQU4uLiMJvN7Nq1C4CVK1eSkpKC0WhkxIgRrFmzBoC0tDRSUlIu+euKjsG+d60nMAAch77EVZ6nYUVCtE/NhobBYKCyspKkpCT27dsHgMvluuQvZDabefnll3nggQe48cYb6dGjB1OmTAHgtddeY/HixUyZMoW6ujrmz2/4G+IzzzzDsmXLuPHGG9m5cye//vWvL/nrio5BrW+687JaJ7sxC9HSmj1PY/ny5Sxfvpx3332XWbNmERUVRWxsrN9Mf5UxjY7BmfMt9eve8FwroRaCf7RYtkkX4hK1yCFMdXV1BAUFUVhYyP79+7nmmmsICAho0UJbi4RGx+HM/RbHsQyUoAhMg6egC4nSuiQh/M4VD4S7XC5WrFjBli1b0Ov1TJgwwW8CQ3QshsQhGBKHaF2GEO1as6HxwgsvcOLECWbOnImqqnz88cfk5uby4IMPtkV9QgghfEiz3VOpqal89tlnGI1GAGw2GzfddBPr1q1rkwKvlHRPCSGE9674jPCoqKhGs6UURSEsLKxlqhNCCOFXmu2e6tevH3PnzmX27Nno9XrWrFlDZGQk77//PoBX+1AJIURH5rbZqNq2FVdlJaEjR2GK7aJ1SZet2dCw2Wz07duXgwcPAhAfHw/AsWPHWrcyIUSb+SJ3I9tP7ybMFMr0Hql0C0vUuqR2Q1VVTv3+FawnTwBQ9tkqEh5bSED3HhpXdnm8mnLrz2RMQ4iLy8jfwQdH/uO5DjQE8uK4JwgwyCzJllB//Djf/a7xBq5h464m9u6faVTRxV3xlNtt27bx3nvvUVnZeHXt8uXLr7y6Dq6m3kHG/gKsDhdjB8RiiQjUuiTRAR0sPdzout5Zz4mKbKrs1WRV5dA9vBujY4ehU5odAhXnozvPRqs6fdvX0UKaDY2nnnqKO++8k8REeVxtSXaHixf/sZOi8obdg9dtz2XRgpHERAVpXJnoaLoEx/Bt8QHPtYLCt8UHyCjYDsDX+ds5XVvIzb2maVWiXwvs2Yug/ldRd/gQAIo5gMiJkzSu6vI1GxqdOnXy7AUlWs63mSWewACot7nYvK+AOdf11LCq9kN1WLHvWY2r+CT62L6YhtyIojdqXZZPmpiYwsnKHI6WZ2LSGZnRcwqrTzaeUr8lb5uExhWI+9VD1OzZjbOykpBhwzFG+e9uBc2GxvXXX88HH3zAtdde6zkQCaBr166tWlh7Z9A3fdQ36C/9vJCOzF1TCijn3S7E+tVfcGbtBMCVdwi1vlK2Sr+AQEMgvxz6v1TYKgnQBxBgMPPf3E3YXHbPPQEGs4YV+j/FYCB05City2gRzYZGeXk5r7/+OoGB5/rbFUVh9+7drVpYeze4ZyeSYkLJKawGICzYREqyBLE3VLcL64b3cJ7YBoCh11gCrvsZik535nUnzuxdjd7jOLFNQqMZEeZwz7/P6DGZfx3+DyoqCgo39ZiiYWX+xW21UrOn4c9fyLAR6MztK3CbDY0NGzawZcsWOnfu3Bb1dBgGvY6Fdw5j17FibHYXw/tGExIo3SfecGbt9AQGgDNzK87uwzF2H4G7qgjFHIwSFIlaW+a5RxfSSYtS/daYLiPoEd6N7KpcuoUlEh0k//97w1VXS+6Lz+Moajia2rh6FYlPPYOropzSz1bhqqwibNw4wsZerXGll8+rMY0oP+5/82VGg54xV8VqXYbfcVc2PSveXfodtfvScRdmgt6AocconFm7wGkDUxDmsXdoUKl/iw7qLGFxiaq3feMJDABH4Wmqtn5N2WercJ2ZgVp3+CCK0UjoCP/srmo2NPr06cPcuXOZMGECJpPJ0y4rwYVWDElDse9KA9Xd0KDoUevKGwIDwOXEmbmVoJufA6cNXadEFGP76iIQvkl1Opu02fPzPYFxVvXOHe03NKxWK927dyc7O7sNyhGiefpOCQRO/hX2/esABdPgKTgOf9X4JlVFtVZhiB+oRYk+z+FycLQ8k2BjEN3Dk7Qup90IHTWGsjWf4aquAkAfFkbYmHFUfvVlo/uMnfz3CU5WhIt2wXHsa6xf/dlzrQSEEjz3NRSZ9dNEha2S3+9aSpm1HIChlkH8dNCdGlfVfjjKy6n6ejMoCuFXX4MhIpKST5ZTtvYzUFXMCQnEPfgIBh/d+PWyT+771a9+xZtvvsmMGTPO+8ZVq1a1TIWtTEKj47Af/ALHsQx0QeGYRtyMvpMsSD2ftMw1fJ77VaO23wy/nx7yxNGqHOXluKqrMCckoii+O73+srcR+dnPGvZFefrpp1u+KiGukKs8D+fxDDAFYeo3HiUgBNOAGzANuEHr0nxeraPuPJQ432QAABkOSURBVG21V/SZLnfD8Qn6M9tjuFU3B0uPUGWvZnDnAYSaLvxLqKMwRkZijIzUuowrJt1Twu+4SnOpS3sBXA4AlPAYgue8KCu+vZRZkcWbe/6E+8xEgkhzBJOTJlBpr2Zo9CDiQhq27a6yV5NbdYrEsHjCTKEA5NUUsKdoP1EBkYyMHYpB0bPixGdsOpWBTtExKXECU7tP5J29f+VA6RGgYfHgb4b/nC7BMdp8w35Edbmo3r4Ne2EBwYOHEtij7XfCvezuqfZCQqP9sX79LxwHv2jUFpD6S9wVp3Eez0AJCsc88hb00f659XRbyKzI4puCnQQbg8iqzOFEZTYAOkXHA0N+Sq2jnr8d/BCn6sKg6LlrwFwCDYG8vfcvnrDpH9WHa+PG8N7+fzT67Dv63sK/j37cqO3qrqOZ2++WNvnefIEtLw/VYcec1O2CXVF1Rw5T+mka7vo6wlMmEDHhegree4fq7WfWICkKXe69n9DhI9qw8hbY5VYIX6MYTE3aXAVHcew/s19S+SnqirMImfc6ilG29z6fXhHd6RXRnbyaAr7I3ehpd6tuvjqVwanqPJxqQ5eTU3WxIvMzuoZ08QQGwOGyY0QFNO1uyavJb9LmcDta4bvwPaqqUvCnd6jZ2bDZY0CPnsQ/9Ai6gADc1npc9VaMkZE4KyvIe/N1VEfDz6Xog3+gGPVU79j+/Q+j4ov1bR4azZHQEH7HeNX1OI5tQa0/M62xS1/cVcWNb7LX4SrMlCm3Z5TWl7HsWBo51afoE9GT2/rOIsQYjP48253rFB1V9ppGbVWOGuLPc2/3sES+zj+3Ol9BYWyXkWRV5fJddR4AekXPtXFjWvg78k11hw56AgPAevIElVs2gdtNSdonqHY7gf36Ezp6jCcwzqo9fPiHHwd639tCXUJDY2VVVuxON7GyJbrXdKGdCb5tMc6cb1FMQegTB2PflYYr53v7oSkKugj/PVKzpb1/8EOyqnIB2FW0F4C7B84jNjiGZMtA9p7ZGt2oMzIx4VqCjUFsyfvG8/4xscMZGTuMg6WHcbgbFrANsQxibNeR2Fx2vvxuE3pFT2q360kMi+dXQ+/hm4KdVNmrGR6dTHxox9hXzVle1qTNnp9H5eZNcGYkoP7IYYzR0U3uC+zWHb3ZTOWmM09+ej1RU25s1Xovh4SGhv6x7igb9+ShAv2TIvnlnMGYjb73NwtfpJiDMfY5t3+PKXkqrqKTuPIOgsGMedQcz35Tqq0WTIEoHfQQIbvL7gmMs46WZ+JW3WRWnGRiQgqjY4dRZq1gcOer6BQYRWJoPDGBnTlZmUOPiG6MjxuHXqfn6dEPs6/kEFEBkQzq3B+AlPixxIV0Qa/T0T2sYdpuoCGACQnXtPn3qrXgwUPQBQbirj9z7IFOhykuzhMYZ7lra+k082bKPluF6nQSPHQY4ddNQDEYCRk6HPvpAoIHDfbJs8QlNDRy7LsKvtqT57k+nFPOpr35TBqRoGFV/kNV3bjLC9CFRKGYAlFMQQRNewR3XQWKMRDFaMZdU0b9F3/EXXQCJTiKgPF3d8juKpPeREyQhcK6c114XYJjeXnHm+TVFAAwsFM/7hn8E8/pfIqiYNKbMOvNmHUmFEXBrbr55vQudhftIyoggk4BkViCOvPmnj+RU/UdAH0ienL/kP/BoOuYv1oMYWEkPPoEZevTUe12Iq67noDu3SldmYa77txU5+DByYRffS2Rk1Jx2x2NFvoFDxpM8KDBWpTvFf2zzz77rNZFtKb6evsPQ94nHM4pZ8/xkkZtsZFBDOopu7E2x115mrq0F3Hs+RT7gc/RBUd5FvKptjoUvQFFb8S66X1cp86cSOeox3XqAMZBkzrkE0e3sESOl5+k1llHYmgcvSK6Nzqtr6i+hG5hCUQHWQD4OHMVq06u41RNPvtLD2Nz2ThdW8SnJ9OpcdRSXF/KvpKDBBoC+Dr/XB9+qbWcrsExdA3pWBtxum02bN/logsKwhjVidBhwwkdOQqjxYJiMBLUrz+OslJ0AWaiptxIxISJACgGo89tna4oCkFBTSebnKXJXwfuvPNOysrKPIc6Pf/88+Tm5vLOO+/gdDpZsGAB8+bNAyAjI4PFixdjs9mYOnUqDz74oBYlt7hBPaIwG/XYHA0zVBRgeF+LtkX5Cdv25ahVZ3YSddqxfv1P9F37Y93wHq6CI6A3YR55C+7Sxl0yan0Val0lSgfcJj0pLIFFYx7B6rIRaAhgReZnTe6p/t4Cv+8HAUBG/nYSwxo/BVfZqzl1nplS1fYrWyjob2oP7KfgT0tx19ejCwqm689/QVC//qhuN/bTpzF26kRA9x7EP/iw1qW2iDYPDVVVyc7OZsOGDZ7QKCws5MEHH+STTz7BZDJx++23M3r0aOLj41m4cCH//Oc/6dKlC/fccw8bN25k/PjxbV32Fau3OfnmUCE2u4vRV8UQGWrm4TuGsGZrDjaHiwlD4+ib6P+rRduCu6qocYPDin3PqobAAHDZsW37Pwy9xuGuPO25TRfRBSW4427zrygKgYaGKcgjYoay4bstuM5Mqw02BNEtNAGby45ZbyLIEIj9eyf3BRoC6RIcw7HyTE+bQWdgbJdRbD+92zM4HqA3MyS6Y3UBFn3wT88YhruulqJ/f0CXe+4j/603cJQUowsMJOYn/+NzU2cvV5t3T508eZK0tDQ2bdrE+++/j8vlorCwEJ1Ox7Rp0zAajZSVlZGZmYnBYODYsWMsWLAAnU6Hqqps3LiRG27wfqsIX+iesjtcPP+3nWQcOM3B7DK+3l/AqP7RxFtCGH1VDOMGdqFr52Bti/Qjan3VuYAAdFEJKHoD7oqCRveZRtyMYgpEra9CH9OTgOt+ii7QNzeJa2vh5lD6RvXGrbpJCk1AReXTk2vZ8N1mzHozgzpfxf6SQ6io6BQdP+ozi3FdRpJVmUO5rYIAfQC39pnJYMtVXBXVt+FzwhKY2+8WLB3oDA5VVSn+z0eNBrrddjv206exZWc13ON0Un/0MJE3pHpOl/RlPtc9VVVVxdixY3n66adxOBzMnz+fqVOnYrGc65qJjo5m3759FBUVNWkvLGx6AM/FXGxlY1vZ/G0ep8vODYLVWp3syiylU3ggH6Yfod7u5MZx3blr+lU+vZGZr1BT76AiNJC6Yzswdo4javwd1J3cS8n3jnjVBYYQmzwK3Uj/eyptKxbLQEb3GshH+1eRUdDQHWV3O/g4cxVvT3uBkT0GklmaRe9O3bEEN3TpLY57jLL6CkKMQZjOLLK0WPozvGd/zb4PrZWPHUPJlq8915ZrxlF18FCje1zV1UQE6jBFhHranHX1nF6zlvr8AjqNGUXUqJFtVvOVaPPQGDp0KEOHDvVcz5kzh8WLF3Pfffd52lRVbZit4XY3+iV6tv1S+MI2IlVV9U3aCktq+OjzY57rFV9l0jnUxNgBHWsA8bL1ScXUJxWAchsQNwrzmDIcZ7YR0Xfpx3f/fB70RkxDpmGIu0rben3YieIfjP2oKruyDrG3+CAnK3PoGdGNOb1vItQUwtd521ib/V9cqosJCdeQmjRBo6p9R8Ttd+IKDsN68gSBvfsQNv0mHLqPsZ4+1zUa0L0Huf/djLOqktARozDFxvLdqy9Tf7Thibnov18Sc9dPCb9a+2nKPreNyM6dO3E4HIwdOxZo+AMaFxdHcfG56YDFxcVER0cTGxt73nZ/M7R3Z7p2Dia/pGGAMCTQSOfwwCb3ncirlNC4TO6aMgw9R2MaPAVXYSZ1K18CziymKjhK8G2/RRfmf3922sKATv0azaQKNASSUbCDI2XHAdhZ+C1Wp41pPSbx4ff2lFp5Yi1dg2MZ2LnjPmUA6AICsNz6o0ZtnWbPQTEaqT14AFNcPPa8UxT+430Ayj5bRezP7vUExllVWzb5RGg0p8072Kqrq3nllVew2WzU1NSwYsUKXn31VbZu3UpZWRn19fWsX7+elJQUkpOTycrKIicnB5fLxerVq0lJSWnrkq+Y0aDnqfnD+cnUftxxQ29e+OlohvbuzA+fmfokRGhSnz9T3S7q//sutR8+RO0HD1G/4c84sndzNjAAcDtxfrdPsxp93dguI5nV80a6BsfSL7I39yf/D8fKTzS653DZMU5UZDd5b2ZFVhtV6V90RiOdZ88h6elnibh2PLacbM9rqsPRsCnhD7YI0QX5x64Qbf6kMWHCBPbu3cusWbNwu93MnTuX4cOH8+CDDzJ//nwcDgdz5sxh8OCGxS0vv/wyDzzwADabjfHjxzNlypS2LrlFBJgMpCSf20ohPNjE3dP6k7b5JDaHmwlD4xjVX7aOvlTOrJ04T5zd7kLFefxrDAMmNrlPFy5PcOezp2g/R8qOER/alcdH/spzHkZcSBfP3lFnr7uFNT3Uqlu4HHTVLF3TLnVdQACRqVMoX9sw9VkXGEin6Te1dWWXRbZGF37NtvtT7Ds/adRmHHELamkOzqydoCgY+12H+Zr5MsngB77I3dhovcbo2OHMv6qhm+W76jz+euADiupLsAR2JtkygGp7DS7VxaHSozhVF9fFX83MnlO1Kt+nuW02ytenY806SWDfftTu2+vpjlLMASQ+/iTmhARs3+ViLywkqF9/9CHaT9oBOU9DQqOdc5XkULfi2XNTHhU9gZN/BTodSkgnFFMguiDp9jufZ7f+juL6Us+1TtHx6jXPklWVS2FdMVdF9SXQGMDyY6vYWbTHc9/07qlM6TZRQvgiGp2LAYRffwNBvXrjrKokZNhwjFG+u8DU5wbChWhJ+s5JBKb+Cvv+daAoKIER1KcvAVSUwDACpz0GEhrnFWBofNaIUWdgeeYqthbsAMCg6PnpwDvZXby30X1bC3Ywtbscq3shqtNJ9c4djdpqdmwjZu6PNaqoZfn+ShMhmmFIGkLQ9McIuP5enCe2cXYQXK2vwr57pbbF+bBp3SehV84Nxl6fcC3fFOz0XDtVFxu+24JZ33ihV5DRPwZsNaPXow9tvIjUEBGBs7ICW94pjYpqOfKkIdoN1VoDZ7bF8LTVVWhUje8b1Pkqnh37KMfLTxIf2pUgQyBrs//b6B4XLqb3mMzyY5+iomLQGZjRY7JGFfsHRVGIvn0up//fe6hOJ4o5AFNcPCcfeQjcbsyJScQ/+DD60NDmP8wHSWiIdkF1uxrGLyzdcRefmwZq+N6ZG6KpqIBIRncZ7rkeYhnEt8X7gYZT+CbEX8OQ6EH0j+rDqep8ekV0J9wsW7E0J3TkKAL79cN+6hS6oCByX3jW85otN4eydWuxzLlNuwKvgISG8HvOvENYN7yHWleBEh6Dsd94VFsthu4jMPbqGMeMtpS7BtzBzsL+fFedR0FtIcuPr2J30T5u7TOT4THJWpfnVwyhYRj6X0XtwQNNXnOUFJ/nHf5BQkP4NVV1Y/3qL55uKLWyEHeohaAb28c21G3NoDMwpssIdhft4+iZHW13FVVgc9m5L/kujavzT4G9+6APDcNVXeVpCx3mvzveSmgI/2avR61tfC6zuzwP53f7sO1aCQ4rxv4TMA2U2T6X4nDZsYteC+/pTCbiH3mMstWf4qysJGzs1YSOGq11WZdNQkP4NcUcjC66J+6ic9te6GP7UL/uLThzxoMt41/oQjph6Db0Qh8jfiA+pAu5P1gRLi6fuWscsXf/DLfdjt5Ptgu5EAkN4ZccJ7fjOLQBxRiAacg0nJlbcZXkYogfgC4q4czU23Ocp/ZLaFyCuf1u5a8H/kVRfQnRQZ2Z12+O1iX5tYqNX1Gy/CPcVivByUPo8rN7fe6YV29JaAi/48w7hPWLpY2ug29/BV1QONCwSvyHdFHxbVZfe5AQ2pVFYx6hxlFLiDFYVn9fAUdZGUUf/APcbgBqv91D+fp0Os2YqXFll0dCQ/gd5/cOW2posOHM3Yu76ITnPA1DzzE4s3eD24Ghx2iMff1vd2StKYpCqMk39kPyZ/b8U57AOMv2Xe4F7vZ9EhrC75zvXAxXcTbOIxsBUKtLcNZVEHTrYnTmIBSzHKUrtBPQoyeKOQDVZvW0BQ3w33PUZRsR4XeM/a5Df/YkPkXBOGAi1P9g5bfLiVp5WgJDaE4fFEzcrx4koGcvjBYLnW6aRXjKdVqXddlkl1vht9wVp8FoRhcciX3vWmzbPjr3ot5I8LzX0QX451YNov1y1ddTvf0bVJud0NGjMYRHoLrd1OzZhf30aUIGJ2NO0O6cEtkaXUKjQ1DdTmwZH+I49jVKUAQBY26X2VLC57jtdnKffwb76QIA9CGhJD79LKUrP6Eq4+uGm3Q6uv78AUKGaPPnV7ZGFx2CojMQcM18Aq6Zr3UpQlC9exfWE8cJ6Nmb0GHn9vaq/XaPJzAAXDXVlH+xnqqtGefe7HZTvj5ds9BojoSGEEK0oJKVKyhbdXZL/nRsM2bSeebNmtbUkmQgXAghWlDFf79ofP3FegCsuTkoJgPG2HOr6/UhoUTekErY2HHn3qDTEZk6pU1qvRzypCGEEC1IMRp+cG3k9N/+H1VbNgOgj4yk0+w56IxGQkc1DITH/OR/CB4yDMfpAoIHJWNOSNCidK/Ik4YQQrSgTtMbr/QOvybFExgArvJyXJWVRE6ajCG84ShiRacjdNhwom6c7tOBAfKkIYQQLSpiwvUE9OyJNfM4AT174a6ro2zN6kb3fH+bdH8joSGEEC0sIDGJgMQkAFSnE6PFgqP4zMFLikLYOP89UVLWaQghRCtzVpRT/vk6z3kawT68jYgs7pPQEEIIrzUXGjIQLoQQwmsypiGEEC3Iba2nevt23A47oSNGYQgP17qkFiXdU0II0ULcdju5Lz6LPT8fAH1oKIlPP4cxKoq6Y0dxVVURPHAguoBAbQu9iHax99SqVat45513cDqdLFiwgHnz5mldkhBCNFH77R5PYAC4qqup2rIJe0E+1Tu2A6APDSPhiacwRTc9F8Yf+PyYRmFhIUuWLOHDDz8kLS2Njz76iMzMTK3LEkKIJlSa9mo4Kio8gQENazTKP09vy7JalM+HRkZGBmPGjCEiIoKgoCAmT55Merr//sCFEO1XSPJQjLGxnmtdSAhB/fo3uc9dV9eWZbUon++eKioqwmKxeK6jo6PZt2+fhhUJIcT56cxmEhcuonrbN7jtNsJGj0EfGkZZl67YC850WykK4deO17bQK+DzoeF2u1EUxXOtqmqj6+ZcbEBHCCFaXiixSY33n4p85SUK1qTjKC/Hct14wvr306i2K+fzoREbG8vOnTs918XFxURfwgCSzJ4SQmhPIXDiVAIBG1BcXK11QRfk94v7xo0bx9atWykrK6O+vp7169eTkpKidVlCCNEh+fyTRkxMDA8++CDz58/H4XAwZ84cBg8erHVZQgjRIcniPiGEEB5+3z0lhBDCd0hoCCGE8JqEhhBCCK9JaAghhPCaz8+eulI6nfcLAYUQoqNr7ndmu589JYQQouVI95QQQgivSWgIIYTwmoSGEEIIr0loCCGE8JqEhhBCCK9JaAghhPCahIYQQgivSWgIIYTwmoSGEEIIr0lo+IC+fftqXUK78P2fo/xML9+Ffo7yM710Z39md9xxB5999lmj1+rq6hg9ejRlZWValHbZJDSEEKKV3XLLLaxatapR2/r16xk9ejRRUVEaVXV5JDSEEKKVTZ06ld27d1NRUeFp+/TTT7nllls0rOrySGgIIUQrCw4OZuLEiaSnpwNQWFhIVlYW11xzjcaVXToJDSGEaAOzZ89m9erVAKxatYqbbroJvV6vcVWXTkJDCCHawMiRIykuLqagoMBvu6ZAQkMIIdrMrFmzeOeddwgPDycxMVHrci6LhIYQQrSR2bNn8/HHH/vtUwbIyX1CCCEugTxpCCGE8JqEhhBCCK9JaAghhPCahIYQQgivSWgIIYTwmoSGEC3smWee4frrr2fJkiValyJEi5Mpt0K0sH79+vHVV18RGxurdSlCtDgJDSFa0Ny5c9m1axd9+vQhMzOTyZMnc/ToUR566CG6devG888/T0VFBYqicPfddzNr1iwA3nvvPZYvX05wcDAjRozgv//9L19++aXG340QTRm0LkCI9uTDDz+kb9++/P3vf2fOnDn07t2bN954A6fTyZQpU3j00UdJTU2lsLCQW2+9laSkJGpqavjkk09Yvnw5oaGhPPnkk1p/G0JckIxpCNGKRowYAUB2djY2m43U1FQAYmJiSE1NZfPmzWzcuJEpU6YQFhaGoijMmzdPy5KFuCgJDSFaUVBQEAAulwtFURq9pqoqTqcTg8HA93uJ/XG7bNFxSGgI0QZ69OiBwWBg/fr1QMMhPOvWrWPcuHGMHz+e9evXU11dDcDy5cu1LFWIi5IxDSHagNFoZOnSpbz44ov84Q9/wOVycf/99zNmzBgAbrvtNn70ox8REBBA7969CQwM1LhiIc5PZk8JobH9+/ezZ88e5s+fD8D777/P3r17eeONNzSuTIimJDSE0FhNTQ0LFy7k5MmTKIpCly5deOGFF4iJidG6NCGakNAQQgjhNRkIF0II4TUJDSGEEF6T0BBCCOE1CQ0hhBBek9AQQgjhNQkNIYQQXvv/vFa+oCnWosQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Make bee swarm plot\n", "_ = sns.swarmplot(x='ID', y='impact force (mN)', data=df)\n", "\n", "# Label axes\n", "_ = plt.xlabel('frog')\n", "_ = plt.ylabel('impact force (mN)')\n", "plt.savefig('../images/frog-swarmplot.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Permutation test on frog data\n", "The average strike force of Frog A was 0.71 Newtons (N), and that of Frog B was 0.42 N for a difference of 0.29 N. It is possible the frogs strike with the same force and this observed difference was by chance. You will compute the probability of getting at least a 0.29 N difference in mean strike force under the hypothesis that the distributions of strike forces for the two frogs are identical. We use a permutation test with a test statistic of the difference of means to test this hypothesis.\n", "\n", "For your convenience, the data has been stored in the arrays ```force_a``` and ```force_b```." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "force_a = np.array(df[df['ID'] == 'II']['impact force (mN)'])\n", "force_b = np.array(df[df['ID'] == 'IV']['impact force (mN)'])" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value = 0.0049\n" ] } ], "source": [ "def diff_of_means(data_1, data_2):\n", " \"\"\"Difference in means of two arrays.\"\"\"\n", " \n", " # The difference of means of data_1, data_2: diff\n", " diff = np.mean(data_1) - np.mean(data_2)\n", " \n", " return diff\n", "\n", "# Compute difference of mean impact force from experiment: empirical_diff_means\n", "empirical_diff_means = diff_of_means(force_a, force_b)\n", "\n", "# Draw 10,000 permutation replicates: perm_replicates\n", "perm_replicates = draw_perm_reps(force_a, force_b, diff_of_means, size=10000)\n", "\n", "# Compute p-value: p\n", "p = np.sum(perm_replicates >= empirical_diff_means) / len(perm_replicates)\n", "\n", "# Print the result\n", "print('p-value =', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bootstrap hypothesis tests\n", "- Pipeline for hypothesis testing\n", " - Clearly state the null hypothesis\n", " - Define your test statistics\n", " - Generate many sets of simulated data assuming the null hypothesis is true\n", " - Compute the test statistics for each simulated data set\n", " - The p-value is the fraction of your simulated data sets for which the test statistic is at least as extreme as for the real data\n", "- One sample test\n", " - Compare one set of data to a single number\n", "- Two sample test\n", " - Compare two sets of data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A one-sample bootstrap hypothesis test\n", "Another juvenile frog was studied, Frog C, and you want to see if Frog B and Frog C have similar impact forces. Unfortunately, you do not have Frog C's impact forces available, but you know they have a mean of 0.55 N. Because you don't have the original data, you cannot do a permutation test, and you cannot assess the hypothesis that the forces from Frog B and Frog C come from the same distribution. You will therefore test another, less restrictive hypothesis: The mean strike force of Frog B is equal to that of Frog C.\n", "\n", "To set up the bootstrap hypothesis test, you will take the mean as our test statistic. Remember, your goal is to calculate the probability of getting a mean impact force less than or equal to what was observed for Frog B if the hypothesis that the true mean of Frog B's impact forces is equal to that of Frog C is true. You first translate all of the data of Frog B such that the mean is 0.55 N. This involves adding the mean force of Frog C and subtracting the mean force of Frog B from each measurement of Frog B. This leaves other properties of Frog B's distribution, such as the variance, unchanged." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "def bootstrap_replicate_1d(data, func):\n", " \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n", " bs_sample = np.random.choice(data, len(data))\n", " return func(bs_sample)\n", "\n", "def draw_bs_reps(data, func, size=1):\n", " \"\"\"Draw bootstrap replicates.\"\"\"\n", " \n", " # Initialize array of replicates: bs_replicates\n", " bs_replicates = np.empty(size)\n", " \n", " # Generate replicates\n", " for i in range(size):\n", " bs_replicates[i] = bootstrap_replicate_1d(data, func)\n", " \n", " return bs_replicates" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 0.0053\n" ] } ], "source": [ "# Make an array of translated impact forces: translated_force_b\n", "translated_force_b = force_b - np.mean(force_b) + 550\n", "\n", "# Take bootstrap replicates of Frog B`s translated impact forces: bs_replicates\n", "bs_replicates = draw_bs_reps(translated_force_b, np.mean, 10000)\n", "\n", "# Compute fraction of replicates that are less than the observed Frog B force: p\n", "p = np.sum(bs_replicates <= np.mean(force_b)) / 10000\n", "\n", "# Print the p-value\n", "print('p = ', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A two-sample bootstrap hypothesis test for difference of means\n", "We now want to test the hypothesis that Frog A and Frog B have the same mean impact force, but not necessarily the same distribution, which is also impossible with a permutation test.\n", "\n", "To do the two-sample bootstrap test, we shift both arrays to have the same mean, since we are simulating the hypothesis that their means are, in fact, equal. We then draw bootstrap samples out of the shifted arrays and compute the difference in means. This constitutes a bootstrap replicate, and we generate many of them. The p-value is the fraction of replicates with a difference in means greater than or equal to what was observed." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "forces_concat = np.concatenate((force_a, force_b))" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value = 0.0036\n" ] } ], "source": [ "# Compute mean of all forces: mean_force\n", "mean_force = np.mean(forces_concat)\n", "\n", "# Generate shifted arrays\n", "force_a_shifted = force_a - np.mean(force_a) + mean_force\n", "force_b_shifted = force_b - np.mean(force_b) + mean_force\n", "\n", "# Compute 10,000 bootstrap replicates from shifted arrays\n", "bs_replicates_a = draw_bs_reps(force_a_shifted, np.mean, 10000)\n", "bs_replicates_b = draw_bs_reps(force_b_shifted, np.mean, 10000)\n", "\n", "# Get replicates of difference of means: bs_replicates\n", "bs_replicates = bs_replicates_a - bs_replicates_b\n", "\n", "# Compute and print p-value: p\n", "p = np.sum(bs_replicates >= empirical_diff_means) / 10000\n", "print('p-value =', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A/B testing\n", "- Used by organizations to see if a strategy change gives a better result\n", "- Null hypothesis of an A/B test\n", " - The test statistics is impervious to the change" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The vote for the Civil Rights Act in 1964\n", "The Civil Rights Act of 1964 was one of the most important pieces of legislation ever passed in the USA. Excluding \"present\" and \"abstain\" votes, 153 House Democrats and 136 Republicans voted yea. However, 91 Democrats and 35 Republicans voted nay. Did party affiliation make a difference in the vote?\n", "\n", "To answer this question, you will evaluate the hypothesis that the party of a House member has no bearing on his or her vote. You will use the fraction of Democrats voting in favor as your test statistic and evaluate the probability of observing a fraction of Democrats voting in favor at least as small as the observed fraction of 153/244. (That's right, at least as small as. In 1964, it was the Democrats who were less progressive on civil rights issues.) To do this, permute the party labels of the House voters and then arbitrarily divide them into \"Democrats\" and \"Republicans\" and compute the fraction of Democrats voting yea." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value = 0.0002\n" ] } ], "source": [ "# Construct arrays of data: dems, reps\n", "dems = np.array([True] * 153 + [False] * 91)\n", "reps = np.array([True] * 136 + [False] * 35)\n", "\n", "def frac_yea_dems(dems, reps):\n", " \"\"\"Compute fraction of Democrat yea votes.\"\"\"\n", " frac = np.sum(dems) / len(dems)\n", " return frac\n", "\n", "# Acquire permutation samples: perm_replicates\n", "perm_replicates = draw_perm_reps(dems, reps, frac_yea_dems, 10000)\n", "\n", "# Compute and print p-value: p\n", "p = np.sum(perm_replicates <= 153/244) / len(perm_replicates)\n", "print('p-value =', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A time-on-website analog\n", "It turns out that you already did a hypothesis test analogous to an A/B test where you are interested in how much time is spent on the website before and after an ad campaign. The frog tongue force (a continuous quantity like time on the website) is an analog. \"Before\" = Frog A and \"after\" = Frog B. Let's practice this again with something that actually is a before/after scenario.\n", "\n", "We return to the no-hitter data set. In 1920, Major League Baseball implemented important rule changes that ended the so-called dead ball era. Importantly, the pitcher was no longer allowed to spit on or scuff the ball, an activity that greatly favors pitchers. In this problem you will perform an A/B test to determine if these rule changes resulted in a slower rate of no-hitters (i.e., longer average time between no-hitters) using the difference in mean inter-no-hitter time as your test statistic. The inter-no-hitter times for the respective eras are stored in the arrays ```nht_dead``` and ```nht_live```, where \"nht\" is meant to stand for \"no-hitter time.\"" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "nht_dead = np.array([ -1, 894, 10, 130, 1, 934, 29, 6, 485, 254, 372,\n", " 81, 191, 355, 180, 286, 47, 269, 361, 173, 246, 492,\n", " 462, 1319, 58, 297, 31, 2970, 640, 237, 434, 570, 77,\n", " 271, 563, 3365, 89, 0, 379, 221, 479, 367, 628, 843,\n", " 1613, 1101, 215, 684, 814, 278, 324, 161, 219, 545, 715,\n", " 966, 624, 29, 450, 107, 20, 91, 1325, 124, 1468, 104,\n", " 1309, 429, 62, 1878, 1104, 123, 251, 93, 188, 983, 166,\n", " 96, 702, 23, 524, 26, 299, 59, 39, 12, 2, 308,\n", " 1114, 813, 887])\n", "\n", "nht_live = np.array([ 645, 2088, 42, 2090, 11, 886, 1665, 1084, 2900, 2432, 750,\n", " 4021, 1070, 1765, 1322, 26, 548, 1525, 77, 2181, 2752, 127,\n", " 2147, 211, 41, 1575, 151, 479, 697, 557, 2267, 542, 392,\n", " 73, 603, 233, 255, 528, 397, 1529, 1023, 1194, 462, 583,\n", " 37, 943, 996, 480, 1497, 717, 224, 219, 1531, 498, 44,\n", " 288, 267, 600, 52, 269, 1086, 386, 176, 2199, 216, 54,\n", " 675, 1243, 463, 650, 171, 327, 110, 774, 509, 8, 197,\n", " 136, 12, 1124, 64, 380, 811, 232, 192, 731, 715, 226,\n", " 605, 539, 1491, 323, 240, 179, 702, 156, 82, 1397, 354,\n", " 778, 603, 1001, 385, 986, 203, 149, 576, 445, 180, 1403,\n", " 252, 675, 1351, 2983, 1568, 45, 899, 3260, 1025, 31, 100,\n", " 2055, 4043, 79, 238, 3931, 2351, 595, 110, 215, 0, 563,\n", " 206, 660, 242, 577, 179, 157, 192, 192, 1848, 792, 1693,\n", " 55, 388, 225, 1134, 1172, 1555, 31, 1582, 1044, 378, 1687,\n", " 2915, 280, 765, 2819, 511, 1521, 745, 2491, 580, 2072, 6450,\n", " 578, 745, 1075, 1103, 1549, 1520, 138, 1202, 296, 277, 351,\n", " 391, 950, 459, 62, 1056, 1128, 139, 420, 87, 71, 814,\n", " 603, 1349, 162, 1027, 783, 326, 101, 876, 381, 905, 156,\n", " 419, 239, 119, 129, 467])" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-val = 0.0002\n" ] } ], "source": [ "# Compute the observed difference in mean inter-no-hitter times: nht_diff_obs\n", "nht_diff_obs = diff_of_means(nht_dead, nht_live)\n", "\n", "# Acquire 10,000 permutation replicates of difference in mean no-hitter time: perm_replicates\n", "perm_replicates = draw_perm_reps(nht_dead, nht_live, diff_of_means, 10000)\n", "\n", "# Compute and print the p-value: p\n", "p = np.sum(perm_replicates <= nht_diff_obs) / 10000\n", "print('p-val =', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test of correlation\n", "- Hypothesis test of correlation\n", " - Posit null hypothesis: the two variables are completely uncorrelated\n", " - Simulate data assumning null hypothesis is true\n", " - Use Pearson correlation, $\\rho$, as test statistics\n", " - Compute p-value as fraction of replicates that have $\\rho$ as least as large as observed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hypothesis test on Pearson correlation\n", "The observed correlation between female illiteracy and fertility may just be by chance; the fertility of a given country may actually be totally independent of its illiteracy. You will test this hypothesis. To do so, permute the illiteracy values but leave the fertility values fixed. This simulates the hypothesis that they are totally independent of each other. For each permutation, compute the Pearson correlation coefficient and assess how many of your permutation replicates have a Pearson correlation coefficient greater than the observed one." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "df = pd.read_csv('./dataset/female_literacy_fertility.csv')\n", "fertility = np.array(df['fertility'])\n", "illiteracy = np.array(100 - df['female literacy'])" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "def pearson_r(x, y):\n", " \"\"\"Compute Pearson correlation coefficient between two arrays\n", " \n", " Args:\n", " x: arrays\n", " y: arrays\n", " \n", " returns:\n", " r: int\n", " \"\"\"\n", " # Compute correlation matrix: corr_mat\n", " corr_mat = np.corrcoef(x, y)\n", " \n", " # Return entry[0, 1]\n", " return corr_mat[0, 1]" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-val = 0.0\n" ] } ], "source": [ "# Compute observed correlation: r_obs\n", "r_obs = pearson_r(illiteracy, fertility)\n", "\n", "# Initialize permutation replicates: perm_replicates\n", "perm_replicates = np.empty(10000)\n", "\n", "# Draw replicates\n", "for i in range(10000):\n", " # Permute illiteracy measurements: illiteracy_permuted\n", " illiteracy_permuted = np.random.permutation(illiteracy)\n", " \n", " # Compute Pearson correlation\n", " perm_replicates[i] = pearson_r(illiteracy_permuted, fertility)\n", " \n", "# Compute p-value: p\n", "p = np.sum(perm_replicates >= r_obs) / 10000\n", "print('p-val =', p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Do neonicotinoid insecticides have unintended consequences?\n", "As a final exercise in hypothesis testing before we put everything together in our case study in the next chapter, you will investigate the effects of neonicotinoid insecticides on bee reproduction. These insecticides are very widely used in the United States to combat aphids and other pests that damage plants.\n", "\n", "In a recent study, Straub, et al. ([Proc. Roy. Soc. B, 2016](http://dx.doi.org/10.1098/rspb.2016.0506)) investigated the effects of neonicotinoids on the sperm of pollinating bees. In this and the next exercise, you will study how the pesticide treatment affected the count of live sperm per half milliliter of semen." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "control = np.array([ 4.159234, 4.408002, 0.172812, 3.498278, 3.104912, 5.164174,\n", " 6.615262, 4.633066, 0.170408, 2.65 , 0.0875 , 1.997148,\n", " 6.92668 , 4.574932, 3.896466, 5.209814, 3.70625 , 0. ,\n", " 4.62545 , 3.01444 , 0.732652, 0.4 , 6.518382, 5.225 ,\n", " 6.218742, 6.840358, 1.211308, 0.368252, 3.59937 , 4.212158,\n", " 6.052364, 2.115532, 6.60413 , 5.26074 , 6.05695 , 6.481172,\n", " 3.171522, 3.057228, 0.218808, 5.215112, 4.465168, 2.28909 ,\n", " 3.732572, 2.17087 , 1.834326, 6.074862, 5.841978, 8.524892,\n", " 4.698492, 2.965624, 2.324206, 3.409412, 4.830726, 0.1 ,\n", " 0. , 4.101432, 3.478162, 1.009688, 4.999296, 4.32196 ,\n", " 0.299592, 3.606032, 7.54026 , 4.284024, 0.057494, 6.036668,\n", " 2.924084, 4.150144, 1.256926, 4.666502, 4.806594, 2.52478 ,\n", " 2.027654, 2.52283 , 4.735598, 2.033236, 0. , 6.177294,\n", " 2.601834, 3.544408, 3.6045 , 5.520346, 4.80698 , 3.002478,\n", " 3.559816, 7.075844, 10. , 0.139772, 6.17171 , 3.201232,\n", " 8.459546, 0.17857 , 7.088276, 5.496662, 5.415086, 1.932282,\n", " 3.02838 , 7.47996 , 1.86259 , 7.838498, 2.242718, 3.292958,\n", " 6.363644, 4.386898, 8.47533 , 4.156304, 1.463956, 4.533628,\n", " 5.573922, 1.29454 , 7.547504, 3.92466 , 5.820258, 4.118522,\n", " 4.125 , 2.286698, 0.591882, 1.273124, 0. , 0. ,\n", " 0. , 12.22502 , 7.601604, 5.56798 , 1.679914, 8.77096 ,\n", " 5.823942, 0.258374, 0. , 5.899236, 5.486354, 2.053148,\n", " 3.25541 , 2.72564 , 3.364066, 2.43427 , 5.282548, 3.963666,\n", " 0.24851 , 0.347916, 4.046862, 5.461436, 4.066104, 0. ,\n", " 0.065 ])\n", "\n", "treated = np.array([1.342686, 1.058476, 3.793784, 0.40428 , 4.528388, 2.142966,\n", " 3.937742, 0.1375 , 6.919164, 0. , 3.597812, 5.196538,\n", " 2.78955 , 2.3229 , 1.090636, 5.323916, 1.021618, 0.931836,\n", " 2.78 , 0.412202, 1.180934, 2.8674 , 0. , 0.064354,\n", " 3.008348, 0.876634, 0. , 4.971712, 7.280658, 4.79732 ,\n", " 2.084956, 3.251514, 1.9405 , 1.566192, 0.58894 , 5.219658,\n", " 0.977976, 3.124584, 1.297564, 1.433328, 4.24337 , 0.880964,\n", " 2.376566, 3.763658, 1.918426, 3.74 , 3.841726, 4.69964 ,\n", " 4.386876, 0. , 1.127432, 1.845452, 0.690314, 4.185602,\n", " 2.284732, 7.237594, 2.185148, 2.799124, 3.43218 , 0.63354 ,\n", " 1.142496, 0.586 , 2.372858, 1.80032 , 3.329306, 4.028804,\n", " 3.474156, 7.508752, 2.032824, 1.336556, 1.906496, 1.396046,\n", " 2.488104, 4.759114, 1.07853 , 3.19927 , 3.814252, 4.275962,\n", " 2.817056, 0.552198, 3.27194 , 5.11525 , 2.064628, 0. ,\n", " 3.34101 , 6.177322, 0. , 3.66415 , 2.352582, 1.531696])" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'ECDF')" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhTZdr48e85J92p3SilpSICUioUqIILgii7UoTxHWWoy7iA844I6rwuyKiIgCPqiA6yKCOMiDq+/GRACgN1HREE8RVpEQErSLEtha7SQmlzcn5/hISkTTeaJml6f65rruEkJ8nzUMzdZ7tvxTAMAyGEEMIF1dsNEEII4bskSAghhKiXBAkhhBD1kiAhhBCiXhIkhBBC1EuChBBCiHp5JEgsWLCA4cOHk5SUxMGDB13eo+s6c+bMYeTIkYwaNYo1a9Z4omlCCCEaYPLEh4wYMYI777yT2267rd57NmzYQG5uLpmZmZSVlTFx4kSuvvpqEhMTm/VZpaWVWCx1j37ExHSguLii2W33Vf7WH5A+tQX+1h+QPqmqQlRUWL3PeyRIDBw4sNF7Nm3axC233IKqqkRHRzNy5Eg2b97MlClTmvVZFovhMkjYnvMn/tYfkD55il6Ygzl/P6aE3mhxPZv1Wl/sT0v5e59y8so5kFtKUtcoenaJaNb7eCRINEVBQQEJCQn26/j4eI4dO+bFFgnhn/TCHE5lvACWGqrVAELTHmt2oPAXOXnlfJ5VQGJMaLO/PNuKnLxyXnxvN2bdgklTeXRyarP66jNBwl1iYjrU+1xsbLgHW9L6/K0/IH1yt6pfDnD6yPeEXNSH4MQkAEoPHuaUxQyGARYzQeWHieqb2uT39Jef0f6fS3jpn7sxmy2YTCrz//saeneL9naz3Mb2c/o8qwBdt2AYoOsWfik+xdUDmj6N7zNBIj4+nvz8fPr16wfUHVk0VXFxhcuhY2xsOCdOnGxxO32Fv/UHpE/u5jhiKHUYMegRF4NqAosZVBNnIi5uchvd2Z+WTIG4w46sPGrM1i9Ps9nCjqw8YsICPN6O1uD4c0qMCUXTVNAtaJpKYkyo089QVZUGf7n2mSAxduxY1qxZw+jRoykrK+Pjjz/mnXfe8XazhGizzPn7wVJjHzGY8/ejxfVEi+tJaNpj570m0VyugkFLp0DcIalrFCZNRT/75ZnUNcqjn+8pPbtE8OjkVN9ek5g3bx6ZmZkUFRVx9913ExkZycaNG5k6dSozZswgJSWFCRMmsGfPHkaPHg3AtGnTuPDCCz3RPCH8khLcAesudwuoJkwJve3P2YJFa6svGBzILcXsMAVyILfU40HC9uX5S/Epv16TAGtfz7d/HgkSTz75JE8++WSdx5cvX27/s6ZpzJkzxxPNEcLv6YU5nNn+Lhg6KBpBg9O9sji9PbuAGrPF2iaHYOArv8X37BLB1QMS/W6a0518ZrpJCNE8DW1jtU81AWBgVHn+HMDn3+XxxZ58+7WqKvZg0NIpEOE5EiSEaIMa28ZqSuhNtRpgX5x2nGryhJy8clZnHsRxD8mQlHinYNCSKRDhORIkhPCQlhxgq62+RWkbdyxOt2T30YHcUqddhpqqMDglvtltEN4nQUIID3D3AbamjBSaujjdGruPkrpGEWBSMZstKKrCbaN7yaihjZIgIYQHNPabf3O5axtra+0+kjUH/yFBQoha3DktZNMaawQt3caak1fO+i8PYTZbMHD/7iNZc/APEiSEcNBaeY08fYCtMTl55bzw7reYdeu6gaLgFAxkJCBsJEgIwbnRg1FR7NZpIUeeOsDWFNuzC+wBAqBb53Amj+wlu49EHRIkRLvnOHpA0az/c3FK2V/k5JVzpND58NhFceESEIRLEiREu1dzcBvo1WevLJh6D0PtEOMT00LN0ZQtq7aFatspaABNk+2pon4SJES7phfmUHNg67kHFJXAXte0qeAAzmsMJk3hsfTLXAYK264lAAW49OIoJgzpLqMIUS+P1LgWwleZ8/db8xsBoGBKGtrmAgQ4rzGYdYPt2QUu77PtWlIVMJlUCRCiUTKSEO2WXpiDpaLYaQ0isNc13m5WkzlOLzWV7FoSzSVBQrRLtRerA3oPI6CNTDPl5JWzPbuAL7ML0C0GJk1l8shL0DQFi26gNrLGILuWRHNIkBDthuMhOacT0FhQOsS0mQBRe+FZ1y1Unq7h8fTLZIQg3E6ChGgXah+SCxqcDl7MktoUrnYrOS4829gOwckIQbQGCRLCb9U7crCYMaoqfOoEdG2ff5fHO5kHsVgMTKZzOZUc02UoqsLQlHgG10rBLYQ7SZAQfqkpIwdfOgHtyF6LwWLbrXQup5IsPAtPkyAh/FJbGznY2JLuOdZiUBXFaQeTTCsJT5IgIfySEtwB6zEgi8+PHGxsU0y6Y4CQWgzCyyRICL+jF+ZwZvu71kNyikbQ4HSfDg5Qd4oJoI+chhY+QIKE8Ct6YQ5nvlkHes3ZRwyMqgqvtqkhth1MxeVVdcp9SoAQvkCChPAb9sVqe7I+xee2t9Y+JW2rCqeqCiZNQdcNKfcpfIoECdGmudzmCoCC1qUPQQMn+sxUU+1Sodf07WwvEWpYDIb2TyAmIlh2LQmfIkFCtFmNbXP1pQAB1iR8tpPS+tkDcY4lQuW8g/BFEiREm9WWtrnm5JWzNSvffq2q1vxKg1Pi5cyD8GkSJESb4ji9ZEroTbWPH5Db/3MJO7LyOJT/K47ZNFK6x9iDggQH4cskSIg2o/b0UmjaYz47cgDr6OGlfzon47OJCAv0QouEaD4JEsLn2UYPRkWx0/SSOX8/QalpPhccbA7kllJjti5MA6iKtemNpfIWwpdIkBA+rXbdB8cCQb60tdWVsJAAe4AAGHNFV0KDTbL+INoUCRLCp9Uc3OZw7sGCqfcw1A4xPjm9VFvusZNO16fPmLnlet9usxC1SZAQPksvzKHmwNZzDygqgW2kehxAeWV14zcJ4eMkSAifZc7fb82/BICCKWmozwcI24nqsJAAsg4V2x/XVGQdQrRJHgsShw8fZubMmZSVlREZGcmCBQvo1q2b0z3FxcU88cQTFBQUUFNTw1VXXcWTTz6JySSxrD1yyuSqBRDY6xpvN6lBjieqFUXBOJuLSQGG9kuQdQjRJqme+qDZs2eTnp7Oli1bSE9P5+mnn65zz7Jly+jRowcbNmxgw4YNfP/992RmZnqqicKHOGdyVX0+k6utDoTZfC7NhqoqqAqYTKqMIkSb5ZFf0YuLi9m3bx8rV64EIC0tjblz51JSUkJ0dLT9PkVRqKysxGKxUF1dTU1NDXFxcZ5oovABVb8c4My+b13kYfL9TK4vvPstZv3syEGxptuYPPISUFUSY0JlFCHaLI8EiYKCAuLi4tA0DQBN0+jUqRMFBQVOQeL+++9n+vTpDBkyhNOnT3Pbbbdx+eWXe6KJwsv0whwKNr6AoddfbtRXbc8usAcIgG6dw5k80prFNTY2nBMnTjbwaiF8m09N9m/evJmkpCTeeustKisrmTp1Kps3b2bs2LFNfo+YmA71PhcbG+6OZvoMf+hP1S8HOH3ke/j1BIZuth+UC9VqiLn9GU4f+Z6Qi/oQnJjk7abWKyQkwOm6d7dorh6QaL/2h5+TI3/rD0ifGuKRIBEfH09hYSG6rqNpGrquc/z4ceLjnedpV69ezXPPPYeqqoSHhzN8+HB27tzZrCBRXFzhVLzFxt9+o/OH/tQ5KKdqYAFUE2ciLsYclAC9EjgJnPSBvjrWgnCcPoqNCHa6LzYi2P6z8YefkyN/6w9In1RVafCXa48sXMfExJCcnExGRgYAGRkZJCcnO001ASQmJvLFF18AUF1dzVdffcUll1ziiSYKL3DK4mpYCO93PYEDbyY07TGfW6S27Vxa+8UhXnxvNzl55fbnKk/X2P+s1LoWoq3z2O6mZ555htWrVzNmzBhWr17NnDlzAJg6dSrZ2dkAzJo1i//7v/9j/PjxTJw4kW7dunHrrbd6qonCg/TCHCwVxWdTbaigmgjvd53P5mI6kFtq37lk1i0cyC21P5fUNYoAk2rfyWSrOieEP/DYmkSPHj1Ys2ZNnceXL19u/3PXrl3tO6CE/6o9zRTQexgBva4hODHJJ6aVXDlVZcY2iWkY1rxMNj27RPDo5FSpCyH8kk8tXIv2oXY+JqVDjE+OHmw+/y6PzV/nOj1We0qpZ5cICQ7CL3lsukkIcJ2PyZe3t+bklbM686BTNldNVWRKSbQbEiSER7W1fEwHckuddsspCtw2upeMGkS7IUFCeNS5fExKm8jHdKrK7HQ99oquXDegi5daI4TnSZAQHtMW8zFt2XXU6bHQYFnGE+2LBAnhMeb8/aCfXfA1LD6djwms6TYcp5pkLUK0R/JrkfAY61STfSPp2Wvf4lgPYmt2gf1xVdYiRDslQUJ4jKXoiMOV4nMjiYbqQVzbP0HWIkS7JEFCeESdra+q5nNbX7dnF1BjtlgvDGs9CMMw0DSpByHaLwkSwiN8fevr59/l8cWefPu1pimkj+pF5ekaOUUt2jUJEqLVOeVpwgKqyae2vtoOzDkmDx6SEi/TS0IgQUK0svryNPnSKKL2gTlNVWR6SYizJEiIVqMX5nDmm3Vnt70a+GqeprBaRYNGD7pQppeEOEuChGgV9hGEPZGf4rNlSLN/Kna6Pn3GXM+dQrQ/EiSEW+mFOZjz92NUFFunmABQ0Lr0IWjgRJ8bReTklfNdTpG3myGEz5IgIdymTjlSh4VqXwkQtUuQHsgtdcrwqirIeoQQDiRICLdxKkeKBVPvYagdYjAl9PZagHAMCoD9sJxJU3l0cqq9qpzZbEFRFW6XU9VCOJEgIdzmXIbXc9tcvTl6cDxBbdJUrunbGbNuLUGqny1BOu7qblJVTogGSJAQbuGc4VXziQyvjieodd36/yZNRdctaNq5WtRSVU6I+kmQEC3ieqHa8Hpeppy8crZmnTtBrZ49+zA4JV5GDUI0gwQJcd4aWqj29lbX7dkFnB08AJDSPcYeFCQ4CNF0EiTEeTtXH8IAQ8eUfJ3XF6rBOoo4UnjS6bGIsEAvtUaItk2ChDhvtetDaB0vIjD5Oi+2yJqo753Mg+iOaTY0SbMhxPmSICHOm/O6g/frQ9gT9TkEiD4XRzFhSHeZYhLiPEn5UnHelOAOZ9chFNACfGIdonaiPgkQQrSMBAlxXpy3vKpe3/JaZzeTlBsVwi0kSIjzUnNwm0PyPu9vea29m6l/z45SD0IIN5AgIZqtTilSRfX6VFN5ZbXTtexmEsI9JEiIZqs5uA0stnTa3i9F+vl3eexxyOSqqZKkTwh3kd1NolnqjCJUzaulSD//Lo+3txxwyuQ6tF+CrEUI4SYykhDNYs7fb12sBrw9irBteXUMEFJ6VAj3kpGEaDK9MAdLRbFT+g1vjSJy8spZ/+Uhpy2viuxoEsLtJEiIJqmdpymg9zACvJQK3NWpavVsLQjZ0SSEe0mQEA2qk+X1bEEhpUOMVwKEnKoWwrM8FiQOHz7MzJkzKSsrIzIykgULFtCtW7c6923atImlS5diGAaKorBy5Uo6duzoqWYKB76Y5VVOVQvhWR4LErNnzyY9PZ0JEyawfv16nn76aVatWuV0T3Z2Nq+99hpvvfUWsbGxnDx5ksBA2e/uSbaRgymht0+VI83JK2d7doGcqhbCwzwSJIqLi9m3bx8rV64EIC0tjblz51JSUkJ0dLT9vn/84x/cc889xMbGAhAeHu6J5omzHEcO1WoAQYPTQQ2wnonwYjlSWxlSW5U5AAW4tn+CrEEI0co8EiQKCgqIi4tD0zQANE2jU6dOFBQUOAWJn376icTERG677TZOnTrFqFGj+OMf/4iiKE3+rJiYDvU+FxvrX0HH3f0pPXiYUxazdeRgMROq1RBz+zOcPvI9IRf1ITgxya2f54qrPq35z6E6ASIgQGXc0B5t4mfaFtrYHP7WH5A+NcSnFq51XefAgQOsXLmS6upqpkyZQkJCAhMnTmzyexQXVzjNWdvExoZz4sRJF69om1qjP9V6ANavYAVUE2ciLsYclAC9EjgJnGzlv7/afXI1xaSp1sNyg1PiiQkL8Pmfqfy7833tvU+qqjT4y7VHDtPFx8dTWFiIrlsPYem6zvHjx4mPdz70lJCQwNixYwkMDKRDhw6MGDGCrKwsTzSx3fPFrK4vvrebz7/LtyfuU7AGiDvH9pZ1CCE8xCNBIiYmhuTkZDIyMgDIyMggOTnZaaoJrGsVX375JYZhUFNTw44dO+jd27uJ49oLX8zq6jjFBGAyqXKaWggP81hajmeeeYbVq1czZswYVq9ezZw5cwCYOnUq2dnZAIwbN46YmBhuvPFGJk6cSM+ePfntb3/rqSa2W76W1bV2bQhNhesGJPDo5FQZQQjhYR5bk+jRowdr1qyp8/jy5cvtf1ZVlSeeeIInnnjCU80SnM3HZNHt197O6lq7NkS/Hh25c6yMKIXwBknwJ6xlSHE4oNbxIu81xgWpDSGE90iQELXWHxSvr0eEBDkPcLt29r/tiUK0FY0GiT179niiHcKLTAm9QT37xaxqXl+P2PJ1rtNjladrvNQaIUSjQeLuu+92ur755ptbrTFCbM8uwPGYi6pAUtco7zVIiHau0SBhGM4H044ePdpqjRGtRy/M4czuDPTCnDrPORUSMizWax/Rv2dH2dEkhBc1urupdkqM5qTIEL6hdk6m0LTHnHYvmRJ6U+2Qo8mb00211yNSesR4qSVCCGhCkDCbzXzwwQf2EUVNTQ3/7//9P6d75CyDb3JZC8Jixpy/3ylIaHE9CU17zJ791VvbX/f/XMKWXc4jVVmPEMK7Gg0S/fv3Z926dfbrvn37sn79evu1oigSJHxQc2tBaHE9vXo2AuDTb47WqRUh6xFCeFejQeLtt9/2RDuEm/lSLYjG2BP5ZRfYH5NaEUL4hiafuC4vLycrK4vy8nIiIyNJSUkhIkL+A/ZVpoTeVCsaGGZQVK/VgmiM1IoQwrc1KUgsXryY119/HV3XiYqKoqSkBJPJxH333ccDDzzQ2m0UfkwS+Qnh2xoNEps2bWL16tW8+OKLjBgxApPJhNls5uOPP2bu3Ll0796dG2+80RNtFc1Qc3CbdbcS2Le1+tpI4vPv8vhij0MiP01haEo8g1PiZZpJCB/RaJBYs2YNM2fOZMyYMedeZDIxduxYqquref/99yVI+Bhfy+rqSk5eOaszDzodnBt9RVduGdbDe40SQtTR6GG6H374gWHDhrl8btiwYezf7zsHr4SVr2V1tcnJK2fjVz+Tk1fOgdzSOjuZhg/s6r3GCSFcanQkUV1dTWRkpMvnIiIiqKmRfey+xhezun7+XR7vZB7EYjEwmVQmj7yEAJOK2WxBURVuG92L3t2i/a6MpBBtXaNBwjCMBlNx1E7bIbzPV7K62kYMYSEB1qmlsyMHs26h8nQNj05O5UBuKUldo2QNQggf1WiQOH36NKNHj643GEiaDt+jBHewHp4zLKAFeGU9wra11axbUBTFaWpJVRR7YJDgIIRvazRIyJpD26IX5nBm+7vWhH2KRtDgdK+sRxzILcVstmBgHW1qqjVQ2KaWJDgI0TY0GiSqqqrIzc2lV69edZ47ePAgF110EUFBQa3SONF0dfI0AWB4baopLCQAx7Hn6EEXEhpskqklIdqYRoPE3//+d3799VdmzZpV57m1a9cSHh7OtGnTWqVxommam6fJE3KPOS9Anz5j5pbrvb/DSgjRPI1ugd20aRP33nuvy+fuvvtuNm7c6PZGieapObgN9GprnibDgilpKIEDb66TEtyTyiurvfK5Qgj3anQkUVhYSFxcnMvn4uLiKCwsdHujRNO5Ojjn7TxNOXnlZP1UZL/WVCTNhhBtVKMjiZCQEAoKClw+l5+fT0hIiNsbJZrOqaocik8cnNueXYDukI6pXw+pLidEW9VokBg2bBgvv/yyy+deffXVek9ji9anF+ZgqSg+uw6hghZAYK9rvN2sOlNNEWGBXmqJEKKlGp1ueuihh5g0aRI33XQTo0ePJjY2lhMnTvDRRx9RUVHBP//5T0+0U9RS9csBp8XqgN7DCPCBdOCff5fHnhyZahLCXzQaJGJjY/nXv/7FihUr2Lp1K2VlZURGRnL99ddz9913S00JLzl95HvQawADDB2lQ4zXA4SrpH1D+yXIVJMQbVijQWLevHk8+eSTPPzww4A1K+wtt9xif3769OksWrSo9VooXNJCwjmXn8k4m6/Je3Lyyln/5aE6SftkFCFE29bomsTatWudrl988UWn623btrm3RaJJzhQecrjyXn4msE4xLXjnW74/XGp/TJWT1UL4hSYl+GvoWnieXpjDqT2fnntA1bx2aM4+xeQwguhzcRQThnSXACGEH2h0JFE7gZ8k9PM+a70I2x5T7217rW+KSQKEEP6j0ZGEruvs2LHDPoIwm81O1xaLpaGXCzezb3tVNbAAqskr215t9SF0x+yuMsUkhN9pNEjExMQ45W2KjIx0uo6Ojm6dlok6nHI0qd7b9vr5d3m8veUAjjOPMsUkhH9qNEh8+umnjd0iPMSeownAYvHKtlfbGoRjgJApJiH8V6NrEu5y+PBhJk2axJgxY5g0aRI///xzvfceOnSI/v37s2DBAk81z+fVydHkpcXq7dkFTmsQioJMMQnhxzwWJGbPnk16ejpbtmwhPT2dp59+2uV9uq4ze/ZsRo4c6ammtQnWxWrdfh3e73qPjiJy8spZtXk/W7Py7Y+pCtwxJonrBnTxWDuEEJ7lkSBRXFzMvn37SEtLAyAtLY19+/ZRUlJS59433niD6667jm7dunmiaW2G9bDcud/ggzp399hn20qRfv5dvj1xnwJc2z9BAoQQfq7RNQl3KCgoIC4uDk3TANA0jU6dOlFQUOC08L1//36+/PJLVq1axZIlS87rs2Ji6j95HBsbfl7v6QtKD9ZwBgVroFDQT5/0WH/W/OcQNeZzu9gUICBAZdzQHm5vQ1v+GdXH3/rkb/0B6VNDPBIkmqKmpoannnqKv/zlL/Zgcj6Kiyuc5sxtYmPDOXHipItX+D69MIfqwnzrtlfDWnEu5KI+HunP59/lsWXHz/ZrTbXmYxqcEk9MWIBb29CWf0b18bc++Vt/QPqkqkqDv1x7JEjEx8dTWFiIrutomoau6xw/fpz4+HN5fU6cOEFubi733XcfAL/++iuGYVBRUcHcuXM90UyfVLs0qW3ba3BiEidb+R92fQn77hzrndPdQgjP80iQiImJITk5mYyMDCZMmEBGRgbJyclOU00JCQns3LnTfr1o0SJOnTrF448/7okm+iynba94dtvrgdxSSdgnRDvnsd1NzzzzDKtXr2bMmDGsXr2aOXPmADB16lSys7M91Yw2xVVpUk9uew0LCXC6Hj3oQtnqKkQ747E1iR49erBmzZo6jy9fvtzl/dOnT2/tJvkUvTAHc/5+TAm97SMFb5cmzT3mPJ11+ozZY58thPANPrNw3Z45rjtUqwGEpj2GFtcTU0JvqtUAsJi9kqOpdhlSIUT7I0HCB1gPytWAYYDFjDl/P1pcT7S4noSmPVZnhOEJUoZUCAESJHyC9aCcCli3tzquO9iChSdJGVIhhI3HFq6Fa3phDme2v2tde1BUgganezzdxsavfiYnr9z+mOxqEkLYyEjCS2wL1UZFsXWqCQDDY2VIc/LK2Z5dwJfZBegWA5Om8ujkVHp2iSCpaxQBJhWz2YIiNSKEaNckSHhB7QNyKBqupppai6uCQbpu4UBuKT27RNCzSwSPTk7lQG4pSV2jJEAI0Y5JkPCC2gfkTL2HoXaI8cjitKua1ACappLUNcp+bQsWQoj2TYKEh7k6IBfooepyrmpSq4o1m+vglHgJCkKIOiRIeJi3Dsjl5JXzwrvfYtada1LfPrqXpPsWQtRLgoSHOW131QI8dkBue3aBU4C4OD6cySNlQVoI0TDZAutB3trumpNXzpFC5xQbF8WFS4AQQjRKRhIeZM7fD/rZ7a6GpVW3u+bklXMgt5SwkADe+/hHp6JBmibnHoQQTSNBwoOcS5AaZ6/dz1Zu1KxbUBQF4+xCtQJcenEUE4Z0l1GEEKJJJEh4kKXoiMOV0mojie3ZBedGDoaBqioYhoGmqRIghBDNIkHCQ+psfVW1Vjk4l5NXztasfPu1pimkj+pF5ekaORgnhGg2CRIe4qmtr9uzC9DPLT+Q0j1GtrgKIc6bBIlWZsvRpAR3AA/UhqhdAyIiLLBVPkcI0T5IkGhFTjma1ACCBqdjVFW0SvoNW8K+rJ+kBoQQwn0kSLQipxxNFjNGVQVBqWlu/QzHbK6Oh+UUpAaEEKLlJEi0Elc5mty9UG3b6up4BsLGZFJlFCGEaDEJEq2k5uA26/oD0FoL1U5bXc/SNIWhKfGSsE8I4RYSJFqBq+2u7l6o3vzVz3yxx2Grq2qdXpLgIIRwJwkSraC1t7vm5JWzbG1WnRrUd45t/YJFQoj2RRL8tYJzmV6VVsn0eiC31KmqnNSgFkK0FgkSbuaJTK+nqsxO16MHXShTTEKIViFBws1aO9NrTl45W3YddXosNFhmDYUQrUOChJu1ZqZXV+VHNVVxqk0thBDuJL+CullrZXr9/Ls83sk86LQWoaoKt42W6nJCiNYjQcKNWivTa05eOaszDzqNIFJ7xXLDlV0lQAghWpUECTdqra2v27ML6kwxpY/pTUxYQIvfWwghGiJBwo3ObX21tGjrq630qG2tYWt2gf05VYHbRveid7doTpw4Wd9bCCGEW0iQcBPnra/aeW99dSw9atJUrunb2T6KUIBr+ydIfQghhMfI7iY3ccr4itGkBeucvHI2fvUzOXnl9scO5JZi1i0YBuhnqweZNBVVkaR9QgjP89hI4vDhw8ycOZOysjIiIyNZsGAB3bp1c7pn8eLFbNq0CU3TMJlMPPzwwwwdOtRTTTxvTc34WnsayXHE8OjkVHp2iSCpaxQmTUXXLWiaNSgMTom3v04WqoUQnuSxIDF79mzS09OZMGEC69ev5+mnn2bVqlVO9/Tr14977rmHkJAQ9u/fz+23386XX35JcHCwp5p5XpqS8dW2hYo3ZiAAABuzSURBVNViMTCZrNNIjiOGA7ml9OwSQc8uETw6ObVOUJDgIITwBo9MNxUXF7Nv3z7S0qwFd9LS0ti3bx8lJSVO9w0dOpSQkBAAkpKSMAyDsrIyTzTxvDUl46ttC6tuMTAAc61pJE1TnQ7E9ewSwbiru0lgEEJ4nUdGEgUFBcTFxaFpGgCaptGpUycKCgqIjo52+Zp169bRtWtXOnfu3KzPiomp/4RzbGx4s96rKUoPHuaUYavpoBDefzixfVOd7vk8y3kLq6oojBvag3FDe5D9UxEpPTrSu5vrv4eGtEZ/vE365Pv8rT8gfWqIT+5u+vrrr3n11VdZsWJFs19bXFzh9IVsExsb3ipbRqv1AKz7jqwZX/WuVzh9Tk5eOdk/nnB6zehBF9rPOFzXz7oQ3dy2tVZ/vEn65Pv8rT8gfVJVpcFfrj0SJOLj4yksLETXdTRNQ9d1jh8/Tnx83Z06u3fv5tFHH2XJkiV0797dE807b41te3VVXlRBEvIJIdoOj6xJxMTEkJycTEZGBgAZGRkkJyfXmWrKysri4Ycf5m9/+xt9+vTxRNNapHbG16LjRU5bWm3bWR2ZTKok5BNCtBke+5X2mWeeYebMmSxZsoQLLriABQsWADB16lRmzJhBSkoKc+bMoaqqiqefftr+uhdeeIGkpCRPNbNRjttYuzpkfDUw2PxdMduqDtm3tDpuZ1VUqT0thGh7PBYkevTowZo1a+o8vnz5cvufP/jgA08157wc2fsd//fpFxysjuPDbXHMHlSEbWnIAEKoctrSOu7qbi63swohRFshk+ONsI0c+oaXEfHVa4wJMjMySGPJydHk1CSQqgWCxYyiaBy2xNfZ0mo7+yCEaDldN1NaegKzubrxm5vo+HEVi8XS+I1tiKs+mUyBREXFomnN+9qXIOHAcSqpZ5cIpzxKFSF7uSFYR1EAw0KvwEISkscSmvIY5vz9mBJ6c6s5VkYNQrSi0tITBAeHEhbWGUVR3PKeJpOK2exfQaJ2nwzDoLLyV0pLT9CxY/NS+0iQOKt2Yj3bNJFZt3CRdoILlAoMRQXDAFXlquHXclGXCCDCvqOpJ3IyWojWZDZXuzVAtBeKohAWdgEVFc0/nCxB4qzt2QX2raq2NYWkrlH0CCziD6GZaOgoaAQmDyOg1zVEuqFOhBCi+SRAnJ/z/XuTLLBYRxFbs/Lt1+rZutE9u0Rwb3IZAYqOpoCCgdIhxi2FhIQQwtHJkyd55523zvv1mzZt4MknH3Nji6wkSGAdRTgeZ0jpHkPPLhHohTmE/vI19vhbT3ZXIYRoqYqKk7z77qp6nzebzfU+15ra/XRTTl45Rwqdj69HhAUCZw/LWXT74+4qRyqE8Kzam1Lcae/eLBYvfpVTp04BMG3ag4SHh/PKKy9RVXWa4OAQHnroEZKT+1BQkM+UKXdw0003s2PHNqqqqpg582n69x/Ayy8voKKigrvuSic4OJhly1bwwAP3kZLSn3379hIYGMiLL77Kv/+dwXvvvY2iKCQkJPLYY7OIimp+7rematdBwpa+W3esH60p9sI+isNhOQCt40WebqIQooVcbUpxV6D49ddyZs16lPnzXyAlpT+6rlNeXsaUKXfyxBNPM2jQlXzzzdf8+c+P8f776wAoLy+nb99+/OEP08jM/DfLlv2NpUtX8Kc/Pc6UKXfwj3+86/QZhw7l8Ne/LsJkMnHoUA7Llr3Gm2+upmPHjixfvpSFC1/k2Wf/4pb+uNJup5ty8spZveWAU4Doc3EUj6dfZv8H5FxdTmlStTkhhG+pXe3xQG6p2957795sunW7mJSU/oA1w3VpaSkBAQEMGnQlAAMHXkFAQAC5uUcACAkJ5ZprrMXU+vRJIS8vr8HPGDVqLCaT9ff5b7/9hquvvoaOHTsCMGHCzXzzzddu648r7TZIbM8uwDFZrKrAhCHdnX7DUII7gKJhy/Aq6xFCtD229Diuare0lGHUzThtGIbLnUS2hwIDA+yPqaqKrje81hASEurw3nV3KbX2Zq92GyTKK51PbPbv2dEpQDhneFXrZHgVQrQNtmqPv7m2u1unmgBSUvrx88+H2bs3CwBd14mOjqa6uppvv/0GsP72bzabufDChqerw8LCqKqqanCB+vLLB/HVV9soLi4CYMOGdQwceIWbeuNau1yT+Py7PPbkFNmvNRVuuOoi9MIc++lp66L12QyvGDLVJEQb1lrpcS64IIL5819g0aKFVFWdRlFUpk17kPnzX3BauJ43bwEBAQGNvtfo0Tfw+9//jvDwC1i2rG49ne7de/CHP0zj4YennV247sKjj85ye78cKYar8VIb1ljRoZy8cp5/51une64bkMBtqSZOZbxgDQxqAEGD060jCYsZVBOhaY/51EiivRdKaSv8rU/e7s+xY0fo3Nm9G0jaQ1oOG1d/fz5RdMiXHMgtdQoQmmrdzWTO32oNEIYBFjNGVQWhaefyMvlSgBBCCE9pd0EiLMR5yDd60IXWg3Om3lQrGhhm+6E5La6nBAchRLvW7hauc485D5VPn/HOKUYhhGgL2l2QCCz7mZHB2XQznXB63Jy/37qTCcCwWK+FEKKda1fTTbu2fsWo8jVoITo6GksrRjM45XIATAm9qVYD7AvVciZCCCHaUZA4sGsXRtYGTCYd9WzhoJGJp+zb4rS4nrJQLYQQtbSLIKEX5qB9/Fd6mWpQAIsBOirxfVKd7pOFaiGEcNYu1iTM+ftRLGY0BSzAgZp4vr3wDi7qO8DbTRNCtGFvvvk6NTU1jd/YDD/+eIBPPvnovF5bUJDPuHEj3NqedhEkiqpULIZtBKGxuWoA1VHdvN0sIUQbt3LlcpdBoiW1H3788SCffXZ+QaI1+P10k16YQ9jetYCBBYW1pwZx1NKJdDcm+RJC+DbHlDvumlL+618XAPDHP96DoqjEx8fTqVMcR48epayslBUrVvP993tZtmwRlZWVAEyZ8t8MHjwEs9nMY489RHl5OWfOnOHSS/vw6KOzOHWqkr//fRmnTlVy113pDBiQykMPPVrv+wB88MH/8r//+y4xMR1JTb3cLX1z5PdBwpy/H9UwoyhgMQzClDP2A3RCCP+nF+bYU+5UqwFuS7HzP//zOP/61xqWLl1BaGgo8+c/w9692bz22huEhIRw8uRJXnrpOV588W907NiRoqIipk69k1Wr3qdDhw7Mnj2PiIhIDMNg3rzZbNy4nokTf8uUKf/N9u1bmTfvBYAG36ew8BirVq1g5cp3iI6O4aWXnm9xv2rz+yBRVKUSblhLBylApRFEQrDfd1sIcZY9WefZlDvm/P2ttkHluutGEBISAsDevXsoKMjnkUdm2J9XFIW8vKNcckkS7723mh07tmOx6Jw8eZLg4GCX79nQ+2RnZzF48BCio2MAmDDhN26fqvL7b8ujuYUkY60XYTGgg3rGrfnkhRC+zZNnoEJDQ+x/Ngzo0eMSFi9eXue+zZs3kpX1HUuWLCc0NIxVq1Zw9Giuy/ds6H2ysva4r/H18PuF62o1BIWzxTqAqJgYmWoSoh2xnYEKHHiz27M5h4aGUVnpuoxA3779+OWXXHtdCYAffvgewzCoqDhJREQkoaFhVFRU8NFHm+33hIVZH2vK+1x22UC++mobpaUlAGRkrHdb32z8fiTRKcyCpQQ0BXQDukX7fVwUQtTSWmegfve725gx478JCgomPj7e6bkLLriA559/mcWLX+XVV/+K2VxDQkIXFixYyNixaWzd+gW3334rsbGx9O+fypkzZwC4/PIreO+91fz+95NJTb2Mhx56tN736dnzEu64427++Md7iY6O4eqrh7i9j35fTyJj/ScMPvYOGhZ0VLZ3vo20Ce7dR+wN3s7r3xqkT77P2/2RehJN4856EvJrtRBCiHr5/XRTz4BCNCz2fE09Awq93SQhhGgz/H4kcbxSdVq4Pl7p910WQgi38ftvzEDLaesZCcV6ViLQctrbTRJCtICfLaN6zPn+vXksSBw+fJhJkyYxZswYJk2axM8//1znHl3XmTNnDiNHjmTUqFGsWbOmxZ8bFhnpNJIIi4xs8XsKIbzDZAqksvJXCRTNZBgGlZW/YjIFNvu1HluTmD17Nunp6UyYMIH169fz9NNPs2rVKqd7NmzYQG5uLpmZmZSVlTFx4kSuvvpqEhMTz/tzK8vKsHBuC2xlWVkLeyKE8JaoqFhKS09QUeG+/45VVcVi8a/dTa76ZDIFEhUV2+z38kiQKC4uZt++faxcuRKAtLQ05s6dS0lJCdHR0fb7Nm3axC233IKqqkRHRzNy5Eg2b97MlClTzvuzi4K7oqOBYd0CWxTctcX9EUJ4h6aZ6NgxvvEbm8Hb23pbgzv75JEgUVBQQFxcHJqmAaBpGp06daKgoMApSBQUFJCQkGC/jo+P59ixYy367N4DB7L0/SJ6qMf4ydKZ3w0c2KL3E0KI9sTvtsDWPhQSGxtOVOTNZP9UxOgeHendLbqeV7Y9sbHh3m6C20mffJ+/9QekTw3xSJCIj4+nsLAQXdfRNA1d1zl+/HidY+zx8fHk5+fTr18/oO7IoilKSyudTlwDxIYHcsuIXhQXV1Bc7DrPSlsTE9PBb/piI33yff7WH5A+qapCVFRYvc97JEjExMSQnJxMRkYGEyZMICMjg+TkZKepJoCxY8eyZs0aRo8eTVlZGR9//DHvvPNOsz6roc42dPS8LfK3/oD0qS3wt/6A9KkhHsvd9NNPPzFz5kx+/fVXLrjgAhYsWED37t2ZOnUqM2bMICUlBV3XefbZZ9m2bRsAU6dOZdKkSZ5onhBCCBf8LsGfEEII9/H7E9dCCCHOnwQJIYQQ9ZIgIYQQol4SJIQQQtRLgoQQQoh6SZAQQghRLwkSQggh6uX3QaIpdSzaktLSUqZOncqYMWMYP348DzzwACUlJd5ullu89tprJCUlcfDgQW83pcXOnDnD7NmzGT16NOPHj+epp57ydpNa5LPPPmPixIlMmDCB8ePHk5mZ6e0mNduCBQsYPnx4nX9jbfk7wlWf3P4dYfi5O+64w1i3bp1hGIaxbt0644477vByi1qmtLTU2LFjh/36+eefN5544gkvtsg99u7da9x7773GddddZxw4cMDbzWmxuXPnGvPnzzcsFothGIZx4sQJL7fo/FksFmPgwIH2n8sPP/xgDBgwwNB13csta55du3YZ+fn5xvXXX+/0b6wtf0e46pO7vyP8eiRhq2ORlpYGWOtY7Nu3r03/5h0ZGcmVV15pvx4wYAD5+flebFHLVVdX8+yzzzJ79mwURfF2c1qssrKSdevW8eCDD9r707FjRy+3qmVUVeXkSWt9gpMnT9KpUydUtW19fQwcOLBOUtG2/h3hqk/u/o7wu1Thjppax6KtslgsvPfeewwfPtzbTWmRV199lZtuuokLL7zQ201xi6NHjxIZGclrr73Gzp07CQsL48EHH2RgG61loigKr7zyCvfffz+hoaFUVlby+uuve7tZbiHfEY1rW78KCCdz584lNDSU22+/3dtNOW+7d+8mOzub9PR0bzfFbcxmM0ePHuXSSy9l7dq1PPLII0yfPp2KiraZjtpsNvP666+zZMkSPvvsM5YuXcrDDz9MZWWlt5smGuGO7wi/DhKOdSyAeutYtEULFizgyJEjvPLKK21u2O9o165dHDp0iBEjRjB8+HCOHTvGvffey5dffuntpp23hIQETCaTfQqjf//+REVFcfjwYS+37Pz88MMPHD9+nMsvvxyAyy+/nJCQEH766Scvt6zl5DuicW3326UJHOtYAPXWsWhrFi5cyN69e1m8eDGBgYHebk6L3HfffXz55Zd8+umnfPrpp3Tu3Jk333yTIUOGeLtp5y06Oporr7zSnvL+8OHDFBcXc9FFF3m5Zeenc+fOHDt2jEOHDgHWtP9FRUV07dr268XLd0Tj/D5VeH11LNqqH3/8kbS0NLp160ZwcDAAiYmJLF682Mstc4/hw4ezbNkyevXq5e2mtMjRo0eZNWsWZWVlmEwmHnroIYYNG+btZp23Dz/8kOXLl9sX4mfMmMHIkSO93KrmmTdvHpmZmRQVFREVFUVkZCQbN25s098Rrvr0yiuvuPU7wu+DhBBCiPPn19NNQgghWkaChBBCiHpJkBBCCFEvCRJCCCHqJUFCCCFEvSRIiCbLz88nNTXVfvDojjvuYM2aNQCsXbuWyZMn2+9NTU3l6NGjXmmnK++++y6DBw8mNTWV0tLSZr125syZLFy4EIBvvvmGMWPGtEYThfBJEiREkyUkJLB79257npuG7N6922dyMdXU1PD888+zYsUKdu/eTVRU1Hm/18CBA9myZYsbWydaYvjw4Wzfvt3bzfBrEiSE3ysuLubMmTP07NnT203xKrPZ7O0mOPG19gjXJEi0c8OHD+fvf/8748ePZ8CAAcyaNYuioiKmTJlCamoqd911F+Xl5QD88ssvJCUlNek/7qSkJI4cOQJYU0s/9thjXHXVVVx//fUsWbIEi8UCnJumWrBgAYMGDWL48OH85z//sb/P2rVrGTFiBKmpqQwfPpwPP/zQ5edVV1czf/58hgwZwpAhQ5g/fz7V1dUcPnyYsWPHAjBo0CDuvPNOl6+fMWMG11xzDZdffjm33XYbP/74o8v7du7cybXXXgvAG2+8wYwZM5yenzdvHvPmzbP3e9asWQwZMoShQ4eycOFC+1RdbVlZWdx8881cdtllDB48mL/85S/Aub/z999/3963FStW2F9nsVh44403GDlyJFdeeSUPPvggZWVlTq9ds2YN1113Hb///e/tj33wwQcMGzaMQYMG8d5775GVlcX48eMZOHAgzz77rMs2AixatIgZM2bw0EMPkZqaym9+8xv2799vf76wsJDp06dz1VVXMXz4cFatWlXntY888giXXXYZ//rXv+q8/8yZM3nmmWfs//5+97vfceLECebPn8+gQYMYO3Ys+/btq7d9wv0kSAgyMzNZuXIlW7Zs4bPPPmPq1Kn86U9/YufOnVgsFt5+++0Wvf/cuXM5efIkH3/8MW+//Tbr16/ngw8+sD+flZXFxRdfzI4dO5gyZQp//vOfMQyDU6dOMW/ePJYvX87u3bv55z//SXJyssvPWLp0KXv27GH9+vV8+OGHZGdns2TJEi6++GJ7Xp5du3Y5fWk5uvbaa9myZQtfffUVl156KY888kij/Ro3bhz/+c9/7NlddV1n8+bN9sR+jz/+OCaTiczMTNatW8e2bdvsazi1zZ8/nzvvvJNvv/2Wjz76iBtuuMHp+Z07d5KZmcmbb77JG2+8YZ9iWbVqFR9//DGrV69m69atRERE1PmS37VrF5s2beLNN9+0P7Znzx4yMzNZuHAhzz33HMuWLeMf//gHGzdu5N///jdff/11vf3+5JNPGDt2LF9//TVpaWncf//91NTUYLFY+OMf/0hSUhJffPEFb731Fm+99RZbt26t89pvvvmG8ePHu3z/f//73zz00EPs2LGDwMBAJk2aRJ8+fdixYwdjxoyxB1DhGRIkBLfffjsdO3YkLi6OgQMH0q9fPy699FICAwMZNWpUi35z03WdTZs28T//8z906NCBxMRE7r77bqcRQUJCArfeeiuapvGb3/yGEydOUFRUBFiL3fz4449UVVXRqVMnLrnkEpefs2HDBqZNm0ZMTAzR0dFMmzat3lGHK7/97W/p0KEDgYGBTJ8+nf3799uL7NSnS5cuXHrppXz88ccA7Nixg+DgYAYMGEBRURFffPEFs2bNIjQ0lJiYGO666y42btzo8r1MJhO5ubmUlJQQFhbGgAEDnJ6fNm0aoaGhJCUlcfPNN9sD3/vvv8/DDz9M586dCQwM5IEHHmDLli1Oo73p06cTGhpqz+Nje7+goCCGDBlCaGgoaWlpxMTE2P8NNPQz79OnD2PHjiUgIIC7776b6upq9uzZQ3Z2NiUlJTzwwAMEBgZy4YUXcuutt7Jp0yb7awcMGMDIkSNRVdWpPY5GjRpF3759CQoKYtSoUQQFBTFx4kQ0TePGG2/khx9+aPDnItzLr4sOiaZxrJoWFBTkdB0cHMypU6fO+71LS0upqakhISHB/lhCQgKFhYUuPz8kJASAU6dOERsby8KFC1mxYgV//vOfueyyy3j88cfp0aNHnc85fvx4nc84fvx4k9qo6zoLFy5k8+bNlJSU2NMql5aWEh4e3uBr09LSyMjIYOLEiWRkZNhHEfn5+ZjNZqdsthaLpd4U1PPnz+dvf/sbN9xwA4mJiTzwwANcf/319ucdX9elSxd7PeP8/HymTZvmlApaVVWKi4vt1507d67zeTExMfY/BwUF1blu6Gfu+H6qqhIXF2f/uz5+/LhTcSVd152uXbWlobYFBwe79d+jaD4JEqJVRUVFERAQQH5+vn3h2FYNrCmGDh3K0KFDqaqq4pVXXuGpp57i3XffrXNfp06dyM/Pt480CgoK6NSpU5M+Y8OGDXzyySesXLmSxMRETp48yaBBg2hK7ssbbriBBQsWcOzYMT766CPef/99APtv9jt27MBkavw/s27duvHyyy9jsVjIzMxkxowZ7Ny50/58QUGBPTjm5+fb+9a5c2eee+45e60HR7/88guA20vCHjt2zP5ni8VCYWEhnTp1QtM0EhMTyczMrPe1/lCetr2R6SbRqjRNY+zYsSxcuJCKigry8vJYuXIlN910U6OvLSoq4pNPPuHUqVMEBgYSGhpa7/bbcePGsXTpUkpKSigpKWHx4sX1znnXVllZSWBgIFFRUZw+fZqXX365yf2Ljo7miiuu4IknniAxMdH+Rd6pUyeuueYann/+eSoqKrBYLOTm5tY7179+/Xr7KOaCCy4AcOrrkiVLOH36ND/++CNr167lxhtvBGDy5Mm88sor5OXlAVBSUmKf/mot33//PZmZmZjNZt566y0CAwPp378//fr1o0OHDrzxxhtUVVWh6zoHDx4kKyurVdtTU1PDmTNn7P+TXVPuJUFCtLqnnnqKkJAQRo4cSXp6OmlpafzXf/1Xo6+zWCysXLmSoUOHcsUVV7Br1y5mz57t8t7777+fvn37ctNNN3HTTTfRp08f7r///ia1b+LEiSQkJDB06FDGjRtXZz2gMWlpaWzfvt0+1WTzwgsvUFNTw4033sigQYOYMWMGJ06ccPkeW7duZdy4caSmpjJ//nwWLlxIUFCQ/fkrrriCUaNGcdddd3HPPffYp7HuvPNOhg8fzj333ENqaiq33nprq38pjxgxgk2bNjFo0CDWr1/PokWLCAgIQNM0li5dyv79+xkxYgRXXXUVTz75ZKuXbb3vvvvo16+f/X+LFi1q1c9rb6SehBA+7JdffmHEiBF8//33TZq2am2LFi3iyJEjvPTSS95uivAQGUkIIYSolwQJIYQQ9ZLpJiGEEPWSkYQQQoh6SZAQQghRLwkSQggh6iVBQgghRL0kSAghhKiXBAkhhBD1+v+lT1Ui0UVYRQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute x, y values for ECDFs\n", "x_control, y_control = ecdf(control)\n", "x_treated, y_treated = ecdf(treated)\n", "\n", "# Plot the ECDFs\n", "_ = plt.plot(x_control, y_control, marker='.', linestyle='none')\n", "_ = plt.plot(x_treated, y_treated, marker='.', linestyle='none')\n", "\n", "# Set the margins\n", "plt.margins(0.02)\n", "\n", "# Add a legend\n", "plt.legend(('control', 'treated'), loc='lower right')\n", "\n", "# Label axes\n", "plt.xlabel('millions of alive sperm per mL')\n", "plt.ylabel('ECDF')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bootstrap hypothesis test on bee sperm counts\n", "Now, you will test the following hypothesis: On average, male bees treated with neonicotinoid insecticide have the same number of active sperm per milliliter of semen than do untreated male bees. You will use the difference of means as your test statistic." ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value = 0.0001\n" ] } ], "source": [ "# Compute the difference in mean sperm count: diff_means\n", "diff_means = np.mean(control) - np.mean(treated)\n", "\n", "# Compute mean of pooled data: mean_count\n", "mean_count = np.mean(np.concatenate((control, treated)))\n", "\n", "# Generate shifted data sets\n", "control_shifted = control - np.mean(control) + mean_count\n", "treated_shifted = treated - np.mean(treated) + mean_count\n", "\n", "# Generate bootstrap replicates\n", "bs_reps_control = draw_bs_reps(control_shifted, np.mean, size=10000)\n", "bs_reps_treated = draw_bs_reps(treated_shifted, np.mean, size=10000)\n", "\n", "# Get replicates of difference of means: bs_replicates\n", "bs_replicates = bs_reps_control - bs_reps_treated\n", "\n", "# Compute and print p-value: p\n", "p = np.sum(bs_replicates >= np.mean(control) - np.mean(treated)) / len(bs_replicates)\n", "print('p-value =', p)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }