[OctDev] octave, splines package, csaps, csaps_sel, nonequidistant sampling points

[Nir K. <nkr...@cc...>]

Dear Sebastian,

Thank you for letting me know of these problems. I think I fixed the
extrapolation now, and have posted corrected versions of csaps and
csaps_sel in the Subversion repository.

> (i) My usual call of
> "y_interp=csaps(t, y, -1, t_interp);"
> , i.e. automatic(?) relaxation between cubic spline interpolation and
> regression, produced different results for the "old" and "new" function.
> Indeed, this became a minor issue, when I noticed that the reproduction
> of the "old" results became better by applying
> "y_interp=csaps_sel(t, y, t_interp, [], "gcv");"

If I get the chance, I will take a look at the implementation in
spline-gcvspl to try and understand the reason for the difference.

> (Also, I am not sure if the documentation "p<0 or not given: an
> intermediate amount of smoothing is chosen (the smoothing term and the
> interpolation term are scaled to be of the same magnitude)" is not
> misleading.)

Why is it misleading?

> Additionally (but this is not an issue for me but might belong to any
> possible solution), there is a different method of extrapolation (csaps:
> constant, csaps_sel: spline). I am not sure if this should be mentioned
> in the documentation.

The extrapolation should be linear in both cases. I found some
mistakes in how I implemented this originally. I used the following
code (slightly non-equally-spaced points) to test the extrapolation:

x = ([1:10 10.5 11.3])'; y = sin(x);
xx = 0:0.1:12;
[yy, p] = csaps(x,y,1,xx);
plot(x, y, 's', xx, yy)

> (ii) In my point of view, there is a more severe issue: application of
> non-equidistant sampling points. In the latter case, the solution
> appears to be alright, but the first derivatives have discontinuities
> (at least for the difference quotient). In my testing file SmoothingSpline, the different
> behaviour (equidistant/non-equidistant) can be switched by
> commenting/uncommenting the line below "MAJOR ISSUE".

The first derivative of a cubic spline should be continuous and
piecewise quadratic. To get the derivative, you can use ppder. For the
test case above,

[pp, p] = csaps(x,y,1);
ppd = ppder(pp, 1);
yyd = ppval(ppd, xx);
plot(xx, yyd)

which looks to be continuous and, inside the domain, reasonably in
agreement with the derivative of sin(x).

Using csaps_sel in this case returns p=0.41 so there is some smoothing
compared with p = 1, but the derivative is still continuous:

x = ([1:10 10.5 11.3])'; y = sin(x);
xx = 0:0.1:12;
[pp, p] = csaps_sel(x,y);
yy = ppval(pp, xx);
ppd = ppder(pp, 1);
yyd = ppval(ppd, xx);
plot(xx, yyd)