<jupyter_start><jupyter_text>Time-related feature engineeringThis notebook introduces different strategies to leverage time-related featuresfor a bike sharing demand regression task that is highly dependent on businesscycles (days, weeks, months) and yearly season cycles.In the process, we introduce how to perform periodic feature engineering usingthe :class:`sklearn.preprocessing.SplineTransformer` class and its`extrapolation="periodic"` option.  Data exploration on the Bike Sharing Demand datasetWe start by loading the data from the OpenML repository.<jupyter_code>from sklearn.datasets import fetch_openml

bike_sharing = fetch_openml("Bike_Sharing_Demand", version=2, as_frame=True)
df = bike_sharing.frame<jupyter_output><empty_output><jupyter_text>To get a quick understanding of the periodic patterns of the data, let ushave a look at the average demand per hour during a week.Note that the week starts on a Sunday, during the weekend. We can clearlydistinguish the commute patterns in the morning and evenings of the work daysand the leisure use of the bikes on the weekends with a more spread peakdemand around the middle of the days:<jupyter_code>import matplotlib.pyplot as plt


fig, ax = plt.subplots(figsize=(12, 4))
average_week_demand = df.groupby(["weekday", "hour"]).mean()["count"]
average_week_demand.plot(ax=ax)
_ = ax.set(
    title="Average hourly bike demand during the week",
    xticks=[i * 24 for i in range(7)],
    xticklabels=["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"],
    xlabel="Time of the week",
    ylabel="Number of bike rentals",
)<jupyter_output><empty_output><jupyter_text>The target of the prediction problem is the absolute count of bike rentals ona hourly basis:<jupyter_code>df["count"].max()<jupyter_output><empty_output><jupyter_text>Let us rescale the target variable (number of hourly bike rentals) to predicta relative demand so that the mean absolute error is more easily interpretedas a fraction of the maximum demand.NoteThe fit method of the models used in this notebook all minimize the    mean squared error to estimate the conditional mean instead of the mean    absolute error that would fit an estimator of the conditional median.    When reporting performance measure on the test set in the discussion, we    instead choose to focus on the mean absolute error that is more    intuitive than the (root) mean squared error. Note, however, that the    best models for one metric are also the best for the other in this    study.<jupyter_code>y = df["count"] / df["count"].max()
fig, ax = plt.subplots(figsize=(12, 4))
y.hist(bins=30, ax=ax)
_ = ax.set(
    xlabel="Fraction of rented fleet demand",
    ylabel="Number of hours",
)<jupyter_output><empty_output><jupyter_text>The input feature data frame is a time annotated hourly log of variablesdescribing the weather conditions. It includes both numerical and categoricalvariables. Note that the time information has already been expanded intoseveral complementary columns.<jupyter_code>X = df.drop("count", axis="columns")
X<jupyter_output><empty_output><jupyter_text>NoteIf the time information was only present as a date or datetime column, we   could have expanded it into hour-in-the-day, day-in-the-week,   day-in-the-month, month-in-the-year using pandas:   https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.htmltime-date-componentsWe now introspect the distribution of the categorical variables, startingwith `"weather"`:<jupyter_code>X["weather"].value_counts()<jupyter_output><empty_output><jupyter_text>Since there are only 3 `"heavy_rain"` events, we cannot use this category totrain machine learning models with cross validation. Instead, we simplify therepresentation by collapsing those into the `"rain"` category.<jupyter_code>X["weather"].replace(to_replace="heavy_rain", value="rain", inplace=True)
X["weather"].value_counts()<jupyter_output><empty_output><jupyter_text>As expected, the `"season"` variable is well balanced:<jupyter_code>X["season"].value_counts()<jupyter_output><empty_output><jupyter_text>Time-based cross-validationSince the dataset is a time-ordered event log (hourly demand), we will use atime-sensitive cross-validation splitter to evaluate our demand forecastingmodel as realistically as possible. We use a gap of 2 days between the trainand test side of the splits. We also limit the training set size to make theperformance of the CV folds more stable.1000 test datapoints should be enough to quantify the performance of themodel. This represents a bit less than a month and a half of contiguous testdata:<jupyter_code>from sklearn.model_selection import TimeSeriesSplit

