module Statistics.Distribution.Hypergeometric
(
HypergeometricDistribution
, fromParams
, hdM
, hdL
, hdK
) where
import Control.Exception (assert)
import Data.Array.Vector
import Data.Typeable (Typeable)
import Statistics.Math (choose, logFactorial)
import Statistics.Constants (m_max_exp)
import qualified Statistics.Distribution as D
data HypergeometricDistribution = HD {
hdM :: !Int
, hdL :: !Int
, hdK :: !Int
} deriving (Eq, Read, Show, Typeable)
instance D.Distribution HypergeometricDistribution where
probability = probability
cumulative = cumulative
inverse = inverse
instance D.Variance HypergeometricDistribution where
variance = variance
instance D.Mean HypergeometricDistribution where
mean = mean
variance :: HypergeometricDistribution -> Double
variance (HD m l k) = (k' * ml) * (1 ml) * (l' k') / (l' 1)
where m' = fromIntegral m
l' = fromIntegral l
k' = fromIntegral k
ml = m' / l'
mean :: HypergeometricDistribution -> Double
mean (HD m l k) = fromIntegral k * fromIntegral m / fromIntegral l
fromParams :: Int
-> Int
-> Int
-> HypergeometricDistribution
fromParams m l k =
assert (m > 0 && m <= l) .
assert (l > 0) .
assert (k > 0 && k <= l) $
HD m l k
probability :: HypergeometricDistribution -> Double -> Double
probability (HD mi li ki) x
| l <= 70 = (mi <> xi) * ((li mi) <> (ki xi)) / (li <> ki)
| r > maxVal = 1/0
| otherwise = exp r
where
a <> b = fromIntegral (a `choose` b)
r = f m + f (lm) f l f xi f (kxi) + f k
f (mxi) f (lmk+xi) + f (lk)
f = logFactorial
maxVal = fromIntegral (m_max_exp 1) * log 2
xi = floor x
m = fromIntegral mi
l = fromIntegral li
k = fromIntegral ki
cumulative :: HypergeometricDistribution -> Double -> Double
cumulative d@(HD m l k) x
| x < fromIntegral imin = 0
| x >= fromIntegral imax = 1
| otherwise = min r 1
where
imin = max 0 (k l + m)
imax = min k m
r = sumU . mapU (probability d . fromIntegral) . enumFromToU imin . floor $ x
inverse :: HypergeometricDistribution -> Double -> Double
inverse = error "Statistics.Distribution.Hypergeometric.inverse: not yet implemented"