통계/머신러닝 기반의 시계열 이상치(Anomaly) 탐지 조사

  이상치 탐지 솔루션은 공정, 금융 등의 분야에서 비교적 흔하게 사용되고 있는 기술이다. 예컨대, 공정에서 공기청정기를 생산한다고 할 때, 이상치 탐지 기술을 이용하면 (그 안에 들어가는) 모터에서 발생하는 소음의 "이상" 여부를 분석하여 해당 부품의 불량 여부를 체크할 수 있고, 금융권의 경우, 은행에서 특정 사용자의 계좌에서 평소와 다른 "이상" 패턴이 탐지된다면 그 사용자의 계좌에 대해 동결 조치를 취할 수 있다.

  이러한 이상치 탐지 기술은 통계적인 접근과 머신러닝 접근으로 구분될 수 있는데, 전자는 오래 전부터 사용된 고전적인 방법인 반면 후자는 비교적 최근 들어 더 인기를 얻게 되었다. 이번 포스팅에서는 이 둘에 대한 대표적인 방법을 소개한다. 통계적 접근은 개발이 쉽고 결과 해석이 용이하다는 이점에 있고, 머신러닝 접근은 탐지 성능 측면에서 더 우수할 수 있다는 이점이 있어 보인다.

1. 통계 기반 이상치 탐지 방법론

1.1 Shewhart Control Charts

  통계적인 접근의 이상치 탐지 방법은 Walter A. Shewhart (1891 - 1967)에 의해 처음으로 제안되었다. 이 방법은 관리도(Control Chart)라고 불리는데, 이것은 시간에 따라 변화하는 데이터의 변동을 추적 관찰(Monitoring)하면서 해당 데이터가 특정 시점에 이상이 발생하는지 여부를 모니터링하는 그래프이다. 이 관리도는  관리 통계량(Control statistics, points)과 관리 한계선(Control limits, lines)으로 구성된다. 여기서 관리 통계량이란 관리도가 추적 관찰하는 대상의 통계량 (평균, 표준편차 등)을 일컫는다. 이것은 이상 여부의 임계치 역할을 하는 관리 한계선을 벗어나는지를 추적 관찰하기 위한 대상이다.

그림1 Control Chart

  관리도에서 관리 한계선을 찾을 때에는 정상 데이터만으로 구성된 데이터세트의 통계량으로부터 계산되어야 한다. 더욱이, 이 과정은 Phase1에서만 진행되어야 한다. 예를 들어, 그림2와 같이 15 시점을 기준으로 과거 기간인 phase1 동안에 관리 한계선이 생성되면, 15 이후의 기간인 phase2에 이 생성된 관리 한계선이 이상 여부를 탐지하기 위해 사용된다.

[그림2] phase1,2 구분

1.1.1 X bar Chart

  먼저, 대표적인 관리도 방법인 X bar Chart에 대해 알아보자. 이것은 연속적인 데이터의 이상 여부를 판별할 때 적합한 방법론이다. 이것은 변수가 하나인 단변량 데이터를 다루는 데 적합하다. 여기서 X bar ($\overline{x}$)가 의미하는 것은 하나의 샘플에 대한 변량의 평균값이다. 아래 [그림3]을 참고할 때, 하나의 서브그룹에는 5번의 반복된 관찰을 통해 얻은 5 개의 변량이 있다. 여기서 예컨대 서브그룹1에 대한 $\overline{x}$는 1행에 있는  5 개 값의 평균인 102.9이다. 이러한 $\overline{x}$를 기반으로 관리 한계선을 구하는 방법을 X bar 차트라고 일컫는다.

[그림3] example

 위 [그림3]의 $\overline{\overline{x}}$는 이 테이블 내 데이터를 모두 요약하고 있는 각 평균값에 대한 평균이다. 이 값은 관리 한계선의 중심선(Center line)으로 간주된다. 하지만 여기서 아래 [그림4]와 같이 Upper Control Limit(UCL)과 Lower Conter Limit(LCL)을 구하는 것이 중요하다. 이들을 통해 변동하는 시계열 데이터가 특정 시점에 이상인지 여부를 판정할 수 있기 때문이다.

[그림4] x bar chart 예

