{
 "metadata": {
  "name": "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",
      "\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": [
        "Nubmer of test significant at the 0.050000 level: 5 \n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "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: %d, corrected: %d  5 0\n",
        "Uncorrected threshold: 0.050000, corrected threshold: 0.000513\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "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": 4
    },
    {
     "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": 5
    },
    {
     "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": 6,
       "text": [
        "[<matplotlib.lines.Line2D at 0x3afdbd0>]"
       ]
      },
      {
       "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 0x399b550>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "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": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEMCAYAAADEXsFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1YVGXeB/DvwIw+KCgipjJDosIyY8oAgW+BjrWGscn6\nxoq7m2Uo5OambrX51B/B7mbR9nrJVuSqZQbi24qVjC7mFKzKmBj2hBeiiY1UlgmC0QqM5/ljYGCE\ncWbgzAvw/VwXl3Nm7rnPb45wvnPe7iMRBEEAERH1e17uLoCIiDwDA4GIiAAwEIiIqBUDgYiIADAQ\niIioFQOBiIgAOCEQHn74YYwcORKTJk2y2uaxxx5DWFgY1Go1Tp48KXYJRETUDaIHwrJly6DVaq2+\nvn//fpw9exZVVVV4++23sXLlSrFLICKibhA9EOLj4zFs2DCrr+/btw8PPvggAGDKlCmoq6vDpUuX\nxC6DiIgcJHX1DGtqahAcHGyeVigUuHjxIkaOHGnRTiKRuLo0IqI+obsDULjloPLNxVpb+QuCwB9B\nwLPPPuv2Gjzlh8uCy6InywLoDz/d5/JAkMvlMBgM5umLFy9CLpe7ugwSgVIJSCTtP15eltPO+MnM\ndP482n5kMvH68vbuWX8DBji2LKRS1y2n7v6I+fti7+9FX/fSSz17v8sDISkpCVu3bgUAHDt2DP7+\n/p12F1Hv8N13ltNCHxsmsaVFvL5u3OhZf83NjrU3Grs/L1fxxN+Xp582/RsZCdx7LzBihClcZbL2\nfwHTv7Gx7e+bMKE9+NsMGAC88Ybp8RtvAP/zP0B8vKnPjuE0axawc6dlHW3za5vnwIFASYmpJmnr\njn4vL8v5FRYCjz/ewwUgiCwlJUUYPXq0IJPJBIVCIWzatEl46623hLfeesvc5tFHHxXGjx8vRERE\nCCdOnOiyHyeU1msdPnzY3SV0smKFIEilgmD6s3blz2GXzUsiccfn88xl4fk/9i2LUaPaHw8YIAjv\nvCMIXl6CMGKEIFRXu/uvShw9WXd2/51OxkDwXOHhguDtbfmHtnOn6Y8tMND0hzZ1qunfmx9Pm2b5\n3qgoQZg1SxBkMtP0kCGmdoGBgjBypCDExQnCwIGmP9i4ONNrXQWRv78g3Hab6f1tzw0Z0r5Sl0gE\nYdgwUz8jRpj6njat/TV/f0EYOrT9vcOHC8KpU6ba2lYYJSWmPgBBiIw0za/jZ5HJTP20TXt5CUJA\ngCCMHm16b1CQIMTHmz7PyJGCUFho+VkCAkyfsW1ZtH2me+8VhMREQYiNNT3v62v52YcMaa9j+PD2\neU2davpsUqmpr7ZlERfX/tk7/shkpuUbF9f+Ods+B2Dqx9vbss+2+mUy07KaOtU07ednauftbXqt\nrb+hQ039xcWZPldgYHsfISGmfydMaO9HKjW99+b/847/z1KpqV+JxLQMAwJMnyMgwPL/Jy5OEGpr\n3f3X43w9WXdKWjvwOBKJBB5aWr+lVAJVVabdH21kMuD77wF/f/fVRUTterLu5NAVZJe0NODsWcsw\nAIATJxgGRH0FA4HscuZM5wOVJSXALUYoIaJexuUXplHvolQClZXt097ewPDhgF4PjBnjvrqISHw8\nhkC35O8PXL3aPu3l1TtOaSTqr3gMgUSRlgYMGmR5IU/HMACATz91T21E5HzcZdTPDRhg30VPXl6m\nMLjrLufXRETuwS2Efsrf37SStycMSkpMu4kYBkR9G7cQ+hl7tgi8vEynl0okQHExg4Cov+AWQj9j\nLQzaxkSZNQv48UfTtZ03bjAMiPoTbiH0E0pl58HoACAuDvjgA15cRkQMhF6pqyEkHDViBHD8OK8l\nIKJ23GXUC333Xc/CIDHRNP4Qw4CIOmIg9EJtY7J3x5gxwPvvi1cLEfUdDIRexN/fdHOMujogIMA0\nhMSAAe03zAAsb+LR1r7t9bg44PPPebyAiLrGoSt6Eam0fdgIHx+gsdG99RCR5+HQFX1c2/1028JA\nIgFKS91bExH1PTzLyIPYO4zE7NkcdpqIxMddRh6k4423rVEqgaNHeRyAiLrmcbuMtFotlEolwsLC\nkJWV1en12tpazJ8/H2q1GlOmTMGXX37pjDJ6BX//9pFFbYmKYhgQkfOIHghGoxGrVq2CVqtFRUUF\n8vLycPr0aYs269evR3R0NMrLy7F161asXr1a7DJ6BaWy8/DSQPvZQRKJaVwhLy/TAHNlZQwDInIe\n0Y8h6PV6hIaGIiQkBACQkpKCgoICqFQqc5vTp09j3bp1AIDw8HBUV1fjhx9+wIgRIyz6ysjIMD/W\naDTQaDRil+sWaWnA5s1d32impITjBxGR/XQ6HXQ6nSh9iR4INTU1CA4ONk8rFAqU3nRKjFqtxp49\nexAXFwe9Xo8LFy7g4sWLtwyEvuJWB44ZBkTkqJu/LGdmZna7L9F3GUns2Bm+bt061NXVISoqCtnZ\n2YiKioJ323CbfVhaWtdhUFhoGl2UYUBE7iT6FoJcLofBYDBPGwwGKBQKizZ+fn7YvHmzeXrs2LEY\nN26c2KV4nDNnOj/HrQIi8hSibyHExMSgqqoK1dXVaGpqQn5+PpKSkizaXL16FU1NTQCAjRs3YubM\nmfD19RW7FI8zaFD747g4oLaWYUBEnkP0LQSpVIrs7GwkJCTAaDQiNTUVKpUKOTk5AID09HRUVFTg\noYcegkQiwcSJE7Fp0yaxy/A4aWlAfT0wahRw7BhHGiUiz8ML01xEowE++cT0ODkZ2LHDreUQUR/l\ncRemkaW0NOCLL0yPo6KAt992bz1ERF3hWEZOdvNppiNG8OIyIvJM3EJwsptPM/2//3NPHUREtjAQ\nRJSWZjqTqG1sopsvyRg4EDhyxD21ERHZwkAQ0ZkzwM8/d/3agAFAZSXPLiIiz8VAEFHH6ww6GjMG\nuHSJYUBEno2BIKLcXOBXvwICA9vvYxwTw/sYE1HvwLOMRJKWZtplNGgQUFXFACCi3ocXponAy8s0\nOF0bXnhGRO7Sk3UnA0EEN59NVFvLLQQicg9eqexBdu5kGBBR78RAEEFhYfu/ixa5txYiou7iLiMi\noj6Eu4zcJC3NNIppYiJQV+fuaoiIeoaB0ANnzpiGtC4sNIUDEVFv5tG7jFasELBtm/XhIDzFHXeY\nboXJg8lE5G59dpfRrcYG8iTnzjEMiKj38+hAsDY2kCeRSAC93t1VEBH1nFMCQavVQqlUIiwsDFlZ\nWZ1ev3z5MubMmYPIyEhMnDgR77zzTpf95OaaDtjedhsgk3X9LXzYMNOw0rGxppWzt7fp+aFD28cT\nAoDIyPYLyLxaP7Wvr2VfkZGm14YNM93IJiDANF9rZDKgvByYNMl6GyKi3kL0YwhGoxHh4eEoKiqC\nXC5HbGws8vLyoFKpzG0yMjJw/fp1PP/887h8+TLCw8Nx6dIlSDuswXnaKRGR4zzqGIJer0doaChC\nQkIgk8mQkpKCgoICizajR49GfX09AKC+vh7Dhw+3CAMiInI90dfCNTU1CA4ONk8rFAqUlpZatFmx\nYgXuvvtuBAUFoaGhATusjASXkZFhfqzRaKDRaMQul4ioV9PpdNDpdKL0JXogSG4e6a0L69evR2Rk\nJHQ6Hc6dO4fZs2ejvLwcfn5+Fu06BgIREXV285flzMzMbvcl+i4juVwOg8FgnjYYDFAoFBZtjhw5\nguTkZADA+PHjMXbsWFRWVopdChEROUD0QIiJiUFVVRWqq6vR1NSE/Px8JCUlWbRRKpUoKioCAFy6\ndAmVlZUYN26c2KUQEZEDRN9lJJVKkZ2djYSEBBiNRqSmpkKlUiEnJwcAkJ6ejqeffhrLli2DWq3G\njRs38OKLLyIgIEDsUoiIyAEePXSFh5ZGROSxPOq0UyIi6p0YCEREBICBQERErRgIREQEgIFARESt\nGAhERASAgUBERK0YCEREBICBQERErRgIREQEgIFAREStPDoQ/P1N9za+cMHdlRAR9X0ePbgdYCpN\noQA63GKBiIis6NOD2w0aBJSUuLsKIqK+z6MDQaEAKiqAMWPcXQkRUd/n0buMPLQ0IiKP5bRdRkaj\nEa+++mq3OiYiot7lloHg7e2N3NxcV9VCRERuZHOX0dq1a9Hc3IzFixdj8ODB5uejo6OdWxh3GRER\nOawn606bgaDRaFpPAbV0+PBhq+/RarVYs2YNjEYjli9fjqeeesri9Zdeegnvv/8+AKClpQWnT5/G\n5cuX4e/v314YA4GIyGFODQRHGY1GhIeHo6ioCHK5HLGxscjLy4NKpeqy/YcffojXXnsNRUVFloUx\nEIiIHNaTdafU2gsvv/yyufOu/OlPf+ryeb1ej9DQUISEhAAAUlJSUFBQYDUQcnNzsWTJEkdqJiIi\nJ7AaCA0NDZBIJKisrMTx48eRlJQEQRDw4YcfYvLkyVY7rKmpQXBwsHlaoVCgtLS0y7aNjY04cOAA\n3njjjS5fz8jIMD/WaDTQaDQ2Pg4RUf+i0+mg0+lE6ctqILStjOPj41FWVgY/Pz8AQGZmJhITE612\naG2LoisffPAB4uLiLI4ddFUDERF17eYvy5mZmd3uy+aVyt9//z1kMpl5WiaT4fvvv7faXi6Xw9Bh\n4CGDwQCFQtFl2+3bt3N3ERGRh7C6hdBm6dKlmDx5MhYsWABBELB37148+OCDVtvHxMSgqqoK1dXV\nCAoKQn5+PvLy8jq1u3r1Kj799FNe50BE5CGsBsJXX32FcePG4ZlnnsGcOXNQXFwMiUSCd955B1FR\nUdY7lEqRnZ2NhIQEGI1GpKamQqVSIScnBwCQnp4OANi7dy8SEhLg4+Mj8kciIqLusHra6Z133okT\nJ07gnnvuwaFDh1xdF087JSLqBqecdmo0GvHcc8+hsrISr7zyisUMJBKJ1dNOiYiod7J6UHn79u3w\n9vaG0WhEQ0MDrl27Zv5paGhwZY1EROQCNq9U3r9//y1PM3UW7jIiInKcRw1dIRYGAhGR4/r0LTSJ\niMg1GAhERATgFmcZ7d6927zp0dVwFAsWLHBqYURE5FpWA+GDDz6ARCLB999/jyNHjuDuu+8GYLoP\nwvTp0xkIRER9jNVAeOeddwAAs2fPRkVFBUaPHg0A+Pbbb285dAUREfVONo8hGAwGjBo1yjw9cuRI\nfP31104tioiIXM/m4Ha//OUvkZCQgN/+9rcQBAH5+fmYPXu2K2ojIiIXsnkdgiAI+Ne//oXi4mIA\nwIwZMzB//nznF8brEIiIHOaUsYw6dh4dHQ0/Pz/Mnj0bjY2NaGhoMN8wh4iI+gabxxDefvttJCcn\n45FHHgEAXLx4EfPmzXN6YURE5Fo2A+Ef//gHSkpKMGTIEADAL37xi1veMY2IiHonm4EwcOBADBw4\n0Dzd0tLi0H2TiYiod7AZCDNnzsRzzz2HxsZG/Pvf/0ZycjLmzp3ritqIiMiFbJ5ldOPGDfzzn//E\nwYMHAQAJCQlYvny507cSeJYREZHjnDr89euvv47Vq1fbfE5sDAQiIsc5dfjrtiEsOtqyZcst36PV\naqFUKhEWFoasrKwu2+h0OkRFRWHixInQaDR2FUtERM5jdQshLy8Pubm5KC4uRnx8vPn5hoYGeHt7\n49ChQ112aDQaER4ejqKiIsjlcsTGxiIvLw8qlcrcpq6uDnfddRcOHDgAhUKBy5cvIzAw0LIwbiEQ\nETnMKRemTZ8+HaNHj8YPP/yAJ554wjwDPz8/qNVqqx3q9XqEhoYiJCQEAJCSkoKCggKLQMjNzcXC\nhQuhUCgAoFMYEBGR61kNhDFjxmDMmDE4duyYQx3W1NQgODjYPK1QKFBaWmrRpqqqCs3NzZg1axYa\nGhqwevVqPPDAA536ysjIMD/WaDTctUREdBOdTgedTidKXzaHrjh69Cgee+wxnD59GtevX4fRaISv\nry/q6+u7bG/P2UfNzc0oKyvDoUOH0NjYiGnTpmHq1KkICwuzaNcxEIiIqLObvyxnZmZ2uy+bgbBq\n1Sps374dv/nNb/DZZ59h69atqKystNpeLpfDYDCYpw0Gg3nXUJvg4GAEBgbCx8cHPj4+mDFjBsrL\nyzsFAhERuY5d91QOCwuD0WiEt7c3li1bBq1Wa7VtTEwMqqqqUF1djaamJuTn5yMpKcmiza9//WuU\nlJTAaDSisbERpaWlmDBhQs8+CRER9YjNLYTBgwfj+vXrUKvV+POf/4xRo0bd8gi2VCpFdnY2EhIS\nYDQakZqaCpVKhZycHABAeno6lEol5syZg4iICHh5eWHFihUMBCIiN7N5YVp1dTVGjhyJpqYmvPrq\nq6ivr8cf/vAHhIaGOrcwnnZKROQwp16p7C4MBCIixznlOoTk5GTs3LkTkyZN6nKGp06d6tYMiYjI\nM1ndQvjmm28QFBSE6urqLt/YduGZs3ALgYjIcU7ZQggKCgJguqfyqFGj4OPjAwD4+eefcenSpW7N\njIiIPJfN004XLVoEb2/v9jd4eWHRokVOLYqIiFzPZiAYjUYMGDDAPD1w4EA0Nzc7tSgiInI9m4EQ\nGBiIgoIC83RBQQEHoyMi6oNsnnZ69uxZ/O53v8M333wDwDRY3XvvvcfrEIiIPJBLrkO4du0aAMDX\n17dbM3IUA4GIyHFOvWNaXV0d1q5di5kzZ2LmzJl4/PHHcfXq1W7NzFEaDZCYCNTVuWR2RET9ms1A\nePjhhzFkyBDs3LkTO3bsgJ+fH5YtW+aK2vDJJ0BhIZCW5pLZERH1azZ3GanVapSXl9t8TvTCJBIA\nAmJjgYMHAX9/p86OiKhPcOouIx8fHxQXF5unS0pKMGjQoG7NzFHJyQwDIiJXsbmF8Pnnn2Pp0qXm\n4wbDhg3Du+++e8v7KotSGA8qExE5zCVnGV29ehUSiQRDhgzp1owcxUAgInKcU3cZvfbaa6ivr8eQ\nIUOwdu1aREdH48CBA92aGREReS6bgbB582YMGTIEBw8exJUrV7B161asW7fOFbUREZEL2QyEtk2P\njz76CA888AAmTpzo9KKIiMj1bAbCnXfeiXvvvRf79+9HQkIC6uvr4eVl821ERNTL2LXL6Pnnn8dn\nn32GwYMHo7m5GVu2bLnle7RaLZRKJcLCwpCVldXpdZ1Oh6FDhyIqKgpRUVH429/+1v1PQEREorB6\ng5w2R48ehVqthq+vL9577z2UlZVhzZo1VtsbjUasWrUKRUVFkMvliI2NRVJSElQqlUW7mTNnYt++\nfT3/BEREJAqbWwiPPPIIBg8ejPLycrzyyisIDQ3F0qVLrbbX6/UIDQ1FSEgIZDIZUlJSLIbPbsNT\nSomIPIvNLQSpVAqJRIK9e/fi0UcfxfLly7Fp0yar7WtqahAcHGyeVigUKC0ttWgjkUhw5MgRqNVq\nyOVyvPTSS5gwYUKnvjIyMsyPNRoNNBqNHR+JiKj/0Ol00Ol0ovRlMxD8/Pywfv16bNu2DcXFxTAa\njbe8Y5ppDKJbi46OhsFgwKBBg1BYWIh58+bhzJkzndp1DAQiIurs5i/LmZmZ3e7L5i6j/Px8DBw4\nEJs3b8aoUaNQU1ODJ5980mp7uVwOg8FgnjYYDFAoFBZt/Pz8zOMh3XfffWhubsaVK1e6+xmIiEgE\nNgNh9OjRePzxxxEfHw8AuHDhAo4dO2a1fUxMDKqqqlBdXY2mpibk5+cjKSnJos2lS5fMxxD0ej0E\nQUBAQEBPPgcREfWQzV1GAFBWVoa8vDzs2LEDY8eOxcKFC613KJUiOzsbCQkJMBqNSE1NhUqlQk5O\nDgAgPT0du3btwptvvgmpVIpBgwZh+/btXfaVmAjk5nK0UyIiV7A6uF1lZSXy8vKQn5+PESNGIDk5\nGX//+9/x9ddfu6aw1vshJCcDO3a4ZJZERL2eU0Y79fLywv3334/s7GzcfvvtAICxY8fi/Pnz3a/U\nkcIkEsTGCrwfAhGRA5wy2umePXvg4+ODGTNm4JFHHsGhQ4dcfu0Aw4CIyHVs3g/h2rVrKCgoQF5e\nHg4fPoylS5di/vz5uPfee51bGO+HQETkMJfcIAcArly5gl27dmH79u34+OOPuzVDezEQiIgc57JA\ncCUGAhGR45x6xzQiIuofGAhERASAgUBERK0YCEREBICBQERErRgIREQEgIFAREStGAhERASAgUBE\nRK0YCEREBICBQERErRgIREQEgIFAREStGAhERATASYGg1WqhVCoRFhaGrKwsq+2OHz8OqVSKPXv2\nOKMMIiJygOiBYDQasWrVKmi1WlRUVCAvLw+nT5/ust1TTz2FOXPm8L4HREQeQPRA0Ov1CA0NRUhI\nCGQyGVJSUlBQUNCp3YYNG7Bo0SKMGDHCal+JiUBdndgVEhFRV6Rid1hTU4Pg4GDztEKhQGlpaac2\nBQUF+Pjjj3H8+HFIJJIu+yoszMBddwHJyYBGo4FGoxG7XCKiXk2n00Gn04nSl+iBYG3l3tGaNWvw\nwgsvmG/1Zm2XUWxsBg4eBPz9xa6SiKhvuPnLcmZmZrf7Ej0Q5HI5DAaDedpgMEChUFi0OXHiBFJS\nUgAAly9fRmFhIWQyGZKSkizaMQyIiFxHIoh8RLelpQXh4eE4dOgQgoKCMHnyZOTl5UGlUnXZftmy\nZZg7dy4WLFhgWVgPbhRNRNRf9WTdKfoWglQqRXZ2NhISEmA0GpGamgqVSoWcnBwAQHp6utizJCIi\nEYi+hSAWbiEQETmuJ+tOXqlMREQAGAhERNSKgUBERAAYCERE1IqBQEREABgIRETUioFAREQAGAhE\nRNSKgUBERAAYCERE1IqBQEREABgIRETUioFAREQAGAhERNSKgUBERAAYCERE1IqBQEREABgIRETU\nioFAREQAnBQIWq0WSqUSYWFhyMrK6vR6QUEB1Go1oqKicOedd+Ljjz92RhlEROQAiSDyneyNRiPC\nw8NRVFQEuVyO2NhY5OXlQaVSmdv89NNPGDx4MADgiy++wPz583H27FnLwnpwo2giov6qJ+tO0bcQ\n9Ho9QkNDERISAplMhpSUFBQUFFi0aQsDALh27RoCAwPFLoOIiBwkFbvDmpoaBAcHm6cVCgVKS0s7\ntdu7dy/+93//F99++y0OHjzYZV8ZGRnmxxqNBhqNRuxyiYh6NZ1OB51OJ0pfou8y2r17N7RaLTZu\n3AgA2LZtG0pLS7Fhw4Yu2xcXF2P58uWorKy0LIy7jIiIHOZRu4zkcjkMBoN52mAwQKFQWG0fHx+P\nlpYW/Pjjj2KXQkREDhA9EGJiYlBVVYXq6mo0NTUhPz8fSUlJFm3OnTtnTrCysjIAwPDhw8UuhYiI\nHCD6MQSpVIrs7GwkJCTAaDQiNTUVKpUKOTk5AID09HTs3r0bW7duhUwmg6+vL7Zv3y52GURE5CDR\njyGIhccQiIgc51HHEIiIqHdiIBAREQAGAhERtWIgEBERAAYCERG1YiAQEREABgIREbViIBAREQAG\nAhERtfLoQEhMBOrq3F0FEVH/4NGBUFgIpKW5uwoiov7Bo8cyio0VcPAg4O/v7mqIiHqHnoxl5NGB\nUFsrMAyIiBzQZwPBQ0sjIvJYHO2UiIh6jIFAREQAGAhERNSKgUBERAAYCL2CTqdzdwkeg8uiHZdF\nOy4LcTglELRaLZRKJcLCwpCVldXp9ffffx9qtRoRERG46667cOrUKWeU0Wfwl70dl0U7Lot2XBbi\nkIrdodFoxKpVq1BUVAS5XI7Y2FgkJSVBpVKZ24wbNw6ffvophg4dCq1Wi7S0NBw7dkzsUoiIyAGi\nbyHo9XqEhoYiJCQEMpkMKSkpKCgosGgzbdo0DB06FAAwZcoUXLx4UewyiIjIQaJfmLZr1y4cOHAA\nGzduBABs27YNpaWl2LBhQ5ftX3rpJZw5cwZvv/22ZWESiZhlERH1G91drYu+y8iRFfnhw4exefNm\n/Oc//+n0Gq9SJiJyLdEDQS6Xw2AwmKcNBgMUCkWndqdOncKKFSug1WoxbNgwscsgIiIHiX4MISYm\nBlVVVaiurkZTUxPy8/ORlJRk0ebrr7/GggULsG3bNoSGhopdAhERdYPoWwhSqRTZ2dlISEiA0WhE\namoqVCoVcnJyAADp6en4y1/+gtraWqxcuRIAIJPJoNfrxS6FiIgcIbhZYWGhEB4eLoSGhgovvPBC\nl23++Mc/CqGhoUJERIRQVlbm4gpdx9ay2LZtmxARESFMmjRJmD59ulBeXu6GKl3Dnt8LQRAEvV4v\neHt7C7t373Zhda5lz7I4fPiwEBkZKdxxxx3CzJkzXVugC9laFj/88IOQkJAgqNVq4Y477hC2bNni\n+iJdYNmyZcJtt90mTJw40Wqb7qw33RoILS0twvjx44Xz588LTU1NglqtFioqKizafPTRR8J9990n\nCIIgHDt2TJgyZYo7SnU6e5bFkSNHhLq6OkEQTH8Y/XlZtLWbNWuW8Ktf/UrYtWuXGyp1PnuWRW1t\nrTBhwgTBYDAIgmBaKfZF9iyLZ599Vli3bp0gCKblEBAQIDQ3N7ujXKf69NNPhbKyMquB0N31pluH\nrrDnmoV9+/bhwQcfBGC6ZqGurg6XLl1yR7lOxes32tmzLABgw4YNWLRoEUaMGOGGKl3DnmWRm5uL\nhQsXmk/eCAwMdEepTmfPshg9ejTq6+sBAPX19Rg+fDikUtH3jLtdfHz8LU/G6e56062BUFNTg+Dg\nYPO0QqFATU2NzTZ9cUVoz7LoaNOmTUhMTHRFaS5n7+9FQUGB+ThUX71uxZ5lUVVVhStXrmDWrFmI\niYnBe++95+oyXcKeZbFixQp8+eWXCAoKglqtxuuvv+7qMj1Cd9ebbo1Oe/+IhZuuSeiLf/xiXb/R\nF9izLNasWYMXXnjBfHeom39H+gp7lkVzczPKyspw6NAhNDY2Ytq0aZg6dSrCwsJcUKHr2LMs1q9f\nj8jISOh0Opw7dw6zZ89GeXk5/Pz8XFChZ+nOetOtgWDPNQs3t7l48SLkcrnLanQVXr/Rzp5lceLE\nCaSkpAAALl++jMLCQshksk6nOPd29iyL4OBgBAYGwsfHBz4+PpgxYwbKy8v7XCDYsyyOHDmCZ555\nBgAwfvx4jB07FpWVlYiJiXFpre7W7fWmKEc4uqm5uVkYN26ccP78eeH69es2DyofPXq0zx5ItWdZ\nXLhwQRi8cy0jAAAE3UlEQVQ/frxw9OhRN1XpGvYsi44eeuihPnuWkT3L4vTp08I999wjtLS0CD/9\n9JMwceJE4csvv3RTxc5jz7JYu3atkJGRIQiCIHz33XeCXC4XfvzxR3eU63Tnz5+366CyI+tNt24h\n2HPNQmJiIvbv34/Q0FAMHjwYW7ZscWfJTsPrN9rZsyz6C3uWhVKpxJw5cxAREQEvLy+sWLECEyZM\ncHPl4rNnWTz99NNYtmwZ1Go1bty4gRdffBEBAQFurlx8S5YswSeffILLly8jODgYmZmZaG5uBtCz\n9abog9sREVHvxDumERERAAYCERG1YiAQEREABgIREbViIFCv5OvrK3qfFy5cQF5enuj93kpGRgZe\nfvlll86TyBoGAvVKzrha/fz588jNzRW931vpi1fdU+/FQKBeTafTQaPRIDk5GSqVCr///e/Nr4WE\nhOCpp55CREQEpkyZgnPnzgEAHnroIezevdvcrm1Yg3Xr1qG4uBhRUVGdxsDR6XSYMWMG7r//fiiV\nSqxcubLT0ABXr15FSEiIefqnn37C7bffjpaWFmzcuBGTJ09GZGQkFi1ahJ9//tncri0UNBoNTpw4\nAcB09fXYsWMBAEajEU8++SQmT54MtVrd6f7jRGJhIFCv9/nnn+P1119HRUUFvvrqKxw5cgSAaUXr\n7++PU6dOYdWqVVizZo35+a5kZWUhPj4eJ0+exOrVqzu9fvz4cWRnZ6OiogLnzp3Dnj17LF4fOnSo\neRwdAPjwww8xZ84cSKVSLFy4EHq9Hp9//jlUKhU2bdrUqX+JRNJlbZs2bYK/vz/0ej30ej02btyI\n6upqRxYRkV0YCNTrTZ48GUFBQZBIJIiMjLRYWS5ZsgQAkJKSgqNHj96yH1vXaE6ePBkhISHw8vLC\nkiVLUFJS0qnN4sWLkZ+fDwDYvn07Fi9eDAD44osvEB8fj4iICLz//vuoqKiw+/MdPHgQW7duRVRU\nFKZOnYorV67g7Nmzdr+fyF59b6Bw6ncGDhxofuzt7Y2WlpYu27V9+5ZKpbhx4wYA4MaNG2hqarJr\nPh2/vQuCAIlEgr179yIzMxOA6Zv83Llz8fTTT6O2thZlZWW4++67AZh2U+3btw+TJk3Cu+++a96K\n6KhjXf/9738tXsvOzsbs2bPtqpOou7iFQH1a27f1/Px8TJ8+HYDp2ELbvvp9+/aZx4Dx8/NDQ0OD\n1b70ej2qq6tx48YN7NixA/Hx8Zg3bx5OnjyJkydPIjo6Gr6+voiNjcVjjz2GuXPnmkPk2rVrGDVq\nFJqbm7Ft2zbz80KHobtDQkLw2WefAQB27dplnm9CQgLeeOMNc9CdOXMGjY2Noi0jojYMBOqVOn5b\nv9WZOrW1tVCr1diwYQNeffVVAKabqHzyySeIjIzEsWPHzKewqtVqeHt7IzIystNBZYlEgtjYWKxa\ntQoTJkzAuHHjMG/evC7nuXjxYuTm5pp3FwHAX//6V0yZMgVxcXFQqVQW/bbV/8QTT+DNN99EdHQ0\nfvzxR/Pzy5cvx4QJExAdHY1JkyZh5cqVVreCiHqCg9tRnzV27FicOHFClNEudTodXn75ZXzwwQci\nVEbkmbiFQH2WmOf4WzsDiKgv4RYCEREB4BYCERG1YiAQEREABgIREbViIBAREQAGAhERtWIgEBER\nAOD/Ad2okmpWJQqYAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3b07350>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alph = 0.05\n",
      "df = 14\n",
      "t14 = tdist(df)\n",
      "n = 200\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:  9\n",
        "Bonferroni:  0\n",
        "FDR:  0\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Adding signal :\n",
      "z = st.norm.rvs(0,1, size=n)\n",
      "print sum(z>1.5)\n",
      "z = z[z>1.5]+1\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": [
        "13\n",
        "number of true positive:  16\n",
        "Uncorrected:  0.05 16\n",
        "Bonferroni:  0.000512801416262 0\n",
        "FDR: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.00463758728239 15 16\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEMCAYAAADEXsFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQ1XX+x/HX4ZLjBUHTVeFQqLCCKaCJpKVSreGySRdl\npN3NIhNz1y3drc1p/wj7zVY2U20T2yatl8wgLNuwUnI1cTUvmCg02SheKKRJM0UtK/D4/f1xLlyP\n54DnAvh8zDic8/1+zvf75kN9X+f7+d5MhmEYAgBc8QL8XQAAoGMgEAAAkggEAIANgQAAkEQgAABs\nCAQAgCQvBMIDDzygAQMGaOTIkU7bPPzww4qJiVFCQoL27t3r6RIAAO3g8UDIyspScXGx0/nr1q3T\noUOHVFlZqby8PM2dO9fTJQAA2sHjgTBhwgT16dPH6fy1a9fqvvvukyQlJyertrZWx48f93QZAIA2\nCvL1CmtqahQZGel4bzabdezYMQ0YMKBJO5PJ5OvSAKBLaO8NKPxyULl5sc42/oZh8M8w9OSTT/q9\nho7yryv2xezZhiZNMvTrXxu6776G16dPt97GPr21vmitXfP5Awcauuqqpv/Cwgz96letf8ZVzTNn\nOl9na/UMG2YoNNRQv36Ghg5teF1V1fC5xm0aT3f2+zXui9BQQ4GBhgICDF1/fUNb+zK7dTOUnGzI\nbDZ0443O+6rx+hq3tf+NIiOt01pbjrM67f3frVtDn7v6m0dGGhowwPnfq/m6Lnfj7HFHjx41RowY\n0eq8OXPmGAUFBY73w4YNM7755psW7bxUWqf05JNP+ruEDqMr9sWkSYYhWf/179/wOiOj9Tb26a31\nRWvtnM1v7V9rn2lvzc7qCQ1tmGYyNbw2mxs+17hN4+nOfr/GfREY2Prv1XiZ7v7erfVX49/X2XKc\n1elqea7atra+5uu6nG2nz/cQ0tPTtXLlSknSzp07FRYW1mK4COhKsrOllBQpLU2qrW05v0cP68+k\nJCkhoeF1Xl7rbRpPv9SyWmtnn9+aUaMuvey21uysnuDghnlhYQ2vt21r+FzjNo2nu9MPAc22ava2\n9mXa54eGul6WfX2N29p/X/u01pbjrM7m/T9qlOu/uX35zdn/Xu7+t+GWdkeJE5mZmcagQYOM4OBg\nw2w2G0uXLjVeffVV49VXX3W0+eMf/2gMHTrUiI+PN/bs2dPqcrxQWqe1efNmf5fQYXTGvnD1rf30\naev006ebvnbWxq61vnD2+cbz77jDMH7zG+u/tDTrzzvucP4ZZ8txVbOzeqqqrN/6q6qavm7M2XRn\n62rcFxUVhtG9u2Fs29a0rX2ZFRXW6VVVl+6rxutr3Lb5tNaWc6m/o73/7X3u6m9eVWUYd97Z8LdK\nS7O+t7dv/vnL2XaabAvocEwmkzpoaeigsrOlgwet35jy8xu+fV6q7eHD0rXXSr17N/1Mdrb0/vvS\niRPWb5ShodKePVJqqvTNN9JPP1m/ofXp0/RzsbHW+cHB0qefWpedliatX2/9Brdhw6XrAi7X5Ww7\nCQR0GSkp0pYt1tcZGdLq1e61tWv8mdbmm83SuXPSmTPOPxcW1jDfbJaqq63DRNnZ1t15wgDedjnb\nTm5dgS6jLWOprY0NX2qs1z6W3XwcuvnnWhv7DguzBgZhgI6OQECn4ergbH6+9du6O8My9rbl5a1/\nJj9fuuMO6bbbpPBwaf9+6/DPp59av/nv29f65+zz7e2BzoQhI3QabRkSAq5UDBmh02v+7T872/rN\nvG9fafJk6zSPnl4HoAX2ENAhNP/2f+JE04O6GRnWEODgLHBp7CGg02v+7b/xQV37BTgcnAW8iz0E\ndAjNT82srZWysqyXc61YQQgA7uI6BACAJIaMAAAeQCAAACQRCAAAGwIBACDJD4/QxJWh8Z1Hf/EL\nqarKvbuQAvAfAgFecfBgw4Vl/ftL335rfZ2dzS0ngI6KISN4hbtP1ALQcXAdAryi8YVmErecAHyF\nC9MAAJK4MA0A4AEEAgBAEoEAALAhEAAAkggEAIANgQAAkEQgAABsCAQAgCQCAQBgQyAAACQRCAAA\nGwIBACCJQAAA2BAIAABJBAIAwIZAAABIIhAAADZeCYTi4mLFxsYqJiZGixcvbjH/5MmTmjJlihIT\nEzVixAitWLHCG2UAANrA44/QtFgsGjZsmDZu3KiIiAglJSWpoKBAcXFxjjY5OTn6+eef9cwzz+jk\nyZMaNmyYjh8/rqCgoIbCeIQmALRZh3qEZmlpqaKjoxUVFaXg4GBlZmaqqKioSZtBgwbp7NmzkqSz\nZ8/q6quvbhIGAADf8/hWuKamRpGRkY73ZrNZu3btatJm9uzZuuWWWxQeHq5z585p9erVrS4rJyfH\n8TolJUUpKSmeLhcAOrWSkhKVlJR4ZFkeDwSTyeSyzdNPP63ExESVlJTo8OHDmjx5ssrLyxUSEtKk\nXeNAAAC01PzL8qJFi9q9LI8PGUVERKi6utrxvrq6WmazuUmb7du3KyMjQ5I0dOhQDR48WAcOHPB0\nKQCANvB4IIwZM0aVlZWqqqpSXV2dCgsLlZ6e3qRNbGysNm7cKEk6fvy4Dhw4oCFDhni6FABAG3h8\nyCgoKEi5ublKTU2VxWLRrFmzFBcXpyVLlkiS5syZoyeeeEJZWVlKSEjQxYsX9dxzz6lv376eLgUA\n0AYeP+3UUzjtFADa7nK2nZzreYXLzpYOHpSOHJGuuUbq3VvKz5fCwvxdGQBfYw/hCpeSIm3Z0nRa\nRobk5ExgAB1ch7owDZ1Ljx7Wn6Gh1p9JSVJenv/qAeA/BMIVLj/fukdQXm79uWEDw0XAlYohIwDo\nQjioDKc4aAzAXQRCF2Df6PfoIR06JJ04IQUHS59+ap1uP2hsv4A8O5uDxgBaIhC6gMYb/eBgqb7e\n+vqmm6SRI62vQ0OlM2c4aAzAOQ4qdyLZ2VJ4uNS3rzR5slRba51uP1MoKck6JGSftm0bB40BuI+D\nyp1I82sG7NcL1NZawyIvz7oXcNNN1jC49lq/lQrATy5n20kgdCJpadL69dbXo0ZJH3/Mt30ATREI\nV4jaWikrSzIMacUKwgBASwQCAEASt64AAHgAgQAAkEQgAABsCAQAgCQCAQBgQyAAACQRCAAAGwIB\nACCJQAAA2BAIAABJBAIAwIZAAABIIhAAADYEAgBAkotAsFgsevHFF31VCwDAjy4ZCIGBgcrPz/dV\nLQAAP3L5gJwFCxaovr5eM2bMUM+ePR3TR48e7d3CeEAOALSZV5+YlpKSIpPJ1GL65s2b27VCdxEI\nANB2PEITACDp8radQc5mPP/8846Ft+bPf/5zu1YIAOiYnAbCuXPnZDKZdODAAe3evVvp6ekyDEMf\nfPCBxo4d68saAQA+4HLIaMKECVq3bp1CQkIkWYMiLS1NW7du9W5hDBkBQJtdzrbT5YVpJ06cUHBw\nsON9cHCwTpw40a6VAQA6LpeBMHPmTI0dO1Y5OTl68sknlZycrPvuu++SnykuLlZsbKxiYmK0ePHi\nVtuUlJRo1KhRGjFihFJSUtpVPADAc5wOGR05ckRDhgyRJO3Zs0dbt26VyWTSxIkTNWrUKKcLtFgs\nGjZsmDZu3KiIiAglJSWpoKBAcXFxjja1tbW68cYb9dFHH8lsNuvkyZPq169f08IYMgKANvPKWUYZ\nGRnas2ePbr31Vm3atEnXX3+9WwssLS1VdHS0oqKiJEmZmZkqKipqEgj5+fmaNm2azGazJLUIAwCA\n7zkNBIvFor///e86cOCAXnjhhSaJYzKZnJ52WlNTo8jISMd7s9msXbt2NWlTWVmp+vp63XzzzTp3\n7pweeeQR3XvvvS2WlZOT43idkpLC0BIANFNSUqKSkhKPLMtpILz11lt67733ZLFYdO7cObcX6Oy6\nhcbq6+tVVlamTZs26fz58xo3bpxuuOEGxcTENGnXOBAAAC01/7K8aNGidi/LaSDExsZq4cKFio+P\nV1pamtsLjIiIUHV1teN9dXW1Y2jILjIyUv369VP37t3VvXt3TZw4UeXl5S0CAQDgOy7PMmpLGEjS\nmDFjVFlZqaqqKtXV1amwsFDp6elN2txxxx3atm2bLBaLzp8/r127dmn48OFtqxwA4FFO9xDavcCg\nIOXm5io1NVUWi0WzZs1SXFyclixZIkmaM2eOYmNjNWXKFMXHxysgIECzZ88mEADAz7i5nRdkZ0sH\nD0qHD0vXXit99ZV0zTVS795Sfr4UFubvCgF0VV652+maNWscC27tQPHdd9/drhW6XVgnDoSUFGnL\nltbnZWRIq1f7tBwAVxCvXIfw/vvvy2Qy6cSJE9q+fbtuueUWSdbnIIwfP97rgdDZ2PcKevSQ7Hf6\nCA2Vzpxp+JmUJOXl+bdOAHDG6UHlFStWaPny5aqrq9P+/fu1Zs0arVmzRp9//rnq6up8WWOncPCg\nda9g/XqpVy/rnkB5edOfGzYwXASg43J5ULm6uloDBw50vB8wYIC++uorrxbVGfXoYf2ZlCQtX96w\n4bcPDzFMBKCjcxkIv/rVr5Samqrf/va3MgxDhYWFmjx5si9q61Ty863DRnl57AUA6JxcnmVkGIb+\n85//OJ5/MHHiRN11113eL6wTH1QGAH/xykHlxgsfPXq0QkJCNHnyZJ0/f17nzp1zPDDnStb4QDKn\nkwLo7FxeqZyXl6eMjAw99NBDkqRjx47pzjvv9HphnUHjA8nZ2f6uBgAuj8tA+Oc//6lt27apd+/e\nkqRf/vKXPDHNpvGBZE4nBdDZuQyEbt26qVu3bo73Fy5ccOuOpleC/HxOJwXQdbgMhEmTJunvf/+7\nzp8/r//+97/KyMjQ1KlTfVFbhxcWZj2dlDAA0BW4PMvo4sWL+ve//60NGzZIklJTU/Xggw96fS+B\ns4wAoO28ci8ju5deekmPPPKIy2meRiAAQNtdzrbT5ZDRihUrWkxbvnx5u1YGAOi4nF6HUFBQoPz8\nfB09erTJMYNz587p6quv9klxAADfcRoI48eP16BBg/Ttt9/q0UcfdeyChISEKCEhwWcFAgB8gwfk\nAEAX4tVjCDt27FBSUpJ69eql4OBgBQQEOC5SAwB0HS4DYd68ecrPz1dMTIx++uknLV26VH/4wx98\nURsAwIdcBoIkxcTEyGKxKDAwUFlZWSouLvZ2XQAAH3N5t9OePXvq559/VkJCgv76179q4MCBjO0D\nQBfkcg9h5cqVunjxonJzc9WjRw8dO3ZMa9as8UVtAAAf4iwjAOhCvPKAnIyMDL399tsaOXJkqyus\nqKho1woBAB2T0z2Er7/+WuHh4aqqqmr1g1FRUV4siz0EAGgPr+whhIeHS7I+U3ngwIHq3r27JOnH\nH3/U8ePH27UyAEDH5fKg8vTp0xUYGNjwgYAATZ8+3atFAQB8z2UgWCwWXXXVVY733bp1U319vVeL\n6giys6WUFCktTaqt9Xc1AOB9LgOhX79+KioqcrwvKipSv379vFpUR3DwoLRli7R+vTUcAKCrc3na\n6aFDh/S73/1OX3/9tSTJbDbrjTfeUHR0tHcL8/NB5bQ0axgkJfHMZACdh1efmGb3/fffS5J69erV\nrhW1lb8DobbWumeQl0cYAOg8vHq309raWi1YsECTJk3SpEmT9Je//EVnzpxp18o6Gvtxgl69pKuu\nkgICpD59pMmTrfNXryYMAFw5XAbCAw88oN69e+vtt9/W6tWrFRISoqysLF/U5nX24wQ//CDV10uG\nYd0z2LiR4wYArjwuh4wSEhJUXl7ucprHC/PBkJH9OEFQkHThQsP0UaOkjz9m7wBA5+PVIaPu3btr\n69atjvfbtm1Tjx492rWyjiY/X8rIkMrKpPBw6bbbpDvuIAwAXJlc7iHs27dPM2fOdBw36NOnj15/\n/XWvP1fZ3weVAaAz8uoeQmJioioqKlRRUaHPPvtM+/btcxkGxcXFio2NVUxMjBYvXuy03e7duxUU\nFKR333237ZUDADzKZSD84x//0NmzZ9W7d28tWLBAo0eP1kcffeS0vcVi0bx581RcXKz9+/eroKBA\nX3zxRavtHn/8cU2ZMoU9AQDoAFwGwrJly9S7d29t2LBBp06d0sqVK7Vw4UKn7UtLSxUdHa2oqCgF\nBwcrMzOzyZXOdi+//LKmT5+u/v37X95vAADwCJeP0LR/e//www917733asSIEZdsX1NTo8jISMd7\ns9msXbt2tWhTVFSkjz/+WLt375bJZGp1WTk5OY7XKSkpSklJcVUuAFxRSkpKVFJS4pFluQyE66+/\nXrfddpuOHDmiZ555RmfPnlVAgPMdC2cb98bmz5+vZ5991nHww9mQUeNAAAC01PzL8qJFi9q9LJeB\nsGzZMu3du1dDhw5Vz5499d1332n58uVO20dERKi6utrxvrq6WmazuUmbPXv2KDMzU5J08uRJrV+/\nXsHBwUpPT2/v7wEAuEwuA2HHjh1KSEhQr1699MYbb6isrEzz58932n7MmDGqrKxUVVWVwsPDVVhY\nqIKCgiZtjhw54nidlZWlqVOnEgYA4GcuDyo/9NBD6tmzp8rLy/XCCy8oOjpaM2fOdNo+KChIubm5\nSk1N1fDhwzVjxgzFxcVpyZIlWrJkiUeLBwB4jssL00aNGqW9e/dq0aJFioiI0IMPPqjRo0errKzM\nu4VxYRoAtJlXnqlsFxISoqefflqrVq3S1q1bZbFYrognpgHAlcblkFFhYaG6deumZcuWaeDAgaqp\nqdFjjz3mi9oAAD7k9gNy7LZu3aqCggK98sor3qpJEkNGANAeXh0ykqSysjIVFBRo9erVGjx4sKZN\nm9aulQEAOi6ngXDgwAEVFBSosLBQ/fv3V0ZGhgzD8NgVcQCAjsXpkFFAQIBuv/125ebm6pprrpEk\nDR48WEePHvVNYQwZAUCbeeX21++++666d++uiRMn6qGHHtKmTZvYQANAF+byoPL333+voqIiFRQU\naPPmzZo5c6buuusu3Xbbbd4tjD0EAGizy9l2tukso1OnTumdd97RW2+9pY8//rhdK3QXgQAAbeez\nQPAlAgEA2s6rj9AEAFwZCAQAgCQCAQBgQyAAACQRCAAAGwIBACCJQAAA2BAIAABJBAIAwIZAAABI\nIhAAADYEAgBAEoEAALAhEAAAkggEAIANgQAAkEQgAABsCAQAgCQCAQBgQyAAACQRCAAAGwIBACCJ\nQAAA2BAIAABJBAIAwMYrgVBcXKzY2FjFxMRo8eLFLea/+eabSkhIUHx8vG688UZVVFR4owwAQBuY\nDMMwPLlAi8WiYcOGaePGjYqIiFBSUpIKCgoUFxfnaLNjxw4NHz5coaGhKi4uVk5Ojnbu3Nm0MJNJ\nHi4NALq8y9l2enwPobS0VNHR0YqKilJwcLAyMzNVVFTUpM24ceMUGhoqSUpOTtaxY8c8XQYAoI2C\nPL3AmpoaRUZGOt6bzWbt2rXLafulS5cqLS2t1Xk5OTmO1ykpKUpJSfFUmQDQJZSUlKikpMQjy/J4\nIJhMJrfbbt68WcuWLdMnn3zS6vzGgQAAaKn5l+VFixa1e1keD4SIiAhVV1c73ldXV8tsNrdoV1FR\nodmzZ6u4uFh9+vTxdBkAgDby+DGEMWPGqLKyUlVVVaqrq1NhYaHS09ObtPnqq6909913a9WqVYqO\njvZ0CQCAdvD4HkJQUJByc3OVmpoqi8WiWbNmKS4uTkuWLJEkzZkzR0899ZROnz6tuXPnSpKCg4NV\nWlrq6VIAAG3g8dNOPYXTTgGg7TrUaacAgM6JQAAASCIQAAA2BAIAQBKBAACwIRAAAJIIBACADYEA\nAJBEIAAAbAgEAIAkAgEAYEMgAAAkEQgAABsCAQAgiUAAANgQCAAASQQCAMCGQAAASCIQAAA2BAIA\nQBKBAACwIRAAAJIIBACADYEAAJBEIAAAbAgEAIAkAgEAYEMgAAAkEQgAABsCAQAgiUAAANh06EDI\nzpZSUqS0NKm21t/VAEDXZjIMw/B3Ea0xmUwKDDRksVjfZ2RIq1f7tyYA6OhMJpPau1nv0IEgWUsL\nDpZOnJDCwvxbEwB0dJcTCB16yMguOZkwAABv6xSBcOSIvyvwr5KSEn+X0GHQFw3oiwb0hWd4JRCK\ni4sVGxurmJgYLV68uNU2Dz/8sGJiYpSQkKC9e/c6LzBA2r7dG1V2HvzH3oC+aEBfNKAvPMPjgWCx\nWDRv3jwVFxdr//79Kigo0BdffNGkzbp163To0CFVVlYqLy9Pc+fObXVZISHWvYNrr/V0lQCA5jwe\nCKWlpYqOjlZUVJSCg4OVmZmpoqKiJm3Wrl2r++67T5KUnJys2tpaHT9+vMWyzp4lDADAV4I8vcCa\nmhpFRkY63pvNZu3atctlm2PHjmnAgAFN2lnPNIIkLVq0yN8ldBj0RQP6ogF9cfk8HgjubsSbnxbV\n/HMd9GxYAOiyPD5kFBERoerqasf76upqmc3mS7Y5duyYIiIiPF0KAKANPB4IY8aMUWVlpaqqqlRX\nV6fCwkKlp6c3aZOenq6VK1dKknbu3KmwsLAWw0UAAN/y+JBRUFCQcnNzlZqaKovFolmzZikuLk5L\nliyRJM2ZM0dpaWlat26doqOj1bNnTy1fvtzTZQAA2srws/Xr1xvDhg0zoqOjjWeffbbVNn/605+M\n6OhoIz4+3igrK/Nxhb7jqi9WrVplxMfHGyNHjjTGjx9vlJeX+6FK33DnvwvDMIzS0lIjMDDQWLNm\njQ+r8y13+mLz5s1GYmKicd111xmTJk3ybYE+5Kovvv32WyM1NdVISEgwrrvuOmP58uW+L9IHsrKy\njF/84hfGiBEjnLZpz3bTr4Fw4cIFY+jQocbRo0eNuro6IyEhwdi/f3+TNh9++KHx61//2jAMw9i5\nc6eRnJzsj1K9zp2+2L59u1FbW2sYhvV/jCu5L+ztbr75ZuM3v/mN8c477/ihUu9zpy9Onz5tDB8+\n3KiurjYMw7pR7Irc6Ysnn3zSWLhwoWEY1n7o27evUV9f749yvep///ufUVZW5jQQ2rvd9OutKzx5\nzUJn505fjBs3TqGhoZKsfXHs2DF/lOp17vSFJL388suaPn26+vfv74cqfcOdvsjPz9e0adMcJ2/0\n69fPH6V6nTt9MWjQIJ09e1aSdPbsWV199dUKCvL4yLjfTZgwQX369HE6v73bTb8GQmvXI9TU1Lhs\n0xU3hO70RWNLly5VWlqaL0rzOXf/uygqKnJc5d5Vr1lxpy8qKyt16tQp3XzzzRozZozeeOMNX5fp\nE+70xezZs/X5558rPDxcCQkJeumll3xdZofQ3u2mX6PTU9csdAVt+Z02b96sZcuW6ZNPPvFiRf7j\nTl/Mnz9fzz77rONWv83/G+kq3OmL+vp6lZWVadOmTTp//rzGjRunG264QTExMT6o0Hfc6Yunn35a\niYmJKikp0eHDhzV58mSVl5crJCTEBxV2LO3Zbvo1ELhmoYE7fSFJFRUVmj17toqLiy+5y9iZudMX\ne/bsUWZmpiTp5MmTWr9+vYKDg1uc4tzZudMXkZGR6tevn7p3767u3btr4sSJKi8v73KB4E5fbN++\nXX/7298kSUOHDtXgwYN14MABjRkzxqe1+lu7t5seOcLRTvX19caQIUOMo0ePGj///LPLg8o7duzo\nsgdS3emLL7/80hg6dKixY8cOP1XpG+70RWP3339/lz3LyJ2++OKLL4xbb73VuHDhgvHDDz8YI0aM\nMD7//HM/Vew97vTFggULjJycHMMwDOObb74xIiIijO+++84f5Xrd0aNH3Tqo3Jbtpl/3ELhmoYE7\nffHUU0/p9OnTjnHz4OBglZaW+rNsr3CnL64U7vRFbGyspkyZovj4eAUEBGj27NkaPny4nyv3PHf6\n4oknnlBWVpYSEhJ08eJFPffcc+rbt6+fK/e8e+65R1u2bNHJkycVGRmpRYsWqb6+XtLlbTc77CM0\nAQC+1SmemAYA8D4CAQAgiUAAANgQCAAASQQCOqlevXp5fJlffvmlCgoKPL7cS8nJydHzzz/v03UC\nzhAI6JS8cbX60aNHlZ+f7/HlXkpXvOoenReBgE6tpKREKSkpysjIUFxcnH7/+9875kVFRenxxx9X\nfHy8kpOTdfjwYUnS/fffrzVr1jja2W9rsHDhQm3dulWjRo1qcQ+ckpISTZw4UbfffrtiY2M1d+7c\nFrcGOHPmjKKiohzvf/jhB11zzTW6cOGCXnvtNY0dO1aJiYmaPn26fvzxR0c7eyikpKRoz549kqxX\nXw8ePFiSZLFY9Nhjj2ns2LFKSEhQXl7e5XYb0CoCAZ3evn379NJLL2n//v06cuSItm/fLsm6oQ0L\nC1NFRYXmzZun+fPnO6a3ZvHixZowYYL27t2rRx55pMX83bt3Kzc3V/v379fhw4f17rvvNpkfGhrq\nuI+OJH3wwQeaMmWKgoKCNG3aNJWWlmrfvn2Ki4vT0qVLWyzfZDK1WtvSpUsVFham0tJSlZaW6rXX\nXlNVVVVbughwC4GATm/s2LEKDw+XyWRSYmJik43lPffcI0nKzMzUjh07LrkcV9dojh07VlFRUQoI\nCNA999yjbdu2tWgzY8YMFRYWSpLeeustzZgxQ5L02WefacKECYqPj9ebb76p/fv3u/37bdiwQStX\nrtSoUaN0ww036NSpUzp06JDbnwfc1fVuFI4rTrdu3RyvAwMDdeHChVbb2b99BwUF6eLFi5Kkixcv\nqq6uzq31NP72bhiGTCaT3nvvPS1atEiS9Zv81KlT9cQTT+j06dMqKyvTLbfcIsk6TLV27VqNHDlS\nr7/+umMvorHGdf30009N5uXm5mry5Mlu1Qm0F3sI6NLs39YLCws1fvx4SdZjC/ax+rVr1zruARMS\nEqJz5845XVZpaamqqqp08eJFrV69WhMmTNCdd96pvXv3au/evRo9erR69eqlpKQkPfzww5o6daoj\nRL7//nsNHDhQ9fX1WrVqlWO60ejW3VFRUfr0008lSe+8845jvampqXrllVccQXfw4EGdP3/eY30E\n2BEI6JQaf1u/1Jk6p0+fVkJCgl5++WW9+OKLkqwPUdmyZYsSExO1c+dOxymsCQkJCgwMVGJiYouD\nyiaTSUmaLV/zAAAAtElEQVRJSZo3b56GDx+uIUOG6M4772x1nTNmzFB+fr5juEiS/u///k/Jycm6\n6aabFBcX12S59vofffRR/etf/9Lo0aP13XffOaY/+OCDGj58uEaPHq2RI0dq7ty5TveCgMvBze3Q\nZQ0ePFh79uzxyN0uS0pK9Pzzz+v999/3QGVAx8QeArosT57j7+wMIKArYQ8BACCJPQQAgA2BAACQ\nRCAAAGwIBACAJAIBAGBDIAAAJEn/D2V29zbOX7dsAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3d7fb50>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "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.44980773  0.3049634   0.23894538  0.95284122  0.85829427  2.60215513\n",
        "  1.47310698  0.4147775  -0.09161498  0.37855162  1.17150395  0.69070465\n",
        "  3.3456306   0.89502638  0.68098054]\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "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.44980773 -0.3049634   0.23894538 -0.95284122 -0.85829427  2.60215513\n",
        "  1.47310698  0.4147775   0.09161498  0.37855162  1.17150395  0.69070465\n",
        "  3.3456306   0.89502638  0.68098054]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "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": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d_H0 = distrib_under_H0(z, my_stat)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size(d_H0)\n",
      "hist(d_H0)\n",
      "mean(z)\n",
      "empirical_p = 1.0 * sum(d_H0 > mean(z))/size(d_H0)\n",
      "print my_stat(z), empirical_p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.957711624634 0.0\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD9CAYAAABZVQdHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEMlJREFUeJzt3X9MVfXjx/HXSdinNSFlyZVd/O42gSGISDm0P9hoci3d\nvMOPG4mt7tJay1pr9a30r6CV4h/9UTpn6+vada3UtgK2ktEfXWuWXZe6tV0n1LDgCnfqhYm6Rdr5\n/mG7C4Hr5XJ/cN8+Hxsb3Hvufb89nvvkcjicY9m2bQsAYKR7Mj0BAEDqEHkAMBiRBwCDEXkAMBiR\nBwCDEXkAMFjMyPf39+vRRx9VZWWlli5dqg8++ECSFIlE5Ha7VVZWpjVr1mhkZCT6mF27dqm0tFTl\n5eXq7u5O7ewBADFZsY6THxoa0tDQkJYvX66rV6/q4YcfVnt7uz7++GM98MADeuONN7R7924NDw+r\nra1NwWBQmzdv1smTJxUKhdTQ0KCenh7dcw8/MABAJsSs78KFC7V8+XJJ0ty5c7VkyRKFQiF1dnbK\n6/VKkrxer9rb2yVJHR0dam5uVm5urlwul0pKShQIBFL8TwAATCUn3gXPnz+v06dPa+XKlQqHw3I4\nHJIkh8OhcDgsSbpw4YJWrVoVfUxxcbFCodC457EsKxnzBoC7TiInKIhrP8rVq1e1ceNGvf/++8rL\nyxt3n2VZMcM92X22bfORpI+33nor43Mw5YN1yfqczR+JumPk//rrL23cuFFPPfWUGhsbJd169z40\nNCRJGhwcVGFhoSTJ6XSqv78/+tiBgQE5nc6EJwcAmJmYkbdtW1u3blVFRYVeeeWV6O0ej0c+n0+S\n5PP5ovH3eDw6dOiQxsbG1NfXp97eXtXW1qZw+gCAWGLukz9+/Lg++eQTLVu2TDU1NZJuHSK5fft2\nNTU16cCBA3K5XDpy5IgkqaKiQk1NTaqoqFBOTo727dvHPvgUq6+vz/QUjMG6TC7W5+wQ8xDKlAxo\nWTPavwQAd6NE28kB7ABgMCIPAAYj8gBgMCIPAAYj8gBgMCKPu0J+fkH0r7PT9ZGfX5DpfzbAIZS4\nO9z6e410b3ds60geDqEEAExA5AHAYEQeAAxG5AHAYEQeAAxG5AHAYEQeAAxG5AHAYEQeAAxG5AHA\nYEQeAAxG5AHAYDEv5A2kQn5+gUZHhzM9DeCuwFkokXaZOiMkZ6FENuMslACACYg8ABiMyAOAwYg8\nABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMs1ACKZPzz8nY0icvb76uXImk\ndUzMbpyFEml3N52FkjNfIlk4CyUAYAIiDwAGI/IAYDAiDwAGI/IAYDAiDwAGI/IAYDAiDwAGI/IA\nYDAiDwAGI/IAYDAiDwAGixn5LVu2yOFwqKqqKnpbS0uLiouLVVNTo5qaGh09ejR6365du1RaWqry\n8nJ1d3enbtYAgLjEPAvl999/r7lz5+rpp5/WL7/8IklqbW1VXl6eXn311XHLBoNBbd68WSdPnlQo\nFFJDQ4N6enp0zz3jv49wFkpwFsrUjsnry0yJtjPm+eTr6up0/vz5CbdPNlBHR4eam5uVm5srl8ul\nkpISBQIBrVq1asKyLS0t0c/r6+tVX18/7YkDgMn8fr/8fv+Mnyehi4bs2bNHBw8e1IoVK/Tee+9p\n3rx5unDhwrigFxcXKxQKTfr4f0ceADDR7W+AW1tbE3qeaf/i9YUXXlBfX5/OnDmjoqIivfbaa1Mu\nm+6r4gAAxpt25AsLC2VZlizL0rPPPqtAICBJcjqd6u/vjy43MDAgp9OZvJkCAKZt2pEfHByMfv7l\nl19Gj7zxeDw6dOiQxsbG1NfXp97eXtXW1iZvpgCAaYu5T765uVnHjh3TpUuXtGjRIrW2tsrv9+vM\nmTOyLEsPPvigPvzwQ0lSRUWFmpqaVFFRoZycHO3bt4/dNQCQYVzIG2nHIZSpHZPXl5m4kDcAYAIi\nDwAGI/IAYDAiDwAGI/IAYDAiDwAGI/IAYDAiDwAGI/IAYDAiDwAGI/IAYLCELhoCYLbKSfuJAfPy\n5uvKlUhax0T8OEEZ0o4TlJk3Jq/p1OMEZQCACYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8\nABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwbj8310uP79Ao6PDmZ4G\ngBTh8n93OS7Fx5jJGJPXdOpx+T8AwAREHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQB\nwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMFjPyW7ZskcPhUFVVVfS2SCQit9utsrIyrVmzRiMjI9H7\ndu3apdLSUpWXl6u7uzt1swYAxCVm5J955hl1dXWNu62trU1ut1s9PT1avXq12traJEnBYFCHDx9W\nMBhUV1eXtm3bpr///jt1MwcA3FHMyNfV1Wn+/Pnjbuvs7JTX65Ukeb1etbe3S5I6OjrU3Nys3Nxc\nuVwulZSUKBAIpGjaAIB4TPvyf+FwWA6HQ5LkcDgUDoclSRcuXNCqVauiyxUXFysUCk36HC0tLdHP\n6+vrVV9fP91pAIDR/H6//H7/jJ9nRtd4tSzrn8vHTX3/ZP4deQDARLe/AW5tbU3oeaZ9dI3D4dDQ\n0JAkaXBwUIWFhZIkp9Op/v7+6HIDAwNyOp0JTQoAkBzTjrzH45HP55Mk+Xw+NTY2Rm8/dOiQxsbG\n1NfXp97eXtXW1iZ3tgCAaYm5u6a5uVnHjh3TpUuXtGjRIr399tvavn27mpqadODAAblcLh05ckSS\nVFFRoaamJlVUVCgnJ0f79u2LuSsHAJB6lm3bdloHtCyleUjEcOsbcbr/PxjTtDF5Tadeou3kL14B\nwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGBE\nHgAMRuQBwGBEHgAMRuQBwGBEHgAMRuQBwGAxL+SN9MrPL9Do6HCmpwHAIFzIexbhotqMma1j8ppO\nPS7kDQCYgMgDgMGIPAAYjMgDgME4ugbADOX8c9BA+uTlzdeVK5G0jpmtiDyAGbqhdB/RMzqa3m8q\n2YzdNQBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAYj\n8gBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAYj8gBgMCIPAAZL+ELeLpdL+fn5mjNnjnJzcxUIBBSJ\nRPTEE0/o999/l8vl0pEjRzRv3rxkzhcAMA0Jv5O3LEt+v1+nT59WIBCQJLW1tcntdqunp0erV69W\nW1tb0iYKAJi+Ge2usW173NednZ3yer2SJK/Xq/b29pk8PQBghhLeXWNZlhoaGjRnzhw9//zzeu65\n5xQOh+VwOCRJDodD4XB40se2tLREP6+vr1d9fX2i0wAAI/n9fvn9/hk/j2Xf/nY8ToODgyoqKtLF\nixfldru1Z88eeTweDQ8PR5cpKChQJBIZP6BlTfgJALdYliUp3euGMRkzO8e82zqSaDsT3l1TVFQk\nSVqwYIE2bNigQCAgh8OhoaEhSbe+CRQWFib69ACAJEgo8tevX9fo6Kgk6dq1a+ru7lZVVZU8Ho98\nPp8kyefzqbGxMXkzBQBMW0K7a/r6+rRhwwZJ0o0bN/Tkk09qx44dikQiampq0h9//DHlIZTsrpka\nu2sYkzHjH/Nu60ii7Ux4n3yiiPzUiDxjMmb8Y95tHUn7PnkAwOxH5AHAYEQeAAxG5AHAYAn/xavp\n8vMLNDo6fOcFAWAW4+iaKXCkC2My5uweMxs6kkwcXQMAmIDIA4DBiDwAGIzIA4DBiDwAGIzIA4DB\niDwAGIzIA4DBiDwAGIzIA4DBiDwAGIzIA4DBiDwAGIzIA4DBiDwAGIzIA4DBiDwAGIzIA4DBiDwA\nGIzIA4DBiDwAGCwn0xMAgOnLkWVZaR0xL2++rlyJpHXMZCDyALLQDUl2WkccHU3vN5VkYXcNABiM\nyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOAwYg8ABiMyAOA\nwbLiVMN//vmnPv30U924cSPTUwGArJIVkQ8Gg3rhhf/VnDn/zfRUACCrZEXkJek///kfXbnyURpH\n/L80jgUAqcE++aznz/QEDOLP9AQM48/0BKAURL6rq0vl5eUqLS3V7t27k/30mMCf6QkYxJ/pCRjG\nn+kJQEmO/M2bN/XSSy+pq6tLwWBQn332mc6ePZvMIQAA05DUyAcCAZWUlMjlcik3N1ebNm1SR0dH\nMocAAExDUn/xGgqFtGjRoujXxcXF+umnnyYsZ1mJXvU83VdLz8TV2RMZszUDY87UbB1zpusykTGT\nbTaNmez1Gc+YKRwx4XZlTlIjH88KsG07mUMCAGJI6u4ap9Op/v7+6Nf9/f0qLi5O5hAAgGlIauRX\nrFih3t5enT9/XmNjYzp8+LA8Hk8yhwAATENSd9fk5ORo7969euyxx3Tz5k1t3bpVS5YsSeYQAIBp\nSPpx8mvXrtW5c+f066+/aseOHfr8889VWVmpOXPm6NSpU1M+juPr4xOJROR2u1VWVqY1a9ZoZGRk\n0uVcLpeWLVummpoa1dbWpnmWs1s829rLL7+s0tJSVVdX6/Tp02meYXa50/r0+/26//77VVNTo5qa\nGr3zzjsZmGV22LJlixwOh6qqqqZcZtrbpp1iZ8+etc+dO2fX19fbP//886TL3Lhxw168eLHd19dn\nj42N2dXV1XYwGEz11LLS66+/bu/evdu2bdtua2uz33zzzUmXc7lc9uXLl9M5tawQz7b21Vdf2WvX\nrrVt27ZPnDhhr1y5MhNTzQrxrM9vv/3WXr9+fYZmmF2+++47+9SpU/bSpUsnvT+RbTPlpzUoLy9X\nWVlZzGU4vj5+nZ2d8nq9kiSv16v29vYpl7U5kmmCeLa1f6/jlStXamRkROFwOBPTnfXife2yLcan\nrq5O8+fPn/L+RLbNWXHumsmOrw+FQhmc0ewVDoflcDgkSQ6HY8r/YMuy1NDQoBUrVuijj9J5YrfZ\nLZ5tbbJlBgYG0jbHbBLP+rQsSz/88IOqq6u1bt06BYPBdE/TGIlsm0n5xavb7dbQ0NCE23fu3Kn1\n69ff8fHZ+AcGqTTV+nz33XfHfW1Z1pTr7vjx4yoqKtLFixfldrtVXl6uurq6lMw3m8S7rd3+zpNt\ndHLxrJeHHnpI/f39uu+++3T06FE1Njaqp6cnDbMz03S3zaRE/ptvvpnR4zm+frxY69PhcGhoaEgL\nFy7U4OCgCgsLJ12uqKhIkrRgwQJt2LBBgUCAyCu+be32ZQYGBuR0OtM2x2wSz/rMy8uLfr527Vpt\n27ZNkUhEBQUFaZunKRLZNtO6u2aq/XIcXx8/j8cjn88nSfL5fGpsbJywzPXr1zU6OipJunbtmrq7\nu2P+tv5uEs+25vF4dPDgQUnSiRMnNG/evOguMowXz/oMh8PR134gEJBt2wQ+QQltm8n5nfDUvvji\nC7u4uNi+9957bYfDYT/++OO2bdt2KBSy161bF13u66+/tsvKyuzFixfbO3fuTPW0stbly5ft1atX\n26Wlpbbb7baHh4dt2x6/Pn/77Te7urrarq6utisrK1mft5lsW9u/f7+9f//+6DIvvviivXjxYnvZ\nsmVTHhWGW+60Pvfu3WtXVlba1dXV9iOPPGL/+OOPmZzurLZp0ya7qKjIzs3NtYuLi+0DBw7MeNu0\nbJtfewOAqWbF0TUAgNQg8gBgMCIPAAYj8gBgMCIPAAYj8gBgsP8HiFS09+HJC44AAAAASUVORK5C\nYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3fa1e90>"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "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": {}
  }
 ]
}