Randomness in Haskell01 Oct 2016
Randomness is a constant nuisance point for Haskell beginners who may be coming from a language like Python or R. While in Python you can just get away with something like:
In : numpy.random.rand(3) Out: array([ 0.61426175, 0.05309224, 0.38861597])
or in R:
> runif(3)  0.49473012 0.68436352 0.04135914
In Haskell, the situation is more complicated. It’s not too much worse when you get the hang of things, but it’s certainly one of those things that throws beginners for a loop - and for good reason.
In this article I want to provide a simple guide, with examples, for getting started and becoming comfortable with randomness in Haskell. Hopefully it helps!
I’m writing this from a hotel during my partner’s birthday, so it’s being slapped together very rapidly with a kind of get-it-done attitude. If anything is unclear or you have any questions, feel free to shoot me a ping and I’ll try to improve it when I get a chance.
Randomness on Computers in General
Check out the R code I posted previously. If you just open R and type
runif(3) on your machine, then odds are you’ll get a different triple of
numbers than what I got above.
These numbers are being generated based on R’s global random number generator
(RNG), which, absent any fiddling by the user, is initialized as needed based
on the system time and ID of the R process. So: if you open up the R
interpreter and call
runif(3), then behind the scenes R will initialize the
RNG based on the time and process ID, and then use a particular algorithm to
generate random numbers based on that initialized value (called the ‘seed’).
These numbers aren’t truly random - they’re pseudo-random, which means they’re generated by a deterministic algorithm such that the resulting values appear random over time. The default algorithm used by R, for example, is the famous Mersenne Twister, which you can verify as follows:
> RNGkind()  "Mersenne-Twister" "Inversion"
You can also set the seed yourself in R, using the
set.seed function. Then
if you type something like
runif(3), R will use this initialized RNG rather
than coming up with its own seed based on the time and process ID. Setting
the seed allows you to reproduce operations involving pseudo-random numbers;
just re-set the seed and perform the same operations again:
> set.seed(42) > runif(3)  0.9148060 0.9370754 0.2861395 > set.seed(42) > runif(3)  0.9148060 0.9370754 0.2861395
(It’s good practice to always initialize the RNG using some known seed before running an experiment, simulation, and so on.)
So the big thing to notice here, in any case, is that R uses a global RNG.
It maintains the state of this RNG implicitly and behind the scenes. When you
runif(3), R consults this implicit RNG, gives you your pseudo-random
numbers based on its value, and updates the global RNG without you needing to
worry about any of this plumbing yourself. The same is generally true for
randomness in most programming languages - Python, C, Ruby, and so on.
Explicit RNG Management
But let’s come back to Haskell. Haskell, unlike R or Python, is purely-functional. State, or effects in general, are never implicit in the same way that R updates its global RNG. We need to either explicitly pass around a RNG ourselves, or at least allow some explicit monad to do it for us.
Passing around a RNG manually is annoying, so in practice this means everyone uses a monad to handle RNG state. This means that one needs to be comfortable working with monadic code in order to practically use random numbers in Haskell, which presents a big hurdle for beginners who may have been able to ignore monads thus far on their Haskell journey.
Let’s see what I mean by all of this by going through a few examples. Make sure you have stack installed, and then grab a few libraries that we’ll make use of in the remainder of this post:
$ stack install random mwc-random primitive
The Really Annoying Method - Manual RNG Management
Let me demonstrate the simplest conceptual method for dealing with random numbers: manually grabbing and passing around a RNG without involving any monads whatsoever.
First, open up GHCi:
$ stack ghci
And let’s also get some quick preliminaries out of the way:
Prelude> :set prompt "> " > import System.Random > import Control.Monad > let runif_pure = randomR (0 :: Double, 1) > let runif n = replicateM n (randomRIO (0 :: Double, 1)) > let set_seed = setStdGen . mkStdGen
> let rng = mkStdGen 42 > rng 43 1
We can use this thing to generate random numbers. A simple function to do that
randomR, which will generate pseudo-random values for some ordered
type in a given range. We’ll use the
runif_pure alias for it that we defined
previously, just to make things look similar to the previous R example and also
emphasize that this one is a pure function:
> runif_pure rng (1.0663729393723398e-2,2060101257 2103410263)
You can see that we got back a pair of values, the first element of which is
our random number
1.0663729393723398e-2. Cool. Let’s try to generate
> runif_pure rng (1.0663729393723398e-2,2060101257 2103410263)
Hmm. We generated the same number again. This is because the value of
hasn’t changed - it’s still the same value we made via
mkStdGen 42. Since
we’re using the same random number generator to generate a pseudo-random value,
we get the same pseudo-random value.
If we want to make new random numbers, then we need to use a different
generator. And the second element of the pair returned from our call to
runif_pure is exactly that - an updated RNG that we can use to generate
additional random numbers.
Let’s try that all again, using the generator we get back from the first function call as an input to the second:
> let (x, rng1) = runif_pure rng > x 1.0663729393723398e-2 > let (y, rng2) = runif_pure rng1 > y 0.9827538369038856
I mean.. sort of. It works and all, and it does constitute a general-purpose solution. But manually binding updated RNG states to names and swapping those in for new values is still pretty annoying.
You could also generate an infinite list of random numbers using the
randomRs function and just take from it as needed, but you still
probably need to manage that list to make sure you don’t re-use any numbers.
You kind of trade off managing the RNG for managing an infinite list of random
numbers, which isn’t much better.
The Less-Annoying Method - Get A Monad To Do It
The good news is that we can offload the job of managing the RNG state to a monad. I won’t actually explain how that works in detail here - I think most people facing this problem are initially more concerned with getting something working, rather than deeply grokking monads off the bat - so I’ll just claim that we can get a monad to handle the RNG state for us, and that will hopefully (mostly) suffice for now.
Still rolling with the
System.Random module for the time being, we’ll use the
runif alias for the
randomRIO function that we defined previously to
generate some new random numbers:
> runif 3 [0.9873934690803106,0.3794382930121829,0.2285653405908732] > runif 3 [0.7651878964537555,0.2623159001635825,0.7683468476766804]
Simpler! Notice we haven’t had to do anything with a generator manually - we
just ask for random numbers and then get them, just like in R. And if we want
to set the value of the RNG being used here, we can use the
function with an RNG that we’ve already created. Here let’s just use the
set_seed alias we defined earlier, to mimic R’s
> set_seed 42 > runif 3 [1.0663729393723398e-2,0.9827538369038856,0.7042944187434987] > set_seed 42 > runif 3 [1.0663729393723398e-2,0.9827538369038856,0.7042944187434987]
So things are similar to how they work in R here - we have a global RNG of
sorts, and we can set its state as desired using the
set_seed function. But
since this is Haskell, the effects of creating and updating the generator state
must still be explicit. And they are explicit - it’s just that they’re
explicit in the type of
> :t runif runif :: Int -> IO [Double]
runif returns a value that’s wrapped up in
IO. This is how we
indicate explicitly - at the type level - that something is being done with the
generator in the background.
IO is a monad, and it happens to be the thing
that’s dealing with the generator for us here.
What this means for you, the practitioner, is that you can’t just mix values of
a with values of type
IO a willy-nilly. You may be writing a
f with type
[Double] -> Double, where the input list of doubles is
intended to be randomly-generated. But if you just go ahead and generate a
xs of random numbers, they’ll have type
IO [Double], and you’ll stare
in confusion at some type error from GHC when you try to apply
Here’s what I mean. Take the example of just generating some random numbers and then summing them up. First, in R:
> xs = runif(3) > sum(xs)  1.20353
And now in Haskell, using the same mechanism we tried earlier:
> let xs = runif 3 > :t xs xs :: IO [Double] > sum xs <interactive>:16:1: No instance for (Num [Double]) arising from a use of ‘sum’ In the expression: sum xs In an equation for ‘it’: it = sum xs
This means that to deal with the numbers we generate, we have to treat them a little differently than we would in R, or compared to the situation where we were managing the RNG explicitly in Haskell. Concretely: if we use a monad to manage the RNG for us, then the numbers we generate will be ‘tagged’ by the monad. So we need to do something or other to make those tagged numbers work with ‘untagged’ numbers, or functions designed to work with ‘untagged’ numbers.
This is where things get confusing for beginners. Here’s how we could add up some random numbers in GHCi:
> xs <- runif 3 > sum xs 1.512024272587933
We’ve used the
<- symbol to bind the result of
runif 3 to the name
let xs = .... But this is sort of particular to running code in
GHCi; if you try to do this in a generic Haskell function, you’ll possibly wind
up with some more weird type errors. To do this in regular ol’ Haskell code,
you need to both use
<--style binding and also acknowledge the ‘tagged’
nature of randomly-generated values.
The crux is that, when you’re using a monad to generate random numbers in Haskell, you need to separate generating them from using them. Rather than try to explain what I mean here precisely, let’s rely on example, and implement a simple Metropolis sampler for illustration.
A Metropolis Sampler
The Metropolis algorithm will help you approximate expectations over certain probability spaces. Here’s how it works. Picture yourself strolling around some bumpy landscape; you want to walk around it in such a fashion that you visit regions of it with probability proportional to their altitude. To do that, you can repeatedly:
- Pick a random point near your current location.
- Compare your present altitude to the altitude of that point you picked. Calculate a probability based on their ratio.
- Flip a coin where the chance of observing a head is equal to that probability. If you get a head, move to the location you picked. Otherwise, stay put.
Let’s implement it in Haskell, using a monadic random number generator to do
so. This time we’re going to use
mwc-random - a more industrial-strength
randomness library that you can confidently use in production code.
mwc-random uses Marsaglia’s multiply-with-carry algorithm to generate
pseudo-random numbers. It requires you to explicitly create and pass a RNG to
functions that need to generate random numbers, but it uses a monad to update
the RNG state itself. This winds up being pretty nice; let’s dive in to see.
Create a module called
Metropolis.hs and get some imports out of the way:
module Metropolis where import Control.Monad import Control.Monad.Primitive import System.Random.MWC as MWC import System.Random.MWC.Distributions as MWC
The first thing we want to do is implement is point (1) from above:
Pick a random point near your current location.
We’ll just use a standard normal distribution of the appropriate dimension to do this - we just want to take a location, perturb it, and return the perturbed location.
propose :: [Double] -> Gen RealWorld -> IO [Double] propose location rng = traverse (perturb rng) location where perturb gen x = MWC.normal x 1 gen
So at finer detail: we’re walking over the coordinates of the current location
and generating a normally-distributed value centered at each coordinate. The
MWC.normal function will do this for a given mean and standard
deviation, and we can use the
traverse function to walk over each coordinate.
Note that we pass a
mwc-random RNG - the value with type
Gen RealWorld - to
propose function. We need to supply this generator anywhere we want to
generate random numbers, but we don’t need to manually worry about tracking and
updating its state. The
IO monad will do that for us. The resulting
randomly-generated values will be tagged with
IO, so we’ll need to deal with
Now let’s implement point (2):
Compare your present altitude to the altitude of that point you picked. Calculate a probability based on their ratio.
So, we need a function that will compare the altitude of our current point to the altitude of a proposed point and compute a probability from that. The following will do: it takes a function that will compute a (log-scale) altitude for us, as well as the current and proposed locations, and returns a probability.
moveProbability :: ([Double] -> Double) -> [Double] -> [Double] -> Double moveProbability altitude current proposed = whenNaN 0 (exp (min 0 (altitude proposed - altitude current))) where whenNaN val x | isNaN x = val | otherwise = x
Finally, the third step of the algorithm:
Flip a coin where the chance of observing a head is equal to that probability. If you get a head, move to the location you picked. Otherwise stay put.
So let’s get to it:
decide :: [Double] -> [Double] -> Double -> Gen RealWorld -> IO [Double] decide current proposed prob rng = do accept <- MWC.bernoulli prob rng return $ if accept then proposed else current
Here we need to flip a coin, so we require a source of randomness again. The
decide function thus takes another generator of type
Gen RealWorld that we
then supply to the
MWC.bernoulli function, and the result - the final
location - is once again wrapped in
This function clearly demonstrates the typical way that you’ll deal with random
numbers in Haskell code.
decide is a monadic function, so it proceeds using
do-notation. When you need to generate a random value - here we generate a
False value according to a Bernoulli distribution - you bind
the result to a name using the
<- symbol. Then afterwards, in the scope of
the function, you can use the bound value as if it were pure. But the entire
function must still return a ‘wrapped-up’ value that makes the effect of
passing the generator explicit at the type level; right here, that means that
the value will be wrapped up in
Putting Everything Together
The final Metropolis transition is a combination of steps one through three. We can put them together like so:
metropolis :: ([Double] -> Double) -> [Double] -> Gen RealWorld -> IO [Double] metropolis altitude current rng = do proposed <- propose current rng let prob = moveProbability altitude current proposed decide current proposed prob rng
metropolis is monadic, so we start off with a
do to make monadic
programming easy on us. Whenever we need a random value, we bind the result of
a random number-returning function using the
propose function returns a random location, so we bind its result to the
proposed using the
<- symbol. The
moveProbability function, on the
other hand, is pure - so we bind that using a
let prob = ... expression. The
decide function returns a random value, so we can just plop it right on the
end here. The entire result of the
metropolis function is random, so it is
wrapped up in
The result of
metropolis is just a single transition of the Metropolis
algorithm, which involves doing this kind of thing over and over. If we do
that, we observe a bunch of points that trace out a particular realization of a
Markov chain, which we can generate as follows:
chain :: Int -> ([Double] -> Double) -> [Double] -> Gen RealWorld -> IO [[Double]] chain epochs altitude origin rng = loop epochs [origin] where loop n history@(current:_) | n <= 0 = return history | otherwise = do next <- metropolis altitude current rng loop (n - 1) (next:history)
Now that we have our
chain function, we can use it to trace out a collection
of points visited on a realization of a Markov chain. Remember that we’re
supposed to be wandering over some particular abstract landscape; here, let’s
stroll over the one defined by the following function:
landscape :: [Double] -> Double landscape [x0, x1] = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1)
What we’ll now do is pick an origin to start from, wander over the landscape
for some number of steps, and then print the resulting realization of the
Markov chain to stdout. We’ll do all that through the following
main :: IO () main = do rng <- MWC.createSystemRandom let origin = [-0.2, 0.3] trace <- chain 1000 landscape origin rng mapM_ print trace
Running that will dump a trace to stdout. If you clean it up and plot it, you’ll see that the visited points have traced out a rough approximation of the landscape:
Hopefully this gives a broad idea of how to go about using random numbers in Haskell. I’ve talked about:
- Why randomness in Haskell isn’t as simple as randomness in (say) Python or R.
- How to handle randomness in Haskell, either by manual generator management or by offloading that job to a monad.
- How to get thing done when a monad manages the generator for you - separating random number generation from random number processing.
- Doing all the above with an industrial-strength RNG, using a simple Metropolis algorithm as an example.
Hopefully the example gives you an idea of how to work with random numbers in practice.
I’ll be the first to admit that randomness in Haskell requires more work than randomness in a language like R, which to this day remains my go-to interactive data analysis language of choice. Using randomness effectively in Haskell requires a decent understanding of how to work with monadic code, even if one doesn’t quite understand monads entirely yet.
What I can say is that when one has developed some intuition for monads - acquiring a ‘feel’ for how to work with monadic functions and values - the difficulty and awkwardness drop off a bit, and working with randomness feels no different than working with any other effect.
Happy generating! I’ve dumped the code for the Metropolis example into a gist.