45

第3章

モデリングの基礎

3.1 行動データマイニング(e)

行動データマイニングの概念

図 3.1 に示すように，行動調査の調査法は従来の紙ベースから GPS ベースの時間的・空間的に細かな情報を用いた調査に変わりつつある．それに伴い，計画の単位もマクロからメゾ，ミクロとより細かな単位での評価が可能となっている．ただし，GPS 機能を用いた行動調査の場合，多量の位置デー
タは取得できるものの被験者に付加的に入力を行ってもらわない限り移動手
段やトリップ目的などの情報を得ることはできない．またデータの粒度が細かくなるに従い，観測誤差や欠損といった現象をモデル中で扱うことの重要性も高まっている．
このような課題に対し，データの特徴量（位置，時刻，速度，加速度）を用
いて移動手段を自動的に判別する問題や，誤差や欠損を補正しながらデータを扱うための手法として，データマイニングや機械学習の手法が有用となる．
データマイニングは多量のデータから隠れた法則や有用なパターンを見出したり行動データを正規化して従前のモデルに持ち込むための前処理を行った
りする際に用いられる．本章では，代表的な行動データ分析の対象とその手法について概説する．

(f)

代表的な分析手法

クラスタリング
クラスタリングは，統計・パターン認識・データベース・データマイニン
グなどの分野で広く用いられている古典的な手法である．非常に多くのデー

46

第 3 章 モデリングの基礎

図 3.1

行動データ種別 (羽藤ら (2005))

3.1

行動データマイニング

47

タが集まったとき，これを自動的にいくつかの似ているもの同士グループにわけることができれば，データの概要を効率的に概観しセグメンテーションに基づく細かな情報提供や施策を行える可能性がある．クラスタリング手法
は大きく，最短距離法などの階層的手法と k-means などの分割最適化手法に分けられる．ここでは，階層的手法の代表例である凝集型クラスタリングとk-means 法，さらに確率的なクラスタリングについて述べる．

• 凝集型クラスタリング 今，対象とするデータが N 個あるとする．凝集型クラスタリングでは，1 個の対象だけを含む N 個のクラスタがある初期状態から，クラスタ間の距離（非類似度）関数に基づき最も距離の近い二つのクラスタを逐次的に統合していく手法である．この統合
をすべての対象が一つのクラスタに含まれるまで繰り返すことで階層構造を獲得する．下（バラバラの状態）から上（統合された状態）に順に統合していくのでボトムアップクラスタリングとも呼ばれる．クラスタリングの結果はデンドログラムと呼ばれる樹形図によって表示されることが一般的である．
凝集型クラスタリングをアルゴリズムの形で表すと次のようになる．
ただし，sim(C1 , C2 ) は 2 つのクラスタ間の距離を表す．
ステップ 0
入力：事例集合 D = {x(1) , x(2) , · · · , x(N ) }クラスタ：C = {c1 , c2 , · · · , cN }ステップ 1

1 つのクラスタに 1 つの事例を割り当てる．
c1 = {x(1) }, c1 = {x(2) }, · · · , cN = {x(N ) }ステップ 2
最も似ているクラスタ対を見つける．

(cm , cn ) = arg max sim(ci , cj )ci ,cj ∈Cステップ 3見つかったクラスタ対を融合する．

48

第 3 章 モデリングの基礎

merge(cm , cn )ステップ 4|C| = 1 になるまでステップ 2，3 を繰り返す．
ここで，クラスタ C1 とクラスタ C2 の距離関数 sim(C1 , C2 ) の違いにより以下のような手法がある．

– 最短距離法または単連結法
この手法は，2 つのクラスタが与えられたときその中で最も近い
事例対の類似度をその 2 つのクラスタの類似度とする方法である．

sim(C1 , C2 ) =min

x1 ∈C1 ,x2 ∈C2sim(x1 , x2 )(3.1)

– 最長距離法または完全連結法
この手法は，2 つのクラスタが与えられたときその中で最も遠い事例対の類似度をその 2 つのクラスタの類似度とする方法である．

sim(C1 , C2 ) =max

x1 ∈C1 ,x2 ∈C2sim(x1 , x2 )(3.2)

– 群平均法
群平均法では，与えられた 2 つのクラスタの重心間の類似度を 2つのクラスタの類似度とする方法である．

sim(C1 , C2 ) =

∑
1
n1 n2

∑

sim(x1 , x2 )(3.3)

x1 ∈C1 x2 ∈C2– ウォード法
ウォード法は各対象からその対象を含むクラスタのセントロイドまでの距離の二乗の総和を最小化する手法である．

sim(C1 , C2 ) = E(C1 ∪ C2 ) − E(C1 ) − E(C2 )ただし，

E(Ci ) =∑x∈Ci(sim(x, ci ))2

(3.4)

(3.5)

3.1

行動データマイニング

49

最短距離法，最長距離法，および群平均法は任意の対象完の距離

sim(C1 , C2 ) が与えられている場合に適用できる．対象が属性ベクトルで記述されている場合には属性ベクトル完のユークリッド
距離などを求めて適用する．ウォード法は対象が属性ベクトルで与えられている場合にのみ適用できる．

• k-means 法 k-means 法は，分割の良さに関する評価関数を設定しその評価関数を最適にする分割を探索する手法である．クラスタ数 k は事前に分析者が決定する．各クラスタは平均ベクトルなどの代表ベクトルで表現される．

k-means 法をアルゴリズムの形で表すと次のようになる．
ステップ 0
入力：事例集合 D = {x(1) , x(2) , · · · , x(N ) }クラスタ数：k
無作為に m1 , m2 , · · · , mk を決定．
ステップ 1
すべての x(i) ∈ D に対して以下の操作を行う．

∀c, cm ax = arg max sim(x(i) , mc )cx(i) を cmax に加える．
ステップ 2
各クラスタの代表ベクトルを再計算する．
1∀c, mc = |c|∑(i)x(i) ∈c x

ステップ 3
収束するまでステップ 1,2 を繰り返す．

k-means 法は初期値によって結果が変化するため，適切な初期値を設定することが計算時間短縮，よい解を求めることの双方に重要となる．

50

第 3 章 モデリングの基礎

例えば，凝集型クラスタリングの結果を初期値として利用する方法などが使われることもある．

• 確率的クラスタリング確率モデルに基づいたクラスタリング手法は，
データマイニング，統計，パターン認識の分野で多く研究されている．
これは，i 番目のクラスタについて，パラメータ θi をもつ対象の確率（密度）分布 fi (x|θi ) と，クラスタの混合比 ∀αi > 0,∑ki αi = 1 であらわされる次の混合分布を用いる手法である．

P r[x|θ1 , . . . , θk ] =k∑αi fi (x|θi )

(3.6)

i

fi (x|θi ) には，対象の属性が実数地の場合は正規分布，カテゴリ値の場合は多項分布が用いられることが多い．与えられた対象集合 X が，
この混合分布に基づき発生する場合の尤度を最大にする最尤推定によりパラメータを求めて混合比と確率分布の積が最大になるクラスタへ
各対象を分類する．最尤推定には一般に各クラスタへの確率的割り当てと各クラスタの分布のパラメータ推定を交互に行う EM アルゴリズ
ムが用いられる．
この手法ですべてのクラスタと属性で標準偏差が等しく共分散が単位行列である正規分布を用いた場合と k-means 法では類似した結果が得られ，混合分布を用いた手法は k-means の制約を緩めた手法であると考えることができる．
以下では，fi (x|θi ) に混合正規分布を用いた例について具体的に記述する．
【混合正規分布を用いたクラスタリング】
いま，fi (x|θi ) は正規分布であると仮定する．ただし，分散 σ 2 は既知であり，かつクラスタによって変化しないものとする．よって，簡単のためクラスタ c として fi (x|θi ) = fi (x|c) と書くことにする．mc をク
ラスタ c の平均ベクトルとすると，

3.1

51

行動データマイニング

1|x(i) − mc |2
fi (x(i) |c) = √
exp(−
)
sσ 2
(2πσ 2 )d

(3.7)

ここで，

