NORMSINV function in PL/SQL

 NORMSINV Formula in Excel, calculates the inverse of the normal cumulative distribution for a given probability. Below is the PL/SQL function for the same :

FUNCTION NORMSINV(p BINARY_DOUBLE)
RETURN BINARY_DOUBLE
IS
  q      BINARY_DOUBLE;
  r      BINARY_DOUBLE;
  a1     BINARY_DOUBLE;
  a2     BINARY_DOUBLE;
  a3     BINARY_DOUBLE;
  a4     BINARY_DOUBLE;
  a5     BINARY_DOUBLE;
  a6     BINARY_DOUBLE;
  b1     BINARY_DOUBLE;
  b2     BINARY_DOUBLE;
  b3     BINARY_DOUBLE;
  b4     BINARY_DOUBLE;
  b5     BINARY_DOUBLE;
  c1     BINARY_DOUBLE;
  c2     BINARY_DOUBLE;
  c3     BINARY_DOUBLE;
  c4     BINARY_DOUBLE;
  c5     BINARY_DOUBLE;
  c6     BINARY_DOUBLE;
  d1     BINARY_DOUBLE;
  d2     BINARY_DOUBLE;
  d3     BINARY_DOUBLE;
  d4     BINARY_DOUBLE;
  p_low  BINARY_DOUBLE;
  p_high BINARY_DOUBLE;
BEGIN
    /* coefficients in rational approximations */
    a1 := -39.696830286653757;
    a2 := 220.9460984245205;
    a3 := -275.92851044696869;
    a4 := 138.357751867269;
    a5 := -30.66479806614716;
    a6 := 2.5066282774592392;
    b1 := -54.476098798224058;
    b2 := 161.58583685804089;
    b3 := -155.69897985988661;
    b4 := 66.80131188771972;
    b5 := -13.280681552885721;
    c1 := -0.0077848940024302926;
    c2 := -0.32239645804113648;
    c3 := -2.4007582771618381;
    c4 := -2.5497325393437338;
    c5 := 4.3746641414649678;
    c6 := 2.9381639826987831;
    d1 := 0.0077846957090414622;
    d2 := 0.32246712907003983;
    d3 := 2.445134137142996;
    d4 := 3.7544086619074162;

    /* define break points */
    p_low := 0.02425;
    p_high := 1 - p_low;

    IF ( p > 0
         AND p < p_low ) THEN
      /* rational approximation for lower region */
      q := Sqrt (-2 * Ln (p));

      RETURN ( ( ( ( ( c1 * q + c2 ) * q + c3 ) * q + c4 ) * q + c5 ) * q + c6 )
             /
             ( ( ( ( d1 * q + d2 ) * q + d3 ) * q + d4 ) * q + 1 );
    ELSIF ( p >= p_low
            AND p <= p_high ) THEN
      /* rational approximation for central region */
      q := p - 0.5;
      r := q * q;

      RETURN ( ( ( ( ( a1 * r + a2 ) * r + a3 ) * r + a4 ) * r + a5 ) * r + a6 )
             * q
             / ( ( (
             ( ( b1 * r + b2 ) * r + b3 ) * r + b4 ) * r + b5 ) * r + 1 );
    ELSIF ( p > p_high
            AND p < 1 ) THEN
      /* rational approximation for upper region */
      q := Sqrt (-2 * Ln (1 - p));

      RETURN -( ( ( ( ( c1 * q + c2 ) * q + c3 ) * q + c4 ) * q + c5 ) * q + c6
              ) /
             ( ( ( (
             d1 * q + d2 ) * q + d3 ) * q + d4 ) * q + 1 );
    /* on error returning 0 */
    ELSE
      RETURN 0;
    END IF;
END normsinv;