https://colab.research.google.com/drive/1CsnK5YrMDtSiAeQezPmw_7YPjZlg_8uS?hl=ko 

dacon의 따릉이 자전거 이용 예측 프로젝트입니다. 

현재 기준으로 데이콘 리더보드에서 6위를 기록하고 있습니다.

감사합니다 :)

https://colab.research.google.com/drive/1eesj-6EBRPpZn6AebXDIOhdULzh7x_BU?hl=ko 

Mask Detection을 수행해 봤습니다. dataset은 캐글에서 받았습니다.

아래 그림은 유튜브 영상 중 일부를 편집하여 모델로 돌린 결과로써 threshold 0.1로 설정했음을 알려드립니다. 

dataset이 좀 더 많고 다양하다면 지금의 성능보다 훨씬 좋아 질것 같습니다. 자세한 내용은 위 링크에서 보실수 있습니다.

본 내용은 Dacon의 동서발전 태양광 발전량 예측 AI 경진대회에 참가한 프로젝트 내용입니다.

주어진 데이터를 통해 태양광 발전량 예측 모델을 만들어 봤습니다.

아래는 Data& Target Data 일부분을  캡처한 그림입니다. 

1. Import and Libraries

!pip install tsfresh

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
from tqdm import tqdm

# Ignore the warnings
import warnings 
warnings.filterwarnings('ignore') # 경고 뜨지 않게 설정

# System related and data input controls
import os

# Data manipulation and visualization
import pandas as pd
pd.options.display.float_format = '{:,.2f}'.format
pd.options.display.max_rows = 10 
pd.options.display.max_columns = 20
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import preprocessing

# Modeling algorithms
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import RFECV
from sklearn import linear_model
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm
from scipy import stats
import tsfresh

# Model selection
from sklearn.model_selection import train_test_split

# Evaluation metrics
from sklearn.metrics import mean_squared_log_error, mean_squared_error, r2_score, mean_absolute_error

from google.colab import drive

2. 선형보건법 fcst/ obs 비교 시각화

데이콘에서 제시한 선형보건법 함수와 주어진 데이터인 예보데이터(fcst),실측치(obs)를 활용하여 가장 오차가 작은 예보 시간대를 선택합니다.

for number in [11,14,17,20]:
    inter_obs = obs_to_interpolate(dangjin_obs,False)
    inter_fcst= fcst_to_interpolate(dangjin_fcst,number,False)
    inter_obs2=inter_obs[24:].reset_index(drop=True)
    
    pd.concat([inter_obs2[['Temperature']],
         inter_fcst[0:25608][['Temperature']]], axis=1).plot(kind='line',figsize=(20,10),
                                                   linewidth=3, fontsize=20,
                                                   ylim=(-15,40))
    print('%s일 경우의 그래프'%number)
    plt.title('Temparature plot', fontsize=20)
    plt.xlabel('label', fontsize=15)
    plt.ylabel('Temparature', fontsize=15)
    plt.show()
    
    #잔차
    Resid=((inter_obs2[['Temperature']].values.flatten()-inter_fcst[0:25608][['Temperature']].values.flatten())**2).mean()
    print('Resid',Resid)

*11시 예보 데이터 사용 그래프 => Resid 2.899278434194671

*14시 예보 데이터 사용 그래프=>Resid 2.899264549644704

*17시 예보 데이터 사용 그래프=>Resid 2.81235024211184

*20시 예보 데이터 사용 그래프=>Resid 2.811606550904232 (20시 예보 데이터가 실측치와 가장 근접)

3. Exploratory Data Analysis

3.1 연도별

for col in concat_df.columns[:4]:
    plt.figure(figsize=(15,6))
    g = sns.boxplot(x='year', y=col, data=concat_df, showfliers=False)
    g.set_title('box plot of %s'%col, size=20)
    g.set_xticklabels(g.get_xticklabels(), rotation=90)
    plt.show()

연도별 시각화 결과, dangjin_floating 수치만 감소세를 보입니다. 수치로 확인하기 위해 아래처럼 작업을 진행했습니다.

concat_df.groupby('year').agg({'dangjin_floating':['count','sum','mean']})

2018년도는 3월부터 기록되어 있으므로 Count 수가 다릅니다. 하지만 태양 에너지 평균은 가장 높습니다.

dangjin_floating 만 감소세를 보이는 원인에 대해서는 파악이 되지 않으나 dangjin_floating 에너지 수급 패턴을 짐작할수 있습니다.

3.2 월별

Y_colname = 'dangjin_floating','dangjin_warehouse','dangjin','ulsan'
for col in Y_colname:
    concat_df2 = concat_df[concat_df[col]!=0] 
    result = concat_df2.groupby('month').agg({col:['sum','mean','count']}).reset_index(drop=False)
    # print(result.iloc[:,3])
    plt.figure(figsize=(15,6))
    plt.subplot(1,2,1)
    plt.plot(result.iloc[:,0], result.iloc[:,2])
    plt.title('%s mean'%col)
    plt.subplot(1,2,2)
    plt.plot(result.iloc[:,0], result.iloc[:,1])
    plt.title('%s sum'%col)
    plt.show()

concat_df.boxplot(column ='ulsan', by ='month', grid =True, figsize=(12,5))
concat_df.boxplot(column ='dangjin_floating', by ='month', grid =True, figsize=(12,5))

여름에 가장 많은 태양광 에너지를 얻을거라는 예상과 달리 3~4월 기점으로 에너지 수급이 7월까지 급격한 감소 추세를 보입니다. 7월 이후로 반등하는 모습을 보이다가 10월 이후, 겨울이 본격적으로 시작되면서 다시 감소 추세입니다. 가장 수급이 좋은 시기는 초봄인 3월~4월과 더위가 점차 사그라드는 9~10월로 보입니다. 즉 적당한 환경일때 효율적인 수급이 가능 한걸로 추측됩니다.

이후 조사 결과 태양광 에너지 수급에 있어서 중요한 핵심은 태양광의 모듈 온도 이라고 합니다. 이 모듈 온도가 높아지면 효율이 떨어진다고 합니다. 그래서 실제로 이 온도를 낮추기 위한 시설이나 장치를 설치하는 경우가 있다고 합니다. 

3.3 시간대별

for col in Y_colname:
  concat_df2 = concat_df[concat_df[col]!=0]     
  result=concat_df2.groupby(['hour']).agg({col:['count','sum','mean']})
  result.iloc[:,2].plot(label = '%s_energy'%col, legend=True, figsize=(20,10))

시간대는 약 13시30분에 가장 높은 수급을 보입니다.

정리를 하면, 3~4월의 13시30분에 가장 높은 에너지를 확보 할수 있습니다.

데이터 보기전에는 무조건 온도가 높으면 좋을 줄 알았는데 조금 의외의 결과였습니다.

3.4 날씨별

*Cloud : 하늘상태(1-맑음, 2-구름보통, 3-구름많음, 4-흐림)

concat_df['Cloud'] = concat_df['Cloud'].astype('int') #일기예보의 구름 양 예측에 따른 에너지 양
Y_colname = ['dangjin_floating','dangjin_warehouse','dangjin','ulsan']
cloud = concat_df.groupby(['Cloud'])[Y_colname].sum()
cloud.plot(kind = 'bar')

대체로 구름이 없을수록 태양열 에너지 수급이 원활함을 알수 있음 특이한 점은 구름이 많은 상태인 3에서 소폭 상승을 보였다는 점이다. 자세한 파악을 위해서 상관성을 시각화 해보았습니다.

import seaborn as sns
plt.figure(figsize=(16,10))
sns.heatmap(concat_df.corr(), annot=True)
plt.show()

Corr 결과만 봤을때는 음의 상관관계이므로 1->4로 갈수록 수급이 떨어져야 정상입니다.  추측으로는 다중공선성의 영향이 있을 거라 생각됩니다. 

3.5 일자별

#Date 설정
concat_df['date'] = concat_df['Forecast_time'].dt.date
concat_df['date'] = pd.to_datetime(concat_df['date'])

#2019/01일 부터 시각화
date_df = concat_df[concat_df['date'] >= '2019-01-01']
date_df =date_df.groupby('date')[Y_colname].sum().reset_index()
plt.figure(figsize=(20,9))
plt.plot('date', 'dangjin_floating', 'g-',label='dangjin_floating',data=date_df)
plt.gcf().autofmt_xdate()
plt.legend(loc='center left', bbox_to_anchor=(1,0.5), fontsize=20)
plt.title('date_df')
plt.show()

월별 시각화와 일맥상통으로 봄에서 겨울 사이에는 감소를 보이다가 다시 상승하는 패턴을 보입니다.

상승이나 하강의 추세는 아직까진 없어 보이는 대신 계절성 패턴이 명확하게 보입니다. 

4. Data Preprocessing & Feature Engineering

# 선형 보건법
def to_date(x):
    return pd.DateOffset(hours=x)

def fcst_to_interpolate(dangjin_fcst,date_number,graph=True):

    #날짜타입으로 변경
    dangjin_fcst['Forecast_time'] = pd.to_datetime(dangjin_fcst['Forecast time'])

    #적용할 date_number
    fcst = dangjin_fcst[dangjin_fcst['Forecast_time'].dt.hour==date_number]
    fcst = fcst[(fcst['forecast']>=24-date_number) & (fcst['forecast']<=47-date_number)] 

    #forecast 시간 합치기
    fcst['Forecast_time'] = fcst['Forecast_time'] + fcst['forecast'].map(to_date)

    #사용될 컬럼 추출
    fcst = fcst[['Forecast_time', 'Temperature', 'Humidity', 'WindSpeed', 'WindDirection', 'Cloud']]
    
    #선형보건법이 적용된 데이터 프레임 생성 
    fcst_= pd.DataFrame()
    fcst_['Forecast_time'] = pd.date_range(start='2018-03-02 00:00:00', end='2021-03-01 23:00:00', freq='H')
    fcst_merge = pd.merge(fcst_, fcst, on='Forecast_time', how='outer')
    inter_fcst = fcst_merge.interpolate()

    if graph:
        plt.figure(figsize=(20,5))
        days = 5
        plt.plot(inter_fcst.loc[:24*days, 'Forecast_time'], inter_fcst.loc[:24*days, 'Temperature'], '.-')
        plt.plot(fcst_merge.loc[:24*days, 'Forecast_time'], fcst_merge.loc[:24*days, 'Temperature'], 'o')
    
    return inter_fcst

def obs_to_interpolate(dangjin_obs,graph=True):

    #날짜타입으로 변경
    dangjin_obs['Forecast_time'] = pd.to_datetime(dangjin_obs['time'])

    #사용될 컬럼 추출
    obs = dangjin_obs[['Forecast_time', 'Temperature', 'Humidity', 'WindSpeed', 'WindDirection', 'Cloud']]
    
    #선형보건법이 적용된 데이터 프레임 생성 
    obs_= pd.DataFrame()
    obs_['Forecast_time'] = pd.date_range(start='2018-03-01 00:00:00', end='2021-01-31 23:00:00', freq='H')
    obs_merge = pd.merge(obs_, obs, on='Forecast_time', how='outer')
    inter_obs = obs_merge.interpolate()

    if graph:
        plt.figure(figsize=(20,5))
        days = 5
        plt.plot(inter_obs.loc[:24*days, 'Forecast_time'], inter_obs.loc[:24*days, 'Temperature'], '.-')
        plt.plot(obs_merge.loc[:24*days, 'Forecast_time'], obs_merge.loc[:24*days, 'Temperature'], 'o')
    
    return inter_obs

def feature_engineering(concat_df, target):
    if 'Forecast_time' in concat_df.columns:
        concat_df['Forecast_time']=pd.to_datetime(concat_df['Forecast_time'])
        concat_df.index =concat_df['Forecast_time']

    concat_df=concat_df.asfreq('H')

    #0으로 되어 있는 값을 작년 수치 & 변수 대체
    concat_df.loc['2020-06-19 11:00:00':"2020-06-25 23:00:00",'dangjin_floating'] = concat_df.loc['2019-06-19 11:00:00':'2019-06-25 23:00:00','dangjin_floating'].values
    concat_df.loc['2020-06-19 11:00:00':"2020-06-25 23:00:00",'Temperature':]  = concat_df.loc['2020-06-19 11:00:00':"2020-06-25 23:00:00",'Temperature':].values 

    #nan값 채우기
    concat_df.loc['2020-06-26 01:00:00':"2020-06-27 00:00:00",'dangjin_floating'] = concat_df.loc['2020-06-28 01:00:00':'2020-06-29 00:00:00','dangjin_floating'].values
    concat_df.loc['2018-03-17 01:00:00':"2018-03-19 00:00:00",'dangjin_warehouse'] = concat_df.loc['2018-03-15 01:00:00':'2018-03-17 00:00:00','dangjin_warehouse'].values
    
    #날짜 형식 추가
    concat_df['hour']=concat_df['Forecast_time'].dt.hour
    concat_df['year']=concat_df['Forecast_time'].dt.year
    concat_df['month']=concat_df['Forecast_time'].dt.month
    concat_df['day']=concat_df['Forecast_time'].dt.day

    result =sm.tsa.seasonal_decompose(concat_df[target], model='additive')
    Y_seasonal = pd.DataFrame(result.seasonal)
    Y_seasonal.fillna(method='bfill', inplace=True)
    Y_seasonal.fillna(method='ffill', inplace=True)
    Y_seasonal.columns= ['%s_seasonal'%target]

    concat_df = pd.concat([concat_df,Y_seasonal], axis=1)
    return concat_df

def feature_engineering_lag(concat_fe,x_colname,number,multi_number=1):
  for col in x_colname:
    for mul_number in range(1,multi_number+1):
      shift_number = concat_fe[[col]].shift(number*mul_number)
      shift_number.fillna(method='bfill', inplace=True)
      shift_number.fillna(method='ffill', inplace=True)
      shift_number.columns =['{}_lag{}'.format(col,number*mul_number)]
      concat_fe = pd.concat([concat_fe,shift_number], axis=1)
  return concat_fe
  
def feature_engineering_diff(concat_fe, x_colname, diffs):
  diff_multi_number =2
  for col in x_colname:
    for number in range(1,diffs+1):
      diff_number = concat_fe[[col]].diff(number)
      diff_number.fillna(method='bfill', inplace=True)
      diff_number.fillna(method='ffill', inplace=True)
      diff_number.columns =['{}_diff{}'.format(col,number)]
      concat_fe =pd.concat([concat_fe, diff_number], axis=1)
  return concat_fe

def feature_engineering_duplicated(x_train,x_val,target):
  target = ['%s_seasonal'%target]  # '%s_trend'%target '%s_diff'%target,'%s_lag_1'%target, '%s_lag_2'%target,'Y_%s_week'%target, 'Y_%s_day'%target
  for col in target:
    x_val.loc['2021-01-01 00:00:00':'2021-01-31 23:00:00',col] = x_train.loc['2020-01-01 00:00:00':'2020-01-31 23:00:00',col].values
  return x_val
  
#Datasplit
def train_dataset(concat_df,criteria, X_colname, target):
    if 'time' in concat_df.columns:
      del concat_df['time']
    if 'Forecast_time' in concat_df.columns:
      del concat_df['Forecast_time']

    train_df = concat_df.loc[concat_df.index<criteria,:]
    val_df = concat_df.loc[concat_df.index>=criteria,:]
    
    x_train = train_df[X_colname]
    y_train = train_df[target]
    
    x_val = val_df[X_colname]
    y_val = val_df[target]
    print('x_train:',x_train.shape, 'y_train:',y_train.shape)
    print('x_val',x_val.shape,'y_val',y_val.shape)
    return x_train,x_val, y_train,y_val

#Scaling
def Scaling(x_train,x_val,scaler):
  from sklearn import preprocessing
  scaler = scaler
  scaler_fit = scaler.fit(x_train)
  x_train_feRS = pd.DataFrame(scaler_fit.transform(x_train), index =x_train.index, columns=x_train.columns)
  x_val_feRS = pd.DataFrame(scaler_fit.transform(x_val),index =x_val.index, columns=x_val.columns)
  return x_train_feRS,x_val_feRS

#다중공선성
def feature_engineering_VIF(x_train):
  vif= pd.DataFrame()
  vif['VIF_Factor'] = [variance_inflation_factor(x_train.values,i) for i in range(x_train.shape[1])]
  vif['Feature'] = x_train.columns
  # x_colname_vif = vif.sort_values(by='VIF_Factor', ascending=True)['Feature'][:num_variables].values
  return  vif

#evalution
def evaluation(y_train, pred, graph=True):
    loss_length = len(y_train.values.flatten()) - len(pred)
    if loss_length !=0:
        print('길이가 맞지 않습니다.')
    if graph == True:
        pd.concat([y_train, pd.DataFrame(pred, index=y_train.index, columns=['prediction'])], axis=1).plot(kind='line', figsize=(20,6),
                                                                                                        xlim=(y_train.index.min(), y_train.index.max()),
                                                                                                        linewidth=3, fontsize=20)
        plt.title('Time Series of Target', fontsize=20)
        plt.xlabel('index', fontsize=15)
        plt.ylabel('Target value', fontsize=15)
    Residual = pd.DataFrame(y_train.values.flatten() - pred, index= y_train.index, columns=['Error'])
    return Residual

def evaluation_trte(y_train, y_train_pred, y_val, y_val_pred, graph=True ):
    Residual_train = evaluation(y_train, y_train_pred,graph=graph)
    Residual_val = evaluation(y_val, y_val_pred, graph=graph)
    return Residual_train,Residual_val

def nmae_10(y_pred, dataset,capacity=1000):
    y_true = dataset.get_label()
    absolute_error = abs(y_true - y_pred)
    absolute_error /= capacity
    target_idx = np.where(y_true>=capacity*0.1)
    nmae = 100 * absolute_error[target_idx].mean()
    return 'score', nmae, False

def sola_nmae(answer, pred,capacity=1000):
    absolute_error = np.abs(answer - pred)
    absolute_error /= capacity
    target_idx = np.where(answer>=capacity*0.1)
    nmae = 100 * absolute_error[target_idx].mean()
    return nmae
    
#feature engineering
concat_fe=feature_engineering(concat_df,"dangjin_floating")

#X's diff 생성 
x_colname=['Temperature','WindSpeed','Humidity'] # 'year', 'month', 'day'  
concat_fe = feature_engineering_diff(concat_fe,x_colname,6) 

#X's lag 생성 
# x_colname2=['WindDirection']
# concat_fe = feature_engineering_lag(concat_fe,x_colname2,4)

#datasplit
Y_colname = ['dangjin_floating','dangjin_warehouse','dangjin','ulsan']
X_colname =[col for col in concat_fe.columns if col not in Y_colname+['time','Forecast_time',"year"]]
x_train_feR,x_val_feR, y_train_feR,y_val_feR =train_dataset(concat_fe,'2021-01-01 00:00:00',X_colname,'dangjin_floating')

# #feature_engineering_duplicated
x_val_feR = feature_engineering_duplicated(x_train_feR,x_val_feR,'dangjin_floating')

#설명변수 스케일링
x_train_feRS, x_val_feRS = Scaling(x_train_feR,x_val_feR,preprocessing.StandardScaler())

#다중공선성
vif = feature_engineering_VIF(x_train_feRS)
vif.sort_values(by=['VIF_Factor'], inplace=True)
vif.reset_index(drop=True, inplace=True)
VIF_colname =[]
for idx,x in enumerate(vif['VIF_Factor'].values): 
  if x <= 15 : 
    VIF_colname.append(vif['Feature'][idx])

#REV
ols = linear_model.LinearRegression()
refecv = RFECV(estimator=ols, step=1, scoring='neg_mean_squared_error')
refecv.fit(x_train_feRS, y_train_feR)
refecv.transform(x_train_feRS)
refecv_colname =[colname for idx,colname in enumerate(x_train_feRS.columns) if refecv.ranking_[idx]==1]

5. modeling

*target의 정규화를 각 모델에 적용한 결과, 비 정규화가 더 높은 성능을 보였습니다. 정규화 결과는 간결한 진행을 위해서 제외 시켰습니다. 

5.1 auto_arima

!pip install pmdarima

from pmdarima.arima import auto_arima

arima_model = auto_arima(y_train_feR, start_p=2, d=0, start_q=2,
                         max_p=3, max_d=1, max_q=3,
                          start_P=1, D=1, start_Q=1,
                         max_P=7, max_D=3, max_Q=7, m=24, seasonal=True,
                         error_action='ignore', trace=True,
                         supress_warnings=True, stepwise=True,
                         random_state=20) #, n_fits=50
arima_model.fit(y_train_feR)

auto_arima로 파리미터를 찾으려고 했으나 computing 문제로 초반부터 중단이 되었습니다. SARIMAX의 acf,pacf 결과로 order 추정을 해보았습니다. 

5.2 SARIMAX

위 설정값은 2020년도 값만 적용했을때의 결과값입니다. 이 order 값을 전체 데이터에 적용해 보았지만 computing 문제로 실행되지 않았습니다. 따라서 order 값을 적용했을때의 모습을 보여 드리기 위하여 아래 임의의 order값을 적용했습니다.

#SARIMAX

trend_order = (0,0,0)
seasonal_order = (2,1,2,24)
model = sm.tsa.SARIMAX(y_train_feR, trend='c', order =trend_order, seasonal_order= seasonal_order, exog = x_train_feRS) 
result = model.fit()
print('SARIMAX {}*{} and AIC: {}'.format(trend_order,seasonal_order,result.aic))
display(result.summary())

#예측
SARIMAX_train_pred = result.predict(start=0, end=len(y_train_feR)-1, exog = x_train_feRS)
SARIMAX_val_pred = result.forecast(len(y_val_feR), exog = x_val_feRS) #1월달 예측

display('train_CV Score : ', sola_nmae(y_train_feR.to_numpy(), SARIMAX_train_pred.to_numpy()))
display('val_CV Score : ', sola_nmae(y_val_feR.to_numpy(), SARIMAX_val_pred.to_numpy()))

#시각화
plt.figure(figsize=(20,10))
plt.subplot(2,1,1)
plt.plot(SARIMAX_train_pred, label = 'train_pred')
plt.plot(y_train_feR, label ='true')
plt.legend()
plt.title('train plot')

plt.subplot(2,1,2)
plt.plot(SARIMAX_val_pred, label = 'val_pred')
plt.plot(y_val_feR, label ='true')
plt.legend()
plt.title('val plot')
plt.show()

#잔차진단
result.plot_diagnostics(figsize=(10,8))
plt.tight_layout()
plt.show()

res= result.resid
sm.tsa.graphics.plot_acf(res, lags=100)
sm.tsa.graphics.plot_pacf(res, lags=100)

분석결과 

자기상관성 존재하며 정규분포는 아니며 등분산은 아니라는 결과를 받았습니다. 정상성을 만들지 못했습니다. 총 62번의 실험을 해보았지만 정상성을 만족하는 order값은 차지 못했습니다. 또한 데이콘에서 제시한 베이스 코드의 점수인 8.3에 못미치는 점수가 나왔으므로 SARIMAX로 모델링 하는건 시간적으로 무리가 있다고 생각했습니다. 

5.3 지수평활

CV_Score=[]
numbers=range(1,100)
# log = True
for number in numbers:
  fit_exp = ExponentialSmoothing(y_train_feR, seasonal_periods=24*number, trend=None, seasonal='additive').fit(use_boxcox=False)
  forecast_exp = fit_exp.forecast(24*31)

  #예측점수 계산
  CV_Score.append(sola_nmae(y_val_feR.to_numpy(),forecast_exp.to_numpy()))
  idxs.append(number)
pd.options.display.max_rows=100
pd.DataFrame(CV_Score, index= numbers, columns=['Score'])

아래는 가장 점수 성능이 좋은 number=12 일 때의 시각화입니다

seasonal_periods의 number 값을 100까지 실행한 결과 9.5 점수를 얻은 number=12일때 가장 좋은 성능을 보였습니다. 

*extract features

개선된 모델을 위하여 중간에 tsfresh 라이브러리를 사용하여 변수를 늘려보았습니다.

!pip install tsfresh

import tsfresh
# extracted_features=tsfresh.extract_features(concat_x, column_id='time')

#데이터 불러오기
extracted_features = pd.read_csv('/content/drive/MyDrive/extracted_features.csv', sep=',')

#rename
extracted_features.rename(columns = {'Unnamed: 0' : 'Forecast_time'}, inplace = True)
extracted_features.set_index('Forecast_time',inplace=True)

extracted_features_df= extracted_features.copy()
extracted_features_df.dropna(axis=1, inplace=True) #결측값 제거
extracted_features_df.isnull().sum().sum()

#열의 평균이 0 or 1이면 삭제
del_list=[]
columns = extracted_features_df.columns
for columns in columns:
  mean_calc = extracted_features_df[columns].mean()
  if mean_calc == 0.0 or mean_calc==1.0:
    del_list.append(columns)
extracted_features_df.drop(columns=del_list, axis=1,inplace=True

#datasplit
from sklearn.model_selection import train_test_split
target = concat_fe['dangjin_floating']
x_train_feR,x_val_feR, y_train_feR,y_val_feR = train_test_split(extracted_features_df, 
                                                    target, 
                                                    test_size=0.1, 
                                                    shuffle=False, 
                                                    random_state=1004)

# 다중공선성
vif = feature_engineering_VIF(x_train_feR)
vif.sort_values(by=['VIF_Factor'], inplace=True)
vif.reset_index(drop=True, inplace=True)
VIF_colname =[]
for idx,x in enumerate(vif['VIF_Factor'].values): 
  if x <= 10: 
    VIF_colname.append(vif['Feature'][idx])
extracted_features_df[VIF_colname]

총 생성된 변수 1300여개 중 전처리 후 다중공선성을 제거하기 위해 vif를 실행한 결과 입니다. 

VIF_colname

위 9개 변수를 포함하여 가장 성능이 좋았던 random forest에 적용 해본 결과, 성능 다운이 발생했습니다. 변수가 많아져서 과적합이 발생한 걸로 추측해 볼수 있습니다. 아래에 진행 될 모델은 위 9개 변수를 제외하여 진행한 내용입니다.

5.4 randomForeest and LGB

def model(target,fcst,energy_df, n_components=3, y_log =False):

  #선형보건법
  #20시간으로 한 선형보건법 적용
  target_fcst= fcst_to_interpolate(fcst,20,False)
  target_fcst = target_fcst.loc[1:25607].reset_index(drop=True)

  #energy + dangjin_fcst
  energy=energy_df[24:-1].reset_index(drop=True)
  concat_df = pd.concat([energy,target_fcst], axis=1)
  concat_df.index = concat_df['Forecast_time']

  #feature engineering
  concat_fe=feature_engineering(concat_df,target)

  #X's diff 생성 
  x_colname=['Temperature','WindSpeed','Humidity'] # 'year', 'month', 'day'  
  concat_fe = feature_engineering_diff(concat_fe,x_colname,6) 

  #X's lag 생성 
  x_colname2=['day']
  concat_fe = feature_engineering_lag(concat_fe,x_colname2,24,2)

  #datasplit
  Y_colname = ['dangjin_floating','dangjin_warehouse','dangjin','ulsan']
  X_colname =[col for col in concat_fe.columns if col not in Y_colname+['time','Forecast_time',"year"]]
  x_train_feR,x_val_feR, y_train_feR,y_val_feR =train_dataset(concat_fe,'2021-01-01 00:00:00',X_colname,target)

  # #feature_engineering_duplicated
  x_val_feR = feature_engineering_duplicated(x_train_feR,x_val_feR,target)

  #설명변수 스케일링
  x_train_feRS, x_val_feRS = Scaling(x_train_feR,x_val_feR,preprocessing.StandardScaler())
  # 1. Normal, 2.Robust = 3.Standard  4. Minmax

  # target log -> 예측범위를 좁히는데 용이
  if y_log:
    y_train_feR = np.log1p(y_train_feR).copy()
    y_val_feR = np.log1p(y_val_feR).copy()

  #다중공선성
  vif = feature_engineering_VIF(x_train_feRS)
  vif.sort_values(by=['VIF_Factor'], inplace=True)
  vif.reset_index(drop=True, inplace=True)
  VIF_colname =[]
  for idx,x in enumerate(vif['VIF_Factor'].values): 
    if x <= 15 : 
      VIF_colname.append(vif['Feature'][idx])

  #PCA 변환
  pca_col =[col for col in x_train_feRS.columns if col not in VIF_colname]
  pca = PCA(n_components=n_components)

  #pca fit
  pca.fit(x_train_feRS[pca_col])
  train_pca_colname = pca.transform(x_train_feRS[pca_col])
  val_pca_colname= pca.transform(x_val_feRS[pca_col])
  train_pca_df = pd.DataFrame(train_pca_colname, columns=['pca_components%s'%x for x in range(n_components)]) #['pca_components1','pca_components2','pca_components3']
  val_pca_df = pd.DataFrame(val_pca_colname, columns=['pca_components%s'%x for x in range(n_components)])
  x_train_feRS = x_train_feRS[VIF_colname]
  x_val_feRS = x_val_feRS[VIF_colname]
  
  for x in range(n_components):
    x_train_feRS.loc[:,'pca_components{}'.format(x)]=train_pca_df.iloc[:,x].values
    x_val_feRS.loc[:,'pca_components{}'.format(x)]=val_pca_df.iloc[:,x].values

  #REV
  ols = linear_model.LinearRegression()
  refecv = RFECV(estimator=ols, step=1, scoring='neg_mean_squared_error')
  refecv.fit(x_train_feRS, y_train_feR)
  refecv.transform(x_train_feRS)
  refecv_colname =[colname for idx,colname in enumerate(x_train_feRS.columns) if refecv.ranking_[idx]==1]

  # 1.RandomForest
  # fit_ranfor = RandomForestRegressor(n_estimators=600, random_state=123,min_samples_leaf=6, min_samples_split=14).fit(x_train_feRS, y_train_feR)
  fit_ranfor = RandomForestRegressor(n_estimators=600, random_state=123).fit(x_train_feRS, y_train_feR)
  random_train_pred = fit_ranfor.predict(x_train_feRS)
  random_val_pred = fit_ranfor.predict(x_val_feRS)

  # 2.LGB
  import lightgbm as lgb
  params = {
      'learning_rate': 0.01,
      'objective': 'regression',
      'metric':'mae',
      'seed':42,
  }

  train_datasets = lgb.Dataset(x_train_feRS.to_numpy(), y_train_feR.to_numpy())
  val_dataset = lgb.Dataset(x_val_feRS.to_numpy(), y_val_feR.to_numpy())
  fit_lgb = lgb.train(params, train_datasets, 10000, val_dataset, feval=nmae_10, verbose_eval=500, early_stopping_rounds=200)

  #lgb 예측
  lgb_train_pred = fit_lgb.predict(x_train_feRS)
  lgb_val_pred = fit_lgb.predict(x_val_feRS)

  #y_log 값 환원 
  if y_log:
    y_train_feR = np.expm1(y_train_feR).copy()
    y_val_feR = np.expm1(y_val_feR).copy()
    random_train_pred = np.expm1(random_train_pred).copy()
    random_val_pred = np.expm1(random_val_pred).copy()
    lgb_train_pred = np.expm1(lgb_train_pred).copy()
    lgb_val_pred = np.expm1(lgb_val_pred).copy()

  # 잔차 진단
  rf_Residual_train,rf_Residual_val = evaluation_trte(y_train_feR,random_train_pred, y_val_feR, random_val_pred, graph=True)
  display('rf_CV Score : ', sola_nmae(y_val_feR.to_numpy(), random_val_pred))
  rf_score = sola_nmae(y_val_feR.to_numpy(), random_val_pred)

  lgb_Residual_train,lgb_Residual_val = evaluation_trte(y_train_feR,lgb_train_pred, y_val_feR, lgb_val_pred, graph=True)
  display('lgb_CV Score : ', sola_nmae(y_val_feR.to_numpy(), lgb_val_pred))
  lgb_score = sola_nmae(y_val_feR.to_numpy(), lgb_val_pred)
  
  return rf_score,lgb_score,rf_Residual_train,lgb_Residual_train,fit_ranfor,fit_lgb,concat_fe,x_train_feRS,y_train_feR,x_val_feRS,y_val_feR,VIF_colname
  
rf_score,lgb_score,rf_Residual_train,lgb_Residual_train,fit_ranfor,fit_lgb,concat_fe,x_train_feRS,y_train_feR,x_val_feRS,y_val_feR,VIF_colname = model('dangjin_floating',dangjin_fcst,energy)

# 랜덤 포레스트 중요도 
importances_values = fit_ranfor.feature_importances_
importances_df = pd.Series(importances_values, index = x_train_feRS.columns).sort_values(ascending=False)
plt.figure(figsize=(20,9))
plt.title('Importancsvariables')
sns.barplot(x = importances_df, y = importances_df.index)
plt.show()

sns.distplot(rf_Residual_train['Error'].iloc[1:], norm_hist='True', fit=stats.norm)
figure, axes = plt.subplots(1, 4, figsize=(30,5))
pd.plotting.lag_plot(rf_Residual_train['Error'].iloc[1:], lag=1, ax=axes[0])
pd.plotting.lag_plot(rf_Residual_train['Error'].iloc[1:], lag=5, ax=axes[1])
pd.plotting.lag_plot(rf_Residual_train['Error'].iloc[1:], lag=10, ax=axes[2])
pd.plotting.lag_plot(rf_Residual_train['Error'].iloc[1:], lag=50, ax=axes[3])

figure, axes = plt.subplots(2,1,figsize=(12,5))
figure = sm.graphics.tsa.plot_acf(rf_Residual_train['Error'].iloc[1:], lags=100, use_vlines=True, ax=axes[0])
figure = sm.graphics.tsa.plot_pacf(rf_Residual_train['Error'].iloc[1:], lags=100, use_vlines=True, ax=axes[1])

rf(랜덤포레스트)의 정규분포/자기상관성 상태를 시각화 해보았습니다. 

ACF를 보면 초반에 튀는 값이 약 7개가 존재하며 중간에 24주기를 중심으로 살짝 튀는 값들도 존재합니다.

제가 이번 프로젝트를 하면서 보았던 자기상관성 그래프들 중 가장 안정적인 모습을 보입니다.

하지만 튀는 값은 존재하므로 제가 잡지 못한 패턴은 존재합니다. 

6. Test

def test_dataset(target,concat_fe,start='2021-02-01 00:00:00', end = '2021-02-28 23:00:00', duplicated_start='2020-02-01 00:00:00', duplicated_end = '2020-02-28 23:00:00',n_components=3):
  from sklearn import preprocessing
  from sklearn.decomposition import PCA

  # data loading
  dangjin_fcst = pd.read_csv('/content/drive/MyDrive/dangjin_fcst_data.csv', sep=',')
  ulsan_fcst = pd.read_csv('/content/drive/MyDrive/ulsan_fcst_data.csv', sep=',')

  #당진
  ##선형보건법
  dangjin_fcst= fcst_to_interpolate(dangjin_fcst,20,False)
  dangjin_start_index =dangjin_fcst[dangjin_fcst['Forecast_time']==start].index[0]
  dangjin_end_index =dangjin_fcst[dangjin_fcst['Forecast_time']==end].index[0]
  dangjin_fcst_test = dangjin_fcst.iloc[dangjin_start_index:dangjin_end_index+1,:]

  #울산
  ##선형보건법
  ulsan_fcst= fcst_to_interpolate(ulsan_fcst,20,False)
  ulsan_start_index =ulsan_fcst[ulsan_fcst['Forecast_time']==start].index[0]
  ulsan_end_index =dangjin_fcst[ulsan_fcst['Forecast_time']==end].index[0]
  ulsan_fcst_test = dangjin_fcst.iloc[ulsan_start_index:ulsan_end_index,:]

  ##feature_engineering
  def test_feature_engineering(fcst_test,target):
      if 'Forecast_time' in fcst_test.columns:
          fcst_test['Forecast_time']=pd.to_datetime(fcst_test['Forecast_time'])
          fcst_test.set_index(fcst_test['Forecast_time'], drop=True, inplace=True)

      fcst_test=fcst_test.asfreq('H')

      #날짜 형식 추가
      fcst_test['hour']=fcst_test['Forecast_time'].dt.hour
      fcst_test['month']=fcst_test['Forecast_time'].dt.month
      fcst_test['day']=fcst_test['Forecast_time'].dt.day
      del fcst_test['Forecast_time']

      ##X's diff 생성 
      x_colname=['Temperature','WindSpeed','Humidity'] # 'year', 'month', 'day'  
      test_fe = feature_engineering_diff(fcst_test,x_colname,6) 

      #X's lag 생성 
      x_colname2=['day']
      test_fe = feature_engineering_lag(test_fe,x_colname2,24,2)

      #seasonal_duplicated
      target2 = ['%s_seasonal'%target]
      test_fe.loc[:,target2] = concat_fe.loc[duplicated_start:duplicated_end,target2].values

      #스케일링
      scaler = preprocessing.StandardScaler()
      test_feRS = pd.DataFrame(scaler.fit_transform(test_fe), index = test_fe.index, columns=test_fe.columns)
      
      #PCA 변환
      
      
      pca_col =[col for col in test_feRS.columns if col not in VIF_colname]
      pca = PCA(n_components=n_components)

      #pca fit
      pca.fit(test_feRS[pca_col])
      test_pca_colname = pca.transform(test_feRS[pca_col])
      test_pca_df = pd.DataFrame(test_pca_colname, columns=['pca_components%s'%x for x in range(n_components)])
      test_feRS = test_feRS[VIF_colname]
      for x in range(n_components):
        test_feRS.loc[:,'pca_components{}'.format(x)]=test_pca_df.iloc[:,x].values
      return test_feRS

  #전처리 작업
  if target =='ulsan':
    test_feRS = test_feature_engineering(ulsan_fcst_test,target)
  else:
    test_feRS = test_feature_engineering(dangjin_fcst_test,target)

  return test_feRS
test_feRS = test_dataset('dangjin_floating',concat_fe)
test_feRS.head()

# 2월 test 예측
rf_pred = fit_ranfor.predict(test_feRS)
pd.DataFrame(rf_pred,index = test_feRS.index, columns =['pred']).plot()

* dangjin_floating

* dangjin_warehouse 

* dangjin_warehouse 

* ulsan

울산 1월 예측 점수는 다른 target과는 달리 큰 폭으로 오차가 줄어들었습니다. lgb가 randomforest보다 더 좋은 성능을 보이고 있습니다. 하지만 일관성을 위해서 randomforest로 2월달을 예측하였습니다.

*프로젝트를 끝내며

데이콘 첫 참가라서 꽤나 어리둥절 했던 것같습니다 프로젝트를 거의 다 만들고

본격적인 평가를 할려고 보니까 6월10일부터 누적 평가라고 하더라구요 이를 확인했을때는 이미 

6월 말쯤이라서.. 뒤늦게 제출해도 지나간 날에는 0점으로 매긴다고 합니다. 그래서 프라이빗 평가 대신 대회 끝난 후 퍼블릭 평가를 해볼 생각입니다. 

상관성을 보면 추세부분에도 미세하게 연관성이 있어 보입니다. 시계열 모델로 하면 ARMA로 추세부분을 어떻게 해볼순있겠지만 가장 좋은 성능을 보인 랜덤에서는 추세를 어떻게 해야 할까 고민을 좀 해봤습니다.

그래서 생각해 본게 시계열 분해에서 나온 추세를 그냥 넣지 말고 차분이나 lag 등 전처리를 하면 조금이라도 잡을수 있지 않을까 생각해 보면서 마무리를 지어 봅니다.

이 대회 첫참가 하고 데이터를 딱 받아 봤을때 뭘 먼저 해야하는지 감도 잡히지 않아서 일주일동안 답답해 했던 날들이 생각납니다. 이번 대회 덕분에 자신감도 생기고 뿌듯했던 경험이었습니다

뉴욕 택시수요예측

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

목차

1. EDA

1. 1 EDA(Region)

1.1.1 퍼센트 별 Trip_cnt

1.1.2 zip_code별 Trip_cnt

1.2 EDA(Time)

1.2.1 시간

           1.2.2 요일

           1.2.3 주말

2. 데이터전처리

    2.1 target 분포 확인 / log 적용

    2.2 날씨변수 추가

    2.3 Scaling

    2.4 PCA 생성

    2.5 원핫인코딩

    Train/Test Dataset

3. 모델 구축과 검증

    3.1 XGBoost

           3.1.1 GridSearch

           3.1.2 최적 파라미터 적용

    3.2 LightGBM

            3.2.1 GridSearch

            3.2.2 최적 파라미터 적용

            3.2.3 LightGBM_PCA삭제 한 모델 검증

     3.3 스태킹 앙상블(zip_code 7개로 축소한 모델)

            3.3.1 개별 Model 학습&평가

            3.3.2 스태킹 최종 모델 결과

     3.4 CV기반의 스태킹

            3.4.1 CV 함수정의

            3.4.2 개별 Model 학습&평가

     3.4 딥러닝

            3.4.1 TimeStep 설정

            3.4.2 Scaling

            3.4.3 Train, Test Dataset

            3.4.4. Model 설정

                      3.4.4.1 기본모델

                      3.4.4.2 LSTM

                      3.4.4.3 LSTM(layter1)

                      3.4.4.3 LSTM(layter2+ DropOut)

                      3.4.4.4 LSTM(layter3+ DropOut)

                      3.4.4.5 Gru

                      3.4.4.6 Gru(layter1)

                      3.4.4.7 Gru (layter2+ DropOut)

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

국민청원글을 대상으로 토픽모델링을 진행했습니다.

첫번째로 LDA를 활용하여 토픽모델링 함으로써 공통적인 속성을 알고자 했습니다

두번째로 ATM은 카테고리별 문서주제 분포를 파악하여 유사도를 분석해 보고자 했습니다 

corab url

1.  import/ 파일 불러오기

2. 텍스트 전처리 전 시각화

>각 청원글을 공백 기준으로 단어로 만듭니다

>공백으로 잘라준 total_tokens의 빈도를 확인합니다

> 전처리 전 청원글의 빈도 그래프 입니다. 보다시피 유의미한 의미를 찾기 힘들기 때문에

명사만을 추출하여 진행했습니다.

3. 텍스트 전처리 후 시각화

전처리 과정에 앞서서 먼저 author2doc를 만들어 줍니다 (전처리 후에 해도 상관없습니다) 

> 카테고리 별로 group by를 해줍니다. 총 882개의 청원글의 인덱스가 카테고리별로 이쁘게 모입니다.

cls2doc은 ATM의 author2doc에 해당됩니다.

이제 카테고리 중 "행정" 청원글을 가지고 텍스트 전처리/시각화를 진행하겠습니다.  

>cls2doc['행정']을 실행하면 문서 인덱스가 list로 나오기 때문에 for문을 이용할수 있습니다.

>행정 문서내용이 담긴 cls1을 doc으로 넘겨준 뒤 mecab.nouns(doc) 통해 명사만 추출하는데 

길이가 1이상인 것만 가져옵니다. 한글자 단어는 이 청원글에서 큰 의미를 주지 못하기 때문에 제외 시켜줬습니다.

각 명사별 빈도 top3가 "국민,정부,코로나" 인것을 확인할수 있습니다.

시각화를 하기 앞서서 불용어 처리를 위해 불용어 사전을 업로드 후 그래프 및 워드 클라우드를 시각화 합니다.

> 텍스트 전처리 후 상위 100개 단어를 가지고 그래피와 워드클라우 시각화 결과입니다.

ngram을 하지 않아서 각 단어가 어떤 연결성이 존재하는지는 알수 없지만 일자리에 관련된 청원글의 

키워드를 대략 파악할수 있었습니다.

4. topic modeling

lda/atm 모델에는 dictionary ,corpus가 필요합니다.  dictionary는 영어 단어장 처럼 청원글에 제시된 단어를 말합니다.

저는 mecab.nouns를 사용해서 dictionary를 만들었습니다. corpus는 각 단어를 (단어의 인덱스,빈도) 형식으로 수치화 시켜준 것을 말합니다. 이해를 돋기 위해서 아래 그림을 참고해주세요

>왼쪽이미자가 corpus, 오른쪽이미지는 인덱스를 단어로 변경했을때 모습입니다.

그러면 dictionary와 corpus를 어떻게 만드는지 아래 코드를 보면서 설명드리겠습니다.

>왼쪽그림에서 문서를 명사 추출 후 text_cleaning통해 한번더 걸러줍니다. 

명사를 추출하면서 특수문자나 url와 같은 문자들은 추출이 되지 않지만 한 글자만 남아 있는 것들 제거 하기 위해서 실행해 줬습니다.

왼쪽그림에서 한 글자을 제거하는 코드를 여러 방면으로 시도해 봤지만 모델을 돌릴때 인덱스 오류가 지속적으로 발생했습니다. 제 생각에는 한 글자들이 빠지면서 생긴 공백 때문이지 않을까 생각됩니다. 

좀더 쉽고 간단하게 할수 있는 방법을 찾다가 실습때 배웠던 text_cleaning 함수에 len(doc)>1를 추가 했는데 에러가 잡히지 않았습니다.

>명사화된 token들을 Dictionary,doc2bow를 통해 dictionary,corpus를 만들어줍니다.

4-1 LDA Model

>토픽 수를 결정하기 위해서 최소 4개~20개까지 각각 계산해 합니다.

>corpus, dictionary를 넣어주고 num_topic를 이용해서 4개일때,5개일때, ... , 20개일때의 lda를 각각 계산합니다.

주석처리된 부분은 모델을 저장 할때 사용합니다. 

>계산 결과 topic 8개 일때 coherence가 0.45로 가장 높은결과를 얻었습니다(passes=100일때)

coherence의 기준을 보면 0.45가 낮은 수준임을 알수가 있습니다. 즉 유사한 단어끼리 뭉쳤다고는 좋게 볼수가 없습니다. passes를 1000회 이상 높이면 지금보다 더 좋은 결과를 얻을수 있을 거라 생각됩니다. 

이어서 topic_num=8 일때의 lda결과를 보겠습니다.

> num_topic를 8로 설정하고, passes= 100 으로 모델을 실행합니다.

> topic label를 정해줍니다.

>각 토픽별 어떤 단어끼리 뭉쳤는지 출력합니다. 

>LDA model 결과로는 총 8개의 topic으로 구분되었습니다. 가장 많은 비중을 차지한 topic1은 코로나와 관련된 키워드 "코로나,방역,정부" 등이 있으며 topic3은 "주택,부동산,지역" 등 부동산과 관련된 topic임을 짐작 할수 있었습니다.

coherence가 낮기때문에 양질의 분류는 아닐지라도 어느정도 유추 할수 있는 topic들이 있었습니다. 

여기까지 전체 청원글에 대한 LDA를 이용하여 topic modeling을 진행해 보았습니다.

4-2 ATM_model

> author2doc이 추가 되었고 나머지는 LDA와 같습니다.

>실행결과 num_topic=8 일때 수치가 가장 높았지만 lda와 같은 0.45로 낮은 수치를 기록했습니다.

그러면 num_topic=8로 다시 설정하여 ATM을 실행하겠습니다.

>각 토픽별 word를 출력한 결과입니다. topic0은 범죄와 관련되어 있고, topic1은 사회적 거리두기 방역지침과 밀접해 보입니다. LDA와 마찬가지로 같은 num_topic(8), coherence를 보였습니다. 전반적으로 토픽의 성격도 비슷하게 묶였음을 확인할수 있었습니다. 그러면 author 간의 거리를 측정해서 유사도를 보겠습니다.

> 특정 저자가 어떠한 저자와 유사한지를 구하는 함수입니다.

즉, author-topic space에서 가장 가까운 벡터를 구해야 합니다. 거리를 재는 방법 중 hellinger distance를

사용합니다. hellinger distance는 두 개의 확률분포가 가우시안 분포를 따를 때 그 사이의 거리를 재는 방법입니다.

get_author_topics에는 해당 저자(카테고리)가 7개 토픽 중 차지하는 비율 정보가 있습니다.

이 정보로 가지고 저자 간의 거리를 구합니다.

(해당 함수는 강의 실습 코드를 참고했습니다.)

>일자리 author와 가장 유사도가 높은 것은 경제 민주화입니다.

따라서 이 두개의 author는 청원글의 주제가 겹칠수 있다는걸 알수 있습니다.

토픽모델링을 공부하면서 느낌점

부트캠프에서 크롤링->R의 워드클라우드로 프로젝트를 한적이 있었습니다. 그 당시에도 들었던 생각은 키워드만 가지고 뭔가 인사이트를 얻는게 충분치 않다 라는 느낌을 받았습니다. 이번계기로 NLP가 파고들수록 어려운 분야 라는걸 깨달았습니다. 문맥,뉘앙스 처럼 사람만이 이해하는 것들을 학습시키는게 상당히 어려울 거란걸 초보인 저도 가늠할수 있었습니다.  이처럼 어려운 분야의 전문가가 되는 제 모습을 상상하면서 더욱 매진해 볼 생각입니다 :) 

만일 틀린거라든가 보완해야 할 점이 있다면 언제점 말씀해주세요 ~!  

시작하기에 앞서

첫 블로그 내용으로는 BTS 이미지를 분류해봤습니다. BTS를 선택한 이유는 팬이기도 하고 유명하니까 자료 얻기도 쉽지 않을까 하는 생각에 선택했습니다! 그래서 BTS 멤버 이미지를 가지고 분류하는 코드를 작성해봤는데, 틀렸다거나 다른 이견이 있으시다면 언제든 말씀해 주시면 감사하겠습니다 ^^

자료는 다음 이미지를 크롤링했습니다. 한명 당 1300건을 수집했고, 훈련에 부적합한 이미지를 거른 후 총 7649장을 사용했습니다. 

업로드할 코드 이미지는 추가적인 설명 또는 막혔었던 부분을 설명하겠습니다. 

참고서적은 케라스 창시자에게 배우는 딥러닝, 핸즈온 머신러닝입니다. 

CNN 딥러닝 분류 모델

성능을 높히기 위해서 데이터 증식과 VGG16 사전 학습 모델을 적용시켜봤습니다.

큰 틀의 내용은 아래 목차를 올리겠습니다.

colab.research.google.com/drive/1n5m6i-4EuGOOWRnpgaMW1AziiijyaiDf?usp=sharing

필요한 모듈 import

데이터 전처리

채널이 RGB가 아닌 사진들을 제거합니다. 삭제 전 후 이미지 수가 같은 이유는 사전에 미리 코드를 돌렸습니다. 

대략 30장이 삭제되었는데 거의 15분 정도가 걸렸습니다. 

각 멤버별로 인덱스를 enumerate로 넣어줍니다. class 2_idx는 tfrecord를 통해서 모델의 label로 사용됩니다.

첫 번째 고민했던 부분

실습 때와 다르게 각 멤버별로 이미지 수가 달라서 train, test를 어떻게 나눌까 였습니다. 

"사실 처음 시도했을 때는 각 멤버별 이미지 수가 다르다는 걸 생각지 못했었습니다 ㅠ "

위 코드로 각 멤버별 이미지 수를 파악한 후에 7:3 비율로 분리했습니다.  

7:3으로 분리하는 사용자 정의 함수를 만든 후, 각 멤버 별로 적용했습니다. "이제 보니까 for문으로 좀 더 깔끔하게 할 수 있지 않았을까 생각이 듭니다."

train_dir, validation_dir으로 분리된 이미지의 수를 파악해봅니다. train 이미지는 5358, validation 이미지는 2291장입니다. ""각 멤버 별 train 수가 1000장이 안되기 때문에 기대에 맞는 성능은 나오지 못할 거란 걸 알고 있었지만 크롤링한 이미지를 모두 전수 검사를 하다 보니까 (3시간 소요) 힘이 들어서 양을 늘릴 엄두를 못 냈습니다""

