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.