diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,7 +1,26 @@ +Changes in 0.13.2.0 + + * Generators for beta, Bernoully, Dirichlet and categorical distributions + added. + + * Functions for generating random shuffles added. + + +Changes in 0.13.1.2 + + * GHC 7.9 support + + +Changes in 0.13.1.1 + + * Long standing performance problem in normal distribution fixed (#16) + + Changes in 0.13.1.0 * createSystemRandom added + Changes in 0.13.0.0 * Workaround for GHC bug 8072 (bug 25). GHC 7.6 on 32-bit platrofms is @@ -9,6 +28,7 @@ * Generators for truncated exponential and geometric distributions added. + Changes in 0.12.0.0 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,4 @@ -{-# LANGUAGE BangPatterns, GADTs #-} +{-# LANGUAGE BangPatterns, GADTs, FlexibleContexts, ScopedTypeVariables #-} -- | -- Module : System.Random.MWC.Distributions -- Copyright : (c) 2012 Bryan O'Sullivan @@ -13,25 +13,40 @@ module System.Random.MWC.Distributions ( -- * Variates: non-uniformly distributed values + -- ** Continuous distributions normal , standard , exponential , truncatedExp , gamma , chiSquare + , beta + -- ** Discrete distribution + , categorical , geometric0 , geometric1 + , bernoulli + -- ** Multivariate + , dirichlet + -- * Permutations + , uniformPermutation + , uniformShuffle -- * References -- $references ) where -import Control.Monad (liftM) +import Prelude hiding (mapM) +import Control.Monad (liftM,when) import Control.Monad.Primitive (PrimMonad, PrimState) import Data.Bits ((.&.)) +import Data.Foldable (Foldable,foldl') +import Data.Traversable (Traversable,mapM) import Data.Word (Word32) -import System.Random.MWC (Gen, uniform) -import qualified Data.Vector.Unboxed as I +import System.Random.MWC (Gen, uniform, uniformR) +import qualified Data.Vector.Unboxed as I +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as M -- Unboxed 2-tuple data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double @@ -117,9 +132,9 @@ -> Gen (PrimState m) -- ^ Generator -> m Double {-# INLINE exponential #-} -exponential beta gen = do +exponential b gen = do x <- uniform gen - return $! - log x / beta + return $! - log x / b -- | Generate truncated exponentially distributed random variate. @@ -130,12 +145,12 @@ -> Gen (PrimState m) -- ^ Generator. -> m Double {-# INLINE truncatedExp #-} -truncatedExp beta (a,b) gen = do +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 - return $! a - log ( (1 - p) + p*exp(-beta*delta)) / beta + return $! a - log ( (1 - p) + p*exp(-scale*delta)) / scale -- | Random variate generator for gamma distribution. gamma :: PrimMonad m @@ -208,7 +223,92 @@ geometric1 p gen = do n <- geometric0 p gen return $! n + 1 +-- | Random variate generator for Beta distribution +beta :: PrimMonad m + => Double -- ^ alpha (>0) + -> Double -- ^ beta (>0) + -> Gen (PrimState m) -- ^ Generator + -> m Double +{-# INLINE beta #-} +beta a b gen = do + x <- gamma a 1 gen + y <- gamma b 1 gen + return $! x / (x+y) +-- | Random variate generator for Dirichlet distribution +dirichlet :: (PrimMonad m, Traversable t) + => t Double -- ^ container of parameters + -> Gen (PrimState m) -- ^ Generator + -> m (t Double) +{-# INLINE dirichlet #-} +dirichlet t gen = do + t' <- mapM (\x -> gamma x 1 gen) t + let total = foldl' (+) 0 t' + return $ fmap (/total) t' + +-- | Random variate generator for Bernoulli distribution +bernoulli :: PrimMonad m + => Double -- ^ Probability of success (returning True) + -> Gen (PrimState m) -- ^ Generator + -> m Bool +{-# INLINE bernoulli #-} +bernoulli p gen = ( v Double -- ^ List of weights [>0] + -> Gen (PrimState m) -- ^ 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 + return $! case G.findIndex (>=p) cv of + Just i -> i + Nothing -> pkgError "categorical" "bad weights!" + +-- | Random variate generator for uniformly distributed permutations. +-- 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) + => Int + -> Gen (PrimState m) + -> m (v Int) +{-# INLINE uniformPermutation #-} +uniformPermutation n gen = do + when (n<=0) (pkgError "uniformPermutation" "size must be >0") + v <- G.unsafeThaw (G.generate n id :: v Int) + let lst = n-1 + loop i | i == lst = G.unsafeFreeze v + | otherwise = do + j <- uniformR (i,lst) gen + M.unsafeSwap v i j + loop (i+1) + loop 0 + + +-- | Random variate generator for a uniformly distributed shuffle of a +-- vector. +uniformShuffle :: (PrimMonad m, G.Vector v a, G.Vector v Int) + => v a + -> Gen (PrimState m) + -> m (v a) +{-# INLINE uniformShuffle #-} +uniformShuffle xs gen + | G.length xs <= 1 = return xs + | otherwise = do + idx <- uniformPermutation (G.length xs) gen + return $! G.backpermute xs idx + + sqr :: Double -> Double sqr x = x * x {-# INLINE sqr #-} @@ -216,6 +316,8 @@ pkgError :: String -> String -> a pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func ++ ": " ++ msg + + -- $references -- diff --git a/benchmarks/tsts.hs b/benchmarks/tsts.hs new file mode 100644 --- /dev/null +++ b/benchmarks/tsts.hs @@ -0,0 +1,65 @@ +{-# LANGUAGE BangPatterns #-} +import Control.Monad +import System.Random.MWC +import System.Random.MWC.Distributions + +main = do + withSystemRandom $ \g -> replicateM_ (300*1000) $ do + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + return () :: IO () + diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,5 +1,5 @@ name: mwc-random -version: 0.13.1.2 +version: 0.13.2.0 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random diff --git a/test/KS.hs b/test/KS.hs --- a/test/KS.hs +++ b/test/KS.hs @@ -10,10 +10,12 @@ import Statistics.Test.KolmogorovSmirnov import Statistics.Distribution +import Statistics.Distribution.Binomial import Statistics.Distribution.Exponential import Statistics.Distribution.Gamma import Statistics.Distribution.Normal import Statistics.Distribution.Uniform +import Statistics.Distribution.Beta import qualified System.Random.MWC as MWC import qualified System.Random.MWC.Distributions as MWC @@ -38,6 +40,10 @@ -- Exponential , testCase "exponential l=1" $ testKS (exponential 1) (MWC.exponential 1) g , testCase "exponential l=3" $ testKS (exponential 3) (MWC.exponential 3) g + -- Beta + , testCase "beta a=0.3,b=0.5" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g + , testCase "beta a=0.1,b=0.8" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g + , testCase "beta a=0.8,b=0.1" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g ] testKS :: (Distribution d) => d -> (MWC.GenIO -> IO Double) -> MWC.GenIO -> IO ()