DEIM Forum 2014 D4-2

大規模時系列データからの特徴自動抽出松原 靖子†

櫻井保志†

ChristosFaloutsos††† 熊本大学†† Carnegie Mellon University
E-mail: †{yasuko,yasushi}@cs.kumamoto-u.ac.jp, ††christos@cs.cmu.edu
あらまし 本論文では，大規模時系列データのための特徴自動抽出手法である AUTO P LAIT について述べる．AUTO -P LAIT は，様々な時系列パターンを含む複雑なシーケンスが与えられたときに，そのシーケンスデータの中から重要な特徴を発見し，それらの情報を統計的に要約，表現する．実データを用いた実験では，AUTO P LAIT が様々な時系列
データの中から有用なパターンを正確に発見することを確認し，さらに，最新の既存手法と比較を行い提案手法が大幅な精度，性能向上を達成していることを明らかにした．
時系列データ，特徴自動抽出

らの大規模な時系列シーケンスの中から，典型的なパターンや異常値を発見することは非常に重要な課題である．本論文では，大規模時系列データを対象とし，重要な時系列パターンの200抽出を自動的に行なうことを目的とする．そして，大規模時系

600

800

1000

1200

1400

(a) AUTO P LAIT の出力結果

列データの中から，異なるトレンドを発見し，すべての時系列
パターンを表現する手法として，AUTO P LAIT を提案する．（注 1）4001

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

0.5

left/right
0
arms
2

時系列シーケンスは，センサデータや Web アクセス履歴等，

tail feathers
tail feathers
beaks
left/right beaks
claps
claps
wings
legs 1wings

3

1. ま え が き

4

キーワード

beaks

!"wings

tail feathers

claps#

より具体的には，(a) X の中のパターンの変化点を発見し，部
分シーケンス集合（セグメント）に分割し，(b) それらのセグ
メントをグループ化し，類似時系列パターン（本論文では「レジーム (regime)」と呼ぶ） を発見する．さらに重要な点として，これらの処理は (c) 高速かつ，自動で行う．
具体例． 図 1 (a) は，MoCap データにおける「チキンダンス(chicken dance)」の時系列シーケンスデータと，AUTO P LAIT

（注 2）の出力結果例である．このモーションは，4 次元のシーケンスで構成され，それぞれの次元が，左右の腕と足の加速度を表現している．チキンダンスは，図 1 (b) に示す通り，beaks, wings,(b) 「チキンダンス (chicken dance)」における 4 つの代表的なステップ
図1

MoCap データにおける AUTO P LAIT の出力例．

索 [2], [9], [13] 等は重要な課題である．しかし，これらの先行研
究は，基本的にすべて，パラメータの設定やチューニングが必要である．例えば，文献 [3], [15] は，セグメントの個数や，エラーの閾値等，いくつかのユーザ指定のパラメータ入力が必要
であり，これらのパラメータが出力結果に大きな影響を与える

tail feathers, claps の 4 つの代表的なステップ から構成される． （5. 章を参照）．したがって，理想的にはこれらのパラメータ設
定やユーザの介入を必要としない手法が望ましい．
図 1 (a) の下の段は，AUTO P LAIT が自動抽出した 4 つのレジー
ムを示している．提案手法は，ダンスに含まれる 4 つのステップを抽出し，そして各ステップの切れ目も正しく発見することができる．ここで最も重要なこととして，AUTO P LAIT は，これらの 4 つのステップに関する事前知識を必要とせず，適切な数のレジームとその位置を自動的に把握することができる．

1. 1 関連研究と自動抽出手法の重要性
時系列データを対象とした研究課題は，数多く存在する．パ
ターン発見 [12], [8], [11], [14]，情報要約 [1], [6], [7]，クラスタ
リング [4]，セグメンテーション [3], [15] や類似シーケンス探（注 1）：ソースコード： http://www.cs.kumamoto-u.ac.jp/∼yasuko/software.html（注 2）：Chicken dance: http://www.youtube.com/watch?v=6UV3kRV46Zs&t=49s

また，さらに重要な問題がビックデータの解析である．時系

列データは様々なシステムにおいて大量に発生している．大規
模なデータを解析するにあたり，ユーザの手を介したパラメー
タ設定を行なうことは，多くの時間的コストが必要となるため，現実的ではない．すなわち，ビックデータの解析において，自動処理技術は必要不可欠な重要な要素である．

1. 2 本論文の貢献
AUTO P LAIT は以下の特長がある．
（ 1 ） 時系列パターン（レジーム） の個数と種類を把握し，それぞれの適切な変化点を発見する．さらに AUTO P LAIT は，
パターン変化点の最適解の検出を保証する．
（ 2 ） 3. 章で述べる提案モデル (MLCM) により，ユーザの

直感に合致した時系列パターンの抽出を行なう．

列（定義 4）とコスト関数（式 (5)）についての詳細は，3. 章に（ 3 ） AUTO P LAIT はパラメータ設定を必要としない．ユーおいて示す．

ザの介入を必要とせず，適切なレジームの数，変化点の数を，問 題 1 で 示 し た 通 り，本 論 文 の 目 的 は ，X の 特 徴 を 抽

自動的に発見することができる．

出し，すべての時系列パターンを表現するパラメータ集合（ 4 ） 計算コストは入力データの長さに対して線形である．

{m, r, S, Θ, F } を発見することである．ここで，この全パラ

2. 問 題 設 定

メータ集合を候補解 C と呼ぶ．

