サンプリングによる近似ベイズ推論その1（モンテカルロ積分, 棄却サンプリング）

【概要】

ベイズ推論について実装して理解するシリーズの第２弾
今回からサンプリングに基づいた近似ベイズ推論を扱います
ここでは、モンテカルロ積分と棄却サンプリングを実装してみます

【目次】

はじめに
ベイズ推論の概要
サンプリングに基づいた近似ベイズ推論
モンテカルロ積分
- 問題設定
- 実装
  - 解析解
  - 近似解
棄却サンプリング
実装コード全体
終わりに
参考文献

はじめに

ベイズ推論についての書籍をいくつか読んでいて、なんとなく理解はできても具体的なイメージってつきにくくないですか？
ということで、実装して理解を深めていきたいと思います。

複雑なモデルや共役関係にない確率モデルでは、確率推論を解析的に実施するのが困難なので近似計算を利用する必要があります。ベイズ推論でよく用いられる近似推論手法には大きく分けて「MCMC(Markov Chain Monte Carlo)法」と呼ばれるアルゴリズム群と「変分推論」という手法があります。

今回から数回に渡って、近似推論手法の中でもMCMC法として知られているサンプリングに基づいたアルゴリズム群を実装して理解していこうと思います。

本記事は、「ベイズ深層学習」の「4章近似ベイズ推論」を参考にしています。

ベイズ深層学習 (機械学習プロフェッショナルシリーズ)

作者: 須山敦志
出版社/メーカー: 講談社
発売日: 2019/08/08
メディア: 単行本
この商品を含むブログを見る

間違いなどあったら指摘していただけると助かります。

【トップに戻る】

ベイズ推論の概要

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

パラメータの推論としては最尤推定がよく用いられます。最尤推定では、未知のパラメータは固定の値として推定していきます。一方ベイズ推論では、この未知のパラメータにも確率分布を仮定し、モデル全体を確率分布の組み合わせとします。推定対象となるパラメータは全て確率分布を持ちますので、ベイズ推論による推論結果は、観測データ $\mathcal{D}$ に基づいたパラメータ $\Theta$ の事後確率分布 $p(\Theta | \mathcal{D})$ が得られます。

事後確率分布 $p(\Theta | \mathcal{D})$ が得られたなら、未知の入力データ $x_{*}$ に対しての予測値は以下のように予測分布 $p(y_{*} | x_{*}, \mathcal{D})$ として得ることができます。

${ \displaystyle p(y_{*} | x_{*}, \mathcal{D}) = \int p(y_{*} | x_{*}, \Theta) p(\Theta | \mathcal{D}) d\Theta }$

詳しくは、ベイズ深層学習やPRML、また、当ブログの以前の記事「ベイズ線形回帰（解析解）の実装 - 機械と学習する」を参照してください。

事後確率分布 $p(\Theta | \mathcal{D})$ は、確率の性質にしたがって以下のように展開できます*1。

${ \displaystyle p(\Theta | \mathcal{D}) = \frac{p(\mathcal{D} | \Theta)p(\Theta)}{p(\mathcal{D})} }$

【トップに戻る】

サンプリングに基づいた近似ベイズ推論

事後確率分布 $p(\Theta | \mathcal{D})$ は上記の式に従うことで計算できます。

しかし、モデルのパラメータ $\Theta$ が多数あるような複雑なモデルであったり、 $p(\mathcal{D} | \Theta), p(\Theta)$ が共役関係という特別な関係にない場合には、 $p(\Theta | \mathcal{D})$ を解析的に求めることができない（または困難）です。そこで、近似推論が必要になります。

近似推論には、二つの方法がよく使われます。

サンプリングに基づく手法：求めたい分布 $p(\Theta | \mathcal{D})$ から何らかの方法でサンプルを多数取得し、そのサンプルの統計的性質から求めたい分布 $p(\Theta | \mathcal{D})$ の特性を知る
変分推論：数値最適化に基づき、ある簡単な分布 $q(\Theta)$ を求めたい分布 $p(\Theta | \mathcal{D})$ に近づけていく

まずは、今回から数回にわたって、MCMCなどサンプリングに基づく手法を実装して理解を進めていきます。

