{
 "metadata": {
  "name": "010-Multiple_comparison"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Multiple comparisons"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# first set up inline plotting\n",
      "%pylab inline "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].\n",
        "For more information, type 'help(pylab)'.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The problem occurs when you permform several statistical tests: by chance, some will be significant even if the data has been generated under the null hypothesis.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "FWER : Bonferroni"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Imaging you are doing $n$ t tests each one with a risk of type one error of $\\alpha$, and we denote by $t_\\alpha$ the t-value such that for a variable t drawn from a Student distribution, we have $P(t < t_\\alpha) = \\alpha$,   the probability of getting at least one significant test is these n tests is : "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "\\alpha_{F} = P(\\textrm{ one or more } t_i > t_\\alpha) = 1 - \\prod_{i} P(t_i < t_\\alpha) = 1 - (1 - \\alpha)^n \n",
      "$$ \n",
      "\n",
      "$$\n",
      "\\alpha = 1 -(1 - \\alpha_{F})^{\\frac{1}{n}} = 1 - \\exp\\left(\\frac{1}{n} \\log(1 - \\alpha_{F})\\right)\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# code demo: a t-distribution\n",
      "from scipy.stats import t as tdist\n",
      "df = 14\n",
      "n = 100\n",
      "alph = .05\n",
      "\n",
      "# Often easier to compute the static distribution of t with df : t14 = tdist(df)\n",
      "# rvs : draw one or several random variables\n",
      "ts = tdist(df).rvs(size=n)\n",
      "\n",
      "# The threshold as a t value: \n",
      "t_alph = tdist.isf(alph,df)\n",
      "print t_alph\n",
      "print \"Nubmer of test significant at the %f level: %d \" %(alph, sum(ts >= t_alph))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.76131013577\n",
        "Nubmer of test significant at the 0.050000 level: 4 \n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To correct for the number of tests with the Bonferroni correction :"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# if we have p values : \n",
      "p = tdist(df).sf(ts)\n",
      "p_corr = (1 - (1 - p)**n)\n",
      "print \"Number of test significant with uncorrected threshold: %d, corrected: %d \" \\\n",
      "            % (sum(p <= alph), sum(p_corr <= alph))\n",
      "\n",
      "# This can be put in a little function:\n",
      "def p_corrected(ps):\n",
      "    return( 1 - (1 - p)**n )\n",
      "\n",
      "# alternatively : compare orginal p with : 1.0 - exp(1./n * log(1 - .05))\n",
      "\n",
      "def alph_corr(alph,n):\n",
      "    return 1.0 - exp(1./n * log(1 - alph))\n",
      "\n",
      "print \"Uncorrected threshold: %f, corrected threshold: %f\" % (alph, alph_corr(alph,n))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Number of test significant with uncorrected threshold: 4, corrected: 0 \n",
        "Uncorrected threshold: 0.050000, corrected threshold: 0.000513\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "On average we would have to run about 20 times the previous code to get a (spurious) significant effect.\n",
      "\n",
      "Note that when the tests are NOT INDEPENDENT, then this procedure is conservative (on average you will get less false positive than the specified risk or error), and therefore not the most sensitive. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "False Discovery Rate "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The false detection rate (introduced by Benjamini and Hochberg (1995), see [http://en.wikipedia.org/wiki/False_discovery_rate] for an history of this measure."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's denote : \n",
      "\n",
      "* $S$ the number of true positive\n",
      "* $V$ the number of false positives (Type I error)\n",
      "* $R$ is the number of rejected null hypotheses (also called \"discoveries\")\n",
      "* $Q$ the proportion of false discoveries amongst the discoveries = $\\frac{V}{S+V}$\n",
      "\n",
      "We want to control such that $Q = \\frac{V}{S+V} = \\frac{V}{R}$ is small. \n",
      "\n",
      "The false discovery rate is the expectation of the false discovery proportion $\\textrm{FDR} = E(Q)$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The procedure : step by step.\n",
      "\n",
      "First, we order the p-values by increasing order. We start to test the smallest p-value (the highest t or z), and declare the test significant if its p-value is less than $\\alpha/n$, for instance $0.05/n$. If this test is significant, we test the following pOur new threshold to declare significance is $\\alpha * 2 /n)$. This is iterated, until a test is not significant. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "T_{BH} = \\textrm{max} \\left \\{ { P_{(i)} : P_{(i)} < \\frac{\\alpha i }{n} } \\right  \\}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An intuition : $i  \\frac{\\alpha}{n} $ is the expected value if I draw $i$ times in a bernouilli distribution of probability $ \\frac{\\alpha}{n} $ "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's implement this and compare to Bonferroni or uncorrected thresolds:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats as st\n",
      "norm01 = st.norm(0,1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# adapted from nipy/algorithms/statistics/empirical_pvalue\n",
      "#q = ( (range(n)+1)/n) * 0.05 ) \n",
      "\n",
      "def fdr(p_values=None, verbose=0):\n",
      "    \"\"\"Returns the FDR associated with each p value\n",
      "\n",
      "    Parameters\n",
      "    -----------\n",
      "    p_values : ndarray of shape (n)\n",
      "        The samples p-value\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    q : array of shape(n)\n",
      "        The corresponding fdr values\n",
      "    \"\"\"\n",
      "    #p_values = check_p_values(p_values)\n",
      "    n_samples = p_values.size\n",
      "    order = p_values.argsort()\n",
      "    sp_values = p_values[order]\n",
      "\n",
      "    # compute q while in ascending order\n",
      "    q = np.minimum(1, n_samples * sp_values / np.arange(1, n_samples + 1))\n",
      "    for i in range(n_samples - 1, 0, - 1):\n",
      "        q[i - 1] = min(q[i], q[i - 1])\n",
      "\n",
      "    # reorder the results\n",
      "    inverse_order = np.arange(n_samples)\n",
      "    inverse_order[order] = np.arange(n_samples)\n",
      "    q = q[inverse_order]\n",
      "\n",
      "    if verbose:\n",
      "        import matplotlib.pylab as mp\n",
      "        mp.figure()\n",
      "        mp.xlabel('Input p-value')\n",
      "        mp.plot(p_values, q, '.')\n",
      "        mp.ylabel('Associated fdr')\n",
      "    return q\n",
      "\n",
      "\n",
      "def fdr_threshold(p_values, alpha=0.05):\n",
      "    \"\"\"Return FDR threshold given p values\n",
      "\n",
      "    Parameters\n",
      "    -----------\n",
      "    p_values : array of shape (n), optional\n",
      "        The samples p-value\n",
      "    alpha : float, optional\n",
      "        The desired FDR significance\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    critical_p_value: float\n",
      "        The p value corresponding to the FDR alpha\n",
      "    \"\"\"\n",
      "    n_samples = np.size(p_values)\n",
      "    sp_values = np.sort(p_values)\n",
      "    p_corr = alpha / n_samples\n",
      "    i = np.arange(1, n_samples + 1) # i from 1 to n_sample\n",
      "    critical_set = sp_values[sp_values < p_corr * i]\n",
      "    if len(critical_set) > 0:\n",
      "        critical_p_value = critical_set.max()\n",
      "    else:\n",
      "        critical_p_value = p_corr\n",
      "    return critical_p_value\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 1000\n",
      "unif = st.uniform()\n",
      "p = unif.rvs(size=n)\n",
      "\n",
      "q = []\n",
      "alph = arange(0.01,0.1,0.01)\n",
      "for a in alph:\n",
      "    q.append(fdr_threshold(p,a))\n",
      "\n",
      "    \n",
      "#more elegant: \n",
      "q = [fdr_threshold(p,a) for a in alph]\n",
      "    \n",
      "print size(alph), size(q)   \n",
      "plot(alph, q,'--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9 9\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "[<matplotlib.lines.Line2D at 0x3b30a10>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAD9CAYAAAB5lZr/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UVPeZP/D3JIztqpuSWBk3M2THzgwM+GOkJYynuxpa\nAwSbTmxNCHhs5kRi3W2UjVqrPU03+kcEdrcnyR5rtPm6FpKeQNL2CBakbJuS6uqAP0i7hyENTYd0\nmAFchWmik+by4/n+QZgA4gAycAd4v86Zc7j3fj53njsx95l77/P5jEZEBERERJNwm9oBEBHRzMdk\nQkREk8ZkQkREk8ZkQkREk8ZkQkREk8ZkQkREkzZmMqmpqYHVaoXFYkFxcfGobQoKCmCxWGCz2dDY\n2Dhm366uLmRkZCAhIQGZmZkIBAKhbYWFhbBYLLBaraitrQ2tLy8vh81mw/Lly7Fv375bOlgiIpoi\nEkZvb6+YTCbxeDyiKIrYbDZxu93D2lRVVUl2draIiLhcLrHb7WP23bNnjxQXF4uISFFRkezdu1dE\nRJqamsRms4miKOLxeMRkMkl/f79cuXJF7rnnHrly5YqIiDidTvn1r38dLnQiIppGYa9MGhoaYDab\nYTQaodVqkZubi4qKimFtKisr4XQ6AQB2ux2BQAAdHR1h+w7t43Q6ceLECQBARUUF8vLyoNVqYTQa\nYTabUV9fjz/96U+wWCxYtGgRAGDdunX42c9+FtmsSkREtywm3Eafz4f4+PjQssFgQH19/ZhtfD4f\n/H7/Tft2dnZCp9MBAHQ6HTo7OwEAfr8fq1evHtbH7/fjy1/+Mv7whz/gvffeg16vx4kTJ9DT0zMs\nDo1GM6EDJyKiARKBiVDCXpmM9wQ9nkBEZNT9aTSaMd8nNjYWL774Ih599FGsXbsWS5cuxe233z7q\ne0Tb65lnnlE9BsbEmOZiXIxpfK9ICXtlotfr4fV6Q8terxcGgyFsm7a2NhgMBvT09NywXq/XAxi4\nGuno6MCSJUvQ3t6OuLi4m+5rsM+DDz6IBx98EADwox/9CDExYUMnIqJpFPbKJDU1FS0tLWhtbYWi\nKCgvL4fD4RjWxuFwoLS0FADgcrkQGxsLnU4Xtq/D4UBJSQkAoKSkBBs2bAitLysrg6Io8Hg8aGlp\nQVpaGgDg8uXLAIDu7m68+OKLeOKJJyL4MRARRVYwCPT3qx3FNJIxVFdXS0JCgphMJjl48KCIiBw5\nckSOHDkSavPkk0+KyWSSlStXysWLF8P2FRG5evWqrFu3TiwWi2RkZEh3d3do27PPPismk0kSExOl\npqYmtD4vL0+Sk5MlOTlZysvLb4hzHIeiit/85jdqh3ADxjQ+jGn8ojEuNWNqahJZtkzkxInh66Px\nc4rUuVPz8c5mPI1GE9H7f0REt6K0FNi9GygqArZsAaK9NihS504+eCAiioBgENi+HTh3DnjjDWDF\nCrUjml5MJkREEbBjB9DTA5w/DyxcqHY004+3uYiIIuD6dWD+/Oi/rTVSpM6dTCZERHNYpM6dnDWY\niIgmjcmEiGgCSkuBzZvVjiL68AE8EdE4BIMDD9nPngVee03taKIPr0yIiMbQ3AykpQGKMlCtNdfK\nfseDyYSIKIyLF4G1a4FduwZucc3Fst/xYDUXEVEYvb3Au+8CiYlqRzI1WBo8ApMJEdHEsTSYiIii\nBpMJEREGqrV27wY+/rULmiAmEyKa89zugWqty5cHpkShiWMyIaI5rbQUuO8+YOdOVmtNBgctEtGc\nJAJs3QqcOTM3p4yPNF6ZENGcpNEAWVnAhQtMJJHA0mAiojmMpcFERBQ1xkwmNTU1sFqtsFgsKC4u\nHrVNQUEBLBYLbDYbGhsbx+zb1dWFjIwMJCQkIDMzE4FAILStsLAQFosFVqsVtbW1ofXHjx/HihUr\nYLPZkJ2djatXr97SARPR3ON2Ay6X2lHMchJGb2+vmEwm8Xg8oiiK2Gw2cbvdw9pUVVVJdna2iIi4\nXC6x2+1j9t2zZ48UFxeLiEhRUZHs3btXRESamprEZrOJoiji8XjEZDJJf3+/fPTRR3LXXXfJ1atX\nRUTkO9/5juzfv39YHGMcChHNUSUlIp/9rMirr6odSXSK1LkzbDVXQ0MDzGYzjEYjACA3NxcVFRVI\nSkoKtamsrITT6QQA2O12BAIBdHR0wOPx3LRvZWUl3nzzTQCA0+lEeno6ioqKUFFRgby8PGi1WhiN\nRpjNZjQ0NODee+/FnXfeiWvXruHOO+/E+++/D4vFckO8+/fvD/2dnp6O9PT0W06yRDSzBYPA9u3A\nuXOs1hqqrq4OdXV1Ed9v2GTi8/kQHx8fWjYYDKivrx+zjc/ng9/vv2nfzs5O6HQ6AIBOp0NnZycA\nwO/3Y/Xq1cP6tLW1wW6344UXXsDy5cuxcOFCJCQk4Ic//OEN8Q5NJkQ0dzU3A488AqSkDEwZz7Ej\nnxj5RfvAgQMR2W/YZyYajWZcO5FxVAKIyKj702g0Yd9Ho9Hg/fffR0FBAX73u9/B7/djxYoVKCws\nHFdsRDT3vPce8NRTHIQ4ncJemej1eni93tCy1+uFwWAI26atrQ0GgwE9PT03rNfr9QAGrkY6Ojqw\nZMkStLe3Iy4u7qb70uv1aG5uxtKlS7F06VIAwCOPPHLTYgAiogceUDuCuSfslUlqaipaWlrQ2toK\nRVFQXl4Oh8MxrI3D4UBpaSkAwOVyITY2FjqdLmxfh8OBkpISAEBJSQk2bNgQWl9WVgZFUeDxeNDS\n0oK0tDR87nOfw9tvv40rV64AAP77v/8bycnJkf0kiIjoloW9MomJicGhQ4eQlZWFvr4+5OfnIykp\nCUePHgUAbNu2DevXr0d1dTXMZjMWLFiA48ePh+0LAPv27UNOTg6OHTsGo9GI1z7+QeXk5GTk5OQg\nOTkZMTExOHz4MDQaDRYvXoyDBw/iS1/6Em677TYYjUb8+Mc/nsKPhYhmiqYmYNkytaMgjoAnohkp\nGAR27AAaGgZ+WnfePLUjmpk4Ap6I5qzBKeMVZaD0l4lEfUwmRDSjDE4Zv2sXq7WiCaegJ6IZo7sb\neOklDkKMRnxmQkQ0h/GZCRERRQ0mEyKKSsEg0N+vdhQ0XkwmRBR1Bqu1Tp5UOxIaLyYTIooqg9Va\nO3cCIybcoCjGai4iigqcMn5mYzIhoqiwYwfQ08Mp42cqlgYTUVS4fh2YPx8Y5y9fUIRE6tzJZEJE\nNIdxnAkREUUNJhMimlalpcDmzWpHQZHGB/BENC0Gp4w/exb4+CeMaBbhlQkRTbnm5k+mjD9/nmW/\nsxGTCRFNqYsXgbVrOWX8bMdqLiKaUr29wLvvAomJakdCo2Fp8AhMJkREEzdtpcE1NTWwWq2wWCwo\nLi4etU1BQQEsFgtsNhsaGxvH7NvV1YWMjAwkJCQgMzMTgUAgtK2wsBAWiwVWqxW1tbUAgA8++AAp\nKSmh1+LFi7Fz585bPmgiIoowCaO3t1dMJpN4PB5RFEVsNpu43e5hbaqqqiQ7O1tERFwul9jt9jH7\n7tmzR4qLi0VEpKioSPbu3SsiIk1NTWKz2URRFPF4PGIymaSvr++GuL7whS/I6dOnh60b41CIaIpd\nvy6ya5dIZ6fakdBEROrcGbY0uKGhAWazGUajEQCQm5uLiooKJCUlhdpUVlbC6XQCAOx2OwKBADo6\nOuDxeG7at7KyEm+++SYAwOl0Ij09HUVFRaioqEBeXh60Wi2MRiPMZjMaGhqwevXq0Pu98847uHz5\nMv7xH//xhnj3798f+js9PR3p6ekTz65ENGHNzcAjjwApKQNTolD0qqurQ11dXcT3GzaZ+Hw+xMfH\nh5YNBgPq6+vHbOPz+eD3+2/at7OzEzqdDgCg0+nQ2dkJAPD7/cMSx+C+hiorK0Nubu6o8Q5NJkQ0\nPUpLgd27gaIiYMsWzq0V7UZ+0T5w4EBE9hs2mWjG+a9CxvHwRkRG3Z9Gown7PiO3lZeX45VXXhlX\nXEQ0dUSArVuBM2c4ZTyN8QBer9fD6/WGlr1eLwwGQ9g2bW1tMBgMo67X6/UABq5GOjo6AADt7e2I\ni4u76b4G+wDA7373O/T29iIlJWXCB0pEkaXRAFlZwIULTCQ0RjJJTU1FS0sLWltboSgKysvL4Rjx\n02cOhwOlpaUAAJfLhdjYWOh0urB9HQ4HSkpKAAAlJSXYsGFDaH1ZWRkURYHH40FLSwvS0tJC7/Xq\nq69i06ZNkTt6IpqURx7hIEQaEPY2V0xMDA4dOoSsrCz09fUhPz8fSUlJOHr0KABg27ZtWL9+Paqr\nq2E2m7FgwQIcP348bF8A2LdvH3JycnDs2DEYjUa89vFEPcnJycjJyUFycjJiYmJw+PDhYbe5Xn/9\ndZw6dWpKPggiIrp1HLRIRGNyu4H33weG1MfQLMHfMyGiaVFaCtx3H9DaqnYkFM04BT0RjSoYBLZv\nB86dY7UWjY1XJkR0g8Ep43t6OGU8jQ+fmRDRDWpqAJ+PgxDnAs4aPAKTCRHRxPEBPBERRQ0mE6I5\nrqlJ7QhoNmAyIZqjgkEgPx/IzR34bXaiyWAyIZqD3O6Bai1FGSj9nTdP7YhopmMyIZpjBgch7to1\n8Dfn1qJI4KBFojmkuxt46SUOQqTIY2kwEdEcxtJgIiKKGkwmRLNUMAj096sdBc0VTCZEs9BgtdbJ\nk2pHQnMFkwnRLDNYrbVzJzDih1GJpgyruYhmCU4ZT2piMiGaJXbs+GTKeI4doek25m2umpoaWK1W\nWCwWFBcXj9qmoKAAFosFNpsNjY2NY/bt6upCRkYGEhISkJmZiUAgENpWWFgIi8UCq9WK2tra0HpF\nUfDNb34TiYmJSEpKws9//vNbOmCi2eo//5ODEElFEkZvb6+YTCbxeDyiKIrYbDZxu93D2lRVVUl2\ndraIiLhcLrHb7WP23bNnjxQXF4uISFFRkezdu1dERJqamsRms4miKOLxeMRkMkl/f7+IiPzrv/6r\nfP/73w+975UrV4bFMcahEBHRKCJ17gx7ZdLQ0ACz2Qyj0QitVovc3FxUVFQMa1NZWQmn0wkAsNvt\nCAQC6OjoCNt3aB+n04kTJ04AACoqKpCXlwetVguj0Qiz2YyGhgYAwPHjx/Hd73439L6LFi2KSDIl\nIqLJC/vMxOfzIT4+PrRsMBhQX18/Zhufzwe/33/Tvp2dndDpdAAAnU6Hzs5OAIDf78fq1atv2Nfg\nbbCnn34adXV1MJlMOHToEOLi4obFsn///tDf6enpSE9PH/MDIJppSkuB2lrglVfUjoRmorq6OtTV\n1UV8v2GTiWacv9cp4xiKLyKj7k+j0Yz5Pr29vWhra8M//MM/4Ac/+AGee+45fPvb30ZpaemwdkOT\nCdFsM1itdfYs8PrrakdDM9XIL9oHDhyIyH7D3ubS6/Xwer2hZa/XC4PBELZNW1sbDAbDqOv1ej2A\ngauRjo4OAEB7e3voCuNmfRYtWoT58+fj61//OgDg4YcfxqVLl27pgIlmosFBiD09wIULLPul6BM2\nmaSmpqKlpQWtra1QFAXl5eVwjBgF5XA4QlcILpcLsbGx0Ol0Yfs6HA6UlJQAAEpKSrBhw4bQ+rKy\nMiiKAo/Hg5aWFqSlpUGj0eCrX/0qfvOb3wAAfv3rX2PZsmWR/SSIotTFi5wynmaAsZ7QV1dXS0JC\ngphMJjl48KCIiBw5ckSOHDkSavPkk0+KyWSSlStXysWLF8P2FRG5evWqrFu3TiwWi2RkZEh3d3do\n27PPPismk0kSExOlpqYmtP69996TtWvXysqVK+X+++8Xr9c7LM5xHArRjNTTI/L222pHQbNVpM6d\nnIKeiGgO4xT0REQUNZhMiKJEMAjs3g1cvqx2JEQTx2RCFAUGq7UuXwbmz1c7GqKJYzIhUtnQKeNZ\nrUUzFWcNJlKJCLB1K3DmDKeMp5mPVyZEKtFogKwsDkKk2YGlwUREcxhLg4mIKGowmRBNA7cbcLnU\njoJo6jCZEE2xwWqt1la1IyGaOqzmIpoig1PGnzvHai2a/XhlQjQFmps/mTL+/HkmEpr9WM1FNAVq\nagCfD9iyZaAEmChaRercyWRCRDSHsTSYiIiiBpMJ0SQ1NakdAZH6mEyIblEwCOTnA7m5gKKoHQ2R\nuphMiG7B4JTxijJQ+jtvntoREamLyYRoggYHIe7axSnjiQaNmUxqampgtVphsVhQXFw8apuCggJY\nLBbYbDY0NjaO2berqwsZGRlISEhAZmYmAoFAaFthYSEsFgusVitqa2tD69PT02G1WpGSkoKUlBRc\nuXLllg6YaDK6u4GXXhoYhMiyX6IhJIze3l4xmUzi8XhEURSx2WzidruHtamqqpLs7GwREXG5XGK3\n28fsu2fPHikuLhYRkaKiItm7d6+IiDQ1NYnNZhNFUcTj8YjJZJL+/n4REUlPT5eLFy/eNNYxDoWI\niEYRqXNn2CuThoYGmM1mGI1GaLVa5ObmoqKiYlibyspKOJ1OAIDdbkcgEEBHR0fYvkP7OJ1OnDhx\nAgBQUVGBvLw8aLVaGI1GmM1m1NfXD018kcqhREQUQWHn5vL5fIiPjw8tGwyGYSf3m7Xx+Xzw+/03\n7dvZ2QmdTgcA0Ol06OzsBAD4/X6sXr16WB+/3x9adjqd0Gq12LhxI55++ukb4t2/f3/o7/T0dKSn\np4c7PKKwgkHg058GbuOTRZpF6urqUFdXF/H9hk0mmnHeEB7PFYOIjLo/jUYzrvf5yU9+grvvvhvX\nrl3Dxo0b8fLLL+Mb3/jGsDZDkwnRZLjdQE4O8OyzwEMPqR0NUeSM/KJ94MCBiOw37HcuvV4Pr9cb\nWvZ6vTAYDGHbtLW1wWAwjLper9cDGLga6ejoAAC0t7cjLi7upvsa7HP33XcDABYuXIhNmzahoaFh\n4kdLNA6D1Vo7dwIOh9rREM0MYZNJamoqWlpa0NraCkVRUF5eDseI/7scDgdKS0sBAC6XC7GxsdDp\ndGH7OhwOlJSUAABKSkqwYcOG0PqysjIoigKPx4OWlhakpaWhr68vVL3V09ODkydPYgWnYaUICwYH\nKrQKCweqtfLzWa1FNF5hb3PFxMTg0KFDyMrKQl9fH/Lz85GUlISjR48CALZt24b169ejuroaZrMZ\nCxYswPHjx8P2BYB9+/YhJycHx44dg9FoxGuvvQYASE5ORk5ODpKTkxETE4PDhw9Do9Hgr3/9Kx54\n4AH09PSgr68PGRkZ2Lp161R+LjQH7djxyZTxHDtCNDGcNZjoY9evA/Pn82qE5hZOQT8CkwkR0cRx\nCnoiIooaTCY055SWAps3qx0F0ewS9gE80WwSDA48ZD97Fvi45oOIIoRXJjQnNDd/MmX8+fMAK8uJ\nIovJhGa9ixeBtWs5ZTzRVGI1F816vb3Au+8CiYlqR0IUfVgaPAKTCRHRxLE0mIiIogaTCc0awSCw\nezdw+bLakRDNPUwmNCu43QPVWpcvD0yJQkTTi8mEZryhU8azWotIHRy0SDOWCLB1K3DmzMCU8Rw7\nQqQeXpnQjKXRAFlZwIULTCREamNpMBHRHMbSYCIiihpMJjQjuN2Ay6V2FER0M0wmFPUGq7VaW9WO\nhIhuhtVcFLWCQWD7duDcOVZrEUW7Ma9MampqYLVaYbFYUFxcPGqbgoICWCwW2Gw2NDY2jtm3q6sL\nGRkZSEhIQGZmJgKBQGhbYWEhLBYLrFYramtrb3gvh8OBFTyrzHqDU8b39HDKeKIZQcLo7e0Vk8kk\nHo9HFEURm80mbrd7WJuqqirJzs4WERGXyyV2u33Mvnv27JHi4mIRESkqKpK9e/eKiEhTU5PYbDZR\nFEU8Ho+YTCbp6+sLvdfPfvYz2bRpk6xYseKGWMc4FJphTp0S+X//T6S/X+1IiGa3SJ07w16ZNDQ0\nwGw2w2g0QqvVIjc3FxUVFcPaVFZWwul0AgDsdjsCgQA6OjrC9h3ax+l04sSJEwCAiooK5OXlQavV\nwmg0wmw2o6GhAQBw7do1PPfcc3j66adZAjwHPPAAkJ8/MJaEiKJf2GcmPp8P8fHxoWWDwYD6+vox\n2/h8Pvj9/pv27ezshE6nAwDodDp0dnYCAPx+P1avXj2sj9/vBwB8//vfx7e//W3MDzPx0v79+0N/\np6enIz09PdzhERHNOXV1dairq4v4fsMmE804vxaO50pBREbdn0ajCfs+IoK33noLf/rTn/Dcc8+h\nNUxJz9BkQjNHUxOwbJnaURDNDSO/aB84cCAi+w17m0uv18Pr9YaWvV4vDAZD2DZtbW0wGAyjrtfr\n9QAGrkY6OjoAAO3t7YiLiwu7L5fLhQsXLmDp0qVYs2YN3nnnHXz5y1++1WOmKBEMDtzKys0d+G12\nIpq5wiaT1NRUtLS0oLW1FYqioLy8HA6HY1gbh8OB0tJSAIDL5UJsbCx0Ol3Yvg6HAyUlJQCAkpIS\nbNiwIbS+rKwMiqLA4/GgpaUFaWlp+Kd/+if4fD54PB6cOXMGCQkJeOONNyL+YdD0GZwyXlEGSn/n\nzVM7IiKajLC3uWJiYnDo0CFkZWWhr68P+fn5SEpKwtGjRwEA27Ztw/r161FdXQ2z2YwFCxbg+PHj\nYfsCwL59+5CTk4Njx47BaDTitddeAwAkJycjJycHycnJiImJweHDh2+4BXaz22U0c5SWDvyIVXEx\n8PjjfMhONBtwokeaVt3dgMMBHD7MsSNE0SBS504mEyKiOYyzBhMRUdRgMqEpEwwC/f1qR0FE04HJ\nhKbEYLXWyZNqR0JE04HJhCJucMr4nTsHHrYT0ezHKegpYjhlPNHcxWRCEbNjxydTxi9cqHY0RDSd\nWBpMEXP9OjB/PgchEs0kHGcyApMJEdHEcZwJERFFDSYTmrDSUmDzZrWjIKJowgfwNG7B4MBD9rNn\ngY/n5iQiAsArExqn5uZPpow/f55lv0Q0HJMJjeniRWDtWmDXroFbXCz7JaKRWM1FY+rtBd59F0hM\nVDsSIoo0lgaPwGRCRDRxLA0mIqKowWRCIcHgwM/pXr6sdiRENNMwmRCAT6aMv3x5YEoUIqKJGDOZ\n1NTUwGq1wmKxoLi4eNQ2BQUFsFgssNlsaGxsHLNvV1cXMjIykJCQgMzMTAQCgdC2wsJCWCwWWK1W\n1NbWhtY/8MADWLVqFZYtW4b8/Hz09PTc0gHTjYZOGc9qLSK6JRJGb2+vmEwm8Xg8oiiK2Gw2cbvd\nw9pUVVVJdna2iIi4XC6x2+1j9t2zZ48UFxeLiEhRUZHs3btXRESamprEZrOJoiji8XjEZDJJf3+/\niIh88MEHoffcuHGjvPzyy8PiGONQaBT9/SL5+SKJiSK//73a0RCRGiJ17gx7ZdLQ0ACz2Qyj0Qit\nVovc3FxUVFQMa1NZWQmn0wkAsNvtCAQC6OjoCNt3aB+n04kTJ04AACoqKpCXlwetVguj0Qiz2Yz6\n+noAwMKPvy739PRAURR89rOfjVhCnas0GiArC7hwgYMQiWhywk6n4vP5EB8fH1o2GAyhk3u4Nj6f\nD36//6Z9Ozs7odPpAAA6nQ6dnZ0AAL/fj9WrV9+wr0FZWVk4f/48MjIy8MADD9wQ7/79+0N/p6en\nIz09PdzhEYBHHlE7AiKaTnV1dairq4v4fsMmE804f5hCxlGjLCKj7k+j0YR9n6HbfvnLX+Kjjz7C\no48+ipKSktDVzaChyYSIiG408ov2gQMHIrLfsLe59Ho9vF5vaNnr9cJgMIRt09bWBoPBMOp6vV4P\nYOBqpKOjAwDQ3t6OuLi4m+5rsM+gT33qU9i4cSPOnz8/oQOd69xuwOVSOwoimq3CJpPU1FS0tLSg\ntbUViqKgvLwcDodjWBuHw4HS0lIAgMvlQmxsLHQ6Xdi+DocDJSUlAICSkhJs2LAhtL6srAyKosDj\n8aClpQVpaWm4fv062tvbAQC9vb34xS9+gZSUlMh+ErPYYLVWa6vakRDRbBX2NldMTAwOHTqErKws\n9PX1IT8/H0lJSTh69CgAYNu2bVi/fj2qq6thNpuxYMECHD9+PGxfANi3bx9ycnJw7NgxGI1GvPbx\nfObJycnIyclBcnIyYmJicPjwYWg0Gly/fh0PPfQQPvroI4gIsrKysGXLlqn8XGaFYBDYvh04dw54\n4w0+ZCeiqcO5uWap5uaBh+spKcCLL3LsCBGNjhM9jsBkMlxNDeDzAVu2DJQAExGNhslkBCYTIqKJ\n46zBREQUNZhMZoGmJrUjIKK5jslkBgsGgfx8IDd34LfZiYjUwmQyQw1OGa8oA6W/8+apHRERzWVM\nJjPQ4CDEXbs4ZTwRRYewgxYp+nR3Ay+9xEGIRBRdWBpMRDSHsTSYiIiiBpNJFAsGgf5+taMgIhob\nk0mUGqzWOnlS7UiIiMbGZBKFBqu1du4ERsz4T0QUlVjNFUU4ZTwRzVRMJlFkxw6gpwc4f55jR4ho\nZmFpcBS5fh2YP59TxhPR9OEU9CPMhmRCRDTdOM6EiIiiBpOJCkpLgc2b1Y6CiChyxkwmNTU1sFqt\nsFgsKC4uHrVNQUEBLBYLbDYbGhsbx+zb1dWFjIwMJCQkIDMzE4FAILStsLAQFosFVqsVtbW1AIAP\nP/wQX/nKV5CUlITly5fju9/97i0fsJoGp4wvLAT27lU7GiKiCJIwent7xWQyicfjEUVRxGazidvt\nHtamqqpKsrOzRUTE5XKJ3W4fs++ePXukuLhYRESKiopk7969IiLS1NQkNptNFEURj8cjJpNJ+vv7\nJRgMSl1dnYiIKIoia9askVOnTg2LY4xDUZ3bLbJsmcjmzSIffKB2NEREAyJ17gx7ZdLQ0ACz2Qyj\n0QitVovc3FxUVFQMa1NZWQmn0wkAsNvtCAQC6OjoCNt3aB+n04kTJ04AACoqKpCXlwetVguj0Qiz\n2Yz6+nr8zd/8De677z4AgFarxec//3n4fL5I5tQpdfEisHYtp4wnotkr7DgTn8+H+Pj40LLBYEB9\nff2YbXw+H/x+/037dnZ2QqfTAQB0Oh06OzsBAH6/H6tXr75hX0MFAgGcPHkSTz311A3x7t+/P/R3\neno60tPTwx3etLHZgDNngMREtSMhormurq4OdXV1Ed9v2GSiGeeABxlHWZmIjLo/jUYT9n2Gbuvt\n7UVeXh7FoXN3AAALmklEQVT+5V/+BUaj8Ya2Q5NJNImJYSIhougw8ov2gQMHIrLfsLe59Ho9vF5v\naNnr9cJgMIRt09bWBoPBMOp6vV4PYOBqpKOjAwDQ3t6OuLi4m+5rsA8AfPOb30RiYiIKCgomfKBE\nRDR1wiaT1NRUtLS0oLW1FYqioLy8HI4RMw86HA6UlpYCAFwuF2JjY6HT6cL2dTgcKCkpAQCUlJRg\nw4YNofVlZWVQFAUejwctLS1IS0sDADz99NN4//338dxzz0X2E4igYBDYvRu4fFntSIiIptlYT+ir\nq6slISFBTCaTHDx4UEREjhw5IkeOHAm1efLJJ8VkMsnKlSvl4sWLYfuKiFy9elXWrVsnFotFMjIy\npLu7O7Tt2WefFZPJJImJiVJTUyMiIl6vVzQajSQnJ8uqVatk1apVcuzYsWFxjuNQplRTE6u1iGjm\nidS5k9OpREBp6cAVSVERsGUL59YiopkjUudOzho8CSLA1q0DlVqcMp6I5jJOpzIJGg2QlQVcuMBE\nQkRzG29zERHNYZw1mIiIogaTyTi53YDLpXYURETRiclkHEpLgfvuA1pb1Y6EiCg6sZorjGAQ2L4d\nOHeO1VpEROHwyuQmmpuBtDSgpwc4f56JhIgoHFZz3URNDeDzcRAiEc1ukTp3MpkQEc1hLA0mIqKo\nwWQCoKlJ7QiIiGa2OZ1MgkEgPx/IzQUURe1oiIhmrjmbTNzugWotRRko/Z03T+2IiIhmrjmZTAYH\nIe7aNfD3woVqR0RENLPNuUGL3d3ASy9xECIRUSSxNJiIaA5jaTAREUWNWZ1MPvwQ6O9XN4a6ujp1\nAxgFYxofxjR+0RgXY5peYyaTmpoaWK1WWCwWFBcXj9qmoKAAFosFNpsNjY2NY/bt6upCRkYGEhIS\nkJmZiUAgENpWWFgIi8UCq9WK2tra0Prvfe97uOeee/C3f/u34zowtxu4917g5MlxNZ8y0fiPhzGN\nD2Mav2iMizFNr7DJpK+vD9u3b0dNTQ3cbjdeffVVNDc3D2tTXV2NP/7xj2hpacGPfvQj/PM///OY\nfYuKipCRkYF33nkH69atQ1FREQDA7XajvLwcbrcbNTU1+Na3vhW6l/fQQw+hoaFhXAc1WK21cyfg\ncEzsAyEiookLm0waGhpgNpthNBqh1WqRm5uLioqKYW0qKyvhdDoBAHa7HYFAAB0dHWH7Du3jdDpx\n4sQJAEBFRQXy8vKg1WphNBphNptRX18PAEhLS8OSJUvCHszgIMSDBweqtfLzOUkjEdG0kDBef/11\neeKJJ0LLL7/8smzfvn1YmwcffFD+53/+J7S8bt06uXDhgvz0pz+9ad/Y2NjQ+v7+/tDy9u3b5ZVX\nXglty8/Pl5/+9KfD3m/hwoWjxgqAL7744ouvW3hFQthxJppxfq2XcZSVicio+9NoNGHfJ5IxEBHR\n1Ah7m0uv18Pr9YaWvV4vDAZD2DZtbW0wGAyjrtfr9QAAnU6Hjo4OAEB7ezvi4uJuuq/BPkREFL3C\nJpPU1FS0tLSgtbUViqKgvLwcjhFPtB0OB0pLSwEALpcLsbGx0Ol0Yfs6HA6UlJQAAEpKSrBhw4bQ\n+rKyMiiKAo/Hg5aWFqSlpUX8oImIKMLGug9WXV0tCQkJYjKZ5ODBgyIicuTIETly5EiozZNPPikm\nk0lWrlwpFy9eDNtXROTq1auybt06sVgskpGRId3d3aFtzz77rJhMJklMTJSamprQ+j179ojBYJDb\nb79dDAaDHDhwYBJ394iIKJIi8+RlGpw6dUoSExPFbDZLUVHRqG127NghZrNZVq5cKZcuXQqtf/zx\nxyUuLk6WL18eFTH9+c9/lvT0dElOTpZly5bJCy+8oHpMH374oaSlpYnNZpOkpCTZt2+f6jEN6u3t\nlVWrVsmDDz4YFTH9/d//vaxYsUJWrVol9957b1TE1N3dLRs3bhSr1SpJSUly7tw51eN6++23ZdWq\nVaHXHXfcEbF/65P5rA4ePCjJycmyfPlyycvLk7/+9a+qx/T888/L8uXLZdmyZfL8889HJJ7xxNTc\n3CyrV6+WT33qU/If//EfEz6eoWZEMunt7RWTySQej0cURRGbzSZut3tYm6qqKsnOzhYREZfLJXa7\nPbTtt7/9rVy6dCmiyWQyMbW3t0tjY6OIiHzwwQeSkJBwQ9/pjklE5Pr16yIi0tPTI3a7XU6fPq16\nTCIiP/jBD2TTpk3y1a9+ddLxRCImo9EoV69ejUgskYrpsccek2PHjonIwH+/QCAQFXEN6uvrkyVL\nlsif//xnVWPyeDyydOnSUALJycmRH//4x6rG9L//+7+yfPly+fDDD6W3t1fuv/9++eMf/zgtMV2+\nfFnOnz8v3/ve94Ylk/H0HWlGTKcymfEuALBmzRrceeedURFTZ2cnlixZglWrVgEAFi5ciKSkJPj9\nflVjAoD58+cDABRFQV9fH+666y7VY2pra0N1dTWeeOKJiFXsTTYmIPLVg5OJ6S9/+QtOnz6NLVu2\nAABiYmLwmc98RvW4hvrVr34Fk8mE+Ph4VWO64447oNVqEQwG0dvbi2AwGJEin8mco5qbm2G32/Hp\nT38at99+O+677z78/Oc/n5aYFi9ejNTUVGi12gn3HWlGJBOfzzfsH6HBYIDP55twm2iIqa2tbVib\n1tZWNDY2wm63qx5TX18fVq1aBZ1Ohy996UtITk5WLabBNjt37sS///u/47bbIvdPdbIxaTQa3H//\n/UhNTcVLL72kakxtbW3weDxYvHgxHn/8cXz+85/H1q1bEQwGVY9rqLKyMmzatEnVmHw+H+666y7s\n3r0b99xzD+6++27Exsbi/vvvVy0mv9+PFStW4PTp0+jq6kIwGERVVdUNn99UxRTJvjMimdzqWJPx\n9rsVkYjp2rVrePjhh/HCCy9gYQR+oWuyMd1+++1466230NbWht/+9rcRmUfoVmMSEfziF79AXFwc\nUlJSInolMNmxS2fOnEFjYyNOnTqFH/7whzh9+rRqMWk0GvT29uLSpUv41re+hUuXLmHBggWhKYrU\njGuQoig4efIkHnnkEVVjAoB3330Xzz//PFpbW+H3+3Ht2jX85Cc/UTUmq9WKvXv3IjMzE9nZ2UhJ\nSYnIl6fJnP9upe+MSCa3Ot5lKseoTDamnp4ebNy4EZs3bw6VRqsd06DPfOYz+MpXvoILFy6oGtPZ\ns2dRWVmJpUuXIi8vD2+88QYee+wxVWMCgLvvvhvAwC2Cr33ta+OeM26qYjIYDDAYDLj33nsBAA8/\n/DAuXbo06ZgmG9egU6dO4Qtf+AIWL16sekwXLlzAF7/4RSxatAgxMTH4+te/jrNnz6oaEwBs2bIF\nFy5cwJtvvonY2FgkJiZOS0wR7TvppzzToKenRz73uc+Jx+ORjz76aMyHW+fOnbvhIaDH44noA/jJ\nxNTf3y/f+MY35KmnnopYPJON6f/+7/9CJdrBYFDWrFkjv/rVr1SNaai6urqIVXNNJqbr16/L+++/\nLyIi165dky9+8Yvyy1/+UtWYRETWrFkjf/jDH0RE5JlnnpHvfOc7k44pEnGJiDz66KMRecgdiZga\nGxtl2bJlEgwGpb+/Xx577DE5dOiQqjGJiHR2doqIyHvvvSdWq1X+8pe/TEtMg5555plhD+An0nfQ\njEgmIpMb75Kbmyt/93d/J/PmzRODwSD/9V//pWpMp0+fFo1GIzabLVQ2eerUKVVj+v3vfy8pKSli\ns9lkxYoV8m//9m8RiWcyMQ1VV1cXsWquycT07rvvis1mE5vNJsuWLRs2fkqtmERE3nrrLUlNTZWV\nK1fK1772tYhVc002rmvXrsmiRYtCCTgaYiouLg6VBj/22GOiKIrqMa1Zs0aSk5PFZrPJG2+8EZF4\nxhNTe3u7GAwGueOOOyQ2Nlbi4+Plgw8+uGnfcGbNz/YSEZF6ZsQzEyIiim5MJkRENGlMJkRENGlM\nJkRENGlMJkRENGlMJkRENGn/H7zpkxCk0EBEAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x39c8890>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 1000\n",
      "z = norm01.rvs(size=n)\n",
      "p = norm01.sf(z)\n",
      "\n",
      "print fdr_threshold(p)\n",
      "print sum(p < fdr_threshold(p))\n",
      "\n",
      "q = fdr(p, verbose=1)\n",
      "\n",
      "if len(z) <= 20:\n",
      "    order = z.argsort()    \n",
      "    print \"z \\n\", z \n",
      "    print \"p \\n\", p \n",
      "    print \"q \\n\", q\n",
      "    print \"sorted z \\n\", z[order]\n",
      "    print \"sorted fdr \\n\", q[order]\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "5e-05\n",
        "0\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEMCAYAAADal/HVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUlHX+B/D3w8yIF0BELYGhQGEZTBgwLuVPFE3E3GJj\nkyNurWkq5uZmbrp53NNJts2y0/XEtkdctbSF0HTDSsmldRTSmkTFjniQSFykvCUKSQoMz++Ppxmu\ncwHmPu/XOXNm5rl+5jn6fPg+35sgiqIIIiIiI7wcHQARETk3JgoiIjKJiYKIiExioiAiIpOYKIiI\nyCQmCiIiMsnqieLxxx/H7bffjujoaKPbPPXUU4iIiIBarcbx48cNy4uLi6FSqRAREYENGzZYOzQi\nIuoHqyeKhQsXori42Oj6vXv34ttvv0V1dTXy8vKwbNkyAIBOp8Py5ctRXFyMyspKFBQU4PTp09YO\nj4iI+sjqiSI5ORkjRowwun7Pnj147LHHAABJSUm4du0aLly4AK1Wi/DwcISGhkKhUCArKwtFRUXW\nDo+IiPpIbu8T1tfXIyQkxPBdqVSivr4e33//fY/lX331VY/9BUGwS5xERO6mvwNxOKQye6Cjhoii\nyJco4vnnn3d4DM7y4rXgteC16PqKjBQhk4kA9K/+s3uiCA4ORl1dneH7+fPnoVQqeyyvq6uDUqm0\nd3hERC4nOxuQywFB6HhVVQE6nXWOb/dEkZ6ejm3btgEAvvzyS/j7++P2229HfHw8qqurUVtbi5aW\nFhQWFiI9Pd3e4RERuRR/f2DTJtNJYdCggZ3D6nUU8+bNw8GDB3HlyhWEhIQgJycHra2tAIClS5di\n9uzZ2Lt3L8LDwzFs2DBs3bpVCkQuR25uLtLS0qDT6bBo0SJERUVZOzy3kpKS4ugQnAavRQdeiw7u\nfi38/YHr101vs28fMGuWVMroL0EURZcaZlwQBLhYyERENiGXdy1JKBTAsGHAjRuAlxfw9deAvkvb\nQO6ddm/1REREA5OdDZw5A7S3dyzbuROYM8c252OiICJyMWfOAAcPSp+9vIATJzpKDrbAREFE5EJU\nKqCmRvocHQ0cOiTVVdgSBwUkInIhFy4AbW3S5x9/tH2SAJgoiIicUm99IwSho5XT4MHA4cP2iYWJ\ngojICZ05Y7pvhCgCd95pn1iYKIiInNDQocbXCYLU9NVe2I+CiMiJ6Ju+KhRSiaK0VKqTEATAzw+4\ndQvQavveyon9KIiIXFx2NrBlS9fHTZmZwC8DWzgUHz0RETmB7nUSggDk5Tkuns5YoiAicjB/f6Cx\nseuy0lL7NH21BEsUREQ2YqyJa29NXjtXH5w8Cfzf/zku7u5YoiAispBKJc3zYEsnT9p2OI7+YImC\niMhCFy7Y9vhlZc6XJAAmCiLyYNnZQFCQ+UdD3XtF95VMBowaJb0CAqRj+fh0rJ88GWhocK7HTZ0x\nURCRR8rOBnbsAH74oe/7CgIwfLj0HhfXMSlQdLQ0muurr3ZsW1Ym9YO4fFl6/fijNDx4U5NULyGK\nzlVx3Rt2uCMijxQU1L8kYct5H2xpIPdOliiIyCPdvNnxOSAA8PYGRoyQ5pfWv/v6dt3HVZPEQLHV\nExG5NXPzSttrTgdXxhIFEbm1n34yvX7cOCYJc1iiICK31Nu80t3dcQewdav9YnJVLFEQkVvSzyst\nilKrJH9/qUWSSiWtj48HKipYmrAESxRE5Jb08zkkJAD79zMhDASbxxKRW7p2TXr8lJfHJAEM7N7J\nREFE5AHYj4KI6BcqlVSCGD0aOHfO0dG4B5YoiMgt9NZfQqkE6uocE4+zYYmCiDxeb/0lysrsH4c7\nYqIgIrfg1e1utnMncOedjonF3fDRExG5NH3HOp0O+PJLqYPdoUPOO2S3owzk3sl+FETk0vQd6wAg\nM1MaOpysi4+eiMhlZWdLU4cC0rwQeXmOjcddsURB5MGys4GPP7b9FJ/2EBjIjnW2wkRB5KH0M7z1\nd3pPZzNokKMjcF989ETkoc6ccZ8kMXYsR4G1JZskiuLiYqhUKkRERGDDhg091jc0NCAjIwNqtRpJ\nSUk4deqUYV1oaChiYmIQFxeHxMREW4RHROgYNG/4cCA5GVAoOr6PHCnN+JaQ0DHy6qhR0rIXXpC2\n8/eXmqB6eUkzxI0YIW3r4yOtHz++45h6I0dKxwSAYcOk9+hoaX8AkMmkz51nmRs8WFp3xx1SbHr6\n5rCTJwPl5XzsZEtWbx6r0+kQGRmJkpISBAcHIyEhAQUFBYiKijJss3r1avj5+eG5555DVVUVnnzy\nSZSUlAAAwsLCUF5ejgD9v5zuAbN5LJFVcNA8z+JUPbO1Wi3Cw8MRGhoKhUKBrKwsFBUVddnm9OnT\nmDZtGgAgMjIStbW1uHz5smE9EwGR7fn7S3UUTBJkjtUrs+vr6xESEmL4rlQq8dVXX3XZRq1WY/fu\n3Zg8eTK0Wi3OnTuH8+fPY/To0RAEATNmzIBMJsPSpUuxZMmSHudYt26d4XNKSgpSUlKs/TOIiFya\nRqOBRqOxyrGsnigEQTC7zZo1a7BixQrExcUhOjoacXFxkMlkAICysjIEBQXh8uXLSE1NhUqlQnJy\ncpf9OycKIuobfU/moUOB/HyWKNxV9z+ic3Jy+n0sqyeK4OBg1HUarrGurg5KpbLLNr6+vtiyZYvh\ne1hYGMaOHQsACAoKAgCMHj0aGRkZ0Gq1PRIFEVkmOxvYskUa3sLYevZkJnOsXkcRHx+P6upq1NbW\noqWlBYWFhUhPT++yzfXr19HS0gIA2LRpE6ZOnQofHx80NzejqakJAHDjxg3s378f0dHR1g6RyGPo\nx0Ayhj2ZyRJWL1HI5XLk5uYiLS0NOp0OixYtQlRUFDZu3AgAWLp0KSorK7FgwQIIgoAJEyZg8+bN\nAICLFy8iIyMDANDW1oZHHnkEM2fOtHaIRB4hOxv45hvj6199lY+dyDIcPZbITaWkdAyWp+/rcNtt\nQGUl8M47wLJlDg2P7IyjxxIRsrOB998Hfv6563KVCjhyhKUH6j+WKIicXG9TfPaFl5fpegryDE7V\n4Y6IrKu3KT774tAh68RBnouJgsjJdZ/isy/KyjjTGw0cEwWRE/L3lwbYEwRpak+FQhoQT9+f9Z57\npAH29N8FQRpQTxCk5CCK0otJgqyBdRRETkgu71qvMGQI0NzsuHjI9bGOgsiFZWdLj5f0JQhB6Jok\nBAHoNlwakV2xREHkYEFBwA8/GF9/8qQ0ZwPRQLBEQeTCbt7sfbkgMEmQc2CHOyIHys6WKqv19LPE\nNTcDWi2TBDkHJgoiB1GpgG+/7aiP+M1vgI8+cmxMRL1hHQWRg3Tuca1QAJcucZgNsh3WURC5IIVC\nepfJgPJyJglyXkwURA5y9CigVAI1NayLIOfGOgoiB/D3l8Zw8vICGhsdHQ2RaayjIHKAzj2v2eua\n7IF1FEQuRj/Qn5cXe12T82OiIHKA8nKpJHHiBOsnyPmxjoLIjjrXTZSXM0mQa2CJgsgOsrOBoUOl\nfhM6HdDaCiQlOToqIsswURDZwSef9JzLmnUT5CqYKIjsoPvAf/v28bETuQ4mCiI7uPtu6V0ul0aE\nnTXLsfEQ9QUTBZEd7NwJZGYCly+zJEGuhx3uiIg8wEDunWwe2weDBkmtVVyZTNZ1mk2yH4WCTWLJ\nNbFE0adzO+S05EY4XAc5is2G8NDpdHjjjTf6dWB3xERBAyEIbBJLrslkopDJZMjPz7dXLE6vtFT6\nz+7vD3h7A9OmSa/Bg4H4eGndqFHAvfdK2ysU0vfbbpPWA1KrF/27IEjLvb2BkSOBgADp8VZvhg+X\ntps8WdoOkPbvnLzi4jrmOPDz61guk0nLU1OBsjLpOApFRyzdt5fLpd+VmgqMHi39vp07O44NAPfc\nA9x+u/T7Ro6U4u48fpH+2CNGdCz38emI199f2nfQIOlYgwZ1XC/98Xx9O453221Sk9JBg6RjREZK\n64YO7XqNBg+WttNfI/1yQFp2223S+Tpft97+AIiNla6Tft/OfHx6LpPLpd/a/VoqFNLywYOBigo+\ndiLXZPbR08qVK9Ha2oq5c+di2LBhhuUTJ060eXC9YWU2EVHfDeTeaTZRpKSkQOjlT64DBw7064QD\nxURBRNR3Nk0UzoaJgoio72zSPPa1114zHLw3f/rTn/p1QleUnQ2cOSM9D8/P59zGRORZjCaKpqYm\nCIKAqqoqfP3110hPT4coivjkk0+QmJhozxgd7swZ4OBB6XN2NrBjh2PjISKyJ7OPnpKTk7F37174\n/tIEpampCbNnz0ZpaaldAuzOEY+eZs+WWtIkJAD797NEQUSux6ZToV66dAmKTu0iFQoFLl261K+T\nuar8fGmcHiYJIvJEZhPF/PnzkZiYiHXr1uH5559HUlISHnvsMZP7FBcXQ6VSISIiAhs2bOixvqGh\nARkZGVCr1UhKSsKpU6cs3tcR/P2lx01MEkTkiYw+evruu+8wduxYAEB5eTlKS0shCAKmTJmCuLg4\nowfU6XSIjIxESUkJgoODkZCQgIKCAkRFRRm2Wb16Nfz8/PDcc8+hqqoKTz75JEpKSizal62eiIj6\nziatnjIzM1FeXo777rsPn3/+Oe7WD6hvhlarRXh4OEJDQwEAWVlZKCoq6nKzP336NNasWQMAiIyM\nRG1tLS5duoSamhqz+xIRkX0ZTRQ6nQ4vvvgiqqqq8Prrr3fJRIIgGG0eW19fj5CQEMN3pVKJr7oN\ncKNWq7F7925MnjwZWq0W586dw/nz5y3aFwDWrVtn+JySkoKUlBSzP5SIyJNoNBpoNBqrHMtoovjg\ngw/w0UcfQafToampyeIDGut30dmaNWuwYsUKxMXFITo6GnFxcZDJZBbtC3RNFERE1FP3P6JzcnL6\nfSyjiUKlUmHNmjWIiYnB7NmzLT5gcHAw6urqDN/r6uqgVCq7bOPr64stW7YYvoeFhWHcuHH4+eef\nze5LRET2ZbbVU1+SBADEx8ejuroatbW1aGlpQWFhIdLT07tsc/36dbS0tAAANm3ahKlTp8LHx8ei\nfYmIyL6sPsOdXC5Hbm4u0tLSoNPpsGjRIkRFRWHjxo0AgKVLl6KyshILFiyAIAiYMGECNm/ebHJf\nIiJyHA4KSETkAWzSPHbXrl2GA/dWyfzb3/62XyckIiLXYjRRfPzxxxAEAZcuXcLhw4cxffp0ANI8\nFJMmTWKiICLyEEYTxbvvvgsASE1NRWVlJQIDAwEAP/zwg9khPIiIyH2YbfVUV1eHMWPGGL7ffvvt\n+N///mfToIiIyHmYbfU0Y8YMpKWl4Xe/+x1EUURhYSFSU1PtERsRETkBs62eRFHEv//9b8P8E1Om\nTEFGRoZdgusNWz0REfWdTVo9dT74xIkT4evri9TUVDQ3N6OpqckwkZG7U6mACxcAhQI4ehS4805H\nR0REZF9m6yjy8vKQmZmJJ554AgBw/vx5PPTQQzYPzFlcuABcvw5cuQJMnuzoaIiI7M9sovj73/+O\nsrIy+Pn5AQB+9atfedQMd7duSe9eXsDevY6NhYjIEcwmCm9vb3h7exu+t7W1WTzKqztQq6X39nbg\nhRccGwsRkSOYTRRTp07Fiy++iObmZvznP/9BZmYmHnzwQXvE5hQCAqT3hAQgL8+xsRAROYLZVk/t\n7e345z//if379wMA0tLSsHjxYoeVKuzd6unaNSA7W0oSnDObiFzVQO6dZhPFW2+9hRUrVphdZi9s\nHktE1HcDuXeaffSkH8qjs61bt/brZERE5HqM9qMoKChAfn4+zp4926VOoqmpCSNHjrRLcERE5HhG\nE8WkSZMQGBiIy5cvY9WqVYYii6+vL9T6pkBEROT2OHEREZEHsGkdxZEjR5CQkAAfHx8oFAp4eXkZ\nOt8REZH7M5soli9fjvz8fERERODmzZvYvHkz/vCHP9gjNiIicgJmEwUAREREQKfTQSaTYeHChSgu\nLrZ1XERE5CTMjh47bNgw3Lp1C2q1Gn/+858xZswY1hEQEXkQsyWKbdu2ob29Hbm5uRg6dCjOnz+P\nXbt22SM2IiJyAmz1RETkAWwycVFmZiZ27tyJ6OjoXk948uTJfp2QiIhci9ESxffff4+goCDU1tb2\numNoaKgNwzKOJQoior6zSYkiKCgIgDRn9pgxYzBkyBAAwM8//4yLFy/262REROR6zFZmz5kzBzKZ\nrGMHLy/MmTPHpkEREZHzMJsodDodBg0aZPju7e2N1tZWmwZFRETOw2yiGDVqFIqKigzfi4qKMGrU\nKJsGRUREzsNs89hvv/0WjzzyCL7//nsAgFKpxPbt2xEeHm6XALtjZTYRUd/ZdIY7vZ9++gkA4OPj\n068TWQsTBRFR39l09Nhr165h5cqVmDp1KqZOnYpnnnkG169f79fJiIjI9ZhNFI8//jj8/Pywc+dO\n7NixA76+vli4cKE9YiMiIidg9tGTWq1GRUWF2WX2Yq9HT9nZwJkzwNChQH4+4O9v81MSEdmMTTrc\n6Q0ZMgSlpaVITk4GAJSVlWHo0KH9Opkz8/cHjD1Re/BBoLTUvvEQETkLsyWKEydOYP78+YZ6iREj\nRuC9995z2LzZtipRyOWATmd8PevPiciV2aXV0/Xr1yEIgkXToBYXF+Ppp5+GTqfD4sWL8eyzz3ZZ\nf+XKFTz66KO4cOEC2trasGrVKixYsACANIaUn58fZDIZFAoFtFpt14BtlCgGDQKM9SPctw+YNcvq\npyQishubtnp688030djYCD8/P6xcuRITJ07EZ599ZnR7nU6H5cuXo7i4GJWVlSgoKMDp06e7bJOb\nm4u4uDicOHECGo0GzzzzDNra2gw/RqPR4Pjx4z2ShC2VlwPe3oBCAYwYIb0DTBJERGYTxZYtW+Dn\n54f9+/fj6tWr2LZtG9asWWN0e61Wi/DwcISGhkKhUCArK6tLz24ACAwMRGNjIwCgsbERI0eOhFze\nUV3iiH4S0dHAzZtASwtw9ar0LopMEkREZiuz9TftTz/9FL///e8xYcIEk9vX19cjJCTE8F2pVOKr\nr77qss2SJUswffp0BAUFoampCTt27DCsEwQBM2bMgEwmw9KlS7FkyZIe51i3bp3hc0pKClJSUsz9\nDCIij6LRaKDRaKxyLLOJ4u6778bMmTPx3Xff4aWXXkJjYyO8vIwXRARBMHvS9evXIzY2FhqNBjU1\nNUhNTUVFRQV8fX3xxRdfIDAwEJcvX0ZqaipUKpWhxZVe50RBREQ9df8jOicnp9/HsujR00svvYSj\nR49i2LBhaG1txdatW41uHxwcjLq6OsP3uro6KJXKLtscPnwYmZmZAIBx48YhLCwMVVVVAKTHUgAw\nevRoZGRk2LWegoiIejKbKI4cOYLIyEj4+/tj+/bt+Nvf/obhw4cb3T4+Ph7V1dWora1FS0sLCgsL\nkZ6e3mUblUqFkpISAMDFixdRVVWFsWPHorm5GU1NTQCAGzduYP/+/b1OxUpERPZjNlE88cQTGDZs\nGCoqKvD6668jPDwc8+fPN7q9XC5Hbm4u0tLSMH78eMydOxdRUVHYuHEjNm7cCABYu3Ytjh49CrVa\njRkzZuCVV15BQEAALly4gOTkZMTGxiIpKQkPPPAAZs6cab1fS0REfWa2H0VcXByOHz+OnJwcBAcH\nY/HixZg4cSKOHTtmrxi74OixRER9Z9MhPHx9fbF+/Xq8//77KC0thU6n4wx3REQexOyjp8LCQnh7\ne2PLli0YM2YM6uvrsXr1anvEZjf+/tIQHoMGAd984+hoiIici8VDeOiVlpaioKAA77zzjq1iMskW\nj546j/M0ZAjQ3GzVwxMROZxNHz0BwLFjx1BQUIAdO3YgLCwMDz/8cL9O5qy8vKRE4eUFdOsbSETk\n8YwmiqqqKhQUFKCwsBCjR49GZmYmRFG0Wk8/Z1JeDiQlSUmCrXGJiLoymiiioqLwwAMP4LPPPsMd\nd9wBAHj99dftFpgteXn1Pmx4VRUTBRFRd0Yrs3fv3o0hQ4ZgypQpeOKJJ/D555+7TbNUYz/jl87i\nRETUidnK7J9++glFRUUoKCjAgQMHMH/+fGRkZDisI5w1KrONDUe1cycwZ86ADk1E5JTsMnERAFy9\nehUffvghPvjgA/z3v//t1wkHyhqJorgYuP9+6bOfH9DYyCRBRO7NbonCGbBnNhFR39l0hjsiIvJs\nTBRERGQSEwUREZnEREFERCYxURARkUlMFEREZBITBRERmcREQUREJjFREBGRSUwURERkksckiuxs\naSY7Ly9OeUpE1BcekyjOnJFmsRNFoLVVmqiIiIjM85hEMXRox2dB4JSnRESW8phEkZ8PzJwJDB4M\nVFRwJjsiIktxmHEiIg/AYcaJiMhmmCiIiMgkJgoiIjKJiYKIiExioiAiIpOYKIiIyCQmCiIiMomJ\ngoiITGKiICIik5goiIjIJJskiuLiYqhUKkRERGDDhg091l+5cgWzZs1CbGwsJkyYgHfffdfifYmI\nyL6sPtaTTqdDZGQkSkpKEBwcjISEBBQUFCAqKsqwzbp163Dr1i289NJLuHLlCiIjI3Hx4kUIgmB2\nX471RETUd0411pNWq0V4eDhCQ0OhUCiQlZWFoqKiLtsEBgaisbERANDY2IiRI0dCLpdbtC8REdmX\n3NoHrK+vR0hIiOG7UqnEV90mf1iyZAmmT5+OoKAgNDU1YceOHRbvC0glEr2UlBSkpKRY90cQEbk4\njUYDjUZjlWNZPVEIgmB2m/Xr1yM2NhYajQY1NTVITU1FRUWFxefonCiIiKin7n9E5+Tk9PtYVn/0\nFBwcjLq6OsP3uro6KJXKLtscPnwYmZmZAIBx48YhLCwMVVVVUCqVZvclIiL7snqiiI+PR3V1NWpr\na9HS0oLCwkKkp6d32UalUqGkpAQAcPHiRVRVVWHs2LEW7UtERPZl9UdPcrkcubm5SEtLg06nw6JF\nixAVFYWNGzcCAJYuXYq1a9di4cKFUKvVaG9vxyuvvIKAgAAA6HVfIiJyHE6FSkTkAZyqeSwREbkX\nJgoiIjKJiYKIiExioiAiIpOYKIiIyCQmCiIiMomJgoiITGKiICIik5goiIjIJCYKIiIyiYmCiIhM\nYqIgIiKTmCiIiMgkJgoiIjKJiYKIiExioiAiIpM8IlFkZwMpKcDs2cC1a46OhojItXhEojhzBjh4\nENi3T0oaRERkOY9IFEOHSu8JCUBenmNjISJyNR4xZ/a1a1JJIi8P8Pe3UWBERE5sIHNme0SiICLy\ndAO5d3rEoyciIuo/l0wUgtC/19ChwLlzjo6eiMi1uOSjJ6D/ISuVQF2d9eIhInIFfPRkIZkMKCtz\ndBRERK7FJROFl5fUeikwUOob4e0NyOXSOpms9318fYGaGuDOO+0XJxGRO3DJR08uFjIRkcPx0RMR\nEdkMEwUREZnEREFERCYxURARkUlMFEREZBITBRERmcREQUREJjFREBGRSUwULkyj0Tg6BKfBa9GB\n16IDr4V12CRRFBcXQ6VSISIiAhs2bOix/tVXX0VcXBzi4uIQHR0NuVyOa79MZh0aGoqYmBjExcUh\nMTHRFuG5Df4n6MBr0YHXogOvhXXIrX1AnU6H5cuXo6SkBMHBwUhISEB6ejqioqIM26xatQqrVq0C\nAHzyySd488034f/L1HOCIECj0SAgIMDaoRERUT9YvUSh1WoRHh6O0NBQKBQKZGVloaioyOj2+fn5\nmDdvXpdlHMuJiMiJiFa2c+dOcfHixYbv27dvF5cvX97rtjdu3BADAgLEhoYGw7KwsDAxNjZWvPvu\nu8W8vLwe+0CajIIvvvjii68+vvrL6o+epImFLPPxxx9j8uTJhsdOAPDFF18gMDAQly9fRmpqKlQq\nFZKTkw3rRZY2iIjsyuqPnoKDg1HXaQq5uro6KJXKXrf94IMPejx2CgwMBACMHj0aGRkZ0Gq11g6R\niIj6wOqJIj4+HtXV1aitrUVLSwsKCwuRnp7eY7vr16/j0KFD+M1vfmNY1tzcjKamJgDAjRs3sH//\nfkRHR1s7RCIi6gOrP3qSy+XIzc1FWloadDodFi1ahKioKGzcuBEAsHTpUgDARx99hLS0NAwZMsSw\n78WLF5GRkQEAaGtrwyOPPIKZM2daO0QiIuqLftdu2Ni+ffvEyMhIMTw8XHz55Zd73eaPf/yjGB4e\nLsbExIjHjh2zc4T2Y+5avP/++2JMTIwYHR0tTpo0SayoqHBAlPZhyb8LURRFrVYrymQycdeuXXaM\nzr4suRYHDhwQY2NjxbvuukucOnWqfQO0I3PX4vLly2JaWpqoVqvFu+66S9y6dav9g7SDhQsXirfd\ndps4YcIEo9v0577plImira1NHDdunHj27FmxpaVFVKvVYmVlZZdtPv30U/H+++8XRVEUv/zySzEp\nKckRodqcJdfi8OHD4rVr10RRlP7DePK10G83bdo08de//rX44YcfOiBS27PkWjQ0NIjjx48X6+rq\nRFGUbpbuyJJr8fzzz4tr1qwRRVG6DgEBAWJra6sjwrWpQ4cOiceOHTOaKPp733TKITws6YuxZ88e\nPPbYYwCApKQkXLt2DRcvXnREuDZlybW49957MXz4cADStTh//rwjQrU5S/vovP3225gzZw5Gjx7t\ngCjtw5JrkZ+fj4cfftjQmGTUqFGOCNXmLLkWgYGBaGxsBAA0NjZi5MiRkMut/uTd4ZKTkzFixAij\n6/t733TKRFFfX4+QkBDDd6VSifr6erPbuOMN0pJr0dnmzZsxe/Zse4Rmd5b+uygqKsKyZcsA9K25\ntiux5FpUV1fj6tWrmDZtGuLj47F9+3Z7h2kXllyLJUuW4NSpUwgKCoJarcZbb71l7zCdQn/vm06Z\nUi39zy1261PhjjeFvvymAwcOYMuWLfjiiy9sGJHjWHItnn76abz88ssQBAGi9GjVDpHZnyXXorW1\nFceOHcPnn3+O5uZm3HvvvbjnnnsQERFhhwjtx5JrsX79esTGxkKj0aCmpgapqamoqKiAr6+vHSJ0\nLv25bzplorCkL0b3bc6fP4/g4GC7xWgvlvZLOXnyJJYsWYLi4mKTRU9XZsm1KC8vR1ZWFgDgypUr\n2LdvHxSBbli+AAAFT0lEQVQKRa9NtF2ZJdciJCQEo0aNwpAhQzBkyBBMmTIFFRUVbpcoLLkWhw8f\nxl/+8hcAwLhx4xAWFoaqqirEx8fbNVZH6/d90yo1KFbW2toqjh07Vjx79qx469Yts5XZR44ccdsK\nXEuuxblz58Rx48aJR44ccVCU9mHJtehswYIFbtvqyZJrcfr0afG+++4T29raxBs3bogTJkwQT506\n5aCIbceSa7Fy5Upx3bp1oiiK4oULF8Tg4GDxxx9/dES4Nnf27FmLKrP7ct90yhKFJX0xZs+ejb17\n9yI8PBzDhg3D1q1bHRy1bVhyLf7617+ioaHB8FxeoVC4ZY92S/voeAJLroVKpcKsWbMQExMDLy8v\nLFmyBOPHj3dw5NZnybVYu3YtFi5cCLVajfb2drzyyituOUL1vHnzcPDgQVy5cgUhISHIyclBa2sr\ngIHdNwVRdNOHuEREZBVO2eqJiIicBxMFERGZxERBREQmMVEQEZFJTBTkVnx8fKx+zHPnzqGgoMDq\nxzVl3bp1eO211+x6TiJjmCjIrdiid/7Zs2eRn59v9eOa4o6jDJDrYqIgt6TRaJCSkoLMzExERUXh\n0UcfNawLDQ3Fs88+i5iYGCQlJaGmpgYAsGDBAuzatcuwnX54hzVr1qC0tBRxcXE9xgjSaDSYMmUK\nHnjgAahUKixbtqzHEAnXr19HaGio4fuNGzdwxx13oK2tDZs2bUJiYiJiY2MxZ84c/Pzzz4bt9Mki\nJSUF5eXlAKTe5mFhYQAAnU6H1atXIzExEWq1Gnl5eQO9bES9YqIgt3XixAm89dZbqKysxHfffYfD\nhw8DkG7A/v7+OHnyJJYvX46nn37asLw3GzZsQHJyMo4fP44VK1b0WP/1118jNzcXlZWVqKmpwe7d\nu7usHz58uGGcIQD45JNPMGvWLMjlcjz88MPQarU4ceIEoqKisHnz5h7HFwSh19g2b94Mf39/aLVa\naLVabNq0CbW1tX25REQWYaIgt5WYmIigoCAIgoDY2NguN1H9XO1ZWVk4cuSIyeOY65OamJiI0NBQ\neHl5Yd68eSgrK+uxzdy5c1FYWAhAmit+7ty5AIBvvvkGycnJiImJwb/+9S9UVlZa/Pv279+Pbdu2\nIS4uDvfccw+uXr2Kb7/91uL9iSzllEN4EFmDt7e34bNMJkNbW1uv2+n/WpfL5WhvbwcAtLe3o6Wl\nxaLzdP5rXxRFCIKAjz76CDk5OQCkv/wffPBBrF27Fg0NDTh27BimT58OQHrctWfPHkRHR+O9994z\nlDo66xzXzZs3u6zLzc1FamqqRXES9RdLFOSR9H/dFxYWYtKkSQCkugt9XcCePXsMY+T4+vqiqanJ\n6LG0Wi1qa2vR3t6OHTt2IDk5GQ899BCOHz+O48ePY+LEifDx8UFCQgKeeuopPPjgg4bk8tNPP2HM\nmDFobW3F+++/b1gudhoiPTQ0FEePHgUAfPjhh4bzpqWl4Z133jEkwDNnzqC5udlq14hIj4mC3Ern\nv+5NtRxqaGiAWq3G22+/jTfeeAOANLnNwYMHERsbiy+//NLQ1FatVkMmkyE2NrZHZbYgCEhISMDy\n5csxfvx4jB07Fg899FCv55w7dy7y8/MNj50A4IUXXkBSUhImT56MqKioLsfVx79q1Sr84x//wMSJ\nE/Hjjz8ali9evBjjx4/HxIkTER0djWXLlhktNRENBAcFJI8TFhaG8vJyq4weqtFo8Nprr+Hjjz+2\nQmREzoklCvI41uyjYKxFEpE7YYmCiIhMYomCiIhMYqIgIiKTmCiIiMgkJgoiIjKJiYKIiExioiAi\nIpP+H1QWxOtRE1nQAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x39f47d0>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alph = 0.05\n",
      "df = 14\n",
      "t14 = tdist(df)\n",
      "n = 1000\n",
      "ts = t14.rvs(size=n)\n",
      "# ts = np.sort(ts)           # careful of ts = ts.sort() .... you get None !!!\n",
      "ps = t14.sf(ts)\n",
      "\n",
      "#uncorrected :\n",
      "print  'Uncorrected: ',sum(ps < alph)\n",
      "\n",
      "#Bonferroni : \n",
      "print 'Bonferroni: ',sum(ps < alph_corr(alph,n))\n",
      "\n",
      "#FDR : (suppose z)\n",
      "fdr_alph = fdr_threshold(ps)\n",
      "print 'FDR: ',sum(ps < fdr_alph)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Uncorrected:  55\n",
        "Bonferroni:  0\n",
        "FDR:  0\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Adding signal :\n",
      "z = st.norm.rvs(0,1, size=1000)\n",
      "print sum(z>1.5)\n",
      "z = z[z>1.]+ 2\n",
      "\n",
      "alph, df, n = 0.05, 14, 100\n",
      "t14 = tdist(df)\n",
      "ts = t14.rvs(size=n)\n",
      "# add some signal : \n",
      "print \"number of true positive: \", sum(ts > 1.5)\n",
      "ts[ts > 1.5] += 1.5\n",
      "ps = t14.sf(ts)\n",
      "\n",
      "#uncorrected :\n",
      "print  'Uncorrected: ',alph, sum(ps < alph)\n",
      "#Bonferroni : \n",
      "print 'Bonferroni: ', alph_corr(alph,n), sum(ps < alph_corr(alph,n))\n",
      "#FDR :\n",
      "fdr_alph = fdr_threshold(ps)\n",
      "q = fdr(ps, verbose=1)\n",
      "print 'FDR: ',fdr_alph, sum(ps < fdr_alph), sum(q <= alph)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "60\n",
        "number of true positive:  11\n",
        "Uncorrected:  0.05 11\n",
        "Bonferroni:  0.000512801416262 0\n",
        "FDR: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.00473460436807 10 11\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEMCAYAAADEXsFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9QVXX+x/Hn5UcuKoqmqfxoUWEFEy6aaNqqVGu4lGQp\nI+5uP8jU3HFTd/vhtDMb7HcrbabaJrYdbNUyg7CssFJyLSnNH5goONkY/qCAWs0EtazA6/n+ceFe\nfl3vBe69/PD1mGG455zPPefNZ/S8z/mcz/l8TIZhGIiIyGXPp6MDEBGRzkEJQUREACUEERGpo4Qg\nIiKAEoKIiNRRQhAREcADCeHee+9l0KBBxMTEOCzzwAMPEBkZidlsZv/+/e4OQURE2sDtCSEtLY38\n/HyH2zdt2sSRI0coLS1l5cqVLFy40N0hiIhIG7g9IUyaNIl+/fo53L5x40buvvtuAMaPH091dTUn\nTpxwdxgiItJKft4+YGVlJWFhYbbl0NBQKioqGDRoUKNyJpPJ26GJiHQLbR2AokMeKjcN1tHJ3zAM\n/RgGjz32WIfH0Fl+VBeqC9XFpX/aw+sJISQkhPLycttyRUUFISEh3g5DRESa8HpCSE5OZu3atQDs\n3r2boKCgZs1FIiLifW5PCHPmzGHixIkcPnyYsLAwVq9eTVZWFllZWQAkJSUxbNgwIiIiWLBgAS+8\n8IK7Q+h2EhISOjqETkN1Yae6sOuqdTF/PvTsCf7+MGAAfPmldV1wMPTvD1OnQnW1dd2QIXDllfZ1\nnmAy2tvo5CEmk6nd7WEiIp1ZQgJ89JF9OTQUhg9vvC4lBU6ebL5u/fqW99mec6feVBYRaaf5860n\n96Sk1l299+zZ+POOHY3XjR4NK1e2vM4TdIcgItJK8+fDF19YT9TZ2TBjhv0K/lJX701VV8Pvfw8H\nDsDOnfDLX1rXpaWBYcBLL0FQkHXdPfeAyQRr1ljXOdKec6cSgoiIA01P/PUn4oZNPSkp8P33sHkz\nxMfDli2XPmF7WnvOnV5/MU1EpKv44gv7iX/+fPuVf30TTny8vflm/nzr545MBu2lOwSRy0D9le6x\nY3D11dCnDwwcaO3V0vTqNyrKWg7g+uvhrbdcO8k1vJoeONB6pVxTAwEBEB4OX31lP3bD4zX9/tGj\n1qaThuVbinX+fHj3XfjuO+v6nj1h6FD7/h9+uOWr+5YEBcGZM9bP/frB/v3WGJKSWr7yr+/50xkT\nQLvOnUYn1YlDE+lypkwxDGurtP1n4ED755QUe9m+fRuXa7jN1WM03HdLPy3ts6UYLxXrpcqnpDTe\n7uxv8PVt/P3QUOv6qirrd6uqXKuDzqA95071MhK5DNQ3cfTta/0dHw9ms/1zw14r/v72z7Gxrvdo\nadiMUr9vsF6xNz12S/tsGqOzWBv2vGl4nPoyLTXrOOLT4Ez4i19Ye/uA9ep//frOdxfgMW5MTG7V\niUMT6XLqr3TLyuxXvI6ufsvKDCM42DBuuaV1V8YN91dVZRi33WYYM2bYj9nw2K7E6CzWqirr/pOS\nGh+nvkxrru5LSgyjRw/DGDTIup+urD3nTj1DEBHpRvRimkgX0dILTI7WNR2+oOE+goOhR4+Wt4u0\nlRKCiBfVd2PcvNl6Yr/Uum++gaoq2LrVvr7htpqalreLtJXeQ5DLXsPuklddBWVll+6q2LS7Y69e\nLXervO46+N//rA9pP/3U2o2xpQedl1oHzYcqaPow1ZNDGcjlRc8Q5LLX8K3TgQPh22+tnx0NQdB0\nQLKWpKRY+63X920PDYXy8pb7rzta13T4gnr122pq4IornA9lIJcXDV0h0g4NXz7q29faBHOpIQjq\ny9cLDIRz56zfPXPG/t3ISDh1ynpFf+iQ9Q5BxNP0UFmkHbKz7Vf0r79u/+zoqrt+MLOkJOvvgwet\n3ykubvzdTz+13hkoGUhXoTsEEZFuRHcIIiLSbkoIIiICKCGIB7gye1RUlPXFKn9/a9nqauu6oCD7\nyJauHqd3b+sD3YED4Xe/s67z9wc/P2svnIMH3fe3iXRneoYgbtd08pCWum42HG64vlxL3TRdPU69\nK66wdsdsKCAAzp9vzV8g0nXpGYJ0Kq6MMtlwRM2YGGu5+nX1c8u6ehw/P/vyhAmNy/j4wJ49rscu\ncjlTQhC3a9iN01HXzU8/tY7Hk5QEH3/ctm6a9ccpKrJ/7+23ret27LDeGRw4YE04IuKcmoxERLoR\nNRmJiEi7KSGIiAighCAiInWUEEREBFBCEBGROkoI4jJX3kAWka5LCUFc1tJUjyLSfWgKTXGo4dSS\n2dmuvYEsIl2XEoI4VH9HANbkkJ3dfKpHEek+lBDEoaZ3BEFBLQ9UJyLdg4auEIdamvxdRDq39pw7\nlRBERLqRTjeWUX5+PlFRUURGRrJixYpm20+dOsW0adOIi4tj1KhRvPTSS54IQ0REWsHtdwgWi4UR\nI0awdetWQkJCiI+PJycnh+joaFuZ9PR0fv75Z5588klOnTrFiBEjOHHiBH5+9kcaukMQEWm9TnWH\nUFhYSEREBOHh4fj7+5OamkpeXl6jMkOGDOHs2bMAnD17liuvvLJRMhAREe9z+1m4srKSsLAw23Jo\naCh7mkxZNW/ePG688UaCg4M5d+4c6x10XUlPT7d9TkhIICEhwd3hioh0aQUFBRQUFLhlX25PCCaT\nyWmZJ554gri4OAoKCjh69ChTp06luLiYwMDARuUaJgQREWmu6cVyRkZGm/fl9iajkJAQyhvMjl5e\nXk5oaGijMjt37iQlJQWA4cOHM3ToUA4fPuzuUEREpBXcnhDGjh1LaWkpZWVl1NTUkJubS3JycqMy\nUVFRbN26FYATJ05w+PBhhg0b5u5QRESkFdzeZOTn50dmZiaJiYlYLBbmzp1LdHQ0WVlZACxYsIBH\nH32UtLQ0zGYzFy9e5KmnnqJ///7uDkVERFpBL6aJiHQjnarbqYiIdE1KCCIiAighXHY065mIOKKE\ncJnRrGci4ogSwmVGs56JiCPqZXSZ0RwHIt2b5kMQERFA3U5FRMQNlBBERARQQhARkTpKCCIiAigh\niIhIHSUEEREBPDD8tXjP/PnWN4979oSrroKyMjh6FH75S+jTB7Kz9a6BiLhO7yF0YQkJ1mEoAAYO\nhG+/bbw9JQUcTFctIt2U3kO4DMyfD8HB0L8/TJ1qfeO44TAUZrP1c9++9nUamkJEWkMJoYv44gv4\n5huoqoKtW60JIjvbehewZQu8/rr1c3GxfZ2ai0SkNdRk1EUkJVlHKAUYPRo+/FAnfBFpTmMZXQaq\nqyEtDQwDXnpJyUBEWqaEICIigB4qi4iIGyghiIgIoIQgIiJ1lBBERARQQhARkTpKCCIiAighiIhI\nnUsmBIvFwrPPPuutWEREpANdMiH4+vqSnZ3trVhERKQDOX1TeenSpdTW1jJ79mx69eplWz9mzBjP\nBqY3lUVEWs2jQ1ckJCRgMpmard+2bVubDugqJQQRkdbTWEYiIgK079zpcArNp59+2rbzlvz5z39u\n0wHFruEUmJruUkQ6msOEcO7cOUwmE4cPH2bv3r0kJydjGAbvvvsu48aN82aMXVrTk/7DD8O6dVBb\nax3K2mKxl9N0lyLSkZw2GU2aNIlNmzYRGBgIWBNFUlIS27dv92xg3aTJqOG8xykpcPKkfblefLxm\nOBMR9/Do8NcnT57E39/ftuzv78/JkyfbdLDObP5868k7Kck6GU29oCDw84MrroCDB63LPj7Wn8GD\n4de/hrAw6++rr7auczTv8cqV9mWAgADr8ZQMRKRTMJz4xz/+YcTExBiPPfaY8be//c2IjY01Hn/8\n8Ut+Z/PmzcaIESOMiIgIY/ny5S2W2bZtmxEXF2dcc801xpQpU5ptdyE0t5oyxTCsjTiGkZJiX+/r\na18fENB42dlPSophVFXZfxuG9fcttxhGcLBhlJV59U8UkctAe86dDpuMjh07xrBhwwDYt28f27dv\nx2QyMXnyZEaPHu0wwVgsFkaMGMHWrVsJCQkhPj6enJwcoqOjbWWqq6u5/vrref/99wkNDeXUqVMM\nGDCg0X683WRUP2dx0+abK66wtvf7+MCBA3DttdZlgF694IcfoG9fOHPG/hs077GIdAyPNBmlpKQA\ncNNNN3HttdeyZMkSFi9efMlkAFBYWEhERATh4eH4+/uTmppKXl5eozLZ2dnMnDmT0NBQgGbJoCNk\nZ1vb+Js23+zbZ23aOXAAYmKsy7/4Bdx8M3z2mfU7xcX23zNmwG23KRmISNfjsJeRxWLh8ccf5/Dh\nwzzzzDONMo7JZHLY7bSyspKwsDDbcmhoKHv27GlUprS0lNraWm644QbOnTvH4sWLufPOO5vtKz09\n3fY5ISGBhIQEV/+uVgsKarmXT0wMnD/fePnHH+3L9d+p//3WWx4LUUSkmYKCAgoKCtyyL4cJ4bXX\nXuPtt9/GYrFw7tw5l3fo6L2FhmpraykqKuKDDz7g/PnzTJgwgeuuu47IyMhG5RomBBERaa7pxXJG\nRkab9+UwIURFRbFs2TJiY2NJSkpyeYchISGUl5fblsvLy21NQ/XCwsIYMGAAAQEBBAQEMHnyZIqL\ni5slBBER8R6n3U5bkwwAxo4dS2lpKWVlZdTU1JCbm0tycnKjMrfddhs7duzAYrFw/vx59uzZw8iR\nI1sXuYiIuJXDO4Q279DPj8zMTBITE7FYLMydO5fo6GiysrIAWLBgAVFRUUybNo3Y2Fh8fHyYN2+e\nEoKISAfT4HYiIt2IRwa327Bhg23HLT0ovuOOO9p0QBER6ZwcJoR33nkHk8nEyZMn2blzJzfeeCNg\nnQdh4sSJSggiIt2Mw4Tw0ksvATB16lQOHTrEkCFDAPjmm2+4++67vRKciIh4j9NeRuXl5QwePNi2\nPGjQIL766iuPBiUiIt7ntJfRb37zGxITE/nd736HYRjk5uYydepUb8QmIiJe5LSXkWEYvPXWW7b5\nDyZPnsztt9/u+cDUy0hEpNU80suo4c7HjBlDYGAgU6dO5fz585w7d842YY6IiHQPTp8hrFy5kpSU\nFO6//34AKioqmDFjhscDExER73KaEP71r3+xY8cO+vTpA8CvfvWrbjljmojI5c5pQujRowc9evSw\nLV+4cMGlEU1FRKRrcZoQpkyZwuOPP8758+f573//S0pKCtOnT/dGbCIi4kVOexldvHiR//znP2zZ\nsgWAxMRE7rvvPo/fJaiXkYhI67Xn3Ok0ITz33HMsXrzY6Tp3U0IQEWk9j8ypXK9+CIuG1qxZ06aD\niYhI5+XwPYScnByys7M5fvx4o2cG586d48orr/RKcCIi4j0OE8LEiRMZMmQI3377LQ8++KDtFiQw\nMBCz2ey1AEVExDs0QY6ISDfi0WcIu3btIj4+nt69e+Pv74+Pj4/tJTUREek+nCaERYsWkZ2dTWRk\nJD/99BOrVq3ij3/8ozdiExERL3KaEAAiIyOxWCz4+vqSlpZGfn6+p+MSEREvczraaa9evfj5558x\nm808/PDDDB48WG37IiLdkNM7hLVr13Lx4kUyMzPp2bMnFRUVbNiwwRuxiYiIF6mXkYhIN+KRCXJS\nUlJ4/fXXiYmJafGAJSUlbTqgiIh0Tg7vEL7++muCg4MpKytr8Yvh4eEeDEt3CCIibeGRO4Tg4GDA\nOqfy4MGDCQgIAODHH3/kxIkTbTqYiIh0Xk4fKs+aNQtfX1/7F3x8mDVrlkeDEhER73OaECwWC1dc\ncYVtuUePHtTW1no0KBER8T6nCWHAgAHk5eXZlvPy8hgwYIBHgxIREe9z2u30yJEj/P73v+frr78G\nIDQ0lFdeeYWIiAjPBqaHyiIirebRGdPqff/99wD07t27TQdqLSUEEZHW8+hop9XV1SxdupQpU6Yw\nZcoU/vKXv3DmzJk2HUxERDovpwnh3nvvpU+fPrz++uusX7+ewMBA0tLSvBGbiIh4kdMmI7PZTHFx\nsdN1bg9MTUYiIq3m0SajgIAAtm/fblvesWMHPXv2bNPBRESk83J6h3DgwAHuuusu23ODfv368fLL\nL3t8XmXdIYiItJ5H7xDi4uIoKSmhpKSEgwcPcuDAAafJID8/n6ioKCIjI1mxYoXDcnv37sXPz483\n33yz9ZGLiIhbOU0I//znPzl79ix9+vRh6dKljBkzhvfff99heYvFwqJFi8jPz+fQoUPk5OTw+eef\nt1jukUceYdq0aboTEBHpBJwmhNWrV9OnTx+2bNnC6dOnWbt2LcuWLXNYvrCwkIiICMLDw/H39yc1\nNbXRm871nn/+eWbNmsXAgQPb9xeIiIhbOJ1Cs/7q/b333uPOO+9k1KhRlyxfWVlJWFiYbTk0NJQ9\ne/Y0K5OXl8eHH37I3r17MZlMLe4rPT3d9jkhIYGEhARn4YqIXFYKCgooKChwy76cJoRrr72Wm2++\nmWPHjvHkk09y9uxZfHwc31g4Ork3tGTJEpYvX257+OGoyahhQhARkeaaXixnZGS0eV9OE8Lq1avZ\nv38/w4cPp1evXnz33XesWbPGYfmQkBDKy8tty+Xl5YSGhjYqs2/fPlJTUwE4deoUmzdvxt/fn+Tk\n5Lb+HSIi0k5OE8KuXbswm8307t2bV155haKiIpYsWeKw/NixYyktLaWsrIzg4GByc3PJyclpVObY\nsWO2z2lpaUyfPl3JQESkgzl9qHz//ffTq1cviouLeeaZZ4iIiOCuu+5yWN7Pz4/MzEwSExMZOXIk\ns2fPJjo6mqysLLKystwavIiIuI/TF9NGjx7N/v37ycjIICQkhPvuu48xY8ZQVFTk2cD0YpqISKt5\nZE7leoGBgTzxxBOsW7eO7du3Y7FYNGOaiEg35LTJKDc3lx49erB69WoGDx5MZWUlDz30kDdiExER\nL3J5gpx627dvJycnhxdeeMFTMQFqMhIRaQuPNhkBFBUVkZOTw/r16xk6dCgzZ85s08FERKTzcpgQ\nDh8+TE5ODrm5uQwcOJCUlBQMw3DbG3EiItK5OGwy8vHx4dZbbyUzM5Orr74agKFDh3L8+HHvBKYm\nIxGRVvPI8NdvvvkmAQEBTJ48mfvvv58PPvhAJ2gRkW7M6UPl77//nry8PHJycti2bRt33XUXt99+\nOzfffLNnA9MdgohIq7Xn3NmqXkanT5/mjTfe4LXXXuPDDz9s0wFdpYQgItJ6XksI3qSEICLSeh6d\nQlNERC4PSggiIgIoIYiISB0lBBERAZQQRESkjhKCiIgASggiIlJHCUFERAAlBBERqaOEICIigBKC\niIjUUUIQERFACUFEROooIYiICKCEICIidZQQREQEUEIQEZE6SggiIgIoIYiISB0lBBERAZQQRESk\njhKCiIgASggiIlJHCUFERAAlBBERqeORhJCfn09UVBSRkZGsWLGi2fZXX30Vs9lMbGws119/PSUl\nJZ4IQ0REWsFkGIbhzh1aLBZGjBjB1q1bCQkJIT4+npycHKKjo21ldu3axciRI+nbty/5+fmkp6ez\ne/fuxoGZTLg5NBGRbq8950633yEUFhYSERFBeHg4/v7+pKamkpeX16jMhAkT6Nu3LwDjx4+noqLC\n3WGIiEgr+bl7h5WVlYSFhdmWQ0ND2bNnj8Pyq1atIikpqcVt6enpts8JCQkkJCS4K0wRkW6hoKCA\ngoICt+zL7QnBZDK5XHbbtm2sXr2aTz75pMXtDROCiIg01/RiOSMjo837cntCCAkJoby83LZcXl5O\naGhos3IlJSXMmzeP/Px8+vXr5+4wRESkldz+DGHs2LGUlpZSVlZGTU0Nubm5JCcnNyrz1Vdfcccd\nd7Bu3ToiIiLcHYKIiLSB2+8Q/Pz8yMzMJDExEYvFwty5c4mOjiYrKwuABQsW8Pe//52qqioWLlwI\ngL+/P4WFhe4ORUREWsHt3U7dRd1ORURar1N1OxURka5JCUFERAAlBBERqaOEICIigBKCiIjUUUIQ\nERFACUFEROooIYiICKCEICIidZQQREQEUEIQEZE6SggiIgIoIYiISB0lBBERAZQQRESkjhKCiIgA\nSggiIlJHCUFERAAlBBERqaOEICIigBKCiIjUUUIQERFACUFEROooIYiICKCEICIidZQQREQEUEIQ\nEZE6SggiIgJ08oQwfz4kJEBSElRXd3Q0IiLdm8kwDKOjg2iJyWTC19fAYrEup6TA+vUdG5OISGdn\nMplo62m9UycEsIbm7w8nT0JQUMfGJCLS2bUnIXTqJiMAX1/Yt0/JQETE0zp9QrBY4OzZjo5CRKT7\n6xJNRj4+2J4liIiIY926yQjg4487OgIRke6vUycEkwl27IDrr+/oSDpWQUFBR4fQaagu7FQXdqoL\n9/BIQsjPzycqKorIyEhWrFjRYpkHHniAyMhIzGYz+/fvb7HMxYtKBqB/7A2pLuxUF3aqC/dwe0Kw\nWCwsWrSI/Px8Dh06RE5ODp9//nmjMps2beLIkSOUlpaycuVKFi5c6O4wRESkldyeEAoLC4mIiCA8\nPBx/f39SU1PJy8trVGbjxo3cfffdAIwfP57q6mpOnDjh7lBERKQV/Ny9w8rKSsLCwmzLoaGh7Nmz\nx2mZiooKBg0a1KictaeRAGRkZHR0CJ2G6sJOdWGnumg/tycEV0/iTbtFNf1eJ+0NKyLSbbm9ySgk\nJITy8nLbcnl5OaGhoZcsU1FRQUhIiLtDERGRVnB7Qhg7diylpaWUlZVRU1NDbm4uycnJjcokJyez\ndu1aAHbv3k1QUFCz5iIREfEutzcZ+fn5kZmZSWJiIhaLhblz5xIdHU1WVhYACxYsICkpiU2bNhER\nEUGvXr1Ys2aNu8MQEZHWMjrY5s2bjREjRhgRERHG8uXLWyzzpz/9yYiIiDBiY2ONoqIiL0foPc7q\nYt26dUZsbKwRExNjTJw40SguLu6AKL3DlX8XhmEYhYWFhq+vr7FhwwYvRuddrtTFtm3bjLi4OOOa\na64xpkyZ4t0AvchZXXz77bdGYmKiYTabjWuuucZYs2aN94P0grS0NOOqq64yRo0a5bBMW86bHZoQ\nLly4YAwfPtw4fvy4UVNTY5jNZuPQoUONyrz33nvGb3/7W8MwDGP37t3G+PHjOyJUj3OlLnbu3GlU\nV1cbhmH9j3E510V9uRtuuMG45ZZbjDfeeKMDIvU8V+qiqqrKGDlypFFeXm4YhvWk2B25UhePPfaY\nsWzZMsMwrPXQv39/o7a2tiPC9aiPP/7YKCoqcpgQ2nre7NChK/TOgp0rdTFhwgT69u0LWOuioqKi\nI0L1OFfqAuD5559n1qxZDBw4sAOi9A5X6iI7O5uZM2faOm8MGDCgI0L1OFfqYsiQIZytGx757Nmz\nXHnllfj5ub1lvMNNmjSJfv36Odze1vNmhyaElt5HqKysdFqmO54IXamLhlatWkVSUpI3QvM6V/9d\n5OXl2d5y767vrLhSF6WlpZw+fZobbriBsWPH8sorr3g7TK9wpS7mzZvHZ599RnBwMGazmeeee87b\nYXYKbT1vdmjqdNc7C91Ba/6mbdu2sXr1aj755BMPRtRxXKmLJUuWsHz5cttQv03/jXQXrtRFbW0t\nRUVFfPDBB5w/f54JEyZw3XXXERkZ6YUIvceVunjiiSeIi4ujoKCAo0ePMnXqVIqLiwkMDPRChJ1L\nW86bHZoQ9M6CnSt1AVBSUsK8efPIz8+/5C1jV+ZKXezbt4/U1FQATp06xebNm/H392/Wxbmrc6Uu\nwsLCGDBgAAEBAQQEBDB58mSKi4u7XUJwpS527tzJX//6VwCGDx/O0KFDOXz4MGPHjvVqrB2tzedN\ntzzhaKPa2lpj2LBhxvHjx42ff/7Z6UPlXbt2ddsHqa7UxZdffmkMHz7c2LVrVwdF6R2u1EVD99xz\nT7ftZeRKXXz++efGTTfdZFy4cMH44YcfjFGjRhmfffZZB0XsOa7UxdKlS4309HTDMAzjf//7nxES\nEmJ89913HRGuxx0/ftylh8qtOW926B2C3lmwc6Uu/v73v1NVVWVrN/f396ewsLAjw/YIV+ricuFK\nXURFRTFt2jRiY2Px8fFh3rx5jBw5soMjdz9X6uLRRx8lLS0Ns9nMxYsXeeqpp+jfv38HR+5+c+bM\n4aOPPuLUqVOEhYWRkZFBbW0t0L7zZqedQlNERLyrU8+YJiIi3qOEICIigBKCiIjUUUIQERFACUG6\nqN69e7t9n19++SU5OTlu3++lpKen8/TTT3v1mCKOKCFIl+SJt9WPHz9Odna22/d7Kd3xrXvpupQQ\npEsrKCggISGBlJQUoqOj+cMf/mDbFh4eziOPPEJsbCzjx4/n6NGjANxzzz1s2LDBVq5+WINly5ax\nfft2Ro8e3WwMnIKCAiZPnsytt95KVFQUCxcubDY0wJkzZwgPD7ct//DDD1x99dVcuHCBF198kXHj\nxhEXF8esWbP48ccfbeXqk0JCQgL79u0DrG9fDx06FACLxcJDDz3EuHHjMJvNrFy5sr3VJtIiJQTp\n8g4cOMBzzz3HoUOHOHbsGDt37gSsJ9qgoCBKSkpYtGgRS5Yssa1vyYoVK5g0aRL79+9n8eLFzbbv\n3buXzMxMDh06xNGjR3nzzTcbbe/bt69tHB2Ad999l2nTpuHn58fMmTMpLCzkwIEDREdHs2rVqmb7\nN5lMLca2atUqgoKCKCwspLCwkBdffJGysrLWVJGIS5QQpMsbN24cwcHBmEwm4uLiGp0s58yZA0Bq\naiq7du265H6cvaM5btw4wsPD8fHxYc6cOezYsaNZmdmzZ5ObmwvAa6+9xuzZswE4ePAgkyZNIjY2\nlldffZVDhw65/Pdt2bKFtWvXMnr0aK677jpOnz7NkSNHXP6+iKu630Dhctnp0aOH7bOvry8XLlxo\nsVz91befnx8XL14E4OLFi9TU1Lh0nIZX74ZhYDKZePvtt8nIyACsV/LTp0/n0UcfpaqqiqKiIm68\n8UbA2ky1ceNGYmJiePnll213EQ01jOunn35qtC0zM5OpU6e6FKdIW+kOQbq1+qv13NxcJk6cCFif\nLdS31W/cuNE2BkxgYCDnzp1zuK/CwkLKysq4ePEi69evZ9KkScyYMYP9+/ezf/9+xowZQ+/evYmP\nj+eBBx5Q2zhdAAABMElEQVRg+vTptiTy/fffM3jwYGpra1m3bp1tvdFg6O7w8HA+/fRTAN544w3b\ncRMTE3nhhRdsie6LL77g/PnzbqsjkXpKCNIlNbxav1RPnaqqKsxmM88//zzPPvssYJ1E5aOPPiIu\nLo7du3fburCazWZ8fX2Ji4tr9lDZZDIRHx/PokWLGDlyJMOGDWPGjBktHnP27NlkZ2fbmosA/u//\n/o/x48fz61//mujo6Eb7rY//wQcf5N///jdjxozhu+++s62/7777GDlyJGPGjCEmJoaFCxc6vAsS\naQ8Nbifd1tChQ9m3b59bRrssKCjg6aef5p133nFDZCKdk+4QpNtyZx9/Rz2ARLoT3SGIiAigOwQR\nEamjhCAiIoASgoiI1FFCEBERQAlBRETqKCGIiAgA/w+8X6gsFA3t9wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3d9cc10>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Permutations : a quick introduction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The principle of the permutation is simple: we construct a distribution of the statistics we are interested in under the assumption that there is no signal. For instance, let's take the one sample t-test scenario."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#N subjects, z[i] = effects observed for subject i.\n",
      "\n",
      "N = 15\n",
      "mu, sigma = 0.75, 1.0\n",
      "normH1 = st.norm(mu, sigma)\n",
      "z = normH1.rvs(size=N) # a series of z values under H1\n",
      "\n",
      "print z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-0.21557732  1.08753659  0.81284445  0.65051999 -0.48500761  0.71379641\n",
        " -0.93627376  0.91966905  1.44485769  0.79344869 -0.59309978  0.13385896\n",
        "  1.3633538   0.06926551  0.28891933]\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The idea here is the following: if the z are drawn from a normal distribution of mean zero, then there there is as much chance to get a negative value than a positive values (the distribution is symmetric). We can therefore arbitrarily swap the sign of our values and construct the distribution of their mean under this hypothesis. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psign = np.random.randint(2, size=len(z))\n",
      "psign[psign==0] = -1\n",
      "print psign\n",
      "print z * psign"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1  1 -1  1  1 -1 -1 -1 -1 -1  1  1 -1 -1  1]\n",
        "[-0.21557732  1.08753659 -0.81284445  0.65051999 -0.48500761 -0.71379641\n",
        "  0.93627376 -0.91966905 -1.44485769 -0.79344869 -0.59309978  0.13385896\n",
        " -1.3633538  -0.06926551  0.28891933]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def distrib_under_H0(x, stat, nperm = 1000):\n",
      "    dist = []\n",
      "    for i in xrange(nperm):\n",
      "        psign = np.random.randint(2, size=len(x))\n",
      "        psign[psign==0] = -1\n",
      "        dist.append(stat(x * psign))\n",
      "    return dist\n",
      "\n",
      "def my_stat(x):\n",
      "    return(mean(x))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d_H0 = distrib_under_H0(z, my_stat)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size(d_H0)\n",
      "hist(d_H0)\n",
      "mean(z)\n",
      "empirical_p = float( sum(d_H0 > my_stat(z)))/size(d_H0)\n",
      "print my_stat(z), empirical_p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.40320746564 0.028\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD9CAYAAACyYrxEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFlJJREFUeJzt3X9M1Pfhx/HXR+GPJkCVKqc9XHCKgRMQVlS6xeaIYteu\nJTgcK7aV1R9Z3JbF1vRXlq3QZIJpTIYa0/5hW9zcrNlWMJsa22XXr9N2uBbbpSzKmrOFE4gWmFjS\nAvr+/uG8+QPwOO4OfPt8JCZ497nP++XJ53Xwufe9P44xxggAYK1J4x0AABBdFD0AWI6iBwDLUfQA\nYDmKHgAsR9EDgOVGLPrW1lYVFhZq/vz5ysrK0rZt2yRJlZWVSk1NVV5envLy8nTw4MHgY6qrq5We\nnq6MjAwdPnw4uukBADfljDSPvqOjQx0dHcrNzdWFCxd0zz33qL6+Xvv27VNiYqKeeuqpa7Zvbm7W\nqlWrdPz4cQUCAS1btkynTp3SpEn84gAA42XEBp4xY4Zyc3MlSQkJCcrMzFQgEJAkDfX60NDQoPLy\ncsXHxystLU1z585VY2NjFGIDAEIVF+qGp0+fVlNTkwoKCnT06FFt375du3fvVn5+vrZu3aopU6bo\nzJkzKigoCD4mNTU1+MJwheM4kUsPALeRcBcyCOmcyoULF7Ry5UrV1tYqISFBGzZskN/v14kTJzRz\n5kxt2rRp2McOVezGmAn/54UXXhj3DLbkvBUykpOcE/3PWNy06AcGBlRaWqrHHntMJSUlkqSUlBQ5\njiPHcbRu3brg6Rm3263W1tbgY9va2uR2u8cUEAAwNiMWvTFGa9eulcfj0caNG4O3t7e3B79+8803\nlZ2dLUkqLi7W3r171d/fL7/fr5aWFi1atChK0QEAoRjxHP3Ro0f1m9/8Rjk5OcrLy5Mkbd68Wb/7\n3e904sQJOY6j2bNn65VXXpEkeTwelZWVyePxKC4uTjt37rxlz8l7vd7xjhCSWyHnrZBRImekkXPi\nGHF6ZVQGdJwxn28CgNvNWLqTCe4AYDmKHgAsR9EDgOUoegCwHEUPAJaj6AHAchQ9AFiOogcAy1H0\nAGA5ih4ALEfRA4DlKHoAsBxFDwCWo+gBwHIUPQBYjqIHAMtR9ABgOYoeACxH0QOA5Sh6ALAcRQ8A\nlqPoAcByFD0AWI6iBwDLUfS4rSUlJctxnKj/SUpKHu9/Km5jjjHGxHRAx1GMhwSG5TiOpFh8P/J9\nj7EZS3fyEz0AWI6iBwDLUfQAYDmKHgAsR9EDgOUoegCwHEUPAJaj6AHAchQ9AFhuxKJvbW1VYWGh\n5s+fr6ysLG3btk2S1NXVpaKiIs2bN0/Lly9XT09P8DHV1dVKT09XRkaGDh8+HN30AICbGnEJhI6O\nDnV0dCg3N1cXLlzQPffco/r6er322muaNm2annnmGW3ZskXd3d2qqalRc3OzVq1apePHjysQCGjZ\nsmU6deqUJk363+sJSyBgImEJBNwqorYEwowZM5SbmytJSkhIUGZmpgKBgPbv36+KigpJUkVFherr\n6yVJDQ0NKi8vV3x8vNLS0jR37lw1NjaGFQwAEBlxoW54+vRpNTU1afHixers7JTL5ZIkuVwudXZ2\nSpLOnDmjgoKC4GNSU1MVCARu2FdlZWXwa6/XK6/XG2Z8ALCTz+eTz+eLyL5CKvoLFy6otLRUtbW1\nSkxMvOa+K8uwDmeo+64uegDAja7/Ibiqqirsfd101s3AwIBKS0v1+OOPq6SkRNLln+I7OjokSe3t\n7UpJSZEkud1utba2Bh/b1tYmt9sddjgAwNiNWPTGGK1du1Yej0cbN24M3l5cXKy6ujpJUl1dXfAF\noLi4WHv37lV/f7/8fr9aWlq0aNGiKMaHrWJ1QRDgdjDirJu//e1vuu+++5STkxM8KKqrq7Vo0SKV\nlZXps88+U1pamvbt26cpU6ZIkjZv3qxXX31VcXFxqq2t1f3333/tgMy6QQhiORuGWTe4FYylO7nC\nFCYkih64FleYAgAMi6IHAMtR9ABgOYoeACxH0QOA5Sh6ALAcRQ8AlqPoAcByFD0AWI6iBwDLUfQA\nYDmKHgAsR9EDgOUoegCwXMjXjAUwFnExudBJYuJUnT/fFfVxcGthPXpMSDauR8+69xgL1qMHAAyL\nogcAy1H0AGA5ih4ALEfRA4DlmF6JUUlKSlZvb/d4xwAwCkyvxKgw7XHij8PxZSemVwIAhkXRA4Dl\nKHoAsBxFDwCWo+gBwHIUPQBYjqIHAMtR9ABgOYoeACxH0QOA5Sh6ALAcRQ8AlqPoAcByIxb9mjVr\n5HK5lJ2dHbytsrJSqampysvLU15eng4ePBi8r7q6Wunp6crIyNDhw4ejlxoAELIRlyk+cuSIEhIS\ntHr1av3zn/+UJFVVVSkxMVFPPfXUNds2Nzdr1apVOn78uAKBgJYtW6ZTp05p0qRrX0tYpvjWxjLF\nE38cji87RW2Z4iVLlmjq1Kk33D7UYA0NDSovL1d8fLzS0tI0d+5cNTY2hhUKABA5YV1havv27dq9\ne7fy8/O1detWTZkyRWfOnFFBQUFwm9TUVAUCgSEfX1lZGfza6/XK6/WGEwMArOXz+eTz+SKyr1EX\n/YYNG/SLX/xCkvTzn/9cmzZt0q5du4bc9vKv+Te6uugBADe6/ofgqqqqsPc16lk3KSkpchxHjuNo\n3bp1wdMzbrdbra2twe3a2trkdrvDDgYAiIxRF317e3vw6zfffDM4I6e4uFh79+5Vf3+//H6/Wlpa\ntGjRosglBQCEZcRTN+Xl5XrnnXd07tw5zZo1S1VVVfL5fDpx4oQcx9Hs2bP1yiuvSJI8Ho/Kysrk\n8XgUFxennTt3DnvqBgAQOyNOr4zKgEyvvKUxvXLij8PxZaeoTa8EANz6KHoAsBxFDwCWo+gBwHIU\nPQBYjqIHAMtR9ABgOYoeACxH0QOA5Sh6ALAcRQ8AlqPoAcByFD0AWI6iBwDLUfQAYDmKHgAsR9ED\ngOUoegCwHEUPAJaj6AHAchQ9AFiOogcAy1H0AGA5ih4ALEfRA4DlKHoAsBxFDwCWo+gBwHIUPQBY\njqIHAMtR9ABgOYoeACxH0QOA5Sh6ALAcRQ8AlqPoAcByIxb9mjVr5HK5lJ2dHbytq6tLRUVFmjdv\nnpYvX66enp7gfdXV1UpPT1dGRoYOHz4cvdQAgJCNWPRPPPGEDh06dM1tNTU1Kioq0qlTp7R06VLV\n1NRIkpqbm/XGG2+oublZhw4d0o9+9CNdunQpeskBACEZseiXLFmiqVOnXnPb/v37VVFRIUmqqKhQ\nfX29JKmhoUHl5eWKj49XWlqa5s6dq8bGxijFBgCEKm60D+js7JTL5ZIkuVwudXZ2SpLOnDmjgoKC\n4HapqakKBAJD7qOysjL4tdfrldfrHW0MALCaz+eTz+eLyL5GXfRXcxxHjuOMeP9Qri56AMCNrv8h\nuKqqKux9jXrWjcvlUkdHhySpvb1dKSkpkiS3263W1tbgdm1tbXK73WEHAwBExqiLvri4WHV1dZKk\nuro6lZSUBG/fu3ev+vv75ff71dLSokWLFkU2LQBg1EY8dVNeXq533nlH586d06xZs/Tiiy/queee\nU1lZmXbt2qW0tDTt27dPkuTxeFRWViaPx6O4uDjt3LlzxNM6AIDYcIwxJqYDOo5iPCQi6PKLdyz+\n/xgn3HE4vuw0lu7kk7EAYDmKHgAsR9EDgOXGNI8ewEQTF5NJEImJU3X+fFfUx0Fk8GYsRoU3Yxnn\nyjgcx7HFm7EAgGFR9ABgOYoeACxH0QOA5Sh6ALAcRQ8AlqPoAcByFD0AWI6iBwDLUfQAYDmKHgAs\nx6JmlkhKSlZvb/d4xwAwAbGomSVYbIxxYj0Ox3FssagZAGBYFD0AWI6iBwDLUfQAYDmKHgAsR9ED\ngOUoegCwHEUPAJaj6AHAchQ9AFiOogcAy1H0AGA5ih4ALEfRA4DlKHoAsBxFDwCWo+gBwHJhX0ow\nLS1NSUlJmjx5suLj49XY2Kiuri59//vf16effqq0tDTt27dPU6ZMiWReAMAohf0TveM48vl8ampq\nUmNjoySppqZGRUVFOnXqlJYuXaqampqIBQUAhGdMp26uv37h/v37VVFRIUmqqKhQfX39WHYPAIiA\nsE/dOI6jZcuWafLkyfrhD3+o9evXq7OzUy6XS5LkcrnU2dk55GMrKyuDX3u9Xnm93nBjAICVfD6f\nfD5fRPblmDAvK97e3q6ZM2fq7NmzKioq0vbt21VcXKzu7u7gNsnJyerq6rp2wDFcyRzDcxxHUiye\nV8ZhnMvjcBzH1li6M+xTNzNnzpQkTZ8+XStWrFBjY6NcLpc6OjokXX4hSElJCXf3AIAICavo+/r6\n1NvbK0n64osvdPjwYWVnZ6u4uFh1dXWSpLq6OpWUlEQuKQAgLGGduvH7/VqxYoUkaXBwUI8++qie\nf/55dXV1qaysTJ999tmw0ys5dRMdnLphnFiPw3EcW2PpzrDP0YeLoo8Oip5xYj0Ox3Fsjcs5egDA\nrYGiBwDLUfQAYDmKHgAsR9EDgOUoegCwHEUPAJaj6AHAchQ9AFgu7GWKAdzO4v77aezoSkycqvPn\nu26+IUZE0QMIw6BisdRCb2/0X0xuB5y6AQDLUfQAYDmKHgAsR9EDgOUoegCwHLNuoiwpKVm9vd03\n3xAAooQrTEUZV35iHMYZ2zi3U1+MhCtMAQCGRdEDgOUoegCwHEUPAJaj6AHAchQ9AFiOogcAy1H0\nAGA5ih4ALMcSCAAmMK5kFQkUPYAJjCtZRQKnbgDAchQ9AFiOogcAy1H0AGA5ih4ALEfRA4DlKPph\n+Hy+8Y4QIt94BwiBb7wDhMg33gFC5BvvACHyjXeAEPnGO0DURXwe/aFDh7Rx40ZdvHhR69at07PP\nPhvpISKip6dH9fX1w16aq6GhQX6/P8apwuGT5B3nDDfj08TPKJEz0ny6dXLaLaJFf/HiRf3kJz/R\n22+/LbfbrYULF6q4uFiZmZmRHCYi3nrrLW3YUKnJkwuHvL+//1O9/fb/xTgVAEReRIu+sbFRc+fO\nVVpamiTpkUceUUNDw4QsekmKj1+o3t7Xhrm3UgMDlREY5fUI7AMAwhfRog8EApo1a1bw76mpqfr7\n3/9+w3axWLsidCNlqYrBGJFwJWesntdwxgnnuRyPf0+k/s9vNs5YjZRzIn0fROL5jM2/Z2L1UmRF\ntOhDeaKGOycOAIiOiM66cbvdam1tDf69tbVVqampkRwCADBKES36/Px8tbS06PTp0+rv79cbb7yh\n4uLiSA4BABiliJ66iYuL044dO3T//ffr4sWLWrt27YR9IxYAbhcR/8DUAw88oJMnT+rf//63nn/+\neXV1damoqEjz5s3T8uXL1dPTM+TjqqurNX/+fGVnZ2vVqlX66quvIh1tRKHm7Onp0cqVK5WZmSmP\nx6P33ntvQuaULk93zcvL08MPPxzDhJeFkrO1tVWFhYWaP3++srKytG3btphkO3TokDIyMpSenq4t\nW7YMuc1Pf/pTpaena8GCBWpqaopJruvdLOeePXu0YMEC5eTk6Fvf+pY++uijcUgZ2vMpScePH1dc\nXJz++Mc/xjDd/4SS0+fzKS8vT1lZWfJ6vbEN+F83y3nu3Dl9+9vfVm5urrKysvT666/ffKcmyp5+\n+mmzZcsWY4wxNTU15tlnn71hG7/fb2bPnm2+/PJLY4wxZWVl5vXXX492tFHnNMaY1atXm127dhlj\njBkYGDA9PT0xy2hM6DmNMWbr1q1m1apV5uGHH45VvKBQcra3t5umpiZjjDG9vb1m3rx5prm5Oaq5\nBgcHzZw5c4zf7zf9/f1mwYIFN4z55z//2TzwwAPGGGPee+89s3jx4qhmCjfnsWPHgt9/Bw8enLA5\nr2xXWFhovvOd75jf//73EzJnd3e38Xg8prW11RhjzNmzZydkzhdeeME899xzwYzJyclmYGBgxP1G\nfQmE/fv3q6KiQpJUUVGh+vr6G7ZJSkpSfHy8+vr6NDg4qL6+Prnd7mhHG3XO//znPzpy5IjWrFkj\n6fKpqjvvvHPC5ZSktrY2HThwQOvWrRuXmU6h5JwxY4Zyc3MlSQkJCcrMzNSZM2eimuvqz3rEx8cH\nP+sxXPbFixerp6dHnZ2dUc0VTs577703+P23ePFitbW1xTRjqDklafv27Vq5cqWmT58e84xSaDl/\n+9vfqrS0NDiBZNq0aRMy58yZM3X+/HlJ0vnz53XXXXcpLm7ks/BRL/rOzk65XC5JksvlGvKASU5O\n1qZNm/S1r31Nd999t6ZMmaJly5ZFO9qoc/r9fk2fPl1PPPGEvvGNb2j9+vXq6+ubcDkl6cknn9RL\nL72kSZPGZzmjUHNecfr0aTU1NWnx4sVRzTXUZz0CgcBNt4l1iYaS82q7du3Sgw8+GIto1wj1+Wxo\naNCGDRskjc989VBytrS0qKurS4WFhcrPz9evf/3rWMcMKef69ev18ccf6+6779aCBQtUW1t70/1G\n5M3YoqIidXR03HD7L3/5y2v+7jjOkP/Jn3zyiX71q1/p9OnTuvPOO/W9731Pe/bs0aOPPhqJeBHL\nOTg4qA8++EA7duzQwoULtXHjRtXU1OjFF1+cUDn/9Kc/KSUlRXl5eVFdnG2sOa+4cOGCVq5cqdra\nWiUkJEQ85/VZQnH9b0GxLqfRjPfXv/5Vr776qo4ePRrFREMLJeeV48RxHBljxuU3zFByDgwM6IMP\nPtBf/vIX9fX16d5771VBQYHS09NjkPCyUHJu3rxZubm58vl8+uSTT1RUVKQPP/xQiYmJwz4mIkX/\n1ltvDXufy+VSR0eHZsyYofb2dqWkpNywzT/+8Q9985vf1F133SVJ+u53v6tjx45FvOjHmjM1NVWp\nqalauHChJGnlypWqqamJaMZI5Dx27Jj279+vAwcO6Msvv9T58+e1evVq7d69e0LllC4fXKWlpXrs\nscdUUlIS0XxDCeWzHtdv09bWFvNTiaF+JuWjjz7S+vXrdejQIU2dOjWWESWFlvP999/XI488Iuny\nG4kHDx5UfHx8TKdeh5Jz1qxZmjZtmu644w7dcccduu+++/Thhx/GtOhDyXns2DH97Gc/kyTNmTNH\ns2fP1smTJ5Wfnz/8jqPyjsJVnn76aVNTU2OMMaa6unrIN+VOnDhh5s+fb/r6+sylS5fM6tWrzY4d\nO6IdbdQ5jTFmyZIl5uTJk8aYy2+KPPPMMzHLaEzoOa/w+XzmoYceikW0a4SS89KlS+bxxx83Gzdu\njFmugYEB8/Wvf934/X7z1Vdf3fTN2HfffXdc3uQMJeenn35q5syZY959992Y57silJxX+8EPfmD+\n8Ic/xDDhZaHk/Ne//mWWLl1qBgcHzRdffGGysrLMxx9/POFyPvnkk6aystIYY0xHR4dxu93m888/\nH3G/US/6zz//3CxdutSkp6eboqIi093dbYwxJhAImAcffDC43ZYtW4zH4zFZWVlm9erVpr+/P9rR\nwsp54sQJk5+fb3JycsyKFStiPusm1JxX+Hy+cZl1E0rOI0eOGMdxzIIFC0xubq7Jzc01Bw8ejHq2\nAwcOmHnz5pk5c+aYzZs3G2OMefnll83LL78c3ObHP/6xmTNnjsnJyTHvv/9+1DOFk3Pt2rUmOTk5\n+NwtXLhwQua82ngVvTGh5XzppZeCPVRbWzshc549e9Y89NBDJicnx2RlZZk9e/bcdJ+OMSw+AwA2\n4wpTAGA5ih4ALEfRA4DlKHoAsBxFDwCWo+gBwHL/D72ii/VLOiE8AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3da6250>"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this instance, we have a significant p value at the 5 percent risk of error because under the null hypothesis, the actual data statistics is observed less than 5 percent of the times. This framework is \n",
      "\n",
      "* powerful : any statistics can be tested\n",
      "* tricky \n",
      "    - non exangeable data  \n",
      "    - what exactly is the null hypothesis tested\n",
      "* demanding in computing time : make sure we can get a good estimation of the probability"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "How do we deal with multiple comparison using permutations ?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We simply take the distribution of the maximum statistics across tests !\n",
      "\n",
      "An example ?"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "FWER : RFT"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    }
   ],
   "metadata": {}
  }
 ]
}