```help mata _bcsf_bracketing
-------------------------------------------------------------------------------

Title

_bcsf_bracketing() -- Bracketing and bracket convergence functions used
by censlope

Syntax

real scalar    _bcsf_bracketing(yfun, y, yfarleft, yfarright,
mbracket, scalefactor, bracketmat, leftindex, rightindex)

real scalar    _bcsf_bisect(objfun, x0, x1, y0, y1, itcount, iterate,
tolerance [, log])

real scalar    _bcsf_regula(objfun, x0, x1, y0, y1, itcount, iterate,
tolerance [, log])

real scalar    _bcsf_ridders(objfun, x0, x1, y0, y1, itcount,
iterate, tolerance [, log])

where

yfun:  pointer (real scalar function) scalar
y:  real scalar
yfarleft:  real scalar
yfarright:  real scalar
mbracket:  real scalar
scalefactor:  real scalar
bracketmat:  real matrix
leftindex:  real scalar
rightindex:  real scalar
objfun:  pointer (real scalar function) scalar
x0:  real scalar
x1:  real scalar
y0:  real scalar
y1:  real scalar
itcount:  real scalar
iterate:  real scalar
tolerance:  real scalar
log:  real scalar

Description

The prefix _bcsf stands for "bracket convergence for step functions".
The program censlope uses these Mata functions for the numerical solution
of inequalities involving step functions.  If a scalar object function
G() is a step function (instead of being a continuous function), then
there may not be a unique solution to the equation G(x)=0.  Instead,
there may be either no solution or a nonempty interval of solutions.
However, if G() is monotonically nonincreasing, then there may be a
supremum for the set of values of x such that G(x) is positive (or
nonnegative), which is also the infimum for values of x such that G(x) is
nonpositive (or negative).  Similarly, if G() is monotonically
nondecreasing, then there may be a supremum for the set of values of x
such that G(x) is negative (or nonpositive), which is also the infimum
for values of x such that G(x) is nonnegative (or positive).  Examples of
such step functions G() include functions derived as continuous monotonic
transformation of cumulative distribution functions (CDFs). The suprema
and infima will then be percentiles, percentile differences, percentile
ratios, or percentile slopes.  The function _bcsf_bracketing() is used to
find pairs of X values that bracket the suprema and infima.  When such a
pair of X values is found, the functions _bcsf_bisect(), _bcsf_regula(),
and _bcsf_ridders() can be used to redefine these X brackets
progressively, until the relative difference between the upper and lower
x-brackets has been reduced to a level at or below a user-defined
tolerance level. The brackets are then said to have converged.  Bracket
convergence is achieved using variants of standard numerical methods used
to find solutions for scalar equations involving continuous functions,
such as the methods described in chapter 9 of Press et al. (1992).  All
the _bcsf functions return an integer-valued scalar return code,
indicating the outcome of their attempts to carry out their bracketing or

_bcsf_bracketing(yfun, y, yfarleft, yfarright, mbracket, scalefactor,
bracketmat, leftindex, rightindex) inputs a Y value and finds a pair of X
values whose corresponding Y values bracket the input Y value.  The
arguments yfun, y, yfarleft, yfarright, mbracket, and scalefactor are
input, the arguments leftindex and rightindex are output, and the
argument bracketmat is a two-column bracket matrix, which must have at
least two rows on input and may be extended with additional rows on
output.  A bracket matrix is a two-column matrix whose first column is
assumed to contain a set of X values in ascending order and whose second
column is assumed to contain the corresponding Y values.  The argument
yfun contains a pointer to a scalar function, which is used to calculate
a Y value from the corresponding X value, and is assumed to be bounded
and monotonically nonincreasing or nondecreasing.  The argument y
contains the Y value to be bracketed.  The arguments yfarleft and
yfarright are assumed to contain the limits of (*yfun)(x) as x tends to
minus infinity and to plus infinity, respectively.  The argument mbracket
contains the maximum number of rows that the bracket matrix is allowed to
have at the end of execution.  The argument scalefactor contains a scale
factor, which must be nonmissing and strictly greater than 1, and is used
to calculate X values for additional rows of the bracket matrix.  It is
assumed that the lowest X value in the bracket matrix is negative and
that the highest X value in the bracket matrix is positive, so that an
additional X value on the left (or right) can be calculated by
multiplying the lowest (or highest) existing X value by scalefactor.  The
output arguments leftindex and rightindex are set during execution to
contain the highest index of the bracket matrix that brackets y on the
left and the lowest index of the bracket matrix that brackets y on the
right.  This bracketing must be strict, except if y == yfarleft (in which
case the left bracket may be partial) or if y == yfarleft (in which case
the right bracket may be partial).  The assumptions mentioned are assumed
and are not necessarily tested.  The terminology is detailed in
Definitions and methods below.

_bcsf_bisect(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log)
carries out bracket convergence by using a revised version (for step
functions) of the bisection method, described in Press et al. (1992).
The argument objfun is input and contains a pointer to a scalar object
function, assumed to be a monotonically nonincreasing or nondecreasing
step function.  The arguments x0 and x1 contain two X values on entry,
which are assumed to be boundaries for an interval containing the
supremum or infimum X value that we need to estimate, and will be revised
inward during each iteration to contain boundaries for another smaller
interval containing the same supremum or infimum.  The arguments y0 and
y1 are assumed to contain, on entry, the values &objfun(x0) and
&objfun(x1), respectively, and, if so, then they will still be equal to
&objfun(x0) and &objfun(x1) on output.  The argument itcount contains a
running total of iterations and will be increased on exit by the number
of iterations carried out by bcsf_bisect().  The argument iterate
contains the maximum number of iterations that bcsf_bisect() may carry
out.  The argument tolerance is input and specifies that no more
iterations are to be carried out after the value of reldif(x0,x1) is
equal to or less than the value of tolerance.  The argument log is input,
and, if it is not equal to zero, then bcsf_bisect() will output an
iteration log, giving the values of x0, x1, y0, and y1 at the end of each
iteration.  The arguments y0 and y1 must partially bracket a Y value of
zero on entry if iterations are to proceed, with y0 as the partial
bracket and y1 as the strict bracket.  If this is the case on entry, then
it will also be the case for the revised values of y0 and y1 on exit.
The terminology is detailed in Definitions and methods below.

_bcsf_regula(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log)
carries out bracket convergence by using a revised version (for step
functions) of the regula falsi (or false position) method, described in
Press et al. (1992).  The arguments have the same names and functions as
those for bcsf_bisect().

_bcsf_ridders(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log)
carries out bracket convergence by using a revised version (for step
functions) of the Ridders method, described in Press et al. (1992).  The
arguments have the same names and functions as those for bcsf_bisect().

Definitions and methods

A real number value y is said to be strictly bracketed by two other
values y0 and y1 if and only if

(y0 - y) * (y1 - y) < 0

or (in other words) if (y0 - y) and (y1 - y) have opposite nonzero signs.
If the Y values are defined using a step function, then we may use a more
relaxed definition of bracketing.  If we define

sign0 = sign(y0 - y)
sign1 = sign(y1 - y)

then we say that y is partially bracketed by y0 and y1 if and only if

sign1!=0 & sign1!=sign0

This definition implies that y0 may be equal to y but y1 may not. We then
say that y1 is the strict bracket and y0 is the partial bracket.

Whether the bracketing is strict or partial, there may be a function F()
such that F(x0) = y0 and F(x1) = y1 for some x0 and x1.  We then say that
the y-bracket associated with the lower X value is the left bracket and
that the y-bracket associated with the higher X value is the right
bracket.

The function _bcsf_bracketing() finds a pair of y-brackets for a Y value
using a bracket matrix, which is a matrix bracketmat with two columns and
at least two rows, sorted in ascending order of the first column, with
the property that, for each index i,

bracketmat[i,2] = (*yfun)(bracketmat[i,1])

In other words, the first column contains an ascending series of X
values, and the second column contains the corresponding Y values.
_bcsf_bracketing() aims to identify a pair of indices leftindex and
rightindex, having the property that bracketmat[leftindex,2] and
bracketmat[rightindex,2] bracket the input value y. This bracketing must
be strict, except if y==yfarleft (in which case the left bracket may be
partial) or if y==yfarright (in which case the right bracket may be
partial).  If the first row of bracketmat does not contain a potential
before the first row, until the first row of the extended matrix contains
a potential left bracket.  If the last row of bracketmat does not contain
a potential right bracket, then bracketmat is extended by adding
additional rows after the last row, until the last row of the extended
matrix contains a potential right bracket.  In both cases, an additional
row is defined by multiplying the X value in the existing end row by the
factor scalefactor to derive a new X value xnew, and by defining the
corresponding Y value as

ynew=(*yfun)(xnew)

and defining the new row as (xnew,ynew). The argument scalefactor is
usually set to 2, implying successive doubling of the outermost X value
in bracketmat.  This method requires that the first X value and the last
X value in bracketmat are initialized to have opposite signs.  When
bracketmat is known to contain at least one potential left bracket and
one potential right bracket, _bcsf_bracketing() returns the highest
possible index for a left bracket in leftindex and returns the lowest
possible index for a right bracket in rightindex.

The functions _bcsf_bisect(), _bcsf_regula(), and _bcsf_ridders() input a
pointer objfun, assumed to point to a nonmonotonically increasing or
decreasing function, and aim to estimate the supremum (or infimum) of the
set of X values with positive (or negative) values of (*objfun)(x).  The
functions start with a pair of X values x0 and x1, with corresponding Y
values y0=(*objfun)(x0) and y1=(*objfun)(x1), which partially bracket a Y
value of zero, with y1 as the strict bracket and y0 as the partial
bracket.  The functions then carry out a sequence of iterations whose
maximum number is given by the argument iterate, in which one or the
other of the (x,y) pairs is replaced with a new (x,y) pair, reducing the
difference abs(x0-x1) but preserving the bracketing properties.  The aim
of the iterations is to reach a point where the relative difference
reldif(x0,x1) is no more than the value of the tolerance argument, and
therefore either of the X values can be used as the estimate of the
supremum or infimum that we aim to estimate.  The methods used by the
functions are versions of those described in Press et al. (1992),
modified for use with step functions.  The function _bcsf_bisect() uses a
bisection method.  The function _bcsf_regula() uses the bisection method
if the value of y0 is zero at the start of the iteration and otherwise
uses the regula falsi method.  The function _bcsf_ridders() uses the
bisection method if the value of y0 is zero at the start of the iteration
and otherwise uses the method of Ridders (1979).

Return codes

_bcsf_bracketing() returns one of the following return codes:

0:  Bracketing successful
>
1:  Invalid row number for bracketmat on input
>
2:  Missing values in input arguments or in bracketmat on input
>
3:  First and last X values in bracketmat do not have opposite signs
>
4:  scalefactor <= 1 on input
>
5:  Left-extension of bracketmat unsuccessful
>
6:  Right-extension of bracketmat unsuccessful
>
7:  Left bracket index could not be located
>
8:  Right bracket index could not be located
>

_bcsf_bisect(), _bcsf_regula(), and _bcsf_ridders() return one of the
following return codes:

0:  Convergence successful
>
1:  Maximum iterations completed without convergence
>
2:  Missing X value calculated
>
3:  Missing Y value calculated
>
4:  Missing value for iterate on input
>
5:  y0 and y1 do not partially bracket zero on input
>

Conformability

_bcsf_bracketing(yfun, y, yfarleft, yfarright, mbracket, scalefactor,
bracketmat, leftindex, rightindex):
yfun:  1 x 1
y:  1 x 1
yfarleft:  1 x 1
yfarright:  1 x 1
mbracket:  1 x 1
scalefactor:  1 x 1
bracketmat:  M x 2 where M >= 2
leftindex:  1 x 1
rightindex:  1 x 1

_bcsf_bisect(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log),
_bcsf_regula(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log),
_bcsf_ridders(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log):
objfun:  1 x 1
x0:  1 x 1
x1:  1 x 1
y0:  1 x 1
y1:  1 x 1
itcount:  1 x 1
iterate:  1 x 1
tolerance:  1 x 1
log:  1 x 1

Diagnostics

_bcsf_bracketing(yfun, y, yfarleft, yfarright, mbracket, scalefactor,
bracketmat, leftindex, rightindex) can fail only if the function
indicated by yfun fails.

_bcsf_bisect(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log),
_bcsf_regula(objfun, x0, x1, y0, y1, itcount, iterate, tolerance , log)
and _bcsf_ridders(objfun, x0, x1, y0, y1, itcount, iterate, tolerance ,
log) can fail only if the function indicated by objfun fails.

Source code

_bcsf_bracketing.mata, _bcsf_bisect.mata, _bcsf_regula.mata,
_bcsf_ridders.mata

Author

Roger Newson, Imperial College London, UK.
Email: r.newson@imperial.ac.uk

References

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.
1992.  Numerical Recipes in C: The Art of Scientific Computing. 2nd
ed.  Cambridge, UK: Cambridge University Press.

Ridders, C. J. F. 1979. A new algorithm for computing a single root of a
real continuous function.  IEEE Transactions on Circuits and Systems,
vol.  CAS-26(11): 979-980.

Also see

Manual:  [M-0] intro

Online:  mata,
somersd (if installed)
```