Gauche:MetropolisProceduralModeling:MPM

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

ここまで来たらあとはTaltonさんの原論文の絵を見ればだいたいやってることはわかると思うんだ。って説明する気力がなくなってきた。

確率的文脈自由文法で列を生成する。例えば描画コマンド列みたいな。その時、最終結果だけでなく、ルートからすべての導出の木を保存しとく。
コマンドの中にはパラメータがある。枝の長さとか。
木のトポロジーが変わると、パラメータ数も変わっちゃう。つまり次元間のジャンプだ。
そこでRJMCMCを使おう。
- 次元内でのマルコフ過程は、パラメータのひとつをランダムに置き換える。
- 次元間でのジャンプは、導出木の非終端記号をひとつ選んで、そっから下の部分木を再導出する。

これで、「次元の選択も含めて、最適解を探索する」って目的は果たせる。

今回の実験では簡単に結果を見たいので、ごく簡単な2Dのブランチモデルを使ってみる。 Taltonさんの論文にあるやつに、ちょっとだけひねりを加えたもの (角度とか)。コードはmcmc5.scm。

生成文法は次のとおり。Fが枝でパラメータλが長さ。Lは(α+β)度だけ左に曲がる。 Rは(α-β)度だけ右に曲がる。Xは非終端記号。Zは枝の先っぽで、今回は花を描く。

(ところでTaltonさんの論文の"Grammar 1"のキャプションでは、Fのパラメータを N(l, .2)でサンプルするように書いてあるんだけど、ここはlに依存しないんじゃないのかなあ。 Fig.2の例とか見ると平均値lにはなっていないような?)

