standard_assessment_metrics_function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining a model evaluation metrics function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-15T14:59:15.375466Z",
     "start_time": "2019-11-15T14:59:15.358761Z"
    }
   },
   "outputs": [],
   "source": [
    "def model_eval(obs,pred):\n",
    "\n",
    "    import numpy as np\n",
    "    from sklearn import linear_model\n",
    "    \n",
    "    if obs.shape != pred.shape:\n",
    "        print('flattening arrays')\n",
    "        obs = obs.flatten()\n",
    "        pred = pred.flatten()\n",
    "        \n",
    "#     from sklearn import metrics\n",
    "    MAE = 1./len(obs) * (np.sum( np.abs(obs - pred) ))\n",
    "#     MAE = metrics.mean_absolute_error(obs,pred)\n",
    "    MSE = 1./len(obs) * (np.sum( (obs - pred)**2 ))\n",
    "#     MSE = metrics.mean_squared_error(obs,pred)\n",
    "    RMSE = np.sqrt(MSE)\n",
    "    RAE = (np.sum( np.abs(obs - pred) )) / (np.sum( np.abs(obs - np.mean(obs)) ))\n",
    "    RSE = (np.sum( (obs - pred)**2 )) / (np.sum( (obs - np.mean(obs))**2 ))\n",
    "    RSQ = 1-RSE\n",
    "#     RSQ = metrics.r2_score(obs,pred)\n",
    "#     NOTE: Prediction efficiency and R-squared are the same thing\n",
    "#     PE = 1 - ( (np.sum(pred - obs))**2 ) / ( (np.sum(obs - np.nanmean(obs)))**2 )\n",
    "    regr = linear_model.LinearRegression()\n",
    "    regr.fit(obs.reshape(-1, 1),pred)\n",
    "    lin_slope = regr.coef_[0]\n",
    "    lin_intercept = regr.intercept_\n",
    "    \n",
    "\n",
    "    print('Evaluation: \\n')\n",
    "    print('   MAE  = {}'.format(MAE))\n",
    "    print('   MSE  = {}'.format(MSE))\n",
    "    print('   RMSE = {}'.format(RMSE))\n",
    "    print('   RAE  = {}'.format(RAE))\n",
    "    print('   RSE  = {}'.format(RSE))\n",
    "    print('   RSQ  = {}'.format(RSQ))\n",
    "#     print('   PE  = {}'.format(PE))\n",
    "    print('   linear slope  = {}'.format(lin_slope))\n",
    "    print('   linear intercept  = {}'.format(lin_intercept))\n",
    "    \n",
    "    return MAE, MSE, RMSE, RAE, RSE, RSQ, lin_slope, lin_intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-15T14:54:39.312545Z",
     "start_time": "2019-11-15T14:54:39.111968Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients: \n",
      " [938.23786125]\n",
      "Mean squared error: 2548.07\n",
      "Variance score: 0.47\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX98FPWd/1+zmyXZBWGDxCqrSbxaoUaECP64i1VD74oVsRG11OZbrz2v+dZvbSvV1HjYgj1ywEWL16PtHWftnW2uRUVTFQtej3ienNqCAWk8aGsl0QU1CAsNu0k2yXz/CLPMzn4+n/nM7Mzu7O77+Xi0D9nMzH52k3nN+/P+qaiqCoIgCCL/+PK9AIIgCGICEmSCIAiPQIJMEAThEUiQCYIgPAIJMkEQhEcgQSYIgvAIJMgEQRAegQSZIAjCI5AgEwRBeIQyKwfPmDFDra2tdWkpBEEQxcmuXbsOq6paZXacJUGura3Fzp077a+KIAiiBFEUpU/mOHJZEARBeAQSZIIgCI9AgkwQBOERSJAJgiA8AgkyQRCER7CUZUEQBOFlunqi6Ni2HwdjCcwMB9G6aBaa6iP5XpY0JMgEQRQFXT1R3PvkXiSSYwCAaCyBe5/cCwAFI8rksiAIoijo2LY/JcYaieQYOrbtz9OKrEMWMkEQlvCqW+BgLGHpdS9CFjJBENJoboFoLAEVp9wCXT3RfC8NM8NBS697ERJkgiCk8bJboHXRLAQD/rTXggE/WhfNytOKrEMuC4IoQey6HbzsFtDW70V3iiwkyARRYmSTjTAzHESUIb5ecQs01UcKSoCNkMuCIFykqyeKhrXbcW7bFjSs3e4JX2s2bodicAt4GbKQCcIlvJoXm43boRjcAl6GBJkgXEJkieZTwLJ1OxS6W8DLkMuCIFzCqwEwcjt4FxJkgnAJr+bFNtVHsGbpHETCQSgAIuEg1iydQ1avByCXBUG4ROuiWWk+ZMA7lii5HbwJCTJBuAQFwAirkCAThIvk0xL1as8Jgg8JMkEUIV5NuSPEUFCPIIoQL/ecIPiQhUwUNaW6bfdqyh0hhixkomjxcqtIt/Fqyh0hhgSZKFpKedtOxR+FCbksCCkKcetfytt2SrkrTEiQCVO8GrE3e0hY6dngpQdOZ2cnVqxYgf7+flRXV6O9vR3Nzc2Wr0PFH4UHuSwIU7y49ZfxD8tu273ka+7s7ERLSwv6+vqgqir6+vrQ0tKCzs7OnK/FDC+2Fi10SJAJU7y49Zd5SMj2bPDCA0cTt7/8f8sRj8fTfhaPx7FixYqcrUUGLz3EiglyWRCmeHFKhOxDQmbbnu8Hjt4lNHb8MPOY/v7+nKxFFq+2Fi10yEImTPFixN7JtK58p4jpxc0/dQbzmOrq6pysRZZ8P8SKFRJkwpRctWu04pN08iHROLvK0utOoxex8JW3QikrT/t5KBRCe3t7TtYiS74fYsUKuSwIKdyO2FvN5HAyrat734Cl151G7xKaUtcIAIi9+CjGjh9GTY39LAs38XJr0UKGBJnwBHZ8kk49JPK9/TaK25S6RlTN+3NPN42nPGd3IEEmHMVuPm8+RTHfQctCFTfKc3YeEmTCMbIpIOGJok9RcG7bFldFyontd7aFJSRuBEBBPcJBssnnZQXpAGBMVV3Lc9WCiMs37UZFwIdwMGAraEk5uYRTkIVcpNix2LK18rJxOxi37T5FwZiqph3jZJ6r0Zo/Gk8iGPBj/bJ5lq9PObmEU5CFXITYsdicsPKyTYVqqo9gR9tCvLV2McYNYqzhlE/Zyeq8fAcFieKBBLkIsSM2TgiUk7nBbue5Oimi4VCA+Trl5BJWIUEuQuyIjRMC5WQBidvVgU4JfldPFINDoxmvB/wK5eQSliEfchFiJ43LqdQvp7IF3E4Fc6qwoWPbfiTHM90rkyeVkf+YsAwJchFiR2y8WHnlZiqYU4LP20EcSySzXiNRepAgFyFmYiPKpshHcUK+msM7Ifj5LiohigtF5USzWSxYsEDduXOni8shZMhGwIzpXsCEJZyLMl3WugHkbD1upALm8/vMNV6aqlJoKIqyS1XVBabHkSAXFtkKQMPa7UyLLhIOYkfbQkfXqoe37vIyH2KM7b3T67Hzvcme47RQeVH4SunB4waygkxZFgVGtulp2WZT2B3bw1s3S4ytrEcWt1IB3RBjL1b9eWGqSilAPuQCI1tBzcbnyepV0frEHqx6uhfHEkmhIFkVWKd9sLLfm15geXtH7Rxe746dfUew5fVDOBqfeNiEgwGsur5OSqhlq/5ybUVT8UtuIAu5wMg2fzab/F6WWCTHVMQSSVNrjre+ylAgJ9NIZL43o3Vqdi2eeP7klf6UGANALJFE6+N7pKxcGeG7r2svlm/anVMrmhrS5wYS5AIj24KJbIo3ZKwh1ja2qyeK+Ehm8UQw4MfKJXWuTCMxulYaZ1cJv7eunijuemxPhsCy1qydY8U6TI6rUtt7M+Hr6omi85X+jAeG2+4DL47xKkbIZVFgOJGeZjfdi+fuMKIXKlYwCMjcxju53Wa5EjbviuLG+RF07xvI+N60443NjPQoQMZ3Lft9aMgIuFk+eMe2/aauFDco1J7NhQYJcgGSr965LLFgobfyWNt6AJhcLl/JZtVfynMldO8bYGZu8Naowcv4YH0fCsAVTJntvZnwiUTXeH2n/czUs9l9SJAJaYxiEQ4FMDg0mlY6rFlzmhjwLEgrWR1Wm95bDUCJ1qKAP+yUJZ6Ns6uw6ddvIzmWLssBn3xvC5Hw8axyBUi7fjbDAoj8QYJMWMIoFrLFHkZkg0F2eg1bzSQRuR5UAJt3RbGgZjp32Krx9QU103H/M722sizM4FnlzZdXp12fejQXJiTIRFawBKlh7XahGPt9Co6cGEZt2xYAE5kWK5ewBcuOVWu1L4eZK8aqkGW7tRe5GmR9uZSmVpiQIJcgbuewmt30Y+MqEjo3x9F4Eq1P7AGQuZ22kzdtNQClPz5bF0u2iFwN2hq1zySabkI9NgoTEuQSIxe+RavZB8BEPjPLCrXbhc6qlaodP+/+55nVg9OC7Cb0TiHyuSeSY7j/mV4MJcelf29e7N5HmEOCXGLI+BaztaBlszGMsKzQXKdbKYr8607tNHipgXr0hSYaIldKLr83L/beKFRIkEsMM98iy4Jevmk3dvYdweqmOVLvIeMCYMHbTmfjk7UqFjGG8LFed3KnYZZ2J0LkSuF9b04KKGVzOAtV6pUIWuWaWY4sSxxUAJ2v9FsqzdUGlj60bB44Rmcabow8stOoR7ZE2G6zHVZzJjP/dDDgR5jjMrEzcsrJ5kXUdMhZSJBLAP1NyEKmHFgFbN1kTfURNF9eLRTlylAAHTfNddyisiMWsiXCdrIYeGIo8k9rpeSrrq9zpHTZaQGlbA5nIZdFjsinn020JY5YKAe2e5OtbpqDBTXTc/757YiFrO/VThYDTwwrAj4EA36pXsNujZyy+7ulbA5nIUFm4FaPWzt+NifWwrvZFCCjJLh10Sws37Sb6drI5ibLR9mtXbGQWSuvQCMaS6Bh7XZLucGxeBLrl80z/T17ceQUZXM4C7ksDLjRIDwbf6MTa7HSOpHnYgj4FZwYHrXcmD6fuNmhTN81D0jvYcH7PYl+D5rP/a21i7GjbaFrDy+nv5NsugcSmZAgG3AjSGF3m+jUWlg3IQDER0aZwrq6aQ7WL5uXuskqQwFAhVTfYxZ2p4xkex23xUIT0Ug4KNUO0wstLN34TnL1MCkFaKaegXPbtjC36wqAt9YutnVNu3PsnFxLV08Uq57uzSh6kJmLls0cPtn2mzLrN14n4FMwpaIMsbh4WonbWPk9Uc5uaSI7U498yAbcCFLY9bM5uZam+gg6tu3PEGRRcYETHdt4AcVYImkpX5U5rWRcTRVMOJ3/akU4rfyeqIUlIYJcFgbc2Fba3SY6vRYrrhOzVDlA7sEgEm2R+8XonpApMHEq/7WrJ4rWx/ek+e5FI5jy6Ypwyh1EeAOykA24VXJqxzISrcXO1teKJWdWPSYSHP3afIoinMQhehjos1JEjd/NrmeVVU/3pvV4Bias8eWbdgPItMDzNU2DquTs41XXEQkyAy9tK1lrsXsjWnGdiITNmLssWptIjAH5h4EK8TQO0fWswmoupK2B9z3n42+Geh7bw8sPMnJZFCB2sy+suE54wqYF8nh/uDzL2sco1bP6MFBPvr+CiaBgwJ9+0Vy4CbxUFkxVcvbwcrk3WcgFSDY3oqwlZzcQyVvDuDqRPiczRYPnWvErSsZkEje2nfp1svCK4DkV9PXq9t0tvPwgI0EuQHJRrirjF9XfyOFQAKrKdykoSG8hOTw6zn1vXvtOzf2hbTHXLJ1jmnZnh5VL6tD6xJ6MuXgaXikLdqJKzsvbd7fwcrk3uSwKkNZFsxAw+ACsDNGURZTwb6wiPBpPcn2vLN+vaItodK34Gc2I3dxiNtVH0HHT3ImCGAN23CKdnZ2ora2Fz+dDbW0tOjs7ucfG3h3C7kmXAoqCg2Xn4Pc/7xWuM9sij1xt372UDeKFAh0eZCEXKkaNkulx6SCyPXwjgmZFor4PetfKuSdn7xkx66yWzTZce/9sr9PZ2YmWlhbE43EAQF9fH1paWgAAzc3NqeOGBkfxnbn/hs/9YRXm4R0AwMyxd/DCw1tx3qfqTNdpl1xs371mhecrK0YGqtQrQLKpnNOTjdjwqtP0aJVqZnnEZtWCVj8vq6pP/x659JnW1tair68v4/WamhocOHAA46Pj+O6VT+Cal7+F2ci0St/fcwhnXHSmK2sDnPtbyvd7eB3ZSj1yWRQgTlg1Zo2LzLaYMv427RheLw0Np3oUa4i24W40jxLR39/PfL2vrw8vXNIKX8CPO19elibG7+EMbJzehsFDf3RVjIHcbN+9HETzGiTIBYiV7m08shUtM5HV39TGzmgszHoUW/GVigQg1ylP1dXVGa9dipuhArh65wNprx/DVKwKrIby5pto+WANppw5xZU16clFtzYn/l5LBXJZFCBmW3IZRA1xeFFo4xaTlWVxLDHR6KdxdhW69w1kuAV429fKUAChSWWOuBFEW+SDJx8yrM9tt3mUCL0P+QK0YDXexw3oyjjuv3Alzn71SXz40tMdX0O+ceLvtdCh5kJFjBNBCVHqj+wWUzREkxfEYaVqBfwKBodGHWsUJEoH4zVLcstaa25uxsgj+/GF7X+LMTwMPzLT/XofeRVXfeFSV97fC3g5iOY1yEIuUURWC0+0ZIMwZkEcY1DtxPAoM2XOrygYV1VbNzAvcJdLa63nwe2ov/vjzJ/9DMtw1sZv46ovnu/oexLehCxkQoiZ1WK34KCrJ2rartNoWfPS2oyFIPp1AxPugBUrVqC/vx/V1dVob29PpZLxeoBoPmT/yaZHor4cMuhblGrXnH14DFt/+CnUM47/Ba7BoTv+Dn/1j6yfEqUOCXIBYzd9y3je+mXz0s6zu8XUrE8eouCOWXtNY9Mc2fxe49r0TY+0h0w2Yqy/ZsXAMJ7/4R2I4CDz+P8+Yyk++d5mW+9FlAbksihQ7G693Zy8Ico3Fq2NN1XEiD7wZpbfK7s2WTcM6+GnWcYVgyN46nt/g49iH/Pcd31n4UPJKBRWhyWiJCCXRZFjlr7Fs27dnLwhSl0TPSiMFjmvh7Lewubl9/Jet5oLa8wgGRwaTfVI1r6j4UQSBx68nnk+AGxHIxqO/QJnTi3nHkMQekpSkL3U3cruWnhCookFr0xVJhnfbk9dUbqc2bX0Pl+e9a/3YVdXVzMtZFber2htLDeK8f2Nnd/UcRUtD2zB8qENzPfqxQVYens7ZlRPxQ4SY8ICJVcYkutKLbfWEmY0vgEm+g6LLGfZ9C47VVRWp1vzkClWaG9vRygUSjsvFAqhvb1dem28QKWoT8en/3Un+jqWcMV47l8+gcX3/D18p1d4olkNUViUnIXspSkL2ayF5/of57yuCSyvtaURGeFmWfdrls7JmG59NG5toKl2nOhYLXDHy7JgXQ+QC1SyHka3d/4C97zzPe56rrl2I353UcSRzA2idCk5QfZSXb1oLWaujGOcVpc8NIE1CtO0YAAnRkbTev/KpLjxij/WLJ2DyeVllqZb26W5uZkrwCxkO6Pp3RvTX1Hw2n/xK/hefOBXuPKuS7BVehUEwafkBNlLzal5a5kWDJi2K+SdGw4GMDw6LvS/GoXJjh9bZN07+dDLh7+/ddEsfPtv9+LzTz2HL+GfmMf8w82P42uP3YQrXV0JUWqUXNpbvurqWcICsAswKgI+5gghfYqW6HMASMsQ0PeYcErQnOiFYUY+fle/faoX+5bei4XYjik4kfHzF+qX4+rXvuPKexPFC7Xf5JCL7lZGeME7AMy1xDjz3PQWpuhzaJM+1i+bh6HkOGKJpOMBTFGRh1MtHXPZme3d1w4CioLzl16I6/FMmhi/hAZ8a/4WQFVJjAlXKTmXBZD7ke0iYWFNcJZtgCP6HF09Udz12J6MfF6jL9euS0DUwMeJZjIyJdhm58u8/0B/AlU1IbC6Du/BRXio6u/wyLvX4goq6iByQMlZyPnAqk81WwtTs8hZxRX6980m7U5koWfr97Vbgm08X/S5huNj2KzciKqaEPMar+JSfDTRgx+9v5gq7IicUXI+ZDcwEyBRSTEvRUrftMannEpnCwcDWHV9nVDgzEYmab5c2XJiKwLrhN/Xbgm22fmRcBAvfaMRil9shxx76wim1VZKrTVXeKmYibAOlU7nCJkBjqLcX16pcupcwzj6WCKJ1sf3ZByvx2xLH4uPoKsnKmW5Wx1QabekW3b9FYFMMTWKFU/Md9z7ceBe7qWx99k+zFlcjWn8Q/KC14aEEu5BFnKWWLUyRZayMQPBzLLmZSyYWcgATrpEVCSSmQ3TNSvc6noB8fDTYMBvmhUyMxxEfGSUmWViPI/X31gB0tbwfx96AfcOP5BxHY3n/vFNXHvHn3B/nm9oSGjhQ1kWOcLKdI0dbQvB80ayriOyFEU/M5t3B0xYrcOjmWIMAMmx8ZQP1ur78/y7fkVhWs6rnu7N8PcODo0i4Of7bY0Wt/G6mhgvfvQNHFh3HVeMn77+YUBVPS3GgLeKmQh3IUHOEqsDHK28LgpeiX6mBdz8ijgYxSuzPjEyZru0mheQ5AUYY4kks/vc5EllUkNRWaK0bPMOHFh3Hb536BvMc185qwlQVVz/89u41/cSNCS0dCBBzhKrGRFWjm9dNItpKQZ8imnGRVN9BA9+eq7QUjYTbB5m719edurPqjIUMJ04zeJYIokdbQu552lipBelK17cjwPrrsO6369hnvOvwS8BqorLDz4FYMKN1LB2O85t24KGtdvz0mBKBqfyugnvQ4JsA/2N3LFtP26cH5EuNNGniwGntvId2/ZnCEJTfQQdN81Fpa6zWzgYQMfNc6WCOU31Edw4P8J0kwQDftxy2TnMGz0cZHeSS8HRcc2fq+9jMXTSR80TlUpO1zpNaM3EqHXRLNQcGMSBddfhJy/fxbzWM7gOUFV8Pv6DjLV6oeufGfkoZiLyAwX1LOJUOW8uyoJ5kzgqQwGsXFLHzRkGMku6jVgJQvKGm4rKx/XfAy/l6+2d7+GcS1glHaeoaX0GZ08PWV4rQTgJpb25hFPtO3PRBpTX1zc0qSz1HqJqP1GWhZUgJG+4qfG9eOlwxvNOfDAEKArOYV5pgo985SkkQwGEOFt7s7VS3i+RD0iQLWIl4i26qbMpC87VWkXFI7wgpJ1OerKl7Oq4CsXvw2TBMZ2P7cX333wfo7GEsC9xOBRgptaFQwHK+yXyBvmQOXR2dqK2thY+nw+1tbXo7OwEIB/xFvkou3qi3PQ3JyPnTqzVahDSreDTz5VPCSvsfrOlD1BVNN98IXa0LcRbaxcz+4Ro8Dx1qprbpkYEoYcsZAa8EfM7DxzFidHzM45niY7ZTc1rXelk5Lx10Sy0Pr4nNZwTmHgCx0dGcW7blrTpyaLmR9rn0bfzXL5pNzq27U+zQO00FTJ1DZzMBPkU5/wffmIDbtv2ZVxo8bvhNfg/lkhyf0Z5v4TbUFCPAW/EfNm0MxD50iNpr+kDZHpE/YIBtiADwIG1/OkUVunqiWaUXhsxVs/pUQC8pVuP04FI0fV8l30N1yc3c8+9/UMP4Befn41wMIDJ5WWWfb2ioB7Adik5FfAj/3TpQUG9LOCNkh89NpDxmhYgM95kPB/lTJMbHnDuhu3Ytl8oxsCEJexXFGbhhtG1kU0gsrOzMzX/btK0Kpx2xecw9cKFGe974Y/eRdPqT3Kv8w3/Gjx295zUv2OJZCrNzoqvV9Q+dGffEXS+0p/20HTK9UL+aUIECTID3oh5/9QZGa9p8++MNxnP29k4uwoLaqZzxcDJG1Z2iz2mqsw+E0YBslvCa3QBDcfex8jWDVABTKlrBADMeX4Az/R8gXuNGKZh3j2dpp9F9gHBc68AwOZd0TQxVgDcON+ZHtpeGrJLeA8SZAbt7e1pAgIAvkA5wlfemnHszHCQeZOxu0QA3fsGsLops6GOZgU3rN3u2A0r6nymJ6LzJYuscrtZFCtWrEj7LgFAHR1G7MVHMXdkDv7z+c+LF6iqeKEniqDEtGxAbkgswM7uYH3/KiZ+b05AfSkIESTIDFgj5m/44t3YOnQe04pcvmm39LXNcnKdvGFbF83CnSZr0wKJMqlnvDai8ZFRdPVEuT2d+/oyXUBn4xy8ffxtQCDG6th4qjk8y6LldYVjDYldvmk3dvYdST0M9ejdKb7TZiB85a0py13DKcH00pBdwnuQIHNgjZi/imN1iQoojPBuPE28eB5fOzdsU30Eq57uTStlNqJC3hWiHWe85tF4MsOtone9+KfOwNjxCQuzAiEkEAfwNvd9ThwdwX+89T46/r4747s2TstmuX4UBUwrt/OVfiyomZ52DaM7Zez4AI5s3QAAaaLslGCKfNcEUbJZFk5GulnCEPApgIK0oBovI4FX4qyh9fc1my7C+ixm17aTOSBTdqw/ZrC3Gx88twHq+LDwup3/ugvfPxRDNJbI6Gks+u6Mn335pt3cB5vx8/IyavxTq3D27T8SvrddKMui9KAsCwFOR7pFASKzG483jFRDL0ysdXb1RNNyjaOxRNpEkab6yETWwKv9GcUQdi0zGbeK/r8PP/ug8Hq7NryMt//sHKzW/U6M3wbPj85ytVgp+eZl1IwdPwwFcEUwcz1klygcSlKQ3Yh0824ys6II0TBSwFyYVj3dm1b4AUz0E171dG8qHW/zrmiGGMvM5uMh4wedGQ5OjEwS8ORXurH0u1djPoCvMoJpRoyjpXgPO5GVbHQ98DJqamqq03KwCSIXlGTptFci3bzmP2bo18nzD2uv895jcnmZ7YePaYm0ogjF+Kfz1gGqiqXfvTr1msx3r4kpq9S79Yk9mHf/86mWqH/24ekZ5emsHUF7eztCofTJ06FQCO3t7abrIQinKUlB9soEBpEIifoSTzPrVyzxHtFYwnZjdl5/3tGPrUiVOrP4H/wpoKq4pSdzkofZd68XU9ZDJjmmIpZIpgT6tf5jaL68mtlDuKsninn3P4/ati1YsTeMGdd8BVVnRaAoCmpqarBx48aMgC5B5IKSdFmYRbpzFXThbf39ipIa/mnsRQEAJ3RpZpWcikCt8bsoF1nfSAiw5j/Xu2g6L/0HNJm4J6Cq+DPBj1m/E14wU8aaTiTH0L1vICNgafS5A4DykY9h2uwrsfEmucb/BOEWJWkhiyYw5HKSBGvrH/ArOK2iLNW8hzXCKTmmppoUrVxSl3FMwK9g8UVnpTIdzAY12e1k1vWlrYCioPnXd/IPUlV+azUdrN/J+mXzcIDRtU12J8MS7o5t+zMecED6d0oQ+aIkLWSAH4TLZWmrMTsjHApgcGg0rTcDD32Bif4aM8NBNM6uwuZdUW7Gguh6MuzafADzbzoXTYJj9EUdsshmH/AKVIywhNvuJG+CyAVFL8hW3Q+5DvjpRahh7Xam+4GFXmyMQsYq/7VyPR6Hf3cUM86fjvmCYwbfO4EpZ4RMrfJs4D3I9JYvL6VP5MKhajki3xS1y8KO+4EXMLMSSLOLrOgrmPgsvICc1YeHWT7yWHIcUBTMOH8695g3tvYDqoopZ4S4xzhJU30EO9oWYv2yeQhNKkNyXE1N0RYNAW1dNGuiaMdAwG8+yZsg3KaoBdnO5AdekoAgecARunqi8HHepDIUSLXmZBWKGEXZqqUn7GSmKPBP8rN/BuDZ+14BVBUXLBJNuHMH/QMXONW1TrQLaqqPoOPmuWkZLJWhADoooEd4gKIWZFHKV8Pa7cy0rxjHZcB73QzeKCg9ogIRBUi5MSpDAW6hiB5esJAHs5OZogifQo9dtBpQVVz3t5dxj3Ebu6OWmuoj2L3yEziwdjEOrF2Mnm99gsSY8ARF7UPm+Qu1LT+QmfblZDcu3igoAGl5rqICEb01zEN7wERjiVSz+XAwgIqAD7F4MuU753V+S3twmWwFXsPFuFjdhU8DuK9rL3766tsYUyfcBbdcdg6zm1q28OIAufD3U98JIpcUtYXMshSNTWuAdKvKyUGdrD7A8XgcK1asSHstWwHRP2A0KzuWSGIoOY71y+al0sYiooIYE4v4DXwUUFVcrO4CMCHGP3mlP/V+Y6qKn7zSj/u69mb1WYyI4gBuF/jkMgWSIIAiF2RWbisvBUyfRsbLUbYKr3GN8fVsBYT3mYzbd9bD5qvf+Q/TnhNQVVygvpH20k9fZbfP5L1uF5FbghWgC/icC87R9Gki1xS1ywJgp4SZuSSc6sbFa1wz48yZaFi7nZs3DLAteTvorW99utjcjb/G999qE58sKOjgNUQSNUqyg6lbwmjUOxh89UrPE6J0KGoLmYWTLgkzWI1ryiuCKLvss2nb4M27orhxfiTNKm++vDpjnXYwWt9n/vIP2HHvx4VirI6Np4lxV080IwjqF7g37PTH4L2PyC3BGuIqqrhjXV+EV3qeEKVDUTaoNwvEWA3UZBPY0Y8Hqq6uxqTLPouRcxsyjmM1ite/b0XAh0SSN6mPjb6xevQ3RxGZw88jBoBj7yYw7UMVae/NaxZ/cfU07HjziPB6laEAVi72jI3dAAAaXElEQVSRa/HJm/5x4/xIxu5B+1y8FpsKkNE6kzdEYEpFWVrgU2YiiZPN6onSQLZBfdEJstM3kdPXO7dti7SI6OG5WoxoWRZaQ57rLjgTZRViz9TGh3bi3xLHhC4UI5GTx2lZFjxkvyvRFBLeAFbeOZWhAEKTytKOlxmzxVqr7MOYsjEIESUryDLjhbx4Pb+i4MFP84sTeEKuwRQ+kxS2X9z2BIa/fDmzo5wZ+gdIbdsW0+N546c07DyomFavXwFUZJRRy5aS2/m95sKSJsEvbEpmhJPxD1V2dI8sMoEdKzcLrzHOmKoK22CKPluG2JkI8TNn344lb38fnwQw7/7nLYsxYL2U3KzNp538b32QUsvBNvqUgYnMCG3nYIadvxO3G1I5PXKM8C4FHdRj5YnypMhuIMYssGM1V1VLq2MFxUQpVbxg5EO6POMPlNOFYvwyLgdUFUve/n7qNdFEahFaT2bgVO9lM+x8PrNga1N9JHWuSHC1smozpgUDlgJ/gPvZGJR+VzoUtCCz/lBZt6TdLIqunijiI6PC69m9WXjiwbuJRfnRDyh3A4qC08EOsg2hHFBV/Kn6cupzaaJjF7OezDzsfD4zZEZhadfTrh8OBjL7SPsUnBgZtVwI4nY2BqXflQ4F7bIQ/UGGgwEcS7Cj5zKw/ILadfXDQa3eLNp1eZht0fWfY9NZXwPe/S7u5p4BQFVRwXh/MwGbPMmPcGgSDp4UJxasnsz68m0WVj6fLGbCpG84JArYxUdGM9qfyrgezCbQZIuT5fyEtyloQRb5VSeXl2H3yk/YvrbscFCrN4vImpO9ib/3Fw/jy7/8IpYJjqmtqUF7ezv0k+G6eqK467E9pr7UgF9B+w2nrFNeINKnKKlRUiyxc1Ok9Fjyr+swrpm3Y4jGEqnPyUI/AVxzAVUEnNt8ui34hHcoaJeF6A8y2+2crOVrxffZ1RMVpl6ZbdFf/vHvAUXBl3/5Re4xysn/aY2MtO5yoo5y+nMj4WBaK0qe2wY4FYhkbemdLEE3g9ezxCoii1PGdTE8eipP/Gg86Vjfi1x+l0R+KWhB1oZ8ssh2OycTzGtYux3LN+1GRcCHcDAgvFm6eqJofWIP9/0i4SD3Bjv0mw8ARcGf3voR7vmnTa3JECF9IyMzP2skHMRbhvl1moiLpph4IbikFyxArmc0C5awa5h9Tgq8EU5Q0IIMTASU3CiFFlm+xsyKo/EkhkfTO6sZWfHUXmZKlmi9yXgSUBScNWcGd52XLfspau95FoPHxY2MzHYMcV3WhIZMsAw4taXXk+tOadoEEVYDKVlh1ISdh515fE4E3qjrXOlQ8ILs1nZOdF2r1tB9XXtxYoQvbEx/o6IgMHkS95zrG/4Ftfc8i/dqTwMAlIermMfNOHMmAPMdA2uLbUVMjOfmy2LMVhhN25RycDPTgqzv0qGgg3oaTnVnk72ulZu+qyeKzlfY1quGJoYA0HTx2cJj1172CP7tL87KCPDcfvc38YPV38Dw0Kk1KGXlKLvssxPuEolJzcaMgmnBgHSesvHcfKVqOZGRYCeI5mbgjdLeSoeiEORcY+Wm79i2X6qN5v+u/iSwmv/zHYGr0TDSjTYAs5mVgddgy+sH8YetD2Ps+GH4p85A+MpbMWnWVejYtj9VDqydZ5bK1tUTxQlGMM8HgNfiSC8QZt+RvnmRsf9GU33EdqmwE8JonGrNe3/jGm+cH0H3vgHHy5sp7a10KLpeFrnASu8CUQ+Kwd5uHH72QfM3lPwdWekHYdajw6xxj1l/D9F3BIBrrZt1eJPtHOd234dcdoKjrnOFT8n0ssgHshYUwLdu/nr9f+C+kX8Qvk/Dmv9E4+wqdOua2YvExYolZWZJ8rbDsXgSK5fUcc/Vi2E4FEB5mS+jQKdh7Xau6ySRHGN2kLPSG8ItF5Yet/tX6LHy90YUNgUvyPnqgiV70xuF7+Kn38OT/3ub8Jzae56d+I9YAj/R+Z/NmspY2a431Uews+9I2pDSG+ef+kwicecJBJBu+R6NJxEM+LF+2by09Zr5Pq2WleeDXPt1c/GQIfJPQQuyE12w3BZ07Vo/XvUSNj/9GeGxysn/rxEcI7LCrFhSXT1RbN4VTRtSuulXb2PL64cQiycx7WSvB32qnl7cWQLBsnxZ6xVV1gHgll57yWdKfl3CDQpakLPdNjrZ1pAn7Ef6/oimi89Gk+DcCpRjGMMAAP9Ufs6xhsgKk7WkWN9dclxNFYHEEkn4MOEz5k3UkF1XNJYwnSGoIfIhe6lUmMqZCTcoaEHOdtvohB+wqyea1sMAmBCgv3lsD5ouPhuioUkf8n8E74/9Djgpxr5AOcJX3mr6njPDwawte5nvaBwT8UTRJBP9iKpJ06ow5YrPYUpdY9oxCpCyJvUzBLv3DXCzLBbUTPe0z5T8uoQbFLQgZ7ttzFbQeZ3Ttq37Ombht9zzrlr4U/RdchrGertR/tKPMXJsANXV1bjhi3dj69B5wlxhBUDt6cGsLXszt4GGKA+5s7MTLS0tiMfjAIDh2PsY2bYBAFKizJqenUiOoXvfgHAyh1d8pqIHn1fWSBQPBS3I2W4bZXNleRaQ0cI+sO464ft9/oLv4IUl56f+PaWuEafVNaZZoJO79qYF2v6kKoTfv38iJWoqgP958wi3PFhWIGQKRVjov5ND//x1DJ8UYw01OYzBl36M0+oaXZngkktYLq3lm3bjzk27TUdSEYQdClqQs902igRdxr+sicqmdd/GZfgV930evq4LP2o4zdSaZwXa9GKsYVbUIYOxhzGPylCAO4F6ODbAPGfk2EDqIcPLZy6E4JdoAAKNUSLcoCh6WexoW5jRqUz23Gz6VXzmn/fgwLrruGL8tbMeBVQVf/3Mp6TadMpOQOFhVeT0DXl4LL7orFRjG+N6eAHI6urq1H/bHc3kBcwecNRPgnAaT1nI+cgpttOv4olrH8FNv7gNazjXfBWX4tBrT6IRSMsu4JXW6i1QuygAqt77FSoqb8FwbADl4Srcfvc3sX7FHabnioSne98A160RvvJWHNm6AerocOq1UCiE9vb21L9zHfxy8m9Ixs9eCK4XonDwTOm018pDWVvtP//lXjy8617hefX3b0vl8Z4YGc3I49VKh/XVbINDo8LJz6zAmJHB3m4c2bYBavKUOCqBcnxt5QOmoiwqoxb1vdDed1AXmGxvb0dzc7PgDPdw+m9IZtyVvlycIHjIlk57RpBFotC6aFaa1dM4u8qVJi569Ddj9W+P4MWnxOloXa+9IxUkCwcDGB4dlw6maXm52uf1cYom3vnBFzB2PNOnWx4+A0NH30t9Jpb1KBIyM8vdSz0VzPpz2IHnPwe89dkJb1NwvSxERQXG4JqVcmIZeEIVf3sQn/3UbOG5yXgSgWAZOgT9GfTItrMEkCpn1j98eEUVY8cPM6+hBd5kgpS8rb5RrDVh8lqmgVlhip216l1a+SrTJ0oHz1vIognGeuxaQUzr0OfD9jV/hbPwLve8gf89jKrZp6f+Lerqlg3BgD/Dcr1xfgTP7jmUJu48Czkw7Qw81v0ad7ipzPdWKELE+xvSIIuWyBcF57LgbZtlt/asFpOi97r/mV7mrDizXOJ9//4aZt9Sn3atbANyPHgPI5bbY7C3OyPAppSVY8a1X8G0ukahj/qA5PcmS74E3G2fb6E8mAjvUXAuCwAoL/OlbqbKUAArl9RJi51sypc2bNQ43676CRUvvrmEe95/fvMFfPzbV0HvwJARABZ6X+RgbzeOvfgoRnVN5afUNQofRiy3h1YZF3vx0bQG9aGPXm0aMBSNuLeKk/1BrCKTW203KyKfn4soHTyRh6z9seuFZig5MZdCNAlYw0pea8e2/Wli3PDsWziw7jquGG9p+29AVfHxb1/FvBZPNP0KfxC9iokHzonebhzdtgGjxwcAqBg7PoAjWzdg0ls70qYoyzKlrhFn3/4j1NzzDM6+/UcZPSV4a7lz0240rN3uyNDMfM9/M8uttluQku/PRZQGnhBksz92/RDQcDCA/3N5te2hppqFVPfSe3h93S3o7P0K87hH6tsBVcXiNVeYXouF3tUw2NuNd37wBfStW4J3fvAFDPZ2Yyg5DmXXzzCuS1MDAHV0GP3P/xCA3MPIKZyaZJyLUumuniga1m7HuW1buA+SbApSWNenuXZELvCEy0I2wwIAhkfHsaBmOlY38ce1izj7YAW+/uMNuAFdGT8bhR9XXPpz/Pm91VLXlykcMPp2NSsYAD44dJB5znBsAPc+uRdrls7BmqVzMrrJuUW2Ey+6eqLcnGmnSqVlXQd2C1J41+cNfC2EEnCicPCEIPOEza8ojo3J+f2LB9F91Uo8h5/hNAxm/PyrZ/8dnm6+CADSGrWzbmRRbqqR2IuPpgXagAkrOPbio/BPncHMjPBPnYFEcgx3btqNyZP8GBnljRU9ebxkJooM2QQneQNdFcCxUmkrLVPtdGPjXb8i4GNmvBRCCThROHjCZcHbXvJERssrldleD/S+DygKzrsqgi/i4Qwx/kq4A7X3PJsSY+BUo3YVmVt5zYLS93bge4v5+cFjxw/jT675a5RXpFtYSll6T+QTI2PCoFwkHMS4BTEO+BRUhgLcn2tBPkDONaCHt9NR4VzgS7SbcsIHLpolyOt7QhBO4QlB5jX5EQW1zHye8cNxQFFQdeGHMn7WiwvQHHoK46PjePb/ftR0fXp/Nq8BEE+UA9OqmK+XTavC37d9GR+69qvwT60CoMA/tQrTr7lDKhiHk+/ZOLtKetvsVxR03DwXPd/6BB5aNo+5ZhUTn1H/4GE9mFjw1mE1OMlCeziIHj1O+MDDnIdVOBTIqpEVQcjgCUEG2F3bzIJarCj3WHIcP/M3I1Q1mXlOJz6LDw++jptfugQf6+iWLubQrHLell7FhPWpJ+BT8OXWbyIUCqW97guU447WbwIAlI98zHJmhP49N++KonF2lVQmyoOfnpvWXF3UxtNOVoFbnd2MuxIeTmQ98DYbDnmECEKIJ3zIPKzmlarjKnZNuhyfwa+Zx37w1nE0155mK39YP4aIBavnxoTv+VosqK1MjTnSN+BpWLtd+v15aNM39H0njH5tBcDF1dOw6ule3LlpN4CJtLuwIFBlJ6vArc5uovRCK+uT4RgneMp7nSCcxNOCDJwKzMg0Oh+KDeFShhj3/fJ3qPn4eTgdE9YWr4xYhOhozQrkBZGam5uZHdCcSpkyXofV0H7Hm0fSXjsaT8LvUxDwKWk+au2z8B6CZu4RN8YaWfmess16oGnSRD7xvCBryIxrCk4P4oUFd+OSnd/HZMSx5DNP4Dc1FZj56360Tp+4oVofF4uxMZJulkWhWcZAeu9jvWXIK7mVSZvzKYAgpgcACE3yY/mm3ZZ7aYyNq5gaCiA0qUyqqVC+sgp43xOr+1q265MdC0Zl1IQbeKaXhQyyNwGvLwagIpHkp5Dx3A48a1HriyBqXwmwhU37WevjezKyKAJ+BR03zeW2xjQeaywDt4KoB4hXRIf3/fIa/jvxfqLP7bXe3YT3KbjmQk5i1vWLheiG4oliOBjAquv5/TYqQwEcT4wKu6zVf/t5ZpMjvdjrC0MmT/Ij4PfhWGIiR/rE8Khp0YjIyi+UButeeTgA7vRdJoqbgmwuJIPMjWnHN7tm6RyceOMF1N6QGXzTrm/sEBdLJIXWK0tojWuMcY45eDKv1nj9cRVYdX1dak3ntm0Rfq7IyR7Km371NtMSN9vie0UI3fBN24XKqAm38EzamwyyubFWAzCVoQBOvPECWlpa0NfXB1VV0dfXh5aWFnR2dgKYEITQpMznl9VOb8Y18tY6Mxzkpp7duWk3Pnzvc7ivay/3fAXAQ8vmYUfbQqxumoOOm+ciHDyVY1sZCqTcIjzs5CKXAqLfGUFkQ0G5LGS3ilbS2jR/7Z03NKCvry/j5zU1NThw4AAAZ5vQP7RsHtdHLDNDT6Phw9PxWv+xjPObL5frxyG8dpZbc7es63xb7eRDJqxSlC4L2a2iTP4ykD6CaGl/P/OYvr7+1PgfmawIWfQFGvq1WhFjAHjlD0fx4Kfn2hYokbhlszV3q3+wG9e1KvC5nqRNlA4FJchWckTN8peNVt6kaVUYjr2fcZx/6ozUTX/j/EjGPDurAsqiqT6CE2+8gNvWfR3DsYG0RvVmjKmqbf+qmbjpv+/B3u5U8/uyqTNQ1Xsrzm+4litEVpoAWcHp69oVeC/5tInioaB8yHZKc0Xn6JvnTLnic1DKytOO0zf60VfE6XtuNF9ebflzhHT9nbt6ojh/2d/gc5+/7eQD4VSj+sHe7rTzWH2VRY3wzTArj9a+O62F6NjJRvqjxwfwwdYN2L/jOa5P2a3Al9PXpcbzhJcoKAvZ7laxInBqNJSWqgak5wfzRiDprdRoLMG0jLa8fkiYUWEkOaamROzeJ/fiD1sf5rbo1N6f11f5ivNmALhW+r31mHVO0z7nZzb8pXB9LOvUrYo3p69LGROElygoQQasbRVZwZfhk72FWZbRlLpGoZuAN39u5ZI6S70xkuNqygJLJMeELTo1eH2Vdz/1T8D6NuH72akU1G/bR45l9mzWr48lXrIVb1Zx+rpUKk14iYJyWViFtx39+mO7bQXnVAB3PbYnY4uubx8qy8FYIiVk/qkzmMeUh0+17uSJdj8nGKnBSl27c9NuNP/Ly2iczW4NCqRv22ecOZN5jLZung/fjf7BTl/XrQ51BGGHgrOQrcDbdpr1hhAxpqrckUGiIKIRTcSisQTCV96a5o4AgFAohI0bvoMH3gxMNALiTBeprhb7sHmd0na8eQQ9/THhuVpxinLJLVCe+8e09Wn+dZF4uRX4cvK6lDFBeImitJA7OztRW1uLA7rglwzGfsY8REEfK1OytWOn1DVi+jV3pBrVV50VwcaNG9Hc3IzFF50FAAhfeWtG0DEUCqG9vV34XiJfaFzQ1wM4VZwS/OjVaevTGulPu3BhUeTeUuN5wisUjYWs+Ul/u+M5HNm2ITXNWT9U1CyNrOPmudjZdwSdr/SbprLxhI5lcTXOrhI2wenYth/RukZMu3AhxlQVkXAQky+YsDq79w2krV0LOpaHq7Bxw3eYbT312M2d1h4ay0/2T2b518dPptwR+S9WIYqDohBkffDu6IuPpsRYw5ixwKNj2360LpqFBTXTTYtKREEfmS21cVAqgFQTIn0urF749aKoAGhuTu/SxhKF1kWzUk3pjSgKUFHmz3BpaJkoTfUR4fdAga8J3CqCIUqPgnJZ8IZu6v2kMhkLPPQ30o62hTiwdjEeWjbP8aBPV08UrY/vSRuUakRzi/BETwXSvgNe3wlgoryaRbDMhxvnR9ICZA8tm4fdKz8xcZ7AHy7TmEgWq8NUnT4/WyiXmXCKgrGQRVaI3orkBb94mQxGjFVfbgR97n3ydeEkaY2DsQTWL5vHTanTfwciUdjRthD3de1F56v9abPh4slxbN4VzfADm/UCqQwFsHJJnSPWX7bWpResU8plJpyiYARZJDian3SwtxvjI0MZ5+or7mRg9cZw6ubu6okKm+TrmRkOmvbl0L4DM1FY3TQH3fsGMq7BKjvu2LYfA7t/ySySiTjsH822FNqtEm0rUC4z4RQF47IQCU7rolkY2f9fE6ljQ39M+7kveBqmX3NHhv84GPCjkjPy3c6NJLttlt3GKkDKJaBlAfCIxhLc8fX6zyJryf12x3NppdL6Um6nW3Bma116wTqlXGbCKQpGkEU9aJvqIxh99d8zKtkAQAlUZIixX1GwZukcrFxSl5HqFvBZ941a6RssKxTNl1dnWHiivhWDQ6MI+NN/bhQF2T6+f3zpx9xSacBZ/2i2vYW90JvYrSIYovQoGEE2s0IOv3uQeR4rmJeWrmXUOBu9eqwEdWSEojIUwIKazECcaDhrclzF5EllQlGQbbQ0HBOXSgPOWaDZWpdesU4pl5lwgoLxIZsF16qrq5kN5lnBPE0UO7btzxgQmhxTM/yPZjmmom2z8dzG2VUZLTyB9OnSR+NJtD6xJ+1zAxMiK0rFO5ZIpjIkWPC+QyC90ZJMYNQpCzTboClV2hHFREFNDBHR2dmJlpYWxOPx1GtKWXmG/zjgU9Bx88ToIt4EEP0kZpnpELz0sMpQAEPJ8QzxNQ4qjcVHcGIkM6OhMhRAz7dOCaxZ9oPdIZvG9Rs7ywHp3yVNxyAIa8hODCkYl4UZzc3N2LhxI2pqaqAoCmpqavC1VQ/gnEtOCVo4GEiJMWDuf+zqieKux/aYuiN422ZVZc/cOzEyhuHRcaw/OfOOJcZA5pBUzVepn42nfz+ZgaWswKPRwteXcivKRCn3eUu/jtNOZlmQGBOEOxSNhWyH+7r2ZpRJa9YfAKE1qreiAbZbY/mm3cISbM2irRVMjj6wdjHzdauluiJLn5dSR2PtCcIZinKmnpN09USxeVc0TTAVADfOP9W1TdTf2Ghds3KVzcqvNcs0HAwglshscM+yhEXvJ0IUeHSrdzFBENYoGpeFVVgCpeJUMx9RFoGsWJl1ftNEfdX17PQ7bbKJE4gCj5S2RRDeoGQtZLOCAl71lZbDLCNW2jGrnu7NsID1op6LTAGzajIa2kkQ+adkBdlMoHjbeKuWoyZ0Zj5ftwWxddEstD6+J62Hhp0iGIIg3KNkBdnMb+q01eoJC9SBIhiCINyjZAVZRnA9IaIOIVsEQxBE/ihZQQaKS3BFdPVEudke1CKSILxDyWZZlApa/jEPahFJEN6BBLnI4U2dBijXmCC8BglykSNySVCuMUF4CxLkIofnkojoppEQBOENSJCLHK/0CyYIwpySzrIoBahfMEEUDiTIJUCppPcRRKFDLguCIAiPQIJMEAThEUiQCYIgPAIJMkEQhEcgQSYIgvAIlmbqKYoyAKDPveUQBEEUJTWqqlaZHWRJkAmCIAj3IJcFQRCERyBBJgiC8AgkyARBEB6BBJkgCMIjkCATBEF4BBJkgiAIj0CCTBAE4RFIkAmCIDwCCTJBEIRH+P+W7vfXVBUaCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# # Testing the linear fit functionality of the function above\n",
    "\n",
    "# from sklearn import datasets, linear_model\n",
    "# import matplotlib.pyplot as plt\n",
    "# import numpy as np\n",
    "# from sklearn.metrics import mean_squared_error, r2_score\n",
    "    \n",
    "# # Load the diabetes dataset\n",
    "# diabetes = datasets.load_diabetes()\n",
    "\n",
    "\n",
    "# # Use only one feature\n",
    "# diabetes_X = diabetes.data[:, np.newaxis, 2]\n",
    "\n",
    "# # Split the data into training/testing sets\n",
    "# diabetes_X_train = diabetes_X[:-20]\n",
    "# diabetes_X_test = diabetes_X[-20:]\n",
    "\n",
    "# # Split the targets into training/testing sets\n",
    "# diabetes_y_train = diabetes.target[:-20]\n",
    "# diabetes_y_test = diabetes.target[-20:]\n",
    "\n",
    "# # Create linear regression object\n",
    "# regr = linear_model.LinearRegression()\n",
    "\n",
    "# plt.figure()\n",
    "# plt.scatter(diabetes_X_train,diabetes_y_train)\n",
    "\n",
    "# # Train the model using the training sets\n",
    "# regr.fit(diabetes_X_train, diabetes_y_train)\n",
    "\n",
    "# # Make predictions using the testing set\n",
    "# diabetes_y_pred = regr.predict(diabetes_X_test)\n",
    "\n",
    "# # The coefficients\n",
    "# print('Coefficients: \\n', regr.coef_)\n",
    "# # The mean squared error\n",
    "# print(\"Mean squared error: %.2f\"\n",
    "#       % mean_squared_error(diabetes_y_test, diabetes_y_pred))\n",
    "# # Explained variance score: 1 is perfect prediction\n",
    "# print('Variance score: %.2f' % r2_score(diabetes_y_test, diabetes_y_pred))\n",
    "\n",
    "# # Plot outputs\n",
    "# plt.scatter(diabetes_X_test, diabetes_y_test,  color='black')\n",
    "# plt.plot(diabetes_X_test, diabetes_y_pred, color='blue', linewidth=3)\n",
    "# plt.plot(diabetes_X_test, diabetes_X_test*regr.coef_ + regr.intercept_, '-.', color='red', linewidth=3)\n",
    "\n",
    "# plt.xticks(())\n",
    "# plt.yticks(())\n",
    "\n",
    "# plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-15T14:59:18.482723Z",
     "start_time": "2019-11-15T14:59:18.474388Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluation: \n",
      "\n",
      "   MAE  = 41.22709128976146\n",
      "   MSE  = 2548.07239872597\n",
      "   RMSE = 50.478434986892864\n",
      "   RAE  = 0.7215101730794795\n",
      "   RSE  = 0.5274245520177285\n",
      "   RSQ  = 0.47257544798227147\n",
      "   PE  = -6.630079274986685e+31\n",
      "   linear slope  = 0.4598381320767749\n",
      "   linear intercept  = 85.52527409006767\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(41.22709128976146,\n",
       " 2548.07239872597,\n",
       " 50.478434986892864,\n",
       " 0.7215101730794795,\n",
       " 0.5274245520177285,\n",
       " 0.47257544798227147,\n",
       " -6.630079274986685e+31)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# model_eval(diabetes_y_test,diabetes_y_pred)"
   ]
  },
  {
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}