「Pythonによる異常検知」を寄り道写経 ～ 第２章2.4節「高度な特徴抽出による異常検知」③主成分分析(PCA)
第２章「非時系列データにおける異常検知」
書籍の著者 曽我部東馬 先生、監修 曽我部完 先生
この記事は、テキスト「Pythonによる異常検知」第２章「非時系列データにおける異常検知」2.4節「高度な特徴抽出による異常検知」の通称「寄り道写経」を取り扱います。
５回連続で高度な特徴抽出技術を利用した異常検知を寄り道写経します！
今回は 主成分分析（PCA）による異常検知です。
ではテキストを開いて異常検知の旅に出発です🚀
はじめに
テキスト「Pythonによる異常検知」のご紹介
このシリーズは書籍「Pythonによる異常検知」（オーム社、「テキスト」と呼びます）の寄り道写経です。
テキストは、2021年2月に発売され、機械学習等の誤差関数から異常検知を解明して、異常検知に関する実践的なPythonコードを提供する素晴らしい書籍です。
とにかく「誤差関数」と「異常度」の強い結びつきを堪能できる１冊です。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「Pythonによる異常検知」第１版第６刷、著者 曽我部東馬、監修者 曽我部完、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
2.4 高度な特徴抽出による異常検知
主に Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
今回の寄り道ポイントです。
主成分分析で異常検知を行います。
scikit-learn の PCA で主成分分析を実行します。
学習データで異常検知モデルを作成して、未知データの異常検知を行います。
学習データと未知データは前回から少々変更しました。
次元削減が必要なので3次元にしました。
この記事で用いるライブラリをインポートします。
### インポート
# 数値計算、確率計算
import numpy as np
import pandas as pd
### 機械学習
from sklearn.decomposition import PCA
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
主成分分析（PCA）
主成分分析を異常検知に用いるアイデアは以下のとおりです。
入力データ（高次元）に対して主成分分析を実行して次元削減する
次元削減後のデータ（低次元）を高次元に戻す
（入力データに再構成される）入力データと高次元に戻したデータの誤差のユークリッド距離を異常度とする
次元削減時に入力データの重要な特徴が確保されていると再高次元化での差異は小さくなり、一方で、重要な特徴が失われていると再高次元化での差異は大きくなります。
各データ点の重要の特徴の違いが正常・異常の判別に使われる感じです。
なお、この差異を「再構成誤差」と呼びます。
誤差関数と異常検知のマリアージュ💞
テキストの異常度の数式表現をお借りします。
$$
\alpha(X) = ||F^{\top}FX-X||^2
$$
$${X}$$は入力データ、$${F}$$は高次元から低次元への変換や、低次元から高次元への変換に用いる変換行列です。
テキストによると最適な変換行列$${F}$$は「共分散行列の固有関数行列」です。たぶん、入力データの分散共分散行列の固有ベクトルだと思われます。
異常度は入力データと再高次元化データのユークリッド距離（L2ノルム）です。
なおこの記事の高次元・低次元の変換は、変換行列を用いず、scikit-learn の PCA で実行します。
準備
１．学習データの作成と可視化
前回までの2次元データに1次元追加した3次元データを作成します。
最初の2次元は前回までと同じ300個の乱数データを作成します。
200個はグループ0、100個はグループ1を想定して、別々の２変量正規分布で乱数を作成します。
3次元目は連続一様分布乱数で作成します。
### データの作成 ※３変数
# 乱数生成器の設定
rng = np.random.default_rng(seed=1)
# 乱数生成器の設定
rng = np.random.default_rng(seed=1) # 23 or 24
# グループ0の200個のデータ作成：２変量正規分布乱数
group0 = rng.multivariate_normal(mean=[1, 1], cov=[[0.3, 0.081], [0.081, 0.3]],
size=200)
# グループ1の100個のデータ作成：２変量正規分布乱数
group1 = rng.multivariate_normal(mean=[-1, -1], cov=[[0.1, 0.009], [0.009, 0.1]],
size=100)
# 3次元目を追加（一様分布乱数）
data_3rd = rng.uniform(low=-1, high=1, size=300).reshape(-1, 1)
# データの統合
data = np.vstack([group0, group1])
data = np.hstack([data, data_3rd])
print('data.shape:', data.shape)
print(data[:5, :], '…')
# グループラベルの作成
label_true = np.hstack([np.zeros(200), np.ones(100)]).astype(int)
print('\nグループラベル:', label_true[:10], '…')【実行結果】
学習データの散布図を描画して、分布などを確認しましょう。
### データの散布図の描画 ※変数1,2の組み合わせ（以後同様）
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 散布図の描画
sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=label_true, alpha=0.7,
palette=['tab:blue', 'tab:orange'], ax=ax)
# 凡例
handles, _ = ax.get_legend_handles_labels()
ax.legend(handles=handles, labels=['グループ0', 'グループ1'], title='正解ラベル')
# 修飾
ax.set(xlabel='変数1', ylabel='変数2', title='データの散布図')
ax.grid(lw=0.5);【実行結果】
1次元目と2次元目を用いた散布図は前回までと同じ様相を見せてくれます。
異常検知モデルの構築
テキストの STEP に沿って、学習データによる異常検知モデルの構築に取り組みます。
１．STEP 1：PCAの実行
PCAを実行して学習データの低次元化データ（２次元）を取得します。
主成分数は２です。
### STEP1 PCAの実行
# 主成分の数の設定
n_components = 2
# PCAの実行、低次元化データ（主成分得点）の取得
pca = PCA(n_components=n_components).fit(data)
data_trans = pca.transform(data)
# クラスタ中心点の取得
centers = pca.mean_【実行結果】なし
２．STEP 2：異常度の算出
PCA の inverse_transform で元の３次元に再高次元化します。
さらに numpy の norm（ノルム）を用いて、再高次元化データと学習データのユークリッド距離を算出します。
### STEP2 異常度の算出
# データの再高次元化
data_inv = pca.inverse_transform(data_trans)
# 異常度αの算出 : Pとデータの差（つまり誤差）のユークリッド距離
anomaly_score = np.linalg.norm(data_inv - data, ord=2, axis=1)
print('anomaly_score.shape:', anomaly_score.shape)
print('anomaly_score[:5] :',anomaly_score[:5], '…')【実行結果】
【ちなみに情報】再高次元化のギモン
テキストの最適な変換行列$${F}$$を用いた再高次元化と、PCA の inverse_transform を用いた再高次元化は、結果が異なりました。
理由は不明です。なぜかな・・・？
低次元化を PCA で実行しており、処理の一貫性を重視して、再高次元化も PCA で行います。
３．STEP 3：異常度の閾値の算出
分位点法による閾値を算出します。
分位点は 1% です。
### STEP3 異常度の閾値の算出 by 分位点法
# 分位%点の設定
q = 0.01
# 閾値の算出
num_quantile = int(round(len(anomaly_score) * q, 0))
anomaly_thres = sorted(anomaly_score)[::-1][num_quantile - 1]
# 閾値の表示
print('分位点数:', num_quantile)
print('閾値:', anomaly_thres)【実行結果】
閾値はおよそ 1.28 です。
４．STEP 3-2：学習データで異常検知を確認
学習データの異常度に閾値を適用して異常検知を実行し、結果を可視化します。
### STEP3-2 学習データの異常検知
# 正常/異常ラベルの算出
anomaly_label = np.array((anomaly_score >= anomaly_thres) + 0)
# 算出結果の表示
print('anomaly_label.shape:', anomaly_label.shape)
print('anomaly_label[:10] :',anomaly_label[:10], '…')
print('異常データ数:', anomaly_label.sum())【実行結果】
学習データのうち、３点が異常判定となるモデルになりました。
異常検知の結果を可視化します。
### STEP3-2 異常度と閾値を可視化して確認
## 異常度と閾値の可視化関数の定義
def plot_anomaly1(a_score, a_threshold, a_label, q):
# 設定
palette = ['tab:blue', 'tab:red']
sizes = [40, 80]
labels = ['正常', '異常']
# 正常データがない場合
if sum(a_label == 0) == 0:
palette = ['tab:red']
sizes = [80]
labels = ['異常']
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# データの散布図の描画
sns.scatterplot(x=range(len(a_score)), y=a_score,
hue=a_label, palette=palette,
size=a_label, sizes=sizes)
# 閾値の水平線の描画
ax.axhline(a_threshold, color='tab:red', ls='--')
handles, _ = ax.get_legend_handles_labels()
ax.legend(handles=handles, labels=labels, bbox_to_anchor=(1, 1))
ax.set(xlabel='データ番号', ylabel='異常度 $\\alpha$',
title=f'PCAによる異常度\n閾値 {a_threshold:.3f} (q={q:.0%})')
ax.grid(lw=0.5)
plt.show()
## 異常度の可視化
plot_anomaly1(anomaly_score, anomaly_thres, anomaly_label, q)【実行結果】
いい感じに異常値を判別していると思います。
変数１、変数２の散布図の形式で異常データを見てみましょう。
### 異常データの散布図の描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 学習データの散布図の描画
sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=anomaly_label,
palette=['tab:blue', 'tab:red'], ec='white', alpha=0.7,
size=anomaly_label, sizes=[30, 80], ax=ax)
# クラスタ中心点の描画
ax.scatter(centers[0], centers[1], marker='x', color='tab:red', s=200,
lw=3)
handles, _ = ax.get_legend_handles_labels()
ax.legend(handles=handles, labels=['正常', '異常'])
# 修飾
ax.set(xlabel='変数1', ylabel='変数2', title='未知データの散布図')
ax.grid(lw=0.5);【実行結果】
赤いデータ点を異常と捉えたようです。
横軸３，縦軸２付近の正常なデータ点と、３つのデータ点の特徴の違いはどこにあるのでしょう・・・？
未知データの準備
未知データを作成して、さきほど構築した異常検知モデルに適用して異常検知しましょう！
１．未知データの作成
### 未知データの作成
data_new = np.array([[-1, 2, -0.5], [-1, 0.3, 0.8], [0, -1.5, 0.2],
[0, -1, -0.4], [1, -1, 1], [2, 0, -1]])
print('未知データ:')
print(data_new)【実行結果】
６個のデータを作成しました。
１次元・２次元の値は前回までの未知データと同じです。
未知データを可視化します。
### 未知データの可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 学習データの散布図の描画
sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=anomaly_label,
palette=['tab:blue', 'tab:red'], ec='white', alpha=0.7,
size=anomaly_label, sizes=[30, 30], legend=False, ax=ax)
# 中心点の描画
ax.scatter(centers[0], centers[1], marker='x', color='tab:red', s=200,
lw=3)
# 未知データの散布図の描画
ax.scatter(data_new[:, 0], data_new[:, 1], color='lightpink', ec='tab:red',
s=80, label='未知データ')
# 修飾
ax.set(xlabel='変数1', ylabel='変数2', title='未知データの散布図')
ax.legend()
ax.grid(lw=0.5);【実行結果】
未知データの異常検知
いよいよ本番です！
テキストの STEP に沿って、未知データの異常検知に取り組みます。
１．STEP 4：未知データの異常度の算出
学習済みの主成分分析器 pca を用いて、未知データの低次元化・再高次元化を実行し、異常度を計算します。
### STEP4 未知データの異常度の算出
# 未知データをPCAで低次元化
data_trans_new = pca.transform(data_new)
# 低次元データを再高次元化
data_inv_new = pca.inverse_transform(data_trans_new)
# 異常度の算出
anomaly_score_new = np.linalg.norm(data_inv_new - data_new, ord=2, axis=1)
print('未知データの異常度 :', anomaly_score_new)【実行結果】
２．STEP 5 未知データの異常検知
異常度計算・閾値を適用して異常データを炙り出します！
### STEP5 未知データの異常検知
anomaly_label_new = np.array((anomaly_score_new >= anomaly_thres) + 0)
print('未知データの異常検知:', anomaly_label_new)【実行結果】
未知データの１番目、５番目、６番目が異常判定されました！
未知データの異常検知結果をデータフレーム化しましょう。
### 未知データの異常検知結果のデータフレームの作成
data_new_df = pd.DataFrame(data_new, columns=['変数1', '変数2', '変数3'])
data_new_df['異常度'] = anomaly_score_new
data_new_df['閾値'] = anomaly_thres
data_new_df['異常検知'] = anomaly_label_new
display(data_new_df)【実行結果】
３．異常検知結果を可視化
未知データの異常検知結果をさまざまな可視化手法で堪能しましょう。
① 異常度の散布図
### 異常度の可視化
plot_anomaly1(anomaly_score_new, anomaly_thres, anomaly_label_new, q)【実行結果】
未知データの異常度はさまざまな値をとっていることが分かります。
② ２変数の散布図に正常・異常をマップ
学習データを含む２変数の散布図に未知データを追加してみましょう。
未知データは正常・異常を識別してプロットします。
### 未知データの異常検知の可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 学習データの散布図の描画
sns.scatterplot(x=data[:, 0], y=data[:, 1], hue=anomaly_label,
palette=['tab:blue', 'tab:red'], ec='white', alpha=0.7,
size=anomaly_label, sizes=[30, 30], legend=False, ax=ax)
# 中心点の描画
ax.scatter(centers[0], centers[1], marker='x', color='tab:red', s=200,
lw=3)
# 未知データの散布図の描画
sns.scatterplot(x=data_new[:, 0], y=data_new[:, 1],
hue=anomaly_label_new, palette=['tab:blue', 'tab:red'],
alpha=0.7, ec='wheat',
style=anomaly_label_new, markers=['o', 'X'],
size=anomaly_label_new, sizes=[200, 250], legend=False,
ax=ax)
for i in range(len(data_new)):
ax.annotate(text=i, xy=data_new[i, :2]+0.12, fontsize=14)
# 凡例
ax.scatter([], [], marker='o', c='tab:blue', ec='wheat', s=70, label='正常')
ax.scatter([], [], marker='X', c='tab:red', ec='wheat', s=100, label='異常')
ax.legend(bbox_to_anchor=(1.2, 1))
# 修飾
ax.set(xlabel='変数1', ylabel='変数2', title='未知データの異常検知')
ax.grid(lw=0.5);【実行結果】
左下から右上方向のデータの向き（分散最大化の方向）を想像し、その向きの外側に位置する未知データが異常判定されるように感じます。
EM法などのクラスタを識別する方法と比べて、主成分分析による異常検知はクラスタを識別しないため、学習データ全体の特徴から外れるデータ点が異常になるような感じです。
（参考：EM法の異常検知）
今回の寄り道写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？