誰？

著書(共著)

著書 2

R ユーザのための RStudio[実践]入門
− tidyverse によるモダンな分析フローの世界−

この発表の目的

• ベイズ実践編の方々(LT, 応用セッション)に繋げる
• 繋げる: 全部は理解できなくても、何をしているのかくらいはだいたい分かるようにする
• ベイズ推定で頻繁に用いられる MCMC の基礎を理解
• 最近人気の Stan の基礎を理解

おことわり

• 数学的な解説は最小限に留めます(教科書に書いてあるから)。
• 「何をやっているのか」を直感的に理解することを優先

このスライドを見るには

https://github.com/ymattu/TokyoR70

• このスライドはrevealjs パッケージで作成しています

なぜ MCMC か

$$ \begin{aligned} p(\alpha \mid x) &= \frac{p(x \mid \alpha) \times p(\alpha)}{p(x)} \end{aligned} $$

• 知りたいのは事後分布p(α ∣ x)
• 例: y = βx + e で βの事後分布を求めたい
• でも p(x)=∫p(x ∣ α)×p(α)dα は一般に計算が大変
  → MCMC をつかう

MCMC(Markov Chain Monte Carlo)

• すごくざっくり:「サンプルをたくさん発生させて事後分布のかわりにしよう」
• ざっくり: データxが与えられたとき、事後分布p(α ∣ x)からパラメータαをサンプリングする手法

細かいことは抜きにして:

• Markov Chain: 前の状態で次が決まる
• Monte Carlo: 確率的なアルゴリズム

細かいことは抜きにして。。。

1. 初期点を決める
2. パラメータをしらみつぶしに変化させて尤度を計算していく
3. パラメータをプラス方向かマイナス方向に変えてみて、尤度が大きくなるようならその方向に続けて変えていく
4. 逆に尤度が小さくなるようなら、反対方向にパラメータを変えていく
5. これを収束するまで繰り返す
6. MCMC サンプルが完成→これをパラメータの分布にしよう!

よく見る「あれ」が出来上がります

• これを traceplot といいます
• なんか 1 回じゃ心もとないから複数回(chain)やったほうがいいよね

ここで疑問

1. なんでこうなるの？
2. 収束ってどういうこと？

仮定

• マルコフ連鎖では任意の確率変数からスタートした時に、最終的な確率変数の分布が不変になることを仮定
• この分布を不変分布という

お天気の例

https://www.slideshare.net/teramonagi/ss-5190440

MCMC のアイディア

• サンプルを得たい確率密度関数が p(x)
• これを不変分布に持つようなマルコフ連鎖を生成(推移確率を与える)することで、実際のサンプリングをしよう
• サンプルが不変分布に近づく十分条件に詳細釣り合い条件というものがある

詳細釣り合い条件

ある状態 x^t から x^t + 1 に移るとき、次の式が成り立つ

p(x_t)q(x_t + 1 ∣ x_t)=p(x_t + 1)q(x_t ∣ x_t + 1)

この q(x ∣ y) は提案分布と呼ばれ、現在の状態 x^t から次の状態 x^t + 1 を得るために用いる確率分布

MCMC アルゴリズムの種類

• メトロポリス・ヘイスティングス法(M-H 法)
• ギブスサンプリング

アイディア

• 必ずしも詳細つり合い条件を満たさない提案分布q(x ∣ y)に対して、詳細釣り合い条件を成り立たせるための重み関数$a(x y) $が満たす条件を考える

p(x^tt)q(x^t + 1 ∣ x^t)w(x^t + 1 ∣ x^t)=p(x^tt)q(x^t + 1 ∣ x^t)

• 提案分布q(x ∣ y)によりサンプルされたxをa(x ∣ y)の確率で採択することで、詳細つり合い条件を満たすようにサンプリングをしよう

このときの採択確率

$$ a(x^{t+1} \mid x^t) = \min \left[1, \frac{p(x^t)q(x^{t+1} \mid x^t)}{p(x^{t+1})q(x^t \mid x^{t+1})}\right] $$

