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 **scikit-learn** implementation of the *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."
   ]
  },
  {
   "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(normalize=False).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 difference between predicted and actual, squared. 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!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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 feaures 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', LinearRegression(normalize=False))])\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": [
    "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])\n",
    "\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",
    "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": [
    "It's a little better - great!\n",
    "\n",
    "### Try a Different Algorithm\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",
    "Let's try training our regression model by using an ensemble algorithm called **Gradient Boosting**. We can do this by just changing the last step of the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "\n",
    "pipeline = Pipeline(steps=[('preprocessor', preprocessor),\n",
    "                            ('regressor', GradientBoostingRegressor())])\n",
    "\n",
    "\n",
    "# train a logistic regression model on the training set\n",
    "model = pipeline.fit(X_train, (y_train))\n",
    "print (model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, let's see how our ensemble model performs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = model.predict(X_test)\n",
    "\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",
    "plt.scatter(y_test, predictions)\n",
    "plt.xlabel('Actual Labels')\n",
    "plt.ylabel('Predicted Labels')\n",
    "plt.title('Daily Bike Share Predictions - Ensemble')\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 certainly some improvement!\n",
    "\n",
    "In a real project, you'd likely try a few more algorithms along with more preprocessing; 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": [
    "## 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",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}