Bike Sharing Prediction Project Lukas Tapparello

In this notebook I will build a regression model to predict the number of bikes rented each day using data from the Bike Sharing Dataset.

I'll start by loading the dataset and doing some basic exploration to understand what type of data im working with

Project Objective

The aim of this project is to predict the number of bike rentals based on various factors such as weather, calendar variables, and time.

By accurately forecasting demand, bike sharing companies can optimize bike distribution, plan staffing, and improve customer satisfaction.

This project explores and compares different machine learning model including Linear Regression, Decision tree, random forest, and a Neural Network to find the best performing approach.

In [1]:
# First, I'm going to import the main libraries I will use in the project

import pandas as pd  

import matplotlib.pyplot as plt 

import seaborn as sns 

# This will make plots show directly in notebook

%matplotlib inline

# Now I will load the dataset which contains daily data of bike rent

df = pd.read_csv('day.csv')

#I will checkthe first 5 rows of the dataset

df.head()
Out[1]:
instant dteday season yr mnth holiday weekday workingday weathersit temp atemp hum windspeed casual registered cnt
0 1 2011-01-01 1 0 1 0 6 0 2 0.344167 0.363625 0.805833 0.160446 331 654 985
1 2 2011-01-02 1 0 1 0 0 0 2 0.363478 0.353739 0.696087 0.248539 131 670 801
2 3 2011-01-03 1 0 1 0 1 1 1 0.196364 0.189405 0.437273 0.248309 120 1229 1349
3 4 2011-01-04 1 0 1 0 2 1 1 0.200000 0.212122 0.590435 0.160296 108 1454 1562
4 5 2011-01-05 1 0 1 0 3 1 1 0.226957 0.229270 0.436957 0.186900 82 1518 1600

Now I’ll see how big the dataset is, what data types the columns have, and if there is any missing data.

In [2]:
# check how many rows and columns the dataset has
df.shape
Out[2]:
(731, 16)
In [3]:
# now let's see the types of each column
df.dtypes
Out[3]:
instant         int64
dteday            str
season          int64
yr              int64
mnth            int64
holiday         int64
weekday         int64
workingday      int64
weathersit      int64
temp          float64
atemp         float64
hum           float64
windspeed     float64
casual          int64
registered      int64
cnt             int64
dtype: object
In [4]:
# last quick thing, checking if there's any missing values

df.isnull().sum()
Out[4]:
instant       0
dteday        0
season        0
yr            0
mnth          0
holiday       0
weekday       0
workingday    0
weathersit    0
temp          0
atemp         0
hum           0
windspeed     0
casual        0
registered    0
cnt           0
dtype: int64

Next thing I will check some basic statistics to see the distribution of values in the dataset

In [5]:
# this gives me a quick summary
df.describe()
Out[5]:
instant season yr mnth holiday weekday workingday weathersit temp atemp hum windspeed casual registered cnt
count 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000 731.000000
mean 366.000000 2.496580 0.500684 6.519836 0.028728 2.997264 0.683995 1.395349 0.495385 0.474354 0.627894 0.190486 848.176471 3656.172367 4504.348837
std 211.165812 1.110807 0.500342 3.451913 0.167155 2.004787 0.465233 0.544894 0.183051 0.162961 0.142429 0.077498 686.622488 1560.256377 1937.211452
min 1.000000 1.000000 0.000000 1.000000 0.000000 0.000000 0.000000 1.000000 0.059130 0.079070 0.000000 0.022392 2.000000 20.000000 22.000000
25% 183.500000 2.000000 0.000000 4.000000 0.000000 1.000000 0.000000 1.000000 0.337083 0.337842 0.520000 0.134950 315.500000 2497.000000 3152.000000
50% 366.000000 3.000000 1.000000 7.000000 0.000000 3.000000 1.000000 1.000000 0.498333 0.486733 0.626667 0.180975 713.000000 3662.000000 4548.000000
75% 548.500000 3.000000 1.000000 10.000000 0.000000 5.000000 1.000000 2.000000 0.655417 0.608602 0.730209 0.233214 1096.000000 4776.500000 5956.000000
max 731.000000 4.000000 1.000000 12.000000 1.000000 6.000000 1.000000 3.000000 0.861667 0.840896 0.972500 0.507463 3410.000000 6946.000000 8714.000000

Now I will check if people rent more bikes in certain seasons.

I will group the data by season and look at how the number of rentals change

In [6]:
# plotting average number of rentals for each season

sns.barplot(x='season', y='cnt', data=df)

# add title labels

plt.title('Average Bike Rentals per Season')
plt.xlabel('Season')
plt.ylabel('Total Rentals')

plt.show()
No description has been provided for this image

From this plot, I can see that the highest number of bike rentals happens during season 3 (fall), followed by season 2 (summer).

Season 1 (spring) has the lowest average rentals. This suggests that weather and seasonal conditions affect bike usage quite a lot.

Here I want to check how bike rentals change in the year

I will group the data by month and look at the average number of rents

In [7]:
#showing how average rentals change month by month

sns.barplot(x='mnth', y='cnt', data=df)

plt.title('Average Bike Rentals per Month')  # title
plt.xlabel('Month')  # x axis is months
plt.ylabel('Total Rentals')  # y axis is total bike rentals

plt.show()
No description has been provided for this image

Rentals increase steadily from January and peak around June–September, then drop again towards December.

This confirms that warmer months have more activity, which makes sense as people prefer biking when the weather is nicer

Bike rentals on working days vs non-working days

Now I will check if people rent bikes more during the week or on weekends and holidays.

The workingday column tells me that:
1 means it’s a working day
0 means it’s a weekend or holiday

In [8]:
# barplot to compare working vs non-working days
sns.barplot(x='workingday', y='cnt', data=df)

plt.title('Average Bike Rentals Working Day vs Non Working Day')  # adding a title
plt.xlabel('Working Day (1 = Yes, 0 = No)')  # clarifying what the numbers mean
plt.ylabel('Total Rentals')  # y-axis as usual
plt.show()
No description has been provided for this image

The average number of rentals is slightly higher on working days than on nonn working days, but the difference is not very big.

This might be because people also use bikes for leisure during weekend

Now I will look at the correlations between the variables.

This helps me see which features are more strongly connected to the target variable "cnt".

In [9]:
# calculating correlation matrix but dropping the date column ( not numeric)
corr_matrix = df.drop('dteday', axis=1).corr()

# now plotting the heatmap
plt.figure(figsize=(12, 8)) #size
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap='coolwarm') #show values inside
plt.title('Correlation Matrix')

plt.show()
No description has been provided for this image

The strongest positive correlation with the target "cnt" is with "registered" (0.95) and "casual" (0.67), which makes sense since both make up the total count

Also "temp" and "atemp" have a good correlation with "cnt" meaning weather is a strong factor

Features like "weekday" or "holiday" have very weak or no correlation

Distribution of total bike rentals

Now I want to look at how the total number of rentals is distributed

This is useful to see if most days have similar demand or if there are many very high low values

In [10]:
                                   # histogram of total rentals
plt.figure(figsize=(8, 5) )

sns.histplot(df['cnt'], bins=30, kde=True)

plt.title('Distribution of Total Bike Rentals (cnt)')

plt.xlabel('Number  of Rents')

plt.ylabel(' Frequency ')

plt.show()
No description has been provided for this image
In [11]:
# Plotting trend of bike rentals over time

plt.figure(figsize=(10, 4))

plt.plot (df['cnt'].values)
plt.title('Bike Rentals Over Time')

plt.xlabel( 'Days')
plt.ylabel('Total Rentals')

plt.show()
No description has been provided for this image

The total number of bike rentals seems to follow a roughly normal distribution, but slightly skewed to the right.

Most days have between 3000 and 5000 rentals, with a few days going as high as 8000

This tells me that: The data is fairly balanced, no extreme outliers or massive gap is seen

Predicting "cnt" should be easier for models because there’s a clear center

I dont need to apply any transformation (like log scale) since the distribution is already clean

Relationship between temperature and bike rentals:

Now I will create a scatter plot between temperature and total bike rentals (cnt).

My goal here is to check if there's a clear trend for example "do more people rent bikes when it’s warmer? "

Since temperature is already normalized in the dataset (between 0 and 1) we can use it directly.

In [12]:
plt.figure(figsize=(8, 5))
sns.scatterplot(x='temp', y='cnt', data=df)
plt.title ('Bike Rentals vs Temperature')
plt.xlabel ('Normalized Temperature')
plt.ylabel ('Total Rentals')


plt.show()

# scatter plot: cnt vs temp
No description has been provided for this image

This scatter plot shows a clear positive trend between temperature and the number of bike rentals

When the temperature is low, people rent fewer bikes. As temperature increases around 0.6 to 0.7 on the normalized scale, rentals go up

However there’s a slight drop at the highest temperatures because extremely hot days discourage biking

So in general warmer weather seems to increase bike usage but there's a point where it stops being helpful.

Rentals by day of the week

Now I want to see if there's a difference in how many bikes are rented depending on the day of the week.

The "weekday" column goes from 0 (Sunday) to 6 (Saturday)

Using a boxplot will help me spot patterns and see the spread of rentals on each day.

In [13]:
# boxplot: cnt by weekday

plt.figure(figsize=(8, 5))
sns.boxplot( x='weekday',  y='cnt', data=df )
plt.title( 'Bike Rentals by Weekday' )
plt.xlabel('Weekday (0 = Sunday)')
plt.ylabel ('Total Rentals')


plt.show()
No description has been provided for this image

What I observe from this plot?

The number of bike rentals is fairly consistent across all days of the week

There is a slight increase ondays 5 and 6 (which are Friday and Saturday), suggesting higher bike usage during the weekend.

The spread (difference between high and low usage) is quite wide every day so the number of rentals can vary a lot but the average stays prettystable.

This tell me that weekday is not a very strong predictor but weekend behavior might still have a small effect

What about weather and rents ?

Now i want to understand how weather affects the number of bike rentals, this is useful to see if people tend to rent less bikes when the weather gets worse.

In [14]:
weather_avg = df.groupby('weathersit')['cnt'].mean( )    # average rentals for each weather situation

plt.figure(figsize=(8, 5))
plt.bar(weather_avg.index, weather_avg.values, color='skyblue')

plt.title ('Average Bike Rentals by Weather Situation')
plt.xlabel ('Weather Situation')
plt.ylabel ('Average Rentals')

plt.xticks(ticks=[1, 2, 3, 4 ], labels=['Clear', 'Mist Cloudy', 'Light Rain Snow', 'Heavy Rain'], rotation=10)
plt.grid(True)

plt.show()
No description has been provided for this image

The bar chart shows a clear relationship between the weather and the number of bike rentals

The best weather condition (clear ) has the highest average rentals.

As the weather worsens from cloudy to light rain or snow the average rentals drop a lot.

When the weather is very bad rentals are almost zero.

This shows how sensitive bike demand is to the weather conditions and why its an important factor for prediction.

Here I create some new features that could help the model learn better patterns:

-is_weekend : marks if the day is Saturday or Sunday

-cnt_lag1: the number of rentals from the previous day so the model can see short term trends

-temp_diff: difference between temperature and “feels like” temperature

These engineered features will be added to the dataset and used in the prediction models

In [15]:
# Is the day a weekend?

df['is_weekend'] = df['weekday'].apply(lambda x: 1 if x == 0 or x == 6 else 0)

# Lag feature: bike count of the previous day

df['cnt_lag1'] = df['cnt'].shift(1)

# Temperature change from previous day

df['temp_diff'] = df['temp'].diff()

# Drop the first row since it has NaN from shift and diff

df = df.dropna().reset_index(drop=True)
In [16]:
X = df[['season', 'yr', 'mnth', 'holiday', 'weekday', 'workingday',
        'weathersit', 'temp', 'atemp', 'hum', 'windspeed',
        'is_weekend', 'cnt_lag1', 'temp_diff']]

Now I will scale the feature data using StandardScaler

This is useful because some features like humidity and temperature are on different scales

Now we will be preparing the data before scaling

Earlier I tried to scale the feature data using "StandardScaler", but I got an error because the variable "X" was not defined

This happened because I had not yet separated the features and the target variable.

So now, I will:

Remove the columns that are not useful for the predictio

Define X as the input features

Define y as the target variable (cnt)

In [17]:
# Is the day a weekend?

df['is_weekend'] = df['weekday'].apply(lambda x: 1 if x == 0 or x == 6 else 0)

# Lag feature: bike count of the previous day

df['cnt_lag1'] = df['cnt'].shift(1)

# Temperature change from previous day

df['temp_diff'] = df['temp'].diff()

# Drop the first row since it has NaN from shift and diff

df = df.dropna().reset_index(drop=True)

# Define the input features

X = df[['season', 'yr', 'mnth', 'holiday', 'weekday', 'workingday',
        'weathersit', 'temp', 'atemp', 'hum', 'windspeed',
        'is_weekend', 'cnt_lag1', 'temp_diff']]

# Define the target

y = df['cnt']

Now that X is defined, I will scale it using StandardScaler.

I am doing this because many machine learning models perform better when the input features are on the same scale

For example, temperature values are between 0 and 1, but humidity might be very different

Scaling makes sure no variable dominates just because of its range

In [18]:
from sklearn.preprocessing import StandardScaler

# creating the scaler
scaler = StandardScaler()

# applying it to X (the features)

X_scaled = scaler.fit_transform(X)  # now it works because X is defined

Training a Linear Regression model

Now I will train a basic Linear Regression model using the scaled feature data.

I will first split the data into training and test sets, so I can check how well the model performs Then I will fit the model and evaluate it using standard metrics for regression tasks.

In [19]:
#splitting the data

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# split the updated feature set


X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)

# create and train the mode

lr_model = LinearRegression()
lr_model.fit(X_train, y_train)

# make predictions

lr_preds = lr_model.predict(X_test)

# evaluate performance

lr_r2 = r2_score(y_test, lr_preds)
lr_rmse = mean_squared_error(y_test, lr_preds) ** 0.5

# print results

print("Linear Regression R2:", round(lr_r2, 3))

print("Linear Regression RMSE:", round(lr_rmse, 2))
Linear Regression R2: 0.87
Linear Regression RMSE: 670.61

The model gives me two metrics:

R2 score

This shows how well the model explains the variation in the target variable (cnt)

A value close to 1 means the model explains most of the variation in bike rentals A value closer to 0 means the model does not explain much and performs poorly

RMSE (Root Mean Squared Error)

This tells me how far off the models predictions are from the actual values on average

The value is in the same unit as the target (cnt) so I can interpret it as the average error in number of rented bikes

Now I will try a different model (Decision tree) and compare the results.

Now I will try a Decision Tree model

Decision trees work by splitting the data into smaller and smaller groups based on the features trying to reduce prediction error at each step

This model usually fits the training data well but it can sometimes overfit

Later I will compare its performance with the Linear Regression model to see which one works better

In [20]:
from sklearn.tree import DecisionTreeRegressor

tree = DecisionTreeRegressor(random_state=42)  # create and train the decision tree model
tree.fit(X_train, y_train)

#make predictions
y_pred_tree = tree.predict(X_test)

#evaluate the model
r2_tree = r2_score(y_test, y_pred_tree)
rmse_tree = mean_squared_error(y_test, y_pred_tree) ** 0.5

#print results

print("Decision Tree R2:", round(r2_tree, 3))

print("Decision Tree RMSE:", round(rmse_tree, 2))
Decision Tree R2: 0.767
Decision Tree RMSE: 898.44

Here I used a Decision Tree Regressor to predict bike rentals

The R*2 score tells me how well the model explains the variation in the target variable

The RMSE shows the average error in the predicted number of rentals

If the R2 is higher and the RMSE is lower than the values from the Linear Regression model
then the Decision Tree is performing better

If not, it may be overfitting or not generalizing well

In the next step, I will try a Random forest model and compare all results

Now I will train a Random forest model.

This model combines the results of several decision trees and takes the average to make the final prediction

Random Forest usually gives better performance than a single decision tree because it reduces overfitting and makes the model more stable

In [21]:
from sklearn.ensemble import RandomForestRegressor


rf = RandomForestRegressor(random_state=42) #create and train the random forest model
rf.fit(X_train, y_train)

# make predictions

y_pred_rf = rf.predict(X_test)

# evaluate the model

r2_rf = r2_score(y_test, y_pred_rf)
rmse_rf = mean_squared_error(y_test, y_pred_rf) ** 0.5

# print results

print ("Random Forest R2:", round(r2_rf, 3))

print("Random Forest RMSE:", round(rmse_rf, 2))
Random Forest R2: 0.854
Random Forest RMSE: 711.9

Gradient Boosting Regressor

This model builds trees sequentially, improving each new tree based on the errors of the previous ones.
It usually performs very well for structured data like this.

In [22]:
from sklearn.ensemble import GradientBoostingRegressor

# Create and train the model

gb = GradientBoostingRegressor(random_state=42)
gb.fit(X_train, y_train)

# Make predictions

y_pred_gb = gb.predict(X_test)

r2_gb = r2_score(y_test, y_pred_gb)
rmse_gb = mean_squared_error(y_test, y_pred_gb) ** 0.5


print ("Gradient Boosting R2:", round(r2_gb, 3))

print ("Gradient Boosting RMSE:", round(rmse_gb, 2))
Gradient Boosting R2: 0.868
Gradient Boosting RMSE: 678.0

cross-validation of models

To make our evaluation more robust, we now apply cross validation.

This splits the training set into 5 parts, trains on 4 and tests on 1, repeating the process.

It gives a better estimate of model generalization compared to a single train-test split.

In [23]:
!pip install xgboost
error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install
    python3-xyz, where xyz is the package you are trying to
    install.
    
    If you wish to install a non-Debian-packaged Python package,
    create a virtual environment using python3 -m venv path/to/venv.
    Then use path/to/venv/bin/python and path/to/venv/bin/pip. Make
    sure you have python3-full installed.
    
    If you wish to install a non-Debian packaged Python application,
    it may be easiest to use pipx install xyz, which will manage a
    virtual environment for you. Make sure you have pipx installed.
    
    See /usr/share/doc/python3.12/README.venv for more information.

note: If you believe this is a mistake, please contact your Python installation or OS distribution provider. You can override this, at the risk of breaking your Python installation or OS, by passing --break-system-packages.
hint: See PEP 668 for the detailed specification.
In [24]:
from xgboost import XGBRegressor

from sklearn.model_selection import cross_val_score

# Define models

models = {
    "Linear Regression": LinearRegression(),
    "Decision Tree": DecisionTreeRegressor(random_state=42),
    "Random Forest": RandomForestRegressor(random_state=42),
    "XGBoost": XGBRegressor(random_state=42, verbosity=0)
}

# apply 5 fold cross validation to each model

cv_results = {}

for name, model in models.items():
    scores = cross_val_score(model, X_train, y_train, cv=5, scoring='neg_mean_squared_error')
    rmse_scores = (-scores) ** 0.5  # Convert negative MSE to RMSE
    cv_results[name] = rmse_scores
    
    print(f"{name} CV RMSE: Mean = {rmse_scores.mean():.2f}, Std = {rmse_scores.std():.2f}")
Linear Regression CV RMSE: Mean = 821.34, Std = 43.99
Decision Tree CV RMSE: Mean = 1108.39, Std = 121.36
Random Forest CV RMSE: Mean = 764.15, Std = 58.86
XGBoost CV RMSE: Mean = 757.74, Std = 26.91

Cross-Validation of All Models

Here in the code cell above I just applied 5 fold cross validation to all the models I trained (Linear Regression, Decision Tree, Random Forest, and XGBoost).

This means the training data is split into 5 parts: each model is trained on 4 parts and validated on the remaining one, rotating 5 times.

This gives a more reliable estimate of each model performance and reduces the risk of overfitting compared to a single train-test split.

