'Python 분석과 프로그래밍/Python 통계분석'에 해당되는 글 13건

  1. 2021.10.31 [Python] 분산 안정화 변환과 차분으로 정상확률과정으로 변환(variance stabilization transformation and differencing for stationarity)
  2. 2021.10.10 [Python] 통계적 가설 검정을 통한 시계열 정상성 여부 확인 (checking stationarity using statistical hypothesis test: ADF test, KPSS test)
  3. 2021.10.04 [Python] 백색잡음과정 vs. 확률보행과정 vs. 정상확률과정 (White Noise Process vs. Random Walk Process vs. Stationary Process)
  4. 2021.10.03 [Python] 스펙트럼 분석 (spectral analysis, spectrum analysis)
  5. 2021.09.07 [Python] scipy를 이용한 시계열 데이터 보간 (interpolation of time series data using scipy) (3)
  6. 2021.07.29 [Python] 지수평활법 모형 훈련 및 예측, 모델평가 (Exponential Smoothing in Python) (19)
  7. 2021.07.18 [Python 시계열 자료 분석] 시계열 모형의 예측 성능 평가 지표 (7)
  8. 2021.05.08 [Python] DataFrame에서 여러개의 변수에 대해 일원분산분석 검정하기 (ANOVA test for multiple numeric variables in pandas DataFrame)
  9. 2021.05.07 [Python] 샘플 크기가 다른 2개 이상 그룹간 일원분산분석 (one-way ANOVA with different sized samples)
  10. 2020.01.12 [Python] 선형회귀모형, 로지스틱 회귀모형에 대한 각 관측치 별 변수별 기여도(민감도) 분석 (Sensitivity analysis of linear regression & Logistic regression per each variables and each observations) (1)

이전 포스팅에서는 

 (i) 정상확률과정(stationary process)의 정의 (https://rfriend.tistory.com/691)

 (ii) 통계적 가설 검증을 통한 시계열 정상성(stationarity test) 여부 확인 (https://rfriend.tistory.com/694)

하는 방법을 소개하였습니다. 

 

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으로 만들고, 시계열 그래프를 그려보겠습니다. 

 

import numpy as np
import pandas as pd
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)

#               value
# date	
# 	            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('Non-Stationary Process: Increasing Variance + Trend + Seasonality', 
          fontsize=16) 
plt.show()

non-stationary process

 

위의 시계열 그래프에서 볼 수 있는 것처럼, (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)을 선택해서 정상확률과정으로 변환해줄 수 있습니다. 아래의 예제 시도표를 참고하세요. 

variance stabilizing transformation (power transfortion)

 

이번 포스팅에서 사용하는 예제는 시간이 흐릴수록 분산이 점점 커지는 형태를 띠고 있으므로 로그 변환(log transformation) 이나 제곱근 변환 (root transformation) 을 해주면 정상 시계열로 변환이 되겠네요. 아래 코드에서는 자연로그를 취해서 로그 변환을 해주었습니다. 

 

## Variance Stabilizing Transformation (VST) by Taking Logarithm
df_vst = np.log(df.value)

df_vst.head()

# date
# 1991-07-01    1.260332
# 1991-08-01    1.157161
# 1991-09-01    1.179338
# 1991-10-01    1.283986
# 1991-11-01    1.271408
# Name: value, dtype: float64


## plotting
df_vst.plot(figsize=(12, 8))
plt.title("Variance Stabilizing Transformation by taking Logarithm", 
          fontsize=16)
plt.show()

variance stabilizing transformation (VST)

 

 

위의 시도표를 보면 시간이 경과해도 분산이 안정화되었음을 알 수 있습니다.  KPSS 검정을 한번 더 해주면 아직 추세(trend)와 계절성(seasonality)가 남아있으므로 여전히 비정상확률과정을 따른다고 나옵니다. 

 

## 귀무가설 (H0): 정상 시계열이다  
## 대립가설 (H1): 정상 시계열이 아니다  <-- p-value 0.01 

kpss_test(df_vst)

# Results of KPSS Test:
# Test Statistic           2.118189
# 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

 

 

 

(2) 추세가 있는 경우 차분을 통한 추세 제거 (de-trend by differencing)

 

차분(differencing)은 현재의 시계열 값에서 시차 t 만큼의 이전 값을 빼주는 것입니다. 

 

1차 차분 = Delta1_Z(t) = Z(t) - Z(t-1) 

