計量行動分析　第 7 回

計量行動分析　第7回 離散選択モデルロジットモデル

連続的選択と離散的選択 • 標準的消費者行動モデル(連続的選択) • 一定期間に購入する　財・サービスの量を説明　（お金・時間をいくらずつ割り振るか？） • 離散的（質的）選択モデル • どの種類の財を選択するか？　（どこへ行くか？何を買うか？） • 離散的(discrete)な　　選択肢(alternative)からの選択 5万円

離散的選択のモデル化 • 個人は，採りうる選択肢alternativeを列挙する • それぞれの選択肢の特徴と費用に基づいて，評価得点utilityをつける • 評価点が高いものを選ぶフランス旅行　40点アメリカ旅行50点中国旅行　60点

確定的選択モデル • 個人は「評価得点が少しでも高いほうの選択肢を必ず選ぶ」と考えるモデル • 全員が同じ考え方で評価 • 同じ状況に直面する人は，全員同じ選択肢を選択 • 異なる条件の人の選択結果から効用関数を推定 • 判別関数モデル，数量化理論II類モデル • 同じ状況に対する個人間での考え方の違いを考慮 • 例：高齢者は着席可能性を，若者は低運賃を重視 • 犠牲量(最小化）モデル • 比較における「あいまいさ」を認める • ファジィ選択モデル

確率的選択：評価点の差と選択率 • 実際には • ほとんど評価点が同じときは，どちらも選択される可能性がある • 評価点の差が大きいときは，片方しか選ばれない．選択肢Aが選ばれる可能性 1 Aが圧倒的に良いほとんどAだけが選ばれる 2つは同じ魅力 50％ずつ Aが圧倒的に劣る Aが選ばれることはほとんどない 0 選択肢Aの得点－選択肢Bの得点

ロジットモデル • S字型の曲線として，　　という式で表わされる曲線を使うと， • いろいろな計算が簡単にできる • 3つ以上の選択肢からの選択も同じ形になる • 2000年ノーベル経済学賞　 McFadden（1937-）が提案 • 各自の評価点が安定している部分と確率的に変動する部分の和である場合の選択から理論的に導いた。（ランダム効用モデル）

ロジットモデル • 選択肢が2つの2項ロジットモデル Binary （Binomial) Logit • 選択肢が多数（ｎ個）ある場合の　多項ロジットモデル Multinomial Logit

モデルの使用手順 選択モデルの定式化（得点からSカーブで選択率を計算）パラメータη，βの推定（実際の選択結果から、どういう得点付けだったかを推測）将来の選択率の予測（将来の選択肢の状況から得点を計算し、選択率を計算）

実際の観測結果から、もとの事象の発生確率を推定する実際の観測結果から、もとの事象の発生確率を推定する • 手品師がイカサマかも知れないコインを6回振ったら、表が5回、裏が1回出た。 • このコインは正しいコインだろうか？表が出やすいようなイカサマのコインだろうか？ • 表5回、裏1回という観測結果から、このコインを1回振ったときに表が出る確率pの値を推測したい。 • 比率を用いてそのままｐ＝5/6と推測する方法 • さいゆうほう（最尤法）

比率によるロジットモデルの推定 • 集計的方法(集団の選択率にあわせる) • ある選択肢の状況下で観測された集団の選択確率ｐを用いる • ロジットモデルから、その時の2つの選択肢の魅力度VAとVBの差が逆算できる • 魅力度の差がうまく一致するように、魅力度の関数の形を調整する先の例：6回のうち5回表だから、p=5/6と推測した。

回帰分析 Linear Regression 平面をうまく決める 2つ(以上)の変数を用いて、目的の変数yの値をうまく予測できるような平面を決める方法