うん、イメージがわかない

• 先人の手を借りましょう
• 可視化で理解するマルコフ連鎖モンテカルロ法

アイデア

• あるパラメータ以外を固定した条件付き分布からサンプリング
• 採択確率は常に1
• 提案分布と不変分布の関係を仮定できる

$$ q(x \mid x^t) = \begin{cases} p(x_i \mid x_{-i}^t) \; &x_{-i} = x_{-i}^t \\ 0 \; &otherwise \end{cases} $$

• これも可視化を見たほうが分かりやすい

バーンイン

• 最初の方は初期値依存があるので事後分布の形成から省くことが多い

MCMCでの収束

• Geweke の方法
• Geweke統計量Z値を計算し、MCMCのburn-in期間（前半10%）と後半50%を比較し、|Z|値が十分に小さければMCMCが収束したと考える

Stan とは

• Hamiltonian Monte Carlo 法(HMC)とそれを応用した No-U-Turn Sampler(NUTS)用いたサンプラー
  • 注: マルコフ連鎖を使っていないので厳密には MCMC ではない
• 要は「モデル式を書けば勝手に MC(Monte Carlo)してサンプリングしてくれるもの」
  • モデル式: 尤度関数と事前分布
  • C++で動くので速い
• 公式マニュアルが丁寧(PDF 注意)
  • 和訳もあるよ

おことわり

• ここでは「使ってみる程度」の説明しかできません。
• いろいろな書き方は以下の資料が詳しいです。
  • https://www.slideshare.net/simizu706/stan-62042940
  • https://www.slideshare.net/simizu706/stan-64926504

Stan のインストール

• OS ごとに違うので以下を見てください
• https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started-(Japanese)

rstan パッケージのインストール

• CRAN から

install.packages("rstan")

R で Stan を使う流れ

1. Stan コードを書く(sample.stan)
2. C++でコンパイル

stanmodel <- rstan::stan_model("sample.stan")

3. サンプリング

datlist <- list(...)
fit <- rstan::sampling(stanmodel, datlist, iter=10000, ...)

4. 結果を出力

とりあえず重回帰

• y = b₀ + b₁x₁ + b₂x₂ + e,  e ∼ N(μ, σ)
• サンプルデータ
  • b0=2, b1=5, b2=-3 のデータを 400 個生成

n <- 400
b0 <- 2
b1 <- 5
b2 <- -3
sigma <- 2

set.seed(71)
xx1 <- rnorm(n, 0, 1)
xx2 <- rnorm(n, 2, 2)
bx <- b0 + b1*xx1 + b2*xx2
yy <- rnorm(n, bx, sigma)

df <- data.frame(y = yy, x1 = xx1, x2 = xx2)

サンプルデータ

Stan の文法

• function{}
• data{}
• transformed data{}
• parameters{}
• transformed parameters{}
• model{}
• generated quantities{}

注: 順番はこの順でなければならない

重回帰を書いてみる

これを sample.stan として保存

data {
  int N ; // num of obs
  real x1[N] ; // independent variable 1
  real x2[N] ; // independent variable 2
  real y[N] ; // dependent variable
}

parameters {
  real b0 ;
  real b1 ;
  real b2 ;
  real<lower=0> sigma ;
}

model {
  for (i in 1:N) {
    y[i] ~ normal(b0 + b1 * x1[i] + b2 * x2[i], sigma) ;
  }
}

実行はRで

library(rstan)

# コンパイル
stanmodel <- stan_model("sample.stan")

datlist <- list(y = df$y,
                x1 = df$x1,
                x2 = df$y,
                N = nrow(df))

fit <- sampling(stanmodel, datlist)

収束を確認1

トレースプロットを見る

traceplot(fit)

収束を確認2

• Rhat という指標を見る
  • このあたりのパラメータの取り出し方は北條さん資料が詳しい
• これが 1.05(1.1とも)以下ならOK

summary(fit)$summary

             mean      se_mean         sd        2.5%         25%