X bar chart에서 관리 한계선은 다음 [그림5] 같이 구한다. 여기서 $\overline{R}$는 각 샘플의 (최댓값 - 최솟값)의 평균을 의미한다. $A_2$는 특정한 상수이다. 여기서 $A_2\overline{R}$가 의미하는 것은 중심선인 $\overline{\overline{x}}$로부터 허용 가능한 정상 범위이다.

[그림5] X bar 차트의 관리 한계선 산식

  수학적인 접근으로 UCL에 대해 더 자세히 이해해보자. $\overline{x}$ chart의 UCL을 이해하기 위해, 식 (1) $UCL=E(\overline{x})+3sd(\overline{x})$ 에서 시작해보자. 여기서 3은 경험적으로 얻은 계수이고, sd는 표준편차(standard deviation)이다. (1)에서 $E(\overline{x})$는 $\overline{\overline{x}}$로 근사되고 후항은 (2) $sd(\overline{x})=\frac{\sigma }{\sqrt{n}}$로 달리 표현될 수 있다. 여기서 관건은 우리가 (2)에서 $\sigma$를 알지 못한다는 점이다. 그래서 이것을 우리가 구해내야 할 estimator로서 간주한다는 의미에서 $\hat{\sigma }$으로 대체해보자. 그러면 (1)은 (3) $\overline{\overline{x}}+\frac{\hat{\sigma }}{\sqrt{n}}$으로 표현될 수 있다.

  이제, 이 (3)이 (4) $UCL=\overline{\overline{x}}+A_2\overline{R}$이 되는 과정에 대해 알아보자. 먼저, Relative range (5) $w=\frac{R}{\sigma }$를  정의해보자. 여기서 $R$은 샘플의 max value - min value이다. 이 $w$의 기댓값은 $E(w)=d_2$이고 표준편차는  $sd(w)=d_3$이다. 이 두 값은 이미 계산되어 있는 상수다. 예컨대, 샘플의 개수가 20개일 때 $d_2$는 보통 5 내외이고 $d_3$는 0.4 내외이다. 우리가 구하고자 하는 것은 estimator인데, 이것은 (6) $\hat{\sigma }=\frac{\overline{R}}{d_2}$이다. 즉, $\frac{\overline{R}}{d_2}$이 우리가 알고 있지 못했던 $\sigma$를 추정하는 좋은 estimator라는 것이다. 이 말은 즉 $\hat{\sigma }$이 unbiased estimator라는 것을 뜻한다.

하지만 unbiased estimator를 만족하기 위한 조건이 있다. 이것은 "만약 $E(\hat{\sigma })=\sigma$ 이면,$\hat{\sigma }$은 unbiased estimator다"에 의해 성립한다. 그러면, $E(\hat{\sigma })$가 과연 이 조건을 만족하는지 알아보자면, (7) $E(\hat{\sigma })=E(\frac{\overline{R}}{d_2})=\frac{1}{d_2}E(\overline{R})$인데, $R$에 대한 평균을 구하기 위해 (5)를 변형하면 (8) $R=w\sigma$가 된다. 이 (8)에 Expectation을 취하면 $E(R)=\sigma \ast E(w)$이고, 이것은 전술한 $E(w)=d_2$에 의해 $\sigma \ast E(w)=\sigma \ast d_2$가 된다. 이제, $R=\overline{R}$이라고 가정할 때, (7)에서 $E(\hat{\sigma })=\frac{1}{d_2}E(\overline{R})$에 (8)을 대입하면 $E(\hat{\sigma })=\frac{1}{d_2}\sigma d_2$가 되는데, 그러므로 이는 $E(\overline{\hat{\sigma }})=\sigma$가 됨으로써 $\hat{\sigma }$이 unbiased estimator라는 것이 증명되는 것이다.

이제 (6)을 (1)에 대입하면, (9) $\overline{\overline{x}}+3\frac{\frac{\overline{R}}{d_2}}{\sqrt{n}}$이 되는데, 그러면 우리는 이 식에서 n은 반복수 ([그림3]의 테이블을 참고할 때 여기서의 반복수는 5이다.)이고, $d_2$는 샘플 개수가 20개인 경우 5 내외이고, $\overline{R}$은 샘플들의 max-min 값의 평균일 뿐이므로 모두 구할 수 있는 값이다. 이것은 [그림5] 식에서의 $A_2$는 (9)를 참고할 때 $A_2=\frac{3}{d_2\sqrt{n}}$인 것이다. 이 $A_2$값 역시 $d_2,d_3$처럼 이미 계산된 값이고 샘플 수가 20개인 경우 이 값은 대략 0.18이다. LCL 역시도 동일한 방법으로 구할 수 있다.

