{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hypothesis Testing\n", "\n", "In this section we will talk about the process of using data to test a hypothesis. \n", "\n", "We can roughly break down the process of hypothesis testing into the following steps:\n", "\n", "1. State relevant null and alternative hypotheses\n", "2. State the test statistic\n", "3. Derive the distribution of the test statistic under the null hypothesis\n", "4. Select significance level ($\\alpha$) at which the null hypothesis will be rejected\n", "5. Compute an observed test statistic from the data\n", "6. Calculate a p-value and decide whether to reject the appropriate hypothesis\n", "\n", "Let's go through these steps one by one, and look at what each of these steps looks like for a simple coin flip example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. State relevant null and alternative hypotheses\n", "When constructing an analysis, we should be careful to state the _null_ and _alternative_ hypotheses. Typically, the _null hypothesis_ ($H_o$) refers to a statement or default position that there is no relationship between two phenomena, that there is no signal, or that our current models describe the data well. The _alternative_ ($H_a$) hypothesis usually describes the opposite: that there is a relationship between two phenomena, that there is a signal, or that our models to NOT describe the data well. Here are some example pairs of hypotheses you might encounter in IceCube:\n", "\n", "- $H_o$: Neutrino event arrival directions are correlated with the positions of blazars\n", "\n", "- $H_a$: Neutrino event arrival directions are not correlated with the positions of blazars\n", "\n", "\n", "\n", "- $H_o$: A 3-flavor oscillation model describes our data well\n", "\n", "- $H_a$: The data is inconsistent with the 3-flavor model\n", "\n", "\n", "\n", "- $H_o$: The prompt component of the atmospheric neutrino spectrum is zero\n", "\n", "- $H_a$: The prompt component of the atmospheric neutrino spectrum is nonzero\n", "\n", "\n", "For this talk, let's consider a simple example: We want to know if a coin flip is \"fair\" in the sense that the $P(heads)=P(tails)=1/2$. In this case, we can set up our hypotheses like so:\n", "\n", "- $H_o$: The coin is fair, $P(heads)=1/2$\n", "\n", "- $H_a$: The coin is not fair, $P(heads) < 1/2$\n", "\n", "This is a _one tailed_ hypothesis tests, where we are only interested in the case that $P(heads) < 1/2$. If $P(heads) > 1/2$, then we won't be able to reject the null hypothesis. You could also set up your hypothesis to be _two tailed_:\n", "\n", "- $H_o$: The coin is fair, $P(heads)=1/2$\n", "\n", "- $H_a$: The coin is not fair, $P(heads) \\neq 1/2$\n", "\n", "...where this would allow you to reject the null hypothesis both if $P(heads) < 1/2$ or if $P(heads) > 1/2$. For now let's use the one tailed test, and we'll come back and look at this two tailed test later. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. State the test statistic\n", "\n", "A _test statistic_ is a quantity derived from the sample, intended to be a numerical summary of a dataset that reduces the data to one value that can be used for hypothesis testing. In general, the test statistic is selected in such a way as to quantify, within observed data, behaviors that would distinguish the null and alternative hypotheses. These are very similar to estimators, and we can oftentimes use estimators as test statistics. \n", "\n", "An important property of the test statistic is that its sampling distribution under the null hypothesis must be calculable either exactly or approximately. There are many types of test statistics that you can use, and what the exact form of the test statistic is will usually be determined by your analysis. Here are some examples of test statistics used in IceCube:\n", "\n", "- $\\chi^2 = \\sum_{i \\in {bins}} \\frac{(n_i^{\\nu+\\mu_{atm}} - n_i^{data})^2} {(\\sigma_i^{data})^2+(\\sigma_{\\nu+\\mu_{atm,i}})^2}$ (Chi squared TS used for oscillation analyses, usually with an additional systematic term)\n", "\n", "- $TS = -2 \\log(\\frac{\\mathcal{L}(n_s=0)}{\\mathcal{L}(\\hat{n}_s)})$ (Likelihood ratio TS used in point source analyses)\n", "\n", "The *likelihood ratio* is a very powerful test statistic. Given two hypotheses $H_o$ and $H_a(a)$ where $H_o = H_a(\\tilde{a})$ (they are _nested_ models), then the Neyman-Pearson lemma states that the likelihood ratio test has the strongest statistical power (i.e. correctly reject $H_o$ when $H_1$ is true). This is the case for many point source searches, where the null hypothesis (\"There are zero signal events, $n_s$ = 0) is nested in the alternative hypothesis (\"There are a nonzero number of signal events, $n_s$ > 0). \n", "\n", "You'll hear more about the likelihood ratio test later this week in the point source example session. For now, let's go back to our simple example of a coin flip. For our coin flip example, let's use just use the fraction of heads obtained as our test statistic:\n", "\n", "- $TS = (N_{heads})/(N_{flips})$\n", "\n", "This may not be the _best_ test statistic that we could possibly use, but it does satisfy the property of the distribution being calculable under the null hypothesis. This does have an analytic result (the binomial distribution), but let's practice generating a distribution from monte carlo." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "%matplotlib inline\n", "\n", "\n", "def flip_coin(pheads=0.5):\n", " #This function flips a single coin, returning heads with probability pheads, and tails with probability 1-pheads\n", " flipfloat = np.random.uniform(0,1)\n", " if flipfloat <= pheads:\n", " flipresult = 'H'\n", " else:\n", " flipresult = 'T'\n", " return flipresult\n", "\n", "def do_one_trial(pheads=0.5, Ncoins=20):\n", " #This function performs 1 trial of flipping many coins (set by Ncoins). \n", " #You can adjust the value of pheads for these trials with the pheads arg\n", " trialdata = []\n", " for i in range(0,Ncoins):\n", " trialdata.append(flip_coin(pheads=pheads))\n", " heads = [trial for trial in trialdata if trial=='H']\n", " fracheads = float(len(heads))/Ncoins\n", " return fracheads\n", "\n", "def do_many_trials(pheads=0.5, Ntrials=1000):\n", " #This does many trials, with each trial flipping 1000 coins.\n", " #we use this (with pheads=0.5) to calculate our null hypothesis test statistic distribution\n", " resultarr = []\n", " for i in range(0, Ntrials):\n", " resultarr.append(do_one_trial(pheads=0.5, Ncoins=20))\n", " return resultarr\n", "\n", "#fair_results is an array of the results of 1000 trials,\n", "#with each element identifying the fraction of heads obtained in that trial\n", "fair_results = do_many_trials(pheads=0.5, Ntrials=1000)\n" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFr9JREFUeJzt3Xm0ZWV55/Hvj8FgC+VEaSpgWU6gFYcSrxKHqGBinAJqO6FoNK6UremoaExLDxHtlU67bEy6O9pCqw1qMDg2ONBBGURZgBZSMpRibEQaUt0UUQEnJp/+Y+9LHcpb55576+59bt39/ax11tnT2fu577r3PPd9373fN1WFJGl49ph2AJKk6TABSNJAmQAkaaBMAJI0UCYASRooE4AkDZQJQJIGygQgSQNlApCkgdqry5MnuRq4GbgDuL2qZsYdv//++9e6deu6DEmSVpSLL774hqpavZjPdpoAWodV1Q2THLhu3To2bdrUdTyStGIk+cFiP2sTkCQNVNcJoIAzk1ycZONcByTZmGRTkk3btm3rOBxJ0qyuE8BTquoQ4NnAHyd56o4HVNWJVTVTVTOrVy+qGUuStAidJoCquq59vx74LPCELq8nSZpcZwkgyT2S7De7DDwTuLyr60mSFqbLu4DuD3w2yex1Tqmq/9Xh9SRJC9BZAqiqq4DHdHV+SdKu6eM5AGm3dspF13Da5uumHQZHbjiAlx+6dtphaAXxOQBpHqdtvo4tW2+aagxbtt60LJKQVhZrANIE1q9Zxamve+LUrv/SEy6Y2rW1clkDkKSBMgFI0kCZACRpoEwAkjRQJgBJGigTgCQNlAlAkgbKBCBJA2UCkKSBMgFI0kCZACRpoEwAkjRQJgBJGigTgCQNlAlAkgbKBCBJA2UCkKSBMgFI0kCZACRpoEwAkjRQJgBJGigTgCQNlAlAkgbKBCBJA2UCkKSBMgFI0kCZACRpoEwAkjRQJgBJGqjOE0CSPZNckuTzXV9LkjS5PmoAbwK+3cN1JEkL0GkCSHIg8Fzgg11eR5K0cF3XAP4a+DPglzs7IMnGJJuSbNq2bVvH4UiSZnWWAJI8D7i+qi4ed1xVnVhVM1U1s3r16q7CkSTtoMsawJOBI5JcDfwdcHiSj3V4PUnSAnSWAKrq2Ko6sKrWAS8Dzq6qo7u6niRpYXwOQJIGaq8+LlJV5wLn9nEtSdJkrAFI0kCZACRpoEwAkjRQE/UBJDkAeODo8VV1XldBSZK6N28CSPJu4KXAFuCOdnMBJgBJ2o1NUgN4PnBwVd3SdTCSpP5M0gdwFbB314FIkvo1SQ3gZ8DmJGcBd9YCquqNnUUlSercJAng9PYlSVpB5k0AVXVykrsBB7Wbrqyq27oNS5LUtUnuAno6cDJwNRDgAUn+wNtAJWn3NkkT0PHAM6vqSoAkBwEfBx7XZWCSpG5NchfQ3rNf/gBV9V28K0iSdnuT1AA2JfkgMDuZyyuATd2FJEnqwyQJ4PXAHwOzt31+FXh/ZxFJknoxyV1AtwDvbV+SpBVipwkgySeq6iVJLqMZ++cuqurRnUYmSerUuBrAm9r35/URiCSpXzu9C6iqtraLb6iqH4y+gDf0E54kqSuT3Ab6u3Nse/ZSByJJ6te4PoDX0/yn/5Akl47s2g84v+vAJEndGtcHcApwBvCXwNtHtt9cVT/sNCpJUud2mgCq6kbgxiT/aodd+ybZt6qu6TY0SVKXJnkQ7As0t4EG2Ad4EHAl8JsdxiVJ6tgkD4I9anQ9ySF4F5Ak7fYmuQvoLqrqm8ChHcQiSerRJPMBvGVkdQ/gEOAfO4tIktSLSfoA9htZvp2mT+DT3YQjSerLJH0A7wRIsqpZrZs7j0qS1Ll5+wCSzLQDwl0KXJbkW0lmug9NktSlSZqAPkwzHtBXAZI8pd3maKCStBub5C6gO2a//AGq6ms0fQFjJdknydfbGsMVSd65K4FKkpbWuLGADmkXv5LkBJqJ4At4KXDuBOe+BTi8qn6SZG/ga0nOqKoLdzFmSdISGNcEdPwO6+8YWf6VCWJ2VFUF/KRd3bt9zfs5SVI/xo0FdNiunjzJnsDFwEOB91XVRbt6TknS0hjXBHR0VX1shwfB7lRV884RXFV3ABuS3Av4bJJHVtXlO1xnI7ARYO3atQsKXpK0eOM6ge/Rvu+3k9fEqurHwDnAs+bYd2JVzVTVzOrVqxdyWknSLhjXBHRC24RzU1X91UJPnGQ1cFtV/TjJ3WlmFnv34kOVJC2lsbeBtk04Ry3y3GuAc9rZxL4BfKmqPr/Ic0mSltgkD4Kdn+RvgFOBn85ubEcF3amquhR47K6FJ0nqyiQJYEP7/q6RbQUcvvThSJL6MkkCeG1VXTW6IcmDO4pHktSTSYaC+NQc2z651IFIkvo17jmAh9PM+3vPJC8c2bWKZm5gSdJubFwT0MHA84B7Ab8/sv1m4I+6DEqS1L1xzwGcBpyW5IlVdUGPMUkAnHLRNZy2+bpph8GWrTexfs2qaYchLblJ+gBekGRVkr2TnJVkW5KjO49Mg3fa5uvYsvWmaYfB+jWrOHLDAdMOQ1pyk9wF9Myq+rMkLwCuBl4InAd8rMvAJGi+fE993ROnHYa0Ik1SA9i7fX8u8MmqurHDeCRJPZmkBvC5JN8Bfg68vh3j5xfdhiVJ6tq8NYCqejvwJGCmqm6jGQ7iyK4DkyR1a9xzAIdX1dmjzwAkGT3kM10GJknq1rgmoKcBZ3PXZwBmFSYASdqtjXsO4B3t+2v6C0eS1JdxTUBzTgU5a5IpISVJy9e4JqDZaR8PBh4PnN6u/z7w9S6DkiR1b1wT0DsBkpwHHFJVN7frxwFf6CU6SVJnJnkQ7P7ArSPrt7bbJEm7sUkeBPsI8PUkn23Xnw+c1FlEkqRezJsAquovkpwB/Ha76TVVdUm3YUmSujZJDWB2Avixk8BLknYvk/QBSJJWIBOAJA2UCUCSBmrck8A304z58yu7gKoq58iTpN3YuAfB9tvZPknS7m9cDeA+4z5YVT9c+nAkSX0ZdxvoxTRNQJljXwEP7iQiSVIvxjUBPajPQCRJ/Zr3QbAkT51re1Wdt/ThSJL6MsmTwG8bWd4HeAJN89DhnUQkSerFJGMB3WVKyCQPAP66s4gkSb1YzINg1wKPmO+gJA9Ick6SLUmuSPKmRVxLktSRSfoA/ivbHwjbA9jAZAPD3Q68taq+mWQ/4OIkX6qqLYuOVpK0ZCbpA9g0snw78PGqOn++D1XVVmBru3xzkm8DBwAmAElaBibpAzh5Vy+SZB3wWOCiXT2XJGlpTNIE9GTgOOCB7fGzYwFN9CBYkn2BTwNvrqqb5ti/EdgIsHbt2okDlyTtmkmagD4EHENz6+cdCzl5kr1pvvz/tqo+M9cxVXUicCLAzMzMXIPPSZI6MEkCuLGqzljoiZOEJnl8u6reu+DIJEmdmiQBnJPkPcBngFtmN7bTRI7zZOCVwGVJNrfb/nVVfXFRkUqSltQkCeDQ9n1mZFsxz5PAVfU15h5ITpK0DExyF9BhfQQiSeqXU0JK0kCZACRpoHaaAJK8uH13XgBJWoHG1QCObd8/3UcgkqR+jesE/qckZwIPSnL6jjur6ojuwpIkdW1cAngucAjwUeD4fsKRJPVl3JzAtwIXJnlSVW1rx/Shqn7SW3SSpM5MchfQ/ZNcAlwBbElycZJHdhyXJKljkySAE4G3VNUDq2ot8NZ2myRpNzZJArhHVZ0zu1JV5wL36CwiSVIvJhkL6Kok/46mMxjgaOCq7kKSJPVhkhrAHwKraUYD/TSwf7tNkrQbm2QwuB8Bb+whFklSjxwLSJIGygQgSQNlApCkgRo3GuiZI8vH7uw4SdLuaVwNYPXI8ou7DkSS1K9xCaB6i0KS1Ltxt4E+uB0GOiPLd3I4aEnavY1LAEeOLP+nrgORJPVr3HDQX5ldTrK63batj6AkSd0bdxdQkrwjyQ3AlcB3k2xL8uf9hSdJ6sq4TuBjgKcAj6+q+1TVvYFDgScnOaaX6CRJnRmXAF4JHFVV35/dUFVX0YwG+qquA5MkdWtcAti7qm7YcWPbD7B3dyFJkvowLgHcush9kqTdwLjbQB+T5Caa5wBg+4NhAfbpNCpJUufG3Qa6Z5+BSJL6tdMEkGQf4F8ADwUuBT5cVbf3FZgkqVvj+gBOBmaAy4DnAMf3EpEkqRfj+gDWV9WjAJJ8CPj6Qk6c5MPA84Drq+qRiw9RktSFcTWA22YXFtn0cxLwrEV8TpLUg0nuAoLmzp+7j9wVVFW1atyJq+q8JOuWJEpJ0pLzLiBpN7Fl60289IQLph0GR244gJcfunbaYWgJjKsB9CLJRmAjwNq1/lJJczlywwHTDgFokhBgAlghpp4AqupE4ESAmZkZZyGT5vDyQ9cuiy/d5VAD0dIZ1wksSVrBOksAST4OXAAcnOTaJK/t6lqSpIXrrAmoqo7q6tySpF1nE5AkDZQJQJIGygQgSQNlApCkgTIBSNJAmQAkaaCm/iSwlqdTLrqG0zZfN9UYtmy9ifVrxo45KGkXWAPQnE7bfN2d475My/o1q5bNGDjSSmQNQDu1fs0qTn3dE6cdhqSOWAOQpIEyAUjSQNkEtAzZASupD9YAliE7YCX1wRrAMmUHrKSuWQOQpIEyAUjSQJkAJGmgTACSNFAmAEkaKBOAJA2UCUCSBsoEIEkDZQKQpIEyAUjSQJkAJGmgTACSNFAmAEkaKBOAJA2UCUCSBsr5ACQtyJatN/HSEy6YagxHbjiAlx+6dqoxrAQmAEkTWw6zxM3OlmcC2HUmAEkTe/mha6f+xTvt2sdK0mkfQJJnJbkyyfeSvL3La0mSFqazBJBkT+B9wLOB9cBRSdZ3dT1J0sJ0WQN4AvC9qrqqqm4F/g44ssPrSZIWoMs+gAOA/zOyfi1w6LgPXLXtp7bv0XRyrV+zatphSMvWcrgTaSWYeidwko3ARoB91zxkytEsD+vXrFoWd1tIy5F/G0snVdXNiZMnAsdV1e+168cCVNVf7uwzMzMztWnTpk7ikaSVKMnFVTWzmM922QfwDeBhSR6U5G7Ay4DTO7yeJGkBOmsCqqrbk/xL4O+BPYEPV9UVXV1PkrQwnfYBVNUXgS92eQ1J0uI4GJwkDZQJQJIGygQgSQNlApCkgTIBSNJAdfYg2GIkuRm4ctpxLBP7AzdMO4hlwHLYzrLYzrLY7uCq2m8xH5z6UBA7uHKxT7StNEk2WRaWwyjLYjvLYrskix4+wSYgSRooE4AkDdRySwAnTjuAZcSyaFgO21kW21kW2y26LJZVJ7AkqT/LrQYgSepJ7wlgvonik/xaklPb/RclWdd3jH2ZoCzekmRLkkuTnJXkgdOIsw/zlcXIcf88SSVZsXeATFIWSV7S/m5ckeSUvmPsywR/I2uTnJPkkvbv5DnTiLMPST6c5Pokl+9kf5L8l7asLk1yyLwnrareXjTDQv9v4MHA3YBvAet3OOYNwAfa5ZcBp/YZ4zIri8OAf9Yuv37IZdEetx9wHnAhMDPtuKf4e/Ew4BLg3u36/aYd9xTL4kTg9e3yeuDqacfdYXk8FTgEuHwn+58DnAEE+C3govnO2XcNYJKJ4o8ETm6XPwU8I0l6jLEv85ZFVZ1TVT9rVy8EDuw5xr5M8nsB8O+BdwO/6DO4nk1SFn8EvK+qfgRQVdf3HGNfJimLAmYn0L4n8I89xterqjoP+OGYQ44EPlKNC4F7JVkz7px9J4C5JorfcYLPO4+pqtuBG4H79hJdvyYpi1GvpcnuK9G8ZdFWZx9QVV/oM7ApmOT34iDgoCTnJ7kwybN6i65fk5TFccDRSa6lmXvkT/oJbVla6HfKsnsSWHNIcjQwAzxt2rFMQ5I9gPcCr55yKMvFXjTNQE+nqRWel+RRVfXjqUY1HUcBJ1XV8e085B9N8siq+uW0A9sd9F0DuA54wMj6ge22OY9JshdNte6feomuX5OUBUl+B/g3wBFVdUtPsfVtvrLYD3gkcG6Sq2naN09foR3Bk/xeXAucXlW3VdX3ge/SJISVZpKyeC3wCYCqugDYh2acoCGa6DtlVN8JYJKJ4k8H/qBdfhFwdrU9HCvMvGWR5LHACTRf/iu1nRfmKYuqurGq9q+qdVW1jqY/5IiqWvQYKMvYJH8j/5Pmv3+S7E/TJHRVn0H2ZJKyuAZ4BkCSR9AkgG29Rrl8nA68qr0b6LeAG6tq67gP9NoEVDuZKD7Ju4BNVXU68CGaatz3aDo8XtZnjH2ZsCzeA+wLfLLtB7+mqo6YWtAdmbAsBmHCsvh74JlJtgB3AG+rqhVXS56wLN4K/Pckx9B0CL96hf7DSJKP0yT+/ds+j3cAewNU1Qdo+kCeA3wP+BnwmnnPuULLSpI0D58ElqSBMgFI0kCZACRpoEwAkjRQJgBJGigTgCQNlAlASybJHUk2j7zWLcE5n59k/cj6u9qnozuT5OPtcLrH7LD9uCR/usO2q9uHsXb1mk9P8vkFHH/3JF9J8piR8v5hku+3y19Oskc7PPDlSS5L8o0kD2o//+Uk997VuLV7cywgLaWfV9WGne1Mslc7wN9CPB/4PLAFoKr+fBfim1eSXwceX1UP7fI6S+APgc9U1beADQBJTgI+X1WfatePAn4DeHRV/TLJgcBP289/lGbo9b/oO3AtH9YA1Kkkr05yepKzgbOS7Jtmcptvtv+VHjly7Kva/7y/leSjSZ4EHAG8p/2v9iFJTkryovb4Z6SZCOSydrKMX2u3X53knSPXePgcce2T5H+0+y9Jcli760zggPZ6v73An/XoJF9vP3tCkj3b7f8tyaY0k7e8c+T4ZyX5TpJvAi8c2f60kf/qL0my3xyXewVw2jwhrQG2zg6MVlXXzg4hTTNswFEL+fm0Ak17kgNfK+dFMyzB5vb12Xbbq2kGL7tPu74XsKpd3p/msfUAv0kzqNn+7b7Z408CXjRyjZNoxojah2bo24Pa7R8B3twuXw38Sbv8BuCDc8T6VpqhBQAeTjOmzD7AOnY+4cZxNINrbR553dr+HI8APgfs3R77fuBVO/wsewLnAo8eif9h7c//CZr/3mnP8+R2eV9grx3iuBvwf+eIb8eyOrAti83A8cBjdzj+H4D7Tvv3xtf0XtYAtJR+XlUb2tcLRrZ/qapmJ7II8B+SXAp8mWa88vsDhwOfrKobAEaO35mDge9X1Xfb9ZNpZkya9Zn2/WKaL/UdPQX4WHut7wA/oBlUbT5/NfIzbmD7BCTPAB4HfCPJ5nb9we2+l7T/5V9Ck+jW0ySd71fVP1RVzcbSOh94b5I3AveqX2022x+Yd+jnqrqWppyOBX5JUwN7xsgh19M0EWmg7ANQH346svwKYDXwuKq6Lc3wzvt0cM3ZobPvoJ/f8wAnV9Wxd9nYdLr+KU2/wo/advqxP29V/cckX6AZ2Ov8JL/XJqlZP5/vHCPnuoVmIqEzkvw/mj6Vs9rd+7Tn0kBZA1Df7glc3375HwbMTnR/NvDiJPcFSHKfdvvNNPMB7OhKYF2S2c7aVwJfWUAcX6VJRiQ5CFjbnnOxzgJelOR+7Tnvk+SBNNMV/hS4Mcn9gWe3x3+njf8h7fqd7fFJHlJVl1XVu2mGRL5LH0Y17fh7JhmbBJIckuQ32uU9aJqeftCuB/h1miYiDZQJQH37W2AmyWXAq2i+CKmqK2juSPlKkm/RzAAGzTywb2s7Q2e/LKmqX9AMd/vJ9ly/BD6wgDjeD+zRfvZUmmGEFz3hTlVtAf4tcGbbvPUlYE01d+lc0v6cp9A078zGvxH4Qts8NDrfw5vbWzcvBW5j7qlAz6RpxhrnfsDnklwOXArcDvxNu+9xwIVzNC9pQBwOWtoNpZkj+ZiqeuUiP/+faWYVO2veg7ViWQOQdkNV9U3gnNlbTRfhcr/8ZQ1AkgbKGoAkDZQJQJIGygQgSQNlApCkgTIBSNJA/X/ecAenUOxumQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "#Now let's plot the PDF\n", "hist, bins = np.histogram(fair_results, bins = np.linspace(0,1,11), density=True)\n", "print(bins)\n", "#remember to set density=True for a PDF (integral must be 1)!\n", "#np.histogram returns bin edges :-/\n", "hist = np.concatenate([hist,[0.]])\n", "plt.plot(bins, hist, drawstyle='steps-post')\n", "plt.xlabel('Fraction of Heads (TS)')\n", "plt.ylabel('PDF of null distribution')\n", "plt.xlim(0,1)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great! So we've shown that we can use the fraction of heads in 1000 coin flips as a test statistic, since we've calculated the PDF under the null hypothesis. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Select a significance level at which the null hypothesis will be rejected\n", "\n", "### Type I and Type II Errors\n", "Since in testing our hypothesis, we have two options for our result (accept/reject the null hypothesis), and there are two options for the \"actual\" state of our null hypothesis (either it is true or it isn't), there are four possible outcomes describing whether we were correct or not. \n", "\n", "|* |$H_o$ rejected |Accept (Fail to reject) $H_o$ |\n", "|------------|-------------------------------|------------------------------------|\n", "|$H_o$ False |Correct |Type II error (false negative) |\n", "|$H_o$ True |Type I error (false positive) |Correct |\n", "\n", "If you reject a true null hypothesis, then this is a type II error. If you accept a false null hypothesis, then this is a type I error. \"Type I\" and \"type II\" are needlessly confusing names, so let's instead refer to these with the more intuitive names \"false positive\" and \"false negative\". \n", "\n", "For our coin flip example, if we reject the null hypothesis and state that the coin is unfair, but the coin was actually fair, then we got a _false positive_. \n", "\n", "If we were to accept the null hypthesis and state that the coin was fair, but the coin was actually unfair, then we got a _false negative_. \n", "\n", "Both types of errors are bound to happen sometimes, but we're typically more concerned with false positives: We don't want to make any claims of discovery when in fact it was simply a background fluctuation. We can control the rate at which this happens with the _significance_ of the test ($\\alpha$). The _significance_ refers to the probability of obtaining a false positive. \n", "\n", "Given a hypothesis test with some test statistic (where we can calculate the distribution of this test statistic under the null hypothesis, see above!), we can divide the range of our test statistic into the _acceptance region_, and the _rejection region_. If our observed test statistic falls in the _acceptance region_, we will accept the null hypothesis to be true. If our observed test statistic falls into the _rejection region_, we will reject the null hypothesis. \n", "\n", "The probability of a false positive is then simply the integral of our null hypothesis test statistic distribution over the rejection region. This is the _significance_ of the test. In principle, we can set the significance of our test to be whatever we want, simply by changing the acceptance and rejection regions. In practice, we typically want to set a high bar for announcing a result (this corresponds to a _low_ significance). In particle physics, the threshold for rejecting the null hypothesis is typically set such that the significance is $2.87 \\times 10 ^{-7}$ (the _one-sided five sigma threshold_). In other disciplines, the convention is different, with social scientists and biologists typically rejecting the null hypothesis at the $.01$ or $.05$ significance level. \n", "\n", "You might sometimes hear significance refered to in units of \"sigma\" (e.g. \"$5 \\sigma$\" or \"$3 \\sigma$\"). This can be obtained by the process shown in the figure below, by mapping the significance to a normal distribution. \n", "\n", "\n", "\n", "For our coin flip example, let's pretend to be social scientists and set our significance to be 0.1. We will be able to reject our null hypothesis at the 10% significance level. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFv1JREFUeJzt3Xm0ZWV55/Hvj0lUKKcqTQUtywm04lDiVYKaKNhtOwXUVhFFonGljKajojGtPcQhK93tsjHp7mhLtdqgBoNjgwMdlEGUBWghyFCKsRFtCB0KB0BUJp/+Y+8Lh7Lq3HNv3b3PvXd/P2uddfbeZ5+9n/uuc89z3vfd+31TVUiShme3aQcgSZoOE4AkDZQJQJIGygQgSQNlApCkgTIBSNJAmQAkaaBMAJI0UCYASRqoPbo8eJIrgRuB24Hbqmpm3P6rV6+u9evXdxmSlpGrfnQTAA+83z2nHIm0dF1wwQXXVdWahby30wTQOqSqrptkx/Xr17Nly5au49Ey8ZaPnAvAe44+eMqRSEtXkh8s9L02AUnSQHWdAAo4LckFSTbtaIckm5JsSbJl27ZtHYcjSZrVdQJ4alUdCDwb+OMkv7v9DlW1uapmqmpmzZoFNWNJkhag0wRQVVe3z9cCnwWe1OX5JEmT6ywBJLlnkn1nl4FnApd2dT5J0vx0eRXQA4DPJpk9z4lV9b87PJ8kaR46SwBVdQXwuK6OL0naNX3cByAtayee/0NOvujqaYfB4Rv342UHrZt2GFpBvA9AmsPJF13N1mtumGoMW6+5YUkkIa0s1gCkCWxYu4qTXjO9O5KPOO7cqZ1bK5c1AEkaKBOAJA2UCUCSBsoEIEkDZQKQpIEyAUjSQJkAJGmgTACSNFAmAEkaKBOAJA2UCUCSBsoEIEkDZQKQpIEyAUjSQJkAJGmgTACSNFAmAEkaKBOAJA2UCUCSBsoEIEkDZQKQpIEyAUjSQJkAJGmgTACSNFAmAEkaKBOAJA2UCUCSBsoEIEkDZQKQpIHqPAEk2T3JhUk+3/W5JEmT66MG8Abg2z2cR5I0D50mgCQPBJ4LfLDL80iS5q/rGsBfA38G/GpnOyTZlGRLki3btm3rOBxJ0qzOEkCS5wHXVtUF4/arqs1VNVNVM2vWrOkqHEnSdrqsATwFOCzJlcDfAYcm+ViH55MkzUNnCaCq3lZVD6yq9cBLgTOq6qiuzidJmh/vA5Ckgdqjj5NU1VnAWX2cS5I0GWsAkjRQJgBJGigTgCQN1ER9AEn2Ax48un9Vnd1VUJKk7s2ZAJK8GzgC2Arc3m4uwAQgScvYJDWA5wMHVNXNXQcjSerPJH0AVwB7dh2IJKlfk9QAfg5clOR04I5aQFW9vrOoJEmdmyQBnNI+JEkryJwJoKpOSLIXsH+76fKqurXbsCRJXZvkKqCnAycAVwIBHpTk970MVJKWt0magI4FnllVlwMk2R/4OPCELgOTJHVrkquA9pz98geoqu/iVUGStOxNUgPYkuSDwOxkLi8HtnQXkiSpD5MkgNcCfwzMXvb5VeD9nUUkSerFJFcB3Qy8t31IklaInSaAJJ+oqpckuYRm7J+7qKrHdhqZJKlT42oAb2ifn9dHIJKkfu30KqCquqZdfF1V/WD0Abyun/AkSV2Z5DLQf76Dbc9e7EAkSf0a1wfwWppf+g9LcvHIS/sC53QdmCSpW+P6AE4ETgX+I/DWke03VtWPO41KktS5nSaAqroeuD7Jv97upX2S7FNVP+w2NElSlya5EewLNJeBBtgbeAhwOfBbHcYlSerYJDeCPWZ0PcmBeBWQJC17k1wFdBdV9U3goA5ikST1aJL5AN40srobcCDwj51FJEnqxSR9APuOLN9G0yfw6W7CkST1ZZI+gHcCJFnVrNaNnUclSercnH0ASWbaAeEuBi5J8q0kM92HJknq0iRNQB+mGQ/oqwBJntpuczRQSVrGJrkK6PbZL3+AqvoaTV/AWEn2TvL1tsZwWZJ37kqgkqTFNW4soAPbxa8kOY5mIvgCjgDOmuDYNwOHVtXPkuwJfC3JqVV13i7GLElaBOOagI7dbv3tI8u/NkHM9qqqgJ+1q3u2jznfJ0nqx7ixgA7Z1YMn2R24AHg48L6qOn9XjylJWhzjmoCOqqqPbXcj2B2qas45gqvqdmBjknsDn03y6Kq6dLvzbAI2Aaxbt25ewUuSFm5cJ/A92+d9d/KYWFX9FDgTeNYOXttcVTNVNbNmzZr5HFaStAvGNQEd1zbh3FBVfzXfAydZA9xaVT9NcneamcXevfBQJUmLaexloG0TzpELPPZa4Mx2NrFvAF+qqs8v8FiSpEU2yY1g5yT5G+Ak4KbZje2ooDtVVRcDj9+18CRJXZkkAWxsn981sq2AQxc/HElSXyZJAK+uqitGNyR5aEfxSJJ6MslQEJ/awbZPLnYgkqR+jbsP4JE08/7eK8kLR15aRTM3sCRpGRvXBHQA8Dzg3sDvjWy/EfjDLoOSJHVv3H0AJwMnJzm4qs7tMSYJgH+64Wau+9nNHHHcdD9+W6+5gQ1rV001BqkLk/QBvCDJqiR7Jjk9ybYkR3UemQbvup/dzM9vmXPk8c5tWLuKwzfuN+0wpEU3yVVAz6yqP0vyAuBK4IXA2cDHugxMArjHXntw0msOnnYY0oo0SQ1gz/b5ucAnq+r6DuORJPVkkhrA55J8B/gF8Np2jJ9fdhuWJKlrc9YAquqtwJOBmaq6lWY4iMO7DkyS1K1x9wEcWlVnjN4DkGR0l890GZgkqVvjmoCeBpzBXe8BmFWYACRpWRt3H8Db2+dX9ReOJKkv45qAdjgV5KxJpoSUJC1d45qAZqd9PAB4InBKu/57wNe7DEqS1L1xTUDvBEhyNnBgVd3Yrr8D+EIv0UmSOjPJjWAPAG4ZWb+l3SZJWsYmuRHsI8DXk3y2XX8+cHxnEUmSejFnAqiqv0xyKvA77aZXVdWF3YYlSeraJDWA2Qngx04CL0laXibpA5AkrUAmAEkaKBOAJA3UuDuBb6QZ8+fXXgKqqpwjT5KWsXE3gu27s9ckScvfuBrAfce9sap+vPjhSJL6Mu4y0AtomoCyg9cKeGgnEUmSejGuCeghfQYiSerXnDeCJfndHW2vqrMXPxxJUl8muRP4LSPLewNPomkeOrSTiCRJvZhkLKC7TAmZ5EHAX3cWkSSpFwu5Eewq4FFz7ZTkQUnOTLI1yWVJ3rCAc0mSOjJJH8B/484bwnYDNjLZwHC3AW+uqm8m2Re4IMmXqmrrgqOVJC2aSfoAtows3wZ8vKrOmetNVXUNcE27fGOSbwP7ASYASVoCJukDOGFXT5JkPfB44PxdPZYkaXFM0gT0FOAdwIPb/WfHAproRrAk+wCfBt5YVTfs4PVNwCaAdevWTRy4JGnXTNIE9CHgGJpLP2+fz8GT7Enz5f+3VfWZHe1TVZuBzQAzMzM7GnxOktSBSRLA9VV16nwPnCQ0yePbVfXeeUcmSerUJAngzCTvAT4D3Dy7sZ0mcpynAK8ALklyUbvt31TVFxcUqSRpUU2SAA5qn2dGthVz3AlcVV9jxwPJSZKWgEmuAjqkj0AkSf1ySkhJGigTgCQN1E4TQJIXt8/OCyBJK9C4GsDb2udP9xGIJKlf4zqBf5TkNOAhSU7Z/sWqOqy7sCRJXRuXAJ4LHAh8FDi2n3AkSX0ZNyfwLcB5SZ5cVdvaMX2oqp/1Fp0kqTOTXAX0gCQXApcBW5NckOTRHcclSerYJAlgM/CmqnpwVa0D3txukyQtY5MkgHtW1ZmzK1V1FnDPziKSJPVikrGArkjy72k6gwGOAq7oLiRJUh8mqQH8AbCGZjTQTwOr222SpGVsksHgfgK8vodYJEk9ciwgSRooE4AkDZQJQJIGatxooKeNLL9tZ/tJkpancTWANSPLL+46EElSv8YlgOotCklS78ZdBvrQdhjojCzfweGgJWl5G5cADh9Z/s9dByJJ6te44aC/MrucZE27bVsfQUmSujfuKqAkeXuS64DLge8m2Zbkz/sLT5LUlXGdwMcATwWeWFX3rar7AAcBT0lyTC/RSZI6My4BvAI4sqq+P7uhqq6gGQ306K4DkyR1a1wC2LOqrtt+Y9sPsGd3IUmS+jAuAdyywNckScvAuMtAH5fkBpr7AODOG8MC7N1pVJKkzo27DHT3PgORJPVrpwkgyd7AHwEPBy4GPlxVt/UVmCSpW+P6AE4AZoBLgOcAx/YSkSSpF+P6ADZU1WMAknwI+Pp8Dpzkw8DzgGur6tELD1GS1IVxNYBbZxcW2PRzPPCsBbxPktSDSa4CgubKn7uPXBVUVbVq3IGr6uwk6xclSknSovMqIGmZ2HrNDRxx3LnTDoPDN+7Hyw5aN+0wtAjG1QB6kWQTsAlg3To/VNKOHL5xv2mHADRJCDABrBBTTwBVtRnYDDAzM+MsZNIOvOygdUviS3cp1EC0eMZ1AkuSVrDOEkCSjwPnAgckuSrJq7s6lyRp/jprAqqqI7s6tiRp19kEJEkDZQKQpIEyAUjSQJkAJGmgTACSNFAmAEkaqKnfCayl6cTzf8jJF1091Rh+fstt3GMvP6JSV6wBaIdOvujqO8Z9mZZ77LUHq/e521RjkFYyf15ppzasXcVJrzl4aud/y0ccd0bqkjUASRooE4AkDZRNQEvQUuiA3XrNDWxYO3bSN0nLnDWAJWgpdMBuWLtqyUxCIqkb1gCWqGl3wEpa+awBSNJAmQAkaaBMAJI0UCYASRooE4AkDZQJQJIGygQgSQNlApCkgTIBSNJAmQAkaaBMAJI0UCYASRooE4AkDZQJQJIGygQgSQPlfACS5mXrNTdwxHHnTjWGwzfux8sOWjfVGFYCE4CkiS2FWeJmZ8szAew6E4Ckib3soHVT/+Kddu1jJem0DyDJs5JcnuR7Sd7a5bkkSfPTWQJIsjvwPuDZwAbgyCQbujqfJGl+uqwBPAn4XlVdUVW3AH8HHN7h+SRJ89BlH8B+wP8dWb8KOGjcG67YdpPtezSdXBvWrpp2GNKStRSuRFoJpt4JnGQTsAlgn7UPm3I0S8OGtauWxNUW0lLk/8bi6TIBXA08aGT9ge22u6iqzcBmgJmZmTrpNQd3GJKk5W4pXIm0lHzijxb+3i77AL4BPCLJQ5LsBbwUOKXD80mS5qGzGkBV3ZbkXwF/D+wOfLiqLuvqfJKk+em0D6Cqvgh8sctzSJIWxsHgJGmgTACSNFAmAEkaKBOAJA2UCUCSBipVNe0Y7pDkRuDyacexRKwGrpt2EEuA5XAny+JOlsWdDqiqfRfyxqkPBbGdy6tqZtpBLAVJtlgWlsMoy+JOlsWdkmxZ6HttApKkgTIBSNJALbUEsHnaASwhlkXDcriTZXEny+JOCy6LJdUJLEnqz1KrAUiSetJ7Aphrovgkd0tyUvv6+UnW9x1jXyYoizcl2Zrk4iSnJ3nwNOLsw1xlMbLfv0xSSVbsFSCTlEWSl7SfjcuSnNh3jH2Z4H9kXZIzk1zY/p88Zxpx9iHJh5Ncm+TSnbyeJP+1LauLkxw450GrqrcHzbDQ/wd4KLAX8C1gw3b7vA74QLv8UuCkPmNcYmVxCHCPdvm1Qy6Ldr99gbOB84CZacc9xc/FI4ALgfu06/efdtxTLIvNwGvb5Q3AldOOu8Py+F3gQODSnbz+HOBUIMBvA+fPdcy+awCTTBR/OHBCu/wp4BlJ0mOMfZmzLKrqzKr6ebt6Hs2saivRJJ8LgL8A3g38ss/gejZJWfwh8L6q+glAVV3bc4x9maQsCpidQPtewD/2GF+vqups4Mdjdjkc+Eg1zgPunWTtuGP2nQB2NFH89hN83rFPVd0GXA/cr5fo+jVJWYx6NU12X4nmLIu2OvugqvpCn4FNwSSfi/2B/ZOck+S8JM/qLbp+TVIW7wCOSnIVzdwjf9JPaEvSfL9TltydwNqBJEcBM8DTph3LNCTZDXgv8Moph7JU7EHTDPR0mlrh2UkeU1U/nWpU03EkcHxVHZvkYOCjSR5dVb+admDLQd81gEkmir9jnyR70FTrftRLdP2apCxI8s+AfwscVlU39xRb3+Yqi32BRwNnJbmSpn3zlBXaETzJ5+Iq4JSqurWqvg98lyYhrDSTlMWrgU8AVNW5wN404wQN0UTfKaP6TgCTTBR/CvD77fKLgDOq7eFYYeYsiySPB46j+fJfqe28MEdZVNX1VbW6qtZX1Xqa/pDDqmrBY6AsYZP8j/wvml//JFlN0yR0RZ9B9mSSsvgh8AyAJI+iSQDbeo1y6TgFOLq9Gui3geur6ppxb+i1Cah2MlF8kncBW6rqFOBDNNW479F0eLy0zxj7MmFZvAfYB/hk2w/+w6o6bGpBd2TCshiECcvi74FnJtkK3A68papWXC15wrJ4M/A/khxD0yH8yhX6g5EkH6dJ/KvbPo+3A3sCVNUHaPpAngN8D/g58Ko5j7lCy0qSNAfvBJakgTIBSNJAmQAkaaBMAJI0UCYASRooE4AkDZQJQIsmye1JLhp5rF+EYz4/yYaR9Xe1d0d3JsnH2+F0j9lu+zuS/Ol2265sb8ba1XM+Pcnn57H/3ZN8JcnjRsr7x0m+3y5/Oclu7fDAlya5JMk3kjykff+Xk9xnV+PW8uZYQFpMv6iqjTt7Mcke7QB/8/F84PPAVoCq+vNdiG9OSX4DeGJVPbzL8yyCPwA+U1XfAjYCJDke+HxVfapdPxL4TeCxVfWrJA8Ebmrf/1Gaodf/su/AtXRYA1CnkrwyySlJzgBOT7JPmsltvtn+Kj18ZN+j21/e30ry0SRPBg4D3tP+qn1YkuOTvKjd/xlpJgK5pJ0s427t9iuTvHPkHI/cQVx7J/mf7esXJjmkfek0YL/2fL8zz7/1qCRfb997XJLd2+3/PcmWNJO3vHNk/2cl+U6SbwIvHNn+tJFf9Rcm2XcHp3s5cPIcIa0FrpkdGK2qrpodQppm2IAj5/P3aQWa9iQHPlbOg2ZYgovax2fbba+kGbzsvu36HsCqdnk1zW3rAX6LZlCz1e1rs/sfD7xo5BzH04wRtTfN0Lf7t9s/AryxXb4S+JN2+XXAB3cQ65tphhYAeCTNmDJ7A+vZ+YQb76AZXOuikcct7d/xKOBzwJ7tvu8Hjt7ub9kdOAt47Ej8j2j//k/Q/HqnPc5T2uV9gD22i2Mv4P/tIL7ty+qBbVlcBBwLPH67/f8BuN+0Pzc+pvewBqDF9Iuq2tg+XjCy/UtVNTuRRYD/kORi4Ms045U/ADgU+GRVXQcwsv/OHAB8v6q+266fQDNj0qzPtM8X0Hypb++pwMfac30H+AHNoGpz+auRv3Ejd05A8gzgCcA3klzUrj+0fe0l7a/8C2kS3QaapPP9qvqHqqrZWFrnAO9N8nrg3vXrzWargTmHfq6qq2jK6W3Ar2hqYM8Y2eVamiYiDZR9AOrDTSPLLwfWAE+oqlvTDO+8dwfnnB06+3b6+ZwHOKGq3naXjU2n65/S9Cv8pG2nH/v3VtV/SvIFmoG9zknyL9okNesXcx1j5Fg300wkdGqSf6LpUzm9fXnv9lgaKGsA6tu9gGvbL/9DgNmJ7s8AXpzkfgBJ7ttuv5FmPoDtXQ6sTzLbWfsK4CvziOOrNMmIJPsD69pjLtTpwIuS3L895n2TPJhmusKbgOuTPAB4drv/d9r4H9au39Een+RhVXVJVb2bZkjku/RhVNOOv3uSsUkgyYFJfrNd3o2m6ekH7XqA36BpItJAmQDUt78FZpJcAhxN80VIVV1Gc0XKV5J8i2YGMGjmgX1L2xk6+2VJVf2SZrjbT7bH+hXwgXnE8X5gt/a9J9EMI7zgCXeqaivw74DT2uatLwFrq7lK58L27zyRpnlnNv5NwBfa5qHR+R7e2F66eTFwKzueCvQ0mmasce4PfC7JpcDFwG3A37SvPQE4bwfNSxoQh4OWlqE0cyQfU1WvWOD7/wvNrGKnz7mzVixrANIyVFXfBM6cvdR0AS71y1/WACRpoKwBSNJAmQAkaaBMAJI0UCYASRooE4AkDdT/B3InCoOBauOJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def get_thresh(bgdata, sig=.1):\n", " #figure out the .05 significance threshold on our TS\n", " bgdata = sorted(bgdata)\n", " bgdata = np.array(bgdata)\n", " for val in bgdata:\n", " integral = float(len(bgdata[bgdata sig:\n", " lowerbound = val\n", " break\n", " return lowerbound\n", "\n", "lowerbound = get_thresh(fair_results, sig=0.1)\n", "#replot our pdf\n", "plt.plot(bins, hist, drawstyle='steps-post')\n", "plt.xlabel('Fraction of Heads (TS)')\n", "plt.ylabel('PDF of null distribution')\n", "\n", "#Plot the 5% significance thresholds as vertical lines.\n", "#This defines our acceptance region\n", "#If our observed TS value isn't above the vertical line, we can reject H0 at the 5% significance level. \n", "plt.axvline(lowerbound, color='steelblue')\n", "plt.xlim(0,1)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Compute the observed TS from data\n", "\n", "Now that we've defined our set of hypotheses, a test statistic, and a significance level at which we plan to be able to reject the null hypothesis, let's actually compute a TS for some observed data. " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "observed data array: ['H' 'T' 'T' 'T' 'T' 'T' 'T' 'T' 'H' 'T' 'H' 'H' 'H' 'T' 'T' 'T' 'T' 'T'\n", " 'H' 'T']\n", "observed test statistic: 0.3\n" ] } ], "source": [ "#use our functions to simulate some data\n", "obsdata = []\n", "for i in range(0,20):\n", " oneflip = flip_coin(pheads=.5)\n", " obsdata.append(oneflip)\n", " \n", "#...or actually flip a coin yourself! I used a biased coin for this data\n", "obsdata = ['H', 'T', 'T', 'T', 'T', 'T', 'T', 'T', 'H', 'T', 'H', 'H', 'H', 'T', 'T', 'T', 'T', 'T', 'H', 'T']\n", "obsdata = np.array(obsdata)\n", "fracheads = len(obsdata[obsdata=='H'])/len(obsdata)\n", "print(\"observed data array:\", obsdata)\n", "print(\"observed test statistic:\", fracheads)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Decide whether to accept/reject the null hypothesis\n", "\n", "Let's see where this observed test statistic lies relative to the bounds we set up in step 4:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHHtJREFUeJzt3Xu8VWW97/HPVy4uEdACKgW5pGCCwEqWuo3aoualdEtqRpaa7l4bt+4KTd1Wenbgzs7p5bHLOduOcpKDpiVpGeblbE28ixggKqCUeWsZJwELBJSbv/PHGAsmuJhrrMUaY661xvf9es3XHGPMMcf4rYfJ/M3necZ4HkUEZmZWPrvVOgAzM6sNJwAzs5JyAjAzKyknADOzknICMDMrKScAM7OScgIwMyspJwAzs5JyAjAzK6nueR5c0ivAW8AWYHNENFTbv3///jF06NA8Q7Ks1ixLnvseWLMQGletA2BQvz1rFoNZR7dgwYKVETGgLe/NNQGkjoqIlVl2HDp0KPPnz887HsvitxOS508+VLMQLr1pLgBXn31EzWIw6+gkvdrW97oJyMyspPJOAAHcJ2mBpMnN7SBpsqT5kuavWLEi53DMzKxJ3gng4xFxCPAp4F8k/f2OO0TE9IhoiIiGAQPa1IxlZmZtkGsfQES8nj6/IekO4DDgkTzPaWbta9OmTTQ2NvLOO+/UOpRSq6urY9CgQfTo0aPdjplbApC0J7BbRLyVLh8HXJnX+cwsH42NjfTp04ehQ4ciqdbhlFJEsGrVKhobGxk2bFi7HTfPJqAPAo9JegZ4Crg7Iv5vjuczsxy888479OvXz1/+NSSJfv36tXstLLcaQES8BIzN6/hmVhx/+ddeHv8GRdwHYNap/Wzea8xe9Hqtw2Bi/UC+cPjgWodhXYjvAzBrwexFr7N0+ZqaxrB0+ZoOkYRqpbGxkYkTJzJ8+HD2339/pkyZwsaNGwGYOXMmX/nKV2oc4Xv17t17u/VVq1ZRX19PfX09H/rQhxg4cODW9Y0bN3LVVVcxatQoxowZQ319PfPmzcs9RtcAzDIYuU9fZp1XuzuSJ10/t2bnrrWI4NRTT+X8889n9uzZbNmyhcmTJ3P55Zdz9dVX53LOzZs30717+3499uvXj0WLFgEwdepUevfuzSWXXALA3Llzueuuu1i4cCG77747K1eu3Jrg8uQagJl1aHPmzKGuro5zzz0XgG7duvGDH/yAGTNmsH79egD+9Kc/MWHCBIYPH860adMAWLduHSeeeCJjx47l4IMPZtasWQAsWLCAI488knHjxnH88cezfPlyACZMmMCFF15IQ0MDV111FUOGDOHdd9/deqz99tuPTZs28cc//pETTjiBcePG8YlPfIIXXngBgJdffpkjjjiC0aNHc8UVV7Tqb1y+fDn9+/dn9913B6B///7su+++u1hyLXMNwMyyW3Ah/HVR+x7zffUw7oc7fXnJkiWMGzduu219+/Zl8ODBvPjiiwA89dRTLF68mF69enHooYdy4okn8uqrr7Lvvvty9913A7B69Wo2bdrEV7/6VWbPns2AAQOYNWsWl19+OTNmzABg48aNW8cjW7hwIQ8//DBHHXUUd911F8cffzw9evRg8uTJXHfddQwfPpx58+ZxwQUXMGfOHKZMmcL555/P2WefzbXXXtuqIjjuuOO48sorGTFiBJ/85CeZNGkSRx55ZKuO0RauAZhZp3fsscfSr18/9thjD0499VQee+wxRo8ezf33389ll13Go48+yl577cWyZctYvHgxxx57LPX19XznO9+hsbFx63EmTZq03XJTreHWW29l0qRJrF27lieeeILTTz+d+vp6zjvvvK01iMcff5wzzjgDgLPOOqtV8ffu3ZsFCxYwffp0BgwYwKRJk5g5c+YulkrLXAMws+yq/FLPy8iRI7n99tu327ZmzRpee+01DjjgABYuXPieSyQlMWLECBYuXMg999zDFVdcwTHHHMMpp5zCqFGjmDu3+T6VPffcNvT4ySefzLe+9S3efPNNFixYwNFHH826devYe++9t7bl72hXLtXs1q0bEyZMYMKECYwePZobb7yRc845p83Hy8I1ADPr0I455hjWr1/PTTfdBMCWLVu4+OKLOeecc+jVqxcA999/P2+++SZvv/02v/71rxk/fjx//vOf6dWrF2eeeSaXXnopCxcu5MADD2TFihVbE8CmTZtYsmRJs+ft3bs3hx56KFOmTOGkk06iW7du9O3bl2HDhnHbbbcBSQf1M888A8D48eO59dZbAbjlllta9TcuW7aMP/zhD1vXFy1axJAhQ1p1jLZwAjCzDk0Sd9xxB7fddhvDhw9nxIgR1NXV8d3vfnfrPocddhinnXYaY8aM4bTTTqOhoYHnnnuOww47jPr6eqZNm8YVV1xBz549uf3227nssssYO3Ys9fX1PPHEEzs996RJk7j55pu3axq65ZZbuOGGGxg7diyjRo1i9uzZAPzoRz/i2muvZfTo0bz+eusu2V27di1f+tKXGDlyJGPGjGHp0qVMnTq1dQXVBoqI3E+SVUNDQ3hCmA7CE8Js1XQJZke4DLQWMTz//PMcdNBBhZ/X3qu5fwtJC1qabXFnXAMwMyspJwAzs5JyAjCzFnWkpuKyyuPfwAnAzKqqq6tj1apVTgI11DQfQF1dXbse1/cBmFlVgwYNorGxEc/ZXVtNM4K1JycAM6uqR48e7ToLlXUcbgIyMyspJwAzs5JyAjAzKyknADOzknICMDMrKScAM7OScgIwMyspJwAzs5JyAjAzKyknADOzknICMDMrKScAM7OScgIwMyspJwAzs5JyAjAzKyknADOzkso9AUjqJulpSXflfS4zM8uuiBrAFOD5As5jZmatkGsCkDQIOBH4SZ7nMTOz1su7BvBD4F+Bd3e2g6TJkuZLmu9Jp83MipNbApB0EvBGRCyotl9ETI+IhohoGDBgQF7hmJnZDvKsAYwHTpb0CnArcLSkm3M8n5mZtUJuCSAivhkRgyJiKPB5YE5EnJnX+czMrHV8H4CZWUl1L+IkEfEQ8FAR5zIzs2xcAzAzKyknADOzknICMDMrqUx9AJIGAkMq94+IR/IKyszM8tdiApD0PWASsBTYkm4OwAnAzKwTy1ID+AxwYERsyDsYMzMrTpY+gJeAHnkHYmZmxcpSA1gPLJL0ALC1FhARX8stKjMzy12WBHBn+jAzsy6kxQQQETdK6gmMSDcti4hN+YZlZmZ5y3IV0ATgRuAVQMB+kr7ky0DNzDq3LE1A1wDHRcQyAEkjgJ8D4/IMzMzM8pXlKqAeTV/+ABHxe3xVkJlZp5elBjBf0k+ApslcvgjMzy8kMzMrQpYEcD7wL0DTZZ+PAj/OLSIzMytElquANgDfTx9mZtZF7DQBSPpFRHxO0nMkY/9sJyLG5BqZmZnlqloNYEr6fFIRgZiZWbF2ehVQRCxPFy+IiFcrH8AFxYRnZmZ5yXIZ6LHNbPtUewdiZmbFqtYHcD7JL/39JT1b8VIf4PG8AzMzs3xV6wP4GXAv8F+Bb1Rsfysi3sw1KjMzy91OE0BErAZWS7psh5d6S+odEa/lG5qZmeUpy41gd5NcBiqgDhgGLANG5RiXmZnlLMuNYKMr1yUdgq8CMjPr9LJcBbSdiFgIHJ5DLGZmVqAs8wF8vWJ1N+AQ4M+5RWRmZoXI0gfQp2J5M0mfwC/zCcfMzIqSpQ9gGoCkvslqvJV7VGZmlrsW+wAkNaQDwj0LPCfpGUkN+YdmZmZ5ytIENINkPKBHASR9PN3m0UDNzDqxLFcBbWn68geIiMdI+gKqklQn6am0xrBE0rRdCdTMzNpXtbGADkkXH5Z0PclE8AFMAh7KcOwNwNERsVZSD+AxSfdGxJO7GLOZmbWDak1A1+yw/u2K5fdMELOjiAhgbbraI320+D4zMytGtbGAjtrVg0vqBiwADgCujYh5u3pMMzNrH9WagM6MiJt3uBFsq4hocY7giNgC1EvaG7hD0sERsXiH80wGJgMMHjy4VcGbmVnbVesE3jN97rOTR2YR8TfgQeCEZl6bHhENEdEwYMCA1hzWzMx2QbUmoOvTJpw1EfGD1h5Y0gBgU0T8TdIeJDOLfa/toZqZWXuqehlo2oRzRhuPvQ/wYDqb2O+A+yPirjYey8zM2lmWG8Eel/QfwCxgXdPGdFTQnYqIZ4GP7lp4ZmaWlywJoD59vrJiWwBHt384ZmZWlCwJ4MsR8VLlBkkfzikeMzMrSJahIG5vZttt7R2ImZkVq9p9AB8hmfd3L0mnVrzUl2RuYDMz68SqNQEdCJwE7A38Q8X2t4B/yjMoMzPLX7X7AGYDsyUdERFzC4zJDIC/rNnAyrUbmHR9bT9+S5evYeQ+fWsag1kesvQBnCKpr6Qekh6QtELSmblHZqW3cu0G1m9sceTx3I3cpy8T6wfWOgyzdpflKqDjIuJfJZ0CvAKcCjwC3JxnYGYAvXp2Z9Z5R9Q6DLMuKUsNoEf6fCJwW0SszjEeMzMrSJYawG8kvQC8DZyfjvHzTr5hmZlZ3lqsAUTEN4CPAQ0RsYlkOIiJeQdmZmb5qnYfwNERMafyHgBJlbv8Ks/AzMwsX9WagI4E5rD9PQBNAicAM7NOrdp9AN9On88tLhwzMytKtSagZqeCbJJlSkgzM+u4qjUBNU37eCBwKHBnuv4PwFN5BmVmZvmr1gQ0DUDSI8AhEfFWuj4VuLuQ6MzMLDdZbgT7ILCxYn1jus3MzDqxLDeC3QQ8JemOdP0zwMzcIjIzs0K0mAAi4ipJ9wKfSDedGxFP5xuWmZnlLUsNoGkC+KqTwJuZWeeSpQ/AzMy6ICcAM7OScgIwMyupancCv0Uy5s97XgIiIjxHnplZJ1btRrA+O3vNzMw6v2o1gPdXe2NEvNn+4ZiZWVGqXQa6gKQJSM28FsCHc4nIzMwKUa0JaFiRgZiZWbFavBFM0t83tz0iHmn/cMzMrChZ7gS+tGK5DjiMpHno6FwiMjOzQmQZC2i7KSEl7Qf8MLeIzMysEG25EawROKilnSTtJ+lBSUslLZE0pQ3nMjOznGTpA/ifbLshbDegnmwDw20GLo6IhZL6AAsk3R8RS9scrZmZtZssfQDzK5Y3Az+PiMdbelNELAeWp8tvSXoeGAg4AZiZdQBZ+gBu3NWTSBoKfBSYt6vHMjOz9pGlCWg8MBUYku7fNBZQphvBJPUGfglcGBFrmnl9MjAZYPDgwZkDNzOzXZOlCegG4CKSSz+3tObgknqQfPnfEhG/am6fiJgOTAdoaGhobvA5MzPLQZYEsDoi7m3tgSWJJHk8HxHfb3VkZmaWqywJ4EFJVwO/AjY0bUyniaxmPHAW8JykRem2b0XEPW2K1MzM2lWWBHB4+txQsS1o4U7giHiM5geSMzOzDiDLVUBHFRGImZkVy1NCmpmVlBOAmVlJ7TQBSDo9ffa8AGZmXVC1GsA30+dfFhGImZkVq1on8CpJ9wHDJN2544sRcXJ+YZmZWd6qJYATgUOAnwLXFBOOmZkVpdqcwBuBJyV9LCJWpGP6EBFrC4vOzMxyk+UqoA9KehpYAiyVtEDSwTnHZWZmOcuSAKYDX4+IIRExGLg43WZmZp1YlgSwZ0Q82LQSEQ8Be+YWkZmZFSLLWEAvSfovJJ3BAGcCL+UXkpmZFSFLDeAfgQEko4H+EuifbjMzs04sy2BwfwW+VkAsZmZWII8FZGZWUk4AZmYl5QRgZlZS1UYDva9i+Zs728/MzDqnajWAARXLp+cdiJmZFataAojCojAzs8JVuwz0w+kw0KpY3srDQZuZdW7VEsDEiuX/nncgZmZWrGrDQT/ctCxpQLptRRFBmZlZ/qpdBSRJ35a0ElgG/F7SCkn/Vlx4ZmaWl2qdwBcBHwcOjYj3R8T7gMOB8ZIuKiQ6MzPLTbUEcBZwRkS83LQhIl4iGQ307LwDMzOzfFVLAD0iYuWOG9N+gB75hWRmZkWolgA2tvE1MzPrBKpdBjpW0hqS+wBg241hAupyjcrMzHJX7TLQbkUGYmZmxdppApBUB/wzcADwLDAjIjYXFZiZmeWrWh/AjUAD8BzwaeCaQiIyM7NCVOsDGBkRowEk3QA81ZoDS5oBnAS8EREHtz1EMzPLQ7UawKamhTY2/cwETmjD+8zMrABZrgKC5MqfPSquCoqI6FvtwBHxiKSh7RKlmZm1O18FZNZJLF2+hknXz611GEysH8gXDh9c6zCsHVSrARRC0mRgMsDgwf5QmTVnYv3AWocAJEkIcALoImqeACJiOjAdoKGhwbOQmTXjC4cP7hBfuh2hBmLtp1onsJmZdWG5JQBJPwfmAgdKapT05bzOZWZmrZdbE1BEnJHXsc3MbNe5CcjMrKScAMzMSsoJwMyspJwAzMxKygnAzKyknADMzEqq5ncCW8f0lzUbWLl2A1fW8M7P9Rs306unP6JmeXENwJq1cu0G1m+s7QRwvXp2p3/v3Wsag1lX5p9XtlO9enZn1tlH1Oz8l97kcWfM8uQagJlZSTkBmJmVlJuAOqCfzXuN2Yter2kMl3RzB6xZV+caQAc0e9HrWyfeqBV3wJp1ff6J10GN3Kcvs86rXQcsv6065bOZdQGuAZiZlZQTgJlZSTkBmJmVlBOAmVlJOQGYmZWUE4CZWUk5AZiZlZQTgJlZSTkBmJmVlBOAmVlJOQGYmZWUE4CZWUk5AZiZlZQTgJlZSTkBmJmVlOcDMLNWWbp8DZOun1vTGCbWD+QLhw+uaQxdgROAmWU2sX5grUPYOlueE8CucwIws8y+cPjgmn/x1rr20ZXk2gcg6QRJyyS9KOkbeZ7LzMxaJ7cEIKkbcC3wKWAkcIakkXmdz8zMWifPGsBhwIsR8VJEbARuBSbmeD4zM2uFPPsABgJ/qlhvBA6v9oaXVqxz+x5JJ9fIffrWOgyzDqsjXInUFdS8E1jSZGAyQO999q9xNB3DyH36doirLcw6Iv/faD95JoDXgf0q1gel27YTEdOB6QANDQ0x67wjcgzJzDq7jnAlUkfyi39u+3vz7AP4HTBc0jBJPYHPA3fmeD4zM2uF3GoAEbFZ0leA/wS6ATMiYkle5zMzs9bJtQ8gIu4B7snzHGZm1jYeDM7MrKScAMzMSsoJwMyspJwAzMxKygnAzKykFBG1jmErSW8By2odRwfRH1hZ6yA6AJfDNi6LbVwW2xwYEX3a8saaDwWxg2UR0VDrIDoCSfNdFi6HSi6LbVwW20ia39b3ugnIzKyknADMzEqqoyWA6bUOoANxWSRcDtu4LLZxWWzT5rLoUJ3AZmZWnI5WAzAzs4IUngBamihe0u6SZqWvz5M0tOgYi5KhLL4uaamkZyU9IGlILeIsQktlUbHfaZJCUpe9AiRLWUj6XPrZWCLpZ0XHWJQM/0cGS3pQ0tPp/5NP1yLOIkiaIekNSYt38rok/Y+0rJ6VdEiLB42Iwh4kw0L/Efgw0BN4Bhi5wz4XANely58HZhUZYwcri6OAXuny+WUui3S/PsAjwJNAQ63jruHnYjjwNPC+dP0DtY67hmUxHTg/XR4JvFLruHMsj78HDgEW7+T1TwP3AgL+DpjX0jGLrgFkmSh+InBjunw7cIwkFRhjUVosi4h4MCLWp6tPksyq1hVl+VwA/DvwPeCdIoMrWJay+Cfg2oj4K0BEvFFwjEXJUhYBNE2gvRfw5wLjK1REPAK8WWWXicBNkXgS2FvSPtWOWXQCaG6i+B0n+Ny6T0RsBlYD/QqJrlhZyqLSl0mye1fUYlmk1dn9IuLuIgOrgSyfixHACEmPS3pS0gmFRVesLGUxFThTUiPJ3CNfLSa0Dqm13ykd7k5ga4akM4EG4Mhax1ILknYDvg+cU+NQOoruJM1AE0hqhY9IGh0Rf6tpVLVxBjAzIq6RdATwU0kHR8S7tQ6sMyi6BpBlovit+0jqTlKtW1VIdMXKUhZI+iRwOXByRGwoKLaitVQWfYCDgYckvULSvnlnF+0IzvK5aATujIhNEfEy8HuShNDVZCmLLwO/AIiIuUAdyThBZZTpO6VS0Qkgy0TxdwJfSpc/C8yJtIeji2mxLCR9FLie5Mu/q7bzQgtlERGrI6J/RAyNiKEk/SEnR0Sbx0DpwLL8H/k1ya9/JPUnaRJ6qcggC5KlLF4DjgGQdBBJAlhRaJQdx53A2enVQH8HrI6I5dXeUGgTUOxkonhJVwLzI+JO4AaSatyLJB0eny8yxqJkLIurgd7AbWk/+GsRcXLNgs5JxrIohYxl8Z/AcZKWAluASyOiy9WSM5bFxcD/lnQRSYfwOV30ByOSfk6S+PunfR7fBnoARMR1JH0gnwZeBNYD57Z4zC5aVmZm1gLfCWxmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiXlBGDtRtIWSYsqHkPb4ZifkTSyYv3K9O7o3Ej6eTqc7kU7bJ8q6ZIdtr2S3oy1q+ecIOmuVuy/h6SHJY2tKO83Jb2cLv9W0m7p8MCLJT0n6XeShqXv/62k9+1q3Na5eSwga09vR0T9zl6U1D0d4K81PgPcBSwFiIh/24X4WiTpQ8ChEXFAnudpB/8I/CoingHqASTNBO6KiNvT9TOAfYExEfGupEHAuvT9PyUZev2qogO3jsM1AMuVpHMk3SlpDvCApN5KJrdZmP4qnVix79npL+9nJP1U0seAk4Gr01+1+0uaKemz6f7HKJkI5Ll0sozd0+2vSJpWcY6PNBNXnaT/k77+tKSj0pfuAwam5/tEK//WMyU9lb73eknd0u3/S9J8JZO3TKvY/wRJL0haCJxasf3Iil/1T0vq08zpvgjMbiGkfYDlTQOjRURj0xDSJMMGnNGav8+6oFpPcuBH13mQDEuwKH3ckW47h2Twsven692Bvulyf5Lb1gWMIhnUrH/6WtP+M4HPVpxjJskYUXUkQ9+OSLffBFyYLr8CfDVdvgD4STOxXkwytADAR0jGlKkDhrLzCTemkgyutajisTH9Ow4CfgP0SPf9MXD2Dn9LN+AhYExF/MPTv/8XJL/eSY8zPl3uDXTfIY6ewP9rJr4dy2pQWhaLgGuAj+6w/x+AfrX+3PhRu4drANae3o6I+vRxSsX2+yOiaSILAd+V9CzwW5Lxyj8IHA3cFhErASr235kDgZcj4vfp+o0kMyY1+VX6vIDkS31HHwduTs/1AvAqyaBqLflBxd9Yz7YJSI4BxgG/k7QoXf9w+trn0l/5T5MkupEkSefliPhDRERTLKnHge9L+hqwd7y32aw/0OLQzxHRSFJO3wTeJamBHVOxyxskTURWUu4DsCKsq1j+IjAAGBcRm5QM71yXwzmbhs7eQjGfcwE3RsQ3t9uYdLpeQtKv8Ne0nb7q3xsR/03S3SQDez0u6fg0STV5u6VjVBxrA8lEQvdK+gtJn8oD6ct16bGspFwDsKLtBbyRfvkfBTRNdD8HOF1SPwBJ70+3v0UyH8COlgFDJTV11p4FPNyKOB4lSUZIGgEMTo/ZVg8An5X0gfSY75c0hGS6wnXAakkfBD6V7v9CGv/+6frW9nhJ+0fEcxHxPZIhkbfrw4ikHb+bpKpJQNIhkvZNl3cjaXp6NV0X8CGSJiIrKScAK9otQIOk54CzSb4IiYglJFekPCzpGZIZwCCZB/bStDO06cuSiHiHZLjb29JjvQtc14o4fgzslr53Fskwwm2ecCcilgJXAPelzVv3A/tEcpXO0+nf+TOS5p2m+CcDd6fNQ5XzPVyYXrr5LLCJ5qcCvY+kGauaDwC/kbQYeBbYDPxH+to44MlmmpesRDwctFknpGSO5Isi4qw2vv9HJLOKPdDiztZluQZg1glFxELgwaZLTdtgsb/8zTUAM7OScg3AzKyknADMzErKCcDMrKScAMzMSsoJwMyspP4/WrT0XFYxDi8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#replot our pdf\n", "plt.plot(bins, hist, drawstyle='steps-post')\n", "plt.xlabel('Fraction of Heads (TS)')\n", "plt.ylabel('PDF of null distribution')\n", "\n", "#Plot the 5% significance thresholds as vertical lines.\n", "#If our observed TS value isn't between these vertical lines, we can reject H0 at the 5% significance level. \n", "plt.axvline(lowerbound, color='steelblue')\n", "plt.axvline(fracheads, color='orange', label='Observed TS')\n", "plt.xlim(0, 1)\n", "plt.legend(loc=0)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since our observed TS was not within our acceptance region, we can reject the null hypothesis at the 5% level. \n", "(If you got that the TS was in the acceptance region, then you would instead fail to reject the null hypothesis)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### p-values\n", "\n", "An alternative way of doing this construction involves the calculation of a _p-value_. A _p-value_ is the probability, under the null hypothesis, of obtaining a test statistic at least as extreme as the one we have observed from data. If our p-value is smaller than our significance, then we can reject the null hypothesis. The p-value can be calculated by simply integrating our null hypothesis test statistic distribution. For the one-tailed case (above), this is: \n", "\n", "p-value = $\\int_0^{TS_{obs}} PDF_{H_o}$" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "one-tailed p-value: 0.019\n" ] } ], "source": [ "def get_p_value(obsval, bgdata):\n", " #calculate a 1 sided p-value for an observed value and a background distribution\n", " bgdata = np.array(bgdata)\n", " lowertail = float(len(bgdata[bgdata1/2$, we still can't reject the null hypothesis. \n", "\n", "We could have instead set up our hypothesis to be _two-tailed_:\n", "\n", "- $H_o$: The coin is fair, $P(heads)=1/2$\n", "\n", "- $H_a$: The coin is not fair, $P(heads) \\neq 1/2$\n", "\n", "Where now if either $P(heads) < 1/2$ OR $P(heads) > 1/2$, we will reject the null hypothesis. The procedure is similar to that of above, but we should be careful when setting up our acceptance and rejection regions. Since our hypothesis is now two tailed, our rejection region is _also_ two tailed: If the test statistic is either below the lower bound, or above the upper bound, we can reject the null hypothesis. " ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHJxJREFUeJzt3Xu8VWW97/HPV0CXCGgBFYJcUjBRYCVL3WYlal5Kt6RmZKnp7rVx467Q1G2lZwfu7JxeHrucs+0oJzloahKWYahna+JdxABRLkqZt5ZxtogFihcu/s4fYyyY0GKusS5jzLXW+L5fr/maY4w55hi/9ay15m8+zzPG8ygiMDOz8tml1gGYmVltOAGYmZWUE4CZWUk5AZiZlZQTgJlZSTkBmJmVlBOAmVlJOQGYmZWUE4CZWUn1zPPgkl4E3gC2AJsjoqHa/gMGDIjhw4fnGZJltX5V8txv/zYfonHtBgCG9N+jIyKydvLvo3tavHjxaxExsC3vzTUBpI6KiNey7Dh8+HAWLVqUdzyWxW8nJM+feqDNh7jkxgUAXHX24e2Px9rNv4/uSdJLbX2vm4DMzEoq7wQQwD2SFkua3NwOkiZLWiRp0Zo1a3IOx8zMmuSdAD4eEQcDnwb+WdInd9whImZERENENAwc2KZmLDMza4Nc+wAi4pX0+VVJtwOHAg/leU4z61ibNm2isbGRd955p9ahlFpdXR1DhgyhV69eHXbM3BKApD2AXSLijXT5OOCKvM5nZvlobGykb9++DB8+HEm1DqeUIoK1a9fS2NjIiBEjOuy4eTYBfRB4RNJTwBPAnRHxf3M8n5nl4J133qF///7+8K8hSfTv37/Da2G51QAi4nlgXF7HN7Pi+MO/9vL4HRRxH4BZl3bLwpeZu/SVWofBxPrBfPGwobUOw7oR3wdg1oK5S19h5er1NY1h5er1nSIJ1UpjYyMTJ05k5MiR7LvvvkydOpWNGzcCMGvWLL761a/WOMK/1adPn+3W165dS319PfX19XzoQx9i8ODBW9c3btzIlVdeyYEHHsjYsWOpr69n4cKFucfoGoBZBqMH9WP2ebW7g3bSdQtqdu5aiwhOPfVUpkyZwty5c9myZQuTJ0/msssu46qrrsrlnJs3b6Znz479eOzfvz9Lly4FYNq0afTp04eLL74YgAULFjBv3jyWLFnCbrvtxmuvvbY1weXJNQAz69Tmz59PXV0d5557LgA9evTghz/8ITNnzuStt94C4E9/+hMTJkxg5MiRTJ8+HYANGzZw4oknMm7cOA466CBmz54NwOLFiznyyCMZP348xx9/PKtXrwZgwoQJXHDBBTQ0NHDllVcybNgw3nvvva3H2meffdi0aRN//OMfOeGEExg/fjyf+MQnePbZZwF44YUXOPzwwxkzZgyXX355q37G1atXM2DAAHbbbTcABgwYwN57793OkmuZawBmlt3iC+AvSzv2mO+rh/E/2unLK1asYPz48dtt69evH0OHDuW5554D4IknnmD58uX07t2bQw45hBNPPJGXXnqJvffemzvvvBOAdevWsWnTJr72ta8xd+5cBg4cyOzZs7nsssuYOXMmABs3btw6HtmSJUt48MEHOeqoo5g3bx7HH388vXr1YvLkyVx77bWMHDmShQsXcv755zN//nymTp3KlClTOPvss7nmmmtaVQTHHXccV1xxBaNGjeJTn/oUkyZN4sgjj2zVMdrCNQAz6/KOPfZY+vfvz+67786pp57KI488wpgxY7j33nu59NJLefjhh9lzzz1ZtWoVy5cv59hjj6W+vp7vfve7NDY2bj3OpEmTtltuqjXceuutTJo0iTfffJPHHnuM008/nfr6es4777ytNYhHH32UM844A4CzzjqrVfH36dOHxYsXM2PGDAYOHMikSZOYNWtWO0ulZa4BmFl2Vb6p52X06NHcdttt221bv349L7/8Mvvttx9Lliz5m0skJTFq1CiWLFnCXXfdxeWXX84xxxzDKaecwoEHHsiCBc33qeyxx7ahsk8++WS+/e1v8/rrr7N48WKOPvpoNmzYwF577bW1LX9H7blUs0ePHkyYMIEJEyYwZswYbrjhBs4555w2Hy8L1wDMrFM75phjeOutt7jxxhsB2LJlCxdddBHnnHMOvXv3BuDee+/l9ddf5+233+bXv/41RxxxBH/+85/p3bs3Z555JpdccglLlixh//33Z82aNVsTwKZNm1ixYkWz5+3Tpw+HHHIIU6dO5aSTTqJHjx7069ePESNGMGfOHCDpoH7qqacAOOKII7j11lsBuPnmm1v1M65atYo//OEPW9eXLl3KsGHDWnWMtnACMLNOTRK33347c+bMYeTIkYwaNYq6ujq+973vbd3n0EMP5bTTTmPs2LGcdtppNDQ0sGzZMg499FDq6+uZPn06l19+Obvuuiu33XYbl156KePGjaO+vp7HHntsp+eeNGkSN91003ZNQzfffDPXX38948aN48ADD2Tu3LkA/PjHP+aaa65hzJgxvPJK6y7ZffPNN/nyl7/M6NGjGTt2LCtXrmTatGmtK6g2UETkfpKsGhoawhPCdBKeEGarpkswO8NloO2Joa2/j2eeeYYDDjigzee1jtPc70LS4pZmW9wZ1wDMzErKCcDMrKScAMysRZ2pqbis8vgdOAGYWVV1dXWsXbvWSaCGmuYDqKur69Dj+j4AM6tqyJAhNDY24jm7a6tpRrCO5ARgZlX16tWrQ2ehss7DTUBmZiXlBGBmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiXlBGBmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiXlBGBmVlJOAGZmJeUEYGZWUrknAEk9JD0paV7e5zIzs+yKqAFMBZ4p4DxmZtYKuSYASUOAE4Gf5nkeMzNrvbxrAD8C/gV4b2c7SJosaZGkRZ502sysOLklAEknAa9GxOJq+0XEjIhoiIiGgQMH5hWOmZntIM8awBHAyZJeBG4FjpZ0U47nMzOzVsgtAUTEtyJiSEQMB74AzI+IM/M6n5mZtY7vAzAzK6meRZwkIh4AHijiXGZmlo1rAGZmJeUEYGZWUk4AZmYllakPQNJgYFjl/hHxUF5BmZlZ/lpMAJK+D0wCVgJb0s0BOAGYmXVhWWoAnwX2j4h38w7GzMyKk6UP4HmgV96BmJlZsbLUAN4Clkq6D9haC4iIr+cWlZmZ5S5LArgjfZiZWTfSYgKIiBsk7QqMSjetiohN+YZlZmZ5y3IV0ATgBuBFQMA+kr7sy0DNzLq2LE1AVwPHRcQqAEmjgJ8D4/MMzMzM8pXlKqBeTR/+ABHxe3xVkJlZl5elBrBI0k+BpslcvgQsyi8kMzMrQpYEMAX4Z6Dpss+HgZ/kFpGZmRUiy1VA7wI/SB9mZtZN7DQBSPpFRHxe0jKSsX+2ExFjc43MzMxyVa0GMDV9PqmIQMzMrFg7vQooIlani+dHxEuVD+D8YsIzM7O8ZLkM9Nhmtn26owMxM7NiVesDmELyTX9fSU9XvNQXeDTvwMzMLF/V+gBuAe4G/ivwzYrtb0TE67lGZWZmudtpAoiIdcA6SZfu8FIfSX0i4uV8QzMzszxluRHsTpLLQAXUASOAVcCBOcZlZmY5y3Ij2JjKdUkH46uAzMy6vCxXAW0nIpYAh+UQi5mZFSjLfADfqFjdBTgY+HNuEZmZWSGy9AH0rVjeTNIn8Mt8wjEzs6Jk6QOYDiCpX7Iab+QelZmZ5a7FPgBJDemAcE8DyyQ9Jakh/9DMzCxPWZqAZpKMB/QwgKSPp9s8GqiZWReW5SqgLU0f/gAR8QhJX0BVkuokPZHWGFZImt6eQM3MrGNVGwvo4HTxQUnXkUwEH8Ak4IEMx34XODoi3pTUC3hE0t0R8Xg7YzYzsw5QrQno6h3Wv1Ox/DcTxOwoIgJ4M13tlT5afJ+ZmRWj2lhAR7X34JJ6AIuB/YBrImJhe49pZmYdo1oT0JkRcdMON4JtFREtzhEcEVuAekl7AbdLOigilu9wnsnAZIChQ4e2KngzM2u7ap3Ae6TPfXfyyCwi/grcD5zQzGszIqIhIhoGDhzYmsOamVk7VGsCui5twlkfET9s7YElDQQ2RcRfJe1OMrPY99seqpmZdaSql4GmTThntPHYg4D709nEfgfcGxHz2ngsMzPrYFluBHtU0r8Ds4ENTRvTUUF3KiKeBj7avvDMzCwvWRJAffp8RcW2AI7u+HDMzKwoWRLAVyLi+coNkj6cUzxmZlaQLENB3NbMtjkdHYiZmRWr2n0AHyGZ93dPSadWvNSPZG5gMzPrwqo1Ae0PnATsBfx9xfY3gH/MMygzM8tftfsA5gJzJR0eEQsKjMkMgFsWvszcpa/UOgxWrl7P6EH9ah2GWYfL0gdwiqR+knpJuk/SGkln5h6Zld7cpa+wcvX6WofB6EH9mFg/uNZhmHW4LFcBHRcR/yLpFOBF4FTgIeCmPAMzg+TDd/Z5h9c6DLNuKUsNoFf6fCIwJyLW5RiPmZkVJEsN4DeSngXeBqakY/y8k29YZmaWtxZrABHxTeBjQENEbCIZDmJi3oGZmVm+qt0HcHREzK+8B0BS5S6/yjMwMzPLV7UmoCOB+Wx/D0CTwAnAzKxLq3YfwHfS53OLC8fMzIpSrQmo2akgm2SZEtLMzDqvak1ATdM+7g8cAtyRrv898ESeQZmZWf6qNQFNB5D0EHBwRLyRrk8D7iwkOjMzy02WG8E+CGysWN+YbjMzsy4sy41gNwJPSLo9Xf8sMCu3iMzMrBAtJoCIuFLS3cAn0k3nRsST+YZlZmZ5y1IDaJoAvuok8GZm1rVk6QMwM7NuyAnAzKyknADMzEqq2p3Ab5CM+fM3LwEREZ4jz8ysC6t2I1jfnb1mZmZdX7UawPurvTEiXu/4cMzMrCjVLgNdTNIEpGZeC+DDuURkZmaFqNYENKLIQMzMrFgt3ggm6ZPNbY+Ihzo+HDMzK0qWO4EvqViuAw4laR46OpeIzMysEFnGAtpuSkhJ+wA/yi0iMzMrRFtuBGsEDmhpJ0n7SLpf0kpJKyRNbcO5zMwsJ1n6AP4n224I2wWoJ9vAcJuBiyJiiaS+wGJJ90bEyjZHa2ZmHSZLH8CiiuXNwM8j4tGW3hQRq4HV6fIbkp4BBgNOAGZmnUCWPoAb2nsSScOBjwIL23ssMzPrGFmagI4ApgHD0v2bxgLKdCOYpD7AL4ELImJ9M69PBiYDDB06NHPgZmbWPlmagK4HLiS59HNLaw4uqRfJh//NEfGr5vaJiBnADICGhobmBp8zM7McZEkA6yLi7tYeWJJIksczEfGDVkdmZma5ypIA7pd0FfAr4N2mjek0kdUcAZwFLJO0NN327Yi4q02RmplZh8qSAA5LnxsqtgUt3AkcEY/Q/EByZmbWCWS5CuioIgIxM7NieUpIM7OScgIwMyupnSYASaenz54XwMysG6pWA/hW+vzLIgIxM7NiVesEXivpHmCEpDt2fDEiTs4vLDMzy1u1BHAicDDwM+DqYsIxM7OiVJsTeCPwuKSPRcSadEwfIuLNwqIzM7PcZLkK6IOSngRWACslLZZ0UM5xmZlZzrIkgBnANyJiWEQMBS5Kt5mZWReWJQHsERH3N61ExAPAHrlFZGZmhcgyFtDzkv4LSWcwwJnA8/mFZGZmRchSA/gHYCDJaKC/BAak28zMrAvLMhjcX4CvFxCLmZkVyGMBmZmVlBOAmVlJOQGYmZVUtdFA76lY/tbO9jMzs66pWg1gYMXy6XkHYmZmxaqWAKKwKMzMrHDVLgP9cDoMtCqWt/Jw0GZmXVu1BDCxYvm/5x2ImZkVq9pw0A82LUsamG5bU0RQZmaWv2pXAUnSdyS9BqwCfi9pjaR/LS48MzPLS7VO4AuBjwOHRMT7I+J9wGHAEZIuLCQ6MzPLTbUEcBZwRkS80LQhIp4nGQ307LwDMzOzfFVLAL0i4rUdN6b9AL3yC8nMzIpQLQFsbONrZmbWBVS7DHScpPUk9wHAthvDBNTlGpWZmeWu2mWgPYoMxMzMirXTBCCpDvgnYD/gaWBmRGwuKjAzM8tXtT6AG4AGYBnwGeDqQiIyM7NCVOsDGB0RYwAkXQ880ZoDS5oJnAS8GhEHtT1EMzPLQ7UawKamhTY2/cwCTmjD+8zMrABZrgKC5Mqf3SuuCoqI6FftwBHxkKThHRKlmZl1OF8FZNZFrFy9nknXLWjz+9/4S/J9rj3HAJhYP5gvHja0XcewzqFaDaAQkiYDkwGGDvUflVlzJtYPrnUIQJKEACeAbqLmCSAiZgAzABoaGjwLmVkzvnjY0HZ/6F5yY/LN/6qzD2/zMdpbe7DOpVonsJmZdWO5JQBJPwcWAPtLapT0lbzOZWZmrZdbE1BEnJHXsc3MrP3cBGRmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiVV8zuBrXP6z/Xv8tqb73JFDceeWbl6PaMHVR1z0MzawTUAa9Zrb77LWxtrOwHc6EH9Os0YOGbdkWsAtlO9d+3J7HaMG9MRY8+YWX5cAzAzKyknADOzknITUCd0y8KXmbv0lZrGcHGPzfTe1X8eZt2ZawCd0Nylr2ydeKNWeu/akwF9dqtpDGaWL3/F66RGD+rH7PNq2Hn6W19+adbduQZgZlZSTgBmZiXlBGBmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiXlBGBmVlJOAGZmJeUEYGZWUk4AZmYl5QRgZlZSTgBmZiXlBGBmVlKeD8DMWmXl6vVMum5BTWOYWD+YLx42tKYxdAdOAGaW2cT6wbUOYetseU4A7ecEYGaZffGwoTX/4K117aM7ybUPQNIJklZJek7SN/M8l5mZtU5uCUBSD+Aa4NPAaOAMSaPzOp+ZmbVOnjWAQ4HnIuL5iNgI3ApMzPF8ZmbWCnn2AQwG/lSx3ggcVu0Nz6/Z4PY9kk6u0YP61ToMs06rM1yJ1B3UvBNY0mRgMkCfQfvWOJrOYfSgfp3iaguzzsj/Gx0nzwTwCrBPxfqQdNt2ImIGMAOgoaEhZp93eI4hmVlX1xmuROpMfvFPbX9vnn0AvwNGShohaVfgC8AdOZ7PzMxaIbcaQERslvRV4D+AHsDMiFiR1/nMzKx1cu0DiIi7gLvyPIeZmbWNB4MzMyspJwAzs5JyAjAzKyknADOzknICMDMrKUVErWPYStIbwKpax9FJDABeq3UQnYDLYRuXxTYui232j4i+bXljzYeC2MGqiGiodRCdgaRFLguXQyWXxTYui20kLWrre90EZGZWUk4AZmYl1dkSwIxaB9CJuCwSLodtXBbbuCy2aXNZdKpOYDMzK05nqwGYmVlBCk8ALU0UL2k3SbPT1xdKGl50jEXJUBbfkLRS0tOS7pM0rBZxFqGlsqjY7zRJIanbXgGSpSwkfT7921gh6ZaiYyxKhv+RoZLul/Rk+n/ymVrEWQRJMyW9Kmn5Tl6XpP+RltXTkg5u8aARUdiDZFjoPwIfBnYFngJG77DP+cC16fIXgNlFxtjJyuIooHe6PKXMZZHu1xd4CHgcaKh13DX8uxgJPAm8L13/QK3jrmFZzACmpMujgRdrHXeO5fFJ4GBg+U5e/wxwNyDg74CFLR2z6BpAloniJwI3pMu3AcdIUoExFqXFsoiI+yPirXT1cZJZ1bqjLH8XAP8GfB94p8jgCpalLP4RuCYi/gIQEa8WHGNRspRFAE0TaO8J/LnA+AoVEQ8Br1fZZSJwYyQeB/aSNKjaMYtOAM1NFL/jBJ9b94mIzcA6oH8h0RUrS1lU+gpJdu+OWiyLtDq7T0TcWWRgNZDl72IUMErSo5Iel3RCYdEVK0tZTAPOlNRIMvfI14oJrVNq7WdKp7sT2Joh6UygATiy1rHUgqRdgB8A59Q4lM6iJ0kz0ASSWuFDksZExF9rGlVtnAHMioirJR0O/EzSQRHxXq0D6wqKrgFkmSh+6z6SepJU69YWEl2xspQFkj4FXAacHBHvFhRb0Voqi77AQcADkl4kad+8o5t2BGf5u2gE7oiITRHxAvB7koTQ3WQpi68AvwCIiAVAHck4QWWU6TOlUtEJIMtE8XcAX06XPwfMj7SHo5tpsSwkfRS4juTDv7u280ILZRER6yJiQEQMj4jhJP0hJ0dEm8dA6cSy/I/8muTbP5IGkDQJPV9kkAXJUhYvA8cASDqAJAGsKTTKzuMO4Oz0aqC/A9ZFxOpqbyi0CSh2MlG8pCuARRFxB3A9STXuOZIOjy8UGWNRMpbFVUAfYE7aD/5yRJxcs6BzkrEsSiFjWfwHcJyklcAW4JKI6Ha15IxlcRHwvyVdSNIhfE43/cKIpJ+TJP4BaZ/Hd4BeABFxLUkfyGeA54C3gHNbPGY3LSszM2uB7wQ2MyspJwAzs5JyAjAzKyknADOzknICMDMrKScAM7OScgKwDiNpi6SlFY/hHXDMz0oaXbF+RXp3dG4k/TwdTvfCHbZPk3TxDtteTG/Gau85J0ia14r9d5f0oKRxFeX9uqQX0uXfStolHR54uaRlkn4naUT6/t9Kel9747auzWMBWUd6OyLqd/aipJ7pAH+t8VlgHrASICL+tR3xtUjSh4BDImK/PM/TAf4B+FVEPAXUA0iaBcyLiNvS9TOAvYGxEfGepCHAhvT9PyMZev3KogO3zsM1AMuVpHMk3SFpPnCfpD5KJrdZkn4rnVix79npN++nJP1M0seAk4Gr0m+1+0qaJelz6f7HKJkIZFk6WcZu6fYXJU2vOMdHmomrTtL/SV9/UtJR6Uv3AIPT832ilT/rmZKeSN97naQe6fb/JWmRkslbplfsf4KkZyUtAU6t2H5kxbf6JyX1beZ0XwLmthDSIGB108BoEdHYNIQ0ybABZ7Tm57NuqNaTHPjRfR4kwxIsTR+3p9vOIRm87P3pek+gX7o8gOS2dQEHkgxqNiB9rWn/WcDnKs4xi2SMqDqSoW9HpdtvBC5Il18EvpYunw/8tJlYLyIZWgDgIyRjytQBw9n5hBvTSAbXWlrx2Jj+HAcAvwF6pfv+BDh7h5+lB/AAMLYi/pHpz/8Lkm/vpMc5Il3uA/TcIY5dgf/XTHw7ltWQtCyWAlcDH91h/z8A/Wv9d+NH7R6uAVhHejsi6tPHKRXb742IpoksBHxP0tPAb0nGK/8gcDQwJyJeA6jYf2f2B16IiN+n6zeQzJjU5Ffp82KSD/UdfRy4KT3Xs8BLJIOqteSHFT9jPdsmIDkGGA/8TtLSdP3D6WufT7/lP0mS6EaTJJ0XIuIPERFNsaQeBX4g6evAXvG3zWYDgBaHfo6IRpJy+hbwHkkN7JiKXV4laSKyknIfgBVhQ8Xyl4CBwPiI2KRkeOe6HM7ZNHT2For5OxdwQ0R8a7uNSafrxST9Cn9J2+mr/rwR8d8k3UkysNejko5Pk1STt1s6RsWx3iWZSOhuSf9J0qdyX/pyXXosKynXAKxoewKvph/+RwFNE93PB06X1B9A0vvT7W+QzAewo1XAcElNnbVnAQ+2Io6HSZIRkkYBQ9NjttV9wOckfSA95vslDSOZrnADsE7SB4FPp/s/m8a/b7q+tT1e0r4RsSwivk8yJPJ2fRiRtOP3kFQ1CUg6WNLe6fIuJE1PL6XrAj5E0kRkJeUEYEW7GWiQtAw4m+SDkIhYQXJFyoOSniKZAQySeWAvSTtDmz4siYh3SIa7nZMe6z3g2lbE8RNgl/S9s0mGEW7zhDsRsRK4HLgnbd66FxgUyVU6T6Y/5y0kzTtN8U8G7kybhyrne7ggvXTzaWATzU8Feg9JM1Y1HwB+I2k58DSwGfj39LXxwOPNNC9ZiXg4aLMuSMkcyRdGxFltfP+PSWYVu6/Fna3bcg3ArAuKiCXA/U2XmrbBcn/4m2sAZmYl5RqAmVlJOQGYmZWUE4CZWUk5AZiZlZQTgJlZSf1/KtUIM5l2DcQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def get_thresh_2tail(bgdata, sig=.1):\n", " #figure out the .05 significance threshold on our TS\n", " bgdata = sorted(bgdata)\n", " bgdata = np.array(bgdata)\n", " #make sure you adjust for the fact that you are now doing a two-tailed test!\n", " #half the significance goes into each tail\n", " twotailedsig = sig/2.\n", " for val in bgdata:\n", " integral = float(len(bgdata[bgdata twotailedsig:\n", " lowerbound = val\n", " break\n", " for val in bgdata[::-1]:\n", " integral = float(len(bgdata[bgdata>val]))/len(bgdata)\n", " if integral > twotailedsig:\n", " upperbound = val\n", " break\n", " return lowerbound, upperbound\n", "\n", "lowerbound_2t, upperbound_2t = get_thresh_2tail(fair_results, sig=0.1)\n", "#replot our pdf\n", "plt.plot(bins, hist, drawstyle='steps-post')\n", "plt.xlabel('Fraction of Heads (TS)')\n", "plt.ylabel('PDF of null distribution')\n", "\n", "#Plot the 5% significance thresholds as vertical lines.\n", "#If our observed TS value isn't between these vertical lines, we can reject H0 at the 5% significance level. \n", "plt.axvline(lowerbound_2t, color='steelblue')\n", "plt.axvline(upperbound_2t, color='steelblue')\n", "plt.axvline(fracheads, color='orange', label='Observed TS')\n", "plt.xlim(0, 1)\n", "plt.legend(loc=0)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So even if we had done the two-tailed test, we still would have concluded that the coin was unfair. Also, we should be careful to take both tails into account when calculating a p-value:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "two-tailed p-value: 0.034\n" ] } ], "source": [ "def get_p_value_2t(obsval, bgdata):\n", " #calculate a 2 sided p-value for an observed value and a background distribution\n", " bgdata = np.array(bgdata)\n", " upperval = np.median(bgdata)+abs(np.median(bgdata)-obsval)\n", " lowertail = float(len(bgdata[bgdataupperval]))/len(bgdata)\n", " return lowertail + uppertail\n", " \n", "print(\"two-tailed p-value:\", get_p_value_2t(fracheads, fair_results))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that our p-value for testing the two-tailed hypothesis is different from the one-tailed hypothesis. It's always important to keep track of the hypothesis that the p-value corresponds to. If we _only_ wanted to test the hypothesis that coin was biased to produce tails more often, then we should use the one-tailed test. If we instead were interested in either the case that the coin is biased in the favor of tails _or_ in favor of heads, then we should use the two-tailed test. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Trial Factors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each time you obtain the answer of an analysis from real data, you have done one _experimental trial_. It is common to perform multiple trials when working with a specific analysis technique.\n", "\n", "**Example**: In the neutrino sources group, we might test the probability that _one_ location is consistent with background expectation. If we were to apply this framework to a list of 30 source candidate locations, then this would be 30 trials.\n", "\n", "**Example**: Repeating an analysis with new cuts. We might count the number of neutrinos with E > 1 TeV in 1 year, then count the number of neutrinos with E > 10 TeV in 1 year. This would be 2 trials. \n", "\n", "**Example**: Repeating with different signal hypotheses. You might look for emission from blazars using both the hypothesis of steady emission and a flare hypothesis (2 trials)\n", "\n", "Trials can be both independent, or correlated:\n", "\n", "**Independent**: Search for a new source/signal in 1 year of data, repeat the search in the next year of data (separately). This would be 2 independent trials\n", "\n", "**Correlated**: Search for a new source/signal in 1 year of data, repeat the search with 2 years of data, including the first year. This would be 2 trials, but they would be correlated.\n", "\n", "Let's go back to our coin flip example to see why trials matter." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TS | p-value (pre-trial)\n", "0.35 | 0.058\n", "0.45 | 0.257\n", "0.65 | 0.886\n", "0.55 | 0.593\n", "0.45 | 0.257\n", "0.5 | 0.403\n", "0.5 | 0.403\n", "0.5 | 0.403\n", "0.55 | 0.593\n", "0.6 | 0.764\n", "0.5 | 0.403\n", "0.55 | 0.593\n", "0.35 | 0.058\n", "0.3 | 0.019\n", "0.5 | 0.403\n", "0.6 | 0.764\n", "0.4 | 0.129\n", "0.45 | 0.257\n", "0.5 | 0.403\n", "0.6 | 0.764\n", "0.35 | 0.058\n", "0.5 | 0.403\n", "0.65 | 0.886\n", "0.55 | 0.593\n", "0.55 | 0.593\n", "0.5 | 0.403\n", "0.55 | 0.593\n", "0.45 | 0.257\n", "0.55 | 0.593\n", "0.45 | 0.257\n" ] } ], "source": [ "#Simulate 30 trials, or ask the class to do it experimentally\n", "#We'll go back to using a 1-tailed test here. \n", "print(\"TS\", \"|\" ,\"p-value (pre-trial)\")\n", "N_trials_pvals = []\n", "Ntrials = 30\n", "for i in range(0,Ntrials):\n", " one_trial = do_one_trial(pheads=0.5, Ncoins=20)\n", " print(one_trial, \"|\", get_p_value(one_trial, fair_results))\n", " N_trials_pvals.append(get_p_value(one_trial, fair_results))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we simulate enough trials, we'll _eventually_ obtain a p-value < 0.1, and will thus _eventually_ reject the null hypothesis, even if the null hypothesis is true. This is bad. We adjust for this with **trial factors**. Let's call the probability of seeing a result in a single trial is the \"pre-trial probability\" ($P_{pre}$), and the probability of seeing a result after many trials the \"post-trial probability\" ($P_{post}$). For N independent trials, $P_{post}$ can be written as:\n", "\n", "$P_{post} = 1-(1-P_{pre})^N$ (The Sidak Correction)\n", "\n", "Let's look at how this affects our results above." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "best pre-trial p-value: 0.019\n" ] } ], "source": [ "#Find our best p-value from 30 trials. \n", "bestp = min(N_trials_pvals)\n", "print(\"best pre-trial p-value:\", bestp)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Trial corrected p-value: 0.4375677420703794\n" ] } ], "source": [ "def trialcorrect(pvalue, Ntrials):\n", " return 1.-(1.-pvalue)**Ntrials\n", "\n", "print(\"Trial corrected p-value:\", trialcorrect(bestp, Ntrials))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fact that we did 30 trials had a pretty significant effect here. After trial correcting, our p-value isn't less than .1 anymore, so we can't reject the null hypothesis. This makes sense, since our trials were all generated from a fair coin. \n", "\n", "Suppose that one of our trials had turned up a p-value of $2.8 \\times 10^{-7}$ (so a $5\\sigma$ result). How might our 30 trials have affected this?" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Trial corrected 5 sigma result: 8.399965894967387e-06\n", "Trial corrected significance: 4.303660134506316\n" ] } ], "source": [ "correctedp = trialcorrect(2.8*10**-7, Ntrials)\n", "\n", "#What's the new significance?\n", "import scipy.stats as stats\n", "newsig = stats.norm.ppf(1-correctedp)\n", "\n", "print(\"Trial corrected 5 sigma result:\", correctedp)\n", "print(\"Trial corrected significance:\", newsig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So our 30 independent trials bumped our result down from $5\\sigma$ to $4.3\\sigma$. This is still a pretty significant result! The probability of a false positive is $8.3\\times10^{-6}$. This is why we choose the 5-sigma threshold in particle physics. If a result is $>5\\sigma$, then we can be pretty sure we're properly rejecting the null hypothesis, even with order ~10 trials. \n", "\n", "In practice, we oftentimes don't have fully independent trials. A major challenge for a lot of a analyses is figuring out how to handle correlated trials in an efficient way. This is a difficult problem, and is often analysis-dependent, so we won't discuss it here. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "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.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }