Data Analysis Using Statistical Packages: Discrete Regression Analysis

このページに対応するRmdファイル：GitHub

このページでは次を扱う：

• 二値回帰モデル：アウトカムが二値 {0, 1} の場合のプロビットモデルとロジットモデル．
• ランダム効用モデルに基づく離散選択モデル：2つ以上のカテゴリー（離散的なアウトカム）から1つを選択する行動モデル．ロジットモデルのみ扱う．

1 Binary regression model / 二値回帰モデル

Titanic データの読み込み．

titanic <- read.csv("https://raw.githubusercontent.com/kurodaecon/bs/main/data/titanic3_csv.csv")

library(tidyverse)

1.1 Linear probability model / 線形確率モデル

生存したかどうかを表す Survived 変数を，性別，年齢，客室等級変数に回帰．

\[ \mbox{Survived} = \beta_0 + \beta_1 \mbox{Male} + \beta_2 \mbox{Age} + \beta_3 \mbox{2nd-class} + \beta_4 \mbox{3rd-class} \]

summary(lm(survived ~ sex + age + factor(pclass), data = titanic))

## 
## Call:
## lm(formula = survived ~ sex + age + factor(pclass), data = titanic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09442 -0.24408 -0.08375  0.23289  0.99397 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.1049549  0.0438211  25.215  < 2e-16 ***
## sexmale         -0.4914132  0.0255520 -19.232  < 2e-16 ***
## age             -0.0052695  0.0009316  -5.656 2.00e-08 ***
## factor(pclass)2 -0.2113738  0.0348568  -6.064 1.85e-09 ***
## factor(pclass)3 -0.3703874  0.0325039 -11.395  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3913 on 1041 degrees of freedom
##   (263 observations deleted due to missingness)
## Multiple R-squared:  0.3691, Adjusted R-squared:  0.3667 
## F-statistic: 152.3 on 4 and 1041 DF,  p-value: < 2.2e-16

解釈

• 男性の生存率は49%ポイント低い
• 年齢が1歳高い乗客の生存率は0.5%ポイント低い
• 2等客室の乗客の生存率は1等と比べて21%ポイント低い

1.2 Binary probit model / 二値プロビット・モデル

\[ \Pr (Y = 1 \mid X) = \Phi (\beta_0 + \beta_1 X) \]

glm 関数で推定．

probit_titanic <- glm(formula = survived ~ sex + age + factor(pclass), data = titanic, 
                      family = binomial(link = "probit"))
summary(probit_titanic)

## 
## Call:
## glm(formula = survived ~ sex + age + factor(pclass), family = binomial(link = "probit"), 
##     data = titanic)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.055101   0.181290  11.336  < 2e-16 ***
## sexmale         -1.485640   0.093813 -15.836  < 2e-16 ***
## age             -0.019426   0.003606  -5.387 7.16e-08 ***
## factor(pclass)2 -0.760170   0.129921  -5.851 4.89e-09 ***
## factor(pclass)3 -1.303160   0.126633 -10.291  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1414.62  on 1045  degrees of freedom
## Residual deviance:  984.55  on 1041  degrees of freedom
##   (263 observations deleted due to missingness)
## AIC: 994.55
## 
## Number of Fisher Scoring iterations: 5

1.2.1 Estimation by numerical calculation / 数値計算による推定

対数尤度関数を汎用の最適化関数 optim で最大化することで推定する．

選択確率は以下で与えられる．

\[ \Pr (Y=1 \mid X) = \Phi (\beta_0 + \beta_1 X) = \frac{1}{2 \pi} \int_{-\infty}^{\beta_0 + \beta_1 X} \exp \left( - \frac{t^2}{2} \right) dt \]

よって，尤度関数は以下のようになる．

\[ L (\beta_0, \beta_1) = \prod_i \left[ \{ \Pr (Y_i = 1 \mid X_i) \}^{Y_i} \times \{ \Pr (Y_i = 0 \mid X_i) \}^{1 - Y_i} \right] \] \[ = \prod_i \left[ \{ \Phi (\beta_0 + \beta_1 X_i) \}^{Y_i} \times \{ 1 - \Phi (\beta_0 + \beta_1 X_i) \}^{1 - Y_i} \right] \]

対数を取ると対数尤度関数が得られる．

\[ \log L (\beta_0, \beta_1) = \sum_i \left[ Y_i \times \log \{ \Phi (\beta_0 + \beta_1 X_i) \} + (1 - Y_i) \times \log \{ 1 - \Phi (\beta_0 + \beta_1 X_i) \} \right] \]

これを，\((\beta_0, \beta_1)\) を引数とする関数として定義すればよい．

• 標準正規分布の分布関数は pnorm で与えられる．
• 上の式では一つの説明変数を用いて \(\beta_0 + \beta_1 X\) と書いているが，一般には複数の説明変数を用いて \(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_K X_K = \mathbf{X} \boldsymbol{\beta}\) と書ける．

ml_probit <- function (y_probit, x_probit) {
  LL_probit <- function (param) {
    pnorm_xbeta <- pnorm(q = x_probit %*% param, mean = 0, sd = 1)
    sum(y_probit %*% log(pnorm_xbeta) + (1-y_probit) %*% log(1 - pnorm_xbeta))
  }
  optim(par = rep(0, ncol(x_probit)), fn = LL_probit, control = list(fnscale = -1))
}

上と同じ説明変数を使って推定する．

まずは，計算のために character/factor 型に対応する数値を計算しておく．

titanic_numeric <- titanic %>% 
  filter(!is.na(age)) %>% 
  mutate(male = 1*(sex == "male"),
         pclass2 = 1*(pclass == 2),
         pclass3 = 1*(pclass == 3))

titanic_x <- titanic_numeric %>% 
  dplyr::select(male, age, pclass2, pclass3) %>% 
  as.matrix
ml_probit(y_probit = titanic_numeric$survived, x_probit = cbind(1, titanic_x))

## $par
## [1]  2.05418759 -1.48531017 -0.01941115 -0.76010160 -1.30246558
## 
## $value
## [1] -492.2765
## 
## $counts
## function gradient 
##      468       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

1.2.2 Marginal effect / 限界効果

定義

\[ \frac{d \Pr(Y_i = 1 \mid X_i)}{d X_i} = \frac{d \Phi(\beta_0 + \beta_1 X_i)}{d X_i} = \phi(\beta_0 + \beta_1 X_i) \cdot \beta_1 \]

説明変数の平均値 \(\bar{X}\) を用いる方法．

\[ \phi (\hat{\beta_0} + \hat{\beta_1} \bar{X}) \hat{\beta_1} \]

xb_mean <- probit_titanic$coef %*% c(1, mean(titanic_numeric$male), mean(titanic$age, na.rm = TRUE), 
                                     mean(titanic$pclass == 2), mean(titanic$pclass == 3))
as.numeric(dnorm(x = xb_mean)) * probit_titanic$coef

##     (Intercept)         sexmale             age factor(pclass)2 factor(pclass)3 
##      0.77727888     -0.56189779     -0.00734726     -0.28751081     -0.49288013

Individual ごとの限界効果の平均値を計算する方法．

\[ \frac{1}{n} \sum_i \phi (\hat{\beta_0} + \hat{\beta}_1 X_i) \hat{\beta}_1 \]

xb <- as.matrix(cbind(1, titanic_numeric[, c("male", "age", "pclass2", "pclass3")])) %*% probit_titanic$coef
c(male = mean(dnorm(x = xb) * probit_titanic$coef[2], na.rm = TRUE), 
  age = mean(dnorm(x = xb) * probit_titanic$coef[3], na.rm = TRUE), 
  class2 = mean(dnorm(x = xb) * probit_titanic$coef[4], na.rm = TRUE), 
  class3 = mean(dnorm(x = xb) * probit_titanic$coef[5], na.rm = TRUE))

##         male          age       class2       class3 
## -0.392874167 -0.005137142 -0.201025116 -0.344617610

mfx パッケージを使う場合．

# install.packages("mfx")
library(mfx)
probitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)

## Call:
## probitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)
## 
## Marginal Effects:
##                      dF/dx  Std. Err.        z     P>|z|    
## sexmale         -0.5408979  0.0287920 -18.7864 < 2.2e-16 ***
## age             -0.0074648  0.0013850  -5.3898 7.052e-08 ***
## factor(pclass)2 -0.2672290  0.0404958  -6.5989 4.141e-11 ***
## factor(pclass)3 -0.4666238  0.0395420 -11.8007 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "sexmale"         "factor(pclass)2" "factor(pclass)3"

1.3 Binary logit model / 二値ロジット・モデル

\[ \Pr (Y = 1) = \Lambda (\beta_0 + \beta_1 X) \]

glm 関数で推定．

logit_titanic <- glm(formula = survived ~ sex + age + factor(pclass), data = titanic, 
                     family = binomial(link = "logit"))
summary(logit_titanic)

## 
## Call:
## glm(formula = survived ~ sex + age + factor(pclass), family = binomial(link = "logit"), 
##     data = titanic)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.522075   0.326702  10.781  < 2e-16 ***
## sexmale         -2.497845   0.166037 -15.044  < 2e-16 ***
## age             -0.034393   0.006331  -5.433 5.56e-08 ***
## factor(pclass)2 -1.280571   0.225538  -5.678 1.36e-08 ***
## factor(pclass)3 -2.289661   0.225802 -10.140  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1414.62  on 1045  degrees of freedom
## Residual deviance:  982.45  on 1041  degrees of freedom
##   (263 observations deleted due to missingness)
## AIC: 992.45
## 
## Number of Fisher Scoring iterations: 4

余裕がある履修者向けの宿題：プロビットと同様に対数尤度関数を書いて optim で最適化計算する．

1.3.1 Marginal effect / 限界効果

logitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)

## Call:
## logitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)
## 
## Marginal Effects:
##                      dF/dx  Std. Err.        z     P>|z|    
## sexmale         -0.5514348  0.0293172 -18.8092 < 2.2e-16 ***
## age             -0.0080960  0.0014881  -5.4405 5.312e-08 ***
## factor(pclass)2 -0.2673472  0.0404682  -6.6064 3.939e-11 ***
## factor(pclass)3 -0.4901785  0.0400797 -12.2301 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "sexmale"         "factor(pclass)2" "factor(pclass)3"

余裕がある履修者向けの宿題：プロビットと同様に（mfx::logitmfx 関数を使わずに）計算する．

2 Discrete choice based on random utility model / ランダム効用モデルに基づく離散選択

経済主体（消費者など）の効用最大化を前提とすれば，選択肢が複数（3つ以上でもよい）あるときにそのうち一つを選択するということは，選択した選択肢から得られる効用（≒ 嬉しさ・満足度）が他の選択肢から得られる効用と比較してより大きいことを意味するはず．

選択肢 \(j (\in \{ 1,\ldots,J \} )\) を選択・消費することから得る間接効用を \(u_j\) で表すとき，実際に選択される選択肢 \(j^*\) は次のように決まる．

\[ j^* = \arg \max_{j} u_j \]

\(u_j\) を分析者にとって観察できる要因によって説明できる部分 \(V_j\) と観察できない要因による部分 \(\epsilon_j\) に分解する．

\[ u_j = V_j + \epsilon_j \]

これをランダム効用モデルと呼ぶ．

離散選択モデルの選択確率はこのランダム効用モデルに基づいて次のように記述できる．

\[ \Pr (y_j = 1) = \Pr (u_j \ge u_k, \forall k \ne j) \\ = \Pr (V_j + \epsilon_j \ge V_k + \epsilon_k, \forall k \ne j) = \Pr (\epsilon_j - \epsilon_k \ge - V_j + V_k, \forall k \ne j) \]

ここで \(y_j\) は選択肢 \(j\) を選択した場合に1をとる indicator である．

\(\epsilon_j\) に第一種極値分布を仮定すると（このとき，\(\epsilon_j - \epsilon_k\) がロジスティック分布に従う），この確率はロジット選択確率で表される（証明は Train (2009) §3.10 などを参照）．

\[ \Pr (y_j = 1) = \frac{\exp(V_j)}{\sum_k \exp(V_k)} \]

推定に用いるデータの形式によって \(V_{ij}\) は異なる関数形を取る（\(i\): 消費者など）．

• たとえば，\(x\) が \(j\) 間で変分を持つときは \(V_{ij} = \sum_{m} \beta_{m} x_{jm}\) のように定式化できる．このとき \(\beta\) は \(j\) 間で共通のものが推定される．これは conditional logit と呼ばれる．\(x_{ijm}\) の場合もこの枠組みで推定できる．
• \(x\) が \(j\) 間で変分を持たないとき，すなわち \(x\) は \(i\) 間でのみ変分を持つとき，\(V_{ij} = \sum_{m} \beta_{jm} x_{im}\) のように定式化される．このとき \(\beta\) は \(j\) ごとに推定される．これは multinomial logit と呼ばれる．

2.1 Data for simulation

間接効用関数は次のように設定する（\(\beta = -2, j = 0\): outside option）：

\[ u_j = - 2 p_j + \epsilon_j, \\ u_0 = \epsilon_0 \]

# install.packages("evd")
library(evd)  # 第一種極値分布からの乱数ドローに必要
N <- 1000  # i = 1, ..., 1000
J <- 2  # j = 0 (outside option = not purchase), 1, 2
param_beta <- -2
set.seed(0)
data <- tibble::tibble(
  id = rep(1:N, each = 1+J),
  j = rep(0:J, times = N),
  p = rep(c(0, 0.2, 0.3), times = N),
  e = evd::rgumbel(N * (1+J)),
  u = param_beta * p + e,
  y = 0
)
for (i in 1:N) {  # j* = arg max u_j
  index_i <- data$id == i
  index_max_ui <- which.max(data$u[index_i])[1]
  data$y[index_i][index_max_ui] <- 1
}
data %>% group_by(j) %>% summarise(n_chosen = sum(y))

## # A tibble: 3 × 2
##       j n_chosen
##   <int>    <dbl>
## 1     0      427
## 2     1      313
## 3     2      260

2.2 ML estimation / 最尤推定

2.2.1 survival::clogit function

# install.packages("survival")
library(survival)
cl0 <- survival::clogit(y ~ p + strata(id), data)
cl0

## Call:
## survival::clogit(y ~ p + strata(id), data)
## 
##      coef exp(coef) se(coef)      z        p
## p -1.6307    0.1958   0.2482 -6.569 5.05e-11
## 
## Likelihood ratio test=42.81  on 1 df, p=6.043e-11
## n= 3000, number of events= 1000

logLik(cl0)  # maximized (ln L)

## 'log Lik.' -1077.209 (df=1)

2.2.2 mlogit::mlogit function

# install.packages("mlogit")
library(mlogit)
mlogit_data <- mlogit.data(data, shape = "long", choice = "y", chid.var = "id", alt.var = "j")
summary(mlogit(y ~ p | 0, data = mlogit_data))

## 
## Call:
## mlogit(formula = y ~ p | 0, data = mlogit_data, method = "nr")
## 
## Frequencies of alternatives:choice
##     0     1     2 
## 0.427 0.313 0.260 
## 
## nr method
## 3 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.0329 
## successive function values within tolerance limits 
## 
## Coefficients :
##   Estimate Std. Error z-value  Pr(>|z|)    
## p -1.63068    0.24823 -6.5693 5.054e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1077.2

2.2.3 Handmade likelihood function + general-purpose optimizer

尤度関数と対数尤度関数

\[ L (\beta) = \prod_i \prod_j \Pr(y_{ij}=1)^{y_{ij}}, \\ \ln L (\beta) = \sum_i \sum_j \left[ y_{ij} \times \ln \Pr(y_{ij}=1) \right] \]

• exp の計算で数値計算上の問題が生じにくいように次のような工夫をする（\(\max(V) = \max \{ V_0, V_1, V_2 \}\)）：

\[ \frac{\exp(V_j)}{\sum_k \exp(V_k)} = \frac{\exp(V_j) \times \exp(-\max(V))}{\sum_k \exp(V_k) \times \exp(-\max(V))} = \frac{\exp(V_j - \max(V))}{\sum_k \exp(V_k - \max(V))} \]

lnL <- function (param, d) { # 対数尤度関数
  d <- d %>% 
    mutate(v = param * p) %>% 
    group_by(id) %>% 
    mutate(max_v = max(v),  # max of {V_i0, V_i1, V_i2}
           v_minus_max = exp(v - max_v),  # exp V-max
           sum_v = sum(v_minus_max)) %>%  # sum(exp V-max)
    ungroup() %>% 
    mutate(prob = v_minus_max / sum_v,  # exp V-max / sum(exp V-max)
           y_lnprob = y * log(prob))
  sum(d$y_lnprob)
}

2.2.3.1 optim function

デフォルトでは目的関数を最小化するので control = list(fnscale = -1) で最大化問題として指定．

opt_logit <- optim(par = -1, fn = lnL, d = data,
 control = list(fnscale = -1), method = "BFGS", hessian = TRUE)
opt_logit

## $par
## [1] -1.630674
## 
## $value
## [1] -1077.209
## 
## $counts
## function gradient 
##       19        3 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
##           [,1]
## [1,] -16.22952

• $par: 推定されたパラメタ \(\hat{\beta}\)
• $value: 最大化された対数尤度関数の値
• $hessian: Hessian = \(\nabla^2 \ln L(\beta)\) （対数尤度関数の二階微分）

推定されたパラメタの標準誤差は Hessian を利用して次のように計算できる．

\[ \mbox{SE} (\beta) = \sqrt{ \mbox{diag} \left( I^{-1} \right) } \quad \mbox{where} \quad I = - \nabla^2 \ln L (\beta) \]

• \(I\): フィッシャー情報行列

opt_se <- - opt_logit$hessian |> solve() |> diag() |> sqrt()  # std. error
opt_se

## [1] 0.2482259

opt_logit$par / opt_se  # z value 

## [1] -6.569315

2 * (1 - pnorm(q = abs(opt_logit$par / opt_se)))  # p value 

## [1] 5.054734e-11

2.2.3.2 Note: Counterfactual simulation

元の変数（価格変数のつもり）は \(p = (0, 0.2, 0.3)\) だが，これを \(p = (0, 0.1, 0.3)\) に変更したときに，各選択肢の選択確率がどのように変化するかを計算してみる．

このシミュレーションは，「元の価格変数（変数名の suffix：_act）」と「反実仮想の価格変数（変数名の suffix：_cf）」それぞれについて次を計算することによって実行される．

• \(V\) を計算
• \(i\) ごとに \(\max V, \exp(V_j - \max V), \sum_k (\exp(V_k - \max V))\) を計算
• \(\Pr(y_{ij} = 1)\) を計算
• \(\Pr(y_{ij} = 1)\) の平均値を \(j\) ごとに計算（市場シェアとして解釈可能）

data %>%
  mutate(
    v_act = opt_logit$par * p,  # act = actual 
    p_cf = rep(c(0, 0.1, 0.3), times = N),  # cf = counterfactual 
    v_cf = opt_logit$par * p_cf
  ) %>% 
  group_by(id) %>% 
  mutate(
    # actual 
    max_v_act = max(v_act),
    v_minus_max_act = exp(v_act - max_v_act),
    sum_v_act = sum(v_minus_max_act),
    # counterfactual 
    max_v_cf = max(v_cf),
    v_minus_max_cf = exp(v_cf - max_v_cf),
    sum_v_cf = sum(v_minus_max_cf)
  ) %>% 
  ungroup() %>% 
  mutate(
    prob_act = v_minus_max_act / sum_v_act,
    prob_cf = v_minus_max_cf / sum_v_cf
  ) %>% 
  group_by(j) %>%
  summarise(
    prob_act = mean(prob_act),
    prob_cf = mean(prob_cf)
  )

