DEIM Forum 2016 G1-6

自動パターン検出のためのストリームアルゴリズム
川畑光希†

松原 靖子††

櫻井保志††† 熊本大学工学部情報電気電子工学科†† 熊本大学大学院自然科学研究科情報電気電子工学専攻
E-mail: †kouki@dm.cs.kumamoto-u.ac.jp, ††{yasuko,yasushi}@cs.kumamoto-u.ac.jp
あらまし

近年，データストリーム処理に関する研究が盛んに行なわれている．本論文は，様々な時系列パターンを含むデータストリームの中から，(a) 重要な特徴を自動的に発見，それらの情報を統計的に要約し，(b) 得られた特徴に類似した部分シーケンスを検出することを目的とする．時系列ストリームは，センサデータや Web アクセス履歴
等，様々なアプリケーションにおいて大量に生成されており，これらの大規模な時系列ストリームの中から典型的なパターンや異常値を発見することは非常に重要な課題である．しかし，ネットワーク分析，センサ監視等，データ量が多く緊急性が要求されるような近年のアプリケーションでは，すべてのデータを蓄積してから処理することが困難である．本研究では，このような問題を解決するアルゴリズムを提案する．実データを用いた実験では，提案手法がデータストリームから有用なパターンを正確に発見，かつそれに類似した部分シーケンスを検出できることを確認し，オンライン処理であるにもかかわらず，従来の手法と同等の精度を保ち，計算時間について性能向上を達成していることを明らかにした．
キーワード

ストリームデータ処理，特徴自動抽出

1. は じ め に

時系列ストリームは，センサデータや Web アクセス履歴等，なっている．したがって，すべてのデータをメモリ空間に格納することなく，高速に処理することができるストリームアルゴ
リズムは非常に重要である．

様々なアプリケーションにおいて大量に生成されている．一般

1. 2 本論文の貢献に，実際に生成されるデータストリームは複数のトレンドやパ

S TREAM S COPE は以下の特長がある．
（ 1 ） データストリームから時系列パターン（レジーム） のターンを持つことが多い．また，ネットワーク通信のモニタリ
ングシステムから発生するデータについては，正常と異常のパターンが考えられる．本論文では，このような時系列ストリー
ムを対象とした自動パターン検出の手法として S TREAM S COPEを提案する．S TREAM S COPE は監視するデータストリームに関する事前知識を必要とせず，その中からユーザの直感にあった
特徴を自動的に抽出し，モデルパラメータとしてその情報を圧縮する．このとき，S TREAM S COPE は刻一刻と生成されるストリームが，(1) 抽出した特徴の中のどれに該当するかをリアルタイムに把握し，(2) それらの特徴を表すモデルを最新の情報に更新しながら (3) 特徴間の変化点を発見するだけでなく，それらとは異なる (4) 新たなパターンを発見することが可能である．
より具体的には，データストリーム X が与えられたとき，(a) X の中のパターンの変化点を発見し，部分シーケンス集合（セグメント）に分割し，(b) それらのセグメントをグループ化
し，類似時系列パターン（本論文では「レジーム (regime)」と呼ぶ） を発見する．さらに重要な点として，これらの処理は(c) 高速かつ，リアルタイムで行う．
1. 1 自動特徴抽出とストリーム処理の重要性
時系列データを対象とした研究課題は，数多く存在する．パ
ターン発見 [20], [14], [19], [24]，情報要約 [5], [12], [13]，クラス
タリング [10]，セグメンテーション [9], [25] や類似シーケンス
探索 [6], [17], [21] 等は重要な課題である．しかし，これらの先行研究はすべてのデータを蓄積してから処理することが前提と個数と種類をオンライン処理によって把握し，それぞれの適切な変化点をリアルタイムに発見する．さらに S TREAM S COPE
は，パターン変化点の最適解の検出を保証する．
（ 2 ） 4. 章で述べるモデルにより，ユーザの直感に合致した
時系列パターンの抽出を行う．
（ 3 ） S TREAM S COPE は一度検出した時系列パターンの部分
シーケンスに該当するデータを必要とせず，適切なレジームの数，変化点の数を，自動的に発見することができる．
（ 4 ） 従来のオフライン手法と比較して，クラスタリングの精度を落とさず計算時間について性能向上を達成した．

2. 関 連 研 究
関連研究は以下の 3 つに分類される
パターン検出．時系列データの解析に関する研究は様々な分野で進められている．[1], [3], [16]．自己回帰モデル (AR: autoregres-sive model), 線形動的システム (LDS: linear dynamical systems),
カルマンフィルタ (KF: Kalman filters) は代表的な技術であり，これらに基づく時系列の解析と予測手法が数多く提案されている [10], [22]．著者らの先行研究 [12], [13], [19], [20] では，主にセンサデータや Web 関連のデータを用いて，大規模時系列デー
タ集合のための探索，モデル推定，予測手法を提案している．
時系列データからの情報抽出に関しては，Li らが文献 [9] において，欠損を含む大規模時系列シーケンス集合のためのアルゴリズムである DynaMMo を提案している．DynaMMo は LDS

検出及び，レジームの発見を高速かつリアルタイムで行うことに基づき，時系列データのパターンを発見し，シーケンスのセである．本論文で取り組む問題を以下のように定義する．

グメント化の能力を持つ．また，データストリームにおけるパ［ 問題 1］ 多次元時系列ストリーム X が与えられたとき，X

ターン検出の研究としては，Lp 距離に基づくストリームのトレンド検出や相関検出，予測を始めとする様々な手法が提案さを表現するような以下の情報を抽出する．
（ 1 ） セグメントの総数 m と各セグメントの位置：れている． SPIRIT はデータストリームから相関とトレンドにS = {s1 , . . . , sm }相当する隠れ値を検出する問題に取り組んだものである [16]．

