random-fu 0.2.6.0 → 0.2.6.1

raw patch · 8 files changed

+132/−18 lines, 8 filesdep +logfloatdep ~basenew-uploaderPVP: minor bump suggested

API additions: PVP suggests at least a minor version bump

Dependencies added: logfloat

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• changelog.md +23/−0
• random-fu.cabal +5/−17
• src/Data/Random.hs +1/−1
• src/Data/Random/Distribution.hs +7/−0
• src/Data/Random/Distribution/Beta.hs +20/−0
• src/Data/Random/Distribution/Binomial.hs +47/−0
• src/Data/Random/Distribution/Normal.hs +19/−0
• src/Data/Random/Distribution/Uniform.hs +10/−0

@@ -0,0 +1,23 @@+* Changes in 0.2.4.0: Added a Lift instance that resolves a common+  overlapping-instance issue in user code.++* Changes in 0.2.3.1: Should now build on GHC 7.6++* Changes in 0.2.3.0: Added stretched exponential distribution,+  contributed by Ben Gamari.++* Changes in 0.2.2.0: Bug fixes in+  Data.Random.Distribution.Categorical.++* Changes in 0.2.1.1: Changed some one-field data types to newtypes,+  updated types for GHC 7.4's removal of Eq and Show from the context+  of Num, and added RVarT versions of random variables in+  Data.Random.List++* Changes in 2.6.1: now supports probability density functions and log+  probability density functions via the PDF class, similar to R and+  initially just for the Beta, Binomial, Normal and Uniform+  distributions. The log Binomial probability density function uses+  *Fast and Accurate Computation of Binomial Probabilities* by+  Catherine Loader (this is what is implemented in R and Octave) to+  minimize the occurrence of underflow.

@@ -1,5 +1,5 @@ name:                   random-fu-version:                0.2.6.0+version:                0.2.6.1 stability:              provisional  cabal-version:          >= 1.6@@ -28,24 +28,11 @@                         comparable to other Haskell libraries, but still                         a fair bit slower than straight C implementations of                          the same algorithms.-                        .-                        Changes in 0.2.4.0: Added a Lift instance that resolves-                        a common overlapping-instance issue in user code.-                        .-                        Changes in 0.2.3.1: Should now build on GHC 7.6-                        .-                        Changes in 0.2.3.0: Added stretched exponential distribution,-                        contributed by Ben Gamari.-                        .-                        Changes in 0.2.2.0: Bug fixes in Data.Random.Distribution.Categorical.-                        .-                        Changes in 0.2.1.1: Changed some one-field data types-                        to newtypes, updated types for GHC 7.4's removal of Eq -                        and Show from the context of Num, and added RVarT versions-                        of random variables in Data.Random.List  tested-with:            GHC == 7.4.2, GHC == 7.6.1 +extra-source-files:     changelog.md+ source-repository head   type:                 git   location:             https://github.com/mokus0/random-fu.git@@ -109,7 +96,8 @@                         syb,                         template-haskell,                         transformers,-                        vector >= 0.7+                        vector >= 0.7,+                        logfloat >=0.12 && <0.13    if os(Windows)     cpp-options:        -Dwindows

@@ -42,7 +42,7 @@       runRVar, runRVarT, runRVarTWith,        -- ** Concrete ('Distribution')-      Distribution(..), CDF(..),+      Distribution(..), CDF(..), PDF(..),              -- * Sampling random variables       Sampleable(..), sample, sampleState, sampleStateT,

@@ -60,6 +60,13 @@     rvarT :: d t -> RVarT n t     rvarT d = lift (rvar d) +-- FIXME: I am not sure about giving default instances+class Distribution d t => PDF d t where+    pdf :: d t -> t -> Double+    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.

