지난번 포스팅에서는 독립된 2개 표본의 평균 차이를 검정하는 t-test 에 대해서 알아보았습니다. 이번 포스팅에서는 Python을 이용한 짝을 이룬 표본의 평균 차이를 검정하는 paired t-test 에 대해서 소개하겠습니다. (R을 이용한 짝을 이룬 표본에 대한 평균 차이 검정은 https://rfriend.tistory.com/128 를 참고하세요)
짝을 이룬 t-test(paired t-test)는 짝을 이룬 측정치(paired measurements)의 평균 차이가 있는지 없는지를 검정하는 방법입니다. 데이터 값이 짝을 이루어서 측정되었다는 것에 대해 예를 들어보면, 사람들의 그룹에 대해 신약 효과를 알아보기 위해 (즉, 신약 투약 전과 후의 평균의 차이가 있는지) 신약 투약 전과 후(before-and-after)를 측정한 데이터를 생각해 볼 수 있습니다.
[ One sample t-test vs. Independent samples t-test vs. Paired samples t-test ]
짝을 이룬 t-test(paired t-test)는 종속 표본 t-test (the dependent samples t-test), 짝을 이룬 차이 t-test (the paired-difference t-test), 매칭된 짝 t-test (the matched pairs t-test), 반복측정 표본 t-test (the repeated-sample t-teset) 등의 이름으로도 알려져있습니다. 동일한 객체에 대해서 전과 후로 나누어서 반복 측정을 합니다.
아래의 도표에는 짝을 이룬 두 표본의 대응 비교 데이터셋에 대한 모습입니다. (Xi, Yi) 가 동일한 대상에 대해서 before-after 로 반복 측정되어서, 각 동일 객체의 전-후의 차이 (즉, Di = Xi - Yi) 에 대해서 검정을 진행하게 됩니다.
paired t-test
[ 가정사항 (Assumptions) ]
(1) 측정 대상이 독립적(independent)이어야 합니다. 하나의 객체에 대한 측정은 어떤 다른 객체의 측정에 영향을 끼치지 않아야 합니다. (2) 각 짝을 이룬 측정치는 동일한 객체로 부터 얻어야 합니다. 예를 들면, 신약의 효과를 알아보기 위해 투약 전-후(before-after) 를 측정할 때 동일한 환자에 대해서 측정해야 합니다. (3) 짝을 이뤄 측정된 전-후의 차이 값은 정규분포를 따라야한다는 정규성(normality) 가정이 있습니다. 만약 정규성 가정을 충족시키지 못하면 비모수 검정(nonparametric test) 방법을 사용해야 합니다.
[ (예제) 신약 치료 효과 여부 검정 ]
새로운 당뇨병 치료제를 개발한 제약사의 예를 계속 들자면, 치료에 지대한 영향을 주는 외부요인을 통제하기 위해 10명의 당뇨병 환자를 선별하여 1달 동안 '위약(placebo)'을 투여한 기간의 혈당 (Xi)과 동일 환자에게 '신약(new medicine)'을 투여한 1달 기간 동안의 혈당 수치(Yi)를 측정하여 짝을 이루어 혈당 차이를 유의수준 5%에서 비교하는 방법이 짝을 이룬 표본에 대한 검정이 되겠습니다. (palacebo 와 신약 투여 순서는 무작위로 선정. 아래 예는 그냥 예시로 아무 숫자나 입력해본 것임. 혈당 수치 이런거 전 잘 몰라요. ^^;)
scipy 모듈의 scipy.stats.ttest_rel 메소드를 사용해서 쌍을 이룬 t-test 를 실행합니다. "Calculate the t-test on TWO RELATED samples of scores, a and b." 라는 설명처럼 TWO RELATED samples 에서 rel 을 타서 메소드 이름을 지었습니다. (저라면 ttest_paired 라고 메소드 이름 지었을 듯요...)
alternative='two-sided' 가 디폴트 설정인데요, 이번 예제에서는 'H1: 신약이 효과가 있다. (즉, 신약 먹기 전보다 신약 먹은 후에 혈당이 떨어진다)' 는 가설을 검정하기 위한 것이므로 alternative='greater' 로 설정을 해주었습니다.
## -- Paired t-test
import numpy as np
from scipy import stats
## sample data-set (repeated measurements of Before vs. After for the same objects)
bef = np.array([51.4, 52.0, 45.5, 54.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1])
aft = np.array([50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0])
## paired t-test using Python scipy module
# H0: New medicine is not effective (i.e., no difference b/w before and after)
# H1: New medicine is effective (i.e., there is difference b/w before and after)
stat, p_val = stats.ttest_rel(bef, aft, alternative='greater')
print('statistic:', stat, ' p-value:', p_val)
# statistic: 3.550688262985491 p-value: 0.003104595950799298
분석 결과 p-value 가 0.003 으로서 유의수준 0.05 하에서 귀무가설 (H0: 신약은 효과가 없다. 즉, before와 after의 차이가 없다) 을 기각(reject)하고, 대립가설(H1: 신약은 효과가 있다. 즉, before 보다 after의 혈당 수치가 낮아졌다)을 채택(accept) 합니다.
모집단의 평균과 분산에 대해서는 알지 못하는 경우가 많으므로, 보통은 모집단에서 무작위로 표본을 추출(random sampling)해서 모집단의 평균과 분산을 추정합니다. 표본의 크기가 작은 집단 간 평균의 차이가 있는지를 검정할 때 t-분포에 기반한 t-통계량(t-statistics)을 사용하여 검정을 합니다.
t-검정은 대상 표본 집단이 1개인지 2개인지에 따라서 아래와 같이 구분할 수 있습니다.
* One-sample t-test : 모집단의 평균이 귀무가설의 특정 평균 값과 같은지를 검정
* Two-sample t-test: 두 모집단의 평균이 같다는 귀무가설을 검정
One-sample t-test와 Two-sample t-test에서 사용하는 통계량에 대해서는 아래에 정리해보았습니다.
여기서부터는 독립된 두 표본 간의 평균 차이에 대한 t-검정 (independent two-sample t-test) 에 대해서만 자세하게 소개하도록 하겠습니다.
t-test 를 통해 나온 p-value 가 유의수준보다 작으면 귀모가설을 기각하고 대립가설을 채택(즉, 두 모집단의 평균이 차이가 있다)하게 됩니다.
(2) Two-sample t-test 의 가정사항 (Assumptions)
Two-sample t-test 의 결과가 유효하기 위해서는 아래의 가정사항을 충족시켜야 합니다.
(a) 한 표본의 관측치는 다른 표본의 관측치와 독립이다. (independent)
(b) 데이터는 정규분포를 따른다. (normally distributed)
(c) 두 집단의 표본은 동일한 분산을 가진다. (the same variance).
(--> 이 가설을 만족하지 못하면 Welch's t-test 를 실행합니다.)
(d) 두 집단의 표본은 무작위 표본추출법을 사용합니다. (random sampling)
정규성 검정(normality test)을 위해서 Kolmogorov-Smirnov test, Shapiro-Wilk test, Anderson-Darling test 등을 사용합니다. 등분산성 검정(Equal-Variance test) 을 위해서 Bartlett test, Fligner test, Levene test 등을 사용합니다.
(3) Python을 이용한 Two-sample t-test 실행
(3-1) 샘플 데이터 생성
먼저 numpy 모듈을 사용해서 정규분포로 부터 각 관측치 30개를 가지는 표본을 3개 무작위 추출해보겠습니다. 이중 표본집단 2개는 평균과 분산이 동일한 정규분포로 부터 무작위 추출하였으며, 나머지 1개 집단은 평균이 다른 정규분포로 부터 무작위 추출하였습니다.
다음으로 집단 x1과 x2, 집단 x1과 x3에 대한 등분산 가정 검정 결과, p-value 가 모두 유의수준 0.05 보다 크므로 두 집단 간 분산이 같다고 할 수 있습니다. (귀무가설 H0: 두 집단 간 분산이 같다.)
## (2) Equal variance test using Bartlett's tes
var_test_stat_x1x2, var_test_p_val_x1x2 = stats.bartlett(x1, x2)
var_test_stat_x1x3, var_test_p_val_x1x3 = stats.bartlett(x1, x3)
print('[x1 vs. x2]', 'statistic:', var_test_stat_x1x2, ' p-value:', var_test_p_val_x1x2)
print('[x1 vs. x3]', 'statistic:', var_test_stat_x1x3, ' p-value:', var_test_p_val_x1x3)
# [x1 vs. x2] statistic: 0.4546474955289549 p-value: 0.5001361557169177
# [x1 vs. x3] statistic: 0.029962346601998174 p-value: 0.8625756934286083
처음에 샘플 데이터를 생성할 때 정규분포로 부터 분산을 동일하게 했었으므로 예상한 결과대로 잘 나왔네요.
(3-3) 독립된 두 표본에 대한 t-test 평균 동질성 여부 검정
이제 독립된 두 표본에 대해 t-test 를 실행해서 두 표본의 평균이 같은지 다른지 검정을 해보겠습니다.
- (귀무가설 H0) Mu1 = Mu2 (두 집단의 평균이 같다)
- (대립가설 H1) Mu1 ≠ Mu2 (두 집단의 평균이 다르다)
분산은 서로 같으므로 equal_var = True 라고 매개변수를 설정해주었습니다.
그리고 양측검정(two-sided test) 을 할 것이므로 alternative='two-sided' 를 설정(default)해주면 됩니다. (왜그런지 자꾸 에러가 나서 일단 코멘트 부호 # 로 막아놨어요. scipy 버전 문제인거 같은데요... 흠... 'two-sided'가 default 설정이므로 # 로 막아놔도 문제는 없습니다.)
## (3) Identification test using Independent 2 sample t-test
## x1 vs. x2
import scipy.stats as stats
t_stat, p_val = stats.ttest_ind(x1, x2,
#alternative='two-sided', #‘less’, ‘greater’
equal_var=True)
print('t-statistic:', t_stat, ' p-value:', p_val)
#t-statistic: -0.737991822907993 p-value: 0.46349499774375136
#==> equal mean
## x1 vs. x3
import scipy.stats as stats
t_stat, p_val = stats.ttest_ind(x1, x3,
#alternative='two-sided', #‘less’, ‘greater’
equal_var=True)
print('t-statistic:', t_stat, ' p-value:', p_val)
#t-statistic: -15.34800563666855 p-value: 4.370531118607397e-22
#==> different mean
(3-1)에서 샘플 데이터를 만들 때 x1, x2 는 동일한 평균과 분산의 정규분포에서 무작위 추출을 하였으며, x3만 평균이 다른 정규분포에서 무작위 추출을 하였습니다.
위의 (3-3) t-test 결과를 보면 x1, x2 간 t-test 에서는 p-value 가 0.46으로서 유의수준 0.05 하에서 귀무가설(H0)을 채택하여 두 집단 x1, x2 의 평균은 같다고 판단할 수 있습니다.
x1, x3 집단 간 t-test 결과를 보면 p-value 가 4.37e-22 로서 유의수준 0.05 하에서 귀무가설(H0)을 기각(reject)하고 대립가설(H1)을 채택(accept)하여 x1, x3 의 평균이 다르다고 판단할 수 있습니다.
ARIMA 모형과 같은 통계적 시계열 예측 모델의 경우 시계열데이터의 정상성 가정을 충족시켜야 합니다. 따라서 만약 시계열 데이터가 비정상 확률 과정 (non-stationary process) 이라면, 먼저 시계열 데이터 변환을 통해서 정상성(stationarity)을 충족시켜주어야 ARIMA 모형을 적합할 수 있습니다.
이번 포스팅에서는 Python을 사용하여
(1) 분산이 고정적이지 않은 경우 분산 안정화 변환 (variance stabilizing transformation, VST)
(2) 추세가 있는 경우 차분을 통한 추세 제거 (de-trend by differencing)
(3) 계절성이 있는 경우 계절 차분을 통한 계절성 제거 (de-seasonality by seaanl differencing)
하는 방법을 소개하겠습니다.
[ 비정상확률과정을 정상확률과정으로 변환하기 (Transforming non-stationary to stationary process) ]
먼저 예제로 사용할 약 판매량 (drug sales) 시계열 데이터를 가져와서 pandas DataFrame으로 만들고, 시계열 그래프를 그려보겠습니다.
위의 시계열 그래프에서 볼 수 있는 것처럼, (a) 분산이 시간의 흐름에 따라 증가 하고 (분산이 고정이 아님), (b) 추세(trend)가 있으며, (c) 1년 주기의 계절성(seasonality)이 있으므로, 비정상확률과정(non-stationary process)입니다.
KPSS 검정을 통해서 확인해봐도 p-value가 0.01 이므로 유의수준 5% 하에서 귀무가설 (H0: 정상 시계열이다)을 기각하고, 대립가설(H1: 정상 시계열이 아니다)을 채택합니다.
## 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)
## 귀무가설 (H0): 정상 시계열이다
## 대립가설 (H1): 정상 시계열이 아니다 <-- p-value 0.01
kpss_test(df)
# Results of KPSS Test:
# Test Statistic 2.013126
# 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
(1) 분산이 고정적이지 않은 경우 분산 안정화 변환 (variance stabilizing transformation, VST)
분산이 고정적이지 않은 경우 멱 변환(Power Transformation)을 통해서 분산을 안정화(variance stabilization) 시켜줍니다. 분산이 고정적이지 않고 추세가 있는 경우 분산 안정화를 추세 제거보다 먼저 해줍니다. 왜냐하면 추세를 제거하기 위해 차분(differencing)을 해줄 때 음수(-)가 생길 수 있기 때문입니다.
power transformation
원래의 시계열 데이터의 분산 형태에 따라서 적합한 멱 변환(power transformation)을 선택해서 정상확률과정으로 변환해줄 수 있습니다. 아래의 예제 시도표를 참고하세요.
이번 포스팅에서 사용하는 예제는 시간이 흐릴수록 분산이 점점 커지는 형태를 띠고 있으므로 로그 변환(log transformation) 이나 제곱근 변환 (root transformation) 을 해주면 정상 시계열로 변환이 되겠네요. 아래 코드에서는 자연로그를 취해서 로그 변환을 해주었습니다.
Python의 diff() 메소드를 사용해서 차분을 해줄 수 있습니다. 이때 차분의 차수 만큼 결측값이 생기는 데요, dropna() 메소드를 사용해서 결측값은 제거해주었습니다.
## De-trend by Differencing
df_vst_diff1 = df_vst.diff(1).dropna()
df_vst_diff1.plot(figsize=(12, 8))
plt.title("De-trend by 1st order Differencing", fontsize=16)
plt.show()
de-trending by 1st order differencing
위의 시도표를 보면 1차 차분(1st order differencing)을 통해서 이제 추세(trend)도 제거되었음을 알 수 있습니다. 하지만 아직 계절성(seasonality)이 남아있어서 정상성 조건은 만족하지 않겠네요. 그런데 아래에 KPSS 검정을 해보니 p-value가 0.10 으로서 유의수준 5% 하에서 정상성을 만족한다고 나왔네요. ^^;
## 귀무가설 (H0): 정상 시계열이다 <-- p-value 0.10
## 대립가설 (H1): 정상 시계열이 아니다
kpss_test(df_vst_diff1)
# Results of KPSS Test:
# Test Statistic 0.121364
# p-value 0.100000
# Lags Used 37.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) 계절성이 있는 경우 계절 차분을 통한 계절성 제거 (de-seasonality by seaanl differencing)
아직 남아있는 계절성(seasonality)을 계절 차분(seasonal differencing)을 사용해서 제거해보겠습니다. 1년 12개월 주기의 계절성을 띠고 있으므로 diff(12) 함수로 계절 차분을 실시하고, 12개의 결측값이 생기는데요 dropna() 로 결측값은 제거해주었습니다.
## Stationary Process: De-seasonality by Seasonal Differencing
df_vst_diff1_diff12 = df_vst_diff1.diff(12).dropna()
## plotting
df_vst_diff1_diff12.plot(figsize=(12, 8))
plt.title("De-seasonality by Seasonal Differencing",
fontsize=16)
plt.show()
de-seasonality by seasonal differencing
위의 시도표를 보면 이제 계절성도 제거가 되어서 정상 시계열처럼 보이네요. 아래에 KPSS 검정을 해보니 p-value 가 0.10 으로서, 유의수준 5% 하에서 귀무가설(H0: 정상 시계열이다)을 채택할 수 있겠네요.
## 귀무가설 (H0): 정상 시계열이다 <-- p-value 0.10
## 대립가설 (H1): 정상 시계열이 아니다
kpss_test(df_vst_diff1_diff12)
# Results of KPSS Test:
# Test Statistic 0.08535
# p-value 0.10000
# Lags Used 8.00000
# Critical Value (10%) 0.34700
# Critical Value (5%) 0.46300
# Critical Value (2.5%) 0.57400
# Critical Value (1%) 0.73900
# dtype: float64
이제 비정상 시계열(non-stationary process)이었던 원래 데이터를 (1) log transformation을 통한 분산 안정화, (2) 차분(differencing)을 통한 추세 제거, (3) 계절 차분(seasonal differencing)을 통한 계절성 제거를 모두 마쳐서 정상 시계열(stationary process) 로 변환을 마쳤으므로, ARIMA 통계 모형을 적합할 수 있게 되었습니다.
지난번 포스팅에서는 백색잡음과정(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) 여부를 확인하는 방법을 소개하겠습니다.
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
## 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_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)을 꼭 그려보고 눈으로도 반드시 확인해보는 습관을 들이시기 바랍니다.
이전 포스팅에서 시계열 데이터를 예측하는 분석 방법 중의 하나로서 지수평활법(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번의 정상확률과정이 눈으로만 봐서는 구분이 잘 안가지요?
백색잡음과정, 확률보행과정, 정상확률과장이 서로 관련이 있고, 얼핏 봐서는 비슷하기도 해서 참 헷갈립니다. 정확하게 이 셋을 구분해서 이해하시기 바래요.
이번 포스팅에서는 주기성을 가지는 센서 시계열 데이터, 음향 데이터, 신호 처리 데이터 등에서 많이 사용되는 스펙트럼 분석 (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
예제로 사용할 샘플데이터로서, 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개만 가져다가 시도표로 그려보면 아래처럼 코사인 파동 모양의 주기성을 띠는 데이터임을 눈으로 확인할 수 있습니다.
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 에서 스파이크가 엄청 높게 튀는군요.
시계열 데이터를 분석할 때 꼭 확인하고 처리해야 하게 있는데요, 바로 결측값 여부 확인과 결측값 처리입니다.
시계열 데이터의 결측값을 처리하는 방법에는
(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()
그럼, 실제로 다양한 지수평활법으로 위의 약 판매량 시계열 데이터에 대한 예측 모델을 만들어 예측을 해보고, 다양한 평가지표로 모델 간 비교를 해봐서 가장 우수한 모델이 무엇인지, 과연 앞서 말한 "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를 제공합니다. 아래의 코드는 순서대로
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() 함수를 사용하며, 괄호 안에는 예측하고자 하는 기간을 정수로 입력해줍니다.
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보다 크게 되었기 때문입니다.)
시계열 모델 성능 지표별로 보면 전반적으로 7번째 지수평활법 시계열 예측 모델인 "Multiplicative Winters' method with Linear Trend" 모델이 상대적으로 가장 우수함을 알 수 있습니다. 역시 시도표를 보고서 우리가 예상했던대로 결과가 나왔습니다.
MAPE (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'] 으로 적합된 모델의 속성값에 접근할 수 있습니다.
이번 포스팅에서는 시계열 자료 예측 모형의 성능, 모델 적합도 (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)
지난번 포스팅에서는 샘플 크기가 다른 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 의 칼럼 이름으로 선택해오면 됩니다.
Greenplum DB에서 PL/Python (또는 PL/R)을 사용하여 여러개의 숫자형 변수에 대해 일원분산분석을 분산병렬처리하는 방법(one-way ANOVA in parallel using PL/Python on Greenplum DB)은 rfriend.tistory.com/640 를 참고하세요.
