02 - Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression\n",
    "\n",
    "*Supervised* machine learning techniques involve training a model to operate on a set of *features* and predict a *label* using a dataset that includes some already-known label values. The training process *fits* the features to the known labels to define a general function that can be applied to new features for which the labels are unknown, and predict them. You can think of this function like this, in which ***y*** represents the label we want to predict and ***x*** represents the features the model uses to predict it.\n",
    "\n",
    "$$y = f(x)$$\n",
    "\n",
    "In most cases, *x* is actually a *vector* that consists of multiple feature values, so to be a little more precise, the function could be expressed like this:\n",
    "\n",
    "$$y = f([x_1, x_2, x_3, ...])$$\n",
    "\n",
    "The goal of training the model is to find a function that performs some kind of calcuation to the *x* values that produces the result *y*. We do this by applying a machine learning *algorithm* that tries to fit the *x* values to a calculation that produces *y* reasonably accurately for all of the cases in the training dataset.\n",
    "\n",
    "There are lots of machine learning algorithms for supervised learning, and they can be broadly divided into two types:\n",
    "\n",
    "- **_Regression_ algorithms**: Algorithms that predict a *y* value that is a numeric value, such as the price of a house or the number of sales transactions.\n",
    "- **_Classification_ algorithms**: Algorithms that predict to which category, or *class*, an observation belongs. The *y* value in a classification model is a vector of probability values between 0 and 1, one for each class, indicating the probability of the observation belonging to each class.\n",
    "\n",
    "In this notebook, we'll focus on *regression*, using an example based on a real study in which data for a bicycle sharing scheme was collected and used to predict the number of rentals based on seasonality and weather conditions. We'll use a simplified version of the dataset from that study.\n",
    "\n",
    "## Explore the Data\n",
    "\n",
    "The first step in any machine learning project is to explore the data that you will use to train a model. The goal of this exploration is to try to understand the relationships between its attributes; in particular, any apparent correlation between the *features* and the *label* your model will try to predict. This may require some work to detect and fix issues in the data (such as dealing with missing values, errors, or outlier values), deriving new feature columns by transforming or combining existing features (a process known as *feature engineering*), *normalizing* numeric features (values you can measure or count) so they're on a similar scale, and *encoding* categorical features (values that represent discrete categories) as numeric indicators.\n",
    "\n",
    "Let's start by loading the bicycle sharing data as a **Pandas** DataFrame and viewing the first few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# load the training dataset\n",
    "bike_data = pd.read_csv('data/daily-bike-share.csv')\n",
    "bike_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data consists of the following columns:\n",
    "\n",
    "- **instant**: A unique row identifier\n",
    "- **dteday**: The date on which the data was observed - in this case, the data was collected daily; so there's one row per date.\n",
    "- **season**: A numerically encoded value indicating the season (1:spring, 2:summer, 3:fall, 4:winter)\n",
    "- **yr**: The year of the study in which the observation was made (the study took place over two years - year 0 represents 2011, and year 1 represents 2012)\n",
    "- **mnth**: The calendar month in which the observation was made (1:January ... 12:December)\n",
    "- **holiday**: A binary value indicating whether or not the observation was made on a public holiday)\n",
    "- **weekday**: The day of the week on which the observation was made (0:Sunday ... 6:Saturday)\n",
    "- **workingday**: A binary value indicating whether or not the day is a working day (not a weekend or holiday)\n",
    "- **weathersit**: A categorical value indicating the weather situation (1:clear, 2:mist/cloud, 3:light rain/snow, 4:heavy rain/hail/snow/fog)\n",
    "- **temp**: The temperature in celsius (normalized)\n",
    "- **atemp**: The apparent (\"feels-like\") temperature in celsius (normalized)\n",
    "- **hum**: The humidity level (normalized)\n",
    "- **windspeed**: The windspeed (normalized)\n",
    "- **rentals**: The number of bicycle rentals recorded.\n",
    "\n",
    "In this dataset, **rentals** represents the label (the *y* value) our model must be trained to predict. The other columns are potential features (*x* values).\n",
    "\n",
    "As mentioned previously, you can perform some *feature engineering* to combine or derive new features. For example, let's add a new column named **day** to the dataframe by extracting the day component from the existing **dteday** column. The new column represents the day of the month from 1 to 31."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bike_data['day'] = pd.DatetimeIndex(bike_data['dteday']).day\n",
    "bike_data.head(32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, let's start our analysis of the data by examining a few key descriptive statistics. We can use the dataframe's **describe** method to generate these for the numeric features as well as the **rentals** label column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numeric_features = ['temp', 'atemp', 'hum', 'windspeed']\n",
    "bike_data[numeric_features + ['rentals']].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statistics reveal some information about the distribution of the data in each of the numeric fields, including the number of observations (there are 731 records), the mean, standard deviation, minimum and maximum values, and the quartile values (the threshold values for 25%, 50% - which is also the median, and 75% of the data). From this, we can see that the mean number of daily rentals is around 848; but there's a comparatively large standard deviation, indicating a lot of variance in the number of rentals per day.\n",
    "\n",
    "We might get a clearer idea of the distribution of rentals values by visualizing the data. Common plot types for visualizing numeric data distributions are *histograms* and *box plots*, so let's use Python's **matplotlib** library to create one of each of these for the **rentals** column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# This ensures plots are displayed inline in the Jupyter notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Get the label column\n",
    "label = bike_data['rentals']\n",
    "\n",
    "\n",
    "# Create a figure for 2 subplots (2 rows, 1 column)\n",
    "fig, ax = plt.subplots(2, 1, figsize = (9,12))\n",
    "\n",
    "# Plot the histogram   \n",
    "ax[0].hist(label, bins=100)\n",
    "ax[0].set_ylabel('Frequency')\n",
    "\n",
    "# Add lines for the mean, median, and mode\n",
    "ax[0].axvline(label.mean(), color='magenta', linestyle='dashed', linewidth=2)\n",
    "ax[0].axvline(label.median(), color='cyan', linestyle='dashed', linewidth=2)\n",
    "\n",
    "# Plot the boxplot   \n",
    "ax[1].boxplot(label, vert=False)\n",
    "ax[1].set_xlabel('Rentals')\n",
    "\n",
    "# Add a title to the Figure\n",
    "fig.suptitle('Rental Distribution')\n",
    "\n",
    "# Show the figure\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show that the number of daily rentals ranges from 0 to just over 3,400. However, the mean (and median) number of daily rentals is closer to the low end of that range, with most of the data between 0 and around 2,200 rentals. The few values above this are shown in the box plot as small circles, indicating that they are *outliers* - in other words, unusually high or low values beyond the typical range of most of the data.\n",
    "\n",
    "We can do the same kind of visual exploration of the numeric features. Let's create a histogram for each of these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot a histogram for each numeric feature\n",
    "for col in numeric_features:\n",
    "    fig = plt.figure(figsize=(9, 6))\n",
    "    ax = fig.gca()\n",
    "    feature = bike_data[col]\n",
    "    feature.hist(bins=100, ax = ax)\n",
    "    ax.axvline(feature.mean(), color='magenta', linestyle='dashed', linewidth=2)\n",
    "    ax.axvline(feature.median(), color='cyan', linestyle='dashed', linewidth=2)\n",
    "    ax.set_title(col)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numeric features seem to be more *normally* distributed, with the mean and median nearer the middle of the range of values, coinciding with where the most commonly occurring values are.\n",
    "\n",
    "> **Note**: The distributions are not truly *normal* in the statistical sense, which would result in a smooth, symmetric \"bell-curve\" histogram with the mean and mode (the most common value) in the center; but they do generally indicate that most of the observations have a value somewhere near the middle.\n",
    "\n",
    "We've explored the distribution of the numeric values in the dataset, but what about the categorical features? These aren't continuous numbers on a scale, so we can't use histograms; but we can plot a bar chart showing the count of each discrete value for each category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# plot a bar plot for each categorical feature count\n",
    "categorical_features = ['season','mnth','holiday','weekday','workingday','weathersit', 'day']\n",
    "\n",
    "for col in categorical_features:\n",
    "    counts = bike_data[col].value_counts().sort_index()\n",
    "    fig = plt.figure(figsize=(9, 6))\n",
    "    ax = fig.gca()\n",
    "    counts.plot.bar(ax = ax, color='steelblue')\n",
    "    ax.set_title(col + ' counts')\n",
    "    ax.set_xlabel(col) \n",
    "    ax.set_ylabel(\"Frequency\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many of the categorical features show a more or less *uniform* distribution (meaning there's roughly the same number of rows for each category). Exceptions to this include:\n",
    "\n",
    "- **holiday**: There are many fewer days that are holidays than days that aren't.\n",
    "- **workingday**: There are more working days than non-working days.\n",
    "- **weathersit**: Most days are category *1* (clear), with category *2* (mist and cloud) the next most common. There are comparitively few category *3* (light rain or snow) days, and no category *4* (heavy rain, hail, or fog) days at all.\n",
    "\n",
    "Now that we know something about the distribution of the data in our columns, we can start to look for relationships between the features and the **rentals** label we want to be able to predict.\n",
    "\n",
    "For the numeric features, we can create scatter plots that show the intersection of feature and label values. We can also calculate the *correlation* statistic to quantify the apparent relationship.. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in numeric_features:\n",
    "    fig = plt.figure(figsize=(9, 6))\n",
    "    ax = fig.gca()\n",
    "    feature = bike_data[col]\n",
    "    label = bike_data['rentals']\n",
    "    correlation = feature.corr(label)\n",
    "    plt.scatter(x=feature, y=label)\n",
    "    plt.xlabel(col)\n",
    "    plt.ylabel('Bike Rentals')\n",
    "    ax.set_title('rentals vs ' + col + '- correlation: ' + str(correlation))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results aren't conclusive, but if you look closely at the scatter plots for **temp** and **atemp**, you can see a vague diagonal trend showing that higher rental counts tend to coincide with higher temperatures; and a correlation value of just over 0.5 for both of these features supports this observation. Conversely, the plots for **hum** and **windspeed** show a slightly negative correlation, indicating that there are fewer rentals on days with high humidity or windspeed.\n",
    "\n",
    "Now let's compare the categorical features to the label. We'll do this by creating box plots that show the distribution of rental counts for each category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# plot a boxplot for the label by each categorical feature\n",
    "for col in categorical_features:\n",
    "    fig = plt.figure(figsize=(9, 6))\n",
    "    ax = fig.gca()\n",
    "    bike_data.boxplot(column = 'rentals', by = col, ax = ax)\n",
    "    ax.set_title('Label by ' + col)\n",
    "    ax.set_ylabel(\"Bike Rentals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show some variance in the relationship between some category values and rentals. For example, there's a clear difference in the distribution of rentals on weekends (**weekday** 0 or 6) and those during the working week (**weekday** 1 to 5). Similarly, there are notable differences for **holiday** and **workingday** categories. There's a noticable trend that shows different rental distributions in summer and fall months compared to spring and winter months. The **weathersit** category also seems to make a difference in rental distribution. The **day** feature we created for the day of the month shows little variation, indicating that it's probably not predictive of the number of rentals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train a Regression Model\n",
    "\n",
    "Now that we've explored the data, it's time to use it to train a regression model that uses the features we've identified as potentially predictive to predict the **rentals** label.  The first thing we need to do is to separate the features we want to use to train the model from the label we want it to predict."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate features and labels\n",
    "X, y = bike_data[['season','mnth', 'holiday','weekday','workingday','weathersit','temp', 'atemp', 'hum', 'windspeed']].values, bike_data['rentals'].values\n",
    "print('Features:',X[:10], '\\nLabels:', y[:10], sep='\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After separating the dataset, we now have numpy arrays named **X** containing the features, and **y** containing the labels.\n",
    "\n",
    "We *could* train a model using all of the data; but it's common practice in supervised learning to split the data into two subsets; a (typically larger) set with which to train the model, and a smaller \"hold-back\" set with which to validate the trained model. This enables us to evaluate how well the model performs when used with the validation dataset by comparing the predicted labels to the known labels.\n",
    "\n",
    "To randomly split the data, we'll use the **train_test_split** function in the **scikit-learn** library. This library is one of the most widely used machine learning packages for Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Split data 70%-30% into training set and test set\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=0)\n",
    "\n",
    "print ('Training Set: %d, rows\\nTest Set: %d rows' % (X_train.size, X_test.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the following four datasets:\n",
    "\n",
    "- **X_train**: The feature values we'll use to train the model\n",
    "- **y_train**: The corresponding labels we'll use to train the model\n",
    "- **X_test**: The feature values we'll use to validate the model\n",
    "- **y_test**: The corresponding labels we'll use to validate the model\n",
    "\n",
    "Now we're ready to train a model by fitting a suitable regression algorithm to the training data. We'll use the a *linear regression* algorithm, a common starting point for regression that works by trying to find a linear relationship between the *X* values and the *y* label. The resulting model is a function that conceptually defines a line where every possible X and y value combination intersect.\n",
    "\n",
    "In Scikit-Learn, training algorithms are encapsulated in *estimators*, and in this case we'll use the **LinearRegression** estimator to train a linear regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the model\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# Fit a linear regression model on the training set\n",
    "model = LinearRegression().fit(X_train, y_train)\n",
    "print (model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate the Trained Model\n",
    "\n",
    "Now that we've trained the model, we can use it to predict rental counts for the features we held back in our validation dataset. Then we can compare these predictions to the actual label values to evaluate how well (or not!) the model is working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "predictions = model.predict(X_test)\n",
    "np.set_printoptions(suppress=True)\n",
    "print('Predicted labels: ', np.round(predictions)[:10])\n",
    "print('Actual labels   : ' ,y_test[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing each prediction with its corresponding \"ground truth\" actual value isn't a very efficient way to determine how well the model is predicting. Let's see if we can get a better indication by visualizing a scatter plot that compares the predictions to the actual labels. We'll also overlay a trend line to get a general sense for how well the predicted labels align with the true labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's a definite diagonal trend, and the intersections of the predicted and actual values are generally following the path of the trend line; but there's a fair amount of difference between the ideal function represented by the line and the results. This variance represents the *residuals* of the model - in other words, the difference between the label predicted when the model applies the coefficients it learned during training to the validation data, and the actual value of the validation label. These residuals when evaluated from the validation data indicate the expected level of *error* when the model is used with new data for which the label is unknown.\n",
    "\n",
    "You can quantify the residuals by calculating a number of commonly used evaluation metrics. We'll focus on the following three:\n",
    "\n",
    "- **Mean Square Error (MSE)**: The mean of the squared differences between predicted and actual values. This yields a relative metric in which the smaller the value, the better the fit of the model\n",
    "- **Root Mean Square Error (RMSE)**: The square root of the MSE. This yields an absolute metric in the same unit as the label (in this case, numbers of rentals). The smaller the value, the better the model (in a simplistic sense, it represents the average number of rentals by which the predictions are wrong!)\n",
    "- **Coefficient of Determination (usually known as *R-squared* or R<sup>2</sup>**: A relative metric in which the higher the value, the better the fit of the model. In essence, this metric represents how much of the variance between predicted and actual label values the model is able to explain.\n",
    "\n",
    "> **Note**: You can find out more about these and other metrics for evaluating regression models in the [Scikit-Learn documentation](https://scikit-learn.org/stable/modules/model_evaluation.html#regression-metrics)\n",
    "\n",
    "Let's use Scikit-Learn to calculate these metrics for our model, based on the predictions it generated for the validation data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we've quantified the ability of our model to predict the number of rentals. It definitely has *some* predictive power, but we can probably do better!\n",
    "\n",
    "## Experiment with Algorithms\n",
    "\n",
    "The linear regression algorithm we used to train the model has some predictive capability, but there are many kinds of regression algorithm we could try, including:\n",
    "\n",
    "- **Linear algorithms**: Not just the Linear Regression algorithm we used above (which is technically an *Ordinary Least Squares* algorithm), but other variants such as *Lasso* and *Ridge*.\n",
    "- **Tree-based algorithms**: Algorithms that build a decision tree to reach a prediction.\n",
    "- **Ensemble algorithms**: Algorithms that combine the outputs of multiple base algorithms to improve generalizability.\n",
    "\n",
    "> **Note**: For a full list of Scikit-Learn estimators that encapsulate algorithms for supervised machine learning, see the [Scikit-Learn documentation](https://scikit-learn.org/stable/supervised_learning.html). There are many algorithms to choose from, but for most real-world scenarios, the [Scikit-Learn estimator cheat sheet](https://scikit-learn.org/stable/tutorial/machine_learning_map/index.html) can help you find a suitable starting point. \n",
    "\n",
    "### Try Another Linear Algorithm\n",
    "\n",
    "Let's try training our regression model by using a **Lasso** algorithm. We can do this by just changing the estimator in the training code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "# Fit a lasso model on the training set\n",
    "model = Lasso().fit(X_train, y_train)\n",
    "print (model, \"\\n\")\n",
    "\n",
    "# Evaluate the model using the test data\n",
    "predictions = model.predict(X_test)\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Try a Decision Tree Algorithm\n",
    "\n",
    "As an alternative to a linear model, there's a category of algorithms for machine learning that uses a tree-based approach in which the features in the dataset are examined in a series of evaluations, each of which results in a *branch* in a *decision tree* based on the feature value. At the end of each series of branches are leaf-nodes with the predicted label value based on the feature values.\n",
    "\n",
    "It's easiest to see how this works with an example. Let's train a Decision Tree regression model using the bike rental data. After training the model, the code below will print the model definition and a text representation of the tree it uses to predict label values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.tree import export_text\n",
    "\n",
    "# Train the model\n",
    "model = DecisionTreeRegressor().fit(X_train, y_train)\n",
    "print (model, \"\\n\")\n",
    "\n",
    "# Visualize the model tree\n",
    "tree = export_text(model)\n",
    "print(tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we have a tree-based model; but is it any good? Let's evaluate it with the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the model using the test data\n",
    "predictions = model.predict(X_test)\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tree-based model doesn't seem to have improved over the linear model, so what else could we try?\n",
    "\n",
    "### Try an Ensemble Algorithm\n",
    "\n",
    "Ensemble algorithms work by combining multiple base estimators to produce an optimal model, either by appying an aggregate function to a collection of base models (sometimes referred to a *bagging*) or by building a sequence of models that build on one another to improve predictive performance (referred to as *boosting*).\n",
    "\n",
    "For example, let's try a Random Forest model, which applies an averaging function to multiple Decision Tree models for a better overall model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "# Train the model\n",
    "model = RandomForestRegressor().fit(X_train, y_train)\n",
    "print (model, \"\\n\")\n",
    "\n",
    "# Evaluate the model using the test data\n",
    "predictions = model.predict(X_test)\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For good measure, let's also try a *boosting* ensemble algorithm. We'll use a Gradient Boosting estimator, which like a Random Forest algorithm builds multiple trees, but instead of building them all independently and taking the average result, each tree is built on the outputs of the previous one in an attempt to incrementally reduce the *loss* (error) in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the model\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "\n",
    "# Fit a lasso model on the training set\n",
    "model = GradientBoostingRegressor().fit(X_train, y_train)\n",
    "print (model, \"\\n\")\n",
    "\n",
    "# Evaluate the model using the test data\n",
    "predictions = model.predict(X_test)\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimize Hyperparameters\n",
    "\n",
    "Take a look at the **GradientBoostingRegressor** estimator definition in the output above, and note that it, like the other estimators we tried previously, includes a large number of parameters that control the way the model is trained. In machine learning, the term *parameters* refers to values that can be determined from data; values that you specify to affect the behavior of a training algorithm are more correctly referred to as *hyperparameters*.\n",
    "\n",
    "The specific hyperparameters for an estimator vary based on the algorithm that the estimator encapsulates. In the case of the **GradientBoostingRegressor** estimator, the algorithm is an ensemble that combines multiple decision trees to create an overall predictive model. You can learn about the hyperparameters for this estimator in the [Scikit-Learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html).\n",
    "\n",
    "We won't go into the details of each hyperparameter here, but they work together to affect the way the algorithm trains a model. In many cases, the default values provided by Scikit-Learn will work well; but there may be some advantage in modifying hyperparameters to get better predictive performance or reduce training time.\n",
    "\n",
    "So how do you know what hyperparameter values you should use? Well, in the absence of a deep understanding of how the underlying algorithm works, you'll need to experiment. Fortunately, SciKit-Learn provides a way to *tune* hyperparameters by trying multiple combinations and finding the best result for a given performance metric.\n",
    "\n",
    "Let's try using a *grid search* approach to try combinations from a grid of possible values for the **learning_rate** and **n_estimators** hyperparameters of the **GradientBoostingRegressor** estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import make_scorer, r2_score\n",
    "\n",
    "# Use a Gradient Boosting algorithm\n",
    "alg = GradientBoostingRegressor()\n",
    "\n",
    "# Try these hyperparameter values\n",
    "params = {\n",
    " 'learning_rate': [0.1, 0.5, 1.0],\n",
    " 'n_estimators' : [50, 100, 150]\n",
    " }\n",
    "\n",
    "# Find the best hyperparameter combination to optimize the R2 metric\n",
    "score = make_scorer(r2_score)\n",
    "gridsearch = GridSearchCV(alg, params, scoring=score, cv=3, return_train_score=True)\n",
    "gridsearch.fit(X_train, y_train)\n",
    "print(\"Best parameter combination:\", gridsearch.best_params_, \"\\n\")\n",
    "\n",
    "# Get the best model\n",
    "model=gridsearch.best_estimator_\n",
    "print(model, \"\\n\")\n",
    "\n",
    "# Evaluate the model using the test data\n",
    "predictions = model.predict(X_test)\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "# overlay the regression line\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note**: The use of random values in the Gradient Boosting algorithm results in slightly different metrics each time. In this case, the best model produced by hyperparameter tuning is unlikely to be significantly better than one trained with the default hyperparameter values; but it's still useful to know about the hyperparameter tuning technique!\n",
    "\n",
    "## Preprocess the Data\n",
    "\n",
    "We trained a model with data that was loaded straight from a source file, with only moderately successful results.\n",
    "\n",
    "In practice, it's common to perform some preprocessing of the data to make it easier for the algorithm to fit a model to it. There's a huge range of preprocessing trasformations you can perform to get your data ready for modeling, but we'll limit ourselves to a few common techniques:\n",
    "\n",
    "- Scaling numeric features so they're on the same scale. This prevents features with large values from producing coefficients that disproportionately affect the predictions.\n",
    "- Encoding categorical variables. For example, by using a *one hot encoding* technique you can create individual binary (true/false) features for each possible category value.\n",
    "\n",
    "To apply these preprocessing transformations, we'll make use of a Scikit-Learn feature named *pipelines*. These enable us to define a set of preprocessing steps that end with an algorithm. You can then fit the entire pipeline to the data, so that the model encapsulates all of the preprocessing steps as well as the regression algorithm. This is useful, because when we want to use the model to predict values from new data, we need to apply the same transformations (based on the same statistical distributions and catagory encodings used with the training data).\n",
    "\n",
    ">**Note**: The term *pipeline* is used extensively in machine learning, often to mean very different things! In this context, we're using it to refer to pipeline objects in Scikit-Learn, but you may see it used elsewhere to mean something else."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the model\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.preprocessing import StandardScaler, OneHotEncoder\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import numpy as np\n",
    "\n",
    "# Define preprocessing for numeric columns (scale them)\n",
    "numeric_features = [6,7,8,9]\n",
    "numeric_transformer = Pipeline(steps=[\n",
    "    ('scaler', StandardScaler())])\n",
    "\n",
    "# Define preprocessing for categorical features (encode them)\n",
    "categorical_features = [0,1,2,3,4,5]\n",
    "categorical_transformer = Pipeline(steps=[\n",
    "    ('onehot', OneHotEncoder(handle_unknown='ignore'))])\n",
    "\n",
    "# Combine preprocessing steps\n",
    "preprocessor = ColumnTransformer(\n",
    "    transformers=[\n",
    "        ('num', numeric_transformer, numeric_features),\n",
    "        ('cat', categorical_transformer, categorical_features)])\n",
    "\n",
    "# Create preprocessing and training pipeline\n",
    "pipeline = Pipeline(steps=[('preprocessor', preprocessor),\n",
    "                           ('regressor', GradientBoostingRegressor())])\n",
    "\n",
    "\n",
    "# fit the pipeline to train a linear regression model on the training set\n",
    "model = pipeline.fit(X_train, (y_train))\n",
    "print (model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, the model is trained, including the preprocessing steps. Let's see how it performs with the validation data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get predictions\n",
    "predictions = model.predict(X_test)\n",
    "\n",
    "# Display metrics\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions')\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pipeline is composed of the transformations and the algorithm used to train the model. To try an alternative algorithm you can just change that step to a different kind of estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a different estimator in the pipeline\n",
    "pipeline = Pipeline(steps=[('preprocessor', preprocessor),\n",
    "                           ('regressor', RandomForestRegressor())])\n",
    "\n",
    "\n",
    "# fit the pipeline to train a linear regression model on the training set\n",
    "model = pipeline.fit(X_train, (y_train))\n",
    "print (model, \"\\n\")\n",
    "\n",
    "# Get predictions\n",
    "predictions = model.predict(X_test)\n",
    "\n",
    "# Display metrics\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print(\"MSE:\", mse)\n",
    "rmse = np.sqrt(mse)\n",
    "print(\"RMSE:\", rmse)\n",
    "r2 = r2_score(y_test, predictions)\n",
    "print(\"R2:\", r2)\n",
    "\n",
    "# Plot predicted vs actual\n",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions - Preprocessed')\n",
    "z = np.polyfit(y_test, predictions, 1)\n",
    "p = np.poly1d(z)\n",
    "plt.plot(y_test,p(y_test), color='magenta')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've now seen a number of common techniques used to train predictive models for regression. In a real project, you'd likely try a few more algorithms, hyperparameters, and preprocessing transformations; but by now you should have got the general idea. Let's explore how you can use the trained model with new data.\n",
    "\n",
    "### Use the Trained Model\n",
    "\n",
    "First, let's save the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import joblib\n",
    "\n",
    "# Save the model as a pickle file\n",
    "filename = './models/bike-share.pkl'\n",
    "joblib.dump(model, filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can load it whenever we need it, and use it to predict labels for new data. This is often called *scoring* or *inferencing*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the model from the file\n",
    "loaded_model = joblib.load(filename)\n",
    "\n",
    "# Create a numpy array containing a new observation (for example tomorrow's seasonal and weather forecast information)\n",
    "X_new = np.array([[1,1,0,3,1,1,0.226957,0.22927,0.436957,0.1869]]).astype('float64')\n",
    "print ('New sample: {}'.format(list(X_new[0])))\n",
    "\n",
    "# Use the model to predict tomorrow's rentals\n",
    "result = loaded_model.predict(X_new)\n",
    "print('Prediction: {:.0f} rentals'.format(np.round(result[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model's **predict** method accepts an array of observations, so you can use it to generate multiple predictions as a batch. For example, suppose you have a weather forecast for the next five days; you could use the model to predict bike rentals for each day based on the expected weather conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An array of features based on five-day weather forecast\n",
    "X_new = np.array([[0,1,1,0,0,1,0.344167,0.363625,0.805833,0.160446],\n",
    "                  [0,1,0,1,0,1,0.363478,0.353739,0.696087,0.248539],\n",
    "                  [0,1,0,2,0,1,0.196364,0.189405,0.437273,0.248309],\n",
    "                  [0,1,0,3,0,1,0.2,0.212122,0.590435,0.160296],\n",
    "                  [0,1,0,4,0,1,0.226957,0.22927,0.436957,0.1869]])\n",
    "\n",
    "# Use the model to predict rentals\n",
    "results = loaded_model.predict(X_new)\n",
    "print('5-day rental predictions:')\n",
    "for prediction in results:\n",
    "    print(np.round(prediction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learn More\n",
    "\n",
    "To learn more about Scikit-Learn, see the [Scikit-Learn documentation](https://scikit-learn.org/stable/user_guide.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Citation\n",
    "\n",
    "The data used in this exercise is derived from [Capital Bikeshare](https://www.capitalbikeshare.com/system-data) and is used in accordance with the published [license agreement](https://www.capitalbikeshare.com/data-license-agreement)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6 - AzureML",
   "language": "python",
   "name": "python3-azureml"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}