こ こ で は 本 論 文 で 必 要 な 概 念 に つ い て 定 義 を 行 な う．

X = {x1 , . . . , xn } を d 次元の時系列シーケンスとし，xt を時
刻 t における d 次元ベクトルとする．シーケンス X が与えられたとき，本研究は X を m 個のセグメント集合 S= {s1 , . . . , sm }に分割することを目的とする．si は i 番目のセグメントの開始点，終了点で構成され（つまり，si = {ts , te }），各セグメントは重複がないものとする．本研究ではさらに，発見したセグメ
ント集合を類似セグメントのグループ（レジーム: regime ）に分類する．
［定義 1］（レジーム） r を最適なセグメントグループの個数とする．それぞれのセグメント s はセグメントグループの１つに割り当てられる．これらグループをレジーム (regime) と呼び，それぞれのレジームは統計モデル θi (i = 1, . . . , r) として表現［定義 3］ X を表現する全パラメータ集合 C = {m, r, S, Θ, F }を候補解と呼ぶ．候補解 C は，セグメント集合，各セグメントのレジームへの割当て，レジームを表現する確率モデル，これらすべてを表現する．
結論として，本論文の目的は最適な解 C を発見することである．
ここで非常に重要な課題は，(a) どのようにセグメントおよびレジームの数を推定するか，(b) どのようにレジームを表現し，セグメントの割当てを行なうかである．本研究では，ユーザによるパラメータ設定を介せず，自動処理によって最適解を求めるための新手法を提案する．

3. データ圧縮と情報要約
本章では，問題 1 を解決するためのモデルを提案する．

される．

提案モデルは以下の 2 つのアイデアに基づく．
（ 1 ） 多階層連鎖モデル (MLCM: multi-level chain model):例えば，図 1 において，モーションは m = 8 個のセグメント

複数のレジーム間の時系列パターンとその遷移を表現するためから構成され，各々のセグメントは r = 4 個のレジーム (beaks,に，多層的な連鎖モデル (MLCM) を提案する．図 2 は提案モ

winds, tail feathers, claps) のうちの１つに割り当てられる．

デルの概念図である．青いセル (state 1,2,3) はレジーム 1 に所

［定義 2］（セグメントメンバーシップ） F = {f1 , . . . , fm } を，

属する隠れマルコフモデルの状態遷移を表現し，赤いセルはレ

m 個の整数列とし，fi を i 番目のセグメントが所属するレジー
ムの番号とする (1 <= fi <= r)．

列 A を持つ．例えばレジームはそれぞれ図 1 のダンスにおけ

ジーム 2 に該当する．各レジーム内部においてそれぞれ遷移行

図 1 では，1 番目のセグメントは 4 番目のレジーム (beaks) に，る beaks と wings のステップに相当する．ここで，提案モデル

2 番目のセグメントは 3 番目のレジーム (wings) にそれぞれの重要な要素として，それぞれのレジーム間のレジーム遷移確

所属する．つまり，この場合のセグメントメンバーシップは率 (∆r×r ) の概念を導入する．レジーム遷移確率は r × r の行

F ={4, 3, 2, 1, 4, 3, 2, 1} となる．

列として表現され，r 個のレジーム間の遷移を表現する（図 2

本研究の目的は，大規模時系列シーケンスが与えられたときの例では，r = 2 となる）．

に，そのシーケンスのセグメント化と分割位置の検出及び，レ（ 2 ） モデル表現コスト: セグメントとレジームの発見のた

ジームの発見を高速かつ自動で行なうことである．本論文で取めに，最小記述長 (MDL: minimum description length) の概念をり組む問題を以下のように定義する．

用いる．MDL は情報理論に基づくモデル選択基準のひとつで［ 問題 1］ 多次元時系列シーケンス X が与えられたとき，Xあり，可逆圧縮を行なうことができるが，そのものの概念だけを表現するような以下の情報を抽出する．
（ 1 ） セグメントの総数 m と各セグメントの位置：では本論文の目的を直接解決することはできない．そこで，与

S = {s1 , . . . , sm }えられたシーケンス X を適切に表現するモデルを見つけるた
めに，新しい符号体系を定義する．具体的には，(a) 新たな関数（ 2 ） レジームの総数 r とセグメントメンバーシップ：（式 (5)）を用いて候補解 C のモデルコストを推定し，(b) 最適

F = {f1 , . . . , fm }（ 3 ） r 個のレジームを表現するモデルのパラメータ集合：解を発見するための効果的なアルゴリズムを提案する．
3. 1 MLCM: 多階層連鎖モデル

Θ = {θ1 , . . . , θr , ∆r×r }これらの情報はコスト関数（式 (5)）を最小化するものを選ぶ．

にあるように，隠れマルコフの状態遷移をレジームにグルー

本論文では，レジームを表現するモデルパラメータ集合 Θを，r 個の隠れマルコフモデル (HMM: hidden Markov model)，{θ1 , . . . , θr }，として表現する．（注 3）ここでさらに，新たな概念として，レジーム遷移行列 ∆r×r を導入する．レジーム遷移行（注 3）：本論文で提案する枠組みは，HMM 以外の時系列モデルに適用することも可能である．

多階層連鎖モデル (MLCM: multi-level chain model) は，図 2
プ化し，階層的な時系列パターンの遷移を表現する．本論文では，これ以降，主に 2 層の遷移について焦点を当てるが，2 層以上の多層遷移を表現することも可能である．ここで，図 2 を用いて MLCM の説明を行なう．図は，総計 5 つの状態 (state)から構成される連鎖モデルであるが，ここでは従来の HMM のように，5 × 5 の遷移行列を用いるのではなく，上位層の状態