fi (c|x(i) ) =P (x, c)P (c, x(i) )P (c)P (x(i) |c)
=∑
=∑
(i)
(i)
P (x)
)
|c)
c P (c, x
c P (c)P (x(i)
P (c)exp(− |x 2σ−2mc | )
=
(i)
∑
|x
−mc |
)
c P (c)exp(−
2σ 2

(3.8)

クラスタ c の代表ベクトル mc は確率の重みづけで計算することになるから，以下のようになる．

∑
mc =fi (c|x(i) )x(i)
(i)
x
∈D
∑
f (c|x(i) )x(i) ∈D i(3.9)

以上の計算をアルゴリズムにまとめると以下のようになる．
ステップ 0
入力：ベクトル集合 D = {x(1) , x(2) , · · · , x(N ) }クラスタ数：k
無作為に m1 , m2 , · · · , mk を決定．
ステップ 1
すべての x(i) ∈ D，c に対して以下の操作を行う．
′

(i)
|x
−m c |)
P (c)exp(−(i)
2σ 2
fi (c|x ; m ) =
′
(i)
∑
|x
−mc |
)
c P (c)exp(−
2σ 2′(3.10)

ステップ 2
各クラスタの代表ベクトルを再計算する．

∑′
fi (c|x(i) ; m )x(i)(i)
x∈Dmc = ∑′
f (c|x(i) ; m )x(i) ∈D i(3.11)

52

第 3 章 モデリングの基礎

ステップ 3
収束するまでステップ 1,2 を繰り返す．
今回は簡単のため σ 2 を既知としたが，σ 2 も未知の場合を扱うために′は，代表ベクトルm の再計算に加え，σ 2 も再計算を行う．
相関ルール
相関ルールとは相関ルールは 1993 年に Agrawal らが提起しデータマイニングの発展の一端を担った．元来はスーパーマーケットでの買い物の傾向を分析して商品の配置計画や販売促進に役立てようとするバスケット分析を指向して提起された手法論であるが，その後一般的なデータ解析にも適用可能な柔軟な手法として広く用いられている．相関ルールは次のような指標を用いて評価される．
アイテム集合 X ⊆ L とする．トランザクションとは，（購買データであれば）一人の顧客が購買したアイテムのリスト D を表す．count(X) は，トランザ
クション D の中で X を含むものの個数を表す．

• 支持度 (support)
全トランザクション数に対する，X と Y とを共に含むトランザクショ
ン数の比．P r[X, Y ] に相当する．

count(X)|D|∪
count(X Y )support(X ⇒ Y ) =|D|
support(X) =(3.12)(3.13)• 確信度 (conﬁdence)
X を満たすトランザクション数に対する，X と Y を共に含むトランザ
クション数の比．P r[Y |X] に相当する．

conf idence(X ⇒ Y ) =∪∪count(X Y )support(X Y )=(3.14)
count(X)support(X)【例】「A を買った人の 80%が B, C も買っておりその 3 種類の商品すべてを買った人は全顧客の 6%である」というとき，各指標は以下のようになる．

3.1

行動データマイニング

53[A] ⇒ [B, C]:sup = 6%, conf = 80%
最低支持度（minsup）と最低確信度（minsconf）を指定して，データベース
からすべての相関ルールを抽出するといった方法が用いられることが多い．
相関ルールの問題点
相関ルールの問題点として以下のような点が挙げられている．

• 相関ルールには統計的な相関がないため抽出したルールは無意味である• アイテムが密になると，アイテムの組み合わせ数が膨大となりルールの抽出の計算が困難である．また，大量のルールが抽出された場合，それらの視察が実質的に困難になる．

統計的な相関がないという問題に関しては，条件部を空とした場合との確信度の比であるリフト値という新たな指標が考案された．

list(X ⇒ Y ) =∪conf idence(X Y )support(Y )(3.15)

また，膨大なルールの計算方法についても効率的な計算方法の工夫が提案されている．

(g) 系列パターンの類似性分析
パタン認識，音声認識，画像認識，ゲノム解析など膨大なデータを扱う多くの問題においてデータ文字列をデータベース中の文字列にマッチングさせる問題は非常に重要な問題である．2 つのパターンの類似度は編集距離や相同性スコアと呼ばれる．
行動データにおいても，活動目的や活動パタンを離散的に表す際にこのような概念を用いて類似度を計算することが試みられている．（活動データの分析で具体的な例を示す．
）この方法は，系列データの類似度を計算できるという点で行動データの生成過程のモデリングやシミュレーションにおいても適用可能性は高いと考えられる．
類似度計算におけるマッチングの概念
ここでの類似度の概念は，ある離散的なパターン列 s と t が存在するとき，それらのパターン列類似度が以下の 3 つの操作によって s を t に変換するのに必要な最小コストとして定義される．

54

第 3 章 モデリングの基礎• パターン列中の記号 si が別の文字 tj に置換される．
• パターン列中の記号 si が脱落する．
• パターン列中の記号 tj が挿入される．
ここで，上記の置換，削除，挿入の各操作のコストをそれぞれ cs (i, j), cd (i), ci (j)と置くと，位置 (i, j) における類似度 d(i, j) は以下のような再帰式で記述することができる．

d(i, j) = min{d(i − 1, j − 1) + cs (i, j), d(i, j − 1) + ci (j) + d(i − 1, j) + (3.16)
cd (i)}
このような再帰式を動的計画法を用いて解くことで，最適な操作パターンと類似度を求めることができる．

(h) 因果構造，関係性のモデル (確率的グラフィカルモデル)
確率的グラフィカルモデルとは
本項では，確率分布の図式的な表現である確率的グラフィカルモデルを扱う．確率的グラフィカルモデルでは，確率変数間の関係性をモデル化するため，画像のノイズ除去や修復問題など部分的な観測から全体を復元する問題
や，全体を精度よく表現するための効率的な観測などに関して画像処理や行動データの因果推論など広く研究が進められている．われわれが分析対象とする行動データに関しても，確率変数同士関係を持つ場合の推論問題は多く存在し，これらの手法の適用性は高いと考えられる．
確率的グラフィカルモデルとは以下のような特徴を持っている．

• 確率モデルの構造を視覚化する簡単な方法を提供し，新しいモデルの設計方針を決めるのに役立つ．

• グラフ構造を調べることで，条件付き独立性などモデルの性質に関する知見が得られる．

• 推論や学習を行う際，数学的な表現を暗に伴うグラフ上の操作として表現することができる．
グラフはリンクによって接続されたノードの集まりからなる．確率的グラ
フィカルモデルでは，各ノードが確率変数，リンクがこれらの変数間の確率

3.1

55

行動データマイニング

Ă


ď

図 3.2

グラフィカルモデルの例

的関係を表現し，グラフは「全確率変数上の同時分布が一部の変数のみに依存する因子の積としてどのように分解可能か」という情報を表現する．
代表的なグラフィカルモデルには，グラフのリンクが特定の方向性をもつ
ベイジアンネットワーク（有向グラフィカルモデル）とリンクが方向性を持たないマルコフ確率場（無向グラフィカルモデル）が存在する．有向グラフ
は確率変数間の因果関係を表現でき，無向グラフは確率変数間の緩い束縛関係を表現することができる．また，推論問題を解く際には有効グラフや無向グラフを因子グラフと呼ばれる別の表現に変換することが多い．次項からは，先に述べた代表的な 2 つのモデルと推論問題を解く手法について述べる．
ベイジアンネットワーク
まず，簡単な a, b, c 上の任意の同時分布 p(a, b, c) を考えよう．確率の乗法定理を用いることにより，同時分布はp(a, b, c) = p(c|a, b)p(a, b)= p(c|a, b)p(a|b)p(a)(3.17)となる．このような分解は任意の同時分布に対して可能である．
このような分解の例を図 3.2 に示した．ノード a からノード b へのリンクが存在するとき，a は b の親ノード，b は a の子ノードと呼ぶ．この例で興味深いのは，左辺は 3 変数 a, b, c に関して対象であるのに右辺は対象でないことである．
次に，K 変数の同時分布 p(x1 , . . . , xK ) の場合を考えよう．確率の乗法定
理を繰り返し適用することにより，同時分布は各変数上の条件付き確率分布ϭ

56

第 3 章 モデリングの基礎の積として書ける．

p(x1 , . . . , xK ) = p(xK |x1 , . . . , xK ) · · · p(x2 |x1 )p(x1 )(3.18)

K の値を決めれば，この同時分布は K 個のノードをもつ有効グラフとして表現される．有向グラフと対応する（変数の）分布の間の一般的な関係について以下に述べる．グラフによって定義される同時分布はグラフ上で親に対応する変数によって条件づけられた，各ノード変数上の条件付き分布の（全ノード分の）積によって与えられる．よって，K 個のノードを持つグラフに対応する同時分布は以下の式で与えられる．

p(x) =K∏p(xk |pak )(3.19)

k=1

ただし，pak は xk の親ノード集合を表し，x=x1 , . . . , xK を表す．この等式
は，与えられた有効グラフィカルモデルに対応する同時分布の分解特性を表現している．
多項式曲線フィッティングの例
ここでは，ベイズ多項式回帰モデルについて考える．このモデルの確率変数は，多項式係数ベクトル w および観測データ t=(t1 , . . . , tN )T である．このとき，同時分布は事前分布 p(w) と N 個の条件付き分布 p(tn | w) (n =1, . . . , N ) の積によって与えられる．
p(t, w) = p(w)N∏

p(tn |w)(3.20)

n=1

モデルパラメータを陽に書く場合は次のようになる．
2p(t, w|x, α, σ ) = p(w|α)

N∏p(tn |w, xn , σ 2 )(3.21)n=1

ただし，α はガウス事前分布の精度を表すハイパーパラメータ，σ 2 はノイズの分散である．
マルコフ確率場
有向グラフィカルモデルは，全変数集合上の同時分布を局所的な条件付き

3.1

57

行動データマイニング

分布の積に分解することが可能なモデルであった．無向グラフであるマルコフ確率場においても，因数分解および条件付き独立性が規定される．
ここで，1 つのリンクによって直接接続されない 2 つのノード xi およびxj について考える．これらの変数は，グラフ上の他のすべてのノードが与えられたもとで，条件付き独立でなければならない．この条件付き独立性は以下のように表現される．

p(xi , xj |x\{i,j} ) = p(xi |x\{i,j} )p(xj |x\{i,j} )(3.22)

ここで，x\{i,j} はすべての変数集合 x から xi および xj を取り除いたものである．したがって，この条件付き独立性がグラフに属するすべての可能な分布について成り立つためには，xi と xj とが同じ因子に含まれないように同時分布が因数分解される必要がある．
このようなことから，クリークと呼ばれるグラフ上の概念を導入すると便利である．クリークの定義は，すべてのノードの組にリンクが存在するようなグラフの部分集合である．さらに，極大クリークとは，もう 1 つノードを加えるとクリークでなくなってしまうようなクリークとして定義される．
クリークを C と書き，クリーク内の変数集合をxC と書く．すると同時分
布は，グラフの極大クリーク上のポテンシャル関数の積の形で書ける．

p(x) =1 ∏Z C

ψc (xC )(3.23)

ここで，Z は規格化定数（分配関数）であり，以下で与えられる．

Z=

∑∏
x

ψc (xC )(3.24)

C

ポテンシャル関数は，周辺分布や条件付き確率分布のように確率的解釈が可能なものに限定されない．したがって，一般にそれらの積は正しく規格化されない．よって規格化定数が陽に組み込まれているのである．
有効グラフの場合は，因子分解の各因子である条件付き確率分布が規格化されていたため，同時分布は自動的に規格化されていた．無向グラフでは，M個の K 状態離散変数ノードをもつモデルの場合，規格化因子を計算するた
めには KM 個の状態に対する足し算が必要であり，計算量は最悪の場合モ

58

第 3 章 モデリングの基礎

デルのサイズに対して指数的に増大する．ただし，分配関数はポテンシャル
関数 ψc (xC ) を支配する関数でありパラメータを学習する際には計算が必要となるが，局所的な条件付き分布の計算には不要である．なぜなら，条件付
き分布は 2 つの周辺確率の比であり比を計算する際には分配関数はキャンセ
ルされるからである．
ここで，無向グラフの条件付き独立性と因数分解との関係について考えよう．これを定式化するためには，ポテンシャル関数 ψc (xC ) が狭義に正（任意の xC に対して 0 でも負でもない）場合に制限する必要がある．この制限のもとでは，因数分解と条件付き独立性との間の精密な関係を定義することができる．
ポテンシャル関数は協議に正に限ったので，指数関数を用いて以下のように表す．

ψc (xC ) = exp{−E(xC )}(3.25)

ここで，E(xC ) はエネルギー関数と呼ばれこの指数関数表現はボルツマ
ン分布と呼ばれる．同時分布はポテンシャル関数の積で与えられるので全エ
ネルギーは極大クリークごとのエネルギーの和として与えられる．先にも述べたように，有向グラフの同時分布は確率分布の積に分解されたのに対し無向グラフの因子であるポテンシャル関数は確率的解釈を持たない．規格化制
約がないのでポテンシャル関数を自由に選ぶことができる．
有向グラフモデルと無向グラフモデルの関係
まず，有効フラフで記述されるモデルを無向グラフに変換する問題を考えよう．単純なグラフの場合であればこの変換は容易に実行できる．図の例を考えてみる．有向グラフの同時分布は，条件付き分布の積として以下の形に分解できる．

p(x) = p(x1 )p(x2 |x1 )p(x3 |x2 ) · · · p(xN |xN −1 )(3.26)

これを無向グラフによる表現に変換してみよう．この無向グラフにおける極大クリークは単なる隣接ノード対であるので同時分布は以下のように書ける．

p(x) =

1Z

ψ1,2 (x1 , x2 )ψ2,3 (x2 , x3 ) · · · ψN −1,N (xN −1 , xN )(3.27)

3.1

行動データマイニング

59

したがって，容易に以下のような対応付けができる．

ψ1,2 (x1 , x2 ) = p(x1 )p(x2 |x1 )
ψ2,3 (x2 , x3 ) = p(x3 |x2 )..
.

(3.28)

ψN −1,N (xN −1 , xN ) = p(xN |xN −1 )(3.29)
一般に，有効グラフを無向グラフに変換するには以下のようにすればよい．

1. グラフの各ノードに対してそのすべての親同士の対に無向リンクを付加する．さらにもともとのリンクから矢印の方向性を取り除いてモラルグ
ラフを作る．

2. モラルグラフのすべてのクリークポテンシャル関数を 1 に初期化する．
3. もともとの有向グラフの条件付き分布因子を 1 つ取ってきて，対応する
クリークポテンシャルの 1 つにかける．モラル化を行ったため，各条件付き分布因子の変数をすべて含む極大クリークが少なくとも一つ必ず存在する．
グラフィカルモデルの枠組みは，有効リンクと無向リンクを両方持つグラ
フにも矛盾なく拡張でき，そのようなグラフを連鎖グラフと呼ぶ．
推論問題
ここでは，グラフィカルモデルにおける推論の問題を考える．これは，グラフのいくつかのノードが観測値に固定されたとき，残ったノード（のうちいくつか）に関する事後分布を計算したい，という問題である．
効率の良い推論アルゴリズムを見つけたりそれらのアルゴリズムの構造の見通しを良くする際にグラフの「構造」を利用することができる．具体的には，多くのアルゴリズムが局所的なメッセージのグラフ全体にわたる伝播として表現できる．本節では，木構造グラフにおける厳密推論のためのアルゴ
リズムについて述べる．

60

第 3 章 モデリングの基礎

積和アルゴリズム
積和アルゴリズムとは，ノードあるいはノード部分集合上の局所的な周辺分布を計算する手法である．次節で述べる max-sum アルゴリズムはこの手法の修正したものであり，最も確からしい状態を見つけるのに用いられる．
以下の議論では，モデルの持つすべての変数は離散的であるとする．よって周辺化は和演算に対応するが，ここで得られる枠組みは周辺化が積分演算によって得られる線形ガウスモデルにも同様に適用可能である．ループのない有向グラフにおける厳密推論のためのアルゴリズムとして確率伝播法が存在するが，これは積和アルゴリズムの特別な場合と考えてよい．ここでの目
的は以下の 2 つである．

1. 周辺分布を求めるための効率よい推論アルゴリズムを得る．
2. 複数の周辺分布を計算したい場合に，計算の重複をなくして効率化する．
初めに，ある特定の変数ノード x 上の周辺分布 p(x) を求める問題から考える．今，すべての変数が隠れている（観測されていない）と仮定する．現
実には，一部の変数が観測されたもとでそれらの観測値に条件づけられた事
後分布を計算した場合が多いが，拡張は容易である．
ステップ 1
変数ノード x を因子グラフの根ノードとみなし，以下の式を用いてグ
ラフの葉ノード f におけるメッセージを初期化する．

µx−→f (x) = 1(3.30)µf −→x (x) = f (x)(3.31)

ステップ 2
メッセージパッシングが以下の式で再帰的に適用される．

3.1µfs −→x (x)
∑
∑
=
···
fs (x, x1 , . . . , xM )
x1

=

∑
x1

xM

···

∑

行動データマイニング

∏

[

∑

61

Gm (xm , Xs m)]m∈ne(fs )\x Xs m∏

fs (x, x1 , . . . , xM )xM

µxm −→fs (xm ) (3.32)m∈ne(fs )\x∏µxm −→fs (xm ) =[∑Fl (xm , Xm l)]l∈ne(fs )\x Xm l∏=µfl −→xm (xm )(3.33)l∈ne(fs )\x(3.34)
ステップ 3
根ノードがすべての隣接ノードからのメッセージを受け取ったら，求める周辺分布は以下の式を用いて計算される．

p(x) =∑Fs (x, Xs )(3.35)s∈ne(x)max-sum アルゴリズム
積和アルゴリズムを用いれば，因子グラフによって表現される同時分布 p(x)に関して，その成分変数上の周辺分布を効率よくとくことができる．これ以外に 1) 確率が最大となる変数の組を見つける，2) その時の確率値を求める，という 2 つのタスクが必要になる場合がある．このタスクを実現するための手法が max-sum アルゴリズムであり，積和アルゴリズムと深く関係するア
ルゴリズムによって実現される．なお，max-sum アルゴリズムは動的計画法のグラフィカルモデルへの応用と考えることもできる．
ここで考えたいタスクは，同時分布を最大にするベクトル xmax を求める問題として以下のように定義できる．

