
  0 | module Statistics.Norm.Precise
  1 | 
  2 | import Statistics.Erf
  3 | import public Statistics.Norm
  4 | import Statistics.Norm.Rough
  5 | 
  6 | %default total
  7 | 
  8 | export
  9 | NormCDF where
 10 |   normcdf x = P $ erfc (-x / sqrt 2) / 2
 11 | 
 12 | -- theoretically, it can get rough `InvNormCDF` as an auto-implicit, i.e. to improve any implementation
 13 | export
 14 | InvNormCDF where
 15 |   -- Code below is based on taken from
 16 |   -- https://hackage.haskell.org/package/erf-2.0.0.0/docs/src/Data-Number-Erf.html
 17 |   -- Licensed with BSD-3-Clause, copyright Lennart Augustsson
 18 |   invnormcdf p =
 19 |     if      p == 0 then NegInf
 20 |     else if p == 1 then PosInf
 21 |     else let -- Do one iteration with Halley's root finder to get a more accurate result.
 22 |       x = invnormcdf @{Rough} p
 23 |       e = (normcdf x).asDouble - p.asDouble
 24 |       x = x.asDouble
 25 |       u = e * sqrt (2*pi) * exp (x*x / 2)
 26 |       x' = x - u / (1 + x * u / 2)
 27 |       in roughlyFit x'
 28 |