Contents Previous Next Subchapters

Smoothing Splines Of Arbitrary Order And Dimension
Syntax smosplc(mytsiglevel)    
See Also interp2 , cubespl

The routine smosplc returns coefficients that define a smoothing spline that interpolates a set of measurement values. The routine smosplv uses these coefficients to evaluate the smoothing spline and its derivatives. Both smosplc and smosplv are defined in the file SMOSPL.OMS (not in SMOSPLC.OMS or SMOSPLV.OMS). A spline passes though the specified values while a smoothing spline only approximates the specified values. Smoothing splines are described in the article Surface fitting with scattered noisy data on Euclidean d-space and on the sphere, (1984) by G. Wahba, Volume 14, Number 1, of Rocky Mountain Journal of Mathematics.

The integer scalar m specifies the order of the partial derivatives that are minimized by the smoothing spline. The real or double-precision column vector y contains the measurement values that the spline is smoothing. The matrix t has the same type and number of rows as y. The i-th row of t specifies the value of the independent variable at which the value y(i) is measured. The column dimension of t is referred to as d and is the dimension of the space on which the spline is defined. If d is even, m - d must be greater than 0. If d is odd, m - d must be greater than or equal to 0. The vector sig has the same type and length as the vector y and each of the elements of sig is greater than 0. The value sig(i) is the standard deviation in the measurement y(i). The integer scalar level specifies the level of tracing. If level > 1, the values lam, R(f), the derivative of R(f) with respect to lam, and the change in lam are traced at each iteration. (a discussion of R(f) and lam follows below).

The derivative penalty term is defined by
            ------  / [   n(1)     n(2)            n(m)         ] 2
            \       | [  d        d               d             ]
     J(f) =  >      | [  -------  -------  . . .  -------  f(x) ]   dx
            /       | [    n(1)     n(2)            n(m)        ]
            ------  / [ d x      d x             d x            ]
              n    /       1        2               m           ]
where the integral is over Euclidean space of dimension d and the summation is over all vectors of integers n that have m elements each of which is between 1 and d. The average residual term is
            1 -----                   2      2
     R(f) = - >     { y  - f[t(i,:)] }  / sig  
            N -----    i                     i 
              i = 1 
where N is the number of rows in t and t(i,:) denotes the i-th row of t. The smoothing spline s(x) solves the problem

     minimize R(f) + lam J(f) with respect to f
where lam is adjusted so that the average residual, R(s), one equal to 1. If this cannot be accomplished, the message smosplc: average squared normalized residual is not 1 is printed in the Command window.

A monomial on the space of dimension d is a function of the form
              n(1)  n(2)        n(d)
     P (x) = x     x    . . .  x
      i       1     2           d
where n is a nonnegative integer vector of length d. The degree of the monomial is the sum of the element of n. The index i above is used to denote a single index counting of the different possible values for n.

The smoothing spline, s(x), is specified by the return value of smosplc. If c is the return value,
              N                               M
            -----                           -----
     s(x) = >      c  E[ |x - t(i,:)| ]  +  >      c    P (x)
            -----   i                       -----   i+N  i
            i = 1                           i = 1
where M is the number of monomials on the space of dimension d and of degree less than m, and the function E(tau) is defined by
                /             (2 m - d)
               /  log(tautau            if d is even
     E(tau) = {   
               \     (2 m - d)
                \ tau                     if d is even

Returns a column vector containing the value at the points specified by p of the derivative specified by k of the smoothing spline specified by m, t and c. The vector c is the coefficient vector returned by smosplc. The arguments m and t must be the same as used in the call to smosplc that generated the coefficient vector c. The matrix p has column dimension d and each row of p specifies a point at which to evaluate the smoothing spline. The integer row vector k has column dimension d and all of its elements are nonnegative. It specifies the derivative of the spline that we are evaluating. If all of the elements of k are 0, the spline is evaluated. In general, the i-th element of the return value is the function
      k(1)     k(2)            k(m)
     d        d               d 
     -------  -------  . . .  -------  s(x)
        k(1)     k(2)            k(m)
     d x      d x             d x 
        1        2               m 
evaluated at the point p(i,:). If one of the elements of k is nonzero. If a row of p is equal to a row of t the corresponding element of the return value is not be finite.

The following program fits a smoothing spline to 20 data points of the form
     y  = sin(t ) +  w
      i        i      i
where w is a vector of independently distributed normal random variables with mean zero and standard deviation .1. It plots the data and the smoothing spline on a 100 point grid. It then plots the derivative of the smoothing and the analytic derivative of sin(t). Note that the points where the derivative of the smoothing spline is not finite are not included in this plot (see the last sentence in the discussion of smosplv above).

include smospl.oms
include isfinite.oms
N      = 20
sig    = fill(.1, N, 1)
t      = 2 * pi * seq(N) / real(N)
y      = sin(t) + sig % snormal(N, 1)
m      = 2
level  = 1
c = smosplc(m, y, t, sig, level)
p      = 2 * pi * seq(100) / 100.
k      = 0
s      = smosplv(m, t, c, p, k)
gaddview(0., 0., 1., .5)
gplot(t, y, "plus")
gplot(p, s, "solid")
k      = 1
ds     = smosplv(m, t, c, p, k)
ok     = isfinite(ds)
p      = p.row(ok)
ds     = ds.row(ok)
gaddview(0., .5, 1., .5)
gplot(p, ds)
gplot(p, cos(p))