xmax = arg max p(x)xこのときの同時確率値は，次の式で与えられる．

(3.36)

62

第 3 章 モデリングの基礎

p(xmax ) = max p(x)x(3.37)

まずは，単純に成分に関する最大値演算を書き出してみると次のようになる．

max p(x) = max · · · max p(X)
x

x1

xM

(3.38)

ただし，M は変数の総数である．次に，同時分布を因子の積で表現して p(x)に代入する．ここで，最大値演算に関する以下の法則を用いると，積演算と和演算が交換可能になる．

max(ab, ac) = a max(b, c)(3.39)

対数をとって，以下のように書ける．

ln(max p(x)) = max lnp(x)xx(3.40)

3.1.1 行動データの正規化
人の行動を分析する際トリップに区分してから考えることが一般的である.PP 調査でもパーソントリップ調査と比較するためにトリップ単位で集計することが多い. 移動の状態を分析する上で基礎となる情報としては、1) 移動・滞在状態を分けること、2) 移動時にどの交通機関を利用しているか、3)
移動時どの経路を移動しているかなどが大きな要素となる. 例えば図 3.3 の時間と位置の生データが得られているとする. 図 3.4 は図 3.3 を平面上にプロットしたものである. この状態ではラベルづけされていないため、移動・滞在、交通機関がわからない. また取得データには誤差が含まれているため、正確にどの経路を移動したのかはわからない. 従って図 3.5 のようにまず移動、滞在識別を行うことで、トリップを作成する. 次に図 3.6 のように移動中の交通機関を識別する. またトリップ全体でどの経路を利用したかを図 3.7 のように特定する.

(a) 移動滞在識別
GPS データを移動滞在識別に用いる場合、そのままでは移動中の点 (移動点) なのか、滞在中の点 (滞在点) なのか不明である. そこで 1) 前の時刻までの測位店重心からの距離を計算して一定距離以上離れたら移動と判定する方

3.1

図 3.3

タイムスペース図

図 3.5

移動滞在識別

行動データマイニング

63

図 3.4 平面上にプロットした図

図 3.6

交通機関識別64

第 3 章 モデリングの基礎

図 3.7

交通機関識別法や、2) 測位点間の速度を計算して一定速度以上であれば移動と判定する手法によって移動滞在の識別を行うことができる. ここでは距離を用いた移動滞在識別手法について簡単に説明する [1].

図 3.8

移動滞在識別取得された GPS データが 2 次元座標上にあり、図 3.8 のように時間的に連続する 2 つの位置データ Pt 、Pt+1 があるとする. このときデータ間の距
離は式 (3.41) のようになる.

3.1

d=

行動データマイニング√
2
2(xt+1 − xt ) + (yt+1 − yt )65

(3.41)

この距離があらかじめ与えた閾値 D 内にあるとき、時刻 t t + 1 の間に移動はないものとみなす. 時刻 t までの識別が終わっているとして、時刻 t + 1の識別を考える. 時刻 t + 1 の点識別は、時刻 t の点が滞在点であるか、移動点であるかによって異なる. また時刻 t が滞在点であることは時刻 t − 1 と時
刻 t の条件によって決められるが、移動点であるかどうかは時刻 t と t + 1との関係を計算しない限り決まらない. 従って、移動点には暫定的な移動点と、確定した移動点の 2 種類があることになる. 以下にそれぞれの場合のアルゴリズムを示す.
時刻 t の点が滞在の場合ステップ 1
時刻 t 以前の時間的に連続する N 個の点 (時刻 (i = t − N + 1, .., .t))は滞在点であるとする. それらの点の重心位置 (x̄, ȳ) と識別対象である時刻 t + 1 の点 (xt+1 , yt+1 ) との距離を式 (3.42) により求める.
1x̄ = N1ȳ = Nt∑i=t−N +1
t
∑

xi
yi

(3.42)

i=t−N +1

ステップ 2-1
距離が閾値 D 以内なら時刻 t + 1 の点は時刻 t までの点と同じ位置での滞在点とする. 平均値、及び、連続する点の数を更新し N = N + 1とする.
ステップ 2-2
距離が閾地を越えているなら、時刻 t + 1 の点を暫定的な移動点とする.
時刻 t の点が暫定的移動点の場合
ステップ 1
時刻 t の点と時刻 t + 1 の距離を式 (3.41) により求める.

66

第 3 章 モデリングの基礎

ステップ 2-1
距離が閾値 D 以内なら時刻 t の点と時刻 t + 1 の点はいずれも同一地
点での滞在点である. このとき 2 点の平均位置座標を求めておく.
ステップ 2-2
距離が閾値 D 以上なら時刻 t の点と時刻 t + 1 の点は異なる位置の点であり、時刻 t の点は真の移動点とし、時刻 t + 1 の点は暫定的な移動点とする.
補足制約
上記のルール以外にも時間によって移動滞在を分ける以下の制約条件を設ける.

1: 滞在時間制約
短時間の滞在を除外する. 同一地点に 5 分以上滞在しない限り、滞在とはみなさない.

2: 移動時間制約
短時間トリップを除外する. 移動時間が 5 分以上でなければトリップとみなさない.
以上のような流れで測位データを下に移動滞在判別をすることでトリップ候
補を自動的に作成することができる.

(b) 交通機関識別
移動滞在を識別し、トリップ候補の作成が完了したのち交通機関の自動識別を行う. 交通機関自動識別の手法するために扱うデータとしては 1)GPS
データによる識別 2) 加速度データによる識別などがある.1) の GPS による交通機関識別では移動滞在識別のときと同様に閾値を各交通機関ごとに閾値を設けることで識別を行う.2) の、加速度データによる識別では交通機関ごとに現れる特徴を利用して識別を行う. それぞれの識別手法についてについて詳しく説明していく. またここでは”滞在”、”徒歩”、”自転車”、”車”の 4 交通
機関識別について扱うこととする.

GPS データによる識別GPS データによる識別では GPS データから各点の間の移動速度に交通機関

3.1

行動データマイニング

67に閾値を設けて交通機関識別を行う. 各交通機関ごとに速度の発生確率に特徴があるのでそれをもとに識別を行う. 以下に具体的なアルゴリズムを示す.
ステップ 1
最も識別しやすい徒歩トリップから識別を開始し、ステップ 1 へ戻ってきたら、自転車→車の順でトリップを抽出する.
ステップ 2
移動滞在識別により作成された各トリップ候補について各測位データ
間の移動速度が一定値以下になるまでトリップを再分割.
ステップ 3ステップ 2 終了時にトリップ開始時刻、終了時刻、交通機関を確定する. 残りのデータでステップ 1 へ戻り、車まで終われば終了.
加速度データによる識別
加速度による識別も各交通機関ごとの加速度の特徴をもとに識別を行う。図3.16、3.17、3.18 は交通機関ごとの 3 軸方向 (進行方向、左右方向、鉛直方向) 加速度の推移を表わしたものである. 徒歩、自転車、車の順に 3 軸加速度の分散が大きいのがわかる. このように交通機関ごとに異なった加速度の特徴をもつので、これを利用して識別を行う. ここでは SVM、HMM を用いた交通機関識別の流れについて説明する.

1: SVM(サポートベクターマシン)SVM は 2 クラス識別問題を扱うための手法であり、文字認識や音声認識などに幅広く用いられている.SVM の特徴としては主に以下の 3 点が挙げられる.

1. マージン最大化を基準とすることで、高い判別性能を実現.
2. カーネル関数を用いた柔軟なモデリング．
3. パラメータ推定が凸 2 次最適化問題として定式化されるため、高速な最適化アルゴリズムが利用可能．
式 (3.43) で表わされる線形モデルの 2 値分類問題を解くことでその特徴について説明を行う.

y (x) = wT ϕ (x) + b(3.43)

68

第 3 章 モデリングの基礎

ここで ϕ(x) はある固定された特徴空間変換関数である。訓練データは x1 , ..., xNの N 個の入力ベクトルとそれに対応した出力値 t1 , ..., tN ∈ {−1, 1} からなり、未知データ x は y(x) が正か負のどちらかによって分類が行われる. つ
まり式 (3.44) が成立する.

tn = +1 → y (xn ) > 0
tn = −1 → y (xn ) < 0(3.44)

したがって必ず式 (3.45) が成立する.

tn y (xn ) > 0(3.45)

今、図 (3.9) のように訓練データを正確に分類する解が y1 、y2 のように複数存在する場合、汎化誤差が最も小さくなるような解を求めることが最も望ましい.SVM では、マージンという概念を用いてそのような解を求める. ここでマージンとは境界面と訓練データの間の最短距離を表わす. このマージンを最大化する分類境界を選択する. 分類境界面 y(x) = 0 から点 x までの距離
は |y (x)|/∥w∥ で表わされる. また今式 (3.45) がすべての n について成り立つので、分類境界から点 xn までの距離は式 (3.46) のように表わされる.

