diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -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. 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.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 diff --git a/src/Data/Random.hs b/src/Data/Random.hs --- a/src/Data/Random.hs +++ b/src/Data/Random.hs @@ -42,7 +42,7 @@ runRVar, runRVarT, runRVarTWith, -- ** Concrete ('Distribution') - Distribution(..), CDF(..), + Distribution(..), CDF(..), PDF(..), -- * Sampling random variables Sampleable(..), sample, sampleState, sampleStateT, 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 @@ -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. 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 @@ -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 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 @@ -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 |]) 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 @@ -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 #-} 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 @@ -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