単一のレジームのモデル θ は，状態数 k（log∗ (k)）と確率モデ

Θ

δ11

a1;11

θ1

a1;12

state 2

a1;31 state 1state 3

a1:22
a1;23

a1;33

δ12

a2;11

state 1

a2;12

θ2
a2;22

a2;21

δ21

Regime 1 (e.g., “beaks”)
図2

ル (θ = {π, A, B}) の表現コストが必要となる．まとめると，

δ22

Regimeswitchstate 2

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

(2)

ここで，cF は浮動小数点のコストを示す．（注 6） 同様にして，レジーム遷移行列には，CostM (∆) = cF · r2 のコストを要する．
3. 2. 2 時系列シーケンスの符号化コスト
先述の通り，本論文では MLCM モデルを用いてシーケンス

Regime 2 (e.g., “wings”)

X の時系列パターンを表現するが，ここで重要なのは，推定したモデルが X を正しく表現しているかを判断する指標の導入

多階層連鎖モデル Θ （ここでは r = 2）．

である．ハフマン符号を用いた情報圧縮では，モデル θ が与え(super-state) の概念を導入することによって，パターンのグルーられた際の X の符号化コストを負の対数尤度を用いて次のよ

プ化を行なう．ここで，このグループを「レジーム」と呼ぶ．

うに表現することができる．

例えば，図 2 では，3 つの青い状態を用いて beaks のステップを表現し，残りの 2 つの赤い状態を用いて wings のステップを表現する．ここで，各レジームは内部に遷移行列 (a1;ji ∈ A1 ,a2;ji ∈ A2 ) を保持し，それらは黒い矢印で表される．さらに，上位層ではレジーム間の遷移行列 (δvu ∈ ∆2×2 ) を保有し，これらは緑の矢印で示されている．
［定義 4］（レジーム遷移行列） ∆r×r を r 個のレジーム群の遷

CostC (X|θ) = log2

行 列 ∆ 内 部 の 要 素 δi,j は 確 率 を 表 し ，0 <= δij <
= 1,
P
δ
=
1
という条件を持つ．そして提案モデルは
r
個の
j ij
レジーム集合 Θ = {θ1 , . . . , θr , ∆r×r } で表現され，θi は i 番
目のレジームのモデルパラメータを表現する．ここで，θi はHMM に基づき，初期確率，遷移確率，出力確率の三つ組で次のように表現される：θi = {πi , Ai , Bi }. （注 4）

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

CostC (X|Θ) =次に，大規模時系列データを表現するための符号化スキームを導入する．直感的には，データが与えられたときのモデルのよさは次の式で表現できる：CostT = Cost(M) + Cost(X|M)．
ここで，Cost(M) はモデル M を表現するためのコストを示
し，Cost(X|M) は，M が与えられたときのデータ X の符号化のコストを示す．

log (d) ビット

（注 5）

•

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

バーシップ）：m log(r) ビットmX− 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 } が与えられたときの X の符号長は次のように表現される．

CostT (X; C) = CostT (X; m, r, S, Θ, F )= log∗ (n) + log∗ (d) + log∗ (m) + log∗ (r) + m log(r)m−1Xlog∗ |si | + CostM (Θ) + CostC (X|Θ)(5)

各セグメントの長さ s：•

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

Pm−1
i=1

log∗ |si | ビット

CostM (θi ) + CostM (∆).

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

4. 最適化アルゴリズム
前章では，候補解 C = {m, r, S, Θ, F } が与えられた上でシー
ケンス X を表現するためのコスト関数として，式 (5) を示し

•

(1)

i=1

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

=i=1

提案モデルの表現コストは以下の要素から構成される．
• 多次元シーケンスデータの長さ n と次元数 d：log∗ (n) +CostM (Θ) =CostC (X[si ]|Θ)i=1+

3. 2. 1 モデル表現コスト

r
X

m
X

3. 2. 3 符号化コスト関数

3. 2 特徴抽出とデータ圧縮∗(3)

ここで，P (X|θ) は X の尤度を示す．シーケンス X と r 個の移行列と呼ぶ．ここで，要素 δij ∈ ∆ は i 番目のレジームから

j 番目のレジームへの遷移確率を示す．

1= − ln P (X|θ).
P (X|θ)た．続いて本章では，式 (5) に基づき，最適な解 C を発見するためのアルゴリズムを提案する．

4. 1 概要本研究では，前章で述べたコストモデルに基づき，セグメ
ントおよびレジームの個数を自動的に選択する．直感的には，データの圧縮率が高ければ，そのモデルはデータに含まれるパ（注 6）：本論文では 4 × 8 ビットとする．

記号
シーケンス
n
d
X
セグメント
m
S
F
レジーム
r
Θ

θi
ki
∆r×r

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

表1

定義

主な記号と定義．

f1 = 2

θ1
θ2

X に含まれるセグメントの総数
X に含まれるセグメント集合：S = {s1 , . . . , sm }セグメントメンバーシップ：F = {f1 , . . . , fm }(a) Iteration 1 (レジーム: r = 2, セグメント: m = 4)

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

f1 = 2

f1 = 2 f2 = 4 f3 = 3 f4 = 1 f5 = 2 f6 = 4 f7 = 3 f8 = 1

X

ンバーシップ F から構成される候補解 C に対し，X の符号化

