뉴욕 택시수요예측

2021. 2. 16. 21:07

728x90

뉴욕 택시수요예측

2015년 5월의 yellow 뉴욕택시 데이트를 사용해서 수요량을 예측합니다. 시간은 30분 단위로 진행되었으며 주 평가척도는 mae입니다. 외부데이터는 날씨 데이터를 참고했으며 LSTM, Gru의 timestep은 메모리 크기에 맞게 유동적으로 변경하면서 진행했습니다.

Import

!pip install chart_studio

import chart_studio.plotly as py
import cufflinks as cf
import pandas as pd
import numpy as np
import seaborn as sns
import math
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

%%time
base_query="""
WITH base_data AS 
(
  SELECT nyc_taxi.*, gis.* EXCEPT (zip_code_geom)
  FROM (
    SELECT *
    FROM `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2015`
    WHERE 
        EXTRACT(MONTH from pickup_datetime) = 5
        and pickup_latitude  <= 90 and pickup_latitude >= -90
    ) AS nyc_taxi
  JOIN (
    SELECT zip_code, state_code, state_name, city, county, zip_code_geom
    FROM `bigquery-public-data.geo_us_boundaries.zip_codes`
    WHERE state_code='NY'
    ) AS gis 
  ON ST_CONTAINS(zip_code_geom, st_geogpoint(pickup_longitude, pickup_latitude))
)

SELECT 
    zip_code,
    DATETIME_TRUNC(pickup_datetime, hour) as pickup_hour,
    concat (format_datetime("%F %H:",pickup_datetime),
    case floor(extract(minute from pickup_datetime)/30) when 1.0 then 30 else 00 end) as minute, 
    EXTRACT(MONTH FROM pickup_datetime) AS month,
    EXTRACT(DAY FROM pickup_datetime) AS day,
    CAST(format_datetime('%u', pickup_datetime) AS INT64) -1 AS weekday,
    EXTRACT(HOUR FROM pickup_datetime) AS hour,
    CASE WHEN CAST(FORMAT_DATETIME('%u', pickup_datetime) AS INT64) IN (6, 7) THEN 1 ELSE 0 END AS is_weekend,
    round(sum(passenger_count)) AS passenger_cnt, #승객수
    round(sum(trip_distance)) AS trip_cnt,  #운행거리
    round(sum(total_amount)) AS total_amount_cnt, #승객에게 부과된 택시비
    COUNT(*) AS cnt	#콜 수
FROM base_data 
GROUP BY zip_code, pickup_hour, month, day, weekday, hour, is_weekend,minute
ORDER BY zip_code,pickup_hour

"""
base_df = pd.read_gbq(query=base_query, dialect='standard', project_id='abcd')

30분 단위로 나눕니다 (minute 변수로 지정)

EDA(Region)_퍼센트 별 trip_cnt

#전체 콜 대비 각 zip_code 비율
zipcode_sum=base_df.groupby(['zip_code'])[['cnt']].sum()
zipcode_sum

temp=(zipcode_sum/base_df[['cnt']].sum()).rename(columns={'cnt':'percent'})*100
temp.reset_index(inplce=True)
temp

temp['percent'].plot()

#log 값으로 변환
temp['log_cnt']=np.log10(temp['percent'])
plt.figure(figsize=(16,6));
sns.barplot(x='zip_code', y='log_cnt', data=temp);
sns.scatterplot(x='zip_code', y='log_cnt', data=temp);

log10(percent value) -> 0퍼센트 미만 값은 음수
log10(percent value) -> 0보다 큰 값은 양수
정리, 전체 비율에서 1%를 기준으로 분류

# 퍼센트 비율이 1%미만을 제거 하기 위한 필터
# 1%이상인 zip_code는 26개
filter_temp=temp[temp['log_cnt']>0]  
filter_temp.head()
# len(filter_temp)


zip_code	percent	log_cnt
0	10001	5.002529	0.699190
1	10002	1.962145	0.292731
2	10003	5.366801	0.729715
7	10009	2.003886	0.301873
8	10010	3.014356	0.479195

# 인덱스 초기화
#총 zip_code 383개
zipcode_sum.reset_index(inplace=True)
#zipcode_sum.head()
len(zipcode_sum)

383

# 1% 이상의 26개 지역 trip_cnt 보기
temp_cnt=[]
temp_zip_code=[]
for zip_code,cnt in np.array(zipcode_sum):
    for zip_code2 in np.array(filter_temp['zip_code']):
      if zip_code == zip_code2:
        temp_cnt.append(cnt)
        temp_zip_code.append(zip_code)
                              
