Instantly share code, notes, and snippets.

# sir-deenicus/#Probabilistic Programming in the Browser.md Last active Dec 29, 2018

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

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!

 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)
 //////////////////////// 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)) //////////////////////////
 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))