統計モデリング実習 2022 TMDU

1. 導入、直線回帰
2. 確率分布、擬似乱数生成
3. 尤度、最尤推定
4. 一般化線形モデル (GLM)
5. 個体差、一般化線形混合モデル (GLMM)
6. ベイズの定理、事後分布、MCMC
7. StanでGLM
8. 階層ベイズモデル (HBM)

コイントス4回、たまたま表が1回だったら

コイントスの回数が増えていったら

最尤推定: 推定値が真の値に近づいていく

ベイズ推定: 確率分布がどんどん尖り、確信が強まる

確率おさらい

同時分布/結合確率: AかつBの確率

$\text{Prob}(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B}) = \text{Prob}(\textcolor{#E69F00}{A} \cap \textcolor{#0072B2}{B}) = \text{Prob}(\textcolor{#E69F00}{A}) \text{Prob}(\textcolor{#0072B2}{B})$

周辺確率: BによらずAになる確率

$\text{Prob}(\textcolor{#E69F00}{A}) = \sum_{\textcolor{#0072B2}{B}} \text{Prob}(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B})$

条件付き確率: Bである条件の下でAになる確率。重要。

$\text{Prob}(\textcolor{#E69F00}{A} \mid \textcolor{#0072B2}{B}) = \frac {\text{Prob}(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B})} {\text{Prob}(\textcolor{#0072B2}{B})}$

条件付き確率がよくわかる具体例

B BreweryのビールがAwesomeな確率

$\text{Prob}(\textcolor{#E69F00}{\text{Awesome}} \mid \textcolor{#0072B2}{\text{B Brewery}}) = \frac {\text{Prob}(\textcolor{#E69F00}{\text{Awesome}},~\textcolor{#0072B2}{\text{B Brewery}})} {\text{Prob}(\textcolor{#0072B2}{\text{B Brewery}})}$

かなり高い確率。良い醸造所。

AwesomeなビールがB Breweryのものである確率

$\text{Prob}(\textcolor{#0072B2}{\text{B Brewery}} \mid \textcolor{#E69F00}{\text{Awesome}}) = \frac {\text{Prob}(\textcolor{#E69F00}{\text{Awesome}},~\textcolor{#0072B2}{\text{B Brewery}})} {\text{Prob}(\textcolor{#E69F00}{\text{Awesome}})}$

かなり低い確率。Awesomeなビールはほかにもたっくさんある。

ベイズの定理

乗法定理

$\text{Prob}(\textcolor{#E69F00}{A},~\textcolor{#0072B2}{B}) = \text{Prob}(\textcolor{#E69F00}{A} \mid \textcolor{#0072B2}{B})~\text{Prob}(\textcolor{#0072B2}{B}) = \text{Prob}(\textcolor{#0072B2}{B} \mid \textcolor{#E69F00}{A})~\text{Prob}(\textcolor{#E69F00}{A})$

移項するだけでベイズの定理:

宴会場にビールが運ばれてきた。これはどこのブルワリーの？

事前確率: $\text{Prob}(\textcolor{#0072B2}{B})$

飲む前、手元のビールがB Breweryのである確率。

↓ 🍻 飲んでみて更新

事後確率: $\text{Prob}(\textcolor{#0072B2}{B} \mid \textcolor{#E69F00}{A})$

飲んでみてAwesomeだったビールが B Breweryのである確率。

ベイズの定理 in 感染症検査

• 有病率 $\text{Prob}(I)$ : 0.3% (この地域の感染者の割合; 事前確率)
• 感度 $\text{Prob}(P \mid I)$ : 90% (感染してる人に陽性判定が出る)
• 特異度 $\text{Prob}(\lnot P \mid \lnot I)$: 99% (感染してない人に陰性判定が出る)

さて、陽性適中率(検査陽性の人が実際に感染者である確率)は？

感染者を隔離するスクリーニング目的では使いものにならない性能。

🔰 同様に $\text{Prob}(\lnot I \mid \lnot P)$ 陰性的中率を計算してみよう
🔰 計算結果が検査性能だけでなく有病率にも依存することを確認しよう

ベイズの定理 in 統計モデリング

モデル$M$に対する確信度合いをデータ$D$に基づいて更新する。
モデル$M$を仮説$H$やパラメータ$\theta$に置き換えてもいい。

周辺尤度は「確率分布の積分は1」を満たすための正規化定数とみなせる。
比例関係だけ抜き出してこう書くことが多い:

表が出る確率のベイズ推定: 1. 事前分布

コイントスを繰り返して、表が出る確率pをベイズ推定したい。

事前分布にはベータ分布を採用(理由は後で分かる):

分布の形は $a,~b$ によって決まる。
ガンマ関数の部分は厳つく見えるけどただの正規化定数。
投げる前なのでとりあえず真っ平らを仮定 $\text{Beta}(p \mid a = 1, b = 1)$:

表が出る確率のベイズ推定: 2. 尤度関数

4回投げて表が1回だった、というデータで尤度を計算(二項分布):

これに事前分布を掛けて正規化したら事後分布になるはず。

 ⨉

表が出る確率のベイズ推定: 3. 事後分布

なんと、事後分布もベータ分布になる。

ベータ分布の形パラメータ$a$が表、$b$が裏の回数分だけ増加。

 $\propto$

 ⨉

表が出る確率のベイズ推定: 4. 逐次学習

さっきの事後分布を事前分布として、さらにデータを集める。

コイントス4回のうち表1回、に基づく事前分布: $\text{Beta}(p \mid 2,~4)$

さらに16回投げたら表が7回、の尤度: $\text{Binomial}(7 \mid 16,~p)$

事後分布はまた事前分布と同じ形になる:

データを加えるたびに更新していける:

共役事前分布

事後分布が事前分布と同じ形なので計算しやすい、という組み合わせ。

共役事前分布を使うことが常に最善とは限らない。
計算コストがかかっても無情報事前分布を使う風潮。

事後分布を用いた推定

ベイズ推定の中間まとめ

• 推定結果は事後分布 ∝ 尤度関数。
  • 広がり具合によって不確実性も表現できる。
  • 逐次学習で尖っていき、確信が強まる。

コイン投げモデルのベータ分布は美しい例。
→ 解析的(数学的)に解ける。

実践的なモデル・事後分布はもっと複雑。
→ コンピュータに頼って数値計算: MCMC

MCMC: Marcov Chain Monte Carlo

マルコフ連鎖

次の時点の挙動が現在の値だけで決定されるような確率過程。

$\ldots \to X_{t - 2} \to X_{t - 1} \to X_{t} \to X_{t + 1}$

$\text{Prob}(X_{t+1} \mid X_{t}, X_{t-1}, X_{t-2}, \ldots) = \text{Prob}(X_{t+1} \mid X_{t})$

e.g., すごろく

モンテカルロ法

乱数を用いた計算方法の総称。

モンテカルロ法は乱数を用いた計算方法

e.g., 半径1の円の面積

数学を知っていれば $\pi r ^ 2 \approx 3.14159$

面積4の正方形に400個の一様乱数を打ち込んだら318個が円に乗った:
$4 \times \frac {318} {400} = 3.18$

変な分布もモンテカルロ法で扱えそう

e.g., 確率密度分布に従って変数Xを集める(棄却サンプリング)。

$\;\sim\;$

でも、ハズレの値もけっこう引いてしまう。

次元の呪い: 高次元になるほど当たりにくくなる

(N次元球の体積 / N次元の立方体) はゼロに近づいていく。

• 2次元: $\frac {\pi r ^ 2} {(2r) ^ 2} = \frac \pi 4 \approx 0.79$
• 3次元: $\frac {\frac 4 3 \pi r ^ 3} {(2r) ^ 3} = \frac \pi 6 \approx 0.52$
• N次元: $\frac {\frac {\pi ^ {N/2}} {\Gamma (N/2 + 1)} r ^ N} {(2r) ^ N} = \frac {\pi ^ {N/2}} {2^N \Gamma (N/2 + 1)} \to 0$

パラメータが増えると計算量(≈乱数の無駄撃ち)が急増。

密度の高い「当たり」付近を効率よく探索したい。
「当たり」は「当たり」の近くにありがちだろう。
→ マルコフ連鎖が使えそう

Metropolis–Hastings法 (MH法)

0. パラメータ $\theta$ の初期値を選ぶ
1. $\theta$ をちょっと増減させて $\theta_\text{new}$ を作る
2. それぞれ尤度を計算し、比較。
   • $L(\theta_\text{new}) \ge L(\theta)$ なら $\theta_\text{new}$ を即採択
   • $L(\theta_\text{new}) < L(\theta)$ でも 確率 $r = \frac {L(\theta_\text{new})} {L(\theta)}$ で $\theta_\text{new}$ を採択
3. $\theta_\text{new}$ が採択されたら $\theta$ を更新。手順1に戻る。

採択されたパラメータ値の軌跡

尤度が高い方にただ向かうだけでなく、結構うろつく。
通ったパラメータ値を集めるといい感じの分布が得られる。

尤度に比例する事後分布からサンプルしたのと等価

全体にばら撒く棄却サンプリングよりも効率よく集められる。
が、パラメータ1つの1次元ではご利益はわかりにくい。

$\;\sim\;$ $\;\propto\;$

パラメータが複数ある場合は？

Gibbs Sampling

パラメータが複数の場合「ほかを固定してひとつ更新」を繰り返す。

e.g., 二次元正規分布。(-2, 2) からスタート。

何回やっても似たような結果になってほしい

乱数や初期値によって偶々、じゃないことを確認したい。

e.g., chains = 3 。ほぼ同じところをうろうろ:

収束(convergence)の判定については後ほど。

初期値の影響が消えるまでウォーミングアップ

定常分布の山に到達してからが本番。

e.g., iter_warmup = 200, iter_sampling = 600 で灰色の部分を捨てる:

どれくらい長く捨てるべきかは場合による。

適度に間引いて自己相関を軽減したい

直前の値と似すぎていたら独立サンプルとして扱えないので。

e.g., thin = 5 で5回に1回だけサンプルする:

間引かなくても大丈夫な場合も、間引いても解決しない場合もある。

収束判定

• 複数chainsで異なる初期値から実行し、軌跡を可視化(traceplot)
• Gelman-Rubin統計量 $\hat R < 1.05$
• Effective Sample Size (ESS) $N_\text{eff} > 100$ per chain

diagnose() みたいな機能が提供されていれば利用する。

実行時に警告してもらえることもある。

収束・自己相関が悪い場合にどう改善するか

• 小手先の対処
  • iteration (warmup + sampling) をもっと長く
  • thinを大きくして間引く
• ちょっと大掛かり
  • プログラムの書き方を改善する
  • モデルの構造を見直す
  • アルゴリズム・ソフトウェアを変える

MCMCの方法いろいろ

採択率を高め、早く収束するように改良されてきている。

• Metropolis–Hastings法
  • Gibbs Sampling
  • Hamiltonian Monte Carlo (HMC)
    • No-U-Turn Sampler (NUTS)

MCMCソフトウェア

• BUGS
  • クローズドソースで、ほぼWindows専用。
• JAGS
  • オープンソースで、さまざまなOS・言語から利用可能。
  • マニュアルや用例が不足。
• Stan 👈
  • オープンソースで、さまざまなOS・言語から利用可能。
  • 開発も利用も活発。マニュアルや用例も充実。
  • HMC/NUTSにより早く収束。
• PyMC3
• NumPyro
• TensorFlow Probability

ベイズ推定とMCMCまとめ

• 推定結果は事後分布 ∝ 尤度関数。
  • 広がり具合によって不確実性も表現できる。
  • 逐次学習で尖っていき、確信が強まる。
• コンピュータに頼った計算方法としてのモンテカルロ法。
  • 高次元・変な形の分布でも効率よくサンプルできるかが肝。
  • ソフトウェアの代表はStan → これから見ていく

Stan

• Stan言語でモデルを柔軟に記述できる。
• C++で書かれていて高速に動作。
• RやPythonなどから呼び出して使うのが便利。

R Interface

RStan

Rcppを介してStanHeadersを取り込んだパッケージ。

CmdStanR 👈 今後の主流

CmdStan を呼び出し、書き出されたCSVを読み取る。

CmdStanR: インストールがちょっと特殊

実行の前後にRを再起動してまっさらにすることを推奨。

1. C++言語を扱うためのツールを用意。
   • macOS: Command Line Tools (xcode-select --install)
   • Windows: RTools (Rのバージョンに合わせる)
2. CmdStanRパッケージをインストール。 (まだCRANに登録されていない):

   install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
   

3. CmdStan本体と可視化パッケージのインストール:

   library(cmdstanr)
   check_cmdstan_toolchain(fix = TRUE)
   install_cmdstan()
   install.packages("bayesplot")
   

https://mc-stan.org/cmdstanr/articles/cmdstanr.html

🔰 とりあえずStanを動かしてみよう

library(cmdstanr)
library(bayesplot)

This is cmdstanr version 0.5.3
- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
- CmdStan path: /Users/watal/.cmdstan/cmdstan-2.31.0
- CmdStan version: 2.31.0

This is bayesplot version 1.10.0

おおまかな流れ:

1. データ準備
2. Stan言語でモデルを書く
3. モデルをコンパイルして機械語に翻訳 → 実行ファイル
4. 実行ファイルにデータを渡してMCMCサンプリング
5. 結果を見る

説明変数なしのベイズ推定: データ準備

表が出る確率 $p=0.7$ のイカサマコインをN回投げたデータを作る。
この $p$ をStanで推定してみよう。

true_p = 0.7
N = 40L
coin_data = list(N = N, x = rbinom(N, 1, true_p))
print(coin_data)

$N
[1] 40

$x
 [1] 1 0 0 1 1 1 1 0 1 0 1 1 0 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 1 0 1 1 0 1 1 0 1 1 1 1

Rならlist型、Pythonならdict型にまとめてStanに渡す。

説明変数なしのベイズ推定: Stan言語でモデル定義

別ファイルに書いておく。 e.g., coin.stan:

data {
  int<lower=0> N;
  array[N] int x;
}
parameters {
  real<lower=0,upper=1> p;
}
model {
  x ~ binomial(1, p);
}

• いくつかのブロックに分けて記述する:
  R/Pythonから受け取る data, 推定する parameter, 本体の model.
• 変数には型や制約を設定できる
• 関数もたくさん用意されている

Stan言語の7種のブロック

順番厳守。よく使うのは太字のやつ。

1. functions {...}
2. data {...}
3. transformed data {...}
4. parameters {...}
5. transformed parameters {...}
6. model {...}
7. generated quantities {...}

https://mc-stan.org/docs/reference-manual/overview-of-stans-program-blocks.html

説明変数なしのベイズ推定: MCMCサンプル

予め実行速度の速い機械語に翻訳(コンパイル):

model = cmdstanr::cmdstan_model("stan/coin.stan")

モデルとデータを使ってMCMCサンプリング:

fit = model$sample(coin_data, seed = 24601L)

いろいろオプションはあるけど、ここではデフォルトに任せる:
chains, inits, iter_warmup, iter_samples, thin, …

問題があったら警告してくれるのでちゃんと読む。

説明変数なしのベイズ推定: 結果を眺める

parameters ブロックに書いた変数の情報が出てくる。
乱数を使った計算なので(乱数シードを固定しない限り)毎回変わる。

print(fit)

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
     lp__ -25.62 -25.35 0.70 0.30 -27.00 -25.13 1.00     1949     2646
     p      0.72   0.72 0.07 0.07   0.60   0.82 1.00     1586     2132

真の値に近い $p \approx 0.7$ が得られた (0.6 から 0.82 である確率が90%)。
$\hat R$ もほぼ1で $N_\text{eff}$ も大きいのでよさそう。

lp__ はlog posterior(対数事後確率)。後述。

念のため trace plot も確認しておこう→

説明変数なしのベイズ推定: trace plot 確認

どのchainも似た範囲を動いていて、しっかり毛虫っぽい:

draws = fit$draws()
params = names(model$variables()$parameters)
bayesplot::mcmc_trace(draws, pars = params)

説明変数なしのベイズ推定: 自己相関の確認

2–3ステップくらいで自己相関がほぼ消えるので問題なし:

bayesplot::mcmc_acf_bar(draws, pars = params)

説明変数なしのベイズ推定: 推定結果確認

サンプルサイズNが小さいせいか裾野の広い推定結果。
真の$p$の値も含まれている:

bayesplot::mcmc_hist(draws, bins = 20, pars = params)

lp__: log posterior とは?

model ブロックに次のように書いてあると:

model {
  mu ~ normal(0.0, 10.0);  // prior
  x ~ normal(mu, 1.0);     // likelihood
}

内部的には次のような処理が行われている:

target += normal_lpdf(theta | 0.0, 10.0)  // prior
target += normal_lpdf(x | theta, 1.0);    // likelihood

つまり、事前確率と尤度の対数の和を取っている。
ベイズの定理により、事後確率はこれに比例する。
lp__ はこの target 変数を記録しておいたようなもの。

ベイズ推定とMCMC

• 推定結果は事後分布 ∝ 尤度関数。
  • 広がり具合によって不確実性も表現できる。
  • 逐次学習で尖っていき、確信が強まる。
• コンピュータに頼った計算方法としてのモンテカルロ法。
  • 高次元・変な形の分布でも効率よくサンプルできるかが肝。
  • ソフトウェアの代表はStan。
  • Stan言語でモデルを書き、RやPythonから動かす。

参考文献

• データ解析のための統計モデリング入門 久保拓弥 2012
• StanとRでベイズ統計モデリング 松浦健太郎 2016
• RとStanではじめる ベイズ統計モデリングによるデータ分析入門 馬場真哉 2019
• データ分析のための数理モデル入門 江崎貴裕 2020
• 分析者のためのデータ解釈学入門 江崎貴裕 2020
• 統計学を哲学する 大塚淳 2020