今回は二つのアルゴリズムを実装します。
最初に、確率推論の世界では予測分布や周辺尤度などを求めるために積分が非常に重要になります。その積分を数値的に近似計算するための手法として「モンテカルロ法」という単純なアルゴリズムを実装してみます。
次に、モンテカルロ法では確率分布からサンプルを得られることが前提ですが、複雑なモデルの事後確率分布 $p(\Theta | \mathcal{D})$ は解析的に求められなかったのでした。そこで、そのような確率分布からサンプルを得るための最も単純な方法である「棄却サンプリング」というアルゴリズムを実装してみます。

【トップに戻る】

モンテカルロ 積分

分布 $p(\mathbf{z})$ に関して、ある関数 $f(\mathbf{z})$ の期待値 $\int f(\mathbf{z})p(\mathbf{z})d\mathbf{z}$ を計算することを考えます。このような積分は、パラメータの期待値を単純に知りたい場合やベイズの公式の分母（周辺尤度）を求める場合など頻出する計算です。

ここで、分布 $p(\mathbf{z})$ からT個のサンプルを得られたとすると、以下のように積分は総和の形で近似することができます。

${ \displaystyle \int f(\mathbf{z})p(\mathbf{z})d\mathbf{z} \approx \frac{1}{T}\sum^{T}_{t=1} f(\mathbf{z}^{(t)}) }$

ここで、 $\mathbf{z}^{(t)}$ は分布 $p(\mathbf{z})$ からのt番目のサンプルを表します。確率の高い領域はそれだけ多くのサンプル $\mathbf{z}^{(t)}$ が得られるので、サンプルの量（密度）が確率の大きさを表現しています。期待値の式を見ると、確率で重みつけられた平均のように見えます。そのため、サンプルの量が重みとなり、 $f(\mathbf{z}^{(t)})$ の平均が期待値を近似するというイメージです。

問題設定

ここでは、モンテカルロ積分を使った具体例として、周辺尤度（モデルエビデンス）を近似計算してみます。

N個のデータ $X=\{x_1, x_2, \cdots, x_N\}$ が観測されており、パラメータ $\theta$ を持つモデル $p(\mathbf{X}, \theta)$ を考えます。周辺尤度は以下のように近似できます。

${ \displaystyle p(\mathbf{X}) = \int p(\theta) p(\mathbf{X} | \theta) d\theta \approx \frac{1}{T} \sum^{T}_{t=1} \prod^{N}_{n=1} p(\mathbf{X} | \theta^{(t)}) }$

ここでは、データ $x_n$ が従う分布としてベルヌーイ分布を設定し、パラメータ $\theta$ が従う分布としてベータ分布を設定することにします。

${ \displaystyle \begin{eqnarray} p(x_n | \theta) = \mathrm{Bern}(x_n | \theta) \\ p(\theta) = \mathrm{Beta}(\theta | a, b) \end{eqnarray} }$

ベルヌーイ分布とベータ分布は共役な関係なので、厳密な解を求めることも可能です。

$\begin{eqnarray} p(X) = \int \mathrm{Beta}(\theta | a, b) \prod^{N}_{n=1} \mathrm{Bern}(x_n | \theta) d\theta \\ \ \ = C_{B} \int^{1}_{0} \theta^{m+a-1} (1-\theta)^{N-m+b-1} d\theta \end{eqnarray}$

ここで、 $C_B$ はベータ分布の正規化定数です。また、 $m=\sum x_i$ を表しています。上記の積分はベータ関数の積分公式を利用して次のように展開できます。

${ \displaystyle \begin{align} p(X) &= C_{B} \frac{m^{\prime}!n^{\prime}!}{(m^{\prime}+n^{\prime}+1)!} \\ &= \frac{\Gamma(a+b)\Gamma(m+a)\Gamma(N-m+b)}{\Gamma(a)\Gamma(b)\Gamma(a+b+N)} \end{align} }$

途中で出てくる $m^{\prime}, n^{\prime}$ は以下の通りです。

$m^{\prime} = m+a-1, \\ n^{\prime} = N-m+b-1$

実装

では上記の問題設定にしたがって、モンテカルロ積分を実装してみます。コードの全体は、実装コード全体にnotebookを貼っていますので参照ください。

