時系列異常検知（折線の屈折点の特定）

Question

傾き違い二本直線が連接しています。その屈折点（連接点）を特定したかったです。時系列異常検知ライブラリ（BANPEI番兵)を使いました。
一次方程式で作ったデータに対して、成功しましたが（test1.py),　実際な実験データ（時系列データ）に応用したら、失敗しました。（test2.py)。お願いですが、もっと精度いい,有効な方法を教えていただきたいです。

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
import banpei # 番兵(banpei)を利用して、異常検知します。

################
# 変化点がある折線を作ります。
raw_data = []
nx =np.linspace(0,30)
for i in nx:
    if i <= 20:
        a = 0.5*(i-10)
    else:
        a = i-15
    raw_data.append (a)
################
# 内挿でデータ点数を増えます。
xp=np.arange(nx[0], nx[-1], 0.1)
rp=scipy.interpolate.splrep(nx, raw_data, s=0)
yp=scipy.interpolate.splev(xp,rp,der=0)
################
# 番兵(banpei)を利用して、異常検知します。
model = banpei.SST(w=8)
results = model.detect(yp)
###############
#グラフを表示します。
fig=plt.figure()
ax1 = fig.add_subplot(111)
ax2 = ax1.twinx()
ax1.plot(nx,raw_data) # 元データです。
ax2.plot(xp,results,color="red") #異常検知係数です。
plt.show()

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
import banpei # 番兵(banpei)を利用して、異常検知します。

################
# 変化点がある時系列データを入力する。

nx = [0,50,100,150,200,250,300,350,400,450,500,
        550,600,650,700,750,800,850,900,950,1000,
        1050,1100,1150,1200,1250,1300,1350,1400,1450,1500,
        1550,1600,1650,1700,1750,1800,1850,1900,1950,2000,
        2050,2100,2150,2200,2250,2300,2350,2400,2450,2500,2550,2600,
        2650,2700,2750,2800,2850,2900,2950,3000,3050,3100,3150,3200,
        3250,3300,3350,3400,3450,3500,3550,3600,3650,3700,3750,3800,
        3850,3900,3950,4000,4050,4100,4150,4200,4250,4300,4350,4400,
        4450,4500,4550,4600,4650,4700,4750,4800,4850,4900,4950,5000,
        5050,5100,5150,5200,5250]

raw_data = [373.766,373.796,373.826,373.86,373.891,373.923,373.954,373.985,
            374.015,374.047,374.078,374.11,374.14,374.173,374.204,374.236,
            374.267,374.299,374.33,374.362,374.393,374.425,374.458,374.49,
            374.522,374.554,374.587,374.62,374.651,374.681,374.711,374.743,
            374.776,374.807,374.837,374.867,374.9,374.932,374.964,374.996,
            375.028,375.06,375.092,375.124,375.155,375.187,375.219,375.25,
            375.281,375.313,375.346,375.378,375.41,375.441,375.472,375.504,
            375.537,375.568,375.6,375.632,375.663,375.694,375.724,375.757,
            375.788,375.821,375.851,375.882,375.914,375.944,375.974,376,
            376.022,376.041,376.059,376.078,376.098,376.114,376.128,376.143,
            376.158,376.172,376.19,376.205,376.22,376.234,376.245,376.257,
            376.27,376.282,376.297,376.309,376.323,376.34,376.352,376.366,
            376.381,376.396,376.408,376.425,376.438,376.451,376.464,376.476,
            376.487,376.501]


################
# 内挿でデータ点数を増えます。
xp=np.arange(nx[0], nx[-1], 1)
rp=scipy.interpolate.splrep(nx, raw_data, s=0)
yp=scipy.interpolate.splev(xp,rp,der=0)
################
# 番兵(banpei)を利用して、異常検知します。
model = banpei.SST(w=8)
results = model.detect(yp)
###############
#グラフを表示します。
fig=plt.figure()
ax1 = fig.add_subplot(111)
ax2 = ax1.twinx()

ax1.plot(nx,raw_data) # 元データです。
ax2.plot(xp,results,color="red",alpha = 0.5) #異常検知係数です。
plt.show()

Answer 1

時系列データが、ある点を境にしてそれぞれほぼ直線上に近い点であることが確実なら、以下のアプローチで解いた方が確実でしょう。

全体のデータを傾きの境の近くを除外したうえで、前半と後半に分離する。
ノイズを除去するために、raw_data[i+1]-raw_data[i-1]あるいは、raw_data[i+2]-raw_data[i-2]を調べ、変化する点を境とし、その前後の5個なり10個なりの点を除外すればできる。

前半の点の回帰直線、後半の点の回帰直線を求める。numpy.polyfitを使えばすぐに求めることができる。

前半の回帰直線と後半の回帰直線の交点を数値計算で求める。これは簡単。

ノイズ除去のストライドをどうするか、境界点の前後の丸くなっていうところをどれぐらい取り除けばよいかは試行錯誤が必要でしょう。

Comments

[Deng]

ご解答をありがとうございました。勉強になりました。助かりました。

[Deng]

Paulさんの解答になり、ソースを作りました。大成功です。
import numpy as np
import matplotlib.pyplot as plt

################
# 変化点がある時系列データを入力する。

nx = [0,50,100,150,200,250,300,350,400,450,500,
550,600,650,700,750,800,850,900,950,1000,
1050,1100,1150,1200,1250,1300,1350,1400,1450,1500,
1550,1600,1650,1700,1750,1800,1850,1900,1950,2000,
2050,2100,2150,2200,2250,2300,2350,2400,2450,2500,2550,2600,
2650,2700,2750,2800,2850,2900,2950,3000,3050,3100,3150,3200,
3250,3300,3350,3400,3450,3500,3550,3600,3650,3700,3750,3800,
3850,3900,3950,4000,4050,4100,4150,4200,4250,4300,4350,4400,
4450,4500,4550,4600,4650,4700,4750,4800,4850,4900,4950,5000,
5050,5100,5150,5200,5250]

raw_data = [373.766,373.796,373.826,373.86,373.891,373.923,373.954,373.985,
374.015,374.047,374.078,374.11,374.14,374.173,374.204,374.236,
374.267,374.299,374.33,374.362,374.393,374.425,374.458,374.49,
374.522,374.554,374.587,374.62,374.651,374.681,374.711,374.743,
374.776,374.807,374.837,374.867,374.9,374.932,374.964,374.996,
375.028,375.06,375.092,375.124,375.155,375.187,375.219,375.25,
375.281,375.313,375.346,375.378,375.41,375.441,375.472,375.504,
375.537,375.568,375.6,375.632,375.663,375.694,375.724,375.757,
375.788,375.821,375.851,375.882,375.914,375.944,375.974,376,
376.022,376.041,376.059,376.078,376.098,376.114,376.128,376.143,
376.158,376.172,376.19,376.205,376.22,376.234,376.245,376.257,
376.27,376.282,376.297,376.309,376.323,376.34,376.352,376.366,
376.381,376.396,376.408,376.425,376.438,376.451,376.464,376.476,
376.487,376.501]

#差分を取ります。
raw_dif = np.diff(raw_data, n=1)

#分界点を探します。
d= []
for i in range(2,len(raw_dif)-2):
if abs(raw_dif[i+2]-raw_dif[i])/abs(raw_dif[i]) > 0.15 :
d.append(nx.index(nx[i]))

#回帰直線
f = np.polyfit(nx[3:d[0]],raw_data[3:d[0]],1) #前半部分
f2 =np.polyfit(nx[d[10]:],raw_data[d[10]:],1) #後半部分

#連立方程式を解きます
m=np.matrix([[f[0],-1.0],[f2[0],-1.0]])
n=np.matrix([[-f[1]],[-f2[1]]])

x=np.linalg.solve(m,n) #解

x = np.array(x)

plt.plot(nx, raw_data)

plt.scatter(x[0],x[1],marker="x", c = "red")

plt.show()