Re: [Haskell-cafe] Splittable random numbers

12 Nov 2010

      On 11/11/10 21:34, Luke Palmer wrote:
...
On Thu, Nov 11, 2010 at 3:13 AM, Richard Senington  wrote:
...
I got hold of, and looked through the paper suggested in the root of this
thread “Pseudo random trees in Monte-Carlo", and based upon this
I have thrown together a version of the binary tree based random number
generator suggested.
I would like to point out that I do not know very much about random number
generators, the underlying mathematics or any subsequent papers on this
subject, this is just a very naive implementation based upon this one paper.
As a question, the following code actually generates a stream of numbers
that is more random than I was expecting, if anyone can explain why I would
be very interested.
What do you mean more random than you were expecting?  Shouldn't they
be "maximally random"?
My issue is how it should react, given how the underlying data structure 
works.
If you just use this tree to generate numbers, you are taking 
Left,Left,Left .....
If you split the tree, you get Left and Right.

So, in my test code at the bottom I have taken a generator and then 
generated 10 numbers from it.
Then we split. The left hand branch (g1 in the test) should generate 
numbers that are just "tail.randoms $ gen"
but this is not what happens, at least not for raw integers.

I have been doing some more testing, and if you limit the range (0-1000 
seems to be stable) then it works as described above, however
if you use wider ranges, or even the maximum range, then the sequences 
do not match as expected.
This worries me, as one advantage of PRNGs (see paper as I am 
paraphrasing) is repeatability, or certain expected properties.

The underlying system is working, so I probably have the range or data 
type conversion wrong somewhere.
You can test the underlying tree like so.

rawTest :: LehmerTree->IO()
rawTest t = do print $ myTake 10 t
            print $ myTake 10 $ leftBranch t
            print $ myTake 10 $ rightBranch t
   where
     myTake 0 _ = []
     myTake x t = nextInt t : (myTake (x-1) (leftBranch t))
...
BTW, nice module.  Do you want to hackage it up?  If not, I will.
I would be happy to hackage it up, but I think this is a bit premature. 
I started to read a bit more about PRNGs, and I came across tests for
randomness. It seemed that a library of test systems for RandomGens 
would be quite cool, so I started coding last night. That is far too
premature to even post up here, but in short, this system gives some 
very odd results.

For example, mean averages (I tried medians but that did not tell me 
much, I am going to look at modals some time this weekend).

mean :: [Float]->Double
mean [] = error "empty list has no mean?"
mean xs = ((sum.(map realToFrac)) xs) / (fromIntegral.length $ xs)

rangedMeanTest :: RandomGen g=>g->Int->(Int,Int)->Double
rangedMeanTest g count range = let p = take count $ randomRs range g
                                in mean (map fromIntegral p)

So, I am testing discrete randomness, ints. We take a range we wish to 
generate (0-3 for example), and generate some number of test values (I 
used 1000).
This list is converted into floats, and averaged. We can then predict 
what we think the average should be, given that this is an unbiased 
uniform (or nearly uniform) system.

It does not give the results you would want. This may have something to 
do with picking "good" parameters for the mkLehmerTree function.
For example, using a random setup, I just got these results
result   expected              range
16.814  expected = 16.0  (1,31)
16.191  expected = 16.5  (1,32)
16.576  expected = 17.0  (1,33)
17.081  expected = 17.5  (1,34)
17.543  expected = 18.0  (1,35)

In short, I am worried by the properties of this random number 
generator. I propose improving the testing system, and then posting both 
the test suite and this random generator to
Hackage, unless you really want it up now in a very very preliminary form.

RS
...
...
import System.Random
data LehmerTree = LehmerTree {nextInt :: Int,
                               leftBranch :: LehmerTree,
                               rightBranch :: LehmerTree}
instance Show LehmerTree where
   show g = "LehmerTree, current root = "++(show $ nextInt g)
mkLehmerTree :: Int->Int->Int->Int->Int->Int->LehmerTree
mkLehmerTree aL aR cL cR m x0 = innerMkTree x0
   where
     mkLeft x = (aL * x + cL) `mod` m
     mkRight x = (aR * x + cR) `mod` m
     innerMkTree x = let l = innerMkTree (mkLeft x)
                         r = innerMkTree (mkRight x)
                     in LehmerTree x l r
mkLehmerTreeFromRandom :: IO LehmerTree
mkLehmerTreeFromRandom = do gen<-getStdGen
                             let a:b:c:d:e:f:_ = randoms gen
                             return $ mkLehmerTree a b c d e f
This can be pure:
mkLehmerTreeFromRandom :: (RandomGen g) =>  g ->  LehmerTree
...
instance RandomGen LehmerTree where
   next g = (fromIntegral.nextInt $ g, leftBranch g)
   split g = (leftBranch g, rightBranch g)
   genRange _ = (0, 2147483562) -- duplicate of stdRange
test :: IO()
test = do gen<-mkLehmerTreeFromRandom
           print gen
           let (g1,g2) = split gen
           let p = take 10 $ randoms gen :: [Int]
           let p' = take 10 $ randoms g1 :: [Int]
           -- let p'' = take 10 $ randoms g2 :: [Float]
           print p
           print p'
           -- print p''
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe