Discrete Probability Monad with F# computation Expressions that runs in Fable.

#Probabilistic Programming in the Browser.md

Discrete Probability Monad with F# computation Expressions that runs in Fable

Copy and paste each part in turn into http://fable.io/repl/. Probabilistic programming part should start around line 220. Should be compatible with almost everything in http://greenteapress.com/wp/think-bayes/

Also provides text histogram plots!

Common.fs

open Fable.Core.JsInterop
open Fable.Import 
open Fable.Core
open System

let flip f x y = f y x

let konst x _ = x

let keepLeft f (x,y) = x , f y

let keepRight f (x,y) = f x , y

let round (n:int) (x:float) = System.Math.Round(x,n)

let stringr n x = string (round n x)

let stringr100 n x = string (round n (x * 100.))
///////////////////

module Map  =
       let sum m = Map.fold (fun sum _ x -> sum + x) 0. m

module Array =
   let inline normalizeWeights (a:_[]) =
       let sum = Array.sumBy snd a 
       Array.map (fun (x,w) -> x, w/sum) a

   let indexed a = Array.mapi (fun i x -> (i,x)) a
   
   let selectColumn c a = 
       Array.map (indexed >> Array.filter (fst >> (=) c) >> Array.map snd) a  

module String = 
   let pad padlen (s:string) = s + String.replicate (max 0 (padlen - s.Length)) " "
   let inline toupper (str:string) = str.ToUpper()
   

