角度データの変化検知を可視化するプロット

概要

角度データの解析を勉強しています。

作者: 清水邦夫
出版社/メーカー: 近代科学社
発売日: 2018/01/30
メディア: 単行本
この商品を含むブログを見る

この書籍中ではP95で、時系列角度データの変化点の検出について1節割かれています。しかし実際に可視化した際にどのような挙動を示すのか図表がなく分からなかったので、可視化してみました。

データの生成

いつも通りjupyterで実験しました。簡単にvon mises分布を2つ容易して、それを結合します。以下の例ですと系列店700の段階で傾向が変わったということですね。可視化した際にもなんとなくわかりますが、von mises分布自体のノイズが大きいので判断が難しいです。

import matplotlib.pyplot as plt
from scipy.stats import vonmises
%matplotlib inline

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
I = 1000
i1 = 700
i2 = 300
x = np.arange(0, I)
y1 = vonmises.rvs(kappa=0.8, loc=0, size=i1)
y2 = vonmises.rvs(kappa=0.4, loc=np.pi/2, size=i2)
y = np.concatenate([y1, y2])
ax.plot(x, y)

f:id:ajhjhaf:20180418201352p:plain

累積和プロット(Cumulative sum plot; CUSUM plot)

CUMSUMプロットでは、生成された角度を正弦と余弦に分解した時に、それぞれの成分の累積値をプロットします。同じ分布から生成された角度の累積値は一定の傾向を持ち続けますが、分布が変化した時にこの傾向が変わる、といった点から判断できるようです実際にノイズが取り除かれて、可視化結果も分かりやすくなっています。

fc = np.vectorize(lambda x: np.cos(x / 2 / np.pi))
cosy = fc(y)
c = np.cumsum(cosy)

fs = np.vectorize(lambda x: np.sin(x / 2 / np.pi))
siny = fs(y)
s = np.cumsum(siny)

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(c, s)
ax.set_xlabel("Cumulative cosine value")
ax.set_ylabel("Cumulative sine value")

f:id:ajhjhaf:20180418201929p:plain

累積平均方向プロット(cumulative mean directional plot)

このプロットでは、系列における時点jまでの累積平均方向をプロットします。分布が変化すると、平均の方向も変化するのでやはり判定できるといったものです。テキストではcos, sinの累積成分 $C_j$ , $S_j$ とそこから求められる合成ベクトルの長さ $R_j$ に対して

$\cos{\bar{\theta_j}} = \frac{C_j}{R_j}$

あるいは

$\sin{\bar{\theta_j}} = \frac{S_j}{R_j}$

を満たす $\theta_j$ をプロットすると書いてあります。ですが、これを求める際には

$\tan{\bar{\theta_j}} = \frac{\sin{\bar{\theta_j}}}{\cos{\bar{\theta_j}}} = \frac{\frac{S_j}{R}}{\frac{C_j}{R}} = \frac{S_j}{C_j}$

なので、別に $R_j$ は必要ありません。また系列の最初の方ではnが少なくて値が安定しません。なので、最初100を削ってスケールを合わせたものも出してみました。

y_bar = np.arctan(s/c)

fig = plt.figure()

ax = fig.add_subplot(2,1,1)
ax.plot(x, y_bar)
ax.set_xlabel("series")
ax.set_ylabel("average direction")

ax = fig.add_subplot(2,1,2)
ax.plot(x[100:], y_bar[100:])
ax.set_xlabel("series")
ax.set_ylabel("average direction")
ax.set_xlim(0, I)

f:id:ajhjhaf:20180418202829p:plain

ランク累積和プロット(rank CUSUM plot)

ランク化して計算します。今回扱うプロットの中では一番分かりやすいですね。元々の分布がなんだったかにも影響しそうですが、計算が面倒くさいだけあっていい感じです。

from scipy.stats import rankdata
gamma = 2 * np.pi * rankdata(y) / y.size
U = np.sqrt(2 / y.size) * np.cumsum(np.cos(gamma))
V = np.sqrt(2 / y.size) * np.cumsum(np.sin(gamma))
B = np.sqrt(np.power(U, 2) + np.power(V, 2))

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(x / I, B)
ax.set_xlabel("relative series")
ax.set_ylabel("B")

f:id:ajhjhaf:20180418203106p:plain

累積円周分散プロット(cumulative circular variance plot)

テキスト中では

各j=1,...,nに対して、円周分散1-R_jを計算し…

と書いてありますが、正しくは合成ベクトルの大きさ $R_j$ ではなく、平均合成ベクトル長 $\bar{R_j}$ を使って計算する必要があります。傾向が累積平均方向ベクトルと変わらないですね。

r = np.sqrt(np.power(c, 2) + np.power(s, 2))
r_bar = r / np.arange(1, I+1)
v = 1 - r_bar

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(x, v)
ax.set_xlabel("series")
ax.set_ylabel("cumulative circular variance")

f:id:ajhjhaf:20180418203433p:plain