diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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. 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.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 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 @@ -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 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 @@ -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\", +-- 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 diff --git a/src/Data/Random/Distribution/Simplex.hs b/src/Data/Random/Distribution/Simplex.hs new file mode 100644 --- /dev/null +++ b/src/Data/Random/Distribution/Simplex.hs @@ -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')