0 | module Statistics.Confidence
  1 |
  2 | import public Data.Functor.TraverseSt
  3 |
  4 | import Data.DPair
  5 | import public Data.Nat
  6 | import Data.String
  7 | import public Data.Vect
  8 |
  9 | import public Statistics.Norm.Rough
 10 | import public Statistics.Probability
 11 |
 12 | -- This module reexports the rough implementation for `InvNormCDF` for compatility reasons, this may be changed in the future
 13 |
 14 | %default total
 15 |
 16 | -- In order to get an accurate measurement with small sample sizes, we're using the Wilson score interval
 17 | -- https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Wilson_score_interval
 18 | export
 19 | wilsonBounds :
 20 |   InvNormCDF =>
 21 |   (confidence : Probability) ->
 22 |   (count : Nat) ->
 23 |   (0 _ : IsSucc count) =>
 24 |   (successes : Probability) ->
 25 |   (Probability, Probability)
 26 | wilsonBounds confidence count successes =
 27 |   let
 28 |     n = cast count
 29 |     p = successes.asDouble
 30 |     z = cast $ invnormcdf $ inv $ confidence / 2
 31 |     zz_n = z * z / n
 32 |
 33 |     midpoint = p + zz_n / 2
 34 |
 35 |     offset = z / (1 + zz_n) * sqrt (p * (1 - p) / n + zz_n / (4 * n))
 36 |
 37 |     denominator = 1 + zz_n
 38 |
 39 |     low  = (midpoint - offset) / denominator
 40 |     high = (midpoint + offset) / denominator
 41 |
 42 |   in if low == low && high == high
 43 |        then mapHom (P . roughlyFit) (low, high)
 44 |        else (0, 1) -- we've gone too close to infinite `z`
 45 |
 46 | --- Stuff for expressing what and how should be tested ---
 47 |
 48 | public export
 49 | record CoverageTest a where
 50 |   constructor Cover
 51 |   checkedProbability : (Probability, Probability)
 52 |   {auto 0 minMaxCorrect : So $ fst checkedProbability <= snd checkedProbability}
 53 |   successCondition : a -> Bool
 54 |
 55 | namespace CoverageTest
 56 |
 57 |   public export %inline
 58 |   minimalProbability : CoverageTest a -> Probability
 59 |   minimalProbability = fst . checkedProbability
 60 |
 61 |   public export %inline
 62 |   (.minimalProbability) : CoverageTest a -> Probability
 63 |   (.minimalProbability) = minimalProbability
 64 |
 65 |   public export %inline
 66 |   maximalProbability : CoverageTest a -> Probability
 67 |   maximalProbability = snd . checkedProbability
 68 |
 69 |   public export %inline
 70 |   (.maximalProbability) : CoverageTest a -> Probability
 71 |   (.maximalProbability) = maximalProbability
 72 |
 73 | namespace CoverageTest
 74 |
 75 |   public export %inline
 76 |   DefaultTolerance : Probability
 77 |   DefaultTolerance = 95.percent
 78 |
 79 |   export
 80 |   coverBetween : {default DefaultTolerance tolerance : Probability} ->
 81 |                  (minP, maxP : Probability) ->
 82 |                  (successCondition : a -> Bool) ->
 83 |                  CoverageTest a
 84 |   coverBetween minP maxP = do
 85 |     let (minP, maxP) = (min minP maxP, max minP maxP)
 86 |     let minP = minP * tolerance
 87 |     let maxP = inv $ maxP.inv * tolerance
 88 |     let maxP = if maxP <= minP then minP else maxP -- just in case `maxP` is a non-normalised very small value
 89 |     Cover (minP, maxP) @{believe_me Oh}
 90 |
 91 |   export
 92 |   coverMin : {default DefaultTolerance tolerance : Probability} ->
 93 |              (minP : Probability) ->
 94 |              (successCondition : a -> Bool) ->
 95 |              CoverageTest a
 96 |   coverMin minP = coverBetween {tolerance} minP 1
 97 |
 98 |   export
 99 |   coverWith : {default DefaultTolerance tolerance : Probability} ->