（ 2 ） レジームの総数 r とセグメントメンバーシップ：

BRAID はデータストリーム間の遅延相関を検出するための手

F = {f1 , . . . , fm }法である [20]．

（ 3 ） r 個のレジームを表現するモデルのパラメータ集合：

確率モデル．隠れマルコフモデル (HMM: Hidden Markov model)Θ = {θ1 , . . . , θr , ∆r×r }これらの情報はコスト関数（式 (5)）を最小化するものを選ぶ．

は音声認識を含む様々な分野において，時系列解析手法として広く利用されている．HMM に基づく大規模時系列シーケンスのための研究として，文献 [6] では大規模 HMM データ集合のための高速探索アルゴリズムが扱われている．Wang ら [25] は
文献 [7] を改良し，pHMM (pattern-based hidden Markov model)を提案している．pHMM は時系列のセグメント化とクラスタリングのための動的モデルであり，時系列シーケンスをマルコフモデルに基づいて線形のセグメントに分割する能力をもつ．
階層的確率モデルとして，Fine ら [4] は階層的 HMM (HHMM:

hierarchical HMM) を提案し，Fox ら [5] はベータ過程に基づくモデルとして，BP-AR-HMM (beta process autoregressive HMM)を提案している．
情報抽出とクラスタリング． 情報抽出とクラスタリングの手法については，CLARANS [15], BIRCH [26], TRACLUS [8] を含
め，様々なものが提案されている．文献 [2], [23] においては，

MDL の概念を用いて情報要約とクラスタリング問題を扱っている．また，著者らは文献 [11] において，大規模時系列データを対象とした自動特徴抽出手法を提案している．

3. 問 題 設 定
ここでは本論文で必要な概念について定義を行う．データ
ストリーム X は，x1 , x2 , . . . , xn , . . . の値からなる半無限長のシーケンスである．xn は最も新しい値であり，時刻が進むごとに n は増加する．xt は時刻 t における d 次元ベクトルとする．このようなデータストリーム X が与えられたとき，本研
究は X を m 個のセグメント集合 S= {s1 , . . . , sm } に分割することを目的とする．si は i 番目のセグメントの開始点，終了点で構成され（つまり，si = {ts , te }），各セグメントは重複がないものとする．本研究ではさらに，発見したセグメント集合を
類似セグメントのグループ（レジーム: regime ）に分類する．
［定義 1］（レジーム） r を最適なセグメントグループの個数とする．それぞれのセグメント s はセグメントグループの１つに割り当てられる．これらグループをレジーム (regime) と呼び，それぞれのレジームは統計モデル θi (i = 1, . . . , r) として表現される．
［定義 2］（セグメントメンバーシップ） F = {f1 , . . . , fm } を，