let inline joinToStringWith sep (s:'a seq) = String.Join(sep, s)           
////////////
let buildTableRow (collens:_[]) (row:string[]) =
       row |> Array.mapi (fun i s ->  String.pad collens.[i] s) |> joinToStringWith " | "


let makeTable newline headers title (table:string[][]) =
       let hlen = Array.map String.length headers

       let lens = table |> Array.map (Array.map (String.length))

       let longest = [|for c in 0..headers.Length - 1 -> max hlen.[c] (Array.selectColumn c lens |> Array.map Seq.head |> Array.max)|]

       let t0 = table |> Array.map (buildTableRow longest) |> joinToStringWith newline

       let hrow = [|headers; [|for i in 0..headers.Length - 1 -> String.replicate longest.[i] "-"|]|] 
                  |> Array.map (buildTableRow longest) 
                  |> joinToStringWith newline
       String.Format("{0}{1}{1}{2}{1}{3}", (String.toupper title), newline, hrow, t0)

DiscreteProbabilityMonad.fs

////////////////////////

let rnd = System.Random()                

//from Expert F# code by Don Syme
//added some helpers to make the code more easily express bayesian formalisms

type Distribution<'T when 'T : comparison> =
    abstract Sample : 'T
    abstract Support : Set<'T>
    abstract Expectation: ('T -> float) -> float
                                                          
let always x =
    { new Distribution<'T> with
    member d.Sample = x
    member d.Support = Set.singleton x
    member d.Expectation(H) = H(x) }       


let coinFlip (p:float) (d1:Distribution<'T>) (d2:Distribution<'T>) =
    if p < 0.0 || p > 1.0 then failwith "invalid probability in coinFlip"
    { new Distribution<'T> with
        member d.Sample =         
            if rnd.NextDouble() < p then d1.Sample else d2.Sample
        member d.Support = Set.union d1.Support d2.Support
        member d.Expectation(H) =     
          p * d1.Expectation(H) + (1.0-p) * d2.Expectation(H)}


let bind (dist:Distribution<'T>) (k: 'T -> Distribution<'U>) =
    { new Distribution<'U> with
        member d.Sample = (k dist.Sample).Sample
        member d.Support = Set.unionMany (dist.Support |> Set.map (fun d -> (k d).Support))
        member d.Expectation H = dist.Expectation(fun x -> (k x).Expectation H) }

/////////////////////////////////////

type DistributionBuilder() =
    member x.Delay f = bind (always ()) f
    member x.Bind (d, f) = bind d f
    member x.Return v = always v
    member x.ReturnFrom vs = vs
 
let distr = DistributionBuilder()  

///////////////////

let weightedCases (inp: ('T * float) list) =
    let rec coinFlips w l =
        match l with
            | [] -> failwith "no coinFlips"
            | [(d,_)] -> always d
            | (d,p)::rest -> coinFlip (p/(1.0-w)) (always d) (coinFlips (w+p) rest)
    coinFlips 0.0 inp

let countedCases inp =
    let total = Seq.sumBy (fun (_,v) -> v) inp
    weightedCases (inp |> List.map (fun (x,v) -> (x, float v / float total)))

/////////Helpers
let dualCase ((_,p) as o1) o2 = weightedCases [o1; o2 , 1. - p]

let bernoulli p = dualCase (true, p) false

let uniform (items:'a list) = 
    let num = float items.Length
    List.map (fun item -> item, 1./num) items |> weightedCases

let categorical distr = weightedCases distr    

////////////////////////////// 

let conditionalProb eventx condition  (t:Distribution<'a>) =
   let n = t.Expectation(fun c -> if eventx c && condition c then 1. else  0.)  

   let z = t.Expectation(fun c -> if condition c then 1. else  0.) 
   n / z


let probabilityOf eventx (t:Distribution<'a>) =
    t.Expectation(fun c -> if eventx c then 1. else  0.) 

let rec update likelihood state obsv = 
        Map.map (fun item (prob:float) -> prob * likelihood obsv item) state    

let updates likelihood state obsvs =
    List.fold (update likelihood) state obsvs

let probabilityOfWeighted eventx (weights:Map<_,_>)  =
    let n = Map.filter (fun k _ -> eventx k) weights |> Map.sum   
    let z = Map.sum weights  
    n / z

let probabilityOfWeighted2 eventx (weights:('a * float)[])  =
    let n = Array.filter (fun (k, _) -> eventx k) weights |> Array.sumBy snd
    let z = Array.sumBy snd weights  
    n / z
    
let pmfWeighted (weights:Map<_,_>)  =    
    let t = Map.sum weights  
    Map.map (fun _ p -> p/t) weights 

let pmf (dist:Distribution<_>) = 
    dist.Support 
    |> Set.toArray 
    |> Array.map (fun c -> c, probabilityOf ((=) c) dist)  


let histogram len (d:_[]) =
    let _, maxp = Array.maxBy snd d
    Array.map (fun (x,p:float)-> 
          [|sprintf "%A" x ;
            stringr100 2 p + "%"; 
            String.replicate (int(round 0 (p/maxp * len))) "#" |]) d
    |> makeTable "\n" [|"item";"p"; ""|] ""  

/////////////
    
let toBits x = x / log 2. 

let inline log0 x = if x = 0. then 0. else log x

let inline entropy dist = -Seq.sumBy (fun (_,p) -> p * log0 p) dist

let entropyDistr d = entropy (pmf d)

let inline conditionWith projectWith f d = Array.filter (fst >> projectWith >> f) d |> Array.normalizeWeights

let inline conditionEntropyOn projectWith x d = conditionWith projectWith ((=) x) d |> entropy


let conditionalEntropy projectWith (joint:Distribution<_>) =
    Set.map projectWith joint.Support 
    |> Seq.sumBy (fun x -> 
          let p = probabilityOf (projectWith >> (=) x) joint
          let h = conditionEntropyOn projectWith x (pmf joint)
          p * h)

let mutualInformation (joint:Distribution<_>) =
    joint.Support |> Seq.sumBy (fun(x,y) ->
        let px = probabilityOf (fst >> (=) x) joint
        let py = probabilityOf (snd >> (=) y) joint 
        let pxy = probabilityOf ((=) (x,y)) joint
        pxy * log0(pxy/(px * py)))


let kldivergence (pA:Distribution<_>) (pB:Distribution<_>) =
    pA.Support |> Seq.sumBy (fun x ->
        let p_a = probabilityOf ((=) x) pA
        let p_b = probabilityOf ((=) x) pB
        p_a * log0(p_a/ p_b))

let kldivergence2 pA pB =
    pA |> Array.sumBy (fun ((x,_)) ->
        let p_a = probabilityOfWeighted2 ((=) x) pA
        let p_b = probabilityOfWeighted2 ((=) x) pB
        p_a * log0(p_a/ p_b))
//////////////////////////

Sample.fs

let coin = distr { let! n = uniform [0.01..0.05..1.0]
                   return n}

printfn "Uniform Prior"
printfn "p >= 50%% = %s%%" (stringr100 0 (probabilityOf ((<=) 0.5) coin))

printfn "%s" (histogram 20. (pmf coin))

let fair = [true; false; false; true; true; false; true; false]

let biased = [ false; false; false; false; false; false; true; false]

let fairPosterior = fair |> updates (fun c p -> if c then p else 1. - p) (Map.ofArray (pmf coin))

printfn "\nFair tosses Posterior:"
printfn "p >= 50%% = %s%%" (stringr100 1 (probabilityOfWeighted ((<=) 0.5) fairPosterior))
kldivergence2 (Map.toArray fairPosterior) (pmf coin) |> toBits |> round 2 |> printfn "Surprise: %g bits"
printfn "%s" (histogram 30. (fairPosterior |> pmfWeighted |> Map.toArray))

printfn "\nBiased tosses Posterior:"
let biasedPosterior = biased |> updates (fun c p -> if c then p else 1. - p) (Map.ofArray (pmf coin))

printfn "p >= 50%% = %s%%" (stringr100 1 (probabilityOfWeighted ((<=) 0.5) biasedPosterior))
kldivergence2 (Map.toArray biasedPosterior) (pmf coin) |> toBits |> round 2 |> printfn "Surprise: %g bits"
printfn "%s" (histogram 30. (biasedPosterior |> pmfWeighted |> Map.toArray))