パラメータ推定の理論と実践

山梨大学
佐々木邦明
BEhavior Study for Transportation Graduate school, Univ. of Yamanashi

モデルのパラメータ(母数)とは• 離散型選択モデルの母数

ランダム効用モデル– 確率分布を特徴付ける数– 例えば，期待値，分散等期待値と分散が母数– 変動要因を規定する少数の要素に分解する
V1 = b x

U ~ N (V1 , s )• パラメータ推定– 観察されたデータに基づいて母数を統計的に推定すること

V1 = b x2

推定アルゴリズム

3最尤推定法Ln (q | x ) = Õ f (xi | q )n• 点推定量を求める最もポピュラーな方法

i =1• 右上の式をθの関数とみなしたものが尤度関数

データ（a，b）が得られたとき，全体の平均がいくつとするのがよいか＝平均がいくつだったら（a，b）が得られやすいか？

• 尤度関数を最大化するθの値を最尤推定量とするのが最尤推定法

選択モデルの場合，f が選択確率個人の選択確率を全員で掛け合わせる

MaxLikelihood = max Õ Pi (xi | q )q ÎQi4

確率モデルによる個人の選択記述
U = U1 - U 2• 確率的選択モデルは効用がある分布に従うと仮定• ある選択肢の選択確率は確率分布を積分する• ロジットモデルやNL等のGEVモデルはこの積分した結果が積分のない形式で表現可能２項モデル

P(1)

０

MaxLikelihood = max Õ ò f (xi | q )q ÎQiÕ ò f (x | q ) をどうやって最大化するか？
ii5

最尤推定のアルゴリズムの考え方•

