|
|
Contents | Previous | Next | Subchapters |
| Syntax |
smosplc(m, y, t, sig, level) m, t, c, p, k) |
| See Also | interp2 , cubespl |
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.
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).
/
------ / [ 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 ]
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
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
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.
d is a function
of the form
n(1) n(2) n(d)
P (x) = x x . . . x
i 1 2 d
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.
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
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(, 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
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.
y = sin(t ) + w
i i i
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))