blog_notebooks/regression/regression_blog_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import cudf \n",
    "import cuml\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This data description, and the data that this example uses are available at \n",
    "[the UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "creating bike data directory\n"
     ]
    }
   ],
   "source": [
    "import urllib.request\n",
    "\n",
    "data_dir = '../../data/bike/'\n",
    "if not os.path.exists(data_dir):\n",
    "    print('creating bike data directory')\n",
    "    os.system('mkdir ../../data/bike')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading https://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip to ../../data/bike/Bike-Sharing-Dataset.zip\n"
     ]
    }
   ],
   "source": [
    "base_url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00275/'\n",
    "fn = 'Bike-Sharing-Dataset.zip'\n",
    "if not os.path.isfile(data_dir+fn):\n",
    "        print(f'Downloading {base_url+fn} to {data_dir+fn}')\n",
    "        urllib.request.urlretrieve(base_url+fn, data_dir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Files extracted: ['Readme.txt', 'day.csv', 'hour.csv']\n"
     ]
    }
   ],
   "source": [
    "from zipfile import ZipFile \n",
    "files = []\n",
    "with ZipFile(data_dir+fn) as myzip:\n",
    "    files = myzip.namelist()\n",
    "    print(\"Files extracted: \"+ str(files))\n",
    "    myzip.extractall(data_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Description: \n",
    "#         \n",
    "#         - instant: record index\n",
    "#         - dteday : date\n",
    "#         - season : season (1:springer, 2:summer, 3:fall, 4:winter)\n",
    "#         - yr : year (0: 2011, 1:2012)\n",
    "#         - mnth : month ( 1 to 12)\n",
    "#         - hr : hour (0 to 23)\n",
    "#         - holiday : weather day is holiday or not (extracted from http://dchr.dc.gov/page/holiday-schedule)\n",
    "#         - weekday : day of the week\n",
    "#         - workingday : if day is neither weekend nor holiday is 1, otherwise is 0.\n",
    "#         + weathersit : \n",
    "#                 - 1: Clear, Few clouds, Partly cloudy, Partly cloudy\n",
    "#                 - 2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist\n",
    "#                 - 3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds\n",
    "#                 - 4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog\n",
    "#         - temp : Normalized temperature in Celsius. The values are divided to 41 (max)\n",
    "#         - atemp: Normalized feeling temperature in Celsius. The values are divided to 50 (max)\n",
    "#         - hum: Normalized humidity. The values are divided to 100 (max)\n",
    "#         - This is a good example of variables we might not have a good theoretical intuation around\n",
    "#         - windspeed: Normalized wind speed. The values are divided to 67 (max)\n",
    "#         - casual: count of casual users\n",
    "#         - registered: count of registered users\n",
    "#         - cnt: count of total rental bikes including both casual and registered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use `files[2]` and read in the hour.csv dataset.  Valid datasets are `files[1]` (day.csv) and `files[2]` (hour.csv).   Please do not use `files[0]`, which is a README file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = cudf.read_csv(data_dir+files[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We drop the index, the timestamp (because they have broken it down for us), \n",
    "and the individual counts that make up our target variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "drop_list = ['instant', 'dteday', 'casual', 'registered']\n",
    "gdf = gdf.drop(drop_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to create one-hot encoded variables, also known as dummy variables, for each of the time variables as well as the weather situation. We're going to drop one of each of these dummy variables so we don't create colinearity. \n",
    "\n",
    "The next data munging step we take is to convert all of our data into the same type, because that is what the cuML algorithms are expecting. \n",
    "\n",
    "Last, we split our data into test and train sets, training on 2011 data, and testing on 2012. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dummies_list = ['season','yr', 'mnth', 'hr', 'weekday', 'weathersit']\n",
    "\n",
    "for item in dummies_list:\n",
    "    codes = gdf[item].unique()\n",
    "    gdf = gdf.one_hot_encoding(item, '{}_dummy'.format(item), codes)\n",
    "    gdf = gdf.drop('{}_dummy_1'.format(item))\n",
    "\n",
    "#cuML current requires all data be of the same type, so this loop converts all values into floats\n",
    "for col in gdf.columns.tolist():\n",
    "    gdf[col] = gdf[col].astype('float32')\n",
    "    \n",
    "test = gdf.query('yr == 1').drop(dummies_list)\n",
    "train = gdf.query('yr == 0').drop(dummies_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I am going to test out how well a small variable does against all the variables available. I select \"weathersit_dummy_4\" because as we see above, that seems like the worst weather for bike riding. Also, based on personal experience, bike riding in high wind is not fun either. I add workday because I'm sure it has some impact, but I'm not sure what. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = train['cnt']\n",
    "y_test = test['cnt']\n",
    "\n",
    "#Some of the variables, chosen by theory\n",
    "X_train_1 = train[['weathersit_dummy_4', 'windspeed', 'workingday']]\n",
    "X_test_1 = test[['weathersit_dummy_4', 'windspeed', 'workingday']]\n",
    "\n",
    "#all of the varibles.\n",
    "X_train_2 = train.drop('cnt')\n",
    "X_test_2 = test.drop('cnt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, I run two regressions. The first is based on small set of variable I think will be most impactful to bike ridership. The second regression inclueds all the variables availale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "OLS_1 = cuml.LinearRegression()\n",
    "fit_1 = OLS_1.fit(X_train_1, y_train)\n",
    "y_hat_1 = fit_1.predict(X_test_1)\n",
    "MSE_1 = ((y_test - y_hat_1)**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "OLS_2 = cuml.LinearRegression()\n",
    "fit_2 = OLS_2.fit(X_train_2, y_train)\n",
    "y_hat_2 = fit_2.predict(X_test_2)\n",
    "MSE_2 = ((y_test - y_hat_2)**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = {'MSE_OLS_THEORY':MSE_1,\n",
    "         'MSE_OLS_ALL': MSE_2}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like I was wrong, and the model with everything performs better on the out-of-sample data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE_OLS_THEORY:     827127.25\n",
      "MSE_OLS_EVERYTHING: 1321845.875\n"
     ]
    }
   ],
   "source": [
    "print('MSE_OLS_THEORY:     {}'.format(MSE_1))\n",
    "print('MSE_OLS_EVERYTHING: {}'.format(MSE_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"everything\" model outperformed the small model, but let's see if we can do better by doing a Ridge regression. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to do a small hyperparameter search for alpha, checking 100 different values. This is fast to do with RAPIDS. Also notice that I am appending the results of each Ridge model onto the dictionary containing our earlier results, so I can more easily see which model is the best at the end. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "for alpha in np.arange(0.01, 1, 0.01): #alpha value has to be positive\n",
    "    \n",
    "    Ridge = cuml.Ridge(alpha=alpha, fit_intercept=True)\n",
    "    fit_3 = Ridge.fit(X_train_2, y_train)\n",
    "    y_hat_3 = fit_3.predict(X_test_2)\n",
    "    MSE_3 = ((y_test - y_hat_3)**2).sum()\n",
    "    output['MSE_RIDGE_{}'.format(alpha)] = MSE_3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that our regulaized model does better than the rest, include OLS with all the variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min MSE: MSE_OLS_THEORY\n"
     ]
    }
   ],
   "source": [
    "print('Min MSE: {}'.format(min(output, key=output.get)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'MSE_OLS_ALL': 1321845.875,\n",
      " 'MSE_OLS_THEORY': 827127.25,\n",
      " 'MSE_RIDGE_0.01': 1321725.25,\n",
      " 'MSE_RIDGE_0.02': 1321610.25,\n",
      " 'MSE_RIDGE_0.03': 1321495.25,\n",
      " 'MSE_RIDGE_0.04': 1321380.375,\n",
      " 'MSE_RIDGE_0.05': 1321266.0,\n",
      " 'MSE_RIDGE_0.060000000000000005': 1321151.5,\n",
      " 'MSE_RIDGE_0.06999999999999999': 1321037.25,\n",
      " 'MSE_RIDGE_0.08': 1320923.5,\n",
      " 'MSE_RIDGE_0.09': 1320809.75,\n",
      " 'MSE_RIDGE_0.09999999999999999': 1320696.0,\n",
      " 'MSE_RIDGE_0.11': 1320582.625,\n",
      " 'MSE_RIDGE_0.12': 1320469.5,\n",
      " 'MSE_RIDGE_0.13': 1320356.375,\n",
      " 'MSE_RIDGE_0.14': 1320243.5,\n",
      " 'MSE_RIDGE_0.15000000000000002': 1320130.875,\n",
      " 'MSE_RIDGE_0.16': 1320018.5,\n",
      " 'MSE_RIDGE_0.17': 1319906.0,\n",
      " 'MSE_RIDGE_0.18000000000000002': 1319794.0,\n",
      " 'MSE_RIDGE_0.19': 1319681.875,\n",
      " 'MSE_RIDGE_0.2': 1319570.25,\n",
      " 'MSE_RIDGE_0.21000000000000002': 1319458.625,\n",
      " 'MSE_RIDGE_0.22': 1319347.25,\n",
      " 'MSE_RIDGE_0.23': 1319235.875,\n",
      " 'MSE_RIDGE_0.24000000000000002': 1319124.875,\n",
      " 'MSE_RIDGE_0.25': 1319013.875,\n",
      " 'MSE_RIDGE_0.26': 1318903.0,\n",
      " 'MSE_RIDGE_0.27': 1318792.5,\n",
      " 'MSE_RIDGE_0.28': 1318682.0,\n",
      " 'MSE_RIDGE_0.29000000000000004': 1318571.75,\n",
      " 'MSE_RIDGE_0.3': 1318461.75,\n",
      " 'MSE_RIDGE_0.31': 1318351.75,\n",
      " 'MSE_RIDGE_0.32': 1318242.0,\n",
      " 'MSE_RIDGE_0.33': 1318132.375,\n",
      " 'MSE_RIDGE_0.34': 1318022.875,\n",
      " 'MSE_RIDGE_0.35000000000000003': 1317913.375,\n",
      " 'MSE_RIDGE_0.36000000000000004': 1317804.5,\n",
      " 'MSE_RIDGE_0.37': 1317695.375,\n",
      " 'MSE_RIDGE_0.38': 1317586.5,\n",
      " 'MSE_RIDGE_0.39': 1317477.75,\n",
      " 'MSE_RIDGE_0.4': 1317369.25,\n",
      " 'MSE_RIDGE_0.41000000000000003': 1317261.0,\n",
      " 'MSE_RIDGE_0.42000000000000004': 1317152.75,\n",
      " 'MSE_RIDGE_0.43': 1317044.875,\n",
      " 'MSE_RIDGE_0.44': 1316936.75,\n",
      " 'MSE_RIDGE_0.45': 1316829.0,\n",
      " 'MSE_RIDGE_0.46': 1316721.25,\n",
      " 'MSE_RIDGE_0.47000000000000003': 1316613.875,\n",
      " 'MSE_RIDGE_0.48000000000000004': 1316506.625,\n",
      " 'MSE_RIDGE_0.49': 1316399.5,\n",
      " 'MSE_RIDGE_0.5': 1316292.25,\n",
      " 'MSE_RIDGE_0.51': 1316185.5,\n",
      " 'MSE_RIDGE_0.52': 1316078.75,\n",
      " 'MSE_RIDGE_0.53': 1315972.25,\n",
      " 'MSE_RIDGE_0.54': 1315865.625,\n",
      " 'MSE_RIDGE_0.55': 1315759.5,\n",
      " 'MSE_RIDGE_0.56': 1315653.0,\n",
      " 'MSE_RIDGE_0.5700000000000001': 1315547.25,\n",
      " 'MSE_RIDGE_0.5800000000000001': 1315441.25,\n",
      " 'MSE_RIDGE_0.59': 1315335.5,\n",
      " 'MSE_RIDGE_0.6': 1315230.0,\n",
      " 'MSE_RIDGE_0.61': 1315124.5,\n",
      " 'MSE_RIDGE_0.62': 1315019.25,\n",
      " 'MSE_RIDGE_0.63': 1314913.75,\n",
      " 'MSE_RIDGE_0.64': 1314808.875,\n",
      " 'MSE_RIDGE_0.65': 1314704.0,\n",
      " 'MSE_RIDGE_0.66': 1314599.125,\n",
      " 'MSE_RIDGE_0.67': 1314494.5,\n",
      " 'MSE_RIDGE_0.68': 1314390.0,\n",
      " 'MSE_RIDGE_0.6900000000000001': 1314285.75,\n",
      " 'MSE_RIDGE_0.7000000000000001': 1314181.5,\n",
      " 'MSE_RIDGE_0.7100000000000001': 1314077.25,\n",
      " 'MSE_RIDGE_0.72': 1313973.25,\n",
      " 'MSE_RIDGE_0.73': 1313869.5,\n",
      " 'MSE_RIDGE_0.74': 1313765.875,\n",
      " 'MSE_RIDGE_0.75': 1313662.125,\n",
      " 'MSE_RIDGE_0.76': 1313558.75,\n",
      " 'MSE_RIDGE_0.77': 1313455.5,\n",
      " 'MSE_RIDGE_0.78': 1313352.125,\n",
      " 'MSE_RIDGE_0.79': 1313249.125,\n",
      " 'MSE_RIDGE_0.8': 1313146.25,\n",
      " 'MSE_RIDGE_0.81': 1313043.375,\n",
      " 'MSE_RIDGE_0.8200000000000001': 1312940.875,\n",
      " 'MSE_RIDGE_0.8300000000000001': 1312838.25,\n",
      " 'MSE_RIDGE_0.8400000000000001': 1312735.875,\n",
      " 'MSE_RIDGE_0.85': 1312633.5,\n",
      " 'MSE_RIDGE_0.86': 1312531.5,\n",
      " 'MSE_RIDGE_0.87': 1312429.375,\n",
      " 'MSE_RIDGE_0.88': 1312327.375,\n",
      " 'MSE_RIDGE_0.89': 1312225.625,\n",
      " 'MSE_RIDGE_0.9': 1312124.0,\n",
      " 'MSE_RIDGE_0.91': 1312022.5,\n",
      " 'MSE_RIDGE_0.92': 1311921.0,\n",
      " 'MSE_RIDGE_0.93': 1311819.75,\n",
      " 'MSE_RIDGE_0.9400000000000001': 1311718.5,\n",
      " 'MSE_RIDGE_0.9500000000000001': 1311617.5,\n",
      " 'MSE_RIDGE_0.9600000000000001': 1311516.625,\n",
      " 'MSE_RIDGE_0.97': 1311415.875,\n",
      " 'MSE_RIDGE_0.98': 1311315.0,\n",
      " 'MSE_RIDGE_0.99': 1311214.625}\n"
     ]
    }
   ],
   "source": [
    "#You can uncomment the code below to see all 100 model MSEs\n",
    "#Notice in particular that MSE for OLS with everything and Ridge with alpha = 0 are essentially the same. \n",
    "import pprint\n",
    "pprint.pprint(output)"
   ]
  }
 ],
 "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": 4
}