파이썬 z-score 이상치 제거 - paisseon z-score isangchi jegeo

Question

이상치 (anomaly)란 주어진 데이터 분포 중심에서 멀리 떨어진 데이터를 말합니다. 말 그대로 정상 데이터가 아니라 비정상 데이터인 것이죠. 주어진 데이터에서 이상치를 찾는 가장 간단한 방법은 Z-score 입니다.

목차 Show

Modified Z-score
Interquantile range

Z-score

Z-score 는 평균과 표준오차가 정의되어 있을 떄 해당 데이터가 얼마나 벗어나 있는지 측정하는 지표로

$Z-score = \frac{x_i-\mu}{\sigma}$

와 같이 정의됩니다. 평균에서 얼마나 떨어져 있는지 계산하고 표준 오차로 나눠줌으로써 평균에서 어느 방향으로 얼마나 떨어져 있는지 계산합니다. Z-score의 절대값이 클수록 이상치라고 생각할 수 있습니다.

Example

기본적인 라이브러리를 iport 하고 1950년부터 열린 월드컵 최다 득점자 정보를 가진 데이터를 생성합니다.

import numpy as np
import pandas as pd

import scipy.stats as ss
import matplotlib.pyplot as plt

year = list(range(1950, 2019, 4))
players = list(map(chr, range(97, 97+len(year))))
goals = [8,11,13,4,9,10,7,6,6,6,6,6,6,8,5,5,6,6]                      
assert len(year) == len(players) == len(goals)

wc_goals = pd.DataFrame(list(zip(year, players, goals)), columns=['Year', 'Player', 'Goals'])

Z-score를 위의 수식대로 계산하고 (scipy.stats의 zscore 함수)

wc_goals['Zscore'] = ss.zscore(wc_goals['Goals'])

Z-score의 절대값이 2보다 크면 이상치라고 가정했을 때, 이상치가 얼마나 나오는지 살펴보면,

def plot_anomaly_goals(score_data, threshold, title='Z-score'):
    score_data = score_data.sort_values(ascending=False).values
    ranks = np.linspace(1, len(score_data), len(score_data))
    mask_outlier = (score_data > threshold)

    plt.figure()
    plt.plot(ranks[mask_outlier], score_data[mask_outlier], 'o', color='r', label='anomalies')
    plt.plot(ranks[~mask_outlier], score_data[~mask_outlier], 'o', color='b', label='typical') 
    plt.axhline(threshold, color='r', label='threshold', alpha=0.5)
    plt.legend(loc='upper right')
    plt.title(title, fontweight='bold')
    plt.xticks(np.arange(0, len(score_data), step=2))
    plt.xlabel('Player Rank')
    plt.ylabel('Z-score')
    plt.show()   
    
plot_anomaly_goals(wc_goals['Zscore'], 2)

하나의 데이터만 이상치가 됩니다. 하지만 위의 그림을 보면 Z-score 대부분이 0 근방에 위치한 것으로 보아 이상치는 더 발견되어야 할 것 같습니다. 이는 Z-score 계산에 사용된 평균과 표준 오차의 함정인데, 평균은 이상치에 의해 굉장히 크게 영향을 받기 떄문에 평균이 지나치게 증가하여 하나만 이상치로 추출되었다고 생각할 수 있습니다.

Modified Z-score

이를 해결하기 위해 modified Z-score를 사용합니다. Modified Z-score는 Median Absolute Deviation (MAD) 를 사용하는 것으로 다음과 같이 정의됩니다.

$\frac{x_i - \tilde{X}}{MAD}$

$\tilde{X}$는 데이터 $X$의 중간값 (median)이고 MAD는 $|x_i - \tilde{X}|$의 중간값으로 정의됩니다. 이때, MAD를 표준 오차를 대신해 사용하기 위해 consistent scale factor $k$를 다음과 같이 정의합니다.

$\Phi^{-1}$은 정상 분포의 CDF의 역함수로 파라미터로 들어간 3/4는 $\pm MAD$가 정상 분포의 CDF가 50%만큼 포함할 수 있도록 한 것입니다.

따라서 $\Phi(MAD/\sigma) - \Phi(-MAD/\sigma) = 1/2$가 성립해야 하고 $\Phi(-MAD/\sigma)$ =$1-\Phi(MAD/\sigma)$ (대칭 확률 분포의 CDF 특징) 임을 이용하면 $MAD/\sigma = \Phi^{-1}(3/4)$ 임을 알 수 있습니다. 최종적으로 modified Z-score는 다음과 같이 정의됩니다.

$\frac{x_i - \tilde{X}}{k*MAD}$

또한, 다음과 같이 쉽게 구현할 수 있습니다.

def modified_z_score(my_data):
    ## First Calculate Median
    median_my_data = np.median(my_data)
    
    ## Median Absolute Deviation 
    ## Median of | X_i - median of X| for all X_i
    mad = np.median(my_data.map(lambda x: np.abs(x - median_my_data)))
    
    ## Modified Z score
    ## 0.6745 * (X_i - median of X)/Median Absolute Deviation
    modified_z_score = list(my_data.map(lambda x: 0.6745* (x - median_my_data)/mad))
    return modified_z_score

Modified Z-score를 사용한 결과 4명의 이상 데이터를 발견했습니다. 표준 오차는 2이고 $k$*MAD는 1.45 정도로 이상치는 표준 오차에 큰 영향을 미친 것을 알 수 있으며, 중앙값을 이용하여 이상치가 미치는 영향을 감소시킨 것을 알 수 있습니다.

Interquantile range

Modified Z-score 이외에 interquantile range를 통해서 이상치를 찾을 수 있습니다. 먼저, quantile이란 데이터를 크기 별로 줄세운 이후 해당 퍼센트에 따른 값으로 interquantile range는 데이터의 상위 25% 값에서 하위 25% 값을 뺀 범위를 말합니다. q1 (하위 25%), q3 (상위 25%) 에서 interquantile range의 1.5배가 넘어서는 것을 이상치로 판단합니다.

def get_lower_upper_bound(my_data):
    # Get first and third quartile
    q1 = np.percentile(my_data, 25)
    q3 = np.percentile(my_data, 75)
    
    # Calculate Interquartile range
    iqr = q3 - q1
    
    # Compute lower and upper bound
    lower_bound = q1 - (iqr * 1.5)
    upper_bound = q3 + (iqr * 1.5)
    
    return lower_bound, upper_bound
    
def get_outliers_iqr(my_data):
    lower_bound, upper_bound = get_lower_upper_bound(my_data)
    # Filter data less than lower bound and more than upper bound
    return my_data[np.where((my_data > upper_bound) |
                            (my_data < lower_bound))]