「異常検知と変化検知」は異常検知の本です。アルゴリズム部分をpythonで実装していきたいと思います。たくさんの方が同じ内容をより良い記事でとして投稿しています。
個人的な勉強のメモ書きとなります。
私には難しい内容であり、コードがかなり汚いです。申し訳ありません。
まとめることが目的ではないので本文について参考書とほぼ同じ表現となっています。問題があればお教えください。
興味を持ったり、導出や詳細などを知りたい方は「異常検知と変化検知」を読んでいただければと思います。
参考
前回
ガウス過程回帰による異常検知
部分空間法による変化検知
累積和法: 変化検知の古典技術
ある量$\xi$が、多少のばらつきがあったとしても安定して$\mu$という値を取ることが分かっていれば、正常時のモデルとして正規分布$N(\xi|\mu,\sigma^2)$を考えることは自然である。ここで、$\sigma^2$は$\xi$の分散である。
さらに過去の異常事例の分析から、「正常時から上に$\nu_+$振れるようなら異常である」というような事前知識があるとする。
正常状態と異常状態は次のようなモデルで表現できることになる。
- 正常状態：$N(\xi|\mu,\sigma^2)$
- 異常状態：$N(\xi|\mu+\nu_+,\sigma^2)$
現在の状態がどちらに近いかを調べることにより、正常状態から異常状態への変化を検知できることになる。
時刻$t$における変化度として、
$$
\begin{align}
a(\xi^{(t)})&=\ln\frac{\mathcal{N}(\xi^{(t)}|\mu+\nu_+,\sigma^2)}{\mathcal{N}(\xi^{(t)}|\mu,\sigma^2)}\
&=\bigg(\frac{\nu_+}{\sigma} \bigg)\frac{\xi^{(t)}-\mu-(\nu_+/2)}{\sigma}
\end{align}
$$
と定義できる。$\xi^{(t)}$は時刻$t$における観測値である。
継続的に異常状態が生じているときは、上記の量は継続的に正の値を取る。
したがって、時刻$t$における変化度を漸化式の形で
\begin{align}
a_+^{(t)}&=\bigg[a_+^{(t-1)}+a(\xi^{(t)}) \bigg]_+\\
&=\bigg[a_+^{(t-1)}+\bigg(\frac{\nu_+}{\sigma} \bigg)\frac{\xi^{(t)}-\mu-(\nu_+/2)}{\sigma}) \bigg]_+
\end{align}
のように定義できる。$[・]_+$はカッコの中が正なら何もせず、負なら0に変える演算を意味する。
このように定義される変化度を累積和(CUMSUM)あるいは特に上側累積和(upper CUMSUM)と呼ぶ。累積和統計量(CUMSUM statistic)と呼ばれることもある。
下振れも監視するためには、
\begin{align}
a_-^{(t)}&=\bigg[a_-^{(t-1)}-a(\xi^{(t)}) \bigg]_+\\
&=\bigg[a_-^{(t-1)}-\bigg(\frac{\nu_-}{\sigma} \bigg)\frac{\xi^{(t)}-\mu+(\nu_-/2)}{\sigma}) \bigg]_+
\end{align}
という量を定義する。これは下側累積和(lower CUMSUM)と呼ぶ。
累積和ないしそれに類する統計量に上限や下限を付して時系列的に監視するための図を管理図ないしシューハート管理図などと呼ぶ。
アルゴリズム(累積和法による変化検知)
- 上振れと下振れの目安$\nu_+,\nu_-$を与える。上側累積和、下側累積和のそれぞれについて閾値$a_+^{th},a_-^{th}$を与える。
- 各時刻において、上側・下側の累積和を
a_+^{(t)}=\bigg[a_+^{(t-1)}+\bigg(\frac{\nu_+}{\sigma} \bigg)\frac{\xi^{(t)}-\mu-(\nu_+/2)}{\sigma}) \bigg]_+\\
a_-^{(t)}=\bigg[a_-^{(t-1)}-\bigg(\frac{\nu_-}{\sigma} \bigg)\frac{\xi^{(t)}-\mu+(\nu_-/2)}{\sigma}) \bigg]_+
にて計算する。
3. 上側・下側累積和のいずれかが閾値を超えていたら警報を出す。
import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0,100)
y1 = np.concatenate([np.random.normal(loc=5,scale=0.1, size=60),np.random.normal(loc=6,scale=0.2, size=60)])
y2 = np.concatenate([np.random.normal(loc=5,scale=0.1, size=60),np.random.normal(loc=4,scale=0.2, size=60)])
plt.plot(y1)
plt.plot(y2)
# 上振れと下振れの目安
nu_plus = 1
nu_minus = 1
# パラメータ
mu = 5
sigma = 0.1
# 変化度
LogLR_upper1 = (nu_plus/sigma)*((y1-mu-(nu_plus/2))/sigma)
LogLR_lower1 = (nu_minus/sigma)*((y1-mu+(nu_minus/2))/sigma)
# 上側・下側累積和の計算
UpperCumsum1 = [0]
LowerCumsum1 = [0]
for i in range(1, len(LogLR_upper1)):
if UpperCumsum1[i-1] + LogLR_upper1[i] < 0:
UpperCumsum1.append(0)
else:
UpperCumsum1.append(UpperCumsum1[i-1] + LogLR_upper1[i])
if LowerCumsum1[i-1] - LogLR_lower1[i] < 0:
LowerCumsum1.append(0)
else:
LowerCumsum1.append(LowerCumsum2[i-1] - LogLR_lower1[i])
# 管理図
fig, ax = plt.subplots(3,1,figsize=(8,5))
ax[0].plot(y1)
ax[1].plot(LogLR_upper1)
ax[1].plot(LogLR_lower1)
ax[2].plot(UpperCumsum1)
ax[2].plot(LowerCumsum1)
# 変化度
LogLR_upper2 = (nu_plus/sigma)*((y2-mu-(nu_plus/2))/sigma)
LogLR_lower2 = (nu_minus/sigma)*((y2-mu+(nu_minus/2))/sigma)
# 上側・下側累積和の計算
UpperCumsum2 = [0]
LowerCumsum2 = [0]
for i in range(1, len(LogLR_upper2)):
if UpperCumsum2[i-1] + LogLR_upper2[i] < 0:
UpperCumsum2.append(0)
else:
UpperCumsum2.append(UpperCumsum2[i-1] + LogLR_upper2[i])
if LowerCumsum2[i-1] - LogLR_lower2[i] < 0:
LowerCumsum2.append(0)
else:
LowerCumsum2.append(LowerCumsum2[i-1] - LogLR_lower2[i])
# 管理図
fig, ax = plt.subplots(3,1,figsize=(8,5))
ax[0].plot(y2)
ax[1].plot(LogLR_upper2)
ax[1].plot(LogLR_lower2)
ax[2].plot(UpperCumsum2)
ax[2].plot(LowerCumsum2)
近傍法による異常部位検出
観測値として長さ$T$の時系列が
$$
D={\xi^{(1)},\xi^{(2)},\cdots,\xi^{(T)} }
$$
のように与えられていると考える。
$M$個の隣接した観測値をまとめて
$$
\boldsymbol{x}^{(1)}=\begin{pmatrix} \xi^{(1)}\\xi^{(2)}\\vdots\\xi^{(M)} \end{pmatrix},\ \boldsymbol{x}^{(2)}=\begin{pmatrix} \xi^{(2)}\\xi^{(3)}\\vdots\\xi^{(M+1)} \end{pmatrix},\cdots
$$
のように、データを$M$次元ベクトルの集まりとして表すことにする。
これにより、長さ$T$の観測値からなる時系列データは、
$$
N=T-M+1
$$
本の$M$次元ベクトルに変換される。長さ$M$の窓のことをスライド窓ないし滑走窓と呼ぶ。
スライド窓により生成したベクトルのことを、部分時系列と呼ぶことがある。
異常部位検出では、$k=1$とした近傍法による外れ値検知の方法が多くみられる。
ある標本$\boldsymbol{x}^{(n)}$から最近傍標本までの距離を$\epsilon^{(n)}$とすれば、
$$
a(\boldsymbol{x}^{(n)})=M\ln\epsilon^{(n)}
$$
により異常度を計算できる。
from scipy.spatial.distance import cdist
# 窓幅
M = 3
# 部分時系列
X = np.array([y1[i:i+M] for i in range(len(y1)-M+1)])
# 距離の計算
dist = cdist(X,X)
dist[dist==0]=100 # 近傍点に自分が入らないようにする
# 異常度
a = M*np.log(np.min(dist, axis=1))
# 可視化
fig, ax = plt.subplots(2,1,figsize=(8,4))
ax[0].plot(y1)
ax[1].plot(a)
変化検知問題と密度比
長さ$T$の1次元時系列が
$$
D={\xi^{(1)},\xi^{(2)},\cdots,\xi^{(T)} }
$$
のように与えられており、窓幅$M$のスライド窓を使い$M$次元ベクトルの
$$
D={\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},\cdots,\boldsymbol{x}^{(N)} }
$$
に変換したとする。ここで、$N=T-M+1$である。
時刻$t$において過去側と現在側に2つの領域を設定し、その領域内に入るベクトルを使って、それぞれの領域で確率分布を推定する。
過去側の分布を$p^{(t)}(\boldsymbol{x})$とし、現在側の分布を$p'^{(t)}(\boldsymbol{x})$とする。
変化検知問題とは、本質的には、この2つの確率分布の相違度を計算する問題として定義できる。
ネイマン・ピアソン則を拡張させると、時刻$t$における変化度は
$$
a^{(t)}=\int d\boldsymbol{x} p^{(t)}(\boldsymbol{x})\ln\frac{p^{(t)}(\boldsymbol{x})}{p'^{(t)}(\boldsymbol{x})}
$$
のように定義することができる。したがって、変化検知問題とは、本質的には密度比の逐次推定問題である。
特異スペクトル変換法
特異値分解による特徴的なパターンの自動抽出
次の問題を考える:$D$に含まれる部分時系列において、もっとも典型的な波形は何か。
ここで、
$$
\boldsymbol{x}^{(1)}v_1+\boldsymbol{x}^{(2)}v_2+\cdots+\boldsymbol{x}^{(N)}v_N=X\boldsymbol{v}
$$
のような一般の1次結合を考える。ただし、$X=[\boldsymbol{x}^{(1)},\cdots,\boldsymbol{x}^{(N)}]$および$\boldsymbol{v}=[v_1,\cdots,v_N]^T$と置いた。
直感的にはこの問題は
$$
maximize\ |X\boldsymbol{v}|^2\ subject\ to\ \boldsymbol{v}^T\boldsymbol{v}=1
$$
により求めればよいことがわかる。ラグランジュ法により、条件式
$$
X^TX\boldsymbol{v}=\gamma\boldsymbol{v}
$$
が得られる。$\boldsymbol{v}$として、行列$X^TX$の規格化された固有ベクトルを選べばよいことを意味する。
上式に左から$X$をかけて
$$
\boldsymbol{u}=X\boldsymbol{v}/\sqrt{\gamma}
$$
とおくと、$\boldsymbol{u}$が別の固有値方程式
$$
XX^T\boldsymbol{u}=\gamma\boldsymbol{u}
$$
を満たすことがただちにわかる。$\boldsymbol{u}$と$\boldsymbol{v}$の立場を入れ替えれば
$$
\boldsymbol{v}=X^T\boldsymbol{u}/\sqrt{\gamma}
$$
も得られる。
$M<N$とすると、固有値方程式の解$(\gamma_1,\boldsymbol{v}_1),(\gamma_2,\boldsymbol{v}_2),\cdots$に対応して、正規直行ベクトルの2つの集合${\boldsymbol{v}_i}$と${\boldsymbol{u}_i}$が得られ、これらは、
$$
X\boldsymbol{v}_i=\sqrt{\gamma_i}\boldsymbol{u}_i\
X^T\boldsymbol{u}_i=\sqrt{\gamma_i}\boldsymbol{v}_i\
$$
を満たす。
いま、$U=[\boldsymbol{u}_1,\cdots,\boldsymbol{u}_M]$と$V=[\boldsymbol{v}_1,\cdots,\boldsymbol{v}_M]$とおくと、
$$
X^TU=V\Gamma^{1/2}
$$
となる。これから、
$$
X=U\Gamma^{1/2}V^T
$$
という関係式が得られる。
これを$X$の特異値分解(SVD)と呼ぶ。
${\boldsymbol{u}_i}$を左特異ベクトル、${\boldsymbol{u}_i}$を右特異ベクトル、${\gamma_i}$を特異値と呼ぶ。
特異値分解による$X$の低階数近似はノイズ除去をやっているのと同じである。
変化度の定義
過去側で$n$本、現在側で$k$本の部分時系列を使って2つのデータ行列
$$
X^{(t)}=[\boldsymbol{x}^{(t-n-M+1)},\cdots,\boldsymbol{x}^{(t-M-1)},\boldsymbol{x}^{(t-M)}]\
Z^{(t)}=[\boldsymbol{x}^{(t-k+L-M+1)},\cdots,\boldsymbol{x}^{(t-M+L-1)},\boldsymbol{x}^{(t-M+L)}]
$$
を定義する。
$X^{(t)}$は$\xi^{(t-n-M+1)}$から$\xi^{(t-1)}$までのデータで構成され、履歴行列と呼ぶ。
対して、$Z$はテスト行列と呼ぶ。$L$は、ラグと呼ばれる。
2つの行列の特異値分解を行い、$X^{(t)}$において$r$本、$Z^{(t)}$において$m$本の左特異ベクトルを選ぶことにすれば、
$$
U_r^{(t)}=[\boldsymbol{u}^{(t,1)},\boldsymbol{u}^{(t,2)},\cdots,\boldsymbol{u}^{(t,r)}]\
Q_m^{(t)}=[\boldsymbol{q}^{(t,1)},\boldsymbol{q}^{(t,2)},\cdots,\boldsymbol{q}^{(t,m)}]
$$
となり、それぞれの列ベクトルで張られる空間を主部分空間と呼ぶ。
異常度として、
$$
\begin{align}
a(t)&=1-|U_r^{(t)T}Q_m^{(t)}|_2\
&=1-\big(U_r^{(t)T}Q_m^{(t)}の最大特異値 \big)
\end{align}
$$
という定義が考えられる。ここで、$||_2$は行列2ノルムを意味する。
特異値分解により時系列データの特徴パターンを求め、それに基づき変化度を求める手法を**特異スペクトル変換(特異スペクトル解析)**と呼ぶ。
アルゴリズム(特異スペクトル変換)
時系列$D={\xi^{(1)},\xi^{(2)},\cdots,\xi^{(T)} }$を用意する。窓幅$M$、履歴行列の列数$n$とパターン数$r$、テスト行列の列数$k$とパターン数$m$、ラグ$L$を決める。$t=(M+k),\cdots,(T-L+1)$において次の計算を行う。
- 履歴行列とテスト行列の構成
$$
X^{(t)}=[\boldsymbol{x}^{(t-n-M+1)},\cdots,\boldsymbol{x}^{(t-M-1)},\boldsymbol{x}^{(t-M)}]\
Z^{(t)}=[\boldsymbol{x}^{(t-k+L-M+1)},\cdots,\boldsymbol{x}^{(t-M+L-1)},\boldsymbol{x}^{(t-M+L)}]
$$
を作る。- 特異値分解
$X^{(t)},Z^{(t)}$を特異値分解し、左特異ベクトルの行列$U_r^{(t)},Q_m^{(t)}$を求める。- スコアの計算
$U_r^{(t)T}Q_m^{(t)}$の最大特異値を計算し、
$$
a(t)=1-\big(U_r^{(t)T}Q_m^{(t)}の最大特異値 \big)
$$
により異常度を計算する。
def singular_spectrum_analysys(X, t, n, k, r, m, L):
# 過去側のデータ行列(履歴行列)
Xt = X[:,(t-n-M):(t-M)]
# 現在のデータ行列(テスト行列)
Zt = X[:,(t-k+L-M+1):(t-M+L+1)]
# 特異値分解
U = np.linalg.svd(Xt)[0][:,:r]
Q = np.linalg.svd(Zt)[0][:,:m]
# UQの最大固有値
max_eig_value = np.linalg.svd(U.T@Q)[1][0]
a = 1 - max_eig_value
return a
M = 5# 窓幅
n = 4# 履歴行列の列数
r = 2# 履歴行列のパターン数
k = 3# テスト行列の列数
m = 3# テスト行列のパターン数
L = 4# ラグ
# 部分時系列の作成
X = np.array([y1[i:i+M] for i in range(len(y1)-M+1)]).T
T = len(X.T)
a = [0]*(M+n+L-1)
for t in range(M+n, T+M-L):
a.append(singular_spectrum_analysys(X, t, n, k, r, m, L))
fig, ax = plt.subplots(2,1,figsize=(8,4))
ax[0].plot(y1)
ax[1].plot(a)
ランチョス法による特異スペクトル変換の高速化
次の問題を考える。
固有射影計算問題
$M$次元実対称行列$C$と、$\boldsymbol{a}^T\boldsymbol{a}=1$を満たす$M$次元ベクトル$\boldsymbol{a}$が与えられている。
$C$の固有ベクトルを$\boldsymbol{u}^{(1)},\boldsymbol{u}^{(2)},\cdots$とするとき、これらと$\boldsymbol{a}$の内積
$$
h_i(C,\boldsymbol{a})=\boldsymbol{a}^T\boldsymbol{u}^{(i)}
$$
を計算せよ
特異スペクトル変換の変化度で$m=1$とすれば
$$
a(t)=1-\sqrt{\sum_{i=1}^rh_i(XX^T,\boldsymbol{q}^{(t,1)})^2}
$$
のように、$\boldsymbol{q}^{(t,1)}$の計算と内積の計算に帰着される。
解法は次のように与えられる。
アルゴリズム(ランチョス法)
- $C$と$\boldsymbol{a}$を与える。3重対角行列の次元$s$を決める。
$\boldsymbol{a}_0=\boldsymbol{0},\beta_0=1,\boldsymbol{r}=\boldsymbol{a}$とおく。- $j=1,2,\cdots,s$について次式を反復する。
\boldsymbol{a}_j←\boldsymbol{r}/\beta_{j-1}\\
\alpha_j←\boldsymbol{a}_j^TC\boldsymbol{a}_j\\
\boldsymbol{r}←C\boldsymbol{a}_j-\alpha_j\boldsymbol{a}_j-\beta_{j-1}\boldsymbol{a}_{j-1}\\
\beta_j←|\boldsymbol{r}|
- $s×s$の3重対角行列$T_s$を結果として返す。
アルゴリズム(固有射影計算問題の解法)
- $C$と$\boldsymbol{a}$を与える。必要な固有射影の個数$s$を与える。
- ランチョス法の計算を行い、3重対角行列$T_s$を求める。
- $T_s$の固有値分解を行い、$s$個の固有ベクトル$\boldsymbol{b}_1,\cdots,\boldsymbol{b}_s$を求める。
- $h_i(C,\boldsymbol{a})=(\boldsymbol{b}_iの第1成分)$
def eig_value_lanczos(X, a, M, s, r):
# Cの計算
C = X@X.T
#a=np.ones((M,1))/np.sqrt(M)
# 結果の初期化
aj = np.zeros((1+s,M,1))
betaj = np.ones((1+s,1))
_r = a
# パラメータの更新
for i in range(1,s+1):
aj[i] = _r/betaj[i-1]
alpha = aj[i].T@C@aj[i]
_r = C@aj[i]-alpha*aj[i]-betaj[i-1]*aj[i-1]
betaj[i] = np.sqrt(np.sum(_r**2))
# 3重対角行列の計算
Ts = aj[1:,:,0]@C@aj[1:,:,0].T
return np.sqrt(np.sum(np.linalg.eig(Ts)[1][0,:r]**2))
def singular_spectrum_analysys(X, t, n, k, r, m, L):
# 過去側のデータ行列(履歴行列)
Xt = X[:,(t-n-M):(t-M)]
# 現在のデータ行列(テスト行列)
Zt = X[:,(t-k+L-M+1):(t-M+L+1)]
# 特異値分解
#U = np.linalg.svd(Xt)[0][:,:r]
Q = np.linalg.svd(Zt)[0][:,:m]
# 異常度
a=1-eig_value_lanczos(X, Q[:,0:1], M, 3, r)
return a
M = 5# 窓幅
n = 4# 履歴行列の列数
r = 2# 履歴行列のパターン数
k = 4# テスト行列の列数
m = 2# テスト行列のパターン数
L = 4# ラグ
# 部分時系列の作成
X = np.array([y1[i:i+M] for i in range(len(y1)-M+1)]).T
T = len(X.T)
a = [0]*(M+n+L-1)
for t in range(M+n, T+M-L):
a.append(singular_spectrum_analysys(X, t, n, k, r, m, L))
fig, ax = plt.subplots(2,1,figsize=(8,4))
ax[0].plot(y1)
ax[1].plot(a)
アメリカのインフルエンザの患者数データで試す。データはアメリカの週ごとの患者数を表している。
2020年以降はご存じの通りcovid19の影響で患者数が激減しているが、それ以前では患者数が大きく増加していたことがわかる。
増加が始まった領域から毎年のように状況が変わっているが、変化スコアにもその様子が表れている。
import pandas as pd
data = pd.read_csv('ilidata.csv')
data['YEAR_WEEK'] = data['YEAR'].astype('str')+'-'+data['WEEK'].astype('str')
data = data.set_index(data['YEAR_WEEK'])
data = data['ILITOTAL']
M = 7# 窓幅
n = 5# 履歴行列の列数
r = 3# 履歴行列のパターン数
k = 5# テスト行列の列数
m = 3# テスト行列のパターン数
L = 52# ラグ
# 部分時系列の作成
X = np.array([data[i:i+M] for i in range(len(data)-M+1)]).T
T = len(X.T)
a = [0]*(M+n+L-1)
for t in range(M+n, T+M-L):
a.append(singular_spectrum_analysys(X, t, n, k, r, m, L))
fig, ax = plt.subplots(figsize=(10,4))
ax2 = ax.twinx()
ax.plot(data)
ax2.plot(data.index, a, color='red', alpha=.5)
ax.set_xticks(np.arange(0,701,50));
ax.tick_params(axis="x", labelsize=6)
ax.set_xlabel('week')
以上となります。
次回
疎構造学習による異常検知