728x90
*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다
저번 포스팅에서 시계열 모델링을 수동으로 자동화 하는법을 포스팅했습니다.
이번 포스팅에서도 이어서 진도를 나갔는데요
원래 자동화 해주는 라이브러리가 있지만 사용하기에는 제약이 있습니다. 그래서 어떤 제약이 있는지에 대한 실습과
prophet(시계열 자동화 라이브러리) 사용방법과 알고리즘을 실습을 통해 알아 보겠습니다.
Auto-ARIMA Library
먼저 제약이 있는 자동화 라이브러리는 소규모 데이터에는 작동이 가능합니다.
import pmdarima as pm
from pmdarima import model_selection
#Data loading and Split
data = pm.datasets.load_wineind() #와인의 생산량 데이터
train,, test = model_selectiom.train_test_split(data,train_size=150) #데이터 분리
#모델링
autoarima = pm.auro_arima(train, stationarity=False, with_intercept=True,
start_p=0, d=None, start_q=0,
max_p=5, max_d=2, max_1=5, #p,d,q의 범위를 지정
seasonal_P=0, D=None, max_Q=5,
max_P=5, max_D=2, max_Q=0,
max_order = 100, maxiter =10,
information_criterion = 'bic'#어떤값으로 성능을 평가할지에 대한 기준선택
#경고문 패스: suppress_warnings
trace = True, suppress_warnings=True) #trace: 파라미터가 업데이트되는 걸 보고자할때
display(autoarima.summary())
#predict_in_sample: train data를 예측
pred_tr_Ts_autoarima = autoarima.predict_in_sample()
#test 예측: predict 사용
pred_tr_ts_autoarima = autoarima.predict(n_periods= len(train))#점추정만 출력
pred_te_ts_autoarima = autoarima.predict(n_periods= len(test),
return_conf_inf = True)[0] #0번째: 점추정
pred_te_ts_autoarima_ci = autoarima.predict(n_periods = len(test),
return_conf_int=True)[1] #1번째: 구간추정
# 검증
Score_ts_autoarima, Resid_tr_ts_autoarima, Resid_te_ts_autoarima = evaluation_trte(pd.DataFrame(train),pred_tr_ts_autoarima,
pd.DataFrame(test), pred_te_ts_autoarima,
graph_on=True)
display(Score_ts_autoarima)
ax = pd.DataFrame(test).plot(figsize=12,4))
pd.DataFrame(pred_te__ts_autoarima, columns =['prediction']).plot(kind='line', linewidth=3, fontsize=20,
ax=ax)
ax.fill_between(pd.DataFrame(pred_te_ts_autoarima_ci).index,
pd.DataFrame(pred_te_ts_autoarima_ci).iloc[:,0],
pd.DataFrame(pred_te_ts_autoarima_ci).iloc[:,1], color='k', alpha=0.15)
plt.show()
# 잔차진단
error_analysis(Resid_tr_ts_autoarima, ['Error'], pd.DataFrame(train), graph_on = True)
가장 BIC가 낮은 파라미터는 (2,1,1)(1,0,0)[12]이며 아래는 model summary 결과입니다.
자기상관이 존재 하지 않으며 정규분포가 아니라고 하지만 Skew, Kurtosis는 정규분포과 비슷하게 나왔다
분산은 0.05이므로 유의수준 5%라고 할때 등분산 이라고 말하고 있다.
아래는 성능 평가와 예측과 실측치를 시각화한 그래프이며 ACF,PACF를 통해 정말 자기상관이 0인지를 살펴본다
ACF, PACF 모두 자기 상관이 존재 한다. PACF같은 현상은 저번 포스팅에서 말씀드린 것처럼
뒤쪽으로 갈수록 unstable 현상은 계절성 패턴이 있음을 의심해 봐야 합니다.
auto_arima를 자전거 수요 데이터로 실행하면 메모리 부족 현상으로 인해 에러가 걸립니다. 위 데이터는
데이터가 적어서 실행이 가능했던 것입니다. 아래 실습은 데이터만 다르고 모두 동일합니다.
import pmdarima as pm
# auto_arima를 자전거 수요에 적용
# Data Loading
# location = 'https://raw.githubusercontent.com/cheonbi/DataScience/master/Data/Bike_Sharing_Demand_Full.csv'
location = './Data/BikeSharingDemand/Bike_Sharing_Demand_Full.csv'
raw_all = pd.read_csv(location)
# Feature Engineering
raw_fe = feature_engineering(raw_all)
###Reality###
target = ['count_trend', 'count_seasonal', 'count_Day', 'count_Week', 'count_diff']
raw_feR = feature_engineering_year_duplicated(raw_fe, target)
#data split
Y_colname = ['count']
X_remove = ['datetime', 'DateTime', 'temp_group', 'casual', 'registered']
X_colname = [x for x in raw_fe.columns if x not in Y_colname+X_remove]
X_train_feR, X_test_feR, Y_train_feR, Y_test_feR = datasplit_ts(raw_feR, Y_colname, X_colname, '2012-07-01')
### Reality ###
target = ['count_lag1', 'count_lag2']
X_test_feR = feature_engineering_lag_modified(Y_test_feR, X_test_feR, target)
###Scaling###
X_train_feRS, X_test_feRS = feature_engineering_scaling(preprocessing.Normalizer(), x_train_feR, X_test_feR)
###############
### multicollinearity ###
print('Number_of_Total_X: ', len(X_train_feRS.columns))
X_colname_vif = feature_engineering_XbyVIF(X_train_feRS,12)
print('Number_of_Selected_X: ', len(X_colname_vif))
X_train_feRSM, X_test_feRSM = X_train_feRS[X_colname_vif].copy(), X_test_feRS[X_colname_vif].copy()
# 모델링
## Parameter Setting
trend_diff_order =0
seasonal_diff_order, seasonal_order =0,24
#SARIMAX
#auto_arima가 안되는 이유: 정보를 저장하는데에 메모리 부족
fit_ts_autoarima = pm..auto_arima(Y_train_feR, satationary= False,
stationary=False,
with_intercept=True,
start_p=0, d=None, start_q=0,
max_p=2, max_d=1, max_q=2,
seasonal=True, m=24,
start_P=0, D=None, start_Q=0,
max_P=2, max_D=1, max_Q=2,
max_order=30, maxiter=3,
stepwise=False,
exogenous=X_train_feRSM,
information_criterion='aic',
trace=True, suppress_warnings=True)
display(fit_ts_autoarima.summary())
pred_tr_ts_autoarima = fit_ts_autoarima.predict_in_sample(exogenous=X_train_feRSM)
pred_tr_ts_autoarima = fit_ts_autoarima.predict(n_periods=len(Y_train_feR), exogenous=X_train_feRSM)
pred_te_ts_autoarima = fit_ts_autoarima.predict(n_periods=len(Y_test_feR),
exogenous=X_test_feRSM,
return_conf_int=True)[0]
pred_te_ts_autoarima_ci = fit_ts_autoarima.predict(n_periods=len(Y_test_feR),
exogenous=X_test_feRSM,
return_conf_int=True)[1]
# 검증
Score_ts_autoarima, Resid_tr_ts_autoarima, Resid_te_ts_autoarima = evaluation_trte(Y_train_feR, pred_tr_ts_autoarima,
Y_test_feR, pred_te_ts_autoarima, graph_on=True)
display(Score_ts_autoarima)
ax = pd.DataFrame(Y_test_feR).plot(figsize=(12,4))
pd.DataFrame(pred_te_ts_autoarima,
index=Y_test_feR.index, columns=['prediction']).plot(kind='line',
xlim=(Y_test_feR.index.min(),Y_test_feR.index.max()),
linewidth=3, fontsize=20, ax=ax)
ax.fill_between(pd.DataFrame(pred_te_ts_autoarima_ci, index=Y_test_feR.index).index,
pd.DataFrame(pred_te_ts_autoarima_ci, index=Y_test_feR.index).iloc[:,0],
pd.DataFrame(pred_te_ts_autoarima_ci, index=Y_test_feR.index).iloc[:,1], color='k', alpha=0.15)
plt.show()
# 잔차진단
error_analysis(Resid_tr_ts_autoarima, ['Error'], Y_train_feR, graph_on=True) 
Prophet Library
이 라이브러리의 특징은 현실성을 상당히 반영할수 있다는 점입니다. 도메인 지식을 접목 시킬수 있다는 것이 장점입니다.
# Data Loading
# location = 'https://raw.githubusercontent.com/cheonbi/DataScience/master/Data/Bike_Sharing_Demand_Full.csv'
location = './Data/BikeSharingDemand/Bike_Sharing_Demand_Full.csv'
raw_all = pd.read_csv(location)
# Feature Engineering
raw_fe = feature_engineering(raw_all)
### Reality ###
target = ['count_trend', 'count_seasonal', 'count_Day', 'count_Week', 'count_diff']
raw_feR = feature_engineering_year_duplicated(raw_fe, target)
###############
# Data Split
# Confirm of input and output
Y_colname = ['count']
X_remove = ['datetime', 'DateTime', 'temp_group', 'casual', 'registered']
X_colname = [x for x in raw_fe.columns if x not in Y_colname+X_remove]
X_train_feR, X_test_feR, Y_train_feR, Y_test_feR = datasplit_ts(raw_feR, Y_colname, X_colname, '2012-07-01')
### Reality ###
target = ['count_lag1', 'count_lag2']
X_test_feR = feature_engineering_lag_modified(Y_test_feR, X_test_feR, target)
###############
### Scaling ###
X_train_feRS, X_test_feRS = feature_engineering_scaling(preprocessing.Normalizer(), X_train_feR, X_test_feR)
###############
### Multicollinearity ###
print('Number_of_Total_X: ', len(X_train_feRS.columns))
X_colname_vif = feature_engineering_XbyVIF(X_train_feRS, 12)
print('Number_of_Selected_X: ', len(X_colname_vif))
X_train_feRSM, X_test_feRSM = X_train_feRS[X_colname_vif].copy(), X_test_feRS[X_colname_vif].copy()위 자전거 수요 예측 과정을 선행 후 모델을 진행합니다.
Default
import timeit
start = timeit.default_timer()
from fbprophet import Prophet
from fbprophet.plot import plot_plotly, plot_components_plotly
from fbprophet.plot import add_changepoints_to_plot
# Rearrange of Data
Y_train_feR_prophet = Y_train_feR.reset_index() #시간축을 인덱스로 가지고 있으면 안된다.
#Y값과 시간정보가 들어간 컬럼이 있어야 한다.
#컬럼명도 아래처럼 입력해야 작동된다.
Y_train_feR_prophet.columns =['ds','y'] #Prophet실행 시 무조건 해야 하는것
#하나의 데이터프레임으로 만들어서 fit 함수에 적용
YX_train_prophet = pd.concat([Y_train_feR_prophet, X_train_feRSM.reset_index().iloc[:,1:]],
axis=1)
# Applying Prophet Model
#change point: 그래프에 변화를 보일때 시점(CP)
#선형적인 효과로 추세가 있다.
fit_reg1_feRSM_prophet = Prophet(growth='linear',#추세가 선형인지, 2차함수인지 로그함수인지 임의로 지정
'''
여기서 말하는 trend를 이해하기 쉽게 설명.
한 회사의 사장님이 구속되는 일이 발생하면 주식(trend)은 상향에서 하향으로 변화 할것이다.
변화하는 시점을 CP로 본다
'''
'''
trend라는 건 잘개 쪼개진 것들의 합으로 구성된 것
잘개 쪼깨는 기준이 CP가 된다
이 알고리즘이 원하는 것: CP 정보만 주면 예측하겠다는 것이다.
이때, CP의 정보가 아래이다.
'''
#미래든 과거든 이벤트가 발생한 날을 적용
changepoints=None, # CP가 발생하는 시점들의 list ['2012-01-01']
#예상되는 cp를 넣으면 해당 범위에서 분석한다(max 값)
n_changepoints=25, # CP의 수
#CP값을 비율로 나타낸것
# 1000개 중 0.8, 즉 최대로 나올수 있는 CP는 800개 라는 뜻
changepoint_range=0.8, # CP의 기존 데이터 수 대비 최대 비율
# 퍼센트 개념x, 많을수록 CP를 많이 추정하려고 한다.
changepoint_prior_scale=0.05, # CP 추정 민감도로 높을수록 민감
'''
Seasonality
'''
#계절성이 y값에 더하기로 표현되는지 곱하기로 표현되는지 지정해준다.
seasonality_mode='additive', # 계절성 모델: 'additive' or 'multiplicative'
#계절성 가중치를 결정
seasonality_prior_scale=10.0, # 계절성 추정 민감도로 높을수록 민감
#데이터를 보고 자동으로 일~연계절성을 자동으로 추정해준다
#없으면 추정 안한다.
#auto: 있으면 추정 없으면 추정 안함
yearly_seasonality='auto', # 연계절성
weekly_seasonality='auto', # 주계절성
daily_seasonality='auto', #일계절성
# Holiday
# 2021-12-25 크리스마스 <-이런형식으로 데이터프레임으로 입력
holidays=None, # 휴일 또는 이벤트 시점 dataframe
holidays_prior_scale=10.0, # 휴일 추정 민감도로 높을수록 민감
'''
Others
'''
#미래에 나올 예측구간을 설정한다. -> 기존의 범위(train)에서 80% 반영한 범위를 사용하겠다
interval_width=0.8, # 추세 예측 정확도 구간범위
mcmc_samples=0) # 계절성 예측 정확도 제어(몬테카를로) changepoints =None
fit_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.fit(YX_train_prophet)
#예측할 데이터의 freq가 뭔지, 몇개 시점을 예측해야 하는지를 입력
#실제 예측한다는 용도보다는 freq가 적당한지를 검증하기 위한 정도->도메인 전문가에게 검증을 받는 용도
forecast = fit_reg1_feRSM_prophet.make_future_dataframe(freq='H',
periods=Y_test_feR.shape[0])
pred_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.predict(forecast) #train,test 예측값을 함께 출력한다.
#붙어있는 train,test를 분리한다
pred_tr_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[:Y_train_feR.shape[0]-1,['yhat']])
pred_te_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[Y_train_feR.shape[0]:, ['yhat']])
# Result
display(pred_reg1_feRSM_prophet[['ds'.'yhat','yhat_lower','yhat_upper']])
#trend 시각화
fig =fit_reg1_feRSM_prophet.plot(pred_reg1_feRSM_prophet) #fit한 것과 예측한 것을 같이 시각화fit한 것과 예측한 것을 같이 시각화
#gca: fig안에 라벨값, 폰트 값 등 셋팅 값을 가지고있음
add_chnagepoints_to_plot(fig.gca(), fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet)
plt.xticks(rotation=90)
plt.show()
#seasonal 시각화
# plot_components: 일,연,달 계절성, 예측값과 관련된 값이 나온다
fig = fig_reg1_feRSM_prophet.plot_components(pred_reg1_feRSM_prophet)
ax = fig.get_axes()
#아래처럼 하나씩 작업해야 시각화가 됨
for i in range(len(ax)):
ax[i].tick_params(axis='x', labelsize=15, rotation=90)
plt.show()
## using plotly
# display(plot_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
#plot_components_plotly: 핏팅 값과 예측값을 각각 시각화
# display(plot_components_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# Evaluation
Score_reg1_feRSM_prophet, \
Resid_tr_reg1_feRSM_prophet, \
Resid_te_reg1_feRSM_prophet = evaluation_trte(Y_train_feR,
pred_tr_reg1_feRSM_prophet,
Y_test_feR,
pred_te_reg1_feRSM_prophet,
graph_on=True)
display(Score_reg1_feRSM_prophet)
# Error Analysis
error_analysis(Resid_tr_reg1_feRSM_prophet, ['Error'], X_train_feRSM, graph_on=True)
stop = timeit.default_timer()
print('Time: ', stop - start) 
검정색이 실제 train값, 파란색이 train,test 예측값
주황색: 점추정값
빨간색 점선: CP
각 해당되는 빈도로 계절성을 뽑은 패턴입니다. 이를 통해 각 빈도의 흐름을 파악할수 있습니다.
본 강의에서는 trend와 seasonal 파라미터 조정 시 어떻게 변화하는지에 대한 실습이 많이 있지만 그중에서 최종적으로
적용된 최종 실습 하나만 올리도록 하겠습니다
Final
import timeit
start = timeit.default_timer()
from fbprophet import Prophet
from fbprophet.plot import plot_plotly, plot_components_plotly
from fbprophet.plot import add_changepoints_to_plot
# Rearrange of Data
Y_train_feR_prophet = Y_train_feR.reset_index()
Y_train_feR_prophet.columns = ['ds', 'y']
YX_train_prophet = pd.concat([Y_train_feR_prophet,
X_train_feRSM.reset_index().iloc[:,1:]],
axis=1)
# Applying Prophet Model_비선형적 추세를 넣기 위한 작업
### saturating ###
saturating = 'logistic' #logistic를 쓸때는 Y값의 상한과 하한을 컬럼으로 생성해야한다.
if saturating == 'logistic':
#예상되는 train 상한과 하한을 정해준다
YX_train_prophet['cap'] = 900
YX_train_prophet['floor'] = 0
### holiday ### -> 수동으로 넣을때 방법
newyear = pd.DataFrame({'holiday': 'newyear',
'ds': pd.to_datetime(['2011-01-01', '2012-01-01']),
'lower_window': 0, # 이전날짜 포함범위
'upper_window': 1}) # 이후날짜 포함범위
christmas = pd.DataFrame({'holiday': 'christmas',
'ds': pd.to_datetime(['2011-12-25', '2012-12-25']), # 미래 시점도 포함
'lower_window': 0, # 이전날짜 포함범위
'upper_window': 1}) # 이후날짜 포함범위
holidays = pd.concat([newyear, christmas])
fit_reg1_feRSM_prophet = Prophet(growth=saturating,
# Trend
changepoints=None, # CP가 발생하는 시점들의 list ['2012-01-01']
n_changepoints=25, # CP의 수
changepoint_range=0.9, # CP의 기존 데이터 수 대비 최대 비율
changepoint_prior_scale=0.5, # CP 추정 민감도로 높을수록 민감
# Seasonality
seasonality_mode='additive', # 계절성 모델: 'additive' or 'multiplicative'
seasonality_prior_scale=20.0, # 계절성 추정 민감도로 높을수록 민감
yearly_seasonality='auto', # 연계절성
weekly_seasonality='auto', # 주계절성
daily_seasonality='auto', #일계절성
# Holiday
holidays=holidays, # 휴일 또는 이벤트 시점 dataframe
holidays_prior_scale=5.0, # 휴일 추정 민감도로 높을수록 민감 (이전결과참고반영)
# Others
interval_width=0.8, # 추세 예측 정확도 구간범위
mcmc_samples=0) # 계절성 예측 정확도 제어
### holiday ###
#US를 적으면 미국의 휴일,공유일이 자동으로 인식된다.(한국은 없음)
fit_reg1_feRSM_prophet.add_country_holidays(country_name='US')
### monthly seasonality(계절성 추가하기)###
#성능이 안좋아져서 주석처리 후 위에 holiday 수동으로 직접 넣음
# fit_reg1_feRSM_prophet.add_seasonality(name='monthly', period=24*7*4, fourier_order=1)
### extra feature ###
# X(설명변수)들을 고려하여 예측하기
target_name = X_train_feRSM.columns[:5]
for col in target_name:
fit_reg1_feRSM_prophet.add_regressor(col)
fit_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.fit(YX_train_prophet)
forecast = fit_reg1_feRSM_prophet.make_future_dataframe(freq='H',
periods=Y_test_feR.shape[0])
### saturating ###
# test 상하한 범위 지정
if saturating == 'logistic':
forecast['cap'] = 700
forecast['floor'] = 0
### extra feature ###
forecast = pd.concat([forecast,
pd.concat([X_train_feRSM[target_name],
X_test_feRSM[target_name]], axis=0).reset_index().iloc[:,1:]],
axis=1)
pred_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.predict(forecast)
pred_tr_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[:Y_train_feR.shape[0]-1, ['yhat']])
pred_te_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[Y_train_feR.shape[0]:, ['yhat']])
# Result
display(pred_reg1_feRSM_prophet[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
fig = fit_reg1_feRSM_prophet.plot(pred_reg1_feRSM_prophet)
add_changepoints_to_plot(fig.gca(), fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet)
plt.xticks(rotation=90)
plt.show()
fig = fit_reg1_feRSM_prophet.plot_components(pred_reg1_feRSM_prophet)
ax = fig.get_axes()
for i in range(len(ax)):
ax[i].tick_params(axis="x", labelsize=15, rotation=90)
plt.show()
## using plotly
# display(plot_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# display(plot_components_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# Evaluation
Score_reg1_feRSM_prophet, \
Resid_tr_reg1_feRSM_prophet, \
Resid_te_reg1_feRSM_prophet = evaluation_trte(Y_train_feR,
pred_tr_reg1_feRSM_prophet,
Y_test_feR,
pred_te_reg1_feRSM_prophet,
graph_on=True)
display(Score_reg1_feRSM_prophet)
# Error Analysis
error_analysis(Resid_tr_reg1_feRSM_prophet, ['Error'], X_train_feRSM, graph_on=True)
stop = timeit.default_timer()
print('Time: ', stop - start) 
위 실습은 어떻게 적용하는지에 대한 교육 초점을 뒀기 때문에 성능 개선을 위한 정밀한 작업은 없었음을 알려드립니다.
728x90
'Data Diary' 카테고리의 다른 글
| 2021-04-07(시계열데이터 심화14_공적분 모형) (0) | 2021.04.07 |
|---|---|
| 2021-04-06(시계열데이터 심화13_VAR&Granger Causality ) (0) | 2021.04.06 |
| 2021-04-03 기록 (0) | 2021.04.03 |
| 2021-04-02 기록(1차 방문) (0) | 2021.04.03 |
| 2021-04-01(시계열데이터 심화11_선형확률과정 분석싸이클1) (0) | 2021.04.01 |