random-fu 0.2.7.7 → 0.3.0.0
raw patch · 25 files changed
+539/−459 lines, 25 filesdep −random-sourcedep ~basedep ~rvar
Dependencies removed: random-source
Dependency ranges changed: base, rvar
Files
- changelog.md +4/−0
- random-fu.cabal +11/−13
- src/Data/Random.hs +25/−30
- src/Data/Random/Distribution.hs +20/−20
- src/Data/Random/Distribution/Bernoulli.hs +58/−27
- src/Data/Random/Distribution/Beta.hs +8/−8
- src/Data/Random/Distribution/Binomial.hs +95/−30
- src/Data/Random/Distribution/Categorical.hs +26/−28
- src/Data/Random/Distribution/Dirichlet.hs +1/−1
- src/Data/Random/Distribution/Gamma.hs +6/−7
- src/Data/Random/Distribution/Multinomial.hs +2/−2
- src/Data/Random/Distribution/Normal.hs +55/−26
- src/Data/Random/Distribution/Poisson.hs +82/−24
- src/Data/Random/Distribution/Rayleigh.hs +2/−2
- src/Data/Random/Distribution/Triangular.hs +1/−1
- src/Data/Random/Distribution/Uniform.hs +54/−74
- src/Data/Random/Distribution/Weibull.hs +1/−1
- src/Data/Random/Distribution/Ziggurat.hs +50/−50
- src/Data/Random/Internal/Find.hs +8/−8
- src/Data/Random/Internal/Fixed.hs +1/−1
- src/Data/Random/Internal/TH.hs +0/−79
- src/Data/Random/Lift.hs +4/−6
- src/Data/Random/List.hs +2/−3
- src/Data/Random/RVar.hs +4/−3
- src/Data/Random/Sample.hs +19/−15
changelog.md view
@@ -1,3 +1,7 @@+* Chnages in 0.3.0.0:++ * Drop usage of `random-source` in favor of `random`+ * Changes in 0.2.7.7: Update to random-1.2. Revert 0.2.7.6 changes (which added an extra constraint to `Data.Random.Sample.sampleState` and `Data.Random.Sample.sampleStateT`). * Changes in 0.2.7.4: Compatibility with ghc 8.8.
random-fu.cabal view
@@ -1,5 +1,5 @@ name: random-fu-version: 0.2.7.7+version: 0.3.0.0 stability: provisional cabal-version: >= 1.10@@ -12,30 +12,30 @@ category: Math synopsis: Random number generation-description: Random number generation based on modeling random +description: Random number generation based on modeling random variables in two complementary ways: first, by the parameters of standard mathematical distributions and, second, by an abstract type ('RVar') which can be composed and manipulated monadically and sampled in either monadic or \"pure\" styles. .- The primary purpose of this library is to support + The primary purpose of this library is to support defining and sampling a wide variety of high quality random variables. Quality is prioritized over speed, but performance is an important goal too. .- In my testing, I have found it capable of speed + In my testing, I have found it capable of speed comparable to other Haskell libraries, but still- a fair bit slower than straight C implementations of + a fair bit slower than straight C implementations of the same algorithms. -tested-with: GHC == 7.10.3+tested-with: GHC == 8.10.7 extra-source-files: changelog.md source-repository head type: git- location: https://github.com/mokus0/random-fu.git+ location: https://github.com/haskell-numerics/random-fu subdir: random-fu Flag base4_2@@ -72,7 +72,6 @@ Data.Random.Distribution.Ziggurat Data.Random.Internal.Find Data.Random.Internal.Fixed- Data.Random.Internal.TH Data.Random.Lift Data.Random.List Data.Random.RVar@@ -83,25 +82,24 @@ else cpp-options: -Dold_Fixed build-depends: base >= 4 && <4.2- + if flag(mtl2) build-depends: mtl == 2.* cpp-options: -DMTL2 else build-depends: mtl == 1.*- + build-depends: math-functions, monad-loops >= 0.3.0.1, random >= 1.2 && < 1.3, random-shuffle,- random-source == 0.3.*,- rvar == 0.2.*,+ rvar >= 0.3, syb, template-haskell, transformers, vector >= 0.7, erf- + if impl(ghc == 7.2.1) -- Doesn't work under GHC 7.2.1 due to -- http://hackage.haskell.org/trac/ghc/ticket/5410
src/Data/Random.hs view
@@ -1,39 +1,39 @@ -- |Flexible modeling and sampling of random variables. ----- The central abstraction in this library is the concept of a random --- variable. It is not fully formalized in the standard measure-theoretic --- language, but rather is informally defined as a \"thing you can get random --- values out of\". Different random variables may have different types of +-- The central abstraction in this library is the concept of a random+-- variable. It is not fully formalized in the standard measure-theoretic+-- language, but rather is informally defined as a \"thing you can get random+-- values out of\". Different random variables may have different types of -- values they can return or the same types but different probabilities for -- each value they can return. The random values you get out of them are -- traditionally called \"random variates\".--- --- Most imperative-language random number libraries are all about obtaining --- and manipulating random variates. This one is about defining, manipulating --- and sampling random variables. Computationally, the distinction is small --- and mostly just a matter of perspective, but from a program design +--+-- Most imperative-language random number libraries are all about obtaining+-- and manipulating random variates. This one is about defining, manipulating+-- and sampling random variables. Computationally, the distinction is small+-- and mostly just a matter of perspective, but from a program design -- perspective it provides both a powerfully composable abstraction and a -- very useful separation of concerns.--- +-- -- Abstract random variables as implemented by 'RVar' are composable. They can -- be defined in a monadic / \"imperative\" style that amounts to manipulating -- variates, but with strict type-level isolation. Concrete random variables -- are also provided, but they do not compose as generically. The 'Distribution'--- type class allows concrete random variables to \"forget\" their concreteness --- so that they can be composed. For examples of both, see the documentation --- for 'RVar' and 'Distribution', as well as the code for any of the concrete +-- type class allows concrete random variables to \"forget\" their concreteness+-- so that they can be composed. For examples of both, see the documentation+-- for 'RVar' and 'Distribution', as well as the code for any of the concrete -- distributions such as 'Uniform', 'Gamma', etc.--- +-- -- Both abstract and concrete random variables can be sampled (despite the -- types GHCi may list for the functions) by the functions in "Data.Random.Sample".--- +-- -- Random variable sampling is done with regard to a generic basis of primitive--- random variables defined in "Data.Random.Internal.Primitives". This basis +-- random variables defined in "Data.Random.Internal.Primitives". This basis -- is very low-level and the actual set of primitives is still fairly experimental, -- which is why it is in the \"Internal\" sub-heirarchy. User-defined variables -- should use the existing high-level variables such as 'Uniform' and 'Normal' -- rather than these basis variables. "Data.Random.Source" defines classes for--- entropy sources that provide implementations of these primitive variables. +-- entropy sources that provide implementations of these primitive variables. -- Several implementations are available in the Data.Random.Source.* modules. module Data.Random ( -- * Random variables@@ -43,32 +43,26 @@ -- ** Concrete ('Distribution') Distribution(..), CDF(..), PDF(..),- + -- * Sampling random variables- Sampleable(..), sample, sampleState, sampleStateT,- + Sampleable(..), sample, sampleState, samplePure,+ -- * A few very common distributions Uniform(..), uniform, uniformT, StdUniform(..), stdUniform, stdUniformT, Normal(..), normal, stdNormal, normalT, stdNormalT, Gamma(..), gamma, gammaT,- + -- * Entropy Sources- MonadRandom, RandomSource, StdRandom(..),- + StatefulGen, RandomGen,+ -- * Useful list-based operations randomElement, shuffle, shuffleN, shuffleNofM- + ) where import Data.Random.Sample-import Data.Random.Source (MonadRandom, RandomSource)-import Data.Random.Source.IO ()-import Data.Random.Source.MWC ()-import Data.Random.Source.StdGen ()-import Data.Random.Source.PureMT ()-import Data.Random.Source.Std import Data.Random.Distribution import Data.Random.Distribution.Gamma import Data.Random.Distribution.Normal@@ -78,3 +72,4 @@ import Data.Random.List import Data.Random.RVar +import System.Random.Stateful (StatefulGen, RandomGen)
src/Data/Random/Distribution.hs view
@@ -13,7 +13,7 @@ -- > data Normal a -- > = StdNormal -- > | Normal a a--- +-- -- Where the two parameters of the 'Normal' data constructor are the mean and -- standard deviation of the random variable, respectively. To make use of -- the 'Normal' type, one can convert it to an 'rvar' and manipulate it or@@ -21,39 +21,39 @@ -- -- > x <- sample (rvar (Normal 10 2)) -- > x <- sample (Normal 10 2)--- +-- -- A 'Distribution' is typically more transparent than an 'RVar'--- but less composable (precisely because of that transparency). There are +-- but less composable (precisely because of that transparency). There are -- several practical uses for types implementing 'Distribution':--- --- * Typically, a 'Distribution' will expose several parameters of a standard +--+-- * Typically, a 'Distribution' will expose several parameters of a standard -- mathematical model of a probability distribution, such as mean and std deviation for -- the normal distribution. Thus, they can be manipulated analytically using -- mathematical insights about the distributions they represent. For example, -- a collection of bernoulli variables could be simplified into a (hopefully) smaller -- collection of binomial variables.--- +-- -- * Because they are generally just containers for parameters, they can be--- easily serialized to persistent storage or read from user-supplied +-- easily serialized to persistent storage or read from user-supplied -- configurations (eg, initialization data for a simulation).--- +-- -- * If a type additionally implements the 'CDF' subclass, which extends -- 'Distribution' with a cumulative density function, an arbitrary random -- variable 'x' can be tested against the distribution by testing -- @fmap (cdf dist) x@ for uniformity.--- +-- -- On the other hand, most 'Distribution's will not be closed under all the -- same operations as 'RVar' (which, being a monad, has a fully turing-complete--- internal computational model). The sum of two uniformly-distributed --- variables, for example, is not uniformly distributed. To support general --- composition, the 'Distribution' class defines a function 'rvar' to --- construct the more-abstract and more-composable 'RVar' representation +-- internal computational model). The sum of two uniformly-distributed+-- variables, for example, is not uniformly distributed. To support general+-- composition, the 'Distribution' class defines a function 'rvar' to+-- construct the more-abstract and more-composable 'RVar' representation -- of a random variable. class Distribution d t where -- |Return a random variable with this distribution. rvar :: d t -> RVar t rvar = rvarT- + -- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'. -- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using -- either 'lift' or 'sample'.@@ -66,8 +66,8 @@ pdf d = exp . logPdf d logPdf :: d t -> t -> Double logPdf d = log . pdf d- + 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@@ -76,19 +76,19 @@ -- -- In the case where 't' is an instance of Ord, 'cdf' should correspond -- 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)@ -- must be uniformly distributed over (0,1). Inclusion of either endpoint is optional, -- though the preferred range is (0,1].- -- - -- Note that this definition requires that 'cdf' for a product type - -- should _not_ be a joint CDF as commonly defined, as that definition + --+ -- Note that this definition requires that '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. That is, -- it should represent the CDF with respect to the lexicographic order -- of the product.- -- + -- -- The present specification is probably only really useful for testing -- conformance of a variable to its target distribution, and I am open to -- suggestions for more-useful specifications (especially with regard to
src/Data/Random/Distribution/Bernoulli.hs view
@@ -1,22 +1,21 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,- UndecidableInstances,- TemplateHaskell+ UndecidableInstances #-} {-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-} module Data.Random.Distribution.Bernoulli where -import Data.Random.Internal.TH- import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Uniform import Data.Ratio import Data.Complex+import Data.Int+import Data.Word -- |Generate a Bernoulli variate with the given probability. For @Bool@ results, -- @bernoulli p@ will return True (p*100)% of the time and False otherwise.@@ -57,7 +56,7 @@ newtype Bernoulli b a = Bernoulli b -instance (Fractional b, Ord b, Distribution StdUniform b) +instance (Fractional b, Ord b, Distribution StdUniform b) => Distribution (Bernoulli b) Bool where rvarT (Bernoulli p) = boolBernoulli p@@ -66,34 +65,66 @@ where cdf (Bernoulli p) = boolBernoulliCDF p -$( replicateInstances ''Int integralTypes [d|- instance Distribution (Bernoulli b) Bool - => Distribution (Bernoulli b) Int- where- rvarT (Bernoulli p) = generalBernoulli 0 1 p- instance CDF (Bernoulli b) Bool- => CDF (Bernoulli b) Int- where- cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p- |] )+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Integer where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Integer where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Int where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int8 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Int8 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int16 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Int16 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int32 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Int32 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int64 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Int64 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Word where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word8 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Word8 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word16 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Word16 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word32 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Word32 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word64 where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Word64 where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p -$( replicateInstances ''Float realFloatTypes [d|- instance Distribution (Bernoulli b) Bool - => Distribution (Bernoulli b) Float- where- rvarT (Bernoulli p) = generalBernoulli 0 1 p- instance CDF (Bernoulli b) Bool- => CDF (Bernoulli b) Float- where- cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p- |] )+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Float where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Float where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Double where+ rvarT (Bernoulli p) = generalBernoulli 0 1 p+instance CDF (Bernoulli b) Bool => CDF (Bernoulli b) Double where+ cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p instance (Distribution (Bernoulli b) Bool, Integral a)- => Distribution (Bernoulli b) (Ratio a) + => Distribution (Bernoulli b) (Ratio a) where rvarT (Bernoulli p) = generalBernoulli 0 1 p instance (CDF (Bernoulli b) Bool, Integral a)- => CDF (Bernoulli b) (Ratio a) + => CDF (Bernoulli b) (Ratio a) where cdf (Bernoulli p) = generalBernoulliCDF (>=) 0 1 p instance (Distribution (Bernoulli b) Bool, RealFloat a)
src/Data/Random/Distribution/Beta.hs view
@@ -1,16 +1,13 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,- UndecidableInstances,- TemplateHaskell+ UndecidableInstances #-} {-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-} module Data.Random.Distribution.Beta where -import Data.Random.Internal.TH- import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Gamma@@ -57,7 +54,10 @@ where pdf (Beta a b) = realToFrac . exp . logBetaPdf (realToFrac a) (realToFrac b) . realToFrac -$( replicateInstances ''Float realFloatTypes [d|- instance Distribution Beta Float- where rvarT (Beta a b) = fractionalBeta a b- |])+instance Distribution Beta Float+ where+ rvarT (Beta a b) = fractionalBeta a b++instance Distribution Beta Double+ where+ rvarT (Beta a b) = fractionalBeta a b
src/Data/Random/Distribution/Binomial.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,- UndecidableInstances, TemplateHaskell,+ UndecidableInstances, BangPatterns #-} @@ -9,13 +9,14 @@ module Data.Random.Distribution.Binomial where -import Data.Random.Internal.TH- import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Beta import Data.Random.Distribution.Uniform +import Data.Int+import Data.Word+ import Numeric.SpecFunctions ( stirlingError ) import Numeric.SpecFunctions.Extra ( bd0 ) import Numeric ( log1p )@@ -131,31 +132,95 @@ data Binomial b a = Binomial a b -$( replicateInstances ''Int integralTypes [d|- instance ( Floating b, Ord b- , Distribution Beta b- , Distribution StdUniform b- ) => Distribution (Binomial b) Int- where- rvarT (Binomial t p) = integralBinomial t p- instance ( Real b , Distribution (Binomial b) Int- ) => CDF (Binomial b) Int- where cdf (Binomial t p) = integralBinomialCDF t p- instance ( Real b , Distribution (Binomial b) Int- ) => PDF (Binomial b) Int- where pdf (Binomial t p) = integralBinomialPDF t p- logPdf (Binomial t p) = integralBinomialLogPdf t p- |])+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Integer where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Integer) => CDF (Binomial b) Integer where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Integer) => PDF (Binomial b) Integer where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Int) => CDF (Binomial b) Int where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Int) => PDF (Binomial b) Int where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int8 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Int8) => CDF (Binomial b) Int8 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Int8) => PDF (Binomial b) Int8 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int16 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Int16) => CDF (Binomial b) Int16 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Int16) => PDF (Binomial b) Int16 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int32 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Int32) => CDF (Binomial b) Int32 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Int32) => PDF (Binomial b) Int32 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int64 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Int64) => CDF (Binomial b) Int64 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Int64) => PDF (Binomial b) Int64 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Word) => CDF (Binomial b) Word where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Word) => PDF (Binomial b) Word where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word8 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Word8) => CDF (Binomial b) Word8 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Word8) => PDF (Binomial b) Word8 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word16 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Word16) => CDF (Binomial b) Word16 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Word16) => PDF (Binomial b) Word16 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word32 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Word32) => CDF (Binomial b) Word32 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Word32) => PDF (Binomial b) Word32 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p+instance (Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word64 where+ rvarT (Binomial t p) = integralBinomial t p+instance (Real b, Distribution (Binomial b) Word64) => CDF (Binomial b) Word64 where+ cdf (Binomial t p) = integralBinomialCDF t p+instance (Real b, Distribution (Binomial b) Word64) => PDF (Binomial b) Word64 where+ pdf (Binomial t p) = integralBinomialPDF t p+ logPdf (Binomial t p) = integralBinomialLogPdf t p -$( replicateInstances ''Float realFloatTypes [d|- instance Distribution (Binomial b) Integer- => Distribution (Binomial b) Float- where rvar (Binomial t p) = floatingBinomial t p- instance CDF (Binomial b) Integer- => CDF (Binomial b) Float- where cdf (Binomial t p) = floatingBinomialCDF t p- instance PDF (Binomial b) Integer- => PDF (Binomial b) Float- where pdf (Binomial t p) = floatingBinomialPDF t p- logPdf (Binomial t p) = floatingBinomialLogPDF t p- |])+instance Distribution (Binomial b) Integer => Distribution (Binomial b) Float where+ rvar (Binomial t p) = floatingBinomial t p+instance CDF (Binomial b) Integer => CDF (Binomial b) Float where+ cdf (Binomial t p) = floatingBinomialCDF t p+instance PDF (Binomial b) Integer => PDF (Binomial b) Float where+ pdf (Binomial t p) = floatingBinomialPDF t p+ logPdf (Binomial t p) = floatingBinomialLogPDF t p+instance Distribution (Binomial b) Integer => Distribution (Binomial b) Double where+ rvar (Binomial t p) = floatingBinomial t p+instance CDF (Binomial b) Integer => CDF (Binomial b) Double where+ cdf (Binomial t p) = floatingBinomialCDF t p+instance PDF (Binomial b) Integer => PDF (Binomial b) Double where+ pdf (Binomial t p) = floatingBinomialPDF t p+ logPdf (Binomial t p) = floatingBinomialLogPDF t p
src/Data/Random/Distribution/Categorical.hs view
@@ -23,9 +23,7 @@ import Control.Arrow import Control.Monad import Control.Monad.ST-import Data.Foldable (Foldable(foldMap)) import Data.STRef-import Data.Traversable (Traversable(traverse, sequenceA)) import Data.List import Data.Function@@ -37,7 +35,7 @@ categorical :: (Num p, Distribution (Categorical p) a) => [(p,a)] -> RVar a categorical = rvar . fromList --- |Construct a 'Categorical' random process from a list of probabilities +-- |Construct a 'Categorical' random process from a list of probabilities -- and categories, where the probabilities all sum to 1. categoricalT :: (Num p, Distribution (Categorical p) a) => [(p,a)] -> RVarT m a categoricalT = rvarT . fromList@@ -47,7 +45,7 @@ weightedCategorical :: (Fractional p, Eq p, Distribution (Categorical p) a) => [(p,a)] -> RVar a weightedCategorical = rvar . fromWeightedList --- |Construct a 'Categorical' random process from a list of weights +-- |Construct a 'Categorical' random process from a list of weights -- and categories. The weights do /not/ have to sum to 1. weightedCategoricalT :: (Fractional p, Eq p, Distribution (Categorical p) a) => [(p,a)] -> RVarT m a weightedCategoricalT = rvarT . fromWeightedList@@ -73,14 +71,14 @@ numEvents :: Categorical p a -> Int numEvents (Categorical ds) = V.length ds --- |Construct a 'Categorical' distribution from a list of weighted categories, +-- |Construct a 'Categorical' distribution from a list of weighted categories, -- where the weights do not necessarily sum to 1. fromWeightedList :: (Fractional p, Eq p) => [(p,a)] -> Categorical p a fromWeightedList = normalizeCategoricalPs . fromList -- |Construct a 'Categorical' distribution from a list of observed outcomes. -- Equivalent events will be grouped and counted, and the probabilities of each--- event in the returned distribution will be proportional to the number of +-- event in the returned distribution will be proportional to the number of -- occurrences of that event. fromObservations :: (Fractional p, Eq p, Ord a) => [a] -> Categorical p a fromObservations = fromWeightedList . map (genericLength &&& head) . group . sort@@ -91,10 +89,10 @@ -- binary search. -- |Categorical distribution; a list of events with corresponding probabilities.--- The sum of the probabilities must be 1, and no event should have a zero +-- The sum of the probabilities must be 1, and no event should have a zero -- or negative probability (at least, at time of sampling; very clever users--- can do what they want with the numbers before sampling, just make sure --- that if you're one of those clever ones, you at least eliminate negative +-- can do what they want with the numbers before sampling, just make sure+-- that if you're one of those clever ones, you at least eliminate negative -- weights before sampling). newtype Categorical p a = Categorical (V.Vector (p, a)) deriving Eq@@ -117,19 +115,19 @@ | n == 1 = return (snd (V.head ds)) | otherwise = do u <- uniformT 0 (fst (V.last ds))- + let -- by construction, p is monotone; (i < j) ==> (p i <= p j) p i = fst (ds V.! i) x i = snd (ds V.! i)- + -- findEvent -- =========== -- invariants: (i <= j), (u <= p j), ((i == 0) || (p i < u)) -- (the last one means 'i' does not increase unless it bounds 'p' below 'u') -- variant: either i increases or j decreases. -- upon termination: ∀ k. if (k < j) then (p k < u) else (u <= p k)- -- (that is, the chosen event 'x j' is the first one whose - -- associated cumulative probability 'p j' is greater than + -- (that is, the chosen event 'x j' is the first one whose+ -- associated cumulative probability 'p j' is greater than -- or equal to 'u') findEvent i j | j <= i = x j@@ -139,7 +137,7 @@ -- midpoint rounding down -- (i < j) ==> (m < j) m = (i + j) `div` 2- + return $! if u <= 0 then x 0 else findEvent 0 (n-1) where n = V.length ds @@ -156,22 +154,22 @@ instance Fractional p => Monad (Categorical p) where return x = Categorical (V.singleton (1, x))- + -- I'm not entirely sure whether this is a valid form of failure; see next -- set of comments. #if __GLASGOW_HASKELL__ < 808 fail _ = Categorical V.empty #endif- + -- Should the normalize step be included here, or should normalization -- be assumed? It seems like there is (at least) 1 valid situation where- -- non-normal results would arise: the distribution being modeled is - -- "conditional" and some event arose that contradicted the assumed - -- condition and thus was eliminated ('f' returned an empty or + -- non-normal results would arise: the distribution being modeled is+ -- "conditional" and some event arose that contradicted the assumed+ -- condition and thus was eliminated ('f' returned an empty or -- zero-probability consequent, possibly by 'fail'ing).- -- + -- -- It seems reasonable to continue in such circumstances, but should there- -- be any renormalization? If so, does it make a difference when that + -- be any renormalization? If so, does it make a difference when that -- renormalization is done? I'm pretty sure it does, actually. So, the -- normalization will be omitted here for now, as it's easier for the -- user (who really better know what they mean if they're returning@@ -180,7 +178,7 @@ xs >>= f = {- normalizeCategoricalPs . -} fromList $ do (p, x) <- toList xs (q, y) <- toList (f x)- + return (p * q, y) instance Fractional p => Applicative (Categorical p) where@@ -191,7 +189,7 @@ mapCategoricalPs :: (Num p, Num q) => (p -> q) -> Categorical p e -> Categorical q e mapCategoricalPs f = fromList . map (first f) . toList --- |Adjust all the weights of a categorical distribution so that they +-- |Adjust all the weights of a categorical distribution so that they -- sum to unity and remove all events whose probability is zero. normalizeCategoricalPs :: (Fractional p, Eq p) => Categorical p e -> Categorical p e normalizeCategoricalPs orig@(Categorical ds)@@ -200,13 +198,13 @@ lastP <- newSTRef 0 nDups <- newSTRef 0 normalized <- V.thaw ds- + let n = V.length ds skip = modifySTRef' nDups (1+) save i p x = do d <- readSTRef nDups MV.write normalized (i-d) (p, x)- + sequence_ [ do let (p,x) = ds V.! i@@ -218,7 +216,7 @@ writeSTRef lastP $! p | i <- [0..n-1] ]- + -- force last element to 1 d <- readSTRef nDups let n' = n-d@@ -242,14 +240,14 @@ -- event will have a probability equal to the sum of all the originals). collectEvents :: (Ord e, Num p, Ord p) => Categorical p e -> Categorical p e collectEvents = collectEventsBy compare ((sum *** head) . unzip)- + -- |Simplify a categorical distribution by combining equivalent events (the new -- event will have a weight equal to the sum of all the originals). -- The comparator function is used to identify events to combine. Once chosen, -- the events and their weights are combined by the provided probability and -- event aggregation function. collectEventsBy :: Num p => (e -> e -> Ordering) -> ([(p,e)] -> (p,e))-> Categorical p e -> Categorical p e-collectEventsBy compareE combine = +collectEventsBy compareE combine = fromList . map combine . groupEvents . sortEvents . toList where groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)
src/Data/Random/Distribution/Dirichlet.hs view
@@ -18,7 +18,7 @@ fractionalDirichlet as = do xs <- sequence [gammaT a 1 | a <- as] let total = foldl1' (+) xs- + return (map (* recip total) xs) dirichlet :: Distribution Dirichlet [a] => [a] -> RVar [a]
src/Data/Random/Distribution/Gamma.hs view
@@ -9,10 +9,10 @@ module Data.Random.Distribution.Gamma ( Gamma(..) , gamma, gammaT- + , Erlang(..) , erlang, erlangT- + , mtGamma ) where @@ -31,10 +31,10 @@ {-# SPECIALIZE mtGamma :: Float -> Float -> RVarT m Float #-} mtGamma :: (Floating a, Ord a,- Distribution StdUniform a, + Distribution StdUniform a, Distribution Normal a) => a -> a -> RVarT m a-mtGamma a b +mtGamma a b | a < 1 = do u <- stdUniformT mtGamma (1+a) $! (b * u ** recip a)@@ -42,11 +42,11 @@ where !d = a - fromRational (1%3) !c = recip (sqrt (9*d))- + go = do x <- stdNormalT let !v = 1 + c*x- + if v <= 0 then go else do@@ -89,4 +89,3 @@ instance (Integral a, Real b, Distribution (Erlang a) b) => CDF (Erlang a) b where cdf (Erlang a) x = incompleteGamma (fromIntegral a) (realToFrac x)-
src/Data/Random/Distribution/Multinomial.hs view
@@ -24,9 +24,9 @@ go n (p:ps) (psum:psums) f = do x <- binomialT n (p / psum) go (n-x) ps psums (f . (x:))- + go _ _ _ _ = error "rvar/Multinomial: programming error! this case should be impossible!"- + -- less wasteful version of (map sum . tails) tailSums [] = [0] tailSums (x:xs) = case tailSums xs of
src/Data/Random/Distribution/Normal.hs view
@@ -1,6 +1,6 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,- UndecidableInstances, ForeignFunctionInterface, BangPatterns, + UndecidableInstances, ForeignFunctionInterface, BangPatterns, RankNTypes #-} @@ -10,26 +10,25 @@ ( Normal(..) , normal, normalT , stdNormal, stdNormalT- + , doubleStdNormal , floatStdNormal , realFloatStdNormal- + , normalTail- + , normalPair , boxMullerNormalPair , knuthPolarNormalPair ) where -import Data.Random.Internal.Words import Data.Bits -import Data.Random.Source import Data.Random.Distribution import Data.Random.Distribution.Uniform import Data.Random.Distribution.Ziggurat import Data.Random.RVar+import Data.Word import Data.Vector.Generic (Vector) import qualified Data.Vector as V@@ -37,6 +36,8 @@ import Data.Number.Erf +import qualified System.Random.Stateful as Random+ -- |A random variable that produces a pair of independent -- normally-distributed values. normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)@@ -44,7 +45,7 @@ -- |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 +-- 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).@@ -55,27 +56,27 @@ t <- stdUniform let r = sqrt (-2 * log u) theta = (2 * pi) * t- + x = r * cos theta 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 +-- 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 v1 <- uniform (-1) 1 v2 <- uniform (-1) 1- + let s = v1*v1 + v2*v2 if s >= 1 then knuthPolarNormalPair else return $ if s == 0 then (0,0)- else let scale = sqrt (-2 * log s / s) + 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)@@ -110,7 +111,7 @@ normalFInv y = sqrt ((-2) * log y) -- | integral of 'normalF' normalFInt :: (Floating a, Erf a, Ord a) => a -> a-normalFInt x +normalFInt x | x <= 0 = 0 | otherwise = normalFVol * erf (x * sqrt 0.5) -- | volume of 'normalF'@@ -120,8 +121,8 @@ -- |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'--- class). --- +-- class).+-- -- 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@@ -135,13 +136,13 @@ -- @Distribution Normal@ instance declaration. realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVarT m a realFloatStdNormal = runZiggurat (normalZ p getIU `asTypeOf` (undefined :: Ziggurat V.Vector a))- where + where p :: Int p = 6- + getIU :: (Num a, Distribution Uniform a) => RVarT m (Int, a) getIU = do- i <- getRandomWord8+ i <- Random.uniformWord8 RGen u <- uniformT (-1) 1 return (fromIntegral i .&. (2^p-1), u) @@ -159,18 +160,32 @@ {-# NOINLINE doubleStdNormalZ #-} doubleStdNormalZ :: Ziggurat UV.Vector Double-doubleStdNormalZ = mkZiggurat_ True - normalF normalFInv - doubleStdNormalC doubleStdNormalR doubleStdNormalV +doubleStdNormalZ = mkZiggurat_ True+ normalF normalFInv+ doubleStdNormalC doubleStdNormalR doubleStdNormalV getIU (normalTail doubleStdNormalR)- where + where getIU :: RVarT m (Int, Double) getIU = do- !w <- getRandomWord64+ !w <- Random.uniformWord64 RGen let (u,i) = wordToDoubleWithExcess w return $! (fromIntegral i .&. (doubleStdNormalC-1), u+u-1) +-- NOTE: inlined from random-source+{-# INLINE wordToDouble #-}+-- |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)++{-# INLINE wordToDoubleWithExcess #-}+-- |Same as wordToDouble, but also return the unused bits (as the 12+-- least significant bits of a 'Word64')+wordToDoubleWithExcess :: Word64 -> (Double, Word64)+wordToDoubleWithExcess x = (wordToDouble x, x `shiftR` 52)++ -- |A random variable sampling from the standard normal distribution -- over the 'Float' type. floatStdNormal :: RVarT m Float@@ -185,17 +200,31 @@ {-# NOINLINE floatStdNormalZ #-} floatStdNormalZ :: Ziggurat UV.Vector Float-floatStdNormalZ = mkZiggurat_ True - normalF normalFInv - floatStdNormalC floatStdNormalR floatStdNormalV +floatStdNormalZ = mkZiggurat_ True+ normalF normalFInv+ floatStdNormalC floatStdNormalR floatStdNormalV getIU (normalTail floatStdNormalR) where getIU :: RVarT m (Int, Float) getIU = do- !w <- getRandomWord32+ !w <- Random.uniformWord32 RGen let (u,i) = word32ToFloatWithExcess w return (fromIntegral i .&. (floatStdNormalC-1), u+u-1)++-- NOTE: inlined from random-source+{-# INLINE word32ToFloat #-}+-- |Pack the low 23 bits from a 'Word32' into a 'Float' in the range [0,1).+-- Used to convert a 'stdUniform' 'Word32' to a 'stdUniform' 'Double'.+word32ToFloat :: Word32 -> Float+word32ToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)++{-# INLINE word32ToFloatWithExcess #-}+-- |Same as word32ToFloat, but also return the unused bits (as the 9+-- least significant bits of a 'Word32')+word32ToFloatWithExcess :: Word32 -> (Float, Word32)+word32ToFloatWithExcess x = (word32ToFloat x, x `shiftR` 23)+ normalCdf :: (Real a) => a -> a -> a -> Double normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s)
src/Data/Random/Distribution/Poisson.hs view
@@ -1,15 +1,12 @@ {-# LANGUAGE MultiParamTypeClasses,- FlexibleInstances, FlexibleContexts, UndecidableInstances,- TemplateHaskell+ FlexibleInstances, FlexibleContexts, UndecidableInstances #-} {-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-} module Data.Random.Distribution.Poisson where -import Data.Random.Internal.TH- import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Uniform@@ -18,6 +15,9 @@ import Control.Monad +import Data.Int+import Data.Word+ -- from Knuth, with interpretation help from gsl sources integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVarT m a integralPoisson = psn 0@@ -87,24 +87,82 @@ newtype Poisson b a = Poisson b -$( replicateInstances ''Int integralTypes [d|- instance ( RealFloat b- , Distribution StdUniform b- , Distribution (Erlang Int) b- , Distribution (Binomial b) Int- ) => Distribution (Poisson b) Int where- rvarT (Poisson mu) = integralPoisson mu- instance (Real b, Distribution (Poisson b) Int) => CDF (Poisson b) Int where- cdf (Poisson mu) = integralPoissonCDF mu- instance (Real b, Distribution (Poisson b) Int) => PDF (Poisson b) Int where- pdf (Poisson mu) = integralPoissonPDF mu- |] )+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Integer) b, Distribution (Binomial b) Integer) => Distribution (Poisson b) Integer where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Integer) => CDF (Poisson b) Integer where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Integer) => PDF (Poisson b) Integer where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int) b, Distribution (Binomial b) Int) => Distribution (Poisson b) Int where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Int) => CDF (Poisson b) Int where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Int) => PDF (Poisson b) Int where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int8) b, Distribution (Binomial b) Int8) => Distribution (Poisson b) Int8 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Int8) => CDF (Poisson b) Int8 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Int8) => PDF (Poisson b) Int8 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int16) b, Distribution (Binomial b) Int16) => Distribution (Poisson b) Int16 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Int16) => CDF (Poisson b) Int16 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Int16) => PDF (Poisson b) Int16 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int32) b, Distribution (Binomial b) Int32) => Distribution (Poisson b) Int32 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Int32) => CDF (Poisson b) Int32 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Int32) => PDF (Poisson b) Int32 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int64) b, Distribution (Binomial b) Int64) => Distribution (Poisson b) Int64 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Int64) => CDF (Poisson b) Int64 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Int64) => PDF (Poisson b) Int64 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word) b, Distribution (Binomial b) Word) => Distribution (Poisson b) Word where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Word) => CDF (Poisson b) Word where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Word) => PDF (Poisson b) Word where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word8) b, Distribution (Binomial b) Word8) => Distribution (Poisson b) Word8 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Word8) => CDF (Poisson b) Word8 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Word8) => PDF (Poisson b) Word8 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word16) b, Distribution (Binomial b) Word16) => Distribution (Poisson b) Word16 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Word16) => CDF (Poisson b) Word16 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Word16) => PDF (Poisson b) Word16 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word32) b, Distribution (Binomial b) Word32) => Distribution (Poisson b) Word32 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Word32) => CDF (Poisson b) Word32 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Word32) => PDF (Poisson b) Word32 where+ pdf (Poisson mu) = integralPoissonPDF mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word64) b, Distribution (Binomial b) Word64) => Distribution (Poisson b) Word64 where+ rvarT (Poisson mu) = integralPoisson mu+instance (Real b, Distribution (Poisson b) Word64) => CDF (Poisson b) Word64 where+ cdf (Poisson mu) = integralPoissonCDF mu+instance (Real b, Distribution (Poisson b) Word64) => PDF (Poisson b) Word64 where+ pdf (Poisson mu) = integralPoissonPDF mu -$( replicateInstances ''Float realFloatTypes [d|- instance (Distribution (Poisson b) Integer) => Distribution (Poisson b) Float where- rvarT (Poisson mu) = fractionalPoisson mu- instance (CDF (Poisson b) Integer) => CDF (Poisson b) Float where- cdf (Poisson mu) = fractionalPoissonCDF mu- instance (PDF (Poisson b) Integer) => PDF (Poisson b) Float where- pdf (Poisson mu) = fractionalPoissonPDF mu- |])+instance Distribution (Poisson b) Integer => Distribution (Poisson b) Float where+ rvarT (Poisson mu) = fractionalPoisson mu+instance CDF (Poisson b) Integer => CDF (Poisson b) Float where+ cdf (Poisson mu) = fractionalPoissonCDF mu+instance PDF (Poisson b) Integer => PDF (Poisson b) Float where+ pdf (Poisson mu) = fractionalPoissonPDF mu+instance Distribution (Poisson b) Integer => Distribution (Poisson b) Double where+ rvarT (Poisson mu) = fractionalPoisson mu+instance CDF (Poisson b) Integer => CDF (Poisson b) Double where+ cdf (Poisson mu) = fractionalPoissonCDF mu+instance PDF (Poisson b) Integer => PDF (Poisson b) Double where+ pdf (Poisson mu) = fractionalPoissonPDF mu
src/Data/Random/Distribution/Rayleigh.hs view
@@ -1,5 +1,5 @@ {-# LANGUAGE- MultiParamTypeClasses, + MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, UndecidableInstances #-}@@ -19,7 +19,7 @@ -- |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
src/Data/Random/Distribution/Triangular.hs view
@@ -52,7 +52,7 @@ = realToFrac (1 - (c - x)^(2 :: Int) / ((c - a) * (c - b))) | otherwise = 1- + instance (RealFloat a, Ord a, Distribution StdUniform a) => Distribution Triangular a where rvarT (Triangular a b c) = floatingTriangular a b c instance (RealFrac a, Distribution Triangular a) => CDF Triangular a where
src/Data/Random/Distribution/Uniform.hs view
@@ -37,11 +37,9 @@ , enumUniformCDF ) where -import Data.Random.Internal.TH-import Data.Random.Internal.Words+ import Data.Random.Internal.Fixed -import Data.Random.Source import Data.Random.Distribution import Data.Random.RVar @@ -51,39 +49,14 @@ import Control.Monad.Loops +import qualified System.Random.Stateful as Random+ -- |Compute a random 'Integral' value between the 2 values provided (inclusive). {-# INLINE integralUniform #-}-integralUniform :: (Integral a) => a -> a -> RVarT m a-integralUniform !x !y = if x < y then integralUniform' x y else integralUniform' y x--{-# SPECIALIZE integralUniform' :: Int -> Int -> RVarT m Int #-}-{-# SPECIALIZE integralUniform' :: Int8 -> Int8 -> RVarT m Int8 #-}-{-# SPECIALIZE integralUniform' :: Int16 -> Int16 -> RVarT m Int16 #-}-{-# SPECIALIZE integralUniform' :: Int32 -> Int32 -> RVarT m Int32 #-}-{-# SPECIALIZE integralUniform' :: Int64 -> Int64 -> RVarT m Int64 #-}-{-# SPECIALIZE integralUniform' :: Word -> Word -> RVarT m Word #-}-{-# SPECIALIZE integralUniform' :: Word8 -> Word8 -> RVarT m Word8 #-}-{-# SPECIALIZE integralUniform' :: Word16 -> Word16 -> RVarT m Word16 #-}-{-# SPECIALIZE integralUniform' :: Word32 -> Word32 -> RVarT m Word32 #-}-{-# SPECIALIZE integralUniform' :: Word64 -> Word64 -> RVarT m Word64 #-}-{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVarT m Integer #-}-integralUniform' :: (Integral a) => a -> a -> RVarT m a-integralUniform' !l !u- | nReject == 0 = fmap shift prim- | otherwise = fmap shift loop- where- m = 1 + toInteger u - toInteger l- (bytes, nPossible) = bytesNeeded m- nReject = nPossible `mod` m-- !prim = getRandomNByteInteger bytes- !shift = \(!z) -> l + (fromInteger $! (z `mod` m))-- loop = do- z <- prim- if z < nReject- then loop- else return z+integralUniform :: Random.UniformRange a => a -> a -> RVarT m a+integralUniform !x !y = Random.uniformRM (x, y) RGen+ -- Maybe switch to uniformIntegralM (requires exposing from `random` internals):+ -- Random.uniformIntegralM (x, y) RGen integralUniformCDF :: (Integral a, Fractional b) => a -> a -> a -> b integralUniformCDF a b x@@ -92,15 +65,6 @@ | x > b = 1 | otherwise = (fromIntegral x - fromIntegral a) / (fromIntegral b - fromIntegral a) --- TODO: come up with a decent, fast heuristic to decide whether to return an extra--- byte. May involve moving calculation of nReject into this function, and then--- accepting first if 4*nReject < nPossible or something similar.-bytesNeeded :: Integer -> (Int, Integer)-bytesNeeded x = head (dropWhile ((<= x).snd) powersOf256)--powersOf256 :: [(Int, Integer)]-powersOf256 = zip [0..] (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@@ -120,13 +84,17 @@ -- |Compute a uniform random 'Float' value in the range [0,1) floatStdUniform :: RVarT m Float floatStdUniform = do- x <- getRandomWord32- return (word32ToFloat x)+ x <- uniformRangeRVarT (0, 1)+ -- exclude 1. TODO: come up with something smarter+ if x == 1 then floatStdUniform else pure x -- |Compute a uniform random 'Double' value in the range [0,1) {-# INLINE doubleStdUniform #-} doubleStdUniform :: RVarT m Double-doubleStdUniform = getRandomDouble+doubleStdUniform = do+ x <- uniformRangeRVarT (0, 1)+ -- exclude 1. TODO: come up with something smarter+ if x == 1 then doubleStdUniform else pure x -- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type realFloatStdUniform :: RealFloat a => RVarT m a@@ -284,32 +252,42 @@ -- (that is, 0 to 1 including 0 but not including 1). data StdUniform t = StdUniform -$( replicateInstances ''Int integralTypes [d|- instance Distribution Uniform Int where rvarT (Uniform a b) = integralUniform a b- instance CDF Uniform Int where cdf (Uniform a b) = integralUniformCDF a b- |])+instance Distribution Uniform Integer where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Integer where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Int where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Int where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Int8 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Int8 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Int16 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Int16 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Int32 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Int32 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Int64 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Int64 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Word where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Word where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Word8 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Word8 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Word16 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Word16 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Word32 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Word32 where cdf (Uniform a b) = integralUniformCDF a b+instance Distribution Uniform Word64 where rvarT (Uniform a b) = integralUniform a b+instance CDF Uniform Word64 where cdf (Uniform a b) = integralUniformCDF a b -instance Distribution StdUniform Word8 where rvarT _ = getRandomWord8-instance Distribution StdUniform Word16 where rvarT _ = getRandomWord16-instance Distribution StdUniform Word32 where rvarT _ = getRandomWord32-instance Distribution StdUniform Word64 where rvarT _ = getRandomWord64+instance Distribution StdUniform Word8 where rvarT _ = Random.uniformWord8 RGen+instance Distribution StdUniform Word16 where rvarT _ = Random.uniformWord16 RGen+instance Distribution StdUniform Word32 where rvarT _ = Random.uniformWord32 RGen+instance Distribution StdUniform Word64 where rvarT _ = Random.uniformWord64 RGen+instance Distribution StdUniform Word where rvarT _ = uniformRVarT -instance Distribution StdUniform Int8 where rvarT _ = fromIntegral `fmap` getRandomWord8-instance Distribution StdUniform Int16 where rvarT _ = fromIntegral `fmap` getRandomWord16-instance Distribution StdUniform Int32 where rvarT _ = fromIntegral `fmap` getRandomWord32-instance Distribution StdUniform Int64 where rvarT _ = fromIntegral `fmap` getRandomWord64+instance Distribution StdUniform Int8 where rvarT _ = uniformRVarT+instance Distribution StdUniform Int16 where rvarT _ = uniformRVarT+instance Distribution StdUniform Int32 where rvarT _ = uniformRVarT+instance Distribution StdUniform Int64 where rvarT _ = uniformRVarT -instance Distribution StdUniform Int where- rvar _ =- $(if toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)- then [|fromIntegral `fmap` getRandomWord64 :: RVar Int|]- else [|fromIntegral `fmap` getRandomWord32 :: RVar Int|])+instance Distribution StdUniform Int where rvarT _ = uniformRVarT -instance Distribution StdUniform Word where- rvar _ =- $(if toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)- then [|fromIntegral `fmap` getRandomWord64 :: RVar Word|]- else [|fromIntegral `fmap` getRandomWord32 :: RVar Word|]) -- Integer has no StdUniform... @@ -331,7 +309,7 @@ instance CDF Uniform Double where cdf (Uniform a b) = realUniformCDF a b instance Distribution StdUniform Float where rvarT _ = floatStdUniform-instance Distribution StdUniform Double where rvarT _ = getRandomDouble+instance Distribution StdUniform Double where rvarT _ = uniformRangeRVarT (0, 1) instance CDF StdUniform Float where cdf _ = realStdUniformCDF instance CDF StdUniform Double where cdf _ = realStdUniformCDF instance PDF StdUniform Float where pdf _ = realStdUniformPDF@@ -349,15 +327,17 @@ instance Distribution Uniform () where rvarT (Uniform _ _) = return () instance CDF Uniform () where cdf (Uniform _ _) = return 1-$( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|- instance Distribution Uniform Char where rvarT (Uniform a b) = enumUniform a b- instance CDF Uniform Char where cdf (Uniform a b) = enumUniformCDF a b - |])+instance Distribution Uniform Char where rvarT (Uniform a b) = enumUniform a b+instance CDF Uniform Char where cdf (Uniform a b) = enumUniformCDF a b+instance Distribution Uniform Bool where rvarT (Uniform a b) = enumUniform a b+instance CDF Uniform Bool where cdf (Uniform a b) = enumUniformCDF a b+instance Distribution Uniform Ordering where rvarT (Uniform a b) = enumUniform a b+instance CDF Uniform Ordering where cdf (Uniform a b) = enumUniformCDF a b instance Distribution StdUniform () where rvarT ~StdUniform = return () instance CDF StdUniform () where cdf ~StdUniform = return 1-instance Distribution StdUniform Bool where rvarT ~StdUniform = fmap even (getRandomWord8)+instance Distribution StdUniform Bool where rvarT ~StdUniform = uniformRVarT instance CDF StdUniform Bool where cdf ~StdUniform = boundedEnumStdUniformCDF instance Distribution StdUniform Char where rvarT ~StdUniform = boundedEnumStdUniform
src/Data/Random/Distribution/Weibull.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, FlexibleContexts #-} module Data.Random.Distribution.Weibull where import Data.Random.Distribution
src/Data/Random/Distribution/Ziggurat.hs view
@@ -7,16 +7,16 @@ -- |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 +-- '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(..)@@ -48,10 +48,10 @@ data Ziggurat v 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 + -- 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.@@ -64,8 +64,8 @@ -- -- * 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 + --+ -- * 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@@ -74,21 +74,21 @@ -- a double (using 52 bits) and a bin number (using up to 12 bits), -- for example. zGetIU :: !(forall m. RVarT m (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 :: (forall m. RVarT m 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 :: !(forall m. t -> t -> RVarT m 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 in -- its native form only samples from one-sided distributions.@@ -113,7 +113,7 @@ -- (or 0 to 1 if not mirroring the distribution). -- Let X be U scaled to the size of the selected bin. (!i,!u) <- zGetIU- + -- 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@@ -123,7 +123,7 @@ else if i == 0 then sampleTail u else sampleGreyArea i $! (u * zTable_xs ! i)- + -- 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.@@ -133,7 +133,7 @@ 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)@@ -143,28 +143,28 @@ | 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) 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 the 'zGetIU' random tuple--- +-- -- * an RVar sampling from the tail (the region where x > R)--- +-- {-# INLINE mkZiggurat_ #-} {-# SPECIALIZE mkZiggurat_ :: Bool -> (Float -> Float) -> (Float -> Float) -> Int -> Float -> Float -> (forall m. RVarT m (Int, Float)) -> (forall m. RVarT m Float ) -> Ziggurat UV.Vector Float #-} {-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. RVarT m (Int, Double)) -> (forall m. RVarT m Double) -> Ziggurat UV.Vector Double #-}@@ -191,13 +191,13 @@ , zTailDist = tailDist , zMirror = m }- where + where xs = zigguratTable f fInv c r v --- |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 -- argument for the tail distribution as a function of the selected -- R value.@@ -213,15 +213,15 @@ -> (forall m. t -> RVarT m t) -> Ziggurat v t mkZiggurat m f fInv fInt fVol c getIU tailDist =- mkZiggurat_ m f fInv c r v getIU (tailDist r) + mkZiggurat_ m f fInv c r v getIU (tailDist r) where (r,v) = findBin0 c f fInv fInt fVol -- |Build a lazy recursive ziggurat. Uses a lazily-constructed ziggurat -- as its tail distribution (with another as its tail, ad nauseam).--- +-- -- Arguments:--- +-- -- * flag indicating whether to mirror the distribution -- -- * the (one-sided antitone) PDF, not necessarily normalized@@ -254,7 +254,7 @@ fix g = g (fix g) z = mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU z)) -mkTail :: +mkTail :: (RealFloat a, Vector v a, Distribution Uniform a) => Bool -> (a -> a) -> (a -> a) -> (a -> a)@@ -269,16 +269,16 @@ return (x + r * signum x) where fIntR = fInt r- + f' x | x < 0 = f r | otherwise = f (x+r) fInv' = subtract r . fInv fInt' x | x < 0 = 0 | otherwise = fInt (x+r) - fIntR- + fVol' = fVol - fIntR- + zigguratTable :: (Fractional a, Vector v a, Ord a) => (a -> a) -> (a -> a) -> Int -> a -> a -> v a zigguratTable f fInv c r v = case zigguratXs f fInv c r v of@@ -292,19 +292,19 @@ where xs = Prelude.map x [0..c] -- sample c x ys = Prelude.map f xs- + x 0 = v / f r x 1 = r x i | i == c = 0 x i | i > 1 = next (i-1) x _ = error "zigguratXs: programming error! this case should be impossible!"- + next i = let x_i = xs!!i in if x_i <= 0 then -1 else fInv (ys!!i + (v / x_i))- - excess = xs!!(c-1) * (f 0 - ys !! (c-1)) - v + excess = xs!!(c-1) * (f 0 - ys !! (c-1)) - v + precomputeRatios :: (Vector v a, Fractional a) => v a -> v a precomputeRatios zTable_xs = generate (c-1) $ \i -> zTable_xs!(i+1) / zTable_xs!i where@@ -314,7 +314,7 @@ -- Search the distribution for an appropriate R and V. -- -- Arguments:--- +-- -- * Number of bins -- -- * target function (one-sided antitone PDF, not necessarily normalized)@@ -326,20 +326,20 @@ -- * estimate of total volume under function (integral from 0 to infinity) -- -- Result: (R,V)-findBin0 :: (RealFloat b) => +findBin0 :: (RealFloat b) => Int -> (b -> b) -> (b -> b) -> (b -> b) -> b -> (b, b) findBin0 cInt f fInv fInt fVol = (rMin,v rMin) where c = fromIntegral cInt v r = r * f r + fVol - fInt r- + -- initial R guess: r0 = findMin (\r -> v r <= fVol / c) -- find a better R:- rMin = findMinFrom r0 1 $ \r -> - let e = exc r + rMin = findMinFrom r0 1 $ \r ->+ let e = exc r in e >= 0 && not (isNaN e)- + exc x = zigguratExcess f fInv cInt x (v x) instance (Num t, Ord t, Vector v t) => Distribution (Ziggurat v) t where
src/Data/Random/Internal/Find.hs view
@@ -19,7 +19,7 @@ -- specified point with the specified stepsize, performs an exponential -- search out from there until it finds an interval bracketing the -- change-point of the predicate, and then performs a bisection search--- to isolate the change point. Note that infinitely-divisible domains +-- to isolate the change point. Note that infinitely-divisible domains -- such as 'Rational' cannot be searched by this function because it does -- not terminate until it reaches a point where further subdivision of the -- interval has no effect.@@ -33,31 +33,31 @@ -- a feasible answer fixZero 0 = 0 fixZero z = z- + -- preconditions: -- not (p l) -- 0 <= l < x- ascend l x + ascend l x | p x = bisect l x | otherwise = ascend x $! 2*x-z0- + -- preconditions: -- p h -- x < h <= 0- descend x h + descend x h | p x = (descend $! 2*x-z0) x | otherwise = bisect x h- + -- preconditions: -- not (p l) -- p h -- l <= h- bisect l h + bisect l h | l /< h = h | l /< mid || mid /< h = if p mid then mid else h | p mid = bisect l mid | otherwise = bisect mid h- where + where a /< b = not (a < b) mid = (l+h)*0.5
src/Data/Random/Internal/Fixed.hs view
@@ -35,7 +35,7 @@ -- |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 +-- 'unMkFixed' as defined here will work. Both are implemented using -- 'unsafeCoerce'. mkFixed :: Integer -> Fixed r mkFixed = unsafeCoerce
− src/Data/Random/Internal/TH.hs
@@ -1,79 +0,0 @@-{-# LANGUAGE- TemplateHaskell- #-}---- |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--import Data.Word-import Data.Int-import Control.Monad---- |Names of standard 'Integral' types-integralTypes :: [Name]-integralTypes = - [ ''Integer- , ''Int, ''Int8, ''Int16, ''Int32, ''Int64- , ''Word, ''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 getDecls = liftM concat $ sequence- [ do- decls <- getDecls- sequence- [ everywhereM (mkM (return . replaceName standin t)) dec- | dec <- decls- ]- | t <- types- ]-
src/Data/Random/Lift.hs view
@@ -5,7 +5,6 @@ import Data.RVar import qualified Data.Functor.Identity as T import qualified Control.Monad.Trans.Class as T-import Data.Random.Source.Std #ifndef MTL2 import qualified Control.Monad.Identity as MTL@@ -19,10 +18,10 @@ -- For instances where 'm' and 'n' have 'return'/'pure' defined, -- these instances must satisfy -- @lift (return x) == return x@.--- +-- -- This form of 'lift' has an extremely general type and is used primarily to -- support 'sample'. Its excessive generality is the main reason it's not--- exported from "Data.Random". 'RVarT' is, however, an instance of +-- exported from "Data.Random". 'RVarT' is, however, an instance of -- 'T.MonadTrans', which in most cases is the preferred way -- to do the lifting. class Lift m n where@@ -41,7 +40,7 @@ lift = return . T.runIdentity instance Lift (RVarT T.Identity) (RVarT m) where- lift x = runRVar x StdRandom+ lift x = runRVar x RGen -- | This instance is again incoherent with the others, but provides a -- more-specific instance to resolve the overlap between the@@ -58,7 +57,7 @@ lift = return . MTL.runIdentity instance Lift (RVarT MTL.Identity) (RVarT m) where- lift x = runRVarTWith (return . MTL.runIdentity) x StdRandom+ lift x = runRVarTWith (return . MTL.runIdentity) x RGen -- | This instance is again incoherent with the others, but provides a -- more-specific instance to resolve the overlap between the@@ -67,4 +66,3 @@ lift = T.lift #endif-
src/Data/Random/List.hs view
@@ -28,11 +28,11 @@ shuffleT [] = return [] shuffleT xs = do is <- zipWithM (\_ i -> uniformT 0 i) (tail xs) [1..]- + return (SRS.shuffle xs (reverse is)) -- | A random variable that shuffles a list of a known length (or a list--- prefix of the specified length). Useful for shuffling large lists when +-- prefix of the specified 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. shuffleN :: Int -> [a] -> RVar [a]@@ -53,4 +53,3 @@ is <- sequence [uniformT 0 i | i <- take n [m-1, m-2 ..1]] return (take n $ SRS.shuffle (take m xs) is) shuffleNofMT _ _ _ = error "shuffleNofMT: negative length specified"-
src/Data/Random/RVar.hs view
@@ -2,13 +2,14 @@ module Data.Random.RVar ( RVar, runRVar , RVarT, runRVarT, runRVarTWith+ , RGen(..), uniformRVarT, uniformRangeRVarT ) where import Data.Random.Lift-import Data.Random.Internal.Source import Data.RVar hiding (runRVarT)+import System.Random.Stateful --- |Like 'runRVarTWith', but using an implicit lifting (provided by the +-- |Like 'runRVarTWith', but using an implicit lifting (provided by the -- 'Lift' class)-runRVarT :: (Lift n m, RandomSource m s) => RVarT n a -> s -> m a+runRVarT :: (Lift n m, StatefulGen g m) => RVarT n a -> g -> m a runRVarT = runRVarTWith lift
src/Data/Random/Sample.hs view
@@ -1,6 +1,6 @@ {-# LANGUAGE MultiParamTypeClasses,- FlexibleInstances, FlexibleContexts, + FlexibleInstances, FlexibleContexts, IncoherentInstances #-} @@ -8,39 +8,43 @@ module Data.Random.Sample where -import Control.Monad.State +import Control.Monad.State+import Control.Monad.Reader import Data.Random.Distribution import Data.Random.Lift import Data.Random.RVar-import Data.Random.Source-import Data.Random.Source.Std +import System.Random.Stateful+ -- |A typeclass allowing 'Distribution's and 'RVar's to be sampled. Both may -- also be sampled via 'runRVar' or 'runRVarT', but I find it psychologically -- pleasing to be able to sample both using this function, as they are two -- separate abstractions for one base concept: a random variable. class Sampleable d m t where -- |Directly sample from a distribution or random variable, using the given source of entropy.- sampleFrom :: RandomSource m s => s -> d t -> m t+ sampleFrom :: StatefulGen g m => g -> d t -> m t instance Distribution d t => Sampleable d m t where- sampleFrom src d = runRVarT (rvar d) src+ sampleFrom gen d = runRVarT (rvar d) gen -- 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+ sampleFrom gen x = runRVarT x gen -- |Sample a random variable using the default source of entropy for the -- monad in which the sampling occurs.-sample :: (Sampleable d m t, MonadRandom m) => d t -> m t-sample = sampleFrom StdRandom+sample :: (Sampleable d m t, StatefulGen g m, MonadReader g m) => d t -> m t+sample thing = ask >>= \gen -> sampleFrom gen thing -- |Sample a random variable in a \"functional\" style. Typical instantiations -- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.-sampleState :: (Sampleable d (State s) t, MonadRandom (State s)) => d t -> s -> (t, s)-sampleState thing = runState (sample thing)+-- sample :: (Distribution d a, StatefulGen g m, MonadReader g m) => d t -> m t+-- sample thing gen = runStateGen gen (\stateGen -> sampleFrom stateGen thing) --- |Sample a random variable in a \"semi-functional\" style. Typical instantiations--- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.-sampleStateT :: (Sampleable d (StateT s m) t, MonadRandom (StateT s m)) => d t -> s -> m (t, s)-sampleStateT thing = runStateT (sample thing)+sampleState :: (Distribution d t, RandomGen g, MonadState g m) => d t -> m t+sampleState thing = sampleFrom StateGenM thing++-- |Sample a random variable in a \"functional\" style. Typical instantiations+-- of @g@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.+samplePure :: (Distribution d t, RandomGen g) => d t -> g -> (t, g)+samplePure thing gen = runStateGen gen (\stateGen -> sampleFrom stateGen thing)