すごいH本の素朴な確率モナド

年末年始の連休から中五日あっての三連休で、正月ボケをぶり返してしまいそうな今日この頃ですが、いかがお過ごしでしょうか。

すごいH本こと、書籍「すごいHaskellたのしく学ぼう！」の最後のほう、第１４章「もうちょっとだけモナド」の 14.8 (P356)にて、素朴な確率モナドが紹介されています。

すごいHaskellたのしく学ぼう!

普通、モナドは作りたいと思って作るものではありません。むしろ、とある問題のある側面をモデル化した型を作り、後からその型が文脈付きの値を表現していてモナドのように振る舞うと分かった場合に、Monad のインスタンスを与える場合が多いです。

というのが印象的で。ふむふむ確かになるほどなあという感じです。

ぼけーっと、ただ連休をだらだらと過ごすだけなのもなんなので、正月ボケのリハビリを兼ねて何か書いておこうかなということで、これを F# で書いてみようと思います。

確率を表現するための有理数を表す型

数学では通常、確率はパーセント（％）ではなく、0 から 1 までの実数で表します。確率が 0 ということは絶対にありえないということであり、確率が 1 というのは確実に起こるということを意味します。確率を浮動小数点で表すのも間違いではないのですが、どうしても精度が落ちてしまう。そこで Haskell では、Rationalという分数を表すために最適な有理数を表す型があり、例えば ４分の１は、1%4 のように、分子と分母は % で区切って表現することができる。

では、F# はどうでしょう。標準には用意されていませんが、F# では、F# PowerPack という追加ライブラリにて数学に関する様々な機能が提供されています。これを導入することで分数の表現に対応することができます(NuGetで簡単に導入することもできます)。有理数を表すことができる BigRational という型が定義されているので、それを使えます。BigRational は、Nリテラルを用いて表現することができ、４分の１は、1N/4N というように表せます。

F# で素朴な確率モナド

Haskellでの実装例は書籍や(Learn You a Haskell for Great Good! - For a Few Monads More)に出ている通りなので、そちらを参照されたい。

BigRational型と FSharpx を使って、F# で素朴な確率モナドをとりえず実装してみる。

namespace FSharpx.Monad

// 素朴な確率モナド
module Probability =
  let probMap f m = List.map (fun (x,p) -> (f x, p)) m

  type ProbBuilder() =
    member this.ReturnFrom(x) = x
    member this.Return(x) = [x,1N/1N]
    member this.Bind(m, f) = 
      let flatten xs = 
        let concatMap f m = List.concat( List.map (fun x -> f x) m )
        let multAll (innerxs,p) = List.map (fun (x,r) -> (x, p*r)) innerxs
        concatMap multAll xs
      flatten (probMap f m) 
        
    member this.Zero () = []

  let prob = new ProbBuilder()

  open FSharpx
  open Operators 
  let inline returnM x = returnM prob x 
  let inline (>>=) m f = bindM prob m f
  let inline (=<<) f m = bindM prob m f
  let inline (<*>) f m = applyM prob prob f m
  let inline ap m f = f <*> m
  let inline map f m = liftM prob f m
  let inline (<!>) f m = map f m
  let inline lift2 f a b = returnM f <*> a <*> b
  let inline lift3 f a b c = returnM f <*> a <*> b <*> c
  let inline ( *>) x y = lift2 (fun _ z -> z) x y
  let inline ( <*) x y = lift2 (fun z _ -> z) x y
  let inline (>>.) m f = bindM prob m (fun _ -> f)
  let inline (>=>) f g = fun x -> f x >>= g
  let inline (<=<) x = flip (>=>) x

使ってみる。３枚のコイン（イカサマコインが１つ混入している）がすべて裏が出る確率を出す。

module Program =
  open System
  open FSharpx.Monad.Probability

  type Coin = Heads | Tails 
  
  let coin = [(Heads,1N/2N); (Tails,1N/2N)]
  let loadedCoin = [(Heads,1N/10N); (Tails,9N/10N)]

  let flipThree = prob {
    let! a = coin
    let! b = coin
    let! c = loadedCoin
    return List.forall (function |Tails->true |_->false) [a;b;c]
  }

  flipThree |> printfn "%A"

実行結果

[(false, 1/40N); (false, 9/40N); (false, 1/40N); (false, 9/40N); (false, 1/40N); (false, 9/40N); (false, 1/40N); (true, 9/40N)]

確率モナドによって、３枚とも裏が出る確率は、４０分の９であると導きだすことができた。すごいH本と同じ結果になりましたね。めでたしめでたし。

続いて、6面のサイコロを２回振ったとき、その出目の合計値ごとの確率を出してみる。

  let d sides = [for i in [1 .. sides] -> (i, 1N/ BigRational.FromInt(sides))]
  let dice = d 6

  let diceTwoSum = prob {
    let! a = dice
    let! b = dice
    return a+b
  }
  diceTwoSum |> printfn "%A"

実行結果

[(2, 1/36N); (3, 1/36N); (4, 1/36N); (5, 1/36N); (6, 1/36N); (7, 1/36N); (3, 1/36N); (4, 1/36N); (5, 1/36N); (6, 1/36N); (7, 1/36N); (8, 1/36N); (4, 1/36N); (5, 1/36N); (6, 1/36N); (7, 1/36N); (8, 1/36N); (9, 1/36N); (5, 1/36N); (6, 1/36N); (7, 1/36N); (8, 1/36N); (9, 1/36N); (10, 1/36N); (6, 1/36N); (7, 1/36N); (8, 1/36N); (9, 1/36N); (10, 1/36N); (11, 1/36N); (7, 1/36N); (8, 1/36N); (9, 1/36N); (10, 1/36N); (11, 1/36N); (12, 1/36N)]

ここまでがすごいH本で書かれている範囲でできること。これから先については、読者への演習問題としている。

上記の実行結果を見てわかるように、確率の結果がまとまっておらず、バラバラに出力されていて結果の内容がわかりにくい。これはひどい。

できれば、

[(false, 31/40N); (true, 9/40N)]

とか

[(2, 1/36N); (3, 1/18N); (4, 1/12N); (5, 1/9N); (6, 5/36N); (7, 1/6N); (8, 5/36N); (9, 1/9N); (10, 1/12N); (11, 1/18N); (12, 1/36N)]

というように、結果が一致する事象の確率については１つにまとめて出力してくれるのが分かり易くて理想だよね、と。せっかくなので、この演習問題をやってみましょう。

とりあえず、結果が一致する事象の確率を１つにまとめてみる

あまり何も考えずに、とりあえず実装してみた版。

  let rec merge (k,p) xs = xs |> function
    | []  -> []
    | (k,p)::kps -> kps |> function
      | [] -> [(k,p)]
      | (k',p')::kps' ->
        if k = k' then (k,p+p')::(merge (k,p) kps')
        else (k,p)::(merge (k',p') kps)

  let agglomerate f pd = 
    let xs : ('b * BigRational) list = (probMap f pd) |> List.sort 
    List.foldBack merge pd xs

  let agg pd = agglomerate id pd 

使ってみる。

  let flipThree = prob {
    let! a = coin
    let! b = coin
    let! c = loadedCoin
    return List.forall (function |Tails->true |_->false) [a;b;c]
  }

  flipThree |> agg |> printfn "%A"

  //let flipThree2 = agg <| lift3 (fun a b c -> List.forall (function |Tails->true |_->false) [a;b;c]) coin coin loadedCoin
  //flipThree2 |> printfn "%A"

実行結果

[(false, 31/40N); (true, 9/40N)]

使ってみる。

  let d sides = [for i in [1 .. sides] -> (i, 1N/ BigRational.FromInt(sides))]
  let dice = d 6

  let diceTwoSum = prob {
    let! a = dice
    let! b = dice
    return a+b
  }

  diceTwoSum |> agg |> printfn "%A"

  //let diceTwoSum2 = agg <| lift2 (+) dice dice
  //diceTwoSum2 |> printfn "%A"

実行結果

[(2, 1/36N); (3, 1/18N); (4, 1/12N); (5, 1/9N); (6, 5/36N); (7, 1/6N); (8, 5/36N); (9, 1/9N); (10, 1/12N); (11, 1/18N); (12, 1/36N)]

うん。とりあえず動いているね。これで一応目的は達成できているのだけど、なんだか冗長な感じがするしカッコ悪い。俺が欲しいのコレジャナイ感がぱない。もっとシンプルに行きたい。

結果が一致する事象の確率を１つにまとめる（改訂版）

どのあたりがコレジャナイ感を出しているのか。落ち着いて先ほどの実装をよく見てみてみよう。

  let rec merge (k,p) xs = xs |> function
    | []  -> []
    | (k,p)::kps -> kps |> function
      | [] -> [(k,p)]
      | (k',p')::kps' ->
        if k = k' then (k,p+p')::(merge (k,p) kps')
        else (k,p)::(merge (k',p') kps)

  let agglomerate f pd = 
    let xs : ('b * BigRational) list = (probMap f pd) |> List.sort 
    List.foldBack merge pd xs

  let agg pd = agglomerate id pd 

List.sort でソートしたリストと、ソートする前のリストとを比較して、再帰でマージしながら結果をまとめあげる実装となっている。

そもそもここでやりたいことは、集合として確率の結果をまとめ上げたいということ。集合を扱いたい場合、F# では set が使える。また、コレジャナイ実装では、List.foldBack で入力要素を順々に受け取りながら marge 関数で確率の和を求めながら結果の状態を順次更新していっているが、set を使って集合化することができれば、集合の要素ごとの確率の和をそれぞれ算出してゆくだけでよいことになる。あ、それって List.reduce 使えばいんじゃね？ となる。

ところで、List.reduce とはなんだったのか。例えば、List.foldを用いてintのリストの和を求める場合を思い出してみよう。

List.fold (+) 0 [1;2;3] 

のように書けるのでした。育てる種となる初期値の 0 を与えて、次々にリストを畳み込むことにより、結果 6 が得られる。

ここで、育てる種となる初期値の 0 を与えずにリストの和を求めるには、育てる種の初期値としてリストの最初の要素を採用すればよい。最初の要素と次の要素によって演算を開始するという処理を行えばよいことなる。

こう書く事ができる。

List.reduce (+) [1;2;3]

そう、List.fold が簡易化されたものが List.reduce ということだった。

ということで集合を扱える set と 育てる種を必要としない List.reduce を用いて実装すると次のように書ける。

  let merge f (a,x) (b,y) : 'a * BigRational = f a b, x + y

  let agglomerate f pd =
    let d = pd |> List.map(fun (x,_) -> x) |> set |> Set.toList 
    List.map (fun x -> List.filter(fun (y,b) -> x = y) pd |> List.reduce (merge f)) d

  let agg pd = agglomerate (fun _ x -> x) pd

ほむ。だいぶシンプルになりました。これはリハビリをして正解でしたね。