CA Reward

Tech Blog

Tech Blogトップに戻る

確率モナドでMCMCサンプラーを実装してみる

2017.04.10

  • このエントリーをはてなブックマークに追加
  • Pocket
しらい
しらいデータサイエンティスト

 データマイニングエンジニアの白井です。今回は個人的に興味を持っている確率モナドを使った話題について紹介したいと思います。

はじめに

 関数型言語において参照透過性が保たれるためには、ある値を代入された関数はいかなる場合でも同じ値を返さなければなりません。しかし、確率的な事象を扱う場合には毎回異なる値を返すため、その性質が壊されてしまいます。その解決策として確率をモナドにして、参照透過な実装とサンプリングの操作を完全に切り離すことが考えられます。今回は連続的な確率においてどのように確率モナドを構築するのかを説明し、統計モデリングで使われる Metropolis-Hastings 法の実装をしてみます。全体を通して scala を使って実装してみたいと思います。
 今回扱う内容は確率モナドの有用性を理解することを目的としており、実用性は重視していません。もし、業務で統計モデリングが必要な場合は Stan などを使う方がよいかと思います。

ベイズの定理を表現する

 ベイズ統計ではデータだけでなく、その背後にあるパラメータも確率的に生成されていると考えます。 例えば、コインを投げたときに表が出る確率を知りたいとします。 正常なコインであればその確率は 1/2 になりますが、それを実際にコインを繰り返し投げることで調べることができます。この場合、ベイズの定理を用いて事後分布を計算することで、表が出る確率自体の分布を得ることができます。
 ベイズの定理は以下のように書かれ、ここではDを実際に観測したデータ、θをコインの表が出る確率を表します。(規格化因子は省いています。)

01bayes-theorem.png

このベイズの定理により、p(θ) という事前分布を与え、データの分布からコインの表が出る確率分布である事後分布 p(θ | D) を得ることができます。ここで尤度である p(D | θ) はパラメータが与えられたときのデータの分布を表します。 もし、これを関数で表現することができたら f: Coin => Distribution[Data] のように書かかれるでしょう。 つまり、Distribution が以下のような flatMap を持っているとすれば、ベイズの定理の枠組みを表現できるかもしれません。

このように flatMap が与えられるということは、Distribution がモナドとして振る舞うことが期待できます。 もちろん、モナド則など満たすべき条件は他にもありますが、モナドを構築する糸口としては十分だと思います。

