Werner Pamler

2017-04-03 16:32:25 UTC

Permalink

In addition to the incomplete gamma function I am planning to add aRaw Message

series of other special functions useful for statistical calculations,

among them the incomplete beta function and its inverse.

Calculation of the inverse incomplete beta function usually is performed

in the literature by numerical root finding procedures. Numlib is

well-equipped for this purpose, it has the unit "roo" with the procedure

"roof1r" which does a bisection search of the root lying between

positive and negative guess values. The function for which the zero

point is to be found must be specified as a procedure variable of type

rfunc1r (defined in "typ", probably meaning "real function of one real

variable"); in case of the inverse incomplete beta function, this

function is the complete beta function, of course, and this depends on

three parameters, a, b, and x.

In the current state of numlib, parameters must be declared globally so

that the function passed as a parameter gets access to them:

var

_a, _b, _y: ArbFloat;

function _betai(x: ArbFloat): ArbFloat;

begin

Result := betai(_a, _b, x) - _y;

end;

function invbetai(a, b, y: ArbFloat): ArbFloat;

const

EPS = 1e-7;

var

term: ArbInt = 0;

begin

_a := a;

_b := b;

_y := y;

roof1r(@_betai, 0, 1, EPS, EPS, Result, term);

end;

This is not very nice. It would be better if the test function could be

declared as a nested function inside "invbetai":

function invbetai(a, b, y: ArbFloat): ArbFloat;

function _invbetai(x: ArbFloat): ArbFloat;

begin

Result := betai(a, b, x) - y;

end;

const

EPS = 1e-7;

var

term: ArbInt = 0;

begin

roof1r(@_betai(0, 1, EPS, EPS, Result, term);

end;

In order to make this construction pass the compiler I would like to

- declare another rfuncf1r type, now as "is nested"

type

rfuncf1rn = function(x: ArbFloat): ArbFloat is nested; //

trailing "n" means "nested"

- modify roof1r to use the nested function

Procedure roof1r(f: rfunc1rn; a, b, ae, re: ArbFloat; Var x:

ArbFloat; var term: ArbInt); overload;

- to avoid breaking existing code and to avoid forcing users of the

original version to add a {$modeswitch nestedprocvars} to their code, I

want to overload this function with a version using the non-nested

function, but calling the nested version:

Procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; Var x: ArbFloat;

Var term: ArbInt); overload;

function _f(x: ArbFloat): ArbFloat;

begin

Result := f(x);

end;

begin

roof1r(@_f, a, b, ae, re, x, term);

end;

Is there a chance that such a patch would be accepted?

_______________________________________________

fpc-devel maillist - fpc-***@lists.freepascal.org

http://lists.freepascal.org/cgi-bin/mailman/listinfo