書籍-入門機械学習による異常検知―Rによる実践ガイド-の7章
「状態空間モデルによる異常検知」で、カルマンフィルタによる異常検知を
実行すると、結果がどうなるのか気になっていました。

そこで、Pythonのnumpyとscipyを使ってゴリゴリ実装してみました。

#理論と実装
本稿では、結果の式だけを載せておきます。
式の意味や導出は、書籍をご覧ください。

書籍自体は大変分かりやすく、異常検知を始める人におススメです。
ただ、後半（7章は最後の方）は少し数式が多くなり、難しくなってきます。

##データロード
書籍でも使われている心拍数のデータを使います。

import urllib.request
import numpy as np
from numpy.linalg import svd, inv, matrix_rank
import matplotlib.pyplot as plt
from scipy.linalg import sqrtm

url =  "http://www.cs.ucr.edu/~eamonn/discords/qtdbsel102.txt"
f_name = "data.txt"
 
# ダウンロードを実行
urllib.request.urlretrieve(url, f_name)

# データ抽出
data = np.loadtxt(f_name, delimiter="\t")
data = data[2000:5000,1]

plt.plot(data)
plt.show()

このデータは長い時系列のデータとなっていますが、ここでは
時間：2000～5000（3000個）のデータを使います。

##データ作成

$X_p,X_f$を作ります。
$X_p$は時間的に前のデータを格納した行列、$X_f$は後のデータを格納した行列です。

w = 50 # 抽出窓の幅
m = 2  # zの次元
T = 1000 # 学習区間

x, x_after = [], []
    
for i in range(w, T):
    x.append(data[i-w: i])
    x_after.append(data[i: i+w])
    
x_p = np.array(x).T
x_f = np.array(x_after).T

##カルマンフィルタ
異常検知で使うカルマンフィルタは以下の式で定義されます。

\begin{align}
    p(x^{(t)}|z^{(t)}) &=N(x^{(t)}|Cz^{(t)},R) \\
    p(z^{(t)}|z^{(t-1)}) &=N(z^{(t)}|Az^{(t-1)},Q)
\end{align}

制御などで使うカルマンフィルタとはちょっと形が違います。

ただし、xは観測値を表すベクトル、zは内部状態を表すベクトル、Rは観測値のばらつきを表す共分散行列、Qは内部状態のばらつきを表す共分散行列です。状態変数系列のzを更新しながら、異常検知を行います。

行列A,C,Q,Rは部分空間同定法で求めます。

##部分空間同定法
カルマンフィルタで使うA,Q,C,Rを部分空間同定法で求めます。
まずは、Wを求めます。

\begin{align}
    S_{pf}=X_pX_f^T\\
    S_{pp}=X_pX_p^T\\
    S_{ff}=X_fX_p^T
\end{align}

W=S_{pp}^{-1/2}S_{pf}S_{ff}^{-1/2}

$X^{-1/2}$は$X^{1/2}$の逆行列を表します。

s_pf = x_p @ x_f.T
s_pp = x_p @ x_p.T
s_ff = x_f @ x_f.T

# 行列W
W = inv(sqrtm(s_pp)) @ s_pf @ inv(sqrtm(s_ff))

次に、Wについて特異値を求めます。
そして、特異値が大きい順にm本の左右特異ベクトルを求めます。

# 左右特異ベクトル
left, _, right = svd(W)

左特異ベクトルを特異値の大きい順に$\tilde{\alpha}^1,...,\tilde{\alpha}^m$、
右特異ベクトルを特異値の大きい順に$\tilde{\beta}^1,...,\tilde{\beta}^m$とし、$\tilde{\alpha}^i$、$\tilde{\beta}^i$を求めます。

\begin{align}
    \alpha^i=S_{pp}^{-1/2}\tilde{\alpha}^i\\
    \beta^i=S_{ff}^{-1/2}\tilde{\beta}^i\\
\end{align}

# 特異ベクトルm本抽出
alpha, beta = [], []
for i in range(m):
    alpha.append(inv(sqrtm(s_pp)) @ left[:,i])
    beta.append(inv(sqrtm(s_ff)) @ right[i,:])

次に$\hat{Z}_p$ , $\hat{Z}_f$を求めます。

\begin{align}
    \hat{Z}_p=[\alpha^1,...,\alpha^m]^TX_p\\
    \hat{Z}_f=[\beta^1,...,\beta^m]^TX_f\\
\end{align}

# zを算出
z_p = np.array(alpha) @ x_p
z_f = np.array(beta) @ x_f

$\omega$を検出窓の大きさとして、以下の値を求めます。

