Weighted least-squares fitting by B-splines to discrete One-Dimensional data is performed. Constraints on the spline or its derivatives are optional. The spline function
its derivatives, or the square root of its variance function are evaluated after the fitting.
Required Arguments
DATA= DATA(1:3,:) (Input/Output) An assumed-shape array with size(data,1) = 3. The data are placed in the array: data(1,i) = xi , data(2,i) = yi, and data(3,i) = σi, i = 1,…,ndata. If the variances are not known but are proportional to an unknown value, users may set data(3,i) = 1, i = 1,…,ndata.
KNOTS= KNOTS (Input) A derived type, ?_spline_knots, that defines the degree of the spline and the breakpoints for the data fitting interval.
Optional Arguments
CONSTRAINTS= SPLINE_CONSTRAINTS (Input) A rank-1 array of derived type ?_spline_constraints that give constraints the output spline is to satisfy.
COVARIANCE = G (Output) An assumed-shape rank-2 array of the same precision as the data. This output is the covariance matrix of the coefficients. It is optionally used to evaluate the square root of the variance function.
IOPT= IOPT(:) (Input/Output) Derived type array with the same precision as the input array; used for passing optional data to SPLINE_FITTING. The options are as follows:
Packaged Options for SPLINE_FITTING
Prefix = None
Option Name
Option Value
SPLINE_FITTING_TOL_EQUAL
1
SPLINE_FITTING_TOL_LEAST
2
IOPT(IO)= ?_OPTIONS(SPLINE_FITTING_TOL_EQUAL, ?_VALUE) This resets the value for determining that equality constraint equations are rank-deficient. The default is ?_value = 10-4.
IOPT(IO)= ?_OPTIONS(SPLINE_FITTING_TOL_LEAST, ?_VALUE) This resets the value for determining that least-squares equations are rank-deficient. The default is ?_value = 10-4.
FORTRAN 90 Interface
Generic: CALLSPLINE_FITTING (DATA, KNOTS[, …])
Specific: The specific interface names are S_SPLINE_FITTING and D_SPLINE_FITTING.
Description
This routine has similar scope to CONFT found in IMSL (2003, pp 734-743). We provide the square root of the variance function, but we do not provide for constraints on the integral of the spline. The least-squares matrix problem for the coefficients is banded, with band-width equal to the spline order. This fact is used to obtain an efficient solution algorithm when there are no constraints. When constraints are present the routine solves a linear-least squares problem with equality and inequality constraints. The processed least-squares equations result in a banded and upper triangular matrix, following accumulation of the spline fitting equations. The algorithm used for solving the constrained least-squares system will handle rank-deficient problems. A set of reference are available in Hanson (1995) and Lawson and Hanson (1995). The CONFT routine uses QPROG (loc cit., p. 959), which requires that the least-squares equations be of full rank.
Fatal and Terminal Error Messages
See the messages.gls file for error messages for SPLINE_FITTING. These error messages are numbered 1340–1367.
Examples
Example 1: Natural Cubic Spline Interpolation to Data
The function
is interpolated by cubic splines on the grid of points
Those natural conditions are
Our program checks the term const. appearing in the maximum truncation error term
error≈const.×Δx4
at a finer grid.
USE spline_fitting_int
USE show_int
USE norm_int
implicit none
! This is Example 1 for SPLINE_FITTING, Natural Spline
! Interpolation using cubic splines. Use the function
write(*,*) 'Example 1 for SPLINE_FITTING is correct.'
end if
end
Output
Example 1 for SPLINE_FITTING is correct.
Example 2: Shaping a Curve and its Derivatives
The function
is fit by cubic splines on the grid of equally spaced points
The term noise is uniform random numbers from the normalized interval [-,] where = 0.01. The spline curve is constrained to be convex down for 0 ≤x≤ 1 convex upward for 1 < x≤ 4, and have the second derivative exactly equal to the value zero at x = 1. The first derivative is constrained with the value zero at x = 0 and is non-negative at the right and of the interval, x = 4. A sample table of independent variables, second derivatives and square root of variance function values is printed.
use spline_fitting_int
use show_int
use rand_int
use norm_int
implicit none
! This is Example 2 for SPLINE_FITTING. Use 1st and 2nd derivative
"The x values, 2-nd derivatives, and square root of variance.")
! See that differences are relatively small and the curve has
! the right shape and signs.
diffs=norm(values-ydata)/norm(ydata)
if (all(values > zero) .and. all(derivat1 < epsilon(zero))&
.and. diffs <= tol) then
write(*,*) 'Example 2 for SPLINE_FITTING is correct.'
end if
end
Output
Example 2 for SPLINE_FITTING is correct.
Example 3: Splines Model a Random Number Generator
The function
is an unnormalized probability distribution. This function is similar to the standard Normal distribution, with specific choices for the mean and variance, except that it is truncated. Our algorithm interpolates g(x) with a natural cubic spline, f(x). The cumulative distribution is approximated by precise evaluation of the function
Gauss-Legendre quadrature formulas, IMSL (1994, pp. 621-626), of order two are used on each polynomial piece of f(t) to evaluate q(x) cheaply. After normalizing the cubic spline so that q(1) = 1, we may then generate random numbers according to the distribution f(x≅g(x. The values of x are evaluated by solving q(x) = u, ‑1 < x < 1. Here u is a uniform random sample. Newton’s method, for a vector of unknowns, is used for the solution algorithm. Recalling the relation
we believe this illustrates a method for generating a vector of random numbers according to a continuous distribution function having finite support.
use spline_fitting_int
use linear_operators
use Numerical_Libraries
implicit none
! This is Example 3 for SPLINE_FITTING. Use splines to
! generate random (almost normal) numbers. The normal distribution
! function has support (-1,+1), and is zero outside this interval.
! Constrain the values so they fall back into the interval.
! Newton's method may give approximates outside the interval.
where(x <= -one .or. x >= one) x=zero
if(norm(fn,1) <= sqrt(epsilon(one))*norm(x,1))&
exit solve_equation
end do solve_equation
! Check that Newton's method converges.
if (niterat <= limit) then
write (*,*) 'Example 3 for SPLINE_FITTING is correct.'
end if
end
Output
Example 3 for SPLINE_FITTING is correct.
Example 4: Represent a Periodic Curve
The curve tracing the edge of a rectangular box, traversed in a counter-clockwise direction, is parameterized with a spline representation for each coordinate function, (x(t), y(t)). The functions are constrained to be periodic at the ends of the parameter interval. Since the perimeter arcs are piece-wise linear functions, the degree of the splines is the value one. Some breakpoints are chosen so they correspond to corners of the box, where the derivatives of the coordinate functions are discontinuous. The value of this representation is that for each t the splines representing (x(t), y(t)) are points on the perimeter of the box. This “eases” the complexity of evaluating the edge of the box. This example illustrates a method for representing the edge of a domain in two dimensions, bounded by a periodic curve.
use spline_fitting_int
use norm_int
implicit none
! This is Example 4 for SPLINE_FITTING. Use piecewise-linear
! splines to represent the perimeter of a rectangular box.