Failed uncertainty estimate experiments

Failed_estimate_uncertainty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try two methods to predict uncertainty.\n",
    "\n",
    "Ensemble method `forestci`:\n",
    "- https://github.com/scikit-learn-contrib/forest-confidence-interval\n",
    "- Usage: http://contrib.scikit-learn.org/forest-confidence-interval/auto_examples/plot_mpg.html#sphx-glr-auto-examples-plot-mpg-py\n",
    "\n",
    "\n",
    "Conformal prediction `nonconformist`:\n",
    "- https://github.com/donlnz/nonconformist\n",
    "- Usage: https://github.com/donlnz/nonconformist/blob/master/README.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zhuangjw/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/sklearn/ensemble/weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.\n",
      "  from numpy.core.umath_tests import inner1d\n",
      "Failed to import duecredit due to No module named 'duecredit'\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "import forestci as fci\n",
    "\n",
    "from nonconformist.cp import IcpRegressor\n",
    "from nonconformist.nc import NcFactory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Frozen(SortedKeysDict(OrderedDict([('lon', 144), ('lat', 91), ('lev', 25), ('time', 20)])))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = xr.open_mfdataset('./nc/*.nc4').drop('KPP_RFactive')\n",
    "ds.dims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = ds.isel(time=0, lev=0, drop=True).to_dataframe().reset_index(drop=True).rename(columns=lambda s: s.replace('KPP_', ''))\n",
    "\n",
    "df_y = df.iloc[:, df.columns.str.startswith('AFTER_CHEM')].rename(columns=lambda s: s.replace('AFTER_CHEM_', ''))\n",
    "\n",
    "df_x = df.iloc[:, df.columns.str.startswith('BEFORE_CHEM')].rename(columns=lambda s: s.replace('BEFORE_CHEM_', ''))\n",
    "\n",
    "df_in = df.iloc[:, ~df.columns.str.startswith('AFTER_CHEM')].rename(columns=lambda s: s.replace('BEFORE_CHEM_', ''))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(13104, 1)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# scaling has no effect for decision tree\n",
    "# do this just for neural nets\n",
    "X_all= StandardScaler().fit_transform(df_in)\n",
    "X_all.shape\n",
    "\n",
    "# predict the difference\n",
    "# pick up a single variable for now\n",
    "Y_all = StandardScaler().fit_transform( (df_y-df_x)[['O3']]) \n",
    "Y_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, Y_train, Y_test = train_test_split(X_all, Y_all, test_size=0.2, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10483, 127)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic random forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_rf = RandomForestRegressor(n_estimators=20, max_leaf_nodes=1000, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zhuangjw/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/ipykernel_launcher.py:1: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building tree 1 of 20\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    1.0s remaining:    0.0s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building tree 2 of 20\n",
      "building tree 3 of 20\n",
      "building tree 4 of 20\n",
      "building tree 5 of 20\n",
      "building tree 6 of 20\n",
      "building tree 7 of 20\n",
      "building tree 8 of 20\n",
      "building tree 9 of 20\n",
      "building tree 10 of 20\n",
      "building tree 11 of 20\n",
      "building tree 12 of 20\n",
      "building tree 13 of 20\n",
      "building tree 14 of 20\n",
      "building tree 15 of 20\n",
      "building tree 16 of 20\n",
      "building tree 17 of 20\n",
      "building tree 18 of 20\n",
      "building tree 19 of 20\n",
      "building tree 20 of 20\n",
      "CPU times: user 18.9 s, sys: 206 ms, total: 19.1 s\n",
      "Wall time: 19.4 s\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done  20 out of  20 | elapsed:   19.4s finished\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,\n",
       "           max_features='auto', max_leaf_nodes=1000,\n",
       "           min_impurity_decrease=0.0, min_impurity_split=None,\n",
       "           min_samples_leaf=1, min_samples_split=2,\n",
       "           min_weight_fraction_leaf=0.0, n_estimators=20, n_jobs=1,\n",
       "           oob_score=False, random_state=None, verbose=2, warm_start=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time model_rf.fit(X_train, Y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1c180c77b8>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAD8CAYAAAC/1zkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm0HVWZ9/HvjzAkl0AiMxkkCLRAAga4gt2NItBgumVUkNAi4cXutAivNCovqGhHXGhkEFEwGJTJhQEHplaw6RXAlgbBmxAyMBgiEULCFKaEEDTJ8/5R+ySVc8+995zkDHXv/X3WOitVu3bV2acW5MmuXfvZigjMzMzqYZNWN8DMzPoOBxUzM6sbBxUzM6sbBxUzM6sbBxUzM6sbBxUzM6sbBxUzM6sbBxUzM6sbBxUzM6ubTVvdgGbbbrvtYtSoUa1uhplZrzJjxoxXImL7nur1u6AyatQoOjo6Wt0MM7NeRdKfq6nXtMdfkpaX7Z8m6cq0/XlJj0uaLWm6pF1y9UZLulfSHyXNl/RVSUrHPpnOmS3pQUnva9bvMTOzzorSU3kUaI+IFZLOAC4GTpI0CLgTOCMi7pHUBvwS+CxwFfAMcEhEvCbpH4GpwEHdfdGyZXOYfu9u65UdftiCuv8gM7P+qBAD9RFxX0SsSLu/B0ak7X8G/jci7kn1VgBnAeen/Qcj4rUK55mZWQs0s6cySNKs3P42ZL2Qcp8G7k7bo4EZ+YMRsUDSYElbR8SbXZxnZmYt0Myg8nZEjC3tSDoNaM9XkHRKKjukVAR0teDL2nJJh5IFlYMrVZQ0EZgIsMMORXniZ2bW9xTi8ReApH8AvgIcExHvpOJ5dA487wGWR8SytL8v8CPg2IhYWunaETE1Itojon3o0ML8ZDOzPqcQf8NK2g/4IVlAeSl36Cbg4BRwSAP33yMbyEfSu4FbgU9FxB+b22ozMytXlGdBlwCDgZ+nt4WfjYhjIuJtSccC35d0FTAA+AlwZTrva8C2wA/Seasior3T1XO22mofDj/M81TMzBpB/W2N+vb29vDkRzOz2kia0dM/2qEgj7/MzKxvcFAxM7O6cVAxM7O6qWqgPuXtehkYFxFP5cq/CyyOiNLbWFcAJwAjI2JNKjuNLAXLWVV8zwHA9cAg4C7g7IgISdsAtwCjgIXAJ1Jqlj2B64D9ga9ExKU9fcfixYuZNGlSp/JKZWZmVptaeio3A+NLO5I2IQsgt+T2jweeAz60ge2ZQjZJcY/0GZfKzwemR8QewPS0D/Aq8Dmgx2BiZmaNV0tQmUYuqJAFjoURUUqHfCgwlywwnFxrQyTtDGwdEQ9F9krajcBx6fCxwA1p+4ZSeUS8FBF/AP5a6/eZmVn9VR1UImI2sCaXXn48WaApOTnt3wYcJWmzGtsyHFiU21+UygB2jIglqR1LgB1qubCkiZI6JHWsWLGi5xPMzGyD1DpQPw0YL2lTst7DzwEkbQ78E3B7SvL4MHBkjddWhbK6TKLJp2lpa2urxyXNzKyCWmfUTwPuAX4LzM6lVBkHDAHmpJntbcAK4Nc1XHsR66euHwEsTtsvSto5Ipakx2QvdTq7SsOGDfOgvJlZg9TUU4mIBcBSYDKdH339S0SMiohRwK7AkWlRrWqvvQRYJukDaWXHU4E70uE7gQlpe0Ku3MzMCqTHoJIedb2TK5oG7Ek2dkIKHB8h1yuJiLeAB4CjU9FpkhblPl0tpnUGWcbhp4EFrFsfZTJwhKT5wBFpH0k7SVoEfB64IF17655/tpmZNUKPub/SwPw1EXFgc5rUWM79ZWZWu7rk/pL0GbKeyQX1apiZmfVd3QaViLg6IvYurRG/odKM/Pz+AkkvS5ol6QVJayQdUqm+pBGS7pA0P513RXrbLH+9d0taLumLG9NOMzPbOK1aT+UbpNQtkiYBp5O9kvzbfKU0YH8rMCUijpU0AJgKXAScm6t6OVWuT/+X55ez6PzfdSofMfmDG/AzzMwsrygJJa8FTko5vvIOA1ZGxHUAEbEaOAc4vfRmmaTjgD+RLT1sZmYt1KygMig96polaRZwYdnx5WSB5eyy8tHAjHxBmlz5LLC7pC2B84CvN6bZZmZWi2Y9/no7IsaWdkqZi8vqfA+YJemyXJmoPKu+VP514PKIWJ4mXVYkaSJZokqGb73jhrTfzMyqUJQ16omI1yX9FPhsrnge8PF8vTQPZSTZPJaDgBMkXQwMJctNtjIirsyfExFTycZi2HfnPfvX+slmZk1UmKCSfAf4A+vaNR2YLOnUiLgxDdRfBlwfESuAtaPracB/eXlAKbf58MEelDcza5CiDNQDEBGvkM3U3yLtB9kaLSem2fR/BFYCX25ZI83MrEs9zqjvazyj3sysdnWZUW9mZlYLBxUzM6ubQgQVSSHpJ7n9TVMal1+V1btD0kNlZZMkPZ/mwMyVdEyz2m1mZusryttfbwFjJA2KiLfJ0ts/n68gaSiwP7Bc0q4R8Uzu8OURcamkvYDfSdohItZU+qIX//Q0l510VKfyL9zyqwq1zcysFoXoqSR3Ax9N26X17vM+DvwncDMwvtIFIuIJYBWwXYPaaGZm3ShSULkZGC9pILAv2Tr3eaVAMy1tdyLpIGAN8HJZ+URJHZI63nrnL3VvuJmZZYry+IuImC1pFFnAuCt/TNKOwO7AAxERklZJGhMRc1OVcySdAiwDToqy96TzM+pHbjO0f71DbWbWREXqqUC2Fv2ldH70dRLwLuAZSQuBUaz/COzyiBgbER+MiM557c3MrCkK01NJrgXeiIg5kj6cKz8ZGBcRDwFI2hX4bzZgRcod37O7B+XNzBqkUD2ViFgUEVfky9IjsXcDv8/VewZ4M42hmJlZQRSipxIRgyuU3Q/cn3aHVzi+f9osH9A3M7MWKVRPxczMejcHFTMzq5vCBRVJq1PKlcckzZT0d93UHZFSt8yXtEDSFZI2b2Z7zcxsnUKMqZRZu/SwpI8A3wIOyVdIi3WtAW4FpkTEsalsKnARcG5XF3/pz8u46jP3Vjx25tWH1eUHmJn1V4XrqZTZGngNQNKHJd2XlhyeAxwGrIyI6wAiYjVwDnC6pLZWNdjMrD8rYk9lkKRZwEBgZ7LgUXIgMCYinpH0OWBG/sSIeFPSs2Sz72eXyiVNBCYCvGvwDg1uvplZ/1XEnsrbaXb8nsA44EZJSsceyWUnFlAp5Uqn8oiYGhHtEdE+eODQhjXczKy/K2JQWSvNoN8O2D4VvZU7PA9Yb2lLSVsDI4EFTWmgmZmtp4iPv9aStCcwAFha4fB0YLKkUyPixjRQfxlwfUSs6OqaO+yylQfkzcwapIg9lUHpleJZwC3AhDQIv56Uifh44ERJ84E/AiuBLze1tWZmtlbheioRMaCL8vtZl7alVPYccHTjW2VmZtUoYk/FzMx6qYYHFUnL05+bSPqepLmS5kj6Q0phj6Qhkm5Ms+IXpO0h6dgoSW+nR2KPS7o6XWuUpJD0f3PfdaWk0xr9m8zMrLJmPv46CRgG7BsRaySNYN3bXD8G5kbEqQCSvg78CDgxHV8QEWMlbQrcCxwHzAReAs6W9MOIqGqd4JVz5/HEnnt1W2evJ5+o7ZeZmRnQ3MdfOwNLImINrF075TVJuwMHAN/I1b0QaJe0W/4CEbEKeJBsciNka9FPByY0uvFmZtazZgaVnwFHp8dYl0naL5XvDczKv+GVtmcBo/MXSOlXDidL01IyGfhCeqXYzMxaqGlBJSIWAe8FvkSWDHK6pMOpbmb8bukV4/8Ffh0Rd+eu+wzwCPDPXX23pImSOiR1vLp6VV1+j5mZddbUV4oj4h3gbuBuSS+SjY1cAewnaZPSozFJmwDvA0qDGwtKmYu78E3gF8D/dPG9U8kyGDNm4KBKAczMzOqgaT0VSftLGpa2NwH2Bf4cEU8DjwIX5KpfAMxMx3oUEU8CjwNH1bfVZmZWi2b2VHYArpG0Rdp/BLgybX8a+L6kp8keez2UympxEVlw6tbAMaPZq6OjxkubmVk1lGU76T/a29ujw0HFzKwmkmZERHtP9Tyj3szM6sZBxczM6qaZA/XblrIPS3pB0vO5/chtz5J0fjrnfknP5hbpQtLtudQvFVO4NOs3mZnZ+po2UB8RS4GxAJImAcsj4tK0v7ybV4ZfB/4eeEDSULKZ+XmVUrjc2lU75i2dxz437NNtW+dMmNPtcTMzq6w3/Kv+ZmB82v4YXQSMCilczMysyYoSVAaVPf46KXdsOvChlIZlPNnCXZ10kcLFzMyaqCiLdL3dzeOv1cADZFmOB0XEwtwQC6xL4RLAHfkULiWSJgITATbbdrO6NtzMzNYpSlDpyc3AbcCkCsd6SuGyXpqWQbs6TYuZWaMU5fFXT34HfAuY1uqGmJlZ14rSUxmUHmGV/CYizi/tRDbt/9J6fNHobUfTMcEz6s3MGqElQSUiJpXtV1wLJSI+3EX54PTnQmBMfVtnZmYbqrc8/jIzs17AQcXMzOqmcEGllIIlt3+apCvT9qRcepcnJU1xWhYzs+IoykB9LS6PiEtTMPkf4BDgvnwFSZumGfadLX4UJg3p+VsmvbHxLTUz62d687/yNwcGAq/B2uST35T0W+DslrbMzKyfKmJPpfz14m2AO3P750g6BdgFuDsi8nWHRsQhzWikmZl1VsSeytsRMbb0Ab5WdvzyVL4DsKWk8bljXeUFmyipQ1LHyys8od7MrFGKGFSqEhF/BX4DfChX/FYXdadGRHtEtG/fpkpVzMysDnptUEkLd/0dsKDVbTEzs0wRx1R6UhpT2QyYDfygprOH7QeTnKbFzKwRlKXV6j/a29ujo8NBxcysFpJmRER7T/V67eMvMzMrHgcVMzOrGwcVMzOrm6oG6lM+rpeBcRHxVK78u8DiiLg47V8BnACMjIg1qew0oD0izqriew4ArgcGAXcBZ0dESNqGbA7KKGAh8ImIeE3SJ4Hz0unLgTMi4rHuvmPO828w6vxf9/ibF07+aI91zMxsfbX0VG4G1k40TLm3TiBNOEz7xwPPsf7ckVpMIVtLfo/0GZfKzwemR8QewPS0D/AMcEhE7At8g7RksJmZtUYtQWUauaBCFjgWRsSf0/6hwFyywHByrQ2RtDOwdUQ8lFZ6vBE4Lh0+Frghbd9QKo+IByPitVT+e2BErd9rZmb1U3VQiYjZwBpJ70tF41l/zfiT0/5twFGSNquxLcOBRbn9RakMYMeIWJLasYQsRUu5TwN3V7pwPk3L6hXOPmxm1ii1DtRPA8ZL2pSs9/BzAEmbA/8E3B4RbwIPA0fWeO1K+VOqmkQj6VCyoHJepeP5NC0D2qpIe29mZhuk1hn104B7gN8CsyPipVQ+DhgCzMmyp9AGrAB6HhFfZxHrP74aASxO2y9K2jkilqTHZKXvRdK+wI+Af4yIpT19yT7Dh9DhQXgzs4aoqacSEQuApcBkOj/6+peIGBURo4BdgSMltdVw7SXAMkkfSHm9TgXuSIfvBCak7QmlcknvBm4FPhURf6zlt5iZWf31GFTSo653ckXTgD3Jxk5IgeMj5HolEfEW8ABwdCo6TdKi3KerAfUzyHodT5MliiyNkUwGjpA0Hzgi7UOWFn9b4AdpiWHnXzEza6Eec3+lgflrIuLA5jSpsZz7y8ysdnXJ/SXpM2Q9kwvq1TAzM+u7ug0qEXF1ROwdEffU6wslLZf0cHpcNUvSs5JelrRPOj5R0pPp84ikg3Pn3i/pqXTeE5Im1qtdZma28VqynkpEHFTazqVxmSPpKODfgIMj4hVJ+wO3SzowIl5Ip3wyIjpS6pYFkq6PiL/kry9pQESsrvTd1aZpAadqMTOrVdESSp4HnBsRrwBExEyyGfRnVqg7mGz54NWwtgd0oaSHgb9tUnvNzCynFT2VQZJm5fa3IXtlGGA0MKOsfgfrXicGuEnSO2S5wf491yPZEpgbEV9rQJvNzKwKrQgqb0fE2NJO6fFXN/XF+jPrS4+/tgcelPSblH9sNfDLihfIxl4mAgzYevuNbL6ZmXWlaI+/HgcOKCvbP5WvJyJeBmYCpfGZlV2NozhNi5lZc7RkoL4bFwPfljQuIpZKGgucxrrAsVaadLlfOqdqTtNiZtY4hQoqEXGnpOFkj7UCWAacUspQnNwk6W1gC+D6iCgfgzEzsxbpcUZ9X+MZ9WZmtavLjHozM7NaOKiYmVndNC2oSNpJ0s2SFkh6XNJdkv5G0tyyepMkfTFtX5LStcyWdJukobl6B6c0LqWULk7ZYmbWYk0ZqE/ro9wG3BAR41PZWGDHHk79b+BLEbFK0reBLwHnSdoJ+ClwXETMlLQd8F+Sno+IbnOwPLZsBTvdN6u7Kt164dCxPVcyM+unmtVTORT4a0RcXSqIiFnAc92dFBH3RMSqtPt71q0MeSbZm18zU71XgP8HnF/vhpuZWfWa9UrxGDqnXynZrSxty07ApRXqnQ7ckrZHk+UEy+tI5WZm1iJFmKeyoCxty6TyCpK+AqwCbioVsX7qlpKK70fn07RssuPOG9lcMzPrSrMef82jc/qVqkiaABxFlvOrFDTm0Tlf2AFUSOcC66dp2WTI0EpVzMysDprVU7kX+Kakf42IawAkvR9o6+4kSePI0uEfEhErcoeuAh6WdGtEzJK0LfBt4MKeGvK+rdro8GC7mVlDNKWnknoYxwNHpFeK5wGTgMU9nHolsBXw32m1x6vT9ZYApwDXSHoSeBC4NiL+s1G/wczMeuY0LWZm1iOnaTEzs6ZzUDEzs7opfFCRtDqNp8yV9J+lVC2SRnWX4sXMzJqvCPNUerJ2+WFJN5DNpr9oQy+2bNkcpt+720Y16PDDFmzU+WZmfVXheyplHgKGt7oRZmZWWW/oqQAgaQBwOPDjXHFVKV7yM+p32KHX/GQzs16nN/RUBqXAsRTYhixzccmCiBhb+gBXV7pAfkb90KG94SebmfVOveFv2NKYyi7A5mRjKmZmVkC95llQRLwh6XPAHZKmbOh1ttpqHw4/zJMfzcwaoTf0VNaKiEeBx4DxrW6LmZl1VvieSkQMLts/Orc7puzYpGa0yczMKutVPRUzMys2BxUzM6ubhgYVSTtJujmlu39c0l2S/qa79CqSrpf0TErN8qSk/+jm+t+QNDvVvUfSsEb+HjMz617DxlQkCbgNuCEixqeyscCOVZx+bkT8QtJA4HFJN0bEM2XXHwBcEhFfTfufA74GfKa7Cy9evJhJkybV/HvK1eMaZmZ9TSN7KocCf42ItRMSI2IW8FwN1xiY/nwLQNJCSV+T9ABwYkS8mau7JV2sUW9mZs3RyLe/xgAzujjWU3qVSyRdAOwOfC8iXsodWxkRB5d2JF0EnAq8QRbIOsmnaRkyZEitv8PMzKrUqoH6ntKrnJvKdwIOl/R3uWO35CtGxFciYiRwE3BWpS/Lp2lpa2ur488wM7O8RgaVecABG3OBiFgO3A8cnCt+q4vqPwU+vjHfZ2ZmG6eRj7/uBb4p6V8j4hoASe8Hqu4qSNoUOAj4fhfH94iI+Wn3GODJnq45bNgwD7KbmTVIw3oqERHA8cAR6ZXiecAkYHEVp1+SxlxmA3OAW7uoNzmtCDkbOBI4e+NbbmZmG0rZ3/39R3t7e3R0OKGkmVktJM2IiPae6nlGvZmZ1U1Lg4qkEZLukDRf0p8kXSlpC0kflvSrCvXvl9SetodIujE9Wlsg6SZJ72r+rzAzs5KWZSlOM+5vBaZExLFphvxU4GKymfg9+TEwNyJOTdf7OnA9cGx3J/3l+eUsOv93G9P0tUZM/mBdrmNm1le0MvX9YWQTGa8DiIjVks4B/sz6SwZ3Iml3steVT8oVXwgskPTeiHiqQW02M7NutPLx12jKZtyntCsLyWbSd2dvYFZErM6duxp4FNirvs00M7NqtTKoiMq5ulTvcyVNlNQhqePVFa/X0EQzM6tFK4PKPGC919MkbU2Wxbinx1fzgP0krW1/2t4XmFleOZ+mZZu2oRvdcDMzq6yVQWU60CapNNA+ALgMuBJ4u7sTI+JpskddF+SKLwCmR8SzjWmumZn1pGUD9RERko4HrpL0VWB74JaIuEjSh8kSSS7KnXJi2SVOB74v6WlgCPAH4Gh6sPnwwX5ry8ysQVr59hcR8RxZzi5SJuJpkg6IiPuBQRVO+XDu3NeBT6Vz3wvcBXwk/WlmZi3Q0qCSFxEPArts4LlPAbvVt0VmZlYrp2kxM7O6KWRQkbRa0ixJj0maWbZIV77eSEn3SXpC0jxJzlJsZtZChcxSLGl5RAxO2x8BvhwRh5TVGQDsAOwcETMlbUU2mfK4iHi8q2uP3GZo/PsRB3d1uGZfuKVTijIzsz6nL2Up3hp4DSAlmrxP0k+BORGxJCJmAkTEMuAJYHjrmmpm1r8VZqC+zKC0SNdAYGeyPGElBwJjIuKZ/AmSRgH7AQ83qY1mZlamqEHl7YgYCyDpb4EbJY1Jxx6pEFAGA78E/j3lD6Ps+ERgIsC72iq9qWxmZvVQ+MdfEfEQsB3Z5EiAt/LHJW1GFlBuioiKyw7n07RsucXmDW2vmVl/VvigImlPYACwtMIxka2r8kREfKfZbTMzs/UV9fFXaUwFsszDE9J6K+X1/p5sVv2cXP0vR0SXs+p3fM/ufmPLzKxBChlUImJAF+X3A/fn9h+gulT5ZmbWBIV//GVmZr2Hg4qZmdVNrwgqkkZIukPSfEkLJF0haXNJB6Z0LqWULse3uq1mZv1ZIdO05KU3vB4GpkTEdSk9y1TgVeA/gL9ExCpJOwOPAcMiYlVX13v39u+N8z4+pa5tPPPqw3quZGbWi/WlNC2HASsj4jqAiFgNnEO2SBe5ADKQyuvWm5lZkxTy7a8yo8kSRa4VEW9KehbYXdIg4FqytVg+1V0vxczMGqs39FRE5R6IyFYlfjgiRgPvB74kaWCnitJESR2SOpavfL3BzTUz6796Q1CZB6z3HE/S1sBIYEGpLCKeIEvhMoYy+TQtgwcObXBzzcz6r94QVKYDbZJOhbXrqFwGXA/sKGnTVL4L8F5gYWuaaWZmhR9TiYhIrwr/QNJXyQLhXcCXgU8A50v6K7AG+GxEvNLd9XbYZSu/rWVm1iCFDyoAEfEccHSFQz9JHzMzK4De8PjLzMx6CQcVMzOrm0IEFUkh6bLc/hclTcrtT5T0ZPo8Iung3LEfpxQtsyX9Iq0CaWZmLVCUMZV3gI9J+lb5QLuko4B/Aw6OiFck7Q/cLunAiHgBOKe0hLCk7wBnAZO7+qKVc+fxxJ571f0H7PXkE3W/pplZb1OIngqwiiyf1zkVjp0HnFsKNhExE7gBODPtlwKKgEE4VYuZWcsUJagAXAV8UtKQsvJOaVqAjlQOgKTrgBeAPYHvN7KRZmbWtcIEldTjuBH4XBXV10vdEhH/BxgGPAGc1KlyLk3Lq6udGszMrFEKE1SS7wKfBrbMlT0OHFBWb/9UvlbKXnwL8PHyi+bTtGwzoCjDSGZmfU+h/oaNiFcl/YwssFybii8Gvi1pXEQslTQWOA04KI2j7BYRT6fto4Enu/uOgWNGs1dHR+N+hJlZP1aooJJcRvYGFwARcaek4cCDkgJYBpwSEUskbQLckBJMimyRrjNa0WgzMytIUImIwbntF4G2suNTgE7LNUbEGuDvG95AMzOrStHGVMzMrBdzUDEzs7pxUDEzs7qpakxF0nLgZWBcRDyVK/8usDgiLk77VwAnACPTeAeSTgPaI+KsThfu/D0HkC2+NYhszZSz03oqt5AtwAUwFHg9IsZK2hb4BdlSwtdX8x3zls5jnxv2qeZn18WcCXOa9l1mZq1WS0/lZmB8aSe9eXUC2dyQ0v7xwHPAhzawPVOAicAe6TMOICJOioixETEW+CVwa6q/Evgq8MUN/D4zM6ujWoLKNHJBhSxwLIyIP6f9Q4G5ZIHh5FobImlnYOuIeCgigmx2/XFldUS22uM0gIh4KyIeIAsuZmbWYlUHlYiYDayR9L5UNJ70l3tyctq/DThK0mY1tmU4sCi3vyiV5X0QeDEi5tdy4XyaltXLVtfYLDMzq1atA/XTgPGSNgWOBX4OIGlz4J+A21MOr4eBI2u8tiqUlWccLgWumuTTtAzYakCtp5uZWZVqnfw4DbgH+C0wOyJeSuXjgCHAnOwJFW3ACuDXNVx7ETAitz8CWFzaSYHsY3TOA1aT0duOpmOC07SYmTVCTT2ViFgALCVbBKv80de/RMSoiBgF7AocKamt81W6vPYSYJmkD6Sxk1OBO3JV/gF4MiIWVbyAmZm1XI9BJfUQ3skVTSNbt+S2dLwN+Ai5XklEvAU8QJbgEeA0SYtyn3yPJO8M4EfA08AC4O7csfIxnFL7FgLfyX3H3j39JjMzawxlL1p1UyEbmL8mIg5sTpMaq729PTqcpdjMrCaSZkREe0/1uu2pSPoMWe/ggno1zMzM+q5ug0pEXB0Re0fEPRvzJWlGfn5/gaSXJc2S9IKkNZIOqVRf0ghJd0ian867Ir1thqQjJM2QNCf9edjGtNPMzDZOq1Lff4OUukXSJOB0sleSf5uvlAbsbwWmRMSxkgYAU4GLgHOBV4CjI2KxpDHAf9F5bsv6Fj8Kk4bU+ecU0KQ3Wt0CM+uHipJQ8lrgJEnblJUfBqyMiOtg7ZLB5wCnS2qLiEcjovTa8TxgoKQtmtZqMzNbT7OCyqD0qGuWpFnAhWXHl5MFlrPLykcDM/IFaXLls8DuZXU/DjwaEe+Ula83o/7lFd2/mGBmZhuuWY+/3k7JIIF1mYvL6nwPmCXpslyZ6DyrvlO5pNHAt+liFn9ETCV7bEb7sAGOKmZmDVKUx19ExOvAT4HP5ornURZ80nr0I8nmsZDmvNwGnJomZ5qZWYsUYo36nO8Af2Bdu6YDkyWdGhE3poH6y8jWTlkhaSjZpMsvRcT/VvUNw/aDSZ6nYmbWCIXpqQBExCtkvY4t0n6QrdFyoqT5wB/J0tx/OZ1yFtnYyldzYzY7NL/lZmYGVcyo72s8o97MrHZ1mVFvZmZWCwcVMzOrm0IHFUmRf8VY0hfTDPzS/kRJT6bPI5IObklDzcwMKN7bX+XeAT4m6VtpEH8tSUcB/wYcHBGvSNofuF3SgRHxQlcXnPP8G4xoBIKgAAAEIUlEQVQ6v5a1w8zMer+Fkz/alO8pdE8FWEU2afGcCsfOA84tBZuImAncAJzZvOaZmVle0YMKwFXAJyWVZ4HslMIF6Ejl68mnaVm9wokWzcwapfBBJeX6uhH4XBXVK6Z1iYipEdEeEe0D2vpBhmIzsxYpfFBJvgt8GtgyV/Y4cEBZvf1TuZmZtUDRB+oBiIhXJf2MLLBcm4ovBr4taVxELJU0FjgNOKi7a+0zfAgdTRqwMjPrb3pFUEkuI0vLAkBE3ClpOPCgpACWAadExJJWNdDMrL8rdFCJiMG57ReBtrLjU4ApzW6XmZlV1u9yf0laBjzV6nYU2HZkyzRbZ743XfO96VpfuTe7RMT2PVUqdE+lQZ6qJilafyWpw/enMt+brvnedK2/3Zve8vaXmZn1Ag4qZmZWN/0xqExtdQMKzvena743XfO96Vq/ujf9bqDezMwapz/2VMzMrEH6bFCRNE7SU5KelnR+heNbSLolHX9Y0qjmt7I1qrg3H5I0U9IqSSe0oo2tUsW9+bykxyXNljRd0i6taGerVHF/PiNpjqRZkh6QtHcr2tkKPd2bXL0T0lpRffONsIjocx9gALAAeA+wOfAYsHdZnc8CV6ft8cAtrW53ge7NKGBfskSeJ7S6zQW7N4cCbWn7jP7y300N92fr3PYxwG9a3e6i3JtUbyvgf4DfA+2tbncjPn21p3Ig8HRE/Cki/gLcDBxbVudYsvVXAH4BHC5JTWxjq/R4byJiYUTMBta0ooEtVM29uS8iVqTd3wMjmtzGVqrm/ryZ292SClnD+6hq/s4B+AZZ3sKVzWxcM/XVoDIceC63vyiVVawTEauAN4Btm9K61qrm3vRXtd6bTwN3N7RFxVLV/ZF0pqQFZH95VrNkRV/Q472RtB8wMiJ+1cyGNVtfDSqVehzl/2Kqpk5f1F9/dzWqvjeSTgHagUsa2qJiqer+RMRVEbEb2eqsFzS8VcXQ7b2RtAlwOfCFprWoRfpqUFkEjMztjwAWd1VH0qbAEODVprSutaq5N/1VVfdG0j8AXwGOiYh3mtS2Iqj1v52bgeMa2qLi6OnebAWMAe6XtBD4AHBnXxys76tB5Q/AHpJ2lbQ52UD8nWV17gQmpO0TgHsjjaT1cdXcm/6qx3uTHmH8kCygvNSCNrZSNfdnj9zuR4H5TWxfK3V7byLijYjYLiJGRcQosvG4YyKiozXNbZw+GVTSGMlZwH8BTwA/i4h5ki6UdEyq9mNgW0lPA58HunwFsC+p5t5Ier+kRcCJwA8lzWtdi5unyv9uLgEGAz9Pr832m4Bc5f05S9I8SbPI/r+a0MXl+pQq702/4Bn1ZmZWN32yp2JmZq3hoGJmZnXjoGJmZnXjoGJmZnXjoGJmZnXjoGJmZnXjoGJmZnXjoGJmZnXz/wGAmaUOOYAikwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# feature importance diagnostics\n",
    "feat_importances = pd.Series(model_rf.feature_importances_, index=df_in.columns)\n",
    "feat_importances.nlargest(20).plot(kind='barh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done  20 out of  20 | elapsed:    0.0s finished\n",
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done  20 out of  20 | elapsed:    0.0s finished\n"
     ]
    }
   ],
   "source": [
    "Y_pred_train = model_rf.predict(X_train)\n",
    "Y_pred_test = model_rf.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.979036008792889, 0.8272216700069426)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r2_score(Y_pred_train, Y_train), r2_score(Y_pred_test, Y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1c180c71d0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl0XOWZ5/HvU6XSLnnRZuPdgLEcYxyk2IFsQOQQ21kAH3oC6Q5Mp+Oc6YTpTkKGZNI5k9M5Wchkme7MTE+cDk16EshmEphgIBINDYmNQXIc20jCOHZsyzZleZVkrVX1zh9VJZfl0l6lWvT7nKNTpdru64v46dVzn/tec84hIiLZw5PqAYiISGIp2EVEsoyCXUQkyyjYRUSyjIJdRCTLKNhFRLKMgl1EJMso2EVEsoyCXUQky+SkYqPl5eVu8eLFqdi0iEjGampqOuWcqxjtdSkJ9sWLF9PY2JiKTYuIZCwzOzyW16kUIyKSZRTsIiJZRsEuIpJlFOwiIllGwS4ikmUU7CIiWUbBLiKSZRTsIpL2Gpr9PLB1Dw3N/lQPJSMo2EUk7dW3+Cn0ealvUbCPhYJdRNLeuuoqugeCrKuuSvVQMkJKlhQQERmPuhVV1K1QqI+VZuwiIllGwS4ikmUU7CIiWWbMwW5mD5nZSTPbF/PYl8zsmJntjnxtSM4wRURkrMYzY38YeG+cx7/jnFsd+dqWmGGJiMhEjTnYnXMvAGeSOBYREUmARNTYP2lmeyKlmlnDvcjMNptZo5k1tre3J2CzIiISz2SD/Z+AK4HVwAngW8O90Dm3xTlX65yrragY9ZJ9IiIyQZMKduec3zkXdM6FgO8DaxIzLBGZbrQeTOJMKtjNbG7Mt7cD+4Z7rYhkp0QFstaDSZzxtDs+CuwArjGzNjP7KPANM9trZnuAm4FPJWmcIpKmEhXIWg8mcca8Voxz7q44D/8ggWMRkQy0rrqK+hb/pANZ68EkjhYBE5FJUSCnHy0pICKX0YHMzKZgF5HL6EBmZlOwi8hldCAzs6nGLiKXUd08s2nGLiIja90GT9wXvpWMoGAXkZHtfwp8heFbyQgKdhEZ2bL1MNAdvpWMoBq7iIxs+Ybwl2QMzdhFRLKMgl1EJMso2EVEsoyCXUQkyyjYRUSmwhSeD6BgFxGZClN4PoCCXURkKkzh+QDqYxeR8WvdFp55LluvHvexmsLzATRjF8kwabFWupYZSGsKdpEMEBvmabFWupYZSGsqxYhkgNgwT9Q1RidFywykNQW7SAaIDXOtlS6jUbCLTJVJHHBUmMt4qMYuMlWm6IBjWhxclZRSsItMlSk64JgWB1clpVSKEZkqCT7gGO2Qidbdo9Li4KqklIJdZAoMF8KTETszj/1M1eNFpRiRKZCM8si66iq6B4KamctlNGMXmQLJKI9oZi7DUbCLTIFEhnAyyjqSXVSKEUlDI7UsqutFRjPmYDezh8zspJnti3lstpnVm9nrkdtZyRmmyPQyUnirti6jGc+M/WHgvUMe+xzwrHPuauDZyPciMkkjhXfdiioe3LRKZRgZ1phr7M65F8xs8ZCHPwjcFLn/Q+B54IEEjEtkWtOBUZmMydbYq5xzJwAit5WTH5KIiEzGlB08NbPNZtZoZo3t7e1TtVkRkWlnssHuN7O5AJHbk8O90Dm3xTlX65yrraiomORmRURkOJMN9ieAeyL37wEen+TniYjIJI2n3fFRYAdwjZm1mdlHga8D68zsdWBd5HsREUmh8XTF3DXMU+9O0FhERCQBdOapiEiWUbCLTDFd4UiSTcEu6aF1GzxxX/g2y2mtF0k2Bbukhym6Hmg60FovkmxatlfSw7L14VBP8vVA04GWC5BkU7BLekjw9UBFpjOVYkREsoyCXUQkyyjYRUSyjIJdJN1No1ZQSQwFu0i6m0atoJIY6ooRGaOGZj/1Lf7B/vPo/aS3Lk6jVlBJDAW7yBgNPWM0ej/pwa5WUBknlWJExqiiOJfn97dTUZyrs0clrWnGLjJG7V393LSsgvaufp09KmlNM3aRMdIsXTKFZuySkWIPZE7VzFmzdMkUmrFLRpqqpW+1drpkIgW7ZKSpKoto7XTJRCrFSEaairJIQ7OfY2e7Abj3xiVJ3ZZIIinYZVob6aSj+hY/V1eW0D0QVG1dMopKMZK2pqK+HVtqGVp2UReMZCrN2CVtxQZtQmbMrdsunpofOZNzXXUVf2h4hHf2vsTZBXXsZO1gkKsLRjKVZuySthI+Y46zmFbdiire7hopKCxh3snneXDTKoW5ZDzN2CVtJXTG3LoNzh0N31+z+ZKn8lZspK/5SfJWbAw/0PBlaPkVVN8GdV9MzPZFppCCXaaH/U9BxTUw0H3Zglqr190N6+6++EDLryAnP3yrYJcMpFKMZL2GZj8Pn66m7eRpdhfcMPoB2erbINAbvhXJQOacm/KN1tbWusbGxinfrmSQOAc6x2Jo++LD2w/xur+LqpI8OvuDzJuZP9jC+OCmVckavUhSmFmTc652tNepFCPpKfZA5ziC/eHth3j1WAeP//4YS8qLCIQcHoODp7tZu2QWfYGQWhgl6ynYJT1N8KpBp7v6Odc7gM8DZy70c3VVMVWleaxeMJP2rv4pXTRMJFUU7JKeJnjVoLLiXGaez6FnIMSdtfO5/9blSRicSHpLSLCb2Z+ATiAIBMZSAxIZr3hL9Q597N4bl0z5cr4i6SaRXTE3O+dWK9QlWeKttFjf4ue6C9spqf90+ICriKjdUTJHvDNRK4pzKT3SQE5eEXuefZS//eluWo6d1zK7Mq0lqsbugN+YmQO+55zbkqDPFRkU70zU3UfPcd63lls6GqkP1jCzIIcFp1/gvhmH2F1/C492rFRZRqadRAX725xzx82sEqg3s1bn3AuxLzCzzcBmgIULFyZoszIdffOZVp7c+wbVc4rZdeQc2wdW8NyM67nt+nn49r7B/TOeY35XJxd2H6Lwmn9M3CJiIhkiIaUY59zxyO1J4JfAmjiv2eKcq3XO1VZUVCRiszIN7a5/hIW//RzLzr3IM6/6CQSDeD3GQNBx/63Lee7+m1hSXgwOyovy1LMu09KkZ+xmVgR4nHOdkfvvAf5+0iOTaWW4jpeHtx8CLl7BqHfnY/RaHje5Rn7jricUhBwv1CyaefHD1myG/U9Rtmw9Dy7X2aUy/SSiFFMF/NLMop/3iHPu6QR8rkwjl6293rqNkvofM7fzWv5f33W87u9idlEuK6llDS/xdKgGB+QYlBfl0vJGFw3NkfdOsAdeJFtMOtidcweB6xIwFpnG1lVXXbLGC/ufYm7ZbN56didbA9fS3tnHme5+zhXV8kJeLWcC/RiOkIOQg6Vlhaqli0TozFNJvdZt1B14iro3rYflqwbXTi/s6uM3gRsJRdapK/AYIQc+rxFyDp8Xrqos4dPrrrn0l4LINKdgl9SLWfCrIVRDSf2Pycmr5Gj7KZ4JXD/4sv5giNlFuQRCjsrSPPJ9OXx63TW6hJ3IEDpBSVJv2frwBTCWrae+xc8fZ72Do/5T/Lpv9SUvK873UVacS1VpHldWFPOFDdUKdJE4NGOXlGsI1VA/MJ/Oxn5+d+AET9sSevs/Sl8o/LzHYE5JPgHnuPfGJcOH+QTXcBfJNgp2Sbm2nVt594nn+GnXtXS5GszgXTRRl9NEQ6iG10rfxtuurhj9DNIJruEukm0U7JI8o8ygo73rH+7eQVOPhzpPE88GaijyebnXPUMlZ1mQc4reD3x8bCWXCa7hLpJtFOySPKPMoOtb/Jzq7OUnHdfyZtvO08EaPEBlaT7zPYXk93ZRPnMGZWOto6t/XQTQwVNJpuhB0cJKeOK+y5bVrSjO5XcHTvN497X8XXAzz4ZqMI8xuzCHx3I20JqznHOzV8V9r4gMT8EuybN8A3zgu9B98uLMPcbuo+cIOceFgRD9QYcHKM718pr/Ak8HruerOX/NKf/RuO8VkeEp2CX5YtoZgfDs+4n7mOd/nv5g+OwjB2Dhp6+pKuJ8zwA+r1FetQAONIRn/SIyJqqxS/It3xBuaXzVz4F/+y13Hv8e3eRxPUf4KdcC4Uyfke/jPSvnALBq/iy6B4Jc6bsAV9WFZ/0iMiaascuUeHj7IV547SRNR85TH6qhgD78bgZfy9nCuz1NFPo8OBzHznZTUZx7cbndobN9ERmVZuySdA3NfgoO/ob/7Gmi3TuD6zx/xIB5dooDbh63endxcMY7mTczn6srS2jv6ufBTdHldtXpIjJeCnZJmGhf+l2l+1jdswOWrachVMNjj27h73N+QCE9FDLAeQo5xQwCzsu7PKd53bOIfyj4AaFF6wcvZSciE6dgl4SJrql+9veP82JBMTOOPspPjjXyOe8j5BCgiH4GzEuB6+ekmwVAd8lV1AR20Rqazzt6drB6090p/leIZD4FuyRMdE31/GAXK0//OyeCM/lWTj0BjFJ66XR59JBPU2gZXRRyys1gpevlD6XvYnlxr+roIgmiYJeEiS6f2/GVVjyhAa6yEzjCHS+nKKGDYr4WuJu/y32EPzGXmbk+gu/7X1qhUSTBFOySENH6+gdPP0Rt3xlyCGEWDnUHFLk+drpqvp73LxRYgFmew8y847OwXKEukmgKdkmIb/3mNeadfJ7rcn6Ex8KXrHPO8OAI4OGgm0PAV0xRXi6uv5/A7GvU7SKSJAp2SYj3n3mID+c8jY8gHhxBjKdDb6HWs58AHjpsBgdnvZPCmeEzTYuuqAmvAaO100USTsEu49fwZbp2b8UfLKIi8Aa9ff18nA4sWncBPDje42mijTLAy0vlm7iu7m5Y8anwC564T2uniySJzjyVcevavZW2LsfCnmby+09Tbh2D9XQit+H7jnmcwe8pp6wk79KDpDqjVCRpNGOXcWlo9nO8t5a3h35Lr/NSZAODz7nIjD0UWczLZ44e89JRdCV1nl2XfpDWThdJGs3YJa6GZj8PbN1DQ7P/ksfrW/zk5hjlnA6HemR6HtsB87vQSno8eXh8BRTNmMN7ryll/to7UvHPEJmWFOwSV/Qs0vqWi8He0Oznd6+3c1v3Lyi2/sEwN8A5CDm4QC5bFn2L3soamLkQypaG12TX7FxkyijYJa511VUXV1gE/vpHjfQ88uc8130HuYQuCfOQgy6XSwAvjTm1/Oiv3kr5LffBgjWwZnNq/yEi05Bq7BJX9EBn286tdDzzQ759dj+5nvBJR4MMel0u+918ltoJjuUt5pbSk+ELaaiGLpIymrHLYD19d/0jl1xftL7Fz6pjP6X4XCt5NiTUARzsd/OZY2c4lTOXpfmdMPsqXcZOJMU0Y5dLVmV82vK5Yt9/54qZ3+UT5/zM7f/j5YFOuAzT43IooYe+3Nkcn3E9V844DSWVw7cwtm4Lh75OShJJKgV7momuubKuumpii2ONMzwbmv0cO9vNm7t3UNJ3gjnB88zmPDknA8xznXiGhLpzDAb9BU8hR2bcwPzgEZaUWriePnSbsePZ/5ROShKZAirFpJl43SjjEhueMUZqX7y6soT39/2aRaFjLKWNHBegJNSJucs/3iwc7s6gzDq55apilr3/M8y/55/jh3XseHRSksiUSEiwm9l7zew1MztgZp9LxGdOV0O7UcYtGp6FlZfVy+P9wohub3ZRLjPceXwuwCzXGW5jjDNbDzkIRtbiPUfJxZl4ZDvDjif6F4RaH0WSzpyLMy0bzweYeYH9wDqgDXgFuMs51zzce2pra11jY+OktiujiK7FMtANH/hu3BJPQ7Ofh7cfAuD+RQdZ+eJ/whunng6RmrrHR745Ds54K7nnD9G07FPcUbTnku2ISPKYWZNzrna01yWixr4GOOCcOxjZ8E+ADwLDBrtMgehMOlL2iF4EI1bbzq3cebyBE8EZ+I40jvjnWwjIN4d31hKuXrwUln2CRcs3XFpDF5G0kIhgnwccjfm+DVibgM+ViYgN2lFm0O/v+TW+YBPFrjdu50uUcxDwFpC/5h5of+3SJ9WvLpJ2EhHs8SLhsvqOmW0GNgMsXLgwAZuVuKIHK1/ecnl3TDT0CyvheBNlp14hRP+IoQ7Q78nhdO4VzGt/DdpbwJunzhaRNJaIYG8DFsR8Px84PvRFzrktwBYI19gTsN3pYzwtjNESDECnH+q/ePG5+i9yLphP/rnX8dKHN3T5AdKo6H+gfpfD/rxVhGo/xryeHeDNhTMHoOaeSf+zRCQ5EhHsrwBXm9kS4BjwIeDuBHyuRI2n/ztaGmn4MrzyPY5757PrFw+xcFYBq3JLKD7xB8yF4na9QLjsgkHQGV2eIlrm3cn/9nwIDsP9i2A1O6DmXs3WRdLYpIPdORcws08CzwBe4CHn3KuTHplcNORA6Jh0n4SFN9L/+l4e799AqBsezPstRSEjf4TSS8jgghWyJ7SUr878GqvKZuI/fBaARztWsnqTfmeLpLuEnHnqnNsGDNPILJM2kQOUPR1wZDtQygPeRwgBs4JnwC49ezTKOegilz7LJ4iX133XUFacy7rqKo6d7QaYeG+9iEwpLSmQwaK96XeV7mN1z46LNfjWbXDwWQgMsJBj4dPQYsM8JtSjpzEccHP5gy3n3bmtFAx08J7cP7DwxiXUeZqoq4zW99861f9EEZkALSmQCVq3XXIWaVT0bNK+5icv1uBbt0H9F+lxeYSCvZhjxHp6v4OHQ7fi88LSKyopGjhDri+H+b6ucN/7y1vg6Mvw3FfijkFE0o+CPVWGCeu4hln/JbocQN6KjRdP29//FCd9C3D9HeE1Xbg4K49yDkIh+FjgM3wi+BlW5b1BQUExJZ1/5My8m8ktrYDr/jzmDUD36bhjEJH0o1JMqgzX6RKvtXEcB08fu7CKqhO/p9r5yKP/suedg/MU8rRbw668tbxp3gxqK9t58XA3hfTxg9n/hQc3rbr4hjWbL/a+d5/UGaYiGUDBnirDhfXQwI8X9K3b4OUtLPB3cl3VnfS1vciRstn88Ykfcrwznw3WQo4FwmWYSAkmOmnvwUsuAW60fXzK+3Nu8QGFi6j2bufMhX7uKt0HxAS7ziwVyTgK9lQZLjCHBn68mf3+p6DjBFd4Apw7+yLBggqCB56mK7CAv/DsxWsBPMSEuguv9dJDHgPkYAZl3h7uyG+kuPJW6D5J+fyrKPcVQs8OdBqCSGZTjT3dDF3aNt4a5svWQ+lcSioWkLdiI2dOtvFcYBXVdpSDbm54sfRofR3AoA8f/xJ8L9vdSkpyghSXllG8etPFz9Za6SJZQ8GebOM5SDpEQ7OfB16dT8NVf3exLPPEfeEnP/IrGmq/x6f+cAW/89TyFmuh0k6zyPz0kAsW6WqMBHy+DfBhbz193iKOL/kzuPImmF9zcWNaK10kayjYk22YjpZRtW6jpP7TXHdhe/jiGJE2RjpPDn5WfYufpeWF/KTjWgoYIIAPDw4jcpWjyIFSzDCDfE+IYxXv4sTpM+wuuIGdz/yYIx2o00Uky6jGnkyt2+BcZEXjNZsveWrYa5tGD5aeO8rcsgVw+kUqa2+H/d+H2VfSdeI1Gg92Ur1rMR+llH/23c0tnn568ZHjAgQsh3OuhGNWTsB5ec17FTf7mpnt7SGQN4djc26hs/ouHmvxc92sDjj9IgtrPzyFO0VEkk3Bnkz7n4KKa8K16yFL57adrqaw7F207dwKB1outhOeOxp+D7CwlHDoLq8CT/ig6i/6b+Idh/4Rn/VT6dp5c88O8MAJylhg7TS6ZeQzwPFQOblveh+b7vrY4DZnLFvPg7TB/u9TXnoDj3Jj+JfG8qpLxjbWC2GLSHpSKSaZ4h2QjJRm6jy76B4IUufZFS7VtPwqfAvh96zZHLfmXZrv48ngGvqcj6OunIZQDe1uBmuslVOulCt5g596NtD2jge5466Phd8UWz+PbH91zw4e3LTq0r8WJlo2EpG0ohl7MsVraYy0M85fewcPLl8FrXeEg7T6tvCMfc3m+LPlSOgWHnqG7wT/I98O/ofBp+pymnjZLWcJb/B4xce5ve7uyy6DN3T7cbtfJrKKpIiknUlfzHoidDHr8dtd/wiBl/4P3f1Bfhi8lWdDNbzb08RHvM9QRgdnKOWntp7i6z5w6ZmjIpI1xnoxa5ViMkBDs5+fNbYxe+Akffio8zRhwMbc37M0r5OqIi/BkvkcnP1Ojp3tpqHZn+ohi0gKqRSTRmI7ZSDczlhRnMvPG9v4bN92/ujmcKW9wc+Ct1Dg8/BbzxruWBj+i+vmNfdweMc+3n7mMYJPOnYf+yiPdqy8vOtGRLKegj2N1Lf4OdXZyxd+uZe+QJCCXC+hEHg88Ezgetbl7OJfSzbycvcqXCBE64y3wUc+P/j+up1/hbevneK8HJqbn6RwyXXUt/gV7CLTjII9DURn6hXFuTzfepKO3gECwRB9gRBm4Use/ZurYTtvoe9sCC8BAg4WlxVe8jnz194BL58CIG/uRro7grrqkcg0pGBPA/Utfk519PLyoTOUFefy1oGXucHtpD5Yw/OuhtwcL5UleXT2BSnJ93Kue4CCXA8tb3Rd+kExXTirI18iMv0o2NPAuuoqvrKthdsLdnP1+d9SHDzB68yjztMEIXgPu2jJeRvFNe9n99FzHD59gf6AY+O1c1I9dBFJQwr2FIuWYTZeO4e1ra9w1ldERX8nc4N72BZaQ52nCU9uMbX9O8lf8BHuv3X55R+iM0ZFJIbaHVOoodnPV7a14HltG4u3f57DfUX0dndy3mayy3c9czwd7Cm8AQt0c3DW28OLgcWjM0ZFJIZm7CkQnaUfO9vN0rJCVh98iZCvENd5ku+VfZo397zEV990lN0F76OpYyV/Kr6N9q7++AdCR1hoTESmJwX7JA27SuMIHt5+iIPtF7jQF+CKmQU05b+VtQM7aeB69vs7KVxwI3zg7WM7ABpvoTERmdZUipmk+hY/hT7v8GWSYZzp7qdnIMihUxc4Pe8W/mfx3/A771qK83I4eraXB7buYXf9I6NfpENXPhKRIRTsk7SuuorugbH1i3/zmVZu/ubzlObnkJ/jIRB0eD1QdeJ5/qHgB9xRuId8n5fy4lwKfV76mp8cvXauKx+JyBAqxUxS3YqRSzCxpZon975Bfk64/3zlvBm8eqyD7oEgHyrdy6qlV7Bq4Ahf+sBnB9+Tt2Jj+OLSmo2LyDgo2JMstlSz8do5/LyxjdlFuZTm5zAQclxTVcyuwhtYNdAyGOAXf1msAu5O6fhFJPMo2JOooTnc+QJw741LqFtRRXtXP4U+L8/vb2ft4lkcPN3N/HdvAq3nIiIJomBPgMs6Y2Iuf3d15bvCV0paUXVJ0G+8dg7tXf3ctXaRFukSkYRSsCdAbLmlbkUV7H+KIx2wqP3faQhez703Lhk8GWlpeSHlJfnxzyAVEUmASXXFmNmXzOyYme2OfE3L1ozLOmOWrefE6TPsLbqBY+d6gXD4Ly0r5OCpbiqKc3lg6x5dEENEkiIR7Y7fcc6tjnyN0HCdXRqa/YPhXLeiigc3rWL30bPc/M3n+ebhpXSu+za/7FnN0rLCwTJNeWk+X9hQPVhnH2/vu4jIWKgUM0GXlV9gsJ3x541ttHf1D9bRo7X32Fp67JWSREQSKRHB/kkz+wjQCHzGOXc2AZ+Z9tZVV10WztF2xt6BIKc6wyWYeBeWHq33XURkMsw5N/ILzBqAeAt/fwF4CTgFOODLwFzn3F8O8zmbgc0ACxcurDl8+PAkhp2+Hti6h1MdvRw83c0XNlQrwEUkYcysyTlXO+rrRgv2cWxwMfBr59zK0V5bW1vrGhsbE7LddBDvItS6iLSIJNpYg32yXTFzY769Hdg3mc/LBLEHTaOG1tuj5Rd1vohIKky2K+YbZrbXzPYANwOfSsCY0lq81RzjLQQ20VUfRUQma1IHT51zf5GogWSKeAdN4x0Mjfc6EZGpoHbHcapbURW+yPT+74Nn+GuMqvNFRFJFwT4RL2+BzhPhS9JpHXQRSTO60MZEJaaZSEQk4RTsQ8TrernMms2wcI0uHi0iaUmlmCHiLRVwmeUbVIIRkbSlGfsQsa2LY5q9i4ikGQX7ENETjOpWVKkXXUQykoJ9BPFOPBIRSXeqsY9AvegikokU7DGii3lVFOdeso66iEgmUSkmRrSm/uTeN1RbF5GMpWCPEa2pb7x2jmrrIpKxVIqJoZq6iGQDzdhFRLKMgl1EJMtkbbDrrFERma6yNth11qiITFdZG+w6a1REpqus7YpRh4uITFdZO2MXEZmuFOwiIllGwS4ikmUU7CIiWSajgl296SIio8uoYFdvuojI6DIq2NWbLiIyuozqY1dvuojI6DJqxi4iIqNTsIuIZBkFu4hIllGwi4hkGQW7iEiWUbCLiGQZBbuISJYx59zUb9SsHTg8wbeXA6cSOJxk03iTS+NNLo03ucY73kXOuYrRXpSSYJ8MM2t0ztWmehxjpfEml8abXBpvciVrvCrFiIhkGQW7iEiWycRg35LqAYyTxptcGm9yabzJlZTxZlyNXURERpaJM3YRERlBRgS7md1pZq+aWcjMaoc893kzO2Bmr5nZraka43DM7EtmdszMdke+NqR6TPGY2Xsj+/CAmX0u1eMZjZn9ycz2RvZpY6rHM5SZPWRmJ81sX8xjs82s3sxej9zOSuUYYw0z3rT92TWzBWb2nJm1RLLhbyKPp+U+HmG8SdnHGVGKMbNqIAR8D7jfOdcYeXwF8CiwBrgCaACWOeeCqRrrUGb2JaDLOffNVI9lOGbmBfYD64A24BXgLudcc0oHNgIz+xNQ65xLy55lM3sn0AX8q3NuZeSxbwBnnHNfj/zynOWceyCV44waZrxfIk1/ds1sLjDXObfLzEqAJuA24F7ScB+PMN4/Iwn7OCNm7M65Fufca3Ge+iDwE+dcn3PuEHCAcMjL+KwBDjjnDjrn+oGfEN63MkHOuReAM0Me/iDww8j9HxL+HzstDDPetOWcO+Gc2xW53wm0APNI0308wniTIiOCfQTzgKMx37eRxJ01CZ80sz2RP3fT4k/DITJlP8ZywG/MrMnMNqd6MGNU5Zw7AeH/0YHKFI9nLNL9ZxczWwy8GdhJBuzjIeOFJOzjtAl2M2sws31xvkaaOVqcx6a8tjTyohf2AAABwElEQVTK2P8JuBJYDZwAvjXV4xuDtNiP4/Q259z1wHrgE5FSgiRW2v/smlkxsBX4W+dcR6rHM5o4403KPk6ba5465+om8LY2YEHM9/OB44kZ0diNdexm9n3g10kezkSkxX4cD+fc8cjtSTP7JeFy0gupHdWo/GY21zl3IlJzPZnqAY3EOeeP3k/Hn10z8xEOyR875x6LPJy2+zjeeJO1j9Nmxj5BTwAfMrM8M1sCXA28nOIxXSLywxV1O7BvuNem0CvA1Wa2xMxygQ8R3rdpycyKIgegMLMi4D2k534d6gngnsj9e4DHUziWUaXzz66ZGfADoMU59+2Yp9JyHw833mTt40zpirkd+C5QAZwDdjvnbo089wXgL4EA4T9vnkrZQOMws/9L+M8sB/wJ+Hi0BphOIm1W/wPwAg85576S4iENy8yWAr+MfJsDPJJu4zWzR4GbCK/e5wf+G/Ar4GfAQuAIcKdzLi0OWA4z3ptI059dM3s78CKwl3DHHMB/JVy3Trt9PMJ47yIJ+zgjgl1ERMYu00sxIiIyhIJdRCTLKNhFRLKMgl1EJMso2EVEsoyCXUQkyyjYRUSyjIJdRCTL/H8ZoPIsKcmdxAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(Y_train, Y_pred_train, alpha=0.5, s=4)\n",
    "plt.scatter(Y_test, Y_pred_test, alpha=0.5, s=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ensemble method for error estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zhuangjw/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:1044: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  output = mkl_fft.rfftn_numpy(a, s, axes)\n"
     ]
    }
   ],
   "source": [
    "# variance (calibrated) as an indicator of error\n",
    "rf_error = fci.random_forest_error(model_rf, X_train, X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(69.44043213133598, 67.16287773993044)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_error.max(), rf_error.mean() # this is huge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1a9e6ba8>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8lOWZ//HP9cwpkwNJIBhCgAIKEomKHFSqlaKArK7VKtittvXQV7H18MJzrdqK3fWw1mrb3ZYWK1u7+rNVsV1dLR7Qql1tJWiRQAoIRCCEACHnzPm5f3/MJCaYhBwmmUOu96tpkieTmWsy5pub67mf+xZjDEoppdKHlegClFJKxZcGu1JKpRkNdqWUSjMa7EoplWY02JVSKs1osCulVJrRYFdKqTSjwa6UUmlGg10ppdKMMxEPWlBQYCZOnJiIh1ZKqZS1YcOGQ8aY0Ue7XUKCfeLEiZSVlSXioZVSKmWJyCe9uZ22YpRSKs1osCulVJrRYFdKqTSjwa6UUmlGg10ppdJMr4NdRMaLyJsiUiEim0Vkeez4ChGpEpG/x97OG7xylVJKHU1fpjuGgVuMMR+ISA6wQURei33tUWPMw/EvTymlVF/1OtiNMdVAdezjJhGpAIoHqzCllEonhw4dIhwOM2bMmEF/rH712EVkInAK8LfYoetF5CMRWS0i+d18zzIRKRORsoMHD/arWKWUSjXGGJ544gmmTZvGddddNySP2edgF5FsYA1wozGmEVgJHAvMIDqi/3FX32eMWWWMmW2MmT169FGviFVKqZS3fft2FixYwJVXXsnxxx/PvffeOySP26clBUTERTTUnzLGPA9gjKnp8PXHgP+Na4VKKZWCgsEg8+bNo6WlhZUrV7Js2TIsa2gmIvY62EVEgMeBCmPMIx2OF8X67wBfBsrjW6JSSqWODz74gBkzZuB2u/ntb3/L9OnTKSoqGtIa+vLn4wzg68DZR0xtfEhENonIR8B84KbBKFQppZJZfX0911xzDbNmzeKJJ54AYMGCBUMe6tC3WTF/AaSLL70cv3KUUiq1GGN45plnWL58OQcPHuTmm29m6dKlCa0pIcv2KqVUurjuuutYuXIls2bN4uWXX2bmzJmJLkmDXSml+iocDhOJRPB4PFx88cVMnTqV66+/HqczOSJV14pRSqk+WL9+PXPmzGmfurhgwQJuvPHGpAl10GBXSqleaWpqYvny5Zx++unU1NQwZ86cRJfUreT5E6OUUknqrbfe4mtf+xpVVVVce+213HfffeTm5ia6rG5psCul1FEUFBRQWFjIs88+y+mnn57oco5Kg10ppY4QiURYuXIlmzZt4le/+hXTp09n/fr1RK/TTH7aY1dKqQ42btzI5z//eW644QYqKyvx+/0AKRPqoMGulFIAtLa2cvvttzNr1iwqKyt56qmnWLt2LRkZGYkurc802JVSiuisl1//+tdcddVVVFRUcNlll6XUKL0jDXal1LBVU1PDihUrsG2bwsJCtm3bxmOPPcbIkSMTXdqAaLArpYYd27Z57LHHmDZtGg888AAffvghEJ39kg402JVSw8qWLVuYN28ey5Yt4+STT+ajjz5i1qxZiS4rrnS6o1Jq2LBtmyVLlrB//35Wr17NlVdembJ99J5osCul0t7bb7/N7NmzyczM5KmnnqK4uJhjjjkm0WUNGm3FKKXSVm1tLVdddRXz5s3jpz/9KQCnnHJKWoc66IhdKZWGjDE8+eST3HzzzdTX13PHHXewfPnyRJc1ZHTErpRKO7fffjvf+MY3mDJlCh988AEPPPAAmZmZiS5ryOiIXSmVFoLBID6fj9zcXK688komT57MNddcg2UNv/Hr8HvGSqm08+677zJz5ky+/e1vAzB9+nS+853vDMtQBw12pVQKq6+v5zvf+Q5nnHEGjY2NXHbZZYkuKSn0OthFZLyIvCkiFSKyWUSWx46PFJHXRGR77H3+4JWrlFJR7733HiUlJaxatYqbbrqJLVu2cMEFFyS6rKTQlxF7GLjFGFMCnA5cJyInAHcA64wxU4B1sc+VUmpQGGMAmDx5MtOnT+f999/nkUceITs7O8GVJY9eB7sxptoY80Hs4yagAigGLgSeiN3sCeCieBeplFLhcJgf//jHnHvuue2Ldr3++utptxxAPPSrxy4iE4FTgL8BhcaYaoiGP5DeM/+VUkOurKyMU089lVtvvRWPx0NTU1OiS0pqfQ52EckG1gA3GmMa+/B9y0SkTETKDh482NeHVUoNQy0tLdx4442cdtpp7N+/n+eee44XXnghqTeSTgZ9CnYRcREN9aeMMc/HDteISFHs60XAga6+1xizyhgz2xgze/To0QOpWSk1TFiWxcsvv8y3v/1tKioquOSSS9Jy0a5468usGAEeByqMMY90+NILwBWxj68A/id+5Smlhpuqqiquu+46Wltb8Xq9fPjhh/z85z/XUXof9GXEfgbwdeBsEfl77O084EFgoYhsBxbGPldKqT6JRCL8/Oc/p6SkhNWrV/P+++8DkJWVleDKUk+vlxQwxvwF6O7fQOfEpxyl1HD00UcfsWzZMv72t7+xcOFCVq5cybHHHpvoslKWrhWjlEq45cuXs3PnTp588smU3kQ6WWiwK6US4tVXX+Wkk05izJgxrF69mhEjRjBq1KhEl5UWdK0YpdSQqqmp4bLLLuPcc8/lRz/6EQCTJk3SUI8jDXal1JCwbZtf//rXlJSUsGbNGu655x7uv//+RJeVljTYlVJD4v777+db3/oWJ554Ihs3bmTFihV4PJ5El5WWtMeulBo0gUCA2tpaxo4dy7e+9S3Gjh3LlVdeOWzXSR8q+tNVSg2KP//5z5x00kksXboUYwyFhYVcffXVGupDQH/CSqm4qq2t5eqrr2b+/PmEw2Huuecenb44xLQVo5SKmw8//JBFixZRX1/PHXfcwfe///1htYl0stBgV0oNWDgcxul0Mm3aNBYsWMCdd97JiSeemOiyhi1txSil+i0UCvHAAw8wY8aM9kW7nn76aQ31BNNgV0r1y3vvvcfMmTO58847mTZtGj6fL9ElqRgNdqVUn/j9fq699lrOOOMMGhoaeOGFF3juuef0ytEkosGulOoTt9vN5s2bWb58OVu2bOGCCy5IdEnqCBrsSqmj+uSTT/jqV7/K/v37sSyLdevW8eijj5KdnZ3o0lQXNNiVUt0Kh8M88sgjnHDCCbzwwgt88MEHADidOqEumWmwK6W6tGHDBk477TRuueUW5s+fz5YtWzjvvPMSXZbqBf2zq5Tq0qOPPsq+fft45plnWLJkiV49mkI02JVS7V588UUmT57M9OnT+elPf4rD4SAvLy/RZak+0laMUop9+/axZMkSvvSlL/Hwww8DMGrUKA31FKXBrtQwZts2v/jFLygpKeGll17ivvvu41e/+lWiy1IDpK0YpYaxlStXcv3117NgwQJWrlzJcccdl+iSVBz0esQuIqtF5ICIlHc4tkJEqkTk77E3PWWuVJLz+Xz84x//AODqq6/m97//Pa+++qqGehrpSyvmN8DiLo4/aoyZEXt7OT5lKaUGw2uvvUZpaSnnn38+oVAIr9fLpZdeqjNe0kyvg90Y8zZweBBrUUoNkgMHDnD55ZezaNEinE4njz/+OC6XK9FlqUESjx779SLyDaAMuMUYUxeH+1RKxcm2bds4/fTTaW5u5gc/+AHf+973yMjISHRZahANdFbMSuBYYAZQDfy4uxuKyDIRKRORsoMHDw7wYZVSR9Pa2grAcccdxze+8Q02btzIvffeq6E+DAwo2I0xNcaYiDHGBh4DTu3htquMMbONMbNHjx49kIdVSvUgEAiwYsUKJk+e3L5o109+8hNKSkoSXZoaIgMKdhEp6vDpl4Hy7m6rlBp8b731FieffDL33nsv55xzDg6HI9ElqQTodY9dRJ4GvggUiMhe4B7giyIyAzBAJXDNINSolDqKSCTCsmXLWL16NZMmTWLt2rWce+65iS5LJUivg90Y89UuDj8ex1qUUv3kcDgIh8N897vf5Qc/+AGZmZnd3vbZsj0ALJ09fqjKU0NMlxRQKkXt2LGD888/n02bNgHwm9/8hgcffLDHUFfDgwa7UikmFArx4IMPUlpayjvvvMP27dsBenWRUUV1A69u3s/v1+/h0de2UVHdMNjlqgTQYFcqhfz1r39l1qxZfO973+Of/umfqKio4OKLL+7V91ZUN7Dq7V34ghHyvE4afCFWvb1Lwz0N6SJgSqWQF198kcOHD/PHP/6RCy+8sE/fu7a8hlyvi2Z/CIBcr6v9eElRbtxrVYmjI3alkpgxhjVr1vDGG28AcPfdd1NRUdHnUAeoqveRk9F5LJeT4aSq3heXWlXy0GBXKknt3r2bCy+8kCVLlvCf//mfAHi9XnJycvp1f8V5Xpr84U7HmvxhivO8A65VfdazZXvaZyANNQ12pZJMJBLhJz/5CSeccALr1q3jRz/6Ec8888yA73dxaSENvhC+YARjDA2+EA2+EItLC+NQtUomGuxKJZk1a9Zw0003cdZZZ7F582ZuvfVWnM6Bnw4rKcpl2VmT8Lod1PvC5HpdLDtrEiVFuQkdXaajRM8+0pOnSiWB5uZmNm3axNy5c1myZAlr165l0aJFcV8nvaQol0XTxwB6gdJg6W72Udsf0aGgwa5Ugr300ktce+21NDU1sXv3brKzswd1OYAjA71tdFnXGmJvnY/FpYU6S2YAkmH2kbZilEqQ6upqli5dyj//8z+TnZ3Niy++SHZ29pDWoHPb4y8ZZh/piF2poxiMtVWqq6spKSnB7/fzb//2b9x222243e643X9vJcPoMt0U53lp8IU6HRvq2Uc6YldqCNXW1gJQVFTEXXfdxaZNm7jrrrsSEuqQHKPLdJMMs4802JXqQbxmN/h8Pu68807Gjx/fvmjXbbfdxpQpU+JZbp/p3Pb462n20VDRYFeqG/HqP7/++uuceOKJPPDAA3zlK1+hqKjo6N80RJJhdJmO2mYffWXOeG5aOHXI21oa7Ep1o63/7HU7EBFyvS5yvS7Wltf06vuNMXzzm99k4cKFWJbFG2+8wX/9139RUFAwyJX3XjKMLlX86clTpbpRVe+jKLfzxs+96T8bYxARRIQJEybw/e9/nzvvvDNpN5HWue2DI5E/Sx2xK9WN/vSft27dyvz581m7di0A99xzDz/84Q+TNtRVetJgV6obfek/BwIB7r33Xk466SQ2btxIY2NjAiruv6Wzx+toPY1osCvVjd72n//yl78wY8YMVqxYwcUXX0xFRQWXXnppgqpWSnvsSvWoN/3niooK/H4/f/rTn1i8ePFQlqdUl8QYM+QPOnv2bFNWVjbkj6tUPBhjePrppzHGcPnll2PbNn6/XzeRVoNORDYYY2Yf7Xa9bsWIyGoROSAi5R2OjRSR10Rke+x9fn8LVioV7Ny5k8WLF3P55Zfz29/+FmMMlmVpqKuk0pce+2+AI/+deQewzhgzBVgX+1yptBMKhfj3f/93SktLeffdd/nZz37Gyy+/HPdldZWKh14HuzHmbeDwEYcvBJ6IffwEcFGc6lIqqbz77rvccccdnHvuuVRUVHDDDTfgcDgSXZZSXRrorJhCY0w1QOz9MQMvSank0NjYyIsvvgjAvHnzeP/99/nDH/7AuHHjElyZUj0bsumOIrJMRMpEpOzgwYND9bAqTobT1mnGGJ5//nlKSkpYunQpNTXRJQTmzJmT4MqU6p2BBnuNiBQBxN4f6O6GxphVxpjZxpjZo0ePHuDDKjU49uzZw0UXXcQll1zC6NGjefvttyks1AWxVGoZ6Dz2F4ArgAdj7/9nwBWppDNctk5rbGzk5JNPxu/389BDD3HjjTficrkSXZZSfdaX6Y5PA+8Bx4vIXhH5JtFAXygi24GFsc9VGhkOW6dVVlYCMGLECH72s5+xefNmbrvtNg11lbL6Mivmq8aYImOMyxgzzhjzuDGm1hhzjjFmSuz9kbNmVIob6NK1yay5uZlbbrmF4447jldeeQWAr33ta0yaNCmujzOczk+o5KBLCqge9Xfp2mT30ksvce2117J7926WLVvGqaeemuiSlIobDXbVo2TYmDferrnmGlatWsUJJ5zAO++8w5lnnpnokpSKKw121aPFpYXtPfYMl9W+dO1X5qTWXG7btgGwLIu5c+cyfvx4br/99oRtIq3UYNJle1WP0mHrtPLycr7whS/w2GOPAXDllVdy9913D0mox2szbKX6QoNdHVWiN+btL5/Px1133cUpp5zC1q1byc0d2rqHw4wilZy0FaPS0jvvvMNVV13Fjh07uOKKK3j44YeHfBPpthlFzf7oOYpcr6v9eKr8cVSpSYNd9UqqbZvm8/mwLIt169Zx9tlnJ6SGdJ1RpJKfBrtKC8YYnnjiCfbv388dd9zBokWL2Lx5c0IvMkrHGUUqNWiPXaW8bdu2cfbZZ3PVVVfxyiuvEIlEABJ+5WhfNsNWKp402FXS6elKzY5fCwQC/PCHP+TEE0/kww8/5Je//CXr1q1LmnXS02FGkUpN2opRSeXZsj2srzzMnIkjj3rbnTt3cu8P/5Upp57Nm2t+Q1FR0RBU2De92QxbqXjTYFcppbmxnvffXMvS2bdRUlLC1Y8+T37RhKQMdaUSRVsxKiUYY/jd737HTUvPYdUD32Pbtm0A5BdNSHBlSiUfHbGrpLKv3kd5VQM7D7a0r/2e4T/Mtddey7qyCoq/cDmnz13AS59AJCc1LvTRFowaahrsql/aTmDGM7Re+qiKp9/fTUsgTEG2m8pDzfzqz008d993aG5q5pzrH2Bk0US8bkf7VZwtgRB5mbrei1IdabCrpFBR3cB/rNuBAB6nRUurj601Do4vzOHLN9xLQUEBxpVJVV0rEL2Kc9fBJt6vrMNhCZW1rVwxdwLnn1Sc2CeiVBLQHrsaEkfbbGJteQ1h2+CyDHU1e9m9dRONdYepbvDjGTmWJttNTsan45Bt+xv5cHcD4YiNxyE0+kI8+KetvPRR1VA8HaWSmo7YVZ8Nxh6oVfU+bF8jO3dWEva3kD96DGOOKaC2JcjcY6NrvHS8inNjVQOWJTgdFmJZ7euwPPHebh21q2FPR+yqT+K1YuGzZXu4/bmN7aP4TX/9Mx/99c843R4KJk5jzLgJtIYMLofF4tLCz1zF2eQLY4nB4/z0P+Ecj4OaRn9cn69SqUiDXfVJb/ZA7e0en3YkQigYAOBfzpzGrJmncNIJ08jwZtIStMHADWcfS0lR7meu4szxOsnNdON0fPqfcFMgQuGIjO4eTqlhQ1sxqk/itWJh5dbNPHXPrew4cx6Xff7nXPXlczn99AYeeXUbjYEIU47J5uZFndd+73gVZ6bb4sE/bSUYtnFZ0TZNSyDMDWcfO/AnqVSK02BXfdJxxcK2MC/Oz+xxxcKOPfldBxrY/uZzrH7038jIzmXi8dPbb9cW3PlZ0SUFeurbt/XR73lhC/5QhBFeFzecfaz215UiTsEuIpVAExABwsaY2fG4X5V4R85XX1xayF1/KKfRF8LlEHzBSI97oHbsyYcbD/KLx9bQKpmc8M2HyJtYSvPYfCqqG/p18vX8k4p5a9shAB5acnI/n6FS6SeePfb5xpgZGurpraQol3lTR+N2WvhCNl63o8cVCzv25N0eD05vJseeMpfMCaX4wrBxbwM3/W5jp2mKcyaO1Ks1lRoAbcWoPhub56W0OBrkR7ZMOrZd9hxu5U/vfYRdt5fJp5xJbv4oxs88mxZ/iCZ/hAyXgxy3BQL/8cYOJo/OPupja+ArdXTxCnYDvCoiBviVMWbVkTcQkWXAMoAJE3ThplSyvvIw0DlUy6saqG8NdVrTBWhvu4i/gcf/+xUOOwsYaYX5XDiMw+kkELLxhWycluCwBBFhRIaT2pYga8truGnh1D7Vpi0YpT4rXsF+hjFmn4gcA7wmIv8wxrzd8QaxsF8FMHv2bBOnx1X9NJC1XvbV+/ikthVLIM+b0z6X3euyyHYJf/3ofbZuXI/H42HmvBOJZB9D0BYyjMEhQjBiGJXlIhiJ/mcQCNuMynLrXqBKxUlcgt0Ysy/2/oCI/AE4FXi75+9Sya6trbL9QDN1LUFKi0cA8Me/VxGxDeIQmv3h9qs+/7arlrkTstix5e+MP/Z4Lr3ofDK8mVRUNyJAXWuI4rwMfKEIxgDGEIoYAmGbiaN6nlmjlOq9AQe7iGQBljGmKfbxIuCHA65MDbquWixtOs5m8bosfMEID63diiVCcyCMQ8A2hoqaJoJBP/aBnVAwFduRwblLryDDm0lmZhYNvhDTx+YyLv/T0N64p47XKw4QCNvkZMDxhdlYlqV7gSoVJ/EYsRcCfxCRtvv7f8aYtXG4X5UAbS2avXU+cr0umv0hpFXwuh0caIpeJZrjcdISCEdDvv4wf/14L4H9O7j4kok0+EIYh6fT5s1fmTOO8qrG9rXW61pDTByZyYHmACOzPEwsyI7LejNKqagBB7sxZiegZ7DSTNsVpk2+EAeaAtS3hghGbLI9DnIyXPiDYcIBH7bfhzO3kMVnzuHmL50CwCOvbqOuNcR0r4uvzBlHSVEur26u4a1tB8l0OcjzOvGHbCI2fGFKAcsX9O2EqVKqZzrdcRg6sne+r973mROXW/Y18HFNI/sbA0Rsg8cjBMJQ2xykriVEqLUZcbpw5ozC5XRyydwp7SPurjZvtg2cNC6vfT11r9uByyGUV6XGLkhKpRIN9mHk2bI97Kv3UVnbysGmAPvqW9lb18r7lYcZmemiICeDSQVZ7bffdaiVTLeDUNim3hfGbUL4bIvWcARnZhYWFojgtIT/+7i2x8v5u1pjxuUQ6lpD3XyHUqq/NNiHmd+v30NNg59I7HMBAuEI+xttAiGbQ00BCnI8ZHmcZHmceJxCXUsEX+1e6rd/QP7MRYjlxiFOwGAMhG3Dq1tq+Prc7pcG6LjGDMAJY3OjJ2bdjkF/zkoNNxrsw8SzZXv4xZsfs7/Rj93huAC2bbANNPpDOEMRmgJh3A4Lj1OQxhoO/XUdkaCf0+eeTrXXS8Q2hG1DMBz9PttAsz/Mw69s49Zzp3Y5y2ZxaWH7LJsMl0WDL8ToERksO2vSkP0MlBoudD32NNe2ocW//+kf7DncGp0/3kHbp4KhNWTjC0UIhSM0+UNUfbKD8s1bcOWM4sRzLuG0M+eTn+nGsoRg2CZ2fRECZLod7DrUwpPvfdJlHUeup57rdfW4xoxSqv90xJ4mbn9uIxC9xH7pL99tP76hsq79Y8OnQd7xWNiOjroBAmFDMBwGBEf2GEYdl0dm9ggsr5tcr4tbz53Kf76xgwpfEwJYApYIY/O8WAIf7un+ZGjH9dR1zRelBo8GexradbCFJn8Il8Pq1HbpTluoY2yMAWNsxOEkgoXtyCTkC+N2Rnvhk0dns7h0DNsPNGFJdL0Xt9Miy+PEHwxjPvOnQyk11DTY08SuQy1U7Gvk5Y+qaQ5GT41G7MhRvutTxhhAov+TaIgbQEx0kf3WYLh9TZhlZ03ijx/u5VBzEKfDwhLwhyI0ByKcNnlk/J+cUqpPNNjTwPwfvcmu2uj8cOlw3O7F4Nkc0XSPXkD86b04HUIoYvCH7fY1YdaW1zClMAdfqIFAyCYcu48JozL5+tzPDeSpKKXiQIM9xbRd8r++8jBvVBzgcEuwU7vFdPPxkYwxiAgSu100z+Uzt7ONicV89P/b9jfNy3Rz0rhc9tb5aA1G+OLxx/RqWQDtrSs1+DTYU9Av3vy4fYTek66CvdMI3RiIhXt3wjZ4nEJBlhuAJn+Y4jxv+6JebQuJ9XUddaXU4NFgT1LPlu1hfeXh9m3ibn9uI7sOtfDRnnoCkf6doDyy7dI+Uu+BAG6Hxbj8jE6LerWNzNuCXSmVPDTYB9FANrNo0xacz2/YS7ifE06ODHTovvVypAkjvYQjhogRcjss6tVmzkQ9WapUstFgHwL9Cfi2Hro/FOGZsr1xrUd6EegA2W6LL80o5r0dtXxlzvgu69eeuVLJR4M9CT1btoddh1o41BLs9310HKWLSPvJ0t5wxM6oHjMigyZ/mPxMV7/rUEoNPQ32JHL7cxs7zUPvj67aLtD7Ubol4HFZiIleUdrgC1FarJf9K5VKdK2YQdJ28vNot2lr09z+3EaeLdsb91DvTaA7AK9T8Lospo3JIc/rxutxMqkgi2VnTWKs7kWqVErREfsQWF95mF2HWoBPe9JLf/lu+xK5z5TtYX2HNV3ipbej9IIcF/6QAYH61hClxSNYeMKY9lp1oS6lUosGey/19QRo2/6eOw+2sK/ex5HLKh5qCnAwtodob+akd6dtlN4xxHsb6G2aAxFGZrkJRwxh22br/mamj/Ud/RuVUklJg30QVFQ38Na2gwTDNmNznWzdH6I1GOHRV7fyTNkeDjUF2HO4lbCB5n6Gen/bLkfKcFpYCC6HRcS2KchyMzLTrVvWKZXCNNgHwdryGhp9IQ41B6j3hfCFIkQM7GsMcKAp0O/56ND9ydG+cgh4XQ4cDiHH46TRH+aU8XmcMDaXzVX1umWdUiksLidPRWSxiGwVkY9F5I543GeyWV95uNcnQ6vqfbgcQjhiaAmEiXRYzGUgod4VEenXSD0/y80IrwuXw2JUlhuX49P/FPwhW6c4KpXCBjxil+garz8HFgJ7gfUi8oIxZstA73swxeOq0CO1Bb/HIexvDNAajID0vBhXb3TVR+/q86OR2JvTgmNyPACM8Dg52BzEaUUfR7esUyr1xaMVcyrwsTFmJ4CI/A64EEjqYO9Od4HfNqvlaOpbg4QihnDEjq66OMBU79h66ctFRkfyOC1cVnRjjOh+pTYuh+B0WEwYlYkvGKHeF2Z6F8sGKKVSSzyCvRjY0+HzvcBpcbjfpFFR3UBVnY9gxObR17a1L0975B+BXYdaqKrzsWj6GDbtracpMPRz0jvdnugFRwKcNaWgfVu6thk72w8043U7uHnRVMqrGjs9F6VU6opHsHeVNp9JJRFZBiwDmDBhQhwedmhUVDew6u1dhG0blyWddhGCT9svbYEYjNjkZDhjW8n17wRkvGa8OAQcDsFpWZ2uHh2b52VsnpdF0z+keTDSAAAPn0lEQVStuy3YlVKpLx7BvhfoOMwbB+w78kbGmFXAKoDZs2enzMaYa8tryPW6cFrRk4sddxFqW5O8o3DE5vfrd9PkD/f5sbpqtfQn0LPdFhEDhSMymHJMNqXFue1Xj7YFedu/NpRS6Scewb4emCIik4Aq4F+Ay+JwvwM20BOkz5bt4f8+PsTcY0e1H9uyrwFjDM1Znk7B3tauCYQiNPnD9GXJ9Hj10dtMGJXFgcYAp04a2b6srrZYlBo+BhzsxpiwiFwPvEJ02ZHVxpjNA65sEFVUN/Dq5v3UtYbYW+f7zJZuHac15me62F3bQqM/RMQ2fLinHq/L4uKZeTxTtoddB1uobw3y7o5amvwhgmFDXzrr8ZqX3ibLFW27rA92bhEppYaPuFygZIx5GXg5HvcVLx13IOqorWfuC0bI8zrZUHmYd7Yf5L4vl/Lq5hpe2byf1mCEupYgpcW5FI7w8Na2Q4RtGwdCMGzjC0aYWpjFuzsOAbC3zsdxox0EI3avQ32gqzB2Jdttcd5JY9ufc182wdA/AEqlj2F35ekjr27DF4zgdTsA2t8/+d4nVOxvIhi28bosfMEIb207iNsSJhdksr/BT8gYrEiEDKfwfx/XsutQC3WtQepbg+yt8xEcwNVHA22/ZDilfUek3oS0BrlS6WvYBXtda4g8b+enneGy+HBPAyMzXfhD0TF3W+D/Y38TBdkeXM7ops8GCERsXtlcQ3MgFF3bS6A5ED7qlPWBbH7RE49TGOF1cfK4PB5acvKA708pldrSNtg7rq7YsY+en+nqNGKH6CX0BkOGq/MKCxkuC384evWoJUIkYhO2IWRDXUsQt1OwLLBN9K07g9F26SjT7cTrcnLzoqntx3RErtTwlZYbbXRcXTHP62yfe15R3UBpcS6toQi+YARjDL5ghNZQhFPG5+EPfbqoS5M/zNaaJkIRw8EmP/5QhEDEEIoluA0EIyYa6H0M9XgEutsh5Hmd5HqdFOVmcPHMYr1aVCkFpOmI/cn3PuFgk59GfxjbGApzM8j1utrnns+bOpryqgaq6v34QmG8zujovbYlSHMgRKMvxMdhG6cl5Hgscrwuals+e7FRT6P07gwk1C1g1sR8AC6NjcjbZvDoLkdKqTZpF+wV1Q2883EtxoDLEkIRw4ZP6jllQi5bqlvZsq+ButYQTomOpotzvWS4LFoCYfbV+2gORC8s8jgtcjxOwraB6OZCvV/2pW2UPoDNL45UkOnk7BPGMGfiyE7TMedMHKltF6VUJ2nXillbXkN+piuaqSK4HILHabGhsp49tb72aY5V9X4ONgcI24bmQIStNc1EbLs9i8O2wR+2yXQ7yfI46VUuG4Mxho7dmf4uq9uR04IvHF/IVWdMHND9KKWGh7QL9qp6HyVFOYQihogdDVqMYX+Tn6mF2XjdDkSEiDFkOC0ONPk50OTH47RiYRwN4UhsBUSnFb1tL6a8dLrJQLvoTis6L90h4HJY7a0kpZQ6mrQL9uI8LxkuJxNHZuKwhEDYgAijMt18riCr/Xae2AwYX8imtiVIoy+EJdGGixVL5YhtCEVs3A6r/ViXOoR69B8KQu+G+J9lSbSFNC7Py4RRWRSOyGBSQRY5GU6q6nUfUqXU0aVdsC8uLaTBF8LpsBid7WZyQRaTR2fz+WNHdVqYqzDbExuRg9MS/KEIHqcjusytCE4rOsWxNWRTOnYEWR5HFz+sT3vpghlQoC8/ZwrHjc5izIgMphRm44ztaJSX6aI4z0uTP0xx7ASp9tWVUj1Ju5OnJUW5LDtrUvQK03obr9vRvsRu21ICGS6LaWNz8Xqc+IIRdh9uxTaGcfkZ7KkzNPlDgFCU6+GiU8axfMFUXtlSAxJBDJhIEBMKgjsTaRvKS9//Rjokep41J8NBgy9EKGIoLR7B9LG5PP9BFWAzcVQmxfmZNPhCn9kAQ8NdKdWVtAt2iIb7ouljyM+KrhXTFoZtgV/XGmK618Xti4+nvKqRZ8r2kJ/porE1RNg2iERbL/6QzSub9wOQ7XESDIZoaGwkEg5hYXC4PRgsMpwWgsHXyyUFBJg9MZ/8TBfv7zqMM9ZDP35MNnmZbpYvmMqi6YWdatVdjZRSvZWWwd6dtsCHrjeYCNqGDKfQ7I9eqFTbEiRi26wt38+hhmbqDh1AnB4sgUnjxlLbEqIlaCMY7F6s0ijAqCw3EJ2HvnT2eJb+8l0Ablo4tdMa6V3VqpRSvTGsgr0r++p9VNX52FbTRIZTONwS7cO3tcrrWkMEIoasDA9Bbwa2J5sITpoCNg7LomhE9Ee4vynY5f0L0bVcJo+ObnjxRsWBTl+f1OGErlJKxUPaBXvHUe/RTjK2LT0QtqMj9NagTdg2OGL7hIYDrYRtC49DcLg9zDx+EuXVjfhDNrmZLqYck01da4ja5gAZrjCBkI0lELajp1WdFozP9/KLr836zNZzbXV1XFpXR+ZKqXhIu2Bv05uQXFteQ6bLgdOysB3QGokggG1swk312CI4vSPwhQ3OYHQkP2FkJnUtQUZlecjLdLO3zkdNY6B9TRgRIcMlZLodTCrI5rnvfB6gvT+u4a2UGmxpG+zd6RisVfU+xo308o+aJnzBMMaAHQ5E14BxuHB7sxCsWDvFovJwK16XRcQY8jNdPLTkZB59bRvPb9hLbWsApwNyMlx4XQ6aA2HyM119qqcvX1NKqe6k3Tz2vnA7hPW76ojYhgyXg4JsN0h0LrvLm4PTYeF2WWR6HLidDpyWEI6A1+WktDg6Aq+q9zE2z4Mx4LCEHI8DY6JXrbbdRimlhtKwG7F3JEAoHKKu+hNyRxUyMn8kjf5wtL9uG0ZnuynMyQARDjT5aQ1Ge+hLZkbntkP0Stf99T7GjPDQFIgQCBscIhyT49YVF5VSCZFWwX60Tao7CofDfPDRJraufx8rbyyBnHwyXA6mFeZQ0xQg0+1gZKa7fUOOnIxsfMEIsyaObA91iF7puvtwKwC1LQFGZXkYPSKDZWdN0nnnSqmESJtWzJGbVHfcXONIGzZs4LTTTuO911+iuKiQE0pLmVx8DKdPHoXLYZHpdjAu39vlhhyLSws73Vfbla5etwNf6NMrXTXUlVKJkjYj9rXlNeR6XTT7oxti5Hpd7cePDNnnn3+effv2ce+D36My4zg27a3HGEODL0RrKMK4fC95mW6mj82lvCq6fnt+pos5k0Z2GdjdXemqlFKJMKBgF5EVwLeAg7FDdxpjXh5oUf1RVe+jKDej07GOKyK++OKLZGdnM3/+fO6++25uu+028vLyqKhu4JZnNtIajJDrdXHfl0vb55wvnT2+07x4pZRKBfEYsT9qjHk4DvdzVG0h29U0wOI8Lw2+ztvXNfnDjHCEWLJkCWvWrOGCCy5g/vz5eL1evN7oic2Sotz22Ss3LZzafuxIvZl6qKsuKqWSQcq1Ytq2hTsyQBeXFnZavbG+NcD6jVt477/+Ff/+Hdx///3ccsstXd5nx6s/lVIq1cUj2K8XkW8AZcAtxpi6ONxnn3VcrreuNYQ07OD1//gu82ZMYeVrf+S4444b1MfXkbpSKlkcNdhF5HVgTBdfugtYCfwr0aVR/hX4MXB1N/ezDFgGMGHChH6W27OOKyJeMvN0LjrOxXnnnTegPUc1sJVSqeaowW6MWdCbOxKRx4D/7eF+VgGrAGbPnt27hcsHwLIszj///F7dVsNbKZVOBjorpsgYUx379MtA+cBL6lpFdQNryvZQUdPEX7Yf4u+76/ja3M/p1EKllDrCQHvsD4nIDKKtmErgmgFX1IWK6gYeWruVnYdasIhu+PzezsPsbwxw67lTNdyVUqqDAQW7Mebr8SqkJ2vLazjcEsTjsrABt9MiJ8PJoebAZy5A0raKUmq4S4klBarqfQTDNk7r05OgHqdFIBxpvwBJKaVUVEoEe3GeF7fTImx/es41ELbxOB0U6wqKSinVSUoE++LSQkZmuQmEbCIRm2DYpskfpiDb85lFuZRSarhLiWAvKcrl9sXHM7kgCxuwDcydPFJPnCqlVBdSZkmBkqJcLpk9nvWVh3VNFqWU6kFKjNiVUkr1nga7UkqlmZRpxYDOUVdKqd5IqWAHDXellDoabcUopVSa0WBXSqk0o8GulFJpRoNdKaXSjAa7UkqlGQ12pZRKMxrsSimVZjTYlVIqzWiwK6VUmhFjzNFvFe8HFTkIfDIED1UAHBqCxxlM+hySgz6H5DDcn8PnjDGjj3ajhAT7UBGRMmPM7ETXMRD6HJKDPofkoM+hd7QVo5RSaUaDXSml0ky6B/uqRBcQB/ockoM+h+Sgz6EX0rrHrpRSw1G6j9iVUmrYSftgF5EVIlIlIn+PvZ2X6Jp6Q0QWi8hWEflYRO5IdD39ISKVIrIp9nMvS3Q9vSUiq0XkgIiUdzg2UkReE5Htsff5iayxJ93Un1K/ByIyXkTeFJEKEdksIstjx1PpdejuOQz6a5H2rRgRWQE0G2MeTnQtvSUiDmAbsBDYC6wHvmqM2ZLQwvpIRCqB2caYlJp3LCJnAc3Ab40xpbFjDwGHjTEPxv7Q5htjvpvIOrvTTf0rSKHfAxEpAoqMMR+ISA6wAbgIuJLUeR26ew6XMsivRdqP2FPUqcDHxpidxpgg8DvgwgTXNGwYY94GDh9x+ELgidjHTxD9BU1K3dSfUowx1caYD2IfNwEVQDGp9Tp09xwG3XAJ9utF5KPYP1GT9p9uHRQDezp8vpch+g8izgzwqohsEJFliS5mgAqNMdUQ/YUFjklwPf2Rar8HAIjIROAU4G+k6OtwxHOAQX4t0iLYReR1ESnv4u1CYCVwLDADqAZ+nNBie0e6OJaKPbMzjDEzgX8Crou1CFRipOLvASKSDawBbjTGNCa6nv7o4jkM+mvhjPcdJoIxZkFvbicijwH/O8jlxMNeYHyHz8cB+xJUS78ZY/bF3h8QkT8QbTG9ndiq+q1GRIqMMdWx3umBRBfUF8aYmraPU+X3QERcRAPxKWPM87HDKfU6dPUchuK1SIsRe09iL36bLwPl3d02iawHpojIJBFxA/8CvJDgmvpERLJiJ4wQkSxgEanxs+/OC8AVsY+vAP4ngbX0War9HoiIAI8DFcaYRzp8KWVeh+6ew1C8FsNhVsx/E/0njwEqgWvaenTJLDYF6ieAA1htjLkvwSX1iYhMBv4Q+9QJ/L9UeQ4i8jTwRaKr8NUA9wB/BJ4BJgC7gaXGmKQ8QdlN/V8khX4PRORM4B1gE2DHDt9JtEedKq9Dd8/hqwzya5H2wa6UUsNN2rdilFJquNFgV0qpNKPBrpRSaUaDXSml0owGu1JKpRkNdqWUSjMa7EoplWY02JVSKs38fzi+jWNaPnkyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(Y_test, Y_pred_test, yerr=np.sqrt(rf_error)/5, fmt='o', alpha=0.5)\n",
    "plt.plot([-5, 25], [-5, 25], 'k--')\n",
    "# almost the same for all points..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1c1939bd68>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmQnPd93/n393n66XvuGQwAAiAG4E3JulAyJYJM1hIt2VYk59iitbHKTlxFpyqVw9FW1t6tbGprK1WuxE6iqmwlYsnXKorMjWzLLjtWJEaJCVAkRfCwSAIkwZnBPYPpufp+up/jt388PQMMOAPM9PSgex58X1UsYBo93b95OP15fs/3dzxijEEppdTuZ3W7AUoppTpDA10ppWJCA10ppWJCA10ppWJCA10ppWJCA10ppWJCA10ppWJCA10ppWJCA10ppWIicTvfbHR01Bw+fPh2vqVSSu16r7zyyrwxZuxWz7utgX748GFOnTp1O99SKaV2PRE5v5nnaclFKaViQgNdKaViQgNdKaViQgNdKaViQgNdKaVi4rbOcuk104UKJycXuFp0GR9Ic/zoCBNj+W43Syml2nLLHrqI/LaIzInIm9c99q9E5G0R+ZGI/JGIDO5sMztvulDhmVOXqLo+ewfSVF2fZ05dYrpQ6XbTlFKqLZspufwu8NkbHvse8AFjzI8B7wK/1uF27biTkwsMZhz6Mw6WCP0Zh8GMw8nJhW43TSml2nLLQDfGPAcs3vDYd40xfuvLF4EDO9C2HXW16JJPr6045dMJrhbdLrVIKaW2pxODon8X+PMOvM5tNT6QpuL6ax6ruD7jA+kutUgppbZnW4EuIv8H4APfuMlznhKRUyJyqlAobOftOur40RGW6x6lukdoDKW6x3Ld4/jRkW43TSml2tJ2oIvILwCfA/62McZs9DxjzNPGmGPGmGNjY7fcW+a2mRjL8+SxA+TSCWaLLrl0giePHdBZLkqpXautaYsi8lngfwP+ijGm1tkm3T4TY3kNcKVUbGxm2uI3gReA+0Xkkoj8EvDvgD7geyLyuoj8hx1up1JKqVu4ZQ/dGPPFdR7+rR1oi1JKqW3Qpf9KKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTtwx0EfltEZkTkTeve2xYRL4nImdbfw7tbDOVUkrdymZ66L8LfPaGx34V+G/GmHuB/9b6WimlVBfdMtCNMc8Bizc8/AXg91p//z3gZzvcLqWUUlvUbg193BgzA9D6c0/nmqSUUqodOz4oKiJPicgpETlVKBR2+u2UUuqO1W6gXxWRfQCtP+c2eqIx5mljzDFjzLGxsbE2304ppdSttBvofwL8QuvvvwD8cWeao5RSql2bmbb4TeAF4H4RuSQivwT8OvCEiJwFnmh9rZRSqosSt3qCMeaLG/zTpzrcFqWUUtugK0WVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomNNCVUiomthXoIvIrIvKWiLwpIt8UkXSnGqaUUmpr2g50EbkL+IfAMWPMBwAb+LlONUwppdTWbLfkkgAyIpIAssCV7TdJKaVUO9oOdGPMZeA3gAvADFA0xnz3xueJyFMickpEThUKhfZbqpRS6qa2U3IZAr4ATAD7gZyI/PyNzzPGPG2MOWaMOTY2NtZ+S5VSSt3UdkounwamjTEFY4wH/CHwyc40Syml1FZtJ9AvAI+ISFZEBPgUcKYzzVJKKbVV26mhvwR8C3gVeKP1Wk93qF1KKaW2KLGdbzbG/HPgn3eoLUoppbZBV4oqpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMaKArpVRMbCvQRWRQRL4lIm+LyBkR+USnGqaUUmprEtv8/q8A3zHG/C0RSQLZDrRJKaVUG9oOdBHpBx4HfhHAGNMEmp1pllJKqa3aTsnlCFAAfkdEXhORr4lIrkPtUkoptUXbCfQE8FHg3xtjPgJUgV+98Uki8pSInBKRU4VCYRtvp5RS6ma2E+iXgEvGmJdaX3+LKODXMMY8bYw5Zow5NjY2to23U0opdTNtB7oxZha4KCL3tx76FHC6I61SSim1Zdud5fIPgG+0ZrhMAX9n+01SSinVjm0FujHmdeBYh9qilFJqG3SlqFJKxcR2Sy6qQ6YLFU5OLnC16DI+kOb40REmxvLdbpZSahfRHnoPmC5UeObUJaquz96BNFXX55lTl5guVLrdNKXULqKB3gNOTi4wmHHozzhYIvRnHAYzDicnF7rdNKXULqKB3gOuFl3y6bXVr3w6wdWi26UWKaV2Iw30HjA+kKbi+mseq7g+4wPpLrVIKbUbaaD3gONHR1iue5TqHqExlOoey3WP40dHut00pdQuooHeAybG8jx57AC5dILZoksuneDJYwd0lotSakt02mKPmBjLa4ArpbZFe+hKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUTGuhKKRUT2w50EbFF5DUR+dNONEgppVR7OtFD/0fAmQ68jlJKqW3YVqCLyAHgZ4CvdaY5Siml2rXdHvq/Bf4pEHagLUoppbah7XuKisjngDljzCsi8ldv8ryngKcADh061O7bxcp0ocLJyQWuFl3GB9IcPzqi9xNVSm3bdnrojwKfF5FzwO8DPyEi//HGJxljnjbGHDPGHBsbG9vG28XDdKHCM6cuUXV99g6kqbo+z5y6xHSh0u2mKaV2ubYD3Rjza8aYA8aYw8DPAd83xvx8x1oWUycnFxjMOPRnHCwR+jMOgxmHk5ML3W6aUmqX03not9nVoks+vbbSlU8nuFp0u9QipVRcdCTQjTH/wxjzuU68VtyND6SpuP6axyquz/hAukstUkrFRduDoqo9x4+O8MypS0DUM6+4Pst1j888PN7llu0OOqCs1Ma05HKbTYzlefLYAXLpBLNFl1w6wZPHDmgobYIOKCt1c9pD74KJsbwGeBuuH1AGVv88Obmgx1MptIeudhEdUFbq5jTQ1a6hA8pK3ZyWXNSuoQPK29PpAWUdoO492kNXu4YOKLev0wPKOkDdm7SHrnYVHVBuT6cHlHWAujdpoCu1w3qhNHG16LL3hrGGfOtKpxdeT3WGllyU2kG9Upro9ICyDlD3Jg10pXZQr2zGdvzoCMt1j1LdIzSGUt1jue5x/OhIT7ye6gwNdKV2UK/Mne/0gLIOUPcmraErtYNWShMrg4bQvdJEpweUdYC692gPXakdpKUJdTtpD12pHbRSmjg5ucBsa5bLZx4e157tOnphNtBup4Gu1A7T0sStrcwGGsw47G2VqZ45dUnr8lukJRelVNf1ymyg3U576EqprltvoZLr+bwwNa8lmC3QHrpSqutuXKg0X3Z5cWqRlGPrXjFboD10pXaQDvRtzo07ab55pYQBPrC/f7UEA7pXzK1ooCu1Q3Sg7+ZuPNk9emSYc0t1ZosuDT/kk0eHGc1fK8PoXjG3poGu1A7RHQk3tt7J7vmpxdWT3ddfPE9V94rZMq2hK7VDemXZfy+61awWXZDVnrZ76CJyEPh/gb1ACDxtjPlKpxqm1G63nWX/ca+932r73Z1ckBXnY7udkosPfNkY86qI9AGviMj3jDGnO9Q2pXa1dm+Z18u19+lChW+/fpnXLy4DwkcODvCFD9+15XZt5mS3EwuyevnYdkLbJRdjzIwx5tXW38vAGeCuTjVMqd2u3R0Je3WRzXShwlefm+Kl1nTCdMLihalFnj4xteXphN0qqfTqse2UjgyKishh4CPAS514PbW7xPkSdrva6WX26t2ATk4usFhtMpBNknFsAESEhWpzywO93drjplePbadsO9BFJA/8AfCPjTGldf79KeApgEOHDrX1HhoYvSvul7Dd0Etb7l7vatGl6YcMXNeulGOxXGu2NdDbjT1uevXYdsq2ZrmIiEMU5t8wxvzhes8xxjxtjDlmjDk2Nja25ffolVt4qfXF/RK2G253OWK6UOHrL57nN/7rO3z9xfMbfrbGB9IkExauH64+1vBCUo69awJxo2N7eCizqWPQ69oOdBER4LeAM8aYf925Jq2lgdHbdGpe593OuwFtpcN0/OgIw7kkxVqTmudTb/os1ZqM5JK3bTrhZk8+G5kYy3N0JMOf/ugK/+77Z/nTH12hPyk8P7UYi07jdkoujwJfAt4Qkddbj/3vxpj/sv1mXRP3mtduF/dL2G7ZTjliKzNRtrL4aWIszy8/fmTNa3/iyHBbs1za0Yny3ol35vjGy5cYzac4MpajVPf5Ty9f5tF7hjk0nAV29wKwtgPdGHMSkA62ZV0aGL2t3al5cdQLYz0rM1HOzVfpzzqIgRemFrlabvDUY0fe156tdpgmxvL8yhP371j7b6YTK2+feSU6IQxmkwAMZpNYFrx1pcyHDg6vPm+3dhp7fqWorhjrbXqz4EivjPVcPxMl6yTIJBMMZZOrM1FudOMuh9C7HaZOlPdminX6M2tfYyDjsFRrrnmsV4/BrfT8Xi6dmt7UC72nuNI78qztPc6XXd4rVJkvN7i0VOPLT9x3247PZmaiXP9ZsC2YKze4ezi3Y1dYnfrsdeJqfd9AhlLdX+2hA2QSNmnHplT3dv1VZs/30CEKjC89cjf/62fu50uP3N1WmPdC70nF10rvcb7scur8Eg0/ZLQvyUK1cVt/1241E+XGz0LGSWBMdDOJnbjC6uRnrxNX609+7ADLdY/lWpPQhCzXmjRDw997bGLTV5nbHZjdST3fQ+8E3fVO7aTpQoVzC1VePb9EueExkHHIODauFzDWlyYMQ37ze+9yeCS341eHx4+O8NblIufmq5hWDX257nFkLMfxoyPrfhYOj+TIpRN86ZG7O96eTn72OnG1/tj9e4Coln5lORo/+KVHD68+fiu9vu7ijgh0nSmjNnLinTmeeeUSM8U6+wYyPPmxA5v+cMO1D3guaXGm5FKoNEgmLA4N++RTSfYNpDgzUyYwhkeOjOx4ANxqJsofvHr5tn4WNvvZ22xZphPlvcfu37Ol/8fX6/XO4R0R6DpTRq3nxDtz/OazZxnMONw1GNVWf/PZswCb/sCfnFwgDENmik32D6apewG1hseFxTp/7ccGWKh52JYwlEnetjvvrMxEuT4kVwZEt/JZaKf2feP32Ba3fL92e73dGBfr9c7hrqihb5fOlFHruX4KmyUWg9kkgxmHZ165tOnXuFp0mSm5ZJM2I/k0D+3vZ7QvTTZpM1/xKJRdAmO4Zyy3+j23Y+HVRrXrw0OZTX0WVr7/4nyV8wtV/vxHM/yzP36LE+/Mbfo9L85Xef69ef7kLy/zP96ZY65UX/f92lk82K1xsV6fFXRH9NC7tRGQ6m0zxTp3DWbWPNafSXBlefNhOz6Q5pULi+zpiz7Q+ZTDvoE0S3WPqUIFEcGxhfcKVQBG+9IdD4D1eqoblQbOLdU39Vk4OblAGIS8fbVMJplgrC9Fse7x9MlpDgxnN1ykFIYhp2dLzC7XWaw1Gc4m2dufAoHn3p3nsftG39fzbqfX207poxM9+l5fd3FHBDro1Dr1futNYSvV/feFy80cHsowV2pweanGQCbq4YsID433MVNu8ODePs7MlCnVPV4+v8SD431YttWxANioXFGuezy4v3/Nc1dCcjOfhZUrj0wysbqz4kDWYa7sbhia78yUOL9QI5dK4PoBFsJi1UNE+MJD+yjVPXLpxPu+dzNloBvD+O0rpQ1/vq0cp62OZfR65/COCXSlbvTkxw6s1sz7MwlK9ai39UuPHt7U908XKjw/tcjHDg7yysVlinWPWiPgY3cPcqXY4KF9fRwezZNPJ3ivUKVQdrlScjs6L32jnuqlpdq2xo3GB9K8en6Jsb7U6mMNL2Q0n9qwXLRcj8YL0o6N2wzJOHY0ptAMgI0D91a93vXC+MJSjVzS5vB1x/FmP9/KFcfpmRIl16M/7bC3L8W3X7/MaF96S73260+IKyeaP3j1ck+sb9FAV3es7U5hWwnTbNJivtbk4mIdPwwpVBqMZB3OzpX57+8UADg0nOGRiWG8KNv4N997Z9t3/YGNyxVDGYflurf69c1KAyuh9PaVEkXXYzDjMJpP0vADinWPgaxDwwupNQMeGM5vGJoDaYdi3aPuBaQdi5oXYGC1h3994N7Y4370yDCnLizxwtT86jG58Thff9J6cG8fp2fKDOdTmyp9vH2lxIWlGtlkgoG0g+uHvHJhiXoz4PMfvqutXnsvTmHUQFd3tK1OYbs+iF6aWsCyDOfnayAWaUfIOAkuLNTIJBMU602SCYtqI+DUeZe3Z8t8/PDQlvZauZWNyhX37etfraXPtlaEphPyvp7kSiiFQciFpRq2CMWaRzZps3cgzdWSSyMIGM2neGA4j2VZG04meGB/P7mkzWy5QdqxqTUDhnNJxgdSq4Ohn3l4fN0g/ObLF5kruyRtm5G8Q60RrIbjeietQyM56s1gdTHQRj/fiqLrYYusnlwyjk3J9bEtWbcOv/LnzXruvTiFUQNdqU26PogcGwqVBguVBpZAyhGqjRAMlBo+rhcgIizXfBxbsBEWq02+d/oqe/rT7B/MkHWij1+7d/2BqIb/1RPThKFhJJ9kX396tUa/Uhq4vt0rvdmVsFwJpdMzJbKternrBcyWG3zwrkHu2eNvuiRx/OgIzyzVeWhfPx+fGObCfJXTs2UG0w65dGK1TV9/8fyaIGz6AadnSji28NC+LK4fDcY+MN7HyckFbAuee7dAMwjpTzuM5BJcLjZoeAH37evnxw8P8fzU4ro/30pbBzMOxZqH6wWkHIuGF+I2A/YNppmvRFs1lOoefekEFnB5qb7mhPPV56bY25/CD1k9Dr04hVEDXalNWgm/ZhDw7JkCQWhwvQDHtuhLWzSDkGozIOvYVJs+Y/kUrhdScj2aQYgtgi1QrPvY4pJK2PSlnbbv+rNSw39oXx8zJZf5SoPlusdTxyfWhO7NepIroVRyPQbS0eMpx6K4sq9Jw9/0CtIbBwwPjub44o8fWtOW6UKFZ0/PIkQ943vGcrxXqLbmTwtyXS96puRSbQaEGIr1qBS0WHF5ebrKnr4U/9MDe6i6Pl89Mc1D+/po+gE/nI5q5Enb4neen+bevf1cLbos1z32DaRwg+i1+jMOd49kAcOpc8tkkzYDrdC/uFRjvC+9ui/P6xeXefNKEUG4eyTLYC7J6StF9vSlem59iwa6Upt0teji2PDK+SLVpsdQLsl82aUZhJRdb3UzaUsgNFBrBohAwrKwEMSCMISG57Nch+XLRfYNpOlPJ8imElsOguuD+vBoFJqluse5pTqP3dDujXqSqyWbVl0549gsVhos133+7I0rjORSTBcqm64p36xMsXKlkHJsxEDDD6N9b1q9Zi8ImZqvRL3ohEXCEu4ayvKB/QM0mj4vTi+xWG2SsIQ9fUn29EdTTsPQcHauTBBCplUjn6+6/NkbizxYqGDbFkEQMln1eOyeUT5+eJiK63POqfLeXJl8ysbzA84vVqk0fMIw5OxcGS8IOPHeAvOVBk0/xBKYnq+SLrrYVlTnHx/Irh7PXpjCqIGu7hjbnYc8PpDmxNkC2aRNPuUQBIZc2qHienihIZ+ysUQAwbGg5EblFi8ICUKDbQTbgqZvcPwQsYV606dQcvno3UM3Xeh2fdsTFhjg5Hvz7B/McO+ePKP51jz4dS75V0K7GQSrpYVkwuKhvX0cPzrCv/qv7/D2bIlCuUmrc8w9e/pIWBZZx+L//JO3ODSU5YH9/Rses40GCB89Msy5pTpXiy7nFqrsH0jzgf39q73idNJmrtwAoNIIsMQinbSp1D2ageHQUJbZYo0fnlsmn0rgBdGmY1PzNc7OFrl37wAj+SRvXily33j/td79sktool5+LuVQaXg0vJDXLy2TbG1U9suPH+FrJ6Y5O1fi7FwFSyxGckn8UHjzSonJ+SqphE3DD2l4ISHgetG4SCph8eyZAn/jI/tJJYSLi7XVAeWVGnw36uga6OqO0IkZCcePjvBHr1xitD/FWD7JVKFKOmHRsAQ/NDT9kJF8kuFskonRHCfem6dS9/ACQzJhkUvaVBt+VKLJONT9ABCOjuW4d8+1eveNJx1gTe3+B5OLCNCfTlCu+5w6t8yxw4OM5tdftHT86AhPn5hiqlBlMOPgBQHn5qtcWqxxbr7CO7PlaHuCbIKFmkcYGhq+jyXw4vRiq7cMuaTNP3tzlsMjWe5vDbpCdKXw7OlZUo7NB/b3Y0l01bBYbUTlkL1RSei1C8ucSVg8ft8oxw4P8l6hSrHWpC+VwGDoTzvUvYCq62FZwifvHqLmh7w4vUQ2aZNNJqg2A7wgJO0kOHVhmXv3DrCvP83rF5epuk2m55ssVJoU603SCWG+bBCEihtQdpvMV6K1AX/zo9GsotG+JH/xjkvCjs6SJdfHGEPDD/FDQ6bPXg1zAN+A74XUvRDHhu+/PYcXRM/tSztMjGbJJm2eaS3gut2hroGu7gidmJEwMZbnsftGOTNTJjQw0VrO/+blIrZY7OlLtnpxAadnymCEPf1pmmFIxQ2wbYnC0bZIOjZ/69gBRvNpQmOYLbobnnRSCVlt+4vTJYaySQzg+QGBMVginJ2rkLTtdS/5J8by7OlLsVBpslRrslhtMppLUm36fPd0gVTC4sF9fezpz3B6pkjT87m83GBPH+RT0fa6Z66UqLg+uZRD0fWouj5Pn5jCGDg8kkMQxLB6cgH4weQCc2WX2VKdff1pRnIOFdfnxHsLfP5D+3hkYmR1sdHbV6Lad6Xh0zea556xHMP5FC9NznN5qU7SFqrJKGzLrk/ShktLIecKFUoNn6FMgtcvFmkGAQnLwhioNg0Jy6dQJrpCsiwkMPzFu/M0Q8NPP7yXFycXuFp2STnRCcMLQqqN6GQWhHBxsYZZ53fBEP17oeziGzCtJ711pcxkocpQNnXb98IHDXR1h+jUjISf/fBdNPxrM0YuLFQ5v1CLAqcZ0J9JUK4HBKEhYUO5EZBN2jwwnmHZ9ZltuvSlEgxnk6tlkpVe9UYnnRem5nniob1AVCMfyDhgotkhx+4e4uxchZlll49PjKyG+ddfPL+mlx+E8Ph9Y/xwepGUbTFbbuDYQmBCErbF2bkK2WSCtGNTrDXxW73OjGPjBwbEYqHapNGqcwvCQtUll3T4sQOD9GccGn5INmnzlxeX8YKorSvjB7OlBqP5JHUvoOH7656AqjcMMJ4rVFh2fcb6kpRdn7LrU28GJBMWgTEQGl69uMR4f5qEZRGEhiAEP7jWo/ZCKLo+QjS2kUvaWLZwYaHKV09MM1dpMJCJxg+KtSahAcHgGXAswQ3Wi/NIYMAEYFlCYAzzVW/13+bLdRYqDb78n/+Sx+8dvW33XdVAV3eEdpaXr1cvvn4mx8vTC7xxuUQqIRTK0SV+0fXIOjbJhM1wPsm7s2UWqw3OL9TIOBYHh9KkkwkGcw6hMWsG0jba2hZkte39mWiRjyEag32vUGWh0mTvOuWZG3v5Fden5Hosux6ObSFA2rERhBCYqzTY05diej4K91TCot6MFgdZhCxVfRKWxUArvC8s1NnXHwXePWM5TpwtUGkGzC7XGe2L9m+xJJrz7YeGajNgvD/Nct1jZtnlyGhude54woLZUoP+VIKzc2XOL9ZYrDY5Mpbjwwf6+cHkEl4Q4tiCSBTcYoW8daXI6ZkSxhCNVbTyV2BNz9oQBXCpEVApVLm8UGUwlyI0hkzSIjTg+iEmNBggBIKbhPmKkGhQ9kZVD6rLdeqez8uOheub21KC0UBXd4R2lpdvVGNf+fq7b81yaDjLQMah0ljE9UPqDR/Hssgl4Y3LRRp+9GEXosCZnI+WrB8Zzb1vL5CNTjofOTiwuurz6GiWH0wu4jZ9AmO4uFTDbQaIGJ4+McV4X2rdXr7rRT9v0raouj7ppE0zMEyMZJmar0Xh3fDxs0lyyQT7BlIkEzb1ZsjegRRnXQ9jDEs1j3zKZqZYxxZDqdHaeVCiWT3z5QY1L6Dk+vQ5Fk1LqHkBmaRNpeExmEnyySMjZJI2Dd+sudJ54/Iys8t1vBDG8knyqQSL1WhBUNoR5iohxkThbAOBTWvlrXlfWeRmURwCjRBqDY/BXJLlmhfVwTcR4Fs1X/WYn1xkZtml3vD4v372gx1/j+tpoKs7wq02Vdpqjf3k5AJ+GDLcl0YQhnIp3GaAY1kYY5heqOMHZrWnaIgu/23AIMxXPf7hp9Zehm900nny2IHV96wUfT5xZJgXJhe4uBwtfjm0N4tj20wVqrx5aYnDo33Mld0oSBPRjI7BtMPfPT7Bt1+/zJmZEmHDQxAuVpsEQYAnUGn65NMJvvzpe5lcqLNYdlmqLfPG5SJV1ydhC4M5m76Mg9sMqHshCTukVPd47cISJdcjn3YYzDoIghcaBlMWNS+61Vt/2uGBvdFqU4HVvVWulurMFBuU602afnS7vPlq9HcvNFxZqhEYSCcsgjCkEUAABMH2fidcP8SxLLwgoOnf+vnbMb1QY658iU8/ON72zTU2QwNd3TFutsvgyqBcuRHNyV4ZlNuoxn616DKaT9HwQtKOzZ6+VHTbN0w0UyK8dvl/PSNQb/pUXH/dJeYr0/zWO+lc3/bXLv6QD+4fIJMPAegEAAAMCklEQVS89hF2LLi41CCTdFisNrFEcJsBGFhyoh7+rzxxPwcG0vzms2epNjxSCSGZiFaHph0LG3BD8D2P7565SmgM+ZSDCQ2+MYhAwwtIJ21SiTSuH/DC1Dw/urTMUGt2jwi8e7VMvRmwXDM8vLefmh/w8L5+Do5Et8L72onp1b1VLi3WmCs3WLkNajMICVs9cWNgJWvr3rX7pHaCF8K5hdrqCXcnGaDSDPmX3zmjga7UrWxnjvl0obK6j8nKANmp80s8MN7HwdHcut8zPpCm6Qe8PRvdUCGXSjCci2af1Jo+goUEwfuCIjRgW0LZ9Xj7Sul9S8yfn1rcsNY6Xajwx69f5rWLRX50aYnhTJI9/WmqXoDbDCi6UUAv1jwsETLJaD+VxVqTn7h/z+rOgs+eniXj2CRtodwIMBhCYLnu8xdnC5ydK/PO1UprgBCaflSzdyyh3gz5sYlBFisNphdqjPY5fPTQIG9eLjFfaeDYFsPZ5MoaKzAGIwZLoFhv8tJUtAnYy+cWqDQCqk2fYj2aVSJE5ZCNcnsnQrezp4hbe2NmZ2/AIca0f5hE5LPAV4iuJL9mjPn1mz3/2LFj5tSpU22/n1LrWW+vkpVSxWZuYfbs6Vn8IKQZGEZyqWjpe83DDw3/9xceBt6/URNEg49hGK4uu09YFp//4F5+78ULXFisUW/664ZT1oHhXIYPHxrkA/sH3rfgZ19fintaS9ZX3u/SYo2vfP8sV5ZdBjIOZTeafmhZFuN9Ub370lKNoWx0g+pyI1jtdY/lo2Xyz707z08+vJcXpxY4v1Cl3PDJOhaLNR8/CPACw8r43sqgK1wL2qQN6YTNxw4PM19pkEnYWBZcWW5QaK2YFYlOWLlkgkwywaGhDE7CxvV8FipNDg5nKZTrTBZqNFonijvNZx8Y4T/84iNb+h4RecUYc+xWz2v7FnQiYgP/D/BTwEPAF0XkoXZfT6l2tXsLs68+N8WJdwtMz1dZrDapNwO8INoyti+TYGIkWta93q3OAJ48doCDIznuHs7x0x/cz1PHJ5hcqPPw/j7C0JCw3v/xSgiknQQi0YZRDT9aGNT0QgYyDjXX48/fusrF+erq+331uSm+8v2zq/uZWCLYloUlQmgM5Ua0xcBQLsVgNslizSOfSnBwOEsu5VBtBrx8bpFMyuL0TInLy3XqXrRAp1j3afo+rm9WS0QrIXtj2DYDSCQsPnRggL6UQ2AM5xdqlBseaUcwJpqb3fQN9aZPEEZns2wy2hcdorCfKTUIwvBaL/4O8523N/693K7tlFw+DrxnjJkCEJHfB74AnO5Ew5TarHbmmH/79cucm68y0LqPaN0LcT0fQ4qffGjv6oKXmw2WfumRu9dcAazsInhoOMuFhTo/ulzEBpCoPVHghQQh/MT9Y9yzt391K4F0a8n6Qi2aZz5bbnB4LN9acRmtfkzYFk7CQhDyKYdya/vXjGNzz3gfI9kEz51dwLEtQmOoNwPAMJxNMlmocGAoS8MPOTiU4d2rPg0vxPdDrh9btIQ1tX8DUa+b6PGJ4RyvXVxmOJ9kudqkGRjSjkW9CekkhK2xA7GEe/fkmS012DOQptqIBlwLlcbqQLHcOLdQbdt2bhJ9F3Dxuq8vtR5T6rZq58a9r19cpr9Vnhhr3Q80YQsXFmprbmR8tei25oJfs9FNnq9/7ifuGWHfQIaDw2ny6WhOty0wmk9y73ieX3x0guNHR5gvNQgxGGOoewGVhs/BoXS02VdL0w+xrejepCtT6xKWYIkwmHX40KFBHpkY4d7xAfrSNsO5JCP5JKmkxVA2yfhgmlQi2n8m49j0Z5Lcv7efoaxDSJSpVmvO+HoV2JWpggNpm7/+sQN88K5BvvzEfa0bWJho4DIMsQUGMg65pE3KtsmkEhhjKNaipfz96Wh2TMKW1p43YG+hm24LJG0hm7Tu2N79rWynh77eMX3fr4OIPAU8BXDo0KFtvJ1S62vvxr3RUnWAvnTUq76yXKfa9Nfs3b2ZBUkrrn/uaD7NY/cO89/fnSefSjCUTeLYFvl0Ys32titbCZRcj760w7178oQha94vmbAY70/RDAzz5QaZpI0XhCQTQsKy2NuXWl2k1JdO8uB435pbs5XqHleL7up2vynHwrEt7tnTRyphUSg38MIQEYkW05hrg4UJCxzLQizh0w+Or/7sE2N5Hr93lGKtyXxrRk0qYWEQbBHuGc+DgaFcEj80fPTgIDPFBoZo98mVRUtJJ7qacL2bd9Wt1nHoTzvkUtGGXvVmsO5MojvZdgL9EnDwuq8PAFdufJIx5mngaYgGRbfxfkqtq50b937k4AAvTC0iIqsB159x+MzD42v2/97KyeLG5473Z3lkYuR9N0a4vl3rbSXww3NL3LcnvxrSKz3uvlSC9+bKTM/X8EPDpx8c54kHx9dMc3zq+ATPTy1SWtnPvNXeTxwdod6Mblyxsh/4oeE+BrIOYRjyg8lFLDEkbIumZ6j7Pik7mkven3Z49Ogw9473r/nZv/Dhu7habnD6SomyG5WFRIS7BtN86K5BLNtaM4c+lShhWzBXarC3L82FxSpl18cIZB2h4Ufz9m0rWjAUEgX5cN5hYiTPg/v6EODlc4vkkjZNL8CYm89USbS2LL7ds1m6pe1ZLiKSAN4FPgVcBl4G/hdjzFsbfY/OclG9YrpQ4ekTU9H+JF5AyrEZySXXvQ3cVqZEtjN98sbvOTyUWd1y9vpZNdtpA7DuTKBHjwzz/NQiS1WXt66UWao1STs2f++xCf7nj9+9qT3OV6ZS1ps+2aTNXYO33mr35OQC786UuLRco+JGVfylWoOGH1JxfcQS9g2k+dXPPPC+edvThQrffv0y33ljlgtLNYIgGpdYKR8BpBPg2BaphE3asSjXPIrN3on1c7/+M1t6/mZnuWx32uJPA/+WaMzkt40x/+Jmz9dAV71ku/uj7zYb/by9chw6cTK8/uT37kyJi8s1qq5PNuVweDjTKgHxvp//269fXnPT7rsG0vzZm7O8PVum1gwIQ4MXBDRaI8hRCUjoS0dz/avNzefoVsMcblOgb5UGulJKbd2Oz0NXSinVWzTQlVIqJjTQlVIqJjTQlVIqJjTQlVIqJm7rLBcRKQDnb8NbjQLzt+F9dgs9HtfosVhLj8c1vXws7jbGjN3qSbc10G8XETm1mSk+dwo9HtfosVhLj8c1cTgWWnJRSqmY0EBXSqmYiGugP93tBvQYPR7X6LFYS4/HNbv+WMSyhq6UUneiuPbQlVLqjhO7QBcRW0ReE5E/7XZbuk1EzonIGyLyuojc8buiicigiHxLRN4WkTMi8olut6kbROT+1u/Eyn8lEfnH3W5XN4nIr4jIWyLypoh8U0Q2vt1VD4tdyUVE/glwDOg3xnyu2+3pJhE5BxwzxvTq3NrbSkR+DzhhjPmaiCSBrDFmudvt6qbWzd4vAz9ujLkda0R6jojcBZwEHjLG1EXk/wP+izHmd7vbsq2LVQ9dRA4APwN8rdttUb1FRPqBx4HfAjDGNO/0MG/5FDB5p4b5dRJApnXjnizr3H1tN4hVoBPdbOOfcufccepWDPBdEXmldW/XO9kRoAD8Tqsk9zURyXW7UT3g54BvdrsR3WSMuQz8BnABmAGKxpjvdrdV7YlNoIvI54A5Y8wr3W5LD3nUGPNR4KeAvy8ij3e7QV2UAD4K/HtjzEeAKvCr3W1Sd7XKTp8H/nO329JNIjIEfAGYAPYDORH5+e62qj2xCXTgUeDzrbrx7wM/ISL/sbtN6i5jzJXWn3PAHwEf726LuuoScMkY81Lr628RBfyd7KeAV40xV7vdkC77NDBtjCkYYzzgD4FPdrlNbYlNoBtjfs0Yc8AYc5joMvL7xphdeZbtBBHJiUjfyt+BnwTe7G6ruscYMwtcFJH7Ww99CjjdxSb1gi9yh5dbWi4Aj4hIVkSE6HfjTJfb1JZEtxugdsw48EfR7ycJ4D8ZY77T3SZ13T8AvtEqNUwBf6fL7ekaEckCTwC/3O22dJsx5iUR+RbwKuADr7FLV43GbtqiUkrdqWJTclFKqTudBrpSSsWEBrpSSsWEBrpSSsWEBrpSSsWEBrpSSsWEBrpSSsWEBrpSSsXE/w/nEH6BhhQFhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# doesn't look useful...\n",
    "plt.scatter(np.sqrt(rf_error), np.abs(Y_test.ravel() - Y_pred_test), alpha=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conformal prediction for error estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc = NcFactory.create_nc(RandomForestRegressor(n_estimators=10, max_leaf_nodes=1000))\n",
    "model_icp = IcpRegressor(nc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10483, 127)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zhuangjw/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/nonconformist/base.py:60: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
      "  self.model.fit(x, y, **self.fit_params)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.08 s, sys: 31.3 ms, total: 4.11 s\n",
      "Wall time: 4.14 s\n"
     ]
    }
   ],
   "source": [
    "%time model_icp.fit(X_train[0:5000], Y_train[0:5000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.04 s, sys: 331 ms, total: 2.37 s\n",
      "Wall time: 1.95 s\n"
     ]
    }
   ],
   "source": [
    "%time model_icp.calibrate(X_train[5000:], Y_train[5000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 1503163 is out of bounds for axis 0 with size 5483",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m~/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/nonconformist/icp.py\u001b[0m in \u001b[0;36mpredict\u001b[0;34m(self, x, significance)\u001b[0m\n\u001b[1;32m    387\u001b[0m \t\t\t\tp = self.nc_function.predict(x[idx, :],\n\u001b[1;32m    388\u001b[0m                                                              \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcal_scores\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcondition\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 389\u001b[0;31m \t\t\t\t                             significance)\n\u001b[0m\u001b[1;32m    390\u001b[0m                                 \u001b[0;32mif\u001b[0m \u001b[0mn_significance\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    391\u001b[0m                                         \u001b[0mprediction\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0midx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/nonconformist/nc.py\u001b[0m in \u001b[0;36mpredict\u001b[0;34m(self, x, nc, significance)\u001b[0m\n\u001b[1;32m    501\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0msignificance\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    502\u001b[0m                         \u001b[0mintervals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 503\u001b[0;31m                         \u001b[0merr_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merr_func\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_inverse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignificance\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    504\u001b[0m                         \u001b[0merr_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0merr_dist\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mn_test\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    505\u001b[0m                         \u001b[0merr_dist\u001b[0m \u001b[0;34m*=\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/Computing/miniconda3/envs/sci/lib/python3.6/site-packages/nonconformist/nc.py\u001b[0m in \u001b[0;36mapply_inverse\u001b[0;34m(self, nc, significance)\u001b[0m\n\u001b[1;32m    163\u001b[0m                 \u001b[0;31m# TODO: should probably warn against too few calibration examples\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    164\u001b[0m                 \u001b[0mborder\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mborder\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 165\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mborder\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mborder\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    166\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    167\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1503163 is out of bounds for axis 0 with size 5483"
     ]
    }
   ],
   "source": [
    "%time y_pred_icp = model_icp.predict(X_test, significance=0.05)\n",
    "# requires a terribly huge \"calibration set\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}