SVDによる変化検知能力 耐ノイズ 2

前回に引き続いて、変化点検知へのノイズの影響を追いかけてみる。

これが、地の部分の正弦波に対して46%のホワイトノイズを載せた例。

ノイズによる変化度が変化点での変化度の半分になっている。

続いて、40% と 10%の結果。

 

このように比べてみると、この数値実験の条件では地の正弦波に対して46%程度までのホワイトノイズであれば、位相変化を検知できると言えるだろう。

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

SVDによる変化検知能力 耐ノイズ 1

前回までの変化検知実験では理想的なデータを使っていたが、実際にはノイズが入る。前回の実験で変化検知していたのは、地の部分の信号の位相の変化だった。それではどれくらいのノイズが載ったら変化検知ができなくなるのだろう。

一連のテストデータでは、地の部分に p-p 0.02の正弦波がある。ここに、乱数でホワイトノイズの載せて変化度への影響を見る数値実験を試す。

まず、ノイズのない状態がこれ

地の部分に対して、48%のノイズを載せた結果は、こうなる。変化点と同じだけの変化度がノイズで生まれてしまって、変化検知ができていない。

ちなみに、このときのノイズ交じりのソースデータはこのとおり。目視では地の正弦波がしっかり見えてはいるのだが。

<pre>

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

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

n_amp = 0.01 * 0.1
bg_amp = 0.01
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)+bg_amp) * (1-math.cos(math.pi * p / 60)) + n_amp *(2*np.random.rand()-1)

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)

print (n_amp, n_amp / 0.01)
plt.plot(a[:pat_len])

</pre>

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

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

違いを生んだ、データの生成部分だけを示す。

3パルス目以降の offsetの値だけが前回と異なる

pat_len = 3000

pat = np.empty(3600)

for i in range (5):

    if (i<2):

        offset = 720*i

    else :

        offset = 720*i - 360

    for p in range (720):

        pat[offset+p] = (math.exp(-p/60)+0.01) * (1-math.cos(math.pi * p / 60))

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

前回は、より大きな模擬データを使って特異スペクトル変換による変化検知の実験をおこない、オフビートとなっている部分で変化度に顕著なピークを得られた。すばらしい...のか？

もう一度実験をしてみよう。前回とは、オフビートのタイミングが僅かに異なるデータを用意した。目視では違いがわかるだろうか。

その結果は、がらりと変わった。オフビートの所の変化度は前回の1/20程度で、オフビートのほうが低いとなった。目盛に注意して欲しい。

この結果は何を表しているのだろう。前回のデータと変化度のグラフを見ると、周期的パルスのタイミングの乱れを検知したと思ってしまいそうだ。
ところが、特異ベクトルで見ているのはパルス間の地の部分の波であり、周期的パルスはノイズ扱いされていると見てとれる。前回のデータはベースとなるサイン波が不連続となるタイミングでオフビートさせたデータ、今回はベースのサイン波が連続となるタイミングとした。それが、結果となって現れている。

SVDを使った変化検知では、どこの部分をみて定常/変化を見ているのかに注意しなければならない。

「異常検知と変化検知」[1]のp122 図9.5もしかり。p3の心電図データで赤い部分が他の収縮時波形と異なっているから変化度に現れたと考えるのは不適当だろう。






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

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

変化検知実験につかった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