100 |               (singleP : Probability) ->
101 |               (successCondition : a -> Bool) ->
102 |               CoverageTest a
103 |   coverWith singleP = coverBetween {tolerance} singleP singleP
104 |
105 | public export
106 | DefaultConfidence : Probability
107 | DefaultConfidence = 1/1000000000
108 |
109 | --- Coverage condition checker with low-level intermediate results ---
110 |
111 | public export
112 | record CoverageTestState where
113 |   constructor MkCoverageTestState
114 |   testedBounds : (Probability, Probability)
115 |   wilsonBounds : (Probability, Probability)
116 |
117 | export
118 | [NoAnalysis] Interpolation CoverageTestState where
119 |   interpolate $ MkCoverageTestState (minP, maxP) (wLow, wHigh) = do
120 |     let ss = if minP == maxP && wLow /= minP && wHigh /= maxP
121 |                then [ (wLow , "[\{show wLow}")
122 |                     , (minP , "(\{show minP})")
123 |                     , (wHigh, "\{show wHigh}]")
124 |                     ]
125 |                else [ (minP , "(\{show minP}")
126 |                     , (wLow , "[\{show wLow}")
127 |                     , (wHigh, "\{show wHigh}]")
128 |                     , (maxP , "\{show maxP})")
129 |                     ]
130 |     unwords $ map snd $ sortBy (comparing fst) ss
131 |
132 | export
133 | [Means] Interpolation CoverageTestState where
134 |   interpolate $ MkCoverageTestState expected wilson = "expected: \{interpolatePair expected}, wilson: \{interpolatePair wilson}" where
135 |     interpolatePair : (Probability, Probability) -> String
136 |     interpolatePair (lo, hi) = do
137 |       let lo = lo.value
138 |           hi = hi.value
139 |       if lo == hi then "\{show lo}" else with Bounded.(+) with Bounded.(/) with Bounded.(-) do
140 |         let mi = (lo + hi) / 2
141 |         let di = hi - mi
142 |         "\{show mi} ± \{show di}"
143 |
144 | export
145 | checkCoverageConditions' :
146 |   TraversableSt t =>
147 |   InvNormCDF =>
148 |   {default DefaultConfidence confidence : Probability} ->
149 |   Vect n (CoverageTest a) ->
150 |   t a ->
151 |   t $ Vect n CoverageTestState
152 |
153 | checkCoverageConditions' coverageTests = traverseSt initialResults checkCoverageOnce where
154 |
155 |   data PastResults : Type where
156 |     R : (attempts : Nat) -> (successes : Vect n $ Subset Nat (`LTE` attempts)) -> PastResults
157 |
158 |   initialResults : PastResults
159 |   initialResults = R 0 $ coverageTests <&> const (0 `Element` LTEZero)
160 |
161 |   checkCoverageOnce : PastResults -> a -> (PastResults, Vect n CoverageTestState)
162 |   checkCoverageOnce (R prevAttempts prevResults) x = do
163 |     let %inline currAttempts : NatcurrAttempts = S prevAttempts
164 |     mapFst (R currAttempts) $ unzip $ coverageTests `zip` prevResults <&>
165 |       \(Cover checkedP@(minP, maxP) cond, Element prevSucc _) => do
166 |         let pr@(Element currSucc _) = if cond x
167 |                                         then S prevSucc `Element` LTESucc %search
168 |                                         else prevSucc   `Element` lteSuccRight %search
169 |         (pr, MkCoverageTestState checkedP $ wilsonBounds confidence currAttempts $ ratio currSucc currAttempts)
170 |
171 | --- Coverage condition checkers with simpler results ---
172 |
173 | public export
174 | data SignificantBounds
175 |   = BoundsOk
176 |   | UpperBoundViolated Probability
177 |   | LowerBoundViolated Probability
178 |
179 | public export
180 | isOk : SignificantBounds -> Bool
181 | isOk BoundsOk = True
182 | isOk _        = False
183 |
184 | export
185 | Interpolation SignificantBounds where
186 |   interpolate BoundsOk               = "ok"
187 |   interpolate $ UpperBoundViolated x = "(..., ...) --> \{show x}"
188 |   interpolate $ LowerBoundViolated x = "\{show x} <-- (..., ...)"
189 |
190 | -- `Just` means "statistical significance", `Nothing` means "no significance yet".
191 | public export %inline
192 | CoverageTestResult : Type
193 | CoverageTestResult = Maybe SignificantBounds
194 |
195 | export
196 | coverageTestResult : CoverageTestState -> CoverageTestResult
197 | coverageTestResult $ MkCoverageTestState (minP, maxP) (wLow, wHigh) =
198 |   if      wLow >= minP && wHigh <= maxP then Just BoundsOk
199 |   else if wLow > maxP                   then Just $ UpperBoundViolated wLow
200 |   else if                 wHigh < minP  then Just $ LowerBoundViolated wHigh
201 |   else                                       Nothing
202 |
203 | export
204 | Interpolation CoverageTestState where
205 |   interpolate cts@(MkCoverageTestState (minP, maxP) (wLow, wHigh)) = case coverageTestResult cts of
206 |     Just BoundsOk               => "ok"
207 |     Just (UpperBoundViolated _) => "(..., \{show maxP}) --> \{show wLow}"
208 |     Just (LowerBoundViolated _) => "\{show wLow} <-- (\{show minP}, ...)"
209 |     Nothing                     => "\{show wLow} .. \{show wHigh} still covers (\{show minP}, \{show maxP})"
210 |
211 | export %inline
212 | checkCoverageConditions :
213 |   TraversableSt t =>
214 |   InvNormCDF =>
215 |   {default DefaultConfidence confidence : Probability} ->
216 |   Vect n (CoverageTest a) ->
217 |   t a ->
218 |   t $ Vect n CoverageTestResult
219 | checkCoverageConditions = map @{Compose} coverageTestResult .: checkCoverageConditions' {confidence}
220 |
221 | export %inline
222 | checkCoverageCondition :
223 |   TraversableSt t =>
224 |   InvNormCDF =>
225 |   {default DefaultConfidence confidence : Probability} ->
226 |   CoverageTest a ->
227 |   t a ->
228 |   t CoverageTestResult
229 | checkCoverageCondition ct = map head . checkCoverageConditions {confidence} [ct]
230 |