statistics 0.3.6 → 0.4.0

raw patch · 5 files changed

+75/−15 lines, 5 files

Files

• Statistics/Distribution/Binomial.hs +66/−8
• Statistics/Distribution/Hypergeometric.hs +1/−1
• Statistics/Math.hs +5/−3
• Statistics/Sample/Powers.hs +1/−1
• statistics.cabal +2/−2

@@ -27,8 +27,10 @@ import Data.Array.Vector import Data.Int (Int64) import Data.Typeable (Typeable)+import Statistics.Constants (m_epsilon) import qualified Statistics.Distribution as D-import Statistics.Math (choose)+import Statistics.Distribution.Normal (standard)+import Statistics.Math (choose, logFactorial)  -- | The binomial distribution. data BinomialDistribution = BD {@@ -50,17 +52,69 @@     mean = mean  density :: BinomialDistribution -> Double -> Double-density (BD n p) x =-    fromIntegral (n `choose` floor x) * p ** x * (1-p) ** (fromIntegral n-x)+density (BD n p) x+    | not (isIntegral x) = integralError "density"+    | n == 0             = 1+    | x < 0 || x > n'    = 0+    | n <= 50 || x < 2   = sign * p'' ** x' * (n `choose` fx) * q'' ** nx'+    | otherwise          = sign * exp (x' * log p'' + nx' * log q'' + lf)+  where sign = oddX * oddNX+        (x',p',q') | x > n' / 2 = (n'-x, q, p)+                   | otherwise  = (x,    p, q)+        oddX | p' < 0 && odd fx     = -1+             | otherwise            = 1+        oddNX | q' < 0 && odd nx    = -1+              | otherwise           = 1+        p'' = abs p'+        q'' = abs q'+        q   = 1 - p+        nx  = n - fx+        nx' = fromIntegral nx+        fx  = floor x'+        n'  = fromIntegral n+        lf  = logFactorial n - logFactorial nx - logFactorial fx  cumulative :: BinomialDistribution -> Double -> Double-cumulative d =-    sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor+cumulative d x+  | isIntegral x = sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor $ x+  | otherwise    = integralError "cumulative" +isIntegral :: Double -> Bool+isIntegral x = x == floorf x++floorf :: Double -> Double+floorf = fromIntegral . (floor :: Double -> Int64)+ quantile :: BinomialDistribution -> Double -> Double-quantile d@(BD n _p) p = fromIntegral . r64 $ D.findRoot d p (n'/2) 0 n'-    where n'  = fromIntegral n-          r64 = round :: Double -> Int64+quantile dist@(BD n p) prob+    | isNaN prob = prob+    | p == 1     = n'+    | n' < 1e5   = fst (search 1 y0 z0)+    | otherwise  = let dy = floorf (n' / 1000)+                   in  narrow dy (search dy y0 z0)+  where q  = 1 - p+        n' = fromIntegral n+        y0 = n' `min` floorf (µ + σ * (d + γ * (d * d - 1) / 6) + 0.5)+          where µ  = n' * p+                σ  = sqrt (n' * p * q)+                d = D.quantile standard prob+                γ  = (q - p) / σ+        z0 = cumulative dist y0+        search dy y1 z1 | z0 >= prob' = left y1 z1+                        | otherwise   = right y1+          where+            prob' = prob * (1 - 64 * m_epsilon)+            left y oldZ | y == 0 || z < prob' = (y, oldZ)+                        | otherwise           = left (max 0 y') z+                where z  = cumulative dist y'+                      y' = y - dy+            right y | y' >= n' || z >= prob' = (y', z)+                    | otherwise              = right y'+                where z  = cumulative dist y'+                      y' = y + dy+        narrow dy (y,z) | dy <= 1 || dy' <= n'/1e15 = y+                        | otherwise                 = narrow dy' (search dy y z)+            where dy' = floorf (dy / 100)  mean :: BinomialDistribution -> Double mean (BD n p) = fromIntegral n * p@@ -78,3 +132,7 @@     assert (p > 0 && p < 1) $     BD n p {-# INLINE binomial #-}++integralError :: String -> a+integralError f = error ("Statistics.Distribution.Binomial." ++ f +++                         ": non-integer-valued input")

@@ -80,7 +80,7 @@     | r > maxVal = 1/0     | otherwise  = exp r   where-    a <> b = fromIntegral (a `choose` b)+    a <> b = a `choose` b     r = f m + f (l-m) - f l - f xi - f (k-xi) + f k -         f (m-xi) - f (l-m-k+xi) + f (l-k)     f = logFactorial

@@ -49,15 +49,17 @@ -- | The binomial coefficient. -- -- > 7 `choose` 3 == 35-choose :: Int -> Int -> Int+choose :: Int -> Int -> Double n `choose` k-    | k > n = 0-    | otherwise = ceiling . foldlU go 1 . enumFromToU 1 $ k'+    | k > n     = 0+    | k < 30    = foldlU go 1 . enumFromToU 1 $ k'+    | otherwise = exp $ lg (n+1) - lg (k+1) - lg (n-k+1)     where go a i = a * (nk + j) / j               where j = fromIntegral i :: Double           k' | k > n `div` 2 = n - k              | otherwise     = k           nk = fromIntegral (n - k')+          lg = logGamma . fromIntegral {-# INLINE choose #-}  data F = F {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64

@@ -107,7 +107,7 @@     | k == 0    = 1     | otherwise = (/n) . sumU . mapU go . indexedU . takeU (k+1) $ pa   where-    go (i :*: e) = fromIntegral (k `choose` i) * ((-m) ^ (k-i)) * e+    go (i :*: e) = (k `choose` i) * ((-m) ^ (k-i)) * e     n = indexU pa 0     m = mean p {-# INLINE centralMoment #-}

@@ -1,5 +1,5 @@ name:           statistics-version:        0.3.6+version:        0.4.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -65,6 +65,6 @@   -- gather extensive profiling data for now   ghc-prof-options: -auto-all -  ghc-options: -Wall -funbox-strict-fields -O2+  ghc-options: -Wall -funbox-strict-fields   if impl(ghc >= 6.8)     ghc-options: -fwarn-tabs