最尤法による推定(個々の事象の組み合わせを考える）最尤法による推定(個々の事象の組み合わせを考える） • コインを続けて6回振ったところ、そのうち5回が表であった．このコインの表の出る確率ｑはいくらか？ • 母数（ここでは表の出る確率ｑ)を変えたときに，観察された事象が実現する確率を求める(尤度関数L(ｑ)） • 観察された事象の発生確率が最大になるｑの値を求める • 尤度関数を求めてみよう • 確率ｑの事象が5回，（１－ｑ）の事象が1回観察されたのだから、何回目が裏かが6通りあることを踏まえると、

最尤法による母数の推定の例 尤度関数を最大化するあるいは対数をとったものをｑについて最大化してもよいこれは，6回のうち5回という割合に一致

最尤法によるロジットモデルの推定 • 非集計的方法(各個人の選択を考える) • 魅力度の関数形を仮定する • 各個人が直面した選択肢の状況をもとに、ロジットモデルで各自の選択確率を求める。 • それらを掛け合わせて、観測された事象の発生確率(尤度関数)を求める • 尤度関数の値が最大になるように、魅力度の関数の形を調整する

交通手段選択モデル • あるODを移動する消費者が,複数の交通手段の所要時間や費用を考えて手段を選択 • 2項ロジットモデル選択肢jの魅力度が他の選択肢よりも高い確率

2項ロジットモデルの図解 V* 1 V 0 X2 X1

集計データ(比率)を用いたロジットモデルの推定集計データ(比率)を用いたロジットモデルの推定

集計データ(比率)を用いた線形回帰分析による推定集計データ(比率)を用いた線形回帰分析による推定ロジットモデル式から、2つの選択率の比は以下のようになる

集計データ(比率)を用いた線形回帰分析による推定集計データ(比率)を用いた線形回帰分析による推定バスの分担率 VBUS-VCAR

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

最尤法によるロジットモデルの推定

最尤法によるロジットモデルの推定 尤度関数最大になるように値を調整

最尤法によるロジットモデルの推定 ソルバー機能：表計算ソフト：Excel　の中では、いくつかの数値が別のセルの数値に影響を与える場合、目的のセルの数値を最大にするように、元の数値を決定する計算ができる。実際には、この尤度関数は、あまりに値が小さいため、数値計算誤差の影響でうまく計算できない。尤度関数の対数(log)を取ったものを考え、それを最大化する。

最尤法によるロジットモデルの推定 最大になるように値を調整

作成されたロジットモデル

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) Null deviance: 0.02932541 on 8 degrees of freedom Residual deviance: 0.00060564 on 6 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 4

各個人の選択を説明する場合Rにより，個人ごとのデータを作成各個人の選択を説明する場合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

理論的補足ランダム項がある場合の選択

ランダム項の確率分布の仮定 • 集団全体における各自の効用の誤差εの頻度分布を連続的な確率密度関数で表現

選択率の導出

V1-V2 （＞０） V1-V2（＜０）の場合選択率の導出(図解) 選択肢１が選ばれる確率 Prob(U1>U2) =Prob(V1+ε1>V2＋ε２）＝Prob(ε2<V1-V2+ε1) ε2 V1-V2+ε1 ε1 誤差項をどのような分布で考えるか？ (まんじゅうのふくらみ方) ε2 釣鐘（つりがね）まんじゅうを，中央からずれた位置で包丁で切る右下の方の切れ端の体積が選択率を表わす． ε1

具体的なランダム効用モデル • ロジット（LOGIT)モデル • 各選択肢の誤差項が独立で同一の Gumbel分布に従う　と仮定したモデル • 第k選択肢の選択確率は， • プロビット（PROBIT)モデル • 誤差項が多変量正規分布に従うと仮定したモデル • 選択確率を解析的に表示することができない • モンテカルロシミュレーションなどを用い近似値を計算 η：誤差項の分散に　　反比例するパラメータ

非集計データによるロジットモデルのパラメータ推定非集計データによるロジットモデルのパラメータ推定尤度関数尤度関数の値を最大にするようにβを定める

パラメータ推定（２）

ニュートン・ラプソン法

最尤法の計算例 通勤に2 種類の交通手段(1：バス、2：鉄道) があり、10 人の両交通機関での所要時間を調べたところ、表のようであった。ロジットモデルを適用して、所要時間のパラメータとバスの選択し定数項のパラメータを計算せよ。さらに、バス専用レーンの設置により、全てのサンプルのパス所要時間が1 分短縮された場合に、バス選択確率がどのように変化するかを推測せよ。

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

尤度関数の数値計算例

尤度関数の数値計算例(ソルバー) =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)) )

得られたパラメータを用いた予測 現況再現政策実施後

多項ロジットモデルMultinomial Logit Model

モデル作成の手順

多項ロジットモデルMultinomial Logit Model 個人によって全ての選択肢が利用できない場合もありうる

推定結果の評価

計量行動分析　第 7 回