Haskell package for sampling from standard probability distributions

Question

I would like to do some Monte Carlo analysis in Haskell. I would like to be able to write code like this:

do n <- poisson lambda
   xs <- replicateM n $ normal mu sigma
   return $ maximum xs

which corresponds to the stochastic model

n ~ Poisson(lambda)
for (i in 1:n)
  x[i] ~ Normal(mu, sigma)
y = max_{i=1}^n x[i]

I can see how to create the necessary random-sampling monad pretty easily. However, I would prefer not to have to implement samplers for all of the standard probability distributions. Is there a Haskell package that already has these implemented?

I have looked at package random-fu, which has been stalled at version 0.2.7 for three years, but I can't make sense of it; it depends on typeclasses MonadRandom and RandomSource, which aren't well explained.

I've also looked at package mwc-probability, but I can't make sense of it either -- it seems you have to already understand the PrimMonad and PrimState typeclasses.

Both of these packages strike me as having overly complex APIs, and seem to have entirely abandoned the standard random-number-generation framework of Haskell as found in System.Random.

Any advice would be appreciated.

Answer 1

Well, if you want to be able to write code like this:

do n <- poisson lambda
   xs <- replicateM n $ normal mu sigma
   return $ maximum xs

then you presumably want to use random-fu:

import Control.Monad
import Data.Random
import Data.Random.Distribution.Poisson
import Data.Random.Distribution.Normal

foo :: RVar Double
foo = do
  n <- poisson lambda
  xs <- replicateM (n+1) $ normal mu sigma
  return $ maximum xs

  where lambda = 10 :: Double
        mu = 0
        sigma = 6

main :: IO ()
main = print =<< replicateM 10 (sample foo)

I'm not sure that lack of updates over the past three years should be a deciding factor. Have there really been that many exciting advances in the world of gamma distributions?

Actually, it looks like mwc-probability works about the same:

import Control.Monad
import System.Random.MWC.Probability

foo :: Prob IO Double
foo = do
  n <- poisson lambda
  xs <- replicateM (n+1) $ normal mu sigma
  return $ maximum xs

  where lambda = 10 :: Double
        mu = 0
        sigma = 6

main :: IO ()
main = do
  gen <- createSystemRandom
  print =<< replicateM 10 (sample foo gen)