temp_26=pd.concat([pd.DataFrame(temp_zip_code),pd.DataFrame(temp_cnt)], axis=1)
temp_26.columns=["zip_code",'cnt']
temp_26.tail()


zip_code	cnt
21	10075	224972
22	10119	157338
23	10128	371243
24	11371	324979
25	11430	282593

# 상위26개 지역 trip_cnt
plt.xticks(rotation =90)
sns.barplot(x="zip_code",y='cnt',data=temp_26)

EDA(Region)_zip_code 별 trip_cnt

#base_df에서 상위 26개 zip_code 정보 출력
filter_zip_df = base_df[base_df['zip_code'].isin(temp_zip_code)]
filter_zip_df.head()

#zip_code,date로 gropby 시킨 후 count 지표들을 비교
filter_zip_df['date']= filter_zip_df['pickup_hour'].dt.date
filter_zip_daily_df=filter_zip_df.groupby(['zip_code',"date"])[['cnt','passenger_cnt','trip_cnt','total_amount_cnt']].sum().reset_index()
filter_zip_daily_df

#히트맵
plt.figure(figsize=(16,10))
sns.heatmap(filter_zip_daily_df.pivot_table(values='cnt', index='date', columns='zip_code'), lw=.1, cmap='GnBu');

filter_zip_daily_df.pivot_table('cnt', index='date', columns='zip_code').iplot(kind='box', x='zip_code', y='cnt')

상위 26개 zip_code 중 10019가 가장 색이 짙은 걸로 보아 Trip 수가 많았음을 알 수 있다.
인덱스를 사용해서 수치로 보겠습니다

#콜수 최대값의 인덱스
filter_zip_daily_df['cnt'].idxmax()

63

filter_zip_daily_df.iloc[63,:]

zip_code                 10003
date                2015-05-02
cnt                      32042
passenger_cnt            55439
trip_cnt                 76000
total_amount_cnt        437920
Name: 63, dtype: object

5월 2일은 토요일이며, 이 날이 5월 달 중 가장 많은 이동을 보였습니다.
매 주말 토요일마다 이동 수가 급증하는지를 알아볼 필요가 있어 보입니다.
시간대 별 시각화를 통해 전반적인 흐름을 파악합니다.

EDA(Time)_시간

#시간대별 시각화
base_df['minute']=pd.to_datetime(base_df['minute'])
minute_group=base_df.groupby(["minute"])[["passenger_cnt","trip_cnt","total_amount_cnt","cnt"]].sum().reset_index()
minute_group=minute_group.set_index('minute')
minute_group.tail(5)

chart_studio.tools.set_credentials_file(username='name', api_key='abcde')
minute_group.iplot()

#위 링크 걸어둔 corab에서는 그래프가 없습니다(그래프 수 제한 초과)

크기가 맞지 않으므로 따로 분리해서 자세히 보겠습니다.

#운행거리
minute_group[['trip_cnt']].iplot()

운행거리를 나타내는 시각화
5월 2일(토요일)과 5월 31일(일요일), 5월 초, 말에 급증세를 보입니다.

#승객수, 콜 수
minute_group[['passenger_cnt',"cnt"]].iplot()

급증했던 5월 2일과 31일의 콜 수를 볼 때 특이점은 없어 보입니다.
경기 또는 축제처럼 이벤트가 발생한 날인지 알아보기 위해서 조사했지만 특별한 정보는 얻을 수 없었습니다.
가까운 지역 간의 이동이 아니라 장거리 손님이 많았음을 추론이 가능합니다
월 초 주말(2일)에 다른 지역으로 장거리 이동 후, 월 말(31일)에 다시 복귀하는 현상일 수도 있습니다

cnt 지표를 보면, 2일 대비 31일 콜수가 현저히 낮습니다
아마도 2일에 떠나간 사람들이 특정 지역으로 집중되었다면, 그 지역의 31일 골수는 급증세를 보일 거라 예상됩니다.

minute_group.reset_index(inplace=True)
sns.pairplot(minute_group, 
             diag_kind='kde',
             palette='bright') # pastel, bright, deep, muted, colorblind, dark
plt.show()

plt.figure(figsize=(16,10))
sns.heatmap(minute_group.corr(),annot=True)
plt.show()

trip_cnt(운행거리)를 제외하면 모두 강한 상관관계를 보입니다
상관관계를 갖는 변수는 삭제 대신 PCA를 활용해서 살려 보도록 했습니다

EDA(Time)_요일

