23rd Sep 2022行動モデル夏の学校MNL, NL, RLコードの説明小林里瑳Multinomial Logit Model
### Create the data array from the existing dataset Data <- read.csv("C:/R/Data_Clean_English.csv",header=TRUE) ## Count the number of data rows hh <- nrow(Data) ## Set the initial values of the parameters (the number in parenthesis corresponds to the number of parameters to be estimated) b0 <- numeric(5) 1.まずはデータを読み込む2.データ数を数える3.パラメータ数の初期値(=0)を代入したベクトルを作る今回はパラメータ数を5にするので、要素が5の0ベクトルを作る2
#### Define the log-likelihood function of the logit model#### fr <- function(x) { ### declare the parameters### ## Alternative specific constants b1 <- x[1] b2 <- x[2] b3 <- x[3] b4 <- x[4] ## Travel time to destination d1 <- x[5] ## declare the log-likelihood variable, set value to 0 LL = 0 ### For this choice problem, we consider the following 5 modes??? ##鉄道(train) ## バス(bus) ##自動車(car) ##自転車(bike) ## 徒歩(walk) ility function: :introduce the desired explanatory variables in the fu年time # constant it fare <- Data$ModeAvailableTrain*exp(d1*Data$TotalTimeTrain/100 +b1*matrix(1,nrow =hh,ncol=1)) ous +b2*matrix(1,nrow =hh,ncol=1)) <- Data$ModeAvailableBus *exp(d1*Data$TotalTimeBus/100 <- Data$ModeAvailableCar  *exp(d1*Data$TimeCar/100 +b3*matrix(1,nrow =hh,ncol=1)) car bike <- Data$ModeAvailableBike *exp(d1*Data$TimeBike/100 +b4*matrix(1,nrow =hh,ncol=1)) <- Data$ModeAvailableWalk *exp(d1*Data$TimeWalk/100 walkベクトルxを入れると尤度を出す関数frを作る今回の効用関数Vtrain = b1 + d1 *鉄道所要時間/100 Vbus = b2 + d1 *バス所要時間/100 Vcar = b3 + d1 *自動車所要時間/100 Vbike = b4 + d1 *自転車所要時間/100 +d1 *徒歩所要時間/100 Vwalk = MNLの確率exp( μV ; ) P(i) = Σ j ∈ A exp(μVi) W
### Calculate the choice probabilities ## calculate the Inclusive Value (the denominator of the choice probabilities equation) deno <- (car + train + bus + bike + walk) ## Calculate indiviudal choice probabilities Ptrain <- Data$ModeAvailableTrain*(train / deno) <- Data$ModeAvailableTrain  *(bus Pbus / deno) <- Data$ModeAvailableCar  *(car / deno) Pcar Pbike <- Data$ModeAvailableBike *(bike  / deno) <- Data$ModeAvailableWalk *(walk  / deno) Pwalk ## Avoid problems stemming from choice probabilities becoming zero. Ptrain <- (Ptrain!=0)*Ptrain + (Ptrain==0) + (Pbus==0) <- (Pbus!=0)*Pbus Pbus <- (Pcar  !=0)*Pcar   + (Pcar  ==0) Pcar Pbike <- (Pbike  !=0)*Pbike + (Pbike  !==0) Pwalk  <- (Pwalk!=0)*Pwalk   + (Pwalk  ==0) ## Choice results Ctrain Cbus Ccar <- Data$MainModeENG =="Rail" <- Data$MainModeENG =="Bus" <- Data$MainModeENG =="Car" <- Data$MainModeENG =="Bicycle" <- Data$MainModeENG =="Walk" Cbike Cwalk ## Calculate the Log-likelihood function LL <- colSums(Ctrain*log(Ptrain) + Cbus*log(Pbus) + Ccar  *log(Pcar)   + Cbike  *log(Pbike) +Cwalk *log(Pwalk)) exp( μV; ) P(i) = L j ∈ A exp( µ V i ) exp(μVi) P(i) ll Σ j ∈ A exp(μVi)選択確率が0だと尤度も0になるするとln 0 = −∞になってしまいエラーが出る選択確率0の時，1にする処理をするN Yin log Pn(i|0) 4
##Parameter optimization res <= optim(b0,fr,gr=NULL, method = "Nelder-Mead", hessian = TRUE, control=list(fnscale=-1))
●最適化関数optim
●  optim(par, fn, gr = NULL, method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"), lower = -Inf, upper = Inf, control = list(), hessian = FALSE, ...) ● Parというパラメータを初期値としてfnを最小化するようにパラメータを動かしながら反復して探索する（最大化の場合はfnscale =-1） ●探索方法をmethodから選ぶ-   “Nelder-Mead”法：関数値だけを用い頑健(初期値の選択に敏感)だが遅い- “BFGS”法：準ニュートン法．関数値と勾配関数を関数の曲面近似に使う“CG”法：共役勾配法．破綻しやすいがメモリ使用量が少ない- “L-BFGS-B”法：変数の上限・下限を設定した準ニュートン法- “SANN”法：確率的手法である焼きなまし法.関数値だけを用い遅い.
## Parmeter estimation、Hessian matrix calculation b    <- res$par hhh <- res$hessian ## Calculate the t-statistic tval <- b/sqrt(-diag(solve(hhh))) ## L(0), Log-Likelihood when all parameters are 0 L0 <- fr(b0) ## LL, maximium likelihood LL <- res$value #### Output #### print(res) ## L(0) print(L0) ## LL print(LL) ##rho-square print((L0-LL)/L0) ## adjusted rho-square print((L0-(LL-length(b)))/L0) ##estimated parameter values print(b) ## t-statistic print(tval) ○ t値●母集団の平均値の仮説検定に利用● |t| >1.96で5%有意●ヘッセ行列の逆行列の対角成分の平方根が分母●修正済み尤度比LO-(LL-パラメータ数) LO
fr <- function(x) { ### declare the parameters### ## Alternative specific constants b1 <- x[1] b2 <- x[2] b3 <- x[3] b4 <- x[4] ## Travel time to destination d1 <- x[5] ## Fare f1 <- x[6] ##Scale parameter pa <- x[7] MNL同様にベクトルxを入れると尤度を出す関数frを作るx[1]~x[4]は定数項（MNLと同じ） x[5]は時間：上位ネストと下位ネストの組み合わせで決まる効用x[6]は運賃：上位ネストと下位ネストの組み合わせで決まる効用（今回は考慮しない） x[7]はスケールパラメータNLではスケールパラメータを推定スケールパラメータが1ならばMNLと同等

●各選択肢の効用の定義
train bus car bike <- Data$ModeAvailableBike *exp(d1*Data$TimeBike/100 walk <- Data$ModeAvailableWalk *exp(d1*Data$TimeWalk/100 <- Data$ModeAvailableTrain*exp(d1*Data$TotalTimeTrain/100 <- Data$ModeAvailableBus *exp(d1*Data$TotalTimeBus/100 <- Data$ModeAvailableCar  *exp(d1*Data$TimeCar/100 +f1*Data$FareTrain/1000 +f1*Data$FareBus/1000 +b1*matrix(1,nrow =hh,ncol=1)) +b2*matrix(1,nrow =hh,ncol=1)) +b3*matrix(1,nrow =hh,ncol=1)) +b4*matrix(1,nrow =hh,ncol=1)) ●今回はf1で公共交通の運賃を定義しています
ネストの作成：上位ネストの計算V in = - P(m) 48 49 50 51 52 ### Construct nests.First, calculate public or private choice probabilityP(LV1) ## setting of log-sum variables logsum.public  <- log( ( (train+bus)!=0) * (train + bus) + ((train+bus)==0)) logsum.private <- log( ( (car + bike + walk)!=0) * (car + bike + walk) + ((car + bike + walk)==0))
公共／私有の選択についてログサム変数を設定
53 54 ##Calculate public or private choice probability with these log-sum variables nume.public   <- (logsum.public != 0)*exp(pa*logsum.public) + (logsum.public == 0) nume.private  <- (logsum.private != 0)*exp(pa*logsum.private) + (logsum.private == 0) 55ログサム変数を使って公共／私有の選択確率の分子を計算57 58 59 deno <- nume.public + nume.private P.public  <- nume.public  / deno P.private <- nume.private / deno分母を計算し(deno)，公共／私有の選択確率を計算
ネストの作成：上位ネスト内の条件付確率exp( μ m V i ) P(i | m) = F exp(μ mV j) 61 62 63 64 65 ### Next, calculate conditional probability.P(LV2/LV1) ## conditional probability under public choice deno.public      <- train + bus P.train.public  <= train / ((deno.public!=0)*deno.public + (deno.public==0)) P.bus.public    <- bus   / ((deno.public!=0)*deno.public + (deno.public==0))公共一鉄道／バスというネスト構造公共を条件とした鉄道,公共を条件としたバスの選択確率をそれぞれ計算67 ## conditional probability under private choice 68 deno.private  <- car + bike + walk 69 70 71. P.car.private <- car   / ((deno.private!=0)*deno.private + (deno.private==0)) P.bike.private  <- bike  / ((deno.private!=0)*deno.private + (deno.private==0)) P.walk.private  <- walk  / ((deno.private!=0)*deno.private + (deno.private==0))私有一車／自転車／徒歩　というネスト構造私有を条件とした車，自転車，徒歩の選択確率をそれぞれ計算
11

同時確率の計算### Finally, calculate joint probability. P(LV1, LV2) = P(LV2/LV1)*P(LV1) P.train <- P.train.public * P.public <- P.bus.public   * P.public P.bus <= P.car.private  * P.private P . car <- P.bike.private * P.private P.bike <- P.walk.private * P.private P.walk exp( μ m V i ) P(i,m) = P(i | m)P(m) = − − − exp(μ m V i)    exp(μ m V i)ネストの作り方を含めたNLモデルの詳細については2日目の倉内先生の資料を参照ください
●最適化関数optimによる最尤推定res <- optim(b0,fr, method = "L-BFGS-B", lower = c(-Inf,-Inf,-Inf,-Inf,-Inf(,0) upper = c(Inf,Inf,Inf,Inf,Inf(1), hessian = TRUE, control=list(fnscate=-1)) ●ここでは7番目のパラメータがスケールパラメータなので,変数に下限０上限１を設定して，準ニュートン法でパラメータを探索●スケールパラメータが(0,1)の範囲に収まらない場合はネスト構造が不適切●同様に１に貼り付いている場合は不適切であると同時にMNLと同義
13

Fosgerau, Mogens, Emma Frejinger, and Anders Karlstrom. "A link based network route choice model with unrestricted choice set." Transportation Research Part B: Methodological 56 (2013): 70-80. ●自宅から最寄り駅に行くルートをランダム効用最大化理論で考える● MNLモデルでは経路の羅列が不可能（組み合わせ爆発を起こすため）出発

●解決方法
●代表的なルートのみ抽出する：明示的な経路集合を定める●逐次的な行動モデルを考えることで明示的な経路列挙を必要としない
●即時効用，誤差項，目的地までの期待最大効用に従って各ノードで選択●選択そのものはMNLモデル●逐次的かつ再帰的な方法で経路が選択される●逐次的：リンク選択の積み重ねが経路全体の選択になると考える●再帰的：次の状態だけでなく将来の状態まで考えた効用の式...Bellman方程式の利用●マルコフ過程に従うと仮定する●選択肢は無限集合であり，ループを含む経路を選択する可能性がある●道路ネットワークが与えられている場合，経路集合を設定しなくとも良い-アクティビティパス選択などネットワークを限定する必要がある場合はまた別●有向連結グラフ（ネットワーク）G = (S,A)を考える● S:ノードの集合，A:リンクの集合●リンクk,a ∈ A ●リンクkの終点から出るリンクの集合A(k),  a ∈ A(k) ∈ A ● OからDの経路はリンクの列[k o ...k i ]で表される∀i < I,k i+1 ∈ A(k i ), ●効用関数●即時効用un(a|k) = νn(a|k) + μεn(a) ●期待効用   Vnd(a)  ●終点のダミーリンクd ●目的地ノード(D)からダミーリンクdを追加し，吸収状態を定義●全リンクの集合をÃď = A ∪ d ●目的地ノードを終点とする全てのリンクkに対してv,(d|k) = 0を仮定する●主体はマルコフ過程に従って次の状態を選択する●各リンクkでランダム効用εn(a)を観測●即時効用と期待効用の和で表される効用を最大化する  νη(α|k) + μεη(α) + ναα)
●効用関数
●即時効用   un(a|k) = vn(a|k) + μεn(a) -  旅行者nがリンクkにいる状態で選択肢集合A(k)からリンクαを選択する- 確定項  vn(a|k) = vn(xn,a|k；β)誤差項 μεπ(α)：平均値0のガンベル分布を仮定●期待効用   Vnd (a) -リンクαを選択した場合の目的地Dから伸びるダミーリンクdまでの期待効用●効用関数●即時効用un(a|k) = νn(a|k) + μεn(a) ●期待効用Vn (a) O価値関数V đ (k): Bellman方程式と呼ぶリンクkにおいて次にリンクaを選んだ場合の効用max aEA(k νn(a|k) + Vn (a) + μεn(a) Ak E A

21

●以下に示す要素Mから構成される行列M (|A|×|A|)を定義する

22

●  Bellman方程式を解く  Vg (k) = E |    max ( vn(a|k) + Vn(a) + μεn(a) | aEA(k) ●さらに以下の要素からなるベクトルb (|A|×1)とz (|A|×1)を定義する
●すると式(RL 5)は以下のように簡潔な線形方程式で書き下すことができる」は単位行列
● zが解を持つ条件はI − Mが逆行列を持つこと-  逆行列を持たない可能性もある．持つかどうかはパスの数や瞬間効用のサイズに依存する

23

●選択確率に対応する行列Pを作る●目的地に対するリンクの選択確率は出発点に依存しないことに注意●式(3)をA内のリンクについて定義した行列Pに整理する●状態k に対応する行は以下の通り。 ○は要素ごとの積， MkはMのk列目を表す
行列Mの要素
otherwise
●経路選択確率●経路σ = {ki}i=0 = {ko,k1,..,k1}について考える  koは出発リンク,  k1 = d  とする●経路σの選択確率は
●経路選択確率● ν(σ) = ∑ i=1 ν(k i+1|k i )と置くと
●ここでΩは全選択肢の集合●式(10)の分子は経路σの確定項，分母は可能な経路の確定項の和であるため，このモデルは無限の選択肢を持つMNLモデルと同等●従ってIIA特性を持つ
● Bellman方程式に時間割引率βを導入
Sj = d
27

●時間割引率を含んだ価値関数の計算方法
RLモデルと異なり，価値関数の解は非線形方程式系に従うことになる

28

●時間割引率を含んだ価値関数の計算方法
X(z)は|A|×|A|の行列で要素X(z)s j = (zs j ) β式(dRL 2)を解くには価値関数が不動点に達するまで反復計算する必要がある
29

●時間割引率を含んだ価値関数の計算方法→不動点までの反復計算z = MX(z) + b (dRL 2) X(z)は|A|×|A|の行列で要素X(z)s j = (zsj)β反復計算の手順（Mai, Fosgerau, and Frejinger, 2015 ） 1.ベクトルz(0)として初期化する2.     z(1) = MX(z(0) + bとして更新する3.    |z(n+1) — z(n)| < γとなるまで2を繰り返す
データのインプット価値関数（V, z）の初期価の設定対数元度関数の定義•選択確率の計算  P k = MkOzT M k z •対数尤度の計算価値関数の定量•時間割引率を導入しているため価値関数zを閾値まで求解最充街定収長するまで繰り返し価値関数の計算
●インプットデータとして最低限必要なもの●経路データ...トリップID／発リンク／着リンク／終点リンクの列からなる経路の集合●ネットワークデータ...リンクの接続関係を表す●リンクデータ...リンクidごとの説明変数を格納　モデルの説明変数となる
k a 2 6
配布コード：RL_Estimation_Assignment_Node.Rを参照ください！ fr <- function(x) LL <- 0 for(d in 1:D){ dn <- (1:L)[linkid == ddata[d][ dd <- ddata[d] Id <- Ilist[[d]] #  仙街関参z <- array(exp(Theta* V[,d]), dim = c(L, 1)) ZD <-  array(rep(z), dim = c(L, L)) ZD <- t(ZD) #  既問効用M <- Id * Mset(x) ME, dn] <- exp(0) #選択確率を計算Mz <- (M %*% z != 0) * (M %*% z) + (M %*% z == 0) MZ <- array(rep(Mz), dim = c(L, L)) p <- (M * ZD) / MZ p <- (p == 0) * 1 + (p != 0) * p ←終点リンクdごとに計算している選択確率の計算
newV <- function(x) { z <- matrix(1, nrow = L, ncol = D) #exp(Vd) V <- matrix(0, nrow = L, ncol = D) i v- 0 for (d in 1:D){ kd <- (1:L)[linkid == ddata[d]] z[kd, d] <- 1 M <- matrix(0, nrow = L, ncol = L) B <- matrix(0, nrow = L, ncol = 1) B[kd, 1] <- 1 Id <- Ilist[[d]] # 瞬間効用Mを更新for (k in 1:L){ for (a in 1:L){ Ika <- Id[k, a] if(Ika == 1){ if(a == kd){ M[k, a] <- 1 }else{ M[k, a] <- Mset(x)[k, a]関数名：newVの中身←終点リンクdごとに計算している← Mとbを定義←即時効用Mを更新
z = MX(z) + b (dRL 2) X(z)は|A|×|A|の行列で要素X(z)s, = (zs,)β反復計算の手順（Mai, Fosgerau, and Frejinger, 2015 ） 1.ベクトルz(0)として初期化する2.    z(1) = MX(z(0) + bとして更新する3. |z(n+1) - z(n)| < γとなるまで2を繰り返す
35

while(dL>=0.01){ n <- n + 1 b0 <- b Z <- Z0 V <- V0 res <= optim(b0, fr, method = "Nelder-Mead", hessian = TRUE, control=list(fnscale=-1)) b <- res$par cat("b =") print(b) cat("lnL=") print(res$value) #updating value funciton z0 <- newV(b) #convergence condition zz <- (z == 0) * 1 + (z != 0) * z V <- log(zz) zz0 <- (z0 == 0) * 1 + (z0 != 0) * z0 VO <- log(zz0) dL <- sum(abs(b - b0)) # NFXPアルゴリズムcat("dL =") print(dL) ←最尤推定の実施←パラメータの更新←価値関数の求解(newVの実行) ←価値関数の更新
●経路選択モデル(RLモデル) ●  Fosgerau, Mogens, Emma Frejinger, and Anders Karlstrom. "A link based network route choice model with unrestricted choice set." Transportation Research Part B: Methodological 56 (2013): 70- 80. ●災害時におけるドライバーの経路選択モデル(DRLモデル) ●  Oyama, Yuki, and Eiji Hato. "A discounted recursive logit model for dynamic gridlock network analysis." Transportation Research Part C: Emerging Technologies 85 (2017): 509-527. ●  Activity-travelスケジューリングの時空間経路選択モデル（mixed RLモデル） ●  Zimmermann, Maëlle, et al. "Capturing correlation with a mixed recursive logit model for activity- travel scheduling." Transportation Research Part C: Emerging Technologies 93 (2018): 273–291. ●歩行者環境におけるアクセシビリティ分析のための歩行者経路選択モデル（RLモデル） ●  Zheng LIANG*, Ka Fai NG, Hong K. LO “A Link-Based Approach For Measuring Walking Accessibility”, 2022