Archives

Archives / 2010 / August
  • 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.

    Read more...

  • F#, the ACM, and the SEC

    It all started with a twit from @mulambda: “Phil Wadler lists #fsharp as a candidate for SEC regulation spec language: http://tinyurl.com/2edfxka” I downloaded the, nonetheless, Association for Computer Machinery answer to the Securities and Exchange Commission proposal (and ask for comments) on requiring Python programs to be provided to explain contractual cash flow provisions. I quickly skimmed the ACM document and twitted “Java, C#, and F# recommended by the #ACM for SEC regulation spec language http://is.gd/e49Vt #fsharp /via @mulambda”. Later, I read with more care the ACM answer and I found that I really should clarify my twit:

    Read more...