t-SNEとPCAの実行

想定として、ある時間に200~800nmまでの波長データで何かしらの値を取得できるセンサーがあり、そのログデータから異常値を予測したいものとします。

正常データとして使用するデータと異常データとして使用するデータを読み込み
t-SNEによる次元削減を行って可視化してみます。

これでデータの分布に違いがあるか判り、使用する正常データと異常データの選択の指標の一つになります。

データは訓練用データ、検証用(異常値)データ、テスト用(異常値)データ、検証用(異常値)データです。
データの内容 - ./data/train:学習に使用するノイズの少ないsin波データ - ./data/test/normal:検証とテスト兼用の正常sin波データ - ./data/test/anomaly/test:異常値のテストデータ - ./data/test/anomaly/valid:異常値の検証用データ

正常データを1、異常データを-1としてラベルを付与して可視化します。

import os
import glob
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
plt.style.use('seaborn')

# 警告の非表示
import warnings
warnings.filterwarnings('ignore')

分析用疑似データ作成

# 保存用フォルダ作成
if not(os.path.isdir('./data')):
    os.mkdir('./data')
if not(os.path.isdir('./data/train')):
    os.mkdir('./data/train')

if not(os.path.isdir('./data/test/normal')):
    os.mkdir('./data/test/normal')
if not(os.path.isdir('./data/test')):
    os.mkdir('./data/test')

if not(os.path.isdir('data/test/anomaly')):
    os.mkdir('data/test/anomaly')
if not(os.path.isdir('data/test/anomaly/test')):
    os.mkdir('data/test/anomaly/test')
if not(os.path.isdir('data/test/anomaly/valid')):
    os.mkdir('data/test/anomaly/valid')

データ内容の確認

それぞれどのようなデータが入るかを可視化して確認してみます。

