[Data/Data Analysis] Kaggle : 자전거 수요 예측 (Bike Sharing Demand)

---

 

Bike Sharing Demand 

 

https://www.kaggle.com/c/bike-sharing-demand

 

https://www.kaggle.com/c/bike-sharing-demand

Bike Sharing Demand | Kaggle

 

자전거 수요 예측 모델을 만들어보자.

만들긴 만드는데 내가 만든 건 아니고 다른분들 것 필사하며 공부하기 ^^~

이 competition은 7-8년전이라 이미 끝난 대회이다.

자전거 공유 시스템의 수요를 예측하는 것이 목적.

그럼 미국에는 이 때도 자전거 공유 시스템이 있었다는 건데 ..

우린 카카오바이크 같은 게 상용화된 지 얼마 안 된 것 같은데 신기하네 

 

 

데이터 살펴보기

---

Data Dictionary 를 먼저 살펴보고 들어가보자.

 

---

• Datetime : Hourly Date + timestamp
• season : 1 = spring, 2 = summer, 3 = fall, 4 = winter
• holiday : holiday
• workingday : whether the day is neither a weekend nor holiday
• weather : 1 = clear, few clouds, partly cloudy, 대체로 좋은 기후
  2 = Mist + Cloudy, Mist, 살짝 흐린 기후
  3 = Light snow, Light Rain 약하게 안 좋은 기후
  4 = Heavy rain, Ice Pallets + Thunderstorm + Mist, Snow + Fog 대충 악기후 말하는듯.
• temp : temperature in Celsius
• atemp : "feels like" temperature in Clesius, 체감온도
• humidity 
• windspeed
• casual : number of non-registered user rentals initiated 미등록 유저
• registered : 등록유저
• count : 우리의 Target 변수. number of total rentals

---

 

먼저 필요한  library  들을 import 하고 들어가자

 

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import os
from scipy import stats

import missingno as msno
plt.style.use('seaborn')

import warnings
warnings.filterwarnings("ignore")

mpl.rcParams['axes.unicode_minus'] = False

%matplotlib inline

mlp.rcParams['axes.unicode_minus'] = False 는 그래프에서 마이너스가 깨지는 것을 방지하기 위한 코드이다.

 

# 데이터 불러오기
df_train = pd.read_csv("./train.csv", parse_dates = ["datetime"])
df_test = pd.read_csv("./test.csv", parse_dates = ["datetime"])

parse_dates 를 이용해서 날짜와 시간 변수를 datetime 으로 변경하며 가져왔다.

 

registered 변수는 train 데이터셋에만 존재하는 것을 확인했다.

count는 우리의 target 변수

 

train 데이터 확인.

 

msno.matrix(df_train, figsize = (12,5))

missingno를 사용하여 결측치를 시각화해주었다.

꽉꽉 들어차있는 것을 보아 결측치는 없는 것으로 보인다.

 

Target 변수인 count의 왜도와 첨도를 확인해보자.

f,ax = plt.subplots(1, 1, figsize = (10,6))
g = sns.distplot(df_train["count"], color = "b", label = 
                 "Skewness : {:2f}".format(df_train["count"].skew()), ax=ax)
g = g.legend(loc = "best")

print("Skewness : %f" % df_train['count'].skew())
print("Kurtosis : %f" % df_train['count'].kurt())

Skewness : 1.242066
Kurtosis : 1.300093

왜도가 양수인 것을 보아 왼쪽으로 치우쳐져있다는 것을 알 수 있다. 

또한 0에 값들이 집중되어있는 것이 확인됐다.

 

 

EDA

---

Datetime 변수를 먼저 살펴보자.

Datetime 을 연,월,일,시,분,초 단위로 쪼개주었다.

# Datetime 먼저 살펴보기
# 연,월,일,시,분,초 로 나누어준다
df_train['year'] = df_train['datetime'].dt.year
df_train['month'] = df_train['datetime'].dt.month
df_train['day'] = df_train['datetime'].dt.day
df_train['hour'] = df_train['datetime'].dt.hour
df_train['minute'] = df_train['datetime'].dt.minute
df_train['second'] = df_train['datetime'].dt.second

