diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,7 @@ +* Changes in 0.2.7.3: Remove dependence on log-domain. Raise lower bound for base to 4.9. + +* Changes in 0.2.7.1: Add PDF instance for Poisson. + * Changes in 0.2.7.0: Add Simplex, fix logBetaPdf, fix binomialPdf and binomialCdf to actually use the numerically stable method! 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.0 +version: 0.2.7.3 stability: provisional cabal-version: >= 1.6 @@ -78,7 +78,7 @@ Data.Random.Sample Data.Random.Vector if flag(base4_2) - build-depends: base >= 4.2 && <5 + build-depends: base >= 4.9 && <5 else cpp-options: -Dold_Fixed build-depends: base >= 4 && <4.2 @@ -98,13 +98,7 @@ template-haskell, transformers, vector >= 0.7, - log-domain >=0.9 && <1.0 - - if os(Windows) - cpp-options: -Dwindows - build-depends: erf-native - else - build-depends: erf + erf if impl(ghc == 7.2.1) -- Doesn't work under GHC 7.2.1 due to 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 @@ -16,7 +16,7 @@ import Numeric.SpecFunctions ( stirlingError ) import Numeric.SpecFunctions.Extra ( bd0 ) -import Numeric.Log ( log1p ) +import Numeric ( log1p ) -- algorithm from Knuth's TAOCP, 3rd ed., p 136 -- specific choice of cutoff size taken from gsl source 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 @@ -41,13 +41,13 @@ categoricalT :: (Num p, Distribution (Categorical p) a) => [(p,a)] -> RVarT m a categoricalT = rvarT . fromList --- |Construct a 'Categorical' random variable from a list of probabilities --- and categories, where the probabilities all sum to 1. +-- |Construct a 'Categorical' random variable from a list of weights +-- and categories. The weights do /not/ have to sum to 1. 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 probabilities --- and categories, where the probabilities all sum to 1. +-- |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 diff --git a/src/Data/Random/Distribution/Exponential.hs b/src/Data/Random/Distribution/Exponential.hs --- a/src/Data/Random/Distribution/Exponential.hs +++ b/src/Data/Random/Distribution/Exponential.hs @@ -10,6 +10,15 @@ import Data.Random.Distribution import Data.Random.Distribution.Uniform +{-| +A definition of the exponential distribution over the type @a@. + +@'Exp' mu@ models an exponential distribution with mean @mu@. This can +alternatively be viewed as an exponential distribution with parameter @lambda = +1 / mu@. + +See also 'exponential'. +-} newtype Exponential a = Exp a floatingExponential :: (Floating a, Distribution StdUniform a) => a -> RVarT m a @@ -20,9 +29,23 @@ floatingExponentialCDF :: Real a => a -> a -> Double floatingExponentialCDF lambdaRecip x = 1 - exp (negate (realToFrac x) / realToFrac lambdaRecip) +{-| +A random variable which samples from the exponential distribution. + +@'exponential' mu@ is an exponential random variable with mean @mu@. This can +alternatively be viewed as an exponential random variable with parameter @lambda += 1 / mu@. +-} exponential :: Distribution Exponential a => a -> RVar a exponential = rvar . Exp +{-| +A random variable transformer which samples from the exponential distribution. + +@'exponentialT' mu@ is an exponential random variable with mean @mu@. This can +alternatively be viewed as an exponential random variable with parameter @lambda += 1 / mu@. +-} exponentialT :: Distribution Exponential a => a -> RVarT m a exponentialT = rvarT . Exp 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 @@ -24,18 +24,18 @@ psn j mu | mu > 10 = do let m = floor (mu * (7/8)) - + x <- erlangT m if x >= mu then do b <- binomialT (m - 1) (mu / x) return (j + b) else psn (j + m) (mu - x) - + | otherwise = prod 1 j where emu = exp (-mu) - + prod p k = do u <- stdUniformT if p * u > emu @@ -47,15 +47,36 @@ [ exp (fromIntegral i * log lambda - i_fac_ln) | (i, i_fac_ln) <- zip [0..k] (scanl (+) 0 (map log [1..])) ] - + where lambda = realToFrac mu +-- | The probability of getting exactly k successes is +-- given by the probability mass function: +-- +-- \[ +-- f(k;\lambda) = \Pr(X = k) = \frac{\lambda^k e^{-\lambda}}{k!} +-- \] +-- +-- Note that in `integralPoissonPDF` the parameter of the mass +-- function are given first and the range of the random variable +-- distributed according to the Poisson distribution is given +-- last. That is, \(f(2;0.5)\) is calculated by @integralPoissonPDF 0.5 2@. +integralPoissonPDF :: (Integral a, Real b) => b -> a -> Double +integralPoissonPDF mu k = exp (negate lambda) * + exp (fromIntegral k * log lambda - k_fac_ln) + where + k_fac_ln = foldl (+) 0 (map (log . fromIntegral) [1..k]) + lambda = realToFrac mu + fractionalPoisson :: (Num a, Distribution (Poisson b) Integer) => b -> RVarT m a fractionalPoisson mu = liftM fromInteger (poissonT mu) fractionalPoissonCDF :: (CDF (Poisson b) Integer, RealFrac a) => b -> a -> Double fractionalPoissonCDF mu k = cdf (Poisson mu) (floor k :: Integer) +fractionalPoissonPDF :: (PDF (Poisson b) Integer, RealFrac a) => b -> a -> Double +fractionalPoissonPDF mu k = pdf (Poisson mu) (floor k :: Integer) + poisson :: (Distribution (Poisson b) a) => b -> RVar a poisson mu = rvar (Poisson mu) @@ -65,7 +86,7 @@ newtype Poisson b a = Poisson b $( replicateInstances ''Int integralTypes [d| - instance ( RealFloat b + instance ( RealFloat b , Distribution StdUniform b , Distribution (Erlang Int) b , Distribution (Binomial b) Int @@ -73,6 +94,8 @@ 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 |] ) $( replicateInstances ''Float realFloatTypes [d| @@ -80,4 +103,6 @@ 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 |])