Pythonでデータ分析：主成分分析（PCA）による異常検知

導入

データ分析の種類の一つとして、教師なし学習による異常検知というものがあります。ほとんどが正常なデータでまれに異常なデータが混じっている、その異常発生のパターンや異常と他の要因との紐付きがいまいちつかみきれていないというような場合、教師あり学習による2値分類がうまくワークしない、といった状況がありえます。そういった場合には、正常パターンを教師なし学習で学び、その正常パターンから外れているものを異常とする、という方法が有効です。

異常検知の方法ですが、正常パターンを正規分布でフィットするようなものから状態空間モデルに即したものまで多くのバリエーションがあります。その中の一つに、主成分分析（PCA）をベースにした異常検知の方法があります。正常データがある程度相関を持っているような多変数正規分布に従っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています（正規分布から外れている場合は、Box.Cox変換などで無理やり正規分布に近づける、他の手法を使うなどを検討します）。この方法ですが、Rでは実行例は井出さんの異常検知の本
www.coronasha.co.jp
にありますが、Pythonでの実行例というのが書籍やネットを探してもいまいち見つかりません。そのため、メモがてらこちらに簡単に残しておこうと思います。

多変数正規分布による人工データ

異常検知のデモ用に、人工データを用意します。PCAによる異常検知の場合、データは正規分布に従い、かついくつかの変数で相関が高いような場合にうまくはまります。というわけで、多変数正規分布で変数間にいくらか相関のあるようなデータを生成します。

# 必要なライブラリ
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# 多変数正規分布に従う乱数を生成
mean = (4, 2, 3)
sigma = [[1, 0.7, 0.75], [0.7, 1, 0.3], [0.75, 0.3, 1]]
df = np.random.multivariate_normal(mean, sigma, 500)

主成分分析（PCA）

※以下のコードは、無駄にDataFrameが使われていますが、もっとすっきり書ける可能性があります。

主成分分析は、次元圧縮の方法の一つです。データセットに含まれる変数が非常に多いがいくらか相関が強いものがあるとき、適切に基底の回転・取り替えを行うことによってより少ない軸でデータを近似することが期待できます。主成分分析（PCA）は今のデータセットから新しい基底（主成分）を求め、その基底のうちデータの情報を失わない（軸に射影しても分散が小さくならない）ものでいくらか代表して近似してしまおう、といった手法になっています。

PCAは、sklearnのメソッドで簡単に行うことができます。

# 異常度算出の都合上、正規化を行う
std = StandardScaler()
X = DataFrame(std.fit_transform(df))

from sklearn.decomposition import PCA
# インスタンス化
pca = PCA()
# 最初の1/5を学習に用いる
pca.fit(X[0:round(X.shape[0] / 5)])

pca.fitでレコード数×特徴量数の行列に関して特異値分解（もしくはこの行列と転置をかけた共分散行列というものの固有値分解）を行い、新しい基底を求めています。今回、変換前の軸は3つあったのですが、変換後は2つの軸でデータを近似することにします（こちらに関しては、本来は軸の寄与率などでどれくらい軸を採用するかを決めるなどします）。新しい基底はcomponents_というメソッドで引き出すことができ、その行方向を制限します。

Um = DataFrame(pca.components_[0:2, :]).T

ここで、後の都合上転置をしています。また、共分散行列というもの（ $S^2 = {I}_2 - U_m U_m^T$ ）を以下で定義しておきます。

S2 = np.identity(X.shape[1]) - Um.dot(Um.T)

異常度算出

PCAで次元圧縮した際に用いたデータは、正常なデータのみ、もしくはほとんど異常なデータは存在しないとします。すると、正常データは近似した軸で貼られる部分空間上でうまく近似でき、逆に異常なデータはその部分空間から大きく離れるのではないかと期待されます。この部分空間は正常部分空間と呼ばれます。異常度を評価する一つの方法が、データと正常部分空間からの距離を用いる方法です。今、正規化（中心化でもよい）されたデータを $\vec{x}$ （次元は特徴量の数）とすると、異常度 $a(\vec{x})$ は次のように書かれます。
$$a(\vec{x}) = \vec{x}^T S^2 \vec{x}^T$$
導出は井出さんの異常検知の本などを参照してください。

Pythonでは、おそらくライブラリで異常度を出力してくれるメソッドはないので、手で書く必要があります。ここでは、1レコードに対して異常度を計算する関数を作っておき、それを検証用データに対してループさせて使う、としています。

1レコードに対して異常度算出
def CalcAnomaly(x, S2):
    return x.T.dot(S2).dot(x)

# 検証データ（全データ）に対して異常度算出
a = [CalcAnomaly(X.iloc[i, :], S2) for i in range(X.shape[0])]

検証・可視化

各レコードに対し、異常度を可視化してみます。

import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(figsize = (15, 8))
plt.plot(a)
plt.xlabel('n_step')
plt.ylabel('anomaly score')

f:id:tekenuko:20171016205305p:plain

概ね1より小さい値になっています。まれに1程度の相対的に異常度が高い値が出力されていますが、こういった値は乱数の関係上発生しうるものだと思われます。本当に異常なデータがある場合は、今のスケール感でいうと数十といった値が出てくることもあるので、今回はほぼ正常データしかないだろうと思われるデータになっています。

また、もとの空間での散布図と異常度の関係も可視化しておきましょう。3軸なので、散布図も3次元でプロットしています。

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (12, 10))
ax = Axes3D(fig)
p = ax.scatter(df[:, 0], df[:, 1], df[:, 2], marker = 'o', c = a, cmap = 'jet', s = 50)
fig.colorbar(p)

f:id:tekenuko:20171016205724p:plain

新しい軸はデータの相関のある方向（図で言うと斜め方向）などになっているのですが、その軸たちから離れているデータほど異常度が大きくなっている、といった様子が図から読み取れます。このように、PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。

Next Step

今回は、主成分分析（PCA）をベースにした異常検知の方法があります。PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。異常検知の方法は他にも色々あるので、機会があればまた別の方法も紹介したいと思います。