DACON - 상점 신용카드 매출 예측(모델 구축과 검증)

Data

DACON - 상점 신용카드 매출 예측(모델 구축과 검증)

HC-Kang 2021. 5. 30. 15:54

드디어 데이콘 우승작 4장 신용카드 매출 예측의 끝이 보이기 시작한다.

오늘은 모델 구축과 검증부터 성능향상방법까지 끝을 볼 수 있기를.

아 이전에 먼저 다뤘어야겠지만, 이 우승자분께서는 특이하게도 R studio를 사용하셨다. 그런데 책에서는 이걸 rpy2 라이브러리를 활용해서 파이썬에 적용하는 방향으로 작성이 되어있었다.

처음에는 문제를 풀 때, 왜 굳이 R코드를 가상환경을 만들고 굳이 호환되는 버전까지 찾아가며 파이썬으로 옮겨서 작업을 하셨을까 싶은 의문이 있었는데, 알고보니 독자들이 뭐 하나라도 더 배울 수 있게 해 주겠다는 저자분들의 노고가 들어간 것이라고 생각 하려고 한다.

근데 그래도,,, 이렇게는 안할것같다. 차라리 R을 따로 돌리고 말지....

from rpy2.robjects.packages import importr # rpy2 내의 패키지를 불러올 importr클래스

utils = importr('utils') # utils 패키지를 임포트
utils.install_packages('forecast') # r의 forecast 패키지 설치
utils.install_packages('forecastHybrid') # r의 forecastHybrid 패키지 설치

forecast 패키지는 시계열 모델링을 제공하고, forecastHybrid는 이름과 같이 앙상블 예측을 위한 패키지이다.

위 코드를 실행하면 Secure CRAN mirror창에서 R 패키지를 받기 위한 창이 열린다.

그리고 여기부터 또 신기한 게 나온다.

import rpy2.robjects as robjects     # r함수를 파이썬에서 사용 가능하게 변환하는 모듈 
from rpy2.robjects import pandas2ri  # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈

pandas2ri.activate()  # pandas2ri를 활성화

다만 아무래도 어거지로 두개를 이어붙이는 것인 만큼, 내가보기엔 다소 좀 지저분하다고 느껴졌다...

아무래도 하나의 모듈로 뿅! 하는 것도 아니고, robjects 라는 함수 변환 모듈과, 판다스와 r을 이어주기 위한 pandas2ri라는 모듈을 별도로 설치해야 하고, pandas2ri는 또 별도로 activate까지 해줘야 한다는 점 때문이겠지.

auto_arima = """
	function(ts){
		library(forecast)  # forecast 패키지 로드
		d_params = ndiffs(ts)  # ndiffs() 시계열 자료의 차분 횟수 계산
		model = auto.arima(ts, max.p=2, d=d_params)  # p부분차수, d차분차수         
        	auto.arima 모델 생성 forecasted_data = forecast(model, h=3)  # 이후 3개월(h=3)을 예측
		out_df = data.frame(forecasted_data$mean)  # 예측값(3개월의 평균)을 R의 데이터프레임으로 변환 
        	colnames(out_df) = c('amount')  # amount라는 열로 이름을 지정
		out_df
	} 
"""

이부분이 신기했다. rpy2를 가지고 파이썬에서 r 코드를 돌릴 수 있다고만 알고있었는데, 그 방식이 색달랐다.

먼저 auto_arima라는 이름으로 str형태로 r코드를 통째로 짜넣는다. 모델은 이름과 같이 ARIMA 모델을 활용한다.

r 코드의 경우에는 우승자분께서 주석을 워낙 잘 달아주셔서 생략.

# r() 함수로 r 자료형을 파이썬에서 사용 가능하게 변경
auto_arima = robjects.r(auto_arima)
ts = robjects.r('ts')# r의 time series 자료형으로 만들어주는 함수 
c = robjects.r('c') # r의 벡터로 만들어주는 함수

auto_arima = robjects.r(auto_arima)의 형태로 robjects 함수를 이용해서 활용가능한 r코드 함수로 변환해준다.

