サンプリングによる近似ベイズ推論その3（MCMC:ギブスサンプリング）

【概要】

ベイズ推論について実装して理解するシリーズ
今回は、MCMC アルゴリズムの一つであるギブスサンプリングです
ギブスサンプリングによって線形回帰（ベイズ線形回帰）を近似推論してみました

【目次】

はじめに
近似ベイズ推論
ギブスサンプリング (Gibbs Sampling)
ギブスサンプリングのアルゴリズム
ギブスサンプリングによる確率分布の近似推論の実装
- 2次元ガウス分布の近似推論
  - 近似解
- ベイズ線形回帰の近似推論
実装コード全体
おわりに
参考文献

はじめに

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

本記事ではベイズ推論における近似推論について扱います。ベイズ推論では「MCMC(Markov Chain Monte Carlo)法」などのサンプリングに基づくアルゴリズムと「変分推論」という代表的な近似推論手法があります。

サンプリングに基づく近似ベイズ推論(MCMC)について、過去2回と合わせて今回が3回目となります。
その1では、基本的な積分の近似計算である「モンテカルロ積分」と最も単純なサンプリングアルゴリズムである「棄却サンプリング」を実装しました。その2では、MCMC アルゴリズムの一種である「メトロポリス法」を実装しました。その3となる本記事では、MCMC アルゴリズムの一種であるギブスサンプリングを実装してみます。

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

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

作者:須山敦志
出版社/メーカー: 講談社
発売日: 2019/08/08
メディア: 単行本

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

【トップに戻る】

近似ベイズ推論

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

ベイズ推論においては、推定対象のパラメータを含めて全て確率変数と考えます。パラメータ $\Theta$ が従う確率分布 $p(\Theta)$ を観測データ $\mathcal{D}$ を使って更新することが、ベイズ推論における「学習（パラメータ推論）」になります。そのようにして得られる、データ $\mathcal{D}$ を観測した後でのパラメータ $\Theta$ の確率分布 $p(\Theta | \mathcal{D})$ は事後確率分布と呼ばれています。

事後確率分布 $p(\Theta | \mathcal{D})$ はモデル（確率変数の同時分布）を展開することで導出できます*1。しかし、複雑なモデルであったり、モデルを構成する確率分布（ $p(\mathcal{D} | \Theta), p(\Theta)$ ）が共役関係にない場合には、 $p(\Theta | \mathcal{D})$ を解析的に求めることができません。
そこで、近似推論が必要になります。

【トップに戻る】

ギブスサンプリング (Gibbs Sampling)

ギブスサンプリングは、メトロポリス法(前回記事)と同様に、性質のわからない確率分布 $p(\mathbf{Z})$ からサンプル $\{z^{(1)}, z^{(2)}, \cdots, z^{(L)}\}$ を得るためのアルゴリズムです(Lはサンプル数)。 $p(\mathbf{Z})$ からサンプルを得ることができれば、そのサンプルの統計的性質を調べることで $p(\mathbf{Z})$ の性質がわかるのでした。

ギブスサンプリングは、以下の条件を満たす問題の場合に適用でき、メトロポリス法などと比較して高次元のモデルから効率的にサンプリングを行うことができる手法です。

複数の確率変数()を持ったモデルである
- $\mathbf{Z}_{i}$ は $\mathbf{Z}$ をM個の部分集合( $\mathbf{Z}_{1}, \mathbf{Z}_{2}, \cdots , \mathbf{Z}_M$ )に分解したもの
をサンプルすることが容易な条件付き分布を構築することができる
- $\mathbf{Z}_{-i}$ は $\mathbf{Z}_{i}$ を除いた確率変数の集合とする

適用できるか否かには条件がありますが、確率変数 $\mathbf{Z}$ が高次元な場合などで有効な手法ということです*2。

ひとつわからないところがありまして、条件付き分布が既知の分布ではない場合にもサンプルできるものなんでしょうか？？

ギブスサンプリングのアルゴリズム

サンプルを得たい目標分布 :
- 確率変数 $\mathbf{Z}$ : $\mathbf{Z} = \{\mathbf{Z}_{1}, \mathbf{Z}_{2}, \cdots , \mathbf{Z}_M\}$
確率変数 $\mathbf{Z}$ の部分集合 $\mathbf{Z}_i$ の条件付き分布 : $p(\mathbf{Z}_{i} | \mathbf{Z}_{-i})$ (全てのiについて)