tn y (xn )
∥w∥

=

(
)
tn wT ϕ (xn ) + b(3.46)∥w∥

今求めたいのはマージンを最大化するパラメータ w と b である. 従って、解は式 (3.47) を解くことで得られる.

{
arg max
w,b

1
∥w∥

[ (
)]
min tn wT ϕ (xn ) + b}(3.47)

次に直接この問題を解くのは複雑なので、より簡単な形へ変形することを考える. まずパラメータを w → κw、b → κb のように定数倍しても、点xn から分類境界までの距離は変化しない. 従って定数倍することで、境界に最も近い点について式 (3.48) が成立する.

()
tn wT ϕ (xn ) + b = 1(3.48)

3.1

図 3.9

マージンの概念

行動データマイニング

69

図 3.10 訓練データが線形分離不可能な場合

このスケーリングの下ではすべてのデータについて制約式 (3.49) が成り立つ.

()
tn wT ϕ (xn ) + b ≥ 1(3.49)

このとき式 (3.47) は式 (3.50) の最適化問題に帰着する.

arg min
w,b

1
2

∥w∥

2

(3.50)

ここまでの議論では、訓練データが特徴空間 ϕ (x) において線形分離可能であり、そこから得られる SVM はもともとの入力空間 x において訓練デー
タを完全に分離すると仮定してきた. しかし、図 3.10 のように実際にはそのような仮定が満たされないことも多い. このような状況に対応するためには、
一部の訓練データの誤分類を許すように SVM を修正する必要がある. 誤って分類したデータについては無限大のペナルティを与え、正しく分類したデー
タにはペナルティを与えない誤算関数を用いる. この誤差関数を、データ点がマージン内に侵入した場合にはマージン境界からの距離に応じたペナルティを与えることで、マージン境界の間違った側に存在することを許すように定式化を修正する. 具体的にはまずスラック変数 ξn ≥ 0 (n = 1, ..., N ) を導入する. スラック変数は各訓練データごとに定義される変数で、データが正
しく分類されかつマージン境界の上または外側に存在する場合は ξ = 0、それ以外の場合には ξn = |tn − y (xn )| で与えられる. 従って、ちょうど分

70

第 3 章 モデリングの基礎

類境界 y = (x) 上にあるデータについては ξn = 1、誤分類されたデータについては ξ > 1 が成り立つ. このスラック変数を用いて識別関数式 (3.49) を式 (3.51) のように修正する.

tn y (xn ) ≥ 1 − ξn(3.51)

ただし、ξ ≥ 0 である. 図 3.11 のように ξ = 0 が成り立つデータは正しく分類されており、かつマージン境界上、あるいは正しい側に存在する. また

0 < ξn ≤ 1 となるデータ点はマージン内部にあるが正しく分類されている.
そして ξ > 1 となるデータ点は分類境界に対して誤った側にあり、従って誤分類されている. これはハードマージンの制約のソフトマージンへの緩和といわれることがある. このとき式 (3.50) に対する最小化すべき最適化問題は式(3.52) のようになる.
N
∑

C

n=1

図 3.11 ソフトマージン

ξn +

1
2

∥w∥

2

図 3.12(3.52)

カーネルトリック

ここで C > 0 はスラック変数を用いて表わされるペナルティとマージンの大きさ間のトレードオフを制御するパラメータである.C → ∞ の極限においては 分離可能なデータに対するハードマージン問題と同じになる. 以上のようなソフトマージンを用いることで、線形分離不可能な訓練データに対

3.1

行動データマイニング

71

しての識別も可能になる. しかし、本質的に非線形で複雑な場合、このソフト
マージンを用いてもうまく訓練データを識別することができない. そこでカー
ネルトリックという手法を用いることでこの問題を解決する. 非線形な写像

Φ を用いて入力空間をより高次の空間に写像し、写像先の空間で線形分離を行うことで分離を容易にする. 写像 Φ を用いた計算は本来であれば元の次元
より高次元でのベクトル計算を行う必要があるが、カーネルトリックを用いることで高次元でのベクトル計算を避けることが可能となる. カーネルトリックとは写像 Φ に関する計算 Φ (x) · Φ (x′ ) をカーネル関数 K (x, x′ ) を用いて置き換えることで写像 Φ の球解を避け、カーネル関数のみの計算に変換することである. カーネル関数は複数ある. ここでは代表的なカーネルの一つである RBF (RadialBasisF unction) カーネルを式 (3.53) に示す.

()
2
K (x, x′ ) = exp −γ × |x − x′ |(3.53)

ここで、x、x′ はベクトル、γ はパラメータを表わし、RBF カーネルを用いる場合はパラメータ γ に適切な値を与える必要がある.

SVM を用いた交通機関識別
SVM は訓練あり学習を用いる識別手法の一つである. すなわち訓練データから識別クラスごとの特徴量を学習することで、未知データの識別行う. 具体
的には交通機関ごとの 3 軸方向加速度データ (それぞれの方向について最大、最小、合成、平均加速度データが取得されている) から交通機関ごとの特徴を学習することで識別を行う. 表 3.2 に、入出力データを、図 3.13 に交通機
関識別のイメージを示す. 2: HMM(隠れマルコフモデル)
表 3.1
入力データ

出力データ

交通機関識別入出力データ

3 軸方向加速度 (最大、最小、合成、平均)
運動強度
歩数
交通機関 (滞在、徒歩、自転車、車など)

HMM(隠れマルコフモデル) は不確定な時系列データをモデル化するための有効な統計的手法である. ここでは HMM の概要とモデルの推定、最適状態系列の推定の説明を行う [2]、[3]。HMM は出力シンボルによって一意に状態

72

第 3 章 モデリングの基礎

図 3.13 サポートベクターマシン交通機関識別イメージ

図 3.14

HMM イメージ図

遷移先が決まらないという意味での非決定性確率有限オートマトンとして定義される. 出力シンボル系列が与えられても状態遷移系列は唯一には決まらな
い. 観測できるのはシンボル系列だけであることから,”隠れ”マルコフモデルと呼ばれる. 図 3.14 に HMM のイメージを示す.

HMM はパラメータとして状態遷移確率、シンボル出力確率、初期状態確率を持つ. シンボル出力確率の計算方法によって離散型 HMM と連続分布型

HMM に分かれる. またシンボル出力確率が状態で出力される Moore マシン

3.1

行動データマイニング

73と状態遷移で出力される Mealy マシンに分類できる. 以下では Mealy タイプの離散型 HMM について述べる. ただし、Moore タイプと Mealy タイプは相互に変換可能である. 説明で用いる変数を表 3.2 に定義する. ここで HMM の
表 3.2

To1 , o2 , ..., oT
NL
S = {s}sti, jυ = {υ1 , υ2 , ..., υL }変数の定義
観測系列の長さ
観測系列
状態数
観測シンボル数
状態集合
時刻 t における状態 (番号)
状態番号
出力可能なシンボル集合

要素を正確に定式化し、観測系列がどのように生成されるか説明を行う.HMM
は以下のように特徴づけられる.

1. モデル (実際の状態) の数 N ．それぞれの隠れた実際の状態をそれぞれS = {S1 , S2 , · · ·, SN } と表わし、時刻 t における状態を qt と表わす.
2. 観測シンボルの数 M . これはモデル化されたシステムが出力する状態を表わす. それぞれの観測シンボルを V = {v1 , v2 , · · ·, vN } と表わす.

3. 状態遷移確率分布 A = {aij }．ここで aij は式 (3.54) のようになる.
aij = P [qt+1 = Sj |qt = Si ]1 ≤ i, j ≤ N(3.54)

ある状態からすべての状態へ遷移できるモデルを想定する場合すべてのi、j について aij である. ある状態から特定の状態にしか遷移できないモデルであれば、aij = 0 となる場合がある.

4. 状態 j のときの観測シンボル確率分布 B = {bj (k)} ただし式 (3.55)が成り立つ.

bj (k) = P [vk at t|qt = Sj ]

1≤j≤N
1≤k≤M

5. 初期状態確率 π = {πi } ただし、式 (3.56) が成り立つ.

(3.55)

74

第 3 章 モデリングの基礎

πi = P [q1 = Si ] ,1≤i≤N(3.56)

N 、M 、A、B と π が与えられたとき、HMM が式 (3.57) のような観測系
列を生成したとする.

O = O1 O2 · · · OT

(3.57)

それぞれの観測シンボル Ot は V の一つであり、T は観測系列の数を表わす. このとき式 (3.57) で表わされる観測系列が生成される過程は以下のよう
になる.
過程 1
初期状態確率 π に従って初期状態 q1 = Si を選択する
過程 2

t = 1 とする
過程 3
状態 Si のときの観測シンボル確率分布 bi (k) に従って Ot = vk とする.
過程 4
状態 Si のときの状態遷移確率分布 aij に従って新しい状態 qt+1 = Sjに遷移する.
過程 5
もし、t < T ならば t = t + 1 とし、過程 3 へ、そうでなければ過程を終了する
以上の過程は観測系列の生成や観測系列が HMM によって以下に生成されたかに利用される.HMM を規定するためには実際の状態と観測状態の数 N 、M 、観測系列、確率値 A、B 、π が必要であることがわかる. そこで簡単た
めに HMM を式 (3.58) のように略記する.

λ = (A, B, π)(3.58)

3.1

行動データマイニング

75

図 3.15 3 状態 left to right HMM

ここで HMM の例として状態遷移が一定方向に進む left to right モデルの例を図 3.15 に示す. 図 3.15 において HMM は 3 つの状態から構成され、2 種類のラベル a、b のみからなるラベル系列を出力する. 初期状態確率は
π1 = 1.0、π2 = 0、π3 = 0、最終状態を S3 とし、図 3.15 のような遷移のみを行うとする. 図 3.15 において 0.4 など遷移アークに添えられている数
字は状態遷移確率を意味し、[] 内の数字の上段はラベル a の出力確率、下段
はラベル b の出力確率を意味する. 例えば状態 S1 を例にとると、S1 から S1
自身に 0.4 の確率で遷移し、遷移の際に確率 0.7 でラベル a を出力し、0.3 の確率で b を出力する. 他の状態、遷移についてもどうようである. ここでラベル系列が ”aab” を出力する確率を考える. この HMM で許される状態系列において ”aab” を出力する可能性のあるのものは、S1 → S1 → S2 → S3 、

S1 → S2 → S2 → S3 、S1 → S1 → S1 → S3 である. それぞれについて ”aab” を出力する確率は以下のようになる.

S1 → S1 → S2 → S3

0.4 × 0.7 × 0.3 × 1.0 × 0.7 × 0.5 = 0.0294

S1 → S2 → S2 → S3

0.3 × 1.0 × 0.3 × 0.3 × 0.7 × 0.5 = 0.00945

S1 → S1 → S1 → S3

0.4 × 0.7 × 0.4 × 0.7 × 0.3 × 0.2 = 0.004704

従ってこの HMM が ”aab” を出力する確率は 3 つの合計で以下のように
なる.

76

第 3 章 モデリングの基礎

0.0294 + 0.00945 + 0.004704 = 0.038154
HMM では状態系列に意味を持たないが、最尤の経路を推定することはできる. この例では ”aab” を出力する可能性が最も高い状態系列は、先ほどの計算より、S1 → Sd1 → S2 → S3 とわかる.

