DATA602 Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load all the necessary libraries \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from pandas import Series, DataFrame\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "# Load the file\n",
    "airport=pd.read_csv('D:/Machine Learning/airport.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(19, 13)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check the number of rows and columns \n",
    "airport.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    year     dep     arr  dep_dem  arr_dem  eff_aar     adr  cap_aar     cap  \\\n",
      "0   2000  172170  158570   215157   175647   339821  368792   339821  708613   \n",
      "1   2001  147010  133806   176997   144244   320077  336135   320077  656212   \n",
      "2   2002  142779  129089   169072   138061   317088  343905   317088  660993   \n",
      "3   2003  140023  126800   168626   138507   303844  337738   303844  641582   \n",
      "4   2004  161351  160474   203993   179710   302256  337324   302256  639580   \n",
      "5   2005  176738  176139   249869   209219   322701  321895   322836  644731   \n",
      "6   2006  191852  191288   326131   237164   335796  331485   337419  668904   \n",
      "7   2007  222946  222505   407468   332520   365119  343131   370455  713586   \n",
      "8   2008  221203  221458   352677   325057   365838  359354   372714  732068   \n",
      "9   2009  209595  209710   312026   288945   373893  350917   380245  731162   \n",
      "10  2010  200902  200595   277480   253438   349015  333458   351570  685028   \n",
      "11  2011  205231  204685   266431   263065   375289  339264   379989  719253   \n",
      "12  2012  202430  202162   254435   233863   409216  333894   412665  746559   \n",
      "13  2013  203798  203443   264611   252442   410030  327404   416274  743678   \n",
      "14  2014  213113  212793   269706   259825   416622  326557   422429  748986   \n",
      "15  2015  220122  220306   285469   275114   409928  323723   417113  740836   \n",
      "16  2016  225989  226432   296871   271421   421207  348720   426426  775146   \n",
      "17  2017  223126  223879   317699   299227   412406  332844   418746  751590   \n",
      "18  2018  227425  227342   312875   270243   430599  327024   436230  763254   \n",
      "\n",
      "    caputil  eff_dep  eff_arr   saer  \n",
      "0     46.67    91.93    94.31  93.00  \n",
      "1     42.79    93.41    95.48  94.34  \n",
      "2     41.13    94.81    95.87  95.29  \n",
      "3     41.59    94.26    94.73  94.47  \n",
      "4     50.32    93.22    94.93  94.02  \n",
      "5     54.73    91.55    93.92  92.63  \n",
      "6     57.28    91.98    94.20  92.92  \n",
      "7     62.42    90.42    91.82  91.05  \n",
      "8     60.47    89.35    92.10  90.67  \n",
      "9     57.35    89.91    92.47  91.14  \n",
      "10    58.61    92.05    92.89  92.45  \n",
      "11    56.99    91.84    92.76  92.30  \n",
      "12    54.19    93.28    94.91  94.06  \n",
      "13    54.76    93.27    94.07  93.66  \n",
      "14    56.86    93.58    93.62  93.60  \n",
      "15    59.45    93.68    93.92  93.80  \n",
      "16    58.37    93.35    94.09  93.71  \n",
      "17    59.47    92.18    93.47  92.81  \n",
      "18    59.58    92.29    94.00  93.08  \n"
     ]
    }
   ],
   "source": [
    "# View the dataset\n",
    "print(airport)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a data frame\n",
    "airport=pd.DataFrame(airport, columns=['caputil', 'dep'])\n",
    "# Define X and Y\n",
    "X=pd.DataFrame(airport['caputil'])\n",
    "Y=pd.DataFrame(airport['dep'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the model\n",
    "reg = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine three folds\n",
    "scores=[]\n",
    "kfold = KFold(n_splits=3, shuffle=True, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.9029142706220207, 0.9048352664783202, 0.5119143545760048]\n"
     ]
    }
   ],
   "source": [
    "# Create the iteration for the three folds\n",
    "for i, (train, test) in enumerate(kfold.split(X, Y)):\n",
    "    reg.fit(X.iloc[train,:], Y.iloc[train,:])\n",
    "    score = reg.score(X.iloc[test,:], Y.iloc[test,:])\n",
    "    scores.append(score)\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of Determination: 0.8897612752328166\n"
     ]
    }
   ],
   "source": [
    "# Determine the coefficient of determination with sklearn\n",
    "reg = LinearRegression().fit(X, Y)\n",
    "print('Coefficient of Determination:', reg.score(X, Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients: [[4202.71508638]]\n"
     ]
    }
   ],
   "source": [
    "# Determine the coefficient for X\n",
    "print('Coefficients:', reg.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept: [-33354.0929305]\n"
     ]
    }
   ],
   "source": [
    "# Determine the intercept\n",
    "print('Intercept:', reg.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict the Y values\n",
    "Y_pred=reg.predict(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted Values:  [[7.23548102e+08]\n",
      " [6.17807791e+08]\n",
      " [6.00026103e+08]\n",
      " [5.88443420e+08]\n",
      " [6.78078928e+08]\n",
      " [7.42746105e+08]\n",
      " [8.06265941e+08]\n",
      " [9.36945164e+08]\n",
      " [9.29619831e+08]\n",
      " [8.80834714e+08]\n",
      " [8.44300512e+08]\n",
      " [8.62494066e+08]\n",
      " [8.50722261e+08]\n",
      " [8.56471575e+08]\n",
      " [8.95619866e+08]\n",
      " [9.25076696e+08]\n",
      " [9.49734026e+08]\n",
      " [9.37701652e+08]\n",
      " [9.55769124e+08]]\n"
     ]
    }
   ],
   "source": [
    "# Determine the values\n",
    "print('Predicted Values: ',Y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize predicted errors\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "predicted = cross_val_predict(reg, X, Y, cv=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted Errors:  [[160835.55508042]\n",
      " [143928.21704009]\n",
      " [138393.80845549]\n",
      " [140358.20860634]\n",
      " [180662.82450213]\n",
      " [198806.56295092]\n",
      " [209181.70395083]\n",
      " [231558.86084681]\n",
      " [220561.15973044]\n",
      " [207488.73277439]\n",
      " [214087.71393934]\n",
      " [207154.97715331]\n",
      " [193506.03241906]\n",
      " [195900.52219425]\n",
      " [204795.19251424]\n",
      " [215536.69965656]\n",
      " [210137.99688364]\n",
      " [214613.01298902]\n",
      " [216056.74708452]]\n"
     ]
    }
   ],
   "source": [
    "# Determine the predicted errors\n",
    "print('Predicted Errors: ', predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXgUVdbA4d8hYRUQEFA0EEAjISgiBEbHmRFFcUER9yVK3D/DMgjKJgghoiKIIKJxGRi2oKLi6DgygqDiqCgRNcjiBMFIBBFkESQhLOf7o26YTrqzAN3pLOd9nn5SdW8tpytJn657q26JqmKMMcYEU7VwB2CMMabyseRijDEm6Cy5GGOMCTpLLsYYY4LOkosxxpigs+RijDEm6Cy5mHJLRFqKiIpIpJtfICKJR7GdFiKyR0QiQhCjishpwd5uZSEiP4jIRW76IRH5W2mWPYr9/FlEvjvaOE3wWXIxx8R9IOS4D+8tIvJ3Eakbin2p6mWqOrOUMR3+kFLVH1W1rqoeDEVcwSIizURkmohsFpHdIrJWRMaIyHFhime4iCwNUN5YRPJE5Iwj2Z6qPqaqdwcptgJJXVU/VtU2wdi2CQ5LLiYYrlTVukBHoDMwsvAC4rG/tyKISCPgM6A2cK6q1gMuBhoApwZYPrIMwpoN/FFEWhUqvwlYqarflkEMpoKyf3YTNKr6E7AAOANARD4UkUdF5BNgL9BaRI73+Xb+k4iMzW+uEpEIEXlSRLaJyHqgh+/23fbu9pm/R0TWuG/5q0Wko4jMBloA/3RnU0MCNK+dLCJvi8h2EVknIvf4bDNZROaJyCy33VUiEl/CW79cRNa7uCeISDURqem2f6bPtpu6s7wmAbYxCNgN3KqqP7jjuVFVB6hqhltfRaSviGQCma7sjyKyXER2uZ9/9Nnf7S6u3SKyQUQSXPlpIvKRW2ebiLxaxO8zG1gC3Faoqjcw023rVBFZIiK/um2liUiDQNtzx3aOz/xtIpLl1h1RaNkuIvKZiOx0fytTRaSGq8s/m/rG/Y5vFJGuIpLts35b9/ey0/0Oe/rUzRCRZ0XkX+7YfC4ifgncHCNVtZe9jvoF/ABc5KabA6uAR9z8h8CPQDsgEqgO/AN4ATgOaAp8AfyfW/4+YK3bTiPgA0CBSJ/t3e2mrwd+wjtTEuA0ILpwTG6+ZaHtfAQ8B9QCOgBbgW6uLhnIBS4HIoDHgWXFvH91cTbCS2r/9YnxOeAJn2UHAP8sYjvLgDElHGsFFrl91XY/d+B9+EcCN7v5E9zx/Q1o49ZtBrRz0y8DI/C+XNYC/lTMPhOATJ/5NkAe0MTNn4Z3hlUTaAIsBSYX8feRDMxx03HAHuAvbt2ngAM+y3YCznHvqyWwBri/0LE4zWe+K5DtpqsD64CHgBrAhXiJO/9YzAC2A13c9tOAV8L9v1TZXmEPwF4V++U+PPYAO4Es94Fa29V9CKT4LHsisC+/3pXdDHzgppcA9/nUdafo5PIeMKCYmAImF7zEdRCo51P/ODDDTScD7/vUxQE5xbx/BS71me8DLHbTfwA2AtXcfDpwQxHbyfR978Xs60Kf+duALwot8xlwO15y2Qlc63u83TKzgBeBqFL8fuvgJak/uvlHgbeKWb4X8FWg3wUFk8so3w90F2+e7++t0HbvB94sdCyKSi5/Bn7OP+6u7GUg2U3PAP7mU3c5sDbc/0uV7WXNYiYYeqlqA1WNVtU+qprjU7fRZzoa71vlZtdcsRPvLKapqz+50PJZxeyzOfD9UcR6MrBdVXcX2s8pPvM/+0zvBWqV0MdROOaTAVT1c+B34HwRicX7lv92Edv4Fe/soiS++zoZ/2OUBZyiqr8DN+KdDW52TUCxbpkheGd7X7gmozvh8NVce9zrefce9gKvAb1FRPDOZA5fVOGa+l5xTZy/AXOAxqV4HwV+1y7eX322e7qIvCMiP7vtPlbK7R7etqoeKnxcfOYL/45DchFKVWbJxYSa77DbG/HOXBq7ZNRAVeurajtXvxkvaeRrUcx2NxKgozvAPgvbBDQSkXqF9vNTMeuUpHDMm3zmZwK34p1lvK6quUVs433g6lJc9OD73jbhJWxfh9+Lqr6nqhfjJa21wEuu/GdVvUdVTwb+D3hORE5T72quuu51X6H3cANe81c94B2fusddTO1Vtb57r1LCe4BCv2sRqYPXnJcv1cUc47b7UCm3C95xaV7oWB7r79gcIUsupsyo6mZgITBRROq7ju9TReR8t8g84K8iEiUiDYFhxWzub8CDItJJPKeJSP4H7RagdRExbAQ+BR4XkVoi0h64C6/d/WgNFpGGItIcr1/Ft4N8NnA13ofurGK28RRQH5iZ/z5E5BQRecrFGMi7wOkicouIRIrIjXjNeO+IyIki0lO8y5j34TVdHnTbvV5Eotw2duAlh+Iu0/4Yr4ntRbymrDyfunpu2ztF5BRgcDHb8fU6cIWI/Ml11KdQ8POoHl5z3B53xpVUaP0if8dA/hnjEBGpLiJdgSuBV0oZmwkCSy6mrPXG62RdjffB9jr/aw56Ca8v5RtgBTC/qI2o6mt47f9z8Tpr/4HXwQ3et+mRruntwQCr34zXD7MJeBMYraqLjuE9vQV8CXwN/AuY5hNntnsvivchXdT72Q78EdgPfC4iu4HFwC68zulA6/wKXAE8gNekNAS4QlW34f1vP+De43bgfLz+IPAugvhcRPbgNdMNUNUNxcSmeIkxGv8EOQbvEvRd7r0X+TsrtM1VQF+8399mvL+FbJ9FHgRuwfvdvkTBhA1e/81M9zu+odC284CewGXANrx+wN6qurY0sZngENehZYwJERGZDmxSVb/7f4yprMriRixjqiwRaQlcA5wd3kiMKVvWLGZMiIjII8C3wITimp2MqYysWcwYY0zQ2ZmLMcaYoLM+F6dx48basmXLcIdhjDEVypdffrlNVf3Gy7Pk4rRs2ZL09PRwh2GMMRWKiAQcScOaxYwxxgSdJRdjjDFBZ8nFGGNM0FlyMcYYE3SWXIwxxgSdJRdjjKmg0tLm0iomlmoREbSKiSUtbW64QzrMLkU2xpgKKC1tLkkDB1OnWz+a94ojJ3s1SQO9Jx4kJNwS5uhs+JfD4uPj1e5zMcZUFK1iYsmJT6RW9P8e95OblUHt9JlsyCy7pwuIyJeqGl+43JrFjDGmAspan0nNqLgCZTWj4shanxmmiAqy5GKMMRVQdOsY9mWvLlC2L3s10a1jwhRRQZZcjDGmAhqbPIq9i6eSm5WBHjxAblYGexdPZWzyqHCHBliHvjHGVEj5nfYjk1PImpdJdOsYJk6aUC4688E69A+zDn1jjDly1qFvjDGmzFhyMcYYE3SWXIwxxgSdJRdjjKlC9u3bx4oVK0K+H0suxhhTRSxatIgzzzyTiy66iK1bt4Z0X5ZcjDGmktu0aRM33XQT3bt3JzMzkx07djBs2LCQ7tOSizHGVFIHDhxg8uTJxMbG8uqrrxaomz59Op988knI9m03URpjTCX06aefkpSUREZGhl+diHDvvffStm3bkO3fkosxxlQi27ZtY9iwYUybNi1gfceOHUlNTaVLly4hjcOaxYwxphI4dOgQf/vb32jTpk3AxFK/fn2mTp3KF198EfLEAnbmYowxFd7XX39NUlISy5YtC1h/6623MmHCBE466aQyi8mSizHGVGB79uyha9eu7Nq1y6+ubdu2PPfcc3Tt2rXM47JmMWOMqcDq1q3LyJEjC5TVqVOHcePG8fXXX4clsYAlF2OMqfAGDBhAu3btAOjVqxerV69m6NCh1KhRI2wxWbOYMcZUAHv37uXHH38kNjbWr6569eq89NJL/Prrr1xxxRVhiM6fJRdjjCnn3nnnHfr3709kZCQrV66kVq1afsuce+65YYisaNYsZowx5VRWVha9evXiyiuv5IcffmDdunU0OyWKtLS54Q6tRCFLLiLSXEQ+EJE1IrJKRAa48gkislZEMkTkTRFp4LPOcBFZJyLficglPuWXurJ1IjLMp7yViHwuIpki8qqI1HDlNd38OlffMlTv0xhjgi0vL49x48bRtm1b3nrrrQJ1O3fu4v/+OrDcJ5hQnrkcAB5Q1bbAOUBfEYkDFgFnqGp74L/AcABXdxPQDrgUeE5EIkQkAngWuAyIA252ywI8AUxS1RhgB3CXK78L2KGqpwGT3HLGGFPuffDBB5x11lkMHz6cnJycgpVSjXode1DngvsYmZwSngBLKWTJRVU3q+oKN70bWAOcoqoLVfWAW2wZEOWmrwJeUdV9qroBWAd0ca91qrpeVfOAV4CrRESAC4HX3fozgV4+25rppl8HurnljTGmXPr555+59dZbufDCC1m7dq1ffc2TY2mWOJlG3e6hdquOZK3PDEOUpVcmHfquWeps4PNCVXcC+UN1noKXbPJluzKAjYXK/wCcAOz0SVS+y5+Sv46qHhCRXW75bYXiuhe4F6BFixZH/saMMeYYHTx4kNTUVEaMGMFvv/3mV1+tWjXqdupFgwtuR8Q7H9iXvZro1jFlHeoRCXmHvojUBd4A7lfV33zKR+A1naXlFwVYXY+ivLhtFSxQfVFV41U1vkmTJkW/CWOMCYH8cb769+8fMLHcfffdPPvsc+gPy9j347fowQPkZmWwd/FUxiaPCkPEpRfSMxcRqY6XWNJUdb5PeSJwBdBNVfM/9LOB5j6rRwGb3HSg8m1AAxGJdGcvvsvnbytbRCKB44HtwXxvxhhzLLKzsznvvPM4cOCAX91ZZ51Famrq4cuL69Wrx8jkFLLmZRLdOoaJkyaQkHBLWYd8REJ5tZgA04A1qvqUT/mlwFCgp6ru9VnlbeAmd6VXKyAG+AJYDsS4K8Nq4HX6v+2S0gfAdW79ROAtn20luunrgCU+ScwYY8IuKiqK++67r0BZvXr1mDx5Munp6QXuW0lIuIUNmWs5dPAgGzLXlvvEAiCh+swVkT8BHwMrgUOu+CFgClAT+NWVLVPV+9w6I/D6YQ7gNaMtcOWXA5OBCGC6qj7qylvjdfA3Ar4CblXVfSJSC5iN18+zHbhJVdcXF298fLymp6cH460bY0yp7Ny5k9jYWLZs2cJNN93ExIkTOfnkk8Md1hERkS9VNd6v3L7Qeyy5GGNCYffu3ezYsaPIi4b++c9/Urt2bS666KIyjgzS0uZ6zW3rvea2scmjjvisqKjkYnfoG2NMCKgqr732GrGxsSQkJFDUF/krr7zSL7Gkpc2lVUws1SIiaBUTG5IbJtPS5pI0cDA58Yk0HzSfnPhEkgYODtq+7MzFsTMXY0ywZGZm0q9fPxYuXHi4bMaMGSQmJhazlif/Q79Ot37UjIpjX/Zq9i6eSmqQO/FbxcSSE59Irej2h8tyszKonT6TDZn+99kUxc5cjDEmxHJychg9ejRnnHFGgcQCMHjwYHbs2FHiNkYmp1CnWz9qRbdHIiKpFd2eOt36FXlH/tGe5WStz6RmVFyBsppRcUG7OdOSizGmXCiLpqBQWrBgAWeccQYpKSnk5eUVqIuMjOTOO+8s1fNVjuRD/1iatqJbx7Ave3WBsmDenGnJxRgTdqFu/w+ljRs3cu2113L55Zezfr3/Raldu3blm2++Ydy4cRx33HElbu9IPvSP9CzH19jkUexdPJXcrIyQ3JxpycUYE3bH8iEZLvv372fChAm0bduW+fPn+9WfeOKJzJkzhyVLlhAXFxdgC4EdyYf+sTRtJSTcQuqkCdROn8nGp66hdvrMoPbr2MPCjDFhl7U+k+a9AnxIziufgzMuXbqUpKQkVq9e7VdXrVo1+vTpwyOPPEKDBg0CrF28/A/30tyRH906hpzs1QU65Y+kaSsh4ZaQ3ZBpycUYE3bH+iFZllasWMH5558fsK5Lly4899xzdOrU6Zj2UdoP/bHJo0gaOBgKXVk2cdKEY9p/MFizmDEm7ELd/h9MZ599tt9z6hs2bMjzzz/PZ599dsyJ5UiEumnrWNh9Lo7d52JMeAXjbvGysmHDBuLi4sjNzeX222/niSeeoGnTpuEOKyxs+JcSWHIxxvjasWMHBw8epHHjxgHr//73v3Paaafx5z//uYwjK1/sJkpjjCkFVWX27NnExsYycODAIpe74447qnxiKY4lF2OMcVatWkXXrl3p3bs3v/zyC3PmzOHDDz8Md1hBUdY3qVpyMcZUeXv27GHo0KF06NCBpUuXFqjr06eP3x33FU04blK15GKMOayiD8FypFSVN998k7i4OMaPH+/3VMiaNWty0003FTmicUURjptU7T4XYwxQcDTe5r3iyMle7d1DAeX2qq1jsX79evr378+7774bsP7SSy9l6tSpnHrqqWUcWfCF4yZVO3MxxgAVcwiWo5Gbm8sjjzxCu3btAiaWqKgo3njjDd59991KkVgg9INUBmLJxRgDhH4I9vJg4cKFtG/fnlGjRpGbm+tXX79BQ5KTx3DNNdcgImGIMDTCcZOqJRdjDBCeb7fBUpq+on/9619ccsklZGb6J8vqTVpyUuLT1L50MAOHP1zp+prCcSe/3UTp2E2UpqorqycgBltp4z5w4ACdO3fm66+/PlxWLSKCuvFX0+D8xMNnKkfzNMaqzG6iNMYUqzyPU1Wc0vYVRUZGkpqaCoCIkJSUxKFDh2jw51sLNIFVtqbAcLEzF8fOXIypmKpFRNB80Hwkwrv49eDeXVAtkp+euZlDBw/6LT9+/HguuOACOnfuHLTnyFdlduZijKmU8vuKVA+x++t/s+ml+9ix6Pki+4qGDBlC586dgYo1GnOFo6r2UqVTp05qjKl45sxJ0zoNmmj1E5orcPg1fvyEUq/f8rQ2KtWqacvT2uicOWkhWaeyAtI1wGeq3URpjKmwdu3axeefLyNn1za/u+gXLnyPBx98oMRLio/0aYxV7WbTo2XNYsaYCkdVmTt3LrGxsTzzzDN+iaVOnTpcfPHFHDp0KOj7rio3mx4rSy7GmApl7dq1XHTRRSQkJPDzzz/71V999dWsWbOGIUOGEBEREfT9V4WbTYPBkosxpkLYu3cvI0aMoH379ixZssSvvlWrVrzzzjvMnz+fFi1ahCyOinyzaVmy5GKMKffefvtt4uLieOyxx9i/f3+Buho1avDwww+zatUqevToEfJY7Aqz0rEOfWNMuTZt2jTuvvvugHUXXXQRzz77LKeffnqZxZPfaT8yOYWseZlEt45hYgW42bSs2U2Ujt1EaUz5tHv3bmJjY9m0adPhsmbNmjF58mSuv/76SjXAZEVkN1EaYyqkevXqMXnyZAAiIiIYOHAga9eu5YYbbrDEUo5Zs5gxplz45ZdfaNy4MdWq+X/nve666xg8eDAJCQmcddZZYYjOHKmQnbmISHMR+UBE1ojIKhEZ4MobicgiEcl0Pxu6chGRKSKyTkQyRKSjz7YS3fKZIpLoU95JRFa6daaI+xpT1D6MMeXPgQMHmDJlCjExMUyfPj3gMiLC+PHjLbFUIKFsFjsAPKCqbYFzgL4iEgcMAxaragyw2M0DXAbEuNe9QCp4iQIYDfwB6AKM9kkWqW7Z/PUudeVF7cMYU44sW7aMzp07M2DAAH777TeGDh3Ktm3bwh2WCYKQJRdV3ayqK9z0bmANcApwFTDTLTYT6OWmrwJmueFqlgENRKQZcAmwSFW3q+oOYBFwqaurr6qfufFtZhXaVqB9GGPKgV9//ZV7772Xc889t8DzVbZv386wYfZdsDIokw59EWkJnA18DpyoqpvBS0BAU7fYKcBGn9WyXVlx5dkByilmH4XjuldE0kUkfevWrUf79kwVVZqnH5qCDh06xPTp02nTpg0vvfSSX339+vWt6auSCHmHvojUBd4A7lfV34q5uiNQhR5Feamp6ovAi+Bdinwk65qqzQYvPHIZGRkkJSXx6aefBqy/5ZZbePLJJ2nWrFkZR2ZCIaRnLiJSHS+xpKnqfFe8xTVp4X7+4sqzgeY+q0cBm0oojwpQXtw+jAkKG7yw9Hbv3s2gQYPo2LFjwMTSpk0bFi9eTFpamiWWSiSUV4sJMA1Yo6pP+VS9DeRf8ZUIvOVT3ttdNXYOsMs1ab0HdBeRhq4jvzvwnqvbLSLnuH31LrStQPswJihs8MKSqSrz5s0jNjaWSZMmcbDQUyFr167NY489xjfffMOFF14YpihNqISyWew84DZgpYjk99g9BIwD5onIXcCPwPWu7l3gcmAdsBe4A0BVt4vII8Byt1yKqm5300nADKA2sMC9KGYfxgRFdOsYcrJXF3g8rg1eWNCjjz7Kww8/HLCuZ8+ePP3007Rs2bJsgzJlxoZ/cWz4F3MkfPtcakbFsS97NXsXTyXVxpg6bOPGjbRt25bff//9cFl0dDRTpkyhZ8+eYYzMBJMN/2JMECUk3ELqpAnUTp/JxqeuoXb6TEsshTRv3pzRo0cDUL16dR566CFWr159zInFrtKrGOzMxbEzF2OOzrZt22jcuHHAuv3799OvXz8GDhxIbGzsMe/LzhjLn6LOXCy5OJZcjDkyeXl5TJo0iZSUFN588026d+8e8n22ioklJz6xQF9XblYGtdNnsiFzbcj3b/xZs5gxJmg+/PBDOnTowLBhw9i7dy99+/YlNzc35Pu1q/QqDksuxphS27JlC7fddhsXXHABa9asOVy+bt06JkyYEPL92yOGKw5LLsaYEh08eJBnn32WNm3aMGfOHL/6hg0b0rx58wBrHpvCnfeXdb/IHjFcQRR7n4uIDCquvtDNkcaYSmj58uUkJSXx5ZdfBqy/8847GTduHE2aNAnqfgMNsTPntancev3VLFg40x4xXM6VdBNlPfezDdAZ7853gCuBpaEKyhgTfjt27GDEiBE8//zzBLrw58wzzyQ1NZXzzjsvJPv3HWIH8H5268eChdZ5XxEUm1xUdQyAiCwEOrqh8xGRZOC1kEdnjClzqsrs2bN58MEHCTRaeN26dUlJSaF///5ERoZukI+s9Zk07xWg836edd5XBKXtc2kB5PnM5wEtgx6NMSbs/vrXv5KYmBgwsdxwww2sXbuWgQMHhjSxgHXeV3SlTS6zgS9EJFlERuM9l2VW6MIyxoTL7bffTuFHY8TExLBw4UJeffVVTjnllCLWDK6xyaOs874CK1VyUdVH8QaS3AHsBO5Q1cdCGZgxJjw6depEnz59AKhVqxaPPPIIK1eu5OKLLy7TOGyInYqt1Hfoi8ifgBhV/buINAHqquqGkEZXhuwOfVPV7Nq1i+OPPz5g3c6dO+nfvz9jxoyhdevWZRyZqUiO6Q591xQ2FBjuiqoD/he7G2PKvdzcXMaMGUNUVBSrVq0KuEyDBg2YPXu2JRZz1Erb53I10BP4HUBVN/G/y5SNMRXEv//9b8444wySk5PZs2cPffr0CXiZsTHHqrTJJU+9v0AFEJHjQheSMSbYsrOzuf7667nsssv4/vvvD5cvXbo04B33xhyr0iaXeSLyAtBARO4B3gf+FrqwjDHBsH//fiZOnEhsbCyvv/66X33Tpk2pXbt2GCIzlV2pLlRX1SdF5GLgN7y79Uep6qKQRmaMOSb/+c9/6NOnDytXrvSrExGSkpIYO3YsDRs2DEN0prIrVXIRkSdUdSiwKECZMaYc2bp1K0OHDuXvf/97wPr4+HhSU1OJj/e7wMeYoClts1igC9wvC2Ygxphjc+jQIV588UXatGkTMLE0aNCA1NRUli1bZonFhFxJoyInAX2AU0Ukw6eqHvBpKAMzxhyZG2+8MWC/CkBiYiLjx4+nadOmZRyVqapKOnOZizcC8lvuZ/6rk6omhDg2Y8wRSEjw/5ds164dH330ETNmzLDEYspUsclFVXep6g/A08B2Vc1S1Sxgv4j8oSwCNMaUzlVXXUWPHj0AqFOnDuPHj+err77iL3/5S5gjM1VRaYc1TQU6+sz/HqDMGFMG9uzZQ926df3KRYQpU6ZQr149xo8fH5InQxpTWqXt0Bf1uY1XVQ9R+sRkjAmC33//neHDh9O6dWu2bNkScJnWrVvz8ssvW2IxYVfa5LJeRP4qItXdawCwPpSBGWM8qspbb71FXFwc48aNY+vWrQwePDjcYRlTrNIml/uAPwI/AdnAH4B7QxWUMcazYcMGevbsSa9evfjxxx8Pl8+ePZuPPvoojJEZU7zS3qH/C3BTiGMxxjj79u3jySefZOzYseTm5vrVn3zyyeTk5IQhMmNKp6T7XIao6ngReQY3aKUvVf1ryCIzpop6//336du3L//973/96iIiIhgwYADJycnUq2cDk5vyq6QzlzXupz1Fy5gQ27RpEw888ACvvPJKwPrzzjuP1NRUzjzzzDKOzJgjV2xyUdV/up8zyyYcY6qeAwcO8Oyzz/Lwww+ze/duv/rGjRszYcIEevfuTbVqpe0mNSa8SmoW+ycBmsPyqWrPoEdkTBWiqlx44YV8/PHHfnUiwj333MPjjz9Oo0aNwhCdMUevpK9BTwITgQ1ADvCSe+0Bvg1taKaySEubS6uYWKpFRNAqJpa0tLnhDqncEBGuueYav/Kzzz6bzz77jBdeeMESi6mQShr+5SNV/Qg4W1VvVNV/utctwJ+KW1dEpovILyLyrU9ZBxFZJiJfi0i6iHRx5SIiU0RknYhkiEhHn3USRSTTvRJ9yjuJyEq3zhQREVfeSEQWueUXiYg9rCKM0tLmkjRwMDnxiTQfNJ+c+ESSBg62BOOjX79+nHXWWQDUr1+fZ555huXLl/OHP9gIS6biKm0DbhMRaZ0/IyKtgCYlrDMDuLRQ2XhgjKp2AEa5efCG749xr3vxhpZBRBoBo/Huq+kCjPZJFqlu2fz18vc1DFisqjHAYjdvwmRkcgp1uvWjVnR7JCKSWtHtqdOtHyOTU8IdWpkLdEkxQGRkJKmpqSQkJPDdd9/Rr18/IiIiyjg6Y4KrtMllIPChiHwoIh8CHwD3F7eCqi4FthcuBuq76eOBTW76KmCWepbhPU65GXAJsEhVt6vqDryHlV3q6uqr6mduWJpZQC+fbeVfgDDTp9yEQdb6TGpGxRUoqxkVR9b6zDBFVPZ+++037r//ftq3b1/kvSnnnnsuc+bM4aSTTirj6IwJjdLeRPlvEYkBYl3RWlXddxT7ux94T0SexKHf6f8AABrVSURBVEtsf3TlpwAbfZbLdmXFlWcHKAc4UVU3u7g3i4iNMx5G0a1jyMleTa3o9ofL9mWvJrp1TBijKhuqyquvvsqgQYPYvHkzAI8//jgpKVXvrM1UPaU6cxGROsBgoJ+qfgO0EJErjmJ/ScBAVW2OdzY0LX8XAZbVoyg/IiJyr+v7Sd+6deuRrm5KYWzyKPYunkpuVgZ68AC5WRnsXTyVscmjwh1aSH333XdcfPHF3HzzzYcTC8ATTzwR8OZIYyqb0jaL/R3IA85189nA2KPYXyIw302/htePkr8932Fco/CazIorjwpQDrDFNZvhfv5SVDCq+qKqxqtqfJMmJXUhmaORkHALqZMmUDt9Jhufuoba6TNJnTSBhIRbwh1aSOzdu5eRI0dy5plnsnjxYr/6Zs2aYV9kTFVQ2uRyqqqOB/YDqGoOgc8eSrIJON9NXwjkN7y/DfR2V42dA+xyTVvvAd1FpKHryO8OvOfqdovIOe4qsd54T8vM31b+VWWJPuUmTBISbmFD5loOHTzIhsy1lTaxvPPOO7Rr145HH32U/fv3F6irXr06I0aMYPXq1Zx33nlhitCYslPaZ7LkiUhtXNOTiJwKFNvnIiIvA12BxiKSjXfV1z3A0yISCeTyv5GV3wUuB9YBe4E7AFR1u4g8Aix3y6Woav5FAkl4V6TVBha4F8A4YJ6I3AX8CFxfyvdozFHJyspiwIABvPVW4O8x3bp1Y+rUqcTGxgasN6YyEp9ngBW9kMjFwEggDlgInAfcrqofhjS6MhQfH6/p6TaEmim9vLw8nnrqKVJSUgJeBXbSSScxadIkbrzxRtxtWMZUOiLyparGFy4v8czFNTutBa4BzsFrDhugqtuCHqUxFcS+ffuIj4/n22/9B6qoVq0a/fv3Z8yYMRx//PFhiM6Y8Csxuaiqisg/VLUT8K8yiMmYcq9mzZpcdNFFfsnlnHPOITU1lQ4dOoQpMmPKh9J26C8Tkc4hjcSYCmbMmDE0a9YMgEaNGvG3v/2NTz75xBKLMZS+Q/8C4D4R+QH4Ha9pTFW1fbFrGVMJ5OXlUaNGDb/y+vXrM3nyZBYuXMi4ceNo3LhxGKIzpnwqbYd+dKByVc0KekRhYh36prDt27fz0EMPsWrVKj766CN7looxARxVh76I1ALuA04DVgLTVPVAaEI0pnw4dOgQs2bNYvDgwWzb5l23MmPGDO68884wR2ZMxVHSV7GZQDxeYrkM79kuxlRaK1eu5Pzzz+eOO+44nFgAhgwZwq+//hrGyIypWErqc4lT1TMBRGQa8EXoQzKm7O3evZsxY8YwefJkDh486Fd/wgkn8NNPP3HCCSeEITpjKp6SzlwOj2FhzWGmMlJVXn/9ddq2bcvEiRP9EkutWrUYO3YsGRkZtG9v168YU1olnbmcJSK/uWkBarv5/KvF6he9qjHl27p16+jXrx/vvfdewPoePXrwzDPP0KpVqzKOzJiKr9jkoqr2ODxT6eTk5PDEE08wbtw49u3zHyKvRYsWTJkyhZ49e9qwLcYcpdLe52JMpbB9+3a6dOnC999/71cXGRnJgw8+yMiRIznuuOPCEJ0xlYdduG+qlEaNGgXsO+natSvffPMNjz/+uCUWY4LAkoupciZPnkydOnUAaNq0KbNnz2bJkiXExcWFOTJjKg9rFjOV1sGDB4mI8O82bNGiBWPGjOGHH35g7NixNGjQIAzRGVO5WXIxlc4vv/zCkCFDyMvLY+7cuQGXefDBB8s4KmOqFmsWM5XGwYMHef7552nTpg0zZ87k5Zdf5v333w93WMZUSZZcTKXw5Zdfcu6555KUlMTOnTsPl/ft2zfg5cbGmNCy5GIqtJ07d9KvXz86d+7M8uXL/epr1KjBpk2bwhCZMVWbJRdTIakqc+bMoU2bNjz77LMUfnTEcccdx5NPPsmKFSvsDntjwsA69E2Fs3r1avr06cNHH30UsP66665j0qRJREVFlXFkxph8duZiKozff/+dYcOGcdZZZwVMLKeeeioLFizgtddes8RiTJjZmYupENavX0/Xrl3ZuHGjX13NmjUZPnw4Q4cOpVatWmGIzhhTmCUXUyFER0fTtGlTv+RyySWXMHXqVE477bQwRWaMCcSaxUyFEBERQWpq6uFRik855RRef/11FixYYInFmHLIzlxMuXPo0CGqVfP/3tO5c2f69etHzZo1GTVqFPXq1QtDdMaY0rDkYsqNTZs2MWjQIFq2bMm4ceMCLvP000/bM1aMqQAsuZiwO3DgAM888wyjR49m9+7dREZG0rt374CjFFtiMaZisD4XE1affPIJnTp1YtCgQezevRvwkk2fPn38bow0xlQcllxMWGzbto277rqLP/3pT2RkZPjV7969m19//TUMkRljgsGSiylThw4d4qWXXqJNmzZMnz7dr/74449n6tSpfPHFFzRu3DgMERpjgsH6XEyZ+eqrr0hKSuLzzz8PWH/bbbcxYcIETjzxxDKOzBgTbHbmYkJu165dDBgwgPj4+ICJpW3btnzwwQfMmjWrTBNLWtpcWsXEUi0iglYxsaSlBX6wmDHmyNmZiwmpFStW0KNHD37++We/ujp16jB69Gjuv/9+atSoUaZxpaXNJWngYOp060fzXnHkZK8maeBgABISbinTWIypjEJ25iIi00XkFxH5tlB5fxH5TkRWich4n/LhIrLO1V3iU36pK1snIsN8yluJyOcikikir4pIDVde082vc/UtQ/UeTclOP/10qlev7ld+9dVXs2bNGoYMGVLmiQVgZHIKdbr1o1Z0eyQiklrR7anTrR8jk1PKPBZjKqNQNovNAC71LRCRC4CrgPaq2g540pXHATcB7dw6z4lIhIhEAM8ClwFxwM1uWYAngEmqGgPsAO5y5XcBO1T1NGCSW86ESd26dXn66acPz7dq1Yp33nmH+fPn06JFi7DFlbU+k5pRBe+jqRkVR9b6zDBFZEzlErLkoqpLge2FipOAcaq6zy3ziyu/CnhFVfep6gZgHdDFvdap6npVzQNeAa4S7066C4HX3fozgV4+25rppl8HuondeVcmirovpVevXlx99dU8/PDDrFq1ih49epRxZP6iW8ewL3t1gbJ92auJbh0TpoiMqVzKukP/dODPrrnqIxHp7MpPAXyHu812ZUWVnwDsVNUDhcoLbMvV73LL+xGRe0UkXUTSt27desxvrqr64Ycf6NmzJy+88ELAehHhjTfeICUlhdq1a5dxdIGNTR7F3sVTyc3KQA8eIDcrg72LpzI2eVS4QzOmUijrDv1IoCFwDtAZmCcirYFAZxZK4OSnxSxPCXUFC1VfBF4EiI+Pt9vBj1BeXh5PPvkkY8eOJScnh48//phrrrmGpk2b+i1b3k4e8zvtRyankDUvk+jWMUycNME6840JkrJOLtnAfPXaT74QkUNAY1fe3Ge5KGCTmw5Uvg1oICKR7uzEd/n8bWWLSCRwPP7Nc+YYLVmyhL59+7J27drDZTt37mTIkCHMmDEjfIEdgYSEWyyZGBMiZd0s9g+8vhJE5HSgBl6ieBu4yV3p1QqIAb4AlgMx7sqwGnid/m+75PQBcJ3bbiLwlpt+283j6peoDVIVNJs3byYhIYFu3boVSCz5MjMzyc3NDUNkxpjyJJSXIr8MfAa0EZFsEbkLmA60dpcnvwIkqmcVMA9YDfwb6KuqB91ZST/gPWANMM8tCzAUGCQi6/D6VKa58mnACa58EHD48mVz9PJHLo6NjWXuXP+bDU844QSmTZvGxx9/bI8aNsYg9qXeEx8fr+np6eEOo1z6/PPPSUpK4quvvgpYf8899/D4449zwgkBr5swxlRiIvKlqsYXLrc79E2Rtm/fzvDhw3nppZcCXmbcoUMHnnvuOc4999wwRGeMKc8suZiAPvjgA2644Qa2bdvmV1evXj3Gjh1Lnz59iIy0PyFjjD/7ZDABtWnThn379vmV33zzzUycOJFmzZqFISpjTEVhoyKbgE4++WRSUv43zlabNm14//33mTt3riUWY0yJ7MylilPVIm9w7NevH6+++ipXXnklDzzwADVr1izj6IwxFZWduVRh//3vf7nkkkt49913A9ZHRkbyySef8NBDD1liMcYcETtzqYJycnJ4/PHHeeKJJ8jLy+P777/nggsuCDjuV7Vq9v3DGHPk7JOjinn33Xdp164djzzyCHl5eQCsX7+ecePGhTkyY0xlYsmlivjxxx+55ppr6NGjBxs2bPCrX758eZFD5htjzJGy5FLJ7d+/n/Hjx9O2bVvefPNNv/oTTzyRtLQ0/vWvf5W7kYuNMRWX9blUYkuXLiUpKYnVq1f71VWrVo2+ffvyyCOPcPzxx4chOmNMZWbJpRLasmULQ4YMYdasWQHru3TpQmpqKh07dizjyIwxVYUll0rmH//4B3fccQc7d+70q2vYsCHjxo3j7rvvtqvAjDEhZcmlkmnZsiW//fabX/mdd97JuHHjaNKkSRiiMsZUNfb1tZLp0KED/fv3Pzx/5pln8p///Idp06ZZYjHGlBlLLpVQSkoKp59+OhMnTuTLL7/kvPPOC3dIxpgqxpJLBbRq1Sq6detGRkZGwPr69euzatUqBg0aRPXq1cs4OmOMseRSoezZs4chQ4bQoUMHlixZQlJSEocOHQq4rD1nxRgTTpZcKgBVZf78+bRt25YJEyZw4MABAD799FNmzpwZ5uiMMcafJZdy7vvvv6dHjx5ce+21ZGdn+9UvWbIkDFEZY0zxLLmUU7m5uaSkpNCuXTsWLFjgVx8VFcUbb7xR5I2SxhgTTtYwXw4tXLiQvn37sm7dOr+6yMhIBg4cyKhRo6hbt24YojPGmJJZcilHfvrpJwYOHMhrr70WsP4vf/kLzz33HO3atSvjyIwx5shYs1g5MWPGDGJjYwMmliZNmjBr1iw+/PBDSyzGmArBkks50bhxY/bs2VOgTETo06cP3333HbfddpsNiW+MqTAsuZQTV1xxBVddddXh+fj4eJKTx/DuwsWc0LgxrWJiSUubG8YIjTGm9KzPpRx5+umnWb58OSNHjqRu3Xr0fWAodbr1o3mvOHKyV5M0cDAACQm3hDlSY4wpntijbT3x8fGanp4e0n2sWLGCYcOGMXPmTJo1axZwmby8PGrUqEGrmFhy4hOpFd3+cF1uVga102eyIXNtSOM0xpjSEpEvVTW+cLk1i5WBnTt30r9/fzp37syiRYt44IEHily2Ro0aAGStz6RmVFyBuppRcWStzwxprMYYEwyWXEJIVUlLSyM2NpapU6ceHgfs5ZdfLvHO+ujWMezLLvh44n3Zq4luHROyeI0xJlgsuYTImjVruPDCC7n11lvZsmWLX/27775b7Ppjk0exd/FUcrMy0IMHyM3KYO/iqYxNHhWqkI0xJmisQz/Ifv/9d8aOHcvEiRPZv3+/X33r1q155plnuPzyy4vdTn6n/cjkFLLmZRLdOoaJkyZYZ74xpkKwDn0nGB36b731Fn/961/58ccf/epq1KjBsGHDGDZsGLVr1z6m/RhjTHlR5h36IjJdRH4RkW8D1D0oIioijd28iMgUEVknIhki0tFn2UQRyXSvRJ/yTiKy0q0zRdwdhiLSSEQWueUXiUjDUL3HtLS5tIqJRapVo85xdenVq1fAxNK9e3e+/fZbxowZY4nFGFMlhLLPZQZwaeFCEWkOXAz4fgpfBsS4171Aqlu2ETAa+APQBRjtkyxS3bL56+XvaxiwWFVjgMVuPujS0uaSNHAwv9Y7DYmoTs7e3/2WOfnkk5k3bx7//ve/iYmxjnhjTNURsuSiqkuB7QGqJgFDAN/2uKuAWepZBjQQkWbAJcAiVd2uqjuARcClrq6+qn6mXrveLKCXz7byn6A106c8qEYmp1CnWz+qN26OHsgrUBcREcGgQYNYu3Yt119/vQ3bYoypcsq0Q19EegI/qeo3hT5wTwE2+sxnu7LiyrMDlAOcqKqbAVR1s4g0DeqbcLLWZ9K8Vxw1m7djT8Yi8rZ8f7huxYoVtG/fvpi1jTGmciuzS5FFpA4wAgh0LW2gr/Z6FOVHGtO9IpIuIulbt249onXz70ORahE06t6HanUaUL/LtUSferpfYsnvm6kWEWFjhBljqoSyvM/lVKAV8I2I/ABEAStE5CS8M4/mPstGAZtKKI8KUA6wxTWb4X7+UlRAqvqiqsaranyTJk2O6M343odS48RTaXz5QHTDZzw6ZnSB5fL7ZnLiE2k+aD458YkkDRxsCcYYU6mVWXJR1ZWq2lRVW6pqS7wE0VFVfwbeBnq7q8bOAXa5pq33gO4i0tB15HcH3nN1u0XkHHeVWG/gLbert4H8q8oSfcqDKiHhFlInTaB2+kw2PnUNdb6eS2qA+1Dy+2ZqRbdHIiKpFd2eOt36MTI5JRRhGWNMuRCy+1xE5GWgK9AY2AKMVtVpPvU/APGqus0liKl4V3ztBe5Q1XS33J3AQ261R1X17648Hu+KtNrAAqC/qqqInADMA1rgXZF2vaoGurCggFANXFktIoLmg+YjEf/r3tKDB9j41DUcOngw6PszxpiyVNR9LiHr0FfVm0uob+kzrUDfIpabDkwPUJ4OnBGg/Feg2xGGGzLRrWPIyV5dYHRjGyPMGFPZ2dhiIWZjhBljqiIbWyzEbIwwY0xVZGOLOWXxsDBjjKls7GFhxhhjyowlF2OMMUFnycUYY0zQWXIxxhgTdJZcjDHGBJ1dLeaIyFYg6yhXbwxsC2I4lYUdF392TAKz4xJYRTgu0arqNzijJZcgEJH0QJfiVXV2XPzZMQnMjktgFfm4WLOYMcaYoLPkYowxJugsuQTHi+EOoJyy4+LPjklgdlwCq7DHxfpcjDHGBJ2duRhjjAk6Sy7GGGOCrkonFxGZLiK/iMi3AeoeFBEVkcZuXkRkioisE5EMEenos2yiiGS6V6JPeScRWenWmeKeuImINBKRRW75Re4RzuVGUcdFRPqLyHciskpExvuUD3fv8TsRucSn/FJXtk5EhvmUtxKRz937f1VEarjymm5+natvGfp3W3qBjouIdBCRZSLytYiki0gXV14l/l5EpLmIfCAia9zfxQBXHjBmOy4yQUTWuvf+pog08Fmncv0fqWqVfQF/AToC3xYqbw68h3dTZWNXdjne45QFOAf43JU3Ata7nw3ddENX9wVwrltnAXCZKx8PDHPTw4Anwn0sSjouwAXA+0BNN9/U/YwDvgFqAq2A74EI9/oeaA3UcMvEuXXmATe56eeBJDfdB3jeTd8EvBruY1GK47LQ5/d6OfBhVfp7AZoBHd10PeC/7m8iYMx2XOgORLryJ3yOS6X7Pwr7LyHcL6Al/snldeAs4Af+l1xeAG72WeY79wd0M/CCT/kLrqwZsNan/PBy+ev6/BF+F+7jUNJxcX/IFwVYbjgw3Gf+PfdBcC7wXuHl3AfENp9/sMPL5a/rpiPdchLuY1HCcXkPuNHndzy3Kv69+MT9FnBxUTFX9eNSqOxqIM1NV7r/oyrdLBaIiPQEflLVbwpVnQJs9JnPdmXFlWcHKAc4UVU3A7ifTYP2BkLndODP7jT7IxHp7MqP9LicAOxU1QOFygtsy9XvcsuXZ/cDE0RkI/Ak3j8+VMG/F9f8cjbwOUXHXNWPi6878c7EoBL+H1ly8SEidYARQKAH3EuAMj2K8ooqEq+54hxgMDDPtX0H87hUxGOWBAxU1ebAQGCaK69Sfy8iUhd4A7hfVX8rbtEAZVXuuIjICOAAkJZfFGD1Cv1/ZMmloFPx2ju/EZEfgChghYichPfNoLnPslHAphLKowKUA2wRkWYA7ucvQX8nwZcNzFfPF8AhvEH1jvS4bAMaiEhkoXJ813H1xwPbQ/JugicRmO+mXwO6uOkq8/ciItXxPkDTVDX/WBQVc1U/LriLFa4AEtS1XVEJ/48sufhQ1ZWq2lRVW6pqS7xfUkdV/Rl4G+jtrnY5B9jlTsXfA7qLSEN3tUp3vLbPzcBuETnHfcPvjdfuittW/tUwiT7l5dk/gAsBROR0vM7FbXjv5SZ3hUorIAavA3Y5EOOuaKmB17H4tvtn+gC4zm3X9/37HpfrgCU+/3zl1SbgfDd9IZDppqvE34uLdRqwRlWf8qkqKuYqfVxE5FJgKNBTVff6rFL5/o/C3dEVzhfwMrAZ2I+XSO4qVP8D/+vQF+BZvCs3VgLxPsvdCaxzrzt8yuOBb906U/nfiAgnAIvxPogWA43CfSxKOi54yWSOez8rgAt9lh/h3uN3uCt5XPnleFfJfA+M8ClvjfePsw7v237+FWi13Pw6V9863MeiFMflT8CXeFfxfA50qkp/L+79K5ABfO1elxcVsx0X1uH1h+SXPV9Z/49s+BdjjDFBZ81ixhhjgs6SizHGmKCz5GKMMSboLLkYY4wJOksuxhhjgs6SizFBIN4I2rN95iNFZKuIvBPOuEoiIh+KSHy44zCVjyUXY4Ljd+AMEant5i8GfgpHID53bRsTNpZcjAmeBUAPN30z3k2XAIjIceI9D2a5iHwlIle58pYi8rGIrHCvP7ryZiKyVLznxHwrIn925Xt8tnmdiMxw0zNE5CkR+QB4opj91RaRV8R7nsirQH4yNCao7BuOMcHzCjDKNYW1B6YDf3Z1I/CG4bhTvAdEfSEi7+ONh3WxquaKSAxeQooHbsEb/uRREYkA6pRi/6fjPRbhoIg8VsT+/g/Yq6rtRaQ93mgLxgSdJRdjgkRVM9zw6jcD7xaq7g70FJEH3XwtoAXe2GRTRaQDcBAvQYA3ptR0N/jhP1T161KE8JqqHixhf38BpvjEm3Fk79KY0rHkYkxwvY33XJeuFHyGhgDXqup3vguLSDKwBe/hdNWAXABVXSoif8FrZpstIhNUdRYFh06vVWjfv5dif1CBhqw3FZf1uRgTXNOBFFVdWaj8PaC/Gy0XETnblR8PbFbVQ8BteI+1RUSigV9U9SW80XXznzW/RUTaikg1vCcZFqWo/S0FElzZGXjNd8YEnSUXY4JIVbNV9ekAVY8A1YEMEfnWzQM8BySKyDK8JrH8s4+uwNci8hVwLZC/zWHAO8ASvBGai1LU/lKBuq45bAjeqLnGBJ2NimyMMSbo7MzFGGNM0FlyMcYYE3SWXIwxxgSdJRdjjDFBZ8nFGGNM0FlyMcYYE3SWXIwxxgTd/wNmMpF742M2rgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a graph of the predicted versus cross-validated errors\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(Y, predicted, edgecolors=(0, 0, 0))\n",
    "ax.plot([Y.min(), Y.max()], [Y.min(), Y.max()], 'k--', lw=4)\n",
    "ax.set_xlabel('Measured')\n",
    "ax.set_ylabel('Predicted')\n",
    "plt.title('Prediction by Cross-Validation', fontdict=None, loc='center', pad=None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method 2 Using Statsmodels\n",
    "import statsmodels.api as sm\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "summary()\n",
      "                             OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                    dep   R-squared:                       0.890\n",
      "Model:                            OLS   Adj. R-squared:                  0.883\n",
      "Method:                 Least Squares   F-statistic:                     137.2\n",
      "Date:                Wed, 17 Jul 2019   Prob (F-statistic):           1.45e-09\n",
      "Time:                        09:08:05   Log-Likelihood:                -201.11\n",
      "No. Observations:                  19   AIC:                             406.2\n",
      "Df Residuals:                      17   BIC:                             408.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const      -3.335e+04   1.96e+04     -1.698      0.108   -7.48e+04    8092.445\n",
      "caputil     4202.7151    358.786     11.714      0.000    3445.742    4959.688\n",
      "==============================================================================\n",
      "Omnibus:                        2.344   Durbin-Watson:                   0.634\n",
      "Prob(Omnibus):                  0.310   Jarque-Bera (JB):                1.898\n",
      "Skew:                          -0.721   Prob(JB):                        0.387\n",
      "Kurtosis:                       2.437   Cond. No.                         464.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "pvalues\n",
      " const      1.077598e-01\n",
      "caputil    1.454602e-09\n",
      "dtype: float64\n",
      "tvalues\n",
      " const      -1.697874\n",
      "caputil    11.713701\n",
      "dtype: float64\n",
      "rsquared\n",
      " 0.8897612752328166\n",
      "rsquared_adj\n",
      " 0.8832766443641588\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Tony Diana\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.\n",
      "  return ptp(axis=axis, out=out, **kwargs)\n",
      "C:\\Users\\Tony Diana\\anaconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1416: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=19\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "# Add a constant\n",
    "X2 = sm.add_constant(X) \n",
    "# Define the model\n",
    "est = sm.OLS(Y, X2)\n",
    "# Fit the model\n",
    "est2 = est.fit()\n",
    "print(\"summary()\\n\",est2.summary())\n",
    "print(\"pvalues\\n\",est2.pvalues)\n",
    "print(\"tvalues\\n\",est2.tvalues)\n",
    "print(\"rsquared\\n\",est2.rsquared)\n",
    "print(\"rsquared_adj\\n\",est2.rsquared_adj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0     162786.620151\n",
      "1     146480.085616\n",
      "2     139503.578572\n",
      "3     141436.827512\n",
      "4     178126.530216\n",
      "5     196660.503747\n",
      "6     207377.427217\n",
      "7     228979.382761\n",
      "8     220784.088343\n",
      "9     207671.617273\n",
      "10    212967.038282\n",
      "11    206158.639842\n",
      "12    194391.037600\n",
      "13    196786.585199\n",
      "14    205612.286881\n",
      "15    216497.318955\n",
      "16    211958.386661\n",
      "17    216581.373256\n",
      "18    217043.671916\n",
      "dtype: float64\n",
      "r2_score 0.8897612752328166\n"
     ]
    }
   ],
   "source": [
    "# Determine predicted values and R-squared with scikitlearn\n",
    "from sklearn.metrics import r2_score\n",
    "predictions = est2.predict(X2)\n",
    "\n",
    "print(est2.predict(X2))\n",
    "\n",
    "print(\"r2_score\",r2_score(Y,predictions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "scale=StandardScaler()\n",
    "# Load the file\n",
    "airport=pd.read_csv('D:/Machine Learning/airport.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a data frame\n",
    "airport=pd.DataFrame(airport, columns=['caputil', 'dep'])\n",
    "# Define X and Y\n",
    "X=pd.DataFrame(airport['caputil'])\n",
    "Y=pd.DataFrame(airport['dep'])"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}