{ "metadata": { "name": "", "signature": "sha256:945d8418bcc971d99c631cb77619afddcd7f50dc53b103c5edb0d3391c603a55" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Thresholding with false discovery rate" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Our normal warm up\n", "from __future__ import print_function, division\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Jean-Baptiste Poline and Matthew Brett**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The false discovery rate is a different *type* of correction than family-wise correction. Instead of controlling for the risk of *any tests* falsely being declared significant under the null hypothesis, FDR will control the *number of tests falsely declared significant as a proportion of the number of all tests declared significant*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A basic idea on how the FDR works is the following.\n", "\n", "We have got a large number of p values from a set of individual tests. These might be p values from tests on a set of brain voxels.\n", "\n", "We are trying to a find a p value threshold $\\theta$ to do a reasonable job of distinguishing true positive tests from true negatives. p values that are less than or equal to $\\theta$ are *detections* and $\\theta$ is a *detection threshold*. \n", "\n", "We want to choose a detection threshold that will only allow a small number of false positive detections.\n", "\n", "A *detection* can also be called a *discovery*; hence false discovery rate.\n", "\n", "For the FDR, we will try to find a p value within the family of tests (the set of p values), that we can use as a detection threshold.\n", "\n", "Let's look at the p value for a particular test. Let's say there are $N$ tests, indexed with $i \\in 1 .. N$. We look at a test $i$, and consider using p value from this test as a detection threshold; $\\theta = p(i)$. The expected number of false positives (FP) in N tests at this detection threshold would be:\n", "\n", "$$\n", "E(FP) = N p(i)\n", "$$\n", "\n", "For example, if we had 100 tests, and the particular p value $p(i)$ was 0.1, then the expected number of false positive detections, thresholding at 0.1, is 0.1 * 100 = 10." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take some data from a random normal distribution to illustrate:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.random.seed(42) # so we always get the same random numbers\n", "N = 100\n", "z_values = np.random.normal(size=N)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Turn the Z values into p values:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import scipy.stats as sst\n", "normal_distribution = sst.norm(loc=0,scale=1.) #loc is the mean, scale is the variance.\n", "# The normal CDF\n", "p_values = normal_distribution.cdf(z_values)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make it easier to show, we sort the p values from smallest to largest:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "p_values = np.sort(p_values)\n", "i = np.arange(1, N+1) # the 1-based i index of the p values, as in p(i)\n", "plt.plot(i, p_values, '.')\n", "plt.xlabel('$i$')\n", "plt.ylabel('p value')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 4, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAERCAYAAACQIWsgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHBJJREFUeJzt3X1wVNX9x/HP2sRBrTUqmEg2EiUxDwQCNEjVsb8gpWEy\nECwCDU6tRZowONTScWKn7bSF6ciD/tHapnWgI2CVhIi2A+OQrYl0O1UJDwXFGqqRErpEgxMxTUEg\nsN7fH9skm91sHjZ7d+/ufb9mmNy7ubl7cgf2y/me8z3HYRiGIQAA/FwR6wYAAKyH4AAACEJwAAAE\nITgAAIIQHAAAQQgOAIAgpgaHhx9+WKmpqZo8eXLIax599FFlZ2ersLBQR44cMbM5AIBhMjU4LFu2\nTC6XK+T39+zZow8++EAtLS3avHmzVq5caWZzAADDZGpwuOeee3T99deH/P7u3bv10EMPSZJmzpyp\nzs5OnT592swmAQCGIaZjDm1tbcrIyOg9dzqdOnXqVAxbBACQLDAgHbh6h8PhiFFLAAA9kmL55unp\n6fJ4PL3np06dUnp6etB1BAwACE+4y+fFtOdQVlamP/zhD5KkpqYmpaSkKDU1dcBrDcPgj2Ho5z//\neczbYJU/PAueBc9i8D+jYWrPYenSpfrrX/+qjo4OZWRkaO3atbp06ZIkacWKFSotLdWePXuUlZWl\na665Rlu3bjWzOQCAYTI1ONTW1g55TXV1tZlNAACEIeYD0hiZ4uLiWDfBMngWfXgWfXgWkeEwRpuY\nigKHwzHq/BkA2M1oPjvpOQBAAqqsHN3PExwAIAG9//7ofp7gAAAJ6OqrR/fzjDkAQALq7JSuvz78\nz06CAwAkKAakAcDmKiul4mKptNTXaxgteg4AEEcqK32DzVdfLY0bJ5086Tvu6pLeeMN3zeLF0osv\nju6zM6YL7wEABucfDGpqfMd//avve2PHSh0dvuO0NN/XoiJp8+bRvy/BAQAsxj8g+PcIKiv7ZiEV\nFUkpKVJjo+/4pZekqipfYEhJGX0bSCsBgAWECghpaVJ7uy8ANDT0XdvTO+g5HiggjOazk+AAABZQ\nXNyXLvIPCKPpETBbCQDinH+6qKnJN6jc0CBNmOAbXI5Eqmgk6DkAQAwEDjT3vBapMQOJtBIAxB3/\nNFLP1NNII60EAHHGP40UiamnkUbPAQBioLMz8mmkQKSVAABBSCsBACKK4AAACEJwAAAEITgAAIIQ\nHAAAQQgOAIAgBAcAQBCCAwAgCMEBAKIk0vs8m4kKaQAwUahNfMxabM8fe0gDgEX57/kc6X2ezURa\nCQAiKDB1FGoTn2hv3jNSpJUAYJQGSx1t3mz+6quhsCorAERZqIDgv/9zrHsIrMoKAFHWM5ZQXy8d\nP+57Ld5SR4MhOADAMAx3LGHCBN8spHgODBJpJQAIyapjCcNl2bSSy+VSbm6usrOztXHjxqDvd3R0\naO7cuZo6daoKCgq0bds2M5sDACMSKnXUExASoYcQimnBwev1atWqVXK5XGpublZtba2OHTvW75rq\n6mpNmzZNb731ltxutx577DFdvnzZrCYBwIjE6zTUSDAtOBw4cEBZWVnKzMxUcnKyysvLtWvXrn7X\n3Hzzzerq6pIkdXV16cYbb1RSEnV5AKyhpibxxhKGy7RP4ra2NmVkZPSeO51O7d+/v981FRUVuvfe\nezV+/Hj997//1Ytm15IDwAj0pI7syLTg4HA4hrxm3bp1mjp1qtxut44fP645c+bo7bff1rXXXht0\n7Zo1a3qPi4uLVVxcHMHWArAr/0Hnmhrp8cf7n8dTT8HtdsvtdkfkXqYFh/T0dHk8nt5zj8cjp9PZ\n75o333xTP/nJTyRJEydO1K233qr33ntPRUVFQffzDw4AMBqhZiFVVkoff9y3FlJlZXz1HAL/47x2\n7dqw72XamENRUZFaWlrU2tqq7u5u1dXVqaysrN81ubm5amxslCSdPn1a7733nm677TazmgQAkgaf\nheQ/CG31xfHMZFpwSEpKUnV1tUpKSpSfn69vfvObysvL06ZNm7Rp0yZJ0o9//GMdOnRIhYWF+trX\nvqYnn3xSN9xwg1lNAgBJg89C8h+EjqeUUqRRBAcg4QWOK/S8ZuUCtkhg4T0AGERxcd84QjQ22bEK\nNvsBgAD+vYXkZN9rdh9HGAl6DgASkn9vYcEC6corEz+NFIieAwAE8B903rbNXkEhEug5AEhInZ32\nGHQeDAPSAIAgll2yGwAQnwgOABJG4G5tCB/BAUDC8F8Wo7Iy1q2JbwQHAAmDdZEihwFpAHHLrsti\nDBezlQDYkl2XxRguZisBsCXSSOah5wAgblHoNjjSSgCAIKSVANgGtQzRQXAAEFeoZYgOVmUFYHns\nzRB9jDkAsDz2ZggP+zkASCiBxW3szRB9BAcAltMzriD5AkVNDVNWo43gAMByAovbUlKofo42xhwA\nWA7FbZFBERwAIAhFcACAiCI4AACCEBwAAEEIDgCAIAQHAEAQggMAS2C1VWthKisAU/kvhTFunHTy\nZN+yGI8/3ve9ri7pjTd8P8OWn5FBnQMAS/EPCP4f+mPHSh0dvuPFi6WPP+5bJiMtTWpv91VFNzRQ\n/BYJLLwHIOoG6xH4r42Ulub7WlTk+8BvbOxbFuOBB/q+99JLUlUVVdFWQc8BwLANt0dw9qxvM57A\nD/2ee/QEAJbJMBdpJQCmCFw6+777Bk4D+fcIGhr6fpYP/diy7PIZLpdLubm5ys7O1saNGwe8xu12\na9q0aSooKFBxcbGZzQEwDP6zhpqb+2/J6b9aalOTr5fQ0CDt3Nl3nJLSt4oqgSF+mdZz8Hq9ysnJ\nUWNjo9LT0zVjxgzV1tYqLy+v95rOzk7dfffd+vOf/yyn06mOjg6NHTs2uJH0HICo8d91LXCQWKJH\nEE8s2XM4cOCAsrKylJmZqeTkZJWXl2vXrl39rqmpqdH9998vp9MpSQMGBgDRFap3QI/AXkwLDm1t\nbcrIyOg9dzqdamtr63dNS0uLzpw5o1mzZqmoqEjPP/+8Wc0BMAj/VNIzz/QFhAkTCAZ2ZdpUVofD\nMeQ1ly5d0uHDh/Xaa6/ps88+05133qmvfOUrys7ONqtZAAbgP/W0qooCNJgYHNLT0+XxeHrPPR5P\nb/qoR0ZGhsaOHaurrrpKV111lb761a/q7bffHjA4rFmzpve4uLiYwWtgFAJnIQVuy4n45Ha75Xa7\nI3MzwySXLl0ybrvtNuPEiRPGxYsXjcLCQqO5ubnfNceOHTNmz55tXL582Th37pxRUFBgvPvuu0H3\nMrGZgC393/8ZhuT7s3ixYXz6ad9XJI7RfHaa1nNISkpSdXW1SkpK5PV6tXz5cuXl5WnTpk2SpBUr\nVig3N1dz587VlClTdMUVV6iiokL5+flmNQnA/wT2FHoGmoEeFMEBNuGfSnrmGZaqsAMqpAEMKNRy\nF6x6ag+mL7z32WefyePxKCcnJ6w3ARAboRbAY9AZQxmyzmH37t2aNm2aSkpKJElHjhxRWVmZ6Q0D\nMHKBG+YMVtAGDGbItNL06dO1d+9ezZo1S0eOHJEkFRQU6B//+EdUGiiRVgKGy3/pi8WLfT0Elruw\nL1PTSsnJyUoJ+Ft1xRXsLgpYEbOQEClDfspPmjRJ27dv1+XLl9XS0qLvfe97uuuuu6LRNgAjVFND\n6giRMWRa6dy5c3riiSf06quvSpJKSkr005/+VGPGjIlKAyXSSgAQDqayAjYWuBQGPQb0MHXMYdas\nWQO+4d69e8N6QwCR5T9dtbKSMQZExpDB4amnnuo9vnDhgl5++WUlJZm26gaA//HvEYwbJ508GXzM\nonkwS1hppRkzZujgwYNmtGdApJVgR/7TUseOlTo6go+ZrorBmJpWOnPmTO/x559/rkOHDqmrqyus\nNwMQ2mDLaKekSI2NwcdMV4VZhgwO06dP7924JykpSZmZmXr22WdNbxhgN4FjBzU1fT2CntcCj+kp\nwCzMVgIsorRUqq/39QioU0AkmDKV9eWXXx50q8+FCxeG9YbhIDjADjo76REgskwJDt/5zncGDQ5b\nt24N6w3DQXAAgJGjCA4AEMT0/RxeeeUVNTc368KFC72v/exnPwvrDQH0oboZVjVkcFixYoXOnz+v\nvXv3qqKiQjt37tTMmTOj0TYgIYXanY3qZljJkGmlyZMn65133tGUKVN09OhRnT17VnPnztXrr78e\nrTaSVkJC8S9uS0uT2tuZoQRzjOazc8glu6+66ipJ0tVXX622tjYlJSWpvb09rDcD7Ijd2RCPhkwr\nzZs3T59++qmqqqr05S9/WZJUUVFhesOARDFYcRvVzbCqEc1WunDhgi5cuBC0M5zZSCshnlHchlgx\nNa00ZcoUrVu3TsePH9eYMWOiHhiAeMfubIhHQ/YcWltbVVdXpxdffFEOh0Pl5eVasmSJbrnllmi1\nkZ4D4g5TVGEFUSuCa2lp0S9+8Qtt375dXq83rDcMB8EB8cZ/RtLixYwrIDZML4Lz7z184Qtf0JNP\nPhnWmwF2wQY8iHdDBoeZM2equ7tbS5Ys0c6dO3XbbbdFo11A3PFPJT3zjFRVxSJ6iF9DppX++c9/\nKjc3N1rtGRBpJcQDUkmwGlPTSrEODIBVDbZzG6kkxLthjTkACDZUcRsQzwgOQJgCewpUOyORDDnm\ncP78ef3ud7/T66+/LofDoXvuuUcrV67UmDFjotVGxhxgSezcBqsztc5h8eLF+tKXvqRvfetbMgxD\nNTU1+s9//qOdO3eG9YbhIDgAwMiZGhzy8/PV3Nw85GtmIjjAKqh8RjwxdW2l6dOna9++fb3nTU1N\nvauzAnbTMwhdX+8LFECiGjI4HDp0SHfffbcmTJigzMxM3XXXXTp06JAmT56sKVOmDPqzLpdLubm5\nys7O1saNG0Ned/DgQSUlJemPf/zjyH8DIIqYrgq7GHK2ksvlCuvGXq9Xq1atUmNjo9LT0zVjxgyV\nlZUpLy8v6Lof/vCHmjt3LqkjWEJg6ujxx6l8hv0MGRwyMzPDuvGBAweUlZXV+/Pl5eXatWtXUHD4\nzW9+o0WLFungwYNhvQ8wXIN96I8bJ508OfC+zh9/3FfPUFXFdFXYg2l1Dm1tbcrIyOg9dzqd2r9/\nf9A1u3bt0t69e3Xw4EE5HA6zmgMEFa35f+iPHSt1dPiO09J8X3tSRw880P8csIMhxxzCNZwP+tWr\nV2vDhg29I+qklWCmwPEC//OpU/uOA/d1ZrMe2JFpPYf09HR5PJ7ec4/HI6fT2e+av//97yovL5ck\ndXR0qL6+XsnJySorKwu635o1a3qPi4uLVVxcbEq7kbgCl7fwP5dC7+tM5TPihdvtltvtjsi9RrTZ\nz0hcvnxZOTk5eu211zR+/Hjdcccdqq2tDRpz6LFs2TLNnz9fCxcuDG4kdQ4AMGKmb/YT1o2TklRd\nXa2SkhJ5vV4tX75ceXl52rRpkyRpxYoVZr01bM5/4Nl/oJmiNWD4TOs5RBI9BwQaLADcd9/AA83s\nsQC7sWTPATCT/8wj/wBQWdl/oDklRWpsZKYRMFIEB8QN/95CcrLvtVABINRAM4DhIa2EuOG/DeeC\nBdKVVxIAgMGQVoIt+KeLtm3rHwgYSwAii54D4gab6wAjY+p+DlZAcACAkSOthIRBjQJgDQQHWMpg\nU1QZVwCix7SF94BwhFoMjxoFILoYc4Cl+A86SwxAA6PBgDQAIMhoPjtJKwEAgjAgjagINQuJGUmA\nNREcEDGDTUMNNQuJGUmANREcEDHhrJTKqqmANREcMCqjXSnV/5iUEmAdzFbCiPgHg8CNdVgpFbAW\nprLCVP4BoatLeuMN3+uLF0tnz0r19b4eQkMDgQCwEoIDIi5UQEhLk9rb+4JBz7X0EADrITgg4vw3\n1vEPCC+9JFVVEQyAeEARHEatstIXEEpLfUtY+M8uamrypZAaGqQJE3zTTQkMQGKj5wBJ/XsKixf7\negaki4D4xn4OGDX/nkJPQKAgDbAvgoPNhKpifuYZxhIA9CE42EyoKuaqKnoKAPowIG0zbKYDYDgY\nkE5wgRXNPa9RxQwkPuocEFLgLCRSR4B9MFsJ/YRaDI/UEYDhoueQgPx7C/6L4ZE6AuyFngNC9ha2\nbSMoABg5eg5xarhLZxMYAPui52BD/vUKgTut0VsAMFoEhzgVuNyFxLRUAJFDWilOdXYSDAAMjjoH\nAEAQS+/n4HK5lJubq+zsbG3cuDHo+9u3b1dhYaGmTJmiu+++W0ePHjW7SQCAIZjac/B6vcrJyVFj\nY6PS09M1Y8YM1dbWKi8vr/eaffv2KT8/X9ddd51cLpfWrFmjpqam/o2k5yApeIYS6SQAg7Fsz+HA\ngQPKyspSZmamkpOTVV5erl27dvW75s4779R1110nSZo5c6ZOnTplZpPiWs8Mpfp6X6AAALOYOlup\nra1NGRkZvedOp1P79+8Pef2zzz6r0tJSM5sUd1gKA0AsmBocHA7HsK/9y1/+oi1btuiNN94Y8Ptr\n1qzpPS4uLlZxcfEoWxcf/OsZFizo28KTlBKAQG63W263OyL3MjU4pKeny+Px9J57PB45nc6g644e\nPaqKigq5XC5df/31A97LPzjYCcVtAIYr8D/Oa9euDftepo45FBUVqaWlRa2treru7lZdXZ3Kysr6\nXfPvf/9bCxcu1AsvvKCsrCwzmxOXamp8vYWGBgIDgOgxvc6hvr5eq1evltfr1fLly/WjH/1ImzZt\nkiStWLFC3/3ud/WnP/1Jt9xyiyQpOTlZBw4c6N9IG81WYkYSgEihCC6BsDkPgEix7FRWjNxAayYB\nQLTRc7AY1kwCECmklQAAQUgrAQAiiuBgAZWVvoHo0lJfWgkAYo3gYAGsmQTAatgJLkZYMwmAlTEg\nHSP+9QwLFkhXXskMJQCRNZrPTnoOMcKaSQCsjJ5DlAQui9HzGr0FAGahziEOsCwGgGgjrWRRDDoD\niFf0HEzEoDOAWKLnYBGB4woMOgOIVwSHUfIPCF1dUs8up5WVvgDBoDOAeERaaZT8U0dpaVJ7u6+n\nwM5tAGKNhfdiyD911NTElp4AEgM9h1Fi/wUAVkWdQ5SxzzOAeEBaKcpYRRVAoiM4hIF9ngEkOtJK\nw8C6SADiEWMOJghVv8C6SADiBRXSYQrsETz++MABIS3N95U0EgC7sHVw6BlYlnyB4uOP+xe0Sb6A\n8NJLUlUVaSQA9mHrAenAgeVQBW0TJvhSSQQGAHZh6zGHwAI2CtoAJBIGpAEAQSiCAwBElO2CQ2Wl\nbyXV0lJfGgkAEMwWaSVqFgDYEXUOQ/CfskrNAgAMzRZpJfZcAICRSZi0kn/qaNw46eRJ1kICYG+2\nncoaaixh7Fipo8N3zLgCALuy7FRWl8ul3NxcZWdna+PGjQNe8+ijjyo7O1uFhYU6cuTIiO7vv6/C\n8eO+14qKpKlT+44ZVwCAkTMtOHi9Xq1atUoul0vNzc2qra3VsWPH+l2zZ88effDBB2ppadHmzZu1\ncuXKEb1HqLGEnTsTd1zB7XbHugmWwbPow7Pow7OIDNOCw4EDB5SVlaXMzEwlJyervLxcu3bt6nfN\n7t279dBDD0mSZs6cqc7OTp0+fXrY71FTM/D6RykpibsWEn/x+/As+vAs+vAsIsO04NDW1qaMjIze\nc6fTqba2tiGvOXXqVMh7BhawJXIQAIBYMi04OByOYV0XOFgy2M+xdzMARIdpRXDp6enyeDy95x6P\nR06nc9BrTp06pfT09AHvFxg0du6Uhhl/Es7atWtj3QTL4Fn04Vn04VmMnmnBoaioSC0tLWptbdX4\n8eNVV1en2trafteUlZWpurpa5eXlampqUkpKilJTU4PuFQezbQEgoZgWHJKSklRdXa2SkhJ5vV4t\nX75ceXl52rRpkyRpxYoVKi0t1Z49e5SVlaVrrrlGW7duNas5AIARiIsiOABAdFl6baXhFNElKo/H\no1mzZmnSpEkqKCjQr3/9a0nSmTNnNGfOHN1+++36+te/rk4brTvu9Xo1bdo0zZ8/X5J9n0VnZ6cW\nLVqkvLw85efna//+/bZ9FuvXr9ekSZM0efJkPfDAA7p48aJtnsXDDz+s1NRUTZ48ufe1wX739evX\nKzs7W7m5uXr11VeHvL9lg8NwiugSWXJysn75y1/q3XffVVNTk37729/q2LFj2rBhg+bMmaP3339f\ns2fP1oYNG2Ld1Kh5+umnlZ+f3zs5wa7P4vvf/75KS0t17NgxHT16VLm5ubZ8Fq2trfr973+vw4cP\n65133pHX69WOHTts8yyWLVsml8vV77VQv3tzc7Pq6urU3Nwsl8ulRx55RJ9//vngb2BY1JtvvmmU\nlJT0nq9fv95Yv359DFsUWwsWLDAaGhqMnJwco7293TAMw/joo4+MnJycGLcsOjwejzF79mxj7969\nxrx58wzDMGz5LDo7O41bb7016HU7PotPPvnEuP32240zZ84Yly5dMubNm2e8+uqrtnoWJ06cMAoK\nCnrPQ/3u69atMzZs2NB7XUlJibFv375B723ZnsNwiujsorW1VUeOHNHMmTN1+vTp3hldqampI6oo\nj2c/+MEP9NRTT+mKK/r+ytrxWZw4cULjxo3TsmXLNH36dFVUVOjcuXO2fBY33HCDHnvsMd1yyy0a\nP368UlJSNGfOHFs+ix6hfvcPP/ywXynBcD5PLRschltEl+jOnj2r+++/X08//bSuvfbaft9zOBy2\neE6vvPKKbrrpJk2bNi3ktGa7PIvLly/r8OHDeuSRR3T48GFdc801QWkTuzyL48eP61e/+pVaW1v1\n4Ycf6uzZs3rhhRf6XWOXZzGQoX73oZ6LZYPDcIroEt2lS5d0//3368EHH9R9990nyfe/gfb2dknS\nRx99pJtuuimWTYyKN998U7t379att96qpUuXau/evXrwwQdt+SycTqecTqdmzJghSVq0aJEOHz6s\ntLQ02z2LQ4cO6a677tKNN96opKQkLVy4UPv27bPls+gR6t/ESAqOe1g2OPgX0XV3d6uurk5lZWWx\nblbUGIah5cuXKz8/X6tXr+59vaysTM8995wk6bnnnusNGols3bp18ng8OnHihHbs2KF7771Xzz//\nvC2fRVpamjIyMvT+++9LkhobGzVp0iTNnz/fds8iNzdXTU1NOn/+vAzDUGNjo/Lz8235LHqE+jdR\nVlamHTt2qLu7WydOnFBLS4vuuOOOwW8W6QGSSNqzZ49x++23GxMnTjTWrVsX6+ZE1d/+9jfD4XAY\nhYWFxtSpU42pU6ca9fX1xieffGLMnj3byM7ONubMmWN8+umnsW5qVLndbmP+/PmGYRi2fRZvvfWW\nUVRUZEyZMsX4xje+YXR2dtr2WWzcuNHIz883CgoKjG9/+9tGd3e3bZ5FeXm5cfPNNxvJycmG0+k0\ntmzZMujv/sQTTxgTJ040cnJyDJfLNeT9KYIDAASxbFoJABA7BAcAQBCCAwAgCMEBABCE4AAACEJw\nAAAEITgAAIIQHAAAQQgOwChcunRJS5cujXUzgIijQhoAEISeAwAgSFKsGwDEq3/961965ZVXNH78\neC1atCjWzQEiip4DEKb29nbdeOONunjxYqybAkQcYw7AKCxZskRbtmzRF7/4xVg3BYgoeg5AmLq6\nuuRwOPTOO+/EuilAxBEcgDB5vV6lpqaqu7s71k0BIo60EgAgCD0HAEAQggMAIAjBAQAQhOAAAAhC\ncAAABCE4AACCEBwAAEEIDgCAIP8Ph4+/nHvYTNsAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the (more or less) straight line of p value against $i$ index in this case, where there is no signal in the random noise." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to find a p value threshold $p(i)$ where there is only a small *proportion* of false positives among the detections. For example, we might accept a threshold such that 5% of all detections (discoveries) are likely to be false positives. If $d$ is the number of discoveries at threshold $\\theta$, and $q$ is the proportion of false positives we will accept (e.g. 0.05), then we want a threshold $\\theta$ such that $E(FP) / d < q$ where $E(x)$ is the expectation of $x$, here the number of FP I would get _on average_ if I was to repeat my experiment many times.\n", "\n", "So - what is $d$ in the plot above? Now that we have ordered the p values, for any index $i$, if we threshold at $\\theta \\le p(i)$ we will have $i$ detections ($d = i$). Therefore we want to find the largest $p(i)$ such that $E(FP) / i < q$. We know $E(FP) = N p(i)$ so we want the largest $p(i)$ such that:\n", "\n", "$$\n", " N p(i) / i < q \\implies p(i) < q i / N\n", "$$\n", "\n", "Let's take $q$ (the proportion of false discoveries = detections) as 0.05. We plot $q i / N$ (in red) on the same graph as $p(i)$ (in blue):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "q = 0.05\n", "plt.plot(i, p_values, 'b.', label='$p(i)$')\n", "plt.plot(i, q * i / N, 'r', label='$q i / N$')\n", "plt.xlabel('$i$')\n", "plt.ylabel('$p$')\n", "plt.legend()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAERCAYAAABhKjCtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHy1JREFUeJzt3Xlw1PX9x/FXmAQLaA2XiWSjqIk5iBxtgCpTjUYIk4Gg\nHBqdeiACo+MgVRFr/xA6lsPW1iMdB7WK1RLgp3biYNiGiMt0hCRoFBwTATHYJXJMxIAcmmT5/v7Y\nZrO72c21+90j+3zM7Ozudz/55pOP8nnv544zDMMQAABeBoQ7AwCAyESAAAD4RIAAAPhEgAAA+ESA\nAAD4RIAAAPhkeoC47777lJSUpGuuucZvmiVLlig9PV3jxo3Tp59+anaWAAA9YHqAmD9/vqxWq9/P\ny8vL9dVXX+nAgQN6+eWX9cADD5idJQBAD5geIH79619r6NChfj9/7733dM8990iSJk+erObmZh07\ndszsbAEAuhH2MYjGxkalpqa63lssFh0+fDiMOQIASBEQICTJe7ePuLi4MOUEANAuPtwZSElJkd1u\nd70/fPiwUlJSOqUjaABA3/R1y72wtyCKior0j3/8Q5JUVVWlxMREJSUl+UxrGAYPw9BTTz0V9jxE\nyoOyoCwoi64fgTC9BXHHHXdox44dampqUmpqqlauXKnW1lZJ0uLFi1VYWKjy8nKlpaVpyJAhev31\n183OEgCgB0wPEKWlpd2mKSkpMTsbAIBeCnsXE3ovLy8v3FmIGJRFB8qiA2URHHFGoJ1UIRIXFxdw\nfxoAxJpA6s6wz2ICgJ5gJmP3vAPBokWB3Y8AASBq0Ivgn68Aun9/YPdkDAIA+qnBgwP7ecYgAEQF\n6oCu+Sqf5mZp6NC+lxsBAkBUoA7omr/yCaTc6GICgH5i0SIpL08qLHS2HgJFgACAfqKuTtqxQ9q6\nNfAZTBIBAgBCpqGhwef1I0eO6OzZswHf/+BB53NurvTyywHfjgABAKHw9ddfq6qqyudnI0eO1DPP\nPBPw76iqkubNk7ZtkxITA74dg9QAokO01wHLly/X2rVr/X6+e/du1dfX6+677+7T/RmkBoAotGfP\nHlkslk7Xn3jiCVVUVEiSJk6cqMrKylBnrUuspAYAk23ZskW33HJLp+tr1qzxeD9y5Eh99dVXSktL\nC1XWukQLAkC/EIwpnsGeJtpu9+7dys7O7jbduHHj9MknnwTvFweIFgSAfmH/fucUT8lZ0W/eHNp7\n7NmzR5988on27dun6667TsePH9cFF1ygu+++W2fPnvXYK+nEiRPatm2b3nnnHW12+yVDhw7V/kA3\nUAoiWhAA+oX2fYcCmeIZyD2OHTumjIwMHTp0SLNmzdKdd96pp59+WpLkcDg80tbW1qqgoKDTtNdB\ngwappaWlb5k3AQECQL+wYUPgUzwDuce0adNUUVGhmTNnSpI+/fRTjRgxQpIUH+/ZWXPzzTdr/fr1\nuvfeez2unzx5UsOGDetb5k1AgADQLyQmOruEApn/H+g9KisrdcMNN0iS3njjDT322GOSpOTkZJ0+\nfdojbWlpqX7zm9+ovLzcde3IkSMRM0AtESAAIChOnjypEydOaPv27XrllVc0efJkzZ49W5J0ww03\nqKamxiP9lVdeqffff1+TJk1yXfvss880ZcqUkOa7KwxSA0AQbN++XUVFRbrnnns6fTZ79mz9+c9/\n1k033eS6Vlpa6pHmxx9/1M9//nP97Gc/Mz2vPUULAgAC9OWXX+ovf/mLjh8/rlOnTnX6PDExUSNG\njFBTU5Pfe2zcuFGLFy82M5u9xlYbAKJCtNcBhmHo1Vdf1cKFCzt9ZrfbVVtbq1mzZvX5/mZstUGA\nABAVqAO6xl5MAICQIUAAAHwiQAAAfCJAAAB8IkAAAHwiQAAAfCJAAAB8IkAAAHwiQAAAfCJAAECY\n/PDDD9q3b1+4s+EXAQIATPbEE0+ooqKi0/XNmzfrwgsvdL2vrq7WrbfeKovFora2NknOk+qKi4s1\nY8YM7dy5M2R5ltjuGwBMt2bNGp/X7Xa7UlJSXO8nT56s6dOn6+TJk3rnnXd0++23KykpSTNmzNCc\nOXM0aNCgUGVZUghaEFarVZmZmUpPT9fatWs7fd7U1KTp06dr/PjxysnJ0fr1683OEgCE3ZdffqnM\nzEyPa+fPn1dCQoKWLFmiF154wXX9zJkzIQ8OkskBwuFw6KGHHpLValVdXZ1KS0tVX1/vkaakpEQT\nJkzQZ599JpvNpkcffdTVtAKAXomLC86jj/bv368nn3xSFRUVevrpp7VhwwZt2rRJt912W6e0ZWVl\nrhPn2tXW1io3N1dFRUU6cuSIamtr//dn9T1PgTA1QNTU1CgtLU2jR49WQkKCiouLVVZW5pHm0ksv\ndR2wcerUKQ0fPrzTAd8A0COGEZxHH5w5c0Zz587VY489pmnTpumjjz7SuXPnVFBQoIaGBo+0DodD\nra2tGjhwoMf1vXv3auzYsRowYIAefPBBvfjii9q3b58yMjL6XCSBMLUmbmxsVGpqquu9xWJRdXW1\nR5qFCxfqpptu0qhRo/TDDz9o8+bNZmYJAEzx7rvvKicnR8OGDVNLS4uOHz+uBQsW6LnnntO9997r\nkbayslLTpk3rdI/z58+7Xt9///1KS0tTdna2Hn74YbOz75OpAaInzaJVq1Zp/PjxstlsOnjwoKZO\nnao9e/booosu6pR2xYoVrtd5eXnKy8sLYm4BoO+ampo0YcIESc7zqadMmSLJefZ0RUWFysvLVVhY\nKEmqqqrSU0895fHz3i2KxMREzZ07Vx9++KGWLVvW43zYbDbZbLYA/xonUwNESkqK7Ha7673dbpfF\nYvFIs3PnTv3+97+XJF111VW64oortG/fPuXm5na6n3uAAIBIUlxcrNWrV+v999/Xs88+q0ceeUSS\ndOWVV+r99993tRiam5s1dOhQj5/dvXu3Vq9ercGDBys/P981s2nJkiWduuW74/3leeXKlX3+m0w9\ncrStrU0ZGRn64IMPNGrUKE2aNEmlpaXKyspypXnkkUd08cUX66mnntKxY8f0y1/+Unv37tWwYcM8\nM8pxg0BMi6Y6YNKkSfrwww81ZMiQTp+98sormjlzppKTk4P6O6PuyNH4+HiVlJSooKBA2dnZuv32\n25WVlaV169Zp3bp1kqQnn3xSH3/8scaNG6ebb75ZzzzzTKfgAADR4MyZM1qxYoXsdrt27drlM823\n334b9OBgFlNbEMEUTd8eAARff6gDvv76a+3Zs0e33npr0O9tRguCAAEgKlAHdC3qupgAANGLAAEA\n8IkAAQDwiQABAPCJAAEA8Ild8QBEjXDtahqrCBAAokJPp2ouWiTt3y8NHixt2CAlJpqcsX6MLiYA\n/cr+/dKOHdLWrc5ggb4jQADoVwYPdj7n5kovvxzevEQ7VlIDiGreXUrt115+me4lia02AMSwvDxn\nl5IkzZsnceaYJ7baABCz6FIyDy0IAFGtuZkupa7QxQQA8IkuJgAxZdEi59hDYaGzBQFzECAARB3W\nOoQGK6kBRAX36awJCc5rDEybizEIAFHBfTrrrFnSwIEMTPdEIHUnLQgAEcl7AZz7dNb16wkMoUCA\nABCR2scZJGew2LCB6ayhRoAAEJG8F8AlJrJKOtQYgwAQkVgAFxwslAMA+MRCOQBA0BEgAAA+ESAA\nAD4RIAAAPhEgAAA+ESAARAx2aY0sTHMFYDr3bTNGjpS++aZjC43HH+/47NQp6aOPnD/D8aHBwToI\nABHHPSi4V/wjRkhNTc7X8+ZJx493bKmRnCwdPepcPb1tGwvkgoHN+gCERVctA/e9lJKTnc+5uc5K\nv7KyYwuNO+/s+Oztt6Vly1g9HSloQQDolZ62DE6fdh7o413xt9+jPQiwpYa56GICYBrvbbdvucV3\nl5B7y2Dbto6fpeIPr4jeasNqtSozM1Pp6elau3atzzQ2m00TJkxQTk6O8vLyzM4SgG64zyaqq/M8\n3tN9l9WqKmdrYds26f/+r+N1YmLH7qsEh+hlagvC4XAoIyNDlZWVSklJ0cSJE1VaWqqsrCxXmubm\nZk2ZMkX//ve/ZbFY1NTUpBEjRnTOKC0IIGTcT2/zHjiWaBlEk4htQdTU1CgtLU2jR49WQkKCiouL\nVVZW5pFmw4YNmjNnjiwWiyT5DA4AQstfK4GWQWwxNUA0NjYqNTXV9d5isaixsdEjzYEDB3TixAnd\neOONys3N1ZtvvmlmlgD44d6t9NJLHUHh8ssJCLHK1GmucXFx3aZpbW1VbW2tPvjgA509e1bXXnut\nfvWrXyk9Pd3MrAHw4j4tddkyFqnB5ACRkpIiu93uem+3211dSe1SU1M1YsQIDRo0SIMGDdL111+v\nPXv2+AwQK1ascL3Oy8tjQBsIgPfsJO8jPhGdbDabbDZbcG5mmKi1tdW48sorjYaGBuOnn34yxo0b\nZ9TV1Xmkqa+vN/Lz8422tjbjzJkzRk5OjvHFF190upfJWQVizg03GIbkfMybZxjff9/xjP4jkLrT\n1BZEfHy8SkpKVFBQIIfDoQULFigrK0vr1q2TJC1evFiZmZmaPn26xo4dqwEDBmjhwoXKzs42M1sA\n1LnF0D74DLRjoRwQQ9y7lV56iW0tYgErqQH45W9rDHZLjQ1s1gfAL3+b5jEQje5wYBDQz3gfutPV\nojegK3QxAf2M+zYZ8+Y5WwpsjRG76GIC4MLsJAQLLQign+F8BbhjFhMAwKeI3c0VQGh4D0wDwUCA\nAPqB9qms7Yf6AMHAIDUQwdwXuY0cKX3zTefXbLQHsxAggAjmvshtxAipqanz60WLnEGCgWkEG11M\nQATpapHb+PG+X7tPZSU4IJhoQQARxL3F4N0yaL/m/ZqgALMwzRWIIIWFzoHm3Fy2w0BwsA4C6CdY\n5IZgI0AAAHxioRwAIOi6DRDPPvus8vPzNWbMGD355JNqbW0NRb6AmMEqaESqbruYtmzZohkzZsgw\nDG3fvl07duzQH/7wh1Dlz4UuJvQnnPKGUDF1u++jR4+qvLxc119/vfLz83X27Nk+/SIAHTjlDdGg\n2wBht9vV3Nys119/Xd99953a2tp08uRJNTY2avny5aHIIxD13FsM3ltjvP22tGwZM5cQebrtYqqt\nrdW5c+c0ZcoUSdLBgwe1c+dOvfrqq9rR/hUoBOhiQjTjlDeES1imuR45ckSXXnppn35pXxAgEM1Y\nAIdwYR0EEOFYAIdwIUAAEch73IHAgHBgoRwQgTjEB9GOAAGYhEN8EO3Y7hsIIvdupZdeYvoqohsB\nAggi9wVwy5axKhrRjQABBKCrBXB0KyHaESCAAHR1AhzdSoh2BAggAN4thvazoYH+gHUQQABYAIdI\nx0I5AIBPLJQDQogDfhArCBBAL7FCGrHC9ABhtVqVmZmp9PR0rV271m+63bt3Kz4+Xu+++67ZWQIC\nwlRWxApTA4TD4dBDDz0kq9Wquro6lZaWqr6+3me65cuXa/r06YwzICJ4dyO5v3/pJeeZDmzbjf7O\n1GmuNTU1SktL0+jRoyVJxcXFKisrU1ZWlke6F198UXPnztXu3bvNzA7QaWHb4493vB85Uvrmm87n\nRC9aJB0/zgppxB5TA0RjY6NSU1Nd7y0Wi6qrqzulKSsr0/bt27V7927FxcWZmSXEOO+Fbe4V/4gR\nUlOT87X3OdF33un5HogFpnYx9aSyX7p0qdasWeOaikUXE8zkPX7g/n78+I7XVVWe3UgbNtCthNhj\nagsiJSVFdrvd9d5ut8tisXik+eSTT1RcXCxJampq0tatW5WQkKCioqJO91uxYoXrdV5envLy8kzJ\nN/ov760w3N9Lnp+5dyOxQhrRwmazyWazBeVepi6Ua2trU0ZGhj744AONGjVKkyZNUmlpaacxiHbz\n58/XzJkzNXv27M4ZZaEcAPRaIHWnqS2I+Ph4lZSUqKCgQA6HQwsWLFBWVpbWrVsnSVq8eLGZvx4x\nzH0w2n3wmaM/gZ5jqw1Era6CwC23+B58njePriLElohtQQBmcp+R5B4EFi3yHHxOTJQqK5mBBPQW\nAQJRxb3VkJDgvOYvCPgbfAbQM3QxIark5XW0GmbNkgYOJAgAXaGLCTHDveto/XrPYMDYAhBctCAQ\nVTigB+gdDgwCAPhEFxP6FdYwAJGBAIGI09X0VcYZgNDhRDlEHH8b6LGGAQgtxiAQcdwHoiUGpYFA\nMEgNAPApkLqTLiYAgE8MUiNk/M1OYqYSEJkIEAiqrqao+pudxEwlIDIRIBBUfdlhld1WgchEgEDA\nAt1h1f013UtA5GAWE3rNPSB4H87DDqtAZGGaK0znHhROnZI++sh5fd486fRpaetWZ0th2zaCARBJ\nCBAwhb+gkJwsHT3aERDa09JSACIPAQKmcD+cxz0ovP22tGwZAQGIBiyUQ1AsWuQMCoWFzu0u3Gcd\nVVU5u5O2bZMuv9w5FZXgAPRvtCDg4t5imDfP2UKg6wiIbpwHgaBwbzG0BwUWrQGxiwARg/ytdn7p\nJcYWAHQgQMQgf6udly2jxQCgA4PUMYgDeQD0BIPUMcB75XP7NVY7A/0f6yDQJe/ZSXQjAbGDWUzo\nxN8GenQjAegpWhD9lHurwX0DPbqRgNhCCwKS/Lca1q8nMADoPVoQUayn224THIDYRQsiRrmvZ/A+\nsY1WA4BAESCimPfWGBJTVgEED11MUay5mYAAoGusgwAA+BTx50FYrVZlZmYqPT1da9eu7fT5P//5\nT40bN05jx47VlClTtHfv3lBkCwDQBdNbEA6HQxkZGaqsrFRKSoomTpyo0tJSZWVludLs2rVL2dnZ\nuvjii2W1WrVixQpVVVV5ZpQWhKTOM5foWgLQlYhuQdTU1CgtLU2jR49WQkKCiouLVVZW5pHm2muv\n1cUXXyxJmjx5sg4fPmx2tqJW+8ylrVudwQIAzGL6LKbGxkalpqa63lssFlVXV/tN//e//12FhYVm\nZyuqsG0GgHAwPUDExcX1OO2HH36o1157TR999JHPz1esWOF6nZeXp7y8vABzFx3c1zvMmtVxHCjd\nSwC82Ww22Wy2oNzL9ACRkpIiu93uem+322WxWDql27t3rxYuXCir1aqhQ4f6vJd7gIglLIAD0FPe\nX55XrlzZ53uZPgaRm5urAwcO6NChQ2ppadGmTZtUVFTkkea///2vZs+erbfeektpaWlmZynqbNjg\nbDVs20ZwABA6IVkHsXXrVi1dulQOh0MLFizQ7373O61bt06StHjxYt1///3617/+pcsuu0ySlJCQ\noJqaGs+MxtAsJmYqAQgWFsr1MxzwAyBYInqaK3rP1x5LABBqtCAiEHssAQgWupgAAD7RxQQACDoC\nRIRYtMg5OF1Y6OxiAoBwI0BECPZYAhBpOFEujNhjCUAkY5A6jNzXO8yaJQ0cyMwlAMEVSN1JCyKM\n2GMJQCSjBRFC3ltotF+j1QDALKyDiBJsoQEg1OhiimAMRAOIVrQgTMZANIBwogURQbzHGRiIBhCt\nCBBB4B4UTp2S2k9MXbTIGSQYiAYQjehiCgL3bqTkZOnoUWeLgRPgAIQbm/WFmXs3UlUVx4MC6B9o\nQQQB5zcAiFSsgwgDzo0GEA3oYgoDdl8F0N8RIPqIc6MB9Hd0MfUQ+ygBiEaMQZjE3/oG9lECEC1Y\nSR0A75bB44/7DgrJyc5nupQAxIqYDxDtg82SM1gcP+656E1yBoW335aWLaNLCUDsiPlBau/BZn+L\n3i6/3NmtRHAAECtifgzCe5Ebi94A9CcMUgMAfGKQGgCiiWFIra3Sjz9K5851/9z+2t91f2l//DGg\nbMZkC4JtMgBIclbUbW2dK9vuKuTuKmpf173TDBgg/exnzsegQc6H9+uuPvP37JUuLiODLqbusKYB\niGC+vlH35Ftyd2m7ukf7c3cV9QUXdLzupjLutnJ3f44PTQcOXUw94D6dlTUNgB+9rai7qrh7+y28\nvaLuTSXr/jxsWM/ShamijkYxUzLu01dZ04CI19s+6r50hfSkou5JJetemQ8eLA0f3rtv1O2vqagj\nTr/qYnLvRho5UvrmG/ZOQoDc+6j78i05kG/hvflG3dOKvCfdJVTU/UpMT3P1N7YwYoTU1OR8zThD\nP9BVRd2Tb8l9qaj70vXRmy6S7gYZqagRBBE9BmG1WrV06VI5HA7df//9Wr58eac0S5Ys0datWzV4\n8GCtX79eEyZM6PH9/Y0tJCZKlZWMMwSdv1kffe0C6aqi9/WN2v0bsHtF211FPWRI564Pf5W59zMV\nNWKUqf/nOxwOPfTQQ6qsrFRKSoomTpyooqIiZWVludKUl5frq6++0oEDB1RdXa0HHnhAVVVVPf4d\n/sYWpP7bpWSz2ZR3ww3BG0zsbXeJrz5q90q2q37q9j7qvnwL91FR22w25eXlhf4/QgSiLDpQFsFh\naoCoqalRWlqaRo8eLUkqLi5WWVmZR4B47733dM8990iSJk+erObmZh07dkxJSUk9+h0bNngGAveu\npJB0K3XXR92bCtrXs3v6/12znTypvPPnO1fUPZnx0dU36gAq6nChIuhAWXSgLILD1H/pjY2NSk1N\ndb23WCyqrq7uNs3hw4f9Bghfi9w2b5azom7pYtZHbyvk3tzHu4/aVzdIdxV1V90l3pX0X/8q/eEP\nEVVRA+h/TK1h4uLiepTOewClq5/76dM6vfBxsQbpnNpG/SgN7mbWR08XtPiqqHs6yBjqinrgQIID\nANOZWsukpKTIbre73tvtdlksli7THD58WCkpKT7v1x44/tF+4dz/Hu3On5dOn3Y++rmVK1eGOwsR\ng7LoQFl0oCwCZ2qAyM3N1YEDB3To0CGNGjVKmzZtUmlpqUeaoqIilZSUqLi4WFVVVUpMTPTZvRQl\ns3EBoN8wNUDEx8erpKREBQUFcjgcWrBggbKysrRu3TpJ0uLFi1VYWKjy8nKlpaVpyJAhev31183M\nEgCgh6JmoRwAILQi/shRq9WqzMxMpaena+3ateHOTkjZ7XbdeOONGjNmjHJycvTCCy9Ikk6cOKGp\nU6fq6quv1rRp09Tc3BzmnIaOw+HQhAkTNHPmTEmxWxbNzc2aO3eusrKylJ2drerq6pgti9WrV2vM\nmDG65pprdOedd+qnn36KmbK47777lJSUpGuuucZ1rau/ffXq1UpPT1dmZqYqKiq6vX9EB4j2hXZW\nq1V1dXUqLS1VfX19uLMVMgkJCfrrX/+qL774QlVVVfrb3/6m+vp6rVmzRlOnTtX+/fuVn5+vNWvW\nhDurIfP8888rOzvbNWEhVsvi4YcfVmFhoerr67V3715lZmbGZFkcOnRIr7zyimpra/X555/L4XBo\n48aNMVMW8+fPl9Vq9bjm72+vq6vTpk2bVFdXJ6vVqgcffFDnz5/v+hcYEWznzp1GQUGB6/3q1auN\n1atXhzFH4TVr1ixj27ZtRkZGhnH06FHDMAzjyJEjRkZGRphzFhp2u93Iz883tm/fbsyYMcMwDCMm\ny6K5udm44oorOl2PxbL47rvvjKuvvto4ceKE0draasyYMcOoqKiIqbJoaGgwcnJyXO/9/e2rVq0y\n1qxZ40pXUFBg7Nq1q8t7R3QLwtciusbGxjDmKHwOHTqkTz/9VJMnT/ZYaZ6UlKRjx46FOXeh8dvf\n/lZ/+tOfNGBAx/+2sVgWDQ0NGjlypObPn69f/OIXWrhwoc6cOROTZTFs2DA9+uijuuyyyzRq1Cgl\nJiZq6tSpMVkW7fz97d9++63HMoOe1KcRHSB6utCuvzt9+rTmzJmj559/XhdddJHHZ3FxcTFRTlu2\nbNEll1yiCRMm+J3yHCtl0dbWptraWj344IOqra3VkCFDOnWhxEpZHDx4UM8995wOHTqkb7/9VqdP\nn9Zbb73lkSZWysKX7v727sologNETxba9Xetra2aM2eO7rrrLt1yyy2SnN8Kjh49Kkk6cuSILrnk\nknBmMSR27typ9957T1dccYXuuOMObd++XXfddVdMloXFYpHFYtHEiRMlSXPnzlVtba2Sk5Njriw+\n/vhjXXfddRo+fLji4+M1e/Zs7dq1KybLop2/fxO9WZTcLqIDhPtCu5aWFm3atElFRUXhzlbIGIah\nBQsWKDs7W0uXLnVdLyoq0htvvCFJeuONN1yBoz9btWqV7Ha7GhoatHHjRt1000168803Y7IskpOT\nlZqaqv3790uSKisrNWbMGM2cOTPmyiIzM1NVVVU6d+6cDMNQZWWlsrOzY7Is2vn7N1FUVKSNGzeq\npaVFDQ0NOnDggCZNmtT1zYI9YBJs5eXlxtVXX21cddVVxqpVq8KdnZD6z3/+Y8TFxRnjxo0zxo8f\nb4wfP97YunWr8d133xn5+flGenq6MXXqVOP7778Pd1ZDymazGTNnzjQMw4jZsvjss8+M3NxcY+zY\nscatt95qNDc3x2xZrF271sjOzjZycnKMu+++22hpaYmZsiguLjYuvfRSIyEhwbBYLMZrr73W5d/+\nxz/+0bjqqquMjIwMw2q1dnt/FsoBAHyK6C4mAED4ECAAAD4RIAAAPhEgAAA+ESAAAD4RIAAAPhEg\nAAA+ESAAAD4RIIAAtba26o477gh3NoCgYyU1AMAnWhAAAJ/iw50BIJp9/fXX2rJli0aNGqW5c+eG\nOztAUNGCAAJw9OhRDR8+XD/99FO4swIEHWMQQIBuu+02vfbaa7rwwgvDnRUgqGhBAAE4deqU4uLi\n9Pnnn4c7K0DQESCAADgcDiUlJamlpSXcWQGCji4mAIBPtCAAAD4RIAAAPhEgAAA+ESAAAD4RIAAA\nPhEgAAA+ESAAAD4RIAAAPv0/5iac+Ip1LMMAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our job is to look for the largest $p(i)$ value (blue dot) that is still underneath $q i / N$ (the red line)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The red line $q i / N$ is the acceptable number of false positives $q i$ as a proportion of all the tests $N$. Further to the right on the red line corresponds to a larger acceptable number of false positives. For example, for $i = 1$, the acceptable number of false positives $q * i$ is $0.05 * 1$, but at $i = 50$, the acceptable number of expected false positives $q * i$ is $0.05 * 50 = 2.5$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that, if only the first p value passes threshold, then $p(1) < q \\space 1 \\space / \\space N$. So, if $q = 0.05$, $p(1) < 0.05 / N$. This is the Bonferroni correction for $N$ tests." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The FDR becomes more interesting when there is signal in the noise. In this case there will be p values that are smaller than expected on the null hypothesis. This causes the p value line to start below the diagonal on the ordered plot, because of the high density of low p values." ] }, { "cell_type": "code", "collapsed": false, "input": [ "N_signal = 20\n", "N_noise = N - N_signal\n", "noise_z_values = np.random.normal(size=N_noise)\n", "# Add some signal with very low z scores / p values\n", "signal_z_values = np.random.normal(loc=-2.5, size=N_signal)\n", "mixed_z_values = np.sort(np.concatenate((noise_z_values, signal_z_values)))\n", "mixed_p_values = normal_distribution.cdf(mixed_z_values)\n", "plt.plot(i, mixed_p_values, 'b.', label='$p(i)$')\n", "plt.plot(i, q * i / N, 'r', label='$q i / N$')\n", "plt.xlabel('$i$')\n", "plt.ylabel('$p$')\n", "plt.legend()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAERCAYAAABhKjCtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH81JREFUeJzt3Xtw1NX9//FXaKIJaI0QJCSbFjUxF5BLG4LKrzUSITSF\noFw0OPWCSBgcvogXxLEzX7G1BGxtvaTjoFax2gT4qjM4GlYIGNsRQtQo2IKACHaJXCZiQCCBZPn8\n/liz7CafXEj2s9fnY2Ynu589bE4+ree957zPJcowDEMAALTRJ9AVAAAEJwIEAMAUAQIAYIoAAQAw\nRYAAAJgiQAAATFkeIO6++24NGjRIV199dYdlFixYoLS0NI0YMUKffvqp1VUCAHSD5QFi1qxZstvt\nHb5fUVGhL7/8Unv27NELL7ygefPmWV0lAEA3WB4gfvGLX+jSSy/t8P23335bd955pyRpzJgxamho\n0OHDh62uFgCgCwHPQdTV1SklJcX92maz6cCBAwGsEQBACoIAIUltd/uIiooKUE0AAK2iA12B5ORk\nORwO9+sDBw4oOTm5XTmCBgD0TE+33At4D6KwsFB///vfJUnV1dWKj4/XoEGDTMsahsHDMPTYY48F\nvA7B8uBecC+4F+ce119vSHI9Zsww9Ktf9W4vVst7EDNnztQHH3yg+vp6paSk6PHHH1dzc7Mkae7c\nuSooKFBFRYVSU1PVr18/vfLKK1ZXCQDCRnGxtHu31LevFBPjupadLb3wgut5J3OEumR5gCgvL++y\nTGlpqdXVAICwtHu39MEHrudTpkgzZriCQ3x87z874DkInL/c3NxAVyFocC/O4V6cE0n3om9f18/s\nbGnlSt8EhlZRhmGExIFBUVFRCpGqAoBPeQ4jlZVJDz987vXzz0uLFnXca+hN20mAABASmMnYNbM2\nsjdtJ0NMAEIGXxI7ZkUADfg0VwBAcCJAAABMESAAAKYIEAAAUwQIAIApAgQA+Mm+fftMrx88eFCn\nTp3yc226RoAAAD/46quvVF1dbfrewIED9eSTT/q5Rl0jQACAH6xYsUIzZ840fS86Olq//vWv3Ttb\nBwsCBABYbNu2bbLZbO2uP/LII1q/fr0kafTo0aqsrPR31TrFSmoAsNg777yjm266qd31ZcuWeb0e\nOHCgvvzyS6Wmpvqrap2iBwEgLBQXS7m5UkGB1NAQuM8w89FHHykrK6vLciNGjNAnn3ziu1/cS/Qg\nAIQFz3MRioulNWv8+xnbtm3TJ598ol27dum6667TkSNHdOGFF+qOO+7QqVOnvPZKOnr0qDZs2KA3\n33xTazx+yaWXXqrdu3eff8UtQg8CQFjwPBeh9TQ1f37G4cOHlZ6erv3792vKlCm67bbb9MQTT0iS\nnE6nV9na2lrl5+e3m/YaFxenM2fO9KzyFiBAAAgLZWWu09Q2bOj5oTm9+YwJEyZo/fr1mjx5siTp\n008/VUJCgiTXLCVPN954o1auXKm77rrL6/qxY8fUv3//nlXeAgQIAGEhPt41JNSbE9V6+xmVlZW6\n/vrrJUmvvvqqHnroIUlSYmKiTpw44VW2vLxcv/nNb1RRUeG+dvDgwaBJUEsECADwiWPHjuno0aPa\ntGmTXnzxRY0ZM0ZTp06VJF1//fWqqanxKn/FFVfo3XffVU5OjvvaZ599prFjx/q13p0hSQ0APrBp\n0yYVFhbqzjvvbPfe1KlT9ac//Unjxo1zXysvL/cq09TUpB//+MeKjY21vK7dRQ8CAHrpiy++0J//\n/GcdOXJEx48fb/d+fHy8EhISVF9f3+FnrFq1SnPnzrWymueNM6kBhIRQbwMMw9BLL72kOXPmtHvP\n4XCotrZWU6ZM6fHnd3R/enPfCBAAQgJtQOesCBAMMQEATBEgAACmCBAAAFMECACAKQIEAMAUAQIA\nYIoAAQAwRYAAAJgiQAAATBEgACBAvv/+e+3atSvQ1egQAQIALPbII49o/fr17a6vWbNGF110kfv1\n1q1bdfPNN8tms6mlpUWS66S6oqIiTZo0SZs3b/ZbnSW2+wYAyy1btsz0usPhUHJysvv1mDFjNHHi\nRB07dkxvvvmmbr31Vg0aNEiTJk3StGnTFBcX568qS/JDD8JutysjI0NpaWlavnx5u/fr6+s1ceJE\njRw5UsOGDdPKlSutrhIABNwXX3yhjIwMr2tnz55VTEyMFixYoGeffdZ9/eTJk34PDpLFAcLpdGr+\n/Pmy2+3asWOHysvLtXPnTq8ypaWlGjVqlD777DNVVVXpwQcfdHetAOC8REX55tFDu3fv1qOPPqr1\n69friSeeUFlZmVavXq1bbrmlXdm1a9e6T5xrVVtbq+zsbBUWFurgwYOqra394c/qeZ16w9IAUVNT\no9TUVA0ZMkQxMTEqKirS2rVrvcoMHjzYfcDG8ePHNWDAgHYHfANAtxiGbx49cPLkSU2fPl0PPfSQ\nJkyYoA8//FCNjY3Kz8/Xvn37vMo6nU41Nzfrggsu8Lq+fft2DR8+XH369NG9996r5557Trt27VJ6\nenqPb0lvWNoS19XVKSUlxf3aZrNp69atXmXmzJmjcePGKSkpSd9//73WrFljZZUAwBJvvfWWhg0b\npv79++vMmTM6cuSIZs+eraefflp33XWXV9nKykpNmDCh3WecPXvW/fyee+5RamqqsrKydN9991ld\nfVOWBojudIuWLl2qkSNHqqqqSnv37tX48eO1bds2XXzxxe3KLlmyxP08NzdXubm5PqwtAPRcfX29\nRo0aJcl1PvXYsWMluc6eXr9+vSoqKlRQUCBJqq6u1mOPPeb179v2KOLj4zV9+nS9//77WrRoUbfr\nUVVVpaqqql7+NS6WBojk5GQ5HA73a4fDIZvN5lVm8+bN+u1vfytJuvLKK3X55Zdr165dys7Obvd5\nngECAIJJUVGRSkpK9O677+qpp57SAw88IEm64oor9O6777p7DA0NDbr00ku9/u1HH32kkpIS9e3b\nV3l5ee6ZTQsWLGg3LN+Vtl+eH3/88R7/TZYeOdrS0qL09HRt3LhRSUlJysnJUXl5uTIzM91lHnjg\nAV1yySV67LHHdPjwYf385z/X9u3b1b9/f++KctwgENFCqQ3IycnR+++/r379+rV778UXX9TkyZOV\nmJjo098ZckeORkdHq7S0VPn5+crKytKtt96qzMxMrVixQitWrJAkPfroo/r44481YsQI3XjjjXry\nySfbBQcACAUnT57UkiVL5HA4tGXLFtMy33zzjc+Dg1Us7UH4Uih9ewDge+HQBnz11Vfatm2bbr75\nZp9/thU9CAIEgJBAG9C5kBtiAgCELgIEAMAUAQIAYIoAAQAwRYAAAJhiVzwAISNQu5pGKgIEgJDQ\n3amaxcXS7t1S375SWZkUH29xxcIYQ0wAwsru3dIHH0jr1rmCBXqOAAEgrPTt6/qZnS298EJg6xLq\nWEkNIKw0NLh6Di+8wPCSxFYbAIAOsNUGgIhVXCzl5koFBa7eA3yHAAEgpJGUtg4BAkBIIyltHXIQ\nAEIaSenOkaQGAJgiSQ0gopCY9g8CBICQQ2LaPwgQAEIOiWn/IAcBIOSQmO4+ktQAwk7bXVkffphd\nWnuCAAEg7OTmuvIMkjRjhnTkiPfrNWsCVrWQwiwmAGGnbZ6BvIP/0YMAEJTa5hnIO/QMQ0wAwgKn\nwfkeAQJAyPIMCsePSx9+6LpOnsE3etN2ciY1gIBqXfQmSYmJrp/kGYIDSWoAAeWZfK6udvUcNmxg\neCkYMMQEIKBIPluLHAQAwBQ5CABBzTMRPXCg9PXXzFQKBQQIAJbzTEQnJEj19a7nxcXMVApmJKkB\nWM4zET1y5LnnzFQKbuQgAFjOMxEtkZT2p6BOUtvtdi1cuFBOp1P33HOPFi9e3K5MVVWV7r//fjU3\nNyshIUFVVVXtK0qAAIJOZ7kFdl8NDkEbIJxOp9LT01VZWank5GSNHj1a5eXlyszMdJdpaGjQ2LFj\n9d5778lms6m+vl4JCQntK0qAAIKO546rnrkFdl8NHkG7m2tNTY1SU1M1ZMgQxcTEqKioSGvXrvUq\nU1ZWpmnTpslms0mSaXAAEDw8z4OOiXFdM8stsPtq6LM0QNTV1SklJcX92mazqa6uzqvMnj17dPTo\nUd1www3Kzs7Wa6+9ZmWVAJwnz4DQ0OB9HnS/fudWPv/f/3mvgi4rY1V0qLN0mmtUVFSXZZqbm1Vb\nW6uNGzfq1KlTuvbaa3XNNdcoLS3NyqoBaKOjfILnBnrFxd49g5UrvRt/z2Gk+HiGlUKdpQEiOTlZ\nDofD/drhcLiHklqlpKQoISFBcXFxiouL0y9/+Utt27bNNEAsWbLE/Tw3N1e5ublWVR2IOB2tVTDb\nQI9ZSMGrqqrKdKJPT1iapG5paVF6ero2btyopKQk5eTktEtSf/HFF5o/f77ee+89nT59WmPGjNHq\n1auVlZXlXVGS1ECvdTbr6LbbXMNG2dmuhr+y0vX8jTekRYsICKEqaLfaiI6OVmlpqfLz8+V0OjV7\n9mxlZmZqxYoVkqS5c+cqIyNDEydO1PDhw9WnTx/NmTOnXXAA4BudrWguK+t4rQJDRZGJhXJABCko\nMO8lkEgOX0G7DsKXCBBA77GiOfIQIAB0iHOeI1vQ5iAA+F/bgOCZd2D3VJwPAgQQZtoGBFY0o6cI\nEECYMQsI5BrQE+QggDDDGc/wRJIaAGAqaHdzBQCELgIEEAba7rgK+AIBAggDnltwFxcHujYIFwQI\nIAwwlRVWIEkNhKC2i+FarzFzCW0xiwmIAJ5BwfMQH857RmfYagOIAJ4rpM0O8QF8jRwEECI88wzV\n1Zz3DOsxxASECFZIoyfIQQBhiq260VuWrqR+6qmnlJeXp6FDh+rRRx9Vc3Nzj34RgPPH+gYEUpcB\nIj09XRs3btS///1v5eXl6fe//70/6gVArG9AYHUZIA4dOqSKigqdPHlSeXl5Gj16tD/qBUCuYSWS\n0QiULqe5OhwONTQ06JVXXtG3336rlpYWHTt2THV1dVq8eLE/6giENc88w8CB0tdfe+ccWOOAQOky\nSV1bW6vGxkaNHTtWkrR3715t3rxZL730kj5onZTtBySpEa5yc8+tb0hIkOrrXc9ZAAdfCMgspoMH\nD2rw4ME9+qU9QYBAuCoocCWhs7NdPYbKStdzhpXgC0xzBUKY5/oGibUO8C0CBBDk2q5nePhh1jfA\nPwgQQJDzzDPMmCEdOeL9mlwDrMKRo0CQa7uegfUNCAX0IAA/aLuPEvsqwV8YYgIAmGKICQDgcwQI\nAIApAgRgkeJi1+ylggJXzgEINQQIwCJs1Y1QR4AALMJUVoS6LndzBdB9niumn39eWrSIqawIXQQI\nwIdah5UkV3BghTRCmeVDTHa7XRkZGUpLS9Py5cs7LPfRRx8pOjpab731ltVVAizDsBLCiaUBwul0\nav78+bLb7dqxY4fKy8u1c+dO03KLFy/WxIkTWQyHkMYJcAgnlgaImpoapaamasiQIYqJiVFRUZHW\nrl3brtxzzz2n6dOna+DAgVZWB7Bc6wlwBAeEA0sDRF1dnVJSUtyvbTab6urq2pVZu3at5s2bJ8m1\nLBwIFax1QDizNEndncZ+4cKFWrZsmXu/EIaYEOw8ZyodPy59+OG56ySlEU4sDRDJyclyOBzu1w6H\nQzabzavMJ598oqKiIklSfX291q1bp5iYGBUWFrb7vCVLlrif5+bmKjc315J6A53xnKmUmOj6SVIa\nwaKqqkpVVVU++SxLd3NtaWlRenq6Nm7cqKSkJOXk5Ki8vFyZmZmm5WfNmqXJkydr6tSp7SvKbq4I\nEp5nSL/xBmsdENx603Za2oOIjo5WaWmp8vPz5XQ6NXv2bGVmZmrFihWSpLlz51r56wFLlJV5n+XA\nsBLCFedBAEAY4zwIwGLMVkIkIkAA3cDOrIhE7MUEmPCcylpWxhYaiEwECMCE51TW4uL2iWkgEhAg\nABNtewzMVkIkYhYTYKKhgR4DwkNv2k4CBACEMaa5AgB8jgABADBFgAB+wGI4wBsBAvgBi+EAbwQI\n4AcshgO8MYsJ+AFTWxGOmOYKADDFNFcAgM8RIBDRmLkEdIwAgYjGzCWgYwQIRDRmLgEdI0mNiMbM\nJYQ7ZjEBAEwxiwnoJpLSQPcRIBBRSEoD3UeAQEQhKQ10HzkIRBSS0og0JKkBAKZIUgMAfI4AAQAw\nRYAAAJgiQCDssfYB6BkCBMIeax+AniFAIOyx9gHomehAVwCwQnGxq+fQt6/0/PPSokWsfQDOFwEC\nYal1WElyBYc1awJbHyAUESAQFjx7DGVlDCsBvkCAQFjw7DEUF7uCBFtqAL3DVhsIWZ69huZmqbLS\n1WPYsIGgALQK+q027Ha7MjIylJaWpuXLl7d7/x//+IdGjBih4cOHa+zYsdq+fbs/qoUQ5zl9tV8/\nacYMggPgS5YPMTmdTs2fP1+VlZVKTk7W6NGjVVhYqMzMTHeZK664Qv/85z91ySWXyG63q7i4WNXV\n1VZXDSHIs9cQE+O6lp0trVxJYAB8zfIAUVNTo9TUVA0ZMkSSVFRUpLVr13oFiGuvvdb9fMyYMTpw\n4IDV1UIQ8wwCAwdKX399LvnsmWuYMsXVayDPAFjD8gBRV1enlJQU92ubzaatW7d2WP5vf/ubCgoK\nrK4WgphnEEhIkOrrXc+Li71nJ9FrAKxleYCIiorqdtn3339fL7/8sj788EPT95csWeJ+npubq9zc\n3F7WDsHIMwjEx59LPrdOV2V2EtCxqqoqVVVV+eSzLJ/FVF1drSVLlshut0uSSkpK1KdPHy1evNir\n3Pbt2zV16lTZ7Xalpqa2ryizmCKG56lvEgEB6I2gPlGupaVF6enp2rhxo5KSkpSTk6Py8nKvHMR/\n//tfjRs3Tq+//rquueYa84oSIMJW20VuBALAd3rTdlo+xBQdHa3S0lLl5+fL6XRq9uzZyszM1IoV\nKyRJc+fO1e9+9zt99913mjdvniQpJiZGNTU1VlcNAeQZFI4fl1pHFYuL2RYDCBYslIPfdBQUEhOl\nQ4dY5AZYIegXygGS98K2vXtd17KzpepqFrkBwYgAAb/xnJ3kGRR++lPXsBLBAQguDDHBMm2Tz63X\nmJEE+E9Qz2LyFQJEaOgozzBjBslnIBCCehYTIovnKujERNdPzmQAQhM5CPhUR3kGhpSA0MMQE3zK\ncxU0QQEIPHIQCChWQgPBi3UQCCjP9Q3FxYGuDQBfIUCg1zzzDiSjgfDBEBN6jbwDELzIQQAATJGD\ngF8VF0u5uVJBgav3ACA80YNAt7BCGghNrKSGJTrbnlsiKQ2EOwIE3NquZ+ho24w33pAWLSIpDYQ7\nAgTcPANCcbH39NW2QYFhJSD8ESDgZraewXP6KkEBiCwkqSOc57DS888zdASEG5LU6DHPYaVFi+gl\nADiHdRARjm0yAHSEHkQEaDs76eGHGVYC0DVyEBEgN/fcMNKMGdKRI96vGVYCwhdbbaBTbYeRGFYC\n0B30ICJA291W2X0ViBzs5op2OOUNgMQQE0xwyhuA3mIWUxjx7DXExLiukWcA0FMEiBDW2eZ6U6a4\nZiiRZwCCkGFIzc1SU5PU2Nj1z9bnHV3vqGxTU6+qSYAIYZ1trrdyJYEB6JJhSC0t7Rvbrhrkrhpq\ns+tty/TpI8XGuh5xca5H2+cdvde3rzRgwLlrnj/b/pv09B7fHgJECOtqcz0gZJh9o+7Ot+Suynb2\nGa0/u2qoL7zw3PO2jXFcnNS/f/cbd8+f0cHf/DKLKcSwuR4sdb4NdWcN9/l+C29tqM+nkT3fcmZl\nQqCh7g2muYY5jvuMQOc7Rt2ToZDuNNTdaWQ9G+nzHS7xvB7mDXWgsJtrmOnuyW7MTvIDzzHqnnxL\n7s238PP5Rt22AW8doz6f4RIaarTB/xMCyDMQDBwoff11+15CVye7RYzOGurufEvuSUPdk6GPttfN\nGuquhkNoqBEkLB9istvtWrhwoZxOp+655x4tXry4XZkFCxZo3bp16tu3r1auXKlRo0a1r2gYDDG1\n7RncdNO5nkFCglRf73qemCgdOuQKCBs2nPu3QREUOpr10dMhkM4aerNv1J7fgD0bWl+MSXf0Pg01\nQljQDjE5nU7Nnz9flZWVSk5O1ujRo1VYWKjMzEx3mYqKCn355Zfas2ePtm7dqnnz5qm6utrKavlV\nR/mDtj2D+HipsrLjXoJnrqGqqkq511/vu2Ti+Q6XmI1Rd3cc2nN6Xne/hXfSUFdVVSk3N9fa/xFD\nBPfiHO6Fb1gaIGpqapSamqohQ4ZIkoqKirR27VqvAPH222/rzjvvlCSNGTNGDQ0NOnz4sAYNGmRl\n1XrN3fDHGUpMaNGhfY265MImDY5vVL2jUT++oEkl/9uovlub1Hd7o2LVpMsvaVS6mpT+k0b9z7BG\n6XST7Psa9aucJkU1NmrzkCaNTWnUBfc1aU1jozSlTaP/QwNddeyYcs+ebd9Qd+dbc+vzfv3Mx6h7\n2FAHCg3BOdyLc7gXvmHpf+l1dXVKSUlxv7bZbNq6dWuXZQ4cOHD+AaKrWR/dHP6osjep8Wij4vo0\n6eLoRjlPNClOjer7I1eZWKNRg+Ob9McjjfpRi+u9s+qjJsWqUXE63SdOjWcvVKPidHBmnO5WrOoU\np9j4WOX8Mk4122N1XV6cLjzuanRv/k0/Kdb1jfqG/9fJXGvPRvovf5F+97ugaqgBhB9LW5ioqKhu\nlWs7PtbZv3tsxg7dsa7I3VCf+rZRP2pu1AVnm2RE9dGZPrE686M4tcTE6dTZWDVHx+mn6bHaWxer\nhqY4tVwQpz5xsTraFKezF8Qp+qJY1Z+I1dkL++mWOwbos+ZYbT8cp0bF6YKLY3Xke1fDf+ElcTp0\n0vU874ZYfdcYp3c2xWnYz2N18aXR7uEhz6Gi1vzBEz/kD/rFSzf0+G56uOACggMA6xkW2rJli5Gf\nn+9+vXTpUmPZsmVeZebOnWuUl5e7X6enpxuHDh1q91mSePDgwYNHDx49ZenX0OzsbO3Zs0f79+9X\nUlKSVq9erfLycq8yhYWFKi0tVVFRkaqrqxUfH286vGSE+AwmAAg1lgaI6OholZaWKj8/X06nU7Nn\nz1ZmZqZWrFghSZo7d64KCgpUUVGh1NRU9evXT6+88oqVVQIAdFPIbLUBAPCvoD9Rzm63KyMjQ2lp\naVq+fHmgq+NXDodDN9xwg4YOHaphw4bp2WeflSQdPXpU48eP11VXXaUJEyaooaEhwDX1H6fTqVGj\nRmny5MmSIvdeNDQ0aPr06crMzFRWVpa2bt0asfeipKREQ4cO1dVXX63bbrtNp0+fjph7cffdd2vQ\noEG6+uqr3dc6+9tLSkqUlpamjIwMrV+/vsvPD+oA0brQzm63a8eOHSovL9fOnTsDXS2/iYmJ0V/+\n8hf95z//UXV1tf76179q586dWrZsmcaPH6/du3crLy9Py5YtC3RV/eaZZ55RVlaWe6ZbpN6L++67\nTwUFBdq5c6e2b9+ujIyMiLwX+/fv14svvqja2lp9/vnncjqdWrVqVcTci1mzZslut3td6+hv37Fj\nh1avXq0dO3bIbrfr3nvv1dmzZzv/BT1Ob/vB5s2bvWZBlZSUGCUlJQGsUWBNmTLF2LBhg9dMr4MH\nDxrp6ekBrpl/OBwOIy8vz9i0aZMxadIkwzCMiLwXDQ0NxuWXX97ueiTei2+//da46qqrjKNHjxrN\nzc3GpEmTjPXr10fUvdi3b58xbNgw9+uO/va2s0jz8/ONLVu2dPrZQd2DMFtEV1dXF8AaBc7+/fv1\n6aefasyYMV4rzQcNGqTDhw8HuHb+cf/99+uPf/yj+vQ593/bSLwX+/bt08CBAzVr1iz97Gc/05w5\nc3Ty5MmIvBf9+/fXgw8+qJ/85CdKSkpSfHy8xo8fH5H3olVHf/s333wjm83mLted9jSoA0R3F9qF\nuxMnTmjatGl65plndPHFF3u9FxUVFRH36Z133tFll12mUaNGdTjlOVLuRUtLi2pra3XvvfeqtrZW\n/fr1azeEEin3Yu/evXr66ae1f/9+ffPNNzpx4oRef/11rzKRci/MdPW3d3VfgjpAJCcny+FwuF87\nHA6vCBgJmpubNW3aNN1+++266aabJLm+FRw6dEiSdPDgQV122WWBrKJfbN68WW+//bYuv/xyzZw5\nU5s2bdLtt98ekffCZrPJZrNp9OjRkqTp06ertrZWiYmJEXcvPv74Y1133XUaMGCAoqOjNXXqVG3Z\nsiUi70Wrjv6baNueHjhwQMnJyZ1+VlAHCM+FdmfOnNHq1atVWFgY6Gr5jWEYmj17trKysrRw4UL3\n9cLCQr366quSpFdffdUdOMLZ0qVL5XA4tG/fPq1atUrjxo3Ta6+9FpH3IjExUSkpKdq9e7ckqbKy\nUkOHDtXkyZMj7l5kZGSourpajY2NMgxDlZWVysrKish70aqj/yYKCwu1atUqnTlzRvv27dOePXuU\nk5PT+Yf5OmHiaxUVFcZVV11lXHnllcbSpUsDXR2/+te//mVERUUZI0aMMEaOHGmMHDnSWLdunfHt\nt98aeXl5RlpamjF+/Hjju+++C3RV/aqqqsqYPHmyYRhGxN6Lzz77zMjOzjaGDx9u3HzzzUZDQ0PE\n3ovly5cbWVlZxrBhw4w77rjDOHPmTMTci6KiImPw4MFGTEyMYbPZjJdffrnTv/0Pf/iDceWVVxrp\n6emG3W7v8vNZKAcAMBXUQ0wAgMAhQAAATBEgAACmCBAAAFMECACAKQIEAMAUAQIAYIoAAQAwRYAA\neqm5uVkzZ84MdDUAn2MlNQDAFD0IAICp6EBXAAhlX331ld555x0lJSVp+vTpga4O4FP0IIBeOHTo\nkAYMGKDTp08HuiqAz5GDAHrplltu0csvv6yLLroo0FUBfIoeBNALx48fV1RUlD7//PNAVwXwOQIE\n0AtOp1ODBg3SmTNnAl0VwOcYYgIAmKIHAQAwRYAAAJgiQAAATBEgAACmCBAAAFMECACAKQIEAMAU\nAQIAYOr/AyYSIM1eXUY7AAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interesting part is the beginning of the graph, where the blue p values stay below the red line:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "first_i = i[:30]\n", "plt.plot(first_i, mixed_p_values[:30], 'b.', label='$p(i)$')\n", "plt.plot(first_i, q * first_i / N, 'r', label='$q i / N$')\n", "plt.xlabel('$i$')\n", "plt.ylabel('$p$')\n", "plt.legend()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAERCAYAAACKHYuuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wk/XhB/B3arMhVAnU0mLSXU8SaFqgrUvJkHMEUSpV\n4g/qrIrUGaTHjqs43an7Z2XbIdWxWej+aHGifDdKuWOjPUhz0GmUA0s9W9RdQcpovKehrRMM5acp\n4fn+0ZE1bdqnSfPkV9+vu5z58XnSz+Ojn3c+n+f5fB6FKIoiiIiIRpEQ6QoQEVH0Y1gQEZEkhgUR\nEUliWBARkSSGBRERSWJYEBGRJNnDwmazITMzEzqdDhUVFcM+P3nyJBYuXIhJkyZhy5Ytwz73eDzI\ny8vDihUr5K4qERGNIFHOL/d4PFi/fj2ampqgVquRn58Ps9kMvV7vLZOcnIxt27Zh3759fr+jsrIS\nWVlZuHjxopxVJSKiUcjas2hpaYFWq0VGRgaUSiWKi4tRX1/vUyYlJQUGgwFKpXLY9l1dXbBarViz\nZg04d5CIKHJkDQun04n09HTva41GA6fTOebtX3rpJbz11ltISOCpFSKiSJK1FVYoFEFvu3//fsyY\nMQN5eXnsVRARRZis5yzUajUEQfC+FgQBGo1mTNsePXoUDQ0NsFqtuHbtGvr6+rB69Wrs3LnTp9x4\nAomIaCIL5Ie4rD0Lg8GAjo4OOBwOuN1u1NXVwWw2+y07tNKbNm2CIAjo7OzE7t27cd999w0LisHb\nxuvjN7/5TcTrwP3j/k3E/YvnfRPFwEdrZO1ZJCYmoqqqCgUFBfB4PLBYLNDr9aiurgYAlJaWoqen\nB/n5+ejr60NCQgIqKyvR3t6OpKQkn+9iD4KIKHJkDQsAWL58OZYvX+7zXmlpqfd5Wlqaz1CVP4sX\nL8bixYtlqR8REUnjZUZRzmQyRboKsuL+xbZ43r943rdgKMRgBq+iiEKhCGr8jYhoIgu07ZR9GIqI\nSG48pzm6UPygZlgQUVzgCIN/oQpSnrMgIiJJDAsiIpLEsCAiIkkMCyIiksSwICIiSQwLIqII6+zs\n9Pt+d3c3rly5Euba+MewICKKoDNnzqC5udnvZykpKXjzzTfDXCP/GBZERBFUXV2Np556yu9niYmJ\neOihh0ZccTucGBZERBHy+eef+73Hz2uvvYaDBw8CAPLz89HU1BTuqg3DGdxERBGyf/9+PProo8Pe\n37x5s8/rlJQUnD59GlqtNlxVG4Y9CyKKa2vXAiYTUFgIuFyR+w5/Pv30U2RlZUmWy8nJwWeffRa6\nPxwE9iyIKK6dOgV89NHA87VrgT17wvsdn3/+OT777DN89dVXuOeee/DNN9/ghz/8IVavXo0rV674\nrN10/vx5HDp0CHv37sWeQX9k2rRpOHXqVOAVDyH2LIgork2ePPBPgwGoqQn/d/T29mLOnDlwOBx4\n5JFH8PTTT+P3v/89AMDj8fiUbW1tRUFBwbBLaW+99Va43e7gKh8iDAsiimu7dgFPPAEcOgSoVOH/\njmXLluHgwYNYsWIFAKCtrQ133HEHgIGrnQa7//778d577+G5557zef/ChQuYPn16cJUPEYYFEcU1\nlWpg2CjYoAjFdzQ1NXlvDf3+++/jlVdeATBwW+lLly75lK2trcWqVatgtVq973V3d0f05DbAsCAi\nktWFCxdw/vx5fPDBB9i+fTuMRiMef/xxAMDixYvR0tLiU/6uu+7CgQMHsGDBAu97x48fx6JFi8Ja\n76F4gpuISEYffPABzGYzSkpKhn32+OOP4w9/+APuu+8+73u1tbU+Za5du4bbb78dkyZNkr2uo2HP\ngohIJidPnsQf//hHfPPNN+jr6xv2uUqlwh133IFvv/12xO/YvXs3SktL5azmmCjEGL8XYaA3HSei\n+BPL7YAoinjnnXfwwgsvDPtMEAS0trbikUceCfr7R/p3E+i/s7D0LGw2GzIzM6HT6VBRUTHs85Mn\nT2LhwoWYNGkStmzZ4n1fEAQsWbIE2dnZmDt3LrZu3RqO6hIRhY1CofAbFACQnp4+rqAIJdnPWXg8\nHqxfvx5NTU1Qq9XIz8+H2WyGXq/3lklOTsa2bduwb98+n22VSiX+9Kc/ITc3F5cuXcKPf/xjPPDA\nAz7bEhFRYNauDXwb2XsWLS0t0Gq1yMjIgFKpRHFxMerr633KpKSkwGAwQKlU+ryflpaG3NxcAEBS\nUhL0ej3Onj0rd5WJiOJaMJPBZQ8Lp9OJ9PR072uNRgOn0xnw9zgcDrS1tcFoNIayekREE87NGemB\nkH0YavC6J8G6dOkSioqKUFlZiaSkpGGfl5eXe5+bTCaYTKZx/00ionhit9tht9sBAPPnA42NgW0v\ne1io1WoIguB9LQiC3/XbR9Lf34+VK1di1apVfpfyBXzDgoiIhhv6Q7qiYmNA28s+DGUwGNDR0QGH\nwwG32426ujqYzWa/ZYdexiWKIiwWC7KysrBhwwa5q0pERCMIyzyLxsZGbNiwAR6PBxaLBa+//jqq\nq6sBAKWlpejp6UF+fj76+vqQkJCA2267De3t7Th+/Dh++tOfYv78+d7hrDfeeAMPPvjg/3Yghq+v\nJqLQYDswslDNs+CkPCKKeWwHRhZTk/KIiCi2MSyIiKLMxYsX8dVXX0W6Gj4YFkREEfLaa6/h4MGD\nw97fs2ePzzSBY8eO4bHHHoNGo8H169cBDNyBr7i4GA8//DCOHj0qe125RDkRUYRs3rzZ7/uCIECt\nVntfG41GPPjgg7hw4QL27t2LJ598EqmpqXj44YexcuVK3HrrrbLXlT0LIqIocvLkSWRmZvq8d+PG\nDSiVSpSVlfksqHr58uWwBAXAngURTQQhWEkCABDkFVenTp3Ce++9B5PJhJaWFtx111245ZZbsHfv\nXuzZs8enbH19PV566SWf91pbW2EwGDB37lz88pe/RGtrK+6+++6QrJAxVuxZEFH8E8XQPIJw+fJl\nFBUV4ZVXXsGyZctw5MgRXL16FQUFBejs7PQp6/F40N/fjx/84Ac+73/xxReYP38+EhIS8Itf/ALb\ntm3DV199hTlz5gT9ryRQ7FkQEcno73//O+bOnYvp06fD7Xbjm2++gcViwdtvv43nnnvOp2xTUxOW\nLVs27Dtu3Ljhfb5mzRpotVpkZWXhxRdflLv6XuxZEBHJ6Ntvv0VeXh6AgftxL1q0CMDAvbZXrVoF\nq9XqLdvc3IwFCxb4bD+0p6FSqVBUVIQPP/xwWA9ETgwLIiIZFRcXQxAEHDhwAG+++aa353DXXXfh\nwIED3nBwuVyYNm2az7affvopnnzySRw8eNDn1g5lZWW49957w7cT4HIfRBQHYqUdWLBgAT788ENM\nmTJl2Gfbt2/HihUrkJaWFtK/yeU+iIhixOXLl1FeXg5BEPDJJ5/4LXP27NmQB0UosWdBRDEv1tuB\nM2fO4PPPP8djjz0W8u/mqrP/Fev/kRDR+LEdGBmHoYiIKGwYFkREJIlhQUREkhgWREQkiWFBRESS\nuDYUEcWFcK7AOhExLIgo5vGyWflxGIqIiCQxLIiISJLsYWGz2ZCZmQmdToeKiophn588eRILFy7E\npEmTsGXLloC2JSKi8JB1uQ+Px4M5c+agqakJarUa+fn5qK2thV6v95b5z3/+g6+//hr79u3DtGnT\n8PLLL495W4DT/ImIghFVy320tLRAq9UiIyMDSqUSxcXFqK+v9ymTkpICg8EApVIZ8LZERBQesoaF\n0+lEenq697VGo/G5gYdc2xIRUWjJeunseK57DmTb8vJy73OTyQSTyRT03yUiikd2ux12uz3o7WUN\nC7VaDUEQvK8FQYBGown5toPDgoiIhhv6Q3rjxo0BbS/rMJTBYEBHRwccDgfcbjfq6upgNpv9lh16\noiWQbYmISF6y9iwSExNRVVWFgoICeDweWCwW6PV6VFdXAwBKS0vR09OD/Px89PX1ISEhAZWVlWhv\nb0dSUpLfbYmIKPx4pzwiogkoqi6dJSKi+MCwICIiSQwLIiKSxLAgIiJJDAsiIpLEsCAiGoe1awGT\nCSgsBFyuSNdGPgwLIqJxOHUK+OgjoLFxIDjiFcOCiGgcJk8e+KfBANTURLYucuKkPCKicXC5BnoU\nNTWAShXp2oxdoG0nw4KIaALiDG4iIgo5hgUREUliWBARkSSGBRERSWJYEBGRJIYFERFJYlgQEZEk\nhgUREUliWBARkSSGBRERSWJYEBGRJIYFERFJYlgQEZEk2cPCZrMhMzMTOp0OFRUVfsuUlZVBp9Mh\nJycHbW1t3vffeOMNZGdnY968eXj66afx/fffy11dIiLyQ9aw8Hg8WL9+PWw2G9rb21FbW4sTJ074\nlLFarTh9+jQ6OjpQU1ODdevWAQAcDge2b9+O1tZWfPnll/B4PNi9e7ec1SUiohHIGhYtLS3QarXI\nyMiAUqlEcXEx6uvrfco0NDSgpKQEAGA0GuFyudDb24vbb78dSqUSV65cwfXr13HlyhWo1Wo5q0tE\nBGDi3Fc7ELKGhdPpRHp6uve1RqOB0+kcU5np06fj5Zdfxo9+9CPceeedUKlUuP/+++WsLhERgIlz\nX+1AJMr55QqFYkzl/N2t6d///jfefvttOBwOTJ06FU888QT+9re/4ZlnnhlWtry83PvcZDLBZDIF\nW2Uiori8r7bdbofdbg96e1nDQq1WQxAE72tBEKDRaEYt09XVBbVaDbvdjnvuuQfJyckAgMcffxxH\njx6VDAsiovHatSs276s9mqE/pDdu3BjQ9rIOQxkMBnR0dMDhcMDtdqOurg5ms9mnjNlsxs6dOwEA\nzc3NUKlUSE1NxZw5c9Dc3IyrV69CFEU0NTUhKytLzuoSURwL5DyESgXs2RM/QREKsvYsEhMTUVVV\nhYKCAng8HlgsFuj1elRXVwMASktLUVhYCKvVCq1WiylTpmDHjh0AgNzcXKxevRoGgwEJCQm4++67\nsZaDh0QUpJvnIYCB4NizJ7L1iTUK0d8JgxiiUCj8nvMgIhqssHDghLXBABw6xF5DoG0nw4KIJgSX\nK/7OQ4wHw4KIiCQF2nZybSgiIpLEsCAiIkkMCyIiksSwICIiSQwLIiKSxLAgIiJJDAsiIpLEsCAi\nIkmSYbFlyxYsXboU2dnZ+PWvf43+/v5w1IuIiKKIZFjMmTMH//znP/Gvf/0LS5cuxe9+97tw1IuI\niKKIZFj09PTAarXi8uXLWLp0KfLz88NRLyIiiiKSS5QLggCXy4UdO3bg3LlzuH79Oi5cuACn04lX\nX301HHUkIqIIk1xIsLW1FVevXsWiRYsADNzu9OjRo3jnnXfw0c3F4SOICwkSEQUubKvOdnd3Y+bM\nmcFsGlIMCyKiwHGJciIiksQlyomIKOQYFkREJIlhQUREkhgWREQkiWFBRESSGBZERCRJ9rCw2WzI\nzMyETqdDRUWF3zJlZWXQ6XTIyclBW1ub932Xy4WioiLo9XpkZWWhublZ7uoSEZEfsoaFx+PB+vXr\nYbPZ0N7ejtraWpw4ccKnjNVqxenTp9HR0YGamhqsW7fO+9mLL76IwsJCnDhxAl988QX0er2c1SUi\nohHIGhYtLS3QarXIyMiAUqlEcXEx6uvrfco0NDSgpKQEAGA0GuFyudDb24sLFy7g8OHDeP755wEA\niYmJmDp1qpzVJaIYs3YtYDIBhYWAyxXp2sQ3WcPC6XQiPT3d+1qj0cDpdEqW6erqQmdnJ1JSUvDz\nn/8cd999N1544QVcuXJFzuoSUYw5dQr46COgsXEgOEg+kqvOjodCoRhTuaFTzhUKBa5fv47W1lZU\nVVUhPz8fGzZswObNm/Hb3/522Pbl5eXe5yaTCSaTaTzVJqIYMXnywD8NBqCmJrJ1iXZ2ux12uz3o\n7WUNC7VaDUEQvK8FQYBGoxm1TFdXF9RqNURRhEaj8d4/o6ioCJs3b/b7dwaHBRFNHLt2DfQoamoA\nlSrStYluQ39Ib9y4MaDtZR2GMhgM6OjogMPhgNvtRl1dHcxms08Zs9mMnTt3AgCam5uhUqmQmpqK\ntLQ0pKen49SpUwCApqYmZGdny1ldIooxKhWwZw+DIhxk7VkkJiaiqqoKBQUF8Hg8sFgs0Ov1qK6u\nBgCUlpaisLAQVqsVWq0WU6ZMwY4dO7zbb9u2Dc888wzcbjdmzZrl8xkREYUPlygnIpqAuEQ5ERGF\nHMOCiIgkMSyIiEgSw4KIiCQxLIgoqnAJj+jEsCCiqMIlPKITw4KIogqX8IhOnGdBRFHF5eISHuEQ\naNvJsCAimoA4KY+IiEKOYUFEsuMVTrGPYUFEsuMVTrGPYUFEsuMVTrGPJ7iJSHa8win68GooIiKS\nxKuhiIgo5BgWREQkiWFBRESSGBZERCSJYUFERJIYFkREJIlhQUREkhgWREQkSfawsNlsyMzMhE6n\nQ0VFhd8yZWVl0Ol0yMnJQVtbm89nHo8HeXl5WLFihdxVJSKiEcgaFh6PB+vXr4fNZkN7eztqa2tx\n4sQJnzJWqxWnT59GR0cHampqsG7dOp/PKysrkZWVBYVCIWdViYhoFLKGRUtLC7RaLTIyMqBUKlFc\nXIz6+nqfMg0NDSgpKQEAGI1GuFwu9Pb2AgC6urpgtVqxZs0aLulBRBRBsoaF0+lEenq697VGo4HT\n6RxzmZdeeglvvfUWEhJ4aoWIKJIS5fzysQ4dDe01iKKI/fv3Y8aMGcjLy4Pdbh91+/Lycu9zk8kE\nk8kUYE2JiOKb3W6XbEtHI2tYqNVqCILgfS0IAjQazahlurq6oFarsXfvXjQ0NMBqteLatWvo6+vD\n6tWrsXPnzmF/Z3BYEBHRcEN/SG/cuDGg7WUd3zEYDOjo6IDD4YDb7UZdXR3MZrNPGbPZ7A2A5uZm\nqFQqpKWlYdOmTRAEAZ2dndi9ezfuu+8+v0FBRETyk7VnkZiYiKqqKhQUFMDj8cBisUCv16O6uhoA\nUFpaisLCQlitVmi1WkyZMgU7duzw+128GoqIKHJ48yMiCsratQP31p48Gdi1i3fAizW8+RERhcWp\nU8BHHwGNjQPBQfGNYUFEQZk8eeCfBsPAvbUpvnEYioiC4nIN9ChqajgEFYsCbTsZFkREExDPWRAR\nUcgxLIiISBLDgoiIJDEsiIhIEsOCiIgkMSyICMDAZbAmE1BYOHBZLNFgDAsiAsAZ2TQ6hgURAeCM\nbBodJ+UREQDOyJ5oOIObiIgkcQY3EXnxpDWFCsOCKI7xpDWFCsOCKI7xpDWFCs9ZEMUxnrSmkfAE\nNxERSeIJbiIiCjmGBRERSWJYEBGRJIYFERFJCktY2Gw2ZGZmQqfToaKiwm+ZsrIy6HQ65OTkoK2t\nDQAgCAKWLFmC7OxszJ07F1u3bg1HdYmiGifaUUSIMrt+/bo4a9YssbOzU3S73WJOTo7Y3t7uU+bA\ngQPi8uXLRVEUxebmZtFoNIqiKIrd3d1iW1ubKIqiePHiRXH27NnDtg3DLhDJ7oUXRHHxYlFcvlwU\nv/tu9LKLF4siMPB44olw1I7iUaBtp+w9i5aWFmi1WmRkZECpVKK4uBj19fU+ZRoaGlBSUgIAMBqN\ncLlc6O3tRVpaGnJzcwEASUlJ0Ov1OHv2rNxVJgq7QGZac6IdRYLsYeF0OpGenu59rdFo4HQ6Jct0\ndXX5lHE4HGhra4PRaJS3wkQREEgA7NoFPPEEcOgQJ9pR+CTK/QcUCsWYyolDJocM3u7SpUsoKipC\nZWUlkpKShm1bXl7ufW4ymWAymYKqK1Gk7No19pnWKhWwZ0946kXxw263w263B7297GGhVqshCIL3\ntSAI0Gg0o5bp6uqCWq0GAPT392PlypVYtWoVHn30Ub9/Y3BYEMUiBgDJbegP6Y0bNwa0vezDUAaD\nAR0dHXA4HHC73airq4PZbPYpYzabsXPnTgBAc3MzVCoVUlNTIYoiLBYLsrKysGHDBrmrSkREI5C9\nZ5GYmIiqqioUFBTA4/HAYrFAr9ejuroaAFBaWorCwkJYrVZotVpMmTIFO3bsAAAcOXIEf/3rXzF/\n/nzk5eUBAN544w08+OCDclebaNzWrh04cT158sAwE88vUCzjQoJEMjGZBq5wAgZOSHOYiaIJFxIk\nihK8xJXiCXsWRDLhvSQomgXadsp+zoJoouIVThQxHs/Ar5Xz5wce58797/nNR4AYFkRE0Wpwoz+0\nwR/8euhnfX3A7bcDycnA9On/e9x8rdUGXBUOQxERyS3YRv/CBd9Gf3DjP9rzqVOBW24ZtUq8rSoR\nkVxGavSlAqCvb6ABnz4dmDZtoEEf6Vf/4OcqlWSjHyyGBRGRlJuN/miNvFSjP7hhlwoAGRv9YDEs\niGji8HgGhmpuNuoj/dqXavT9NfI3AyDKG/1gMSyIZMRZ2TIZ2uiPdWy/rw+47Tb/Y/r+hnZuBkAc\nNfrBYlgQyYizsiUMvWQzkEb/9tulT94Ofc1GP2icZ0EUgEB7ChNmVrZUoz9SAEg1+lotsGABG/0Y\nxJ4FxZ1AAiDQnkLMzcqWmpw1UgAMbvRH+mXv76QuG/2YwZ4FTXg3b1EKDDTsowVAoD2FiM3KHk+j\nf3NM31+jr9P5DwA2+jQEw4LiTqC3KA1rT2G0Rn+0ABj6S3/or/qhjX4Ak7OIxoLDUBQTAhlaCstQ\n0Xgb/aG/8kMwI5coELwaiuKSbFchjXftnbE09Gz0KQrxnAXFJcmhpVCtvTO0gddq/Y/1c0yfJhj2\nLCi6jNDoX+k6hwP/dx7mRefxw8sBLMMQwbV3iKIZh6EoOoxhTP9Y43kozp/D1BvnoZ1+Hrd8N4ZG\nP4bW3iGKZhyGotDytwzDWE7ojjam/9+rd5qapuHId8k4j+lYsDAZW//KRp8oWrFnMVFE4do7hYVA\nY+PAeYhDh2JkkhtRnOAwVLwL1do7gxv6EK69E3WXuBKRX1EXFjabDRs2bIDH48GaNWvw6quvDitT\nVlaGxsZGTJ48Ge+99x7y8vLGvG3MhkUo1t4J5Fr9MA3vcKE9otgQVecsPB4P1q9fj6amJqjVauTn\n58NsNkOv13vLWK1WnD59Gh0dHTh27BjWrVuH5ubmMW0bFeRae+e/M3LtXV0wLVkS0bV3AuktBLp8\nht1uh8lkCkk9oxH3L3bF874FQ9awaGlpgVarRUZGBgCguLgY9fX1Pg1+Q0MDSkpKAABGoxEulws9\nPT3o7OyU3DakQr32zmjLMATQ6NvLy2EyGuXZ5zEKZK2lQJfPiPf/Ibl/sSue9y0YsoaF0+lEenq6\n97VGo8GxY8ckyzidTpw9e1Zy24D09wO///3YZuRGydo7a9cCBw8CLS3Sv+gD+fUv57LcEVtoj4hk\nJWtYKBSKMZULxzmHtb9IxD2HgO8na7HqxWRMSR99TN/boH4XuYb61Cng668HHlK/6AP59R9IWSAC\ni+0RUfQRZfTJJ5+IBQUF3tebNm0SN2/e7FOmtLRUrK2t9b6eM2eO2NPTM6ZtRVEUAfDBBx988BHE\nIxCy9iwMBgM6OjrgcDhw5513oq6uDrW1tT5lzGYzqqqqUFxcjObmZqhUKqSmpiI5OVlyWwzsrZy7\nQEREkHkYKjExEVVVVSgoKIDH44HFYoFer0d1dTUAoLS0FIWFhbBardBqtZgyZQp27Ngx6rZERBR+\nMT8pj4iI5JcQ6QqMh81mQ2ZmJnQ6HSoqKiJdnZDLyMjA/PnzkZeXhwULFkS6OuPy/PPPIzU1FfPm\nzfO+d/78eTzwwAOYPXs2li1bBpfLFcEajo+//SsvL4dGo0FeXh7y8vJgs9kiWMPxEQQBS5YsQXZ2\nNubOnYutW7cCiJ9jONL+xcMxvHbtGoxGI3Jzc5GVlYXXX38dQBDHLqAzHFHk+vXr4qxZs8TOzk7R\n7XaLOTk5Ynt7e6SrFVIZGRniuXPnIl2NkPj444/F1tZWce7cud73fvWrX4kVFRWiKIri5s2bxVdf\nfTVS1Rs3f/tXXl4ubtmyJYK1Cp3u7m6xra1NFEVRvHjxojh79myxvb09bo7hSPsXL8fw8uXLoiiK\nYn9/v2g0GsXDhw8HfOxitmcxeMKfUqn0TtqLN2KcjBLee++9mDZtms97gydklpSUYN++fZGoWkj4\n2z8gfo5fWloacnNzAQBJSUnQ6/VwOp1xcwxH2j8gPo7h5P9OlnK73fB4PJg2bVrAxy5mw2KkyXzx\nRKFQ4P7774fBYMD27dsjXZ2Q6+3tRWpqKgAgNTUVvb29Ea5R6G3btg05OTmwWCwxO0QzlMPhQFtb\nG4xGY1wew5v795Of/ARAfBzDGzduIDc3F6mpqd7htkCPXcyGxVgn/MWyI0eOoK2tDY2Njfjzn/+M\nw4cPR7pKslEoFHF3TNetW4fOzk4cP34cM2fOxMsvvxzpKo3bpUuXsHLlSlRWVuK2227z+SwejuGl\nS5dQVFSEyspKJCUlxc0xTEhIwPHjx9HV1YWPP/4YH374oc/nYzl2MRsWarUagiB4XwuCAI1GE8Ea\nhd7MmTMBACkpKXjsscfQ0tIS4RqFVmpqKnp6egAA3d3dmDFjRoRrFFozZszw/k+4Zs2amD9+/f39\nWLlyJZ599lk8+uijAOLrGN7cv1WrVnn3L96O4dSpU/HQQw/hs88+C/jYxWxYDJ7w53a7UVdXB7PZ\nHOlqhcyVK1dw8eJFAMDly5dx8OBBnytt4oHZbMb7778PAHj//fe9/4PGi+7ubu/zf/zjHzF9/ERR\nhMViQVZWFjZs2OB9P16O4Uj7Fw/H8Ntvv/UOn129ehWHDh1CXl5e4MdOzjPwcrNareLs2bPFWbNm\niZs2bYp0dULqzJkzYk5OjpiTkyNmZ2fH/P4VFxeLM2fOFJVKpajRaMR3331XPHfunLh06VJRp9OJ\nDzzwgPjdd99FuppBG7p/f/nLX8Rnn31WnDdvnjh//nzxkUceEXt6eiJdzaAdPnxYVCgUYk5Ojpib\nmyvm5uaKjY2NcXMM/e2f1WqNi2P4xRdfiHl5eWJOTo44b9488c033xRFUQz42HFSHhERSYrZYSgi\nIgofhgUREUliWBARkSSGBRERSWJYEBGRJIYFERFJYlgQEZEkhgUREUliWBDJpL+/H0899VSkq0EU\nEpzBTUSgzW14AAAAq0lEQVREktizICIiSYmRrgBRPDpz5gz279+PO++8E0VFRZGuDtG4sWdBJIOe\nnh4kJyfj+++/j3RViEKC5yyIZPKzn/0M7777LpKSkiJdFaJxY8+CSAZ9fX1QKBT48ssvI10VopBg\nWBDJwOPxIDU1FW63O9JVIQoJDkMREZEk9iyIiEgSw4KIiCQxLIiISBLDgoiIJDEsiIhIEsOCiIgk\nMSyIiEgSw4KIiCT9P82tOnSwlrNHAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are looking for the largest $p(i) < qi/N$, which corresponds to the last blue point below the red line." ] }, { "cell_type": "code", "collapsed": false, "input": [ "below = mixed_p_values < (q * i / N) # True where p(i)