マルコフ連鎖モンテカルロ法を実装してみよう

"Chain Bridge, Budapest" by szeke

約2ヶ月ぶりということで、ご無沙汰しております。書きたいネタというのは結構あって書こう書こうとは思っているんですが、なにより書くとなるとあれもこれも伝えなきゃいけないみたいな思いになって、結局分量の多さから諦めてしまうというのが結構続いています。もう少し気を張らずに更新していきたいものです。

さて、最近の自分はマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo, MCMC)についておりました。何にも知らない方ににとってはよく分からないかもしれませんが、この手法はマルコフ連鎖が持つ簡明さとモンテカルロ法が持つ実用性が合わさった凄まじい手法なんです。そしてなによりとてもエレガント。

とりあえず知らない人のためにてきとーに解説しますと、与えられた適当な関数から確率変数をサンプリングするための公式です。べつにこれじゃなくても確率変数のサンプリングにはいろんなアルゴリズムがあるわけですが、これを使うと関数の計算コストがやけに高い場合だとか、ピーク部分を集中的にサンプリングしてくれたりとかで結構便利で、金融や工学、物理の分野でも非常に役立っております。計算物理をやる人にとっては必須に近い手法といっても過言ではありません。

そんなMCMCですが、実際に提案されている手法には様々なものがあります。不変分布から遷移確率を逆算するため、自由度が増えるから遷移確率は一意じゃなくなってしまうんです。今回は最も有名なMetropolis-Hastingsアルゴリズム(*1)を用いて、任意の関数に従って乱数のサンプリングを行うコードを練習がてらpythonで書いてみました。

#!/usr/bin/env python
#-*- coding:utf-8 -*-

import sys
import math
import random

def main():
   f = lambda x: math.exp(-x*x / 10.0)
   for x, fx in metropolis(f):
       print x, fx

def metropolis(func, start=0.0, delta=2.0, burn_in=100, seed=1):
   metro_iter = _metropolis(func, start, delta, seed)
   for i, value in enumerate(metro_iter):
       if i > burn_in:
           yield value
           break
   for value in metro_iter:
       yield value

def _metropolis(func, start, delta, seed):
   random.seed(seed)
   initial_x = start
   initial_fx = func(initial_x)
   proposed_x, proposed_fx = 0.0, 0.0
   while True:
       proposed_x = initial_x + random.uniform(-delta, delta)
       proposed_fx = func(proposed_x)
       ratio = proposed_fx / initial_fx
       if ratio > 1.0 or ratio > random.random():
           initial_x, initial_fx = proposed_x, proposed_fx
           yield proposed_x, proposed_fx
       else:
           yield initial_x, initial_fx

if __name__ == '__main__':
   sys.exit(main())

うーん・・・きれいだ(*2)・・・

こいつのすごいところは、アルゴリズム自体は非常にシンプルで実装しやすいというところにあります。行数も上を見ると分かるのですがとっても少ないし、コア部分はたった25行くらいで構成されています。

なのに使っている学問は結構複雑で、定常分布から遷移行列を導出するために詳細釣り合いの条件からMetropolis的制約を用いて導出、サンプリングを行っております。このアルゴリズム自体は結構知ってる人が多いかもしれませんが、実際に理論方面で理解している人はどうなんでしょう。

んで、正規分布のガウス関数からサンプリングを行った結果がこれ。

赤線が元の確率分布。きちんとサンプリングが行われているようです。今回は正規分布で計算を行いましたが、非負の関数であればなんでも構いません。

あとはこれを研究テーマへと落とし込むことができれば、とりあえず一区切りつけるかなって感じです。

*1: 提案確率関数が対称であるという制約を持つのがMetropolis。それを非対称の場合にも拡張したのがHastings。個人的にはMetropolis-HastingよりMetropolisのほうが綺麗で好きです。というか提案確率密度関数が非対称じゃないといけないようなシチュエーションがよく分からない
*2: yieldとジェネレータって便利ですよね。計算速度も効率的だし、慣れてくるとすべて遅延評価で行わせたいくらい