2차 차분 = Delta2_Z(t) = Z(t) - Z(t-1) - (Z(t-1) - Z(t-2)) =  Z(t) - 2Z(t-1) + Z(t-2)

 

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 통계 모형을 적합할 수 있게 되었습니다. 

 

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

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

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 백색잡음과정(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/

 

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

 

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

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

 

반응형
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 

 

 

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

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

 

반응형
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 를 참고하시기 바랍니다. 

 

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

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

 

반응형
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

 

 

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

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

 

반응형
Posted by Rfriend

댓글을 달아 주세요

  1. huns 2022.01.17 15:52  댓글주소  수정/삭제  댓글쓰기

    매번 도움을 많이 받습니다.~
    질문이 있어서요~~
    xnew 값은 무엇인가요?

    • Rfriend 2022.01.17 23:56 신고  댓글주소  수정/삭제

      안녕하세요.
      scipy.interpolate.interp1d() 함수는 기존의 데이터(포스팅 본문에 있는 그래프의 점)로부터 "보간을 할 수 있는 함수"를 적합합니다.

      그리고 "보간을 할 수 있는 함수"에 보간을 하고 싶은 x_new 값을 대입해주면, 이 함수에 계산이 된 보간 후의 값을 반환하게 됩니다.

      가령, 기존 값이 x=1 일때 y = 10, x=5 일때 y=50 인 값으로 선형보간(linear interpolation) 하는 함수(y= 10 * x)를 적합한 후에, x_new = [1, 2, 3, 4, 5] 를 대입해주면 y_new = [10, 20, 30, 40, 50] 을 반환하는 식입니다.

  2. huns 2022.01.18 23:12  댓글주소  수정/삭제  댓글쓰기

    답변 말씀 감사합니다. ~
    그런데 한가지 더 질문이 있습니다. python을 이용해서 이동평균을 통한 데이터 보간 하는 방법이 궁금합니다. 계속 구글 검색하니깐 나오지가 않아서요~

    R 같은 경우는 imputeTS::na_ma 를 사용하면 되더라구요 파이썬같은경우는 패키지가 있을까요?
    이동평균이나 arima를 통한 보간법이 파이썬에도 있는지 아니면 직접 구현을 해야 하는지 궁금합니다.

지난번 포스팅에서는 시계열 패턴 별 지수평활법 모델을 설정하는 방법 (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

 

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

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

 

반응형
Posted by Rfriend

댓글을 달아 주세요

  1. 최새나 2021.08.05 20:11  댓글주소  수정/삭제  댓글쓰기

    rfriend님 정말 감사합니다! 백신맞고 좀 아팠다가 괜찮아져서 오늘 출근했습니다ㅠㅠ 올려주신글 읽고 싶어서 집에서 누워서 노트북키구 여러번 읽어보고 그랬습니다ㅠㅠ jupyter Notebook으로 직접해봤는데 너무 재미있어요! 알려주신내용 익혀서 회사 Sensor data로 적용시켜보려구요!! python 안에 지수 평활법을 구현할 수 있는 statsmodel 패키지가 있다는 것을 덕분에 알게되었습니다! 정말정말 감사합니다^^

    • Rfriend 2021.08.06 11:09 신고  댓글주소  수정/삭제

      백신 맞고서 아프셨다니 고생 많으셨겠어요. 그래도 백신 맞아서 한스름 놓으셨을듯 해요.
      이반 포스팅이 도움이 되었기를 바래요.

  2. 최새나 2021.08.06 19:06  댓글주소  수정/삭제  댓글쓰기

    forcast는 month 단위로만 가능한가요?ㅠㅠ 저는 일별 예측을 하고싶은데, raw data를 월단위로 끊어서만 넣어줘야 분석이 가능한거죠?

    • Rfriend 2021.08.06 19:48 신고  댓글주소  수정/삭제

      훈련 데이터 (raw data)가 일별 간위이면 forecast 도 일별 단위입니다.
      월별은 일 단위 예측을 월단위 합계 내거나 아니면 raw data를 월별로 집계 후 그걸로 모델 훈련 시킨 후 월별 예측해도 되구요.

  3. 엉클창 2021.08.30 16:46 신고  댓글주소  수정/삭제  댓글쓰기

    다음과 같은 에러 메시지가 뜨는데 이유가 뭘까요? 다른 데이터로 공부해 보려고 엑셀파일을 올려서 분석하는데 이 메시지가 뜹니다
    TypeError: __init__() got an unexpected keyword argument 'initialization_method'

  4. 익명 2021.08.31 09:24  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • Rfriend 2021.08.31 13:39 신고  댓글주소  수정/삭제

      안녕하세요 엉클창님. 반갑습니다.

      독학으로 Python 공부하신다니 힘드실테지만 응원 보내드립니다!

      (1) statsmodels 모듈의 버전을 확인해보세요. initialization_method 옵션은 버전 0.12 부터 추가되었습니다. 저는 version 0.12.2 를 쓰고 있네요. 혹시 0.12 이전 버전 사용하고 계신다면 최신 버전으로 업그레이드 하세요.

      -- 커맨드라인 터미널에서 모듈 버전 확인
      (base) lhongdon@Hongui-MacBookPro ~ % pip show statsmodels
      Name: statsmodels
      Version: 0.12.2
      Summary: Statistical computations and models for Python
      Home-page: https://www.statsmodels.org/
      Author:
      Author-email:
      License: BSD License
      Location: /Users/lhongdon/opt/anaconda3/lib/python3.8/site-packages
      Requires: patsy, pandas, scipy, numpy
      Required-by: pmdarima


      (2) statsmodels 의 메소드는 인풋으로 시계열 형태의 index 와 숫자형 칼럼으로 구성된 데이터셋을 사용합니다. 올려주신 코드를 보니 df_A2에 (a) index 는 Month 로 잘 하셨구요, (b) 칼럼이 숫자형, 범주형으로 해서 여러개가 있는거 같아요. 모델 훈련/예측하려는 숫자형 칼럼 단 1개만 남겨두셔야 이후 에러가 안나요.

      (3) 제가 본문에 올려놓은 샘플 데이터로 한 사이클 분석해보시기 바랍니다. 동일한 에러가 발생하면 statsmodels 모듈의 version 문제일 가능성이 높구요, 동일 에러가 발생하지 않는다면 지금 사용하고 있는 데이터 포맷이나 아니면 코딩에 문제가 있을 수 있겠네요.

    • 익명 2021.08.31 14:22  댓글주소  수정/삭제

      비밀댓글입니다

  5. 손지용 2021.10.12 13:34  댓글주소  수정/삭제  댓글쓰기

    eval_all_df = pd.DataFrame(
    {'SES': eval_all(y_test, forecast_1, T, fit1),
    이 부분에서 이러한 에러가 뜨는데 이유를 알 수 있을까요?
    name 'r2' is not defined

    • Rfriend 2021.10.12 13:37 신고  댓글주소  수정/삭제

      안녕하세요.

      에러 메시지("name 'r2' is not defined")를 보니
      def ADJ_R2() 를 정의할 때
      ==> r2 = 1 - sse/sst
      이 라인을 빼먹은게 아닌가 싶습니다.

    • 손지용 2021.10.12 16:28  댓글주소  수정/삭제

      감사합니다!
      혹시 하나 더 질문해도 될까요?
      fit1 = SimpleExpSmoothing(df_train, initialization_method="estimated").fit()
      이 부분에서
      에러메시지 : Pandas data cast to numpy dtype of object. Check input data with np.asarray(data).
      다음과 같이 에러가 발생하는 이유를 알 수 있을까요?

    • Rfriend 2021.10.12 16:52 신고  댓글주소  수정/삭제

      pandas DataFrame을 numpy 로 변환하는 과정에서 에러가 발생했나 보네요.

      fit1 = SimpleExpSmoothing(np.asarray(df_train), initialization_method="estimated").fit()

      로 직접 numpy array 로 변환해서 해보실래요?

  6. 파이팅 2022.02.22 16:25  댓글주소  수정/삭제  댓글쓰기

    안녕하세요 R프랜드님! 정말 감사히 잘 읽었고, 실제로 잘 적용해보았습니다.
    추가로 여쭤보고 싶은 것이 있는데요!
    가장 적합한 모델을 찾아낸 이후에, 해당 모델을 통해 향후 10년을 예측해보고 싶은데,
    코드를 어떻게 작성해야할까요?
    제가 독학 입문자라 많이 부족합니다 ㅠㅠ 죄송합니다.

    • Rfriend 2022.02.22 17:26 신고  댓글주소  수정/삭제

      안녕하세요.

      model_name.forcast(예측하려는 기간) 함수를 사용하시면 됩니다.

      예측하려는 기간에는, 가령 월 단위 학습을 했고, 1년을 예측하고 싶으면 forcast(12)를 입력해주면 됩니다.

      단, 만약 학습에 사용한 데이터가 5년치 밖에 안되는데 향이후 10년치를 예측하려고 하면 아마 동일한 값이 쭈욱 이측되는갈로 나올거예요. 시간이 뒤로 갈수록 예측값ㅇ 분산도 엄청 커지구요. 학습에 사용한 데이터 기간 대비 너무 장기 예측은 무리일거예요.

    • 데이터마케터 2022.02.22 19:28 신고  댓글주소  수정/삭제

      답변 정말 감사합니다. R프랜드님!!
      잘 되는지 확인하고 싶어, 위에 작성해주신 코드에 그대로 이어서
      forecast_7.forcast(12)라고 작성하였는데
      AttributeError: 'Series' object has no attribute 'forcast' 이란 경고 문구가 나오면서 실행이 되지 않고 있습니다.
      혹시나 모델명을 잘못 알고 있나 싶어 fit7로도 바꿔보고 하였지만, 여전히 작동을 안합니다.
      이런 것까지 여쭤봐 정말 죄송합니다만, 어떻게 적어야 할까요? ㅜㅜ

    • Rfriend 2022.02.22 19:32 신고  댓글주소  수정/삭제

      Training set에 대해 fit() 함수를 사용해서 적합한 모델에 대해서 forecast() 함수를 적용허면 됩니다.

      에러 메시지를 보니 forecast_7 은 작합된 모델 객체가 아니라 그냥 (아마도 예측한?) 숫자가 들어있는 pandas Series 이네요.

    • 데이터마케터 2022.02.22 19:48 신고  댓글주소  수정/삭제

      감사합니다!! 해결했습니다!!

    • Rfriend 2022.02.22 19:49 신고  댓글주소  수정/삭제

      아, 다행이네요! :-)

지난번 포스팅에서는 시계열 패턴별 지수 평활법 (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)

 

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

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

 

반응형
Posted by Rfriend

댓글을 달아 주세요

  1. 최새나 2021.07.20 10:54  댓글주소  수정/삭제  댓글쓰기

    댓글로 소통할 수 있고, 포스팅 올라오길 기다리고, 포스팅 올라왔을 때 순간 순간들이 참 설레네요! 감사합니다! 회사에서 Sensor data들의 자기 자신의 과거 data(성능좋았던) 비교했을때 대비 얼마나 벗어나있는지 정도를 예측하여 관리해보고자 하거든요^^ 앞으로의 포스팅하나하나가 너무 기대됩니다 포스팅해주시는대로 따라해보고있어요! 정말 큰 도움을 받고있습니다. 감사합니다~~

    • Rfriend 2021.07.20 15:32 신고  댓글주소  수정/삭제

      안녕하세요 최새나님.

      댓글에 제가 약속해놓은게 있는지라... 주말에 시골 부모님댁 다녀와서는, 일요일 오전에 커피숍 가셔 블로그 포스팅 먼저 했어요. 저도 덕분에 부담감이 생기면서 정신차리고 포스팅하게 되네요.

      다변량분석 중에 mahalanobis distance (https://rfriend.tistory.com/199), Hotelling's t-square 통계량도 한번 찾아보시면 지금 하시는 일에 도움이 될 듯 합니다.

    • 최새나 2021.07.20 18:37  댓글주소  수정/삭제

      감동입니다ㅠㅠ 좋은정보 알려주셔서 감사합니다!

  2. 최새나 2021.07.28 16:40  댓글주소  수정/삭제  댓글쓰기

    R Friend님 안녕하세요^^
    다음번 포스팅 손꼽아 기다리구있습니다!! 이들 지표를 사용해서 하나의 시계열 자료에 대해 여러 개의 지수 평활법 기법들을 적용해서 가장 모형 적합도가 높은 모델을 찾아보는 코드 너무 궁금합니다ㅠ_ㅠ 실제 Sensor 데이터 그래프별로 어떤 시계열 모형 타입인지 지수 평활법 기법들을 모형적합도따져가면서 구분해보고싶습니다!! 바쁘신와중에 항상 감사합니다~~!!

    • Rfriend 2021.07.28 18:08 신고  댓글주소  수정/삭제

      아, 네. 제가 내일 저녁까진 포스팅 해볼게요. 주말에 영화감상문 포스팅 쓰다가 시간이 다 가서 더 못썼어요. 기다리실것 같긴 했는데 늦어져서 죄송해요. ^^;;;

    • 최새나 2021.07.28 18:31  댓글주소  수정/삭제

      포스팅해주신 영화 여주인공 제가 좋아하는배우라서 저도 주말에 챙겨보려구요^^ 앗 내일저녁까지요 ㅠㅠ 너무 감사합니다!!

    • Rfriend 2021.07.30 17:25 신고  댓글주소  수정/삭제

      안녕하세요 최새나님.

      어제 밤에 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 의 칼럼 이름으로 선택해오면 됩니다. 

 

# 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

 

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

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

 

반응형
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

반응형
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)




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

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



반응형
Posted by Rfriend

댓글을 달아 주세요

  1. Stat_in_KNU 2020.06.08 16:21 신고  댓글주소  수정/삭제  댓글쓰기

    선생님 안녕하십니까? 데이터분석전문가 시험을 준비중인 학생입니다. 시험 관련 자료를 찾던 중 우연히 들어왔는데 예전에도 꽤 방문한 흔적이 있던 블로그 더군요 매번 도움 받을 수 있어서 감사한 마음입니다. 다름이 아니라, 제가 공부한 것들을 정리하는 용도로 블로그를 운영하고있는데 선생님의 게시글을 참고해서 올리고 싶습니다. 블로그를 통한 수익은 전혀 없으며 조회수 또한 거의 없습니다. 혹시 괜찮으시다면 제 블로그에 포스팅해도될까요?