なぜベイズ推定にMCMCなのか?

MCMC (マルコフ連鎖モンテカルロ法) なにも分からん。ベイズ推定を勉強しているといつだって「計算が困難なので MCMC で事後分布を近似しよう!」という話に遭遇する。しかしながら MCMC が何を計算するものなのかよく分からないし、なにが嬉しいのかも分からない....。

そんなあなたに向けて、本記事ではベイズ推定と MCMC について解説します。以下の質問に答えられるようになるのが目標です。

  • MCMC は何を入力として何を出力するものなのか?
  • なぜ MCMCベイズ推論と相性がよいのか?
  • まともにやると困難な計算を、なぜ MCMC では回避できるのか?

なお、本記事では以下の知識があることを想定します。

目次:

例によって、私自身が勉強しながらまとめたものなので間違ってる可能性が大いにあります。特に「すべての MCMC がこうである」って話をする自信はないので、MCMC として代表的な メトロポリス・ヘイスティングス法 の話をします。つまりこの記事の MCMC とはメトロポリスヘイスティングス法のことです。

ベイジアンネットワークと同時分布

まずはじめに、確率モデルの同時分布を、ベイジアンネットワークに基づいて条件付き確率の積として表現できることを説明します。

ベイジアンネットワークとは、確率変数の依存関係を示した有向グラフの一種です。以下の Wikipedia のページ から借りてきたベイジアンネットワークの例を用いて説明します。

f:id:t-keita:20210420221705p:plain:w500

この図には3つの確率変数 RAIN, SPRINKLER, GLASS WET が存在します。RAIN は雨が降ったかどうか、 SPRINKLERスプリンクラーが作動したかどうか、GLASS WET は芝生が濡れているかどうかを示す確率変数です。それぞれの頭文字をとって確率変数  R,  S,  G と示すことにしましょう。

図中に書いてある表は、それぞれの確率変数が取る値とその確率を示しています。たとえば、雨が降ったとき(RAIN = T のとき)にスプリンクラーが作動しない(SPRINKLER = F)確率は 0.99 です。条件付き確率として書くと  \text{Pr}(S = \text{F} | R = \text{T} ) = 0.99 という感じです。これは確率変数 S が確率変数 R に依存しているとも言えます。

ベイジアンネットワークは、このような確率変数同士の依存関係を表現したグラフです。グラフの各ノード(頂点)が確率変数に対応しており、ノード A からノード B へのエッジ(矢印)は確率変数 B が確率変数 A に条件付けられていることを表します。確率変数同士の関係性をグラフとして可視化したものであるとも言えます。

ベイジアンネットワークは、人間が問題を理解しやすくなるだけでなく、確率計算の方法を提供するという点で価値があります。たとえば、先ほどの確率モデルの同時分布  \text{Pr}(R, S, G) は以下のように積の形で計算されます。

 \text{Pr}(R, S, G) = \text{Pr}(G | S, R) \times \text{Pr}(S | R)  \times \text{Pr}(R)

どの確率変数を条件に含めるかは、ベイジアンネットワークのエッジの有無に基づいています。たとえば、確率変数 R から S へのエッジあるので、 \text{Pr}(S | R) となっています。

同時分布を積の形に分解することで、具体的な値の組み合わせが現れる確率を計算できるようになります。たとえば、「雨は降っておらず(RAIN = F)、スプリンクラーが作動し(SPRINKLER = T)、芝が濡れている(GRASS WET=T)」となる確率  \text{Pr}(S = \text{T}, R = \text{F} , G = \text{T}) は以下のように積の形を用いて求められます。


\text{Pr}(R = \text{F}, S = \text{T}, G = \text{T})

=  \text{Pr}(G = \text{T}| S= \text{T}, R= \text{F}) \times \text{Pr}(S= \text{T} | R= \text{F})  \times \text{Pr}(R= \text{F})
 = 0.9 \times 0.4 \times 0.8
 = 0.288

ここまでの話で重要なポイントをおさえておきます。それは「確率モデルがベイジアンネットワークとして表現できるとき、同時分布において具体的な組み合わせが起こる確率を計算できる」ということです。そして、ベイズ推定で使用する確率モデルは一般にベイジアンネットワークとして表現できます。つまり、ベイズ推定では、同時分布において具体的な実現値が起こる確率を簡単に計算できます。 詳しくはあとで述べますが、MCMC ではこの性質を活用します。

ベイズ推定と事後分布の計算

ここからは本題のひとつであるベイズ推定について説明します。ベイズ推定におけるベイジアンネットワークでは、確率変数の値が観測済みであるかどうかが区別されるものとします。観測とは、芝生が濡れているとか、サイコロの 3 の目が出たとか、未知だった情報が明らかになったことを言います。たとえば、以下のようなベイジアンネットワークを考えましょう。観測済みのノードは灰色に塗られています。確率変数  X_1,  X_2 は具体的な実現値  x_1,  x_2 が観測された状況です。

f:id:t-keita:20210421013413p:plain:w400

一般に、ベイズ推定におけるベイジアンネットワークには 2 種類の確率変数が存在することになります。1つは未観測の確率変数  \theta_1 , \dots, \theta_n であり、もう1つは観測済みの確率変数  X_1 , \dots, X_m です。ベイズ推定では、それぞれの確率変数について以下のように確率分布を与えます。

  • 未観測の確率変数  \theta_1 , \dots, \theta_n それぞれに対して事前分布 \text{Pr}(\theta_1) , \dots, \text{Pr}(\theta_n) を人間が与える
  • 観測済みの確率変数  X_1 , \dots, X_m それぞれに対して、確率1で観測値  x_1 , \dots, x_m を返すような確率分布を与える(厳密には ディラックのデルタ関数 を用いると確率の公理を満たすように確率分布を定義できます)

そしてベイズ推定では、これらの確率分布をもとに未観測の確率変数の事後分布 \text{Pr}(\theta_1 , \dots, \theta_n | X_1 = x_1 , \dots, X_m = x_m) を求めます。具体的には、条件付き確率の定義から以下のように事後分布を計算できます。


\text{Pr}(\theta_1 , \dots, \theta_n | X_1 = x_1 , \dots, X_m = x_m)
= \frac{\text{Pr}(\theta_1 , \dots, \theta_n, X_1 = x_1 , \dots, X_m = x_m)}{\text{Pr}(X_1 = x_1 , \dots, X_m = x_m)}

ここで、右辺分子の  \text{Pr}(\theta_1 , \dots, \theta_n, X_1 = x_1 , \dots, X_m = x_m) はすべての確率変数の同時分布であり、ベイジアンネットワークをもとに条件付き確率の積へ変形できるのでした。そして、具体的な実現値が起こる確率を簡単に計算できるのがポイントでした。

一方で、右辺分母の  \text{Pr}(X_1 = x_1 , \dots, X_m = x_m) を計算するには、未観測の確率変数  \theta_1 , \dots, \theta_n の組み合わせごとに確率を求め、その合計値を計算する必要があります(要するに 周辺化 の計算です)。たとえばスプリンクラーの例において「芝生が濡れている(GLASS WET = T)」確率を計算しようとすると、スプリンクラーと雨の状態のすべての組み合わせ 2 × 2 = 4 通りについて芝生が濡れている確率を求める必要があります。4 通りだけなら大丈夫ですが、未観測の確率変数の数 n が大きいときにはその組み合わせの数が膨大になり、コンピュータを使っても私たちが生きている間には計算が終わりません(つまり 組み合わせ爆発 を起こすということです)。それゆえ、事後分布を直接計算することは現実的ではありません。さらに、確率変数が連続値を取る場合は積分計算が必要になり、離散値の場合と同等かそれ以上に計算が難しいです。いずれにせよ、未観測の確率変数の取りうる組み合わせを調べつくすことが困難なのです。

ここでひとつ重要なポイントがあります。それは、分母の  \text{Pr}(X_1 = x_1 , \dots, X_m = x_m) の値を求めることは難しくても、この値が定数であることは確定しているということです。そのため、事後分布(左辺)は同時分布(右辺の分子)に比例します。 詳しくはあとで述べますが、MCMC ではこの性質を活用します。

MCMC は何を入力として何を出力するのか?

いったんベイズ推定のことは忘れて、MCMC が何者なのかを説明します。MCMC の目的は確率分布  \text{Pr}(Z) をコンピュータシミュレーションにより近似することです。

MCMC の入出力を以下に示します。これが明示的に書かれない解説が多いので理解が曖昧になりがちです。

  • 入力:関数  f(z)。ただし、以下の要件を満たすもの。
    1.  f(z) の値が確率  \text{Pr}(Z = z) に比例する。
    2.  f(z) の値を簡単に計算できる。
  • 出力:確率分布  \text{Pr}(Z) に従う実現値(サンプル)の列。

たとえば、近似したい確率分布  \text{Pr}(Z) がサイコロの出目であるとすると、MCMC の出力は 1 から 6 までの実現値の列(たとえば [2, 4, 1, 3, 6, 5, 1, 4,...])となります。サンプル数を 600 とすると、1 から 6 までのそれぞれの値が約 100 回ずつ出現することが期待できるわけです。そして、100 / 600 = 1 / 6 なので、 \text{Pr}(Z) が 1 から 6 までの値を確率 1 / 6 で取るような一様分布であることが推測できるわけです。

重要なポイントは MCMC の入力が関数  f(z) であることです。MCMC の入力は  \text{Pr}(Z) ではないことに注意してください。確率分布  \text{Pr}(Z) そのものが分からなくても、代わりに計算可能な  f(z) さえ与えられれば  \text{Pr}(Z) を近似できるのです。