条件付き分布 $p(\mathbf{Z}_{i} | \mathbf{Z}_{-i})$ が肝です。これがガウス分布などの既知の分布として構築できれば、アルゴリズムとしては単純です。
サンプル $\mathbf{z}^{(l)}$ は以下の手順で取得できます。

$\mathbf{Z} = \{z_{1}^{(0)}, \cdots, z_{M}^{(0)} \}$ を適当に初期化
for l=1 to L :
1. for i=1 to M :
  1. $\mathbf{Z}_{*} \sim p(\mathbf{Z}_{i} | \mathbf{Z}_{-i})$
  2. $z_{i}^{(l)} = \mathbf{Z}_{*}$

メトロポリス法と異なりサンプルした点が必ず受容されるため、ギブスサンプリングが使えればメトロポリス法などと比較すると効率良くサンプルを得ることができます。詳しくは、「ベイズ深層学習」や「PRML(下)」などを参照ください。

【トップに戻る】

ギブスサンプリングによる確率分布の近似推論の実装

ここでは、ギブスサンプリングの実装例として二つの例を取り挙げます。

2次元のガウス分布からのサンプル
ベイズ線形回帰の重みパラメータの事後分布からのサンプル

実装コードの全体は、実装コード全体にnotebookを貼っていますので参照ください。

2次元ガウス分布の近似推論

2次元ガウス分布自体はよく知られた確率分布で、直接サンプルを得ることは容易です*3。ここではあくまで適用例として、ギブスサンプリングを用いて2次元ガウス分布からのサンプルを1次元ガウス分布から取得してみます。

2次元ガウス分布。

f:id:hippy-hikky:20191203140751p:plain:w300

ここで $\mu_i$ はi次元の平均値、 $\Sigma_{ij}$ はi次元とj次元の共分散を表します。

ギブスサンプリングを実行するためには、確率変数を部分集合に分けた条件付き分布を導出する必要があります。

${ \displaystyle p \left( x_{1} \mid x_{2}=x_{2}^{(l)} \right), p \left( x_{2} \mid x_{1}=x_{1}^{(l)} \right) }$

はじめに $x_1$ についての条件付き分布 $p \left( x_{1} | x_{2}=x_{2}^{(l)} \right)$ を導出します。確率の基本性質を使って以下の様に条件付き分布を考えます。

${ \displaystyle \begin{eqnarray} p \left( x_{1} | x_{2} \right) = \frac{p(x_1, x_2)}{p(x_2)} \propto p(x_1, x_2) \end{eqnarray} }$

ここで、 $x_2$ は固定値なので、展開して地道にまとめていくと以下の式となります。計算の詳細については省略します。詳細は「PRML(上)」や「ガウス過程と機械学習」などを参照ください。

${ \displaystyle \begin{eqnarray} p \left( x_{1} | x_{2} \right) = \mathcal{N} \left( \mu_{1} + \Sigma_{12} \Sigma_{22}^{-1} \left( x_{2} - \mu_{2} \right), \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21} \right) \end{eqnarray} }$

$x_2$ についての条件付き分布 $p(x_2 | x_1)$ については、上式を対称にして以下の通りとなります。

${ \displaystyle \begin{eqnarray} p \left( x_{2} | x_{1} \right) = \mathcal{N} \left( \mu_{2} + \Sigma_{21} \Sigma_{11}^{-1} \left( x_{1} - \mu_{1} \right), \Sigma_{22} - \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12} \right) \end{eqnarray} }$

これで条件付き分布を導出することができたので、あとはギブスサンプリングのアルゴリズムに従って、 $x_1$ と $x_2$ をそれぞれサンプルします。

近似解

ギブスサンプリングにより2次元ガウス分布からサンプルを得るためのコードは実装コード全体のnotebookのセル番号5にあります。なおMCMC アルゴリズムでは、定常分布に収束するまでの期間のサンプルは破棄するのが普通ですが、今回その部分は書いていません。

サンプルが2次元ガウス分布からのサンプルになっているか確認します。

f:id:hippy-hikky:20191203223901p:plain — ギブスサンプリングによるサンプル結果．オレンジの点がギブスサンプルによるりサンプルした点．青点がnumpy.randomからサンプルした点．

numpy.randomでサンプルした点とギブスサンプルによってサンプルした点群が重なっていることがわかります。1次元ガウス分布を使って2次元ガウス分布からのサンプルを得ることに成功したようです。
また、最初の50点の動きを赤線で示しています。ギブスサンプリングでは、変数を一つずつ固定するため、カクカクした動きになっていることがわかります。

