이번 포스팅에서는 데이터 속에 이상치가 있는지 없는지 테스트해볼 수 있는 통계적 방법과 판별 기준에 대해서 소개한다. 여기서 다루는 내용은 다음과 같다.
5. Generalized Extreme Studentized Deviation Test
1. Grubbs's Test
- 정의 -
Grubbs's Test는 정규분포를 따르는 데이터에서 하나의 이상치를 발견할 수 있는 검정 방법이다.
Grubbs's Test는 관측된 데이터가 정규분포로 추출되거나 데이터의 분포가 근사적으로 정규분포를 따른다고 가정한다. 따라서 이 검정법을 수행하기 위해서는 사전에 데이터의 분포가 정규분포를 따른다고 볼 수 있는지 확인해야 한다.
- 가설 -
Grubbs's Test가 검정하고자 하는 가설은 다음과 같다.
$H_0$ : 데이터에 이상치가 하나도 없다.
$H_a$ : 이상치가 하나 있다.
- 검정 방법 -
$y_1, y_2, \ldots, y_n$을 현재 관측 데이터라고 할 때 Grubbs's Test의 검정 통계량은 다음과 같다.
$$G = \frac{\max_{i=1, \ldots, n}|y_i-\bar{y}|}{s}$$
여기서 $\bar{y}$는 표본평균, $s$는 표본 표준편차이다.
1) 양측 검정
만약 유의 수준 $\alpha$에 대하여 $G$가 다음을 만족한다면 귀무가설을 기각한다.
$$G > \frac{n-1}{\sqrt{n}}\sqrt{\frac{t_{\alpha/2n, n-2}^2}{n-2+t_{\alpha/2n, n-2}^2}}\tag{1}$$
여기서 $t_{\alpha/2n, n-2}^2$는 자유도가 $n-2$인 $t$분포에서 오른쪽 꼬리 확률이 $\alpha/2n$이 되게 하는 값이다.
2) 단측 검정
만약 데이터의 최소값이 이상치인지 아닌지 확인해보고자 한다면 검정 통계량 $G$를 다음과 같이 설정한다.
$$G = \frac{\bar{y}-\min_{i=1, \ldots, n}y_i}{s}$$
만약 데이터의 최대값이 이상치인지 아닌지 확인해보고자 한다면 검정 통계량 $G$를 다음과 같이 설정한다.
$$G = \frac{\max_{i=1, \ldots, n}y_i-\bar{y}}{s}$$
이때 $G$가 다음을 만족한다면 귀무가설 $H_0$를 기각한다.
$$G > \frac{n-1}{\sqrt{n}}\sqrt{\frac{t_{\alpha/n, n-2}^2}{n-2+t_{\alpha/n, n-2}^2}}$$
- 장단점 -
1) 장점
Grubbs's Test는 직관적이고 구현 과정이 간단하다.
2) 단점
Grubbs's Test는 한 번에 하나의 이상치만을 발견할 수 있다. 그렇다면 여러 개의 이상치를 확인하기 위하여 이상치를 하나 제거하고 Grubbs's Test를 수행하고 이 과정을 반복해볼 수 있을 것이다. 하지만 이렇게 하는 경우 매 스텝마다 데이터가 하나씩 제거되며 이에 따라 원 데이터 관점에서는 이상치가 아니었지만 제거된 데이터의 입장에서는 이상치가 될 수 있다. 다시 말하면 이러한 과정은 불필요한 이상치를 많이 생산하게 된다.
데이터의 분포가 정규분포로 나왔다는 뚜렷한 증거가 없다면 사용할 수 없다.
또한 데이터 개수가 적으면 사용할 수 없다. 위키에 따르면 데이터가 6개 이하라면 사용할 수 없다고 한다(사실 이 정도의 데이터 개수라면 어떤 테스트도 수행할 수 없을 것 같다).
마지막으로 이상치가 여러 개 있다면 $G$값이 작아져 이상치가 있음에도 불구하고 귀무가설을 기각 못하게 되는 상황이 발생한다. 따라서 이상치를 제대로 탐지하지 못하게 된다.
- 구현 코드 -
필요한 모듈을 임포트하자.
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import t
from scipy.stats import probplot
데이터를 만들어준다.
data = [199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57]
data = np.array(data)
분포를 확인해보자.
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
정규분포 가정을 확인해보기 위하여 Q-Q plot을 그려본다. 정규분포 가정을 테스트하는 통계적인 방법이 있지만 추천하지 않는다. 왜냐하면 이상치가 포함되어있을 경우 정규분포라는 귀무가설을 기각할 가능성이 높아지고 이에 따라 당연히 해당 데이터 분포가 정규분포가 아니라는 결과가 나오기 때문이다.
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
ax = fig.add_subplot()
probplot(data, plot=ax)
plt.show()
맨 마지막 이상치로 판단되는 점을 제외하면 직선 형태이므로 이는 주어진 데이터가 정규분포라고 생각해도 괜찮다. 다음으로 Grubbs's Test를 수행하는 함수를 만들었다. 이 함수는 이상치라고 판단되는 데이터가 있다면 해당 데이터 값과 인덱스, 통계량, critical value를 출력한다.
def grubbs_test(data,two_sided=True,alpha=0.05):
mean = np.mean(data)
std = np.std(data,ddof=1)
constant = (len(data)-1)/np.sqrt(len(data))
if two_sided:
diff = np.abs(data-mean)
max_diff = np.max(diff)
G = max_diff/std
sig_alpha = alpha/(2*len(data))
t_critical = t.ppf(1-sig_alpha,len(data)-2)
critical_val = constant*np.sqrt(np.square(t_critical)/(len(data)-2+np.square(t_critical)))
if G > critical_val:
idx = np.where(diff==max_diff)[0][0]
suspected_outlier = data[idx]
outlier_idx = idx
return (outlier_idx, suspected_outlier, G, critical_val)
else:
return (G, critical_val)
else: ## one sided
## first test whether max point is outlier
result = []
G = (np.max(data)-mean)/std
sig_alpha = alpha/(len(data))
t_critical = t.ppf(1-sig_alpha,len(data)-2)
critical_val = constant*np.sqrt(np.square(t_critical)/(len(data)-2+np.square(t_critical)))
if G > critical_val:
idx = np.where(x==np.max(data))[0][0]
suspected_outlier = data[idx]
outlier_idx = idx
result.append((outlier_idx, suspected_outlier, G, critical_val))
## next, test whether min point is outlier
G = (mean-np.min(data))/std
if G > critical_val:
idx = np.where(x==np.min(data))[0][0]
suspected_outlier = data[idx]
outlier_idx = idx
result.append((outlier_idx, suspected_outlier, G))
return result
Grubbs's Test를 이용하여 이상치를 판단해보자. 아래 코드는 검정 결과와 이상치라 판단되는 점을 표시한 것이다.
res = grubbs_test(data,two_sided=True,alpha=0.05)
if len(res) == 4:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print(res)
print('There is no outlier')
이상치를 잘 탐지해낸 것 같다.
지금까지는 이상치가 하나 있는 데이터를 다뤘지만 이번엔 이상치가 여러 개 있는 데이터를 고려해보자. 다음과 같이 데이터를 바꾸고 동일하게 실험해보자.
data = [9.1, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
data = np.array(data)
이 데이터에는 이상치가 3개 있다.
Q-Q plot을 확인해보니 이상치를 제외하면 직선 형태이므로 데이터가 정규분포로부터 추출되었다고 볼 수 있다.
하지만 아래 코드를 수행하면 이상치가 하나도 없다고 나온다. 그 이유는 이상치가 여러 개가 있어서 표준편차 값이 커지게 되고 이에 따라 검정 통계량 값이 작아졌기 때문이다.
res = grubbs_test(data,two_sided=True,alpha=0.05)
if len(res) == 4:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print(res)
print('There is no outlier')
검정통계량 값이 2.3891 정도 나왔다. 그렇다면 이상치 2개를 제외하고 하나만 남겨놓은 상태에서 Grubbs's Test를 수행해보자.
temp_data = [9.1, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
temp_data = np.array(temp_data)
res = grubbs_test(temp_data,two_sided=True,alpha=0.1)
if len(res) == 4:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(temp_data)), temp_data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(temp_data)), range(1,len(temp_data)+1))
plt.show()
else:
print(res)
print('There is no outlier')
위와 같이 이상치를 잘 잡아낼 수 있다. 이때 검정통계량 값을 확인해보자.
print(res[2])
아까 전에 2.3891보다 더 큰 값인 3.8608이다. Grubbs's Test는 이상치가 여러 개 있을 때는 사용할 수 없다고 생각하자.
2. Chauvenet's Criterion
- 정의 -
Chauvenet's Criterion은 표준 정규분포를 따르는 데이터에서 이상치를 판단하기 위한 하나의 기준이다. 이 기준은 아래 그림과 같이 평균을 중심으로 면적 $1-1/2n$을 넘어가는 데이터를 이상치로 판단한다.
여기서 $1/2n$이라는 숫자의 의미를 곰곰이 생각해보았다. 먼저 $n=1$인 경우 $1-1/2n=1/2$가 된다. 이때에는 $1/2$의 확률로 이상치를 판단하겠다는 뜻이다. 이는 굉장히 합리적이다. 왜냐하면 하나의 데이터는 평균(여기서는 평균이 0)으로부터 가깝거나 멀거나 둘 중에 하나이므로 1/2보다 작다면 이상치로 판단하는 것이 합리적이기 때문이다.
이제 데이터가 많은 경우 $n>1$인 경우를 살펴보자. 데이터가 많다면 이상치 판단을 위한 구간은 점점 넓어져야 할 것이다. 즉, 데이터 개수 $n$이 커짐에 따라 구간은 넓어져야할 것이며 이는 평균을 중심으로 한 면적이 넓어져야 한다는 뜻이다. 이러한 원리를 반영하기 위해 $1/n$이라는 인자를 넣었고 면적 1-$1/2n$을 넘어가는 데이터를 이상치로 판단하겠다는 것이 Chauvenet's Criterion에 담긴 배경이라 할 수 있다.
- 판별법 -
데이터 $y_1, y_2, \ldots, y_n$ 이 주어졌다고 했을 때 Chauvenet's Criterion은 다음과 같은 프로세스를 통해 이상치를 판단한다.
1 단계) 데이터를 다음과 같이 변환한다.
$$y_1^*, y_2^*, \ldots, y_n^*$$
여기서 $y_i^*=|y_i-\bar{y}|/s$, $\bar{y}$는 표본 평균, $s$는 표본 표준편차이다.
2 단계) 각 데이터에 대하여 표준 정규분포 누적확률을 계산한다.
$$p_i=P(Z\leq y_i^*) \;\;\forall i=1,2,\ldots,n$$
여기서 $Z$는 표준 정규분포를 따르는 확률 변수이다.
3 단계) 만약 $p_i > 1-1/2n$이라면 해당 데이터를 이상치로 판단한다.
4 단계) 이상치로 판단되었다면 제거한 후 1~3 단계)를 반복한다.
어떤 분들은 3 단계)에서 끝내는 것을 권장하지만 필요에 따라 반복해도 괜찮다고 한다.
- 장단점 -
1) 장점
Grubbs's Test에서 최대값 관련 통계량을 사용하지 않으므로 하나 이상의 이상치를 탐지할 가능성이 있으며 이상치가 여러 개인 경우에도 이상치를 발견할 수 있다. 구현 과정이 간단하다.
2) 단점
데이터의 분포가 정규분포가 아니면 사용하기 어렵다. 또한 bimodal 형태의 데이터 분포를 갖고 있다면 굉장히 안 좋다고 한다. 또한 면적 $1-1/2n$에 대한 이론적인 타당성에 대하여 지적을 많이 받고 있다.
- 구현 코드 -
먼저 필요한 모듈을 임포트해준다.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
데이터는 앞에서 다룬 이상치 3개를 포함하는 것을 사용한다.
data = [9.1, 79.5, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 78.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
data = np.array(data)
이제 Chauvenet's Criterion에 따라서 이상치로 판단된 데이터와 해당 인덱스를 출력하는 함수를 정의한다.
def chauvenet_criterion(data):
mean = np.mean(data)
std = np.std(data,ddof=1)
standardized_data = np.abs((data-mean)/std)
p_val = (1-norm.cdf(standardized_data)) < 1/(2*len(data))
if len(np.where(p_val==True)[0]) > 0:
idx = np.where(p_val==True)[0]
value = data[idx]
return (idx, value)
else:
return False
이제 이상치를 탐지해보자.
res = chauvenet_criterion(data)
if res:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print('There is no outlier')
위 코드를 실행하면 아래 그림과 같이 하나의 이상치가 탐지되었다는 것을 알 수 있다. 앞에서 Grubbs's Test로는 탐지할 수 없었던 이상치를 Chauvenet's Criterion에 의해서 탐지해냈다.
물론 탐지는 했지만 하나밖에 탐지 못하였으므로 반복 수행해준다. 나는 4번 반복 수행하여 이상치를 잘 찾아내는지 확인해보았다. 반복수행을 4번으로 한 이유는 3번을 반복하여 우리가 원하는 3개를 찾고 나머지 이상치가 아닌 데이터에 대해서 이상치라고 판단하는지 확인해보기 위함이다.
temp_data = data.copy()
outlier_idx = []
max_iter = 4
iter_num = 1
while iter_num <= max_iter:
res = chauvenet_criterion(temp_data)
if res:
so = res[1][0]
idx = np.where(data == so)[0][0]
outlier_idx.append(idx)
temp_data = np.delete(temp_data,res[0][0])
else:
break
iter_num += 1
outlier_idx = np.array(outlier_idx)
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(outlier_idx,data[outlier_idx],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
정확하게 3개를 찾아내었다.
3. Peirce's Criterion
- 정의 -
Peirce's Criterion는 Chauvenet's Criterion 보다 더 엄밀하게(rigorous) 이상치를 탐지할 수 있는 방법이라고 한다. 이 기준도 데이터가 정규분포를 따른다는 가정을 하고 있으며 표준화된 데이터의 절대값이 특정 기준값보다 크면 기각한다. Peirce's Criterion은 이상치로 의심되는 데이터가 $m$개 있는 경우 이를 제외하지 않았을 경우의 에러 분포 확률과 제외했을 경우의 에러 분포 확률을 비교하여 전자가 작다면 기각한다.
- 판별법 -
데이터 $y_1, y_2, \ldots, y_n$ 이 주어졌다고 했을 때 Peirce's Criterion은 다음과 같은 프로세스를 통해 이상치를 판단한다.
1 단계) 데이터를 다음과 같이 변환한다.
$$y_1^*, y_2^*, \ldots, y_n^*$$
여기서 $y_i^*=(y_i-\bar{y})/s$, $\bar{y}$는 표본 평균, $s$는 표본 표준편차이다.
2 단계) 특정 기준값 $R$을 찾는다. $R$을 찾기 위해 필요한 변수는 기각할 이상치의 개수 $m$이며 처음에는 $m=1$로 설정한다. $R$값은 아래에 나오는 코드로 찾을 수 있다. 수학적으로 어떻게 찾는지 알고 싶다면 이 논문을 참고하기 바란다.
3 단계) $y_1^*, y_2^*, \ldots, y_n^*$에서 $R$보다 큰 $y_i^*$를 찾고 찾은 개수를 $m'$이라 하고 4 단계를 진행한다. 이때 $y_i^*$에 해당하는 $y_i$가 이상치라고 판단한다. 만약 하나도 찾지 못했다면 프로세스를 종료한다.
4 단계) $m = 1+m'$로 설정하고 $R$을 업데이트한다.
5 단계) 3~4 단계를 반복한다.
- 장단점 -
1) 장점
Peirce's Criterion은 Chauvenet's Criterion에 비하여 엄밀하다고 한다. 또한 여러 개의 이상치를 잡아낼 수 있다. 구현도 어렵지 않다.
2) 단점
데이터의 분포가 정규분포가 아니면 사용하기 어렵다.
- 구현 코드 -
먼저 필요한 모듈을 임포트한다.
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc
데이터는 앞에서 다룬 이상치 3개를 포함하는 데이터를 활용했다.
data = [9.1, 79.5, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 78.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
data = np.array(data)
다음으로 특정 기준값 $R$을 찾는 함수가 필요하다. 나는 위키에서 제공해준 코드를 그대로 가져왔다. 이 코드는 $R$이 아닌 $R^2$을 찾는 함수이므로 실제로 사용할 때에는 제곱근을 취해준다.
def peirce_dev(N: int, n: int, m: int) -> float:
"""Peirce's criterion
Returns the squared threshold error deviation for outlier identification
using Peirce's criterion based on Gould's methodology.
Arguments:
- int, total number of observations (N)
- int, number of outliers to be removed (n)
- int, number of model unknowns (m)
Returns:
float, squared error threshold (x2)
"""
# Assign floats to input variables:
N = float(N)
n = float(n)
m = float(m)
# Check number of observations:
if N > 1:
# Calculate Q (Nth root of Gould's equation B):
Q = (n ** (n / N) * (N - n) ** ((N - n) / N)) / N
#
# Initialize R values (as floats)
r_new = 1.0
r_old = 0.0 # <- Necessary to prompt while loop
#
# Start iteration to converge on R:
while abs(r_new - r_old) > (N * 2.0e-16):
# Calculate Lamda
# (1/(N-n)th root of Gould's equation A'):
ldiv = r_new ** n
if ldiv == 0:
ldiv = 1.0e-6
Lamda = ((Q ** N) / (ldiv)) ** (1.0 / (N - n))
# Calculate x-squared (Gould's equation C):
x2 = 1.0 + (N - m - n) / n * (1.0 - Lamda ** 2.0)
# If x2 goes negative, return 0:
if x2 < 0:
x2 = 0.0
r_old = r_new
else:
# Use x-squared to update R (Gould's equation D):
r_old = r_new
r_new = np.exp((x2 - 1) / 2.0) * erfc(
np.sqrt(x2) / np.sqrt(2.0)
)
else:
x2 = 0.0
return x2
이제 Peirce's Criterion을 이용하여 이상치를 판단하는 함수를 만들어준다. 데이터와 테이블을 인자로 받고 이상치와 해당 인덱스를 출력한다.
def peirce_criterion(data):
mean = np.mean(data)
std = np.std(data,ddof=1)
N = len(data)
temp_data = data.copy()
outlier_idx = []
n = 1 ## 탐지할 이상치 개수
while n <= N-1:
dev = np.abs((temp_data - mean)/std)
R = np.sqrt(peirce_dev(N,n,1))
idx = np.where(dev>R)[0]
if len(idx) > 0:
temp_data_points = temp_data[idx]
temp_idx = []
for x in temp_data_points:
temp_idx += list(np.where(data == x)[0])
outlier_idx += temp_idx
temp_data = np.delete(temp_data,idx)
n += len(idx)
else:
break
outlier_idx = np.array(outlier_idx)
value = data[outlier_idx]
if len(outlier_idx) > 0:
return (outlier_idx, value)
else:
return False
이제 이상치를 탐지해보자.
res = peirce_criterion(data)
if res:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print('There is no outlier')
위 코드를 실행하면 아래와 같이 3개의 이상치를 한 번에 찾을 수 있다. Chauvenet's Criterion은 한번에 하나의 이상치만을 찾았다는 것을 상기하자. 물론 Chauvenet's Criterion을 반복 적용하면 많이 찾을 수 있긴 하다.
4. Dixon's Q-Test
- 정의 -
Dixon's Q-Test는 정규분포로 추출된 데이터의 이상치를 순서 통계량을 이용하여 검정하는 방법이다.
- 가설 -
Dixon's Q-Test는 다음의 가설을 검정한다.
$H_0$ : 데이터에 이상치가 하나도 없다.
$H_a$ : 데이터에 이상치가 1개 있다.
- 검정 방법 -
데이터 $y_1, y_2, \ldots, y_n$ 이 주어졌다고 했을 때 Dixon's Q-Test는 다음과 같이 진행한다. 검정은 단측 검정과 양측 검정으로 나뉠 수 있다. 단측 검정은 최소값 또는 최대값 하나에 대해서만 이상치인지 아닌지 테스트를 수행하는 것이며 양측 검정은 최소값, 최대값 모두 이상치인지 아닌지를 검정하는 것이다.
1 단계) 데이터를 오름차순으로 정렬한다. 정렬된 데이터를 $y_{(i)}, i=1, \ldots, n$라 하자. 즉,
$$y_{(1)} \leq y_{(2)} \leq \cdots\leq y_{(n-1)}\leq y_{(n)}$$
이다.
2 단계) 검정 통계량을 구한다.
$$r = \frac{y_{(2)}-y_{(1)}}{y_{(n)}-y_{(1)}} \text{ or } \frac{y_{(n)}-y_{(n-1)}}{y_{(n)}-y_{(1)}}$$
3 단계) 유의 수준 $\alpha$와 데이터 개수 $n$에 대하여 단측 검정이라면 기각값 $c_{\alpha, n}$을 구한다. 기각값은 테이블을 이용하여 찾을 수 있다(난 기각값을 계산하는 과정을 구현하였다).
4 단계) 검정 통계량 $r$과 $c_{\alpha, n}$을 비교한다. $r > c_{\alpha, n}$이면 귀무가설 $H_0$를 기각한다. 즉, $y_{(1)}$, $y_{(n)}$을 이상치라고 판단한다.
예를 들어 데이터가 1, 2, 3, 4, 5가 있다고 하자. 그렇다면 $r = \frac{2-1}{5-1} = 0.25$다. 유의 수준 5%($\alpha=0.05$)에서 양측 검정을 수행하기 위한 기각값을 구해야 한다. 데이터의 개수가 5개 이므로 아래와 같이 $n=5$인 열, $\alpha=0.05$=($Q_{95\text{%}}$)인 행에 해당하는 값 0.710이 기각값이 된다(아래 표는 위키에서 가져온 것이다).
이때 0.25 < 0.710이므로 귀무가설을 기각할 수 없고 1은 이상치가 아니라고 판단할 수 있다. 마찬가지로 최대값인 5 또한 이상치가 아니라고 판단할 수 있다.
- 장단점 -
1) 장점
검정 통계량 계산이 간단하다. 필요로 데이터 개수가 적을 때(3~10) 사용 가능하다.
2) 단점
Dixon's Q-Test 역시 데이터의 분포가 정규분포가 아니면 유효하지 않다.
만약 인접 데이터($y_{(1)}, y_{(2)}$ 또는 $y_{(n)}, y_{(n-1)}$)가 이상치임이 확실한 상황에서 이 둘의 차이가 작다면 귀무가설을 기각하지 못할 수 있다. 이를 보완하기 위해 다음과 같은 통계량을 제안하기도 했다.
$$r_{j,i-1} = \frac{y_{(n)}-y_{(n-j)}}{y_{(n)}-y_{(i)}} \;\;\text{or } \frac{y_{(1+j)}-y_{(1)}}{y_{(n-i+1)}-y_{(1)}}$$
Rorabacher는 $j=1,2, \; i=1,2$인 경우의 검정 통계량을 사용하는 것을 추천했다. 또한 Rorabacher는 $r_{10}$은 데이터의 개수가 3개 이상 7개 이하일 때, $r_{11}$은 8개 이상 10개 이하, $r_{21}$은 11개 이상 13개 이하, $r_{22}$는 14개 이상일 때 사용하는 것을 권장했다. 하지만 그렇다 해도 각 통계량에 대한 기각값 테이블을 따로 봐야하는 불편함이 있다(참고로 위에 첨부한 테이블은 $r_{10}$에 대한 기각값을 나타낸 것이다). 그리고 보완했다고 하지만 역시 이상치가 여러개 있는 경우 $H_0$를 기각하지 못할 가능성이 있다. 참고로 위에서 살펴보았던 검정통계량 $r$은 사실 $r_{10}$이며 위키에 따르면 $r_{10}$ 이외의 다른 값은 잘 사용하지 않는다고 한다.
기각값 $c_{\alpha, n}$의 계산이 상당히 어렵다. 따라서 웬만하면 테이블을 사용하자.
- 구현 코드 -
먼저 필요한 모듈을 임포트 하자.
import numpy as np
import matplotlib.pyplot as plt
from math import fabs
from collections import deque
from scipy.stats import norm
from scipy.special import erf
from scipy.optimize import root
여기서 사용할 데이터를 만들어 준다.
data = [9.1, 79.5, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 78.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
data = np.array(data)
데이터의 분포를 살펴보자(사실 앞에서 본 데이터와 같다).
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
이상치는 3개인 것이 확실해 보이지만 이 3개의 차이값은 작으므로 기각이 안될 것이 눈에 선하다.
필요한 함수를 정의한다. 먼저 이번 주제와는 상관없지만 곱셉을 수행하는 함수를 정의했다. 이 함수는 배열에 있는 값들을 모두 곱한다. 이렇게 해준 이유는 아주 작은 값을 곱할 때도 꽤나 정확하게 곱해준다고 해서 가져왔다(사실 그냥 곱셈 연산자 * 와 큰 차이는 느끼지 못했다).
def prod(seq):
s = deque()
for x in seq:
s.appendleft(x) if fabs(x) < 1.0 else s.append(x)
while len(s) > 1:
x = s.popleft() * s.pop()
s.appendleft(x) if fabs(x) < 1.0 else s.append(x)
return s[0] if s else 1.0
그리고 아래의 함수는 semi-infinite Gauss quadrature의 가중치와 특정 지점들을 구하기 위한 함수이다. 이 함수에 대한 설명은 여기를 참고 하자(참고로 해당 링크에는 prod 함수를 쓰지 않고 직접 곱셈을 했다).
def get_exp_xpower_int(k, b=None):
'''
if b == None, b == infty
'''
if b:
if k==0:
return (np.sqrt(np.pi)/2)*erf(b)
elif k==1:
return 0.5*(1-np.exp(-np.square(b)))
else:
# return ((k-1)*0.5)*get_exp_xpower_int(k-2,b)-0.5*np.power(b,k-1)*np.exp(-np.square(b))
return ((k-1)*0.5)*get_exp_xpower_int(k-2,b)-0.5*np.exp(-np.square(b)+(k-1)*np.log(b))
else:
if k==0:
return np.sqrt(np.pi)/2
elif k==1:
return 1/2
else:
return ((k-1)*0.5)*get_exp_xpower_int(k-2)
def get_exp_poly_int(coef,b=None):
res = 0
for i in range(len(coef)):
# res += coef[i]*get_exp_xpower_int(i,b)
res += prod([coef[i],get_exp_xpower_int(i,b)])
return res
def get_polynomial(n,b=None):
if n == 0:
return np.array([1])
elif n == 1:
if b:
return np.array([1, -(1-np.exp(-np.square(b)))/(np.sqrt(np.pi)*erf(b))])
else:
return np.array([1, -(1/np.sqrt(np.pi))])
else:
first = get_polynomial(n-1,b)
## calculate gamma
p = np.poly1d(first)
p2 = p*p ## square polynomial
p2_coef = np.flip(p2.coef) ## reverse coefficient
gamma_first = get_exp_poly_int(p2_coef,b) ## gamma[n-1]
p2_coef_shifted = np.insert(p2_coef,0,0) ## x*polynomial
integral_component = get_exp_poly_int(p2_coef_shifted,b)
# alpha = -((1/gamma_first))*integral_component
alpha = -prod([1/gamma_first,integral_component])
first_plus_order = np.insert(first, len(first), 0) ## x*p(x)
first_plus_order = np.poly1d(first_plus_order)
first_multiplied_const = alpha*first ## alpha*p(x)
first_multiplied_const = np.poly1d(first_multiplied_const)
second = get_polynomial(n-2,b) ## second term
ps = np.poly1d(second)
ps2 = ps*ps
ps2_coef = np.flip(ps2.coef)
gamma_second = get_exp_poly_int(ps2_coef,b)
# beta = -gamma_first/gamma_second
beta = -prod([gamma_first,1/gamma_second])
second_multiplied_const = beta*second
second_multiplied_const = np.poly1d(second_multiplied_const)
return (first_plus_order+first_multiplied_const+second_multiplied_const).coef
def get_weight_and_poins(n,b=None):
assert n > 0 and isinstance(n,int), 'n must be integer greater than 0'
coef = get_polynomial(n,b)
coef_prev = get_polynomial(n-1,b)
points = np.roots(coef)
p = np.poly1d(coef)
p_prev = np.poly1d(coef_prev)
p2_prev = p_prev*p_prev
p2_prev_coef = np.flip(p2_prev.coef)
gamma = get_exp_poly_int(p2_prev_coef,b)
first = p.deriv()(points)
second = p_prev(points)
denominator = first * second
weights = gamma/denominator
weights
return (weights, points)
다음으로 Dixon's Q-Test의 기각값을 구하는 함수를 만든다. 기각값을 구할때 시간이 꽤 걸린다(대략 25~30초). $i$와 $j$값을 조정하여 유의 수준 $\alpha$에서의 $r_{j, i-1}$에 대한 기각값을 구할 수 있다.
def get_cv_Dixon_test(n,i=1,j=1,alpha=0.05):
fac = np.math.factorial
cdf = norm.cdf
N = ((2*np.pi)**(-1.5))*(fac(n)/(fac(i-1)*fac(n-j-i-1)*fac(j-1)))
fh_res = np.polynomial.hermite.hermgauss(30)
fh_weights = fh_res[1]
fh_points = fh_res[0]
hh_res = get_weight_and_poins(15)
hh_weights = hh_res[0]
hh_points = hh_res[1]
gl_res = np.polynomial.legendre.leggauss(12)
gl_weights = gl_res[1]
gl_points = gl_res[0]
def J(x,v,r):
first_elem = np.power(cdf(x-v),i-1)
second_elem = np.power((cdf(x-r*v)-cdf(x-v)),n-j-i-1)
third_elem = np.power((cdf(x)-cdf(x-r*v)),j-1)
return prod([first_elem,second_elem,third_elem])
def second_func(t,u,r):
x = u*np.sqrt(2/3)
v = t*np.sqrt(2/(1+r**2))
last_term = np.exp(2*u*t*(1+r)/np.sqrt(3*(1+r**2)))
return prod([J(x,v,r),last_term,t])
def p(r):
const = N*np.sqrt(2/3)*(2/(1+r**2))
summand = 0
for k, u in enumerate(fh_points):
for l, t in enumerate(hh_points):
summand += prod([fh_weights[k],hh_weights[l],second_func(t,u,r)])
return const*summand
def g(R):
# return quad(p,0,R)[0]
## or
summand = 0
for w, y in zip(gl_weights,gl_points):
summand += prod([w,p(prod([0.5,R,y+1]))])
return prod([0.5,R,summand])
def objective(R,alpha):
return 1-alpha-g(R)
## cv = brentq(objective,a,b,args=(alpha))
## Brent method는 초기 구간을 잡는 것이 어렵다.
cv = root(objective,0.5,args=(alpha))
return cv.x[0]
다음으로 Dixon's Q-Test를 수행하는 함수를 정의한다. 아래 함수는 주어진 데이터(data)에 대하여 유의수준 $\alpha$(alpha)에서 검정 통계량 $r_{j,i-1}$을 이용한 Dixon's Q-Test를 수행한다. two_sided 인자는 True인 경우 양측 검정을 False인 경우 단측 검정을 수행한다. 단측 검정인 경우, lowest=True라면 최소값이 이상치인지 테스트하고 False 라면 최대값에 대해서 테스트한다.
def dixon_q_test(data,alpha=0.05,i=1,j=1,two_sided=True,lowest=True):
temp_data = data.copy()
temp_data = np.sort(temp_data)
test_statistics1 = (temp_data[j] - temp_data[0])/(temp_data[-i]-temp_data[0])
test_statistics2 = (temp_data[-1] - temp_data[-1-j])/(temp_data[-1]-temp_data[i-1])
n=len(data)
if two_sided:
cv = get_cv_Dixon_test(n,i=i,j=j,alpha=alpha*0.5)
outlier_idx = []
if test_statistics1 > cv:
idx = np.where(data == temp_data[0])[0][0]
outlier_idx.append(idx)
if test_statistics2 > cv:
idx = np.where(data == temp_data[-1])[0][0]
outlier_idx.append(idx)
if outlier_idx:
outlier_idx = np.array(outlier_idx)
values = data[outlier_idx]
return (outlier_idx, values, test_statistics1, test_statistics1, cv)
else:
return (test_statistics1, test_statistics1, cv)
else:
cv = get_cv_Dixon_test(n,i=i,j=j,alpha=alpha)
if lowest: ## test whether minimum of data is outlier or not
if test_statistics1 > cv:
idx = np.where(data == temp_data[0])[0][0]
return (idx, data[idx], test_statistics1, cv)
else:
return (test_statistics1, cv)
else: ## test whether maxmum of data is outlier or not
if test_statistics2 > cv:
idx = np.where(data == temp_data[-1])[0][0]
return (idx, data[idx], test_statistics2, cv)
else:
return (test_statistics2, cv)
이제 Dixon's Q-Test를 해보자. 시간이 좀 걸린다. 나 같은 경우 20초 걸렸다 ㄷㄷ;;
res = dixon_q_test(data,alpha=0.05,i=1,j=1,two_sided=True,lowest=True)
if len(res)==5:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print('There is no outlier')
역시 예상했던 대로 이상치가 발견되지 않았다. 이번엔 이상치가 하나 있는 데이터로 바꿔서 해본다. 데이터만 아래 걸로 바꾸고 위 코드를 실행하자.
data = [199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57]
data = np.array(data)
이상치를 잘 탐지해냈다.
사실 Dixon's Q-Test에서 기각값을 계산하는데 많이 애를 먹었다. Guass Quadrature, semi-infinite Gauss Quadrature, Gauss-Legendre 적분 근사와 해를 찾기위한 수치적 방법들이 필요하기 때문이다. R에서는 Dixon 검정의 기각값을 계산해주는 라이브러리가 있지만 파이썬에서는 없는 것 같다. 그렇다고 내가 만든 코드는 라이브러리로 만드는 것까지는 무리인 듯하다. 시간이 많이 걸리기 때문이다. 그래도 파이썬으로 기각값을 계산하는 과정을 완성했다는 것에 의미가 있다고 생각하며 나름 뿌듯하다 ㅎㅎ.
5. Generalized Extreme Studentized Deviation Test
- 정의 -
Generalized Extreme Studentized Deviation(G-ESD) Test는 Grubbs's Test가 이상치를 한 번에 하나밖에 찾지 못하는 단점을 보완한 검정방법이다. G-ESD Test는 사전에 이상치의 개수를 $r$로 정해줌으로써 데이터에 여러 개의 이상치가 있는지 통계적으로 검정해볼 수 있다.
- 가설 -
G-ESD Test에서 검정하고자 하는 가설은 다음과 같다.
$H_0$ : 데이터에 이상치가 하나도 없다.
$H_a$ : 데이터에 이상치가 최대 $r$개 있다.
- 검정 방법 -
G-ESD Test는 총 $r$번의 반복 과정이 필요하다. 먼저 데이터 $y_1, y_2, \ldots, y_n$이 주어졌다고 하자. 유의 수준 $\alpha$에서 G-ESD Test는 다음과 같이 수행한다.
1 단계)
$i$ 번째 스텝에서 $R_i$와 $\lambda_i$를 다음과 같이 계산한다($i=1, 2, \ldots, r$).
$$R_i = \max_j|x_j-\bar{x}|/s, \;\; \lambda_i = \frac{(n-i)t_{n-i-1, \alpha/2(n-i+1)}}{\sqrt{(n-i-1+t_{n-i-1,\alpha/2(n-i+1)}^2)(n-i+1)}}$$
이때 $\bar{x}$와 $s$는 $i=1$인 경우 원 데이터에 대한 표본 평균과 표본표준편차이며 $i>1$인 경우 원 데이터에서 R_{i-1}에 대응하는 데이터를 제거한 후 새로 계산한 표본평균과 표본 표준편차이다. 그리고 $t_{n,p}$는 자유도가 $n$인 $t$분포에서 오른쪽 꼬리 확률이 $p$가 되게 하는 값이다.
2 단계)
$m = \max_i\{i : R_i>\lambda_i\}$라 할 때 데이터에는 총 $m$개의 이상치가 있다고 판단하며 $m$이 존재하지 않는 경우에는 이상치가 하나도 없다고 판단하게 된다.
만약 $r=1$인 경우 Grubbs's Test의 양측 검정과 같다.
- 장단점 -
1) 장점
Grubbs's Test에서는 한 번의 검정으로 하나의 이상치만을 탐지해낼 수 있었던 것과는 달리 G-ESD Test는 한 번의 검정으로 여러 개 이상치를 탐지할 수 있다. 비교적 구현이 간단하다.
2) 단점
데이터의 분포가 정규분포와 다른 경우에는 G-ESD Test의 결과가 유효하지 않다. 또한 사전에 이상치의 개수를 미리 지정해야 한다는 점도 단점으로 꼽힌다. 마지막으로 G-ESD Test는 양측 검정 용도로만 제한된다고 한다(즉, 단측 검정은 안된다).
- 구현 코드 -
먼저 필요한 모듈을 임포트한다.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t
다음으로 앞에서 계속 사용했던 이상치가 3개 있는 데이터를 만든다.
data = [9.1, 79.5, 26.8, 81.5, 19.1, 15.2, 22.6, 28.8, 24.1, 23.6,
18.6, 17.3, 25.8, 78.8, 23.1, 11.9, 20.1, 20.3, 14.1, 26.5]
data = np.array(data)
이제 G-ESD Test를 수행하는 함수를 만든다. 이 함수는 검정을 수행할 데이터(data), 이상치 최대 개수 $r$(num_outliers) 그리고 유의 수준(alpha)을 인자로 받으며 이상치 인덱스, 이상치 값, $R_i$, $\lambda_i$를 리턴한다.
def generalized_esd(data,num_outliers,alpha=0.05):
assert len(data)-num_outliers > 0, 'invalid num_outliers'
n = len(data)
temp_data = data.copy()
res = []
for i in range(num_outliers):
mean = np.mean(temp_data)
std = np.std(temp_data, ddof=1)
diff = np.abs(temp_data-mean)
R = np.max(diff)/std
t_val = t.ppf(1-alpha/(2*(n-i)),n-i-2)
lambda_val = (n-i-1)*t_val/np.sqrt((n-i-2+t_val**2)*(n-i))
temp_idx = np.where(diff==np.max(diff))[0][0]
temp_data_point = temp_data[temp_idx]
idx = np.where(data == temp_data_point)[0][0] ## index of suspected outlier
value = data[idx] ## suspected outlier
flag = R > lambda_val
res.append((idx,value,flag,R,lambda_val))
temp_data = np.delete(temp_data,temp_idx)
if res:
idx_suspected_outlier = []
for i, r in enumerate(res):
if r[2] == True:
idx_suspected_outlier.append(i)
num_suspected_outlier = max(idx_suspected_outlier)+1
outlier_idx = [res[i][0] for i in range(num_suspected_outlier)]
outlier_idx = np.array(outlier_idx)
values = data[outlier_idx]
Rs = [res[i][3] for i in range(num_suspected_outlier)]
lambdas = [res[i][4] for i in range(num_suspected_outlier)]
return (outlier_idx, values, Rs, lambdas)
else:
return False ## no outlier detected
이제 데이터에 이상치가 있는지 검정해보자(Grubbs's Test는 이상치를 하나도 잡아내지 못했음을 기억하자). 먼저 3개가 있음을 시각적으로 확인했으니 num_outliers = 3으로 설정했다.
res = generalized_esd(data,num_outliers=3,alpha=0.05)
if len(res) == 4:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print(res)
print('There is no outlier')
이제 정말 G-ESD Test가 좋은지 확인해봐야 한다. num_outliers = 10으로 설정했을 때에도 정확히 3개만 잡아내는지 확인해본다.
res = generalized_esd(data,num_outliers=10,alpha=0.05)
if len(res) == 4:
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.scatter(range(len(data)), data)
plt.scatter(res[0],res[1],s=200,facecolor='none',edgecolors='r')
plt.xticks(range(len(data)), range(1,len(data)+1))
plt.show()
else:
print(res)
print('There is no outlier')
정확하게 3개 만을 잡아내었다.
참고자료
Grubbs's Test - en.wikipedia.org/wiki/Grubbs%27s_test
Grubbs's Test for Outliers - www.itl.nist.gov/div898/handbook/eda/section3/eda35h1.htm
Lin et al - Cleaning Data the Chauvenet Way
Chuavenet's Criterion: Definition & Example - www.statology.org/chauvenets-criterion/
Peirce's Criterion - en.wikipedia.org/wiki/Peirce%27s_criterion
Dixon's Q Test - en.wikipedia.org/wiki/Dixon%27s_Q_test
Statistical Treatment for Rejection of Deviant Values: Critical Values of Dixon’s “Q” Parameter and Related Subrange Ratios at the 95% Confidence Level
Extreme Studentized Deviate Test - www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/esd.htm
댓글