統計モデリング概論 DSHC 2024

岩嵜 航 (Watal M. Iwasaki, PhD)
東北大学 生命科学研究科 進化ゲノミクス分野 特任助教
(Graduate School of Life Sciences, Tohoku University)
  1. 導入
  2. 直線回帰、確率分布、擬似乱数生成
  3. 尤度、最尤推定
  4. 一般化線形モデル(GLM)
  5. 個体差、一般化線形混合モデル(GLMM)
  6. ベイズの定理、事後分布、MCMC
  7. StanでGLM
  8. 階層ベイズモデル(HBM)
2024-08-28 東京海上 Data Science Hill Climb
https://heavywatal.github.io/slides/tokiomarine2024/

ちょっとずつ線形モデルを発展させていく

データ解析のための統計モデリング入門 久保拓弥 2012

久保先生の"緑本"こと
「データ解析のための統計モデリング入門」
をベースに回帰分析の概要を紹介。

線形モデル LM (単純な直線あてはめ)

    ↓ いろんな確率分布を扱いたい

一般化線形モデル GLM

    ↓ 個体差などの変量効果を扱いたい

一般化線形混合モデル GLMM

    ↓ もっと自由なモデリングを！

階層ベイズモデル HBM

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

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

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

plot of chunk freq-vs-bayesplot of chunk freq-vs-bayes

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

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

plot of chunk coin-frequentist

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

plot of chunk coin-bayesian

確率おさらい

同時分布/結合確率: AかつBの確率
$\Pr(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B}) = \Pr(\textcolor{#E69F00}{A} \cap \textcolor{#0072B2}{B}) = \Pr(\textcolor{#E69F00}{A}) \Pr(\textcolor{#0072B2}{B})$
周辺確率: BによらずAになる確率
$\Pr(\textcolor{#E69F00}{A}) = \sum_{\textcolor{#0072B2}{B}} \Pr(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B})$
条件付き確率: Bである条件の下でAになる確率。重要。
$\Pr(\textcolor{#E69F00}{A} \mid \textcolor{#0072B2}{B}) = \frac {\Pr(\textcolor{#E69F00}{A}, \textcolor{#0072B2}{B})} {\Pr(\textcolor{#0072B2}{B})}$

plot of chunk venn

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

B BreweryのビールがAwesomeな確率
$\Pr(\textcolor{#E69F00}{\text{Awesome}} \mid \textcolor{#0072B2}{\text{B Brewery}}) = \frac {\Pr(\textcolor{#E69F00}{\text{Awesome}},~\textcolor{#0072B2}{\text{B Brewery}})} {\Pr(\textcolor{#0072B2}{\text{B Brewery}})}$
かなり高い確率。良い醸造所。
AwesomeなビールがB Breweryのものである確率
$\Pr(\textcolor{#0072B2}{\text{B Brewery}} \mid \textcolor{#E69F00}{\text{Awesome}}) = \frac {\Pr(\textcolor{#E69F00}{\text{Awesome}},~\textcolor{#0072B2}{\text{B Brewery}})} {\Pr(\textcolor{#E69F00}{\text{Awesome}})}$
かなり低い確率。Awesomeなビールはほかにもたっくさんある。
plot of chunk venn

ベイズの定理

乗法定理
$\Pr(\textcolor{#E69F00}{A},~\textcolor{#0072B2}{B}) = \Pr(\textcolor{#E69F00}{A} \mid \textcolor{#0072B2}{B})~\Pr(\textcolor{#0072B2}{B}) = \Pr(\textcolor{#0072B2}{B} \mid \textcolor{#E69F00}{A})~\Pr(\textcolor{#E69F00}{A})$

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

事後確率 事前確率
\[ \Pr(\textcolor{#0072B2}{B} \mid \textcolor{#E69F00}{A}) = \frac {\Pr(\textcolor{#E69F00}{A} \mid \textcolor{#0072B2}{B})~\Pr(\textcolor{#0072B2}{B})} {\Pr(\textcolor{#E69F00}{A})} \]

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

事前確率: $\Pr(\textcolor{#0072B2}{B})$
飲む前、手元のビールがB Breweryのである確率。
↓ 🍻 飲んでみて更新
事後確率: $\Pr(\textcolor{#0072B2}{B} \mid \textcolor{#E69F00}{A})$
飲んでみてAwesomeだったビールが B Breweryのである確率。

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

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

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

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

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

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

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

事後分布 尤度 事前分布
\[ \Pr(M \mid D) = \frac {\Pr(D \mid M)~\Pr(M)} {\Pr(D)} \]
周辺尤度

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

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

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

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

plot of chunk posterior-beta  $\propto$
plot of chunk likelihood-binom  ⨉
plot of chunk prior-beta

コイントスを繰り返して、表が出る確率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)$:

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

plot of chunk posterior-beta  $\propto$
plot of chunk likelihood-binom  ⨉
plot of chunk prior-beta

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

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

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

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

plot of chunk posterior-beta  $\propto$
plot of chunk likelihood-binom  ⨉
plot of chunk prior-beta

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

\[\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$が裏の回数分だけ増加。

表が出る確率のベイズ推定: 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}\]

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

plot of chunk coin-bayesian

共役事前分布

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

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

共役事前分布を使うことが常に最善とは限らない。
無情報事前分布を使って計算機に頑張らせる風潮。

事後分布を用いた推定

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

plot of chunk integrate

ベイズ推定の中間まとめ

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

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

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

MCMC: Marcov Chain Monte Carlo

マルコフ連鎖
次の時点の挙動が現在の値だけで決定されるような確率過程。
$\ldots \to X_{t - 2} \to X_{t - 1} \to X_{t} \to X_{t + 1}$
$\Pr(X_{t+1} \mid X_{t}, X_{t-1}, X_{t-2}, \ldots) = \Pr(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 $\;\sim\;$ plot of chunk mcpdf

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

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

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

plot of chunk circle
  • 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法)

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

plot of chunk metropolis

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

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

plot of chunk metropolis-trajectory

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

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

plot of chunk propto-lik $\;\sim\;$ plot of chunk propto-lik $\;\propto\;$ 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 = 600 で灰色の部分を捨てる:

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サンプル増やす vs データ増やす

plot of chunk propto-lik $\;\sim\;$ plot of chunk propto-lik $\;\propto\;$ plot of chunk propto-lik
  • MCMCサンプルを増やす → 事後分布・尤度関数をより良く近似
  • データを増やす → 分布の裾野が狭まり、確信が強まる
plot of chunk coin-bayesian

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

Stan

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

Python Interface

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

参考文献

7. StanでGLM