ベイズ線形回帰の近似推論

前節では2次元ガウス分布からサンプルを得ることができたわけですが、この例だけでは何も面白くないですよね*4。ガウス分布はよく知られた分布ですし、その平均や分散のパラメータは予め与えられており、推論のプロセスは入っていません。

ギブスサンプリングを使うモチベーションは、確率モデルのパラメータをベイズ推論する際の近似解法として使いたいというものでした。そこで本節では、線形回帰（いわゆる「ベイズ線形回帰」）のパラメータ推論にギブスサンプリングを適用してみます*5。

モデル

線形回帰モデルのグラフィカルモデルは以下の通りです。 f:id:hippy-hikky:20191204105503p:plain

N個のデータ $(X,Y)$ のペアが与えられ、変数Xと重みパラメータWに変数Yが依存する構造です。重みパラメータWが未知のため、このパラメータをデータ $(X,Y)$ から推論したいわけです。
このモデルの同時分布は以下の通り。

${ \displaystyle \begin{eqnarray} p(\mathbf{w}, \mathbf{Y} | \mathbf{X}) = p(\mathbf{w})p(\mathbf{Y} | \mathbf{X}, \mathbf{w}) = p(\mathbf{w}) \prod^{N}_{n=1} p(y_n | \mathbf{x}_{n}, \mathbf{w}) \end{eqnarray} }$

今回は、変数Xを基底関数 $\mathbf{\phi}(x)=\{\phi_1(x), \phi_2(x), \cdots, \phi_M(x)\}^{T}$ で変換し*6、重みWと線型結合するモデルを考えます*7。

${ \displaystyle \begin{eqnarray} p(y_n | \mathbf{x}_{n}, \mathbf{w}) = \mathcal{N}(y_{n} | \mathbf{w}^{T}\phi(\mathbf{x}_n), \sigma^{2}_{y}) \end{eqnarray} }$

ここで簡単のために、分散 $\sigma^{2}_{y}$ は固定のパラメータとします。

重み $\mathbf{w}$ は未知の確率変数であり、これにも何らかの確率分布を仮定する必要があります。重み $\mathbf{w}$ は正負の実数をとりますが、極端に大きい値でないことが好ましいです*8。そこで、重み $\mathbf{w}$ はガウス分布に従うとし、事前分布として平均0で分散 $\sigma^{2}_{w}$ を与えることにします。なお、分散 $\sigma^{2}_{w}$ は固定のパラメータとします。

${ \displaystyle \begin{eqnarray} p(\mathbf{w}) = \mathcal{N}(\mathbf{w} | \mathcal{0}, \sigma^{2}_{w}\mathbf{I}) \end{eqnarray} }$

サンプルデータ

今回は、2次の多項式(2次関数)を真の関数とします。重みは[0.0, -2.0, 0.5]としました*9。この関数から20点をランダムにサンプルします。真の関数(青の実線)とサンプルデータ(オレンジ点)は以下の通りです。

f:id:hippy-hikky:20191204233212p:plain — 真の関数(青実線)とサンプルデータ(オレンジ点)．横軸がxで縦軸がy．

解析解

ベイズ線形回帰の解析解については、以前記事を書きました。

learning-with-machine.hatenablog.com

知りたいのは、未知の重みパラメータwの事後分布 $p(\mathbf{w} | Y, X)$ です。この事後分布を確率の基本性質に従って展開します。

${ \displaystyle \begin{eqnarray} p(\mathbf{w} | Y, X) \propto p(\mathbf{w}) p(\mathbf{Y} | \mathbf{X}, \mathbf{w}) = \mathcal{N}(\mathbf{w} | \mathcal{0}, \sigma^{2}_{w}\mathbf{I}) \prod^{N}_{n=1} \mathcal{N}(y_n | \mathbf{w}^{T}\phi\left(\mathbf{x}_{n}\right), \sigma^{2}_{y}) \end{eqnarray} }$

計算は省略しますが、最終的には以下のガウス分布になることがわかっています。

${ \displaystyle p(\mathbf{w} | \mathbf{Y}, \mathbf{X}) = \mathcal{N}(\mathbf{w} | \hat{\boldsymbol{\mu}}, \hat{\mathbf{\Sigma}}) }$