이어서 'ts'와 'c'도 동일하게 r함수를 활용 가능하도록 정의해준다.

store_0 = resampling_data[resampling_data['store_id']==0] # 0번 상점 
start_year = int(min(store_0['year_month'])[:4]) # 영업 시작 년도 
start_month = int(min(store_0['year_month'])[5:]) # 영업 시작 월

맨 위 부터, resampling_data중에서 store_id가 0인 상점에 대해 뽑아내어 store_0에 저장하고,

min() 함수를 통해 각각 시작 년도와 월을 뽑아내어 저장한다.

train = ts(store_0['amount'], start=c(start_year, start_month), frequency=12)

이후 train data로 만들어 주기 위해 train 이라는 변수에 (start_year, start_month) 부터 time series 형태로 저장하고, 그 주기는 1년이므로 12로 지정한다.

#ensemble model
forecast = auto_arima(train)
np.sum(pandas2ri.ri2py(forecast).values) # 3개월 매출을 합산
	# 결과 : 2007226.9722140946

이후 train 자료를 가지고 이전에 정의한 auto_arima 모델 학습하고, numpy를 통해 forecast의 결과 합산.

대략적으로 이런식으로 모델을 운용할거다~ 라고 이해하면 될 것 같다.

이제 본격적으로 시계열 모델을 활용 해보고, 선택과 검증 과정을 시작한다.

우승자분은 자기회귀누적이동평균(ARIMA), 지수평활법, STL분해를 적용한 지수평활법 등 세가지 모델을 활용하였고, 여기에 로그 정규화를 조건부로 활용한 것을 토대로 앙상블해서 최종 결과를 제출했다.

1. ARIMA 모델

일단 당장 ARIMA부터 좀 알아보기 위해 내가 공부하면서 가장 도움 받았던 채널을 소개하면

https://www.youtube.com/watch?v=ma_L2YRWMHI

김성범 인공지능연구소 소장님 Youtube 채널

위 자료에서 정말 알아듣기 쉽게 설명 해 주신다.

간단하게 정리해보면 ARIMA 모델은 AR모델과 MA모델을 섞은 것이고, AR모델은 AutoRegression, MA 모델은 Moving Average모델로 각각 자기회귀, 이동평균모델이다. 처음엔 뭐 이름이 이래 복잡하냐 했는데. 지금 보니까 이만큼 정직한 이름도 없을듯 싶다.

AR 모델을 수식으로 나타내면 아래와 같다

전체적인 형태는 이름답게 자기 자신의 함수값을 활용한 회귀모델이다. 간단하게 설명하면, y의 값을 이전의 y값에 계수를 곱한값들의 합으로 표현한다는 것이다.

MA 모델의 수식은

위와 같이 y의 값을 현재와 과거 오차들의 가중합으로 표현하겠다는 내용이다.

이걸 쓰면서도 실시간으로 궁금한점, 이해가 안되는 점이 솟구치는데, 그걸 다 찾아보다가는 결국 포스팅을 못 할 것 같아서, 나중에 다시 다뤄보기로 하고 넘어가자.

위 두가지를 합친 ARMA 모델도 있는데, 일단 여기서는 넘어가고, 차분을 포함한 ARIMA 모델로 바로 들어가자.

지금 당장 드는 생각으로는 아마 차분 없이 ARMA모델만 사용하는 경우가 상당히 적어서 아예 다루지 않은것인가 싶기도 하다.

ARIMA 모델의 수식은

단순히 AR 모델과 MA 모델의 합이라고 볼 수 있다.

이런 모델 선정 과정에서 원래 ACF, PACF 등의 활용도 필요하겠지만, 여기서는 우선 다 넘어가도록 하자.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈

# pandas2ri를 활성화 
pandas2ri.activate()

auto_arima = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        d_params = ndiffs(ts) # 시계열 자료의 차분 횟수 계산
        model = auto.arima(ts, max.p=2, d=d_params) # auto.arima 모델 생성
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""
# r() 함수로 r 자료형을 파이썬에서 사용 가능
auto_arima = robjects.r(auto_arima)# str 형식으로 정의된 auto_arima
ts = robjects.r('ts')# r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수