1.2 Time-Weighted Control Charts

 전술한 Shewhart control 방법은 사실 오직 최근 정보만을 고려하여 관리 한계선을 생성한다는 한계가 있다. 이 경우, 모니터링 시스템은 관리 통계량의 작업 변화를 포착하지 못할 수 있다. 예컨대, Shewhart control 방법에 따라 개발된 시스템은 t 시점까지의 평균이 20이고 그 이후부터 평균이 22인데 t 이후의 모든 시점에서의 통계량이 UCL, LCL을 넘지 않는다면 이상anomaly에 대한 알람이 울리지 않을 것이다. 이러한 한계를 극복하기 위해, 이번 섹션에서는 과거의 정보도 고려하여 관리 한계선을 생성하는 대표적인 방법인 Exponentially Weighted Moving Average (EWMA)에 대해 알아본다.

1.2.1 EMWA 

  먼저, EWMA statistics는 (10) $z_i=\lambda x_i+(1-\lambda )z_{i-1}$이다. 즉, 관리도에서 하나의 타점은 (10)에 의해 계산된다. 여기서 $i$는 특정 시점을 일컫고 $\lambda$는 0과 1 사이의 가중치 역할을 하는 실수인 hyperparameter이고, $x_i$는 $i$ 시점에서의 입력값이다. 식 (10)을 자세히 보면 전항은 현재 시점의 정보이고 후항은 과거 시점에 가중합된 정보임을 알 수 있다. 참고로, 첫 시점인 $z_0={\mu }_0=\overline{x}$이다. 이제 (10)을 예를 들어 자세히 알아보자.

세 개의 시점만을 고려한다고 가정할 때, 아래 [그림6]같이 $z_3$가 최종적인 관리 통계량이 될 것이다. 이 $z_3$를 자세히 들여다보면 $z_2$가 있고, 그 안에 $z_1$, $z_0$도 있음을 [그림6]의 마지막 유도식을 통해 확인할 수 있다. 이것은 EMWA statistics가 과거의 모든 정보도 고려하고 있음을 함축한다. 또 다른 주목할 점은 $(1-\lambda {)}^n$ 가중치를 관리 통계량에 곱함으로써 현재에서 더 가까이 있는 항일수록 상대적으로 높은 가중치가 부여되고, 더 과거에 있는 항일수록 상대적으로 낮은 가중치가 부여된다는 점이다.

[그림6 EWMA 식]

나아가, 위 식에 따르면, $\lambda$가 작을수록 과거 관측치에 더 높은 가중을 주는 것을 함의한다. 이는 작은 변화를 탐지하는 데 유익하게 해주는 반면, $\lambda$가 클수록 현재 관측치에 더 높은 가중을 줌으로써 큰 변화를 탐지하는 데 더 유익하게 해준다.

이제, EWMA의 Control Limit을 알아보자. UCL의 일반식은 $UCL=E(z)+3sd(z)$이었다. 이것을 알기 위해서는 EWMA statistics로부터 $E[z_i]$와 $SD[z_i]$를 구해야 한다. 지금부터 이 둘을 구해보자.

가정: $x_i\sim N(\mu ,{\sigma }^2)$, for all i

$E[z_i]$를 풀어 쓰면 아래 [그림7]의 (a)와 같고, (a)는 위 정규분포 가정에 의해 (b)로 표현될 수 있다. 여기서 공통 인수를 묶으면 (c)가 되고, (c)는 등비수열의 합공식에 의해 (d)가 된다. 마지막으로, 이 식을 정리하면 (e)가 되는데, 위에서 참고로 이야기했던 $z_0=\mu$이기 때문에, 식(e)가 정리되어 결과적으로 $E[z_i]=\mu$가 된다.

[그림7] E[z_i]=mu