df_test['year'] = df_test['datetime'].dt.year
df_test['month'] = df_test['datetime'].dt.month
df_test['day'] = df_test['datetime'].dt.day
df_test['hour'] = df_test['datetime'].dt.hour
df_test['minute'] = df_test['datetime'].dt.minute
df_test['second'] = df_test['datetime'].dt.second

각각 year month day hour minute second 칼럼이 새로 생성된 것을 확인할 수 있다.

 

barplot 을 이용하여 각각의 칼럼을 시각화해서 비교했다.

# barplot 을 이용하여 시각화

figure, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows = 2, ncols = 3)
figure.set_size_inches(18,10)

sns.barplot(data = df_train, x = "year", y="count", ax = ax1)
sns.barplot(data = df_train, x = "month", y="count", ax = ax2)
sns.barplot(data = df_train, x = "day", y="count", ax = ax3)
sns.barplot(data = df_train, x = "hour", y="count", ax = ax4)
sns.barplot(data = df_train, x = "minute", y="count", ax = ax5)
sns.barplot(data = df_train, x = "second", y="count", ax = ax6)

ax1.set(ylabel = "count", title = "Rental amount by year")
ax2.set(ylabel = "count", title = "Rental amount by month")
ax3.set(ylabel = "count", title = "Rental amount by day")
ax4.set(ylabel = "count", title = "Rental amount by hour")

그래프를 확인해봤을 때, 2012년에 수요량이 더 상승했다.

겨울보다는 5-10월, 봄부터 가을까지 수요량이 많다.

날짜 별로는 거의 차이가 없다.

시간대별로는 8시와 17-18시까지 수요가 많은 것을 보아 출퇴근시간에 수요가 늘어난다는 것을 알 수 있다.

 

요일변수를 추가로 생성해서 요일별로 확인해보자.

# 요일변수를 추가로 생성
df_train['dayofweek'] = df_train['datetime'].dt.dayofweek
df_test['dayofweek'] = df_test['datetime'].dt.dayofweek
df_train.head()

dayofweek 칼럼이 새로 생겼다.

0은 월요일 6은 일요일이다.

 

df_train['dayofweek'].value_counts()

5    1584
6    1579
3    1553
0    1551
2    1551
1    1539
4    1529
Name: dayofweek, dtype: int64

크게 차이는 없어보이나 5, 6에 조금 더 많은 수요가 있다. 아마도 주말에 수요가 더 많은 것 같다.

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(nrows = 5)
fig.set_size_inches(18,25)

sns.pointplot(data = df_train, x = "hour", y="count", ax = ax1)
sns.pointplot(data = df_train, x = "hour", y="count", hue = 'workingday', ax= ax2)
sns.pointplot(data = df_train, x = "hour", y="count", hue = 'dayofweek', ax = ax3)
sns.pointplot(data = df_train, x = "hour", y="count", hue = 'weather', ax = ax4)
sns.pointplot(data = df_train, x = "hour", y="count", hue = 'season', ax = ax5)

시간의 흐름에 따라서 pointplot 을 그려주었다.

주말에는 오후대에 수요가 가장 많아지는 것을 확인할 수 있다.

나머지는 우리가 위에서 확인한 것과 비슷하다.

악천후일때는 거의 자전거를 타지 않고 날씨가 좋을수록 수요도 높은 것이 눈에 보인다.

계절로 보았을 땐 의외로 봄에 제일 수요가 적고 가을,여름,겨울순인 것을 알 수 있다.

변수들 간 상관관계를 확인해보았다.

corr_data = df_train[['temp','atemp','casual','registered','humidity','windspeed','count']]

colormap = plt.cm.PuBu

f, ax = plt.subplots(figsize = (12,10))
plt.title("Correlation of Numeric Features with Rental Count", y=1, size=18)
sns.heatmap(corr_data.corr(), vmax=0.8, linewidths=0.1, square=True, annot = True,
            cmap=colormap, linecolor="white", annot_kws={'size':14})

