Posts Tagged ‘random’

Programming Praxis – Two Random Exercises

August 17, 2012

In today’s Programming Praxis exercise, our goal is to use convert two random number generators to a different range. Specifically, we have to convert an rng that produces numbers in the closed interval [1,3] to one that has the interval [1,9] and one with interval [1,5] to [1,7]. In both cases the resulting rngs must maintain an equal distribution across the interval. Let’s get started, shall we?

Some imports:

import Control.Applicative
import Control.Monad
import Data.List
import System.Random

These are the two functions we can use as imports.

rand3, rand5 :: IO Int
rand3 = randomRIO (1,3)
rand5 = randomRIO (1,5)

To test, we use a distribution function to see how many of each element occur in a sequence.

dist :: (Eq a, Ord a) => [a] -> [(a, Int)]
dist = map (\x -> (head x, length x)) . group . sort

Let’s start with the easiest one. Since 9 = 3 * 3, all we need to do scale up one input by a factor of 3 and “interpolate” using a second one. The reason for taking a source parameter is to facilitate testing: you can either feed it an IO Int for a random result, or you can give it a [Int] to see all possible results. We’ll get to why the testing is important in a bit.

rand9 :: (Applicative f, Num a, Enum a) => f a -> f a
rand9 src = (+) . (*3) . pred <$> src <*> src

The second one is more complicated, since 5*5 is not a multiple of 7. The correct solution is to do the same thing as above using a subrange that is a multiple of 7, i.e. 21 of the 25 numbers and rerolling when one of the four values outside the accepted range is chosen. The problem with this approach is that it isn’t guaranteed to terminate, since there is nothing preventing a true random number generator from generating the number 23 a few trillion times in a row. Granted, we’re only using pseudo-random number generators, so it probably will never come up. Here’s the implementation:

rand7 :: (Applicative m, Monad m, Integral a) => m a -> m a
rand7 src = (+) . (*5) . pred <$> src <*> src >>= \n -> 
    if n > 21 then rand7 src else return (mod n 7 + 1)

Since I didn’t like the chance of a non-terminating algorithm I came up with the following algorithm: sum seven random numbers in the range 1 to 5, and give the result modulo 7. Or in code:

rand7approx :: (Functor f, Monad f, Integral a) => f a -> f a
rand7approx src = succ . (`mod` 7) . sum <$> replicateM 7 src

Printing the distribution looked perfectly fine, with all of them getting a little over 14%. But something felt off to me. The code was shorter than the version with rerolling and didn’t suffer the inifinite loop problem, so there had to be reason it wasn’t the provided solution. I changed all the functions to take a source argument rather than calling rand3 and rand5 directly and printed the theoretical distributions as well (except for rand7, since the distribution for that one is non-terminating). This revealed that there is a 0.3% difference between the most and least common values. In retrospect, this should’ve been obvious, as 5^7 mod 7 isn’t equal to 0, which means it’s impossible for all 7 values to get the same distribution. Some more reading up revealed that there is indeed no solution that both provides an equal distribution and is guaranteed to terminate, since 7 is an infinitely repeating fraction in base 5.

Here are the tests I used:

main :: IO ()
main = do print $ dist (rand9 [1..3]) == zip [1..9] [1,1..]
          print . dist =<< replicateM 100000 (rand9 rand3)
          print . dist =<< replicateM 100000 (rand7 rand5)
          print $ dist (rand7approx [1..5])
          print . dist =<< replicateM 100000 (rand7approx rand5)

Programming Praxis – Rule 30 RNG

April 29, 2011

In today’s Programming Praxis exercise, our goal is to write a random number generator based on the Rule 30 cellular automaton. Let’s get started, shall we?

Some imports:

import Control.Monad.State
import Data.Bits
import Data.Word

First we define the amount of bits in our RNG state.

size :: Int
size = 7 * 43

We’re going to be needing a function to convert a list of bits to a number.

fromBits :: (Bits a, Integral a) => [Bool] -> a
fromBits = foldl (\a x -> shift a 1 .|. fromIntegral (fromEnum x)) 0

The step function replaces the state with the next generation of cells for a given rule and returns the middle bit.