θ1
θ3
θ4
θ2

コスト CostT (X; m, r, S, Θ, F ) が最小となるとき，C は最適なモデルになる．
次に，具体的な最適化手法を示す．ここでは問題を簡略化す(c) Iteration 4 (レジーム: r = 4, セグメント: m = 8)

るため，次に挙げる 3 つの部分問題に分割する．(1) レジームの個数を r = 2 に固定し，各レジームのモデルパラメータも与

9000
CostT

えられる場合を考える．(2) r = 2 を固定した状態で，モデルパ
ラメータの推定を行う．(3) 最適なレジームの個数を推定する．
より具体的には各部分問題に対応し，以下の 3 つのアルゴリズ（ 3 ） AutoPlait: 最適なレジームの個数 (r = 2, 3, 4, . . .) を求める．
図 3 は，AUTO P LAIT の処理のながれである．AUTO P LAIT は

r = 2 から開始し，セグメントとレジームを分割しながらシー
ケンス X を適切に表現する解 C を発見する．レジームとセグ
メントを分割するとき，コスト関数（図 3 (d)）が下がる．例外
として，図中における赤丸の箇所 (Iteration 3) については，レジーム θ1 が分割しないため，コストが下がらない．

4. 2 Cut-point search
まず最も単純な部分問題として，シーケンス X と，2 つのレジームのモデルパラメータ {θ1 , θ2 , ∆} が与えられている場合を考える．CutPointSearch は レジームのモデルパラメータに基づき，X のパターンの変化点（つまりセグメントの分割位置）を検出することができる．ここで重要な点として，提案ア
ルゴリズムは探索漏れがないことを保証しながらも，高速かつ単一の走査によって，最適なレジーム変化点の個数と位置を検出することができる．図 4 はレジームの変化点が 1 つの場合における CutPointSearch の処理の様子を示している．これは青い(b) Splitting (r=3)(c) Splitting (r=4)

8000
7000 (a) Splitting (r=2)

ムを提案する．
（ 1 ） CutPointSearch: レジームの個数 (r = 2) とモデルパ推定する．

f6 = 1(b) Iteration 2 (レジーム: r = 3, セグメント: m = 6)

数 m，レジームの個数 r，モデルパラメータ集合 Θ，そしてメに，2 つのレジームを表現するモデルパラメータ (θ1 , θ2 , ∆) を

f3 = 1 f4 = 2 f5 = 3

θ1
θ3
θ2

ターンをよく表現しているといえる．つまり，セグメントの個（ 2 ） RegimeSplit: レジームの個数 r = 2 が与えられたとき

f2 = 3

X

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

それぞれのセグメントの分割位置を検出する．

f4 = 1

f3 = 2

X

時系列の長さ時系列の次元数d 次元の時系列シーケンス

ラメータが与えられたときに，X を 2 つのレジームに分割し，

f2 = 1

0

図3

1

2

3
4
Iteration

567(d) 総コストの変化．
AUTO P LAIT の概要図：AUTO P LAIT は X が与えられたとき，反復処理により適切なセグメント/レジームの個数を求める．

レジーム θ1 から赤いレジーム θ2 へ切り替わる例である．
より具体的には，ここでは X が与えられたとき，コスト関数（式 (4)）を最小とするような，複数のレジーム変化点を発見し，それらを 2 つのセグメント集合 (S1 , S2 ) に分割することを考える．ここで，2 つのセグメント集合は，偶数番目のセグ
メントは 1 つ目のレジームに，奇数番目は 2 つ目のレジームに所属する．例えば，図 3(a) では S1 = {s2 , s4 } は θ1 に，S2 ={s1 , s3 } は θ2 にそれぞれ所属する．
ここで，複数の変化点を検出するために，多階層連鎖モデル(MLCM) の概念を用いる．図 4 のように 2 つのレジーム θ1 , θ2が与えられているとする．ここで，この 2 つのレジームは相互に遷移確率を持つ．モデルが与えられた上での符号化コスト

CostC (X|Θ) = − ln P (X|Θ) を計算するために，本論文では動的計画法に基づくアルゴリズムを提案する．

4. 2. 1 アルゴリズム
シーケンス X と 2 つの レジーム θ1 = {π1 , A1 , B1 }, θ2 =
{π2 , A2 , B2 }, およびレジーム遷移行列 ∆ ={δ11 , δ12 , δ21 , δ22 }が与えられたとき，X の尤度 P (X|Θ) は次のように計算さ
れる．

...

state 1
L1;1 (1) = φ

θ1

L1;1 (3) = φ...

state 2...

state 3
L1;3 (2) = φ

δ12

δ12

...

state 1

θ2

L2;1 (3) = {3} L (4) = {3}2;1

L2;1 (5) = {4}...

state 2

L2;2 (4) = {4} L2;2 (5) = {3} L2;2 (6) = {4}

X

t =1

t=2

t =3

t=4

t=5

t=6

Algorithm 1 CutPointSearch (X, θ1 , θ2 , ∆)1: Input: Sequence X, model parameters of two regimes {θ1 , θ2 , ∆}2: Output: (a) Number of segments assigned to each regime, m1 , m2(b) Segment sets of two regimes S1 , S2

3:

