r_nsdのブログ

r_nsdのブログ

勉強したこと・調べたこと・思ったことを残しておくためのブログ

MENU

離散選択モデリング入門（多項ロジットモデルの理論とPythonによる実装）

離散選択モデル 統計モデリング

はじめに

日常には選択行動が溢れています．例えば，移動するときにはバスや鉄道，タクシーなどの選択肢から一つの交通手段を選択します．このような離散的な選択肢の中から選択をする行動をモデル化することを「離散選択モデリング」と呼びます．

本記事では，離散選択モデリングの研究の流れと，最も基礎的なモデルである多項ロジットモデルについて書きました．またPythonによる実装コードも載せています．

離散選択モデリングの研究の流れ

ランダム効用理論

人々は選択するとき，効用が最大となる選択肢を選択すると仮定します．効用とは得られる満足感のようなものであり，交通手段選択の場合は料金，移動時間などの要素からなります．人々はたとえ同じ状況であっても，同じ選択を繰り返すとは限りません．つまり選択はランダム性がある行動です．
そこで効用に確率項を入れてランダム性を考慮して選択行動を扱えるようにしたのが「ランダム効用理論」です．

ランダム効用理論が提唱されたのは，Thurstone の食べ物の嗜好の研究 (1920年代) が初めてだと言われています．しかし，彼は確率項に正規分布を使用していたため，選択確率を簡単に計算できる形に定式化できず，30年以上この研究は衰退していました．

離散選択モデリングの発展

しかし，1960年代〜1980年代の研究で，確率項にガンベル分布を仮定することで選択確率を計算できるようになりました．その結果，離散選択モデリングの研究は発展していきました．一番の功労者は，経済学者であるDaniel Little McFadden（2000年にノーベル経済学賞を受賞）です．彼の研究 [McFadden, 1973] が現在の離散選択モデル（ランダム効用モデル）の基礎となっています．

その後，静的な離散選択モデルを動的に拡張し，一回の選択ではなく，選択の系列をモデル化する「動的離散選択モデル」の研究が1980年代から発展しました．動的離散選択モデルの研究としてベースになっているのはRust の研究 [Rust, 1987a,b] です．

これらの離散選択モデルは，経済学だけでなく，交通工学の分野でも発展していますし，情報科学の分野では逆強化学習という手法として発展しています．

近年では，各分野で提案されたモデルが実は同じものであることがわかったり*1，他分野の手法を取り入れることでモデルの性能を上げたりしながら研究が進んでいます．

多項ロジットモデル

ここからは，離散選択モデルのなかでも最も基礎的なモデルである「多項ロジットモデル」について説明します．

多項ロジットモデルの概要

離散選択モデルでは，効用関数を確定項（決定項）と確率項（誤差項）の和として定義します． 効用は選択に関わる要素から算出されますが，どの要素をどれくらい重要視するかは異なるため，確定項を重みと要素の値の線形和として表します．確率項は，以下の点を包含しています．

  • 選択行動のランダム性
  • 確定項を重みと要素の線形和としたことによる誤差
  • 確定項に含まれなかった要素
  • 変数の測定誤差

これより，ある個人nのある選択肢iに対する効用関数は，次式のようになります．

\begin{aligned} U_{in} &= V_{in} + \varepsilon_{in} \\ &= \sum_k \beta_k x_{ikn}+\varepsilon_{in} \end{aligned}

ここで，V_{in}が確定項，\varepsilon_{in}が確率項，\betaが重み，xが要素の値，kが要素の数です．

確率項にガンベル分布を仮定したものを「ロジットモデル」と呼びます．（正規分布を仮定したものを「プロビットモデル」と呼びます．）

ガンベル分布

ガンベル分布とは連続確率分布の一つで，極値分布のタイプI型に相当します．母数は\mu \in \mathbb{R}, \eta>0です．

平均：E(x)=\mu+\gamma\eta （\gamma=0.577\cdots, オイラー定数）
分散：V(x)=\frac{\pi^ 2 \eta^ 2}{6}
最頻値：\mu

確率密度関数

\displaystyle{ f(x) = \frac{1}{\eta} \left\{ -\left( \frac{x-\mu}{\eta} \right) \right\} \exp \left[ -\exp \left\{ -\left( \frac{x-\mu}{\eta} \right) \right\} \right] }

累積分布関数

\displaystyle{ F(x)=\exp\left[-\exp \left\{-\left(\frac{x-\mu}{\eta} \right) \right\} \right] }