最大値をいきなり求めるのは無理な場合•対数尤度関数の段階的な最大化
1. 初期値を与える
2. 初期値周りで勾配(1次微分）等を用いて次の推定値の方向を決める
3. 初期値付近のステップサイズを１次微分，２次微分を用いて適切に決めて次の推定値を決める
4. 収束基準(尤度関数の一時微分ベクトル)を判定し，収束していない場合は，現在の値を初期値として2に進む

Ln (q | x )

qˆML
q2
q ST
qq1

ロジットモデルのみ全域で単峰性あり

6代表的な繰り返し計算法
尤度関数を最大化:尤度関数の一階微分＝０を解く

g(q)• Newton-Raphson法– テイラー展開の１次近似を利用して進める– 解の収束が早い(ステップ数が少ない)– ヘッセ行列を計算するので時間がかかる• 準Newton法（BFGS法）– ヘッセ行列を逐次近似する．
– 二階微分をする必要がないので計算が早い

q2

qq1

q0

q n +1 = q n - H -1 g1
H:尤度関数の二階微分
ヘッセ行列
g: 尤度関数の一階微分7パラメータ推定がうまくいかない• 収束するとはθn+1とθnが同じになる– g1が0になる• 収束しない– 無限に繰り返す– θ2が計算不能• 局所最適解– 見かけ上の最大化

H-1が存在しない(計算できない）
変数が完全相関
変数が効用関数に影響していない関数の近似状況
2次関数近似初期値の問題そもそも推定不可能
推定プログラムに誤り

gg(q()q)

q1q1

qq
q0q0
q28

そもそも推定不能• 最大値において唯一解が求まらない可能性がある（最大値となるパラメータベクトルが無数にある）→Identification Problem
例
U1

U3

U2
効用に適当な数を足しても選択確率に変化はない
定数項は相対的数値
09

ちょっと考える• 例：3カ所の買物目的地選択
A:駐車場有料，B：駐車場有料，C:駐車場無料
このとき，駐車場無料は大きな選択要因
ならば駐車場無料ダミー変数を選択肢に導入

U1i = a1 ´1 + b1 ´ 0 + K + e 1i
U 2i = a 2 ´1 + b1 ´ 0 + K + e 2i
U 3i = a 3 ´ 0 + b1 ´1 + K + e 3i
これは区別できない
10

駐車場の台数や料金無料は効果が計測不能？
• 例えば自動車保有ダミーddを使ってみる
U1i = (a1 + d di ´ b1 ´100 )´1 + K + e 1i

U 2i = (a 2 + d di ´ b1 ´ 200 )´ 1 + K + e 2i

U 3i = (d di ´ b1 ´1000 )´1 + K + e 3i

11

積分（INTEGRATION）

12

数値積分
•

個人ごとの確率を計算するためには定積分することが必要だが，閉じた形で表現できないものは数値積分を行う必要がある．
9•

数値積分：関数を多項式近似などを用いて積分可能な形にして，ある間隔の関数の和に直して積分を行う方法．
– Gaussの求積法：近似誤差を最も小さくなるように区間と幅を定める方法．
– Gauss-Legendre法：多項式近似ができる関数においては，多項式の次元にあわせた分割を行うと非常に精度よく積分を行える
方法

87654123

複雑な関数を簡単な１次関数で近似した例．

13

シミュレーション(乱数)による積分• 面積を求める別の方法を
考える

987

• 2次元乱数を[３×９]の範囲に散らばらせて，それが積分したい関数より上か(青丸)下か(赤丸) を判定する．
• ドットをうまく大量に発生させれば精度よく積分を近似できる．

654123

この図だと点の合計は20．そのうち赤丸１１個，青丸９個なのでこれより関数の積分結果は14.85となる

14

乱数の生成• 何を積分したいのか？
– 乱数の発生ルーチンを条件に合わせて改良• 制約つきの乱数• 多変量正規分布の乱数• 乱数と同等の性質を持ち，同じ値が出ず，必要な全域がカバーされる数値の生成法があるとより効率的な積分が可能
15

制約つきの乱数• 例えば平均０で分散1の正規分布の-3から1までの間の乱数．

10.80.6• 対応する分布関数の間の乱数を発生させてそれを密度に対応させる

0.40.2• この変換によって，密度の高いところはたくさん，密度の薄いところは少しという密度関数に応じた乱数になる．

0.6
0.4
0.2-4-22416

多変量正規分布の発生• 例えば標準正規分布μを平均b標準偏差sに変換するためには e = b + sm と線形変換をすればいい• 同じ事を多変量正規分布に応用する．標準正規分
布の乱数ベクトルηを e = b + Lh と変換する．ただしLはLL’=Wとなる行列．するとç はMVN(b,W)になる．
例えば，2次元の例を示すと

æ e 1 ö æ b1 ö æ l11 0 öæ h1 ö÷÷çç ÷÷çç ÷÷ = çç ÷÷ + ççè e 2 ø è b2 ø è l12 l22 øèh 2 ø
e 1 = b1 + l11h1

e 2 = b2 + l12h1 + l22h 217

準乱数• 高次積分において乱数による積分の精度を上げるためには，計算速度が多大に犠牲になる321• 対処方法として系統的な数列（準乱数）がある

イメージ12318

準乱数の例o 代表的例としてHalton数列がある．その計算方法は素数pに対して
ttt{st +1 = st , st + 1 p , st + 2 p L, st + ( p - 1) p}

例えばp=3ならば，初期値0として１/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9・・・• 多次元化– 数列の異なる素数pを決めて，それぞれに応じて数列を作り多次元化する．

• 正規分布化– 数列を制約つきの乱数発生と同様の変換をして正規分布化
19

シミュレーションによる尤度計算
•手順
1.
2.
3.
4.
5.

誤差項の密度関数から選択肢の数の次元の(準)乱数を発生させる
この乱数を誤差の値として，各代替案の効用値を計算（積分）する
代替案iの効用値とその他の代替案の効用との値を比較し，それらの大小関係を1-0の変数Gで記述する．
1～3のステップを繰り返す．その反復回数をRとする．
1 R r
P=G
シミュレーションされた確率は i R år =1となり，この値は不偏推定量である．

効用を確定値にする

確定的に選択を決定

比率を確率に置き換える

これを尤度として最大になるようにパラメータをアップデートする
20

連続型シミュレーション
1. 誤差項の密度関数から選択肢の数の次元の乱数を発生させる
2. この乱数を誤差の値として，各代替案の効用値を計算する
3. この効用をロジット変換を行い，連続的な0-1変数に変換する．この変換後の値をSとする
4. 1～3のステップを繰り返す．その反復回数をRとする．
5. シミュレーションされた確率はは不偏推定量である．

1 R r
Pi = å S
R r =1となり，この値21シミュレーションベースのパラメータ推定法• シミュレーション尤度最大化（MSL）– シミュレーションによって計算された確率を尤度として，最大化を行う．

• 特性– サンプル数と乱数発生回数に依存する．
– 乱数発生回数が十分大きいと一致性や漸近的有効性を持ち解析積分と一緒の特性を持つ．
– 乱数発生回数がサンプル数に対して小さく固定されると一致性もない．

22

まとめ23

統計的な視点• 推定尤度に基づくモデル構造の選択– 投入する情報量(データ)が一緒なら尤度関数の大小で判定– 情報量の異なるモデル間ではAIC等の統計量による比較• 推定パラメータに基づく変数の選択– ｔ値による取捨選択– 推定値の符号による判定• 費用と所要時間のパラメータは負24モデル構造はどうやって選ぶ？
• 行動に関するアプリオリな知識による（仮説）– 手段選択には費用は影響する• 線形効用関数に費用を導入する．でもタクシーは選択肢に入らない• 費用は非補償型のモデル構造とする

検証– バスとパークアンドライドの選択において，誤差項に相関があるだろう• NLやCNLの採用を考える• 推定結果からの後天的な知識による（仮説の修正）– モデルの推定結果が仮説を支持しないときに仮説が間違っている可能性を考える• 線形効用関数で費用のパラメータが負にならない→キャプティブな選択があるのでは？
• ネスティッドロジットモデルのログサムの係数が0～1にならない→ネスト構造が違っている？

検証25