plt.show()

우리의 Target변수인 count와 비교했을 때 가장 높은 상관관계를 보인 것은 registered. 그 다음은 casual이다.

temp와 atemp 둘 또한 다중공선성이 있는 것을 확인했다. 어쩌면 당연함 ..

count변수와 다른 습도, 풍속, 온도와도 상관관계가 없는 것을 확인했다.

 

우리가 유의하지 않다고 판단한 변수들에 대해서 scatter plot 을 그려보자.

# Heatmap을 통해 유의하지 않다고 판단된 변수에 대해서 scatter plot을 그려본다.
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(12,5))

temp_scatter_plot = pd.concat([df_train['count'],df_train['temp']], axis=1)
sns.regplot(x='temp', y='count', data = temp_scatter_plot, scatter=True, fit_reg=True, ax=ax1)
windspeed_scatter_plot = pd.concat([df_train['count'],df_train['windspeed']], axis=1)
sns.regplot(x='windspeed', y='count', data= windspeed_scatter_plot, scatter=True, 
            color = 'green', fit_reg=True, ax=ax2)
humidity_scatter_plot = pd.concat([df_train['count'], df_train['humidity']], axis=1)
sns.regplot(x='humidity', y='count', data=humidity_scatter_plot, scatter=True, 
            color = 'red', fit_reg=True, ax=ax3)

풍속변수를 보니 0에 값들이 몰려있는 것을 확인할 수 있다.

windspeed 변수를 좀 더 구체화해서 살펴보자.

fig, axes = plt.subplots(nrows=2, figsize= (18,14))

plt.sca(axes[0])
plt.xticks(rotation=30, ha="right")
axes[0].set(ylabel="count", title = "train windspeed")
sns.countplot(data = df_train, x="windspeed", ax=axes[0])

plt.sca(axes[1])
plt.xticks(rotation = 30, ha='right')
axes[1].set(ylabel="count", title='test windspeed')
sns.countplot(data=df_test, x= 'windspeed', ax=axes[1])

plt.show()

0에 있는 값들이 너무 많다.

나중에 수정해주도록 하자.

 

연도별 데이터의 대여 변화를 좀 더 자세하게 보기 위해서 year와 month 데이터를 붙여준다.

def concatenate_year_month(datetime):
    return "{0}-{1}".format(datetime.year, datetime.month)

df_train['year_month'] = df_train['datetime'].apply(concatenate_year_month)
df_test['year_month'] = df_test['datetime'].apply(concatenate_year_month)

print(df_train.shape)
df_train[['datetime','year_month']].head()

(10886, 20)

 

fig, ax = plt.subplots(figsize = (18,4))
sns.barplot(data =df_train, y="count", x="year_month")

확실히 2012년이 count가 훨씬 많다.

detect_outlier 함수를 만들어서 이상치를 제거해준다.

이상치는 IQR 을 이용하여 그 이상이나 이하는 삭제해주는 방법을 사용했다.

from collections import Counter 

def detect_outliers(df, n, features):
    outlier_indices = []
    for col in features:
        Q1 = np.percentile(df[col], 25)
        Q3 = np.percentile(df[col], 75)
        IQR = Q3 - Q1
        
        outlier_step = 1.5 * IQR
        
        outlier_list_col = df[(df[col] < Q1 - outlier_step) | (df[col] > Q3 + outlier_step)].index
        outlier_indices.extend(outlier_list_col)
    outlier_indices = Counter(outlier_indices)
    multiple_outliers = list(k for k, v in outlier_indices.items() if v > n)
        
    return multiple_outliers
        
Outliers_to_drop = detect_outliers(df_train, 2, ["temp", "atemp", "casual", 
                                                 "registered", "humidity", "windspeed", "count"])

df_train = df_train.drop(Outliers_to_drop, axis = 0).reset_index(drop=True)
df_train.shape

(10846, 20)

삭제해주었다.

행의 개수가 좀 줄어든 것을 확인할 수 있다.

40개정도 ..

 

다시 한 번 왜도와 첨도를 확인한다.