4: /* Compute p1;i (t) and p2;u (t) */5: for t = 1 : n do
6:Compute p1;i (t) for state i = 1, . . . , k1 ; /* Equations 7 and 9 */7:Compute p2;u (t) for state u = 1, . . . , k2 ; /* Equations 8 and 9 */8:Update L1;i (t) for state i = 1, . . . , k1 ;9:Update L2;u (t) for state u = 1, . . . , k2 ;10: end for
11: /* Divide into two sets of segments S1 , S2 */12: Choose the best cut-point set Lbest ;...

図 4 CutPointSearch の様子．X と θ1 , θ2 が与えられたとき，単一の走査でレジームの変化点 (t = 4) を検出する．

13: ts = 1; /* Starting position of first segment */14: for each cut point li in Lbest do
15:

Create segment si = {ts , li };16:

if i is odd thenAdd si into S1 ; m1 = m1 + 1;

17:

else

18:
P (X|Θ) = max

Add si into S2 ; m2 = m2 + 1;19:{p1;i (n)}
max1<
= k1
= i<

// regime θ1

max1<{p2;u (n)} // regime θ2
=u<
= k2

(6)

20:

end if21:

ts = l i ;22: end for

p1;i (t) = maxδ · maxv {p2;v (t − 1)} · π1;i · b1;i (xt ) 21// regime switch from θ2 to θ1 δ11 · maxj {p1;j (t − 1) · a1;ji } · b1;i (xt )23: return {m1 , m2 , S1 , S2 };(7)

して L2;2 (4) = {4} を保持する. 時刻 t = 6 において，アル// staying at regime θ1p2;u (t) = maxδ · maxj {p1;j (t − 1)} · π2;u · b2;u (xt ) 12// regime switch from θ1 to θ2 δ22 · maxv {p2;v (t − 1) · a2;vu } · b2;u (xt )ゴリズムは p2;2 (6) が確率の最大値を持つことを明らかにし，(8)

ここで，p1;i (t) は，時刻 t におけるレジーム θ1 の状態 i の確率の最大値であり，p2;u (t) は，時刻 t における θ2 の状態 u の確率の最大値を示す．時刻 t = 1 においては，各レジームの確率は次のように計算する．

p2;u (1) = δ22 · π2;u · b2;u (x1 )4. 2. 2 理論的な分析［補助定理 1］ CutPointSearch は O(ndk 2 ) の計算量を要する.

// staying at regime θ2

p1;i (1) = δ11 · π1;i · b1;i (x1 )L2;2 (6) = {4} を最適解として出力する．

(9)式 (7) の上段は，レジーム間の切り替え（θ2 から θ1 ）の確率，下段はレジーム θ1 の内部の状態遷移の確率を示す．π1,i , b1;i (xt ),a1;ji はそれぞれ θ1 内における状態 i の初期確率，出力確率，
状態 j から状態 i への遷移確率を示す．
アルゴリズム 1 は CutPointSearch の具体的な処理を示す．ここで L1;i (t) は時刻 t における θ1 の状態 i の変化点集合であり，

L2;u (t) は θ2 の状態 u の変化点集合である．これらは式 (7), (8)の尤度計算に基づき更新される．変化点は，レジームが他方に切り替わった場合，その時刻 t を変化点の候補として集合に加える（詳細は文献 [5] 参照）．CutPointSearch は時刻 t = n において L1;i (n) と L2;u (n) 内のすべての状態 i,u の中から，尤度

P (X|Θ) を最大化するような最適解 Lbest を選出する．
［ 例 1］ ここでは図 4 を用いてアルゴリズムの具体的な例を示
す．時刻 t = 1 と t = 2 において，θ1 内の確率 p1;1 (1), p1;3 (2)はそれぞれ確率の最大値を持つ．時刻 t = 3 において Cut-

PointSearch は p1;3 (2) から p2;1 (3) への変化点の候補を発見す
る（図では θ1 から θ2 への緑の矢印）．同時に，L2;1 (3) = {3}を候補集合として更新する．同様にして，2 番目の候補点と［補助定理 2］ 与えられたモデル Θ = {θ1 , . . . , θr , ∆} に対し，

CutPointSearch は最適な変化点集合を検出することを保証する．
本論文ではスペースの都合で証明は省略する（文献 [5] 参照）．

4. 3 RegimeSplit
次にモデルパラメータの推定について述べる．具体的には，(a) 2 つのレジームのモデルパラメータを推定し，同時に，(b)すべてのレジーム変化点を検出する．本研究では，式 (5) を用いてシーケンス X の表現コストを最小にするようなモデルパ
ラメータの推定を行なう．アルゴリズム 2 は RegimeSplit の処
理を示す．提案アルゴリズムは以下に示す 2 つのステップから構成される反復処理によって，モデルパラメータを求める．
• ステップ 1: CutPointSearch （アルゴリズム 1）を利用し，符号化コストが最小となるレジーム変化点を検出し，セグ
メント集合を 2 つのグループ {S1 , S2 } に分割する．

•

ステップ 2: ステップ 1 で得られたセグメント集合に基づき，2 つのレジームのモデルパラメータ {θ1 , θ2 , ∆} を推定する．ここで，HMM のパラメータの学習には，Baum-Welch ア
ルゴリズムを用いる．
モデルパラメータの推定． HMM のモデルパラメータの推定については，モデル θ に対し，隠れ状態の数 k を与える必要があるが，この k を手動で設定するのは非常に難しい．もし k の値を小さくすれば，データの表現能力が低くなり，適切なセグ
メントおよびレジームを求めることが困難となる．一方で，もし k を大幅に上げてしまうと，オーバーフィッティングを招く．

