統計モデリング概論 DSHC 2022

岩嵜航 (Watal M. Iwasaki, PhD)

東北大学生命科学研究科進化ゲノミクス分野特任助教
(Graduate School of Life Sciences, Tohoku University)

2022-08-24 東京海上 Data Science Hill Climb https://heavywatal.github.io/slides/tokiomarine2022/

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

最尤推定: 推定結果は最も尤もらしい1点。; データが少ないとき過剰適合気味。; 表が出る確率 p = 0.25 のコインだろう。
(信じ難いけどデータはこう言っている)

ベイズ推定: 推定結果は確率分布そのもの。; データが少ないなりの不確実さも表現。; p = 0.25 らへんである確率は高いが、
p = 0.6 とかである可能性もまあある。

plot of chunk freq-vs-bayes

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

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

plot of chunk coin-frequentist

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

plot of chunk coin-bayesian

確率おさらい

同時分布/結合確率: 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})}$

plot of chunk venn

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

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% (感染してない人に陰性判定が出る)

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

\[\begin{split} \text{Prob}(I \mid P) &= \frac {\text{Prob}(P \mid I)~\text{Prob}(I)} {\text{Prob}(P)} \\ &= \frac {\text{Prob}(P \mid I)~\text{Prob}(I)} {\text{Prob}(P \mid I)~\text{Prob}(I) + \text{Prob}(P \mid \lnot I)~\text{Prob}(\lnot I)} \\ &= \frac {0.9 \times 0.003} {0.9 \times 0.003 + 0.01 \times 0.997} \approx 0.21 \end{split}\]

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

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

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

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

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

\[\begin{align} \text{Prob}(M \mid D) &\propto \text{Prob}(D \mid M)~\text{Prob}(M) \tag{Model}\\ \text{Prob}(H \mid D) &\propto \text{Prob}(D \mid H)~\text{Prob}(H) \tag{Hypothesis} \\ \text{Prob}(\theta \mid D) &\propto \text{Prob}(D \mid \theta)~\text{Prob}(\theta) \tag{Parameter} \end{align}\]

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

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

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

\[\begin{split} \text{Beta}(p \mid a, b) = \frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} p^{a-1} (1 - p)^{b-1} \end{split}\]

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

plot of chunk prior-beta

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

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

\[\begin{split} \text{Binom}(1 \mid 4,~p) = \binom {1} {4} p^{1} (1 - p)^{3} \end{split}\]

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

plot of chunk likelihood-binom

⨉

plot of chunk prior-beta

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

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

\[\begin{split} \text{Posterior} &\propto \text{Binom}(1 \mid 4,~p) \times \text{Beta}(p \mid 1, 1)\\ &= \binom {1} {4} p^{1} (1 - p)^{3} \times \frac{\Gamma(1 + 1)}{\Gamma(1) \Gamma(1)} p^{1-1} (1 - p)^{1-1} \\ &= C p^{2-1} (1 - p)^{4-1} \\ &= \text{Beta}(p \mid 2, 4) \end{split}\]

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

plot of chunk posterior-beta

$\propto$

plot of chunk likelihood-binom

⨉

plot of chunk prior-beta

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

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

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

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

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

\[\begin{split} \text{Beta}(p \mid 9, 13) \propto \text{Binomial}(7 \mid 16,~p) \times \text{Beta}(p \mid 2, 4) \end{split}\]

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

共役事前分布

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

尤度関数	共役事前分布
二項分布	ベータ分布
ポアソン分布	ガンマ分布
正規分布	ガンマ分布
正規分布 (分散既知)	正規分布

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

事後分布を用いた推定

区間推定: 幅のある推定値を提示; e.g., 95%ベイズ確信区間:
等裾事後確信区間 (Equal-Tailed Interval)
最高密度区間 (Highest Density Interval)
点推定: 値を1点だけ提示; e.g.,
事後確率最大値 (Maximum A Posteriori)
事後期待値 (Expected A Posteriori)
事後中央値 (Posterior Median)

plot of chunk integrate

ベイズ推定の中間まとめ

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

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

実践的なモデル・事後分布はもっと複雑。
→ コンピュータに頼って数値計算: 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$

plot of chunk circle

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

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

plot of chunk mcpdf

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

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

(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法)

パラメータ $\theta$ の初期値を選ぶ
$\theta$ をちょっと増減させて $\theta_\text{new}$ を作る
それぞれ尤度を計算し、比較。
- $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}$ を採択
$\theta_\text{new}$ が採択されたら $\theta$ を更新。手順1に戻る。

plot of chunk metropolis

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

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

plot of chunk metropolis-trajectory

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

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

plot of chunk propto-lik

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

Gibbs Sampling

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

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

plot of chunk gibbs

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

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

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

plot of chunk chains

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

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

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

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

plot of chunk warmup

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

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

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

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

plot of chunk thin

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

収束判定

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

plot of chunk convergence

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などから呼び出して使うのが便利。

Python Interface

PyStan: WindowsやJupyterで使うには難あり。
CmdStanPy 👈 今後の主流: CmdStan を呼び出し、書き出されたCSVを読み取る。

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

おおまかな流れ:

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

🔰 4-bayesian.ipynb を開き、スライド説明に沿って実行しよう。

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

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

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

$N
[1] 40

$x
 [1] 1 1 1 0 0 1 1 1 0 1 1 1 0 0 1 1 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 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種のブロック

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

functions {...}
data {...}
transformed data {...}
parameters {...}
transformed parameters {...}
model {...}
generated quantities {...}

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

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

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

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

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

fit = model$sample(mydata)

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

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

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

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

print(fit)

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
     lp__ -23.59 -23.30 0.77 0.34 -25.17 -23.05 1.00     1760     2102
     p      0.76   0.77 0.07 0.07   0.65   0.86 1.00     1462     1749

真の値に近い $p \approx 0.7$ が得られた (0.65 から 0.86 である確率が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)

plot of chunk stan-binom-traceplot

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

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

bayesplot::mcmc_acf_bar(draws, pars = params)

plot of chunk stan-binom-ac

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

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

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

Warning: The following arguments were unrecognized and ignored: bins

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot of chunk stan-binom-hist

`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

5. StanでGLM

統計モデリング概論 DSHC 2022

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

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

確率おさらい

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

ベイズの定理

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

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

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

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

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

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

共役事前分布

事後分布を用いた推定

ベイズ推定の中間まとめ

MCMC: Marcov Chain Monte Carlo

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

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

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

Metropolis–Hastings法 (MH法)

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

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

Gibbs Sampling

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

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

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

収束判定

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

MCMCの方法いろいろ

MCMCソフトウェア

ベイズ推定とMCMCまとめ

Stan

Python Interface

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

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

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

Stan言語の7種のブロック

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

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

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

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

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

lp__: log posterior とは?

ベイズ推定とMCMC

参考文献

`lp__`: log posterior とは?