# ガンベル分布の確率密度関数
x = np.arange(-5,20,0.1)
pdf1 = stats.gumbel_r.pdf(x, loc=0.5, scale=2.0)
pdf2 = stats.gumbel_r.pdf(x, loc=1.0, scale=2.0)
pdf3 = stats.gumbel_r.pdf(x, loc=1.5, scale=3.0)
pdf4 = stats.gumbel_r.pdf(x, loc=3.0, scale=4.0)
    
plt.plot(x, pdf1, label=(r'$\mu=0.5, \eta=2.0$'))
plt.plot(x, pdf2, label=(r'$\mu=1.0, \eta=2.0$'))
plt.plot(x, pdf3, label=(r'$\mu=1.5, \eta=3.0$'))
plt.plot(x, pdf4, label=(r'$\mu=3.0, \eta=4.0$'))
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Probability Density Function")
plt.legend()

f:id:ryonsd:20191216135019p:plain

# ガンベル分布の累積分布関数
x = np.arange(-5,20,0.1)
cdf1 = stats.gumbel_r.cdf(x, loc=0.5, scale=2.0)
cdf2 = stats.gumbel_r.cdf(x, loc=1.0, scale=2.0)
cdf3 = stats.gumbel_r.cdf(x, loc=1.5, scale=3.0)
cdf4 = stats.gumbel_r.cdf(x, loc=3.0, scale=4.0)
    
plt.plot(x, cdf1, label=(r'$\mu=0.5, \eta=2.0$'))
plt.plot(x, cdf2, label=(r'$\mu=1.0, \eta=2.0$'))
plt.plot(x, cdf3, label=(r'$\mu=1.5, \eta=3.0$'))
plt.plot(x, cdf4, label=(r'$\mu=3.0, \eta=4.0$'))
plt.xlabel("x")
plt.ylabel("F(x)")
plt.title("Cumulative Distribution Function")
plt.legend()

f:id:ryonsd:20191216135109p:plain

ガンベル分布には以下の性質があります．

性質①
x_1, x_2がそれぞれガンベル分布(\mu_1, \eta), (\mu_2, \eta)に従うとき，x=x_1-x_2の累積分布関数はロジスティック分布の形になる．

\displaystyle{ F(x)=\frac{1}{1+\exp\{\eta(\mu_1-\mu_2-x)\}} }

性質②
x_1, x_2, \cdots, x_{\tau}が，それぞれガンベル分布(\mu_1, \eta), (\mu_2, \eta), \cdots, (\mu_{\tau}, \eta)に従うとき，\max(x_1, x_2, \cdots, x_{\tau})もガンベル分布に従い，その母数は次のようになる．

\displaystyle{ \left(\frac{1}{\eta}\ln\sum_{j=1}^{J}\exp\left(\eta \mu_j\right),\,\, \eta \right) }

選択確率の導出

ある選択肢iが選択される確率は，その選択肢iの効用関数値が他の選択肢の効用関数値よりも高くなる確率であり，この選択確率をもとに尤度を最大化する重みを推定することでモデルを構築します．

\begin{aligned} P_n (i) &= Pr[U_{in} > U_{jn}\,\,\,\, ,{\rm for} \forall j, j \neq i ] \\ &= Pr[U_{in} > \max (U_{jn}, \,\,{\rm for} \forall j, j \neq i)] \\ &= Pr[V_{in}+\varepsilon_{in} > \max (V_{jn}+\varepsilon_{jn}, \,\,{\rm for} \forall j, j \neq i)] \end{aligned}


ここで，\epsilon_{jn}はガンベル分布に従い，V_{jn}は定数であるため，{V_{jn}+\varepsilon_{jn}}はガンベル分布に従います．つまり，\varepsilon_{jn}が(\mu=0, \eta)のガンベル分布に従うと仮定すると，{V_{jn}+\varepsilon_{jn}}は(\mu=V_{jn}, \eta)のガンベル分布に従います．

\displaystyle{ \max (V_{jn}+\epsilon_{jn}, \,\,{\rm for} \forall j, j \neq i) = V_n^* + \epsilon_n^* }

とすると，ガンベル分布の性質②より，

\begin{aligned} V_n^* &= \frac{1}{\eta}\ln\sum_{j=1}^{J}\exp\left(\eta V_{jn}\right) \\ \varepsilon_n^* &\sim gumbel \left(0, \eta \right) \end{aligned}


これより，ガンベル分布の性質①を使うことで，次のように式変形でき選択確率が算出できます．

