서울시 따릉이 데이터 분석 및 예측하기

지난 2021년, 빅데이터를 융합전공하면서 했던 첫 미니프로젝트를 인턴십에서 발표하게 됐다.

솔직히 2년이나 지난 일이라 잘 기억나지 않기도 하고...

지금보니 이상치 제거와 연속 & 범주 데이터가 같이 있는데 정규화도 되어있지 않은 걸 보니 너무 끔찍(?)해서

정리도 할겸, 다시 EDA 해보기로 결정했다.

Problem

예측: Regression

종속 변수 : 따릉이 대여 수(count, integer)

사용된 예측 변수

Hour : 시간, 0부터 23까지의 integer
날씨 변수 : temperature, precipitation(mm) -> float
공휴일 여부(0/1)
1시간 전 이용건수
2시간 전 이용건수
3시간 전 이용건수

측정 척도: Root Mean Squared Error(RMSE)

Preivew

'서울시 공공 데이터’에서 제공하는 자전거 대여 장소와 시간대별 자전거 대여 수(數) 자료를 분석하여 ‘뚝섬유원지역 1번 출구 앞’에 위치한 자전거 대여소의 대여 건수를 예측한다.

해당 대여소의 자전거 대여 수를 예측하여 이용자들의 편의와 자전거의 이용률을 높여 이윤을 극대화 하기 위해 본 프로젝트를 진행하게 되었음.

Data Load

자전거 이용건수 데이터는 2020.07 ~ 2021.06 기간, 총 12개월을 사용했다.

아무래도 예측하기에는 최소 1년 정도의 데이터는 필요하다고 판단했기 때문이다.

또한, 이 데이터에는 서울시의 자전거 대여소 전체를 제공하기 때문에 '뚝섬 유원지역 1번 출구'를 기준으로

예측하기로 했다.

df = pd.concat(df_frame, ignore_index = True)
df = df[df.대여소번호 == 502].drop(['대여소번호','대여소명','대여구분코드','성별','연령대코드','운동량','탄소량','이동거리','사용시간'],axis=1)

# column 전체의 type을 바꿔줌!
df = df.astype({'이용건수':'float'})

그 외 다른 컬럼은 필요가 없으므로 제거했다.

Temper data concat

강수량 데이터가 필요했다. 데이터는 기상청에서 확보할 수 있었다.

특히 온도 데이터는 결측치가 한 곳이 비어있어 평균 온도로 수기 입력했다.

강수량 데이터는 비가 오지 않은 날은 Nan이여서, 그 또한 0.0으로 변경했다.

df = df.groupby(['대여일자','대여시간']).sum().reset_index()

w = pd.read_csv('weather.csv',encoding = "cp949").drop(['지점','지점명'],axis=1)

# 결측치 처리
new_data = {
    '일시' : '2020-08-26 13:00',
    '기온(°C)' :33.7,
    '강수량(mm)':0.0
}

idx = 1357 # 추가하려는 곳의 index 번호
 
temp1 = w[w.index < idx]
temp2 = w[w.index >= idx]

w = temp1.append(new_data,ignore_index=True).append(temp2, ignore_index=True)


w1 = w.iloc[:,:1]['일시'].str.split(' ',expand=True) #일시
w1.columns = ['대여일자','대여시간']
w2 = w.iloc[:,1:] #기온, 강수량

a = w1['대여시간'].str.split(':',expand=True).iloc[:,:1].astype('int64')
a.columns = ['대여시간']

b = pd.concat( [w1.drop(['대여시간'],axis=1), a], axis=1 )
c = pd.concat([b,w2],axis=1)

total = pd.merge(df,c,how='right') #자전거, 날씨 합친데이터
#total.tail(50)

Data Shift

1/2/3시간 전 자전거 이용건수를 변수로 채택했기 때문에, 자전거 이용건수를 1시간 씩 shift했다.

이때 생기는 6개의 결측치는 직접 이전 데이터를 확인하고 수기로 입력했다.

total[['이용건수']] = total[['이용건수']].fillna(0)
total[['강수량(mm)']] = total[['강수량(mm)']].fillna(0)

total['이용건수1'] = total[['이용건수']].shift(periods=1)
total['이용건수2'] = total[['이용건수']].shift(periods=2)
total['이용건수3'] = total[['이용건수']].shift(periods=3)

total['강수량(mm)1'] = total[['강수량(mm)']].shift(periods=1)

공휴일 여부

우선, 데이터가 12개월 뿐이라 공휴일 자체도 그리 많지 않았다. 따라서 직접 공휴일 리스트를 만들고 이를 반복해 공휴일은 1, 공휴일이 아닌 것은 0으로 처리했다.

...처럼 하려고 했는데, 예전 우리 팀원의 노력이 무색하게도 날짜를 쪼갤 수 있는 아주 좋은 방법이 있었다.

temp = total
temp["year"] = temp["대여일자"].dt.year
temp["month"] = temp["대여일자"].dt.month
temp["day"] = temp["대여일자"].dt.day
temp["weekday"] = temp["대여일자"].dt.weekday

최종 데이터

이렇게 최종 데이터가 완성됐다. 이를 heatmap으로 확인하면,

The number of uses(이용건수) 라인을 보면 이전 시간의 사용 건수에 대해 가장 상관관계가 깊다고 알 수 있다.

아무래도 21년도의 여름 데이터는 빠져서 그런지 preciptation(강수량)은 음의 상관관계를 보이게 된다. 조금 아쉽긴하다.

아마 온도보다 가장 영향을 많이 받지 않았을까.

그밖에도 온도와 대여 시간 또한 0.3, 0.4 정도의 상관관계를 보이고 있다.

