特異スペクトル変換による変化検知 - プログラム

変化検知実験につかったPythonプログラムは以下のとおり。

定数のM, n, k, L, r, m の意味は文献[1] 9.4章 (p122 心電図データの変化検知)と同じ。

import numpy as np
import matplotlib.pyplot as plt
import math

def normV(v):
    return v / np.linalg.norm(v,2)

pat_len = 3000
pat = np.empty(3600)
for i in range (5):
    if (i<2):
        offset = 720*i
    else :
        offset = 720*i - 321
    for p in range (720):
        pat[offset+p] = (math.exp(-p/60)+0.01) * (1-math.cos(math.pi * p / 60))

plt.plot(pat[:pat_len])

M = 50
n = M // 2
k = M // 2
L = k // 2
r = 3
m = 2

a = np.zeros(pat_len)
for i in range (pat_len-L-k-M):
    X = np.empty((M,0))
    for j in range(n):
        X = np.hstack((X, np.array([normV(pat[i+j:i+j+M])]).T) )
    U,Sx,Vx = np.linalg.svd(X, full_matrices=True)
    Z = np.empty((M,0))
    for j in range(k):
        Z = np.hstack((Z, np.array([normV(pat[i+L+j:i+L+j+M])]).T) )
    Q,Sz,Vz = np.linalg.svd(Z, full_matrices=True)
    UQ = np.dot(U.T[0:r,:], Q.T[0:m,:].T)
    a[i+L+k+M-2]=1-np.linalg.norm(UQ,2)
  
plt.plot(a[:pat_len])

ちなみに実行環境は、OSが 64bit Windows 7, pythonは、anaconda3 の jupyter notebook.

[1] 井出 剛、杉山 将, 異常検知と変化検知, 講談社, 2017 ISBN978-4-06-152908-3

特異スペクトル変換による変化検知 - その3

特異スペクトル変換の基礎実験ができたので、もう少し大きなデータで実験をしてみる。

例えば、次のような長さ3000のデータについて解析してみよう。
履歴行列、テスト行列の取り方や特異値分解、スコアの計算は、参考文献のp122にある心電図の分析と同じパラメータとした。

対象のデータ


変化度


元データでパルスのリズムに異常のある近傍で著しいピークがあることがわかります。

パルスのリズムの変調を捉えられた...のでしょうか？


[1] 井出 剛、杉山 将, 異常検知と変化検知, 講談社, 2017 ISBN978-4-06-152908-3