\begin{aligned} P_n(i) &= Pr[V_{in}+\varepsilon_{in} > V_n^* + \varepsilon_n^*] \\ &= Pr[\epsilon_n^*-\varepsilon_{in} \leq V_{in}-V_n^*] \\ &= \frac{1}{1+\exp\left\{\eta (0-0-(V_{in}-V_n^*)) \right\}} \\ &= \frac{1}{1+\exp\left\{\eta \left(-V_{in}+\frac{1}{\eta} \ln \sum_{j \neq i}\exp \left(\eta V_{jn}\right) \right) \right\}} \\ &= \frac{1}{1+\exp\left\{-\eta V_{in}+\ln \sum_{j \neq i}\exp \left(\eta V_{jn}\right) \right\}} \\ &= \frac{\exp(\eta V_{in})}{\sum_j \exp(\eta V_{jn})} \\ &= \frac{\exp(V_{in})}{\sum_j \exp(V_{jn})}\,\,\,\,(\eta=1) \end{aligned}


選択確率はソフトマックスの形になり，\etaが温度パラメータになっていますが，一般的には\eta=1に固定します．

重みの推定

全ての個人の選択に対して尤度を最大化する重み（パラメータ）\betaを推定します（最尤推定）．尤度および対数尤度は次のようになります．ここで，y_{in}は，ある個人nが選択肢iを選択したときに1を，それ以外で0をとります．

尤度

\displaystyle{ L(\beta) = \prod_n \prod_i P_n(i)^{y_{in}} }

対数尤度

\begin{aligned} LL(\beta) &= \sum_n \sum_i y_{in}\log P_n(i) \\ &= \sum_n \sum_i y_{in}\log \frac{\exp(V_{in})}{\sum_j \exp(V_{jn})} \\ &= \sum_n \sum_i y_{in}\log\frac{\exp(\sum_k \beta_k x_{ikn})}{\sum_j\exp(\sum_k \beta_k x_{jkn})} \\ &= \sum_n \sum_i y_{in}\left\{\log \exp\left(\sum_k \beta_k x_{ikn}\right)-\log\sum_j \exp\left(\sum_k \beta_k x_{jkn}\right)\right\} \\ &= \sum_n \sum_i y_{in} \left\{\sum_k \beta_k x_{ikn}-\log\sum_j \exp\left(\sum_k \beta_k x_{ikn}\right) \right\} \end{aligned}


この対数尤度LL(\beta)を最大化する\beta_{k}を求めます．つまり\frac{\partial LL(\beta_k)}{\partial \beta_k}=0となる\beta_{k}を求めます．ニュートン法などの最適化手法を用いることで解が求まります．

モデルの評価

評価指標

ロジットモデルの評価指標に関しては，こちらの記事がとても参考になるので，リンクを貼らさせていただきます．

[離散選択] ロジットモデルの決定係数 - ill-identified diary

一般的にモデルの評価には，尤度や情報量基準が使われますが，これらは複数のモデルのうちどれが真のモデルに近いかの比較に用いられます．これらの指標だとモデルがどれくらい選択行動を表現できるかがわかりません．
そこで，ロジットモデルの評価には，回帰分析に使われる決定係数のような指標として，McFadden の疑似決定係数が用いられます [McFadden, 1973]．これは対数尤度から求まる値で，0~1に収まります．1に近いほうが良いモデルです．

\displaystyle{ \rho^ 2 = 1-\frac{\ln L(\beta)}{\ln L(0)} }

ここで，L(\beta)は推定した重みでのモデルの尤度，L(0)は重みを0にしたときのモデルの尤度です． また，この指標ですと，説明変数xの数を増やせばL(\beta)の値が小さくなり，指標の値が大きくなるため，推定する重み（パラメータ）の数hに応じて調整を入れた自由度調整済み疑似決定係数という指標も用いられます．

\displaystyle{ \bar{\rho}^ 2 = 1-\frac{\ln L(\beta)-h}{\ln L(0)} }

評価指標の解釈

これらの評価指標の値をどう解釈するかについては，McFaddenの論文中で触れられています． [McFadden, 1975] では，決定係数R^ 2とMcFadden の疑似決定係数\rho^ 2の関係を示したグラフが示されています．

f:id:ryonsd:20191216213852p:plain

また，[McFadden, 1977] のp.35の文章中では，以下のよう書かれています．

Those unfamiliar with the \rho^ 2 index should be forewarned that its values tend to be considerably lower than those of the R^ 2 index and should not be judged by the standards for a "good fit" in ordinary regression analysis. For example, values of .2 to .4 for \rho^ 2 represent an excellent fit.

