지난번 포스팅에서는 백색잡음과정(White noise process), 확률보행과정(Random walk process), 정상확률과정(Stationary process)에 대해서 소개하였습니다. (https://rfriend.tistory.com/691

 

지난번 포스팅에서 특히 ARIMA와 같은 시계열 통계 분석 모형이 정상확률과정(Stationary Process)을 가정한다고 했습니다. 따라서 시계열 통계 모형을 적합하기 전에 분석 대상이 되는 시계열이 정상성 가정(Stationarity assumption)을 만족하는지 확인을 해야 합니다. 

 

[ 정상확률과정 (stationary process) 정의 ]

1. 평균이 일정하다. 
2. 분산이 존재하며, 상수이다. 
3. 두 시점 사이의 자기공분산은 시차(時差, time lag)에만 의존한다. 

 

정상 시계열 (stationary time series) 여부를 확인하는 방법에는  3가지가 있습니다. 

  [1] 시계열 그래프 (time series plot)

  [2] 통계적 가설 검정 (statistical hypothesis test)

  [3] 자기상관함수(ACF), 편자기상관함수(PACF)

 

정상성 확인하는 방법 (how to check the stationarity)

 

이번 포스팅에서는  이중에서 통계적 가설 검정 (Statistical Hypothesis Test) 을 이용해 정상성(stationarity) 여부를 확인하는 방법을 소개하겠습니다. 

 

- (a) Augmented Dickey-Fuller (“ADF”) test

- (b) Kwiatkowski-Phillips-Schmidt-Shin (“KPSS”) test

 

 

ADF 검정은 시계열에 단위근(unit root)이 존재하는지의 여부를 검정함으로써 정상 시계열인지 여부를 판단합니다. 단위근이 존재하면 정상 시계열이 아닙니다. 

 

단위근(unit root)이란 확률론의 데이터 검정에서 쓰이는 개념입니다. 주로 ‘단위근 검정’의 형식으로 등장합니다. 일반적으로 시계열 데이터는 시간에 따라 일정한 규칙을 가짐을 가정합니다. 따라서 매우 복잡한 형태를 갖는 시계열 데이터라도 다음과 같은 식으로 어떻게든 단순화시킬 수 있을 것이라 생각해볼 수 있습니다.

즉, ‘t시점의 확률변수는 t-1, t-2 시점의 확률변수와 관계를 가지면서 거기에 에러가 포함된 것’이라는 의미입니다. 여기서 편의를 위해 y0=0이라 가정합니다.  이제 아래의 방정식을 볼까요.

여기서 m=1이 위 식의 근이 된다면 이때의 시계열 과정을 단위근을 가진다고 합니다. 단위근 모형은 주로 복잡한 시계열 데이터를 단순하게나마 계산하려 할 때 사용됩니다.

 

 

Python으로 가상의 시계열 데이터셋(1개의 정상 시계열, 3개의 비정상 시계열)을 만들어서 위의 ADF test, KPSS test 를 각각 해보겠습니다. 위의 정상확률과정의 정의에 따라서, 추세(trend)가 있거나, 분산(variance)이 시간의 흐름에 따라 변하거나(증가 또는 감소), 계절성(seasonality)을 가지고 있는 시계열은 정상성(stationarity) 가정을 만족하지 않습니다. 

 

- (1) 정상 시계열 데이터 (stationary time series)

- (2) 추세를 가진 비정상 시계열 데이터 (non-stationary time series with trend)

- (3) 분산이 변하는 비정상 시계열 데이터 (non-stationary time series with changing variance)

- (4) 계절성을 가지는 비정상 시계열 데이터 (non-stationary time series with seasonality)

 

 

이제 (1)~(4)번 데이터별로 (a) ADF test, (b) KPSS test 로 정상성 여부를 차례대로 가설 검정해보겠습니다. 

 

 

(1) 정상 시계열 데이터 (stationary time series)

 

먼저, 자기회귀과정(auto-regressive process)을 따르는 AR(1) 과정의 시계열 데이터를 임의로 만들어보겠습니다. 

 

import numpy as np
import matplotlib.pyplot as plt

## Stationary Process 
## exmaple: AR(1) process 
np.random.seed(1) # for reproducibility
z_0 = 0 
rho = 0.6 # correlation b/w z(t) and z(t-1) 
z_all = [z_0] 

for i in range(200): 
    z_t = rho*z_all[i] + np.random.normal(0, 0.1, 1) 
    z_all.append(z_t) 
    
## plotting 
plt.rcParams['figure.figsize'] = (10, 6)
plt.plot(z_all)
plt.title("stationary time series : AR(1) process", fontsize=16)

## adding horizonal line at mean position 0 
plt.axhline(0, 0, 200, color='red', linestyle='--', linewidth=2) 
plt.show()

stationary process time series, 정상 시계열 

 

 

 

(1-a) ADF 검정 (Augmented Dickey-Fuller test)

 

Python의 statsmodels 모듈에 있는 adfuller 메소드를 사용해서 ADF 검정을 위한 사용자 정의함수를 정의해보겠습니다. 

 

#! pip install statsmodels

## UDF for ADF test
from statsmodels.tsa.stattools import adfuller
import pandas as pd

def adf_test(timeseries):
    print("Results of Dickey-Fuller Test:")
    dftest = adfuller(timeseries, autolag="AIC")
    dfoutput = pd.Series(
        dftest[0:4],
        index=[
            "Test Statistic",
            "p-value",
            "#Lags Used",
            "Number of Observations Used",
        ],
    )
    for key, value in dftest[4].items():
        dfoutput["Critical Value (%s)" % key] = value
    print(dfoutput)

 

 

이제 위에서 만든 ADF 검정 사용자정의함수 adf_test() 를 사용해서 정상시계열 z_all 에 대해서 ADF 검정을 해보겠습니다. p-value 가 8.74e-11 이므로 유의수준 5% 하에서 귀무가설 (H0: 단위근(unit root)이 존재한다. 즉, 정상 시계열이 아니다)을 기각하고 대립가설 (H1: 단위근이 없다. 즉, 정상 시계열이다.) 을 채택합니다. (맞음 ^_^)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(z_all)  # p-value 8.740428e-11 => stationary

# Results of Dickey-Fuller Test:
# Test Statistic                -7.375580e+00
# p-value                        8.740428e-11
# #Lags Used                     1.000000e+00
# Number of Observations Used    1.990000e+02
# Critical Value (1%)           -3.463645e+00
# Critical Value (5%)           -2.876176e+00
# Critical Value (10%)          -2.574572e+00
# dtype: float64

 

 

이번에는 Python의 statsmodels 모듈을 사용해서 KPSS 검정 (Kwiatowski-Phillips-Schmidt-Shin Test) 을 위한 사용자 정의함수를 정의해보겠습니다. 

 

## UDF for KPSS test
from statsmodels.tsa.stattools import kpss
import pandas as pd

def kpss_test(timeseries):
    print("Results of KPSS Test:")
    kpsstest = kpss(timeseries, regression="c", nlags="auto")
    kpss_output = pd.Series(
        kpsstest[0:3], index=["Test Statistic", "p-value", "Lags Used"]
    )
    for key, value in kpsstest[3].items():
        kpss_output["Critical Value (%s)" % key] = value
    print(kpss_output)

 

 

이제 정상 시계열 z_all 에 대해서 위에서 정의한 kpss_test() 함수를 사용해서 정상성 여부를 확인해보겠습니다. p-value 가 0.065 로서 유의수준 10% 하에서 귀무가설 (H0: 정상 시계열이다) 를 채택하고, 대립가설 (H1: 정상 시계열이 아니다) 를 기각합니다. (유의수준 5% 하에서는 대립가설 채택). (맞음 ^_^)

* 이때, KPSS 검정은 ADF 검정과는 귀무가설과 대립가설이 정반대이므로 해석에 조심하시기 바랍니다.  

 

kpss_test(z_all) # p-value 0.065035 => stationary at significance level 10%

# Results of KPSS Test:
# Test Statistic           0.428118
# p-value                  0.065035
# Lags Used                5.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(2) 추세를 가진 비정상 시계열 데이터 (non-stationary time series with trend)

 

다음으로, 추세(trend)를 가지는 가상의 비정상(non-stationary) 시계열 데이터를 만들어보겠습니다. 평균이 일정하지 않으므로 비정상 시계열이 되겠습니다. 

 

## time series with trend
np.random.seed(1) # for reproducibility
ts_trend = 0.01*np.arange(200) + np.random.normal(0, 0.2, 200)

## plotting
plt.plot(ts_trend)
plt.title("time series with trend", fontsize=16)

plt.show()

time series with trend (non-stationary process)

 

 

(2-a) ADF 검정 (ADF test)

 

위의 추세를 가지는 비정상 시계열 데이터에 대해 ADF 검정 (ADF test)를 실시해보면, p-value가 0.96 으로서 유의수준 5% 하에서 귀무가설 (H0: 시계열이 단위근을 가진다. 즉, 정상 시계열이 아니다) 을 채택하고, 귀무가설 (H1: 시계열이 단위근을 가지지 않는다. 즉, 정상 시계열이다.) 을 기각합니다. (맞음 ^_^)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(ts_trend) # p-value 0.96 => non-stationary

# Results of Dickey-Fuller Test:
# Test Statistic                   0.125812
# p-value                          0.967780
# #Lags Used                      10.000000
# Number of Observations Used    189.000000
# Critical Value (1%)             -3.465431
# Critical Value (5%)             -2.876957
# Critical Value (10%)            -2.574988
# dtype: float64

 

 

 

(2-b) KPSS 검정 (KPSS test)

 

추세(trend)를 가지는 시계열에 대해 KPSS 검정(KPSS test)을 실시해보면, p-value 가 0.01 로서 귀무가설 (H0: 정상 시계열이다) 를 기각하고, 대립가설 (H1: 정상 시계열이 아니다) 를 채택합니다. (맞음 ^_^)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_trend) # p-valie 0.01 => non-stationary

# Results of KPSS Test:
# Test Statistic           2.082141
# p-value                  0.010000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(3) 분산이 변하는 비정상 시계열 데이터 (non-stationary time series with changing variance)

 

다음으로, 추세는 없지만 시간이 흐름에 따라서 분산이 점점 커지는 가상의 비정상 시계열을 만들어보겠습니다. 분산이 일정하지 않기 때문에 정상 시계열이 아닙니다. 

 

## time series with increasing variance
np.random.seed(1) # for reproducibility
ts_variance = []

for i in range(200):
    ts_new = np.random.normal(0, 0.001*i, 200).astype(np.float32)[0]
    ts_variance.append(ts_new)

## plotting    
plt.plot(ts_variance)
plt.title("time series with increasing variance", fontsize=16)

plt.show()

time series with increasing variance (non-stationary process)

 

 

(3-a) ADF 검정 (ADF test)

 

위의 시간이 흐름에 따라 분산이 커지는 비정상 시계열에 대해 ADF 검정(ADF test)을 실시하면, p-value가 5.07e-19 로서 유의수준 5% 하에서 귀무가설(H0: 단위근을 가진다. 즉, 정상 시계열이 아니다)를 기각하고, 대립가설(H1: 단위근을 가지지 않는다. 즉, 정상 시계열이다)를 채택합니다. (땡~! 틀림 -_-;;;)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root. (not stationary)
# -- Alternate Hypothesis: The series has no unit root. (stationary)

adf_test(ts_variance) # p-vaue 5.07e-19 => stationary (Opps, wrong result -_-;;;)

# Results of Dickey-Fuller Test:
# Test Statistic                -1.063582e+01
# p-value                        5.075820e-19
# #Lags Used                     0.000000e+00
# Number of Observations Used    1.990000e+02
# Critical Value (1%)           -3.463645e+00
# Critical Value (5%)           -2.876176e+00
# Critical Value (10%)          -2.574572e+00
# dtype: float64

 

 

(3-b) KPSS 검정 (KPSS test)

 

위의 시간이 흐름에 따라 분산이 커지는 비정상 시계열에 대해 KPSS 검정(KPSS test)을 실시하면, p-value가 0.035 로서 유의수준 5% 하에서 귀무가설(H0: 정상 시계열이다)를 기각하고, 대립가설(H1: 정상 시계열이 아니다)를 채택합니다. (딩동댕! 맞음 ^_^)  (ADF test 는 틀렸고, KPSS test 는 맞았어요.)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_variance) # p-value 0.035 => not stationary

# Results of KPSS Test:
# Test Statistic           0.52605
# p-value                  0.03580
# Lags Used                3.00000
# Critical Value (10%)     0.34700
# Critical Value (5%)      0.46300
# Critical Value (2.5%)    0.57400
# Critical Value (1%)      0.73900
# dtype: float64

 

 

 

(4) 계절성을 가지는 비정상 시계열 데이터 (non-stationary time series with seasonality)

 

마지막으로, 코사인 주기의 계절성(seasonality)을 가지는 가상의 비정상(non-stationary) 시계열 데이터를 만들어보겠습니다. 

 

## time series with seasonality
## generating x from 0 to 4pi and y using numpy
np.random.seed(1) # for reproducibility
x = np.arange(0, 50, 0.2) # start, stop, step 
ts_seasonal = np.cos(x) + np.random.normal(0, 0.2, 250)

## ploting 
plt.plot(ts_seasonal) 
plt.title("time series with seasonality", fontsize=16)

plt.show()

time series with seasonality (non-stationary process)

 

 

(4-a) ADF 검정 (ADF test)

 

위의 계절성(seasonality)을 가지는 비정상 시계열에 대해서, ADF 검정(ADF test)을 실시하면, p-value가 3.14e-16 로서 유의수준 5% 하에서 귀무가설(H0: 단위근을 가진다. 즉, 정상 시계열이 아니다)를 기각하고, 대립가설(H1: 단위근을 가지지 않는다. 즉, 정상 시계열이다)를 채택합니다. (땡~! 틀림 -_-;;;)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(ts_seasonal) # p-value 3.142783e-16 => stationary (Wrong result. -_-;;;)

# Results of Dickey-Fuller Test:
# Test Statistic                -9.516720e+00
# p-value                        3.142783e-16
# #Lags Used                     1.600000e+01
# Number of Observations Used    2.330000e+02
# Critical Value (1%)           -3.458731e+00
# Critical Value (5%)           -2.874026e+00
# Critical Value (10%)          -2.573424e+00
# dtype: float64

 

 

 

(4-b) KPSS 검정 (KPSS test)

 

위의 계절성(seasonality)을 가지는 비정상 시계열에 대해서, KPSS 검정(KPSS test)을 실시하면, p-value가 0.1 로서 유의수준 10% 하에서 귀무가설(H0: 정상 시계열이다)를 기각하고, 대립가설(H1: 정상 시계열이 아니다)를 채택합니다. (딩동댕! 맞음 ^_^)  (ADF test 는 틀렸고, KPSS test 는 맞았어요. 유의수준 5% 하에서는 둘 다 틀림.)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_seasonal) # p-value 0.1 => non-stationary (at 10% significance level)

# Results of KPSS Test:
# Test Statistic           0.016014
# p-value                  0.100000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

이상의 정상성 여부에 대한 통계적 가설 검정 결과를 보면,

(1) 정상 시계열에 대해 ADF test, KPSS test 모두 정상 시계열로 가설 검정 (모두 맞음 ^_^)

(2) 추세가 있는 시계열에 대해 ADF test, KPSS test 가 모두 정상 시계열이 아니라고 정확하게 가설 검정 (모두 맞음 ^_^)

