MCMC法の概要

なぜMCMCが必要か¶

積分計算の大変さ¶

ベイズ統計学では事後分布を推定してEAP推定量などの推定量を得るが、その際は積分の計算が出てくる（例えばEAPのための期待値計算での積分）

微分なら演算が閉じている（例えば $d x^2 / dx = 2x$ や $d \log x/dx = 1/x$ のように初等的な演算に還元される）が、積分は閉じていないので評価のコストが高い。

例えば、 $d$ 次元母数ベクトル $\boldsymbol{\theta}=\left(\theta_1, \theta_2, \cdots, \theta_d\right)$ があったとして、1番目の母数 $\theta_1$ のEAP推定量を求めるためには

\begin{aligned} \hat{\theta}_{1 \text { eap }} & =\int \theta_1 f\left(\theta_1 \mid \boldsymbol{x}\right) d \theta_1\\ & =\int \theta_1 \int \cdots \int \underbrace{ f\left(\theta_1, \theta_2, \cdots, \theta_d \mid \boldsymbol{x}\right) }_{ f(\boldsymbol{\theta} \mid \boldsymbol{x}) } d \theta_d d \theta_{d-1} \cdots d \theta_1 \\ \end{aligned}

(1)

のように、事後分布 $f(\boldsymbol{\theta} \mid \boldsymbol{x})$ の $d$ 次元の積分を評価しなければならない。この積分は一般的に、解析的に解くことができない。

事後分布から乱数を発生させて解く¶

そこで、事後分布 $f(\theta \mid x) \propto f(x \mid \theta)\, f(\theta)$ から母数の乱数 $\theta^{(1)}, \theta^{(2)}, \dots$ を発生させて、母数の実現値を取得することにする。母数の実現値たちが事後分布を十分に近似できていれば、その標本平均をEAPの推定値とすることができる。

ただし、このアプローチをとるにあたっては課題が2つある

どうやって乱数を発生させるか？
- 疑似乱数を発生させるプログラムでは、所与の分布のもとでの実現値 $f(x \mid \boldsymbol{\theta})$ を発生させるようにできており、事後分布 $f(\boldsymbol{\theta} \mid \boldsymbol{x})$ の生成のために作られていない
事後分布の正規化定数は評価できないことが多い
- 正規化定数 $\int_{\mathbb{R}^d} f(x)\,dx$ は $d$ 次元の積分計算を含み、簡単には評価できない。
- 正規化定数が不明だと事後確率密度そのものは計算できないが、その状態で カーネルだけで事後分布をどう計算するか？ という課題がある。

マルコフ連鎖¶

MCMC法では マルコフ連鎖（Markov chain） という確率過程を使う。マルコフ連鎖は

ある状態から、別の状態に推移する確率だけが決まっている
1期前だけが次の状態に影響する
連鎖はずっと続く（0%/100%がない）

という確率過程。

$t$ 期の確率変数を $Y_t$ とし、その確率過程を $\{Y_t\}^\infty_{t=0}$ と表すことにする。また確率変数の実現値を $y_t$ と表すことにする。

一般の確率過程においては、ある時点 $T$ の確率変数 $Y_T$ の確率密度関数は

f(y_T|y_0,y_1,\cdots,y_{T-1})

(2)

と、過去の $\{y_t\}^{T-1}_{t=0}$ に依存した条件付き確率密度となる。しかし、マルコフ連鎖においては1期前の確率変数の実現値にのみ依存する。

マルコフ連鎖（Markov chain）

ある時点 $T$ の確率変数 $Y_T$ の確率密度関数が、一期前の実現値 $y_{t-1}$ にのみ依存する、すなわち

f(y_t|y_0,y_1,\cdots,y_{t-1}) = f(y_t|y_{t-1})

(3)

を満たす確率過程を マルコフ連鎖（Markov chain） という。

マルコフ連鎖の重要な性質¶

MCMC法に関わる性質

既約性（irreducibility）
- どこから開始しても、有限回の連鎖で、確率密度が正であるすべての状態に必ずたどり着く（全空間を探索できる）
非周期性（aperiodicity）
- 一定に間隔で必ず訪れる状態がない（収束する）
正再帰性（positive recurrence）
- 平均的に、有限時間で戻って来る ⇔ 定常分布が存在する
- 遷移核が既約性と非周期性を満たすとき、 $t\to\infty$ で $y_t$ の定常分布（不変分布）に収束する。そしてその定常分布は連鎖の均衡分布である。（→この性質で事後分布を近似する）

遷移核¶

$\{Y_t\}^T_{t=0}$ の同時確率密度関数は、 $Y_0$ の周辺確率密度関数を $f_0$ とすると、

\begin{aligned} f(y_0, y_1, \ldots, y_T) & =f_0(y_0) f(y_1 \mid y_0) f(y_2 \mid y_0, y_1) \times \ldots \times f(y_T \mid y_0, \ldots, y_{T-1}) \\ & =f_0(y_0) f(y_1 \mid y_0) f(y_2 \mid y_1) \times \ldots \times f(y_T \mid y_{T-1}) \\ & =f_0(y_0) \prod_{t=1}^T f(y_t \mid y_{t-1}) \end{aligned}

(4)

となり、 $f_0$ と $f(y_t \mid y_{t-1})$ が与えられれば $f(y_0, y_1, \ldots, y_T)=f_0(y_0) \prod_{t=1}^T f(y_t \mid y_{t-1})$ で一意に決まる。

この $f(y_t \mid y_{t-1})$ のことを 遷移核（transition kernel, 推移核とも） という。

定常分布¶

マルコフ連鎖において、各状態に遷移していっても同じ分布のまま変わらなくなったときの分布を 定常分布という。

定常分布

ある特定の確率密度関数を $f$ で表し、状態 $\theta$ から状態 $\theta'$ へ移る確率（遷移核）を $f(\theta' \mid \theta)$ で表すとする。このとき

f(\theta^{\prime})=\int^{+\infty}_{-\infty} f(\theta) f(\theta^{\prime} \mid \theta) d \theta

(5)

が成り立つとき、確率密度関数 $f$ をもつ確率分布をマルコフ連鎖の 定常分布 (stationary distribution) あるいは 不変分布 (invariant distribution) という。

$t\to\infty$ において $f$ が定常分布に収束するとき、そのようなマルコフ連鎖は エルゴード性（ergodicity） を持つという。

エルゴード性（ergodicity）¶

既約性（irreducibility）
非周期性（aperiodicity）
正再帰性（positive recurrence）

を満たす連鎖においては、次の性質を持つ

定常分布 $\pi$ が一意に存在する
初期値によらず、 $\operatorname{Law}(X_t) \to \pi$ （ $\operatorname{Law}(X_t)$ は確率変数 $X_t$ の従う分布の意味）
連鎖の時間についての平均が分布の平均と等しくなる

エルゴード定理

連鎖の時間についての平均が分布の平均と等しくなる

\frac{1}{T} \sum_{t=1}^T f(X_t) \xrightarrow{\text { a.s. }} \mathbb{E}_\pi[f(X)]

(6)

→ MCMCでサンプリングした連鎖の実現値の標本期待値で事後期待値（EAP推定値）が推定できる

詳細つり合い条件¶

マルコフ連鎖が定常分布に収束するための十分条件に詳細釣り合い条件というものがある。

詳細釣り合い条件の式の両辺を $\theta'$ で積分すると

\begin{aligned} \int f(\theta \mid \theta') f(\theta') \, d\theta' &= \int f(\theta' \mid \theta) f(\theta) \, d\theta' \\ &= f(\theta) \int f(\theta' \mid \theta) \, d\theta' \\ % &= f(\theta) \int \frac{f(\theta',\theta)}{f(\theta)}\, d\theta'\\ &= f(\theta) \frac{1}{f(\theta)} \int f(\theta',\theta)\, d\theta'\\ &= f(\theta) \frac{1}{f(\theta)} f(\theta)\\ % &= f(\theta) \cdot 1 \\ &= f(\theta) \end{aligned}

(8)

となる。言葉で説明すると、 $\theta'$ から $\theta$ に遷移する確率密度の、あらゆる点に関する平均確率密度（左辺）は、 $\theta$ の確率密度 $f(\theta)$ そのものである（右辺） ということ。そしてこれが標本空間の任意の2点（ $\theta, \theta'$ ）で上記の条件が成り立つ、というのが詳細釣り合い条件。

MCMCのアルゴリズムは、この条件を満たすように遷移確率を作り、事後分布をサンプリングする。

MCMCのアルゴリズム¶

メトロポリス・ヘイスティング（MH）法¶

最も単純な方法でわかりやすいが、実用上はうまく収束しないことが多い方法。

基本アイデア

遷移核 $f( \mid)$ に正規分布など適当な分布 $q(\mid)$ を置く（ 提案分布 という）
提案分布が詳細釣り合い条件を満たさない場合、釣り合うようにズレたぶんだけ補正をしていく

アルゴリズムの概要：

遷移核にまず適当な確率分布を置く
- 正規分布がよく使われる
- $T(\theta'|\theta)=N(\theta, \sigma)$
その提案分布を補正する処置を施す
- $r N(\theta, \sigma)$
- 提案分布からの採用率を、尤度と事前分布から計算する

r=\frac{ N\left(\theta^{\prime} \mid \theta, \sigma\right) P(\theta \mid y) }{ N\left(\theta \mid \theta^{\prime}, \sigma\right) P\left(\theta^{\prime} \mid y\right) }

(9)

ハミルトニアン・モンテカルロ（HMC）法¶

ハミルトニアン方程式を使って推移確率行列を構成
勾配情報をつかい、運動エネルギーと位置を更新していく

NUTS¶

not u-turn
HMCの改良版
- 更新頻度とハイパーパラメータを自動調整
- Stan, NumPyro, PyMCなど多くのパッケージがNUTSを利用