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:

1. For each event e in the distribution
2. 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
3. Take the sequence of tuple sequences produced in 2 and flatten them
4. Isolate actual events and add up their individual probabilities, to compute the actual probability of each distinct event
5. 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 11-fold 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.