def Make_data(index, eta=0.2, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 40個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 40 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    plt.subplots(figsize=(16, 9))                                            # 表示サイズ指定
    plt.scatter(X, y)                                                        # 散布図
    plt.show()
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    
    return df
# train
Make_data(1, 0.05).head()

f:id:sekihan_0290:20200918215431p:plain

1
0 -0.012546
1 0.050312
2 0.033680
3 0.025586
4 -0.013438
# valid_anomaly
Make_data(1, 0.15, outlier=2.5, mode='test').head()

f:id:sekihan_0290:20200918215443p:plain

1
0 -0.037638
1 0.140455
2 0.080079
3 0.045318
4 -0.082235
# test_anomaly
Make_data(1, 0.2, outlier=3, mode='test').head()

f:id:sekihan_0290:20200918215511p:plain

1
0 -0.050184
1 0.185526
2 0.103278
3 0.055184
4 -0.116633
# 可視化の内容を消去して再度宣言
def Make_data(index, eta=0.05, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 100個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 100 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    return df

訓練用データを作成

関数Make_data()を使用して訓練用データを作成します。
ノイズはデフォルトの0.2でデータ数160個作成します。

df_train = pd.DataFrame([])
for i in range(1, 161):
    df_base =  Make_data(i).T
    df_train = pd.concat([df_train, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_train.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_train.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

訓練データの保存

# index不要のとき
df_train.to_csv('data/train/X_train.csv', index=False)

検証とテスト兼用の正常sin波データ

ノイズはデフォルトの0.2でデータ数80個作成します。

df_normal = pd.DataFrame([])
for i in range(1, 81):
    df_base =  Make_data(i).T
    df_normal = pd.concat([df_normal, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_normal.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_normal.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

検証とテスト兼用の正常sin波データの保存

# index不要のとき
df_normal.to_csv('./data/test/normal/X_normal.csv', index=False)

検証用(異常値)データを作成

ノイズは0.15で、外れ値は2.5、データ数30個作成します。

df_valid = pd.DataFrame([])
for i in range(1, 31):
    df_base =  Make_data(i, 0.15, outlier=2.5, mode='test').T
    df_valid = pd.concat([df_valid, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_valid.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_valid.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
2 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
3 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
4 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
5 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974

5 rows × 1200 columns

検証(異常値)データの保存

# index不要のとき
df_valid.to_csv('./data/test/anomaly/valid/X_valid.csv', index=False)

テスト用(異常値)データを作成

ノイズはデフォルトの0.2で外れ値は3、データ数48個作成します。

df_test = pd.DataFrame([])
for i in range(1, 49):
    df_base =  Make_data(i, 0.2, outlier=3, mode='test').T
    df_test = pd.concat([df_test, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_test.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_test.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
2 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
3 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
4 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
5 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966

5 rows × 1200 columns

テスト用(異常値)データの保存

# index不要のとき
df_test.to_csv('./data/test/anomaly/test/X_test.csv', index=False)

データの読み込みとラベル付与

start = time.time()             # 実行開始時間の取得
"""
正常データの読み込みとラベル付与
"""
# フォルダ内のファイル一覧を取得
check_train_normal_files = glob.glob('./data/train/*')
check_test_normal_files = glob.glob('./data/test/normal/*')

X_train = pd.DataFrame([])
y_train = []
for file_name in check_train_normal_files:
    csv = pd.read_csv(file_name)
    X_train = pd.concat([X_train, csv])

    for i in range(0, len(csv)):
        y_train.append(1)

for file_name in check_test_normal_files:
    csv = pd.read_csv(file_name)
    X_train = pd.concat([X_train, csv])

    for i in range(0, len(csv)):
        y_train.append(1)
"""
異常データの読み込みとラベル付与
"""
check_anomaly1_files = glob.glob('./data/test/anomaly/test/*')
check_anomaly2_files = glob.glob('./data/test/anomaly/valid/*')

X_test = pd.DataFrame([])
y_test = []

for file_name in check_anomaly1_files:
    csv = pd.read_csv(file_name)
    X_test = pd.concat([X_test, csv])

    for i in range(0, len(csv)):
        y_test.append(-1)

for file_name in check_anomaly1_files:
    csv = pd.read_csv(file_name)
    X_test = pd.concat([X_test, csv])

    for i in range(0, len(csv)):
        y_test.append(-1)
# trainとtestのデータを一つにまとめる
df = pd.concat([X_train, X_test])
targets = y_train
targets.extend(y_test)
targets = np.array(targets)

t-SNEの実行

データ全体を使用して学習を実行してみます。

tsne = TSNE(n_components=2, random_state=42)
df_embedded = tsne.fit_transform(df)
# xとyはt-SNE分解からの2つのコンポーネントで、targetsは実際の数
tsne_df = pd.DataFrame(
    # データを縦(カラム方向)に結合、学習結果と説明変数を結合
    np.column_stack((df_embedded, targets)),
    columns=['x', 'y', 'targets']
)
# 念のため正解ラベルのint型への型変換
tsne_df.loc[:, 'targets'] = tsne_df['targets'].astype(int)
grid = sns.FacetGrid(tsne_df, hue='targets', size=8)
grid.map(plt.scatter, 'x', 'y')
# grid.set(xlim=(-20, 20), ylim=(-60, 60), xticks=[-20, -10, 0, 10, 20])
grid.add_legend()
plt.title('t-SNE scatter prot')
plt.tight_layout()
# plt.savefig('./image/t-SNE_data_check.png')
plt.show()

f:id:sekihan_0290:20200918215522p:plain

PCAの実行

データ全体を使用して学習を実行してみます。

pca = PCA(n_components=2, random_state=42)
df_embedded = pca.fit_transform(df)
# xとyはt-SNE分解からの2つのコンポーネントで、targetsは実際の数
tsne_df = pd.DataFrame(
    # データを縦(カラム方向)に結合、学習結果と説明変数を結合
    np.column_stack((df_embedded, targets)),
    columns=['x', 'y', 'targets']
)
# 念のため正解ラベルのint型への型変換
tsne_df.loc[:, 'targets'] = tsne_df['targets'].astype(int)
grid = sns.FacetGrid(tsne_df, hue='targets', size=8)
grid.map(plt.scatter, 'x', 'y')
# grid.set(xlim=(-20, 20), ylim=(-60, 60), xticks=[-20, -10, 0, 10, 20])

# seabornの凡例を使用すると結果と被るので変更
# grid.add_legend()
plt.legend(loc="best", frameon=True, edgecolor="blue")
plt.title('PCA scatter prot')
plt.tight_layout()
plt.savefig('./image/PCA_data_check.png')
plt.show()

f:id:sekihan_0290:20200918215536p:plain

普段は細かくデータの選別をするならt-SNE
よりクラスタ間の距離感をみるならPCAを重視して確認を行っています。

カラムごとのOne-class-SVM(OCSVM)の実行

二次元データで外れ値検知を実行している例はよく見ますが、複数列を持つデータフレームに対し、各列で学習と予測を行い異常カラムを検出するみたいな例は見なかったので実践に使う前のプロト版として作ってみたものを公開してみます。(あまりいい例じゃないですが)

想定として、ある時間に200~800nmまでの波長データで何かしらの値を取得できるセンサーがあり、そのログデータから異常値を予測したいものとします。

具体的には、ノイズを加えたsin波を例に異常検知を実装します。
カラムは、200~800まで0.5刻みで記録されています。
100%分類できてるみたいで面白くないけどあくまで例ということで
データは訓練用データ、検証用(異常値)データ、テスト用(異常値)データ、検証用(異常値)データです。
データの内容 - ./data/train:学習に使用するノイズの少ないsin波データ - ./data/test/normal:検証とテスト兼用の正常sin波データ - ./data/test/anomaly/test:異常値のテストデータ - ./data/test/anomaly/valid:異常値の検証用データ

まずデータ全体を見て異常検知を行ったあと、
カラムごとに異常検知を行い、カラムごとに複数回異常と判定されたカラムとその異常スコアの最大値をDataFrameとして出力してみます。

import os
import glob
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('seaborn')
import itertools
from sklearn.svm import OneClassSVM
from sklearn.metrics import precision_recall_fscore_support,confusion_matrix
from sklearn. metrics import roc_curve, roc_auc_score
from sklearn.metrics import precision_recall_fscore_support

import warnings
warnings.filterwarnings('ignore')

分析用疑似データ作成

# 保存用フォルダ作成
if not(os.path.isdir('./data')):
    os.mkdir('./data')
if not(os.path.isdir('./data/train')):
    os.mkdir('./data/train')

if not(os.path.isdir('./data/test/normal')):
    os.mkdir('./data/test/normal')
if not(os.path.isdir('./data/test')):
    os.mkdir('./data/test')

if not(os.path.isdir('data/test/anomaly')):
    os.mkdir('data/test/anomaly')
if not(os.path.isdir('data/test/anomaly/test')):
    os.mkdir('data/test/anomaly/test')
if not(os.path.isdir('data/test/anomaly/valid')):
    os.mkdir('data/test/anomaly/valid')

データ内容の確認

それぞれどのようなデータが入るかを可視化して確認してみます。

def Make_data(index, eta=0.2, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 40個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 40 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    plt.subplots(figsize=(16, 9))                                            # 表示サイズ指定
    plt.scatter(X, y)                                                        # 散布図
    plt.show()
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    
    return df
# train
Make_data(1, 0.05).head()

f:id:sekihan_0290:20200917221439p:plain

1
0 -0.012546
1 0.050312
2 0.033680
3 0.025586
4 -0.013438
# valid_anomaly
Make_data(1, 0.15, outlier=2.5, mode='test').head()

f:id:sekihan_0290:20200917221451p:plain

1
0 -0.037638
1 0.140455
2 0.080079
3 0.045318
4 -0.082235
# test_anomaly
Make_data(1, 0.2, outlier=3, mode='test').head()

f:id:sekihan_0290:20200917221501p:plain

1
0 -0.050184
1 0.185526
2 0.103278
3 0.055184
4 -0.116633
# 可視化の内容を消去して再度宣言
def Make_data(index, eta=0.05, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 100個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 100 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    return df

訓練用データを作成

関数Make_data()を使用して訓練用データを作成します。
ノイズはデフォルトの0.2でデータ数160個作成します。

df_train = pd.DataFrame([])
for i in range(1, 161):
    df_base =  Make_data(i).T
    df_train = pd.concat([df_train, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_train.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_train.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

訓練データの保存

# index不要のとき
df_train.to_csv('data/train/X_train.csv', index=False)

検証とテスト兼用の正常sin波データ

ノイズはデフォルトの0.2でデータ数80個作成します。

df_normal = pd.DataFrame([])
for i in range(1, 81):
    df_base =  Make_data(i).T
    df_normal = pd.concat([df_normal, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_normal.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_normal.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

検証とテスト兼用の正常sin波データの保存

# index不要のとき
df_normal.to_csv('./data/test/normal/X_normal.csv', index=False)

検証用(異常値)データを作成

ノイズは0.15で、外れ値は2.5、データ数30個作成します。

df_valid = pd.DataFrame([])
for i in range(1, 31):
    df_base =  Make_data(i, 0.15, outlier=2.5, mode='test').T
    df_valid = pd.concat([df_valid, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_valid.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_valid.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
2 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
3 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
4 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
5 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974

5 rows × 1200 columns

検証(異常値)データの保存

# index不要のとき
df_valid.to_csv('./data/test/anomaly/valid/X_valid.csv', index=False)

テスト用(異常値)データを作成

ノイズはデフォルトの0.2で外れ値は3、データ数48個作成します。

df_test = pd.DataFrame([])
for i in range(1, 49):
    df_base =  Make_data(i, 0.2, outlier=3, mode='test').T
    df_test = pd.concat([df_test, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_test.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_test.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
2 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
3 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
4 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
5 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966

5 rows × 1200 columns

テスト用(異常値)データの保存

# index不要のとき
df_test.to_csv('./data/test/anomaly/test/X_test.csv', index=False)

データの読み込みとラベル付与

start = time.time()             # 実行開始時間の取得
# 学習用フォルダ内のファイル一覧を取得
files_train = glob.glob('./data/train/*')
# 検証とテスト兼用の正常sin波データァイル一覧
files_normal  = glob.glob('./data/test/normal/*')
# 異常値のテストデータ
files_anomaly1 = glob.glob('./data/test/anomaly/test/*')
# 異常値の検証用データ
files_anomaly2 = glob.glob('./data/test/anomaly/valid/*')

# CSV形式のファイルを読み込み、学習データを全てx_trainに格納
X_train = pd.DataFrame([])
for file_name in files_train:
    csv = pd.read_csv(filepath_or_buffer=file_name)
    X_train = pd.concat([X_train, csv])
# StandardScalerでz標準化
# sc = StandardScaler()
# sc.fit(X_train)
# X_train = sc.transform(X_train)

正解ラベルの付与

# 正常データのラベルを1,異常データのラベルを-1として
# テストデータ用y_test_true,検証用データ用y_valid_trueに格納
x_test_normal, x_test_anomaly1, x_test_anomaly2, X_test, X_valid = pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([])
y_test_true, y_valid_true = [], []

# 検証とテスト兼用の正常sin波データを読み込み
# 正常データのラベル1を付与
for file_name in files_normal:
    csv = pd.read_csv(file_name)
    x_test_normal = pd.concat([x_test_normal, csv])
    for i in range(0, len(csv)):
        y_test_true.append(1)
        y_valid_true.append(1)
# 異常値の検証用データを読み込み
# 異常データのラベル-1を付与
for file_name in files_anomaly1:
    csv = pd.read_csv(file_name)
    x_test_anomaly1 = pd.concat([x_test_anomaly1, csv])
    for i in range(0, len(csv)):
        y_test_true.append(-1)
# 異常値のテスト用データを読み込み
# 異常データのラベル-1を付与
for file_name in files_anomaly2:
    csv = pd.read_csv(file_name)
    x_test_anomaly2 = pd.concat([x_test_anomaly2, csv])
    for i in range(0, len(csv)):
        y_valid_true.append(-1)
        
# テストデータx_test,検証用データx_validを正常データと異常データを組み合わせて用意
X_test = pd.concat([x_test_normal, x_test_anomaly1])
X_valid = pd.concat([x_test_normal, x_test_anomaly2])
# z標準化
# X_test = sc.transform(X_test)
# X_valid = sc.transform(X_valid)
# 正常データ数,異常データ数(テストデータ),テストデータ総数,検証用データ総数を確認
print('data size')
print('訓練データ:{}'.format(X_train.shape[0]))
print('テスト・検証用正常データ:{}'.format(x_test_normal.shape[0]))
print('異常データ数(テストデータ):{}'.format(x_test_anomaly1.shape[0]))
print('異常データ数(検証データ):{}'.format(x_test_anomaly2.shape[0]))
print('テストデータ総数:{}'.format(X_test.shape[0]))
print('検証用データ総数:{}'.format(X_valid.shape[0]))
data size
訓練データ:160
テスト・検証用正常データ:80
異常データ数(テストデータ):48
異常データ数(検証データ):30
テストデータ総数:128
検証用データ総数:110

One-class-SVM(OCSVM)の実行

データ全体を使用して学習を実行してみます。

# ハイパーパラメータチューニングするリスト
gamma = [0.001, 0.005, 0.01]                # 係数
coef0 = [0.1, 1.0, 5.0]                     # 定数項
degree = [1, 2, 3]                          # 次数

# RBFカーネルのバンド幅パラメータγを変化させたときの検証用データに対するF値の取得
idx, f_score = [], []
for r in gamma:
    ocsvm_rbf = OneClassSVM(kernel='rbf', gamma=r)
    ocsvm_rbf.fit(X_train)
    prec_rec_f = precision_recall_fscore_support(y_valid_true, ocsvm_rbf.predict(X_valid))
    f_score.append(np.average(prec_rec_f[2]))
    idx.append(r)
# F値が最大になるバンド幅γを取得し、再度RBFカーネルで学習
plt.plot(idx, f_score, 'b-')
plt.xlabel('gamma')
plt.ylabel('F-score')
plt.title('ALL_Spectrim_F-score')
# plt.savefig('./alldata_image/ALL_Spectrim_F-score.png')

best_rbf_gamma = gamma[np.argmax(f_score)]
print('RBF best_kernel: gamma:{:2.4}, f-score:{:.4}'.format(best_rbf_gamma, np.max(f_score)))
ocsvm_rbf = OneClassSVM(kernel='rbf', gamma=best_rbf_gamma)
ocsvm_rbf.fit(X_train)

# 多項式カーネルパラメータ(次数d, 係数γ, 定数項c)を変化させたときの検証用データに対するF値を取得
idx, f_score = [], []
for d, r, c in itertools.product(degree, gamma, coef0):
    ocscvm_poly = OneClassSVM(kernel='poly', degree=d, gamma=r, coef0=c)
    ocscvm_poly.fit(X_train)
    prec_rec_f = precision_recall_fscore_support(y_valid_true, ocscvm_poly.predict(X_valid))
    f_score.append(np.average(prec_rec_f[2]))
    idx.append([d, r, c])
# F値が最大になるパラメーターの組み合わせを取得し、再度多項式カーネルで学習
best_idx = idx[np.argmax(f_score)]
print('Polynomial best_kernel: degree:{:1}, gamma:{:.4}, coef0:{:3.2}, f_score:{:.4}'.format(best_idx[0], best_idx[1], best_idx[2], np.max(f_score)))
ocscvm_poly = OneClassSVM(kernel='poly', degree=best_idx[0], gamma=best_idx[1], coef0=best_idx[2])
ocscvm_poly.fit(X_train)
RBF best_kernel: gamma:0.001, f-score:0.2143
Polynomial best_kernel: degree:1, gamma:0.001, coef0:1.0, f_score:0.4211





OneClassSVM(coef0=1.0, degree=1, gamma=0.001, kernel='poly')

[f:id:sekihan_0290:20200917221515p:plain]

# 多項式カーネルパラメータ(次数d, 係数γ, 定数項c)を変化させたときの検証用データに対するF値を取得
idx, f_score = [], []
for d, r, c in itertools.product(degree, gamma, coef0):
    ocscvm_poly = OneClassSVM(kernel='poly', degree=d, gamma=r, coef0=c)
    ocscvm_poly.fit(X_train)
    prec_rec_f = precision_recall_fscore_support(y_valid_true, ocscvm_poly.predict(X_valid))
    f_score.append(np.average(prec_rec_f[2])) 
    idx.append([d, r, c])

# F値が最大になるパラメーターの組み合わせを取得し、再度多項式カーネルで学習
best_idx = idx[np.argmax(f_score)]
print('Polynomial best_kernel: degree:{:1}, gamma:{:.4}, coef0:{:3.2}, f_score:{:.4}'.format(best_idx[0], best_idx[1], best_idx[2], np.max(f_score)))
ocscvm_poly = OneClassSVM(kernel='poly', degree=best_idx[0], gamma=best_idx[1], coef0=best_idx[2])
ocscvm_poly.fit(X_train)
Polynomial best_kernel: degree:1, gamma:0.001, coef0:1.0, f_score:0.4211





OneClassSVM(coef0=1.0, degree=1, gamma=0.001, kernel='poly')
# シグモイドカーネル係数γを変化させたときの検証用データに対するF値を取得
idx, f_score = [], []
for r, c in itertools.product(gamma, coef0):
    ocsvm_smd = OneClassSVM(kernel='sigmoid', gamma=r, coef0=c)
    ocsvm_smd.fit(X_train)
    prec_rec_f = precision_recall_fscore_support(y_valid_true, ocscvm_poly.predict(X_valid))
    f_score.append(np.average(prec_rec_f[2]))
    idx.append([r, c])
# F値が最大になる係数γを取得し、再度シグモイドカーネルで学習
best_idx = idx[np.argmax(f_score)]
print('Sigmoid best_kernel: gamma:{:.4}, coef0:{:2.2}, f_score:{:.4}'.format(best_idx[0], best_idx[1], np.max(f_score)))
ocscvm_smd = OneClassSVM(kernel='sigmoid', gamma=best_idx[0], coef0=best_idx[1])
ocscvm_smd.fit(X_train)
Sigmoid best_kernel: gamma:0.001, coef0:0.1, f_score:0.4211





OneClassSVM(coef0=0.1, gamma=0.001, kernel='sigmoid')
print('One-class SVM result')
One-class SVM result
# 最適なパラメータを使用したRBFカーネル、多項式カーネル、シグモイドカーネルのテストデータに対する結果表示
pred_rbf = ocsvm_rbf.predict(X_test)
pred_poly = ocscvm_poly.predict(X_test)
pred_smd = ocscvm_smd.predict(X_test)
print('rbf kernel')
# 平均精度・再現率・F値と混同行列の表示
prec_rec_f = precision_recall_fscore_support(y_test_true, pred_rbf)
print('Ave. Precision {:0.4f}, Ave. Recall {:0.4f}, Ave. F-score {:0.4f}'.format(np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
print('Confusion Matrix')
df = pd.DataFrame(confusion_matrix(y_test_true, pred_rbf))
df.columns = [u'anomaly', u'normal']
df
rbf kernel
Ave. Precision 0.1875, Ave. Recall 0.5000, Ave. F-score 0.2727
Confusion Matrix
anomaly normal
0 48 0
1 80 0
print('polynomial kernel')
# 平均精度・再現率・F値と混同行列の表示
prec_rec_f = precision_recall_fscore_support(y_test_true, pred_poly)
print('Ave. Precision {:0.4f}, Ave. Recall {:0.4f}, Ave. F-score {:0.4f}'.format(np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
print('Confusion Matrix')
df = pd.DataFrame(confusion_matrix(y_test_true, pred_poly))
df.columns = [u'anomaly', u'normal']
df
polynomial kernel
Ave. Precision 0.3125, Ave. Recall 0.5000, Ave. F-score 0.3846
Confusion Matrix
anomaly normal
0 0 48
1 0 80
print('sigmoid kernel')
# 平均精度・再現率・F値と混同行列の表示
prec_rec_f = precision_recall_fscore_support(y_test_true, pred_rbf)
print('Ave. Precision {:0.4f}, Ave. Recall {:0.4f}, Ave. F-score {:0.4f}'.format(np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
print('Confusion Matrix')
df = pd.DataFrame(confusion_matrix(y_test_true, pred_smd))
df.columns = [u'anomaly', u'normal']
df
sigmoid kernel
Ave. Precision 0.1875, Ave. Recall 0.5000, Ave. F-score 0.2727
Confusion Matrix
anomaly normal
0 0 48
1 80 0
# ファイル保存ディレクトリの作成
# if not(os.path.isdir('alldata_image')):
#     os.mkdir('alldata_image')
#  正解ラベル(y_test_true)と識別関数の出力値(decision_function)を受け取って,ROC曲線を描画(SVMのカーネル比較用)
fpr, tpr, thresholds = roc_curve(y_test_true, pred_rbf, pos_label=1)
roc_auc = roc_auc_score(y_test_true, pred_rbf)
plt.plot(fpr, tpr, 'k--',label='rbf kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="r")

fpr, tpr, thresholds = roc_curve(y_test_true, pred_poly, pos_label=1)
roc_auc = roc_auc_score(y_test_true, pred_poly)
plt.plot(fpr, tpr, 'k--',label='poly kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="g")

fpr, tpr, thresholds = roc_curve(y_test_true, pred_smd, pos_label=1)
roc_auc = roc_auc_score(y_test_true, pred_smd)
plt.plot(fpr, tpr, 'k--',label='sigmoid kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="b")

plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ALL_Spectrim_ROC_curve')
plt.legend(loc='lower right')
# plt.savefig('./alldata_image/ALL_Spectrim_ROC_curve.png')
plt.show()

f:id:sekihan_0290:20200917221524p:plain

列ごとに学習を行い、異常発生時のカラムを検出する

カラムごとに学習を行い、テストデータで予測を行って42データ以上異常が予測されたものを異常値として抽出します。

ROC曲線、混同行列の表示は冗長なのでコメントアウトしてます。

# 異常と判断された波長一覧とその異常値スコア最大値、AUCの保存先
anomaly_spectrim = []
anomaly_score = []
anomaly_roc_auc_score = []
# ファイル保存ディレクトリの作成
# if not(os.path.isdir('image')):
#     os.mkdir('image')
# 平均精度,平均再現率,平均F値,ならびに混同行列を表示する関数
def print_precision_recall_fscore(y_true, y_pred):
    prec_rec_f = precision_recall_fscore_support(y_true, y_pred)
    print('Ave. Precision {:.4}, Ave. Recall {:.4}, Ave. F-score {:.4}'.format(np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
    print('Confusion Matrix')
    df = pd.DataFrame(confusion_matrix(y_true, y_pred))
    df.columns = [u'anomaly', u'normal']
    print(df)
for column_name, item in X_train.iteritems():
        # ハイパーパラメータチューニングするリスト
        gamma = [0.001, 0.005, 0.01]                # 係数
        coef0 = [0.1, 1.0, 5.0]                     # 定数項
        degree = [1, 2, 3]                          # 次数
        # RBFカーネルのバンド幅パラメータγを変化させたときの検証用データに対するF値の取得
        idx, f_score = [], []
        for r in gamma:
            ocsvm_rbf = OneClassSVM(kernel='rbf', gamma=r)
            ocsvm_rbf.fit(item.values.reshape(-1, 1))
            prec_rec_f = precision_recall_fscore_support(y_valid_true, ocsvm_rbf.predict(X_valid[column_name].values.reshape(-1, 1)))
            f_score.append(np.average(prec_rec_f[2]))
            idx.append(r)

        # F値が最大になるバンド幅γを取得し、再度RBFカーネルで学習
        # plt.plot(gamma, f_score)
        # plt.xlabel('gamma')
        # plt.ylabel('F-score')
        # plt.title('{}Spectrim_OCSVM_F-score'.format(column_name))
        # plt.savefig('./image/{}Spectrim_OCSVM_F-score.png'.format(column_name))
        # plt.show()
        
        best_rbf_gamma = gamma[np.argmax(f_score)]
        # print('RBF best_kernel: gamma:{:2.4}, f-score:{:.4}'.format(best_rbf_gamma, np.max(f_score)))
        ocsvm_rbf = OneClassSVM(kernel='rbf', gamma=best_rbf_gamma)
        ocsvm_rbf.fit(item.values.reshape(-1, 1))

        # 多項式カーネルパラメータ(次数d, 係数γ, 定数項c)を変化させたときの検証用データに対するF値を取得
        idx, f_score = [], []
        for d, r, c in itertools.product(degree, gamma, coef0):
            ocscvm_poly = OneClassSVM(kernel='poly', degree=d, gamma=r, coef0=c)
            ocscvm_poly.fit(item.values.reshape(-1, 1))
            prec_rec_f = precision_recall_fscore_support(y_valid_true, ocscvm_poly.predict(X_valid[column_name].values.reshape(-1, 1)))
            f_score.append(np.average(prec_rec_f[2]))
            
            idx.append([d, r, c])

        # F値が最大になるパラメーターの組み合わせを取得し、再度多項式カーネルで学習
        best_idx = idx[np.argmax(f_score)]
        # print('Polynomial best_kernel: degree:{:1}, gamma:{:.4}, coef0:{:3.2}, f_score:{:.4}'.format(best_idx[0], best_idx[1], best_idx[2], np.max(f_score)))
        ocscvm_poly = OneClassSVM(kernel='poly', degree=best_idx[0], gamma=best_idx[1], coef0=best_idx[2])
        ocscvm_poly.fit(item.values.reshape(-1, 1))

        # シグモイドカーネル係数γを変化させたときの検証用データに対するF値を取得
        idx, f_score = [], []
        for r, c in itertools.product(gamma, coef0):
            ocsvm_smd = OneClassSVM(kernel='sigmoid', gamma=r, coef0=c)
            ocsvm_smd.fit(item.values.reshape(-1, 1))
            prec_rec_f = precision_recall_fscore_support(y_valid_true, ocsvm_smd.predict(X_valid[column_name].values.reshape(-1, 1)))
            f_score.append(np.average(prec_rec_f[2]))
            
            idx.append([r, c])

        # F値が最大になる係数γを取得し、再度シグモイドカーネルで学習
        best_idx = idx[np.argmax(f_score)]
        # print('Sigmoid best_kernel: gamma:{:.4}, coef0:{:2.2}, f_score:{:.4}'.format(best_idx[0], best_idx[1], np.max(f_score)))
        ocscvm_smd = OneClassSVM(kernel='sigmoid', gamma=best_idx[0], coef0=best_idx[1])
        ocscvm_smd.fit(item.values.reshape(-1, 1))

        # 最適なパラメータを使用したRBFカーネル、多項式カーネル、シグモイドカーネルのテストデータに対する結果表示
        pred_rbf = ocsvm_rbf.predict(X_test[column_name].values.reshape(-1, 1))
        pred_poly = ocscvm_poly.predict(X_test[column_name].values.reshape(-1, 1))
        pred_smd = ocscvm_smd.predict(X_test[column_name].values.reshape(-1, 1))
        # print('--------------------')
        # print('One-class SVM result')
        # print('--------------------')
        # print('rbf kernel')
        # print_precision_recall_fscore(y_test_true, pred_rbf)
        # print('--------------------')
        # print('polynomial kernel')
        # print_precision_recall_fscore(y_test_true, pred_poly)
        # print('--------------------')
        # print('sigmoid kernel')
        # print_precision_recall_fscore(y_test_true, pred_smd)
        # print('--------------------')
        # fpr, tpr, thresholds = roc_curve(y_test_true, pred_rbf, pos_label=1)
        # roc_auc = roc_auc_score(y_test_true, pred_rbf)
        # plt.figure()
        # plt.plot(fpr, tpr, 'k--',label='rbf kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="r")

        # fpr, tpr, thresholds = roc_curve(y_test_true, pred_poly, pos_label=1)
        # roc_auc = roc_auc_score(y_test_true, pred_poly)
        # plt.plot(fpr, tpr, 'k--',label='poly kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="g")

        # fpr, tpr, thresholds = roc_curve(y_test_true, pred_smd, pos_label=1)
        # roc_auc = roc_auc_score(y_test_true, pred_smd)
        # plt.plot(fpr, tpr, 'k--',label='sigmoid kernel (AUC = %0.2f)' % roc_auc, lw=2, linestyle="-", color="b")

        # plt.xlim([-0.05, 1.05])
        # plt.ylim([-0.05, 1.05])
        # plt.xlabel('False Positive Rate')
        # plt.ylabel('True Positive Rate')
        # plt.title('{}Spectrim_ROC_curve'.format(column_name))
        # plt.legend(loc='lower right')
        # plt.savefig('./image/{}Spectrim__OCSVM_ROC_curve.png'.format(column_name))
        # plt.show()
        # print('--------------------')

        """
        # 最もAUCが高いモデルを使用して異常値を抽出する
        """
        roc_auc_rbf = roc_auc_score(y_test_true, pred_rbf)
        roc_auc_poly = roc_auc_score(y_test_true, pred_poly)
        roc_auc_smd = roc_auc_score(y_test_true, pred_smd)
        # 最もAUCが高い結果を取得
        roc_auc_max = max(roc_auc_rbf, roc_auc_poly, roc_auc_smd)
        # 最もAUCの高いモデルの予測結果と異常スコアの最大値の取得
        if roc_auc_rbf == roc_auc_max:
            pred = roc_auc_rbf
            maximum_anomaly_score = ocsvm_rbf.score_samples(X_test[column_name].values.reshape(-1, 1)).min()
            #print('{}nm Spectrim best AUC kernel'.format(column_name))
        elif roc_auc_poly == roc_auc_max:
            pred = roc_auc_rbf
            maximum_anomaly_score = ocscvm_poly.score_samples(X_test[column_name].values.reshape(-1, 1)).min()
            # print('{}nm Spectrim best AUC kernel'.format(column_name))
        else:
            pred = roc_auc_smd
            maximum_anomaly_score = ocscvm_smd.score_samples(X_test[column_name].values.reshape(-1, 1)).min()
            # print('{}nm Spectrim best AUC kernel'.format(column_name))
        # 1回以上、異常と判断した時カラム名を表示
        if np.count_nonzero(pred == -1) >= 1:
            print('異常と判断した波長:{}nm:'.format(column_name))
            print('異常スコアの最高値:{}:'.format(maximum_anomaly_score))
            # 異常と判断された波長一覧とその異常値スコア最大値の保存先
            anomaly_spectrim.append(column_name)
            anomaly_score.append(maximum_anomaly_score)
            anomaly_roc_auc_score.append(roc_auc_score(y_test_true, pred))

コメントアウトして可視化した結果のイメージ

f:id:sekihan_0290:20200917221539p:plain

f:id:sekihan_0290:20200917221553p:plain

# 異常波長とその異常スコアの最大値をDataFrame化
result_df = pd.DataFrame({
    'anomaly_spectrim':anomaly_spectrim,
    'anomaly_score':anomaly_score,
    'anomaly_roc_auc_score':anomaly_roc_auc_score
})
# 異常値スコアでソート
result_df = result_df.sort_values('anomaly_score')
# 異常度スコアが高いほど異常度が高いと読める用に値を反転
result_df['anomaly_score'] = result_df['anomaly_score'] * (-1)
# 異常度スコアの先頭5データを表示
print(result_df.head())

# ファイル保存ディレクトリの作成
# if not(os.path.isdir('result_anomaly_data')):
#     os.mkdir('result_anomaly_data')
# result_df.to_csv('./result_anomaly_data/result_anomaly.csv', index=False)

# 処理に要した時間を出力
print("Computation time:{0:.3f} sec".format(time.time() - start))
Empty DataFrame
Columns: [anomaly_spectrim, anomaly_score, anomaly_roc_auc_score]
Index: []
Computation time:48.589 sec

カラムごとのIsolationForest(iForest)の実行

二次元データで外れ値検知を実行している例はよく見ますが、複数列を持つデータフレームに対し、各列で学習と予測を行い異常カラムを検出するみたいな例は見なかったので実践に使う前のプロト版として作ってみたものを公開してみます。(あまりいい例じゃないですが)

想定として、ある時間に200~800nmまでの波長データで何かしらの値を取得できるセンサーがあり、そのログデータから異常値を予測したいものとします。

具体的には、ノイズを加えたsin波を例に異常検知を実装します。
カラムは、200~800まで0.5刻みで記録されています。
100%分類できてるみたいで面白くないけどあくまで例ということで
データは訓練用データ、検証用(異常値)データ、テスト用(異常値)データ、検証用(異常値)データです。
データの内容 - ./data/train:学習に使用するノイズの少ないsin波データ - ./data/test/normal:検証とテスト兼用の正常sin波データ - ./data/test/anomaly/test:異常値のテストデータ - ./data/test/anomaly/valid:異常値の検証用データ

まずデータ全体を見て異常検知を行ったあと、
カラムごとに異常検知を行い、カラムごとに複数回異常と判定されたカラムとその異常スコアの最大値をDataFrameとして出力してみます。

import os
import glob
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('seaborn')
# from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import IsolationForest
from sklearn.metrics import precision_recall_fscore_support,confusion_matrix
from sklearn. metrics import roc_curve, roc_auc_score
from sklearn.metrics import precision_recall_fscore_support

import warnings
warnings.filterwarnings('ignore')

分析用疑似データ作成

# 保存用フォルダ作成
if not(os.path.isdir('./data')):
    os.mkdir('./data')
if not(os.path.isdir('./data/train')):
    os.mkdir('./data/train')

if not(os.path.isdir('./data/test/normal')):
    os.mkdir('./data/test/normal')
if not(os.path.isdir('./data/test')):
    os.mkdir('./data/test')

if not(os.path.isdir('data/test/anomaly')):
    os.mkdir('data/test/anomaly')
if not(os.path.isdir('data/test/anomaly/test')):
    os.mkdir('data/test/anomaly/test')
if not(os.path.isdir('data/test/anomaly/valid')):
    os.mkdir('data/test/anomaly/valid')

データ内容の確認

それぞれどのようなデータが入るかを可視化して確認してみます。

def Make_data(index, eta=0.2, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 40個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 40 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    plt.subplots(figsize=(16, 9))                                            # 表示サイズ指定
    plt.scatter(X, y)                                                        # 散布図
    plt.show()
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    
    return df
# train
Make_data(1, 0.05).head()

f:id:sekihan_0290:20200917220919p:plain

1
0 -0.012546
1 0.050312
2 0.033680
3 0.025586
4 -0.013438
# valid_anomaly
Make_data(1, 0.15, outlier=2.5, mode='test').head()

f:id:sekihan_0290:20200917221004p:plain

1
0 -0.037638
1 0.140455
2 0.080079
3 0.045318
4 -0.082235
# test_anomaly
Make_data(1, 0.2, outlier=3, mode='test').head()

f:id:sekihan_0290:20200917221016p:plain

1
0 -0.050184
1 0.185526
2 0.103278
3 0.055184
4 -0.116633
# 可視化の内容を消去して再度宣言
def Make_data(index, eta=0.05, outlier=3, mode='train'):
    '''
    sin波にノイズを加えたデータを作成する
    '''
    np.random.seed(42)
    data_size = 1200                                                         # データ数
    X = np.linspace(0,1, data_size)                                          # 0~1まで20個データを作成する
    noise = np.random.uniform(low= -1.0, high=1.0, size=data_size) * eta    # -1~1までの値をデータサイズ個作成し、引数倍する
    y = np.sin(2.0 * np.pi * X) + noise                                      # sin波にノイズを追加する
    # 外れ値の設定
    # 100個ごとにランダムで0~引数outlierまでの1個の値を加算
    if mode == 'test':
        for cnt in range(data_size):
            if cnt % 100 == 0:
                y[cnt] += np.random.randint(0, outlier, 1)
    # DataFrameを引数をカラム名として作成
    df = pd.DataFrame({
        index:y
    })
    return df

訓練用データを作成

関数Make_data()を使用して訓練用データを作成します。
ノイズはデフォルトの0.2でデータ数160個作成します。

df_train = pd.DataFrame([])
for i in range(1, 161):
    df_base =  Make_data(i).T
    df_train = pd.concat([df_train, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_train.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_train.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

訓練データの保存

# index不要のとき
df_train.to_csv('data/train/X_train.csv', index=False)

検証とテスト兼用の正常sin波データ

ノイズはデフォルトの0.2でデータ数80個作成します。

df_normal = pd.DataFrame([])
for i in range(1, 81):
    df_base =  Make_data(i).T
    df_normal = pd.concat([df_normal, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_normal.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_normal.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
2 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
3 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
4 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991
5 -0.012546 0.050312 0.03368 0.025586 -0.013438 -0.008202 -0.012755 0.073292 0.052022 0.067953 ... -0.090411 -0.033694 -0.052086 -0.019345 -0.071625 0.016194 0.031628 0.036407 0.019725 -0.036991

5 rows × 1200 columns

検証とテスト兼用の正常sin波データの保存

# index不要のとき
df_normal.to_csv('./data/test/normal/X_normal.csv', index=False)

検証用(異常値)データを作成

ノイズは0.15で、外れ値は2.5、データ数30個作成します。

df_valid = pd.DataFrame([])
for i in range(1, 31):
    df_base =  Make_data(i, 0.15, outlier=2.5, mode='test').T
    df_valid = pd.concat([df_valid, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_valid.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_valid.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
2 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
3 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
4 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974
5 -0.037638 0.140455 0.080079 0.045318 -0.082235 -0.077003 -0.101138 0.146527 0.072245 0.109567 ... -0.176941 -0.017259 -0.082909 0.004838 -0.162476 0.090501 0.126326 0.130183 0.069655 -0.110974

5 rows × 1200 columns

検証(異常値)データの保存

# index不要のとき
df_valid.to_csv('./data/test/anomaly/valid/X_valid.csv', index=False)

テスト用(異常値)データを作成

ノイズはデフォルトの0.2で外れ値は3、データ数48個作成します。

df_test = pd.DataFrame([])
for i in range(1, 49):
    df_base =  Make_data(i, 0.2, outlier=3, mode='test').T
    df_test = pd.concat([df_test, df_base])
# 200nm~800nmをカラムとして1桁で丸めて設定
df_test.columns =  np.round(np.linspace(200, 800, 1200), 1)
df_test.head()
200.0 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 ... 795.5 796.0 796.5 797.0 797.5 798.0 798.5 799.0 799.5 800.0
1 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
2 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
3 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
4 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966
5 -0.050184 0.185526 0.103278 0.055184 -0.116633 -0.111403 -0.14533 0.183145 0.082357 0.130375 ... -0.220205 -0.009042 -0.098321 0.016929 -0.207902 0.127655 0.173675 0.177071 0.09462 -0.147966

5 rows × 1200 columns

テスト用(異常値)データの保存

# index不要のとき
df_test.to_csv('./data/test/anomaly/test/X_test.csv', index=False)

データの読み込みとラベル付与

start = time.time()             # 実行開始時間の取得
# 学習用フォルダ内のファイル一覧を取得
files_train = glob.glob('./data/train/*')
# 検証とテスト兼用の正常sin波データァイル一覧
files_normal  = glob.glob('./data/test/normal/*')
# 異常値のテストデータ
files_anomaly1 = glob.glob('./data/test/anomaly/test/*')
# 異常値の検証用データ
files_anomaly2 = glob.glob('./data/test/anomaly/valid/*')

# CSV形式のファイルを読み込み、学習データを全てx_trainに格納
X_train = pd.DataFrame([])
for file_name in files_train:
    csv = pd.read_csv(filepath_or_buffer=file_name)
    X_train = pd.concat([X_train, csv])
# StandardScalerでz標準化
# sc = StandardScaler()
# sc.fit(X_train)
# X_train = sc.transform(X_train)

正解ラベルの付与

# 正常データのラベルを1,異常データのラベルを-1として
# テストデータ用y_test_true,検証用データ用y_valid_trueに格納
x_test_normal, x_test_anomaly1, x_test_anomaly2, X_test, X_valid = pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([]), pd.DataFrame([])
y_test_true, y_valid_true = [], []

# 検証とテスト兼用の正常sin波データを読み込み
# 正常データのラベル1を付与
for file_name in files_normal:
    csv = pd.read_csv(file_name)
    x_test_normal = pd.concat([x_test_normal, csv])
    for i in range(0, len(csv)):
        y_test_true.append(1)
        y_valid_true.append(1)
# 異常値の検証用データを読み込み
# 異常データのラベル-1を付与
for file_name in files_anomaly1:
    csv = pd.read_csv(file_name)
    x_test_anomaly1 = pd.concat([x_test_anomaly1, csv])
    for i in range(0, len(csv)):
        y_test_true.append(-1)
# 異常値のテスト用データを読み込み
# 異常データのラベル-1を付与
for file_name in files_anomaly2:
    csv = pd.read_csv(file_name)
    x_test_anomaly2 = pd.concat([x_test_anomaly2, csv])
    for i in range(0, len(csv)):
        y_valid_true.append(-1)
        
# テストデータx_test,検証用データx_validを正常データと異常データを組み合わせて用意
X_test = pd.concat([x_test_normal, x_test_anomaly1])
X_valid = pd.concat([x_test_normal, x_test_anomaly2])
# z標準化
# X_test = sc.transform(X_test)
# X_valid = sc.transform(X_valid)
# 正常データ数,異常データ数(テストデータ),テストデータ総数,検証用データ総数を確認
print('data size')
print('訓練データ:{}'.format(X_train.shape[0]))
print('テスト・検証用正常データ:{}'.format(x_test_normal.shape[0]))
print('異常データ数(テストデータ):{}'.format(x_test_anomaly1.shape[0]))
print('異常データ数(検証データ):{}'.format(x_test_anomaly2.shape[0]))
print('テストデータ総数:{}'.format(X_test.shape[0]))
print('検証用データ総数:{}'.format(X_valid.shape[0]))
data size
訓練データ:160
テスト・検証用正常データ:80
異常データ数(テストデータ):48
異常データ数(検証データ):30
テストデータ総数:128
検証用データ総数:110

IsolationForest(iForest)の実行

データ全体を使用して学習を実行してみます。

# 探索するハイパーパラメータのリスト
estimators_params = [50, 100, 150, 200]

# アンサンブルする識別器の数n_estimatorsを変化させて検証用データに対するF値を取得
idx, f_score = [], []
for k in estimators_params:
    IF = IsolationForest(n_estimators=k, random_state=2)
    IF.fit(X_train)
    # 検証データの平均F値を追加
    # _predictで局所異常因子を異常度としてある閾値で切った時の2値ラベルの取得
    prec_rec_f = precision_recall_fscore_support(y_valid_true, IF.predict(X_valid))
    f_score.append(np.average(prec_rec_f[2]))
    idx.append(k)
# F値が最大となるアンサンブル数を取得し,iForestを再適合
plt.plot(idx, f_score, 'b-')
plt.xlabel('n_estimators')
plt.ylabel('F-score')
plt.title('ALL_Spectrim_iForest_F-score')
plt.show()

best_k = idx[np.argmax(f_score)]
print('Local Outlier Factor result (n_estimators={}):'.format(best_k))
IF = IsolationForest(n_estimators=best_k, random_state=42)
IF.fit(X_train)

f:id:sekihan_0290:20200917221047p:plain

Local Outlier Factor result (n_estimators=50):





IsolationForest(n_estimators=50, random_state=42)
# 平均精度・再現率・F値と混同行列の表示
y_pred = IF.predict(X_test)
prec_rec_f = precision_recall_fscore_support(y_test_true, y_pred)
print('Ave. Precision {:0.4f}, Ave. Recall {:0.4f}, Ave. F-score {:0.4f}'.format(np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
print('Confusion Matrix')
df = pd.DataFrame(confusion_matrix(y_test_true, y_pred))
df.columns = [u'anomaly', u'normal']
df
Ave. Precision 0.3125, Ave. Recall 0.5000, Ave. F-score 0.3846
Confusion Matrix
anomaly normal
0 0 48
1 0 80
fpr, tpr, thresholds = roc_curve(y_test_true, IF.decision_function(X_test), pos_label=1)
roc_auc = roc_auc_score(y_test_true, IF.decision_function(X_test))
plt.plot(fpr, tpr, 'b--',label='ROC for test data (AUC = {:0.2f})'.format(roc_auc), lw=2, linestyle='-')
plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC curve')
plt.legend(loc="lower right")
<matplotlib.legend.Legend at 0x23f347ceb50>

[f:id:sekihan_0290:20200917221102p:plain]

列ごとに学習を行い、異常発生時のカラムを検出する

カラムごとに学習を行い、テストデータで予測を行って42データ以上異常が予測されたものを異常値として抽出します。

ROC曲線、混同行列の表示は冗長なのでコメントアウトしてます。

# 異常と判断された波長一覧とその異常値スコア最大値、AUCの保存先
anomaly_spectrim = []
anomaly_score = []
anomaly_roc_auc_score = []
# ファイル保存ディレクトリの作成
# if not(os.path.isdir('image')):
#     os.mkdir('image')
for column_name, item in X_train.iteritems():
    # 探索ハイパーパラメータのリスト
    estimators_param = [50, 100, 150, 200]

    # アンサンブルする識別器の数を変化させ検証用データに対するF値を取得
    idx, f_score = [], []
    for k in estimators_param:
        IF = IsolationForest(n_estimators=k, random_state=42)
        # pd.Series.valuesでnumpy.ndarray型にし、np.reshape(-1, 1)でn行1列に変形
        IF.fit(item.values.reshape(-1, 1))
        # 正常と異常のクラスに対する平均F値を取得
        prec_rec_f = precision_recall_fscore_support(y_valid_true, IF.predict(X_valid[column_name].values.reshape(-1, 1)))
        f_score.append(np.average(prec_rec_f[2]))
        idx.append(k)
        
    # F値が最大になるアンサンブル数を取得し、iForestを再結合
    # plt.figure()
    # plt.plot(idx, f_score)
    # plt.xlabel('n_estimators')
    # plt.ylabel('F-score')
    # plt.title('{}Spectrim_iForest_F-score'.format(column_name))
    # plt.savefig('./image/{}Spectrim_iForest_F-score.png'.format(column_name))
    # plt.show()

    best_k = idx[np.argmax(f_score)]
    IF = IsolationForest(n_estimators=best_k, random_state=42)
    IF.fit(item.values.reshape(-1, 1))
    pred = IF.predict(X_test[column_name].values.reshape(-1, 1))
    # 最適なアンサンブル数を使用して、テストデータに対する結果を表示
    # print('--------------------')
    # print('Local Outlier Factor result (n_estimators={})'.format(best_k))
    # print('--------------------')
    # 平均精度・再現率・F値と混同行列の表示
    prec_rec_f = precision_recall_fscore_support(y_test_true, y_pred)
    # print('Ave. Precision %.4f, Ave. Recall %.4f, Ave. F-score %.4f'% (np.average(prec_rec_f[0]), np.average(prec_rec_f[1]), np.average(prec_rec_f[2])))
    # print('Confusion Matrix')
    df = pd.DataFrame(confusion_matrix(y_test_true, y_pred))
    df.columns = [u'anomaly', u'normal']
    # print(df)
    
    # print('--------------------')
    # fpr, tpr, thresholds = roc_curve(y_test_true, IF.decision_function(X_test[column_name].values.reshape(-1, 1)), pos_label=1)
    # roc_auc = roc_auc_score(y_test_true, IF.decision_function(X_test[column_name].values.reshape(-1, 1)))
    # plt.figure()
    # plt.plot(fpr, tpr, 'k--',label='ROC for test data (AUC = {:0.2f})'.format(roc_auc, lw=2, linestyle='-'))
    # plt.xlim([-0.05, 1.05])
    # plt.ylim([-0.05, 1.05])
    # plt.xlabel('False Positive Rate')
    # plt.ylabel('True Positive Rate')
    # plt.title('{}Spectrim_ROC_curve'.format(column_name))
    # plt.legend(loc='lower right')
    # plt.savefig('./image/{}Spectrim_iForest_ROC_curve.png'.format(column_name))
    # plt.show()
    # print('--------------------')
    # 42データ以上、異常と判断した時カラム名を表示
    if np.count_nonzero(pred == -1) >= 42:
        print('異常と判断した波長:{}nm:'.format(column_name))
        print('異常スコアの最高値:{}:'.format(IF.score_samples(X_test[column_name].values.reshape(-1, 1)).min()))
        # 異常と判断された波長一覧とその異常値スコア最大値の保存先
        anomaly_spectrim.append(column_name)
        anomaly_score.append(IF.score_samples(X_test[column_name].values.reshape(-1, 1)).min())
        anomaly_roc_auc_score.append(roc_auc_score(y_test_true, IF.predict(X_test[column_name].values.reshape(-1, 1)) ))

コメントアウトして可視化した結果のイメージ

f:id:sekihan_0290:20200917221115p:plain

f:id:sekihan_0290:20200917221127p:plain

# 異常波長とその異常スコアの最大値をDataFrame化
result_df = pd.DataFrame({
    'anomaly_spectrim':anomaly_spectrim,
    'anomaly_score':anomaly_score,
    'anomaly_roc_auc_score':anomaly_roc_auc_score
})
# 異常値スコアでソート
result_df = result_df.sort_values('anomaly_score')
# 異常度スコアが高いほど異常度が高いと読める用に値を反転
result_df['anomaly_score'] = result_df['anomaly_score'] * (-1)
# 異常度スコアの先頭5データを表示
print(result_df.head())

# ファイル保存ディレクトリの作成
# if not(os.path.isdir('result_anomaly_data')):
#     os.mkdir('result_anomaly_data')
# result_df.to_csv('./result_anomaly_data/result_anomaly.csv', index=False)

# 処理に要した時間を出力
print("Computation time:{0:.3f} sec".format(time.time() - start))
Empty DataFrame
Columns: [anomaly_spectrim, anomaly_score, anomaly_roc_auc_score]
Index: []
Computation time:838.087 sec

ディレクトリ構造の表示 tree

treeコマンドにパスを指定することで、そのパスのディレクトリ構造をツリー形式で表示できます。

>tree data
フォルダー パスの一覧:  ボリューム Windows
ボリューム シリアル番号は 22A2-82FD です
C:\DATA
├─test
│  ├─anomaly
│  │  ├─test
│  │  │  └─.ipynb_checkpoints
│  │  └─valid
│  └─normal
└─train

ラズパイ3をサーバー化してファイル共有する

RasPi 3B+ のファイルServer化
NAS[Network Attached Storage]構築手順書
2019年6月版

2019年6月にラズパイ3をサーバー化した時の手順メモです。 古い情報ですが、2020の9月でもできたので保存場所の統一も兼ねて投稿しました。


Windowsにインストールしたソフト

【サーバー構築に必要な Windowsソフト 】

最低限必要なのはSDカードにOSイメージを書き込むソフトとターミナルソフト 他は適宜

【ターミナルソフト】

TeraTerm

ネットワーク上に設置したRaspberry Piと接続して操作する

【SDカード操作関連】

Win32DiskImager

Raspberri PiのOSイメージをSDカードに書き込む。 直接OSを入れたほうがSDカードのパーティションも複雑に分割されないので

SD Formatter

SDカードを再フォーマットするのに使う。

【ファイル転送】

WinSCP

Raspberry Piとファイルのやり取りをする サーバーにFTPを入れないので必要


ラズパイの設定

1. OSのダウンロード

以下のリンクからRaspbian Stretch with desktop and recommended softwareをダウンロードする(ファイル名:2019-04-08-raspbian-stretch-full.zip) RaspbianDWPage Raspbian_Downloads_Page

2. SDへの書き込み・起動

SD Card Formatter Win32DiskImager(balenaEtcher等)等を使い書き込む ※Win32DiskImagerはWriteをクリックすること

3.ソフトウェアの更新

terminal上で以下のコマンドを入力する

sudo apt-get update
sudo apt-get upgrade

4.ファームウェアの更新

terminal上で以下のコマンドを入力する

sudo rpi-update

5.ラズパイの設定

ラズパイの設定からSSHを有効にする 念のためVNCも有効にする

再起動

sudo reboot

※この段階でVNCSSHが有効になっているので VNC Viewer または Tera Termが使用できる。

6.rootパスワードの設定

初期状態ではrootユーザのパスワードが設定されていない rootは管理ユーザーなので一応鍵をかけておく、 コマンドプロンプトで以下のコマンドを入力する

sudo passwd root

入力後、続けてパスワードを入力する 入力した文字列は表示されないが入力できている。 最後に以下の表示が出れば正常にパスが変更される。

passwd: password updated successfully

7.ネットワーク設定(IPアドレスの固定)

※意外とミスしやすいので注意

今後のSSH接続に備え、LAN上でのRaspberry PiIPアドレスを固定する 設定ファイルは/etc/network/interfaces まず、初期ファイルをコピーしてバックアップします。

sudo cp /etc/network/interfaces /etc/network/interfaces.org

以下コマンドで編集に入ります。(leafpadとしてテキストファイルで編集してますがvimでもnanoでもOK)

sudo leafpad /etc/dhcpcd.conf

末尾に以下の意味の内容を追加してください(環境によって適宜変える箇所があるので注意) - interface eth0 - interface eth0だと優先 - static ip_address=[設定したい固定IPアドレス]/24 - static routers=[デフォルトゲートウェイIPアドレス] - static domain_name_servers=[DNSサーバーのIPアドレス]

固定したいIPアドレスは 例として192.168.100.105とする

interface wlan0
static ip_address=192.168.100.105/24
static routers=192.168.100.1
static domain_name_servers=192.168.100.1

アドレスの設定は※以下で調べて自分の環境に合わせてください。

ネットワーク設定ができたので一度シャットダウンのコマンドは以下

sudo shutdown -h now

ラズパイの電源を再導入し、 コマンドプロンプト(端末など)を起動して

7.固定したIPアドレスの確認

ラズパイを再起動後、PC側で以下のコマンドでネットワークの状態を確認

ping 192.168.100.105

ping以降のコマンドは固定化したIPアドレス

8.SSHで接続

PCからTeraTermを立ち上げ、ラズパイのIPアドレスTCP/IPのホスト(T)に入力する

OKを入力すると以下の画面になれば接続に成功している

ユーザー名にpi、パスフレーズに初期設定で設定したパスワードを入力してOKを押す ログインできれば、TeraTerm上にコマンドプロンプトが表示


Sambaのインストール・設定

1. Sambaをインストール

以下のコマンドでインストール

sudo apt-get install samba

2. vimをインストール

Tera Term上で作業したいのでvimをインストール

sudo apt-get install vim

3. Sambaの設定

sudo vim /etc/samba/smb.conf

以下の意味をもつ内容をShare Definitionsセクションに記載する - [pi] ⇒ 共有名。Windows上で表示されるフォルダ名 - path = /home/pi ⇒ 共有するフォルダのパス - read only = No ⇒ 読み込みも書き込みも可能 - guest ok = yes ⇒ ゲスト権限で接続 - force user = pi ⇒ ファイル操作は、ユーザー:piにより実行

#======================= Share Definitions =======================
・
・
[pi]
path = /home/pi
read only = No
guest ok = Yes
force user = pi

4. 動作確認

Windows

エクスプローラーを開いて、左下のネットワークをクリック。\の後に以下の固定IPアドレスを入力

¥¥192.168.100.105

Mac

ファインダーを開いてCommand+Kを入力 サーバーアドレスを聞かれるので smb://[ラズパイのIPアドレス]を入力 ゲストを選択して接続を行う マウントするボリュームはpublicを選択


ファイルの転送

1. WinSCPの設定

WinSCPを起動しホスト名にラズパイのIPアドレス、ユーザー名とパスワードを入力


参考にしたリンク一覧

ちょっとLinuxのサーバーログを調べるときに使うコマンド

Linuxのログファイルで主になるのが/var/log/messagesになります。

リアルタイムでログを監視

tail -fを使用することで確認できます。

tail -f /var/log/messages

終了はCtrl+cで行います。

ログファイルから特定のメッセージのみを調査したい

grepを使用することで確認できます。

下記の例はoom-killerによりプロセスが強制終了したときのログのみを表示させています。

[root@server ~]# grep oom-killer /var/log/messages
Sep  7 09:05:47 server kernel: SysMonMgr invoked oom-killer: gfp_mask=0x280da, order=0, oom_score_adj=0
Sep  7 09:35:57 server kernel: SysMonMgr invoked oom-killer: gfp_mask=0x280da, order=0, oom_score_adj=0
Sep  7 09:45:42 server kernel: SysMonMgr invoked oom-killer: gfp_mask=0x280da, order=0, oom_score_adj=0
Sep  7 09:58:31 server kernel: SysMonMgr invoked oom-killer: gfp_mask=0x280da, order=0, oom_score_adj=0
Sep  7 09:58:33 server kernel: postgres invoked oom-killer: gfp_mask=0x200da, order=0, oom_score_adj=0

oom-killerにより強制終了した09:05:47のログを抽出して見ます。

[root@server ~]# grep 09:05:47 /var/log/messages
Sep  7 09:05:47 server kernel: SysMonMgr invoked oom-killer: gfp_mask=0x280da, order=0, oom_score_adj=0
Sep  7 09:05:47 server kernel: SysMonMgr cpuset=/ mems_allowed=0
Sep  7 09:05:47 server kernel: CPU: 2 PID: 230184 Comm: SysMonMgr Not tainted 3.10.0-327.el7.x86_64 #1
Sep  7 09:05:47 server kernel: Hardware name: Microsoft Corporation Virtual Machine/Virtual Machine, BIOS Hyper-V UEFI Release v4.0 12/17/2019
Sep  7 09:05:47 server kernel: ffff880b78b42280 00000000fdacee5d ffff88070ccdbb00 ffffffff816351f1
Sep  7 09:05:47 server kernel: ffff88070ccdbb90 ffffffff81630191 ffff880a372069f0 ffff880a37206a08
Sep  7 09:05:47 server kernel: ffffffff00000202 fffeefff00000000 0000000000000001 ffffffff81128803
Sep  7 09:05:47 server kernel: Call Trace:
Sep  7 09:05:47 server kernel: [<ffffffff816351f1>] dump_stack+0x19/0x1b
Sep  7 09:05:47 server kernel: [<ffffffff81630191>] dump_header+0x8e/0x214
Sep  7 09:05:47 server kernel: [<ffffffff81128803>] ? delayacct_end+0x63/0xb0

強制終了の原因は、SysMonMgr cpuset=/ mems_allowed=0より、メモリ不足によるものだとわかります。