Fortran software is developed that calculates a best piecewise monotonic approximation to n univariate data contaminated by random errors. The underlying method minimizes the weighted sum of the squares of the errors by requiring k − 1 sign changes in the first divided differences of the approximation, where k is a given positive integer. Hence, the piecewise linear interpolant to the fit consists of k monotonic sections, alternately increasing and decreasing. This calculation can have about O ( n k ) local minima, because the positions of the turning points of the fit are integer variables of the problem. The method, however, by employing a dynamic programming technique divides the data into at most k disjoint sets of adjacent data and solves a k = 1 problem (monotonic fit or isotonic regression) for each set. So it calculates efficiently a global solution in only O ( n σ + k σ 2 ) computer operations when k ≥ 3, where σ is the number of local minima of the data, always bounded by n /2. This complexity reduces to only O ( n ) when k = 1 or k = 2 (unimodal case). At the end of the calculation a spline representation of the solution and the corresponding Lagrange multipliers are provided. The software package has been tested on a variety of data sets showing a performance that does provide in practice shorter computation times than the complexity indicates in theory. An application of the method on identifying turning points and monotonic trends of data from 1947--1996 on the U.K. pound over the U.S. dollar exchange rate is presented. Generally, the method may have useful applications as, for example, in estimating the turning points of a function from some noisy measurements of its values, or in image and signal processing, or in providing a preliminary or complementary smoothing phase to further analyses of the data.

A Fortran subroutine applies the method of Demetriou and Powell [1991] to restore convexity in n measurements of a convex function contaminated by random errors. The method minimizes the sum of the squares of the errors, subject to nonnegativity of second divided differences, in two phases. First, an approximation close to the optimum is derived in O(n) operations. Then, this approximation is used as the starting point of a dual-feasible quadratic programming algorithm that completes the calculation of the optimum. The constraints allow B-splines to be used, which reduce the problem to an equivalent one with fewer variables where the knots of the splines are determined automatically from the data points due to the constraint equations. The subroutine benefits from this reduction, since common submatrices that occur during the calculation are updated suitably. Iterative refinement improves the accuracy of some calculations when round-off errors accumulate. The subroutine has been applied to a variety of data having substantial differences and has performed fast and stably even for small data spacing, large n , and single-precision arithmetic. Driver programs and examples with output are provided to demonstrate the use of the subroutine.