여기까지 코드는 위에서 살펴본 모델 소개와 동일하다.

final_pred = []

for i in tqdm(resampling_data.store_id.unique()):
    store = resampling_data[resampling_data['store_id']==i]
    start_year = int(min(store['year_month'])[:4]) ## 영업 시작 년도
    start_month = int(min(store['year_month'])[5:]) ## 영업 시작 월
    # R의 ts 함수로 time series 데이터로 변환
    train = ts(store['amount'], start=c(start_year, start_month), frequency=12) 
    # 자동회귀누적이동평균 model
    forecast = auto_arima(train)
    # 3개월 매출을 합산, final_pred에 추가
    final_pred.append(np.sum(pandas2ri.ri2py(forecast).values))

이전 모델 설정에서는 단순히 모델의 형태 설명을 위해 store_0만 예로 실행하였는데, 우리의 자료에서는 약 2천개의 상점이 예측을 위해 대기중이다. 따라서 위의 코드를 통해 해당 상점들을 대상으로 매출을 예측한다.

먼저 final_pred라는 빈 리스트를 만들어서 우리가 예측한 3개월 매출 합을 저장할 수 있도록 한다.

이어서 for문을 활용해서 모든 상점에 반복시행하며 이 때 이전과 같이 unique() 함수와 tqdm라이브러리를 활용한다.

그 외 store 지정, start 지정과 훈련, 결과예측 과정은 이전과 동일하다.

submission = pd.read_csv('./submission.csv')
submission['amount'] = final_pred
submission.to_csv('submission.csv', index=False)

이후 final_pred에 저장된 결과를 submission.csv 파일에 작성하여 저장한다.

위에서 작성한 ARIMA 모델의 제출 결과는 844,384점으로, 18위를 기록했다고 한다.

아 참고로, 맨 첫 포스트에서 다루려다가 간단히 언급만 하고 넘어갔는데, 대회의 채점 기준은 MAE(평균절대오차)이다. 따라서 숫자가 작을수록 높은 점수이다.

2. 지수평활법

우선 단순지수평활법의 식은 아래와 같다.

다시금 개념적으로만 간단히 설명하면, 현재의 값에 가장 큰 가중치를 두고, 과거의 자료일수록 점점 더 비중을 줄여나간다는 것으로, 이런 정의로 인해서 단순지수평활법은 추세나 계절성이 있는 데이터를 예측하는데에는 적합하지 못하다.

물론 이런 점을 보완하기 위해 홀트의 선형추세 기법도 있다.

이 식에서는 단순지수평활법의 단점을 보완하기 예측식과 세 개의 평활식으로 구성되어 있고, 세 개의 평활식은 각각 수준식(ℓ)과 추세식(b), 계절식(s)으로 구성되어있다.

이후에도 또 홀트-윈터스의 덧셈, 곱셈 기법으로 나누어 진다는데... 여기도 일단 그만 알아보자. 정말 공부에는 끝이 없구나.....

여하튼 이 코드에서 사용하는 지수평활법은, 정확히 뭐라고는 나와있지 않은 것 같다. 다만 언급되어있는 것은 AIC가 최소화되는 모델을 추정한다고 하는데, 이 부분도 정확히 이해하지 못했다.

일단 AIC가 뭔지를 좀 알아보면, 풀 네임은 Akaike Information Criterion, 일단 Information Criterion을 먼저 알아봐야 할 것 같다.

이부분을 이해하는데에 우주먼지 님의 포스트에서 많은 도움을 받았다.

내가 이해한 개념으로, 우리가 모델을 구성 할 때 보통 그 성능의 지표로 정확도 혹은 가능성과, 모델에 필요한 변수의 수를 기준으로 판단하는데 이를 지수화 한 것으로, 보통은 아래와 같은 trade-off 관계를 보인다.