#0:월요일 ~ 6:일요일
weekday_df=base_df.groupby(['weekday'])[['cnt','passenger_cnt','trip_cnt','total_amount_cnt']].sum()
weekday_df.reset_index(inplace=True)

                              #콜수, 손님 수
weekday_df.iplot(x='weekday',y=['cnt','passenger_cnt'],kind='line')

				# 운행 거리, 택시 비
weekday_df.iplot(x='weekday',y=['trip_cnt','total_amount_cnt'],kind='line')

weekday_df.iplot(x='weekday',kind='bar')

base_df.groupby(['hour','weekday']).mean()[['cnt']].unstack().iplot()

일요일을 제외한 날들은 저녁 7시에 콜수가 가장 많았음을 알 수 있습니다.
새벽 5시 기점으로 증가합니다.
주말만 별도로 자세히 살펴보겠습니다

EDA(Time)_주말

#0:평일,1:주말
weekend_df=base_df.groupby(['is_weekend'])[['cnt','passenger_cnt','trip_cnt','total_amount_cnt']].sum()
weekend_df.iplot(kind='bar')

base_df.groupby(['hour','is_weekend'])['cnt'].mean().unstack().iplot()

#시간대별 평일VS주말 비교
plt.figure(figsize=(16,10))
sns.heatmap(base_df.groupby(['hour','is_weekend'])['cnt'].mean().unstack().rename(columns={0:'weekday', 1:'weekend'}), lw=.5, annot=True, cmap='GnBu',
            fmt='g',annot_kws={'size':10});

평일은 러시아워 시간대에 콜 수가 집중되어 있습니다.
반면 주말은 비교적 고르게 분포되어 있습니다.

EDA(Time)_뉴욕의 전반적인 콜 수 흐름 보기

base_df2=base_df
base_df2['year']=2015
base_df2['date_except_hour']=pd.to_datetime(base_df2[['year','month','day']])
base_df2.groupby('date_except_hour').sum()[['cnt']].iplot()

5월 25일 월요일에 콜수가 급감을 보입니다.
구글링으로 공휴일, 이벤트 등 조사해봤지만 특이사항은 없었습니다.
현재까지 어떤 이유로 급감을 보였는지는 알 수가 없었습니다.

Model 순서

1. XGBoost

2. LightGBM

3. 스태킹 앙상블

4. LSTM

5. GRU

1. 데이터 전처리

Import

import numpy as np
from lightgbm import LGBMRegressor
from xgboost import XGBRegressor
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold
from sklearn.metrics import mean_absolute_error
import pandas as pd
import numpy as np
import seaborn as sns
import math, time
import matplotlib.pyplot as plt

1.1 target 분포 확인 및 log 적용

#target분포 확인
y_target= base_df['cnt']
y_target.hist()

# target log화(정규분포와 최대한 비슷하게 만들기 위함)
y_log=np.log1p(y_target)

#log_cnt 삽입
base_df['log_cnt']=y_log

y_log.hist()

1.2 날씨 변수 추가

#날씨 정보를 추가
#bad: 비가 내린 날, normal: 비가 내리지 않은 날
base_df["wheather"] = base_df["day"].apply(lambda x: 'bad' if x==16 or x==27 else "normal")

1.3 Scaling

from sklearn.preprocessing import StandardScaler

#객체 생성
scaler=StandardScaler() #기존 변수의 범위를 정규 분포로 변환

#강한 상관관계를 보인 passenger_cnt, total_amount_cnt는 차원축소로 변환할 예정
scaled_df=scaler.fit_transform(base_df[['passenger_cnt','trip_cnt','total_amount_cnt','cnt']])

#PCA 할 변수와 target 변수 분리
scaled_df=pd.DataFrame(scaled_df).iloc[:,0:3]

#기존 passenger_cnt, total_amount_cnt 변수 삭제
base_df.drop(['passenger_cnt','total_amount_cnt'], axis=1, inplace=True)

del scaled_df[1]  #trip_cnt 삭제
scaled_df=pd.DataFrame(scaled_df).rename(columns={0:"passenger_cnt",2:'total_amount_cnt'})
scaled_df

1.4 PCA

# 2개의 변수를 하나로 저차원 축소
from sklearn.decomposition import PCA
pca= PCA(n_components=1)

pca.fit(scaled_df)
pca_result= pca.transform(scaled_df)
print(pca_result,pca_result.shape)

[[ 3.87046796]
 [ 3.2337383 ]
 [ 1.69928273]
 ...
 [-0.70746083]
 [-0.70026516]
 [-0.7066679 ]] (158514, 1)

