diff --git a/README.markdown b/README.markdown deleted file mode 100644 --- a/README.markdown +++ /dev/null @@ -1,29 +0,0 @@ -# Efficient, general purpose pseudo-random number generation -[![Build Status](https://travis-ci.org/Shimuuar/mwc-random.png?branch=master)](https://travis-ci.org/Shimuuar/mwc-random) -[![Build status](https://ci.appveyor.com/api/projects/status/4228vkxje4as3nhw/branch/master)](https://ci.appveyor.com/project/Shimuuar/mwc-random) - -This package provides the System.Random.MWC module, a Haskell library -for generating high-quality pseudo-random numbers in a space- and -time-efficient way. - - -# Get involved! - -Please report bugs via the -[github issue tracker](http://github.com/bos/mwc-random). - -Master git [git repository](http://github.com/bos/mwc-random): - -* `git clone git://github.com/bos/mwc-random.git` - -There's also a [Mercurial mirror](http://bitbucket.org/bos/mwc-random): - -* `hg clone http://bitbucket.org/bos/mwc-random` - -(You can create and contribute changes using either Mercurial or git.) - - -# Authors - -This library is written and maintained by Bryan O'Sullivan, -. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,23 @@ +# Efficient, general purpose pseudo-random number generation +[![Build +Status](https://github.com/haskell/mwc-random/workflows/Haskell-CI/badge.svg)](https://github.com/haskell/mwc-random/actions) + +This package provides the System.Random.MWC module, a Haskell library +for generating high-quality pseudo-random numbers in a space- and +time-efficient way. + + +# Get involved! + +Please report bugs via the +[github issue tracker](https://github.com/haskell/mwc-random/issues). + +Master [git repository](https://github.com/haskell/mwc-random): + +* `git clone https://github.com/haskell/mwc-random.git` + + +# Authors + +This library is written and maintained by Bryan O'Sullivan, +. diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -1,5 +1,6 @@ -{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts, - MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples +{-# LANGUAGE BangPatterns, CPP, DataKinds, DeriveDataTypeable, FlexibleContexts, + FlexibleInstances, MultiParamTypeClasses, MagicHash, Rank2Types, + ScopedTypeVariables, TypeFamilies, UnboxedTuples, TypeOperators #-} -- | -- Module : System.Random.MWC @@ -10,68 +11,124 @@ -- Stability : experimental -- Portability : portable -- --- Pseudo-random number generation. This module contains code for --- generating high quality random numbers that follow a uniform --- distribution. +-- Pseudo-random number generation using Marsaglia's MWC256, (also +-- known as MWC8222) multiply-with-carry generator, which has a period +-- of \(2^{8222}\) and fares well in tests of randomness. It is also +-- extremely fast, between 2 and 3 times faster than the Mersenne +-- Twister. There are two representation of generator: 'Gen' which is +-- generator that uses in-place mutation and 'Seed' which is immutable +-- snapshot of generator's state. -- --- For non-uniform distributions, see the --- 'System.Random.MWC.Distributions' module. -- --- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) --- multiply-with-carry generator, which has a period of 2^8222 and --- fares well in tests of randomness. It is also extremely fast, --- between 2 and 3 times faster than the Mersenne Twister. +-- == Initialization -- --- The generator state is stored in the 'Gen' data type. It can be --- created in several ways: +-- Generator could be initialized in several ways. One is to obtain +-- randomness from operating system using 'createSystemRandom', +-- 'createSystemSeed' or 'withSystemRandomST' (All examples assume +-- that @System.Random.Stateful@ is imported) -- --- 1. Using the 'withSystemRandom' call, which creates a random state. +-- >>> g <- createSystemRandom +-- >>> uniformM g :: IO Int +-- ... -- --- 2. Supply your own seed to 'initialize' function. +-- >>> withSystemRandomST $ \g -> uniformM g :: IO Int +-- ... -- --- 3. Finally, 'create' makes a generator from a fixed seed. --- Generators created in this way aren't really random. +-- Deterministically create generator from given seed using +-- 'initialize' function: -- --- For repeatability, the state of the generator can be snapshotted --- and replayed using the 'save' and 'restore' functions. +-- >>> import Data.Int +-- >>> import qualified Data.Vector.Unboxed as U +-- >>> import System.Random.Stateful +-- >>> g <- initialize $ U.fromList [1,2,3] +-- >>> uniformRM (1,200) g :: IO Int64 +-- 101 -- --- The simplest use is to generate a vector of uniformly distributed values: +-- Last way is to create generator with fixed seed which could be +-- useful in testing -- --- @ --- vs \<- 'withSystemRandom' . 'asGenST' $ \\gen -> 'uniformVector' gen 100 --- @ +-- >>> g <- create +-- >>> uniformM g :: IO Int +-- -8765701622605876598 -- --- These values can be of any type which is an instance of the class --- 'Variate'. -- --- To generate random values on demand, first 'create' a random number --- generator. +-- == Generation of random numbers -- --- @ --- gen <- 'create' --- @ +-- Recommended way of generating random numbers in simple cases like +-- generating uniformly distributed random number in range or value +-- uniformly distributed in complete type domain is to use +-- 'UniformRange' and 'Uniform' type classes. Note that while small +-- self-contained examples usually require explicit annotations +-- usually result type could be inferred. -- --- Hold onto this generator and use it wherever random values are --- required (creating a new generator is expensive compared to --- generating a random number, so you don't want to throw them --- away). Get a random value using 'uniform' or 'uniformR': +-- This example simulates 20 throws of fair 6-sided dice: -- --- @ --- v <- 'uniform' gen --- @ +-- >>> g <- create +-- >>> replicateM 20 $ uniformRM (1, 6::Integer) g +-- [3,4,3,1,4,6,1,6,1,4,2,2,3,2,4,2,5,1,3,5] -- --- @ --- v <- 'uniformR' (1, 52) gen --- @ +-- For generating full range of possible values one could use +-- 'uniformM'. This example generates 10 random bytes, or equivalently +-- 10 throws of 256-sided dice: +-- +-- >>> g <- create +-- >>> replicateM 10 $ uniformM g :: IO [Word8] +-- [209,138,126,150,165,15,69,203,155,146] +-- +-- There are special functions for generation of @Doubles@ and @Float +-- in unit interval: 'Random.uniformDouble01M', +-- 'Random.uniformDoublePositive01M', 'Random.uniformFloat01M', +-- 'Random.uniformFloatPositive01M': +-- +-- >>> uniformDouble01M =<< create +-- 0.5248103628705498 +-- >>> uniformFloat01M =<< create +-- 0.5248104 +-- +-- For normal distribution and others see modules +-- "System.Random.MWC.Distributions" and +-- "System.Random.MWC.CondensedTable". Note that they could be used +-- with any other generator implementing 'Random.StatefulGen' API +-- +-- There're special cases for generating random vectors and +-- bytestrings. For example in order to generate random 10-byte +-- sequences as unboxed vector or bytestring: +-- +-- >>> g <- create +-- >>> uniformVector g 10 :: IO (U.Vector Word8) +-- [209,138,126,150,165,15,69,203,155,146] +-- +-- >>> import qualified Data.ByteString as BS +-- >>> g <- create +-- >>> BS.unpack <$> uniformByteStringM 10 g +-- [138,242,130,33,209,248,89,134,150,180] +-- +-- Note that 'Random.uniformByteStringM' produces different result +-- from 'uniformVector' since it uses PRNG's output more efficiently. +-- +-- +-- == State handling +-- +-- For repeatability, the state of the generator can be snapshotted +-- and replayed using the 'save' and 'restore' functions. Following +-- example shows how to save and restore generator: +-- +-- >>> g <- create +-- >>> replicateM_ 10 (uniformM g :: IO Word64) +-- >>> s <- save g +-- >>> uniformM g :: IO Word32 +-- 1771812561 +-- >>> uniformM =<< restore s :: IO Word32 +-- 1771812561 module System.Random.MWC ( -- * Gen: Pseudo-Random Number Generators Gen , create , initialize - , withSystemRandom + , createSystemSeed , createSystemRandom - + , withSystemRandomST -- ** Type helpers -- $typehelp , GenIO @@ -80,6 +137,8 @@ , asGenST -- * Variates: uniformly distributed values + , Random.Uniform(..) + , Random.UniformRange(..) , Variate(..) , uniformVector @@ -89,7 +148,8 @@ , toSeed , save , restore - + -- * Deprecated + , withSystemRandom -- * References -- $references ) where @@ -98,29 +158,33 @@ #include "MachDeps.h" #endif -import Control.Monad (ap, liftM, unless) -import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO) -import Control.Monad.ST (ST) +import Control.Monad (unless) +import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO, stToPrim) +import Control.Monad.ST (ST,runST) import Data.Bits ((.&.), (.|.), shiftL, shiftR, xor) import Data.Int (Int8, Int16, Int32, Int64) -import Data.IORef (atomicModifyIORef, newIORef) +import Data.IORef (IORef, atomicModifyIORef, newIORef) import Data.Typeable (Typeable) import Data.Vector.Generic (Vector) import Data.Word +import Data.Kind import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Unboxed as I import qualified Data.Vector.Unboxed.Mutable as M import System.IO (hPutStrLn, stderr) import System.IO.Unsafe (unsafePerformIO) import qualified Control.Exception as E -#if defined(mingw32_HOST_OS) -import Foreign.Ptr -import Foreign.C.Types -#endif import System.Random.MWC.SeedSource - +import qualified System.Random.Stateful as Random +#if MIN_VERSION_random(1,3,0) +import Data.List.NonEmpty (NonEmpty(..), toList) +#endif --- | The class of types for which we can generate uniformly +-- | NOTE: Consider use of more principled type classes +-- 'Random.Uniform' and 'Random.UniformRange' instead. +-- +-- The class of types for which we can generate uniformly -- distributed random variates. -- -- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) @@ -225,52 +289,48 @@ {-# INLINE uniformR #-} instance Variate Int where -#if WORD_SIZE_IN_BITS < 64 +#if WORD_SIZE_IN_BITS == 32 uniform = uniform1 fromIntegral -#else +#elif WORD_SIZE_IN_BITS == 64 uniform = uniform2 wordsTo64Bit +#else +#error "Word size is not 32 nor 64" #endif uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word where -#if WORD_SIZE_IN_BITS < 64 +#if WORD_SIZE_IN_BITS == 32 uniform = uniform1 fromIntegral -#else +#elif WORD_SIZE_IN_BITS == 64 uniform = uniform2 wordsTo64Bit +#else +#error "Word size is not 32 nor 64" #endif uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} -{- -instance Variate Integer where - uniform g = do - u <- uniform g - return $! fromIntegral (u :: Int) - {-# INLINE uniform #-} --} - instance (Variate a, Variate b) => Variate (a,b) where - uniform g = (,) `liftM` uniform g `ap` uniform g - uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g + uniform g = (,) <$> uniform g <*> uniform g + uniformR ((x1,y1),(x2,y2)) g = (,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where - uniform g = (,,) `liftM` uniform g `ap` uniform g `ap` uniform g + uniform g = (,,) <$> uniform g <*> uniform g <*> uniform g uniformR ((x1,y1,z1),(x2,y2,z2)) g = - (,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` uniformR (z1,z2) g + (,,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g <*> uniformR (z1,z2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where - uniform g = (,,,) `liftM` uniform g `ap` uniform g `ap` uniform g - `ap` uniform g + uniform g = (,,,) <$> uniform g <*> uniform g <*> uniform g + <*> uniform g uniformR ((x1,y1,z1,t1),(x2,y2,z2,t2)) g = - (,,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` - uniformR (z1,z2) g `ap` uniformR (t1,t2) g + (,,,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g <*> + uniformR (z1,z2) g <*> uniformR (t1,t2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} @@ -385,13 +445,73 @@ {-# INLINE initialize #-} -- | An immutable snapshot of the state of a 'Gen'. -newtype Seed = Seed { - -- | Convert seed into vector. - fromSeed :: I.Vector Word32 - } - deriving (Eq, Show, Typeable) +newtype Seed = Seed (I.Vector Word32) + deriving (Eq, Show, Typeable) --- | Convert vector to 'Seed'. It acts similarily to 'initialize' and +-- | Convert seed into vector. +fromSeed :: Seed -> I.Vector Word32 +fromSeed (Seed v) = v + +-- | @since 0.15.0.0 +instance (s ~ PrimState m, PrimMonad m) => Random.StatefulGen (Gen s) m where + uniformWord32R u = uniformR (0, u) + {-# INLINE uniformWord32R #-} + uniformWord64R u = uniformR (0, u) + {-# INLINE uniformWord64R #-} + uniformWord8 = uniform + {-# INLINE uniformWord8 #-} + uniformWord16 = uniform + {-# INLINE uniformWord16 #-} + uniformWord32 = uniform + {-# INLINE uniformWord32 #-} + uniformWord64 = uniform + {-# INLINE uniformWord64 #-} +#if MIN_VERSION_random(1,3,0) + uniformByteArrayM isPinned n g = stToPrim (Random.fillByteArrayST isPinned n (uniform g)) + {-# INLINE uniformByteArrayM #-} +#else + uniformShortByteString n g = stToPrim (Random.genShortByteStringST n (uniform g)) + {-# INLINE uniformShortByteString #-} +#endif + +-- | @since 0.15.0.0 +instance PrimMonad m => Random.FrozenGen Seed m where + type MutableGen Seed m = Gen (PrimState m) + freezeGen = save +#if MIN_VERSION_random(1,3,0) + modifyGen gen@(Gen mv) f = do + seed <- save gen + case f seed of + (a, Seed v) -> a <$ G.copy mv v + overwriteGen (Gen mv) (Seed v) = G.copy mv v + +instance PrimMonad m => Random.ThawedGen Seed m where +#endif + thawGen = restore + +#if MIN_VERSION_random(1,3,0) +instance Random.SeedGen Seed where + type SeedSize Seed = 1032 -- == 4 * 258 + fromSeed64 seed64 = toSeed $ I.fromListN 258 + [ w32 + | !w64 <- toList seed64 + , !w32 <- [ fromIntegral (w64 `shiftR` 32) + , fromIntegral w64 ] + ] + toSeed64 vSeed = + let w32sToW64 :: Word32 -> Word32 -> Word64 + w32sToW64 w32u w32l = + (fromIntegral w32u `shiftL` 32) .|. fromIntegral w32l + v = fromSeed vSeed + evens = I.ifilter (\i _ -> even i) v + odds = I.ifilter (\i _ -> odd i) v + in case I.toList $ I.zipWith w32sToW64 evens odds of + [] -> + error $ "Impossible: Seed had an unexpected length of: " ++ show (I.length v) + x:xs -> x :| xs +#endif + +-- | Convert vector to 'Seed'. It acts similarly to 'initialize' and -- will accept any vector. If you want to pass seed immediately to -- restore you better call initialize directly since following law holds: -- @@ -401,42 +521,83 @@ -- | Save the state of a 'Gen', for later use by 'restore'. save :: PrimMonad m => Gen (PrimState m) -> m Seed -save (Gen q) = Seed `liftM` G.freeze q +save (Gen q) = Seed <$> G.freeze q {-# INLINE save #-} -- | Create a new 'Gen' that mirrors the state of a saved 'Seed'. restore :: PrimMonad m => Seed -> m (Gen (PrimState m)) -restore (Seed s) = Gen `liftM` G.thaw s +restore (Seed s) = Gen <$> G.thaw s {-# INLINE restore #-} --- | Seed a PRNG with data from the system's fast source of --- pseudo-random numbers (\"@\/dev\/urandom@\" on Unix-like systems or --- @RtlGenRandom@ on Windows), then run the given action. +-- $seeding -- --- This is a somewhat expensive function, and is intended to be called --- only occasionally (e.g. once per thread). You should use the `Gen` --- it creates to generate many random numbers. -withSystemRandom :: PrimBase m - => (Gen (PrimState m) -> m a) -> IO a -withSystemRandom act = do - seed <- acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do - seen <- atomicModifyIORef warned ((,) True) +-- Library provides several functions allowing to intialize generator +-- using OS-provided randomness: \"@\/dev\/urandom@\" on Unix-like +-- systems or @RtlGenRandom@ on Windows. This is a somewhat expensive +-- function, and is intended to be called only occasionally (e.g. once +-- per thread). You should use the `Gen` it creates to generate many +-- random numbers. + +createSystemRandomList :: IO [Word32] +createSystemRandomList = do + acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do + seen <- atomicModifyIORef seedCreatetionWarned ((,) True) unless seen $ E.handle (\(_::E.IOException) -> return ()) $ do hPutStrLn stderr $ "Warning: Couldn't use randomness source " ++ randomSourceName hPutStrLn stderr ("Warning: using system clock for seed instead " ++ "(quality will be lower)") acquireSeedTime - unsafePrimToIO $ initialize (I.fromList seed) >>= act - where - warned = unsafePerformIO $ newIORef False - {-# NOINLINE warned #-} --- | Seed a PRNG with data from the system's fast source of pseudo-random --- numbers. All the caveats of 'withSystemRandom' apply here as well. +seedCreatetionWarned :: IORef Bool +seedCreatetionWarned = unsafePerformIO $ newIORef False +{-# NOINLINE seedCreatetionWarned #-} + + + +-- | Generate random seed for generator using system's fast source of +-- pseudo-random numbers. +-- +-- @since 0.15.0.0 +createSystemSeed :: IO Seed +createSystemSeed = do + seed <- createSystemRandomList + return $! toSeed $ I.fromList seed + +-- | Seed a PRNG with data from the system's fast source of +-- pseudo-random numbers. createSystemRandom :: IO GenIO -createSystemRandom = withSystemRandom (return :: GenIO -> IO GenIO) +createSystemRandom = initialize . I.fromList =<< createSystemRandomList + +-- | Seed PRNG with data from the system's fast source of +-- pseudo-random numbers and execute computation in ST monad. +-- +-- @since 0.15.0.0 +withSystemRandomST :: (forall s. Gen s -> ST s a) -> IO a +withSystemRandomST act = do + seed <- createSystemSeed + return $! runST $ act =<< restore seed + +-- | Seed a PRNG with data from the system's fast source of +-- pseudo-random numbers, then run the given action. +-- +-- This function is unsafe and for example allows STRefs or any +-- other mutable data structure to escape scope: +-- +-- >>> ref <- withSystemRandom $ \_ -> newSTRef 1 +-- >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref +-- 2 +-- >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref +-- 3 +withSystemRandom :: PrimBase m + => (Gen (PrimState m) -> m a) -> IO a +withSystemRandom act = do + seed <- createSystemSeed + unsafePrimToIO $ act =<< restore seed +{-# DEPRECATED withSystemRandom "Use withSystemRandomST or createSystemSeed or createSystemRandom instead" #-} + + -- | Compute the next index into the state pool. This is simply -- addition modulo 256. nextIndex :: Integral a => a -> Int @@ -453,37 +614,19 @@ {-# INLINE aa #-} uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32 +-- NOTE [Carry value] uniformWord32 (Gen q) = do - i <- nextIndex `liftM` M.unsafeRead q ioff - c <- fromIntegral `liftM` M.unsafeRead q coff - qi <- fromIntegral `liftM` M.unsafeRead q i + i <- nextIndex <$> M.unsafeRead q ioff + c <- fromIntegral <$> M.unsafeRead q coff + qi <- fromIntegral <$> M.unsafeRead q i let t = aa * qi + c - -- The comments in this function are a proof that: - -- "if the carry value is strictly smaller than the multiplicator, - -- the next carry value is also strictly smaller than the multiplicator." - -- Eventhough the proof is written in terms of the actual value of the multiplicator, - -- it holds for any multiplicator value /not/ greater than maxBound 'Word32' - -- - -- (In the code, the multiplicator is aa, the carry value is c, - -- the next carry value is c''.) - -- - -- So we'll assume that c < aa, and show that c'' < aa : - -- - -- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32') - -- hence aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2. - -- - -- hence t < 0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2) - -- hence t <= 0x5BCF5AB1FFFFFFFF c' = fromIntegral (t `shiftR` 32) - -- c' < 0x5BCF5AB1 x = fromIntegral t + c' (# x', c'' #) | x < c' = (# x + 1, c' + 1 #) | otherwise = (# x, c' #) - -- hence c'' < 0x5BCF5AB2, - -- hence c'' < aa, which is what we wanted to prove. M.unsafeWrite q i x' M.unsafeWrite q ioff (fromIntegral i) - M.unsafeWrite q coff (fromIntegral c'') + M.unsafeWrite q coff c'' return x' {-# INLINE uniformWord32 #-} @@ -495,11 +638,11 @@ uniform2 :: PrimMonad m => (Word32 -> Word32 -> a) -> Gen (PrimState m) -> m a uniform2 f (Gen q) = do - i <- nextIndex `liftM` M.unsafeRead q ioff + i <- nextIndex <$> M.unsafeRead q ioff let j = nextIndex i - c <- fromIntegral `liftM` M.unsafeRead q coff - qi <- fromIntegral `liftM` M.unsafeRead q i - qj <- fromIntegral `liftM` M.unsafeRead q j + c <- fromIntegral <$> M.unsafeRead q coff + qi <- fromIntegral <$> M.unsafeRead q i + qj <- fromIntegral <$> M.unsafeRead q j let t = aa * qi + c c' = fromIntegral (t `shiftR` 32) x = fromIntegral t + c' @@ -513,14 +656,14 @@ M.unsafeWrite q i x' M.unsafeWrite q j y' M.unsafeWrite q ioff (fromIntegral j) - M.unsafeWrite q coff (fromIntegral d'') + M.unsafeWrite q coff d'' return $! f x' y' {-# INLINE uniform2 #-} -- Type family for fixed size integrals. For signed data types it's --- its unsigned couterpart with same size and for unsigned data types +-- its unsigned counterpart with same size and for unsigned data types -- it's same type -type family Unsigned a :: * +type family Unsigned a :: Type type instance Unsigned Int8 = Word8 type instance Unsigned Int16 = Word16 @@ -532,24 +675,8 @@ type instance Unsigned Word32 = Word32 type instance Unsigned Word64 = Word64 --- This is workaround for bug #25. --- --- GHC-7.6 has a bug (#8072) which results in calculation of wrong --- number of buckets in function `uniformRange'. Consequently uniformR --- generates values in wrong range. --- --- Bug only affects 32-bit systems and Int/Word data types. Word32 --- works just fine. So we set Word32 as unsigned counterpart for Int --- and Word on 32-bit systems. It's done only for GHC-7.6 because --- other versions are unaffected by the bug and we expect that GHC may --- optimise code which uses Word better. -#if (WORD_SIZE_IN_BITS < 64) && (__GLASGOW_HASKELL__ == 706) -type instance Unsigned Int = Word32 -type instance Unsigned Word = Word32 -#else type instance Unsigned Int = Word type instance Unsigned Word = Word -#endif -- Subtract two numbers under assumption that x>=y and store result in @@ -591,12 +718,25 @@ -- | Generate a vector of pseudo-random variates. This is not -- necessarily faster than invoking 'uniform' repeatedly in a loop, -- but it may be more convenient to use in some situations. -uniformVector :: (PrimMonad m, Variate a, Vector v a) - => Gen (PrimState m) -> Int -> m (v a) -uniformVector gen n = G.replicateM n (uniform gen) +uniformVector + :: (PrimMonad m, Random.StatefulGen g m, Random.Uniform a, Vector v a) + => g -> Int -> m (v a) +-- NOTE: We use in-place mutation in order to generate vector instead +-- of generateM because latter will go though intermediate list until +-- we're working in IO/ST monad +-- +-- See: https://github.com/haskell/vector/issues/208 for details +uniformVector gen n = do + mu <- GM.unsafeNew n + let go !i | i < n = Random.uniformM gen >>= GM.unsafeWrite mu i >> go (i+1) + | otherwise = G.unsafeFreeze mu + go 0 {-# INLINE uniformVector #-} +-- This is default seed for the generator and used when no seed is +-- specified or seed is only partial. It's not known how it was +-- generated but it looks random enough defaultSeed :: I.Vector Word32 defaultSeed = I.fromList [ 0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184, @@ -695,3 +835,40 @@ -- -- This is almost as compact as the original code that the compiler -- rejected. + + + +-- $setup +-- +-- >>> import Control.Monad +-- >>> import Data.Word +-- >>> import Data.STRef +-- >>> :set -Wno-deprecations + + +-- NOTE [Carry value] +-- ------------------ +-- This is proof of statement: +-- +-- > if the carry value is strictly smaller than the multiplicator, +-- > the next carry value is also strictly smaller than the multiplicator. +-- +-- Even though the proof is written in terms of the actual value of the +-- multiplicator, it holds for any multiplicator value /not/ greater +-- than maxBound 'Word32' +-- +-- (In the code, the multiplicator is aa, the carry value is c, +-- the next carry value is c''.) +-- +-- So we'll assume that c < aa, and show that c'' < aa : +-- +-- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32') +-- +-- Then we get following: +-- +-- aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2. +-- t < 0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2) +-- t <= 0x5BCF5AB1FFFFFFFF +-- c' < 0x5BCF5AB1 +-- c'' < 0x5BCF5AB2, +-- c'' < aa, which is what we wanted to prove. diff --git a/System/Random/MWC/CondensedTable.hs b/System/Random/MWC/CondensedTable.hs --- a/System/Random/MWC/CondensedTable.hs +++ b/System/Random/MWC/CondensedTable.hs @@ -29,7 +29,6 @@ ) where import Control.Arrow (second,(***)) -import Control.Monad.Primitive (PrimMonad(..)) import Data.Word import Data.Int @@ -41,13 +40,12 @@ import qualified Data.Vector as V import Data.Vector.Generic (Vector) import Numeric.SpecFunctions (logFactorial) +import System.Random.Stateful import Prelude hiding ((++)) -import System.Random.MWC - -- | A lookup table for arbitrary discrete distributions. It allows -- the generation of random variates in /O(1)/. Note that probability -- is quantized in units of @1/2^32@, and all distributions with @@ -77,12 +75,11 @@ -- | Generate a random value using a condensed table. -genFromTable :: (PrimMonad m, Vector v a) => - CondensedTable v a -> Gen (PrimState m) -> m a +genFromTable :: (StatefulGen g m, Vector v a) => CondensedTable v a -> g -> m a {-# INLINE genFromTable #-} genFromTable table gen = do - w <- uniform gen - return $ lookupTable table $ fromIntegral (w :: Word32) + w <- uniformM gen + return $! lookupTable table $ fromIntegral (w :: Word32) lookupTable :: Vector v a => CondensedTable v a -> Word64 -> a {-# INLINE lookupTable #-} @@ -101,7 +98,7 @@ -- | Generate a condensed lookup table from a list of outcomes with -- given probabilities. The vector should be non-empty and the --- probabilites should be non-negative and sum to 1. If this is not +-- probabilities should be non-negative and sum to 1. If this is not -- the case, this algorithm will construct a table for some -- distribution that may bear no resemblance to what you intended. tableFromProbabilities @@ -151,7 +148,7 @@ {-# INLINE tableFromIntWeights #-} tableFromIntWeights v | n == 0 = pkgError "tableFromIntWeights" "empty table" - -- Single element tables should be treated sepately. Otherwise + -- Single element tables should be treated separately. Otherwise -- they will confuse correctWeights | n == 1 = let m = 2^(32::Int) - 1 -- Works for both Word32 & Word64 in CondensedTable @@ -224,8 +221,8 @@ return arr --- | Create a lookup table for the Poisson distibution. Note that --- table construction may have significant cost. For λ < 100 it +-- | Create a lookup table for the Poisson distribution. Note that +-- table construction may have significant cost. For λ < 100 it -- takes as much time to build table as generation of 1000-30000 -- variates. tablePoisson :: Double -> CondensedTableU Int diff --git a/System/Random/MWC/Distributions.hs b/System/Random/MWC/Distributions.hs --- a/System/Random/MWC/Distributions.hs +++ b/System/Random/MWC/Distributions.hs @@ -1,4 +1,5 @@ -{-# LANGUAGE BangPatterns, CPP, GADTs, FlexibleContexts, ScopedTypeVariables #-} +{-# LANGUAGE MultiWayIf #-} +{-# LANGUAGE BangPatterns, GADTs, FlexibleContexts, ScopedTypeVariables #-} -- | -- Module : System.Random.MWC.Distributions -- Copyright : (c) 2012 Bryan O'Sullivan @@ -27,6 +28,8 @@ , geometric0 , geometric1 , bernoulli + , binomial + , poisson -- ** Multivariate , dirichlet -- * Permutations @@ -38,19 +41,16 @@ ) where import Prelude hiding (mapM) -import Control.Monad (liftM) import Control.Monad.Primitive (PrimMonad, PrimState) import Data.Bits ((.&.)) import Data.Foldable (foldl') -#if !MIN_VERSION_base(4,8,0) -import Data.Traversable (Traversable) -#endif import Data.Traversable (mapM) import Data.Word (Word32) -import System.Random.MWC (Gen, uniform, uniformR) +import System.Random.Stateful (StatefulGen(..),Uniform(..),UniformRange(..),uniformDoublePositive01M, uniformDouble01M) import qualified Data.Vector.Unboxed as I import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as M +import Numeric.SpecFunctions (logFactorial) -- Unboxed 2-tuple data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double @@ -58,10 +58,10 @@ -- | Generate a normally distributed random variate with given mean -- and standard deviation. -normal :: PrimMonad m +normal :: StatefulGen g m => Double -- ^ Mean -> Double -- ^ Standard deviation - -> Gen (PrimState m) + -> g -> m Double {-# INLINE normal #-} normal m s gen = do @@ -75,13 +75,13 @@ -- Compared to the ziggurat algorithm usually used, this is slower, -- but generates more independent variates that pass stringent tests -- of randomness. -standard :: PrimMonad m => Gen (PrimState m) -> m Double +standard :: StatefulGen g m => g -> m Double {-# INLINE standard #-} standard gen = loop where loop = do - u <- (subtract 1 . (*2)) `liftM` uniform gen - ri <- uniform gen + u <- subtract 1 . (*2) <$> uniformDoublePositive01M gen + ri <- uniformM gen let i = fromIntegral ((ri :: Word32) .&. 127) bi = I.unsafeIndex blocks i bj = I.unsafeIndex blocks (i+1) @@ -93,14 +93,14 @@ xx = x * x d = exp (-0.5 * (bi * bi - xx)) e = exp (-0.5 * (bj * bj - xx)) - c <- uniform gen + c <- uniformDoublePositive01M gen if e + c * (d - e) < 1 then return x else loop normalTail neg = tailing where tailing = do - x <- ((/rNorm) . log) `liftM` uniform gen - y <- log `liftM` uniform gen + x <- (/ rNorm) . log <$> uniformDoublePositive01M gen + y <- log <$> uniformDoublePositive01M gen if y * (-2) < x * x then tailing else return $! if neg then x - rNorm else rNorm - x @@ -129,36 +129,36 @@ -- | Generate an exponentially distributed random variate. -exponential :: PrimMonad m +exponential :: StatefulGen g m => Double -- ^ Scale parameter - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Double {-# INLINE exponential #-} exponential b gen = do - x <- uniform gen + x <- uniformDoublePositive01M gen return $! - log x / b -- | Generate truncated exponentially distributed random variate. -truncatedExp :: PrimMonad m +truncatedExp :: StatefulGen g m => Double -- ^ Scale parameter -> (Double,Double) -- ^ Range to which distribution is -- truncated. Values may be negative. - -> Gen (PrimState m) -- ^ Generator. + -> g -- ^ Generator. -> m Double {-# INLINE truncatedExp #-} truncatedExp scale (a,b) gen = do -- We shift a to 0 and then generate distribution truncated to [0,b-a] -- It's easier let delta = b - a - p <- uniform gen + p <- uniformDoublePositive01M gen return $! a - log ( (1 - p) + p*exp(-scale*delta)) / scale -- | Random variate generator for gamma distribution. -gamma :: PrimMonad m +gamma :: (StatefulGen g m) => Double -- ^ Shape parameter -> Double -- ^ Scale parameter - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Double {-# INLINE gamma #-} gamma a b gen @@ -167,13 +167,13 @@ where mainloop = do T x v <- innerloop - u <- uniform gen + u <- uniformDoublePositive01M gen let cont = u > 1 - 0.331 * sqr (sqr x) && log u > 0.5 * sqr x + a1 * (1 - v + log v) -- Rarely evaluated case () of _| cont -> mainloop | a >= 1 -> return $! a1 * v * b - | otherwise -> do y <- uniform gen + | otherwise -> do y <- uniformDoublePositive01M gen return $! y ** (1 / a) * a1 * v * b -- inner loop innerloop = do @@ -188,9 +188,9 @@ -- | Random variate generator for the chi square distribution. -chiSquare :: PrimMonad m +chiSquare :: StatefulGen g m => Int -- ^ Number of degrees of freedom - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Double {-# INLINE chiSquare #-} chiSquare n gen @@ -201,14 +201,14 @@ -- | Random variate generator for the geometric distribution, -- computing the number of failures before success. Distribution's -- support is [0..]. -geometric0 :: PrimMonad m +geometric0 :: StatefulGen g m => Double -- ^ /p/ success probability lies in (0,1] - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Int {-# INLINE geometric0 #-} geometric0 p gen | p == 1 = return 0 - | p > 0 && p < 1 = do q <- uniform gen + | p > 0 && p < 1 = do q <- uniformDoublePositive01M gen -- FIXME: We want to use log1p here but it will -- introduce dependency on math-functions. return $! floor $ log q / log (1 - p) @@ -217,19 +217,19 @@ -- | Random variate generator for geometric distribution for number of -- trials. Distribution's support is [1..] (i.e. just 'geometric0' -- shifted by 1). -geometric1 :: PrimMonad m +geometric1 :: StatefulGen g m => Double -- ^ /p/ success probability lies in (0,1] - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Int {-# INLINE geometric1 #-} geometric1 p gen = do n <- geometric0 p gen return $! n + 1 -- | Random variate generator for Beta distribution -beta :: PrimMonad m +beta :: StatefulGen g m => Double -- ^ alpha (>0) -> Double -- ^ beta (>0) - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Double {-# INLINE beta #-} beta a b gen = do @@ -238,9 +238,9 @@ return $! x / (x+y) -- | Random variate generator for Dirichlet distribution -dirichlet :: (PrimMonad m, Traversable t) +dirichlet :: (StatefulGen g m, Traversable t) => t Double -- ^ container of parameters - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m (t Double) {-# INLINE dirichlet #-} dirichlet t gen = do @@ -249,12 +249,12 @@ return $ fmap (/total) t' -- | Random variate generator for Bernoulli distribution -bernoulli :: PrimMonad m +bernoulli :: StatefulGen g m => Double -- ^ Probability of success (returning True) - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Bool {-# INLINE bernoulli #-} -bernoulli p gen = ( uniformDoublePositive01M gen -- | Random variate generator for categorical distribution. -- @@ -262,16 +262,16 @@ -- "System.Random.MWC.CondensedTable" will offer better -- performance. If only few is needed this function will faster -- since it avoids costs of setting up table. -categorical :: (PrimMonad m, G.Vector v Double) +categorical :: (StatefulGen g m, G.Vector v Double) => v Double -- ^ List of weights [>0] - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Int {-# INLINE categorical #-} categorical v gen | G.null v = pkgError "categorical" "empty weights!" | otherwise = do let cv = G.scanl1' (+) v - p <- (G.last cv *) `liftM` uniform gen + p <- (G.last cv *) <$> uniformDoublePositive01M gen return $! case G.findIndex (>=p) cv of Just i -> i Nothing -> pkgError "categorical" "bad weights!" @@ -279,9 +279,9 @@ -- | Random variate generator for categorical distribution where the -- weights are in the log domain. It's implemented in terms of -- 'categorical'. -logCategorical :: (PrimMonad m, G.Vector v Double) +logCategorical :: (StatefulGen g m, G.Vector v Double) => v Double -- ^ List of logarithms of weights - -> Gen (PrimState m) -- ^ Generator + -> g -- ^ Generator -> m Int {-# INLINE logCategorical #-} logCategorical v gen @@ -294,9 +294,9 @@ -- It returns random permutation of vector /[0 .. n-1]/. -- -- This is the Fisher-Yates shuffle -uniformPermutation :: forall m v. (PrimMonad m, G.Vector v Int) +uniformPermutation :: forall g m v. (StatefulGen g m, PrimMonad m, G.Vector v Int) => Int - -> Gen (PrimState m) + -> g -> m (v Int) {-# INLINE uniformPermutation #-} uniformPermutation n gen @@ -306,9 +306,9 @@ -- | Random variate generator for a uniformly distributed shuffle (all -- shuffles are equiprobable) of a vector. It uses Fisher-Yates -- shuffle algorithm. -uniformShuffle :: (PrimMonad m, G.Vector v a) +uniformShuffle :: (StatefulGen g m, PrimMonad m, G.Vector v a) => v a - -> Gen (PrimState m) + -> g -> m (v a) {-# INLINE uniformShuffle #-} uniformShuffle vec gen @@ -320,9 +320,9 @@ -- | In-place uniformly distributed shuffle (all shuffles are -- equiprobable)of a vector. -uniformShuffleM :: (PrimMonad m, M.MVector v a) +uniformShuffleM :: (StatefulGen g m, PrimMonad m, M.MVector v a) => v (PrimState m) a - -> Gen (PrimState m) + -> g -> m () {-# INLINE uniformShuffleM #-} uniformShuffleM vec gen @@ -332,7 +332,7 @@ n = M.length vec lst = n-1 loop i | i == lst = return () - | otherwise = do j <- uniformR (i,lst) gen + | otherwise = do j <- uniformRM (i,lst) gen M.unsafeSwap vec i j loop (i+1) @@ -345,8 +345,223 @@ pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func ++ ": " ++ msg +-- | Random variate generator for Binomial distribution. Will throw +-- exception when parameters are out range. +-- +-- The probability of getting exactly k successes in n trials is +-- given by the probability mass function: +-- +-- \[ +-- f(k;n,p) = \Pr(X = k) = \binom n k p^k(1-p)^{n-k} +-- \] +-- +-- @since 0.15.1.0 +binomial :: forall g m . StatefulGen g m + => Int -- ^ Number of trials, must be positive. + -> Double -- ^ Probability of success \(p \in [0,1]\) + -> g -- ^ Generator + -> m Int +{-# INLINE binomial #-} +binomial n p gen + | n <= 0 = pkgError "binomial" "number of trials must be positive" + | p < 0.0 || p > 1.0 = pkgError "binomial" "probability must be >= 0 and <= 1" + | p == 0.0 = return 0 + | p == 1.0 = return n + | p <= 0.5 = if + | fromIntegral n * p < inv_thr -> binomialInv n p gen + | otherwise -> binomialTPE n p gen + | p > 0.5 = do + ix <- case 1 - p of + p' | fromIntegral n * p' < inv_thr -> binomialInv n p' gen + | otherwise -> binomialTPE n p' gen + pure $! n - ix + -- Reachable when p is NaN + | otherwise = pkgError "binomial" "probability must be >= 0 and <= 1" + where + -- Threshold for preferring the BINV algorithm / inverse cdf + -- logic. The paper suggests 10, Ranlib uses 30, R uses 30, Rust uses + -- 10 and GSL uses 14. + inv_thr = 10 +-- Binomial-Triangle-Parallelogram-Exponential algorithm (BTPE) +-- described in Kachitvichyanukul1988 +binomialTPE :: forall g m . StatefulGen g m => Int -> Double -> g -> m Int +{-# INLINE binomialTPE #-} +binomialTPE n p g = loop + where + -- Main accept/reject loop + loop = do + u <- uniformRM (0.0, p4) g + v <- uniformDoublePositive01M g + selectArea u v + -- Acceptance / rejection of sample [step 5] + acceptReject :: Int -> Double -> m Int + acceptReject !ix !v + | var <= accept = return ix + | otherwise = loop + where + var = log v + accept = logFactorial bigM + logFactorial (n - bigM) - + logFactorial ix - logFactorial (n - ix) + + fromIntegral (ix - bigM) * log (p / q) + -- Select area to be used [Steps 1-4] + selectArea :: Double -> Double -> m Int + selectArea !u !v + -- Triangular region + | u <= p1 = return $! floor $ xm - p1 * v + u + -- Parallelogram region + | u <= p2 = do let x = xl + (u - p1) / c + w = v * c + 1.0 - abs (x - xm) / p1 + if w > 1 || w <= 0 + then loop + else do let ix = floor x + acceptReject ix w + -- Left tail + | u <= p3 = case floor $ xl + log v / lambdaL of + ix | ix < 0 -> loop + | otherwise -> do let w = v * (u - p2) * lambdaL + acceptReject ix w + -- Right tail + | otherwise = case floor $ xr - log v / lambdaR of + ix | ix > n -> loop + | otherwise -> do let w = v * (u - p3) * lambdaR + acceptReject ix w + ---------------------------------------- + -- Constants used in algorithm. See [Step 0] + q = 1 - p + np = fromIntegral n * p + ffm = np + p + bigM = floor ffm + -- Half integer mean (tip of triangle) + xm = fromIntegral bigM + 0.5 + -- p1: the distance to the left and right edges of the triangle + -- region below the target distribution; since height=1, also: + -- area of region (half base * height) + !p1 = let npq = np * q + in fromIntegral (floor (2.195 * sqrt npq - 4.6 * q) :: Int) + 0.5 + xl = xm - p1 -- Left edge of triangle + xr = xm + p1 -- Right edge of triangle + c = 0.134 + 20.5 / (15.3 + fromIntegral bigM) + -- p1 + area of parallelogram region + !p2 = p1 * (1.0 + c + c) + -- + lambdaL = let al = (ffm - xl) / (ffm - xl * p) + in al * (1.0 + 0.5 * al) + lambdaR = let ar = (xr - ffm) / (xr * q) + in ar * (1.0 + 0.5 * ar) + -- p2 + area of left tail + !p3 = p2 + c / lambdaL + -- p3 + area of right tail + !p4 = p3 + c / lambdaR + +-- Compute binomial variate using inversion method (BINV in +-- Kachitvichyanukul1988) +binomialInv :: StatefulGen g m => Int -> Double -> g -> m Int +{-# INLINE binomialInv #-} +binomialInv n p g = do + u <- uniformDoublePositive01M g + return $! invertBinomial n p u + +-- This function is defined on top level to avoid inlining it since it's rather +-- large and we don't need specializations since it's monomorphic anyway +invertBinomial + :: Int -- N of trials + -> Double -- probability of success + -> Double -- Output of PRNG + -> Int +invertBinomial !n !p !u0 = invert (q^n) u0 0 + where + -- We forcing s&a in order to avoid allocating thunks. Those are + -- more expensive than computing them unconditionally + q = 1 - p + !s = p / q + !a = fromIntegral (n + 1) * s + -- + invert !r !u !x + | u <= r = x + | otherwise = invert r' u' x' + where + u' = u - r + x' = x + 1 + r' = r * ((a / fromIntegral x') - s) + + +-- | Random variate generate for Poisson distribution. +-- +-- If parameter λ is within 10 σ or greater than @maxBound :: Int@ +-- error is raised since result may not be representable as @Int@ +-- +-- @since 0.15.3.0 +poisson + :: StatefulGen g m + => Double -- ^ Rate parameter, also known as \( \lambda \) + -> g -- ^ Generator + -> m Int +{-# INLINE poisson #-} +poisson lambda gen + | lambda > maxPoissonLam + = pkgError "poisson" + $ "lambda is too large: "++show lambda++" Result won't fit into Int" + | lambda > 10 = poissonAtkinson lambda gen + | lambda >= 0 = poissonInterArrival lambda gen + | otherwise + = pkgError "poisson" "Lambda parameter must be greater than zero" + +maxPoissonLam :: Double +maxPoissonLam = m - 10 * sqrt m where + m = fromIntegral (maxBound :: Int) + + +-- This uses the fact that if N(t) is a Poisson process +-- with rate lambda, then the counting process N(t) can +-- be represented as interarrival times X[1], X[2],... with +-- X[i] ~ Exp(lambda). +poissonInterArrival :: StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonInterArrival #-} +poissonInterArrival lambda gen = do + loop 0 1.0 + where + loop !k !p = do + p' <- (*p) <$> uniformDouble01M gen + if p' <= bigL then return $! k else loop (k+1) p' + bigL = exp (negate lambda) + +-- Attributed to Atkinson, via Casella. Uses a rejection +-- method that uses logistic distribution as the envelope +-- distribution. +poissonAtkinson :: forall g m . StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonAtkinson #-} +poissonAtkinson lambda gen = loop + where loop :: m Int + loop = do + bigU <- uniformDouble01M gen + let x = (alpha - log ((1.0 - bigU) / bigU)) / bbeta + if x < (-0.5) + then loop + else do + bigV <- uniformDouble01M gen + let n = floor (x + 0.5) :: Int + y = alpha - bbeta * x + logFacN = logFactorial n + lhs = y + log (bigV / (1.0 + exp y)**2) + rhs = bigK + fromIntegral n * logLambda - logFacN + if lhs <= rhs + then return n + else loop + bigC,alpha,bbeta,bigK,logLambda :: Double + bigC = 0.767 - 3.36 / lambda + bbeta = pi / sqrt (3.0 * lambda) + alpha = bbeta * lambda + bigK = log bigC - lambda - log bbeta + logLambda = log lambda + -- $references -- -- * Doornik, J.A. (2005) An improved ziggurat method to generate @@ -357,3 +572,18 @@ -- (2007). Gaussian random number generators. -- /ACM Computing Surveys/ 39(4). -- +-- +-- * Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random +-- Variate Generation. Communications of the ACM, 31, 2 (February, +-- 1988) 216. +-- Here's an example of how the algorithm's sampling regions look +-- ![Something](docs/RecreateFigure.svg) +-- +-- * Devroye, L. (1986) Non-uniform Random Variate Generation. +-- Springer Verlag. Chapter 10: Discrete Univariate Distributions. +-- +-- +-- * Robert, C.P. & Casella, G. Monte Carlo Statistical Methods. +-- Springer Texts in Statistics. Algorithm A6: Atkinson's Method +-- for Generating Poisson Random Variables. +-- diff --git a/System/Random/MWC/SeedSource.hs b/System/Random/MWC/SeedSource.hs --- a/System/Random/MWC/SeedSource.hs +++ b/System/Random/MWC/SeedSource.hs @@ -10,7 +10,6 @@ , randomSourceName ) where -import Control.Monad (liftM) import Data.Word (Word32,Word64) import Data.Bits (shiftR) import Data.Ratio ((%), numerator) @@ -22,16 +21,17 @@ #if defined(mingw32_HOST_OS) import Foreign.Ptr import Foreign.C.Types +#else +import System.IO (IOMode(..), hGetBuf, withBinaryFile) #endif import System.CPUTime (cpuTimePrecision, getCPUTime) -import System.IO (IOMode(..), hGetBuf, withBinaryFile) -- Acquire seed from current time. This is horrible fallback for -- Windows system. acquireSeedTime :: IO [Word32] acquireSeedTime = do - c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime - t <- toRational `liftM` getPOSIXTime + c <- numerator . (% cpuTimePrecision) <$> getCPUTime + t <- toRational <$> getPOSIXTime let n = fromIntegral (numerator t) :: Word64 return [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)] diff --git a/bench-papi/Bench.hs b/bench-papi/Bench.hs new file mode 100644 --- /dev/null +++ b/bench-papi/Bench.hs @@ -0,0 +1,14 @@ +-- | +-- Here we reexport definitions of tasty-bench +module Bench + ( whnf + , nf + , nfIO + , whnfIO + , bench + , bgroup + , defaultMain + , benchIngredients + ) where + +import Test.Tasty.PAPI diff --git a/bench-time/Bench.hs b/bench-time/Bench.hs new file mode 100644 --- /dev/null +++ b/bench-time/Bench.hs @@ -0,0 +1,14 @@ +-- | +-- Here we reexport definitions of tasty-bench +module Bench + ( whnf + , nf + , nfIO + , whnfIO + , bench + , bgroup + , defaultMain + , benchIngredients + ) where + +import Test.Tasty.Bench diff --git a/bench/Benchmark.hs b/bench/Benchmark.hs new file mode 100644 --- /dev/null +++ b/bench/Benchmark.hs @@ -0,0 +1,161 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +module Main(main) where + +import Control.Exception +import Data.Int +import Data.Word +import Data.Proxy +import qualified Data.Vector.Unboxed as U +import qualified System.Random as R +import System.Random.Stateful (StatefulGen) +import System.Random.MWC +import System.Random.MWC.Distributions +import System.Random.MWC.CondensedTable +import qualified System.Random.Mersenne as M + +import Test.Tasty.Options +import Test.Tasty.Runners +import Test.Tasty (includingOptions) +import Bench + +-- | Size of vector used in benchmarks +newtype Iterations = Iterations Int + +instance IsOption Iterations where + defaultValue = Iterations 10000 + parseValue = fmap Iterations . safeRead + optionName = pure "iter" + optionHelp = pure "Number of iteration in sampling benchmarks" + + +loop :: Iterations -> IO a -> IO () +{-# INLINE loop #-} +loop (Iterations n) act = go n where + go i | i <= 0 = pure () + | otherwise = do _ <- evaluate =<< act + go (i - 1) + +makeTableUniform :: Int -> CondensedTable U.Vector Int +makeTableUniform n = + tableFromProbabilities $ U.zip (U.enumFromN 0 n) (U.replicate n (1 / fromIntegral n)) +{-# INLINE makeTableUniform #-} + +main :: IO () +main = do + -- Set up tasty + let tasty_opts = [Option (Proxy :: Proxy Iterations)] + ingredients = includingOptions tasty_opts : benchIngredients + opts <- parseOptions ingredients (bench "Fake" (nf id ())) + let iter = lookupOption opts + -- Set up RNG + mwc <- create + seed <- save mwc + mtg <- M.newMTGen . Just =<< uniform mwc + defaultMainWithIngredients ingredients $ bgroup "All" + [ bgroup "mwc" + -- One letter group names are used so they will fit on the plot. + -- + -- U - uniform + -- R - uniformR + -- D - distribution + [ bgroup "U" + [ bench "Double" $ whnfIO $ loop iter (uniform mwc :: IO Double) + , bench "Int" $ whnfIO $ loop iter (uniform mwc :: IO Int) + , bench "Int8" $ whnfIO $ loop iter (uniform mwc :: IO Int8) + , bench "Int16" $ whnfIO $ loop iter (uniform mwc :: IO Int16) + , bench "Int32" $ whnfIO $ loop iter (uniform mwc :: IO Int32) + , bench "Int64" $ whnfIO $ loop iter (uniform mwc :: IO Int64) + , bench "Word" $ whnfIO $ loop iter (uniform mwc :: IO Word) + , bench "Word8" $ whnfIO $ loop iter (uniform mwc :: IO Word8) + , bench "Word16" $ whnfIO $ loop iter (uniform mwc :: IO Word16) + , bench "Word32" $ whnfIO $ loop iter (uniform mwc :: IO Word32) + , bench "Word64" $ whnfIO $ loop iter (uniform mwc :: IO Word64) + ] + , bgroup "R" + -- I'm not entirely convinced that this is right way to test + -- uniformR. /A.Khudyakov/ + [ bench "Double" $ whnfIO $ loop iter (uniformR (-3.21,26) mwc :: IO Double) + , bench "Int" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int) + , bench "Int8" $ whnfIO $ loop iter (uniformR (-12,4) mwc :: IO Int8) + , bench "Int16" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int16) + , bench "Int32" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int32) + , bench "Int64" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int64) + , bench "Word" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word) + , bench "Word8" $ whnfIO $ loop iter (uniformR (34,63) mwc :: IO Word8) + , bench "Word16" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word16) + , bench "Word32" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word32) + , bench "Word64" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word64) + ] + , bgroup "D" + [ bench "standard" $ whnfIO $ loop iter (standard mwc :: IO Double) + , bench "normal" $ whnfIO $ loop iter (normal 1 3 mwc :: IO Double) + , bench "exponential" $ whnfIO $ loop iter (exponential 3 mwc :: IO Double) + , bench "gamma,a<1" $ whnfIO $ loop iter (gamma 0.5 1 mwc :: IO Double) + , bench "gamma,a>1" $ whnfIO $ loop iter (gamma 2 1 mwc :: IO Double) + , bench "beta" $ whnfIO $ loop iter (beta 2 3 mwc :: IO Double) + , bench "chiSquare" $ whnfIO $ loop iter (chiSquare 4 mwc :: IO Double) + -- NOTE: We switch between algorithms when Np=10 + , bgroup "binomial" + [ bench (show p ++ " " ++ show n) $ whnfIO $ loop iter (binomial n p mwc) + | (n,p) <- [ (2, 0.2), (2, 0.5), (2, 0.8) + , (10, 0.1), (10, 0.9) + , (20, 0.2), (20, 0.8) + -- + , (60, 0.2), (60, 0.8) + , (600, 0.2), (600, 0.8) + , (6000, 0.2), (6000, 0.8) + ] + ] + ] + , bgroup "poisson" + [ bench (show lam) $ whnfIO $ loop iter (poisson lam mwc) + | lam <- [0.1, 1, 8, 12, 100, 1000, 10000] + ] + -- Test sampling performance. Table creation must be floated out! + , bgroup "CT/gen" $ concat + [ [ bench ("uniform "++show i) $ whnfIO $ loop iter (genFromTable tbl mwc) + | i <- [2..10] + , let tbl = makeTableUniform i + ] + , [ bench ("poisson " ++ show l) $ whnfIO $ loop iter (genFromTable tbl mwc) + | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000] + , let tbl = tablePoisson l + ] + , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnfIO $ loop iter (genFromTable tbl mwc) + | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)] + , let tbl = tableBinomial n p + ] + ] + -- Benchmarking of setting up table (no need to use iterations + -- here!). Setting up is rather expensive + , bgroup "CT/table" $ concat + [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i + | i <- [2..30] + ] + , [ bench ("poisson " ++ show l) $ whnf tablePoisson l + | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000] + ] + , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p + | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)] + ] + ] + ] + , bgroup "random" + [ + bench "Double" $ whnfIO $ loop iter (R.randomIO :: IO Double) + , bench "Int" $ whnfIO $ loop iter (R.randomIO :: IO Int) + ] + , bgroup "mersenne" + [ + bench "Double" $ whnfIO $ loop iter (M.random mtg :: IO Double) + , bench "Int" $ whnfIO $ loop iter (M.random mtg :: IO Int) + ] +#if MIN_VERSION_random(1,3,0) + , bgroup "seed" + [ bench "SeedGen.fromSeed" $ let rseed = R.toSeed seed :: R.Seed Seed + in whnf R.fromSeed rseed + , bench "SeedGen.toSeed" $ whnf R.toSeed seed + ] +#endif + ] diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,50 @@ +## Changes in 0.15.3.0 [2025-12-30] + + * `poisson` samples for Poisson distribution is added + +## Changes in 0.15.2.0 [2025-01-18] + + * Support for `random-1.3`. + + +## Changes in 0.15.1.0 [2024-07-09] + + * Additon of binomial sampler using the rejection sampling method in + Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random + Variate Generation. Communications of the ACM, 31, 2 (February, + 1988) 216. . A more + efficient basis for e.g. the beta binomial distribution: + `beta a b g >>= \p -> binomial n p g`. + +## Changes in 0.15.0.2 [2021-08-15] + + * Doctests on 32-bit platforms are fixed. (#79) + + +## Changes in 0.15.0.1 [2020-08-08] + + * Bug in generation of Int/Word in both uniform and uniformR is fixed. (#75) + + +## Changes in 0.15.0.0 [2020-07-31] + + * `withSystemRandomST` and `createSystemSeed` are added. + + * `withSystemRandom` is deprecated. + + * `random>=1.2` is dependency of `mwc-random`. + + * Instances for type classes `StatefulGen` & `FrozenGen` defined in random-1.2 + are added for `Gen`. + + * Functions in `System.Random.MWC.Distributions` and + `System.Random.MWC.CondensedTable` now work with arbitrary `StatefulGen` + + * `System.Random.MWC.uniformVector` now works with arbitrary `StatefulGen` as + well and uses in-place initialization instead of `generateM`. It should be + faster for anything but IO and ST (those shoud remain same). + + ## Changes in 0.14.0.0 * Low level functions for acquiring random data for initialization diff --git a/docs/RecreateFigure.svg b/docs/RecreateFigure.svg new file mode 100644 --- /dev/null +++ b/docs/RecreateFigure.svg @@ -0,0 +1,3 @@ + +0.51.01.568.6231.8859.540.5FourThreeTwoOnep = 0.25n = 200Target PDFMinorizing FunctionMajorizing FunctionProbability Density (Unnormalised)Number of Successes diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,5 +1,17 @@ +cabal-version: 3.0 +build-type: Simple name: mwc-random -version: 0.14.0.0 +version: 0.15.3.0 +license: BSD-2-Clause +license-file: LICENSE +copyright: 2009, 2010, 2011 Bryan O'Sullivan + +author: Bryan O'Sullivan +maintainer: Alexey Khudyakov +homepage: https://github.com/haskell/mwc-random +bug-reports: https://github.com/haskell/mwc-random/issues + +category: Math, Statistics synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random @@ -15,38 +27,124 @@ convenient API, is faster, and supports more statistical distributions. -license: BSD3 -license-file: LICENSE -homepage: https://github.com/bos/mwc-random -bug-reports: https://github.com/bos/mwc-random/issues -author: Bryan O'Sullivan -maintainer: Bryan O'Sullivan -copyright: 2009, 2010, 2011 Bryan O'Sullivan -category: Math, Statistics -build-type: Simple -cabal-version: >= 1.8.0.4 + extra-source-files: + README.md + +extra-doc-files: + docs/*.svg changelog.md - README.markdown +tested-with: + GHC ==8.0.2 + || ==8.2.2 + || ==8.4.4 + || ==8.6.5 + || ==8.8.4 + || ==8.10.7 + || ==9.0.2 + || ==9.2.8 + || ==9.4.8 + || ==9.6.7 + || ==9.8.4 + || ==9.10.2 + || ==9.12.2 + + +source-repository head + type: git + location: https://github.com/haskell/mwc-random.git + +flag BenchPAPI + Description: Enable building of benchmarks which use instruction counters. + It requires libpapi and only works on Linux so it's protected by flag + Default: False + Manual: True + library + default-language: Haskell2010 exposed-modules: System.Random.MWC System.Random.MWC.Distributions System.Random.MWC.CondensedTable System.Random.MWC.SeedSource - build-depends: base >= 4.5 && < 5 - , primitive >= 0.6 + build-depends: base >= 4.9 && < 5 + , primitive >= 0.6.2 + , random >= 1.2 , time , vector >= 0.7 , math-functions >= 0.2.1.0 - ghc-options: -Wall -funbox-strict-fields -fwarn-tabs + ghc-options: -O2 -Wall -funbox-strict-fields -fwarn-tabs -source-repository head - type: git - location: git://github.com/bos/mwc-random +-- We want to be able to build benchmarks using both tasty-bench and tasty-papi. +-- They have similar API so we just create two shim modules which reexport +-- definitions from corresponding library and pick one in cabal file. +common bench-stanza + ghc-options: -Wall + default-language: Haskell2010 + build-depends: base < 5 + , vector >= 0.11 + , mersenne-random + , mwc-random + , random + , tasty >=1.3.1 -source-repository head - type: mercurial - location: https://bitbucket.org/bos/mwc-random +benchmark mwc-bench + import: bench-stanza + type: exitcode-stdio-1.0 + hs-source-dirs: bench bench-time + main-is: Benchmark.hs + Other-modules: Bench + build-depends: tasty-bench >= 0.3 + +benchmark mwc-bench-papi + import: bench-stanza + type: exitcode-stdio-1.0 + if impl(ghcjs) || !flag(BenchPAPI) || impl(ghc < 8.2) + buildable: False + hs-source-dirs: bench bench-papi + main-is: Benchmark.hs + Other-modules: Bench + build-depends: tasty-papi >= 0.1.2 + +test-suite mwc-prop-tests + type: exitcode-stdio-1.0 + hs-source-dirs: tests + main-is: props.hs + default-language: Haskell2010 + ghc-options: + -Wall -threaded -rtsopts + + build-depends: base + , mwc-random + , QuickCheck >=2.2 + , vector >=0.12.1 + , tasty >=1.3.1 + , tasty-quickcheck >=0.10.2 + , tasty-hunit + , random >=1.2 + , mtl + , math-functions >=0.3.4 + +test-suite mwc-doctests + type: exitcode-stdio-1.0 + main-is: doctests.hs + hs-source-dirs: tests + default-language: Haskell2010 + if impl(ghcjs) || impl(ghc < 8.0) + Buildable: False + -- Linker on macos prints warnings to console which confuses doctests. + -- We simply disable doctests on ma for older GHC + -- > warning: -single_module is obsolete + if os(darwin) && impl(ghc < 9.6) + buildable: False + build-depends: + base -any + , mwc-random -any + , doctest >=0.15 && <0.24 + -- + , bytestring + , primitive + , vector >=0.11 + , random >=1.2 diff --git a/tests/doctests.hs b/tests/doctests.hs new file mode 100644 --- /dev/null +++ b/tests/doctests.hs @@ -0,0 +1,4 @@ +import Test.DocTest (doctest) + +main :: IO () +main = doctest ["-fobject-code", "System/Random/MWC.hs"] diff --git a/tests/props.hs b/tests/props.hs new file mode 100644 --- /dev/null +++ b/tests/props.hs @@ -0,0 +1,316 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ViewPatterns #-} +import Control.Monad +import Data.Word +import Data.Proxy +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as MVU +import Numeric.SpecFunctions (logChoose,incompleteGamma,log1p,logFactorial) + +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.Tasty.Runners +import Test.Tasty.Options +import Test.Tasty.HUnit +import Test.QuickCheck.Monadic + +import System.Random.MWC +import System.Random.MWC.Distributions +import System.Random.Stateful (StatefulGen) +#if MIN_VERSION_random(1,3,0) +import qualified System.Random.Stateful as Random (SeedGen(..)) +#endif + +---------------------------------------------------------------- +-- +---------------------------------------------------------------- + +-- | Average number of events per bin for binned statistical tests +newtype NPerBin = NPerBin Int + +instance IsOption NPerBin where + defaultValue = NPerBin 100 + parseValue = fmap NPerBin . safeRead + optionName = pure "n-per-bin" + optionHelp = pure "Average number of events per bin" + + +-- | P-value for statistical test. +newtype PValue = PValue Double + +instance IsOption PValue where + defaultValue = PValue 1e-9 + parseValue = fmap PValue . safeRead + optionName = pure "pvalue" + optionHelp = pure "P-value for statistical test" + +---------------------------------------------------------------- +-- +---------------------------------------------------------------- + +main :: IO () +main = do + -- Set up tasty + let tasty_opts = [ Option (Proxy :: Proxy NPerBin) + , Option (Proxy :: Proxy PValue) + , Option (Proxy :: Proxy QuickCheckTests) + , Option (Proxy :: Proxy QuickCheckReplay) + , Option (Proxy :: Proxy QuickCheckShowReplay) + , Option (Proxy :: Proxy QuickCheckMaxSize) + , Option (Proxy :: Proxy QuickCheckMaxRatio) + , Option (Proxy :: Proxy QuickCheckVerbose) + , Option (Proxy :: Proxy QuickCheckMaxShrinks) + ] + ingredients = includingOptions tasty_opts : defaultIngredients + opts <- parseOptions ingredients (testCase "Fake" (pure ())) + let n_per_bin = lookupOption opts :: NPerBin + p_val = lookupOption opts + -- + g0 <- createSystemRandom + defaultMainWithIngredients ingredients $ testGroup "mwc" + [ testProperty "save/restore" $ prop_SeedSaveRestore g0 +#if MIN_VERSION_random(1,3,0) + , testProperty "SeedGen" $ prop_SeedGen g0 +#endif + , testCase "user save/restore" $ saveRestoreUserSeed + , testCase "empty seed data" $ emptySeed + , testCase "output correct" $ do + g <- create + xs <- replicateM 513 (uniform g) + assertEqual "[Word32]" xs golden + , testCase "beta binomial mean" $ prop_betaBinomialMean + , testProperty "binomial is binomial" $ prop_binomial_PMF n_per_bin p_val g0 + , testProperty "Poisson is Poisson" $ prop_Poisson_PMF n_per_bin p_val g0 + , testCase "poisson mean and variance of 1" $ prop_poissonMeanAndVar 1 + , testCase "poisson mean and variance of 5" $ prop_poissonMeanAndVar 5 + , testCase "poisson mean and variance of 40" $ prop_poissonMeanAndVar 40 + , testCase "poisson mean and variance of 150" $ prop_poissonMeanAndVar 150 + , testCase "poisson 0" $ do + n <- poisson 0 g0 + 0 @=? n + ] + +updateGenState :: GenIO -> IO () +updateGenState g = replicateM_ 250 (uniform g :: IO Word32) + +prop_SeedSaveRestore :: GenIO -> Property +prop_SeedSaveRestore g = monadicIO $ do + run $ updateGenState g + seed <- run $ save g + seed' <- run $ save =<< restore seed + return $ seed == seed' + +#if MIN_VERSION_random(1,3,0) +prop_SeedGen :: GenIO -> Property +prop_SeedGen g = monadicIO $ do + run $ updateGenState g + seed <- run $ save g + return $ seed == (Random.fromSeed . Random.toSeed) seed +#endif + +saveRestoreUserSeed :: IO () +saveRestoreUserSeed = do + let seed = toSeed $ U.replicate 258 0 + seed' <- save =<< restore seed + assertEqual "Seeds must be equal" seed' seed + +emptySeed :: IO () +emptySeed = do + let seed = toSeed U.empty + seed' <- save =<< create + assertEqual "Seeds must be equal" seed' seed + +-- First 513 values generated from seed made using create +golden :: [Word32] +golden = + [ 2254043345, 562229898, 1034503294, 2470032534, 2831944869, 3042560015, 838672965, 715056843 + , 3122641307, 2300516242, 4079538318, 3722020688, 98524204, 1450170923, 2669500465, 2890402829 + , 114212910, 1914313000, 2389251496, 116282477, 1771812561, 1606473512, 1086420617, 3652430775 + , 1165083752, 3599954795, 3006722175, 341614641, 3000394300, 1378097585, 1551512487, 81211762 + , 604209599, 3949866361, 77745071, 3170410267, 752447516, 1213023833, 1624321744, 3251868348 + , 1584957570, 2296897736, 3305840056, 1158966242, 2458014362, 1919777052, 3203159823, 3230279656 + , 755741068, 3005087942, 2478156967, 410224731, 1196248614, 3302310440, 3295868805, 108051054 + , 1010042411, 2725695484, 2201528637, 667561409, 79601486, 50029770, 566202616, 3217300833 + , 2162817014, 925506837, 1527015413, 3079491438, 927252446, 118306579, 499811870, 2973454232 + , 2979271640, 4078978924, 1864075883, 197741457, 296365782, 1784247291, 236572186, 464208268 + , 1769568958, 827682258, 4247376295, 2959098022, 1183860331, 2475064236, 3952901213, 1953014945 + , 393081236, 1616500498, 2201176136, 1663813362, 2167124739, 630903810, 113470040, 924745892 + , 1081531735, 4039388931, 4118728223, 107819176, 2212875141, 1941653033, 3660517172, 192973521 + , 3653156164, 1878601439, 3028195526, 2545631291, 3882334975, 456082861, 2775938704, 3813508885 + , 1758481462, 3332769695, 3595846251, 3745876981, 152488869, 2555728588, 3058747945, 39382408 + , 520595021, 2185388418, 3502636573, 2650173199, 1077668433, 3548643646, 71562049, 2726649517 + , 494210825, 1208915815, 620990806, 2877290965, 3253243521, 804166732, 2481889113, 623399529 + , 44880343, 183645859, 3283683418, 2214754452, 419328482, 4224066437, 1102669380, 1997964721 + , 2437245376, 985749802, 858381069, 116806511, 1771295365, 97352549, 341972923, 2971905841 + , 110707773, 950500868, 1237119233, 691764764, 896381812, 1528998276, 1269357470, 2567094423 + , 52141189, 2722993417, 80628658, 3919817965, 3615946076, 899371181, 46940285, 4010779728 + , 318101834, 30736609, 3577200709, 971882724, 1478800972, 3769640027, 3706909300, 3300631811 + , 4057825972, 4285058790, 2329759553, 2967563409, 4080096760, 2762613004, 2518395275, 295718526 + , 598435593, 2385852565, 2608425408, 604857293, 2246982455, 919156819, 1721573814, 2502545603 + , 643962859, 587823425, 3508582012, 1777595823, 4119929334, 2833342174, 414044876, 2469473258 + , 289159600, 3715175415, 966867024, 788102818, 3197534326, 3571396978, 3508903890, 570753009 + , 4273926277, 3301521986, 1411959102, 2766249515, 4071012597, 959442028, 1962463990, 1098904190 + , 714719899, 562204808, 1658783410, 1471669042, 2565780129, 1616648894, 4236521717, 1788863789 + , 3068674883, 191936470, 253084644, 1915647866, 276372665, 2117183118, 3704675319, 218791054 + , 3680045802, 406662689, 3844864229, 91140313, 3834015630, 25116147, 904830493, 3152559113 + , 820358622, 1301896358, 296152699, 2202014455, 4256659428, 1175171414, 3287520873, 2028006499 + , 327448717, 2095642873, 3798661296, 58567008, 3907537112, 3691259011, 1730142328, 2373011713 + , 3387040741, 3189417655, 2949233059, 1238379614, 1813238023, 1064726446, 1339055235, 1744523609 + , 279811576, 2934103599, 283542302, 994488448, 418691747, 1062780152, 102211875, 4071713296 + , 1790834038, 1035092527, 2374272359, 3558280982, 1927663822, 3645417844, 3481790745, 3566282546 + , 2000290859, 505518126, 363501589, 4075468679, 3247300709, 3705242654, 2731103609, 2836871038 + , 589640144, 2546495106, 84767518, 1376911639, 2400770705, 527489676, 3804134352, 150084021 + , 240070593, 3807594859, 3518576690, 659503830, 2239678479, 1273668921, 4271050554, 3090482972 + , 401956859, 1772128561, 4438455, 1989666158, 2521484677, 3960178700, 4220196277, 1033999035 + , 2214785840, 3428469341, 428564336, 2517446784, 3935757188, 3294001677, 1037971963, 3590324170 + , 1220969729, 1719719817, 807688972, 77076422, 4251553858, 3963852375, 326128795, 3277818295 + , 3671513069, 549617771, 1683950556, 3352913781, 409318429, 2456264774, 4036950639, 1162718475 + , 83888874, 5578966, 172866494, 1542278848, 455546979, 1296511553, 4263636440, 2450589064 + , 372411483, 211216338, 2632256495, 2393754408, 1336054289, 4087203071, 3159642437, 1933346856 + , 2914152714, 3805541979, 2769740793, 1161287028, 2289749561, 4124509890, 2128452935, 210531695 + , 4250709834, 390950534, 1421430300, 3030519715, 3228987297, 3086837053, 2866915453, 2335948692 + , 1684378991, 2575634059, 4153427304, 2426048796, 4197556954, 2605152326, 2909410733, 2424889219 + , 654577921, 811955499, 118126602, 504071559, 1278756230, 3896458168, 4105558075, 750276169 + , 1120805572, 1762689330, 993728154, 1104363215, 774344996, 4077568952, 2183487324, 994724370 + , 3323036885, 3880704963, 746305447, 961608310, 2030117337, 453935768, 800490463, 1034636 + , 2323633564, 602565693, 806061242, 1899269713, 162686347, 467541008, 1529175313, 282891502 + , 2529616339, 2930657178, 464272784, 2878535316, 807165854, 3209080518, 4080120278, 347748171 + , 3972126063, 284174728, 2498328933, 1723872460, 143845955, 4223866687, 1761495357, 1544646770 + , 4206103283, 3771574626, 642165282, 1119501013, 3514063332, 1443320304, 4056369796, 3602131475 + , 1422908288, 804093687, 431176780, 40108717, 2998264213, 3705835674, 169805085, 454593842 + , 2781536994, 2385225212, 4137367775, 2631435125, 2347082354, 629238010, 3283635219, 3815791831 + , 1340400558, 4061846985, 3803921868, 3196119096, 718610843, 3694290834, 2169960411, 2407155570 + , 2557480499, 16164105, 480957288, 2155919829, 2490067282, 2356287132, 511737296, 1602800634 + , 1802275249, 3316832299, 50286484, 2106622541, 2352302834, 2538374315, 344766394, 2777260569 + , 1215135803, 2229011963, 114632277, 1645499402, 1111617833, 3833259754, 928611385, 686744723 + , 1898396834, 2461932251, 2665457318, 3797019621, 868313114, 2366635205, 481934875, 1170532970 + , 642610859, 3150733309, 3508548582, 666714469, 711663449, 2436617656, 2681476315, 1637296693 + , 2487349478, 4174144946, 2793869557, 559398604, 1898140528, 991962870, 864792875, 3861665129 + , 4024051364, 3383200293, 773730975, 33517291, 2660126073, 689133464, 2248134097, 3874737781 + , 3358012678] + +-- We can test two for the price of one +betaBinomial :: StatefulGen g m => + Double -> Double -> Int -> g -> m Int +betaBinomial a b n g = do + p <- beta a b g + binomial n p g + +nSamples :: Int +nSamples = 10000 + +alpha, delta :: Double +alpha = 600.0 +delta = 400.0 + +nTrials :: Int +nTrials = 10 + +prop_betaBinomialMean :: IO () +prop_betaBinomialMean = do + g <- create + ss <- replicateM nSamples $ betaBinomial alpha delta nTrials g + let m = fromIntegral (sum ss) / fromIntegral nSamples + let x1 = fromIntegral nTrials * alpha / (alpha + delta) + assertBool ("Mean is " ++ show x1 ++ " but estimated as " ++ show m) (abs (m - x1) < 0.001) + +prop_poissonMeanAndVar :: Double -> IO () +prop_poissonMeanAndVar lambda = do + gen <- create + ss <- replicateM nSamples $ poisson lambda gen + let m = fromIntegral (sum ss) / fromIntegral nSamples :: Double + let v = (fromIntegral (sum (map (^ 2) ss)) / fromIntegral nSamples) - (m ** 2) + assertBool + ("True mean and var: " ++ show lambda ++ " but estimated as " ++ show (m, v) ++ " respectively") + ((abs (lambda - m) / lambda < 0.025) && (abs (lambda - v) / lambda < 0.025)) + + +-- Test that `binomial` really samples from binomial distribution. +-- +-- If we have binomial random variate with number of trials N and +-- sample it M times. Then number of events with K successes is +-- described by multinomial distribution and we can test whether +-- experimental distribution is described using likelihood ratio test +prop_binomial_PMF :: NPerBin -> PValue -> GenIO -> Property +prop_binomial_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do + p <- choose (0, 1.0) -- Success probability + n_trial <- choose (2, 100) -- Number of trials in binomial distribution + -- Number of binomial samples to generate + let n_samples = n_trial * n_per_bin + n_samples' = fromIntegral n_samples + -- Compute number of outcomes + pure $ ioProperty $ do + hist <- do + buf <- MVU.new (n_trial + 1) + replicateM_ n_samples $ + MVU.modify buf (+(1::Int)) =<< binomial n_trial p g + U.unsafeFreeze buf + -- Here we compute twice log of likelihood ratio. Alternative + -- hypothesis is some distribution which fits data perfectly + -- + -- Asymtotically it's ditributed as χ² with n_trial-1 degrees of + -- freedom + let likelihood _ 0 + = 0 + likelihood k (fromIntegral -> n_obs) + = n_obs * (log (n_obs / n_samples') - logProbBinomial n_trial p k) + let logL = 2 * U.sum (U.imap likelihood hist) + let significance = 1 - cumulativeChi2 (n_trial - 1) logL + pure $ counterexample ("p = " ++ show p) + $ counterexample ("N = " ++ show n_trial) + $ counterexample ("p-val = " ++ show significance) + $ counterexample ("chi2 = " ++ show logL) + $ significance > p_val + +-- Similar test for Poisson distribution. We build histogram by moving +-- all values >λ+2σ into overflow bin. +prop_Poisson_PMF :: NPerBin -> PValue -> GenIO -> Property +prop_Poisson_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do + lam <- choose (0, 100.0) -- Poisson rate + let max_n = ceiling $ lam + 2 * sqrt lam -- Maximum histogram bin + n_samples = max_n * n_per_bin -- Number of samples to generate + n_samples' = fromIntegral n_samples + pure $ ioProperty $ do + hist <- do + buf <- MVU.new (max_n + 1) + replicateM_ n_samples $ do + k <- poisson lam g + MVU.modify buf (+(1::Int)) (min max_n k) + U.unsafeFreeze buf + let likelihood _ 0 + = 0 + likelihood k (fromIntegral -> n_obs) + = n_obs * (log (n_obs / n_samples') - p) + where + p | k == max_n = logComlCdfPoisson lam k + | otherwise = logProbPoisson lam k + let logL = 2 * U.sum (U.imap likelihood hist) + let significance = 1 - cumulativeChi2 max_n logL + pure $ counterexample ("lambda = " ++ show lam) + $ counterexample ("p-val = " ++ show significance) + $ counterexample ("chi2 = " ++ show logL) + $ significance > p_val + +---------------------------------------------------------------- +-- Statistical helpers +---------------------------------------------------------------- + +-- Logarithm of probability for binomial distribution +logProbBinomial :: Int -> Double -> Int -> Double +logProbBinomial n p k + = logChoose n k + log p * k' + log1p (-p) * nk' + where + k' = fromIntegral k + nk' = fromIntegral $ n - k + +logProbPoisson :: Double -> Int -> Double +logProbPoisson lam i = log lam * fromIntegral i - logFactorial i - lam + +logComlCdfPoisson :: Double -> Int -> Double +logComlCdfPoisson lam i = log $ incompleteGamma (fromIntegral (i :: Int)) lam + + +cumulativeChi2 :: Int -> Double -> Double +cumulativeChi2 (fromIntegral -> ndf) x + | x <= 0 = 0 + | otherwise = incompleteGamma (ndf/2) (x/2)