HANSEI in F#

HANSEI.fs

//HANSEIとはprobabilistic programmingを行うためのtaglessなインタプリタのこと。
//確率モデルをプログラムの形で書くとパラメータの推論を自動で行なってくれる。
//oleg氏の論文だとreify/reflectを使ったdirect styleのモナドで記述されている。
//http://okmij.org/ftp/kakuritu

//要はこちらのスライドで紹介されているものをmonadやHOASを使って内部DSLとして実装したもの。
//http://www.slideshare.net/tsubosaka/infernetlda

//こちらの方の記事も参考になる。
//http://d.hatena.ne.jp/rst76/20100706/1278430517

//importance samplingについてはこちらの記事が参考になる。
//http://d.hatena.ne.jp/teramonagi/20120107/1325927802

//パラメータ推論の戦略は色々ある。
//http://okmij.org/ftp/kakuritu/inference.ml


//continuation monad
type Cont<'R, 'A> = (('A -> 'R) -> 'R)
type ContBuilder() =
  member this.ReturnFrom(x) = x
  member this.Return(x) = fun k -> k x 
  member this.Bind(m, f) = fun k -> m (fun a -> (f a) k)
let cont = new ContBuilder()


//probability monad
//基本的には以下の確率モナドと似たものになっている。
//http://www.randomhacks.net/articles/2007/02/21/refactoring-probability-distributions
type Prob = float
type VC<'T> = V of 'T | C of Lazy<PV<'T>>
and PV<'T> = list<Prob * VC<'T>>
type PM<'T> = Cont<PV<bool>, 'T>
type Arr<'A, 'B> = 'A -> PM<'B>
type ProbBuilder() =
  member this.ReturnFrom(x) = x
  member this.Return(x) = [(1.0, V(x))]
  member this.Bind(m, f) =
    let g = function
      | (p, V x) -> (p, C(lazy(f x)))
      | (p, C(Lazy(t))) -> (p, C(lazy(this.Bind(t, f))))
    List.map g m
let prob = new ProbBuilder()


//tagless interpreter
let b = prob.Return
let dist ch = fun k -> List.map (fun (p,v) -> (p, C(lazy(k v)))) ch
let flip_ p = dist [(p, true); (1.0-p, false)]
let neg m = cont { let! x = m in return not x }
let con m1 m2 = cont { let! x = m1 in let! y = m2 in return x && y }
let dis m1 m2 = cont { let! x = m1 in let! y = m2 in return x || y }
let if_ et e1 e2 = cont { let! t = et in return! if t then e1 else e2 }
let lam e = cont.Return(e << cont.Return)
let app e1 e2 = cont { let! f = e1 in let! x = e2 in return! f x }
let let_ e f = app (lam f) e
let reify0 m = m prob.Return
let exact_reify model = explore None (reify0 model)
let fail () = dist []

//構文木の例。
let grassModel1 () =
  let_ (flip_ 0.3) (fun rain ->
  let_ (flip_ 0.5) (fun sprinkler ->
  let_  (dis (con (flip_ 0.9) rain)
             (dis (con (flip_ 0.8) sprinkler)
                  (flip_ 0.1))) (fun grassIsWet ->
  if_ grassIsWet rain (fail ()))))

//実行例。
exact_reify (grassModel1())

//continuation monadとしても使える。
let grassModel2 () = cont {
  let! rain           = flip_ 0.3
  let! sprinkler      = flip_ 0.5
  let! wetInRain      = flip_ 0.9
  let! wetInSprinkler = flip_ 0.8
  let! wetInOther     = flip_ 0.1
  let grassIsWet = rain && wetInRain || sprinkler && wetInSprinkler || wetInOther
  if grassIsWet then
    return rain
  else
    return! fail()
}

exact_reify (grassModel2())

//monadとhigher-order abstract syntaxは似ている。
//両方ともオブジェクトレベルの変数を表現するのにメタレベルの仕組みを流用している。


let explore (maxdepth : int option) (choices : PV<'T>) =
  let insertWith fn k v m =
    let v' = Map.tryPick (fun k' v' -> if k = k' then Some(v') else None) m
    match v' with
    | Some(v') -> Map.add k (fn v v') m
    | None -> Map.add k v m
  let rec loop p depth down susp answers =
    match (down, susp, answers) with
    | (_, [], answers) -> answers
    | (_, (pt, V v) :: rest, (ans, susp)) ->
      loop p depth down rest (insertWith (+) v (pt*p) ans, susp)
    | (true, (pt,C(Lazy(t)))::rest, answers) ->
      let down' = match maxdepth with Some x -> depth < x | None -> true
      loop p depth true rest <| loop (pt*p) (depth+1) down' t answers
    | (down, (pt,c)::rest, (ans,susp)) ->
      loop p depth down rest (ans, (pt*p,c)::susp)
  let (ans, susp) = loop 1.0 0 true choices (Map.empty, [])
  Map.fold (fun a v p -> (p, V v)::a) susp ans

let normalize (choices : PV<'T>) =
  let denominator =
    List.map fst choices
    |> List.sum
  List.map (fun (p, v) -> (p/denominator, v)) choices


//もっと簡単な例。
let coin_flip () = dist [(0.5, true); (0.5, false)]

let two_coin_toss () = cont {
  let! x = coin_flip ()
  let! y = coin_flip ()
  if x || y then
    return (x, y)
  else
    return! fail ()
}

two_coin_toss ()
|> exact_reify


//モンティ・ホール問題。
let uniform xs =
  let p = 1.0 / float(List.length xs)
  xs
  |> List.map (fun x -> (p, x))
  |> dist

type Door = A | B | C
let doors = [A; B; C]

let monty_hall1 () = cont {
  let! answer = uniform doors
  let! selected = uniform doors
  return (answer, selected)
}

let monty_hall2 () = cont {
  let! answer = uniform doors
  let doors' = List.filter ((<>) answer) doors
  let! chosen = uniform doors'
  let doors'' = List.filter ((<>) chosen) doors
  let! selected = uniform doors''
  return (answer, selected)
}

//1より2の方が確率が高くなっている。
exact_reify (monty_hall1())
exact_reify (monty_hall2())


//http://www.slideshare.net/gorn708/infernetlda-12845396の「簡単な確率モデルを作ってみる」の例をやってみる。
let bernoulli p = dist [(p, true); (1.0 - p, false)] //true/falseじゃないほうがいいかも。

let sample_model () = cont {
  let! firstCoin = bernoulli 0.5
  let! secondCoin = bernoulli 0.5
  let bothHeads = firstCoin && secondCoin //どちらかあるいは両方のコインが、
  if bothHeads = false then //裏になったときの、
    return firstCoin //最初のコインの分布を求めてみる。
    //疑問点:Infer.NETのinfer関数呼び出しが複数あるプログラムをどう翻訳する？ モデル自体を複数回コピペするしかないのでは…。
  else
    return! fail ()
}

let dist = exact_reify (sample_model ())
normalize dist
|> printfn "%A" //表になる確率は0.3333333...だった。合ってた。

memo.txt

PRML memo

chapter 11.
変分ベイズ法やEP法も推論アルゴリズムだが、これらは決定論的近似。11章でやるのは確率的近似。
数値的サンプリング（モンテカルロサンプリング）の目的：
f(z)の確率p(z)の下での期待値 E[f] = ∫f(z)p(z)dz を求めたいが、厳密にやるのは無理なので近似値を出したい。
#モンテカルロ法は解が確率的でラスベガス法は計算リソースが確率的。
これはp(z)から独立に摘出されたサンプル集合z_lがあれば f' = (1/L)* Σf(z_l) のように計算できる。
しかしこれだとp(z)が小さくf(z)が大きいような部分があると精度が悪くなるのでサンプル数を多くしないと駄目。

importance sampling
E[f] = ∫f(z)*(p(z)/q(z))*q(z)dz ≒ (1/L)*Σ(p(z_l)/q(z_l))*f(z_l)
q(z)がproposal distribution。r_l = p(z_l)/q(z_l) がimportance weight。
#これって結局「proposal distributionはだいたいこういう分布でいいだろ」みたいな勘が必要ということ？

Box-Muller法のような変換法は逆関数を求める必要があるが、現実には逆関数を求められない分布が多い。
なのでrejection samplingやimportance samplingが必要。
しかしこれらは低次元でのsamplingにしか使えない。高次元ではMCMCのような手法を使う。
logic samplingはimportance samplingの特殊な場合。


http://www.cs.rutgers.edu/~ccshan/rational/ml-slides.pdf
Depth-first enumeration = exact inference
Random dive = rejection sampling
Dive with look-ahead = importance sampling
微妙に違う意味で言葉を使っている…？


Noisy-ORについて。参考になる。
http://d.hatena.ne.jp/yosshi71jp/20120610/1339327581


http://dickie.ees.kyushu-u.ac.jp/HP/wadalab/wadalab2005/graduate/tagawa.pdf
変分ベイズ法は、事後分布近似を計算する必要があるため、EMアルゴリズムよりも計算量は多くなってしまうが、
事前知識を導入し学習データ数が少なくても過学習を起こしにくいという利点があり、十分に学習ができない条件下でも、高い推定精度を保つことができる。
変分ベイズ法はEMアルゴリズムのベイズ推定への拡張とみなすことができる。

reifyとreflectを使ってるとメモ化できるらしい。CPSでもできるかも？
exact inferenceはprobability monadのtreeを全部traverseする。しかし現実のsearch treeはtoo large to traverse completely。
なので一部だけtraverseするapproximate inferenceが必要。
どのnodeをtraverseするのかについて色々な戦略がある。
重みに従ってランダムに選ぶ→このrejection samplingはpoorly approximate。分布のtailとかとれないし、長すぎるパスを取ったり？して実用的でない。
Importance samplingっていうrejection samplingを改良したやつが良いらしい。
IBALというprobabilistic programming languageはevidence pushingという手法を使っている。HANSEIは内部DSLなのでソースコードに手を入れられないからこの技法は使えない。

Importance samplingのようなapproximate inferenceならsearch treeの一部しか観測しないので、lazyにしておけば高速化できるはず。

EP法
http://www.ohno.mgmt.waseda.ac.jp/Thesis/06yang.pdf
変分ベイズ
http://www.kecl.ntt.co.jp/as/members/ueda/pdf/2001_sairyo.pdf

疑問点：PHOASを使えばtermの分解ができるはずだから、それを使ってevidence pushingが適用できたりしないか？　PHOASのためにhigher-kinded polymorphismが必要？

Infer.NETで使われている推論アルゴリズム。cf.http://www.slideshare.net/tsubosaka/tokyowebmininginfernet
Expectation Propagation for Approximate Bayesian Inference
http://research.microsoft.com/en-us/um/people/minka/papers/ep/minka-ep-uai.pdf
Variational Message Passing
http://johnwinn.org/Publications/papers/VMP2004.pdf
Gibbs Sampling
http://www-users.cs.york.ac.uk/jc/teaching/agm/lectures/lect13/lect13.pdf

Infer.NETの競合プロジェクト。
BUGS
http://www.mrc-bsu.cam.ac.uk/bugs/
Factorie
http://factorie.cs.umass.edu/

こんなものもあるらしい。
Hierarchical Bayes Compiler
http://www.cs.utah.edu/~hal/HBC/

多分importance samplingとrejection samplingはあまり関係ない（改良したものではない。それぞれ長所短所がありそう）。

1.連続確率分布が扱え 2.custom確率モデルと 3.確率密度関数が自然に表現できる ような確率プログラミング言語はまだ無い（無かった）らしい。
でもそれを作ることに挑戦してみたよという論文がこれ。今度じっくり読む。
http://ashishagarwal.org/wp-content/uploads/2011/10/POPL2012-pdf-type-theory-preprint.pdf

PRISM
http://sato-www.cs.titech.ac.jp/prism/
PRISMはBayesian networkやPCFGなどをモデリングするための言語。explanation graph上で確率計算を行う。Prologのような論理型言語に近いらしい。

HANSEI as a Declarative Logic Programming Language
http://okmij.org/ftp/kakuritu/logic-programming.html
HANSEIを宣言的な論理型言語として解釈する話（だと思う）。HANSEIもバックトラックするしPRISMと似ているのかも。

A General MCMC Method for Bayesian Inference in Logic-Based Probabilistic Modeling
http://ijcai.org/papers11/Papers/IJCAI11-248.pdf
PRISMでMCMCを行う手法。PCFGでのMetropolis Hastings samplingを一般化してexplanation graph上でサンプリングを行う。