;; Grammar
;;  2-D simple branching model (a modified version of the example in
;;  Talton's paper)
;;
;;  V = {X()}
;;  T = {F(),L(),R(),Z,[,]}
;;  Σ = {l}
;;  ω = {X(1)}
;;  P =
;;    X(l) : l <= M ---(1-l/M)--> F(λ)[L(α,β)X(l+1)][R(α,β)X(l+1)]
;;    X(l) : l <= M ----(l/M)---> F(λ) Z
;;     where λ 〜 N(μ_λ, σ_λ), α 〜 N(μ_α, σ_α), β 〜 N(0, σ_β)

これでランダムに描画するとこんな木が描ける。センスが無いのは勘弁。

[image]

目的は、ターゲットとする形が与えられた時に、花の分布がその形に沿うようにすること。ターゲットは32x32のイメージIで受け取り、花を中心にしたガウス関数と積分する。

(define (likelihood tree I)             ;I :: <f32vector>
  (define buf (likelihood-buf tree))
  (define penalty 0)  ; if we're outside of canvas, we penalize it a lot.
  (define (roundI a min max)
    (round->exact (clamp (* (- I_size 1) (/. (- a min) (- max min)))
                         0 (- I_size 1))))
  (define (add-Z x y)
    (if (or (< x canvas-x0) (< canvas-x1 x)
            (< y canvas-y0) (< canvas-y1 y))
      (inc! penalty 5.0)
      (let ([ix (roundI x canvas-x0 canvas-x1)]
            [iy (roundI y canvas-y0 canvas-y1)])
        (do-ec (: y kernel-size)
               (: x kernel-size)
               (let ([jy (+ iy y (- kernel-center))]
                     [jx (+ ix x (- kernel-center))])
                 (when (and (< -1 jy I_size)
                            (< -1 jx I_size))
                   (f32vector-set! buf (+ jx (* jy I_size))
                                   (+ (f32vector-ref kernel
                                                     (+ x (* y kernel-size)))
                                      (f32vector-ref buf
                                                     (+ jx (* jy I_size))))))))
        )))
  
  (f32vector-fill! buf 0)
  (walk-node (Tree-ω tree) #f (^[x y level] (add-Z x y)))
  (f32vector-clamp! buf 0 1)
  (f32vector-sub! buf I)
  ;; Just to calculate likelihood, we can do (exp (- (f32vector-dot buf buf)))
  ;; but we want the squared image for visualization.
  (f32vector-mul! buf buf)
  (exp (- (+ penalty (f32vector-dot buf I0)))))

RJMCMCのステップの肝はここ。

(define (step tree I)
  (let1 u (random-real)
    (if (< u 0.95)
      (diffuse! tree I)
      (jump! tree I))
    (when (> (Tree-p tree) (Tree-max-p tree))
      (Tree-max-p-set! tree (Tree-p tree))
      (Tree-max-tree-set! tree (save-tree (Tree-ω tree))))))

95%の確率でパラメータをひとつだけ変える次元内の遷移であるdiffuse!を、 5%の確率でトポロジーを変えるjump!を試みる。なお、ステップの度に新しい木を作ってると重いので、各動きはacceptされるとtreeを破壊的に変更する。ステップ後に尤度を見て、より良いフィッティングになってたらその時の木を保存しとく。

diffuseは骨組みだけ見ればむしろmcmc4より簡単。resample!というやつが木のパラメータを破壊してしまうのだけれど、rejectされた時に元に戻すクロージャを返す。

(define (diffuse! tree I)
  (let* ([p (Tree-p tree)]
         [rollback! (resample! tree)]
         [p* (likelihood tree I)]
         [u (random-real)])
    (if (< u (min 1 (/ p* p)))
      (Tree-p-set! tree p*)             ;accepted
      (rollback!))))                    ;rejected

jumpの方は、だいたい論文通り。論文でτ'などと書いてあるのはτ*になってる。 Optimizationのところに書いてあるdelayed rejectは実装してない。

(define (jump! tree I)
  (receive (X w) (pick-variable tree)
    (let ([rollback! (tree-restorer tree X)]
          [p (Tree-p tree)]
          [qτ (/. w (Tree-Σw tree))]
          [Πφτ (subtree-probability tree X)])
      (derive! tree X)
      (let ([p* (likelihood tree I)]
            [qτ* (/. (find-weight tree X) (Tree-Σw tree))]
            [Πφτ* (subtree-probability tree X)])
        (let ([u (random-real)]
              [A (min 1 (/ (* qτ* p* Πφτ*) (* qτ p Πφτ)))])
          (if (< u A)
            (Tree-p-set! tree p*)       ;accepted
            (rollback!)))))))

コードを走らせるには、mcmc5.scmをロードしてrunにターゲットイメージを与える。今のところ一度runしちゃったらプログラム終了以外にターゲットを切り替える方法はない。いくつかサンプルのターゲットイメージがコードの下の方 (I-なんとか) に定義してある。

gosh> (run I-triangle2)
#<thread #f runnable 0x1a84400>

で窓が開いて試行が始まるはず。左下は尤度計算の様子。現在の木は *tree* に束縛されてるのでREPLから覗ける。

gosh> (describe *tree*)
#<Tree 0x246cd90> is an instance of class Tree
slots:
  ω         : (X 1 (Tree (F (0.5293490846909199)) (Tree (L (0.303594163981
  λs        : #((0.7263214279480569) (0.6567707105361102) (0.0856684682868
  αs        : #((0.5744914550716996) (0.5744914550716996) (0.9220321619166
  βs        : #((0.19125564015560367) (0.19125564015560367) (0.25929242561
  num-params: 367
  Xs        : ((X 5 (Tree (F (0.7263214279480569)) Z)) (X 5 (Tree (F (0.65
  ws        : (32 32 64 4 2 2 4 8 8 16 16 32 32 64 128 16 8 8 16 32 4 2 1 
  Σw        : 3062
  Lbuf      : #f32(0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  p         : 1.981363325244679e-43
  max-tree  : (X 1 (Tree (F (0.5293490846909199)) (Tree (L (0.303594163981
  max-p     : 2.695746077484929e-32

GLウィンドウ上でスペースで停止/動作切り替え、'm'で今までの最適形状表示、'i'でターゲットイメージ、尤度計算の表示など、rで木を新たに生成してフィッティングやりなおし。

案外収束が速いなって感じだが、2Dで自由度が少ないからってだけかもしれん。動作の様子: http://www.youtube.com/watch?v=W-Q6K4FZXJI

いくつか結果を貼っとく。