計量行動分析 第7回
離散選択モデル
ロジットモデル
連続的選択と離散的選択
標準的消費者行動モデル(連続的選択)
一定期間に購入する
財・サービスの量を説明
（お金・時間をいくらずつ割り振るか？）
離散的（質的）選択モデル
どの種類の財を選択するか？
（どこへ行くか？何を買うか？）
離散的(discrete)な
選択肢(alternative)からの選択
5万円
離散的選択のモデル化
個人は，採りうる選択肢alternativeを列挙する
それぞれの選択肢の特徴と費用に基づいて，評
価得点utilityをつける
評価点が高いものを選ぶ
中国旅行 60点
フランス旅行 40点
アメリカ旅行50点
確定的選択モデル
個人は「評価得点が少しでも高いほうの選択肢を
必ず選ぶ」と考えるモデル
全員が同じ考え方で評価
同じ状況に直面する人は，全員同じ選択肢を選択
異なる条件の人の選択結果から効用関数を推定
判別関数モデル，数量化理論II類モデル
同じ状況に対する個人間での考え方の違いを考慮
例：高齢者は着席可能性を，若者は低運賃を重視
犠牲量(最小化）モデル
比較における「あいまいさ」を認める
ファジィ選択モデル
確率的選択：評価点の差と選択率
実際には
ほとんど評価点が同じときは，どちらも選択される可
能性がある
評価点の差が大きいときは，片方しか選ばれない．
選択肢Aが選ばれる可能性
1
2つは同じ魅力
50％ずつ
Aが圧倒的に劣る
Aが選ばれること
はほとんどない
0
Aが圧倒的に良い
ほとんどAだけが
選ばれる
選択肢Aの得点－選択肢Bの得点
ロジットモデル
S字型の曲線として，
という式で表わされる曲線を使うと，
いろいろな計算が簡単にできる
3つ以上の選択肢からの選択も同じ形になる
2000年ノーベル経済学賞
McFadden（1937-）が提案
各自の評価点が安定している部分と確率的に変
動する部分の和である場合の選択から理論的に
導いた。（ランダム効用モデル）
ロジットモデル
選択肢が2つの2項ロジットモデル
Binary （Binomial) Logit
exp(V1 )
exp(V2 )
P1
, P2
1 P1
exp(V1 ) exp(V2 )
exp(V1 ) exp(V2 )
選択肢が多数（ｎ個）ある場合の
多項ロジットモデル Multinomial Logit
Pj
exp(V j )
N
k 1
exp(Vk )
モデルの使用手順
選択モデルの定式化（得点からSカーブで選択率を計算）
パラメータη，βの推定
（実際の選択結果から、どういう得点付けだったかを推測）
将来の選択率の予測
（将来の選択肢の状況から得点を計算し、選択率を計算）
実際の観測結果から、
もとの事象の発生確率を推定する
手品師がイカサマかも知れないコインを6
回振ったら、表が5回、裏が1回出た。
このコインは正しいコインだろうか？表が出や
すいようなイカサマのコインだろうか？
表5回、裏1回という観測結果から、このコ
インを1回振ったときに表が出る確率pの値
を推測したい。
比率を用いてそのままｐ＝5/6と推測する方法
さいゆうほう（最尤法）
比率によるロジットモデルの推定
集計的方法(集団の選択率にあわせる)
ある選択肢の状況下で観測された集団の選
択確率ｐを用いる
ロジットモデルから、その時の2つの選択肢の
魅力度VAとVBの差が逆算できる
魅力度の差がうまく一致するように、魅力度の
関数の形を調整する
先の例：6回のうち5回表だから、p=5/6と推測した。
回帰分析 Linear Regression
y ln[PijCAR / PijBUS ]
y
平面
y x1 x2
をうまく決める
x1 (tijCAR tijBUS )
x2 (cijCAR cijBUS )
x1
x2
2つ(以上)の変数を用いて、目的の変数yの値をうまく予
測できるような平面を決める方法
最尤法による推定
(個々の事象の組み合わせを考える）
コインを続けて6回振ったところ、そのうち5回が表であっ
た．このコインの表の出る確率ｑはいくらか？
母数（ここでは表の出る確率ｑ)を変えたときに，観察され
た事象が実現する確率を求める(尤度関数L(ｑ)）
観察された事象の発生確率が最大になるｑの値を求める
尤度関数を求めてみよう
確率ｑの事象が5回，（１－ｑ）の事象が1回観察されたの
だから、何回目が裏かが6通りあることを踏まえると、
6! 5
L(q) 6 C5 q (1 q)
q (1 q) 6q 5 (1 q )
5!1!
5
1
最尤法による母数の推定の例
尤度関数を最大化する
L(q) 6q (1 q)
5
L(θ)
0.45
0.4
dL(q)
6 5q 4 (1 q) q 5 6q 4 (5 5q q)
dq
0.35
0.3
0.25
0.2
6q 4 (5 6q) 0
0.15
あるいは対数をとったものをｑに
ついて最大化してもよい
0.05
0.1
0
0.00
0.20
0.40
0.60
0.80
1.00
L * (q) ln L(q) ln 6 5 ln q ln(1 q)
dL * 5
1
5(1 q) q 5 6q
0
dq q (1 q)
q(1 q)
q(1 q)
q 5 / 6 8.333 これは，6回のうち5回という割合に一致
最尤法によるロジットモデルの推定
非集計的方法(各個人の選択を考える)
魅力度の関数形を仮定する
各個人が直面した選択肢の状況をもとに、ロ
ジットモデルで各自の選択確率を求める。
それらを掛け合わせて、観測された事象の発
生確率(尤度関数)を求める
L n1 j 1 Pnj
N
J
nj
尤度関数の値が最大になるように、魅力度の
関数の形を調整する
交通手段選択モデル
あるODを移動する消費者が,複数の交通手段
の所要時間や費用を考えて手段を選択
U i Vi j
BUS
ij
V
t
BUS
ij
c
BUS
ij
VijCAR tijCAR cijCAR
2項ロジットモデル
exp( Vi )
Pi Prob[ Ui Uk : k ( i) K ]
k 1,2 exp( Vk )
選択肢jの魅力度が他の選択肢よりも高い確率
2項ロジットモデルの図解
V*
1
V
0
X2
X1
集計データ(比率)を用いた
ロジットモデルの推定
表3.15 バスの所要時間
表3.17 バスの所要費用
ｔBij
cBij
1
2
3
1
5
11
13
2
10
12
3
14
16
表3.19 バスの分担率
1
2
3
1
130
140
180
12
2
140
130
7
3
180
220
PBij
1
2
3
1
0.273
0.265
0.253
220
2
0.282
0.248
0.255
130
3
0.239
0.192
0.244
単位：分
単位：円
表3.16 自動車の所要時
間
表3.18 自動車の所要費用
表3.20 自動車の分担率
cCij
PCij
tCij
1
2
3
1
2
3
1
2
3
1
21
45
58
1
0.727
0.735
0.747
1
3
8
10
2
45
42
60
2
0.718
0.752
0.745
2
8
7
11
3
58
60
19
3
0.761
0.808
0.756
3
10
11
3
単位：分
単位：円
集計データ(比率)を用いた
線形回帰分析による推定
ロジットモデル式から、2つの選択率の比は以下のようになる
ln[PijCAR / PijBUS ] (tijCAR tijBUS ) (cijCAR cijBUS )
ｔ
B ij ｔ
C ij
5
10
14
11
12
16
13
12
7
3
8
10
8
7
11
10
11
3
ｃB ij cC ij PB ij PC ij ln(PC /PB ) ｔ
C -tB
130
21 0.273 0.727
0.979
-2
140
45 0.282 0.718
0.935
-2
180
58 0.239 0.761
1.158
-4
140
45 0.265 0.735
1.020
-3
130
42 0.248 0.752
1.109
-5
220
60 0.192 0.808
1.437
-5
180
58 0.253 0.747
1.083
-3
220
60 0.255 0.745
1.072
-1
130
19 0.244 0.756
1.131
-4
ｃC -cB
-109
-95
-122
-95
-88
-160
-122
-160
-111
集計データ(比率)を用いた
線形回帰分析による推定
係数 標準誤差 t
切片 0.3898 0.045 8.731
X 値 1-0.0796 0.006 -13.2
X 値 2-0.0039 3E-04 -12.2
回帰統計
重相関0.
R9899
重決定0.
R9799
2
補正 R 02.9732
標準誤差
0.0237
観測数
9
V
0.0796t
CAR
ij
-1
P-値 下限 95%上限 95%下限 95.
上限
0% 95.0%
1E-04 0.281 0.499055 0.281 0.499
1E-05 -0.09 -0.0648 -0.09 -0.06
2E-05
-0 -0.00309
-0
-0
バスの分担率
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
VijBUS 0.0796tijBUS 0.00387cijBUS
CAR
ij
-2
0.00387c
CAR
ij
0.1
0
0.390
-2
-1
0
1
VBUS-VCAR
2
Rによる線形回帰
tb <- c(5,10,14,11,12,16,13,12,7)
tc <- c(3,8,10,8,7,11,10,11,3)
cb <- c(130,140,180,140,130,220,180,220,130)
cc <- c(21,45,58,45,42,60,58,60,19)
pb <- c(0.273,0.282,0.239,0.265,0.248,0.192,0.253,0.255,0.244)
pc <- rep(1,9)-pb
lnpbpc <- log(pc/pb)
tsa <- tc-tb
csa <- cc-cb
ans <- lm(lnpbpc ~ tsa+csa)
summary(ans)
Rによる線形回帰の推定結果
Call:
lm(formula = lnpbpc ~ tsa + csa)
Residuals:
Min
1Q Median
3Q
Max
-0.021913 -0.017819 -0.006659 0.018098 0.030403
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3898049 0.0446483 8.731 0.000125 ***
tsa
-0.0795878 0.0060416 -13.173 1.18e-05 ***
csa
-0.0038682 0.0003171 -12.200 1.85e-05 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0237 on 6 degrees of freedom
Multiple R-squared: 0.9799,
Adjusted R-squared: 0.9732
F-statistic: 146 on 2 and 6 DF, p-value: 8.165e-06
最尤法による
ロジットモデルの推定
表3.15 バスの所要時間
表3.17 バスの所要費用
ｔBij
cBij
1
2
3
1
2
3
1
5
11
13
1
130
140
180
2
10
12
12
2
140
130
220
3
14
16
7
3
180
220
130
単位：分
単位：円
表3.16 自動車の所要時
間
表3.18 自動車の所要費用
tCij
cCij
1
2
3
1
2
3
1
21
45
58
1
3
8
10
2
45
42
60
2
8
7
11
3
58
60
19
3
10
11
3
単位：分
単位：円
表 バスの利用者数
NBij
1
2
39
22
1
11
31
2
16
15
3
3
21
25
50
表 自動車の利用者数
NCij
1
2
3
104
61
62
1
28
94
73
2
51
63 155
3
最尤法による
ロジットモデルの推定
時間 費用
選択人数 魅力度(仮定値) 個人の選択率 事象発生確率
ｔBij ｔCij ｃBij cCij NBij NCij Vbij
Vcij
Pbij Pcij
Pb^Nb Pc^Nc
5 3 130 21 39 104 -0.8878 0.0831 0.27 0.72529 1.3E-22 3.1E-15
10 8 140 45 11
28 -1.3163 -0.399 0.29 0.71449 1E-06 8.2E-05
14 10 180 58 16
51 -1.7816 -0.605 0.24 0.76436 9E-11 1.1E-06
11 8 140 45 22
61 -1.3944 -0.399 0.27 0.73014 3.1E-13 4.7E-09
12 7 130 42 31
94 -1.4341
-0.31 0.25 0.75484 1.2E-19 3.3E-12
16 11 220 60 15
63 -2.0908 -0.691 0.2 0.80222 2.8E-11 9.3E-07
13 10 180 58 21
62 -1.7035 -0.605 0.25 0.75001 2.3E-13 1.8E-08
12 11 220 60 25
73 -1.7786 -0.691 0.25 0.74801 1.1E-15 6.2E-10
7 3 130 19 50 155 -1.0439 0.0907 0.24 0.75669 2E-31 1.7E-19
魅力度関数の係数の仮定値
-0.078045
α
6E-139 7.8E-87
-0.003828
β
尤度関数 尤度関数 5E-225
0.397569
γ
最大になるように値を調整
最尤法による
ロジットモデルの推定
ソルバー機能：
表計算ソフト：Excel の中では、いくつかの数値が
別のセルの数値に影響を与える場合、
目的のセルの数値を最大にするように、元の数値
を決定する計算ができる。
実際には、この尤度関数は、あまりに値が小さい
ため、数値計算誤差の影響でうまく計算できない。
尤度関数の対数(log)を取ったものを考え、それを
最大化する。
ln n1 j 1 Pnj
N
J
nj
N
n 1
ln
P
nj
nj
j 1
J
最尤法による
ロジットモデルの推定
時間 費用
選択人数 魅力度(仮定値) 個人の選択率選択率の対数 人数分を合計
ｔBij ｔCij ｃBij cCij NBij NCij Vbij
Vcij
Pbij Pcij
lnPbij lnPcij NblnPb NclnPc
5 3 130 21 39 104 -0.8878 0.0831 0.27 0.725 -1.292 -0.321 -50.39 -33.4
10 8 140 45 11
28 -1.3163 -0.399 0.29 0.714 -1.253 -0.336 -13.79 -9.413
14 10 180 58 16
51 -1.7816 -0.605 0.24 0.764 -1.445 -0.269 -23.13 -13.7
11 8 140 45 22
61 -1.3944 -0.399 0.27
0.73 -1.31 -0.315 -28.82 -19.19
12 7 130 42 31
94 -1.4341
-0.31 0.25 0.755 -1.406 -0.281 -43.58 -26.44
16 11 220 60 15
63 -2.0908 -0.691 0.2 0.802 -1.621 -0.22 -24.31 -13.88
13 10 180 58 21
62 -1.7035 -0.605 0.25
0.75 -1.386 -0.288 -29.11 -17.84
12 11 220 60 25
73 -1.7786 -0.691 0.25 0.748 -1.378 -0.29 -34.46 -21.19
7 3 130 19 50 155 -1.0439 0.0907 0.24 0.757 -1.413 -0.279 -70.67 -43.21
魅力度関数の係数の仮定値
-0.078045
α
-318.3 -198.3
-0.003828
β
対数尤度関数 -516.5
0.397569
γ
最大になるように値を調整
作成されたロジットモデル
V BUS 0.078t BUS 0.0038c BUS
V CAR 0.078t CAR 0.0038c CAR 0.3975
PBUS
exp(VBUS )
exp(VBUS ) exp(VCAR )
PCAR 1 PBUS
Rにおける
尤度関数の定義と最大化
#対数尤度関数の定義 係数ベクトルｘの関数と見なす．
LL <- function (x){
vbus <- x[1] * tb + x[2] * cb
vcar <- x[1] * tc + x[2] * cc + x[3]
ppb <- 1/(1+exp(vcar - vbus))
ppc <- 1- ppb
return(sum(pb*log(ppb)+pc*log(ppc)))
}
#Optim関数で最大化を行い，結果をresに代入する．初期値は(0,0,0)
b0=c(0,0,0)
res<-optim(b0,LL, method = "BFGS", hessian = TRUE, control=list(fnscale=-1))
Rによる推定結果の表示
> res
$par
[1] -0.081037584 -0.004007811 0.369193890
$value
[1] -5.047779
$counts
function gradient
62
18
$convergence
[1] 0
$message
NULL
$hessian
[,1]
[,2]
[,3]
[1,] -19.628306 -613.3939 5.313007
[2,] -613.393875 -23980.3173 196.530577
[3,] 5.313007 196.5306 -1.681972
Rによる非集計最尤法の推定
完全非集計データ(1サンプル1行）の時
glm（従属変数～独立変数１＋・・・＋独立変
数n,family=binomial(link=“logit”))
従属変数が比率の場合
glm（従属変数～独立変数１＋・・・＋独立変
数n,family=quasibinomial(link=“logit”))
Rによる最尤法での推定
比率を説明する場合
ans2 <- glm(pc ~ tsa+csa, family=quasibinomial(link="logit"))
summary(ans2)
Call:
glm(formula = pc ~ tsa + csa, family = quasibinomial(link = "logit"))
Deviance Residuals:
Min
1Q
Median
3Q
Max
-0.008856 -0.007615 -0.002659 0.007188 0.013667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3991916 0.0446219 8.946 0.000109 ***
tsa
-0.0787542 0.0060030 -13.119 1.21e-05 ***
csa
-0.0038095 0.0003174 -12.001 2.03e-05 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasibinomial family taken to be 0.0001007740)
各個人の選択を説明する場合
Rにより，個人ごとのデータを作成
nb <- c(39,11,16,22,31,15,21,25,50)
nc <- c(104,28,51,61,94,63,62,73,155)
nn <- nb + nc
nnf <-numeric(length=10)
for (i in 1:9){
nnf[i+1]=nnf[i]+nn[i]
}
chc <- numeric(length=nnf[10]-1)
tim <- numeric(length=nnf[10]-1)
cst <- numeric(length=nnf[10]-1)
for (i in 1:9){
chc[nnf[i]:(nnf[i]+nn[i]-1)] <- c(rep(1,nb[i]), rep(0,nc[i]))
tim[nnf[i]:(nnf[i]+nn[i]-1)] <- rep((tb[i]-tc[i]),nn[i])
cst[nnf[i]:(nnf[i]+nn[i]-1)] <- rep((cb[i]-cc[i]),nn[i])
各個人の選択を説明する場合
ans3 <- glm(chc ~ tim+cst, family=binomial(link="logit"))
summary(ans3)
Call:
glm(formula = chc ~ tim + cst, family = binomial(link = "logit"))
Deviance Residuals:
Min
1Q Median
3Q
Max
-0.8215 -0.7630 -0.7470 -0.1019 1.8016
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.385889 0.507795 -0.760 0.447
tim
-0.079514 0.061803 -1.287 0.198
cst
-0.003873 0.003465 -1.118 0.264
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1034.7 on 919 degrees of freedom
Residual deviance: 1032.4 on 917 degrees of freedom
AIC: 1038.4
Number of Fisher Scoring iterations: 4
理論的補足
ランダム項がある場合の選択
ランダム項の確率分布の仮定
集団全体における各自の効用の誤差εの頻度分
布を連続的な確率密度関数で表現
選択率の導出
選択率の導出(図解)
ε2
V1-V2+ε1
V1-V2
（＞０）
ε1
V1-V2（＜０）の場合
ε2
ε1
選択肢１が選ばれる確率
Prob(U1>U2)
=Prob(V1+ε1>V2＋ε２）
＝Prob(ε2<V1-V2+ε1)
誤差項をどのような分布で考えるか？
(まんじゅうのふくらみ方)
釣鐘（つりがね）まんじゅうを，
中央からずれた位置で包丁
で切る
右下の方の切れ端の体積
が選択率を表わす．
具体的なランダム効用モデル
ロジット（LOGIT)モデル
各選択肢の誤差項が独立で同一の
Gumbel分布に従う
F ( x) Prob( x) exp exp(x)
f ( x) F ' ( x) exp(x) F ( x)
と仮定したモデル
第k選択肢の選択確率は，
Pi exp(Vi ) / exp(Vk )
k
η：誤差項の分散に
反比例するパラメータ
プロビット（PROBIT)モデル
誤差項が多変量正規分布に従うと仮定したモデル
選択確率を解析的に表示することができない
モンテカルロシミュレーションなどを用い近似値を計算
非集計データによる
ロジットモデルのパラメータ推定
尤度関数
尤度関数の値を最大にするようにβを定める
パラメータ推定（２）
ニュートン・ラプソン法
最尤法の計算例
通勤に2 種類の交通手段(1：バ
ス、2：鉄道) があり、10 人の両交
通機関での所要時間を調べたと
ころ、表のようであった。
ロジットモデルを適用して、所要
時間のパラメータとバスの選択し
定数項のパラメータを計算せよ。
さらに、バス専用レーンの設置に
より、全てのサンプルのパス所要
時間が1 分短縮された場合に、
バス選択確率がどのように変化
するかを推測せよ。
サンプル 選択 所要時間
番号
バス バス 鉄道
1
0
22
22
2
0
24
20
3
0
23
19
4
0
21
20
5
0
22
20
6
0
24
21
7
1
21
20
8
1
20
20
9
1
23
25
10
1
20
23
合計→
4
22
21
時間差
Ｚ1-Ｚ2
0
4
4
1
2
3
1
0
-2
-3
←平均
Rによる計算
use <- c(0,0,0,0,0,0,1,1,1,1)
tbus <- c(22,24,23,21,22,24,21,20,23,20)
trail <- c(22,20,19,20,20,21,20,20,25,23)
tdiff <- tbus - trail
ans4 <- glm(use ~ tdiff, family=binomial(link="logit"))
summary(ans4)
Rによる計算結果
Call:
glm(formula = use ~ tdiff, family = binomial(link = "logit"))
Deviance Residuals:
Min
1Q Median
3Q
Max
-1.4457 -0.3896 -0.1086 0.2146 1.5412
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6117
1.1632 0.526 0.599
tdiff
-1.4357
1.0207 -1.406 0.160
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 13.460 on 9 degrees of freedom
Residual deviance: 6.406 on 8 degrees of freedom
AIC: 10.406
Number of Fisher Scoring iterations: 6
尤度関数の数値計算例
alpha/beta
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
alpha/beta
0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
0.0044
0.01169
0.02091
0.02873
0.03337
0.03485
0.03397
0.03163
0.02853
0.02517
0.02184
0.01872
0.00413
0.00378
0.01141
0.01091
0.02099
0.02071
0.0294
0.02964
0.03461
0.03544
0.03651
0.03784
0.03588
0.03754
0.03364
0.03546
0.03052
0.03239
0.02706
0.02887
0.02358
0.02527
0.02028
0.02182
0.00338
0.01022
0.02008
0.02945
0.03585
0.03878
0.03888
0.03704
0.03407
0.03055
0.02688
0.02332
0.00296
0.00252
0.00937
0.00841
0.01913
0.01789
0.02883
0.0278
0.03579
0.03527
0.03929
0.03936
0.03986
0.04044
0.03834
0.03931
0.03555
0.03676
0.03209
0.03342
0.02839
0.02976
0.02474
0.02606
0.0021
0.0074
0.01644
0.0264
0.0343
0.03898
0.04061
0.03993
0.03769
0.03454
0.03095
0.02726
0.00171
0.00136
0.00637
0.00537
0.01484
0.01315
0.02469
0.02273
0.03292
0.03118
0.03815
0.0369
0.04036
0.03968
0.04016
0.04001
0.0383
0.03857
0.03539
0.03597
0.03195
0.03272
0.02831
0.02918
0.00106
0.00443
0.01145
0.02061
0.02914
0.03527
0.03859
0.03947
0.03849
0.03625
0.03325
0.02987
0.00081
0.00359
0.00979
0.0184
0.02687
0.03331
0.03714
0.03856
0.03807
0.03623
0.03353
0.03034
1
1.05
1.1
1.15
1.2
1.25
1.3
1.35
1.4
1.45
1.5
1.55
0.03527
0.03661
0.03773
0.03865
0.03936
0.03989
0.04023
0.04041
0.04044
0.04033
0.04010
0.03976
0.03484
0.03625
0.03745
0.03844
0.03923
0.03982
0.04024
0.04049
0.04058
0.04053
0.04036
0.04006
0.03430
0.03578
0.03705
0.03812
0.03898
0.03965
0.04014
0.04045
0.04061
0.04063
0.04051
0.04027
0.03366
0.03521
0.03655
0.03768
0.03862
0.03936
0.03992
0.04031
0.04054
0.04062
0.04056
0.04038
0.03292
0.03453
0.03594
0.03714
0.03815
0.03896
0.03960
0.04006
0.04036
0.04050
0.04051
0.04040
0.03210
0.03376
0.03522
0.03650
0.03757
0.03846
0.03917
0.03970
0.04007
0.04029
0.04036
0.04031
0.03118
0.03289
0.03442
0.03575
0.03690
0.03785
0.03863
0.03924
0.03968
0.03996
0.04011
0.04012
0.03019
0.03195
0.03353
0.03492
0.03613
0.03715
0.03800
0.03867
0.03918
0.03954
0.03975
0.03983
尤度関数の
数値計算例(ソルバー)
alpha/beta 1.436
0.611735 0.04064
=1/(
(1+EXP(-B$24*0+$A25))*
(1+EXP(-B$24*4+$A25))*
(1+EXP(-B$24*4+$A25))*
(1+EXP(-B$24*1+$A25))*
(1+EXP(-B$24*2+$A25))*
(1+EXP(-B$24*3+$A25))*
(1+EXP(B$24*1-$A25))*
(1+EXP(B$24*0-$A25))*
(1+EXP(B$24*(-2)-$A25))*
(1+EXP(B$24*(-3)-$A25))
)
得られたパラメータを用いた予測
現況再現
α
0.612 β -1.44
Ｖ１
Ｖ２ Ｐ１ Ｐ２ 予測
-30.97 -31.6 0.65 0.352
-33.84 -28.7 0.01 0.994
-32.41 -27.3 0.01 0.994
-29.54 -28.7 0.30 0.695
-30.97 -28.7 0.09 0.905
-33.84 -30.1 0.02 0.976
-29.54 -28.7 0.30 0.695
-28.10 -28.7 0.65 0.352
-32.41 -35.9 0.97 0.03
-28.10
-33 0.99 0.007
バスを利用するサンプルの予測値
政策実施後
1
0
0
0
0
0
0
1
1
1
4
（レーン設置後）
Ｖ１
Ｖ２ Ｐ１ Ｐ２ 予測
-29.54 -31.6 0.89 0.114
-32.41 -28.7 0.02 0.976
-30.97 -27.3 0.02 0.976
-28.10 -28.7 0.65 0.352
-29.54 -28.7 0.30 0.695
-32.41 -30.1 0.09 0.905
-28.10 -28.7 0.65 0.352
-26.67 -28.7 0.89 0.114
-30.97 -35.9 0.99 0.007
-26.67
-33 1.00 0.002
バスを利用するサンプルの予測値
1
0
0
1
0
0
1
1
1
1
6
多項ロジットモデル
Multinomial Logit Model
Pi exp(Vi ) / exp(Vk )
k
モデル作成の手順
多項ロジットモデル
Multinomial Logit Model
個人によって
全ての選択肢が
利用できない
場合もありうる
Pi exp(Vi ) / exp(Vk )
k
推定結果の評価
的中率
モデル上の選択結果と，実際の選択結果との当てはまりを示す指標
多項ロジットモデル用のデータ
(愛媛大学作成：)ehime.csv
SEQ,選択,年齢,保有台数,鉄道可能性,時間,費用,バス可能性,時間,費用,４輪可能性,時間,費用
3869,1,7,5,1,37.2 ,250,1,129.0 ,850,1,25.61,120.61
7447,1,7,5,1,37.2 ,250,1,121.3 ,1174.33,1,25.61,120.61
7924,1,7,5,1,42.1 ,230,1,148.3 ,960,1,29.85,142.34
2460,1,5,5,1,39.0 ,230,1,106.6 ,970,1,37.06,177.84
3800,1,5,5,1,39.0 ,230,1,91.3 ,850,1,37.06,177.84
3347,1,3,5,1,22.0 ,140,1,77.9 ,300,1,9.41,34.04
4143,1,3,5,1,22.0 ,140,1,77.9 ,300,1,9.41,34.04
2529,1,1,5,1,45.0 ,320,1,94.8 ,910,1,34.42,165.4
4985,1,1,5,1,45.0 ,320,1,128.5 ,1050,1,34.42,165.39
6424,1,6,4,1,83.3 ,570,1,179.8 ,1604.61,1,62.26,511.51
7791,1,6,4,1,83.3 ,570,1,172.0 ,1720,1,62.26,511.51
1498,1,5,4,1,43.6 ,650,1,169.3 ,1932.33,1,77.59,400.49
1962,1,5,4,1,54.1 ,609,1,136.0 ,1182.28,1,61.78,313.43
2367,1,5,4,1,120.3 ,1100,1,112.2 ,1449.06,1,60.44,306.36
2369,1,5,4,1,120.3 ,1100,1,112.2 ,1449.06,1,60.44,306.36
3985,1,5,4,1,54.1 ,609,1,134.2 ,1182.28,1,61.78,313.43
4307,1,5,4,1,43.6 ,650,1,183.8 ,2082.33,1,77.59,400.49
4931,1,5,4,1,155.0 ,1100,1,143.1 ,1449.06,1,60.44,306.36
4932,1,5,4,1,155.0 ,1100,1,143.1 ,1449.06,1,60.44,306.36
4701,1,1,4,1,54.7 ,320,1,105.9 ,1020,1,31.65,160.45
5774,1,1,4,1,29.7 ,180,1,67.4 ,450,1,14.63,61.37
6840,1,1,4,1,29.7 ,180,1,97.8 ,450,1,14.63,61.37
7242,1,1,4,1,54.7 ,320,1,98.9 ,882.17,1,31.65,160.45
2368,1,7,3,1,120.3 ,1100,1,112.2 ,1449.06,1,60.44,306.36
3309,1,7,3,1,75.7 ,570,1,167.2 ,1460,1,58.17,497.65
4662,1,7,3,1,18.0 ,140,1,70.4 ,160,1,6.81,20.01
4830,1,7,3,1,35.3 ,296,1,95.3 ,440,1,20.37,92.23
4933,1,7,3,1,155.0 ,1100,1,143.1 ,1449.06,1,60.44,306.36
4980,1,7,3,1,33.7 ,190,1,125.1 ,800,1,22.24,102.38
多項ロジットモデルの最尤法プ
ログラム(愛媛大学作成：)
### Multi Nomial Logit (MNL) estimation program (Original code by EHIME University)
###データファイルの読み込み
Data<-read.csv("ehime.csv",header=T)
hh<-nrow(Data) ##データ数:Data の行数を数える
print(hh)
ch<- 3 ##今回用いる選択肢の数
b0<-c(0, 0, 0, 0, 0, 0)
Srail <- sum(Data[,14]== 1); Sbus <- sum(Data[,14]== 2); Scar <- sum(Data[,14]== 3)
cat("rail:",Srail," bus:",Sbus," car:",Scar,"\n")
## Logit model の対数尤度関数の定義
fr <- function(x) {
LL=0
##効用の計算
rail <- x[1]*Data[, 6]/100 + x[2]*Data[, 7]/100 + x[5]*matrix(1,nrow =hh,ncol=1)
bus <- x[1]*Data[, 9]/100 + x[2]*Data[,10]/100 + x[3]*(Data[, 3]>=6) + x[6]*matrix(1,nrow=hh,ncol=1)
car <- x[1]*Data[,12]/100 + x[2]*Data[,13]/100 + x[4]*(Data[, 4]>=2)
##効用の指数化
Erail <- exp(rail)*Data[, 5]
Ebus <- exp(bus )*Data[, 8]
Ecar <- exp(car )*Data[,11]
PPrail <- Erail/(Erail+Ebus+Ecar)
PPbus <- Ebus /(Erail+Ebus+Ecar)
PPcar <- Ecar /(Erail+Ebus+Ecar)
##選択結果の確率のみを有効化
Prail <- (PPrail!=0)*PPrail + (PPrail==0)
Pbus <- (PPbus !=0)*PPbus + (PPbus ==0)
Pcar <- (PPcar !=0)*PPcar + (PPcar ==0)
##選択結果
Crail <- Data[,14]== 1
多項ロジットモデルの
最尤法による推定結果
rail: 493 bus: 432 car: 708
roh = 0.2432912
rohbar= 0.2398755
$par
[1] -1.7411817 -0.1757195 1.7302273 2.5890795 2.0644240 2.3457336
$value
[1] -1329.232
$counts
function gradient
38 14
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -80.15815 -538.9584 -67.42988 65.52220 27.40143 -114.7780
[2,] -538.95843 -4708.3230 -493.74059 503.37025 123.55640 -800.5406
[3,] -67.42988 -493.7406 -127.96682 31.64655 79.67259 -127.9668
[4,] 65.52220 503.3703 31.64655 -214.17601 136.44173 77.7343
[5,] 27.40143 123.5564 79.67259 136.44173 -293.24442 123.6202
[6,] -114.77804 -800.5406 -127.96682 77.73429 123.62021 -224.7471
[1] -5.857365 -5.520454 12.528884 17.136179 13.928912 10.259022
多項ロジットモデルの
最尤法による推定結果
rail: 493 bus: 432 car: 708
roh = 0.2432912
rohbar= 0.2398755
$par
時間
費用
高齢B
台数C
鉄道
バス
[1] -1.7411817 -0.1757195 1.7302273 2.5890795 2.0644240 2.3457336
[1] -5.857365
-5.520454 12.528884 17.136179 13.928912 10.259022
$value
[1] -1329.232
所要時間は -0.0174［／分］，費用は -0.00176［／円］でパラメータ比から得られる時間
評価値は10［円／分］程度と，若干低い
多項ロジットモデル （パッケージ）
Multinomial Logit Model
install.packages("mlogit")
library("mlogit")
Dehime<-read.csv("ehime.csv",header=T)
Ehime <-mlogit.data(Dehime,varying=c(5:13),
shape="wide",choice="choice")
#MNL without constant term
summary(mlogit(choice~time+cost-1,data=Ehime))
#MNL with constant term
summary(mlogit(choice~time+cost,data=Ehime))
多項ロジットモデル(定数項なし)
Multinomial Logit Model
Call:
mlogit(formula = Choice ~ time + cost - 1, data = Ehime)
Frequencies of alternatives:
bus
car rail
0.26454 0.43356 0.30190
Newton-Raphson maximisation
gradient close to zero. May be a solution
5 iterations, 0h:0m:0s
g'(-H)^-1g = 7.47E-31
Coefficients :
Estimate
Std. Error
t-value Pr(>|t|)
time -0.00329405 0.00197725 -1.6660 0.0957185 .
cost -0.00096882 0.00025323 -3.8259 0.0001303 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Log-Likelihood: -1715.7
多項ロジットモデル(定数項あり)
Multinomial Logit Model
Call:
mlogit(formula = Choice ~ time + cost, data = Ehime)
Frequencies of alternatives:
bus
car rail
0.26454 0.43356 0.30190
Newton-Raphson maximisation
gradient close to zero. May be a solution
5 iterations, 0h:0m:0s
g'(-H)^-1g = 9.17E-26
Coefficients :
Estimate
Std. Error
t-value Pr(>|t|)
altcar
-1.07985435 0.14856671 -7.2685 3.635e-13 ***
altrail
-0.95903599 0.11625041 -8.2497 2.220e-16 ***
time
-0.01607098 0.00261456 -6.1467 7.910e-10 ***
cost
-0.00119568 0.00027125 -4.4081 1.043e-05 ***
Log-Likelihood: -1680.5
McFadden R^2: 0.68776
Likelihood ratio test : chisq = 152.13 (p.value=< 2.22e-16
MNLモデルの限界と注意点
青バス・赤バス問題（I.I.A特性)
例：バスと自動車の交通機関選択モデル
バスも自動車も共に一般化交通費用が等しい
→バス,自動車の選択率は1/2ずつになる
バスの半数を赤く,半数を青く塗って区別
→青バス,赤バス,自動車の選択率は1/3ずつ
色を変えたらバスのシェアが1/2から2/3に上昇？
各選択肢の確率的効用の間に相関がある場合
には，非現実的な選択率を与える
ネスティッド(入れ子)ロジットモデル
Nested Logit Model （NL)
下位の選択肢のうち
大きいほうを取ったと
きの期待値
ネスティッド(入れ子)ロジットモデル
Nested Logit Model （NL)
下位の選択肢の間
の相関のほうが,上
位の選択肢間の相
関より大きい必要が
ある
ネスティッド(入れ子)ロジットモデル
Nested Logit Model （NL)
•段階推定法：下の階層から順次推定して，ログサ
ム変数値を計算して上位階層を推定する
•階層間でパラメータ値が大きく異なるなど,解釈が困難
•同時推定法：全体の対数尤度関数を作り，全パラメータ
値とスケールパラメータを一気に推定
ネスティッドロジットモデルの最尤
法プログラム(愛媛大学作成：)
### Nested Logit estimation program (Original code by EHIME University)
###データファイルの読み込み
Data<-read.csv(“ehime.csv",header=T)
hh<-nrow(Data) ##データ数:Data の行数を数える
print(hh)
ch<- 3 ##今回用いる選択肢の数
b0<-c(0, 0, 0, 0, 0, 0, 1) ##初期パラメータ値(全て0)
## サンプルにおける各選択肢のシェア
Srail <- sum(Data[,2]==“rail”); Sbus <- sum(Data[,2]==“ bus”); Scar <- sum(Data[,2]==“car”)
cat("rail:",Srail," bus:",Sbus," car:",Scar,"\n")
## Logit model の対数尤度関数の定義
fr <- function(x) {
LL=0
##効用の計算
rail <- x[1]*Data[, 6]/100 + x[2]*Data[, 7]/100 + x[5]*matrix(1,nrow =hh,ncol=1)
bus <- x[1]*Data[, 9]/100 + x[2]*Data[,10]/100 + x[3]*(Data[, 3]>=6) + x[6]*matrix(1,nrow
=hh,ncol=1)
car <- x[1]*Data[,12]/100 + x[2]*Data[,13]/100 + x[4]*(Data[, 4]>=2)
##効用の指数化
Erail <- exp(rail)*Data[, 5]
Ebus <- exp(bus )*Data[, 8]
Ecar <- exp(car )*Data[,11]
NLモデルの計算結果(ehime)
ケースA：rohbar= 0.2652 (最も大きい)
パラメータ：
t 値：
-0.8978 0.01277 1.178 0.5003 0.03899 -0.43394 4.873
-5.62
1.41
8.43 6.39 0.42
-3.40
7.04
時間 費用 高齢B 台数C 鉄道 バス スケール
ケースB：rohbar= 0.2416
パラメータ：
t 値：
-1.327 -0.1613 1.244 2.413 1.877 2.382 1.354
-4.86
-6.01 7.20 14.97 11.86 12.29 9.82
時間 費用 高齢B 台数C 鉄道 バス スケール
ケースC：rohbar= 0.2470
パラメータ：
t 値：
-2.357 -0.1976 2.297 3.008 1.992 2.814 0.5830
-5.69 -4.70 11.29 15.79 9.53 9.81 8.72
時間
費用 高齢B 台数C 鉄道 バス スケール
スケー
ルパラ
メータが
1以上
© Copyright 2024 ExpyDoc
