3行で

・井出・杉山の異常検知と変化検知で説明されてるアルゴリズムを実問題(に近いなにか)に適用したい
・まずはホテリングのT2法をkaggleのCredit Card Fraud Detectionをやってみたぞ
https://www.kaggle.com/mlg-ulb/creditcardfraud

理論の説明

詳細は本をご購入ください。

アルゴリズムをざっと説明すると

  1. (パラメータ決定)異常度の分布であるカイ2乗分布から異常と判定する閾値を決定する
  2. (訓練)(正常標本が圧倒的に多いデータから)標本平均、標本共分散行列を計算する
  3. (推論1)2で求めた標本共分散行列で定義される標本平均と推論対象のデータのマハラノビス距離を計算し、その推論対象の異常度とする
  4. (推論2)計算された異常度が閾値を超えたものを異常と判定する

データのロード

In [18]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split
import scipy.stats

df = pd.read_csv("creditcard.csv")

あとで使うモジュールもここでimportしておく

データ概観

In [17]:
df.head(3)
Out[17]:
Time V1 V2 V3 V4 V5 V6 V7 V8 V9 ... V21 V22 V23 V24 V25 V26 V27 V28 Amount Class
0 0.0 -1.359807 -0.072781 2.536347 1.378155 -0.338321 0.462388 0.239599 0.098698 0.363787 ... -0.018307 0.277838 -0.110474 0.066928 0.128539 -0.189115 0.133558 -0.021053 149.62 0
1 0.0 1.191857 0.266151 0.166480 0.448154 0.060018 -0.082361 -0.078803 0.085102 -0.255425 ... -0.225775 -0.638672 0.101288 -0.339846 0.167170 0.125895 -0.008983 0.014724 2.69 0
2 1.0 -1.358354 -1.340163 1.773209 0.379780 -0.503198 1.800499 0.791461 0.247676 -1.514654 ... 0.247998 0.771679 0.909412 -0.689281 -0.327642 -0.139097 -0.055353 -0.059752 378.66 0

3 rows × 31 columns

V1からV28はカード会社の持っているデータをPCAしたもののようです。
生データは機密だからみせられないよ!とのこと。
PCAされているので共分散行列は対角化されていると期待できる。

In [16]:
df.describe()
Out[16]:
Time V1 V2 V3 V4 V5 V6 V7 V8 V9 ... V21 V22 V23 V24 V25 V26 V27 V28 Amount Class
count 284807.000000 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 ... 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 2.848070e+05 284807.000000 284807.000000
mean 94813.859575 1.165980e-15 3.416908e-16 -1.373150e-15 2.086869e-15 9.604066e-16 1.490107e-15 -5.556467e-16 1.177556e-16 -2.406455e-15 ... 1.656562e-16 -3.444850e-16 2.578648e-16 4.471968e-15 5.340915e-16 1.687098e-15 -3.666453e-16 -1.220404e-16 88.349619 0.001727
std 47488.145955 1.958696e+00 1.651309e+00 1.516255e+00 1.415869e+00 1.380247e+00 1.332271e+00 1.237094e+00 1.194353e+00 1.098632e+00 ... 7.345240e-01 7.257016e-01 6.244603e-01 6.056471e-01 5.212781e-01 4.822270e-01 4.036325e-01 3.300833e-01 250.120109 0.041527
min 0.000000 -5.640751e+01 -7.271573e+01 -4.832559e+01 -5.683171e+00 -1.137433e+02 -2.616051e+01 -4.355724e+01 -7.321672e+01 -1.343407e+01 ... -3.483038e+01 -1.093314e+01 -4.480774e+01 -2.836627e+00 -1.029540e+01 -2.604551e+00 -2.256568e+01 -1.543008e+01 0.000000 0.000000
25% 54201.500000 -9.203734e-01 -5.985499e-01 -8.903648e-01 -8.486401e-01 -6.915971e-01 -7.682956e-01 -5.540759e-01 -2.086297e-01 -6.430976e-01 ... -2.283949e-01 -5.423504e-01 -1.618463e-01 -3.545861e-01 -3.171451e-01 -3.269839e-01 -7.083953e-02 -5.295979e-02 5.600000 0.000000
50% 84692.000000 1.810880e-02 6.548556e-02 1.798463e-01 -1.984653e-02 -5.433583e-02 -2.741871e-01 4.010308e-02 2.235804e-02 -5.142873e-02 ... -2.945017e-02 6.781943e-03 -1.119293e-02 4.097606e-02 1.659350e-02 -5.213911e-02 1.342146e-03 1.124383e-02 22.000000 0.000000
75% 139320.500000 1.315642e+00 8.037239e-01 1.027196e+00 7.433413e-01 6.119264e-01 3.985649e-01 5.704361e-01 3.273459e-01 5.971390e-01 ... 1.863772e-01 5.285536e-01 1.476421e-01 4.395266e-01 3.507156e-01 2.409522e-01 9.104512e-02 7.827995e-02 77.165000 0.000000
max 172792.000000 2.454930e+00 2.205773e+01 9.382558e+00 1.687534e+01 3.480167e+01 7.330163e+01 1.205895e+02 2.000721e+01 1.559499e+01 ... 2.720284e+01 1.050309e+01 2.252841e+01 4.584549e+00 7.519589e+00 3.517346e+00 3.161220e+01 3.384781e+01 25691.160000 1.000000

