Pythonデータ分析勉強会　第3回資料

今回は第3回目です。

[Pythonデータ分析勉強会　第1回資料](https://qiita.com/shinmura0/items/1a1549e832a3664b0eac)
[Pythonデータ分析勉強会　第2回資料](https://qiita.com/shinmura0/items/57957b1c757d3cb44be9)
Pythonデータ分析勉強会　第3回資料

第1回、第2回の資料を見れば、Pythonを使ってデータ分析が進められると思います。

今回は、MT法の説明がメインです。機械学習に興味がある方は、本稿を飛ばして「scikit-learn」を、
ディープラーニングに興味がある方は「Keras」を勉強すると、理解が進むと思います。

※　塩尻市内で[Pythonデータ分析勉強会](https://pythondata.connpass.com/event/101107/)を開催します。**初心者の方も是非参加してください！**

![figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/264781/88c0640a-8dc9-5dc1-9470-d6fcd0f68861.png)

#第2回　演習の答え
##その1

```py
import numpy as np
import matplotlib.pyplot as plt

A = np.sin(2*np.pi*np.array(range(100))/100)

plus, minus = [], []

for i in A:
    if i > 0:
        plus.append(i)
    else:
        minus.append(i)

np.savetxt("All.csv", A, delimiter=",")
np.savetxt("plus.csv", plus, delimiter=",")
np.savetxt("minus.csv", minus, delimiter=",")

plt.figure()
plt.title("All")
plt.plot(A)
plt.show()

plt.figure()
plt.title("Plus")
plt.plot(plus)
plt.show()

plt.figure()
plt.title("Minus")
plt.plot(minus, c="red")
plt.show()
```

##その2

```py
import numpy as np
import matplotlib.pyplot as plt

def plot_csv(name):
    x = np.loadtxt(name + ".csv", delimiter=",")
    
    plt.figure()
    plt.title(name)
    plt.plot(x)
    plt.show()
      
plot_csv("All")
plot_csv("Plus")
plot_csv("Minus")
```

#MT法

##概要

MT法は、一言でいうと「変数間の相関を考慮した距離を出し、その距離を見て異常かどうか判定する」
手法といえます。図にすると以下のようになります。

![figure_5.png](https://qiita-image-store.s3.amazonaws.com/0/264781/e9ab714a-91f7-fa0f-53cb-50379435d4f0.png)

正常データの平均値を出して、平均値からの距離で異常検知することを考えます。
通常の距離（ユークリッド距離）を使うと左図にあるように、同心円状に広がっていきます。
ところが、変数間の相関を考慮した距離（マハラノビス距離）を使うと、右図のように楕円上に
広がる距離となります。

MT法では、マハラノビス距離が大きければ「異常」、小さければ「正常」と判定します。

##理論

MT法の理論は、線形代数の知識が必要になります。ただ、理論を理解しなくても使うことはできるので、
お急ぎの方は、コード全文をご覧ください。MT法の数式は以下の資料が参考になります。

・http://kaz7227.art.coocan.jp/mt-system.pdf
・[入門MTシステム](https://www.amazon.co.jp/%E5%85%A5%E9%96%80-MT%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0-%E7%AB%8B%E6%9E%97-%E5%92%8C%E5%A4%AB/dp/4817192844)

MT法は以下の流れで行っています。

![無題.jpg](https://qiita-image-store.s3.amazonaws.com/0/264781/077292b9-fe52-4ca6-791c-50e2be581854.jpeg)

まずは、左図にあるように正常データで平均値、標準偏差、相関行列を求め、単位空間を作成します。
そして、テストデータでマハラノビス距離を求め、その距離で異常かどうか判断します。

ここで、テストデータXが以下のように与えられたとします。

```math
X=[X_{1},X_{2},・・・,X_{n}]
```

Xを使って、以下のようにマハラノビス距離が求められます。

```math
\tilde{X} = [\frac{X_{1}-\bar{X_{1}}}{\sigma_{1}},\frac{X_{2}-\bar{X_{2}}}{\sigma_{2}},・・・,\frac{X_{n}-\bar{X_{n}}}{\sigma_{n}}]\\

D^2 = \frac{1}{n}\tilde{X}R^{-1}\tilde{X}^{T}
```

ただし、$\bar{X}$は各状態量の平均値、$\sigma$は各状態量の標準偏差、$R$は相関行列、$D$はマハラノビス距離です。

###例題

ここでは、以下のデータが与えられたとして単位空間を作ってみます。

||X|Y|
|---|---|---|
|データ1|10|20|
|データ2|11|18|
|データ3|9.5|22|

XとYは状態量です。例えば、Xは機械の温度、Yは機械にかかる力と見なしても良いでしょう。

###平均値を求める

最初に各状態量の平均値を求め、その平均値で引くという作業が必要になります。
表にすると、以下のとおりです。

・平均値を求める

||X|Y|
|---|---|---|
|データ1|10|20|
|データ2|11|18|
|データ3|9.5|22|
|**平均値(avg)**|10.2|20|

・平均値で引く

||XX|YY|
|---|---|---|
|データ1|10-10.2=**-0.2**|20-20=**0**|
|データ2|11-10.2=**0.8**|18-20=**-2**|
|データ3|9.5-10.2=**-0.7**|22-20=**2**|
上の表で、新たな変数XX、YYを作っています。

コードにすると以下のとおりです。

```py
import numpy as np

xx = np.copy(x)
avg = np.zeros(2)

#各状態量から平均値を引く
for i in range(2):
    avg[i] = np.mean(x[:,i])
    for j in range(3):
        xx[j,i] = x[j,i] - avg[i]
```

・xというのは、データが入ったnp配列（3行、2列）です。
・3行目でxをコピーしています。
・4行目で平均値配列(avg)の形を宣言しています。
・7行目以降はfor文です。rangeに2を与えているのは、xの列の数が「2」だからです。
・8行目で各列（状態量）毎に平均値を求めています。
・9行目以降もfor文になっています。7行目で、rangeに3を与えているのは、xの行数が「3」だからです。
・最後の行で、xから各平均値で引き、XXとYYを求めています。

###標準偏差を求める
次に、前項で求めた値を「各状態量の標準偏差で割る」という作業が必要です。
表にすると、以下のとおりです。

・標準偏差を求める

||X|Y|
|---|---|---|
|データ1|10|20|
|データ2|11|18|
|データ3|9.5|22|
|**標準偏差(std)**|0.62|1.63|
Xの標準偏差

```math
\sqrt{\frac{(10-10.2)^2+(11-10.2)^2+(9.5-10.2)^2}{3}}=0.62
```

Yの標準偏差

```math
\sqrt{\frac{(20-20)^2+(18-20)^2+(22-20)^2}{3}}=1.63
```

・標準偏差で割る

||XX|YY|
|---|---|---|
|データ1|-0.2/0.62=**-0.32**|0/1.63=**0**|
|データ2|0.8/0.62=**1.29**|-2/1.63=**-1.23**|
|データ3|-0.7/0.62=**-1.13**|2/1.63=**1.23**|

コードにすると、以下のとおりです。

```py
std = np.zeros(2)

#標準偏差で割る
for i in range(2):
    std[i] = np.std(x[:,i])
    for j in range(3):
        xx[j,i] = xx[j,i] / std[i]
```

・やっていることは、平均値の作業と同じです。「np.mean」が「np.std」になっただけです。
・標準偏差はnp.std()の1行で求められます。これは便利！



###相関行列を求める
相関行列は、以下の式で与えられます。

```math
R = \begin{pmatrix}
1&r_{yx}\\
r_{xy}&1
\end{pmatrix}
```

ただし、$r_{xy}$と$r_{yx}$はxとｙの相関係数です。マハラノビス距離を求める際は、逆行列$R^{-1}$を使います。
コードにすると以下のとおりです。

```py
R = np.corrcoef(x.transpose())
invR = np.linalg.inv(R)
```

ただし、np.corrcoef()に渡す配列は、横長の形にしないといけないため、transpose()で
転置しています。$R^{-1}$は以下のように、求まります。

```
array([[28.        , 27.49545417],
       [27.49545417, 28.        ]])
```

これで、マハラノビス距離を求める準備が整いました。
次に、正常データでマハラノビス距離を求めてみます。

###マハラノビス距離を求める
試しに、データ1のマハラノビス距離を求めてみます。

||XX|YY|
|---|---|---|
|データ1|-0.32|0|

マハラノビス距離は前述したとおり、以下の式で与えられます。

```math
D^2 = \frac{1}{n}\tilde{X}R^{-1}\tilde{X}^{T}
```

データ1の数値を入れると、以下のとおりです。

```math
D^2=\frac{1}{2}

\begin{pmatrix}
-0.32&0
\end{pmatrix}

\begin{pmatrix}
28&27.495\\
27.495&28\\
\end{pmatrix}

\begin{pmatrix}
-0.32\\
0
\end{pmatrix}

=1.4336
```


コードで書くと、以下のとおりです。

```py
#MD^2の計算
d0 = np.array([-0.32, 0])
d1 = np.dot(d0,invR)
d2 = np.dot(d1,d0)/2
```

###等確率楕円

グラフで出てくる楕円は「等確率楕円」と呼ばれ、以下の資料（P2）が参考になります。
http://civilyarou.web.fc2.com/WANtaroHP_html5_win/f90_STAT/dir_SREG/TeX_ellipse.pdf

コードは以下のとおりです。

```py
curve_c = np.zeros((2,51))

low = np.corrcoef(data[:,0],data[:,1])[0,1]

for i in range(div+1):
    r = (-2*(1-low**2)*np.log(1-0.95)/(1-2*low*np.sin(i*2*np.pi/50)*np.cos(i*2*np.pi/50)))**0.5
    curve_c[0,i] = avg[0] + std[0]*r*np.cos(i*2*np.pi/50)
    curve_c[1,i] = avg[1] + std[1]*r*np.sin(i*2*np.pi/50)
```

・1行目で、楕円カーブの座標が入る箱を用意しています。
・3行目で相関係数を求め、for文以降は参考資料のとおりに計算しています。

#Python上達のコツ
・Python上達の近道は、実装して動かしてみることです。ネットには、先人達が書いた素晴らしい
　コードがたくさん落ちています。そのコードをいじって**データを使い**練習しましょう。
　使うデータは、乱数でもフリーのデータでも何でもOKです。

・Pythonのコードを書いて、初めて動かすときは、必ずバグがあると思ってください。
　最初は、エラーが出るのが当たり前です。エラーが出なくなっても、print文で
　「データ数」や「出力値」を出して、**論理的に間違っていないか**確認をとってください。

・論理的な確認をとったら、**数種類のダミーデータで**確認をとってください。
　コードは完璧！と思っていても思わぬ見落としがあったりします。

#コード全文

ちなみに、コード中の「if __name__ == '__main__':」というのは、ここからコードが始まりますよ！という
宣言です。無くても良いですが、見やすさのために書いています。

```py
import numpy as np
import matplotlib.pyplot as plt

divide = 2 #状態量の数
R = np.zeros((divide,divide))    #相関行列
invR = np.zeros((divide,divide)) #相関行列の逆行列
avg = np.zeros(divide)           #平均値
std = np.zeros(divide)           #標準偏差
make = 0

p = 0.95        #マハラノビス距離p=0.95で2σ
md_sikii = 2.448#MDの閾値95%で2.448
div = 50        #Mt楕円の分割数
  
def maha(x):
    global make, R, invR, avg, std
    
    N, _ = x.shape #Nはデータ数
    xx = np.copy(x)
    xx = np.array(xx,dtype="float32")
    x_return = []
    
    #各状態量から平均値を引く
    for i in range(divide):
        if make == 0:
            avg[i] = np.mean(x[:,i])
        for j in range(N):
            xx[j,i] = xx[j,i] - avg[i]

    #各状態量を標準偏差で割る
    for i in range(divide):
        if make == 0:
            std[i] = np.std(x[:,i])
        for j in range(N):
            xx[j,i] = xx[j,i] / std[i]
    
    #make=0のときだけ計算
    if make == 0:
        R = np.corrcoef(xx.transpose())
        invR = np.linalg.inv(R)
        make = 1
        
    #MD^2の計算
    for i in range(N):
        d0 = xx[i,:]
        d1 = np.dot(d0,invR)
        d2 = np.dot(d1,d0)/divide
        x_return.append(d2)
        
    return x_return

if __name__ == '__main__':
    
    curve_c = np.zeros((2,div+1))

    #正常データの作成
    x1 = np.random.normal(1, 0.3, (1, 50))
    y1 = np.random.normal(1, 0.3, (1, 50))
    x2 = np.random.normal(1.5, 0.3, (1, 50))
    y2 = np.random.normal(1.5, 0.3, (1, 50))
    
    #テストデータの作成
    test1 = np.array([1.02,1.5])
    test1 = test1.reshape((1,2))
    test2 = np.array([0.5,2])
    test2 = test2.reshape((1,2))

    #正常データの形を整える
    data = []
    data.append(x1)
    data.append(x2)
    data.append(y1)
    data.append(y2)
    data = np.array(data)
    data = data.reshape(2,100)
    data = data.transpose()
    print(data.shape)
    
    #単位空間の作成
    _ = maha(data)

    #テストデータのマハラノビス距離
    md_1 = maha(test1)
    md_2 = maha(test2)
    
    #楕円のデータ
    low = np.corrcoef(data[:,0],data[:,1])[0,1]

    for i in range(div+1):
        r = (-2*(1-low**2)*np.log(1-p)/(1-2*low*np.sin(i*2*np.pi/div)*np.cos(i*2*np.pi/div)))**0.5
        curve_c[0,i] = avg[0] + std[0]*r*np.cos(i*2*np.pi/div)
        curve_c[1,i] = avg[1] + std[1]*r*np.sin(i*2*np.pi/div)

    #可視化
    plt.figure()

    plt.subplot(1,2,1)
    plt.scatter(x1, y1, c="green", s=50)
    plt.scatter(x2, y2, c="green", s=50)
    plt.scatter(test1[:,0],test1[:,1],c="m", s=50)
    plt.scatter(test2[:,0],test2[:,1],c="red", s=50)
    plt.xlabel("T")
    plt.ylabel("KN")
    
    plt.plot(curve_c[0],curve_c[1],c="c",label="MD^2=2.448")
    plt.legend()
    
    plt.subplot(1,2,2)
    plt.bar([1,2], [md_1[0], md_2[0]], align="center")
    plt.grid(True)
    plt.xticks([1,2], ["Purple", "Red"])
    plt.ylabel("MD^2")
    plt.grid(True)

    plt.show()
```

コードを実行すると、冒頭で出た散布図と棒グラフが出力されます。