HMM の基本問題
HMM を現実の問題に適用することを想定したとき重要な基本問題として以下の 3 つが挙げられる.
問題 1 観測系列 O = O1 O2 · · · OT とモデル λ = (A, B, π) が与えられたとき、観測系列が与えられた HMM モデルでどれくらいの確率で観測されるか問題 2 観測系列 O = O1 O2 · · · OT とモデル λ = (A, B, π) が与えられたとき、どのよううに最も可能性の高い状態遷移系列を推定するか.
問題 3 モデルパラメータ λ(π, A, B) をどのように調整するのか.
問題 1 は評価の問題で、つまりモデルと観測系列が与えられたときモデルによって観測系列が生成された確率をどのように計算するかという問題である.
あるいは、与えられたモデルが観測された系列にどれくらいマッチしているかという問題ととらえることができる. 問題 2 はモデルの隠れた部分を推定する問題で、つまり正しい系列を探索する問題である. 問題 3 はモデルが最もマッチするよなパラメータを求める問題である. モデルパラメータを修正するための系列を訓練系列と呼ぶこととする. 訓練問題は観測された訓練データから最も適したパラメータを採用する問題なので、HMM を適用する上で重要になってくる.
モデル評価問題与えられたモデル λ において観測系列 O = O1 O2 · · · OT 、の観測確率、

P (O|λ) を計算する. 最も単純な方法は、図 3.15 の例で示したように長さ Tの系列の状態を数え挙げる方法がある. 今、ある一つの状態系列が式 (3.59)のように表わされるとする.

Q = q1 q2 · · · qT

(3.59)

3.1

77

行動データマイニング

ここで q1 は初期状態である. 式 (3.59) で表わされる状態系列のとき観測系列

O が観測される確率は式 (??) のようになる.
P (O|Q, λ) =T∏

P (Ot |qt , λ)(3.60)

t=1

ただし、各観測同士は互いに独立だと仮定する. 式 (??) を書き下すと式 (3.61)が成り立つ.

P (O|Q, λ) = bq1 (O1 ) · bq2 (O2 ) · · · bqT (OT )(3.61)

状態系列 Q の確率は式 (3.62) のようになる.

P (Q|λ) = πq1 aq1 q2 aq2 q3 · · · aqT −1 qT(3.62)O と Q が同時に起こる確率は式 (3.61)、(3.62) を使って式 (3.63) のように表わされる.

P (O, Q|λ) = P (O|Q, λ) P (Q|λ)(3.63)

O が得られる確率は式 (3.63) をすべてのあり得る状態の系列について足し合
わせたものなので式 (3.64) のようになる.

∑

P (O|λ) =P (O|Q, λ) P (Q|λ)∑πq1 bq1 (O1 ) aq1 q2 bq2 (O2 )

all Q

=

q1 ,q2 ,···,qT(3.64)· · · · aqT −1 qT bqT (QT )以上の式展開の解釈は次のようになる. まず時刻 t = 1 では確率 πq1 で状態

q1 になり、確率 bq1 (O1 ) でシンボル O1 を生成する.t = t+1 にし (t = 2)、
確率 aq1 q2 で状態 q1 から状態 q2 へ推移する. また確率 bq2 (O2 ) でシンボル O2 を生成する. このような過程を確率 aqT −1 qT で状態 qT −1 から状態 qTに遷移し、シンボル OT が確率 bqT (OT ) で生成されるまで繰り返す.
隠れ最適状態推定問題
この問題は観測系列 O = {O1 O2 · · · OT } が与えられたときに最適状態

78

第 3 章 モデリングの基礎

Q = {q1 q2 · · · qT } を探索する問題である. この隠れ最適状態を推定するために利用されうる手法の一つに Viterbi アルゴリズムがある. ここで Viterbi
アルゴリズムに必要な量を式 (3.65) により定義する.

δt (i) =max

q1 ,q2 ,···,qt−1

P [q1 q2 · · · qt = i, O1 O2 · · · Ot |λ](3.65)

δt (i) は状態 Si で終わる時刻 t までの観測を説明する最もその確率が大きい
パス (遷移系列) である. このとき式 (3.66) が成り立つ.

[]

δt+1 (i) = max δt (i) aij · bj (Ot+1 )i(3.66)

実際に状態系列を得るためには式 (3.66) を最大化する状態を各 i、j について把握していなければいけない. これを配列 ψt (j) を経由して行う. 具体的なアルゴリズムを以下に示す.
ステップ 1: 初期化

δ1 (i) = πi bi (O1 ) ,1≤i≤Nψ1 (i) = 0(3.67)(3.68)

ステップ 2: 反復

δt (j) = max [δt−1 (i) aij ] bj (Ot ) , 2 ≤ t ≤ T
1≤i≤N

1≤j≤N

ψt (j) = arg max [δt−1 (i) aij ] ,(3.69)

2≤t≤T

1≤i≤N

1≤j≤N

ステップ 3: 終了(3.70)

3.1

P ∗ = max

行動データマイニング

79[δT (i)](3.71)

qT∗ = arg max [δT (i)](3.72)

1≤i≤N

1≤i≤N

ステップ 4: 状態系列のバックトラッキング()∗,
qt∗ = ψt+1 qt+1

t = T − 1, T − 2, · · ·, 1

(3.73)

ステップ 4 で求めたものが最適な経路 (状態系列) となる.
パラメータ調整問題
パラメータ調整問題は HMM の中でも最も難しい問題である. パラメータを調整することで、観測系列が観測される確率が最も大きくなるようにする.Baum-

Welch アルゴリズムを用い訓練データから局所的に観測系列の観測確率を最大にするパラメータを選択することができる.Baum-Welch アルゴリズムを説明するに当たりまず変数の定義を行っておく.ξt (i, j) を観測系列において、与えられたモデルで時刻 t のときに状態 Si でかつ時刻 t + 1 のときに状態 Sj
である確率であるとする. つまり式 refdmeq51) のように表わされるとする.

ξt (i, j) = P (qt = Si , qt+1 = Sj |O, λ)(3.74)

この ξt (i, j) は式 (3.75)、(3.76) であるとき式 (3.77) のように表わされる.

αt (i) = P (O1 O2 · · · Ot , qt = Si |λ)(3.75)βt (i) = P (Ot+1 Ot+2 · · · OT |qt = Si , λ)(3.76)

ξt (i, j) =αt (i)aij bj (Ot+1 )βt+1 (j)P (O|λ)αt (i)aij bj (Ot+1 )βt+1 (j)
= ∑
N ∑
N
αt (i)aij bj (Ot+1 )βt+1 (j)(3.77)

i=1 j=1

ここで時刻 t で状態 Si である確率 γt (i) を式 (3.78) で定義する.

γt (i) = P (qt = Si |O, λ)このとき ξt (i, j) と γt (i) は式 (3.79) のような関係がある.

(3.78)

80

第 3 章 モデリングの基礎

N∑

γt (i) =ξt (i, j)(3.79)

j=1

時刻 t に対して γt (i) を足し合わせたT∑−1

γt (i) は状態 Si から遷移する

t=1と期待される数になり、また ξt (i, j) を同様に時刻について足し合わせたT∑−1

ξt (i, j) 状態 Si から状態 Sj に遷移すると期待される数になる. 以上のt=1ような関係式を用い、パラメータ π 、A、B の再推定を行う. 再推定された()

パラメータを λ̄ = Ā, B̄, π̄ とすると式 (3.80)〜(3.82) のようになる.

π̄i = γ1 (i)T∑−1(3.80)

ξt (i, j)

t=1
T∑
−1

(3.81)
γt (i)t=1
T∑

γt (j)

t=1
s.t Ot =vk
T∑(3.82)

γt (j)

t=1

再推定された λ̄ の評価は以下のように行う.

1. λ̄ = λ → (局所的な収束状態)()
2. P O|λ̄ > P (O|λ) →シンボル系列 O を出力するより最適なモデル
λ を推定
Baum-Welch アルゴリズムは、学習データの尤度を最大にするようにパラメータを学習する. ただし、基本的には gradient 学習によるパラメータ収束の学習方法であるためｎ、local minimum の方向にしか学習が進まない. そのため初期値が重要になる.

HMM を用いた交通機関識別HMM ではシステムがパラメータ未知のマルコフ過程 (未来の挙動が、過去の挙動とは無関係に現在の値だけで決定される性質) であると仮定そ観測可能な3.1

行動データマイニング

図 3.16 3 軸加速度推移: 徒歩 (Hato[4] から引用)

図 3.17 3 軸加速度推移: 自転車 (Hato[4] から引用)

81

82

第 3 章 モデリングの基礎

図 3.18 3 軸加速度推移: 車 (Hato[4] から引用)

情報から未知パラメータを推定する. モデル中の状態数は”滞在”、”徒歩”、”
自転車”、”車”の 4 状態であり、各状態において出力シンボルの数は”移動速度”、”単位時間あたりの歩数”、”鉛直方向加速度振幅”の 3 つ (各時点における状態数を表わす) から 3 × 4 × 3 = 48 シンボルとした. 出力値の閾値の設定は、サンプルデータから各交通機関における各盲目の確率密度関数を重ね合わせ、境界となる値に設定する. 図 3.19 は交通機関別速度の発生確率であり、車の閾値に 3.4m/s、自転車と徒歩の閾値に 1.6m/s、徒歩と滞在の閾値に 0.6m/s を設定し 4 クラスにした. 図 3.20 は交通機関別単位時間あたり歩数の発生確率であり、徒歩と自転車の閾値に 1.7step/sec、自転車と徒歩及び滞在の閾値に 0.6step/sec を設定し 3 クラスにした. 図 3.21 は交通機関別
鉛直方向加速度振幅の発生確率であり、徒歩と自転車の閾値に 0.7G, 自転車と自動車の閾値に 0.43G、自動車と滞在の閾値に 0.10G を設定し 4 クラスにした.
次に Baum-Welch アルゴリズムの適用について説明する. このアルゴリズ
ムの目的は状態遷移確率 A を求めることである. 推定の際に必要となる初期状態確率 π をランダムに与え、初期状態遷移確率 A0 を等状態遷移の場合は0.99999、他状態遷移の場合は 0.00001 とした. これは交通機関が同じという状態が連続して起こりやすいからである. 初期シンボル出力確率 B0 は真値
データとの照合に得られる出力確率を与えた. この Baum-Welch アルゴリズ
ムは gradient 学習によるパラメータ収束を用いているため、局所最大にしか3.1

行動データマイニング

図 3.19 交通機関別、速度の発生確率 (大村 [2] から引用)

図 3.20 交通機関別、単位時間あたり歩数の発生確率 (大村 [2] から引用)

83

84

第 3 章 モデリングの基礎

図 3.21 交通機関別、鉛直方向加速度振幅の発生確率 (大村 [2] から引用)

学習は進まない. ここで扱っている HMM はすべての状態がすべての状態に推移する可能性のあるモデルであるので推定が困難である. 従って以下で説明する Viterbi アルゴリズムでは状態遷移確率 A を等状態遷移の場合は 0.99999、
他状態遷移の場合は 0.00001、シンボル出力確率 B を真値データとの照合により得られる確率とした. また Viterbi アルゴリズムの性質上確率 0 は 10−10とした.Viterbi アルゴリズムではモデル λ がシンボル系列 O を出力する最も可能性の高い状態遷移系列を推定し、その系列に対する尤度を求めるアルゴリズムである. モデル λ を決定するパラメータ群は Baum-Whelch アルゴ
リズムぬよるパラメータ推定が行えなかったので、初期状態確率 π は等確率で π = (0.25, 0.25, 0.25, 0.25) とする. このとき入力する系列データは一人の一日分のデータであるが、識別精度を考慮し、ある個人の調査期間の全データとする.

(c)

マップマッチング

マップマッチング基礎概念 高度道路交通システム (ITS) の分野において、
経路ガイダンス、道路利用者への課金、交通流制御などを行うために、あるいは利用者に事故やバス到着時刻情報を知らせるためには自動車の位置情

3.1

行動データマイニング

85