ts_cv = TimeSeriesSplit(
    n_splits=5,
    gap=48,
    max_train_size=10000,
    test_size=1000,
)<jupyter_output><empty_output><jupyter_text>Let us manually inspect the various splits to check that the`TimeSeriesSplit` works as we expect, starting with the first split:<jupyter_code>all_splits = list(ts_cv.split(X, y))
train_0, test_0 = all_splits[0]
X.iloc[test_0]
X.iloc[train_0]<jupyter_output><empty_output><jupyter_text>We now inspect the last split:<jupyter_code>train_4, test_4 = all_splits[4]
X.iloc[test_4]
X.iloc[train_4]<jupyter_output><empty_output><jupyter_text>All is well. We are now ready to do some predictive modeling! Gradient BoostingGradient Boosting Regression with decision trees is often flexible enough toefficiently handle heteorogenous tabular data with a mix of categorical andnumerical features as long as the number of samples is large enough.Here, we do minimal ordinal encoding for the categorical variables and thenlet the model know that it should treat those as categorical variables byusing a dedicated tree splitting rule. Since we use an ordinal encoder, wepass the list of categorical values explicitly to use a logical order whenencoding the categories as integers instead of the lexicographical order.This also has the added benefit of preventing any issue with unknowncategories when using cross-validation.The numerical variables need no preprocessing and, for the sake of simplicity,we only try the default hyper-parameters for this model:<jupyter_code>from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OrdinalEncoder
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import cross_validate


categorical_columns = [
    "weather",
    "season",
    "holiday",
    "workingday",
]
categories = [
    ["clear", "misty", "rain"],
    ["spring", "summer", "fall", "winter"],
    ["False", "True"],
    ["False", "True"],
]
ordinal_encoder = OrdinalEncoder(categories=categories)


gbrt_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", ordinal_encoder, categorical_columns),
        ],
        remainder="passthrough",
    ),
    HistGradientBoostingRegressor(
        categorical_features=range(4),
    ),
)<jupyter_output><empty_output><jupyter_text>Lets evaluate our gradient boosting model with the mean absolute error of therelative demand averaged across our 5 time-based cross-validation splits:<jupyter_code>def evaluate(model, X, y, cv):
    cv_results = cross_validate(
        model,
        X,
        y,
        cv=cv,
        scoring=["neg_mean_absolute_error", "neg_root_mean_squared_error"],
    )
    mae = -cv_results["test_neg_mean_absolute_error"]
    rmse = -cv_results["test_neg_root_mean_squared_error"]
    print(
        f"Mean Absolute Error:     {mae.mean():.3f} +/- {mae.std():.3f}\n"
        f"Root Mean Squared Error: {rmse.mean():.3f} +/- {rmse.std():.3f}"
    )


evaluate(gbrt_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.044 +/- 0.003
Root Mean Squared Error: 0.068 +/- 0.005<jupyter_text>This model has an average error around 4 to 5% of the maximum demand. This isquite good for a first trial without any hyper-parameter tuning! We just hadto make the categorical variables explicit. Note that the time relatedfeatures are passed as is, i.e. without processing them. But this is not muchof a problem for tree-based models as they can learn a non-monotonicrelationship between ordinal input features and the target.This is not the case for linear regression models as we will see in thefollowing. Naive linear regressionAs usual for linear models, categorical variables need to be one-hot encoded.For consistency, we scale the numerical features to the same 0-1 range usingclass:`sklearn.preprocessing.MinMaxScaler`, although in this case it does notimpact the results much because they are already on comparable scales:<jupyter_code>from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import MinMaxScaler
from sklearn.linear_model import RidgeCV
import numpy as np


one_hot_encoder = OneHotEncoder(handle_unknown="ignore", sparse=False)
alphas = np.logspace(-6, 6, 25)
naive_linear_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
        ],
        remainder=MinMaxScaler(),
    ),
    RidgeCV(alphas=alphas),
)


evaluate(naive_linear_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.142 +/- 0.014
Root Mean Squared Error: 0.184 +/- 0.020<jupyter_text>The performance is not good: the average error is around 14% of the maximumdemand. This is more than three times higher than the average error of thegradient boosting model. We can suspect that the naive original encoding(merely min-max scaled) of the periodic time-related features might preventthe linear regression model to properly leverage the time information: linearregression does not automatically model non-monotonic relationships betweenthe input features and the target. Non-linear terms have to be engineered inthe input.For example, the raw numerical encoding of the `"hour"` feature prevents thelinear model from recognizing that an increase of hour in the morning from 6to 8 should have a strong positive impact on the number of bike rentals whilean increase of similar magnitude in the evening from 18 to 20 should have astrong negative impact on the predicted number of bike rentals. Time-steps as categoriesSince the time features are encoded in a discrete manner using integers (24unique values in the "hours" feature), we could decide to treat those ascategorical variables using a one-hot encoding and thereby ignore anyassumption implied by the ordering of the hour values.Using one-hot encoding for the time features gives the linear model a lotmore flexibility as we introduce one additional feature per discrete timelevel.<jupyter_code>one_hot_linear_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
            ("one_hot_time", one_hot_encoder, ["hour", "weekday", "month"]),
        ],
        remainder=MinMaxScaler(),
    ),
    RidgeCV(alphas=alphas),
)

