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 convergence tasks.

_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 left bracket, then bracketmat is extended by adding
additional rows 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)