Implementation of NORMSINV function in MySQL

I need to use the inverse normal cumulative distribution function (what in Excel is NORMSINV function) in MySQL, but there is no such function there. Maybe any of you has an implementation of it for MySQL?

Many thanks for your tim开发者_高级运维e.

Well, I've finally found this. It's not perfect but a quite good aproximation. The code isn't mine, its author is Geoffrey C. Barnes. I've just transformed it from VB.NET to MySQL.

DROP FUNCTION IF EXISTS NORMSINV;
DELIMITER //

CREATE FUNCTION NORMSINV (p DOUBLE) RETURNS DOUBLE
BEGIN

DECLARE q, r DOUBLE;
DECLARE A1, A2, A3, A4, A5, A6 DOUBLE;
DECLARE B1, B2, B3, B4, B5 DOUBLE;
DECLARE C1, C2, C3, C4, C5, C6 DOUBLE;
DECLARE D1, D2, D3, D4 DOUBLE;
DECLARE P_LOW, P_HIGH DOUBLE;

/* coefficients in rational approximations */
SET A1 = -39.696830286653757;
SET A2 = 220.9460984245205;
SET A3 = -275.92851044696869;
SET A4 = 138.357751867269;
SET A5 = -30.66479806614716;
SET A6 = 2.5066282774592392;

SET B1 = -54.476098798224058;
SET B2 = 161.58583685804089;
SET B3 = -155.69897985988661;
SET B4 = 66.80131188771972;
SET B5 = -13.280681552885721;

SET C1 = -0.0077848940024302926;
SET C2 = -0.32239645804113648;
SET C3 = -2.4007582771618381;
SET C4 = -2.5497325393437338;
SET C5 = 4.3746641414649678;
SET C6 = 2.9381639826987831;

SET D1 = 0.0077846957090414622;
SET D2 = 0.32246712907003983;
SET D3 = 2.445134137142996;
SET D4 = 3.7544086619074162;

/* define break points */
SET P_LOW = 0.02425;
SET P_HIGH = 1 - P_LOW;

IF (p > 0 AND p < P_LOW) THEN
/* rational approximation for lower region */
SET q = SQRT(-2 * LOG(p));
RETURN (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) /  
       ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);

ELSEIF (p >= P_LOW AND p <= P_HIGH) THEN
/* rational approximation for central region */
SET q = p - 0.5;
SET 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);

ELSEIF (p > P_HIGH AND p < 1) THEN
/* rational approximation for upper region */
SET q = SQRT(-2 * LOG(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//
DELIMITER ;

Below is an equivalence of Excel's NORMINV function.

CREATE FUNCTION NORMINV(p double, mu double, sigma double)
RETURNS decimal(20,6)
begin
    declare q   double;
    declare r   double;
    declare val double;
BEGIN
    IF p < 0 OR p > 1 THEN
        signal sqlstate '20100' set message_text= 'The probality p must be bigger than 0 and smaller than 1';
    END IF;
    IF sigma < 0 THEN
        signal sqlstate '20100' set message_text=  'The standard deviation sigma must be positive';
    END IF;
    IF p = 0 THEN 
        RETURN to_binary_double('-INF');    
    END IF;
    IF p = 1 THEN
        RETURN to_binary_double('INF');
    END IF;
    IF sigma = 0 THEN
        RETURN mu;
    END IF;

    set q:= p - 0.5;

    IF(ABS(q) <= .425) THEN
            set r:= .180625 - q * q;
            set val :=
                 q * (((((((r * 2509.0809287301226727 +
                            33430.575583588128105) * r + 67265.770927008700853) * r +
                          45921.953931549871457) * r + 13731.693765509461125) * r +
                        1971.5909503065514427) * r + 133.14166789178437745) * r +
                      3.387132872796366608)
                 / (((((((r * 5226.495278852854561 +
                          28729.085735721942674) * r + 39307.89580009271061) * r +
                        21213.794301586595867) * r + 5394.1960214247511077) * r +
                      687.1870074920579083) * r + 42.313330701600911252) * r + 1);
    ELSE
        /* r = min(p, 1-p) < 0.075 */
        IF q > 0 THEN
            set r:= 1 - p;
        ELSE
            set r:= p;
        END IF;

        set r:= SQRT(-LN(r));
        /* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */

        IF (r <= 5) THEN
            set r:= r - 1.6;
            set val := (((((((r * 7.7454501427834140764e-4 +
                         .0227238449892691845833) * r + .24178072517745061177) *
                       r + 1.27045825245236838258) * r +
                      3.64784832476320460504) * r + 5.7694972214606914055) *
                    r + 4.6303378461565452959) * r +
                   1.42343711074968357734)
                  / (((((((r *
                           1.05075007164441684324e-9 + 5.475938084995344946e-4) *
                          r + .0151986665636164571966) * r +
                         .14810397642748007459) * r + .68976733498510000455) *
                       r + 1.6763848301838038494) * r +
                      2.05319162663775882187) * r + 1);
        ELSE /* very close to  0 or 1 */
            set r := r - 5;
            set val := (((((((r * 2.01033439929228813265e-7 +
                         2.71155556874348757815e-5) * r +
                        .0012426609473880784386) * r + .026532189526576123093) *
                      r + .29656057182850489123) * r +
                     1.7848265399172913358) * r + 5.4637849111641143699) *
                   r + 6.6579046435011037772)
                  / (((((((r *
                           2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
                          r + 1.8463183175100546818e-5) * r +
                         7.868691311456132591e-4) * r + .0148753612908506148525)
                       * r + .13692988092273580531) * r +
                      .59983220655588793769) * r + 1);
        END IF;

        IF q < 0.0 THEN
            set val := -val;
        END IF;
    END IF;

    RETURN mu + sigma * val;
END;
END

It has been adapted from [https://forums.mysql.com/read.php?98,411441,411458#msg-411458] which is Oracle equivalence. Have used it, so far so good. Hope it assists someone.