또한, 시간 대비 이용 건수를 확인했을 때 17~18시 사이가 가장 급증하는 형태로 보이고 있다.

아무래도 퇴근 시간이라 그런지... 이것을 어떻게 활용하면 좋을까?

이상치

프로젝트 당시 우리 팀원이 그려줬던 plot이다.

대충 봐도 알 수 있듯이, 저 멀리 200대가 넘어가는 데이터가 포착된다.

그 땐 이상치 탐지에 대한 이론만 배웠지 활용할 생각을 못했다. ~~왜 그랬을까? 과거의 나..~~

이용 건수를 보면 평균은 13대인 것에 반해 max값은 213대다.

뭐 사실 그래프를 보면 알 수 있겠지만... 좀 커다란 이상치는 두어개 정도다.

plt.boxplot(final_df["The number of uses"])

아까 max값의 200대의 영향인지 동 떨어져 있는 데이터가 몇몇 개 보인다.

from scipy import stats
import numpy as np 

# sns.distplot(final_df['The number of uses'])
# sns.distplot(final_df['The number of uses'])
trainWithoutOutliers = final_df[np.abs(final_df['The number of uses'] - final_df['The number of uses'].mean()) <= (3 * final_df['The number of uses'].std())]

#print(final_df.shape)
#print(trainWithoutOutliers.shape)

plt.boxplot(trainWithoutOutliers["The number of uses"])

IQR의 값에 1.5를 곱해서 최대, 최소 값 너머의 값들을 이상치로 판단하고 제거했다.

그리고 나서 다시 boxplot을 보니,

나름 예쁘게 된 것 같다. 데이터는 8760개에서 7932로 줄었다.

로그변환

기계학습을 진행할 때 정확도를 위해 정규분포를 따르는 데이터를 활용하는 게 좋다.

그런데 대부분의 데이터들은 그게 지켜지진 않는 것 같다.

실제로 밀집도를 확인해보니 왼쪽으로 치우쳐져 있는 것이 보인다.

이걸 어떻게 해결해야하나, 했는데 로그 변환을 통해 어느 정도는 해결할 수 있었다.

sns.distplot(final_df['The number of uses'])

그래서 로그 변환을 진행했다. (더 자세한 설명은 https://suppppppp.github.io/posts/Why-Series-MDM-1/)

대강 설명하자면, x(0,1)일때는 기울기가 매우 가파르지만, y의 경우는 -무한대부터 0까지 매우 범위가 넓기 때문에 로그로 변환해주는 것.

그리고 log1p를 사용하는 이유는 0의 경우 로그로 변환하면 -inf와 같은 매우 작은 음수값이 되어버린다.

df_log = temp
df_log = np.log1p(df_log["이용건수"])

완벽하진 않지만 어느정도 쏠림 현상을 해결할 수 있었다. 이용건수와 1/2/3시간 전 이용건수에 적용해서 로그 변환 또한 진행했다.

이렇게!

어느 정도 이제 분석과 회귀 예측을 할 준비가 끝난 것 같다. 그리고 대여일자의 경우 year, month, day, weekday로 쪼갰으니 필요없다. 바로 drop!

분석

import seaborn as sns 
plt.rc("font", family = "Malgun Gothic")
sns.set(font="Malgun Gothic", 
rc={"axes.unicode_minus":False}, style='white')

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

sns.barplot(data=temp, x="year", y="이용건수", ax=ax1)
sns.barplot(data=temp, x="month", y="이용건수", ax=ax2)
sns.barplot(data=temp, x="day", y="이용건수", ax=ax3)
sns.barplot(data=temp, x="weekday", y="이용건수", hue = "weekday", ax=ax4,)

ax1.set(ylabel='대여량',xlabel='년도', title="연도별 대여량",)
ax2.set(xlabel='월',ylabel='대여량',title="월별 대여량")
ax3.set(xlabel='일',ylabel='대여량',title="일별 대여량")
ax4.set(xlabel='요일',ylabel='대여량', title="요일별 대여량")

연도/월/일/요일 별로 대여량을 시각화해봤다. 아무래도 휴일에 자전거를 많이 빌렸고, 상대적으로 추운 1월과 12월은 대여량이 적었으며 5월에 가장 많이 대여한 걸 알 수 있다.

연도별 차이는 크게 없는 게 당연할지도 모르겠다.

상관관계에서 아주 큰 영향을 미쳤던 3가지를 산포도로 나타내보았다.

분포를 봤을 때, 1시간 전 이용건수가 대여 건수에 가장 큰 영향을 미쳤음을 알 수 있다.

당연히 강수량은 0일 때 대여량이 가장 많으며, 온도가 높을 때 대여량이 훨씬 많다.

예측

드디어 예측이다.

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error,mean_absolute_error
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge

#df = pd.read_csv('./bike_and_weather_team_drpm.csv')

y_target = temp["이용건수"]
def get_data():
    X_train, X_test, y_train, y_test = train_test_split(x, y_target, test_size=0.3, random_state=40)
    
    return [X_train, X_test, y_train, y_test]


X_train, X_test, y_train, y_test = get_data()

clf = LinearRegression()
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)

r2 = r2_score(y_test, y_pred)
print('r2 score is = ' + str(r2))


mse = mean_squared_error(y_test, y_pred)
print('mse is = ' + str(mse))

mae = mean_absolute_error(y_test, y_pred)
print('mae is = ' + str(mse))

mse = mean_squared_error(y_test, y_pred)
print('rmse is = ' + str(np.sqrt(mse)))

x는 이용건수를 제외한 모든 컬럼, y_target은 이용건수이다.