haskell / mwc-random

A very fast Haskell library for generating high quality pseudo-random numbers.
http://hackage.haskell.org/package/mwc-random
BSD 2-Clause "Simplified" License
54 stars 25 forks source link

Something like `uniformR` for Natural and Integer. #51

Closed Ericson2314 closed 4 years ago

Ericson2314 commented 8 years ago

Of course we cannot uniformly sample from infinite types, but we can sample from a finite portion of them.

Ericson2314 commented 8 years ago

Here is a rather verbose implementation:

uniformNatR :: forall m. PrimMonad m => Natural -> Gen (PrimState m) -> m Natural
uniformNatR upperBound gen = do
  let numBits = length $ takeWhile (/= zeroBits) $ iterate (\n -> shiftR n 1) upperBound
  let loop = do
        bits <- help numBits
        if bits > upperBound
          then loop
          else return bits
  loop
  where help :: Int -> m Natural
        help numBits =
          if numBits == 0
          then return 0
          else do
            r <- uniform gen :: m Word64
            if numBits > 64
              then do
                upper <- help $ numBits - 64
                return $ fromIntegral r .|. shiftL upper 64
              else return $ fromIntegral $ r .&. (shiftR (complement zeroBits) $ 64 - numBits)
Shimuuar commented 8 years ago

I agree that such function would be useful. I haven't understood your implementation but I would expect it to be buggy. Sampling from non-power of two ranges is not entirely trivial and involves rejection sampling. One cannot fit such range integral number of times into PRNG output.

It also raises question whether Variate type class should be split in two.

Ericson2314 commented 8 years ago

I would not bet money on my implementation, but it does use rejection sampling for that. numBits is the power of two least upper bound, help samples [0, 2^numBits), loop does the rejection sampling. Sorry about the lack of clarity / unclear identifiers / etc, it was just something I through together quickly.

I do think the Variate class should be split. I suppose writing a polymorphic instance of it from two new classes would ease the transition.

Shimuuar commented 8 years ago

I could be wrong and your function may be correct after all I only spend few minutes looking at it.

One thing that should be decided is how to add such function. It basically asks to be abstracted in type class. Question is how to arrange these type classes. Currently we have:

class Variate a where
    uniform :: (PrimMonad m) => Gen (PrimState m) -> m a
    uniformR :: (PrimMonad m) => (a,a) -> Gen (PrimState m) -> m a

uniform assumes that a is finite and uniformR that a is ordered. Integer is good example where we can implement uniformR but not uniform. But are there any good example where uniform could be defined but uniformR could not? If so neither type class for uniform nor uniformR could be superclass for other.

Ericson2314 commented 8 years ago

No worries, I don't like the style of my code above nor find it very clear.

Every finite type can be ordered, but yes the ordering may not be natural. I agree two classes with neither being a superclass of the other is the way to go.

Would be nice to add an Ord bound for the second class, perhaps.

class Foo a where
    uniform :: (PrimMonad m) => Gen (PrimState m) -> m a

class Ord a => Bar a where
    uniformR :: (PrimMonad m) => (a,a) -> Gen (PrimState m) -> m a

-- backwards compat
{-# DEPRECATED "This is inflexible" #-}
instance Foo a => Bar a => Variate a
    ....
Shimuuar commented 4 years ago

random-1.2 now has type class UniformRange which provides exactly this functionality. mwc-random will provide necessary StatefulGen instance in next release thus closing