How to Haskell – developing a Markov chain

One of my favorite things about Tidal is that it’s inside of Haskell, a well-developed full-blown programming language, so all sorts of tricks are available.  After some discussion on Slack I developed code for making Markov chains to use in Tidal, and here I’ll try to describe the process.

First – what’s a Markov chain? Basically it’s some kind of random sequence where each step depends only on the step immediately before it. A common way to express this is to imagine several possible “states” for each step, and then write a “transition matrix” which gives the odds of going from one state to another.

How should we represent all this in Tidal/Haskell?  We ultimately want to have the chain as a Pattern, but it’ll probably be easier to start out with an ordinary list, and then use the listToPat function to turn it into a Pattern.  Let’s number the states with integers, starting at zero.  Then the chain itself will be a list of integers.  The transition matrix is a list of lists of floating point numbers – we’ll use doubles. Then one way of writing things is to have a function that takes the transition matrix and a chain as arguments, and returns a new chain which is just the old one with one new element added. Here’s just the signature:

markovStep :: [[Double]] -> [Int] -> [Int]

Next let’s start the basic structure of what the function will have to do. Let’s call the transition matrix tp and the input chain xs. We need to look at the top element of the list to see what the “from” state is, and find that row of the transition matrix: tp!!(head xs). Next we need to generate a random number r (we’ll worry about how to actually do this later) and somehow use it to choose what the “to” state for the next step of the chain should be.

For example, if the row of the matrix were [0.20, 0.30, 0.50], this would mean 20% chance of transition to state zero, 30% chance to state one, and 50% chance to state two. So if our random number (between zero and one) were between 0 and 0.2, we’d pick state zero, between 0.2 and 0.5 state one, etc. This is easier to write if we take the cumulative sum of the row, and there’s a handy Haskell function we can use for that: scanl1 (+). That’ll turn the row into [0.20, 0.50, 1.00], and now when we have a random number we just need to find the index of the first element in the row larger than our random number. That’ll be the state the chain jumps to next.

There’s a handy Haskell function for finding things in lists too: findIndex. While it doesn’t quite work yet, we have most of our function:

markovStep tp xs = (findIndex (r <=) $ scanl1 (+) (tp!!(head xs))) : xs

(the : xs at the end adds the existing chain after the new element we just made)

One problem is that findIndex doesn’t actually return an Int. It returns a Maybe Int, because it might not find anything.  If our transition probabilities for each row add to one we’ll always be safe, but if we mess up or make a typo things might go wrong.  We have two options, (1) choose a safe default or (2) throw an exception and halt.  For the first option we could use fromMaybe 0 (to default to state zero), and for the second we’d use fromJust. Both require us to import the Data.Maybe module. Since I like to immediately know when I’ve made a math mistake and don’t mind if my performance grinds to a halt, I’m choosing the second option:

import Data.Maybe

markovStep tp xs = (fromJust $ findIndex (r <=) $ scanl1 (+) (tp!!(head xs))) : xs

The final thing to take care of is the random number. While it’s not quite the intended use, we can use the timeToRand function, which takes a rational number as an argument and hashes it to a pseudorandom result. We just have to make sure we get a different number for each step of the chain by using the length of the chain, and do type conversion as necessary:

import Data.Maybe

markovStep tp xs = (fromJust $ findIndex (r <=) $ scanl1 (+) (tp!!(head xs))) : xs 
  where r = timeToRand $ fromIntegral $ length xs

The last touch to make this a bit more usable is to wrap markovStep in a function that will repeat the step process to build a chain n numbers long from an initial state xi, flip it around so the most recent element is at the end, and turn the whole thing into a pattern with one number per cycle.

import Data.Maybe

markovStep :: [[Double]] -> [Int] -> [Int]
markovStep tp xs = (fromJust $ findIndex (r <=) $ scanl1 (+) (tp!!(head xs))) : xs 
  where r = timeToRand $ fromIntegral $ length xs
markovPat n xi tp = slow (fromIntegral n) $ listToPat 
  $ reverse $ (iterate (markovStep tp) [xi]) !! (n-1)

Now it’s ready to use in Tidal:

d1 
 $ fast 4 
 $ n (markovPat 64 0 [[0.3, 0.4, 0.3], [0.3, 0, 0.7], [0.7, 0.3, 0]])
 # s "arpy"

The chain states are always integers starting at zero, and that seems rather limiting.  But we can treat the chain as an index to some other array of our choice, which needn’t even be numerical:

d1
 $ fast 4
 $ n (fmap ([0,3,5]!!) $ markovPat 64 0 [[0.3, 0.4, 0.3], [0.3, 0, 0.7], [0.7, 0.3, 0]])
 # s "arpy"

d2
 $ fast 8
 $ s (fmap (["hc:2", "cr"]!!) $ markovPat 128 1 [[0.88, 0.12], [0.7, 0.3]])
 # cut 1

I’m sure there are still improvements we could make, but this is a pretty fun start!

Leave a Reply

Your email address will not be published. Required fields are marked *