報が必要となる. 近年位置を正確に取得するための技術として GPS(Global

Positioning System) が目覚ましい技術発展を遂げてきた. また、ストリート
キャニオンなどの遮蔽物の多いところや地下トンネルのような GPS だけでは十分に位置を特定できない場所でも、方位情報や加速度情報を用いたデッドレコニング (自律航法) と呼ばれる技術によって正確な位置情報を取得できるようになってきた. このような技術を用いて取得された位置データを用い、ある時点で移動者がどの経路を利用しているのか、またその経路上のどこにいるのかを特定する技術がマップマッチングである. ここではマップマッチン
グ手法のレビューを行う.
マップマッチング基本ステップ
まずマップマッチングの基本的なステップについて説明を行う.
ステップ 1
道路ネットワークデータを準備する (図 3.22)
ステップ 2
道路ネットワーク上に測位位置データをプロットする (図 3.23)
ステップ 3
道路ネットワーク上のリンクやノードと測位位置データとの関係 (距離、角度) を定量化し通過リンクを特定する (図 3.24)
ステップ 4
リンク単位の所要時間、走行速度などパフォーマンス指標を各リンク
について算出 (図 3.25)
マップマッチングの基本的なステップは以上の通りである. リンクを特定する
際の定量化の方法や、測位方法の補正などによって様々なマップマッチング手法がある.
マップマッチング手法
この節では様々なマップマッチングを行う際の手法を紹介する. マップマッ
チングの精度向上のために様々なアプローチが考案されてきた. アプローチの仕方を分類すると 1) 幾何解析マップマッチング手法、2) 位相幾何解析マッ

86

第 3 章 モデリングの基礎

図 3.22 道路ネットワークの準備

図 3.23 測位データプロット

図 3.24 通過リンク特定

図 3.25 リンクごとの所要時間、走行速度
算出

3.1

行動データマイニング

87

プマッチング手法、3) 確率的マップマッチング手法、4) 高度な技術を利用したマップマッチング手法の 4 つに分類される. この節では各手法ごとに簡単な紹介を行う.

1): 幾何解析マップマッチング手法
幾何解析マップマッチング手法はリンクの道路ネットワークデータの形状のみを考慮した手法である. 幾何解析手法で最も利用される手法は point topoint といわれるマッチング手法である [6]. この手法は測位点をネットワー
ク上の最も近いノードにマッチングするという簡単な手法であり、実行も容易で計算もとても速い. 以下にステップを示す.
ステップ 1
測位点と道路ネットワーク上のすべてのノード (交差点など) の距離を計算する
ステップ 2
測位点を最も距離が小さいノードへマッチングする
ステップ 3
すべての測位点について計算すれば終了、そうでなければステップ 1 へ
しかしながら、この手法は道路ネットワークの形状に大きく影響を受ける. 例
えば図 3.26 をみると、距離の観点からだと、P はリンク A 上へマッチング
されることが適切であるにも関わらず、一番近いノードは B2 なので、リンク B 上にマッチングされる. したがって point to point 手法では、図 3.27のように多くのノードがあるネットワークの方が正確にマッチングされることになる. 別の幾何解析アプローチに point to curve といわれるマッチン
グ手法がある [6],[7]. これは先ほどの point to point マッチングと類似しているが、測位点をネットワーク上の最も近いリンクにマッチングする手法である. 図 3.28 を見ると、P とリンク A、B との距離は dA 、dB で表され、今 dA > dB なので、測位データ P はリンク A にマッチングされる. 以下にステップを示す.
ステップ 1
測位点とネットワーク上のすべてのリンク (道路) の距離を計算する

88

第 3 章 モデリングの基礎

ステップ 2
測位点を最も距離が小さいリンクへマッチングする
ステップ 3
すべての測位点について計算すれば終了、そうでなければステップ 1 へ

point to curve は point to point マッチングよりも正確にマッチングが可能であるが、図 3.29 において、移動者はリンク A を利用している可能性が大きいが、P2 はリンク B の方がその距離が近いのでリンク A ではなく
リンク B にマッチングされてしまう. これは実際の適用を考えた場合、密度の高いネットワークでは誤ったマッチングが増えてします可能性があることを示している. これら二つのマッチング手法に近いが、それより精度がよい手法として curve to curve マッチングがある [6],[7],[8]. 以下にステップを示す.
ステップ 1point to point マッチングにより、測位点をマッチングする候補となるノードを探索する
ステップ 2
候補ノードをつなぎ候補リンクを作成する
ステップ 3
候補ノードを結び候補経路を作成する
ステップ 4
測位点を結びその軌跡との距離の和を計算するステップ 5すべての候補経路のうち最も距離が小さい経路を真の経路に決定する
図 3.33 では水色の経路にマッチングされる.

最後に RRFroad reduction

ﬁlter) といわれる手法を紹介する [9]. この手法は、まず候補となる道路をすべて列挙して、間違っている道路を機械的に取り除いていく手法である. まず

curve to curve マッチングによって候補経路を探索する. その後測位位置

3.1

行動データマイニング

89

データと比較することでマッチングを行う. ここでは差分 GPS(VDGPS) が用いられる. これは正確に位置の分かっている場所を参照点として受信機をおき、そこでの観測データから未知点と参照点に共通な誤差から位置の補正を行うための手段である. ステップを以下に示す.
ステップ 1
測位点と道路ネットワークの上のリンクとの距離を計算し最も近い道路ネットワーク上に投影するステップ 2VDGPS を用いて位置を修正するステップ 3次の時点の測位点とステップ 2 で生成した修正点によって新たな点を生成するステップ 4ステップ 3 で生成した点を道路ネットワーク上に投影するステップ 5道路ネットワークが 1 つに定まれば終了、そうでなければステップ 1からステップ 4 を繰り返す

2): 位相幾何解析マップマッチング手法
位相解析マップマッチング手法は幾何的な情報に加え、リンクの接続性、連続性をもとにマッチングを行う [10]. このマッチングでは、連続する測位位置
データが道路ネットワーク上のどのリンク上にマッピングされるかを、距離だけでなく、リンクデータと位置データの交差角度、相対角度距離などから決定する (リンク決定フェイズとする). その後、速度と進行方向角度によってリンク上のどの位置にあるのかを決定する (位置決定フェイズとする). リン
クとリンク上の位置を決定したのち、測位データが現在のリンク上にあるのか、つまり交差点に入ったかどうかを確認する (リンク確認フェイズとする).
リンク決定フェイズ
リンクはいくつかの評価指標の重みつき総和で決定される. 利用する評90

第 3 章 モデリングの基礎

図 3.26 point to point: 誤ったマッチン 図 3.27 point to point: 正しいマッチン
グ例
グ例

図 3.28 point to curve マッチング

図 3.29 point to curve: 誤ったマッチン
グ例

3.1

行動データマイニング

91

図 3.30 curve to curve: ステップ 1 候補 図 3.31 curve to curve: ステップ 2 候補リ
ノード探索
ンク探索

図 3.32 curve to curve: ステップ 3 候補 図 3.33 curve to curve: ステップ 4 位置
ルート探索
データ軌跡と候補ルートの距離計算

92

第 3 章 モデリングの基礎

図 3.34 RRF: ステップ 1 測位データを道 図 3.35 RRF: ステップ 2 VDGPS を用い
路ネットワークに投影て位置を修正図 3.36 RRF: ステップ 3 次の時点での測 図 3.37 RRF: ステップ 4 再修正したデー
位データと修正データで新たな位置を生成 タを道路ネットワークに投影3.1

行動データマイニング

93

図 3.38 RRF: ステップ 5 繰り返し

価指標は 1) 測位位置データの進行方向とリンク方位角の差、2) 測位位
置データとリンクの近接性、3) 測位位置データのリンクに対する相対位置、の 3 つである.

1) 進行方向とリンク方位角の差
図 3.39 において P1 〜P3 がリンク 1 にあると決定されているとき、P4 がどのリンクにあるかを決定する. このとき候補となるのはリンク 2,3,4 である. 図中北 (上) 方向と進行方向がなす角度を

β 、北 (上) 方向と各リンクがなす角度を βi (i = 2, 3, 4) とする.
このとき ∆β = β − βi とし、β ′ を式 3.83 と定義する. このとき式 3.84 で表わされる W SH を進行方向とリンク方位角差の重みつきスコアとする. ただし、AH は正のパラメータである.∆β ′が小さくなれば、つまり、進行方向とリンク方位角の差が小さければ小さいほどこのスコアは大きくなり、そのリンクにマッチングされる確率が大きくなる. ここで W SH は進行方向とリンクと方位角の差の重みスコア、AH (> 0) は重みパラメータである.

94

第 3 章 モデリングの基礎

図 3.39 測位位置データの進行方向とリンク方位角差

図 3.40 リンクとの近接性: 最短距離◦◦
∆β, if − 180 ≤ ∆β ≤ 180 360◦ − ∆β, if ∆β > 180◦
′
∆β =
◦
◦

360 + ∆β, if ∆β < −180W SH = AH cos(∆β ′ )(3.83)(3.84)

3.1

行動データマイニング

952) 測位位置データとリンクの近接性
測位位置データのリンクに対する相対位置には二つの評価指標がある.1
つは測位位置データと候補リンクの最短距離 (位置データから候補リン
クへの垂線の長さ) で、もう 1 つが連続する 2 つの位置データを結んだ
ときの線と、候補リンクが交差する場合の角度である.

a) 最短距離
図 3.40 において P は測位位置データ、AB は道路ネットワーク
上のリンク、C は P からリンク AB に下ろした垂線の足を表わす. このとき、P C の大きさは式 (3.85) のようになる.

D=

x3 (y1 − y2 ) − y3 (x1 − x2 ) + (x1 y2 − x2 y1 )√22
(x1 − x2 ) + (y1 − y2 )(3.85)

この距離 D が小さくなれば測位データとリンクの近接性は大きくなるので、最短距離についての評価指標は式 (3.86) のようになる.

W SP D = AP/D

(3.86)

ここで、W SP D は測位データとリンクの近接性の最短距離についての評価指標、D(> 0) は測位データから道路ネットワーク上のリンクへの垂直距離、AP は重みづけパラメータを表わす.

b) 連続する 2 点を結ぶ線と候補リンクが交差する角度
連続する測位データ、P (xi−1 , yi−1 ) と P (xi , yi ) と道路ネット
ワーク上の候補リンク l が交差しているときその角度が θ(鋭角)のとき、連続する 2 点の測位データと候補リンクの角度に関する評価指標は式 (3.87) のようになる.

W SP I = AP cos (θ)(3.87)

ここで W SP I は交差角度に関する評価指標を表わし、もし連続する 2 点を結ぶ線と候補リンクが交差していなければ 0 とする.

96

第 3 章 モデリングの基礎

図 3.41 候補リンクに対する相対位置

3) 測位位置データのリンクに対する相対位置
図 3.41 において P (xi , yi ) は測位データ、α1 〜α4 はリンク 1〜リンク 4 からの測位データの相対位置を表わす角度とする. 図

3.41 において測位データはリンク 2 かリンク 3 にマッピングされるのが適当であると考えられるので、角度が小さければ小さいほどそのリンクにマッチングされる確率が高くなる. したがって候
補リンクからの相対位置についての評価指標は式 (3.88) のようになる.

W SRP = ARP cos (α)(3.88)W SRP は候補リンクに対する相対位置に関する評価指標、α ≤180 は測位データと最も近いノードを結んだ線と候補リンクがな
す角度、ARP は重みづけパラメータである.
すべてのの重み付きの評価指標を足し合わせたものを T W S とすると
式 (3.89) のようになる.

T W S = W SH + (W SP D + W SP I ) + W SRP(3.89)

3.1

97

行動データマイニング

図 3.42 位置決定フェイズ: 速度と進行方向から推定

この T W S によって道路ネットワークのどのリンク上にあるのかを決定する.
位置決定フェイズ
リンク上のどの位置にあるのかを進行方向と速度 v によって決定する.
図 3.42 において P t と P t+1 はあるリンク A 上の時刻 t、t + 1 の自動車の位置を表すとする.P t における速度が v であるから式 (3.90) が成り立つ.