I used the Root Mean Squared Error (RMSE ) as the evaluation metric since it’s more interpretable in the context of regression problems like this one.

Hyperparameter Tuning for Random Forest

Here I will use GridSearchcv to find the best hyperparameters for the Random Forest model.

It tests multiple combinations of parameters using 5 fold cross validation and selects the one with the lowest RMSE.

This helps to fine-tune the model for better generalization and performance.

In [25]:
from sklearn.model_selection import GridSearchCV

# Define the parameter grid to search
param_grid = {
    'n_estimators': [100, 200],
    'max_depth': [10, 20, None],
    'min_samples_split': [2, 5]
}

# Create Griddsearchcv object
rf_grid = GridSearchCV(
    estimator=RandomForestRegressor(random_state=42),
    param_grid=param_grid,
    scoring='neg_mean_squared_error',
    cv=5,
    n_jobs=-1,         # use all cores
    verbose=1
)

# Fitting  grid search to the training data
rf_grid.fit(X_train, y_train)

#  show the best parameters  
print("✅ Best Parameters found:")
print(rf_grid.best_params_)

best_rmse = (-rf_grid.best_score_) ** 0.5
print(f"Best Cross-Validated RMSE: {best_rmse:.2f}")
Fitting 5 folds for each of 12 candidates, totalling 60 fits
✅ Best Parameters found:
{'max_depth': 20, 'min_samples_split': 2, 'n_estimators': 200}
Best Cross-Validated RMSE: 763.65
In [26]:
# Evaluate the tuned Random Forest model

from sklearn.metrics import r2_score, mean_squared_error

y_pred_best = rf_grid.best_estimator_.predict(X_test)
r2_rf_tuned = r2_score(y_test, y_pred_best)
rmse_rf_tuned = mean_squared_error(y_test, y_pred_best) ** 0.5


# Store all scores, including the tuned RF model

r2_scores = {
    'Linear Regression': lr_r2,
    'Decision Tree': r2_tree,
    'Random Forest (Default)': r2_rf,
    'Random Forest (Tuned)': r2_rf_tuned,
    'Gradient Boosting': r2_gb  }


rmse_scores = {
    'Linear Regression': lr_rmse,
    'Decision Tree': rmse_tree,
    'Random Forest (Default)': rmse_rf,
    'Random Forest (Tuned)': rmse_rf_tuned,
    'Gradient Boosting': rmse_gb  }


# Bar plot for R2

plt.figure(figsize=(9, 5))
sns.barplot(x=list(r2_scores.keys()), y=list(r2_scores.values()))

plt.title("R² Score Comparison")
plt.ylabel("R² Score")
plt.ylim(0, 1)
plt.xticks(rotation=15)


plt.show()

# bar  plot for RMSE

plt.figure(figsize=(9, 5))
sns.barplot(x=list(rmse_scores.keys()), y=list(rmse_scores.values()))

plt.title("RMSE Comparison")
plt.ylabel("RMSE")
plt.xticks(rotation=15)


plt.show()
No description has been provided for this image
No description has been provided for this image

final model performance comparison

I just compared the performance of all the models I trained using two metrics:

R2 Score, which shows how well the model explains the variability of the data (higher is better), and
RMSE (Root Mean Squared Error), which shows the average prediction error in the same units as the target variable (lower is better).

As shown in the bar plots below, the tuned Random forest model performs better than the default one, confirming that hyperparameter tuning improves prediction accuracy.

Gradient Boosting also performs very well, with slightly better RMSE but similar R2.

Linear Regression and Decision tree models performed worse in comparison.

Now I will try a Neural Network model to see how well it can learn the pattern.

I need to scale the target variable"y" because neural networks work best when both input and output are normalized

Let’s see how it performs.

In [27]:
# importing everything we need
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import r2_score, mean_squared_error

#creating and train the neural network model

nn = MLPRegressor(hidden_layer_sizes=(100,100), max_iter=3000, random_state=42)

nn.fit(X_train, y_train)

# make predictions
y_pred_nn = nn.predict(X_test)

