Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
ME3263_Lab-0/ME3263_lab-00.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
561 lines (561 sloc)
85.1 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "from scipy.stats import norm, t\n", | |
| "\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "from ipywidgets import interact\n", | |
| "import ipywidgets as widgets\n", | |
| "\n", | |
| "import check_lab00 as p" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import pretty_plots # script to set up LaTex and increase line-width and font size" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# ME 3263\n", | |
| "## Laboratory # 0\n", | |
| "## Introduction to the Student t-test" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We use statistics to draw conclusions from limited data. No measurement is exact. Every measurement you make has two types of uncertainties, *systematic* and *random*. *Systematic* uncertainties come from faults in your assumptions or equipment. \n", | |
| "*Random* uncertainties are associated with unpredictable (or unforeseen at the\n", | |
| "time) experimental conditions. These can also be due to simplifications of your\n", | |
| "model. Here are some examples for caliper measurements:\n", | |
| "\n", | |
| "In theory, all uncertainies could be accounted for by factoring in all physics\n", | |
| "in your readings. In reality, there is a diminishing return on investment\n", | |
| "for this practice. So we use some statistical insights to draw conclusions. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Mean and standard deviation\n", | |
| "\n", | |
| "If the same measurement is taken multiple times, then there should be some average value, the mean, $\\mu$. If there is an average value, then that means there is also a measure of deviation from the average value, standard deviation, $\\sigma$. The definitions for mean and standard deviation are as such\n", | |
| "\n", | |
| "$\\mu = \\sum_{i=1}^{N}\\frac{x_{i}}{N}$\n", | |
| "\n", | |
| "and\n", | |
| "\n", | |
| "$\\sigma^2 = \\frac{\\sum_{i=1}^{N}(x_{i}-\\mu)^2}{N}$,\n", | |
| "\n", | |
| "where $x_i$ is the $i^{th}$ measurement in a dataset called $x$ and $N$ is the number of data points. \n", | |
| "\n", | |
| "If you know the mean and standard deviation of a normally distributed data set, then you can predict the probability a given measurement will occur. The [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) for the difference between measurement $x_{i}$ and $\\mu$ is shown below for $\\sigma=1$. \n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, 'Probability of\\\\\\\\ Measurement (\\\\%)')" | |
| ] | |
| }, | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEbCAYAAACLGcAmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU5b0/8M83YYeQBRKQNSRhExFIIrKIKCRYaxeXBFzqVhWwtr3WWpb23t+9vddbCnaxVqUEW4taFRKwva2tmrAoFRSSsFhBlkwCYQ9kYQkh2/f3x/PMzMkwk5lJZubM8n2/XvNKzpkz5zwhYb7znOf7fB9iZgghhBDBKsrsBgghhBDtkUAlhBAiqEmgEkIIEdQkUAkhhAhqEqiEEEIENQlUQgghgloXsxsQjvr378/JyclmN0MIIUJKSUnJWWZOdNwfkYGKiBYBmMfMGf44f3JyMoqLi/1xaiGECFtEdMTZ/oi79UdEWQCWA0hv55h0Isonohr9yCcil8cLIYTwn4gKVEQUByDfzTFZAEoA5ACo1o8cACVElOP3RgohhGgjogIVVJCKc/WkDmSFejObmVOZORVAtvX1+hghhBABEjGBSo9LZQFYAMDi4rD5+usKZi6y7tTf5zkcI4QQIgAiIlDp8aXlAIqYOa+dQ+fpr2udPLfK4RghhBABEBGBCuqWXy2AXDfHpQAAM5c6PmHYJ0kVQggRQGGfnk5Eq6ACUDYz17o5XMafRNg4fOYC8j62YO+xOjhbzad71yjMGpOER6YlI65Xt8A3UAgPhXWg0ll68+Ew5uRGe8GsFi6CGRHN19fCsGHDvGmmED71+bE6vLz5MD7Yd8ppgDLae6wOqz+24FtThuOxGSOQFNMjMI0UwgsUrgsn6uy8cgDVOnPP+FwZgBRmJof9DACO+z193iozM5Nlwq8ItB3l1Xhp82F8fLCqQ6/v1iUK8zKHYsHMFAyJ7+Xj1gnhHhGVMHOm4/5w7lHNher9VBNRocNzKQBg2L9c97gs1ufa4e72oRABdfp8A55+Zze2W85d9dzsMUn49k0j0K/P1bf29p88j1c2l+HQmYsAgMbmVrzx6RG8veMoHpw6HD/56lh0iY6UYWwRzMI5UFmlwHXwydJfrRl9FgApRJTumFBhqEzhKrVdiIA7eq4eD/z+U1RWX7btIwLuGH8NvnNLGq4d1Nfla8cM7ItvThiMwv2n8fLmw9h7rA4A0NzKeO2TChyvuYwX75uEHl2j/f5zCNGesP24xMx5zEzOHtDBxrCvQL/MWrXCWQq6dd8qJ88JEXAHTl1Azu+22YJUdBRhbuYQbHxmJl66P73dIGUVFUW4bdxA/OWp6XjjscmYPCLB9tyH+07j23/ciUtXmv32MwjhibAdo2qPqzEq/Zz1HyTD2qvSvakSwP34FCBjVML/dlfW4pHXdqC2vgkA0L1LFF55IB2zxw7o1HmZGcv+8SXyPrbfOJg4NA5/fPQGyQwUfudqjCpse1SdYJ1rVUJEhXocq8ThOSFMs63sLB5Y/aktSPXp3gVrvj2500EKAIgIS28fgx/dNtq2b3dlLeat+hRnzjd0+vxCdIQEKgf6NmAGgAIAmfpRANXDKmjvtUL4W+G+03jktZ241NgCAIjv1RVvPXEjpqT089k1iAhP3ZqG//nmONu+A6cvIHfVdlRW1/vsOkJ4KiJv/fmb3PoT/lC47zQWvlmCllb1f3ZA3+5487EbMXJAjN+u+e6uY3g2f6/tmgP79sCG70zDoLiefrumiFxy60+IEFZZXY9n1u22BYxhCb1QsHCaX4MUANw1aQhWPpCObjpN/dT5Bnz/7V1obmn163WFMJJAJUSQa2ppxffe3oULDSr7bnBcTxQsnIqhCYGZlDtn3ED8/pFMREepPKLiIzV4oehQQK4tBCCBSoig94sPDmB3pZpn3iWK8NL9k5DUN7CljmaMTMQPskbatl/echj/PHQ2oG0QkUsClRBBbPOBM1hlSBX/0W2jMWlYvCltefKWNExPU0kbzMDTa3ej6sIVU9oiIosEKiGC1Km6Bvxw3R7b9i2jE/HEDHcVvvwnOorw63kT0V+XYzp78QqeWbcbra2SkCX8SwKVEEGopZXx9NpdqL7UCEBl+P0ydwKiotzON/erpJge+PW8iSDdjK2HzmLlR2WmtkmEPwlUQgSh3246hE8t1QCAKAJemDcJ/fp0N7lVyoyRiXhypn1Bgl8VHkRxRbWJLRLhLuiK0hLRLKgisqlQ1c/LoCqWFzPzbjPbJkQgfGo5hxc32rPqvjdrJKam+m5Cry88kz0Kn5VXo+RIDVpaGd9/exf+/m8zpMyS8Iug6FER0eNEtJOIWgAUAsgDsBjAAgAr9HYJEZ0joneI6FYTmyuE3zQ0teDZ/D2wDvvcOCIB3589sv0XmaBLdBRevG8SYnt2BQCcqGvAc+/tN7lVIlyZFqiIqC8RLdPBKQ9APID1ABZC1dTL0I9svf08VM29uQCKiOgQEd1lSuOF8JPVH1twrEZVQ4/r1RW/uXeSbf5SsBkc1xMrcq63bReUHLOl0QvhS6bc+iOiuwG8qjefB7CKmcvdvGy94fVZUL2t9UT0IYC5zHzeL40VIkBO1F7Gy1sO27afnTMaA2ODe2n428YNRNbYASjafxoA8F//9wU2PDnN9KQPEV7M6lG9CmAxMycw8xIPglQbzFzEzLkAEgBUAFjqhzYKEVDL/vElGppUaaKx1/TFfZOHmdwiz/zH18baSiztrqzFhl3HTW6RCDemBCodoFb74Dy1zLyQmSVQiZD2meUc/rrnhG37p98YF7S3/BwN79cbj88YYdte/v6XuNDQZGKLRLgJimQKVwzjWB8Q0c+IyP2SpUKEmJZWxn/9dZ9t++sTBrVZaTcUPHVrGgb0VenzVReu4KVNh928QgjPBW2gIqJYAOVQ2X/9ACwBYJFgJcLN2zuOYv9JNcTao2sUlt4+xuQWea939y5YevtY2/YfPimHpeqiiS0S4SRoAxWA1VDzqFL1+iRpUGNSeaa2Sggfqq1vxC8/PGDbfuqWtJBd6+mbEwchY7iqQ9jUwvifv+1z8wohPBPMgSoLQJ410YKZLVDBK9vUVgnhQ78uPIgavaT8kPieeOJm82r5dRYR4affGGcrr7T5QBU2fXna3EaJsGBKoCKiWUSU7MGhjmWipfqlCBtfnjqPNz87atv+9zuuRY+u0Sa2qPOuGxyLe28Yatv+n7/tx5XmFhNbJMKBWT2qVKjxplfaGXPKB5BLRBMAgIgmAZgPYF2A2iiE3zAz/vuv+2wr9t6U1h+3jRtgcqt849k5oxHTQ03RLD97Ca99UmFug0TIMys9fTWAOfpRQ0Q/c3LYIqg5UqVEdAhAMQALVHKFECHtk8PnsK3sHAC1fMZ/fv1aEIVGOro7/fp0xw+yRtm2V24pw3lJVxedYNoYlZ60mwbgSQALdR2/xwzP1zFzKtRk3l0AljBzmlSgEKGOmfFC0UHb9r03DMXIATEmtsj3Hpw6HMP79QIA1F1uwhrpVYlOMD2ZgpnzmDkBwHIAq3UNv1sNz69g5rnM/Lx5rRTCdz45fA7FR2oAAF2jCU/dmmZyi3yva3QUvjfLXkz31X+WS69KdJjpgcqKmVdApZ9vArCRiN4nouEmN0sIn3LsTc27YWjIpqO7c+fEQdKrEj4RNIEKsJVEWgA1Zyoa7hMuhAgpjr2p79wSfr0pqy7SqxI+Ymqg0mnqh4iohYgOWm/5MbOFmbMB3AZgMlwnXAgRMiKpN2UlvSrhC2auRzUbapHEOqjySBeg1pkyjk8V6aoU8+Ak4UKIUBJJvSkr6VUJXzCzR7UcgIWZM5n5eWbOgEpHX+54IDMXGBIupISSCDmR2Juykl6V6CwzA1U6gCKHfYVQq/o6ZUi4ECKkRGJvykp6VaKzzAxUFgCOhc1S9H6XmLnOby0Swg8iuTdlJb0q0RlmBqo8ANl6namJRLQMqhDtKhPbJITPRXJvykp6VaIzzKxMsQKqGvoSACVQpZHymPkXZrVJCF+T3pSd9KpER5manq7nTKUCmAu17tRCM9sjhK9tt0hvyspZr+rSlWYTWyRChVnLfNgm8DJzOTOvt6471dnzCRFMVn9sH3LNzYzc3pTVnRMHYViCvVeVX1xpcotEKDCrR1Wrx6Y6FWD0hOGdUIVrhQgqh05fwOYDVQAAIuCJGaG7KKKvdImOwuMzRti2f/9JuW2pEyFcMStQpcG+xMc7xkm+7hBRMhE9q5f+KARQxMwSqETQeXWr/SbBnGsHYET/3ia2JnjkZAxBXK+uAIDK6sv44ItTJrdIBLsuZlxULyufSUTzodadmktEDKAUat2pMgC1+vA4AP2gUtez9DYBKAAwpzO3DIXwlzMXGvDuruO27fkhvMS8r/Xq1gUPThmO3246DABY9bEFt183MGzW4xK+Z0qgsmLmPAB5RJQFIBfAbAAL2nlJEVQvKk/mU4lg9vq2I2hsaQUATBoWh4zhMk/d6KGpyVj1kQWNLa3YU1mL4iM1uCFZ/o2Ec6YGKitmLoKuUkFEsVC9pxSoKhQWANXMvMu8FgrhufrGZrzx6RHb9nwZm7pKYkx33DVpMNbqZIq8jy0SqIRL7Y5REdFKInrcYXuWPxukV/bdpTMBVzPzRglSIpQUlBxD3WU1mXVYQi/MGTfQ5BYFJ2NSRdH+07BUXTSxNSKYuUumWAAgx2F7kv+aI0Roa2nlNkkUj900AtFRMvbizMgBMZg1JgkAwAz8/p8y3Cycc3frbxdUmaOdAKr1voVENMeDczMzf6VTrRMixHz4xSkcra4HAMT27IrczCEmtyi4PTEjBZu+PANA9USfyR6Ffn26m9wqEWzcBarHobLrrBXNGaqSRKoH55bJESKiMDNWGSb4PjhlOHp1C4ph4KA1JSUB4wfH4vPjdbjS3Io3Pj2Cp7NGmd0sEWTavfWnx4pSmTkKKrGBoNaEivfgISOjIqKUHKnB7ko1q6JbdBQemjbc5BYFPyJqM1b1+vYjaGhqMbFFIhh5POGXmWuhKp4X6oQHtw//NVuI4LN6q703deekQUiK6WFia0LHV8dfg8G6tFT1pUZsKD3u5hUi0nhVmYKZFzLzJsf9UmtPRLrys5fw4b7Ttu3HJSXdY12jo/Do9GTb9qtbLWiVskrCoEMllPT6UR8Q0TkiaoEqhXSOiN73phySEOHij5+Ug/V7662jEzFqQIy5DQox904ehpgeajzPcvYSNh84Y3KLRDDxOlAR0c+h1o/KhhqLKtePeKj6fUVE9IqX51zmyfwsIrqbiJ71ts1C+NOFhiYUlByzbT92k/SmvNWnexfcN3mYbfuP2yrMa4wIOl4FKiK6B6o2Xx2AXGaOYuY0/YiCWlfqPIAFRHSXF6deDM/mZ90LlcwhRNAoKDmGS40qAWBkUh9MT+tncotC04NThsNa7m/robMokwnAQvM2d3YBVNp5OjNXOD7JzAVEVArgsD72XWcn0WWS8hx2LySiye1cOw6qKG1tO8cIEVCtrYw1hk//D09LluKqHTQ0oReyxg5AoR7re31bBX76zetMbpUIBt4GqkwApc6ClBUzW3SwuqGd8yRAFaG1vQyez89yDHBCmOajQ1WoOKcm+Mb06IK7Jg02uUWh7ZFpybZAVVByDM/eNhoxPbqa3CphNrOW+SgnImtQIqge2CoAK9y8tFrS3kUwMfam5mUORe/uMsG3M6al9sPIpD44dOYiLjW2oKDkGB6dPsL9C0VY8zaZogRAOhG5nMmob+ulQ60r5ZJegr5cr01VBDU/q9zNQ4KUCBqWqovYYljB96GpyeY2KAwQER6elmzbXrOtQlLVhdeBahVUD6iQiCY4PklEE6ECFAPI9/SkzDyHmTd42RYhTPX6dvtSHrNGJ2FYv14mtiZ83DVpsC1VveJcPT46VGVyi4TZvLpPoZMl1gO4B0ApEdXC3nPKhH313XxmftWbcxNRMlQChrvcXmbme705txC+dvFKc5uU9EcME1ZF5/Tu3gXzMofiVV1Nfc22Ctw6OsnkVgkzeX1DnZlziSgHKk18BNR8KisLgMXMvN6bcxLRJNgDnruUKYZKUxfCNOtLjuHilWYAQGpib9yU1t/kFoWXh6Ym4/d6EvWWA1WwVF1ESmIfs5slTNKhyhTMXGAoVpsKINUwp8qrIKUthwpQ66ECX2o7j7SOtFkIX3FMSX9EUtJ9bli/Xpg9xt6LMt5mFZGn0ylKzOyL1c4yAZQx81wfnEsIv9p6+CwsZy8BAGK6d8Hd6bLmlD88PC0ZRfvta1U9e9to9JGsyojUoR6VH8QBKDW7EUJ44o+f2D+b5WQOkZR0P7kprT9SE3sDUGOC6w1jgiKyBEug2giV0i5EUCs/ewmbDSnpD0tKut8QER4xpqpvl1T1SBUsgWoxgFQi+qHZDRGiPW8YxkpuGZWI5P69TWxN+Ls7fQhidI/VUnUJWw+fNblFwgzBcs+iFSpYrSCieQDWQWUQOiVzroQZ6hubkV9Sads2TkwV/tG7exfkZA7Ba59UAADe2F6BmaMSTW2TCLxgCVSlUGnnBJVYkam3HZHeHx24pgmh/HnXCVxoUCnpyf164eaR8oYZCA9NTbYFqo1fnkFldT2GJsjk6kjiVaAiolnOVvj1gSVwHpiECArMjNe3V9i2vzVlOKKiJCU9EEb0740ZI/tj66GzYAb+9NlRLLl9jNnNEgHkbY+qiIiqoSqYr2Pm3b5oBDO7K0YrhKmKj9Tgy1MXAAA9ukYhN2OoyS2KLA9NTcbWQ2p8au3Oo3g6ayR6dJUbK5HC22SKjVBLdCwBUEJEB4noh7r8kc8QUV8iSiaivr48rxAdZZxweufEwYjtJUtPBNKsMUkYHNcTAFBT34T39p40uUUikLwKVMxsXX5+IYBNUFUingdQRkQ7iOixzgQXInqWiM4BqAFQBmC+3v8BET3W0fMK0RlnLjTg/X/Z3xgfnOpy8QDhJ9FRhAem2Jeqf/1TqVQRSbxOT2fmOmbOcwhau6ESIFYDqNGBxZul6EFEH0CVUooHUI62Nf8mA8gjove9ba8QnfXOjko0tagh1Mzh8Rg3KNbkFkWmeZlD0S1avWXtqazF3mOy2Hek6NQ8KkPQyoAKMCugAkwWgAIiaiGid4jo1vbOQ0RPQNX42wUghZkd6/mNgOrBZUvPSgRSU0sr3vrsqG1belPm6denO752/TW2ban/Fzk6PeFXjyc9DjX36UfW3VC9ojoAc6GSMHa0c5oFUFl/s5wtc8/MtboHVwfVgxMiIAr3ncap8w0AgP59uuP2665x8wrhT8YPCv+35wSqLzWa2BoRKB0KVDrR4Vki2gk1nrQKqkdUDtWrStWV1BMAzIG6NZhBRD9zccp0AKXMfN7NpYvhfr0qIXzGmJJ+3+Sh6NYlWIq5RKaJQ+MwfrC69drY3Ip1xZVuXiHCgVf/64hoGREdgkp0WAEgA+p23RLYg9MSY0V1Zi4CMBuql5Xt5LRAO1UoHEiQEgFz8PQFfGqpBqAG8++/cZibVwh/I6I2vao3Pz2CFqn/F/a8/Xi4GGpNqF36+3hmzmTm59tb7oOZa6F6W0UuDtkFIJ2IXA4AENEIqEBV7OoYIXzJWNcve+wAXBPb08TWCKtvTBiEOD094FjNZWw5cMbkFgl/8zZQLUDb4FTn6Qv1QotLXTy9DKrHVUhEE4wvAwAimggVoBgqM1AIv7rQ0IQNpfZlJR6SJIqg0aNrNOZm2idcS1JF+PM2UKVA3e5rFxHdTUTPenpSZi6Fun2YBqBUz6ViAD/W35dAZxX6qYSTEG28u+s4LjW2AADSkvpgamo/k1skjL5143BYF1X+6GAVKvRCliI8deTW3yQPjrsXXvZ8dBmlNKg0dNKPeP11I4CMdnpkQvgMM7e57ffglOGy1HyQGdavF24xVFF/UyYAh7V2a/0RUSxUXT+jhUQ0uZ2XxUHNo/J6Nh4zW2BIuCCiWG9uLwrhC9vLzuHQmYsAgN7donF3+mCTWySceWhasm0Ry3XFlXhmzij06hYsC0IIX3L3W00AkGvYZqhkilQPzu0Y4LwmQUqYYc32Ctv392QMQUwPqesXjGaOTMTwfr1w5Fw9zjc04y+7T+C+yZKZGY7avfWnM/msgSkN6jbcKsM+V494uU0nQtHx2sso3Hfati1JFMErKorw4BT772fNtgowS6p6OHLbTzamnRNREYDC9lLRO0pn9s2D+7lSzMz3+vr6QgDAnz49Auu0nOlp/ZCWFGNug0S7cjOG4pcfHsTlphZ8eeoCdlbUYPKIBLObJXzMqxu6zDzHH40gonugSjABbYvROm0GVLKGED7V0NSCd3baKx08NDXZvMYIj8T26oo7Jw3G2ztUPcY12yskUIUhd8kUs/S3xcx8Xq875XF1CC9SyZdCBag8APnwvFKFED7z3t6Tttpxg+N6YvaYJJNbJDzx8LThtkD1wb9O4VRdAwbG9jC5VcKX3PWoiqB6MNlQaeO5UGnnnt4I9nQJznQAJcwsBWeFaRyXmu8SLXX9QsGYgX1x44gEfFZejeZWxls7juKZ7FFmN0v4kLtA9TxUULL2cEqhavz5mgXSixIm2l1Ziz3HVJJpty5RmHeDLDUfSh6elozPylVdxrc+O4rv3pomBYTDSLuBipkXO2xvhJp862sbAeQQUQwzX/DD+YVo1+vbKmzff2PCICT07mZeY4TXsq8dgIF9e+DU+QacvXgF//jXSXxzosx/CxfB8pFjEdRyIRsdav0J4XdnL17B3/bal5p/WJIoQk7X6Cg8YKhuv8bwwUOEPk+TKTrE02QKZq7Tqe/zoWr9uTmcvcpWJKIcqISNFKjKGaUA1uqyTc6OT9fHZ+ldRQCW6ZqEIsy8s+MoGltaAQCThsVh/BBZaj4U3Tt5GF7cdAhNLYzSo7X4/Fid/C7DhKfJFB3lUTIFEf0cwBN6sw5AdSeu6XjuRbDXHbSOhaVDLSuygJlTHY7PAlBoOB4AcqBuTeYyc4Gv2ibM19zSijc/tS81L72p0JUY0x13jL8Gf959AoBKjnk+V27QhANPkyn8LUd/zfJldXQddJZD1R2cbe0REVEcVBp8FhGtYuYFhv3WIJWtF300Bq98IorX62uJMOC41PxXx8tS86Hs4WnJtkD1lz0n8OOvjkW8jDeGPK+SKfwoAUCRH5bwWKC/PmG8bcfMtUSUCzUuNt9w3Hz9dYU1SOnji4goTz8/H/7JfBQmMNb1u1+Wmg95E4fG4fohsdh7rA6Nza1YW1yJhTM9KU0qglmw/K/016q96frrVSsL616RBbD1pABVwgkA1jo51yqHY0SI+/LUeYel5qWuX6gjojYVRd7YfgTNevxRhK5gqUyRB2AdEU1g5j2ent8DubodV92q08EpxeF56/ZVSRPMbE3ySHd8ToSmP/zTXrLyK+MGSjWDMPG166/Bsr/vx7lLjTheexkf7jstt3RDXLBUpjgMFaxKiWgVVFaey4QKZt7gyUndZOnl66/G23hxzg70BBFZbwti2DBZaiDYnb14xTaWAQDfvmmEia0RvtSjazQeuHEYXtx0GID6QCKBKrQFS2WKUn0dArAQrgMh6ec8DYBXn0AlRqyC6j0VORmHay9RohYughkz50GvwZWZmSlrDQS5tz47isZmdUtowtA4pA/r8GcUEYS+NWU4Vn5UhqYWRvGRGuyprMWEofI7DlXBUpliCfycXUhEKVAByjo3aoWLZJH2/prlLz0MXGluwRuGpcu/PT1ZlpoPM0l9e+Dr1w/Chl3HAQCvfVKOF+6dZHKrREd1eN1mIuoL1SuxjllZmHl3R87lauKtrzjMpSoAsFgve+/IAvdjcJKaHuLe23sSVReuAAAG9JWU9HD16PQRtkD1t70nseT2sTIOGaK8zvojor5EtBIqtbsEaqwnH0AJETUT0Q872yh9jWQdDDt7rlVQQaoUQCoz57oIUoA9C/CqhAnDPimeG8KYGb83JFE8NDUZXaVKelgaPyQWk5PV2lTNrYw3Pq0wt0GiwzryP7QUKmmAoG4DPq8fm/T5VhDR+x1pDBE9S0TnoIJgmb4OiOgDInqsA+ezJjgUMXNGOwHKyppg4SwF3bpvlZPnRIjYUV6NL06cBwB07xKF+ydL4ks4MybJvPXZUVxubDGxNaKjvApUutRRClSwimfmOcy8RD+yoSbu7gaQ7W3Piog+gOr5xAMoR9uVficDyOtAALSOQeV6crBOiACARcZelf5+kcMxIgT94RN7b+ru9CFStSDMZV87AEPiewIAauqb8Ofdx01ukegIb3tUWVBJD7OZuc7xST0faTZUkPF4YiwRPQGVAr8LQAozpzkcMgKqx5btac/KOE8KQDkR1bh6OLzUGtRKiKiQiAqhbnEanxMh6Oi5eny477Rt+9vTk81rjAiI6CjCI9OSbdt/+Gc5mCUpN9R4G6jSAZQy83lXB+hgVQogw4vzLoAKgLOYucLZOXWPrQ4qfd0TxqSIODcP47UKdNsLAGTqRwGADClIG9rWbK+A9T1qxsj+GDkgxtT2iMCYe8NQ9O6mZrQcOnMRWw+dNblFwlveBqpd8KwyhXUpDU+5DYBasYfXBzOXMjN58nDx2lxmjtePXFniI7RdaGjC2p2Vtu3HZIJvxOjboytyM+0rNhtv/4rQ4G2gWgsgnojucnUAEd0DPZnWi/N6mknncfkmIYwKSo7h4pVmAEBqYm/cPDLR5BaJQHp0ejKsU+W2HKjC4TMXzW2Q8IpXgUrPdyoCUEBEP9O1/wAAOp18GYB1AEqYeakXp94FtT6Uy6qgRDQCKlD5q4CtCFMtrYzXPqmwbT86fQSiomSCbyQZ3q83ssYOsG2/Jr2qkNJuoCKiFscH7MkSiwGUGfaXQWXGEYA4Inrci3Ys068rdFiKnnU7JkIFKIZ94q4QHincdxpHq+sBALE9u+Lu9MEmt0iY4dvT7bd715ceQ/WlRhNbI7zhrkdVAZUq7viwuNhvfc6rj6t6/GcJgDSowrTnoILSj/X3JVBp6yv8sGaVCGPMjN99VGbbvm/yMPTq1uGCLCKETUlJwLXXqBoCDU2tWLOtwtwGCY+5q/UXsBXHmHkFERVATdwnSOEAACAASURBVKjNgAp28VAlizZClT3aFaj2iPDwWXk1dleqqlfdoqMkJT2CEREWzEzBv72jKr2t2V6BBTNT5INLCPBL7RgiepyIXvH2dcxsYeZsZk5g5iioScUJemKxBCnhNWNv6p6MwUjqK7XeItkd46+xTQCurW+bCSqCV4cClU6cmOXicTfUbbwF7s5jOF+Ls6oTziYVC+GpfSfOY8uBKgAAEfDEDEkajXRdoqMw/2b738GrW8vRJCsABz2v+rxEFAuV9edulVuCd+npFQACdptRRIZVH9t7U7dfNxApiX1MbI0IFrkZQ/FC0SFU6xWA/7rnBO5OH2J2s0Q7vO1RLYUaPyqHKkRbofdbC9NugApS+cw8x4vz5gJI9UXldSEAoLK6Hn/dY1/Bd+FM+RwklJ7dovGooazSqo8sUlYpyHk7ipgDgK21+IgoD2oZ+Q+t2Xi6YvkyIoph5gsenrcaKljlEdG9UL2xna4O9nQpehG5Vm+1oFW/90xL7Yfrh8ial8LuwalqBeD6xhYcOH0Bmw+cwawxA9y/UJjC20CVAqDQusHMFr0y6iSoorFg5jzDQoXf8fC8FtiXos/QD2cfcTq9FL0If2cvXmkzSP7kLdKbEm3F9eqG+yYPs61NtnJLmQSqIOZtoHK2uq0FahkOo1Koauie8vtS9CJyrNlWgSvNaoB83KC+uCmtv8ktEsHosZtGYM22CjS3MnZW1KDkSDUyhieY3SzhhLeBygJVTdyoFFcnVxiX2HDL30vRi8hx6UozXt9+xLa9cGYqiKRckrjaoLieuHPSYBSUHAMArNxiwasPS6AKRt4mU6yDKkr7vi5rBOiK5tZCtbomXxZkyXZhgrd3HEXd5SYAwPB+vXD7dQNNbpEIZgtn2j9PF+0/jYOnPR1WF4HkVY9KV4+YB2AO1NpQ86AqSSyFKlRrgepJeVWTTxez9VShlFESzjQ2t+LVrfZio0/MSEGXaL/MaRdhIi0pBlljB6Bov1pQc9VHFvxy7gQ3rxKB5nXtEGbO0CvyVuvtOiKaDSAfai5ULYA8Zn7Vi9MudvO8cfzqHHTihhBGf959HKfONwAA+vfphpwMmRsj3HvyllRboPrL7uN4Zs4oDI7raXKrhFGHilwx82qH7VJ0bsJue4kXKQDmQlVtX8TMv+jEdUSYamppxUubDtu2H50+Aj26SnKocC9jeDwmJydgR0U1mlsZL28+jJ/dNd7sZgmDDldjJKK+UEHEepPXwsy7O3IuZt7YztMbAazW87NWElEJM2/uyHVE+NpQeqzNUh4PTXW5tJkQV/n+7JH41u8/AwCs21mJJ2emYmhCL5NbJay8voFPRH2JaCWAGqjlN/L1o0TX7PNLdQlmzoNaYFHWoxJtNDa34sWN9t7U/JtTENOjq4ktEqFmelo/TE5WGX/Nrdymdy7M15GR5lIA86Em326EvXzSJr1vhbMCsz5igZoMLIRNQckxHK+9DABI6N0NDxvK4wjhCSLCD7JH2bYLSo/hyLlLJrZIGHkVqIjo51C3+kqhluCYw8xL9CMbQAKA3QCy/dSzSofzScciQl1pbsFLmw7ZthfcnII+3WV9IeG9qan9MDWlHwCgpZXxW+lVBQ1ve1RZUBl4s50twcHMtbAvVT/P05MS0UQ3j1lEtBbACKh5W0IAUOMJJ+rsmX4PytiU6ARjr2pD6TGUn5VeVTDw9qNnOoASZj7v6gBmriWiUnh3i64U7ksoWcsLuEtlFxGioakFL222f+pdODNVVmsVnTJ5RAJmjOyPrYfOopWBFzcewq/nTXT/QuFX3v6v3gXPSiPFQQUfT62G+0BVBqCAmcvdHCcixNs7juL0+SsAgMSY7vjWFOlNic57OmsUth46C0DNq3rq1lSkJcWY3KrI5m2gWgvg50R0FzO/6+wAIroHKph5nJ3HzB6vBiwEAFxubMErW+wLIz51S6rMmxI+kTE8HreMTsSWA1VoZeA3Gw/jt/dNMrtZEc2rMSpdPLYIqlzSz4go2fqcXp5+GVQ9wBJmXurLhgph9KfPjqDqgupNDezbA/dOHmZyi0Q4+UGWfazqb3tP4MApqQFopnYDlZ4X1eYBe7LEYgBlhv1lABbp5+KI6HFvG0NEdxuK3VqTLA7pa+wgIinCJVDf2IyVxt7UrDTpTQmfmjA0DlljkwAAzMBvNh40uUWRzd2tvwp0bJ0or9ZVIKJYqMnDI6AC4G69bxPUeBeglhcpIaI0Zq7oQJtEmHh9+xGcu9QIABgc1xNzM6Wmn/C9p7NGoWj/GQDA3z8/hX0nzuPaQX1NblVkardHxcypzJzWwYc3RWmXQ41rbYS6tQio+n5xAJYzc5TejoLqtYkIVXOpEa8YMv2+OysN3btIb0r43nWDY3HbOPuqv8v+sR/Msr6rGXyyBoKu+9cZWQBq9ARia73AXKje3DIAYOYCqMoU3qwcLMLMC0UHcb6hGQAwon9v3JMuvSnhP89kj0aUvj+09dBZbDlQZW6DIlSHApUeO/qAiM7p8aka/f37RHRrB06Zgqsn8mZCFbo1ztmyrnclItDhMxfw5mdHbdtLbx+Dbl1kvSnhP6MHxrRJ1HnuvX1oamk1sUWRqSNFaZdBjSdlA4gHUK4f8VALKhYR0StentYCVX7Jeo1JULf9ihyOS4GUUIpYP/v7l2hpVbdepqb0Q/a1A9y8QojOeyZ7lK0sV1nVJby946ibVwhf87bW3z1QyQ51AHKZOcowJmUdRzoPYIF1aXoP7QKQbuiNLYW67ZdvuPZsqEAlS9xHoK2HqrDpSzWwTQT8+9fGgsirnB0hOqR/n+546tY02/avCw+irr7JxBZFHm97VAugAkg6M693fFKPI2VAZf15M4l3mX5Nkb6VeA/Ubb9NAKCXFfkQhjErETmaW1rx3N/227ZzM4Zg3KBYE1skIs2j05MxJF6t+ltT34TfGgohC//zNlBlAihtLz2cmS1Q5ZNu8PSkeoXgTKjK6wTVw5pjOGSO3r+EmTd42WYR4tYVH8OB02rCZa9u0Xh2zmiTWyQiTY+u0Vh6+1jb9prtFaiQgrUBEzQj0cxcyswZ+nZipkNNvyy9/3nTGihMcaGhCb8qPGDbfnJmKpL69jCxRSJSfXX8QGQOjwcANLUwlv1jv5tXCF/xNlCVQI0luaz+qSfqpsOHy3FYgxYRPa5vA4oI8fLmMpy9qCb3DortgSdulqRPYQ4iwn987Vrb9gdfnMb2snMmtihyeBuoVkHdgit0Vs5Ilz8qhkMihKd0vcBZLh53A1gCtbqwiACV1fX4wz/tHevFt4+RUknCVBOGxuGuSYNt28+9tw+trTIJ2N+8qp7OzAVEtB4q2aGUiGph7zllQqWUE4B8bypT6F5YEVRPrN1DcXXKughTP//Hl2jUc1YmDI3D168fZHKLhAB+dNto/ONfJ9HQ1IovTpxHQekxzM0canazwprXY1TMnAuVhl4BNXcqG23nVOUys8er+2pLobIFywE8r88N/f3zADbAHgDnODuBCC+bvzyD9z4/adv+f18bi6goSUcX5hsU1xPzZ9hvQS/7+36cu3jFxBaFvw4lUzBzga4DGAUgFUCqYU7VVWnrHshRp+U0Zl4CFfgIwIfMvEQHx4UAsohIVjALcxevNOMn735u275r0mBkDE9o5xVCBNaCmakYHGdPV//vv+0zuUXhzdsJvy1E9L5xHzOX+2DV3RQYbunpFHcAmGTYlwegBl4syChC0y8+OIATdQ0AgITe3doMYAsRDHp374Ln7rrOtv2X3SewWU9IF77nbY+qAqoH5WvOyiJZAEx22FcKKUob1kqO1GDN9grb9v/72rVI6N3NtPYI4cqto5Nw50T7uOlP3v0cF680m9ii8OVtoMoFkEpEP/RxOyxQyRhGpbg6uSIOUpQ2bF1pbsHi9XthXUnhltGJ+OZESaAQwes/vnYt4nt1BQCcqGvALz444OYVoiO8DVTVUMHqx0S0k4iW6VV5nT68OO86APG6+rp1hd9iACnWmoFENAJqORCp9RemVm4pw+EzFwGoChTP3Xmd1PMTQa1fn+74z6+Ps22v2V6BkiM15jUoTJE3C4ERUSvUHCnju4ezExBUcoTHk16IqARqTCqfmefplPUKAH3RdnmPBV4uyhhwmZmZXFzss/nOEeHg6Qu448WtaGpRf07/+fVr8ej0ESa3Sgj3mBmP/nGnba2qkUl98Lfv3yQLenYAEZUws+PdNe/mUUFNuPXL7DZmziCiJ6B6bWDmOl0xPR9qXKwWQF6wBynhvZZWxuL1e21BatKwODw0NdncRgnhISLCc3dehzm//hj1jS04dOYiVm4pw9NZo8xuWtjwdsLvCn81RJ9/tcN2KfyTvCGCyBvbK7DrqMqn6RpNWH7P9YiWOVMihAyJ74Uf3TYaP/2rSlN/efNhfHX8NRg1QGbT+ELQFKU1IqK+upxSZ5e4F0Gu4uwlrDAMQD95S5r85xYh6aGpyZg4NA6AKlr7o/w9aGyW1YB9waNApQPHLF0Udpa/AggRPUtE56DmS5VB1/XTy94/5o9rCvM0NLXgqbdKUd/YAgBIS+qDp26VDrQITdFRhBU516NrtLobsOdYHZ7/4EuTWxUe3AYqnb1XDqAQqihtIYBqL1fwdYuIPoCazGstxWS89zMZQJ7jZGMR2n729/344sR5AEC36Cj8eu5EGYAWIW3UgBj86Db7emmrt5Zj4/7TJrYoPLQbqIhoEoAC2IPHRv01CkCBswrqHaGTKLKhFkxMYeY0h0NGANgEIFt6VuHhH5+fxOvbj9i2f3LHWIwfIqv2itD3+E0pmDUmybb9w/w9OFF72cQWhT53Paql+usqXYdvjg4iq6F6PL4qZ2Rd4n6Ws9WDmbmWmbMB1EHV/BMh7Oi5eixav9e2/ZVxA/HQVJdLnAkRUqKiCL/InYCBeoHP2vomfP/tXWhqkfGqjnIXqNKh5kM9adzJzAv0tx4vN+/BdUqZ+byb44ohlSlCWmNzK773dikuNKhSM0Pie2J5zvUysVeElYTe3fDb+yfZsleLj9TgV4UHTW5V6HIXqFLguhKEBaqkkS94Wm1CglSIW/7+l9hzrA6ASkV/6f50xPbsanKrhPC9G5IT8Ey2fS7Vyi1l+OhglYktCl2eZP05Kxjb3v6O2AX3S9yPgApUUvIhRBXuO43fG1fs/coYWzqvEOHoyZmpmDGyv237mbW7cfp8g4ktCk3BMo9qGZwvcc/AVUvcyzIfIejouXo8m7/Htp01NgmP3SQlkkR4i4oi/HreRCTFdAcAnLvUiO+9vUvmV3kpKAKVrkCxBEAa1BL356CC0o/19yVQmYcrmHmTeS0VHVFzqRGPvLYDdZebAACDYnvgF7kTZFxKRIT+fbrjN/dOgrXYyo7yaixZvxfe1FmNdEERqABbeaY0qDR00o94/XUjgAxmXur6DCIYNTS14PHXi2E5ewkA0K1LFH57fzrieskaUyJyTE3thx/Osc+v2rDrOH75oSRXeMqTWn/pRNTi6sl2nmNm9raWoAWGhRGJKJaZ67w5hwgeLa2Mp9/ZbVv2gAh4Yd5EZAyPN7llQgTed25JRWV1Pd7ZWQkAeGnzYQyK64n7bxxmcsuCnyc9Kurgo9O9NQlSoe259/bh/S9O2bb//Y5r8dXx15jYIiHMY62yfuvoRNu+f//z59j0pVSucKfdHg8z++XWIBGt7MTLmZm/47PGCL94dasFr31SYdv+9vQRkjwhIl6X6Ci8dH867s37FJ8fr0MrA0/9aRfWLpiC64dIBqwrXi2c6LOL2hdgBNrW9POEVwsymiHSF058b+9JPPVWqW379usG4uX70xElS3cIAQA4c6EBd7+yDcdqVGml/n26YcOT0zGsXy+TW2YuXy2c6EsEFawKoRZHLG//cBEKth0+ix+s223bzhwej1/PmyhBSgiDpJge+OOjk3HPym2ou9yEsxdVZuw786cgSZdeEnZm9ahyoJbwyNK7jAEr34NSSkEtUntURftO4ztvldrmiKQk9sb6hdMQ31sy/IRwZmdFNR549TPb/5nh/XrhzcduxNCEyOxZuepRmZKezswFusBtFIC5UCnpc6CK3dYQ0ftE9JgsnBg6/m/PCSx8s8T2H25A3+5Y8+hkCVJCtOOG5AT8Zt5EW03AI+fqMXfVdpRVXTS5ZcHF9HlUOmhlS9AKXW99dhT/9s4uNLeq3vmwhF4oWDgtYj8VCuGN28dfg5UPpKNbtHo7PlnXgLm/244vTkjSs5XpgcrISdDagLZBS1b6DTJ5H5fhx+9+Dusd5JFJfZC/cKoEKSG8MGfcQPzhkRvQs6vKEzt3qRH35n2KkiPVJrcsOARVoDLSQStXB61cqOoU2VCrDAuTMTN+9eEB/Ozv9qW2xw+OxdoFUzFABoOF8NpNI/vjzcdvRN8eKsftQkMzvvXqDmw9JBXXgzZQWRHR3VC9q9l6V9C3Odw1NLVg6YbP8eKmw7Z9k5MT8NYTNyJBxqSE6LCM4fF4e/4U9NP/jy43teCxPxZjXXGlyS0zV1C+6RPR3US0VpdnyofqUe2CWglY6u+YqLK6Hjm/22YrAwMAM0clYs23JyOmh6wrJURnjRsUi3ULp+KaWHVnorGlFYsK9mJxwV40NLmsZhfWgiZQEdEsF8FpCYB4Zs5k5tVSVsk8RftO444Xt+Jfx+2zB745cRBWP5SJnt2Ceg62ECElNVGN9Y5M6mPbt7a4Ene/sg1Hzl0ysWXmMGUele3iRLOgekk51l0ASgGsBZAXqkEp3OZRNbe04peFB7FyS5ltX9dowr/fcS0emjpclusQwk/qG5vx4w2f48+7T9j2xfTogl/kTsBt4waa2DL/cDWPyqwJvyuhxp3i0DY4FTBzyFeoCKdAdeZCA77/9i58arFnHw2K7YGXH0jHpGFyF1YIf2Nm/Omzo/jvv+5DY4t9wcX5N6fgR7eNRtfooLkx1mnBFqiMtf6KoAKVp5iZf+z7VvlOOASq1lbGuuJK/Pz9L1Fb32Tbf/OoRLwwb6IkTQgRYHsqa/GdP5XieO1l274xA2Pwv3eND5ulc4IxUHWUFKX1s/0nz+Mn736O0qO1tn1EwNOzR+G7s9Jss+iFEIFVc6kRz6zbjc0H2qas3zd5KBZ/ZUzIL0gabIHqns68npnX+6ot/hCqgerSlWa8UHQQf/ikAi2t9r+LIfE98fO7r8dNI/ub2DohBKDudrz6Twt+VXgQDU32z/wJvbvhx18di3vSB4fsuHFQBapwF2qBqrWV8fd/ncT/vrcfJ+sabPu7RhPm35yC7946UrL6hAgyldX1+Olfv0DR/jNt9k8ekYD//Pq1GDco1qSWdZwEqgAKlUDV3NKKv+09iZc3H8ahM22LYE5JScBzd16HtKQYk1onhPDEh1+cwn/93xc4YfiQCQCzxyThqVlpSA+hpCcJVAEU7IGqsbkV7+46hle2lOHIufo2z/Xr3Q0/uWMs7poUurcPhIg09Y3N+M3GQ/j91nJbcWir6Wn98N1bR2JKSkLQ/5+WQBVAwRqoausbsaH0OF7darnq01fvbtH41tTh+M7MNMT2kgoTQoSig6cv4IWig/jHv07B8a09Y3g85t+cglljkoI2pV0CVQAFU6BqamnFxwersL70GIr2nWkzDwMAYnt2xaPTk/HItOSQzxgSQiiHz1zAK5vL8Jc9J9okRgEq6eKbEwfhnvQhGDeob1D1siRQBZDZgYqZsf/kBWwoPYY/7z6Osxcbrzqmf59ueHxGCr41ZTj6dO9iQiuFEP525Nwl/O6jMhSUHENTy9Xv9WMGxiAnYwi+MWEQkoJg1QMJVAFkRqC60NCEbWXnsOVAFT4+WNVmUqDRhCGxyMkcipz0IZLJJ0SEOFl3GW9sP4J3dx1vk9lrdN3gvrhlVBJmjk7EpKFx6GLC7UEJVF4gonQASwFk6V1FAJYxs0cVNAIRqC43tuCLE3XYWVGDjw6eQXFFzVWDqFZJMd1xV/pg5KQPwcgBksUnRKRqaWVsLzuHgpJKvP/FqTbzsIxienTBjJH9MWNkItKHxSMtqU9AJvpLoPIQEWUBKNSbFv01RX/NZeYCd+fwdaBqamlFWdVF7Kmsxe7KOuyprMWB0xeuuvdsFNO9C2aOTkROxhDclNbflE9HQojgdaGhCX///CTe3XUcOytq2n0/6dUtGuMHx2Li0DhMGBqH64fEYnBcT5+Pb0mg8gARxQGo0ZvZzFyk9xuDVzwz1zp7vVVHA1VtfSPKqi6irOqS+nrmEixnL+LouXqXvSWjcYP6YuaoRMwclYj04fFBm9kjhAgu5xuasO3wWWw5UIWPDla5vD1o1KtbNFISeyM1sQ9S+vdBapL6fkT/3ujRtWPDCq4ClYyitzVff11hDVIAwMxFRJSnn58PYIWvL/zqVguee2+/x8cTqTVrJgyJw5SUBMwclRgUg6FCiNDTt0dXfOW6a/CV664BM+Pg6Yu2IYXdlbU4c+HKVa+pb2zBv46fb7M+HaDem97/t5sxeqDvhhkkULU1T39d6+S5VVBBah78EKiGxPds9/lBsT0wfkgsJgyNw8QhcbhuSCz6yoq6QggfIyKMHhiD0QNjMP9mte9UXQN2V9Ziz7Fa7KmsxRcnzqPucpPT1zMDQxPafz/zlgSqtlIAwFnSBDOX6vux6f64cGpiH3TvEoWUxD5IText+5qa2Acpib3Rq5v8qoQQ5hgY2wNfiR2Ir1ynFmtkZlRfakRZ1SVYqi7ahiwsVRfRwuzz9ysZozIgIgYAZnY6QujueauOjFG16jGoKFlCQwgRwlpaucMZgq7GqGS0/WrtJUq4fI6I5hNRMREVV1VVuTrMpagokiAlhAh5/khjl0B1tbiOPMfMecycycyZiYmJfmiWEEJEJglUbVncH9Juj0sIIYSPSaBqywLYKlO0YdjnSTATQgjhIxKo2srXX+c5ec66b1WA2iKEEAISqNpg5jz97SJjr0p/v8jhGCGEEAEggepqufprCREVElEhgBKH54QQQgSIzKNyorPV04moCsARPzUvlPQHcNbsRog25HcSnOT3ogxn5qvSpiVQCb8homJnk/eEeeR3Epzk99I+ufUnhBAiqEmgEkIIEdQkUAl/kgzJ4CO/k+Akv5d2yBiVEEKIoCY9KiGEEEFNApUQQoigJoFKCCFMQEQ5RFRCRDVExPr7RWa3KxhJoBIizBFROhHl6zfEGv29X1aqFp7RASkfasXwagCl+vvlRFRmZtuCkQQqETDyCTLwiCgLqgRYDtQbYrX+voSIcsxsW6TSv5PlUEsGZTBzKjNnAIiHqoKTQkRS/NpAsv5EQOiAtFxvWqD+k9qWTmHmVFMaFsaIKA5Ajd7MZuYivT8LQKHeH8/MssZaABFRPtSHhVxmLnB4zvY7Y2ZZ8luTHpXwO/kEaZr5+usKa5ACAP19nsMxInCsH9CKHJ/QHxqs6+K1t9p4RJEelfA7+QRpDiIqgXpTzHAsqKzHqEoAlOoPDSJArOODzopcy/8H5yRQCb/Tg8MpcHGbyd3zomOIqAZAnKs3PCJiQN4Qg4leVigLqhe82Oz2BAsJVMLv5BOkOdwFIglUwUPfHl8F9YGtiJmzTW5SUJExKuF3zFzazlpe+frrikC1J8K010OV3qvJiChF96IKoYLUCglSV5MelTCFfIL0P+lRBTeHTNgCAIuZ2WJik4KW9KhEQMknyIDy5E1PelUm0Fmuy6Em+qYyc64EKde6mN0AEfw6MFO+1lkmmXyCDDgLVOp/uousP+sxIoCIaD7UtAC5k+AhCVTCLV9MxtWfIOdDfYKUT4+BkQ+VQTYP6t/daJ7+KvPXAs+azZdraitCiIxRCb/TnyBXQT5BBpx1HAqGuVSGOVQyPhVgDtVC2r3tyszx/m9RaJBAJfxO5kmZR9fzs2ZWWishZOmvV03AFv5l/JDgjnyIsJNAJfxKPkGaT785LoU9QBUBWNbOlAEhgooEKuFX8glSCNFZEqiEEEIENZlHJYQQIqhJoBJCCBHUJFAJIYQIahKohBBCBDUJVEIIIYKaBCrhN0SUQ0SsH25L9Xh7vAhuugCxLKcuOk0ClQiUuR4cs8DvrRABQUQpAMoArDa7LSL0SaASgRKn16Bqj7vnhRARSAKVCARrjTmXPSZdk854rBBCAJBAJQKjVD9y2jnGGsTyXR1AROlElE9EZXocq0RXZnd1/CJ9TI0+voyIljsbN3Fy7jIiWuV4rN7HhvWcjM9l6eeWG/Yt0vtS9DXKDBXNvfq5DOdK1+N5hfpnK7Nek4jidBuN53L67+7Nv6fDtdP1tVlfv9D476EXxrSuYWYdd1zu7Lze0O1jZz3z9p4ToU8ClQiUtYBtCXpnsqB6U07XqdKLLpZABbtaqMCXDmCVfmN0PL4EapHGFADF+twpABYB2OhwrLUeofXcRQASoNbP8qhOoQeMNQ9tvUZvfy5tAewBvRj659IJKCVQ44HWDwfpAPIdA2sHrwuo31OJvmYBgGrrPj0uBaglXfL09xao9ZfaO6enlumvi407iShft30BM0uPPBwxszzk4ZcH1Jsgwx4wGEB+O8fNh3rTYwCrDM9bX1sGIMWwPw7qTZMBLDLsn6/3FTq5Vpl+znge6zmyHI4t1PtzDPtW6X3pTs5tbftyw75Fel+NsY0d/Lms52pzfcN1Wb8uzkl7l3f0uk6u7fhcfjvXuOr33cm/qRrj78/QrlW+vI4P29vm39jw71xj/D3Jo/2H9KhEQLBa0dfV7T/rbb91Ll5uTVVvszIwq7WtrKukLjUcXw291L2Tc1nXX0ox7EvX53P8NL5YP3yxGrGFmVc47PP257LKY8MSHbrd1iVUnuC2a35Ze17Gn7ej1wWA0nZ+jhTHg/3A1qvSvfPlUAtythn/1LcjTZ3ioG+5pkD9PRrNhQpSsjabh2QpehFIawGkE1EOt12wLwvqzaaWyOlKMXhr3QAABFhJREFUH5kAatnJ+knMbCEiC4AUIopj5lp97qsWBNS3ppwFylLdrnwY1mnSX321ZpOzW1Je/VyGp53djqyGevNzPJezINvR6wL6Fq4H13BK33Ls58GhO9n5oo55UMFpvn6UsvNVo0tdtDWQsqE+oDj+G2ZDkoa8IoFKBJL1TWYedCAxDPS7TKKAulViXFbdlQTonoVOgpgP4AaoHlN7n/ZzoW7z5UAN/tdCjf3kA1jno0++O53s8/rn0hw/oXuro9cFOt+7XGq9vhtOP2zoDzN5UL/bWgCznb2YmZ31pgPNOu7qbH+ek/3CBQlUImD0m4zj7b95+qvT236GrLta2G/7uFKtX5MOlTARB/XGWgp1e6oU6tPsIod2WQCk6ltJuVBvJNbHciKa7az30Rkd+blC+bpW3MlVnHX7bQk5zj5E6N/jAmbOdXzOybFl7o5xUMvMGR62MwUOSSS6Vx/nuF+0TwKVCDTH23850Lf9nB2sg1ut/t5xbMSVfKg3g1zH20dE5Ow2kfVaRdCfgPUbzVKooJYPINWD63o8RtPBn6vTzLquD22E+ne23p50vI0MqA8jHv0umNmT32tHWIOpY48qS19Xbv15QZIpRKBZb3nM8/C2H6Buw8W5mLsUZ53LY9idAjU24GyMo80bmJ7fVOY48K7HuhY7e42W4GSf20/wDrz9uXzFrOuC7HPa3D2u+pswpqGj/aSPdKif0UyuxqcyoMc99d9eIBJQQp70qERAOdz+s96GcpXtZ7UYKoFgIxFlGDPVYO89Gd/YaqE+bacYj9WTWa3BMU63x6LfLOYTUb7xk64hkBpv+1lvFS1A2/lQ1tR6b3j7c/lKIK/rOB61DB4mUxg3dBJGDlTGY57eZ02CSXe4NZsJeyaiWbJw9fgeoDL+rH/vOSHaqw08s/Pj5RG+DxjmUTnst85zumquE5zMo9L7lxteUwZ1j7/GxTmMxxZCvemW6ePzDeeYr49f5OTcZYZ9xjlLKQ7H5sMw98jx5zXsy3Hxb+TNz+XyXNb2OtnvdD6TN9f14NpXXQMqQFnPn+/q5/fy78jx3yOnneumdPR6Pvi7t/3sTv4Ny6DnoTn+v5CH64fc+hNmMPagPPrkzuo2nDWtNwHqU7MFatA828mxC/TzWVC3gooAjGA1wG6tUpGhj18BdSupVJ/b2jPKA5DKbecsWfTrrO2w9rqy9Xm8qsLgzc/lS/6+LqtbXiugehXW+URe07cn83Xb2txaZXVrtxYqU9Pac7OOAfli7ltHWf9+VuhyVot0jzAP6u/yBr1tdq8vZJCO7EIIEfJ0TcEs9iAzz49tWKXb4K9EjYgjPSohRDgJhkQKV/OnRAdJoBJChJNM+K6QsNdczZ8SnSNZf0KIsMGdnFDsA5n6q/SofEh6VEII4TvFUIk1UnDWhySZQgghRFCTHpUQQoigJoFKCCFEUJNAJYQQIqhJoBJCCBHUJFAJIYQIahKohBBCBLX/D7VoCGLdKshCAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "x=np.linspace(-3,3)\n", | |
| "y=norm.pdf(x,0,1)*100 # convert fraction to percent\n", | |
| "plt.plot(x,y)\n", | |
| "plt.xlabel('Measurement=$x_i-\\mu$')\n", | |
| "plt.ylabel(r'Probability of\\\\ Measurement (\\%)')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Half of the measured values are above the mean, $\\mu$, and half below. To determine the number of times any given range of values would be measured, you can integrate the probability density function between two measurements, for example $\\mu\\pm\\sigma$ the mean plus or minus the standard deviation. This distribution is the [cumulative distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function) shown below." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, 'Probability \\\\\\\\ Measurement $<x_i$ (\\\\%)')" | |
| ] | |
| }, | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAEcCAYAAACrjpH4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3xU5Z0/8M83N8IlIYkgIHILKGi5JsHe14rBWteuigngb+3arQquv273t61ysa27blu5WLvbvWgTtNdttxLQWltXJWi3a1urSbhYFJGEgMqdJBAIuX9/fzzPmUyGmcxMMplzZvJ5v17zmpwzZ855kkC+8zzn+3wfUVUQERElixS3G0BERBRLDGxERJRUGNiIiCipMLAREVFSYWAjIqKkwsBGRERJhYGNiIiSSprbDSBgzJgxOnXqVLebQUSUUKqrq0+q6tjA/QxsHjB16lRUVVW53QwiooQiIgeD7fdUYBORRQDyAUwHkAOgFkATgCpV3TmA864AUKeqlTFpKBEReZbrgU1E7gKwEkCBsyvIYSoiTQC2AShT1VeiOH8OgDIA5QBCBjYRKQCwFkCx3VUJYJ2q1sTieCIiig9XApuIZMMEhVUwgawOwFaYwNVgtwEgD6bnthAm8C0FUCoidQBWqeozEVxuUwTtKbbXht+1SwCUiEipqm4ZyPFERBQ/cQ9sIrIEwBN28xGYHtiBMG/b6vf+Ypge3lYReQnAUlU9E3CNEgCLYQJhTpj25KAnSC12hiv9gleFiOSqalN/jiciovhyI93/CQCrVTVPVddEENR6UdVKVS2F6c3Vw/T8Aq0FsAJhgpq1wj5v9L8HZ78uDzimP8cTEVEcxT2w2YAWdngwgvM0qeo9qnpBYFPVQlUVVRUApWFOtcw+PxXktbKAY/pzPBERxZHrySOh+N2HKwBQDWB94JBjjOQDQLCkD1WtERGgJ7GlP8cTESWsjq5unG3tRHNrJ860duBsWyfOt3ehpb0LLe2dON/R5dtu7exCe2c32jq7/Z7Nvs5u9T13dnWjvcs8Tx87Ct/7XGFM2+zJwCYio2GSMnIB1ABYA2CFiOQPQnCLZLhyIMcHZacgrACAyZMnx+KURERhtXZ04diZVhw704ajZ1pxorkNjefa0dDSbp7PtaOxpR2NLR1obu1Aa0f3oLYnNSVYIvzAeDKwwWQy5gCYrqoHRCQfwH6Ye1jLB+F6fSV6NOHCYBbt8RdQ1XLYe3JFRUVcxpyIYqKtswvvNZzHoYZzOHiqxT7O4XBTK441t6KppcPtJvbS2R37P39eDWzFAMqdxBJVrRORTTAp9YOhr0AU7LVojyciiqnOrm7UnzqHt480452jzdh79Az2Hm3GB03noTGMFSkCZGWmY9SwNGRlmseIjDSMyEjF8PRUDM9I9X09LD0Vw9JSzHNqCoalpyAjNQUZaSlIS01BeqogPTUFaSnmOT01BcPTU2PXWMuNdP9FMFVA6sMcmhuwPVi9mjrY+2Z98O+hRXs8EdGAqCrebzyPmkONqKpvxI73GrHv2Fm0d0Y/TJiWIrg4axjGjc7EuKxMXJw9DHkjM3DRyAzkjsxA3gjznDsiA9nD0zA8PRU2dyBhuNFjmw6gUkS+B2BNiHtmFQDuEpF1qrpLRBbA3I8qC3LsQNUByBeRgsCEEFtdxDmmv8cTEUVFVbHv2Fm8uv8kag42oupgA46daYvovSLAJaOHY8pFI+xjJKbkjcCkvBEYl52Ji0ZmIGUQ7mt5SdwDm6puEpEDAL4HYKWIbFDVBwIOWwUzHFljq4zkwwSL1YPQpAp7rWUwiSr+nLT9sgEcT0QU1rm2Tvxu/0n8Zt8J/GbvcRw+3Rr2PeOzMzFrQhZmjc/GrPFZmDUhC9PGjMSwtNgP7yUS0VgOxkZ7cZMZuB5mmHGVqj4Z8PoqAEUA3lDVR/p5jRKYYFSuqitDHOP8EAqdXpjtfVUDgJ0P1+/jwykqKlJW9ycaek6dbcOvdh9B5dvH8Me6BrR3hR5aHDUsDQsm56Bgci6KpuZizsTRyBmREcfWeo+IVKtqUeB+V5NHnMxAG8A2icgaACucIsequjFOTSmFCX7VIuJUEyn2e22gxxMRAQDOt3dh29vH8IsdH+B/9p1AV4iswKzMNPzZZWPxkfw8FEzJxazx2YOSGp+MPJEVqaobRaQcwAYA220NyJWqGnStnUG4/hYRKUTvav1bEKJaf7THE9HQpqr4Q+0pbK35AC/86QjOtXcFPW7W+CxcM+tiXDPzYhRMzkFaqhtVDxOfq0ORwdg5a2UAFtnnUAkmSYNDkUTJqbWjC8/s+ADff/UA3j1+NugxRVNycdP8S1B85ThMGD08zi1MbJ4birRp/2UwiSG1MD20V1S1DsBiWy1/PYDGEAkmRESedPxMK378h4P46R8PojHIhOj8sSOxZMFE3DR/IibljXChhcnNrfXYrgXwEoAdMOWylsNMASj2u79WCaDIJn+Ui8hKBEkwISLyiv3Hz+I/XtmPX+0+jI6u3qNhIzNSUVJ4KUoKJ2H2xOyEmxuWSFwZihSRKgCjVfUyv321AE6p6lUh3rMK5h5W0uWxciiSKLEdbjqP71a+i4rq9xCYC3Jp7nB8/mNTsXThJGRnprvTwCTltaHIAlw412sbgLtDvcEmmHB+GBF5RsO5djz2yn78+LWDF1QBWTg1F3d+YhqKrxjHJJA4cyuwBStL5UzCDklVTw9ai4iIInSurRNPvnoA5b+tw9m2zl6vffKyMfjy4suxYHJgVUCKF7cCWzmA9SLyMIDNMBU7imEqjhARedb2t4/hwWf34IOm8732z5uUg9WfnomPzRjjUsvI4Upgs8OK02ESR1YDEJjKIN92oz1EROEcO9OKh57bg+ffPNpr/4yLR+G+62bi0x8ax4QQj3At3V9VV4rIepj7bTXOEjVERF7S3a346euHsPG/96LZb9gxb2QGVl8/E7cWXMp7aB7jxrI12c6EaxvMBhTQ/M9HRBRL+441Y83W3ag51HslqpLCS/HVG65A7sihXavRq9zosTXZntr6gQQkO8F7A4BKmNJWREQxoar4z9cO4hu/frtXtuO0MSPxrVtm42PTeR/Ny9zoP88AcB1MRZGfi8g1kb5RRKaKyH0i8i7M9IBKVWVQI6KYOd3SgXt/WoOvP7vHF9TSUwVfWjQD//13n2RQSwBurMdWB1NRZAVMFuRSuwxMDYAqmPJaTr8/B8BFMFMBiu22wBQcvo735YgolqoPNuJL/7WjV8bjrPFZ+NfbFuDycVkutoyi4WbyiLNkTTHMUi/XAgi6XppVCdNLK+d8NiKKpe5uRdlv6/Dtl97ptYzMX310Ch644QpkpiddwaOk5vqyNbYmZCUAiMhomN5ZPoA8mAnbDaq6w70WElEyO93Sgb/9+Q78dt8J377szDRsLJmL62dPcLFl1F+uBzZ/tie2wz6IiAZV/clz+MKP3kDdiXO+fQWTc/Dd5QtYdT+BeSqwERHFy+sHGrDyJ1W9lpX5m09Nx5cXX450zktLaAxsRDTkPF3zPlZv3e1bWmZYWgoeXToPN869xOWWUSwwsBHRkNHdrfjOtn3491f2+/aNGZWBTX9VxKLFSYSBjYiGhNaOLnylYhd+vfuIb9/McVl48vNFuDSX99OSCQMbESW9lvZO3P3jKvxu/ynfvqsvH4t//z8LkMXFP5MOAxsRJbXm1g7c+cMqvF7f4Nt3x0en4Os3XsnixUnKc79VEXlcROa73Q4iSnynWzpw+5Ov9wpq9396Jh66aTaDWhLzYo9tOszCoztFZDOAaTCltrbBLG9T72LbiChBNJxrx+1P/BFvHemptf61P78Cd30y38VWUTx4LrCp6nV+m2UAFsOs2bYMQE5AXckqVX0y/q0kIi873tyK25/4I/YdO+vb942bZ+NzH5niYqsoXlwPbCKySFVfDvaaqm4HsN3v2NEwga4IJtgtBcDARkQ+R06fx19u+iPqTppqIiLAhlvnYmnRJJdbRvHiemADsE1EVkTS87Ilt7bYBxFRL6fOtuEvn+gJaqkpgu8snYeb5k90uWUUT164e7oTpsr/Y+EOFJHsOLSHiBJQc2sHPv+DnrqP6amC//g/CxjUhiDXA5uqFgJ4GsA9IvJCsGNEZL6IPAWgIdjrRDS0tXZ04e4fV+HND8yKVikCfHf5AlbnH6JcD2wAoKqlANYAuE5EXheRLAAQkSUi8gaAapg12870cRoiGoI6u7rxxZ/twGt1PZ97H75lDm6Yw6A2VHnhHhsAQFU3ikg1TFp/vYgAPStm1wBYp6pbXWwiEXlMd7di9dY3Ufn2Md++NZ+ZheVXTXaxVeQ2zwQ2EZkKoMRuOtVIawGUqupON9pERN6lqvjmr9/G1pr3fftWXp2Pe66e7mKryAtcH4oUkUUi8iJMEFtpd5fDrJ6dD6DQrbYRkXc99ptafP93B3zbyxdOwprrZ7nYIvIKL/TYKu1zE4B1AMpV9bSds/YyTMZkoare61oLichTntt1GI+8+I5v+zOzx+Nbt8yBvYVBQ5zrPTaY+2elqpqnqo/YuWpQ1dOBGZNM9yeine814b6KXb7tj02/CP+yfD5SUxjUyHA9sKlqUV9JIf4ZkzBltIhoiDrcdB53/agKbZ3dAID8sSPx+O2FGJaW6nLLyEtcD2yRUNWNMOWzZrjdFiJyx7m2Ttz5oyqcPNsGAMgZkY7v37EQo4dzPTXqLSECGwCo6hYwkYRoSOrqVvzdz3fibVupPz1V8L3bCzF1zEiXW0ZelDCBDQBUdYfbbSCi+Nv4wt5ec9W+dfMcfCT/IhdbRF6WUIGNiIaezVXvoey3db7tFX+Wj6ULWamfQmNgIyLP2nGoEV995k3fdvEV47Cac9UoDAY2IvKkU2fbcO9Pa9DRpQCAWeOz8F2m9VMEGNiIyHOcZJEjp1sBAKOHp2PTXxVh5DAv1JQgr2NgIyLP+ZfKfXh1/8me7WXzMSlvhIstokTCwEZEnvLy3mP4t5f3+7a/tGgGrpl1sYstokTDwEZEnnHoVAv+3897FvP45GVj8HfFl7vYIkpErgc2EblWRJZEeOwSEVk02G0iovhr7ejC3/y0GmdaOwEAl4zOxHeXL2CyCEXN9cAGU/9xbYTHrgXrRRIlpX94dg/2HO6pLPLY7YXIG5nhcqsoEbke2Gw1/+3hem0ici2AalU9E5+WEVG8bK1+H09VvefbfvCzH8L8STkutogSmeuBzSpD+F7bKgDr49AWIoqj+pPn8OCzf/Jt37JgIm7/8GQXW0SJzhOBTVUPADggIvODvS4i0wCIqtbHtWFENKjaO7vxpZ/vwLn2LgBA/piR+NYts7lgKA2IJwKbtQ6he23srREloUdfege73z8NwNxX+9fbFmBEBidh08B4JrDZyv25IjLVf7+IjAawUFVfdqNdRDQ4/vfdE72KG6++fhZmTxztYosoWXgmsFkbAKwO2LcGpjdHREni1Nk2fHnzLt/21ZePxRc+Ps3FFlEy8VRgU9XtABaKSLbf7sWqutWtNhFRbKkq7t+yGyeazUrYY0YNw7dL5yGF89UoRjwV2KwyACsAQETuAvCUu80holj64e/r8fLe477tR5fOw9isYS62iJKN5+7SquomEXkXwLcBrFTVhW63iYhi463DZ7Du+b2+7bs+MQ1XXz7WxRZRMvJijw0AtorIUwC2u90QIoqNts4u/P1TO9He1Q0AmD0xG/dfP9PlVlEy8lyPzVoPoAEASw8QJYnvbNuHd441AwAy01Pw3eULMCwt1eVWUTLyZI9NVZsA5LJ8FlFyeKO+AeV+qf1fveEKTB87ysUWUTLzZGADfDUkiSjBnW3rxJc374Sq2f7kZWNw+0emuNsoSmqeDWxElBy+9eu38V7DeQBAVmYaNpbMZcksGlQMbEQ0aF7Zexz/9foh3/Y3bpqNCaOHu9giGgoY2IhoUDSea8fqrbt92zfMGY+b5l/iYotoqGBgI6JB8fVn/4TjftVFvnnzHA5BUlwwsBFRzP1y12H8avcR3/b6JXO4GjbFDQMbEcXU8ebWXguHLiuahOIrx7nYIhpqGNiIKGZUFQ/+Yg+aWjoAABNzhuNrN17hcqtoqHE9sInIOhFZFMFxS0Tkvni0iYj659dvHsELe476tjfcOhdZmekutoiGItcDG8z6awsiOG45zHptRORBp8624cFn9/i2b7tqEj5x2RgXW0RDVdxrRdoVscsDdt8jIlf18bYcAMUAmgatYUQ0IA/+cg8azrUDACaMzsTaGzgESe5wowhyHoBSv20FMN0+wgkMiETkAf/95hH82i8Lct2SOcjmECS5JO6BTVUPiIgTxATAfpjFRTeGeWsD60cSeU/DuXZ83S8LsrTwUnxq5sUutoiGOleWrVHVA87XIlIJYJv/PiJKHA89twcnz5ohyHHZw/C1G690uUU01Lm+HpuqXud2G4iof17acxTP7jzs2374ljkYPZxDkOQu1wMbAIjIVAArAeSHOVRVdfmgN4iIwjrd0oGv/qJnCHLJgom49gpOxCb3uR7YRGQBgCpnM8zhCpP2T0Que/j5t3HCrxbkg5/lECR5g+uBDWZumgDYApNEUtf34UTktt/tP4mnqt7zbX/z5g8hZwRrQZI3eCGwFQGoVdWlbjeEiMJrae/Emqd7lqP5zOzxuH72BBdbRNSbFyqP5ACocbsRRBSZ77y0z7cidnZmGh666UMut4ioNy8Etu0ACtxuBBGFt/O9Jnz/dz0zc75245W4OCvTxRYRXcgLgW01gOki8hW3G0JEobV3dmP1lt3oVrP9iRljUFp4qbuNIgrCC/fYumGC20YRWQZgM/pIIFHVp+PVMCLq8fhvavHOsWYAwPD0VKxbwhWxyZu8ENhqYNL4BSaRpMhuBxK7PzV+TSMiANh3rBn//sq7vu37Pj0Tk/JGuNgiotCiCmwiskhVX45xG9YgeCCLOxGpAFDSxyErVbVXIWYRKQCwFmb1AQCoBLBOVZkQQ0mhq1uxeutudHSZ/6bzJ+Xg8x+b6m6jiPoQbY+tUkQaYKrsb1bVnQNtgKqGK34cT04SS6jlcRr8N0SkGMA2u+kMn5YAKBGRUlXdEvsmEsXXj/9Qjx2HzH+J9FTBxpK5SE3hECR5V7TJI9thlp1ZA6BaRPaJyFdsSayYEJFsEZkqItmxOmcU8gHUqWpuiIcvUIlIDnqC2mJVna6q0wEstvsq7DFECev9xhY88uI7vu17PzUDl4/LcrFFROFFFdhUdTGAXAD3AHgZwAwAjwCoFZHXReTO/gYkEblPRE4BaARQC2CF3f+iiNzZn3NGeX0nCEVa+WSFfd6oqpXOTvt1ecAxRAlHVfHVZ/6ElvYuAMBlF4/CvddEsmwikbuiTvdX1dOqWh4Q5HbCJH1sAtBog9EtkZ5TRF6EKa2VC+AAeteMvApAuYi8EG1bo1RknyO9N7bMPj8V5LWygGOIEs6zOw/jf/adAACIAOtvnYthaczdIu8b0Dw2vyBXCBOUNsIEpWIAW0SkS0R+LiLXhDqHiNwNM3y3A0C+qs4IOGQaTO9w8SD33HwrC4hIhYg0ioiKSLWIrAp1fLAkEb99nHhOCenU2TY89Nwe3/YdH52Kwim5LraIKHIDnqBt74ndBTP/7H5nN0zP6zSApTBJJ6+HOMVKmKzIRapaH/iiqjbZ3uFpmN7hYHHGWFbBJIDU2UcBgA0iUh1wPO+fUdL6xq/eQmNLBwDgktGZuO/TM11uEVHk+hXYbHLHfSLyBsw9sTKYXtcBmF7bdFWdoap5AK6DGaosFJGHg5yuAECNqp4Jc9kqhF+vbSCc3tUWVRVVLbTJINNhhicLRGRDwHtCZU+Gew0iskJEqkSk6sSJE/1vNVGMvfLOcfzCb/HQby2Zg1HDvDDllSgyUQU2EVknIu/CJHdsBFAIM4S4Bj3BbI2q+orJ2WSKa2F6cYuDnDbSZI3BDGpQ1cU2oJUG7K8D4OwLHJLsq9fWZ4/ODuEWqWrR2LFjo28w0SA419aJrz3Ts3joTfMvwTUzL3axRUTRi7bHthqmB7PDfp1r/zg/4h/MAqlqE0xvrjLIyztgekNTQr1fRKbBBLaqUMcMJhvc6mxbnAAbSUDus9dG5DWPvPgOPmgylftzR6TjwRu5eCglnmgD20r0DmanI32jnee1NshL62B6c9tEZJ7/WwBARObDBDSFyZx0mzNJ2wl0FySI+O3joqmUMGoONeJHf6j3bT/42Stx0ahhrrWHqL+iDWz5MMOPfRKRJSJyXyQntBmEa2DmxNXYuWwK4AH7dTVsxuUglPNy2ptjMyAbQ70O87032d4nAFTY52Ap/c6+siCvEXlOW2cXVm/ZDbXF7f7s8rG4ef5EdxtF1E/9GYpcEMFxyxFF78qW1ZoBk9Yv9pFrn7cDKAzR24sJG6wqAeSISK9gZIPadru52u89ziTsVf69Nvv1qoBjiDztsVdq8e7xswCAERmp+NbNs1m5nxJWn6lOIjIaPVU0HPeIyFV9vC0HZh5bVPeX7H0sX3KJiIyOZqgzBkphAtgKEVkKM/yZh97ZkoE/i1KYnlu1iDj3D4v9XiPyvHeONuOx3+z3bd/Pyv2U4MLl8Oah9x9ohUkeiaSuzoB6K3EOak6vrdBOxl6MnkoklQA2+JfN8nvPFhEpRO/q/lvA6v6UILq6Fav8KvcXTM7BX310qruNIhqgPgObqh4QESeICYD9MPeNwlXkb4h3YIoVOywa8YoDNoCxd0YJ6Qe/O4Bd75nBlYzUFGy4lZX7KfGFnXXpn8Zvh9u29ZXa3x8283EZws9VU1VdHstrEw1Vh0614NGX9vm2v7hoBi5j5X5KAlGVE1DV62LdABG5FaYcF9C7+HHQJsAkphDRAKgqHnjmTZzvMJX7Z43Pwj1Xs3I/JYdwySOL7JdVqnrGrrsWcQWQCNPz18IEtHKYRAzO/SIaZBXV7+PV/ScBACkCbLh1LjLSBlw6lsgTwvXYKmF6SYthUvFLYdL4NcLzR7LGRQGAalUdzALHRGQdP9OKb/7qLd/2nZ+YhnmTWNObkke4wPYITBBzelE1iCKxIkK+clVENLhUFV9/9k8409oJAJicNwJfXszK/ZRcwmVFrg7Y3o6eycqxsh1AiYhkqWpzjM9NRH5+tfsIXtxzzLe9bskcDM/g4qGUXLwwqL4KZumb7QG1Iokohk6ebcODz/ZU7r/tqsn4+IwxLraIaHBEmjzSL5Ekj6jqaTuNYAVMrcgwhysXhiLqh394dk+vxUMfuGGWyy0iGhyRJo/0V9gxDhFZD+Buu3kaPdXziShGnn/zCH795hHf9rpb5yIrM93FFhENnkiTRwZTiX0uHqzq/URDWcO5dnz9Fz1DkEuLLsXVl3NxW0peUSWPDJI8AJUMakSD4x9/uQenzrUDAMZlD8NX/5yLh1Jy80LyiCurYhMNBS/uOYpf7jrs2163ZA5GD+cQJCU3L1QeKQewWUTmqequSM9NRH1ramnHV5/pGYJcUjARi2aNc7FFRPHhhcoj+2GCW41d5LMGfSSQqOrTEV6baEh76Lm3cPJsGwBgbNYwPHgjhyBpaPBC5ZEaew0BcA9CB02xr3E2KVEYL/zpCJ7Z8YFv+1s3z0bOiAwXW0QUP16oPLIGg595STRknGhuwwN+Q5A3z78E131ovIstIoqvfk92FpFsmPttzj23OlXdGe157MKeRBQDqoq1T+9Gg82CHJ+diYdumu1yq4jiK+rAZgPaBphKIYGvKYDVqvpofxtkz58Hswr3mf6eh2goqqh+H5VvH/dtP1I6l1mQNOT0J92/BiaoCcyw5CP28bI930YReSHak4rIfSJyCqZuZK29BkTkRRG5sx/tJBpS3m9swT8917MczR0fnYJPXsaJ2DT0RNVjs+Wv8gFUw1QKOR3weg5MsFssIl+JtOcmIi8CKIYJlnXoPaXgKgDFIlKqqtdH016ioaK7W3FfxS6cbTPL0UwbMxJrPnOFy60icke0PbZimESPawODGgCoahOAa2EC1LJITigid8NMJ9gBIF9VZwQcMg2mN7iYPTei4H7w+3q8VmdmyaQI8OjSeVyOhoasaANbAYCavu592eBWA6AwwnOuhAmWi1S1Ptj5VHUxTIFkrrJNFGD/8WZsfGGvb/veT81AweRcF1tE5K5oA9sORFZ5JAcmuEUibLC0qiK8NtGQ0d7ZjS9v3oW2zm4AwJUTsvGlay9zuVVE7oo2sD0FIFdEbgl1gIjcChOAKiM8Z134QwAwqBFd4Dvb9mH3++auQEZqCv552XxkpHmhBCyRe6L6H2DnnFUC2CIiD9vakQAAEZkqIusAbAZQraprIzztDgAFIjIl1AEiMg0msLFgMpH16rsnUfbbWt/2qutnYub4LBdbROQNfQY2EekKfKAnOWQ1gFq//bUAVtnXckTkrgjbsM6+Z5uIzPPbr7YN82ECmsLMnyMa8k6dbcPfb94JtTV7/uzysfjCx6e52ygijwiX7l+P/pW7kkgPVNUaEVkDYD1MIeQme80HROQBmPt1AmAD12wjMtVF7t+yGyeaTYHjMaMy8GjpPKSkRPzfjiiphasVOT0ejVDVjSKyBUAZTDalAMgF0AQzL261qu6IR1uIvO5Hv6/Hy3t7qot8u3QexmYNc7FFRN7S71qRfbHDkAWqem+k71HVOpj5bM45RgebK0c0lL11+Awefr4ntf+uT0zDp2Ze7GKLiLynX4EtzIKjOTAV+6cBCBvY7P25bYFVRRjUiHo7396Fv/2vGrR3mdT+D12Sjfuvn+lyq4i8J9qSWqNhsiILwh2KyNP96wHEZciTKJH906/eQu2JcwCA4emp+NfbFmBYGquLEAWKdsLLWph7YAdgCh/X2/1OIeSnYYJahapeF+E5SwFMF5GvRNkWoiHjl7sO479eP+TbfugvPoTpY0e52CIi74p2KLIEgDr1HEWkHMB+AC85GYsisgLAOhHJUtXmCM7ZABPcykVkOUxP741QB6vq01G2mSih7TvWjNVbdvu2/3zuBJQWXepii4i8LdrAlg9gm7OhqnUiAgALYAoVQ1XLRWQVzJyzSJJH6mDS+wWmN1iI4FMMxO7n2AsNGWfbOnHPf1bjfOvVBqwAABpCSURBVEcXAFO1f/2SObD/74goiGgDW1OQfXUwS8v4q4FfhmMYa9C/uXJESU1VsWrLLtT53Vf73u2FyMrkwqFEfYk2sNUBKArYV4MLk0lyEGFtR1umi4gCPPnqATz/5lHf9rolc1gyiygC0SaPbIYpgvyCLXUF2Kr7TmFkW9exGJEXNyaiAK8faMC6/+6Zr/a5j0zBzQsmutgiosQRVY/NVghZBuA6mPXRlsFUC1kLUxjZWf064rqOtnBypLaxrBYlu+PNrfjiz2rQ1W1G6OdPysHXbuRq2ESRinqCtqoW2lWvG+z2aRG5FkAFzHy0JgDlqvpEhKdcHe6Sfl+fgk1SIUpGnV3d+Nuf7cBxWwcyd0Q6HvvLAs5XI4pCvyqPqOqmgO0a9H+SdV9JJvkAlsKsKLBKVb/dz2sQJYRv/vpt/PFAAwBABPjX2xbgkpzhLreKKLH0u1akiGTDBB4nSaROVXdGex5V3d7Hy9sBbLJz4x4XkWpVfSX61hJ5309eO4gf/r7et/33xZfjk5eNda9BRAkq6qV2RSRbRB4H0AigGmYIsgJAtV2bLeYVRFS1HGZBUq7HRknpf989gX/85R7f9g1zxuOL18xwsUVEias/PbYamALHTj3IGru/EGbIcKOILA4sahwDdQBujfE5iVy3//hZ3PvTnmSRORNH49HS+Vxfjaifoi2CvB5m6LEaQHFgBX4RyYEZPlwsIl9R1Udj1lIzVy7YBHGihNV4rh13/ugNNLd2AgDGZ2fiiTuKMDyDySJE/RVtj60YJkvxWlU9E/iiqjbZDMkGmKkAYQOb33y4UPIArITpJUa6YgCR57V3duOe/6zGwVMtAExlkSfuKMK47EyXW0aU2KINbAUAqoMFNYcNbjUwQ5ORqEH4klrOmEy4qQFECUFV8bVfvOnLgASAf142D7MnjnaxVUTJIdrAtgORlcrKQc+9t3A2IXxgqwWwRVUPRHhOIk97/H9qsbnqfd/2/Z+eietnT3CxRUTJI9rA9hSA9SJyi6o+E+wAEbkVJvhFlMGoqiujbANRQvv564ew8YV3fNtLCibi3k9xrV2iWIkq3d8WLK6EKZ/1sIhMdV4Tkam2PNZmmOHKtbFsKFEyeP7NI3jgmTd92x+elod1XIaGKKb6DGx2XlqvB0xKv8Dc76r1218LYJV9LUdE7oqmISKyxD+RRETmi8i79vyvi8i8aL85Ii/533dP4O9+vgM2qx+zJ2bjiTuKWC6LKMbCDUXWo39rpUX88VNERsNMH5gGEyx32n0vw9yrA8xSOdUiMkNV6/vRHiJX7TjUiJU/qUZHl/nvlD9mJH7411dxbTWiQdBnYFPVeAz8b4C5J1eJnnT+pTBBbYOqrhWREpghzlWIbFVuIs/Yd6wZn//BG2hpN6tgTxidiZ/c9WGMGTXM5ZYRJad+14r0JyLZfU0BCKMYQKOqXue3rxSmp7gOAFTVWRIn0lW5iTzhvYYWfO7JP+L0+Q4AQN7IDPzkzg9jIgsbEw2aqGtFAr77Xy+KyCl7f63Rfv2CiFwT5enyYRYr9VcEU1TZP1g6a70RJYT3Glpw26bXcOyMWYJmZEYqfvjXCzHj4lEut4woufWnCPI6mHtiiwHkAjhgH7kwC5BWishjUZyyDqa6iHP+BTDDkIFVRvLBklqUIOpOnMXSsj/g/cbzAICMtBRsuqMIcy/NCfNOIhqoqAKbnaO2Gmb17FJVTVHVGfaRAnNv7AyAlSJyS4Sn3QGgwK+ntxZmGLLC77rXwgS2umjaS+SGfceasbTsNRw53QrABLWy2wvxseljXG4Z0dAQbY9tJUzQKVDVrYEvquoWmFJaYo+NxDp7fKUd1rwVZhjyZQCwS+S8BL97bkRetefwaSwvfw0nz5rhx+HpqfjB5xfimlkXu9wyoqEj2sBWBKCmr5R7Va2DKae1MJIT2tW3iwDshAlwO2CGNB3X2f1rVPXpKNtLFDc732vCbeWvoeFcOwBg1LA0/OgLV+HjM9hTI4qnmGRFDpQNbqGKJhezRiR53Rv1DfjrH7yBs21m+ZnsTBPUFkzOdbllRENPtD22apj7YVNCHWAnVxfgwkzHfnGCmojcZYcliTzl+TeP4PYn/ugLarkj0vGzuz/CoEbkkmh7bGUwJbW2iUipqu7yf9GWxKpAQPJHJGzdyVDp/DkA1sBUJ/mb6JpMNDhUFY/9phaPvNhT0HjMqGH46V0fxszxWS62jGhoiyqw2YnSW2ESPGpEpAk9PbMimAAkACpU9YlIzml7eJUwvbw+DwUXGiWPaO/sxtqn38TWmp6lZ6aNGYnvf34hpo0Z6WLLiCjqe2yqWmpLXG2A6UH5VwOpA7A6WMZkH9bC3F+rA7AFpurIVACP2NenwwTSClVdFm17iWKt8Vw7Vv5nNV73WyT0w9PyUPa5QuSMyHCxZUQE9DN5xKb1bwEAEZlm9/U3waPEvF1n2POVA9gP4CW/lP8VANaJSJaqNvfzOkQDVnfiLL7wwzdQf6rFt6+k8FI8fMscZKT1q5APEcVYtBO0u0TkBf99qnpggFmLTgFk53zOJOwFfvvKATQiwsVLiQbDS3uO4pbHft8rqK26fiYeKZnLoEbkIdH22OphhgZjKViZrDoAVwXsqwGLIJML2jq7sO75vfjh7+t9+zLTU/DPS+fjM3MmuNcwIgoq2sBWCqBKRL6iqo/GqA11MIkn/mpwYTJJDlgEmeLswMlz+OLParDncE897ktGZ+Lx2wsxbxLrPhJ5UbSBrQEmuJWLyHKYIcQ3Qh0cYaWQzQDW2yHONaq6EybT8lYRuUVVn7H38YphVukmiotf7PgAX33mTZyz66gBwOIrx+GRkrlMEiHysGgDWx3MHDWByWQsRPAVtsXuD7vmvapuFJFlMKWzTgNYBjNfbi0AZx22fHs+3mOjQdfc2oGHnnsLW6p7UvkzUlPwwA2zcMfHpkIk4gXiicgF0Qa2NQgeyAZEVQtF5G6YHiFU9bSt6F8Bc0+vCUB5pHPjiPrrxT1H8Q/P7sHRM62+fdPGjMS/3bYAsyeOdrFlRBSpaCdobxyshqjqpoDtGsQ+UYUoqCOnz+Mfnt2Dl9461mv/LQsm4hs3z8aoYZ4oq0pEEfDc/1YRyYZZeLQhYAVtopjr6lb85A/1+PZL+3y1HgFgzKgMPPjZD+Gzcydw6JEowUQU2GywKULPYp9VsQ46InIfzH01J9VsNYBvi8iLADar6pOxvB7RjkON+Mfn3sKu93rPOFm+cBLWfGYWE0SIElTYwCYiSwBsQk/AAQC1RZCfiUUjbPAqhkk6cZJFHFcBKLbXuz4W16Oh7d1jzXjkxXcuGHacPnYk1i2Zi6um5bnUMiKKhT4Dm4gsgC2dBRNwnKCTD5OxWBBY4T9aNmlkMczctRJVrReRbr9DpsEkkSwWkTvZc6P+er+xBf9S+S6ernkf3X4pUBmpKfi/18zAPZ/Kx7C0sIm8RORx4Xpsa+1zmar6losRkTIAd8Ok3w+0F7USJtNyUbDhTVVtgglqDQDuAcDARlE53tyKx39Ti5++dgjtXd29Xrtx7gR85bqZrMhPlETCBbYCmALFvdZAU9WVtqe1MAZtKABQHcE9uyqEXmWb6AJ7Dp/G91+tx3O7Dl8Q0K6+fCzu//RMpvATJaFwgS0foat91MEMEw5UXfhDfG0h6lNXt+Llvcfx5Kt1eK2u4YLXCybnYNX1s/CR/ItcaB0RxUMkWZHBihT3tT9aOwAsEZEpqnow2AG2pFY+gG0xuiYlmeNnWvHLXYfxk9cO4qBf9X3Hgsk5+Jurp2PxleOYvk+U5Lwwj20dzEKi22zmo5OMogAgIvMBbAdLalGAs22dePFPR/GLnR/gd/tP9koIAYDUFMFnZo/HFz4xDQWTc91pJBHFneuBTVVrRGQNgPUAakSkCSaIPSAiD8BMMxAAG5yFR2noamnvxO/2n8Jzuw7jpbeOorWj+4JjsjPTcNuHJ+OOj07FJTnDXWglEbnJ9cAG+Aohb4EpflwIE8hyYYY7twNYrao7XGwiuaj+5Dm88s5xvPLOCbxWdwrtnRcGMxHgI9Muws0LLsFn512CERme+KdNRC6I5H9/gYh0hXqxj9dUVSP+62JXzvYtJCoio1X1dKTvp+Rx9HQrqg824o36BvzPvhM4cPJcyGNnjsvCLQUT8RfzLmHvjIgARBbY+nunfUB36BMpqIlIAcycv2K7qxLAOlvImfrQ3tmNfceaUXOoEVX1jag+2IgPms73+Z6Z47JwzayLcdP8S3DFhOw4tZSIEkWfgU1VU2J9QRF5fABvV1W9N2aNiQERKUZPtqYzdaEEQIlNhtkS/J1Di6ri6JlW7D3ajL1HmrH36BnsPdKM2hNn0RmY9RFgeHoqPj5jDK6ZNRafmnkxJrJnRkR9cONGhFNpBIi+V6cAPBPYRCQHPUFtsapW2v1OsKsQkVxbPSXpnWvrxLEzrfig6TwOnmrBoYYW1J88h0MN5uuW9pAj2r1kpqdg/qQcFE3Jw4fz83DVtDyWuiKiiLl1h91ZYXsbTB3IAy61Y6BW2OeNTlADAFWtFJFy+/oKAIO2jt1g6uzqxtm2TjS1dKChpR2N59rRcK4djS3taDjXgePNrTh+pg1Hz7Ti2OlWNPst+xKNS3OHY96lOSickouiqbm4YkI20lNjPlhAREOEG4FtKcwf+2KYZBFf7wZARYKtwbbMPj8V5LUymO9zGQYhsDXaANPZrejo6kZHl6LTee7uRnunebT5nrvQZrfPt3ehpb0L5zs6/b7uQnNrJ5pbO9Dc2omzbZ0R97AilZWZhpnjsjBrQhZmjc/GFROycPm4LGRlpsf0OkQ0tIlq3/c3BvXiIiUwQ5PX2l3+vTjPBzkRaQSQo6pBh1RFRAEg1OuOoqIiraqqiura//jLPfjh7+ujes9gy0hNwcXZwzBhdCYm543ElItG2MdITMkbgZwR6az6QUQxIyLVqloUuN/VyT42sWIL0CvIXWcf5SLi9SCXE/6QwZGRNvhDdSLAqGFpyBmRjrwRGcgdmeF7zh2RjotGDcP47EyMy87E+NGZyGXgIiIP8Mws1iBBbhlMqS0nyFXCmytp95UY0oQQwU9EnPtvmDx5ctQXvWhkBqZeNAJpqSlISxGkp6YgPVWQZp/TU1MwLC0Fw9JSMSwtBRnO1+kpGJ6eihEZqRieker7OjM9FVmZacjKTPc9j0hPRUoKAxURJRZXhyIjISK3wvTkigF0RzPpe7CFG2oczKFIIqKhLtRQpKdTz0RkCUyyiXMPzmvtjWTJnSGR6k9E5BVeCxQQkSUi8pQt1VUBoBRmaZuVMPUjvaQO8FUe6cVvX6TrzRERUQx4YlhPRBbBBK4SZxeAGpg0+nIPl9eqgBkiXQbTXn/OVICyuLaIiGiIc+0eWwIHs16c+2gACp3akLa3Vg2Ev78G8B4bEVF/eCbd39aKXIqeddacYLZFVROxAkkpTM+t2mZuAj3FkEvdaRIR0dDldq3IbTCBLQ/AigjmQKmqPjCIbYuaqm4RkUL0ru6/BazuT0TkCjdrRQKmpNbivg4MoAA8FdgAswo42DsjIvKEuN9js/PS+k1Vt8aqLV4hIicAHHS7HR4wBsBJtxtBvfB34k38vRhTVHVs4E7PT9CmoUNEqoLdCCb38HfiTfy99M1z89iIiIgGgoGNiIiSCgMbeUm52w2gC/B34k38vfSB99iIiCipsMdGRERJhYGNiIiSCgMbEVECEJESEakWkUYRUfv1Krfb5UUMbETUi4gUiEiF/QPaaL++YGkmih8bwCoAFABogClFWABgg4jUutk2L2JgI8/iJ9T4E5FimJUpSmD+gDbYr6tFpKSv99LgsL+TDTCLFheq6nRVLYRZn7ISQL6IcHksP8yKJE+yAWyD3ayD+U/tW7xVVae70rAkJiI5ABrt5mJVrbT7i2EKlgNArqpyVfg4EpEKmA8Xpaq6JeA13+8skiWyhgr22Mhz+AnVNSvs80YnqAGA/bo84BiKH+cDXWXgC/ZDRh3gC3IE9tjIg/gJ1R0iUg3zR7QwcMklv8Vza+yHDIoT5/5msGWw+P8hOAY28hx7MzwfIYa9wr1O/SMijQByQv2BdFaL5x9Q7xCRbTDrQG5U1dVut8crGNjIc/gJ1R3hAhcDm3fY4foymA94laoazbqWSY/32MhzVLWmj9XHK+zzxni1Z4jpqwfM3rHLRCTf9tK2wQS1jQxqF2KPjRICP6EOPvbYvC0gU3gLgNWqWudikzyLPTbyNH5CjatI/kiy1+YCmwW8AWZi9nRVLWVQCy3N7QZQ8ulHJYSmYJl2/IQad3UwUykKQmRFOsdQHInICphpFhypiBADG8VcLCZP20+oK2A+ofLTaXxUwGTYLYP5uftbZp85fzD+nGzHUldbkUB4j408x35CLQM/ocadcx8NfnPZ/Oaw8f5anAVUg+lzGFhVcwe/RYmBgY08h/PU3GPrQTqZp06li2L7fMGEeRpc/h8qwuGHjh4MbOQp/ITqPvvHdC16AlolgHV9TMEg8hQGNvIUfkIlooFiYCMioqTCeWxERJRUGNiIiCipMLAREVFSYWAjIqKkwsBGRERJhYGNPENESkRE7SNs6aZojydvswWvc9xuByU+BjbyqqURHLNy0FtBcSEi+QBqAWxyuy2U+BjYyKty7BpsfQn3OhENQQxs5EVOjcKQPTJb09D/WCIiAAxs5E019lHSxzFO0KsIdYCIFIhIhYjU2vtw1XblgFDHr7LHNNrja0VkQ7D7PkHOXSsiZYHH2n3qt56Z/2vF9rUNfvtW2X359hq1fhX3o/q+/M5VYO9HbrPfW61zTRHJsW30P1fQn3s0P8+AaxfYa6u9/jb/n4ddSNZZw8+5b7oh2HmjYdunwXr+fb1GiY+BjbzqKcD88Q/xejFMby3oOm12kdJqmODYBBMoCwCU2T+kgcdXwyxqmg+gyp47H8AqANsDjnXqWTrnrgSQB7N+XER1LiPgXzPT1yuN9vuyVqLnA0AV7PdlE26qYe5nOh8mCgBUBAbifl4XML+nanvNLQAanH32vhpgligqt1/Xwaw/1tc5I7XOPq/23ykiFbbtK1WVPf5kpKp88OGJB8wfTUVPgFEAFX0ctwLmj6QCKPN73XlvLYB8v/05MH9kFcAqv/0r7L5tQa5Va1/zP49zjuKAY7fZ/SV++8rsvoIg53bavsFv3yq7r9G/jf38vpxz9bq+33XVvi8nSHs39Pe6Qa4d+FpFH9e44Pc9wH9Tjf6/P792lcXyOjFsb6+fsd/PudH/98RH3w/22MiT1KyYHWo40hmG3Bzi7U7qf6+Vt9Ws7easQrzW7/gGmN5Er0/2lrP+WL7fvgJ7vsBP+6vtIxarfdep6saAfdF+X45y9VtyxrbbWRLobu295p3Ts/P/fvt7XQCo6eP7yA88eBD4em22978BZgHbXvdv7fCoq1NG7BBwPsy/R39LYYIa1yaMUJrbDSDqw1MACkSkRHsvcFkM88epSSToyjVFAJo0yPphqlonInUA8kUkR1Wb7LkvWEDTDpUFC6w1tl0V8FunzD7Has2yYENkUX1ffi8HGx5tgPljGXiuYEG5v9cF7JByBNcIyg6BXhTBoW9o8EVQy2GC2Qr7qNHgq7LXhGhrPC2G+UAT+DNcDCZJRYWBjbzM+aO0DDbw+CU2hEwagRm6QWDSRRB5sD0Xm/SxAsBCmB5ZX72JUphhxxKYZIcmmHtXFQA2x+iT9RtB9kX9fVmBPYBo9fe6wMB7r2ud64cR9MOJ/fBTDvO7bQJwbbA3q2qw3nq8OfeNg+0vD7KfQmBgI8+yf5QChyOX2eegw5B+WYlN6BmGCqXBvqcAJkEkB+YPcQ3McFkNzKflVQHtqgMw3Q5tlcL84XEeG0Tk2mC9m4Hoz/eVyNd16ABXSbft9yUgBfvQYX+PK1W1NPC1IMfWhjsmQJOqFkbYznwEJM3YUYOcwP3UNwY28rrA4cgS2GHIYAfbYNhkvw68txNKBcwfj9LA4SwRCTZs5VyrEvYTtv3DtBYmCFYAmB7BdSO+x9TP72vA3LpuDG2H+Tk7w6WBw9qA+fAS0e9CVSP5vfaHE3wDe2zF9rociowCk0fI65whmGURDkMCZlgwJ8TcsRxnLpXf7nyYexvB7tH0+oNn55fVBiYa2Ht1q4O9x8oLsi9sDyFAtN9XrLh1XUjPnMJwjwv+Tfin9aPvJJcCmO/RTaHurxXC3re1//bikXCT8NhjI08LGI50hsVCZUM6VsMkTGwXkUL/TD709M78/xA2wXyaz/c/1k4+doJpjm1Pnf3jskJEKvw/SfsFXv9hSGfoaiV6z0dzpipEI9rvK1bied3A+2nrEGHyiP+GTTopgckILbf7nKSfgoCh4iL0ZGq6pRgX3p8ETEak8++9JEF7zfHn9nwDPvhwHvCbxxaw35lndsFcMwSZx2b3b/B7Ty3MPYrGEOfwP3YbzB/pWnt8hd85VtjjVwU5d63fPv85Y/kBx1bAb+5X4Pfrt68kxM8omu8r5Lmc9gbZH3Q+WTTXjeDaF1wDJqA5568I9f1H+e8o8OdR0sd18/t7vRj8u/d970F+hrWw8wAD/1/wEfrBoUhKBP49tIh6BmqGBZ006TyYT+V1MEkCi4Mcu9K+XgwzNFUJYJqahAKnCkmhPX4jzNBWjT230/MqBzBde88Zq7Pvc9rh9OoW2/NEVWUjmu8rlgb7umqG4DbC9Fqc+VxRs8OlFbZtvYZ61Qw1N8Fksjo9Q+ceVizmHvaX8+9noy1vtsr2OMth/l0utNtu9yoThthPAkREQ46tSVmsEWQuDmIbymwbBisxZchhj42IhjIvJI6Emr9G/cTARkRDWRFiV7g6aqHmr9HAMCuSiIYsHeAE8Bgoss/sscUQe2xERO6pgkkkYoHjGGLyCBERJRX22IiIKKkwsBERUVJhYCMioqTCwEZEREmFgY2IiJIKAxsRESWV/w+ESmVHzpm4YwAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "x=np.linspace(-3,3)\n", | |
| "y=norm.cdf(x,0,1)*100 # convert fraction to percent\n", | |
| "plt.plot(x,y)\n", | |
| "plt.xlabel('Measurement=$x_i-\\mu$')\n", | |
| "plt.ylabel(r'Probability \\\\ Measurement $<x_i$ (\\%)')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## The problem\n", | |
| "\n", | |
| "If we take anything less than an infinite number of measurements, then the mean and standard deviation, $\\mu$ and $\\sigma$, are subject to random noise. Take the following five data sets of 20 normally distributed random numbers with $\\mu=10$ and $\\sigma=1$. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "data mean=10.207, standard deviation=0.945\n", | |
| "data mean=9.369, standard deviation=1.253\n", | |
| "data mean=10.015, standard deviation=1.038\n", | |
| "data mean=9.877, standard deviation=0.938\n", | |
| "data mean=10.186, standard deviation=1.194\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "for i in range(0,5):\n", | |
| " data=np.random.normal(10,1,20) # generate 20 random numbers with true mean of 10 and true std of 1\n", | |
| " # print the mean and standard deviation\n", | |
| " print('data mean=%1.3f, standard deviation=%1.3f'%(np.mean(data),np.std(data)))\n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Each average and standard deviation is subject to random noise (here the noise is artificially created by python). " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Problem #1\n", | |
| "\n", | |
| "Try creating 100 normally disributed random numbers with an average of 10 and standard deviation of 1. How close is the average to 10? How close is the standard deviation to 1?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "data # enter your work here\n", | |
| "pts=p.check_p01(data)\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## The solution - the Student t-test\n", | |
| "\n", | |
| "[William S. Gossett](https://en.wikipedia.org/wiki/William_Sealy_Gosset) was working for Guinness in 1908 when he found a solution to the small sample size problem. Guinness insisted his work be kept secret, so he published on the pseudonym, student. Now we refer to his solution as the \"Student t-test\". \n", | |
| "\n", | |
| "We can use the Student's t-distribution to:\n", | |
| "1. Establish confidence limits for the mean estimated from smaller sample sizes. This will give us error bars\n", | |
| "\n", | |
| "2. Test the statistical significance of the difference between means from two independent samples. \n", | |
| "\n", | |
| "### 1 - Confidence intervals\n", | |
| "\n", | |
| "The confidence interval for a given data set is calculated with the mean, $\\mu$, standard deviation, $\\sigma$, number of data points $N$, and t-statistic $t=f(\\sigma, N)$. The equation for the confidence interval of a small data set is as such\n", | |
| "\n", | |
| "x=$\\mu \\pm t(C.I.,df)\\frac{\\sigma}{\\sqrt{N}}$,\n", | |
| "\n", | |
| "where x is the measured quantity, C.I. is the chosen confidence interval (95% here), df=N-1 is the degrees of freedom, and $t\\sigma/\\sqrt{N}$ is an estimated true standard deviation. For example, lets look at the previous data set of $\\mu=10$ and $\\sigma=1$. Notice the confidence interval increases. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "data confidence interval is 9.687 +/- 0.427\n", | |
| "data confidence interval is 9.986 +/- 0.417\n", | |
| "data confidence interval is 9.796 +/- 0.579\n", | |
| "data confidence interval is 9.995 +/- 0.475\n", | |
| "data confidence interval is 10.086 +/- 0.545\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "N=20\n", | |
| "tstat=t.interval(0.95, N-1) # calculate the 95% interval t-statistic\n", | |
| "for i in range(0,5):\n", | |
| " data=np.random.normal(10,1,N) # generate N random numbers with true mean of 10 and true std of 1\n", | |
| " # print the mean and standard deviation\n", | |
| " mu=np.mean(data)\n", | |
| " sigma=np.std(data)\n", | |
| " N=len(data) \n", | |
| " print('data confidence interval is %1.3f +/- %1.3f'%(mu,tstat[1]*sigma/np.sqrt(N)))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEbCAYAAABUV7o5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9dX48c+ZmbCJEHABBRECCFoXCBG3Witbtfqrj8qibW31UcC2VvtUC2j71OVpVdDW1lor4FZ3CS51V8Cl1J0E3FlMQBQV2QLKlszM+f3xvTMZwiSZSSa5M5nzfr3mFe6dO3dOCOTM997zPV9RVYwxxphcFPA7AGOMMaapLIkZY4zJWZbEjDHG5CxLYsYYY3KWJTFjjDE5y5KYMcaYnBXyO4B8svfee2vfvn39DsMYY3JKWVnZelXdJ9lzlsRaUd++fVm0aJHfYRhjTE4RkU/qe84uJxpjjMlZlsSMMcbkLEtixhhjcpYlMWOMMTnLkpgxxpicZUnMGGNMzrIkZowxJmdZEjPGGJOzLIkZY4zJWZbEjDHG5CxLYsYYY3KWJTFjjDE5y5KYMcaYnGVJzBhjTM6yJGaMMSZnWRIzxhiTsyyJGWOMyVmWxIwxxuQsS2LGGGNyliUxY4wxOcuSmDHGmJxlScwYY0zOsiRmjDEmZ4X8DqAuERkBFAH9gUKgAqgCFqnqEj9jM8YYk12yIomJyAXAZKA4tivJYSoiVcA8YKaqvtRa8RljjMlOviUxEekCXA5MwSWtSuARXJLa6G0DdMeNyI7EJbnxwDgRqQSmqOpjrRy6McaYLOFLEhORM4Dbvc0bcCOrlY287JGE14/CjdweEZEXgPGquqVFgjXGGJO1/CrsuB2YqqrdVXVaCglsF6o6X1XH4UZpq3AjOmOMMXnGl5GYqnbP0HmqgAszcS5jjDG5J6tL7EWki4hcJyLPi8i13n00Y4wxBsjiJCYiXYGVwFRgL2AaUGmJzBhjTEzWJjFgNq4qsb+qlgADcPfAZvkalTHGmKyRzUlsFDArVvShqpW4xDba16iMMcZkDV+SmIiMEJG+KRzarc62Zj4aY4wxucqvkVh/3P2tWxu4x1WKm9R8BICIDAUmAXNaKcZGicgUESnzOw5jjMlXviQxVZ0NjPEem0Tk2iSHTcHNASsXkRXAIlwXj6npvp+IFItIqYhs8h6lIlLc+CsbPOcoYDq1rbKMMca0Mt/uiXkTlgcAPwMuFJENInJ+wvObVbU/biLzYmCaqg5ItzOHl2zKgLG4dlYbvT+XicjYpsQuIoW4kaIxxhgf+V7YoaqzvMnP04HZIrJCRE5MeH6Gqo5X1RvSPbeXbOZ5m6NVtb+XGGPFIaXeMekqxVVOGmOM8ZHvSSxGVWfgSuhfBBaIyHMicmAzTzvJ+zpDVecnvNd8akv1J+32qgaIyBRc5eRkapsUG2OM8UHWJDFwbaRUdTJuTliQxos/GjPB+/pwkudm1jmmUd59tOnAfFW1+WrGGOMzX9cT8xbAnIlbBLMCmKyqL3lzwkZ797OuxxV/TFfVK9J8iyIAVS2v+4SqlosIpFeYUYpboHNcmnEYkzV21ER4edk61lRtr/eYor334NsD96YgmFWfc43ZjZ/riY0EXsAr2gDOAuaLyKjYgpfeZb8SrwBjlohMxq0hdkeKb5Ox+1YiEku2o73Gw6m+bhLeJcs+ffpkKhxj0qKq/PLKG1kR2ZuVke5Up/BfvwM19A9uYEBwPX+9Ju2iYGNahaj6M39YRBYBXVV1YMK+CmCDqg6v5zVTgOtUNZjieyhQpap1J03Hnt8EFKpqspWkE48bixuFzVDVqQn7K4Cixl4fU1JSoosWLUrlUGMy4quvd/D44jXMLfuM5Wu/afJ5Du/dlbHDevODI/ansFO7DEZoTONEpMxrP7j7cz4msShuMcyfJey7DZjYUJISka6qujnF91CA+pJMY897xxTiGhFv9CobE5+zJGay0sat1Vz5xAc8894XRKK7/x/v070TJxy0D6Hg7v90d4ajLPhoLWu37NztuXbBAGcO681vTzmYzu19vRth8khDSczPf4WVePesEhTRSMVfqgmsgfeoq7FLg+NxlyU3isi8Os8VASTsn55YBWmMH96o3MAlDy3eLQlFq7ezbdl/+Oa9BVSuepdAoP7PXpHTDqXzgBI6HzaKTgOPQUIFAFRHojz41mreqNzA384eyqG9urbo92JMY/xMYrOA671uHXNwVYKjcJ06MqUSKBKR4rrFHQkdO1Itky+i/oQ4yvs6s57njWlxkajytxdXcPOCFSQOvob3687YYb0569iD0JodAA0mMIBgQNixspwdK8sJdOjM3S++R2nZZ7z7mfsMuXL9Vs649TWu+P5gfnpsX7wiKWNanW+XEyFeLDER19hXcF3rM7ZSs1dUMZM697K856bjEubkppbL2+VEky3WbtnBJQ8t5o3KjfF93ToV8KfxRzBicA+AXRJNKv/vkx3/2OLP+O1j77OtOhJ/bswhPZgx9nC7V2ZaTFbeE4sHINIPV+ZeHlt2JcPnj32Dw2KjMW8UVga73w/zyvpJ5bKgJTGTDV5a9hWXznmHjVur4/uO6tedv541lJ5dO8T3ZSKJAVSu+4aLHljMh1/UdoDrVdiRm88ewrADuzf5+zCmPg0lMb+WYolPXlbVlar6SHMSWCOToWNzuspEZJ53/6qsznOx88TaVM0TkcbupRnjK1XlhueXct5db8cTWEDgV6MG8sDEo3dJYJlUtE9nHvvFsZx7bN/4vjVV2xk/8w3u+E/GP4ca0yC/ZjJWici1zejEAcTXJXsb1yQ4KVWdCwwD5gIl3mMubmQ2tznvb4yfbnxhGX9/qSK+3aNLex6YeDS/GnUQwUbueTVX+1CQq37wLWadM4yuHV3RRySq/N9TH3L3q5bITOvxK4kNoHYZlocSG/42RkT6ishl3vIs83AtoOpNYuC6c6jqOFXt5j3G1dPFo0pVxXs0WvDhNRS2O9qm1d32SsUuCeyEg/bhmYuP5+iivVo1jjHf6skzlxxPcZ/avgJXPfkhj5R91qpxmPzlS3WilyBKvMKLKcB4795VOW7dsApqS98Lgb1wlYGjvG3BjabGtMR9NGOy2f1vfsL1zy6Nb48cvC+3nTPMtxZRvQo7ct8FR/Hj29+kfLX7b/ubue+wR/sQJx3a05eYTP7wvbAD4sUU44CRNDyvaz5u9DUrzfliWcEKO0xz/WvJGn718BJi/22PLurO3ecNp0NB401sMlXYUZ/N22s4a9YbfOQVfLQLBrjj3BKOH7hPo681piFZXZ1Yl4h0pXZOVnfcPK6NqrrY18AywJKYaY4FH61l8r1lhL1JYEf07sr9E49OuXNGSycxgHVf72T8zNdZuX4rAB0Lgtx3wVEMOzBp5zdjUtLk6kQR+YeIXFBne0SmA0zkrei82KtYnK2qC9pCAjOmOV6v2MDP7y+PJ7CB+3bm7vOGZ13rp332bM99FxzF/l5l5PaaCOfd9RYffp7WguzGpKyxi+iTgbF1toe2XDjGmLre+bSKC/75NjvDUQAO6O7uQXXbIzsnF/cq7Mi9FxzFXl58W3aE+cmdb1K5rukNiI2pT4OXE0WkDBiCK7jYCIzGFV2k0qpJVfWkTATZVtjlRJOudV/v5Ps3L2Td164P4r57tmfuhcfSZ69OaZ+rNS4nJnp/zWbOnv0GX+8IA26Nsid/+W32yLLRo8l+zWkAfAHenCpvW4H+3qMx2XWzzZgcE40qv56zJJ7ACjsVcN8FR6WWwK5qpDFvY8835/irXM3Vob26cte5R/LjO95kR02UyvVb+d3j7/Pn8UdYr0WTMQ1eTvTuTfVX1QCuyEKA6UC3FB7Wf8aYZvjHKxUsXLEeABH429lDOajHnj5HlZ6Svt259vTD4tuPeWubGZMpKU8s8VYzngXM84ovGn20XNjGtG1vr9rIn+ctj2//7IT+OVuqfkZxb8YO6x3f/v2/PuDjr772MSLTlqR1cbq+DvMi0kVVrfzImAzYtLWaix9cHF/MsuTAbvx69EE+R5WGJJcer9H2LJY/UKG92F4T4Rc33ce/2v0vHaSmBd7fPj/nkyZN8ReRISLyvIhsEJEIrn3UBhF5Lp0WUsaYXakql5W+wxeb3bpfhZ0KuPnsoYR86saRKZ1kJ38vuJn2uEbFy7QPV4d/4nNUpi1I+3+GiFyP6wI/Gnfva6X36IbrhzhfRG5N85zXpTL/TETOEJHL0o3ZmFxxx39WsmDpV/HtG8cewf6FHX2MKHMGBz7lqtA/49sPRkbyROQYHyMybUFaSUxEzsT1OtwMjFPVgKoO8B4BYDywBZgsIqenceqppDb/7CxcYYkxbc47n1Yx/bnanojnf7sfow7p4WNEmXdW8CX+X+C1+PYVNeezKtq2vkfTutKdsDEZVzpfrKqr6j6pqnNFpBz42Dv2sWQn8VpL1V1N+UIRGd7AexfiGgBXNXCMMTlpy44aLnqwnJqIuw92eO+uTD1psM9RZZ4IXFtwB+9WF/GJ9uQbOnFRzcU80u5K2kvY7/BMDko3iZXgVmBeVd8BqlrpJbIjGzhPd3ZdkDKd+Wd1k58xOe/3j7/Ppxu3A7Bn+xC3nF1Mu1Bu3werz56ynVsK/sYZ1VdTQ4j3tR9/Co/nioIH/A7N5CBf/pd4y6fEktYA3PyzmQn76nt0a2ztMGNyzUvLvuLxJZ/Ht68/8/AmdeTIJYcFVnJF6P749u2R7/NetJ+PEZlclW4SKwOKReTA+g7wLhUW49YFq5eqrvQelXhLrCTsq+9htbOmTdlWHeZ3j70f3/6vIftzyuH7+RhR6zk3+DzHBtz3HiXAtJoLCGvbHH2alpPuv5iZuFHTPBE5ou6TIjIEl7wUKE31pKo6RlUfTTMWY3Len19YzpoqdxmxW6cC/vfUQ3yOqPWIwLWhO+Jl9x9oP+6MnOxzVCbXpDvZea6IPAKcCZSLSBW1I64SalddLlXV29M5t4j0xRWDNLQopheGnpXOuY3JRu99tpk7X61dmPx3pxzCXp3b+xhR6+sbWMsloUeYET4bgD+Hx3Jy4C0OCKzzOTKTK9JuJ62q40RkLK7UvR9uvlhMJTBVVR9J55wiMpTaZNhYZ1DFldobk7PCkSjTHn0XrykH3x6wN2cU9/I3KJ9MDD7DE5FjWaoHsoP2XBE+n3sKrsd6BJtUNOkCtKrOTWgM3B/onzBnLK0E5pmOS16P4JJiQ8UdA5oSszHZ5M5XV/KBt1Bk+1CAP55+aN52di+QCNcX3I7g1ktbGD2cf0WP8zkqkyuavbCPV2nYXCVAhaqOz8C5jMkO9Sxf8ml0H/5cPQNwlw5/xf0c+Lf8vrgwJFDBT4MvcHfELUF4Tc05nBB4h25iC2mahmVLKVAhbuFNY9o0VbgifD47vAR2sKziguAzPkeVHS4LzWF/3NIzG+nCH2p+5HNEJhdkSxJbgCvLN6ZN+1f0OBZGDwdAiHJ9we0USMTnqLJDZ9nBHwrujG8/Ej2BVyPf8jEikwuyJYlNBfqLyKV+B2JMS9moe3JNzTnx7XODz3NEoNLHiLLPiOASTgm8Ht++Inw+O7TAx4hMtmv2PbEMieIS2QwRmQDMwVU6JmVzykwumhGewEa6ANCLdVwWmuNzRNnpyoJ7WLjzcLawB59oT24Nn8avC+b6HZbJUtmSxMpxpfOCK/Io8bbrEm9/sPVCM6b53o8eyMOR78a3rym4mz1kp2/xZLN9ZTNXhO5nWngSADMjpzI+9DK9Zb3PkZlslFYSE5ERqvpiC8QxjeRJy5icpwrX1PwE9a7enxhYzMjgYp+jym7jg69wf2QU72kRO2nH9TVnc0u7v/kdlslC6Y7E5ovIRlwn+TmquiQTQajqjEycx5hs9Ez0KN7SgwEIEeZ3oft8jij7BUT5fcG9jKu+EoCnosfwk+gLDA8s8zkyk23SLexYgFtGZRpQJiLLReRSr2VUxohIFxHpKyJdMnleY1rbDi3g2pofxrd/EnyB/oEvfIwodxwZWMapCUUe19T8hKjm54RwU7+0kpiqjga6ARcCL+K6Z9wAVIjIWyJyfnMSj4hcJiIbgE1ABTDJ2/+8iJzf1PMa45fbI99nDfsA0I2vuSSUdJ1YU4/LCx6INwh+X/sxN/IdnyMy2SbtEntV3ayqs+oktCW4YozZwCYv6ZyeznlF5Hlc+6luwEp27aE4HJglIs+lG68xflmrhdwaPi2+/etQKV1lq48R5Z5esoHJwafi2zPCE/haO/oYkck2zZonlpDQhuGSzwxc8hkFzBWRiIg8JCInNnQeEZmI65m4GChS1br9EfvhRn6jbURmcsX0mrPYRgcABslqzg62RE1U23dh6El6sgGA9RTy94QPBsY0e7Kzd//qAtzcrt/EduNGU5uB8biCkLcaOM1kXHXiCFVdVfdJVa3yRn6bcSM/Y7Lakk+reDRae+nr96F7CUnUx4hyVyfZydSCh+Lbd0ZO5pPovj5GZLJJk5KYV3RxmYi8jbt/NRM3klqJG4319zradwfG4C43DhORa+s5ZTFQrqpbGnnrRTS+3pgxvlJVrnnyg/j26MAijgt+0MArTGNOC7zGEFkBQDUFXBv+YSOvMPkirSQmIteJyApc0cUMYBjuEuA0ahPXtMTO9qo6HxiJG52NTnJaaKA7Rx2WwEzWe+KdzylfXQVAAWF+G7rf54hyX0CUKwvuiW8/Hx3Oa5H8WQXb1C/dkdhU3Jpei70/d1PVElW9oaElWVS1CjdKm1/PIYuBYhE5sL5ziEg/XBJbVN8xxvhtW3WY659dGt/+7+Cz9A2s9TGitmNooILTAwvj29eEzyFiJfd5L90kNpldE9fmVF/oLaJ5eT1PX4cbqc0TkSMSXwYgIkNwyUtxFYzGZKXbF67ki807ANibzVwUetzniNqWqQUP0RH397tUD6Q0coLPERm/pZvEinCXEBskImeIyGWpnlRVy3GXJAcA5d5cMQWu8P5chlf92EJtr4xptvXf7GTmKxXx7UtDc9hTtvsYUdvTUzZxYejJ+PZN4bFs13Y+RmT81pTLiUNTOO4s0hwxea2nBuBK6cV7dPO+LgCGNTCSM8Z3t7z4MVur3dpgA/btzLjgKz5H1DZNDD7DPmwCYC3dudNbDdrkpwZ7J4pIV1yfxEQXisjwBl5WiJsnVpVuMKpaSULxh4h0TeeSpTF++WTDVu5/85P49tSTBhOaYyX1LaGT7ORXoUf4bfgCAG4L/4Czgy/RXb72OTLjh8YaAHcHxiVsK66wo38K566b/NJmCczkihueX0ZNxC3EcGTfbow62OYxtaQJwZe5I/J9KnV/vqYTt4RP4/cF1lg5HzV4OdGrOIwlrQG4S3szE/bV9+hml/5Mvnjn0yqeere2qe+0kw9GxKrmWlJIokwJ1U6Avjcyhk+j+/gYkfFLo0uxJJbOi8h8YF5D5fRN5VUgTqDxuWCqqmdl+v2NaQpV3aWk/qRv9WTYgd18jCh/fC+wiGGyjDIdRA0hbgyP56/t/u53WKaVpbWemKqOaYkgRORMXNsq2LXxb9IwcIUjxvju5eXreL3S9fULBoTfnDTI54jyhwhcXvAgY6uvAuBf0eOYGH2aQ/0Ny7Syxgo7Rnh/XKSqW7x1w1LumpFGOfzluOQ1Cygl9Q4exvgmElWmJ4zCzjryAPrv09nHiPJPSWA5owOLmBctAeD68NnYnbH80thIbD5u5DMaV/o+Dlc6rymeP5jiccVAmapac1+TMx5fvIalX7qKuI4FQS4ZOdDniPLTlNDDLKguJkqA/0QPY+GKdRw/0O6P5YvGktgNuIQVGxmV43omZlolNvoyOWRHTYQ/z1se3554fD/27dLBx4jy18DAGsYHX+ahiLtwdP2zSzmu/94EAlZckw8aTGKqOrXO9gLcxONMWwCMFZE9VdUme5isd8/rq1hT5bpx7LVHOyadkMqsE9NS/ic0l8cjx7GD9nzw+RaeeOdz/mtoL7/DMq2g2euJZcgU3JIuC+r0TjQm62zeVsPfX6ptL3XxyIF0bp9WjZTJsB5SxfnBZ+PbNzy/jJ3hiI8RmdaSamFHk6Ra2KGqm73y/Um43omNHK72G8P45rZ/V7B5ew0AB+7VibOH9/E5IgMwOfQkD0RGsok9WVO1nQfeXM15x/XzOyzTwlIt7GiqlAo7ROR6YKK3uRnY2Iz3NKbFrN2yg7terZ0meemYQbQLZcsFjfzWRbZzUegx/i/8E8D1shxXcoCNktu4VAs7WtpY7+so61JvstnNC1awo8b1RPzW/l049bD9fI7IJPpRcAF3dp7MmqrtbNhazR0LV3LJKKsabcvSKuxoQd2B+ZbATDZbuX4rD739aXx7ykmDrQIuy3SQGn41aiC/mfsuALMXVvLjo/uwV+f2PkdmWkq2XAex1ZpN1vvzvOVEou7CxNFF3fnOwL19jsgkc0Zxbwbu6yadf7MzzK0vVzTyCpPLsqVjxyxgjogcoarvpHp+Y1rL+2s28+Q7n8e3p5w02Jr8ZqlgQLjse4OYfG8ZAPe+/gn//e1+9Crs6HNkpiVkS8eOj3GJrFxEZuImVddb3KGqj6Z4XpPvruqakdPcUD0FGALAmMDbFN/5w4yc17SMMYf0YMgBhSz5tIrqSJS/zFvODeNs9k5blC0dO8q99xHgQupPkuI9l2pyNKbZXo8czCtRl8ACRLksNKeRVxi/iQhTTxrM2bPfAOCR8s+Y9J0iBvbY0+fITKZlS8eOabROFaQxaVGFGeHaRRPOCC7koMAaHyMyqTqm/15856B9+PfydUQVbnxhGTPPKfE7LJNhTZ5AISJdcPfHYvfIKlV1SVPOpaotMbozptnmRYexWF2Jdjtq+FXoEZ8jMumY8r1B/Hv5OgCe/2Ati1dvYmgfW++tLUm7OlFEuojIP3BtospwS6eUAmUiEhaRS5sblPcefb1EaYwvIircEJ4Q3/5xcB69Zb2PEZl0HdqrK6ceXjuXb/pzS1G1iz5tSVNK7Mtx7aEEd2nxBu/xone+GSLyXFOCEZHLRGQDLkFWeO+DiDwvIuc35ZzGNNVj0W+zQnsDsAfb+UXoXz5HZJri0jGDCHrz+d6o3MjCFfZBpC1JK4l57aGKcImsm6qOUdVp3mM0btLyEmB0uiMyEXkeV/nYDVjJris8DwdmNTU5GpOuHVrATTVj49sTQ0+zl9gCC7mo3957MOHIA+Lb059bSjRqo7G2It2R2ChcAcZIVd1c90lVrQJG4hLQhLrP10dEJuLK+BcDRao6oM4h/XAjvdE2IjOt4b7IKNbgFlbszhYuCD7jc0SmOS4ZOZAOBe7X3Qefb+Gp977wOSKTKekmsWKgXFW31HeAl8jKgWFpnHcyLjmOUNVVyc7pjfQ240rwjWkxW7Qjfw//V3z7l6HH6Cw7fIzINFePLh126Wj/pxeWUR2O+hiRyZR0k9hiUuvYUYhLZKlqNDl6FqX4/sY02ezwqWzCzSfqLV/xw2BLzCoxre3CE/rTtWMBAJ9s2MbDb6/2OSKTCekmsYeBbiJyen0HiMiZuEQzP43zVjZ+CGAJzLSwr7SQ2yMnx7cvDZXSXsI+RmQypWvHAn7+3doVuP+6YAVbd9rPNtellcS8+Vzzgbkicq3XSxEAryT+OmAOUKaql6dx6sVAsYgcWN8BItIPl8SsWbBpMTeHT2c7HQA4WFZxWuA1nyMymfTTY/uyX1f3813/TTV3/GdlI68w2a7BJCYikboPags3pgIVCfsrgCnec4UickEacVznvW6eiCQ2OFMvjiG45KW4CkZjMm5VtAcPRU6Mb08JPUxArIqtLelQEOR/Rh0U357170o2fLPTx4hMczU2EluFK3ev+6isZ3/subTae6tqOa711ABcE+ANuIR1hffnMlzp/Qxbc8y0lBvD4wl7TWyOkg/5bsAWVGiLzijuxYCEpVr+/pIt1ZLLGuud2L+h5zNJVWeIyFxgJq6yUXCJqwo3qXqqqi5urXhMfnkv2o+nosfEt6cVPIittNI2hYIBpnxvEJO8pVrue+MTzjuuLwd07+RzZKYpWmRRTBG5QERuTfd1qlqpqqNVtbuqBnATqrt7k6qbnMBEpFhESkVkk/coFZHiJpxnrIiUeedQ789TmhqXyR7TE5r8nhR4i6EB+3Telo0+pAfDDnQ9FKsjUW6at9zniExTNSmJeUUcI+p5nIG7NDg5jfNFknXjSDahugmxjsJdjhyLW6Nso/fnMhEZ29Br65xnCq5HZLF3jnLvz9NFxH7j5bCFkUP5T/QwwJZayRexpVpiHluyho++aGyGj8lG6bad6ioib+OKOObV8yjFVRGmM7lmFZDxS5ciUujFBDBaVft7l0hHe/tKvWMaO88oXEFJFTDMO88w3OXO+UCRt5inyTFRFaaHz45vjw++zIDA5w28wrQVw/t1Z8TgfQFvyZ3nlvockWmKdEdil+PuV63ENf1d5e2PNQF+FHcvq1RVx6Rx3nFA/0x0wK9jkvd1hqrG5615f55V55iGxEaVE70ilNh5qnCxp3oek2Wejh7F++o6ObSn2pZayTNTThoUv/f50rJ1vFm5wd+ATNrSTWJjAVXVAao6DTeiEeAFrwnwOFxbqFEiks4SqhtxyeAKEXlbRK4TkTPqe6Rx3lj/xoeTPDezzjENid0/220Ct5fIKiE+8jM5olqD3BgeH98+L/gcPWWTjxGZ1ja4ZxdOH9orvn3ds7ZUS65Jd1HMImovz6GqleI+xgzFNehFVWd594+mAz9P8byVuJJ6wY30hpF8pWfx9gfTiJfE0VNC7OVe7KkUeIzzXlO1W0AucRXV97zJXvdFRvOJ9gSgK9/ws9CTPkdk/PDr0Qfx1DtfUB2JsuTTKp5+7wtOPXx/v8MyKUo3iSX7JV2JWyolUTm1951SMY3kSau5MjIySpYEE5R6X5OuTi0ik/AuNfbp0ycT4ZgM2Kx7cHO4tnvaL0OP0VW2+hiR8Uvvbp0477i+zPy36343/bmljD6kB+1DqX5WNn5KN4lVAiV19sWq9BLFRyep8NpZtZSGRkdVNDHRecUeM/H6RKrq1GTHqeosvPtvJSUldp0iS9wSPo0qr8lvH1nLOcF5jbzC5KhLJZsAACAASURBVIyruqb9kp9rJx7mJqrYk083bueeq37MxFATlt+5qtkF1SZN6d4Tm4NrAPyc1woKvM7ysabAXo/DUaTe1LelNZSk0k5gIlIkIrFKzCJc0Ug6o07js9XRffln5Hvx7amhh6zJb57rKtu4JPRofPtv4dPZpJ19jMikKq2RmNdVYwIwBre21wTcaORyXFPgStwv9rR6HHqNg1M1L43WU7F4GpLyfayEe30Ac3FdRLIlWZsUTQ9PoBq3JEexLOf7gTd9jshkgx8F53NPZAwrdT+24C43X1lwr99hmUakezkRVR3mrcS80dveLCIjcfeG+uOSwixVvT2N0ya9FJf4tgl/3oBXRJKCStwosbjufa2Ejh0pJSFvHtgk3OXTcZa8clNZdCBPJ7SX+m3B/dZeygDQTiJMDT3IhTW/BuDeyGh+EpxHv8CXPkdmGpJ2EgNQ1dl1tstp3mTlhi7HFQHjcd3zp6jqjWmctxR3aXMCuy/SGSutb3SSckJxxny7dJi7VOGPNT+Kb58SeJ1hgRU+RmSyzfcCixguH/GWHkyYENPDZ3Fbu7/4HZZpQJOSGICIdMElmNjlukpVXdKUc6lqQ909FgCzvUTyDxEpU9WXUjzvLG8ENUVEHo6NxrxR2JTYMYmv8Qo2YhOiY2IjxXGYnPVsdDjl6pbhKCDM1FCy6YMmn4nAFQUP8F/V/wfAc9HhvB0dxJGBZT5HZuqTdu9EEekiIv8ANuF6EpZ6jzKvB2Kmu24A8WSzmPTXE4slnjIRmecVZZTVeQ7YpU3VPBEpStgXS9QrE5oI7/ZoyvdlWke1Brk+ob3UT4PP0yfwlY8RmWw1JFDBDwKvxrf/UPMjbP5z9mpKA+By3KU1wY2SYi2nXvT2zUjWzDdDKnEToVOmqnO918zFTQ8o8f48zHuuMYmFIYWNPEyWuicyhtXaA3ATm38ZetzniEw2+03oYdpRDcA7OoAnE+6jmuyS1uVEEbke90u9DBhVt8u8N2pZAIwWkUtV9U8Zi9QpJo1qwhjvMmKjlwK9jhtSZ1953X0mt1TpHvwtYWLzxaFHbWKzadABgfWcF3yOmZEfADC95izGBBbRQWp8jszUle5IbBSuUnBksmVSvCQwEvdLP5WehACIyJBGHiNE5GGgH25emjEpuzl8Optxc34OlC9tYrNJyS9C/6IbXwOwhn24O2Fuocke6RZ2FANlqlrvwjuqWiUi5aR32a+cxttOxUZDjZXjGxO3ItqLeyK1CypMCz1IO4n4GJHJFV1kO78KPcKV4XMBNwH69OB/6CHWIjWbpJvEFpNaO6lCdi9pb8hsGk9iFcBcVV2ZxnlNHlNVrgr/lLD3z/wo+ZCTAm/7HJXJJT8MLuC+yChWaG+20pHra87mpnb/8DsskyDdJPYwcL2InK6qjyU7QETOxCW6lKsIVTXlVaCNSdWz73/Jq9FDAQgS4eqCf9rEZpOWAolwdehufljzOwAeix7PD6MvWsl9FknrnpjXqHc+rsXUtSLSN/aciPT12kfNwV1yvDyTgRqTju3VEf7w1Ifx7XOC8xgc+NTHiEyuOjb4IacE3ohv/77mp0TUPg1liwaTmDfva5cHtYUbU4GKhP0VuAnEAhSKyAXpBuMtejkkYXuIiKzw3uMtETki3XOa/HTryx/z+eYdAOzFZv4nlMpsCmOSu6LgfjqwE4CPtC8PREb6HJGJaexy4iqats5XWh9TRKQrrmy/Hy45LvH2vUjt/KsS3ITlAaq6qgkxmTzxyYat8bWhwHWp7yrbfIzI5LpesoGLQo9zY9gVXd8YHs8pwTfpLl/7HJlpcCSmqv1VdUATH+k0AJ6Ou4+2AHe5Ely/xEJguqoGvO0AXrsoY+rzf099SHU4CsAR8jFjg//2OSLTFlwQfIY+shaAzXTmhvB4nyMy0LSOHbvx+ig2xyhgk6qOSei/OA43CrwO4p03KklvxWiTZ15a+hXzP3LtpETgmoK7CYj1DDLN10FquDJ0T3z7ociJvBvt52NEBpqYxLx7Vc+LyAbvftgm78/PiciJTThlEbtPYi7BNRVOnJOWyvpgJk/tDEe4+skP4tsTSg7giICtmGMyZ2RwMSMCbvaQEuDKmnOJWpGHr5rSAPg63P2r0UA3YKX36IZbLHO+iNya5mkrge4J7zEUdylxfp3jimhC2ymTH25fuJJVG9y9ry4dQvzme4N8jsi0Rb8P3Us7XPupxTqQRyLH+xxRfksriXlzwKbiVnUep6qBhHtgsftWW4DJInJ6Q+eqYzFQnDCKuxx3KbE04b1H4pKYfbQ2u/li83ZuefHj+PalYwaxV+f2PkZk2qq+gbVMDD4d354ePpst2tHHiPJbuiOxybjkUqyqj9R9MqFjvHjHpuo67zXzvcuTZ+IuJb4I4C398gIJ98iMSXTVEx+wvca1kxrcc09+dFQfnyMybdkvQv9iPzYAsJ6u3BBOuVWsybB0k1gJUN5QibuqVuJaTh2Z6km9TvElwBJcMluMuzQZM8bbP01VH00zZtPGPfveFzz/wdr49tU/+BahYEZqloxJqpPs5LcF98W3742M4e2oXb72Q9b8T1fVclUd5l2iLKnTI3GUt/8G3wI0WWnzthp+/0RtMcdZRx7AUUV7+RiRyRenBN6MF3kATK2ZyI4aay7d2tJNYmW4e1cH1neAN0m5mAwumRJLaCJygXdp0RgA/vjMh6z72nVS2HfP9lz+/YN9jsjkCxH4Q8GddMYVE1Xq/rvclzWtI90kNhN3WW9eshZQXsuoRdQpykiV139xRD2PM4BpuFWljeHVj9czZ9Fn8e1rTjuUrh0LfIzI5Jv9ZSNTQw/Ft297pYKPvqh3pSrTAtLqYq+qc0XkEVzhRbmIVFE74irBlcULUJpOxw5v9DYfN4Jr8FB2L7s3eWh7dYTLH30vvn3yoT056dCePkZk8tWPggt4InIsb+tgwlFl6iPv8ujPjrX7sq0k7b9lVR2HK6VfhZsbNppd54yNU9V0S3Uux1U1rgRu8M6N9+cbgEepTY5jkp3A5Jc/z1vG6o21c8KuPu1bPkdk8lVAlOsKbo/PHXv3s83c9eoqf4PKI036qKCqc72+igGgP9A/Yc7YbqX3KRjrTqsDVHUaLikK8IKqTvMS54XAKBHZsykxm7bjnU+ruOM/tXU/vzvlEPbds4OPEZl8NyDwOReHagun/zRvGZ9s2OpjRPkj3cnOERF5LnGfqq7MwGrLRSRcJvTK9AGGJuybBWwijcU2TdtTE4ky9ZF3iXrtEI8bsBfjSnr7G5QxwOTgUwzu6T5j76iJcvmj76FqfTtbWrojsVW4kVemJWslVQkMr7OvHGsAnNdmvlLB0i/d8hcdCgJcd/rhiC3XbLJAgUSYMfZwAt4/x9cqNlCaUHhkWka6SWwc0F9ELs1wHJW4wpBE5exe6FGINQDOWx9/9TU3L0hoLTV6EH326uRjRMbs6vDehZz/7drO9n94+kPWbtnhY0RtX7pJbCMukV0hIm+LyHXeasxJH2mcdw7QzeuCH1vZeRFQFOvBKCL9cEu2WO/EPLQzHOGXDy6hOuLWCTu8d1fOO66vv0EZk8SvRw+iT3f34WrLjjC/nrOEaNQuK7aUtErscQlEcUUXw7xHsp+OePuDqZxUVWeIyARce6nNwATcnLTLgbkiEluCRbF7YnnphueWxefftAsFmDH2cCthNlmpY7sgM8Yeztmz30AVXv14A7f/p5JJ32mJOzEm3SQ2jeRJq9lUdZiITMSN9lDVzV7n+lLcfbgqYFaaK0abNuCV5eu4PaEa8YqTBzO4Z3PXYTWmBVzVFYCjgZ8FJnBr5DQAbnjmfY6ZfyaHBZpbA9fQe29uuXNnsXQnO89oqUC888+us11OyxSSmByx/pudXDrnnfj2iYP24afH9vUvIGNS9D+hubwaPZR3tD81hLik5hc81e63dJKdfofWpmTl9RgR6eK1oLKP23lMVZky913Wf+P+0+/duT03jDvCqhFNTiiQCH8tuIVOuMKOSt2fa8Ln+BxV25NSEvOSygivAe+IlkouInKZiGzAzQerwOuTKCLPi8j5LfGeJnvd8/onvLj0q/j2jeMOZ29b6NLkkL6BtVxTcFd8+6HICJ6J1J05ZJqj0STmVRmuBObhii3mARvTXLm5USLyPK5oI9a+KvHj9nBgVt2J1qbtWvrlFv74zEfx7fO/3Y/vDtrXx4iMaZozAws5NfB6fHtazUQ+1+4+RtS2NJjERGQoMJfaxLLA+xrAVQ3u1sm+KbyCjtG4xTCLVHVAnUP6AS8Co21E1vbtqIlwyYNLqA67cvqD9+vClJNswUGTm0TgjwV30It1AGxhD/6n+udE1C6LZ0JjI7HLva8zvb6GY7wEMxs3UspUuftkXNXjiGSrRqtqlaqOxpXfX5ih9zRZ6rpnPmLZ2tquHH87ewjtQynN1jAmK3WVbfyl3a0EcB/M3tRDuC3yA5+jahsaS2LFuMa8P0vcqaqTvT8emaE4ioFyVW1sIZ5FWMeONu3Jdz7nn69/Et/+31MPYcC+1vPZ5L4jA8u4KPhYfPvP4bG8FjnEx4jahsaSWBH1d8ioxLWByoRUu3BYAmvD3l+zmd/MrS2nH3NID344vI+PERmTWReHHmOYLAMgQpBf1FzCp9F9fI4qt6VSnZisOW9D+5tiMVAsIgfWd4DXdqqI2kU4TRuy/pudTLpnETtq3OWWor33sHJ60+aEJMrf293MPmwCYBN7MrHmUraqVd02VbbME7sOd49tXp1iEQXw+ikuwtpOtUnV4Sg/u6+Mzze7+TR7tg8x+6cldO1Y4HNkxmReT9nEbe3+El9Ec6n24dc1PyNqhR5NkhVJzOvMMQ0YAJR7c8UU12h4A1CGq5Ccoaov+hepyTRV5conPuDtVe6TqQjcfPZQ+u/T2efIjGk5wwIr+EPozvj289Hh3BzJ6KylvJEVSQziLa0G4ErpxXt0874uAIap6uX1n8HkovveXM2Db62Ob089aTAnDrb5YKbtGx96hfOCz8a3/xIey3ORTNXK5Y9UeicWi0ikvicbeE5VNd3ejJUkLHopIl1VNT+7WuaBNyo3cPUTH8S3TxuyP5O/Y7U7Jn/8NnQ/K7Q3/4keBsCva35GX/mSwYFPfY4sd6QyEpMmPpo9yrME1nZ9unEbP7+/nLC3ztJhvboy/Uxbpdnkl5BEuaXgZvrIWgC20YGJNZeyUW1aSaoaHCmpaotcbhSRfzTj5aqqP89YMKZleUtTJNqsezCx+n/ZqK58fm+qmLX+Ijr8cWNrR2eM7wplK7cX3Mjp1dewlY58qvtyYfWvuKfd9XSQGr/Dy3rprieWKbEOHbBrj8RUKGBJLEdt1facWz2FpV4CKyDMbe3+wn5iCczkr4MCa7ip4FYm1VwKwFt6MD+vuYSZBTdRUP/dHIN/SQxqV3+eh1v4sgVXizPZYIcWMLHmUhbrwPi+6wtmUxJY7mNUxmSHMcEyrtD7uTb8IwBejBbzPzU/568FtxCUFlmLuE3wK4mNxy2zMgpXyDGK2mRWmkL7KZNjajTIRTUX81r00Pi+q0N3c2ZwoY9RGZNdJoWeZot24hav3P6p6DHsEd7BdaHbCVgiS8qXEntVnes1Ew7gEtqLwBhcY+FNIvKciJxvi2K2DREVLq25kPnRYfF9vwk9xE9DL/gYlTHZ6dJQKecGa1edejhyIv8X/jFqOSwp3+eJeQlttCW0tkkVfhc+nyeix8X3/Sz4L34ResLHqIzJXiLw+9C9jA2+Et93V+Rkbgqf6WNU2cv3JJYoSUJ7lF0Tmq3wnENUlT+Gf8SDkRHxfecEX2BK6GEfozIm+wVEuT40m+8H3ozvuzlyJrPCp/gYVXbKqiSWyEto47yENg7XtWM0bnVpk+VUlT+9sJzbI7X/6c4ILOTq0D+xqWDGNC4kUf5ScAsnBJbE910b/hH/DI/xMarsk7VJLEZEzsCNykZ6u7I+5nwXjkS54rH3uOWlj+P7Tgq8xYyCmXZz2pg0tJMItxX8heHyUXzfleFzubFmnN0j82RlQhCRM0TkYa+lVSluJLYYN7+sm6/BmQZtr45w4X1lPPhWbduc7waW8NeCWwhJ1MfIjMlNHaWaO9rdyBBZEd93S+R0poQnUaO24nnWJDERGVFP4poGdFPVElWdba2ostfGrdX88PY3mP/RV/F9ZwQWMqvgT7SXsI+RGZPb9pTtPNDuWk4MLI7vK418l4k1l7Itz9ci8zWJ1Ulc80ieuG6wxJX9Pt24jbG3vcbi1bVrpV54Qn/+VPAP2lnHAWOarZPsZFbBnxkXfDm+7+XoEM6u/h0b8rjXoi+Tnb3eieOBQlznjnLgYWCuqlrnjhzzweebOfeut1n39U7AlQhfeeohnHtcP3izkRcbY1JWIBFmhGbRg03xCdHvaH/OrL6aezZso89enXyOsPVlQ+/Eebgk1h2YlEIXc1XVK1owNpOGl5Z9xS8fWMw3O93lwnbBADdNGMIph+/nc2TGtE0icFlBKT1kE1eGzyVKgFXakzP+8Sozzylh2IH5VTbgd+9EcGXzoxs6sA4FLIn5rCYS5U8vLOe2Vyri+/bsEGL2T0o4umgvHyMzJj+cE5rPvlLFxTUXsZN2rP+mmvEzX+c33xvEpOOLCATyYy6LX0lsnE/vazLgs03buPjBxZQn3P/q2aUDd//3kQzuaY1VjGkt3wsu4n65lguqL6WKPYlEleufXcrrFRv48/gj2Ktz2y/68CWJqeojfryvab4XPviS38x9l83ba9c5OuGgffLmP4wx2aYksJyn2v+WX/a4J15Y9crydXz/5oX89ayhbf7KiJ+XE01rSbIwZbp2aojrw2dzV+Tk+L4gEX4TephJnzxN4EabeWmMX3rLeuZMPoYbX1jGzFcqAVi7ZSc/nP0Gl4w8iItGDCDYRi8vZs08MZO9lkYPYGz1VbsksF6sY067a7gw9JR14TAmCxQEA1x+8sHcdd6RdN+jHQBRhZvmL+fHt7/J6g3bfI6wZVgSM/X6Rjvwh5ofcUr1tbynRfH9YwJv83T7KxgWWNHAq40xfjhx0L48c/HxDO/XPb7v9coNjL7pFf46fwU7atrWvE1LYmY3qvBE5BhG7ryR2yOnEMG1tikgzFWhu5lZcBOFstXnKI0x9enZtQMPXHAUF48cGG+4vTMc5ab5y/neX/7Ny8u+avgEOcSSmNnFx9H9+XHNFVxc80vWUvtJ7ujABzzd7nLODb1gXeiNyQGhYIBfjz6Ix35+HIf1qr0v/smGbZx719tceG8Za6q2+xhhZlhhhwFgk3bmtvCp3Bn5PjUJ/yz2YRO/K7ifHwRes+RlTA4ackAhj//iOB54azU3PLeULTtcY4LnPviSV5av4xcn9ufc4/rRuX1upoPcjNpkzFfalTvC3+feyGi20SG+P0iEnwaf51ehR+giuf9pzZh8FgwI5xx9ICcf2pPrn13K3LLPANheE+HGF5Yze+FK/vu4fpx7bF+6dirwOdr0WBLLU59rd2aFT+XByAh20m6X50pkGdcU3MUhgdU+RWeMaQl7d27PjeOO4KwjD+B3j7/P0i+/BmDz9hpumr+c2Qsr+ckxB3L+t/vlzLxPS2J5ZmW0J7MipzI38p1dLhsCDJLV/CL0L04NvGFl88a0YSV9u/PUL7/N3LLPuPXlClZvdOX33+wMc+vLFdz56kp+OPxALji+H/sXdvQ52oZZEssDW7U9z0SOojRyAm/pwbs9f7hUcFHocUYFyi15GZMnQsEAZw3vw9hhvXny3c/5+0sVfPzVNwDsqIly56srufu1lRw/cB/GlfRm9CE9aB/KvkU4LYm1UapK2SebmLPoU57eeStb2f3T1JGylItCj/OdwLtWtGFMngoFA5w+tDenHdGL5z/4kr+9+DEffrEFcJOlX1m+jleWr6NrxwL+a8j+jCs5gG/t34UUVhxpFZbEWtPnizPSAqo+URXe1SJejAzhqegxVOr+3jO1CSxIhO8G3mFi6CmODixtsViMMbklEBBOPmw/Tjq0Jy8t+4rbF67ktYoN8ec3b6/hn69/wj9f/4SD9+vCqYfvx4mD9uXg/fb0NaFZEstxW7QjC6OH82JkCK9Eh7Ce5Emyv6xhfPAVTg/+h32lKukxxhgjIowY3IMRg3vw6cZtPFL+GXPLPuOzTbVVyh99sYWPvtjCDc8vY7+uHfjuoH0ZMXhfjhuwF53atW5asSSWY7ZoR5ZEB1CuA3kjejCLooMI1/Nj7Mw2/l/wdcYFX2GofGyXDI1py1rgKs8BwK+Ai1V4o+BgSiMn8Gx0ODuorVz8YvMOHnxrNQ++tZp2oQBHF+3FUf26M7RPIUf0LmSPFp5/lhdJTESKgcuBUd6u+cB1qlrux3lSFdYAK3U/FntJqzw6kBXaC22g0crebOa7wSWMCCzmxMASOkp1S4RmjMkjAVGODX7IscEPuVrvZkG0mBcjQ3klegRb2CN+XHU4yr+Xr+Pfy9cBbn7aoB57UnxgIcMO7MaQA7rRp3unjHbUF9W2XY0mIqOAed5mpfc11s12nKrOba3zlOwf1EWTOu+2P6xuefEV2osV2pvl0d6s0F5U6v67lcEnc5hUcmJgMSOCSzhcKq3C0CQlV2+J/1mvbHzx0nSPN/knrAHKdSAvHnMPLy39imVrv270Ne1DAQbs25mB+3ZmYI89GbhvZw7qsScHNJDcRKRMVUuSPdemR2IiUkht4hmtqvO9/bGEVCoi3VS1wZtEmTrPN3TkofB3WaP7sEb35jPdmzW6N1/SPd5ktzFBIhwsqykOrKA4sIJjAx/YPS5jjC9CEmW4LGP4yYOZdvJgPt24jdcrN7B49SbKP6li+VdfU3ectDMc5YPPt/DB51t22V8QFPbr2pFehR3p1a32a+9uDc9Ta9NJDJjkfZ0RSzwAqjpfRGZ5z08CZrTGeVbqfkwLT2rokN3sxwa+FVjF0MAKimUFRwQq6SQ70zqHMca0hgO6d+KA7p0YX3IAAFt21PDOp1WUf1JF+epNfPjFFtZ9nfz3V01EWb1xW3zidaraehKb4H19OMlzM3GJZwKNJ7FMnade+7GBgYHPGChrOEg+Y2DgMwbIGutbaIzJWV06FHD8wH04fuA+8X2btlbz8bpvWL72a1as/YYVX33N8rXf1JvcGtPWk1gRQLLCC1Ut9+Y2FLfWefZgB6cHFtJL1tNb1tFL1tNL1rO/bKCD1KQQhjHG5LZue7TjyD26c2Tf7rvs31Yd5vOq7Xy2aTtrYl+9P3/SwPnadGGHiKtwUNWkdwsbez4T5xGR2KVGgEHAshTDz0V7A+v9DsI0if3scltb//kdqKr7JHsiH5JYlap2q+f5TUBhikms2edp60RkUX0VRCa72c8ut+Xzzy8fVnYubOJzLXUeY4wxGdLWk1hl44eQSn16ps5jjDEmg/IiiXmdNnaRsC+VBJWp87R1s/wOwDSZ/exyW97+/Np6Eiv1vk5I8lxs38xWPE+bpqp5+x8p19nPLrfl88+vTRd2QG3lIDAsViLvjZ7KYPeKQq8LB4mTmptyHmOMMS2vrY/EAMZ5X8tEZJ6IzMNLPAnPAbu0l5onIkXsKuXzGGOMaR1tPol5jXmHAXOBEu8xFzeiSqn5bybPY4wxmSAiY0WkTEQ2iYh6f57id1ytrc1fTjTGJNfaSwuZzPGS1XRvsxJXHR0vMlPV/r4E5oM2PxIz/rFPitnLu/dbBowFNnqPsbjL5WP9jM00zPvZTcclrmGq2l9VhwHdcB9EikQkbwrNbCRmWoR9Usxe3r3fTd5msqWFABpdWsj4Q0RKcR84dlvHMPFnmy/FZjYSMxlnnxSzXr1LC1E73yi9NYNMa4p9GJxf9wnvg0dsXmtedBKykZjJOPukmN1EpAz3i3BY3ftfCdNGyr0PHibLxBosJLt3mY//vyyJmYwTkQrc8jVJL0k19rxpWY01rE51dQeTfbypP6Nwo+ypfsfTGiyJmYyzT4rZLVNLFJns4V3Cn4n7cDhfVUf7HFKrsXtiJuNUtbyBMu1YC68mr4JtMqKhEbCNjnOEiBR5o695uAQ2I58SGNhIzLSSfP6kmG1sJNY21KkAngtMVdW8a0RuIzHTouyTYlaypYVynFfdOx0oB/qr6rh8TGAAIb8DMNnHK7xIR1WySjb7pJi1KnHTHIrrqU6MHWOykIhMwk2BsCsaWBIzSWRiIrL3SXES7pNi3n5KzFKluAq2CbifTyJbWij7xaoOrfE4dk/MtADvk+JM7JNi1rKlhXJTnW4rDV7yVdVuLR+R/yyJmYyzeWDZz+uPGKsUjXV+iDUC3m2SuskOiR80GpMvH0QsiZmMsk+KucO62Ju2wJKYySj7pGiMaU2WxIwxxuQsmydmjDEmZ1kSM8YYk7MsiRljjMlZlsSMMcbkLEtixhhjcpYlMdPqRGSsiKj3aLS9UbrHm+zmNYUu9DsO0zZYEjN+G5/CMZNbPArTKkSkCKgAZvsdi2kbLIkZvxV6a401pLHnjTF5ypKY8VOsZ1+9Iy2vx1/iscYYE2dJzPip3HuMbeCYWIIrre8AESkWkVIRqfDum5V5nfTrO36Kd8wm7/gKEZme7D5NknNXiMjMusd6+zRhPa7E50Z5z01P2DfF21fkvUdFQmf5tL6vhHMVe/cP53nfW0XsPUWk0Isx8VxJ/97T+fus897F3nur9/7zEv8+vMVRY2vVxe5zTk923nR48WmyEX1Dz5m2wZKY8dvD4H7R1/P8KNwoLOl6ZN7Cm2W4RFiFS4rFwEzvl2bd48twC3UWAYu8cxcBU4AFdY6N9YGMnXs+0B23TlpK/SFTkNhrMj7aTPf78kymNtkvwvu+vGKYMtz9x9gHh2KgtG7SbeL7gvs5lXnvORfYGNvn3QcDtzzPLO/Plbh1sRo6Z6qu875OTdwpIqVe7JNV1UbybZWq2sMerfrA/YJUapOJBcN70AAABZNJREFUAqUNHDcJ9wtRgZkJz8deWwEUJewvxP1CVWBKwv5J3r55Sd6rwnsu8Tyxc4yqc+w8b//YhH0zvX3FSc4di316wr4p3r5NiTE28fuKnWuX9094X/VeV5gk3ulNfd8k7133udIG3mO3n3cz/01tSvz5JcQ1M5Pvk8F4d/k7Tvh73pT4c7JH4w8biRlfqVvxub5LirFLiXPqeXms3H6XlaPVrWEWW/X28oTjN+JGCbt8YvfE1s8qSthX7J2v7qf4qd4jE6tVV6rqjDr70v2+YmZpwjIqXtyx5XAm6q5ru8VGbInfb1PfF6C8ge+jqO7BLSA+GvNG9dNxi7Lucr/Vu8Tp6zQN7zJuEe7fY6LxuARma/ClIeR3AMbgLikWi8hY3XUxxlG4X0RVIklXbSkBqjTJ+leqWikilUCRiBSqapV37t0We/QudyVLouVeXKUkrLPlfc3UmlvJLnOl9X0lPJ3sEudG3C/GuudKloCb+r7gXRZO4T2S8i5j7pXCoW9r8gU7Z+ES1yTvUa7JVxUvryfW1jQa9+Gl7t/haKyAKW2WxEw2iP0CmoCXZBKKDuot6MBdfqFuQUQS3fFGJF5BxiTgSNxIq6FRwjjcpcOxuEKEKty9plJgToY+Mb+dZF/a35en7if7dDX1faH5o9LLY+/fiKQfRLwPOrNwP9sqYGSyF6tqslF4a4vd5022f1aS/aYBlsSM77xfQHUvKU7wvia9lJhQHVhF7aWk+mz0XlOMK94oxP3SLcdd8irHfQqeUieuSqC/d3lqHO6XTOwxXURGJhu1NEdTvq9cft8YbeYq31788eKgZB8wvJ/jZFUdV/e5JMdWNHZMHVWqOizFOIuoU9DiXQ0orLvfNM6SmMkWdS8pjsW7lJjsYC/xVXl/rnsvpj6luF8U4+pekhKRZJeeYu81H++Ts/dL6HJcwisF+qfwvinfE2ri99Vsfr1vBi3A/T3HLnnWvTQN7oNKSj8LVU3l59oUsURbdyQ2yntfu5yYJivsMNkidhllQoqXEsFd2iusZ25WYWyuUsLuIty9iGT3VHb55ebN36qoWwTg3Vubmuw1nu5J9jX6yb+OdL+vTPHrfZHaOXuNPXb7N5FYSk/DBSjFuO/RT/XdDxuGd5/V+7fXGsUwbYKNxExWqHNJMXZpq76qxJipuGKGBSIyLLGijtpRV+IvvSrcp/SixGO9ibyxxFnoxVPp/SKZJCKliZ+QE5Js4qXE2OWnyew63ys2PSAd6X5fmdKa71v3/td1pFjYkbjhFYSMxVVmzvL2xQpyiutc7i2htmLSL6PY/X4iuMrE2L/3sTk6GvaH3zX+9si/BwnzxOrsj83j2m0uF0nmiXn7pye8pgJ3T2FTPedIPHYe7hdyhXd8acI5JnnHT0ly7oqEfYlzsorqHFtKwtyqut9vwr6x9fwdpfN91XuuWLxJ9iedr5XO+6bw3ru9By55xc5fWt/3n+a/o7p/H2MbeN+ipr5fBv7dx7/3JH+HFXjz7Or+v7BHww+7nGiySeLIK6VP/Oou7cVKk7vjPm1X4m7gj05y7GTv+VG4y0vzgX7qbvbHuncM846fgbs8Ve6dOzaimgX0113nZFV6r4vFERutjfbOk1Z3inS+r0xq6fdVdxltBm40EpsvlTbvkmepF9sul2vVXS6uwlWUxkZ8sXtOmZjb11Sxfz8zvBZgU7yR5Czcv8sjvW2/R4s5RbzMb4wxbZbXo3GUplBB2IIxzPRiaKmikbxkIzFjTD7IhqKO+uaHmWawJGaMyQclZK5pc9rqmx9mms+qE40xbZ42czJ1BpR4X20klmE2EjPGmJa3CFfkY819M8wKO4wxxuQsG4kZY4zJWZbEjDHG5CxLYsYYY3KWJTFjjDE5y5KYMcaYnGVJzBhjTM76/7MMBrfVj/KsAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| }, | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Black bars show confidence interval based upon N measurements\n", | |
| "Use the slider to change the number of samples, N\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "def convergence_of_mean(N):\n", | |
| " x=np.linspace(-3,3)\n", | |
| " y=norm.pdf(x,0,1) # convert fraction to percent\n", | |
| " plt.plot(x,y)\n", | |
| " y_sample=np.random.normal(0,1,N)\n", | |
| " tstat=t.interval(0.95, N-1)\n", | |
| " CI=tstat[1]*sigma/np.sqrt(N)\n", | |
| " plt.hist(y_sample,density=True)\n", | |
| " plt.vlines([np.mean(y_sample)-CI,np.mean(y_sample)+CI],0,0.4)\n", | |
| " plt.xlabel('Measurement=$x_i-\\mu$')\n", | |
| " plt.ylabel(r'Probability of\\\\ Measurement (\\%)')\n", | |
| " plt.axis([-3,3,0,0.5])\n", | |
| " \n", | |
| "interact(convergence_of_mean,N=(10,1000))\n", | |
| "print('Black bars show confidence interval based upon N measurements')\n", | |
| "print('Use the slider to change the number of samples, N')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The interactive plot above generates `N` normally distributed numbers and plots the histogram with two vertical lines to indicate the confidence interval of the measured mean. More measurements, larger `N`, produce a smaller interval, even though the standard deviation is constant. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Problem #2\n", | |
| "\n", | |
| "Try creating 100 normally disributed random numbers with a true average of 10 and true standard deviation of 1. What is the confidence interval for the measured average? Is it higher or lower than the confidence interval for 20 samples?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "N=100\n", | |
| "data # =... enter your work here\n", | |
| "tstat# =... enter your work here\n", | |
| "avg # =... enter your work here\n", | |
| "std# =... enter your work here\n", | |
| "conf_interval = tstat[1]*std/np.sqrt(N)\n", | |
| "\n", | |
| "p.check_p02(avg,conf_interval)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### 2 - Difference between means\n", | |
| "We will use the t test for independent samples. Each group can have different number of measurements. We assume that the population has a normal distribution if we could make infinite measurements. The student t-test tests the null hypothesis. \n", | |
| "\n", | |
| "*null hypothesis for given data: there is no difference between the two sample means*\n", | |
| "\n", | |
| "The t-test is based upon the means $\\mu_1$ and $\\mu_2$, standard deviations, $\\sigma_1$ and $\\sigma_2$, and number of data points, $N_1$ and $N_2$. The t variable for sets of independent variables based upon the difference in means between sample sets as such\n", | |
| "\n", | |
| "$t=\\frac{|\\mu_1-\\mu_2|}{\\sqrt{AB}}$\n", | |
| "\n", | |
| "$A=\\frac{N_1+N_2}{N_1 N_2}$\n", | |
| "\n", | |
| "$B=\\frac{(N_1-1)\\sigma_1^2+(N_2-1)\\sigma_2^2}{N_1+N_2-2}$\n", | |
| "\n", | |
| "Let's use the same random numbers from the previous examples and examine two sets of data and compare to the t-distribution table.\n", | |
| "\n", | |
| "**df=18**\n", | |
| "\n", | |
| "|p=0.05 | p=0.025 | p=0.01 | p=0.005 |\n", | |
| "|--- | --- | --- | --- |\n", | |
| "|2.10|2.45|2.88|3.20|" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "t=1.38\n", | |
| "df=200\n", | |
| "| p=0.05 | p=0.025 | p=0.01 | p=0.005 |\n", | |
| "| --- | --- | --- | --- |\n", | |
| "| 1.97 | 2.25 | 2.59 | 2.82 |\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "N=200\n", | |
| "data1=np.random.normal(10,1,N)\n", | |
| "data2=np.random.normal(10,1,N)\n", | |
| "\n", | |
| "mu1=np.mean(data1); mu2=np.mean(data2)\n", | |
| "sigma1=np.std(data1); sigma2=np.std(data2); \n", | |
| "N1=N2=N\n", | |
| "\n", | |
| "A=np.abs(N1+N2)/(N1*N2*1.0)\n", | |
| "B=((N1-1)*sigma1**2+(N2-1)*sigma2**2)/(N1+N2-2)\n", | |
| "\n", | |
| "tstat=np.abs(mu1-mu2)/np.sqrt(A*B)\n", | |
| "\n", | |
| "print('t=%1.2f'%tstat)\n", | |
| "\n", | |
| "df=2*N-2\n", | |
| "\n", | |
| "# Print out the table for df degrees of freedom (N1+N2-2)\n", | |
| "print('df=%i'%N )\n", | |
| "print('| p=0.05 | p=0.025 | p=0.01 | p=0.005 |')\n", | |
| "print('| --- | --- | --- | --- |')\n", | |
| "print(\"| %1.2f | %1.2f | %1.2f | %1.2f |\"%(t.interval(0.95, df)[1],\\\n", | |
| " t.interval(0.975, df)[1],\\\n", | |
| " t.interval(0.99, df)[1],\\\n", | |
| " t.interval(0.995, df)[1]))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Our t-statistic is lower than even the p=0.05 confidence interval. This indicates that we cannot reject the null hypothesis. Based upon the current data set, there is no measured difference between the two sample averages. This makes sense because the two sample means should be the same. We just used two random data sets with the same average and same standard deviation. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Problem 3\n", | |
| "\n", | |
| "Try changing the mean between samples. What happens to the value of tstat?\n", | |
| "\n", | |
| "At what point does the t-statistic become statistically significant? Keep the standard deviations, $\\sigma=1$, and vary the averages for `data1` and `data2`, shown above as `10` and `10`. With 200 samples, can you reliably measure the difference between 10 and 11? or 10 and 10.1? \n", | |
| "\n", | |
| "Create two datasets, `data1` and `data2` that *are* different, statistically with p$\\le 0.05$, but have means within 10% of each set. \n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "data1 #=\n", | |
| "data2 #=\n", | |
| "p.check_p03(data1,data2)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Data to analyze - Marble Madness Design Problem\n", | |
| "\n", | |
| "Your company manufactures marbles with glass. There was a mix-up with the glass delivery and all of the affected marbles are at higher risk of shattering. The faulty glass has a slightly different density, $\\rho_{faulty}$=2.80 g/cm$^3$, than your normal manufacturing glass, $\\rho_{mfg}$=2.52 g/cm$^3$. \n", | |
| "\n", | |
| "You take 15 samples from three bins: A, B, and C. You know that bin A is the correct glass, but bins B and C might have the faulty glass. Write a report that clarifies whether bins B and C are safe to ship using the student t-test and creating confidence intervals for the three bins' samples. \n", | |
| "\n", | |
| "|mass(g) A |mass(g) B|mass(g) C|\n", | |
| "|---|---|---|\n", | |
| "|34.6| 40.3| 56.1|\n", | |
| "|40.7| 35.0| 59.3|\n", | |
| "|37.5| 43.7| 51.4|\n", | |
| "|45.8| 43.2| 50.2|\n", | |
| "|41.4| 41.1| 46.4|\n", | |
| "|44.2| 42.3| 56.1|\n", | |
| "|44.5| 39.8| 43.1|\n", | |
| "|51.8| 39.2| 55.1|\n", | |
| "|47.5| 49.4| 46.9|\n", | |
| "|45.4| 43.2| 39.9|\n", | |
| "|36.4| 44.4| 51.8|\n", | |
| "|46.2| 44.8| 46.2|\n", | |
| "|43.0| 44.5| 49.7|\n", | |
| "|43.3| 50.1| 54.9|\n", | |
| "|42.0| 47.1| 64.5|" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Your Report \n", | |
| "\n", | |
| "1. Introduction\n", | |
| "\n", | |
| " a. What is the problem and how can you solve it? \n", | |
| " \n", | |
| " b. What is the student t-test?\n", | |
| "\n", | |
| " c. What assumptions are used in the student t-test?\n", | |
| "\n", | |
| " d. What do the results of the t-test tell you?\n", | |
| "\n", | |
| " e. Why would you use the student t-test?\n", | |
| " \n", | |
| "2. Methods\n", | |
| "\n", | |
| " a. Description of materials and measurements\n", | |
| " \n", | |
| " b. Description of t-test procedure\n", | |
| "\n", | |
| " c. If someone repeated the experiment he/she would read only this section\n", | |
| "\n", | |
| "3. Results and Discussion\n", | |
| "\n", | |
| " a. Present experimental results succinctly\n", | |
| "\n", | |
| " b. Present propagation of error clearly\n", | |
| "\n", | |
| " c. Interpret results clearly\n", | |
| "\n", | |
| " d. Compare the averages and confidence intervals for bins A, B, and C\n", | |
| " \n", | |
| " e. Describe the student t-test results, how confident is you rejection/non-rejection of the null hypothesis\n", | |
| "\n", | |
| "4. Conclusion\n", | |
| "\n", | |
| " a. Is this a good method to use for this application?\n", | |
| " \n", | |
| " b. How could we make it better?\n", | |
| "\n", | |
| " c. What other applications would benefit from this technique/these results?\n", | |
| "\n", | |
| "\n", | |
| "\n", | |
| "**Links of Interest and References**\n", | |
| "\n", | |
| "\\[0\\] Student. (1908). [The probable error of a mean](./student_error-of-mean.pdf). Biometrika, 1-25.\n", | |
| "\n", | |
| "\\[1\\] ['Student's' t Test (For Independent Samples)\n", | |
| "](https://www.ruf.rice.edu/~bioslabs/tools/stats/ttest.html)\n", | |
| "\n", | |
| "\\[2\\] [t-based Confidence Interval for the Mean](http://www.stat.wmich.edu/s216/book/node79.html)\n", | |
| "\n", | |
| "\\[3\\] [Confidence intervals,\n", | |
| "t tests,\n", | |
| "P values](http://evolution.gs.washington.edu/gs560/2011/lecture3.pdf)\n", | |
| "\n", | |
| "\\[4\\] [Student's t-test](https://serc.carleton.edu/introgeo/teachingwdata/Ttest.html)\n", | |
| "\n", | |
| "\\[5\\] Sehgal, J., & Ito, S. (1999). [Brittleness of glass](./brittleness_of_glass.pdf). Journal of non-crystalline solids, 253(1-3), 126-132." | |
| ] | |
| } | |
| ], | |
| "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.8.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |