機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan：単回帰分析

Stan

はじめに


統計モデリングの勉強の一環に「Stan」を触ってみる。
書籍は「アヒル本」（StnaとRでベイズ統計モデリング）を参考とする。

初回は、統計モデリングの「Hello world」と言われる単回帰分析を扱う。（アヒル本のChapter.4）

前処理

データは、社員の年齢と年収のデータ。
（Xが年齢、Yが年収）

> data <- read.csv("data-salary.txt")
> head(data)
   X   Y
1 24 472
2 24 403
3 26 454
4 32 575
5 33 546
6 35 781

プロットするとこんな感じ。

> plot(data$X, data$Y)

f:id:kento1109:20180412121652p:plain

lm関数で推定

まずは、Stanを利用せずにlm関数で最尤推定を行う。

> res_lm <- lm(Y~X, data)
> res_lm

Call:
lm(formula = Y ~ X, data = data)

Coefficients:
(Intercept)            X  
     -119.7         21.9

描画するとこんな感じ。
f:id:kento1109:20180412153607p:plain

シンプルな単回帰分析では、これで十分だと思う。

Stanでモデリング

今回はこのデータをモデリングする。

今回の単回帰分析は以下のように書ける。
Y=a+Xb
これだと、全ての点が直線上にのることになる。
実際はこの直線を中心に各点ごとにノイズが発生しているので、
Y=a+Xb+\varepsilon
となる。
ノイズ\varepsilonは、以下の分布に従って発生するとする。
\varepsilon\sim{\rm Normal}(0,\sigma)
\varepsilonの平均は、ノイズが全くない状態なので、
aX+bである。
つまり、２つの式をまとめて、
Y\sim{\rm Normal}(a+Xb,\sigma)
と書ける。

これをStanの文法に沿って書いてみる。

data {
  int N;
  real X[N];
  real Y[N];
}
parameters {
  real a;
  real b;
  real<lower=0> sigma;
}
model {
  for (n in 1:N){
    Y[n] ~ normal(a + X[n] * b, sigma);
  }
}

簡単であるが、各ブロックを説明する。

  • dataブロック

　実際に使うデータを格納するための変数と型を宣言する。

  • parametersブロック

　推定したいパラメータを宣言する。

  • modelブロック

　モデリングのメインブロック。
　dataブロック 、parametersブロックで宣言した変数を用いてモデルを記述する。

Stanの実行

実行させるには、上記のStanコードを読み込んで実行する。
方法としては、

  1. 直接コードを変数に格納する
  2. Stanファイルにコードを記述して呼び出す

の２種類がある。
簡単なコードの場合、直接コードを変数に格納する方が楽。

1.の場合、下記のようにして実行する。

> my_code = '
+ data {
+     int N;
+     real X[N];
+     real Y[N];
+ }
+ parameters {
+     real a;
+     real b;
+     real<lower=0> sigma;
+ }
+ model {
+     for (n in 1:N){
+         Y[n] ~ normal(a + X[n] * b, sigma);
+     }
+ }
+ '

> model <- stan(model_code = my_code,
+ data = list(N=nrow(data), X=data$X, Y=data$Y))

stan関数の引数である、dataでは、dataブロックで定義した変数に格納するデータを渡す。

2.のファイルを読み込む場合は下記の通り。

> model <- stan(file='model.stan',
+ data = list(N=nrow(data), X=data$X, Y=data$Y))
実行結果

model に格納された結果を確認する。

> model
Inference for Stan model: be6814ec3b164fbfad2c7ec7333e3157.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean    sd    2.5%     25%     50%    75%  97.5% n_eff Rhat
a     -117.47    1.90 72.10 -257.79 -165.45 -116.53 -70.70  23.21  1443    1
b       21.84    0.04  1.61   18.69   20.79   21.87  22.89  25.00  1466    1
sigma   85.27    0.38 15.08   61.98   74.15   83.29  94.10 120.76  1612    1
lp__   -93.61    0.03  1.21  -96.60  -94.20  -93.33 -92.72 -92.16  1481    1

Samples were drawn using NUTS(diag_e) at Thu Apr 12 16:14:09 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

１列目はMCMCサンプルの平均値。
（最尤推定の値と近似していることが分かる。）

次に

  • トレースの確認
> stan_trace(model)

f:id:kento1109:20180412165301p:plain

stan_trace(model,pars = "a")

で変数を選択できる。
f:id:kento1109:20180412165451p:plain

stan_hist(model)

f:id:kento1109:20180412165642p:plain
これも同様に変数選択できる。

stan_hist(model,pars ="a")

f:id:kento1109:20180412165656p:plain

stan_dens(model,pars ="a")

f:id:kento1109:20180412165833p:plain

  • 自己相関
stan_ac(model, pars="a")

f:id:kento1109:20180412170339p:plain

収束の判断

目安の１つとして下記の指標が紹介ある。

  • Rhat：全ての変数が「1.05」以下
  • 密度関数：全てのチェインが１つの分布に収束しているか
  • 自己相関：１つ前のサンプルと相関が生じていないか

（参考）

www.slideshare.net

※自己相関の基準はよくわかってないので、別の機会に調べる。

データの取り出し

MCMCサンプル結果をR上の変数とする場合、以下で取り出す。

ms <- rstan::extract(model)
head(ms$a)
[1] -191.37256  -81.52679  -32.57347 -106.43938 -155.87620  -38.85424

> length(ms$a)
[1] 4000

> summary(ms$a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -409.9  -165.4  -116.5  -117.5   -70.7   169.2

事前分布の考え方については下記が勉強になる。

www.slideshare.net