step :: Monad m => Int -> StateT Integer m Bool
step r = modify (\s -> let b i = testBit s (mod i size) in
    fromBits [ testBit r $ fromBits [b (n+1), b n, b (n - 1)]
             | n <- [size - 1, size - 2..0]]) >> gets (`testBit` div size 2)

To produce a random number, we take the middle bit of the desired number of generations and convert that to an unsigned integer.

randomBits :: Monad m => Int -> Int -> StateT Integer m Word
randomBits n r = return . fromBits =<< replicateM n (step r)

To test the algorithm, we see if we get the correct result when we only set the bit in the middle. Afterwards, we produce a bunch of 8-bit numbers.

main :: IO ()
main = evalStateT (do r <- randomBits 32 30
                      liftIO $ print (r == 3112904540)
                      liftIO . print =<< replicateM 100 (randomBits 8 30)
                  ) $ bit (div size 2)

Looks like everything is working properly.

Programming Praxis – Fortune

April 5, 2011

In today’s Programming Praxis exercise, our goal is to implement the fortune Unix command line tool. Let’s get started, shall we?

Some imports:

import Control.Monad
import System.Environment
import System.Random

Since we wrote some code in a previous exercise to select a random item from a list we can just use that here again:

chance :: Int -> Int -> a -> a -> IO a
chance x y a b = fmap (\r -> if r < x then a else b) $ randomRIO (0, y-1)

fortune :: [a] -> IO a
fortune = foldM (\a (n,x) -> chance 1 n x a) undefined . zip [1..]

All that’s left to do to implement the fortune program is to read the appropriate file and choose a random line.

main :: IO ()
main = putStrLn =<< fortune . lines =<<
       readFile . head . (++ ["fortunes.txt"]) =<< getArgs

Four lines versus the 23 of the C version. Not bad.

Programming Praxis – Two Random Selections

December 10, 2010

In today’s Programming Praxis exercise we have to implement two algorithms that select random items from a list in linear time. Let’s get started, shall we?

Some imports:

import Control.Monad
import Data.List
import System.Random

I found myself doing the same thing in both functions so I factored it out. This function gives you an x in y chance of choosing a instead of b.

chance :: Int -> Int -> a -> a -> IO a
chance x y a b = fmap (\r -> if r < x then a else b) $ randomRIO (0, y - 1)

The first algorithm (selecting one item at random from a list) can be done by folding over the list, with each item having a decreasing chance of becoming the new choice.

fortune :: [a] -> IO a
fortune = foldM (\a (n, x) -> chance 1 n x a) undefined . zip [1..]

The second algorithm (selecting m items from a list of integers) is pretty much the same, except now we also have to keep track of how many items we selected. This version always goes through the entire list rather than stopping when m items have been selected, but since it still runs in O(n) and the resulting code is cleaner I went with this version.

sample :: Int -> Int -> IO [Int]
sample m n = fmap snd $ foldM (\(m', a) x -> chance m' x
    (m' - 1, x:a) (m', a)) (m, []) [n, n-1..1]

With random algorithms it’s always a good idea to check the distribution of the results, as was proven again today because it revealed a bug in my code.

main :: IO ()
main = do let dist n f = mapM_ (\x -> print (length x, head x)) .
                         group . sort . concat =<< replicateM n f
          dist 10000 . fmap return $ fortune ["rock", "paper", "scissors"]
          dist 10000 $ sample 6 43

The frequency distribution is pretty much equal and they sum up to the correct amount, so everything seems to be working correctly.

Programming Praxis – E

August 13, 2010

In today’s Programming Praxis exercise our task is to determine how many random numbers between 0 and 1 we have to sum on average to exceed 1. Let’s get started, shall we?

Some imports:

import Control.Monad
import System.Random

First, the algorithm itself, which counts the number of required additions.

f :: Float -> IO Int
f s = if s > 1 then return 0 else fmap succ . f . (s +) =<< randomRIO (0, 1)

To get an average, we run the algorithm a number of times and take the average.

test :: Int -> IO Float
test n = fmap ((/ toEnum n) . toEnum . sum) $ replicateM n (f 0)

As usual, a test to see if everything works:

main :: IO ()
main = print =<< test 100000

Three sequential test produce 2.72085, 2.71691 and 2.71838. Certainly in the right ballpark of the expected result of e, with the last one accurate to three places behind the comma, so it looks like everything’s okay.