8 rows × 31 columns

stdをみると、V1からV28まで降順で並んでいるのがわかる
V1からV28はPCAしたときの共分散行列の固有値の大きい順に並んでいるようだ。
また、meanを見るとほとんど0に近い値が並んでいる
PCAしているので平均が0になるようにシフトされている

In [23]:
df.cov()
Out[23]:
Time V1 V2 V3 V4 V5 V6 V7 V8 V9 ... V21 V22 V23 V24 V25 V26 V27 V28 Amount Class
Time 2.255124e+09 1.091960e+04 -8.307031e+02 -3.021425e+04 -7.077378e+03 1.134407e+04 -3.986868e+03 4.976739e+03 -2.095683e+03 -4.518322e+02 ... 1.560435e+03 4.964595e+03 1.516599e+03 -4.654076e+02 -5.769855e+03 -9.482254e+02 -9.841860e+01 -1.475443e+02 -125860.970747 -24.300717
V1 1.091960e+04 3.836489e+00 -3.668399e-16 -1.591003e-15 -7.327817e-16 8.378639e-16 3.333094e-16 9.924436e-17 -5.408793e-17 8.332734e-17 ... -1.691534e-16 1.559269e-16 2.421982e-16 -5.106918e-17 -3.270723e-16 -1.284838e-16 9.270791e-17 2.291751e-16 -111.556566 -0.008244
V2 -8.307031e+02 -3.668399e-16 2.726820e+00 1.722431e-16 -5.294779e-16 1.705466e-16 8.324751e-16 -2.573168e-16 4.291109e-18 -1.894075e-16 ... 8.721928e-17 2.080689e-16 9.431707e-17 -1.081634e-16 1.314776e-16 2.006343e-16 -3.456338e-16 -1.968671e-16 -219.485433 0.006260
V3 -3.021425e+04 -1.591003e-15 1.722431e-16 2.299029e+00 -4.621923e-16 -1.297100e-15 2.939634e-15 5.082718e-16 -5.758069e-17 2.412002e-16 ... -1.889959e-16 -2.348384e-16 -7.394678e-17 3.089224e-17 5.925223e-17 -1.715297e-16 3.409186e-16 3.703327e-16 -79.975549 -0.012150
V4 -7.077378e+03 -7.327817e-16 -5.294779e-16 -4.621923e-16 2.004684e+00 -3.576902e-15 -7.963000e-16 -1.949710e-16 1.121726e-15 9.433454e-16 ... -5.786136e-17 2.599239e-16 1.700352e-16 1.347708e-16 4.773297e-16 -2.689739e-16 -6.540199e-17 -4.341005e-18 34.964556 0.007846
V5 1.134407e+04 8.378639e-16 1.705466e-16 -1.297100e-15 -3.576902e-15 1.905081e+00 1.156005e-15 -1.397105e-18 8.149115e-16 8.041837e-16 ... -7.223938e-17 4.573025e-17 9.701149e-17 -8.131900e-16 -6.929392e-17 2.398717e-16 2.521276e-16 -8.280967e-17 -133.380790 -0.005444
V6 -3.986868e+03 3.333094e-16 8.324751e-16 2.939634e-15 -7.963000e-16 1.156005e-15 1.774946e+00 -5.304010e-17 -5.904266e-16 -1.975345e-16 ... -9.924436e-17 -1.017080e-16 3.371062e-17 -8.535315e-16 3.820521e-16 -1.645278e-16 -6.873882e-17 1.897350e-16 71.970931 -0.002415
V7 4.976739e+03 9.924436e-17 -2.573168e-16 5.082718e-16 -1.949710e-16 -1.397105e-18 -5.304010e-17 1.530401e+00 -5.867842e-17 3.567608e-17 ... 3.142395e-17 -5.758880e-16 -2.145430e-16 -3.648690e-18 2.201688e-17 -4.844150e-16 -5.993831e-17 3.056168e-17 122.936845 -0.009620
V8 -2.095683e+03 -5.408793e-17 4.291109e-18 -5.758069e-17 1.121726e-15 8.149115e-16 -5.904266e-16 -5.867842e-17 1.426479e+00 5.726634e-16 ... 3.613762e-17 1.493156e-17 1.519976e-16 -1.545859e-16 -1.134399e-16 -7.983458e-19 1.586213e-16 -2.243975e-16 -30.792991 0.000986
V9 -4.518322e+02 8.332734e-17 -1.894075e-16 2.412002e-16 9.433454e-16 8.041837e-16 -1.975345e-16 3.567608e-17 5.726634e-16 1.206992e+00 ... 1.899876e-16 -1.223621e-16 -6.412962e-17 -1.877797e-16 1.144129e-16 -7.271184e-17 -7.746449e-17 2.986687e-16 -12.158248 -0.004459
V10 1.583108e+03 1.150616e-16 -2.581651e-16 3.983746e-16 -1.549290e-16 1.886092e-16 1.885094e-16 4.369446e-16 5.214196e-17 -3.406442e-16 ... 8.314273e-16 -2.321544e-16 1.570122e-16 -6.012542e-17 -1.976904e-16 -2.011457e-16 -1.489913e-16 8.499888e-17 -27.643420 -0.009807
V11 -1.200595e+04 6.134539e-16 5.897780e-16 1.855406e-16 -3.985991e-16 1.059954e-15 1.252530e-15 -4.654107e-16 1.648335e-16 3.325110e-16 ... 6.411715e-18 7.422121e-18 7.518562e-17 1.014404e-15 -3.322101e-16 -5.055229e-17 -6.875130e-17 -1.178776e-16 0.026545 0.006565
V12 5.900343e+03 2.879035e-16 -5.113904e-16 2.192457e-16 -2.237614e-16 4.985046e-16 3.857382e-16 8.852408e-16 6.361818e-17 -1.452939e-15 ... 4.281192e-16 -4.191939e-17 1.751683e-16 2.662171e-16 5.036439e-18 7.535948e-17 -1.251532e-16 2.307594e-16 -2.384691 -0.010813
V13 -3.114775e+03 -1.150366e-16 -2.421233e-17 -1.552408e-17 6.985526e-18 -3.766882e-16 -2.110751e-16 -9.082431e-17 -3.522666e-16 1.031550e-15 ... 1.101780e-16 -3.090471e-17 -4.245672e-16 -3.867362e-16 -3.703108e-17 -6.879807e-17 -1.954029e-16 3.569869e-16 1.317731 -0.000189
V14 -4.495601e+03 6.929642e-16 -5.958903e-16 9.764268e-16 -1.340223e-16 3.198373e-16 4.424582e-16 3.789648e-17 -2.466140e-16 1.055974e-15 ... -1.480713e-16 4.201232e-16 1.360837e-16 5.972001e-19 -1.113942e-17 -1.394026e-17 8.638351e-18 7.615907e-16 8.092317 -0.012044
V15 -7.974101e+03 -1.809002e-16 1.848171e-16 -9.166008e-17 2.706080e-16 1.742967e-16 -1.507127e-16 -7.203823e-17 1.270056e-16 -9.079687e-16 ... 3.721040e-17 -2.108506e-16 6.455862e-17 -2.440054e-16 1.078079e-16 4.659096e-17 -4.597287e-16 -3.402512e-16 -0.683577 -0.000161
V16 4.952977e+02 5.285673e-16 3.155961e-17 7.528650e-16 -9.837117e-17 7.070849e-16 -1.168329e-16 4.634148e-16 1.289328e-16 -4.782091e-16 ... -2.574385e-16 1.641723e-16 3.941458e-16 -1.879231e-16 -1.675341e-16 -1.998359e-16 2.856098e-16 2.027221e-16 -0.856845 -0.007152
V17 -2.956329e+03 -9.979323e-20 -8.274355e-16 2.238861e-16 -4.450404e-16 5.438107e-16 1.515578e-16 5.852623e-16 -3.612764e-16 7.218917e-16 ... -5.982261e-16 -1.709583e-16 2.318197e-16 -8.313399e-17 5.065130e-17 1.111198e-16 2.314862e-16 -2.508123e-17 1.552706 -0.011515
V18 3.599748e+03 2.506182e-16 3.313010e-16 3.611642e-16 -1.999607e-17 4.699014e-16 2.479862e-17 1.761132e-16 -3.275299e-16 1.275108e-16 ... -7.473702e-16 -3.439810e-16 -1.539934e-16 -9.610399e-17 -1.062299e-16 1.221656e-16 7.003925e-17 2.214466e-16 7.473906 -0.003880
V19 1.120106e+03 2.628616e-16 1.322260e-17 4.161378e-16 -3.186273e-16 -1.326579e-16 4.128321e-17 -6.023769e-17 -3.161449e-16 9.824799e-17 ... 3.482472e-16 -6.066227e-16 3.415922e-16 -4.382794e-17 3.477669e-16 2.153475e-16 -5.407545e-17 -3.659698e-16 -11.432744 0.001176
V20 -1.862195e+03 1.167581e-16 8.602176e-17 1.501888e-17 -2.057175e-16 -1.774199e-16 1.268996e-16 1.786174e-16 1.124420e-16 -2.917704e-16 ... -6.614669e-16 5.314426e-16 6.533961e-17 7.554659e-17 4.378428e-18 -1.297577e-16 -3.006146e-16 -6.298824e-17 65.445072 0.000643
V21 1.560435e+03 -1.691534e-16 8.721928e-17 -1.889959e-16 -5.786136e-17 -7.223938e-17 -9.924436e-17 3.142395e-17 3.613762e-17 1.899876e-16 ... 5.395255e-01 1.865772e-15 2.918952e-16 6.076160e-17 -4.707745e-17 -1.772897e-16 -4.174631e-16 3.679563e-17 19.474041 0.001233
V22 4.964595e+03 1.559269e-16 2.080689e-16 -2.348384e-16 2.599239e-16 4.573025e-17 -1.017080e-16 -5.758880e-16 1.493156e-17 -1.223621e-16 ... 1.865772e-15 5.266428e-01 1.313092e-16 1.332240e-17 -3.658763e-16 -8.804511e-18 4.925419e-17 -1.296150e-16 -11.762131 0.000024
V23 1.516599e+03 2.421982e-16 9.431707e-17 -7.394678e-17 1.700352e-16 9.701149e-17 3.371062e-17 -2.145430e-16 1.519976e-16 -6.412962e-17 ... 2.918952e-16 1.313092e-16 3.899507e-01 2.487034e-17 -2.341586e-16 3.876437e-16 1.088027e-16 2.850687e-16 -17.592087 -0.000070
V24 -4.654076e+02 -5.106918e-17 -1.081634e-16 3.089224e-17 1.347708e-16 -8.131900e-16 -8.535315e-16 -3.648690e-18 -1.545859e-16 -1.877797e-16 ... 6.076160e-17 1.332240e-17 2.487034e-17 3.668084e-01 3.937607e-16 5.488754e-17 -7.522538e-17 -5.594736e-17 0.779572 -0.000182
V25 -5.769855e+03 -3.270723e-16 1.314776e-16 5.925223e-17 4.773297e-16 -6.929392e-17 3.820521e-16 2.201688e-17 -1.134399e-16 1.144129e-16 ... -4.707745e-17 -3.658763e-16 -2.341586e-16 3.937607e-16 2.717308e-01 6.118510e-16 -1.328653e-16 6.097678e-17 -6.237072 0.000072
V26 -9.482254e+02 -1.284838e-16 2.006343e-16 -1.715297e-16 -2.689739e-16 2.398717e-16 -1.645278e-16 -4.844150e-16 -7.983458e-19 -7.271184e-17 ... -1.772897e-16 -8.804511e-18 3.876437e-16 5.488754e-17 6.118510e-16 2.325429e-01 -5.454791e-17 -4.690671e-17 -0.386936 0.000089
V27 -9.841860e+01 9.270791e-17 -3.456338e-16 3.409186e-16 -6.540199e-17 2.521276e-16 -6.873882e-17 -5.993831e-17 1.586213e-16 -7.746449e-17 ... -4.174631e-16 4.925419e-17 1.088027e-16 -7.522538e-17 -1.328653e-16 -5.454791e-17 1.629192e-01 4.927291e-19 2.910121 0.000295
V28 -1.475443e+02 2.291751e-16 -1.968671e-16 3.703327e-16 -4.341005e-18 -8.280967e-17 1.897350e-16 3.056168e-17 -2.243975e-16 2.986687e-16 ... 3.679563e-17 -1.296150e-16 2.850687e-16 -5.594736e-17 6.097678e-17 -4.690671e-17 4.927291e-19 1.089550e-01 0.846923 0.000131
Amount -1.258610e+05 -1.115566e+02 -2.194854e+02 -7.997555e+01 3.496456e+01 -1.333808e+02 7.197093e+01 1.229368e+02 -3.079299e+01 -1.215825e+01 ... 1.947404e+01 -1.176213e+01 -1.759209e+01 7.795722e-01 -6.237072e+00 -3.869364e-01 2.910121e+00 8.469230e-01 62560.069046 0.058496
Class -2.430072e+01 -8.243501e-03 6.260047e-03 -1.214993e-02 7.846318e-03 -5.443715e-03 -2.414579e-03 -9.619937e-03 9.857688e-04 -4.458868e-03 ... 1.232718e-03 2.426933e-05 -6.963168e-05 -1.816117e-04 7.160261e-05 8.922171e-05 2.946665e-04 1.307146e-04 0.058496 0.001725

