Bike Sharing Demand Problem, part 1 – Sklearn regression model, scaling, transformation chain

import warnings
warnings.filterwarnings("ignore")

Libraries

from sklearn import model_selection, linear_model, metrics
import numpy as np
import pandas as pd

%pylab inline

Dataset load

We load data into dataframe.

raw_data = pd.read_csv("bike_sharing_demand.csv", header = 0, sep = ",")

raw_data.head()

datetime – hourly date + timestamp

season – 1 = spring, 2 = summer, 3 = fall, 4 = winter

holiday – whether the day is considered a holiday

workingday – whether the day is neither a weekend nor holiday

weather – 1: Clear, Few clouds, Partly cloudy, Partly cloudy 2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist 3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds 4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog

temp – temperature in Celsius

atemp – “feels like” temperature in Celsius

humidity – relative humidity

windspeed – wind speed

casual – number of non-registered user rentals initiated

registered – number of registered user rentals initiated

count – number of total bike rentals, the target label

print(raw_data.shape)

(10886, 12)

Let us check is any values that are missing in dataset. If so, we need to process to fill them in.

raw_data.isnull().values.any()

The negative result shows there are no missed data cells/ values.

Data preprocessing

Features’ types

raw_data.info()

<class "pandas.core.frame.DataFrame">
RangeIndex: 10886 entries, 0 to 10885
Data columns (total 12 columns):
 #   Column      Non-Null Count  Dtype  
---  ------      --------------  -----  
 0   datetime    10886 non-null  object 
 1   season      10886 non-null  int64  
 2   holiday     10886 non-null  int64  
 3   workingday  10886 non-null  int64  
 4   weather     10886 non-null  int64  
 5   temp        10886 non-null  float64
 6   atemp       10886 non-null  float64
 7   humidity    10886 non-null  int64  
 8   windspeed   10886 non-null  float64
 9   casual      10886 non-null  int64  
 10  registered  10886 non-null  int64  
 11  count       10886 non-null  int64  
dtypes: float64(3), int64(8), object(1)
memory usage: 1020.7+ KB

We may see here that almost all the data are presented as numbers, whether integers or floating point ones.

Transform into the datetime type

So, let’s transform datetime feature (of a type object) into the datetime type in the dataframe:

raw_data.datetime = raw_data.datetime.apply(pd.to_datetime)
raw_data.sample(3)

2 new features from datetime

Based on that transformed feature we calculate 2 new features: month and hour of the event (rental).

raw_data["month"] = raw_data.datetime.apply(lambda x : x.month)
raw_data["hour"] = raw_data.datetime.apply(lambda x : x.hour)

raw_data.sample(3)

Splitting the model into the train and test sets.

Note that the data are time dependent, time spread. So we might get the train data of earlier period and evaluate the built estimator on the test set of later period. Since the data are sorted by time, we just cut of the last 1000 rows and use them as a test set.

train_data = raw_data.iloc[:-1000, :]
hold_out_test_data = raw_data.iloc[-1000:, :]

print(raw_data.shape, train_data.shape, hold_out_test_data.shape)

(10886, 14) (9886, 14) (1000, 14)

print("train period from {} to {}".format(train_data.datetime.min(), train_data.datetime.max()))
print("evaluation period from {} to {}".format(hold_out_test_data.datetime.min(), hold_out_test_data.datetime.max()))

train period from 2011-01-01 00:00:00 to 2012-10-16 06:00:00
evaluation period from 2012-10-16 07:00:00 to 2012-12-19 23:00:00

Data and target function

We cut off the target, column count from the rest of dataset.

We also cut off the datetime data column since it’s now only an object identifier.

# training dataset
train_labels = train_data["count"].values
train_data = train_data.drop(["datetime", "count"], axis = 1)

# test dataset
test_labels = hold_out_test_data["count"].values
test_data = hold_out_test_data.drop(["datetime", "count"], axis = 1)

Let’s visualize the target function (label) at the train and test sets.

Remember that target label (count) shows how many bikes were rented at that time spot (hour).

pylab.figure(figsize = (16, 6))
pylab.subplot(1,2,1)
pylab.hist(train_labels)
pylab.title("train data")
pylab.subplot(1,2,2)
pylab.hist(test_labels)
pylab.title("test data")

Text(0.5, 1.0, "test data")

Numeric features

In this post we work only with numeric features. Let’s separate them.

