Computing gamma, beta and beta pdf functions in postgresql

betapdf.sql

-- numberical recipes 6.1
CREATE OR REPLACE FUNCTION gammaln(xx double precision) RETURNS double precision AS $$
DECLARE
    x double precision;
    tmp double precision;
    ser double precision;
    i integer;
    cof double precision[] := ARRAY[
    57.1562356658629235,
    -59.5979603554754912, 
    14.1360979747417471,
    -0.491913816097620199,
    .339946499848118887e-4, 
    .465236289270485756e-4,
    -.983744753048795646e-4,
    .158088703224912494e-3, 
    -.210264441724104883e-3,
    .217439618115212643e-3,
    -.164318106536763890e-3, 
    .844182239838527433e-4,
    -.261908384015814087e-4,
    .368991826595316234e-5];
BEGIN
    IF xx < 0 THEN
        RAISE EXCEPTION 'negative argument';
    END IF;
    x := xx;
    tmp := x+5.24218750000000000;
    tmp := (x+0.5)*ln(tmp)-tmp;
    ser := 0.999999999999997092;
    
    FOR i IN 1..14 LOOP
        ser := ser + cof[i]/(xx + i);
    END LOOP;

    RETURN tmp+ln(2.5066282746310005*ser/x);
END;
$$ LANGUAGE plpgsql;

CREATE OR REPLACE FUNCTION _err(a double precision, b double precision) RETURNS double precision AS $$
BEGIN
    RETURN abs(a-b)/abs(b);
END;
$$ LANGUAGE plpgsql;


-- expected values were computed in Mathematica with 16 digits precision
SELECT _err(gammaln(5), 3.178053830347946) < 1e-14,
    _err(gammaln(50), 144.56574394634488600891844306296897157498) < 1e-14,
    _err(gammaln(500), 2605.11585036173389265867427356379081031402) < 1e-14,
    _err(gammaln(5000), 37582.62631568535033174656614769685808280445) < 1e-14,
    _err(gammaln(50000), 490984.42327157182173146714600978963141888836) < 1e-14,
    _err(gammaln(500000), 6061176.0464591755665203694) < 1e-14
    ;

-- numberical recipes 6.1
CREATE OR REPLACE FUNCTION beta(a double precision, b double precision) RETURNS double precision AS $$
BEGIN
    RETURN exp(gammaln(a)+gammaln(b)-gammaln(a+b));
END;
$$ LANGUAGE plpgsql;

-- expected values were computed in Mathematica with 16 digits precision
SELECT _err(beta(2,3), 0.08333333333333333) < 1e-14,
       _err(beta(5,3), 0.009523809523809524) < 1e-14,
       _err(beta(20,8),5.630440413049109e-8) < 1e-13
;


-- exp(log(wikipedia formula))
CREATE OR REPLACE FUNCTION betapdf(a double precision, b double precision, x double precision) RETURNS double precision AS $$
BEGIN
    IF a <= 0 THEN
        RAISE EXCEPTION 'Bad a';
    END IF;
    IF b <= 0 THEN
        RAISE EXCEPTION 'Bad b';
    END IF;
    IF x < 0 OR x > 1 THEN
        RAISE EXCEPTION 'Bad x';
    END IF;

    RETURN exp(gammaln(a+b)-gammaln(a)-gammaln(b)+(a-1)*ln(x)+(b-1)*ln(1-x));
END;
$$ LANGUAGE plpgsql;

-- expected values were computed in Mathematica with 16 digits precision
SELECT _err(betapdf(600,1288,.3),9.39848) < 1e-6,
       _err(betapdf(8,20,.01),1.46733e-7) < 1e-5
;

DROP FUNCTION _betacf;

-- numerical methods 6.2
CREATE OR REPLACE FUNCTION _betacf(a double precision, b double precision, x double precision) RETURNS double precision AS $$
DECLARE
    m integer;
    m2 integer;
    aa double precision;
    c double precision;
    d double precision;
    del double precision;
    h double precision;
    qab double precision;
    qam double precision;
    qap double precision;
    EPS double precision = 2.22044604925031308084726333618164e-16;
    DMIN double precision = 2.22507385850720138309023271733240e-308;
    FPMIN double precision = DMIN / EPS;
BEGIN
    qab=a+b;
    qap=a+1.0;
    qam=a-1.0;
    c=1.0;
    d=1.0-qab*x/qap;
    IF abs(d) < FPMIN THEN
        d=FPMIN;
    END IF;
    
    d=1.0/d;
    h=d;

    FOR m IN 1..10000 LOOP
        m2=2*m;
        aa=m*(b-m)*x/((qam+m2)*(a+m2)); 
        d=1.0+aa*d;
        IF abs(d) < FPMIN THEN
            d=FPMIN; 
        END IF;
        c=1.0+aa/c;
        if abs(c) < FPMIN THEN
            c=FPMIN;
        END IF;

        d=1.0/d;
        h = h*d*c;
        aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
        d=1.0+aa*d;
        if abs(d) < FPMIN THEN
            d=FPMIN; 
        END IF;
        c=1.0+aa/c;
        if abs(c) < FPMIN THEN
            c=FPMIN;
        END IF;
        d=1.0/d;
        del=d*c;
        h = h*del;
        if abs(del-1.0) <= EPS THEN
            RETURN h;
        END IF;
    END LOOP;

    RETURN h;