방향이 정확히 맞는지는 모르겠지만, 변수의 수가 늘어날수록 지나치게 늘어난 Variance로 인해 과적합되는 것이 떠오르는데, 이부분을 바탕으로 이해하면 아래의 수식이 좀 더 이해가 되는것 같다.

따라서 위 식에서 AIC는 모형의 적합도가 높을수록 그 수치가 낮아지지만, 변수가 늘어나면 그만큼 점수가 늘어나는 패널티를 주는 것으로, 개념적으로 AIC가 최소화되는 지점이 대체로 가장 적합한 모델일 것이라는 점을 이해 할 수 있다.

이와 비슷한 개념으로 BIC(Bayesian Information Criterion)도 있는데, 그 식은 아래와 같다.

위 식과 비교해서 본다면 쉽게 이해 할 수 있다. AIC와는 다르게 변수의 수에 log를 씌워 그 강도를 약화시켰다. 따라서 변수의 수가 조금 많더라도, 더욱 정확한 모델을 찾는 방향으로 집중하는구나~ 하고 대략적으로 이해 할 수 있다.

다시 원문으로 돌아와서, 지수평활법 코드에 대해서 알아보자.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈

# pandas2ri를 활성화 
pandas2ri.activate()

ets = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        model = ets(ts) # AIC가 낮은 지수평활 모델을 찾음 
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""
# r() 함수로 r 자료형을 파이썬에서 사용 가능
ets = robjects.r(ets)# str 형식으로 정의된 ets
ts = robjects.r('ts')# r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수

이전의 코드와 대부분 동일하지만, 함수의 종류가 auto arima에서 ets 로 변경되었다.

final_pred = []

for i in tqdm(resampling_data.store_id.unique()):
    store = resampling_data[resampling_data['store_id']==i]
    start_year = int(min(store['year_month'])[:4]) # 영업 시작 년도
    start_month = int(min(store['year_month'])[5:]) # 영업 시작 월
    # R의 ts 함수로 time series 데이터로 변환
    train = ts(store['amount'], start=c(start_year, start_month), frequency=12) 
    # 지수평활법l
    forecast = ets(train)
    # 3개월 매출을 합산, final_pred에 추가
    final_pred.append(np.sum(pandas2ri.ri2py(forecast).values))
    
    submission = pd.read_csv('./submission.csv')
    submission['amount'] = final_pred
    submission.to_csv('submission.csv', index=False)

역시나 이전 코드와 동일하게 final_pred 리스트에 결과를 담아준 뒤, submission 파일에 그 값을 저장한다.

지수평활법 모델의 경우 제출 결과는 794,263점을 기록했고, 당시 10위였다고 한다. 약 5만점의 향상을 기록했다.

3. STL 분해를 적용한 지수평활법

이번엔 당연히 STL분해가 무엇인지 부터 알아봐야겠다.

STL(Seasonal and Trend decomposition using Loess)은 시게열 분해 기법으로, 데이터를 EDA에서 다루었던 계절성과 추세, 그리고 잔차로 분해해서 분석하는 기법이다.

코드를 먼저 살펴보자.

from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.pyplot as plt

store_0 = time_series(resampling_data, 0)
# STL 분해
stl = seasonal_decompose(store_0.values, freq=12)
stl.plot()
plt.show()

STL 분해는 r이 아닌 파이썬에서 제공하는 패키지인 statsmodels의 seasonal_decompose 함수를 활용한다. 또한 그래프를 시각화하기 위해 matplotlib도 같이 불러와 준다.

이후 예시로 0번 상점의 자료를 시각화 해본다.

그래프에서 볼 수 있다시피, 맨 위에는 데이터의 관측값 총계가 보여지고, 그 아래로 순서대로 추세와 계절성, 잔차를 보여준다.

STL분해를 적용한 지수평활법은 일반적인 지수평활법에서 이런 STL분해를 통해 얻은 추세와 계절성이라는 정보를 추가적으로 적용하여 예측을 보정하는 방식으로 R의 stlm()함수를 이용한다.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈

# pandas2ri를 활성화 
pandas2ri.activate()
stlm = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        model = stlm(ts, s.window="periodic", method='ets') # STL 분해 후 지수평활법을 통한 예측 
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""
ets = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        model = ets(ts) # AIC가 낮은 지수평활 모델을 찾음 
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""
# r() 함수로 r을 파이썬에서 사용 가능
stlm = robjects.r(stlm)# str 형식으로 정의된 stlm
ets = robjects.r(ets)# str 형식으로 정의된 ets
ts = robjects.r('ts')# r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수

코드를 보니 뭐 이전과 크게 다른점은 없어보인다. 다만, 이전에 사용한 ets가 남아있다. 이게 왜 남아있는거지? 하는 궁금증이 있었는데, 다음 코드와 해설에서 그 부분이 해결됐다.

final_pred = []

for i in tqdm(resampling_data.store_id.unique()):
    store = resampling_data[resampling_data['store_id']==i]
    data_len = len(store)
    
    start_year = int(min(store['year_month'])[:4]) # 영업 시작 년도
    start_month = int(min(store['year_month'])[5:]) # 영업 시작 월
    # R의 ts 함수로 time series 데이터로 변환
    train = ts(store['amount'], start=c(start_year, start_month), frequency=12) 
    # STL 분해를 적용한 지수평활 model
    if data_len > 24:
        forecast = stlm(train)
    # 지수평활 model
    else:
        forecast = ets(train)    # 3개월 매출을 합산, final_pred에 추가
    final_pred.append(np.sum(pandas2ri.ri2py(forecast).values))
    
    submission = pd.read_csv('./submission.csv')
    submission['amount'] = final_pred
    submission.to_csv('submission.csv', index=False)

이전의 코드들과 다르게 조건문이 하나 추가되어있다. 데이터의 길이를 data_len에 입력하고 그 길이가 24이상인 경우에 한해서 stlm() 함수를 적용한다. 이는 stlm 함수를 사용하기 위해서는 2개 시즌 이상의 자료가 필요하기 때문이고, 우리는 주기를 12로 설정했으니 최소 24개 이상의 자료가 있어야만 함수를 사용 할 수 있기에 별도 조건문을 두어서 영업일이 24개월 미만인 상점에 대해서는 일반적인 지수평활법인 ets() 함수를 적용해준다.

그리고 동일하게 submission 파일에 결과를 저장하여 제출한 결과는 869,178점으로 21등을 기록했다고 한다. 여기서 우리는 분석기법을 추가로 적용한다고 해서 무조건적으로 좋은 결과를 기대하는 건 맞지 않는다는 점을 알 수 있다.

이후 추가적인 성능 향상을 위해 우승자분은 두 가지 방법을 이용했다. 첫째는 데이터 전처리에 로그 정규화를 실시하는 것이고, 둘째는 R의 forecastHybrid 패키지를 이용해 앞의 세 가지 예측모델을 앙상블하는 방법이다. 먼저 로그 정규화 코드부터 살펴보자.

시계열 데이터에서 로그 정규화를 활용하는 이유는, 로그의 개념을 이해하면 좀 더 쉽게 이해 할 수 있다.

일단 로그함수는 위 그래프에서 보다시피 1보다 큰 양수에 대해서는 그 변화량을 큰 폭으로 줄여준다. 다시 말해, 어지간히 큰 변화가 아니라면 거의 변화가 없는 것으로 만들어준다는 뜻이다. 일종의 댐퍼가 된다고 보면 되겠다. 이런 로그 정규화의 효과를 알아보기 위해서 우승자분은 0번 상점에 대해서 train set, test set을 분할하여 성능을 검증했다.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈
import numpy as np

# pandas2ri를 활성화 
pandas2ri.activate()

auto_arima = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        d_params = ndiffs(ts) # 시계열 자료의 차분 횟수 계산
        model = auto.arima(ts, max.p=2, d=d_params) # auto.arima 모델 생성
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""

여기까지는 이전 코드와 거의 동일하다.