1.5 원핫인코딩

#원핫인코딩

#방법1
# get_dummies 
ohe_df=pd.get_dummies(base_df, columns=["zip_code",'day','hour','wheather'])
ohe_df.shape

(158514, 447)

#방법2
# #객체 생성
from sklearn.preprocessing import OneHotEncoder

# ohe= OneHotEncoder()
# ohe.fit(base_df[["zip_code",'day','hour','wheather']])
# ohe_output= ohe.transform(base_df[["zip_code",'day','hour','wheather']]).toarray()
# ohe_df=pd.concat([base_df, pd.DataFrame(ohe_output, columns=np.concatenate([ohe.categories_[0],ohe.categories_[1],ohe.categories_[2],ohe.categories_[3]], axis=0))],axis=1)
# ohe_df

#PCA 값과 합치기
df=pd.concat([ohe_df, pd.DataFrame(pca_result)], axis =1)
ohe_df=df.rename(columns={df.columns[-1]:"pca_value"})
ohe_df.tail()

1.6 Train/Test Dataset

def train_test_split(df,date):
  train_df=df[df['minute']<date]
  test_df=df[df['minute']>=date]
  return train_df,test_df

train_df,test_df=train_test_split(ohe_df, '2015-05-24')

#변수 삭제
del train_df['minute']
del test_df['minute']


# train/test 나누기
y_train_log=train_df.pop('log_cnt')
y_train_raw=train_df.pop('cnt')
x_train=train_df.copy()

y_test_log=test_df.pop('log_cnt')
y_test_raw=test_df.pop('cnt')
x_test=test_df.copy()

#시간 데이터 삭제
del x_train['pickup_hour']
del x_test['pickup_hour']

print('x_train shape:',x_train.shape, "y_train_log:",y_train_log.shape)

x_train shape: (117883, 447) y_train_log: (117883,)

XGBoost

def evaluate(y_true, y_pred):
  from sklearn.metrics import mean_absolute_error
  from sklearn.metrics import mean_squared_error
  y_true, y_pred = np.array(y_true), np.array(y_pred)
  mape= np.mean(np.abs((y_true-y_pred)/y_true))*100
  mae= mean_absolute_error(y_true,y_pred)
  mse= mean_squared_error(y_true,y_pred)
  score = pd.DataFrame([mape, mae, mse], index=['mape','mae','mse'], columns=['score'])
  return score

from xgboost import XGBRegressor
from sklearn.model_selection import GridSearchCV

xgb_reg= XGBRegressor()

#GridSearch
params={
    "max_depth":[5,7,10],
    "n_estimators":[100,300,500],
    "learning_rate":[0.001,0.01]
}
gridcv= GridSearchCV(xgb_reg, cv=2, param_grid= params, n_jobs=-1, verbose=True)
gridcv.fit(x_train,y_train_log)

print("GridSearchCV 최적 파라미터:", gridcv.best_params_)
print("GridSearchCV 최적 성능:", gridcv.best_score_)

GridSearchCV 최적 파라미터: {'learning_rate': 0.01, 'max_depth': 10, 'n_estimators': 500}
GridSearchCV 최적 성능: 0.9458285690503948

xgb_reg= XGBRegressor(n_estimators=500, random_state= 150, learning_rate=0.01, max_depth=10, objective='reg:squarederror')
xgb_reg.fit(x_train,y_train_log)
pred=xgb_reg.predict(x_test)
evaluate(y_test_log,pred)

plt.figure(figsize=(12,10))
plt.plot(pred[:200],  color='blue', label='pred')
plt.plot(y_test_log[:200].reset_index(drop=True), color='red', label='y_test')
plt.title('y_test VS pred')
plt.xlabel('index')
plt.ylabel('value')
plt.show()

LightGBM

from lightgbm import LGBMRegressor
from sklearn.model_selection import GridSearchCV

lgbm_reg= LGBMRegressor()

params={
    "max_depth":[5,7,10],
    "n_estimators":[100,300,500],
    "learning_rate":[0.001,0.01]
}
gridcv= GridSearchCV(lgbm_reg, cv=2, param_grid= params, n_jobs=-1, verbose=True)
gridcv.fit(x_train,y_train_log)

print("GridSearchCV 최적 파라미터:", gridcv.best_params_)
print("GridSearchCV 최적 성능:", gridcv.best_score_)

GridSearchCV 최적 파라미터: {'learning_rate': 0.01, 'max_depth': 7, 'n_estimators': 500}
GridSearchCV 최적 성능: 0.939007168663416

from lightgbm import LGBMRegressor
lgbm_reg= LGBMRegressor(n_estimators=500, random_state= 150, learning_rate=0.01, max_depth=7)
lgbm_reg.fit(x_train,y_train_log)
pred=lgbm_reg.predict(x_test)
score=evaluate(y_test_log, pred)
score

plt.figure(figsize=(12,10))
plt.plot(pred[:200],  color='blue', label='pred')
plt.plot(y_test_log[:200].reset_index(drop=True), color='red', label='y_test')
plt.title('y_test VS pred')
plt.xlabel('index')
plt.ylabel('value')
plt.legend()
plt.show()

두 모델 중 XGBoost가 좀 더 좋은 성능을 보입니다.
하지만 학습시간은 131.8 min 소요되었습니다. 반면 LightGBM은 1.2 min이 걸렸습니다
PCA 변수가 성능에 악영향을 주는지 알아보기 위해 전 후 차이를 보겠습니다.

LightGBM_PCA 삭제

#pca 삭제
del x_train['pca_value']
del x_test['pca_value']

lgbm_reg= LGBMRegressor(n_estimators=500, random_state= 150, learning_rate=0.01, max_depth=7)
lgbm_reg.fit(x_train,y_train_log)
pred2=lgbm_reg.predict(x_test)
score2=evaluate(y_test_log, pred2)
score2

plt.figure(figsize=(12,10))
plt.plot(pred2[:100], color='blue',label='pred')
plt.plot(y_test_log[:100].reset_index(drop=True), color='red', label='y_test')
plt.title('y_test VS pred')
plt.xlabel('index')
plt.ylabel('value')
plt.legend()
plt.show()

pca 전후 차이 오류 값에 차이를 보입니다.
pca 삭제 후 다른 모델들을 실행해 보겠습니다.

스태킹 앙상블(zip_code 7개 샘플링)

전체 데이터로는 Corab에서 실행 불가능하므로 데이터 수가 많은 상위 7개의 zip_code를 대상으로 분석함

EDA 시각화 중에서 heatmap 참고하여 데이터가 많은 zip_code을 골랐습니다.

개별 Model 설정&평가

#개별 ML 모델 생성
svr_reg= SVR()
rf_reg= RandomForestRegressor(n_estimators=100,random_state=0)
dt_reg= DecisionTreeRegressor()
ada_reg= AdaBoostRegressor(n_estimators=100)

#스태킹으로 만들어진 데이터세트를 학습,예측할 최종모델
lgbm_reg= LGBMRegressor(n_estimators=300, random_state= 150, learning_rate=0.01, max_depth=5)

#개별 모델 학습
svr_reg.fit(x_train,y_train_log)
rf_reg.fit(x_train,y_train_log)
dt_reg.fit(x_train,y_train_log)
ada_reg.fit(x_train,y_train_log)

# 학습된 개별 모델들이 각자 반환하는 예측 데이터 세트를 생성하고 개별모델의 정확도 측정
svr_pred=svr_reg.predict(x_test)
rf_pred=rf_reg.predict(x_test)
dt_pred=dt_reg.predict(x_test)
ada_pred=ada_reg.predict(x_test)

#평가
svr_score=evaluate(y_test_log, svr_pred)
rf_score=evaluate(y_test_log, rf_pred)
dt_score=evaluate(y_test_log, dt_pred)
ada_score=evaluate(y_test_log, ada_pred)
print(pd.concat([svr_score,rf_score,dt_score,ada_score],axis=1))

#각 개별 분류기의 예측값은 1차원 형태이므로 이걸 행 형태로 붙인 뒤 transpose를 취하여 데이터 세트로 만든다

#4행 형태로 붙인다
pred= np.array([svr_pred, rf_pred,dt_pred,ada_pred])
print(pred.shape)

# transpose로 행과 열의 위치 교환하여 각 컬럼끼리 위치가 맞도록 해준다
final_pred=np.transpose(pred)
print(final_pred.shape)

(4, 2688)
(2688, 4)

스태킹 최종 모델 결과

# 최종 메타 모델인 LightGBM를 실시
lgbm_reg.fit(final_pred, y_test_log)
final=lgbm_reg.predict(final_pred)

evaluate(y_test_log,final)

가장 낮은 에러 성능 수치를 보입니다.
이번에는 CV를 적용하여 스태킹을 재 실행해 보겠습니다.

CV 기반의 스태킹