(3) 분산이 변하는 시계열에 대해 ADF test 는 정상 시계열로 가설 검정 (틀렸음! -_-;;;), KPSS test 는 정상 시계열이 아니라고 가설 검정 (맞음 ^_^)

(4) 계절성이 있는 시계열에 대해 ADF test 는 정상 시계열로 가설 검정 (틀렸음! -_-;;;), KPSS test 는 정상 시계열이 아니라고 가설 검정 (맞음 ^_^)

 

ADF test 는 추세가 있는 비정상 시계열에 대해서는 정상 시계열이 아님을 잘 검정하지만, 분산이 변하거나 계절성이 있는 시계열에 대해서는 정상성 여부를 제대로 검정해내지 못하고 있습니다.

 

반면에 KPSS test 는 위의 4가지의 모든 경우에 정상성 여부를 잘 검정해내고 있습니다. 

 

통계적 가설 검정 외에 시계열 도표 (time series plot)을 꼭 그려보고 눈으로도 반드시 확인해보는 습관을 들이시기 바랍니다. 

 

 

[ Reference ]

- ADF test: https://en.wikipedia.org/wiki/Augmented_Dickey%E2%80%93Fuller_test 

- KPSS test: https://en.wikipedia.org/wiki/KPSS_test

- ADF test and KPSS test in Python: https://www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html

- 단위근 (unit root) : https://www.scienceall.com/%EB%8B%A8%EC%9C%84%EA%B7%BCunit-root-%E5%96%AE%E4%BD%8D%E6%A0%B9/

 

다음번 포스팅에서는 정상성 가정을 만족하지 않는 시계열 데이터에 대해 데이터 변환을 통해 정상화 시키는 방법을 소개하겠습니다. 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

이전 포스팅에서 시계열 데이터를 예측하는 분석 방법 중의 하나로서 지수평활법(exponential smoothing)을 소개하였습니다. 지수평활법은 통계적인 이론 배경이 있다기 보다는 경험적으로 해보니 예측이 잘 맞더라는 경험에 기반한 분석기법으로서, 데이터가 어떤 확률과정(stochastic process)을 따라야 한다는 가정사항이 없습니다. 

 

반면에, ARIMA (AutoRegressive Integreated Moving Average) 모형은 확률모형을 기반으로 한 시계열 분석 기법으로서, 정상확률과정(Stationary Process)을 가정합니다. 만약 시계열 데이터가 정상성(stationarity) 가정을 만족하지 않을 경우 데이터 전처리를 통해서 정상시계열(stationary time series)로 변환해 준 다음에야 ARIMA 모형을 적용할 수 있습니다. 

 

이번 포스팅에서는 

   (1) 백색잡음과정 (White Noise Process)

   (2) 확률보행과정 (Random Walk Process)

   (3) 정상확률과정 (Stationary Process)

에 대해서 소개하겠습니다.  특히, 백색잡음과정, 확률보행과정과 비교해서, ARIMA 모형이 가정하는 "정상확률과정"에 대해서 유심히 살펴보시기 바랍니다. 

 

 

[ 백색잡음과정(white noise) vs. 확률보행과정 (random walk) vs. 정상확률과정 (stationary process) ]

white noise vs. random walk vs. stationary process

 

(1) 백색잡음과정 (White Noise Process)

 

신호처리 분야에서 백색잡음(white noise)은 다른 주파수에서 동일한 강도를 가지고, 항상 일정한 스펙트럼 밀도를 발생시키는 무작위한 신호를 의미합니다. 백색잡음 개념은 물리학, 음향공학, 통신, 통계 예측 등 다양한 과학과 기술 분야에서 광범윟게 사용되고 있습니다. 

시계열 데이터 통계 분석에서 백색잡음과정(white noise process)은 평균이 0, 분산은 유한한 상수(sigma^2)인 확률분포(random variables)로 부터 서로 상관되지 않게(uncorrelated) 무작위로 샘플을 추출한 이산 신호(discrete signal)를 말합니다.  

 

white noise process

 

특히, 평균이 0 인 동일한 정규분포로부터 서로 독립적으로 샘플을 추출한 백색잡음을 가법 백색 가우시언 잡음(additive white Gaussian noise) 라고 합니다. 아래 Python 코드는 평균이 0, 분산이 1인 정규분포로부터 200개의 백색잡음 샘플을 추출하고 시각화 해본 것입니다. 

 

import numpy as np
import matplotlib.pyplot as plt

## generating white noise from normal distribution
x1 = np.random.normal(0, 1, 200) # additive white Gaussian noise

plt.plot(x1, linestyle='-') # marker='o'

## addidng horizontal line at mean '0'
plt.axhline(0, 0, 200, color='red', linestyle='--', linewidth=2)
plt.show()

white Gaussian noise

 

 

 

(2) 확률보행과정 (Random Walk Process)

 

수학에서 확률보행과정(random walk process)은 어떤 수학적인 공간 (예: 정수)에서 연속적인 무작위 스텝의 이동경로를 묘사하는 수학적인 객체를 말합니다. 

 

가장 기본적인 확률보행과정의 예를 들자면, 동전을 던져서 앞면이 나오면(0.5의 확률) 위(up)로 +1 계단을 올라가고, 동전의 뒷면이 나오면(0.5의 확률) 아래(down)로 -1 계단을 내려가는 보행과정을 생각해 볼 수 있습니다.  액체나 기체 속에서 분자가 이동하는 경로를 추적한 브라운 운동(Brownian motion), 수렵채집을 하는 동물의 탐색경로, 요동치는 주식의 가격, 도박꾼의 재무상태 등을 확률보행과정으로 설명할 수 있습니다. 

 

 

아래의 Python 코드는 위로(up)로 0.3의 확률, 아래(down)로 0.7의 발생확률을 가지고 계단을 오르락 내리락 하는 확률보행과정을 만들어보고, 이를 시각화해 본 것입니다. 마치 잠깐 상승하다가 하락장으로 돌아선 주식시장 차트 같이 생기지 않았나요? 

 

# Probability to move up or down
prob = [0.3, 0.7]  
  
# starting position
start = 0
rand_walks = [start]

# creating the random points
rr = np.random.random(200)
downp = rr < prob[0]
upp = rr >= prob[1]

# random walk process
# z(t) = z(t-1) + a(t), where a(t) is white noise
for idown, iup in zip(downp, upp):
    rand_walks.append(rand_walks[-1] - idown + iup)

# plot 
plt.plot(rand_walks)
plt.show()

random walk process (example)

 

 

 

(3) 정상확률과정 (Stationary Process)

 

서두에서 ARIMA 시계열 통계모형은 시계열 데이터의 정상성(stationarity) 조건을 만족하는 경우에 사용할 수 있다고 했습니다. 정상확률과정은 다음의 세 조건을 만족하는 확률과정을 말합니다. 

 

[ 정상시계열의 조건 (conditions for stationary time series)]

  (a) 평균이 일정하다. 

  (b) 분산이 존재하며 상수이다. 

  (c) 두 시점 사이의 자기공분산은 시차(time lag)에만 의존한다. 

 

만약 추세(trend)가 있거나 시간이 갈 수록 분산이 커지거나 작아지는 시계열 데이터라면 정상확률과정 시계열이 아닙니다. 추세가 있거나 분산이 달라지게 되면 ARIMA 모형을 적합하기 전에, 먼저 추세를 없애 주고 분산을 안정화 시켜주는 변환을 통해 정상확률과정으로 만들어주어야 합니다. 

 

stationary process 정상확률과정

 

정상확률과정 중에서 대표적인 모수적 확률모형은 ARMA(p, q) 로 표기되는 자기회귀이동평균과정(AutoRegressive Moving Average Process) 입니다. 아래의 Python 코드는 AR(1) 인 자기회귀과정을 생성해서 시각화해본 것입니다. 

 

AR(1) 과정:  z(t) = Phi * z(t-1) + a(t)

                  (이때 Phi 는 절대값이 1 미만이고, a(t) 는 평균이 0이고 분산이 sigma^2 인 백색잡음과정임.)

 

## Stationary Process
## exmaple: AR(1) process
z_0 = 0
rho = 0.65 # correlation b/w z(t) and z(t-1)
z_all = [z_0]

for i in range(200):
    z_t = rho*z_all[i] + np.random.normal(0, 0.1, 1)
    z_all.append(z_t)
    

## plotting
plt.plot(z_all)

## adding horizonal line at mean position 0
plt.axhline(0, 0, 200, color='red', linestyle='--', linewidth=2)
plt.show()

stationary process plot (example)

 

위의 1번의 백색잡음과정과 3번의 정상확률과정이 눈으로만 봐서는 구분이 잘 안가지요? 

백색잡음과정, 확률보행과정, 정상확률과장이 서로 관련이 있고, 얼핏 봐서는 비슷하기도 해서 참 헷갈립니다. 정확하게 이 셋을 구분해서 이해하시기 바래요. 

 

 

[ Reference ]

1. 박유성, 김기환 저 , "SAS/ETS를 이용한 시계열자료분석 I", 자유아카데미

2. "White Noise" from Wikipedia: https://en.wikipedia.org/wiki/White_noise

3. "Random Walk" from Wikipedia: https://en.wikipedia.org/wiki/Random_walk

4. "Stationary Process" from Wikipedia: https://en.wikipedia.org/wiki/Stationary_process 

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 주기성을 가지는 센서 시계열 데이터, 음향 데이터, 신호 처리 데이터 등에서 많이 사용되는 스펙트럼 분석 (spectral analysis) 에 대해서 소개하겠습니다. 

 

많은 시계열 데이터는 주기적인 패턴 (periodic behavior) 을 보입니다. 스펙트럼 분석 은 시계열 데이터에 내재되어 있는 주기성(periodicity) 을 찾아내는 분석기법으로서, spectral analysis, spectrum analysis, frequency domain analysis 등의 이름으로 불립니다. 

 

스펙트럼 분석을 할 때 우리는 먼저 시계열 데이터를 시간 영역(time domain)에서 주파수 영역(frequency domain)로 변환을 합니다. 시계열의 공분산스펙트럼 확률밀도함수 (spectral density function) 로 알려진 함수로 재표현될 수 있습니다. 스펙트럼 밀도함수는 시계열과 다른 주파수의 사인/코사인 파동 간의 제곱 상관계수로 구하는 주기도(periodogram) 를 사용하여 추정할 수 있습니다. (Venables & Ripley, 2002)

 

아래의 그림은 시계열(time series)과 스펙트럼(spectrum, frequency) 간의 관계를 직관적으로 이해할 수 있게 도와줍니다. 

 

[ Time Series vs. Spectrum/ Frequency, Fourier Transform(FT) vs. Inverse FT ]

 

time-series vs. spectrum

 

 

참고로, 주파수(frequency)는 단위 시간 당 반복되는 이벤트가 발생하는 회수를 의미합니다. 주파수를 측정하는데 사용되는 단위인 Hz (Hertz, 헤르쯔) 는 전자기파의 존재를 처음으로 증명한 Heinrich Rudolf Hertz 의 이름에서 따온 것인데요, 가령 2 Hz 는 1초에 2번의 반복 사이클이 발생하는 시계열의 경우, 4 Hz 는 1초에 4번의 반복 사이클이 발행하는 시계열의 주파수에 해당합니다. (아래 그림 예시 참조)

 