低ければ，エントリをスタックから取り除き，{θ0 , m0 , S0 } を

Algorithm 2 RegimeSplit (X)

出力する．これらの処理をスタックが空になるまで繰り返す．

1: Input: Sequence X
2: Output: (a) Number of segments assigned to each regime, m1 , m2
3:(b) Segment sets of two regimes, S1 , S2

Algorithm 3 AUTO P LAIT (X)4:(c) Model parameters of two regimes {θ1 , θ2 , ∆}1: Input: Sequence X

5: Initialize models θ1 , θ2 ;2: Output: Complete set of parameters C, i.e.,

6: while improving the cost do

3:(a) Number of segments, m

7:/* Find segments (phase 1) */4:(b) Number of regimes, r

8:{m1 , m2 , S1 , S2 } = CutPointSearch (X, θ1 , θ2 , ∆);5:(c) Segment set S = {s1 , . . . , sm }9:/* Update model parameters (phase 2) */6:(d) Model parameters of regimes Θ = {θ1 , . . . , θr ; ∆}10:

θ1 = BaumWelch (X[S1 ]);7:(e) Segment membership F = {f1 , . . . , fm }11:

θ2 = BaumWelch (X[S2 ]);8: Q = ∅; /* Q: stack for number of segments, segment set, regime */12:

Update regime transitions ∆; /* Equation 10 */9: S = ∅; m = 0;

13: end whiler = 0;

S0 = {1, n}; m0 = 1;

10: θ0 = BaumWelch (X[S0 ]); /* Estimate model θ0 of S0 */11: Push an entry {m0 , S0 , θ0 } into Q;

14: return {m1 , m2 , S1 , S2 , θ1 , θ2 , ∆};12: while stack Q =| ∅ do

そこで本研究では，隠れ状態の個数を k = 1, 2, 3, . . . のように

13:

Pop an entry {θ0 , m0 , S0 } from Q;

14:/* Try to refine a regime */変化させながら，コスト関数 CostM (θ) + CostC (X[S]|θ) が15:{m1 , m2 , S1 , S2 , θ1 , θ2 , ∆} = RegimeSplit (X[S0 ]);最小となるような k を求める．

16:/* Compare single regime θ0 v.s. regime pair θ1 and θ2 */17:

if CostT (X; S0 , θ0 ) > CostT (X; S1 , S2 , θ1 , θ2 ) then

さらに，本研究ではレジーム遷移確率 ∆ についても符号化コストを最小にする必要がある．そこで，セグメントの変化点

18:/* Regime pair win - split regime */Push entries {m1 , S1 , θ1 }, {m2 , S2 , θ2 } into Q;

19:

集合 {S1 , S2 } とモデルパラメータ {θ1 , θ2 } が与えられたとき

20:のレジーム遷移確率 ∆ = {δ11 , δ12 , δ21 , δ22 } をラグランジュ乗

21:/* Single regime win - no more split, leave it out of the stack */22:

S = S ∪ S0 ; Θ = Θ ∪ θ0 ; r = r + 1;

数法に基づき次のように計算する．

δ11 =
ここで

P

P

s∈S1

P

s∈S1

|s| − N12

s∈S1

|s|

,

δ12 = P

N12
s∈S1

|s|

,

(10)

|s| はレジーム θ1 に所属するセグメントの長さ

else

23:

Update ∆r×r ; /* 10 */24:

fi = r(i = m + 1, . . . , m0 ); m = m + m0 ;end if25:

26: end while
27: return C = {m, r, S, Θ, F };の総和を示し，N12 は θ1 から θ2 へのレジームの切替回数を示
す．δ21 , δ22 についても同様に計算できる．
また本研究ではモデルパラメータの適切な初期値を求めるた
めに，サンプリングに基づく手法 [5] を用いている．

4. 4 AutoPlait
本論文の最終目標は，問題 1 で述べた通り，大規模時系列
データ の中から任意の数のパターンを自動的に抽出することである．ここで解決すべき問題は，(a) 最適なセグメントおよびレ
ジームの個数 m, r はどのように決定すればよいか，(b) それぞれのセグメントを適切なレジームに割り当てるにはどうしたらいいか，の 2 点である．

4. 4. 1 アルゴリズム
大規模時系列シーケンスの中から自動的にパターンを取り出す手法として AUTO P LAIT を提案する．AUTO P LAIT はスタックを用いた手法であり．貪欲法に基づくアルゴリズムである．

AUTO P LAIT は与えられたシーケンスをセグメントに分割し，新たなレジームを生成し，コスト関数である式 (5) を減少させていく．アルゴリズム 3 は AUTO P LAIT の処理の流れを示している．各ステップにおいて，AUTO P LAIT はスタック Q の中か4. 4. 2 理論的な分析［補助定理 3］ AUTO P LAIT の計算量はシーケンスの長さ n に対し線形である．

5. 評 価 実 験
本論文では AUTO P LAIT の有効性を検証するため，実データを用いた実験を行なった．実験は 32GB のメモリ，Intel Core 2

Duo 1.86GHz の CPU を搭載した Linux のマシン上で実施した．
実験で用いたデータは以下の通りであり，各々は平均値と分散値で正規化 (z-normalization) して使用した．
• MoCap: MoCap は，1 秒 120 フレームでヒトの動きを計測したモーションキャプチャのデータセットである．（注 7） 本実
験ではデータの中から左右の腕と足の 4 次元から構成される加速度の値を使用した．

•