31 rows × 31 columns

共分散行列の非対角成分は確かにほぼ0になっており、対角化されている。

In [4]:
df.isnull().sum()
Out[4]:
Time      0
V1        0
V2        0
V3        0
V4        0
V5        0
V6        0
V7        0
V8        0
V9        0
V10       0
V11       0
V12       0
V13       0
V14       0
V15       0
V16       0
V17       0
V18       0
V19       0
V20       0
V21       0
V22       0
V23       0
V24       0
V25       0
V26       0
V27       0
V28       0
Amount    0
Class     0
dtype: int64

欠損はない。
なので今回の解析では補完しなくてよい。
このデータを作るときに誰かが前処理してPCAしてるはずだが、その際補完したのかもしれない。

In [33]:
for label in df.columns:
    df.plot(kind='hist', y=label , bins=int(np.sqrt(len(df))))
    plt.show()
    plt.clf()
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>

標本のヒストグラムを片っ端から作ってみた。
見てみると多峰的だったり片方にテールを引いていたりするものもある。
現実は綺麗な正規分布になるとは限らないのである
だが、今回はあえて無視して正規分布を仮定してモデリングする。

次に異常と判定されるべきクラスを概観する

In [10]:
group_by_class = df.groupby("Class")
df_fraud = group_by_class.get_group(1)
df_fraud.describe()
Out[10]:
Time V1 V2 V3 V4 V5 V6 V7 V8 V9 ... V21 V22 V23 V24 V25 V26 V27 V28 Amount Class
count 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 ... 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.000000 492.0
mean 80746.806911 -4.771948 3.623778 -7.033281 4.542029 -3.151225 -1.397737 -5.568731 0.570636 -2.581123 ... 0.713588 0.014049 -0.040308 -0.105130 0.041449 0.051648 0.170575 0.075667 122.211321 1.0
std 47835.365138 6.783687 4.291216 7.110937 2.873318 5.372468 1.858124 7.206773 6.797831 2.500896 ... 3.869304 1.494602 1.579642 0.515577 0.797205 0.471679 1.376766 0.547291 256.683288 0.0
min 406.000000 -30.552380 -8.402154 -31.103685 -1.313275 -22.105532 -6.406267 -43.557242 -41.044261 -13.434066 ... -22.797604 -8.887017 -19.254328 -2.028024 -4.781606 -1.152671 -7.263482 -1.869290 0.000000 1.0
25% 41241.500000 -6.036063 1.188226 -8.643489 2.373050 -4.792835 -2.501511 -7.965295 -0.195336 -3.872383 ... 0.041787 -0.533764 -0.342175 -0.436809 -0.314348 -0.259416 -0.020025 -0.108868 1.000000 1.0
50% 75568.500000 -2.342497 2.717869 -5.075257 4.177147 -1.522962 -1.424616 -3.034402 0.621508 -2.208768 ... 0.592146 0.048434 -0.073135 -0.060795 0.088371 0.004321 0.394926 0.146344 9.250000 1.0
75% 128483.000000 -0.419200 4.971257 -2.276185 6.348729 0.214562 -0.413216 -0.945954 1.764879 -0.787850 ... 1.244611 0.617474 0.308378 0.285328 0.456515 0.396733 0.826029 0.381152 105.890000 1.0
max 170348.000000 2.132386 22.057729 2.250210 12.114672 11.095089 6.474115 5.802537 20.007208 3.353525 ... 27.202839 8.361985 5.466230 1.091435 2.208209 2.745261 3.052358 1.779364 2125.870000 1.0

