robust_linear_regression_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is for unsderstanding Robust Linear Regression of statsmodels\n",
    "\n",
    "Reference : https://www.statsmodels.org/stable/generated/statsmodels.robust.robust_linear_model.RLM.html\n",
    "\n",
    "versions :\n",
    "\n",
    "    python 3.6.2\n",
    "    statsmodels 0.8.0\n",
    "    scikit-learn 0.19.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(action='ignore')\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.tools as st\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load data to use\n",
    "data = sm.datasets.stackloss.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = pd.DataFrame(data.exog, columns=data.exog_name)\n",
    "Y = pd.DataFrame(data.endog, columns=[data.endog_name])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AIRFLOW</th>\n",
       "      <th>WATERTEMP</th>\n",
       "      <th>ACIDCONC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>80.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>80.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>88.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>75.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>90.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>62.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>87.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>62.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>87.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>62.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>87.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>62.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>93.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>62.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>93.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>58.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>87.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>58.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>58.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>58.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>88.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>58.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>82.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>58.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>93.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>50.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>50.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>86.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>50.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>72.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>50.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>79.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>50.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>56.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>82.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>70.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>91.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    AIRFLOW  WATERTEMP  ACIDCONC\n",
       "0      80.0       27.0      89.0\n",
       "1      80.0       27.0      88.0\n",
       "2      75.0       25.0      90.0\n",
       "3      62.0       24.0      87.0\n",
       "4      62.0       22.0      87.0\n",
       "5      62.0       23.0      87.0\n",
       "6      62.0       24.0      93.0\n",
       "7      62.0       24.0      93.0\n",
       "8      58.0       23.0      87.0\n",
       "9      58.0       18.0      80.0\n",
       "10     58.0       18.0      89.0\n",
       "11     58.0       17.0      88.0\n",
       "12     58.0       18.0      82.0\n",
       "13     58.0       19.0      93.0\n",
       "14     50.0       18.0      89.0\n",
       "15     50.0       18.0      86.0\n",
       "16     50.0       19.0      72.0\n",
       "17     50.0       19.0      79.0\n",
       "18     50.0       20.0      80.0\n",
       "19     56.0       20.0      82.0\n",
       "20     70.0       20.0      91.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# view the content of data as DataFrame\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x10b294710>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEICAYAAAB/Dx7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEhVJREFUeJzt3XuQZGddxvHv4yZIzMRASJjCbHTR\nwohmNZghYEHpbLwFguAFS1KAiaLrDQt1QRerFPEa1HCR8hYREzUyRi5CJaKGhCGgEtwlwU0MCMKq\nCWHXEAhMKkYXfv7RZ804Oz19ZrI93S/7/VR1zXT3e04//VbPM2feOdOTqkKS1I7Pm3QASdL6WNyS\n1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3Nk2SJyf5+yT3JLk7yd8leUmSpe7yX0k+s+z6rcu2\nTZIPJ/nnIfv+1iQ3JPl0kv9M8o4kT+/uuzjJu5aN/cLusd+Q5Pgklyf55SH7TZIXJflgkvuS/HuS\nS5J8/rIxW7t93dU9t31JLl52//OSvL/LdiDJNUlOOgpTqmOUxa1NkeQLgauBVwOnAKcDLwXeVFUz\nVTUD/DDwD4evV9VXLdvF1wOPBL40yeNX7PuZwF8AfwxsBWaBnwe+bZUcDwfeBvwb8D1V9T8jov8W\nsBP4XuAk4CnAecBVy8b8CfAfwJcAj+jGHuge7xuAXwUurKqTgMeu2FZat+MmHUDHjC8HqKrXddfv\nA/52HdtfBLwZOKH7/B9hcEQMvBz4pap6zbLx7+gu/yfJqcC1wM3A86rqs2s9YJLHAD8KfF1Vvae7\n+dYk3wV8KMl5VXU98HjgJ6vq3m7MTct283gG34xuAqiqu4Er1vG8pSN4xK3N8i/AZ5JckeQp3ZFv\nL0m+AHgmcGV3eVaSh3R3nwmcAbx+xG5OYVDkNwLfP6q0O98I3L6stAGoqv8A3g18c3fTu4HfTvKs\nJF+8Yh83At+a5KVJnrR8iUXaKItbm6KqPgU8GSjgD4D/TPKWJLM9Nv9O4H4GR+hXM/hJ8YLuvkd0\nH+8csY8zGBz1/1H1f4OeU9fY753d/QDfDbwT+DngI0luPrycU1Xv7PJ/LXAN8PEkL0+ypWcG6QgW\ntzZNVd1WVRdX1VbgLOCLgFf22PQi4KqqOlRV9wNv7G4D+Hj38VEj9vE+4IXAW5M8rmfku9bY76O6\n+6mqT1TV7m5NfpbBUsxfdss4VNVbq+rbGBz1PwO4GPiBnhmkI1jcmoiqej9wOYMCHyrJVga/DHxO\nko8l+RiDZZOndmvWH2Dwi8Hv6vGYrwIuAa5Nsubjdq4Hzkhy7opMZwBPBK5b5THuAn6TwTelU1bc\n99mquq7bb5/Hl1ZlcWtTJPmKJLu6Ij5cfhcyWB9ey3MZrI+fCZzdXb4cuJ3BmRoF/BTwc0m+rzvV\n7/O6Uw8vW7mzqvp14FXA25KcueyuLUkeuuzykKr6F+D3gCuTPDHJliRfBbwBeFtVva17Li9LclaS\n47rT/H4E+FBVfTzJM7q174d3pxaeC3xDj+ctDWVxa7N8GngCcGOSexkU1y3ArhHbXQT8TlV9bPmF\nQaFeBFBVrwe+B/h+4KMMTsX7ZQZnoRyhqn4JeA1wXZIv627ezeBMl8OX67vbn9+N/VNgCfhrYJH/\nf4T/BcCbgE8CH2ZwWuDTu/s+Afwg8EHgU91+fqOqrhzxvKWh4j9SkKS2eMQtSY2xuCWpMRa3JDXG\n4pakxozlvUpOPfXU2rZt2zh2fdTce++9nHjiiZOO0UtLWaGtvC1lhbbytpQVJp937969d1XVaX3G\njqW4t23bxp49e8ax66NmcXGR+fn5ScfopaWs0FbelrJCW3lbygqTz5vk3/qOdalEkhpjcUtSYyxu\nSWqMxS1JjbG4JakxFrckNabX6YBJ9jN4d7fPAIeqam6coSRJw63nPO4d3ZvES5ImyKUSSWpMr/fj\nTvIRBm8IX8DvV9UR/1kkyU5gJ8Ds7Ow5CwsLGwq07457NrTdes2eAAfue+D69tNP3pTHXc2o57wy\n69Eyrue8tLTEzMzMWPZ9tLWUFdrK21JWmHzeHTt27O27DN23uL+oqj6a5JHAtcCPV9UNw8bPzc3V\nRv/kfdvuaza03Xrt2n6IS/c9sFK0/5IL1hg9XqOe88qsR8u4nvOk/3R4PVrKCm3lbSkrTD5vkt7F\n3WuppKo+2n08yOBfNJ279haSpHEZWdxJTuz+ASpJTgS+hcH/CpQkTUCfn79ngTclOTz+z6rqr8ea\nSpI01MjirqoPA1+zCVkkST14OqAkNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3\nJDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtS\nYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMb2LO8mW\nJDcluXqcgSRJa1vPEfcLgNvGFUSS1E+v4k6yFbgAeM1440iSRklVjR6UvB74NeAk4IVV9bRVxuwE\ndgLMzs6es7CwsKFA++64Z0PbrdfsCXDgvgeubz/95E153NWMes4rs067PnknOd/LLS0tMTMzM+kY\nvbWUt6WsMPm8O3bs2FtVc33GHjdqQJKnAQeram+S+WHjquoy4DKAubm5mp8fOnRNF+++ZkPbrdeu\n7Ye4dN8DT3//s+c35XFXM+o5r8w67frkneR8L7e4uMhGX6uT0FLelrJCW3n7LJU8CXh6kv3AAnBe\nkj8daypJ0lAji7uqXlxVW6tqG/As4Pqqes7Yk0mSVuV53JLUmHUtnFbVIrA4liSSpF484pakxljc\nktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1J\njbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQY\ni1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1ZmRxJ3lokvckeV+SW5O8dDOCSZJWd1yPMfcD51XV\nUpLjgXcleWtVvXvM2SRJqxhZ3FVVwFJ39fjuUuMMJUkartcad5ItSW4GDgLXVtWN440lSRomgwPq\nnoOThwFvAn68qm5Zcd9OYCfA7OzsOQsLCxsKtO+Oeza03XrNngAH7tuUh3rQWsoK/fJuP/3kzQkz\nwtLSEjMzM5OO0VtLeVvKCpPPu2PHjr1VNddn7LqKGyDJS4B7q+o3h42Zm5urPXv2rGu/h23bfc2G\ntluvXdsPcem+Pkv8k9dSVuiXd/8lF2xSmrUtLi4yPz8/6Ri9tZS3paww+bxJehd3n7NKTuuOtEly\nAvBNwPsfXERJ0kb1OYx7FHBFki0Miv6qqrp6vLEkScP0Oavkn4DHbUIWSVIP/uWkJDXG4pakxljc\nktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1J\njbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQY\ni1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1ZmRxJzkjyduT3Jbk1iQv2IxgkqTVHddjzCFgV1W9\nN8lJwN4k11bVP485myRpFSOPuKvqzqp6b/f5p4HbgNPHHUyStLpUVf/ByTbgBuCsqvrUivt2AjsB\nZmdnz1lYWNhQoH133LOh7dZr9gQ4cN+mPNSD1lJW6Jd3++knb06YEZaWlpiZmZl0jN5aynvw7nsm\n9rrdyOtr0nO7Y8eOvVU112ds7+JOMgO8A/iVqnrjWmPn5uZqz549vfa70rbd12xou/Xatf0Ql+7r\ns1I0eS1lhX55919ywSalWdvi4iLz8/OTjtFbS3lffeWbJ/a63cjra9Jzm6R3cfc6qyTJ8cAbgCtH\nlbYkabz6nFUS4A+B26rq5eOPJElaS58j7icBzwXOS3Jzd3nqmHNJkoYYuQBVVe8CsglZJEk9+JeT\nktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1J\njbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQY\ni1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMaMLO4kr01yMMktmxFIkrS2PkfclwPnjzmH\nJKmnkcVdVTcAd29CFklSD6mq0YOSbcDVVXXWGmN2AjsBZmdnz1lYWNhQoH133LOh7dZr9gQ4cN+m\nPNSD1lJW6Jd3++knb06YEZaWlpiZmRn74xyt1/VGXguTmuuDd9/zOfe6HeXBzPWOHTv2VtVcn7HH\nbfhRVqiqy4DLAObm5mp+fn5D+7l49zVHK9Kadm0/xKX7jtrTH6uWskK/vPufPb85YUZYXFxko6/V\n9Thar+uNvBYmNdevvvLNn3Ov21E2a649q0SSGmNxS1Jj+pwO+DrgH4Azk9ye5HnjjyVJGmbkgk5V\nXbgZQSRJ/bhUIkmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrck\nNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1Jj\nLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjelV3EnOT/KBJB9K\nsnvcoSRJw40s7iRbgN8GngJ8JXBhkq8cdzBJ0ur6HHGfC3yoqj5cVf8NLADPGG8sSdIwqaq1ByTP\nBM6vqh/orj8XeEJVPX/FuJ3Azu7qmcAHjn7co+pU4K5Jh+ippazQVt6WskJbeVvKCpPP+yVVdVqf\ngcf1GJNVbjui7avqMuCyPg86DZLsqaq5Sefoo6Ws0FbelrJCW3lbygpt5e2zVHI7cMay61uBj44n\njiRplD7F/Y/AY5I8OslDgGcBbxlvLEnSMCOXSqrqUJLnA38DbAFeW1W3jj3Z+DWzrENbWaGtvC1l\nhbbytpQVGso78peTkqTp4l9OSlJjLG5JaswxUdxJXpvkYJJblt12SpJrk3yw+/jwSWY8bEjWX0hy\nR5Kbu8tTJ5nxsCRnJHl7ktuS3JrkBd3t0zq3w/JO3fwmeWiS9yR5X5f1pd3tj05yYze3f96dMDBx\na+S9PMlHls3t2ZPOeliSLUluSnJ1d30q53Y1x0RxA5cD56+4bTdwXVU9Briuuz4NLufIrACvqKqz\nu8tfbXKmYQ4Bu6rqscATgR/r3g5hWud2WF6Yvvm9Hzivqr4GOBs4P8kTgZcxyPoY4BPA8yaYcblh\neQFetGxub55cxCO8ALht2fVpndsjHBPFXVU3AHevuPkZwBXd51cA376poYYYknUqVdWdVfXe7vNP\nM/giOJ3pndtheadODSx1V4/vLgWcB7y+u32a5nZY3qmUZCtwAfCa7nqY0rldzTFR3EPMVtWdMPiC\nBh454TyjPD/JP3VLKVOx9LBckm3A44AbaWBuV+SFKZzf7kf5m4GDwLXAvwKfrKpD3ZDbmaJvPCvz\nVtXhuf2Vbm5fkeTzJxhxuVcCPw18trv+CKZ4blc6lou7Jb8LfBmDH0HvBC6dbJz/L8kM8AbgJ6rq\nU5POM8oqeadyfqvqM1V1NoO/Vj4XeOxqwzY31XAr8yY5C3gx8BXA44FTgJ+ZYEQAkjwNOFhVe5ff\nvMrQqZnblY7l4j6Q5FEA3ceDE84zVFUd6L4oPgv8AYMv4qmQ5HgGJXhlVb2xu3lq53a1vNM8vwBV\n9UlgkcG6/MOSHP7Dual8+4llec/vlqeqqu4H/ojpmNsnAU9Psp/Bu52ex+AIfOrn9rBjubjfAlzU\nfX4R8OYJZlnT4RLsfAdwy7Cxm6lbF/xD4Laqevmyu6Zyboflncb5TXJakod1n58AfBODNfm3A8/s\nhk3T3K6W9/3LvoGHwZrxxOe2ql5cVVurahuDt/C4vqqezZTO7WqOib+cTPI6YJ7B2zYeAF4C/CVw\nFfDFwL8D311VE/+l4JCs8wx+jC9gP/BDh9eQJynJk4F3Avt4YK3wZxmsG0/j3A7LeyFTNr9JvprB\nL8i2MDjAuqqqfjHJlzI4SjwFuAl4Tnc0O1Fr5L0eOI3BUsTNwA8v+yXmxCWZB15YVU+b1rldzTFR\n3JL0ueRYXiqRpCZZ3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4Jakx/wsyHs4JhZ7CvAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x124890e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check histgram of STACKLOSS. There seems to be some possible outliers away from the main cluter. \n",
    "Y.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x126abd9e8>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGl5JREFUeJzt3X+U3XV95/HnKxOoDGMPIDhNkcyA\nWPYgPYITrYgLmYn2KKUFu/QoTmtU6nTXH8df7UpNu8K6UylLpevZxXYKQty9y8gCLh4OtmIy4ZcR\nN8EIwYgITmKQXyk/ZBgOQvLeP77fIXcm9873zp353l/f1+Oce+79vO/3x+edb3Lf+f76fBURmJlZ\ncS1rdgfMzKy5XAjMzArOhcDMrOBcCMzMCs6FwMys4FwIzMwKzoXAzKzgXAjMzArOhcDMrOCWN7sD\ntTjyyCOjv7+/rnmfe+45Dj300KXtUJM4l9bTKXmAc2lVi8ll69ateyLiqKzp2qIQ9Pf3s2XLlrrm\n3bRpE6tXr17aDjWJc2k9nZIHOJdWtZhcJO2sZTofGjIzKzgXAjOzgnMhMDMrOBcCM7OCcyEwMys4\nFwIzsxZUKkF/PwwNnUF/f9LOS1tcPmpmViSlEoyMwPQ0gNi5M2kDDA8v/fq8R2Bm1mLWrZspAvtN\nTyfxPLgQmJm1mF27FhZfLBcCM7MWs3LlwuKL5UJgZtZiRkehu3t2rLs7iefBhcDMrMUMD8PYGPT1\ngRT09SXtPE4UgwuBmVlLGh6GyUnYuPFWJifzKwLgQmBmVnguBGZmBedCYGZWcLkXAkldkn4g6aa0\nfaykuyQ9IOnrkg7Ouw9mZlZdI/YIPgHsKGv/LXBZRLwOeAo4vwF9MDOzKnItBJJeA/wecEXaFjAE\nXJdOsh44J88+mJnZ/PLeI/h74D8C+9L2q4CnI+KltL0bODrnPpiZ2TwUEfksWDoLODMiPiJpNfDn\nwAeBzRFxfDrNMcDNEfHbFeYfAUYAent7B8bHx+vqx9TUFD09PfUl0WKcS+vplDzAubSqxeQyODi4\nNSJWZU4YEbm8gC+S/I9/EngUmAZKwB5geTrNqcC/ZC1rYGAg6jUxMVH3vK3GubSeTskjwrm0qsXk\nAmyJGn6vczs0FBF/GRGviYh+4L3AxogYBiaAc9PJ1gI35tUHMzPL1oz7CD4LfFrST0nOGVzZhD6Y\nmVmqIU8oi4hNwKb080PAmxuxXjMzy+Y7i83MCs6FwMys4FwIzMwKzoXAzKzgXAjMzArOhcDMrOBc\nCMzMCs6FwMys4FwIzKxjlErQ3w9DQ2fQ35+021Ujc2nIncVmZnkrlWBkBKanAcTOnUkbYHi4mT1b\nuEbn4j0CM+sI69bN/HDuNz2dxNtNo3NxITCzjrBr18LirazRubgQmFlHWLlyYfFW1uhcXAjMrCOM\njkJ39+xYd3cSbzeNzsWFwMw6wvAwjI1BXx9IQV9f0m63E8XQ+FxcCMysYwwPw+QkbNx4K5OT7VkE\nZjQyFxcCM7OCcyEwMys4FwIzs4LLrRBIeoWk70v6oaT7JF2Uxq+W9DNJ29LXyXn1wczMsuU5xMQL\nwFBETEk6CLhD0rfS7/4iIq7Lcd1mZlaj3ApBRAQwlTYPSl+R1/rMzKw+uZ4jkNQlaRvwOHBLRNyV\nfjUq6R5Jl0n6tTz7YGZm81PyH/ecVyIdBnwD+Djwr8CjwMHAGPBgRPznCvOMACMAvb29A+Pj43Wt\ne2pqip6enjp73lqcS+vplDzAubSqxeQyODi4NSJWZU4YEQ15AZ8H/nxObDVwU9a8AwMDUa+JiYm6\n5201zqX1dEoeEc6lVS0mF2BL1PD7nOdVQ0elewJIOgR4O/BjSSvSmIBzgO159cHMzLLledXQCmC9\npC6ScxHXRsRNkjZKOgoQsA349zn2wczMMuR51dA9wCkV4kN5rdPMzBbOdxabmRWcC4GZWcG5EJiZ\nFZwLgZlZwbkQmJkVnAuBmVnBuRCYmRWcC4GZWcG5EJiZFZwLgZlZwbkQmJkVnAuBmVnBuRCYmRWc\nC4GZWcG5EJiZFZwLgZlZwbkQmJkVnAuBmVnBuRCYmRVcboVA0iskfV/SDyXdJ+miNH6spLskPSDp\n65IOzqsPZmaWLc89gheAoYh4A3Ay8E5JbwH+FrgsIl4HPAWcn2MfzMwsQ26FIBJTafOg9BXAEHBd\nGl8PnJNXH8zMLFuu5wgkdUnaBjwO3AI8CDwdES+lk+wGjs6zD2ZmNj9FRP4rkQ4DvgH8J+CqiDg+\njR8D3BwRv11hnhFgBKC3t3dgfHy8rnVPTU3R09NTb9dbinNpPZ2SBziXVrWYXAYHB7dGxKrMCSOi\nIS/g88BfAHuA5WnsVOBfsuYdGBiIek1MTNQ9b6txLq2nU/KIcC6tajG5AFuiht/nPK8aOirdE0DS\nIcDbgR3ABHBuOtla4Ma8+mBmZtmW57jsFcB6SV0k5yKujYibJP0IGJf0X4AfAFfm2AczM8uQWyGI\niHuAUyrEHwLenNd6zcxsYeY9NCTpTZJ+o6z9fkk3SvqypCPy756ZmeUt6xzBPwK/ApB0OnAx8DXg\nGWAs366ZmVkjZB0a6oqIJ9PP7wHGIuJ64Pr0/gAzM2tzWXsEXZJmisUaYGPZd3meaDarSakE/f2w\nbFnyXio1u0dm7Sfrx/wa4FZJe4DngdsBJB1PcnjIrGlKJRgZgenppL1zZ9IGGB5uXr/M2s28ewQR\nMQp8BrgaeFt6g8LMfB/Pt2tm81u3bn8RmDE9ncTNrHbz7hFI6ga2RsSLafsE4ExgZ0Tc0ID+mVW1\na9fC4mZWWdY5gn8G+uHlw0GbgeOAj0r6Yr5dM5vfypULi5tZZVmF4PCIeCD9vBa4JiI+DrwLOCvX\nnpllGB2F7u7Zse7uJG5mtcsqBOVDkw6RDCVNRPwK2JdXp8xqMTwMY2PQ1wdS8j425hPFZguVddXQ\nPZIuBR4Gjge+DS8PK23WdMPD/uE3W6ysPYIPkwwb3Q/8bkTMXKNxInBpjv0yM7MGmXePICKeBy6W\n9ArgeEmvBx6MiO8C321EB83MLF9Zg84tl3QJ8HOS5wv/L+Dnki6RdFAjOmhmZvnKOjT0X4EjgOMi\nYiAiTgFeCxyGDw2ZmXWErEJwFvDhiHh2JhARvwT+A8mNZWZm1uYyLx8tG1aiPLiX2ZeWmplZm8oq\nBD+S9P65QUl/DPw4ny6ZmVkjZd1H8FHgBkkfAraS7AW8CTgE+MOc+2ZmZg2Qdfnow8DvSBoCXg8I\n+FZEbGhE58zMLH81PVwmIjYy+6E0SNoVEVWH95J0DMljLX+DZDiKsYj4b5IuJLlR7Yl00s9FxM11\n9N3MzJbAYp4ypozvXwI+ExF3S3olsFXSLel3l0WELz81M2sBiykE8141FBGPAI+kn5+VtAM4ehHr\nMzOzHKjC1aH7v5Q+Xe0rYF1EHFHTSqR+4DbgJODTwAeAXwJbSPYanqowzwgwAtDb2zswPj5ey6oO\nMDU1RU9PT13zthrn0no6JQ9wLq1qMbkMDg5ujYhVWdNlFYLPzzdzRFyUuQKpB7gVGI2IGyT1kgxk\nF8AXgBUR8aH5lrFq1arYsmVL1qoq2rRpE6tXr65r3lbjXFpPp+QBzqVVLSYXSTUVgqxDQ1dGxO4q\nK/j9GjpxEHA9UJp5tGVEPFb2/T8BN2Utx8zM8pN1Q9mG9LDOLJI+CPz9fDNKEnAlsCMivlQWX1E2\n2buB7bV21szMll7WHsGngFsknTnzyEpJfwm8DzgjY97TgD8B7pW0LY19DjhP0skkh4YmgT+rs+9m\nZrYEsm4ou1nSC8C3JJ0D/CnJncWnVzrBO2feO6h8ianvGTAzayFZh4ZI7yL+ALAJOA5Yk1UEzMys\nfcy7RyDpWZJDOAJ+DVgDPJ4e/4+I+PX8u2hmZnnKOjT0ykZ1xMzMmiPrUZVvkvSuCvHflzSQX7fM\nzKxRanlU5Y4K8R3pd2Zm1uayCsGrImJybjAifgq8KpcemZlZQ2UVgkPm+e7QpeyImZk1R1Yh+I6k\n0fQqoZdJuog5zycwM7P2lHVn8WeAK4Cflt0d/AaSUUM/nGfHzMysMbIuH32OZEiI40geVQlwX0Q8\nlHvPzMysIWp9VOVDwEMAkl4r6a+A90bESXl2zszM8pc5xAQkI4ZK+qSk7wP3AV3Aebn2zMzMGiLr\nhrIPS9pI8mCZI0kGnXskIi6KiHsb0UEzM8tX1qGh/wFsBt4XEVsAJM37rGIzM2svWYXgN4E/Ar6U\nPmLyWuCg3HtlZmYNM++hoYjYExFfiYjTSUYefYZk9NEdkv6mIT00M7NcZZ0jeMvM54jYHRGXRsQA\ncA7wQt6dMzOz/GVdNXR5pWBE3B8RF+XQHzMza7CaLh81M7POlVUIjpP0zWqv+WaUdIykifR8wn2S\nPpHGj5B0i6QH0vfDlywbK5xSCfr7Ydmy5L1UanaPzNpP1lVDTwB/V+eyXwI+ExF3S3olsFXSLSTP\nP94QERdLugC4APhsneuwAiuVYGQEpqeT9s6dSRtgeLh5/TJrN1mFYCoibq1nwRHxCPBI+vlZSTuA\no4GzgdXpZOuBTbgQWB3WrdtfBGZMTydxFwKz2imi+v1h6V3F74uIR9P2+4F/B+wELoyIJ2taidQP\n3AacBOyKiMPKvnsqIg44PCRpBBgB6O3tHRgfH68xpdmmpqbo6empa95W41xmGxo6gwgdEJeCjRvr\n+v/LgnmbtCbnkhgcHNwaEasyJ4yIqi/gbuCI9PPpwC9ICsEXgOvmm7dsGT3AVuAP0/bTc75/KmsZ\nAwMDUa+JiYm65201zmW2vr4IOPDV17foRdfM26Q1OZcEsCVq+J3OOlm8LPb/r/89wFhEXB8Rfw0c\nn1VkJB0EXA+UIuKGNPyYpBXp9yuAxzOrlVkFo6PQ3T071t2dxM2sdlmFYLmkmfMIa5j9VLJ5zy+k\nTzW7EtgREV8q++qbwNr081rgxtq7a7bf8DCsXQtdXUm7qytp+/yA2cJkFYJrgFsl3Qg8D9wOIOl4\nkuEm5nMa8CfAkKRt6etM4GLgHZIeAN6Rts0WrFSC9eth796kvXdv0vYlpGYLk/WEslFJG4AVwLfT\nY06QFJCPZ8x7B3DgmbzEmoV21GwuXzVktjQyn1AWEd+rEPtJPt0xq92uXQuLm1llHmLC2tbKlQuL\nm1llLgTWtnzVkNnScCGwtjU8DKeeOjt26qk+P2C2UC4E1rY+8hHYsGF2bMOGJG5mtXMhsLY1Nraw\nuJlV5kJgbWvm/oFa42ZWmQuBta2ZO4prjZtZZS4E1rZmnj1Qa9zMKsu8ocysVV2ePlF7bCw5HNTV\nlRSByys+advMqnEhsLZ2+eX+4TdbLB8aMjMrOBcCM7OCcyEwMys4FwIzs4JzITAzKzgXAjOzgnMh\nMDMrOBcCM7OCy60QSPqqpMclbS+LXSjp4TkPszczsybKc4/gauCdFeKXRcTJ6evmHNdfOKUS9PfD\nsmXJe6nUWsuz1jSznYeGzvB2LqjchpiIiNsk9ee1fJutVErG2ZmeTto7d+4ffK2eJ3Yt9fKsNc3e\nzvJ2LqhmnCP4mKR70kNHhzdh/R1p3br9P9ozpqeTeCssz1qTt7MBKCLyW3iyR3BTRJyUtnuBPUAA\nXwBWRMSHqsw7AowA9Pb2DoyPj9fVh6mpKXp6euqat9XMl8vQ0BlE6IC4FGzceOuC17XUy5urU7ZL\nu+eR93ZulnbfLuUWk8vg4ODWiFiVOWFE5PYC+oHtC/1u7mtgYCDqNTExUfe8rWa+XPr6IuDAV19f\nfeta6uXN1Snbpd3zyHs7N0u7b5dyi8kF2BI1/MY29NCQpBVlzXcD26tNawszOgrd3bNj3d1JvN7l\nLZvzt2PZsvqXZ61pqf/eWHvK8/LRa4DNwAmSdks6H7hE0r2S7gEGgU/ltf6iGR5OHtDS1wdS8j42\nVv8JvzvvhH37Zsf27Uvi1jlm/72JRf+9sfaU51VD51UIX5nX+iz5x7tU/4DHxqrH/SCYzjLz92bT\npltZvXp1s7tjTeA7i62ivXsXFjez9uVCYBV1dS0sbmbty4XAKpq5qajWuJm1LxcCq+jyy2HNmtmx\nNWt8fsCsE7kQWEWlEmzePDu2ebPHoTHrRC4EVpGHHjArDhcCq2jXroXFzax9uRBYRStXLixuZu3L\nhcAq8tADZsXhQmAVDQ/D2rX77xvo6kraHnrArPO4EFhFpRKsX7//TuK9e5O2rxoy6zwuBFaRrxoy\nKw4XAqvIVw2ZFYcLgVXkq4bMisOFwCoaHYWDD54dO/hgXzVk1olcCKyquY+zzvHx1mbWRC4EVtG6\ndfDii7NjL77ok8VmnciFwCryyWKz4nAhsIp8stisOPJ8eP1XJT0uaXtZ7AhJt0h6IH0/PK/12+J4\niAmz4shzj+Bq4J1zYhcAGyLidcCGtG0taHg4eVB9Xx9IyfvYmIeYMOtEy/NacETcJql/TvhsYHX6\neT2wCfhsXn2wxRke9g+/WRE0+hxBb0Q8ApC+v7rB6zczszkUOV4cnu4R3BQRJ6XtpyPisLLvn4qI\niucJJI0AIwC9vb0D4+PjdfVhamqKnp6euuZtNc6l9XRKHuBcWtVichkcHNwaEasyJ4yI3F5AP7C9\nrH0/sCL9vAK4v5blDAwMRL0mJibqnrfVOJfW0yl5RDiXVrWYXIAtUcNvbKMPDX0TWJt+Xgvc2OD1\nm5nZHHlePnoNsBk4QdJuSecDFwPvkPQA8I60nYtSCfr7YWjoDPr7PY6+mVk1eV41dF6Vr9bktc4Z\npRKMjMyMpy927kza4KtgzMzm6sg7i/1QFTOz2nVkIfA4OWZmtevIQuBxcszMateRhWB0FJbNyWzZ\nMo+TY2ZWSUcWgquugn37Zsf27UviZmY2W0cWgg0bFhY3MyuyjiwEZmZWOxcCM7OC68hCsKbKLWvV\n4mZmRdaRheA734ETT5xpJaOrnnhiEjczs9k6shCUSjA5OdMSkLQ93pCZ2YE6shB4iAkzs9p1ZCHw\nEBNmZrXryELgISbMzGrXkYVgdBS6u2fHurs9xISZWSUdWQiGh2FsDPr6QAr6+pK2n0VgZnagjiwE\nkPzoT07Cxo23MjnpImBmVk3HFgIzM6uNC4GZWcG5EJiZFZwLgZlZwbkQmJkVnCKi2X3IJOkJYGed\nsx8J7FnC7jSTc2k9nZIHOJdWtZhc+iLiqKyJ2qIQLIakLRGxqtn9WArOpfV0Sh7gXFpVI3LxoSEz\ns4JzITAzK7giFIKxZndgCTmX1tMpeYBzaVW559Lx5wjMzGx+RdgjMDOzeXRcIZA0KeleSdskbUlj\nR0i6RdID6fvhze5nlip5XCjp4TS2TdKZze5nLSQdJuk6ST+WtEPSqe24TaBqLm23XSSdUNbfbZJ+\nKemT7bZd5smj7bYJgKRPSbpP0nZJ10h6haRjJd2VbpOvSzp4ydfbaYeGJE0CqyJiT1nsEuDJiLhY\n0gXA4RHx2Wb1sRZV8rgQmIqIS5vVr3pIWg/cHhFXpH+Ju4HP0WbbBKrm8knacLvMkNQFPAz8DvBR\n2nC7wAF5fJA22yaSjgbuAE6MiOclXQvcDJwJ3BAR45L+AfhhRHxlKdfdcXsEVZwNrE8/rwfOaWJf\nCkXSrwOnA1cCRMSvIuJp2nCbzJNLu1sDPBgRO2nD7VKmPI92tRw4RNJykv9kPAIMAdel3+eyTTqx\nEATwbUlbJY2ksd6IeAQgfX9103pXu0p5AHxM0j2Svtrqu+2p44AngKsk/UDSFZIOpT23SbVcoP22\nS7n3Atekn9txu8wozwPabJtExMPApcAukgLwDLAVeDoiXkon2w0cvdTr7sRCcFpEvBF4F/BRSac3\nu0N1qpTHV4DXAieT/EX5uyb2r1bLgTcCX4mIU4DngAua26W6VculHbcLAOnhrT8A/k+z+7IYFfJo\nu22SFquzgWOB3wQOJfn3P9eSH8/vuEIQEb9I3x8HvgG8GXhM0gqA9P3x5vWwNpXyiIjHImJvROwD\n/okkt1a3G9gdEXel7etIfkzbbptQJZc23S4z3gXcHRGPpe123C4wJ4823SZvB34WEU9ExIvADcBb\ngcPSQ0UArwF+sdQr7qhCIOlQSa+c+Qz8LrAd+CawNp1sLXBjc3pYm2p5zPwDTb2bJLeWFhGPAj+X\ndEIaWgP8iDbbJlA9l3bcLmXOY/bhlLbbLqlZebTpNtkFvEVStySx/9/KBHBuOk0u26SjrhqSdBzJ\n/54h2Y3/3xExKulVwLXASpI/7D+KiCeb1M1M8+TxP0l2dQOYBP5s5nhuK5N0MnAFcDDwEMkVHcto\no20yo0ouX6Y9t0s38HPguIh4Jo211b8VqJpHu/5buQh4D/AS8APgT0nOCYwDR6SxP46IF5Z0vZ1U\nCMzMbOE66tCQmZktnAuBmVnBuRCYmRWcC4GZWcG5EJiZFZwLgRWKpHdLCkn/Jm33S9qefl4t6Zl0\n+IgfS7q0bL4PSHqibDTLr6XxqyWdW2E9r5e0UdJP0lEj/1r77ZkZ8kDSirQ/byub94n0Mk6zhnAh\nsKI5j2SEx/dW+f72dPiIU4CzJJ1W9t3XI+Lk9PX+aiuQdAjJjVkXR8RvAW8guUP0I5Fcr30XcGo6\n+VtJrg1/azrvCcCeiPjXujM0WyAXAisMST3AacD5VC8EAETE88A26hvg633AnRHx7XRZ08DH2D/G\n0p2kP/zp+5eYXRi+W8c6zermQmBFcg7wzxHxE+BJSW+sNmF66OZ1wG1l4feUHRr64DzreT3JqJEv\ni4gHgZ50KOvvsr8QvBn4v8AxafutJIXCrGFcCKxIziO5VZ/0/bwK0/xbSfcAjwI3peMLzSg/NHTV\nPOsR1UeIDOD7wCnpOFIHRcQU8JCk4/EegTXB8uxJzNpfevJ1CDhJUgBdJD/Kl8+Z9PaIOEvSbwF3\nSPpGRGxb4OruI3mATfn6jyN5YtazafunwIeAu9NJvkfyJKpXA/cvcH1mi+I9AiuKc4GvRURfRPRH\nxDHAz0iG9T1Aevjoi0A9j2ksAW+T9HZ4+eTxl4FLyqa5k+QRl5vT9mbgE8D3wgOAWYO5EFhRnMf+\nEV1nXE/y7ORq/gE4XdKxGcv+R0m709fm9ETz2cBfSbofuBf4f8B/L5vnTpInns0UgrtJipIPC1nD\nefRRM7OC8x6BmVnBuRCYmRWcC4GZWcG5EJiZFZwLgZlZwbkQmJkVnAuBmVnBuRCYmRXc/wdsh6bS\nYziJpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x124986390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGnhJREFUeJzt3X+cXXV95/HXOzOohMEGBMcUyIyA\ndV1xGzaj62olM4PuWoqKXX0Izta4uk5LXOrP+ivdIttNi0rB9vFoaKdiiY+dByPLj8XygNVIuAEV\ncCcQfhldlCYRiUQaQIZQ3Ayf/eOcSybDzJyZyT3n3HvP+/l4zOPe7/ece8/nmzu5nznf7/d8jyIC\nMzOrriVlB2BmZuVyIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOziuss\nO4D5OOaYY6K3t7fsMBbsqaee4ogjjig7jMJUrb3gNldFq7Z569atj0bEsVn7tUQi6O3tZXx8vOww\nFqxWq9Hf3192GIWpWnvBba6KVm2zpJ3z2c9dQ2ZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBm\n1oRGR6G3F5YsSR5HR/M7VktMHzUzq5LRURgehn37kvLOnUkZYGio8cfzGYGZWZNZt+5AEqjbty+p\nz4MTgZlZk9m1a2H1h8qJwMysyaxYsbD6Q+VEYGbWZNavh6VLD65bujSpz4MTgZlZkxkagpER6OkB\nKXkcGclnoBg8a8jMrCkNDeX3xT+dzwjMzCrOicDMrOKcCMzMKi73RCCpQ9Jdkq5Pyy+XdIekByR9\nXdIL8o7BzMxmV8QZwUeA7VPKXwAuiYhXAI8BHywgBjMzm0WuiUDS8cDvAF9JywIGgavSXTYCZ+UZ\ng5mZzS3vM4IvA58Cnk3LLwEej4j9afkh4LicYzAzsznkdh2BpDOBPRGxVVJ/vXqGXWOW1w8DwwDd\n3d3UarU8wszVxMRES8a9WFVrL7jNVdHubc7zgrI3Am+XdAbwIuDFJGcIyyR1pmcFxwMPz/TiiBgB\nRgD6+vqiv78/x1DzUavVaMW4F6tq7QW3uSravc25dQ1FxGcj4viI6AXOBjZHxBBwM/CudLc1wHV5\nxWBmZtnKuI7g08DHJf2YZMzgshJiMDOzVCFrDUVEDailzx8EXlfEcc3MLJuvLDYzqzgnAjOzinMi\nMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMLOmNzoKvb2wZEnyODpadkT5K7LNhVxZ\nbGa2WKOjMDwM+/Yl5Z07kzLA0FB5ceWp6Db7jMDMmtq6dQe+EOv27Uvq21XRbXYiMLOmtmvXwurb\nQdFtdiIws6a2YsXC6ttB0W12IjCzprZ+PSxdenDd0qVJfbsqus1OBGbW1IaGYGQEenpASh5HRtp3\noBiKb7NnDZlZ0xsaau8v/pkU2WafEZiZVZwTgZlZxTkRmJlVXG6JQNKLJH1f0t2S7pd0QVp/uaR/\nlLQt/VmZVwxmZpYtz8HiZ4DBiJiQdBjwHUk3ptv+KCKuyvHYZmY2T7klgogIYCItHpb+RF7HMzOz\nxcl1jEBSh6RtwB5gU0TckW5aL+keSZdIemGeMZiZ2dyU/OGe80GkZcC1wHnAPwE/B14AjAA/iYj/\nNsNrhoFhgO7u7lVjY2O5x9loExMTdHV1lR1GYarWXnCbq6JV2zwwMLA1Ivqy9iskEQBIOh94KiIu\nmlLXD3wyIs6c67V9fX0xPj6ec4SNV6vV6O/vLzuMwlStveA2V0WrtlnSvBJBnrOGjk3PBJB0OPBm\n4IeSlqd1As4C7ssrBjMzy5bnrKHlwEZJHSQJ58qIuF7SZknHAgK2AX+QYwxmZpYhz1lD9wCnzlA/\nmNcxzcxs4XxlsZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnF\nORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxeWW\nCCS9SNL3Jd0t6X5JF6T1L5d0h6QHJH1d0gvyisHMzLLleUbwDDAYEb8JrATeKun1wBeASyLiFcBj\nwAdzjMHMzDLklggiMZEWD0t/AhgErkrrNwJn5RWDmZlly3WMQFKHpG3AHmAT8BPg8YjYn+7yEHBc\nnjGYmdncFBH5H0RaBlwL/Anw9xFxclp/AnBDRLxmhtcMA8MA3d3dq8bGxnKPs9EmJibo6uoqO4zC\nVK294DZXRau2eWBgYGtE9GXt11lEMBHxuKQa8HpgmaTO9KzgeODhWV4zAowA9PX1RX9/fxGhNlSt\nVqMV416sqrUX3OaqaPc25zlr6Nj0TABJhwNvBrYDNwPvSndbA1yXVwxmZpYtzzOC5cBGSR0kCefK\niLhe0g+AMUn/HbgLuCzHGMzMLENuiSAi7gFOnaH+QeB1eR3XzMwWZs6uIUmvlfSyKeX3SbpO0l9J\nOjr/8MzMLG9ZYwR/C/wKQNJpwIXA14AnSAdyzcystWV1DXVExN70+XuAkYi4Grg6vT7AzMxaXNYZ\nQYekerI4Hdg8ZVshU0/NrDmMjkJvLwwOrqa3Nylbe8j6Mr8C2CLpUeBp4FYASSeTdA+ZWQWMjsLw\nMOzbByB27kzKAENDZUZmjTDnGUFErAc+AVwO/FYcuAx5CXBevqGZWbNYt66eBA7Yty+pt9Y35xmB\npKXA1oj4f2n5lcAZwM6IuKaA+MysCezatbB6ay1ZYwT/G+iF57qDbgNOBD4s6c/zDc3MmsWKFQur\nt9aSlQiOiogH0udrgCsi4jzgt4Ezc43MzJrG+vWwdOnBdUuXJvXW+rISwdSlSQdJlpImIn4FPJtX\nUGbWXIaGYGQEenpACnp6krIHittD1qyheyRdBPwMOBn4Fjy3rLSZVcjQUPJTq21p65U4qyjrjOBD\nwKMk4wT/LiLq8wb+JXBRjnGZmVlB5jwjiIingQslvQg4WdKrgZ9ExPeA7xURoJmZ5Str0blOSV8E\nfkpyf+H/AfxU0hclHVZEgGZmlq+srqEvAUcDJ0bEqog4FTgJWIa7hszM2kJWIjgT+FBEPFmviIhf\nAueSXFhmZmYtLnP6aMxwd/uImOTgqaVmZtaishLBDyS9b3qlpP8I/DCfkMzMrEhZ1xF8GLhG0geA\nrSRnAa8FDgd+N+fYzMysAFnTR38G/BtJg8CrAQE3RsRNRQRnZmb5m9fNZSJiMwfflAZJuyJi1iWn\nJJ1AclvLl5EsRzESEX8p6fMkF6r9It31cxFxwyJiNzOzBjiUu4wpY/t+4BMRcaekI4Gtkjal2y6J\nCE8/NTNrAoeSCOacNRQRu4Hd6fMnJW0HjjuE45mZWQ40w+zQAxulj8+2CVgXEUfP6yBSL3ALcArw\nceD9wC+BcZKzhsdmeM0wMAzQ3d29amxsbD6HaioTExN0dXWVHUZhqtZecJurolXbPDAwsDUi+rL2\ny0oE58/14oi4IPMAUhewBVgfEddI6iZZyC6APwWWR8QH5nqPvr6+GB8fzzpU06nVapVapbFq7QW3\nuSpatc2S5pUIsrqGLouIh2Y5wNvmEcRhwNXAaP3WlhHxyJTtfwdcn/U+ZmaWn6wLym5Ku3UOIuk/\nAV+e64WSBFwGbI+Ii6fUL5+y2zuB++YbrJmZNV7WGcHHgE2SzqjfslLSZ4H3AqszXvtG4PeAeyVt\nS+s+B5wjaSVJ19AO4PcXGbuZmTVA1gVlN0h6BrhR0lnAfya5svi0mQZ4p732O8w8xdTXDJiZNZGs\nriHSq4jfD9SAE4HTs5KAmZm1jjnPCCQ9SdKFI+CFwOnAnrT/PyLixfmHaGZmecrqGjqyqEDMzKwc\nWbeqfK2k356h/m2SVuUXlpmZFWU+t6rcPkP99nSbmZm1uKxE8JKI2DG9MiJ+DLwkl4jMzKxQWYng\n8Dm2HdHIQMzMrBxZieDbktans4SeI+kCpt2fwMzMWlPWlcWfAL4C/HjK1cG/SbJq6IfyDMzMzIqR\nNX30KZIlIU4kuVUlwP0R8WDukZmZWSHme6vKB4EHASSdJOmPgbMj4pQ8gzMzs/xlLjEByYqhkj4q\n6fvA/UAHcE6ukZmZWSGyLij7kKTNJDeWOYZk0bndEXFBRNxbRIBmZpavrK6hvwZuA94bEeMAkua8\nV7GZmbWWrETw68C7gYvTW0xeCRyWe1RmZlaYObuGIuLRiLg0Ik4jWXn0CZLVR7dL+rNCIjQzs1xl\njRG8vv48Ih6KiIsiYhVwFvBM3sGZmVn+smYNbZipMiJ+FBEX5BCPmZkVbF7TR83MrH1lJYITJX1j\ntp+5XijpBEk3p+MJ90v6SFp/tKRNkh5IH49qWGvMLDejo9DbC4ODq+ntTcrWHrJmDf0C+ItFvvd+\n4BMRcaekI4GtkjaR3P/4poi4UNJngM8An17kMcysAKOjMDwM+/YBiJ07kzLA0FCZkVkjZCWCiYjY\nspg3jojdwO70+ZOStgPHAe8A+tPdNgI1nAjMmtq6dfUkcMC+fUm9E0HrU8Ts14elVxW/NyJ+npbf\nB/wHYCfw+YjYO6+DSL3ALcApwK6IWDZl22MR8bzuIUnDwDBAd3f3qrGxsXk2qXlMTEzQ1dVVdhiF\nqVp7oTptHhxcTYSeVy8Fmzcv6m/FltKqn/PAwMDWiOjL2i8rEdwJvDki9ko6DRgDzgNWAq+KiHdl\nHkDqIlmiYn1EXCPp8fkkgqn6+vpifHw861BNp1ar0d/fX3YYhalae6E6be7thZ07n1/f0wM7dhQd\nTfFa9XOWNK9EkDVYvGTKX/3vAUYi4uqI+K/AyfMI4jDgamA0Iq5Jqx+RtDzdvhzYk/U+Zlau9eth\n6dKD65YuTeqt9WUlgk5J9XGE0zn4rmRzji+kdzW7DNgeERdP2fQNYE36fA1w3fzDNbMyDA3BmjXQ\n0QEQdHQkZY8PtIesRHAFsEXSdcDTwK0Akk4mWW5iLm8Efg8YlLQt/TkDuBB4i6QHgLekZTNrYqOj\nsHEjTE4CiMnJpOwppO0h6w5l6yXdBCwHvhUHBhSWkIwVzPXa7wDPH11KnL7QQM2sPJ411N4y71AW\nEbfPUPd/8wnHzJrRrl0Lq7fW4iUmzCzTihULq7fW4kRgZpk8a6i9ORGYLUKZ6+7Uj71kCYUde2gI\nRkaS6wakoKcnKXt8oD1kjhGY2cHKXHfn4GNT6LGHhpKfWm1LS15cZbPzGYHZAs01g6adj23ty4nA\nbIHKnEHj2TuWBycCswUqcwaNZ+9YHpwIzBaozBk0nr1jeXAiMFugMmfQHHxsPHvHGsKzhswWocwZ\nNPVjmzWKzwjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKzinAjMzCrOicDMrOJySwSSvipp\nj6T7ptR9XtLPpt3M3szMSpTnGcHlwFtnqL8kIlamPzfkeHwryNq10NkJAwOr6exMypafMm5MY+0t\ntyUmIuIWSb15vb81h7Vr4dJL6yUxOXmgvGFDWVG1rzJvTGPtq4wxgv8i6Z606+ioEo5vDTQysrB6\nOzS+MY3lQRGR35snZwTXR8QpabkbeBQI4E+B5RHxgVleOwwMA3R3d68aGxvLLc68TExM0NXVVXYY\nuRoYWA1ohi3BzTdvKTqcwhX9GQ8Oribi+f/eUrB5czH/3lX4vZ6uVds8MDCwNSL6svYrNBHMd9t0\nfX19MT4+3ujwcler1dr+3q6dnTA5+fz6jg7Yv7/4eIpW9Gfc25t0B03X0wM7dhQTQxV+r6dr1TZL\nmlciKLRrSNLyKcV3AvfNtm8rq9Lgab1/er71dmh8YxrLQ26DxZKuAPqBYyQ9BJwP9EtaSdI1tAP4\n/byOX5aqDZ7W2zQyApOTQUeHGB5uz7Y2g/qA8Lp1yX2KV6xIkoAHiu1Q5Dlr6JwZqi/L63jNYq7B\n03b9ctywIfkp4yYtVeQb01ij+criBpupv3yuejOzsjkRNFhHx8LqzczK5kTQYB48NbNW07aJoKzL\n8DdsgHPPrZ8BBB0dSbldxwfMrPW1ZSKoX4a/cydEHLgMv8hksH8/3HzzFvbvdxIws+bWlonAl+Gb\nmc1fWyaCXbsWVm9mVmVtmQhWrFhYvZlZlbVlIvBl+GZm89eWiWBoKLmSt6cHpORxZKT9r8b0DUvM\nbDFyW2KibFW7DN83LDGzxWrLM4Iq8kwpM1ssJ4I24ZlSZrZYTgRtwjOlzGyxnAhyUB+0HRxcXdig\nrWdKmdliORE02MHLW6iw5S2qOlPKzA6dE0GDlTloOzSU3Lf22WeTRycBM5sPJ4IG86CtmbUaJ4IG\n86CtmbWa3BKBpK9K2iPpvil1R0vaJOmB9PGovI5fFg/amlmryfOM4HLgrdPqPgPcFBGvAG5Ky23l\n4EHb8KCtmTW93BJBRNwC7J1W/Q5gY/p8I3BWXscvU33QdvPmLR60NbOmV/QYQXdE7AZIH19a8PHN\nzGwaRUR+by71AtdHxClp+fGIWDZl+2MRMeM4gaRhYBigu7t71djYWG5x5mViYoKurq6ywyhM1doL\nbnNVtGqbBwYGtkZEX9Z+Ra8++oik5RGxW9JyYM9sO0bECDAC0NfXF/39/QWF2Di1Wo1WjHuxqtZe\ncJurot3bXHTX0DeANenzNcB1BR/fzMymyXP66BXAbcArJT0k6YPAhcBbJD0AvCUtW4srY20lM2uc\n3LqGIuKcWTadntcxrXgH3xBHviGOWQvylcV2SHxDHLPW50Rgh8RrK5m1PicCOyReW8ms9TkR5KCs\nwdO1a6GzM7kfQWdnUs6b11Yya31OBA1W1o1p1q6FSy+FycmkPDmZlPNOBl5byaz1ORE0WFmDpyMj\nC6tvJK+tZNbanAgarKzB0/qZwHzrzczqnAgarKzB046OhdWbmdU5ETRYWYOn9Yu45ltvZlbnRNBg\nQ0OwZk39L/GgoyMp591vvmEDnHvugTOAjo6kvGFDvsc1s9bnRNBgo6OwcWO9b15MTiblIqaQbtgA\n+/dDRPLoJGBm8+FE0GBecsHMWo0TQYN5yQUzazVOBA3mJRfMrNU4ETSYl1wws1bjRNBgXnLBzFqN\nE0EOvOSCmbUSJwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKU0SUHUMmSb8AdpYdxyIcAzxadhAF\nqlp7wW2uilZtc09EHJu1U0skglYlaTwi+sqOoyhVay+4zVXR7m1215CZWcU5EZiZVZwTQb4KuHV8\nU6lae8Ftroq2brPHCMzMKs5nBGZmFedE0ACSvippj6T7ptStlHS7pG2SxiW9rswYG03SCZJulrRd\n0v2SPpLWHy1pk6QH0sejyo61UeZo85ck/VDSPZKulbSs7FgbZbY2T9n+SUkh6ZiyYmykudor6TxJ\nP0rrv1hmnI3mrqEGkHQaMAF8LSJOSeu+BVwSETdKOgP4VET0lxhmQ0laDiyPiDslHQlsBc4C3g/s\njYgLJX0GOCoiPl1iqA0zR5uPBzZHxH5JXwBo9zZHxA8knQB8BfgXwKqIaMV59geZ4zPuBtYBvxMR\nz0h6aUTsKTPWRvIZQQNExC3A3unVwIvT578GPFxoUDmLiN0RcWf6/ElgO3Ac8A5gY7rbRpL/RG1h\ntjZHxLciYn+62+0kiaEtzPE5A1wCfIrkd70tzNHec4ELI+KZdFvbJAFwIsjTR4EvSfopcBHw2ZLj\nyY2kXuBU4A6gOyJ2Q/KfCnhpeZHlZ1qbp/oAcGPR8RRhapslvR34WUTcXWpQOZr2Gf8G8CZJd0ja\nIum1ZcbWaE4E+TkX+FhEnAB8DLis5HhyIakLuBr4aET8sux4ijBbmyWtA/YDo2XFlpepbSZp4zrg\nT0oNKkczfMadwFHA64E/Aq6UpBJDbCgngvysAa5Jn/9PoK0GiwEkHUbyn2U0IuptfSTtZ633t7bV\nKfQsbUbSGuBMYCjabOBthjafBLwcuFvSDpKusDslvay8KBtnls/4IeCaSHwfeJZk/aG24ESQn4eB\n1enzQeCBEmNpuPSvocuA7RFx8ZRN3yBJgqSP1xUdW15ma7OktwKfBt4eEfvKii8PM7U5Iu6NiJdG\nRG9E9JJ8Sf7riPh5iaE2xBy/1/+L5P8xkn4DeAGtuQjdjDxrqAEkXQH0k/yF8AhwPvAj4C9JTin/\nGVgbEVvLirHRJP0WcCtwL8lfRwCfI+lPvRJYAewC3h0R0wfSW9Icbf4r4IXAP6V1t0fEHxQfYePN\n1uaIuGHKPjuAvjaZNTTbZ/xt4KvASuBXwCcjYnMpQebAicDMrOLcNWRmVnFOBGZmFedEYGZWcU4E\nZmYV50RgZlZxTgTW8iRdIumjU8rflPSVKeW/kPTx9PnHJP2zpF9Ly/8+XSF2m6SJdHXJbZK+Jqlf\n0hNTtm+T9Ob0dZNp+T5J/yBpmaTXTNlvr6R/TJ9/W1KvpKenvdf70vfaIenWaW3apnQ12ylx3JWu\ninl+/v+qViWdZQdg1gDfA94NfFnSEpLrOV48ZfsbSJZGADgH+D/AO4HLI+KbwDcBJNVI5oePp+V+\n4NaIOHOGYz4dESvT/TYCH46I9STzzJF0OXB9RFyVlnuBn9RfM4MjJZ0QET+V9KoZtt8aEWdKOgLY\nJun6drouxcrlMwJrB98l+bIHeDVwH/CkpKMkvRB4FXCXpJOALuCPSRJCo9zGgRU5F+tK4D3p83OA\nK2baKSKeIlka+aRDPJ7Zc5wIrOVFxMPAfkkrSBLCbSRXOP9boA+4JyJ+xYEv2FuBV0qaz8qob5rW\nnXPQF7CkDuB0kqU1spw07b3eNGXbVcDvps/fBvzDTG8g6SUkC5/dP4/jmc2Lu4asXdTPCt4AXEzy\nF/obgCdIuo4AzgbeGRHPSrqGpDvprzPed7auocMlbQN6Sf5C3zSPGOfqGtoLPCbpbJI18KevWfQm\nSXeRLHtwYUQ4EVjD+IzA2sX3SL74X0PSNXQ7yRnBG4DvSvpXwCuATenaOGdzaN1D9TGCHpIFyD58\nCO9V93WSxDRTt9CtEXFqRKyKiL9pwLHMnuNEYO3iuyTLQO+NiMl0obtlJMngNpIv/c/XV8yMiF8H\njpPUcygHjYgngD8EPpkuX3worgW+SDp4bVYUJwJrF/eSzBa6fVrdE+mqmGeTfNFOdW1aP5fpYwTv\nmr5DRNwF3D2P95o+RvCH097nyYj4QjqeYVYYrz5qZlZxPiMwM6s4JwIzs4pzIjAzqzgnAjOzinMi\nMDOrOCcCM7OKcyIwM6s4JwIzs4r7/+KEfs8Fjz7XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x126a2ce80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGf5JREFUeJzt3X+UJWV95/H3Z2ZAaAbDT9sJON0i\nRonsCTit8ege6GkkRwiJYPAotmFUQseQsP7MKs7uCWy2EzQoZs9RdlsxjLEPIyu4sASNE2buqCdG\nt4dfAuMygMwITvixCNq0QRm++0fVde70dN/qX/Xj3vq8zrnndj236tb3Pl19v13P89RTigjMzKy+\nlpUdgJmZlcuJwMys5pwIzMxqzonAzKzmnAjMzGrOicDMrOacCMzMas6JwMys5pwIzMxqbkXZAczF\nUUcdFf39/WWHUYhnnnmGQw45pOwwKsv1k8111F6d6mfbtm1PRMTRWet1RCLo7+9nYmKi7DAK0Wg0\nGBwcLDuMynL9ZHMdtVen+pG0cy7ruWnIzKzmnAjMzGrOicDMrOacCMzMas6JwMys5pwIzKyjjI9D\nfz8sW5Y8j4+XHVHn64jho2ZmkHzpj4zA1FSyvHNnsgwwPFxeXJ3OZwRm1jHWr9+bBJqmppJyWzgn\nAjPrGLt2za/c5saJwMw6xurV8yu3uXEiMLOOMToKPT37lvX0JOW2cE4EZtYxhodhbAz6+kBKnsfG\n3FG8WB41ZGYdZXjYX/xLzWcEZmY150RgZlZzTgRmZjWXeyKQtFzS7ZJuTpdfKum7knZI+rKkA/OO\nwczMZlfEGcH7gO0tyx8HroyIlwM/AS4oIAYzM5tFrolA0rHA7wKfT5cFDAFfSVfZAJydZwxmZtZe\n3mcEnwb+I/B8unwk8FREPJcuPwwck3MMZmbWRm7XEUg6C3gsIrZJGmwWz7BqzLL9CDAC0NvbS6PR\nyCPMypmcnKzNZ10I108211F7rp/95XlB2RuA35d0JnAQ8EKSM4TDJK1IzwqOBX4808YRMQaMAQwM\nDMTg4GCOoVZHo9GgLp91IVw/2VxH7bl+9pdb01BEXBIRx0ZEP/B2YHNEDANbgHPT1dYBN+YVg5mZ\nZSvjOoKPAB+UdD9Jn8HVJcRgZmapQuYaiogG0Eh/fhB4bRH7NTOzbL6y2Mys5pwIzMxqzonAzKzm\nnAjMzGrOicDMrOacCMzMas6JwMys5pwIzMxqzonAzAo3Pg79/bBsWfI8Pl7cPoeGTi1sn52ikCuL\nzcyaxsdhZASmppLlnTuTZYDh4SL2qUL22Ul8RmBmhVq/fm8SaJqaSsq7aZ+dxInAzAq1a9f8yjt1\nn53EicDMCrV69fzKO3WfncSJwMwKNToKPT37lvX0JOXdtM9O4kRgZoUaHoaxMejrAyl5HhvLt9N2\n331GIfvsJB41ZGaFGx4u/ku4uc9GY6tvVTmNzwjMzGrOicDMrOacCMzMai63RCDpIEnfk3SnpHsk\nXZaWXyPph5LuSB8n5RWDmZlly7Oz+FlgKCImJR0AfFvS19LX/jwivpLjvs3MbI5ySwQREcBkunhA\n+oi89mdmZguTax+BpOWS7gAeAzZFxHfTl0Yl3SXpSkkvyDMGMzNrT8k/7jnvRDoM+CpwMfD/gH8F\nDgTGgAci4r/MsM0IMALQ29u7ZuPGjbnHWQWTk5OsXLmy7DAqy/WTzXXUXp3qZ+3atdsiYiBrvUIS\nAYCkvwCeiYgrWsoGgQ9HxFntth0YGIiJiYmcI6yGRqPhi13acP1kcx21V6f6kTSnRJDnqKGj0zMB\nJB0MvBH4gaRVaZmAs4G784rBzMyy5TlqaBWwQdJykoRzXUTcLGmzpKMBAXcA780xBjMzy5DnqKG7\ngJNnKB/Ka59mZjZ/vrLYzKzmnAjMzGrOicDMrOacCMzMas6JwMys5pwIzMxqzonAzKzmnAjMzGrO\nicDMrOacCMzMas6JwMys5pwIzMxqzonAzKzmnAjMzGrOicDMrOacCMzMas6JwMys5pwIzMxqzonA\nzKzmcksEkg6S9D1Jd0q6R9JlaflLJX1X0g5JX5Z0YF4xmJlZtjzPCJ4FhiLit4CTgDdJeh3wceDK\niHg58BPgghxjMDOzDLklgkhMposHpI8AhoCvpOUbgLPzisHMzLLl2kcgabmkO4DHgE3AA8BTEfFc\nusrDwDF5xmBmZu2tyPPNI2IPcJKkw4CvAifMtNpM20oaAUYAent7aTQaeYVZKZOTk7X5rAvh+snm\nOmrP9bO/XBNBU0Q8JakBvA44TNKK9KzgWODHs2wzBowBDAwMxODgYBGhlq7RaFCXz7oQrp9srqP2\nXD/7y3PU0NHpmQCSDgbeCGwHtgDnpqutA27MKwYzM8uW5xnBKmCDpOUkCee6iLhZ0r3ARkn/Fbgd\nuDrHGMzMLENuiSAi7gJOnqH8QeC1ee3XzMzmp23TkKTXSHpxy/L5km6U9N8kHZF/eGZmlresPoL/\nAfwCQNIpwOXAF4GnSTtyzcyss2U1DS2PiCfTn98GjEXE9cD16fUBZmbW4bLOCJZLaiaL04DNLa8V\nMvTUzKprfBz6+2HZsuR5fLzsiGwhsr7MrwW2SnoC+DnwLQBJx5M0D5lZTY2Pw8gITE0lyzt3JssA\nw8PlxWXz1/aMICJGgQ8B1wD/PiKaVwEvAy7ONzQzq7L16/cmgaapqaTcOkvbMwJJPcC2iPhluvwK\n4ExgZ0TcUEB8ZlZRu3bNr9yqK6uP4OtAP/yqOeg7wHHAn0r663xDM7MqW716fuVWXVmJ4PCI2JH+\nvA64NiIuBs4Azso1MjOrtNFR6OnZt6ynJym3zpKVCFpnBh0imUqaiPgF8HxeQZlZ9Q0Pw9gY9PWB\nlDyPjbmjuBNljRq6S9IVwCPA8cA3AJqTyZlZvQ0P+4u/G2SdEVwIPEHST/A7EdEcI/CbwBU5xmVm\nZgVpe0YQET8HLpd0EHC8pFcBD0TEPwP/XESAZmaWr6xJ51ZI+gTwI5L7C38J+JGkT0g6oIgAzcws\nX1lNQ38DHAEcFxFrIuJk4GXAYbhpyMysK2QlgrOACyPiZ82CiPgp8CckF5aZmVmHyxw+2jKtRGvh\nHma56byZmXWWrERwr6TzpxdKeifwg3xCMjOzImVdR/CnwA2S3gNsIzkLeA1wMPCWnGMzM7MCZA0f\nfQT4bUlDwKsAAV+LiFuLCM7MzPI3p5vLRMRm9r0pDZJ2RcSs00tJegnJbS1fTDIdxVhE/K2kS0ku\nVHs8XfVjEXHLAmI3M7MlsJi7jCnj9eeAD0XEbZIOBbZJ2pS+dmVEePipmVkFLCYRtB01FBG7gd3p\nzz+TtB04ZhH7MzOzHGiG0aF7X5Q+ONtLwPqIOGJOO5H6gW8CJwIfBN4F/BSYIDlr+MkM24wAIwC9\nvb1rNm7cOJdddbzJyUlWrlxZdhiV5frJ5jpqr071s3bt2m0RMZC1XlYi+It2G0fEZZk7kFYCW4HR\niLhBUi/JRHYB/CWwKiLe0+49BgYGYmJiImtXXaHRaDA4OFh2GJXl+snmOmqvTvUjaU6JIKtp6OqI\neHiWHfzeHII4ALgeGG/e2jIiHm15/XPAzVnvY2Zm+cm6oOzWtFlnH5LeDXy63YaSBFwNbI+IT7WU\nr2pZ7Rzg7rkGa2ZmSy/rjOADwCZJZzZvWSnpEuAdwKkZ274B+EPg+5LuSMs+Bpwn6SSSpqGHgD9e\nYOxmZrYEsi4ou0XSs8DXJJ0N/BHJlcWnzNTBO23bbzPzEFNfM2BmViFZTUOkVxG/C2gAxwGnZSUB\nMzPrHG3PCCT9jKQJR8ALgNOAx9L2/4iIF+YfopmZ5SmraejQogIxM7NyZN2q8jWSzpih/Pckrckv\nLDMzK8pcblW5fYby7elrZmbW4bISwZER8dD0woi4Hzgyl4jMzKxQWYng4DavHbKUgZiZWTmyEsE/\nSRpNRwn9iqTLmHZ/AjMz60xZVxZ/CPg8cH/L1cG/RTJr6IV5BmZmZsXIGj76DMmUEMeR3KoS4J6I\neDD3yMzMrBBzvVXlg8CDAJJeJuk/AW+PiBPzDM7MzPKXOcUEJDOGSnq/pO8B9wDLgfNyjczMzAqR\ndUHZhZI2k9xY5iiSSed2R8RlEfH9IgI0M7N8ZTUNfQb4DvCOiJgAkNT2XsVmZtZZshLBrwNvBT6V\n3mLyOuCA3KMyM7PCtG0aiognIuKqiDiFZObRp0lmH90u6a8KidDMzHKV1UfwuubPEfFwRFwREWuA\ns4Fn8w7OzMzylzVq6LMzFUbE/42Iy3KIx8zMCjan4aNmZta9shLBcZJumu3RbkNJL5G0Je1PuEfS\n+9LyIyRtkrQjfT58yT6NmVmXGB+H/n5Ytix5Hh/Pb19Zo4YeBz65wPd+DvhQRNwm6VBgm6RNJPc/\nvjUiLpf0UeCjwEcWuA8zs64zPg4jIzA1lSzv3JksAwwPL/3+shLBZERsXcgbR8RuYHf6888kbQeO\nAd4MDKarbQAaOBGYmf3K+vV7k0DT1FRSnkciUMTs14elVxW/IyL+NV0+H/gDYCdwaUQ8OaedSP3A\nN4ETgV0RcVjLaz+JiP2ahySNACMAvb29azZu3DjHj9TZJicnWblyZdlhVJbrJ5vrqL1OqJ+hoVOJ\n0H7lUrB589z/N1+7du22iBjIWi8rEdwGvDEinpR0CrARuBg4CTghIs7N3IG0kmSKitGIuEHSU3NJ\nBK0GBgZiYmIia1ddodFoMDg4WHYYleX6yeY6aq8T6qe/P2kOmq6vDx56aO7vI2lOiSCrs3hZy3/9\nbwPGIuL6iPjPwPFzCOIA4HpgPCJuSIsflbQqfX0V8FjW+5iZ1cnoKPT07FvW05OU5yErEayQ1OxH\nOI1970rWtn8hvavZ1cD2iPhUy0s3AevSn9cBN849XDOz7jc8DOvWwfLlyfLy5clyHv0DkJ0IrgW2\nSroR+DnwLQBJx5NMN9HOG4A/BIYk3ZE+zgQuB06XtAM4PV02M7PU+Dhs2AB79iTLe/Yky3kNIc26\nQ9mopFuBVcA3Ym+HwjKSvoJ2234b2L+3I3HafAM1M6uLokcNZd6hLCL+ZYay+5Y+FDMzA9i1a37l\ni+UpJszMKmb16vmVL5YTgZlZxVRt1JCZ2awuughWrAApeb7oorIj6g7DwzA2llw3ICXPY2P5jRrK\n7CMwM5vJRRfBVVftXd6zZ+/yZ2ecwN7mY3g4vy/+6XxGYGYLMjY2v3KrLicCM1uQ5hj3uZZbdTkR\nmNmCNK96nWu5VZcTgZktSHN+/LmWW3W5s9jMFqTZITw2ljQHLV+eJAF3FHceJwIzW7DPftZf/N3A\nTUNmZjXnRGBmVnNOBGZmNedEYGZWc04EZmY150RgZlZzTgRmZjXnRGBmVnO5JQJJX5D0mKS7W8ou\nlfTItJvZm5lZifI8I7gGeNMM5VdGxEnp45Yc92+2KOPj0N8Py5Ylz+PjZUdUPXWpo27/nLlNMRER\n35TUn9f7m+VpfDyZN2dqKlneuXPvZGpF3Syk6upSR3X4nGX0EfyZpLvSpqPDS9i/Wab16/f+4TdN\nTSXllqhLHdXhcyoi8nvz5Izg5og4MV3uBZ4AAvhLYFVEvGeWbUeAEYDe3t41GzduzC3OKpmcnGTl\nypVlh1FZRdXP0NCpRGi/cinYvHlr7vtfDNdRe/Otn079nABr167dFhEDmStGRG4PoB+4e76vTX+s\nWbMm6mLLli1lh1BpRdVPX18E7P/o6ytk94viOmpvvvXTqZ8zIgKYiDl8xxbaNCRpVcviOcDds61r\n1qrZWTc0dGohnXWjo9DTs29ZT09Sbomy6qjojttaHAtzyRYLeQDXAruBXwIPAxcAfw98H7gLuImk\nachnBC18RrC/L30poqdn3//GenqS8rz329cXISXPee9vqRR5DBVdR0txLCykfjr1WGCOZwS59hEs\nlYGBgZiYmCg7jEI0Gg0GBwfLDqNS+vuTkRrT9fXBQw8VHU31dfMxtBTHQjfXz3SS5tRH4CuLrfJ2\n7ZpfuXUvHwv5cCKwylu9en7l1r18LOTDicAqrxaddTYnPhby4URglTc8DGNjSTuwFPT1JcvdclVn\nHS105M++xwI+FpZIblNMmC2l4eHk0WhsrU1HX7da7JQNzWPBlo7PCMysUHWYsqHTOBGYWaE88qd6\nnAjMrFAe+VM9TgRmViiP/KkeJwIzK5RH/sxNkXMqedSQmRXOI3/aK/pmOD4jMDOrmKJHVjkRmJlV\nTNEjq5wIzMwqpuiRVU4EZrMo+gYoZk1Fj6xyIjCbQbOzbufO5PYnzc46JwMrQtEjq5wIzGbgaRCs\nbMPDyc12nn8+ec5zlJUTgdkMPA2C1YkTgdkMPA2C1UluiUDSFyQ9JunulrIjJG2StCN9Pjyv/Zst\nhqdBsDrJ84zgGuBN08o+CtwaES8Hbk2XzSrH0yBYneQ2xUREfFNS/7TiNwOD6c8bgAbwkbxiMFsM\nT4NgdVF0H0FvROwGSJ9fVPD+zcxsmspOOidpBBgB6O3tpdFolBtQQSYnJ2vzWRfC9ZPNddSe62d/\nRSeCRyWtiojdklYBj822YkSMAWMAAwMDUZf71DYaDd+Ttw3XTzbXUXuun/0V3TR0E7Au/XkdcGPB\n+zczs2nyHD56LfAd4BWSHpZ0AXA5cLqkHcDp6bJZpua8P0NDp3reH7MllueoofNmeem0vPZp3Wnf\nm3Qo95t0mNWNryy2yvO8P2b5ciKwyvO8P2b5ciKwyvO8P2b56tpE4JuKVNNCfi+e98da+W976VX2\ngrLF2LdzEXcuVsRCfy/N19avh127gtWrxeiof5d15L/tfHTlGYE7F6tpMb+X5k06Nm/emvtNOqy6\n/Ledj65MBO5crCb/XmyxfAzloysTgTsXq8m/F1ssH0P56MpE4M7FavLvxRbLx1A+ujIR+KYi1eTf\nSzV10igcH0P56MpRQ+CbilSVfy/V0omjcHwMLb2uPCMws7nxKBwDJwKzWvMoHAMnArNa8ygcAycC\ns1rzKBwDJwKzWvMoHIMuHjVkZnPjUTjmMwIzs5pzIjAzqzknAjOzmnMiMDOrOScCM7OaU0SUHUMm\nSY8DO8uOoyBHAU+UHUSFuX6yuY7aq1P99EXE0VkrdUQiqBNJExExUHYcVeX6yeY6as/1sz83DZmZ\n1ZwTgZlZzTkRVM9Y2QFUnOsnm+uoPdfPNO4jMDOrOZ8RmJnVnBNBSSS9QtIdLY+fSnq/pEslPdJS\nfmbZsZZF0gck3SPpbknXSjpI0kslfVfSDklflnRg2XGWaZY6ukbSD1uOoZPKjrMskt6X1s09kt6f\nlh0haVN6DG2SdHjZcZbNTUMVIGk58Ajw28C7gcmIuKLcqMol6Rjg28BvRsTPJV0H3AKcCdwQERsl\n/Xfgzoi4qsxYy9KmjgaBmyPiK2XGVzZJJwIbgdcCvwC+DvwJcCHwZERcLumjwOER8ZHyIi2fzwiq\n4TTggYioy0Vzc7UCOFjSCqAH2A0MAc0vuA3A2SXFVhXT6+jHJcdTJScA/xIRUxHxHLAVOAd4M8mx\nAz6GACeCqng7cG3L8p9JukvSF+p62hoRjwBXALtIEsDTwDbgqfSPGuBh4JhyIizfTHUUEd9IXx5N\nj6ErJb2gtCDLdTdwiqQjJfWQnE2+BOiNiN0A6fOLSoyxEpwISpa2cf8+8D/ToquAlwEnkfxxf7Kk\n0EqVJsA3Ay8Ffh04BDhjhlVr27Y5Ux1JeidwCfBK4DXAEUAtmz0iYjvwcWATSbPQncBzbTeqKSeC\n8p0B3BYRjwJExKMRsScingc+R9K+WUdvBH4YEY9HxC+BG4DXA4elzSAAx1LvppAZ6ygidkfiWeDv\nqO8xRERcHRGvjohTgCeBHcCjklYBpM+PlRljFTgRlO88WpqFmgdo6hyS09s62gW8TlKPJJH0o9wL\nbAHOTddZB9xYUnxVMFMdbW/5khNJ+3ddjyEkvSh9Xg28heRv7SaSYwd8DAEeNVSqtN3yR8BxEfF0\nWvb3JM1CATwE/HGzPbNuJF0GvI3kdP524I9I+gQ2kjR53A68M/3Pt5ZmqaOvAUcDAu4A3hsRk6UF\nWSJJ3wKOBH4JfDAibpV0JHAdsJokmb41Ip4sMczSORGYmdWcm4bMzGrOicDMrOacCMzMas6JwMys\n5pwIzMxqzonAupqkcySFpFe2lP2GpFsk3S9pu6TrJPVKGpR0c7rOuyQ9Lun2dJbKf5T0+mnv/WFJ\nP0hnt7xT0vlp+YGSPi3pgXTbGyUd27JdSPrktPe5tGX5/JYZM++V9OEcq8jMicC63nkkM3S+HUDS\nQcA/AFdFxPERcQLJtB5Hz7DtlyPi5Ih4OXA5cIOkE9L3eS9wOvDaiDgROIVk3D7AXwGHAr+Rbvu/\n0m2brz8LvEXSUdN3KOkM4P3A70TEq4BXk8yzZJYbJwLrWpJWAm8ALiBNBMA7gO9ExP9urhcRWyKi\n7dW3EbGF5BaHI2nRx4CLIuKn6etPR8SG9CLBdwMfiIg96Wt/R/LlP5Ru+1z6Xh+YYVeXAB+OiB+n\n2/5bRHxufp/cbH6cCKybnQ18PSLuA56U9GrgRJJZTBfiNuCVkg4FDo2IB2ZY53hgVzNBtJgAXtWy\n/BlgWNKvTVtvMfGZLYgTgXWz80imoyB9Pm+R76eW59kuyZ/ttX3K00TxReA/LDIms0VzIrCulM4n\nMwR8XtJDwJ+TzMlzD7BmgW97MrA9/RJ/RtJxM6xzP9CXnjW0ejXJpHmtPk3SbHVIS9li4jNbECcC\n61bnAl+MiL6I6I+IlwA/BO4DXi/pd5srSnqTpH/X7s0knUrSP9Bsr/9r4DOSXpi+/kJJIxHxDMld\nrz6V3oKUdDRRD7C59T3Tic6uI0kGtLzvJyS9ON32BZJ81mC5ciKwbnUe8NVpZdeTdBafBVycDu28\nF3gXM89J/7b05u/3kXQO/0F6sxNIRhptAf6PpLtJboM4lb52CfBvwH2SdgBvBc6JmWd4/CTwq9FD\nEXELSf/BP0m6h6S/YMUM25ktGc8+amZWcz4jMDOrOScCM7OacyIwM6s5JwIzs5pzIjAzqzknAjOz\nmnMiMDOrOScCM7Oa+/+acMNre0bv8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x126a79c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# distribution of STACKLOSS vs other variables\n",
    "\n",
    "plt.figure(1)\n",
    "plt.xlabel('AIRFLOW')\n",
    "plt.ylabel('STACKLOSS')\n",
    "plt.grid(True)\n",
    "plt.plot(X['AIRFLOW'], Y, 'o', color='blue')\n",
    "\n",
    "plt.figure(2)\n",
    "plt.xlabel('WATERTEMP')\n",
    "plt.ylabel('STACKLOSS')\n",
    "plt.grid(True)\n",
    "plt.plot(X['WATERTEMP'], Y, 'o', color='blue')\n",
    "\n",
    "plt.figure(3)\n",
    "plt.xlabel('ACIDCONC')\n",
    "plt.ylabel('STACKLOSS')\n",
    "plt.grid(True)\n",
    "plt.plot(X['ACIDCONC'], Y, 'o', color='blue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>STACKLOSS</th>\n",
       "      <th>AIRFLOW</th>\n",
       "      <th>WATERTEMP</th>\n",
       "      <th>ACIDCONC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>STACKLOSS</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.919663</td>\n",
       "      <td>0.875504</td>\n",
       "      <td>0.399830</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AIRFLOW</th>\n",
       "      <td>0.919663</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.781852</td>\n",
       "      <td>0.500143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>WATERTEMP</th>\n",
       "      <td>0.875504</td>\n",
       "      <td>0.781852</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.390940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACIDCONC</th>\n",
       "      <td>0.399830</td>\n",
       "      <td>0.500143</td>\n",
       "      <td>0.390940</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           STACKLOSS   AIRFLOW  WATERTEMP  ACIDCONC\n",
       "STACKLOSS   1.000000  0.919663   0.875504  0.399830\n",
       "AIRFLOW     0.919663  1.000000   0.781852  0.500143\n",
       "WATERTEMP   0.875504  0.781852   1.000000  0.390940\n",
       "ACIDCONC    0.399830  0.500143   0.390940  1.000000"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Correlation coefficient\n",
    "Z = pd.concat([Y, X], axis=1)\n",
    "Z.corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Estimated Parameters ---\n",
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "--- Standard Error of Estimated Parameters ---\n",
      "[9.79189854 0.11100521 0.30293016 0.12864961]\n",
      "--- Result ---\n",
      "                    Robust linear Model Regression Results                    \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                   21\n",
      "Model:                            RLM   Df Residuals:                       17\n",
      "Method:                          IRLS   Df Model:                            3\n",
      "Norm:                          HuberT                                         \n",
      "Scale Est.:                       mad                                         \n",
      "Cov Type:                          H1                                         \n",
      "Date:                Thu, 02 May 2019                                         \n",
      "Time:                        01:46:19                                         \n",
      "No. Iterations:                    19                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "var_0        -41.0265      9.792     -4.190      0.000     -60.218     -21.835\n",
      "var_1          0.8294      0.111      7.472      0.000       0.612       1.047\n",
      "var_2          0.9261      0.303      3.057      0.002       0.332       1.520\n",
      "var_3         -0.1278      0.129     -0.994      0.320      -0.380       0.124\n",
      "==============================================================================\n",
      "\n",
      "If the model instance has been used for another fit with different fit\n",
      "parameters, then the fit options might not be the correct ones anymore .\n"
     ]
    }
   ],
   "source": [
    "# add constant term to exogenous variable as well as the case of OLS \n",
    "data.exog = sm.add_constant(data.exog)\n",
    "\n",
    "# M : robust criterion function for downweighting outliers. \n",
    "# There are other options : LeastSquares, HuberT, RamsayE, AndrewWave, TrimmedMean, Hampel, and TukeyBiweight \n",
    "# See the shapes of these weight fucntions : https://www.statsmodels.org/stable/rlm_techn1.html\n",
    "\n",
    "huber_t = sm.RLM(data.endog, data.exog, M=sm.robust.norms.HuberT()) # use Huber weight function\n",
    "hub_results = huber_t.fit()\n",
    "\n",
    "print('--- Estimated Parameters ---')\n",
    "print(hub_results.params)\n",
    "\n",
    "print('--- Standard Error of Estimated Parameters ---')\n",
    "print(hub_results.bse)\n",
    "\n",
    "print('--- Result ---')\n",
    "print(hub_results.summary(yname='y', xname=['var_%d' % i for i in range(len(hub_results.params))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# prediction of training data(not test data)\n",
    "huber_predicted = hub_results.predict(data.exog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x126cc5128>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VOW9x/HPD4QqglaK4oKArd5K\niLKKiKBAFahLtVZbvHBFDQRRFC/uIq0bynJVBESJhIISQWQpyK5AkH0HA2HRIqAipQqKAdmS3/1j\nBsuSQEgyc5KZ7/v1mteZ88xZfs8rw/x4nnPO85i7IyIi8atU0AGIiEiwlAhEROKcEoGISJxTIhAR\niXNKBCIicU6JQEQkzikRiIjEOSUCEZE4p0QgIhLnTgk6gPyoVKmSV69evUD77t69m9NPP71oAyrm\nVOf4oDrHvsLWd9myZd+6+9kn2q5EJILq1auzdOnSAu2bnp5O06ZNizagYk51jg+qc+wrbH3NbHN+\ntlPXkIhInFMiEBGJc0oEIiJxTolARCTOKRGIiMS5iN41ZGabgB+BbOCgu9c3s4rA+0B1YBPwZ3ff\nGck4REQkb9FoETRz99ruXj+8/iQww90vAWaE10VEJCBBdA3dAgwLvx8G3BpADCIixdqWLVsYMGAA\nBw8ejPi5LJJzFpvZF8BOwIFB7p5iZt+7+y8P22anu5+Vy77JQDJA5cqV640cObJAMWRlZVG+fPkC\n7VtSqc7xQXWOTTk5OUyYMIGUlBRycnJ4/fXX+e1vf1ugYzVr1mzZYb0xeXP3iL2A88PLc4BVwDXA\n90dts/NEx6lXr54X1KxZswq8b0mlOscH1Tn2rF+/3ps0aeKAX3/99T5ixIhCHQ9Y6vn4rY5o15C7\nbw0vtwPjgAbAv8zsPIDwcnskYxARKe4OHjxIr169uPzyy8nIyODvf/8706ZN49xzz43K+SOWCMzs\ndDOrcOg90AJYDUwA2oU3aweMj1QMIiLF3cqVK7nyyit58sknufHGG1m7di133303Zha1GCJ5+2hl\nYFy4MqcA77n7VDNbAowysyRgC3BHBGMQESmW9u7dywsvvECvXr2oVKkSo0eP5k9/+lMgsUQsEbj7\nRqBWLuXfAb+L1HlFRIq7+fPnk5SUxLp162jXrh2vvvoqFStWDCwePVksIhIlWVlZPPTQQzRu3Jg9\ne/YwdepUhg4dGmgSACUCEZGomD59OomJiQwYMIDOnTuzevVqWrZsGXRYgBKBiEhE7dixg3vuuYeW\nLVty6qmnMmfOHPr160eFChWCDu1nSgQiIhEyZswYEhISePfdd3n66adZuXIlV199ddBhHaNETFUp\nIlKSbNu2jc6dOzNmzBjq1KnD1KlTqV27dtBh5UktAhGRIuLuDB06lISEBCZOnMjLL7/MokWLinUS\nALUIRESKxKZNm+jYsSPTp0+ncePGDB48uMBjBEWbWgQiIoWQk5ND//79SUxMZP78+bzxxhvMnj27\nxCQBUItARKTA1q5dS/v27Zk/fz6tWrXirbfeolq1akGHddLUIhAROUkHDhzgpZdeonbt2qxbt453\n3nmHyZMnl8gkAGoRiIiclOXLl5OUlMTKlSu544476N+/P5UrVw46rEJRi0BEJB9++uknnnrqKRo0\naMC2bdsYO3Yso0aNKvFJANQiEBE5oblz55KUlMSGDRtISkqiT58+nHXWMRMrllhqEYiI5OHHH3+k\nc+fONGnShP379/PRRx8xePDgmEoCoEQgIpKrKVOmULNmTQYOHMjDDz/M6tWrue6664IOKyKUCERE\nDvPdd99x1113ccMNN1C+fHnmzZvHa6+9xumnnx50aBGjRCAiQmh4iA8++ICEhARGjBhB9+7dWbFi\nBVdddVXQoUWcLhaLSNz75ptvuP/++/nHP/5BvXr1mD59OrVqHTPBYsxSi0BE4pa7M2TIEGrUqMHU\nqVPp3bs3CxcujKskAGoRiEic2rhxI8nJycyYMYNrrrmGwYMHc8kllwQdViDUIhCRuJKdnU3fvn25\n7LLLWLx4MW+++SazZs2K2yQAahGISBzJzMwkKSmJhQsXcsMNN/DWW29x4YUXBh1W4NQiEJGYt3//\nfl544QXq1KnDZ599xvDhw5k4caKSQJgSgYgUC2lpUL06lCoVWqalFc1xly5dyhVXXMFf//pXbrvt\nNjIzM2nTpg1mVjQniAFKBCISuLQ0SE6GzZvBPbRMTi5cMtizZw+PP/44V155Jd9++y3jx49nxIgR\nnHPOOUUXeIxQIhCRwHXrBnv2HFm2Z0+ovCBmz55NrVq16NOnD0lJSWRmZvKHP/yh8IHGKCUCEQnc\nli0nV56XXbt20alTJ5o2bUpOTg4zZswgJSWFM888s/BBxjAlAhEJXNWqJ1eem0mTJlGzZk1SUlLo\n2rUrGRkZNG/evGgCjHFKBCISuB49oFy5I8vKlQuVn8i3335L27ZtuemmmzjzzDOZP38+r7zyCuWO\nPqDkSYlARALXpg2kpEC1amAWWqakhMrz4u6MHDmSGjVqMGrUKJ599lmWL1/OlVdeGb3AY0TEHygz\ns9LAUuBrd7/JzC4CRgIVgeXA/7j7/kjHISLFW5s2x//hP9zXX39Np06d+PDDD2nQoAGpqakkJiZG\nNsAYFo0WQRdg7WHrvYDX3P0SYCeQFIUYRCQGuDtvv/02CQkJfPzxx7zyyivMnz9fSaCQIpoIzKwK\ncCMwOLxuQHNgdHiTYcCtkYxBRGLDP//5Tx555BGSk5OpV68eGRkZdO3aldKlSwcdWokX6RZBX+Bx\nICe8/ivge3c/GF7/CrggwjGISAmWnZ3Nq6++ymWXXcaGDRtISUlhxowZ/OY3vwk6tJgRsWsEZnYT\nsN3dl5lZ00PFuWzqeeyfDCQDVK5cmfT09ALFkZWVVeB9SyrVOT7EQ52/+OILevfuzbp162jUqBHJ\nyclUq1aN2bNnBx1aVETtb+zuEXkBLxP6H/8mYBuwB0gDvgVOCW9zFTDtRMeqV6+eF9SsWbMKvG9J\npTrHh1iu8759+/xvf/ublylTxitVquQjRozwnJycmK5zbgpbX2Cp5+P3OmJdQ+7+lLtXcffqQGtg\npru3AWYBt4c3aweMj1QMIlLyLF68mLp16/Lcc8/x5z//mbVr19K6dWsNEhdBQTxH8ATQ1cw+J3TN\nIDWAGESkmNmzZw+PPPIIV111FT/88AMTJ05k+PDhVKpUKejQYl5UJqZx93QgPfx+I9AgGucVkZJh\n1qxZtG/fno0bN3LffffRq1cvzjjjjKDDiht6slhEAvPDDz+QnJxM8+bNKVWqFOnp6bz55ptKAlGm\nRCAigZgwYQIJCQmkpqby2GOPsWrVKq699tqgw4pLSgQiElXbt2+ndevW3HLLLfzqV79i0aJF9O7d\nW4PEBUiJQESiwt1JS0sjISGBcePG8cILL7B06VLq168fdGhxT4lAJA5Fan7gvHz55ZfcfPPNtG3b\nlksuuYQVK1bwzDPPULZs2UDjkpCo3DUkIsXHofmBD00NeWh+YMj/6J/5lZOTQ0pKCo8//jjZ2dn0\n7duXzp075zo+UDTjkiOpRSASZ4p6fuC8fPbZZzRv3pxOnTrRoEEDMjIy6NKlS56DxEUrLjmWEoFI\nnCmq+YHzcvDgQfr06cPll1/OypUrSU1N5aOPPuLXv/51oHFJ3pQIROJMUcwPnJdVq1bRsGFDHn/8\ncVq1akVmZib33ntvvoaHiGRccnxKBCJxpjDzA+dl3759dO/enfr16/Pll18yatQoxo4dy/nnnx9o\nXJI/SgQicaYg8wMfz4IFC6hTpw4vvvgi//3f/01mZiZ33HHHSQ8SV9RxSf7priGROHQy8wPnZffu\n3XTr1o1+/fpRpUoVJk+ezO9///vA45KTp0QgIift448/pkOHDmzatIkHHniAl19+mQoVKgQdlhSQ\nuoZEJN927txJUlIS119/PWXLluWTTz5hwIABSgIlnBKBiOTLuHHjSEhIYNiwYTz55JOsWrWKJk2a\nBB2WFAF1DYnIcf3rX//iwQcf5IMPPqB27dpMmjSJunXrBh2WFCG1CEQkV+7OO++8Q40aNRg/fjw9\nevT4eRpJiS1qEYjIMbZs2ULHjh2ZOnUqjRo1IjU1lUsvvTTosCRC1CIQkZ/l5OTwxhtvULNmTebM\nmUO/fv2YM2eOkkCMU4tARABYv3497du3Z+7cubRo0YJBgwZRvXr1oMOSKFCLQCTOHThwgJ49e1Kr\nVi3WrFnD0KFDmTp1qpJAHFGLQCSOrVixgqSkJFasWMGf/vQnBgwYwLnnnht0WBJlahGIxKG9e/fS\nrVs3rrjiCrZu3cro0aMZPXq0kkCcUotAJM7MmzePpKQk1q9fz913380rr7xCxYoVgw5LAqQWgUic\nyMrK4qGHHqJJkybs3buXadOm8fe//11JQJQIROLBtGnTqFmzJgMGDODBBx9k9erVtGjRIuiwpJhQ\nIhCJYTt27ODuu++mVatWlCtXjjlz5vD6669Tvnz5oEOTYkSJQCRGjRkzhoSEBIYPH063bt1YsWIF\nV199ddBhSTGki8UiMeabb76hc+fOjB07ljp16jB16lRq164ddFhSjKlFIBIj3J2hQ4eSkJDApEmT\n6NmzJ4sXL1YSkBNSi0AkBmzatInk5GQ++ugjmjRpwttvv81vf/vboMOSEiJiLQIzO9XMFpvZKjNb\nY2bPhcsvMrNFZvaZmb1vZmUjFYNIrMvOzqZfv34kJiayYMEC3njjDdLT05UE5KREsmtoH9Dc3WsB\ntYFWZtYQ6AW85u6XADuBpAjGIBKzNm/ezDXXXEOXLl1o0qQJa9as4f7776dUKfX4ysmJ2DfGQ7LC\nq2XCLweaA6PD5cOAWyMVg0gsOnDgAD169KBDhw6sW7eOd955h8mTJ1O1atWgQ5MSytw9cgc3Kw0s\nAy4G3gD6AAvd/eLw5xcCU9w9MZd9k4FkgMqVK9cbOXJkgWLIysqKu3umVefYtWHDBnr37s0///lP\nGjduzP/+7//G1ZPB8fJ3PqSw9W3WrNkyd69/wg3dPeIv4JfALKAJ8Plh5RcCGSfav169el5Qs2bN\nKvC+JZXqHHv27NnjTzzxhJcuXdrPPfdcHzduXMzXOTfxVufC1hdY6vn4jY5KZ6K7fw+kAw2BX5rZ\nobuVqgBboxGDSEn1ySefUKtWLXr16sXdd99NZmYmt96qHlUpOpG8a+hsM/tl+P1pwHXAWkItg9vD\nm7UDxkcqBpGSbNeuXTzwwANce+21HDx4kI8//pjBgwdz1llnBR2axJhIPkdwHjAsfJ2gFDDK3Sea\nWSYw0sxeBFYAqRGMQaREmjJlCh07duSrr77i4Ycf5sUXX+T0008POiyJUZG8a+hTd6/j7pe7e6K7\nPx8u3+juDdz9Yne/w933RSoGEYC0NKheHUqVCi3T0oKOKG/fffcdd911FzfccAMVKlRg/vz5vPba\na0oCElG64VhiWloaJCfD5s3gHlomJxe/ZODujBo1iho1ajBixAi6d+/O8uXLadiwYdChSRw4YSIw\ns2OGK8ytTKQ46tYN9uw5smzPnlB5cbF161b++Mc/8pe//IWqVauybNkynn/+eX7xi18EHZrEify0\nCPrns0yk2Nmy5eTKo8ndSU1NJSEhgWnTptGnTx8WLlzI5ZdfHnRoEmfyvFhsZlcBjYCzzazrYR+d\nAZSOdGAiRaFq1VB3UG7lQdq4cSMdOnRg5syZXHvttQwePJiLL7442KAkbh2vRVAWKE8oWVQ47LWL\n/9z+KVKs9egB5codWVauXKg8CNnZ2fTt25fLLruMJUuW8NZbbzFz5kwlAQlUni0Cd58NzDazoe6+\n2cxOd/fdUYxNpNDatAktu3ULdQdVrRpKAofKo2nNmjUkJSWxaNEibrzxRt566y2qVKkS/UBEjpKf\nawTnh+/9XwtgZrXMbGBkwxIpOm3awKZNkJMTWkY7Cezfv5/nn3+eOnXq8Pnnn5OWlsaHH36oJCDF\nRn4eKOsLtAQmALj7KjO7JqJRicSIJUuWkJSUREZGBq1bt6Zfv36cffbZQYclcoR8PUfg7l8eVZQd\ngVhEYsaePXt47LHHaNiwId999x3jx49nxIgRSgJSLOWnRfClmTUCPDyb2EOEu4lE5Fjp6el06NCB\nzz//nOTkZHr37s2ZZ54ZdFgiecpPi+A+4AHgAuArQrONPRDJoERKoh9++IH77ruPZs2a4e7MnDmT\nQYMGKQlIsXfCFoG7fwsEcI+FSMkxadIkOnbsyDfffMMjjzzC888/T7mj71sVKaZOmAjMrF8uxT8Q\nmvBAQ0hLXPv3v//Nww8/zHvvvUdiYiJjx46lQYMGQYclclLy0zV0KqHuoM/Cr8uBikCSmfWNYGwi\nxZa7M2LECBISEvjggw949tlnWbZsmZKAlEj5uVh8MdDc3Q8CmNmbwHTgeiAjgrGJFEtfffUVnTp1\nYuLEiTRo0IDU1FQSE4+ZdlukxMhPi+AC4PDB0E8Hznf3bEBzCUjcyMnJISUlhZo1azJjxgxeffVV\n5s+fryQgJV5+WgS9gZVmlg4YcA3wkpmdDnwcwdhEio3PP/+cDh06kJ6eTrNmzXj77bf5zW9+E3RY\nIkXiuInAzIxQN9BkoAGhRPC0ux+acP6xyIYnEqyDBw/y+uuv0717d8qUKcPbb79NUlISoX8aIrHh\nuInA3d3M/uHu9dAk8xJnMjIySEpKYsmSJdx88828+eabXHDBBUGHJVLk8nONYKGZXRHxSESKiX37\n9vG3v/2NunXrsmnTJkaOHMn48eOVBCRm5ecaQTOgo5ltBnYT6h5yd9c0ShJzFi1aRFJSEmvWrKFt\n27a89tprVKpUKeiwRCIqP4ng9xGPQiRgu3fvpnv37vTt25cLLriAiRMncuONNwYdlkhU5GeIic0A\nZnYOoYfLRGLKzJkz6dChAxs3bqRTp0707NmTM844I+iwRKLmhNcIzOwPZvYZ8AUwG9gETIlwXCIR\n9/3339OhQwd+97vfUapUKdLT0xk4cKCSgMSd/FwsfgFoCGxw94uA3wHzIhqVSISNHz+ehIQEhgwZ\nwuOPP86nn37KtddeG3RYIoHITyI44O7fAaXMrJS7zyI09pBIibN9+3Zat27Nrbfeytlnn82iRYvo\n1asXp512WtChiQQmPxeLvzez8sAnQJqZbQcORDYskaLl7qSlpdGlSxeysrJ44YUXeOKJJyhTpkzQ\noYkELj+JYBWwB/hfQvMSnAmUj2RQIkXpyy+/5L777mPy5Mk0bNiQ1NRUEhISgg5LpNjI13ME7p4D\n5ADDAMzs04hGJVIEcnJyGDRoEE888QTZ2dn07duXzp07U7p06aBDEylW8kwEZtYJuB/4zVE//BXQ\nxWIp5jZs2ED79u2ZM2cO1113HSkpKVx00UVBhyVSLB3vYvF7wM2Exhi6+bBXPXdve6IDm9mFZjbL\nzNaa2Roz6xIur2hmH5nZZ+HlWUVQDxEgNEhc7969qVWrFhkZGQwZMoTp06crCYgcR54tAnf/gdCU\nlHcW8NgHgUfcfbmZVQCWmdlHwN3ADHfvaWZPAk8CTxTwHCI/W7VqFbfcMorNmzsAj/KrX+VQtuwp\naKBQkePLzzWCAnH3b4Bvwu9/NLO1hCa5uQVoGt5sGJCOEoEUwr59+xgyZAjDhzvug4DQpPFff12K\n5OTQNm3aBBefSHGXn+cICs3MqgN1gEVA5XCSOJQszolGDBKbFixYQJ06dXj33XcpV+41DiWBQ/bs\ngW7dgolNpKQwd4/sCULPIMwGerj7WDP73t1/edjnO939mOsEZpYMJANUrly53siRIwt0/qysLMqX\nj6+7XeOhzj/99BOpqamMHTuWc845h06dOvHcc3/F/dh+IDNn5szZAUQZWfHwdz5avNW5sPVt1qzZ\nMnevf8IN3T1iL6AMMA3oeljZeuC88PvzgPUnOk69evW8oGbNmlXgfUuqWK/z9OnTvXr16g74Aw88\n4Lt27fJZs2Z5tWrucOyrWrWgI46MWP875ybe6lzY+gJLPR+/1RHrGgpPc5kKrHX3Vw/7aALQLvy+\nHZr5TPJp586d3HvvvbRo0YKyZcvyySefMGDAACpUqADADTfkvl9e5SISEslrBFcD/wM0N7OV4dcN\nQE/g+vCIpteH10WOa9y4cSQkJPDOO+/w1FNPsWrVKpo0aXLENpMn575vXuUiEhLJu4bmEprNLDe/\ni9R5JbZs27aNBx98kNGjR1O7dm0mTZpE3bp1c912y5bcj5FXuYiEROWuIZGT5e4MGzaMhIQEPvzw\nQ1566SUWL178cxJIS4Pq1aFUqdDy44/PoWLF3I9VtWrUwhYpkSLWIhApqM2bN9OxY0emTZtGo0aN\nSE1N5dJLL/3587Q0SE4O3Roa2h569rw012OVLQs9ekQjapGSSy0CKTZycnIYMGAANWvWZO7cufTv\n3585c+YckQQg9FzAoSRwSHZ2KbKzjz1mhQp6mEzkRNQikGJh/fr1JCUlMW/ePFq2bMmgQYOoVq1a\nrtueTJ//jh1FFKBIDFOLQAJ14MABXn75ZWrVqkVmZiZDhw5lypQpeSYBOLk+f10fEDkxJQIJzIoV\nK7jyyit5+umnufnmm8nMzKRdu3bYCUaJ69EDyh05kgSlS+dQtuyRZeXK6fqASH4oEUjU7d27l6ef\nfporrriCrVu3MmbMGD744APOPffcfO3fpg2kpEC1amAWWj755DqGDDmyLCVF1wdE8kPXCCSq5s6d\nS1JSEhs2bOCee+7hlVde4ayzTn5KijZtjvyRT0/fTtOmCfrhFykAtQgkKn788Uc6d+5MkyZN2Ldv\nH9OmTWPIkCEFSgIiUrSUCCTipk2bRmJiIgMHDuShhx5i9erVtGjRIuiwRCRMiUAiZseOHbRr145W\nrVpRrlw55s6dy+uvvx5XwwiLlARKBBIRo0ePpkaNGrz33nt069aNFStW0KhRo6DDEpFcKBHIcR09\npk9a2vG3/+abb7jtttu44447qFKlCkuWLOHFF1/k1FNPLfSxg1ASYhQpLN01JHnKbUyfvOYAdneG\nDh1K165d+emnn+jZsyePPPIIp5yS+1fsZI4dlJIQo0hRUItA8pTbmD65zQH8xRdf0KJFC+69914u\nu+wyPv30U5544ok8k8DJHDtIJSFGkaKgRCB5OtH4/tnZ2fTr14/ExEQWLlzIwIEDSU9P57/+678K\nfezioCTEKFIUlAgkz37wvMbpqVoV1q5dS5MmTejSpQvXXnsta9asoVOnTpQqlb+v1PGOXVyUhBhF\nioISQZw71A++eXNoqvdD/eBpaXnN9etUrryE2rVrs379et59910mTZpE1ZP8dcxtvKDiNjZQSYhR\npCgoEcS54/WD5z7Xr7F48dnceuutrF27lrZt255wkLjc5DZeUHEbG6gkxChSFHTXUJwrSD+4WTXe\nf//9Qp/76PGCiqOSEKNIYalFEGeOvh5wvHl+8+4jP/kWgIgUX0oEcSS36wE//ghlyhy5Xbly8Mwz\ne7j44iHA7mM+Ux+5SGxRIogjuV0P2L8fzjjjyH7w++9fwfPPX8rMme1p1WocF16Yoz5ykRimawRx\nJK9+/x074Ntv4dtvv+XWW0fxf/93A7CJc889QNu2v9APv0iMU4sgjuTV53/hhc7777/Pr3/9DPPm\ntQOqA6XYtu0XP99KKiKxS4kgjuR2X/xpp+VQqdKrtG7dmn37/gacfsTnGlJBJPYpEcSRI++LdypW\n/BFIJjPzGfr06cOBA7nPGawhFURimxJBnGnTBmbO3EizZtexY8cZNGjwORkZGTz66KN53hbqriGY\nRWKZEkEcyc7O5rXXXiMxMZElS5YwaNAgZs6cycUXXwzk3nV0yOFDT4hIbFEiiBOrV6+mUaNGdO3a\nlebNm5OZmUlycvIRg8Qd3nWUG10vEIlNSgQxbv/+/Tz33HPUrVuXjRs38t577/Hhhx9SpUqVXLdv\n0wY2bQo9U5AbXS8QiT0RSwRmNsTMtpvZ6sPKKprZR2b2WXh5VqTOL7BkyRLq1avHs88+y+23305m\nZiZ33nlnvgaJ0xDMIvEjki2CoUCro8qeBGa4+yXAjPC6FLG9e/fy6KOP0rBhQ3bu3MmECRN47733\nOPvss/N9DA3BLBI/IvZksbt/YmbVjyq+BWgafj8MSAeeiFQM8Sg9PZ2kpCS2bt1Kx44d6dWrF2ee\neeZJH+fQ08TduoW6g6pWDSUBPWUsEnvM3SN38FAimOjuieH17939l4d9vtPdc+0eMrNkIBmgcuXK\n9UaOHFmgGLKysihfvnyB9i1JsrKyGDRoEBMnTuS8887jscceo06dOkGHFTXx8nc+nOoc+wpb32bN\nmi1z9/on3NDdI/YiNFbB6sPWvz/q8535OU69evW8oGbNmlXgfYub4cPdq1VzNwsthw8PlU+YMMHP\nP/98L1WqlD/66KM+ZcqUIMMMRCz9nfNLdY59ha0vsNTz8Rsb7buG/mVm5wGEl9ujfP4SK7chpDt0\ncBo1GsAf/vAHKlasyIIFC+jTpw+nnnpq0OGKSAkS7UQwAWgXft8OGB/l85dYuQ0h/dNPxoIFN/Pc\nc8+xbNkyGjRoEExwIlKiRfL20RHAAuC3ZvaVmSUBPYHrzewz4PrwuuRDXvfvm1Xlr3/9K2XLlo1u\nQCISMyKWCNz9Tnc/z93LuHsVd0919+/c/Xfufkl4uSNS5y9Ojp4eMj/DNBw7pWTuF/UrVtS0kSJS\nOHqyOMJy69s/0Zg9ue3z3XcHgAPHbLtrl8b/EZHCUSKIsNz69k80Zk9u+0BZoPQx2x44oPF/RKRw\nlAgiLK++/eON2ZP3Z7n/uTT+j4gUhhJBhOU1Nk+pUrl36ezbt48zzvg+131KH9sgOO45RETyQ4kg\nwvIa4z87+9hrBQsXLqRu3br88MP9lC6974jty5ULba/xf0SkqCkRRNihMf5z+9/8oWsFu3fvpmvX\nrjRq1Ihdu3YxaVJbhg37RXhKydD8ACkpMHDg4VNN/qdc4/+ISGFEbNA5+Y82beB//if3z7ZscS67\n7DK++OILOnXqRM+ePTnjjDN+3i+3Y+mHX0SKkhJBlFStGroN9GjumznllFOYPXs211xzTfQDE5G4\np66hKMn9WsFubrppAatWrVJu1/3UAAAJ9UlEQVQSEJHAKBFESZs20KfPD5Qr928ghzJltvL88//i\nww/v5LTTTgs6PBGJY0oEUeDuDB8+nO7df83Bg1V48cWX2b37bLp3/3XQoYmI6BpBpG3ZsoX77ruP\nKVOmcNVVV5GamkqNGjWCDktE5GdqEURITk4OAwcOpGbNmsyePZvXX3+dOXPmKAmISLGjFkEEbNiw\ngfbt2zNnzhyuu+46UlJSuOiii4IOS0QkV2oRFKGDBw/Sq1cvLr/8cjIyMhgyZAjTp09XEhCRYk0t\ngiKyatUq7r33XpYvX84f//hH3njjDc4777ygwxIROSG1CApp7969PPPMM9SvX5+vv/6a0aNHM3bs\nWCUBESkx1CIohPnz55OUlMS6deto164dr776KhUrVgw6LBGRk6IWQQFkZWXx0EMP0bhxY/bs2cPU\nqVMZOnSokoCIlEhKBCdp+vTpJCYm0r9/fx544AFWr15Ny5Ytgw5LRKTAlAjyaefOndxzzz20bNmS\nU089lTlz5tC/f38qVKgQdGgiIoWiRJAPY8eOJSEhgXfffZennnqKlStX0rhx46DDEhEpErpYfBzb\ntm2jc+fOjBkzhtq1azN58mTq1KkTdFgiIkVKLYJcuDtDhw4lISGBiRMn8tJLL7F48WIlARGJSWoR\nHGXTpk107NiR6dOnc/XVVzN48GAuvfTSoMMSEYkYtQjCcnJy6N+/P4mJicybN4/+/fvzySefKAmI\nSMxTiwBYt24d7du3Z968ebRs2ZJBgwZRrVq1oMMSEYmKuG4RHDhwgJdeeolatWqRmZnJsGHDmDJl\nipKAiMSVuG0RLF++nKSkJFauXMntt9/OgAEDqFy5ctBhiYhEXdy1CH766SeeeuopGjRowLZt2xgz\nZgwffPCBkoCIxK1AEoGZtTKz9Wb2uZk9Ga3zzp07l9q1a9OzZ0/uuusuMjMzue2226J1ehGRYinq\nicDMSgNvAL8HEoA7zSwhkuf88ccf6dy5M02aNGH//v1Mnz6dIUOGcNZZZ0XytCIiJUIQLYIGwOfu\nvtHd9wMjgVsidbLFixeTmJjIwIED6dKlCxkZGVx//fWROp2ISIkTxMXiC4AvD1v/CrgyEifq2LEj\nKSkp1KhRg3nz5nHVVVdF4jQiIiVaEInAcinzYzYySwaSASpXrkx6evrJn8iMv/zlL9x7773s27ev\nQMcoibKysuKmroeozvEh3uocrfoGkQi+Ai48bL0KsPXojdw9BUgBqF+/vjdt2vSkT9S0aVPS09Mp\nyL4lmeocH1Tn2Bet+gZxjWAJcImZXWRmZYHWwIQA4hAREQJoEbj7QTPrDEwDSgND3H1NtOMQEZGQ\nQJ4sdvfJwOQgzi0iIkeKuyeLRUTkSEoEIiJxTolARCTOKRGIiMQ5JQIRkThn7sc81FvsmNm/gc0F\n3L0S8G0RhlMSqM7xQXWOfYWtbzV3P/tEG5WIRFAYZrbU3esHHUc0qc7xQXWOfdGqr7qGRETinBKB\niEici4dEkBJ0AAFQneOD6hz7olLfmL9GICIixxcPLQIRETmOmE4EZtbKzNab2edm9mTQ8USCmQ0x\ns+1mtvqwsopm9pGZfRZexszkzGZ2oZnNMrO1ZrbGzLqEy2O5zqea2WIzWxWu83Ph8ovMbFG4zu+H\nh3WPKWZW2sxWmNnE8HpM19nMNplZhpmtNLOl4bKIf7djNhGYWWngDeD3QAJwp5klBBtVRAwFWh1V\n9iQww90vAWaE12PFQeARd68BNAQeCP9dY7nO+4Dm7l4LqA20MrOGQC/gtXCddwJJAcYYKV2AtYet\nx0Odm7l77cNuG434dztmEwHQAPjc3Te6+35gJHBLwDEVOXf/BNhxVPEtwLDw+2HArVENKoLc/Rt3\nXx5+/yOhH4kLiO06u7tnhVfLhF8ONAdGh8tjqs4AZlYFuBEYHF43YrzOeYj4dzuWE8EFwJeHrX8V\nLosHld39Gwj9cALnBBxPRJhZdaAOsIgYr3O4i2QlsB34CPgn8L27HwxvEovf777A40BOeP1XxH6d\nHZhuZsvC87ZDFL7bgUxMEyWWS5lukYoRZlYeGAM87O67Qv9ZjF3ung3UNrNfAuOAGrltFt2oIsfM\nbgK2u/syM2t6qDiXTWOmzmFXu/tWMzsH+MjM1kXjpLHcIvgKuPCw9SrA1oBiibZ/mdl5AOHl9oDj\nKVJmVoZQEkhz97Hh4piu8yHu/j2QTuj6yC/N7NB/5mLt+3018Acz20SoW7c5oRZCLNcZd98aXm4n\nlPAbEIXvdiwngiXAJeG7DMoCrYEJAccULROAduH37YDxAcZSpML9xKnAWnd/9bCPYrnOZ4dbApjZ\nacB1hK6NzAJuD28WU3V296fcvYq7Vyf0b3emu7chhutsZqebWYVD74EWwGqi8N2O6QfKzOwGQv+L\nKA0McfceAYdU5MxsBNCU0CiF/wL+BvwDGAVUBbYAd7j70ReUSyQzawzMATL4T9/x04SuE8RqnS8n\ndJGwNKH/vI1y9+fN7NeE/rdcEVgBtHX3fcFFGhnhrqFH3f2mWK5zuG7jwqunAO+5ew8z+xUR/m7H\ndCIQEZETi+WuIRERyQclAhGROKdEICIS55QIRETinBKBiEicUyIQOQlmlhVenm9mo0+w7cNmVu6w\n9cmHngcQKU50+6jEPTMrHR7CIT/bZrl7+Xxuuwmo7+7fFiY+kUhTi0BimplVN7N1ZjbMzD41s9Fm\nVi487vtfzWwucIeZ/cbMpoYH+5pjZpeG97/IzBaY2RIze+Go464Ovy9tZv8XHkf+UzN70MweAs4H\nZpnZrPB2m8ysUvh9VzNbHX49fNgx15rZ2+F5B6aHnyQWiSglAokHvwVS3P1yYBdwf7h8r7s3dveR\nhOaGfdDd6wGPAgPD27wOvOnuVwDb8jh+MnARUCd8jjR370doHJxm7t7s8I3NrB5wD3AloTGDOphZ\nnfDHlwBvuHtN4HvgT4Wsu8gJKRFIPPjS3eeF3w8HGoffvw8/j2TaCPggPNTzIOC88DZXAyPC79/N\n4/jXAW8dGh45H4//NwbGufvu8DwDY4Em4c++cPeV4ffLgOonrp5I4cTyMNQihxx9IezQ+u7wshSh\nce5r53P/o1k+tjl6+7wcPm5ONqCuIYk4tQgkHlQ1s6vC7+8E5h7+obvvAr4wszsgNMKpmdUKfzyP\n0OiXAG3yOP504L5DwyObWcVw+Y9AhVy2/wS4NXyt4nTgj4QG0hMJhBKBxIO1QDsz+5TQqJVv5rJN\nGyDJzFYBa/jPtKZdCM2LvAQ4M4/jDyY0KuSn4f3/O1yeAkw5dLH4kPBUm0OBxYRGTR3s7isKWDeR\nQtPtoxLTwtNZTnT3xIBDESm21CIQEYlzahGIiMQ5tQhEROKcEoGISJxTIhARiXNKBCIicU6JQEQk\nzikRiIjEuf8HdH12aZGujm0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x124989be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize\n",
    "plt.figure()\n",
    "plt.xlabel('prediction')\n",
    "plt.ylabel('target')\n",
    "plt.grid(True)\n",
    "\n",
    "# auxiliary line \n",
    "plt.plot([0, 50], [0, 50], 'k-')\n",
    "\n",
    "plt.plot(huber_predicted, data.endog, 'o', color='blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare to OLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Result ---\n",
      "                    Robust linear Model Regression Results                    \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                   21\n",
      "Model:                            RLM   Df Residuals:                       17\n",
      "Method:                          IRLS   Df Model:                            3\n",
      "Norm:                    LeastSquares                                         \n",
      "Scale Est.:                       mad                                         \n",
      "Cov Type:                          H1                                         \n",
      "Date:                Thu, 02 May 2019                                         \n",
      "Time:                        01:46:20                                         \n",
      "No. Iterations:                     2                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "var_0        -39.9197     11.896     -3.356      0.001     -63.235     -16.604\n",
      "var_1          0.7156      0.135      5.307      0.000       0.451       0.980\n",
      "var_2          1.2953      0.368      3.520      0.000       0.574       2.017\n",
      "var_3         -0.1521      0.156     -0.973      0.330      -0.458       0.154\n",
      "==============================================================================\n",
      "\n",
      "If the model instance has been used for another fit with different fit\n",
      "parameters, then the fit options might not be the correct ones anymore .\n"
     ]
    }
   ],
   "source": [
    "ols_model = sm.RLM(data.endog, data.exog, M=sm.robust.norms.LeastSquares()) # set M to be the criterion of OLS\n",
    "ols_results = ols_model.fit()\n",
    "\n",
    "print('--- Result ---')\n",
    "print(ols_results.summary(yname='y', xname=['var_%d' % i for i in range(len(ols_results.params))]))\n",
    "\n",
    "ols_predicted = ols_results.predict(data.exog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x126d22f60>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4k1Xa+PHvXRa1bIIgLtiWUd6R\nUi2biAgKiIA4joyjDrz1J2qgiCD44ihix3EtAo7sIFRgQAkwrILsCi2r7FDZ0YG2IiCiIJZKoe39\n+yOpFEihlCZpk/tzXbmSnDzLfa6G3JxznuccUVWMMcYErxB/B2CMMca/LBEYY0yQs0RgjDFBzhKB\nMcYEOUsExhgT5CwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+RK+zuAgqhatapGREQUat+TJ09Srly5\nog2omLM6Bwerc+C70vpu2rTpqKpWu9R2JSIRREREsHHjxkLtm5SURPPmzYs2oGLO6hwcrM6B70rr\nKyKpBdnOuoaMMSbIWSIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyHn1qiERSQF+BbKBLFVtKCJVgP8A\nEUAK8KSqHvNmHMYYY/LnixZBC1Wtq6oN3e9fA5aqai1gqfu9McYYP/FH19CjwET364lAez/EYIwx\nxVpaWhojRowgKyvL6+cSb65ZLCL7gWOAAmNUNUFEjqvqtXm2OaaqlT3sGwvEAlSvXr3B1KlTCxVD\neno65cuXL9S+JZXVOThYnQNTTk4Oc+fOJSEhgZycHIYOHcof//jHQh2rRYsWm/L0xuRPVb32AG5y\nP18PJAP3AcfP2+bYpY7ToEEDLazExMRC71tSWZ2Dg9U58OzZs0ebNWumgD744IM6ZcqUKzoesFEL\n8Fvt1a4hVT3ofj4CzAYaAT+IyI0A7ucj3ozBGGOKu6ysLAYMGMCdd97Jtm3b+Pe//83ixYu54YYb\nfHJ+ryUCESknIhVyXwOtge3AXKCTe7NOwBxvxWCMMcXd1q1bufvuu3nttdd4+OGH2bVrF8888wwi\n4rMYvHn5aHVgtrsypYHJqrpIRDYA00TEAaQBT3gxBmOMKZZOnTrFu+++y4ABA6hatSozZszgr3/9\nq19i8VoiUNV9QLSH8p+AB7x1XmOMKe7WrFmDw+Fg9+7ddOrUiUGDBlGlShW/xWN3FhtjjI+kp6fT\ns2dPmjZtSkZGBosWLWLChAl+TQJgicAYY3xiyZIlREVFMWLECHr06MH27dtp06aNv8MCLBEYY4xX\n/fzzzzz77LO0adOGq6++mpUrVzJs2DAqVKjg79B+Z4nAGGO8ZObMmURGRvLpp5/y+uuvs3XrVu69\n915/h3WBErFUpTHGlCSHDx+mR48ezJw5k3r16rFo0SLq1q3r77DyZS0CY4wpIqrKhAkTiIyMZN68\nebz//vusW7euWCcBsBaBMcYUiZSUFLp27cqSJUto2rQpY8eOLfQcQb5mLQJjjLkCOTk5DB8+nKio\nKNasWcPIkSNZvnx5iUkCYC0CY4wptF27dtG5c2fWrFlD27ZtGT16NOHh4f4O67JZi8AYYy7TmTNn\n6NevH3Xr1mX37t188sknLFiwoEQmAbAWgTHGXJbNmzfjcDjYunUrTzzxBMOHD6d69er+DuuKWIvA\nGGMK4LfffqNv3740atSIw4cPM2vWLKZNm1bikwBYi8AYYy5p1apVOBwO9u7di8Ph4IMPPqBy5QsW\nViyxrEVgjDH5+PXXX+nRowfNmjXj9OnTfPHFF4wdOzagkgBYIjDGGI8WLlxInTp1GDVqFC+99BLb\nt2+nVatW/g7LKywRGGNMHj/99BNPP/007dq1o3z58qxevZrBgwdTrlw5f4fmNZYIjDEG1/QQ06dP\nJzIykilTpvDGG2+wZcsW7rnnHn+H5nU2WGyMCXqHDh3ihRde4LPPPqNBgwYsWbKE6OgLFlgMWNYi\nMMYELVVl/Pjx1K5dm0WLFjFw4EDWrl0bVEkArEVgjAlS+/btIzY2lqVLl3LfffcxduxYatWq5e+w\n/MJaBMaYoJKdnc2QIUO44447WL9+PR999BGJiYlBmwTAWgTGmCCyc+dOHA4Ha9eupV27dowePZpb\nbrnF32H5nbUIjDEB7/Tp07z77rvUq1ePb775hkmTJjFv3jxLAm6WCIwx/rffSbozghxnCClDI+j5\nqBOns2gOvXHjRu666y7++c9/8thjj7Fz505iYmIQkaI5QQCwRGCM8a/9TrLWxFJeUgkRJaJaKu+3\nj+XLsVeWDDIyMnj11Ve5++67OXr0KHPmzGHKlClcf/31RRd7gLBEYIzxr+Q4SkvGOUXlrsrgzfZx\nxMUV7pDLly8nOjqaDz74AIfDwc6dO/nzn/9cBMEGJksExhj/ykjzWBx2XRppnj/K14kTJ+jWrRvN\nmzcnJyeHpUuXkpCQQKVKlYog0MBlicAY41+hYR6L034KI8zzRx7Nnz+fOnXqkJCQQO/evdm2bRst\nW7YsoiADmyUCY4x/RceTpaHnFJ3MDOXtz+KJj7/07kePHuWpp57iT3/6E5UqVWLNmjV8+OGHhIaG\nXnpnA1giMMb4W80YSjdJIF3DyVEh5cdw+n6WQKvOMcTE5L+bqjJ16lRq167NtGnTeOutt9i8eTN3\n332372IPEF6/oUxESgEbge9V9U8iUhOYClQBNgP/T1VPezsOY0wxVjOG8jVdv/oRwLBeF9/8+++/\np1u3bnz++ec0atSIcePGERUV5fUwA5UvWgS9gF153g8ABqtqLeAY4PBBDMaYAKCqfPzxx0RGRvLl\nl1/y4YcfsmbNGksCV8iriUBEagAPA2Pd7wVoCcxwbzIRaO/NGIwxgeG///0vL7/8MrGxsTRo0IBt\n27bRu3dvSpUq5e/QSjxvtwiGAK8COe731wHHVTXL/f4AcLOXYzDGlGDZ2dkMGjSIO+64g71795KQ\nkMDSpUu59dZb/R1awPDaGIGI/Ak4oqqbRKR5brGHTTWf/WOBWIDq1auTlJRUqDjS09MLvW9JZXUO\nDsFQ5/379zNw4EB2795NkyZNiI2NJTw8nOXLl/s7NJ/w2d9YVb3yAN7H9T/+FOAwkAE4gaNAafc2\n9wCLL3WsBg0aaGElJiYWet+SyuocHAK5zpmZmfrmm29qmTJltGrVqjplyhTNyckJ6Dp7cqX1BTZq\nAX6vvdY1pKp9VbWGqkYAHYBlqhoDJAKPuzfrBMzxVgzGmJJn/fr11K9fn7fffpsnn3ySXbt20aFD\nB5skzov8cR9BH6C3iHyLa8xgnB9iMMYUMxkZGbz88svcc889/PLLL8ybN49JkyZRtWpVf4cW8Hyy\nMI2qJgFJ7tf7gEa+OK8xpmRITEykc+fO7Nu3j+eff54BAwZQsWJFf4cVNOzOYmOM3/zyyy/ExsbS\nsmVLQkJCSEpK4qOPPrIk4GOWCIwxfjF37lwiIyMZN24cr7zyCsnJydx///3+DisoWSIwxvjUkSNH\n6NChA48++ijXXXcd69atY+DAgTZJnB9ZIjDG+ISq4nQ6iYyMZPbs2bz77rts3LiRhg0b+ju0oGeJ\nwJhgst8Jn0WgzhAOjIwg5l4nEREU2frA+fnuu+945JFHeOqpp6hVqxZbtmzhH//4B2XLlj0nNm+t\nW2wuzhKBMcFivxPWx0JGKiJKjcqpJHSOpcnNTmJjvZMMcnJyGD16NHXq1CExMZEhQ4awatUqIiMj\nL4jNG+sWm4KxRGBMsEiOg+wL1wbu92QcGRkUen3g/HzzzTe0bNmSbt260ahRI7Zt20avXr08TxLn\nhXWLTcFZIjAmWFxkbWDgstcHzk9WVhYffPABd955J1u3bmXcuHF88cUX/OEPfyhUbEUVl8mfJQJj\ngsVF1gYGLmt94PwkJyfTuHFjXn31Vdq2bcvOnTt57rnnLj09RBGtW2wKxxKBMcEiOh5KXbg28OvT\n4gkNpUDrA+cnMzOTN954g4YNG/Ldd98xbdo0Zs2axU033VTg2K5k3WJzZXwyxYQxphhwLwVJchx6\nMo3vj4fRZ3I8a76PISGBi64PfDFfffUVDoeDXbt28fTTTzNo0CCuu+66y46tNJC+Jo5Q0kg7Gsag\nZfGXXLfYFA1LBMYEk5oxUDMGAWoAzu6FP9TJkyeJi4tj2LBh1KhRgwULFvDQQw9dUWyXs26xKTqW\nCIwxl+3LL7+kS5cupKSk0L17d95//30qVKjg77BMIdkYgTGmwI4dO4bD4eDBBx+kbNmyrFixghEj\nRlgSKOEsERhjCmT27NlERkYyceJEXnvtNZKTk2nWrJm/wzJFwLqGjDEX9cMPP/Diiy8yffp06tat\ny/z586lfv76/wzJFyFoExhiPVJVPPvmE2rVrM2fOHOLj439fRtIEFmsRGGMukJaWRteuXVm0aBFN\nmjRh3Lhx3H777f4Oy3iJtQiMMb/Lyclh5MiR1KlTh5UrVzJs2DBWrlxpSSDAWYvAGAPAnj176Ny5\nM6tWraJ169aMGTOGiIgIf4dlfMBaBMYEuTNnztC/f3+io6PZsWMHEyZMYNGiRZYEgoi1CIwJYlu2\nbMHhcLBlyxb++te/MmLECG644QZ/h2V8zFoExgShU6dOERcXx1133cXBgweZMWMGM2bMsCQQpKxF\nYEyQWb16NQ6Hgz179vDMM8/w4YcfUqVKFX+HZfzIWgTGBIn09HR69uxJs2bNOHXqFIsXL+bf//63\nJQFjicCYYLB48WLq1KnDiBEjePHFF9m+fTutW7f2d1immLBEYEwA+/nnn3nmmWdo27YtoaGhrFy5\nkqFDh1K+fHl/h2aKEUsExgSomTNnEhkZyaRJk4iLi2PLli3ce++9/g7LFEM2WGxMgDl06BA9evRg\n1qxZ1KtXj0WLFlG3bl1/h2WKMWsRGBMgVJUJEyYQGRnJ/Pnz6d+/P+vXr7ckYC7JWgTGBICUlBRi\nY2P54osvaNasGR9//DF//OMf/R2WKSG81iIQkatFZL2IJIvIDhF5211eU0TWicg3IvIfESnrrRiM\nCXTZ2dkMGzaMqKgovvrqK0aOHElSUpIlAXNZvNk1lAm0VNVooC7QVkQaAwOAwapaCzgGOLwYgzEB\nKzU1lfvuu49evXrRrFkzduzYwQsvvEBIiPX4msvjtW+MuqS735ZxPxRoCcxwl08E2nsrBmMC0Zkz\nZ4iPj6dLly7s3r2bTz75hAULFhAWFubv0EwJJarqvYOLlAI2AbcBI4EPgLWqepv781uAhaoa5WHf\nWCAWoHr16g2mTp1aqBjS09OD7pppq3Pg2rt3LwMHDuS///0vTZs25f/+7/+C6s7gYPk757rS+rZo\n0WKTqja85Iaq6vUHcC2QCDQDvs1Tfguw7VL7N2jQQAsrMTGx0PuWVFbnwJORkaF9+vTRUqVK6Q03\n3KCzZ88O+Dp7Emx1vtL6Ahu1AL/RPulMVNXjQBLQGLhWRHKvVqoBHPRFDMaUVCtWrCA6OpoBAwbw\nzDPPsHPnTtq3tx5VU3S8edVQNRG51v36GqAVsAtXy+Bx92adgDneisGYkuzEiRN0796d+++/n6ys\nLL788kvGjh1L5cqV/R2aCTDevI/gRmCie5wgBJimqvNEZCcwVUTeA7YA47wYgzEl0sKFC+natSsH\nDhzgpZde4r333qNcuXL+DssEKG9eNfS1qtZT1TtVNUpV33GX71PVRqp6m6o+oaqZ3orBGACnEyIi\nICTE9ex0Avud8FkETA5xPe93+jXGXD/99BNPP/007dq1o0KFCqxZs4bBgwdbEjBeZRccm4DmdEJs\nLKSmgqrr+cuxTrLWxEJGKqCu5/Wxfk0Gqsq0adOoXbs2U6ZM4Y033mDz5s00btzYbzGZ4HHJRCAi\nF0xX6KnMmOIoLg4yMs4te7N9HKXlvMLsDEiO811geRw8eJC//OUv/O1vfyMsLIxNmzbxzjvvcNVV\nV/klHhN8CtIiGF7AMmOKnbS0C8vCqnooBMjIp9xLVJVx48YRGRnJ4sWL+eCDD1i7di133nmnT+Mw\nJt/BYhG5B2gCVBOR3nk+qgiU8nZgxhSFsDBXd1BeaUfDiKiWeuHGob67M3ffvn106dKFZcuWcf/9\n9zN27Fhuu+02n53fmLwu1iIoC5THlSwq5Hmc4Ozln8YUa/HxEBp6btnbn8WTpecVlgqF6Hivx5Od\nnc2QIUO444472LBhA6NHj2bZsmWWBIxf5dsiUNXlwHIRmaCqqSJSTlVP+jA2Y65YTIzrOS7O1U0U\nFgatOsdQugmuMYGMNFdLIDoeasZ4NZYdO3bgcDhYt24dDz/8MKNHj6ZGjRpePacxBVGQ+whuEpGF\nuFoHYSISDXRV1Re8G5oxRSMm5mxCyFPq9R/+XKdPn6Z///689957VKxYEafTSceOHRERn5zfmEsp\nSCIYArQB5gKoarKI3OfVqIwJEBs2bMDhcLBt2zY6dOjAsGHDqFatmr/DMuYcBbqPQFW/O68o2wux\nGBMwMjIyeOWVV2jcuDE//fQTc+bMYcqUKZYETLFUkBbBdyLSBFD3amI9cc0ZZIzxICkpiS5duvDt\nt98SGxvLwIEDqVSpkr/DMiZfBWkRPA90B24GDuBabay7N4MypiT65ZdfeP7552nRogWqyrJlyxgz\nZowlAVPsXbJFoKpHAd+MqhlTQs2fP5+uXbty6NAhXn75Zd555x1Cz79u1Zhi6pKJQESGeSj+BdeC\nBzaFtAlqP/74Iy+99BKTJ08mKiqKWbNm0ahRI3+HZcxlKUjX0NW4uoO+cT/uBKoADhEZ4sXYjCm2\nVJUpU6YQGRnJ9OnTeeutt9i0aZMlAVMiFWSw+DagpapmAYjIR8AS4EFgmxdjM6ZYOnDgAN26dWPe\nvHk0atSIcePGERV1wbLbxpQYBWkR3AzknQy9HHCTqmYDtpaACRo5OTkkJCRQp04dli5dyqBBg1iz\nZo0lAVPiFaRFMBDYKiJJgAD3Af1EpBzwpRdjM6bY+Pbbb+nSpQtJSUm0aNGCjz/+mFtvvdXfYRlT\nJC6aCMR1D/wSYAHQCFcieF1Vcxecf8W74RnjX1lZWQwdOpQ33niDMmXK8PHHH+NwOGx6CBNQLpoI\nVFVF5DNVbYAtMm+CzLZt23A4HGzYsIFHHnmEjz76iJtvvtnfYRlT5AoyRrBWRO7yeiTGFBOZmZm8\n+eab1K9fn5SUFKZOncqcOXMsCZiAVZAxghZAVxFJBU7i6h5SVbVllEzAWbduHQ6Hgx07dvDUU08x\nePBgqlat6u+wjPGqgiSCh7wehTF+dvLkSd544w2GDBnCzTffzLx583j44Yf9HZYxPlGQKSZSAUTk\nelw3lxkTUJYtW0aXLl3Yt28f3bp1o3///lSsWNHfYRnjM5ccIxCRP4vIN8B+YDmQAiz0clzGeN3x\n48fp0qULDzzwACEhISQlJTFq1ChLAiboFGSw+F2gMbBXVWsCDwCrvRqVMV42Z84cIiMjGT9+PK++\n+ipff/01999/v7/DMsYvCpIIzqjqT0CIiISoaiKuuYeMKXGOHDlChw4daN++PdWqVWPdunUMGDCA\na665xt+hGeM3BRksPi4i5YEVgFNEjgBnvBuWMUVLVXE6nfTq1Yv09HTeffdd+vTpQ5kyZfwdmjF+\nV5BEkAxkAP+Ha12CSrgWsjemRPjuu+94/vnnWbBgAY0bN2bcuHFERkb6Oyxjio0C3UegqjlADjAR\nQES+9mpUxhSBnJwcxowZQ58+fcjOzmbIkCH06NGDUqVK+Ts0Y4qVfBOBiHQDXgBuPe+HvwI2WGyK\nub1799K5c2dWrlxJq1atSEhIoGbNmv4Oy5hi6WKDxZOBR3DNMfRInkcDVX3qUgcWkVtEJFFEdonI\nDhHp5S6vIiJfiMg37ufKRVAPYwDXJHEDBw4kOjqabdu2MX78eJYsWWJJwJiLyDcRqOovqpqiqh1V\nNTXP4+cCHjsLeFlVa+O6/LS7iEQCrwFLVbUWsNT93pgrlpyczG23vUmfPk9y6tRJypX7kVplyyJz\nasLkEPgsAvY7/R2mMcVOQcYICkVVDwGH3K9/FZFduBa5eRRo7t5sIpAE9PFWHCbwZWZmMn78eCZN\nUlTHAK5F4+8Ln0K907EgGa4NM1Jhfazrdc0Y/wRrTDFUkPsIrpiIRAD1gHVAdXeSyE0W1/siBhOY\nvvrqK+rVq8enn35KaOhgcpMAQL8n4yh3Vca5O2RnQHKcb4M0ppgTVfXuCVz3ICwH4lV1logcV9Vr\n83x+TFUvGCcQkVggFqB69eoNpk6dWqjzp6enU758cF3tGgx1/u233xg3bhyzZs3i+uuvp1u3brz9\n9j9RPbtgTPakEELkwu+3Iiy/aZkvw/WKYPg7ny/Y6nyl9W3RosUmVW14yQ1V1WsPoAywGOidp2wP\ncKP79Y3Anksdp0GDBlpYiYmJhd63pAr0Oi9ZskQjIiIU0O7du+uJEyc0MTFRw8NV4exj/5BwVScX\nPmaH+7kGRSPQ/86eBFudr7S+wEYtwG+117qG3MtcjgN2qeqgPB/NBTq5X3fCVj4zBXTs2DGee+45\nWrduTdmyZVmxYgUjRoygQoUKALRrd+7287a0w2OD96Z2HgqNCV7eHCO4F/h/QEsR2ep+tAP6Aw+6\nZzR90P3emIuaPXs2kZGRfPLJJ/Tt25fk5GSaNWt2zjYLFpy7z5/qLcDj0sIHF3goNCZ4efOqoVW4\nVjPz5AFvndcElsOHD/Piiy8yY8YM6taty/z586lfv77HbdPSzn0fVjXN43Zk5FNuTJDyyVVDxlwu\nVWXixIlERkby+eef069fP9avX/97EljldHJgZAQ5zhAOjIzgyPr1VKlydv+OTZzk5OTz9Q4N80EN\njCk5LBGYYic1NZWHHnqIZ555htq1a7N161b69u37+0yhq5xO6p2OpUblVEJEqVE5lYervU3b2q6b\nxTo2cfJx51hKl8q+8OClQiE63pfVMabYs0Rgio2cnBxGjBhBnTp1WLVqFcOHD2flypXcfvvt52wX\ncfzC+wPKXZXBe4+77g/weP8AgJSCRgl2M5kx5/HaGIExl2PPnj04HA5Wr15NmzZtGDNmDOHh4R63\nvelaz338Yde5yvMdG9AcSwLGeGAtAuNXZ86c4f333yc6OpqdO3cyYcIEFi5cmG8SADh43HMff9pP\nrvK0o/mMAdjYgDEeWSIwfrNlyxbuvvtuXn/9dR555BF27txJp06dEI/XfJ6Vcm08JzNDzyk7mRnK\nP2e5+v5fn3bh5zY2YEz+LBEYnzt16hSvv/46d911FwcPHmTmzJlMnz6dG264oUD7N42JYUvZBA4c\nCycnRzhwLJz5P75Jm9gYwsNh6lcx9P0sgXQNBwRCw21swJiLsDEC41OrVq3C4XCwd+9enn32WT78\n8EMqV778JSmaxsTgWjkVagDXJyXRvDnE/P5bf/ZzY8zFWYvA+MSvv/5Kjx49aNasGZmZmSxevJjx\n48cXKgkYY4qWJQLjdYsXLyYqKopRo0bRs2dPtm/fTuvWrf0dljHGzRKB8Zqff/6ZTp060bZtW0JD\nQ1m1ahVDhw4NqmmEjSkJLBEYr5gxYwa1a9dm8uTJxMXFsWXLFpo0aeLvsIwxHlgiMBfldELPR52k\nDHXN65PujLjour+HDh3iscce44knnqBGjRps2LCB9957j6uvvtrjsSMiICTE9ex04jr2ZxHeXWP4\nMs7hMUZjAoxdNWTy5XTCl2OdjHg69vcpG8qTStaaWNcXJ8/lmKrKhAkT6N27N7/99hv9+/fn5Zdf\npnRpz18xpxNiYyHDPRNEaqrrXH8jltLeXGN4v9N1zOxLn8NTjLHuTWPsgiQTQKxFYPIVFwdvtr9w\n3p7Scu66v/v376d169Y899xz3HHHHXz99df06dMn3ySQe+yM86YDerN93NkkkKuo1xhOjjubBC5x\nDk8xZmS4yo0JJJYITL7S0i4+p392djbDhg0jKiqKtWvXMmrUKJKSkvif//mfAh37fD5ZPyC/Y3ko\n9xTjxcqNKaksEQSTfPrG8/aD93zU6RoHmBxC2vAIfvq1isdDnS5zA82aNaNXr17cf//97Nixg27d\nuhESUrCvVJiHaX98MkdQfsfyUO4pxouVG1NSWSIIFrl94xmpgP7eN77K6SQ21tX/3eEeJ++3j6W8\nuLapUTmVCtecIPNM2XMO9duZq+k84gh79uzh008/Zf78+YRd5q9jfDyEnjcd0NufxZOlXp4jKDre\ndcwCnMNTjKGhrnJjAoklgmCRT994xPG43/vBPc3jf3WZM5z4rQIpP7rm9Un5MRxHwvtk3vBXdu3a\nxVNPPXXJSeI8iYmBhAQIDwcR13OrzjGUbpLgmhvIW3ME1YxxHbMA5/AUY0KCDRSbwGNXDQWLfPrG\nb6p0tjy/Pvrryv/M9d2O/v5eRMnJeemKQ4qJ8fSjGuP9yeFqFvwcnmM0JrBYiyDQ7HfS+IcOF14j\nn0/f+Pd55vbPr48+d57/XGFhl98CMMYUX5YIAol7HODq7B/IOw7AfidEX9j/fjIzlLjp8biXAvY4\nj//JzGt4fdrZTnHrIzcm8FgiCCQXu0a+Zgy9pyWc09ffZWwCn66IoWJFzpnH/9jpG8lRSPkRRqzv\nxKrvOlofuTEBzMYIAsklrpEf8XkMw+ee/RXv2MTJ/iERhFVNI6RcGCf+8BrtXsihyrNrgDBuuOEM\n//rXVaR95IPYjTF+Yy2CQHKJa+TzXuHZsYmTjzvHElEtlRBxdSOV3vQiYXo1EAGEcPjwVcTG2vw6\nxgQ6SwSB5BLXyOe9Lt7TpaKhV2XR78l3zimzKRWMCXyWCAKJ+xr5U6Wq4+ka+bzXxed3qWjYdQWf\nasEYExgsEQSamjGsrT4V/jcH2qdccL18TAwsW7aPI+llPe4uouwfEkHHJmf7g1RtCmZjApklgiCS\nnZ3N4MGDiYqK4vVpcEYvTAYiEFEtlY87x56TDHKnYLZkYEzgsUQQJLZv306TJk3o3bs3LVu25J2J\n31KmyXj3VAsXKndVBv2ePHdwwMYLjAlMlggC3OnTp3n77bepX78++/btY/LkyXz++efUqFHD1W3U\nPgXwfKewjRcYExy8lghEZLyIHBGR7XnKqojIFyLyjfu5srfOb2DDhg00aNCAt956i8cff5ydO3fS\nsWPHCyeJy+ey0/OnlgCbgtl1hlkfAAAN/UlEQVSYQOTNFsEEoO15Za8BS1W1FrDU/d4UsVOnTvH3\nv/+dxo0bc+zYMebOncvkyZOpVq2a5x08XHaapaG8/dm5c0nY9BLGBCav3VmsqitEJOK84keB5u7X\nE4EkoI+3YghGSUlJOBwODh48SNeuXRkwYACVKlW6+E65VxYlx7nuQg4No3R0PK2IIXG/e6WyMFcS\nsOkljAk8oqreO7grEcxT1Sj3++Oqem2ez4+pqsfuIRGJBWIBqlev3mDq1KmFiiE9PZ3y5csXat+S\nJD09nTFjxjBv3jxuvPFGXnnlFerVq+fvsHwmWP7OeVmdA9+V1rdFixabVLXhJTdUVa89cM1VsD3P\n++PnfX6sIMdp0KCBFlZiYmKh9y1uJk1SDQ9XFXE9T5qkqvsm6cnJ1TR7EpoyFHXGt9OFCxf6OVLf\nC6S/c0FZnQPfldYX2KgF+I319aRzP4jIjap6SERuBI74+PwlltPpuo4/dzWx1FT4YqyTx7KeJbTM\nGRAIrwrh1ZPYmVOPC4dnjDHGM19fPjoX6OR+3QmY4+Pzl1hxcWeTQK632sdxTZkz5xZmZ/CHX8f6\nLjBjTInnzctHpwBfAX8UkQMi4gD6Aw+KyDfAg+73pgA8Xb+f33xBV2VbQ8sYU3BeSwSq2lFVb1TV\nMqpaQ1XHqepPqvqAqtZyP//srfMXG/udpDsjyHGGkDI0gp6POgs0TYPT6ZrfJyTE9dz5QdfaAdmT\nQn6fCyi/pSUzS11fpFUwxgQ2u7PYm/Y7yVoTS3lxzfkfUS2V99vH8uXYiyeD3PGA1FTXhG9NbnYy\nuMPZtQNy5wLac+g2PF30dbRsY+/VyRgTcCwReFNyHKXl3I79cldl8Gb7uIvO2XP+eICntQPKXZXB\nA3WSOP8mYYCqp9deSdTGmCBjicCb8lk6Muy6tIvO2XP+Z/mNBZQKyfZYbmMExpjLYYnAm/KZwyd3\nzv9VHvqHMjMzqVjx+Dll+Y0FZOeU8lhuYwTGmMthicCbouPJ0tALikUgvGoq9U7HnpMM1q5dS/36\n9fnllxcoVSrz9/LXp8VzMvPc45zMDGX00tgLyikVyr4KnYu2HsaYgGaJwJtqxlC6SQLpGu5xULfc\nVRlEHI/j5MmT9O7dmyZNmnDixAnmz3+KiROvIjzclTTWfB/DlrIJEBqOqnDgWDixYxP4V9Ko38vz\nLk15JLSVz6tqjCm5fH1ncfCpGUP5mjHkOEMQLswGN1VK47Y77mD//v1069aN/v37U7FiReD8Cd5i\ngBgEqAE4u59bfo7UpCKuhDEmkFki8JGDx8OoUTn1gvK0n26idOnSLF++nPvuu88PkRljgp11DflI\nyrWe+vmvYcY3j5KcnGxJwBjjN5YIfKRpTAyrs4aSdrQGOTlC6tGbWXw8nr+PGMk111zj7/CMMUHM\nuoZ8QFVxOp30eq0P6enp/POf7/Lqq68SXqaMv0MzxhhLBN6WlpbG888/z8KFC7nnnnsYN24ctWvX\n9ndYxhjzO+sa8pKcnBxGjRpFnTp1WL58OUOHDmXlypWWBIwxxY61CLxg7969dO7cmZUrV9KqVSsS\nEhKoWbOmv8MyxhiPrEVQhLKyshgwYAB33nkn27ZtY/z48SxZssSSgDGmWLMWQRFJTk7mueeeY/Pm\nzfzlL39h5MiR3Hjjjf4OyxhjLslaBFfo1KlT/OMf/6Bhw4Z8//33zJgxg1mzZlkSMMaUGNYiuAJr\n1qzB4XCwe/duOnXqxKBBg6hSpYq/wzLGmMtiLYJCSE9Pp2fPnjRt2pSMjAwWLVrEhAkTLAkYY0ok\nSwSXacmSJURFRTF8+HC6d+/O9u3badOmjb/DMsaYQrNEUEDHjh3j2WefpU2bNlx99dWsXLmS4cOH\nU6FCBX+HZowxV8QSQQHMmjWLyMhIPv30U/r27cvWrVtp2rSpv8MyxpgiYYPFF3H48GF69OjBzJkz\nqVu3LgsWLKBevXr+DssYY4qUtQg8UFUmTJhAZGQk8+bNo1+/fqxfv96SgDEmIFmL4DwpKSl07dqV\nJUuWcO+99zJ27Fhuv/12f4dljDFeYy0Ct5ycHIYPH05UVBSrV69m+PDhrFixwpKAMSbgWYsA2L17\nN507d2b16tW0adOGMWPGEB4e7u+wjDHGJ4K6RXDmzBn69etHdHQ0O3fuZOLEiSxcuNCSgDEmqARt\ni2Dz5s04HA62bt3K448/zogRI6hevbq/wzLGGJ8LuhbBb7/9Rt++fWnUqBGHDx9m5syZTJ8+3ZKA\nMSZo+SURiEhbEdkjIt+KyGu+Ou+qVauoW7cu/fv35+mnn2bnzp089thjvjq9McYUSz5PBCJSChgJ\nPAREAh1FJNKb5/z111/p0aMHzZo14/Tp0yxZsoTx48dTuXJlb57WGGNKBH+0CBoB36rqPlU9DUwF\nHvXWydavX09UVBSjRo2iV69ebNu2jQcffNBbpzPGmBLHH4PFNwPf5Xl/ALjbGyfq2rUrCQkJ1K5d\nm9WrV3PPPfd44zTGGFOi+SMRiIcyvWAjkVggFqB69eokJSVd/olE+Nvf/sZzzz1HZmZmoY5REqWn\npwdNXXNZnYNDsNXZV/X1RyI4ANyS530N4OD5G6lqApAA0LBhQ23evPlln6h58+YkJSVRmH1LMqtz\ncLA6Bz5f1dcfYwQbgFoiUlNEygIdgLl+iMMYYwx+aBGoapaI9AAWA6WA8aq6w9dxGGOMcfHLncWq\nugBY4I9zG2OMOVfQ3VlsjDHmXJYIjDEmyFkiMMaYIGeJwBhjgpwlAmOMCXKiesFNvcWOiPwIpBZy\n96rA0SIMpySwOgcHq3Pgu9L6hqtqtUttVCISwZUQkY2q2tDfcfiS1Tk4WJ0Dn6/qa11DxhgT5CwR\nGGNMkAuGRJDg7wD8wOocHKzOgc8n9Q34MQJjjDEXFwwtAmOMMRcR0IlARNqKyB4R+VZEXvN3PN4g\nIuNF5IiIbM9TVkVEvhCRb9zPAbM4s4jcIiKJIrJLRHaISC93eSDX+WoRWS8iye46v+0uryki69x1\n/o97WveAIiKlRGSLiMxzvw/oOotIiohsE5GtIrLRXeb173bAJgIRKQWMBB4CIoGOIhLp36i8YgLQ\n9ryy14ClqloLWOp+HyiygJdVtTbQGOju/rsGcp0zgZaqGg3UBdqKSGNgADDYXedjgMOPMXpLL2BX\nnvfBUOcWqlo3z2WjXv9uB2wiABoB36rqPlU9DUwFHvVzTEVOVVcAP59X/Cgw0f16ItDep0F5kaoe\nUtXN7te/4vqRuJnArrOqarr7bRn3Q4GWwAx3eUDVGUBEagAPA2Pd74UAr3M+vP7dDuREcDPwXZ73\nB9xlwaC6qh4C1w8ncL2f4/EKEYkA6gHrCPA6u7tItgJHgC+A/wLHVTXLvUkgfr+HAK8COe731xH4\ndVZgiYhscq/bDj74bvtlYRofEQ9ldolUgBCR8sBM4CVVPeH6z2LgUtVsoK6IXAvMBmp72sy3UXmP\niPwJOKKqm0SkeW6xh00Dps5u96rqQRG5HvhCRHb74qSB3CI4ANyS530N4KCfYvG1H0TkRgD38xE/\nx1OkRKQMriTgVNVZ7uKArnMuVT0OJOEaH7lWRHL/Mxdo3+97gT+LSAqubt2WuFoIgVxnVPWg+/kI\nroTfCB98twM5EWwAarmvMigLdADm+jkmX5kLdHK/7gTM8WMsRcrdTzwO2KWqg/J8FMh1ruZuCSAi\n1wCtcI2NJAKPuzcLqDqral9VraGqEbj+7S5T1RgCuM4iUk5EKuS+BloD2/HBdzugbygTkXa4/hdR\nChivqvF+DqnIicgUoDmuWQp/AN4EPgOmAWFAGvCEqp4/oFwiiUhTYCWwjbN9x6/jGicI1DrfiWuQ\nsBSu/7xNU9V3ROQPuP63XAXYAjylqpn+i9Q73F1Df1fVPwVynd11m+1+WxqYrKrxInIdXv5uB3Qi\nMMYYc2mB3DVkjDGmACwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhzGUQk3f18k4jMuMS2L4lI\naJ73C3LvBzCmOLHLR03QE5FS7ikcCrJtuqqWL+C2KUBDVT16JfEZ423WIjABTUQiRGS3iEwUka9F\nZIaIhLrnff+niKwCnhCRW0VkkXuyr5Uicrt7/5oi8pWIbBCRd8877nb361Ii8i/3PPJfi8iLItIT\nuAlIFJFE93YpIlLV/bq3iGx3P17Kc8xdIvKxe92BJe47iY3xKksEJhj8EUhQ1TuBE8AL7vJTqtpU\nVafiWhv2RVVtAPwdGOXeZijwkareBRzO5/ixQE2gnvscTlUdhmsenBaq2iLvxiLSAHgWuBvXnEFd\nRKSe++NawEhVrQMcB/56hXU35pIsEZhg8J2qrna/ngQ0db/+D/w+k2kTYLp7qucxwI3ube4Fprhf\nf5rP8VsBo3OnRy7A7f9NgdmqetK9zsAsoJn7s/2qutX9ehMQcenqGXNlAnkaamNynT8Qlvv+pPs5\nBNc893ULuP/5pADbnL99fvLOm5MNWNeQ8TprEZhgECYi97hfdwRW5f1QVU8A+0XkCXDNcCoi0e6P\nV+Oa/RIgJp/jLwGez50eWUSquMt/BSp42H4F0N49VlEO+AuuifSM8QtLBCYY7AI6icjXuGat/MjD\nNjGAQ0SSgR2cXda0F651kTcAlfI5/lhcs0J+7d7/f93lCcDC3MHiXO6lNicA63HNmjpWVbcUsm7G\nXDG7fNQENPdylvNUNcrPoRhTbFmLwBhjgpy1CIwxJshZi8AYY4KcJQJjjAlylgiMMSbIWSIwxpgg\nZ4nAGGOCnCUCY4wJcv8fUI4G7pdo1zcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x126ccd7b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize\n",
    "plt.figure()\n",
    "plt.xlabel('prediction')\n",
    "plt.ylabel('target')\n",
    "plt.grid(True)\n",
    "\n",
    "# auxiliary line \n",
    "plt.plot([0, 50], [0, 50], 'k-')\n",
    "\n",
    "plt.plot(huber_predicted, data.endog, 'o', color='blue')\n",
    "plt.plot(ols_predicted, data.endog, 'o', color='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE(Robust LR) :  9.073643089346657\n",
      "MSE(OLS) :  8.515712457064692\n"
     ]
    }
   ],
   "source": [
    "# Mean Squared Error of Robust Linear Regression(Huber weight function)\n",
    "huber_mse = st.eval_measures.mse(data.endog, huber_predicted)\n",
    "print('MSE(Robust LR) : ', huber_mse)\n",
    "\n",
    "# Mean Squared Error of OLS\n",
    "ols_mse = st.eval_measures.mse(data.endog, ols_predicted)\n",
    "print('MSE(OLS) : ', ols_mse)\n",
    "\n",
    "# By only checking how training data fit to the model, OLS is more fittable than RLS. Estimation using test data should be done to see its effect."
   ]
  }
 ],
 "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.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}