Gauche:MetropolisProceduralModeling:MCMC-1D

Gauche:MetropolisProceduralModeling:MCMC-1D

(Shiro 2012/01/03 10:04:38 UTC): 「正月にMetropolis Procedural Modelingを基礎から理解してみる」その1。


(この項は主に、Christophe Andriew, Nando De Freitas, Arnaud Douchet, Michael I. Jordan, "An Introduction to MCMC for Machine Learning", Machine Learning, 50, 5--43, 2003. を読んで自分なりに理解したところをメモしている。)

モンテカルロ法で円周率を求めてみたことがある人は多いだろう。 左下が(0,0)、右上が(1,1)である正方形の中の一点をランダムに選び、 それが(0,0)を中心とする半径1の円内に入っているかどうかを調べる。 十分に多くのサンプルを取れば、円内に入った点の割合は、 四分円と正方形の面積比π/4に近づく。

この場合、サンプルは正方形の中から一様に取れば良かった。 けれどももっと一般的な問題では、サンプルの得られる確率が一様であるとは限らない。

例えば偏りのあるサイコロで、 目の出る確率が p(1)=0.5, p(2)=p(3)=p(4)=p(5)=p(6)=0.1 (ひでぇサイコロだ) だとしたら? このケースはまだ簡単だ。{0,1}の範囲で一様にサンプリングして、 0.5以下なら出目は1、0.5〜0.6なら2、0.6〜0.7なら3、等とすればいい。

ではサンプルが離散的じゃなく連続的だったら? ちょいとややこしくなるけど、 同じ原理を適用するなら、xが得られる確率密度関数がp(x)で与えられてるとしたら、 まず{0,1}の範囲内で一様にサンプリングしてある値uを得て、 u = ∫_{-∞}^{X} p(x)dx になるようなX、がサンプル値だ。

でも、p(x)が正規分布のような素直な関数ならいいけれど、一般の関数だとこの積分が 簡単には求まらないかもしれない。xが多次元になったらもっとややこしい。 (何か方法があるのかもしれないけど、Shiroはとりあえずお手上げ。)

で、こういう「確率分布p(x)は計算できるんだけど、 それに従ってランダムなサンプリングが欲しい」って時に使えるのがMCMC。

直感的な理解としては、サンプルを取る場所xをマルコフ過程でもって 試しにふらふら動かして、p(x)に照らしてより「ありそう」な時にその動きを 認める、ってことを繰り返してると、なんとなくxはp(x)の高いところを うろうろしがちになるよね、って感じ。ただ、最急降下法みたいな数値計算と 違って、時々離れたところにぽーんと動いたりするんで、 局所最適にハマっちゃうことが避けられる、と思えばいい。 これがうまくいくための条件とか、一般的な計算式は上の参考文献に出ている。 特に広く応用が効くMCMCアルゴリズムが、Metropolis-Hastingsのもの。以降MHと表記。

プログラマとしては数式を眺めてるよりコードの方が理解が速いんで、 上の参考文献の3.1節に出ている例を実装してみた。

mcmc-util.scmは他でも使う関数をいくつか定義している。今使うのは標準正規分布関数φと、 サンプリング関数dpick。後者は、 平均μ、分散σな正規分布で、minからmaxまでstepごとの離散的なサンプルを取る。 上のサイコロの例のような求め方をしているのであまり効率は良くない。

んで、mcmc1.scmが今回の実験コード。

さて今回、ターゲットとする分布はこんな形。 (この関数は全区間で積分して1にはならないので、「確率密度関数」とは言えないのだけれど、 MHはこういう正規化されてない関数を与えてもokで、これに比例した確率でサンプルを出してくれる。)

(define (p x)
  (+ (* 0.3 (exp (* -0.2 (expt x 2))))
     (* 0.7 (exp (* -0.2 (expt (- x 10) 2))))))

グラフにするとこんな感じ。

[image]

こっからXをサンプリングしたい。こいつは1次元なんで数値積分すりゃ直接求められるんだけど、 MCMCを使ってみよう。

今回使うマルコフ過程は、「現在のサンプル位置がxの時、次のサンプル位置x*は 正規分布N(x, 10^2)に従いますよ」というもの。pick1は現在のサンプル位置を与えると、 この分布に従ってx*の候補を一つ返してくる。

;; Pick a sample x* by proposal distribution - using N(x, 100)
(define (pick1 x) (dpick x 10 -10 20 0.5))

で次がMHの肝。

;; Calculate posteriori probability q(x*|x)
(define (q x* x) (/ (φ (/ (- x* x) 10.)) 10.))

;; Single MH step
(define (step x)
  (let ([u  (random-real)]
        [x* (pick1 x)])
    (if (< u (min 1 (/ (* (p x*) (q x* x)) (* (p x) (q x x*)))))
      x*
      x)))

[0,1]の一様分布からひとつサンプリングして、それが (min 1 (/ (* (p x*) (q x* x)) (* (p x) (q x x*)))) より 小さければx*への移動を認め、そうでなければxに留まる。

これを初期値x0からスタートしてN回繰り返し、xのヒストグラムをハッシュテーブルに記録してく。

;; Start from x0, run N steps of markov chain
(define (run x0 N)
  (let1 results (make-hash-table 'eqv?)
    (do ([i 0 (+ i 1)]
         [x x0 (step x)])
        [(= i N) results]
      (hash-table-update! results x (cut + 1 <>) 0))))

簡便のためにファイルに書き出す関数も定義しとこう。

(define (mcmc1 :key (samples 50000) (file "mcmc1.dat"))
  (with-output-to-file file
    (^[] (let1 result (run 0 samples)
           (do-ec (: x -10 20 0.5)
                  (let1 y (hash-table-get result x 0)
                    (format #t "~a\t~a\n" x (/. y samples))))))))

この関数を走らせるとデフォルトで結果をmcmc1.datに結果を書き出す。 (plot'result) とするとgnuplotを使って結果をプロットできる。

(mcmc1 :samples 1000) として1000サンプル取ってみたものがこれ。 緑のバーグラフがxの出現頻度。

[image]

なんとなくp(x)に近いけれどまだばらつきがある。

50000サンプルにするとこのとおり。

[image]

(なお、上に書いたとおり、MHで得られる分布は「p(x)に比例している」ってだけなので、 絶対値はあまり問題ではない。このグラフは一致度が見やすいように分布の方を スケールしてp(x)に合わせてある。)

1次元のp(x)ではあまりMCMCの有難味がわからないので、 次に2次元でやってみよう。→ Gauche:MetropolisProceduralModeling:MCMC-2D

More ...