<jupyter_start><jupyter_text>Capital Bike Share Demand (Regression)  Import data<jupyter_code>import pandas as pd
data = pd.read_csv('train.csv')
data<jupyter_output><empty_output><jupyter_text>Description of variables- 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 rentals<jupyter_code>data.isnull().sum()<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis  Which factors contribute most to the number of bicycles being checked out over the course of a given day?<jupyter_code># Convert datetime to datetime format
data['datetime']= pd.to_datetime(data['datetime'])
data
data.count()<jupyter_output><empty_output><jupyter_text>How much of an impact does weather have on demand?  - Weather has a negative and week correlation (-0.13) on demand, i.e. if weather is bad (e.g. 4) than the demand goes down - in contrast temp, atempand humidity have higher impact than weather<jupyter_code>import seaborn as sns
sns.barplot(x= 'weather', y= 'count', data = data)<jupyter_output><empty_output><jupyter_text>Datetime conversion<jupyter_code>data['year'] = data['datetime'].dt.year
data['month'] = data['datetime'].dt.month
data['day'] = data['datetime'].dt.weekday
data['hour'] = data['datetime'].dt.hour
data
sns.barplot(x= 'year', y= 'count', data = data)
sns.barplot(x= 'month', y= 'count', data = data)
sns.barplot(x= 'day', y= 'count', data = data, palette= 'Set2')
sns.barplot(x= 'hour', y= 'count', data = data)<jupyter_output><empty_output><jupyter_text>How does the behavior of casual users differ from subscribers?<jupyter_code>month = data.groupby(['month'])['registered','casual'].sum()
month
month.reset_index(inplace=True)
month
sns.lineplot(data=month, y='registered', x = 'month')
sns.lineplot(data=month, y='casual', x = 'month')
import matplotlib.pyplot as plt

y = month['registered']
x = month['month']
w = month['casual']


plt.bar(x,y, align='center')
plt.bar(x,w, align='center')


day = data.groupby(['day'])['registered','casual'].sum()
day.reset_index(inplace=True)

sns.lineplot(data=day, y='registered', x = 'day')
sns.lineplot(data=day, y='casual', x = 'day')
data.drop(labels='datetime', axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>PP Score https://towardsdatascience.com/rip-correlation-introducing-the-predictive-power-score-3d90808b9598  works only with older version of scikit-learn<jupyter_code>#import ppscore as pps

#fig, ax = plt.subplots(figsize=(10,10))  
#a = pps.matrix(data).pivot(columns='x', index='y', values='ppscore')
#sns.heatmap(a, annot=True, ax=ax, cmap='coolwarm')<jupyter_output><empty_output><jupyter_text>Heatmap<jupyter_code>fig, ax = plt.subplots(figsize=(10,10))  
sns.heatmap(data.corr(), annot= True, ax=ax, cmap= 'coolwarm')<jupyter_output><empty_output><jupyter_text>--> Order of correlation, starting with highest: registered, casual,temp and atemp, humidity, season, weather, workingday, windspeed, holiday<jupyter_code>from IPython.display import display, Markdown, Latex<jupyter_output><empty_output><jupyter_text>Checking for multicollinearity with Variance Inflation Factor  VIF measures the linear association between an independent variable and all other independent variables. A VIF for any given independent variable is calculated by:                                    $VIF = \frac{1}{1-R²}$  https://towardsdatascience.com/everything-you-need-to-know-about-multicollinearity-2f21f082d6dc<jupyter_code>import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor<jupyter_output><empty_output><jupyter_text>VIF above 10 indicates high multicollinearity<jupyter_code>data= data[~data.isin([np.nan, np.inf, -np.inf]).any(1)]
X_variables = data[['season', 'holiday', 'workingday', 'weather', 'temp', 'atemp','humidity', 'windspeed', 
                    'year','month', 'day', 'hour']]

vif_data = pd.DataFrame()
vif_data["feature"] = X_variables.columns
vif_data["VIF"] = [variance_inflation_factor(X_variables.values, i) 
                   for i in range(len(X_variables.columns))]
vif_data
#removed season, year and atemp

data= data[~data.isin([np.nan, np.inf, -np.inf]).any(1)]
X_variables = data[[ 'holiday', 'workingday','temp', 'windspeed', 'humidity',
                   'month', 'day', 'hour', 'weather']]

vif_data = pd.DataFrame()
vif_data["feature"] = X_variables.columns
vif_data["VIF"] = [variance_inflation_factor(X_variables.values, i) 
                   for i in range(len(X_variables.columns))]
vif_data<jupyter_output><empty_output><jupyter_text>Plots of y and X variables<jupyter_code>f, axes = plt.subplots(3, 3, figsize=(12, 14), sharex=False)
a = sns.lineplot(data=data, x="holiday", y="count", ax=axes[0, 0])
b = sns.lineplot(data=data, x="workingday", y="count", ax=axes[0, 1])
c = sns.lineplot(data=data, x="weather", y="count", ax=axes[0, 2])
d = sns.lineplot(data=data, x="temp", y="count", ax=axes[1, 0])
e = sns.lineplot(data=data, x="humidity", y="count", ax=axes[1, 1])
g = sns.lineplot(data=data, x="windspeed", y="count", ax=axes[1, 2])
h = sns.lineplot(data=data, x="month", y="count", ax=axes[2, 0])
i = sns.lineplot(data=data, x="day", y="count", ax=axes[2, 1])
j = sns.lineplot(data=data, x="hour", y="count", ax=axes[2, 2])<jupyter_output><empty_output><jupyter_text>Train test split<jupyter_code>y = data['count']
y.shape
X = data[['holiday', 'workingday', 'weather', 'temp','humidity', 'windspeed', 'day', 'hour', 'month', 'season']]
X.shape
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.2, random_state = 42)
X_train.shape, X_test.shape, y_train.shape, y_test.shape
X_train.head()
X_test.head()<jupyter_output><empty_output><jupyter_text>Feature Engineering  Preprocessor for windspeed and humidity<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
poly2_transformer = PolynomialFeatures(degree= 2, include_bias= False)<jupyter_output><empty_output><jupyter_text>Preprocessor for month, hour and weather<jupyter_code>from sklearn.preprocessing import OneHotEncoder, MinMaxScaler
encoder = OneHotEncoder(sparse=False, handle_unknown='ignore')
scaler = MinMaxScaler()<jupyter_output><empty_output><jupyter_text>Preprocessor for temp and humidity<jupyter_code>from sklearn.pipeline import make_pipeline
action = make_pipeline(PolynomialFeatures(degree= 3, include_bias= False), PolynomialFeatures(interaction_only=True, include_bias= False))
from sklearn.preprocessing import KBinsDiscretizer
bin_transformer = KBinsDiscretizer(n_bins= 7, encode = 'onehot-dense', strategy = 'quantile')<jupyter_output><empty_output><jupyter_text>ColumnTransformer<jupyter_code>from sklearn.compose import ColumnTransformer

