packages feed

random-fu 0.0.2.0 → 0.0.2.1

raw patch · 16 files changed

+46/−117 lines, 16 filesdep +erfdep +random-shuffle

Dependencies added: erf, random-shuffle

Files

random-fu.cabal view
@@ -1,5 +1,5 @@ name:                   random-fu-version:                0.0.2.0+version:                0.0.2.1 stability:              experimental  cabal-version:          >= 1.2@@ -16,9 +16,7 @@                         for entropy sources and random variable distributions, all                         served up on a monadic platter.  Aspires to be useful                         in an idiomatic way in both \"pure\" and \"impure\" styles,-                        as well as reasonably fast.  May not yet meet the latter-                        goal, but I think the former is starting to shape up-                        nicely.+                        as well as reasonably fast.  Flag base4 @@ -51,17 +49,19 @@                         Data.Random.Source.PureMT                         Data.Random.Source.Std   if flag(base4)-    build-depends:      base >= 4,+    build-depends:      base >= 4 && <5,                         syb   else     build-depends:      base >= 3 && < 4        build-depends:        array,                         containers,+                        erf,                         mersenne-random-pure64,                         monad-loops >= 0.3.0.1,                         mtl,                         random,+                        random-shuffle,                         stateref,                         storablevector,                         template-haskell
src/Data/Random/Distribution.hs view
@@ -9,9 +9,6 @@  import Data.Random.Lift import Data.Random.RVar-import Data.Random.Source-import Data.Random.Source.Std-import Data.Word  -- |A definition of a random variable's distribution.  From the distribution -- an 'RVar' can be created, or the distribution can be directly sampled using @@ -23,7 +20,7 @@  class Distribution d t => CDF d t where     -- |Return the cumulative distribution function of this distribution.-    -- That is, a function taking 'x :: t' to the probability that the next+    -- That is, a function taking @x :: t@ to the probability that the next     -- sample will return a value less than or equal to x, according to some     -- order or partial order (not necessarily an obvious one).     --@@ -31,7 +28,7 @@     -- to the CDF with respect to that order.     --      -- In other cases, 'cdf' is only required to satisfy the following law:-    -- > fmap (cdf d) (rvar d)+    -- @fmap (cdf d) (rvar d)@     -- must be uniformly distributed over (0,1).  Inclusion of either endpoint is optional,     -- though the preferred range is (0,1].     --
src/Data/Random/Distribution/Bernoulli.hs view
@@ -12,10 +12,8 @@  import Data.Random.Internal.TH -import Data.Random.Source-import Data.Random.Distribution import Data.Random.RVar-+import Data.Random.Distribution import Data.Random.Distribution.Uniform  import Data.Int
src/Data/Random/Distribution/Beta.hs view
@@ -9,7 +9,6 @@  module Data.Random.Distribution.Beta where -import Data.Random.Source import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Gamma@@ -29,12 +28,6 @@     x <- erlang a     y <- erlang b     return (x / (x + y))--fractionalBetaMinusHalf :: (Fractional a, Distribution Gamma a) => a -> a -> RVar a-fractionalBetaMinusHalf a b = do-    x <- gamma a 1-    y <- gamma b 1-    return ((x-y)/(2*(x+y)))  beta :: Distribution Beta a => a -> a -> RVar a beta a b = rvar (Beta a b)
src/Data/Random/Distribution/Binomial.hs view
@@ -11,11 +11,8 @@  import Data.Random.Internal.TH -import Data.Random.Source-import Data.Random.Distribution import Data.Random.RVar--import Data.Random.Distribution.Bernoulli+import Data.Random.Distribution import Data.Random.Distribution.Beta import Data.Random.Distribution.Uniform @@ -43,15 +40,14 @@         bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVar a         bin k t p             | t > 10    = do-                a <- selectA (fromInteger (floatRadix undefined) ^ floatDigits undefined) t-                    -                let b = 1 + t - a-                +                let a = 1 + t `div` 2+                    b = 1 + t - a+                         x <- beta (fromIntegral a) (fromIntegral b)                 if x >= p                     then bin  k      (a - 1) (p / x)                     else bin (k + a) (b - 1) ((p - x) / (1 - x))-            +                     | otherwise = count k t                 where                     count k  0    = return k
src/Data/Random/Distribution/Exponential.hs view
@@ -9,7 +9,6 @@  module Data.Random.Distribution.Exponential where -import Data.Random.Source import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Uniform
src/Data/Random/Distribution/Gamma.hs view
@@ -16,7 +16,6 @@  module Data.Random.Distribution.Gamma where -import Data.Random.Source import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Uniform
src/Data/Random/Distribution/Normal.hs view
@@ -34,11 +34,7 @@ import Control.Monad import Foreign.Storable -foreign import ccall "math.h erf" erf :: Double -> Double-foreign import ccall "math.h erfc" erfc :: Double -> Double-foreign import ccall "math.h erff" erff :: Float -> Float-erfg :: RealFrac a => a -> a-erfg = realToFrac . erf . realToFrac+import Data.Number.Erf  normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a) normalPair = boxMullerNormalPair@@ -86,11 +82,11 @@                 else return (r - x)  -- |Construct a 'Ziggurat' for sampling a normal distribution, given--- a suitable error function, logBase 2 c, and the 'zGetIU' implementation.+-- logBase 2 c, and the 'zGetIU' implementation. normalZ ::-  (RealFloat a, Storable a, Distribution Uniform a, Integral b) =>-  (a -> a) -> b -> RVar (Int, a) -> Ziggurat a-normalZ erf p = mkZigguratRec True normalF normalFInv (normalFInt erf) normalFVol (2^p)+  (RealFloat a, Erf a, Storable a, Distribution Uniform a, Integral b) =>+  b -> RVar (Int, a) -> Ziggurat a+normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)  -- | Ziggurat target function normalF :: (Floating a, Ord a) => a -> a@@ -100,17 +96,17 @@ -- | inverse of 'normalF' normalFInv :: Floating a => a -> a normalFInv y  = sqrt ((-2) * log y)--- | integral of 'normalF', parameterized over 'erf'-normalFInt :: (Floating a, Ord a) => (a -> a) -> a -> a-normalFInt erf x +-- | integral of 'normalF'+normalFInt :: (Floating a, Erf a, Ord a) => a -> a+normalFInt x      | x <= 0    = 0     | otherwise = normalFVol * erf (x * sqrt 0.5) -- | volume of 'normalF' normalFVol :: Floating a => a normalFVol = sqrt (0.5 * pi) -realFloatStdNormal :: (RealFloat a, Storable a, Distribution Uniform a) => RVar a-realFloatStdNormal = rvar (normalZ erfg p getIU)+realFloatStdNormal :: (RealFloat a, Erf a, Storable a, Distribution Uniform a) => RVar a+realFloatStdNormal = runZiggurat (normalZ p getIU)     where          p = 6         @@ -120,7 +116,7 @@             return (fromIntegral i .&. (2^p-1), u)  doubleStdNormal :: RVar Double-doubleStdNormal = rvar doubleStdNormalZ+doubleStdNormal = runZiggurat doubleStdNormalZ  -- doubleStdNormalC must not be over 12 if using wordToDoubleWithExcess doubleStdNormalC :: Int@@ -142,7 +138,7 @@             return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)  floatStdNormal :: RVar Float-floatStdNormal = rvar floatStdNormalZ+floatStdNormal = runZiggurat floatStdNormalZ  floatStdNormalC :: Int floatStdNormalC = 512@@ -168,29 +164,33 @@ normalPdf :: Real a => a -> a -> a -> Double normalPdf m s x = recip (realToFrac s * sqrt (2*pi)) * exp (-0.5 * (realToFrac x - realToFrac m)^2 / (realToFrac s)^2) -normalCdf :: Real a => a -> a -> a -> Double-normalCdf m s x = 0.5 * (1 + erf ((realToFrac x - realToFrac m) / (realToFrac s * sqrt 2)))+normalCdf :: (Real a, Erf a) => a -> a -> a -> Double+normalCdf m s x = realToFrac (normcdf ((x-m) / s))  data Normal a     = StdNormal     | Normal a a -- mean, sd  instance Distribution Normal Double where+    {-# SPECIALIZE instance Distribution Normal Double #-}     rvar StdNormal = doubleStdNormal     rvar (Normal m s) = do         x <- doubleStdNormal         return (x * s + m)  instance Distribution Normal Float where+    {-# SPECIALIZE instance Distribution Normal Float #-}     rvar StdNormal = floatStdNormal     rvar (Normal m s) = do         x <- floatStdNormal         return (x * s + m) -instance (Real a, Distribution Normal a) => CDF Normal a where+instance (Real a, Erf a, Distribution Normal a) => CDF Normal a where     cdf StdNormal    = normalCdf 0 1     cdf (Normal m s) = normalCdf m s +{-# SPECIALIZE stdNormal :: RVar Double #-}+{-# SPECIALIZE stdNormal :: RVar Float #-} stdNormal :: Distribution Normal a => RVar a stdNormal = rvar StdNormal
src/Data/Random/Distribution/Poisson.hs view
@@ -11,17 +11,14 @@  import Data.Random.Internal.TH -import Data.Random.Source-import Data.Random.Distribution import Data.Random.RVar-+import Data.Random.Distribution import Data.Random.Distribution.Uniform import Data.Random.Distribution.Gamma import Data.Random.Distribution.Binomial  import Data.Int import Data.Word- import Control.Monad  -- from Knuth, with interpretation help from gsl sources
src/Data/Random/Distribution/Uniform.hs view
@@ -37,10 +37,8 @@ import Data.Random.Distribution import Data.Random.RVar -import Data.Ratio import Data.Word import Data.Int-import Data.Bits import Data.List  import Control.Monad.Loops
src/Data/Random/Internal/Words.hs view
@@ -10,7 +10,6 @@  import Data.Bits import Data.Word-import Control.Monad  {-# INLINE buildWord #-} buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64
src/Data/Random/List.hs view
@@ -1,19 +1,10 @@-{-- -      ``Data/Random/List''- -}-{-# LANGUAGE -    FlexibleContexts-  #-}- module Data.Random.List where  import Data.Random.RVar-import Data.Random.Source-import Data.Random.Distribution import Data.Random.Distribution.Uniform-import GHC.IOBase -import qualified Data.Sequence as S+import qualified System.Random.Shuffle as SRS+import Control.Monad  randomElement :: [a] -> RVar a randomElement [] = error "randomElement: empty list!"@@ -21,47 +12,16 @@     n <- uniform 0 (length xs - 1)     return (xs !! n) -randomSeqElement :: S.Seq a -> RVar a-randomSeqElement s-    | S.null s  = error "randomSeqElement: empty list!"-    | otherwise = do-        n <- uniform 0 (S.length s - 1)-        return (s `S.index` n)- shuffle :: [a] -> RVar [a]-shuffle = shuffleSeq . S.fromList--shuffleSeq :: S.Seq a -> RVar [a]-shuffleSeq s = shuffle (S.length s) s-    where-        shuffle 0 _ = return []-        shuffle (n+1) s = do-            i <- uniform 0 n-            let (x, xs) = extract i s-            ys <- shuffle n xs-            return (x:ys)-        -        extract n s = case S.splitAt n s of-            (l,r) -> case S.viewl r of-                x S.:< r  -> (x, l S.>< r)---- |Shuffle a list using interleaved IO when extracting elements.-lazyShuffleFrom :: (RandomSource IO s) => s -> [a] -> IO [a]-lazyShuffleFrom src = lazyShuffleSeqFrom src . S.fromList+shuffle [] = return []+shuffle xs = do+    is <- zipWithM (\_ i -> uniform 0 i) (tail xs) [1..]+    +    return (SRS.shuffle xs (reverse is)) --- |Shuffle a 'S.Seq' using interleaved IO when extracting elements.-lazyShuffleSeqFrom :: (RandomSource IO s) => s -> S.Seq a -> IO [a]-lazyShuffleSeqFrom src s = shuffle (S.length s) s-    where-        shuffle 0     _  = return []-        shuffle (n+1) s -            | S.null s = return []-            | otherwise = do-                i <- runRVar (uniform 0 n) src-                let (x, xs) = extract i s-                ys <- unsafeInterleaveIO (shuffle n xs)-                return (x:ys)-        -        extract n s = case S.splitAt n s of-            (l,r) -> case S.viewl r of-                x S.:< r  -> (x, l S.>< r)+shuffleN :: Int -> [a] -> RVar [a]+shuffleN 0 xs = return []+shuffleN (n+1) xs = do+    is <- sequence [uniform 0 i | i <- [n,n-1..1]]+    return (SRS.shuffle xs is)+
src/Data/Random/RVar.hs view
@@ -23,16 +23,13 @@     ) where  -import Data.Random.Internal.Words import Data.Random.Source import Data.Random.Lift as L -import Data.Word import Data.Bits  import qualified Control.Monad.Trans as T import Control.Applicative-import Control.Monad import Control.Monad.Reader import Control.Monad.Identity
src/Data/Random/Source.hs view
@@ -11,8 +11,6 @@     ) where  import Data.Word-import Data.Bits-import Data.List import Control.Monad  import Data.Random.Internal.Words
src/Data/Random/Source/DevRandom.hs view
@@ -11,7 +11,6 @@  import Data.Random.Source -import GHC.IOBase (unsafePerformIO) import System.IO (openBinaryFile, hGetBuf, IOMode(..)) import Foreign
src/Data/Random/Source/StdGen.hs view
@@ -10,7 +10,6 @@ import Data.Random.Internal.Words import Data.Random.Source import System.Random-import Control.Monad import Control.Monad.State import qualified Control.Monad.ST.Strict as S import qualified Control.Monad.State.Strict as S