HP PPL and Statistical distributions

[Albert Chan]

Hi, robmio

Here is the translated quad, in PPL.
I hope this is right, as I am learning PPL as I code.

Compare to Lua or Python, there are some gothas:

1. all local variables have to declare up front.
2. there is no Lua elseif or Python elif equivalent.
3. there is no logic precedence: a OR b AND c translated to ((a OR b) AND c)
4. comparison operators has "sugar" in it, which may translated wrongly.

Cas> a > b > c                         → a > b AND b > c ... looks good
Cas> a > b > c > d                   → ((a>b) AND (b>c)) > d ???
Cas> (a>b) == (mode == 2)     → a>b and b == mode == 2 ???

Code:

#cas
quad6(f, a, b, d, n, eps)
BEGIN

LOCAL k, c, t0, t1, t, p, s, err;
LOCAL fp, fm, u, r, x, y, mode;

fp := fm := k := c := 0.;
t0 := 7.38905609893065; // e^2
a := approx(a);
b := approx(b);
mode := (a-a!=0) + 2*(b-b!=0);

IF mode == 0 THEN       // tanh-sinh
    c,d,s := d, (b-a)/2, f((a+b)/2)
ELSE
    IF mode == 3 THEN   // sinh-sinh
        s := f(c)
    ELSE                // exp-sinh
        c, s := a, b;   // c = finite edge
        IF mode == 1 THEN c, s := b, a END
        IF s < 0 THEN d:=-d END
        s := f(c+d)
    END
END

REPEAT
    p, t, t1 := 0, sqrt(t0), t0;
    IF k==0 THEN t1 := t END
    t0 := t;
    IF mode == 0 THEN           // tanh-sinh
        REPEAT
            u := exp(1/t-t);
            r := 2*u/(1+u);
            x := d*r;
            IF a+x != a THEN
                y := f(a+x);
                IF y-y==0 THEN fp:=y ELSE fp*=c END
            END
            IF b-x != b THEN
                y := f(b-x);
                IF y-y==0 THEN fm:=y ELSE fm*=c END
            END
            y := (t+1/t)*r/(1+u) * (fp+fm);
            p, t := p+y, t*t1;
        UNTIL abs(y) <= eps*abs(p);
    ELSE
        t *= 0.5;
        REPEAT
            r := exp(t-0.25/t);
            u := r; fp := 0;
            IF mode == 3 THEN   // sinh-sinh
                r := 0.5*r - 0.5/r;
                u := 0.5*u + 0.5/u;
                y := f(c - d*r);
                IF y-y==0 THEN fp := y*u END
            ELSE                // exp-sinh
                x := c + d/r;
                IF x == c THEN break END
                y := f(x);
                IF y-y==0 THEN fp := y/u END
            END
            y := f(c + d*r);
            IF y-y == 0 THEN fp += y*u END
            fp *= (t+0.25/t);
            p, t := p+fp, t*t1;
        UNTIL abs(fp) <= eps*abs(p);
    END
    s, err, k := s+p, abs(s-p), k+1;

// print(k, d*ldexp(s,1-k), err/abs(s));
UNTIL err <= 10*eps*abs(s) OR k>n;

IF mode==1 OR (mode==3 and a>b) THEN d:=-d END
return d*ldexp(s,1-k), err/abs(s);
END;

quad(f,a,b) := quad6(f,a,b,1,6,1e-9);
#end

Usage:

Cas> pdf(x) := exp(-x*x/2) / sqrt(2*pi)

Cas> quad(pdf, -1, 1)         → [0.682689492137, 2.08159857137e−14]
Cas> quad(pdf, 1, inf)        → [0.158655253928, 1.81990772641e−10]
Cas> quad(pdf, 1, -inf)       → [-0.841344746, 7.74963867073e−11]
Cas> quad(pdf, -inf, 1)       → [ 0.841344746, 7.74963867073e−11]
Cas> quad(pdf, -inf, inf)     → [1., 8.45545855554e−13]