行動モデル夏の学校 2019

歩行者回遊行動モデリングのための最大エントロピー逆強化学習とEMアルゴリズム

山本俊行（名古屋大学）
日髙健・早川敬一郎・西智樹（豊田中央研究所）
薄井智貴（人間環境大学）

行動モデル夏の学校 2019

歩行者の回遊行動・目的地指向行動
移動=不効用，不効用最小化経路
・回遊行動
移動=正の効用，効用最大化経路
いくらでも効用を大きくできるため時間的な制約が必要

目標
入力：
行動系列：𝜻𝜻 = {𝜁𝜁1 , 𝜁𝜁2 , 𝜁𝜁3 }環境の特徴量：歩道，芝生，桜の木，POI
出力：
各特徴量の重み：𝜃𝜃 = [𝜃𝜃ww , 𝜃𝜃ch , 𝜃𝜃POI ]T
報酬関数：R 𝑠𝑠
R 𝑠𝑠 = 𝜃𝜃ww 𝑓𝑓ww + 𝜃𝜃ch 𝑓𝑓ch + 𝜃𝜃POI 𝑓𝑓POI

報酬関数は各特徴量の線形和

Illustrative example of theproblem (POI: point of interest)

1

行動モデル夏の学校 2019

発表の構成
１．逆強化学習
1-1. 逆強化学習の歴史
1-2. 最大エントロピー逆強化学習の導出
1-2-1. マルコフ決定過程（MDP）
1-2-2. Soft-max MDP
1-2-3. 最大エントロピー逆強化学習
1-2-4. 最大エントロピー逆強化学習の推定
1-2-5. Recursive logitモデルとの比較
1-3. 時空間制約下の方策
２．EMアルゴリズム
2-1. EMアルゴリズムとは？
2-2. EMアルゴリズムの導出
2-3. EMアルゴリズムのイメージ
３．歩行者回遊行動への適用

2

行動モデル夏の学校 2019

１．逆強化学習3

行動モデル夏の学校 2019

1-1. 逆強化学習の歴史

所与の報酬の下で最適な行動を見つける

報酬強化学習逆強化学習

最適な行動最適な行動系列を所与として報酬関数を求める（効用関数）• 逆強化学習はRussell(1998)によって提案• 様々な手法が提案されている
 線形計画法を用いた推定（Ng et al., 2000)
 エキスパートの行動軌跡を所与とし，エキスパートと同じような行動軌跡が得られる報酬関数をmax-margin法で推定（Abbeel and Ng, 2004)
 Abbeel and Ng(2004)の方法は条件を満たす複数の解候補があり，解の不定性を解決するため最大エントロピー原理を導入（Ziebart et al., 2008) その他にも相対エントロピー逆強化学習（Bouralias et al., 2011)，最大エントロピー深層逆強化学習（Wulfmeier et al., 2015)など様々な方法へ発展

4

行動モデル夏の学校 2019

1-1. 逆強化学習の歴史
逆強化学習の適用例

機械学習により「動物の行動戦略」を解読 －動物は何を報酬として行動しているのか－
http://www.kyoto-u.ac.jp/ja/research/research_results/2018/180515_1.html

逆強化学習による人の行動予測
Kitani et al.(2012)5行動モデル夏の学校 2019

1-2-1. マルコフ決定過程（MDP）
報酬関数
遷移モデル
行動の集合
状態の集合

マルコフ決定過程（MDP）はこの4つの組で定義される回遊空間の表現：グリッドs7s8s9s4s5s6s1s2s3𝑅𝑅(𝑠𝑠4 )𝑅𝑅(𝑠𝑠2 )状態の集合：グリッドの全てのセル𝑠𝑠1 ~𝑠𝑠9
行動の集合：上下左右遷移モデル：状態𝑠𝑠から行動𝑎𝑎を取ったときに状態𝑠𝑠𝑠に遷移する確率
報酬関数：新たな状態に遷移した結果，環境から受け取る値6

行動モデル夏の学校 2019

1-2-1. マルコフ決定過程（MDP）
最適な行動

⇒

累積報酬和を最大にする行動の結果として選択されたものと仮定する

図の例の場合，軌跡𝜁𝜁の累積報酬和𝑅𝑅 𝜁𝜁 は𝑅𝑅 𝜁𝜁 = 𝑅𝑅 𝑠𝑠1 + 𝑅𝑅 𝑠𝑠2 + 𝑅𝑅 𝑠𝑠5 + 𝑅𝑅 𝑠𝑠6 + 𝑅𝑅 𝑠𝑠9

