## 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.