for col in df_train_num :
    print('{:15}'.format(col),
         'Skewness : {:05.2f}'.format(df_train[col].skew()),
         '  ',
         "Kurtosis : {:06.2f}".format(df_train[col].kurt()))

count           Skewness : 01.21    Kurtosis : 001.20
temp            Skewness : 00.01    Kurtosis : -00.91
atemp           Skewness : -0.10    Kurtosis : -00.85
casual          Skewness : 02.52    Kurtosis : 007.74
registered      Skewness : 01.51    Kurtosis : 002.61
humidity        Skewness : -0.09    Kurtosis : -00.76
windspeed       Skewness : 00.58    Kurtosis : 000.63

boxplot을 통해 확인해보자.

 

fig, axes = plt.subplots(nrows = 5, ncols = 2, figsize = (16,18))
sns.boxplot(data = df_train, y="count", x = "season", orient="v", ax = axes[0][0])
sns.boxplot(data = df_train, y="count", x = "holiday", orient="v", ax = axes[0][1])
sns.boxplot(data = df_train, y="count", x = "workingday", orient="v", ax = axes[1][0])
sns.boxplot(data = df_train, y="count", x = "weather", orient="v", ax = axes[1][1])
sns.boxplot(data = df_train, y="count", x = "dayofweek", orient="v", ax = axes[2][0])
sns.boxplot(data = df_train, y="count", x = "month", orient="v", ax = axes[2][1])
sns.boxplot(data = df_train, y="count", x = "year", orient="v", ax = axes[3][0])
sns.boxplot(data = df_train, y="count", x = "hour", orient="v", ax = axes[3][1])
sns.boxplot(data = df_train, y="count", x = "minute", orient="v", ax = axes[4][0])

axes[0][0].set(ylabel = "count", title="Rental count by season")
axes[0][1].set(ylabel = "count", title="Rental count by holiday")
axes[1][0].set(ylabel = "count", title="Rental count by workingday")
axes[1][1].set(ylabel = "count", title="Rental count by weather")
axes[2][0].set(ylabel = "count", title="Rental count by dayofweek")
axes[2][1].set(ylabel = "count", title="Rental count by month")
axes[3][0].set(ylabel = "count", title="Rental count by year")
axes[3][1].set(ylabel = "count", title="Rental count by hour")
axes[4][0].set(ylabel = "count", title="Rental count by minute")

통계량을 보았을 때 우리가 위에서 시각화하여 살펴본 것과 유사하다. 굳굳 ..

 

 

 

 

전처리

---

f,ax = plt.subplots(1, 1, figsize = (10,6))
g = sns.distplot(df_train["count"], color = "b", label = 
                 "Skewness : {:2f}".format(df_train["count"].skew()), ax=ax)
g = g.legend(loc = "best")

print("Skewness : %f" % df_train['count'].skew())
print("Kurtosis : %f" % df_train['count'].kurt())

Skewness : 1.210923
Kurtosis : 1.200871

한 번 더 왜도와 첨도를 확인했다.

 

Log 변환을 통해서 스케일링 해주었다.

df_train['count_Log'] = df_train['count'].map(lambda i:np.log(i) if i>0 else 0)

f,ax = plt.subplots(1,1, figsize = (10,6))
g = sns.distplot(df_train["count_Log"], color = "b", label = 
                 "Skewness : {:2f}".format(df_train["count"].skew()), ax=ax)
g = g.legend(loc = "best")

print("Skewness : %f" % df_train['count_Log'].skew())
print("Kurtosis : %f" % df_train['count_Log'].kurt())

Skewness : -0.975198
Kurtosis : 0.247435

df_train.drop('count', axis=1, inplace = True)

다음은 풍속을 Feature Engineering 할 차레이다.

우리는 위에서 풍속이 0값이 많은 것을 확인했다.

상식적으로 풍속이 0일 가능성은 거의 없기 때문에 우리는 0값을 다른 값으로 대체해줄 것이다.

그렇다면 어떻게 ?

