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)