∆Ei = (v ∗ 1) sin θ∆Ni = (v ∗ 1) cos θ(3.90)

したがって P t+1 における自動車の位置は式 (??) にようになる.

Ei+1 = Ei + ∆Ei
Ni+1 = Ni + ∆Ni

このようにリンク上の位置を決定する.

(3.91)

測位位置を決定するのにGPS/DR のセンサデータとデジタルロードマップ上の位置データを用いる手法もある. 図 3.43 において GPS/DR によって測位された点をP S 、測位された点を道路ネットワーク上のリンクに投影したものを

98

第 3 章 モデリングの基礎

図 3.43 位置決定フェイズ: デジタルマップと測位データ平均から推定

C P 、デジタルマップからの点を P map とする. リンクの両端のノード A、B の座標と測位データの座標を用いると C P の座標は式 (3.92)のように表わされる.
1 )+Ns (y2 −y1 )]+(y2 −y1 )(x1 y2 −x2 y1 )P Es = (x2 −x1 )[Es (x2 −x(x−x )2 +(y −y )2
2

1

2

1

2

1

2

1

s (y2 −y1 )]−(x2 −x1 )(x1 y2 −x2 y1 )P Ns = (y2 −y1 )[Es (x2 −x1(x)+N−x )2 +(y −y )2(3.92)
今P

map

とC

P

は独立であるので、自動車の位置は式 (??) のようにそれぞれの重みづけ平均で推定される.

(

Ê =
N̂ =

)

2
σE

S
2 +σ 2
σE
map
S

(

2
σN

S

2
2
+σmap
σN
S

(
)
2
σmap
Emap + σ2 +σP ES
2)( ES map )2
σmap
Nmap + σ2 +σ
P NS
2
NS

(3.93)map

2
2
Ê 、N̂ は推定された位置、σmap
は Emap に関する誤差分散、σEはS
2
ES に関する誤差分散、σNは NS に関する誤差分散を表わす.
S

リンク確認フェイズ
正しいリンクにマッチングしたのち、測位データが現在のリンクにあるのかどうかを確認する. 基準は以下の 2 点である.

3.1

行動データマイニング

99

1. 二つの連続する線 (つまり P (xi , yi ) と P (xi+1 , yi+1 ) を結んだ線と P (xi+1 , yi+1 ) と P (xi+2 , yi+2 ) を結んだ線)の進行方向角度が 45°より大きいか

2. 二つの連続する測位点の進行方向角度が 45°より大きいか、また α(測位店と最も近いノードを結んだ線とリンクのなす角度) が 90°より大きいか
もしこれらの条件が満たされれば、自動車は交差点に入ったとみなし、再びリンク決定フェイズに戻る.
以上のマップマッチングアルゴリズムを以下に示す.
ステップ 1
最初の測位点と最も近いノードを探索する
ステップ 2ステップ 1 で探索した最も近いノードを通るすべてのリンクを選択する
ステップ 3式 (3.89) で表わされる重みづけ式を用いて正しいリンクを決定する (このとき最初の測位点と次の測位点はこのリンクにマッチングされる)ステップ 4式 (3.93) を用いてリンク上の位置を決定する
ステップ 5∆β ′ < 45◦ 、α ≤ 90◦ を確認しもし満たされていれば式 (3.93) によりリンク上の位置を決定することを継続し、もし満たされていなければステップ 1 へ戻る
ステップ 6すべての測位点がマッチングされるまでステップ 5 を継続する

3): 確率的マップマッチング手法
確率的手法マップマッチング手法ではセンサーから取得された測位データ周辺に楕円や円形の信頼域を定義することが必要となる. この技術は DR セン

100

第 3 章 モデリングの基礎

サからの位置をマッチングさせるのに初めて導入された [11]. また GPS セン
サから取得されたデータに対しても利用され、エラー領域は GPS 測位デー
タの誤差分散から導くことが提案された [12]. ここではより精度の高い確率的
マップマッチング手法を紹介する [13]. このマップマッチング手法においても位相幾何マップマッチング手法と同様にマッチングするリンクを決定したのちにリンク上の場所を決定する. 確率的マップマッチングではリンクを決定する段階で確率的手法が用いられている. 図 3.44 において測位データが P だとすると、その周囲の楕円形がエラー域として定義され、このエラー域内にあるリンクが候補リンクとなる. 図 3.44 の場合候補リンクは”リンク 3”となる. エラー域は図 3.45 のように長方形を用いることもできる. このとき楕円の軸の長さ a、b は式 (3.94)、(3.95) のように表わされる.

v
(
)
√(
u
)2
u
t1
2
2
2
2
2
a = σ̂0 /2 σx + σy +σx − σy + 4σxy

(3.94)

v
(
)
√(
u
)2
u
t1
2
2
2
2
2
b = σ̂0 /2 σx + σy −σx − σy + 4σxy

(3.95)

また楕円形が図 3.44 中の上矢印からどれだけ傾いているかを表わす ϕ は
式 (3.96) のように表わされる.

(
ϕ = π/2 − 1/2 arctan2σxy
σx2 − σy2)
(3.96)

ここで σx 2 と σy 2 は測位位置の分散、σx y は共分散、σ̂0 は拡張係数を表わす. 拡張係数は GPS の軌道が安定していないこと、大気伝搬、様々な伝搬経路を通ること、受信者側のノイズなどを補うための定数で 3.03 がよいとされている [12]. 図 3.46 にこのリンク決定アルゴリズムのフローを示す. リンクが決定した後は位相幾何マップマッチングのときと同様にリンク上の位置を決定する.

4): 高度な技術を利用したマップマッチング手法
この節ではカルマンフィルタ、パーティクルフィルタ [15]、DDR などより高

3.1

図 3.44 信頼区域: 楕円形

行動データマイニング

図 3.45 信頼区域: 長方形

図 3.46 リンク決定アルゴリズム フローチャート

101

102

第 3 章 モデリングの基礎

度な概念を用いたマップマッチング手法を紹介する. カルマンフィルタ ここ
ではカルマンフィルタを用いたマップマッチング手法について説明する [16].
カルマンフィルタとはシステムや観測値に様々なランダムなノイズがのっているときにそれぞれのノイズの大きさに応じて適切な重みづけを行い、変化していくシステムの状態を精度よく推定するための手法である. 今、自立航法(DeadReckoning) によって取得したデータを j(i)、実際の走行軌跡の候
補データを m(i) としたとき、m と j を観測モデル式 (3.97) と状態モデル
式 (3.98) でモデル化する.

j (i) = M (i) X (i) + υ (i)(3.97)X (i + 1) = X (i)(3.98)

ただし M (i)、X(i) はそれぞれ式 (3.99) で表わされるとする.

M (i) = [m (i − 2) m (i − 1) m (i) m (i + 1) m (i + 2)]
T

X (i) = [x1 x2 x3 x4 ](3.99)
観測モデル、状態モデルでカルマンフィルタを用いて X(i) の推定値 X̂ (i)を求める. ここではその推定法については割愛する. 推定値と実測値の差であるイノベーション v(i) を式 (3.100) によって計算する.

v (i) = j (i) − X̂ (i)(3.100)カルマンフィルタのイノベーションは平均値 0 のガウス性白色雑音であるという性質がある. したがって複数の走行軌跡の候補が与えられた場合、それぞれの候補に対してイノベーションを求め、その白色度が最も高いものが真の走行軌跡であるとする. イノベーションの白色度を評価するものとして式(3.101) を用いる.
ρv =

1

N
−2
∑

(N − 1) σv2 k=0{v (k) − m} {v (k + 1) − m}

(3.101)

ただし m は v(i) の平均値である. 今状態モデル式 (3.98) は自動車が急に進
行方向を変化させないという仮定を置いている. したがって、急激に進行方向が変わった場合すなわち式 (3.102) が成り立つとき3.1|m (i) − m (i − k)| > ε,

行動データマイニング

k = −2, −1, 1, 2

103

(3.102)m (i − k) = m (i) として X(i) の推定を行う.
パーティクルフィルタ
パーティクルフィルタは、観測モデルとして非線形の任意のモデルを適用できるため、画像などの観測過程が複雑なコンピュータビジョン等で多く用いられている. 前節のカルマンフィルタは誤差が正規分布に従う場合用いることができないのでパーティクルフィルタを利用することができる. 自動車に関する位置データで分析を行う場合、まずマップマッチングを行う. しかし速度が低下している交差点近傍で位置補正を行う際は、移動時よりも誤差によるばらつきが有意になることがあるため、誤ったマップマッチングが発生しやす
い. そこでパーティクルフィルタを用いた位置補正をあらかじめ行うことで、交差点近傍でのマップマッチングの精度を向上させる. パーティクルフィル
タは、時刻 t における観測値 θt から状態ベクトル xt の事後確率である確率密度関数 P (xt |θt ) を推定する. この事後確率 P (xt |θt ) は、ベイズの定理
により以下の式 (3.103) に置換される.

P (xt |θt ) =P (yt |xt ) · P (xt |θt )P (yt |θt−1 )(3.103)∫P (xt |θt−1 ) =P (xt |xt−1 ) · P (xt−1 |θt−1 ) dxt−1(3.104)

まず t − 1 の状態をもとに抽出されたパーティクル群 P (xt−1 |θt−1 ) を選択し、対象領域に拡散する. これを事前分布とする. つぎに、運動モデル

P (xt−1 |θt−1 ) を用い、パーティクル群を一定の規則に基づいて誤差を加えた形で移動させる (拡散: 運動モデルを用いる (後述)). さらに、各パーティ
クルに対しその尤度を計算する (重み付け (後述)). 最後にリサンプリングを実施し尤度の高いパーティクルについては複数個サンプリングされ、低いものに関しては消滅する. このとき最も尤度の高いパーティクルが集中する位置を推定値として位置を確定する.

104

参考文献

運動モデル
パーティクルフィルタでは、推定対象の動作時性をモデル化したもの
で、モデルの精度を左右するものである. ここでは速度が低下している交差点近傍で位置補正を行う際、誤差項が小さくなるように定式化を行う. 現時刻 t における i 番目のパーティクルの状態ベクトル xt の運動モデル式は式 (3.105) のようになる.



at





at−1





α · V at−1 · ε




 





 
 =  bt−1  +  α · V bt−1 · ε 





  β·ξ
 0
 V at 
β·ξ
V bt
0

 bt

xit = 


(3.105)(a, b) は 2 次元平面座標、(Va , Vb ) は速度を表わす.α、β は分散の大()()
きさを表わすパラメータであり、ε∼N 1, σ 2 、ξ∼N 0, σ 2 は正
規乱数を表わしている.
重み付け
重み付けは、観測モデルを用い、パーティクルの重みを評価する. 観測
モデルは、あらかじめ状態空間内で観測される値の分布状況を測定し、これを学習データとしてモデル化したものである. 観測モデルでは、現時刻における観測値とパーティクルが示す状態量を入力する. 入力した
状態量において入力した観測値が観測される確率密度を求める. この確率密度から、そのパーティクルが示す状態量における移動体の存在尤
度を決定し、重みに反映する. ここでは各パーティクルから現時刻 t における観測値 yi までの距離の逆数を尤度とする.
図 3.47 にパーティクルフィルタの概念図を示す. 以上のような流れでパーティ
クルフィルタを用いて測位位置データを補正することでマップマッチングの精度を向上させる.

参考文献

図 3.47105

パーティクルフィルタ概念図 (三谷, 羽藤 (2009)[15] から引用)

