サンプリングによる近似ベイズ推論その2（MCMC:メトロポリス法）

統計をベースとした機械学習では、観測データ $\mathcal{D}$ に基づき、確率分布を組み合わせてモデルを構築します。そして、モデルが観測データ $\mathcal{D}$ をうまく表現できるようにモデルのパラメータ $\Theta$ を推論します。こうして推定するパラメータ $\Theta$ の分布が事後確率分布 $p(\Theta | \mathcal{D})$ と呼ばれています。

事後確率分布 $p(\Theta | \mathcal{D})$ は確率モデル（パラメータの同時分布）に基づいて確率の基本性質に従うことで、導出できます（ベイズの公式）。しかし、複雑なモデルであったり、 $p(\mathcal{D} | \Theta), p(\Theta)$ が共役関係にない場合には、 $p(\Theta | \mathcal{D})$ を解析的に求めることができないのでした。

そこで、近似推論が必要になります。

learning-with-machine.hatenablog.com

前回は、サンプリングに基づく近似推論手法の中でも最もシンプルな「棄却サンプリング(Rejection Sampling)」を実装してみました。
今回は、MCMC法の一つである「メトロポリス法(Metropolis method)」を実装してみます。

【トップに戻る】

メトロポリス・ヘイスティングス法(Metropolis-Hastings method)

メトロポリス法やメトロポリス・ヘイスティングス法とは、MCMCと呼ばれる近似推論手法の一つです。 MCMC法というのは、高次元の空間で効率的にサンプリングを行う手段として知られており、マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo method)の略です。

具体的なアルゴリズムとしては、今回扱うメトロポリス法の他に、ハミルトニアンモンテカルロ法(Hamiltonian Monte Carlo method)やギブスサンプリング(Gibbs sampling)などが知られています。ギブスサンプリングは次回実装してみるつもりです。

メトロポリス・ヘイスティングス法のアルゴリズム

メトロポリス・ヘイスティングス法も棄却サンプリングと同様に、提案分布からサンプルを取得して、そのサンプルを目標分布に比例する分布に従って受容(accept)するか棄却(reject)するかを決めていくというアルゴリズムです。

まず、サンプルを得たい目標分布を $p(z)$ 、これに比例する分布 $\tilde{p}(z)$ を設定します。 $p(z)$ からサンプルを得たいけど、 $p(z)$ はどんな性質かわからないとします。次に、マルコフモデルにおける遷移確率 $p(z | z^{\prime})$ を設計するのですが、直接設計することが難しい場合には提案分布 $q(z | z^{\prime})$ を適当に設計します。（遷移確率を直接設計できる場合とはどのような場合なんでしょうか？解析的に定常分布に収束する分布が分かっているということ？）

そして、以下の手順でサンプルを取得します。

提案分布 $q(z_{*} | z^{t})$ から次のサンプル候補 $z_{*}$ をサンプル
比率r（下記）を計算
確率min(1,r)に従って $z _ *$ を受容し、状態を更新（ $z^{(t+1)} = z_{*}$ ）する。受容されない場合は状態を更新せず（ $z^{(t+1)}=z^{(t)}$ ）に1に戻る

${ \displaystyle \begin{eqnarray} r = \frac{ \tilde{p}(z_*)q(z^{(t)} | z_{*}) }{ \tilde{p}(z^{(t)}) q(z_*|z^{(t)}) } \end{eqnarray} }$

上記の比率rの計算で、提案分布 $q(z | z^{\prime})$ が対称であれば、提案分布に関する項は約分される。この場合を特に「メトロポリス法」と呼ぶそうです。

この比率rは、サンプル候補 $z_*$ の尤度 $\tilde{p}(z_*)$ が高いほど受容されやすくなるということを表しています。比率rに従って受容されるか否かが決まるので、確率が低い場合にも受容されることが有り得ます。そのため、尤度の山を行ったりきたりしながら登っていくようなイメージかなと思っています。

問題設定

前回の棄却サンプリングの例と同じで、ベルヌーイ分布のパラメータ $\theta$ の事後分布 $p(\theta | X)$ を推定する問題を考えます。具体的には、チョコボールの購入結果としてエンゼル有無の列をデータ $X$ とし、エンゼル出現確率 $\theta$ を推定するような問題です。

前節の手順に従って、まずは目標分布 $p(z)$ に比例する分布 $\tilde{p}(z)$ を設計します。目標分布は事後分布 $p(\theta | X)$ なので、ベイズの公式に従って次のように展開できます。

${ \displaystyle \begin{align} p(\theta | X) = \frac{p(X | \theta)p(\theta)}{p(X)} \propto p(X | \theta)p(\theta) \end{align} }$

ということで、目標分布 $p(z)$ に比例する分布 $\tilde{p}(z)$ は $\tilde{p}(z=\theta) = p(X | \theta)p(\theta)$ となります。また、確率の基本性質から、 $p(X, \theta) = p(X | \theta)p(\theta)$ ですので、モデルに登場する確率変数の同時分布が $\tilde{p}(z)$ になると考えてもOKです。

$p(X | \theta)$ はベルヌーイ分布ですので、 $p(\theta)$ はベータ分布を使いましょう。

${ \displaystyle \begin{align} \tilde{p}(\theta) = \mathrm{Bern}(X | \theta) \mathrm{Beta}(\theta | a,b) \end{align} }$

次は提案分布 $q(\theta)$ です。メトロポリス法の提案分布は、ガウス分布 $N(z | z^{\prime}, \sigma^2)$ がよく例として挙げられるようです。しかし今回の問題では、サンプルzの値域は0から1の範囲です（ベルヌーイ分布のパラメータは確率なので）。そこで、提案分布として一様分布を用いることにします。

一様分布は前の状態に依存しないので、 $q(z | z^{\prime}) = q(z^{\prime} | z)$ であり、対称です。なので、メトロポリス・ヘイスティングス法ではなく、メトロポリス法を使って近似推論します。

実装

以上の問題設定に従って、事後分布 $p(\theta | X)$ をメトロポリス法で近似推論してみます。コードの全体は、実装コード全体にnotebookを貼っていますので参照ください。

なお、真の確率 $\theta$ は0.1とします。観測データ $X$ はN=100個観測したとします。データ未観測でのパラメータ $\theta$ の分布（ベータ分布）は $\mathrm{Beta}(\theta | a=1, b=1)$ とします。

解析解

ここで扱う問題設定の場合、ベータ分布とベルヌーイ分布は共役関係にあるため、事後分布 $p(\theta | X)$ は厳密解を求めることができます。

${ \displaystyle \begin{align} p (\theta | X) &\propto p(X | \theta)p(\theta) \\ &= \left\{ \prod^{N}_{n=1} p(x_n | \theta) \right\} p(\theta) \\ &= \mathrm{Beta}(\theta | \hat{a}, \hat{b}) \end{align} }$

ここで、

${ \displaystyle \begin{eqnarray} \hat{a} = \sum^{N}_{n=1} x_n + a \\ \hat{b} = N - \sum^{N}_{n=1} x_n + b \end{eqnarray} }$

これを実装してみます。といっても、単純にベータ分布のパラメータを計算するだけですが。

# 事前分布のパラメータ
a = 1.0
b = 1.0
# 事後分布のパラメータ
a_hat = a + np.sum(X)
b_hat = b + N - np.sum(X)
# 事後分布の計算
analysis_dist = stats.beta(a=a_hat, b=b_hat)
ts = np.linspace(0, 1, 1000)
tt = analysis_dist.pdf(ts)

密度関数を確認してみます。設定値（黒線）をカバーするように事後分布が推定できていますね。

f:id:hippy-hikky:20191122102154p:plain — 事後分布の解析解．

近似推論

メトロポリス法により事後分布 $p(\theta | X)$ からサンプルを得ることを試してみます。応用や効率化など考えずにコード書いたので、その辺りは突っ込まないでいただけたらと思います。

## 提案分布
proposal = stats.uniform(loc=0, scale=1)
## 尤度
p_tilde_lik = stats.bernoulli
# 事前分布
a = 1
b = 1
p_tilde_pri = stats.beta(a=a, b=b)

def lik(X, z):
    lik = np.prod(p_tilde_lik.pmf(k=X, p=z))
    prior = p_tilde_pri.pdf(z)
    p_tilde = lik * prior
    return p_tilde

## メトロポリス法
n_sample = 1000
zs = []
zt = p_tilde_pri.rvs() # 初期状態は事前分布に従ってサンプル
cnt = 0
while len(zs) < n_sample:
    # 提案分布からサンプリング
    z = proposal.rvs()
    # 比率の計算
    r = lik(X, z) / lik(X, zt)
    # 受容判定
    u = stats.uniform.rvs(loc=0, scale=1)
    if u < r:
        zs.append(z)
        zt = z
    cnt+=1
    if cnt>n_sample*100:
        print('over iter')
        break

得られたサンプル（近似解）、先に求めた解析解、設定値を比較してみましょう。解析解を概ね近似していると言えそうですよね（？）

f:id:hippy-hikky:20191122142926p:plain — メトロポリス法による事後分布の推定結果．

次に、サンプル系列をプロットしてみます。下図左は全サンプルの系列。右はバーンインを200とした場合の系列（インデックスが200以降の系列）です。右図から、定常分布からサンプルできている（ランダムな系列）かなということがわかります。また左図から、極短い期間で定常分布に収束しているなということがわかります。

f:id:hippy-hikky:20191122143126p:plain — サンプル系列．左は全サンプルの系列．右はバーンインを200とした場合の系列．

今回の実験では、n_sample=1,000のサンプルを得るまでに10,699回ループが回っており、サンプル受容率は役10%でした。この辺りは、提案分布を工夫することで受容率は改善するだろうなと思います。（ベータ分布を使うなど）

比較

棄却サンプリングと推論にかかる時間やサンプル受容率を比較してみましょう。両者で問題設定やモデルは同じです。計算環境は僕のMacbookProで動かしています(。

CPU : 2.3 GHz Intel Core i5
メモリ : 16GB

比較結果は以下の通りです。今回の設定では受容率が10倍くらい違ってますね。

アルゴリズム	サンプル数	繰り返し回数	受容率	処理時間
棄却サンプリング	1,000	99,066	0.010	57.3sec
メトロポリス法	1,000	10,699	0.093	6.33sec

【トップに戻る】

実装コード全体

本記事で実装したコード全体は以下のnotebookの通りです。

【トップに戻る】

おわりに

ということで、サンプリングをベースにした実用的な近似推定アルゴリズムであるMCMC アルゴリズムの一つ、メトロポリス法を実装してみました。

棄却サンプリングと比較すると同じモデル、同じ提案分布を使っていても計算時間（≒サンプル受容率）が10倍くらい違ってきました。ただ、今回の実験では提案分布を一様分布にしているので、MCMC アルゴリズムの本当に美味しい使い方では無いんですよね。前の状態に対して独立になっていますので。
いずれ、もっと美味い例でアルゴリズム間の比較などをやっていこうと思います。

次回はMCMC アルゴリズムの一つであるギブスサンプリングを実装してみたいと思います。

【トップに戻る】