canineAge_regressionModeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. &nbsp;&nbsp;Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract .gz file content\n",
    "import pandas as pd\n",
    "df = pd.read_csv('canineAge_dataMatrix.txt.gz', compression='gzip', sep = '\\t', low_memory=False)\n",
    "df = df.transpose()\n",
    "\n",
    "# Recode species labels\n",
    "df[0] = df[0].apply(lambda x: 1 if x=='W' else 0 if x=='D' else x )\n",
    "\n",
    "# load dataset\n",
    "y = df[df[1]!='?'][1][1:].astype(float)\n",
    "X = df[df[1]!='?'].drop([1], axis = 1)[1:].astype(float)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data into training / testing sets\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, stratify=X[0], random_state=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ge metrics for regression model assessment\n",
    "from sklearn.metrics import r2_score, mean_squared_error\n",
    "def compute_metrics(y_test, y_pred):\n",
    "    r2 = r2_score(y_test, y_pred)\n",
    "    mse = mean_squared_error(y_test, y_pred)\n",
    "    return (r2, mse)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. &nbsp;&nbsp;PCA Feature Reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'explained_variance')"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAD3CAYAAADWiwWzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnhklEQVR4nO3deZRcdZ338fe3q7t6X9Lp7qSzJyQkICBLIBFQg6ATomNgdBRUHNyYjOCgj3McnfFx1PHMGZfHRxkVBhUBNwYZHfM4KCDKogwhCSYhISshSyedpNPpfV++zx/3dtM0vVSFrq6lP69z6nTVrXurvpeQT9/87m8xd0dERNJDVrILEBGR2Cm0RUTSiEJbRCSNKLRFRNKIQltEJI0otEVE0kh2sguYSBUVFb5gwYJklyEiEpfNmzefdPfKWPbNqNBesGABmzZtSnYZIiJxMbODse6r5hERkTSi0BYRSSMKbRGRNKLQFhFJIwptEZE0MuVDu7evn+bOHnr7+pNdiojIuKZ8aF/7nac47/MPs+Noc7JLEREZ15QP7fxoBIC27t4kVyIiMr4pH9oFYWh3dPcluRIRkfEptMPQbldoi0gamPKhnZ8TjOTXlbaIpIMpH9ovXWmrTVtEUp9CeyC0e3SlLSKpT6EdVfOIiKQPhbZuRIpIGpnyoZ2vNm0RSSNTPrR1pS0i6UShrdAWkTQy5UM7XzciRSSNTPnQVj9tEUknCQ9tM1ttZrvNbJ+ZfXqE983Mbgvf32ZmFw5574CZPWdmW8wsISv2qnlERNJJQldjN7MI8G3gzUANsNHM1rv780N2uxpYEj5WALeHPwdc4e4nE1XjYD9tDa4RkTSQ6CvtS4B97r7f3buB+4C1w/ZZC9zrgaeBMjOrTnBdg3SlLSLpJNGhPRs4POR1Tbgt1n0ceNjMNpvZTYkocLCfdpfatEUk9SW0eQSwEbZ5HPtc5u5HzawKeMTMdrn7Ey87OAjzmwDmzZsXd4EFOS/NPeLumI1UjohIakj0lXYNMHfI6znA0Vj3cfeBnyeAXxA0t7yMu9/p7svdfXllZWXcBWZHsohGsnCHrl6tEykiqS3Rob0RWGJmC80sClwHrB+2z3rg/WEvkpVAk7vXmlmhmRUDmFkh8BZgeyKKzFe7toikiYQ2j7h7r5ndAjwERIC73H2Hma0L378DeBBYA+wD2oEPhIfPAH4RNldkAz9x998kos6CaISmjh7au3spL4wm4itERCZEotu0cfcHCYJ56LY7hjx34OYRjtsPvDbR9YHWiRSR9DHlR0TCS3211TwiIqlOoY3atEUkfSi00fwjIpI+FNpoVKSIpA+FNpCfo+lZRSQ9KLRR84iIpA+FNlCQ+9JQdhGRVKbQBgrUPCIiaUKhjW5Eikj6UGijftoikj4U2uhGpIikD4U2ah4RkfSh0Abyo7oRKSLpQaENFKp5RETShEIb3YgUkfSh0OalqVk7NLhGRFKcQhvdiBSR9KHQ5qXmEd2IFJFUp9AGCnKC0G7r7iVY/UxEJDUptIHsSBbRSBbu0NXbn+xyRERGpdAOqQeJiKQDhXZIfbVFJB0otEO6GSki6SCu0Dazy83sA+HzSjNbmJiyJt9AX201j4hIKos5tM3sn4C/Bz4TbsoBfhTDcavNbLeZ7TOzT4/wvpnZbeH728zswmHvR8zsT2b2q1hrPR1q0xaRdBDPlfa1wNuBNgB3PwoUj3WAmUWAbwNXA2cD15vZ2cN2uxpYEj5uAm4f9v6twM446jwtAwNsOnrUpi0iqSue0O72oBOzA5hZYQzHXALsc/f97t4N3AesHbbPWuBeDzwNlJlZdfgdc4C3At+Lo87TMhDabV260haR1BVPaN9vZv9OEKofAX4LfHecY2YDh4e8rgm3xbrPN4BPAQnvPJ2vdSJFJA1kx7qju3/NzN4MNANLgc+5+yPjHGYjfVQs+5jZ24AT7r7ZzFaN+gVmNxE0qzBv3rxxyhldYa66/IlI6os5tMOeIk8OBLWZ5ZvZAnc/MMZhNcDcIa/nAEdj3OedwNvNbA2QB5SY2Y/c/X1DD3b3O4E7AZYvX37aY9AHb0Rqpj8RSWHxNI/8jJc3U/SF28ayEVhiZgvNLApcB6wfts964P1hL5KVQJO717r7Z9x9jrsvCI/73fDAnkgFah4RkTQQ85U2kB3eTATA3bvDIB6Vu/ea2S3AQ0AEuMvdd5jZuvD9O4AHgTXAPqAd+ECc5zAhND2riKSDeEK7zsze7u7rAcxsLXByvIPc/UGCYB667Y4hzx24eZzPeAx4LI5a46Z+2iKSDuIJ7XXAj83sWwQ3Dw8D709IVUlQoLlHRCQNxNN75AVgpZkVAebuLYkra/KpeURE0kE8vUdygXcAC4Bss6Cnnrt/MSGVTbLBdSIV2iKSwuJpHvkl0ARsBroSU07yqHlERNJBPKE9x91XJ6ySJNONSBFJB/H0037KzM5NWCVJNtg8osE1IpLC4rnSvhy40cxeJGgeMYIee+clpLJJphuRIpIO4gntqxNWRQrQyjUikg7i6fJ3EMDMqgjmAskoBTnh1Kzdvbg7A71jRERSSTwr17zdzPYCLwKPAweAXyeorkmXHckiGsnCHbp6Ez4TrIjIaYnnRuQ/AyuBPe6+ELgS+GNCqkqSGaW5AOw5nlHjhkQkg8QT2j3uXg9kmVmWu/8eOD8xZSXH5YsrAXhsd12SKxERGVk8od0YDmF/gmAOkm8CGTUSZdXSgdA+keRKRERGFk9orwU6gE8AvwFeAP48EUUly2WLK8iJGFsON9LY3j3+ASIikyzm0Hb3Nnfvc/ded7/H3W8Lm0syRlFuNsvnl9Pv8OTecWedFRGZdOOGtpn9IfzZYmbNQx4tZtac+BIn10tNJGrXFpHUM25ou/vl4c9idy8Z8ih295LElzi5Vi2tAuDxPXX095/2kpMiIgkRU/OImWWZ2fZEF5MKzpxRxMySPE62dvF8bcb9Q0JE0lxMoe3u/cBWM5uX4HqSzszUi0REUlY8vUeqgR1m9qiZrR94JKqwZFK7toikqngmjPpCwqpIMZctriAayWLzoQbqWrqoLM5NdkkiIkB8E0Y9nshCUklxXg6vX1LBo7tO8Jsdx7hh5fxklyQiAsQ3YdRKM9toZq1m1m1mfZnY5W/A1edWA/Dr52qTXImIyEviadP+FnA9sBfIBz4cbstIbz5rBjkR4+n99dS3ZtySmCKSpuIJbdx9HxAJR0b+AFiVkKpSQGlBDpctrqDf4aEdx5NdjogIEF9ot5tZFNhiZl8xs08AheMdZGarzWy3me0zs0+P8L6Z2W3h+9vM7MJwe56ZPWNmW81sh5lN+o3QNQNNJNvVRCIiqSGe0L4h3P8WoA2YC7xjrAPMLAJ8m2CpsrOB683s7GG7XQ0sCR83AbeH27uAN7n7awmmgF1tZivjqPdVe8vZM8jOMp56oZ5TbZpASkSSL57QvpBgId9md/+Cu/+vsLlkLJcA+9x9v7t3A/cRzBY41FrgXg88DZSZWXX4ujXcJyd8TOq48rKCKK87Yzp9/c4jzx+bzK8WERlRPKH9dmCPmf3QzN5qZrF0F5wNHB7yuibcFtM+ZhYxsy3ACeARd98w/AvM7CYz22Rmm+rqJn4wzFvDJpJHnle7togkXzxTs34AWAz8DHgP8IKZfW+cw0ZaHXf41fKo+4Q3PM8H5gCXmNk5I9R1p7svd/fllZWV45QTvxWLpgOw42jG9m4UkTQSb++RHoLFfO8DNvPKpo7hagjavgfMAY7Gu4+7NwKPAavjqXcizCsvIC8ni9qmTi2MICJJF8/gmtVmdjewD3gn8D2C+UjGshFYYmYLw54n1wHD5ytZD7w/7EWyEmhy91ozqzSzsvC784GrgF2x1jtRIlnG0hnFAOw6pgV/RSS54pl75EaCK+y/dveYRpu4e6+Z3QI8BESAu9x9h5mtC9+/A3gQWEPwy6Ad+EB4eDVwT9gDJQu4391/FUe9E2bpzGK21jSx+1gLK8PmEhGRZIhn7pHrxnrfzP7H3V83wnEPEgTz0G13DHnuwM0jHLcNuCDW+hJp2cxgrYddx9SuLSLJFVeb9jjyJvCzUsqyajWPiEhqmMjQzti1uQautHcfa9ESZCKSVBMZ2hmrvDBKVXEu7d19HG5oT3Y5IjKFTWRoj9TfOmMsnakmEhFJvokM7Rsm8LNSzlnV4c3IWoW2iCTPuL1HzKyFMdqr3b0k/JnRq7UvC6+0dx9XDxIRSZ5xQ9vdiwHM7IvAMeCHBE0h7wWKE1pdChlsHtGVtogkUTzNI3/m7t9x95Zwpr/bGWdq1kyyuKqISJZxoL6Nju6+ZJcjIlNUPKHdZ2bvDWfeyzKz9wJTJr1ysyMsqiik32HvCV1ti0hyxBPa7wHeBRwPH38ZbpsylulmpIgkWTxTsx5w97XuXuHule5+jbsfSGBtKWfgZuSzhxqSXImITFXxzPJ3ppk9ambbw9fnmdlnE1da6nnz2TMA+O9ttbR39ya5GhGZiuJpHvku8BmgBwYndBpzEqlMc+aMYs6fW0ZLVy8PPqflx0Rk8sUT2gXu/sywbVPucvPdFwfrNdy/8fA4e4qITLx4QvukmZ1BONDGzN4J1CakqhT256+dRUE0wjMHTrG/rnX8A0REJlA8oX0z8O/AMjM7Anwc+JtEFJXKinKzBxf7vX9TTZKrEZGpJp7eI/vd/SqgEljm7pdPtd4jAwaaSB7YXENPX3+SqxGRqSTmlWvMLJdgBOQCINssmNTP3b+YkMpS2EXzp3FGZSEv1LXx8I7jvPW88ZbKFBGZGPE0j/ySYPX1XqBtyGPKMTNuvHQBAF95aBddvVNmYKiIJFk8C/vOcffVCaskzVx3yTzufuoAL9S1ce9TB/nIGxYluyQRmQLiudJ+yszOTVglaSYnksVn33Y2ALf9bi/1rTEtUC8i8qrEE9qXA5vNbLeZbTOz58xsW6IKSwdXLK3ijWdW0tLZy//97Z5klyMiU0A8zSNXJ6yKNPbZt57FH/ad5CcbDlGcl8PNVyymKDee/6wiIrEb90rbzErCpy2jPKa0JTOK+fiVS+h3uP2xF7jia4/xyy1Hkl2WiGSoWJpHfhL+3AxsCn9uHvJ6TGa2OmxS2Wdmnx7hfTOz28L3t5nZheH2uWb2ezPbaWY7zOzWmM9qkn3syiX84qOXcv7cMupaurj1vi189aFduI+6SpuIyGmxRAaLmUWAPcCbgRpgI3C9uz8/ZJ81wMeANcAK4JvuvsLMqoFqd3/WzIoJfklcM/TY4ZYvX+6bNo37eyRh+vudnzxziH9av4O+fuddy+fwL9eeS3ZkItdPFpFMY2ab3X15LPvG1fhqZtOAJUDewDZ3f2KMQy4B9rn7/vD4+wj6eg8N3rXAvR789njazMrMrNrdawnnNnH3FjPbCcwedmxKycoy3rdyPrPK8vjoj5/l/k01tHX38a3rL2BgMJKIyKsRz3zaHwaeAB4CvhD+/Pw4h80Ghk6HVxNui2sfM1sAXABsiLXeZHrTshn8+MMrKc7L5r+31XKfZgQUkQkSz7/bbwUuBg66+xUEIVo3zjEjXV4Ob48Zcx8zKwL+E/i4uze/4gvMbjKzTWa2qa5uvHImz0Xzp/Gla84B4J9/9TwH66fk4FERmWDxhHanu3dCMA+Ju+8Clo5zTA0wd8jrOcDRWPcxsxyCwP6xu/98pC9w9zvdfbm7L6+srIz5ZCbD2vNn87bzqmnv7uOT92+lr183JkXk1YkntGvMrAz4L+ARM/slrwzg4TYCS8xsoZlFCVa6WT9sn/XA+8NeJCuBJnevtaAR+PvATnf/ehx1ppQvXXMOM0py2XSwge8+uT/Z5YhImotnatZr3b3R3T8P/G+CQL1mnGN6gVsI2r93Ave7+w4zW2dm68LdHgT2A/sIljT7aLj9MuAG4E1mtiV8rIn5zFJEWUGUL7/jPAC+9bt9NLR1J7kiEUln43b5M7Pysd5391MTWtGrkOwuf2O54fsbeHLvSf5m1Rn8/eplyS5HRFLIRHf520xwY3C0G4aa3i4Gn3zLUp7ce5K7/3iAD162kMri3GSXJCJpaNzmEXdf6O6Lwp/DHwrsGJ0/t4yrzqqio6ePOx5/IdnliEiaimuonpn9hZl93cz+j5ldk6CaMtYn3nwmAD96+iDHmzuTXI2IpKN4Btd8B1gHPAdsB9aZ2bcTVVgmes2sUq4+ZyZdvf186J6NHD7VnuySRCTNxHOl/Ubgz9z9B+7+A4K5QlYlpKoM9g9rzmJueT7bjzTztn/7A4/uPJ7skkQkjcQT2ruBeUNezwWm9CIIp2NueQG/uuX1XHVWFU0dPXzonk1863d7NSOgiMQkntCeDuw0s8fM7DGCiZuqzGy9mQ0fMCNjKC3I4c4blvOp1Usxg689vIdPPbCN7t7+ZJcmIikunln+PpewKqagrCzjo6sWc0ZlEbfe9yd+trmGg/XtfPDyBbzhzEoKolr9RkReKeb5tM3s7OFzWZvZKnd/LBGFnY5UHlwzlq2HG/nQPZs4GS4OnJudxVVnz+ATV53J4qqiJFcnIokWz+CaeEJ7O3Av8FWC+bS/Aix399edbqETLV1DG+BESyf/ufkIDz9/jD8dagQgkmW8++K5fOKqMzUYRySDJSq0C4EvAxcBxcCPgS+7e8o0xKZzaA9V29TBv/1uH/+x8TB9/c6MklweWHcpc8sLkl2aiCRAPKEdz43IHqADyCe40n4xlQI7k1SX5vMv157LQx9/PRfOK+N4cxc3fH8DdS1dyS5NRJIsntDeSBDay4HLgevN7IGEVCUALK4q5u4PXsJrZpVwoL6dv7rrGZo6epJdlogkUTyh/RFgL/AP7n6MYDHeLYkoSl5SkpfDPR+8hEUVhTxf28yqr/6eL/3qefadaE12aSKSBPG0ad8O9ANvcvezwkV+H3b3ixNZYDwypU17JEcaO/ibH21mW03T4Lbz5pTy1nOrWXNutdq7RdJYom5EPuvuF5rZn9z9gnDbFnc///RLnViZHNoA7s5zR5r46TOHWL/lKG3dfYPvzSsvYMXCci5fUsHV51QTzY5rLjARSaJEhfYG4FJgYxjelQRX2hecfqkTK9NDe6jOnj4e213Hfz9Xy+93naC1q3fwvdll+dx8xWLeedEchbdIGkhUaL8XeDdwIXAP8E7gs+7+s9MtdKJNpdAeqrevn+drm9mw/xT/senwYHv3zJI83nXxXN598Vxml+UnuUoRGU1CQjv84GXAlQSr2Dzq7jtPr8TEmKqhPVRfv/Pgc7Xc9uhe9obhbQYXzZvG+XPLOH9eGa+ZVcq88gIiWSMtRiQiky1hoZ3qFNovcXf+Z389P33mMA9tP0Z338u71OflZHFGZRFlBTnk50SoLM7l5isWM2eabmiKTLaJXiNS0pCZcekZFVx6RgVN7T08e7iBrYcb2Xq4kZ21LRxr7mTH0eaXHfPY7jp++KEVmu9EJIXpSnuKamrv4YWTrbR19dLR3cd3n9zPxgMNlBdGufeDl3DO7NJklygyZah5ROLW0d3Huh9t5vE9deTnRPiLC2dz3cXzOHeOwlsk0RTaclq6e/v51ANb+a8tRwe3zSrNo6I4l7KCKPPLC7hgXhnnzy1jYUUhZrqRKTIRUiq0zWw18E0gAnzP3f912PsWvr8GaAdudPdnw/fuAt4GnHD3c8b7LoX2xNh9rIX7Nh7iF386QmP7yHOdTC+McsnCci5eUM5rZpWwZEYx5YXRSa5UJDOkTGibWQTYA7wZqCGYdOr6oYspmNkagnlM1gArgG+6+4rwvTcArcC9Cu3J193bz5HGDhrau2lo62b38Ra2HGrk2UONgws2DFVRlMtF88tYsXA6KxaVc9bMErLUrVBkXKnUe+QSYJ+77wcws/uAtQTrSw5YSxDKDjxtZmVmVu3ute7+hJktSHCNMopodhYLKwpZSCEAV541Awi6Ex6ob+eZF+vZdKCBPSda2Xe8hZOtXTy04zgP7QhWmC8ryOF1i6bz2rllVBXnUlmcy7SCKEW52RTnZTOtIKpQF4lTokN7NnB4yOsagqvp8faZDdQmtjQ5XWYWhHlFIe++eB4QBPnB+naeefEUG148xdP76znS2MGvtx/j19uPjfg5hdEIS2cWs3RmCctmFnPmjGLOnFFEeWFU7eUio0h0aI/0N294e0ws+4z+BWY3ATcBzJs3L/bKZEKZGQsqCllQUci7Lp6Lu3P4VAdPvXCSfSdaqWvt4kRzF00dPbR29dLU0UNTRw/Phs0tQxXlZjNnWj5zpuUzqyx4zCzJY3pRlPLCKJXFuUwvzNWITpmSEh3aNcDcIa/nAEdPY59RufudwJ0QtGmfXpky0cyMedMLmDd99F+k9a1d7D7Wws5jLew93sLu4y3sO95KS1cvu461sOtYy6jHRrKMquJcqkryqC7JY2ZpHmfOKObc2aUsnVmsibIkYyU6tDcCS8xsIXAEuA54z7B91gO3hO3dK4Amd1fTyBQwvSiXSxfncuniisFt7k5TRw81DR3UNLRztLGTI40dnGjpor61i/rWbupauzjV1k1tUye1TZ1sHfa50UgW58wuYfmCcpbPn8bFC8qZpp4tkiESGtru3mtmtwAPEXT5u8vdd5jZuvD9O4AHCXqO7CPo8veBgePN7KfAKqDCzGqAf3L37yeyZkkuM6OsIEpZQXTMUZldvX2caO7ieHMnx5o7OdLQwY6jzWw/0sT+k22DzS53hvsvm1nMxQvKuWj+NC6aP4050/LVbi5pSYNrJOMEbeUNbD7QwDMHTrHlcCPdvS+fMCuancW0ghymFUSZXhSlsijo3TKzNJ/ZZXnMKsunujSfiiLdFJXES6UufyKTrjQ/hyuWVnHF0iogWDBi6+FGNh1s4NmDDTx7qIGG9h6ON3dxvHnsFe6j2VnMKMmlvDCX6YXRwZAvL4yydGYxr1s0nbycyGSclgig0JYpIC8nwopF01mxaPrgto7uvmDQUHt30E7e0sWJli5qmzo42tjBkcZOaps6aGzv4fCpDg6f6hjxswujEd64tJJFFUXkZmeRm5NFWX4Q6tMKo5TkZVOUl01ZfpT8qMJdXj2FtkxJ+dEI+dGgO+FY2rt7Od4c3PhsaOvmVFs39W1ByG94sZ4dR5t58LmR+6EPV1aQw+yyfKqKcynOy6E4L5vywuhgL5iKoqAtv7wgSml+jgYeyYgU2iJjKIhms7Aim4UVhSO+f/hUO4/vqeNUWzddvX109vTT0N49GPItXb20dvbS2N4z+NgRw/dmZxmV4SjSsoIoZfk5lORnk58TIT8nQkl+8Atg9rR8Zpbmqd/6FKLQFnkV5pYX8L6V88fdr7/fOdnWxdHGTk62dNHa1UtLZw/1bd0cb+7iRHMn9W3dNIaB39zZO9ilMRaRLKOiKMq5s8u48qwqrlxWRVVJ3qs9PUlBCm2RSZCVZVQV51FVHFuQdvb0UdfSxcnWLho7emhq76G5s4eO7j46evo41dbN0cYOaho6ON7cOeTG6nF+uzOY+yU3O4vivOAKvTy8gTq9KJfqkjyqy/KZVRr2kinLIzdb7e3pQqEtkoLyciLMLS9gbnlsa3Z29/ZT29TBH/fV8+jO4zz1Qj0dPX10tQbBv5+2UY81g/KCKBVFuUwvilIQzR68qTqjJC+cUqCA2WV5VJfmU5ir2Egm/dcXyQDR7CzmTy9k/vRC3rNiHu5OZ08/zZ09NHcEzTCnwhuoQbNLB7XhaNNjYdNMfVs3HB//u8oKcpgX/kKZO62AWWGYzyrLY3ZZPqX5OerbnkAKbZEMZGZhD5kIM0ryWDLGvr19/Zxq6+Zkazf1bV20d/fR3dtPR08fx5o6OXyqnZqGjqA7ZFNneEO1iW01TSN+XkE0wuyyfOaVFzBnWj5VJXlMLwyu5GdPC26eluTlJObEpwCFtsgUlx3JoqokL6Ybl+5OXWsXh091cOhUG0cagiCvbewYnCemtauXvSda2XuiddTPKS+McsHcMi6cP43XziljQUUB1aX56gETAw1jF5EJ4+40d/RyuKGdmobgCv1kazf1rcHgpSONwURgnT39rzg2GsmisjiX8sIoZQU5zCjJG+zWWJqfQ2E0GKg0d1p+xs25rmHsIpIUZkZpQQ6lBaWjTvjl7tQ0dLD5YAObDp5i97EWDta3D4b6kcaRR58OVZybzexp+RTlZlOQm01xbjYl+dmU5OdQUZjLjNI8qkvzOKOyKOPWLlVoi8ikMrPBnjHXXDB7cHt7dy91LV00tPfQEE69e6SxndrGTpo7e2nv7g2nFWgfnHM9FrNK8zh7VillBTnkRLLIz4lw3pxSViwqp7p07BGxqUihLSIpoSCazfzp2cyfPvZ+7k5Dew9HGzto7+6jrTsYdTqwGtLJ1i6ONXVytLGDvSdaOdrUydFRBilVFedSlJtNNDtomrl4QTkrF01nWXUxxbnZKdkEo9AWkbRiZpQXRmNq9ujrd1482cbuYy20dffS09dPY3sPmw6cYtOBBk6EE4UB7DrWwpN7Tw4em51lTCuMsvo1M/nMmmUURFMjLlOjChGRBIhkGYurilhcVfSK9/r6ndqmDrp6++ns6eNgfTsb9tez4cVTHDrVTnt3MCr1h08f5I8vnOS26y4Yc2GOyaLeIyIiI+jq7WP3sRb+7mdb2XO8lZyI8bE3LeGv37howof9x9N7RKufioiMIDc7wnlzylh/y+W8/3Xz6elzvv7IHlZ/40me3FuXtLoU2iIiY8jLifDFtefwkw+v4IzKQl482cYN33+G32yPbR71iabQFhGJwaWLK/j1rW/gxksXAPCzTYeTUodCW0QkRtHsLD666gwA/rDvJB3dfZNeg0JbRCQOVSV5vHZuGV29/fxh38nxD5hgCm0RkThdtawKgN8+H8NcthNMoS0iEqerzp4BwKO7TtDfP7ndphXaIiJxWjazmNll+Zxs7WJrTeOkfnfCQ9vMVpvZbjPbZ2afHuF9M7Pbwve3mdmFsR4rIpIMZsZVZwVNJI/uPDGp353Q0DazCPBt4GrgbOB6Mzt72G5XA0vCx03A7XEcKyKSFANNJAMLKU+WRM89cgmwz933A5jZfcBa4Pkh+6wF7vVgPP3TZlZmZtXAghiOFRFJihULp1OUm82uYy38xXf+ODgj4LffcyEzS8dfBeh0Jbp5ZDYwtAd6Tbgtln1iORYzu8nMNpnZprq65A0tFZGpJZqdxdXnzATg2UONbD7YwOaDDXT1JrbvdqKvtEeajHb4rdbR9onlWNz9TuBOCCaMirdAEZHT9aVrz+H6FfNe1oNkRgxrbb4aiQ7tGmDukNdzgKMx7hON4VgRkaTJzY5w4bxpk/qdiW4e2QgsMbOFZhYFrgPWD9tnPfD+sBfJSqDJ3WtjPFZEZEpJ6JW2u/ea2S3AQ0AEuMvdd5jZuvD9O4AHgTXAPqAd+MBYxyayXhGRVKdFEEREkkyLIIiIZCiFtohIGlFoi4ikkYxq0zazOuDgaRxaAUz+xLiTS+eY/jL9/GDqnuN8d6+M5eCMCu3TZWabYr0JkK50jukv088PdI6xUPOIiEgaUWiLiKQRhXbgzmQXMAl0jukv088PdI7jUpu2iEga0ZW2iEgaUWiLiKSRKR/ambYOpZnNNbPfm9lOM9thZreG28vN7BEz2xv+nNz5JBPAzCJm9icz+1X4OqPOMVzF6QEz2xX+eb4uk87RzD4R/j+63cx+amZ56X5+ZnaXmZ0ws+1Dto16Tmb2mTB7dpvZn8XyHVM6tDN0Hcpe4JPufhawErg5PKdPA4+6+xLg0fB1ursV2Dnkdaad4zeB37j7MuC1BOeaEedoZrOBvwWWu/s5BDN5Xkf6n9/dwOph20Y8p/Dv5XXAa8JjvhNm0pimdGgzZA1Ld+8GBtahTFvuXuvuz4bPWwj+os8mOK97wt3uAa5JSoETxMzmAG8Fvjdkc8aco5mVAG8Avg/g7t3u3kgGnSPB1ND5ZpYNFBAscpLW5+fuTwCnhm0e7ZzWAve5e5e7v0gwPfUl433HVA/tmNahTFdmtgC4ANgAzAgXlyD8WZXE0ibCN4BPAf1DtmXSOS4C6oAfhE1A3zOzQjLkHN39CPA14BBQS7D4ycNkyPkNM9o5nVb+TPXQjmkdynRkZkXAfwIfd/fmZNczkczsbcAJd9+c7FoSKBu4ELjd3S8A2ki/poJRhe26a4GFwCyg0Mzel9yqJt1p5c9UD+1Y1rBMO2aWQxDYP3b3n4ebj5tZdfh+NXAiWfVNgMuAt5vZAYImrTeZ2Y/IrHOsAWrcfUP4+gGCEM+Uc7wKeNHd69y9B/g5cCmZc35DjXZOp5U/Uz20M24dSjMzgnbQne7+9SFvrQf+Knz+V8AvJ7u2ieLun3H3Oe6+gODP7Hfu/j4y6xyPAYfNbGm46UrgeTLnHA8BK82sIPx/9kqC+y+Zcn5DjXZO64HrzCzXzBYCS4Bnxv00d5/SD4L1KfcALwD/mOx6JuB8Lif4J9Y2YEv4WANMJ7hzvTf8WZ7sWifofFcBvwqfZ9Q5AucDm8I/y/8CpmXSOQJfAHYB24EfArnpfn7ATwna6HsIrqQ/NNY5Af8YZs9u4OpYvkPD2EVE0shUbx4REUkrCm0RkTSi0BYRSSMKbRGRNKLQFklBZrbKzC5Ndh2SehTaIqlpFcFgE5GXUWhLyjKzBeGUpN8Np/B82MzyR9l3sZn91sy2mtmzZnaGBb4aTv35nJm9O9x3lZk9bmb3m9keM/tXM3uvmT0T7ndGuN/dZnaHmT0Z7ve2cHuemf0g3PdPZnZFuP1GM/u5mf0mnIbzK0Pqe4uZ/U9Y28/CaQYwswNm9oVw+3NmtiycM2Yd8Akz22JmrzezvwzPY6uZPZHQ//CS2pLdGV0PPUZ7AAsIppo9P3x9P/C+UfbdAFwbPs8jmDXuHcAjBNN+ziAYhVdNcBXbGD7PBY4AXwiPvRX4Rvj8buA3BBc3SwgGS+QBnwR+EO6zLPzcPOBGYD9QGr4+SDBMuQJ4AigMj/l74HPh8wPAx8LnHwW+Fz7/PPB3Q87vOWB2+Lws2X82eiTvoSttSXUvuvuW8PlmgiB/GTMrJgi0XwC4e6e7txOMDv2pu/e5+3HgceDi8LCNHkxj20UwIu3hcPtzw77jfnfvd/e9BIG8LPzcH4bftYsgnM8M93/U3ZvcvZNg2Pl8gnnNzwb+aGZbCIYyzx/yHQPzw4x4fqE/Aneb2UcIfgnJFJWd7AJExtE15HkfMFLzyEizpY21ffjn9g953c/L/14MHzLscXxuX/hZBjzi7tePc8zA/q/g7uvMbAXBHOJbzOx8d68fow7JULrSlrTnwdSzNWZ2DUA4AU8BQZPEuy1YlqySYFGB8Sfkebm/NLOssJ17EcEcEU8A7w2/60xgXrh9NE8Dl5nZ4vCYgvC4sbQAxQMvzOwMd9/g7p8DTvLy2eFkClFoS6a4AfhbM9sGPAXMBH5BMNnSVuB3wKc8mD0vHrsJmlV+DawLmz2+A0TM7DngP4Abw2aWEbl7HUF790/D+p4maGYZy/8Drh24EQl8NbxRuZ3gl8bWOM9DMoQmjBIZhZndTTCD4APJrkVkgK60RUTSiK60Ja2Y2bcJVq4Z6pvu/oNk1CMy2RTaIiJpRM0jIiJpRKEtIpJGFNoiImlEoS0ikkYU2iIiaUShLSKSRv4/jg0YdGq55JwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run PCA to address multicollinearity amongst predictor variables in regression model\n",
    "from sklearn import decomposition\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "pca = decomposition.PCA()\n",
    "\n",
    "# Plot & identify # of components to include\n",
    "pca.fit(X)\n",
    "plt.figure(1, figsize=(6, 4))\n",
    "plt.clf()\n",
    "plt.axes([.2, .2, .7, .7])\n",
    "plt.plot(pca.explained_variance_ratio_, linewidth=2)\n",
    "plt.axis('tight')\n",
    "plt.xlabel('n_components')\n",
    "plt.ylabel('explained_variance')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.&nbsp;&nbsp; Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-squared on test data = 0.30.  MSE on test data = 2.90\n"
     ]
    }
   ],
   "source": [
    "# Create linear regression object\n",
    "from sklearn.linear_model import LinearRegression\n",
    "lr = LinearRegression()\n",
    "\n",
    "# Train the model & make predictions\n",
    "pca = decomposition.PCA(n_components = 50, random_state = 0)\n",
    "pca.fit(X_train)\n",
    "y_pred_lin = lr.fit(pca.transform(X_train), y_train).predict(pca.transform(X_test))\n",
    "\n",
    "# Out metrics\n",
    "print(\"R-squared on test data = %.2f.  MSE on test data = %.2f\" % compute_metrics(y_test, y_pred_lin))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.&nbsp;&nbsp;Hyperparameter Tuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# Use GridSearch for 4-fold CV optimization of elastic net model parameters \n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.linear_model import ElasticNet\n",
    "\n",
    "enet = ElasticNet(max_iter = 10000, random_state = 0)\n",
    "param_grid = dict(alpha=[0.01, 0.1, 0.2, 0.5, 1.0], l1_ratio=[0.2, 0.1, 0.05, 0.01])\n",
    "grid_search = GridSearchCV(enet, scoring = 'r2', param_grid = param_grid, verbose = 10, cv = 4)\n",
    "grid_search.fit(X_train, y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ElasticNet(alpha=0.1, l1_ratio=0.01, max_iter=10000, random_state=0)\n"
     ]
    }
   ],
   "source": [
    "# Identify \"optimal\" values of alpha & L1_ratio\n",
    "print(grid_search.best_estimator_)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.&nbsp;&nbsp;Elastic Net Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-squared on test data = 0.41.  MSE on test data = 3.39\n"
     ]
    }
   ],
   "source": [
    "# run elastic net using those parameters\n",
    "enet = grid_search.best_estimator_.fit(X_train, y_train)\n",
    "y_pred_enet = enet.predict(X_test)\n",
    "\n",
    "print(\"R-squared on test data = %.2f.  MSE on test data = %.2f\" % compute_metrics(y_test, y_pred_enet))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.&nbsp;&nbsp;Bootstrapping for Enhanced Variable Selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "# Produce a bootstrap distribution for each\n",
    "from sklearn.linear_model import ElasticNet\n",
    "results = {}\n",
    "\n",
    "for i in range(20):\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, stratify=X[0], random_state=i)\n",
    "    lr = LinearRegression()\n",
    "    pca = decomposition.PCA(n_components = 50, random_state = 0)\n",
    "    enet = ElasticNet(alpha=0.2, l1_ratio=0.01, max_iter = 1000, random_state = 0)\n",
    "    \n",
    "    pca.fit(X_train)\n",
    "    y_pred_lin = lr.fit(pca.transform(X_train), y_train).predict(pca.transform(X_test))\n",
    "    y_pred_enet = enet.fit(X_train, y_train).predict(X_test)\n",
    "    r2_lin, mse_lin = compute_metrics(y_test, y_pred_lin)\n",
    "    r2_enet, mse_enet = compute_metrics(y_test, y_pred_enet)\n",
    "    y_quiz =  enet.predict(X_quiz)\n",
    "\n",
    "    results[i] = {'r2_lin':r2_lin, 'mse_lin': mse_lin,\n",
    "                 'r2_enet':r2_enet, 'mse_enet': mse_enet, 'y_quiz':y_quiz }\n",
    "    \n",
    "    print (i, results[i])\n",
    "      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ENet R-squared on test data = 0.42.  MSE on test data = 4.46\n",
      "LinRegression R-squared on test data = 0.34.  MSE on test data = 5.17\n"
     ]
    }
   ],
   "source": [
    "# Output results\n",
    "print(\"ENet R-squared on test data = %.2f.  MSE on test data = %.2f\" % (df_results['r2_enet'].mean(), df_results['mse_enet'].mean()))\n",
    "print(\"LinRegression R-squared on test data = %.2f.  MSE on test data = %.2f\" % (df_results['r2_lin'].mean(), df_results['mse_lin'].mean()))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.&nbsp;&nbsp;Distance-Weighted Prediction Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quantize the distribution of known ages\n",
    "import numpy as np\n",
    "vc = pd.Series(df[df[1]!='?'][1][1:].astype(float).values).value_counts()\n",
    "\n",
    "# Calculate the inverse square of the pairwise Euclidean distance between known & predicted age \n",
    "def quant_y(x):\n",
    "    pp =  np.argmin([(1.0/float(p))*(p-x)**2 for (p,v) in list(zip(vc.index, vc.values))])\n",
    "    return list(zip(vc.index, vc.values))[pp][0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-squared on test data = 0.58.  MSE on test data = 2.42\n"
     ]
    }
   ],
   "source": [
    "# Use the distance measures as weights to post-process predicted age\n",
    "y_pred_quantized = [quant_y(x) for x in y_pred_enet]\n",
    "\n",
    "print(\"R-squared on test data = %.2f.  MSE on test data = %.2f\" % compute_metrics(y_test, y_pred_quantized))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.&nbsp;&nbsp;Age Estimates Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_quiz = df[df[1]=='?'].drop([1], axis = 1)[1:].astype(float)\n",
    "y_quiz = [quant_y(x) for x in enet.predict(X_quiz)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 canid  age\n",
      "0                 005F  6.0\n",
      "1                 006M  4.0\n",
      "2                 007F  2.0\n",
      "3                 008M  4.0\n",
      "4  AFGH_AKC11_ODL13172  3.0\n",
      "5   AFGH_AKC3_ODL10621  2.0\n",
      "6   AFGH_AKC9_ODL10626  4.0\n",
      "7  AKIT_AKC10_ODL12758  4.0\n",
      "8   AKIT_AKC12_ODL1133  4.0\n"
     ]
    }
   ],
   "source": [
    "df_quiz = pd.DataFrame(list(zip(X_quiz.index, y_quiz)))\n",
    "df_quiz.columns = ['canid', 'age']\n",
    "print(df_quiz)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAEvCAYAAAAgi0SBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYWElEQVR4nO3df4xd5Z3f8fe3OFEdhvCjrG/A0A6tLLossyFh5LKNWs2EJTEGxWm1SUGUhU2i2ayCm628yjqNlK5UbRdpl2yzIgryBgpRU6ZRAgrCTsCimWWRkmxsCgysw+JS72Kb2stm42QCEp3m2z/mjHs93OuZuffce+e5vF/SaO5znufc5/nOHc/H59wzZyIzkSSpVH9n0AuQJKkbBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaOsGvYBWzj///BwdHT3Z/ulPf8qZZ545uAX1wLDVNGz1gDWVYNjqgeGrqa569u/f/0pm/lyrvjUZZKOjo+zbt+9ke2ZmhomJicEtqAeGraZhqwesqQTDVg8MX0111RMRf9muz1OLkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKKtyXst9tvozt2ntA/dft2AViJJWi2PyCRJRTPIJElFM8gkSUUzyCRJRTPIJElFM8gkSUUzyCRJRVs2yCLi4oj4dkQciIjnIuKT1fbzImJvRLxQfT63zf5bIuL5iDgYETvrLkCS9Oa2kiOyeWBHZv48cBXwiYi4DNgJPJaZm4DHqvYpIuIM4AvAtcBlwI3VvpIk1WLZIMvMlzPzyerxT4ADwEZgG3BfNew+4IMtdt8MHMzMFzPzdWC62k+SpFqs6j2yiBgF3gV8D2hk5suwEHbAhha7bAReamofrrZJklSLyMyVDYwYAf4E+N3MfCAifpSZ5zT1/21mnrtknw8B78/Mj1Xtm4HNmbm9xfNPAVMAjUbjyunp6ZN9c3NzjIyMrLa2FZs9cuKU9tjGs3s216Je19Rvw1YPWFMJhq0eGL6a6qpncnJyf2aOt+pb0U2DI+ItwNeBr2TmA9XmYxFxQWa+HBEXAMdb7HoYuLipfRFwtNUcmbkL2AUwPj6eExMTJ/tmZmZobtft1qU3Db6pd3Mt6nVN/TZs9YA1lWDY6oHhq6kf9azkqsUA7gYOZObnmroeAm6pHt8CfKPF7t8HNkXEJRHxVuCGaj9JkmqxkvfI3gPcDLw3Ip6qPrYCtwPXRMQLwDVVm4i4MCL2AGTmPHAb8AgLF4l8NTOf60EdkqQ3qWVPLWbmE0C06b66xfijwNam9h5gT6cLlCTpdLyzhySpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBtkAjO7czeyRE4zu3M3ozt2DXo4kFc0gkyQVbd1yAyLiHuB64HhmXl5t+2/ApdWQc4AfZeYVLfY9BPwE+L/AfGaO17JqSZIqywYZcC9wJ/DlxQ2Z+a8WH0fEHcCJ0+w/mZmvdLpASZJOZ9kgy8zHI2K0VV9EBPBh4L01r0uSpBXp9j2yfwYcy8wX2vQn8GhE7I+IqS7nkiTpDSIzlx+0cET28OJ7ZE3bvwgczMw72ux3YWYejYgNwF5ge2Y+3mbsFDAF0Gg0rpyenj7ZNzc3x8jIyMoq6sDskVPPjI5tPLtncy3O11gPx17rz3z90OvXaBCsae0btnpg+Gqqq57Jycn97a6z6DjIImIdcAS4MjMPr+A5fgeYy8w/WG7s+Ph47tu372R7ZmaGiYmJZdfZqaWXwB+6/bqezbU4346xee6YXdeX+fqh16/RIFjT2jds9cDw1VRXPRHRNsi6ObX4y8AP2oVYRJwZEWctPgbeBzzbxXySJL3BskEWEfcD3wEujYjDEfHRqusG4P4lYy+MiD1VswE8ERFPA38G7M7Mb9W3dEmSVnbV4o1ttt/aYttRYGv1+EXgnV2uT5Kk0/LOHpKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkohlkkqSiGWSSpKIZZJKkoi0bZBFxT0Qcj4hnm7b9TkQciYinqo+tbfbdEhHPR8TBiNhZ58IlSYKVHZHdC2xpsf0PM/OK6mPP0s6IOAP4AnAtcBlwY0Rc1s1iJUlaatkgy8zHgR928NybgYOZ+WJmvg5MA9s6eB5Jktrq5j2y2yLimerU47kt+jcCLzW1D1fbJEmqTWTm8oMiRoGHM/Pyqt0AXgES+A/ABZn5kSX7fAh4f2Z+rGrfDGzOzO1t5pgCpgAajcaV09PTJ/vm5uYYGRlZdXErNXvkxCntsY1n92yuxfka6+HYa/2Zrx96/RoNgjWtfcNWDwxfTXXVMzk5uT8zx1v1revkCTPz2OLjiPhj4OEWww4DFze1LwKOnuY5dwG7AMbHx3NiYuJk38zMDM3tut26c/cp7UM39W6uxfl2jM1zx+y6vszXD71+jQbBmta+YasHhq+mftTT0anFiLigqfkvgGdbDPs+sCkiLomItwI3AA91Mp8kSe0se0QWEfcDE8D5EXEY+PfARERcwcKpxUPAr1djLwS+lJlbM3M+Im4DHgHOAO7JzOd6UYQk6c1r2SDLzBtbbL67zdijwNam9h7gDZfmS5JUF+/sIUkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkq2rJBFhH3RMTxiHi2advvR8QPIuKZiHgwIs5ps++hiJiNiKciYl+N65YkCVjZEdm9wJYl2/YCl2fmLwJ/AXz6NPtPZuYVmTne2RIlSWpv2SDLzMeBHy7Z9mhmzlfN7wIX9WBtkiQtq473yD4CfLNNXwKPRsT+iJiqYS5Jkk4Rmbn8oIhR4OHMvHzJ9s8A48C/zBZPFBEXZubRiNjAwunI7dURXqs5poApgEajceX09PTJvrm5OUZGRlZc1GrNHjlxSnts49k9m2txvsZ6OPZaf+brh16/RoNgTWvfsNUDw1dTXfVMTk7ub/cW1bpOnzQibgGuB65uFWIAmXm0+nw8Ih4ENgMtgywzdwG7AMbHx3NiYuJk38zMDM3tut26c/cp7UM39W6uxfl2jM1zx+y6vszXD71+jQbBmta+YasHhq+mftTT0anFiNgC/Dbwgcx8tc2YMyPirMXHwPuAZ1uNlSSpUyu5/P5+4DvApRFxOCI+CtwJnAXsrS6tv6sae2FE7Kl2bQBPRMTTwJ8BuzPzWz2pQpL0prXsqcXMvLHF5rvbjD0KbK0evwi8s6vVSZK0DO/sIUkqWscXe0ijTRfJ7BibZ2JwS5H0JuYRmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloywZZRNwTEccj4tmmbedFxN6IeKH6fG6bfbdExPMRcTAidta5cEmSYGVHZPcCW5Zs2wk8lpmbgMeq9iki4gzgC8C1wGXAjRFxWVerlSRpiWWDLDMfB364ZPM24L7q8X3AB1vsuhk4mJkvZubrwHS1nyRJten0PbJGZr4MUH3e0GLMRuClpvbhapskSbWJzFx+UMQo8HBmXl61f5SZ5zT1/21mnrtknw8B78/Mj1Xtm4HNmbm9zRxTwBRAo9G4cnp6+mTf3NwcIyMjq6tsFWaPnDilPbbx7J7NtThfYz0ce60/8/VK89etsR42nFdmHe30+vtuEIatpmGrB4avprrqmZyc3J+Z46361nX4nMci4oLMfDkiLgCOtxhzGLi4qX0RcLTdE2bmLmAXwPj4eE5MTJzsm5mZobldt1t37j6lfeim3s21ON+OsXnumF3Xl/l6pfnrtmNsng/38DUahF5/3w3CsNU0bPXA8NXUj3o6PbX4EHBL9fgW4Bstxnwf2BQRl0TEW4Ebqv0kSarNSi6/vx/4DnBpRByOiI8CtwPXRMQLwDVVm4i4MCL2AGTmPHAb8AhwAPhqZj7XmzIkSW9Wy55azMwb23Rd3WLsUWBrU3sPsKfj1UmStAzv7CFJKppBJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKlqnt6jSGje69LZbt183oJVoOc2vla+TtHoekUmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKK1nGQRcSlEfFU08ePI+I3l4yZiIgTTWM+2/WKJUlq0vHfI8vM54ErACLiDOAI8GCLoX+amdd3Oo8kSadT16nFq4H/mZl/WdPzSZK0IpGZ3T9JxD3Ak5l555LtE8DXgcPAUeC3MvO5Ns8xBUwBNBqNK6enp0/2zc3NMTIy0vU625k9cuKU9tjGs3s21+J8jfVw7LXezdePmprnaKyHDef19uvWb73+vlvU/HXs9fdev2rql2GrB4avprrqmZyc3J+Z4636ug6yiHgrCyH1C5l5bEnf24GfZeZcRGwFPp+Zm5Z7zvHx8dy3b9/J9szMDBMTE12t83Sa/9Q89P7PzY/u3M2OsXnumF3Xs/n6UVPzHDvG5tl+07ba5xikXn/fLWr+Ovb6e69fNfXLsNUDw1dTXfVERNsgq+PU4rUsHI0dW9qRmT/OzLnq8R7gLRFxfg1zSpIE1BNkNwL3t+qIiHdERFSPN1fz/U0Nc0qSBHRx1SJARLwNuAb49aZtHwfIzLuAXwF+IyLmgdeAG7KON+UkSap0FWSZ+Srw95Zsu6vp8Z3AnUv3kySpLt7ZQ5JUtK6OyKR+6+cVfpLK4BGZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWhDe6/Ffv/VZ0nSYHhEJkkqmkEmSSqaQSZJKppBJkkqmkEmSSqaQSZJKppBJkkqWldBFhGHImI2Ip6KiH0t+iMi/igiDkbEMxHx7m7mkyRpqTp+IXoyM19p03ctsKn6+CfAF6vPkiTVotenFrcBX84F3wXOiYgLejynJOlNpNsgS+DRiNgfEVMt+jcCLzW1D1fbJEmqRWRm5ztHXJiZRyNiA7AX2J6Zjzf17wZ+LzOfqNqPAZ/KzP0tnmsKmAJoNBpXTk9Pn+ybm5tjZGRkVWubPXLilPbYxrNrGVuH2SMnaKyHY6/1br5+1NQ8R2M9bDivt1+3pXNa0+p18m9pLRu2emD4aqqrnsnJyf2ZOd6qr6v3yDLzaPX5eEQ8CGwGHm8achi4uKl9EXC0zXPtAnYBjI+P58TExMm+mZkZmtsrcevSmwbf1H7/1Yytw607d7NjbJ47Ztf1bL5+1NQ8x46xeT68yteo2zmtafU6+be0lg1bPTB8NfWjno5PLUbEmRFx1uJj4H3As0uGPQT8anX14lXAicx8uePVSpK0RDdHZA3gwYhYfJ7/mpnfioiPA2TmXcAeYCtwEHgV+LXulitJ0qk6DrLMfBF4Z4vtdzU9TuATnc4hSdJyvLOHJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWgGmSSpaAaZJKloBpkkqWjrBr0ASb03unM3ADvG5pkY7FKk2nlEJkkqWsdBFhEXR8S3I+JARDwXEZ9sMWYiIk5ExFPVx2e7W64kSafq5tTiPLAjM5+MiLOA/RGxNzP/fMm4P83M67uYR5Kktjo+IsvMlzPzyerxT4ADwMa6FiZJ0krU8h5ZRIwC7wK+16L7lyLi6Yj4ZkT8Qh3zSZK0KDKzuyeIGAH+BPjdzHxgSd/bgZ9l5lxEbAU+n5mb2jzPFDAF0Gg0rpyenj7ZNzc3x8jIyKrWNXvkxCntsY1n1zK2DrNHTtBYD8de6918/aipeY7GethwXm+/bkvntKbVz9Gvmvqlk58Na92w1VRXPZOTk/szc7xVX1dBFhFvAR4GHsnMz61g/CFgPDNfOd248fHx3Ldv38n2zMwMExMTq1rb4uXGiw7dfl0tY+swunM3O8bmuWN2Xc/m60dNzXPsGJtn+03bap/jdHNa0+rn6FdN/dLJz4a1bthqqqueiGgbZN1ctRjA3cCBdiEWEe+oxhERm6v5/qbTOSVJWqqbqxbfA9wMzEbEU9W2fwf8fYDMvAv4FeA3ImIeeA24Ibs9lylJUpOOgywznwBimTF3And2OockScvxzh6SpKIZZJKKNLpzN7NHTrzhwqa651j80NplkEmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKK1s2fcSmW902T6tXvP047rBb/6O6t1dez1K9j8/fDvVvO7Pl8HpFJkopmkEmSimaQSZKKZpBJkopmkEmSimaQSZKKZpBJkopmkEmSitZVkEXEloh4PiIORsTOFv0REX9U9T8TEe/uZj5JkpbqOMgi4gzgC8C1wGXAjRFx2ZJh1wKbqo8p4IudzidJUivdHJFtBg5m5ouZ+TowDWxbMmYb8OVc8F3gnIi4oIs5JUk6RTdBthF4qal9uNq22jGSJHUsMrOzHSM+BLw/Mz9WtW8GNmfm9qYxu4Hfy8wnqvZjwKcyc3+L55ti4fQjwKXA803d5wOvdLTQtWvYahq2esCaSjBs9cDw1VRXPf8gM3+uVUc3d78/DFzc1L4IONrBGAAycxewq1VfROzLzPHOl7r2DFtNw1YPWFMJhq0eGL6a+lFPN6cWvw9siohLIuKtwA3AQ0vGPAT8anX14lXAicx8uYs5JUk6RcdHZJk5HxG3AY8AZwD3ZOZzEfHxqv8uYA+wFTgIvAr8WvdLliTp/+vqD2tm5h4Wwqp5211NjxP4RDdzVFqecizcsNU0bPWANZVg2OqB4aup5/V0fLGHJElrgbeokiQVbc0H2XK3wSpJRFwcEd+OiAMR8VxEfHLQa6pLRJwREf8jIh4e9FrqEBHnRMTXIuIH1ev1S4NeUzci4t9W33PPRsT9EfF3B72m1YqIeyLieEQ827TtvIjYGxEvVJ/PHeQaV6tNTb9ffd89ExEPRsQ5A1ziqrSqp6nvtyIiI+L8uudd00G2wttglWQe2JGZPw9cBXyi8HqafRI4MOhF1OjzwLcy8x8D76Tg2iJiI/BvgPHMvJyFi7NuGOyqOnIvsGXJtp3AY5m5CXisapfkXt5Y017g8sz8ReAvgE/3e1FduJc31kNEXAxcA/xVLyZd00HGym6DVYzMfDkzn6we/4SFH47F3+kkIi4CrgO+NOi11CEi3g78c+BugMx8PTN/NNBFdW8dsD4i1gFvo83vc65lmfk48MMlm7cB91WP7wM+2M81datVTZn5aGbOV83vsvD7t0Vo8xoB/CHwKaAnF2Ws9SAb2ltcRcQo8C7gewNeSh3+EwvfpD8b8Drq8g+Bvwb+c3W69EsRceagF9WpzDwC/AEL/xt+mYXf53x0sKuqTWPxd1OrzxsGvJ66fQT45qAX0Y2I+ABwJDOf7tUcaz3IosW24i+zjIgR4OvAb2bmjwe9nm5ExPXA8Va3HSvYOuDdwBcz813ATynvlNVJ1ftG24BLgAuBMyPiXw92VVpORHyGhbcjvjLotXQqIt4GfAb4bC/nWetBtuJbXJUiIt7CQoh9JTMfGPR6avAe4AMRcYiFU7/vjYj/Mtglde0wcDgzF4+Wv8ZCsJXql4H/lZl/nZn/B3gA+KcDXlNdji3+RY3q8/EBr6cWEXELcD1wU5b9O1L/iIX/QD1d/Yy4CHgyIt5R5yRrPchWchusYkREsPC+y4HM/Nyg11OHzPx0Zl6UmaMsvD7/PTOL/t9+Zv5v4KWIuLTadDXw5wNcUrf+CrgqIt5WfQ9eTcEXryzxEHBL9fgW4BsDXEstImIL8NvABzLz1UGvpxuZOZuZGzJztPoZcRh4d/VvrDZrOsiqNzwXb4N1APhqZj432FV15T3AzSwctTxVfWwd9KLU0nbgKxHxDHAF8B8Hu5zOVUeWXwOeBGZZ+Hdf3N0jIuJ+4DvApRFxOCI+CtwOXBMRL7BwVdztg1zjarWp6U7gLGBv9TPirtM+yRrSpp7ez1v2Uask6c1uTR+RSZK0HINMklQ0g0ySVDSDTJJUNINMklQ0g0ySVDSDTJJUNINMklS0/we7Qg+wGjhFFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "plt.figure(1, figsize=(7, 5))\n",
    "df[df[1]!='?'][1][1:].astype(float).hist(bins=100)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:mypython3]",
   "language": "python",
   "name": "conda-env-mypython3-py"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}