random-fu 0.0.3 → 0.0.3.2

raw patch · 17 files changed

+386/−95 lines, 17 filesdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• random-fu.cabal +10/−3
• src/Data/Random.hs +3/−3
• src/Data/Random/Distribution.hs +3/−2
• src/Data/Random/Distribution/Normal.hs +40/−4
• src/Data/Random/Distribution/Rayleigh.hs +4/−7
• src/Data/Random/Distribution/Triangular.hs +17/−5
• src/Data/Random/Distribution/Uniform.hs +50/−6
• src/Data/Random/Distribution/Ziggurat.hs +151/−56
• src/Data/Random/Internal/Find.hs +1/−1
• src/Data/Random/Internal/Fixed.hs +22/−0
• src/Data/Random/Internal/TH.hs +43/−1
• src/Data/Random/Internal/Words.hs +7/−2
• src/Data/Random/List.hs +13/−1
• src/Data/Random/RVar.hs +11/−0
• src/Data/Random/Sample.hs +1/−1
• src/Data/Random/Source/PureMT.hs +5/−3
• src/Data/Random/Source/StdGen.hs +5/−0

@@ -1,5 +1,5 @@ name:                   random-fu-version:                0.0.3+version:                0.0.3.2 stability:              experimental  cabal-version:          >= 1.2@@ -19,6 +19,8 @@                         as well as reasonably fast.  Flag base4+Flag base4_2+    Description:        base-4.2 has an incompatible change in Data.Fixed (HasResolution)  Library   hs-source-dirs:       src@@ -50,8 +52,13 @@                         Data.Random.Source.Std                         Data.Random.Source.StdGen   if flag(base4)-    build-depends:      base >= 4 && <5,-                        syb+    build-depends:      syb+    +    if flag(base4_2)+      build-depends:    base >= 4 && <4.2+    else+      cpp-options:      -Dbase_4_2+      build-depends:    base >= 4.2 && <5   else     build-depends:      base >= 3 && < 4     

@@ -7,14 +7,14 @@  -- |Random numbers and stuff... -- --- Data.Random.Source exports the typeclasses for entropy sources, and+-- "Data.Random.Source" exports the typeclasses for entropy sources, and -- Data.Random.Source.* export various instances and/or functions with which -- instances can be defined. -- --- Data.Random.Distribution exports the typeclasses for sampling distributions,+-- "Data.Random.Distribution" exports the typeclasses for sampling distributions, -- and Data.Random.Distribution.* export various specific distributions. ----- Data.Random.RVar exports the 'RVar' type, which is a probability distribution+-- "Data.Random.RVar" exports the 'RVar' type, which is a probability distribution -- monad that allows for concise definitions of random variables, as well as -- a couple handy 'RVar's. 

@@ -16,7 +16,6 @@ class Distribution d t where     -- |Return a random variable with this distribution.     rvar :: d t -> RVar t-    rvar = rvarT  class Distribution d t => CDF d t where     -- |Return the cumulative distribution function of this distribution.@@ -34,7 +33,9 @@     --      -- Thus, 'cdf' for a product type should not be a joint CDF as commonly      -- defined, as that definition violates both conditions.-    -- Instead, it should be a univariate CDF over the product type.+    -- Instead, it should be a univariate CDF over the product type.  That is,+    -- it should represent the CDF with respect to the lexicographic order+    -- of the tuple.     cdf :: d t -> t -> Double  -- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.

