Commits

Bryan O'Sullivan committed c0aca4b Merge

Merge

Comments (0)

Files changed (2)

Statistics/Distribution/Binomial.hs

 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
-import Numeric.SpecFunctions (choose)
+import Numeric.SpecFunctions (choose,incompleteBeta)
 
 
 -- | The binomial distribution.
 
 -- Summation from different sides required to reduce roundoff errors
 cumulative :: BinomialDistribution -> Double -> Double
-cumulative d@(BD n _) x
+cumulative (BD n p) x
   | isNaN x      = error "Statistics.Distribution.Binomial.cumulative: NaN input"
   | isInfinite x = if x > 0 then 1 else 0
   | k <  0       = 0
   | k >= n       = 1
-  | k <  m       = D.sumProbabilities d 0 k
-  | otherwise    = 1 - D.sumProbabilities d (k+1) n
+  | otherwise    = incompleteBeta (fromIntegral (n-k)) (fromIntegral (k+1)) (1 - p)
   where
-    m = floor (mean d)
     k = floor x
 {-# INLINE cumulative #-}
 

tests/Tests/Distribution.hs

 -- CDF for discrete distribution uses <= for comparison
 cdfDiscreteIsCorrect :: (DiscreteDistr d) => T d -> d -> Property
 cdfDiscreteIsCorrect _ d
-  = printTestCase (unlines $ map show badN)
+  = printTestCase (unlines badN)
   $ null badN  
   where
     -- We are checking that:
     -- Apporixmate equality is tricky here. Scale is set by maximum
     -- value of CDF and probability. Case when all proabilities are
     -- zero should be trated specially.
-    badN = [ (i,p,p1,dp, (p1-p-dp) / max p1 dp)
+    badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)
            | i <- [0 .. 100]
            , let p      = cumulative d $ fromIntegral i - 1e-6
                  p1     = cumulative d $ fromIntegral i