「Pythonではじめるベイズ機械学習入門」を読む　～第1章ベイジアンモデリング② MCMC～

はじめに

久保拓弥著「データ解析のための統計モデリング入門 ―一般化線形モデル・階層ベイズモデル・MCMC―」を読んだときに，ベイズモデルや階層ベイズモデルに触れた。本書ではコードがRで書かれていたので，ChatGPTに相談しながらCmdStanPyで実装した。
Pythonによるベイズモデルの実装をきちんと学ぼうと思い，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読むことにした。

本記事は，第1章「ベイジアンモデリング」の第1.4節のうちMCMCに関する読書メモである。

本書の紹介ページ

www.kspub.co.jp

第1章ベイジアン モデリングとは

本書は，確率的プログラミング言語(Probabilistic Programming Language; PPL)を用いたベイジアンモデリングや機械学習のデータ解析手法を構築するための実践書である。

本章では，

ベイジアンモデリングの基礎：同時確率，周辺確率と条件付き確率など
代表的な確率分布
近似推論手法：MCMC，変分推論法 など

が説明されている。

1.4 近似推論法

ベイズ推定で難しいのは，事後分布を求めることである。

たとえば回帰モデル $y = wx + \varepsilon$ の場合，興味があるのはデータ $x$ が得られた時のパラメータ $w$ の分布 $\pi(w | x)$ である。パラメータ $w$ が求められたら，新しい入力 $\hat{x}$ が得られた時の新しい出力 $\hat{y}$ が得られるためである。

事後分布は，ベイズの定理より

$\begin{align} \pi(w|x) = \frac{f(x|w) \pi(w) }{ p(x) } = \frac{f(x|w) \pi(w) }{ \int f(x|w) \pi(w) dw } \\ \end{align}$

となる。尤度 $f(x|w)$ は，回帰式を設計することで求められる。事前分布 $\pi(w)$ は，事前知識によって設計する。
ただし，データの発生機構 $p(x)$ はよくわからないし， $p(x) = \int f(x|w) \pi(w) dw$ には積分が発生する。そのため，事後分布を求めることは困難である。

本節では近似推論の方法として，

マルコフ連鎖モンテカルロ法 (Markov Chain Monte Carlo; MCMC)
変分推論法

を説明している。なお本書では，「ベイジアンモデリングを使ううえで推論手法の理論的な理解は必須というわけではない」という立場から，近似推論の手法はは要点を押さえた説明となっている。

マルコフ連鎖モンテカルロ法(MCMC)

事後分布を近似する方法の1つの手段は，事後分布の式を明示的に求める代わりに，事後分布に従う乱数(サンプル)をたくさん作成することである。
このようなサンプルが得られれば，分布の性質(平均や分散などの統計量)を計算することができる。

MCMCでは，初期値にかかわらず定常分布に収束するというエルゴード性(ergodicity)を満たす(ただし，収束にどの程度の時間がかかるかは別の問題である)。

代表的なMCMCのアルゴリズムとして，

メトロポリス・ヘイスティング法(Metropolis-Hasting algorithm, MH法)
ハミルトニアン・モンテカルロ法(Hamiltonian Monte Carlo method, HMC法)

などが挙げられる。

メトロポリス・ヘイスティング法の実験

本書のサポートページ・ソースコード(GitHub) のコードを参考に，メトロポリス・ヘイスティング法(MH法)を試した。

問題設定は以下の通りである。

データ

$\mathcal{N}(\mu, 1.0^2)$ にしたがう100個の正規乱数 $X = \{\ x_1, \cdots, x_{100} \}$ を発生させる。
なおサンプルコードでは，正解として $\mu = 3$ としている。

統計モデル

尤度は分散を既知として $f(x_n|\mu) \sim \mathcal{N}(x_n | \mu, 1.0^2)$ とする。
事前分布は， $\mu$ に関する事前知識がないので，分散が大きい正規分布として $\pi(\mu) \sim \mathcal{N}(0, 10^2)$ とする。

求めるもの

上記のデータから，パラメータ $\mu$ をベイズ推定する。事後分布 $\pi(\mu | x)$ は，MH法によって求める。

データのヒストグラムは下図のようになる。

なお上記のデータでは，標本平均は $\bar{X} = 3.06$ ，標本標準偏差は $s = 0.89$ である。

MH法の実装のポイントは以下の通りである。

提案分布

$q(\cdot | z^{(t)}) = \mathcal{N}(z^{(t)}, \sigma I)$ として， $\sigma = 0.3$ としている。

またサンプルコードでは，提案された点 $z_{*}$ は， $z_{*} = z^{(t)} + \sigma \cdot \varepsilon, \quad \varepsilon \sim \mathcal{N}(0,1)$ として求めている。
すなわち提案された点は，今あるサンプルからステップ幅 $\sigma$ だけ離れた点であると言える。分布を形成するサンプルを取得するためには，このステップ幅の選び方が重要になりそうである。

正規化されていない事後分布

正規化されていない事後分布は，尤度と事前分布の積としている。
また，正規化されていない事後分布は，そのまま計算すると小さい値になってしまうので， $p(z) = \exp(\log p(z))$ であることを利用して，アンダーフローを避ける。

MCMCのサンプルサイズ・初期値

サンプルサイズは6000とする。また初期値は0とする。

初期値と事後分布サンプル

MCMCの各ステップで得られたサンプルを結んだ折れ線プロットであるトレースプロットは下図のようになる。

真の平均は3である。初期値は0なので，最初は真の平均からずれた値になっているが，だんだんと真の値に近づいていることが分かる。

なお，事後分布のサンプル $x_{posterior}$ に，真の確率密度関数 $\mathcal{N}(3, 1.0^2)$ を当てはめて描画した結果が下図のようになる(なお，可視化用に，当てはめ結果には乱数を加えている)

真値 $\mu=3$ 付近で，事後分布 $\pi(\mu | x)$ のサンプリングができていることが確認できた。

次に，試みに初期値を10にして計算した結果が以下の通りである。

この例も同様に，だんだんと真の値に近づいてきていることが分かる。
MCMC法の特徴として，初期値に関係なく定常分布に近づくというものがあるが，実験的にもこのことは言えそうである。
ただし，初期値に近い部分は定常分布からは離れる可能性があるので，最初のバーンイン期間のサンプルは捨てる必要があることがわかる。

ステップサイズを変えた場合

先ほどの例では，ステップサイズを $\sigma = 0.3$ としたが，ステップサイズを $\sigma = 0.01$ および $\sigma = 10.0$ としてMH法を実行した結果が下図となる。

ステップサイズが小さすぎる，または大きすぎるため，定常分布に収束していないように見える。

原因としては，サンプルの系列の自己相関が大きくなる，すなわち各ステップで得られたサンプルが以前のサンプルの影響を強く受けており，定常分布からの独立なサンプルとみなせなくなっているためである。

まとめと感想

MCMCの特徴を把握できた

MH法やHMC法は，パッケージがあればすぐに実装できると思うが，原理を知っておくことは重要だと考える。
最近ではデータサイエンスについて学ぶハードルがかなり低くなっているが，分析者として差が出るのはトラブルシューティング，すなわち，うまくいかなかったときの対策を練られることだと考える。
対策を練るためには，うまくいかなかった理由が分かる必要があり，そのためにはアルゴリズムの特徴や原理を理解しておく必要がある。

本節では，NumPyなどの基本的なパッケージを用いてMH法の説明をしていた。特にサンプリングがうまくいかないパターンやその原因，確認方法について説明がされていた。
実務上はMCMCのパッケージを用いることになるが，うまくいかなかったときは原理的なところに立ち返ることができるよう，基本的なことは身に着けていきたい。