A quick and dirty implementation of Excel NORMINV in F#

A couple of weeks ago I posted an example implementation of Excel NORMINV function in C#, in there I mentioned that what I actually needed was an F# version, but I used C# as an stepping stone moving from the C++ version I originally found. Well, here you have my attempt at implementing NORMINV in F#:

let probit mu sigma p =
    if p < 0. || p > 1. then failwith "The probability p must be greater than 0 and lower than 1"
    elif sigma < 0. then failwith "The standard deviation sigma must be positive"
    elif sigma = 0. then mu
    else
        let q = p - 0.5
        let value =
            if abs q <= 0.425 then          // 0.075 <= p <= 0.925
                let r = 0.180625 - q * q
                let a0, a1, a2, a3, a4, a5, a6, a7 =
                    3.387132872796366608, 133.14166789178437745, 1971.5909503065514427, 13731.693765509461125,
                    45921.953931549871457, 67265.770927008700853, 33430.575583588128105, 2509.0809287301226727
                let b0, b1, b2, b3, b4, b5, b6, b7 =
                    1., 42.313330701600911252, 687.1870074920579083, 5394.1960214247511077,
                    21213.794301586595867, 39307.89580009271061, 28729.085735721942674, 5226.495278852854561
                q*(((((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2)*r + a1)*r + a0)
                    / (((((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2)*r + b1)*r + b0)
            else                            // closer than 0.075 from {0,1} boundary
                let r = sqrt -(log (if q > 0. then 1. - p else p))
                let val1 =
                    if r <= 5. then         // <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11
                let a0, a1, a2, a3, a4, a5, a6, a7 =
                     1.42343711074968357734, 4.6303378461565452959, 5.7694972214606914055, 3.64784832476320460504,
                     1.27045825245236838258, 0.24178072517745061177, 0.0227238449892691845833, 7.7454501427834140764e-4
                 let b0, b1, b2, b3, b4, b5, b6, b7 =
                     1., 2.05319162663775882187, 1.6763848301838038494, 0.68976733498510000455,
                    0.14810397642748007459, 0.0151986665636164571966, 5.475938084995344946e-4, 1.05075007164441684324e-9
                        let r1 = r - 1.6
                        (((((((a7*r1 + a6)*r1 + a5)*r1 + a4)*r1 + a3)*r1 + a2)*r1 + a1)*r1 + a0)
                            / (((((((b7*r1 + b6)*r1 + b5)*r1 + b4)*r1 + b3)*r1 + b2)*r1 + b1)*r1 + b0)
                    else                    // very close to  0 or 1
                        let a0, a1, a2, a3, a4, a5, a6, a7 =
                            6.6579046435011037772, 5.4637849111641143699, 1.7848265399172913358, 0.29656057182850489123,
                            0.026532189526576123093, 0.0012426609473880784386, 2.71155556874348757815e-5, 2.01033439929228813265e-7
                        let b0, b1, b2, b3, b4, b5, b6, b7 =
                            1., 0.59983220655588793769, 0.13692988092273580531, 0.0148753612908506148525,
                            7.868691311456132591e-4, 1.8463183175100546818e-5, 1.4215117583164458887e-7, 2.04426310338993978564e-15
                        let r1 = r - 5.
                        (((((((a7*r1 + a6)*r1 + a5)*r1 + a4)*r1 + a3)*r1 + a2)*r1 + a1)*r1 + a0)
                            / (((((((b7*r1 + b6)*r1 + b5)*r1 + b4)*r1 + b3)*r1 + b2)*r1 + b1)*r1 + b0)
                if q < 0. then -val1 else val1
        sigma*value + mu

I wouldn’t get angry if some of you think the code is ugly, I’m afraid such is the way of evaluating several polynomials (each pair of polynomials carefully tuned for a specific segment) and then dividing them. At any rate, I invite you to compare this version with the C# version (let alone the C++ version), I hope you find the F# implementation terser. You can download a Visual Studio 2010 sample project here:

Please remember I have made only a handful of tests, so I encourage you to do (a lot of) additional testing if you plan to use this function for anything even slightly serious.

Published Sunday, August 8, 2010 11:24 PM by Edgar Sánchez

Comments

No Comments