@@ -36,9 +36,17 @@  import Data.Number.Erf +-- |A random variable that produces a pair of independent+-- normally-distributed values. normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a) normalPair = boxMullerNormalPair +-- |A random variable that produces a pair of independent+-- normally-distributed values, computed using the Box-Muller method.+-- This algorithm is slightly slower than Knuth's method but using a +-- constant amount of entropy (Knuth's method is a rejection method).+-- It is also slightly more general (Knuth's method require an 'Ord'+-- instance). {-# INLINE boxMullerNormalPair #-} boxMullerNormalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a) boxMullerNormalPair = do@@ -51,6 +59,10 @@         y = r * sin theta     return (x,y) +-- |A random variable that produces a pair of independent+-- normally-distributed values, computed using Knuth's polar method.+-- Slightly faster than 'boxMullerNormalPair' when it accepts on the +-- first try, but does not always do so. {-# INLINE knuthPolarNormalPair #-} knuthPolarNormalPair :: (Floating a, Ord a, Distribution Uniform a) => RVar (a,a) knuthPolarNormalPair = do@@ -65,8 +77,7 @@             else let scale = sqrt (-2 * log s / s)                    in (v1 * scale, v2 * scale) --- |Draw from the tail of a normal distribution (the region beyond the provided value), --- returning a negative value if the Bool parameter is True.+-- |Draw from the tail of a normal distribution (the region beyond the provided value) {-# INLINE normalTail #-} normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>               a -> RVar a@@ -82,13 +93,13 @@                 else return (r - x)  -- |Construct a 'Ziggurat' for sampling a normal distribution, given--- logBase 2 c, and the 'zGetIU' implementation.+-- @logBase 2 c@ and the 'zGetIU' implementation. normalZ ::   (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+-- | Ziggurat target function (upper half of a non-normalized gaussian PDF) normalF :: (Floating a, Ord a) => a -> a normalF x     | x <= 0    = 1@@ -105,6 +116,22 @@ normalFVol :: Floating a => a normalFVol = sqrt (0.5 * pi) +-- |A random variable sampling from the standard normal distribution+-- over any 'RealFloat' type (subject to the rest of the constraints -+-- it builds and uses a 'Ziggurat' internally, which requires the 'Erf'+-- and 'Storable' classes).  +-- +-- Because it computes a 'Ziggurat', it is very expensive to use for+-- just one evaluation, or even for multiple evaluations if not used and+-- reused monomorphically (to enable the ziggurat table to be let-floated+-- out).  If you don't know whether your use case fits this description+-- then you're probably better off using a different algorithm, such as+-- 'boxMullerNormalPair' or 'knuthPolarNormalPair'.  And of course if+-- you don't need the full generality of this definition then you're much+-- better off using 'doubleStdNormal' or 'floatStdNormal'.+--+-- As far as I know, this should be safe to use in any monomorphic+-- @Distribution Normal@ instance declaration. realFloatStdNormal :: (RealFloat a, Erf a, Storable a, Distribution Uniform a) => RVar a realFloatStdNormal = runZiggurat (normalZ p getIU)     where @@ -115,6 +142,8 @@             u <- uniform (-1) 1             return (fromIntegral i .&. (2^p-1), u) +-- |A random variable sampling from the standard normal distribution+-- over the 'Double' type. doubleStdNormal :: RVar Double doubleStdNormal = runZiggurat doubleStdNormalZ @@ -137,6 +166,8 @@             let (u,i) = wordToDoubleWithExcess w             return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1) +-- |A random variable sampling from the standard normal distribution+-- over the 'Float' type. floatStdNormal :: RVar Float floatStdNormal = runZiggurat floatStdNormalZ @@ -165,8 +196,11 @@ normalCdf :: (Real a) => a -> a -> a -> Double normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s) +-- |A specification of a normal distribution over the type 'a'. data Normal a+    -- |The \"standard\" normal distribution - mean 0, stddev 1     = StdNormal+    -- |@Normal m s@ is a normal distribution with mean @m@ and stddev @s@.     | Normal a a -- mean, sd  instance Distribution Normal Double where@@ -189,8 +223,10 @@  {-# SPECIALIZE stdNormal :: RVar Double #-} {-# SPECIALIZE stdNormal :: RVar Float #-}+-- |'stdNormal' is a normal variable with distribution 'StdNormal'. stdNormal :: Distribution Normal a => RVar a stdNormal = rvar StdNormal +-- |@normal m s@ is a random variable with distribution @'Normal' m s@. normal :: Distribution Normal a => a -> a -> RVar a normal m s = rvar (Normal m s)

@@ -17,14 +17,11 @@ floatingRayleigh s = do     u <- stdUniformPos     return (s * sqrt (-2 * log u))-    -    where -        stdUniformPos = do-            u <- stdUniform-            if u == 0-                then stdUniformPos-                else return u +-- |The rayleigh distribution with a specified mode (\"sigma\") parameter.+-- Its mean will be @sigma*sqrt(pi/2)@ and its variance will be @sigma^2*(4-pi)/2@+-- +-- (therefore if you want one with a particular mean @m@, @sigma@ should be @m*sqrt(2/pi)@) newtype Rayleigh a = Rayleigh a  rayleigh :: Distribution Rayleigh a => a -> RVar a

@@ -13,12 +13,23 @@ import Data.Random.Distribution import Data.Random.Distribution.Uniform -data Triangular a = Triangular-    { triLower  :: a-    , triMid    :: a-    , triUpper  :: a-    } deriving (Eq, Show)+-- |A description of a triangular distribution - a distribution whose PDF+-- is a triangle ramping up from a lower bound to a specified midpoint+-- and back down to the upper bound.  This is a very simple distribution+-- that does not generally occur naturally but is used sometimes as an+-- estimate of a true distribution when only the range of the values and+-- an approximate mode of the true distribution are known.+data Triangular a = Triangular {+    -- |The lower bound of the triangle in the PDF (the smallest number the distribution can generate)+    triLower  :: a,+    -- |The midpoint of the triangle (also the mode of the distribution)+    triMid    :: a,+    -- |The upper bound of the triangle (and the largest number the distribution can generate)+    triUpper  :: a}+    deriving (Eq, Show) +-- |Compute a triangular distribution for a 'Fractional' type.  The name is+-- a historical accident and may change in the future. realFloatTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a realFloatTriangular a b c     | a <= b && b <= c@@ -33,6 +44,7 @@ --        x <- stdUniform         return (b - ((1 - sqrt x) * (b-d))) +-- |@realFloatTriangularCDF a b c@ is the CDF of @realFloatTriangular a b c@. realFloatTriangularCDF :: RealFrac a => a -> a -> a -> a -> Double realFloatTriangularCDF a b c x     | x < a

@@ -14,6 +14,8 @@ 	     , StdUniform(..)     , stdUniform+    , stdUniformNonneg+    , stdUniformPos          , integralUniform     , realFloatUniform@@ -47,6 +49,7 @@  import Control.Monad.Loops +-- |Compute a random 'Integral' value between the 2 values provided (inclusive). integralUniform :: (Integral a) => a -> a -> RVar a integralUniform a b     | a > b     = compute b a@@ -71,26 +74,33 @@     Just x -> x powersOf256 = iterate (256 *) 1 +-- |Compute a random value for a 'Bounded' type, between 'minBound' and 'maxBound'+-- (inclusive for 'Integral' or 'Enum' types, in ['minBound', 'maxBound') for Fractional types.) boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a boundedStdUniform = uniform minBound maxBound  boundedStdUniformCDF :: (CDF Uniform a, Bounded a) => a -> Double boundedStdUniformCDF = cdf (Uniform minBound maxBound) +-- |Compute a random value for a 'Bounded' 'Enum' type, between 'minBound' and+-- 'maxBound' (inclusive) boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a boundedEnumStdUniform = enumUniform minBound maxBound  boundedEnumStdUniformCDF :: (Enum a, Bounded a, Ord a) => a -> Double boundedEnumStdUniformCDF = enumUniformCDF minBound maxBound +-- |Compute a uniform random 'Float' value in the range [0,1) floatStdUniform :: RVar Float floatStdUniform = do     x <- getRandomWord     return (wordToFloat x) +-- |Compute a uniform random 'Double' value in the range [0,1) doubleStdUniform :: RVar Double doubleStdUniform = getRandomDouble +-- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type  realFloatStdUniform :: RealFloat a => RVar a realFloatStdUniform = do     let (b, e) = decodeFloat one@@ -102,6 +112,8 @@          where one = 1 +-- |Compute a uniform random 'Fixed' value in the range [0,1), with any+-- desired precision. fixedStdUniform :: HasResolution r => RVar (Fixed r) fixedStdUniform = x     where@@ -110,42 +122,52 @@             u <- uniform 0 (res)             return (mkFixed u) +-- |The CDF of the random variable 'realFloatStdUniform'. realStdUniformCDF :: Real a => a -> Double realStdUniformCDF x     | x <= 0    = 0     | x >= 1    = 1     | otherwise = realToFrac x +-- |@floatUniform a b@ computes a uniform random 'Float' value in the range [a,b) floatUniform :: Float -> Float -> RVar Float floatUniform 0 1 = floatStdUniform floatUniform a b = do     x <- floatStdUniform     return (a + x * (b - a)) +-- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b) doubleUniform :: Double -> Double -> RVar Double doubleUniform 0 1 = doubleStdUniform doubleUniform a b = do     x <- doubleStdUniform     return (a + x * (b - a)) +-- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for+-- any 'RealFloat' type realFloatUniform :: RealFloat a => a -> a -> RVar a realFloatUniform 0 1 = realFloatStdUniform realFloatUniform a b = do     x <- realFloatStdUniform     return (a + x * (b - a)) +-- |@fixedUniform a b@ computes a uniform random 'Fixed' value in the range +-- [a,b), with any desired precision. fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVar (Fixed r) fixedUniform a b = do     u <- integralUniform (unMkFixed a) (unMkFixed b)     return (mkFixed u) -realUniformCDF :: Real a => a -> a -> a -> Double+-- |@realUniformCDF a b@ is the CDF of the random variable @realFloatUniform a b@.+realUniformCDF :: RealFrac a => a -> a -> a -> Double realUniformCDF a b x     | b < a     = realUniformCDF b a x     | x <= a    = 0     | x >= b    = 1-    | otherwise = realToFrac (x-a) / realToFrac (b-a)+    | otherwise = realToFrac ((x-a) / (b-a)) +-- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for+-- any 'Enum' type enumUniform :: Enum a => a -> a -> RVar a enumUniform a b = do     x <- integralUniform (fromEnum a) (fromEnum b)@@ -160,6 +182,9 @@          where e2f = fromIntegral . fromEnum +-- @uniform a b@ computes a uniformly distributed random value in the range+-- [a,b] for 'Integral' or 'Enum' types and in the range [a,b) for 'Fractional'+-- types.  Requires a @Distribution Uniform@ instance for the type. uniform :: Distribution Uniform a => a -> a -> RVar a uniform a b = rvar (Uniform a b) @@ -170,14 +195,33 @@ stdUniform :: (Distribution StdUniform a) => RVar a stdUniform = rvar StdUniform -stdUniformPos :: (Distribution StdUniform a, Ord a, Num a) => RVar a+-- |Like 'stdUniform', but uses 'abs' to return only positive or zero values.+stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVar a+stdUniformNonneg = abs `fmap` stdUniform++-- |Like 'stdUniform' but only returns positive values.+stdUniformPos :: (Distribution StdUniform a, Num a) => RVar a stdUniformPos = do-    x <- stdUniform-    if x > 0+    x <- stdUniformNonneg+    if x /= 0         then return x         else stdUniformPos -data Uniform t = Uniform !t !t+-- |A definition of a uniform distribution over the type @t@.  See also 'uniform'.+data Uniform t = +    -- |A uniform distribution defined by a lower and upper range bound.+    -- For 'Integral' and 'Enum' types, the range is inclusive.  For 'Fractional'+    -- types the range includes the lower bound but not the upper.+    Uniform !t !t++-- |A name for the \"standard\" uniform distribution over the type @t@,+-- if one exists.  See also 'stdUniform'.+--+-- For 'Integral' and 'Enum' types that are also 'Bounded', this is+-- the uniform distribution over the full range of the type.+-- For un-'Bounded' 'Integral' types this is not defined.+-- For 'Fractional' types this is a random variable in the range [0,1)+-- (that is, 0 to 1 including 0 but not including 1). data StdUniform t = StdUniform  $( replicateInstances ''Int integralTypes [d|

@@ -1,24 +1,22 @@-{-- -      ``Data/Random/Distribution/Ziggurat''- -      A generic "ziggurat algorithm" implementation.  Fairly rough right- -      now.- -      - -      There is a lot of room for improvement in 'findBin0' especially.- -      It needs a fair amount of cleanup and elimination of redundant- -      calculation, as well as either a justification for using the simple- -      'findMinFrom' or a proper root-finding algorithm. - -      - -      It would also be nice to add (preferably via its own library)- -      support for numerical integration and differentiation, so that- -      tables can be derived from only a PDF (if the end user is- -      willing to take the performance hit for the convenience).- -} {-# LANGUAGE         MultiParamTypeClasses,         FlexibleInstances, FlexibleContexts,         RecordWildCards   #-} +-- |A generic \"ziggurat algorithm\" implementation.  Fairly rough right+--  now.+--  +--  There is a lot of room for improvement in 'findBin0' especially.+--  It needs a fair amount of cleanup and elimination of redundant+--  calculation, as well as either a justification for using the simple+--  'findMinFrom' or a proper root-finding algorithm. +--  +--  It would also be nice to add (preferably by pulling in an +--  external package) support for numerical integration and +--  differentiation, so that tables can be derived from only a +--  PDF (if the end user is willing to take the performance and +--  accuracy hit for the convenience). module Data.Random.Distribution.Ziggurat     ( Ziggurat(..)     , mkZigguratRec@@ -38,54 +36,129 @@  vec ! i = index vec i -data Ziggurat t = Ziggurat-    { zTable_xs         :: Vector t-    , zTable_x_ratios   :: Vector t-    , zTable_ys         :: Vector t-    , zGetIU            :: RVar (Int, t)-    , zTailDist         :: RVar t-    , zUniform          :: t -> t -> RVar t-    , zFunc             :: t -> t-    , zMirror           :: Bool+-- |A data structure containing all the data that is needed+-- to implement Marsaglia & Tang's \"ziggurat\" algorithm for+-- sampling certain kinds of random distributions.+--+-- The documentation here is probably not sufficient to tell a user exactly+-- how to build one of these from scratch, but it is not really intended to+-- be.  There are several helper functions that will build 'Ziggurat's.+-- The pathologically curious may wish to read the 'runZiggurat' source.+-- That is the ultimate specification of the semantics of all these fields.+data Ziggurat t = Ziggurat {+        -- |The X locations of each bin in the distribution.  Bin 0 is the+        -- 'infinite' one.+        -- +        -- In the case of bin 0, the value given is sort of magical - x[0] is+        -- defined to be V/f(R).  It's not actually the location of any bin, +        -- but a value computed to make the algorithm more concise and slightly +        -- faster by not needing to specially-handle bin 0 quite as often.+        -- If you really need to know why it works, see the 'runZiggurat'+        -- source or \"the literature\" - it's a fairly standard setup.+        zTable_xs         :: Vector t,+        -- |The ratio of each bin's Y value to the next bin's Y value+        zTable_x_ratios   :: Vector t,+        -- |The Y value (zFunc x) of each bin+        zTable_ys         :: Vector t,+        -- |An RVar providing a random tuple consisting of:+        --+        --  * a bin index, uniform over [0,c) :: Int (where @c@ is the+        --    number of bins in the tables)+        -- +        --  * a uniformly distributed fractional value, from -1 to 1 +        --    if not mirrored, from 0 to 1 otherwise.+        --+        -- This is provided as a single 'RVar' because it can be implemented+        -- more efficiently than naively sampling 2 separate values - a+        -- single random word (64 bits) can be efficiently converted to+        -- a double (using 52 bits) and a bin number (using up to 12 bits),+        -- for example.+        zGetIU            :: RVar (Int, t),+        +        -- |The distribution for the final \"virtual\" bin+        -- (the ziggurat algorithm does not handle distributions+        -- that wander off to infinity, so another distribution is needed+        -- to handle the last \"bin\" that stretches to infinity)+        zTailDist         :: RVar t,+        +        -- |A copy of the uniform RVar generator for the base type,+        -- so that @Distribution Uniform t@ is not needed when sampling+        -- from a Ziggurat (makes it a bit more self-contained).+        zUniform          :: t -> t -> RVar t,+        +        -- |The (one-sided antitone) PDF, not necessarily normalized+        zFunc             :: t -> t,+        +        -- |A flag indicating whether the distribution should be+        -- mirrored about the origin (the ziggurat algorithm it+        -- its native form only samples from one-sided distributions.+        -- By mirroring, we can extend it to symmetric distributions+        -- such as the normal distribution)+        zMirror           :: Bool     } +-- |Sample from the distribution encoded in a 'Ziggurat' data structure. {-# INLINE runZiggurat #-} runZiggurat :: (Num a, Ord a, Storable a) =>                Ziggurat a -> RVar a runZiggurat Ziggurat{..} = go     where         go = do+            -- Select a bin (I) and a uniform value (U) from -1 to 1+            -- (or 0 to 1 if not mirroring the distribution).+            -- Let X be U scaled to the size of the selected bin.             (i,u) <- zGetIU             let x = u * zTable_xs ! i             +            -- if the uniform value U falls in the area "clearly inside" the+            -- bin, accept X immediately.+            -- Otherwise, depending on the bin selected, use either the+            -- tail distribution or an accept/reject test.             if abs u < zTable_x_ratios ! i                 then return $! x                 else if i == 0-                    then if x < 0 then fmap negate zTailDist else zTailDist-                    else do-                        v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)-                        if v < zFunc (abs x)-                            then return $! x-                            else go+                    then sampleTail x+                    else sampleGreyArea i x+        +        -- when the sample falls in the "grey area" (the area between+        -- the Y values of the selected bin and the bin after that one),+        -- use an accept/reject method based on the target PDF.+        sampleGreyArea i x = do+            v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)+            if v < zFunc (abs x)+                then return $! x+                else go+        +        -- if the selected bin is the "infinite" one, call it quits and+        -- defer to the tail distribution (mirroring if needed to ensure+        -- the result has the sign already selected by zGetIU)+        sampleTail x+            | x < 0     = fmap negate zTailDist+            | otherwise = zTailDist  --- |Build the tables to implement the "ziggurat algorithm" devised by +-- |Build the tables to implement the \"ziggurat algorithm\" devised by  -- Marsaglia & Tang, attempting to automatically compute the R and V -- values. --  -- Arguments: --  --  * flag indicating whether to mirror the distribution---  * the (one-sided antitone) CDF---  * the inverse of the CDF+-- +--  * the (one-sided antitone) PDF, not necessarily normalized+-- +--  * the inverse of the PDF+--  --  * the number of bins+--  --  * R, the x value of the first bin+--  --  * V, the volume of each bin---  * an RVar providing a random tuple consisting of:---          - a bin index, uniform over [0,c) :: Int---          - a uniformly distributed fractional value, from -1 to 1 if not mirrored, from 0 to 1 otherwise.+-- +--  * an RVar providing the 'zGetIU' random tuple+--  --  * an RVar sampling from the tail (the region where x > R)-+--  mkZiggurat_ :: (RealFloat t, Storable t,                Distribution Uniform t) =>               Bool@@ -99,21 +172,21 @@               -> Ziggurat t mkZiggurat_ m f fInv c r v getIU tailDist = z     where z = Ziggurat-            { zTable_xs            = zigguratTable f fInv c r v-            , zTable_x_ratios    = precomputeRatios (zTable_xs z)-            , zTable_ys      = Vec.map f (zTable_xs z)+            { zTable_xs         = zigguratTable f fInv c r v+            , zTable_x_ratios   = precomputeRatios (zTable_xs z)+            , zTable_ys         = Vec.map f (zTable_xs z)             , zGetIU            = getIU-            , zUniform = uniform-            , zFunc = f-            , zTailDist = tailDist-            , zMirror = m+            , zUniform          = uniform+            , zFunc             = f+            , zTailDist         = tailDist+            , zMirror           = m             } --- |Build the tables to implement the "ziggurat algorithm" devised by +-- |Build the tables to implement the \"ziggurat algorithm\" devised by  -- Marsaglia & Tang, attempting to automatically compute the R and V -- values. -- --- Arguments are the same as for |mkZigguratRec|, with an additional+-- Arguments are the same as for 'mkZigguratRec', with an additional -- argument for the tail distribution as a function of the selected -- R value. mkZiggurat :: (RealFloat t, Storable t,@@ -138,15 +211,31 @@ -- Arguments: --  --  * flag indicating whether to mirror the distribution---  * the (one-sided antitone) CDF---  * the inverse of the CDF---  * the integral of the CDF (definite, from 0)---  * the estimated volume under the CDF (from 0 to +infinity)---  * the chunk size (number of bins).  64 seems to perform well in practice.---  * an RVar providing a random tuple consisting of:---          - a bin index, uniform over [0,c) :: Int---          - a uniformly distributed fractional value, from -1 to 1 if not mirrored, from 0 to 1 otherwise.-+--+--  * the (one-sided antitone) PDF, not necessarily normalized+--+--  * the inverse of the PDF+--+--  * the integral of the PDF (definite, from 0)+--+--  * the estimated volume under the PDF (from 0 to +infinity)+--+--  * the chunk size (number of bins in each layer).  64 seems to+--    perform well in practice.+--+--  * an RVar providing the 'zGetIU' random tuple+--+mkZigguratRec ::+  (RealFloat t, Storable t,+   Distribution Uniform t) =>+  Bool+  -> (t -> t)+  -> (t -> t)+  -> (t -> t)+  -> t+  -> Int+  -> RVar (Int, t)+  -> Ziggurat t mkZigguratRec m f fInv fInt fVol c getIU =      mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU))         where@@ -201,12 +290,18 @@ -- Arguments: --  --  * Number of bins---  * function (one-sided antitone CDF, not necessarily normalized)+--+--  * target function (one-sided antitone PDF, not necessarily normalized)+-- --  * function inverse+-- --  * function definite integral (from 0 to _)+-- --  * estimate of total volume under function (integral from 0 to infinity) -- -- Result: (R,V)+findBin0 :: (RealFloat b) => +    Int -> (b -> b) -> (b -> b) -> (b -> b) -> b -> (b, b) findBin0 cInt f fInv fInt fVol = (r,v r)     where         c = fromIntegral cInt

@@ -1,7 +1,7 @@ {-  -      ``Data/Random/Internal/Find''  -  Utilities for searching fractional domains.  Needs cleanup, testing,- -  and such.+ -  and such.  Used for constructing generic ziggurats.  -}  module Data.Random.Internal.Find where

@@ -1,9 +1,24 @@+{-# LANGUAGE CPP #-} module Data.Random.Internal.Fixed where  import Data.Fixed import Unsafe.Coerce +#ifdef base_4_2+-- So much for backward compatibility through base (>=5) ...+ resolutionOf :: HasResolution r => f r -> Integer+resolutionOf = resolution++resolutionOf2 :: HasResolution r => f (g r) -> Integer+resolutionOf2 x = resolution (res x)+    where+        res :: HasResolution r => f (g r) -> g r+        res = undefined++#else++resolutionOf :: HasResolution r => f r -> Integer resolutionOf x = resolution (res x)     where         res :: HasResolution r => f r -> r@@ -15,6 +30,13 @@         res :: HasResolution r => f (g r) -> r         res = undefined +#endif++-- |The 'Fixed' type doesn't expose its constructors, but I need a way to+-- convert them to and from their raw representation in order to sample+-- them.  As long as 'Fixed' is a newtype wrapping 'Integer', 'mkFixed' and+-- 'unMkFixed' as defined here will work.  Both are implemented using +-- 'unsafeCoerce'. mkFixed :: Integer -> Fixed r mkFixed = unsafeCoerce 

@@ -5,7 +5,21 @@         TemplateHaskell   #-} -module Data.Random.Internal.TH where+-- |Template Haskell utility code to replicate instance declarations+-- to cover large numbers of types.  I'm doing that rather than using+-- class contexts because most Distribution instances need to cover+-- multiple classes (such as Enum, Integral and Fractional) and that+-- can't be done easily because of overlap.  +-- +-- I experimented a bit with a convoluted type-level classification +-- scheme, but I think this is simpler and easier to understand.  It +-- makes the haddock docs more cluttered because of the combinatorial +-- explosion of instances, but overall I think it's just more sane than +-- anything else I've come up with yet.+module Data.Random.Internal.TH+    ( replicateInstances+    , integralTypes, realFloatTypes+    ) where  import Data.Generics import Language.Haskell.TH@@ -13,20 +27,48 @@ import Data.Word import Data.Int +-- |Names of standard 'Integral' types+integralTypes :: [Name] integralTypes =      [ ''Int,   ''Integer     , ''Int8,  ''Int16,  ''Int32,  ''Int64     , ''Word8, ''Word16, ''Word32, ''Word64     ] +-- |Names of standard 'RealFloat' types+realFloatTypes :: [Name] realFloatTypes =     [ ''Float, ''Double ] +-- @replaceName x y@ is a function that will+-- replace @x@ with @y@ whenever it sees it.  That is:+--+-- > replaceName x y x  ==>  y+-- > replaceName x y z  ==>  z+--  (@z /= x@) replaceName :: Name -> Name -> Name -> Name replaceName x y z     | x == z    = y     | otherwise = z +-- | @replicateInstances standin types decls@ will take the template-haskell+-- 'Dec's in @decls@ and substitute every instance of the 'Name' @standin@ with+-- each 'Name' in @types@, producing one copy of the 'Dec's in @decls@ for every+-- 'Name' in @types@.+-- +-- For example, 'Data.Random.Distribution.Uniform' has the following bit of TH code:+-- +-- @ $( replicateInstances ''Int integralTypes [d|                                                  @+-- +-- @       instance Distribution Uniform Int   where rvar (Uniform a b) = integralUniform a b       @+-- +-- @       instance CDF Uniform Int            where cdf  (Uniform a b) = integralUniformCDF a b    @+-- +-- @   |])                                                                                          @+-- +-- This code takes those 2 instance declarations and creates identical ones for+-- every type named in 'integralTypes'.+replicateInstances :: (Monad m, Data t) => Name -> [Name] -> m [t] -> m [t] replicateInstances standin types decls = do     decls <- decls     sequence

@@ -12,6 +12,11 @@ import Data.Word  {-# INLINE buildWord #-}+-- |Build a word out of 8 bytes.  No promises are made regarding the order+-- in which the bytes are stuffed.  Note that this means that a 'RandomSource'+-- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,+-- may return different random values on different platforms when started +-- with the same seed, depending on the platform's endianness. buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64 buildWord b0 b1 b2 b3 b4 b5 b6 b7     = unsafePerformIO . allocaBytes 8 $ \p -> do@@ -26,7 +31,7 @@         peek (castPtr p)  {-# INLINE wordToFloat #-}--- |Pack 23 unspecified bits from a 'Word64' into a 'Float' in the range [0,1).+-- |Pack the low 23 bits from a 'Word64' into a 'Float' in the range [0,1). -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'. wordToFloat :: Word64 -> Float wordToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)@@ -38,7 +43,7 @@ wordToFloatWithExcess x = (wordToFloat x, x `shiftR` 23)  {-# INLINE wordToDouble #-}--- |Pack 52 unspecified bits from a 'Word64' into a 'Double' in the range [0,1).+-- |Pack the low 52 bits from a 'Word64' into a 'Double' in the range [0,1). -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'. wordToDouble :: Word64 -> Double wordToDouble x = (encodeFloat $! toInteger (x .&. 0x000fffffffffffff {- 2^52-1 -})) $ (-52)

@@ -6,12 +6,17 @@ import qualified System.Random.Shuffle as SRS import Control.Monad +-- | A random variable returning an arbitrary element of the given list.+-- Every element has equal probability of being chosen.  Because it is a+-- pure 'RVar' it has no memory - that is, it \"draws with replacement.\" randomElement :: [a] -> RVar a randomElement [] = error "randomElement: empty list!" randomElement xs = do     n <- uniform 0 (length xs - 1)     return (xs !! n) +-- | A random variable that returns the given list in an arbitrary shuffled+-- order.  Every ordering of the list has equal probability. shuffle :: [a] -> RVar [a] shuffle [] = return [] shuffle xs = do@@ -19,9 +24,16 @@          return (SRS.shuffle xs (reverse is)) +-- | A random variable that shuffles a list of a known length. Useful for +-- shuffling large lists when the length is known in advance.+-- Avoids needing to traverse the list to discover its length.  Each ordering+-- has equal probability.+--+-- Throws an error the list is not exactly as long as claimed. shuffleN :: Int -> [a] -> RVar [a] shuffleN 0 xs = return []-shuffleN (n+1) xs = do+shuffleN m@(n+1) xs = do     is <- sequence [uniform 0 i | i <- [n,n-1..1]]     return (SRS.shuffle xs is)+     

@@ -81,6 +81,17 @@     getRandomWord   = RVarT $ \k (RVarDict s) -> getRandomWordFrom   s >>= \a -> k a     getRandomDouble = RVarT $ \k (RVarDict s) -> getRandomDoubleFrom s >>= \a -> k a +-- I would really like to be able to do this, but I can't because of the+-- blasted Eq and Show in Num's class context...+-- instance (Applicative m, Num a) => Num (RVarT m a) where+--     (+) = liftA2 (+)+--     (-) = liftA2 (-)+--     (*) = liftA2 (*)+--     negate = liftA negate+--     signum = liftA signum+--     abs = liftA abs+--     fromInteger = pure . fromInteger+ -- some 'fundamental' RVarTs -- this maybe ought to even be a part of the RandomSource class... {-# INLINE nByteInteger #-}

@@ -25,7 +25,7 @@ instance Distribution d t => Sampleable d m t where     sampleFrom src d = runRVarT (rvar d) src --- This instance conflicts with the other, but because RVarT is not a Distribution there is no conflict.+-- This instance overlaps with the other, but because RVarT is not a Distribution there is no conflict. instance Lift m n => Sampleable (RVarT m) n t where     sampleFrom src x = runRVarT x src 

@@ -1,12 +1,14 @@-{-- -      ``Data/Random/Source/PureMT''- -} {-# LANGUAGE     MultiParamTypeClasses,     FlexibleContexts, FlexibleInstances,     UndecidableInstances   #-} +-- |This module provides functions useful for implementing new 'MonadRandom'+-- and 'RandomSource' instances for state-abstractions containing 'PureMT'+-- values (the pure pseudorandom generator provided by the+-- mersenne-random-pure64 package), as well as instances for some common+-- cases. module Data.Random.Source.PureMT where  import Data.Random.Internal.Words

@@ -5,6 +5,11 @@     MultiParamTypeClasses, FlexibleInstances, UndecidableInstances   #-} +-- |This module provides functions useful for implementing new 'MonadRandom'+-- and 'RandomSource' instances for state-abstractions containing 'StdGen'+-- values (the pure pseudorandom generator provided by the System.Random+-- module in the \"random\" package), as well as instances for some common+-- cases. module Data.Random.Source.StdGen where  import Data.Random.Internal.Words