Contents Previous Next Subchapters Current Chapters-> interp interp1 lagrange polyfit smospl cubespl cubeval interp2 mlmode_interp2 snewton brent linlsq linlsqb nlsq dnlsq nlsqbox dnlsqb conjdir neldermead conjgrad minline lemke Qpos Qbox Qpro sqp relative fordiff cendiff autodiff testder testgrad testhess Parent Chapters-> Omatrix6 fit smospl Search Tools-> contents reference index search

Smoothing Splines Of Arbitrary Order And Dimension
 Syntax `smosplc(`m`, `y`, `t`, `sig`, `level```)     smosplv(```m`, `t`, `c`, `p`, `k`)` See Also interp2 , cubespl

Description
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, `2 m - d` must be greater than 0. If `d` is odd, `2 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 ```                 N              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(tau) tau            if d is even      E(tau) = {                   \     (2 m - d)                 \ tau                     if d is even `````` smosplv(```m`, `t`, `c`, `p`, `k```) ```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.

Example
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). ``` clear 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)) ```