1.はじめに

前回の記事【異常検知】マハラノビス距離を嚙み砕いて理解する (1)の続きです。
今回は、実際のデータセットを利用して、マハラノビス距離を計算してみます。

2.データセット(Davis)

有名なDavisデータセットを使います。元を辿ると1990年の論文のデータです。
C.Davis, "Body image and weight preoccupation: A comparison between exercising and non-exercising women" (1990) Link

データCSVのダウンロードはこちらです。
Davis Dataset

データの集計（平均、分散）

　データを覗くと、Featureは['sex', 'weight', 'height', 'repwt', 'repht']の五つで、データの長さは200です。従ってshapeは(5,200)

　この中で'weight', 'height'を選び、マハラノビス距離を計算します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.spatial import distance

#Davis Dataset
df = pd.read_csv('Davis.csv')
print(df.head())
print(df.columns)
df = df[['weight', 'height']]

data = df.to_numpy()


print('data',data.shape)
print('data.T',data.T.shape)
print(data.T[0])

#Plot
fig1 = plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot()
ax1.scatter(data.T[0], data.T[1])
ax1.set_title('Davis Data')
ax1.set_xlabel('weight')
ax1.set_ylabel('height')
ax1.set_aspect('equal')
ax1.grid()
plt.show()

怪しい奴が二つ見えてますね。
(w,h) = (120,180)
(w,h) = (165,50)

3.マハラノビス距離の計算

それでは、このデータ('weight', 'height'の2組（＝次元）のデータ群）のマハラノビス距離を計算してみましょう。
マハラノビス距離の式はこちらです。

必要なパラメータを計算しましょう。

#統計情報の集計

#mean
mu_mat = np.mean(data, axis=0)
print(mu_mat)

#data - mean
data_m_mat = data - mu_mat

#covariance matrix
cov_mat = np.cov(data.T)
cov_i_mat = np.linalg.pinv(cov_mat)
print(cov_mat)
print(cov_i_mat)

結果

print(mu_mat)
[ 65.8  170.02]

print(cov_mat)
[[227.85929648  34.3758794 ]
 [ 34.3758794  144.19055276]]

print(cov_i_mat)
[[ 0.00455241 -0.00108532]
 [-0.00108532  0.00719401]]

3.1.scipyの関数利用

scipy.spacialのdistanceにマハラノビスを計算する関数があるので、それを利用します。


mahala_result = []
for i in range(len(data_m_mat)):
    mahala_result.append(distance.mahalanobis(data[i],mu_mat, cov_i_mat ))
    print(data[i])

プロットしてみます。

fig2 = plt.figure(figsize=(8,5))
ax1 = fig2.add_subplot()
ax1.plot(mahala_result)
ax1.set_title('Mahalanobis Distance')
ax1.set_xlabel('sample')
ax1.set_ylabel('Mahalanobis Distance')
# ax1.set_aspect('equal')
ax1.grid()

3.2.numpyの利用

同様にnumpyを利用し、直接にマハラノビス距離を計算することも可能です。こちらのほうが早い気がします。

mahala_result = np.sqrt(np.sum(np.dot(data_m_mat,cov_i_mat)*data_m_mat,axis = 1))

4.等高線図の作成

それでは、見やすくするため、マハラノビス距離の等高線を書いてみます。
やり方は、x座標、y座標、そして、(x,y)座標に対応するマハラノビス距離を格納する3次元のデータを作成します。(x座標、y座標、マハラノビス距離)の3組（次元）データ。 それをmatplotのcontourを利用して描画します。
その後、scatterを利用し、一緒に描画するだけです。


#Contour

# weight [30,200]
# height [50,210]


z = np.array(
    [
        [(i,j, distance.mahalanobis([i,j], mu_mat, cov_i_mat))
         for i in np.linspace(30,200,100)]
        for j in np.linspace(50,210,100)
    ]

)

print('z shape',z.shape)
print(z[:,:,0])

fig3 = plt.figure()
ax1 = fig3.add_subplot()
ax1.set_title('Davis Data, Mahalanobis Distance')
CS = ax1.contour(z[:,:,0], z[:,:,1],z[:,:,2])
ax1.clabel(CS)
ax1.set_xlabel('weight')
ax1.set_ylabel('height')
ax1.set_aspect('equal')

ax1.grid()

ax1.scatter(data.T[0], data.T[1])

plt.show()

5.結論

マハラノビス距離、知れば知るほどすごいですね。多次元データの異常検知では、やはりこれを使わない手はないと思います。

【異常検知】マハラノビス距離を嚙み砕いて理解する (2)