1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
{-# 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 (0 ,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.
|