${ \displaystyle \ \ \hat{\mathbf{\Sigma}}^{-1} = \sigma_{y}^{-2} \sum_{n=1}^{N} \phi\left(\mathbf{x}_{n}\right) \phi\left(\mathbf{x}_{n}\right)^{T}+\sigma_{w}^{-2} \mathbf{I} \\ \ \ \hat{\mu} = \hat{\Sigma} \sigma_{y}^{-2} \sum_{n=1}^{N} y_{n} \phi\left(\mathbf{x}_{n}\right) }$

データを当てはめ、重みパラメータwの事後分布 $p(\mathbf{w} | Y, X)$ を推論した結果は以下の通りです。

f:id:hippy-hikky:20191204233546p:plain — 重みパラメータwの事後分布．上から0次, 1次, 2次の項の重みパラメータの事後分布．

この事後分布から関数（重みパラメータ $\mathbf{w}$ ）を20本サンプルしてきた結果を以下に示します。グレーの線がサンプルしてきた関数です。

f:id:hippy-hikky:20191204234106p:plain — ベイズ線形回帰の解析解．グレーの線は事後分布からサンプルした関数．

うまく設定値(真の関数)を推論できていることがわかりますね。

これらの実装例は、実装コード全体に貼ってあるnotebookの9~12セルです。

近似解

次に、重みパラメータ $\mathbf{w}$ の事後分布をギブスサンプリングによって近似推論してみます。

ギブスサンプリングによって近似推論するためには、確率変数の部分集合の条件付き分布が必要でした。ここでは、重みパラメータ $\mathbf{w}$ が未知の確率変数です。また、 $\mathbf{w}$ は0次から2次までの3次元のベクトルでした。そこで、 $\mathbf{w}$ の各次元毎の条件付き分布 $p(w_{i} | w_{-i}, Y, X)$ を導出してみます。なお、 $w_{-i}$ とはi次元を除いたwのベクトルを意味します。

f:id:hippy-hikky:20191205000051p:plain:w450

ここで、 $\left( w_i \phi_i(x_n) + w_{-i} \phi_{-i}(x_n) \right)$ については、ギブスサンプリングにおける毎回のサンプリング手順では $w_{-i}$ が与えられているものとするため、未知の確率変数は $w_i$ だけです。そのため、この事後分布は1次元のガウス分布になります。計算過程について細かくは省略しますが、本節の最後に手書きのメモを貼っておきますので、興味のある方は参考にしてください（汚い字でごめんなさい）。

f:id:hippy-hikky:20191205002448p:plain:w400

1次元のガウス分布のため、解析解の導出と比較してちょっとだけ簡単になったような気がします。

実装コードの詳細は実装コード全体に貼ってあるnotebookを参照ください。セル番号の13,14が該当します。事後分布の導出さえできてしまえば、アルゴリズムは非常に単純です。

サンプルのトレースを確認してみます。ちょっとトレンドを持っているように見えますが、定常分布からサンプルできているとしましょう。

f:id:hippy-hikky:20191205003109p:plain — サンプルのトレース．サンプル数は1000だが、最初の300点は捨てている(burn-in)．

事後分布は次のようになります。解析解の事後分布とだいたい同じようなものであることがわかります。

f:id:hippy-hikky:20191205003412p:plain — ギブスサンプリングによる重みパラメータの事後分布．上から, 0次，1次，2次の項のパラメータ．

最後に、この近似推論した事後分布から関数をサンプルしてみます。これも、解析解と比較すると概ね同じものであることがわかります。

f:id:hippy-hikky:20191205003651p:plain — 近似事後分布からの関数のサンプル．グレーの線がサンプルした関数．

線形回帰問題の条件付き分布の導出過程メモ

f:id:hippy-hikky:20191205001538j:plain:w250 f:id:hippy-hikky:20191205001558j:plain:w250

【トップに戻る】

実装コード全体

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

【トップに戻る】

おわりに

ということで、サンプリングをベースにした実用的な近似推論アルゴリズムの一つ、ギブスサンプリングを実装してみました。

MCMC アルゴリズムの中でも、前回実装したメトロポリス法と比較すると全てのサンプルが受容されるため、かなり効率の良いアルゴリズムであることはわかりました。実際、計算時間もだいぶ速かったです。しかし、条件付き分布を導出するのが厳しい場面は少なくないんじゃないかなというのが素直な感想です。今回のように単純な線形回帰モデルであれば条件付き分布を求めるのも容易ですが。モデルを構築した際に、うまくパラメータが独立になるように工夫することが必要になるんじゃないかなと思います。

この辺りはまだまだ勉強中なので、ご指摘いただけると助かります。

【トップに戻る】