
  0 | module Statistics.Norm.Rough
  1 | 
  2 | import public Statistics.Norm
  3 | 
  4 | %default total
  5 | 
  6 | ||| Calculation of the inverse of CDF for normal distribution with precision only up to 10^9.
  7 | ||| This implementation does not have any external dependencies.
  8 | export %defaulthint
  9 | Rough : InvNormCDF
 10 | Rough = R where
 11 |   -- Code below is based on taken from
 12 |   -- https://hackage.haskell.org/package/erf-2.0.0.0/docs/src/Data-Number-Erf.html
 13 |   -- Licensed with BSD-3-Clause, copyright Lennart Augustsson
 14 |   [R] InvNormCDF where
 15 |     invnormcdf p =
 16 |         let res = if      p == 0       then NegInf
 17 |                   else if p == 1       then PosInf
 18 |                   else if p < pLow     then closeToLowBound p
 19 |                   else if p.inv < pLow then - closeToLowBound p.inv
 20 |                   else                      middling p
 21 |         in roughlyFit res
 22 | 
 23 |       where
 24 |         pLow : Probability
 25 |         pLow = 0.02425
 26 | 
 27 |         closeToLowBound : Probability -> Double
 28 |         closeToLowBound p = let
 29 |           c1 = -7.784894002430293e-03
 30 |           c2 = -3.223964580411365e-01
 31 |           c3 = -2.400758277161838e+00
 32 |           c4 = -2.549732539343734e+00
 33 |           c5 =  4.374664141464968e+00
 34 |           c6 =  2.938163982698783e+00
 35 | 
 36 |           d1 =  7.784695709041462e-03
 37 |           d2 =  3.224671290700398e-01
 38 |           d3 =  2.445134137142996e+00
 39 |           d4 =  3.754408661907416e+00
 40 | 
 41 |           q = sqrt(-2*log(p.asDouble))
 42 | 
 43 |           in (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) /
 44 |              ((((d1*q+d2)*q+d3)*q+d4)*q+1)
 45 | 
 46 |         %inline
 47 |         middling : Probability -> Double
 48 |         middling p = let
 49 |           a1 = -3.969683028665376e+01
 50 |           a2 =  2.209460984245205e+02
 51 |           a3 = -2.759285104469687e+02
 52 |           a4 =  1.383577518672690e+02
 53 |           a5 = -3.066479806614716e+01
 54 |           a6 =  2.506628277459239e+00
 55 | 
 56 |           b1 = -5.447609879822406e+01
 57 |           b2 =  1.615858368580409e+02
 58 |           b3 = -1.556989798598866e+02
 59 |           b4 =  6.680131188771972e+01
 60 |           b5 = -1.328068155288572e+01
 61 | 
 62 |           q = p.asDouble - 0.5
 63 |           r = q*q
 64 | 
 65 |           in (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q /
 66 |              (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1)
 67 |