unit integrate;

(*
   PROJECT          : StatPascal library
   FILE-SPEC.       : integrate.sp
   AUTHOR           : R.-D. Reiss, M. Thomas
   CREATION DATE    : Mar-1-2000
   LAST UPDATE      : Mar-1-2000
   RELATED DOCUMENTS: Klugman, S.A. (1992):	
                      Bayesian Statistics in Actuarial Science. 
		      Kluwer, Boston.        
   VERSION          : 1.0
   FUNCTIONAL DESCR.: Unit implementing Gaussian integration.
   SPECIAL FEATURES : requires Xtremes Version 3.0
*)

interface

type realfunc = function (real): real;

function GaussIntegrate (f: realfunc; a, b: real): real;

(*****************************************************************)

implementation

function GaussIntegrate;

% Adaptive Gaussian integration taken from Klugman, S.A. (1992):	
% Bayesian Statistics in Actuarial Science. Kluwer, Boston.        

    var x, w: realvector;
        p: real;

    function integrate (a, b: real): real;
        begin
            return sum ((b - a) * 0.5 * w * f ((b + a) * 0.5 + (b - a) * 0.5 * x))
        end;

    function improve (a, b, prev, err: real; level: integer): real;
        var d, v1, v2: real;
        begin
            v1 := integrate (a, 0.5 * (a + b));
            v2 := integrate (0.5 * (a + b), b);
            d := prev - v1 - v2;
            if (abs (d) < 1e6 * err) or (level = 7) then 
                return v1 + v2
            else 
                return improve (a, 0.5 * (a + b), v1, 0.5 * err, level + 1) +
                       improve (0.5 * (a + b), b, v2, 0.5 * err, level + 1)
        end;

    begin
        x := combine (-0.973906528517172, -0.865063366688985,
                      -0.679409568299024, -0.433395394129247,
                      -0.148874338981631,  0.148874338981631,
                       0.433395394129247,  0.679409568299024,
                       0.865063366688985,  0.973906528517172);
        w := combine ( 0.066671344308688,  0.149451349050581,
                       0.219086362515982,  0.269266719309996,
                       0.295524334714753,  0.295524334714753,        
                       0.269266719309996,  0.219086362515982,
                       0.149451349050581,  0.066671344308688);
        p := integrate (a, b);
        return improve (a, b, p, p / 1e10, 1)
    end;

end.