m 個の整数列とし，fi を i 番目のセグメントが所属するレジー< fi =< r)．
ムの番号とする (1 =本研究の目的は，長さ n の多次元時系列ストリームが与えられたときに，その部分シーケンスのセグメント化と分割位置の本論文では，レジームを表現するモデルパラメータ集合 Θを，r 個の隠れマルコフモデル (HMM: hidden Markov model)，{θ1 , . . . , θr }，として表現する．（注 1） さらに，レジーム間の関係性を表現するため，レジーム遷移行列 ∆r×r を使用する．レ
ジーム遷移行列（定義 4）とコスト関数（式 (5)）についての詳細は，4. 章で説明する．
問 題 1 で 示 し た 通 り，本 論 文 の 目 的 は ，X の 特 徴 を 抽出し，すべての時系列パターンを表現するパラメータ集合{m, r, S, Θ, F } を発見することである．本論文では，この全パラメータ集合を候補解 C と呼ぶ．
［定義 3］ X を表現する全パラメータ集合 C = {m, r, S, Θ, F }を候補解と呼ぶ．候補解 C は，セグメント集合，各セグメントのレジームへの割当て，レジームを表現する確率モデル，これらすべてを表現する．
結論として，本論文の目的は最適な解 C を発見することである．
ここで非常に重要な課題は，(a) どのようにセグメントおよびレジームの数を推定するか，(b) どのようにレジームを表現し，セグメントの割当てを行うかである．本研究では，リアルタイ
ム処理によって最適解を求めるための手法を提案する．

4. 特徴抽出とデータ圧縮
本章では，問題 1 を解決するために使用するモデルについて説明する．
本研究では，複数のレジーム間の時系列パターンとその遷移を表現するために，多層的なモデルを使用する．具体的には，
隠れマルコフの状態遷移をレジームにグループ化し，階層的な時系列パターンの遷移を表現する．つまり，隠れマルコフの状態 (state) に対して，レジームを上位層の状態 (super-state) と考え，遷移確率が以下のように定義される．
［定義 4］（レジーム遷移行列） ∆r×r を r 個のレジーム群の遷移行列と呼ぶ．ここで，要素 δij ∈ ∆ は i 番目のレジームから

j 番目のレジームへの遷移確率を示す．
行 列 ∆ 内 部 の 要 素 δi,j は 確 率 を 表 し ，0 <= δij <= 1,
δ = 1 という条件を持つ．よって本論文で用いるモデj ij

P

ルは r 個のレジーム集合 Θ = {θ1 , . . . , θr , ∆r×r } で表現され，

θi は i 番目のレジームのモデルパラメータを表現する．ここで，（注 1）：本論文で提案する枠組みは，HMM 以外の時系列モデルに適用することも可能である．

表1

記号
シーケンス
n
d
X
xt
X[ts : te ]
セグメント
m
S
F
レジーム
r
Θ

定義

時系列の長さ時系列の次元数d 次元の時系列シーケンス/ストリーム
X の t 番目の値/要素 (t = 1, . . . , n)
ts から te までの X の部分シーケンス

X に含まれるセグメントの総数
X に含まれるセグメント集合：S = {s1 , . . . , sm }セグメントメンバーシップ：F = {f1 , . . . , fm }

θi
ki
∆r×r

図1

コスト関数
C
CostM (Θ)
CostC (X|Θ)
CostT (X; C)

本論文で使用するモデル Θ （ここでは r = 3）．

θi は HMM に基づき，初期確率，遷移確率，出力確率の三つ組
で次のように表現される：θi = {πi , Ai , Bi }. （注 2）
次に，最小記述長 (MDL: minimum description length) の考え方に基づき，大規模時系列データを表現する符号化スキームについて説明する．直感的には，データが与えられたときのモデルのよさは次の式で表現できる：CostT = Cost(M)+Cost(X|M)．
ここで，Cost(M) はモデル M を表現するためのコストを示
し，Cost(X|M) は，M が与えられたときのデータ X の符号

log∗ (d) ビット（注 3）•セグメントとレジームの個数 m, r：log∗ (m) + log∗ (r)•各セグメントのレジームへの割当て（セグメントメン

バーシップ）：m log(r) ビット
各セグメントの長さ s：•

r 個のレジームのモデルパラメータ集合：CostM (Θ)

CostM (Θ) =

r
X

i=1

log∗ |si | ビット

CostM (θi ) + CostM (∆)

候補解：C = {m, r, S, Θ, F }Θ のモデル表現コスト
Θ による X の符号化コスト
C による X の総コスト

よって，シーケンス X と r 個のレジームのモデルパラメー
タ集合 Θ が与えられたとき，データ圧縮のためのコストの総数は次の通りである．

CostC (X|Θ) =m
X

CostC (X[si ]|Θ)

i=1

=

m
X

− ln(δvu · (δuu )|si |−1 · P (X[si ]|θu ))

(4)

i=1

ここで，i と (i − 1) 番目のセグメントはそれぞれ u と v 番目のレジームに所属し，fi = u, fi−1 = v, f0 = f1 とする．X[si ]はセグメント si の部分シーケンス，P (X[si ]|θu ) はセグメント

si の尤度，θu はセグメント si が所属するレジームである．
まとめると，候補解 C = {m, r, S, Θ, F } が与えられたとき

•

Pm−1

X に含まれるレジームの総数
r 個のレジームのモデルパラメータ集合：Θ = {θ1 , . . . , θr , ∆r×r }i 番目の レジームのモデルパラメータθi の状態数
レジーム遷移行列：∆ = {δij }r
i,j=1

ここで，P (X|θ) は X の尤度を示す．

化のコストを示す．表現コスト Cost(M) は以下の要素から構成される．
• 多次元シーケンスデータの長さ n と次元数 d：log∗ (n) +主な記号と定義．

(1)i=1

単一のレジームのモデル θ は，状態数 k（log (k)）と確率モデの X の符号長は次のように表現される．

CostT (X; C) = CostT (X; m, r, S, Θ, F )= log∗ (n) + log∗ (d) + log∗ (m) + log∗ (r) + m log(r)∗ル (θ = {π, A, B}) の表現コストが必要となる．まとめると，

m−1

+

X

log∗ |si | + CostM (Θ) + CostC (X|Θ)

(5)

i=1

CostM (θ) = log∗ (k) + cF · (k + k2 + 2kd)

(2)

ここで，cF は浮動小数点のコストを示す．（注 4） 同様にして，レジーム遷移行列には，CostM (∆) = cF · r2 のコストを要する．
本論文では確率モデルを用いてシーケンス X の時系列パターンを表現する．そこで，推定したモデルが X を正しく表現しているかを判断する指標が必要である．ハフマン符号を用いた情報圧縮では，モデル θ が与えられた際の X の符号化コストを負の対数尤度を用いて次のように表現することができる．

CostC (X|θ) = log2

1
= − ln P (X|θ)
P (X|θ)(3)（注 2）：本論文では出力確率 B に多次元ガウス分布を仮定する．これにより多次元
ベクトルのシーケンスを確率モデルで表現する（つまり B = {N (µi , σi2 )}k
．
i=1 ）
（注 3）
：こ こ で ， log∗ は 整 数 の ユ ニ バ ー サ ル 符 号 長 を 表 す：log∗ (x) ≈log2 (x) + log2 log2 (x) + . . .． [18].
（注 4）：本論文では 4 × 8 ビットとする．

本論文の最終的な目標は，上記のコスト関数を最小化するようなセグメントおよびレジーム集合を時系列ストリームの中から発見することである．

5. 最適化アルゴリズム
前章では，候補解 C = {m, r, S, Θ, F } が与えられた上でストリーム X を表現するためのコスト関数である式 (5) について述べた．続いて本章では，式 (5) に基づき，最適な解 C を発見するためのストリームアルゴリズムである S TREAM S COPE を提案する．

5. 1 概要本研究では，刻々と生成されるデータストリームから，前章で述べたコストモデルに基づき，セグメントおよびレジームの個数を自動的に選択する．直感的には，ストリームを監視しな複数のレジーム間における変化点を検出する．そして，時刻 t
における入力 xt と，モデルが与えられた上での符号化コスト

CostC (X|Θ) = − ln P (X|Θ) をインクリメンタルに計算するためのアルゴリズムを提案する．

5. 2. 1 パターン変化点の検出=
本 研 究 で は ，入 力 xt と 2 つ 以 上 の レ ジ ー ム θ1<=r=i<{πi , Ai , Bi }, およびレジーム遷移行列 ∆r×r が与えられたとき，xn の尤度 P (xn |Θ) は次のようにインクリメンタルに計算される．
P (xn |Θ) = max { max1<u<ki {pi;u (n)}
= =
1<
=i<
=r

図2

S EGMENTA SSIGNMENT の様子．X と θ1 , θ2 , . . . , θr が与えられたとき，インクリメンタルな走査でレジームの変化点を検出し，

pi;u (t) = max

セグメントの割り当てを行う．

// regime θi }
δ · maxv {pj;v (t − 1)} · πi;u · bi;u (xt ) ji// regime switch from θj to θi δii · maxw {pi;w (t − 1) · ai;wu } · bi;u (xt )

(6)

(7)

// staying at regime θipi;u (t) は，時刻 t におけるレジーム θi の状態 u の確率の最大
がらその部分シーケンスを適切なレジームに割り当てていく．
また，どのレジームにも当てはまらない場合には新たにレジー
ムを生成する．そして，2 つの処理を繰り返しながら，時刻 nでの解 C を最適なものに更新していく．よって，提案手法に必値を示し，式 (7) の上段は，レジーム間の切り替え（θj から

θi ）の確率，下段はレジーム θi の内部の状態遷移の確率を示す．また，δji はレジーム j からレジーム i への遷移確率を示
し，maxv {pj;v(t−1) } は時刻 t − 1 における θj 内の確率の最大要な処理は以下の 2 つである．
（ 1 ） 入力 xn に対する候補解 C 中のモデルとの尤度をイン

値を示す．π1,i , b1;i (xt ), a1;ji はそれぞれ θ1 内における状態 i

クリメンタルに計算し，レジーム変化点を検出する．変化点に

式 (7) より，時刻 t における遷移確率 pi;u (t) は，時刻 t − 1

より分割したセグメントを適切なレジームに割り当てる．

の初期確率，出力確率，状態 j から状態 i への遷移確率を示す．
の遷移確率から計算される．したがって，ストリーム処理にお（ 2 ） 候補解 C 中にはない新たな時系列パターンを検出し，いては，時刻 t − 1 における各レジームの各状態の遷移確率を

適切なレジームを作成する．
本論文では，(1) についてのアルゴリズムである S EGMENTA S -保持しながら，累積的に尤度を計算する．そうすることで，時SIGNMENT ついて説明し，最後に提案手法である S TREAM S COPEてのみであり，探索漏れがないことを保証しながらも，過去のについて説明する．

5. 2 S EGMENTA SSIGNMENT
ここではストリーム X からレジーム変化点を検出し，各セグメントを適切なレジームに割り当てるためのアルゴリズ
ムを提案する．今，ストリーム X のうち，現時刻のデータであるベクトル xt と，2 つ以上のレジームのモデルパラメータ{θ1 , θ2 , . . . , θr , ∆r×r } が与えられていると考える．このとき，
S EGMENTA SSIGNMENT は レジームのモデルパラメータに基づき，X のパターンの変化点（つまりセグメントの分割位置）を検出することができる．ここで重要な点として，提案アルゴリ
ズムは探索漏れがないことを保証しながらも，一度処理を行ったデータを遡ることなくパターン変化点を検出することができる．
図 2 は時系列ストリームからある一つのレジーム変化点を検出するときの S EGMENTA SSIGNMENT の処理の様子を示している．これは青のレジーム θ1 から緑のレジーム θ3 へ切り替わる
例である．
より具体的には，X が与えられたとき，コスト関数（式 (4)）
を最小とするような，レジーム変化点をリアルタイムに発見し，それをもとに分割した部分シーケンスを 1 つのセグメント (S)として適切なレジームに割り当てることを考える．
レジーム遷移行列 ∆ はすべてのレジームの組み合わせに対して遷移確率を持っており，本論文ではその遷移行列を用いて

刻 t における尤度の計算は図 2 中の 1 つのトレリス構造についデータを遡ることなく高速に尤度を計算することができる．
続いて，検出したレジーム変化点の候補集合について述べる．
L = {l1 , l2 , . . . , li , . . .} を変化点集合とし，li は i 番目の変化点を表す．変化点 l は変化点の時刻と遷移先のレジーム ID から構成される．本論文では，これらの候補集合を各レジームの各状態（つまり，式 (8) において， 1 <=i<= r, 1 <=u<= ki ）に対して保持し，その中から最適な変化点集合を発見する．
Li;u (t) =

Lj;v (t − 1) ∪ {l} // switch from θj to θi
Li;u (t − 1)// staying at regime θi(8)

ここで Li;u (t) は時刻 t における θi の状態 u の変化点集合である．これらは式 (7) の尤度計算に基づき更新される．ある時刻 tにおいて，レジーム j からレジーム i に切り替わった場合，変化点 l = (t, i) を変化点の候補として集合に加える．S EGMEN TA SSIGNMENT はこのような候補集合から尤度 P (xn |Θ) を最大化するような最適解 Lbest を選出する．
［補助定理 1］ モデル Θ = {θ1 , . . . , θr , ∆} および長さ n のシーケンス X[1 : n] が与えられたとき，尤度 P (xn |Θ) を最大にする変化点集合は最適である．
［ 証明 1］ S EGMENTA SSIGNMENT は時刻 t − 1 から時刻 t にかけての１つのトレリス構造について尤度を累積的に更新する．このとき，レジーム間の変化点の候補は集合としてすべて記録される．よって，各時刻において，各レジームの各状態のViterbi パスがすべて記録されており，尤度 P (xn |Θ) を最大にする変化点集合は最適 Viterbi パスである．

✷

［補助定理 2］ モデル Θ = {θ1 , . . . , θr , ∆} および，シーケン示す．S EGMENTA SSIGNMENT は時刻 t における遷移確率 p，変ス X[1 : t + γ] が与えられたとき，時刻 1 から t までの部分パ

化点集合 L を更新する．

スから生成されるセグメントは時刻 t + γ において最適である．

次に，レジーム変化点を検出した場合，つまり，Lbest に変［ 証明 2］ 尤度 P (xt+γ |Θ) は時刻 t + γ において確率の最大化点候補が存在するとき (|Lbest | > 0) ，式 (9) より γ を計算し，値を示す．したがって，補助定理 1 より，時刻 1 から t + γ ま

変化点集合が γ 保証されるのを待つ．

での変化点集合は最適である．よって，時刻 1 から t における部分パスは時刻 t + γ まで最適である．

✷

本研究では，時刻 1 から t における部分パスが補助定理 2 を満たすことを，γ 保証と呼ぶ．

γ 保証が完了した場合，Lbest の情報をもとに，部分シーケ
ンスをセグメントとしてレジームに割り当てる．そして，すべての変化点集合を初期化し，カレントウィンドウの開始点 tc に遷移したレジームの開始点を代入する．カレントウィンドウの［定義 5］ 時刻 t において，P (xt |Θ) を出力するレジームがレ

定義については次節で説明する．

ジーム θi への変化した場合の γ を以下のように定義する．

［ 例 1］ ここでは図 2 を用いてアルゴリズムの具体的な例を示

γ∝

l

1
1 − δii

m

(1 <
=i<
= r)

(9)

す．S EGMENTA SSIGNMENT は各時刻において，時刻 t − 1 から
現時刻 t についての尤度を計算する．時刻 1 において，θ1 内の確率 p1;3 (1) は確率の最大値を持つ．時刻 2 において S EGMEN -すなわち，γ は遷移先のセグメント長に基づくものである．変化点の候補を検出した場合，ただちにその変化点によりセグメ
ントの割り当てを行わず，γ 保証された変化点集合をもとにセグメントをレジームに割り当てる．直感的には，あるレジーム

θi への変化点の候補を検出したとき，過去にレジーム θi に割り当てられたセグメントの長さを考慮するために，γ 時刻未来において尤もらしいことが保証された変化点を用いてセグメン
トの割り当てを行う．

TA SSIGNMENT は p1;3 (1) から p2;3 (2) への変化点の候補を発見する（図では θ1 から θ2 への矢印）．同時に，L2;3 (2) = {(2, 2)}を候補集合として更新する．ここで，(2, 2) は時刻 2 において，レジーム 2 へ遷移したことを表している．さらに，式 (9) にしたがって，γ を計算する．同様にして，2 番目の候補点として

L3;1 (4) = {(4, 3)} を保持する. 初めて変化点の候補を検出した
時刻 t = 2 から γ 時刻後，アルゴリズムは p3;1 (2 + γ) が確率の最大値を持つことを明らかにし，L3;1 (2 + γ) = {(4, 3)} を最適解として出力する．

5. 2. 2 アルゴリズム

5. 3 S TREAM S COPE
Algorithm 1 S EGMENTA SSIGNMENT (xt , θ1 , . . . , θr , ∆)本論文の最終目標は，大規模時系列ストリームの中から，直1: Input: Vector xt , model parameters of all regimes {θ1 , . . . , θr , ∆}感にあった時系列パターンを発見し，それらをグループ化する

2: Output: (a) Number of segments, m

ことである．ここで解決すべき問題は，刻一刻と生成される時

3:(b) Segment set S = {s1 , . . . , sm }4:(c) Segment membership F = {f1 , . . . , fm }系列パターンを既知レジームにグループ化しながら，それらと5:(d) Best cut point set Lbest = {l1 , l2 , . . .}は異なる新たな時系列パターンをどのようにして検出するかと

6: /* Compute pi;j (t) */いうことである．

7: for i = 1 to r do
8:Compute pi;j (t) for state j = 1, . . . , ki ; /* Equation (7) */9:

Update Li;j (t) for state j = 1, . . . , ki ; /* Equation (8) */10: end for
11: Choose the best cut-point set Lbest ;/* Equation (6) */12: if |Lbest | > 0 and γ == 0 then13:

Compute γ; /* Equation (9) */14: else15:

if γ > 0 then16:γ = γ − 1;

17:

else

5. 3. 1 アルゴリズム
大規模時系列ストリームの中から自動的に時系列パターンを取り出す手法として S TREAM S COPE を提案する．S TREAM S COPE
はコスト関数である式 (5) をもとに，与えられたストリームの中からこれまでに検出したパターンと同じパターンを検出しな
がら，新たな時系列パターンを検出した場合，自動的に新たなレジームを生成する．
ストリーム処理の場合，すべてのデータを蓄積して使用する

18:/* Assign guaranteed segments into optimal regime */ことができない．そこで，S TREAM S COPE では新たな時系列パ

19:

ts = tc ;20:

for each cut point l = (te , i) in Lbest do

ターンを発見するために必要な部分シーケンスのみを保持し21:

Create a new segment s = {ts , te };22:

Add s into S; m = m + 1; fm = i;ウィンドウと呼ぶ．

23:

ts = te ;［定義 6］ (カレントウィンドウ) 候補解 C のうち，最も新しい24:

end for25:

tc = t e ;initialize L;

26:
27:

end if

28: end if
29: return {m, S, F ,Lbest };ながら処理を行う．この部分シーケンスを本研究ではカレント

セグメント sm = {ts , te } の開始点 ts をカレントウィンドウの開始点 tc とし，tc から現時刻 t によって切り出されるストリー
ム X の部分シーケンス X[tc : t] をカレントウィンドウとする．
アルゴリズム 2 は S TREAM S COPE の処理の流れを示している．S TREAM S COPE は毎時刻 S EGMENTA SSIGNMENT によって既知のレジームによるパターンの変化点を検出する．

アルゴリズム 1 に S EGMENTA SSIGNMENT の具体的な処理を続いて，未知のレジーム，つまり新たな時系列パターンをを推定する．続いて，すべてのモデルのペア {θs1 , θs2 } に対し，Algorithm 2 S TREAM S COPE (xt )1: Input: a new vector xt at time tick t

符号化コストを計算し，最も適切なペア {θ1 , θ2 } を初期モデル

2: Output: Complete set of parameters C, i.e.,として選出する．

3:(a) Number of segments, m4:(b) Number of regimes, r5:(c) Segment set S = {s1 , . . . , sm }6:(d) Model parameters of regimes Θ = {θ1 , . . . , θr ; ∆}7:(e) Segment membership F = {f1 , . . . , fm }{θ1 , θ2 } =CostC (X|θs1 , θs2 )arg min(10)

θs1 ,θs2 |s1 ,s2 ∈X

ここで，X = {s1 , s2 , . . .} は，カレントウィンドウから取り出

8: {m, S, F , L}=S EGMENTA SSIGNMENT (xt , θ1 , . . . , θr , ∆);したサンプルの集合を示す．

9: {m′ , S ′ , θ ′ }=R EGIME G ENERATION(X[tc : t]);

モデルパラメータの推定． HMM のモデルパラメータの推定手

10: if S ′ =| φ then
11:

Compute CostSW and CostCU T/* Equation (11), (12) */法である Baum-Welch アルゴリズムは，モデル θ に対し，隠れ12:

if CostSW > CostCU T or |L| == 0 then状態の数 k を与える必要がある．しかし，この k を手動で設定

13:

S = S ∪ S ′ ; Θ = Θ ∪ θ ′ ; r = r + 1;

するのは非常に難しい．もし k の値を小さくすれば，データの14:

Update tc ;15:

Update ∆r×r ; /* Equation (13) */fi = r16:′

表現能力が低くなり，適切なセグメントおよびレジームを求め′(i = m + 1, . . . , m ); m = m + m ;end if17:

18: else
19:

i = fm ;20:

Update θi ;21: end if

検出するため，R EGIME G ENERATION によってカレントウィン
ドウを分割することを試みる．直感的には，1 つのモデルパラ
メータ θ では表現できない時系列パターンを発見したとき，それは新たな時系列パターンであり，それを表現するためのモデル θ ，およびレジームを新たに作成する．

R EGIME G ENERATION． カレントウィンドウ X[tc : t] が与えられたとき，R EGIME G ENERATION は式 (5) を用いてカレント
ウィンドウ内における表現コストを最小にするようなモデルパ
ラメータの推定を行う．具体的には以下の 2 つのステップから構成される反復処理によって，モデルパラメータの推定を行う．

•

ステップ 1: 定理 1 より，カレントウィンドウの符号化コストが最小となるレジーム変化点を検出し，セグメント集合
を 2 つのグループ {S1 , S2 } に分割する．

•

ステップ 2: ステップ 1 で得られたセグメント集合に基ることが困難となる．一方で，もし k を大幅に上げてしまうと，オーバーフィッティングを招く．そこで本研究では，隠れ状態の個数を k = 1, 2, 3, . . . のように変化させながら，コスト関数

CostM (θ) + CostC (X[S]|θ) が最小となるような k を求める．
R EGIME G ENERATION が処理されるとき，S EGMENTA SSIGN MENT において変化点の候補が検出されていない，もしくは γ

保証をすでに終えている場合，新たに作られたレジームは最適であることが保証されている．なぜならば，既存のレジームの中で最も当てはまりの良いモデルと比較して，それを複数に分割し，新たなレジームを作成した方が，パターンをうまく表現できているという結果が得られるからである．
しかしながら，γ 保証が完了していない状態で R EGIME G EN ERATION により新たなレジームが発見された場合，γ 保証が完了した後に遷移したであろうレジームと，新たに作られたレジームは同じ時系列パターンを表現している可能性がある．
よって，本研究では，時刻 t における Lbest の変化点によるセグメント集合に必要なコストと，R EGIME G ENERATION によって出力されたレジームのカレントウィンドウに必要なコストを比較することで，新たに作られたレジームが解 C 内のレジームと重複しないことを保証する．より具体的には次の 2 つのコストを計算し，CostCU T が CostSW より小さくなるとき，新たなレジームを作成する．

づき，2 つのレジームのモデルパラメータ {θ1 , θ2 } を推定する．

CostSW = CostC (X|, ΘC )(11)

HMM のパラメータの学習には，Baum-Welch アルゴリズムを

CostCU T = CostC (X|θ1 , θ2 , ∆) + CostM (θ)(12)

用いる．分割して得られたモデルパラメータが，分割しない場合のモデルパラメータを用いたときのコスト関数 (式 (5)) をより小さくする場合，R EGIME G ENERATION は新たなパターンを表現するレジームを作成する．
モデルパラメータの初期化． R EGIME G ENERATION では，はじめにモデルパラメータ {θ1 , θ2 } を初期化する必要がある．最も簡易的な方法としては，カレントウィンドウの中に含まれる部
分シーケンスをランダム抽出し，モデルの初期値に設定することである．しかし，この方法を用いる場合，初期値に大きく依存するため局所解へ収束してしまう可能性がある．そこで，本研究ではこの問題を解決するため，サンプリングに基づく手法を提案する．まずカレントウィンドウの中から複数個のセグメ
ント/部分シーケンスをサンプルとして均等に取り出す．次に，それぞれのサンプルセグメント s に対し，モデルパラメータ θs

ここで，ΘC は，Lbest に含まれるすべてのモデルパラメータの集合とする．
最後に，新たな時系列パターンが発見されず，既存のレジー
ムに割り当てられた場合，割り当てられたデータをもとに，パラメータ θ を更新する．
レジーム遷移確率 ∆ の更新． 未知レジームの発見により候補
解 C 中のレジーム数 r が変更された場合，レジーム遷移確率 ∆を以下の式によって更新する．
δii =

ここで

P
P

s∈Si

|s| −

P

s∈Si

Pr

s∈Si

j=1

|s|

Nij
,

δij = P

Nij
s∈Si

|s|

(13)

|s| はレジーム θi に所属するセグメントの長さの総和を示し，Nij は θi から θj へのレジームの切替回数を示す．

beakswings

1

tail feathers
claps

tail feathers
beaks
claps
wings

left/rightlegs

6. 2 計算コスト
本研究では S TREAM S COPE の性能を確認するために実際的な状況で実験を行った．図 4 は S TREAM S COPE と，時系列シー

0.5

ケンスから特徴抽出するためのオフライン手法である Auto-left/rightarms 0

600

800

1000

1200

Plait [11]，pHMM [25] との計算時間の比較である．pHMM は
タ集合には MoCap を用い，シーケンス長 n を変化させている．

3

パラメータを必要とするため，ǫr = 0.1,ǫc = 0.8 とした．デー

S TREAM S COPE の処理はカレントウィンドウ内のデータのみ

4

1

400

2

200

を必要とする．そのため，S TREAM S COPE の計算時間はシーケ
ンス長 n に依存することなく，すべてのデータを必要とするオ(a) S TREAM S COPE の出力結果

beaks!"wings

tail feathers

claps#

フライン処理に比べ高速に動作する．S TREAM S COPE はオフライン手法である AutoPlait と比較して 100 倍以上の性能向上を達成した．

(b) 「チキンダンス (chicken dance)」における 4 つの代表的なステップ

6

10

MoCap データにおける S TREAM S COPE の出力結果．

6. 評 価 実 験
本論文では S TREAM S COPE の有効性を検証するため，実デー
タを用いた実験を行なった．実験は 16GB のメモリ，Intel Corei5 3.3GHz の CPU を搭載した iMac 上で実施した．実験で用いたデータは以下の通りであり，各々は平均値と分散値で正規化

Wall clock time (sec.)

図3(z-normalization) して使用した．
• MoCap: MoCap は，1 秒 120 フレームでヒトの動きを計

StreamScope
AutoPlait
pHMM

4

10

2

y=x
2

10

y=x

0

10

−2

10

−4

10

3

4

10

10

Sequence length

測したモーションキャプチャのデータセットである．（注 5） 本実

図4

験ではデータの中から左右の腕と足の 4 次元から構成される加S TREAM S COPE の計算コスト．

速度の値を使用した．

•

GoogleTrend（注 6）: このデータセットは，Google による検索クエリの頻度を週毎に 12 年間に渡り集計したものである．各
シーケンスは各クエリの出現頻度を表す．
6. 1 時系列データからの特徴抽出
図 3 (a) は，MoCap データにおける「チキンダンス (chicken

7. アプリケーション
本章では，S TREAM S COPE の実用的なアプリケーションとして，GoogleTrend データを使用して，Web 上のユーザ動向の情報を自動検出する例を紹介する．S TREAM S COPE は時系列デーの時系列シーケンスデータと，S TREAM S COPE のタの中から，未知のパターンと任意の数のレジームを自動的に出力結果である．このモーションは，4 次元のシーケンスで構発見することができる．また，さらに解析を進めることで，一成され，それぞれの次元が，左右の腕と足の加速度を表現して度発見したパターンと類似したパターンを発見するこも可能でいる．チキンダンスは，図 3 (b) に示す通り，beaks, wings, tailある．

dance)」（注 7）

feathers, claps の 4 つの代表的なステップ から構成される．図 3

外れ値検出． 図 5 (a) は，スターウォーズに関連するキーワーの下の段は，S TREAM S COPE が自動抽出した 4 つのレジームを

ド 5 つ（“star wars movie”, “star wars video” 等）の 12 年間の検示している．ストリーム処理を進めるうち，S TREAM S COPE は索数を示している．このデータは毎年 12 月に検索数が増加す“beak” ，“wing” と次々にそれぞれのステップの特徴をレジームるという周期性を持っている．しかし，それ以外に，映画の公として抽出し，4 つすべての特徴を捉えることに成功している．

開によって検索数が爆発的に増加するという特徴を持っている．

さらに，S TREAM S COPE は “clap” から再び “beak” にステップ

はじめに，S TREAM S COPE は 2005 年に映画が公開されたことが変化したことを捉えていることから，過去に抽出した特徴をを例外的なパターン (レジーム#2) として検出することに成功表すモデルが正しく学習できていること，そしてそれに類似すしている．その後，通常時の周期性 (レジーム#1) に戻り，2015

るパターンを検出することができたと言える．

年に映画が公開されたとき，過去に映画が公開されたときと同様のパターンとして，レジーム#2 への変化を検出することに成功している．

（注 5）：http://mocap.cs.cmu.edu/

トレンド発見． 図 5 (b) はスマートフォンに関連するキーワー（注 6）：http://www.google.com/insights/search/ド（“iphone”, “galaxy” 等）の時系列データであり，それぞれの（注 7）：Chicken dance: http://www.youtube.com/watch?v=6UV3kRV46Zs&t=49s

(a) スターウォーズ関連ワード （レジーム r = 2）(b) スマートフォン関連ワード （レジーム r = 5）
図5

GoogleTrend データにおけるトレンドの変化点抽出例．

データが新機種の発表や発売，または 12 月に検索数が増加する
複雑なシーケンスになっている．S TREAM S COPE は過去 12 年間のスマートフォンのシェア拡大の様子を 5 つのレジームとしてそれぞれリアルタイムに検出した．具体的には，(1) スマー
トフォンという言葉はまだ普及しておらず，欧米で blackberryがビジネスマンの間で流行し始めた時代．(2)iPhone が発売され，スマートフォンが注目を浴び始めた時代．(3)Android 搭載
端末としての Galaxy シリーズのシェア拡大期．そして，(4)LGエレクトロニクスのスマートフォンが普及し始め，(5)iOS 端末の iPhone，そして Android 端末の Galaxy，LG スマートフォン
が世界シェアのほぼ全体をしめる今日の傾向といった 5 つの特徴を捉えている．

8. む す び
本論文では自動パターン検出のためのストリームアルゴリズ
ムとして S TREAM S COPE を提案した．S TREAM S COPE は，与えられた時系列ストリームに対し，ユーザの直感に合う複雑な時
系列パターン（レジーム）とその変化点を発見することができる．様々な種類の実データを用いて実験を行い，S TREAM S COPE
はオンライン処理でありながら，その有用性を確認することができた．また，計算時間についても，従来の手法と比較し大幅に向上していることを示した．
文

献[1] G. E. Box, G. M. Jenkins, and G. C. Reinsel. Time Series Analysis:
Forecasting and Control. Prentice Hall, Englewood Cliffs, NJ, 3rdedition, 1994.

[2] D. Chakrabarti, S. Papadimitriou, D. S. Modha, and C. Faloutsos.
Fully automatic cross-associations. In KDD, pages 79–88, 2004.
[3] L. Chen and R. T. Ng. On the marriage of lp-norms and edit distance.
In VLDB, pages 792–803, 2004.
[4] S. Fine, Y. Singer, and N. Tishby. The hierarchical hidden markov
model: Analysis and applications. Machine Learning, 32(1):41–62,
1998.
[5] E. B. Fox, E. B. Sudderth, M. I. Jordan, and A. S. Willsky. Sharing features among dynamical systems with beta processes. In NIPS,
pages 549–557, 2009.
[6] Y. Fujiwara, Y. Sakurai, and M. Yamamuro. Spiral: efficient and exact model identification for hidden markov models. In KDD, pages
247–255, 2008.
[7] E. J. Keogh, S. Chu, D. Hart, and M. J. Pazzani. An online algorithm
for segmenting time series. In ICDM, pages 289–296, 2001.
[8] J.-G. Lee, J. Han, and K.-Y. Whang. Trajectory clustering: apartition-and-group framework. In SIGMOD Conference, pages 593–
604, 2007.
[9] L. Li, J. McCann, N. Pollard, and C. Faloutsos. Dynammo: Mining
and summarization of coevolving sequences with missing values. In
KDD, 2009.
[10] L. Li, B. A. Prakash, and C. Faloutsos. Parsimonious linear fingerprinting for time series. PVLDB, 3(1):385–396, 2010.
[11] Y. Matsubara, Y. Sakurai, and C. Faloutsos. Autoplait: Automatic
mining of co-evolving time sequences. In SIGMOD, 2014 (to appear, http://www.cs.kumamoto-u.ac.jp/∼yasuko/PUBLICATIONS/
sigmod14-autoplait.pdf).
[12] Y. Matsubara, Y. Sakurai, C. Faloutsos, T. Iwata, and M. Yoshikawa.
Fast mining and forecasting of complex time-stamped events. In
KDD, pages 271–279, 2012.
[13] Y. Matsubara, Y. Sakurai, B. A. Prakash, L. Li, and C. Faloutsos. Riseand fall patterns of information diffusion: model and implications. In
KDD, pages 6–14, 2012.
[14] A. Mueen and E. J. Keogh. Online discovery and maintenance oftime series motifs. In KDD, pages 1089–1098, 2010.
[15] R. T. Ng and J. Han. Clarans: A method for clustering objects for spatial data mining. IEEE Trans. Knowl. Data Eng., 14(5):1003–1016,
2002.
[16] S. Papadimitriou, J. Sun, and C. Faloutsos. Streaming pattern discovery in multiple time-series. In Proceedings of VLDB, pages 697–708,
Trondheim, Norway, August-September 2005.
[17] T. Rakthanmanon, B. J. L. Campana, A. Mueen, G. E. A. P. A.
Batista, M. B. Westover, Q. Zhu, J. Zakaria, and E. J. Keogh. Searching and mining trillions of time series subsequences under dynamic
time warping. In KDD, pages 262–270, 2012.
[18] J. Rissanen. A Universal Prior for Integers and Estimation by Minimum Description Length. Ann. of Statist., 11(2):416–431, 1983.
[19] Y. Sakurai, C. Faloutsos, and M. Yamamuro. Stream monitoring under the time warping distance. In ICDE, pages 1046–1055, 2007.
[20] Y. Sakurai, S. Papadimitriou, and C. Faloutsos. Braid: Stream mining
through group lag correlations. In SIGMOD, pages 599–610, 2005.
[21] Y. Sakurai, M. Yoshikawa, and C. Faloutsos. Ftw: Fast similarity
search under the time warping distance. In PODS, pages 326–337,
2005.
[22] Y. Tao, C. Faloutsos, D. Papadias, and B. Liu. Prediction and indexing of moving objects with unknown motion patterns. In Proceedingsof ACM SIGMOD, pages 611–622, 2004.
[23] N. Tatti and J. Vreeken. The long and the short of it: summarising
event sequences with serial episodes. In KDD, pages 462–470, 2012.
[24] M. Toyoda, Y. Sakurai, and Y. Ishikawa. Pattern discovery in data
streams under the time warping distance. VLDB J., 22(3):295–318,
2013.
[25] P. Wang, H. Wang, and W. Wang. Finding semantics in time series.
In SIGMOD Conference, pages 385–396, 2011.
[26] T. Zhang, R. Ramakrishnan, and M. Livny. Birch: an efficient dataclustering method for very large databases. In SIGMOD, pages 103–
114. ACM, 1996.