ここでは，単純な累積和の代わりに割引累積報酬和を用いる
𝑅𝑅 𝜁𝜁 = 𝑅𝑅 𝑠𝑠1 + 𝛾𝛾𝑅𝑅 𝑠𝑠2 + 𝛾𝛾 2 𝑅𝑅 𝑠𝑠5 + 𝛾𝛾 3 𝑅𝑅 𝑠𝑠6 + 𝛾𝛾 4 𝑅𝑅 𝑠𝑠9

ただし，𝛾𝛾 ∈ [0,1]は割引率割引率を考慮する理由としては，・値が発散しない・将来の不確実性を考慮

逆強化学習は，このような累積報酬を最大にする行動が観測と合うような報酬関数を見つける方法

s7s8s9s4s5s6s1s2s3

グリッド空間上の行動
ζ = {𝑠𝑠1 , 𝑠𝑠2 , 𝑠𝑠5 , 𝑠𝑠6 , 𝑠𝑠9 }7

行動モデル夏の学校 2019

1-2-1. マルコフ決定過程（MDP）
ここで，状態価値関数V(s)と行動価値関数Q(s, a)を定義する

状態価値関数（状態sにいるときの期待累積報酬和）

行動価値関数（状態sにおいて行動aをとるときの期待累積報酬和）

ベルマンの最適性原理より以下の関係が得られる
次の状態の報酬

次の状態の価値（＝将来の報酬和）𝜋𝜋 ∗このときの最適な行動方策は以下のように書ける

最も行動価値が高い行動𝜋𝜋 ∗𝑅𝑅(𝑠𝑠, 𝑎𝑎)

𝜋𝜋 ∗ここでのモデルは決定論的なモデル

8

行動モデル夏の学校 2019

1-2-2. Soft-max MDP

GoalZiebart et al.(2009)は意思決定に関わる不確実性を考慮するためにMaxの代わりに期待最大効用（ログサム）を導入

𝑉𝑉 ≈ (𝑠𝑠)𝑅𝑅(𝑠𝑠, 𝑎𝑎)

このとき行動方策は以下のロジットの形で得られる次の状態の効用