def get_stakong_base_datasets(model, x_train, y_train, x_test, n_folds):
  from sklearn.model_selection import KFold
  from sklearn.metrics import mean_absolute_error

  kf = KFold(n_splits=n_folds, shuffle=False, random_state=0)

  train_fold_pred=np.zeros((x_train.shape[0],1)) #x_train에서 n_folds수 만큼 학습,검증폴드 분리됨, 즉 검증세트=x_train.shape[0] 될때까지 진행
  test_pred= np.zeros((x_test.shape[0],n_folds)) #x_test의 예측값이 n_folds수 만큼 생성됨 
  print(model.__class__.__name__,'model start')

  for folder_counter, (train_idx, valid_idx) in enumerate(kf.split(x_train)):
    print('\n폴드세트:', folder_counter,'학습/검증폴드 생성')
    train_x= x_train[train_idx] #x_train의 학습폴드
    train_y= y_train[train_idx] #y_train의 학습폴드
    test_x= x_train[valid_idx]  #x_train의 검증폴드

    #학습폴드 인덱드에 해당하는 x_train, y_train를 가지고 fit 진행
    model.fit(train_x,train_y)

    #학습한 model를 검증폴드 예측(최종모델의 학습 데이터로 사용)
    train_fold_pred[valid_idx,:]=model.predict(test_x).reshape(-1,1)

    #학습한 model를 원본 테스트 데이터인 x_test 예측
    test_pred[:,folder_counter]=model.predict(x_test)
  
  #x_test를 예측한 값을 열 기준으로 평균내고 세로로 쌓는다(최종모델의 test 데이터로 사용)
  test_mean=np.mean(test_pred, axis=1).reshape(-1,1)

  return train_fold_pred,test_mean

개별 Model 설정&평가

#개별 ML 모델 생성
svr_reg= SVR()
rf_reg= RandomForestRegressor(n_estimators=100,random_state=0)
dt_reg= DecisionTreeRegressor()
ada_reg= AdaBoostRegressor(n_estimators=100)

#스태킹으로 만들어진 데이터세트를 학습,예측할 최종모델
lgbm_reg= LGBMRegressor(n_estimators=300, random_state= 150, learning_rate=0.01, max_depth=5)

#ndarray로 변경(안하면 get_stakong_base_datasets 적용 x)
x_train=np.array(x_train)
y_train_log=np.array(y_train_log)
x_test=np.array(x_test)

svr_train,svr_test=get_stakong_base_datasets(svr_reg,x_train,y_train_log,x_test,3)
rf_reg_train,rf_reg_test=get_stakong_base_datasets(rf_reg,x_train,y_train_log,x_test,3)
dt_reg_train,dt_reg_test=get_stakong_base_datasets(dt_reg,x_train,y_train_log,x_test,3)
ada_train,ada_test=get_stakong_base_datasets(ada_reg,x_train,y_train_log,x_test,3)

stack_final_x_train=np.concatenate((svr_train,rf_reg_train,dt_reg_train,ada_train),axis=1)
stack_final_x_test=np.concatenate((svr_test,rf_reg_test,dt_reg_test,ada_test), axis=1)
print('원본 학습 피처 데이터 shape:', x_train.shape, '원본 테스트 피처 shape:',x_test.shape)
print('스태킹 학습 피처 데이터 shape:', stack_final_x_train.shape, '스태킹 테스트 피처 데이터 shape:',stack_final_x_test.shape)

원본 학습 피처 데이터 shape: (7728, 68) 원본 테스트 피처 shape: (2688, 68)
스태킹 학습 피처 데이터 shape: (7728, 4) 스태킹 테스트 피처 데이터 shape: (2688, 4)

lgbm_reg.fit(stack_final_x_train,y_train_log) #원본 학습 label과 fit
stack_final=lgbm_reg.predict(stack_final_x_test)
evaluate(y_test_log,stack_final)

CV 적용하기 전에 비해 성능이 안 좋아졌습니다.
CV를 적용한다고 무조건의 성능개선은 이뤄지진 않습니다.
LSTM, GRU를 활용해서 딥러닝에서의 성능은 어떨지 실행해 보겠습니다.

딥러닝

Import

import tensorflow as tf
from tensorflow.keras.layers import Input, LSTM, GRU, Dense,Dropout
from tensorflow.keras.models import Model, Sequential
from tensorflow.keras import layers
from tensorflow.keras.optimizers import SGD, Adam, RMSprop
from tensorflow.keras import models
from tensorflow.keras.wrappers.scikit_learn import KerasRegressor
from sklearn.model_selection import GridSearchCV
from sklearn.datasets import make_regression

