{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Hypothesis Tests (Parametric)\n", "====\n", "\n", "## Unit 10, Lecture 1\n", "\n", "*Numerical Methods and Statistics*\n", "\n", "----\n", "\n", "### Reading\n", "\n", "Langley: Pages 137-189, 199-211, 230-245\n", "\n", "----\n", "\n", "#### Prof. Andrew White, April 3 2018" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Goals:\n", "---\n", "\n", "1. Understand the meaning of a null hypothesis\n", "2. Be able to construct a null hypothesis\n", "3. Be able to calculate a p-value, understand its meaning, and test it against significance level\n", "4. Visualize the intervals used to compute p-values for $z$ and $t$ tests" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from math import sqrt, pi, erf\n", "import scipy.optimize" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Hypothesis Testing\n", "====" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Hypothesis**: Going to college increases your starting salary. How do you prove this?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "It's not really possible to prove this directly. We can, however, disprove the opposite hypothesis. We construct the opposite hypothesis to what we're interested in, called the *null hypothesis*. The null hypothesis is an assumption of no-difference and/or no-correlation." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "In our example, this seems simple at first: *Going to college has no effect on your salary*. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "But, maybe the null hypothesis is: *For People who qualified and were accepted to college, attending college has no effect on their salary*. \n", "\n", "Or it might be *People who can afford, are smart enough, and motivated enough to go to college but did not, have the same salary as those that did.*" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Let's assume we know well enough what our null hypothesis is. Hypothesis testing is the ability to use statistics to disprove the null hypothesis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Null Hypothesis**: The simplest explanation, typically meaning no correlation or that all data is from the same population. Because a simpler explanation is preferred to a more complex one, the null hypothesis should be our default belief. We construct our null hypothesis as sort-of the opposite of what we want to study. For example, if we want to know if a sample is significantly different than our population, our null hypothesis is that it *is not* significantly different and then we aim to disprove the null hypothesis. Hypothesis testing is about showing significance by disproving or rejecting the null hypothesis." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Hypothesis Testing\n", "----\n", "\n", "We construct a null hypothesis and take it to be true. For example, we believe college has no influence on income. This allows us to build a simple probability model. For colelge then, we might take income to be normally distributed. Then we see how likely our *observed* data is according to that null hypothesis model. For example, we check to see if the sample mean of people who graduated from college matches our null hypothesis mean." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Hypothesis Test Example\n", "----\n", "\n", "I open up a cash4gold store and people bring me their jewlery. I know the probability distribution for gold melting is normal with mean 1060 $^{\\circ}$C and my measurements have a standard deviation of 3$^{\\circ} $C. I melt a sample at 1045 $^{\\circ} $ C and want to know if I should be suspicious." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Null Hypothesis**: The sample is gold" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's see if I can disporve this. If the sample is gold, what is the probability of that measurement? Zero, because it's a single point. Let's ask though how big of an interval would we need to include that measurement." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.999999426697\n" ] } ], "source": [ "from scipy import stats as ss\n", "\n", "Z = abs(1045 - 1060.) / 3\n", "interval_p = ss.norm.cdf(Z) - ss.norm.cdf(-Z)\n", "\n", "print(interval_p)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8U1XaB/DfSZo0adI9Tfe9ZRMBARcQtSg67q8K7uK+\njPr6DuM4jjI6476NjI6jIAoKoiIKrogiW9lBWihbKbR03/e9TbPc94+bFFpa6JLk5uQ+38/nGEhu\nkodafj2cnIUJggBCCCH8UkhdACGEkOGhICeEEM5RkBNCCOcoyAkhhHMU5IQQwjkKckII4RwFOSGE\ncI6CnBBCOEdBTgghnPNxx5sYDAYhISHBHW/Vt+aj4m3ASOlqIIRw46g9MkZKHBmZmZm1giCEnek6\ntwR5QkICMjIy3PFWfVufJt7OSJeuBkIIN9LSxNv0dCmrABhjRQO5joZWCCGEc27pkUtu7HNSV0AI\n4chznEWGPII8YobUFRBCODKDs8iQx9BKQ5bYCCFkALKyxMYLefTIM+eIt/RhJyFkAObYI0PqDzsH\nSh49ckII8WIU5IQQwjl5DK0Q0kuryYKNOdUormuDQe+L6aOMCA/QSF0WIUNCQU5k57t9pXjxx2w0\ndpi77/NRMjx4URKeunwEfJT0D1XCF3kE+fjXpK6AeIi31x7F+5vykByux93TExFj8ENtiwkbD1Th\nw/TjOFTWhE/vPRcqCnNZe42zyGCCILj8TSZPnixIukSfEACfbCvAS6uzMWWEATdNiYVSwXo8vjW7\nGt/tLsWMs8Lx8V2TwBjr55UIcQ/GWKYgCJPPdJ08uh01O8RGZGtfcQNeXXME4+KDMLOPEAeAi8YY\nceU5kVh/uAofby+QoEriKXbsEBsv5DG0sn+ueEvzyGWp02zFn1dkIchPhdumxUPRR4g7zBgfgfyq\nVvzrlxzMGGVEkkHvxkqJp5hrjwyaR06Ih1i8rQCFde24dVocNGrlaa9VMIZbL4wHYwx/+/agmyok\nZHgoyIlXq2kx4YNNeTg7LhApkQEDek6wXo3Lx0dgT3491mZXurhCQoaPgpx4tfc25MJkseG6ydGD\net7FY4wI0Kowb90xuGNCACHDQUFOvFZNiwkrMkpwbkoIDIGDW+yj8lHgsnHhOFbRgl+pV048nDw+\n7Jz0rtQVEAks2VEAs8WGS8eGD+n5F4wwYMOBKvx3Yx6uOivSydURT/YuZ5Ehjx558ASxEdloM1mw\nbGcRxiUEDbo37qDyUWDamDBklzUjs7jByRUSTzZhgth4IY8gr1wvNiIbP+4vR3OnBdPPGlpv3OGC\nEaHwUTIs3JrvpMoID9avFxsv5DG0cugV8ZZOCpKNr34vRmSwBrFhfsN6Hb1GhXMSg5F+pBq1rSYY\n9L5OqpB4slfskcHLSUHy6JETWTlS0Yz9pU24YITBKcvsLxwVhi6LDct+H9CB5oS4HQU58Tpf/V4M\nlZJhYnKIU14v1uCH8CANVu+voKmIxCNRkBOvYrJY8d2+MoyLD4LO1zkjh4wxTE4OwfGqVhwob3LK\naxLiTBTkxKtsOVaL5k4Lzk1xTm/cYWJSCBiAL/cUO/V1CXEGeXzYed5CqSsgbvLT/nLoNT5IHuBy\n/IEK1quREumP9YerYL3eBqWC+kDebCFnkSGP78aAkWIjXq2jy4r1R6pwdnxQn9vUDtek5BDUtXRh\n8/Fap7828SwjR4qNF/II8tKfxEa82sacarR3WTEpKdglr39WXCAUDPghq9wlr088x08/iY0X8hha\nyZkn3sZcJ20dxKV+2l+OQD8VEoyu2UNc5+uDlEh/7MithdVGwyvebJ49Mq7jJDLoO5F4hU6zFZuP\n1WBsXOBpD44YrnHxQahpNmF3ES3ZJ56Dgpx4hR3Ha9FhtuLsuCCXvs/YuCAwAN9nlbn0fQgZDApy\n4hXWZVdDo1IgKcK1R7MF+KmQGK7HtmO1tDiIeAwKcsI9m03AxpwqjIwOgI/S9d/SY+MCUd7Qgeyq\nZpe/FyEDIY8PO6csk7oC4kKHyptQ1WzC5ePds2f46JgA/LinDGsOVeKsiEC3vCdxr2WcRYY8euS6\nWLERr7Q+uwqMASNinLsIqD/GQA1C9GpszaX55N4qNlZsvJBHkBetEBvxSuuPVCPJqIde455/YDLG\nMDomEDllzWjsMLvlPYl7rVghNl7II8hzF4iNeJ3q5k5kVzRjdIx7hzhGxwSgy2LDb0foPE9vtGCB\n2HghjyAnXmtbnji8MTLa363vmxLpDx8lw4aj1W59X0L6QkFOuLY1txb+Gh9Ehmjd+r5qHwVSIvyx\nt6ABNpqGSCRGQU64ZbMJ2Jpbg9QofyiccBLQYI2OCUBNswkHy2kaIpEWBTnhVk5lC2pbuzA62j2z\nVXobESUO56zPqZLk/QlxkMc88mkrpa6AuMC2vBoAQHKke8fHHYyBGgRoVdiVXwdcJkkJxEVWchYZ\n8ghyjUHqCogLbM2tRWSwBkE6tSTvzxhDapQ/jpQ1o9NshUallKQO4nwGziJDHkMr+UvERrxGp9mK\n3QX1GBElzbCKw4hIf7R2WrC7qF7SOohzLVkiNl5QkBMu/V5Qjy6LDaMkDvJU+zj5xmM0DdGbUJAT\n4gY78+ugVDAkhOskrSNIp4Yx0BcZBbQ/OZEOBTnh0q7jdYg1+MHXA8alUyP9kVvRgqZOWq5PpEFB\nTrjTZrLgYFkTUly89/hApUaJy/W32mfREOJuFOSEO5lFDbDYBCRHSDPtsLeUCD0YgC20GyKRiDym\nH6atkboC4kS78uugYECCUdrxcQc/Xx9EhWiRVdIodSnESdZwFhny6JH7+ImNeIXdBfWIM+g8Ynzc\nISlCj4KqVtrW1kv4+YmNF/II8mPzxUa4195lwf6SRqRItJqzP8nhepitArbn0/CKN5g/X2y8kEeQ\nF38tNsI9x/h4SoRnDKs4OA593n6cgtwbfP212HghjyAnXmN3fj0UDIgzesaMFQe9RoXwQA32lzRJ\nXQqRIQpywpVd+XWINeg8cl+TpAg9jle2oNVkkboUIjMU5IQbHV1WZJU2esz88d6SI/ToNNuwi/Zd\nIW5GQU64sbe4ARargBQPmT/eW1K4+APGcfwcIe4ij3nkM9KlroA4we6CejAGxHvI/PHegnRqhPqr\nkVVM+67wLj1d6goGh3rkhBuZRfWIDtFCo/a88XGHpHA9cita0GG2Sl0KkRF5BPmRt8VGuGWx2rC3\nuBEJHjZbpbfkCH+0mazILKVeOc/efltsvJBHkJetFhvhVk5lCzq6rEj28CB3jJNvpX1XuLZ6tdh4\nIY8gJ9zLKBRngsR56Pi4Q6i/GgF+KuyjcXLiRhTkhAsZRQ0I1qkQrJfmfM6BYowh0ahDbmUrrDZB\n6nKITFCQE48nCAL2FNZ7/Pi4Q4JRj4bWLhytbpa6FCIT8ghypVZshEtljR2oajZ1jz97Osf2utuO\n10lcCRkqrVZsvJDHPPLpv0hdARmGzCJxvNlT54/3FhPqB5WSIaOoAQ9fKHU1ZCh+4Swy5NEjJ1zL\nKGyAr0qByGA+ukhKBUOsQYec8mYIAo2TE9eTR5AffFlshEsZRfWID9NBqWBSlzJgieE6lNa1o6rV\nJHUpZAhefllsvJBHkFdtEBvhTkunGUcrW7gZH3dIMOphE4CdBTROzqMNG8TGC3kEOeHWvuJG2AQg\nkZPxcYeEMLHe3QW0EyJxPQpy4tEyihrAGBAbxleQ6zQ+MAb64lApHTRBXI+CnHi0zKJ6RAdrPfIg\niTNJNOqRX92K9i46aIK4ljyC3DdUbIQr3RtlcTY+7pAYrke7yYqsMuqV8yY0VGy8kMc88otWSV0B\nGQJeNsrqj2Nh0I78WkxN5CgVCFZxFhny6JETLvGyUVZ/wgJ8ofNVYl9xo9SlEC8njyDPelZshCvi\nRllqj98oqz+MMSQY9citbIGNFgZx5dlnxcYLeQyt1O6UugIyBBlFDd3DE7xKMOpwuKQJBfVtSA7l\nc4hIjnZyFhny6JET7pQ3dqCyqZO7hUC9JdrH97cdp4MmiOtQkBOP1L1RFmfzx3uLMfhBqWDYU0gH\nTRDXoSAnHimzqAFqHwUiQ/jYKKs/ah8FokO1yC6jDbSI68hjjNwvRuoKyCDtLWpAnL03y7tEox47\ncmrQ1GlGkJbPD27lJoazyJBHkE/9XOoKyCB0dFmRXdGM6WPDpS7FKRKMOmw+XI1dhfW4cnSE1OWQ\nAfics8igoRXicQ6UNsJiE7ifseLg+MBzN+2ESFxEHkGeOUdshAuZ9hPoedsoqz8BfiqE6NXYX0JL\n9XkxZ47YeCGPoZWGLKkrIIOwt6gR4YEa6DXe8+2ZYNQhr6IVZqsNKqU8+k88y+IsMug7ingUQRCQ\nWVSP+DA/qUtxqkSjHs0dZhypapG6FOKFKMiJRymsa0dDu5n7hUC9Ocb7t9PCIOICFOTEozgWAsV5\nyfi4Q2SwFr4+iu4/HyHO5D2DkKfjP0LqCsgA7S1ugFathDFII3UpTqVQMMQbdcipEBcGMcb//Hhv\nNoKzyJBHkJ//kdQVkAHKLGpAfJgOCi8MugSjDuv2V6K61YRwf+/6QeVtPuIsMmhohXiM5k4zjlW1\ncHfQ8kAlGPUQBPGgCUKcSR5BvvthsRGPtr+kEYJwYgGNt4kP04EB2E0baHm8hx8WGy/kMbTSckzq\nCsgAZBY1gEHcMdAbadVKRARrcLCUFgZ5umOcRYY8euSEC5lFDYgK0UKjVkpdisskGPU4XtWKLotV\n6lKIF6EgJx7BZhOwr7iR+/3HzyTBqENnlxWZJXSOJ3EeCnLiEXKrW9FqsnjdQqDeHAuDduTTBlrE\neeQxRh48QeoKyBk4FsrEetnS/N4M/r7Qa3ywt5g+8PRkEziLDHkE+aR3pa6AnEFmUQP8NT4w+PtK\nXYpLMcaQYNThaEULLQzyYO9yFhmDHlphjOkYY977aRSRxN7iBsQbdbIItgSjHrXNJhTWt0ldCvES\nZwxyxpiCMXYHY+xnxlg1gBwAFYyxw4yxfzHGUl1f5jDtuEtsxCPVt3WhoLbNa+eP9+ZY8LTtOI2T\ne6q77hIbLwbSI98EIBnAswAiBEGIFQTBCOAiALsAvMEY8+w/cnup2IhH2msfH/eWE4HOJCZUPIt0\nT1G91KWQfpSWio0XAxkjnyEIgrn3nYIg1ANYBWAVY0zl9MqIbOwtboBSwRAd6t0fdDqofBSIDfXD\n4dJmqUshXuKMPXJHiDPG3mX9DGD2FfSEDFRmUQNiQrVQ+8hnNmyCUYei2ja0mOivDhm+wfzNaQXw\nI2NMBwCMsSsYY9tdUxaRC7PVhv2l3r8QqLcEow4Wq4BdBTS8QoZvwNMPBUF4jjF2B4B0xpgJQBuA\nZ1xWmTMZpkhdAenHkYpmdJptXr8QqLcE+we7uwrqcPmocImrIb1N4SwyBhzkjLHLADwEMcAjATwg\nCMJRVxXmVBNel7oC0o/uhUAGefXIA/xUCNGrkUVL9T3S65xFxmCGVv4O4HlBENIAzAKwgjF2qUuq\nIrKRUdiAUL0awXq11KW4XWK4HrkVLbDabFKXQjg34CAXBOFSQRC22X99EMBVAF5xVWFOtXWm2IhH\nEQQBvxfWy2baYW8JRh2aOyzIrmyRuhTSy8yZYuPFQBYE9TdTpQLAZae7xmOY6sRGPEpxfTtqWkxI\njvCXuhRJJHQvDKITgzxNXZ3YeDGgBUGMsScYY3En38kYUwOYwhhbCuAel1RHvNqeQnktBOotMkgL\nXx8FMopoAy0yPAP5sPNKAPcDWM4YSwLQAEAL8YfAbwDeEQQhy3UlEm+1p6AeOl8ljEHyPIhYoWCI\nN+qQU95MG2iRYTljkAuC0AlgPoD59hWcBgAdgiDQx+1kWPYUiePjChkHWIJRh3X7K1HZYkJkgDx/\noJHhG/CHnYyxqwBsBZAO4CPG2AWuKsrpwi8TG/EYda0m5Ne0ITlcnuPjDolGPQQB2J5P4+Se5LLL\nxMaLwexHPh/AXQCyAUwC8DZj7ANBEJa7pDJnOvt5qSsgvTjGxxPD5Tk+7hAXpgMDsLugHrMmxEhd\nDrF7nrPIGEyQVwmC4FiSv54xthPAbgCeH+TE42QU1kOllM9GWf3RqpWICNbgYGmT1KUQjg1mQVAh\nY+wV+2wVADAD4GMC7KarxEY8xp7CesSH6eCjlM9GWf1JCvdHflUrbaDlQa66Smy8GMzfIgHATQBK\nGGPbAORB3HfF8w+WsHaIjXiE9i4LDpc3dx+wIHfJEXp0WWy0gZYH6egQGy8Gs2nW7QDAGNMAGAtg\nvL0tYowlCYIQ65oSibfJKm6ExSbIbqOs/ji+DlvzamkDLTIkgz582T4dMcPeCBm0PYUNYADiZHK0\n25kE+KkQFuCLfcW0MIgMDQ1QErfbU1iP6BAttGo6w9shOUKP3MpWdJqtUpdCODToHjmXoq+VugJi\nZ7HasLe4AROTQqQuxaMkReix61gdfi+qx8UpYVKXI3vXchYZ8gjy0U9JXQGxO1jWhPYuK1IiaVjl\nZI6FUVvyainIPcBTnEUGDa0Qt9qVL87MSKQPOnsI1qsRoldjL22gRYZAHkG+Pk1sRHK78usQGaSB\nv1YldSkeJzlCj6PlLTBb6aAJqaWliY0X8ghy4hHMVhv2FNYjKYJ6431JCtejzWRBZgn1ysngUJAT\nt3GMj4+IDJC6FI/kOGBjS26NxJUQ3lCQE7fZeVw8ciVB5htl9SfUX40APxUdNEEGjYKcuM2u/DpE\nBWtpfLwfjDEkR+iRU0bj5GRw5DH9MO4WqSuQPbPVhozCBkxKofnjp5Mcrse+/AZklTbi3Hj6Wknl\nFs4iQx5BPuIxqSuQvQOljegwWzEiUt4HSZxJapT49dl4tJqCXEKPcRYZ8hhasbSLjQzbvn37oFQq\nceGFFw7qeY754wlumj9edPQInr/rJtw+PgkPXTwRX3/wbwiCcNrn5B8+gBfvvxWzzx2Fe84/Cwue\n/ys62tp6XHNg51bMve063DkxFQ9cNAHL3n4FVoulxzXbf/kRf7lhBm6fkIRHLj0X3y+ef8p7/fLF\np/i/qy/G7eOT8MSV05D+/TcAAIO/L4J1auzMq8ZLL72E5ORkaDQajB8/Hr/++muP12hpacGcOXMQ\nHx8PrVaLqVOnYs+ePUP5cpFe2tvFxgt5BHn61WIjw/bxxx/jsccew6FDh3DkyJHTXms2n9hfe+fx\nOkSHaKHXuP4fge2tLXjpgdsQaAjDm9+swf1/fxk/LF6Anz5d2O9z6qsq8eL9tyE8Jh5vrFiN5xd9\ngZK8o3j/2Tnd1xQezcarD8/G+Asvwdvf/YYn5y3Ano2/4fN5r3Zfs3fLRrz71OO4/Ja78M5Pm/Dw\nP17H6iUfY83nn3Rf8+vypfh83qu4+bE/453Vm3DrE0/h45fmYs/G38AYQ2qUP7Z9/l8s+PBDvPfe\ne8jOzsYf//hH3Hjjjdi3b1/36zz44INYu3Ytli5dioMHD+KKK67AjBkzUFZW5uSvqPxcfbXYeCGP\nICdO0dHRgS+//BIPPfQQZs2ahcWLF3c/VlhYCMYYli9fjksvvRRarRYLF4rBuXnLNnz78oPY/fzV\neOjiiVj4wjNobz1xJsm+rZvw3J034O7zRuOe88fgpQduR+nx3CHXueWnb2Hq6MATb7yLuBGjMOUP\n1+DGhx7HT0s+6rdXnpG+HowxPPTP1xGdlIKUsyfgkRfexK7ffkZFUQEAYPvPPyA2JRW3PvEUIuMT\ncdZ5UzD7qefw65dL0dHaKv5Zf1iJydMvx5V33IuI2HhMSpuBmx7+X3y/6IPu997yw0rMuPlOXHTt\njYiIjce0a27A5bfche8XfQAASI30R+PBjZh532O45pprkJSUhEcffRRXX3015s2b1/3/YtWqVXjj\njTeQlpaGlJQUvPDCC0hJScGCBQuG/LUjfKIgJwO2cuVKxMfHY9y4cZg9ezY+++yzHr1uAHj22Wfx\n2GOPITs7GzfccAMOHjyIq666Eprk8/DEwu/x1/8uQuGRw/hg7pPdz+lsb8c1dz+EN77+GS9+thJ+\n/gF4/dF7YO7q6r7mT9em4c6JKf22P12b1n3tsaxMjJ58Pnw12u77JkxLQ311JarLSvr8s1m6TFCq\nVFAqT+zIqNaIp9ofyfwdAGDu6oJK3fOke7VGgy5TJ44fPnDSNb6nXFNXWYGastLTXpN3MAsWsxmp\nkf4QLGYUNPb82mq1Wmzbtk2s12KB1WqFRqPp9xoiH/L4sJM4xaJFizB79mwAwCWXXAI/Pz/8+OOP\nmDlzZvc1TzzxBGbNmtX9+7lz52LstCtRO/EmTJ40Flq1Eg+/8DqeuvEKNNXVIjDUgCl/uKbH+/zv\na+9g9uQRyDu4D6MnnQ8A+PvCz2G19H8UmtLnxJTGxppqhEZE9ng8MNTQ/Vh4TNwpzx97wTQsefNF\nfPvRf3HdvY/A1NGOz+e91v0cAJgw7RKsXvoR0n9YiYuuuQGNdTX4Zv47AICGmir7NWn49LV/IGtb\nOsZNvRiVRQX40T6k01BTBWNMLCZMS8OGVctxwRVXI3nseBw/dAAbVn4Ji9mMloZ6BBvDETxyMjav\nWoKj/3cHUlNTsWHDBnz77bewWsVtbv39/TFlyhS88sorGDt2LCIiIrB8+XLs3LkTKSkp/X6diHei\nICcDkpeXh+3bt2P5cvGsbcYY7rzzTixatKhHkE+ePLnH8zIzM5FzLBcs/Wc8+J74D0DHEENlcSEC\nQw2oLC7E8v+8hdwD+9BcXwdBsMFms6G2vAyYJL6OMXqQJ8wz1uO3jhEV1ut+h7jUkXji9Xex5M0X\nsfw/b0GhUOLq2fcjyBAGhf1c0QnT0nD30//A4pf/jg/m/hkqtRqzHp2DIxm7obD35C+/5U5UlRTi\nzcfvh8Vihp/eH9fMfgAr3p/Xfc2sx+agsbYac2+/HoIgICg0DGk33IzvF83vvmb6g3/DL/99AWPG\njBHnlycn47777sOnn37aXfOyZctw//33IyYmBkqlEhMnTsTtt9+OvXv3Du5rRbgnjyBPulfqCri3\naNEiWK1WxMWd6M06Armk5MRwhU7Xc9WmxWqFbtwVuGb2g5g2xtjjsZDwCADA64/eg5DwCDzy4psI\nDY+EQqnEnGvTegzb/OnaNNSWl/ZbnyEqBv9ZnQ4ACAozorG25zL35vpaAECgof8tYi+67iZcdN1N\naKytga/WD4wxrF7yEYwn9eCvv+8RXHfvw2ioroIuMBA1ZaX44t+vd/fyGWOY/dRzuOPPz6KxthoB\nwaE4uEsc6jBGi6ch+mq0ePy1d/DIi2+hqa4GQWHhWPf159Dq9PAPFqccjhsZhz03Poe3Z47GlChf\nREVF4ZlnnkFiYmJ3LcnJydi8eTPa2trQ3NyMyMhI3HrrrT2uIUNz771SVzA4FOTkjCwWC5YuXYrX\nX38d1/bacX/27Nn49NNPcffdd/f53KjkMSg6mIvzJp2FyD6OdmtpqEfp8Vw8+PxrOPsCcUpj/uED\np0zpG8zQyogJk/D526+iy9QJta84hrx/+xaEGCO6w/R0guxhv2HVcqh8fTF+6sU9HmeMdf8Q2vrz\ndzBERiFxzNk961EqERouDu9s+/l7jJwwqXt4x8FHpUJoRBQA8YPUSWkzoFCIvf/kCH8oGLCjuAWz\nzk2C2WzGqlWrcEsfK1V0Oh10Oh0aGhqwdu1avPXWW2f8M5LToyD3RJ1ibwwaw+mvI336+eefUVtb\ni4ceegihoaE9HrvtttuwYMEC3HXXXX0+N+myO7D5t7uw5v2X8IfbZkOr06MsPw8Zm9bhjy+9BV1g\nEAKCQ7D+my9giIxCfVUlPvvXy1D69PzWHMzQykXX3ohvPvg33n92Dmb9cQ7KC/Px3cfv45bHn+we\nWtm97hd8/u/X8MKSr7sDd83nn2DUOZOh8dNh/44t+OxfL+OuJ+dCFxDY/drfL56Pc6ZNB1MosHvd\nGnz/8Qd48p0Puz8kbW6ow45fV2PseVNh7jJh47crsPPX1Xhp2aru1ygvOI5jB/ZhxPiJaGtuwk9L\nFqI49yieeOM/3dcUH9kPbUUmtmc2YUtYE1588UXYbDY8/fTT3desXbsWNpsNo0aNQl5eHv76179i\n5MiRuO+++wb8tSJ9q7VHhoGTyJBHkG+zf/g2I13SMni1ePFiTJ8+/ZQQB4Cbb74ZzzzzDNavX3/K\nY4IgIKcrBBc++T5qd3yBf8yeCZvNivCYeJx3+ZUAAIVCgSff+RCLX30ef77uUkTEJ+Cep/+Jt//0\n4JDr1fkH4B+Lv8LHL8/F07Ougi4wUBwSue+R7mvaWppRXnAcVvOJnn/ewSys+O88dLa3ITopBY+8\n+BbS/mdWj9fet2UTVn34HixdXYgfNQZ/++BTTLz40h7XbP5+JZb962UIgoAREybhxc9WInXcOd2P\n22w2rF6yEGUFx+Hjo8JZ50/Fa8t/gDHmxL8WzCYTytcvQWNlKW5YGIBrr7kay5YtQ1BQUPc1TU1N\nePbZZ1FaWoqQkBDMnDkTr776KlQq2stmuByf16enS1rGgLEzrXZzhsmTJwsZGRkuf59+OQ6VoCB3\nq4LaNkx/Ox03T4nFlFF0fNlg5Ve24v1fjuGfN4zBfRfQuLc7OQ6VkDrIGWOZgiBMPtN1NI+cuMw2\n+77aqbS/ypDEG3XQqJXYlEP7k5PToyAnLrM1txahejVCA3zPfDE5hVLBMDLKH1lFjbDaaFtb0j8K\ncuISXRYbth+vxYiogH7nbpMzGx0TiOYOM3YW1EtdCvFg8viwM/VRqSuQnT2F9WgzWTE2LvDMF5N+\njYoWj8Vbe6QS05I5mULhBR7lLDLkEeTxt0pdgexsOFINlZLRQcvDFOCnQnSIFjvz6qQuRVZu5Swy\n5DG00lYiNuI2G3OqkBLpD1+V8swXk9MaHROI/KpWVDZ3Sl2KbJSUiI0X8gjynbPFRtyioLYNhXXt\nOCuWhlWcYXRMAGwCsOZwhdSlyMbs2WLjhTyCnLjVxhxxt0DH+C4ZnrgwHbRqJTYerZa6FOKhKMiJ\n023KqUZkkAYh/jTt0BmUCoZR0QHIKmyEyULTEMmpKMiJU7WaLNhVUIfRsdQbd6axcYFo7bQgPZd6\n5eRUFOTixca4AAAPhklEQVTEqbYcq4HFKuCsGBofd6ZRMYFQKhhWH6RxcnIqeUw/HPUXqSuQjV8P\nVcJf44O4PrasJUOnVSuRGumPnbl1sNls3dvdEtf4C2eRIY8gj7lO6gpkwWSxYkNOFc6OD4JSQas5\nnW1sXCBW7izBnpJGnB8fInU5Xu06ziJDHj/Wm4+KjbjU9rxatJmsmBAfLHUpXmlsXBAYgB/3l0ld\nitc7elRsvJBHj/x3+z7UtI2tS/16qBJatRJJkTSs4goBfirEh+mw9VgtBEGgPWxc6BF7ZEi9je1A\nyaNHTlzOYrVhXXYVxsQGwEdJ31auMjY+EMW17cipapG6FOJB6G8ccYrfC+rR0G7GhAQaVnGls+PF\nE4JW7uv/IGoiPxTkxCl+OVQJtY8CqVE0f9yVwgI0iAn1w/rsKrjjdC/CBwpyMmwWqw0/H6zA6JgA\nqH3oW8rVJiYFo6imHYcqmqUuhXgIeXzYOfY5qSvwatvyalHf1oWbLog988Vk2MYnBOPHPWVYkVGC\ns6+nhVeu8BxnkSGP7lPEDLERl/ghqxx+vkqMiKazOd0hWK9GUrgeG47Q8IqrzJghNl7II8gbssRG\nnK69y4K1hysxPj6IZqu40TlJwaho6MTvRXQEnCtkZYmNF/L4m5c5R2zE6dZlV6G9y4rJybTS0J3G\nJwRBwYAVmTR7xRXmzBEbL+QR5MRlfsgqR4hejfhwWgTkTnqNCiOiArApuwpmK21tK3cU5GTI6lpN\n2HKsBhMSg6GgVYZud15qKBrazFh9iHZElDsKcjJk3+4tg8Um4NwUGlaRwti4QOh8lVj+e7HUpRCJ\nUZCTIREEAcv3FCPRqEN4kFbqcmTJR6nAxOQQZBY0oKKJDmaWM3kE+fjXxEacJqOoAfk1bbhghEHq\nUmTt/FQDrDYBn+0ulLoUr/Laa2LjhTyCPGyq2IjTLP+9GFq1EuMSgqQuRdaiQrSINfjhx6xy2Gz0\noaezTJ0qNl7II8hrdoiNOEVThxlrDlbgnKRg+KqUUpcje+enhqKsvgPpuTVSl+I1duwQGy/kEeT7\n54qNOMU3GSXoNNswhYZVPMLE5BBoVAos2lYgdSleY+5csfFCHkFOnMZqE7BkRyGSw/WIDvWTuhwC\nQKNS4vwRBuzKq0NRfZvU5RAJUJCTQVmXXYnShg6knWWUuhRykmmjwyAIwIIt+VKXQiRAQU4G5ZNt\nhTD4+2J0LO2650lC/X0xJjYQq7PK0W6ySF0OcTMKcjJgB0ub8HthPaaNNkChoJWcnuaiMWFo7bTg\nk100Vi438tiPfNK7UlfgFean50GrVuLclFCpSyF9SI30R6zBD0u3F+GRC5Og8qEZRUP1LmeRIY8e\nefAEsZEhO1rZgl8OVWLa6DBofeXx8583jDFcPi4CNc0mfL6nROpyuDZhgth4IY8gr1wvNjJk72/K\ng69KgUvG0IecnmxMXCAigjT4aPNxWGlXxCFbv15svJBHkB96RWxkSPKqW7H6QDmmjQqDn4Z6455M\nwRguHx+BisZOfJVJvfKheuUVsfFCHkFOhuXd9cegVipwCU055ML4hGCEB2nwnw256LJYpS6HuAEF\nOTmtvcUNWH2gApeMNUKvVUldDhkAhYLh2klRqG4yYeFWmlcuBxTkpF+CIOC1n48gUKvC9LHhUpdD\nBmFMbCCSwvX4aHM+mjvNUpdDXIyCnPRr7eFKZBQ14MpzImlzLM4wxnDdudFo6bTgjbU5UpdDXEwe\nn1ydt1DqCrjT3mXBy6uPICpYi8mpNG+cR/FhOkxMCsaK3SWYfV48RkcGSF0SNxZyFhny6JEHjBQb\nGbB31h1DWWMHZk2JhZJWcXLrf86LgdpHgadW7YcgCFKXw42RI8XGC3kEeelPYiMDcqisCYu3FWDq\nSAMSwvVSl0OGwV+rwjWTonC4tBmf7S6Suhxu/PST2HghjyDPmSc2ckYmixVPrzwAf40YAIR/F4w0\nICFMhzfW5KCkvl3qcrgwb57YeCGPICcD9tavR5Fd0YybL4yjpfheQsEYbr84HlabgEe/3AubjYZY\nvA0FOemWfrQai7cV4KLRYRhD29R6lbAADW48PwaHSpswb/0xqcshTkZBTgAApQ3tePLr/YgO0eKa\nydFSl0Nc4LzUUIxPCMKCTXnYeLRa6nKIE1GQE7SaLHhgaQY6zVbcPT0Rah/6tvBGjDHcOi0exkAN\n/vfLvcivaZW6JOIk8vgbO2WZ2MgpLFYb5ny1D7lVLbhnehLCAjRSl0RcSKNS4oEZyQCA2Z/8jvpW\nk8QVeaZly8TGC3kEuS5WbKQHq03AX77Zj/VHqjHzglikRvlLXRJxg1B/X9x7aRKqmztxy8e7aAl/\nH2JjxcYLeQR50QqxkW5Wm4BnVh3AD1nluHZSFKaMCpO6JOJGKRH+uGd6EvKrW3E7hfkpVqwQGy/k\nEeS5C8RGAACdZise/2IvvsksxZUTInDpuAipSyISGBMbiNlpiThS3owbPtiO6uZOqUvyGAsWiI0X\n8ghy0q2mxYQ7Pt6FtYcrceP5MbjiHFr0I2fjE4Lx0OUpKG3owDXvb8Ph8iapSyJDQEEuI9vzanH1\ne1txqKwZ916aiIvo2DYCYGR0AB6/KhVdFhtu+GA7vthdRPuycIaCXAbaTBa8sjobdy3eDR8fhj9f\nNxJnxwdLXRbxILEGHZ68fhQSw/X4+3eHcN+SPahsoqEWXtAabC9mswn4+WAFXl9zBOVNnZg60oDr\nz42GmvYWJ33w16rw8OUp2Hy4Gr/uK8f0een402UpuHdqIjT0PePRmDv+CTV58mQhIyPD5e/Tr85a\n8VZjkK4GN7LaBKzLrsS763ORU9mC6BAtbp4Sizgj7WRIBqa22YRvd5Ugp6wZBn81/u/SVMyaFAM/\ntTz6frX2yDBIHBmMsUxBECaf8TpZBLlMlDd24Lt9ZfhidxHKGzsRHuiLKyZEYXxCEBS0pzgZgryK\nFqzZW47C6jbofX1w67mxmDUpBqMi/MEYfU+5GgX5yfKXiLdJ90pXgwtYbQKyy5ux/Xgtfj1UiayS\nRgDAiCh/XDhK3PiKDoUgwyUIAgqq27D9SDX2FzbCJgDxoX64dlwkpqWE4Zy4IK8belmyRLy9914p\nq6Ag72l9mng7I126Goapy2JDQW0bciqbcayqBYfLm5FR2IBWkwUAEGfww7j4IIyLD4IhkJbZE9do\n6TDjUFEjsgobkVfZAkEA1EoFxsUGYlx0EEZH+mN0ZAASDTroON4GOS1NvE1Pl7KKgQc5v19pzgiC\nAJPFhi6rDSazeNtlscFksaLNZEVzpxnNHWY0dZjR1G5GY4cZlU2dKG/qQHlDB6pbTHD8yFUqGMID\nNTg7IQgjIv0Rb9QhSKeW9M9H5MFfq8KUUWGYMioMHV1W5Fe1Ir+iBcerWrFsVyHM1hMdwyA/FaKD\ntIgN8UOY3hfBfioE69QI9lMjyE8Fna8PtColNColtGoltCqxqX0UUDDQ0M0geHSQbz5Wg5dXZ3fP\naRW6/9N90+Mxofsx+332378X1ghAwOOvb+jzeQN6zV7/chGEvq8/8diJGsxWW49v8IHw9VEg0E+F\nIL0a8RF6TEgOhjFQg4ggDcICNfBR0sxRIi2tWomzYgNxln3veptNQG2LCeX1HahvMaG+tQv1rV3Y\nW9KA1k4LOkxWDOZvgZIxKBUMCoXYefFRiAGvVIj3M8bgiHpH5p+456T7Tvp54Hi8r58RJ79GVskE\nMAAz/r1/EBX37bZzY/HgRUnDfp3T8egg1/v6YGS4fSMnduLG8ZP61P+Jpz4GBgRYfGAVBMSE+dmv\n7/l/kfV6/ZN/xVjP1+p9Te9viJOvdXzT+CiZvSngo2BQKRVQ2u9TKRVQKRXw8/WBn6+9Z6JW9ghq\nJfVMCAeUSobIIC0ig7Q97rcJYrfGZhPQ3mVFu8mCtk4Luiw2mC0n/nVqtoi3FpsAQRBgswmwCThx\n232f+PvuDlQfPx16d7L6euzknyon7hN/5aNkUDB2In+GwaD3HfZrnIlHB/mk+GBMcsbClfXitLtv\nZkwZ/msRQrxe2lfi7Qd3TpS2kAHy6CB3mrQ1UldACOHIGs4iQx5B7uMndQWEEI74cRYZ8vjE7Nh8\nsRFCyADMny82XsgjyIu/FhshhAzA11+LjRfyCHJCCPFiFOSEEMI5CnJCCOEcBTkhhHDOLZtmMcZq\nABS5/I1OzwCgVuIaCCH88ITMiBcEIexMF7klyD0BYyxjILuIEUIIwFdm0NAKIYRwjoKcEEI4J6cg\n/0jqAgghXOEmM2QzRk4IId5KTj1yQgjxSlwHOWPsE8ZYNWPs0En3hTDG1jHGcu23wb2ecy5jzMoY\nm3XSfVbGWJa9/ejOPwMhxD0GkxeMsb+elAmH7BkRYn+skDF20P6YhIcRn8B1kANYAuDKXvc9A2CD\nIAipADbYfw8AYIwpAbwJYG2v53QIgjDB3q53Yb2EEOkswQDzQhCEfzkyAcCzADYLglB/0vOm2x/3\niOmJXAe5IAhbANT3uvt/ACy1/3opgBtOeuwJAKsAVLu+OkKIJxlCXjjcDmC5C0sbNq6DvB/hgiBU\nAID91ggAjLFoADcC+LCP52gYYxmMsV2Msb7+RxJCvFOfeeHAGPOD2ItfddLdAoDfGGOZjLGH3Vbp\nacjjhCDRuwD+JgiCtffhywDiBEEoZ4wlAdjIGDsoCMJx95dICPEw1wHY3mtY5UJ7XhgBrGOM5dh7\n+5Lxxh55FWMsEgDst45hlMkAvmKMFQKYBWC+o/ctCEK5/TYfQDqAc9xcMyFEGv3lhcNt6DWsclJe\nVAP4DsB5bqjztLwxyH8EcI/91/cA+AEABEFIFAQhQRCEBAArATwmCML3jLFgxpgvADDGDAAuBJDt\n/rIJIRLoMy8AgDEWCOCSXvfpGGP+jl8DuAJA9ywYqXA9tMIYWw4gDYCBMVYK4J8A3gDwNWPsAQDF\nAG4+w8uMBrCQMWaD+IPtDUEQKMgJ8TJDyIsbAfwmCELbSfeFA/jOPjzrA+BLQRB+dUP5p0UrOwkh\nhHPeOLRCCCGyQkFOCCGcoyAnhBDOUZATQgjnKMgJIYRzFOSEEMI5CnJCCOEcBTmRJcbYJsbY5fZf\nv8IYe0/qmggZKq5XdhIyDP8E8JJ946NzANA+9IRbtLKTyBZjbDMAPYA0QRBapK6HkKGioRUiS4yx\nswFEAjBRiBPeUZAT2bFvV/oFxNNh2hhjf5C4JEKGhYKcyIr9xJdvAfxFEIQjAF4G8IKkRREyTDRG\nTgghnKMeOSGEcI6CnBBCOEdBTgghnKMgJ4QQzlGQE0II5yjICSGEcxTkhBDCOQpyQgjh3P8DRzFP\nQDfA9sgAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "sample = 1045\n", "mu = 1060\n", "w = mu - sample\n", "x = np.linspace(mu - w - 2, w + mu + 2, 1000)\n", "y = ss.norm.pdf(x, loc=mu, scale=3)\n", "plt.plot(x,y)\n", "plt.fill_between(x, y, where= np.abs(x - mu) < w, color='lightblue')\n", "plt.text(mu,np.max(y) / 3, 'Area={}'.format(np.round(interval_p,6)), fontdict={'size':14}, horizontalalignment='center')\n", "plt.axvline(mu - w, linestyle='--', color='orange')\n", "plt.axvline(mu + w, linestyle='--', color='blue')\n", "plt.xticks([mu - w, mu + w], [sample, mu + w])\n", "plt.yticks([])\n", "plt.ylabel(r'$p(x)$')\n", "plt.xlabel(r'$x$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We would expect to see values outside or at the boundary of this interval 0.00000001% of the time." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "What is signficiant?\n", "----\n", "\n", "Would we call it significant if it was 0.1%? what about 1%? It turns out the convention is 5%. This is called our $\\alpha$ or significance level." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We saw in the cash4gold example how to compare a single sample with a known population. What about when we don't know the standard deviation of the population?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "I open up a cash4gold store and people bring me their jewlery. I know the probability distribution for gold melting is normal with mean 1060 $^{\\circ}$C. I do not know the standard deviation. If I get a sample that melts at 1045$^\\circ{}$C, should I be suspicious?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "I don't know. We have no way of estimating the standard deviation with one sample, so we can't say anything. If we have at least 3 samples though, we can compute the sample standard deviation" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "This leads us to the beautiful part about hypothesis testing: since we're assuming the null hypothesis, that the samples are gold, we can use what we learned about normal distributions to estimate the population standard deviation from our samples.\n", "\n", "----" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "General Strategy\n", "----\n", "\n", "The general approach for hypothesis tests is as follows:\n", "\n", "1. Start with a null hypothesis, $H_0$, that corresponds to a probability distribution/mass function\n", "2. Compute the probability of an interval which *inculdes* values as extreme as your sample data. Note this may be single- or double-sided depending on if \"extreme\" means both above and below or only above/below the mean. This is your $p$-value.\n", "3. Reject the null hypothesis if the $p$-value is below your significance level ($\\alpha$), which is generally 0.05" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Student's $t$-Test\n", "====" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "I open up a cash4gold store and people bring me their jewlery. I know the probability distribution for gold melting is normal with mean 1060 $^{\\circ}$C. I do not know the standard deviation. Someone brings in 6 samples and they melt at 1035, 1050, 1020, 1055, and 1046 $^{\\circ}$C. Should I reject the null hpyothesis, that these are gold?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "#### Logical Steps: \n", "\n", "1. Assuming the null hypothesis, compute our uncertainty in the true mean confidence interval. This is our probability distribution\n", "2. We happen to have the true mean, so then we see how big the confidence interval has to be to include it. This is us constructing our probability interval\n", "3. Take the area of that interval, our $p$-value, and compare with the significance level" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX5wPHP996bvRfZIQTZYSggy4EbReRXZ1WstI5W\npdpd/aktv9a2WrVabd22aGtBVBQnWqYge0MChJGQTQYZZCf3nt8fJ0FGdu69596T5/3ivG5yc+45\nT4A895vveL5K0zSEEEJ4L4vRAQghhOgbSeRCCOHlJJELIYSXk0QuhBBeThK5EEJ4OUnkQgjh5SSR\nCyGEl5NELoQQXk4SuRBCeDmbO24SHR2tpaamuuNW36o+oD+GDnPvfYUQplJeXg5AVFSU2++9bdu2\nMk3TYro6zy2JPDU1la1bt7rjVt9aPl1/vHy1e+8rhDCVBQsWADB37ly331spdbQ750nXihBCeDm3\ntMgNkf6Y0REIIUzgoosuMjqELpk3kcddbnQEQggTSEtLMzqELpk3kVfs1B8jxhkbhxBO0NzcTH5+\nPg0NDUaH0u/Y7XYArFary+7h7+9PUlISPj4+vXq9eRP5tp/ojzLYKUwgPz+fkJAQUlNTUUoZHU6/\nUlZWBkB0dLRLrq9pGuXl5eTn5zNo0KBeXUMGO4XwAg0NDURFRUkSNyGlFFFRUX36bUsSuRBeQpK4\nefX139a8XStCdELTNFZnlZJZWE1ieABXj47Dz+a6PlAhXEla5KLfKa9p5OZXN/D9f27h6S8P8JN3\ndzL96dVkFlUZHZrXufvuu8nMzOzVa3NyckhPT3dyRM4zffp09y9k7CXztsjH/tHoCIQHqm1s4bbX\nN5JdXsd3LxjI2NRwckpqWbTuKDe9soGP513A4Jhgo8P0Gm+88YbRIbhcaGio0SF0ybwt8pip+iHE\nKR5fupeDJTXcddlgzh8ShZ+PlWGJoTxw9VA04PtvbaGpxW50mB6ntraWmTNnMnbsWNLT03n33XeB\n01utwcHBPProo4wdO5bJkydz7NgxAA4fPszkyZOZOHEiv/nNbwgOPvuN0m6388tf/pKJEycyZswY\nXn311W7H8Lvf/Y6JEyeSnp7Ovffei6ZpJ2P76U9/ykUXXcSIESPYsmUL119/PUOGDOGxx/QFgzk5\nOQwfPpw777yTMWPGcOONN1JXV3fafX19fVm9ejVTpkzhvPPO46abbqKmpsZJf7POYd4Weel6/VGS\nuWi18Ug5S7YXcMWYWIYkhJz2tehQP26ZNpB/rjzCs8uzeGTGCIOi7Nr/fZJBZmG1U685MiGU384a\n1eHXly1bRkJCAp999hkAVVVnd0PV1tYyefJk/vCHP/CrX/2K119/nccee4yHHnqIhx56iFtvvZVX\nXnml3eu/+eabhIWFsWXLFhobG5k2bRpXXnnladPxOoph3rx5/OY3vwHgjjvu4NNPP2XWrFmAnoS/\n/vpr/vrXvzJ79my2bdtGZGQkgwcP5qc//SkABw4c4M0332TatGn84Ac/4KWXXuIXv/jFyfsWFhby\n+9//nuXLlxMUFMRTTz3FX/7yl5P39ATmbZHv+l/9EAJwODTmf5xBVLAvl42Nb/ec0QPDGZUcxtvr\nj1Je0+jmCD3b6NGjWb58Ob/+9a9Zu3YtYWFhZ53j6+vLtddeC8D48ePJyckBYMOGDdx0000A3Hbb\nbe1e/6uvvuLtt99m3LhxTJo0ifLycg4ePNitGFatWsWkSZMYPXo0K1euJCMj4+RrrrvuupOvHTVq\nFPHx8fj5+ZGWlkZeXh4AycnJTJs2DYA5c+awbt260+67evVqMjMzmTZtGuPGjeOtt97i6NFu1bJy\nG/O2yIU4xcr9JewvPsGci1LxtXXcfplxXjzPLt3PMyuy+NPs0W6MsPs6azm7ytChQ9m2bRuff/45\njzzyCFdeeeVZLVIfH5+T0+isVistLS3dvr6mabz44otcddVVPYrhV7/6Fffffz9bt24lOTmZ+fPn\nnzYf28/PDwCLxXLy47bP2+I7c+rfmZ9rmsbFF1/MkiVLuv39uJt5W+RCtNI0jZdWHyIqxJexgyI6\nPTcxMpDRKWF8tK2AyromN0Xo+QoLCwkMDGTOnDn84he/YPv27d1+7eTJk/nggw8AWLRoUbvnXHXV\nVbz88ss0NzcDkJWVRW1tbZcxtCXt6OhoampqeP/993v8veXm5rJhwwYAFi5cyAUXXHDa18ePH8/m\nzZs5dOgQAHV1dWRlZfX4Pq4kLXJhervzq9ieW8kNk5OxWrpeeDE9PZY9uVW8tekoD10yxA0Rer49\ne/bwy1/+EovFgo+PDy+//HK3X/v8888zZ84cnn32WWbOnNlut8zdd99NTk4O5513HpqmERMTw0cf\nfdRlDOHh4dxzzz2MHj2a1NRUJk6c2OPvbcSIEbz11lv88Ic/ZMiQIdx3332nfT06OpoXX3yRW2+9\nlcZGvcvtiSeeYOjQoT2+l6uothFeV5owYYImG0sIozyyZA8fbM/nt7eMJsC360U/mqbx9Ef78LVZ\nWPPz6R6xonLfvn2MGOG5A7CdqaurIyAgAKUUixYtYuHChSxdutTosAB91sq1117L3r17OzzH1bVW\n2rT3b6yU2qZp2oSuXmveFvn4542OQHiAuqYWPt5VwNjU8G4lcdD7SKcNj+GDjXl8c6SMCwZ3udOW\n6MS2bduYN28emqYRHh7OP/7xD6ND6pH2foPwNOZN5FK+VgBfZhRT22hn8tCetabGpUXw4aY8Fm/P\nl0TeRxdeeCG7du0yOox2paamdtoaB3pdWtadzDvYWbxcP0S/9vmeYsKDfBg4IKhHrwvyszE0IZR1\nB8posTtcFJ3wBo2NjSf7xj2VeRP53if0Q/RbNY0trMkqZczAcCy96Oc+d1AEx2uaWHmw1AXRCW9x\n4sQJTpw4YXQYnTJvIhf93sr9JTS1OBiX2vmUw46kDwzHalF8uLPAyZEJ4VySyIVpfbGniLBAH1J6\n2K3SJsDXyrCEEDYdKsfukO4V4bkkkQtTampxsCarlFEpYb3qVmkzMjmM4zVNbMmtdGJ0/VNbsayc\nnBz+85//nHx+69atPPjgg52+Vkredk4SuTClbUcrqGuyMzKxbyVIRyTpU8++yChyRliCsxP5hAkT\neOGFFwyMyPuZN5Gf/6p+iH5pTVYpVosiLT6k65M7ERHsS1y4PxsOlTspMu/UVu717rvvJj09ndtv\nv53ly5czbdo0hgwZwubNmwGYP38+zzzzzMnXpaennyye1ebhhx9m7dq1jBs3jueee47Vq1efLLY1\nf/587rjjDi699FKGDBnC66+/flYs7i55++yzzxIWFtatkregFwBrr+Ttww8/zMiRIxkzZsxp1RWd\nwbzzyEOHGR2BMNDXWaUMGhCEv0/ft28bkRTGmoxjFFc3EBfq74To+uiLh6F4j3OvGTcarn6y01MO\nHTrEe++9x2uvvcbEiRP5z3/+w7p16/j444/54x//eNaS+o48+eSTPPPMM3z66aeAXl3wVLt372bj\nxo3U1tZy7rnnMnPmzNO+bkTJ27bE21XJ27KyMp544omzSt7OmzePDz/8kP3796OUorLSuV115m2R\n53+iH6LfKaluILOo+mS3SF+NSA7FocGX+4qdcj1vNWjQIEaPHo3FYmHUqFFcdtllKKUYPXr0Wa3u\nvpg9ezYBAQFER0dzySWXnGztt3F3ydvU1NSTBbO6Knm7cePGdkvehoaG4u/vz913382SJUsIDAx0\n2t8XmLlFvv9Z/TFplrFxCLf7+qBeG2NoQt+6VdqkxgRhsyq+OVTOnZNSnXLNPumi5ewqZ5aBPbVE\nbFtJWJvNhuOUGT6nlpTtru6UlXVnyVtN0zhx4gSxsbHdiu2KK65g4cKFZ8W0efNmVqxYwaJFi/jb\n3/7GypUru/ib6D7ztshFv7X2YCkhATbiIwOccj2b1UJabDB78itxuKHInDdLTU09WeJ2+/btZGdn\nn3VOSEhIpwtsli5dSkNDA+Xl5axevfqsioaeXPJ28uTJfPPNN2eVvK2pqaGqqoprrrmG559/np07\nd/b43p0xb4tc9EuaprHxSDnnxIX0adrhmc6JC+Hz7YVkl9cyOFo2Z+7IDTfccLLbY+LEie2Weh0z\nZgw2m42xY8cyd+5czj333NO+fv755zNz5kxyc3N5/PHHSUhIOK3rxpNL3sbExLBgwYKzSt6GhIQw\ne/ZsGhoa0DSN5557rsf37oyUsRWmcrS8loufXs1NU1OYMsx5ZUdzSmp44bMsHr1uBPdMTXPadbvL\nm8vY9sT8+fMJDg52+qyOvmgrY1tTU9Nlydu+6EsZW+laEaay6chxAAbF9m41Z0eSo4PwtVnYcLh/\nT0MUnsm8XStT/mV0BMIAm7KPE+xvIzbMudMErRZFWmwwGQXVODTNqd024lvz5883OoSzhIeHA3rf\nuqta431l3hZ5ULJ+iH5lU3Y5abHBLtnV55z4EI5VNnC4rMbp1+4Od3SDirPZbDZsNte2efv6b2ve\nRH70Xf0Q/UZBZT35FfWcE+eawci01u6adQZ0r/j7+1NeXi7J3AD19fXU19e77PqaplFeXo6/f+9/\nizRv18rB1s1hB95ibBzCbTZn6wk2LdY1iTwpKhCrRbEl5zjfn5zqknt0eO+kJPLz8yktldro7ta2\nxL6t6Jcr+Pv7k5SU1OvXmzeRi35na04FAb5W4iKcM3/8TDarheSoQPYVVqNpmls3Zfbx8TltCbpw\nnwULFgAwd+5cQ+PojHm7VkS/szOvkuToQCwW1yXY1AFB5JXXUVHf7LJ7CNFTksiFKdQ1tbC/6ASp\nMc6ddnim1AFBtNg1NuXINEThOSSRC1PYk1+FXdN6vMlyT6UO0PtJN2Yfd+l9hOgJ8/aRX9DzOgrC\ne+3I08uCJke7NpGHBvoQGezLnvwql95HeI6bb77Z6BC6ZN5E7u+85dnC8+3MrSQm1I9gf9f/l04d\nEMTh4hqaWuz42vpe71x4NmeXnHUF83atHFmgH8L0NE1je24FKdHu+YEbNCCYqrpmMo91XMFPmMfO\nnTudXq3Q2SSRC69XVNVAyYlGBrlo/viZUloHVDfnSD95fyCJXAg32JHb1j/unhZ5fIQ/VotiZ55z\nt+sSorckkQuvtzOvAh+rIt5FC4HOZLNaSIgM4OCxGlkyLzyCJHLh9XbkVpIUFYjN6r7/zinRgeSV\n1VLbZHfbPYXoiCRy4dVa7A72Flad7Ld2l+ToIBqaHewuku4VYTzzTj+c/rnREQg3yC6rpaHZQXKU\ne6eIJbXeb2tOJVNTZaqrmd1+++1Gh9Al8yZym+fP/RR9l1lUDeCyQlkdiQ33x9dmYXeBLAwyOx8f\nH6ND6JJ5u1ayXtIPYWqZhdXYrIrYcOfuCNQVq0WRGBnA4WMnZMDT5LZs2cKWLVuMDqNT5k3kuYv1\nQ5haRmE18eEBWF1Y8bAjydFB5JfXUdUglRDNLCMjg4yMDKPD6JR5E7kwPU3TyCisIjHKvd0qbZKj\nA2m2a+yU7hVhMEnkwmsdq26koq6ZpEhjxkPaFiBty5UVnsJYksiF18os0lvCcZHGtMijQ/3w97Gw\nJ7/akPsL0UYSufBamYV6Ak0wKJFblCIxKpDsUlnhKYxl3umHl682OgLhYnsLq1tbxcaVkk2MDGTD\ngVJONLYQ6u/509REz3nyXp1tpEUuvFZGYRVJBrXG2yRGBdBs19hTJAOewjjmTeT7ntEPYUonGprJ\nO15PoptXdJ4psfWNZHuuLNU3q/Xr17N+/Xqjw+iUeRN5waf6IUxpf7G+qYO7Kh52JDY8AJtVkVEo\nA55mlZWVRVZWltFhdMq8iVyYWkbr3G2j5pC3sVoUceEBZJfWGBqH6N8kkQuvlFlUTbC/jdAA4wcY\nEyMDyC+vo665xehQRD8liVx4pb2F1SRGBqCU+5fmnykxKpDaRjsHS6RVLoxh3kRuDdAPYTrNdgcH\nj504WUrWaG0DnjvyZcDTjHx8fDy+AqJ555Ff8oXREQgXOVRSQ7NdM2wh0JkSIgNQwB6puWJK3lCP\n3LwtcmFaRq/oPJOfj5XoUD8OHZOuFWEM8ybyPb/XD2E6mUXV+FgtxIS6twZ5ZxKjAskrr6PZ4TA6\nFOFka9asYc2aNUaH0SnzJvJjK/RDmE5GYRXxEf5YDKhB3pHEyACO1zSRX1FvdCjCybKzs8nOzjY6\njE6ZN5ELU9JrkFefLCHrKdpWmG7LqzA4EtEfSSIXXqWgsp4TDS0e0z/epm3miuzhKYwgiVx4lbaB\nzvgIz2qRhwT4EBrgQ1Zr6QAh3Mm80w/9ooyOQLhAZlE1CoiP8JyBzjaJUQEcLavFoWlYPGChknCO\nwEDPajS0x7yJ/MIPjI5AuEBGYTUxYX74GViDvCOJkYEcKKimvLaRmGDPe6MRvXPzzTcbHUKXpGtF\neBW9BrlntpASowJwaLBLKiEKNzNvIt/5iH4I06iqa6awssHwGuQdaRuA3SVL9U1l+fLlLF++3Ogw\nOmXerpWyDUZHIJwss6h1oNPDZqy0iQrxw89mYV+RtMjNJD8/3+gQumTeFrkwnYzC1hrkHprILUqR\nEBlAdkmt0aGIfkYSufAamUXVhAb6EOIBNcg7khAZQP7xOuqapDa5cB9J5MJrZBRWk2jw1m5dSYgM\npLHZwf4SmU8u3Me8iTwwST+EKTS22DlUUkOShy3NP1Nbt8/OPBnwNIvQ0FBCQ0ONDqNT5h3snPpv\noyMQTnTwWA12h0aCh7fI4yICUAr2FspSfbO4/vrrjQ6hS+ZtkQtTObk030MHOtv42iwMCPPnkAx4\nCjcybyLf9hP9EKaQWVSNn81CdKif0aF0KTEygLyyOuwOzehQhBMsW7aMZcuWGR1Gp8zbtVKx0+gI\nhBNlFFYRHxngFTVMEiID2X6kgvzKegZ66CpU0X3FxcVGh9Al87bIhWk4HBqZRdUeuzT/TG0Dntvz\npTa5cA9J5MLj5VXUUdtoJynKs/vH27Qt1d+TLwOewj0kkQuP1zbQGefhM1batNUm318sS/WFe5i3\njzxkqNERCCfJLKrGorwnkYPeKs8prUPTNJQX9OuLjkVFef7eBuZN5JNeMzoC4SSZhdXEhvnja/Oe\nXyATowLIKqymurGFMH/PLSkgujZr1iyjQ+iS9/xkiH5rb2EVCV7SP94mMTIQhwY7paStcAPzJvJN\n9+qH8GrlNY0cq24kyUNrkHfkZG3yAknk3u6TTz7hk08+MTqMTpm3a+VEltERCCfYV6QXn/L0pfln\nig7xw9dmOTlQK7xXeXm50SF0ybwtcmEKmUX6FL54L5lD3sZiUSREBHC4pMboUEQ/IIlceLSMwmoi\ngnwI9ve+Xx4TIgPIK6+nxe4wOhRhcpLIhUfLLKz2+EJZHUmIDKChyc4BqU0uXMy8iTxinH4Ir9XQ\nbOdwaQ3JXjbQ2aZtk+jtUpvcq8XFxREXF2d0GJ3yvt9Xu2v880ZHIProQPEJHJo+lc8bxbfWJt9T\nIEv1vdmMGTOMDqFLPW6RK6WClFJWVwQjxKkyvKQGeUd8bRZiQv3IOiZdK8K1ukzkSimLUuo2pdRn\nSqkSYD9QpJTKUEo9rZQa4vowe2H9HP0QXiuzqIoAXysRwb5Gh9JriZGB5JbpS/WFd1qyZAlLliwx\nOoxOdadFvgoYDDwCxGmalqxp2gDgQmAj8KRSyvMyZl2+fgivlVlYTUKEd9Qg70hCZADHa5oorm4w\nOhTRS9XV1VRXe/Z6gO70kV+uaVrzmU9qmnYc+AD4QCklxSSEU9kdGvuKTzDxnEijQ+mTthWe2/Mq\nmRnmnV1EwvN12SJvS+JKqedVB2Xc2kv0QvRFTnkt9U12r1uaf6a2gVpZqi9cqSeDnTXAx0qpIACl\n1JVKqW9cE5bo77ytBnlHQgN9CAmwnSw1IIQrdHv6oaZpjymlbgNWK6UagVrgYZdF1lfRU4yOQPTB\n3sIqrBZFXLi/0aH0WUJkINmlslTfWyUlJRkdQpe6nciVUpcB96An8HjgLk3TDrgqsD4b9yejIxB9\nkFlYTXyEPzar969ZS4wMYE3GCWoaWwj2M+/SDbO6/PLLjQ6hSz35KXkUeFzTtOnAjcC7SqlLXRKV\n6Nc0TWNvQRUJXroQ6EyJkQHYHZr0kwuX6XYi1zTtUk3T1rV+vAe4GnjCVYH12dob9EN4neLqBirq\nmkmJNkcib3tD2iWbTHilxYsXs3jxYqPD6FSXv+cppZTWzmoGTdOKWrtbOjzHUI2eX0NYtC+jQB/o\nTPDSFZ1nign1w8eq2Cu1yb1SXV2d0SF0qVsLgpRSP1ZKpZz6pFLKF5iilHoLuNMl0Yl+KaOwGoVe\nq8QMLBZFfEQAh4/JgKdwje6MvMwAfgAsVEqlARVAAPqbwFfAc5qm7XRdiKK/ySisYkCYP34+5inp\nkxAZyO6cClrsDlMM4ArP0mUi1zStAXgJeKl1BWc0UK9pmnT4CZfIKKw2TbdKm8SoADZmlZFVWsPI\nuFCjwxEm0+2mgVLqamAtsBp4TSk12VVBOUXsZfohvEplXRMFlfVev6LzTImtb0w7pDa51xk0aBCD\nBg0yOoxO9WRS60vAHCATGA88o5T6u6ZpC10SWV+NftzoCEQvtK3oNFsij48IQAG7Cyq5fWJKl+cL\nz3HxxRcbHUKXepLIj2ma1rYkf7lSagOwCfDMRC680t7Cts2WzdW14udjJTrUj6xiGfAUzteTUZcc\npdQTrbNVAJoBzy0gsepq/RBeRd9s2dflmy1bWhoIqC0g6EQugTV5WFtcP8UsMTKA3LJaqU3uZd55\n5x3eeecdo8PoVE9+WjTgeuAepdRBIAV4Ryk1RNO0gy6Jri/s9UZHIHrBJQOdmkZ4RQYDir4hpmQL\noVUHCKgvPeu0Jp8QqiKGUxE5huKECymPmYBmcd4bSkJUIDtzKimsricxzFxdR2bW3Oz5xV17UjTr\nVgCllD+QDoxtPd5QSqVpmpbsmhBFf1HfZOdIaQ1XjI13yvV8GysYdHARKTkfEXLiKABVYUMoibuA\n2uBk6gMG4LD4oDQ7/g3lBNYWEl6RweCstxm6/00a/SLITb2OI0Nupzak7/3aCRFttcmrJJELp+px\nc6N1OuLW1kMIp9lXXI1Dg+SovrXI/erLGJ7xMgOPfIDN3kDpgPM5OOJuChMvo8m/640qrC11xBZ+\nTVLuMgZnvcM5B96mIPlKMsc8RE1oWq/jSmz9vnblVzIr3TlvVkJALxK5EK7y7WbLvWutWuxNnHPg\nLYZlvILV3khu6nUcHPEDToSd06Pr2G2BFKbMoDBlBv51x0g7tJBzDrxNQv5ycgbfRMbYn9LsG9bj\n+EIDfAj2t7G/WJbqC+cybyJPvNboCEQPZRZWEeRnJTyo5zsHRpTvYfzGRwitPkRR4iXsGffLPrWe\n2zQExpI55iccHjKH4RkvM+jQIhLyV7BzwmMUJl/Vo2sppUiIDCCnVB/w7GDDLeFhhg4danQIXTJv\nIh/xC6MjED20K6+KpKjAHiU45Whh+N6XGJb5Ko3+0Xxz8ascS3D+vN/GgGh2TXicnME3Mn7To0xe\n9xC5A2exY+J87D5B3b5OUlQgq/ce43hdE1FBfk6PUzjf1KlTjQ6hS1L0QXiEhmY7WcdOkBLd/aTo\n11DOtNV3MyLjJfIGXsvyaz5xSRI/VVXECFZduZjM0Q+SnPsZl355I6GV3d9fJTk6EIcGW3IrXBil\n6G/Mm8iXT9cP4RUyi6ppcWgMjOleIo8o28Wly64nqmwHWyf9kW1TnqLZ1z01TDSLjf3p97P2kn/i\n01zDJV/dTOLRz7v12rY3qu2SyL3GggULWLBggdFhdMq8iVx4ld2tNUgSujFjJSHvv1y08ns4LD6s\nvmIRuWnXuzq8dpXFTmLFjA+piExn0vqfMSzjFehisU94kD7gubegyk1Riv5AErnwCLvzqwgL9CEs\nsPOBzsEH3mbSugepDB/OqisXUxUxwk0Rtq8xIJp1l/yT3NTrGLX7ecZvegTl6HgBiVKKlOhADpfI\nCk/hPOYd7BReZXdBFwOdmkb6zmcYuv9NCpKuYOuUP2O3eUY9FofVl62Tn6ImZCAj97yIT9MJNk/7\nCw5r+4OZydFB7Muv5lhNI3Eh/m6OVpiRtMiF4WoaWzhcWtNx/7imMXbbEwzd/yaHh9zGpmnPe0wS\nP0kp9qc/wM7xj5NQsIIpa+7rsH5LcnQgGrApW7YjFM5h3kSecrN+CI+3J78KTYPk9krXag7GbZ3P\n4IPvkDX8++wa/zhYPHfnoCNDb2frpD8xoGQjF6y6C1vT2XXlUmL073NbrtQm9wajRo1i1KhRRofR\nKfN2rQy93+gIRDftbt1dPjH6jETusHPe5sdJzV7CgZH3kjHmp+AFi2hy075Di08g53/zc6atuZd1\n0984ba55sL8PEcG+J2uvC882ceJEo0Poknlb5C11+iE83u78KqJCzihdq2mcu3U+qdlL2Jf+gNck\n8TaFyVexedqzRJTvZurX92FtOb0aZ0p0INklNThkwNPjNTc3e3wFRPMm8tXX6IfweLvzK0/fEUjT\nSN/5ZwYdfo/9o37EvtE/9qok3qYw+Sq2Tn6S6JItTF47D4u98eTXkqMDKa9p4uhxaWx4Om+oR27e\nRC68wvHaJvIq6k8b6ByW+SpD9/+Tw0NuJ3P0QwZG13f5qbPYPukJYou/YdI3Pzk5NTG5dWHQ5qPH\njQxPmIQkcmGotv7x5Nb+8bSsfzNq9/McTZ3NrvGPemVL/ExH025gx4TfEl+wivM2Pw6aRnJUIArY\nJis8hROYd7BTeIWdeZUoIDEykOTsjxm37QkKky5n+6Q/gDJPOyN7yK34NZQzcu/faPCPIWPcz4kJ\n82OfDHgKJ5BELgy17WgF8ZEBJB/fwvhN/0tJ7GQ2T33WqVuseYr96Q/g31DKsH2v0xAQQ0r0Rewv\nqKbJbsfX6rlTKoXnM99PS5u0uUZHILrgcGjszKtkRnwFk9b9ghOhg9h4wYsdroj0ekqxc/xv8Gs4\nzpjtf2J2qi9bG4aSUVzNuYkRRkcnOjBu3DijQ+iSeX53PVPaXEnmHu5QaQ0BDaU8Wv4b7DZ/1l/8\nGi2+IUa1z+DbAAAb+ElEQVSH5VoWK1umPkN5zHhuO/oEUy17WX9YVnh6snHjxnl8MjdvIm8o0w/h\nsXYdLuAfvk8TaD/B+otepT6of+xj6bD6seGil6gJHcSrPs9RdHCH0SGJTtTV1VFX59nTRM2byNfd\nqB/CM9lbSF//EMMtuWye9hxVkSONjsitmn1DWX/xqzRZ/JmX/wjaiWKjQxIdWLx4MYsXLzY6jE6Z\nN5ELz6Vp8PkvGFGzkZeC76ck0bW7+niq+qAEXh74J0Id1bT8+2ZoqjU6JOGlJJEL9/vmedj2T/7e\nch1ZA/v3b03WlPP4cfM8bMf2wAd3g8NudEjCC0kiF+61531YPp9jKdfyTMvNpA7o/h6dZpQSE8RK\nx3iWJj4IBz6HLx81OiThhSSRC/fJ+QY+ug8GTuM/8Q+jlIXE9krX9iMBvlbiIvx5rfFymHw/bHoZ\nNr5idFjCy5h3HvmQ+4yOQJyqNAsW3QbhA+GWf7PlnSwSIwPx85GFMKkDgtlx5DgN9/wf/hVHYdnD\nEJ4Cw6XomyeYMGGC0SF0ybwt8oG36IcwXk0JvHMjWH1gzvs0+4WzI7eS1Nj+3a3SJnVAEA3NDrYW\nVMMNr0PCOPjgLijYbnRoAkhPTyc9Pd3oMDpl3kRem6cfwlhNtbDwu3oyv+1diEhld34V9c12hsSZ\nfPFPN7WNE6w7XAa+QXDruxAYDf+5BSpzDY5OVFVVUVVVZXQYnTJvIt9wh34I4zjs+kyMgu1w45uQ\nOB74dq9KaZHrokP8CAv0YWt2ayXEkFi4/T1oaYR3boJ62RLOSB9++CEffvih0WF0yryJXBhL0/S+\n3gOfw9V/huEzT35p05HjxEf4E+zvY2CAnkMpxeC4YPYXVtNid+hPDhgOt/wLyg/Du3P0pC5EBySR\nC9fY8HfY/BpMmQeT7j35dIvdwdac4wwaEGxgcJ7nnPgQahpa2J5/Sus77WKY/XfIWQtLHwCHw7gA\nhUeTRC6cL3MpfPUYjLgOrvj9aV/KKKymtsnOkHjpHz/VOa3jBWsOlp7+hbG3wKWPw573YOXvDIhM\neAPzTj8UxsjdBEvuhaSJcP1rYDm9rbA5W9/abFCstMhPFRXiS3igz8m/n9Nc+HOoyoN1z0FYMky8\ny/0BCo9m3kQ+/OdGR9D/lB/WZ6iEJsCti8An4KxTNmWXMyDMj9BA6R8/lVKKwfEh7CvQ+8ltVsup\nX4RrnoXqIvj8FxCaCMNmGBdsPzNlyhSjQ+iSebtWkmbph3CP2jJ9rrhScPv7EBR11il2h8bm7OOk\nSWu8XYPjgqlpaGFbfjv7eFptcOM/IG4MvP99KNjm/gD7qWHDhjFs2DCjw+iUeRN59QH9EK7XVKvP\nea4u1FviUYPbPW13fiXVDS0MSwh1c4Deoa2ffHVWafsn+AXDbYshqHWO+fFsN0bXf5WVlVFW5tl7\nG5g3kW/+oX4I12ppgnfvgMLtcMObkHx+h6euO6j/MAyOlxZ5ezrtJ28TEgu3fwD2Zv03oLpOzhVO\n8emnn/Lpp58aHUanzJvIhes5HHoRrMMrYNZfYcS1nZ6+9lAZydGBMn+8A0ophiTo/eQNzZ2Us40Z\nCrcuhMo8PZk31rgvSOGRJJGL3mlb8LP3fbh8Ppz3vU5Pr2lsYfvRCobKtMNODUsMpa7RztojXfwq\nP3Aq3PRPKNypFyNrbnBPgMIjSSIXvfP1M7D5VX3Bz7SfdHn6piPltDg0hiVKIu/M0IRQFLBi37Gu\nTx4+U18wlL1GL7Jlb3F5fMIzSSIXPbf1H7DqCRjzXX3Bj1JdvmTtwTJ8bRYGxkj/eGeC/W0kRQey\n6Ug3+77H3QoznoL9n8InD8rqz37KvPPI0x8zOgJz2rkQPv0ZDLkKZv/trAU/HVl3qIy02GB8bNJ2\n6MrwxFCW7y6msLqehNCz5+KfZfKPoKESVv8J/MPgqj92681VdM9FF11kdAhdMu9PVdzl+iGcZ8/7\nsPR+GHQR3PyWXl+8Gwoq6zlUUsPwRJl22B3DE0P1IYiM4u6/6OJfw6QfwcaXYPWTrguuH0pLSyMt\nLc3oMDpl3kResVM/hHNkLtWX3qdM6XDVZkfa+ntHJEki746UmCD8fa2s6Wg+eXuUgqv+BOPmwJon\nYfVTrguwnykuLqa4uAdvqgYwb9fKttYBuMtXGxqGKRz4At7/gV5P/LZ3wbdn+2wu31dCbJg/MWH+\nLgrQXKwWxbCEEHYerTx7uX5nLBa47gXQHLC6tXvl4l+5Nth+YNmyZQDMnTvX2EA6Yd4WuXCOrK9g\n8ff0peFz3ge/ns06OdHQzIbDZYxMltZ4T4xMDqOqrpmvu5qGeCaLVR+7GHsrrPoDrHnaNQEKjyKJ\nXHQsc6k+R3nACLhjiT6Q1kNrD5bRbNdITwl3QYDmNTI5DIuCT3cX9fzFFqs+LXHMd/XZRV9LMjc7\n83atiL7Z9S589CNInKBvOxbQu0S8fN8xgv1spMTItm49EeRnY3BcCOsPlqFpGqqns1AsVvifl/SP\nVz6hLxi69DGZzWJS0iIXZ9v6T/jwhzBwGtzxYa+TeIvdwar9JQxPCsVqkQTSU+kpYRRXNrCroJcb\n/7Yl8/O+B2uf0UvgyjxzUzJvi3zsH42OwDtt+Dt8+b9wzhX6npE9mJ1ypq1HK6ioayY9peddMgLS\nU8L5cFM+S3cVMi6pl11TFivMegH8w2H9C9BQBf/zcrenjgq47LLLjA6hS+ZN5DFTjY7Auzgc8N/H\nYcPf9C3abngDbH59uuRnu4vwtVkYJvPHeyUi2JekqABWHyjhtzNH9v5CSsGVv4eACFjxf9BQra8D\n6MObdH+SnJxsdAhdMm/XSul6/RBda2nUa3Vs+Bucfy/ctKDPSbzF7uDzPUWMTArFz8fqnDj7odEp\n4WSX1JJVcqLvF7vwZzDzL3DwK3hrFtT0YJ56P5aXl0deXp7RYXTKvIl81//qh+hcfQX863rIWAKX\n/x9c/Wf91/E+2pR9nPLaJs5Li3RCkP3XuLQIABZtdVIimXgX3Pw2FO+BNy6DUtl8pSsrVqxgxYoV\nRofRKfMmctG149nwjxmQtwmufwMu+InTZjV8ursQfx8LQ6VbpU9iQv1JiQnkq73FaJrmnIuOvA7m\nfgbNdfDmFZD9tXOuKwwjiby/OrIaXr8EThTDnA9gzE1Ou3Sz3cGyvcWMSg7DV4pk9dn4tEjyj9ez\nLa/SeRdNmgB3r4DgOPjXd2D7v5x3beF28lPW32gabHxF704JjoN7V0HaxU69xcr9JVTUNTN+8Nkb\nMIueGzcoAouChVtynXvhiIFw11f6NNOP5+lVLVuanHsP4RaSyPuT5npYOg+W/RqGzoC7/wuRzq/q\n9t7WPMIDfRiSIJtIOENIgA9DE0JZta+EFruT54EHhMOcJTD1Qdj6JiyYCdW9WE0qDGXeRD7+ef0Q\nutIseP0y2PlvuOhXcMu/e1w3pTtKqhtYdaCU8YMjZRGQE40fHMnxmia+yHRBFT6rTZ+eeNMCOJYB\nr14EOeucfx8vNWPGDGbMmGF0GJ0ybyKPGKcfQl9u/9p0qCnWd2C/9NFubwjRU0t2FGB3aJw/VLpV\nnGn0wHACfK38e9NR191k1HfgnhX6G/xbs2DVH2X7OCAuLo64uDijw+iUeRN58XL96M8aquGjB+DD\neyFhHPxoHQxx3WYbDofG4q15pMUGExMqJWudyddm4fwhUWw9UkF+ZZ3rbjRgBPxwjV5wa81TsOAa\nqHDhm4cXOHLkCEeOHDE6jE6ZN5HvfUI/+qvstfDyNNj1H7jol/C9jyE0waW3XHeojCOltUwdFu3S\n+/RXU4ZFY3dovPlNtmtv5BcC33lZn5J6LBNeuRB2L9YHyvuhr7/+mq+/9uwpmuZN5P1Vcz188TC8\nda1eT+MHX+pV76yur8bwz2+yCQvwYUyqlKx1hQFh/pwTF8zHOwqdP+jZnjE3wY/WQswwWHIPLLwV\nqgtdf1/RY5LIzeTIar0Vvullfan9j9ZC8vluuXV2WS2rDpQyeVhU93e0ET12wYgYymuaWLTNTUvG\nIwfBD5bBlX+AI6vg75Nhx7/7bevcU8lPnBnUlMAH98DbswENvrcUrnkafN1XA/yt9TnYLIqpw2Lc\nds/+KD0lnJhQP95cm+28lZ5dsVhh6jy4bz3EpcPSB+Dt66Bkv3vuL7okidyb2VtgyxvwtwmQ+ZG+\nk/p9GyBtulvDKDnRwMLNuZw3OIKQQCmP6koWi2J6eizZpbX8d/8x9948ajDc+aleeKtoN7wyDZb9\nr14aVxjKvGVsz3/V6AhcR9P0CnZfPQ5lByD1Qv2HK2aoIeG8sTabZruDK8Z49hQts5gwOJIvthfy\n4spDXDnCzX/nFoteeGvk/8DK38HGl2DPe3D5fBj7XacUXPM01157rdEhdMm8LfLQYfphNsV79C6U\n/9wMjha45R248xPDknh5TSP/2nCU89IiiZIph27hY7MwPT2WPXlVrDxQYkwQQVEw669wz0oIT4Gl\n9+vjM/s/N13/eXR0NNHRnj0Ty7yJPP8T/TCLYxn6bvavXADFu2HGU3D/RhhxraH7ML6y5jANzXau\nGCutcXe6YEQMoQE+PLlsv/v6ytuTeB7c9V99Vai9CRbdCm9eaaqVoQcOHODAAc8u92veRL7/Wf3w\ndm0J/OWpcGilPif8wR0w+Udg8zU0tJyyWhasz+H8IVHEhElr3J18bRauHBdHVtEJPtpt8JRAi0Vf\nFfrAJr2VXpWn12z5x9VwcLnXt9A3bNjAhg0bjA6jU+ZN5N5M0/SphO/crCfww6v0+ig/2a3PCQ+I\nMDpCAP70xT6sFsU157l2oZFo36Sh0USH+PHUF/tpbLYbHY6+bmH8XL2hMeNJqDwK79yg127Zu0SW\n+7uQJHJP0tKoz9F95QK9H7xwO0x/BB7apddHCfSc3XbWHyrjy4xjXDo6VmaqGMRqUXxnchLFlQ08\ns9yDfvX3CYDJ98GDO2H23/UNLN7/PrwwDtb+BWrLjI7QdMw7a8WblB6A7W/DrkVQVwYDRuo/AOk3\ngo/ndVnUNbXw8JI9DAj146JRsUaH06+NSApjzMBwFqzL4dYJKaTFBBsd0rdsvnDuHBh7Kxz4HDa/\npm/+vPpJSL8eJt4NieMNHeMxC0nkRmms0ed+b39b32rNYoNhV8OEH0DaJR79n/vpLw+Qe7yOH18z\nVHYA8gCzz0/iQEE18xbt4JMHLvC88sEWK4yYpR8l+2HL63qjZddCiB6mT1sccwuEJRodqddS7hjx\nnjBhgrZ161aX3+c0ta1LmIOS3XvfzjTVwaH/6v2FWV9CSz1ED4Vz79D/MwcPMDrCLn1zqIw5b27i\nguExfGeyB/3d9nObD5azaN1R7r9kML+6arjR4XStoRoyPtSTee4GQOk7VY2+WW/QeFA3YlWVvuAp\nLCzM7fdWSm3TNG1Cl+eZNpF7ioZqvUZF5sdw4AtoroWgGBg5G0bfBMmTPLr1faqiqnpmvrAOf18r\nP545FD8f8y3+8FaapvH26mz2HK3knXsmMSXNs+c9n+b4Eb1m/q6F+gCpxaYvchsxC4ZfCyH9t/tO\nEvnRd/XHgbe4974AZYfg4JeQtQyObgBHMwRE6ruXj7pe3yPRDdUInamh2c6tr29kX9EJfjZrGNEy\n3dDj1DW28PwnB2hqsfPZjy8kOTLQ6JB6RtP0Af59n+gNn+OHAQVJE+Gcy2Hwpfq8dTevHt27dy8A\n6enpbr0vSCKH5dP1x8tXu/5eVQX6AoictfpRkaM/HzMChl6p74+ZdL7XJe82LXYH972zneWZx5h7\n6SBGD/SM6Y/ibMcqG3jhswPEhPjxyQPTCAs0dq1Br2kalO7Xk/qBz6FwJ6CBf7jeBTP4Mki9QN9z\n1sW/0S5YsACAuXPnuvQ+7eluIvfOzGIkhwPKD0LBNr1vL2ed/qsh6P/JUi+AKfNgyJX6LuVezu7Q\neHjJHv6beYwbJidJEvdwseH+fG/6IN5YfpibXtvI+z+aQqi/F04PVUrfrWjACLj4V1BbrndRHl4F\nh1dA5lL9vKABkDIZUqbAwCkQO9prG0x90f++457QNKgu0FsDBdv0o3AHNFbrX/cLg9Rp+jSq1Ash\nNt1le2EaobHFzs/e3cVne4qYMS6OaSM8fzBWwLDEUO68ZBBvrcrmxlc28M5dk4gJ8TM6rL4JioLR\nN+pHW2v96HrI3agf+z7Wz/MJgrjR+taG8eMgfqw+ocDkyd3c311PNFRDyT4oydCXxR/L1D9uK9Fp\nsUHsKH2AMnG8fkQPMWW1N4CS6gbmLdzB5uzjzJ6YyMXp/XfAyRulp4Rz5yWD+NfqbGa+uJa3vj+R\nEfHun3XhEqe21ifepT9Xla8n9LzNULRTn9bb/Ir+NVuAXkc9Nl3f7ShmGMQMh5B4r5lo0JX+lcib\n6qAiG8oP6wMp5YfheLb+8Ymib8/zC9UX5aTfoD/GjYH4MfqKtX5gTVYpP1+8k+qGFu64OJVz0zxn\nKpjovvSUcOZdM5Q3lx9h9t/X8/DVw/n+1FSUSZLXacKSvm2xAzjsUHYQinbpib1olz7dsaHy29f4\nheqt9Zjheq31iNRvj4AIr0ry5hns1DR96W9Vnt4dcvyA/lhToQ9GVuXDiTOKCwXFQORg/R8xajAM\nGAWxIyEs2av+EZ3lWHUDT36xnw93FBAX7s+dlwwiNrx/vHmZWXVdM++uO8q+gmrOTQnn97PTSU80\nSeu8JzQNakv1bpnSA61H68e1Z5QD9gvVx7giUqkLGggh8QRGJ0JIAoTGQ3CcW4rWmWPWiqbp76A1\npfpfdE2J/g9RU9L6eevztaVw4hjYG09/vc0fQhP1FWOhSfr+g5FpetKOTAP/fvifuR1FVfW8sTab\nf288it2hcenoWC4bE4ePrNo0DU3T2JhVzufbC6lraOGaMfHcP30woxLkZwDQV1pXHtVnnJ11HD07\nt6AgKFrvnglNgJA4vWEYFKM/f/LjGL1138suWHMk8t3vwZK7z35eWU//ywoeAMGx+q9XoYn6Y9Va\n8A2Gwd/v+zdgQrWNLXydVcoH2wtY2bpl2IRzorhqXBwRwV4+MCY6VN/YwvLdxXyzv4ymFgfnpYRz\nw/gkrk6PJzLIS6cqutjOHTugqZZxA8Ogukj/zf60xyI4UQz1x0FznH2BKfPgqj/06t7mSOTlh/VF\nNUEDIDim9XGAvrimq9kh7pxH7gWq6pvZnV/JjtxKth6tYOORcppaHIQG2Jg4JIopQ6OJ9PaZDaLb\n6htb2JBVxpaD5RyrasSiYFRCGFPPieL81EhGxIcSH+Zvzv70Hur2PHKHHeor9B6Ck0e5PtA6cGqv\n7m2OeeRRg2HKA0ZH4XFa7A4aWhw0NNtbDwf1TXYq6po4XvvtUVhVz9HyOrLLajle2wSAAuIi/Jky\nLJoxKWGkDAj2vCJLwuUC/GxcOjqOS9JjKaqoZ8/RSg4WneCNtdm8ukZfFxHkZ+WcAcEkhQeSEO5P\nXFgA0cG+hPr7EBpgI8Tfh1B/HwJ8rPjYFL5WC1aL6r/J39LWUxANjHDrrT26Rb4mq5Tff5oJcHI7\nq5PRnhJ224enfi/PRf8UgB+X/KXT805+TTvrsj2+57fPnXre6Rc+/frdvUbrlTRobHFgd3T9b6YU\nhAX4EBXqR3SoHwNC/UiKCiQpOogAX3NOmRR919TiIL+sjuLKeooq6jlW2UBVXTOVtU0027vx/w7w\nsVrwsSl8rBasSoH+B4tS6J+2Pp75sVIovn3eHbrznjO2bicAe4LO7dU9vjsxmbsvTOvVa03RIg/2\nszEsNuTbJ9RpD6e983/7XOtrHTbsDo0BEf7tnscZ55/+3NlPqrM+OPXDUz5SZ512Vtztnd/1vfT/\n8Dar/gOi/7BY8LEqbFYLvjYLQf42gvxsBPvbCPCzYmm9uLW/tpBEjwX4WBkSH8KQ+G9/7uyahqZp\n1DXaqW1sob7JTkOTnfrWo7nFQYtDw25vfXRotNg1WhwOHA442ZzR9EbJyQbSmZ+3PueureG6exdr\no8Ki1Om5qAei3TDm5NGJfPzACMb3dkn4cr3A/pLLe9c3JYQQAAsWHARg7u3nGRxJxzw6kffJ9M+N\njkAIYQK333670SF0ybyJ3OZlJTyFEB7Jx8fzi46Zd8VH1kv6IYQQfbBlyxa2bNlidBidMm8iz12s\nH0II0QcZGRlkZGQYHUanzJvIhRCin5BELoQQXk4SuRBCeDlJ5EII4eXcskRfKVUKHHX5jc4WDZQZ\ncF8hhLkYlUsGapoW09VJbknkRlFKbe1OnQIhhOiMp+cS6VoRQggvJ4lcCCG8nNkT+WtGByCEMAWP\nziWm7iMXQoj+wOwtciGEMD2PTuRKqX8opUqUUntPeS5SKfVfpdTB1seIM14zUSllV0rdeMpzy5RS\nlUqpTzu518+UUplKqd1KqRVKqYGu+a6EEO7Wk1yilPqlUmpn67G3NZ9Etn7tp0qpjNbnFyql/Nu5\nV4pSapVSakdrPrnG1d+fRydyYAEw44znHgZWaJo2BFjR+jkASikr8BTw5RmveRq4o4t77QAmaJo2\nBngf+HPvwxZCeJgFdDOXaJr2tKZp4zRNGwc8AqzRNO24UioReBA9T6QDVuC77dzrMWCxpmnntn7d\n5WVYPTqRa5r2NXD8jKdnA2+1fvwW8D+nfO3HwAdAyRnXWQGc6OJeqzRNq2v9dCOQ1MuwhRAephe5\npM2twMJTPrcBAUopGxAIFLZ3OyC09eOwDs5xKo9O5B2I1TStCKD1cQBA67vld4BXnHCPu4AvnHAd\nIYTnajeXtFFKBaK34j9oPacAeAbIBYqAKk3TvmrnuvOBOUqpfOBz9AamS3ljIu/I88CvNU2z9+Ui\nSqk5wAT07hghRP81C/hG07TjAK196LOBQUACENSaL850K7BA07Qk4BrgX0opl+Zab0zkx5RS8QCt\nj23dKBOARUqpHOBG4CWlVHu/KnVIKXU58ChwnaZpjc4LWQjhgTrKJW2+y+ndKpcD2ZqmlWqa1gws\nAdrb3f0uYDGApmkbAH/0Wi0u442J/GPgztaP7wSWAmiaNkjTtFRN01LRByvv1zTto+5eVCl1LvAq\nehI/8x9UCGE+7eYSAKVUGHDxqc+hd6lMVkoFKqUUcBmwr53r5rZ+DaXUCPREXur06E/h0YlcKbUQ\n2AAMU0rlK6XuAp4ErlBKHQSuaP28q+usBd4DLmu9zlWtz/9OKXVd62lPA8HAe63Tjj52wbckhDBA\nL3LJd4CvNE2rbXtC07RN6I3E7cAe9Pz5Wuv1T80lPwfuUUrtQm/Rz9VcvPJSVnYKIYSX8+gWuRBC\niK5JIhdCCC8niVwIIbycJHIhhPByksiFEMLLSSIXQggvJ4lcCCG8nCRy0S+11ou+ovXjJ5RSLxgd\nkxC9ZTM6ACEM8lvgd0qpAcC5wHVdnC+Ex5KVnaLfUkqtQS/LMF3TtE7r1QvhyaRrRfRLSqnRQDzQ\nKElceDtJ5KLfaS1Z+g56benatiJqQngrSeSiX2nd9WUJ8HNN0/YBv0ff0UUIryV95EII4eWkRS6E\nEF5OErkQQng5SeRCCOHlJJELIYSXk0QuhBBeThK5EEJ4OUnkQgjh5SSRCyGEl/t/gUO307TcOrEA\nAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "samples = np.array([1035., 1050., 1020., 1055., 1046.])\n", "sigma = np.sqrt(np.var(samples, ddof=1))\n", "sample_mean = np.mean(samples)\n", "mu = 1060\n", "w = mu - sample_mean\n", "\n", "x = np.linspace(mu - w - 2, w + mu + 2, 1000)\n", "y1 = ss.norm.pdf(x, loc=mu, scale=3)\n", "y2 = ss.t.pdf(x, loc=mu, scale=sigma / np.sqrt(len(samples)), df=len(samples) - 1)\n", "plt.plot(x,y1, label='single sample')\n", "plt.plot(x,y2, label='multiple samples')\n", "plt.fill_between(x, y1, where= np.abs(x - mu) < w, color='lightblue')\n", "\n", "plt.axvline(mu - w, linestyle='--', color='orange')\n", "plt.axvline(mu + w, linestyle='--', color='gray')\n", "plt.xticks([mu - w, mu + w], [sample_mean, mu + w])\n", "plt.yticks([])\n", "plt.ylabel(r'$p(x)$')\n", "plt.xlabel(r'$x$')\n", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "13.953494186 1041.2 1060.0\n", "-3.01272766638\n" ] } ], "source": [ "samples = np.array([1035., 1050., 1020., 1055., 1046.])\n", "sigma = np.sqrt(np.var(samples, ddof=1))\n", "sample_mean = np.mean(samples)\n", "true_mean = 1060.\n", "\n", "print(sigma, sample_mean, true_mean)\n", "\n", "T = (sample_mean - true_mean) / (sigma / np.sqrt(len(samples)))\n", "print(T)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now we have a $T$ and we know $P(T)$. However, just like the $zM$ test, we can't compute $P(T)$ since that's a single point and we're using a continuous distribution. So instead, we build an interval and see how big it must be to catch that $T$. \n", "\n", "Specifically, we want to find $\\int_{-T}^T p(t)\\,dt$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "T = -3.01272766638\n", "0.0197221623052 Is the single sided p-value\n", "0.0394443246103 Is the double sided p-value\n", "notice, just using 2 * the single-sided value gives the same answer\n" ] } ], "source": [ "print('T = ', T)\n", "\n", "p = ss.t.cdf(T, len(samples) - 1) # The integral from -infinity to T\n", "print(p, 'Is the single sided p-value')\n", "\n", "p = 1 - (ss.t.cdf(-T, len(samples) - 1) - ss.t.cdf(T, len(samples) - 1)) # 1- The integral from -T to T\n", "print(p, 'Is the double sided p-value')\n", "\n", "print('notice, just using 2 * the single-sided value gives the same answer')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "What if accidentally reverse our T-value?\n", "---" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.01272766638\n", "CDF gives: 0.980277837695\n", "Recognize that includes from -infty up to a positive T, so we need to find the complementary area\n", "0.0197221623052 Is the single sided p-value\n", "0.0394443246103 Is the double sided p-value\n" ] } ], "source": [ "T = (true_mean - sample_mean) / (sigma / sqrt(len(samples)))\n", "print (T)\n", "\n", "p = (scipy.stats.t.cdf(T, len(samples) - 1))\n", "print ('CDF gives: ', p)\n", "print ('Recognize that includes from -infty up to a positive T, so we need to find the complementary area')\n", "\n", "print (1 - p, 'Is the single sided p-value')\n", "print ((1 - p) * 2, 'Is the double sided p-value')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Summary of Methods for Comparing Single Measurement with Normal Population\n", "====" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "$zM$ Test\n", "====\n", "\n", "**Data Type:** Measurements and Ranks\n", "\n", "**Compares:** A single sample vs a normally distributed population\n", "\n", "**Null Hypothesis:** The sample came from the population\n", "\n", "**Conditions:** Standard deviation of population is known\n", "\n", "**Related Test 1:** Student's $t$-test, which is used when the standard deviation is not known\n", "\n", "**Python:** Integrate an interval with a Z-statistic and `erf` or `scipt.stats.norm.cdf`" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Student's $t$-test\n", "====\n", "\n", "**Data Type:** Measurements and Ranks\n", "\n", "**Compares:** A set of samples vs a normally distributed population\n", "\n", "**Null Hypothesis:** The sample came from the population\n", "\n", "**Conditions:** Standard deviation of population is not known\n", "\n", "**Related Test 1:** $zM$ test, which is used when standard deviation is known\n", "\n", "**Python:** Integrate an interval with a T-stastic and `scipt.stats.t.cdf`" ] } ], "metadata": { "celltoolbar": "Slideshow", "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.1" } }, "nbformat": 4, "nbformat_minor": 1 }