確率論的なモデル(Recursive logitと同じ）次の状態にいるときの期待最大効用

将来を考慮した逐次選択のモデル

9

1-2-3. 最大エントロピー逆強化学習

行動モデル夏の学校 2019

逆強化学習は，< 𝑆𝑆, 𝐴𝐴, 𝑇𝑇 >とエキスパートの行動から報酬関数𝑅𝑅を推測する問題
報酬関数は以下の対数尤度関数が最大となるように決定する，すなわちargmax𝑅𝑅 � log 𝑃𝑃(𝜁𝜁𝑖𝑖 |𝑅𝑅)
𝑖𝑖

最大エントロピー原理より，観測された軌跡𝜁𝜁𝑛𝑛 が得られる確率は以下の式になる

exp{𝑅𝑅(𝜁𝜁𝑖𝑖 )}𝑃𝑃 𝜁𝜁𝑖𝑖 𝑅𝑅 =∑𝜁𝜁∈𝑍𝑍 exp{𝑅𝑅(𝜁𝜁)}

Ziebart et al.(2008)は，報酬関数に線形性を仮定すること，すなわち𝑅𝑅 ≔ � 𝜽𝜽T 𝐟𝐟𝐬𝐬𝑠𝑠∈𝜁𝜁𝑖𝑖により報酬関数が効率的に推定できることを示した．
ただし，𝜃𝜃はパラメータで，𝐟𝐟𝒔𝒔 は𝑘𝑘次元の特徴量ベクトル10

行動モデル夏の学校 2019

1-2-4. 最大エントロピー逆強化学習の推定
パラメータ𝜃𝜃は対数尤度関数の最大化によって決定される

exp{𝑅𝑅(𝜁𝜁𝑖𝑖 )}𝑃𝑃 𝜁𝜁𝑖𝑖 𝑅𝑅 =∑𝜁𝜁∈𝑍𝑍 exp{𝑅𝑅(𝜁𝜁)}

上の式の勾配を計算すると，以下の式が得られる

観測における特徴量の合計の平均値

モデルの特徴量の合計の期待値

𝔼𝔼𝑃𝑃𝜃𝜃 (𝜁𝜁) 𝒇𝒇𝑠𝑠 = � 𝑃𝑃(𝑠𝑠|𝜃𝜃)𝒇𝒇𝑠𝑠
𝑠𝑠

パラメータ𝜃𝜃のモデルで，状態𝑠𝑠を何回訪れるかState Visitation Frequency(SVF)と呼ばれる

11

行動モデル夏の学校 2019

1-2-4. 最大エントロピー逆強化学習の推定
𝑠𝑠1

モデルの特徴量の合計の期待値の計算

𝔼𝔼𝑃𝑃𝜃𝜃 (𝜁𝜁) 𝒇𝒇𝑠𝑠 = � 𝑃𝑃(𝑠𝑠|𝜃𝜃)𝒇𝒇𝑠𝑠
𝑠𝑠𝑠𝑠2
𝑠𝑠3

𝑠𝑠3

環境モデル

パラメータ𝜃𝜃のモデルで，状態𝑠𝑠を何回訪れるかState Visitation Frequency(SVF)と呼ばれる

右上の図のSVFの計算例

𝑠𝑠1

𝑠𝑠2

特徴量

t=1SVF

木

ベンチ

お店

0.3

t=2
0.1

1.4

100

0

0.7

0.6

1.3

11

00

00.3

0.3

00

1

t=0
１

𝑃𝑃(𝑠𝑠1 |𝑠𝑠1 )モデルから計算される存在確率特徴量（歩道）の期待値＝内積(SVF,木)＝2.7
特徴量（ベンチ）の期待値＝内積(SVF,ベンチ)＝1.3
特徴量（お店）の期待値＝内積(SVF,お店)＝0.3

12

行動モデル夏の学校 2019

1-2-4. 最大エントロピー逆強化学習の推定
価値関数の計算：価値反復法
数値反復計算によりBellman最適方程式を求める方法
強化学習分野では一般的な方法の一つ
Step1
全てのsにおいて価値関数𝑉𝑉 ≈ 𝑠𝑠 の値を初期化
これを𝑉𝑉0≈ 𝑠𝑠 と表す．
Step2全てのsに対して以下を実行する．

≈𝑄𝑄𝑘𝑘≈ 𝑠𝑠, 𝑎𝑎 = � 𝑃𝑃(𝑠𝑠|𝑠𝑠 ′ , 𝑎𝑎) 𝑅𝑅 𝑠𝑠𝑠, 𝑎𝑎 + 𝛾𝛾𝑉𝑉𝑘𝑘−1
𝑠𝑠𝑠𝑠𝑠′𝑉𝑉𝑘𝑘≈ 𝑠𝑠 = log � exp 𝑄𝑄𝑘𝑘≈ 𝑠𝑠, 𝑎𝑎

通常のMDPならMax𝑎𝑎

Step3≈𝑉𝑉𝑘𝑘≈ 𝑠𝑠 − 𝑉𝑉𝑘𝑘−1
𝑠𝑠 の値が十分小さくなるまでStep2を繰り返す

𝑃𝑃(𝑠𝑠|𝑠𝑠 ′ , 𝑎𝑎)：状態𝑠𝑠 ′ で行動𝑎𝑎を取ったときに状態𝑠𝑠に遷移する確率

13

行動モデル夏の学校 2019

1-2-4. 最大エントロピー逆強化学習の推定

パラメータ更新（勾配法）
𝜃𝜃’ ← 𝜃𝜃 − 𝛼𝛼𝛻𝛻𝜃𝜃 𝐿𝐿

報酬関数
R 𝑠𝑠 = 𝜽𝜽𝑻𝑻 𝒇𝒇
価値関数の計算（価値反復法）𝑉𝑉 ≈ 𝑠𝑠
方策

繰り返し計算

観測との差特徴量カウントの期待値

𝔼𝔼𝑃𝑃𝜃𝜃 (𝜁𝜁) 𝒇𝒇𝑠𝑠
14

1-2-5. Recursive logit modelとの違い

行動モデル夏の学校 2019

• 価値関数の求解方法
 Recursive logit modelは逆行列計算を利用- 逆行列計算ができるよう，負の効用のみを扱う
 最大エントロピー逆強化学習は反復計算を利用- 正の効用も扱うことができる
ただし，Discounted recursive logit (Oyama and Hato, 2017)ではrecursive logitに時間割引を導入し，反復計算によって価値関数を計算している．

15

行動モデル夏の学校 2019

1-3. 時空間制約下の方策
最大エントロピー逆強化学習で得られる方策：𝜋𝜋(𝑎𝑎𝑡𝑡 |𝑠𝑠𝑡𝑡 )状態間の遷移確率は方策と遷移モデルを用いて以下のように表せる．

𝑃𝑃 𝑠𝑠𝑡𝑡+1 |𝑠𝑠𝑡𝑡 = ∑𝑎𝑎 𝑃𝑃(𝑠𝑠𝑡𝑡+1 |𝑠𝑠𝑡𝑡 , 𝑎𝑎𝑡𝑡 )𝜋𝜋(𝑎𝑎𝑡𝑡 |𝑠𝑠𝑡𝑡 )遷移モデル

ここで，目的地に関する時空間制約を導入したい．
ある目的地に到着必須時刻𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 に到着するという制約の下に，時空間制約がないときの遷移確率と整合的な状態遷移確率𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 , 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 を求めたい．

16

行動モデル夏の学校 2019

1-3. 時空間制約下の方策
以下の後ろ向き確率を導入する．

方策に従って行動したときに目的地に時間通りに到着する確率（後向き確率）𝜷𝜷 𝒔𝒔𝒕𝒕 ≡ 𝑷𝑷 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 |𝒔𝒔𝒕𝒕 = �𝒔𝒔𝒕𝒕+𝟏𝟏=�

𝑷𝑷 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒗𝒗𝒗𝒗 |𝒔𝒔𝒕𝒕+𝟏𝟏 𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕

𝒔𝒔𝒕𝒕+𝟏𝟏𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 𝜷𝜷(𝒔𝒔𝒕𝒕+𝟏𝟏 )ベイズの定理より以下の関係が得られる．
𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 , 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 ===𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 , 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 |𝒔𝒔𝒕𝒕
𝑷𝑷 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 |𝒔𝒔𝒕𝒕

𝑷𝑷 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 |𝒔𝒔𝒕𝒕+𝟏𝟏 𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕
𝑷𝑷 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂 |𝒔𝒔𝒕𝒕

𝜷𝜷 𝒔𝒔𝒕𝒕+𝟏𝟏
𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕
𝜷𝜷 𝒔𝒔𝒕𝒕

後ろ向き確率の比を用いて計算できる

17

行動モデル夏の学校 2019

1-3. 時空間制約下の方策
簡易ネットワーク上の計算例

問題設定：𝑡𝑡 = 3に𝑠𝑠 (3) に到着

経路は3種類

1 3 2
14 4 383 2 1
1= � � =4 3 24𝑃𝑃 𝜁𝜁1 = � � =𝑃𝑃 𝜁𝜁33 1 2
4 3 3𝑃𝑃 𝜁𝜁2 = � � =16

それぞれの経路の選択確率は𝑃𝑃 𝜁𝜁1 |𝑍𝑍 = 3/13, 𝑃𝑃 𝜁𝜁2 |𝑍𝑍 = 4/13, 𝑃𝑃 𝜁𝜁3 |𝑍𝑍 = 6/13

前向き確率の計算（時空間制約なし）

18

行動モデル夏の学校 2019

1-3. 時空間制約下の方策
簡易ネットワーク上の計算例

𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 , 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂

𝜷𝜷(𝒔𝒔𝒕𝒕 )𝜷𝜷(𝒔𝒔𝒕𝒕+𝟏𝟏 )𝟑𝟑
𝟏𝟏⁄𝟐𝟐 𝟏𝟏
� =𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 ==𝜷𝜷(𝒔𝒔𝒕𝒕 )𝟏𝟏𝟏𝟏⁄𝟐𝟐𝟐𝟐 𝟒𝟒 𝟏𝟏𝟏𝟏

𝑷𝑷 𝒔𝒔𝒕𝒕+𝟏𝟏 |𝒔𝒔𝒕𝒕 , 𝒔𝒔𝑻𝑻𝒂𝒂𝒂𝒂𝒂𝒂𝒂𝒂

𝜷𝜷(𝒔𝒔𝒕𝒕+𝟏𝟏 )後向き確率の計算

時空間制約下についても後ろ向き確率を用いることで逐次選択，経路列挙必要なしで整合的な経路を与えることができる

Prism-based path set (Oyama and Hato, 2019) ？

19

行動モデル夏の学校 2019

２．EMアルゴリズム

20行動モデル夏の学校 20192-1. EMアルゴリズムとは？
EMアルゴリズムとは• 潜在変数を持つ確率モデルの最尤解を求めるための一般的な手法• 観測変数𝑿𝑿の尤度関数𝑃𝑃(𝑿𝑿|𝜽𝜽)を最大化するパラメータ𝜽𝜽を見つける際に，潜在変数𝒁𝒁を導入した以下の式の方が扱いやすい場合がある

𝑝𝑝 𝑿𝑿 θ = � 𝑝𝑝(𝑿𝑿, 𝒁𝒁|𝜃𝜃)
𝒁𝒁• 混合分布の推定や潜在クラスモデルの推定などによく用いられる複雑な分布を考えるより𝜋𝜋1 𝒩𝒩 𝜇𝜇1 , 𝜎𝜎1 + 𝜋𝜋2 𝒩𝒩 𝜇𝜇2 , 𝜎𝜎2正規分布の和で表現どちらの正規分布から得られたかが潜在変数

21

行動モデル夏の学校 20192-2. EMアルゴリズムの導出
ln 𝑝𝑝 𝑿𝑿 θ = ln � 𝑝𝑝(𝑿𝑿, 𝒁𝒁|𝜃𝜃)

最大化したい対数尤度𝒁𝒁

𝑝𝑝(𝑿𝑿, 𝒁𝒁|𝜃𝜃)= ln � 𝑞𝑞(𝒁𝒁)
𝑞𝑞(𝒁𝒁)𝒁𝒁𝑝𝑝(𝑿𝑿, 𝒁𝒁|𝜃𝜃)≥ � 𝑞𝑞 𝒁𝒁 ln𝑞𝑞(𝒁𝒁)
𝒁𝒁= ℒ(𝑞𝑞, 𝜃𝜃)EMアルゴリズムの理解には以下のページがおすすめ「EMアルゴリズム徹底解説」
https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb

Jensenの不等式

変分下限と呼ぶ

ln 𝑝𝑝 𝑿𝑿 θ − ℒ(𝑞𝑞, 𝜃𝜃)を計算するとln 𝑝𝑝 𝑿𝑿 θ − ℒ 𝑞𝑞, 𝜃𝜃 =（中略）= 𝐾𝐾𝐾𝐾[𝑞𝑞(𝒁𝒁)||𝑝𝑝(𝒁𝒁|𝑿𝑿, 𝜃𝜃)]KLダイバージェンス(2つの確率分布間の距離）

𝑞𝑞 𝒁𝒁 = 𝑝𝑝(𝒁𝒁|𝑿𝑿, 𝜃𝜃)を選ぶことでln 𝑝𝑝 𝑿𝑿 θ = ℒ(𝑞𝑞, 𝜃𝜃)になる

22

行動モデル夏の学校 20192-2. EMアルゴリズムの導出

EMアルゴリズムでは，ln 𝑝𝑝 𝑿𝑿 θ の代わりに変分下限ℒ(𝑞𝑞, 𝜃𝜃)の最大化を行う．
ℒ(𝑞𝑞, 𝜃𝜃)の引数𝑞𝑞と𝜃𝜃を交互に最適化しながら，ℒ(𝑞𝑞, 𝜃𝜃)の最大化を行う．
E-step: 𝑞𝑞の最適化

𝜃𝜃を固定して， ln 𝑝𝑝 𝑿𝑿 θ 最大となる𝑞𝑞を選ぶ．
すなわち，以下を求める．

𝑞𝑞 𝒁𝒁 = 𝑝𝑝(𝒁𝒁|𝑿𝑿, 𝜃𝜃)M-step: 𝜃𝜃の最適化

𝑞𝑞(𝒁𝒁)を固定して，θを最大化する．
その際に事後分布の算出に用いたθはθold として固定したまま数値計算ℒ 𝑞𝑞, 𝜃𝜃 = � 𝑞𝑞 𝒁𝒁|𝑿𝑿, 𝜃𝜃 old ln𝑝𝑝 𝑿𝑿, 𝒁𝒁 𝜃𝜃 + const
𝒁𝒁= 𝑄𝑄 𝜃𝜃, 𝜃𝜃 old + const

完全データの対数尤度関数

𝑄𝑄 𝜃𝜃, 𝜃𝜃 old を𝜃𝜃に関して最大化すれば良い．
完全データの対数尤度が最適化可能であれば𝑄𝑄 𝜃𝜃, 𝜃𝜃 old は最大化できる23

行動モデル夏の学校 20192-3. EMアルゴリズムのイメージ
E-stepE-stepM-stepln 𝑝𝑝 𝑿𝑿 θℒ(𝑞𝑞, 𝜃𝜃)𝜃𝜃𝜃𝜃ln 𝑝𝑝 𝑿𝑿 θ𝜃𝜃ℒ(𝑞𝑞, 𝜃𝜃)𝜃𝜃𝜃𝜃

M-stepln 𝑝𝑝 𝑿𝑿 θold𝜃𝜃

newℒ(𝑞𝑞, 𝜃𝜃)

最尤解に近づく𝜃𝜃

ln 𝑝𝑝 𝑿𝑿 θℒ(𝑞𝑞, 𝜃𝜃)𝜃𝜃𝜃𝜃 old 𝜃𝜃 new
E-stepとM-stepを収束するまで繰り返す

24

行動モデル夏の学校 2019

３．歩行者回遊行動への適用

25

行動モデル夏の学校 20193-1. 歩行者回遊行動への適用
適用例
入力：
行動系列：𝜻𝜻 = {𝜁𝜁1 , 𝜁𝜁2 , 𝜁𝜁3 }環境の特徴量：歩道，芝生，桜の木，POI
出力：
各特徴量の重み：𝜃𝜃 = [𝜃𝜃ww , 𝜃𝜃ch , 𝜃𝜃POI ]T
報酬関数：R 𝑠𝑠R 𝑠𝑠 = 𝜃𝜃ww 𝑓𝑓ww + 𝜃𝜃ch 𝑓𝑓ch + 𝜃𝜃POI 𝑓𝑓POIIllustrative example of theproblem (POI: point of interest)

報酬関数は各特徴量の線形和推定されたパラメータは，特徴の持つ価値（魅力度）を表す
歩行者の軌跡は効用最大化行動の結果とみなす26

3-1. 歩行者回遊行動への適用

行動モデル夏の学校 2019

正の効用を扱う場合，同じ場所に滞在することで効用をいくらでも大きくできてしまう⇒ 時空間制約が必要
しかしながら，実際の観測データからの推定を考えた場合，時間制約は観測できない
図の例のように最短経路の軌跡が観測された場合，Goal地点の報酬が高いのか，単に急いでいるか区別がつかない．
本来は到着時間の制約を持っているにも関わらず，制約がないとみなして推定を行うと，推定にバイアスが生じる

最短経路

27

行動モデル夏の学校 20193-2. 潜在変数導入の考え方
潜在的な到着時間制約の確率分布（ex. 負の二項分布）

時間制約：強い𝑷𝑷(𝝉𝝉)𝜏𝜏0

時間制約：弱い

この時間制約の条件付きの中でLogitによって経路が選択されている

この時間制約の条件付きの中でLogitによって経路が選択されている

𝝉𝝉𝑃𝑃 𝜁𝜁 = � 𝑃𝑃(𝜏𝜏)𝑃𝑃(𝜁𝜁|𝜏𝜏)

軌跡𝜁𝜁の得られる確率𝜏𝜏

到着時間制約が𝜏𝜏である確率
到着時間制約𝜏𝜏を満たす軌跡の集合の（何らかの確率分布に従うと仮定） 中からの軌跡𝜁𝜁の選択確率

基本的な考え方は，確率的な選択肢集合形成
Manski(1977)と同じ

28

行動モデル夏の学校 20193-3. 定式化

𝜁𝜁𝑛𝑛 = 𝑠𝑠0 , 𝑠𝑠1 , ⋯ , 𝑠𝑠𝑇𝑇

軌跡の観測数：N状態の系列観測された軌跡
到着時間制約（潜在変数）

到着時間制約に対応する𝜏𝜏のみ1になる2値変数

潜在変数の導入により，軌跡の得られる確率𝑝𝑝(𝜁𝜁)は以下のように表される

観測された軌跡集合𝜁𝜁の得られる確率

到着時間制約の確率分布
到着時間制約を満たす軌跡の集合の中からの軌跡𝜁𝜁の選択確率

29

行動モデル夏の学校 20193-3. 定式化
個々の到着時間制約の事前分布は以下の式で表される

負の二項分布(𝜏𝜏0 = 4, 𝜇𝜇 = 0.4)

ここで，到着時間制約が従う確率分布として，以下の負の二項分布を仮定
負の二項分布𝜇𝜇 ：パラメータ（試行の成功確率に相当）𝜏𝜏0 ：最短所要時間（ODから決まる）

また，到着時間制約の下の軌跡の選択は最大エントロピー原理より以下に従うと仮定

30

3-3. 定式化
したがって，完全データの対数尤度関数は行動モデル夏の学校 2019

𝑅𝑅 𝜁𝜁 = 𝜃𝜃 𝑇𝑇 𝐟𝐟𝜁𝜁 を仮定

完全データの対数尤度関数が最適化可能であればEMアルゴリズムの適用が可能

31

行動モデル夏の学校 20193-4. EMアルゴリズムによる推定
E-step𝑧𝑧𝑛𝑛𝜏𝜏 の𝒛𝒛の事後分布による平均を求める．

𝒛𝒛の事後分布はなので，は負担率(responsibility)

各到着時間制約τが経路ζの観測を説明する度合

32

行動モデル夏の学校 20193-4. EMアルゴリズムによる推定
M-stepE-stepで求めた事後分布を用いて対数尤度の期待値を求める

Q関数と呼ばれる

M-stepではQ関数を微分して，最尤解を求めれば良いµの最尤解は

𝜕𝜕𝜕𝜕= 0を解くことで求まる．
𝜕𝜕𝜕𝜕𝜃𝜃の最尤解は

𝜕𝜕𝜕𝜕= 0を解くことで求まる．
𝜕𝜕𝜃𝜃

33

行動モデル夏の学校 20193-4. EMアルゴリズムによる推定
到着時間制約分布のパラメータ：𝜇𝜇の最尤解
以下の式を解けば良い

したがって，𝜇𝜇について解くと= 𝑁𝑁

𝜇𝜇は最短所要時間に対する迂回の度合の逆数として求まる= 𝑁𝑁𝜏𝜏

=最短所要時間
到着時間制約の平均値=１
平均迂回係数34

行動モデル夏の学校 20193-4. EMアルゴリズムによる推定
到着時間制約分布のパラメータ：𝜃𝜃の最尤解
以下の式を解けば良い

したがって，観測の特徴量の平均値負担率による重みづけ

各到着時間制約τの下の特徴量の期待値（求め方については，付録を参照）

強い制約を受けている場合，経路の選択肢が少ないためパラメータの大小に関わらず，モデルと観測が一致しやすい⇒自由に回遊している（と思われる）データほどパラメータ推定に寄与する

35

3-4. EMアルゴリズムによる推定

行動モデル夏の学校 2019

Algorithm

1. 負担率γ(𝑧𝑧𝑛𝑛𝜏𝜏 )を設定

2. 収束するまで以下を繰り返す
1. パラメータの更新
1. µの更新

2. θの更新

2. 負担率の更新

36参考文献

行動モデル夏の学校 2019

関連発表
[1] 日高健，早川敬一郎，西智樹，薄井智貴，山本俊行，逆強化学習を用いた報酬関数推定と時空間制約下における歩行者の行動軌跡生成，
第58回土木計画学研究発表会，CD-ROM，2018.
[2] 日高健，早川敬一郎，西智樹，薄井智貴，山本俊行，歩行者の報酬関数と潜在的な到着時間制約を同時に推定する逆強化学習法，第59回
土木計画学研究発表会，CD-ROM，2019.
参考文献
Abbeel, P., Ng, A. Y., Apprenticeship learning via inverse reinforcement learning, In: Proceedings of the 21st international
conference on Machine learning (ICML), 2004.
Boularias, A., Kober, J., Peters, J. Relative entropy inverse reinforcement learning. In Proceedings of the 14th International
Conference on Artificial Intelligence and Statistics (AISTATS), 182-189, 2011.
Fosgerau, M., Frejinger, E., Karlstrom, A., A link based network route choice model with unrestricted choice set. Transp. Res. Part
B: Methodol. 56, 70–80, 2013.
Kitani, K. M., Ziebart, B. D., Bagnell, J. A., Hebert, M., Activity forecasting. In: European Conference on 418 Computer Vision(ECCV), 201–214, 2012.
Manski, C. F., The structure of random utility models. Theory and decision 8 (3), 229–254, 1977.
Ng, A. Y., Russell, S. J., Algorithms for inverse reinforcement learning. In: Proceedings of the 17th International Conference onMachine Learning (ICML), 663–670, 2000.
Oyama, Y., Hato, E., A discounted recursive logit model for dynamic gridlock network analysis. Transp. Res. Part C: Emerg.
Technol. 85, 509–527, 2017.
Oyama, Y., Hato, E., Prism-based path set restriction for solving Markovian traffic assignment problem. Transp. Res. Part B:Methodol. 122, 528-546, 2019.
Russell, S.: Learning agents for uncertain environments, In: Proceedings of the 11th annual conference on computational learningtheory, 101-103, 1998.
Wulfmeier, M., Ondruska, P., Posner, I. Maximum entropy deep inverse reinforcement learning. arXiv preprint arXiv:1507.04888,2015.
Ziebart, B. D., Maas, A. L., Bagnell, J. A., Dey, A. K., Maximum entropy inverse reinforcement learning. In: Proceedings of theAssociation for the Advancement of Artificial Intelingence (AAAI), 1433–1438, 2008.
Ziebart , B. D., Ratliff, N., Gallagher, G., Mertz, C., Peterson, K., Bagnell, J. A., Hebert, M., Dey, A. K., Srinivasa, S., Planningbased prediction for pedestrians. In IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 3931–3936,
2009.

37

行動モデル夏の学校 20193-5. 数値実験
実験による検証①StartGoal𝜏𝜏0 = 4𝜇𝜇 = 0.4

設定した到着時間制約分布

𝜃𝜃walkway , 𝜃𝜃cherry , 𝜃𝜃PoI = [2,2,4]𝜏𝜏̅ =𝜏𝜏0= 10𝜇𝜇

到着時間制約の平均10ステップ

経路データのみから到着時間制約分布のパラメータµと特徴量の重みθを推定できるか検証

38

行動モデル夏の学校 20193-5. 数値実験

・𝑁𝑁 = 1000
・計算する𝜏𝜏は30ステップまで（データの長さも全て30ステップ分）・負担率の初期値 𝑧𝑧𝑛𝑛,30 = 1

・10回程度の反復で収束・設定パラメータが再現できている・10回程度の反復で収束・設定パラメータが再現できている

時間制約を考慮しない推定値

39

3-5. 数値実験

到着時間制約の比較（推定値は𝜏𝜏𝑒𝑒𝑒𝑒𝑒𝑒 = argmax𝜏𝜏 𝛾𝛾(𝑧𝑧𝑛𝑛𝜏𝜏 )）行動モデル夏の学校 2019

𝛾𝛾(𝑧𝑧𝑛𝑛𝜏𝜏 )の例

40

行動モデル夏の学校 20193-5. 数値実験
実験による検証②StartGoal𝜏𝜏0 = 4𝜇𝜇 = 0.4

設定した到着時間制約分布

𝜃𝜃walkway , 𝜃𝜃cherry , 𝜃𝜃PoI = [2,2,4]𝜏𝜏̅ =𝜏𝜏0= 10𝜇𝜇

到着時間制約の平均
10ステップ

経路データのみから到着時間制約分布のパラメータµと特徴量の重みθを推定できるか検証

41

行動モデル夏の学校 20193-5. 数値実験
・𝑁𝑁 = 1000
・計算する𝜏𝜏は30ステップまで・負担率の初期値 𝑧𝑧𝑛𝑛,30 = 1

・3回程度の反復で収束・設定パラメータが再現できている

・3回程度の反復で収束・設定パラメータが再現できている

時間制約を考慮しない推定値

423-5. 数値実験

到着時間制約の比較（推定値は𝜏𝜏𝑒𝑒𝑒𝑒𝑒𝑒 = argmax𝜏𝜏 𝛾𝛾(𝑧𝑧𝑛𝑛𝜏𝜏 )）行動モデル夏の学校 2019

𝛾𝛾(𝑧𝑧𝑛𝑛𝜏𝜏 )の例

43

3-5. 数値実験

行動モデル夏の学校 2019

どちらの実験でもバイアスのない推定ができた

44

行動モデル夏の学校 2019

（付録）到着時間制約下の特徴量期待値の求め方

時間制約𝜏𝜏の下での特徴量期待値

時間制約𝜏𝜏の下でのSVF

到着時間制約

forward probability

𝑆𝑆𝑆𝑆𝑆𝑆𝜏𝜏 = � 𝑃𝑃(𝑠𝑠𝑡𝑡 |𝑠𝑠0 , 𝑠𝑠𝜏𝜏 )𝑡𝑡

𝑃𝑃 𝑠𝑠𝑡𝑡 𝑠𝑠0 , 𝑠𝑠𝜏𝜏 =𝑃𝑃 𝑠𝑠𝜏𝜏 𝑠𝑠0 , 𝑠𝑠𝑡𝑡 𝑃𝑃(𝑠𝑠0 , 𝑠𝑠𝑡𝑡 )𝑃𝑃(𝑠𝑠0 , 𝑠𝑠𝜏𝜏 )𝑃𝑃 𝑠𝑠𝜏𝜏 𝑠𝑠𝑡𝑡 𝑃𝑃 𝑠𝑠𝑡𝑡 𝑠𝑠0 𝑃𝑃(𝑠𝑠0 )=𝑃𝑃(𝑠𝑠0 , 𝑠𝑠𝜏𝜏 )==backward probability

𝑃𝑃 𝑠𝑠𝜏𝜏 𝑠𝑠𝑡𝑡 𝑃𝑃 𝑠𝑠𝑡𝑡 𝑠𝑠0
𝑃𝑃(𝑠𝑠𝜏𝜏 |𝑠𝑠0 )𝛼𝛼(𝑠𝑠𝑡𝑡 )𝛽𝛽(𝑠𝑠𝑡𝑡 )𝛽𝛽(𝑠𝑠0 )45

（付録）到着時間制約下の特徴量期待値の求め方

行動モデル夏の学校 2019

まとめると
時間制約𝜏𝜏の下でのSVF
𝛼𝛼(𝑠𝑠𝑡𝑡 )𝛽𝛽(𝑠𝑠𝑡𝑡 )𝑆𝑆𝑆𝑆𝑆𝑆𝜏𝜏 = �𝛽𝛽(𝑠𝑠0 )𝑡𝑡

時間制約のない一般のSVF
𝑆𝑆𝑆𝑆𝑆𝑆 = � 𝛼𝛼(𝑠𝑠𝑡𝑡 )𝑡𝑡

46