content: ---
  title: Picking a random element from a frequency list
  description: A mildly interesting algorithm I used in Lorem Markdownum
  tags: haskell
  ---
  
  Introduction
  ============
  
  A week or so ago, I wrote [Lorem Markdownum]: a small webapp to generate random
  text (like the many [lorem ipsum] generators out there), but in [markdown]
  format.
  
  [Lorem Markdownum]: http://jaspervdj.be/lorem-markdownum/
  [lorem ipsum]: http://www.lipsum.com/
  [markdown]: http://daringfireball.net/projects/markdown/
  
  This blogpost explains a mildly interesting algorithm I used to pick a random
  element from a frequency list. It is written in Literate Haskell so you should
  be able to drop it into a file and run it -- the raw version can be found
  [here](https://raw.github.com/jaspervdj/jaspervdj/master/posts/2013-11-21-random-element-frequency-list.lhs).
  
  > import           Data.List       (sortBy)
  > import           Data.Ord        (comparing)
  > import qualified Data.Map.Strict as M
  > import           System.Random   (randomRIO)
  
  The problem
  ===========
  
  Lorem ipsum generators usually create random but realistic-looking text by using
  sample data, on which a model is trained. A very simple example of that is just
  to pick words according to their frequencies. Let us take some sample data from
  a song that [fundamentally changed](http://www.youtube.com/watch?v=EIyixC9NsLI)
  the music industry in the early 2000s:
  
  <blockquote>
  Badger badger badger</br>
  Mushroom mushroom</br>
  Badger badger badger</br>
  Panic, a snake</br>
  Badger badger badger</br>
  Oh, it's a snake!</blockquote>
  
  This gives us the following frequency list:
  
  > badgers :: [(String, Int)]
  > badgers =
  >     [ ("a",        2)
  >     , ("badger",   9)
  >     , ("it's",     1)
  >     , ("mushroom", 2)
  >     , ("oh",       1)
  >     , ("panic",    1)
  >     , ("snake",    2)
  >     ]
  
  The sum of all the frequencies in this list is 18. This means that we will e.g.
  pick "badger" with a chance of 9/18. We can naively implement this by expanding
  the list so it contains the items in the given frequencies and then picking one
  randomly.
  
  > decodeRle :: [(a, Int)] -> [a]
  > decodeRle []            = []
  > decodeRle ((x, f) : xs) = replicate f x ++ decodeRle xs
  
  > sample1 :: [(a, Int)] -> IO a
  > sample1 freqs = do
  >     let expanded = decodeRle freqs
  >     idx <- randomRIO (0, length expanded - 1)
  >     return $ expanded !! idx
  
  This is obviously extremely inefficient, and it is not that hard to come up with
  a better definition: we do not expand the list, and instead use a specialised
  indexing function for frequency lists.
  
  > indexFreqs :: Int -> [(a, Int)] -> a
  > indexFreqs _   [] = error "please reboot computer"
  > indexFreqs idx ((x, f) : xs)
  >     | idx < f     = x
  >     | otherwise   = indexFreqs (idx - f) xs
  
  > sample2 :: [(a, Int)] -> IO a
  > sample2 freqs = do
  >     idx <- randomRIO (0, sum (map snd freqs) - 1)
  >     return $ indexFreqs idx freqs
  
  However, `sample2` is still relatively slow when we our sample data consists of
  a large amount of text (imagine what happens if we have a few thousand different
  words). Can we come up with a better but still elegant solution?
  
  Note that lorem ipsum generators generally employ more complicated strategies
  than just picking a word according to the frequencies in the sample data.
  Usually, algorithms based on [Markov Chains] are used. But even when this is the
  case, picking a word with some given frequencies is still a subproblem that
  needs to be solved.
  
  [Markov Chains]: http://en.wikipedia.org/wiki/Markov_chain
  
  Frequency Trees
  ===============
  
  It is easy to see why `sample2` is relatively slow: indexing in a linked list is
  expensive. Purely functional programming languages usually solve this by using
  trees instead of lists where fast indexing is required. We can use a similar
  approach here.
  
  A leaf in the tree simply holds an item and its frequency. A branch also holds
  a frequency -- namely, the sum of the frequencies of its children. By storing
  this computed value, we will be able to write a fast indexing this method.
  
  > data FreqTree a
  >     = Leaf   !Int !a
  >     | Branch !Int (FreqTree a) (FreqTree a)
  >     deriving (Show)
  
  A quick utility function to get the sum of the frequencies in such a tree:
  
  > sumFreqs :: FreqTree a -> Int
  > sumFreqs (Leaf   f _)   = f
  > sumFreqs (Branch f _ _) = f
  
  Let us look at the tree for `badgers` (we will discuss how this tree is computed
  later):
  
  ![A nicely balanced tree for the badgers example](/images/2013-11-21-badgers-balanced.png)
  
  Once we have this structure, it is not that hard to write a faster indexing
  function, which is basically a search in a binary tree:
  
  > indexFreqTree :: Int -> FreqTree a -> a
  > indexFreqTree idx tree = case tree of
  >     (Leaf _ x)             -> x
  >     (Branch _ l r)
  >         | idx < sumFreqs l -> indexFreqTree idx                l
  >         | otherwise        -> indexFreqTree (idx - sumFreqs l) r
  
  > sample3 :: FreqTree a -> IO a
  > sample3 tree = do
  >     idx <- randomRIO (0, sumFreqs tree - 1)
  >     return $ indexFreqTree idx tree
  
  There we go! We intuitively see this method is faster since we only have to walk
  through a few nodes -- namely, those on the path from the root node to the
  specific leaf node.
  
  But how fast is this, exactly? This depends on how we build the tree.
  
  Well-balanced trees
  ===================
  
  Given a list with frequencies, we can build a nicely balanced tree (i.e., in the
  sense in which binary tries are balanced). This minimizes the longest path from
  the root to any node.
  
  We first have a simple utility function to clean up such a list of frequencies:
  
  > uniqueFrequencies :: Ord a => [(a, Int)] -> [(a, Int)]
  > uniqueFrequencies =
  >     M.toList . M.fromListWith (+) . filter ((> 0) . snd)
  
  And then we have the function that actually builds the tree. For a singleton
  list, we just return a leaf. Otherwise, we simply split the list in half, build
  trees out of those halves, and join them under a new parent node. Computing the
  total frequency of the parent node (`freq`) is done a bit inefficiently, but
  that is not the focus at this point.
  
  > balancedTree :: Ord a => [(a, Int)] -> FreqTree a
  > balancedTree = go . uniqueFrequencies
  >   where
  >     go []       = error "balancedTree: Empty list"
  >     go [(x, f)] = Leaf f x
  >     go xs       =
  >         let half     = length xs `div` 2
  >             (ys, zs) = splitAt half xs
  >             freq     = sum $ map snd xs
  >         in Branch freq (go ys) (go zs)
  
  Huffman-balanced trees
  ======================
  
  However, well-balanced trees might not be the best solution for this problem.
  It is [generally known] that few words in most natural languages are extremely
  commonly used (e.g. "the", "a", or in or case, "badger") while most words are
  rarely used.
  
  [generally known]: http://www.oxforddictionaries.com/us/words/the-oec-facts-about-the-language
  
  For our tree, it would make sense to have the more commonly used words closer to
  the root of the tree -- in that case, it seems intuitive that the *[expected]*
  number of nodes visited to pick a random word will be lower.
  
  [expected]: http://en.wikipedia.org/wiki/Expected_value
  
  It turns out that this idea exactly corresponds to a [Huffman tree]. In a
  Huffman tree, we want to minimize the expected code length, which equals the
  expected path length. Here, we want to minimize the expected number of nodes
  visited during a lookup -- which is precisely the expected path length!
  
  [Huffman tree]: http://en.wikipedia.org/wiki/Huffman_coding
  
  The algorithm to construct such a tree is surprisingly simple. We start out with
  a list of trees: namely, one singleton leaf tree for each element in our
  frequency list.
  
  Then, given this list, we take the two trees which have the lowest total sums of
  frequencies (`sumFreqs`), and join these using a branch node. This new tree is
  then inserted back into the list.
  
  This algorithm is repeated until we are left with only a single tree in the
  list: this is our final frequency tree.
  
  > huffmanTree :: Ord a => [(a, Int)] -> FreqTree a
  > huffmanTree = go . map (\(x, f) -> Leaf f x) . uniqueFrequencies
  >   where
  >     go trees = case sortBy (comparing sumFreqs) trees of
  >         []             -> error "huffmanTree: Empty list"
  >         [ft]           -> ft
  >         (t1 : t2 : ts) ->
  >             go $ Branch (sumFreqs t1 + sumFreqs t2) t1 t2 : ts
  
  This yields the following tree for our example:
  
  ![A tree built using the huffman algorithm](/images/2013-11-21-badgers-huffman.png)
  
  Is the second approach really better?
  =====================================
  
  Although Huffman trees are well-studied, for our example, we only *intuitively*
  explained why the second approach is *probably* better. Let us see if we can
  justify this claim a bit more, and find out *how much* better it is.
  
  The expected path length *L* of an item in a balanced tree can be very easily
  approached, since it is just a binary tree and we all know those (suppose *N* is
  the number of unique words):
  
  ![](/images/2013-11-21-expected-length-bal.png)
  
  However, if we have a tree we built using the `huffmanTree`, it is not that easy
  to calculate the expected path length. We know that for a Huffman tree, the path
  length should approximate the entropy, which, in our case, gives us an
  approximation for the path length for item with a specified frequency *f*:
  
  ![](/images/2013-11-21-length-huf.png)
  
  Where *F* is the total sum of all frequencies. If we assume that we know the
  frequency for every item, the expected path length is simply a weighted mean:
  
  ![](/images/2013-11-21-expected-length-huf.png)
  
  This is where it gets interesting. It turns out that the frequency of words in a
  natural language is a [well-researched](http://planetmath.org/ZipfsLaw)
  [topic](http://en.wikipedia.org/wiki/Zipf%27s_law), and predicted by something
  called *Zipf's law*. This law tells us that the frequency of an item *f* can be
  estimated by:
  
  ![](/images/2013-11-21-zipfs-law.png)
  
  Where *s* characterises the distribution and is typically very close to 1 for
  natural languages. *H* is the generalised [harmonic number]:
  
  [harmonic number]: http://en.wikipedia.org/wiki/Harmonic_number
  
  ![](/images/2013-11-21-harmonic-number.png)
  
  If we substitute in the definition for the frequencies into the formula for the
  expected path length, we get:
  
  ![](/images/2013-11-21-expected-length-huf-expanded.png)
  
  This is something we can work with! If we plot this for *s = 1*, we get:
  
  ![](/images/2013-11-21-graph-bal-huf.png)
  
  It is now clear that the expected path length for a frequency tree built using
  `huffmanTree` is expected to be significantly shorter than a frequency tree
  built using `balancedTree`, even for relatively small *N*. Yay! Since the
  algorithm now works, the conclusion is straightforward.
  
  Conclusion
  ==========
  
  Lorem markdownum constitit foret tibi Phoebes propior poenam. Nostro sub flos
  auctor ventris illa choreas magnis at ille. Haec his et tuorum formae [obstantes
  et viribus](http://www.wtfpl.net/) videret vertit, spoliavit iam quem neptem
  corpora calidis, in. Arcana ut puppis, ad agitur telum conveniant quae ardor?
  Adhuc [arcu acies corpore](http://haskell.org/) amplexans equis non velamina
  buxi gemini est somni.
  
  Thanks to [Simon Meier](https://github.com/meiersi), [Francesco
  Mazzoli](http://mazzo.li/) and some other people at
  [Erudify](http://www.erudify.com/) for the interesting discussions about this
  topic!
  