feature_transform = ColumnTransformer([
    ('do_nothing', 'passthrough', ['workingday', 'day']),
    #('poly2', poly2_transformer, [ 'windspeed', 'humidity']),
    ('inter', action, ['windspeed', 'humidity', 'temp']),
    #('temp', bin_transformer, ['temp']),
    #('scaler', scaler, ['weather']),
    ('encoder', encoder, ['month', 'hour','weather']),
    #('hours', hours_transformer, ['hour']),
    #('drop_features','drop',['day', 'hour']),

])

X_train_trans = feature_transform.fit_transform(X_train)
X_train_trans<jupyter_output><empty_output><jupyter_text>Fit Linear Regression<jupyter_code>import numpy as np
y_train_log = np.log1p(y_train)
y_test_log = np.log1p(y_test)
from sklearn.linear_model import LinearRegression
m = LinearRegression()

#ylog_train = np.log1p(y_train)
m.fit(X_train_trans, y_train_log)
round(m.score(X_train_trans, y_train_log),3)
X_test_trans = feature_transform.transform(X_test)
round(m.score(X_test_trans, y_test_log),3)
y_pred = m.predict(X_test_trans)
y_pred2 = np.exp(y_pred)-1
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score, mean_squared_log_error
print(f'R^2: {round(r2_score(y_test, y_pred2),2)}')
print(f'MSE: {round(mean_squared_error(y_test, y_pred2),2)}')
print(f'MAE: {round(mean_absolute_error(y_test, y_pred2),2)}')
print(f'RMSLE: {round(mean_squared_log_error(y_test, y_pred2),2)}')<jupyter_output>RMSLE: 0.39<jupyter_text>OLS output<jupyter_code>import statsmodels.api as sm
from scipy import stats
import numpy as np
#X_train_trans.shape, y_train.shape
mod = sm.OLS(y_train_log, X_train_trans)
res = mod.fit()
print(res.summary())
#X_train_trans<jupyter_output><empty_output><jupyter_text>Poisson regression<jupyter_code>P = data[['holiday', 'workingday', 'weather', 'temp','humidity', 'windspeed','month', 'day', 'hour']]
P.shape
from sklearn.model_selection import train_test_split

P_train, P_test, y_train, y_test = train_test_split(P,y, test_size = 0.2, random_state = 42)
P_train.shape, P_test.shape, y_train.shape, y_test.shape
P_train.head()
P_test.head()
P_train_trans = feature_transform.fit_transform(P_train)
P_train_trans
P_test_trans = feature_transform.transform(P_test)
from sklearn.linear_model import PoissonRegressor

mp = PoissonRegressor(alpha=0.000001)
mp.fit(P_train_trans, y_train)

y_train.min()
p_pred = mp.predict(P_test_trans)
p_pred
mp.score(P_train_trans, y_train)
from sklearn.metrics import mean_poisson_deviance #Mean Poisson deviance regression loss
print (f'Mean Poisson deviance regression loss: {round(mean_poisson_deviance(y_test, p_pred),2)}')
print(f'RMSE: {round(mean_squared_error(y_test, p_pred),2)}') 
print(f'MSE: {round(mean_squared_error(y_test, p_pred),2)}')
print(f'MAE: {round(mean_absolute_error(y_test, p_pred),2)}')
print(f'RMSLE: {round(mean_squared_log_error(y_test, p_pred),2)}')<jupyter_output>RMSLE: 2.52<jupyter_text>Lasso regression<jupyter_code>from sklearn.linear_model import Lasso
m_lasso = Lasso(alpha = 1)
m_lasso.fit(X_train_trans,y_train)
ypred_lasso = m_lasso.predict(X_test_trans)
ypred_lasso2 = np.exp(y_pred)-1
ypred_lasso2
m_lasso.coef_ 
round(m_lasso.score(X_train_trans, y_train),3)
round(m_lasso.score(X_test_trans, y_test),3)
print(f'R^2: {round(r2_score(y_test, ypred_lasso2),2)}')
print(f'MSE: {round(mean_squared_error(y_test, ypred_lasso),2)}')
print(f'MAE: {round(mean_absolute_error(y_test, ypred_lasso),2)}')
#X_train_trans = pd.DataFrame(X_train_trans, columns= feature_transform.get_feature_names_out())
#X_train_trans
#X_train_trans.columns
#pd.DataFrame({'features': X_train_trans.columns, 'weights': m_lasso.coef_})<jupyter_output><empty_output><jupyter_text>Random Forest<jupyter_code>from sklearn.ensemble import RandomForestRegressor

rf = RandomForestRegressor(n_estimators= 100,random_state = 42) 
rf_log = RandomForestRegressor(n_estimators= 100,random_state = 42) 
rf.fit(X_train_trans,y_train)

rf_log.fit(X_train_trans, y_train_log)

y_pred_rf = rf.predict(X_test_trans)
y_pred_rf_1 = rf_log.predict(X_test_trans)

#print(f'R^2: {round(r2_score(y_test, y_predrf),2)}')
#print(f'MSE: {round(mean_squared_error(y_test, y_predrf),2)}')
#print(f'RMSE: {round(mean_squared_log_error(y_test, y_pred),2)}') ##log before??
round(rf.score(X_train_trans, y_train),2)
round(rf.score(X_test_trans, y_test),2)
round(rf_log.score(X_train_trans, y_train_log),2)
round(rf_log.score(X_test_trans, y_test_log),2)
y_pred_rf_2 = np.exp(y_pred_rf_1)-1
print(f'MAE: {round(mean_absolute_error(y_test,y_pred_rf),2)}')
print(f'MAE: {round(mean_absolute_error(y_test,y_pred_rf_2),2)}')
print(f'RMSLE: {round(mean_squared_log_error(y_test, y_pred_rf),2)}')
print(f'RMSLE: {round(mean_squared_log_error(y_test, y_pred_rf_2),2)}')
print(f'R^2: {round(r2_score(y_test, y_pred_rf_2),2)}')
from sklearn.model_selection import cross_val_score
#For the forest
cross_val_score(rf, X_train_trans, y_train)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1457536031267812, 'fraction_numerical': 0.01335396140379448, 'mean_word_length': 7.7979942693409745, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 88, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '528723', 'license': 'MIT', 'n_tokens_mistral': 4746, 'n_tokens_neox': 4550, 'n_words': 875, 'path': 'Bike_Sharing _Demand_project.ipynb', 'repo_name': 'fasoula/bike_sharing_demand_prediction', 'size': '12281'}