TFRecord File 생성

tfrecord에 대해서 더 공부해야 할 부분이 많지만 장점을 말씀드리자면 대규모 데이터를 효율적으로 처리할 수 있다는 장점을 가지고 있습니다. tfrecord형식은 이진 레코드의 시퀀스를 저장할 수 있는 간단한 형식입니다.

경험상 dataset으로 만드는 과정이 비교적 간단하다고 생각합니다.

tfrecord file을 저장할 경로를 설정해줍니다.

tensorflow 홈페이지에 있는 코드입니다. 이 함수들을 통해서 일반 data를 tfrecord data로 변환합니다.

image는 bite로 직렬화 해서 만들고, 레이블은 int형식으로 지정해줍니다. validation도 동일하게 위 코드처럼 진행해서 이미지는 생략했습니다.

생성이 되었는지 확인해 봅니다.

hyper parameter를 설정합니다. 

tfrecord file를 data로 해석해주는 함수를 정의합니다. 해석된 data로 dataset를 만들 예정입니다. 그리고

IMG_SIZE= 224로 설정했습니다.

loss, acc 그래프를 그려본 결과입니다. loss감소, val_loss 증가 형태는 과적 합의 그래프입니다. epoch 수가 적기 때문에 단지 튀는 과정일 수도 있습니다. 따라서 acc 그래프도 같이 보았습니다.(오른쪽 사진)

acc 그래프에는 val_loss가 일정한 걸 보아서는 epoch를 늘려도 변함이 없어 보임으로 과적합임을 확인할 수 있었습니다. 과적합을 없애기 위해서 BatchNormalization과 데이터 증식을 적용시켜 봤습니다.

train 이미지에 다양한 옵션을 설정했습니다. validation 이미지에는 rescale만 적용시켜 줍니다. 

위 사진처럼 같은 이미지라도 다르게 만들어서 Data의 부족을 해결합니다.

class_mode에서 이진 분류를 한다면 binary를 넣고, 이진이 아닌 다중 분류라면 정수 라벨을 반환하는 sparse로 설정합니다. BTS는 7개의 라벨이 있기 때문에 sparse로 적용시켜줍니다.

파라미터 정리

2D one-hot 부호화된 라벨이 반환됩니다.

binary : 1D 이진 라벨이 반환됩니다

sparse : 1D 정수 라벨이 반환됩니다.

val_accuray: 29.73%

loss, acc 그래프를 본 결과 과적합이 사라 졌음을 볼 수 있었습니다. epoch를 더 늘리면 성능이 좀 더 올라갈 것으로 보입니다. 

이 결과를 통해 새로운 이미지를 test 해보겠습니다.

지민 사진을 넣었는데 예측 값으로 진이 나왔습니다.

과적합을 해결했지만 성능 측면에서 아쉬운 부분이 있습니다. 따라서, 사전학습을 통해 성능을 올려 보도록 하겠습니다.

사용될 사전학습 모델은 VGG16입니다.

모델을 imagenet을 통해 다운로드합니다.

이 외에 나머지 코드는 위에서 했던 코드와 상당수 일치하므로 새로 추가된 코드만 올리겠습니다

정확도가 전 모델에 비해 크게 상승했습니다. 그래프를 통해 과적합 여부를 쉽게 파악합니다.

epoch 10 이상부터는 과적합 전조 현상이 보입니다. v16 model에는 적용하지 않았지만 과적합을 줄이기 위해서 dropout를 시도해 볼 수도 있겠습니다.

마지막으로 아까와 같은 이미지로 테스를 진행해보겠습니다.

그 결과 

"지민" -> "김남준"으로 오답을 예측했습니다.

과적합이 발생했기 때문에 정답을 맞혔더라도 모델 수정은 반드시 이뤄져야 할 것 같습니다.

느낀 점

처음으로 CNN 프로젝트를 짧게나마 진행해봤습니다. 깔끔하게 정리된 이미지를 가지고 실습하다가 실전 이미지를 다루다 보니 많은 시행착오가 있었습니다. 전처리하는데 대부분 시간이 소요되었습니다. 아직 많이 부족하지만 배울게 많다는 게 인공지능 분야의 장점이라고 생각합니다. 감사합니다