\begin{align}
    X &=[x^{(1)}, ...,x^{(T-w)}]\\
    Z &=[z^{(1)}, ...,z^{(T-w)}]\\
    Z_+ &=[z^{(2)}, ...,z^{(T-w)}]\\
    Z_- &=[z^{(1)}, ...,x^{(T-w-1)}]
\end{align}

最後にA,Q,C,Rを求めます。

\begin{align}
    A &=(Z_+Z_-^T)(Z_-Z_-^T)^{-1}\\
    Q &=\frac{1}{T-w-1}(Z_+-AZ_-)(Z_+-AZ_-)^T\\
    C &=(XZ^T)(ZZ^T)^{-1}\\
    R &=\frac{1}{T-w}(X-CZ)(X-CZ)^T
\end{align}

# A, Q, C, R
X = np.copy(x_p)
Z = np.copy(z_p)
Z_plus = z_p[:,1:]
Z_minus = z_p[:,:-1]

A = (Z_plus @ Z_minus.T) @ inv(Z_minus @ Z_minus.T)
Q = ((Z_plus - A @ Z_minus) @ (Z_plus - A @ Z_minus).T)/(T-w-1)
C = (X @ Z.T) @ inv(Z @ Z.T)
R = ((X - C @ Z) @ (X - C @ Z).T)/(T-w)

##漸化式

$z^{(t-1)}$を以下の式で置き換え、

p(z^{(t-1)}|X_{t-1}) = N(z^{(t-1)}|\mu_{t-1},V_{t-1})

次の漸化式で$\mu_t,Q_t$を推定します。

\begin{align}
    K_t &= Q_{t-1}C^T(R+CQ_{t-1}C^T)^{-1}\\
    \mu_t &= A\mu_{t-1}+K_t(x^{(t)}-CA\mu_{t-1})\\
    V_t &= (I_m-K_tC)Q_{t-1}\\
    Q_t &= Q+AV_tA^T
\end{align}

def update(A, Q, C, R, q, x, mu, m):
    k_t = q @ C.T @ inv(R + C @ q @ C.T)
    mu_t = A @ mu + k_t @ (x - C @ A @ mu)
    v_t = (np.eye(m) - k_t @ C) @ q
    q_t = q + A @ v_t @ A.T
    return mu_t, q_t

ただし、$\mu_0,Q_0$の初期値は適当に与えます。

q = np.eye(m)
mu = np.zeros(m)

##異常値
異常値aは以下の式で定義されます。

\begin{align}
    a(x^{(t)}) &=(x^{(t)}-CA\mu_{t-1})^T\Sigma_t^{-1}(x^{(t)}-CA\mu_{t-1})\\
   \Sigma_t &=R+CQ_{t-1}C^T
\end{align}

def anomaly_score(A, C, R, q, x, mu):
    sigma = R + C @ q @ C.T
    predict = C @ A @ mu
    score = (x - predict).T @ sigma @ (x - predict)
    return [score, predict, sigma]

#結果
最後に異常検知を行います。

result_score, result_predict, result_sigma = [], [], []
begin = 950
end = 1000

for i in range(2000):
    # t-1
    x = data[begin + i: end + i]
    mu, q = update(A, Q, C, R, q, x, mu, m)

    # t
    x = data[begin + i + 1: end + i + 1]
    result = anomaly_score(A, C, R, q, x, mu)

    result_score.append(result[0])
    result_predict.append(result[1])
    
    sigma_predict = np.array(result[2])
    sigma_predict = sigma_predict[-1, -1] - sigma_predict[-1,:-1] @ inv(sigma_predict[:-1,:-1]) @ sigma_predict[-1,:-1].T
    result_sigma.append(np.max([sigma_predict, 0]))

グラフを表示します。

result_score = np.array(result_score)
result_predict = np.array(result_predict)
result_sigma = np.array(result_sigma)

#boundary_upper = result_predict[:,-1] + 2*np.sqrt(result_sigma)
#boundary_lower = result_predict[:,-1] - 2*np.sqrt(result_sigma)

plt.figure(figsize=(12,12))
plt.subplot(2,1,1)
plt.plot(data[1000:3000], label="true")
#plt.fill_between(np.arange(2000), boundary_upper, boundary_lower, facecolor='y',alpha=0.3)
plt.plot(result_predict[:,-1], label ="predict")
#plt.xlim(1300,1400)
plt.legend()

plt.subplot(2,1,2)
plt.plot(result_score, label="score")
#plt.xlim(1300,1400)
plt.legend()
plt.show()

ちゃんと異常検知できましたが、K近傍法と同じような感度かなぁという印象。
やはり目的は違いますが、特異スペクトル変換（書籍参照）の方が感度良く検知します。

#参考文献