evaluate(one_hot_linear_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.099 +/- 0.011
Root Mean Squared Error: 0.131 +/- 0.011<jupyter_text>The average error rate of this model is 10% which is much better than usingthe original (ordinal) encoding of the time feature, confirming our intuitionthat the linear regression model benefits from the added flexibility to nottreat time progression in a monotonic manner.However, this introduces a very large number of new features. If the time ofthe day was represented in minutes since the start of the day instead ofhours, one-hot encoding would have introduced 1440 features instead of 24.This could cause some significant overfitting. To avoid this we could use:func:`sklearn.preprocessing.KBinsDiscretizer` instead to re-bin the numberof levels of fine-grained ordinal or numerical variables while stillbenefitting from the non-monotonic expressivity advantages of one-hotencoding.Finally, we also observe that one-hot encoding completely ignores theordering of the hour levels while this could be an interesting inductive biasto preserve to some level. In the following we try to explore smooth,non-monotonic encoding that locally preserves the relative ordering of timefeatures. Trigonometric featuresAs a first attempt, we can try to encode each of those periodic featuresusing a sine and cosine transformation with the matching period.Each ordinal time feature is transformed into 2 features that together encodeequivalent information in a non-monotonic way, and more importantly withoutany jump between the first and the last value of the periodic range.<jupyter_code>from sklearn.preprocessing import FunctionTransformer


def sin_transformer(period):
    return FunctionTransformer(lambda x: np.sin(x / period * 2 * np.pi))


def cos_transformer(period):
    return FunctionTransformer(lambda x: np.cos(x / period * 2 * np.pi))<jupyter_output><empty_output><jupyter_text>Let us visualize the effect of this feature expansion on some synthetic hourdata with a bit of extrapolation beyond hour=23:<jupyter_code>import pandas as pd

hour_df = pd.DataFrame(
    np.arange(26).reshape(-1, 1),
    columns=["hour"],
)
hour_df["hour_sin"] = sin_transformer(24).fit_transform(hour_df)["hour"]
hour_df["hour_cos"] = cos_transformer(24).fit_transform(hour_df)["hour"]
hour_df.plot(x="hour")
_ = plt.title("Trigonometric encoding for the 'hour' feature")<jupyter_output><empty_output><jupyter_text>Let's use a 2D scatter plot with the hours encoded as colors to better seehow this representation maps the 24 hours of the day to a 2D space, akin tosome sort of a 24 hour version of an analog clock. Note that the "25th" houris mapped back to the 1st hour because of the periodic nature of thesine/cosine representation.<jupyter_code>fig, ax = plt.subplots(figsize=(7, 5))
sp = ax.scatter(hour_df["hour_sin"], hour_df["hour_cos"], c=hour_df["hour"])
ax.set(
    xlabel="sin(hour)",
    ylabel="cos(hour)",
)
_ = fig.colorbar(sp)<jupyter_output><empty_output><jupyter_text>We can now build a feature extraction pipeline using this strategy:<jupyter_code>cyclic_cossin_transformer = ColumnTransformer(
    transformers=[
        ("categorical", one_hot_encoder, categorical_columns),
        ("month_sin", sin_transformer(12), ["month"]),
        ("month_cos", cos_transformer(12), ["month"]),
        ("weekday_sin", sin_transformer(7), ["weekday"]),
        ("weekday_cos", cos_transformer(7), ["weekday"]),
        ("hour_sin", sin_transformer(24), ["hour"]),
        ("hour_cos", cos_transformer(24), ["hour"]),
    ],
    remainder=MinMaxScaler(),
)
cyclic_cossin_linear_pipeline = make_pipeline(
    cyclic_cossin_transformer,
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_cossin_linear_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.125 +/- 0.014
Root Mean Squared Error: 0.166 +/- 0.020<jupyter_text>The performance of our linear regression model with this simple featureengineering is a bit better than using the original ordinal time features butworse than using the one-hot encoded time features. We will further analyzepossible reasons for this disappointing outcome at the end of this notebook. Periodic spline featuresWe can try an alternative encoding of the periodic time-related featuresusing spline transformations with a large enough number of splines, and as aresult a larger number of expanded features compared to the sine/cosinetransformation:<jupyter_code>from sklearn.preprocessing import SplineTransformer


def periodic_spline_transformer(period, n_splines=None, degree=3):
    if n_splines is None:
        n_splines = period
    n_knots = n_splines + 1  # periodic and include_bias is True
    return SplineTransformer(
        degree=degree,
        n_knots=n_knots,
        knots=np.linspace(0, period, n_knots).reshape(n_knots, 1),
        extrapolation="periodic",
        include_bias=True,
    )<jupyter_output><empty_output><jupyter_text>Again, let us visualize the effect of this feature expansion on somesynthetic hour data with a bit of extrapolation beyond hour=23:<jupyter_code>hour_df = pd.DataFrame(
    np.linspace(0, 26, 1000).reshape(-1, 1),
    columns=["hour"],
)
splines = periodic_spline_transformer(24, n_splines=12).fit_transform(hour_df)
splines_df = pd.DataFrame(
    splines,
    columns=[f"spline_{i}" for i in range(splines.shape[1])],
)
pd.concat([hour_df, splines_df], axis="columns").plot(x="hour", cmap=plt.cm.tab20b)
_ = plt.title("Periodic spline-based encoding for the 'hour' feature")<jupyter_output><empty_output><jupyter_text>Thanks to the use of the `extrapolation="periodic"` parameter, we observethat the feature encoding stays smooth when extrapolating beyond midnight.We can now build a predictive pipeline using this alternative periodicfeature engineering strategy.It is possible to use fewer splines than discrete levels for those ordinalvalues. This makes spline-based encoding more efficient than one-hot encodingwhile preserving most of the expressivity:<jupyter_code>cyclic_spline_transformer = ColumnTransformer(
    transformers=[
        ("categorical", one_hot_encoder, categorical_columns),
        ("cyclic_month", periodic_spline_transformer(12, n_splines=6), ["month"]),
        ("cyclic_weekday", periodic_spline_transformer(7, n_splines=3), ["weekday"]),
        ("cyclic_hour", periodic_spline_transformer(24, n_splines=12), ["hour"]),
    ],
    remainder=MinMaxScaler(),
)
cyclic_spline_linear_pipeline = make_pipeline(
    cyclic_spline_transformer,
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_linear_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.097 +/- 0.011
Root Mean Squared Error: 0.132 +/- 0.013<jupyter_text>Spline features make it possible for the linear model to successfullyleverage the periodic time-related features and reduce the error from ~14% to~10% of the maximum demand, which is similar to what we observed with theone-hot encoded features. Qualitative analysis of the impact of features on linear model predictionsHere, we want to visualize the impact of the feature engineering choices onthe time related shape of the predictions.To do so we consider an arbitrary time-based split to compare the predictionson a range of held out data points.<jupyter_code>naive_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
naive_linear_predictions = naive_linear_pipeline.predict(X.iloc[test_0])

one_hot_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
one_hot_linear_predictions = one_hot_linear_pipeline.predict(X.iloc[test_0])

cyclic_cossin_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_cossin_linear_predictions = cyclic_cossin_linear_pipeline.predict(X.iloc[test_0])

cyclic_spline_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_spline_linear_predictions = cyclic_spline_linear_pipeline.predict(X.iloc[test_0])<jupyter_output><empty_output><jupyter_text>We visualize those predictions by zooming on the last 96 hours (4 days) ofthe test set to get some qualitative insights:<jupyter_code>last_hours = slice(-96, None)
fig, ax = plt.subplots(figsize=(12, 4))
fig.suptitle("Predictions by linear models")
ax.plot(
    y.iloc[test_0].values[last_hours],
    "x-",
    alpha=0.2,
    label="Actual demand",
    color="black",
)
ax.plot(naive_linear_predictions[last_hours], "x-", label="Ordinal time features")
ax.plot(
    cyclic_cossin_linear_predictions[last_hours],
    "x-",
    label="Trigonometric time features",
)
ax.plot(
    cyclic_spline_linear_predictions[last_hours],
    "x-",
    label="Spline-based time features",
)
ax.plot(
    one_hot_linear_predictions[last_hours],
    "x-",
    label="One-hot time features",
)
_ = ax.legend()<jupyter_output><empty_output><jupyter_text>We can draw the following conclusions from the above plot:- The **raw ordinal time-related features** are problematic because they do  not capture the natural periodicity: we observe a big jump in the  predictions at the end of each day when the hour features goes from 23 back  to 0. We can expect similar artifacts at the end of each week or each year.- As expected, the **trigonometric features** (sine and cosine) do not have  these discontinuities at midnight, but the linear regression model fails to  leverage those features to properly model intra-day variations.  Using trigonometric features for higher harmonics or additional  trigonometric features for the natural period with different phases could  potentially fix this problem.- the **periodic spline-based features** fix those two problems at once: they  give more expressivity to the linear model by making it possible to focus  on specific hours thanks to the use of 12 splines. Furthermore the  `extrapolation="periodic"` option enforces a smooth representation between  `hour=23` and `hour=0`.- The **one-hot encoded features** behave similarly to the periodic  spline-based features but are more spiky: for instance they can better  model the morning peak during the week days since this peak lasts shorter  than an hour. However, we will see in the following that what can be an  advantage for linear models is not necessarily one for more expressive  models. We can also compare the number of features extracted by each featureengineering pipeline:<jupyter_code>naive_linear_pipeline[:-1].transform(X).shape
one_hot_linear_pipeline[:-1].transform(X).shape
cyclic_cossin_linear_pipeline[:-1].transform(X).shape
cyclic_spline_linear_pipeline[:-1].transform(X).shape<jupyter_output><empty_output><jupyter_text>This confirms that the one-hot encoding and the spline encoding strategiescreate a lot more features for the time representation than the alternatives,which in turn gives the downstream linear model more flexibility (degrees offreedom) to avoid underfitting.Finally, we observe that none of the linear models can approximate the truebike rentals demand, especially for the peaks that can be very sharp at rushhours during the working days but much flatter during the week-ends: the mostaccurate linear models based on splines or one-hot encoding tend to forecastpeaks of commuting-related bike rentals even on the week-ends andunder-estimate the commuting-related events during the working days.These systematic prediction errors reveal a form of under-fitting and can beexplained by the lack of interactions terms between features, e.g."workingday" and features derived from "hours". This issue will be addressedin the following section.  Modeling pairwise interactions with splines and polynomial featuresLinear models do not automatically capture interaction effects between inputfeatures. It does not help that some features are marginally non-linear as isthe case with features constructed by `SplineTransformer` (or one-hotencoding or binning).However, it is possible to use the `PolynomialFeatures` class on coarsegrained spline encoded hours to model the "workingday"/"hours" interactionexplicitly without introducing too many new variables:<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import FeatureUnion


hour_workday_interaction = make_pipeline(
    ColumnTransformer(
        [
            ("cyclic_hour", periodic_spline_transformer(24, n_splines=8), ["hour"]),
            ("workingday", FunctionTransformer(lambda x: x == "True"), ["workingday"]),
        ]
    ),
    PolynomialFeatures(degree=2, interaction_only=True, include_bias=False),
)<jupyter_output><empty_output><jupyter_text>Those features are then combined with the ones already computed in theprevious spline-base pipeline. We can observe a nice performance improvemntby modeling this pairwise interaction explicitly:<jupyter_code>cyclic_spline_interactions_pipeline = make_pipeline(
    FeatureUnion(
        [
            ("marginal", cyclic_spline_transformer),
            ("interactions", hour_workday_interaction),
        ]
    ),
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_interactions_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.078 +/- 0.009
Root Mean Squared Error: 0.104 +/- 0.009<jupyter_text>Modeling non-linear feature interactions with kernelsThe previous analysis highlighted the need to model the interactions between`"workingday"` and `"hours"`. Another example of a such a non-linearinteraction that we would like to model could be the impact of the rain thatmight not be the same during the working days and the week-ends and holidaysfor instance.To model all such interactions, we could either use a polynomial expansion onall marginal features at once, after their spline-based expansion. However,this would create a quadratic number of features which can cause overfittingand computational tractability issues.Alternatively, we can use the Nyström method to compute an approximatepolynomial kernel expansion. Let us try the latter:<jupyter_code>from sklearn.kernel_approximation import Nystroem


cyclic_spline_poly_pipeline = make_pipeline(
    cyclic_spline_transformer,
    Nystroem(kernel="poly", degree=2, n_components=300, random_state=0),
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_poly_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.053 +/- 0.002
Root Mean Squared Error: 0.076 +/- 0.004<jupyter_text>We observe that this model can almost rival the performance of the gradientboosted trees with an average error around 5% of the maximum demand.Note that while the final step of this pipeline is a linear regression model,the intermediate steps such as the spline feature extraction and the Nyströmkernel approximation are highly non-linear. As a result the compound pipelineis much more expressive than a simple linear regression model with raw features.For the sake of completeness, we also evaluate the combination of one-hotencoding and kernel approximation:<jupyter_code>one_hot_poly_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
            ("one_hot_time", one_hot_encoder, ["hour", "weekday", "month"]),
        ],
        remainder="passthrough",
    ),
    Nystroem(kernel="poly", degree=2, n_components=300, random_state=0),
    RidgeCV(alphas=alphas),
)
evaluate(one_hot_poly_pipeline, X, y, cv=ts_cv)<jupyter_output>Mean Absolute Error:     0.082 +/- 0.006
Root Mean Squared Error: 0.111 +/- 0.011<jupyter_text>While one-hot encoded features were competitive with spline-based featureswhen using linear models, this is no longer the case when using a low-rankapproximation of a non-linear kernel: this can be explained by the fact thatspline features are smoother and allow the kernel approximation to find amore expressive decision function.Let us now have a qualitative look at the predictions of the kernel modelsand of the gradient boosted trees that should be able to better modelnon-linear interactions between features:<jupyter_code>gbrt_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
gbrt_predictions = gbrt_pipeline.predict(X.iloc[test_0])

one_hot_poly_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
one_hot_poly_predictions = one_hot_poly_pipeline.predict(X.iloc[test_0])

cyclic_spline_poly_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_spline_poly_predictions = cyclic_spline_poly_pipeline.predict(X.iloc[test_0])<jupyter_output><empty_output><jupyter_text>Again we zoom on the last 4 days of the test set:<jupyter_code>last_hours = slice(-96, None)
fig, ax = plt.subplots(figsize=(12, 4))
fig.suptitle("Predictions by non-linear regression models")
ax.plot(
    y.iloc[test_0].values[last_hours],
    "x-",
    alpha=0.2,
    label="Actual demand",
    color="black",
)
ax.plot(
    gbrt_predictions[last_hours],
    "x-",
    label="Gradient Boosted Trees",
)
ax.plot(
    one_hot_poly_predictions[last_hours],
    "x-",
    label="One-hot + polynomial kernel",
)
ax.plot(
    cyclic_spline_poly_predictions[last_hours],
    "x-",
    label="Splines + polynomial kernel",
)
_ = ax.legend()<jupyter_output><empty_output><jupyter_text>First, note that trees can naturally model non-linear feature interactionssince, by default, decision trees are allowed to grow beyond a depth of 2levels.Here, we can observe that the combinations of spline features and non-linearkernels works quite well and can almost rival the accuracy of the gradientboosting regression trees.On the contrary, one-hot encoded time features do not perform that well withthe low rank kernel model. In particular, they significantly over-estimatethe low demand hours more than the competing models.We also observe that none of the models can successfully predict some of thepeak rentals at the rush hours during the working days. It is possible thataccess to additional features would be required to further improve theaccuracy of the predictions. For instance, it could be useful to have accessto the geographical repartition of the fleet at any point in time or thefraction of bikes that are immobilized because they need servicing.Let us finally get a more quantative look at the prediction errors of thosethree models using the true vs predicted demand scatter plots:<jupyter_code>fig, axes = plt.subplots(ncols=3, figsize=(12, 4), sharey=True)
fig.suptitle("Non-linear regression models")
predictions = [
    one_hot_poly_predictions,
    cyclic_spline_poly_predictions,
    gbrt_predictions,
]
labels = [
    "One hot + polynomial kernel",
    "Splines + polynomial kernel",
    "Gradient Boosted Trees",
]
for ax, pred, label in zip(axes, predictions, labels):
    ax.scatter(y.iloc[test_0].values, pred, alpha=0.3, label=label)
    ax.plot([0, 1], [0, 1], "--", label="Perfect model")
    ax.set(
        xlim=(0, 1),
        ylim=(0, 1),
        xlabel="True demand",
        ylabel="Predicted demand",
    )
    ax.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0777158587350611, 'fraction_numerical': 0.011850409560897005, 'mean_word_length': 5.289907094594595, 'pattern_counts': {'":': 0, '<': 121, '<?xml version=': 0, '>': 121, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '65021', 'license': 'MIT', 'n_tokens_mistral': 8842, 'n_tokens_neox': 8315, 'n_words': 3105, 'path': 'testpad/plot_cyclical_feature_engineering.ipynb', 'repo_name': 'furyhawk/kaggle_practice', 'size': '29788'}