# r() 함수로 r 자료형을 파이썬에서 사용 가능
auto_arima = robjects.r(auto_arima)
ts = robjects.r('ts')# r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수
log = robjects.r('log')# 로그 변환 함수
exp = robjects.r('exp')# 로그 역변환 함수

로그 정규화에 필요한 R의 log 함수와 exp 함수도 같이 불러와준다.

# 0번 상점 추출
store_0 = resampling_data[resampling_data['store_id']==0]
start_year = int(min(store_0['year_month'])[:4]) # 영업 시작 년도
start_month = int(min(store_0['year_month'])[5:]) # 영업 시작 월

# train, test 분리
train = store_0[store_0.index <= len(store_0)-4]
test = store_0[store_0.index > len(store_0)-4]

# R의 ts 함수로 r의 time series 자료형으로 변환
train_log = ts(log(train['amount']), start=c(start_year, start_month), frequency=12) # log 정규화 
train = ts(train['amount'], start=c(start_year, start_month), frequency=12) # log 정규화를 하지 않음

# model arima
forecast_log = auto_arima(train_log)
forecast = auto_arima(train)

# pred
pred_log = np.sum(pandas2ri.ri2py(exp(forecast_log)).values) #로그 역변환 후 3개월 합산
pred = np.sum(pandas2ri.ri2py(forecast).values) #3개월 매출을 합산

# test(2018-12~2019-02)
test = np.sum(test['amount'])

# mae
print('log-regularization mae: ', abs(test-pred_log))
print('mae:', abs(test-pred))

0번 상점의 자료를 불러와서 마지막 4달을 테스트셋으로 분리한 후 정규화 자료와 일반 자료를 각각 train_log와 train 변수에 저장한다.

이후 예측값을 각각 pred_log 와 pred에 저장한 후 비교해준다. 다만 여기서 주의할 점으로, 로그 정규화를 해준 예측값의 경우에는 반드시 exp() 함수를 통해 이전값으로 돌려보내 주어야 한다는 것이다.

이후 마지막 4개월치 자료를 test 변수에 저장해주고, 그 값을 각각 pred_log, pred와 비교해준다.

결과는 로그정규화 MAE : 2401.97, 일반 MAE : 5884.67로 0번 상점에서는 로그 정규화된 자료를 활용했을 때 더욱 정확한 예측값을 보인다는 것을 알 수있다.

다만 주의할 점으로는, 앞에서 언급했듯 뭔가 새로운 방법을 적용한다고 해서 무조건 성능이 좋아지는 것은 아니라는 것이다. 따라서 로그 정규화도 무조건적인 성능향상을 기대 할 수는 없다. 그런데 그 조건에 대해서 명확하게 알 수 있는 기준이 없기때문에 우승자분은 여러번 반복적으로 로그 정규화를 실시하였고, 변동계수가 0.3 미만인 경우에만 로그 정규화를 진행하기로 결정했다.

여기서 변동계수라는 개념이 또 등장했는데, 그 식은 아래와 같다.

쉽게 말해서, 평균에 비해서 자료의 흩어진 정도가 0.3 미만인 경우에 로그정규화가 효과가 있다는 뜻으로, 내가 이해한대로 풀자면 그 이상으로 변동폭이 큰 매장의 경우에는 정규화 이전 값으로도 충분하게 변동을 감지 할 수 있으니 정규화 하지 않는 것이 낫고, 매출이 0.3보다 작은 진동폭으로 쓸데없이 흔들리는 상점의 경우에는 로그 정규화로 그 변화폭을 잡아주면 더욱 정확한 예측값을 얻을 수 있다는 뜻인 것 같다. 음.. 아니면 0.3 이상의 변동이 있는 상점만이 의미있는 변화이고, 그 이하는 그저 잡음으로 이해하는 게 맞을 수도 있겠다.

여하튼, 여기서 명확히 알아야 할 점은 변동계수 0.3과 같은 기준은 어디까지나 경험적으로 얻어내야 하고, 이것이 분석가의 기량이라는 점인 것 같다.