8 rows × 31 columns

countを見ると、異常クラスの標本サイズは492個あるようだ
データ全体が284807個あるので、およそ0.2%が異常。
データ全体と比べ、meanは0に近くない。
ホテリングのT2法は標本平均からのズレをマハラノビス距離で評価するものだから
異常クラスの標本平均が標本全体の標本平均と違うことはホテリングのT2法が「全然使えないわけじゃない」ことを予感させる。
ただし、異常クラスのV1のminは標本全体のV1のminより大きい。 これはV1の最小値である-5.640751e+01の値を持つ標本は正常クラスに属していることを示している。
V1以外の特徴量を考慮することでこの点は克服できるだろうか?

訓練

0. データ分割

データを訓練用、交差検証用の2つに分割する

In [12]:
df_use = df.drop(['Time','Amount'],axis=1)
df_train, df_test= train_test_split(df_use, test_size=0.1)
target_train = df_train['Class']
target_test = df_test['Class']
df_train = df_train.drop('Class',axis=1)
df_test = df_test.drop('Class',axis=1)

1. 異常度の分布であるカイ2乗分布から異常と判定する閾値を決定する

下記の式にしたがって異常度の閾値を決定する。
$$ 1-\alpha = \int ^{a_{\mathrm{th}}} _0 dx \chi ^2 (x|M,1) $$ $\alpha$は誤報率。カイ2乗分布に従う確率変数に対して、どれぐらいの確率で起きそうなことが起きてしまったら異常であるとみなすか決めている。
$a_{\mathrm{th}}$が異常度の閾値。 $M$は自由度。 理屈はこうだが、計算しやすいように少し変形する。
カイ二乗分布の累積分布関数は以下のように定義される。
$$ X ^2 (x|k,s) = \int ^{x} _0 dt \chi ^2 (t|k,s) $$ scipy.stats.chi2で定義されている便利な関数にSurvival functionとその逆関数がある。
Survival functionは$1-X ^2 (x|k,s)$で定義される。
Survival functionの逆関数を$isf(x|k,s)$とすれば、求めるべき$a_{\mathrm{th}}$は $$ a_{\mathrm{th}} = isf(\alpha |M,1) $$ となる。

