Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.