ベルヌーイ分布のパラメータ（確率, $\theta$ ）は0.1とし( $\mathrm{Bern}(x | \theta=0.1)$ )、20回試行したとします( $X=\{x_1, \cdots, x_{20}\}$ )。パラメータ $\theta$ の分布であるベータ分布のパラメータは(a=1, b=1)とします（一様分布, $\mathrm{Beta}(\theta | a=1, b=1)$ ）。

解析解

解析解は先に導出した式に従えば得られます。しかし、ガンマ関数（階乗）の積がたくさん出てくるので対数をとった対数周辺尤度を計算しています。

from scipy.special import gammaln

m = sum(X)
ln_px_c = gammaln(a+b)+gammaln(m+a)+gammaln(N-m+b)-gammaln(a)-gammaln(b)-gammaln(a+b+N)
px_c = np.exp(ln_px_c)

近似解

モンテカルロ積分は、分布 $p(\mathbf{z})$ から得たT個のサンプル毎に $f(z^{(t)})$ を計算して平均をすれば良いだけでした。

def approximate_marginal_likelihood(pz, sample_size=100000):
    zs = pz.rvs(size=sample_size)
    px = np.sum(list(map(lambda z:np.prod(stats.bernoulli.pmf(k=X, p=z)), zs)))/sample_size
    return px

a = 1.0
b = 1.0
pz = stats.beta(a=a, b=b)
px = approximate_marginal_likelihood(pz, sample_size=10000)

解析解を比較してモンテカルロ積分によって積分が実行できることを確認してみます。

f:id:hippy-hikky:20191115004012p:plain — モンテカルロ積分による近似解の算出の様子

この図は、上からサンプルサイズを[100, 500, 1,000]として、それぞれ500回づつ近似解を算出した結果のヒストグラムです。赤線が解析解を示しています。この図から、モンテカルロ積分によって解析解を近似できていることがわかります。また、サンプルサイズを増やすことで近似解のばらつきが少なくなることがわかります。

ただ、結構時間がかかるんですよね。ソースコードをもっと効率化すれば良いのかもしれませんが。

【トップに戻る】

棄却サンプリング

確率分布 $p(\mathbf{z})$ からサンプル( $z^{(l)}$ )を得ることができれば、モンテカルロ法で数値的に積分を計算できることがわかりました。

しかし、そもそも近似ベイズ推論を扱いたいというスタート地点は事後確率分布 $p(\Theta | \mathcal{D})$ を解析的に得ることが難しいので近似的に推論したいということでしたよね。

そこで、このように性質のわかっていない確率分布からサンプル $z^{(l)}$ を得るシンプルな方法として「棄却サンプリング(Rejection Sampling)」があります。

アルゴリズム

まず、サンプルを得たい分布を目標分布( $p(z)$ )と呼びます（これはどのような分布かわからないとする）。次に、目標分布に比例する分布 $\tilde{p}(z)$ を設計します（ $p(z) = \frac{1}{Z_p}\tilde{p}(z)$ ）。 $p(z)$ の性質がわからないのに $\tilde{p}(z)$ をどうやって用意するんだと思った方もいると思いますが（僕がそうでした）、あとで具体例で解説します。

そして、適当な分布（ $q(z)$ , 提案分布）を用意し、 $q(z)$ が $\tilde{p}(x)$ を覆うように定数kを設定します（ $kq(z) \ge \tilde{p}(z)$ ）。

あとは以下の手順でサンプルを生成していきます。

サンプル $z^{(t)}$ を提案分布 $q(z)$ から取得する
サンプル $\tilde{u}$ を区間 $[0,kq(z^{(t)})$ ]の一様分布 $\mathrm{Uni}(u|0,kq(z^{(t)}))$ から取得する
$\tilde{u} > \tilde{p}(z^{(t)})$ ならばサンプル $z^{(t)}$ を棄却する（捨てる）。 $\tilde{u} \le \tilde{p}(z^{(t)})$ ならばサンプル $z^{(t)}$ を受容する（ $q(z)$ に従うサンプルとする）。

提案分布からのサンプルzを確率 $\frac{\tilde{p}(z^{(t)})}{kq(z^{(t)})}$ で採択するということです。

このアルゴリズムを図にしてみたものがこちらです。

提案分布（青線）と $\tilde{p}(z)$ （赤線）の間のサンプルを捨てています。これで $\tilde{p}(z)$ （赤線）の範囲のサンプルを積分して1になるように、つまり、目標分布 $p(z)$ に従うサンプルになるということですね。

問題設定

棄却サンプリングの例として、パラメータ $\theta$ の事後分布 $p(\theta | X)$ を推定する問題を考えます。パラメータ $\theta$ はベルヌーイ分布 $\mathrm{Bern}(x|\theta)$ のパラメータであり、 $X=\{x_1, x_2, \cdots, x_N\}$ はこの分布に従って得られたデータとします。

ベルヌーイ分布なので、具体的にはクジの当たり確率を推定するような問題です。

まず、目標分布 $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(X)$ はパラメータ $\theta$ に関係しない（積分消去する）項なので、比例定数と考えると、 $\tilde{p}(z=\theta) = p(X | \theta)p(\theta)$ となります。ここでは、 $p(X | \theta)$ はベルヌーイ分布です。 $\theta$ は確率なので、 $p(\theta)$ はベータ分布を使うのが良いでしょう。

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

次に、提案分布 $q(\theta)$ です。 $\tilde{p}(\theta)$ の式を眺めると、ベータ分布のパラメータを（a=1, b=1）とすれば最大値が1となることがわかります。そこで、 $q(\theta)$ は[0,1]の一様分布を割り当てることにしてみます（k=1）。

なおこの問題設定の場合、ベータ分布とベルヌーイ分布は共役な関係であるため、事後分布は以下のように解析的に求めることができます。そこで、解析解と近似解の比較もしてみましょう。

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

計算は省略しますが、最終的に以下のベータ分布となります。

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

実装

上記の問題設定に従って、棄却サンプリングによって事後分布 $p(\theta | X)$ を推定してみます。こちらもコードの全体は、実装コード全体にnotebookを貼っていますので参照ください。

なお、真の確率 $\theta$ は0.1とします。観測データ $X$ はN=100個観測したとします。データが未観測でのパラメータ $\theta$ の分布はモンテカルロ積分の例と同様に $\mathrm{Beta}(\theta | a=1, b=1)$ とします。

解析解

解析解は上述の通り、ベータ分布になります。

a_hat = a + np.sum(X)
b_hat = b + N - np.sum(X)
analysis_dist = stats.beta(a=a_hat, b=b_hat)

この分布の密度関数を描いてみると以下のようになります。

f:id:hippy-hikky:20191115135311p:plain — パラメータの事後分布の解析解

近似推定

棄却サンプリングのアルゴリズムに従って近似解を求めれば良いのですが、一つ注意することがあります。

$\tilde{p}(z)$ は尤度 $\prod^{N}_{n=1} \mathrm{Bern}(x_n | \theta)$ と事前分布 $p(\theta)$ の積です。尤度はその形からわかるようにデータ数が増えるほど小さくなります。そのため、提案分布を一様分布、 $k=1$ の設定のままではサンプルが棄却されるばかりでほとんど受容されることがなくなってしまいます。そのため、尤度を適切に調整（定数倍）する必要があります（kを調整するのと一緒）。調整する値としては、二項係数が良いでしょう。二項係数をかけることで、尤度が積分して1になるように正規化されます。なおこれは尤度関数に二項分布を用いることと同じです。

proposal_dist = stats.uniform(loc=0, scale=1)
k = 1.0
a=1.0
b=1.0
p_tilde_lik = stats.bernoulli
p_tilde_prior = stats.beta(a=a, b=b)

n_sample = 500
zs_1 = []
cnt=0
while len(zs_1) < n_sample:
    z = proposal_dist.rvs()
    q_z = proposal_dist.pdf(z)
    u = stats.uniform.rvs(loc=0, scale=k*q_z)
    lik = scipy.special.comb(N, sum(X)) * np.prod(p_tilde_lik.pmf(k=X, p=z))
    prior = p_tilde_prior.pdf(z)
    p_tilde = lik * prior
    if u <= p_tilde:
        zs_1.append(z)

先に求めた解析解(下図赤線)と近似解(下図のヒストグラム)、設定値(下図の縦線)を比較してみます。