ほくそ笑む

R言語と統計解析について

実践 統計モデリング入門 【2. 最も簡単なモデルの作成】

はじめに

本シリーズの概要については、下記記事を最初に読んでください。

それでは、モデルの構築を始めていきます。
問題は次の通りでした。

スーパーで定価の同じ 3 種類の牛乳が並べて売られています。
このとき、消費者が 3 種類のブランドのうちどの牛乳を選ぶのか、その行動原理をモデル化して下さい。

しかし、いきなりこの問題に対するモデリングを始めるわけではありません。
統計モデリングのコツ「シンプルなモデルから始める」を思い出して下さい。

最も簡単なモデルの作成

まずは、この問題の最も簡単な場合として、商品が 1 種類しかない場合の消費者の行動選択モデルを考えます。
データとしては、商品の値引き額(x)と買ったか買わなかったか(y)があるとします。

y  x
0 50
1 49
0 48
1 48
0 47

y = 1 が購入、y = 0 が購入しない、x は商品の値引き額(円)です。

このとき、次のようなモデルが考えられます。

  • 消費者は値引き額が大きいほどお得感を感じる
  • お得感が、ある閾値を超えた場合に購入する
  • 感じるお得感は、同じ値引き額に対してもその時の気分によってばらつきがある

これを数式で表すと、次のようになります。
任意の時点  t において、

  • お得感 U_t = x_t + \varepsilon_t ( x_t は値引き額)
  •  \varepsilon_t \sim Norm(0, \sigma)

ある閾値  th に対して、

  •  y_t=1 ( U_t \geq th の場合)
  •  y_t=0 ( U_t < th の場合)

このモデルにおける、推定したいパラメータは、お得感の閾値 th と、気分によるばらつき  \sigma です。
ここまで表現できれば、あとは Stan を使ったモデル記述に落とし込んでいきます。

Stan によるモデルの記述

Stan では、一般的なモデルであれば sampling statement *1 によってモデルを直接的に記述することができます。
しかし、今回のモデルではそうはいきません。
そこで、increment_log_prob() 関数によって対数尤度を足しこんでいくという方法を取ります。

したがって、まずは対数尤度の式を導出します。
まず、 y_t=0 の場合、尤度は  Norm(x_t, \sigma) に従う確率変数が閾値  th を超えない確率と同じなので、

  •  LogLik_t = \log(cdf_t(th))

となります。ただし、 cdf_t Norm(x_t, \sigma) の累積分布関数です。
同様に、 y_t=1 の場合、尤度は  Norm(x_t, \sigma) に従う確率変数が閾値  th を超える確率と同じなので、

  •  LogLik_t = \log(1 - cdf_t(th))

となります。

Stan には、正規分布の対数累積分布関数 normal_cdf_log() およびその complementary 版である normal_ccdf_log() があるため、これを使えば上記の対数尤度を簡単に記述できます。*2

Stan によるモデルの記述は次のようになります。

// model1.stan
data {
  int<lower=0> N;
  real x[N];
  int<lower=0, upper=1> y[N];
}
parameters {
  real<lower=0> sigma;
  real threshold;
}
model {
  for(i in 1:N) {
    if (y[i] == 0)
      increment_log_prob( normal_cdf_log(threshold, x[i], sigma));
    else
      increment_log_prob(normal_ccdf_log(threshold, x[i], sigma));
  }
  sigma ~ cauchy(0, 5);
  threshold ~ normal(0, 1000);
}

Stan では、data ブロックに入力されるデータ(x, y)、parameters ブロックに推定したいパラメータ( \sigma,  th)、model ブロックにモデルを記述します。

以上により、モデルは完成です。

しかし、ここで終わってはいけません。
統計モデリングのコツ「シミュレーションを行う」に従って、このモデルが理想的なデータに対してパラメータ推定できることを確かめてみましょう。

シミュレーションデータの作成

Stan へ入力するためのシミュレーションデータを作成します。

generate_data <- function(N) {
  sigma <- 5      # parameter
  threshold <- 50 # parameter
  x <- rpois(N, lambda = 50)
  u <- x + rnorm(N, 0, sigma)
  y <- ifelse(u >= threshold, 1, 0)
  data.frame(y=y, x=x)
}

data <- generate_data(N = 1000)
head(data, 5)
  y  x
1 0 50
2 0 56
3 1 51
4 1 59
5 0 50

真のパラメータ  \sigma=5 th = 50 のシミュレーションデータができました。
それでは、このデータに対して先程のモデルを当てはめてみましょう。

モデルの当てはめ

Stan コードにより記述されたモデルをデータに当てはめるための R コードは次のようになります。

library(rstan)

standata <- list(N=nrow(data), x=data$x, y=data$y)     # データの整形
model <- stan_model(file = "model1.stan")              # Stan コードのコンパイル
result <- sampling(model, data = standata, chains = 3) # MCMC サンプリングの実行

# 結果の表示
traceplot(result, inc_warm = FALSE)
print(result, digits=2)

これを実行し、パラメータの推定結果を見てみます。


             mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
sigma        4.84    0.01 0.29    4.29    4.63    4.83    5.02    5.44  1743    1
threshold   50.35    0.01 0.24   49.88   50.18   50.35   50.51   50.81  1968    1
lp__      -394.87    0.03 1.03 -397.67 -395.27 -394.54 -394.13 -393.87   966    1

Rhat < 1.1 であり、パラメータはちゃんと定常分布に収束していることがわかります。
sigma = 5 および threshold = 50 がうまく推定できていることが確認できました。

まとめ

最も簡単なモデルに対して Stan によるモデルの記述を行い、シミュレーションデータに対してうまくパラメータ推定できることを確認しました。

このモデルでは、ブランドを 1 種類に限定し、買うか買わないかの判断基準をモデル化しましたが、次はブランドを 2 種類に増やした場合について考えます。

*1:y ~ normal(mu, sigma) みたいな書き方

*2:ccdf = 1 - cdf