GoogleTrend（注 8）: このデータセットは，Google による検索クエリの頻度を週毎に 9 年間に渡り集計したものである．各
シーケンスは各クエリの出現頻度を表す．
5. 1 時系列データからの特徴抽出

らエントリ {θ0 , m0 , S0 } を取り出す．続いて，現在のレジーム提案手法である AUTO P LAIT の情報抽出の効果を検証する

θ0 の分割を試み，新たなレジームの候補ペア {θ1 , θ2 } とそのため，大規模時系列データの解析のための最新の手法として

セグメント集合 {S1 , S2 } を生成する．もし新しいレジームの候

DynaMMo [3], pHMM [15] と比較した．図 5 は MoCap データ

補のコストが現在のレジームのコストより低い場合は（つまり，レジームの候補ペアが勝った場合），AUTO P LAIT は候補ペアをスタック Q に追加する．もし現在のレジームのコストのほうが（注 7）：http://mocap.cs.cmu.edu/（注 8）：http://www.google.com/insights/search/

Recall

Ideal
AutoPlait

0.9

pHMM

0.8

DynaMMo

0.6
DynaMMo:

N/A

0.4
0.2

0.7
0.2

0.4

0.6
0.8
Precision(a) 適合率と再現率
図6(a) AUTO P LAIT（パラメータ不要）
lag

lag

lag

m=9

10 9 8 7 6 5 4 3 2 1(b) DynaMMo （パラメータ：m = 9）

図5

0.8
1

1000 2000 3000 4000 5000 6000 7000 8000

7 6 5 4 3 2 1

0.5
left/rightarms
0

kick jump jump wave hands(left) (right)
walk

Entropy (CE)walk run jump walkleft/right1
legs

5. 2 精10

AutoPlait

pHMM(b) CE スコア

AUTO P LAIT の精度．

度

続いて，与えられたシーケンスに対する提案手法の変化点抽出とクラスタリングの精度について検証する．
変化点抽出の精度． 図 6 (a) は，提案手法と比較手法におけるセグメントの分割点の抽出の精度を適合率と再現率に基づき比較したものである．ここでは MoCap データから合計 20 個のシーケンスを選び使用した．各シーケンスにはそれぞれ平均して 10 個程度のセグメントが含まれており，本実験で使用した
データはおよそ n = 20 × 10, 000 のモーションフレームを含ん(c) pHMM （パラメータ：ǫr = 0.1, ǫc = 0.8）
MoCap における AUTO P LAIT の出力結果のようす．

でいる．図 6 (a) では，AUTO P LAIT は 1 点のみで表現されている．これは，提案手法がパラメータを持たず，出力結果が 1 つに定まるためである．図のように，提案手法は，95% 以上の正におけるモーションの特徴抽出の結果を示している．図 1 にお解率を示している．

いても既にモーションの解析結果を示しているが，ここで使用

一方，DynaMMo は，図 5 で示した通り，ユーザの設定するパしたデータは，異なるユーザによる，より複雑な動きから生成

ラメータ m を要する．ここではパラメータを m = 2, 4, 6, . . . , 30

したデータである．このシーケンスは walking, jumping, kickingのように変化させてそれぞれ精度を比較している．同様に，等の複数のパターンを含み，walking 以外のパターンは一度しかpHMM はモデルの学習精度に関連する閾値のパラメータを要出現しない．従来のクラスタリング手法にはない AUTO P LAIT

する．本実験では ǫr を 0.1 から 10.0 に変化させながら精度をの強みの 1 つとして，重複のないトレンドの検出が挙げられる．

検証した．DynaMMo と pHMM はパラメータによって適合率

実際に，図 5 (a) において AUTO P LAIT は running, jumping のよと再現率が大きく左右されることがわかる．

うな，繰り返し出現しないモーションを抽出し，m = 9 個のセクラスタリング精度． AUTO P LAIT は与えられたシーケンスのグメントと r = 7 個のレジームの抽出に成功している．

中からレジームを発見すると同時に，各クラスタがどのレジー

一方，図 5 (b) は DynaMMo によるセグメントの分割結果でムに所属するかの情報（つまり，セグメントメンバーシップ）

ある．DynaMMo はユーザの指定するパラメータとして，セグを抽出することができる．図 6 (b) は提案手法と比較手法におメントの個数 m を与える必要があるため，ここでは m = 9

けるクラスタリングの精度を示している．より具体的には，レを用いて処理を行っている．図 5 (b) では，点線が正解値を示

ジームの正解ラベルおよび推定したラベルに関する混合行列

し，色のついた実線が DynaMMo によって得られた値である．

(CM: confusion matrix) を作成し，条件付きエントロピー (CE:

DynaMMo は正しいセグメントの個数として m = 9 を与えた

conditional entropy) [5] を計算した．図 6 (b) は，AUTO P LAIT と場合であっても，変化点の位置を正確に特定することができな

pHMM における CE の平均スコアを示している．pHMM と異い．さらに重要なこととして，DynaMMo はクラスタリングのなり，提案手法はほぼすべてのレジームを正しく抽出している．

能力がないため，レジームを発見することができない．

なお，DynaMMo はクラスタ発見の能力を持たない．

図 5 (c) は pHMM の出力結果である．pHMM は，モデル学

5. 3 計算コスト

習の際のエラー値に関連する閾値（すなわちパラメータ）とし

図 7 はシーケンスの長さ n を変化させた際の AUTO P LAIT とて ǫr と ǫc の 2 つを設定しなくてはいけない．文献 [15] を参考