確率モナド

 離散的な確率分布であれば確率変数はいくつかの決まった値しかとらないためリストを用いて表現できそうですが、連続的な確率を扱う場合は無限要素を持つリストを用意するわけにもいかないので、確率分布を用意しておき、必要なときにサンプリングしなければなりません。しかし、サンプリングを行う場合には取り出される値は毎回異なるため、参照透過性が壊れてしまいます。そのためここでは確率をモナドにすることで、サンプリングを行う操作は計算の論理構造と完全に切り離します。(実装方法は Practical Probabilistic Programming with Monads http://mlg.eng.cam.ac.uk/pub/pdf/SciGhaGor15.pdf の中に書かれているコードを参考にしています。)
 確率分布を表現するために Distribution という型を用意します。これからこの Distribution をモナドにしていきます。

モナドは flatMappoint (PureとかReturnと呼ばれることもある)が実装されている必要があります。ここでは具体的な実装を切り離したいので、以下のようなcase classを用意し計算のプレースホルダーとして使います。Conditional は事前分布と尤度関数を扱うために使うもので、 Primitive は組み込み関数としてある分布関数をラップして確率モナドとして扱えるようにするものです。

確率であることがわかるように別名をつけておきます。

これらのcase classを用いて、モナドを以下のように定義します。モナドはファンクターである必要があるので、mapの定義もしておきます。(scalazの力を借りることにしました。)

最低限これだけ用意してあれば、Distribution がモナドとして振る舞うことができ、flatMapの糖衣構文としてfor文が使えるようになります。今のところ具体的な動作を一切定義していないので、計算した結果つくられるものはPoint, FlatMap, Primitiveなどのcase classがネストしたオブジェクトです。最後にこれらの型に依る具体的な動作を定義したインタープリタを実装します。
 以下のように trait Distributionsampleメソッドをもたせます。sampleメソッドの中にそれぞれの型についてどのような動作をするか定義しておきます。FlatMap の処理がネストして複雑になっているのはstacklessにするためで、Haskellのように末尾再帰を気にしなくてよければもっと簡潔に書くことができます。

以上で必要なものはすべて揃いました。

 Primitive の使い方についてですが、例えば、scala.util.Random にある nextGaussian() という関数からサンプリングしたい場合、以下のようなclassを作って、sampleというメソッドを持たせます。しかし、このままでは Normalは確率モナドの文脈で扱うことができません。そこで Primitive(new Normal(mean, stdDev)) のようにラップすることで、確率モナドとして振る舞うことができ、インタプリタで内で Primitivesample が呼ばれた際に下で定義されている sampleメソッドが呼ばれることになります。

Metropolis-Hastings法

 ここでは上で定義した確率モナドを使って実際にMetropolis-Hastings法の実装をしてみたいと思います。Metropolis-Hastings法のアルゴリズムについてはここでは説明しないので、統計モデリングの本などを参照されたいと思います。実装は以下のようになります。基本的に Metropolis-Hastings 法のアルゴリズムを実装しているだけで、確率分布は全て確率モナドになっている以外特別なことはしていません。

bernoulliというメソッドは以前に説明した Primitive を使って、ベルヌーイ分布からのサンプリングを行います。

また、Prior.prior(d) は提案分布を与えます。やっていることは引数でテストデータ受け取って尤度関数を作ることです。愚直に実装すると、 https://gist.github.com/stenoritama/2c7ca0b5a49ec2b98e0ed6f639110823 のように末尾再帰ではなくなってしまいます。scalaの場合は末尾再帰にしなければ最適化されないので、トランポリンと呼ばれる方法を使います。 https://gist.github.com/stenoritama/ee7630ec7d12aef73b19724394c9ca6e (長くなるのでここだけgistに書きました)

 データは Conditional を使って尤度関数と共に与えます。尤度関数はモデルパラメータを引数にとりデータが既に適用された状態、例えば、y = ax + b のような線形回帰をする場合、尤度をガウス分布だとすると、データ (x, y) は既に適用され、パラメータを引数にとって確率密度を返すようにします。

これを尤度関数として、Conditionalに入れて、全データを使って Condtional がネストした構造を作ります。(foldLeftとか使えば簡単に作れる。)

最も内側にある priorDist は事前分布なのでパラメータ a, b の事前分布を与えます。これもガウス分布だとすると、

のように同時確率を作ります。このようにデータと事前分布を与えることで、priorの中で尤度が計算されます。ここまで実装すれば実際に動かすことができます。

prosterior が事後分布を表すオブジェクトになっていますが、この時点ではまだサンプリングは行われてなく、ただcase classがネストしたものになっています。最後の行の sample を実行して初めて必要なサンプリングが全て行われ、事後分布が数値として得られます。あたかも直接事後分布からサンプリングしているかのように見えます。

 今回は実装方法の説明なのであまり重要ではないですが、実際に y = -0.5 x + 0.3 にガウスノイズを加えたデータを使って計算してみたので、その結果を載せておきます。aが-0.5、bが0.3の周りで振動していることがわかります。

02param_a.png 03param_b.png

おわりに

 長々と説明してきましたが、重要なことは確率モナドを使うことで毎回結果の違うサンプリングという操作を他の参照透過な実装から完全に切り離すことが出来ることです。さらに構造だけを先に作ってしまって最後にサンプリングを行っているため、あたかも事後分布から直接サンプリングしているかのように見えるところは確率分布を直感的に扱える点で面白いと思いました。今回の例はそれほど実用性を重視していませんでしたが、確率モナドを理解するという点で非常に勉強になったと思います。

参考文献