표준편차의 제곱인 분산 $V[z_i]$도 [그림8]과 같이 유사한 방식으로 구할 수 있다. (f)에서 상수 $(1-\lambda {)}^i{z}_0$의 분산은 0이고 s$d(x_i)$에서 제곱을 취한 $V[x_i]$에서 계수는 제곱으로 나오니 모두 ${\lambda }^2$, $(1-\lambda {)}^2$로 나오게 된다. 이것을 정리하면 (g)를 거쳐 (h)가 되고, 동일하게 등비수열 합공식에 따라 (i)가 되고 이것을 정리하면 (j)가 된다. 여기에 루트를 씌운 것이 $SD(z_i)$이다.

[그림8] V[zi]

위에서 구한 $E(z_i)$와 $SD(z_i)$를 이용하여 EWMA의 UCL, LCL의 식을 구하면 아래와 같다. 아래에서 $L$은 위 일반식에서 3으로 정의했던 것이다. 이 값은 사용자에 의해 설정될 수 있는 HYPERPARAMETER이다. 또한, ${\mu }_0$는 전술했듯이 $E(z_i)$와 동일한 값이다. 즉, 우리는 위에서 구한 기댓값과 표준편차를 이용하여 EWMA의 관리 한계선을 구할 수 있다.

[그림9] EWMA의 관리 한계선

 

위 [그림9]의 UCL, LCL에서 볼 수 있는 특이한 점은 특정 시점을 나타내는 인덱스 $i$가 식 내부에 있다는 점이다. 이것은 시점마다 관리 한계선이 변할 수 있음을 함의하고, 나아가 이 사실은 $i$가 커짐에 따라 $(1-\lambda {)}^{2}i$는 점점 0으로 수렴한다는 것을 함축한다. 그 결과 특정 시점 이후부터 UCL, LCL은 정상 상태를 유지하게 된다.

[그림10] EWMA chart 예

마지막으로 강조할 점은Time-Weighted Control Charts에서는 Shewhart Control Charts에서와 같이 Phase1에서 정상 데이터만을 사용하지 않는다는 것이다. 즉, EWMA 방법에서는 여기서는 정상/이상 데이터 모두 사용하여 관리 한계선을 생성한다. 또한, 현업에서는 일반적으로$\lambda$를 0.05에서 0.25 사이의 값을 사용하고 $L=3$ 내외에서 설정된다.

2. 머신러닝 기반 이상치 탐지 방법론

  지금까지 고전적인 통계 기반의 이상치 탐지 방법에 대해 알아보았다. 이번 장에서는 최근 들어 활발히 연구되고 있는 머신러닝 기반의 이상치 탐지 방법에 대해 알아보자. 그 중에서도 딥러닝 접근을 이용한 직관적으로 이해하기 쉬운 구조를 가진 DeepAnT [7]에 대해 간략히 살펴본다. 

  먼저, 머신러닝 기반의 이상치 탐지 태스크에 대해 정의해보자.

• Dataset: 시계열 세트 $T=\{x_1,\dots ,x_t\}$가 있고, 임의의 시점 $k$에 대해 각 $x_k$는 $x_k=\{x_{1,k},\dots ,x_{m,k}]\}$로 이루어져 있다. 여기서 $m$은 변수의 개수이고 $t$는 시점의 개수이다. 만약 $m=1$이면 단변량 시계열 데이터이고 $m>1$이면 다변량 시계열 데이터이다.

• Task: 길이가 $s$인 time window $W_t$가 있을 때, 이 태스크는 $W_t=\{x_{t-s+1},\dots ,x_{t-1},x_t\}$를 모델 $\varphi (\theta ,W_t)$에 입력하여 $t$ 시점의 $x_t$가 anomaly인지 여부를 탐지한다.

  그전에 먼저 이상치 탐지에서 다루고 있는 "이상"(anomaly)에 대해 정의할 필요가 있다. 일반적으로 이상은 아래 [그림11] 같이 4 가지로 분류될 수 있다. 우리가 다루는 종류의 anomaly는 1에 해당하는 Point-wise outlier이다. 

1. 연속적인 시계열 안의 특정한 outlier 포인트가 있다.

2. 특정 기간 내 파형이 다른 기간 내 파형과 다른 Shapelet outlier가 있다.

3. 특정 기간 내 파동이 다른 기간 내 파동보다 파장이 긴 Seasonal outlier가 있다.