比較手法における計算コストを示している．ここでは MoCapにして，数多くの閾値を試し，それらの中から最も適したもの

データを用いた．DynaMMo はパラメータとしてモデルの隠を選び，図 5 (c) に示している．図に示すように，pHMM は線

れ状態の個数を必要とするため，本実験では k = 4 とした．

形のパターンをモデル化する手法であるため，複雑な時系列パ

pHMM については ǫr = 0.1, ǫc = 0.8 とした．AUTO P LAIT とターンや長期的なトレンドを抽出することができない．

DynaMMo はデータの長さに対し，線形 O(n) である（対数ス
ケールにおいて傾きは slope = 1.0 である）．一方，pHMM

2

Value

4

10

5

"flu fever""flu symptom" April 2009
"flu headache" "Swine flu"flu medicine"outbreak"Oct 2009"Swine flu
vaccine"

3

10

0

2

10

4

5

10

6

10
10
Length of bundle

図7

50

100

150

200 250 300
Time (per week)

350

400

450

2 1

AutoPlait
DynaMMopHMM(a) インフルエンザ関連ワード （レジーム r = 2）

7

10は O(n2 ) の計算量を要する（slope ≈ 2.0）．AUTO P LAIT は

1. xbox,15
playstation
10
5

達成している．

0

本章では，AUTO P LAIT の実用的なアプリケーションとして，

50

100

"ps2, ps3"

"wii"

2. wiiDec
Nov 2006 2009
wii launches

150

200 250 300
Time (per week)"android"

3. mobilegames
Dec 2011

350

400

450

3 2 1

pHMM と比較し，n = 100, 000 において 472 倍の性能向上を

6. アプリケーション"xbox"

20

AUTO P LAIT の計算コスト．
Value

Wall clock Time (s)

10

y=x

y=x

GoogleTrend データを使用して，Web 上のユーザ動向の情報を(b) ゲーム関連ワード （レジーム r = 3）

自動検出する例を紹介する．AUTO P LAIT は時系列データの中

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

から，未知のパターンと任意の数のレジームを自動的に発見することができる．
異常検出． 図 8 (a) は，インフルエンザに関連するキーワード 4つ（”flu fever”, ”flu symptom”等）の 9 年間の検索数を示している．このデータは年単位の周期性を持ち，毎年 10 月から 2 月にかけて検索数が増加し，次第に春，夏にかけて減少していく．
この傾向は毎年同様であるが，2009 年は異なる傾向を持つ．これは，豚インフルエンザが世界的に大流行したことが原因である．AUTO P LAIT はこの例外的なパターンをレジーム#1 として検出することに成功している．
ドレンド発見． 図 8 (b) はゲーム産業に関連するキーワード（“xbox”, “wii” 等）の時系列データである．毎年 12 月のクリス
マスにかけてピークを持つ．近年，ゲーム産業界は様々な企業の参入により，競争的になりつつある．AUTO P LAIT は過去 9年間のゲーム機戦争における 3 つのレジームを発見している．
具体的には，(1) Xbox と Playstation が主流である時代から，(2)Wii が 2006 年に販売開始して以来の大幅なシェアの獲得，そして (3) 高性能携帯端末の普及によるモバイルゲームやソーシャ
ルゲームの出現への流れを捉えている．

7. む す び
本論文では大規模時系列データのための特徴自動抽出手法として AUTO P LAIT を提案した．AUTO P LAIT は，ユーザによるパラメータ設定や事前知識を要すること無く，与えられた大規模シーケンスに対し，ユーザの直感に合う複雑な時系列パターン（レジーム）とその変化点を発見することができる．様々な種類の実データを用いて実験を行い，AUTO P LAIT は最新の時
系列解析手法と比べてより高い精度と性能を持つことを示した．
文

献[1] 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.
[2] Y. Fujiwara, Y. Sakurai, and M. Yamamuro. Spiral: efficient and exact model identification for hidden markov models. In KDD, pages
247–255, 2008.
[3] L. Li, J. McCann, N. Pollard, and C. Faloutsos. Dynammo: Mining
and summarization of coevolving sequences with missing values. In
KDD, 2009.
[4] L. Li, B. A. Prakash, and C. Faloutsos. Parsimonious linear fingerprinting for time series. PVLDB, 3(1):385–396, 2010.
[5] 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).
[6] 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.
[7] 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.
[8] A. Mueen and E. J. Keogh. Online discovery and maintenance oftime series motifs. In KDD, pages 1089–1098, 2010.
[9] 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.
[10] J. Rissanen. A Universal Prior for Integers and Estimation by Minimum Description Length. Ann. of Statist., 11(2):416–431, 1983.
[11] Y. Sakurai, C. Faloutsos, and M. Yamamuro. Stream monitoring under the time warping distance. In ICDE, pages 1046–1055, 2007.
[12] Y. Sakurai, S. Papadimitriou, and C. Faloutsos. Braid: Stream mining
through group lag correlations. In SIGMOD, pages 599–610, 2005.
[13] Y. Sakurai, M. Yoshikawa, and C. Faloutsos. Ftw: Fast similarity
search under the time warping distance. In PODS, pages 326–337,
2005.
[14] M. Toyoda, Y. Sakurai, and Y. Ishikawa. Pattern discovery in data
streams under the time warping distance. VLDB J., 22(3):295–318,
2013.
[15] P. Wang, H. Wang, and W. Wang. Finding semantics in time series.
In SIGMOD Conference, pages 385–396, 2011.