## # A tibble: 3 × 3
##       j prob_act prob_cf
##   <int>    <dbl>   <dbl>
## 1     0    0.428   0.406
## 2     1    0.309   0.345
## 3     2    0.263   0.249

価格が下がった \(j=1\) の選択確率が上昇し，その分 \(j=0,2\) の選択確率が低下したことが観察できる．

注：

• 本分析では変数 \(p\) の値は \(j\) ごとにすべての \(i\) 間で一律の値を設定しているので，実際には上をすべての \(i\) について計算する必要はなく，特定の \(i\) について計算した選択確率がそのまま市場シェアとして解釈できる．
• 選択確率の平均値の比較ではなく \(\epsilon\) を極値分布からドローして \(j^*\) を求めて比較することも可能（だが，計算結果は \(\epsilon\) のドローによる不確実性を含む）．

2.2.3.3 nlm function

こちらも同様にデフォルトでは目的関数を最小化するので，目的関数を \(- \ln L\) に置き換えてこの関数の最小化問題とする．

注：

• $minimum は \(- \ln L\) の最小値なので，最大化された \(\ln L\) の値は $minimum \(\times (-1)\) となる．
• $hessian は計算上の目的関数 \(- \ln L\) に対応するものなので，本来の目的関数 \(\ln L\) に対応する Hessian は $hessian \(\times (-1)\) となる．

opt_logit_nl <- nlm(f = function (param) -1 * lnL(param, d = data), 
                    p = -1, hessian = TRUE)
opt_logit_nl

## $minimum
## [1] 1077.209
## 
## $estimate
## [1] -1.630676
## 
## $gradient
## [1] 6.971763e-07
## 
## $hessian
##          [,1]
## [1,] 16.22948
## 
## $code
## [1] 1
## 
## $iterations
## [1] 3

opt_logit_nl$hessian |> solve() |> diag() |> sqrt()  # std. error

## [1] 0.2482263

2.3 Bayesian estimation: Metropolis

iter <- 1000
beta_old <- beta_tmp <- 0  # 初期値
beta_history <- numeric(iter)
ll_old <- lnL(param = beta_tmp, d = data)
accept <- numeric(iter)  # initialize (=1 if accepted) 
start_time <- Sys.time()  # 計算時間を確認するため
for (i in 1:iter) {
  beta_tmp <- rnorm(n = 1, mean = beta_tmp, sd = 0.5)
  ll_tmp <- lnL(param = beta_tmp, d = data)
  accept_prob <- min(0, ll_tmp - ll_old)
  if (log(runif(n = 1)) < accept_prob) {
    beta_old <- beta_tmp
    ll_old <- ll_tmp
    accept[i] <- 1
  } else {
    beta_tmp <- beta_old
  }
  beta_history[i] <- beta_tmp
  if (i %% round(iter/2, 0) == 0) {
    elapsed_sec <- round(as.numeric(Sys.time() - start_time, units = "secs"), 1)
    cat("iter=", i, ", lnL=", ll_old, ", beta=", beta_tmp, 
        ", elapsed=", elapsed_sec,"\n", sep = "")
  }
}

## iter=500, lnL=-1077.48, beta=-1.447759, elapsed=3.7
## iter=1000, lnL=-1077.472, beta=-1.450626, elapsed=7.4

mean(accept)

## [1] 0.485

plot(beta_history, type = "l")
abline(h = opt_logit$par, col = 2)  # ML estimate 

Burn-in は最初の20%とする．

beta_after_burnin <- beta_history[round(iter/5, 0):iter]  # 間引かない
hist(beta_after_burnin, breaks = 20)
abline(v = opt_logit$par, col = 2)

c(mean(beta_after_burnin), sd(beta_after_burnin))

## [1] -1.6278531  0.2288725

3 Advanced: Conditional logit

興味がある履修者は 自動車保有の離散選択問題：Conditional logit を参照いただきたい．

ランダム効用モデルに基づいて定式化すれば選択者サイドの厚生 (welfare) が計算できることも分かる．

（注：残念ながらコードがところどころ間違っているので，雰囲気だけ感じてください．）