import seaborn as sns
import numpy as np
from datetime import datetime 

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler

zip_code 7개 중 5개로 줄여서 진행했습니다(메모리 부족)

TimeStep 설정

# ex) 12시~12시30분의 minute 값은 0, 12시 30분~13시의 minute 값은 30 
ohe_df['minute']=pd.to_datetime(ohe_df['minute'])
ohe_df['minute2']=ohe_df['minute'].dt.minute

# 시간대별 정렬
ohe_df=ohe_df.sort_values(by=['hour','minute2'], axis=0)
# ohe_df.to_excel('imsi.xlsx')
ohe_df.head()

timestep은 시간대 별로 정렬했습니다. 4시간 단위로 묶었습니다.
ex) 12시~4시까지의 데이터로 4시~4시 30분의 값을 예측합니다.
각 시간대별 몇 개의 데이터가 있는지 아래 코드를 통해 알아봅니다.

num=0
total_=[]
for idx in range(np.array(ohe_df).shape[0]):
  try:
    if ohe_df.iloc[idx]['minute2'] == ohe_df.iloc[idx+1]['minute2'] and ohe_df.iloc[idx]['hour'] == ohe_df.iloc[idx+1]['hour']:
      num+=1
    else:
      print('zip_code는{0} 이며 {1}의 총 갯수는 {2}입니다'.format(ohe_df.iloc[idx]['zip_code'],ohe_df.iloc[idx]['minute'], num+1))
      total_.append(num+1)  # 1를 더해준 이유: 처음 시작 인덱스는 0부터이다. 코드에서는 1부터 시작한다. 
      num=0  #카운트 초기화
  except:

    try:
      if ohe_df.iloc[idx]['minute2'] == 30 and ohe_df.iloc[idx]['hour'] == 23:  #23시30분~00시 처리부분
        print('zip_code는{0} 이며 {1}의 총 갯수는 {2}입니다'.format(ohe_df.iloc[idx]['zip_code'],ohe_df.iloc[idx]['minute'], num+1))
        total_.append(num)
      else:
        break
    except:
      break
      
      
      
zip_code는10019 이며 2015-05-31 00:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 00:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 01:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 01:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 02:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 02:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 03:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 03:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 04:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 04:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 05:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 05:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 06:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 06:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 07:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 07:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 08:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 08:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 09:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 09:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 10:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 10:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 11:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 11:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 12:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 12:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 13:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 13:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 14:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 14:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 15:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 15:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 16:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 16:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 17:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 17:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 18:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 18:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 19:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 19:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 20:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 20:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 21:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 21:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 22:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 22:30:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 23:00:00의 총 갯수는 155입니다
zip_code는10019 이며 2015-05-31 23:30:00의 총 갯수는 155입니다

같은 시간대를 묶은 결과 각 시간대 별 155개가 있음을 확인했습니다.
4시간을 기준으로 Timestep을 설정해야 하므로 155*8=1240 길이가 됩니다.

#변수 삭제
ohe_df.drop(['pickup_hour','month',"zip_code",'minute','minute2',"wheather",'hour','day'], axis=1, inplace=True)

#feature/label 정의
label_raw=ohe_df['cnt']
label_raw=np.array(label_raw).reshape(-1,1)

label_log=ohe_df['log_cnt']
label_log=np.array(label_log).reshape(-1,1)

input_data=ohe_df.drop(['cnt','log_cnt'],axis=1)
input_data=np.array(input_data)
print("input_data:",input_data.shape, 'label_data_y_log:',label_log.shape)

input_data: (7440, 65) label_data_y_log: (7440, 1)

Scaling

from sklearn.preprocessing import MinMaxScaler

# label_data scaling
label_scaler = MinMaxScaler()
label_raw = label_scaler.fit_transform(label_raw)

다양한 시도를 위해 이번에는 cnt값을 MinMaxScaler를 사용하겠습니다.

Train, Test Dataset

T=620 #2시간 과거 데이터
N=len(input_data)-T 
D=input_data.shape[1] 

#데이터 넣을 공간 만들기
x=np.zeros((N,T,D))
y_log=np.zeros(N)
y_raw=np.zeros(N)

for t in range(N):
  x[t,:,:]= input_data[t:t+T]
  y_log[t]= label_log[t+T]
  y_raw[t]= label_raw[t+T]

print(x.shape, y_log.shape)

(6789, 651, 65) (6789,)

원래 계획은 4시간 길이였지만 메모리 부족으로 인해 2시간 데이터로 변경되었습니다.

from sklearn.model_selection import train_test_split

x_train,x_test,y_train,y_test= train_test_split(x,y_raw,test_size=0.2 ,shuffle=False ,random_state=5)
print("x_train:", x_train.shape, "y_train:",y_train.shape)

x_train: (5431, 651, 65) y_train: (5431,)

신경망 모델

신경망 모델에서는 LSTM과 Gru를 주로 사용했습니다. 성능 개선과 과적합 제거 목적으로 각 모델별로 조금씩 다르게 구성했습니다. 머신러닝과 직관적인 비교를 위해 성능지표는 mae로 통일해서 진행했습니다.

mae (Mean absolute error)-실제 값과 측정(예측) 값과의 차이

n = 오차의 개수

∑ = 합을 나타내는 기호

|xi−x| = 절대 오차

mae는 실제값에서 예측값을 빼준 절대 오차값을 n만큼 나눈 값이기 때문에 낮을수록 좋은 성능을 나타냅니다.

기본 모델 설정

비교를 위해 기본 모델로 실행했습니다.
머신러닝과의 직관적인 성능 차이를 보기 위해서 성능지표를 mae로 통일

model = Sequential()
model.add(layers.Flatten(input_shape=(T,input_data.shape[1])))  #(8,58)
model.add(layers.Dense(32, activation='relu'))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
 history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    validation_data=(x_test, y_test))

#그래프 시각화
plt.plot(history.history['mae'], label='mae')
plt.plot(history.history['val_mae'],label='val_mae')
plt.legend()
plt.show()

LSTM

model = Sequential()
model.add(layers.LSTM(32,input_shape=(T,input_data.shape[1])))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) 
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

LSTM(층 추가)

model = Sequential()
model.add(layers.LSTM(32,return_sequences=True,input_shape=(T,input_data.shape[1])))
model.add(layers.LSTM(64))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) # loss감소가 10번 이상 미발생이면 stop
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

LSTM(층 추가+Dropout)

model = Sequential()
model.add(layers.LSTM(32,return_sequences=True,input_shape=(T,input_data.shape[1])))
model.add(Dropout(0.3))
model.add(layers.LSTM(64))
model.add(Dropout(0.3))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) # loss감소가 10번 이상 미발생이면 stop
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

층을 추가한 LSTM 모델에서 val_mae:0.0771로 가장 좋은 성능을 보였고
층을 추가한 상태에서 Dropout을 적용한 결과 같은 epoch대비 성능을 개선시키지 못했습니다
patience=15로 짧게 설정했기 때문에 개선 과정 중 멈췄을 가능성도 있습니다

GRU

# Gru
model = Sequential()
model.add(layers.GRU(32,input_shape=(T,input_data.shape[1])))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) # loss감소가 10번 이상 미발생이면 stop
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

GRU(층 추가)

# Gru
model = Sequential()
model.add(layers.GRU(32,return_sequences=True,input_shape=(T,input_data.shape[1])))
model.add(layers.GRU(64))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) # loss감소가 10번 이상 미발생이면 stop
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

GRU(층 추가+Dropout)

# Gru
from tensorflow.keras.layers import Dropout

model = Sequential()
model.add(layers.GRU(32,return_sequences=True,input_shape=(T,input_data.shape[1])))
model.add(Dropout(0.3))
model.add(layers.GRU(64))
model.add(Dropout(0.3))
model.add(layers.Dense(1))
model.compile(optimizer=Adam(), loss = 'mse', metrics=['mae'])
earlystopper = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=15, verbose=1) # loss감소가 10번 이상 미발생이면 stop
history= model.fit(x_train,y_train,
                    epochs=100,
                    batch_size=32,
                    callbacks=[earlystopper],
                    validation_data=(x_test, y_test))

Gru 모델도 층 추가한 모델이 val_loss:0.0863으로 가장 좋은 성능을 보였습니다.
추가로 Dropout을 한 결과 오히려 과적합 현상이 의심되는 그래프를 보입니다.

Lightgbm 스태킹 mae : 0.095
LSTM의 val_mae: 0.0771

728x90

'projects' 카테고리의 다른 글

따릉이 자전거 이용 예측_dacon (0)	2021.09.30
Face Mask Detection (0)	2021.08.25
동서발전 태양광 발전량 예측 프로젝트_dacon (0)	2021.06.29
국민청원글 토픽모델링 (0)	2021.01.18
CNN_BTS 멤버 image classification (0)	2021.01.10

H_record