参考文献[1] 高橋厚年, 朝倉康夫, 羽藤英二: 移動体通信システムによる交通行動データ特性に関する基礎的分析, 土木計画学研究 講演集, Vol.22, pp.413-416, 1999
[2] 大村朋之: 位置-加速度情報を用いた交通機関識別と移動負荷モデルの開発, 修士課程
学位論文, 東京大学, 2011
[3] 村上仁一, 確率的言語モデルによる自由発話認識に関する研究, 博士課程学位論文, 豊橋技術科学大学, 2003
[4] Eiji Hato: Development of behavioral context addressable loggers in the shell for
travel-activity analysis, Transportation Research PartC, Vol. 18, pp. 55-67, 2010
[5] C.M. ビショップ : パターン認識と機械学習 下 - ベイズ理論による統計的予測,[6] Bernstein, D., Kornhauser, A. : An introduction to map-matching for personalnavigation assistants,?,1996[7] White, C.E., Bernstein, D., Kornhauser, A.L.: Some map-matching algorithms
for personal navigation assistants. Transportation Research Part C, Vol. 8, pp.91108,200

106

参考文献

[8] Phuyal, B.: Method and use of aggragated dead reckoning sensor and GPS data
for map-matching, In: proceedings of the Institute of Navigation (ION) annual
conference, Portland, OR (20-27 September), 2002
[9] Taylor, G., Blewitt, G., Steup, D., Corbett, S., Car, A.: Road reduction ﬁltering
for GPS-GIS navigation, Transactions in GIS, Vol. 5, pp. 193-207, 2001[10] Quddus, M.A., Ochieng, W.Y., Zhao, L., Noland, R.B.: A general map-matchingalgorithm for transport telematics applications, GPS Solutions, Vol.7 , pp. 157167, 2003[11] Honey, S.K., Zavoli, W.B., Milnes, K.A., Phillips, A.C., White, M.S., Loughmiller, G.E.: Vehicle navigational system and method, United States PatentNo., 4796191, 1989
[12] Zhao, Y.: Vehicle Location and Navigation System. Artech House, Inc., MA.,
1997[13] Ochieng, W.Y., Quddus, M.A., Noland, R.B.: Map-matching in complex urban
road networks, Brazilian Journal of Cartography (Revista Brasileira de Cartograﬁa), Vol. 55(2), pp. 1-18.
[14] Najjar, M.E., Bonnifait, P.: A roadmap matching method for precise vehicle
localization using belief theory and Kalman ﬁltering. In: The 11th International
Conference in Advanced Robotics, Coimbra, Portugal, 2003[15] 三谷卓摩, 羽藤英二:パーティクルフィルタを用いた空間データの自動作成法, 土木計
画学研究・講演集, Vol.40, CD-ROM, 2009.
[16] Jo, T., Haseyamai, M., Kitajima, H.: A map-matching method with the innovation of the Kalman ﬁltering, IEICE Trans. Fund. Electron. Comm. Comput.
Sci. E79-A, pp. 1853-1855, 1996

3.1.2 行動データの分析
本節では，実際の行動データに対する応用例を示す．

(a) 活動パターン分析
sequential aliment method
活動パターンの分類は行動分析における一つの重要なトピックである．活動パターンの分析は主に空間情報や社会経済変数を用いて行われ，活動パター
ン生成シミュレータを用いた需要予測や交通シミュレーションの評価などに用いられる．従来より用いられてきた方法として活動パターン同士のユーク
リッド距離を用いて相互の関係を評価する方法があるが，この方法では活動の連続性を評価することができない．遺伝子配列の分析や言語処理では活動の連続性（系列）を考慮した類似度の指標が用いられているが，一次元の配列
パターンの評価であることが多く，交通分野において対象とする空間や時間，
社会経済変数など（図 3.48）多次元の情報とその相互作用を扱うには不十分でણ৿ঃॱ‫ش‬থ

参考文献

107

୎৿

୎৿

୎৿

୎৿

ᆯ௫

୫হ

ਁୱ$

ঞ५ॺছথ$

開始・終了時刻













ઐৢਃঢ়

ং५

ள੊

െన



ள੊

ং५

ણ৿ৃਚ

図 3.48

活動パターン例ある．Chang-Hyeon et al(2002) では，従来の一次元の sequential aliment

method を多次元の変数に拡張した分析手法を示した．
問題の定義

2 つの多次元アクティビティパターン（ソースパターン，ターゲットパター
ン）を比較する問題を考える．ここで，それぞれのアクティビティパターン
は，K 次元の属性を持ち，それぞれ m,n 次元のシーケンスを持つとする．このとき，ソースパターン s とターゲットパターン g を比較するとは，K×mのマトリクスと K×n のマトリクスを比較することになる．
′

s = s[s1 · · · sk · · · sK ](3.106)

sk = sk [sk0 · · · ski · · · skm ](3.107)

ただし，′

g = g[g1 · · · gk · · · gK ](3.108)

gk = gk [gk0 · · · gki · · · gkn ](3.109)

ただし，

ここで，sk ,gk はそれぞれソースパターン，ターゲットパターンの k 番目
配列の属性ベクトルである．
類似度の定義と計算例

108

参考文献

図 3.49 類似度コスト計算例（Chang-Hyeon et al(2002)）

図 3.50 類似度コスト途中計算例（Chang-Hyeon et al(2002)）

簡単のため，まず一次元の属性におけるアクティビティパターンの比較を考える．ここで，2 つのアクティビティの類似度を，2 つのパターンが同じ配列になるために必要な最小の配列要素の操作（削除，挿入，置換，同一）コストによって 2 つのアクティビティ間の類似度を定義する．この手法を sequential

aliment method と呼ぶ．
図 3.49 には，類似度コスト計算の例を示した．ここでは，s = s[ACB],g =g[ABC] のパターンを考え，s が g と等しくなるための最小コストをもつ操作パターンを求める．ただし，削除，挿入，置換，同一のコストをそれぞれ

1,1,2,0 とする．
具体的な計算は，動的計画法を用いて直前の操作と最小コストを記憶しながら図 3.49 の表を埋めていく．図 3.50 に途中計算例を示した．この例では，最適な操作パターンは（同一，挿入，同一，削除）となり，そのときのコス
トは 2 となる．
多次元への拡張
前節で述べた手法は，2 つのパターンを同一にするための最小オペレーショ
ンコストであり，それぞれのパタンの連続性（順序）を考慮した類似度計算手法といえる．ただし，単一の属性における系列類似度の評価であり多次元の属性間の相互依存性を考慮した系列類似度評価は行えない．

Chang-Hyeon

参考文献

109

図 3.51 2 次元の類似度概念（Chang-Hyeon et al(2002)）

et al(2002) は，この手法を多次元に拡張し，多次元の属性の各操作の set を一つの操作としてみなす方法を提案している．多次元の属性に関する操作の組
み合わせを一つの操作とみなし評価する方法といえる．具体的な例を図 3.51に示した．一次元の操作であれば，D の削除，挿入と 4 の削除，挿入は別々の操作となり，図 3.51 の例では合計 8 回の操作が必要となるが，提案手法では D,4 の削除，挿入をセットで一つの操作としてみなすため，全体で 6 回の操作となる．
このように集合での操作を認める場合，そのコストをどう設定するかが問題になる．従来法では，集合内の 2 つの操作コストの平均，中央値，最大値をとるなどの方法が存在するが，ここでは最大値を用いて以下のように定義
する．

max(β1 , β2 ) × wd

(3.110)

ただし，wd は削除コスト，βk は k 番目の属性の重みとする．このように置くことで，多次元の属性が同時に生じることを表現することができる．現実的には，すべての組み合わせを考慮することは非常に計算コストが高いため，
次に示すようなヒューリスティックな工夫によって計算コストを低減する．
多次元パターンの計算手法
具体的な多次元での計算手法を以下に述べる．

ステップ 1
各属性における一次元での最小コスト（操作パターンを求める）

110

参考文献

以下のように，各属性での操作とその位置を記憶する．

Ok,l = {p|p = d(i, k) ∨ i(j, k) ∨ s(i, j, k)}(3.111)

ただし，i ∈ {1, ..., m},j ∈ {1, ..., n},1 ≤ k ≤ K
ここで，Ok,l は k 番目の属性系列の l 番目の一次元オペレーションを表す．d(i, k) は i 番目の要素の削除，i(j, k) は j 番目への要素の挿入，s(i, j, k) は i 番目から j 番目の要素への置換を表す．
ステップ 2ステップ 1 で求めた一次元の属性での最小コスト操作パターンを多次元の属性間のアライメント操作へ統合する
ここで，Tk を k 番目の一次元のオペレーション集合とすると，その組み合わせ集合は T は以下のように与えられる．

T =

K
∏

Tk

(3.112)k=1

また，O t を t 次元のオペレーションセットとすると，以下のように書ける．

O t = {p|p = d(i, k) ∨ i(j, k) ∨ s(i, j, k)}(3.113)【具体例】
今，次のような 3 つの属性に関する一次元オペレーションセットがあ
るとする．
t
O1,l
′ = {p|p = d(2, 1), i(4, 1), s(5, 6, 1)}(3.114)tO2,l= {p|p = d(3, 2), i(4, 2), s(5, 5, 2)}(3.115)tO3,l′′ = {p|p = d(2, 3), i(4, 3), s(5, 5, 3)}(3.116)

このとき、t 番目の多次元オペレーションセットは以下のようにかける．

O t = {p|p = d(2, {1, 3}), d(3, {2}), i(4, {1, 2, 3}), s(5, 5, {2, 3}), s(5, 6, {1})}(3.117)

参考文献

111

図 3.52

アクティビティパターンへの適用例（Chang-Hyeon et al(2002)）

図 3.53

アクティビティパターンでの計算例（Chang-Hyeon et al(2002)）

ステップ 3
統合されたアライメント操作の総コストを求める．

Ct =

∑

cp

(3.118)

p∈Ot


max


 wd βk
cp =





wi βkmax
ws βkmax

if p = d(i, k)
if p = i(j, k)(3.119)if p = s(i, j, k)

アクティビティパターンへの適用例
ここでは，前節までの手法を用いて実際のアクティビティパターンに適用した例を示す．
まとめ
このように，異なる次元のアクティビティパターンを多次元に拡張された

112

参考文献

sequential aliment method によって統一的に評価することが可能となる．
このような手法は，アクティビティパターンの生成過程を記述したりシミュ
レータのアクティビティパタン生成モデルに利用できる可能性がある．アクティビティを符号化して操作することで本来連続的な属性も扱いやすくすることが可能である．個々の属性の重みや各操作のコスト設定に関してはさま
ざまな設定が可能であり，より記述力の高い重み設定を工夫する余地もあるといえる．

3.2 行動モデル
3.2.1 離散選択モデル(a) 効用と選択確率の定式化
効用の定式化
離散選択モデルでは，選択肢ごとに効用をある関数により与えた上で，選択肢間の効用の比較によって選択結果が導出される．選択肢の効用は直接観測することはできないが，効用は連続変数であり，かつ説明変数の線形和であると仮定する．効用関数は次となる．

Vin = β1 x1in + β2 x2in + · · · + βK xKin

(3.120)

ここで，Vin は個人 n の選択肢 i に対する効用を示す．xkin は個人 n の選
択肢 i に対する k 番目の説明変数，βk は k 番目の未知パラメータを示す．交通手段選択モデルの場合には，説明変数には所要時間・費用・乗り換え回数・交通手段までのアクセス時間等の交通手段別特性と，移動目的・性別・年齢等の個人ごとの移動・個人特性が用いられる．
説明変数 (xkin ) の観測誤差や個人個人の重み付けパラメータ (βk ) の差異
が実際は存在しており，式 (3.120) の効用関数は測定可能な内容による近似に過ぎないと考えられる．そこで，真の効用関数は測定可能な効用と不可能な効用による和で表す．測定不可能な効用を確率変数を用いて表現し，真の効用関数は次となる．