In Part 1, we cloned PFP, a library for computing with probability distributions. PFP represents a distribution as a list of possible values, each with an associated probability.

But in the real world, things aren't always so easy. What if we wanted to pick a random number between 0 and 1? Our previous implementation would break, because there's an infinite number of values between 0 and 1---they don't exactly fit in a list.

As it turns out, Sungwoo Park and colleagues found an elegant solution to this problem. They represented probability distributions as sampling functions, resulting in something called the λ ◯ calculus. (I have no idea how to pronounce this!)

With a little bit of hacking, we can use their sampling functions as a drop-in replacement for PFP.

A common interface

Since we will soon have two ways to represent probability distributions, we need to define a common interface.

type Weight = Float class ( Functor d , Monad d ) => Dist d where weighted :: [( a , Weight )] -> d a uniform :: Dist d => [ a ] -> d a uniform = weighted . map ( \ x -> ( x , 1 ))

The function uniform will create an equally-weighted distribution from a list of values. Using this API, we can represent a two-child family as follows:

data Child = Girl | Boy deriving ( Show , Eq , Ord ) child :: Dist d => d Child child = uniform [ Girl , Boy ] family :: Dist d => d [ Child ] family = do child1 <- child child2 <- child return [ child1 , child2 ]

Now, we need to implement this API two different ways: Once with lists, and a second time with sampling functions.

Finite distibutions, revisted

The monad from part 1 now becomes:

type FDist = PerhapsT ( [] ) instance Dist FDist where weighted [] = error "Empty distribution" weighted xws = PerhapsT ( map weight xws ) where weight ( x , w ) = Perhaps x ( P ( w / sum )) sum = foldl' ( + ) 0 ( map snd xws ) exact :: FDist a -> [ Perhaps a ] exact = runPerhapsT

We can run this as follows:

> exact family [ Perhaps [ Girl , Girl ] 25.0 % , Perhaps [ Girl , Boy ] 25.0 % , Perhaps [ Boy , Girl ] 25.0 % , Perhaps [ Boy , Boy ] 25.0 % ]

Sampling functions

This second part is also easy. First, we define a type Rand a , which represents a value of type a computed using random numbers.

We implement Rand by reduction to Haskell's IO monad. (There's any number of other ways to do this, of course.)

newtype Rand a = Rand { runRand :: IO a } randomFloat :: Rand Float randomFloat = Rand ( getStdRandom ( random ))

Next, we make Rand a trivial, sequential monad:

instance Functor Rand where fmap = liftM instance Monad Rand where return x = Rand ( return x ) r >>= f = Rand ( do x <- runRand r runRand ( f x ))

We can turn any FDist into a Rand by picking an element at random:

instance Dist Rand where weighted = liftF . weighted liftF :: FDist a -> Rand a liftF fdist = do n <- randomFloat pick ( P n ) ( runPerhapsT fdist ) pick :: Monad m => Prob -> [ Perhaps a ] -> m a pick _ [] = fail "No values to pick from" pick n (( Perhaps x p ) : ps ) | n <= p = return x | otherwise = pick ( n - p ) ps

Of course, a single random sample won't do us much good. We need functions for repeated sampling and for displaying the results:

sample :: Rand a -> Int -> Rand [ a ] sample r n = sequence ( replicate n r ) sampleIO r n = runRand ( sample r n ) histogram :: Ord a => [ a ] -> [ Int ] histogram = map length . group . sort

And sure enough, we get a nice, even distribution:

> liftM histogram ( sampleIO family 1000 ) [ 243 , 242 , 254 , 261 ]

That's it! Two monads for the price of one.

Part 3: Coming soon.