module Statistics.Distribution.Binomial
(
BinomialDistribution
, binomial
, bdTrials
, bdProbability
) where
import Data.Binary (Binary)
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D
import qualified Statistics.Distribution.Poisson.Internal as I
import Numeric.SpecFunctions (choose,incompleteBeta)
import Numeric.MathFunctions.Constants (m_epsilon)
import Data.Binary (put, get)
import Control.Applicative ((<$>), (<*>))
data BinomialDistribution = BD {
bdTrials :: !Int
, bdProbability :: !Double
} deriving (Eq, Read, Show, Typeable, Data, Generic)
instance Binary BinomialDistribution where
put (BD x y) = put x >> put y
get = BD <$> get <*> get
instance D.Distribution BinomialDistribution where
cumulative = cumulative
instance D.DiscreteDistr BinomialDistribution where
probability = probability
instance D.Mean BinomialDistribution where
mean = mean
instance D.Variance BinomialDistribution where
variance = variance
instance D.MaybeMean BinomialDistribution where
maybeMean = Just . D.mean
instance D.MaybeVariance BinomialDistribution where
maybeStdDev = Just . D.stdDev
maybeVariance = Just . D.variance
instance D.Entropy BinomialDistribution where
entropy (BD n p)
| n == 0 = 0
| n <= 100 = directEntropy (BD n p)
| otherwise = I.poissonEntropy (fromIntegral n * p)
instance D.MaybeEntropy BinomialDistribution where
maybeEntropy = Just . D.entropy
probability :: BinomialDistribution -> Int -> Double
probability (BD n p) k
| k < 0 || k > n = 0
| n == 0 = 1
| otherwise = choose n k * p^k * (1p)^(nk)
cumulative :: BinomialDistribution -> Double -> Double
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
| otherwise = incompleteBeta (fromIntegral (nk)) (fromIntegral (k+1)) (1 p)
where
k = floor x
mean :: BinomialDistribution -> Double
mean (BD n p) = fromIntegral n * p
variance :: BinomialDistribution -> Double
variance (BD n p) = fromIntegral n * p * (1 p)
directEntropy :: BinomialDistribution -> Double
directEntropy d@(BD n _) =
negate . sum $
takeWhile (< negate m_epsilon) $
dropWhile (not . (< negate m_epsilon)) $
[ let x = probability d k in x * log x | k <- [0..n]]
binomial :: Int
-> Double
-> BinomialDistribution
binomial n p
| n < 0 =
error $ msg ++ "number of trials must be non-negative. Got " ++ show n
| p < 0 || p > 1 =
error $ msg++"probability must be in [0,1] range. Got " ++ show p
| otherwise = BD n p
where msg = "Statistics.Distribution.Binomial.binomial: "