Modelling Stochastically Independent Processes with F# Computation Expressions: Part 2
Code for this entry: here.
In this first part, we started shamelessly plagiarizing creatively reproducing a Clojure monad that does a great job describing stochastically independent processes.
After a few trials I became convinced, that the approach in the blog post mentioned above was indeed optimal: a “natural” approach of modelling an event as a (value, probability) tuple is flawed in a sense that the probability needs to always be carried around and that gets in the way of describing the process in a declarative way.
Auxiliaries
I have defined a couple of auxiliary functions, just like in the Clojure case:
type Microsoft.FSharp.Collections.Map<'Key, 'Value when 'Key : comparison> with static member mapValues (f : 'Value > 'U) x = x > Map.toSeq > Seq.map (fun (k, v) > (k, f v)) let uniform (v : int seq) = let v = v > Seq.toList v > Seq.map (fun i > i, 1N / BigRational.FromInt v.Length) > Map.ofSeq
I have extended the Map type (because I can), to make the code below even more readable.
Computation expression
type IndEventBuilder () = member this.Zero () = Map.empty.Add(0, 0.) member this.Return (e : 'a) : Map<'a, BigRational> = Map.empty.Add(e, 1N) member this.Bind (m : Map<'a, BigRational>, fn : 'a > Map<'b, BigRational>) : Map<'b, BigRational> = m > Map.toSeq > Seq.map (fun (key, value) > (fn key) > Map.mapValues (fun x > x * value)) > Seq.concat > Seq.groupBy(fun (key, value) > key) > Seq.map(fun (key, sq) > (key, sq > Seq.map (fun (key, value) > value) > Seq.sum)) > Map.ofSeq let independent = IndEventBuilder()
Zero and Return functions are necessary: Zero because I want to use “if” statements within the computation expression, and Return is always necessary as a terminator of binding chains: Bind(.., Bind(.., (.., Return))). In that precise sense, Return makes the event marked by its argument occur with the probability 1.
Bind takes a distribution and unwraps it by making each event in the distribution occur together with what comes next:
 For each event e in the distribution
 Take the event and map it to the distribution that comes afterwards by multiplying the probability of this event by the probability of every event in the distribution that follows
 Take the sequence of tuple sequences produced in 2 and flatten them
 Isolate actual events and add up their individual probabilities, to compute the actual probability of each distinct event
 Turn the resulting sequence into a map.
Red Dog
Here is a slightly more interesting application: modelling the Red Dog game. Here is the “naive” version:
let simpleRedDog = independent { let! card1 = cardProb let! card2 = cardProb let! card3 = cardProb let minCard = min card1 card2 let maxCard = max card1 card2 if card1 = card2 then if card2 = card3 then return 10 else return 0 elif minCard + 1 = maxCard then return 0 elif minCard < card3 && card3 < maxCard then return 1 else return 1 }
This is the version from the Clojure blog (the rules are explained there as well). All we need to do is understand the rules and write them out. The only problem is that this implementation assumes cards are being replaced in the deck while the hand is played. Easy to fix, however, it actually gives close to correct results (with probabilities being small enough, rigorously accounting for replacements does not matter much):
val simpleRedDog : Map<int,BigRational> = map [(1, 88/169N); (0, 36/169N); (1, 44/169N); (10, 1/169N)]
Here is the fixed version, modelling the process without replacements:
let redDog = let removeCard (v : BigRational) (i : BigInteger) = if v = 0N then v else let mult = BigInteger range / v.Denominator BigRational.FromBigInt (v.Numerator * mult  i) / BigRational.FromBigInt (v.Denominator * mult  BigInteger.One) let distWithRemoved card dist = dist > Map.map (fun key v > if key <> card then removeCard v BigInteger.Zero else removeCard v BigInteger.One) independent { let! card1 = cardProb let removedCard = distWithRemoved card1 cardProb let! card2 = removedCard let! card3 = distWithRemoved card2 removedCard let minCard = min card1 card2 let maxCard = max card1 card2 if card1 = card2 then if card2 = card3 then return 10 else return 0 elif minCard + 1 = maxCard then return 0 elif minCard < card3 && card3 < maxCard then return 1 else return 1 } let expect dist = dist > Map.toSeq > Seq.sumBy (fun (k, v) > float k * BigRational.ToDouble v)
We simply adjust for the removed card in the distWithRemoved function. (Expressing this in Clojure is harder because monads are not part of the language, and so succinctness characteristic of Clojure will be affected).
The results:
val redDog : Map<int,BigRational> = map [(1, 8624/16575N); (0, 1112/5525N); (1, 352/1275N); (10, 1/425N)]
These are not all that different from the approximation. Even though the chance of an 11fold gain appears to be about 3 times higher in the simplified version, the real difference between 0.59% and 0.23% (actual) is probably not all that great. Expected gains are also pretty close:
> expect redDog;; val it : float = 0.220693816 > expect simpleRedDog;; val it : float = 0.201183432
BUGS/JAGS
BUGS & JAGS are famous enough with Bayesian statisticians to not require detailed description. Usage is based on a very cool principle of declarative programming: simply describe the distributions and let the system do the sampling behind the scenes. While the creators of these packages probably were not thinking in terms of functional programming, the results, in view of everything described above, are if not close, at least strikingly similar: we extend the language to model the process chain in a declarative way.

August 17, 2015 at 8:27 pmModelling Stochastically Independent Processes with F# Computation Expressions: Part 2  Dinesh Ram Kali.