Gauche:MetropolisProceduralModeling:RJMCMC-fit

正月にMetropolis Procedural Modelingを基礎から理解する、のその4。

Gauche:MetropolisProceduralModeling:MCMC-fitでは次数を決めてからMCMCを走らせた。でも、あらかじめ次数がわかってないとしたら?

これまでのMCMCは、サンプルの存在する空間の次元は固定で、その空間の中でサンプルを動かしていた。それを拡張して、ある次元の空間から別の次元の空間へとサンプルがジャンプできるようにしたのが、Reversible Jump MCMCだ。

なんで"Reversible"がつくかっていうと、マルコフ過程でふらふらするサンプルが実際の確率分布に収束するためには、「状態xのもとで状態x*へ遷移する確率」と「状態x*のもとで状態xへ遷移する確率」が常に両方求まらないとならないから。詳しいことは参考文献を見られよ。ここでは論文に書いてある数式をコードに落とすことだけ考える。

コードは mcmc4.scm。尤度関数pはmcmc3と同じ。pickNとqも同じ。

違うのはstepから。確率δで、次数がひとつ上、もしくはひとつ下へとジャンプする。 (但し次数は1が最小、5が最大なので、1では下にジャンプしないかわりに上への移動確率を2δに、5では逆に下への移動確率を2δにしている。) walkがこれまでのstep関数のやってたことで、同じ次数m内でパラメータを動かす。

;; Single MH step.  Returns next m, θ, and error^2 value.
(define (step D m θ)
  (let1 u0 (random-real)
    (case m
      [(1) (if (< u0 (* 2 δ))
             (jump↑ D m (+ m 1) θ)
             (walk D m θ))]
      [(5) (if (< u0 (* 2 δ))
             (jump↓ D m (- m 1) θ)
             (walk D m θ))]
      [else (cond [(< u0 δ)       (jump↓ D m (- m 1) θ)]
                  [(< u0 (* 2 δ)) (jump↑ D m (+ m 1) θ)]
                  [else (walk D m θ)])])))

(define (walk D m θ)
  (let* ([u  (random-real)]
         [θ* (pickN θ)]
         [A  (min 1 (/ (* (p D θ*) (q θ* θ)) (* (p D θ) (q θ θ*))))]
         [θ. (if (< u A) θ* θ)]
         [e  (err D (polynomial θ.))])
    (values m θ. e)))

そしてジャンプする遷移。もうちょい簡略化できるかもしれんが、とりあえず理解したところに忠実にコード書いた。

;; jump from n to m, n > m
(define (jump↓ D n m θ)
  (let* ([u  (random-real)]
         [θ* (pickN (map (cut * 5 <>) θ))]
         [q_n->m (/ (φ (/ (* 5 (car θ)) σ_u)) σ_u)]
         [A  (min 1 (* (/ (p D θ*) (p D θ))
                       (if (= n 5) 2 1) ; q(n|m)/q(m|n)
                       (/ q_n->m)
                       (expt 5 (+ n 1))))]) ;Jacobian
    (if (< u A)
      ;; Jump is accepted.  We discard θ*'s first element.
      (values m (cdr θ*) (err D (polynomial (cdr θ*))))
      ;; Jump is rejected.
      (values n θ (err D (polynomial θ))))))

;; jump from n to m, n < m
(define (jump↑ D n m θ)
  (let* ([u  (random-real)]
         [θadj (map (^t (/ t 5)) (cons (cpick 0 σ_u) θ))]
         [θ* (pickN θadj)]
         [q_m->n (/ (φ (/ (car θadj) σ_u)) σ_u)]
         [A  (min 1 (* (/ (p D θ*) (p D θadj))
                       (if (= n 1) 2 1) ; q(n|m)/q(m|n)
                       q_m->n
                       (/ (expt 5 (+ m 1)))))])   ;Jacobian
    (if (< u A)
      ;; Jump is accepted.
      (values m θ* (err D (polynomial θ*)))
      ;; Jump is rejected.
      (values n θ (err D (polynomial θ))))))

ドライバ。mcmc3から違うところは、次数mも一緒に引き回してるところ。

(define (run D m0 θ0 N)
  (receive (m θ e) (step D m0 θ0)
    (let loop ([i 0] [m m] [θ θ] [m-best m] [θ-best θ] [e-best e])
      (if (= i N)
        (values m-best θ-best e-best)
        (receive (m1 θ1 e1) (step D m θ)
          (print i "\t" e1 "\t" m1)
          (if (< e1 e-best)
            (loop (+ i 1) m1 θ1 m1 θ1 e1)
            (loop (+ i 1) m1 θ1 m-best θ-best e-best)))))))

(define (mcmc4 θ0 :key (samples 50000)
                       (datafile "mcmc3.input")
                       (errfile "/dev/null"))
  (let1 D (map (cut apply cons <>) (slices (file->sexp-list datafile) 2))
    (with-output-to-file errfile
      (^[] (run D (- (length θ0) 1) θ0 samples)))))

さて走らせてみる。

gosh> (mcmc4 '(0 0) :samples 200000 :errfile "mcmc4.err")
3
(0.051415295647450496 -0.15600867185158387 0.15430664534317431 -0.009431577640634155)
28.055289211817964
gosh> (plot'result '(0.051415295647450496 -0.15600867185158387 0.15430664534317431 -0.009431577640634155))
#<undef>
gosh> (plot'error :xmin 0 :ymax 500)
#<undef>

[image]

下の方のグラフは、赤が誤差、緑が次数。ジャンプ直後にぐわっと誤差が大きくなることが多いんだけど、その後しばらく同じ次数内で探索して低いところ探して、またジャンプして、って感じ。

このデータはわりと綺麗に探してるけれど、場合によっては5次まで探しに行ったところで誤差がどかんと跳ね上がって戻ってこないことがあるんで、この実装はあんまり良くないのかもしれん。(間違えてなければ、理論上はいつかは収束するはずだけど。)

さて、RJMCMCの原理までわかったところで、いよいよTalton論文の実装をやってみよう。 → Gauche:MetropolisProceduralModeling:MPM