{-# LANGUAGE ViewPatterns #-}
-- | Calculation of confidence intervals
module Statistics.ConfidenceInt (
    poissonCI
  , poissonNormalCI
  , binomialCI
  , naiveBinomialCI
    -- * References
    -- $references
  ) where

import Statistics.Distribution
import Statistics.Distribution.ChiSquared
import Statistics.Distribution.Beta
import Statistics.Types



-- | Calculate confidence intervals for Poisson-distributed value
-- using normal approximation
poissonNormalCI :: Int -> Estimate NormalErr Double
poissonNormalCI n
  | n < 0     = error "Statistics.ConfidenceInt.poissonNormalCI negative number of trials"
  | otherwise = estimateNormErr n' (sqrt n')
  where
    n' = fromIntegral n

-- | Calculate confidence intervals for Poisson-distributed value for
--   single measurement. These are exact confidence intervals
poissonCI :: CL Double -> Int -> Estimate ConfInt Double
poissonCI cl@(significanceLevel -> p) n
  | n <  0    = error "Statistics.ConfidenceInt.poissonCI: negative number of trials"
  | n == 0    = estimateFromInterval m (m1,m2) cl
  | otherwise = estimateFromInterval m (m1,m2) cl
  where
    m  = fromIntegral n
    m1 = 0.5 * quantile      (chiSquared (2*n  )) (p/2)
    m2 = 0.5 * complQuantile (chiSquared (2*n+2)) (p/2)

-- | Calculate confidence interval using normal approximation. Note
--   that this approximation breaks down when /p/ is either close to 0
--   or to 1. In particular if @np < 5@ or @1 - np < 5@ this
--   approximation shouldn't be used.
naiveBinomialCI :: Int         -- ^ Number of trials
                -> Int         -- ^ Number of successes
                -> Estimate NormalErr Double
naiveBinomialCI n k
  | n <= 0 || k < 0 = error "Statistics.ConfidenceInt.naiveBinomialCI: negative number of events"
  | k > n           = error "Statistics.ConfidenceInt.naiveBinomialCI: more successes than trials"
  | otherwise       = estimateNormErr p σ
  where
    p = fromIntegral k / fromIntegral n
    σ = sqrt $ p * (1 - p) / fromIntegral n


-- | Clopper-Pearson confidence interval also known as exact
--   confidence intervals.
binomialCI :: CL Double
           -> Int               -- ^ Number of trials
           -> Int               -- ^ Number of successes
           -> Estimate ConfInt Double
binomialCI cl@(significanceLevel -> p) ni ki
  | ni <= 0 || ki < 0 = error "Statistics.ConfidenceInt.binomialCI: negative number of events"
  | ki > ni           = error "Statistics.ConfidenceInt.binomialCI: more successes than trials"
  | ki == 0           = estimateFromInterval eff (0, ub) cl
  | ni == ki          = estimateFromInterval eff (lb,0 ) cl
  | otherwise         = estimateFromInterval eff (lb,ub) cl
  where
    k   = fromIntegral ki
    n   = fromIntegral ni
    eff = k / n
    lb  = quantile      (betaDistr  k      (n - k + 1)) (p/2)
    ub  = complQuantile (betaDistr (k + 1) (n - k)    ) (p/2)


-- $references
--
--  * Clopper, C.; Pearson, E. S. (1934). "The use of confidence or
--    fiducial limits illustrated in the case of the
--    binomial". Biometrika 26: 404–413. doi:10.1093/biomet/26.4.404
--
--  * Brown, Lawrence D.; Cai, T. Tony; DasGupta, Anirban
--    (2001). "Interval Estimation for a Binomial Proportion". Statistical
--    Science 16 (2): 101–133. doi:10.1214/ss/1009213286. MR 1861069.
--    Zbl 02068924.