4. 특정 시점에 진폭의 변화가 큰 Trend outlier가 있다.

  머신러닝 기반의 시계열 이상치 탐지(ADTS)의 방법 또한 여러가지이다. 아래 [그림12]에서 볼 수 있듯이 이것은 5 가지로 분류될 수 있다. 먼저, 클러스터링 방법이 있다. 이것은 다변량 시계열 데이터의 피처 공간 $S^n$에 포함되는 정상 시계열 데이터 피처의 서브 공간 $N^n$에  feature extractor가 시계열 $X$를 피처 공간으로 매핑시킨 $f(X)$이 $N^n$에 포함되지 않는 경우 이 $X$를 이상치로 판정하는 방법이다. 다음으로, 밀도 추정 방법이 있다. 이것은 $p_{\theta }$인 정상 시계열의 확률 밀도 함수를 추정하는 것을 목표로 한다. 시계열 $X$가 주어질 때, $L(\theta |x)$가 임계치 $\tau$보다 큰 경우, 이 $X$는 이상치로 판정된다. 거리 기반 방법으로 이상치를 탐지할 수도 있다. 이 방법은 기준이 되는 정상 시계열(Reference time series)을 $R$로, 주어진 시계열을 $X$로 간주하고, 거리 함수를 $D()$로 정의할 때, $D(X,R)>\delta$인 경우 해당 $X$는 이상치로 판정한다.

[그림12] Paradigms for ADTS

  신경망 접근을 활용할 때 가장 인기있는 두 가지 방법은 재구축(Reconstruction) 기반 접근과 예측(forecasting) 기반 접근일 것이다. 이들 중 재구축 방법은 주어진 원본 데이터와 오토인코더 등을 통해 재구축된 시계열 간의 차이가 특정 임계치보다 큰 경우 이상으로 판정한다. 이와 유사하게, 예측 기반 방법도 임의의 시점 $i$에 대해 주어진 시계열 $X_i$와 신경망 모델에 의해 예측된 $\widehat{X_i}$ 간의 차이가 특정 임계치보다 큰 경우 이상치로 판정한다.

  일종의 예측 기반 방법론인 DeepAnT [7]를 간략히 살펴보자. 이것은 주어진 과거 시점의 시계열 데이터를 기반으로 미래 시점의 값을 예측하는 Time Series Predictor와 예측값과 실제값의 차이를 계산하여 해당 시점의 실제값(그림의 기댓값X)이 이상인지 여부를 판정하는 Anomaly Detector로 구성된다. 여기서 전자는 1D-CNN과 Max pooling 레이어를 이용하여 시계열을 학습한다. 이 예측기는 MAE를 목적 함수로 학습된다.

[그림13] DeepAnT 구조

 

학습된 DeepAnT를 이용하여 시계열의 이상치를 검출해보면 아래 [그림14] 같이 Predictor의 예측값과 실제값 간의 차이가 비교적 큰 경우 Anomaly Detector에서 예측값 $y_t^{\prime }$과 실제값 $y_t$ 간의 거리 차가 클 것이고, 이는 이상치로 판정될 가능성이 클 것임을 알 수 있다.

[그림14] DeepAnT를 이용한 이상치 탐지 예

 

 

참고 문헌

[1] 산업공학과 김성범 교수님의 수업 참고

[2] 그림1 https://cdn1.byjus.com/wp-content/uploads/2018/11/maths/2016/06/03072323/Control-Charts.jpg

[3] 그림2 https://www.researchgate.net/publication/346641575/figure/fig3/AS:965150186143746@1607121144354/Phase-I-and-II-control-chart.png

[4] 그림3 https://www.smartersolutions.com/wp-content/uploads/2017/12/Data-set-sample-of-5-subgrouping-jpeg.jpg

[5] 그림 4 https://www.spcforexcel.com/files/images/xbarbowlgif.gif

[6] 그림 10 https://analyse-it.com/docs/user-guide/process-control/images/ewmachart.png

[7] Munir, Mohsin, et al. "DeepAnT: A deep learning approach for unsupervised anomaly detection in time series." Ieee Access 7 (2018): 1991-2005.

[8] Mejri, Nesryne, et al. "Unsupervised Anomaly Detection in Time-series: An Extensive Evaluation and Analysis of State-of-the-art Methods." arXiv preprint arXiv:2212.03637 (2022).