따라서 이후에는 변동계수 0.3 미만인 상점의 데이터에 로그 정규화를 실시하여 예측하는 코드를 살펴보자.

# 매출 변동 계수를 구하는 함수
def coefficient_variation(df, i):
    cv_data = df.groupby(['store_id']).amount.std()/df.groupby(['store_id']).amount.mean()
    cv = cv_data[i]
    return cv

위에서 살펴 본 식 대로, 변동계수 CV를 구하는 함수를 정의한다.

또한 적용할 모델은 가정 우수한 성능을 보였던 지수평활법으로 한다.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈
import numpy as np

# pandas2ri를 활성화 
pandas2ri.activate()

ets = """
    function(ts){
        library(forecast) # forecast 패키지 로드
        model = ets(ts) # AIC가 낮은 지수평활 모델을 찾음 
        forecasted_data = forecast(model, h=3) # 이후 3개월(h=3)을 예측
        out_df = data.frame(forecasted_data$mean) # 예측값을 R의 데이터프레임으로 변환
        colnames(out_df) = c('amount') # amount라는 열로 이름을 지정
        out_df
    }
"""

# r() 함수로 r 자료형을 파이썬에서 사용 가능
ets = robjects.r(ets)
ts = robjects.r('ts') # r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수
log = robjects.r('log') # 로그 변환 함수
exp = robjects.r('exp')# 로그 역변환 함수

동일한 코드이니 넘어가자.

final_pred = []

for i in tqdm(resampling_data.store_id.unique()):
    store = resampling_data[resampling_data['store_id']==i]
    start_year = int(min(store['year_month'])[:4]) # 영업 시작 년도
    start_month = int(min(store['year_month'])[5:]) # 영업 시작 월
    
    cv = coefficient_variation(resampling_data, i)
    # 매출액 변동 계수가 0.3 미만인 경우만 log를 씌움
    if cv < 0.3:
        train_log = ts(log(store['amount']), start=c(start_year,start_month), frequency=12) 
        # ets model
        forecast_log = ets(train_log)
        final_pred.append(np.sum(pandas2ri.ri2py(exp(forecast_log)).values))
    # 매출액 변동 계수가 0.3 이상인 경우
    else:
        train = ts(store['amount'], start=c(start_year,start_month), frequency=12)
        # 지수평활법
        forecast = ets(train)
        final_pred.append(np.sum(pandas2ri.ri2py(forecast).values))
        
submission = pd.read_csv('./submission.csv')
submission['amount'] = final_pred
submission.to_csv('submission.csv', index=False)

일반 지수평활법 모델에서 cv값을 확인하는 조건문이 추가되었다.

for문과 unique() 함수를 통해 각각의 상점에 대해서 동일하게 시작지점을 정한 뒤, cv값을 저장하여 if문에서 판단 지표로 활용한다. 이후 cv값이 0.3 미만이면 로그 정규화를 하여 예측값을 도출하고, 그 외의 경우에는 일반적인 지수평활법을 적용한다. 이후 예측값을 다시 submission 파일에 저장한다.

결과는 793,546점으로 일반 모델에 비해서 약 700점가량 매우 소폭 증가했다. 그러나 최종 결과에서 1, 2등 점수차가 겨우 67점에 불과했다고 하니, 마냥 무시할만한 점수차이는 아닌 것 같다.

마지막으로 forecastHybrid 패키지를 통해 앙상블을 진행한다.

네이버 국어사전에 나오는 앙상블의 정의는 아래와 같다.

이처럼, 예측모델의 앙상블이란 여러 모델의 예측결과를 합치는 것인데, 우승자분은 여러 종류의 앙상블 중 가장 기본적인 평균 앙상블을 사용했다고 한다. 즉, 자기회귀누적이동평균(ARIMA)와 지수평활법, STL분해 적용 지수평활법 세 모델의 예측값을 평균내어 최종 예측치를 구하는 것이다.

