Documentation |
Smoothing spline
sp = spaps(x,y,tol)
[sp,values] = spaps(x,y,tol)
[sp,values,rho] = spaps(x,y,tol)
[...] = spaps(x,y,tol,arg1,arg2,...)
[...] = spaps({x1,...,xr},y,tol,...)
sp = spaps(x,y,tol) returns the B-form of the smoothest function f that lies within the given tolerance tol of the given data points (x(j), y(:,j)), j=1:length(x). The data values y(:,j) may be scalars, vectors, matrices, even ND-arrays. Data points with the same data site are replaced by their weighted average, with its weight the sum of the corresponding weights, and the tolerance tol is reduced accordingly.
[sp,values] = spaps(x,y,tol) also returns the smoothed values, i.e., values is the same as fnval(sp,x).
Here, the distance of the function f from the given data is measured by
$$E(f)={\displaystyle \sum _{j=1}^{n}w(j)|(y(:,j)-f(x(j))){|}^{2}}$$
with the default choice for the weights w making E(f) the composite trapezoidal rule approximation to $$\int}_{\mathrm{min}(x)}^{\mathrm{max}(x)}|y-f{|}^{2$$, and |z|^{2} denoting the sum of squares of the entries of z.
Further, smoothest means that the following roughness measure is minimized:
$$F({D}^{m}f)={\displaystyle \underset{\mathrm{min}(x)}{\overset{\mathrm{max}(x)}{\int}}\lambda (t)|{D}^{m}}f(t)|dt$$
where D^{m}f denotes the mth derivative of f. The default value for m is 2, the default value for the roughness measure weight λ is the constant 1, and this makes f a cubic smoothing spline.
When tol is nonnegative, then the spline f is determined as the unique minimizer of the expression ρE(f) + F(D^{m}f), with the smoothing parameter ρ (optionally returned) so chosen that E(f) equals tol. Hence, when m is 2, then, after conversion to ppform, the result should be the same (up to roundoff) as obtained by csaps(x,y,ρ/(ρ + 1)). Further, when tol is zero, then the "natural" or variational spline interpolant of order 2m is returned. For large enough tol, the least-squares approximation to the data by polynomials of degree <m is returned.
When tol is negative, then ρ is taken to be -tol.
The default value for the weight function λ in the roughness measure is the constant function 1. But you may choose it to be, more generally, a piecewise constant function, with breaks only at the data sites. Assuming the vector x to be strictly increasing, you specify such a piecewise constant λ by inputting tol as a vector of the same size as x. In that case, tol(i) is taken as the constant value of λ on the interval (x(i-1) .. x(i)), i=2:length(x), while tol(1) continues to be used as the specified tolerance.
[sp,values,rho] = spaps(x,y,tol) also returns the actual value of ρ used as the third output argument.
[...] = spaps(x,y,tol,arg1,arg2,...) lets you specify the weight vector w and/or the integer m, by supplying them as an argi. For this, w must be a nonnegative vector of the same size as x; m must be 1 (for a piecewise linear smoothing spline), or 2 (for the default cubic smoothing spline), or 3 (for a quintic smoothing spline).
If the resulting smoothing spline, sp, is to be evaluated outside its basic interval, it should be replaced by fnxtr(sp,m) to ensure that its m-th derivative is zero outside that interval.
[...] = spaps({x1,...,xr},y,tol,...) returns the B-form of an r-variate tensor-product smoothing spline that is roughly within the specified tolerance to the given gridded data. (For scattered data, use tpaps.) Now y is expected to supply the corresponding gridded values, with size(y) equal to [length(x1),...,length(xr)] in case the function is scalar-valued, and equal to [d,length(x1),...,length(xr)] in case the function is d-valued. Further, tol must be a cell array with r entries, with tol{i} the tolerance used during the i-th step when a univariate (but vector-valued) smoothing spline in the i-th variable is being constructed. The optional input for m must be an r-vector (with entries from the set {1,2,3}), and the optional input for w must be a cell array of length r, with w{i} either empty (to indicate that the default choice is wanted) or else a positive vector of the same length as xi.
The statements
w = ones(size(x)); w([1 end]) = 100; sp = spaps(x,y, 1.e-2, w, 3);
give a quintic smoothing spline approximation to the given data that close to interpolates the first and last datum, while being within about 1.e-2 of the rest.
x = -2:.2:2; y=-1:.25:1; [xx,yy] = ndgrid(x,y); rng(39); z = exp(-(xx.^2+yy.^2)) + (rand(size(xx))-.5)/30; sp = spaps({x,y},z,8/(60^2)); fnplt(sp), axis off
produces the figure below, showing a smooth approximant to noisy data from a smooth bivariate function. Note the use of ndgrid here; use of meshgrid would have led to an error.