[Haskell-cafe] Splittable random numbers

[Richard Senington]

On 11/11/10 21:34, Luke Palmer wrote:
> On Thu, Nov 11, 2010 at 3:13 AM, Richard Senington<sc06r2s at leeds.ac.uk>  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 at haskell.org
>> http://www.haskell.org/mailman/listinfo/haskell-cafe
>>
>>
>>      
>