In [13]:
def a_threshold(prob,degree_of_freedom):
    anomality = scipy.stats.chi2.isf(prob, degree_of_freedom)
    return anomality

def test_a_threshold():
    prob = 0.05
    a_th = a_threshold(prob,30)
    print(a_th)

テスト用のコードを実行すると異常度の閾値が表示されるので
これを分布表と照らしあわせて正しい値が出ていることを確認できる。
例えば上記の条件だと、異常度の閾値は43.77297182574217になるが
これは https://www.biwako.shiga-u.ac.jp/sensei/mnaka/ut/chi2disttab.html を見ても
正しい値が計算できている

2. 訓練用データの標本平均、標本共分散行列を計算し、異常度を算出できるようにする

TsqDetectorというクラスを定義する。
このクラスはホテリングのT2法を実行するために必要なパラメータをメンバ変数にもち
与えられたdfに対してホテリングのT2法を実行して異常度と判定結果を返すメソッドを持つ。

実装が正しく行われていることを確認するためのテストのための関数をつけておいた。
「手計算」でホテリングのT2法が実行できる小さなデータセットに対して正しく異常度を計算できるか確認できる。

In [20]:
class TsqDetector(object):
    def __init__(self,mean,cov,prob,dof):
        self.mean = mean
        self.cov = cov
        self.dof = dof
        self.prob = prob
        self.threshold = a_threshold(prob,dof)
        
        
    def set_mean(self,df_train):
        self.mean = df_train.mean()
        
    def set_cov(self,df_train):
        self.cov = df_train.cov()
        
    def set_threshold(self,prob):
        self.threshold = a_threshold(prob,self.dof)
    
    def calc_anomality(self,df_test):
        #calc anomality for each test data
        result = []
        metric = self.cov.values
        metric = np.linalg.inv(metric)
        iter = 0
        for index, row in df_test.iterrows():
            dist = row - self.mean
            dist = np.dot(metric, dist.values.T)
            dist = np.dot(row-self.mean,dist)
            iter += 1
            result.append(dist)
            if iter % 1000 == 999:
                #print(iter+1," th iteration")
                pass
        #return list of anomality
        return result
        
    def examine(self,df_test):
        #calc anomality
        list_anomality = self.calc_anomality(df_test)
        
        #compare a_th
        list_prediction = [0 if anomality < self.threshold else 1 for anomality in list_anomality]

        #return result list
        return list_anomality, list_prediction 
    
    def re_examine(self,list_anomality,prob):
        list_prediction = [0 if anomality < a_threshold(prob,self.dof) else 1 for anomality in list_anomality]
        return list_prediction
        

#debug code
def test_TsqDetector():
    df_check=pd.DataFrame({'V1':[101,102,130,99,98],'V2':[95,105,70,103,101],'Class':[0,0,1,0,0]})
    group_by_class = df_check.groupby("Class")
    df_neg = group_by_class.get_group(0)
    Detector = TsqDetector(mean=df_neg.drop('Class',axis=1).mean(),cov=df_neg.drop('Class',axis=1).cov(),prob=0.005)
    df_check['anomality'], df_check['prediction'] = Detector.examine(df_check)
    no0_anomality = (101-100)*((10/3)**(-1))*(101-100) + (95-101)*(18.6666667**(-1))*(95-101)
    print(no0_anomality)
    print(df_check)

3. 交差検証用データで推論の妥当性を検証する

試しに誤報率を0.2%に設定して分類器を作ってみる。 これは標本全体に占める異常クラスの標本数の割合と等しい。

In [89]:
Detector = TsqDetector(mean=df_train.mean(),cov=df_train.cov(),prob=0.002,dof=len(df_train.columns))
anomality, prediction = Detector.examine(df_test)
1000  th iteration
2000  th iteration
3000  th iteration
4000  th iteration
5000  th iteration
6000  th iteration
7000  th iteration
In [17]:
confusion_matrix(target_test,prediction)

