picomath

Haskell

NormalCdfInverse.hs

module NormalCdfInverse (normalCdfInverse) where

rational_approximation t = let
    -- Abramowitz and Stegun formula 26.2.23.
    -- The absolute value of the error should be less than 4.5 e-4.
    c = [2.515517, 0.802853, 0.010328]
    d = [1.432788, 0.189269, 0.001308]
    numerator = ((c!!2)*t + (c!!1))*t + (c!!0)
    denominator = (((d!!2)*t + (d!!1))*t + (d!!0))*t + 1.0
    in t - numerator / denominator

normalCdfInverse p | p > 0.0 && p < 1 =
    -- See article above for explanation of this section.
    if p < 0.5
        -- F^-1(p) = - G^-1(p)
        then -rational_approximation( sqrt(-2.0*log(p)) )
        -- F^-1(p) = G^-1(1-p)
        else rational_approximation( sqrt(-2.0*log(1.0-p)) )