statfail-simul.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def run(n=20, effect_size=.5, confound_size=.5, runs=10000):\n",
    "    effect = np.random.randn(runs, n) + effect_size\n",
    "    confound = np.random.randn(runs, n) + confound_size\n",
    "    group1 = np.random.randn(runs, n)\n",
    "    group2 = np.random.randn(runs, n)\n",
    "\n",
    "    unconfounded_results = scipy.stats.ttest_1samp(group1 + effect, axis=1)[1]\n",
    "    confounded_results = scipy.stats.ttest_1samp(group1 + effect + confound, axis=1)[1]\n",
    "    stim_tests = scipy.stats.ttest_1samp(confound.T, 0)[1]\n",
    "\n",
    "    sig_stimtests = np.where(stim_tests < .05)\n",
    "    insig_stimtests = np.where(stim_tests > .05)\n",
    "    assert(len(stim_tests) == runs)\n",
    "\n",
    "    # Fake effect (not detected) \n",
    "    # stimuli not classified as confounded\n",
    "    # significant effect with confounds\n",
    "    nonrejected_insign_in_conf = np.intersect1d(insig_stimtests, np.where(confounded_results < .05)[0])\n",
    "    no_alarm = unconfounded_results[np.intersect1d(nonrejected_insign_in_conf, \n",
    "                                                   np.where(unconfounded_results > .05)[0])]\n",
    "\n",
    "    # Fake effect (not detected) \n",
    "    # stimuli classified as confounded\n",
    "    # significant effect with confounds\n",
    "    rejected_sign_in_conf = np.intersect1d(sig_stimtests, np.where(confounded_results < .05)[0])\n",
    "    good_alarm = unconfounded_results[np.intersect1d(rejected_sign_in_conf, \n",
    "                                                     np.where(unconfounded_results > .05)[0])]\n",
    "\n",
    "    # Real effect (falsely rejected) \n",
    "    # stimuli classified as confounded\n",
    "    # significant effect without confounds\n",
    "    false_rejection = unconfounded_results[np.intersect1d(sig_stimtests, \n",
    "                                                          np.where(unconfounded_results < .05)[0])]\n",
    "\n",
    "    return no_alarm, good_alarm, false_rejection\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_hist(no_alarm, good_alarm, false_rejection, runs=100000, \n",
    "              effect_size=.5, n=20, confound_size=.3, ax=None, plty=False, pltx=False):\n",
    "    ax.hist([no_alarm, good_alarm, false_rejection], \n",
    "             label=('Nondetected false results %: {}'.format(int(10*len(no_alarm)/runs*100)/10), \n",
    "                    'Detected false results %: {}'.format(int(10*len(good_alarm)/runs*100)/10), \n",
    "                    '(Falsely) rejected \"real\" results %: {}'.format(int(10*len(false_rejection)/runs*100)/10)),\n",
    "             bins=np.arange(0, 1.05, .05))\n",
    "    ax.legend()\n",
    "    if pltx: ax.set_xlabel(\"p value of unconfounded data\")\n",
    "    else: ax.set_xticks([])\n",
    "    if plty: ax.set_ylabel(\"Counts\")\n",
    "    ax.set_title(\"n: {}. Effect size: {}. Confound size: {}\".format(n, effect_size, confound_size))\n",
    "\n",
    "\n",
    "def run_stim(n=20, stims=40, effect_size=.5, confound_size=.3, runs=100000):\n",
    "    effect = np.random.randn(runs, n) + effect_size\n",
    "    confound = np.random.randn(runs, n) + confound_size\n",
    "#    group1 = np.random.randn(runs, n)\n",
    "#    group2 = np.random.randn(runs, n)\n",
    "\n",
    "    unconfounded_results = scipy.stats.ttest_1samp(effect, 0, axis=1)[1]\n",
    "    confounded_results = scipy.stats.ttest_1samp(effect + confound, 0, axis=1)[1]\n",
    "    stim_tests = scipy.stats.ttest_1samp(confound.T, 0)[1]\n",
    "\n",
    "    sig_stimtests = np.where(stim_tests < .05)\n",
    "    insig_stimtests = np.where(stim_tests > .05)\n",
    "    assert(len(stim_tests) == runs)\n",
    "\n",
    "    # Fake effect (not detected) \n",
    "    # stimuli not classified as confounded\n",
    "    # significant effect with confounds\n",
    "    nonrejected_insign_in_conf = np.intersect1d(insig_stimtests, np.where(confounded_results < .05)[0])\n",
    "    no_alarm = unconfounded_results[np.intersect1d(nonrejected_insign_in_conf, \n",
    "                                                   np.where(unconfounded_results > .05)[0])]\n",
    "\n",
    "    # Fake effect (not detected) \n",
    "    # stimuli classified as confounded\n",
    "    # significant effect with confounds\n",
    "    rejected_sign_in_conf = np.intersect1d(sig_stimtests, np.where(confounded_results < .05)[0])\n",
    "    good_alarm = unconfounded_results[np.intersect1d(rejected_sign_in_conf, \n",
    "                                                     np.where(unconfounded_results > .05)[0])]\n",
    "\n",
    "    # Real effect (falsely rejected) \n",
    "    # stimuli classified as confounded\n",
    "    # significant effect without confounds\n",
    "    false_rejection = unconfounded_results[np.intersect1d(sig_stimtests, \n",
    "                                                          np.where(unconfounded_results < .05)[0])]\n",
    "\n",
    "    return no_alarm, good_alarm, false_rejection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# n=20\n",
    "stims=40\n",
    "effect_size=.5\n",
    "confound_size=.3\n",
    "runs=100000\n",
    "\n",
    "plot_hist(*run_stim(n=n, stims=stims, effect_size=effect_size, confound_size=confound_size, runs=runs), \n",
    "          n=n, effect_size=effect_size, confound_size=confound_size, runs=runs,\n",
    "          ax=plt.gca(), pltx=True)\n",
    "\n",
    "plt.gcf().savefig(\"/Users/jona/dropbox/manuscripts/statfail/statfail_simul.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAFtCAYAAAAu8eIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8FHX+x/HXphHYFIgEAYFQpCkEjcBJFawogoAGTSge\nqCBnEI1SA4RQBCkCggQb+jMg5e4QG3enIFwUEERQEBQl9BYCwZhdyKbs/P7gWAnpgQUmeT8fDx8P\nMvOdz37nuzHvndmZ+VoMwzAQERER0/C41h0QERGRklF4i4iImIzCW0RExGQU3iIiIiaj8BYRETEZ\nhbeIiIjJKLzluvLLL7/wxBNP0KNHD3r27EliYqJr3fr16+nevTtdunRh2LBh2Gy2IuuNGjWKjh07\n0qNHj1z/LV68GIBvvvmGzp07Ex4ezg8//MC9995Lr169OHr0aIn7vmPHDmJjY0u0zaBBg0hKSirx\naxXl3LlzvPTSSzz00EN06dKFNWvWXFa7SzkcDubMmUPPnj3p0aMH3bp14+23376sPttsNp544gm6\ndevGl19+eVm1irJz507uvvvuYrdPTk7miSeecEtfDhw4QGRkJF27diU8PJx9+/bl2+6rr77iL3/5\ni+t3uGfPntjtdrf0SUzAELmOdOvWzVizZo1hGIbx66+/GrfffruRlZVlnD592mjTpo1x8OBBwzAM\nY8aMGcaECROKrDdq1Chj0aJFha6Pj483DMMw5s2bZ8TExJS67//85z+NwYMHl3r7K+nVV181xo0b\nZxiGYRw7dsxo3769ceLEiVK3u5jT6TQGDhxoTJ482XA4HIZhGMaZM2eM3r17G3PmzCl1n7ds2WLc\nd999pd6+JHbs2GF07tz5qrxWUR599FHjs88+MwzDMP773/8aXbt2zbfdzJkzjTfffPNqdk2uYzry\nFrfYvHkzTzzxBCNGjKBnz5507dqVzZs3A+ePenr06JHvdqtWreKee+4B4ODBgwQGBuLp6ck333xD\naGgoderUASAiIoJPP/20WH0xCngO0TvvvMNXX33Fhx9+yAMPPMCyZctYu3Ytw4cPB+Dvf/87vXr1\nomfPngwYMMB1RGS32xk9ejQPPPAAXbt2Zfbs2Zw4cYLXX3+drVu3MmbMmDyv9eGHH/LII4/w2GOP\n0adPH9fR9t13381PP/3EkiVLcp0ZaNasGa+//jpw/oird+/e9OzZk4iICH744Qfg/NFgjx49SElJ\nyfN6a9eupXfv3gDUqFGD9u3b869//avU7S723XffsX//fsaMGYOPjw8AlStXZvr06bRq1QqAEydO\n8Oyzz9KtWze6devGu+++C8CRI0e49957mTx5MuHh4dx///2sXr2a/fv3ExMTQ3JyMj179sThcLBm\nzRp69uxJ9+7diYyMZMeOHQDMmzePSZMmufpz8c/9+vXjtddeo2/fvtx9992MGDHC9f5feJ8fe+wx\nPvzww3z3LTs7m9jYWLp160avXr0YNmwYZ8+e5ciRI9x+++0AjB492vU+PfTQQzRp0oRvv/0WgPj4\neHr16kWPHj147rnnOHnypGucBw0alOf1kpOT2b9/P127dgWgY8eOnD17lt27d+dpu337djZt2kSv\nXr3o06cPW7duLfR9kjLuWn96kLLp22+/NW655Rbj559/NgzDMBYtWmT07du3WNs6nU7jnnvuMZo2\nbWosXrzYMAzDePPNN43x48e72mRlZRmNGzc2bDZbobVGjhxpdOjQwXjkkUdy/ffrr78ahpH7yHze\nvHnGpEmTDMMwjM2bNxt9+vQxzp07ZxiGYXz99dfGQw89ZBiGYbzyyitGdHS04XQ6jczMTKNv377G\n5s2bjZUrV+Z75J2dnW00a9bMSElJMQzDMFatWmWsWLHCMAzD6Ny5s/HTTz/lar9kyRLjscceM86d\nO2fs37/fePjhh43ff//dMIzzZyPatWtnnD17ttD9bt68uXHq1CnXz7NnzzamTp1a6nYXe/fdd40X\nXnih0DZ9+vQx3nvvPcMwDCM9Pd3o3r278fnnnxuHDx82GjdubKxfv94wDMP4z3/+4zoC3rx5s/Hw\nww8bhmEYe/fuNdq1a2ccPnzYMAzD2LRpk9GuXTsjPT3dmDdvnjFx4kTXa138vvXt29fVN5vNZnTo\n0MH49ttvjd27dxtt27Z17WtcXFy+R97fffed8eCDD7p+njFjhrF9+3bj8OHDxm233ZanfXR0tBEX\nF2cYhmF89NFHxosvvmhkZ2cbhmEYy5YtM5555plCx2n79u1Gly5dci2LiIgwvvrqqzxto6KijC+/\n/NIwDMPYunWr0bp16yLPkkjZ5XWtPzxI2VWzZk2aNGkCQNOmTVm5cmWxtrNYLKxZs4YjR47Qp08f\nGjRoUODRs6enZ5G1BgwYwIABA4p8XcMwXK+zfv16Dh48mOt7zrS0NNLS0ti0aROjR4/GYrHg7e1N\nQkICQIH75+npSZcuXXj88cfp1KkT7dq1o1OnTvm2/fLLL3nvvfdYtmwZvr6+bNiwgZSUFJ588slc\n9Q4dOkTjxo0L3Ben05lnmYdH3hNtxW136fqcnJwC1589e5bt27fz3nvvAeDn5+e6fqFFixZ4eXlx\n1113Aed/L37//Xcg9xmSb7/9ljZt2lCrVi0A7rzzTm644QZ27dpVaN8AOnfuDIDVaiUkJIS0tDR2\n7dpF+/btueGGGwB4/PHHWb9+fZ5tGzdujKenJ+Hh4bRv357777+f0NBQjhw5kqft1KlTsdvtzJw5\nE4B169axc+dOHn30UQBycnJwOByF9jW/8Yf834N58+a5/n3HHXdw++23s2HDBnr16lXoa0jZpPAW\nt/H19XX922KxFNk+KyuLL774goceegiLxUKtWrVo27YtP//8MzVq1ODHH390tU1OTiYwMDDXa1xJ\nhmHwyCOP8PLLL7t+Pn78OIGBgXh55f7fJjk5mQoVKhRab8aMGezdu5cNGzbw9ttv849//IMFCxbk\navP9998zceJE3n//fVfIGIZBmzZtmD17tqvdsWPHqF69eqGvV7NmTU6ePOmqk5yczC233FLqdhe7\n7bbb+OCDD3A6nblCZseOHSxevJjY2NhcH4TgfEhlZ2cD4O3t7VpusVgK/GB26fKLa1y8LjMzM1e7\nS38nDMPAw8MjV1AW9KHP39+fjz/+mG3btvHtt9/y4osv0q9fP+69995c7RYtWsTWrVtZsmSJ63fb\nMAwGDRrk+sCXmZnp+mBSkJo1a3Lq1Klcy5KTk/O8v+np6SxZsoRnn302135dPJZSvug7b7lueHt7\nM3fuXD7//HPg/B+xzZs307p1a9q3b8+PP/7IwYMHAVi2bJnru/GiFBQOha1v164dn3/+uev75BUr\nVjBw4EAA2rRpw6pVqzAMg8zMTIYOHcrWrVvx8vIiKysrT63U1FQ6depEYGAgTz75JMOGDWPPnj25\n2uzdu5dhw4bx2muv0aBBA9fyO++8kw0bNri+b09MTKRHjx55AutS99xzD8uXLwfOf//89ddf53u0\nX9x2F7vtttuoV68eU6dOdfXj1KlTTJo0idq1a2O1WmnRooXre+X09HQ+/vhj2rVrV+R7ccFf/vIX\nNmzYwOHDhwHYtGkTycnJtGjRgqCgINcR+NmzZ/nmm29ybXvpa1gsFtq2bcuGDRtITk4GCj5Lsm7d\nOp588kluv/12oqKi6NGjB3v27Mn14fPTTz/lww8/ZOHChbk+KLRv354VK1a47oKYP38+o0aNKnQ/\nq1evTu3atVm9ejUAX3/9NZ6ennnOqlSqVIkPP/yQL774AoDdu3ezc+dOOnToUGh9Kbt05C1X3c6d\nOxk3bhyrVq3Ks27+/PlMnDiRd955Bw8PD0aMGMGtt94KwCuvvMLzzz9PVlYWderUYfr06UXWA3j/\n/ff55JNPci277bbbmDBhAvDnWQGLxeL6d/v27Xn66acZOHAgFosFf39/3njjDQCioqKYMmUK3bt3\nx+l08tBDD3Hvvfdy+PBh5syZw9ChQ3Od4gwKCmLIkCH89a9/pUKFCnh5eTF58mTXaxqGwdSpU8nJ\nyWHatGmuU9LNmzdn0qRJTJw4kejoaAzDwMvLi/j4eHx9fUlOTmbw4MG8/fbbBAcH59q/oUOHMmHC\nBB5++GFycnIYOXIktWvXBmDs2LE0a9aMJ554otjtLjVv3jxee+01evXqhaenJ06nk549e7o+4Myc\nOZOJEyfyz3/+k6ysLLp3707Pnj05cuRInrMw+Z2Vufnmm4mNjWXo0KHk5ORQsWJF4uPj8fPzo3v3\n7iQmJnL//fdz4403EhYWVmS9Ro0aMXz4cJ588kmsViuhoaH5trvrrrv4+uuvefjhh6lUqRKVK1dm\n0qRJOJ1OV/sxY8ZQvXp1Bg8e7Dqaj4iIoHfv3iQnJ/P4449jsVioWbMm06ZNA85fsLZ8+XLeeuut\nPK85e/Zsxo4dS3x8PBUqVGDu3LmudT169GDKlCnceuutLFiwgMmTJzNv3jw8PT2ZM2cOlStXzlNP\nygeLUdyPwiIiInJdKNZp89OnT3PXXXexf/9+du/eTYcOHejXrx/9+vVz3VayYsUKHn300VwXgmRk\nZDB06FD69OnDoEGDSE1NdduOiIiIlBdFHnlnZWXxwgsvkJSUxIIFC/j++++x2Wy5rt5NSUlh4MCB\nrFy5EofDQUREBP/85z9ZsmQJdrudqKgoVq9ezfbt24mJiXH7TomIiJRlRR55T58+nYiICNd3art2\n7WL9+vX07duXmJgY7HY7O3bsICwsDG9vb/z8/AgJCWHPnj1s27aNjh07AtChQwc2bdrk3r0REREp\nBwoN75UrVxIUFET79u1dy0JDQxk5ciSLFy+mdu3azJ8/H7vdjr+/v6uN1WrFZrNhs9mwWq2uZenp\n6W7aDRERkfKjyPDeuHEj/fr145dffnFN8nDhPtD77ruPn3/+GT8/v1wPyL8Q5hcvt9vtBAQEFNkh\nXT8nIiJSuEJvFbsw8xKcf2ZwXFwczz33HDExMYSGhrJx40aaNWtGaGgos2fPJjMzE4fDQVJSEo0a\nNSIsLIzExERCQ0NJTEykZcuWRXbIYrGQkqIjdHcKDvbXGF8FGmf30xi7n8bY/YKD/YtudIkS3edt\nsViIi4sjLi4OLy8vqlWrxsSJE7FarfTv35/IyEicTifR0dH4+PgQERHByJEjiYyMxMfHh1mzZpW4\ngyIiIpLbdXmftz7luZc+SV8dGmf30xi7n8bY/Upz5K3Ho4qIiJiMwltERMRkFN4iIiImo4lJRESk\nVDIzMzl8+GC+62rXDsHHx+cq96j8UHiLyBVV2B/00lIQXJ8OHz7IiE/GY73kgit7SjrTu0+kQYOG\n16hnZZ/CW0SuqMOHDzJsxidUCqx2ReqdTTvJ3OHdCw2Cbdu2MmbMy3zwwXKqVbsRgPj4edStW48H\nH3y41K/94YcJVKlSpcAa+/btJT09nRYtbi92zczMTL74YjUPP9yjWO1fe+1VOne+l9tvv8O17OjR\nIwwfPoxmzUIZMyY2zzbHjx9jwoQY3nzzvWL3q7Sswf741yx8atLHHuvG0qUrSU09zd69v9GuXdHz\nkH/22So+/fRjGjVqwksvjQQgLm4sw4ePplIla4Hb/fHHH0yeHIvNlo6vry8jRoylevXqudps2rSB\n999/B4CmTW/hhReGF9mf642+8xaRK65SYDX8qtx0Rf4r7ocAb28fXnklzvVzfvN1l1RRNdatW8v+\n/ftKVPP06VN8+unHl9WHHTt+oG3bDvkG9/Xowrz133//HTt3/lisbf7979W8+eZ7nDp1kvT0dDZu\n/IYWLW4vNLgBEhLeo3nzFixY8A6Rkf2ZO3dGrvVnz9qJj3+dGTPm8Oab71G1ajBnzpwp9b5dKzry\nFhHTs1gshIW1BAz++c8VPPpo71zrly5dzFdffYGnpxctWtzOkCFDeffdNzlx4jhnzqRy4sQJnn8+\nmtat7yQxcT3vv/8OgYGBWCwW7ruvCwALF85nx44fcDqdPP54JM2bt+Bf//oMHx8fGjduQkZGBm+/\nHY+Hhwc33VSL4cPHkJOTzSuvxJGcnExWVhYvvjiCzz//mAMH9vH+++8QHv4EU6dO5I8//gDghRde\npn79m1m16h988slHVK4cREbGOTp1use1LydOnGDx4vfJyMigVq1ahITU4/3338HpdHLu3DliYyfj\n5fXnn/Y333yDH374nuzsHDp1ups+fZ4kKWkvc+fOxDAMAgMDGT16PFarn2ubd999k59+2kFGxjle\nfXUaX3zxFWvWfIHFAvfccz+PPfYE//3vVyxa9BbH7MepUNVKSPitJK/bj7d/BXxusnLs2FHmzJnJ\nvHlvAuB0Olm8+H0cDgfNmoVy8mQy//7353h4eNCkyS288MLLud4zX19fHA4HOTk5eHhYWL36E+Li\npuZqs2DB63TufA9Nm97qWnbgwD4GDfobAM2bhzJ27Mhc2+zcuYP69W9m3rzZHDt2lIcffoQqVaqU\n6PfteqDwFhHTu/CsqZdeGsUzzzzJnXe2da1LStrLunVrWLjwPTw9PYmJGc7Gjd9gsVjw8fFh5szX\n+e67zSxbtoQ77mjFvHmzeffdBAICAoiLGwucP816/PgxFix4B4fDwbPPDmDevLd46KFu3HBDVZo2\nvZUnnujFwoWLqFy5Mu+8s5B//eszzp61U7NmLeLipnLkyGE2bvyGJ598in37kvjrX59mwYLXadmy\nNT16PMbhw4eYOnUiU6bMYMWKpXzwwXI8PDwYOnRwrqPv6tWr07fvXzl06CA9ejzGRx/9g3HjJlG1\nalUSEt5j3bo13H//g672a9b8h3nz3uKGG25g9epPAXj11cnExEwgJKQun322iiVLPnAFHpz/MFSv\nXn2ef/4l0tKS+eqrNcTHv/u/J2hG0bp1G9as+YKHHurG51n/JevkOZyOHCjkTIWHhwf9+g3g0KGD\ntG/fkWee6c9LL42mSZOmrFr1D3JycvD09HS1799/IFOmTKBjx8588cW/6dq1O0uW/B8nT56kd+8I\n6tQJ4W9/ez7P69x8cyO++SaRhg0b8803iTgcGbnWp6X9zvbtW3n//aX4+lbkueeeplmzUGrXrlOs\n37XrhcJbRMqMgIBAnn/+JSZPjqV58xYAHDp0gFtvbe4KhhYtbmf//iQAGjZsBEC1ajeSmengzJkz\n+PlZXZMoXfgue//+JPbs+YWhQwcDkJOTw/Hjx4DzHxzOnDlDauppxo07f5TncDho1eovpKX97vog\nUatWbXr3jnBtB+e/M9++fStr134JQHr6Hxw9epiQkHquo+fmzVvkO2HThWVVq1ZlzpwZVKpUiZSU\nk4SG3par3fjxk4iPf53U1NOuvhw8uJ+ZM88fxWZnZ+cbXHXqhADw66+/cuLEcZ5//lkAbLZ0jh49\nzNChL/LGG3M4tvtXrLUCCWwaXMS7c77PF/o9enQsy5Yt5tixozRrFppnH0NDbyM09DZsNhszZ07l\njjsG8e23G3nmmWeZO3cW48dPyvc1+vUbwJw5M4iKGkSbNu1c10BcYLVaqV07hNTU0wDUrVuPb775\nL61b32mqCyMV3iJSprRr14HExHX861+f8be/PU9ISF2WLVvyv9OvHvzww3a6dOnK3r2/ArmPFKtU\nqYLNZufMmVSqVAli9+6fuP32O6hTpy5hYXcwYkQM2dnZJCS8x0031cLDw8N16rlatWq8+uprVKpk\nJTFxPf7+/uzd+xs//7yb9u3v4ujRIyxa9CaDBj2H0+kEICSkHk2aNOW++7qQknKSL7/8N7Vq1WH/\n/n04HBn4+FTg55935TqTALlnX5w+/RVWrPiYihUrMmXKBFdtgKysLNatW0Nc3CsYhkG/fr25554H\nqFOnLuPGTaRatRv54YdtpKWl5TOS58emfv361KvXgFmzXgdg2bLF1K9/M5988hE9ejzG/oanSdt6\ngrSfU/Dw8iAr3YEPVg4c2J+nooeHh6t/n366ipdfHo2Pjw/R0UPZtWtnvhf+LV78Pv36DcDhyMDD\n4/wHsHPnzhb4/v/wwza6d+9Js2ahrF+/Nk/NSpUqsfPXHYxfNxUPH0+O/bCH3ZUP8o9PVpvqCnmF\nt4hccWfTTl7VWhaLJdep5WHDXuL7778DoH79m7n77nsZMuQpDMNJaOjtdOzYib17f821jcViwdPT\nk5dfHsXLLw/Dz8+fSpUqAtC+fUe2b/+e5557hnPnztKxY2cqVapE48ZNeOON1wkJqcuwYS/x8svD\nMAwnVqsfY8dO5NZbmzN16kSiogbhdDoZNuxlqlQJIjs7i4UL5/PkkwOZOnUSn3zyEXa7naeeGkzl\nypV58smBDBnyNAEBAXh65v0zffH+3n//gzz33NNUrRpMnTp1OX36lKuNt7c3AQGBDBr0VypUqEDr\n1ndSvXp1Xn55NJMmjScnJweLxcLo0ePzfQ2AJk2acMcdrRgy5CkyMzO59dZmBAdXo2nTW5k9ezrH\nHSfxDqhAQOOqODOyObDiJ9J+PcXB0Bp5xrdBg5v54INFNG7chAYNGvDcc09TqZKV4OBq3HJLszx9\nOHHiODabjQYNbsYwDJKTTzB8+AuuU/z5fecdElKXyZNjAQN//0DXRX3Lly/hpptqU6NGDW5oW4uT\nX5z/cBF0Ww1uaFaD9GO/F/VrVqirfc97mZiYxOl0smHtl5DPrrS4sw0BAYFXqmtlgiYauDrK6zhf\nzfu8y+sYX02FjXFS0m/EbZqR51ax9GO/E9tm+HV5FOuuPicl/Vbqe97dPiXo9SozM5Okd9+iRQXf\nXMt/z8zkl4oVad3hrmvUM5Hyx8fH57r8oy3mYsantxXnnvcrpUyEN4CvlxdWb+9cyzIv+u5HRETM\nQ09vK1yZCW8RESlbruaRrNnoCWsiIiImo/AWERExGZ02F5ErSrOKibifwltErqiCLjQqreJcoLRt\n21bGjx9NvXr1MQyDnJxswsMjufvuewvcpjzNCHY5V27ffffdLF78jxLNCPbHz6eY+OV4QkNvK9GM\nYDk5OcybN5s9e37GZrNxrin416yM7cDvHP9iLwDeVStCm9zb2Ww2Jk4cx9mzdrKzs4mKepFmzZoX\n2U8zU3iLyBV3tS80slgs3HFHK+LiXgHg3LlzREUNonbtOq5HoF5q3bq13HBD1RKF94UZwYob3oXN\nCBYV9UKxX/dyXe6V2xdmBDt06GCxwtv2ayqzpr3NO+8sJD09nZ07fyzWjGD/+c9qcnJyiI9/l61b\nNzNx+fn389i/f6PuE83wqVyRPQu3cPDggVx9Xr58Ca1a/YXw8Cc4dOggEybEsGjR4iL7aWYKbxEx\nvUufNVWxYkUeeaQX69evpWHDRuV2RrBRo8bz3Xff8tlnH/N7+jE876hJ8J21+X33SVK+OYQzO4c3\nDs9l1qx5LFr0FjfcUJUePR7l4MEDzJw5tcgZwTIzHZzy/R3/3rk/qFm8PMjMzCzxjGBbtnxL/foN\nGDHiBex2G9Zbz9dtOKglFg8LOY5snJk5+PpWzFXr8cf74P2/W4Wzs7OpUKFCyX6BTEjhLSJlUlBQ\nEL/++gvffrux3M4Itn//Pr76ag1jx8YRt3EGKV8eIODmG/h950mC29fBs7IPtznDsNvthc5dXtCM\nYN7eXrz0zggMp4HF48/tq4RV5+2347n77vtKNCNYWtrvHD16hOnT5/D5558w+/3XCGpyIxYPC/bD\naRz6+y68AnzyTOHp53f+w8vp06eYPHk8w4a9nKd2WaPwFpEy6fjx4wQHV2Pfvr3ldkawffuSOHHi\nONOmTeb4H4fACY7Us9TscjMnvz6I/XAavzUJwsMjd3AX1OdLZwTbt28v2RUz//do6j9r+Nbw47le\nw7jxxholmhEsICCQtm3bA9CkSVOy0v6cztNaO5Cm0W05uGo3n3/+Sa4jdjg/9euECWOIinqxRF+F\nmJXCW0TKHLvdxmefrWLy5OkcPHig3M4IFhJSl3r1GvC3vw0lbtMMMvb/ge+NfqRuPUb1zvU4l3YW\n9hj897/r8PHxcU1q8uuvv+SpmN+MYIcPH+TZsYOxH/4Dv7p5r3Eo6YxgoaG3sWnTBu66624OHTqI\nl9/5C+n2vvM99fqE4lnRGw9vDzw8ct/lvH//PsaNG8mkSa/SoMHNBdYvSxTeInLF2a/gZCHFqWWx\nWNi2bStDhw7Gw8OTnJxsnnrqWWrXrkPt2nXK7YxgN9/ckDvuaMXkyRM48vth/OtVwdu/ApVqBbB/\n8Q4MDCrf2Jh27Tpit9sYP34UP/ywjcaNmxZrRjCLxYJnRS8q1QrI04dTp1JKPCNY9+49mTlzKoMH\nDyAjI4OqHc+fVQhuX4d9CT9i8fIALwtdnuoKQHR0FNOnz+Gtt94gKyubOXNmAODn58/UqTOL/L0x\ns2LNKnb69Gl69erF+++/j4eHB6NGjcLDw4OGDRsSGxuLxWJhxYoVLF++HC8vL4YMGUKnTp3IyMhg\n+PDhpKamYrVamTZtGkFBQUV2qqSzBGVkZLDq6ScJu+QihjMOB5ZBz2pikktoJqaro7yOs2YVu/5c\nzkxa12pWMXfO/nW91XXLrGJZWVmMHz+eihUrYhgGU6dOJTo6mlatWhEbG8vatWtp0aIFCQkJrFy5\nEofDQUREBG3btmXp0qU0btyYqKgoVq9eTXx8PDExMSXupIiYh2YVE3G/Ih+POn36dCIiIggODgZg\n9+7dtGrVCoCOHTuyceNGdu7cSVhYGN7e3vj5+RESEsKePXvYtm0bHTt2BKBDhw5s2rTJjbsiIiJS\nPhR65L1y5UqCgoJo3749b775Zq6rDQGsVivp6enYbDb8/f1zLbfZbNhsNqxWa662IiKSPzPOYS3X\nRpHhbbFY2LhxI7/88gujRo3izJkzrvU2m42AgAD8/Pyw2+2u5Xa7HX9//1zL7XY7AQF5L2rIT0nP\n/2dkeONRwD2KgYEVS/V9QlmnMbk6NM7uV5bG+Ndffy3wSWjv/PU1brop/6fFFeXMGb8C1wUF+RU5\nhgWtv9y6hXFXbbPVLUih4b148Z+Pl+vXrx9xcXFMnz6dLVu20Lp1axITE2nTpg2hoaHMnj2bzMxM\nHA4HSUm0XhL2AAAgAElEQVRJNGrUiLCwMBITEwkNDSUxMZGWLVsWq1OluWDNWcB1d2lp53RByyV0\nkc/VoXF2v7I2xqmptgIfLZuaaiv1vqam2gpdV1jdwsb4cuoWxV21r8e6brlg7WIWi4VRo0Yxbtw4\nsrKyaNCgAV26dMFisdC/f38iIyNxOp1ER0fj4+NDREQEI0eOJDIyEh8fH2bNmlXiDoqIuWhWMRH3\nK3Z4JyQk5PvvC8LDwwkPD8+1zNfXl7lz515G90TEbA4fPsjGF5+nRqVKV6Te8bNnaTv79SKvYE9L\n+5233lrA8OFjiIoahMPhwNfX17V+9uw3cj03/IJt27by8ccrXZOaFMfq1Z9y6NBBnn02Kt/1U6ZM\n4OWXR5f4Gdtpab9z6utD+D9e/EldSjrT2ZHP9lD51mr41fvzEaNRUYPo2rU7ixa9xd///kmJ+nyx\nE1/tw9u/Ag5bBm//FI/dbmfgwEG5ZlW7ki5+7/7733XcemtzqlatWug2TqeTE/9O4kS2QfV76uPf\nIAhH6jlObTicZ7aygnzwwSKSkvbm+Z1xOnLYv/hHcjJzMHKc1OzSEDxLu3eF00NaROSKq1GpEnX8\nru530W+/Hc+jj/YGzp8lHDduousxoYUp7Jnepd3mvvu68OGHHzBgwDMlqhsYWJmqHfI+7rQwJZ3p\njHz6fuON1bn55obceGP1Er12QbW9rN5UpgqVKlkJDq52eTWL6e9/X0q9evWAwsP70KGDePlXoM7D\njTn00c/4NwjiZOIBqoQVb983bdrApk0b8h2r33ck43dzEMF31sZx6iwH/76LGo+457ZJhbeImJ7d\nbuOXX36mfv2LH42Z9zqYdevW8NFH/yA7OxuLxcIrr8zIdQfNK6/EcfToERwOB+HhT/DAAw+xffv3\neWYcu+CTTz7iyJFD/O1vw8jJyWHgwD68804Cd9zRinnzXssT3o8++jAhIfWoV68evXtHMmPGKzgc\nDipUqMCIETGkpKRw9KM9NHnuL9j2n+HE2n3gYcHD15Oc1jk4HBmlmunsq6++5MjnP1OhckWcmTlw\nS3Cufj3/fDSBgZWZPHl6nn4OGTKIUaPG5OpntWo3snDhfLZv/54jyXux1q5MnZ5NXfV8a/rT9fbu\n1K4dQuXKuc8i9OvXmzp1QvD29mH48NH59veVV+LYu/c3jp45xI0d61KlRXV2v7aRJs/fCcCKFUu5\n7bYwqlevAcCmTd+wd+9vTJ48gTlz3iAubix2ux2HI4NBg/5Gq1Z3/tk3X1+MHCfOLCcePp7YD6VR\n4YZKeFb0drX54480Xn11MlOmzMjV9yNHDvPJJx/x1FOD+fTTVXl+vyqH3oh/rfP7azidWLyLvBu7\n1BTeImJ6u3b9lOcoe9KkWNdp8y5dutK1a3eOHDnMjBlzqFDBlxkzXmHz5m9dz7A4e/YsP/64nbfe\neh84Pz0lwKuvTskz49iF0+/33vsAAwf25dlnh7J58ybCwlq5pqasUiWIpKS9uZ61nZJykvfe+5CA\ngADGjx/NY489wZ13tmXr1i0sXDifBx74c0axw5/8QsOn78DL6sPBj3fz9df/Zfv270s109kXX/yb\nm3o2wf+myiS9t+3iOUSA80f8gCtoL+7nlCnj8vTz5ZdHERAQwIgRY5iwcTrHVu4h6w+Hq57Fw0Kl\nSpXyBDecv8D4r399hoYNG+Xb35kzX+fHH7czZsx4pn83H8OWfb5mIe9/mzbtadiwEcOHj+HEiRP8\n8Ucas2bN48yZMxw6lPv6i+rVa+Bl9ebYF3u5sVNdTqzdT437G3D0P7/xj6PLGTEihoCAwDzBffbs\nWV577VXGjZvI/v378u2HRwVPPLw9yUp3cOifu6n5YKN8PkJeGQpvETG9tLTfqVIl96OX8zttXrly\nFSZPnkDFihU5dOggzZqFutZVqlSJ559/iVdfnYLdbueBBx4scMaxWrVqu7a5/fYwNm/exOrVnzJw\n4J9H2jfcUJW0tN9zvX5gYGXXLbP79iWRkPAeS5b8H4ZhuEIfINueSbYtkwPLfzr/89lMTt90ijNn\nUks101mNGjX53dOJxcOCtU7l/E5KFNjP3377jWPHTuTqp49PBVJTU4mPn8eptEM4M3MwLppUpSh/\nznqWt78X3odFi97m5In9VG1dq9h1AerXb0D37r2YMOH8RDSPPfZEnjZV7qiBf83KnNlxgsAmVUnd\neoyAplWxelvZunULrVr9Jc823323mdTUVMaPH43Nls6pU6dYsuT/6NPnyVztziXbOPT3XdR44Gb8\n6lYm/djveWpdCQpvETG9KlWCsNkuvRUnd0LZbDYWLXqLlSs//99dMVG5TpmfPn2KPXt+5pVXZuBw\nOHj00Ye5//4H851x7MSJ467tunXrweLF/8cff6TlOm2fnv4HQUE35OrDxVNvhoSEEBHRj2bNQtm3\nby+7d//kWudZyRvvgArnZ9Kq4MWJDQe49dZmnD17rlQznR09eoSK2dUxDIOzR/7Av2Hhc0xc3M/6\n9evTq9cTufr57bcbSUlJZsiQoYz/6hWO/P2XIj8Q5K7vUWB/L7wPzz8fTezXr3Jk6S6qtKiOxcuD\n7HQHhmFw6NABbrstLE9Np9PJvn17OXv2LNOnz+HUqVMMGfKUa5rRizmzckjbnULdJ5pzdPWvYDl/\nLcO5c+fy7fNdd3Xmrrs6A7B9+/esWvXPPMGdmXqOo18dJOTxZlS8seD7vq8EhbeIXHHHzxY87WNp\natUros2ttzYnPn7eJUtzn2j18/OjefMWDB48gCpVqlC7dginT5+iRo2aWCwWbrihKqmppxkyZCAe\nHp5ERvbDy8srz4xjMTFxnDhx3HXR2i23NOPo0SOui+Xg/BXNKSkp1K17ac//7NNzz73AzJnTyMx0\n4HA4eOGF4dhs57/7tVgs3PRgI/Yn/IhhgIHBTd1qccstpZvprHv3nvzfR+9zKsAXi2dxLtD7s82I\nESMYM2Zcrn7WqFGD//u/d5g2bRKnso5TqVYAWemZ+Q17obXz6++F92HSpFhO2I8R3K4OFg8L1drX\nYd/iH/Hw9aLFTc3+rPa/96FZs1CmTIll6tTX2L79e9atW4PT6eSZZ57Ntxenvj1C1TvPn0EJCqvB\nwX/sYmc1H55++tkCv/POtRcXXfgXHR3F4MHPkbrlGEaOk6Of/wqAZ0Uvqt5VsgsQi6tYs4pdbZpV\nzL3K2oMtrlfldZyv1axiM2dO5ZFHetGwYeMr+tqXunT/nE4nU6bEMXz4aBo2bIyPjw+bNn3Db7/9\nSv/+A0tUe8OGRGLnjKfxkNa5ll/rGa80q5h767r9IS0iIkW5VrOKPfXUs7z11gJGjnTvzIWHDx90\nPcI06w8Hyf/Zh3+TGxj/xStM951I/fo3s2bNfxgxomT9OHkymYUL38Bar/j3eEv5pfAWkTKhSpUq\nbg/uC1yPMK0JQU1uBHBdmHT+HvNJJa5ZrdqNTJgwmbhNBZ+qFbnAfTehiYiIiFsovEVERExG4S0i\nImIyCm8RERGTUXiLiIiYjMJbRETEZBTeIiIiJqP7vEWkTCrsSW8FPbFNxCwU3iJSJl38JLSL2VPS\nmd594jV5CpzIlaLwFpEyy/UkNJEyRt95i4iImIzCW0RExGQU3iIiIiaj8BYRETEZhbeIiIjJKLxF\nRERMRuEtIiJiMkXe552Tk8PYsWM5cOAAFouFuLg4srKyGDx4MHXr1gUgMjKSBx98kBUrVrB8+XK8\nvLwYMmQInTp1IiMjg+HDh5OamorVamXatGkEBQW5e79ERETKrCLDe926dXh4eLB06VK2bNnC7Nmz\n6dy5MwMHDmTAgAGudikpKSQkJLBy5UocDgcRERG0bduWpUuX0rhxY6Kioli9ejXx8fHExMS4dadE\nRETKsiLD+95776Vz584AHD16lICAAHbt2sX+/ftZu3YtISEhjBkzhh07dhAWFoa3tzfe3t6EhISw\nZ88etm3bxjPPPANAhw4dWLBggXv3SEREpIwr1uNRPT09GTVqFGvWrGHu3LkkJyfTu3dvbrnlFhYu\nXMj8+fNp2rQp/v5/PkPYarVis9mw2WxYrVbXsvT0dPfsiYiISDlR7GebT5s2jVOnTtG7d2+WLl3K\njTfeCMB9993HpEmTaNWqFXa73dXebrfj7++Pn5+fa7ndbicgIKDI1wq+ZCKBomRkeONhseS7LjCw\nYonrlQcak6tD4+x+BY3xmTN+BW4TFORX6vfGXXXdWfty617tMXZnbbPVLUiR4b1q1SqSk5MZPHgw\nvr6+WCwWhg4dytixYwkNDWXjxo00a9aM0NBQZs+eTWZmJg6Hg6SkJBo1akRYWBiJiYmEhoaSmJhI\ny5Yti+xUSkrJjs4zMjJwGka+69LSzpW4XlkXHOyvMbkKNM7FczlTdxY2xqmptgK3S021lfq9cVdd\nd9a+nLrXYozdWft6rFuaYC8yvLt06cKoUaPo27cv2dnZxMTEULNmTeLi4vDy8qJatWpMnDgRq9VK\n//79iYyMxOl0Eh0djY+PDxEREYwcOZLIyEh8fHyYNWtWiTspImWXpu4UKbkiw9vX15c5c+bkWb50\n6dI8y8LDwwkPD8+z/dy5cy+jiyJS1mnqTpGS0XzeIlKkyzm1LSJXnsJbRIqkU9si1xeFt4gUi05t\ni1w/9GxzERERk1F4i4iImIzCW0RExGQU3iIiIiaj8BYRETEZhbeIiIjJKLxFRERMRuEtIiJiMgpv\nERERk1F4i4iImIzCW0RExGQU3iIiIiaj8BYRETEZzSomUoZo3m2R8kHhLVKGaN5tkfJB4S1Sxmje\nbZGyT995i4iImIzCW0RExGQU3iIiIiaj8BYRETEZhbeIiIjJKLxFRERMRuEtIiJiMkXe552Tk8PY\nsWM5cOAAFouFuLg4fHx8GDVqFB4eHjRs2JDY2FgsFgsrVqxg+fLleHl5MWTIEDp16kRGRgbDhw8n\nNTUVq9XKtGnTCAoKuhr7JiIiUiYVGd7r1q3Dw8ODpUuXsmXLFl577TUAoqOjadWqFbGxsaxdu5YW\nLVqQkJDAypUrcTgcRERE0LZtW5YuXUrjxo2Jiopi9erVxMfHExMT4/YdExERKauKDO97772Xzp07\nA3D06FECAwPZuHEjrVq1AqBjx45s2LABDw8PwsLC8Pb2xtvbm5CQEPbs2cO2bdt45plnAOjQoQML\nFixw4+6IiIiUfcX6ztvT05NRo0YxZcoUunXrhmEYrnVWq5X09HRsNhv+/v65lttsNmw2G1arNVdb\nERERKb1iP9t82rRpnDp1ivDwcDIzM13LbTYbAQEB+Pn5YbfbXcvtdjv+/v65ltvtdgICAop8reBL\nJlUoSkaGNx4WS77rAgMrlrheeaAxuTqu9jifOeNX4LqgIL9S98ddda9E7YLWl8excFfdqz3G7qxt\ntroFKTK8V61aRXJyMoMHD8bX1xcPDw+aNWvGli1baN26NYmJibRp04bQ0FBmz55NZmYmDoeDpKQk\nGjVqRFhYGImJiYSGhpKYmEjLli2L7FRKSsmOzjMyMnBedDbgYmlp50pcr6wLDvbXmFwF12KcU1Nt\nha4rbX/cVfdyaxc2xuVtLNxV91qMsTtrX491SxPsRYZ3ly5dGDVqFH379iU7O5uYmBjq16/PuHHj\nyMrKokGDBnTp0gWLxUL//v2JjIzE6XQSHR2Nj48PERERjBw5ksjISHx8fJg1a1aJOykiIiJ/KjK8\nfX19mTNnTp7lCQkJeZaFh4cTHh6eZ/u5c+deRhdFRETkYnpIi4iIiMkU+4I1EblyMjMzOXz4YL7r\natcOwcfH5yr3SETMROEtcg0cPnyQEZ+Mx3rJhSr2lHSmd59IgwYNr1HPRMQMFN4i14g12B//mpWv\ndTdExIQU3iIF0KltEbleKbxFCqBT2yJyvVJ4ixRCp7ZF5HqkW8VERERMRuEtIiJiMgpvERERk1F4\ni4iImIzCW0RExGQU3iIiIiajW8XE9PQwFREpbxTeYnp6mIqIlDcKbykT9DAVESlP9J23iIiIySi8\nRURETEbhLSIiYjIKbxEREZNReIuIiJiMwltERMRkFN4iIiImo/AWERExGYW3iIiIySi8RURETKbQ\nx6NmZWUxZswYjh07RmZmJkOGDKF69eoMHjyYunXrAhAZGcmDDz7IihUrWL58OV5eXgwZMoROnTqR\nkZHB8OHDSU1NxWq1Mm3aNIKCgq7Gfsl1SBOIiIhcGYWG96effkpQUBAzZswgLS2NRx55hOeee46B\nAwcyYMAAV7uUlBQSEhJYuXIlDoeDiIgI2rZty9KlS2ncuDFRUVGsXr2a+Ph4YmJi3L5Tcn3SBCIi\nIldGoeHdpUsXHnjgAQCcTideXl7s2rWL/fv3s3btWkJCQhgzZgw7duwgLCwMb29vvL29CQkJYc+e\nPWzbto1nnnkGgA4dOrBgwQL375Fc1zSBiIjI5Ss0vCtVqgSAzWZj2LBhvPjiizgcDnr37s0tt9zC\nwoULmT9/Pk2bNsXf/8+jKavVis1mw2azYbVaXcvS09PduCsiIiLlQ5FTgh4/fpyoqCj69OlD165d\nSU9PdwX1fffdx6RJk2jVqhV2u921jd1ux9/fHz8/P9dyu91OQEBAsToVfMlp1aJkZHjjYbHkuy4w\nsGKJ65UH12JMzpzxK3BdUJBfqft0PdctqM313OerWfdK1C4rY+zO2mYbY3fWNlvdghQa3qdOnWLg\nwIHExsZy5513AvD0008TExNDaGgoGzdupFmzZoSGhjJ79mwyMzNxOBwkJSXRqFEjwsLCSExMJDQ0\nlMTERFq2bFmsTqWklOwIPSMjA6dh5LsuLe1cieuVdcHB/tdkTFJTbYWuK22frte6hY3z9drnq133\ncmuXpTF2Z22zjbE7a1+PdUsT7IWG98KFC0lPT+eNN97gjTfeAGDMmDFMnToVLy8vqlWrxsSJE7Fa\nrfTv35/IyEicTifR0dH4+PgQERHByJEjiYyMxMfHh1mzZpW4gyIiIpJboeE9duxYxo4dm2f50qVL\n8ywLDw8nPDw81zJfX1/mzp17mV0UERGRixX5nbeUL7oXW0Tk+qfwllx0L7aIyPVP4S156F5sEZHr\nm55tLiIiYjIKbxEREZNReIuIiJiMvvM2KV0VLiJSfim8TUpXhYuIlF8KbxPTVeEiIuWTvvMWEREx\nGYW3iIiIyei0Obr4S0REzEXhjS7+EhERc1F4/4+7Lv7SUb2IiFxpCm8301G9iIhcaQrvq0C3dImI\nyJWkq81FRERMRuEtIiJiMgpvERERk1F4i4iImIzCW0RExGQU3iIiIiaj8BYRETEZhbeIiIjJKLxF\nRERMRuEtIiJiMoU+HjUrK4sxY8Zw7NgxMjMzGTJkCA0aNGDUqFF4eHjQsGFDYmNjsVgsrFixguXL\nl+Pl5cWQIUPo1KkTGRkZDB8+nNTUVKxWK9OmTSMoKOhq7ZuIiEiZVGh4f/rppwQFBTFjxgzS0tJ4\n5JFHaNq0KdHR0bRq1YrY2FjWrl1LixYtSEhIYOXKlTgcDiIiImjbti1Lly6lcePGREVFsXr1auLj\n44mJibla+yYiIlImFXravEuXLjz//PMAOJ1OvLy82L17N61atQKgY8eObNy4kZ07dxIWFoa3tzd+\nfn6EhISwZ88etm3bRseOHQHo0KEDmzZtcvPuiIiIlH2FhnelSpWwWq3YbDaGDRvGCy+8gNPpdK23\nWq2kp6djs9nw9/fPtdxms2Gz2bBarbnaioiIyOUpckrQ48ePExUVRZ8+fXj44YeZMWOGa53NZiMg\nIAA/Pz/sdrtrud1ux9/fP9dyu91OQEBAsToVfMnc10XJyPDGw2LJd11gYMUi650541fguqAgvxL3\n52rUvty6Ba0vj2PhzrpXe5zNVvdK1C4rY+zO2mYbY3fWNlvdghQa3qdOnWLgwIHExsZy5513AtC0\naVO2bNlC69atSUxMpE2bNoSGhjJ79mwyMzNxOBwkJSXRqFEjwsLCSExMJDQ0lMTERFq2bFmsTqWk\nlOwIPSMjA6dh5LsuLe1ckfVSU22Fritpf65G7cupGxzsX+D68jYW7qx7LcbZbHUvt3ZZGmN31jbb\nGLuz9vVYtzTBXmh4L1y4kPT0dN544w3eeOMNAGJiYpgyZQpZWVk0aNCALl26YLFY6N+/P5GRkTid\nTqKjo/Hx8SEiIoKRI0cSGRmJj48Ps2bNKnEHRUREJLdCw3vs2LGMHTs2z/KEhIQ8y8LDwwkPD8+1\nzNfXl7lz515mF0VERORiekiLiIiIySi8RURETEbhLSIiYjIKbxEREZNReIuIiJiMwltERMRkFN4i\nIiImo/AWERExGYW3iIiIySi8RURETEbhLSIiYjIKbxEREZNReIuIiJiMwltERMRkFN4iIiImo/AW\nERExGYW3iIiIySi8RURETEbhLSIiYjIKbxEREZNReIuIiJiMwltERMRkFN4iIiImo/AWERExGYW3\niIiIySi8RURETKZY4f3jjz/Sr18/AHbv3k3Hjh3p168f/fr141//+hcAK1as4NFHH+Xxxx9n/fr1\nAGRkZDB06FD69OnDoEGDSE1Ndc9eiIiIlCNeRTV4++23+eSTT7BarQDs2rWLAQMGMGDAAFeblJQU\nEhISWLlyJQ6Hg4iICNq2bcvSpUtp3LgxUVFRrF69mvj4eGJiYty3NyIiIuVAkUfeISEhzJ8/H8Mw\nAPjpp59Yv349ffv2JSYmBrvdzo4dOwgLC8Pb2xs/Pz9CQkLYs2cP27Zto2PHjgB06NCBTZs2uXdv\nREREyoEiw/v+++/H09PT9XOLFi0YOXIkixcvpnbt2syfPx+73Y6/v7+rjdVqxWazYbPZXEfsVquV\n9PR0N+yCiIhI+VLkafNL3Xfffa6gvu+++5g0aRKtWrXCbre72lwIcz8/P9dyu91OQEBAsV4jONi/\n6EYXycjwxsNiyXddYGDFIuudOeNX4LqgIL8S9+dq1L7cugWtL49j4c66V3uczVb3StQuK2Psztpm\nG2N31jZb3YKUOLyffvppYmJiCA0NZePGjTRr1ozQ0FBmz55NZmYmDoeDpKQkGjVqRFhYGImJiYSG\nhpKYmEjLli2L9RopKSU7Qs/IyMD5v9P6l0pLO1dkvdRUW6HrStqfq1H7cuoGB/sXuL68jYU7616L\ncTZb3cutXZbG2J21zTbG7qx9PdYtTbAXO7wt/zuyjYuLIy4uDi8vL6pVq8bEiROxWq3079+fyMhI\nnE4n0dHR+Pj4EBERwciRI4mMjMTHx4dZs2aVuIMiIiKSW7HCu1atWixbtgyAJk2asHTp0jxtwsPD\nCQ8Pz7XM19eXuXPnXoFuioiIyAV6SIuIiIjJKLxFRERMRuEtIiJiMgpvERERk1F4i4iImIzCW0RE\nxGQU3iIiIiaj8BYRETEZhbeIiIjJKLxFRERMRuEtIiJiMgpvERERk1F4i4iImIzCW0RExGQU3iIi\nIiaj8BYRETEZhbeIiIjJKLxFRERMRuEtIiJiMgpvERERk1F4i4iImIzCW0RExGQU3iIiIiaj8BYR\nETEZhbeIiIjJKLxFRERMpljh/eOPP9KvXz8ADh48SEREBH369GHChAkYhgHAihUrePTRR3n88cdZ\nv349ABkZGQwdOpQ+ffowaNAgUlNT3bMXIiIi5UiR4f32228zduxYsrKyAJg6dSrR0dEsWbIEwzBY\nu3YtKSkpJCQksGzZMt59911mzZpFZmYmS5cupXHjxixZsoQePXoQHx/v9h0SEREp64oM75CQEObP\nn+86wt69ezetWrUCoGPHjmzcuJGdO3cSFhaGt7c3fn5+hISEsGfPHrZt20bHjh0B6NChA5s2bXLj\nroiIiJQPRYb3/fffj6enp+vnCyEOYLVaSU9Px2az4e/vn2u5zWbDZrNhtVpztRUREZHLU+IL1jw8\n/tzEZrMREBCAn58fdrvdtdxut+Pv759rud1uJyAg4Ap0WUREpHzzKukGTZs2ZcuWLbRu3ZrExETa\ntGlDaGgos2fPJjMzE4fDQVJSEo0aNSIsLIzExERCQ0NJTEykZcuWxXqN4GD/ohtdJCPDGw+LJd91\ngYEVi6x35oxfgeuCgvxK3J+rUfty6xa0vjyOhTvrXu1xNlvdK1G7rIyxO2ubbYzdWdtsdQtS7PC2\n/C8cR40axbhx48jKyqJBgwZ06dIFi8VC//79iYyMxOl0Eh0djY+PDxEREYwcOZLIyEh8fHyYNWtW\nsV4rJaVkp9czMjJwXnQ6/2JpaeeKrJeaait0XUn7czVqX07d4GD/AteXt7FwZ91rMc5mq3u5tcvS\nGLuzttnG2J21r8e6pQn2YoV3rVq1WLZsGQB169YlISEhT5vw8HDCw8NzLfP19WXu3Lkl7pSIiIgU\nTA9pERERMRmFt4iIiMkovEVERExG4S0iImIyCm8RERGTUXiLiIiYjMJbRETEZBTeIiIiJqPwFhER\nMRmFt4iIiMkovEVERExG4S0iImIyCm8RERGTUXiLiIiYjMJbRETEZBTeIiIiJuN1rTsgIiJypWVm\nZnL48ME8yw8dyrvMjBTeIiJS5hw+fJBhMz6hUmC1XMtPH/mZmx4sfd3r5UOBwltERAp1vQRWSVUK\nrIZflZtyLTublgycKnVNd30oKCmFt4iIFOp6CazrhTs+FJSUwltERIrkrsAy61H9tabwFhGRa0ZH\n9aWj8BYRkWvqejgNbTYKbxGRMkKnoMsPhbeIyFXmrpDVKejyQ+EtIleFjgr/5M6Q1Sno8qHU4d2z\nZ0/8/PwAqF27NoMHD2bUqFF4eHjQsGFDYmNjsVgsrFixguXLl+Pl5cWQIUPo1KnTleq7iJiI2Y4K\n3f1hQyErl6NU4e1wOABISEhwLXv22WeJjo6mVatWxMbGsnbtWlq0aEFCQgIrV67E4XAQERFB27Zt\n8fHxuTK9FxFTcUdg6RS0lEelCu9ffvmFc+fO8dRTT5Gdnc2LL77I7t27adWqFQAdO3Zkw4YNeHh4\nECJ5YdYAABCpSURBVBYWhre3N97e3oSEhLBnzx6aN29+RXdCRMovnYKW8qhU4V2xYkWeeuopwsPD\nOXDgAE8//XSu9VarlfT0dGw2G/7+/rmW22y2y+uxiLiNWb+XVshKeVOq8K5bty4hISGuf1euXJmf\nf/7Ztd5msxEQEICfnx92u9213G63ExAQcJldFhF30aliEXMoVXivXLmSPXv2EBsbS3JyMna7nXbt\n2rFlyxZat25NYmIibdq0ITQ0lNmzZ5OZmYnD4SApKYmGDRsWWT842L/INhfLyPDGw2LJd11gYMUi\n650541fguqAgvxL352rUvty6Ba0vj2PhzrpXe5wvt+6ZM34lPootbn8L61tBilM7KMg9dd3VX3fW\nNtsYg/nGwp2/FyVRqvB+7LHHGD16NH369AFg6tSpVK5cmXHjxpGVlUWDBg3o0qULFouF/v37ExkZ\nidPpJDo6ulgXq6WkpJeoPxkZGTgNI991aWnniqyXmlrwqfzUVFuJ+3M1al9O3eBg/wLXl7excGdd\nd45zaU5vF6duYf26nLruqh0c7O+2PpttLNxV151jfKGdO2qbqW5pgr1U4e3l5cWMGTPyLL/46vML\nwsPDCQ8PL83LiFwV7vyeV1dCi4g76CEtUu65Mwh1JbSIuIPC+//bu9ugKKv3D+DfBRYckwFEnWZS\nSM3EtBQks8hBDDMFJIQETXzMzBKbQnPVCSllpDFfKDjOaOYLdNRKZBp0fIIpLWEljewBNQ15kFBQ\nnlxWFriv3wv+7t8VWNRY1xu+n1fs3mfPubh2Z6+9z717DhFsWwhZZImos7F4U6dT68+NiIjUgsW7\nm7JlgbXVVDE/FBARtWDx7qZs/YUnW0wV80taREQtWLy7MTVei1VjzEREna3bFe+2pl4f52lXThUT\nEdG9ul3xbmvqtTOmXfl7XiIielS6XfEGWk+9dsa0K3/PS0REj0q3LN62wiJLRESPgoO9AyAiIqIH\nw+JNRESkMizeREREKsPiTUREpDIs3kRERCrD4k1ERKQyLN5EREQqw+JNRESkMizeREREKsPiTURE\npDIs3kRERCrD4k1ERKQyLN5EREQqw+JNRESkMizeREREKsPiTUREpDIs3kRERCrjZOsBFEVBYmIi\nLl68CK1Wi6SkJHh5edl6WCIioi7L5mfex48fR2NjI/bu3Ytly5YhOTnZ1kMSERF1aTYv3mfPnsW4\nceMAACNHjsQff/xh6yGJiIi6NJtPm9+6dQu9evUy33Z0dISiKHBwsP65ITf3FGpqalrd7+CggaKI\nxX0mkwmXmhU0O1r2aXByxIsuLq36qK+5bnHbWHcThoq6Vu0MFXUoLi5qM77Bg4d02G9n9a22ftUY\nM3Nh+37VGDNzYft+1RizLXNxvzQiIh03e3jJyckYOXIkJk+eDAAIDAzEjz/+aMshiYiIujSbT5v7\n+fnhxIkTAID8/HwMHTrU1kMSERF1aTY/8xYRJCYm4sKFCwCA9evXY+DAgbYckoiIqEuzefEmIiKi\nzsVFWoiIiFSGxZuIiEhlWLyJiIhUxi7FW1EUJCQkICYmBrGxsSguLrY4np2djaioKMTExODbb7+1\nR4hdQkd5zszMxPTp0zFjxgysWbMG/PrDg+sox3d8+umn2Lhx4yOOrmvoKMfnzp3D22+/jZkzZ+Kj\njz6CyWSyU6Tq1lGejx07hsjISERFRWHPnj12ilL9fvvtN8TGxra6/4HrntjBkSNHRKfTiYhIfn6+\nLF682HzMZDLJxIkTpba2Vkwmk0RGRkplZaU9wlQ9a3k2Go0SHBwst2/fFhGRjz/+WLKysuwSp5pZ\ny/Ede/bskejoaNm4ceOjDq9LsJZjRVEkPDxciouLRURk3759cvnyZbvEqXYdvZaDgoKkpqbG4j2a\nHsy2bdskNDRUoqOjLe5/mLpnlzNva0umXr58GV5eXnB1dYVWq8Xo0aORl5dnjzBVz1qeXVxcsG/f\nPrj83wp0TU1N6NGjh13iVLOOlv89e/Yszp07h+joaM5sPCRrOS4sLIS7uzt27tyJ2NhY1NbWYtCg\nQfYKVdU6ei1rtVrU1taioaEBIgKNRmOPMFXN29sbqamprd4LHqbu2aV4t7dk6p1jrq6u5mNPPPEE\n6upaLzlHHbOWZ41Gg969ewMA0tLSYDQa8corr9glTjWzluPr169jy5YtSEhIYOH+D6zluKqqCr/+\n+itmzZqFnTt3IicnB7m5ufYKVdWs5RkA5s2bh8jISISGhiIoKMiiLd2f119/HY6Ojq3uf5i6Z5fi\n3atXLxgMBvPtu9c6d3V1tThmMBjg5ub2yGPsCqzl+c7tL774Ajk5OUhJSbFHiKpnLcdHjhxBVVUV\nFi5ciO3btyMzMxMZGRn2ClW1rOXY3d0dXl5eGDRoEJycnDBu3DhufvSQrOW5rKwMu3fvRnZ2NrKz\ns3Hjxg0cPnzYXqF2OQ9T9+xSvK0tmTpo0CAUFRWhpqYGJpMJeXl5GDVqlD3CVL2OlqZNSEiAyWTC\nli1bzNPn9GCs5Tg2Nhbp6elIS0vDu+++i9DQULz55pv2ClW1rOV4wIABqK+vN3+56syZMxgy5OE3\ne+jOrOW5oaEBDg4OcHZ2hoODA3r37s0Z0U70MHXP5ruKtWXixIn4+eefERMTA6BlydTMzEzU19dj\n+vTp0Ol0WLBgARRFQVRUFPr162ePMFXPWp5HjBiB/fv3w9/fH7NnzwYAzJkzB8HBwfYMWXU6ei3f\njdcIH05HOU5KSkJ8fDxEBH5+fggMDLRzxOrUUZ4jIiIQExMDFxcXeHt7IyIiws4Rq9ed94L/Uve4\nPCoREZHKcJEWIiIilWHxJiIiUhkWbyIiIpVh8SYiIlIZFm8iIiKVYfEmIiJSGRZvov8gJSUFqamp\nNuu/rKwMb7zxBiIjIy1WYLKXW7duYdq0aYiIiMCVK1dsMsbhw4excuXK+25/8+ZNTJgwocN2vr6+\nVo+XlJRg9erV9z0ukT3ZZZEWoq7C1guvnD59GsOHD39sthMtKCiAs7Mz9u7da+9QOl1ZWVm7W7oS\nPW5YvKnb0Ov12Lp1KwCgvLwcL7zwAtatWwdnZ2dzm+TkZPTr1w/z588HACxduhRhYWHw9vbG2rVr\nYTQacfPmTcybN6/Vnrw+Pj44f/48ACA9PR15eXlYv349zp07h+TkZNy+fRseHh747LPP0L9/f4vH\nFhYWIiEhATU1NejZsydWr14NrVaLTZs2ob6+HomJiUhMTDS3T0lJgUajwZIlSwAAEyZMQFpaGvR6\nPU6ePIna2lqUlJQgICDAvFf7l19+iePHj8PJyQnR0dGYPXt2m+M+//zz0Ol0cHV1xZ9//ony8nIs\nWbIE48ePx6pVq1BZWYn3338fW7ZsQVJSEnJzc6HRaDB16lQsXLgQer0eqampSEtLAwDodDq89NJL\nGDNmDD744AM8++yzKCgogKenJzZt2gQ3Nzd8//332Lp1K3r27Alvb2/zcr3t5a6goACrV6+GiGDE\niBFtPt9Xr17F8uXLYTAY8Nxzz5k3h7l27RpWrVqFW7duoaKiAiEhIYiPj8e6detQWlqKtWvXYtWq\nVVizZg0uXbqEyspKDBw4EKmpqVxGmB4fnb9jKdHjKTc3V0aNGiVFRUWiKIosXbpUdu7cadHmr7/+\nkmnTpomISF1dnbz66qvS0NAgSUlJkpOTIyIixcXF4uvrKyIimzdvlpSUFBERGTp0qLmf9PR00el0\nYjKZJCwsTP79918RETlx4oTMnTu3VWyRkZFy7NgxEWnZSzkoKEgaGhrM/dwrJSXFPK5Iy17LV69e\nlf3798v48ePFYDCI0WiUwMBAuXDhghw6dEhmzJghJpNJDAaDhIeHS0VFRbvjrlixQuLi4kRE5MKF\nCzJmzBgREdHr9TJr1iwREdm1a5csWbJEFEURo9EoUVFR8sMPP1i0ERHR6XRy4MABKSkpER8fHyko\nKBARkbi4OElLS5Py8nIJCAiQiooKaW5ulsWLF3eYu9DQUPnpp59ERGTHjh0SFBTUKkeLFi2Sffv2\niYjI4cOHzc/Pjh075MCBAyIiUltbK35+flJVVWURd15ennz++eci0rJn+KxZs+TIkSOtxiCyF555\nU7fy8ssvw8vLCwAQHh6Ob775BnPnzjUfHzZsGEwmE4qLi3H27FkEBQXB2dkZOp0OJ06cwLZt23D+\n/HkYjcb7Gu/KlSsoKSnBe++9Z77v3mvXBoMBJSUl5nXlR44cCTc3NxQWFj7QVqJ32vr6+qJnz54A\nWjbuqKmpwS+//IIpU6ZAq9VCq9UiIyPD6rgajQYBAQEAgCFDhqCmpsZiDKBlJiMiIgIajQY9evRA\nWFgYcnJyrF5/9vT0hI+Pj7nf6upq5Ofnw9fXF3369AEATJs2DVlZWe3mrqqqCtevXzfHFxUVhV27\ndrUaS6/Xmy83TJo0ybyF5fz585Gbm4uvv/4aFy9eRGNjI4xGo8X/5u/vD3d3d+zevRv//PMPioqK\nUF9ff9/PBZGtsXhTt+Lk9P8veUVR2txbNywsDAcPHkR+fj4WLlwIAPjwww/h7u6OoKAgTJkyBYcO\nHbI6TmNjIwCgubkZAwYMMG8FqigKKioqLNqKSKsiLSJQFMXqNfW7H9PU1GT++96pXRGBk5OTRfvS\n0lK4ubm1OW5zczMAmC8ntBfDvXErioKmpqZW7e/k4u4+7+5Xo9FY9HNnG0pruWur/b3u7ffOc52c\nnIzS0lKEhYUhODgYOTk5rfKQlZWFlJQUzJkzB5GRkaiurm5zDCJ74bfNqVvR6/WoqKiAoijIyMho\ncweqqVOn4tChQygqKoK/vz8A4NSpU4iLi8OECRNw+vRpAC3F5G4eHh74+++/ISLIzs4G0LLV350z\nXwDYv38/li1bZvG4Xr16YcCAATh27BiAlu0YKysrMWTIkHbPvD08PHDp0iUALdeFKyoqrBb6F198\nEUePHkVTUxOMRiPeeecd3Lhxo91x78fYsWORkZEBRVFgNBqRmZmJsWPHwsPDAyUlJTCZTKiursaZ\nM2es9jN69Gjk5+ejvLwcIoKDBw8CaD93Hh4e6N+/P7KysgC07MzUloCAAKSnpwMATp48aZ49OHXq\nFBYsWIBJkyahrKwM165dQ3NzMxwdHc0fXHJycjB58mRERETA09MTeXl5Fh+QiOyNZ97UrfTr1w/L\nli0zT7u+9dZbrdo8+eST6N27t8V+unFxcZg5cyb69OkDf39/DB48GKWlpRYFMz4+HosWLULfvn3h\n5+eH6upqODs7Y9OmTUhKSkJDQwNcXV2RnJzcaswNGzZgzZo12Lx5M1xcXJCamgonJ6d2C3JISAiO\nHj2KkJAQDB8+3PyFrLbaazQaBAcH4/fff0dERAREBHPnzsXTTz/d5rhardb8uLv7uPfv6OhoFBYW\nIjw8HI2NjQgPDzdPwQcGBiIkJARPPfWU+QOQRqNpMz5PT08kJiZi/vz56NGjh3kfaWu527BhA1au\nXInU1FSMHj26zX4TEhKwfPlyfPfdd/Dx8TFPyy9atAiffPIJPD098cwzz2Ds2LEoLS3FsGHDUFdX\nhxUrVmDBggWIj4/H0aNH0bdvX7z22mu4evVqm88FkT1wS1DqNvR6PbZv346vvvrK3qEQEf0nnDan\nbqO9Mz8iIrXhmTcREZHK8MybiIhIZVi8iYiIVIbFm4iISGVYvImIiFSGxZuIiEhlWLyJiIhU5n/v\nivCOaU6Q8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116d84ef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=30\n",
    "stims=40\n",
    "effect_size=0.\n",
    "confound_size=.5\n",
    "runs=100000\n",
    "\n",
    "plot_hist(*run_stim(n=n, stims=stims, effect_size=effect_size, confound_size=confound_size, runs=runs), \n",
    "          n=n, effect_size=effect_size, confound_size=confound_size, runs=runs,\n",
    "          ax=plt.gca(), pltx=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}