MCMC の仕組みを簡単に説明します。特に、代表的な MCMC アルゴリズムであるメトロポリスヘイスティングス法について述べます。

まずは初期値として適当なサンプル  z^{(1)} から開始します。このとき時刻  t = 1 とします。つぎに、遷移先の候補として z^{(t)} の近くのサンプル  z^{\ast} をランダムに得ます。ここで、実際に  z^{\ast} に遷移するかどうかは確率的に決まります。その確率は  f(z^{\ast}) f(z^{(t)}) の値から計算されます。(これらの値が簡単に計算できることは MCMC の入力の要件「 f(z) の値を簡単に計算できる」によって保証されています。) z^{\ast} へ遷移する場合は  z^{(t+1)} = z^{\ast} として、遷移しない場合は  z^{(t+1)} = z^{(t)} とします。あとは、この遷移の操作を十分な数のサンプル列  z^{(1)}, z^{(2)}, \dots, z^{(N)} が集まるまで繰り返すだけです。このサンプル列が MCMC の出力であり、確率分布  \text{Pr}(Z) に従うことが期待できます。

サンプルが取りうる値の集合をサンプル空間と呼ぶことにすると、MCMC ではサンプル空間を行ったり来たりして探索することになります。MCMC の重要なポイントは、サンプル空間を調べつくすことなく、確率分布を近似できるということです。地味に聞こえますが、MCMC は 20 世紀に発見されたアルゴリズムの中でもっとも重要であると言われるくらいすごいものなのです。

MCMC で事後分布を近似する

話をベイズ推定に戻します。おさらいですが、事後分布は以下のように計算できるのでした。


\text{Pr}(\theta_1 , \dots, \theta_n | X_1 = x_1 , \dots, X_m = x_m)
= \frac{\text{Pr}(\theta_1 , \dots, \theta_n, X_1 = x_1 , \dots, X_m = x_m)}{k}

右辺分母は計算は困難であるものの定数であることは分かっていたので  k で置き換えています。ベイズ推定の目的は左辺の事後分布を求めることです。そして、事後分布は右辺分子の同時分布に比例するのでした。

ここで、本記事の伏線を一気に回収します。ズバリ、ベイズ推定の事後分布を MCMC で近似します。具体的には、MCMC で近似したい確率分布  \text{Pr}(Z) を事後分布  \text{Pr}(\theta_1 , \dots, \theta_n | X_1 = x_1 , \dots, X_m = x_m) として、MCMC の入力となる関数  f(z) を同時分布  \text{Pr}(\theta_1 , \dots, \theta_n, X_1 = x_1 , \dots, X_m = x_m) とします。この同時分布は関数  f(z) の要件を満たしています。以下の2点を満たしているからです。

  1. 同時分布は事後分布に比例する。
  2. 同時分布はベイジアンネットワークによって積の形に分解できるため、具体的な実現値が起こる確率を簡単に計算できる。

これにて MCMC でサンプリングした結果が事後分布の近似になっていることが分かったわけです。めでたしめでたし。

MCMC は素朴にやると困難な計算をなぜ回避できるのかを整理しましょう。事後分布を正確に計算しようとすると、どうしてもサンプル空間を調べつくす必要があったわけです。具体的には、サンプル空間に含まれるそれぞれのサンプルに対して  X_1 = x_1 , \dots, X_m = x_m となるときの確率を求めて合計する必要がありました。

一方で、MCMC ではサンプル空間を調べつくすことなく、決められた個数のサンプルによって事後分布を近似できるのでした。つまり、MCMC によってサンプル空間を調べつくすコストが回避されるわけです。裏を返すと MCMC では計算の手抜きをしているため、事後分布を正確に近似できない可能性もあります。このあたりの詳細については世の中の書籍を参照ください。

まとめ

この記事ではベイズ推定と MCMC の関係について整理しました。私がざっくり調べた感じ、ベイズ推定や MCMC の説明自体は世の中にたくさんあるのですが、ベイズ推定と MCMC の相性がよい理由について説明したものは多くはありませんでした。そこでこの記事を書くに至りました。

整理のため今回の内容を1枚絵にしてみました。

f:id:t-keita:20210422010637p:plain

ベイズ推定の確率モデルは、たった1つの条件付き確率の式として書けます。(多くのベイズ推定の解説ではベイズの定理の式が用いられますが、個人的には同時分布をそのまま残した方がスッキリしてよいと思ってます。)MCMC は「工場」であり、準備として同時分布を受け取ります。MCMC 工場では一定数のサンプルが製造されるのですが、そのサンプルを集めたものが事後分布の近似になっていることが期待できます。

ちなみに、確率モデルと観測値と事前分布を記述するだけで、MCMC によって事後分布を自動的に計算してくれるソフトウェアとして StanPyMC3 などがあります。これらについての記事もそのうち書きたいと思います。

参考文献