많은 방법들이 있지만 우리는 RandomForest를 이용하여 풍속을 Target변수로, 나머지를 독립변수로 사용하여

풍속을 예측하여 값을 채워준 뒤 사용할 것이다.

 

풍속이 0이 아닌 값들을 학습시키고,

그 모델에 풍속이 0인 데이터들을 넣어서 예측할 것이다.

trainwind0 = df_train.loc[df_train['windspeed'] == 0]
trainWindNot0 = df_train.loc[df_train['windspeed'] != 0 ]

trainwind0에 windspeed가 0인 데이터를 넣어준다.

tainwindNot0에는 windspeed가 0이 아닌 데이터들을 넣어준다.

from sklearn.ensemble import RandomForestClassifier

def predict_windspeed(data):
    dataWind0 = data.loc[data['windspeed'] == 0]
    dataWindNot0 = data.loc[data['windspeed'] != 0]
    
    wcol = ['season', 'weather', 'humidity', 'day', 'temp', 'atemp']
    
    dataWindNot0['windspeed'] = dataWindNot0['windspeed'].astype('str') #데이터 형식 바꾸기
    
    rf_wind = RandomForestClassifier()
    rf_wind.fit(dataWindNot0[wcol], dataWindNot0['windspeed']) #독립변수를 wcol애들로, 종속변수를 windspeed로 하여 학습
    wind0 = rf_wind.predict(X=dataWind0[wcol]) # rf_wind 모델에 wind0 wcol데이터를 넣어 wind0에 저장
    
    predictWind0 = dataWind0 # predictWind0에 datawind0를 넣어준다.
    predictWindNot0 = dataWindNot0 
    
    predictWind0['windspeed'] = wind0 #predictWind0의 windspeed에 예측한 wind0을 넣어준다.
    
    data = predictWindNot0.append(predictWind0) #data에 predictWind0를 추가해준다(행추가)
    
    data['windspeed'] = data['windspeed'].astype('float') #데이터타입 바꿔주기
    
    data.reset_index(inplace = True) #index 초기화
    data.drop('index', inplace = True, axis = 1) #이후 index칼럼 제거
    
    return data

train data와 test data를 각각 predict_windspeed 함수에 넣어준 뒤 결과를 시각화시켰다.

df_train = predict_windspeed(df_train)
df_test = predict_windspeed(df_test)

fig, (ax1, ax2) = plt.subplots(nrows = 2, figsize = (18,14))

plt.sca(ax1)
plt.xticks(rotation = 30, ha = "right")
ax1.set(ylabel = "count", title = "train windspeed")
sns.countplot(data = df_train, x = "windspeed", ax = ax1)

plt.sca(ax2)
plt.xticks(rotation = 30, ha = "right")
ax2.set(ylabel = "count", title = "test windspeed")
sns.countplot(data = df_test, x = "windspeed", ax = ax2)

이제 0값은 없다.

상관계수를 간단히 확인해주었다.

# 상관계수 확인
corr_data = df_train[["count_Log", "windspeed"]]
corr_data.corr()

 

범주형 변수들을 원핫인코딩 시켰다.

weather 와 season 을 one-hot encoding

df_train = pd.get_dummies(df_train, columns = ["weather"], prefix = "weather")
df_test = pd.get_dummies(df_test, columns = ["weather"], prefix = "weather")

df_train = pd.get_dummies(df_train, columns = ["season"], prefix = "season")
df_test = pd.get_dummies(df_test, columns = ["season"], prefix = "season")

# train, test 정리하고 마무리 
datetime_test = df_test['datetime']

df_train.drop(["datetime", "registered","casual","holiday", "year_month", "minute", "second"], axis = 1, inplace = True)
df_test.drop(["datetime","holiday", "year_month", "minute", "second"], axis = 1, inplace = True)

이렇게 새로운 칼럼들이 생기고 1,0으로 True, False 표시를 하고있다.

이후 변수를 정리하고 마무리했다.

registered 와 casual 은 count 예측에는 필요없는 변수이기 때문에 날려주었다.

 

 

 

 

모델링

---

