[ML] 마할로비스(Mahalnobis) 거리와 이상치 감지

마할로비스(Mahalnobis) 거리

마할로비스(Mahalnobis) 거리

벡터인 한점과 분포 사이의 거리를 측정하는 효과적인 다변량 거리 메트릭스

이 지표는 다변량의 이상 탐지, 높은 불균형한 데이터 세트의 분류, 단일 클래스 분류, 새로운 데이터들의 예측에 효과적입니다.2차원의 두 점사이의 거리를 나타낼 경우 일반적으로 유클리드거리가 사용됩니다. 그 두점 (p1, p2), (q1, q2)라고 하면 유클리드 거리는 다음과 같이 계산됩니다.

$$\tag{식 1}d(p, q) =\sqrt{(p_1-q_1)^2+(p_1-q_2)^2}$$

식 1을 다차원으로 확장하면 즉, $(p_1, p_2, \cdots, p_n), \;(q_1, q_2, \cdots, q_n)$ 유클리드 거리는 식 2와 같이 계산됩니다.

$$\tag{식 2}d(p, q) =\sqrt{(p_1-q_1)^2+(p_1-q_2)^2+ \cdots +(p_n-q_n)^2}$$

식 2와 같이 다차원의 경우 모든 차원에 대한 가중치는 전혀 고려되지 않습니다. 즉, 차원들 사이에 영향이 없다는 가정이 성립되어야 합니다(모든 차원은 독립적).

import numpy as np
import numpy.linalg as la
import pandas as pd
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style("darkgrid")

np.random.seed(3)
X=np.random.normal(0, 1, 100)
X1=np.linspace(-2,2, 100)
y1=X1+np.random.randn(100)

p1=(-1, 1)
p2=(1, 1)

fig, ax=plt.subplots(1, 2, figsize=(7, 3), sharey=True)
ax[0].scatter(X1, X, s=5)
ax[0].scatter(p1[0], p1[1])
ax[0].scatter(p2[0], p2[1])
ax[0].plot([0,p1[0]], [0, p1[1]], color="orange")
ax[0].plot([0,p2[0]], [0, p2[1]], color="g")
ax[0].set_xlabel("x")
ax[0].set_ylabel("y")
ax[0].set_title("(a)")
ax[1].scatter(X1, y1, s=5)
ax[1].scatter(p1[0], p1[1])
ax[1].scatter(p2[0], p2[1])
ax[1].plot([0,p1[0]], [0, p1[1]], color="orange")
ax[1].plot([0,p2[0]], [0, p2[1]], color="g")
ax[1].set_xlabel("x")
ax[1].set_title("(b)")
plt.show()

위 그림의 (a)는 x와 y의 데이터 사이에 상관성이 없다고 할 수 있지만 (b)의 경우는 상당한 상관성을 추정할 수 있습니다. 이 관계에서 각 그래프의 중점에서 주황색과 녹색의 점까지의 유클리드 거리 같습니다. 그러나 (b)의 경우는 녹색점이 데이터들에 근접해 있음을 나타냅니다. 즉, 유클리드 거리는 분포에 대한 고려없이 산술적인 거리를 계산하므로 분류등에 적용에 적합하지 않습니다. 이 경우 유클리드 거리를 대신해 점과 분포의 거리인 마할로비스 거리를 사용할 수 있습니다. 마할로비스 거리는 다음 단계로 계산됩니다.

1. 열들을 상관성이 없는 변수로 변환합니다.
2. 각 열을 분산 1로 조정합니다.
3. 유클리드 거리를 계산합니다.

마할로비스 거리를 위한 공식은 식 3과 같습니다.

\begin{align} D^2&=(x-m)^T\Sigma^{-1}(x-m)\\ \tag{식 3}D: &\;\text{마힐로비스 거리}\\ x: &\;\text{관찰치의 벡터(데이터 셋의 행)}\\ m: &\;\text{설명변수의 평균벡터(각 열의 평균)}\\ \sigma: &\;\text{공분산 행렬}\end{align}

식 3의 (x-m)은 각 변수의 관찰치와 평균의 거리입니다. 이 결과는 공분산 행렬로 나누어지므로 본질적으로 다변량에 대한 표준화(식 4)와 같습니다.

\begin{align}\tag{식 4}z=\frac{x-\mu}{\sigma}&\quad \text{for univariate}\\ z=\frac{\vec{u}-\overrightarrow{\text{mean}}}{\overrightarrow{\text{covariance}}}&\quad \text{for multivariate}\end{align}

• 변수들 사이의 강한 상관성 ⇒ 높은공분산 ⇒ 거리는 감소됨
• 변수들 사이의 약한 상관성 ⇒ 낮은공분산 ⇒ 거리는 감소효과는 작음