numeric_columns = ["temp", "atemp", "humidity",\
"windspeed", "casual", "registered", "month", "hour"]

train_data = train_data[numeric_columns]
test_data = test_data[numeric_columns]

train_data.head()

We now train a model with SGDRegressor

SGDRegressor develops a model based on Stochastic Gradient Descent.

We supply the default paramets for the regressor model.

regressor = linear_model.SGDRegressor(random_state = 0, max_iter=5)
regressor.get_params()

{"alpha": 0.0001,
 "average": False,
 "early_stopping": False,
 "epsilon": 0.1,
 "eta0": 0.01,
 "fit_intercept": True,
 "l1_ratio": 0.15,
 "learning_rate": "invscaling",
 "loss": "squared_loss",
 "max_iter": 5,
 "n_iter_no_change": 5,
 "penalty": "l2",
 "power_t": 0.25,
 "random_state": 0,
 "shuffle": True,
 "tol": 0.001,
 "validation_fraction": 0.1,
 "verbose": 0,
 "warm_start": False}

1. Train the model regressor.fit()
2. Estimate the model with MAE mean_absolute_error

regressor.fit(train_data, train_labels)
metrics.mean_absolute_error(test_labels, regressor.predict(test_data))

9340385490952.236

The number is too big…

Let’s output the targets and model predictions:

print("Target\n", test_labels[:10])

Target
 [525 835 355 222 228 325 328 308 346 446]

print("Prediction\n", regressor.predict(test_data)[:10])

Prediction
 [1.47320172e+13 2.64787988e+13 1.23495857e+13 1.35429951e+13
 1.25998915e+13 1.74809610e+13 1.75786549e+13 2.11426913e+13
 1.96417926e+13 2.32672606e+13]

The model predicts super big numbers. It should not be this way.

Let’s take a look at the regression coefficients.

regressor.coef_

array([ 1.12866511e+10,  6.68434729e+10, -7.05920006e+10, -2.09051100e+10,
        1.65497847e+11,  2.92967602e+10, -6.71456293e+10,  8.32559786e+10])

The regression model coefficients are too high. This might happen if the features differ in scale.

Scaling the features

from sklearn.preprocessing import StandardScaler

We separately apply StandardScaler to the train and test data.

# create StandardScaler and train/fit it
scaler = StandardScaler()
scaler.fit(train_data, train_labels)
scaled_train_data = scaler.transform(train_data)
scaled_test_data = scaler.transform(test_data)

After we’ve got scaled data we train the model again.

regressor.fit(scaled_train_data, train_labels)
metrics.mean_absolute_error(test_labels, regressor.predict(scaled_test_data))

0.11793762963431514

The error is in this case is very low.

If we output targets and predictions (below), we see that the model has done a good estimation. The distance is less than 1 bike.

print(test_labels[:10])

[525 835 355 222 228 325 328 308 346 446]

list(map(lambda x : round(x, 2), regressor.predict(scaled_test_data)[:10]))

[524.83, 834.83, 354.86, 221.89, 227.84, 324.88, 327.9, 307.93, 345.9, 445.9]

The regression result is suspiciously good

We take a look at the coefficients and see that almost all of the weights are too small.

#print(regressor.coef_)
map_obj = map(lambda x : round(x, 2), regressor.coef_)
print("regressor.coef_\n", list(map_obj))

regressor.coef_
 [1.41, -1.41, 0.01, -0.04, 50.86, 148.0, -0.01, 0.01]

What are the features with non-small weights, 50.86, 148.0 ?

Casual and register features are.

How does the bike rental system work?

The bike rental system works as follows: both registered and unregistered users can use the system. In this case, column registered represent the number of registered users who use the system. The casual column shows the number of users who have not registered but also want to rent a bike.

Let’s output the values of the target function and see the following pattern: if the value of these two columns (registered and casual) is added together, we get our target label.

train_labels[:10]

array([16, 40, 32, 13,  1,  1,  2,  3,  8, 14], dtype=int64)

np.all(train_data.registered + train_data.casual == train_labels)

This was our mistake. We indeed used the features in dataset (casual, registered) that replicate the target, count.

So we cut off these data (columns) from train and test sets to properly build regession model.

train_data.drop(["casual", "registered"], axis = 1, inplace = True)
test_data.drop(["casual", "registered"], axis = 1, inplace = True)