これらより，ロジットモデルの評価指標であるMcFadden の疑似決定係数は，回帰分析に使われるような決定係数R^ 2よりも低い値をとり，0.2~0.4の値を取れば良いモデルと言えるようです．

選択行動の理解と予測

モデルの構築により，各要素xに対する重み\betaを求めました．この重み\betaから選択行動に対する考察を行います．

まずは，各要素が選択に対して影響を与えているのかどうかを確認します．この確認にはt検定を用います．

t検定

t検定は「ある一つの変数とある特定の値や，ある二つの変数の間で，平均値に差があるかどうか」の検定に使われます． 具体的には，

  • 二つの母集団がどちらも正規分布に従うと仮定したとき，平均が等しいかどうか (Welch のt検定と呼ばれる）
  • 正規分布に従う母集団の平均が，特定の値に等しいかどうか
  • 回帰直線の勾配が0と有意に異なるかどうか

の検定に使われます．

ここでは，3つ目の例のように，「推定した重み\betaが0と有意に異なるかどうか」をt検定で調べます．

重みが0と異ならない場合は，その要素は選択に影響を与えていないということになります． 逆に，t検定により重みが0と有意に異なることが証明できれば，その要素が選択に影響を与えていたということがわかります．


t検定により要素が選択に影響を与えていたことが分かると，次は各要素の重みの正負から，要素が選択に対してポジティブorネガティブどちらの影響を与えているかがわかります．

また，重みの比をとることで，ユーザーがどの要素に対してどれくらいの価値を持っているかを示す限界支払意思額 (Marginal willingness to pay: MWTP) を調べることもできます．

さらに，モデルの構築により，未知の条件での選択確率を計算することができます．つまり将来の選択行動の予測が行えます．要素の値を変えたときの選択確率を調べることを感度分析と呼びます．

Python による多項ロジットモデルの実装

実装するにあたりRの離散選択モデリングのパッケージであるmlogit とこちらのスライドを参考にさせていただきました．

使用したデータおよびコードはこちらにあげています． github.com

使用データ

今回使用するデータは，Rのmlogit で用意されている交通手段選択のデータ ("Mode" という名前のデータ) です． 車，カープール，バス，鉄道の四つの交通手段から，コスト（金銭的な？）と時間をもとに一つ選択します．

こちらはwide型と呼ばれるデータです．一行が一個人を表しています．一番上の人は，各交通手段のコストと時間をもとにcarを選択したということです．

f:id:ryonsd:20191225172434p:plain
wide型のデータ

多項ロジットモデルを適用するに当たってデータの形式をwide型からlong型へと変換します．この変換はRのmlogit で用意されているメソッドを使って行うことができます．（ここの部分は実装していません...）

f:id:ryonsd:20191225172824p:plain
long型のデータ

こちらがlong型のデータです．一行が一選択肢を表しています．1~4行目までが一個人のデータです．

モデルの構築および評価

上記の説明の対数尤度を定義し，最適化を行い重みを推定します．最適化の部分はSciPyを用いました．SciPy には関数最小化のためのscipy.optimize.minimize というメソッドが用意されています．今回はL-BFGS-Bを用いて最適化しました．これはヘッセ行列を近似的に推定しながらニュートン法を行う準ニュートン法の一つです．

f:id:ryonsd:20191225173730p:plain

こちらが重みを推定した結果になります．McFaddenの疑似決定係数が0.348となりモデルがデータにフィッティングしていることがわかります．重みのt検定の結果を見ると，どの重みも0とは異なる値になっていることがわかります．

cost に対する重みもtime に対する重みも負の値を取っていることから，コスト・時間がかかればかかるほど，その交通手段は選択しなくなることがわかります．また定数の値はバスが0，自家用車が3.29，カープールが-0.90，鉄道が0.627となっており，自家用車にポジティブな印象を持っているのに対して，カープールには少しネガティブな印象を持ってそうなことがわかります．

まとめ

このような流れで多項ロジットモデルの構築と評価と考察を行います．
モデルを構築することで，重みの考察や要素の値を感度分析などができ，モデル化した選択行動の理解と予測が可能になります．

今回は離散選択モデルの中でも基礎的な多項ロジットモデルについて書きました．他のモデルとしては，ネスティッドロジットモデルや混合ロジットモデル，潜在クラスロジットモデルなどがあります．これらのサーベイ論文としては [Train, 2009] があります．

参考文献

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