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

#1.はじめに
前回の記事[【異常検知】マハラノビス距離を嚙み砕いて理解する (1)](https://qiita.com/kotai2003/items/297930db4466b71f06b0)の続きです。
今回は、実際のデータセットを利用して、マハラノビス距離を計算してみます。

#2.データセット(Davis)
有名なDavisデータセットを使います。元を辿ると1990年の論文のデータです。
C.Davis, "Body image and weight preoccupation: A comparison between exercising and non-exercising women" (1990) [Link](https://www.sciencedirect.com/science/article/abs/pii/019566639090096Q?via%3Dihub)

データCSVのダウンロードはこちらです。
[Davis Dataset](https://www.kaggle.com/ravinduabey/davis-data-set)
#データの集計（平均、分散）
　データを覗くと、Featureは['sex', 'weight', 'height', 'repwt', 'repht']の五つで、データの長さは200です。従ってshapeは(5,200)

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

```python:
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)

![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/1f90d88d-f5b3-cb40-7fa5-3dc5974ffb38.png)


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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/f0862652-c8b9-d720-140e-7ed4cf4bd26d.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/1144305f-aec4-eeb9-6624-48d807e8fb76.png)

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

```python:
#統計情報の集計

#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)

```

```python:結果
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にマハラノビスを計算する関数があるので、それを利用します。

```python:

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])
```
プロットしてみます。

```python:
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()

```

![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/6f1be7e4-df9a-c4ad-d1d8-98988d639c88.png)

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

```python:
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を利用し、一緒に描画するだけです。

```python:

#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()

```
![Figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/2427af32-559e-11f0-6e36-1c249239f106.png)

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

#5. 参考資料
1. [井出剛さんのTwitter ](https://twitter.com/Idesan/status/1469385932742664195
)
1. [入門 機械学習による異常検知―Rによる実践ガイド](https://www.amazon.co.jp/%E5%85%A5%E9%96%80-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%81%AB%E3%82%88%E3%82%8B%E7%95%B0%E5%B8%B8%E6%A4%9C%E7%9F%A5%E2%80%95R%E3%81%AB%E3%82%88%E3%82%8B%E5%AE%9F%E8%B7%B5%E3%82%AC%E3%82%A4%E3%83%89-%E4%BA%95%E6%89%8B-%E5%89%9B/dp/4339024910)
1. [On the generalized distance in statistics](http://bayes.acs.unt.edu:8083/BayesContent/class/Jon/MiscDocs/1936_Mahalanobis.pdf)
1. [マハラノビス距離の意味を２次元の場合で理解する](https://mathwords.net/mahalanobis)
1. [異常検知 マハラノビス距離を実装してみる](https://hituji-ws.com/code/python/md_anomaly/)
1. [Maharanobis Distance at Scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.mahalanobis.html)
1. [【python できること】異常検知 マハラノビス距離を実装してみる](https://hituji-ws.com/code/python/md_anomaly/)