@@ -14,6 +14,8 @@ import Data.Random.Distribution.Gamma import Data.Random.Distribution.Uniform +import Numeric.SpecFunctions+ {-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVarT m Float #-} {-# SPECIALIZE fractionalBeta :: Double -> Double -> RVarT m Double #-} fractionalBeta :: (Fractional a, Eq a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVarT m a@@ -34,6 +36,24 @@ betaT a b = rvarT (Beta a b)  data Beta a = Beta a a++-- FIXME: I am far from convinced that NaNs are a good idea.+logBetaPdf :: Double -> Double -> Double -> Double+logBetaPdf a b x+   | a <= 0 || b <= 0 = nan+   | x <= 0 = 0+   | x >= 1 = 0+   | otherwise = (a-1)*log x + (b-1)*log (1-x) - logBeta a b+  where+    nan = 0.0 / 0.0++instance PDF Beta Double+  where+    pdf (Beta a b) = exp . logBetaPdf a b++instance PDF Beta Float+  where+    pdf (Beta a b) = realToFrac . exp . logBetaPdf (realToFrac a) (realToFrac b) . realToFrac  $( replicateInstances ''Float realFloatTypes [d|         instance Distribution Beta Float

@@ -14,6 +14,10 @@ import Data.Random.Distribution.Beta import Data.Random.Distribution.Uniform +import Numeric.SpecFunctions ( stirlingError )+import Numeric.SpecFunctions.Extra ( bd0 )+import Data.Number.LogFloat ( log1p )+     -- algorithm from Knuth's TAOCP, 3rd ed., p 136     -- specific choice of cutoff size taken from gsl source     -- note that although it's fast enough for large (eg, 2^10000) @@ -52,6 +56,35 @@         p' = realToFrac p         n `c` k = product [n-k+1..n] `div` product [1..k] +-- TODO: improve performance and re-use in CDF+integralBinomialPDF :: (Integral a, Real b) => a -> b -> a -> Double+integralBinomialPDF t p x =+    fromInteger (toInteger t `c` toInteger x) * p' ^^ x * (1-p') ^^ (t-x)+    +    where +        p' = realToFrac p+        n `c` k = product [n-k+1..n] `div` product [1..k]++integralBinomialLogPdf :: (Integral a, Real b) => a -> b -> a -> Double+integralBinomialLogPdf nI pR xI+  | p == 0.0 && xI == 0   = 1.0+  | p == 0.0              = 0.0+  | p == 1.0 && xI == nI  = 1.0+  | p == 1.0              = 0.0+  |             xI == 0   = n * log (1-p)+  |             xI == nI  = n * log p+  | otherwise = lc - 0.5 * lf+  where+    n = fromIntegral nI+    x = fromIntegral xI+    p = realToFrac pR+    lc = stirlingError n -+         stirlingError x -+         stirlingError (n - x) -+         bd0 x (n * p) -+         bd0 (n - x) (n * (1 - p))+    lf = log (2 * pi) + log x + log1p (- x / n)+   -- would it be valid to repeat the above computation using fractional @t@? -- obviously something different would have to be done with @count@ as well... {-# SPECIALIZE floatingBinomial :: Float  -> Float  -> RVar Float  #-}@@ -64,6 +97,12 @@ floatingBinomialCDF :: (CDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double floatingBinomialCDF t p x = cdf (Binomial (truncate t :: Integer) p) (floor x) +floatingBinomialPDF :: (PDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double+floatingBinomialPDF t p x = pdf (Binomial (truncate t :: Integer) p) (floor x)++floatingBinomialLogPDF :: (PDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double+floatingBinomialLogPDF t p x = logPdf (Binomial (truncate t :: Integer) p) (floor x)+ {-# SPECIALIZE binomial :: Int     -> Float  -> RVar Int #-} {-# SPECIALIZE binomial :: Int     -> Double -> RVar Int #-} {-# SPECIALIZE binomial :: Integer -> Float  -> RVar Integer #-}@@ -98,6 +137,10 @@         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     |])  $( replicateInstances ''Float realFloatTypes [d|@@ -107,4 +150,8 @@         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     |])

@@ -197,6 +197,19 @@ normalCdf :: (Real a) => a -> a -> a -> Double normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s) +normalPdf :: (Real a, Floating b) => a -> a -> a -> b+normalPdf mu sigma x =+  (recip (sqrt (2 * pi * sigma2))) * (exp ((-((realToFrac x) - (realToFrac mu))^2) / (2 * sigma2)))+  where+    sigma2 = realToFrac sigma^2++normalLogPdf :: (Real a, Floating b) => a -> a -> a -> b+normalLogPdf mu sigma x =+  log (recip (sqrt (2 * pi * sigma2))) ++  ((-((realToFrac x) - (realToFrac mu))^2) / (2 * sigma2))+  where+    sigma2 = realToFrac sigma^2+ -- |A specification of a normal distribution over the type 'a'. data Normal a     -- |The \"standard\" normal distribution - mean 0, stddev 1@@ -219,6 +232,12 @@ instance (Real a, Distribution Normal a) => CDF Normal a where     cdf StdNormal    = normalCdf 0 1     cdf (Normal m s) = normalCdf m s++instance (Real a, Floating a, Distribution Normal a) => PDF Normal a where+  pdf StdNormal    = normalPdf 0 1+  pdf (Normal m s) = normalPdf m s+  logPdf StdNormal = normalLogPdf 0 1+  logPdf (Normal m s) = normalLogPdf m s  {-# SPECIALIZE stdNormal :: RVar Double #-} {-# SPECIALIZE stdNormal :: RVar Float #-}

@@ -157,6 +157,13 @@     | x >= 1    = 1     | otherwise = realToFrac x +-- |The PDF of the random variable 'realFloatStdUniform'.+realStdUniformPDF :: Real a => a -> Double+realStdUniformPDF x+    | x <= 0    = 0+    | x >= 1    = 0+    | otherwise = 1+ -- |(internal) basic linear interpolation; @lerp x y@ is a linear function whose -- value is @x@ at 0 and @y@ at 1 lerp :: Num a => a -> a -> a -> a@@ -327,6 +334,9 @@ instance Distribution StdUniform Double     where rvarT _ = getRandomDouble instance CDF StdUniform Float               where cdf   _ = realStdUniformCDF instance CDF StdUniform Double              where cdf   _ = realStdUniformCDF+instance PDF StdUniform Float               where pdf   _ = realStdUniformPDF+instance PDF StdUniform Double              where pdf   _ = realStdUniformPDF+  instance HasResolution r =>           Distribution Uniform (Fixed r)     where rvarT (Uniform a b) = fixedUniform  a b