frequency (* https://rfriend.tistory.com)

 

 

 

주기성을 띠는 시계열 데이터나 공간 데이터로 부터 주파수 (spectrum, frequency) 를 계산하거나, 반대로 주파수로 부터 원래의 데이터로 변환할 때 고속 푸리에 변환 (Fast Fourier Transform, 줄여서 FFT) 알고리즘을 많이 사용합니다.  시계열을 주파수로 변환하는 과정 자체를 Fourier Transform 이라고 하며, 이산 푸리에 변환 (Discrete Fourier Transform, DFT) 은 시계열 값을 다른 주파수의 구성요소들로 분해함으로써 구할 수 있습니다. 지금까지 가장 널리 사용되는 고속 푸리에 변환 (FFT) 알고리즘은 Cooley-Tukey 알고리즘 입니다.  이것은 분할-정복 알고리즘 (divide-and-conqer algorithm) 으로서, 반복적으로(recursively) 원래의 시계열 데이터를 훨씬 적은 개수의 주파수로 분할합니다. 

 

[ Fast Fourier Transform(FFT) vs. Inverse FFT]

FFT vs. Inverse FFT

 

 

이제 Python 의 spectrum 모듈을 사용해서 스펙트럼 분석을 해보겠습니다. 터미널에서 $ pip install spectrum 으로 먼저 spectrum 모듈을 설치해줍니다. 

 

## installment of spectrum python library at a terminal
$ pip install spectrum

 

 

spectrum 모듈에서 Periodogramm data_cosine 메소드를 불러옵니다.

예제로 사용할 샘플데이터로서, data_cosine() 메소드를 사용해서 200 Hz 주파수 (freq=200) 의 주기성을 띠는 코사인 신호(cosine signal)에 백색 잡음 (amplitude 0.1) 도 포함되어 있는 이산형 샘플 데이터 1,024개 (N=1024, sampling=1024) 를 생성해보겠습니다. 

 

## -- importing modules
from spectrum import Periodogram, data_cosine

## generating a toy data sets
## data contains a cosine signal with a frequency of 200Hz 
## buried in white noise (amplitude 0.1). 
## The data has a length of N=1024 and the sampling is 1024Hz.

data = data_cosine(N=1024, A=0.1, sampling=1024, freq=200)

## plotting the generated data
plt.rcParams['figure.figsize'] = [18, 8]
plt.plot(data)
plt.show()

 

 

위의 시계열 그래프가 1,024개 데이터 포인트가 모두 포함되다보니 코사인 신호의 주기성이 눈에 잘 안보이는데요, 앞 부분의 50개만 가져다가 시도표로 그려보면 아래처럼 코사인 파동 모양의 주기성을 띠는 데이터임을 눈으로 확인할 수 있습니다. 

 

plt.rcParams['figure.figsize'] = [18, 8]
plt.plot(data[:50])
plt.show()

cosine signal waves

 

 

spectrum 모듈의 Periodogram() 메소드에 위에서 생성한 샘플 데이터를 적용해서 periodogram 객체를 만든 후에 run() 함수를 실행시켜서 스펙트럼을 추정할 수 있고, p.plot() 으로 간단하게 추정한 스펙트럼을 시각화할 수 있습니다. 

 

## the Power Spectrum Estimation method provided in spectrum.
## creating an object Periodogram. 
p = Periodogram(data, sampling=1024)

## running the spectrum estimation
p.run() # or p()

p.plot(marker='o')

 

 

스펙트럼 분석 결과 중 개수(N), 주파수의 세기 (power of spectrum dendity) 속성값을 확인할 수 있습니다.  

 

## some information are stored and can be retrieved later on.
p.N
[Out] 1024

p.psd

# [Out] array([1.06569152e-04, 1.74100855e-03, 3.93335837e-03, 1.10485152e-03,
#        8.28499278e-03, 7.64848301e-04, 2.94332240e-03, 7.35029918e-03,
#        1.21273520e-03, 3.72342384e-05, 1.87058542e-03, 4.35817479e-03,
#        6.60803332e-04, 3.01094523e-03, 3.14914528e-03, 3.40040266e-03,
#        1.70112710e-04, 4.04579597e-04, 1.73102721e-03, 5.17048610e-03,
#        5.89537834e-03, 3.43004928e-03, 2.60894164e-03, 2.29556960e-03,
#        6.23836725e-03, 7.35588289e-03, 3.88793797e-03, 2.25691304e-03,
#        ... 중간 생략함 ...
#        2.57852156e-03, 1.23483225e-03, 1.49139300e-03, 9.09082081e-04,
#        5.24999741e-04])

 

 

우리가 위에서 data_cosine(N=1024, A=0.1, sampling=1024, freq=200) 메소드를 사용해서 200 Hz 의 데이터를 생성하였기 때문에, 이 샘플 데이터로부터 주파수를 분해하면 200 Hz 가 되어야 할텐데요, 스펙트럼 분석 결과도 역시 200 Hz 에서 스파이크가 엄청 높게 튀는군요. 

 

plt.rcParams['figure.figsize'] = [10, 6]
plt.plot(p.psd)
plt.show()

 

 

* R 을 이용한 스펙트럼 분석은 https://rfriend.tistory.com/64 를 참고하시기 바랍니다. 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

시계열 데이터를 분석할 때 꼭 확인하고 처리해야 하게 있는데요, 바로 결측값 여부 확인과 결측값 처리입니다. 

 

시계열 데이터의 결측값을 처리하는 방법에는

   (1) 보간 (Interpolation)

   (2) 이전 값 또는 다음 값 이용 (previous/next value)

   (3) 이동 평균 (Moving average)

등의 여러가지 방법이 있습니다. 

 

 

[ 시계열 데이터 결측값 처리 방법 (How to handle the time series missing data) ]

 

 

아래의 보간(Interpolation)에 대한 내용은 Wikipedia 의 내용을 번역하여 소개합니다. 

 

데이터 분석의 수학 분야에서는 "보간법(Interpolation)을 이미 알려진 데이터 포인트들의 이산형 집합의 범위에 기반해서 새로운 데이터 포인트들을 만들거나 찾는 추정(estimation)의 한 유형"으로 봅니다. 

 

공학과 과학 분야에서는 종종 샘플링이나 실험을 통해서 많은 수의 데이터 포인트들을 획득하는데요, 이들 데이터는 어떤 함수(a function)의 값이나 독립변수(independent variable)의 제한적인 수의 값을 표현한 것입니다. 종종 독립변수의 중간 사이의 값을 위한 함수의 값을 추정(estimate the value of that function for an intermediat value of the independent variable)하는 보간이 필요합니다.  

 

밀접하게 관련된 문제로서 복잡한 함수를 간단한 함수로 근사하게 추정(the approximation of a complicated function by a simple function)하는 것이 있습니다. 어떤 주어진 함수의 공식이 알려져있지만, 너무 복잡해서 효율적으로 평가하기가 어렵다고 가정해봅시다. 원래의 함수로부터 적은 수의 새로운 데이터 포인트는 원래의 값과 상당히 근접한 간단한 함수를 생성해서 보간할 수 있습니다. 단순성(simplicity)으로부터 얻을 수 있는 이득이 보간에 의한 오차라는 손실보다 크고, 연산 프로세스면서도 더 좋은 성능(better performance in calculation process)을 낼 수도 있습니다.   

 

 

이번 포스팅에서는 Python scipy 모듈을 이용해서 시계열 데이터 결측값을 보간(Interpolation)하는 방법을 소개하겠습니다. 

 

1. 이전 값/ 이후 값을 이용하여 결측값 채우기 (Imputation using the previous/next values)

2. Piecewise Constant Interpolation

3. 선형 보간법 (Linear Interpolation)

4. 스플라인 보간법 (Spline Interpolation)

 

 

[ Python scipy 모듈을 이용한 결측값 보간 (Interpolation using Python scipy module)  ]

 

 

먼저 '0.5'로 동일한 간격을 가지는 x 값들에 대한 사인 함수 (sine function) 의 y값을 계산해서  예제 데이터로 사용하겠습니다. 아래 예졔의 점과 점 사이의 값들이 비어있는 결측값이라고 간주하고, 이들 값을 채워보겠습니다. 

 

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

## generating the original data with missing values
x = np.arange(0, 4*np.pi, 0.5)
y = np.sin(x)

plt.plot(x, y, "o")
plt.show()

original data with missing values

 

 

1. 이전 값/ 다음 값을 이용하여 결측값 채우기 (Imputation using the previous/next values)

 

데이터 포인트 사이의 값을 채우는 가장 간단한 방법은 이전 값(previous value) 나 또는 다음 값(next value)을 이용하는 것입니다. 함수를 추정하는 절차가 필요없으므로 연산 상 부담이 적지만, 데이터 추정 오차는 단점이 될 수 있습니다. 

 

## Interpolation using the previous value
f_prev = interpolate.interp1d(
    x, y, kind='previous') # next
y_new_prev = f_prev(xnew)

plt.plot(x, y, "o", xnew, y_new_prev, '-')
plt.show()

interpolation using the previous value

 

 

 

2. Piecewise Constant Interpolation

 

위 1번의 이전 값 또는 다음 값을 이용한 사이값 채우기를 합쳐놓은 방법입니다. Piecewise Constant Interpolation은 특정 데이터 포인트를 기준으로 가장 가까운 값 (nearest value) 을 가져다가 사이값을 보간합니다. ("최근접 이웃 보간"이라고도 함)

 

간단한 문제에서는 아래 3번에서 소개하는 Linear Interpolation 이 주로 사용되고, Piecewise Constant Interpolation 은 잘 사용되지 않는 편입니다. 하지만 다차원의 다변량 보간 (in higher-dimensional multivariate interpolation)의 경우, 속도와 단순성(speed and simplicity) 측면에서 선호하는 선택이 될 수 있습니다. 

 

## Piecewise Constant Interpolation
f_nearest = interpolate.interp1d(
    x, y, kind='nearest')

y_new_nearest = f_nearest(xnew)

plt.plot(x, y, "o", xnew, y_new_nearest)
plt.show()

Piecewise constant interpolation

 

 

 

3. Linear Interpolation

 

선형 보간법은 가장 쉬운 보간법 중의 하나로서, 연산이 빠르고 쉽습니다. 하지만 추정값이 정확한 편은 아니며, 데이터 포인트 Xk 에서 미분 가능하지 않다는 단점도 있습니다. 

 

일반적으로, 선형 보간법은 두 개의 데이터 포인트, 가령 (Xa, Ya)와 (Xb, Yb), 를 사용해서 다음의 공식으로 두 값 사이의 값을 보간합니다. 

 

Y = Ya + (Yb - Ya) * (X - Xa) / (Xb- Xa)    at the point (x, y)

 

## Linear Interpolation
f_linear = interpolate.interp1d(
    x, y, kind='linear')

y_new_linear = f_linear(xnew)

plt.plot(x, y, "o", xnew, y_new_linear, '-')
plt.show()

Linear interpolation

 

 

 

4. Spline Interpolation

 

다항식 보간법(Polynomial Interpolation)은 선형 보간법을 일반화(generalization of linear interpolation)한 것입니다. 선형 보간법에서는 선형 함수를 사용했다면, 다항식 보간법에서는 더 높은 차수의 다항식 함수를 사용해서 보간하는 것으로 대체한 것입니다. 

일반적으로, 만약 우리가 n개의 데이터 포인트를 가지고 있다면 모든 데이터 포인트를 통과하는 n-1 차수의 다차항 함수가 존재합니다. 보간 오차는 데이터 포인트 간의 거리의 n 차승에 비례(interpolation error is proportional to the distance between the data points to the power n)하며, 다차항 함수는 미분가능합니다. 따라서 선형 보간법의 대부분의 문제를 다항식 보간법은 극복합니다. 하지만 다항식 보간법은 선형 보간법에 비해 복잡하고 연산에 많은 비용이 소요됩니다. 그리고 끝 점(end point) 에서는 진동하면서 변동성이 큰 값을 추정하는 문제가 있습니다. 

 

스플라인 보간법은 각 데이터 포인트 구간별로 낮은 수준의 다항식 보간을 사용 (Spline interpolation uses low-degree polynomials in each of the intervals) 합니다. 그리고 이들이 함께 부드럽게 연결되어서 적합될 수 있도록 다항식 항목을 선택(, and chooses the polynomial pieces such that they fit smoothly together)합니다. 이렇게 적합된 함수를 스플라인(Spline) 이라고 합니다. 

 

스플라인 보간법(Spline Interpolation)은 다항식 보간법의 장점은 살리고 단점은 피해간 보간법입니다. 스플라인 보간법은 다항식 보간법처럼 선형 보간법보다 보간 오차가 더 작은 반면에, 고차항의 다항식 보간법보다는 보간 함수가 부드럽고 평가하기가 쉽습니다.  

 

## Spline Interpolation
f_quadr = interpolate.interp1d(
    x, y, kind='quadratic') # cubic

y_new_quadr = f_quadr(xnew)

plt.plot(x, y, "o", xnew, y_new_quadr)
plt.show()

Polynomial interpolation

 

 

[ Reference ]

1. 보간법(interpolation): https://en.wikipedia.org/wiki/Interpolation

2. scipy 모듈: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 시계열 패턴 별 지수평활법 모델을 설정하는 방법 (https://rfriend.tistory.com/511), 그리고 시계열 데이터 예측 모델의 성능, 모델 적합도 (goodness-of-fit of time series model)를 평가하는 지표 (https://rfriend.tistory.com/667) 에 대해서 소개하였습니다.  

 

이번 포스팅에서는 간단한 시계열 예제 데이터에 대해 앞서 소개한 이론들을 Python 코드로 실습하여 

(1) 지수평활법(Exponential Smoothing) 기법으로 시계열 예측 모델을 만들고, 

(2) 모델 적합도를 평가하는 여러 지표들로 성능을 평가해서 최적의 모델을 선택

해보겠습니다. 

 

 

예제로 사용할 시계열 데이터는 1991년 7월부터 2008년 6월까지 월 단위(monthly)로 집계된 약 판매량 시계열 데이터셋 (drug sales time series dataset) 입니다. 

 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

## getting drug sales dataset
file_path = 'https://raw.githubusercontent.com/selva86/datasets/master/a10.csv'
df = pd.read_csv(file_path, 
                 parse_dates=['date'], 
                 index_col='date')
                 
df.head(12)

# [Out]
#             value
# date	
# 1991-07-01	3.526591
# 1991-08-01	3.180891
# 1991-09-01	3.252221
# 1991-10-01	3.611003
# 1991-11-01	3.565869
# 1991-12-01	4.306371
# 1992-01-01	5.088335
# 1992-02-01	2.814520
# 1992-03-01	2.985811
# 1992-04-01	3.204780
# 1992-05-01	3.127578
# 1992-06-01	3.270523


## time series plot
df.plot(figsize=[12, 8])
plt.title('Drug Sales', fontsize=14)
plt.show()

drug sales time series plot

 

위의 약 판매량 시도표 (time series plot of drug sales) 를 보면 세 가지를 알 수 있습니다. 

 

(a) 선형 추세가 존재합니다. (linear or quadratic trend)

(b) 1년 단위의 계절성이 존재합니다. (seasonality of a year, ie. 12 months)

(c) 분산이 시간의 흐름에 따라 점점 커지고 있습니다. (increasing variation)

 

우리가 관찰한 이들 세가지의 시계열 데이터의 패턴으로 부터 "Multiplicative Winters' method with Linear / Quadratic Trend" 모델이 적합할 것이라고 판단해볼 수 있습니다. 

* 참조: 시계열 패턴 별 지수평활법 모델을 설정하는 방법 (https://rfriend.tistory.com/511)

 

그럼, 실제로 다양한 지수평활법으로 위의 약 판매량 시계열 데이터에 대한 예측 모델을 만들어 예측을 해보고, 다양한 평가지표로 모델 간 비교를 해봐서 가장 우수한 모델이 무엇인지, 과연 앞서 말한 "Multiplicative Winters' method with Linear Trend" 모델이 가장 우수할 것인지 확인해보겠습니다. 보는게 믿는거 아니겠습니까!

 

 

이제 월 단위로 시간 순서대로 정렬된 시계열 데이터를 제일 뒤 (최근)의 12개월치는 시계열 예측 모델의 성능을 평가할 때 사용할 test set으로 따로 떼어놓고, 나머지 데이터는 시계열 예측모형 훈련에 사용할 training set 으로 분할 (splitting between the training and the test set)을 해보겠습니다. 

 

## split between the training and the test data sets. 
## The last 12 periods form the test data
df_train = df.iloc[:-12]
df_test = df.iloc[-12:]

 

 

(1) 지수평활법 기법으로 시계열 예측모형 적합하기
       (training the exponential smoothing models in Python)

 

Python의 statsmodels 라이브러리는 시계열 데이터 분석을 위한 API를 제공합니다. 아래의 코드는 순서대로 

 

(a) Simple Exponential Smoothing 

-- 추세 요소의 유형 (type of trend components)에 따라서 

(b) Holt's (trend=None)

(c) Exponential trend (exponential=True)

(d) Additive damped trend (damped_trend=True)

(e) Multiplicative damped trend (exponential=True, damped_trend=True)

 

로 구분할 수 있습니다. 

 

## exponential smoothing in Python
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt

# Simple Exponential Smoothing
fit1 = SimpleExpSmoothing(df_train, 
                          initialization_method="estimated").fit()

# Trend
fit2 = Holt(df_train, 
            initialization_method="estimated").fit()

# Exponential trend
fit3 = Holt(df_train,
            exponential=True, 
            initialization_method="estimated").fit()

# Additive damped trend
fit4 = Holt(df_train,
            damped_trend=True, 
            initialization_method="estimated").fit()

# Multiplicative damped trend
fit5 = Holt(df_train,
            exponential=True, 
            damped_trend=True, 
            initialization_method="estimated").fit()

 

 

 

statsmodels.tsa.api 메소드를 사용해서 적합한 시계열 데이터 예측 모델에 대해서는 model_object.summary() 메소드를 사용해서 적합 결과를 조회할 수 있습니다. 그리고 model_object.params['parameter_name'] 의 방식으로 지수평활법 시계열 적합 모델의 속성값(attributes of the fitted exponential smoothing results)에 대해서 조회할 수 있습니다. 

 

## accessing the results of SimpleExpSmoothing Model
print(fit1.summary())

#                        SimpleExpSmoothing Model Results                       
# ==============================================================================
# Dep. Variable:                  value   No. Observations:                  192
# Model:             SimpleExpSmoothing   SSE                            692.196
# Optimized:                       True   AIC                            250.216
# Trend:                           None   BIC                            256.731
# Seasonal:                        None   AICC                           250.430
# Seasonal Periods:                None   Date:                 Thu, 29 Jul 2021
# Box-Cox:                        False   Time:                         19:26:46
# Box-Cox Coeff.:                  None                                         
# ==============================================================================
#                        coeff                 code              optimized      
# ------------------------------------------------------------------------------
# smoothing_level            0.3062819                alpha                 True
# initial_level              3.4872869                  l.0                 True
# ------------------------------------------------------------------------------


## getting all results by models
params = ['smoothing_level', 'smoothing_trend', 'damping_trend', 'initial_level', 'initial_trend']
results=pd.DataFrame(index=[r"$\alpha$",r"$\beta$",r"$\phi$",r"$l_0$","$b_0$","SSE"],
                     columns=['SES', "Holt's","Exponential", "Additive", "Multiplicative"])

results["SES"] =            [fit1.params[p] for p in params] + [fit1.sse]
results["Holt's"] =         [fit2.params[p] for p in params] + [fit2.sse]
results["Exponential"] =    [fit3.params[p] for p in params] + [fit3.sse]
results["Additive"] =       [fit4.params[p] for p in params] + [fit4.sse]
results["Multiplicative"] = [fit5.params[p] for p in params] + [fit5.sse]

results
# [Out]
#      SES	Holt's	Exponential	Additive	Multiplicative
# 𝛼 	0.306282	0.053290	1.655513e-08	0.064947	0.040053
# 𝛽 	NaN	0.053290	1.187778e-08	0.064947	0.040053
# 𝜙 	NaN	NaN	NaN	0.995000	0.995000
# 𝑙0 	3.487287	3.279317	3.681273e+00	3.234974	3.382898
# 𝑏0 	NaN	0.045952	1.009056e+00	0.052515	1.016242
# SSE  692.195826	631.551871	5.823167e+02	641.058197	621.344891

 

 

계절성(seasonality)을 포함한 지수평활법에 대해서는 statsmodels.tsa.holtwinters 내의 메소드를 사용해서 시계열 예측 모델을 적합해보겠습니다. 

 

seasonal_periods 에는 계절 주기를 입력해주는데요, 이번 예제에서는 월 단위 집계 데이터에서 1년 단위로 계절성이 나타나고 있으므로 seasonal_periods = 12 를 입력해주었습니다. 

trend = 'add' 로 일단 고정을 했으며, 대신에 계절성의 경우 (f) seasonal = 'add' (or 'additive'), (g) seasonal = 'mul' (or 'multiplicative') 로 구분해서 모델을 각각 적합해보았습니다. 

 

위의 시계열 시도표를 보면 시간이 흐름에 따라 분산이 점점 커지고 있으므로  seasonal = 'mul' (or 'multiplicative') 가 더 적합할 것이라고 우리는 예상할 수 있습니다. 실제 그런지 데이터로 확인해보겠습니다. 

 

## Holt's Winters's method for time series data with Seasonality
from statsmodels.tsa.holtwinters import ExponentialSmoothing as HWES

# additive model for fixed seasonal variation
fit6 = HWES(df_train, 
             seasonal_periods=12, 
             trend='add', 
             seasonal='add').fit(optimized=True, use_brute=True)

# multiplicative model for increasing seasonal variation
fit7 = HWES(df_train, 
             seasonal_periods=12, 
             trend='add', 
             seasonal='mul').fit(optimized=True, use_brute=True)

 

 

 

 

(2) 지수평활법 적합 모델로 예측하여 모델 간 성능 평가/비교하고 최적 모델 선택하기

     (forecasting using the exponential smoothing models, model evaluation/ model comparison and best model selection)

 

 

(2-1) 향후 12개월 약 판매량 예측하기 (forecasting the drug sales for the next 12 months) 

 

forecast() 함수를 사용하며, 괄호 안에는 예측하고자 하는 기간을 정수로 입력해줍니다. 

 

## forecasting for 12 months
forecast_1 = fit1.forecast(12)
forecast_2 = fit2.forecast(12)
forecast_3 = fit3.forecast(12)
forecast_4 = fit4.forecast(12)
forecast_5 = fit5.forecast(12)
forecast_6 = fit6.forecast(12)
forecast_7 = fit7.forecast(12


y_test = df_test['value']

t_p = pd.DataFrame({'test': y_test, 
                    'f1': forecast_1, 
                    'f2': forecast_2, 
                    'f3': forecast_3, 
                    'f4': forecast_4, 
                    'f5': forecast_5, 
                    'f6': forecast_6, 
                    'f7': forecast_7})
                    
                    
print(t_p)

#                  test        f1         f2         f3         f4         f5  \
# 2007-07-01  21.834890  20.15245  20.350267  20.970823  20.454211  20.261185   
# 2007-08-01  23.930204  20.15245  20.511013  21.160728  20.616026  20.421627   
# 2007-09-01  22.930357  20.15245  20.671760  21.352352  20.777032  20.582528   
# 2007-10-01  23.263340  20.15245  20.832507  21.545711  20.937233  20.743882   
# 2007-11-01  25.250030  20.15245  20.993254  21.740821  21.096633  20.905686   
# 2007-12-01  25.806090  20.15245  21.154000  21.937699  21.255236  21.067932   
# 2008-01-01  29.665356  20.15245  21.314747  22.136359  21.413046  21.230618   
# 2008-02-01  21.654285  20.15245  21.475494  22.336818  21.570067  21.393736   
# 2008-03-01  18.264945  20.15245  21.636241  22.539092  21.726303  21.557283   
# 2008-04-01  23.107677  20.15245  21.796987  22.743198  21.881758  21.721254   
# 2008-05-01  22.912510  20.15245  21.957734  22.949152  22.036435  21.885642   
# 2008-06-01  19.431740  20.15245  22.118481  23.156972  22.190339  22.050443   

#                    f6         f7  
# 2007-07-01  20.882394  21.580197  
# 2007-08-01  22.592307  22.164134  
# 2007-09-01  21.388974  22.136953  
# 2007-10-01  25.178594  24.202879  
# 2007-11-01  26.776556  25.322929  
# 2007-12-01  26.424404  27.549045  
# 2008-01-01  30.468546  31.194773  
# 2008-02-01  19.065461  18.181477  
# 2008-03-01  21.834584  21.034783  
# 2008-04-01  19.061426  20.085947  
# 2008-05-01  23.452042  23.076423  
# 2008-06-01  22.353486  22.598155

 

 

 

(2-2) 모델 간 모델적합도/성능 평가하고 비교하기

 

Python의 scikit-learn 라이브러리의 metrics 메소드에 이미 정의되어 있는 MSE, MAE 지표는 가져다가 쓰구요, 없는 지표는 모델 적합도 (goodness-of-fit of time series model)를 평가하는 지표 (https://rfriend.tistory.com/667) 포스팅에서 소개했던 수식대로 Python 사용자 정의 함수(User Defined Funtion)을 만들어서 사용하겠습니다. 

 

아래의 모델 성능 평가 지표 중에서 AIC, SBC, APC, Adjsted-R2 지표는 모델 적합에 사용된 모수의 개수 (the number of parameters) 를 알아야지 계산할 수 있으며, 모수의 개수는 지수평활법 모델 종류별로 다르므로, 적합된 모델 객체로부터 모수 개수를 세는 사용자 정의함수를 먼저 정의해보겠습니다. 

 

## UDF for counting the number of parameters in model
def num_params(model):
    n_params = 0

    for p in list(model.params.values()):
        if isinstance(p, np.ndarray):
            n_params += len(p)
            #print(p)
        elif p in [np.nan, False, None]:
            pass
        elif np.isnan(float(p)):
            pass
        else:
            n_params += 1
            #print(p)
    
    return n_params
    

num_params(fit1)
[Out] 2

num_params(fit7)
[Out] 17

 

 

이제 각 지표별로 Python을 사용해서 지수평활법 모델별로 성능을 평가하는 사용자 정의함수를 정의해보겠습니다. MSE와 MAE는 scikit-learn.metrics 에 있는 메소드를 가져다가 사용하며, 나머지는 공식에 맞추어서 정의해주었습니다.

(혹시 sklearn 라이브러리에 모델 성능 평가 지표 추가로 더 있으면 그거 가져다 사용하시면 됩니다. 저는 sklearn 튜토리얼 찾아보고 없는거 같아서 아래처럼 정의해봤어요. 혹시 코드 잘못된거 있으면 댓글에 남겨주시면 감사하겠습니다.)

 

## number of observations in training set
T = df_train.shape[0]
print(T)
[Out] 192


## evaluation metrics
from sklearn.metrics import mean_squared_error as MSE
from sklearn.metrics import mean_absolute_error as MAE

# Mean Absolute Percentage Error
def SSE(y_test, y_pred):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    return np.sum((y_test - y_pred)**2)

def ME(y_test, y_pred):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    return np.mean(y_test - y_pred)

def RMSE(y_test, y_pred):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    return np.sqrt(np.mean((y_test - y_pred)**2))   
    #return np.sqrt(MSE(y_test - y_pred))

def MPE(y_test, y_pred): 
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    return np.mean((y_test - y_pred) / y_test) * 100

def MAPE(y_test, y_pred): 
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    return np.mean(np.abs((y_test - y_pred) / y_test)) * 100

def AIC(y_test, y_pred, T, model):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    sse = np.sum((y_test - y_pred)**2)
    #T = len(y_train) # number of observations
    k = num_params(model) # number of parameters
    return T * np.log(sse/T) + 2*k

def SBC(y_test, y_pred, T, model):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    sse = np.sum((y_test - y_pred)**2)
    #T = len(y_train) # number of observations
    k = num_params(model) # number of parameters
    return T * np.log(sse/T) + k * np.log(T)

def APC(y_test, y_pred, T, model):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    sse = np.sum((y_test - y_pred)**2)
    #T = len(y_train) # number of observations
    k = num_params(model) # number of parameters
    return ((T+k)/(T-k)) * sse / T

def ADJ_R2(y_test, y_pred, T, model):
    y_test, y_pred = np.array(y_test), np.array(y_pred)
    sst = np.sum((y_test - np.mean(y_test))**2)
    sse = np.sum((y_test - y_pred)**2)
    #T = len(y_train) # number of observations
    k = num_params(model) # number of parameters
    r2 = 1 - sse/sst
    return 1 - ((T - 1)/(T - k)) * (1 - r2)
    

## Combining all metrics together
def eval_all(y_test, y_pred, T, model):
    sse = SSE(y_test, y_pred)
    mse = MSE(y_test, y_pred)
    rmse = RMSE(y_test, y_pred)
    me = ME(y_test, y_pred)
    mae = MAE(y_test, y_pred)
    mpe = MPE(y_test, y_pred)
    mape = MAPE(y_test, y_pred)
    aic = AIC(y_test, y_pred, T, model)
    sbc = SBC(y_test, y_pred, T, model)
    apc = APC(y_test, y_pred, T, model)
    adj_r2 = ADJ_R2(y_test, y_pred, T, model)
    
    return [sse, mse, rmse, me, mae, mpe, mape, aic, sbc, apc, adj_r2]

 

 

 

위에서 정의한 시계열 예측 모델 성능평가 지표 사용자 정의 함수를 모델 성능 평가를 위해 따로 떼어놓았던 'test set' 에 적용하여 각 지수평활법 모델별/ 모델 성능평가 지표 별로 비교를 해보겠습니다. 

모델 성능평가 지표 중에서 실제값과 예측값의 차이인 잔차(residual)를 기반으로 한 MSE, RMSE, ME, MAE, MPE, MAPE, AIC, SBC, APC 는 낮을 수록 상대적으로 더 좋은 모델이라고 평가할 수 있으며, 실제값과 평균값의 차이 중에서 시계열 예측 모델이 설명할 수 있는 부분의 비중 (모델의 설명력) 을 나타내는 Adjusted-R2 는 값이 높을 수록 상대적으로 더 좋은 모델이라고 평가합니다. (Ajd.-R2 가 음수가 나온 값이 몇 개 보이는데요, 이는 예측 모델이 불량이어서 예측값이 실제값과 차이가 큼에 따라 SSE가 SST보다 크게 되었기 때문입니다.)

 

eval_all_df = pd.DataFrame(
    {'SES': eval_all(y_test, forecast_1, T, fit1), 
    "Holt's": eval_all(y_test, forecast_2, T, fit2), 
    'Exponential': eval_all(y_test, forecast_3, T, fit3), 
    'Trend_Add': eval_all(y_test, forecast_4, T, fit4), 
    'Trend_Mult': eval_all(y_test, forecast_5, T, fit5), 
    'Trend_Season_Add': eval_all(y_test, forecast_6, T, fit6), 
    'Trend_Season_Mult': eval_all(y_test, forecast_7, T, fit7)}
    , index=['SSE', 'MSE', 'RMSE', 'ME', 'MAE', 'MPE', 'MAPE', 'AIC', 'SBC', 'APC', 'Adj_R2'])


print(eval_all_df)
#                SES      Holt's  Exponential   Trend_Add  Trend_Mult  \
# SSE     205.629692  155.645609   130.567477  150.867568  159.856508   
# MSE      17.135808   12.970467    10.880623   12.572297   13.321376   
# RMSE      4.139542    3.601454     3.298579    3.545744    3.649846   
# ME        3.018502    1.936578     1.123475    1.841425    2.019134   
# MAE       3.453205    2.946251     2.576567    2.878085    3.004307   
# MPE      11.739201    6.900165     3.313095    6.486233    7.259990   
# MAPE     14.079695   12.280881    10.960354   12.010763   12.510302   
# AIC      17.167661  -32.303428   -66.036202  -36.289847  -25.178006   
# SBC      23.682652  -19.273447   -53.006220  -20.002370   -8.890530   
# APC       1.093535    0.845150     0.708977    0.827788    0.877109   
# Adj_R2   -1.146688   -0.642161    -0.377571   -0.600262   -0.695608   

#         Trend_Season_Add  Trend_Season_Mult  
# SSE            56.743149          48.994715  
# MSE             4.728596           4.082893  
# RMSE            2.174533           2.020617  
# ME             -0.118946          -0.089689  
# MAE             1.863421           1.641140  
# MPE            -0.847971          -0.607618  
# MAPE            8.538100           7.461092  
# AIC          -200.040404        -228.230323  
# SBC          -144.662983        -172.852902  
# APC             0.352956           0.304759  
# Adj_R2          0.356850           0.444674

 

 

시계열 모델 성능 지표별로 보면 전반적으로 7번째 지수평활법 시계열 예측 모델인 "Multiplicative Winters' method with Linear Trend" 모델이 상대적으로 가장 우수함을 알 수 있습니다. 역시 시도표를 보고서 우리가 예상했던대로 결과가 나왔습니다. 

MAPE (Mean Absolute Percentage Error)  지표에 대해서만 옆으로 누운 막대그래프로 모델 간 성능을 비교해서 시각화를 해보면 아래와 같습니다. (MAPE 값이 작을수록 상대적으로 더 우수한 모델로 해석함.)

 

# horizontal bar chart
eval_all_df.loc['MAPE', :].plot(kind='barh', figsize=[8, 6])
plt.title('Mean Absolute Percentage Error (MAPE)', fontsize=16)
plt.show()

Mean Absolute Percentage Error (MAPE)

 

 

아래의 시도표는 실제값과 예측값을 시간의 흐름에 따라 선그래프로 시각화해 본 것입니다. 

예상했던 대로 4번째의 "1차 선형 추세만 반영하고 계절성은 없는 이중 지수 평활법 모델" 보다는, 7번째의 "선형추세+계절성까지 반영한 multiplicative holt-winters method exponential smoothing with linear trend 지수평활법 모델" 이 더 실제값을 잘 모델링하여 근접하게 예측하고 있음을 눈으로 확인할 수 있습니다. 

 

# 1차 선형 추세는 있고 계절성은 없는 이중 지수 평활법
plt.rcParams['figure.figsize']=[12, 8]
past, = plt.plot(df_train.index, df_train, 'b.-', label='Sales History')
test, = plt.plot(df_test.index, df_test, 'r.-', label='y_test')
pred, = plt.plot(df_test.index, forecast_4, 'y.-', label='y_pred')
plt.title('Two Parameter Exponential Smoothing', fontsize=14)
plt.legend(handles=[past, test, pred])
plt.show()

 

 

 

# 1차 선형 추세와 확산계절변동이 있는 승법 윈터스 지수평활법
plt.rcParams['figure.figsize']=[12, 8]
past, = plt.plot(df_train.index, df_train, 'b.-', label='Sales History')
test, = plt.plot(df_test.index, df_test, 'r.-', label='y_test')
pred, = plt.plot(df_test.index, forecast_7, 'y.-', label='y_pred')
plt.title('Multiplicative Winters Method Exponential Smoothing with Linear Trend', fontsize=14)
plt.legend(handles=[past, test, pred])
plt.show()

 

 

적합된 statsmodels.tsa 시계열 예측 모델에 대해 summary() 메소드를 사용해서 적합 결과를 조회할 수 있으며, model.params['param_name'] 으로 적합된 모델의 속성값에 접근할 수 있습니다. 

 

#print out the training summary
print(fit7.summary())

#                        ExponentialSmoothing Model Results                       
# ================================================================================
# Dep. Variable:                    value   No. Observations:                  192
# Model:             ExponentialSmoothing   SSE                             75.264
# Optimized:                         True   AIC                           -147.808
# Trend:                         Additive   BIC                            -95.688
# Seasonal:                Multiplicative   AICC                          -143.854
# Seasonal Periods:                    12   Date:                 Thu, 29 Jul 2021
# Box-Cox:                          False   Time:                         19:40:52
# Box-Cox Coeff.:                    None                                         
# =================================================================================
#                           coeff                 code              optimized      
# ---------------------------------------------------------------------------------
# smoothing_level               0.2588475                alpha                 True
# smoothing_trend               0.0511336                 beta                 True
# smoothing_seasonal           3.3446e-16                gamma                 True
# initial_level                 2.1420252                  l.0                 True
# initial_trend                 0.0466926                  b.0                 True
# initial_seasons.0             1.4468847                  s.0                 True
# initial_seasons.1             1.4713123                  s.1                 True
# initial_seasons.2             1.4550909                  s.2                 True
# initial_seasons.3             1.5754307                  s.3                 True
# initial_seasons.4             1.6324776                  s.4                 True
# initial_seasons.5             1.7590615                  s.5                 True
# initial_seasons.6             1.9730448                  s.6                 True
# initial_seasons.7             1.1392096                  s.7                 True
# initial_seasons.8             1.3057795                  s.8                 True
# initial_seasons.9             1.2354317                  s.9                 True
# initial_seasons.10            1.4064560                 s.10                 True
# initial_seasons.11            1.3648905                 s.11                 True
# ---------------------------------------------------------------------------------



print('Smoothing Level: %.3f' %(fit7.params['smoothing_level']))
print('Smoothing Trend: %.3f' %(fit7.params['smoothing_trend']))
print('Smoothing Seasonality: %.3f' %(fit7.params['smoothing_seasonal']))
print('Initial Level: %.3f' %(fit7.params['initial_level']))
print('Initial Trend: %.3f' %(fit7.params['initial_trend']))
print('Initial Seasons:', fit7.params['initial_seasons'])

# Smoothing Level: 0.259
# Smoothing Trend: 0.051
# Smoothing Seasonality: 0.000
# Initial Level: 2.142
# Initial Trend: 0.047
# Initial Seasons: [1.44688473 1.47131227 1.45509092 1.5754307  1.63247755 1.75906148
#  1.97304481 1.13920961 1.30577953 1.23543175 1.40645604 1.36489053]

 

[Reference]

* Exponential Smoothing
https://www.statsmodels.org/stable/examples/notebooks/generated/exponential_smoothing.html
* Holt-Winters 
https://www.statsmodels.org/devel/generated/statsmodels.tsa.holtwinters.ExponentialSmoothing.html

 

Exponential smoothing — statsmodels

Exponential smoothing Let us consider chapter 7 of the excellent treatise on the subject of Exponential Smoothing By Hyndman and Athanasopoulos [1]. We will work through all the examples in the chapter as they unfold. [1] [Hyndman, Rob J., and George Athan

www.statsmodels.org

 

statsmodels.tsa.holtwinters.ExponentialSmoothing — statsmodels

An dictionary containing bounds for the parameters in the model, excluding the initial values if estimated. The keys of the dictionary are the variable names, e.g., smoothing_level or initial_slope. The initial seasonal variables are labeled initial_season

www.statsmodels.org

 

* 시계열 패턴 별 지수평활법 모델을 설정하는 방법
   :
https://rfriend.tistory.com/511

* 시계열 데이터 예측 모델의 성능, 모델 적합도 (goodness-of-fit of time series model)를 평가하는 지표
   :
https://rfriend.tistory.com/667

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 시계열 패턴별 지수 평활법 (exponential smoothing by time series patterns) (https://rfriend.tistory.com/511) 에 대해서 소개하였습니다. 

 

이번 포스팅에서는 시계열 자료 예측 모형의 성능, 모델 적합도 (goodness-of-fit of the time series model) 를 평가할 수 있는 지표, 통계량을 알아보겠습니다. 

 

아래의 모델 적합도 평가 지표들의 리스트를 살펴보시면 선형회귀모형의 모델 적합도 평가 지표와 유사하다는 것을 알 수 있습니다.(실제값과 예측값의 차이 또는 설명비율에 기반한 성능 평가라는 측면에서는 동일하며, 회귀모형은 종속변수와 독립변수간 상관관계에 기반한 반면에 시계열 모형은 종속변수와 자기자신의 과거 데이터와의 자기상관관계에 기반한 다는것이 다른점 입니다.)

 

아래 평가지표들 중에서 전체 분산 중에서 모델이 설명할 수 있는 비율을 나타내는 수정결정계수는 통계량 값이 높을 수록 좋은 모델이며, 나머지 오차에 기반한 평가 지표(통계량)들은 값이 낮을 수록 상대적으로 더 좋은 모델이라고 평가를 합니다. (단, SST는 제외) 

 

이들 각 지표별로 좋은 모델 여부를 평가하는 절대 기준값(threshold)이 있는 것은 아니며, 여러개의 모델 간 성능을 상대적으로 비교 평가할 때 사용합니다. 

 

- 전체제곱합 (SST, total sum of square)

- 오차제곱합 (SSE, error sum of square)

- 평균오차제곱합 (MSE, mean square error)

- 제곱근 평균오차제곱합 (RMSE, root mean square error)

- 평균오차 (ME, mean error)

- 평균절대오차 (MAE, mean absolute error)

- 평균비율오차 (MPE, mean percentage error)

- 평균절대비율오차 (MAPE, mean absolute percentage error)

- 수정결정계수 (Adj. R-square) 

- AIC (Akaike's Information Criterion)

- SBC (Schwarz's Bayesian Criterion)

- APC (Amemiya's Prediction Criterion)

 

 

[ 시계열 데이터 예측 모형 적합도 평가 지표 (goodness-of-fit measures for time series prediction models) ]

goodness-of-fit measures for time series models

 

 

1. 전체제곱합 (SST, total sum of square)

 

total sum of square, SST

SST는 시계열 값에서 시계열의 전체 평균 값을 뺀 값으로, 시계열 예측 모델을 사용하지 않았을 때 (모든 모수들이 '0' 일 때) 의 오차 제곱 합입니다. 나중에 결정계수(R2), 수정결정계수(Adj. R2)를 계산할 때 사용됩니다. (SST 는 모델 성능 평가에서는 제외)

 

 

 

2. 오차제곱합 (SSE, error sum of square)

 

error sum of square, SSE

 

 

 

3. 평균오차제곱합 (MSE, mean square error), 제곱근 평균오차제곱합 (RMSE, root mean square error)

 

root mean squared error

 

MSE는 많은 통계 분석 패키지, 라이브러리에서 모델 훈련 시 비용함수(cost function) 또는 모델 성능 평가시 기본 설정 통계량으로 많이 사용합니다.  

 

RMSE (Root Mean Square Error) 는 MSE에 제곱근을 취해준 값으로서, SSE를 제곱해서 구한 MSE에 역으로 제곱근을 취해주어 척도를 원래 값의 단위로 맞추어 준 값입니다. 

 

 

 

4. 평균오차 (ME, mean error)

 

mean error

ME(Mean Error) 는 MAE(Mean Absolute Error)와 함께 해석하는 것이 필요합니다. 왜냐하면 큰 오차 값들이 존재한다고 하더라도 ME 값만 볼 경우 + 와 - 값이 서로 상쇄되어 매우 작은 값이 나올 수도 있기 때문입니다. 따라서 MAE 값을 통해 실제값과 예측값 간에 오차가 평균적으로 얼마나 큰지를 확인하고, ME 값의 부호를 통해 평균적으로 과다예측(부호가 '+'인 경우) 혹은 과소예측(부호가 '-'인 경우) 인지를 가늠해 볼 수 있습니다. 

 

 

 

5. 평균절대오차 (MAE, mean absolute error)

 

mean absolute error

 

 

 

6. 평균비율오차 (MPE, mean percentage error)

 

mean percentage error

위의 ME와 MAE 는 척도 문제 (scale problem) 을 가지고 있습니다. 반면에 MPE (Mean Percentage Error)와 MAPE (Mean Absolute Percentage Error) 는 0~100%로 표준화를 해주어서 척도 문제가 없다는 특징, 장점이 있습니다. 역시  MAE 와 MAPE 값을 함께 확인해서 해석하는 것이 필요합니다. 0에 근접할 수록 시계열 예측모델이 잘 적합되었다고 평가할 수 있으며, MAE의 부호(+, -)로 과대 혹은 과소예측의 방향을 파악할 수 있습니다. 

 

 

 

7. 평균절대비율오차 (MAPE, mean absolute percentage error)

 

mean absolute percentage error

 

 

 

8. 수정결정계수 (Adj. R-square) 

 

adjusted R2

 

결정계수 R2 는 SST 에서 예측 모델이 설명하는 부분의 비율(R2 = SSR/SST=1-SSE/SST)을 의미합니다. 그런데 결정계수 R2는 모수의 개수가 증가할 수록 이에 비례하여 증가하는 경향이 있습니다. 이러한 문제점을 바로잡기 위해 예측에 기여하지 못하는 모수가 포함될 경우 패널티를 부여해서 결정계수의 값을 낮추어주게 수정한 것이 바로 수정결정계수(Adjusted R2) 입니다. (위의 식에서 k 는 모델에 포함된 모수의 개수를 의미합니다.)

 

위의 2~7번의 통계량들은 SSE(Error Sum of Square)를 기반으로 하다보니 값이 작을 수록 좋은 모델을 의미하였다면, 8번의 수정결정계수(Adj. R2)는 예측모델이 설명력과 관련된 지표로서 값이 클 수록 더 좋은 모델을 의미합니다.(1에 가까울 수록 우수) 

 

 

 

9. AIC, SBC, APC

 

AIC (Akaike's Information Criterion), 

SBC (Schwarz's Bayesian Criterion), 

APC (Amemiya's Prediction Criterion)

AIC, SBC, APC

 

위의 AIC, SBC, APC 지표들도 SSE(Error Sum of Square) 에 기반한 지표들로서, 값이 작을 수록 더 잘 적합된 모델을 의미합니다. 이들 지표 역시 SSE 를 기본으로 해서 여기에 모델에서 사용한 모수의 개수(k, 패널티로 사용됨), 관측치의 개수(T, 관측치가 많을 수록 리워드로 사용됨)를 추가하여 조금씩 변형을 한 통계량들입니다. 

 

 

다음번 포스팅에서는 이들 지표를 사용해서 Python으로 하나의 시계열 자료에 대해 여러 개의 지수 평활법 기법들을 적용해서 가장 모형 적합도가 높은 모델을 찾아보겠습니다.(https://rfriend.tistory.com/671)

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. 

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 샘플 크기가 다른 2개 이상의 집단에 대해 평균의 차이가 존재하는지를 검정하는 일원분산분석(one-way ANOVA)에 대해 scipy 모듈의 scipy.stats.f_oneway() 메소드를 사용해서 분석하는 방법(rfriend.tistory.com/638)을 소개하였습니다. 

 

이번 포스팅에서는 2개 이상의 집단에 대해 pandas DataFrame에 들어있는 여러 개의 숫자형 변수(one-way ANOVA for multiple numeric variables in pandas DataFrame) 별로 일원분산분석 검정(one-way ANOVA test)을 하는 방법을 소개하겠습니다. 

 

숫자형 변수와 집단 변수의 모든 가능한 조합을 MultiIndex 로 만들어서 statsmodels.api 모듈의 stats.anova_lm() 메소드의 모델에 for loop 순환문으로 변수를 바꾸어 가면서 ANOVA 검정을 하도록 작성하였습니다. 

 

 

 

먼저, 3개의 집단('grp 1', 'grp 2', 'grp 3')을 별로 'x1', 'x2', 'x3, 'x4' 의 4개의 숫자형 변수를 각각 30개씩 가지는 가상의 pandas DataFrame을 만들어보겠습니다. 이때 숫자형 변수는 모두 정규분포로 부터 난수를 발생시켜 생성하였으며, 'x3'와 'x4'에 대해서는 집단3 ('grp 3') 의 평균이 다른 2개 집단의 평균과는 다른 정규분포로 부터 난수를 발생시켜 생성하였습니다.  

 

아래의 가상 데이터셋은 결측값이 없이 만들었습니다만, 실제 기업에서 쓰는 데이터셋에는 혹시 결측값이 존재할 수도 있으므로 결측값을 없애거나 또는 결측값을 그룹 별 평균으로 대체한 후에 one-way ANOVA 를 실행하기 바랍니다. 

 

## Creating sample dataset
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# generate 90 IDs
id = np.arange(90) + 1

# Create 3 groups with 30 observations in each group.
from itertools import chain, repeat
grp = list(chain.from_iterable((repeat(number, 30) for number in [1, 2, 3])))

# generate random numbers per each groups from normal distribution
np.random.seed(1004)

# for 'x1' from group 1, 2 and 3
x1_g1 = np.random.normal(0, 1, 30)
x1_g2 = np.random.normal(0, 1, 30)
x1_g3 = np.random.normal(0, 1, 30)

# for 'x2' from group 1, 2 and 3
x2_g1 = np.random.normal(10, 1, 30)
x2_g2 = np.random.normal(10, 1, 30)
x2_g3 = np.random.normal(10, 1, 30)

# for 'x3' from group 1, 2 and 3
x3_g1 = np.random.normal(30, 1, 30)
x3_g2 = np.random.normal(30, 1, 30)
x3_g3 = np.random.normal(50, 1, 30) # different mean

x4_g1 = np.random.normal(50, 1, 30)
x4_g2 = np.random.normal(50, 1, 30)
x4_g3 = np.random.normal(20, 1, 30) # different mean

# make a DataFrame with all together
df = pd.DataFrame({'id': id, 
                   'grp': grp, 
                   'x1': np.concatenate([x1_g1, x1_g2, x1_g3]), 
                   'x2': np.concatenate([x2_g1, x2_g2, x2_g3]), 
                   'x3': np.concatenate([x3_g1, x3_g2, x3_g3]), 
                   'x4': np.concatenate([x4_g1, x4_g2, x4_g3])})
                   
df.head()
[Out] 

id	grp	x1	x2	x3	x4
0	1	1	0.594403	10.910982	29.431739	49.232193
1	2	1	0.402609	9.145831	28.548873	50.434544
2	3	1	-0.805162	9.714561	30.505179	49.459769
3	4	1	0.115126	8.885289	29.218484	50.040593
4	5	1	-0.753065	10.230208	30.072990	49.601211


df[df['grp'] == 3].head()
[Out] 

id	grp	x1	x2	x3	x4
60	61	3	-1.034244	11.751622	49.501195	20.363374
61	62	3	0.159294	10.043206	50.820755	19.800253
62	63	3	0.330536	9.967849	50.461775	20.993187
63	64	3	0.025636	9.430043	50.209187	17.892591
64	65	3	-0.092139	12.543271	51.795920	18.883919

 

 

 

가령, 'x3' 변수에 대해 집단별로 상자 그래프 (Box plot for 'x3' by groups) 를 그려보면, 아래와 같이 집단1과 집단2는 유사한 반면에 집단3은 평균이 차이가 많이 나게 가상의 샘플 데이터가 생성되었음을 알 수 있습니다. 

 

## Boxplot for 'x3' by 'grp'
plt.rcParams['figure.figsize'] = [10, 6]
sns.boxplot(x='grp', y='x3', data=df)
plt.show()

 

 

여러개의 변수에 대해 일원분산분석을 하기 전에, 먼저 이해를 돕기 위해 Python의 statsmodels.api 모듈의 stats.anova_lm() 메소드를 사용해서 'x1' 변수에 대해 집단(집단 1/2/3)별로 평균이 같은지 일원분산분석으로 검정을 해보겠습니다. 

 

    - 귀무가설(H0) : 집단1의 x1 평균 = 집단2의 x1 평균 = 집단3의 x1 평균

    - 대립가설(H1) : 적어도 1개 이상의 집단의 x1 평균이 다른 집단의 평균과 다르다. (Not H0)

 

# ANOVA for x1 and grp
import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols('x1 ~ grp', data=df).fit()
sm.stats.anova_lm(model, typ=1)
[Out]

df	sum_sq	mean_sq	F	PR(>F)
grp	1.0	0.235732	0.235732	0.221365	0.639166
Residual	88.0	93.711314	1.064901	NaN	NaN

 

일원분산분석 결과 F 통계량이 0.221365, p-value가 0.639 로서 유의수준 5% 하에서 귀무가설을 채택합니다. 즉, 3개 집단 간 x1의 평균의 차이는 없다고 판단할 수 있습니다. (정규분포 X ~ N(0, 1) 를 따르는 모집단으로 부터 무작위로 3개 집단의 샘플을 추출했으므로 차이가 없게 나오는게 맞겠습니다.)

 

 

한개의 변수에 대한 일원분산분석하는 방법을 알아보았으니, 다음으로는 3개 집단별로 여러개의 연속형 변수인 'x1', 'x2', 'x3', 'x4' 에 대해서 for loop 순환문으로 돌아가면서 일원분산분석을 하고, 그 결과를 하나의 DataFrame에 모아보도록 하겠습니다. 

 

(1) 먼저, 일원분산분석을 하려는 모든 숫자형 변수와 집단 변수에 대한 가능한 조합의 MultiIndex 를 생성해줍니다. 

 

# make a multiindex for possible combinations of Xs and Group
num_col = ['x1','x2', 'x3', 'x4']
cat_col =  ['grp']
mult_idx = pd.MultiIndex.from_product([num_col, cat_col],
                                   names=['x', 'grp'])

print(mult_idx)
[Out]
MultiIndex([('x1', 'grp'),
            ('x2', 'grp'),
            ('x3', 'grp'),
            ('x4', 'grp')],
           names=['x', 'grp'])
           

 

 

(2) for loop 순환문(for x, grp in mult_idx:)으로 model = ols('{} ~ {}'.format(x, grp) 의 선형모델의  y, x 부분의 변수 이름을 바꾸어가면서 sm.stats.anova_lm(model, typ=1) 로 일원분산분석을 수행합니다. 이렇게 해서 나온 일원분산분석 결과 테이블을 anova_tables.append(anova_table) 로 순차적으로 append 해나가면 됩니다.  

 

# ANOVA test for multiple combinations of X and Group
import statsmodels.api as sm
from statsmodels.formula.api import ols

anova_tables = []
for x, grp in mult_idx:
    model = ols('{} ~ {}'.format(x, grp), data=df).fit()
    anova_table = sm.stats.anova_lm(model, typ=1)
    anova_tables.append(anova_table)

df_anova_tables = pd.concat(anova_tables, keys=mult_idx, axis=0)

df_anova_tables
[Out]

df	sum_sq	mean_sq	F	PR(>F)
x1	grp	grp	1.0	0.235732	0.235732	0.221365	6.391661e-01
Residual	88.0	93.711314	1.064901	NaN	NaN
x2	grp	grp	1.0	0.448662	0.448662	0.415853	5.206912e-01
Residual	88.0	94.942885	1.078896	NaN	NaN
x3	grp	grp	1.0	6375.876120	6375.876120	259.202952	5.779374e-28
Residual	88.0	2164.624651	24.598007	NaN	NaN
x4	grp	grp	1.0	13760.538009	13760.538009	256.515180	8.145953e-28
Residual	88.0	4720.684932	53.644147	NaN	NaN

anova tables

 

 

만약 특정 변수에 대한 일원분산분석 결과만을 조회하고 싶다면, 아래처럼 DataFrame의 MultiIndex 에 대해 인덱싱을 해오면 됩니다. 가령, 'x3' 에 대한 집단별 평균 차이 여부를 검정한 결과는 아래처럼 인덱싱해오면 됩니다. 

 

## Getting values of 'x3' from ANOVA tables
df_anova_tables.loc[('x3', 'grp', 'grp')]
[Out]

df         1.000000e+00
sum_sq     6.375876e+03
mean_sq    6.375876e+03
F          2.592030e+02
PR(>F)     5.779374e-28
Name: (x3, grp, grp), dtype: float64

 

 

F 통계량과 p-value 에 대해서 조회하고 싶으면 위의 결과에서 DataFrame 의 칼럼 이름으로 선택해오면 됩니다. 

 

# F-statistic
df_anova_tables.loc[('x3', 'grp', 'grp')]['F']
[Out]
259.2029515179077


# P-value
df_anova_tables.loc[('x3', 'grp', 'grp')]['PR(>F)']
[Out]
5.7793742588216585e-28

 

 

 

MultiIndex 를 인덱싱해오는게 좀 불편할 수 도 있는데요, 이럴 경우  df_anova_tables.reset_index() 로  MultiIndex 를 칼럼으로 변환해서 사용할 수도 있습니다. 

# resetting index to columns
df_anova_tables_2 = df_anova_tables.reset_index().dropna()


df_anova_tables_2
[Out]

level_0	level_1	level_2	df	sum_sq	mean_sq	F	PR(>F)
0	x1	grp	grp	1.0	0.235732	0.235732	0.221365	6.391661e-01
2	x2	grp	grp	1.0	0.448662	0.448662	0.415853	5.206912e-01
4	x3	grp	grp	1.0	6375.876120	6375.876120	259.202952	5.779374e-28
6	x4	grp	grp	1.0	13760.538009	13760.538009	256.515180	8.145953e-28

 

 

Greenplum DB에서 PL/Python (또는 PL/R)을 사용하여 여러개의 숫자형 변수에 대해 일원분산분석을 분산병렬처리하는 방법(one-way ANOVA in parallel using PL/Python on Greenplum DB)은 rfriend.tistory.com/640 를 참고하세요. 

 

 

[reference] 

* ANOVA test using Python statsmodels
 
: https://www.statsmodels.org/stable/generated/statsmodels.stats.anova.anova_lm.html

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 2개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 2개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  이론적인 부분에 대한 좀더 자세한 내용은 https://rfriend.tistory.com/131 를 참고하세요. 

 

one-way ANOVA  일원분산분석

 

이번 포스팅에서는 Python의  scipy 모듈의 stats.f_oneway() 메소드를 사용하여 샘플의 크기가 서로 다른 3개 그룹 간 평균에 차이가 존재하는지 여부를 일원분산분석(one-way ANOVA)으로 분석하는 방법을 소개하겠습니다. 

 

분산분석(Analysis Of Variance) 검정의 3가지 가정사항을 고려해서, 샘플 크기가 서로 다른 가상의 3개 그룹의 예제 데이터셋을 만들어보겠습니다. 

 

[ 분산분석  검정의 가정사항 (assumptions of ANOVA test) ]

  (1) 독립성: 각 샘플 데이터는 서로 독립이다. 
  (2) 정규성: 각 샘플 데이터는 정규분포를 따르는 모집단으로 부터 추출되었다. 
  (3) 등분산성: 그룹들의 모집단의 분산은 모두 동일하다. 

 

먼저, 아래의 예제 샘플 데이터셋은 그룹1과 그룹2의 평균은 '0'으로 같고, 그룹3의 평균은 '5'로서 나머지 두 그룹과 다르게 난수를 발생시켜 가상으로 만든 것입니다. 

 

# 3 groups of dataset with different sized samples 
import numpy as np
import pandas as pd
np.random.seed(1004)

data1 = np.random.normal(0, 1, 50)
data2 = np.random.normal(0, 1, 40)
data3 = np.random.normal(5, 1, 30) # different mean

data123 = [data1, data2, data3]


print(data123)
[Out]
[array([ 0.59440307,  0.40260871, -0.80516223,  0.1151257 , -0.75306522,
       -0.7841178 ,  1.46157577,  1.57607553, -0.17131776, -0.91448182,
        0.86013945,  0.35880192,  1.72965706, -0.49764822,  1.7618699 ,
        0.16901308, -1.08523701, -0.01065175,  1.11579838, -1.26497153,
       -1.02072516, -0.71342119,  0.57412224, -0.45455422, -1.15656742,
        1.29721355, -1.3833716 ,  0.3205909 , -0.59086187, -1.43420648,
        0.60998011,  0.51266756,  1.9965168 ,  1.42945668,  1.82880165,
       -1.40997132,  0.49433367,  0.9482873 , -0.35274099, -0.15359935,
       -1.18356064, -0.75440273, -0.85981073,  1.14256322, -2.21331694,
        0.90651805,  2.23629   ,  1.00743665,  1.30584548,  0.46669171]), array([-0.49206651, -0.08727244, -0.34919043, -1.11363541, -1.71982966,
       -0.14033817,  0.90928317, -0.60012686,  1.03906073, -0.03332287,
       -1.03424396,  0.15929405,  0.33053582,  0.02563551, -0.09213904,
       -0.91851177,  0.3099129 , -1.24211638, -0.33113027, -1.64086666,
       -0.27539834, -0.05489003,  1.50604364, -1.37756156, -1.25561652,
        0.16120867, -0.42121705,  0.2341905 , -1.20155195,  1.48131392,
        0.29105321,  0.4022031 , -0.41466037,  1.00502917,  1.45376705,
       -0.07038153,  0.52897801, -2.37895295, -0.75054747,  1.10641762]), array([5.91098191, 4.14583073, 4.71456131, 3.88528941, 5.23020779,
       5.12814125, 3.44610618, 5.36308351, 4.69463298, 7.49521024,
       5.41246681, 3.8724271 , 4.60265531, 4.60082925, 4.9074518 ,
       3.8141367 , 6.4457503 , 4.27553455, 3.63173152, 6.25540542,
       3.77536981, 7.19435668, 6.25339789, 6.43469547, 4.27431061,
       5.16694916, 7.21065725, 5.68274021, 4.81732021, 3.81650656])]

 

 

 

위의 3개 그룹의 커널밀도추정 그래프 (Kernel Density Estimates Plot)를 겹쳐서 그려보면, 아래와 같이 2개 집단은 서로 평균이 비슷하고 나머지 1개 집단은 평균이 확연히 다르다는 것을 직관적으로 알 수 있습니다. 

# Kernel Density Estimate Plot
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['figure.figsize'] = [10, 6]

sns.kdeplot(data1)
sns.kdeplot(data2)
sns.kdeplot(data3)
plt.show()

 

 

상자 그래프 (Box Plot) 으로 3개 집단 간 평균의 위치와 퍼짐 정도를 비교해보면, 역시 아래와 같이 그룹1과 그룹2는 서로 중심위치가 서로 비슷하고 그룹3만 중심위치가 확연히 다름을 알 수 있습니다. 

 

# Boxplot
sns.boxplot(data=data123)
plt.xlabel("Group", fontsize=14)
plt.ylabel("Value", fontsize=14)
plt.show()

 

 

이제 scipy 모듈의 scipy.stats.f_oneway() 메소드를 사용해서 서로 다른 개수의 샘플을 가진 3개 집단에 대해 평균의 차이가 존재하는지 여부를 일원분산분석을 통해 검정을 해보겠습니다. 

 

  - 귀무가설(H0): 3 집단의 평균은 모두 같다. (mu1 = mu2 = m3)

  - 대립가설(H1):  3 집단의 평균은 같지 않다.(적어도 1개 집단의 평균은 같지 않다) (Not H0)

 

F통계량은 매우 큰 값이고 p-value가 매우 작은 값이 나왔으므로 유의수준 5% 하에서 귀무가설을 기각하고 대립가설을 채택합니다. 즉, 3개 집단 간 평균의 차이가 존재한다고 평가할 수 있습니다

# ANOVA with different sized samples using scipy
from scipy import stats

stats.f_oneway(data1, data2, data3)
[Out] 
F_onewayResult(statistic=262.7129127080777, pvalue=5.385523527223916e-44)

 

 

F통계량과 p-value 를 각각 따로 따로 반환받을 수도 있습니다. 

# returning f-statistic and p-value
f_val, p_val = stats.f_oneway(*data123)

print('f-statistic:', f_val)
print('p-vale:', p_val)
[Out]
f-statistic: 262.7129127080777
p-vale: 5.385523527223916e-44

 

 

scipy 모듈의  stats.f_oneway() 메소드를 사용할 때 만약 데이터에 결측값(NAN)이 포함되어 있으면  'NAN'을 반환합니다. 위에서 만들었던 'data1'  numpy array 의 첫번째 값을  np.nan 으로 대체한 후에 scipy.stats.f_oneway() 로 일원분산분석 검정을 해보면  'NAN'(Not A Number)을 반환한 것을 볼 수 있습니다. 

 

# if there is 'NAN', then returns 'NAN'
data1[0] = np.nan
print(data1)
[Out]
array([        nan,  0.40260871, -0.80516223,  0.1151257 , -0.75306522,
       -0.7841178 ,  1.46157577,  1.57607553, -0.17131776, -0.91448182,
        0.86013945,  0.35880192,  1.72965706, -0.49764822,  1.7618699 ,
        0.16901308, -1.08523701, -0.01065175,  1.11579838, -1.26497153,
       -1.02072516, -0.71342119,  0.57412224, -0.45455422, -1.15656742,
        1.29721355, -1.3833716 ,  0.3205909 , -0.59086187, -1.43420648,
        0.60998011,  0.51266756,  1.9965168 ,  1.42945668,  1.82880165,
       -1.40997132,  0.49433367,  0.9482873 , -0.35274099, -0.15359935,
       -1.18356064, -0.75440273, -0.85981073,  1.14256322, -2.21331694,
        0.90651805,  2.23629   ,  1.00743665,  1.30584548,  0.46669171])

# returns 'nan'
stats.f_oneway(data1, data2, data3)
[Out] 
F_onewayResult(statistic=nan, pvalue=nan)

 

 

샘플 데이터에 결측값이 포함되어 있는 경우, 결측값을 먼저 제거해주고 일원분산분석 검정을 실시해주시기 바랍니다. 

 

# get rid of the missing values before applying ANOVA test
stats.f_oneway(data1[~np.isnan(data1)], data2, data3)
[Out]
F_onewayResult(statistic=260.766426640122, pvalue=1.1951277551195217e-43)

 

 

위의  일원분산분석(one-way ANOVA) 는 2개 이상의 그룹 간 평균의 차이가 존재하는지만을 검정할 뿐이며, 집단 간 평균의 차이가 존재한다는 대립가설을 채택하게 된 경우 어느 그룹 간 차이가 존재하는지는 사후검정 다중비교(post-hoc pair-wise multiple comparisons)를 통해서 알 수 있습니다. (rfriend.tistory.com/133)

 

pandas DataFrame 데이터셋에서 여러개의 숫자형 변수에 대해 for loop 순환문을 사용하여 집단 간 평균 차이 여부를 검정하는 방법은 rfriend.tistory.com/639 포스팅을 참고하세요. 

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. 

 

[reference]

* scipy.stats.f_oneway : https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html

 

scipy.stats.f_oneway — SciPy v1.6.3 Reference Guide

G.W. Heiman, “Understanding research methods and statistics: An integrated introduction for psychology”, Houghton, Mifflin and Company, 2001.

docs.scipy.org

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 선형회귀모형에 대한 각 관측치별 변수별 기여도 분석 (each variable contribution per each observations)에 대해서 소개하겠습니다. 


변수 중요도 (variable importance, feature importance)가 전체 관측치를 사용해 적합한 모델 단위의 변수별 (상대적) 중요도를 나타내는 것이라면, 이번 포스팅에서 소개하려는 관측치별 변수별 기여도(민감도)는 개별 관측치 단위에서 한개의 칼럼이 모델 예측치에 얼마나 기여를 하는지를 분석해보는 것입니다. 


선형회귀 모형을 예로 들면, 변수 중요도는 회귀모형의 회귀 계수라고 할 수 있겠구요, 관측치별 변수별 기여도는 특정 관측치의 특정 칼럼값만 그대로 사용하고 나머지 변수값에는 '0'을 대입해서 나온 예측치라고 할 수 있겠습니다. 변수 중요도가 Global 하게 적용되는 model weights 라고 한다면, 관측치별 변수별 기여도는 specific variable's weight * value 라고 할 수 있겠습니다. 




변수 중요도를 분석했으면 됐지, 왜 관측치별 변수별 기여도 분석을 할까 궁금할 수도 있을 텐데요, 관측치별 변수별 기여도 분석은 특정 관측치별로 개별적으로 어떤 변수가 예측치에 크게 영향을 미쳤는지 파악하기 위해서 사용합니다. (동어반복인가요? ^^;) 


모델의 변수별 가중치(high, low)와 각 관측치별 변수별 값(high, low)의 조합별로 아래와 같이 관측치별 변수별 예측치에 대한 기여도가 달라지게 됩니다. 가령, 관측치별 변수별 예측치에 대한 기여도가 높으려면 모델의 가중치(중요도)도 높아야 하고 동시에 개별 관측치의 해당 변수의 관측값도 역시 높아야 합니다. 



관측치별 변수별 기여도 분석은 복잡하고 어려운 이론에 기반을 둔 분석이라기 보다는 단순한 산수를 반복 연산 프로그래밍으로 계산하는 것이므로 아래의 예를 따라가다 보면 금방 이해할 수 있을 것이라고 생각합니다. 



예를 들어보기 위해 공개된 전복(abalone) 데이터를 사용하여 전체 무게(whole weight)를 예측하는 선형회귀모형을 적합하고, 관측치별로 각 변수별 예측치에 기여한 정도를 분석해보겠습니다. 



  1. 선형회귀모형에 대한 관측치별 변수별 기여도(민감도) 분석 

    (Sensitivity analysis for linear regression model)



(1) abalone dataset 을 pandas DataFrame으로 불러오기



import numpy as np

import pandas as pd


url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data'

abalone = pd.read_csv(url

                      , sep=','

                      , header=None

                      , names=['sex', 'length', 'diameter', 'height', 

                               'whole_weight', 'shucked_weight', 'viscera_weight',

                               'shell_weight', 'rings']

                     , index_col=None)


abalone.head()

[Out]:

sexlengthdiameterheightwhole_weightshucked_weightviscera_weightshell_weightrings
0M0.4550.3650.0950.51400.22450.10100.15015
1M0.3500.2650.0900.22550.09950.04850.0707
2F0.5300.4200.1350.67700.25650.14150.2109
3M0.4400.3650.1250.51600.21550.11400.15510
4I0.3300.2550.0800.20500.08950.03950.0557




(2) X, y 변수 데이터셋 생성 (creating X and y dataset)


예측하려는 목표변수 y 로는 '전체 무게(whole weight)' 변수를 사용하겠으며, 설명변수 X 에는 'sex', 'length', 'diameter', 'height', 'shell_weight', 'rings' 의 6개 변수를 사용하겠습니다. 이중에서 'sex' 변수는 범주형 변수이므로 가변수(dummy variable)로 변환하였습니다. 



# transformation of categorical variable to dummy variable

abalone['sex'].unique()

[Out]: array(['M', 'F', 'I'], dtype=object)


abalone['sex_m'] = np.where(abalone['sex'] == 'M', 1, 0)

abalone['sex_f'] = np.where(abalone['sex'] == 'F', 1, 0)


# get X variables

X = abalone[["sex_m", "sex_f", "length", "diameter", "height", "shell_weight", "rings"]]


import statsmodels.api as sm

X = sm.add_constant(X) # add a constant to model

print(X)


[Out]:
Index(['const', 'sex_m', 'sex_f', 'length', 'diameter', 'height',
       'shell_weight', 'rings'],
      dtype='object')

      const  sex_m  sex_f  length  diameter  height  shell_weight  rings
0       1.0      1      0   0.455     0.365   0.095        0.1500     15
1       1.0      1      0   0.350     0.265   0.090        0.0700      7
2       1.0      0      1   0.530     0.420   0.135        0.2100      9
3       1.0      1      0   0.440     0.365   0.125        0.1550     10
4       1.0      0      0   0.330     0.255   0.080        0.0550      7
...     ...    ...    ...     ...       ...     ...           ...    ...
4172    1.0      0      1   0.565     0.450   0.165        0.2490     11
4173    1.0      1      0   0.590     0.440   0.135        0.2605     10
4174    1.0      1      0   0.600     0.475   0.205        0.3080      9
4175    1.0      0      1   0.625     0.485   0.150        0.2960     10
4176    1.0      1      0   0.710     0.555   0.195        0.4950     12

[4177 rows x 8 columns]


# get y value

y = abalone["whole_weight"]

 



Train:Test = 8:2 의 비율로 데이터셋을 분할 (Split train and test set)하겠습니다. 



# train, test set split

from sklearn.model_selection import train_test_split


X_train, X_test, y_train, y_test = train_test_split(X, 

                                                    y, 

                                                    test_size=0.2, 

                                                    random_state=2004)

print('Shape of X_train:', X_train.shape

print('Shape of X_test:', X_test.shape

print('Shape of y_train:', y_train.shape

print('Shape of y_test:', y_test.shape)

[Out]:

Shape of X_train: (3341, 8) Shape of X_test: (836, 8) Shape of y_train: (3341,) Shape of y_test: (836,)





(3) 선형회귀모형 적합 (fitting linear regression model)


이번 포스팅이 선형회귀모형에 대해서 설명하는 것이 목적이 아니기 때문에 회귀모형을 적합하는데 필요한 가정사항 진단 및 조치생략하고, 그냥 statsmodels.api 라이브러리를 이용해서 회귀모델을 적합해서 바로 민감도 분석(Sensitivity analysis)으로 넘어가겠습니다


이번 포스팅의 주제인 '관측치별 변수별 기여도 분석'에서 사용하는 모델은 어떤 알고리즘도 전부 가능하므로 개별 관측치별 변수별 영향력을 해석을 하는데 유용하게 사용할 수 있습니다. (가령, 블랙박스 모형인 Random Forest, Deep Learning 등도 가능)



# multivariate linear regression model

import statsmodels.api as sm


lin_reg = sm.OLS(y_train, X_train).fit()


lin_reg.summary()

OLS Regression Results
Dep. Variable:whole_weightR-squared:0.942
Model:OLSAdj. R-squared:0.942
Method:Least SquaresF-statistic:7738.
Date:Fri, 17 Jan 2020Prob (F-statistic):0.00
Time:20:06:12Log-Likelihood:2375.3
No. Observations:3341AIC:-4735.
Df Residuals:3333BIC:-4686.
Df Model:7
Covariance Type:nonrobust
coefstd errtP>|t|[0.0250.975]
const-0.36090.015-23.7950.000-0.391-0.331
sex_m0.03490.0066.0390.0000.0240.046
sex_f0.02170.0063.4930.0000.0100.034
length1.25390.10711.7240.0001.0441.464
diameter0.07490.1350.5540.580-0.1900.340
height0.38800.0884.4100.0000.2160.561
shell_weight2.43840.03864.8020.0002.3652.512
rings-0.01540.001-18.4230.000-0.017-0.014
Omnibus:698.165Durbin-Watson:1.978
Prob(Omnibus):0.000Jarque-Bera (JB):4041.763
Skew:0.866Prob(JB):0.00
Kurtosis:8.102Cond. No.866.


# prediction

predicted = lin_reg.predict(X_test)

actual = y_test


act_pred_df = pd.DataFrame({'actual': actual

                            , 'predicted': predicted

                            , 'error': actual - predicted})


act_pred_df.head()

[Out]:

actualpredictederror
34730.0455-0.0947030.140203
35230.09700.0206040.076396
18620.51850.690349-0.171849
29661.48201.3589500.123050
6590.95851.137433-0.178933



# Scatter Plot: Actual vs. Predicted

import matplotlib.pyplot as plt


plt.plot(act_pred_df['actual'], act_pred_df['predicted'], 'o')

plt.xlabel('actual', fontsize=14)

plt.ylabel('predicted', fontsize=14)

plt.show()


# RMSE (Root Mean Squared Error)

from sklearn.metrics import mean_squared_error

from math import sqrt


rmse = sqrt(mean_squared_error(actual, predicted))

rmse

[Out]: 0.11099248621173345





(4) 관측치별 변수별 예측모델 결과에 대한 기여도 분석 (contribution per each variables from specific observation)


아래 예에서는 전체 836개의 test set 관측치 중에서 '1번 관측치 (1st observation)' 의 8개 변수별 기여도를 분석해보겠습니다. 



X_test.shape

[Out]: (836, 8)


# get 1st observation's value as an example

X_i = X_test.iloc[0, :]

X_i

[Out]:

const           1.000
sex_m           0.000
sex_f           0.000
length          0.210
diameter        0.150
height          0.055
shell_weight    0.013
rings           4.000
Name: 3473, dtype: float64




관측치별 변수별 기여도(민감도, variable's contribution & sensitivity) 분석의 핵심 개념은 전체 데이터를 사용해서 적합된 모델에 특정 관측치의 변수별 값을 변수별로 순서대로 돌아가면서 해당 변수 측정값은 그대로 사용하고 나머지 변수들의 값은 '0'으로 대체해서 예측을 해보는 것입니다. 아래 코드의 for i, j in enumerate(X_i): X_mat[i, i] = j 부분을 유심히 살펴보시기 바랍니다. 



# all zeros matrix with shape of [col_num, col_num]

X_mat = np.zeros(shape=[X_i.shape[0], X_i.shape[0]])

X_mat

[Out]:
array([[0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0.]])


# fill only 1 variable's value and leave '0' for the others

for i, j in enumerate(X_i):

    X_mat[i, i] = j


X_mat

[Out]:
array([[1.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.21 , 0.   , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.15 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.055, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.013, 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 4.   ]])




바로 위에서 만든 X_mat 행렬 (1번 관측치의 각 변수별로 하나씩 실제 값을 사용하고, 나머지는 '0'으로 대체)을 사용해서 (3)에서 적합한 선형회귀모형으로 y값을 예측(lin_reg.predict(X_mat))을 하면 우리가 원하는 '1번 관측치의 개별 변수별 y값 예측에 대한 기여도'를 구할 수 있습니다. 


개별 변수별 y 예측값에 대한 기여도를 기준으로 내림차순 정렬해서 DataFrame을 만들고, 가로로 누운 막대그래프도 그려보았습니다. 


# sensitivity analysis

sensitivity_df = pd.DataFrame({'x': X_test.iloc[0, :]

                               , 'contribution_x': lin_reg.predict(X_mat)}).\

                 sort_values(by='contribution_x', ascending=False)


sensitivity_df

[Out]:

xcontribution_x
length0.2100.263315
shell_weight0.0130.031699
height0.0550.021341
diameter0.1500.011237
sex_m0.0000.000000
sex_f0.0000.000000
rings4.000-0.061437
const1.000-0.360857



horizontal bar plot by column's contribution

sensitivity_df['contribution_x'].plot(kind='barh', figsize=(10, 5))

plt.title('Sensitivity Analysis', fontsize=16)

plt.xlabel('Contribution', fontsize=16)

plt.ylabel('Variable', fontsize=16)

plt.yticks(fontsize=14)

plt.show()


물론, 당연하게 관측치별 변수별 기여도 (민감도) 분석에서 나온 결과를 전부 다 합치면 애초의 모형에 해당 관측치의 전체 변수별 값을 넣어서 예측한 값과 동일한 결과가 나옵니다. 



# result from sum of contribution analysis

sum(sensitivity_df['contribution_x'])

[Out]: -0.09470251191012563


# result from linear regression model's prediction

lin_reg.predict(X_test.iloc[0, :].to_numpy())

[Out]: array([-0.09470251])

 




(5) 관측치별 변수별 예측 기여도 (민감도) 분석을 위한 사용자 정의 함수


관측치별 변수별 기여도 분석 결과를 --> pandas DataFrame으로 저장하고, --> 기여도별로 정렬된 값을 기준으로 옆으로 누운 막대그래프를 그려주는 사용자 정의함수를 만들어보겠습니다. 



# UDF for contribution(sensitivity) analysis per each variables

def sensitivity_analysis(model, X, idx, bar_plot_yn):

    

    import numpy as np

    import pandas as pd

    import matplotlib.pyplot as plt

    import statsmodels.api as sm

    pd.options.mode.chained_assignment = None

    

    # get one object's X values

    X_i = X.iloc[idx, :]

    

    # make a matrix with zeros with shape of [num_cols, num_cols]

    X_mat = np.zeros(shape=[X_i.shape[0], X_i.shape[0]])

    

    # fil X_mat with values from one by one columns, leaving the ohters zeros

    for i, j in enumerate(X_i):

        X_mat[i, i] = j

    

    # data frame with contribution of each X columns in descending order

    sensitivity_df = pd.DataFrame({'idx': idx

                                   , 'x': X_i

                                   , 'contribution_x': model.predict(X_mat)}).\

                    sort_values(by='contribution_x', ascending=True)

    

    # if bar_plot_yn == True then display it

    col_n = X_i.shape[0]

    if bar_plot_yn == True:

        sensitivity_df['contribution_x'].plot(kind='barh', figsize=(10, 0.7*col_n))

        plt.title('Sensitivity Analysis', fontsize=18)

        plt.xlabel('Contribution', fontsize=16)

        plt.ylabel('Variable', fontsize=16)

        plt.yticks(fontsize=14)

        plt.show()

    

    return sensitivity_df



# check UDF

sensitivity_df = sensitivity_analysis(model=lin_reg, X=X_test, idx=0, bar_plot_yn=True)

sensitivity_df




아래는 위에서 정의한 sensitivity_analysis() 사용자 정의 함수에서 bar_plot_yn=False 로 설정을 해서 옆으로 누운 막대그래프는 그리지 말고 기여도 분석 결과 DataFrame만 반환하라고 한 경우입니다. 



# without bar plot (bar_plot_yn=False)

sensitivity_df = sensitivity_analysis(model=lin_reg, X=X_test, idx=0, bar_plot_yn=False)

sensitivity_df

[Out]:

idxxcontribution_x
length00.2100.263315
shell_weight00.0130.031699
height00.0550.021341
diameter00.1500.011237
sex_m00.0000.000000
sex_f00.0000.000000
rings04.000-0.061437
const01.000-0.360857





(6) 다수의 관측치에 대해 '개별 관측치별 변수별 예측 기여도 분석'


위에서 부터 10개의 관측치를 선별해서, 개별 관측치별 각 변수별로 위의 (5)번에서 정의한 sensitivity_analysis() 사용자정의함수와 for loop 반복문을 사용해서 변수 기여도를 분석해보겠습니다. (단, 변수 기여도 막대 그래프는 생략)



# calculate sensitivity of each columns of the first 10 objects using for loop


# blank DataFrame to save the sensitivity results together

sensitivity_df_all = pd.DataFrame()

to_idx = 10


for idx in range(0, to_idx):

    sensitivity_df_idx = sensitivity_analysis(model=lin_reg

                                              , X=X_test

                                              , idx=idx

                                              , bar_plot_yn=False)

    

    sensitivity_df_all = pd.concat([sensitivity_df_all, sensitivity_df_idx], axis=0)

    

    print("[STATUS]", idx+1, "/", to_idx, "(", 100*(idx+1)/to_idx, "%) is completed...")


[STATUS] 1 / 10 ( 10.0 %) is completed...
[STATUS] 2 / 10 ( 20.0 %) is completed...
[STATUS] 3 / 10 ( 30.0 %) is completed...
[STATUS] 4 / 10 ( 40.0 %) is completed...
[STATUS] 5 / 10 ( 50.0 %) is completed...
[STATUS] 6 / 10 ( 60.0 %) is completed...
[STATUS] 7 / 10 ( 70.0 %) is completed...
[STATUS] 8 / 10 ( 80.0 %) is completed...
[STATUS] 9 / 10 ( 90.0 %) is completed...
[STATUS] 10 / 10 ( 100.0 %) is completed...

 



결과를 한번 확인해볼까요?



sensitivity_df_all[:20]

[Out]:

idxxcontribution_x
length00.21000.263315
shell_weight00.01300.031699
height00.05500.021341
diameter00.15000.011237
sex_m00.00000.000000
sex_f00.00000.000000
rings04.0000-0.061437
const01.0000-0.360857
length10.26000.326009
shell_weight10.03050.074372
height10.07000.027161
diameter10.20500.015357
sex_m10.00000.000000
sex_f10.00000.000000
rings14.0000-0.061437
const11.0000-0.360857
length20.52000.652018
shell_weight20.18400.448668
height20.11000.042682
diameter20.41000.030713




 2. 로지스틱 회귀모형에 대한 관측치별 변수별 민감도 분석

   (Sensitivity analysis for Logistic Regression model) 


제가 서두에서 민감도분석에 어떤 모델도 가능하도고 말씀드렸습니다. 그러면 이번에는 같은 abalone 데이터셋에 대해 목표변수로 전체무게(whole weight)가 평균보다 크면 '1', 평균보다 작으면 '0' 으로 범주를 구분한 후에, 이를 이진 분류(binary classification)할 수 있는 로지스틱 회귀모형(Logistic Regression)을 적합한 후에 --> 민감도 분석을 적용해보겠습니다. 


먼저, y_cat 이라는 범주형 변수를 만들고, train/ test set으로 분할을 한 후에, 로지스틱 회귀모형(Logistic Regression)을 적합후, 이진분류를 할 수 있는 확률을 계산(확률이 0.5보다 크면 '1'로 분류)해보겠습니다. 



# make a y_category variable: if y is greater or equal to mean, then 1

y_cat = np.where(abalone["whole_weight"] >= np.mean(abalone["whole_weight"]), 1, 0)

y_cat[:20]

[Out]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])


cat_class, counts = np.unique(y_cat, return_counts=True)

dict(zip(cat_class, counts))

[Out]: {0: 2178, 1: 1999}


# train, test set split

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, 

                                                    y_cat, 

                                                    test_size=0.2, 

                                                    random_state=2004)


# fitting logistic regression

import statsmodels.api as sm

pd.options.mode.chained_assignment = None


logitreg = sm.Logit(y_train, X_train)

logitreg_fit = logitreg.fit()


print(logitreg_fit.summary())

[Out]:

Optimization terminated successfully. Current function value: 0.108606 Iterations 11 Logit Regression Results ============================================================================== Dep. Variable: y No. Observations: 3341 Model: Logit Df Residuals: 3333 Method: MLE Df Model: 7 Date: Fri, 17 Jan 2020 Pseudo R-squ.: 0.8431 Time: 21:00:24 Log-Likelihood: -362.85 converged: True LL-Null: -2313.0 Covariance Type: nonrobust LLR p-value: 0.000 ============================================================================= coef std err z P>|z| [0.025 0.975] -------------------------------------------------------------------------------- const -35.3546 2.219 -15.935 0.000 -39.703 -31.006 sex_m 1.3064 0.250 5.217 0.000 0.816 1.797 sex_f 1.1418 0.258 4.420 0.000 0.635 1.648 length 36.6625 5.348 6.855 0.000 26.180 47.145 diameter 11.7960 6.160 1.915 0.055 -0.277 23.869 height 7.6179 2.011 3.789 0.000 3.677 11.559 shell_weight 38.0930 3.286 11.593 0.000 31.653 44.533 rings -0.1062 0.038 -2.777 0.005 -0.181 -0.031 =============================================================================


# prediction

test_prob_logitreg = logitreg_fit.predict(X_test)

test_prob_logitreg.head()

[Out]:
3473    9.341757e-12
3523    2.440305e-10
1862    1.147617e-02
2966    9.999843e-01
659     9.964952e-01
dtype: float64




다음으로 위의 선형회귀모형에 대한 민감도 분석 사용자정의 함수를 재활용하여 로지스틱 회귀모형에도 사용가능하도록 일부 수정해보았습니다. 


아래의 사용자 정의 함수는 로지스틱 회귀모형 적합 시 statsmodels 라이브러리를 사용한 것으로 가정하고 작성하였습니다. 만약 로지스틱 회귀모형의 model 적합에 from sklearn.linear_model import LogisticRegression 의 라이브러리를 사용하였다면 '#'으로 비활성화해놓은 부분을 해제하여 사용하면 됩니다. (predict_proba(X_mat)[:, 1] 의 부분이 달라서 그렇습니다)



# UDF for contribution(sensitivity) analysis per each variables

# task: "LinearReg" or "LogitReg"

def sensitivity_analysis_LinearReg_LogitReg(task, model, X, idx, bar_plot_yn):

    

    import numpy as np

    import pandas as pd

    import matplotlib.pyplot as plt

    import statsmodels.api as sm

    pd.options.mode.chained_assignment = None

    

    # get one object's X values

    X_i = X.iloc[idx, :]

    

    # make a matrix with zeros with shape of [num_cols, num_cols]

    X_mat = np.zeros(shape=[X_i.shape[0], X_i.shape[0]])

    

    # fil X_mat with values from one by one columns, leaving the ohters zeros

    for i, j in enumerate(X_i):

        X_mat[i, i] = j

        

    # data frame with contribution of each X columns in descending order

    sensitivity_df = pd.DataFrame({

        'idx': idx

        , 'task': task

        , 'x': X_i

        , 'contribution_x': model.predict(X_mat)     

    })

    

#     # ==== Remark =====

#     # if you used LogisticRegressionsklearn from sklearn.linear_model

#     # then use codes below

#     if task == "LinearReg":

#         sensitivity_df = pd.DataFrame({

#             'idx': idx

#             , 'task': task

#             , 'x': X_i

#             , 'contribution_x': model.predict(X_mat) 

#         })

        

#     elif task == "LogitReg":

#         sensitivity_df = pd.DataFrame({

#             'idx': idx

#             , 'task': task

#             , 'x': X_i

#             , 'contribution_x': model.predict_proba(X_mat)[:,1] 

#         })

#     else:

#         print('Please choose task one of "LinearReg" or "LogitReg"...')

    

    

    sensitivity_df = sensitivity_df.sort_values(by='contribution_x', ascending=True)

    

    # if bar_plot_yn == True then display it

    col_n = X_i.shape[0]

    if bar_plot_yn == True:

        sensitivity_df['contribution_x'].plot(kind='barh', figsize=(10, 0.7*col_n))

        plt.title('Sensitivity Analysis', fontsize=18)

        plt.xlabel('Contribution', fontsize=16)

        plt.ylabel('Variable', fontsize=16)

        plt.yticks(fontsize=14)

        plt.show()

    

    return sensitivity_df.sort_values(by='contribution_x', ascending=False)



# apply sensitivity analysis function on 1st observation for Logistic Regression

sensitivity_analysis_LinearReg_LogitReg(task="LogitReg"

                                        , model=logitreg_fit

                                        , X=X_test

                                        , idx=0

                                        , bar_plot_yn=True)




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 

 - 시계열의 구성 요인 (time series component factors) (https://rfriend.tistory.com/509)과 

 - 시계열 분해 (time series decomposition) (https://rfriend.tistory.com/510

에 대해서 소개하였습니다. 

 

이번 포스팅에서는 직관적이어서 이해하기 쉽고 (ARIMA 대비 상대적으로요), 또 시계열 자료의 구성요소가 변동이 느리거나 매우 규칙적(periodic)일 때 예측 정확도도 높아서 시계열 자료 예측에 실무적으로 많이 사용하는 지수 평활법 (exponential smoothing) 에 대해서 소개하겠습니다. 

 

지수 평활법의 기본 개념은 "최근 관측한 값에 높은 가중치를 주고, 먼 과거의 관측치에는 낮은 관측치를 주어 가중 평균값을 구한는 것" 입니다. 이때 가중치는 현재로 부터 과거로 갈 수록 지수적으로 감소(exponential decay)시켜서 차등적으로 주는 평활 함수(smoothing function)을 사용하는 것에서 "지수(exponential)" 이름을 따왔습니다. 그리고 여러개의 관측치를 모아서 가중 평균을 사용하기 때문에 "평활(smoothing)"되는 효과가 있어서, 이 둘을 합쳐서 "지수 평활법(exponential smoothing)"이라고 하는 것입니다. 

 

 

[ 지수 분포의 확률밀도함수 곡선 ]

probability density function of exponential distribution

 

이전 Python 으로 차수 m인 단순이동평균(simple moving average) 구하는 방법을 포스팅한적이 있는데요, 단순이동평균은 가중치가 1/m 로 모두 같은 형태인 (m차 관측치만 이동해서 고려하는) 지수 평활법의 특수한 경우로 볼 수도 있겠습니다. 

 

지수 평활법의 기법을 결정하는데는 시도표 (time series plot)를 그려보고 아래의 3가지를 고려합니다. 

 

 

[ 시계열 특성에 따른 지수평활법 기법 결정 ]

 

(1) 시계열 자료에 추세(Trend)가 있는지, 추세가 있다면 1차 선형인지 아니면 비선형 인가? 

    - 추세 없음 (No Trend)                  --> Simple Exponential Smoothing

    - 1차 선형 추세 (Linear Trend)         --> Two Parameter Exponential Smoothing

    - 2차 비선형 추세 (Quadratic Trend) --> Three Parameter Smoothing

 

(2) 시계열 자료에 계절성(Seasonality)이 있는가? 

    - 계절성 없음 (No Seasonality)

    - 계정성 있음 (with Seasonality)  --> Winters' Method

 

(3) 시계열 자료의 계절성이 시간이 지남에 따라 고정(fixed)되어 있는지 아니면 확산(increasing)되는가? 

    - 고정(상수) 계절 변동 (fixed seasonal variation) --> Additive Model

    - 확산 계절 변동 (increasing seasonal variation)  --> Multiplicative Model

 

 

위의 설명을 각 시계열 자료 패턴별로 시도표 예시와 함께 decision tree 형태로 구분해서 지수 평활법 (exponential smoothing) 기법을 짝지어보면 아래와 같습니다. 

 

 

물론 위와 같이 분석가가 눈으로 시도표를 보고 나서 적합한 지수 평활법 기법을 선택하는 것이 아니라, 분석 시스템 성능이 뒷받침이 되고 분석 시간이 오래 소요되어도 문제가 없다면, 혹은 사람 개입 없이 자동화를 하고자 한다면 모든 지수 평활법 기법을 적용해보고 이들 모델 긴 적합도 평가 통계량(가령, RMSE, MAPE, MAE 등)을 비교해 본 후 적합이 가장 잘 된 것을 선택할 수도 있겠습니다.

 

 

[ 지수 평활법 기법 별 수식 ]

 

를 지수 평활법으로 예측하고자 하는 t 시점의 시계열 값이라고 하고, 

를 t 시점의 계절 요소값, 

를 t 시점에서의 오차라고 했을 때, 각 지수평활법 기법별로 모형의 수식을 표기해보면 아래와 같습니다. 지수평활법 기법 간 수식을 비교해보면 좀더 잘 이해할 수 있을 것입니다. 

 

  • 추세와 계절성이 모두 없는 단순 지수 평활법
    (Simple Exponential Smoothing)



  • 추세는 없고 고정계절변동이 있는 가법 윈터스 방법
     (Additive Winters' Method Exponential Smoothing with No Trend)



  • 추세는 없고 승법 윈터스 방법
    (Multiplicative Winters' Method Exponential Smoothing with No Trend)



  • 1차 선형 추세는 있고 계절성은 없이중 지수 평활법
    (Two Parameter Exponential Smoothing)



  • 1차 선형 추세와 고정계절변동이 있는 가법 윈터스 지수 평활법
    (Additive Winters' Method Exponential Smoothing with Linear Trend)



  • 1차 선형 추세와 확산계절변동이 있는 승법 원터스 지수 평활법
    (Multiplicative Winters' Method Exponential Smoothing with Linear Trend)



  • 2차 비선형 추세는 있고 계절성은 없는 삼중 지수 평활법
    (Three Parameter Exponential Smoothing)



  • 2차 비선형 추세와 고정계절변동이 있는 가법 윈터스 지수 평활법
    (Additive Winters' Method Exponential Smoothing with Quadratic Trend)



  • 2차 비선형 추세와 확산계절변동이 있는 승법 윈터스 지수 평활법
    (Multiplicative Winters' Method Exponential Smoothing with Quadratic Trend)

    : 

 

다음번 포스팅에서는 시계열 예측 모형의 적합도를 평가하는 지표(goodness-of-fit measures for time series models) (https://rfriend.tistory.com/667) 를 살펴보고, 이어서 Python 으로 단순 지수 평활법으로 시계열 자료를 예측하고 모델 성능을 비교평가하는 방법(https://rfriend.tistory.com/671)을 소개하겠습니다. 

 

많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~

'를 꾹 눌러주세요. :-)

 

 

728x90
반응형
Posted by Rfriend
,