karpEtal15.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Re-analyze Karp et al. (2105) data\n",
    "### P.B. Stark\n",
    "\n",
    "Python 3 notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.stats import binom, chi2, hypergeom\n",
    "import numbers\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "import re\n",
    "import pandas as pd\n",
    "import itertools\n",
    "import permute as pt\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "from permute.utils import binom_conf_interval\n",
    "from cryptorandom import cryptorandom as cr\n",
    "from cryptorandom import sample\n",
    "\n",
    "seed = \"statistics 159/259 fall 2018\"\n",
    "prng = cr.SHA256(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = '../Data/karp_farms.csv'\n",
    "with open(fn, newline='') as csvfile:\n",
    "    reader = csv.DictReader(csvfile)\n",
    "    farms = [r for r in reader]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "odict_keys(['Farm ID', 'county', 'town', 'E_Coli_Change', 'Salm_Change', 'Change in Nature', 'Change in Riparian', 'Change in Crop'])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# what are the column names?\n",
    "farms[0].keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecoli = np.array([float(f['E_Coli_Change']) for f in farms])     # the change in E coli \n",
    "nature = np.array([float(f['Change in Nature']) for f in farms]) # the change in natural vegetation\n",
    "\n",
    "# values for regression line plot\n",
    "nature_line = np.linspace(np.amin(nature), np.amax(nature), num=100)\n",
    "\n",
    "# augment data with a vector of 1s for regression with an intercept\n",
    "nature1 = sm.add_constant(nature)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.003921569"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(ecoli)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_regress(y, x):\n",
    "    regress = sm.OLS(y, x).fit()  # fit the linear regression model using OLS\n",
    "    print(regress.summary())\n",
    "    coeffs = regress.params                 # coeffs[0] is the intercept; coeffs[1] is the slope\n",
    "    plt.scatter(x[:,1], y,  color='blue')                   # data\n",
    "    plt.plot(nature_line, coeffs[0] + coeffs[1] * nature_line)  # regression line\n",
    "    plt.xlabel('Relative change in natural vegetation')\n",
    "    plt.ylabel('Relative change in EHEC')\n",
    "    plt.show()\n",
    "    return(coeffs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.299\n",
      "Model:                            OLS   Adj. R-squared:                  0.272\n",
      "Method:                 Least Squares   F-statistic:                     11.08\n",
      "Date:                Sun, 16 Sep 2018   Prob (F-statistic):            0.00262\n",
      "Time:                        20:37:19   Log-Likelihood:                 57.716\n",
      "No. Observations:                  28   AIC:                            -111.4\n",
      "Df Residuals:                      26   BIC:                            -108.8\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.0033      0.007      0.492      0.627      -0.011       0.017\n",
      "x1            -0.3036      0.091     -3.328      0.003      -0.491      -0.116\n",
      "==============================================================================\n",
      "Omnibus:                       27.451   Durbin-Watson:                   2.072\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               78.562\n",
      "Skew:                           1.760   Prob(JB):                     8.72e-18\n",
      "Kurtosis:                      10.413   Cond. No.                         15.1\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl81Xed7/HXh4QkhC2BHCgFkkAh\nVFortOneQh1rtztWZ+x0Eb11dORq3e6MS9tBr7Zal3rvjM7oqOhotcWldRvsqFhrgVZKW6B0L0lY\nm0JLwhJ2yPK5f/y+ORwOyTkncE5ykryfj8d55Lee3+f7O8n55Pf7fn/fr7k7IiIiqQzp6wBERCT/\nKVmIiEhaShYiIpKWkoWIiKSlZCEiImkpWYiISFpKFiIikpaShYiIpKVkISIiaRX2dQDZUlFR4dXV\n1X0dhohIv7J69epmd4+l227AJIvq6mpWrVrV12GIiPQrZrY5k+10G0pERNJSshARkbSULEREJC0l\nCxERSUvJQkRE0lKyEBGRtHKaLMzsKjNbZ2YNZnZbF+v/ycxeNLNnzexhM6tKWHezmdWH1825jFNE\nRFLLWbIwswLgW8DVwEzgJjObmbTZ00Ctu58F/AK4O+w7BvgccD5wHvA5MyvPVawiIpJaLq8szgMa\n3H2Dux8Bfga8PXEDd3/E3Q+E2ZXApDB9JfCQu+90913AQ8BVOYxVRERSyGWymAi8kjDfGJZ15/3A\n709wXxERyaFcdvdhXSzzLjc0ezdQC8ztyb5mNh+YD1BZWXliUYqISFq5vLJoBCYnzE8CtiZvZGaX\nAwuAa939cE/2dfeF7l7r7rWxWNp+sERE5ATlMlk8BUw3sylmVgTcCCxO3MDMZgPfJUoU2xNWLQGu\nMLPyULF9RVgmIiJ9IGe3ody9zcw+QvQlXwD8wN1fMLM7gVXuvhj4GjACeMDMALa4+7XuvtPMvkCU\ncADudPeduYpVRERSM/cuqxH6ndraWlcX5SIiPWNmq929Nt12eoJbRETSUrIQEZG0lCxERCQtJYsB\nbtEiqK6GIUOin4sW9XVEItIfDZgxuOV4ixbB/PlwIHSosnlzNA8wb17fxSUi/Y+uLAawBQuOJopO\nBw5Ey0VEekLJYgDbsqVny0VEuqNkMYB1112WutESkZ5SshjA7roLSkuPXVZaGi0XEekJJYsBbN48\nWLgQqqrALPq5cKEqt0Wk59QaaoCbN0/JQUROnq4sREQkLSULERFJS8lCRETSUrIQEZG0lCxERCQt\nJQsREUlLyUJERNJSshARkbSULEREJC0lCxERSUvJQkRE0lKyEBGRtJQsREQkLSULERFJS8lCRETS\nUrIQEZG0lCxERCQtJQsREUlLyUJERNJSshARkbSULEREJK1uk4WZXWlm13WxfJ6ZvTW3YYmISD5J\ndWVxB7Csi+UPA3fmJhwREclHqZJFqbs3JS9099eA4bkLSURE8k2qZFFiZoXJC81sKDAsdyGJiEi+\nSZUsfgV8z8ziVxFh+jthXVpmdpWZrTOzBjO7rYv1c8xsjZm1JdePmFm7ma0Nr8WZFUdERHIhVbL4\nDPA6sNnMVpvZamAT0BTWpWRmBcC3gKuBmcBNZjYzabMtwHuBn3TxFgfdfVZ4XZvueCIikjvH3Wbq\n5O5twG1mdgcwLSxucPeDGb73eWH7DQBm9jPg7cCLCcfYFNZ19Dx0ERHpLamazn4aICSH0939uc5E\nYWZfyuC9JwKvJMw3hmWZKjGzVWa20sze0YP9REQky1LdhroxYfr2pHVXZfDe1sUyz2C/TpXuXgu8\nC/i6mZ123AHM5oeEsqqp6biGWyIikiWpkoV1M93VfFcagckJ85OArRnGhbtvDT83AEuB2V1ss9Dd\na929NhaLZfrWIiLSQ6mShXcz3dV8V54CppvZFDMrIrpSyahVk5mVm1lxmK4ALiahrkNERHpXtxXc\nwJvMbA/RVcSwME2YL0n3xu7eZmYfAZYABcAP3P0FM7sTWOXui83sXODXQDnwNjO7w93PAN4AfDdU\nfA8BvuLuShYiIn3E3HtSjZC/amtrfdWqVX0dhohIv2Jmq0P9cErdXlmY2ZhUO7r7zhMJTERE+p9U\nt6FWE9VNGDCBqHK6s2Lbgam5DU1ERPJFqofypnROm9nT7n5cayQRERkcMh38aGBUbIiIyAnRSHki\nIpJWqgruf0qYHZc0j7v/S86iEhGRvJKqgntkwvT3kuZFRGQQSVXBfUdvBiIiIvkrVa+z9ydMfzVp\n3R9zGZSIiOSXVBXc0xOm35q0Tr32iYgMIpl2JNiTdSIiMsCkquAuNbPZRAllWJi28BrWG8GJiEh+\nSJUstgGdzWNfS5junBcRkUEiVWuoN/dmICIikr/0BLeIiKSlZCEiImkpWYiISFqpKrjjzGwiUJW4\nvbsvz1VQIiKSX9Imi/D09g3Ai0B7WOyAkoWIyCCRyZXFO4AZ7n4418GIiEh+yqTOYgMwNNeBiIhI\n/srkyuIAsNbMHgbiVxfu/rGcRSUiInklk2SxOLxERGSQSpss3P1HvRGIiIjkr1TDqt7v7teb2XN0\n0cusu5+V08hERCRvpLqy+Hj4+de9EYiIiOSvVB0Jbgs/N/deOCIiko/U3YeIiKSlZCEiImlllCzM\nbJiZzch1MCIikp/SJgszexuwFvhDmJ9lZnruQkRkEMnkyuLzwHnAbgB3XwtU5y4kERHJN5kkizZ3\nb8l5JCIikrcy6e7jeTN7F1BgZtOBjwErchuWiIjkk0yuLD4KnEHUieBPgT3A/85lUCIikl8y6Rvq\nALAgvEREZBDKpDXUb81scdLrXjP7uJmVpNn3KjNbZ2YNZnZbF+vnmNkaM2szs+uS1t1sZvXhdXPP\niyYiItmS6eBH+4Dvhdce4HWgJsx3ycwKgG8BVwMzgZvMbGbSZluA9wI/Sdp3DPA54HyillifM7Py\nDGIVEZEcyKSCe7a7z0mY/62ZLXf3OWb2Qor9zgMa3H0DgJn9DHg70VjeALj7prCuI2nfK4GH3H1n\nWP8QcBVRnYmIiPSyTK4sYmZW2TkTpivC7JEU+00EXkmYbwzLMnEy+4qISJZlcmXxCeAxM1sPGDAF\nuMXMhgOpBkayLpYdNy7GyexrZvOB+QCVlZXH7SAiItmRSWuo34XnK04n+hJ/2d0PhdVfT7FrIzA5\nYX4SsDXDuBqBy5L2XdpFbAuBhQC1tbWZJiIREemhTHudPYfoWYuzgOvN7H9msM9TwHQzm2JmRcCN\nZD6W9xLgCjMrDxXbV4RlIiLSB9JeWZjZvcBpRJ0JtofFDvw41X7u3mZmHyH6ki8AfuDuL5jZncAq\nd19sZucCvwbKgbeZ2R3ufoa77zSzLxAlHIA7Oyu7RUSk95l76rs3ZvYSMNPTbdjHamtrfdWqVX0d\nhohIv2Jmq929Nt12mdyGeh445eRDEhGR/iqT1lAVwItm9iRR/1AAuPu1OYtKRETySibJ4vO5DkJE\nRPJbJk1nl/VGICIikr8y6UjwAjN7ysz2mdkRM2s3sz29EZyIiOSHTCq4vwncBNQDw4B/CMtERGSQ\nyKTOAndvMLMCd28HfmhmGilPRGQQySRZHAhPYK81s7uBbcDw3IYlIiL5JJPbUO8J230E2E/U39M7\ncxmUiIjkl0xaQ20Ok4eAO3IbjoiI5KNM+oa6mOhZi6rE7d19au7CEhGRfJJJncV/Av8IrOZoR4KD\nyrONu/nS715ibs045tbEeMOEkZh1NeSGiMjAlEmyaHH33+c8kjy252AbLQfb+OofXuarf3iZcSOL\nuXR6jDk1FVw6PcaY4UV9HaKISE512+usmZ0dJq8n6mL8VxzbN9SanEfXA73R6+z2PYdYXt/M0nXb\neayhmd0HWjGDsyaOZm5NjDk1MWZNLqOwINNhQkRE+lamvc6mShaPpNjP3f2vTjS4XOjtLsrbO5xn\nG3ezrK6J5XVNrH1lNx0Oo0oKuWR6RTx5TBg9rNdiEhHpqZNOFv1NX49n0XKglccamlle18SyuiZe\n2xONPDt93Ih44jhvyhhKhhb0WYwiiRYtggULYMsWqKyEu+6CefP6OirpbVlLFmb2JeBud98d5suB\nT7j7Z7ISaZb0dbJI5O7Uvb4vnjie3LiTI+0dlAwdwvlTxjK3JsbcGTGmVgxXRbn0iUWLYP58OHDg\n6LLSUli4UAljsMlmsnja3WcnLVvj7md3t09fyKdkkezAkTae2LAzfstqQ/N+ACaWDWPujBhza2Jc\ndNpYRpYM7eNIZbCorobNm49fXlUFmzb1djTSlzJNFpm0hiows2J3PxzeeBhQfLIBDialRYW8+fRx\nvPn0cQC8svMAS0Pi+K+nX+UnT2yhcIhxdlV5dNVRE2PmhFEMGaKrDsmNLVt6tlwkkyuLTwPXAj8E\nHHgfsNjd7859eJnL5yuLVI60dbBmy674LasXtka9v1eMKDqmeW7FCOVnyR5dWUinrFZwm9lVwOWA\nAX909yUnH2J29ddkkaxp72EerY8Sx6P1zezcfwSAMyeOClcd45hdWcZQNc+Vk6A6C+mk1lADQEeH\n8/zWFpati5LH06/spr3DGVlcyEXTxjK3ZhxzaiqYVF7a16FKP6TWUAJKFgNSy8FWVjQ0s7y+iWXr\nmtjaEjXPPS02PJ44Lpg6Vs1zRSRjShYDnLuzvmkfy+qaWVbXxBMbdnC4rYPiwiGcN2VMvKJ82rgR\nap4rIt3Kdp3FMKDS3ddlI7hcGGzJItmh1nae2Lgz3LLazvqmqHnuqaNLmBMSx0XTKhg9TM1zReSo\nbD5n8Tbg/wJF7j7FzGYBd7r7tdkJNTsGe7JI1rjrAMvrmllWt50VDTvYe7iNgiHG7Mll8SfK3zhx\ntJrnigxy2UwWq4G/ApZ2PpxnZs+6+1lZiTRLlCy619rewdpXdrN03XaW1zXz3KstAIwZXsQl06J+\nrC6tqWDcyJI+jlREels2H8prc/cW3ffuv4YWDOHc6jGcWz2GT10JO/Yd5tH6qB+r5fVNLH5mKwAz\nJ4xi7owYc6bHOKeqnKJCNc8VkUgmyeJ5M3sX0ZPc04GPAStyG5bk0tgRxbxj9kTeMXsiHR3Oi9v2\nsLy+iaXrmvje8g18e+l6hhcVcNG0CubUxLisJsbkMWqeKzKYZXIbqhRYAFwRFi0Bvujuh3IcW4/o\nNlR27D3Uyor1O+L9WDXuOgjA1IrhzKmJxZvnlhZl8n+GiOS7bNZZzHb3p7MWWY4oWWSfu7OxeX+8\nK5LHN+zgUGsHRQVDOHdKebyifMZ4DTMr0l9lM1k8AkwAHgB+5u4vZCfE7FKyyL1Dre08tWlnPHnU\nvb4PgPGjiuOJ45JpFZSVaphZkf4i289ZnEI0vOoNwCjg5+7+xZOOMouULHrftpaD8cTxWH0zew61\nMcTgTaF57tyaGGdNKqNAzXNF8lZOnuA2szcCnwZucPe8+vdRyaJvtbV38Ezj7uihwPpmnm3cjTuU\nlQ6NN8+dUxNj/Cg1zxXJJ9m8DfUGoiuK64AdwM+AX7r79mwEmi1KFvll5/4jxwwz27T3MACnnzIy\nftVxTnU5xYXqx0qkL2UzWawEfgo84O5bsxRf1ilZ5C9356Vte+MtrFZt3klru1NaVMCFU8fGuyOp\nrhje16GKDDp50ZFgGAfjG0AB8H13/0rS+mLgx8A5RFctN7j7JjOrBl4COvuiWunuH0x1LCWL/mPf\n4TZWhua5y+qa2LIzGlShamwpc6ZHiePC08YyvFjNc0Vy7aSf4Daz+939ejN7jmiEvPgqwNN192Fm\nBcC3gLcCjcBTZrbY3V9M2Oz9wC53n2ZmNwJfJbrlBbDe3WelK4D0PyOKC7l85ngunzkegE3N++Pd\nrv9yTSP3rtzM0AKjtmpM/InyN0xQ81yRvtTtlYWZTXD3bWZW1dV6d+9iUMZj9r8Q+Ly7Xxnmbw/7\nfTlhmyVhm8fNrBB4DYgBVcCD7n5mpgXRlcXAcLitndWbdsWvOl5+bS8A40YWHzPM7JjhedW+QqTf\nOukrC3ffFiZvcfdbk978q8Ctx+91jInAKwnzjcD53W3j7m1m1gKMDeummNnTwB7gM+7+aJrjyQBQ\nXBh1M3LRtApuv+YNvL7nULyu408vvc4v1zRiBmdNHB1VlM+I8aZJZRRqmFmRnMrkpvBbOT4xXN3F\nsmRd3TNIvozpbpttRONn7DCzc4DfmNkZ7r7nmJ3N5gPzASorK9OEI/3R+FElXF87metrJ9Pe4Tzb\nuDve9fo3H2ng3/7cwKiSQi6ZXsGc6VHz3FPLhvV12CIDTqo6iw8BtwBTzezZhFUjgb9k8N6NwOSE\n+UlAcmuqzm0aw22o0cBOj+6NHQZw99Vmth6oAY65z+TuC4GFEN2GyiAm6ccKhhizK8uZXVnOxy+f\nTsuBVh5riBLH8rpmfvfcawDUjB8Rf67j3OoxGmZWJAtS1VmMBsqBLwO3Jaza6+47075x9OVfB7wF\neBV4CnhXYnchZvZh4I3u/sFQwf23oVI9RpQ02s1sKvBo2K7b46rOYnBzd+pe3xdPHE9u3MmR9g5K\nhg7hgqlj48ljasVwVZSLJMh601kzGwfEH7919y0Z7HMN8HWiprM/cPe7zOxOYJW7LzazEuBeYDaw\nE7jR3TeY2TuBO4E2oB34nLv/NtWxlCwk0YEjbazcsINl65pYXt/MxuZomNlJ5cOi3nOnx7h42lhG\nlmiYWRncsj2s6r8ApwLbiVoqveTuZ2Qj0GxRspBUtuw4wLL6qKJ8RUMz+4+0UzjEOLuqPP5E+cwJ\nozTMrAw62UwWzxANq/ond59tZm8GbnL3+dkJNTuULCRTR9o6WLNlF8vrogGfXtwWtZuoGFHEpeGh\nwEunVzB2RHEfRyqSe9lMFqvcvTYkjdnu3mFmT7r7edkKNhuULOREbd97iMfqm1lW18Sj9c3s3H8E\nMzjz1NHxuo7ZlWUMVfNcGYCymSz+BLyDqKK7guhW1LnuflE2As0WJQvJho4O5/mtLfEOENds2U17\nhzOyuJCLpo1lbs045tRUMKlcw8zKwJDNZDEcOET0TMQ8ouati9x9RzYCzRYlC8mFloOtPL4+uupY\ntq6JrS3RaMKnxYbHE8cFU8eqea70W3nRkWBvUrKQXHN31jftY1ldlDye2LCDw20dFBcO4bwpY5hb\nE+OyGTFOi41Q81zpN046WZjZXrroQLDzp7uPykag2aJkIb3tUGs7T2zcydJ121le18T6pqh57sSy\nYcypiZ4ov2haBaOHqXmu5C9dWYj0ssZdB+Jdkaxo2MHew23RU+edw8zOiHHmqaPVPFfySrbH4L4E\nmO7uPzSzCmCku2/MQpxZo2Qh+aS1vYOnt+yOV5Q/92oLAGOGF3Fp6Mfq0poKxo3UMLPSt7JZwf05\noBaY4e41ZnYq0ah5F2cn1OxQspB81rzvMI/VR8PMLq9vonnfEQBmThgVH7PjnKpyigrVPFd6VzaT\nxVqi7jjWuPvssOzZdIMf9TYlC+kvOjqcF7ftiY/ZsWbzLto6nOFFBVx4WgVzZ8S4rCbG5DFqniu5\nd9LjWSQ44u5uZh7eWAMli5yEIUOMMyeO5syJo/nwm6ex91ArK8Iws53jdgBMqRgeHgqMmueWFmmY\nWek7mfz23W9m3wXKzOwDwPuA7+c2LJHBY2TJUK484xSuPOMU3J0NzfvjdR0/e2oL96zYRFHB0ea5\nc2pi1IxX89wTtWgRLFgAW7ZAZSXcdRfMm9fXUeW/TCu43wpcQdRsdom7P5TrwHpKt6FkIDrU2s5T\nm3bGk0fd6/sAOGVUSdQ8tybGJdMqKCvVMLOZWLQI5s+HAweOListhYULB2/CyFnTWTMrIOpKfNGJ\nBpcLShYyGGxrORhPHI/VN7PnUBtDDGZNLos/UX7WpDIK1Dy3S9XVsHnz8curqmDTpt6OJj9k46G8\nUcCHicbJXgw8FOY/Bax197dnL9yTp2Qhg01bewfPNO5m2bomltU382zjbtyhrHQol0yriHe9Pm6U\nmud2GjIEuvrKM4OOjt6PJx9kI1n8F7ALeJxotLtyoAj4uLuvzWKsWaFkIYPdzv1HeKyhOX7l0bT3\nMACnnzIyXtdRW11OceHg7cdKVxbHy0ayeM7d3ximC4BmoNLd92Y10ixRshA5yt15adteltdHHSCu\n2ryT1nZn2NACLjptLHPCVUd1xeBq3Kg6i+Nlo+lsa+dEGAt7Y74mChE5lpkx89RRzDx1FB+cexr7\nD7fx+PodLK+PBnx6+OXtAFSNLWVOGPDpwtPGMrx4YDfP7UwIag3Vc6muLNqB/Z2zwDDgAOpIUKTf\n29S8P37V8fiGHRw40g7txqHGMZTsivGR62J84v0j1Tx3EFBHgiKSkR/d287Hv7gLO7WJYVOaKBoX\n3UAYWVjMFWdFHSBeOq2C8uFqnjsQKVmISEaSK30LRhyiZEoTFWc0Mfr0ZloOtmIGZ00KvefWVPCm\nSWUUapjZAUHJQkQykqo5aWub80xj1Hvu8rom1r6ymw6HUSWFXDK9It7KasLoYb0fuGSFkoWIZKQn\nzUl3Hzi2ee7re6LmuTXjR8QTx7nVYzIaZlbdbuQHJQuRLBuoX24n2pzU3al7fR/L6razrK6Jpzbu\n4kh7ByVDh3DB1LHx5DG1YvhxFeVqwpo/lCxEsmigf7llIxEeONLGExvCMLP1zWxsjhpTTiofFk8c\nF502lpElQ/VwXB5RshDJIn259dyWHQei5rl1TaxoaGb/kXYKhxhnV5Xzx3tiHNwY48jro4ha40cG\nc7cbfUXJQiSL1KfQyTnS1sHqzbviz3a8uG0PAO37izi4McahjTEObqxg8rhiJd9epmQhkkW6ssiu\n79xziH/+92YKJjZRUt1EQWkr7jCxdDR/e0F0y2p2ZRlD1Tw357I5Up7IoHfXXV3XWdx1V9/F1J99\n8L0ljBw6iQULJrHlv53KN7Vw+bubaBnWxLeXreebjzQwsriQi6dFw8zOqYkxsUzNc/uSrixEMjRQ\nW0Plm5aDraxoaI7fstracgiAaeNGRP1YzYhx/pTMmudKeroNJSL9nrvTsH1fND55fTMrN+zgSFsH\nxYVHh5m9bEaM02IaZvZEKVmIyIBz8Eg7T2zcwfK6ZpbVbWd9U9Q8d2LZsGiY2ekxLp5ewaiSoSd8\njMF2BalkkSWD7RdHpD9p3HUgnjhWNOxg7+E2CoYYZ1eWxZ/tOPPU0QzJcJjZgf48TVeULLJgMP7i\niPRXre0dPL1ld7wrkudebQFg7PCieD9Wl06PERtZzC23RH/H7e1R8+fhw2H//qiJdHv78e89kFu9\nKVlkgZpLSiddYfY/zfsO81h91I/V8vommvcdAWBU+yheeTJ6KPDwq+XQkb557kB+nibTZKFGzCls\n2dKz5TIwdV5hbt4cPZi3eXM0v2hRX0cmEH0O1dXRVUF19dHPpWJEMe+YPZF/uWEWT/7z5Tz40Uv4\n1JUzaNpayKjzNnDKu1Yy+WMPEfubVYyYtZnC0Qe6PUZlZffHGSxyemVhZlcB3wAKgO+7+1eS1hcD\nPwbOAXYAN7j7prDuduD9QDvwMXdfkupYurKQXNHvQf7q6laxWZTUq6qOXgEmXhm6gxW1UlK1g2FT\nogGfCssOAtC6YzgHN4arjlfG4K2FlJbCzTfDj3507HEAxo6Fb3yjf19lZnplgbvn5EWUINYDU4Ei\n4BlgZtI2twDfCdM3Aj8P0zPD9sXAlPA+BamOd84553i23Xefe2mpe/TrFb1KS6PlMniYHfs70Pky\n6+vI8tN997lXVUXnp6oqt38vVVVdfzaJf68f+tDxf8fHvjq8cMxeH3nOBh933RM++Z9+51W3PuiV\nn/idj79hpV/58fVeddYeh45u32Ps2OyXs7fOI7DKM/lOz2SjE3kBFwJLEuZvB25P2mYJcGGYLgSa\niXoVO2bbxO26e+UiWbj37i++5KfuvpCqqvo6svzT2/9gdZfIE18FBem3OeZV0OYl1du97M0v+IT3\nLfWqWx/0qlsf9Im3/MnHXr3WS2ds9SHFR7pMTNkqZ2+ex0yTRc5uQ5nZdcBV7v4PYf49wPnu/pGE\nbZ4P2zSG+fXA+cDngZXufl9Y/p/A7939F90dT89ZSK6oVVzmevuWXXfHy6aCkQcpndpEcXUTw6qb\nGVLShnfA4W1lHNowjoMbKzjyWhm4Za2cvXke86FvqK4aNidnpu62yWRfzGw+MB+gsrKyp/GJZKQz\nIag1VHq93Sikqz67khUUdN0cNlPte4ex95lK2usraT7YQfGpuymZ0sSwKc2MvqSOskvraD84lEOb\nKtixMcbre2KMH1Vy4gckPxvX5LI1VCMwOWF+ErC1u23MrBAYDezMcF/cfaG717p7bSwWy2LoIsea\nNy/6j66jI/qpRNG17v5nO5n/5VK1Qpo3L7rCq6qK5pN7/CgtjZJJaWn379+5bypVVdFxxo4ZwuFX\nx9Dy2Axeu/diGv/9rTQtnsXBhvEUT97J2Gue5fwvPcxVX1/Ol3//Eisamjnc1vNMlYvzeNIyuVd1\nIi+iq5YNRBXUnRXcZyRt82GOreC+P0yfwbEV3BvogwpuEemZbN9r7+n7dVfHeN99USV0d/UMb3lL\n93UYZscer/v36vCvfa/Fv720wW/87uM+7Z//26tufdDf8Nnf+/t++KT/aMVG39i0r0/OYyr0dQV3\nFAPXAHVErZkWhGV3AteG6RLgAaABeBKYmrDvgrDfOuDqdMdSshDJD9lsFJLtxgWpYps5s+tE8aEP\n9fy93N33Hmr1P77wmn/m18/5pV/9c7yifM7df/bP/uY5f+iF13zfodYTijWbMk0WeoJbRPJWb49Q\nmMsn9Tc172dZ6Irk8fU7ONjaztACo7ZqTDRmx/QYb5gwstd7z1V3HyLS7w3UByIPt7WzetMultY1\nsbyuiZdf2wvAuJHFzAkdIF46rYLy4UU5j0XJQkT6vcHSbPm1lkMsr48Sx6P1zbQcbMUM3jSpjDk1\nMT5w6RRGnkS366koWYjIgDDYOnFs73CeaTzae27D6/tY/dm3UlSYm8arShYiIgPAwSPtDCvK3RCy\n6nVWRGQAyGWi6AklCxERSUvJQkRE0lKyEBGRtJQsREQkLSULERFJS8lCRETSUrIQEZG0lCxERCQt\nJQsREUlLyUJERNJSshARkbT7xyZsAAAKq0lEQVSULEREJC0lCxERSUvJQkRE0lKyEBGRtJQsREQk\nrQEzUp6ZNQFdDO2ecxVAcx8cN5tUhvygMuSH/l6GnsZf5e6xdBsNmGTRV8xsVSZDEuYzlSE/qAz5\nob+XIVfx6zaUiIikpWQhIiJpKVmcvIV9HUAWqAz5QWXID/29DDmJX3UWIiKSlq4sREQkLSWLDJjZ\nGDN7yMzqw8/yLraZZWaPm9kLZvasmd2QsO4eM9toZmvDa1bvliArZZhiZk+E/X9uZkW9W4LMyhC2\n+4OZ7TazB5OW94vPIWzXXRn60+dwc9im3sxuTli+1MzWJXwO43op7qvCcRvM7LYu1heHc9oQznF1\nwrrbw/J1ZnZlb8TblRMtg5lVm9nBhHP+nR4f3N31SvMC7gZuC9O3AV/tYpsaYHqYPhXYBpSF+XuA\n6/p5Ge4HbgzT3wE+lI9lCOveArwNeDBpeb/4HNKUoV98DsAYYEP4WR6my8O6pUBtL8dcAKwHpgJF\nwDPAzKRtbgG+E6ZvBH4epmeG7YuBKeF9CvrgvJ9MGaqB50/q+L1d4P74AtYBE8L0BGBdBvs8k/DF\nmw9fUidcBsCIHvIpDMsvBJbkcxmAy/I0WZxwGfrT5wDcBHw3Yf67wE1hui+SxTHnCrgduD1pmyXA\nhWG6MJxrS942cbt+VIaTTha6DZWZ8e6+DSD8THnZbGbnEWX+9QmL7wq3dv7VzIpzF2q3TqYMY4Hd\n7t4WVjcCE3MYa3d6VIZu9KvPIUl/+hwmAq8kzCfH+sNwO+SzZma5CzXjeI7ZJpzjFqJznsm+veFk\nygAwxcyeNrNlZnZpTw9e2PN4ByYz+xNwSherFvTwfSYA9wI3u3tHWHw78BrRl+9C4FbgzhOPtttj\n56QM3fwx56QZXbbK0I1+9Tl09dZdLMvXzyFVrPPc/VUzGwn8EngP8OOeR9kjmZy77rbptfOexsmU\nYRtQ6e47zOwc4Ddmdoa778n04EoWgbtf3t06M3vdzCa4+7bwRbq9m+1GAf8NfMbdVya897YwedjM\nfgh8Mouhx+WwDM1AmZkVhv9WJgFbsxw+kJ0ypHjvfvM5dKM/fQ6NRLfROk0iuv2Eu78afu41s58A\n55H7ZNEITE6KJ/ncdW7TaGaFwGhgZ4b79oYTLoNH96UOA7j7ajNbT1RHuSrTg+s2VGYWA52tOW4G\n/it5g9Aq5dfAj939gaR1E8JPA94BPJ/TaLt2wmUIv2iPANel2r8XpC1DKv3lc+hOP/sclgBXmFl5\naC11BbDEzArNrALAzIYCf03vfA5PAdNDa7IiosrfxUnbJJbrOuDP4ZwvBm4MLY2mENXjPdkLMSc7\n4TKYWczMCgDMbCpRGTb06Oi9XUnTH19E9/weBurDzzFheS3w/TD9bqAVWJvwmhXW/Rl4juiP4j5g\nRD8sw1SiP5AG4AGgOB/LEOYfBZqAg0T/aV3Znz6HNGXoT5/D+0KcDcDfh2XDgdXAs8ALwDfopZZF\nwDVAHVE93IKw7E7g2jBdEs5pQzjHUxP2XRD2Wwdc3dvn/GTLALwznO9ngDXA23p6bD3BLSIiaek2\nlIiIpKVkISIiaSlZiIhIWkoWIiKSlpKFiIikpWQxyJhZe+hm4Xkz+62ZlWWwz74068vM7JaE+VPN\n7BfZiLeb4y01sz4fIznX5QzHmGVm12T5Pfvk/GValuTtzOzarnpYld6lZDH4HHT3We5+JtHTqR/O\nwnuWEfV2CYC7b3X361JsPyD0UjlnEbWtz5hF8vFvO9OyHLOduy9296/kLCrJSD7+QknveZyEjsjM\n7FNm9lToaO+O5I3NbISZPWxma8zsOTN7e1j1FeC0cMXytdB3/vNhnyfM7IyE91hqZueY2XAz+0E4\n3tMJ75V8zE+HYz1jZolfGH9nZk+aWV1np2jhuI+G+NaY2UVh+WXhuL8ws5fNbFFnf1dmdk1Y9piZ\n/ZuF8SMyiS+pnO81s19ZNA5FvZnd3U15NpnZHQnn8PSw/DwzWxGOtcLMZoSndO8Ebgjn9gYz+7yZ\nfTLh/Z4PcVSb2Utm9h9ED11NNrNvm9kqi8YnOe7zTIrrajO7P2H+MjP7bZi+wqJxTtaY2QNmNqKn\n566bsmRa5vea2TfDe1eF38Fnw8/KsPyeEMMKM9tgZgP+n5Ve11dPIurVNy9gX/hZQPSk51Vh/gqi\nzvWM6J+IB4E5SfsUAqPCdAXRU6LHdX+cOA/8I3BHmJ4A1IXpLwHvDtNlRE+lDk+K9WpgBVAa5juf\nFF4K/L8wfQ3wpzBdCpSE6enAqjB9GVHvm5NC2R4HLiF62vUVYErY7qeELsEzjC+xnO8l6j5hdHjf\nzcDkLs7/JuCjYfoWjj49P4qjXY9fDvwy4X2/mbD/54FPJsw/H+KoBjqACxLWdZ6vgnDOzko4f7VJ\ncRUCWzrLCHyb6In+CmB5wvJbgf9zIueui7JkWub4PPBbog4uIXpC/Ddh+h6i3+chRONPNPT139pA\ne+nKYvAZZmZrgR1EA9M8FJZfEV5PE/1nejrRF24iA75kZs8CfyK6Khmf5nj3A38Xpq8n+oPuPN5t\nIZalRF8+lUn7Xg780N0PALj7zoR1vwo/VxN9UQIMBb5nZs+F48xM2P5Jd2/0qCfgtWGf04EN7r4x\nbPPThO0ziS/Zw+7e4u6HgBeBqm626yr20cAD4UrlX4Ezutgvnc2e0IElcL2ZrSH6TM/g2PNxDI86\nJvwD8DaLOqD7H0R9Pl0Q9vtLOBc3h3Jl49ydSJkvBH4Spu8lSvqdfuPuHe7+Iul/L6WH1Ovs4HPQ\n3WeZ2Wiiq4cPA/9GlAi+7O7fTbHvPCAGnOPurWa2ieiLoFsedUW9w8zOAm4A/ldYZcA73X1dit2N\n7ruCPhx+tnP09/gfgdeBNxH9h3moi+0T90k1jkIm8XUXU3Jc3W2XuM0XgEfc/W8sGgpzaTf7tnHs\n7ePE87+/c8KiDu8+CZzr7rvM7B7SfFbAz4l+H3YCT3nUK6wBD7n7TYkbmtnsFO/T5bkzs/OTtsu0\nzKkk/n4knv/eGCNjUNGVxSDl7i3Ax4BPWtT75xLgfQn3oyfa8WMjjwa2h0TxZo7+57wXGJnicD8D\nPg2MdvfnwrIlwEcT6g66+vL5Y4ipNGwzJk2xRgPbwtXDe4huv6TyMjDVjo61fEPCukziy6bRwKth\n+r0Jy5PP7Sbg7BDT2UTDfHZlFFHyaDGz8US39NJZGt77A0SJA2AlcLGZTQvHLDWzGk7s3CWXJdMy\nJ1pB1NsqRP+8PJZBuSQLlCwGMXd/mqgXyhvd/Y9El/ePh9s4v+D4P9hFQK2ZrSL6Q305vM8OotsU\nz5vZ17o41C+I/sDvT1j2BaLbRs+G2xBf6CK+PxB1ubwq3NJIN/7EfwA3m9lKor7696fa2N0PEtUb\n/MHMHiO6KmnJNL4suxv4spn9hWOT3CPAzM7KXqLBgsaE8/EhovqA47j7M0S3n14AfgD8JV0A7t5O\ndLV5dfiJuzcRfZH/NNx+XAmcfoLnLrksmZY50ceAvw+xvAf4eLpySXao11kZ1MxshLvvC/8Ffwuo\nd/d/7eu4+gOdu8FFVxYy2H0g/Jf+AtFtkVR1NnIsnbtBRFcWIiKSlq4sREQkLSULERFJS8lCRETS\nUrIQEZG0lCxERCQtJQsREUnr/wNRcoVuGfXYxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "coeffs= do_regress(ecoli, nature1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the coefficient lines in the regression output:\n",
    "\n",
    ">     ==============================================================================\n",
    ">                    coef    std err          t      P>|t|      [0.025      0.975]\n",
    ">     ------------------------------------------------------------------------------\n",
    ">     const          0.0033      0.007      0.492      0.627      -0.011       0.017\n",
    ">     x1            -0.3036      0.091     -3.328      0.003      -0.491      -0.116\n",
    "\n",
    "We are interested in the slope, x1, which is summarized in the bottom row.\n",
    "\n",
    "The least-squares estimate of x1 is coef=-0.3036. \n",
    "This value can be calculated provided the design matrix $X$ has full rank (so that $X^TX$ is invertible). \n",
    "The estimate is unbiased if the model $Y = X\\beta + \\epsilon$ is true, with $\\mathbb{E} \\epsilon = 0$.\n",
    "\n",
    "The standard error 0.091 can be calculated provided $X$ has full rank and $n > p$ (the number of data is greater than the number of parameters). But it is \"really\" an estimate of the standard error of the estimated slope only if the model is true, $X$ is full rank, $n > p$, $\\mathbb{E}\\epsilon = 0$, and $\\mathrm{cov} \\epsilon = \\sigma^2 I_{n}$.\n",
    "\n",
    "The t statistic -3.328 can be computed whenever coef and std err can be computed.\n",
    "But the p-value, P>|t|, is correct only if the model is true, $X$ has full rank, $n > p$, and $\\epsilon \\sim N(0, \\sigma^2 I_n)$. \n",
    "The same is true for the confidence interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Significance of the slope\n",
    "\n",
    "### Standard regression model\n",
    "\n",
    "The \"usual\" regression model says $Y = X\\beta + \\epsilon$, where \n",
    "\n",
    "\\begin{equation*}  X = \\left [ \\begin{array}{ccccc} X_{11} & X_{12} & \\cdots & X_{1p} \\\\ \n",
    "                X_{21} & X_{22} & \\cdots & X_{2p} \\\\\n",
    "                \\vdots & \\vdots & \\vdots & \\vdots \\\\ \n",
    "                X_{n1} & X_{n2} & \\cdots & X_{np} \\end{array} \\right ]\\end{equation*}\n",
    "\n",
    "is the _design matrix_, $\\beta = [\\beta_1 \\beta_2 \\cdots \\beta_p]^T$ are the _coefficients_, $\\epsilon = [\\epsilon_1 , \\ldots, \\epsilon_n]^T$ is the _noise_ or _random error_, and it is assumed that\n",
    "$\\{ \\epsilon_j \\}_{j=1}^n$ are IID with mean zero and variance $\\sigma^2$.\n",
    "We assume that $X$ is full rank $p$ and either that $X$ is not random or that $X$ and $\\epsilon$ are independent, in which case we will condition on $X$. \n",
    "We also assume $n > p$.\n",
    "\n",
    "The OLS estimate is $\\hat{\\beta} \\equiv (X^TX)^{-1}X^T Y$.\n",
    "\n",
    "#### Theorem: the OLS estimate is unbiased and $\\mathrm{cov}(\\hat{\\beta}) = (X^TX)^{-1} X^T \\mathrm{cov}(\\epsilon | X) X (X^TX)^{-1}$.\n",
    "\n",
    "**Proof.** \n",
    "\n",
    "\\begin{equation*}  (X^TX)^{-1}X^T Y = (X^TX)^{-1}X^T (X\\beta + \\epsilon) = (X^TX)^{-1}X^TX\\beta + (X^TX)^{-1}X^T\\epsilon \\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = \\beta + (X^TX)^{-1}X^T\\epsilon.\\end{equation*}\n",
    "\n",
    "Now, \n",
    "\n",
    "\\begin{equation*} \\mathbb{E} \\left (\\beta + (X^TX)^{-1}X^T\\epsilon | X \\right ) = \\beta + (X^TX)^{-1}X^T \\mathbb{E} (\\epsilon | X)\\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = \\beta + (X^TX)^{-1}X^T \\mathbf{0} = \\beta,\\end{equation*}\n",
    "\n",
    "so $\\hat{\\beta}$ is unbiased.\n",
    "\n",
    "\\begin{equation*}  \\mathrm{cov} \\left (\\hat{\\beta} | X \\right ) \\equiv \\mathbb{E} ((\\hat{\\beta} - \\beta)(\\hat{\\beta} - \\beta)^T | X) \\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = \\mathbb{E} \\left ( (X^TX)^{-1}X^T \\epsilon ((X^TX)^{-1}X^T \\epsilon)^T | X \\right ) \\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = \\mathbb{E} \\left ( (X^TX)^{-1}X^T \\epsilon \\epsilon^T X (X^TX)^{-1} | X \\right ) \\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = (X^TX)^{-1}X^T \\mathbb{E}(\\epsilon \\epsilon^T | X) X (X^TX)^{-1}\\end{equation*}\n",
    "\n",
    "\\begin{equation*}  = (X^TX)^{-1}X^T \\mathrm{cov}(\\epsilon | X) X (X^TX)^{-1}.\\end{equation*}\n",
    "\n",
    "Since $\\mathrm{cov}(\\epsilon | X) = \\sigma^2 I_{n}$, this simplifies to\n",
    "\n",
    "\\begin{equation*}  \\mathrm{cov} (\\hat{\\beta} | X) = \\sigma^2 (X^TX)^{-1}X^T I_n X (X^TX)^{-1} = \\sigma^2 (X^TX)^{-1}.\\end{equation*}\n",
    "\n",
    "Define $e = Y - X \\hat{\\beta}$, the _residuals_.\n",
    "\n",
    "A conditionally unbiased estimate for $\\sigma^2$ is \n",
    "\n",
    "\\begin{equation*}  \\hat{\\sigma}^2 = \\frac{1}{n-p} \\sum_{j=1}^n e_j^2 = \\frac{1}{n-p} \\sum_{j=1}^n (Y_j - X_j \\hat{\\beta})^2.\\end{equation*}\n",
    "\n",
    "#### The t-test\n",
    "\n",
    "The $t$-test for whether a coefficient is significantly different from zero requires the further assumption that\n",
    "$\\{ \\epsilon_j \\}$ have a distribution that is jointly normal (given $X$):  given $X$, $\\{ \\epsilon_j \\}$ are IID $N(0, \\sigma^2)$.\n",
    "\n",
    "Then $\\hat{\\beta} \\sim N(\\beta, \\sigma^2 (X^TX)^{-1})$.\n",
    "\n",
    "The $t$ statistic for the $\\hat{\\beta}_k$ is \n",
    "\n",
    "\\begin{equation*}  t = \\hat{\\beta}_k/\\hat{\\mathrm{SE}}_k, \\end{equation*}\n",
    "\n",
    "where $\\hat{\\mathrm{SE}}_k \\equiv \\hat{\\sigma} \\sqrt{(X^TX)^{-1}_{kk}}$.\n",
    "\n",
    "Under all the assumptions, $t$ has Student's $t$ distribution with $n-p$ degrees of freedom if $\\beta_k = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Permutation test\n",
    "\n",
    "Because the assumptions that lead to the t-test are unreasonable in this application, we will look at other ways of assessing the evidence that the two variables are associated.\n",
    "\n",
    "Let's specialize to the case $p = 2$, with $X_{k1} = 1$ (an intercept term in the regression).\n",
    "\n",
    "Suppose that the components of $\\epsilon$ are IID (but not necessarily zero mean, normal, or anything else). Under the hypothesis that $\\beta_1 = 0$, the data are then\n",
    "\n",
    "\\begin{equation*}  Y_j = \\beta_0 + 0\\times X_j + \\epsilon_j = \\beta_0 + \\epsilon_j, \\;\\; j=1, \\ldots, n.\\end{equation*}\n",
    "\n",
    "Since $\\{ \\epsilon_j \\}$ are IID, so are $\\{ Y_j = \\beta_0 + \\epsilon_j \\}$.\n",
    "\n",
    "__Definition.__ A set of random variables $\\{Z_j \\}_{j=1}^n$ is _exchangeable_ if the joint distribution of $(Z_{\\pi_1}, Z_{\\pi_2}, \\ldots, Z_{\\pi_n})$ is the same for all permutations $\\pi$ of $\\{1, 2, \\ldots, n\\}$.\n",
    "\n",
    "If a set of variables is IID, it is exchangeable. (The converse is not necessarily true.)\n",
    "\n",
    "Thus if $\\beta_1 = 0$, $\\{Y_j \\}_{j=1}^n$ are exchangeable.\n",
    "\n",
    "Suppose $\\{Z_j \\}_{j=1}^n$ are exchangeable, and let $z \\in \\Re^n$. \n",
    "Then for every permutation $\\pi$ of $\\{1, \\ldots, n \\}$,\n",
    "\n",
    "\\begin{equation*} \\Pr \\{ Z_j = z_j, j=1, \\ldots, n \\} = \\Pr \\{ Z_j = z_{\\pi_j}, j = 1, \\ldots, n \\}.\\end{equation*}\n",
    "\n",
    "Let $z = (z_1, \\ldots, z_n)$.\n",
    "The set of all $(z_{\\pi_j})_{j=1}^n$ for all permutations $\\pi$ of $\\{1, \\ldots, n \\}$ is _the orbit of $z$ under the permutation group_.\n",
    "\n",
    "#### What's the justification for *this* model here?\n",
    "\n",
    "Why linear? Why additive errors? Why IID errors?\n",
    "\n",
    "#### Another approach: the Neyman \"ticket\" model\n",
    "\n",
    "| item | response if assigned $X_1$ | $\\cdots$ | response if assigned $X_n$ |\n",
    "|:----:|:--------------------------:|:--------:|:--------------------------:|\n",
    "| 1    |  $Y_{11}$                  | $\\cdots$ | $Y_{1n}$                  |\n",
    "|  $\\vdots$ | $\\vdots$              | $\\vdots$ |  $\\vdots$                 |\n",
    "| $n$  |  $Y_{n1}$                  | $\\cdots$ | $Y_{nn}$                  |\n",
    "\n",
    "Suppose there's no connection between $X$ and $Y$: response of item $j$ would be the same, regardless of\n",
    "the value of $X_i$.\n",
    "Then all \"shufflings\" of $Y$ versus $X$ would be equally likely;\n",
    "treat the data as if the treatment were assigned at random to the farm.\n",
    "\n",
    "Is that more plausible or less plausible than the regression models?\n",
    "\n",
    "Leads to exactly the same permutation test as $Y = X\\beta + \\epsilon$ with \n",
    "the components of $\\epsilon$ IID.\n",
    "\n",
    "Why use $|\\hat{\\beta}_1|$ as the test statistic?\n",
    "\n",
    "Equivalent to a permutation test based on the Pearson correlation:\n",
    "\n",
    "$\\hat{\\beta}_1 = r \\mathrm{SD}(Y)/\\mathrm{SD}(X)$, where $r$ is the (sample) correlation between $X$ and $Y$.\n",
    "\n",
    "One-to-one, monotonic relationship between them, so reject the null for exactly the\n",
    "same data sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to shuffle\n",
    "\n",
    "1. Assign a pseudo-random number to every element, then sort based on that number\n",
    "\n",
    "1. Fisher-Yates-Knuth-Durstenfeld shuffle\n",
    "\n",
    "**Algorithm Fisher-Yates-Knuth-Durstenfeld shuffle (backwards version)**\n",
    "\n",
    "    for i=1, ..., n-1:\n",
    "        J <- random integer uniformly distributed on {i, ..., n}\n",
    "        (a[J], a[i]) <- (a[i], a[J])\n",
    "        \n",
    "This requires the ability to generate independent random integers on various ranges. Doesn't require sorting.\n",
    "\n",
    "There's a version suitable for streaming, i.e., generating a random permutation of a list that has an (initially) unknown number  nn of elements:\n",
    "\n",
    "**Algorithm Fisher-Yates-Knuth-Durstenfeld shuffle (streaming version)**\n",
    "\n",
    "    i <- 0\n",
    "    a = []\n",
    "    while there are records left:\n",
    "        i <- i+1\n",
    "        J <- random integer uniformly distributed on {1, ..., i}\n",
    "        if J < i:\n",
    "            a[i] <- a[J]\n",
    "            a[J] <- next record\n",
    "        else:\n",
    "            a[i] <- next record\n",
    "\n",
    "Again, need to be able to generate independent uniformly distributed random integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fykd(a, gen=np.random):\n",
    "    for i in range(len(a)-1):\n",
    "        J = gen.randint(i,len(a))\n",
    "        a[i], a[J] = a[J], a[i]\n",
    "    return(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof that the streaming Fisher-Yates-Knuth-Durstenfeld algorithm works\n",
    "\n",
    "Induction: \n",
    "\n",
    "For $i=1$, obvious.\n",
    "\n",
    "At stage $i$, suppose all $i!$ permutations are equally likely. For each such permutation, there are\n",
    "$i+1$ distinct, equally likely permutations at stage $i+1$, resulting from swapping the $i+1$st item with one\n",
    "of the first $i$, or putting it in position $i+1$. These possibilities are mutually exclusive of\n",
    "the permutations attainable from a different permutation of the first $i$ items.\n",
    "\n",
    "Thus, this enumerates all $(i+1)i! = (i+1)!$ permutations of $(i+1)$ items, and they are equally likely. &#x25A0;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How not to shuffle\n",
    "\n",
    "Sort using a \"random\" comparison function, e.g.,\n",
    "\n",
    "    def RandomSort(a,b):\n",
    "        return (0.5 - np.random.random())\n",
    "\n",
    "But this fails the basic properties of an ordering, e.g., transitivity, reflexiveness, etc. Doesn't produce random permutation. Output also depends on sorting algorithm. \n",
    "\n",
    "Notoriously used by Microsoft to offer a selection of browsers in the EU. Resulted in IE being 5th of 5 about half the time in IE and about 26% of the time in Chrome, but only 6% of the time in Safari (4th about 40% of the time).\n",
    "\n",
    "See, e.g., http://www.computerworld.com/article/2520190/web-apps/microsoft-s-eu-ballot-fails-to-randomize-browser-order.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to generate random integers\n",
    "\n",
    "Textbook approach:\n",
    "\n",
    "To generate a random integer between 0 and m-1\n",
    "\n",
    "* Generate $U \\sim U[0, 1)$\n",
    "* Return $\\lfloor mU \\rfloor$\n",
    "\n",
    "Works if $U$ _really were_ $U[0, 1)$, but that's not what PRNGs generate. Instead, typically generate a $k$-bit binary fraction (or its floating-point representation).\n",
    "\n",
    "Introduces truncation effects. Even if the PRNG were perfectly uniform on its output, result would be uniformly distributed on $\\{0, 1, \\ldots, m-1\\}$ only if $m$ is a power of 2 (not exceeding the number of bits of output of the PRNG).\n",
    "\n",
    "Better approach (used in python):\n",
    "\n",
    "Let $k = \\lceil \\log_2(m-1) \\rceil$. Generate $k$ random bits and interpret the result $M$ as an integer. If $M \\le m-1$, return $M$. Otherwise, generate $k$ more random bits and try again.\n",
    "\n",
    "See [cryptorandom](https://github.com/statlab/cryptorandom/blob/master/cryptorandom/cryptorandom.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How R generates random integers\n",
    "\n",
    "R uses a variant of the floor approach, but uses floating-point numbers at an intermediate step. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use cryptorandom's permutation\n",
    "def permute_regress(ec, nat, reps=int(10**5)):\n",
    "    '''\n",
    "    Permutation distribution of OLS estimate of slope, regressing ec on nat\n",
    "    '''\n",
    "    vals = np.zeros(reps)                                         # allocate storage for slope estimates\n",
    "    for i in range(reps):\n",
    "        permutation = np.array(sample.fykd_sample(len(ec), len(ec), gen=prng))-1  # shuffle indices\n",
    "        p_ec = ec[permutation]                                    # permute the E coli values\n",
    "        vals[i] = sm.OLS(p_ec, nat).fit().params[1]               # regress\n",
    "    return vals\n",
    "\n",
    "def p_CI(v,t): \n",
    "    '''\n",
    "    estimated probability that |v| >= t based on independent trials,\n",
    "    and upper confidence bound for the true probability\n",
    "    '''\n",
    "    exceed = sum(abs(v) >= t)\n",
    "    print(\"estimated p:\", exceed/len(v), binom_conf_interval(n=len(v),x=exceed, alternative='upper'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated p: 0.03099 (0.0, 0.032082742274273775)\n"
     ]
    }
   ],
   "source": [
    "vals = permute_regress(ecoli, nature1)\n",
    "p_CI(vals, abs(coeffs[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity to farm 31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecoli_edit = np.array([float(f['E_Coli_Change']) for f in farms if f['Farm ID'] != '31'])     # omit farm 31\n",
    "nature_edit = np.array([float(f['Change in Nature']) for f in farms if f['Farm ID'] != '31']) # omit farm 31\n",
    "nature_edit1 = sm.add_constant(nature_edit)    # add intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.059\n",
      "Model:                            OLS   Adj. R-squared:                  0.021\n",
      "Method:                 Least Squares   F-statistic:                     1.562\n",
      "Date:                Sun, 16 Sep 2018   Prob (F-statistic):              0.223\n",
      "Time:                        20:37:57   Log-Likelihood:                 75.457\n",
      "No. Observations:                  27   AIC:                            -146.9\n",
      "Df Residuals:                      25   BIC:                            -144.3\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.0054      0.003      1.645      0.113      -0.001       0.012\n",
      "x1            -0.0635      0.051     -1.250      0.223      -0.168       0.041\n",
      "==============================================================================\n",
      "Omnibus:                       45.099   Durbin-Watson:                   2.175\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              188.313\n",
      "Skew:                           3.230   Prob(JB):                     1.28e-41\n",
      "Kurtosis:                      14.210   Cond. No.                         17.2\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmUHOV57/Hvb2YkzYyWkZBYJUYS\nQYYLMWEZY/vaybUTg4EbjO8xMcKyr7wkJF6zeYGjJAYcHOMkJvYBYiu2MQbFbHFy5dhBxix2DBg0\ngFiEwQiMQAsgoQ2Y0f7cP6pa6ml195RmqrunR7/POX26lre7nqqe6affeqveVxGBmZnZcLU0OgAz\nMxsdnFDMzCwXTihmZpYLJxQzM8uFE4qZmeXCCcXMzHLhhGJmZrlwQjEzs1w4oZiZWS7aGh1APU2b\nNi1mzZrV6DDMzJrGtGnTWLJkyZKIOGOwsgdUQpk1axa9vb2NDsPMrKlImpalnE95mZlZLpxQzMws\nF04oZmaWCycUMzPLhROKmZnlwgnFzMxy4YRiZma5cEIxM7NcOKGYmVkunFDMzCwXTihmZpYLJxQz\nM8uFE4qZmeXCCcXMzHLhhGJmZrlwQjEzs1w4oZiZWS6cUMzMLBdOKGZmlgsnFDMzy4UTipmZ5aKh\nCUXSGZKelLRC0oVl1o+TdGO6/j5Js9LlYyRdK+lRSb+UdFG9Yzczs4EallAktQJXAWcCxwHnSzqu\npNhHgI0RcTRwBXB5uvwPgHER8XrgFOCPC8nGzMwao5E1lFOBFRHxTERsB24Azikpcw5wbTp9C/B7\nkgQEMF5SG9ABbAe21CdsMzMrp5EJZTrwfNH8qnRZ2TIRsRPYDEwlSS6vAWuB54B/iIgN5TYi6QJJ\nvZJ6161bl+8emJnZHo1MKCqzLDKWORXYBRwBzAb+UtJR5TYSEQsjoicieg4++ODhxGtmZlU0MqGs\nAo4smp8BrKlUJj291QVsAN4H3BoROyLiJeBuoKfmEZuZWUWNTChLgTmSZksaC8wFFpeUWQzMT6fP\nBe6IiCA5zfW7SowH3gQ8Uae4zcysjIYllLRN5BPAEuCXwE0RsVzSpZLelRb7FjBV0grgL4DCpcVX\nAROAx0gS0zUR8Uhdd8DMzAZQ8oP/wNDT0xO9vb2NDsPMrKlIeiAiBm1W8J3yZmaWCycUM6to0SKY\nNQtaWpLnRYsaHZGNZG2NDsDMRqZFi+CCC6CvL5lfuTKZB5g3r3Fx2cjlGoqZlbVgwd5kUtDXlyw3\nK8cJxczKeu65/Vtu5oRiZmV1d+/fcjMnFDMr67LLoLNz4LLOzmS5WTlOKGZW1rx5sHAhzJwJUvK8\ncKEb5K0yX+VlZhXNm+cEYtm5hmJmZrlwQjEzs1w4oZiZWS6cUMzMLBdOKGZmlgsnFDMzy4UTipmZ\n5cIJxczMclExoUh6p6RzyyyfJ+m02oZlZmbNploN5RLgp2WW3w5cWptwzMysWVVLKJ0Rsa50YUS8\nAIyvXUhmZtaMqiWUdkn79PUlaQzQUbuQzMysGVVLKN8H/kXSntpIOv31dJ2Zmdke1RLKXwEvAisl\nPSDpAeBZYF26zszMbI+K3ddHxE7gQkmXAEeni1dERH9dIjMzs6ZS7bLhzwKkCeTYiHi0kEwkfbFO\n8ZmZWZOodsprbtH0RSXrzqhBLGZm1sSqJRRVmC43b2ZmB7hqCSUqTJebNzOzA1y1hPJbkrZIegU4\nIZ0uzL8+j41LOkPSk5JWSLqwzPpxkm5M198naVbRuhMk3StpuaRHJbXnEZOZmQ1Ntau8Wmu5YUmt\nwFXAacAqYKmkxRHxeFGxjwAbI+JoSXOBy4Hz0hsurwc+EBEPS5oK7KhlvGZmVl3FhCLpoGovjIgN\nw9z2qSSXIT+Tbu8G4BygOKGcA1ycTt8CXClJwOnAIxHxcBrLy8OMxczMhqliQgEeIGkrEXA4sIa9\njfEBHDXMbU8Hni+aXwW8sVKZiNgpaTMwFXgdEJKWAAcDN0TEl8ttRNIFwAUA3d3dwwzZzMwqqXbK\na3ZhWtJDEXFSztsud6VYaWN/pTJtwFuBNwB9wO2SHoiI2/cpHLEQWAjQ09PjiwnMzGok6wBbtfgi\nXgUcWTQ/g6QWVLZM2m7SBWxIl/80ItZHRB/wI+DkGsRoZmYZNXLExqXAHEmzJY0luZFycUmZxcD8\ndPpc4I6ICGAJyZVnnWmi+V8MbHsxM7M6q9Yo/xdFs4eUzBMRXxnOhtM2kU+QJIdW4NsRsVzSpUBv\nRCwGvgVcJ2kFSc1kbvrajZK+QpKUAvhRRPxwOPGYmdnwKPnBX2aF9PlqL4yIS2oSUQ319PREb29v\no8MwM2sqaRt1z2DlqjXKN13CMDOzxqnW2/BNRdOXl6z7cS2DMjOz5lOtUX5O0fRpJesOrkEsZmbW\nxLJ2Drk/68zM7ABU7U75TkknkSSdjnRa6aOjHsGZmVnzqJZQ1gKFS4NfKJouzJuZme1R7Sqvt9cz\nEDMza26NvFPezMxGEScUMzPLhROKmZnlolqj/B6SpgMzi8tHxM9qFZSZmTWfQRNKepf8eSS9+e5K\nFwfghGJmZntkqaG8GzgmIrbVOhgzM2teWdpQngHG1DoQMzNrbllqKH3AMkm3A3tqKRHxqZpFZWZm\nTSdLQlnMviMpmpmZDTBoQomIa+sRiJmZNbdqQwDfFBHvlfQoZXoXjogTahqZmZk1lWo1lD9Nn3+/\nHoGYmVlzq9Y55Nr0eWX9wjEzs2blrlfMzCwXTihmZpaLTAlFUoekY2odjJmZNa9BE4qks4FlwK3p\n/ImSfF+KmZkNkKWGcjFwKrAJICKWAbNqF5KZmTWjLAllZ0RsrnkkZmbW1LJ0vfKYpPcBrZLmAJ8C\n7qltWGZm1myy1FA+CRxP0jHk94AtwJ/VMigzM2s+gyaUiOiLiAUR8YaI6Emnt+axcUlnSHpS0gpJ\nF5ZZP07Sjen6+yTNKlnfLelVSZ/OIx4zMxu6LCM2/oB9+/LaDPQC3xhqcpHUClwFnAasApZKWhwR\njxcV+wiwMSKOljQXKIweWXAF8F9D2b6ZmeUr6wBbrwL/kj62AC8Cr0vnh+pUYEVEPBMR24EbgHNK\nypwDFHo7vgX4PUkCkPTuNLblw4jBzMxykqVR/qSI+J2i+R9I+llE/I6k4XyZTweeL5pfBbyxUpmI\n2ClpMzBVUj/wOZLajU93mZmNAFlqKAdL6i7MpNPT0tntw9i2yiwrPbVWqcwlwBUR8eqgG5EukNQr\nqXfdunVDCNPMzLLIUkP5S+Dnkp4m+YKfDXxM0nj2no4ailXAkUXzM4A1FcqsktQGdAEbSGoy50r6\nMjAZ2C1pa0RcWbqRiFgILATo6enZZ1wXMzPLR5YRG3+U3n9yLElCeaKoIf6fhrHtpcAcSbOB1cBc\n4H0lZRYD84F7gXOBOyIigN8uFJB0MfBquWRiZmb1k6WGAnAKSXcrbcAJkoiI7w5nw2mbyCeAJUAr\n8O2IWC7pUqA3IhYD3wKuk7SCpGYydzjbNDOz2lHyg79KAek64DdIOojclS6OiPhUjWPLXU9PT/T2\n9jY6DDOzpiLpgYjoGaxclhpKD3BcDJZ5zMzsgJblKq/HgMNqHYiZmTW3LDWUacDjku4n6c8LgIh4\nV82iMjOzppMloVxc6yDMzKz5Zbls+Kf1CMTMzJpbliGA3yRpadqr73ZJuyRtqUdwZmbWPLI0yl8J\nnA88BXQAf5guMzMz2yPTjY0RsUJSa0TsAq6R5BEbzcxsgCwJpU/SWGBZ2nfWWmB8bcMyM7Nmk+WU\n1wfScp8AXiPprPE9tQzKzMyaT5arvFamk1tJuo03MzPbR5YhgN9Cci/KzOLyEXFU7cIyM7Nmk+WU\n17eArwBvBd5Q9DAzG3UWLYJZs6ClJXletKjRETWPLI3ymyPiv2oeiZlZgy1aBBdcAH19yfzKlck8\nwLx5jYurWVTsvl7Syenke0nGK/k+A/vyerDm0eXM3debWTWzZiVJpNTMmfDss/WOZuTIo/v6fyyZ\nL36zAH53KIGZmY1Uzz23f8ttoIoJJSLeXs9AzMwarbu7fA2lu7v+sTSjLH15fVHS5KL5KZL+trZh\nmZnV32WXQWfnwGWdnclyG1yWq7zOjIhNhZmI2AicVbuQzMwaY948WLgwaTORkueFC90gn1WWq7xa\nJY2LiG0AkjqAcbUNy8ysMebNcwIZqiwJ5XrgdknXkDTGfxi4tqZRmZlZ08nS9cqXJT0CvAMQ8IWI\nWFLzyMzMrKlk7b7+VuDWGsdiZmZNLEujvJmZ2aCcUMzMLBeZEoqkDknH1DoYMzNrXllubDwbWEba\nhiLpREmLax2YmY0c7oHXsshSQ7kYOBXYBBARy4BZtQupOfVv39XoEMxqotAD78qVELG3B14nFSuV\n5SqvnRGxWVLuG5d0BvBVkt6MvxkRXypZPw74LnAK8DJwXkQ8K+k04EvAWGA78JmIuCP3ADN6ZesO\nXn/xjzlo/FimT+7giMntTJ/cyfQpHUyf3MGM9Hly5xhqcRzNamnBgr3duRf09SXLfQOgFcuSUB6T\n9D6SO+bnAJ8C7hnuhiW1AlcBpwGrgKWSFkfE40XFPgJsjIijJc0FLgfOA9YDZ0fEGkm/CSwBpg83\npqEK4DPvPIbVm/pZvbGfp9e9xs9+tZ7+HQNrLePHtnLE5I49iWZgwunkkInjaGlxwrGRxT3wWlZZ\nEsongQUkY6H8K8mXdx6dQ54KrIiIZwAk3QCcAxQnlHNITrkB3AJcKUkR8VBRmeVAe3H3MPU2qX0M\nH3/70QOWRQQb+3awemM/qzf1sWpjP6s39bNmU/L88POb2Ni3Y8BrxrSKw7s60lpOknBmFCWewye3\nM66ttZ67ZuYeeC2zLAnlmIhYQJJU8jQdeL5ofhXwxkplImKnpM3AVJIaSsF7gIcalUwqkcRB48dy\n0PixvH5GV9kyr23byZpN/axKazari57vXrGeF1/ZSvH4ZxIcPGHcgBrO3oSTnGKbMC7TvapmmV12\n2cBRDME98Fp5Wb59viLpcOBm4IaIWJ7Ttsud2ykdPrJqGUnHk5wGO73iRqQLgAsAukfYT6rx49qY\nc+hE5hw6sez67Tt388Lmraza1LdPwnl09WZ+vPxFtu/aPeA1k9rbmD6lc8+ptBlT0tpOmnimjh/r\ndhzbL4V2kgULktNc3d1JMnH7iZXK0pfX2yUdRjIU8EJJk4AbI2K4p71WAUcWzc8A1lQos0pSG9AF\nbACQNAP4d+D/RsTTVeJfCCyEZAjgYcZcV2PbWuie2kn31M6y63fvDta/uq1sDef5DX384pmXeXXb\nzgGvaR/TsifBzBjQlpPUcA6dOI62Vt/vagO5B17LImtfXi8AX5N0J/BZ4G8YfjvKUmCOpNnAamAu\n8L6SMouB+cC9wLnAHRER6YBfPwQuioi7hxlH02ppEYdMaueQSe2c3D1ln/URwZb+nWVrOKs39fP4\nmi28/Nr2Aa9pbRGHTWofcNFA6XP7GLfjWO0tWuRaUbMZNKFI+h8kV1adS3Lp7g3AXw53w2mbyCdI\nGvlbgW9HxHJJlwK9EbEY+BZwnaQVJDWTuenLPwEcDfy1pL9Ol50eES8NN67RRBJdnWPo6uzi+CPK\nt+P0b9/Fms0DazirNvaxelM/9/96A2s397O7pF43bcLYgYlmcsee02zTp3TQ1TGmDntno1nh3pdC\nu03h3hdwUhnJFFH9LJCkXwDfA26OiNJTUk2lp6cnent7Gx1GU9m5azcvbNlatoZTmN62c2A7zsRx\nbUwvabspvohg2oTRe3m0f1XnY9as8leWzZwJzz5b72hM0gMR0TNYuSxtKG/KJyRrRm2tLcyY0smM\nKeXbcSKCl1/bXraGs2pjP0uf3cArWwe244xta+GIrvaiGs7Am0AP62pnTBO24/hXdX5870tzqlhD\nkXRTRLxX0qMMvPpKQETECfUIME+uoTTGK1t3DKzdbBx4qfS6VwZe8d0iOLSoHae4plO4TLpz7Mi7\nPNq/qvPjYzmy5FFD+dP0+ffzCckOVBPbx3DsYWM49rBJZddv3bGLtZuT02qF+3JWbUwuJHhg5UZ+\n+MhadpY05EzpHFOxhtOobm78qzo/vvelOVVMKBGxNp38WER8rnidpMuBz+37KrP91z6mldnTxjN7\n2viy63ftDl56ZW87zqqimk6lbm46x7aWv1ItnT5kYjutObfj+I7y/Pjel+aUpVH+wYg4uWTZIz7l\nZSPFwG5u9rbhFF9AsKlKNzeF02rD7eamtA0Fkl/VCxf6i9Ca27BPeUn6KPAx4ChJjxStmggcsPd+\n2Mgz3G5ufv5U9m5uijv3nNg+8PJo/6q2A121RvkuYArwd8CFRateiYgNdYgtd66hWCXVurlZvamf\ntZu27tPNTVfHmLKn09zNjY02w66hRMRmYDNwfvqGhwDtwARJEyLCTY02amTp5mbdq9tYtXFvj9GF\nhPPcy33cs2I9r5UMslatm5sjJrdz2KR2d3Njo0qWO+XPBr4CHAG8BMwEfgkcX9vQzEaOlhZx6KR2\nDp3Uzikzc+7mpmSYguLTau7mxppJlov5/xZ4E/CTiDhJ0ttJay1mlsjazU1xLwPFNZ37fr2Btcvc\nzY01tywJZUdEvCypRVJLRNyZXjZsZvuhY2wrRx8ygaMPmVB2fbVubp5Y+wq3//Klit3clLbfFK5a\nG83d3NjIkyWhbJI0AfgZsEjSS8DOQV5jZvspSzc361/dvk8bTuG+nKXPbmBLaTc3rS0cMXn0dXNj\nI1OW+1DGA1tJulyZRzImyaKIeLn24eXLV3nZaLdl644k4Qyxm5vi58LgbCOxmxurrzw7h3ytaPba\nYUVlZjU1qX0MkzJ0c7Nmn2TTx4PPNU83NzYyVbux8RXKdArJ3s4hy//FmtmIVY9ubo4oGn56+uRO\nDpnodpwDRbX7UMoPdG5mo1ZrS9IlzeFdHZQ7vzGwm5u+AQln9aZ+lj2/qWw3N4d1tQ+o4Qy3mxsb\nmTKdHJX0VmBORFwjaRowMSJ+XdvQzGykydrNTaXB2O5ekb2bm+KhC0q7uWk0D6RWXpYbGz8P9ADH\nANcAY4HrgbfUNjQza0bjx7XxukMn8rpDy5/kKNfNTeHKtcdWb+bHy1/cp5ubSe1te+6/mVHmMul6\ndnPjgdQqy3KV1zLgJODBiDgpXebehs2sJnbvDta/uq1sR56F51e3Dbw8up7d3ByIg3/ldpUXsD0i\nQlKkb1y+Ne8A5uqvWX5aWsQhk9o5ZFI7J3dn7+Zmzebk+ba1W1j/aoVubsqNkbOf3dx4ILXKsiSU\nmyR9A5gs6Y+ADwPfrG1YzcPVX7P6ytLNzdYdu0ouGOhjzabk6rX7f72BF7ZsZVfJ5dGVurl5+J4O\nrry8g+eebqO7Wxx0ELxc5i48D6SW4ZQXgKTTgNNJLhleEhG31TqwWqjFKa8Dsfpr1uyqdXOzemM/\nz73cz84Y2I6ze1sbO7d0sPuVDnZt6WD7xg52bkkeY7d3cPUV4/jA+0fn5dFZT3llSiglb9wKzI2I\nRUMNrlFqkVBaWqDcIZRg9+59l5vZyLZoEXzoQ8GuMdtp6+qnbVI/bV19tE7qp23SVtom9TNmch8a\nV9LNTVsLR3SNzm5u8hixcRLwcWA6sBi4LZ3/DLAMaLqEUgseR9zA7WijyYIFsGOHYMc4tveNY/va\nyfuUkWBz3459LhgoXEhw55PrMnVzc0TJZdLN3s1NteivAzYC9wJ/SJJIxgLnRMSyOsTWFC67rPw4\n4pdd1riYrL7cjja6ZGlc7+6Gie1jOLZKNzfbdu5i7aatZfpUG73d3FQbAvjRiHh9Ot0KrAe6I+KV\nOsaXq1pdNuxfpwc2t6ONLpU+zwIJ/uRP4OqrBy7f3++Bat3cFJ4rdXNTPAjbjCkdvPk3pnLIxPah\n7/Qght2GIunBiDi50nwz8n0oVgtuRxv59ufLPmlDgR07yq+H5CzEwoV736O0llquzP4arJub1Zv6\n93Rzc91HTuW35xw8tA1lkEdC2QUUehoW0AH00cSdQzqhWC24hjKylfuyB5g6Fb761eQLvzThnHUW\n3HRT+cuDC4o/30p/A62tcO21tTtj8dq2nazZ1M8RkzsYP6527S9ZEwoR0bAHcAbwJLACuLDM+nHA\njen6+4BZResuSpc/Cbwzy/ZOOeWUMMvb9ddHdHZGJPWU5NHZmSy3xps5c+BnU/o5ffSj1T8/qfxr\npb3bqFRmtPwtAL2R4Tu2Ydewpe0yVwFnAscB50s6rqTYR4CNEXE0cAVwefra44C5wPEkSenq9P3M\n6m7evOTUxsyZyWmumTOHd6pjtFu0KPlF39KSPC+q8fWi1RrZ+/qSz6q09tLXB/PnJ7FVumIzAqZN\nq16m8F7vf39t9rXex3JQWbJOLR7Am0lukiyucVxUUmYJ8OZ0uo3kwgCVli0uV+3hGopZYzWiNlet\nhjLYo1CDaWmpXGbs2PK1nFrXVup5LBnpNRSS+1ueL5pflS4rWyYidgKbgakZX2tmI8yCBeVrAwsW\n1G6bl12WNJBX0lrl3EZfH/zoR9Xff/v2pMzChdXfq/B+ee1rI47lYBqZUMpdTF16hUClMllem7yB\ndIGkXkm969at288QzSxPjehYsXBKcurUfdd1diYN9tUSznPPDX613nPPJdu59trq71Uom4eR2Ell\nIxPKKuDIovkZwJpKZSS1AV3AhoyvBSAiFkZET0T0HHxw7S6rM7PBVWprqHXPEvPmwfr1cP31+7Z1\nXX119dpFd/fgNY9C/MXtaYOVHa5GHctqGplQlgJzJM2WNJakkX1xSZnFwPx0+lzgjvR83mJgrqRx\nkmYDc4D76xS3mQ1RudNP9exZYt685FLf3buT58KFE5VqF4XYCj0fVFIcf2Eb119f231t9LEsK0tD\nS60ewFnAr4CngQXpskuBd6XT7cDNJJcH3w8cVfTaBenrngTOzLI9N8qbNd711ycN5VLyPJIuqa0W\n2/jx5Rvap04d2vvVOt48kbFRfr97G25mvrHRbPSpV9dHtbgbvllkvbGxOftSNjNj75f8ypVJfaHQ\nMWct7sfw/UaDcw3FzJqWu72pD9dQzGzUG4mXzh7InFDMrGmNxEtnD2ROKGbWtEbkpbMHMCcUM2ta\nbigfWZp7AGMzO+DNm+cEMlK4hmJmZrlwQjEzs1w4oZiZWS6cUMzMLBdOKGZmlgsnFDMzy4UTipmZ\n5cIJxczMcuGEYmZmuXBCMTOzXDihmJlZLpxQzMwsF04oZmaWCycUMzPLhROKmZnlwgnFzMxy4YRi\nZma5cEIxM7NcOKGYmVkunFDMzCwXTihmZpaLhiQUSQdJuk3SU+nzlArl5qdlnpI0P13WKemHkp6Q\ntFzSl+obvZmZldOoGsqFwO0RMQe4PZ0fQNJBwOeBNwKnAp8vSjz/EBHHAicBb5F0Zn3CNjOzShqV\nUM4Brk2nrwXeXabMO4HbImJDRGwEbgPOiIi+iLgTICK2Aw8CM+oQs5mZVdGohHJoRKwFSJ8PKVNm\nOvB80fyqdNkekiYDZ5PUcszMrIHaavXGkn4CHFZm1YKsb1FmWRS9fxvwPeBrEfFMlTguAC4A6O7u\nzrhpMzPbXzVLKBHxjkrrJL0o6fCIWCvpcOClMsVWAW8rmp8B3FU0vxB4KiL+aZA4FqZl6enpiWpl\nzcxs6Bp1ymsxMD+dng/8vzJllgCnS5qSNsafni5D0t8CXcCf1SFWMzPLoFEJ5UvAaZKeAk5L55HU\nI+mbABGxAfgCsDR9XBoRGyTNIDltdhzwoKRlkv6wETthZmZ7KeLAOQvU09MTvb29jQ7DzKypSHog\nInoGK+c75c3MLBcHVA1F0jpgZQM2PQ1Y34Dt5qXZ4wfvw0jhfRgZ9mcf1gNExBmDFTygEkqjSOrN\nUl0cqZo9fvA+jBTeh5GhVvvgU15mZpYLJxQzM8uFE0p9LGx0AMPU7PGD92Gk8D6MDDXZB7ehmJlZ\nLlxDMTOzXDih5CDLgGGSTpR0bzoo2COSzita9x1Jv07v+l8m6cT67kEu+zBb0n3p62+UNLa+e7Bf\nA7fdKmmTpP8sWd4Un0NartI+NNPnsM8AeunyuyQ9WfQ5lOuNvFaxn5Fue4WkcuM0jUuP64r0OM8q\nWndRuvxJSe+sV8wl8Q0pfkmzJPUXHfOvDymAiPBjmA/gy8CF6fSFwOVlyrwOmJNOHwGsBSan898B\nzm3yfbgJmJtOfx346Ejch3Td75EMe/CfJcub4nMYZB+a4nMADgKeSZ+npNNT0nV3AT0NiLsVeBo4\nChgLPAwcV1LmY8DX0+m5wI3p9HFp+XHA7PR9Wpso/lnAY8OOod4f2mh8AE8Ch6fThwNPZnjNw0Vf\nziPhi2zI+0Ay1MB6oC1d/mZgyUjeB5KerEdiQhnyPjTT5wCcD3yjaP4bwPnpdKMSyoDjBVwEXFRS\nZgnw5nS6LT3eKi1bXK5J4s8lofiUVz6yDBi2h6RTSX5BPF20+LL0NNIVksbVLtSKhrMPU4FNEbEz\nXb3PYGh1sl/7UEFTfQ4lmulzGGwAvWvSUy9/Lanc2Ei1MOigfsVl0uO8meS4Z3ltrQ0nfoDZkh6S\n9FNJvz2UAGo2Hspoo+EPGFZ4n8OB64D5EbE7XXwR8ALJF/RC4HPApUOPtuK2a7IPFf7ha3L5YF77\nUEFTfQ7l3rrMspH6OVSLdV5ErJY0Efg34APAd/c/yv2W5fhVKlO3Y1/FcOJfC3RHxMuSTgH+Q9Lx\nEbFlfwJwQskohj9gGJImAT8E/ioiflH03mvTyW2SrgE+nWPoe9RwH9YDkyW1pb96ZgBrcg4fyGcf\nqrx303wOFTTT51BxAL2IWJ0+vyLpX4FTqU9CWQUcWRJT6fErlFmlZNTYLmBDxtfW2pDjj+Qc2DaA\niHhA0tMkbab71T27T3nlY9ABw9Krbf4d+G5E3Fyy7vD0WcC7gcdqGm15Q96H9I/xTuDcaq+vgywD\nt1XULJ9DJU32OZQdQE9Sm6RpAJLGAL9P/T6HpcCc9Eq5sSSN1otLyhTv27nAHelxXwzMTa+imk3S\ntnh/neIuGHL8kg6W1Aog6SiS+CsOrV5RPRuNRuuD5Bzk7cBT6fNB6fIe4Jvp9PuBHcCyoseJ6bo7\ngEdJ/nGuByY04T4cRfIPtALO7drBAAAGfklEQVS4GRg3Evchnf9vYB3QT/KL7Z3N9DkMsg/N9Dl8\nOI1zBfChdNl44AHgEWA58FXqeLUUcBbwK5K2wQXpskuBd6XT7elxXZEe56OKXrsgfd2TwJn1Pu7D\niR94T3q8HwYeBM4eyvZ9p7yZmeXCp7zMzCwXTihmZpYLJxQzM8uFE4qZmeXCCcXMzHLhhGL7kLQr\n7fbiMUk/kDQ5w2teHWT9ZEkfK5o/QtItecRbYXt3SWr4uN+13s90GydKOivn92zI8cu6L6XlJL2r\nXO+6Vl9OKFZOf0ScGBG/SXIX8MdzeM/JJD2dAhARayLi3CrlR4U67eeJJPcfZKbESPz/z7ovA8pF\nxOKI+FLNorJMRuIflI0s91LUwZykz0hamnageElpYUkTJN0u6UFJj0o6J131JeA30prP36fjLzyW\nvuY+SccXvcddkk6RNF7St9PtPVT0XqXb/Gy6rYclFX+p/IGk+yX9qtDZXbrd/07je1DS/0yXvy3d\n7i2SnpC0qNBHmaSz0mU/l/Q1pWOQZImvZD8/KOn7SsYyeUrSlyvsz7OSLik6hsemy0+VdE+6rXsk\nHZPeEX0pcF56bM+TdLGkTxe932NpHLMk/VLS1SQ3rx0p6Z8l9SoZ42afz7MkrjMl3VQ0/zZJP0in\nT1cyVs6Dkm6WNGF/j12Ffcm6zx+UdGX63jPTv8FH0ufudPl30hjukfSMpFH/g6buGnE3px8j+wG8\nmj63ktxVe0Y6fzpJp4ki+THyn8DvlLymDZiUTk8juSN3n+6xi+eBPwcuSacPB36VTn8ReH86PZnk\nDuDxJbGeCdwDdKbzhbuy7wL+MZ0+C/hJOt0JtKfTc4DedPptJD2vzkj37V7grSR3Fj8PzE7LfY+0\ny/iM8RXv5wdJurPoSt93JXBkmeP/LPDJdPpj7O2pYBJ7u6Z/B/BvRe97ZdHrLwY+XTT/WBrHLGA3\n8KaidYXj1ZoesxOKjl9PSVxtwHOFfQT+maT3hGnAz4qWfw74m6EcuzL7knWf98wDPyDpuBSSu/H/\nI53+DsnfcwvJ+CUrGv2/NtoerqFYOR2SlgEvkwyAdFu6/PT08RDJL9xjSb6Uiwn4oqRHgJ+Q1G4O\nHWR7NwF/kE6/l+SfvrC9C9NY7iL5guouee07gGsiog8gIjYUrft++vwAyZcpwBjgXyQ9mm7nuKLy\n90fEqkh6gV6WvuZY4JmI+HVa5ntF5bPEV+r2iNgcEVuBx4GZFcqVi70LuDmt8VwBHF/mdYNZGUUd\nkwLvlfQgyWd6PAOPxwCRdDh5K3C2ko4F/zdJP11vSl93d3os5qf7lcexG8o+vxn413T6OpIfBgX/\nERG7I+JxBv+7tP3k3oatnP6IOFFSF0kt5OPA10iSxd9FxDeqvHYecDBwSkTskPQsyZdFRZF0Vf6y\npBOA84A/TlcJeE9EPFnl5aJyN+Hb0udd7P1b/3PgReC3SH6pbi1Tvvg11cbiyBJfpZhK46pUrrjM\nF4A7I+L/KBm69a4Kr93JwNPZxcf/tcKEkk4MPw28ISI2SvoOg3xWwI0kfw8bgKWR9Ags4LaIOL+4\noKSTqrxP2WMn6Y0l5bLuczXFfx/Fx79e46wcMFxDsYoiYjPwKeDTSnp+XQJ8uOj8+HTtO953F/BS\nmkzezt5f4K8AE6ts7gbgs0BXRDyaLlsCfLKoLaPcF9SP05g60zIHDbJbXcDatBbyAZJTPdU8ARyl\nvWOHn1e0Lkt8eeoCVqfTHyxaXnpsnwVOTmM6mWRI2nImkSSYzZIOJTl9OJi70vf+I5LkAvAL4C2S\njk632SnpdQzt2JXuS9Z9LnYPSU+7kPzA+XmG/bIcOKFYVRHxEEkPpHMj4sckpxLuTU8Z3cK+/9SL\ngB5JvST/zE+k7/MyySmRxyT9fZlN3ULyJXBT0bIvkJyieiQ95fGFMvHdStIld296+mSwMUyuBuZL\n+gXJeA+vVSscEf0k7Ri3Svo5Se1mc9b4cvZl4O8k3c3ARHgncFyhgZpkUKqD0uPxUZL2iX1ExMMk\np7qWA98G7h4sgIjYRVJrPTN9JiLWkXzZfy891fkL4NghHrvSfcm6z8U+BXwojeUDwJ8Otl+WD/c2\nbDYISRMi4tX01/RVwFMRcUWj42oGPnYHFtdQzAb3R+mv/eUkp2CqtSHZQD52BxDXUMzMLBeuoZiZ\nWS6cUMzMLBdOKGZmlgsnFDMzy4UTipmZ5cIJxczMcvH/AVMNZqDcor7CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated p: 0.11594 (0.0, 0.11794043933579963)\n"
     ]
    }
   ],
   "source": [
    "coeffs_edit = do_regress(ecoli_edit, nature_edit1)\n",
    "vals = permute_regress(ecoli_edit, nature_edit1) # permutation distribution of estimated slope\n",
    "p_CI(vals, abs(coeffs_edit[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fourth-root transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.254\n",
      "Model:                            OLS   Adj. R-squared:                  0.226\n",
      "Method:                 Least Squares   F-statistic:                     8.873\n",
      "Date:                Sun, 16 Sep 2018   Prob (F-statistic):            0.00620\n",
      "Time:                        20:38:39   Log-Likelihood:                 9.3891\n",
      "No. Observations:                  28   AIC:                            -14.78\n",
      "Df Residuals:                      26   BIC:                            -12.11\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.1275      0.038      3.348      0.002       0.049       0.206\n",
      "x1            -1.5264      0.512     -2.979      0.006      -2.580      -0.473\n",
      "==============================================================================\n",
      "Omnibus:                        1.056   Durbin-Watson:                   1.957\n",
      "Prob(Omnibus):                  0.590   Jarque-Bera (JB):                0.995\n",
      "Skew:                          -0.407   Prob(JB):                        0.608\n",
      "Kurtosis:                       2.563   Cond. No.                         15.1\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4HHd97/H3V5JlW77bsqSNHVk2\ntuPYWuMQkxCSQhLHjiVDgKdAoC5NyoGcklL6tIe2pKE9NNSF0p62cCgXl8M15pLQAgFLcUISk5K7\nnYtXduJLnPiSrHy/3yV9zx8zsteKLiNrVrvSfl7Po2dnZmd3vjNr72dnfjO/MXdHREQkiqJcFyAi\nIgOHQkNERCJTaIiISGQKDRERiUyhISIikSk0REQkMoWGiIhEptAQEZHIFBoiIhJZSa4LiFt5ebnX\n1NTkugwRkQFl7dq1e919Yk/zDbrQqKmpYc2aNbkuQ0RkQDGzbVHm0+EpERGJTKEhIiKRKTRERCQy\nhYaIiESm0BARkcgUGiIiEplCQ0REIlNoiIhIZAqNArFiBdTUQFFR8LhiRa4rEpGBaNBdES5vtGIF\n3HYbHD8ejG/bFowDLF2au7pEZODRnkYBuPPOc4HR7vjxYLqISG8oNArA9u29my4i0hWFRgGoru7d\ndBGRrig0CsCyZVBWdv60srJguohIbyg0CsDSpbB8OUyZAmbB4/LlagQXkd7T2VMFYulShYSI9J32\nNEREJDKFhoiIRKbQEBGRyBQaIiISWU5Dw8wWm9lGM9tiZp/pYp4PmtkGM1tvZj/s7xpFROScnJ09\nZWbFwL8DC4GdwDNmdp+7b8iYZwZwB3C1ux8ws4rcVCsiIpDbPY0rgC3uvtXdTwM/Bt7TYZ6PA//u\n7gcA3H13P9coIiIZchkak4AdGeM7w2mZZgIzzewxM3vSzBb3W3UiIvIGuby4zzqZ5h3GS4AZwLXA\nZOC/zazW3Q+e90ZmtwG3AVSrQyURkazJ5Z7GTuDijPHJwOudzPMLdz/j7q8AGwlC5Dzuvtzd57v7\n/IkTJ2atYBGRQpfL0HgGmGFmU82sFPgQcF+HeX4OXAdgZuUEh6u29muVIiJyVs5Cw91bgE8Cq4AX\ngXvcfb2Z3WVmN4WzrQL2mdkG4BHgL9x9X24qFhERc+/YjDCwzZ8/39esWZPrMkREBhQzW+vu83ua\nT1eEi4hIZAoNERGJTKEhIiKRKTRERCQyhYaIiESm0BARkcgUGiIiEplCQ0REIlNoiBS4FSugpgaK\nioLHFStyXZHks1z2cisiObZiBdx2Gxw/Hoxv2xaMAyxdmru6JH9pT6MXzrS2cexUS67LEInNnXee\nC4x2x48H00U6oz2NXnhsy17+6O61XDuzgrpkFQsurWTkUG1CGbi2b+/ddBF94/XCpLHDuXn+xTQ2\nNXP/+mZKS4p458yJ1IcBMnrYkFyXKNIr1dXBIanOpot0Rr3cXoC2NmfNtgM0NqVpTDXTfPgkpcVF\n/M6McuqSCRbOrmTMcAWI5L+ObRoAZWWwfLnaNApN1F5uFRp91NbmPLfjII2pNI1Nzbx28ARDio2r\np5dTn0ywaHYlY8tK+60ekd5asSJow9i+PdjDWLZMgVGIFBo54O48v+MgjU3NNKTS7DxwgpIi4+3T\ny6mvrWLRnCrGj1CAiEj+UWjkmLuTeu0QDakgQLbvP05xkXHVtAnUJau4cU4V5SOH5rpMERFAoZHr\nMs7j7qx//TCNTWkaUs28svcYRQZXTp1AfbKKG2urqBg1LNdlikgB63NomNmNwCh3/2mH6UuB3e7+\nYCyVxiwfQyOTu/NS8xEaU2lWptK8vOcYZvDWmvEsSSZYXFtF5WgFiIj0rzhC40ng3e6+p8P0KuBn\n7n5VLJXGLN9DI5O7s2nXURpSaRqb0mzadRQzmD9lHHW1CeqSVSTGDM91mSJSAOIIjXXuPre3z+Xa\nQAqNjrbsPnK2DeSl5iMAvKV6LPXJBHXJBJPGKkBEJDviCI1NwGx3b+kwfQiwwd1nxFJpzAZyaGTa\nuucojU3NrFyXZkP6MABvvngsS5JV1NUmuHh8WY4rFJHBJI7Q+CJQCXzS3Y+F00YAXwH2uvtfxVhv\nbAZLaGR6de8xGpuaaWxKs27nIQDmTh5DXW2C+mQVUyaMyHGFIjLQxREaJcDfAx8D2jsaqAb+H/A3\n7n4mplpjNRhDI9OO/cdpbEqzMtXMCzsOAjDnotHUJxPUJxNMLVeAiEjvxXbKrZkNB6aHo1vc/UQM\n9WXNYA+NTDsPHOf+pmZWptI8tz0IkFlVo1gStoFMrxiZ4wpFZKCIY0/jL939S+HwB9z93ozn/sHd\n/zq2amNUSKGR6fWDJ4KOFJvSrNl2AHeYWTmSutoES+YmmFk5KtclikgeiyM0nnX3t3Qc7mw8nxRq\naGTadfgkjangQsJntu3HHaZXjKS+toq6ZIJZVaMws1yXKSJ5JGpodNc1unUx3Nm45JHK0cO49eqp\n3Hr1VHYfPsmq9c00pJr56iNb+MrDW5hWPoK68CysOReNVoCISGTdhYZ3MdzZuOSpitHD+MhVNXzk\nqhr2HDnFAxuC60C+vvpl/v2Rl5kyoezsWVjJSWMUICLSre4OT7UCxwj2KoYD7T3uGzDM3fPyhhE6\nPBXNvqOneGDDLhpSaZ54eR8tbc7kccPPnoX15skKEJFCog4LJbIDx07z4ItBgDy2ZS9nWp1JY4dT\nF7aBXHbxWIqKFCAig1kcDeHju3uhu++/wNqySqHRN4eOnzkbIL/dvJfTrW0kxgxjcW0V9ckEl1eP\nU4DIoKCbT50vjtB4haDtwoAE8DrnGsDd3afFVGusFBrxOXzyDA+9uIuV65p5dPMeTre0UTFqKHVh\ngMyvGU+xAkQGIN3m9o1iPTxlZs+5+2WxVHb++y4GvgwUA99y9y92Md/7gXuBt7p7t4mg0MiOIyfP\n8PBLu2lIpVm9cQ+nWtooHzmUxbWV1CcTXFEznpLiolyXKRJJTQ1s2/bG6VOmwKuv9nc1+SHu0Ij9\nugwzKwY2AQuBncAzwIfdfUOH+UYBK4FSgn6wFBo5duxUC49sDALk4Zd2c/JMGxNGlLJoThVLkgne\nNk0BIvmtqAg6++ozg7a2/q8nH8RxnUa2XUHQLclWADP7MfAeYEOH+T4PfAn4dP+WJ10ZMbSEd829\niHfNvYjjp1tYvXEPDak0v3j+NX709HbGlQ3hxjnBIayr3jSBIQoQyTPV1Z3vaVRX938tA02XoWFm\nf54xWtFhHHf/lz4uexKwI2N8J3BlhxouAy5291+ZmUIjD5WVlpw9TffE6VZ+s2k3jU3N/PKF1/nx\nMzsYM3wIi2ZXUj83wdVvKqe0RAEiubdsWedtGsuW5a6mgaK7PY3Mzor+o8N4HDprQT27w2hmRcC/\nArf2+EZmtwG3AVTrp0LODC8tZnFtgsW1CU6eaeXRTXu4v6mZ+5uauXftTkYNK2Hh7EqWJBNcM6Oc\noSXFuS5ZOiiUM4ra16kQ1jVuObtOw8yuAj7n7jeG43cAuPsXwvExwMvA0fAlVcB+4Kbu2jXUppF/\nTrW08tvNe2lINfPghmYOn2xh1NASFlxaQX0ywTtmTmTYEAVIrumMosIWxym397j7B8Phf8y86ZKZ\nPeDui/pYYAlBQ/gC4DWChvDfc/f1Xcy/Gvi0GsIHttMtbTz28l4aU2ke2LCLg8fPMKK0mAWXVlKf\nrOLaSyoUIDmiM4oKWxwN4Zm3c10IZN6pb+KFFtbO3VvM7JPAKoJTbr/t7uvN7C5gjbvf19dlSP4p\nLSniuksquO6SCpa1tvHEy/tobEqzav0u7nvhdcpKi7luVgVLkgmuvWQiZaW5PFejsGzf3rvpUpjU\nNbrkhZbWNp56ZT8rU2keWN/M3qOnGTYkCJj6ZILrZ1UwYqgCJJu0p1HY4tjTKAvPXioChofDxrkO\nDEViU1JcxNXTy7l6ejmff08tT7+yn4ZUOrw3ejNDS4q49pKJZwNk1LC87C9zQMvVGUWF0vg+WHS3\np/FIdy909+uyUlEfaU9jcGltc9ZuOxAGSJpdh09RWlLEO2ZMpD5ZxQ2zKxndzwEymL/k+nvd1Pie\nP9TLrQw6bW3Os9sP0JBqprEpTfrQSYYUG78zI9gDWXhpJWPKshsg+pKLlw6J5Q+FhgxqbW3O8zsP\n0rAuOIT12sETlBQZV08vZ0kywcLZlYwbURr7cvUlFy9155E/FBpSMNydF3YeCu6L3pRmx/4TFBcZ\nb3/TBOqTCRbNrmTCyKGxLEtfcvFSCOcPhYYUJHdn/euHWZlK05BKs23fcYqLjLdNG09dbYLFtVWU\n9yFA9CUXLx3uyx9x93I7CZhCxtlW7v5onyrMEoWGtHN3NqQP05gK7ou+de8xigyunDqB+mQVN9ZW\nUTFqWK/eU19y8RvMJxYMJLGFhpn9I3AzQe+zreFkd/eb+lxlFig0pDPuzsZdR2hYl2ZlKs3Le45h\nBm+tGU99bRWLaxNUjYkWIPqSk8EoztDYCMx191NxFZdNCg2JYtOuI8FpvKlmNu46AsDlU8ZRn0xQ\nV1vFRWN1KZIUljhDoxH4gLsf7XbGPKHQkN7asvto2IjezIvpwwDMu3gsS5JBG8jF48tyXKFI9sUZ\nGv8JvBl4CDi7t+Hun+prkdmg0JC+eGXvsbMXEja9FgTImyePoS6ZoL42QfUEBYgMTnGGxi2dTXf3\n711gbVml0JC4bNt3jMamoBF93c5DANROGh3cdKo2QU35iBxXKBIfnXIrEqMd+4/T2JSmIdXM8zsO\nAjA7MZr6ZHBb22kTR+a4QpG+ie1+GmaWIuOOeu3cfW7fy4yfQkOy7bWDJ2gMO1Ncu+0AALOqRoW3\nva1iekXcN7kUyb44QiPh7mkzm9LZ8+7eySVOuafQkP6UPnSC+8NDWGu2HcAdZlSMpC6ZYEkywczK\nkZh1dmdjkfyiw1Mi/WzX4ZOsWt/MynVpnn51P+4wbeII6msT1CcTXJoYpQAZQArtehyFhkgO7T5y\nklXrd9GwLs1Tr+yjzaFmQll4CCvBnItGK0DyWCFe+a/QEMkT+46eCgIkleaJrftobXOqx5dRl6xi\nSTJBctIYBUieKcQ+xuLue2o4UO3uG+MoLpsUGpLP9h87zYMbmmlINfPYlr20tDmTxw0/eyX6vIvH\nKkDyQCH2Zhw1NIoivNG7geeB+8PxeWZ2X99LHFhWrAh+fRQVBY8rVuS6IhmIxo8o5ea3VvO9j17B\nms/ewD+9fy4zK0fxncde4X1fe5yrv/gwd/1yA2te3U9b2+A6CjCQVFf3bnohiXJx31rgemC1u18W\nTltXSKfcFuLxTelfh06c4dcbdtHYlObRTXs53dpG5eih1IWN6POnjKOoSHsg2ZTZ8D1+PBw+DGfO\nnHt+sP+fj/OK8Kfc/Uoze65QQ6MQj29K7hw5eYaHXtxNQyrN6k17ON3SRsWooSyuDS4kfGvNeIoV\nILHq7IdhaSmMGgX79+vsqUwlPc0ANJnZ7wHFZjYD+BTweF8LHEi2b+/ddJG+GDVsCO+9bBLvvWwS\nR0+18PBLu2lMpblnzQ6+/8Q2ykeWcuOcIECunDqekuIejzJLF26/Pdh7aG1943OnT8PIkbB3b//X\nlc+i7GmUAXcCiwADVgGfd/eT2S+v97SnIYPV8dMtPPLSHhpSaR5+aTcnzrQyYUQpi+ZUUZ+s4m3T\nJjBEARLZ7bfD17/e/TyDueG7I51yGyO1aUi+OXG6ld9s2k1DqpmHXtzFsdOtjCsbwqLZVdQlq7h6\nerkCpAclJZ3vYWQqpB+GsR2eMrNf8sa+pw4Ba4Bv5useR5zag6GQrg6VzuXLVcLDS4tZXJtgcW2C\nk2daeXRTsAeyMpXmJ2t2MGb4EBbOrqQ+WcU10ydSWqIA6ainwCgrCz5fOV+Uw1NfBiYCPwon3Qw0\nA8OB0e7+kaxW2Eu6TkOyZSDscZ4808pvN++lIZXmwQ27OHKqhVHDSlh4aSX1yQTXzChn2JDiXJeZ\nF7rb05gypfB+GMZ59tSj7v6OzqaZ2Xp3n9PHWmOl0JBsGWhtW6daWnl8yz4aUmke2LCLQyfOMHJo\nCQsuraA+meCdMycWdIB01abxiU/A177W//XkWpxnT000s2p33x6+cTVQHj53ug81igwoA+0suqEl\nxVw3q4LrZlWwrKWNJ7buo2FdmlUbmvnF868zorSY6y+tpL62imsvqWB4aWEFSHswtJ89VVwc7EkW\nYmD0RpQ9jXrgG8DLBGdPTQVuB1YDH3f3f8tyjb2iPQ3JloG2p9GVM61tPLV1PytTaVatb2b/sdMM\nH1LM9bMqqEtWcf2sCspKo/yelMEk7r6nhgKzCELjpXxu/FZoSLYMhDaN3mppbePpV84FyN6jpxk2\npIhrZwYBsuDSSkYOVYBkypeTIeIWd2i8Hagh43CWu3+/LwVmi0JDsmmwfmEAtLY5z7y6n4bwroR7\njpyitKSId86cyJJkggWXVjBq2JBcl5lTg/GHQ7s4G8J/ALyJoNPC9nMN3N0/1ecqs0ChIdJ3rW3O\ns9sPsHJdmvubmmk+fJLS4iLeMbOc+mSCBZdWMmb44AuQnn4UDJZDlJ2JMzReBGZ7Fq4CNLPFwJeB\nYuBb7v7FDs//OfAxoAXYA3y0p9vMKjRE4tXW5jy342CwB5JK8/qhkwwpNq6ZHgTIwtmVjC0rzXWZ\nfd4LjLIXMZi7TI8zNO4FPuXu6biKC9+3GNgELAR2As8AH3b3DRnzXAc85e7HzewTwLXufnN376vQ\nEMked+f5MEAaUs28dvAEJUXG26eXsyRZxcLZVYwf0f8BEsdhoyh7EdrTiBYajwDzgKeBU+3T3f2m\nPhZ4FfA5d78xHL8jfN8vdDH/ZcBX3f3q7t5XoSHSP9yd1GuHWJlK05hqZvv+4xQXGVdNm0B9MsGi\nOZWUjxzaL7XE8WUeZS9CbRrRrtP4XN/L6dQkYEfG+E7gym7m/x9AY5ZqEZFeMjPmTh7L3Mlj+czi\nWax//XC4B5Lmr3+W4rM/T3Hl1AnUz01w45xKKkYNy1otcVxDU13defBk3nhJXQrlsMNCM/sAcKO7\nfywc/whwhbv/SSfz/j7wSeCd7n6qk+dvA24DqK6uvnxbZ5+8iPQLd+el5iNn+8LauucYZnBFzXiW\nzE2weE4VFaPjDZCoexrdtXsM5r2IKOI8PPU24P8ClwKlBI3Wx9x9dB8LjHR4ysxuCJf/Tnff3dP7\n6vCUSP5wdzbtOnp2D2Tz7qOYwfwp46irTVCXrCIxZnifl9PZF75ZcLipvR8p6DkUBvMp1T2JMzTW\nAB8C7gXmA38AzHD3v+5jgSUEDeELgNcIGsJ/z93XZ8xzGfBTYLG7b47yvgoNkfy1edcRGpuaaUil\nean5CABvqR5LfTJBXTLBpLEXHiDtX/jbtp0LjHZlZTB8OOzb98bXDYZG7DjEGhruPj/zFq9m9ri7\nvz2GIuuBfyPYe/m2uy8zs7uANe5+n5n9GkgC7Wdube+pAV6hITIwvLznKPc3NbNyXZoN6cMAzLt4\nLPXJKupqE1w8vuyC3rerQ1VdGQyny8Yh1l5ugRuAbxF0iZ4GbnX3N8dRaNwUGiIDz6t7j9HQFJyF\nlXrtEADJSWOoTyaoT1YxZcKIyO/V1VlQXbnQPY3BdigrztCYAuwiaM/4M2AM8DV33xJHoXFTaIgM\nbNv3HaexKWgDeWFnECBzLhodBkiCqeXdB0hXexoTJsCJE/E0dA/GRnPd7lVEBrwd+49zf1MzDU1p\nntt+EIBZVaNYEraBTK8Y+YbXdPeFDvHsHQzGi/zi3NO4muBajSmc32HhtD7WmBUKDZHB6bWDJ2hM\nBX1hrdl2AIBLKkdRl6yiPplgZuWos/Nm+9DRYOxOJM7QeIngsNRaznVYiLt3ch5C7ik0RAa/5kMn\naWwKeuN95tX9uMP0ipHU11ZRl0wwq2oUZpa15WtPo/s3esrdu7tSO68oNEQKy+7DJ1m1vpmVqTRP\nv7KfNodp5SPC03irmJ0YHXuAqE2j8zd4Szj4QYJTYv+L8/ueejaGOmOn0BApXHuOnGLV+uA6kCe3\n7qPNoWZCGXXJBPW1CWonxRcgOnvqjW/wSDevc3e//kKLyyaFhogA7Dt6igc37GJlKs3jL++jtc25\nePxw6muDRvQ3Tx6T1UNYA43OnhIRCR04dvpsgDy2ZS8tbc6kscODCwmTCeZNHktRUWEHSJxtGv8A\nfMndD4bj44D/5e6fjaXSmCk0RKQ7h46f4cEXd9GQSvPbzXs53dpGYswwFtdWsSSZ4C3V4woyQOIM\njefc/bIO055197d09ZpcUmiISFSHT57hoRd3sXJdM49u3sPpljYqRw8NOlOsrWJ+zXiKCyRA4ryf\nRrGZDW3vktzMhgP9c2cVEZEsGj1sCO+7bDLvu2wyR06e4eGXdrNyXZofPb2d7z7+KuUjh7K4tpL6\nZIIrp04omADpTpTQuBt4yMy+AzjwUeB7Wa1KRKSfjRo2hPfMm8R75k3i6KkWHnlpNw2pND9du5O7\nn9xO+chSFs0JDmFdOXU8JcVFuS45JyI1hJvZYoJOCw14wN1XZbuwC6XDUyISp+OnW1i9cQ8rU2ke\nfnE3J860Mn5EKTfOqaSuNsFVb5rAkEEQIDp7SkQkZidOt/KbTXtoSKV56MVdHDvdypjhQ1g0u5L6\nuQmuflM5pSUDM0AUGiIiWXTyTCuPbtpDY1Mzv96wiyOnWhg9rISFs6uoT1ZxzYxyhpYU57rMyBQa\nIiL95FRLK7/dvJeGVDMPbmjm8MkWRg0t4YbZldTVVvGOmRMZNiS/AyTW0AjPmKp2941xFJdNCg0R\nyaXTLW08tmUvDak0D2zYxaETZxg5tIQFl1ZQV5vg2kvyM0DivE7j3cA/A6XuPtXM5gF39XTb1VxR\naIhIvjjT2sbjL++jMZVm1fpmDhw/Q1lpMdfPqqA+GQRIWWmUk1izL87QWAtcD6xuv8gv837h+Uah\nISL5qKW1jSe37qehKc2qpmb2HTvN8CHFXDdrInW1Ca6fVcGIobkLkDgv7mtx90Pq2EtE5MKVFBdx\nzYxyrplRzl03zeHpV/bT2NRMY1MzDalmhpYUce0lE6lPBgEyatiQXJfcqSih0WRmv0dwZfgM4FPA\n49ktS0Rk8CopLuLt08t5+/RyPnfTHNa8uj8MjzSr1u+itKSId8yYyJK5VSy4tJLReRQgUQ5PlQF3\nAovCSauAv3f3k1mu7YLo8JSIDFRtbc7a7QdoSKVpTDXTfPgkpcVF/M6McuqSCd41N5G1RvQ42zQu\nc/fnYqssyxQaIjIYtLU5z+04SGMquK3tnqOnePZvFjIyS+0ecYbGI0ACuBf4sbuvj6fE7FBoiMhg\n4+5s23ecmvIRWVtG1NDo8Xp3d78OuBbYAyw3s5SZ5eW9NEREBiMzy2pg9EakTlLcvdndvwL8EfA8\n8LdZrUpERPJSj6FhZpea2efMrAn4KsGZU5OzXpmIiOSdKC0q3wF+BCxy99ezXI+IiOSxHkPD3d/W\nH4WIiEj+6zI0zOwed/+gmaUI7th39inA87UbERERyZ7u9jT+NHx8V38UIiIi+a/LhnB3T4eDt7v7\ntsw/4Pb+KU9ERPJJlFNuF3YyrS7uQkREJP9116bxCYI9imlmti7jqVHAY9kuTERE8k93exo/BN4N\n3Bc+tv9d7u6/H8fCzWyxmW00sy1m9plOnh9qZj8Jn3/KzGriWK6IZN+KFVBTA0VFweOKFbmu6Jzu\narvhBjA7/6+7+rO9nnm3Hd090h9QAVS3/0V9XTfvVwy8DEwDSoEXgNkd5rkd+EY4/CHgJz297+WX\nX+4iklt33+1eVuYO5/7KyoLpudZdbQsWnD+9p/qzvZ79uR2BNR7huzvq7V7/BbgI2A1MAV509zl9\nCSszuwr4nLvfGI7fEYbYFzLmWRXO84SZlQDNwETvpmh1WCiSezU1sG3bG6dPmQKvvtrf1Zyvu9o6\nm95xnsz6s72e/bkdY+uwEPh74G3AJnefCiwgnjaNScCOjPGd4bRO53H3FuAQMKHjG5nZbWa2xszW\n7NmzJ4bSRKQvtm/v3fT+1JfaOs6T7fXMx+0YJTTOuPs+oMjMitz9EWBeDMvu7P6xHfcgosyDuy93\n9/nuPn/ixIkxlCYifVFd3bvp/akvtXWcJ9vrmY/bMUpoHDSzkcCjwAoz+zLQEsOydwIXZ4xPBjr2\nbXV2nvDw1BhgfwzLFpEsWrYMysrOn1ZWFkzPte5qW7Cg69d1Vn+21zMvt2NPjR7ACIJG6xLgFoJ7\nhE+I0mDSw/uWAFuBqZxrCJ/TYZ4/5vyG8Ht6el81hIvkh7vvdp8yxd0seMyHRvB23dXWWWN4d/Vn\nez37azsSV0N4NplZPfBvBKH0bXdfZmZ3hcXfZ2bDgB8AlxHsYXzI3bd2955qCBcR6b2oDeHdXdx3\nhE46KuRch4Wj+1qkuzcADR2m/W3G8EngA31djoiIxKPL0HD3Uf1ZiIiI5L9It3s1s2vM7A/D4XIz\nm5rdskREJB9Fud3r/wb+CrgjnFQK3J3NokREJD9F2dN4H3ATcAzAg1u+6tCViEgBihIap8PTsRzA\nzEZktyQREclXUULjHjP7JjDWzD4O/Br4VnbLEhGRfNTd7V4BcPd/NrOFwGHgEuBv3f3BrFcmIiJ5\np8fQAAhD4kEAMys2s6Xunute3UVEpJ91eXjKzEab2R1m9lUzW2SBTxJ0/fHB/itRRETyRXd7Gj8A\nDgBPAB8D/oLgdNv3uPvz/VCbiIjkme5CY5q7JwHM7FvAXoI79h3pl8pERCTvdHf21Jn2AXdvBV5R\nYIiIFLbu9jTebGaHw2EDhofjsXVYKCIiA0t3HRYW92chIiKS/yJ1WCgiIgIKDRER6QWFhoiIRKbQ\nEBGRyBQaIiISmUJDREQiU2iIiEhkCg0REYlMoSEiIpEpNEREJDKFhoiIRKbQEBGRyBQaIiISmUJD\nREQiU2iIiEhkCg0RkQ5WrICaGigqCh5XrMh1Rfmjuzv3iYgUnBUr4Lbb4PjxYHzbtmAcYOnS3NWV\nL7SnISKS4c47zwVGu+PHg+mi0BAROc/27b2bXmhyEhpmNt7MHjSzzeHjuE7mmWdmT5jZejNbZ2Y3\n56JWESks1dW9m15ocrWn8RncEdWhAAAMaklEQVTgIXefATwUjnd0HPgDd58DLAb+zczG9mONIlKA\nli2DsrLzp5WVBdMld6HxHuB74fD3gPd2nMHdN7n75nD4dWA3MLHfKhSRgrR0KSxfDlOmgFnwuHy5\nGsHbmbv3/0LNDrr72IzxA+7+hkNUGc9fQRAuc9y9rZPnbwNuA6iurr5827ZtWahaRGTwMrO17j6/\np/mydsqtmf0aqOrkqV6dg2BmCeAHwC2dBQaAuy8HlgPMnz+//1NQRKRAZC003P2Grp4zs11mlnD3\ndBgKu7uYbzSwEvisuz+ZpVJFRCSiXLVp3AfcEg7fAvyi4wxmVgr8DPi+u9/bj7WJiEgXchUaXwQW\nmtlmYGE4jpnNN7NvhfN8EHgHcKuZPR/+zctNuSIiAjlqCM+m+fPn+5o1a3JdhojIgBK1IVxXhIuI\nSGQKDRERiUyhISIDgrorzw/qGl1E8p66K88f2tMQkbyn7srzh0JDRPKeuivPHwoNEcl76q48fyg0\nRCTvqbvy/KHQEJG8p+7K84fOnhKRAWHpUoVEPtCehoiIRKbQEBGRyBQaIiISmUJDREQiU2iIiEhk\nCg0REYlMoSEiIpEpNEREJLJBd7tXM9sDbMvBosuBvTlYbpy0DvlB65AfBvo69Lb+Ke4+saeZBl1o\n5IqZrYlyf918pnXID1qH/DDQ1yFb9evwlIiIRKbQEBGRyBQa8Vme6wJioHXID1qH/DDQ1yEr9atN\nQ0REItOehoiIRKbQ6AUzG29mD5rZ5vBxXCfzzDOzJ8xsvZmtM7ObM577rpm9YmbPh3/z+ncNYlmH\nqWb2VPj6n5hZaf+uQbR1COe738wOmtmvOkwfEJ9DOF9X6zCQPodbwnk2m9ktGdNXm9nGjM+hop/q\nXhwud4uZfaaT54eG23RLuI1rMp67I5y+0cxu7I96O3Oh62BmNWZ2ImObf6PXC3d3/UX8A74EfCYc\n/gzwj53MMxOYEQ5fBKSBseH4d4H3D/B1uAf4UDj8DeAT+bgO4XMLgHcDv+owfUB8Dj2sw4D4HIDx\nwNbwcVw4PC58bjUwv59rLgZeBqYBpcALwOwO89wOfCMc/hDwk3B4djj/UGBq+D7FOdjufVmHGqCp\nT8vv7xUeyH/ARiARDieAjRFe80LGF3A+fFld8DoARnCxUEk4/SpgVT6vA3BtnobGBa/DQPocgA8D\n38wY/ybw4XA4F6Fx3rYC7gDu6DDPKuCqcLgk3NbWcd7M+QbQOvQ5NHR4qncq3T0NED52uzttZlcQ\n/BJ4OWPysvCQz7+a2dDsldqlvqzDBOCgu7eET+8EJmWx1q70ah26MKA+hw4G0ucwCdiRMd6x1u+E\nh0n+xswse6VGrue8ecJtfIhgm0d5bX/oyzoATDWz58zsN2b2O71duO4R3oGZ/Rqo6uSpO3v5Pgng\nB8At7t4WTr4DaCb4El4O/BVw14VX2+Wys7IOXfynzsrpd3GtQxcG1OfQ2Vt3Mi1fP4fual3q7q+Z\n2SjgP4GPAN/vfZW9EmXbdTVPv233HvRlHdJAtbvvM7PLgZ+b2Rx3Pxx14QqNDtz9hq6eM7NdZpZw\n93T4hbq7i/lGAyuBz7r7kxnvnQ4HT5nZd4BPx1j6WVlch73AWDMrCX+9TAZej7l8IJ516Oa9B8zn\n0IWB9DnsJDi81m4ywWEp3P218PGImf0QuILsh8ZO4OIO9XTcdu3z7DSzEmAMsD/ia/vDBa+DB8er\nTgG4+1oze5mgDXNN1IXr8FTv3Ae0n/1xC/CLjjOEZ7H8DPi+u9/b4blE+GjAe4GmrFbbuQteh/Af\n3CPA+7t7fT/ocR26M1A+h64MsM9hFbDIzMaFZ1ctAlaZWYmZlQOY2RDgXfTP5/AMMCM8+6yUoJH4\nvg7zZK7X+4GHw21+H/Ch8MykqQTtfE/3Q80dXfA6mNlEMysGMLNpBOuwtVdL7+9GnIH8R3BM8CFg\nc/g4Ppw+H/hWOPz7wBng+Yy/eeFzDwMpgv8cdwMjB+A6TCP4j7IFuBcYmo/rEI7/N7AHOEHwy+vG\ngfQ59LAOA+lz+GhY5xbgD8NpI4C1wDpgPfBl+ulMJKAe2ETQTndnOO0u4KZweFi4TbeE23haxmvv\nDF+3Eajr723e13UAfjfc3i8AzwLv7u2ydUW4iIhEpsNTIiISmUJDREQiU2iIiEhkCg0REYlMoSEi\nIpEpNAqUmbWG3Tc0mdkvzWxshNcc7eH5sWZ2e8b4RWb20zjq7WJ5q80s5/dwzvZ6hsuYZ2b1Mb9n\nTrZf1HXpOJ+Z3dRZj67SvxQaheuEu89z91qCq13/OIb3HEvQuyYA7v66u7+/m/kHhX5az3kE5+ZH\nZoF8/D8edV3Om8/d73P3L2atKokkH/9BSf97gowOz8zsL8zsmbBDv7/rOLOZjTSzh8zsWTNLmdl7\nwqe+CLwp3IP5p7Dv/qbwNU+Z2ZyM91htZpeb2Qgz+3a4vOcy3qvjMv8yXNYLZpb5xfEBM3vazDa1\nd74WLve/w/qeNbO3h9OvDZf7UzN7ycxWtPenZWb14bTfmtlXLLx/RZT6OqznrWb2XxbcB2OzmX2p\ni/V51cz+LmMbzgqnX2Fmj4fLetzMLgmv+r0LuDnctjeb2efM7NMZ79cU1lFjZi+a2dcILt662My+\nbmZrLLg/yhs+zw511ZnZPRnj15rZL8PhRRbcZ+VZM7vXzEb2dtt1sS5R1/lWM/tq+N5Twn+D68LH\n6nD6d8MaHjezrWY26H+09LtcXdGov9z+AUfDx2KCK0cXh+OLCDrxM4IfFb8C3tHhNSXA6HC4nOCq\n0zd0u5w5DvwZ8HfhcALYFA7/A/D74fBYgqtcR3SotQ54HCgLx9uvPF4N/J9wuB74dThcBgwLh2cA\na8Lhawl6+5wcrtsTwDUEV8/uAKaG8/2IsCvyiPVlruetBN0yjAnfdxtwcSfb/1XgT8Lh2zl3Nf5o\nznV5fgPwnxnv+9WM138O+HTGeFNYRw3QBrwt47n27VUcbrO5Gdtvfoe6SoDt7esIfJ2gh4By4NGM\n6X8F/O2FbLtO1iXqOp8dB35J0JEmBFec/zwc/i7Bv+cigvtfbMn1/7XB9qc9jcI13MyeB/YR3CDn\nwXD6ovDvOYJfqrMIvngzGfAPZrYO+DXBXkplD8u7B/hAOPxBgv/Y7cv7TFjLaoIvoeoOr70B+I67\nHwdw9/0Zz/1X+LiW4AsTYAjwH2aWCpczO2P+p919pwc9Dz8fvmYWsNXdXwnn+VHG/FHq6+ghdz/k\n7ieBDcCULubrrPYxwL3hnsu/AnM6eV1PtnlGR5nAB83sWYLPdA7nb4/zeNAB4v3Auy3o6G4JQZ9S\nbwtf91i4LW4J1yuObXch63wV8MNw+AcE4d/u5+7e5u4b6PnfpfSSerktXCfcfZ6ZjSHYm/hj4CsE\ngfAFd/9mN69dCkwELnf3M2b2KsEXQpc86AJ7n5nNBW4G/mf4lAG/6+4bu3m50XUX1KfCx1bO/Xv+\nM2AX8GaCX5wnO5k/8zXd3cchSn1d1dSxrq7my5zn88Aj7v4+C27RubqL17Zw/uHlzO1/rH3Ago71\nPg281d0PmNl36eGzAn5C8O9hP/CMB73QGvCgu384c0Yzu6yb9+l025nZlR3mi7rO3cn895G5/fvj\nHh0FRXsaBc7dDwGfAj5tQW+jq4CPZhyvnmRvvHfzGGB3GBjXce6X9BFgVDeL+zHwl8AYd0+F01YB\nf5LRttDZl9ADYU1l4Tzje1itMUA63Jv4CMFhme68BEyzc/eCvjnjuSj1xWkM8Fo4fGvG9I7b9lXg\nLWFNbyG4/WhnRhOEyCEzqyQ41NeT1eF7f5wgQACeBK42s+nhMsvMbCYXtu06rkvUdc70OEHvrhD8\niPlthPWSGCg0BHd/jqDXyw+5+wMEu/1PhId3fsob/+OuAOab2RqC/7Avhe+zj+DwRZOZ/VMni/op\nwX/0ezKmfZ7gcNK68PDE5zup736Crp7XhIc6err/xdeAW8zsSYJ7BRzrbmZ3P0HQrnC/mf2WYC/l\nUNT6YvYl4Atm9hjnh90jwOz2RmGCmxaND7fHJwjaC97A3V8gOCy1Hvg28FhPBbh7K8HeZ134iLvv\nIfhC/1F4WPJJYNYFbruO6xJ1nTN9CvjDsJaPAH/a03pJPNTLrQjBGWHufjT8VfzvwGZ3/9dc1zUQ\naNsVFu1piAQ+Hv5qX09wuKS7Nh05n7ZdAdGehoiIRKY9DRERiUyhISIikSk0REQkMoWGiIhEptAQ\nEZHIFBoiIhLZ/wfFM8SvGQsfhAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Karp et al. use fourth-root transformation of the contamination data.\n",
    "\n",
    "ecoli4th = np.multiply(np.sign(ecoli), np.power(np.abs(ecoli), 1/4))   # signed fourth root\n",
    "coeffs4th = do_regress(ecoli4th, nature1)                              # reference estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated p: 0.00672 (0.0, 0.007245947095395319)\n"
     ]
    }
   ],
   "source": [
    "# simulated p-value\n",
    "vals = permute_regress(ecoli4th, nature1) # permutation distribution of estimated slope\n",
    "p_CI(vals, abs(coeffs4th[1])) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4th root without farm 31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.126\n",
      "Model:                            OLS   Adj. R-squared:                  0.091\n",
      "Method:                 Least Squares   F-statistic:                     3.599\n",
      "Date:                Sun, 16 Sep 2018   Prob (F-statistic):             0.0694\n",
      "Time:                        20:39:20   Log-Likelihood:                 9.6339\n",
      "No. Observations:                  27   AIC:                            -15.27\n",
      "Df Residuals:                      25   BIC:                            -12.68\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.1311      0.037      3.504      0.002       0.054       0.208\n",
      "x1            -1.1041      0.582     -1.897      0.069      -2.303       0.095\n",
      "==============================================================================\n",
      "Omnibus:                        0.695   Durbin-Watson:                   1.912\n",
      "Prob(Omnibus):                  0.707   Jarque-Bera (JB):                0.587\n",
      "Skew:                          -0.332   Prob(JB):                        0.746\n",
      "Kurtosis:                       2.716   Cond. No.                         17.2\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4HOWV7/HvkWTZlndZsjG2tRjL\nBDMYgxWMWexkzH4HmHlCwuIwMLnENyEJc5OHTGCYmZuQOCHJZJ2sDllIcBYgm5NMYICADTEQxGaW\nYFk2XjFYkvcVSzr3jyq12+2WVHJXq1vq3+d5+umq6re7TlXbOl31Vp3X3B0REZEoinIdgIiI9B9K\nGiIiEpmShoiIRKakISIikSlpiIhIZEoaIiISmZKGiIhEpqQhIiKRKWmIiEhkJbkOIG4VFRVeU1OT\n6zBERPqVZ555psXdK3tqN+CSRk1NDQ0NDbkOQ0SkXzGz9VHa6fSUiIhEpqQhIiKRKWmIiEhkShoi\nIhKZkoaIiESmpCEiIpEpaYiISGRKGiIiEpmShoiIRKakIVLgliyBmhooKgqelyzJdUSSzwZcGRER\niW7JEli4EPbtC+bXrw/mARYsyF1ckr90pCFSwG677XDC6LRvX7BcJB0lDZECtmFD75aLKGmIFLCq\nqt4tF8lp0jCzi8xslZk1mdktaV6/3syazez58HFDLuIUGagWLYKysiOXlZUFy0XSyVnSMLNi4JvA\nxcB04Gozm56m6S/cfWb4uLNPgxQZ4BYsgMWLoboazILnxYvVCS5dy+XVU2cATe6+FsDMfg5cDryS\nw5hECs6CBUoSEl0uT09NBDYmzW8Kl6V6l5mtNLP7zGxyug8ys4Vm1mBmDc3NzdmIVUREyG3SsDTL\nPGX+d0CNu88AHgLuSvdB7r7Y3evdvb6ysschbkVE5BjlMmlsApKPHCYBryc3cPdWdz8Yzn4PmNVH\nsYmISBq5TBpPA3VmVmtmpcBVwNLkBmY2IWn2MuCvfRifiIikyFlHuLu3mdmHgQeAYuAH7v6ymd0O\nNLj7UuAmM7sMaAO2AdfnKl4REQFzT+1G6N/q6+u9oaEh12GIiPQrZvaMu9f31E53hIuISGRKGiIi\nEpmShoiIRKakISIikSlpiIhIZEoaIiISmZKGiIhEpqQhIiKRKWmIiEhkShoiIhKZkoaIiESmpCEi\nIpEpaYiISGRKGiJSkJYsgZoaKCoKnpcsyXVE/YOSRi/sPdhG09bdDLRy8iKFZskSWLgQ1q8H9+B5\n4UIljihymjTM7CIzW2VmTWZ2SzftrjAzN7Mea71n06Ormjnvy8s5+44/ceuvVvLHF7ewc/+hXIYk\nIsfgtttg374jl+3bFyyX7uVsECYzKwYagfMJxgt/Grja3V9JaTcC+ANQCnzY3bsdYSmbgzBt3X2A\nh17ZyrLGraxoamX3wTaKi4zTJo9m3rRK5k6r5JSJoygqsqysX0TiUVQUHGGkMoOOjr6PJx9EHYQp\nZ8O9AmcATe6+FsDMfg5cDryS0u7TwBeAm/s2vKONGzGEa2ZXcc3sKg61d/Dchh0sb2xmWWMzX3qw\nkS892Ej5sFLOmVrB3GmVzK2rYNzIIbkOW0RSVFUFp6TSLZfu5TJpTAQ2Js1vAmYnNzCz04DJ7v57\nM8t50kg2qLiIM2rLOaO2nJsvPJHWPQd5vKmFZauaWb66haUvvA7ASRNGMndaBfOmVVJfXU5pibqR\nRHJt0aKgDyP5FFVZWbBcupfLpJHuHE7igNHMioCvANf3+EFmC4GFAFU5+qkwdvhgLp85kctnTqSj\nw/nrG7tY1tjMslXNfP+x1/jusrWUlRZz1gljw6OQSmoqhuUkVpFCt2BB8HzbbbBhQ3CEsWjR4eXS\ntVz2acwBPunuF4bztwK4++fC+VHAGmBP+JbjgG3AZd31a2SzT+NY7TnYxhNrWlnWuJVljc1s3LYf\ngOqxZcytq2TetErmnDCWYYNzmcNFgquH9Ie0MEXt08hl0igh6AifD2wm6Ai/xt1f7qL9o8DNuewI\nj4O7s651X6Iv5Ik1rew/1M6gYqO+ujw4CplWwfQJIzFTh7r0nc7LUFNP2SxerMRRCPI+aQCY2SXA\nV4Fi4AfuvsjMbgca3H1pSttHGQBJI9XBtnaeWbedZauDU1mvvrEbgMoRgzm3LugLObeukvJhpTmO\nVAa6mpr0ncPV1bBuXV9HI32tXySNbOhvSSPVm7sOsLwx6Ex/fHUz2/cdwgxmTBwVHoVUctrk0ZQU\nq0Nd4qXLUAtbxknDzC4ERrj7fSnLFwBb3f3BWCKNWX9PGsnaO5wXN+8Mr8hq5rkN2+lwGDG4hLOm\njmXetHHMnVbBpDFluQ5VBgAdaRS2OJLGk8Cl7t6csvw44NfuPieWSGM2kJJGqp37D7GiqYXl4ams\n13ceAOCEymGJo5Aza8cytLQ4x5FKf5SrPg11vueHOJLGSnef0dvXcm0gJ41k7s6a5j0sa2xhWWMz\nT61t5WBbB6UlRcyuLU/coV43brg61CWyvv4Drs73/BFH0mgEprt7W8ryQcAr7l4XS6QxK5SkkerA\noXaeem1b4lRW09bgSuUJo4YEl/WeWMnZJ1QwqmxQjiPt//TLOD46JZY/4kgadwDjCeo97Q2XDQO+\nDrS4+ydijDc2hZo0Um3esZ/Hwst6H29qYfeBNooMTqsaw9y64LLeGZNGU6w6Wb2iX8bxUud7/ogj\naZQAnwFuADp/C1QB3wf+3d3zsryrksbR2to7eH5jWCdrdQsrN+3AHUaXDUrUyZo3rZLxqpPVI/0y\njpf2Z/6I7ZJbMxsKTA1nm9x9fwzxZY2SRs+27X0rqU5WM827DwJw4vgRzDsxKHHy9toxDC5Rh3oq\n/TKOl47c8kccRxr/4u5fCKff7e73Jr32WXf/19iijZGSRu+4O3/dspvlq5tZ3tjM0+u2cajdGTqo\nmDOnHO5Qr60Ypg519Ms4G9RHlB/iSBrPuvvpqdPp5vOJkkZm9h5s48m1rSxrDJLIutbgJ+Dk8qFh\nX0glZ50wlhFDCrNDXb+MZaCKYzwN62I63bwMEMMGlzD/pPHMP2k8ABta94WFFlv4zXObWfLUBkqK\njNOrxzAv7AuZPmFkwQw8peqoUuh0pCGRvdXWwTPrtyduLnxlyy4AKoaXcm5Yrfecugoqhg/OcaQi\n0ltxnJ5qB/YSHFUMBToPyA0Y4u55eX5CSaPvbN19gMfCmwsfb2ph2963APibiSODvpC6Sk6vHsMg\n1ckSyXsqWCh9qqPDeen1w3Wynt2wg/YOZ/jgksTAU/OmVTK5XHWyRPJRHEca5d290d23HWNsWaWk\nkR92HTjEiqbWxKmszTuCK7WnVAxLJJDZU8opK9XAUyL5II6k8RrB8KsGTABe53AHuLv7lJhijZWS\nRv4J6mTtDUu+N/Pk2lYOHOqgtLiIt9eOSVzWe+L4EbqsVyRHYj09ZWbPuftpsUSWZUoa+e/AoXae\nXrctMXph45tBnazxIwcnLus9t66C0WUaeEqkr8SdNLJytZSZXQR8jWDkvjvd/Y6U1z8AfAhoJxgr\nfKG7v9LdZypp9D9bdu4PjkIaW3hsdTO7wjpZMyaNThyFzJysOlki2ZT3ScPMignGCD8f2EQwRvjV\nyUnBzEa6+65w+jLgRne/qLvPVdLo39raO3hh087EzYUvhHWyRg0N6mR1JpHjRqlOlkicMr65z8w+\nljQ7LmUed/9yBvEBnEFQy2ptuL6fA5cDiaTRmTBCwwj6WGQAKykuYlb1GGZVj+Fj509je1gnq/NU\n1h9e3ALAtPHDEwnk7TXlDBmkOlkifaG7S1dGJE1/L2U+DhOBjUnzm4DZqY3M7EPAx4BS4G9jjkHy\n3JhhpVx66vFceurxuDur3tydOJV114r1fO+x1xgyqIjZtWMTSeSEStXJEsmWnN2nYWbvBi509xvC\n+WuBM9z9I120vyZsf12a1xYCCwGqqqpmrU9XUU4GnH1vtfHU2m2JU1lrW/YCMHH00MRlvWdNHcvI\nAq2TJZkptEKKcVxye4+7vyec/nzyoEtm9j/ufkGGAc4BPunuF4bztwK4++e6aF8EbHf3Ud19rvo0\nCtfGbfsSCWTFmlb2HGyjuMg4vepwh/rfHD+qYOpkybErxMKUcSSNxGW2aWpPZXwJbjjIUyMwH9hM\n0BF+jbu/nNSmzt1Xh9OXAv+vp41S0hCAQ+0dPNtZJ6uxmZc2B91jY4eVck5d0KF+bl0llSNUJ0uO\nVogl8OOoctvdeauMz2m5e5uZfRh4gOCS2x+4+8tmdjvQ4O5LgQ+b2XnAIWA7cNSpKZF0BhUXMXvK\nWGZPGcvHL3wbLXsO8vjqlsSRyG+ffx2Ak48fydywTtas6jGUlqhOlgSnpHqzvJB0d6TxKnA1UATc\nDVxDcEe4AXe7+0l9FWRvZOtIo9DObw5kHR3OK1t2sSy8IuvZ9dtp63CGlRYz54QK5p1Yyby6SqrG\nqk5WodKRRjftukkaj3T3Rnd/5zHGllXZSBqFeH6zkOw+cIgVa1oTZU42bgvqZNWMLUschcw5YSzD\nBqtOVqEoxP/zqnIbo0L81VGo3J11rfsS94U8saaV/YfaGVRs1FeXJ8ZQP2mC6mQNNKlnEy65BP77\nvwvn7IKSRoyKiiDdbjKDjo5YVyV55mBbOw3rtieSyKtv7AZg3IjBwcBTJ1Zy7tQKxgxTnaz+6MYb\ng6OH9vajXxvoRxaplDRipCMN6fTGzgMsXx10pj/e1MKOfYcwgxkTRx1RJ6tEA0/lvRtvhG9/u/s2\nhfR/XEkjRoV4flN61t7hrNy0I3FF1vMbd9DhMGJICedMrQj6Q6ZVMnH00FyHKmmUlKQ/wkhWSGcT\n4i5YOBGoJukSXXdfnlGEWaKrpyRXdu47dESdrDd2HQCgbtzwRAKZXas6WfkiSpeUjjTStOspaZjZ\n54ErCQoJduZld/fLMo4yC3Rzn2RT1B8P7s7qrXsSCeSp17bxVlsHg0uKOKO2nHlhmZOp44arQz1H\nejrSKLSzCXEmjVXADHc/GFdw2aSkIdmSyWnK/W+189RrrSxvbGFZ41bWNAd1so4fNSRxFHL21ApG\nDVWdrL7SXZ9GdXXhnU2IM2n8EXi3u++JK7hsUtKQbInzgohN2/exvDE4lfXnphZ2h3WyZk4+XCfr\nlImjNPBUliVfPVVcHPwo+Na3ch1VbsSZNH4JnAo8DCSONtz9pkyDzAYlDcmWbF16fai9g+c37kic\nynpx807cYUzZIM6pC05jza2rYNxIDTwl2RNn0khb78nd7zrG2LJKSUOypa8uvW7dc5DHmzrrZLXQ\nsif4rfa240YkSpzMqhnD4BJ1qEt8dMmtSMxycel1R4fz1zd2JfpCnlm/nUPtTllpMXOmjE3coV5T\nMSw7AUjBiG08DTN7kTRVbd19RuZhxk9JQ7Ip15de7znYxhNJdbLWtwYZrKq8jLnTKphbV8lZUysY\nrjpZWZPrfwPZEkfSmODuW8ysOt3r7p6Xw+MpaUghWdeyN3GH+oo1rex7K6iTdXrVmMRRyPQJIzXw\nVEwG8o2+Oj0lUmDeauugYf228FRWM3/dEgw8VTF8MHPrgpLv50ytYOxwDTzVlZ6OIgZySSElDZEC\nt3X3gcRlvY+tbmZ7WCfrlImjmFsXXNZ7WtVoBg2QOlmZnjaKchQxkIuX9oukYWYXAV8jGLnvTne/\nI+X1jwE3AG1AM/C+nk6LKWmIHK29w3lp885EnaznNu6gvcMZMbiEs6aOZd60ccydVsGkMf1z4Kk4\nThtFOYrQkUb02lNDgSp3XxVHcOFnFhOMEX4+sIlgjPCr3f2VpDbvBJ5y931m9kHgHe5+ZXefq6Qh\n0rOd+w+xoqklGEN9VTOv7wzqZE2pHJa4ufDM2rEMLe0fl/XG8cc8ylGE+jS6HyO884MuBf4TKAVq\nzWwmcHsMtafOAJrcfW24np8DlxPUuALA3ZNHD3wSeG+G6xQRYNTQQVx8ygQuPmUC7s6a5j0sC/tC\nfvrUBn7453WUlhQxu7acueG4IXV5XCcrjjG9q6rSJ56qqsPTnYlhIF49FVWUm/ueAf4WeNTdTwuX\nrcz0klszuwK4yN1vCOevBWa7+4e7aP8N4A13/0ya1xYCCwGqqqpmrU/3zYtIJAcOtfOX17YlTmWt\n3hpUEDpu5JDgst5pQYf66LL8GXgqjiONgXwUEUVsRxpAm7vvzMIvjHQfmDaDmdl7gXpgXrrX3X0x\nsBiC01NxBShSiIYMKk4UUQR4fcf+xH0h97/0Bvc0bKLI4NSkOlmnThqd0zpZixal/4O/aNGR7brr\nLNdRRDRRjjS+T1B36hbgXcBNwCB3/0BGKzabA3zS3S8M528FcPfPpbQ7D/gvYJ67b+3pc9WnIZI9\nbe0dvLBpB8tWNbNsdQsrN+3APTjddU5dRaLk+/gc1MlKTgjl5cGybdsO//GHwj6S6EmctafKgNuA\nCwiODh4APu3uBzIMsISgI3w+sJmgI/wad385qc1pwH0Ep7FWR/lcJQ2RvrN971s8ljTwVPPuoE7W\nieNHJG4ufHtt39bJ6uo009Ch0Np6dPuBcOVTHPrLJbeXAF8luOT2B+6+yMxuBxrcfamZPQScAmwJ\n37Khpw54JQ2R3HB3Xn1jdyKBPL1uG4fanaGDijlzSnniVFZtxbCsdqh31b/RlYFwj0Uc4jzS+B1H\n9zXsBBqA72Z6xBE3JQ2R/LD3YBtPru2sk9XCay3BwFOTxgwN+kzqKjl76lhGDIl34KmuLp3tyrEe\naQy0GlRxdoSvBSqBn4XzVwJvAtOA7wHXHmuQIjJwDRtcwvyTxjP/pPEAbGjdx7KwTtZvn9vMT5/a\nQEnRkXWyTj4+8zpZXV06O3Ys7N/fc2d5FKmnwNavD+ahfyeOKKIcaSx397nplpnZy+5+clYj7CUd\naYjkv7faOnh2w/bEqayXXw/qZI0dVsq5dcFlvefWVVI5ovd1srq7dBbiOToYiHeGx3l66q/Ahe6+\nIZyvAu539+lm9lznvRv5QklDpP9p3n2Qx5uCu9MfW91C6963ADj5+JGJvpBZ1WMi18nK9qmjgViD\nKs6kcQnwHWANwdVTtcCNwKPA+939qxlHGyMlDZH+raPDefn1XYkSJ89s2E57hzN8cAlzThibuKx3\ncnnu6mTpSKPnDxsMvI0gabyab53fyZQ0RAaWXQcOsaKpNXGH+uYd+wGoreisk1XBmVPGUlbadwNP\nDcS7x+NOGmcBNSR1nLv7jzMJMFuUNEQGLndnbctelq0K7lB/cm0rBw51UFpcxNtrxyTqZJ04fkTW\n62QV6tVTUU5P/QQ4AXgeaA8Xu7vflHGUWaCkIVI4Dhxq5+l124LLehtbWPXmbgDGjxzMueGYIedO\nrWDMsPypk5Wv4u4In+79ZLQmJQ2RwvXGzgNBX0hjM4+vbmHn/mDgqRmTRod9IRWcOmk0JQNk4Kk4\nxZk07gVucvct3TbME0oaIgLBwFOddbKWr27mhY076HAYOaSEc+oqEqeyJowamutQ80KcSeMRYCbw\nF+Bg5/IYxtPICiUNEUlnx763eDysk7W8sYU3dgXX89SNG564rPeM2nKGDOofA0/FLc6k0VU58mXH\nGFtWKWmISE/cncY39yRuLvzLa9t4q72DIYOKmF07NpFETqjMbp2sfNIvChZmg5KGiPTWvrfaeGrt\n4YGn1oZ1siaODupkzZtWwVlTKxgZc52sfBLnkcaZBONZnEQw5GsxsNfdR8YRaNyUNEQkUxu37Uvc\nXLhiTSt7DrZRXGScXjWaueFVWadMHJVxnax8EmfSaACuAu4lGD3vH4E6d//XOAKNm5KGiMTpUHsH\nz23YkRi9cOWmnQCUDyvlnKkV4SiHFYwb0fcDT8Up1qTh7vXJ44Kb2Qp3PyumWGOlpCEi2dSy5yCP\nr25JJJGWPUGdrJMmjEzcoV5fXU5pSf+6rDfOpLEcOA+4E3iDYECk69391BiCvAj4GsEprzvd/Y6U\n1+cSDNI0A7jK3e/r6TOVNESkr3R0OK9sSaqTtX47bR3OsNJi5pwwNuwPqaR67LBch9qjOJNGNcH4\nGaXAR4FRwLfcvSnDAIsJhns9H9hEMNzr1e7+SlKbGmAkcDOwVElDRPLZnoNtPLGmlWWNW1nW2MzG\nbUGdrOqxZcFRSF0lc04Yy7DBfVcnK6rYBmFy985ajgeAT2UaWJIzgCZ3XwtgZj8HLgcSScPd14Wv\n9dNiwyJSSIYPLuH86eM5f/p43J11rftYtmory1e3cG/DJn78xHoGFRv11eWJgadOmpD9Ollx6jFp\nmNnZwCeBao4sWDglw3VPBDYmzW8CZmf4mSIiecHMqK0YRm1FLdefXcvBtnYa1h0eeOqOP77KHX98\nlcoRgzm3roJ50yo5Z2oFY4f3fuCpvhTlGOn7BKelnuFwwcI4pEutx3TTiJktBBYCVFVVZRKTiEhW\nDC4p5uypFZw9tYJbLzmJrbsOsHx1C8sam3nk1a386tnNmMEpE0clbi48bXL+1cmKkjR2uvsfs7Du\nTcDkpPlJwOvH8kHuvhhYDEGfRuahiYhk17iRQ7hi1iSumDWJ9g7nxc07E0ch33ykif/6UxMjhpRw\n9gkVwamsaZVMHJ37OlldJg0zOz2cfMTMvgj8iiNrTz2b4bqfBurMrBbYTHAvyDUZfqaISL9TXGTM\nnDyamZNHc9P8OnbuO8Sf17Qkii3e//IbADz77+dTnuMy790daXwpZT65V92Bv81kxe7eZmYfBh4g\nuOT2B+7+spndDjS4+1Izezvwa2AMcKmZfcrdT85kvSIi+W5U2SAuOWUCl5wyAXenaeseXti0M+cJ\nA1R7SkREiH7JbY89LGb2WTMbnTQ/xsw+k2mAIiLS/0Tplr/Y3Xd0zrj7duCS7IUkIiL5KkrSKDaz\nxIXDZjYUyO8LiUVEJCuiXHJ7N/Cwmf2QoAP8fcBdWY1KRETyUpQyIl8ws5UERQsN+LS7P5D1yERE\nJO9Eqprl7vcD92c5FhERyXP5dX+6iIjkNSUNERGJLFLSMLOhZnZitoMREZH8FuXmvkuB5wn7NMxs\nppktzXZgIiKSf6IcaXySYMCkHQDu/jxQk72QREQkX0VJGm3uvjPrkYiISN6LcsntS2Z2DcGd4XXA\nTcCK7IYlIiL5KMqRxkeAkwnG0vgpsBP4v9kMSkRE8lOUI40T3f024LZsByMiIvktypHGl83sVTP7\ntJlpACQRkQLWY9Jw93cC7wCagcVm9qKZ/VscKzezi8xslZk1mdktaV4fbGa/CF9/ysxq4liviIgc\nm0g397n7G+7+deADBPds/EemKzazYuCbwMXAdOBqM5ue0ux/A9vdfSrwFeDzma5XRPrGkiVQUwNF\nRcHzkiW5juiw7mI77zwwO/LRXfzZ3s6824/u3u0DOIngXo2XgGXAB4FxPb0vwufOAR5Imr8VuDWl\nzQPAnHC6BGghHKK2q8esWbNcRHLr7rvdy8rc4fCjrCxYnmvdxTZ//pHLe4o/29vZl/sRaPAIf7t7\nHCPczJ4Efgbc6+6vx5WszOwK4CJ3vyGcvxaY7e4fTmrzUthmUzi/JmzT0tXnaoxwkdyrqYH1649e\nXl0N69b1dTRH6i62dMtT2yTHn+3t7Mv9GHWM8CjjaZwZT0hHsXSrO4Y2mNlCYCFAVVVV5pGJSEY2\nbOjd8r6USWypbbK9nfm4H7vs0zCze8LnF81sZdLjxXBQpkxtAiYnzU8CUo9kEm3MrAQYBWxL/SB3\nX+zu9e5eX1lZGUNoIpKJrn675cNvukxiS22T7e3Mx/3YXUf4P4fPfwdcmvTonM/U00CdmdWaWSlw\nFZBaCHEpcF04fQXwJ+/pfJqI5NyiRVBWduSysrJgea51F9v8+V2/L1382d7OvNyPPXV6AJ+PsuxY\nHsAlQCOwBrgtXHY7cFk4PQS4F2gC/gJM6ekz1REukh/uvtu9utrdLHjOh07wTt3Flq4zvLv4s72d\nfbUfibEj/Fl3Pz1l2Up3nxFf6oqPOsJFRHov445wM/sgcCMwJaUPYwTw58xDFBGR/qa7q6d+CvwR\n+ByQfLf2bnc/qjNaREQGvi6ThgdjaOwErgYws3EEfQzDzWy4u+fBxXMiItKXIg33amargdcI7ghf\nR3AEIiIiBSZK7anPAGcCje5eC8xHfRoiIgUpStI45O6tQJGZFbn7I8DMLMclIiJ5KMogTDvMbDiw\nHFhiZluBtuyGJSIi+SjKkcblwH7go8D9BDfixXFHuIiI9DNRChbuTZq9K4uxiIhInuvu5r7dHFlR\n1sJ5A9zdR2Y5NhERyTPd3acxoi8DERGR/BdpuFczO8fM/imcrjCz2uyGJSIi+SjKzX3/D/gEwXCs\nAKXA3dkMSkRE8lOUI41/AC4D9gJ4MOSrTl2JiBSgKEnjrbDWugOY2bDshiQiIvkqStK4x8y+C4w2\ns/cDDwF3ZjcsERHJRz0mDXf/T+A+4JfAicB/uPvXM1mpmZWb2YNmtjp8HtNFu/vNbIeZ/T6T9YmI\nSDwiXT3l7g+6+8fd/WbgT2a2IMP13gI87O51wMMcOV5Hsi8C12a4LhERiUmXScPMRprZrWb2DTO7\nwAIfBtYC78lwvZdz+O7yu4C/T9fI3R8Gdme4LhERiUl3ZUR+AmwHngBuAD5OcLnt5e7+fIbrHe/u\nWwDcfUs4wNMxM7OFwEKAqqqqDEMTEZGudJc0prj7KQBmdifQAlS5e6Rf/mb2EHBcmpdu63WUPXD3\nxcBigPr6eu+huYiIHKPuksahzgl3bzez16ImjPA953X1mpm9aWYTwqOMCcDWqJ8rIiK5011H+Klm\ntit87AZmdE6b2a4M17sUuC6cvg74bYafJyIifaDLpOHuxe4+MnyMcPeSpOlMK9zeAZwfjj1+fjiP\nmdWHp8II5x8D7gXmm9kmM7sww/WKiEgGoozcF7tw+Nj5aZY3EHS6d86f25dxiYhI9yLdpyEiIgJK\nGiIi0gtKGiIiEpmShoiIRKakISIikSlpiIhIZEoaIiISmZKGiIhEpqQhIiKRKWmIiEhkShoiIhKZ\nkoaISIolS6CmBoqKguclS3IdUf7IScFCEZF8tWQJLFwI+/YF8+vXB/MACxbkLq58oSMNEZEkt912\nOGF02rcvWC5KGiIiR9iwoXfv46WcAAANcUlEQVTLC01OkoaZlZvZg2a2Onwek6bNTDN7wsxeNrOV\nZnZlLmIVkcJSVdW75YUmV0catwAPu3sd8HA4n2of8I/ufjJwEfBVMxvdhzGKSAFatAjKyo5cVlYW\nLJfcJY3LgbvC6buAv09t4O6N7r46nH4d2ApU9lmEIlKQFiyAxYuhuhrMgufFi9UJ3ilXV0+Nd/ct\nAO6+xczGddfYzM4ASoE1fRGciBS2BQuUJLqStaRhZg8Bx6V5qVfXIJjZBOAnwHXu3tFFm4XAQoAq\nnXgUEcmarCUNdz+vq9fM7E0zmxAeZUwgOPWUrt1I4A/Av7n7k92sazGwGKC+vt4zi1xERLqSqz6N\npcB14fR1wG9TG5hZKfBr4Mfufm8fxiYiIl3IVdK4AzjfzFYD54fzmFm9md0ZtnkPMBe43syeDx8z\ncxOuiIgAmPvAOptTX1/vDQ0NuQ5DRKRfMbNn3L2+p3a6I1xERCJT0hARkciUNEREJDIlDRERiUxJ\nQ0REIlPSEJF+QaPp5QeN3CcieU+j6eUPHWmISN7TaHr5Q0lDRPKeRtPLH0oaIpL3NJpe/lDSEJG8\np9H08oeShojkPY2mlz909ZSI9AsaTS8/6EhDREQiU9IQEZHIlDRERCSynCQNMys3swfNbHX4PCZN\nm2ozeyYcse9lM/tALmIVEZHDcnWkcQvwsLvXAQ+H86m2AGe5+0xgNnCLmR3fhzGKiEiKXCWNy4G7\nwum7gL9PbeDub7n7wXB2MDqVJiKSc7n6Qzze3bcAhM/j0jUys8lmthLYCHze3V/vwxhFRCRF1u7T\nMLOHgOPSvBS5xJi7bwRmhKelfmNm97n7m2nWtRBYCFClugIiIllj7t73KzVbBbzD3beY2QTgUXc/\nsYf3/BD4g7vf10O7ZmB9fNFGVgG05GC9cdI25AdtQ37o79vQ2/ir3b2yp0a5uiN8KXAdcEf4/NvU\nBmY2CWh19/3h1VVnA1/u6YOjbHQ2mFmDu9fnYt1x0TbkB21Dfujv25Ct+HPVp3EHcL6ZrQbOD+cx\ns3ozuzNscxLwlJm9ACwD/tPdX8xJtCIiAuToSMPdW4H5aZY3ADeE0w8CM/o4NBER6YYuY43P4lwH\nEANtQ37QNuSH/r4NWYk/Jx3hIiLSP+lIQ0REIlPS6IWINbNmmtkTYb2slWZ2ZdJrPzKz18J6Ws+b\n2cy+3YJYtqHWzJ4K3/8LMyvt2y2Itg1hu/vNbIeZ/T5leb/4HsJ2XW1Df/oergvbrDaz65KWP2pm\nq5K+h7Q3+WYh7ovC9TaZ2VEljMxscLhPm8J9XJP02q3h8lVmdmFfxJvOsW6DmdWY2f6kff6dXq/c\n3fWI+AC+ANwSTt9CcJd6aptpQF04fTxBDa3R4fyPgCv6+TbcA1wVTn8H+GA+bkP42nzgUuD3Kcv7\nxffQwzb0i+8BKAfWhs9jwukx4WuPAvV9HHMxsAaYApQCLwDTU9rcCHwnnL4K+EU4PT1sPxioDT+n\nOAf7PZNtqAFeymj9fb3B/fkBrAImhNMTgFUR3vNC0h/gfPhjdczbABjBzUIl4fI5wAP5vA3AO/I0\naRzzNvSn7wG4Gvhu0vx3gavD6VwkjSP2FXArcGtKmweAOeF0SbivLbVtcrt+tA0ZJw2dnuqdSDWz\nOpnZGQS/BNYkLV4UnvL5ipkNzl6oXcpkG8YCO9y9LXx5EzAxi7F2pVfb0IV+9T2k6E/fw0SC2nGd\nUmP9YXia5N/NzLIXauR4jmgT7uOdBPs8ynv7QibbAFBrZs+Z2TIzO7e3K9cY4SkshppZ4edMAH4C\nXOfuHeHiW4E3CP4ILwY+Adx+7NF2ue6sbEMX/6mzcvldXNvQhX71PaT76DTL8vV76C7WBe6+2cxG\nAL8ErgV+3PsoeyXKvuuqTZ/t9x5ksg1bgCp3bzWzWQQ1/U52911RV66kkcLdz+vqNTN708wm+OGa\nWVu7aDcS+APwb+7+ZNJnbwknD1pQS+vmGENPyOI2tACjzawk/PUyCchK5eE4tqGbz+4330MX+tP3\nsIng9FqnSQSnpXD3zeHzbjP7KXAG2U8am4DJKfGk7rvONpvMrAQYBWyL+N6+cMzb4MH5qoMA7v6M\nma0h6MNsiLpynZ7qnc6aWdB1zaxS4NfAj9393pTXJoTPRjCGyEtZjTa9Y96G8B/cI8AV3b2/D/S4\nDd3pL99DV/rZ9/AAcIGZjQmvrroAeMDMSsysAsDMBgF/R998D08DdeHVZ6UEncRLU9okb9cVwJ/C\nfb4UuCq8MqmWoJ/vL30Qc6pj3gYzqzSzYgAzm0KwDWt7tfa+7sTpzw+Cc4IPA6vD5/JweT1wZzj9\nXuAQ8HzSY2b42p+AFwn+c9wNDO+H2zCF4D9KE3AvMDgftyGcfwxoBvYT/PK6sD99Dz1sQ3/6Ht4X\nxtkE/FO4bBjwDLASeBn4Gn10JRJwCdBI0E93W7jsduCycHpIuE+bwn08Jem9t4XvWwVc3Nf7PNNt\nAN4V7u8XgGeBS3u7bt0RLiIiken0lIiIRKakISIikSlpiIhIZEoaIiISmZKGiIhEpqRRoMysPSzf\n8JKZ/c7MRkd4z54eXh9tZjcmzR9vZvfFEW8X63vUzHI+hnO2tzNcx0wzuyTmz8zJ/ou6LantzOyy\ndBVdpW8paRSu/e4+093/huBu1w/F8JmjCaprAuDur7v7Fd20HxD6aDtnElybH5kF8vH/eNRtOaKd\nuy919zuyFpVEko//oKTvPUFSwTMz+7iZPR0W9PtUamMzG25mD5vZs2b2opldHr50B3BCeATzxbB2\n/0vhe54ys5OTPuNRM5tlZsPM7Afh+p5L+qzUdf5LuK4XzCz5D8e7zewvZtbYWXwtXO9jYXzPmtlZ\n4fJ3hOu9z8xeNbMlnfW0zOyScNnjZvZ1C8eviBJfynZeb2a/smAcjNVm9oUutmedmX0qaR++LVx+\nhpmtCNe1wsxODO/6vR24Mty3V5rZJ83s5qTPeymMo8bM/mpm3yK4eWuymX3bzBosGB/lqO8zJa6L\nzeyepPl3mNnvwukLLBhn5Vkzu9fMhvd233WxLVG3+Xoz+0b42dXhv8GV4XNVuPxHYQwrzGytmQ34\nHy19Lld3NOqR2wewJ3wuJrhz9KJw/gKCIn5G8KPi98DclPeUACPD6QqCu06PKrucPA98FPhUOD0B\naAynPwu8N5weTXCX67CUWC8GVgBl4XznncePAl8Kpy8BHgqny4Ah4XQd0BBOv4Og2uekcNueAM4h\nuHt2I1AbtvsZYSnyiPElb+f1BGUZRoWfux6YnGb/rwM+Ek7fyOG78UdyuOT5ecAvkz73G0nv/yRw\nc9L8S2EcNUAHcGbSa537qzjcZzOS9l99SlwlwIbObQS+TVAhoAJYnrT8E8B/HMu+S7MtUbc5MQ/8\njqCQJgR3nP8mnP4Rwb/nIoLxL5py/X9toD10pFG4hprZ80ArwQA5D4bLLwgfzxH8Un0bwR/eZAZ8\n1sxWAg8RHKWM72F99wDvDqffQ/Afu3N9t4SxPErwR6gq5b3nAT90930A7r4t6bVfhc/PEPzBBBgE\nfM/MXgzXMz2p/V/cfZMHlYefD9/zNmCtu78WtvlZUvso8aV62N13uvsB4BWguot26WIfBdwbHrl8\nBTg5zft6st6TCmUC7zGzZwm+05M5cn8cwYMCiPcDl1pQ6O5/EdSUOjN835/DfXFduF1x7Ltj2eY5\nwE/D6Z8QJP9Ov3H3Dnd/hZ7/XUovqcpt4drv7jPNbBTB0cSHgK8TJITPuft3u3nvAqASmOXuh8xs\nHcEfhC55UAK71cxmAFcC/yd8yYB3ufuqbt5udF2C+mD43M7hf88fBd4ETiX4xXkgTfvk93Q3jkOU\n+LqKKTWurtolt/k08Ii7/4MFQ3Q+2sV72zjy9HLy/t/bOWFBYb2bgbe7+3Yz+xE9fFfALwj+PWwD\nnvagCq0BD7r71ckNzey0bj4n7b4zs9kp7aJuc3eS/30k7/++GKOjoOhIo8C5+07gJuBmC6qNPgC8\nL+l89UQ7euzmUcDWMGG8k8O/pHcDI7pZ3c+BfwFGufuL4bIHgI8k9S2k+yP0P2FMZWGb8h42axSw\nJTyauJbgtEx3XgWm2OGxoK9Mei1KfHEaBWwOp69PWp66b9cBp4cxnU4w/Gg6IwmSyE4zG09wqq8n\nj4af/X6CBALwJHC2mU0N11lmZtM4tn2Xui1RtznZCoLqrhD8iHk8wnZJDJQ0BHd/jqDq5VXu/j8E\nh/1PhKd37uPo/7hLgHozayD4D/tq+DmtBKcvXjKzL6ZZ1X0E/9HvSVr2aYLTSSvD0xOfThPf/QSl\nnhvCUx09jX/xLeA6M3uSYKyAvd01dvf9BP0K95vZ4wRHKTujxhezLwCfM7M/c2SyewSY3tkpTDBo\nUXm4Pz5I0F9wFHd/geC01MvAD4A/9xSAu7cTHH1eHD7j7s0Ef9B/Fp6WfBJ42zHuu9RtibrNyW4C\n/imM5Vrgn3vaLomHqtyKEFwR5u57wl/F3wRWu/tXch1Xf6B9V1h0pCESeH/4q/1lgtMl3fXpyJG0\n7wqIjjRERCQyHWmIiEhkShoiIhKZkoaIiESmpCEiIpEpaYiISGRKGiIiEtn/B0YSHOKnue3LAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sensitivity to farm 31\n",
    "ecoli4th_edit = np.array([float(f['E_Coli_Change']) for f in farms if f['Farm ID'] != '31']) # omit farm 31\n",
    "ecoli4th_edit = np.multiply(np.sign(ecoli4th_edit), np.power(np.abs(ecoli4th_edit), 1/4))    # signed 4th root\n",
    "coeffs4th_edit = do_regress(ecoli4th_edit, nature_edit1)           # reference coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated p: 0.06815 (0.0, 0.0697294654933772)\n"
     ]
    }
   ],
   "source": [
    "# simulated p-value\n",
    "vals = permute_regress(ecoli4th_edit, nature_edit1)  # permutation distribution of estimated slope\n",
    "p_CI(vals, abs(coeffs4th_edit[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions:\n",
    "\n",
    "1. How did Karp et al. decide to use the 4th-root transformation?\n",
    "\n",
    "1. What is the effect of using that transformation?\n",
    "\n",
    "1. How many other transformations were considered or tried? Does the article say? \n",
    "\n",
    "1. What effect does testing for normality and trying different transformations have on the p-value?\n",
    "\n",
    "Gelman and Loken, 2013: [\"Garden of Forking Paths\"](http://www.stat.columbia.edu/~gelman/research/unpublished/p_hacking.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Technical issues:\n",
    "\n",
    "1. Formally, $\\hat{\\beta} = (X^TX)^{-1}X^TY$, but that's not how to compute it. Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading extensions from ~/.ipython/extensions is deprecated. We recommend managing extensions like any other Python packages, in site-packages.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda/lib/python3.6/site-packages/IPython/core/formatters.py:839: FormatterWarning: JSON expects JSONable list/dict containers, not JSON strings\n",
      "  FormatterWarning)\n"
     ]
    },
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.6.1 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]"
        },
        {
         "module": "IPython",
         "version": "6.4.0"
        },
        {
         "module": "OS",
         "version": "Darwin 17.7.0 x86_64 i386 64bit"
        },
        {
         "module": "scipy",
         "version": "1.1.0"
        },
        {
         "module": "numpy",
         "version": "1.14.5"
        },
        {
         "module": "pandas",
         "version": "0.23.1"
        },
        {
         "module": "matplotlib",
         "version": "2.2.2"
        },
        {
         "module": "notebook",
         "version": "5.5.0"
        },
        {
         "module": "jupyter",
         "version": "1.0.0"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.6.1 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]</td></tr><tr><td>IPython</td><td>6.4.0</td></tr><tr><td>OS</td><td>Darwin 17.7.0 x86_64 i386 64bit</td></tr><tr><td>scipy</td><td>1.1.0</td></tr><tr><td>numpy</td><td>1.14.5</td></tr><tr><td>pandas</td><td>0.23.1</td></tr><tr><td>matplotlib</td><td>2.2.2</td></tr><tr><td>notebook</td><td>5.5.0</td></tr><tr><td>jupyter</td><td>1.0.0</td></tr><tr><td colspan='2'>Sun Sep 16 20:40:01 2018 PDT</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.6.1 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] \\\\ \\hline\n",
       "IPython & 6.4.0 \\\\ \\hline\n",
       "OS & Darwin 17.7.0 x86\\letterunderscore{}64 i386 64bit \\\\ \\hline\n",
       "scipy & 1.1.0 \\\\ \\hline\n",
       "numpy & 1.14.5 \\\\ \\hline\n",
       "pandas & 0.23.1 \\\\ \\hline\n",
       "matplotlib & 2.2.2 \\\\ \\hline\n",
       "notebook & 5.5.0 \\\\ \\hline\n",
       "jupyter & 1.0.0 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Sun Sep 16 20:40:01 2018 PDT} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.6.1 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]\n",
       "IPython 6.4.0\n",
       "OS Darwin 17.7.0 x86_64 i386 64bit\n",
       "scipy 1.1.0\n",
       "numpy 1.14.5\n",
       "pandas 0.23.1\n",
       "matplotlib 2.2.2\n",
       "notebook 5.5.0\n",
       "jupyter 1.0.0\n",
       "Sun Sep 16 20:40:01 2018 PDT"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%load_ext version_information\n",
    "%version_information scipy, numpy, pandas, matplotlib, notebook, jupyter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}