csaps vs. smooth from PGS
he...@web.de
The MATLAB Spline Toolbox contains a function called csaps. According
to the
documentation this function is an implementation of the fortran
function
SMOOTH from "A Practical Guide to Splines" by Carl de Boor. Now I
wanted to understand how (an if) this is true. I created a set of test
data in MATLAB and called csaps. Then I have written a program that
calls SMOOTH on the same data. Doing this I encountered several
problems, the main one is: I don't understand the documentation in
the Fortran Code.
The SMOOTH function expects parameters x and y (this seam to be the
values to be approximated), npoints (the number of points). Then dy
(estimate of uncertainty in data?) and s (upper bound on the discrete
weighted mean square distance of the approximation from the data).
The parameter s sounds like SMOOTH is more like the MATLAB function
spaps. But according to the documentation spaps uses Reinsch Splines.
I don't know what dy is about - how do I estimate the uncertainty of
my data? I have a parameter p given, such as the MATLAB function csaps
expects.
Furthermore the documentation of SMOOTH explains that p is calculated
so that sf(p) = s where sf(p) = s(f). f is
the spline function - but what is sf, and what is s(f)? The parameter
s is a scalar, so s(f) seems quite useless!
Finally, I filled dy and s with arbitrary values to test the function
(dy = 1 and s = 0). This yields yunk results, and upon closing the
application a message is given that the stack is corrupted.
Can anyone help me with these problems? I'm looking for the link or
relationship between csaps and SMOOTH - obviously csaps
is not equal to SMOOTH.
Best regards,
Torsten
PS: Please excuse my English. This text was written in German and then
translated (quite literally) to English.
Val Schmidt
I happen to be trying to work through both de Boor's book on splines and
Matlab's code at the moment and so can answer a few of your questions.
However I'll caveat my response with the fact that I am a student and am
terribly confused myself about many of the details.
The function 'smooth' is a smoothing spline function in which a cubic curve is
fit piecewise to each segment of 4 data points in the data while
simultaneously fitting a line to all the data points. The cubic curve is of
course a cubic spline which alone would pass through every data point. The
line of course will provide a mean approximation to the data. Combined they
provide both a smoothing and an interpolation of the data. In matlab, you
specify a "smoothing parameter", p that specifies how each of these two
operations is to be weighted. When p = 1, you get normal spline fit, when p
= 0 you get a linear least squares fit. You can try this with csaps in matlab
and you'll see what I mean.
However the fortran code "smooth" doesn't allow you to specify this weighting
factor, p, directly. Instead you specify the upper bound of the mean-squared
error of the fit. [Explicitly, you calculate the difference in the interpolated
function and the original data at each point, square this, weight the result by
the uncertainty in the original measurement (dy) and take the average of
those results over all data points.] If you specify a 0 mean squared error I
would expect you get a cubic spline (in which case the interpolated function
passed directly through the original data).
How do you find the uncertainty in your data? That depends on your problem.
If you've used an instrument (maybe a thermometer) to measure your data,
you might look for manufacturer's specifications. [Or you can calibrate your
thermometer by taking several readings of things at known temperatures and
calculating the standard deviation or variance in your observations.] OR if you
have multiple measurements of the same thing at each time step, you can
take use their mean in the interpolation and their standard deviation as your
uncertainty.
I'm afraid I don't have any experience with fortran to help you with the
smooth code and it's execution.
Best of luck!!!
-Val
he...@web.de wrote in message <f7fd4da6-0983-4064-a0e2-
e13956...@s13g2000prd.googlegroups.com>...
John D'Errico
<fpmlkm$3kq$1...@fred.mathworks.com>...
> Torsten,
>
> I happen to be trying to work through both de Boor's book on splines and
> Matlab's code at the moment and so can answer a few of your questions.
> However I'll caveat my response with the fact that I am a student and am
> terribly confused myself about many of the details.
Ask. I might be able to help.
> The function 'smooth' is a smoothing spline function in which a cubic curve
is
> fit piecewise to each segment of 4 data points in the data while
> simultaneously fitting a line to all the data points.
I would not say this at all.
A better description is to use a mechanical
metaphor.
A spline is a mathematical analogue of a thin,
flexible beam. (Not really exactly so if you use
the not-a-knot end conditions, but it will do
for now.)
Imagine that you have such a flexible beam
that lies in the plane. To that beam, you will
attach springs to each data point. These
springs can move only in the vertical (y)
direction, and they store an elastic potential
energy that is proportional to the square of
their extension. Each spring is hooked to your
(fixed) data points in the plane, and at their
other end, to the beam.
Now, imagine what will happen when you
create this system. The beam itself is flexible,
but it stores an elastic potential energy due to
bending in it. As well, the springs also store
potential energy due to their extension. The
shape of the curve with minimal TOTAL
potential energy is the one that will result.
Now, suppose that you have control of several
things in this beam/spring system. You can
control the relative stiffness of the beam,
compared to the spring constants in your
springs. This is what p does in the smoothing
code. When p == 1, you get an infinitely stiff
beam. P == 0 generates an perfectly flexible
beam, so the springs do not elongate at all,
and you get an interpolating spline.
What else might you do? The dy parameters
(as I remember them without pulling down
my copy of de Boor's book) effectively allow
you to change the spring constants for each
data point, making some of them stiffer than
the rest.
In more mathematical terms, a smoothing
spline is just a spline estimated using linear
regression methodology, but with a bias term
in the analysis for smoothness. If you like to
think of things as ridge regression, its just a
modified ridge estimator. For those of you
with a Bayesian flair, you can think of it in
terms of prior information, that the curve is
smooth. All of these viewpoints are appropriate.
> The cubic curve is of
> course a cubic spline which alone would pass through every data point. The
> line of course will provide a mean approximation to the data. Combined
they
> provide both a smoothing and an interpolation of the data. In matlab, you
> specify a "smoothing parameter", p that specifies how each of these two
> operations is to be weighted. When p = 1, you get normal spline fit, when
p
> = 0 you get a linear least squares fit. You can try this with csaps in matlab
> and you'll see what I mean.
>
> However the fortran code "smooth" doesn't allow you to specify this
weighting
> factor, p, directly. Instead you specify the upper bound of the mean-
squared
> error of the fit.
This is simple. The code merely finds the value
of p that gives you the desired MSE.
> [Explicitly, you calculate the difference in the interpolated
> function and the original data at each point, square this, weight the result
by
> the uncertainty in the original measurement (dy) and take the average of
> those results over all data points.] If you specify a 0 mean squared error I
> would expect you get a cubic spline (in which case the interpolated
function
> passed directly through the original data).
>
> How do you find the uncertainty in your data? That depends on your
problem.
> If you've used an instrument (maybe a thermometer) to measure your data,
> you might look for manufacturer's specifications. [Or you can calibrate your
> thermometer by taking several readings of things at known temperatures
and
> calculating the standard deviation or variance in your observations.] OR if
you
> have multiple measurements of the same thing at each time step, you can
> take use their mean in the interpolation and their standard deviation as
your
> uncertainty.
As stated, you don't always know the
uncertainty in your data. But many times
you may know it. Good experimental
practice would suggest that you know
your apparatus, and keep it under control.
You can also use replication to uncover
an estimate of the noise in your data. A
tool that I provide is estimatenoise, found
on the file exchange. It can work with
whatever data you have to infer the noise
added on top of a smooth curve.
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=16683&objectType=FILE
HTH,
John
Torsten Hensel
> However the fortran code "smooth" doesn't allow you to specify this weighting
> factor, p, directly. Instead you specify the upper bound of the mean-squared
> error of the fit.
So far, so good.
> How do you find the uncertainty in your data? That depends on your problem.
I don't want to find the uncertainty, I need to specify a given
parameter p, like I would use in csaps. This is because I don't want
to compensate for the measurement error (this is only a secondary
issue here), but rather get a smooth approximation specified by p.
> I'm afraid I don't have any experience with fortran to help you with the
> smooth code and it's execution.
That's too bad. Fortran is quite strange if you are used to
"traditional" languages like C, C++, C# or even MATLAB.
> Best of luck!!!
Thanks a lot.
Best regards,
Torsten