Grubb's test 이상치 탐지

금융 데이터를 다루는 일을 하다보니 이상치 데이터를 거의 매번 만나게 됩니다.

보통의 경우 이상치 데이터는 변환해야할 대상이지만 금융쪽 도메인에선 이상치가 굉장히 유의미한 경우가 많습니다. 

그리고 데이터 자체도 정규성을 따르지 않는 경우도 있기 때문에 이상치 추출에 대한 여러가지 방법을 시도하는 도중 정리하기 위해 글을 쓰게 되었습니다.

 

Grubb's test라고 불리는 이 방법은 이상치를 탐지하는데 사용되는데 모집단이 정규분포임을 가정합니다. 즉 정규성 가정에 기반을 두고 있으며 테스트를 들어가기전에 정규 분포를 따르는지 확인해야 합니다.

 

Grubb's test는 한 번에 하나의 이상치만 탐지합니다. 그래서 특이치가 감지되지 않을 때 까지 계속해서 반복해야합니다. 그러나 반복을 하다보면 샘플의 갯수가 작아지고 그에 따라 탐지 확률이 바뀌고 그렇게 되면 이상치가 아닌 경우에도 이상치라고 하기 때문에 표본의 크기가 6이하일 경우에는 종료하는 것이 좋습니다. 

 

Grubb's test의 귀무가설은 "데이터셋에 이상치가 없다"입니다.

Grubb's test의 대립가설은 "데이터셋에 이상치가 하나는 존재한다"입니다.

 

Grubb's test의 Test Statistic 즉 검정 통계량은

https://en.wikipedia.org/wiki/Grubbs's_test_for_outliers

와 같습니다. 

여기서 Y¯는 표본 평균 그리고 s 는 표준 편차를 말합니다. Grubb's test의 검정 통계량은 표본 평균에서 샘플중에서 가장 큰 절대 편차를 가져오는 것입니다.

 

Critical Value 즉 임계치는

https://en.wikipedia.org/wiki/Grubbs's_test_for_outliers

와 같습니다.

 

아래는 Python3으로 구현한 Grubb's test 입니다.

우선 필요한 패키지들을 import 합니다.

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

 Jupter notebook 내에서 실행될 수 있도록 inline 처리를 합니다.  

 

%matplotlib inline

 

이상치를 포함한 배열 Y를 생성합니다. 여기서 이상치에 해당하지 않는 데이터의 경우 정규분포를 따라야 합니다. 

y = np.array([10, 11, 12, 13, 14, 15, 30, 50, 100])
x = np.arange(len(y))

산점도를 찍으면 아래와 같이 나옵니다. 

plt.scatter(x, y)

 

 

데이터셋 Y의 평균을 구합니다.

avg_y = np.mean(y)
print(avg_y)

28.333333333333332

 

Y의 평균에서 Y의 값을뺀 후 각 값들에 대해서 절대값을 취합니다.

abs_val_minus_avg = abs(y - avg_y)
print(abs_val_minus_avg)

 

[18.33333333 17.33333333 16.33333333 15.33333333 14.33333333 13.33333333
  1.66666667 21.66666667 71.66666667]

 

 

구한 값 중에서 가장 큰 값을 가져옵니다.

max_of_deviation = max(abs_val_minus_avg)
print(max_of_deviation)

 

71.66666666666667

 

Y값에서 표준 편차를 구합니다. 

 

std_y = np.std(y)
print(std_y)

 

28.14644244344607

검정 통계량을 구해보면 아래와 같습니다. 

Gcalculated = max_of_deviation / std_y
print(Gcalculated)

 

2.546206925108375

 

임계치를 계산하는 함수힙니다. 

def calculate_critical_value(size, alpha):
    """Calculate the critical value with the formula given for example in
    https://en.wikipedia.org/wiki/Grubbs%27_test_for_outliers#Definition
    Args:
        ts (list or np.array): The timeseries to compute the critical value.
        alpha (float): The significance level.
    Returns:
        float: The critical value for this test.
    """
    t_dist = stats.t.ppf(1 - alpha / (2 * size), size - 2)
    numerator = (size - 1) * np.sqrt(np.square(t_dist))
    denominator = np.sqrt(size) * np.sqrt(size - 2 + np.square(t_dist))
    critical_value = numerator / denominator
    print("Grubbs Critical Value: {}".format(critical_value))
    return critical_value

 

Gcritical = calculate_critical_value(len(y), 0.01)

Grubbs Critical Value: 2.3868098750782827

원본 코드 출처: https://github.com/bhattbhavesh91/outlier-detection-grubbs-test-and-generalized-esd-test-python/blob/master/generalized-esd-test-for-outliers.ipynb