b0       2.025081 0.0028269389 0.14803756    1.740184    1.923768
b1       4.834611 0.0016919143 0.10700606    4.628631    4.760272
b2      -2.933941 0.0009892484 0.05203129   -3.036034   -2.968873
sigma    2.046228 0.0011815352 0.07382550    1.906331    1.995531
lp__  -484.974192 0.0321805712 1.39814965 -488.418567 -485.714889
              50%         75%       97.5%    n_eff      Rhat
b0       2.025092    2.123802    2.321898 2742.275 0.9999702
b1       4.835004    4.906985    5.048263 4000.000 0.9999912
b2      -2.933319   -2.898268   -2.835307 2766.422 0.9996749
sigma    2.044362    2.096117    2.199582 3904.086 1.0002324
lp__  -484.645567 -483.942269 -483.207848 1887.643 1.0013678

Stan の書き方いろいろ

モデルを短く書きたい

for (i in 1:N) {
    y[i] ~ normal(b0 + b1 * x1[i] + b2 * x2[i], sigma) ;
  }

↓

data {
  int N ; // num of obs
  real x1[N] ; // independent variable 1
  real x2[N] ; // independent variable 2
  real y[N] ; // dependent variable
}

parameters {
  real b0 ;
  real b1 ;
  real b2 ;
  real<lower=0> sigma ;
}

transformed parameters {
  real mu[N] ;
  for (i in 1:N) {
    mu[N] = b0 + b1 * x1[i] + b2 * x2[i]
  }
}

model {
  for (i in 1:N) {
    y[i] ~ normal(mu[i], sigma) ;
  }
}

行列で書きたい

• データやパラメータを行列やベクトルで宣言する必要あり

data {
  int N ;
  int M ; \\ num of independent variables
  vector[N] y ; \\ dependent variable
  matrix[N, M] x ; \\ independent variable matrix
}

parameters {
  vector[M] b ; \\ parameter vector
  real<lower=0> sigma ;
}

model {
  y ~ normal(x*b, sigma) ;
}

R側も行列を作ります

datlist2 <- list(N = nrow(df),
                 y = df$y,
                 x = cbind(1, df$x1, df$x2),
                 M = ncol(x))

あれ、そういえば事前分布は？

• 何も指定しなければ無情報事前分布(一様分布)
• でも指定してあげたほうがいい(ベイズだしね！)
• model ブロックに以下を書く(さっきの行列の例)

y ~ normal(x*b, sigma) ;
beta ~ normal(0, 100) ; 
sigma ~ cauchy(0, 5) ;

　サンプリングの他の書き方

• y ~ normal(mu, sigma) は他の書き方もできる
• target += normal_lpdf(y | mu, sigma)
• ***lpdf は対数確率を足し上げる作業
  • 実はこっちがStan内部でやっていることで、 ~ のほうが簡略表記

Stan は楽でいいね！

• Stan を使えばモデルを書くだけでパラメータ推定ができる!
• 尤度さえわかればモデルが書ける

一方で

• モデルが正しく書けないとどうしようもない
  • モデル: 分布、行列etc…
• 発展的なモデル(階層モデル、状態空間モデル)を勉強するときはサボらずに尤度、行列表記くらいまではちゃんと理解できるよう努めましょう。

なにをやったか

• MCMC をざっくり理解
• Stan の書き方をざっくり理解

※正直なところ、30-40分では「ちゃんと」は理解できないと思ってます

初心者セッション 1-3 をもっと詳しく知りたい人は

R で楽しむベイズ統計

初心者セッション 1-3 をもっと詳しく知りたい人は

ベイズ統計モデリング: R,JAGS,Stan によるチュートリアル

初心者セッション 1-3 をもっと詳しく知りたい人は

Stan と R でベイズ統計モデリング

理論面を理解したい人

• 『データ解析のための統計モデリング入門』（通称: 緑本）
• 『計算統計 II－マルコフ連鎖モンテカルロ法とその周辺』（ガチなやつ）
• https://www.youtube.com/watch?v=-H28H1unn0M (MCMC全般)