confusion_matrixのindexは(grand truth, prediction)
非対角成分が予測に失敗している標本の数になる
つまり、
左上がTN
右上がFP
左下がFN
右下がTP

4. ROC曲線の作成

ROC曲線とは、以下の式で定義されるxy平面上の点の集合がなす曲線のことである
$$ (x,y) = (1-r_0(\tau), r_1(\tau)) $$ なお、$r_0(\tau)$、$r_1(\tau)$はそれぞれ閾値$\tau$の時の正常標本精度、異常標本精度で、以下の式で定義される
$$ r_0(\tau) = \frac{TN(\tau)}{TN(\tau)+FP(\tau)} \\ r_1(\tau) = \frac{TP(\tau)}{TP(\tau)+FN(\tau)} $$

In [22]:
Detector = TsqDetector(mean=df_train.mean(),cov=df_train.cov(),prob=0.002,dof=len(df_train.columns))
list_anomality = Detector.calc_anomality(df_test)
In [27]:
def make_ROC_curve(tau_list, list_anomality, Detector):
    x = np.array([1])
    y = np.array([1])
    CMs = np.array([])
    
    for tau in tau_list:
        prediction = Detector.re_examine(list_anomality, tau)
        CM = confusion_matrix(target_test,prediction)
        r0 = CM[0][0]/(CM[0][0]+CM[0][1])
        r1 = CM[1][1]/(CM[1][0]+CM[1][1])
        x=np.append(x,1-r0)
        y=np.append(y,r1)
        CMs=np.append(CMs,CM)
        
    x=np.append(x,0)
    y=np.append(y,0)    
    
    plt.plot(x,y,marker='.')
    plt.xlim(0,1)
    plt.ylim(0,1)
    plt.title("ROC curve")
    plt.xlabel("false positive rate")
    plt.ylabel("recall")
    plt.show()
    
    return x,y,CMs

tau_list = [0.99,0.5,0.3,0.1,0.03,0.01,0.001,0.0001]
x,y,CMs=make_ROC_curve(tau_list, list_anomality, Detector)

ROC曲線が得られたので、AUCを計算してみる。
AUCはROC曲線の下部面積で定義される。
ROC曲線は区分的に直線だと近似して台形公式で面積を求める。

In [28]:
def AUC_calc(x,y):
    #make sure input x, y are sorted
    curve=np.vstack([x,y])
    curve=curve[:,np.argsort(curve[0])]
    
    area = 0
    for i in range(len(x)-1):
        area += 1/2*(curve[0][i+1]-curve[0][i])*(curve[1][i]+curve[1][i+1])
    return area

print("AUC=",AUC_calc(x,y))
AUC= 0.943841317690511

top kernelのAUCが0.95なのでかなり良く出来ているのでしょう。
https://www.kaggle.com/joparga3/in-depth-skewed-data-classif-93-recall-acc-now
でもなーFPがTPの10倍出てるわけで、この分類器使うと鬱陶しいかもしれない。 詐欺検出なんだから、疑わしきは調査というスタンスならこういう運用もあり。
実用上は調査にかかるコストと検出漏れの際の損失を天秤にかけてバランスとるんでしょう。
つまり、1回あたりの調査コストを$f$、1回あたりの検出漏れ損失額を$g$として $$ argmin_{\tau}\left[ \frac{FP(\tau)+TP(\tau)}{TP(\tau)+TN(\tau)+FP(\tau)+FN(\tau)}f+\frac{FN(\tau)}{TP(\tau)+TN(\tau)+FP(\tau)+FN(\tau)}g \right] $$ を満たす$\tau$を取ってくればいいでしょう。

異常度のヒストグラム

交差検証用のデータセットについて異常度が計算できたので、異常度の分布を調べてみましょう。

In [74]:
hist_max=100

df_result = pd.DataFrame({'Grand truth':target_test,"Anomality":list_anomality})
df_result.plot(kind='hist', y="Anomality" , bins=int(np.sqrt(len(df_result))), 
               range=(0,hist_max), density=True)

chi2_x = range(0,hist_max,1)
chi2_y = [scipy.stats.chi2.pdf(x=x,df=28) for x in chi2_x] 
plt.plot(chi2_x,chi2_y,label="chi2 dof=28")

chi2_x = range(0,hist_max,1)
chi2_y = [scipy.stats.chi2.pdf(x=x,df=15) for x in chi2_x] 
plt.plot(chi2_x,chi2_y,label="chi2 dof=15")

plt.legend()
plt.show()

