diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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. diff --git a/random-fu.cabal b/random-fu.cabal --- a/random-fu.cabal +++ b/random-fu.cabal @@ -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 diff --git a/src/Data/Random.hs b/src/Data/Random.hs --- a/src/Data/Random.hs +++ b/src/Data/Random.hs @@ -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) diff --git a/src/Data/Random/Distribution.hs b/src/Data/Random/Distribution.hs --- a/src/Data/Random/Distribution.hs +++ b/src/Data/Random/Distribution.hs @@ -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 diff --git a/src/Data/Random/Distribution/Bernoulli.hs b/src/Data/Random/Distribution/Bernoulli.hs --- a/src/Data/Random/Distribution/Bernoulli.hs +++ b/src/Data/Random/Distribution/Bernoulli.hs @@ -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) diff --git a/src/Data/Random/Distribution/Beta.hs b/src/Data/Random/Distribution/Beta.hs --- a/src/Data/Random/Distribution/Beta.hs +++ b/src/Data/Random/Distribution/Beta.hs @@ -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 diff --git a/src/Data/Random/Distribution/Binomial.hs b/src/Data/Random/Distribution/Binomial.hs --- a/src/Data/Random/Distribution/Binomial.hs +++ b/src/Data/Random/Distribution/Binomial.hs @@ -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 diff --git a/src/Data/Random/Distribution/Categorical.hs b/src/Data/Random/Distribution/Categorical.hs --- a/src/Data/Random/Distribution/Categorical.hs +++ b/src/Data/Random/Distribution/Categorical.hs @@ -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) diff --git a/src/Data/Random/Distribution/Dirichlet.hs b/src/Data/Random/Distribution/Dirichlet.hs --- a/src/Data/Random/Distribution/Dirichlet.hs +++ b/src/Data/Random/Distribution/Dirichlet.hs @@ -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] diff --git a/src/Data/Random/Distribution/Gamma.hs b/src/Data/Random/Distribution/Gamma.hs --- a/src/Data/Random/Distribution/Gamma.hs +++ b/src/Data/Random/Distribution/Gamma.hs @@ -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) - diff --git a/src/Data/Random/Distribution/Multinomial.hs b/src/Data/Random/Distribution/Multinomial.hs --- a/src/Data/Random/Distribution/Multinomial.hs +++ b/src/Data/Random/Distribution/Multinomial.hs @@ -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 diff --git a/src/Data/Random/Distribution/Normal.hs b/src/Data/Random/Distribution/Normal.hs --- a/src/Data/Random/Distribution/Normal.hs +++ b/src/Data/Random/Distribution/Normal.hs @@ -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) diff --git a/src/Data/Random/Distribution/Poisson.hs b/src/Data/Random/Distribution/Poisson.hs --- a/src/Data/Random/Distribution/Poisson.hs +++ b/src/Data/Random/Distribution/Poisson.hs @@ -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 diff --git a/src/Data/Random/Distribution/Rayleigh.hs b/src/Data/Random/Distribution/Rayleigh.hs --- a/src/Data/Random/Distribution/Rayleigh.hs +++ b/src/Data/Random/Distribution/Rayleigh.hs @@ -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 diff --git a/src/Data/Random/Distribution/Triangular.hs b/src/Data/Random/Distribution/Triangular.hs --- a/src/Data/Random/Distribution/Triangular.hs +++ b/src/Data/Random/Distribution/Triangular.hs @@ -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 diff --git a/src/Data/Random/Distribution/Uniform.hs b/src/Data/Random/Distribution/Uniform.hs --- a/src/Data/Random/Distribution/Uniform.hs +++ b/src/Data/Random/Distribution/Uniform.hs @@ -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 diff --git a/src/Data/Random/Distribution/Weibull.hs b/src/Data/Random/Distribution/Weibull.hs --- a/src/Data/Random/Distribution/Weibull.hs +++ b/src/Data/Random/Distribution/Weibull.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, FlexibleContexts #-} module Data.Random.Distribution.Weibull where import Data.Random.Distribution diff --git a/src/Data/Random/Distribution/Ziggurat.hs b/src/Data/Random/Distribution/Ziggurat.hs --- a/src/Data/Random/Distribution/Ziggurat.hs +++ b/src/Data/Random/Distribution/Ziggurat.hs @@ -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 diff --git a/src/Data/Random/Internal/Find.hs b/src/Data/Random/Internal/Find.hs --- a/src/Data/Random/Internal/Find.hs +++ b/src/Data/Random/Internal/Find.hs @@ -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 diff --git a/src/Data/Random/Internal/Fixed.hs b/src/Data/Random/Internal/Fixed.hs --- a/src/Data/Random/Internal/Fixed.hs +++ b/src/Data/Random/Internal/Fixed.hs @@ -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 diff --git a/src/Data/Random/Internal/TH.hs b/src/Data/Random/Internal/TH.hs deleted file mode 100644 --- a/src/Data/Random/Internal/TH.hs +++ /dev/null @@ -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 - ] - diff --git a/src/Data/Random/Lift.hs b/src/Data/Random/Lift.hs --- a/src/Data/Random/Lift.hs +++ b/src/Data/Random/Lift.hs @@ -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 - diff --git a/src/Data/Random/List.hs b/src/Data/Random/List.hs --- a/src/Data/Random/List.hs +++ b/src/Data/Random/List.hs @@ -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" - diff --git a/src/Data/Random/RVar.hs b/src/Data/Random/RVar.hs --- a/src/Data/Random/RVar.hs +++ b/src/Data/Random/RVar.hs @@ -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 diff --git a/src/Data/Random/Sample.hs b/src/Data/Random/Sample.hs --- a/src/Data/Random/Sample.hs +++ b/src/Data/Random/Sample.hs @@ -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)