#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.結論
マハラノビス距離、知れば知るほどすごいですね。多次元データの異常検知では、やはりこれを使わない手はないと思います。
#5. 参考資料