交差検証用データの異常度の分布はこんな感じ。 カイ二乗分布に従う想定をしたので、カイ二乗分布を同時にプロットした。
自由度28のカイ二乗分布から想定されるヒストグラムよりかなりピークが左にずれてしまっている。
書籍によれば、データの有効次元が理論的に想定される値よりかなり小さくなるのはよくあることらしい。
実際、このデータに対して例えば自由度15のカイ二乗分布をあてはめてみると、多少マシになる(それでもぴったりは合わないが・・・)
最適なパラメータの計算方法は書籍7.4章を参照のこと(そのうち実装するかも)

次に異常クラスの分布を見る。

In [95]:
df_result_groupby = df_result.groupby("Grand truth")
df_result_fraud = df_result_groupby.get_group(1)
df_result_fraud.plot(kind='hist', y="Anomality" , bins=int(np.sqrt(len(df_result_fraud))),
                     density=True, title="histgram of anomalies")
                     
df_result_fraud.plot(kind='hist', y="Anomality" , bins=int(np.sqrt(len(df_result_fraud))),
                     density=True, title="zoom up of histgram of anomalies", range=(0,100))

df_result_fraud.plot(kind='hist', y="Anomality" , bins=5000, 
               cumulative=True,histtype='step',density=True, title="cumulative ratio")

df_result_fraud.plot(kind='hist', y="Anomality" , bins=5000, 
               cumulative=True,histtype='step',density=True,range=(0,100),title="zoom up of cumulative ratio")
Out[95]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7558776860>

この累積分布に対して1つ異常度の閾値を決めると、その点より右側は検出できる異常、左側は検出できない異常になります。 大きな異常度のところにある程度集まっているのは想定通りですね。

ROC曲線の平坦な部分は累積分布の20あたりから70あたりの平坦な部分に閾値が来た時に対応していますね。
a_threshold(10e-4,29)は58ぐらいなので、精度を落とさずにもう少し閾値を上げられるように思いますが
安直に実行するのはこのデータセットだけに過学習する行為です。
もし試すなら、検証用データを何度も取り替えてこの平坦な領域が出てくるか確認すべきでしょう。

閾値を20以下にするとすべての異常を拾えますが、FPの割合が増えてしまいます。
これがこのモデルの実力です。
もしFPを増やさずに異常を拾う割合を増やしたければ、モデル自体を変更する必要があるでしょう。

In [99]:
a_threshold(10e-4,29)
Out[99]:
58.30117348979492

考察

df_trainからモデルを作る際、正常クラスに属するデータが圧倒的に多いことが必要だった。 モデルを作る際は正常クラスを表すモデルを作ればいいので、正常クラスのみを使ってモデルを作ると予測精度が上がると思われる。

・こんなシンプルなモデルでもそれなりにいい成績が出る。PCAが偉大なんだと思う。  

・変数間の相関を無視して、各変数は統計的に独立だとすると精度がわるくなるだろうか?
試しに計算してみる。

TsqDetectorを継承してSimpleTsqDetectorクラスを定義する。
ほとんど親クラスと同じだが、異常度の計算の実装が異なっている。
マハラノビス距離の計算の前に共分散行列の非対角項を落としている。
(ここ、多分df.var()で最初から分散のベクトル持ってきた方が計算速いと思う)
ROC_curve_simpleSimpleTsqDetectorを呼び出すように別の関数を定義   ROC_curve_simpleの柔軟性、低いなあ・・・。

追記:曲線のデータ作成機能と異常度の計算機能を分離できることに気づいたので、同じmake_ROC_curveを使えるようになりました

In [26]:
a=np.array([[1,2],[3,4]])
np.diag(np.diag(a))

class SimpleTsqDetector(TsqDetector):
    def __init__(self,mean,cov,prob,dof):
        super().__init__(mean,cov,prob,dof)
    
    def calc_anomality(self,df_test):
        #calc anomality for each test data
        result = []
        metric = self.cov.values
        metric = np.diag(np.diag(metric))
        metric = np.linalg.inv(metric)
        iter = 0
        for index, row in df_test.iterrows():
            dist = row - self.mean
            dist = np.dot(metric, dist.values)
            dist = np.dot(row-self.mean,dist)
            iter += 1
            result.append(dist)
            if iter % 1000 == 999:
                #print(iter+1," th iteration")
                pass
        #return list of anomality
        return result

tau_list = [0.99,0.5,0.3,0.1,0.03,0.01,0.001]
Detector = SimpleTsqDetector(mean=df_train.mean(),cov=df_train.cov(),
                                     prob=tau_list[0],dof=len(df_train.columns))
list_anomality = Detector.calc_anomality(df_test)
x,y,CMs=make_ROC_curve(tau_list,list_anomality, Detector)
print("AUC=",AUC_calc(x,y))
AUC= 0.9395565768488205

AUCはほとんど変わらず。
PCAは偉大である。
この手の異常検知タスクはまずPCAをしてホテリングのT2法できないか考えてみようと思えました。