Square-root LASSO with highdim knockoff data generating mechanism

knockoff_data_generating_with_sqrt_lasso.ipynb

lasso_cv.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division, print_function\n",
    "\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from IPython.display import HTML\n",
    "\n",
    "# this code below is available at https://github.com/jonathan-taylor/selective-inference\n",
    "from selection.algorithms.lasso import lasso\n",
    "from selection.algorithms.sqrt_lasso import choose_lambda\n",
    "\n",
    "%R library(glmnet);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data generating mechanism\n",
    "\n",
    "We use the same data generating mechanism as in [Barber and Candes (2016)](http://arxiv.org/abs/1602.03574)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.3362433974107397e-15"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def cov(p, rho=0.25):\n",
    "    idx = np.arange(p)\n",
    "    return rho**np.fabs(np.subtract.outer(idx, idx))\n",
    "\n",
    "def sqrt_cov(p, rho=0.25):\n",
    "    idx = np.arange(p)\n",
    "    C = rho**np.fabs(np.subtract.outer(idx, idx))\n",
    "    return np.linalg.cholesky(C)\n",
    "\n",
    "# Testing we have the square-root correct\n",
    "p = 2500\n",
    "A = cov(p, rho=0.3)\n",
    "B = sqrt_cov(p, rho=0.3)\n",
    "np.linalg.norm(B.dot(B.T) - A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cholesky_factors = {}\n",
    "for rho in [0, 0.25, 0.5, 0.75]:\n",
    "    cholesky_factors[(rho, 2500)] = sqrt_cov(2500, rho=rho)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def instance(n=2000, k=30, p=2500, signal=4.5, rho=0.25, sigma=1, q=0.2): \n",
    "                                                                       \n",
    "    if (rho, p) in cholesky_factors.keys():\n",
    "        _sqrt_cov = cholesky_factors[(rho, p)]\n",
    "    else:\n",
    "        cholesky_factors[(rho, p)] = sqrt_cov(p, rho=rho)\n",
    "        _sqrt_cov = cholesky_factors[(rho, p)]\n",
    "    X = np.random.standard_normal((n, p)).dot(_sqrt_cov.T)\n",
    "\n",
    "    X /= (np.sqrt((X**2).sum(0))) # like normc\n",
    "    beta = np.zeros(p)\n",
    "    beta[:k] = signal * (2 * np.random.binomial(1, 0.5, size=(k,)) - 1) \n",
    "    np.random.shuffle(beta)\n",
    "    Y = (X.dot(beta) + np.random.standard_normal(n)) * sigma\n",
    "    true_active = np.nonzero(beta != 0)[0]\n",
    "\n",
    "    %R -i X,Y Y = as.matrix(Y)\n",
    "    %R X = as.matrix(X)\n",
    "    %R G_CV = cv.glmnet(X, Y, standardize=FALSE, intercept=FALSE)\n",
    "    lam_1SE = %R G_CV$lambda.1se\n",
    "    lam_minCV = %R G_CV$lambda.min\n",
    "    \n",
    "    lam_1SE *= n\n",
    "    lam_minCV *= n\n",
    "    \n",
    "    P0 = {}\n",
    "    PA = {}\n",
    "    full_model_FDP = {}\n",
    "    full_model_power = {}\n",
    "    directional_FDP = {}\n",
    "    active_size = {}\n",
    "\n",
    "    for method, lam in zip(['minCV', '1SE'], [lam_minCV, lam_1SE]):\n",
    "        L = lasso.gaussian(X, Y, float(lam), sigma=sigma)\n",
    "        L.fit() \n",
    "        S = L.summary('onesided')\n",
    "        if set(true_active).issubset(L.active):        \n",
    "            P0[method] = [p for v, p in zip(S['variable'], S['pval']) if v not in true_active]\n",
    "            PA[method] = [p for v, p in zip(S['variable'], S['pval']) if v in true_active]\n",
    "        else:\n",
    "            P0[method] = []\n",
    "            PA[method] = []\n",
    "\n",
    "        selected = sm.stats.multipletests(S['pval'], q, 'fdr_bh')[0]\n",
    "        type1_errors = selected * np.array([v not in true_active for v in S['variable']])\n",
    "        sign_errors = selected * np.array([(s != np.sign(beta[v])) * (v in true_active) for v, s in \n",
    "                                            zip(S['variable'], L.active_signs)])\n",
    "\n",
    "        full_model_FDP[method] = [np.sum(type1_errors) / max(np.sum(selected), 1)]\n",
    "        directional_FDP[method] = [np.sum(type1_errors + sign_errors) / max(np.sum(selected), 1)]\n",
    "        full_model_power[method] = [np.sum(selected * np.array([v in true_active for v in S['variable']])) / len(true_active)]\n",
    "        active_size[method] = [len(L.active)]\n",
    "   \n",
    "    return P0, PA, full_model_FDP, directional_FDP, full_model_power, active_size\n",
    "\n",
    "P0, PA, full_model_FDP, directional_FDP, full_model_power, active_size = instance(n=200, p=400, k=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_table(nsim, rho, q, n=2000, p=2500, k=30):\n",
    "\n",
    "    P0, PA, full_model_FDP, directional_FDP, full_model_power, active_size = {}, {}, {}, {}, {}, {}\n",
    "\n",
    "    rows = r\"\"\"\n",
    "    <tr><td>Method</td><td>Model FDP</td><td>Directional FDP</td><td>Model power</td><td>Active set</td></tr>\n",
    "    \"\"\"\n",
    "\n",
    "    for _ in range(nsim):\n",
    "        _P0, _PA, _full_model_FDP, _directional_FDP, _full_model_power, _active_size = instance(n=n, p=p, k=k, rho=rho, q=q)\n",
    "        for method in _P0.keys():\n",
    "            for d1, d2 in zip([ P0, PA, full_model_FDP, directional_FDP, full_model_power, active_size],\n",
    "                              [_P0, _PA, _full_model_FDP, _directional_FDP, _full_model_power, _active_size]):\n",
    "                d1.setdefault(method, []).extend(d2[method])\n",
    "\n",
    "    for method in P0.keys():\n",
    "        rows += \"\"\"\n",
    "        <tr><td>%s</td><td>%0.3f</td><td>%0.3f</td><td>%0.3f</td><td>%0.1f</td></tr>\n",
    "        \"\"\" % (method, np.mean(full_model_FDP[method]), np.mean(directional_FDP[method]), np.mean(full_model_power[method]), np.mean(active_size[method]))\n",
    "\n",
    "    table = \"\"\"\n",
    "    <table>\n",
    "    %s\n",
    "    </table>\"\"\" % rows\n",
    "\n",
    "    return table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\rho=0.25$ using `cv.glmnet` two choices of $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <table>\n",
       "    \n",
       "    <tr><td>Method</td><td>Model FDP</td><td>Directional FDP</td><td>Model power</td><td>Active set</td></tr>\n",
       "    \n",
       "        <tr><td>minCV</td><td>0.114</td><td>0.114</td><td>0.106</td><td>152.4</td></tr>\n",
       "        \n",
       "        <tr><td>1SE</td><td>0.085</td><td>0.085</td><td>0.619</td><td>48.2</td></tr>\n",
       "        \n",
       "    </table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTML(make_table(100, 0.25, 0.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1269e8e50>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAJnCAYAAAAXwG6AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xu8lnWd//vXl5McRQE5CIIKIioqppiKyBJJ7WDm6ELN\nkpimnJ3tcTe/ZpqpqWzGpixrZqrdrimdRsusZWba5Dix4AJRUctBQE7KSZSTgAoiymF99x9rQQtY\ni3Wfr/vwej4e6yHrvr/3dX14LMU338PnCjFGJEmSVDk6pV2AJEmSsmOAkyRJqjAGOEmSpApjgJMk\nSaowBjhJkqQKY4CTJEmqMKkHuBDCnSGEjSGEBYcZ850QwgshhPkhhHGlrE+SJKncpB7ggP8ALmvv\nzRDCe4GRMcaTgJuAH5SqMEmSpHKUeoCLMc4FXjvMkCuBu1vGPgX0DSEMKkVtkiRJ5Sj1AJeBocDa\nVt+/0vKaJElSTaqEACdJkqRWuqRdQAZeAY5r9f2wltcOEULwwa6SJKlixBhDLp8rlxm40PLVloeA\nGwFCCOcBr8cYN7Z3oRijXxX69eUvfzn1GvzyZ1eLX/78KvvLn1+WX/feS5w6NbX7r1mzZv+v85F6\ngAsh3As8AYwOIbwUQpgeQrgphPBJgBjj74BVIYQXgR8Cn0qxXEmSVMlWroQTTkjl1kmScM455/DS\nSy/lfa3Ul1BjjB/OYMynS1GLJEmqcqtWwbnnlvy2SZJQX19PQ0MDw4cPz/t6qc/ASfvU1dWlXYJy\n5M+usvnzq2z+/LKUwgxc6/BWqJ9XyHcNtpyEEGI1/X4kSVKBHX88NDbCyJElud2cOXO4+uqr2wxv\nIQRijocYDHCSJKk27N4NvXvDm29C164lueXKlSt5+eWXueiiiw55zwDXwgAnSZLatWIFXHIJrF6d\ndiVAfgHOPXCSJKk2rFoFJ56YdhUFYYCTJEm1IcUWIoVmgJMkSbVh5cqizsAlScI3vvGNol2/NQOc\nJEmqDUVcQt3XKuTcEvWYM8BJkqTaUKQl1GL0eeuIAU6SJNWGIszApRHewDYikiSpFrzxBgwdCtu3\nQ8ipc8ch9u7dy4QJE/j617+eU3jLp41I6s9ClSRJKrpVq5qXTwsU3gA6d+7M448/TufOnQt2zUy5\nhCpJkqpfkQ4wpBHewAAnSZJqQRX1gAMDnCRJqgUFmIFbsWIF5bLX3gAnSZKqX55NfJMk4bzzzmPl\nypUFLCp3BjhJklT98lhCbd0qZOTIkQUuLDe2EZEkSdWtqQl69YItW6Bnz6w+Wsw+b/m0EXEGTpIk\nVbf166Fv36zD29y5c1Np0psJ+8BJkqTqluMBhhNOOIFf//rXXHjhhUUoKj8GOEmSVN1yPMAwdOhQ\nhg4dWoSC8ucSqiRJqm5V1gMODHCSJKnaFekpDGkywEmSpOqWwQxckiTceuutpamnAAxwkiSpunWw\nB25fq5ByO2l6OPaBkyRJ1evtt+Goo2DHDmjjwfPF7PPWEfvASZIktWX1ajjuuLILb/kywEmSpOrV\nzgGGpqYmbr311ooMb2AfOEmSVM3aOcDQqVMnZs6cSadOlTmXVZlVS5IkZeIwBxgqNbyBAU6SJFWz\nKuwBBwY4SZJUzVqWUJcuXUpTU1Pa1RSMAU6SJFWnGGHVKpJ165g4cSLLly9Pu6KCsQ+cJEmqTlu2\nkIwYQX2PHmV52jSfPnCeQpUkSVUpaWig/p13aPjtb8suvOXLGThJklR1nnjiCa68/HIazjqLutmz\n0y6nTc7ASZIktTJ69Gh+c8MNXNC7d9qlFIUBTpIkVZ0BAwYwoKmpKluIgAFOkiRVupkz4QtfOPT1\nZcvg6qtLX08JuAdOkiRVtj//cxg4EK688sDXO3eGs89u80H25SCfPXAGOEmSVLliJBk8mP96//v5\n5l13pV1NVvIJcDbylSRJFSu5917qN2/m/R/9aNqllJQzcJIkqSIlSUL9FVfQcN551P3+92mXkzVn\n4CRJUk1JkoT6+noa3vUu6j7ykbTLKTkDnCRJqigxRm6//XYafvEL6p5/Hi65JO2SSs4lVEmSVHFi\njIT58+G665rbhVQgl1AlSVJNCSFAY2NNzr6BAU6SJFWqGTNgypS0q0iFAU6SJJW1hQsXsnfv3gNf\n3LULnngC6upSqSltBjhJklS2Zs2axeTJk1m8ePGBb8ybB6NHQ79+6RSWMgOcJEkqS7NmzWLq1Kk0\nNDRw+umnH/hmDe9/AwOcJEkqQ63DW11by6Q1vP8NbCMiSZLKzLx587jiiivaD2/bt8OQIbBpE/Ts\nWfL6CiWfNiJdCl2MJElSPk455RQefvhhzjvvvLYHzJkD48dXdHjLl0uokiSprPTt27f98AY1v/8N\nDHCSJKnS1Pj+N3APnCRJqiSbNjW3D9m8GbpU9k4wH6UlSZIqUpIk3HzzzZl/YOZMuOiiig9v+TLA\nSZKkVCRJQn19PfX19Zl/yP1vgKdQJUlSoa1ZA88+e9ghyaJF1N9xBw2f/Sx1r70Gv/51Ztd+9FG4\n5ZYCFFnZ3AMnSZIK65prYP16GDiwzbeTV1+l/plnaBg/nrpjjsnu2v37w49+BCGnrWNlxT5wkiSp\nPOzdC7NmwcKFcOyxh7wdY+S711xDw6OPtt2kVxlxBk6SJBXOH/8IH/0oHPzw+VZijIQqmEHLl6dQ\nJUlSecjgkIHhLX8GOEmSVDg22S0JA5wkSSqMd96BJ5+ESZP2vzR//nx27dqVYlHVyQAnSZIK48kn\n4dRT4aijgOY+b+95z3tYfJj9cMqNAU6SJBVGq/1v+5r0NjQ0MG7cuJQLqz6eQpUkSYVx/vnw1a+S\ndOq0P7zZKqR99oGTJEnp2rYNFi3imW7dqL/qKsNbkRngJElS/mbPhne/m9Pe9S5+97vfMX78+LQr\nqmrugZMkSflr2f/Ws2dPw1sJGOAkSVL+7P9WUh5ikCRJ+dmwobl9yKuvQufOaVdTMXyUliRJKrkk\nSZg2bRrMnNncvNfwVjIGOEmSlLV9fd6mT5+e0fNPVVgGOEmSlJXWTXrrJk1y/1sKDHCSJCljB4S3\nujpYsQL27IGTT067tJpigJMkSRm76667DmzSu2/5NOS0F1858hSqJEnK3dSp8P73w7RpaVdScfI5\nhWqAkyRJuWlqgoEDYf58GDYs7Woqjs9ClSRJhfGf/wm33ZbZ2L17mwOc4a3knIGTJEl/cumlUF8P\ndXX8YeFCTh01ip49erQ/vl8/6N+/dPVVEZdQWxjgJEnKwzvvwIABsHYtyfz51NfX88gjj3DOOeek\nXVlV8kkMkiQpf08+Caeeuj+8NTQ0GN7KlAFOkiQ1mzGD5KSTDuzzprLkEqokSQLgf884g0tfeomG\nBx80vJWAe+BaGOAkScrRtm28c+yxLJ4xg7POOy/tamqCe+AkSVJ+Zs/miPPOM7xVCAOcJEn60yOx\nVBEMcJIkCWbMMMBVEAOcJEk1KEkS6uvrm7/ZsAFeeQXOPjvdopQxH6UlSVKN2RfeGhoaml+YORPq\n6qBz51TrUuacgZMkqYa0Dm/7W4W4/63iGOAkSaoRbYa3GN3/VoEMcJIk1Yj77rvv0CcsrFgBu3fD\nmDGp1aXs2chXkqRa9sMfwuOPw913p11JzbGRryRJyo373yqSM3CSJNWqpiYYOBDmz4dhw9KupuY4\nAydJkg7w1FNPsW3btsMPeu456N/f8FaBDHCSJFWZJEn4wAc+wJIlSw4/0OXTimWAkySpirRuFfLu\nd7/78IMbG2HKlNIUpoJyD5wkSVWizT5v7dm1CwYMgNWroV+/UpSng1T8HrgQwuUhhKUhhOUhhM+1\n8f6RIYSHQgjzQwgLQwgfS6FMSZLK1sKFCzMPbwDz5sHo0Ya3CpX6DFwIoROwHLgEWAc8A1wXY1za\naszfA0fGGP8+hDAAWAYMijHuOehazsBJkmrSnj17WLJkCaeffnpmH/jyl+Htt+H224tbmNpV6TNw\n5wIvxBjXxBh3A/cBVx40JgJ9Wn7dB9hycHiTJKmWdenSJfPwBh5gqHDlEOCGAmtbff9yy2utfQ84\nNYSwDngOuKVEtUmSVH22b2/u/XbhhWlXohx1SbuADF0G/G+McXIIYSTw+xDCGTHGNw8eeOutt+7/\ndV1dXWb7ACRJqjAxRkLIafUN5syB8eOhZ8/CFqXDSpKEJEkKcq1y2AN3HnBrjPHylu//Dogxxttb\njfkt8LUY4+Mt3zcCn4sx/uGga7kHTpJU9ZIk4Zvf/Ca//e1vcwtxf/3XzYcX/uEfCl+cMlbpe+Ce\nAUaFEEaEELoB1wEPHTRmDTAFIIQwCBgNrCxplZIklYF9rUL+5m/+JvcZOPe/VbzUl1BjjHtDCJ8G\n/ofmQHlnjHFJCOGm5rfjvwO3AT8JISxo+djfxhi3plSyJEmpyKrPW3s2bYI1a5qXUFWxUl9CLSSX\nUCVJ1aog4Q3gvvvg3nvhoYMXu1Rqlb6EKkmSOvDwww/nH97A5dMq4QycJEm15MQTm2ffxo5Nu5Ka\n5wycJEnq2MqV8NZbcNppaVeiPBngJEmqFfuWT3M9vaqyYYCTJKnMPPHEE2zevLnwF3b/W9UwwEmS\nVEaSJOHKK69k+fLlhb1wUxPMnGmAqxIGOEmSykTrViEXXHBBYS++aBH07QsjRhT2ukqFAU6SpDJQ\nsD5v7Zkxw9m3KmKAkyQpZUuWLClueAP3v1UZ+8BJkpSypqYmli1bximnnFKcG+zeDQMGwIoVzf9U\nWbAPnCRJFaxTp07FC28ATz8NI0ca3qqIAU6SpGrn/req0yXtAiRJqngLFjQHpL17MxoeYySUspnu\nm2/CI4+U7n4qOvfASZKUr699Ddauhdtu63BoMncu//jNb9L44IOlC3GdOsFRR5XmXspYPnvgnIGT\nJClfjY1wyy3Qr99hhyVJQv3HP05DQwOhf/8SFadq5AycJEn52LkTBg6EV16BI49sd1jR+7yp4ngK\nVZKktDzxBJx+uuFNJWWAkyQpHxk0yG1sbDS8qaBcQpUkKR/vfjfcfjsYzpSlfJZQDXCSJOXq9ddh\n+HB49VU44oi0q1GFcQ+cJElpSBI4/3zDm0rOACdJUq7a2P82d+5c1q9fn1JBqhUGOEmScnVQgEuS\nhKuuuooVK1akWJRqgQFOkqRcrFsHGzfCuHHAga1CLrzwwpSLU7UzwEmSlIvGRrj4Yujc2T5vKjkD\nnCRJuWhZPn3xxRcNbyo524hIkpStGJvbhzQ2Ek86iRUrVjBq1Ki0q1KFsY2IJEml9MILEAKcdBIh\nBMObSs4AJ0lStmbMaD59GnKaPJHyZoCTJCkLe/fuzej5p1IxGeAkScpQkiRcOGECe2fNgsmT0y5H\nNcwAJ0lSBva1Cvna9Ol0HjwYjj027ZJUwwxwkiR14IA+b6+/DlOmpF2SalyXtAuQJFWR7dvhBz+A\nvXvTrqRgkpUrqb/3Xho+/GHq5s2Dn/0Mbrst7bJU4wxwkqTCefhhuOceeN/70q6kYOYtX07DBz9I\nXb9+8Prr8KEPwXvek3ZZqnE28pUkFc7HPw5nnQWf/nTalUhlz0a+kqT0xfin/miSisoAJ0kqjBUr\nYM8eGDMm7UqkqmeAkyQVxr7mthX8dII5c+awatWqtMuQOmSAkyQVRoU/nSBJEq6++mrWrl2bdilS\nhzzEIEnKX1MTDBwI8+fDsGFpV5O1A/q81dWlXY5qhIcYJEnpeu45GDDA8CaViAFOkpS/Cl0+XbNm\nDVOnTjW8qeLYyFeSlL/GRvjEJ9KuImsjRozgmWeeYcSIEWmXImXFPXCSpPzs2tW8fLp6NfTrl3Y1\nUsVwD5wkKT3z5sHJJxvepBIywEmS8lNB+9/27NmTdglSQRjgJEn5qZAAlyQJ5557riFOVcFDDJKk\n3G3f3txCZMKEtCs5rNatQrp08X99qnzOwEmScjdnDowfDz17pl1Ju+zzpmpkgJMk5a7Ml08Nb6pW\nBjhJUu7KPMAtWLDA8KaqZB84SVJuNm2C0aNh82ZwX5mUNfvASZJKb+ZMmDTJ8CalwP/qJKkWxdj8\nlY8yXz6VqpkBTpJq0ZQpzTNoIafVm2Zdu8KCBYWrKU+zZ89m0KBBjBkzJu1SpKJzCVWSas22bfDU\nU/DWW9DUlPvXO+80P0KrDCRJwjXXXMPGjRvTLkUqCQOcJNWa2bPh3e+GHj3SrqQgWrcKmTRpUtrl\nSCVhgJOkWjNjRvMSahWwz5tqlQFOkmpNlRw+WLduHVOnTjW8qSbZB06SasmGDXDKKc292zp3Trua\nvK1bt45jjz027TKknNgHTpKUmX2926ogvAGGN9UsA5wk1ZIq2v8m1TIDnCTVihgrev/brl270i5B\nKhsGOEmqFStWwO7dUIGNbpMk4ayzzjLESS18EoMk1Yp9s2/5PH0hBa1bhXTr1i3tcqSy4AycJNWK\nCtz/Zp83qW22EZGkWtDUBAMHwvz5MGxY2tVkxPCmamcbEUnS4T33HPTvXzHhDWDlypWGN6kdzsBJ\nUi244w5YuRK+//20K5HUwhk4SdLhVeD+N0ntcwZOkqrdrl0wYACsXg39+qVdjaQWzsBJkto3bx6M\nHl3W4S1JEp599tm0y5AqhgFOkqpdmT99Yd9p023btqVdilQxDHCSVO3KeP+brUKk3LgHTpKq2fbt\nMGQIbNoEPXumXc0BDG+qdfnsgTPASVKlefhhmD07s7Hr1sH69TBrVnFrytKrr77K6aefzn333Wd4\nU80ywLUwwEmqCWPHwvve1/xkhUxMnAjvfndxa8rBq6++yjHHHJN2GVJqDHAtDHCSqt6GDXDKKbB5\nM3TunHY1kvJgGxFJqhWNjVBXZ3iTapwBTpIqSZm3BGnPzp070y5BqioGOEmqFDFWZIBLkoQzzjjD\nECcVUJe0C5AkZWjFCti9G8aMSbuSjLVuFdKjR4+0y5GqhjNwklQpGhubG/KGnPY8l5x93qTiMcBJ\nUqWYMaNilk8Nb1Jx2UZEkipBU1Nz37f582HYsLSr6dDPf/5zhgwZYniTDsM+cC0McJKq1v/+L1x3\nHSxblnYlkgrEPnCSVO327X+TJAxwklQZKmj/m6TicwlVksrdrl0wYACsXg39+qVdzSGSJKFr165M\nmDAh7VKkiuISqiRVs3nzYPTosg1v9fX17N69O+1SpJpigJOkclem+99sFSKlxwAnSeWuDPe/Gd6k\ndLkHTpLK2fbtMGQIbNoEPXumXQ0AW7duZezYsdx7772GNykP+eyB81moklTO5syB8ePLJrwB9OvX\nj0WLFtGvDPfkSbXCJVRJKmdluv/N8CalywAnSeWsDPe/SUqfAU6SytWmTfDSS3DOOamWsWPHjlTv\nL+lQBjhJKldPPAEXXABd0tuunCQJp512miFOKjMeYpCkcrViBZx8cmq3b90qpFevXqnVIelQzsBJ\nUrlauRJOOCGVW9vnTSpvBjhJKlcrV8KJJ5b8toY3qfyVRYALIVweQlgaQlgeQvhcO2PqQgj/G0JY\nFEKYVeoaJankVq1KZQZu+/bthjepzKX+JIYQQidgOXAJsA54Brguxri01Zi+wBPApTHGV0IIA2KM\nm9u4lk9ikFQdmpqam/du3VpWTXwlFU4+T2Iohxm4c4EXYoxrYoy7gfuAKw8a82HgVzHGVwDaCm+S\nVFXWr4ejjza8SWpTOQS4ocDaVt+/3PJaa6OBfiGEWSGEZ0IIHy1ZdZKUhhQPMEgqf5XSRqQL8C5g\nMtALeDKE8GSM8cV0y5KkIinRAYYkSdi1axeXXnpp0e8lqXDKIcC9Agxv9f2wltdaexnYHGN8G3g7\nhDAHOBM4JMDdeuut+39dV1fnJlxJlakEBxhanzaVVHxJkpAkSUGuVQ6HGDoDy2g+xLAeeBq4Psa4\npNWYMcB3gcuBI4CngGtjjIsPupaHGCRVhxtvhIsvhunTi3J5W4VI6avoQwwxxr3Ap4H/AZ4H7osx\nLgkh3BRC+GTLmKXAo8ACYB7w7weHN0mqKqtWFW0J1fAmVb7UZ+AKyRk4SVVj6FB48kkYPrzjsVnY\ntm0bY8eO5e677za8SSnLZwbOACdJ5WbnzuYWIjt2QOfOBb/8tm3bOPLIIwt+XUnZqeglVEnSQdas\naZ55K0J4AwxvUhUwwElSubEHnKQOGOAkqdwU8ADDG2+8UZDrSCovBjhJKjcFmoFLkoTTTjvNECdV\nIQOcJJWbAjyFYV+rkJ/+9Kf07du3QIVJKhcGOEkqN3kuodrnTap+BjhJKicx5rWEaniTaoMBTpLK\nydatze1Djj46p483NTUZ3qQaUA4Ps5ck7ZPn/rfJkycXsBhJ5coZOEkqJ/aAk5QBA5wklZMiPsRe\nUvUwwElSOcliBi5JEh588MEiFySpHBngJKmcZLgHbt9p06OOOqoERUkqNwY4SSonGSyh2ipEUogx\n5v7hELoBA4CdMcbXClZV7vXEfH4/kpSqPXugVy/Ytg2OOKLNIYY3qXqEEIgxhlw+m1UbkRBCH+A6\n4D3ARcAxrd7bAywAZgIPxBifyqUgSapZL78Mgwa1G9527NjBn//5nxveJGUW4EIIQ4EvAh8Gere8\n/DqwDNgK9AD6A+OAs4HPhhDmA3fEGH9e6KIlqSp1sP+tV69eLFy4kF69epWwKEnlqMMAF0L4R+Cv\ngSOA3wP3AY/HGFe0MbYXMB64DLgB+FkI4RbgkzHGBYUsXJKqTgYnUA1vkiCzQwyfBf4dGB5jfF+M\n8e62whtAjHFHjDGJMf49MAK4EugKfKhgFUtStbIHnKQMZbKEOirGuC7bC7ecJngYeDiEMDjryiSp\n1qxcCe9///5vt2zZQv/+/VMsSFK56nAGLpfw1sY1NuR7DUmqeq1m4JIkYezYsWzZsiXloiSVIx9m\nL0nlouUQQ+tWIc7ASWpLzo18Qwg3hBAaQwhbQwh7Wv45I4RwQyELlKSa8Oab8OabJEuW2OdNUoey\nbuQbQugK3A98AAjAXmAzzQ19OwMR+C1wTYxxd0Gr7bg2G/lKqkwLF5JccQX1O3YY3qQakU8j31xm\n4P4euAJ4CrgY6B5jHAJ0ByYDT9Mc7j6XS0GSVJNWrqT70KGGN0kZyWUG7kWgCRgbY9zVxvtHAIta\nrj2qIFVmXpszcJIq07/+a/MeuO98J+1KJJVIqWfghgG/aSu8AcQY3wF+AwzNpSBJqkkdPIVBklrL\nJcCto7k57+F0bRknScpEBk9hkKR9cmkjci/wsRDCl2KM2w5+M4RwFHANcGe+xUlSRdiwAV57LauP\nJE8/zSsbN3LDFVc0v7B8uTNwkjKWyx64bsAvgTHAPwJzgI3AIGASzQ+9XwJM9RSqpKoXIwwbBn36\nQMhsK0uyYwf169bRcOyx1O17tmmvXjBnDvTsWcRiJZWTfPbA5TIDt3PffYF72qoHOAl4Oxz4h1mM\nMdo4WFJ1ef556N4dli7NaPj+Jr0zZnjaVFLOcglUj9Hc602S1NgIl1yS0dDWT1gwvEnKR9YBLsZY\nV4Q6JKkyNTbCDR0/gObtt9/mpptuMrxJKois98CVM/fASSqpPXtgwAB44QU45pgOh7/99tt07969\nBIVJqgSl7gMnSQL4wx/g+OMzCm+A4U1SwXQY4EIInw0h5PynTgjhrBDCe3P9vCSVrRkzMt7/JkmF\nlMkM3FeBFSGEz4UQjs3koqHZZSGEXwN/AM7Mp0hJKkuHOcCwcePGEhcjqZZ0uAcuhDAa+DbwPmAv\n8AQwl+Zgth54jeYH2fenuTfcecAlwGBgC/Bl4Icxxqbi/BYOqNU9cJJK4623YNAgWL8eevc+4K0k\nSbj22mtZsGABgwYNSqlASeWuqH3gYozLgQ+EEC4AbgauBibSdiuRfUUsA24H/iPGuD2XwiSprD3+\nOIwb12Z429cqxPAmqVgybiMSY3wCeCKE8JfARcCFwHCaZ952ApuABUASY3y+CLVKUvloY/+bfd4k\nlYptRCQpF+ecA//yLzBxImB4k5S9fJZQDXCSlK2tW5vbh2zeDN26AbBgwQK2bt1qeJOUsaI/CzWE\ncCMwP8a4IJebSFJVSRK48ML94Q3gjDPOSK8eSTUn00a+PwE+1PqFEMK0EMLMglckSeXO/m+SUpbP\nkxiOByYVqA5JqhxZPMBekorBR2lJUjZefplkwwZ+9NRTaVciqYYZ4CQpC8n3vkf9229z0sknp12K\npBpmgJOkDCVJQv2//isNn/60p00lpSqbAGd/Dkk1a3+ft169qPvUp9IuR1KNy6gPXAihiewDXIwx\nZvykh0KwD5ykYti1axdnn3023/3sZ6n78pdh1SoIObVukqT9it7ItyXAZS3GWNIlWgOcpKz98z/D\nk092OGxXUxPdNmxofv7pnXeWoDBJ1c4nMbQwwEnKyp490L8//OhH0KNHZp8ZPx4GDy5uXZJqQtGf\nxCBJVemZZ5ofiTV1atqVSFJWPIUqqXa105D3lVdeSaEYScpc1gEuhHBmCOH2EMKsEMKiEMLCEMLM\nEMLXQginF6NISSqKNgJckiScddZZhjhJZS3jPXAhhM7Ad4FPAqHlq7XY8vV94JY0NqO5B05Sxt56\nCwYOhPXroU8foFWrkIYG+7xJKrpS7YG7A/hLYBfwSyABXqE5yB0LTAauAW4G3gb+NpeCJKkk5s5t\nPlFqeJNUgTJtIzIKWAqsBS6PMS5rZ9wY4L+BocDJMcaVBay1Q87AScrY5z4H3bvDV77C7Nmzueaa\nawxvkkoqnxm4TPfAfZTmmbaPtRfeAGKMS4FpQGfgI7kUJEkl0Wr/25AhQ7j//vsNb5IqRqYzcL8H\nhsQYx2Z00RAWAetijJfmWV9WnIGTlJGtW5vbh2zeDN26pV2NpBpVihm4McBTWVz3qZbPSFL5mTUL\nJkwwvEmqWJkGuKOATVlcdyNwdPblSFIJtNP/TZIqRaYBrhewM4vrvgP0zL4cSSquJEn4l/vvN8BJ\nqmg+iUFSzUiShPqrr+asd96BM89MuxxJylk2feA+FEI4PsOxZ2VfiiQVz/4+b5/4BHUrV0In//4q\nqXJlE+DGtXxlyuOgksrCAU1677zT5VNJFS/TNiLTcrl4jPE/c/lcrmwjIulge/bs4bzzzuOOO+6g\nbtIkGDqUDtAHAAAgAElEQVQU5syBUaPSLk1SjcunjUjGz0KtBAY4SW3Zs2cPXbp0gSVL4PLLYfVq\nCDn9mSlJBVOKPnCSVLG6dGnZLTJjBkyZYniTVPEyDnAhhE+FEP4+hND1MGO6tYz5vwpTniQVkP3f\nJFWJjAJcCOEC4LtAlxjj7vbGxRh3AUcA3wshvLswJUpS5latWkWbWyn27IHZs2Hy5NIXJUkFlukp\n1GnANuBbGYz9JvD/AH9Odo/fkqS87Dtt+vRPfsIJRx554JsrVjQfYBg8OJ3iJKmAMg1wE4GZMca3\nOhoYY9wRQmhs+YwklcT+ViF33cUJV18N48cfOuhTnyp9YZJUBJkGuOHAw1lc90XgsuzLkaTsHdDn\nrU8fGDMGHnss7bIkqWgyPcTQmewa88Ysri1JOXvsscf+FN7q6mDVKjjxxLTLkqSiynQG7lVgZBbX\nHQlszr4cScrOiBEjeOCBB5g4sWXXxsqVcMIJ6RYlSUWWaYB7BnhPCKFvjPGNww0MIfQF3gPMyLc4\nSerI8OHDGT58+J9eWLUKxo5NryBJKoFMlzl/DhwJ/L8ZjP0e0KflM5JUWitXuoQqqeplGuB+BTwB\nXB9CmB1CmBJC6LbvzZYGvlNCCAnwYeDxGOOvCl+uJHXAJVRJNSDjZ6GGEAYBjwJn0HxIYQ+wpeXt\n/jQvxwbgOeCyGOOmglfbcY0+C1WqYkmSMHfuXP7hH/6h7QF790KvXvD669C9e2mLk6QsleRZqDHG\njcD5wBeBtUBXYHDLV9eW1/4BuCCN8Capuu1rFXLhhRe2P+iVV6B/f8ObpKqX6SEGAGKMO4GvAl8N\nIQwDhrS8tT7G+HKhi5MkOKjPW11d+wNtISKpRmQV4FprCWyGNklFlXF4Aw8wSKoZmT7M/qIQwvCO\nR+4ff0YI4cbcy5IkaGpq4gtf+EJm4Q08wCCpZmS6B24W8LHWL4QQPhdC2NL2cK4C/iOPuiSJTp06\n8dhjj2UW3sAlVEk1I9MA19YJie7AUQWsRZIO0alTFk/lcwZOUo3weaWSqoczcJJqhAFOUtlYvnw5\nOfdyfOut5v5vQ4Z0PFaSKpwBTlJZSJKECRMm8MILL+R2gVWrYMQIyGbJVZIqlH/SSUpd61Yho0eP\nzu0iLp9KqiHZBDifUSWp4LLq83Y4HmCQVEOyaeR7awjh1oNfDCHsLVw5kmrJ448/XpjwBjbxlVRT\nspmBC1l+SdJhjRw5kgcffDD/8AYuoUqqKRnNwMUY3SsnqeAGDx7M4MGDC3Mxl1Al1ZCQ85H9MhRC\niNX0+5GUoRihTx945RXo2zftaiQpIyEEYow5rVo6syap8r36KnTvbniTVDMMcJJKIkkSPv/5zxfn\n4i6fSqoxBjhJRbevVcill15anBt4gEFSjTHASSqqgvV5Oxxn4CTVGAOcpKIpSXgDZ+Ak1RwDnKSi\niDFy2223FT+8gU18JdUc24hIKpoYIyGUoK/38cdDYyOMHFn8e0lSgVR8G5EQwuUhhKUhhOUhhM8d\nZtz4EMLuEMKflbI+SbkpSXjbvRvWr4fhw4t/L0kqE6kHuBBCJ+B7wGXAacD1IYQx7Yz7OvBoaSuU\nVNZeegmGDIGuXdOuRJJKJvUAB5wLvBBjXBNj3A3cB1zZxrj/G7gf2FTK4iRl5vnnn6epqan0N/YA\ng6QaVA4BbiiwttX3L7e8tl8I4VjgQzHG/w8owZqMpGwkSUJdXR1Lliwp/c09wCCpBpVDgMvEvwKt\n98YZ4qQy0bpVyGmnnVb6AuwBJ6kGdUm7AOAVoPXu42Etr7V2DnBfaN4RPQB4bwhhd4zxoYMvduut\nt+7/dV1dXfHbF0g1rGR93g5n1Sr40IfSubckZSFJEpIkKci1Um8jEkLoDCwDLgHWA08D18cY21yL\nCSH8B/BwjPGBNt6zjYhUIk8++SQf/OAH0w1vAOPHw3e/C+edl14NkpSDfNqIpD4DF2PcG0L4NPA/\nNC/p3hljXBJCuKn57fjvB3+k5EVKOsSYMWN46KGHOP/889MtxEMMkmpQ6jNwheQMnFRj3ngDhg6F\n7duhFD3nJKmAKr6RryTlZNWq5gMMhjdJNcYAJ6lyuXwqqUYZ4CR1aNasWdxyyy1pl3EoW4hIqlGp\nH2KQVN5mzZrF1KlTaWho+NOLGzbAk0+mV9Q+SQLveU/aVUhSyRngJLWrdXg7oFXIP/0TzJuX/gPk\nu3aFyZPTrUGSUuApVEltaje8AZx8MvziFzBuXCq1SVI18BSqpIKKMfIv//IvbYe3l1+GrVvhjDNS\nqU2S5BKqpDaEEPjNb35DaKs9R2MjXHwxdPLvf5KUFv8EltSmNsMbwIwZMGVKaYuRJB3APXCSMhdj\n85MPHnsMRo5MuxpJqmjugZOUlwULFrBnz56OBy5dCt262TxXklJmgJNqXJIkXHLJJSxevLjjwY2N\ncMklPrpKklJmgJNqWJIk1NfX09DQwBmZnCqdMaM5wEmSUuUeOKlGtQ5vh7QKacuePXDMMbBkCQwe\nXPT6JKnauQdOUlaefvrp7MIbwLPPwrBhhjdJKgP2gZNq0CmnnMJ//dd/ce6552b+oX373yRJqXMG\nTqpBffr0yS68gfvfJKmMuAdOUsd27mze/7ZuHRx5ZNrVSFJVcA+cpOJ64gk4/XTDmySVCQOcVOWS\nJOETn/hEfhdpbPTxWZJURgxwUhXb1yrkhhtuyO9C7n+TpLLiHjipSmXd5609r78Oxx0HmzfDEUcU\nrD5JqnXugZN0gIKFt+aLwfnnG94kqYwY4KQq9IMf/KAw4Q3s/yZJZcglVEmHd8op8NOfwtlnp12J\nJFWVfJZQDXBSJXniCfjYx6BU/57HCG+8ARs2QOfOpbmnJNUIA1wLA5yq3mc/2xyq/vIvS3fPvn1h\n4MDS3U+SakQ+Ac5noUqVpLERvvc9OOmk/S89++yznHrqqXTv3j3FwiRJpeQhBqlSbN4MK1dCq2eY\nJknCZZddxuLFi1MsTJJUagY4qVLMmgUTJ0LXrsCBrULe9a53pVycJKmUDHBSpWj1NISC9nmTJFUc\nDzFIlWLUKHjgAf64ezeXX3654U2SKpynUFsY4FS11qyB8eNhwwbe3rWL559/nrPtyyZJFc1HaUnV\nbt/TEDp1onv37oY3SapxBjipEvg4K0lSKy6hSuUuRhgyBJ58Ek44Ie1qJEkF4hKqVKWSJOHD738/\n9OxpeJMk7WeAk8rUvlYhnzzxRJgyJe1yJEllxAAnlaED+ry99JL73yRJB3APnFRmDghvF14IAwbA\nCy/AMcekXZokqYDcAydVkZ/+9Kd/atL7hz/AiBGGN0nSAZyBk8rZbbfBa6/Bt76VdiWSpAJzBk6q\nVvZ/kyS1wRk4qVy99RYMGgTr10Pv3mlXI0kqMGfgpAr19NNPs2PHjrbffPxxOPNMw5sk6RAGOCkl\nSZLw/ve/n+eff77tAY2N9n+TJLWpS9oFSBVt9WqYMSPrjyXLllH/gx/Q8Jd/ybkLFsCCBYcOevBB\n+NGP8q9RklR13AMn5eMTn4AXX4SRIzP+SLJ+PfUzZ9IweTJ1Q4a0P7BXL7jjDujatQCFSpLKTT57\n4AxwUj5OPBEefhhOOy2j4c899xxTpkz5U583SVLNMsC1MMCppFauhAkTYN06CJn997d7924WL17M\nmWeeWeTiJEnlzlOoUhoaG2Hy5IzDG0DXrl0Nb5KkvBngpFx5SlSSlBKXUKVcNDXB4MHNzyodPrzd\nYTFGQhYzdJKk2uESqlRqixZB376HDW9JknDVVVfhXyokSYVmHzgpFx08ozRJEurr62loaHAGTpJU\ncM7ASbmYMaPdANc6vNkqRJJUDO6Bk7K1ezcMGAArVjT/sxXDmyQpU+6Bk0rp6aebn7xwUHgD+NWv\nfmV4kyQVnTNwUrb+8R9h+3b45jfTrkSSVMGcgZNK6TD73yRJKgVn4KRs7NgBgwbBxo3ND5uXJClH\nzsBJpfLYY3D22dCrF08++SSvvfZa2hVJkmqQAU7KRkv/tyRJ+OAHP8jSpUvTrkiSVIMMcFI2Zswg\nOfro/a1Czj///LQrkiTVIPfASZnavJlkxAjqe/a0VYgkKW/ugZNK4Pm776Z+zx7DmyQpdc7ASRna\n+8lPsvTooznt9tvTLkWSVAXymYHzYfZShjrPnMlpDzyQdhmSJLmEKmVkzRrYtg3Gjk27EkmSDHBS\new5Yjm9shMmToZP/yUiS0uf/jaQ2JEnCZZdd9qcQ19L/TZKkcmCAkw6SJAn19fV8/vOfJ4QAMTYH\nuClT0i5NkiTAACcdYF94O6BVyPPPQ8+ecMIJqdYmSdI+BjipRZvhDVw+lSSVHQOc1OLRRx9tu0mv\nAU6SVGZs5Csdzp49MGAALF8OAwemXY0kqYr4KC2pWJ55BkaMMLxJksqKAU46HJdPJUllyACnmvT4\n44+zcePGjgca4CRJZcgAp5qTJAkf+tCHePHFFw8/8K23mpdQL7qoNIVJkpQhA5xqSutWIRMmTDj8\n4LlzYdw46NOnNMVJkpQhA5xqRrt93trj8qkkqUwZ4FQTli1bll14AwOcJKls2QdONSHGyPLlyzn5\n5JMz+8DWrXD88bB5M3TrVtTaJEm1yT5wUgdCCJmHN4BZs2DCBMObJKksGeCktrh8KkkqYwY4VaWm\npqb8LmCAkySVMQOcqk6SJEyaNCn3ELd2LWzZAmeeWdjCJEkqkC5pFyAVUutWIZ065fj3k8ZGmDwZ\ncv28JElFZoBT1Wi3z1uM8P3vw/btmV3ooYdg2rSi1ChJUiHYRkRV4bBNehctgilT4GMfy+xinTrB\nX/81DBhQ6DIlSdovnzYizsCpKjz22GPtN+ltbIQrroCvf73kdUmSVAzOwKn6XXEF3HADXHdd2pVI\nkrRfPjNwBjhVtz17oH9/eOEFGDgw7WokSdrPJzFI7XnmmeZHYhneJElVxACnivPYY4+xdu3azAbb\nkFeSVIUMcKooSZLwZ3/2Z6xevTqzDxjgJElVyD1wqhiHbRXSlrfeal46Xb8e+vQpen2SJGXDPXCq\nelmHN4C5c2HcOMObJKnqGOBU9lauXJl9eAOXTyVJVcslVJW9GCOrVq3ixBNPzO6D55wD3/42XHRR\ncQqTJCkP9oFrYYDTflu3NrcP2bwZunVLuxpJkg7hHjjpYLNmwYQJhjdJUlUqiwAXQrg8hLA0hLA8\nhPC5Nt7/cAjhuZavuSGE09OoU6Wxd+/e/C/i/jdJUhVLPcCFEDoB3wMuA04Drg8hjDlo2Ergohjj\nmcBtwI9KW6VKJUkSzjvvvPxDnAFOklTFuqRdAHAu8EKMcQ1ACOE+4Epg6b4BMcZ5rcbPA4aWtEKV\nROtWIZ07d879QmvXwpYtcOaZhStOkqQykvoMHM1hrPVzkV7m8AHtL4BHilqRSi6nPm/taWyEyZOh\nUzn86y1JUuGVwwxcxkIIFwPTgQvbG3Prrbfu/3VdXV3+YUBFV9DwBi6fSpLKUpIkJElSkGul3kYk\nhHAecGuM8fKW7/8OiDHG2w8adwbwK+DyGOOKdq5lG5EK9O1vf5t3vetdhQlvMcLQoTBnDowalf/1\nJEkqkoruAxdC6AwsAy4B1gNPA9fHGJe0GjMcaAQ+etB+uIOvZYCrdYsXw/veB6tWQcjpvwlJkkoi\nnwCX+hJqjHFvCOHTwP/QvCfvzhjjkhDCTc1vx38Hvgj0A74fQgjA7hjjuelVrZKKsfkrEzNmNC+f\nGt4kSVUs9Rm4QnIGrgrt2QPDhsGmTZmNDwHuvx+uuqq4dUmSlKeKnoFTbZk9ezbHHnssJ510UmYf\neOYZGDQINmwobmGSJFUQ+yyoZJIk4ZprrmHdunWZf8gTpZIkHcIAp5Jo3Spk0qRJmX+wsRGmTCle\nYZIkVSD3wKnocu7z9tZbMHBg8/Jp795Fq0+SpDTkswfOGTgV1dq1a5k6dWpuTXrnzoVx4wxvkiQd\nxEMMKqrjjjuOP/7xjxx33HHZf9j9b5IktckZOBVdTuEN3P8mSVI73AOn8rR1Kxx/PGzeDN26pV2N\nJEkF5x44lY3du3cX5kKzZsGECYY3SZLaYIBTwSRJwtlnn12YEOf+N0mS2uUhBhVE61YhXbt2zf+C\njY3wyU/mfx1JkqqQM3DKW8593trz8svNe+DOOCP/a0mSVIUMcMpLwcMbNM++XXwxdPJfT0mS2uL/\nIZWXpUuXFja8AcyY4f43SZIOwzYiKi8xwtCh8NhjMHJk2tVIklQ0thFR9Vi6tLl1yIknpl2JJEll\nywCn8rKvfUjI6S8kkiTVBAOcMpYkCQsXLizuTez/JklShwxwysi+06Zbtmwp3k327oUkMcBJktQB\nA5w6VJRWIW354x9h2DAYNKh495AkqQoY4HRYJQtv4PKpJEkZMsCpXRs2bODaa68tTXgDA5wkSRmy\nD5wOa8OGDQwePLj4N9q5EwYOhFdegSOPLP79JElKmX3gVDQlCW8ATzwBp59ueJMkKQMGOJUHl08l\nScqYAU77vfPOO+nd3AAnSVLGDHACmk+bnnnmmemEuNdfh8WL4fzzS39vSZIqUJe0C1D6WrcKOeKI\nI9IoAC64ANK4tyRJFcgZuBpX0j5v7XH5VJKkrBjgalhZhDcwwEmSlCUDXA1bu3Zt+uFt3TrYuBHG\njUuvBkmSKoyNfJWue+6B3/wG7r8/7UokSSopG/mqcrl8KklS1gxwSk+MBjhJknJggKsRSZLw1FNP\npV3GgV54ofmfJ52Ubh2SJFUYA1wN2HfadOfOnWmXcqAZM2DKFAg5Lf9LklSzDHBVrmxahbTF5VNJ\nknLiKdQqVtbhbe9eGDgQFi6EY49NuxpJkkrOU6g6xJYtW7j++uvLM7wBzJ8PgwYZ3iRJyoHPQq1S\n/fv3Z9GiRfTv3z/tUtq2b/+bJEnKmjNwVaxswxu4/02SpDy4B66SbN8O5XaSNBe7d8OYMbB2LRx1\nVNrVSJKUinz2wLmEWil27YLhw6Fr1zbffitGelZSO47Jkw1vkiTlyABXKebNg1Gj4JlnDnkrSRI+\n/vGPs3DhQnr27JlCcZIkqZQMcJWinT1jrVuFGN4kSaoNHmKoFDNmHBLgyrrPmyRJKhoPMVSC7dth\nyBDYtAlaZtkMb5IkVTYb+Va7OXNg/Pj94Q2aG/Ua3iRJqk3ugasEjY2HNL29+uqrUypGkiSlzRm4\nSmDTW0mS1Ip74Mrdpk0wejRs3gxdnDCVJKlauAeums2cSXLqqSRz56ZdiSRJKhMGuDKX3HMP9QsW\npF2GJEkqI67JlbEkSaj/7/+m4cc/9rSpJEnazz1wZSpJEuqvvpqGGKnbsgUq6TmnkiSpQz7Mvsq8\n/vrrfOQjH6Fh+nTq1q83vEmSpAO4B64MHXXUUSxatIi6l18+pP+bJEmSS6jlqqkJBg+GP/wBhg9P\nuxpJklRgthGpRosWQd++hjdJknQIA1wZ2L59+6Ev+vQFSZLUDgNcypIk4bTTTjs0xM2Y4f43SZLU\nJgNcipIkob6+nrvvvps+ffr86Y3du2HuXLj44vSKkyRJZcsAl5J94a2hoeHQJr1PPw0jR0L//qnU\nJkmSypunUEvl61+HpUsBSDZsoD5JaKiro27w4EPHLl0KEyfCN79Z4iIlSVKp5HMK1QBXCtu3w7HH\nwne+AyHw3wsX0r1rV+rGjGn/M5ddBkOGlK5GSZJUUj6JodzNmQPjx8P06QBcnnI5kiSpsrkHrhRs\nCSJJkgrIAFcKBjhJklRABrgiS379a3774otwzjlplyJJkqqEAa6IkiShfto0eo8bB13cbihJkgrD\nVFEk+/u8TZxI3WWXpV2OJEmqIs7AFcEBTXqXLHH/myRJKij7wBXYm2++ydixY/nJT35C3YgRcMEF\nsG4dhJzavEiSpCplH7gy0rt3bxYtWkTv3r3hxz9unn0zvEmSpAJyCbUIevfu3fyLGTNcPpUkSQXn\nEmqxNDXB4MHwhz/A8OFpVyNJkspMPkuozsDl6bXXXmv7jUWLoG9fw5skSSo4A1wekiRh7NixbYe4\nxkaYMqX0RUmSpKrnIYYctW4VcvTRRx86YMaM/Q+vlyRVtuOPP541a9akXYYq1IgRI1i9enVBr+ke\nuBwc0Oetru7QAbt3w4ABsHIl9O9f9HokScXVslcp7TJUodr798c9cCXUYXgDePppGDnS8CZJkorC\nAJelTp06HT68gfvfJElSUbmEWgwXXQRf+AL4DFRJqgouoSofLqFWgh074Nln4cIL065EkqSy06lT\nJ1auXAnA9OnT+dKXvlSU+4wdO5Y5c+YU5drlwACXjxUrYMmSA79++Us4+2zo1Svt6iRJNeCEE05g\n5syZaZeRsVCix0suWrSIiy66KOfPf+c73+H000+nd+/eDB8+nGuvvZZFixZx++23M2nSpEPGb9my\nhSOOOILFixfnU3bGbCNyGEmSsGnTJqZOnXrom4sWwbnnwogRh773f/5P8YuTJCkDe/fupXPnzmmX\nsV8lLEX/1V/9FY888gg//vGPueCCC9i7dy+//vWv+d3vfsdHPvIRvvjFL7JmzRpGtMoAP//5zznj\njDM49dRTS1KjM3Dt2HfadODAgW0P+J//gRtvPHQGbskS+Iu/KG2xkqSadOONN/LSSy9xxRVXcOSR\nR3LHHXewZs0aOnXqxF133cWIESO45JJLmD17Nscdd9wBn209cxdj5Otf/zqjRo3imGOO4brrruP1\n119v8577rvXtb3+bQYMGMXToUH7yk5/sf//iiy/mrrvu2v/9f/7nfzJx4sS8f6/Tp0/n5ptv5n3v\nex99+vRh4sSJbNy4kc985jP069ePU089leeee67N399XvvIVrr32WqZNm8aRRx7J6aefzrPPPtvm\nfV588UW+//3vc9999zFp0iS6du1K9+7duf766/nbv/1bhg4dysUXX8w999xzwOfuuecepk2blvfv\nM1MGuDZk1CrEk6aSpJTdfffdDB8+nN/+9rds27aNz372s/vfmzNnDkuXLuXRRx8FDr90+Z3vfIeH\nHnqIxx57jHXr1nH00UfzqU99qt3xGzZsYPv27axbt44f//jH3Hzzzbzxxhvtji/UsmlDQwP//M//\nzJYtW+jWrRvnn38+55xzDlu2bOHqq6/mM5/5TLufffjhh/nwhz/MG2+8wRVXXMHNN9/c5rjGxkaO\nO+44zj777HavNW3atAMC3LJly3juuee4/vrrc//NZckAd5CMwtvu3TB3Llx8cUlrkySpLQcvS4YQ\n+MpXvkKPHj044ogjOvz8D3/4Q7761a8yZMgQunbtype+9CXuv/9+mpqa2hzfrVs3vvjFL9K5c2fe\n+9730rt3b5YtW1aQ38vhXHXVVYwbN45u3bpx1VVX0aNHD2644QZCCFx77bXMnz+/3c9eeOGFXHbZ\nZYQQ+OhHP8qCBQvaHLdlyxaGDBnSYR0bN25k3rx5QPPs23vf+176l7D/qwGulZ07d/IXf/EXHfd5\ns1GvJKm1EArzVUDDhg3LeOyaNWu46qqr6Nev3/7lyK5du7Jx48Y2x/fv359Onf4UIXr27Mmbb76Z\nd80dGTRo0P5f9+jR45DvD1fD4MGD9/+6Z8+evP32220G1P79+7N+/frD1tGjRw+uueYa7r77bgB+\n9rOflXT5FAxwB+jRowcLFy48fHiD5uXTSy4pSU2SpAoQY2G+ctDe8mTr13v16sVbb721//u9e/fy\n6quv7v9++PDhPPLII2zdupWtW7fy2muvsWPHjg5notpy8L02bNiQ9TXSdMkll/Dyyy+3u0dun2nT\npvHLX/6S3//+97z55pt84AMfKFGFzQxwB+nRo0fHg2bMMMBJksrC4MGD9/dV2+fgJdXRo0fz9ttv\n88gjj7Bnzx5uu+02du3atf/9m266ic9//vO89NJLALz66qs89NBDOdUzbtw4HnjgAXbu3MmLL77I\nnXfemfFnO3XqlFfvtmxOuLY3dtSoUXzqU5/i+uuvZ/bs2ezevZt33nmHX/ziF3zjG9/YP27ixIn0\n7duXT37yk1x33XV06VLaxh4GuGzta9RbgBM1kiTl6+/+7u/4p3/6J/r168e3v/1t4NBZuSOPPJLv\nf//7fPzjH2fYsGH06dPngCXWW265hSuvvJJLL72Uvn37csEFF/D0009nXEPr+33mM5+ha9euDB48\nmOnTp/ORj3yk3bGtrV27dv8J0Y7ukUkdHY0/3Pv/9m//xqc//Wluvvlmjj76aEaNGsWDDz7IFVdc\nccC4faeAb7zxxg5rK7SafpTWq6++yjHHHJPdTf77v+FrX4PZs7OsTpJUqXyUVvH97Gc/Y/HixXz1\nq19Nu5SCK8ajtGo2wCVJwnXXXcfChQuzC3F/8zfQpw8U6dEfkqTyY4BTPnwWaoHsaxVy3333ZT8D\n5/43SZKUspqbgcuoz1t7Nm9ubh+yeTN07Zp7oZKkiuIMnPLhDFye8gpvALNmNR9eMLxJkqQU1VSA\n69OnT+7hDez/JkmSykLNLaHmZdQoeOABOOOM4t1DklR2XEJVPlxCTdOaNbBtG4wdm3YlkiSpxpW2\nbXAla2yEyZOhk5lXkmrNiBEjMmokK7VlxIgRBb9mWQS4EMLlwL/SPCN4Z4zx9jbGfAd4L7AD+FiM\ncf7hrpkkCatWrWL69OmFKbKxEaZMKcy1JEkVZfXq1WmXIB0g9emkEEIn4HvAZcBpwPUhhDEHjXkv\nMDLGeBJwE/CDw11z32nTE044oTBFxugBhhJIkiTtEpQjf3aVzZ9fZfPnV5tSD3DAucALMcY1Mcbd\nwH3AlQeNuRK4GyDG+BTQN4QwqK2L5d0qpC3PPw89e0KhAqHa5B9ClcufXWXz51fZ/PnVpnIIcEOB\nta2+f7nltcONeaWNMQCFD2/g7JskSSorZbEHrpAaRo2i7lvfgm99q3AXfe45+MY3Cnc9SZKkPKTe\nBy6EcB5wa4zx8pbv/w6IrQ8yhBB+AMyKMf6i5fulwKQY48aDrmWTHkmSVDFy7QNXDjNwzwCj/v/2\n7vBnIJMAAAUiSURBVD/U7rqO4/jzReIfik00KKQ2Upgi5Y9hDmKEk8It+sPESQ2UwmCQRv9Z/4h/\nNEgQIcRfLPrxj1KhpPsjmf4xmZMNGs4fpVbqKLQ01BwtGKz19o9z1NPd+d77vXfb93u+3OcDDnf3\ne78XXuzN99wXn3O+n5NkFfAP4BvAN+ecsx24Cfj1uPC9N7e8wdL/EyRJkoak9wJXVUeT3Aw8zkfb\niLyUZMvox7Wtqn6X5KtJXmG0jcgJ2htEkiRpeHp/CVWSJEmLMwt3oS5akg1JXk7y5yQ/aDjnriR/\nSfJskku6zqjpFppdks1Jnhs/dif5fB85NV2ba2983heSHElyTZf5NL+Wz51XJNmf5A9JdnadUdO1\neO78eJLt4795LyT5Vg8xNUWSnyV5K8nz85yz6M4yuAJ3Mjb+VTfazA54DfhSVV0MbAV+2m1KNWk5\nvw/Oux3Y0W1Czaflc+cK4B7ga1X1OWBT50F1jJbX3k3AH6vqEmA9cGeS3t8mJQB+wWh2Uy21swyu\nwHGCN/5VpxacXVXtraqD42/30rDfn3rR5toD+B7wEPDPLsNpQW3mtxl4uKreAKiqtzvOqOnazK6A\nM8b/PgN4p6r+22FGNaiq3cC/5jllSZ1liAXuhG78q061md2k7wCPndREWowF55fkHODqqroP8K7w\n2dLm+lsNnJVkZ5LfJ7m+s3SaT5vZ3Q1cmOTvwHPA9zvKpuO3pM7i8qpmUpL1jO42Xtd3Fi3KT4DJ\n9+dY4oblFGANcCVwOrAnyZ6qeqXfWGrhKmB/VV2Z5DzgiSQXVdWhvoPp5BhigXsDWDnx/afHx+ae\n85kFzlH32syOJBcB24ANVTXfsrO61WZ+lwG/ShLgE8DGJEeqantHGdWszfxeB96uqsPA4SS7gIsB\nC1y/2szu28CPAarq1SQHgAuAfZ0k1PFYUmcZ4kuoH278m+RURhv/zv3jsB24AT78pIepG/+qcwvO\nLslK4GHg+qp6tYeMarbg/Krq3PHjs4zeB/ddy9vMaPPc+SiwLsnHkpwGrAVe6jinjtVmdn8Fvgww\nfv/UakY3hWk2hOZXJJbUWQa3AufGv8PVZnbArcBZwL3jVZwjVXV5f6n1gZbz+79f6TykGrV87nw5\nyQ7geeAosK2qXuwxtmh97W0FfjmxVcUtVfVuT5E1IcmDwBXA2Un+BtwGnMpxdhY38pUkSRqYIb6E\nKkmStKxZ4CRJkgbGAidJkjQwFjhJkqSBscBJkiQNjAVOkiRpYCxwkiRJA2OBkyRJGhgLnCRJ0sBY\n4CRJkgbGAidJc4w/NPx/SX6e5PwkjyR5J8mhJE8l+UrfGSUtbxY4SWp2LrAHOBO4H/gNsAZ4LMmm\nPoNJWt78MHtJmiPJKuAAUMAdVfXDiZ+tAfYC/wZWVdWhflJKWs5cgZOkZgeBH00eqKpngAcYrcp9\nvY9QkmSBk6Rmz1TVf6YcfxIIcGm3cSRpxAInSc3eajj+5vjriq6CSNIkC5wkNftkw/FPjb8e7CqI\nJE2ywElSszVJTp9yfD2jGxz2d5xHkgALnCTNZwVw2+SBJJcBm4H3gN/2EUqSTuk7gCTNsF3AjUnW\nAk8D5wDXMbqBYYtbiEjqiytwktTsAPBF4F1gC3AtsA/YWFUP9RlM0vLmCpwkzaOq/oT7vUmaMa7A\nSZIkDYwFTpIkaWAscJI0XY0fkjRz/DB7SZKkgXEFTpIkaWAscJIkSQNjgZMkSRoYC5wkSdLAWOAk\nSZIGxgInSZI0MO8Dh/G2UkcjAsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115e6e990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "U = np.linspace(0, 1, 100)\n",
    "plt.figure(figsize=(10,10))\n",
    "if(len(P0['1SE'])):\n",
    "    plt.plot(U, sm.distributions.ECDF(P0['1SE'])(U), 'k', label='true null, 1 SE')\n",
    "if(len(P0['minCV'])):\n",
    "    plt.plot(U, sm.distributions.ECDF(P0['minCV'])(U), 'r', label='true null, min CV')\n",
    "plt.plot([0,1], [0,1], 'k--')\n",
    "plt.gca().set_xlabel('p', fontsize=20)\n",
    "plt.gca().set_ylabel('ECDF(p)', fontsize=20)\n",
    "plt.legend(loc='lower right')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

other_multiples.ipynb