random-fu 0.2.6.2 → 0.2.7.0
raw patch · 5 files changed
+89/−33 lines, 5 files
Files
- changelog.md +11/−8
- random-fu.cabal +3/−2
- src/Data/Random/Distribution/Beta.hs +2/−2
- src/Data/Random/Distribution/Binomial.hs +23/−21
- src/Data/Random/Distribution/Simplex.hs +50/−0
changelog.md view
@@ -1,3 +1,14 @@+* Changes in 0.2.7.0: Add Simplex, fix logBetaPdf, fix binomialPdf and+ binomialCdf to actually use the numerically stable method!++* 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.+ * Changes in 0.2.4.0: Added a Lift instance that resolves a common overlapping-instance issue in user code. @@ -13,11 +24,3 @@ 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.
random-fu.cabal view
@@ -1,5 +1,5 @@ name: random-fu-version: 0.2.6.2+version: 0.2.7.0 stability: provisional cabal-version: >= 1.6@@ -29,7 +29,7 @@ a fair bit slower than straight C implementations of the same algorithms. -tested-with: GHC == 7.4.2, GHC == 7.6.1, GHC == 7.8.3+tested-with: GHC == 7.10.3 extra-source-files: changelog.md @@ -63,6 +63,7 @@ Data.Random.Distribution.Pareto Data.Random.Distribution.Poisson Data.Random.Distribution.Rayleigh+ Data.Random.Distribution.Simplex Data.Random.Distribution.T Data.Random.Distribution.Triangular Data.Random.Distribution.Uniform
src/Data/Random/Distribution/Beta.hs view
@@ -41,8 +41,8 @@ logBetaPdf :: Double -> Double -> Double -> Double logBetaPdf a b x | a <= 0 || b <= 0 = nan- | x <= 0 = 0- | x >= 1 = 0+ | x <= 0 = log 0+ | x >= 1 = log 0 | otherwise = (a-1)*log x + (b-1)*log (1-x) - logBeta a b where nan = 0.0 / 0.0
src/Data/Random/Distribution/Binomial.hs view
@@ -20,7 +20,7 @@ -- 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) + -- note that although it's fast enough for large (eg, 2^10000) -- @Integer@s, it's not accurate enough when using @Double@ as -- the @b@ parameter. integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVarT m a@@ -31,12 +31,12 @@ | t > 10 = do let a = 1 + t `div` 2 b = 1 + t - a- + x <- betaT (fromIntegral a) (fromIntegral b) if x >= p then bin k (a - 1) (p / x) else bin (k + a) (b - 1) ((p - x) / (1 - x))- + | otherwise = count k t where count !k' 0 = return k'@@ -45,26 +45,28 @@ count (if x < p then k' + 1 else k') (n-1) count _ _ = error "integralBinomial: negative number of trials specified" --- TODO: improve performance integralBinomialCDF :: (Integral a, Real b) => a -> b -> a -> Double-integralBinomialCDF t p x = sum- [ fromInteger (toInteger t `c` toInteger i) * p' ^^ i * (1-p') ^^ (t-i)- | i <- [0 .. x]- ]- - where - p' = realToFrac p- n `c` k = product [n-k+1..n] `div` product [1..k]+integralBinomialCDF t p x = sum $ map (integralBinomialPDF t p) $ [0 .. x] --- TODO: improve performance and re-use in CDF+-- | The probability of getting exactly k successes in n trials is+-- given by the probability mass function:+--+-- \[+-- f(k;n,p) = \Pr(X = k) = \binom n k p^k(1-p)^{n-k}+-- \]+--+-- Note that in `integralBinomialPDF` the parameters of the mass+-- function are given first and the range of the random variable+-- distributed according to the binomial distribution is given+-- last. That is, \(f(2;4,0.5)\) is calculated by @integralBinomialPDF 4 0.5 2@.+ 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]+ exp $ integralBinomialLogPdf t p x +-- | We use the method given in \"Fast and accurate computation of+-- binomial probabilities, Loader, C\",+-- <http://octave.1599824.n4.nabble.com/attachment/3829107/0/loader2000Fast.pdf> integralBinomialLogPdf :: (Integral a, Real b) => a -> b -> a -> Double integralBinomialLogPdf nI pR xI | p == 0.0 && xI == 0 = 1.0@@ -84,7 +86,7 @@ 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 #-}@@ -132,7 +134,7 @@ , Distribution Beta b , Distribution StdUniform b ) => Distribution (Binomial b) Int- where + where rvarT (Binomial t p) = integralBinomial t p instance ( Real b , Distribution (Binomial b) Int ) => CDF (Binomial b) Int@@ -144,7 +146,7 @@ |]) $( replicateInstances ''Float realFloatTypes [d|- instance Distribution (Binomial b) Integer + instance Distribution (Binomial b) Integer => Distribution (Binomial b) Float where rvar (Binomial t p) = floatingBinomial t p instance CDF (Binomial b) Integer
+ src/Data/Random/Distribution/Simplex.hs view
@@ -0,0 +1,50 @@+{-# LANGUAGE+ MultiParamTypeClasses,+ FlexibleContexts, FlexibleInstances,+ UndecidableInstances, GADTs+ #-}++module Data.Random.Distribution.Simplex+ ( StdSimplex(..)+ , stdSimplex+ , stdSimplexT+ , fractionalStdSimplex+ ) where++import Control.Applicative+import Control.Monad+import Data.List+import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++-- |Uniform distribution over a standard simplex.+newtype StdSimplex as =+ -- | @StdSimplex k@ constructs a standard simplex of dimension @k@+ -- (standard /k/-simplex).+ -- An element of the simplex represents a vector variable @as = (a_0,+ -- a_1, ..., a_k)@. The elements of @as@ are more than or equal to @0@+ -- and @sum as@ is always equal to @1@.+ StdSimplex Int+ deriving (Eq, Show)++instance (Ord a, Fractional a, Distribution StdUniform a) => Distribution StdSimplex [a] where+ rvar (StdSimplex k) = fractionalStdSimplex k++-- |@stdSimplex k@ returns a random variable being uniformly distributed over+-- a standard simplex of dimension @k@.+stdSimplex :: Distribution StdSimplex [a] => Int -> RVar [a]+stdSimplex k = rvar (StdSimplex k)++stdSimplexT :: Distribution StdSimplex [a] => Int -> RVarT m [a]+stdSimplexT k = rvarT (StdSimplex k)++-- |An algorithm proposed by Rubinstein & Melamed (1998).+-- See, /e.g./, S. Onn, I. Weissman.+-- Generating uniform random vectors over a simplex with implications to+-- the volume of a certain polytope and to multivariate extremes.+-- /Ann Oper Res/ (2011) __189__:331-342.+fractionalStdSimplex :: (Ord a, Fractional a, Distribution StdUniform a) => Int -> RVar [a]+fractionalStdSimplex k = do us <- replicateM k stdUniform+ let us' = sort us ++ [1]+ return $ zipWith (-) us' (0 : us')