import rpy2.robjects as robjects # r 함수를 파이썬에서 사용 가능하게 변환하는 모듈
from rpy2.robjects import pandas2ri # 파이썬 자료형과 R 자료형의 호환을 도와주는 모듈
import numpy as np

# pandas2ri를 활성화 
pandas2ri.activate()

hybridModel = """
    function(ts){
        library(forecast)
        library(forecastHybrid)
        d_params=ndiffs(ts)
        hb_mdl<-hybridModel(ts, models="aes", # auto_arima, ets, stlm
                        a.arg=list(max.p=2, d=d_params), # auto_arima parameter
                        weight="equal") # 가중치를 동일하게 줌(평균)
        forecasted_data<-forecast(hb_mdl, h=3) # 이후 3개월(h=3)을 예측
        outdf<-data.frame(forecasted_data$mean)
        colnames(outdf)<-c('amount')
        outdf
    }
""" 

# r() 함수로 r 자료형을 파이썬에서 사용 가능
hybridModel = robjects.r(hybridModel)
ts = robjects.r('ts') # r 자료형 time series 자료형으로 만들어주는 함수
c = robjects.r('c') # r 자료형 벡터를 만들어주는 함수
log = robjects.r('log') # 로그 변환 함수
exp = robjects.r('exp')# 로그 역변환 함수

뭐 대단한 변화는 없는 것 같다. r에 들어가는 코드에서 hybridModel() 함수에 models = 'aes' 라니. 이게 원래 이렇게 간단하게 들어가는건가 싶다. 이런식이라면 이 코드는 평균만 가능하고 가중치는 못넣는건가..? 하는 등의 궁금증이 많이 튀어나온다. 확인해보니 RSA, DFA 등의 모델도 있는가보다. 다만 정확한 자료는 아직 못 찾았다. 이부분은 R을 좀 더 공부 해보아야겠다.

final_pred = []

for i in tqdm(resampling_data.store_id.unique()):
    store = resampling_data[resampling_data['store_id']==i]
    start_year = int(min(store['year_month'])[:4]) # 영업 시작 년도
    start_month = int(min(store['year_month'])[5:]) # 영업 시작 월
    
    cv = coefficient_variation(resampling_data, i)
    # 매출액 변동 계수가 0.3 미만인 경우만 log를 씌움
    if cv < 0.3:
        train_log = ts(log(store['amount']), start=c(start_year,start_month), frequency=12) 
        # 앙상블 예측
        forecast_log = hybridModel(train_log)
        final_pred.append(np.sum(pandas2ri.ri2py(exp(forecast_log)).values)) 
    # 매출액 변동 계수가 0.3 이상인 경우
    else:
        train = ts(store['amount'], start=c(start_year,start_month), frequency=12)
        # 앙상블 예측
        forecast = hybridModel(train)
        final_pred.append(np.sum(pandas2ri.ri2py(forecast).values))

submission = pd.read_csv('./submission.csv')
submission['amount'] = final_pred
submission.to_csv('submission.csv', index=False)

또한 로그 정규화도 동일하게 실시 해주고 파일로 저장한다.

결과는 741,880점으로, 기존 최고점에 비해서 약 5만점이 넘게 상승하는 결과를 보였다. 이 모델을 통해 우승자분은 최종 1위를 기록 할 수 있었다.

하지만 대회 종류 후에도 해커톤 형식으로 계속해서 새로운 모델들이 등장했고 현재는 퍼블릭 점수 기준 585,952으로 박예인 님이 1위를 기록했다.

이번 코드리뷰를 하면서 정말 모르는점이 많아서 단순 리뷰를 하는것임에도 꽤나 많은 시간이 소요되었고 많은 점을 배울 수 있었다. 다만 몇 가지 이상했던 점, 아직 이해가 안되는 점 등이 있었는데, 그런 점을 보충하기 위해 추후에 3위를 기록하신 지성민 님의 코드도 리뷰를 해보아야 겠다.

1, 2, 3, 5장도 각각 리뷰를 해야겠는데, 아직 정말로 갈길이 멀다. 꾸준히 갈 수 있는 사람이 되자.