END;
$$ LANGUAGE plpgsql;

CREATE OR REPLACE FUNCTION _min(a double precision, b double precision) RETURNS double precision AS $$
BEGIN
    IF a < b THEN
        RETURN a;
    END IF;
    RETURN b;
END;
$$ LANGUAGE plpgsql;


CREATE OR REPLACE FUNCTION _max(a double precision, b double precision) RETURNS double precision AS $$
BEGIN
    IF a > b THEN
        RETURN a;
    END IF;
    RETURN b;
END;
$$ LANGUAGE plpgsql;


CREATE OR REPLACE FUNCTION _betaiapprox(a double precision, b double precision, x double precision) RETURNS double precision AS $$
DECLARE
    i integer;
    xu double precision;
    t double precision;
    sum double precision;
    ans double precision;
    a1 double precision = a - 1.0;
    b1 double precision = b - 1.0;
    mu double precision = a / (a + b);
    lnmu double precision = ln(mu);
    lnmuc double precision = ln(1.-mu);
    y double precision[] := ARRAY[
    0.0021695375159141994, 
    0.011413521097787704, 
    0.027972308950302116, 
    0.051727015600492421, 
    0.082502225484340941, 
    0.12007019910960293, 
    0.16415283300752470, 
    0.21442376986779355, 
    0.27051082840644336, 
    0.33199876341447887, 
    0.39843234186401943, 
    0.46931971407375483, 
    0.54413605556657973, 
    0.62232745288031077, 
    0.70331500465597174, 
    0.78649910768313447, 
    0.87126389619061517, 
    0.95698180152629142
    ];
    w double precision[] := ARRAY[
    0.0055657196642445571, 
    0.012915947284065419, 
    0.020181515297735382, 
    0.027298621498568734, 
    0.034213810770299537, 
    0.040875750923643261, 
    0.047235083490265582, 
    0.053244713977759692, 
    0.058860144245324798, 
    0.064039797355015485,
    0.068745323835736408, 
    0.072941885005653087, 
    0.076598410645870640, 
    0.079687828912071670, 
    0.082187266704339706, 
    0.084078218979661945, 
    0.085346685739338721, 
    0.085983275670394821
];
BEGIN
    t = sqrt(a * b / ((a + b)*(a + b) * (a + b + 1.0)));

    IF (x > a / (a + b)) THEN
        if (x >= 1.0) THEN
            return 1.0;
        END IF;
        xu = _min(1., _max(mu + 10. * t, x + 5.0 * t));
    ELSE
        IF (x <= 0.0) THEN
            RETURN 0.0;
        END IF;
        xu = _max(0., _min(mu - 10. * t, x - 5.0 * t));
    END IF;

    sum = 0;
    FOR j IN 1..18 LOOP
        t = x + (xu - x) * y[j];
        sum = sum + w[j] * exp(a1 * (ln(t) - lnmu) + b1 * (ln(1 - t) - lnmuc));
    END LOOP;

    ans = sum * (xu - x) * exp(a1 * lnmu - gammaln(a) + b1 * lnmuc - gammaln(b) + gammaln(a + b));

    IF ans > 0.0 THEN 
        return 1.0 - ans;
    ELSE   
        return -ans;
    END IF;
END;
$$ LANGUAGE plpgsql;

-- numerical methods 6.2
CREATE OR REPLACE FUNCTION betacdf(a double precision, b double precision, x double precision) RETURNS double precision AS $$
DECLARE
    switch double precision = 3000;
    bt double precision;
BEGIN
    IF a <= 0 THEN
        RAISE EXCEPTION 'Bad a';
    END IF;
    IF b <= 0 THEN
        RAISE EXCEPTION 'Bad b';
    END IF;
    IF x < 0 OR x > 1 THEN
        RAISE EXCEPTION 'Bad x';
    END IF;

    IF x = 0.0 OR x = 1.0 THEN
        RETURN x;
    END IF;

    IF a > SWITCH AND b > SWITCH THEN
        RETURN _betaiapprox(a,b,x); 
    END IF;
    
    bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b)+a*ln(x)+b*ln(1.0-x)); 

    IF x < (a+1.0)/(a+b+2.0) THEN
        RETURN bt * _betacf(a,b,x)/a;
    ELSE 
        RETURN 1.0-bt * _betacf(b,a,1.0-x)/b;
    END IF;
END;
$$ LANGUAGE plpgsql;

-- expected values were computed in Mathematica with 16 digits precision
SELECT _err(betacdf(600,1288,.3), 0.0472618) < 1e-6,
       _err(betacdf(8,20,.5), 0.990421) < 1e-6,
     _err(betacdf(6000,12880,.31), 0.0104679) < 1e-5;