random-fu 0.2.7.0 → 0.2.7.3
raw patch · 6 files changed
+65/−19 lines, 6 filesdep −erf-nativedep −log-domaindep ~base
Dependencies removed: erf-native, log-domain
Dependency ranges changed: base
Files
- changelog.md +4/−0
- random-fu.cabal +3/−9
- src/Data/Random/Distribution/Binomial.hs +1/−1
- src/Data/Random/Distribution/Categorical.hs +4/−4
- src/Data/Random/Distribution/Exponential.hs +23/−0
- src/Data/Random/Distribution/Poisson.hs +30/−5
changelog.md view
@@ -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!
random-fu.cabal view
@@ -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
src/Data/Random/Distribution/Binomial.hs view
@@ -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
src/Data/Random/Distribution/Categorical.hs view
@@ -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
src/Data/Random/Distribution/Exponential.hs view
@@ -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
src/Data/Random/Distribution/Poisson.hs view
@@ -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 |])