이제 준비가 끝났으니 모델링에 들어가보자.

from sklearn.model_selection import train_test_split
from sklearn import metrics
X_train = df_train.drop("count_Log", axis=1).values
target_label = df_train["count_Log"].values
X_test = df_test.values
X_tr, X_vld, y_tr, y_vld = train_test_split(X_train, target_label, test_size = 0.2, 
                                           random_state = 2000)

X_train 데이터에 count를 제외한 변수들을 넣어주고, target 에 count_Log(Count를 Log변환한 것) 을 넣어주자.

test_size는 0.2로 8:2의 비율이다. 

x_test는 test 데이터를 넣어주었다.

Gradient Boost를 사용하여 예측모델을 만들었다.

# Gradient Boost 사용
from sklearn.ensemble import GradientBoostingRegressor
regressor = GradientBoostingRegressor(n_estimators = 2000, learning_rate = 0.05,
                                     max_depth = 4,
                                     min_samples_leaf= 15, min_samples_split=10,
                                     random_state = 42)

regressor.fit(X_tr,y_tr)

y_hat = regressor.predict(X_tr)
plt.scatter(y_tr, y_hat, alpha = 0.2)
plt.xlabel('Targets(y_tr)', size= 18)
plt.ylabel("Predictions (y_hat)", size = 18)
plt.show()

오 잘 맞는 것 같다.

y_hat_test = regressor.predict(X_vld)
plt.scatter(y_vld, y_hat_test, alpha=0.2)
plt.xlabel("Targets(y_yld)", size = 18)
plt.ylabel("Predictions(y_hat_test)", size=18)
plt.show()

test데이터로도 눈으로 보기에는 잘 맞는 것 같다.

 

# 정확도 확인
from sklearn.metrics import mean_squared_log_error,mean_squared_error, r2_score,mean_absolute_error # for regression
from sklearn.metrics import accuracy_score,precision_score,recall_score,f1_score  # for classification

models=[GradientBoostingRegressor()]
model_names=['regressor']
rmsle=[]
d={}
for model in range (len(models)):
    clf=models[model]
    clf.fit(X_tr,y_tr)
    test_pred=clf.predict(X_vld)
    rmsle.append(np.sqrt(mean_squared_log_error(test_pred,y_vld)))
d={'Modelling Algo':model_names,'RMSLE':rmsle}   
d

{'Modelling Algo': ['regressor'], 'RMSLE': [0.1342902879875913]}

from sklearn.model_selection import cross_val_score 
accuracies = cross_val_score(estimator = regressor, X = X_tr, y = y_tr, cv = 8)

print(accuracies.mean()) 
print(accuracies.std())

0.95667106187069
0.0035836216554062743

96%정도의 정확성을 보인다.

 

use_logvals = 1

pred_xgb = regressor.predict(X_test)

sub_xgb = pd.DataFrame()
sub_xgb['datetime'] = datetime_test
sub_xgb['count'] = pred_xgb
if use_logvals == 1:
    sub_xgb['count'] = np.exp(sub_xgb['count'])
    
sub_xgb.to_csv("xgb.csv", index=False)

저장하고 제출하면 된다 ~~

 

 

도움받은 분들

---

해당 링크들을 필사하여 공부했습니다.

 

https://www.kaggle.com/code/kongnyooong/bike-sharing-demand-for-korean-beginners/notebook

[Bike Sharing Demand] for Korean Beginners (한글커널)

 

https://hong-yp-ml-records.tistory.com/category/KOREAN%20%EC%BA%90%EA%B8%80%20%ED%8A%9C%ED%86%A0%EB%A6%AC%EC%96%BC%20%ED%94%84%EB%A1%9C%EC%A0%9D%ED%8A%B8/Bike%20Sharing%20Demand

'KOREAN 캐글 튜토리얼 프로젝트/Bike Sharing Demand' 카테고리의 글 목록

 

 

https://www.kaggle.com/code/dogdriip/bike-sharing-demand/notebook

Bike Sharing Demand 선형/트리 기반 회귀 (파이썬 머신러닝 완벽 가이드)