#evaluate the model
r2_nn = r2_score(y_test, y_pred_nn)
rmse_nn = mean_squared_error(y_test, y_pred_nn) ** 0.5

# print results

print("Neural Network R2:", round(r2_nn, 3))

print("Neural Network RMSE:", round(rmse_nn, 2))
Neural Network R2: 0.893
Neural Network RMSE: 609.03
/usr/local/lib/python3.12/dist-packages/sklearn/neural_network/_multilayer_perceptron.py:785: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (3000) reached and the optimization hasn't converged yet.
  warnings.warn(

Residual Analysis for Neural Network

This residual plot shows how well the model fits the data.

Ideally, the residuals should be randomly distributed around 0.

If we see patterns, it might indicate bias or model limitations.

In [28]:
residuals = y_test - y_pred_nn

plt.figure(figsize=(8, 5))
plt.scatter(y_pred_nn, residuals, alpha=0.5)
plt.axhline(y=0, color='red', linestyle='--')
plt.title('Residual Plot - Neural Network')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.grid(True)

plt.show ( )
No description has been provided for this image
In [29]:
r2_scores["Neural Network"] = round(r2_nn, 3)   # I update the dictionaries with Neural Network results

rmse_scores["Neural Network"] = round(rmse_nn, 2)

# Re-plot updated R2 scores

plt.figure(figsize=(8, 4))
sns.barplot(x=list(r2_scores.keys()), y=list(r2_scores.values()))

plt.title("R² Score Comparison")
plt.ylabel("R² Score")
plt.ylim(0, 1)

plt.show()

# plot updated RMSE score

plt.figure(figsize=(8, 4))
sns.barplot(x=list(rmse_scores.keys()), y=list(rmse_scores.values()))

plt.title("RMSE Comparison")
plt.ylabel("RMSE")

plt.show()
No description has been provided for this image
No description has been provided for this image

Error Analysis actual vs Predicted

To better understand where the model performs well or struggles, I plotted the actual vs predicted values and the residuals (prediction errors).

In [30]:
# Predict with the best model on the test set


y_pred_best = nn.predict(X_test)  # the Neural Network was actually the best model overall

# plot Actual vs Predicted

plt.figure(figsize=(8, 6))
sns.scatterplot(x=y_test, y=y_pred_best, alpha=0.5)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--')  # Diagonal line
plt.xlabel("Actual Values")
plt.ylabel("Predicted Values")
plt.title("Actual vs. Predicted Values (Neural Network)")

plt.grid(True)

plt.show()

#plot Residuals

residuals = y_test - y_pred_best

plt.figure(figsize=(8, 5))
sns.histplot(residuals, bins=30, kde=True)
plt.title("Distribution of Residuals")
plt.xlabel("Prediction Error")
plt.ylabel("Frequency")

plt.grid(True)

plt.show()
No description has been provided for this image
No description has been provided for this image

Final conclusion!

Here I just completed the entire machine learning pipeline, starting from exploring the data and building new features all the way to training and evaluating multiple regression model.

Throughout the project I tested different algorithms including linear regression, decision tree, random forest (both default and tuned), Gradient boosting, and a Neural Network, and compared their performance using metrics like R2 and RMSE.

The best overall performance was achieved by the Neural Network, which had the highest R2 and the lowest RMSE among all the models tested. Among the tree based models, the tuned Random Forest was the strongest, benefiting from hyperparameter optimization using GridSearchcv, which shows how useful model tuning can be.

The error analysis confirmed that the model does quite a good job across most demand levels, although predictions can still be off for very high or very low rental days. The residuals were roughly normally distributed, which is a good sign of a stable model.

Key learnings:

Feature engineering (like extracting time variables) adds real value. Cross-validation is essential for fair model comparison. Hyperparameter tuning can boost accuracy significantly. Visualizing prediction errors helps diagnose model behavior.

Final thoughts:

With more time, I would try ensembling models or incorporating more complex features (like lag features or holiday flags).

Overall, this project gave me hands on experience applying regression techniques seen during classes on a real world dataset, it also helped me while studying for the exam to understand better certain concepts.