다음 자료에 대해 마할로비스 거리를 계산해 봅니다.

filepath = 'https://raw.githubusercontent.com/selva86/datasets/master/diamonds.csv'
df = pd.read_csv(filepath).iloc[:, [0,4,6]]
df.head()

	carat	depth	price
0	0.23	61.5	326
1	0.21	59.8	326
2	0.23	56.9	327
3	0.29	62.4	334
4	0.31	63.3	335

데이터를 표준화합니다.

scal=StandardScaler().fit(df)
data=scal.transform(df)
data[:10]

array([[-1.19816781, -0.17409151, -0.90409516],
       [-1.24036129, -1.36073849, -0.90409516],
       [-1.19816781, -3.38501862, -0.9038445 ],
       [-1.07158736,  0.45413336, -0.90208985],
       [-1.02939387,  1.08235823, -0.90183918],
       [-1.17707106,  0.73334442, -0.90158852],
       [-1.17707106,  0.3843306 , -0.90158852],
       [-1.13487758,  0.10511955, -0.90133785],
       [-1.21926455,  2.33880797, -0.90133785],
       [-1.19816781, -1.63994954, -0.90108719]])

각 샘플과 평균과의 거리, 공분산의 역행렬을 계산합니다.

xm=data-data.mean(axis=0)
cov=np.cov(data, rowvar=False)
covi=la.inv(cov)

공분산행렬에서 각 변수에 관게된 분산은 주대각요소들입니다. 그러므로 최종 결과에서 주 대각요소들의 제곱근이 마할로비스 거리가 됩니다.

d1=xm@covi
d2=d1@xm.T
D=np.sqrt(np.diagonal(d2))
D[:10]

array([1.30761617, 1.88151452, 3.56581283, 1.20601365, 1.53207029,
       1.50750889, 1.3504983 , 1.20689511, 2.78642814, 2.03875282])

위 과정을 함수로 작성합니다.

def mahalDist(x):
    #x는 array(각 변수는 열로 이루어짐)
    m=x.mean(axis=0)
    x_m=x-m
    cov=np.cov(x, rowvar=False)
    covi=np.linalg.inv(cov)
    left=np.dot(x_m, covi)
    tot=np.dot(left, x_m.T)
    return np.sqrt(np.diagonal(tot))

d=mahalDist(data)
d[:10]

array([1.30761617, 1.88151452, 3.56581283, 1.20601365, 1.53207029,
       1.50750889, 1.3504983 , 1.20689511, 2.78642814, 2.03875282])

마할로비스 거리를 적용하여 이상치(outlier)의 감지합니다.

데이터가 n 자유도를 가진 카이제곱분포를 따른다고 가정할 수 있습니다. 이 가정에서 유도수준 0.01, 자유도 2에서의 임계값을 계산하면 다음과 같습니다. 이 계산은 scipy.ch2.ppf() 메서드를 적용합니다. ppf(percent point function)는 cdf의 역수로 전달한 면적값(누적확률)에 대응하는 x축의 값을 반환합니다.

from scipy.stats import chi2

criticalPoint=chi2.ppf((1-0.01), df=2)
criticalPoint

9.21034037197618

일반적으로 자유도 2인 Χ² 분포는 다음과 같으며 유의수준 0.01인 점선을 기준으로 오른쪽에 위치한 값들은 극단적인 값으로 간주할 수 있습니다.

x=np.sort(np.linspace(0, 20))
y=chi2.pdf(x, df=2)

plt.figure(figsize=(4,3))
plt.plot(x, y)
plt.vlines(criticalPoint, 0, 0.5, ls="dashed", color="orange")
plt.xlabel("mahalnobis's distance")
plt.ylabel("probability")
plt.show()

각 거리에 대한 p-value를 계산합니다.

df1=df.copy()
pval=1-chi2.cdf(d, df=2)
df1["mahalD"]=d 
df1["p-value"]=pval
df1.head(3)

	carat	depth	price	mahalD	p-value
0	0.23	61.5	326	1.307616	0.520062
1	0.21	59.8	326	1.881515	0.390332
2	0.23	56.9	327	3.565813	0.168149

데이터의 각 샘플에서 계산한 p-value가 0.01보다 작은 경우 유의미한 차이를 보인다고 할 수 있습니다. 즉, 자유도 2인 Χ² 분포에 부합하지 않는다고 할 수 있습니다. 결과적으로 이상치로 간주할 수 있습니다.

np.where(pval<0.01)

(array([ 4518,  6341, 10377, 16283, 17196, 19339, 19346, 21758, 22741,
        23644, 24328, 25998, 25999, 26444, 27130, 27415, 27630, 41918,
        52860, 52861], dtype=int64),)