Skip to content

NormsInv function in F#

July 20, 2011

NORMSINV(p) returns the value z such that, with probability p, a standard normal random variable takes on a value that is less than or equal to z. A standard normal random variable has mean 0 and standard deviation 1 (and also variance 1 because variance = standard deviation squared)

// An algorithm for computing the inverse normal cumulative distribution function
let normsInv(p : float) = 
    //todo: validate that p is between 0 and 1

    //Coefficients in rational approximations      
    let a = [|-3.969683028665376e+01; 2.209460984245205e+02;
             -2.759285104469687e+02; 1.383577518672690e+02;
             -3.066479806614716e+01; 2.506628277459239e+00|]

    let b = [|-5.447609879822406e+01; 1.615858368580409e+02;
             -1.556989798598866e+02; 6.680131188771972e+01;
             -1.328068155288572e+01|]

    let c = [|-7.784894002430293e-03; -3.223964580411365e-01;
             -2.400758277161838e+00; -2.549732539343734e+00;
             4.374664141464968e+00; 2.938163982698783e+00|]

    let d = [|7.784695709041462e-03; 3.224671290700398e-01;
             2.445134137142996e+00; 3.754408661907416e+00|]

    //Define break-points.
    let plow = 0.02425
    let phigh = 1. - plow

    let approxLower(var) : float = 
        let q = Math.Sqrt(-2. * Math.Log(var))
        let r = (((((c.[0] * q + c.[1]) * q + c.[2]) * q + c.[3]) * q + c.[4]) * q + c.[5]) / ((((d.[0] * q + d.[1]) * q + d.[2]) * q + d.[3]) * q + 1.)
        r

    let approxUpper(var) = 
        let q = Math.Sqrt(-2. * Math.Log(1. - var));
        let r = (((((c.[0] * q + c.[1]) * q + c.[2]) * q + c.[3]) * q + c.[4]) * q + c.[5]) /  ((((d.[0] * q + d.[1]) * q + d.[2]) * q + d.[3]) * q + 1.)
        r * -1.

    let approxCentral(var) = 
        let q = var - 0.5;
        let r = q * q
        (((((a.[0] * r + a.[1]) * r + a.[2]) * r + a.[3]) * r + a.[4]) * r + a.[5]) * q / (((((b.[0] * r + b.[1]) * r + b.[2]) * r + b.[3]) * r + b.[4]) * r + 1.)

    match p with 
    | var when p < plow  -> approxLower(var)
    | var when p > phigh -> approxUpper(var)
    | _ -> approxCentral(p)
Advertisements
No comments yet

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: